SUBROUTINE ORBIT (OBLIQD,ECCEN,OMEGVP, DAY, DISTQ,SIND,COSD)
C****
C**** ORBIT receives the orbital parameters and time of year, and
C**** returns the distance from the Sun and its declination angle.
C**** The reference for the following caculations is: V.M.Blanco
C**** and S.W.McCuskey, 1961, "Basic Physics of the Solar System",
C**** pages 135 - 151.
C****
C**** Input: OBLIQD = latitude of Tropic of Cancer in degrees
C****        ECCEN  = eccentricity of the orbital ellipse
C****        OMEGVP = longitude of perihelion =
C****               = spatial angle from vernal equinox to perihelion
C****                 in radians with Sun as angle vertex
C****        DAY    = day of the year in days; 0 = Jan 1, hour 0
C****
C**** Constants: EDAYzY = Earth days per year = 365
C****            VEREQX = occurence of vernal equinox = day 79 = Mar 21
C****
C**** Intermediate quantities:
C****    BSEMI = semi minor axis in units of the semi major axis
C****   TAofVE = TA(VE) = true anomaly of the vernal equinox = -OMEGVP
C****   EAofVE = EA(VE) = eccentric anomaly of the vernal equinox
C****   MAofVE = MA(VE) = mean anomaly of the vernal equinox
C****   PERIHE = perihelion in days
C****       TA = true anomaly = spatial angle from perihelion to
C****            current location with Sun as angle vertex
C****       EA = eccentric anomaly = spatial angle measured along the
C****            eccentric circle of the Earth's orbit from perihelion
C****            to the Earth's absisca (where the absisca is directed
C****            from the center of the eccentric circle to perihelion)
C****       MA = mean anomaly = temporal angle from perihelion to
C****          = current time in units of 2*pi per tropical day
C****     DIST = distance to the Sun in units of the semi major axis
C****
C**** Output: DISTQ = square of the distance to the Sun in units of
C****                 the semi major axis
C****          SIND = sine of the declination angle
C****          COSD = cosine of the declination angle
C****
      IMPLICIT REAL*8 (A-Z)
      PARAMETER (TWOPI=6.283185307179586477d0, EDAYzY=365., VEREQX=79.)
      OBLIQ  = OBLIQD*(TWOPI/360.)
C****
C**** Determine EAofVE from geometry:  tan(EA) = b*sin(TA)/[e+cos(TA)]
C**** Determine MAofVE using Kepler's equation:  MA = EA - e*sin(EA)
C**** Determine MA knowing time from vernal equinox to current day.
C****
      BSEMI  = SQRT(1.-ECCEN*ECCEN)
      TAofVE = - OMEGVP
      EAofVE = ATAN2(SIN(TAofVE)*BSEMI,COS(TAofVE)+ECCEN)
      MAofVE = EAofVE - ECCEN*SIN(EAofVE)
C     PERIHE = VEREQX - MAofVE*EDAYzY/TWOPI
      MA     = (DAY-VEREQX)*TWOPI/EDAYzY + MAofVE
C****
C**** Numerically invert Kepler's equation:  MA = EA - e*sin(EA)
C****
      EA  = MA + ECCEN*(SIN(MA)+ECCEN*SIN(2.*MA)/2.)
   10 dEA = (MA-EA+ECCEN*SIN(EA)) / (1.-ECCEN*COS(EA))
      EA  = EA + dEA
      IF(ABS(dEA).gt.1.d-10)  GO TO 10
C****
C**** Calculate the distance to the Sun and true anomaly
C****
      COSEA = COS(EA)
      SINEA = SIN(EA)
      DISTQ = (1.-ECCEN*COSEA)*(1.-ECCEN*COSEA)
      TA    = ATAN2(SINEA*BSEMI,COSEA-ECCEN)
C****
C**** Change the reference frame to be a nonrotating reference frame
C**** with the Earth at the center, the negative x axis to be the ray
C**** from the Earth to the Sun were the Earth at vernal equinox, and
C**** the x-y plane be the Earth's equatorial plane.  The distance
C**** from the current Sun to this x axis is:  DIST sin(TA-TAofVE).
C**** The Sun is located at:
C****
C**** SUN = (-DIST cos(TA-TAofVE),
C****        -DIST sin(TA-TAofVE) cos(OBLIQ),
C****         DIST sin(TA-TAofVE) sin(OBLIQ))
C****
C**** At Vernal Equinox:  mod(TA-TAofVE,2*pi) = 0,  SUN = (-DIST,0,0)
C****
      SIND   = SIN(TA-TAofVE)*SIN(OBLIQ)
      COSD   = SQRT(1.-SIND*SIND)
      SUNX   = -COS(TA-TAofVE)
      SUNY   = -SIN(TA-TAofVE)*COS(OBLIQ)
C****
      RETURN
      END