program zodiac_modelc 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,A8,A80,A8exp !,cPA,cPA8 real*4 arg1,arg2,arg3,argx,argy,weight 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 c data R,Omega,RRR,REW /5.,70.,6.,6./ !values 16 Nov 2006 data R,Omega,RRR,REW /5.,78.25,6.,0./ !values 14 November 2008: note 0.75*value RRR is used in North, CBZodi value for Omega data A8,A80,A8exp / 8800.,1200.,10./ !values 18 November 2008: PA dependent term; an old name, but now exponential in 2D 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 /2/ !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 c**************************beware the 65 degree cut below, straddled by "c@#$!"**********************************8 data infile /'F:\EngineeringXStars\cNeq2_200Y_DOY_HHMMSS.grd'/ !Bright-stars subtracted data data outfile /'F:\EngineeringXStars\cNeq3_200Y_DOY_HHMMSS.grd'/ !Bright-stars subtracted data c data infile /'F:\EngineeringOrbits\cNeq2_200Y_DOY_HHMMSS.grd'/ !Bright-stars not subtracted data c data outfile /'F:\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,29 !29 engineering-data calibrations are used covering into year 5. c if(n.gt.27.and.ncam.eq.3)go to 99 c ***Note camera 3, n=18 (2005,DOY306) has bad quaternions near G.C. Thought 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 if(ncam.gt.1)gain=1.00+0.01*(float(iyr-2003)+doy/365.25) !gain correction for CCD response loss vs time: note we use 1%/yr here, not 1.6%/yr. if(ncam.eq.1)gain=1.00+0.01*(float(iyr-2003)+doy/365.25) !gain correction 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... c 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 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) aEloP=180.-rEloP c@#$%! if(rEloP.lt.60..and.imax.eq.3600.and.j.gt.300.and.j.le.1500)eqmap(3601-i,j-300)=0. c@#$%! theta1=PA !+R*cosd(rLngSun-Omega) c cPA=sind(theta1) c cPA8=cPA**8 c 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 ****Entering Parameter Hell**** c {"Abandon all understanding, ye who enter here..."} 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 c zod=75. +120.*c !values 15 Nov 2006 c zod=43. +120.*c !values 15 April 2008 zod=65. +120.*c !values 28 April 2008 ****************Icarus Article Values*************** if(c.le.0.)then !antisolar hemisphere zod=zod+154.*c*c+88.*(c**3) !values 18 April 2008; don't forget c**3 is (-) here... else !solar hemisphere c zod=zod-185.*c*c+140.*(c**3) !values 15 Nov zod=zod-185.*c*c+65.*(c**3) !values 18 Nov 2008, reduced to compensate for added PA term below 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 c arg=abs(sind(argLat))/(0.007*(rEloP+40.)) !2 Nov 2006 values arg=abs(sind(argLat))/(0.009*(rEloP+40.)) !28 April 2008 values zod=zod*(10**(-arg)) c c Ad hoc Latitude adjustments c zod=zod+8.*(1.0-cosd(argLat)) !23 June 2008, see also dc offset below c c For EloP < 90 degrees, include another contribution dropping off as cos(rLat). c if(rEloP.lt.90.)then zod=zod+(1./ssq - 1.)*30.*cosd(argLat) !16 Nov 2006 endif 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 zod=zod+ 6.*exp(-rLat*rLat/512.) !gaussian contribution sigma = 16 deg., 19 May 2008, works best for camera 1 c c Here subtract the double-exponential Gegenschein terms, as in our Icarus manuscript c The weight formula here is an empirical match to the observed Gegenschein distribution near the antisolar direction. c It uses 2% of the sine^2 of twice the angle times the antisolar elongation (deg.) to diminish the geometric weighting of the x,y components. c This is ugly, but reproduces the lozenge shape that the Gegenschein has. c This formula gets somewhat crazy for aEloP > 90 degrees, so we reduce Gegenschein to < 0.1 ADU at aEloP > 90 deg., using a Gaussian rolloff c argx=(180.-abs(rLng))/aEloP argy=rLat/aEloP arg1=5.5*exp(-aEloP/4.)+38.5*exp(-aEloP/19.5) !This one for Latitude arg2=5.5*exp(-aEloP/4.)+38.5*exp(-aEloP/27.) !This one for Longitude arg3=argy*argy*arg1+argx*argx*arg2 weight=(argx*argy)**2 weight=1.-0.02*weight*aEloP if(weight.lt.0.)weight=0. !Kill negative values if(aEloP.gt.60.)weight=weight*exp(-((aEloP-60.)**2)/300.) !Roll it down fast for aEloP > 60 degrees zod=zod+weight*arg3 c c Here add a position-angle term which drops off exponentially with distance from the Sun, thus negligible in the antisolar hemisphere. c c zod=zod+A8*cPA8*exp(-rEloP/A8exp) !April 2008. Old. This varies like the 8th power of the sine of the position angle. c c Here use a double exponential, with a 3 degree hyperbola derived from the deviation of PA from 90 degrees. c arg=abs(PA)-90. zod=zod+(A8*exp(-(sqrt(1.+0.111*arg*arg)-1.)/10.)-A80)*exp(-rEloP/A8exp) !18 November 2008 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. c fractNS=s endif if(adhocNS.lt.0.)then !Ad hoc N for dust slab above if(rElop.gt.90.)then RUSE=0.75*RUSE !Use 3/4 as much in the North in antisolar hemisphere else RUSE=(0.75+2.*c)*RUSE !Use 3/4 as much in the North, but more toward Sun. endif 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(rElop.le.90.)RUSE=(1.+2.0*c)*RUSE !More towards Sun in the South. 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 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 c arg=abs(R*adhocEW*sind(rEloP)) !bend it by the R degree tilt... c RUSE=REW*abs(c) !cosine added 9 Nov to prevent E-W asymmetry, ugly... c if(rEloP.gt.90.)then c RUSE=RUSE*s !Drop off the slab brightness in antisolar hemisphere c else c RUSE=RUSE/ssq !Increase slab brightness towards Sun. c 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 c if(adhocEW.lt.0.)then c if(iii.eq.-1)then !West c if(rLat.ge.0.)zslab=zslab-RUSE*adhocEW c if(rLat.lt.0..and.rLng.gt.-180.)zslab=zslab-RUSE*adhocEW*exp(-(180.+rLng)*2./R) !feather NS line using R/2 dropoff c if(rLat.lt.0..and.rLng.lt.0.)zslab=zslab-RUSE*adhocEW*exp(rLng*3./R) !feather NS line using R/3 dropoff c if(rLat.lt.0..and.rLat.gt.-arg)zslab=zslab-RUSE*adhocEW*(arg+rLat)/arg !feather equatorial band c endif c if(iii.eq.+1)then !East c if(rLat.le.0.)zslab=zslab-RUSE*adhocEW c if(rLat.gt.0..and.rLng.lt.-180.)zslab=zslab-RUSE*adhocEW*exp((180.+rLng)*2./R) !feather NS line using R/2 dropoff c if(rLat.gt.0..and.rLng.gt.0.)zslab=zslab-RUSE*adhocEW*exp(-rLng*3./R) !feather NS line using R/3 dropoff c if(rLat.gt.0..and.rLat.lt.arg)zslab=zslab-RUSE*adhocEW*(arg-rLat)/arg !feather equatorial band c endif c else c if(iii.eq.-1)then !West c if(rLat.le.0.)zslab=zslab+RUSE*adhocEW c if(rLat.gt.0..and.rLng.gt.180.)zslab=zslab+RUSE*adhocEW*exp((180.-rLng)*2./R) !feather NS line using R/2 dropoff c if(rLat.gt.0..and.rLng.lt.0.)zslab=zslab+RUSE*adhocEW*exp(rLng*3./R) !feather NS line using R/3 dropoff c if(rLat.gt.0..and.rLat.lt.arg)zslab=zslab+RUSE*adhocEW*(arg-rLat)/arg !feather equatorial band c endif c if(iii.eq.+1)then !East c if(rLat.ge.0.)zslab=zslab+RUSE*adhocEW c if(rLat.lt.0..and.rLng.lt.180.)zslab=zslab+RUSE*adhocEW*exp(-(180.-rLng)*2./R) !feather NS line using R/2 dropoff c if(rLat.lt.0..and.rLng.gt.0.)zslab=zslab+RUSE*adhocEW*exp(-rLng*3./R) !feather NS line using R/3 dropoff c if(rLat.lt.0..and.rLat.gt.-arg)zslab=zslab+RUSE*adhocEW*(arg+rLat)/arg !feather equatorial band c endif c 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 c Nov 2 2006's value is 6 {from the 75 constant above top line for zod, times 10**(-1/(.007*130))}, plus the value below; but see comment just below c 28 April's 2008 value is 9.1 {from the 65 constant above top line for zod, times 10**(-1/(.009*130))}; but see just below c Then we add another 8 for the latitude adjustment, and finally the below number for the comparison, works out OK. c c zod=zod+16. !DC offset, a case can be made that this should be "zod=zod+REW" if(ncam.eq.1)zod=zod+7. !value 7 Nov. 2008 if(ncam.eq.2)zod=zod+7. !value 29 Oct. 2008 if(ncam.eq.3)zod=zod-3. !value 25 Nov. 2008 c c Diminish zodical light a bit near ecliptic and around 90 degrees, this probably needs refinement, especially towards the Sun c This was added to take account of the "dust bands", which were created here by punching down the middle of the ecliptic -- hokey. 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='F:\engineeringxstars\zodiac1.grd') if(n.eq.2)open (10,file='F:\engineeringxstars\zodiac2.grd') if(n.eq.3)open (10,file='F:\engineeringxstars\zodiac3.grd') if(n.eq.4)open (10,file='F:\engineeringxstars\zodiac4.grd') if(n.eq.5)open (10,file='F:\engineeringxstars\zodiac5.grd') if(n.eq.6)open (10,file='F:\engineeringxstars\zodiac6.grd') if(n.eq.7)open (10,file='F:\engineeringxstars\zodiac7.grd') if(n.eq.8)open (10,file='F:\engineeringxstars\zodiac8.grd') if(n.eq.9)open (10,file='F:\engineeringxstars\zodiac9.grd') if(n.eq.10)open (10,file='F:\engineeringxstars\zodiac10.grd') if(n.eq.11)open (10,file='F:\engineeringxstars\zodiac11.grd') if(n.eq.12)open (10,file='F:\engineeringxstars\zodiac12.grd') if(n.eq.13)open (10,file='F:\engineeringxstars\zodiac13.grd') if(n.eq.14)open (10,file='F:\engineeringxstars\zodiac14.grd') if(n.eq.15)open (10,file='F:\engineeringxstars\zodiac15.grd') if(n.eq.16)open (10,file='F:\engineeringxstars\zodiac16.grd') if(n.eq.17)open (10,file='F:\engineeringxstars\zodiac17.grd') if(n.eq.18)open (10,file='F:\engineeringxstars\zodiac18.grd') if(n.eq.19)open (10,file='F:\engineeringxstars\zodiac19.grd') if(n.eq.20)open (10,file='F:\engineeringxstars\zodiac20.grd') if(n.eq.21)open (10,file='F:\engineeringxstars\zodiac21.grd') if(n.eq.22)open (10,file='F:\engineeringxstars\zodiac22.grd') if(n.eq.23)open (10,file='F:\engineeringxstars\zodiac23.grd') if(n.eq.24)open (10,file='F:\engineeringxstars\zodiac24.grd') if(n.eq.25)open (10,file='F:\engineeringxstars\zodiac25.grd') if(n.eq.26)open (10,file='F:\engineeringxstars\zodiac26.grd') if(n.eq.27)open (10,file='F:\engineeringxstars\zodiac27.grd') if(n.eq.28)open (10,file='F:\engineeringxstars\zodiac28.grd') if(n.eq.29)open (10,file='F:\engineeringxstars\zodiac29.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='F:\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(29),ihh1(29),imm1(29),iss1(29),idoy2(29),ihh2(29),imm2(29),iss2(29),idoy3(29),ihh3(29),imm3(29),iss3(29) 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, 06,125,184/ 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, 00, 01, 05/ 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, 56, 13, 44/ 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, 51, 37, 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, 07,127,186/ 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, 00, 00, 05/ 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, 39, 38, 08/ 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, 04, 02, 48/ 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, 08/ 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, 17, 00/ 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, 57, 21/ 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, 47, 18/ 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