      program kepler
C
C     E_TRUE      ECCENTRICITY
C     M_ANOM      MEAN ANOMALY
C     E_ANOM      ECCENTRIC ANOMALY
C     T_ANOM      TRUE ANOMALY
C
      real*8 e_true,m_anom,e_anom,t_anom,twopi
      twopi=2.0d0*3.141592653589793d0

      write(6,*)"Enter e, M!"
      read(5,*)e_true,m_anom

      if(abs(m_anom).le.10.d-14)then
       m_anom=0.D0
       e_anom=0.D0
       t_anom=0.D0
      elseif(abs(abs(m_anom)-pi).le.10.d-14)then
       e_anom=m_anom
       t_anom=e_anom
      else

C GET THE FIRST ESTIMATE OF THE ECCENTRIC ANOMALY

      E_ANOM      = M_ANOM + E_TRUE * SIN(M_ANOM)
     1                   + 0.5D0*E_TRUE**2 * SIN(2.D0*M_ANOM)

C ITERATE TO GET THE ECCENTRIC ANOMALY

      DO J = 1, 6
          E_ANOM  = E_ANOM + (M_ANOM-E_ANOM+E_TRUE*SIN(E_ANOM))/
     1                         (1.D0 - E_TRUE * COS(E_ANOM))
      END DO

C GET THE TRUE ANOMALY

      T_ANOM      = 2.D0*ATAN(
     1                 SQRT((1.D0+E_TRUE)/(1.D0-E_TRUE))
     2                 * TAN(E_ANOM/2.D0)
     3                )
      endif

C SINCE E_ANOM IS ALWAYS POSITIVE, T_ANOM SHOULD BE TOO:

      IF (T_ANOM.LT.0.) THEN
          T_ANOM = TWOPI+DMOD(T_ANOM,TWOPI)
      ELSE
          T_ANOM = DMOD(T_ANOM,TWOPI)
      END IF

      write(6,*)t_anom/twopi*360

      END

