delm.f

Go to the documentation of this file.

      SUBROUTINE delm(RE,GE,AJ2,X,DX)
C  VERSION 4/17/87
C    COMPUTES OSCULATING - MEAN ELEMENTS WHEN GIVEN MEAN ELEMENTS
C  INPUT
C    RE      = RADIUS OF PLANET (KM)
C    GE      = GRAVITATIONAL CONSTANT * MASS OF PLANET (KM**3/SEC**2)
C    AJ2     = J2, = -C20
C    X       = MEAN ELEMENTS
C    X(1)    = A, SEMI-MAJOR AXIS (KM)
C    X(2)    = E, ECCENTRICITY
C    X(3)    = I, INCLINATION (RAD)
C    X(4)    = NODE, LONGITUDE OF ASCENDING NODE (RAD)
C    X(5)    = W, ARGUMENT OF PERIAPSIS (RAD)
C    X(6)    = M, MEAN ANOMALY (RAD)
C  OUTPUT
C    DX      = ARRAY OF 6 VALUES = OSCULATING - MEAN ELEMENTS
C  CALL SUBROUTINES
C    ANOMLY
C  REFERENCES
C    JPL IOM 312/85.2-927, 23 JANUARY 1985, BY C. UPHOFF
C    JPL EM 312/87-153, 20 APRIL 1987
C  ANALYSIS
C    JOHNNY H. KWOK - JPL
C  PROGRAMMER
C    JOHNNY H. KWOK - JPL
C  MODIFICATIONS
C    NONE
C  COMMENTS
C    THIS PROGRAM USES AN ALGORITHM DERIVED BY C. UPHOFF (REF) WHICH
C    USES A COMBINATION OF KOZAI'S AND IZSAK'S THEORY
C
      IMPLICIT DOUBLE PRECISION (a-h,o-z)
      dimension x(6),dx(6)
      DATA half,one,two,three,four,five,six,seven
     1  /0.5d0,1.d0,2.d0,3.d0,4.d0,5.d0,6.d0,7.d0/
      DATA f14,f34,f32,f18/0.25d0,0.75d0,1.5d0,0.125d0/
      DATA tpi/6.283185307179586d0/
      a=x(1)
      e=x(2)
      ai=x(3)
      an=x(4)
      w=x(5)
      am=x(6)
      CALL anomly(2,am,e,f)
C
C *** AM AND F HAVE TO BE IN THE SAME QUADRANT
C
      f=dmod(f+tpi,tpi)
      am=dmod(am+tpi,tpi)
      si=dsin(ai)
      ci=dcos(ai)
      ti=si/ci
      si2=si*si
      ci2=ci*ci
      sf=dsin(f)
      cf=dcos(f)
      s2f=dsin(two*f)
      u=f+w
      e2=e*e
      esf=e*sf
      d1=one-e2
      d2=dsqrt(d1)
      d3=e*cf
      d4=one+d3
      d42=d4*d4
      d5=one+d2
      d6=(three*ci2-one)/d5
      p=a*d1
      d7=dsqrt(ge/p)
      r=p/d4
      rdot=d7*esf
      twou=two*u
      twow=two*w
      s2u=dsin(twou)
      c2u=dcos(twou)
      sf2w=dsin(f+twow)
      d8=three*f+twow
      s3f2w=dsin(d8)
      cf2w=dcos(f+twow)
      c3f2w=dcos(d8)
      q1=aj2*(re/p)**2
      di=f34*q1*si*ci*(c2u+e*cf2w+e/three*c3f2w)
      dp=two*p*ti*di
      dummy1=f-am+esf-half*s2u-half*e*sf2w-e*s3f2w/six
      dn=-f32*q1*ci*dummy1
      dr=-f14*p*q1*((three*ci2-one)*(two*d2/d4+d3/d5+one)-si2*c2u)
      drdot=f14*d7*q1*(d6*esf*(d2*d5+d42)-two*si2*d42*s2u)
      du=-f18*q1*(six*(one-five*ci2)*(f-am)+four*esf*((one-six*ci2)
     1  -d6)-d6*e2*s2f+two*(five*ci2-two)*e*sf2w+(seven*ci2-one)*s2u
     2  +two*ci2*e*s3f2w)
      pnw=p+dp
      ainw=ai+di
      annw=an+dn
      rnw=r+dr
      rdotnw=rdot+drdot
      unw=u+du
      aa=pnw/rnw-one
      bb=dsqrt(pnw/ge)*rdotnw
      enw2=aa*aa+bb*bb
      enw=dsqrt(enw2)
      fnw=datan2(bb,aa)
      anw=pnw/(one-enw2)
      wnw=unw-fnw
      CALL anomly(5,fnw,enw,amnw)
      dx(1)=anw-a
      dx(2)=enw-e
      dx(3)=ainw-ai
      dx(4)=annw-an
      dx(5)=wnw-w
      dx(6)=amnw-am
      RETURN
      END