ekepl1 Function

public pure function ekepl1(em, e)

Solve kepler's equation, em = ekepl - e*sin(ekepl), with legendre-based starter and halley iterator (function has also been used under the name eafkep)

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: em
real(kind=wp), intent(in) :: e

Return Value real(kind=wp)


Source Code

    pure function ekepl1(em, e)

    implicit none

    real(wp) :: ekepl1
    real(wp),intent(in) :: em
    real(wp),intent(in) :: e

    real(wp) :: c,s,psi,xi,eta,fd,fdd,f

    real(wp),parameter :: testsq = 1.0e-8_wp

    c = e*cos(em)
    s = e*sin(em)
    psi = s/sqrt(one - c - c + e*e)

    do

         xi = cos(psi)
         eta = sin(psi)
         fd = (one - c*xi) + s*eta
         fdd = c*eta + s*xi
         f = psi - fdd
         psi = psi - f*fd/(fd*fd - half*f*fdd)
         if (f*f < testsq) exit

     end do

     ekepl1 = em + psi

    end function ekepl1