pv2els Subroutine

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

    implicit none

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

    real(wp) :: esq1,es,eses,ec,ecec,esq,e,v,e1
    real(wp) :: eh,em,ecesq,en,adj,vsq,rtal,d,h,p,alp

    real(wp),parameter   :: sw = 0.25_wp
    logical,parameter    :: l = .false.

    !(all orbits)
    vsq = vr*vr + vt*vt
    al = two*gm/r - vsq
    alp = abs(al)
    rtal = sqrt(alp)
    d = r*vr
    h = r*vt
    p = h*h
    esq1 = p*al
    es = d*rtal
    eses = es*es
    ec = r*vsq - gm
    ecec = ec*ec
    if (al>zero) then
        !(one esq formula superior for the ellipse)
        esq = ecec + eses
    else
        !(different formula superior for the hyperbola)
        esq = gm*gm - esq1
    end if
    e = sqrt(esq)
    q = p/(gm + e)
    if (al==zero) then
        !(parabola)
        tau = d*(two*q + r)/(three*gm)
        v = two*atan2(vr, vt)
    else if (e==zero) then
        !(circle)
        tau = zero
        v = zero
    else
        !(ellipse or hyperbola)
        e1 = al*q
        if (al>zero) then
            !(ellipse)
            eh = atan2(es, ec)
            if (gm*eh*eh/six + e1 >= gm*sw) then
                !(general case)
                em = gm*eh - es
                ecesq = gm*ec - esq
            else
                !(for e1 & eh both near zero)
                em = gm*emkep (e1/gm, eh)
                ecesq = (esq1*ecec - esq*eses)/(esq + gm*ec)
            end if
        else
            !(hyperbola)
            eh = asinh(es/e)
            if (gm*eh*eh/six - e1 >= gm*sw) then
                !(general case)
                em = es - gm*eh
                ecesq = esq - gm*ec
            else
                !(for e1 & eh both near zero)
                em = e*shmkep(-e1/e, es/e)
                ecesq = -(esq1*ecec + esq*eses)/(esq + gm*ec)
            end if
        end if
            !(ellipse or hyperbola still)
            en = alp*rtal
            tau = em/en
            v = atan2(es*h*rtal, ecesq)
    end if

    !(all orbits)
    om = u - v

    !
    !  note: the following is never executed... set l=true and test...
    !

    if (l .and. al>zero) then
        !(for ellipse, adjust revolutions if required (using l))
        adj = twopi*sign(real(int(abs(om/twopi) + half),wp), om)
        om = om - adj
        tau = tau + adj/en
    end if

    end subroutine pv2els