NASA Ocean Color

ocssw V2022
       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
#define f
Definition: l1_czcs.c:696
#define pi
Definition: vincenty.c:23
kepler
subroutine kepler(AM, E, EA, SE, CE)
Definition: kepler.f:2