c from C.T. Russell, Geophysical Coordinate Transformations" Cosmic
c Electrodynamics, 2, 1971, pp 184-196.
c
c Appendix 2. The Calculation of the Position of the Sun
c 
c G.D. Mead (private communication) has written a simple subroutine to calculate the position of the Sun in GEI coordinates.
c It is accurate for years 1901 through 2099, to within 0.006 deg. The input is the year, day of year and seconds of the day in
c UT. The output is Greenwich Mean Sideral Time in degrees, the ecliptic longitude, apparent right ascension and declination
c of the Sun in degrees. The listing of this program follows. We note that the cartesian coordinates of the vector from the Earth
c to the Sun are: 
c 
c   X = cos(SRASN) cos(SDEC)
c   Y = sin(SRASN) cos(SDEC)
c   Z = sin(SDEC)
c 
c 
  SUBROUTINE SUN(IYR, IDAY, SECS, GST, SLONG, SRASN, SDEC)
C PROGRAM TO CALCULATE SIDEREAL TIME AND POSITION OF THE SUN. 
C GOOD FOR YEARS 1901 THROUGH 2099. ACCURACY 0.006 DEGREE.
C INPUT IS IYR, IDAY (INTEGERS), AND SECS, DEFINING UN. TIME. 
C OUTPUT IS GREENWICH MEAN SIDEREAL TIME (GST) IN DEGREES,
C LONGITUDE ALONG ECLIPTIC (SLONG), AND APPARENT RIGHT ASCENSION
C AND DECLINATION (SRASN, SDEC) OF THE SUN, ALL IN DEGREES
  DATA RAD /57.29578/ 
  DOUBLE PRECISION DJ, FDAY 
  IF(IYR. LT. 1901. OR. IYR. GT. 2099) RETURN
  FDAY = SECS/86400
  DJ = 365* (IYR-1900) + (IYR-1901)/4 + IDAY + FDAY -0.5D0 
  T = DJ / 36525 
  VL = DMOD (279.696678 + 0.9856473354*DJ, 360.D0) 
  GST = DMOD (279.690983 + 0.9856473354*DJ + 360.*FDAY + 180., 360.D0)
  G = DMOD (358.475845 + 0.985600267*DJ, 360.D0) / RAD 
  SLONG = VL + (1.91946 -0.004789*T)*SIN(G) + 0.020094*SIN (2.*G) 
  OBLIQ = (23.45229 -0.0130125*T) / RAD 
  SLP = (SLONG -0.005686) / RAD 
  SIND = SIN (OBLIQ)*SIN (SLP) 
  COSD = SQRT(1.-SIND**2)
  SDEC = RAD * ATAN (SIND/COSD) 
  SRASN = 180. -RAD*ATAN2
  (COTAN (OBLIQ)*SIND/COSD, COS (SLP)/COSD) 
  RETURN 
  END