Subroutine ORBIT (ECCEN,OBLIQ,OMEGVP, DAY, * SIND,COSD,SUNDIS,SUNLON,SUNLAT,EQTIME) C**** C**** ORBIT receives orbital parameters and time of year, and returns C**** distance from Sun, declination angle, and Sun's overhead position. C**** Reference for following caculations is: V.M.Blanco and C**** S.W.McCuskey, 1961, "Basic Physics of the Solar System", pages C**** 135 - 151. Existence of Moon and heavenly bodies other than C**** Earth and Sun are ignored. Earth is assumed to be spherical. C**** C**** Program author: Gary L. Russell 2004/11/16 C**** Angles, longitude and latitude are measured in radians. C**** C**** Input: ECCEN = eccentricity of the orbital ellipse C**** OBLIQ = latitude of Tropic of Cancer C**** OMEGVP = longitude of perihelion (sometimes Pi is added) = C**** = spatial angle from vernal equinox to perihelion C**** with Sun as angle vertex C**** DAY = days measured since 2000 January 1, hour 0 C**** C**** Constants: EDAYzY = tropical year = Earth days per year = 365.2425 C**** VE2000 = days from 2000 January 1, hour 0 till vernal C**** equinox of year 2000 = 31 + 29 + 19 + 7.5/24 C**** C**** Intermediate quantities: C**** BSEMI = semi minor axis in units of semi major axis C**** PERIHE = perihelion in days since 2000 January 1, hour 0 C**** in its annual revolution about Sun 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 C**** eccentric circle (that circumscribes Earth's orbit) C**** from perihelion to point above (or below) Earth's C**** absisca (where absisca is directed from center of C**** eccentric circle to perihelion) C**** MA = mean anomaly = temporal angle from perihelion to C**** current time in units of 2*Pi per tropical year C**** TAofVE = TA(VE) = true anomaly of vernal equinox = - OMEGVP C**** EAofVE = EA(VE) = eccentric anomaly of vernal equinox C**** MAofVE = MA(VE) = mean anomaly of vernal equinox C**** SLNORO = longitude of Sun in Earth's nonrotating reference frame C**** VEQLON = longitude of Greenwich Meridion in Earth's nonrotating C**** reference frame at vernal equinox C**** ROTATE = change in longitude in Earth's nonrotating reference C**** frame from point's location on vernal equinox to its C**** current location where point is fixed on rotating Earth C**** SLMEAN = longitude of fictitious mean Sun in Earth's rotating C**** reference frame (normal longitude and latitude) C**** C**** Output: SIND = sine of declination angle = sin(SUNLAT) C**** COSD = cosine of the declination angle = cos(SUNLAT) C**** SUNDIS = distance to Sun in units of semi major axis C**** SUNLON = longitude of point on Earth directly beneath Sun C**** SUNLAT = latitude of point on Earth directly beneath Sun C**** EQTIME = Equation of Time = C**** = longitude of fictitious mean Sun minus SUNLON C**** C**** From the above reference: C**** (4-54): [1 - ECCEN*cos(EA)]*[1 + ECCEN*cos(TA)] = (1 - ECCEN^2) C**** (4-55): tan(TA/2) = sqrt[(1+ECCEN)/(1-ECCEN)]*tan(EA/2) C**** Yield: tan(EA) = sin(TA)*sqrt(1-ECCEN^2) / [cos(TA) + ECCEN] C**** or: tan(TA) = sin(EA)*sqrt(1-ECCEN^2) / [cos(EA) - ECCEN] C**** C Use C540, Only: EDAYzY,VE2000,TWOPI Implicit Real*8 (A-H,M-Z) Parameter (TWOPI=6.283185307179586477d0, * EDAYzY=365.2425d0, VE2000=79.3125) C**** C**** Determine EAofVE from geometry: tan(EA) = b*sin(TA) / [e+cos(TA)] C**** Determine MAofVE from 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 (BSEMI*Sin(TAofVE), ECCEN+Cos(TAofVE)) MAofVE = EAofVE - ECCEN*Sin(EAofVE) C PERIHE = VE2000 - MAofVE*EDAYzY/TWOPI MA = Modulo (TWOPI*(DAY-VE2000)/EDAYzY + MAofVE, TWOPI) 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) > 1d-10) GoTo 10 C**** C**** Calculate distance to Sun and true anomaly C**** SUNDIS = 1 - ECCEN*Cos(EA) TA = ATan2 (BSEMI*Sin(EA), Cos(EA)-ECCEN) C**** C**** Change reference frame to be nonrotating reference frame, angles C**** fixed according to stars, with Earth at center and positive x C**** axis be ray from Earth to Sun were Earth at vernal equinox, and C**** x-y plane be Earth's equatorial plane. Distance from current Sun C**** to this x axis is SUNDIS sin(TA-TAofVE). At vernal equinox, Sun C**** is located at (SUNDIS,0,0). At other times, Sun is located at: C**** C**** SUN = (SUNDIS cos(TA-TAofVE), C**** SUNDIS sin(TA-TAofVE) cos(OBLIQ), C**** SUNDIS sin(TA-TAofVE) sin(OBLIQ)) C**** SIND = Sin(TA-TAofVE) * Sin(OBLIQ) COSD = Sqrt (1 - SIND*SIND) SUNX = Cos(TA-TAofVE) SUNY = Sin(TA-TAofVE) * Cos(OBLIQ) SLNORO = ATan2 (SUNY,SUNX) C**** C**** Determine Sun location in Earth's rotating reference frame C**** (normal longitude and latitude) C**** VEQLON = TWOPI*VE2000 - TWOPI/2 + MAofVE - TAofVE ! modulo 2*Pi ROTATE = TWOPI*(DAY-VE2000)*(EDAYzY+1)/EDAYzY SUNLON = Modulo (SLNORO-ROTATE-VEQLON, TWOPI) If (SUNLON > TWOPI/2) SUNLON = SUNLON - TWOPI SUNLAT = ASin (Sin(TA-TAofVE)*Sin(OBLIQ)) C**** C**** Determine longitude of fictitious mean Sun C**** Calculate Equation of Time C**** SLMEAN = TWOPI/2 - TWOPI*(DAY-Floor(DAY)) EQTIME = Modulo (SLMEAN-SUNLON, TWOPI) If (EQTIME > TWOPI/2) EQTIME = EQTIME - TWOPI C**** Return End