NASA Ocean Color

ocssw V2022
       SUBROUTINE ephem(ISUN,IMOON,ISRP,T,ES,EM)
 C  VERSION OF 1/31/87
 C  PURPOSE
 C    COMPUTES THE ORBITAL ELEMENTS OF A SUN AND/OR A MOON
 C  INPUT ARGUMENTS
 C    ISUN   = SEE SUBROUTINE ASAP
 C    IMOON  = SEE SUBROUTINE ASAP
 C    ISRP   = SEE SUBROUTINE ASAP
 C    T      = CURRENT JULIAN TIME (SEC)
 C    ES(3)  = INCLINATION OF THE SUN (RAD)
 C  OUTPUT ARGUMENTS
 C    ES(6)  = MEAN ANOMALY OF THE SUN AT T (RAD)
 C    ES(7)  = MEAN MOTION OF THE SUN (RAD/SEC)
 C    EM     = SEVEN ORBITAL ELEMENTS OF THE MOON AT T IN EARTH MEAN
 C             EQUATOR AND EQUINOX OF EPOCH
 C      (1)  = SEMI-MAJOR AXIS (KM)
 C      (2)  = ECCENTRICITY
 C      (3)  = INCLINATION (RAD)
 C      (4)  = LONGITUDE OF ASCENDING NODE (RAD)
 C      (5)  = ARGUMENT OF PERIAPSIS (RAD)
 C      (6)  = MEAN ANOMALY AT T (RAD)
 C      (7)  = MEAN MOTION (RAD/SEC)
 C  CALL SUBROUTINES
 C    NONE
 C  REFERENCES
 C    JPL EM 312/87-153, 20 APRIL 1987
 C    EXPLANATORY SUPPLEMENT TO THE ASTRONOMICAL EPHEMERIS AND THE
 C      AMERICAN EPHEMERIS AND NAUTICAL ALMANAC
 C  ANALYSIS
 C    JOHNNY H. KWOK - JPL
 C  PROGRAMMER
 C    JOHNNY H. KWOK - JPL
 C  MODIFICATIONS
 C    NONE
 C  COMMENTS
 C    FOR MERCURY, VENUS, AND MARS, THERE IS NO LUNI-PERTURBATION, USER
 C    SHOULD INPUT THE SUN'S ORBITAL ELEMENTS RELATIVE TO PLANET EQUATOR
 C    OF EPOCH.  THE INERTIAL X-AXIS COULD BE PLANET EQUINOX OR THE
 C    DIRECTION OF THE PRIME MERIDIAN AT SOME EPOCH SUCH AS DEFINED BY
 C    THE IAU.  FOR EARTH, THE ORBITAL ELEMENTS OF THE MOON VARIES TOO
 C    MUCH TO USE TWO-BODY EPHEMERIS.  USER MAY USE THE MEAN LUNAR
 C    EPHEMERIS BUILT INTO THE PROGRAM.  THE SHORT PERIOD VARIATIONS OF
 C    THE MOON MAY INTRODUCE AN ERROR UP TO TWO DEGREES.  AS FOR THE SUN,
 C    THE MEAN ORBITAL ELEMENTS OF THE SUN RELATIVE TO EARTH EQUATOR AND
 C    EQUINOX ARE BUILT IN ALSO.  BUT SINCE THEY DON'T CHANGE AS RAPIDLY
 C    AS THOSE OF THE MOON, THEY ARE INITIALIZE ELSEWHERE AND THEN HELD
 C    FIXED FOR THE DURATION OF THE TRAJECTORY PROPAGATION
 C
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
       dimension es(7),em(7)
       DATA rjd/2.41502d6/
       DATA pi,tpi/3.141592653589793d0,6.283185307179586d0/
       DATA ans,anm/1.720196977d-2,.22997150294101d0/
       DATA dts/8.64d4/
       dt=t/dts-rjd
       dt1=dt*1.d-4
       dt2=dt1*dt1
       dt3=dt1*dt2
 C
 C  SUN RELATIVE TO EARTH IN MEAN EQUATOR OF DATE
 C
       IF (isun.EQ.1.OR.isrp.EQ.1) THEN
       es(6)=dmod(6.256583575d0+ans*dt
      1  -1.9548d-7*dt2-1.22d-9*dt3,tpi)
       es(7)=ans/dts
       ENDIF
 C
 C  MOON IN EARTH MEAN EQUINOX OF DATE
 C
       IF (imoon.EQ.1) THEN
       em(1)=3.844d5
       em(2)=5.4900489d-2
       em(3)=8.980410805d-2
       em(4)=4.523601515d0-9.242202942d-4*dt+2.71748d-6*dt2+8.73d-10*dt3
       em(5)=5.835151539d0+1.944368001d-3*dt-1.35071d-5*dt2-4.538d-9*dt3
      1  -em(4)
       em(6)=4.719966572d0+anm*dt-1.4835d-6*dt2
      1  +6.80678d-10*dt3-em(4)-em(5)
       DO 130 i=4,6
   130 em(i)=dmod(em(i),tpi)
       em(7)=anm/dts
 C
 C  CONVERT TO EARTH MEAN EQUATOR OF DATE
 C
       se=dsin(es(3))
       ce=dcos(es(3))
       si=dsin(em(3))
       ci=dcos(em(3))
       so=dsin(em(4))
       co=dcos(em(4))
       ca=si*se*co-ce*ci
       alpha=dacos(ca)
       em(3)=pi-alpha
       sa=dsin(alpha)
       sw=so*si/sa
       sd=so*se/sa
       sso=dsign(1.d0,so)
       swso=sw*so*sso
       cwso=(sw*ce*co+sd*ci)*sso
       em(4)=datan2(swso,cwso)
       sdso=sd*so*sso
       cdso=(sd*co*ci+sw*ce)*sso
       del=datan2(sdso,cdso)
       em(5)=em(5)+del
       ENDIF
       RETURN
       END
#define pi
Definition: vincenty.c:23
dmod
double dmod(double a, double p)
Description:
Definition: nav.c:23
ephem
subroutine ephem(ISUN, IMOON, ISRP, T, ES, EM)
Definition: ephem.f:2