subroutine sun(y,m,DD,UT,lon,lat,RA,Dec,LST,Az,El,mjd,day)

  implicit none

  integer y                         !Year
  integer m                         !Month
  integer DD                        !Day
  integer mjd                       !Modified Julian Date
  real UT                           !UTC in hours
  real RA,Dec                       !RA and Dec of sun

! NB: Double caps here are single caps in the writeup.

! Orbital elements of the Sun (also N=0, i=0, a=1):
  real w                            !Argument of perihelion
  real e                            !Eccentricity
  real MM                           !Mean anomaly
  real Ls                           !Mean longitude

! Other standard variables:
  real v                            !True anomaly
  real EE                           !Eccentric anomaly
  real ecl                          !Obliquity of the ecliptic
  real d                            !Ephemeris time argument in days
  real r                            !Distance to sun, AU
  real xv,yv                        !x and y coords in ecliptic
  real lonsun                       !Ecliptic long and lat of sun
! Ecliptic coords of sun (geocentric)
  real xs,ys
! Equatorial coords of sun (geocentric)
  real xe,ye,ze
  real lon,lat
  real GMST0,LST,HA
  real xx,yy,zz
  real xhor,yhor,zhor
  real Az,El

  real day
  real rad
  data rad/57.2957795/

! Time in days, with Jan 0, 2000 equal to 0.0:
  d=367*y - 7*(y+(m+9)/12)/4 + 275*m/9 + DD - 730530 + UT/24.0
  mjd=d + 51543
  ecl = 23.4393 - 3.563e-7 * d

! Compute updated orbital elements for Sun:
  w = 282.9404 + 4.70935e-5 * d
  e = 0.016709 - 1.151e-9 * d
  MM = mod(356.0470d0 + 0.9856002585d0 * d + 360000.d0,360.d0)
  Ls = mod(w+MM+720.0,360.0)

  EE = MM + e*rad*sin(MM/rad) * (1.0 + e*cos(M/rad))
  EE = EE - (EE - e*rad*sin(EE/rad)-MM) / (1.0 - e*cos(EE/rad))

  xv = cos(EE/rad) - e
  yv = sqrt(1.0-e*e) * sin(EE/rad)
  v = rad*atan2(yv,xv)
  r = sqrt(xv*xv + yv*yv)
  lonsun = mod(v + w + 720.0,360.0)
! Ecliptic coordinates of sun (rectangular):
  xs = r * cos(lonsun/rad)
  ys = r * sin(lonsun/rad)

! Equatorial coordinates of sun (rectangular):
  xe = xs
  ye = ys * cos(ecl/rad)
  ze = ys * sin(ecl/rad)

! RA and Dec in degrees:
  RA = rad*atan2(ye,xe)
  Dec = rad*atan2(ze,sqrt(xe*xe + ye*ye))

  GMST0 = (Ls + 180.0)/15.0
  LST = mod(GMST0+UT+lon/15.0+48.0,24.0)    !LST in hours
  HA = 15.0*LST - RA                        !HA in degrees
  xx = cos(HA/rad)*cos(Dec/rad)
  yy = sin(HA/rad)*cos(Dec/rad)
  zz =             sin(Dec/rad)
  xhor = xx*sin(lat/rad) - zz*cos(lat/rad)
  yhor = yy
  zhor = xx*cos(lat/rad) + zz*sin(lat/rad)
  Az = mod(rad*atan2(yhor,xhor) + 180.0 + 360.0,360.0)
  El = rad*asin(zhor)
  day=d-1.5

  return
end subroutine sun