C ==================================================================== C Version: 6-Sep-1988 G.P. C 13-Sep-1988 new version C 2-Dec-1993 Lahey Fortran version (S.Z.) C 22-Oct-1994 some bugs corrected C 25-Mar-1996 I filter added C 8-gru-2005 poprawka starych algorymtmow (S.Z) C ==================================================================== PROGRAM FinalRed implicit double precision(a-h,o-z) character*60 line0,line1,line2 character*30 FInp,FOut character*30 SName character *4 filtry(6) character*3 Month,Mnth1 character*1 key,ft1(10) character *1 filter parameter (ndim=5000) dimension time(ndim),flux(ndim),dmag(ndim),phase(ndim),dme(ndim) dimension zv(ndim) c ******************************* c dimension ndni(12),nmr(12) c character*3 mies(12) c data mies/'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP', c $'OCT','NOV','DEC'/ c data /nmr/1,2,3,4,5,6,7,8,9,10,11,12/ c data ndni/31,28,31,30,31,30,31,31,30,31,30,31/ c equivalence (mies(1),nmr(1)),(mies(2),nmr(2)),(mies(3),nmr(3)), c $(mies(4),nmr(4)),(mies(5),nmr(5)),(mies(6),nmr(6)), c $(mies(7),nmr(7)),(mies(8),nmr(8)),(mies(9),nmr(9)), c $(mies(10),nmr(10)),(mies(11),nmr(11)),(mies(12),nmr(12)) c ******************************* write(*,*)' INPUT Data File Name: ' read(*,'(a)')FInp write(*,*)' OUTPUT Data File Name: ' read(*,'(a)')FOut filtry(1)='U' filtry(2)='B' filtry(3)='V' filtry(4)='R' filtry(5)='I' filtry(6)='-' C C***************** Input Data of MAIN & COMP Stars ********************* C call TakeCoor(SName,xm0,p,av,dv,ac,dc) XmType=xm0+2400000. call RAConv(av) call DCConv(dv) call RAConv(ac) call DCConv(dc) open(1,file=FInp) open(2,file=FOut) open(3,file='tst1.fredm') C open(9,file='tst2.fredn') C C******************** General DATA Input ******************************* C write(*,*)'Mean Extinction Coefficients (Y/N) [Y]: ' read(*,'(a)')KEY write(*,'(1(/1x))') read(1,'(a)')line0 77 ntot=1 MODE=1 mft=0 read(1,'(a)',END=999,ERR=1000)line1 write(*,'(1x,a60)') line1 read(1,763)line2 763 format(A60) write(*,'(1x,a60)') line2 read(line1,'(6x,i2,1x,a3,1x,i4,32x,i2)')NDay,Month,NYear,NZone C ================ To jest tylko dla testow !!! =========== ihg0=50 do 345 ihg=1,10 ft1(ihg)=line2(ihg0+ihg:ihg0+ihg) 345 continue do 10 iip=1,10 10 if(ft1(iip).ne.' ') filter=ft1(iip) if(ft1(7).ne.' ') then filter=ft1(7) else filter=ft1(1) end if C write(*,'(10x,a1)') filter C ================ To ^ jest tylko dla testow !!! =========== do 654 iut=1,6 if(filter.eq.filtry(iut)) then mft=1 end if 654 continue if(mft.eq.0) then write(*,*) ' This filter name N/A - Enter new filter name: ' read(*,'(a1)') filter end if call ExtCf(filter,aext,key) write(*,765) filter 765 format(5x,25HReducing data in filter: ,A1) 33 read(1,*,ERR=555,END=555)ndat,time(ntot),flux(ntot) ntot=ntot+1 goto 33 C C************** END of MAIN DATA READING ******************************* 555 ntot=ntot-1 write(*,*) ' Number of data: ',ntot C************** START of REDUCTION ************************************* C NMon=MConv(Month) call JD(NDay,NMon,NYear,hjd,tit) C write(9,'(" tint : ",F12.5)') tint C write(9,'(" hjd : ",F12.5)') hjd if(MODE.eq.1) then BACK=-1 if(time(1).gt.11.) BACK=-BACK HJD1=hjd+BACK call JDConv(HJD1,ND1,Mnth1,NY1) endif STEP=(time(ntot)+time(1))/48. tt=tit+STEP/36525. c call HelCorr(av,dv,tt,hc) write(3,*) 'HJD1 before call hellcorrn, HJD,tt ', HJD1,HJD,tt call HelCorrn(av,dv,hjd1,hc) C write(9,*) 'hjd1', HJD1 do 11 i=1,NTOT time(i)=(time(i)-NZone)/24. if(flux(i).gt.0.) then dmag(i)=-2.5*dlog10(flux(i)) else dmag(i)=77.777 end if 11 continue C tgr=SidGr0(tit) tgr=SidGr0n(HJD,tt) write(3,'(" tgr: ",F14.8)') tgr call MainExt(time,av,dv,ac,dc,dmag,dme,tgr,aext,zv,NTOT) do 200 k=1,ntot C write(9,'(" time:",F12.5," HJD: ",F12.5, " hc: ", F8.5)') C $ time(k),hjd,hc time(k)=time(k)+hjd+hc C write(9,'(" time:",F12.5," HJD: ",F12.5, " hc: ", F8.5)') C $ time(k),hjd,hc C write(9,*) 'hjd: ', hjd C write(9,*) 'hc: ', hc 200 continue call CPhase(xm0,p,time,phase,NTOT) C C************** END of One Filter REDUCTION **************************** C*************** Write Data to OUTPUT File ***************************** C C if(MODE.gt.1) goto 111 c ---------------------------------------------------------------------- c ---------------------------------------------------------------------- if(BACK.eq.1.) then c ---------------------------------------------------------------------- write(2,'(1x,''DATE: '',i2,''-'',a3,''-'',i4,''/'',7x,''OBJECT: '' &,a30)')Nday,Month,NYear,SName IF(MONTH.EQ.'FEB') THEN nday1=NDAY+1 if(Nday.eq.28) nday=1 ELSE IF(MONTH.EQ.'APR'.OR.MONTH.EQ.'JUN'.OR.MONTH.EQ.'SEP'.OR. & MONTH.EQ. 'NOV') THEN nday1=NDAY+1 if(Nday.eq.30) nday1=1 ELSE nday1=NDAY+1 if(Nday.eq.31) nday1=1 END IF END IF write(2,'(6x,''/'',i2,''-'',a3,''-'',i4,2x,''JD0= '',f12.4,2x, &''P= '',f17.10,'' k='',f5.2)')NDay1,Mnth1,NY1,XmType,p,aext else write(2,'(1x,''DATE: '',i2,''-'',a3,''-'',i4,''/'',7x,''OBJECT: '' &,a30)')ND1,Mnth1,NY1,SName write(2,'(6x,''/'',i2,''-'',a3,''-'',i4,2x,''JD0= '',f12.4,2x, &''P= '',f17.10,'' k='',f5.2)')Nday,Month,NYear,XmType,p,aext c ---------------------------------------------------------------------- endif c ---------------------------------------------------------------------- c ---------------------------------------------------------------------- write(2,'(1x,67(''-''))') write(2,'(3x,''JDhel'',8x,''DMag'',5x,''Phase'',4x, &''ExtCorr'',2x,''ZDist'',4x,''FILTER'',2x,''No.'')') write(2,'(1x,5(''-''),'' 2400000.+ '',51(''-''))') 111 MODE=MODE+1 do 88 i=1,NTOT write(2,'(1x,f11.5,2x,f7.3,1x,f9.5,3x,f6.3,3x,f4.1,7x,a1,2x,I5)') &time(i),dme(i),phase(i),dmag(i)-dme(i),zv(i),filter,i 88 continue write(2,'(1x,67(''-''))') C C*************** Next FILTER Reduction ********************************* C backspace(1) goto 77 1000 write(*,*)' ERROR OCCURED READING INPUT FILE !!! ' goto 1001 999 continue close(2) close(1) 1001 STOP end FUNCTION MConv(Month) character*3 Month,mon(12) data mon/'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP', &'OCT','NOV','DEC'/ do 1 i=1,12 if(Month.eq.mon(i)) then mconv=i return endif 1 continue return end SUBROUTINE HelCorrn(Alf,Dec,hjd,HC) implicit double precision(a-h,o-z) C C Sun coordinates calculated according to aproximate forlmulae from the C Astronomical Almanac (S. Zola, May 2005) C pi=4.d0*DATAN(1.d0) xconv = pi/180.d0 Xn = HJD - 51545.d0 C Xn = Xn-1 write(3,*)' in helcor, before call SunCoorn, HJD, Xn: ', HJD, Xn eps = 23.439 -0.0000004 * Xn factor=0.0057755 eps = eps * xconv call SunCoorn(HJD,Slon,RSun) C write(9,*) " Slon, eps, Rsun ", slon,eps,rsun c hc=dcos(alf)*dcos(dec)*dcos(slon) c h0=(dsin(alf)*dcos(dec)*dcos(eps)+dsin(dec)*dsin(eps))*dsin(slon) c hc=hc+h0 c hc=-hc*rsun*factor hc1 = dcos(slon)*dcos(alf)*dcos(dec) hc2 = dsin(eps)*dsin(dec)+dcos(eps)*dsin(alf)*dcos(dec) hc2 = hc2 * dsin(slon) hc = factor*Rsun*(Hc1+hc2) hc = -hc write(3,*) ' ======= In HelCorr =========' write(3,*) ' hjd Lambda_Slon Rsun ' write(3,'(1x,F12.6,1x,F12.7,2x,F10.7)') HJD,Slon,Rsun write(3,'(14x,F12.7,2x,F10.7)') Slon/pi*180.,Rsun write(3,*) ' hel. corr. eps ' write(3,'(1x,F12.6,2x,F10.7)') hc,eps write(3,*) ' RA Dec ' write(3,'(1x,F12.6,1x,F12.6)') alf,dec write(3,'(1x,F12.6,1x,F12.6)')alf/pi*180.d0/360.*24.,dec/pi*180.d0 return end SUBROUTINE RAConv(alfa) C*************** Change HH.MMSS-->RAD ********************************** implicit double precision(a-h,o-z) pi=4.*datan(1.d0)/180. h=dint(alfa) h0=(alfa-h)*100. hm=dint(h0) hs=(h0-hm)*100. alfa=((hs/60.+hm)/60.+h)*15.*pi return end SUBROUTINE DCConv(dec) C*************** Change oo.''""-->RAD ********************************** implicit double precision(a-h,o-z) pi=4.*datan(1.d0)/180. SGN=1. if(dec.lt.0.) SGN=-1 h=dint(SGN*dec) h0=(SGN*dec-h)*100. hm=dint(h0) hs=(h0-hm)*100. dec=((hs/60.+hm)/60.+h)*pi*SGN return end SUBROUTINE JD(NDay,NMonth,NYear,HJD,TInt) implicit double precision(a-h,o-z) if(NMonth.le.2) then qmon=NMonth+12. year=NYear-1. else qmon=NMonth year=NYear endif a=dint(year/100.) aa=dint(a/4.) b=2.-a+aa bb=dint(365.25*year)+dint((30.6001*(qmon+1))) HJD=bb+b+NDay-679005.5 TInt=(hjd-15020.)/36525. return end SUBROUTINE SunCoorn(HJD,SLon,RSun) C C introduced in May 2005 - aproximate formulae from the C Astronomical Almanac (S. Zola) C implicit double precision (A-H,O-Z) pi=4.d0*DATAN(1.d0) pi2= 2.d0* Pi xconv = pi/180.d0 Xn = HJD - 51545.d0 C Xn = Xn -1. eps = 23.439 -0.0000004 * Xn Gdeg = 357.528 + 0.9856003 * Xn XL = 280.460 + 0.9856474 * Xn write(3,*) ' ------------- in SunCoorn ---------------------' do 133 ikl=1,1000 if(XL.lt.0.) Xl = Xl + 360.d0 if(Xl.gt.360.d0) Xl = Xl - 360.d0 if(Gdeg.lt.0.) Gdeg = Gdeg + 360.d0 if(Gdeg.gt.360.d0) Gdeg = Gdeg - 360.d0 133 continue write(3,*) ' Gdeg XL ' write(3,'(1x,F12.6,1x,F12.6)') Gdeg,Xl G = Gdeg * xconv Slon = Xl + 1.915d0*dsin(G) + 0.020d0*dsin(2.d0*G) write(3,*) ' Slon ' write(3,'(1x,F12.6)') slon C if(Slon.lt.0.d0) Slon = Slon + 360. C if(Slon.gt.360.d0) Slon = Slon - 360. Slon = Slon * xconv RSun= 1.00014 - 0.01671 * dcos(G) - 0.00014*dcos(2.d0*G) write(3,*) ' -------------------------------------------------' return end SUBROUTINE SunCoor(T,SLon,RSun) implicit double precision (A-H,O-Z) C proceudra juz zle liczy wspolrzedne Slonca C C pii=3.141592653d+00 pii=4.*datan(1.d0) pi=pii/180.d0 ch=2.*pii ql=2.8122183333d2+1.719174999d0*t+4.527777d-4*t**2+3.33d-6*t**3 qm=3.584758333d2+3.599904975277d4*t-0.00015*t**2-3.333d-6*t**3 ql=ql*pi qm=qm*pi qm=dmod(qm,ch) ql=dmod(ql,ch) em=.01675104-4.18d-5*t-1.26d-7*t**2 cc=(2.*em-.25*em**3+(5./96.)*em**5)*dsin(qm) cc=cc+(1.25*em**2-(11./24.)*em**4)*dsin(2.*qm) cc=cc+((13./12.)*em**3-(43./64.)*em**5)*dsin(3.*qm) cc=cc+((103./96.)*em**4)*dsin(4.*qm) v=qm+cc xl=v+ql ro=((1.-em**2)/(1.+em*dcos(v)))*1.00000023 a=pi*(153.23+2.25187541d4*t) b=pi*(216.57+4.50375082d4*t) c=pi*(312.69+3.29643577d4*t) d=pi*(350.74+4.452671142d5*t-.00144*t**2) e=pi*(231.19+20.20*t) h=pi*(353.4+6.59287155d4*t) a=dmod(a,ch) b=dmod(b,ch) c=dmod(c,ch) d=dmod(d,ch) e=dmod(e,ch) h=dmod(h,ch) xx=pi*(.00134*dcos(a)+.00154*dcos(b)+.002*dcos(c)) xx=xx+pi*(.00179*dsin(d)+.00178*dsin(e)) xy=5.43d-6*dsin(a)+1.575d-5*dsin(b)+1.627d-5*dsin(c) xy=xy+3.076d-5*dcos(d)+9.27d-6*dsin(h) xle=xl+xx SLon=dmod(xle,ch) RSun=ro+xy return end SUBROUTINE ExtCf(Filter,aext,key) implicit double precision(a-h,o-z) character*1 Filter,f(5),KEY dimension COEF(5) data f/'U','B','V','R','I'/ c -------------------------------------------------------- c dopisalem wsp. ekst. dla filtrow R & I (bez wyznaczania) c -------------------------------------------------------- data COEF/0.77,0.40,0.28,0.25,0.20/ if(KEY.eq.'N'.or.KEY.eq.'n') then write(*,'(a)') 'Ext. Coeff. for Filter: ',Filter read(*,*) aext return endif do 3 i=1,5 if(Filter.eq.f(i)) then aext=COEF(i) return endif 3 continue aext=COEF(2) return end SUBROUTINE CPhase(xm0,Per,hjd,phase,NTOT) implicit double precision(a-h,o-z) dimension hjd(ntot),phase(ntot) do 10 i=1,NTOT phase(i)=(hjd(i)-xm0)/Per phase(i)=phase(i)-dint(phase(i)) if(phase(i).lt.0.) phase(i)=phase(i)+1. 10 continue return end SUBROUTINE MainExt(time,av,dv,ac,dc,dmag,dme,tgr,aext,zv,NTOT) implicit double precision(a-h,o-z) dimension time(NTOT),dmag(NTOT),dme(NTOT),zv(NTOT) C C Lattitude and longitude should be in radians C C data SFi,SLag/0.872624557, -0.374877337/ KRK50cm data SFi,SLag/0.86539242,-0.35066573/ Suhora C data SFi,SLag/0.66275966,-0.39482266/ Kryonerion pi=4.*datan(1.d0)/180. Q= 1.00273790935 CONS=0.0012 cons1=0.0018167 cons2=0.002875 cons3=0.0008083 do 100 i=1,NTOT write(3,*) ' TGR, time in MainExt: ',tgr*24.,time(i)*24. tv=(tgr+Q*time(i))*24.*15.*pi-SLag tgv=tv-av tgc=tv-ac szv=1./(dsin(SFi)*dsin(dv)+dcos(SFi)*dcos(dv)*dcos(tgv)) szc=1./(dsin(SFi)*dsin(dc)+dcos(SFi)*dcos(dc)*dcos(tgc)) cosv=1./szv zv(i)=dacos(cosv)/pi C secz â^À^Ó 0.0018167(sec z -1) â^À^Ó 0.002875(sec z -1)Â? â^À^Ó 0.0008083(sec z â^À^Ó1)Â? C stare przyblizenie na airmass C xv=szv*(1.-CONS*(szv*szv-1.)) C xc=szc*(1.-CONS*(szc*szc-1.)) C C New aproximation for airmass C xv=szv-cons1*(szv-1.)-cons2*(szv -1.)**2-cons3*(szv-1.)**3 xc=szc-cons1*(szc-1.)-cons2*(szc -1.)**2-cons3*(szc-1.)**3 dx=xv-xc if(xv.lt.4.0)then de=dx*aext else de=999999. endif dme(i)=dmag(i)-de 100 continue return end FUNCTION SidGr0n(HJD0,TT) implicit double precision(a-h,o-z) C C Procedure introduced in May 2005 after prescription taken C from http://aa.usno.navy.mil/faq/docs/GAST.html C D = HJD0 + TT - 51545.0 D0 = HJD0 - 51545.0 T = D/3652 GMST = 6.697374558 + 0.06570982441908* D0 + 0.000026 * T*T do 10 i = 1,20 if(GMST.GT.24.d0) GMST = GMST - 24. if(GMST.LT.0.d0) GMST = GMST + 24. 10 continue SidGr0n = GMST /24. return end FUNCTION SidGr0(TInt) implicit double precision(a-h,o-z) C old folrmula by G. Pajdosz T=6.6460655+2400.051262*TInt+2.58622d-5*TInt*TInt if(T.lt.0.0) T=T+2400.0 SidGr0=(T-24.0*dint(T/24.0))/24. return end SUBROUTINE JDConv(hjd,NDay,Month,NYear) implicit double precision (a-h,o-z) character*3 Month,mon(12) data mon/'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP', &'OCT','NOV','DEC'/ z=dint(hjd+2400000.5) a=(dint(z-1867216.25))/36524.25 aa=dint(a/4.) b=z+a-aa+1525. c=dint((b-122.1)/365.25) d=dint(365.25*c) f=dint((b-d)/30.6001) ee=dint(30.6001*f) Nday=b-d-ee if(f.lt.13.5) qmon=f-1. if(f.ge.13.5) qmon=f-13. if(qmon.lt.2.5) NYear=c-4715. if(qmon.ge.2.5) NYear=c-4716. NMon=qmon Month=mon(NMon) if(Nday.eq.0) then hjd=hjd-1. z=dint(hjd+2400000.5) a=(dint(z-1867216.25))/36524.25 aa=dint(a/4.) b=z+a-aa+1525. c=dint((b-122.1)/365.25) d=dint(365.25*c) f=dint((b-d)/30.6001) ee=dint(30.6001*f) Nday=b-d-ee+1. if(f.lt.13.5) qmon=f-1. if(f.ge.13.5) qmon=f-13. if(qmon.lt.2.5) NYear=c-4715. if(qmon.ge.2.5) NYear=c-4716. NMon=qmon Month=mon(NMon) endif return end SUBROUTINE TakeCoor(Name,xm0,Per,av,dv,ac,dc) implicit double precision(a-h,o-z) parameter (nst=500) dimension am(nst),p(nst),a1(nst),d1(nst),a2(nst),d2(nst) character*30 OldNam(nst),Name character*73 LINE character*1 key C ================================================================ C a full path to the catalogue with data of stars declared here C ================================================================ open(7,file='./fcatal.str') 19 rewind 7 IFAT=0 write(*,*)' Enter STAR NAME [CK Boo]: ' read(*,'(a)')Name do 9 l=1,nst read(7,'(a)',END=500)OldNam(L) read(7,'(a)')LINE backspace(7) read(7,*)am(l),p(l),a1(l),d1(l),a2(l),d2(l) if(Name.eq.OldNam(l)) then xm0=am(l)-2400000. Per=p(l) av=a1(l) dv=d1(l) ac=a2(l) dc=d2(l) write(*,*)' JD0 Period R.A.var Dec.var & R.A.cmp Dec.cmp' write(*,'(a)')LINE goto 600 endif 9 continue 600 write(*,*)'Edit, Choose Another, Accept (E/C/A) [A]: ' read(*,'(a)')key if(key.eq.'c'.or.key.eq.'C') goto 19 if(key.eq.'e'.or.key.eq.'E') then IFAT=1 goto 700 endif close(7) return 500 continue OldNam(l)=Name write(*,*)' NO SUCH STAR !!!' 700 write(*,*)' Enter JD0 , PERIOD: ' read(*,*)xm0,Per am(l)=xm0 xm0=xm0-2400000. p(l)=Per write(*,*)' MAIN STAR - R.A.(hh.mmss): ' read(*,*)av a1(l)=av write(*,*)' Dec.(oo.``""): ' read(*,*)dv d1(l)=dv write(*,*)' COMP STAR - R.A.(hh.mmss): ' read(*,*)ac a2(l)=ac write(*,*)' Dec.(oo.``""): ' read(*,*)dc d2(l)=dc write(*,*)' SAVE THIS DATA (Y/N) [Y]: ' read(*,'(a)')key if(key.eq.'n'.or.key.eq.'N') goto 19 if(IFAT.eq.1)then do 29 k=l+1,nst read(7,'(a)',END=555)OldNam(k) read(7,*)am(k),p(k),a1(k),d1(k),a2(k),d2(k) 29 continue 555 l=k-1 endif rewind 7 do 99 i=1,l write(7,'(a)')OldNam(i) write(7,'(1x,f13.5,2x,f15.10,2x,f8.5,2x,f9.5,2x,f8.5,2x,f9.5)') &am(i),p(i),a1(i),d1(i),a2(i),d2(i) 99 continue close(7) return end