SUBROUTINE sunpos(jd,ra,dec,elng,dist)

      IMPLICIT real*8 (a-z)
      DATA twopi / 6.2831853071795862d0 /
     

c jd is input julian day (correct, not reduced)
c ra is output right ascension of the sun, in radians
c dec is output declination of the sun, in radians
c elng is output ecliptic longitude of earth, in radians
c dist is output earth-sun distance, in a.u.

c                       Compute components, in radians.

c n: time argument (the beginning of the tabulated ephemeris year)
c L: mean longitude of sun, corrected for aberration
c g: mean anomaly (angular distance from perihelion)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c COMMENT OUT ONE OF THE FOLLOWING TWO DECLARATION STATEMENTS ...
c 
c-This statement uncommented for calculation in degrees.

c     REAL*8
c    &   f   / 1.745329252d-2 / , cir / 360.0 /,
c    &   n0  / 2451545.0 / , 
c    &   L0  / 280.472   / , L1  / 0.9856474 / ,
c    &   g0  / 357.528   / , g1  / 0.9856003 / ,
c    &   el1 / 1.915     / , el2 / 0.020     / ,
c    &   ob0 / 23.439    / , ob1 / 4.d-7     / ,
c    &   R0  / 1.00014   / , R1  / 0.01671   / , R2  / 1.4d-4 /

c-This statement uncommented for calculation in radians.

      REAL*8
     &   f   / 1. / , cir / 6.2831853071795862 /,
     &   n0  / 2451545.0   / , 
     &   L0  / 4.8951599   / , L1  / 1.7202792d-2 / ,
     &   g0  / 6.2400407   / , g1  / 1.7201970d-2 / ,
     &   el1 / 3.3423d-2   / , el2 / 3.4907d-4    / ,
     &   ob0 / 4.090877d-1 / , ob1 / 6.98d-9      / ,
     &   R0  / 1.00014     / , R1  / 0.01671   / , R2  / 1.4d-4 /

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      cirrange(x,y)= mod(mod(x,y)+y,y)

      n = jd - n0
      L = L0 + L1*n 
      L= cirrange(L, cir)
      g = g0 + g1*n
      g = cirrange(g, cir)
      elng = L + el1*sin(g*f) + el2*sin(2.*g*f)
      elng = cirrange(elng, cir)
      oblq = ob0 - ob1*n
      t = tan(oblq * f / 2.) ** 2

c                       Compute celestial coordinates.

      ra = elng - (t * sin(2. * elng * f)) + 
     &                (0.5 * (t**2) * sin(4. * elng * f))
      ra=cirrange(ra, cir)

      dec = asin(sin(oblq*f) * sin(elng*f)) / f

      dist= R0 - (R1 * cos(g * f)) - (R2 * cos(2. * g * f))

      return
      end