eqnox.f

Go to the documentation of this file.

      SUBROUTINE eqnox(X,GM,Y)  
C  VERSION OF 4/1/85
C  PURPOSE   
C    COMPUTES THE EQUINOCTIAL ELEMENTS A,H,K,P,Q,LAMDA AND 
C    CLASSICAL ELEMENTS E,I,NODE,W,M AND OTHER MISC VARIABLES
C  INPUT 
C    X      = 6-D CARTESIAN COORD X,Y,Z,XD,YD,ZD (KM,KM/SEC)
C    GM     = GRAVITATION CONSTANT * MASS OF PLANET (KM**3/SEC**2)
C  OUTPUT
C    Y(1)   = A, SEMI-MAJOR AXIS (KM)
C    Y(2)   = H, E * SIN(W + NODE) 
C    Y(3)   = K, E * COS(W + NODE) 
C    Y(4)   = P, TAN(I/2) * SIN(NODE)  
C    Y(5)   = Q, TAN(I/2) * COS(NODE)  
C    Y(6)   = LAMDA, M + NODE + W (RAD)
C    Y(7)   = E, ECCENTRICITY
C    Y(8)   = I, INCLINATION (RAD)  
C    Y(9)   = NODE, LONGITUDE OF ASCENDING NODE (RAD)
C    Y(10)  = W, ARGUMENT OF PERIAPSIS (RAD)
C    Y(11)  = MA, MEAN ANOMALY (RAD)
C    Y(12)  = TA, TRUE ANOMALY (RAD)
C    Y(13)  = EA, ECCENTRIC ANOMALY (RAD)
C    Y(14)  = L, TRUE LONGITUDE (RAD)
C    Y(15)  = F, ECCENTRIC LONGITUDE (RAD)
C    Y(16)  = R, RADIUS (KM)
C    Y(17)  = V, VELOCITY (KM/SEC)
C    Y(18)  = ORBITAL PERIOD (SEC)
C  CALL SUBROUTINES
C    NONE
C  REFERENCES
C    JPL EM 312/87-153, 20 APRIL 1987
C    AIAA #75-9, ON THE FORMULATION OF THE GRAVITATIONAL POTENTIAL 
C    IN TERMS OF EQUINOCTIAL VARIABLES, CEFOLA AND BROUCKE, 1975
C  ANALYSIS
C    JOHNNY H. KWOK - JPL  
C  PROGRAMMER
C    JOHNNY H. KWOK - JPL  
C  PROGRAM MODIFICATIONS 
C    NONE  
C  COMMENTS  
C    THIS PROGRAM RECYCLES STORAGE DUMMY VARIABLES.  SO BE CAREFUL
C    WHEN MODIFYING THIS ROUTINE.  THIS PROGRAM DOES NOT WORK FOR 
C    PARABOLIC OR HYPERBOLIC OR RETROGRADE EQUATORIAL ORBITS.  
C    IF EITHER INCLINATION OR ECCENTRICITY IS NEAR ZERO, THE   
C    CLASSICAL ORBITAL ELEMENTS MAY NOT BE DEFINED AND THE USER
C    SHOULD INTERPRET ACCORDINGLY. 
C   
      IMPLICIT DOUBLE PRECISION (a-h,o-z)   
      dimension x(6),y(18)
      dimension v1(3),v2(3),v3(3)   
      DATA zero,one,two/0.d0,1.d0,2.d0/
      DATA tpi/6.283185307179586d0/
C *** D2 = V2   
C *** D3 = RDV  
      y(16)=dsqrt(x(1)**2+x(2)**2+x(3)**2)  
      d2=x(4)**2+x(5)**2+x(6)**2
      d3=x(1)*x(4)+x(2)*x(5)+x(3)*x(6)  
      y(17)=dsqrt(d2)   
      y(1)=one/(two/y(16)-d2/gm)
      y(18)=tpi/dsqrt(gm/y(1)**3)
C *** D4 = V2/GM - 1/R  
C *** D5 = RDV/GM   
      d4=d2/gm-one/y(16)
      d5=d3/gm  
C *** V1 = ECCENTRICITY VECTOR  
      DO 10 i=1,3   
   10 v1(i)=d4*x(i)-d5*x(i+3)   
C *** V2 = VECTOR W = R CROSS V, ANGULAR MOMENTUM VECTOR
      v2(1)=x(2)*x(6)-x(3)*x(5) 
      v2(2)=x(3)*x(4)-x(1)*x(6) 
      v2(3)=x(1)*x(5)-x(2)*x(4) 
C *** D1 = MAGNITUDE OF W   
      d1=dsqrt(v2(1)**2+v2(2)**2+v2(3)**2)  
      DO 20 i=1,3   
   20 v2(i)=v2(i)/d1
C *** D1 = 1 + WZ   
      d1 = one+v2(3)
      y(4)=v2(1)/d1 
      y(5)=-v2(2)/d1
C *** D1 = P**2 
C *** D2 = Q**2 
      d1=y(4)*y(4)  
      d2=y(5)*y(5)  
C *** D3
C *** D4
      d3=one+d1+d2  
      d4=two*y(4)*y(5)  
C *** V2 = UNIT VECTOR G
C *** V3 = UNIT VECTOR F
      v2(1)=d4  
      v2(2)=one+d1-d2   
      v2(3)=two*y(5)
      v3(1)=one-d1+d2   
      v3(2)=d4  
      v3(3)=-two*y(4)   
      DO 30 i=1,3   
      v2(i)=v2(i)/d3
   30 v3(i)=v3(i)/d3
      y(2)=v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3)  
      y(3)=v1(1)*v3(1)+v1(2)*v3(2)+v1(3)*v3(3)  
C *** D1 = X1   
C *** D2 = Y1   
      d1=x(1)*v3(1)+x(2)*v3(2)+x(3)*v3(3)   
      d2=x(1)*v2(1)+x(2)*v2(2)+x(3)*v2(3)
      y(14)=datan2(d2,d1)
C *** D3 = DSQRT(1-H*H-K*K) 
C *** D4 = A * DSQRT(1-H*H-K*K) 
C *** D5 = BETA = 1/(1+A*DSQRT(1-H*H-K*K))  
      d3=dsqrt(one-y(2)**2-y(3)**2) 
      d4=y(1)*d3
      d5=one/(one+d3)   
C *** D6 = COS(F)   
C *** D7 = SIN(F)   
      d6=y(3)+((one-y(3)*y(3)*d5)*d1-y(2)*y(3)*d5*d2)/d4
      d7=y(2)+((one-y(2)*y(2)*d5)*d2-y(2)*y(3)*d5*d1)/d4
      y(15)=datan2(d7,d6)
      y(6)=y(15)-y(3)*d7+y(2)*d6
      y(7)=dsqrt(y(2)*y(2)+y(3)*y(3))   
      y(8)=two*datan(dsqrt(y(4)*y(4)+y(5)*y(5)))
      IF (y(8).NE.zero) THEN
         y(9)=datan2(y(4),y(5)) 
      ELSE  
         y(9)=zero  
      ENDIF 
      IF (y(7).NE.zero) THEN
         y(10)=datan2(y(2),y(3))-y(9)   
      ELSE  
         y(10)=zero 
      ENDIF
      y(11)=y(6)-y(9)-y(10)
      y(12)=y(14)-y(9)-y(10)
      y(13)=y(15)-y(9)-y(10)
C
C *** MAKE ALL ANGLES BETWEEN ZERO AND TWO PI
C
      y(6)=dmod(y(6)+two*tpi,tpi)
      DO 40 i=8,15
   40 y(i)=dmod(y(i)+two*tpi,tpi)
      RETURN
      END