      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.
*
*	Routine due to G. D. Mead, privately communicated to C. T. Russell.
*	Published in Russell, C. T., Geophysical Coordinate Transformations,
*	Cosmical Electrodynamics, vol. 2, pp. 184-196, 1971.
*
      REAL*8 DJ, FDAY
*
      DATA RAD /57.29578/
*
      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.0D0)
      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.005686 * T) / RAD
      SLP = (SLONG - 0.005686) / RAD
      SIN_D = SIN (OBLIQ) * SIN (SLP)	! SIN_D, COS_D renamed because of
      COS_D = SQRT (1. - SIN_D ** 2)	! SIND, COSD functions in FORTRAN 77
      SDEC = RAD * ATAN(SIN_D / COS_D)
      SRASN = 180. - RAD * ATAN2 ((1.0 / TAN (OBLIQ)) * SIN_D / COS_D,
     &                            - COS (SLP) / COS_D)
      RETURN
      END
