subroutine SolarEcl( k, jde, gamma, u )
c
c     Compute the general circumstances for a solar eclipse at a
c     given new moon.  Uses algorithm of J. Meeus (Astronomical
c     Algorithms, Willman-Bell, 1991).
c
c     B. Knapp, 94.03.08
c
C
C     RCS DATA
C     
C     $Header$
C     
C     $Log$
C
C
      implicit none
c
c     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).
c
c     Output:
      real*8      jde      !Julian Ephemeris Date of new moon
c
      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.
c
      real*8      u        !Radius of the moon umbral cone in the
                           !fundamental plane in units of the equatorial
                           !radius of the earth.
!     Constants:
      real*8 PI,D2R
      parameter( PI=3.14159265358979324d0, D2R=PI/180.d0 )
!
!     ELP Arguments (polynomial coefficients, arc seconds)
      real*8  M_C(5)
     &/ 1287104.79306d0,  1295965810.4740d0,  -55.29d0,  0.147d0,
     &  0.0d0    /
      real*8 MP_C(5)     
     &/  485868.28096d0, 17179159234.7280d0, 3238.93d0, 51.651d0,
     & -2.4470d0 /
      real*8  F_C(5)
     &/  335779.55755d0, 17395272630.9830d0,-1225.05d0, -1.021d0,
     &  0.0417d0 /
      real*8 OM_C(5)
     &/  450160.39816d0,    -69628902.656d0,  747.42d0,  7.702d0,
     & -0.5939d0 /
c
c     Local variables
      real*8 tk, jdm, t, m, mp, f, om, e, f1, a1, p, q, w
c
c     Determine jdm (mean new moon) and the time argument t
      tk = k/1236.85d0
      jdm = 2451550.09765d0 + 29.530588853d0*k +
     &  (1.337d-4 + (-1.5d-7 + 7.3d-10*tk)*tk)*tk**2
      t = (jdm - 2451545.0d0)/365250.d0
!
!     Compute the moon angles using t
      m  =  M_C(1)+( M_C(2)+( M_C(3)+( M_C(4)+ M_C(5)*t)*t)*t)*t
      m  = mod(  m/3600.d0, 360.d0 )*D2R
      mp = MP_C(1)+(MP_C(2)+(MP_C(3)+(MP_C(4)+MP_C(5)*t)*t)*t)*t
      mp = mod( mp/3600.d0, 360.d0 )*D2R
      f  =  F_C(1)+( F_C(2)+( F_C(3)+( F_C(4)+ F_C(5)*t)*t)*t)*t
      f  = mod(  f/3600.d0, 360.d0 )*D2R
      om = OM_C(1)+(OM_C(2)+(OM_C(3)+(OM_C(4)+OM_C(5)*t)*t)*t)*t
      om = mod( om/3600.d0, 360.d0 )*D2R
c
      e  = 1.0d0 + (-2.516d-2 - 7.4d-4*t)*t
      f1 = f - 0.02665d0*D2R * sin( om )
      a1 = (299.77d0 + 0.107408d0*k - 0.009173d0*tk**2)*D2R
c
      jde = jdm
     &          -0.4075d0     * sin( mp )
     &          +0.1721d0 * e * sin( m )
     &          +0.0161d0     * sin( 2*mp )
     &          -0.0097d0     * sin( 2*f1 )
     &          +0.0073d0 * e * sin( mp-m )
     &          -0.0050d0 * e * sin( mp+m )
     &          -0.0023d0     * sin( mp-2*f1 )
     &          +0.0021d0 * e * sin( 2*m )
     &          +0.0012d0     * sin( mp+2*f1 )
     &          +0.0006d0 * e * sin( 2*mp+m )
     &          -0.0004d0     * sin( 3*mp )
     &          -0.0003d0 * e * sin( m+2*f1 )
     &          +0.0003d0     * sin( a1 )
     &          -0.0002d0 * e * sin( m-2*f1 )
     &          -0.0002d0 * e * sin( 2*mp-m )
     &          -0.0002d0     * sin( om )
c
      p =
     &          +0.2070d0 * e * sin( m )
     &          +0.0024d0 * e * sin( 2*m )
     &          -0.0392d0     * sin( mp )
     &          +0.0116d0     * sin( 2*mp )
     &          -0.0073d0 * e * sin( mp+m )
     &          +0.0067d0 * e * sin( mp-m )
     &          +0.0118d0     * sin( 2*f1 )
c
      q = 5.2207d0
     &          -0.0048d0 * e * cos( m )
     &          +0.0020d0 * e * cos( 2*m )
     &          -0.3299d0     * cos( mp )
     &          -0.0060d0 * e * cos( mp+m )
     &          +0.0041d0 * e * cos( mp-m )
c
      w = abs( cos( f1 ) )
      gamma = (p*cos( f1 ) + q*sin( f1 ))*(1.d0-0.0048d0*w)
      u = 0.0059d0
     &          +0.0046d0 * e * cos( m )
     &          -0.0182d0     * cos( mp )
     &          +0.0004d0     * cos( 2*mp )
     &          -0.0005d0     * cos( m+mp )
c
      return
      end