ekepl2 Function

    pure function ekepl2(em, e)

    implicit none

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

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

    real(wp),parameter :: sw = 0.1_wp
    real(wp),parameter :: a  = (pi-one)**2/(pi+two/three)
    real(wp),parameter :: b  = two*(pi-asixth)**2/(pi+two/three)

    !range-reduce em to line 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)
        !started for e = 1 by cube root of bilinear function
        if (ee<asixth) then
            ee = (six*ee)**athird
        else
            w = pi - ee
            ee = pi - a*w/(b - w)
        end if
        if (emr<zero) ee = -ee

        !interpolate for e
        ee = emr + (ee - emr)*e

        !do two iterations of halley, each followed by newton
        e1 = one - e
        l = (e1 + ee*ee/six) >= sw
        do iter=1,2
            fdd = e*sin(ee)
            fddd = e*cos(ee)
            if (l) then
                f = (ee - fdd) - emr
                fd = one - fddd
            else
                f = emkepl(e,ee) - emr
                fd = e1 + two*e*sin(half*ee)**2
            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
    ekepl2 = ee + (em - emr)

    end function ekepl2