els2pv Subroutine

    pure subroutine els2pv (gm, al, q, om, tau, r, u, vr, vt)

    implicit none

    real(wp),intent(in)  :: gm    !! grav. parameter [km^3/s^2]
    real(wp),intent(in)  :: al    !! alpha [km^2/s^2]
    real(wp),intent(in)  :: q     !! periapsis distance [km]
    real(wp),intent(in)  :: om    !! argument of periapsis relative to assumed reference direction [rad]
    real(wp),intent(in)  :: tau   !! time from periapsis [sec]
    real(wp),intent(out) :: r     !! radial distance [km]
    real(wp),intent(out) :: u     !! angle from reference direction [rad]
    real(wp),intent(out) :: vr    !! radial velocity [km/2]
    real(wp),intent(out) :: vt    !! transverse velocity >=0 [km/s]

    real(wp) :: d,h,v,e1,e,ep1,alp,rtal,em,ee2,s2,c2,emv,s,c

    if (al==zero) then

        !(parabola - gm cannot be zero)

        d = dcbsol(half/gm, q, 1.5_wp*gm*tau)
        r = q + half*d*d/gm
        h = sqrt(two*gm*q)
        v = two*atan2(d,h)

    else

        !(ellipse or hyperbola)

        e1 = al*q
        e = gm - e1
        ep1 = gm + e
        h = sqrt(q*ep1)
        alp = abs(al)
        rtal = sqrt(alp)
        !(last 6 items could be saved if repeating gm, al & q)

        em = tau*alp*rtal
        if (al>zero) then

            !(ellipse - gm cannot be zero)
            ! make sure e1 argument to ekepl is between [0,1]
            ee2 = half*ekepl(em/gm, max(zero,min(one,e1/gm)))
            s2 = sin(ee2)
            c2 = cos(ee2)
            r = q + two*e*s2*s2/al
            d = two*e*s2*c2/rtal
            v = two*atan2(ep1*s2, h*rtal*c2)
            emv = em/gm - v
            v = v + fourpi*sign(real(int(abs(emv/fourpi) + half),wp), emv)

        else

            !(hyperbola)
            s = shkepl(em/e, -e1/e)
            s2 = s*s
            c = sqrt(one + s2)
            s2 = s2/(c+one)
            r = q - e*s2/al
            d = e*s/rtal
            v = atan2(s*h*rtal, -gm*s2 - e1)

        end if

    end if

    !(all orbits)
    u = om + v
    vr = d/r
    vt = h/r

    end subroutine els2pv