subroutine besselts(k, n, jdmax, gamma, u, type, annular_total, & magnitude, jd0, x, y, sind, cosd, mu, dd_dt, dmu_dt, l1, l2, & tanf1, tanf2) ! ! For a given lunation number k, returns the general circumstances ! of the solar eclipse that might occur at the New Moon correspond- ! ing to k, and n Besselian elements of the eclipse tabulated ! at 10-minute intervals centered on the time of the eclipse. (If ! there is no eclipse, type = 0 is returned, and no Besselian ! elements are computed.) ! ! B. Knapp, 2002-08-13 ! ! implicit none ! ! Input integer*4 k !Lunation number (k = 0 for new moon of !6 Jan 2000, and k is approximately !(year-2000)*12.3685, where year is year !and fraction of year). ! integer*4 n !Number of Besselian elements requested !(dimension of output arrays) ! ! Output: real*8 jdmax !Julian Date (UT) of maximum eclipse ! real*8 gamma !Least distance from the axis of the moon's !shadow to the center of the earth, in units !of the equatorial radius of the earth. ! real*8 u !Radius of the moon's umbral cone in the !fundamental plane in units of the !equatorial radius of the earth. ! integer*4 type !Type of eclipse: 0=none, 1=possibly visible !from an orbiting spacecraft, 2=partial, !3=non-central annular or total, 4=central ! integer*4 annular_total !0=neither annular nor total, 1=annular, !2=annular-total, 3=total ! real*8 magnitude !Magnitude of the eclipse ("fraction" of the !Sun's disk obscured by the Moon at the !point of maximum eclipse; magnitudes greater !than 1 indicate how much larger than the !Sun's disk is the area covered by the Moon) ! ! Besselian elements at 10-minute intervals beginning at jd0 (UT) real*8 jd0 !JD (UT) at which table values commence real*8 x(n), y(n) !Coordinates of shadow axis in fundamental !plane real*8 sind(n) !Sine of shadow axis declination real*8 cosd(n) !Cosine of shadow axis declination real*8 mu(n) !Greenwich hour angle of shadow axis real*8 dd_dt !Time rate of change of shadow axis !declination (radians/hour) real*8 dmu_dt !Time rate of change of Greenwich hour angle !of shadow axis (radians/hour) real*8 l1(n) !Penumbra radius (in Earth radii) real*8 l2(n) !Umbra radius (in Earth radii) real*8 tanf1 !Tangent of the penumbral cone vertex !halfangle real*8 tanf2 !Tangent of the umbral cone vertex halfangle ! ! Constants (IAU 1976, WGS84 values) real*8 au !Astronomical Unit (km) parameter (au=149597870.66d0) !2003:(au=149597870.691d0) real*8 r_earth !Equatorial radius of the Earth (km) parameter (r_earth=6378.137d0) real*8 r_sun !Radius of the Sun (Earth radii) parameter (r_sun=109.121087391622d0) !2003:(r_sun=109.122224500d0) real*8 r_moon !Radius of the Moon (Earth radii) parameter (r_moon=0.2725076d0) !IAU76!AApD3:(r_moon=0.272493d0) ! ! Misc. constants real*8 jd2000 parameter (jd2000 = 2441545.0d0) real*8 deg2rad parameter (deg2rad=3.14159265358979324d0/180.d0) ! ! External functions, subroutines (in ~knapp/src/astronomy) real*8 td2ut, ut2td, obliq, sidtim external solarecl, td2ut, ut2td, ephem, obliq, equecl, sphxyz, & sidtim ! ! Local variables integer*4 body, posnid, framid, coorid, status, t0, i, j real*8 jdemax, jdmax, jd, ag, pvsun(6), pvmoon(6), t, eps, & step, idir(3), jdir(3), kdir(3), kmag, sina, cosa, zi, & cscf1, tanf1i, c1, cscf2, tanf2i, c2, dt ! ! ! Get the general circumstances call solarecl(k, jdemax, gamma, u) ag = abs( gamma ) if ( ag .gt. 1.6433d0+u ) then type = 0 elseif ( ag .ge. 1.5433d0+u ) then type = 1 elseif ( ag .ge. 0.9972d0+abs(u) ) then type = 2 elseif ( ag .ge. 0.9972d0 ) then type = 3 else type = 4 endif ! if (type .lt. 3) then annular_total = 0 else if ( u .lt. 0 ) then annular_total = 3 elseif ( u .lt. 0.00464d0*sqrt(1.d0-ag**2) ) then annular_total = 2 else annular_total = 1 endif magnitude = (1.5433d0+u-ag)/(0.5461d0+2.d0*u) ! ! Get UT for maximum eclipse and t0 for the Besselian elements table jdmax = td2ut(jdemax) t0 = int(mod(jdmax,1.d0)*24.d0)-n/12 if (t0 .lt. 0) then jd0 = dint(jdmax-1) + (1.d0+dble(t0)/24.d0) else jd0 = dint(jdmax) + dble(t0)/24.d0 endif ! ! Is there an eclipse? if (type .eq. 0) return ! ! Compute the Besselian elements step = 10.d0/1440.d0 do i=1,n jd = jd0+(i-1)*step ! ! Get Sun's apparent position in ECI rectangular coordinates body = 0 !Sun posnid = 2 !apparent framid = 1 !equatorial coorid = 2 !rectangular call ephem(jd, body, posnid, framid, coorid, pvsun, status) ! ! Get Moon's apparent position in Ecliptic, spherical coordinates body = 10 !Moon posnid = 2 !apparent framid = 2 !ecliptic coorid = 1 !spherical call ephem(jd, body, posnid, framid, coorid, pvmoon, status) ! ! Apply the libration corrections pvmoon(1) = pvmoon(1) + 0.50d0/3600.d0 !dLambda pvmoon(2) = pvmoon(2) - 0.25d0/3600.d0 !dBeta ! ! Get the obliquity of the ecliptic t = (ut2td(jd)-jd2000)/365250.d0 eps = obliq(t) ! ! Transform moon to equatorial, rectangular coordinates do j=1,2 pvmoon(j) = pvmoon(j)*deg2rad pvmoon(j+3) = pvmoon(j+3)*deg2rad enddo call equecl(eps, -1, pvmoon, pvmoon) call sphxyz(1, pvmoon, pvmoon) ! ! Convert from AU, AU/day to km, km/sec, then to earth radii do j=1,3 pvsun(j) = pvsun(j)*au/r_earth pvmoon(j) = pvmoon(j)*au/r_earth pvsun(j+3) = pvsun(j+3)*au/86400.d0/r_earth pvmoon(j+3) = pvmoon(j+3)*au/86400.d0/r_earth enddo ! ! Get the direction vector of the shadow axis kmag = 0.d0 do j=1,3 kdir(j) = pvsun(j)-pvmoon(j) kmag = kmag+kdir(j)**2 enddo kmag = sqrt(kmag) do j=1,3 kdir(j) = kdir(j)/kmag enddo ! ! Need sines, cosines of the the right ascension, declination ! of the shadow axis position sind(i) = kdir(3) cosd(i) = sqrt(1.d0-sind(i)**2) if (cosd(i) .gt. 0.d0) then sina = kdir(2)/cosd(i) cosa = kdir(1)/cosd(i) else sina = 0.d0 cosa = 0.d0 endif ! ! Compute x, y coordinates of the shadow axis idir(1) = -sina idir(2) = cosa idir(3) = 0.d0 jdir(1) = -cosa*sind(i) jdir(2) = -sina*sind(i) jdir(3) = cosd(i) x(i) = 0.d0 y(i) = 0.d0 zi = 0.d0 do j=1,3 x(i) = x(i) + pvmoon(j)*idir(j) y(i) = y(i) + pvmoon(j)*jdir(j) zi = zi + pvmoon(j)*kdir(j) enddo ! ! For the Greenwich hour angle of the shadow axis, we need ! the sidereal time mu(i) = sidtim(jd) - atan2(kdir(2), kdir(1))/deg2rad if (mu(i) .lt. 0.d0) mu(i) = mu(i)+360.d0 if (mu(i) .gt. 360.d0) mu(i) = mu(i)-360.d0 ! ! Now compute the penumbra and umbra vertex angles and radii ! ! Penumbra cscf1 = kmag/(r_sun+r_moon) tanf1i = 1.d0/sqrt(cscf1**2-1.d0) c1 = zi+r_moon*cscf1 l1(i) = c1*tanf1i ! ! Umbra cscf2 = kmag/(r_sun-r_moon) tanf2i = 1.d0/sqrt(cscf2**2-1.d0) c2 = zi-r_moon*cscf2 l2(i) = c2*tanf2i ! ! Return the vertex halfangle tangents at the middle of the ! eclipse if (i .eq. n/2) then tanf1 = tanf1i tanf2 = tanf2i endif enddo ! ! Finally, the time rates of change of mu and d (radians/hour) dt = dble(n-1)/6.d0 dmu_dt = (mu(n)-mu(1))/dt*deg2rad dd_dt = (asin(sind(n))-asin(sind(1)))/dt return end