Root of Kepler's Equation

The equation of eccentric anomaly is given by

$$\begin{array}{}\text{(1)}& \psi -e\sin \psi =\omega t\end{array}$$
Using simple iteration which can be written as

$$\begin{array}{}\text{(2)}& {\psi }_{n+1}=\omega t+e\sin {\psi }_n\end{array}$$
, in which I iterate over and over until the difference between

${\psi }_{n+1}$ and

${\psi }_n$ is less than

$\epsilon$. In this case I let

$\epsilon =1.D0$.

Code

!========================================================
!   Purpose: Solve Kepler eccentric anomaly
!
!   Methods: Simple iteration (E=M+e*sinE)
!       
!       Date        Programer       Description of change
!       ====        =========       =====================
!      9/30/20        Morris        Original Code
!========================================================
PROGRAM kepler
    IMPLICIT NONE
    REAL(16), PARAMETER :: pi = 3.1415926535898D0
    REAL(16) :: phi, e, M, tmp
    INTEGER :: i

    open(1, file='kepler_2.dat', status='unknown')
    e = 0.2D0
    do i = 1, 100
        M = pi * dble(i) / 50.D0
        do 
            tmp = phi
            phi = M + e * sin(phi) 
            write(*, *) abs(tmp-phi)
            if (abs(tmp - phi) <= 1D-8) then
                EXIT
            end if
        end do
        write(1, '(3F20.17)') M, phi
    end do
END PROGRAM kepler

Result

In this program, it tries to find

$\psi$ as a function of

$\omega t$. By letting

$\omega t$ go through

$0$ to

$2\pi$. I can get the following relation. Futhermore, I can adjust different

$e$ to get a different function.