program gegenscheinmovies implicit none c c This program reads Paul's *.ecl files. c Form a daily median map if more than 4 orbits present c Output an antisolar latitude-longitude hemisphere map 180 x 180 degrees square (camera 1), c an antisolar latitude-longitude full sky map 360 x 180 degrees (camera 2), and c a solar-centered latitude-longitude hemisphere map 180 x 180 degrees square (camera 1). c When these are combined in either of the two fisheye.for programs, these are properly combined, although c presently the toward-sun program cannot reach beyond 90 degrees elongation for camera 3. Tough, fix it if you will! c A background-subtraction option is offered, with ecliptic latitude background subtraction. c This is presently thought to be obsolete. c Deviant responses are thrown out with fixed thresholds for cameras 1 & 2, cam 2 different in antisolar hemisphere c For camera 3, an rEloP**-3 threshold is used: this is empirically determined not to remove zodiacal light towards c the Sun, but cleans things up out to radii where this threshold is less than that for camera 2. Beyond this, the c camera 2 threshold value is used. Relaxed by 1.5x for data after DOY261 2008. c character*1 a(21),aa(80),c,asave1,asave2,asv(21) character*46 infile,outfile integer*4 i,j,k,kmax,kk,iflip,iarg,jarg,karg,indx(40),imapcount,ncam,kcount,nout,nbad,iback integer*4 iy,idoy,ih,im,is,icross,itic,imodel,islab,imask integer*1 DOY(3),HH(2),MM(2),SS(2) integer*1 ii(4),jj(4) real*4 r,rr,frame(720,360,40),medfil(720,360),array(40),arraysv(40),sum,back(360),arraycut(3),argcut,zslab,weight,amask real*4 ARRH,Omega,RRR,REW,adhocNS,adhocEW,rLngSun,rEloP,rLat,cc,ssq,s,RUSE,arg,arg1,arg2,arg3,arg4,arg5,arg6,rad,x equivalence (r,ii) equivalence (rr,jj) data c /'c'/ data infile /'G:\GegenscheinXStars\cNecl_yyyy_ddd_hhmmss.fts'/ c data outfile /'G:\GegenscheinXStars\cNecl_yyyy_ddd_hhmmss.grd'/ !get this when iback=0 c data outfile /'G:\GegenscheinXStars\cNantiyyyy_ddd_hhmmss.grd'/ !get this when iback=1 c c Flip maps to agree with Norton's & other's atlas maps? 0 = no, 1 = yes data iflip /1/ c c data ARRH,Omega,RRR,REW /5.,70.,6.,6./ !values 16 Nov 2006 data ARRH,Omega,RRR,REW /5.,90.,6.,6./ !Omega value AQ page 330 c data ARRH,Omega,RRR,REW /5.,90.,0.,6./ !***Try 0 ADU slab 28 July 2008...Omega value AQ page 330 c data ARRH,Omega,RRR,REW /5.,45.,6.,6./ !18 June 2008, undoes E/W Paul did {using Dec 2006} c data ARRH,Omega,RRR,REW /5.,45.,6.,0./ !28 June 2008, no Paul E/W: Paul did no slab at all... data back /360*0./ data imodel /0/ !imodel = 1 means "yes, subtract the model". data islab /0/ !islab = 1 means "yes, subtract the slab ". nout=0 nbad=0 kmax=40 arraycut(1)=12. c arraycut(1)=60. !for gegenschein included, on "F:\", but see the use of "arraysv" below... c arraycut(2)=20. !for antisolar hemisphere, Gegenschein movies arraycut(2)=40. !for solar hemisphere arraycut(3)=3000. icross=0 itic=0 c c iback=0 means no background subtraction, iback=1 means subtract empirical latitude correction, now (7-15-08) thought to be negligible c iback=0 if(iback.eq.0)outfile='G:\GegenscheinXStars\cNecl_yyyy_ddd_hhmmss.grd' if(iback.eq.1)outfile='G:\GegenscheinXStars\cNantiyyyy_ddd_hhmmss.grd' c ncam=3 c c Get background correction versus ecliptic latitude c if(iback.eq.1)then if(ncam.eq.1)open(10, file='G:\GegenscheinXStars\GridStorage\c1avg_2003_037_011820.dat',readonly) if(ncam.eq.2)open(10, file='G:\GegenscheinXStars\GridStorage\c2avg_2003_037_011820.dat',readonly) print *, 'Removing average background for camera ',ncam if(ncam.eq.3)print *,'WARNING -- Camera 3 has no background subtraction available' if(ncam.le.2)then do j=1,360 read(10,10)back(j) 10 format(8x,f8.1) enddo endif endif c c Get list of frames c if(ncam.eq.1)open (20, file='G:\GegenscheinXStars\namelist1.txt',readonly) if(ncam.eq.2)open (20, file='G:\GegenscheinXStars\namelist2.txt',readonly) if(ncam.eq.3)open (20, file='G:\GegenscheinXStars\namelist3.txt',readonly) 94 k=0 c 96 read (20,20)(a(i),i=1,21) 20 format(39x,21a1) c 97 k=k+1 if(k.gt.kmax)then !should never happen print *,' More than 40 orbits in a day!' nbad=nbad+1 go to 99 endif c if(a(1).ne.c)go to 98 !end of input file list if(k.eq.1)then do i=1,21 asv(i)=a(i) enddo endif c c Decode data/time c iy=ichar(a(10))-48+2000 DOY(1)=ichar(a(12))-48 DOY(2)=ichar(a(13))-48 DOY(3)=ichar(a(14))-48 HH(1)=ichar(a(16))-48 HH(2)=ichar(a(17))-48 MM(1)=ichar(a(18))-48 MM(2)=ichar(a(19))-48 SS(1)=ichar(a(20))-48 SS(2)=ichar(a(21))-48 idoy=100*DOY(1)+10*DOY(2)+DOY(3) c c Here determine the mask-in offset to be used with the "arraycut" limit c imask=0 if((iy.eq.2006.and.idoy.ge.331).or.(iy.eq.2007.and.idoy.lt.99))imask=2 if((iy.eq.2007.and.idoy.ge.256).or.(iy.eq.2008.and.idoy.lt.42).or.(iy.eq.2008.and.idoy.gt.53))imask=3 c ih=10*HH(1)+HH(2) im=10*MM(1)+MM(2) is=10*SS(1)+SS(2) rLngSun=(float(idoy)+float(ih)/24.+float(im)/1440.+float(is)/86400.)*360./365.25 rLngSun=rLngSun+270. !Subtracting Solar longitude on Jan 0 adhocNS=cosd(rLngSun-Omega-90.) adhocEW=cosd(rLngSun-Omega) c print *,iy,idoy,ih,im,is,rLngSun,adhocEW c write(infile(22:42),'(21a1)')(a(i),i=1,21) c print *,infile if(k.eq.1)then write(outfile(23:23),'(a1)')a(2) write(outfile(28:42),'(15a1)')(a(i),i=7,21) asave1=a(13) asave2=a(14) c print 1,(a(i),i=1,21) 1 format('Reading ',21a1) endif if(a(13).ne.asave1.or.a(14).ne.asave2)go to 98 open (10, file =infile, form = 'binary', readonly) do i=1,72 !read first header read (10) aa c if(i.le.61)print 100,(aa(j),j=1,78) 100 format(1x,80a1) enddo do j=1,360 !read in ecliptic-centered map read (10) (frame(i,j,k),i=1,720) c c Here swap the bytes to go from UNIX convention to Windows Absoft c do i=1,720 r=frame(i,j,k) jj(4)=ii(1) jj(3)=ii(2) jj(2)=ii(3) jj(1)=ii(4) frame(i,j,k)=rr enddo enddo close(10) if(a(13).eq.asave1.and.a(14).eq.asave2)go to 96 c 98 continue if(a(1).eq.c)then asave1=a(13) asave2=a(14) endif k=k-1 if(k.lt.5) then print *,(asv(i),i=1,21),' skip median filtering, k = ',k k=0 nbad=nbad+1 if(a(1).ne.c)go to 99 go to 97 endif if(k.lt.16)print *,(asv(i),i=1,21),' begin median filtering with',k,' files' if(k.ge.16)print *,(asv(i),i=1,21),' begin median filtering with',k,' files, ***warning, too many files***' c do i=1,720 iarg=i+360 !places antisolar point in center of plot if(ncam.eq.3)iarg=i !except for camera 3 if(iarg.gt.720)iarg=iarg-720 do j=1,360 imapcount=0 jarg=j rLat=(float(j-1)/2.)-89.75 medfil(iarg,jarg)=0. do kk=1,k array(kk)=0. arraysv(kk)=0. zslab=0. if(frame(i,j,kk).ne.0.)then array(kk)=frame(i,j,kk) x=(float(iarg-361))/2.+0.25 cc=cosd(x)*cosd(rLat) !Law of Cosines rEloP=acosd(cc) if(islab.eq.0)go to 95 !skip the slab ssq=1.-cc*cc s=sqrt(ssq) RUSE=RRR RUSE=RUSE*s !Drop off the slab brightness in antisolar hemisphere c if(adhocNS.lt.0.)then !Ad hoc N for dust slab above if(rLat.gt.ARRH)zslab=-RUSE*adhocNS if(rLat.lt.ARRH.and.rLat.gt.-ARRH)zslab=-RUSE*adhocNS*(ARRH+rLat)/(2.*ARRH) !blur edge by +-ARRH endif if(adhocNS.gt.0.)then !Ad hoc S for dust slab below if(rLat.lt.-ARRH)zslab=RUSE*adhocNS if(rLat.lt.ARRH.and.rLat.gt.-ARRH)zslab=RUSE*adhocNS*(ARRH-rLat)/(2.*ARRH) endif c RUSE=REW*abs(cc) c RUSE=RUSE*s !***This drops the EW off in the towards-solar direction, NOT what Paul did*** c if(adhocEW.lt.0..and.j.ge.180.and.iarg.ge.360)array(kk)=array(kk)-RUSE*adhocEW !West c if(adhocEW.lt.0..and.j.lt.180.and.iarg.lt.360)array(kk)=array(kk)-RUSE*adhocEW !East c if(adhocEW.ge.0..and.j.lt.180.and.iarg.ge.360)array(kk)=array(kk)+RUSE*adhocEW !West c if(adhocEW.ge.0..and.j.ge.180.and.iarg.lt.360)array(kk)=array(kk)+RUSE*adhocEW !East if(iback.eq.1)array(kk)=array(kk)-back(j)+2. array(kk)=array(kk)-zslab 95 continue arraysv(kk)=array(kk) if(imodel.eq.1)then c c The Gegenschein parameterization c arg=rLat arg1=x c rad=sqrt(arg*arg+arg1*arg1) !Old approx... c rad=180.-rElop !???Where did this come from, angle from Sun??? rad=rEloP arg2=arg/rad arg3=arg1/rad arg4=7.5*exp(-rad/4.)+39.5*exp(-rad/25.)!This one for Latitude arg5=7.5*exp(-rad/4.)+39.5*exp(-rad/35.)!This one for Longitude arg6=arg2*arg2*arg4+arg3*arg3*arg5 !sine^2 & cosine^2 weighting weight=(arg2*arg3)**2 !sine^2 of twice the angle, ugly! weight=1.-0.02*weight*rad !Want small % of sine^2 less at 30 deg., even more ugly! c*** array(kk)=array(kk)-weight*arg6 !enable this when subtracting Gegenschein from map arraysv(kk)=arraysv(kk)-weight*arg6 endif imapcount=imapcount+1 c endif enddo if(imapcount.ge.5)then karg=k-imapcount/2 call indexx(k,array,indx) sum=0. kcount=0 c c Here form the extended median by summing all below the true median, unless missing or outside of ADUs limit c This latter narrows the "lips" around the Moon and places where the data have gaps. c do kk=1,karg if(ncam.ne.3)then !cameras 1 & 2 c if(array(indx(kk)).ne.0..and.abs(arraysv(indx(kk))).lt.arraycut(ncam))then if(array(indx(kk)).ne.0..and.abs(arraysv(indx(kk))).lt.1.5*arraycut(ncam))then !used after DOY 84, 2009, but... kcount=kcount+1 sum=sum+array(indx(kk)) endif else !camera 3 argcut=arraycut(3)/(rEloP/10.)**3 if(argcut.lt.arraycut(2))argcut=arraycut(2) amask=10.*float(imask) c if(array(indx(kk)).ne.0..and.abs(arraysv(indx(kk))-amask).lt.argcut)then !used with most mask-out data & before Fall 2008 c if(array(indx(kk)).ne.0..and.abs(arraysv(indx(kk))-amask).lt.1.5*argcut)then !mask-out starting 2008 + DOY 44-55, 2009 if(array(indx(kk)).ne.0..and.abs(arraysv(indx(kk))-70.-amask).lt.1.5*argcut)then !used with mask-in data, starting Fall 2008 kcount=kcount+1 sum=sum+array(indx(kk)) endif endif enddo if(kcount.ge.1)medfil(iarg,jarg)=sum/float(kcount) !extended median c c medfil(iarg,jarg)=array(indx(karg)) !use true median c endif if(medfil(iarg,jarg).gt.99999.)medfil(iarg,jarg)=99999. if(medfil(iarg,jarg).lt.-9999.)medfil(iarg,jarg)=-9999. c c put in cross at center and ecliptic tics at sides, if desired (icross, itic = 1) c if(icross.eq.1.and.iabs(iarg-360).le.2.and.iabs(jarg-180).eq.0)medfil(iarg,jarg)=99999. if(icross.eq.1.and.iabs(iarg-360).eq.0.and.iabs(jarg-180).le.2)medfil(iarg,jarg)=99999. c if(ncam.ne.2)then if(itic.eq.1.and.iabs(iarg-181).le.3.and.iabs(jarg-180).eq.0)medfil(iarg,jarg)=99999. if(itic.eq.1.and.iabs(iarg-540).le.3.and.iabs(jarg-180).eq.0)medfil(iarg,jarg)=99999. else if(itic.eq.1.and.iabs(iarg-1) .le.3.and.iabs(jarg-180).eq.0)medfil(iarg,jarg)=99999. if(itic.eq.1.and.iabs(iarg-720).le.3.and.iabs(jarg-180).eq.0)medfil(iarg,jarg)=99999. endif enddo enddo c c write output *.grd files c if(ncam.ne.2)then open (11,file=outfile) c print 2,outfile 2 format('Writing ',a46) if(ncam.eq.1)write(11,111) 111 format('DSAA',/,'360 360',/,'-90 90',/,'-90 90',/,'-10 10') c111 format('DSAA',/,'360 360',/,'-90 90',/,'-90 90',/,'-10 50') !use this one for Gegenschein IN if(ncam.eq.3)write(11,113) 113 format('DSAA',/,'360 360',/,'-90 90',/,'-90 90',/,'-100 100') do j=1,360 if(iflip.eq.0)write(11,11)(medfil(i,j),i=181,540) if(iflip.eq.1)write(11,11)(medfil(721-i,j),i=181,540) 11 format(20f8.1) enddo close(11) else open (11,file=outfile) c print 2,outfile write(11,112) 112 format('DSAA',/,'720 360',/,'-180 180',/,'-90 90',/,'-10 10') do j=1,360 if(iflip.eq.0)write(11,11)(medfil(i,j),i=1,720) if(iflip.eq.1)write(11,11)(medfil(721-i,j),i=1,720) enddo close(11) endif nout=nout+1 k=0 if(a(1).eq.c)go to 97 99 continue print *,nout,' daily average maps written,',nbad,' orbits < 5 maps' stop 'done' end c********************************************************************************* SUBROUTINE INDEXX(N,ARRIN,INDX) !from Numerical Receipes... implicit none integer*4 N integer*4 INDX(N),I,J,L,IR,INDXT real*4 ARRIN(N),Q DO 11 J=1,N INDX(J)=J 11 CONTINUE L=N/2+1 IR=N 10 CONTINUE IF(L.GT.1)THEN L=L-1 INDXT=INDX(L) Q=ARRIN(INDXT) ELSE INDXT=INDX(IR) Q=ARRIN(INDXT) INDX(IR)=INDX(1) IR=IR-1 IF(IR.EQ.1)THEN INDX(1)=INDXT RETURN ENDIF ENDIF I=L J=L+L 20 IF(J.LE.IR)THEN IF(J.LT.IR)THEN IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1 ENDIF IF(Q.LT.ARRIN(INDX(J)))THEN INDX(I)=INDX(J) I=J J=J+J ELSE J=IR+1 ENDIF GO TO 20 ENDIF INDX(I)=INDXT GO TO 10 END