program zodi_distfile_make implicit none c c This program makes a six-year-long data file of the Sun-Earth distance c at one-day intervals c character*3 ADOY real*4 rDis,doy,rLng,rLat,dSMEI real*8 dLngSun,dLatSun,dDisSun integer*4 iyr,idoy c open(11,file='G:\GegenscheinXStars\zodi_dist.txt') c do iyr=2003,2009 c open (10,file='F:\EngineeringXStars\ADOY.txt',readonly) c do idoy=1,367 if(idoy.ge.367.and.(iyr.eq.2004.or.iyr.eq.2008))go to 9 if(idoy.ge.366.and.iyr.ne.2004.and.iyr.ne.2008)go to 9 read (10,10)ADOY 10 format(a3) doy=float(idoy) call SunNewcomb(0,iyr,doy,dLngSun,dLatSun,dDisSun) rDis=sngl(dDisSun) rLng=sngl(dLngSun) rLat=0. call ECLIPTIC_EQUATOR(0,iYr,Doy,rLng,rLat) rLng=rLng-280.-360.*(doy-1.)/365.25 c if(rLng.lt.-180.)rLng=rLng+360. dSMEI=acosd(cosd(rLng)*cosd(rLat+8.8)) write(11,11)adoy,rDis,rLng,rLat,dSMEI,iyr-2003 11 format(1x,a3,f9.6,1x,3f10.4,1x,i3) enddo 9 close(10) 8 enddo c close(11) print *, 'done' stop 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*********************************************************************************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 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