! Implemented as Model III in IDL procedure smei_zld_model.pro (PPH) program zodiac_modelb implicit none c c This program uses an ad hoc parameterization to generate a zodiac sky brightness c This is then subtracted from the selected skymap c character*46 infile,outfile character*10 head(5) character*1 charn,chars,charp real*4 theta1,rEloP,rP,PA,R,temp,tconst,doy,rLng,rLngSun,rLat real*4 EloSun,rLngI,rDis,xx,yy,RUSE real*4 Omega,RRR,REW,adhocNS,adhocEW,c,s,zod,ssq,arg,aEloP,rrrLng,rrrLat,zslab,argLat,fractNS real*8 dLngSun,dLatSun,dDisSun integer*4 i,j,iyr,idoy,hh,mm,ss,iflip,imap,ix,iy,imin,imax,jmin,jmax,jj,rr,izodimap,n,ncam integer*4 ecliptic(3600,1200),iecliptic,imark,islab,iii real*4 heliomap(3600,1800),eqmap(3600,1200),polemap(800,800) real*4 outmap(3600),rrLng,rrLat,gain data charn,chars,charp /'n','s','p'/ c data tconst /6096.312/ c data R,Omega,RRR,REW /3.,90.,3.,12./ !values 6 Nov 2006 !1st 2 values from AQ, 4th, page 330 data R,Omega,RRR,REW /5.,70.,6.,6./ !values 16 Nov 2006 data heliomap /6480000*0./ data eqmap /4320000*0./ data ecliptic /4320000*0/ data polemap /640000*0./ data iflip /1/ !iflip=1 means sky atlas presentation: note this applies to both input & output... c !what this means is that if you ran median filtering with iflip=1, you must do it here too. data imap /1/ !imap =0 means skip the sky map subtraction & output: 1=eq,2=np,3=sp data izodimap /0/ !0 means skip outputting the zodiac map data iecliptic /0/ !0 means skip outputting the ecliptic map data imark /0/ !0 means skip marking the Sun's location with a cross data ncam /1/ !choose which camera to do data islab /1/ !islab=0 means skip doing the ad hoc slab data infile /'G:\EngineeringXStars\cNeq2_200Y_DOY_HHMMSS.grd'/ !Bright-stars subtracted data data outfile /'G:\EngineeringXStars\cNeq3_200Y_DOY_HHMMSS.grd'/ !Bright-stars subtracted data c data infile /'G:\EngineeringOrbits\cNeq2_200Y_DOY_HHMMSS.grd'/ !Bright-stars not subtracted data c data outfile /'G:\EngineeringOrbits\cNeq3_200Y_DOY_HHMMSS.grd'/ !Bright-stars not subtracted data c if(imap.eq.2)then write (infile(24:24),'(a1)')charn write(outfile(24:24),'(a1)')charn write (infile(25:25),'(a1)')charp write(outfile(25:25),'(a1)')charp endif if(imap.eq.3)then write (infile(24:24),'(a1)')chars write(outfile(24:24),'(a1)')chars write (infile(25:25),'(a1)')charp write(outfile(25:25),'(a1)')charp endif c do n=1,26 !26 engineering-data cameras 1 & 2 calibrations are used covering into year 5. if((n.eq.13.or.n.eq.14.or.n.eq.15.or.n.gt.24).and.ncam.eq.3)go to 99 c if(ncam.eq.3.and.n.eq.3)go to 99 !now implicit in subroutine "getum": camera 3: discard bad-quaternions data, 2003_179. c ***Note camera 3, n=19 (2005,DOY206) has bad quaternions near G.C. Appears to be due to nearby Venus c call getum(ncam,n,iyr,idoy,hh,mm,ss,infile,outfile) c print *,'Reading',infile doy=float(idoy)+(3600.*float(hh)+60.*float(mm)+float(ss)+0.5*tconst)/86400. c c shift zodiac calculation to middle of time window for the two day (1st year) or one day calibration durations c if(iyr.eq.2003.or.(iyr.eq.2004.and.idoy.lt.100).or.(iyr.eq.2006.and.idoy.eq.56))then doy=doy+1. else doy=doy+0.5 endif c c Add a bit more for camera 2 2nd file, take away for camera 2 3rd file & camera 1 16th file, and add for camera 3 4th file. c if(ncam.eq.2.and.n.eq.2)doy=doy+0.5 if(ncam.eq.2.and.n.eq.3)doy=doy-0.35 if(ncam.eq.1.and.n.eq.16)doy=doy-0.35 if(ncam.eq.3.and.n.eq.4)doy=doy+0.2 c gain=1.00+0.01*(float(iyr-2003)+doy/365.25) !gain corrects for CCD response loss vs time: note we use 1%/yr here, not 1.6%/yr. if(ncam.eq.1)gain=gain*0.97 !camera relative gains, Chris Eyles feels these #s should be much closer to unity if(ncam.eq.3)gain=gain*0.92 !Value for 16 2006 Feb params... call SunNewcomb(0,iyr,doy,dLngSun,dLatSun,dDisSun) rDis=sngl(dDisSun) rLng=sngl(dLngSun) rLngSun = rLng rrLng=rLng rLat=sngl(dLatSun) rrLat=rLat adhocNS=cosd(rLngSun-Omega-90.) adhocEW=cosd(rLngSun-Omega) c print *,'Gain =',gain,', rDis =',rDis,', tilt =',R*cosd(rrLng-Omega),', adhoc_s = ',adhocNS,adhocEW call ECLIPTIC_EQUATOR(0,iyr,doy,rrLng,rrLat) c print *,'Sun at Ecliptic Longitude',rLng,', DOY',doy,': RA,Dec =',rrLng,rrLat if(imap.eq.1)then open(11,file=infile,readonly) read (11,11)head 11 format(a10) do j=1,1200 read (11,110)(eqmap(i,j),i=1,3600) 110 format(20f8.1) enddo close(11) print *,n,' Finished reading ',infile c c Mark the Sun's location in the eq map c if(imark.ne.0)then ix=int(1.+10.*rrLng) iy=901+int(10.*rrLat) imin=ix-3 imax=ix+3 jmin=iy-3 jmax=iy+3 if(imin.lt.1)imin=1 if(imax.gt.3600)imax=3600 do i=imin,imax do j=jmin,jmax rr=float(i-ix)**2+float(j-iy)**2 if(rr.le.9.)then if(iflip.eq.0)eqmap(i,j-300)=99999. if(iflip.eq.1)eqmap(3601-i,j-300)=99999. endif enddo enddo endif endif if(imap.gt.1)then open(11,file=infile,readonly) read (11,11)head do j=1,800 read (11,110)(polemap(i,j),i=1,800) enddo close(11) print *,n,'Finished reading ',infile endif c imax=3600 jmax=1800 if(imap.gt.1)then imax=800 jmax=800 endif c do i=1,imax if(imap.le.1)then rLngi=float(i)/10.-0.05 else xx=float(i)/10.-40.05 endif do j=1,jmax if(imap.le.1)then !equatorial map rLng=rLngi rLat=-90.05+float(j)/10. elseif(imap.eq.2)then !north polar map yy=float(j)/10.-40.05 rLng=atan2d(yy,xx) rLat=90.-sqrt(xx*xx+yy*yy) elseif(imap.eq.3)then !south polar map yy=float(j)/10.-40.05 rLng=atan2d(yy,xx) rLat=-90.+sqrt(xx*xx+yy*yy) endif call ECLIPTIC_EQUATOR(1,iyr,doy,rLng,rLat) rLng=rLng-rLngSun !Geocentric ecliptic rP=max(0.25,cosd(rLng)*cosd(rLat)) !Geocentric distance to "point P" PA=atan2d(sind(rLng),tand(rLat)) c print 1000,i,j,rLng,rLat,rP,PA 1000 format(2i6,4f10.1) rrrLng=rLng rrrLat=rLat call POINT_ON_LOS(rrrLng,rrrLat,rP,rEloP,EloSun,iii) !gets Elongation, note rrrLng & rrrLat get changed. rEloP=abs(float(iii)*rEloP) theta1=PA !+R*cosd(rLngSun-Omega) c c ****Entering Parameter Hell**** c {"Abandon all understanding, ye who enter here..."} c Exclusive of the choice of functional forms, presently I count 15 + 4 = 19 nonslab + slab parameters. c c Calculate expected zodiacal-light contribution, using the first of Steve Price's two suggested functional forms. c c This is based on Lamy and Perrin (Astr. Ap. 163, 269 (1986)), with a VSF from Ney and Merrill (Science 194, 1051 (1976)), c which integrates to a sum of cosines up through power 3, all divided by (rDis sin(rEloP))**2.3. c We have included an ad hoc 6th power cosine and exponential in the antisolar hemisphere to remove the Gegenschein residue. c aEloP=180.-rEloP c=cosd(rEloP) ssq=1.-c*c s=sqrt(ssq) c c Here alter the angular scale in like fashion to the (rDis**2.3) correction at the very end, see comment just above c ssq=s**2.3 c c zod=75. +120.*c !values 15 Nov 2006 c zod=43. +120.*c !values 15 April 2008 zod=60. +120.*c !values 15 April 2008 when killing gaussian below...*************************** if(c.le.0.)then !antisolar hemisphere c zod=zod+154.*c*c+88.*(c**3)+20.*(c**6) !values 17 Nov; don't forget c**3 is (-) here... zod=zod+154.*c*c+88.*(c**3)+22.*(c**6) !values 18 April 2008; don't forget c**3 is (-) here... c zod=zod+20.*exp(-aEloP/6.) !cutoff at 15 deg removed 27 Nov 2006 else !solar hemisphere zod=zod-185.*c*c+140.*(c**3) !values 15 Nov zod=zod/ssq endif c c Dropoff exponential in sind(rLat), with rEloP slope c The cusp is removed by using a hyperbola in place of rLat in the exponential c This has a flat bottom, pulls the rest of the curve down by 1.5 degrees, with the transition at about 1.5 degrees c This may generate another cusp at +-90 deg., hopefully too small to see...{to fix, multiply by 90/88.5} c The rEloP offset below comes from observed slopes on the rLat - rEloP plot of log10 (Observed/Model) c An offset range between about 30 to 45 is permissable. c argLat=1.5*(sqrt(1.+0.444*rLat*rLat)-1.) !hyperbola arg=abs(sind(argLat))/(0.007*(rEloP+40.)) !2 Nov 2006 values zod=zod*(10**(-arg)) c c Here include the small-scale Gegenschein term outside of the Latitude dropoff, since it is "round in the sky". zod=zod+21.*exp(-aEloP/6.) !value 18 April 2008 c c This gaussian contribution is used with camera 2 when making C2{eq,np,sp}3 grid files to be used in allskymap.for. It is c empirically found to clear out a broad ecliptic band in the eq siderial map, so the siderial residue has "no memory" of the ecliptic. c Here, "no memory" means down below a few ADUs, as best we can see... c This arbitrarily defines the "siderial component" to be subtracted from all maps, and "if it looks good, it's good". c It doesn't strictly constrain the zodiacal-model parameters, although we should probably minimize the subsequent changes. c c zod=zod+17.*exp(-rLat*rLat/400.) !gaussian contribution sigma = 14 deg., 15 April 2008, for camera 2 zod=zod+ 5.*exp(-rLat*rLat/400.) !partial gaussian contribution sigma = 14 deg., 18 April 2008, for camera 1 only c c For EloP < 90 degrees, include another contribution dropping off as cos(rLat) c if(rEloP.lt.90.)zod=zod+(1./ssq - 1.)*30.*cosd(argLat) !16 Nov 06 c if(rEloP.lt.90.)zod=zod+(1./ssq - 1.)*27.*cosd(argLat) !31 Jan 07 c c The ad hoc dust slab... c if(islab.ne.0)then zslab=0. fractNS=1. RUSE=RRR if(rEloP.gt.90.)then RUSE=RUSE*s !Drop off the slab brightness in antisolar hemisphere else RUSE=RUSE/ssq !Increase slab brightness towards Sun. fractNS=s endif c if(adhocNS.lt.0.)then !Ad hoc N for dust slab above if(rLat.gt.R)zslab=-RUSE*adhocNS if(rLat.lt.R.and.rLat.gt.-R)zslab=-RUSE*adhocNS*(R+rLat)/(2.*R) !blur edge by +-R endif if(adhocNS.gt.0.)then !Ad hoc S for dust slab below if(rLat.lt.-R)zslab=RUSE*adhocNS if(rLat.lt.R.and.rLat.gt.-R)zslab=RUSE*adhocNS*(R-rLat)/(2.*R) endif c c Close to the Sun the NS asymmetry is less and we include a symmetric slab c zslab=zslab*fractNS+(1.-fractNS)*RUSE*abs(adhocNS) c c Here's a ugly attempt at the EW/NS component, greatest when Earth passes through rotated dust plane. c There's got to be a better way than this: maybe building and integrating a 3-dimensional model isn't so bad... c arg=abs(R*adhocEW*sind(rEloP)) !bend it by the R degree tilt... RUSE=REW*abs(c) !cosine added 9 Nov to prevent E-W asymmetry, ugly... if(rEloP.gt.90.)then RUSE=RUSE*s !Drop off the slab brightness in antisolar hemisphere else RUSE=RUSE/ssq !Increase slab brightness towards Sun. endif c c Here the slab is placed on the sky. Feathered along the equatorial line by the R line tilted. c Feathered exponentially across the NS line, near zero and 180 longitude, by -longitude/R/2 antisolar and and /3 solar c This latter is the only change to date in Parameter Hell since December 2006 c if(adhocEW.lt.0.)then if(iii.eq.-1)then !West if(rLat.ge.0.)zslab=zslab-RUSE*adhocEW if(rLat.lt.0..and.rLng.gt.-180.)zslab=zslab-RUSE*adhocEW*exp(-(180.+rLng)*2./R) !feather NS line using R/2 dropoff if(rLat.lt.0..and.rLng.lt.0.)zslab=zslab-RUSE*adhocEW*exp(rLng*3./R) !feather NS line using R/3 dropoff if(rLat.lt.0..and.rLat.gt.-arg)zslab=zslab-RUSE*adhocEW*(arg+rLat)/arg !feather equatorial band endif if(iii.eq.+1)then !East if(rLat.le.0.)zslab=zslab-RUSE*adhocEW if(rLat.gt.0..and.rLng.lt.-180.)zslab=zslab-RUSE*adhocEW*exp((180.+rLng)*2./R) !feather NS line using R/2 dropoff if(rLat.gt.0..and.rLng.gt.0.)zslab=zslab-RUSE*adhocEW*exp(-rLng*3./R) !feather NS line using R/3 dropoff if(rLat.gt.0..and.rLat.lt.arg)zslab=zslab-RUSE*adhocEW*(arg-rLat)/arg !feather equatorial band endif else if(iii.eq.-1)then !West if(rLat.le.0.)zslab=zslab+RUSE*adhocEW if(rLat.gt.0..and.rLng.gt.180.)zslab=zslab+RUSE*adhocEW*exp((180.-rLng)*2./R) !feather NS line using R/2 dropoff if(rLat.gt.0..and.rLng.lt.0.)zslab=zslab+RUSE*adhocEW*exp(rLng*3./R) !feather NS line using R/3 dropoff if(rLat.gt.0..and.rLat.lt.arg)zslab=zslab+RUSE*adhocEW*(arg-rLat)/arg !feather equatorial band endif if(iii.eq.+1)then !East if(rLat.ge.0.)zslab=zslab+RUSE*adhocEW if(rLat.lt.0..and.rLng.lt.180.)zslab=zslab+RUSE*adhocEW*exp(-(180.-rLng)*2./R) !feather NS line using R/2 dropoff if(rLat.lt.0..and.rLng.gt.0.)zslab=zslab+RUSE*adhocEW*exp(-rLng*3./R) !feather NS line using R/3 dropoff if(rLat.lt.0..and.rLat.gt.-arg)zslab=zslab+RUSE*adhocEW*(arg+rLat)/arg !feather equatorial band endif endif c zod=zod+zslab endif c c note when rEloP=90 and rLat=90 (ecliptic poles), this + DC offset below should agree with AQ p329 (22.4), or p330 (24) c Nov 2's value is 6 {from the 75 constant above top line for zod, times 10**(-1/(.007*130))}; but see just below c 15 April's value is 3.4 {from the 43 constant above top line for zod, times 10**(-1/(.007*130))}; but see just below c Then we add the below number for the comparison, works out fine. c c zod=zod+16. !DC offset, a case can be made that this should be "zod=zod+REW" zod=zod+19. !value 14 April 2008 c c Diminish zodical light a bit near ecliptic and around 90 degrees, this probably needs refinement, especially towards the Sun c c zod=zod-8.*sind(rEloP)*exp(-rLat*rLat/10.) !added 17 Nov 2006 c c c ****Leaving Parameter Hell**** c temp=zod heliomap(i,j)=temp/(rDis**2.3) !Correction to 1 AU: heliomap is reduced for > 1 AU and v.v. c if(heliomap(i,j).gt.99999.)heliomap(i,j)=99999. if(heliomap(i,j).lt.-9999.)heliomap(i,j)=-9999. if(abs(abs(PA)-90.).lt.0.1.and.j.gt.300.and.j.le.1500)ecliptic(i,j-300)=1 !mark ecliptic enddo enddo print *,'writing output files' if(imap.eq.1.and.izodimap.ne.0)then if(n.eq.1)open (10,file='G:\engineeringxstars\zodiac1.grd') if(n.eq.2)open (10,file='G:\engineeringxstars\zodiac2.grd') if(n.eq.3)open (10,file='G:\engineeringxstars\zodiac3.grd') if(n.eq.4)open (10,file='G:\engineeringxstars\zodiac4.grd') if(n.eq.5)open (10,file='G:\engineeringxstars\zodiac5.grd') if(n.eq.6)open (10,file='G:\engineeringxstars\zodiac6.grd') if(n.eq.7)open (10,file='G:\engineeringxstars\zodiac7.grd') if(n.eq.8)open (10,file='G:\engineeringxstars\zodiac8.grd') if(n.eq.9)open (10,file='G:\engineeringxstars\zodiac9.grd') if(n.eq.10)open (10,file='G:\engineeringxstars\zodiac10.grd') if(n.eq.11)open (10,file='G:\engineeringxstars\zodiac11.grd') if(n.eq.12)open (10,file='G:\engineeringxstars\zodiac12.grd') if(n.eq.13)open (10,file='G:\engineeringxstars\zodiac13.grd') if(n.eq.14)open (10,file='G:\engineeringxstars\zodiac14.grd') if(n.eq.15)open (10,file='G:\engineeringxstars\zodiac15.grd') if(n.eq.16)open (10,file='G:\engineeringxstars\zodiac16.grd') if(n.eq.17)open (10,file='G:\engineeringxstars\zodiac17.grd') if(n.eq.18)open (10,file='G:\engineeringxstars\zodiac18.grd') if(n.eq.19)open (10,file='G:\engineeringxstars\zodiac19.grd') if(n.eq.20)open (10,file='G:\engineeringxstars\zodiac20.grd') if(n.eq.21)open (10,file='G:\engineeringxstars\zodiac21.grd') if(n.eq.22)open (10,file='G:\engineeringxstars\zodiac22.grd') if(n.eq.23)open (10,file='G:\engineeringxstars\zodiac23.grd') if(n.eq.24)open (10,file='G:\engineeringxstars\zodiac24.grd') if(n.eq.25)open (10,file='G:\engineeringxstars\zodiac25.grd') if(n.eq.26)open (10,file='G:\engineeringxstars\zodiac26.grd') write(10,10) endif if(imap.ne.0)open (11,file=outfile) 10 format('DSAA'/'3600 1800'/'0 360'/'-90 90'/'0 1000') head(5)='0 1000 ' if(imap.ge.1)write(11,11)head if(izodimap.ne.0.or.imap.eq.1)then do j=1,1800 jj=j-300 if(imap.eq.1.and.izodimap.ne.0)then if(iflip.eq.0)write(10,110)(heliomap(i,j),i=1,3600) if(iflip.eq.1)write(10,110)(heliomap(3601-i,j),i=1,3600) endif if(imap.eq.1.and.jj.gt.0.and.jj.le.1200)then do i=1,3600 outmap(i)=eqmap(i,jj) if(outmap(i).ne.99999.)outmap(i)=outmap(i)*gain if(iflip.eq.0)then if(outmap(i).ne.0..and.outmap(i).ne.99999.)outmap(i)=outmap(i)-heliomap(i,j) endif if(iflip.eq.1)then if(outmap(i).ne.0..and.outmap(i).ne.99999.)outmap(i)=outmap(i)-heliomap(3601-i,j) endif if(outmap(i).gt.99999.)outmap(i)=99999. if(outmap(i).lt.-9999.)outmap(i)=-9999. enddo write(11,110)(outmap(i),i=1,3600) endif enddo endif if(imap.ge.2)then do j=1,800 do i=1,800 outmap(i)=polemap(i,j) if(outmap(i).ne.99999.)outmap(i)=outmap(i)*gain if(iflip.eq.0)then if(outmap(i).ne.0..and.outmap(i).ne.99999.)outmap(i)=outmap(i)-heliomap(i,j) endif if(iflip.eq.1.and.imap.eq.2)then if(outmap(i).ne.0..and.outmap(i).ne.99999.)outmap(i)=outmap(i)-heliomap(i,801-j) endif if(iflip.eq.1.and.imap.eq.3)then if(outmap(i).ne.0..and.outmap(i).ne.99999.)outmap(i)=outmap(i)-heliomap(801-i,801-j) endif if(outmap(i).gt.99999.)outmap(i)=99999. if(outmap(i).lt.-9999.)outmap(i)=-9999. enddo write(11,110)(outmap(i),i=1,800) enddo endif if(izodimap.ne.0)close(10) if(imap.ne.0)close(11) 99 continue enddo !end of n-loop if(iecliptic.ne.0)then open (12,file='G:\EngineeringOrbits\ecliptic.grd') write(12,100) 100 format('DSAA'/'3600 1200'/'0 360'/'-60 60'/'-1 1') do j=1,1200 if(iflip.eq.0)write(12,111)(ecliptic(i,j),i=1,3600) 111 format(80i2) if(iflip.eq.1)write(12,111)(ecliptic(3601-i,j),i=1,3600) enddo close(12) endif print *, 'done' stop end c********************************************************************************* subroutine getum(ncam,n,iyr,idoy,hh,mm,ss,infile,outfile) c c This messy program implicitly defines the year, days, hours, minutes, and seconds which name each filtered skymap c implicit none character*46 infile,outfile integer*4 ncam,n,iyr,idoy,hh,mm,ss,idum,idum1,idum2,idum3 integer*4 idoy1(26),ihh1(26),imm1(26),iss1(26),idoy2(26),ihh2(26),imm2(26),iss2(26),idoy3(26),ihh3(26),imm3(26),iss3(26) character*1 dum c data idoy1 /56,127,183,242,319,361,55,113,172,238,291,345,41,101,184,247,304,364,61,149,204,264,323,76,138,194/ data ihh1 /16, 15, 11, 15, 15, 01,12, 16, 17, 16, 01, 00,00, 00, 01, 01, 01, 00,00, 00, 00, 00, 01,23, 03, 00/ data imm1 /05, 41, 10, 56, 25, 25,54, 12, 30, 44, 10, 31,54, 10, 06, 10, 17, 30,50, 14, 55, 07, 18,57, 55, 53/ data iss1 /14, 11, 38, 25, 03, 01,03, 13, 21, 44, 02, 41,32, 30, 55, 04, 09, 52,08, 43, 32, 32, 19,58, 30, 23/ c data idoy2 /58,124,182,240,311,358,53,112,171,241,290,346,38,100,185,248,305, 01,56,148,206,266,325,79,137,196/ data ihh2 /15, 16, 04, 16, 16, 15,01, 16, 17, 04, 01, 00,00, 00, 00, 00, 00, 01,00, 00, 00, 01, 00,01, 05, 00/ data imm2 /30, 33, 41, 31, 04, 50,38, 29, 47, 00, 27, 14,05, 28, 49, 52, 59, 37,37, 32, 20, 13, 42,04, 54, 17/ data iss2 /14, 46, 45, 31, 01, 35,03, 50, 58, 40, 41, 02,56, 11, 14, 22, 27, 01,07, 27, 05, 39, 50,03, 50, 51/ c data idoy3 /60,129,194,198,238,316,356,50,111,170,240,296,347,39,97,187,246,306,365,58,147,205,265,324,078,136/!,195,325/ data ihh3 /16, 16, 01, 00, 15, 16, 16,17, 15, 18, 04, 14, 01,01,06, 00, 01, 04, 00,00, 00, 00, 01, 01, 01, 04/!, 00, 16/ data imm3 /36, 47, 11, 01, 25, 17, 25,45, 05, 05, 18, 56, 37,29,25, 13, 27, 04, 13,01, 50, 37, 31, 00, 21, 31/!, 35, 16/ data iss3 /51, 44, 18, 09, 01, 46, 46,15, 51, 35, 18, 13, 59,52,58, 51, 45, 56, 09,41, 10, 48, 23, 35, 48, 01/!, 37, 12/ c c Insert camera number and year c iyr=2003 if(ncam.eq.1)then if(n.gt.6)iyr=iyr+1 if(n.gt.12)iyr=iyr+1 if(n.gt.18)iyr=iyr+1 if(n.gt.23)iyr=iyr+1 if(n.gt.26)iyr=iyr+1 endif if(ncam.eq.2)then if(n.gt.6)iyr=iyr+1 if(n.gt.12)iyr=iyr+1 if(n.gt.17)iyr=iyr+1 if(n.gt.23)iyr=iyr+1 if(n.gt.26)iyr=iyr+1 endif if(ncam.eq.3)then if(n.gt.7)iyr=iyr+1 if(n.gt.13)iyr=iyr+1 if(n.gt.19)iyr=iyr+1 if(n.gt.24)iyr=iyr+1 if(n.gt.28)iyr=iyr+1 endif c dum=char(ncam+48) write (infile(23:23),'(a1)')dum write(outfile(23:23),'(a1)')dum dum=char(iyr-2000+48) write (infile(31:31),'(a1)')dum write(outfile(31:31),'(a1)')dum c c Insert idoy, hour, minutes, seconds c if(ncam.eq.1)then idoy=idoy1(n) hh=ihh1(n) mm=imm1(n) ss=iss1(n) endif if(ncam.eq.2)then idoy=idoy2(n) hh=ihh2(n) mm=imm2(n) ss=iss2(n) endif if(ncam.eq.3)then idoy=idoy3(n) hh=ihh3(n) mm=imm3(n) ss=iss3(n) endif c c Conversion to strings: big pain! c idum=idoy/100 dum=char(idum+48) write (infile(33:33),'(a1)')dum write(outfile(33:33),'(a1)')dum idum1=idoy-100*idum idum2=idum1/10 dum=char(idum2+48) write (infile(34:34),'(a1)')dum write(outfile(34:34),'(a1)')dum idum3=idoy-100*idum-10*idum2 dum=char(idum3+48) write (infile(35:35),'(a1)')dum write(outfile(35:35),'(a1)')dum c idum=hh/10 dum=char(idum+48) write (infile(37:37),'(a1)')dum write(outfile(37:37),'(a1)')dum idum1=hh-10*idum dum=char(idum1+48) write (infile(38:38),'(a1)')dum write(outfile(38:38),'(a1)')dum c idum=mm/10 dum=char(idum+48) write (infile(39:39),'(a1)')dum write(outfile(39:39),'(a1)')dum idum1=mm-10*idum dum=char(idum1+48) write (infile(40:40),'(a1)')dum write(outfile(40:40),'(a1)')dum c idum=ss/10 dum=char(idum+48) write (infile(41:41),'(a1)')dum write(outfile(41:41),'(a1)')dum idum1=ss-10*idum dum=char(idum1+48) write (infile(42:42),'(a1)')dum write(outfile(42:42),'(a1)')dum return end c********************************************************************************* subroutine rotate(ALFA,BETA,GAMMA,PHI,RLAT) implicit none c INPUTS: c ALFA real phase angle of the pole Z(new) in old coordinate system c BETA real polar angle of the pole Z(new) in old coordinate system c GAMMA real phase angle of X(new) in coordinate system (2) c PHI real phase angle in old coordinate system c RLAT real latitude in old coordinate system (-90<=RLAT<=90) c OUTPUTS: c PHI real phase angle in new coordinate system (0<=PHI<360) c RLAT real latitude in new coordinate system real ALFA real BETA real GAMMA real PHI real RLAT double precision SINB double precision COSB double precision SINL double precision COSL double precision SINP double precision COSP double precision COST ! Needed for g77 when using goniod.f ! What happens if these are available as intrinsic functions?? double precision dsind double precision dcosd double precision dacosd double precision datan2d if (abs(RLAT) .gt. 90.0)print *,'Bummer: RLAT = ',RLAT !call Say('rotate','E','invalid','latitude') PHI = PHI-ALFA ! Rotation (1) SINB = dsind(dble(BETA)) ! Rotation (2) COSB = dcosd(dble(BETA)) SINL = dsind(dble(RLAT)) COSL = dcosd(dble(RLAT)) SINP = dsind(dble(PHI )) COSP = dcosd(dble(PHI )) COST = SINL*COSB+COSL*SINB*COSP ! cos(polar angle) if (dabs(COST) .ne. 1.0d0) then ! changed abs --> dabs from Paul's, probably doesn't matter... SINP = COSL*SINP COSP = COSL*COSB*COSP-SINL*SINB else !Point on new Z-axis; PHI not defined !------- ! The following statement is inserted for the benefit of the ! synoptic map drawing program. ! If BETA=0 then SINB=0, COSB=+/-1, i.e. the z-axis stays the same ! or is flipped by 180 degrees. The input phase angle for points on ! the z-axis will remain the same or is changed to 180-(phase angle) ! (just as for points not on the z-axis). if (SINB .eq. 0.0d0) COSP = COSB*COSP end if PHI = datan2d(SINP,COSP) PHI = PHI-GAMMA ! Rotation (3) PHI = amod(amod(PHI,360.0)+360.0,360.0)!0<=PHI<360 RLAT = 90.0d0-dacosd(COST) ! acosd returns between 0 and 180 deg return end c********************************************************************************* subroutine ECLIPTIC_EQUATOR(ID,iYr,Doy,PHI,RLAT) implicit none cINPUTS: c ID integer ID=0: ecliptic ---> equatorial c ID=1: equatorial ---> ecliptic c iYr integer year of current date c Doy real day of year, including a fraction for the time of day c PHI real ecliptic longitude for equinox of current date c RLAT real ecliptic latitude for equinox of current date c OUTPUTS: c PHI real right ascension c RLAT real declination integer ID integer iYr real Doy real PHI,BETA !BETA added to this list to satisfy implicit none real RLAT double precision EPS double precision JD double precision JEpoch call Julian(0,iYr,Doy,JD,JEpoch) !Get date in Julian years EPS = 23.439291d0-0.00013004d0*(JEpoch-2000.0d0) !Angle (ecliptic-equator) BETA = sngl(EPS) ! .. to equatorial coord if (ID .eq. 1) BETA = -BETA ! .. to ecliptic coord call rotate(90.0,BETA,-90.0,PHI,RLAT) return end c********************************************************************************* subroutine Julian(IDIN,iYr,Doy,JDio,JEpoch) implicit none c INPUTS: c IDIN integer ID=0 date --> JD, Julian epoch, Besselian epoch c =1 JD --> Date, Julian epoch, Besselian epoch c =2 Julian epoch --> Date, JD, Besselian epoch c (not active) =3 Besselian epoch --> Date, JD, Julian epoch c Add 10: use modified Julian days (=JD-2400000.5) c Add 20: use days relative to Jan. 1, 2000 (=JD-2451544.5) c INPUTS/OUTPUTS: (depending on value of ID) c iYr integer year; the year xxxBC should be entered as -xxx+1. c Doy real day of year (including fraction for the time of day). c JD double precision Julian day c JEpoch double precision Julian epoch = time in Julian years c*******data type specifications are here added on to satisfy implicit none******* integer IDIN,ID,I,J,Laps,Leap,iDoy,JD0,I365,I1582 integer iYr real Doy double precision JDio double precision JEpoch parameter (I365 = 365) parameter (I1582 = -418) !1582-2000 double precision JD,Frac double precision BEpoch double precision JD2000 JD2000 = 2451545 ! Jan 1.5, 2000 (noon) if (IDIN .ge. 10) JD2000 = 51544.5d0 if (IDIN .ge. 20) JD2000 = 0.5d0 ID = mod(IDIN,10) ! ID=0,1,2,3 if (ID .eq. 0) then ! Date ---> Julian day, etc. Laps = iYr-2000 ! Whole years from 1 Jan 1.0d 2000 to idem of iYr I = Laps if (Laps .gt. 0) I = Laps-1 J = max(I1582,I) Leap = I/4-J/100+J/400 if (Laps .gt. 0) Leap = Leap+1 ! # leap years JD = I365*Laps+Leap+dble(Doy-1.5) ! Rel. to 2000 Jan. 1.5 (noon) if (iYr .le. 1582) JD = JD+10 JEpoch = 2000+JD/365.25d0 BEpoch = 1900+(JD+36524.68648d0)/365.242198781d0 else ! Julian day, etc. ---> Date if (ID .eq. 2) then ! Julian epoch JD = JEpoch-2000 BEpoch = 2000.0012775d0+1.000021359d0*JD JD = JD*365.25d0 else if (ID .eq. 3) then JD = BEpoch-1900 JEpoch = 1900.000858d0+0.999978642d0*JD JD = -36524.68648d0+JD*365.242198781d0 else ! Assume ID=1 JD = JDio-JD2000 JEpoch = 2000+JD/365.25d0 BEpoch = 1900+(JD+36524.68648d0)/365.242198781d0 end if JD = JD+0.5d0 ! Shift to previous midnight JD0 = int(JD) ! Whole days between previous midnight ... if (JD .lt. JD0) JD0 = JD0-1! .. and midnight of 1 Jan 2000 Frac = JD-JD0 ! Day fraction from previous midnight Laps = 0 iDoy = JD0 do while (iabs(iDoy) .gt. 365) Laps = Laps+iDoy/366 ! Underestimate number of full years in iDoy I = Laps if (Laps .gt. 0) I = Laps-1 J = max(I1582,I) Leap = I/4-J/100+J/400 if (Laps .gt. 0) Leap = Leap+1 ! # leap years iDoy = JD0-Laps*I365-Leap if (Laps .le. I1582) iDoy = iDoy-10 end do !----- ! Leave do while loop with -365<=iDoy<=365. Doy <= 0 happens for dates earlier ! than 2000 Jan. 1.0; Doy >=0 for dates later than 2000 Jan 1.0. iYr = 2000+Laps if (iDoy .lt. 0) then iYr = iYr-1 if (iYr .eq. 1582) then iDoy = iDoy+355 ! 355: # days in 1582 if (iDoy .lt. 0) iYr = iYr-1 end if end if Leap = 0 ! Suppose it's not a leap year if (iYr/4*4 .eq. iYr .and. (iYr .lt. 1582 .or. !It's a leap year & iYr/100*100 .ne. iYr .or. iYr/400*400 .eq. iYr)) Leap = 1 if (iDoy .lt. 0) iDoy = 365+Leap+iDoy ! Now iDoy>=0 if (iDoy .eq. 365+Leap) then iDoy = 0 iYr = iYr+1 end if Doy = iDoy+1+Frac JD = JD-0.5D0 end if JDio = JD2000+JD return end c********************************************************************************* subroutine POINT_ON_LOS(XLNG,XLAT,RP,ELOLD,ELNEW,iEorW) implicit none cINPUTS: c (either topocentric or heliocentric) c c XLNG real longitude (0<=XLNG<=360) and ... c XLAT real .. latitude (-90<=XLAT<=90) of line of sight (degrees) c RP real radial distance to point P c (in units of the observer-Sun distance) c OUTPUTS: c (either heliocentric or topocentric) c c XLNG real longitude (0<=XLNG<360) and ... c XLAT real .. latitude (-90<=XLAT<=90) of point P (degrees) c RP real distance of point P in units of observer-Sun distance c ELOLD real elongation in the old coordinate system (deg); 0<=ELOLD<=180 c ELNEW real elongation in the new coordinate system 0<=ELNEW<=180 c iEorW integer +1 or -1; indicates on which side (East or West) of c the Sun-Earth line the point is located real XLNG real XLAT real RP real ELOLD real ELNEW integer iEorW double precision DLNG double precision DLAT double precision DRP double precision R double precision cosE double precision sinP double precision cosP ! Needed for g77 when using goniod.f ! What happens if these are available as intrinsic functions?? double precision dsind double precision dcosd double precision dasind double precision dacosd double precision datan2d DLNG = dble(XLNG) ! Convert to double DLAT = dble(XLAT) DRP = dble(RP ) if (RP .lt. 0.0) then ! Negative RP: DLNG = 180.0d0+DLNG ! .. change to opposite direction DLAT = -DLAT DRP = -DRP end if cosE = DLNG if (dabs(cosE) .gt. 180.0d0) cosE = cosE-dsign(1.0d0,cosE)*360.0d0 iEorW = dnint( dsign(1.0d0,cosE) ) cosE = dcosd(DLNG)*dcosd(DLAT)! Cos(old elongation) ELOLD = dacosd(cosE) ! Old elongation in [0,180] R = 1d0+DRP*(DRP-2d0*cosE)! Cosine rule: (dist Sun-P)**2 if (R .gt. 0d0) then R = dsqrt(R) ! Dist Sun-P RP = R ! Real*4 output sinP = dsind(DLAT)*DRP/R ! Sin(Latitude) sinP = dmin1(1d0,dmax1(-1d0,sinP)) XLAT = sngl(dasind(sinP)) ! Latitude [-90,90] cosP = 1d0-DRP*cosE cosP = dmin1(R,dmax1(-R,cosP)) ! prop. cos(XLNG) ELNEW = dacosd(cosP/R) ! New elongation sinP = -DRP*dcosd(DLAT)*dsind(DLNG) ! prop. sin(XLNG) XLNG = sngl(datan2d(sinP,cosP)) XLNG = amod(amod(XLNG,360.0)+360.0,360.0)! 0<=XLNG<360 else RP = 0.0 ! (just to prevent 'divide by zero' XLAT = 0.0 ! condition) XLNG = 0.0 ELOLD = 0.0 ELNEW = 0.0 end if return end c********************************************************************************* subroutine SunNewcomb(ID,iYr,Doy,rLng,rLat,rDis) implicit none cINPUTS: c ID integer 0 = return two-body solution c 1 = take planetary perturbations into account c iYr integer year c Doy real doy of year, including fraction for time of day c OUTPUTS: c rLng double precision ecliptic longitude in degrees [0,360) c rLat double precision ecliptic latitude in arcsec c rDis double precision Sun-Earth distance in AU integer ID integer iYr real Doy double precision rLng double precision rLat double precision rDis double precision JD double precision T double precision DLP double precision G double precision G2 double precision G4 double precision G5 double precision G6 double precision D double precision A double precision U double precision DLO double precision SLO double precision DL double precision R1 ! Needed for g77 when using goniod.f ! What happens if these are available as intrinsic functions?? double precision dsind double precision dcosd call Julian(0,iYr,Doy,JD,T) T = (T-1900)/100 ! Time since 0.5 Jan 1900 (Julian centuries) DLP = (1.882d0-.016d0*T)*dsind( ( 57.24d0+150.27d0*T) ) & + 6.40d0*dsind( (231.19d0+ 20.20d0*T) ) & + .266d0*dsind( ( 31.80d0+119.00d0*T) ) G = 358.4758333d0 + mod(35999.d0*T,360.d0) + (T*(179.10d0-.54d0*T)+DLP)/3600.d0 !------- ! LO = DLO+SLO DLO = 279.6966778d0 + mod(36000.d0*T,360.d0) SLO = T*(2768.13d0+1.089d0*T) + .202d0*dsind( 315.6d0+893.3d0*T ) + DLP DL = (6910.057d0-17.240d0*T)*dsind( G) & + ( 72.338d0- 0.361d0*T)*dsind(2.d0*G) & + 1.054d0*dsind(3.d0*G) rDis =( .00003057d0-.00000015d0*T) & + (-.00727412d0+.00001814d0*T)*dcosd( G ) & + (-.00009138d0+.00000046d0*T)*dcosd( 2.d0*G ) & + (-.00000145d0) *dcosd( 3.d0*G ) rLat = 0.d0 if (ID .eq. 1) then G2 = 212.45d0+58517.493d0*T G4 = 319.58d0+19139.977d0*T G5 = 225.28d0+ 3034.583d0*T & + 1300.d0*dsind( 133.775d0+39.804d0*T )/3600.d0 G6 = 175.60d0+ 1221.794d0*T D = 350.737486d0+445267.114217d0*T A = 296.104608d0+477198.849108d0*T U = 11.250889d0+483202.025150d0*T !------- ! DL2 (Venus) DL = DL & +4.838d0*dcosd( 299.102d0+ G2- G ) & +0.116d0*dcosd( 148.900d0+2.d0*G2- G ) & +5.526d0*dcosd( 148.313d0+2.d0*G2-2.d0*G ) & +2.497d0*dcosd( 315.943d0+2.d0*G2-3.d0*G ) & +0.666d0*dcosd( 177.710d0+3.d0*G2-3.d0*G ) & +1.559d0*dcosd( 345.253d0+3.d0*G2-4.d0*G ) & +1.024d0*dcosd( 318.150d0+3.d0*G2-5.d0*G ) & +0.210d0*dcosd( 206.200d0+4.d0*G2-4.d0*G ) & +0.144d0*dcosd( 195.400d0+4.d0*G2-5.d0*G ) & +0.152d0*dcosd( 343.800d0+4.d0*G2-6.d0*G ) & +0.123d0*dcosd( 195.300d0+5.d0*G2-7.d0*G ) & +0.154d0*dcosd( 359.600d0+5.d0*G2-8.d0*G ) !------- ! DL4 (Mars) DL = DL & +0.273d0*dcosd( 217.700d0- G4+ G ) & +2.043d0*dcosd( 343.888d0-2.d0*G4+2.d0*G ) & +1.770d0*dcosd( 200.402d0-2.d0*G4+ G ) & +0.129d0*dcosd( 294.200d0-3.d0*G4+3.d0*G ) & +0.425d0*dcosd( 338.880d0-3.d0*G4+2.d0*G ) & +0.500d0*dcosd( 105.180d0-4.d0*G4+3.d0*G ) & +0.585d0*dcosd( 334.060d0-4.d0*G4+2.d0*G ) & +0.204d0*dcosd( 100.800d0-5.d0*G4+3.d0*G ) & +0.154d0*dcosd( 227.400d0-6.d0*G4+4.d0*G ) & +0.101d0*dcosd( 96.300d0-6.d0*G4+3.d0*G ) & +0.106d0*dcosd( 222.700d0-7.d0*G4+4.d0*G ) !------- ! DL5 (Jupiter) DL = DL & +0.163d0*dcosd( 198.600d0- G5+2.d0*G ) & +7.208d0*dcosd( 179.532d0- G5+ G ) & +2.600d0*dcosd( 263.217d0- G5 ) & +2.731d0*dcosd( 87.145d0-2.d0*G5+2.d0*G ) & +1.610d0*dcosd( 109.493d0-2.d0*G5+ G ) & +0.164d0*dcosd( 170.500d0-3.d0*G5+3.d0*G ) & +0.556d0*dcosd( 82.650d0-3.d0*G5+2.d0*G ) & +0.210d0*dcosd( 98.500d0-3.d0*G5+ G ) !------- ! DL6 (Saturn) DL = DL & +0.419d0*dcosd( 100.580d0- G6+ G ) & +0.320d0*dcosd( 269.460d0- G6 ) & +0.108d0*dcosd( 290.600d0-2.d0*G6+2.d0*G ) & +0.112d0*dcosd( 293.600d0-2.d0*G6+ G ) !------- ! DLM (Moon) DL = DL & +6.454d0*dsind( D ) & +0.177d0*dsind( D+A ) & -0.424d0*dsind( D-A ) & +0.172d0*dsind( D-G ) !------- ! DR2 (Venus) R1 = & 2359.d0*dcosd( 209.080d0+ G2- G ) & + 160.d0*dcosd( 58.400d0+2.d0*G2- G ) & +6842.d0*dcosd( 58.318d0+2.d0*G2-2.d0*G ) & + 869.d0*dcosd( 226.700d0+2.d0*G2-3.d0*G ) & +1045.d0*dcosd( 87.570d0+3.d0*G2-3.d0*G ) & +1497.d0*dcosd( 255.250d0+3.d0*G2-4.d0*G ) & + 194.d0*dcosd( 49.500d0+3.d0*G2-5.d0*G ) & + 376.d0*dcosd( 116.280d0+4.d0*G2-4.d0*G ) & + 196.d0*dcosd( 105.200d0+4.d0*G2-5.d0*G ) & + 163.d0*dcosd( 145.400d0+5.d0*G2-5.d0*G ) & + 141.d0*dcosd( 105.400d0+5.d0*G2-7.d0*G ) !------- ! DR4 (Mars) R1 = R1 & + 150.d0*dcosd( 127.700d0- G4+ G ) & +2057.d0*dcosd( 253.828d0-2.d0*G4+2.d0*G ) & + 151.d0*dcosd( 295.000d0-2.d0*G4+ G ) & + 168.d0*dcosd( 203.500d0-3.d0*G4+3.d0*G ) & + 215.d0*dcosd( 249.000d0-3.d0*G4+2.d0*G ) & + 478.d0*dcosd( 15.170d0-4.d0*G4+3.d0*G ) & + 105.d0*dcosd( 65.900d0-4.d0*G4+2.d0*G ) & + 107.d0*dcosd( 324.600d0-5.d0*G4+4.d0*G ) & + 139.d0*dcosd( 137.300d0-6.d0*G4+4.d0*G ) !------- ! DR5 (Jupiter) R1 = R1 & + 208.d0*dcosd( 112.000d0- G5+2.d0*G ) & +7067.d0*dcosd( 89.545d0- G5+ G ) & + 244.d0*dcosd( 338.600d0- G5 ) & + 103.d0*dcosd( 350.500d0-2.d0*G5+3.d0*G ) & +4026.d0*dcosd( 357.108d0-2.d0*G5+2.d0*G ) & +1459.d0*dcosd( 19.467d0-2.d0*G5+ G ) & + 281.d0*dcosd( 81.200d0-3.d0*G5+3.d0*G ) & + 803.d0*dcosd( 352.560d0-3.d0*G5+2.d0*G ) & + 174.d0*dcosd( 8.600d0-3.d0*G5+ G ) & + 113.d0*dcosd( 347.700d0-4.d0*G5+2.d0*G ) !------- ! DR6 (Saturn) R1 = R1 & + 429.d0*dcosd( 10.600d0- G6+ G ) & + 162.d0*dcosd( 200.600d0-2.d0*G6+2.d0*G ) & + 112.d0*dcosd( 203.100d0-2.d0*G6+ G ) !------- ! DRM (Moon) R1 = R1 & +13360.d0*dcosd( D ) & + 370.d0*dcosd( D+A ) & - 1330.d0*dcosd( D-A ) & - 140.d0*dcosd( D+G ) & + 360.d0*dcosd( D-G ) rDis = rDis+R1*1.d-9 rLat = & -0.210d0*dcosd( 151.800d0+3.d0*G2-4.d0*G ) & -0.166d0*dcosd( 265.500d0-2.d0*G5+ G ) & +0.576d0*dsind( U ) endif rLng = DLO+(SLO+DL)/3600.d0 rLng = mod( mod(rLng,360.d0)+360.d0 , 360.d0 ) rDis = 10**rDis return end