ekepl Function

public pure function ekepl(em, e1)

Kepler's equation, em = ekepl - (1 - e1)*sin(ekepl), with e1 in range 1 to 0 inclusive, solved accurately (based on ekepl3, but entering e1, not e)

Arguments

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

Return Value real(kind=wp)


Calls

proc~~ekepl~~CallsGraph proc~ekepl gooding_module::ekepl proc~dcbsol gooding_module::dcbsol proc~ekepl->proc~dcbsol proc~emkep gooding_module::emkep proc~ekepl->proc~emkep proc~dcubrt gooding_module::dcubrt proc~dcbsol->proc~dcubrt

Called by

proc~~ekepl~~CalledByGraph proc~ekepl gooding_module::ekepl proc~els2pv gooding_module::els2pv proc~els2pv->proc~ekepl proc~els3pv gooding_module::els3pv proc~els3pv->proc~els2pv proc~propagate gooding_module::propagate proc~propagate->proc~els3pv

Source Code

    pure function ekepl(em, e1)

    implicit none

    real(wp) :: ekepl
    real(wp),intent(in) :: em
    real(wp),intent(in) :: e1

    real(wp) :: emr,ee,e,w,fdd,fddd,f,fd,dee
    integer :: iter

    real(wp),parameter :: sw = 0.25_wp

    !range-reduce em to lie in range -pi to pi
    emr = mod(em,twopi)
    if (emr<pineg) emr = emr + twopi
    if (emr>pi) emr = emr - twopi
    ee = emr

    if (ee/=zero) then

        if (ee<zero) ee = -ee

        !(emr is range-reduced em & ee is absolute value of emr)
        !starter by first solving cubic equation
        e = one - e1
         w = dcbsol(e,two*e1, three*ee)

        !effectively interpolate in emr (absolute value)
        ee = (ee*ee + (pi - ee)*w)/pi
        if (emr<zero) ee = -ee

        !do two iterations of halley, each followed by newton
        do iter=1,2
            fdd = e*sin(ee)
            fddd = e*cos(ee)
            if (ee*ee/six + e1 >= sw) then
                f = (ee - fdd) - emr
                fd = one - fddd
            else
                f = emkep(e1,ee) - emr
                fd = two*e*sin(half*ee)**2 + e1
            end if
            dee = f*fd/(half*f*fdd - fd*fd)
            f = f + dee*(fd + half*dee*(fdd + athird*dee*fddd))
            !to reduce the danger of underflow replace the last line by
            !    w = fd + half*dee*(fdd + athird*dee*fddd)
            fd = fd + dee*(fdd + half*dee*fddd)
            ee = ee + dee - f/fd
            !if replacing as above, then also replace the last line by
            !ee = ee - (f - dee*(fd - w))/fd
        end do

    end if

    !range-expand
    ekepl = ee + (em - emr)

    end function ekepl