kepler.f

Go to the documentation of this file.

      SUBROUTINE kepler(AM,E,EA,SE,CE)  
C  VERSION OF 1/27/87
C  PURPOSE   
C    SOLVES KEPLER EQUATION
C  INPUT 
C    AM     = MEAN ANOMALY (RAD)
C    E      = ECCENTRICITY 
C  OUTPUT
C    EA     = ECCENTRIC ANOMALY (RAD)
C    SE     = SIN(EA)  
C    CE     = COS(EA)  
C  CALL SUBROUTINES  
C    NONE  
C  REFERENCES
C    JPL EM 312/87-153, 20 APRIL 1987
C    JPL IOM 312/79.4-426, FINAL REPORT JPL CONTRACT #955140   
C    AND USER'S GUIDE FOR MULCON, 1 JAN 79.
C  ANALYSIS  
C    J. H. KWOK - JPL  
C  PROGRAMMER
C    J. H. KWOK - JPL  
C  PROGRAM MODIFICATION  
C    NONE  
C  COMMENTS  
C    THIS PROGRAM IS A STRIPED VERSION OF THE ORGINAL VERSION FOR  
C    MULCON.  THIS VERSION DOES NOT WORK FOR HYPERBOLIC CASES.  
C    THE TOLERANCES ARE SET FOR DOUBLE PRECISION. CHANGE TOL1 AND  
C    TOL2 FOR OTHER PRECISION.   SEE REFERENCE 2 FOR ALGORITHM 
C    EXPLANATION.  TOL1 SHOULD BE SET TO 1/3 OF MACHINE PRECISION.
C    TOL2 SHOULD BE SET TO 1/2 OF MACHINE PRECISION.  FOR EXAMPLE,
C    TOL1=1.D-6 AND TOL2=1.D-9 FOR 18 DIGIT MACHINE PRECISION.
C
      IMPLICIT DOUBLE PRECISION (a-h,o-z)   
      DATA tol1,tol2/1.d-5,1.d-8/
      DATA half,zero,one/0.5d0,0.d0,1.d0/
      DATA pi,tpi/3.141592653589793d0,6.283185307179586d0/  
      k=0   
   10 CONTINUE  
      abm=dabs(am)  
      IF (abm.LE.pi) GO TO 20   
      IF (am.GT.pi) am=am-tpi   
      IF (am.LT.-pi) am=am+tpi  
      GO TO 10  
   20 sm=dsin(abm)  
      ea=abm+e*sm/(one-dsin(abm+e)+sm) 
   30 se=dsin(ea)   
      ce=dcos(ea)   
      k=k+1 
      IF (k.GT.10) stop 
   50 ese=e*se  
      f=abm+ese-ea  
      d=one-e*ce   
      c=f/(d+half*f*ese/d) 
      ea=ea+c   
      IF (dabs(c).GT.tol1) GO TO 30 
      IF (dabs(c).LT.tol2) GO TO 40 
      a=one-half*c*c  
      sen=a*se+c*ce 
      cen=a*ce-c*se 
      se=sen
      ce=cen
      GO TO 50  
   40 sen=se+c*ce   
      cen=ce-c*se   
      se=sen
      ce=cen
      IF (am.GT.zero) RETURN
      se=-se
      ea=-ea
      RETURN
      END