c..this files contains several subroutines. c..organized roughly in functional groups, they are: c..for converting julian dates and calander dates: c..routine juldat2 converts calander dates to julian dates c..routine caldat2 converts julian dates to calander dates c..for converting between universal time and dynamical time: c..routine deltat gets the difference between ut and dt c..routine polint does polynomial interpolation c..for playing with dates and times: c..function mod24 reduces a time argument to the 0.0 to 24.0 interval c..routine dayweek returns the day of the week c..routine numday returns the day number of the year given the month and day c..routine inumday returns the month and day given the day number of the year c..routine daysave computes the dates of daylight savings time c..for converting orbital elements to position and velocity vectors: c..routine kep_to_state converts orbital elements to the state vector c..routine ellip used by kep_to_state for elliptical elements c..function kep_ellip contains keplers equation for elliptical orbits c..routine hyper used by kep_to_state for hyperbolic elements c..function kep_hyper contains keplers equation for hyperbolic orbits c..routine parab used by kep_to_state for parabolic elements c..routine scaled_stumpff evaluates the stumpff functions c..for converting position and velocity vectors to orbital elements: c..routine state_to_kep converts a state vector to orbital elements c..for converting state vectors to orbital plane state vectors: c..routine state_to_orb converts a state vector to orbital vectors c..routine ludcmp does the lu decomposition of a matrix c..routine lubksb does the backsubstitution step in solving linear systems c..for vector coordinate transforms: c..routine rotax is the basic rotation matrix about the x axis c..routine rotay is the basic rotation matrix about the y axis c..routine rotaz is the basic rotation matrix about the z axis c..routine requ_to_recl converts equatorial to ecliptic state vectors c..routine recl_to_requ converts ecliptic to equatorial state vectors c..routine requ_to_rhor converts equatorial to horizon state vectors c..routine rhor_to_requ converts horizon to equatorial state vectors c..routine dh2h converts decimal degrees to hours, minutes, and seconds c..for oblate spheroids: c..for finding the geodetic latitude from the geocentric latitude c..for finding the roots of an equation: c..function zbrent is a general one-dimensional root finder, brent's method c..routine zbrac brackets a one-dimensional root c..for finding the minima of an equation: c..function mnbrent is a general one-dimensional minima finder, brent's method c..routine mnbrak brackets a one-dimensional minima c..for chopping up files that are read in: c..routine getnam is a basic string parser c..function value converts a number in a character string to a real number c..for exact solutions to the two body problem: c..routine twobody solves the two body problem c..routine keprot containes the universal keplers equation c..for orbit integrations: c..routine derv is the right hand side for orbit integrations c..routine cometint is a general ode integration driver c..routine bsstep is a burlish-stoer stepper for odes c..routine mmid kdoes one modified midpoint step for bsstep c..routine pzextr extrapolates a vector of functions to zero c..the next two routines are alternatives to the juldat and caldat routines c..in the novas package, since the novas routines do not take care of the c..switch from julian to gregorian calanders on 15oct1582. hence, c..the novas routines are of no use with the long jpl ephemris de406. c..input calendar dates 1582-oct-15 and after are taken to be expressed c..in the gregorian calendar system. prior dates are assumed to be in the c..julian calendar system. c..historically, not all regions switched calendars at the same time c..(or even in the same century). thus, the user must be aware of c..which calendar was in effect for a particular historical record. c..it should not be assumed this system's calendar automatically c..correlates with a date from an arbitrary historical document. c..here is the progression near the calendar switch point: c c calendar type calendar date julian day number c ------------- ------------- ----------------- c julian 1582-oct-03 2299158.5 c julian 1582-oct-04 2299159.5 ---> c (skipped) "1582-oct-05" | c (skipped) "1582-oct-06" | c (skipped) "1582-oct-07" | c (skipped) "1582-oct-08" | c (skipped) "1582-oct-09" | c (skipped) "1582-oct-10" | c (skipped) "1582-oct-11" | c (skipped) "1582-oct-12" | c (skipped) "1582-oct-13" | c (skipped) "1582-oct-14" | c gregorian 1582-oct-15 2299160.5 <--- c gregorian 1582-oct-16 2299161.5 c gregorian 1582-oct-17 2299162.5 c..in this system there are zero and negative years. c..the progression is as follows: c c julian day number labeling-convention c (jan 1 00:00) bc/ad arithmetical c ----------------- ----- ------------ c 1720327.5 3bc -2 c 1720692.5 2bc -1 c 1721057.5 1bc 0 c 1721423.5 1ad 1 c 1721788.5 2ad 2 subroutine juldat2(iy,im,id,rh,tjd) implicit none save c..converts a calander date to a julian date. c..input time value can be in any ut-like time scale (utc, ut1, tt, etc.) c..and the output julian date will have same basis. c..the astronomical calendar is used. thus the year before 1 ad is 0, c..the year before that is 1 bc. the change to the gregorian calander c..on oct 15, 1582 is accounted for. c..input: c..iy = integer year c..im = integer month c..id = integer day c..rh = number of hours past midnight c..output: c..tjd = julian date c..declare integer im,iy,id,igreg,i,jd parameter (igreg = 15+31*(10+12*1582)) double precision tjd,xy,xm,xa,xb,rh if (im .gt. 2) then xy = iy xm = im + 1 else xy = iy - 1 xm = im + 13 end if i = id + 31 * (im + 12 * iy) if (i .ge. igreg) then xa = int(0.01d0 * xy) xb = 2.0d0 - xa + int(0.25d0 * xa) else xb = 0.0d0 end if jd = int(365.25d0*(xy + 4716.0d0)) + int(30.6001d0*xm) + id tjd = dble(jd) + xb - 1524.5d0 + rh/24.d0 return end subroutine caldat2(tjd,iy,im,id,rh) implicit none save c..converts a julian date to a calander date. c..input time value can be in any ut-like time scale (utc, ut1, tt, etc.) c..and the output calander date will have same basis. c..the astronomical calendar is used. thus the year before 1 ad is 0, c..the year before that is 1 bc. the change to the gregorian calander c..on oct 15, 1582 is accounted for. c..input: c..tjd = julian date c..iy = integer year c..im = integer month c..id = integer day c..rh = number of hours past midnight c..output: c..iy = integer year c..im = integer month c..id = integer day c..rh = number of hours past midnight c..declare integer id,im,iy,igreg parameter (igreg = 2299161) double precision tjd,rh,x1,z,f,x2,xa,xb,xc,xd,xe,rd,c1,c2,c3 parameter (c1 = 1.0d0/36524.25d0, 1 c2 = 1.0d0/365.25d0, 2 c3 = 1.0d0/30.6001d0) x1 = tjd + 0.5d0 z = int(x1) f = x1 - z if (x1 .ge. igreg) then x2 = int((x1-1867216.25d0)*c1) xa = z + 1.0d0 +x2 - int(0.25d0 * x2) else xa = z end if xb = xa + 1524.0d0 xc = int((xb - 122.1d0)*c2) xd = int(365.25d0*xc) xe = int((xb-xd)*c3) rd = xb - xd - int(30.6001d0*xe) + f id = rd rh = 24.0d0*(rd - dble(id)) im = xe - 1 if (im .gt. 12) im = im - 12 iy = xc - 4715 if (im .gt. 2) iy = iy - 1 return end subroutine deltat(tjd,tjde,secdif) implicit none save c..this routine computes the difference between c..universal time and dynamical time. c..input: c..tjd = ut julian day number c..output: c..tjde = tjd + secdif c..secdif = ut - et in seconds c..declare the pass double precision tjd,tjde,secdif c..for storing the table integer tabstart,tabend,tabsize parameter (tabstart = 1620, 1 tabend = 2005, 2 tabsize = tabend - tabstart + 1) double precision dt(tabsize),years(tabsize) c..local variables integer i,j,iy,im,id,iflag,mp,iat parameter (mp = 2) double precision hr,b,xx,dy c..parameter mp above determines the order of the interpolatant c..when interpolating in the table. mp=2=linear mp=3=quadratic and so on. c..don't be too greedy about the order of the interpolant since c..for many years the data is flat, and trying to fit anything other c..that a line will produce unwanted oscillations. thus i choose mp=2. c..these tables of observed and extrapolated data are c..are from the us naval observatory ftp://maia.usno.navy.mil/ser7/ c..years 1620 to 1710 data (dt(i), i=1,90) / 1 124.00d0, 119.00d0, 115.00d0, 110.00d0, 106.00d0, 102.00d0, 2 98.00d0, 95.00d0, 91.00d0, 88.00d0, 85.00d0, 82.00d0, 3 79.00d0, 77.00d0, 74.00d0, 72.00d0, 70.00d0, 67.00d0, 4 65.00d0, 63.00d0, 62.00d0, 60.00d0, 58.00d0, 57.00d0, 5 55.00d0, 54.00d0, 53.00d0, 51.00d0, 50.00d0, 49.00d0, 6 48.00d0, 47.00d0, 46.00d0, 45.00d0, 44.00d0, 43.00d0, 7 42.00d0, 41.00d0, 40.00d0, 38.00d0, 37.00d0, 36.00d0, 8 35.00d0, 34.00d0, 33.00d0, 32.00d0, 31.00d0, 30.00d0, 9 28.00d0, 27.00d0, 26.00d0, 25.00d0, 24.00d0, 23.00d0, & 22.00d0, 21.00d0, 20.00d0, 19.00d0, 18.00d0, 17.00d0, 1 16.00d0, 15.00d0, 14.00d0, 14.00d0, 13.00d0, 12.00d0, 2 12.00d0, 11.00d0, 11.00d0, 10.00d0, 10.00d0, 10.00d0, 3 9.00d0, 9.00d0, 9.00d0, 9.00d0, 9.00d0, 9.00d0, 4 9.00d0, 9.00d0, 9.00d0, 9.00d0, 9.00d0, 9.00d0, 5 9.00d0, 9.00d0, 9.00d0, 9.00d0, 10.00d0, 10.00d0/ c..years 1711 to 1799 data (dt(i), i=91, 180) / 1 10.00d0, 10.00d0, 10.00d0, 10.00d0, 10.00d0, 10.00d0, 2 10.00d0, 11.00d0, 11.00d0, 11.00d0, 11.00d0, 11.00d0, 3 11.00d0, 11.00d0, 11.00d0, 11.00d0, 11.00d0, 11.00d0, 4 11.00d0, 11.00d0, 11.00d0, 11.00d0, 11.00d0, 11.00d0, 5 12.00d0, 12.00d0, 12.00d0, 12.00d0, 12.00d0, 12.00d0, 6 12.00d0, 12.00d0, 12.00d0, 12.00d0, 13.00d0, 13.00d0, 7 13.00d0, 13.00d0, 13.00d0, 13.00d0, 13.00d0, 14.00d0, 8 14.00d0, 14.00d0, 14.00d0, 14.00d0, 14.00d0, 14.00d0, 9 15.00d0, 15.00d0, 15.00d0, 15.00d0, 15.00d0, 15.00d0, & 15.00d0, 16.00d0, 16.00d0, 16.00d0, 16.00d0, 16.00d0, 1 16.00d0, 16.00d0, 16.00d0, 16.00d0, 16.00d0, 17.00d0, 2 17.00d0, 17.00d0, 17.00d0, 17.00d0, 17.00d0, 17.00d0, 3 17.00d0, 17.00d0, 17.00d0, 17.00d0, 17.00d0, 17.00d0, 4 17.00d0, 17.00d0, 17.00d0, 17.00d0, 16.00d0, 16.00d0, 5 16.00d0, 16.00d0, 15.00d0, 15.00d0, 14.00d0, 14.00d0/ c..years 1800 to 1890 data (dt(i), i=181, 270) / 1 13.70d0, 13.40d0, 13.10d0, 12.90d0, 12.70d0, 12.60d0, 2 12.50d0, 12.50d0, 12.50d0, 12.50d0, 12.50d0, 12.50d0, 3 12.50d0, 12.50d0, 12.50d0, 12.50d0, 12.50d0, 12.40d0, 4 12.30d0, 12.20d0, 12.00d0, 11.70d0, 11.40d0, 11.10d0, 5 10.60d0, 10.20d0, 9.60d0, 9.10d0, 8.60d0, 8.00d0, 6 7.50d0, 7.00d0, 6.60d0, 6.30d0, 6.00d0, 5.80d0, 7 5.70d0, 5.60d0, 5.60d0, 5.60d0, 5.70d0, 5.80d0, 8 5.90d0, 6.10d0, 6.20d0, 6.30d0, 6.50d0, 6.60d0, 9 6.80d0, 6.90d0, 7.10d0, 7.20d0, 7.30d0, 7.40d0, & 7.50d0, 7.60d0, 7.70d0, 7.70d0, 7.80d0, 7.80d0, 1 7.88d0, 7.82d0, 7.54d0, 6.97d0, 6.40d0, 6.02d0, 2 5.41d0, 4.10d0, 2.92d0, 1.82d0, 1.61d0, 0.10d0, 3 -1.02d0, -1.28d0, -2.69d0, -3.24d0, -3.64d0, -4.54d0, 4 -4.71d0, -5.11d0, -5.40d0, -5.42d0, -5.20d0, -5.46d0, 5 -5.46d0, -5.79d0, -5.63d0, -5.64d0, -5.80d0, -5.66d0/ c..years 1890 to 1979 data (dt(i), i=271, 360) / 1 -5.87d0, -6.01d0, -6.19d0, -6.64d0, -6.44d0, -6.47d0, 2 -6.09d0, -5.76d0, -4.66d0, -3.74d0, -2.72d0, -1.54d0, 3 -0.2d0, 1.24d0, 2.64d0, 3.86d0, 5.37d0, 6.14d0, 4 7.75d0, 9.13d0, 10.46d0, 11.53d0, 13.36d0, 14.65d0, 5 16.01d0, 17.20d0, 18.24d0, 19.06d0, 20.25d0, 20.95d0, 6 21.16d0, 22.25d0, 22.41d0, 23.03d0, 23.49d0, 23.62d0, 7 23.86d0, 24.49d0, 24.34d0, 24.08d0, 24.02d0, 24.00d0, 8 23.87d0, 23.95d0, 23.86d0, 23.93d0, 23.73d0, 23.92d0, 9 23.96d0, 24.02d0, 24.33d0, 24.83d0, 25.30d0, 25.70d0, & 26.24d0, 26.77d0, 27.28d0, 27.78d0, 28.25d0, 28.71d0, 1 29.15d0, 29.57d0, 29.97d0, 30.36d0, 30.72d0, 31.07d0, 2 31.35d0, 31.68d0, 32.18d0, 32.68d0, 33.15d0, 33.59d0, 3 34.00d0, 34.47d0, 35.03d0, 35.73d0, 36.54d0, 37.43d0, 4 38.29d0, 39.20d0, 40.18d0, 41.17d0, 42.23d0, 43.37d0, 5 44.49d0, 45.48d0, 46.46d0, 47.52d0, 48.53d0, 49.59d0/ c..Using data from the Astronomical Almanac 2004 page K9 (publ 2002) and the c..most recent data from IERS, http://maia.usno.navy.mil/bulletin-a.html, c..david walden suggests these lines: c..years 1980 to 2005 data (dt(i), i=361, 386) / 1 50.54d0, 51.38d0, 52.17d0, 52.96d0, 53.79d0, 54.34d0, 2 54.87d0, 55.32d0, 55.82d0, 56.30d0, 56.86d0, 57.57d0, 3 58.31d0, 59.12d0, 59.98d0, 60.78d0, 61.63d0, 62.29d0, 4 62.97d0, 63.47d0, 63.83d0, 64.09d0, 64.30d0, 64.470d0, 5 64.57d0, 64.63d0/ c..first time flag data iflag/0/ c..if this is the first time, initialize c..put the julain date at the first of the year in the array c..years so that we can interpolate/extrapolate directly on c..the given ut julian date if (iflag .eq. 0) then iflag = 1 do i=1,tabsize j = tabstart + (i-1) call juldat2(j,1,1,0.0d0,xx) years(i) = xx end do endif c..convert the given ut julian date to a ut calander date call caldat2(tjd,iy,im,id,hr) c..if we are outside the table on the low end c..use the stephenson and morrison expression 948 to 1600, c..and the borkowski formula for earlier years if (iy .lt. tabstart) then if (iy .ge. 948) then b = 0.01d0 * float(iy - 2000) secdif = b * (b * 23.58d0 + 100.3d0) + 101.6d0 else b = 0.01d0 * float(iy -2000) + 3.75d0 secdif = 35.0d0 * b * b + 40.0d0 end if c..if we are outside the table on the high end c..use a linear extrapolation into the future else if (iy .gt. tabend) then b = float(iy - tabend) secdif = dt(tabsize) + b*(dt(tabsize) - dt(tabsize-1)) c..otherwise we are in the table c..get the table location and interpolate else iat = iy - tabstart + 1 iat = max(1, min(iat - mp/2 + 1, tabsize - mp + 1)) call polint(years(iat),dt(iat),mp,tjd,secdif,dy) c..the astronomical almanac table is corrected by adding the expression c.. -0.000091 (ndot + 26)(year-1955)^2 seconds c..to entries prior to 1955 (page K8), where ndot is the secular tidal c..term in the mean motion of the moon. entries after 1955 are referred c..to atomic time standards and are not affected by errors in lunar c..or planetary theory. a value of ndot = -25.8 arcsec per century squared c..is the value used in jpl's de403 ephemeris, the earlier de200 ephemeris c..used the value -23.8946. note for years below the table (less than 1620) c..the time difference is not adjusted for small improvements in the c..current estimate of ndot because the formulas were derived from c..studies of ancient eclipses and other historical information, whose c..interpretation depends only partly on ndot. c..here we make the ndot correction. if (iy .lt. 1955) then b = float(iy - 1955) secdif = secdif - 0.000091d0 * (-25.8d0 + 26.0d0)*b*b end if end if c..add the difference to the ut julian date to get the dynamical julian date tjde = tjd + secdif/86400.0d0 return end subroutine polint(xa,ya,n,x,y,dy) implicit none save c..given arrays xa and ya of length n and a value x, this routine returns a c..value y and an error estimate dy. if p(x) is the polynomial of degree n-1 c..such that ya = p(xa) ya then the returned value is y = p(x) c..input: c..xa(1:n) = array of x values c..ya(1:n) = array of y values c..n = order of interpolant, 2=linear, 3=quadratic ... c..x = x value where interpolation is desired c..output: c..y = interpolated y value c..dy = error esimate c..declare integer n,nmax,ns,i,m parameter (nmax=10) double precision xa(n),ya(n),x,y,dy,c(nmax),d(nmax),dif,dift, 1 ho,hp,w,den c..find the index ns of the closest table entry; initialize the c and d tables ns = 1 dif = abs(x - xa(1)) do i=1,n dift = abs(x - xa(i)) if (dift .lt. dif) then ns = i dif = dift end if c(i) = ya(i) d(i) = ya(i) enddo c..first guess for y y = ya(ns) c..for each column of the table, loop over the c's and d's and update them ns = ns - 1 do m=1,n-1 do i=1,n-m ho = xa(i) - x hp = xa(i+m) - x w = c(i+1) - d(i) den = ho - hp if (den .eq. 0.0) stop ' 2 xa entries are the same in polint' den = w/den d(i) = hp * den c(i) = ho * den enddo c..after each column is completed, decide which correction c or d, to add c..to the accumulating value of y, that is, which path to take in the table c..by forking up or down. ns is updated as we go to keep track of where we c..are. the last dy added is the error indicator. if (2*ns .lt. n-m) then dy = c(ns+1) else dy = d(ns) ns = ns - 1 end if y = y + dy enddo return end double precision function mod24(x) implicit none save c..reduce the time argument x to the 0.0 to 24.0 interval c..input: c..x = number of hours c..output: c..mod24 = x between 0 and 24 c..declare and go double precision x if (x .ge. 0.0) then mod24 = x - 24.0d0*int(x/24.0d0) else mod24 = x - 24.0d0*int(x/24.0d0) + 24.0d0 endif return end subroutine dayweek(iy,im,id,cday) implicit none save c..given the month day and year this routine returns the day of the week. c..input: c..iy = integer year c..im = integer month c..id = integer day c..output: c..cdat = day of week as a 3 letter character string c..declare character*3 cday,week(0:6) integer iy,im,id double precision tjd,a,b,rh,seventh parameter (seventh = 1.0d0/7.0d0, rh=0.0d0) data week/'sun','mon','tue','wed','thu','fri','sat'/ call juldat2(iy,im,id,rh,tjd) a = (tjd + 1.5d0) * seventh b = int(a) cday = week(nint(7.0d0*abs(a-b))) return end subroutine numday(iyear,imonth,iday,nday) implicit none save c..given the month day and year, n is returned as number day of the year. c..input: c..iy = integer year c..im = interger month c..id = integer day c..output: c..nday = integer day number of year c..declare integer iyear,imonth,iday,nday,k double precision a,b,c1,c2 parameter (c1 = 275.0d0/9.0d0, c2 = 1.0d0/12.0d0) c..figure out if it was a leap year (k=1) if (mod(iyear,4) .eq. 0) then k = 1 else k = 2 end if c..gregorian modification to the leap year scheme if (iyear .gt. 1582 .and. mod(iyear,100).eq.0 .and. 1 mod(iyear,400).ne.0) k = 2 c..the number day of the year a = float(imonth)*c1 b = (float(imonth) + 9.0d0)*c2 nday = int(a) - k*int(b) + iday - 30 return end subroutine inumday(iyear,imonth,iday,nday) implicit none save c..given the number day n and the year, this routine returns the month and day c..input: c..nday = integer day number of year c..iy = integer year c..output: c..im = interger month c..id = integer day c..declare integer iyear,imonth,iday,nday,k double precision day,c1,c2,c3 parameter (c1 = 1.0d0/9.0d0, c2 = 1.0d0/12.0d0, 1 c3 = 1.0d0/275.0d0) c..figure out if it was a leap year (k=1) if (mod(iyear,4) .eq. 0) then k = 1 else k = 2 end if c..gregorian modification to the leap year if (iyear .gt. 1582 .and. mod(iyear,100).eq.0 .and. 1 mod(iyear,400).ne.0) k = 2 c..the month and day if (nday .lt. 32) then imonth = 1 else imonth = int(9.0d0*(k+nday)*c2 + 0.98d0) end if iday = nday - int(275.0d0*imonth*c1) + 1 k*int((imonth+9.0d0)*c2) + 30.0d0 return end subroutine daysave(iyear,imonth,iday, 1 ismonth,isday,ifmonth,ifday,saving) implicit none save c..this routine returns the months and days of when daylight savings applies. c..input: c..iyear = integer year c..imonth = integer month c..iday = integer day c..output: c..ismonth = integer month when daylight savings turns on in spring c..isday = integer day when daylight savings turns on in spring c..ifmonth = integer month when daylight savings turns off in fall c..ifday = integer day when daylight savings turns off in fall c..saving = value of the daylight savings correction c..declare the pass integer iyear,imonth,iday,ismonth,isday,ifmonth,ifday double precision saving c..local variables character*3 cday integer i,iy double precision tjd,tjdx,tjdspring,tjdfall,rh,x c..get the julian day for jan 1 of the year and its present value call juldat2(iyear,1,1,0.0d0,tjd) call juldat2(iyear,imonth,iday,0.0d0,tjdx) c..find the first sunday in april, start with a reasonable guess tjdspring = tjd + 85.0d0 do i=1,30 tjdspring = tjdspring + 1.0d0 call caldat2(tjdspring,iy,ismonth,isday,rh) call dayweek(iyear,ismonth,isday,cday) if (ismonth.eq.4 .and. isday.ge.1.0 .and. cday.eq.'sun') goto 20 enddo stop 'error: did not find april first in routine daysave' 20 continue c..find the last sunday in october, tstart with a reasonable guess tjdfall = tjd + 265.0d0 do i = 1,50 tjdfall = tjdfall + 1.0d0 call caldat2(tjdfall,iy,ifmonth,ifday,rh) call dayweek(iyear,ifmonth,ifday,cday) if (ifmonth.eq.10 .and. ifday.ge.1.0 .and. cday.eq.'sun') then x = tjdfall else if (ifmonth .eq. 11) then goto 30 end if enddo stop 'error: did not find last october first in routine daysave' 30 continue call caldat2(x,iy,ifmonth,ifday,rh) c..compute the daylight savings correction saving = 0.0d0 if (tjdx .ge. tjdspring .and. tjdx .lt. tjdfall) saving = 1.0d0 return end subroutine kep_to_state(gm, 1 tjdp,tjd,xmanom,qper,xecen,lan,inc,aop, 2 pos,vel) implicit none save c..given the orbital elements, this routine returns the state vector c..for either elliptical, hyperbolic, or parabolic orbits c..reference plane is determined by the orbital elements c..input: c..gm = g*(m1 + m2) = kgrav**2 (1 + m2/m1) c.. where kgrav= kgauss = 0.01720209895d0 for orbits about the sun c.. kgrav = kgauss*(mb/msun) for orbits about mb c.. kgrav= kearth = 0.07436680d0 for orbits about the earth c..qper = perihelion distance in au c..xecen = orbital eccentricity c..lan = longitude of the ascending node in decimal degrees c..inc = inclination in degrees in decimal degrees c..aop = argument of perihelion in decimal degrees c..xmanon = mean anomaly in decimal degrees c..tjdp = julian day number of perhelion c..tjd = julian day number where the state vector os to be computed c..output: c..pos = coordinates of the ecliptic position vector in au at tjd c..vel = coordinates of the ecliptic velocity vector in au/day at tjd c..declare the pass double precision gm,tjdp,tjd,xmanom,qper,xecen,lan,inc,aop, 1 pos(3),vel(3) c..math constants double precision pi,a2rad,rad2a parameter (pi = 3.1415926535897932384d0, 1 a2rad = pi/180.0d0, 2 rad2a = 180.0d0/pi) c..for silly compilers that don't define the trig functions in degrees external cosd,sind,tand,acosd,asind,atan2d double precision cosd,sind,tand,acosd,asind,atan2d c..local variables double precision delta,semi_maj, 1 xorb,yorb,vxorb,vyorb, 2 p(3),q(3),r(3),c1,s1,c2,s2,c3,s3, 3 xmlimit,eclimit parameter (xmlimit = 5.0d0, 1 eclimit = 0.1d0) c..popular format statements 222 format(1x,1p6e16.8) c..test value for the eccentricity delta = abs(1.0d0 - xecen) c..parabolic orbit if (abs(xmanom) .lt. xmlimit .and. delta .lt. eclimit) then c write(6,*) 'parabolic orbit' if (tjdp .ne. 0.0) then call parab(gm,tjdp,tjd,qper,xecen,xorb,yorb,vxorb,vyorb) else stop 'tjdp is zero and tried a parablic orbit' end if c..elliptical orbit else if (xecen .lt. 1.0d0) then c write(6,*) 'elliptical orbit' semi_maj = qper/delta call ellip(gm,xmanom,semi_maj,xecen,xorb,yorb,vxorb,vyorb) c..hyperbolic orbit else c write(6,*) 'hyperbolic orbit' if (tjdp .ne. 0.0) then semi_maj = qper/delta call hyper(gm,tjdp,tjd,semi_maj,xecen,xorb,yorb,vxorb,vyorb) else stop 'tjdp is zero and tried a hyperbolic orbit' end if end if c..we should now have the orbital state vector xorb,yorb,vxorb,vyorb c..now get the gaussian vectors p, q, and r c1 = cosd(aop) s1 = sind(aop) c2 = cosd(inc) s2 = sind(inc) c3 = cosd(lan) s3 = sind(lan) p(1) = c1*c3 - s1*c2*s3 p(2) = c1*s3 + s1*c2*c3 p(3) = s1*s2 q(1) = -s1*c3 - c1*c2*s3 q(2) = -s1*s3 + c1*c2*c3 q(3) = c1*s2 r(1) = s2*s3 r(2) = -s2*c3 r(3) = c2 c..expand orbital state to the heliocentric ecliptic state vector pos(1) = p(1)*xorb + q(1)*yorb pos(2) = p(2)*xorb + q(2)*yorb pos(3) = p(3)*xorb + q(3)*yorb vel(1) = p(1)*vxorb + q(1)*vyorb vel(2) = p(2)*vxorb + q(2)*vyorb vel(3) = p(3)*vxorb + q(3)*vyorb return end subroutine ellip(gm,xmanom,semi_maj,xecen,xorb,yorb,vxorb,vyorb) implicit none save c..this routine computes the orbital state vector for elliptical orbits c..input: c..gm = gravitational constant = kgrav**2 * (1.0d0 + m2/m1) c..xmanom = mean anomoly in decimal degrees c..semi_maj = semi-major axis in au c..xecen = eccentricity c..output: c..xorb = orbital x coordinate in au c..yorb = orbital y coordinate in au c..vxorb = orbital x speed in au/day c..vyorb = orbital y spped in au/day c..declare the pass double precision gm,xmanom,semi_maj,xecen,xorb,yorb,vxorb,vyorb c..math constants double precision pi,a2rad,rad2a parameter (pi = 3.1415926535897932384d0, 1 a2rad = pi/180.0d0, 2 rad2a = 180.0d0/pi) c..for communicating the solution of keplers equation double precision manom_kep,ecen_kep common /eqkep/ manom_kep,ecen_kep c..for solving keplers equation external kep_ellip logical succes integer niter double precision low,high,kep_ellip,tol,zbrent parameter (tol = 1.0d-10) c..local variables double precision eanom,xx,rho,fac,ceanom,seanom c..load the common block for keplers equation manom_kep = xmanom * a2rad ecen_kep = xecen c..bracket the eccentric anomaly in keplers equation if (manom_kep .eq. 0.0) then low = -1.0d0 high = 1.0d0 else low = 0.5d0 * manom_kep high = 2.0d0 * manom_kep end if call zbrac(kep_ellip,low,high,succes) if (.not. succes) then write(6,*) 'could not bracket in routine ellip' write(6,*) 'low =',low,' high=',high end if c..solve keplers equations for the eccentric anomoly in radians eanom = zbrent(kep_ellip,low,high,tol,niter) c..useful factors xx = sqrt(gm/semi_maj) ceanom = cos(eanom) seanom = sin(eanom) rho = 1.0d0 - xecen*ceanom fac = sqrt((1.0d0-xecen) * (1.0d0+xecen)) c..positions and velocities in the orbital plane xorb = semi_maj * (ceanom - xecen) yorb = semi_maj * fac * seanom vxorb = -xx * seanom/rho vyorb = xx*fac*ceanom/rho return end double precision function kep_ellip(eanom) implicit none save c..keplers equation for elliptical orbits. c..input: c..eanom = eccentric anomaly in radians c..output: c..kep_ellip = keplers equation as f(x) = 0 c..declare double precision eanom c..common block communication double precision manom_kep,ecen_kep common /eqkep/ manom_kep,ecen_kep kep_ellip = eanom - ecen_kep * sin(eanom) - manom_kep return end subroutine hyper(gm,tjdp,tjd,semi_maj,xecen,xorb,yorb,vxorb,vyorb) implicit none save c..this routine computes the orbital state vector for hyperbolic orbits c..input: c..gm = gravitational constant = kgrav**2 * (1.0d0 + m2/m1) c..tjdp = julian day number of perhelion passage c..tjd = julian day number for the date of interest c..semi_maj = semi-major axis in au c..xecen = eccentricity c..output: c..xorb = orbital x coordinate in au c..yorb = orbital y coordinate in au c..vxorb = orbital x speed in au/day c..vyorb = orbital y spped in au/day c..declare the pass double precision gm,tjdp,tjd,semi_maj,xecen,xorb,yorb,vxorb,vyorb c..for communicating the solution of keplers equation double precision manom_kep,ecen_kep common /eqkep/ manom_kep,ecen_kep c..for solving keplers equation external kep_hyper logical succes integer niter double precision low,high,kep_hyper,tol,zbrent parameter (tol = 1.0d-10) c..local variables double precision a,xx,xmanom,eanom,ceanom,seanom,fac,rho c..get the mean anomoly in radians a = abs(semi_maj) xx = sqrt(gm/a) xmanom = xx/a * (tjd - tjdp) c..load the common block for keplers equation manom_kep = xmanom ecen_kep = xecen c..bracket the eccentric anomaly in keplers equation if (manom_kep .ge. 0.0) then fac = log(1.8 + 2.0d0*manom_kep/2.7182818d0) else fac = log(1.8 - 2.0d0*manom_kep/2.7182818d0) end if low = 0.5d0 * fac high = 2.0d0 * fac call zbrac(kep_hyper,low,high,succes) if (.not. succes) then write(6,*) 'could not bracket in routine hyper' write(6,*) 'low =',low,' high=',high end if c..solve keplers equations for the eccentric anomoly in radians eanom = zbrent(kep_hyper,low,high,tol,niter) c..useful factors ceanom = cosh(eanom) seanom = sinh(eanom) rho = xecen*ceanom - 1.0d0 fac = sqrt((xecen - 1.0d0) * (xecen + 1.0d0)) c..positions and velocities in the orbital plane xorb = a * (xecen - ceanom) yorb = a * fac * seanom vxorb = -xx * seanom/rho vyorb = xx * fac * ceanom/rho return end double precision function kep_hyper(eanom) implicit none save c..keplers equation for hyperbolic orbits. c..ecentric anomoly eanom, in radians, is iterated on by a root finder. c..declare double precision eanom c..common block communication double precision manom_kep,ecen_kep common /eqkep/ manom_kep,ecen_kep kep_hyper = ecen_kep * sinh(eanom) - eanom - manom_kep return end subroutine parab(gm,tjdp,tjd,qper,xecen,xorb,yorb,vxorb,vyorb) implicit none save c..this routine computes the orbital state vector for parabolic orbits c..input: c..gm = gravitational constant = kgrav**2 * (1.0d0 + m2/m1) c..tjdp = julian day number of perhelion passage c..tjd = julian day number for the date of interest c..qper = perhelion distance in au c..xecen = eccentricity c..output: c..xorb = orbital x coordinate in au c..yorb = orbital y coordinate in au c..vxorb = orbital x speed in au/day c..vyorb = orbital y spped in au/day c..declare the pass double precision gm,tjdp,tjd,qper,xecen,xorb,yorb,vxorb,vyorb c..local variables integer i,maxit parameter (maxit = 40) double precision xx,tau,e2,fac,e20,a,u,u2,c0,c1,c2,c3,r, 1 tol,third parameter (tol = 1.0d-10, 1 third = 1.0d0/3.0d0) c..initialize xx = sqrt(gm / (qper*(1.0d0 + xecen)) ) tau = sqrt(gm) * (tjd - tjdp) e2 = 0.0d0 fac = 0.5d0 * xecen c..iteration on the barker/kepler equation for u do i=1,maxit e20 = e2 a = 1.5d0*sqrt(fac/(qper*qper*qper))*tau a = (sqrt(a*a + 1.0d0) + a)**third u = a - 1.0d0/a u2 = u * u e2 = u2*(1.0d0 - xecen)/fac call scaled_stumpff(e2,c0,c1,c2,c3) c call stumpff(e2,c1,c2,c3) fac = 3.0d0*xecen*c3 if (abs(e2 - e20) .le. tol) goto 20 enddo write(6,*) 'could not get the root in routine parab' stop 'error getting barkers equation solved' 20 continue c..positions and velocities in the orbital plane r = qper * (1.0d0 + u2*c2*xecen/fac) xorb = qper * (1.0d0 - u2*c2/fac) yorb = qper * sqrt((1.0d0 + xecen)/fac) * u * c1 vxorb = -xx * yorb/r vyorb = xx * (xorb/r + xecen) return end subroutine scaled_stumpff(x,c0,c1,c2,c3) implicit none save c..given a positive, zero, or negative x, c..this routine gets the four stumpff functions c0 c1 c2 and c3. c..c0 = cos(x) c..c1 = sin(x)/x c..c2 = (1 - cos(x))/x**2 c..c3 = (x - sin(x))/x**3 c..its done this way to minimize small cancellations that occur c..if evalutated directly when x is small. c..declare integer n double precision x,c0,c1,c2,c3,xx,xmax parameter (xmax = 0.1d0) c..scale the argument if need be, by the half angle formulas n = 0 xx = x 10 if ( abs(xx) .gt. xmax) then n = n + 1 xx = 0.25d0 * xx go to 10 end if c..the first four, which for the scaled argument c..guarantees double precision accuracy for c2 and c3. c..c0 and c1 are had by the recursion relation. c2 = (1.0d0 - xx*(1.0d0 - xx*(1.0d0 - xx*(1.0d0 - xx*(1.0d0 - xx* 1 (1.0d0 - xx/182.0d0)/132.0d0)/90.0d0)/56.0d0)/30.0d0)/ 2 12.0d0)/2.0d0 c3 = (1.0d0 - xx*(1.0d0 - xx*(1.0d0 - xx*(1.0d0 - xx*(1.0d0 - xx* 1 (1.0d0 - xx/210.0d0)/156.0d0)/110.0d0)/72.0d0)/42.0d0)/ 2 20.0d0)/6.0d0 c1 = 1.0d0 - xx * c3 c0 = 1.0d0 - xx * c2 c..unscale, if need be, using the double angle formulas 20 if (n .eq. 0) return n = n - 1 c3 = 0.25d0 * (c2 + c0*c3) c2 = 0.5d0 * c1 * c1 c1 = c0 * c1 c0 = 2.0d0 * c0 * c0 - 1.0d0 go to 20 end subroutine stumpff(x,c1,c2,c3) implicit none save c..calculation of stumpff's functions for parabolic orbits c..c1 = sin(x)/x c..c2 = (1 - cos(x))/x**2 c..c3 = (x - sin(x))/x**3 c..its done this way to minimize small cancellations that occur c..if evalutated directly when x is small. c..declare integer i,imax parameter (imax=50) double precision x,c1,c2,c3,add,xn,eps parameter (eps = 1.0d-14) c1 = 0.0d0 c2 = 0.0d0 c3 = 0.0d0 add = 1.0d0 xn = 1.0d0 do i=1,imax c1 = c1 + add add = add/(2.0d0 * xn) c2 = c2 + add add = add/(2.0d0 * xn + 1.0d0) c3 = c3 + add add = -x*add xn = xn + 1.0d0 if (abs(add) .le. eps) return end do stop 'stumpffs functions did not converge' end subroutine state_to_kep(gm,tjd,pos,vel, 1 semi_maj,qper,xecen,lan,inc,aop,xmanom,period,tjdp) implicit none save c..given state vector, this routine returns the orbital elements for c..elliptical, hyperbolic or parabolic orbits c..input: c..gm = gravitational constant = g*(m1 + m2) = kgrav**2 (1 + m2/m1) c.. where kgrav= kgauss = 0.01720209895d0 for orbits about the sun c.. kgrav= kearth = 0.07436680d0 for orbits about the earth c..tjd = julian day number of the position and velocity vectors c..pos = coordinates of the ecliptic position vector in au c..vel = coordinates of the ecliptic velocity vector in au/day c..output: c..semi_maj = semi-major axis in au c..qper = perihelion distance in au c..xecen = orbital eccentricity c..lan = longitude of the ascending node in decimal degrees c..inc = inclination in degrees in decimal degrees c..aop = argument of perihelion in decimal degrees c..xmanom = mean anomaly in decimal degrees c..period = period in days c..tjdp = julian day number of perhelion c..declare the pass double precision gm,tjd,pos(3),vel(3), 1 semi_maj,qper,xecen,lan,inc,aop,xmanom, 2 period,tjdp c..math constants double precision pi,a2rad,rad2a,twopi parameter (pi = 3.1415926535897932384d0, 1 a2rad = pi/180.0d0, 2 rad2a = 180.0d0/pi, 3 twopi = 2.0d0 * pi) c..for silly compilers that don't define the trig functions in degrees external cosd,sind,tand,acosd,asind,atan2d double precision cosd,sind,tand,acosd,asind,atan2d c..local variables double precision x,y,z,vx,vy,vz,r,v2,rdotv, 1 hx,hy,hz,h,ex,ey,ez,e2,delta,wx,wy,wz,w, 2 c1,s1,eanom,nu,a,fac1,fac2,elimit parameter (elimit = 1.0d-3) c..popular format statements 222 format(1x,1p6e16.8) c..initailize xmanom = 0.0d0 semi_maj = 0.0d0 xecen = 0.0d0 lan = 0.0d0 inc = 0.0d0 aop = 0.0d0 c..unload for easier coding x = pos(1) y = pos(2) z = pos(3) vx = vel(1) vy = vel(2) vz = vel(3) c..distance, velocity squared, and rdotv r = sqrt(x**2 + y**2 + z**2) v2 = vx**2 + vy**2 + vz**2 rdotv = x*vx + y*vy + z*vz c..angular momentum = r cross v components hx = y*vz - z*vy hy = z*vx - x*vz hz = x*vy - y*vx h = sqrt(hx**2 + hy**2 + hz**2) c..eccentricity components, the eccentricity and its difference from unity fac1 = v2/gm - 1.0d0/r fac2 = rdotv/gm ex = fac1*x - fac2*vx ey = fac1*y - fac2*vy ez = fac1*z - fac2*vz e2 = ex**2 + ey**2 + ez**2 xecen = sqrt(e2) delta = abs(1.0d0 - xecen) c..ascending node = zhat cross h components wx = -hy wy = hx wz = 0.0d0 w = sqrt(wx**2 + wy**2 + wz**2) c..perhelion distance qper = h*h/(gm * (1.0d0 + xecen)) c..inclination inc = acosd(hz/h) c..longitude of the ascending node, taking care of i=0 orbits if (w .le. 1.0d-12 .and. abs(wx) .le. 1.0d-12) then lan = 0.0d0 else lan = acosd(wx/w) if (wy .lt. 0.0) lan = 360.0d0 - lan end if c..argument of perihelion, taking care of e=0 or i=0 orbits fac1 = wx*ex + wy*ey + wz*ez if (abs(w*xecen) .lt. 1.0d-12) then aop = 0.0d0 else aop = acosd(fac1/(w*xecen)) if (ez .lt. 0.0) aop = 360.0d0 - aop end if c..for ellipses if (delta .gt. elimit .and. xecen .lt. 1.0d0) then c write(6,*) 'elliptical elemets' c..semimajor axis semi_maj = 1.0d0/(2.0d0/r - v2/gm) c..ecen*cos(eanom) and ecen*sin(eanom) c..hence the eccentric anomaly and mean anomaly c1 = 1.0d0 - r/semi_maj s1 = rdotv/(sqrt(semi_maj * gm)) eanom = atan2d(s1, c1) xmanom = eanom - s1*rad2a if (abs(xmanom) .lt. 1.0d-6) xmanom = 0.0d0 fac1 = xmanom if (xmanom .lt. 0.0) xmanom = xmanom + 360.0d0 c..period and perihelion date fac2 = sqrt((semi_maj*semi_maj*semi_maj)/gm) period = twopi * fac2 tjdp = tjd - fac1*a2rad * fac2 c..for hyperbolas else if (delta .gt. elimit .and. xecen .gt. 1.0d0) then c write(6,*) 'hyperbolic elements' c..semimajor axis is negative semi_maj = 1.0d0/(2.0d0/r - v2/gm) c..period is undefined for hyperbolas period = 1.0e30 c..ecen*cosh(eanom) and ecen*sinh(eanom) c..the eccentric anomaly is an inverse hyperbolic tangent, and mean anomaly c1 = 1.0d0 - r/semi_maj s1 = rdotv/(sqrt(-semi_maj * gm)) a = s1/c1 eanom = rad2a * 0.5d0*log( abs( (a + 1.0d0)/(a - 1.0d0) ) ) xmanom = s1*rad2a - eanom if (abs(xmanom) .lt. 1.0d-6) xmanom = 0.0d0 fac1 = xmanom if (xmanom .lt. 0.0) xmanom = xmanom + 360.0d0 c..period and perihelion date fac2 = sqrt((-semi_maj*semi_maj*semi_maj)/gm) tjdp = tjd - fac1*a2rad * fac2 c..for parabolas else if (delta .le. elimit) then c write(6,*) 'parabolic elements' c..semimajor axis and period are not defined for parbolas semi_maj = 1.0d30 period = 1.0e30 c..barkers equation to get the mean anomaly a = rdotv/sqrt(2.0d0 * gm * qper) xmanom = rad2a * a * (1.0d0 + a*a/3.0d0) if (abs(xmanom) .lt. 1.0d-6) xmanom = 0.0d0 fac1 = xmanom if (xmanom .lt. 0.0) xmanom = xmanom + 360.0d0 c..perihelion date fac2 = sqrt(2.0d0 * (qper*qper*qper)/gm) tjdp = tjd - fac1*a2rad * fac2 end if return end subroutine state_to_orb(gm,pos,vel,porb,vorb) implicit none save c..converts a state vector to an orbital plane state vector c..input: c..gm = gravitational constant = g*(m1 + m2) = kgrav**2 (1 + m2/m1) c.. where kgrav= kgauss = 0.01720209895d0 for orbits about the sun c.. kgrav= kearth = 0.07436680d0 for orbits about the earth c..pos = x y z coordinates c..vel = vx vy vz velocities c.. where the reference plane should be ecliptic coordinates c.. for bodies orbiting the sun, earth equatorial coordinates c.. for bodies orbiting the earth c..output: c..porb = xorb yorb zorb coordinates in the orbital plane c..vorb = vxorb vyorb vzorb velocities in the orbital plane c.. if all was done correctly zorb and vzorb should be zero, c.. to machine limits, since the orbit takes place in the c.. xorb-yorb plane. we include them as a useful check. c..declare the pass double precision gm,pos(3),vel(3),porb(3),vorb(3) c..for silly compilers that don't define the trig functions in degrees external cosd,sind,tand,acosd,asind,atan2d double precision cosd,sind,tand,acosd,asind,atan2d c..local variables integer i,indx(3) double precision x,y,z,vx,vy,vz, 1 r,v2,rdotv,hx,hy,hz,h, 2 fac1,fac2,ex,ey,ez,e2,xecen, 3 wx,wy,wz,w, 4 inc,lan,aop, 5 c1,s1,c2,s2,c3,s3,p(3),q(3),r1(3), 6 a(3,3),d c..unload for easier coding x = pos(1) y = pos(2) z = pos(3) vx = vel(1) vy = vel(2) vz = vel(3) c..distance, velocity squared, and rdotv r = sqrt(x**2 + y**2 + z**2) v2 = vx**2 + vy**2 + vz**2 rdotv = x*vx + y*vy + z*vz c..angular momentum = r cross v components hx = y*vz - z*vy hy = z*vx - x*vz hz = x*vy - y*vx h = sqrt(hx**2 + hy**2 + hz**2) c..eccentricity components, the eccentricity and its difference from unity fac1 = v2/gm - 1.0d0/r fac2 = rdotv/gm ex = fac1*x - fac2*vx ey = fac1*y - fac2*vy ez = fac1*z - fac2*vz e2 = ex**2 + ey**2 + ez**2 xecen = sqrt(e2) c..ascending node = zhat cross h components wx = -hy wy = hx wz = 0.0d0 w = sqrt(wx**2 + wy**2 + wz**2) c..inclination of the orbital plane inc = acosd(hz/h) c..longitude of the ascending node, taking care of i=0 orbits if (w .le. 1.0d-12 .and. abs(wx) .le. 1.0d-12) then lan = 0.0d0 else lan = acosd(wx/w) if (wy .lt. 0.0) lan = 360.0d0 - lan end if c..argument of perihelion, taking care of e=0 or i=0 orbits fac1 = wx*ex + wy*ey + wz*ez if (abs(w*xecen) .lt. 1.0d-12) then aop = 0.0d0 else aop = acosd(fac1/(w*xecen)) if (ez .lt. 0.0) aop = 360.0d0 - aop end if c..now get the gaussian vectors p, q, and r c1 = cosd(aop) s1 = sind(aop) c2 = cosd(inc) s2 = sind(inc) c3 = cosd(lan) s3 = sind(lan) p(1) = c1*c3 - s1*c2*s3 p(2) = c1*s3 + s1*c2*c3 p(3) = s1*s2 q(1) = -s1*c3 - c1*c2*s3 q(2) = -s1*s3 + c1*c2*c3 q(3) = c1*s2 r1(1) = s2*s3 r1(2) = -s2*c3 r1(3) = c2 c..form the conversion matrix from the gaussian vectors a(1,1) = p(1) a(2,1) = p(2) a(3,1) = p(3) a(1,2) = q(1) a(2,2) = q(2) a(3,2) = q(3) a(1,3) = r1(1) a(2,3) = r1(2) a(3,3) = r1(3) c..decompose the matrix, the right hand side is the given state vector, c..and then do the back substitution call ludcmp(a,3,3,indx,d) do i=1,3 porb(i) = pos(i) vorb(i) = vel(i) enddo call lubksb(a,3,3,indx,porb) call lubksb(a,3,3,indx,vorb) return end subroutine ludcmp(a,n,np,indx,d) implicit none save c..given the matrix a(n,n), with physical dimsnsions a(np,ap) this routine c..replaces a by the lu decompostion of a row-wise permutation of itself. c..input are a,n,np. output is a, indx which records the row c..permutations effected by the partial pivoting, and d which is 1 if c..the number of interchanges is even, -1 if odd. c..use routine lubksb to solve a system of linear equations. c..input: c..a(np,np) = matrix with physical dimensions np by np c..n = logical dimension of the matrix n .le. np c..np = physical dimension of the matrix c..output: c..indx(np) = vector of permutations c..d = even or odd number of permutations c..declare integer n,np,indx(np),nmax,i,j,k,imax parameter (nmax=50) double precision a(np,np),d,tiny,vv(nmax),aamax,sum,dum parameter (tiny=1.0d-20) if (np .gt. nmax) stop 'np > nmax in ludcmp' c..vv stores the implicit scaling of each row c..loop over the rows to get the scaling information d = 1.0d0 do i=1,n aamax = 0.0d0 do j=1,n if (abs(a(i,j)) .gt. aamax) aamax = abs(a(i,j)) enddo if (aamax .eq. 0.0) stop 'singular matrix in ludcmp' vv(i) = 1.0d0/aamax enddo c..for each column apply crouts method; see equation 2.3.12 do j=1,n do i=1,j-1 sum = a(i,j) do k=1,i-1 sum = sum - a(i,k)*a(k,j) enddo a(i,j) = sum enddo c..find the largest pivot element aamax = 0.0d0 do i=j,n sum = a(i,j) do k=1,j-1 sum = sum - a(i,k)*a(k,j) enddo a(i,j) = sum dum = vv(i)*abs(sum) if (dum .ge. aamax) then imax = i aamax = dum end if enddo c..if we need to interchange rows if (j .ne. imax) then do k=1,n dum = a(imax,k) a(imax,k) = a(j,k) a(j,k) = dum enddo d = -d vv(imax) = vv(j) end if c..divide by the pivot element indx(j) = imax if (a(j,j) .eq. 0.0) a(j,j) = tiny if (j .ne. n) then dum = 1.0d0/a(j,j) do i=j+1,n a(i,j) = a(i,j)*dum enddo end if c..and go back for another column of crouts method enddo return end subroutine lubksb(a,n,np,indx,b) implicit none save c..solves a set of n linear equations ax=b. a is input in its lu decomposition c..form, determined by the routine above ludcmp. indx is input as the c..permutation vector also returned by ludcmp. b is input as the right hand c..side vector and returns with the solution vector x. c..a,n ans np are not modified by this routine and thus can be left in place c..for successive calls (i.e matrix inversion) c..input: c..a(np,p) = lu decomposed matrix c..n = logical dimension of matrix c..np = physical dimension of matrix c..indx(np) = vector of row permutations c..b = right hand side c..output: c..b = solution c..declare integer n,np,indx(np),i,ii,j,ll double precision a(np,np),b(np),sum c..when ii is > 0, ii becomes the index of the first nonzero element of b c..this is forward substitution of equation 2.3.6, and unscamble in place ii = 0 do i=1,n ll = indx(i) sum = b(ll) b(ll) = b(i) if (ii .ne. 0) then do j=ii,i-1 sum = sum - a(i,j) * b(j) enddo c..nonzero element was found, so dos the sums in the loop above else if (sum .ne. 0.0) then ii = i end if b(i) = sum enddo c..back substitution equation 2.3.7 do i = n,1,-1 sum = b(i) if (i .lt. n) then do j=i+1,n sum = sum - a(i,j) * b(j) enddo end if b(i) = sum/a(i,i) enddo return end subroutine rotax(rin,theta,rout) implicit none save c..the rotation matrix times a vector, rotate about x axis c..the sign is chosen such that a positive angle theta c..corresponds to a positive rotation (i.e. an antilcockwise rotation) c..of the reference frame around the rotational axis. c..on input, the array rout can be the same as array rin double precision rin(3),theta,rout(3),ct,st,a,b,c c..for silly compilers that don't define the trig functions in degrees external cosd,sind,tand,acosd,asind,atan2d double precision cosd,sind,tand,acosd,asind,atan2d ct = cosd(theta) st = sind(theta) a = rin(1) b = ct*rin(2) + st*rin(3) c = -st*rin(2) + ct*rin(3) rout(1) = a rout(2) = b rout(3) = c return end subroutine rotay(rin,theta,rout) implicit none save c..the rotation matrix times a vector, rotate about y axis c..the sign is chosen such that a positive angle theta c..corresponds to a positive rotation (i.e. an antilcockwise rotation) c..of the reference frame around the rotational axis. c..on input, the array rout can be the same as array rin double precision rin(3),theta,rout(3),ct,st,a,b,c c..for silly compilers that don't define the trig functions in degrees external cosd,sind,tand,acosd,asind,atan2d double precision cosd,sind,tand,acosd,asind,atan2d ct = cosd(theta) st = sind(theta) a = ct*rin(1) - st*rin(3) b = rin(2) c = st*rin(1) + ct*rin(3) rout(1) = a rout(2) = b rout(3) = c return end subroutine rotaz(rin,theta,rout) implicit none save c..the rotation matrix times a vector, rotate about z axis c..the sign is chosen such that a positive angle theta c..corresponds to a positive rotation (i.e. an antilcockwise rotation) c..of the reference frame around the rotational axis. c..on input, the array rout can be the same as array rin double precision rin(3),theta,rout(3),ct,st,a,b,c c..for silly compilers that don't define the trig functions in degrees external cosd,sind,tand,acosd,asind,atan2d double precision cosd,sind,tand,acosd,asind,atan2d ct = cosd(theta) st = sind(theta) a = ct*rin(1) + st*rin(2) b = -st*rin(1) + ct*rin(2) c = rin(3) rout(1) = a rout(2) = b rout(3) = c return end subroutine requ_to_recl(xequ,yequ,zequ,epsil,xecl,yecl,zecl) implicit none save c..converts rectangular equatorial coordinates to c..rectangular ecliptic coordinates c..input: c..xequ,yequ,zequ = equatorial coordinates c..epsil = true obliquity of the ecliptic c..output: c..xecl,yecl,zecl = ecliptic coordinates c..declare the pass double precision xecl,yecl,zecl,xequ,yequ,zequ,epsil c..local variables double precision cepsil,sepsil c..for silly compilers that don't define the trig functions in degrees external cosd,sind,tand,acosd,asind,atan2d double precision cosd,sind,tand,acosd,asind,atan2d c..rotate cepsil = cosd(epsil) sepsil = sind(epsil) xecl = xequ yecl = yequ*cepsil + zequ*sepsil zecl = -yequ*sepsil + zequ*cepsil return end subroutine recl_to_requ(xecl,yecl,zecl,epsil,xequ,yequ,zequ) implicit none save c..converts rectangular ecliptic coordinates to c..rectangular equatorial coordinates c..input: c..xecl,yecl,zecl = ecliptic coordinates c..epsil = true obliquity of the ecliptic c..output: c..xequ,yequ,zequ = equatorial coordinates c..declare the pass double precision xecl,yecl,zecl,xequ,yequ,zequ,epsil c..local variables double precision cepsil,sepsil c..for silly compilers that don't define the trig functions in degrees external cosd,sind,tand,acosd,asind,atan2d double precision cosd,sind,tand,acosd,asind,atan2d c..rotate cepsil = cosd(epsil) sepsil = sind(epsil) xequ = xecl yequ = yecl*cepsil - zecl*sepsil zequ = yecl*sepsil + zecl*cepsil return end subroutine requ_to_rhor(xequ,yequ,zequ,tjd,glon,glat, 1 xhor,yhor,zhor) implicit none save c..converts rectangular equatorial coordinates to c..rectangular horizon coordinates c..note that since the sidereal time is used in the conversion, c..that the passed julian date should be the ut julian date, c..not the dynamical julian date. c.. c..the orientation is: c..north is 0 degrees azimuth c..east is 90 degrees azimuth c..south is 180 degrees azimuth c..west is 270 degrees azimuth c..input: c..xequ,yequ,zequ = equatorial coordinates c..tjd = ut julian day number, not the ephemeris julian day number c..glon = geographic longitude (+ west of greenwhich, - east of greenwhich) c..glat = geographic latitude (+ north of equator, - south) c..output: c..xhor,yhor,zhor = alt-az coordinates c..declare the pass double precision xequ,yequ,zequ,tjd,glon,glat,xhor,yhor,zhor c..local variables double precision tjdh,tjdl,gst,lst,cglat,sglat,clst,slst,x1,y1,z1 c..for silly compilers that don't define the trig functions in degrees external cosd,sind,tand,acosd,asind,atan2d double precision cosd,sind,tand,acosd,asind,atan2d c..get the local sidereal time in degrees tjdh = int(tjd) tjdl = tjd - tjdh call sidtim(tjdh,tjdl,1,gst) lst = 15.0d0*gst - glon if (lst .lt. 360.0) lst = lst + 360.0d0 c..rotate about the z axis through local sidereal time clst = cosd(lst) slst = sind(lst) x1 = xequ*clst + yequ*slst y1 = xequ*slst - yequ*clst z1 = zequ c..rotate about the y axis through the geographic latitude cglat = cosd(glat) sglat = sind(glat) xhor = -x1*sglat + z1*cglat yhor = -y1 zhor = x1*cglat + z1*sglat return end subroutine rhor_to_requ(xhor,yhor,zhor,tjd,glon,glat, 1 xequ,yequ,zequ) implicit none save c..converts rectangular horizon coordinates to c..rectangular horizon coordinates. c..note that since the sidereal time is used in the conversion, c..that the passed julian date should be the ut julian date, c..not the dynamical julian date. c..the orientation is: c..north is 0 degrees azimuth c..east is 90 degrees azimuth c..south is 180 degrees azimuth c..west is 270 degrees azimuth c..input: c..xhor,yhor,zhor = alt-az coordinates c..tjd = julian day number c..glon = geographic longitude (+ west of greenwhich, - east of greenwhich) c..glat = geographic latitide (+ north of equator, - south) c..output: c..xequ,yequ,zequ = equatorial coordinates c..declare the pass double precision xequ,yequ,zequ,tjd,glon,glat,xhor,yhor,zhor c..local variables double precision tjdh,tjdl,gst,lst,cglat,sglat,clst,slst,x1,y1,z1 c..for silly compilers that don't define the trig functions in degrees external cosd,sind,tand,acosd,asind,atan2d double precision cosd,sind,tand,acosd,asind,atan2d c..get the local sidereal time in degrees tjdh = int(tjd) tjdl = tjd - tjdh call sidtim(tjdh,tjdl,1,gst) lst = 15.0d0*gst - glon if (lst .lt. 360.0) lst = lst + 360.0d0 c..rotate about the y axis through the geographic latitude cglat = cosd(glat) sglat = sind(glat) x1 = -xhor*sglat + zhor*cglat y1 = -yhor z1 = xhor*cglat + zhor*sglat c..rotate about the z axis through the hour angle clst = cosd(lst) slst = sind(lst) xequ = x1*clst + y1*slst yequ = x1*slst - y1*clst z1 = zequ return end subroutine dh2h(dec,h,m,s) implicit none save c..converts decimal hours/degees into hour/degree min sec format. c..input: c..dec = decimal degrees or decimal hours c..output c..h = whole degrees or hours c..m = whole minutes c..s = seconds double precision dec,h,m,s,x x = dec h = int(x) x = (x - h) * 60.0d0 m = int(x) s = (x - m) * 60.0d0 return end double precision function geodetic(phi) implicit none save c..this function is to be called by a root finder to find the c..geodedtic latitude phi. c..expressions from the explanatory supplement, page 126. c..our only approximation here is to take ht as the height above c..mean sea level, as opposed to its proper definition as the c..height above the reference ellipsoid. c..input: c..phi = geodedtic latitude in decimal degrees c..glat = geocentric latitude in decimal degrees c..ht = height above mean sea level in meters c..output: c..geodetic = geodetic equation as f(x) = 0 c..declare the pass double precision phi c..communicate with higher level routines, c..must be filled by the calling routine double precision glat_com,ht_com common /gedc1/ glat_com,ht_com c..math constants double precision pi,a2rad,rad2a parameter (pi = 3.1415926535897932384d0, 1 a2rad = pi/180.0d0, rad2a = 180.0d0/pi) c..for silly compilers that don't define the trig functions in degrees external cosd,sind,tand,acosd,asind,atan2d double precision cosd,sind,tand,acosd,asind,atan2d c..local variables c..f is the grs80 ellipsoid value (supplement page 223) flattening value c..and rearth is the 1986 equatorial distance from center to surface double precision ct,st,c,s,x,y,f,f2,rearth parameter (f = 1.0d0/298.257222101d0, 1 f2 = (1.0d0 - f)*(1.0d0 - f), 2 rearth = 6378.137d0) c..go ct = cosd(phi) st = sind(phi) c = 1.0d0/sqrt(ct**2 + f2*st**2) s = f2 * c x = (rearth*s + ht_com*1.0d-3)*st y = (rearth*c + ht_com*1.0d-3)*ct geodetic = glat_com - atan2d(x,y) c..here is a series expansion c geodetic = phi*a2rad - glat_com*a2rad c 1 - (f + 0.5d0*f*f)*sind(2.0d0*phi) c 2 + (0.5d0*f*f + 0.5d0*f**3)*sind(4.0d0*phi) c 3 - (1.0d0/3.0d0 * f**3)*sind(6.0d0*phi) return end double precision function zbrent(func,x1,x2,tol,niter) implicit none save c..using brent's method this routine finds the root of a function func c..between the limits x1 and x2. the root is when accuracy is less than tol. c..input: c..func = external function name that contains f(x) = 0 c..x1 = lower bracket c..x2 = upper bracket c..tol = fractional accuracy to find the root to c..output: c..niter = number of iterations taken c..zbrent = value of the root x*, which makes f(x*) = 0 c..declare external func integer niter,itmax,iter parameter (itmax = 100) double precision func,x1,x2,tol,a,b,c,d,e,fa, 1 fb,fc,xm,tol1,p,q,r,s,eps parameter (eps=3.0d-15) c..note: eps is an estimate of the machine floating point precision c..initialize niter = 0 a = x1 b = x2 fa = func(a) fb = func(b) if ( (fa .gt. 0.0 .and. fb .gt. 0.0) .or. 1 (fa .lt. 0.0 .and. fb .lt. 0.0) ) then write(6,100) x1,fa,x2,fb 100 format(1x,' x1=',1pe11.3,' f(x1)=',1pe11.3,/, 1 1x,' x2=',1pe11.3,' f(x2)=',1pe11.3) stop 'root not bracketed in routine zbrent' end if c = b fc = fb c..rename a,b,c and adjusting bound interval d do iter =1,itmax niter = niter + 1 if ( (fb .gt. 0.0 .and. fc .gt. 0.0) .or. 1 (fb .lt. 0.0 .and. fc .lt. 0.0) ) then c = a fc = fa d = b-a e = d end if if (abs(fc) .lt. abs(fb)) then a = b b = c c = a fa = fb fb = fc fc = fa end if tol1 = 2.0d0 * eps * abs(b) + 0.5d0 * tol xm = 0.5d0 * (c-b) c..convergence check if (abs(xm) .le. tol1 .or. fb .eq. 0.0) then zbrent = b return end if c..attempt quadratic interpolation if (abs(e) .ge. tol1 .and. abs(fa) .gt. abs(fb)) then s = fb/fa if (a .eq. c) then p = 2.0d0 * xm * s q = 1.0d0 - s else q = fa/fc r = fb/fc p = s * (2.0d0 * xm * q *(q-r) - (b-a)*(r - 1.0d0)) q = (q - 1.0d0) * (r - 1.0d0) * (s - 1.0d0) end if c..check if in bounds if (p .gt. 0.0) q = -q p = abs(p) c..accept interpolation if (2.0d0*p .lt. min(3.0d0*xm*q - abs(tol1*q),abs(e*q))) then e = d d = p/q c..or bisect else d = xm e = d end if c..bounds decreasing to slowly use bisection else d = xm e = d end if c..move best guess to a a = b fa = fb if (abs(d) .gt. tol1) then b = b + d else b = b + sign(tol1,xm) end if fb = func(b) enddo stop 'too many iterations in routine zbrent' end subroutine zbrac(func,x1,x2,succes) implicit none save c..given a function func and an initial guessed range x1 to x2, the routine c..expands the range geometrically until a root is bracketed by the returned c..values x1 and x2 (in which case succes returns as .true.) or until the c..range becomes unacceptably large (in which succes returns as .false.). c..success guaranteed for a function which has opposite sign for sufficiently c..large and small arguments. c..input: c..func = external function name which contains f(x) = 0 c..x1 = guess at the lower bracket c..x2 = guess at the upper bracket c..output: c..x1 = bracketing lower bound c..x2 = bracketing upper bound c..succes = logical indicating if the bracketing was successful c..declare external func logical succes integer ntry,j parameter (ntry=50) double precision func,x1,x2,factor,f1,f2 parameter (factor = 1.6d0) if (x1 .eq. x2) stop ' x1 = x2 in routine zbrac' f1 = func(x1) f2 = func(x2) succes = .true. do j=1,ntry if (f1*f2 .lt. 0.0) return if (abs(f1) .lt. abs(f2)) then x1 = x1 + factor * (x1-x2) f1 = func(x1) else x2 = x2 + factor * (x2-x1) f2 = func(x2) end if enddo succes = .false. return end double precision function mnbrent(ax,bx,cx,func,tol,xmin,niter) implicit none save c..finds the minima of a one-dimensional function c..input: c..ax bx cx = x-coordinates which bracket the minimum c.. such that fb < fa and fb < fc c..func = function which contains the expression we seek the minima of c..tol = fractional accuaracy with which to find the minima c..output c..xmin = x-coordinate of the minima c..mnbrent = y-coordinate of the minima c..niter = number of function evaluations required to find the minima c..declare the pass external func integer niter double precision ax,bx,cx,func,tol c..local variables c..cgold is the golden ratio c..zeps protects against a minimum that happens to be exactly zero. integer itmax,iter parameter (itmax=100) double precision xmin,cgold,zeps, 1 a,b,v,w,e,d,fx,fu,fv,fw,xm,tol1,tol2,r,p,q, 2 etemp,u,x parameter (cgold = 0.381966011250105d0, zeps=1.0d-14) c..a and b in ascending order; e is distance moved on the step before last a = min(ax,cx) b = max(ax,cx) v = bx w = v x = v e = 0.0d0 fx = func(x) fv = fx fw = fx niter = 1 c..main loop do iter=1,itmax xm = 0.5d0 * (a+b) tol1 = tol*abs(x)+zeps tol2 = 2.0d0 * tol1 c..bail if we are finished if (abs(x-xm) .le. (tol2 - 0.5d0*(b-a))) then xmin = x mnbrent = fx return end if c..construct a trial parabolic fit if (abs(e) .gt. tol1) then r = (x-w)*(fx-fv) q = (x-v)*(fx-fw) p = (x-v)*q - (x-w)*r q = 2.0d0 * (q-r) if (q .gt. 0.0) p=-p q = abs(q) etemp = e e = d c..determine the acceptability of the parabolic fit if (abs(p) .ge. abs(0.5d0*q*etemp) .or. 1 p .le. q * (a-x) .or. 2 p .ge. q * (b-x)) go to 1 c..parabolic fit ok, take the parabolic step and skip the golden mean section d = p/q u = x + d if (u-a .lt. tol2 .or. b-u .lt. tol2) d=sign(tol1,xm-x) go to 2 end if c..golden mean section 1 if (x .ge. xm) then e = a - x else e = b - x end if d = cgold*e c..arrives here with d calculated from either the golden mean or parabolic fit c..note the one function evaluation per iteration 2 if (abs(d) .ge. tol1) then u = x + d else u = x + sign(tol1,d) end if fu = func(u) niter = niter + 1 c..now decide what to do with the function evaluation and the housekeeping if (fu. le. fx) then if (u .ge. x) then a = x else b = x end if v = w fv = fw w = x fw = fx x = u fx = fu else if (u .lt. x) then a = u else b = u end if if (fu .le. fw .or. w .eq. x) then v = w fv = fw w = u fw = fu else if (fu .le. fv .or. v .eq. x .or. v .eq. w) then v = u fv = fu end if end if enddo stop 'exceed maximum iterations in function mnbrent' end subroutine mnbrak(ax,bx,cx,fa,fb,fc,func) implicit none save c..this routine brackets the minimum of a one-dimensional function c..input: c..ax = one of two distinct initial x-coordinates c..bx = one of two distinct initial x-coordinates c..func = function which contains the expression we seek the minima of c..output: c..ax bx cx = x-coordinates which bracket the minimum c..fa fb fc = y-coordinates of ax bx and cx such that fb < fa and fb < fc c..declare the pass external func double precision ax,bx,cx,fa,fb,fc,func c..local variables double precision gold,glimit,tiny, 1 dum,r,q,u,ulim,fu parameter (gold=1.618033398874989d0, 1 glimit=100.0d0, tiny=1.0d-20) c..gold is the inverse of the golden ratio, and is the ratio with c..which succesive intervals are magnified. glimit is the maximum c..magnification allowed for a parabolic fit. c..get the initial function values fa = func(ax) fb = func(bx) c..switch a and b if needed so that we can go downhill from a to b if (fb .gt. fa) then dum = ax ax = bx bx = dum dum = fb fb = fa fa = dum end if c..first guess for c cx = bx + gold*(bx-ax) fc = func(cx) c..a do while construction, keep returning here until bracketed 01 if (fb .ge. fc) then c..compute u by parabolic extrapolation from a, b and c r = (bx-ax)*(fb-fc) q = (bx-cx)*(fb-fa) u = bx-((bx-cx)*q-(bx-ax)*r) / 1 (2.0d0 * sign(max(abs(q-r),tiny),q-r)) ulim = bx + glimit*(cx-bx) c..now test the various possibilities c..for u between b and c if ((bx-u)*(u-cx) .gt. 0.0) then fu = func(u) c..got a minimum between b and c; relabel a, b and c and bail if (fu .lt. fc) then ax = bx fa = fb bx = u fb = fu return c..got a minimum between a and u; relabel and bail else if (fu .gt. fb) then cx = u fc = fu return end if c..got no minimum; parabolic fit was of no use; try the golden mean u = cx + gold*(cx-bx) fu = func(u) c..so u is not between b and c c..if u is between c and its allowed limit, relabel else if ((cx-u)*(u-ulim) .gt. 0.0) then fu = func(u) if (fu .lt. fc) then bx = cx cx = u u = cx + gold*(cx-bx) fb = fc fc = fu fu = func(u) end if c..if u is past the maximum allowed, limit it else if ((u-ulim)*(ulim-cx) .ge. 0.0) then u = ulim fu = func(u) c..complete rejection of parabolic fit; back to golden mean else u = cx + gold*(cx-bx) fu = func(u) end if c..eliminate the oldest point and continue ax = bx bx = cx cx = u fa = fb fb = fc fc = fu goto 1 end if return end integer function getnam(string,word,ipos) implicit none save c..this routine is a basic string parser, only spaces and equal signs are c..used as delimiters. c.. c..input: c..string = character string to chop c..ipos = where in the string to start chopping c.. c..output: c..word = the parsed word of the string c..ipos = new starting position of where to start chopping c..getnam = length of the word c..if no word is found, meaning that there are no more words in the string, c..then getnam is returned as zero and the word is filled with spaces. c..declare character*(*) string,word integer kpos,ipos,k,lend,nbegin c..get the length of the input line, blank word and save ipos lend = len(string) word = ' ' kpos = ipos c..find begining of the word; if nothing found return getnam as zero do k = kpos,lend nbegin = k if (string(k:k).ne.' ' .and. string(k:k).ne.'=') goto 25 enddo getnam = 0 return c..find end of the word 25 continue do k = nbegin,lend ipos = k if( string(k:k).eq.' ' .or. string(k:k).eq.'=') goto 35 enddo c..build the word, set getnam and return 35 continue word = string(nbegin:ipos-1) getnam = ipos - nbegin return end double precision function value(string) implicit none save c..this routine takes a character string and converts it to a real number. c..on error during the conversion, a fortran stop is issued in this version c..input: c..string = character string c..output: c..value = real value of characters in string c..declare logical pflag character*(*) string character*1 plus,minus,decmal,blank,se,sd,se1,sd1 integer noblnk,long,ipoint,power,psign,iten,j,z,i double precision x,sign,factor,rten,temp parameter (plus = '+' , minus = '-' , decmal = '.' , 1 blank = ' ' , se = 'e' , sd = 'd' , 2 se1 = 'E' , sd1 = 'D' , rten = 10.0, 3 iten = 10 ) c..initialize x = 0.0d0 sign = 1.0d0 factor = rten pflag = .false. noblnk = 0 power = 0 psign = 1 long = len(string) c..remove any leading blanks and get the sign of the number do z = 1,7 noblnk = noblnk + 1 if ( string(noblnk:noblnk) .eq. blank) then if (noblnk .gt. 6 ) goto 30 else if (string(noblnk:noblnk) .eq. plus) then noblnk = noblnk + 1 else if (string(noblnk:noblnk) .eq. minus) then noblnk = noblnk + 1 sign = -1.0d0 end if goto 10 end if enddo c..main number conversion loop 10 continue do i = noblnk,long ipoint = i + 1 c..if a blank character then we are done if ( string(i:i) .eq. blank ) then x = x * sign value = x return c..if an exponent character, process the whole exponent, and return else if (string(i:i).eq.se .or. string(i:i).eq.sd .or. 1 string(i:i).eq.se1 .or. string(i:i).eq.sd1 ) then if (x .eq. 0.0 .and. ipoint.eq.2) x = 1.0d0 if (sign .eq. -1.0 .and. ipoint.eq.3) x = 1.0d0 if (string(ipoint:ipoint) .eq. plus) ipoint = ipoint + 1 if (string(ipoint:ipoint) .eq. minus) then ipoint = ipoint + 1 psign = -1 end if do z = ipoint,long if (string(z:z) .eq. blank) then x = sign * x * rten**(power*psign) value = x return else j = ichar(string(z:z)) - 48 if ( (j.lt.0) .or. (j.gt.9) ) goto 30 power= (power * iten) + j end if enddo c..if an ascii number character, process ie else if (string(i:i) .ne. decmal) then j = ichar(string(i:i)) - 48 if ( (j.lt.0) .or. (j.gt.9) ) goto 30 if (.not.(pflag) ) then x = (x*rten) + j else temp = j x = x + (temp/factor) factor = factor * rten goto 20 end if c..must be a decimal point if none of the above c..check that there are not two decimal points! else if (pflag) goto 30 pflag = .true. end if 20 continue end do c..if we got through the do loop ok, then we must be done x = x * sign value = x return c..error processing the number 30 write(6,40) long,string(1:long) 40 format(' error converting the ',i4,' characters ',/, 1 ' >',a,'< ',/, 2 ' into a real number in function value') stop ' error in routine value' end subroutine twobody(mu,t1,x1,y1,z1,vx1,vy1,vz1, 1 t2,x2,y2,z2,vx2,vy2,vz2) implicit none save c..solution of the 2 body initial value problem for newtonian gravity. c..using the universal kepler equation and the f and g functions, c..so its valid for circular, elliptical, hyperbolic and parabolic orbits c..with no special checks for which type of orbit. c..derived from danby, page 179 c..input: c..mu = g * (m1 + m2) in units of your choice c..t1 = time point for initial positions and velocities c..x1,y1,z1 = initial positions c..vx1,vy1,vz1 = initial velocities c..t = time point at which a solution is desired c..output: c..x,y,z = position vector at time t c..vx,vy,vz = velocity vector at time t c..declare the pass integer niter double precision mu,t1,x1,y1,z1,vx1,vy1,vz1, 1 t2,x2,y2,z2,vx2,vy2,vz2 c..for the solution of the universal kepler equation external keprot logical succes double precision keprot,zbrent,slo,shi,s,tol parameter (tol = 1.0d-14) c..communication between routine twobody and function keprot double precision muu,r0,rdotv,alpha,dt,c1,c2,c3,keprim common /bod2com/ muu,r0,rdotv,alpha,dt,c1,c2,c3,keprim c..local variables double precision v0s,f,g,fdot,gdot c..set the common block information muu = mu r0 = sqrt(x1*x1 + y1*y1 + z1*z1) v0s = vx1*vx1 + vy1*vy1 + vz1*vz1 rdotv = x1*vx1 + y1*vy1 + z1*vz1 alpha = 2.0d0 * mu /r0 - v0s dt = t2 - t1 c..guess the root, bracket it and solve kepler's equation slo = dt/r0 shi = 1.1d0 * slo slo = 0.9d0 * slo call zbrac(keprot,slo,shi,succes) if (.not.succes) then write(6,*) 'could not bracket root in routine twobody' write(6,*) 'slo =',slo write(6,*) 'shi =',shi stop 'error: routine zbrac failed in twobody' end if c..get the root, s is the arc length s = zbrent(keprot,slo,shi,tol,niter) c..form the output via the f and g functions f = 1.0d0 - (mu/r0)*c2 g = dt - mu*c3 fdot = -(mu/(keprim*r0))*c1 gdot = 1.0d0 - (mu/keprim)*c2 x2 = x1 * f + vx1 * g y2 = y1 * f + vy1 * g z2 = z1 * f + vz1 * g vx2 = x1 * fdot + vx1 * gdot vy2 = y1 * fdot + vy1 * gdot vz2 = z1 * fdot + vz1 * gdot return end double precision function keprot(s) implicit none save c..sets up the universal kepler equation for a root find. c..derived from danby, page 179 c..input: c..s = path length c..output: c..keprot = value of the universal kepler equation c..declare double precision s,x,c0 c..communication between routine twobody and function keprot double precision muu,r0,rdotv,alpha,dt,c1,c2,c3,keprim common /bod2com/ muu,r0,rdotv,alpha,dt,c1,c2,c3,keprim c..form the stumpff functions x = s * s * alpha call scaled_stumpff(x,c0,c1,c2,c3) c1 = c1 * s c2 = c2 * s * s c3 = c3 * s * s * s c..hence the universal keplers equation and its derivative keprot = r0*c1 + rdotv*c2 + muu*c3 - dt keprim = r0*c0 + rdotv*c1 + muu*c2 return end subroutine derv(t,y,dydt) implicit none save c..this routine evaluates right hand sides c..for the comet integrations c..input is the state vector y, with y(1), y(2), and y(3) c..being the heliocentric ecliptic rectangular position coordinates (x,y,z) c..and y(4), y(5), and y(6) being the heliocentric ecliptic rectangular c..speed coordinates (vx,vy,vz) c..output is dydt, which contains the rates of change. c..in other words this routine evaluates the accelerations on c..a comet taking into account the planetary perturbations c..declare the pass double precision t,y(*),dydt(*) c..local variables integer ifirst,imercury,ivenus,iearth,imars,ijupiter, 1 isaturn,iuranus,ineptune,ipluto,i parameter (imercury=1,ivenus=2,iearth=3,imars=4, 1 ijupiter=5,isaturn=6,iuranus=7,ineptune=8, 2 ipluto=9) double precision xc,yc,zc,vxc,vyc,vzc,ax,ay,az, 1 r2,r,mur3,dx,dy,dz, 2 xp,yp,zp,rrd(6),tjd2000, 3 oblm,oblt,eqeq,dpsi,deps parameter (tjd2000 = 2451545.0d0) c..commuication between routine auxpos and this routine for the origin integer icenter common /auxint1/ icenter c..physical constants double precision kgauss,gm_sun,mu(9) parameter (kgauss = 0.01720209895d0, 1 gm_sun = kgauss*kgauss) c..jpl 1996 reference values c..each is the sum of the parent body + any moons + any rings if (ifirst .eq. 0) then ifirst = 1 mu(imercury) = gm_sun/6023600.0d0 mu(ivenus) = gm_sun/408523.71d0 mu(iearth) = gm_sun/328900.56d0 mu(imars) = gm_sun/3098708.0d0 mu(ijupiter) = gm_sun/1047.3486d0 mu(isaturn) = gm_sun/3497.898d0 mu(iuranus) = gm_sun/22902.98d0 mu(ineptune) = gm_sun/19412.24d0 mu(ipluto) = gm_sun/1.35e8 call etilt (tjd2000,oblm,oblt,eqeq,dpsi,deps) end if c..load the input xc = y(1) yc = y(2) zc = y(3) vxc = y(4) vyc = y(5) vzc = y(6) c..solar attraction r2 = xc*xc + yc*yc + zc*zc r = sqrt(r2) mur3 = gm_sun/(r*r2) ax = -mur3*xc ay = -mur3*yc az = -mur3*zc c..for each planet do i=1,9 c..get heliocentric equatorial rectangular coordinates of planet from pleph c..and convert to heliocentric ecliptic coordinates call pleph (t, i, icenter, rrd ) call requ_to_recl(rrd(1),rrd(2),rrd(3),oblm,xp,yp,zp) c..direct acceleration term dx = xc - xp dy = yc - yp dz = zc - zp r2 = dx*dx + dy*dy + dz*dz r = sqrt(r2) mur3 = mu(i)/(r*r2) ax = ax - mur3*dx ay = ay - mur3*dy az = az - mur3*dz c..indirect acceleration term r2 = xp*xp + yp*yp + zp*zp r = sqrt(r2) mur3 = mu(i)/(r*r2) ax = ax - mur3*xp ay = ay - mur3*yp az = az - mur3*zp end do c..load the output dydt(1) = vxc dydt(2) = vyc dydt(3) = vzc dydt(4) = ax dydt(5) = ay dydt(6) = az return end subroutine jaco(t,y,dfdy,nlog,nphys) implicit none save c..just a dummy stub, in case one ever wants to try an implicit c..integration of the orbit equations. integer nlog,nphys double precision t,y(*),dfdy(nphys,nphys) return end subroutine cometint(start,stptry,stpmin,stopp,bc, 1 eps,dxsav,kmax, 2 xrk,yrk,xphys,yphys,xlogi,ylogi, 3 nok,nbad,kount,odescal,iprint, 4 derivs,jacob,steper) implicit none save c..drives the comet integrations c..declare external derivs,jacob,steper integer nok,nbad,nmax,stpmax,kmax,kount,xphys, 1 yphys,xlogi,ylogi,iprint,i,j,nstp parameter (nmax = 10, stpmax=10000) double precision bc(yphys),stptry,stpmin,eps,stopp,start, 1 yscal(nmax),y(nmax),dydx(nmax), 2 dxsav,xrk(xphys),yrk(yphys,xphys), 3 odescal,x,xsav,h,hdid,hnext,zero,one,tiny,ttiny parameter (zero=0.0d0, one=1.0d0, 1 tiny=1.0d-30, ttiny=1.0d-15) c..here are the format statements for printouts as we integrate 100 format(1x,i4,1pe10.2) 101 format(1x,1p12e10.2) c..initialize if (ylogi .gt. yphys) stop 'ylogi > yphys in routine cometint' if (yphys .gt. nmax) stop 'yphys > nmax in routine cometint' x = start h = sign(stptry,stopp-start) nok = 0 nbad = 0 kount = 0 c..store the first step do i=1,ylogi y(i) = bc(i) enddo xsav = x - 2.0d0 * dxsav c..take at most stpmax steps do nstp=1,stpmax call derivs(x,y,dydx) c..scaling vector used to monitor accuracy do i=1,ylogi c..constant fractional accuracy c yscal(i) = abs(y(i)) + tiny c..const. frac. cept near zero c yscal(i)=abs(y(i)) + abs(h*dydx(i)) + tiny c..step size dependent accuracy c yscal(i) = abs(odescal * h * dydx(i)) + tiny c..for stiffs yscal(i) = max(odescal,y(i)) c..strait scaling (decrease to get more time steps) c yscal(i) = odescal enddo c..store intermediate results if (kmax .gt. 0) then if ( (abs(dxsav) - abs(x-xsav)) .le. ttiny) then if ( kount .lt. (kmax-1) ) then kount = kount+1 xrk(kount) = x do i=1,ylogi yrk(i,kount) = y(i) enddo if (iprint .eq. 1) then write(6,100) kount,xrk(kount) c write(6,101) (yrk(j,kount), j=1,ylogi) end if xsav=x end if end if end if c..if the step can overshoot the stop point or the dxsav increment then cut it if ((x+h-stopp)*(x+h-start).gt.zero) h = stopp - x if (dxsav.ne.zero .and. h.gt.(xsav-x+dxsav)) h = xsav + dxsav-x c..do an integration step call steper(y,dydx,ylogi,x,h,eps,yscal,hdid,hnext,derivs,jacob) if (hdid.eq.h) then nok = nok+1 else nbad = nbad+1 end if c..this is the normal exit point, save the final step if (nstp .eq. stpmax .or. (x-stopp)*(stopp-start) .ge. zero) then do i=1,ylogi bc(i) = y(i) enddo if (kmax.ne.0) then kount = kount+1 xrk(kount) = x do i=1,ylogi yrk(i,kount) = y(i) enddo if (iprint .eq. 1) then write(6,100) kount,xrk(kount) c write(6,101) (yrk(j,kount), j=1,ylogi) end if end if return end if c..set the step size for the next iteration; stay above stpmin h=hnext if (abs(hnext).lt.stpmin) stop 'hnext < stpmin in cometint' c..back for another iteration or death enddo write(6,*) '> than stpmax steps required in cometint' return end subroutine bsstep(y,dydx,nv,x,htry,eps,yscal,hdid,hnext, 1 derivs,jacob) implicit none save c..bulirsch-stoer step with monitoring of local truncation error for c..accuracy. input are the dependent variable vector y(1:nv) and its c..derivative dydx(1:nv) at the independent variable x. also input are c..the step size to be attempted htry, the requested accuracu eps, and c..the scaling vector yscal(1:nv). on output y and x are replaced by their c..new values, hdid is the step that was actually used and hnext is the c..estimated next step to use. derivs is a user supplied routine that returns c..the right hand side of the ode's. plug into odeint. c.. c..note: c..kmaxx is the maximum row number used in the extrapolation c..imax is the next row number; safe1 and safe2 are safety factors c..redmax is the maxium factor a step is reduced; redmin the minimum c..1/sclmax is the maximum factor by which a stepsize can be increased. c.. c..declare external derivs,jacob logical first,reduct integer nv,nmax,kmaxx,imax parameter (nmax = 20, kmaxx=8, imax=kmaxx+1) integer i,iq,k,kk,km,kmax,kopt,nseq(imax) double precision y(nv),dydx(nv),x,htry,eps,yscal(nv),hdid,hnext, 1 eps1,epsold,errmax,fact,h,red,scale,work,wrkmin, 2 xest,xnew,a(imax),alf(kmaxx,kmaxx),err(kmaxx), 3 yerr(nmax),ysav(nmax),yseq(nmax),safe1,safe2, 4 redmax,redmin,tiny,scalmx parameter (safe1 = 0.25d0, safe2 = 0.7d0, redmax=1.0d-5, 1 redmin = 0.7d0, tiny = 1.0d-30, 2 scalmx = 0.1d0) c 2 scalmx = 0.5d0) data first/.true./, epsold/-1.0d0/ data nseq /2, 4, 6, 8, 10, 12, 14, 16, 18/ c..a new tolerance, so reinitialize if (eps .ne. epsold) then hnext = -1.0e29 xnew = -1.0e29 eps1 = safe1 * eps c..compute the work coefficients a_k a(1) = nseq(1) + 1 do k=1,kmaxx a(k+1) = a(k) + nseq(k+1) enddo c..compute alf(k,q) do iq=2,kmaxx do k=1,iq-1 alf(k,iq) = eps1**((a(k+1) - a(iq+1)) / 1 ((a(iq+1) - a(1) + 1.0d0) * (2*k + 1))) enddo enddo epsold = eps c..determine optimal row number for convergence do kopt=2,kmaxx-1 if (a(kopt+1) .gt. a(kopt)*alf(kopt-1,kopt)) go to 01 enddo 01 kmax = kopt end if c..save the starting values h = htry do i=1,nv ysav(i) = y(i) enddo c..a new stepsize or a new integration, re-establish the order window if (h .ne. hnext .or. x .ne. xnew) then first = .true. kopt = kmax end if reduct = .false. c..evaluate the sequence of modified midpoint integrations 02 do k=1,kmax xnew = x + h if (xnew .eq. x) stop 'stepsize too small in routine bsstep' call mmid(ysav,dydx,nv,x,h,nseq(k),yseq,derivs) xest = (h/nseq(k))**2 call pzextr(k,xest,yseq,y,yerr,nv) c..compute normalized error estimate if (k .ne. 1) then errmax = tiny do i=1,nv errmax = max(errmax,abs(yerr(i)/yscal(i))) enddo errmax = errmax/eps km = k - 1 err(km) = (errmax/safe1)**(1.0d0/(2*km+1)) end if c..in order window if (k .ne. 1 .and. (k .ge. kopt-1 .or. first)) then c..converged if (errmax .lt. 1.0) go to 04 c..possible step size reductions if (k .eq. kmax .or. k .eq. kopt + 1) then red = safe2/err(km) go to 3 else if (k .eq. kopt) then if (alf(kopt-1,kopt) .lt. err(km)) then red = 1.0d0/err(km) go to 03 end if else if (kopt .eq. kmax) then if (alf(km,kmax-1) .lt. err(km)) then red = alf(km,kmax-1) * safe2/err(km) go to 03 end if else if (alf(km,kopt) .lt. err(km)) then red = alf(km,kopt-1)/err(km) go to 03 end if end if c..back for another sequence of midpoint rules enddo c..reduce stepsize by at least redmin and at most redmax 03 red = min(red,redmin) red = max(red,redmax) h = h * red reduct = .true. go to 2 c..successful step; get optimal row for convergence and corresponding stepsize 04 x = xnew hdid = h first = .false. wrkmin = 1.0e35 do kk=1,km fact = max(err(kk),scalmx) work = fact * a(kk+1) if (work .lt. wrkmin) then scale = fact wrkmin = work kopt = kk + 1 end if enddo c..check for possible order increase, but not if stepsize was just reduced hnext = h/scale if (kopt .ge. k .and. kopt .ne. kmax .and. .not.reduct) then fact = max(scale/alf(kopt-1,kopt),scalmx) if (a(kopt+1)*fact .le. wrkmin) then hnext = h/fact kopt = kopt + 1 end if end if return end subroutine mmid(y,dydx,nvar,xs,htot,nstep,yout,derivs) implicit none save c..modified midpoint step. dependent variable vector y(1:nvar) and its c..derivative vector dydx(1:nvar) are input at xs. also input is htot, the c..total step to be made, and nstep, the number of substeps to be used. the c..output is returned as yout which need not be distinct from y; if it is c..distinct, y and dydx are returned unchanged. c..declare external derivs integer nvar,nstep,nmax,i,n parameter (nmax = 20) double precision y(nvar),dydx(nvar),xs,htot,yout(nvar), 1 ym(nmax),yn(nmax),h,x,h2,swap c..the step size on this trip h = htot/nstep do i=1,nvar ym(i) = y(i) yn(i) = y(i) + h*dydx(i) enddo c..use yout as temporary storage x = xs + h call derivs(x,yn,yout) h2 = 2.0d0 * h c..for the general step do n=2,nstep do i=1,nvar swap = ym(i) + h2*yout(i) ym(i) = yn(i) yn(i) = swap enddo x = x + h call derivs(x,yn,yout) enddo c..and the last step do i=1,nvar yout(i) = 0.5d0 * (ym(i) + yn(i) + h*yout(i)) enddo return end subroutine pzextr(iest,xest,yest,yz,dy,nv) implicit none save c..use polynomial extrapolation to evaluate nv functions at x=0 by fitting c..a polynomial to a sequence of estimates with progressively smaller values c..x=xest, and corresponding function vectors yest(1:nv). the call is number c..iest in the sequence of calls. extrapolated function values are output as c..yz(1:nv), and their estimated error is output as dy(1:nv) c..declare integer iest,nv,j,k1,nmax,imax parameter (nmax=50, imax=13) double precision xest,dy(nv),yest(nv),yz(nv),delta,f1,f2,q, 1 d(nmax),qcol(nmax,imax),x(imax) c..save current independent variables x(iest) = xest do j=1,nv dy(j) = yest(j) yz(j) = yest(j) enddo c..store first estimate in first column if (iest .eq. 1) then do j=1,nv qcol(j,1) = yest(j) enddo else do j=1,nv d(j) = yest(j) enddo do k1=1,iest-1 delta = 1.0d0/(x(iest-k1) - xest) f1 = xest * delta f2 = x(iest-k1) * delta c..propagate tableu 1 diagonal more do j=1,nv q = qcol(j,k1) qcol(j,k1) = dy(j) delta = d(j) - q dy(j) = f1*delta d(j) = f2*delta yz(j) = yz(j) + dy(j) enddo enddo do j=1,nv qcol(j,iest) = dy(j) enddo end if return end