shkepl Function

    pure function shkepl (el, g1)

    implicit none

    real(wp)            :: shkepl
    real(wp),intent(in) :: el
    real(wp),intent(in) :: g1

    real(wp) :: s,g,cl,al,w,s0,s1,s2,s3,fdd,fddd,f,fd,ds,stemp
    integer  :: iter

    real(wp),parameter :: sw=half

    s = el

    if (el/=zero) then

        !started based on lagrange's theorem
        g = one - g1
        cl = sqrt(one + el**2)
        al = asinh(el)
        w = g**2*al/cl**3
        s = one - g/cl
        s = el + g*al/dcubrt(s**3 + w*el*(1.5_wp - g/0.75_wp))

        !two iterations (at most) of halley-then-newton process
        do iter=1,2
            s0 = s*s
            s1 = s0 + one
            s2 = sqrt(s1)
            s3 = s1*s2
            fdd = g*s/s3
            fddd = g*(one - two*s0)/(s1*s3)
            if (asixth*s0 + g1 >= sw) then
                f = (s - g*asinh(s)) - el
                fd = one - g/s2
            else
                f = shmkep(g1, s) - el
                fd = (s0/(s2 + one) + g1)/s2
            end if
            ds = f*fd/(half*f*fdd - fd*fd)
            stemp = s + ds
            if (stemp==s) exit
            f = f + ds*(fd + half*ds*(fdd + athird*ds*fddd))
            fd = fd + ds*(fdd + half*ds*fddd)
             s = stemp - f/fd
        end do

    end if

    shkepl = s

     end function shkepl