ephem.f

Go to the documentation of this file.

      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