pv3els Subroutine

public pure subroutine pv3els(gm, pv, e)

Algorithm for three-dimensional conversion from position and velocity to orbital elements.

Arguments

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

grav. parameter [km^3/s^2]
real(kind=wp), intent(in), dimension(6) :: pv

[x, y, z, xdot, ydot, zdot]
real(kind=wp), intent(out), dimension(6) :: e

[al, q, ei, bom, om, tau]

Calls

proc~~pv3els~~CallsGraph proc~pv3els gooding_module::pv3els proc~pv2els gooding_module::pv2els proc~pv3els->proc~pv2els proc~emkep gooding_module::emkep proc~pv2els->proc~emkep proc~shmkep gooding_module::shmkep proc~pv2els->proc~shmkep

Called by

proc~~pv3els~~CalledByGraph proc~pv3els gooding_module::pv3els proc~lambert_test lambert_module::lambert_test proc~lambert_test->proc~pv3els proc~propagate gooding_module::propagate proc~propagate->proc~pv3els

Source Code

    pure subroutine pv3els (gm, pv, e)

    implicit none

    real(wp),intent(in)               :: gm   !! grav. parameter [km^3/s^2]
    real(wp),dimension(6),intent(in)  :: pv   !! [x, y, z, xdot, ydot, zdot]
    real(wp),dimension(6),intent(out) :: e    !! [al, q, ei, bom, om, tau]

    real(wp) :: x,y,z,xdot,ydot,zdot,al,q,ei,bom,om,tau,xsqysq,&
                rsq,r,vr,hx,hy,hz,hsq,u,vt,bx,by,bz,w,h

    if (all(pv==zero)) then

        e = zero

    else

        x    = pv(1)
        y    = pv(2)
        z    = pv(3)
        xdot = pv(4)
        ydot = pv(5)
        zdot = pv(6)

        xsqysq = x*x + y*y
        rsq = xsqysq + z*z
        r = sqrt(rsq)
        vr = (x*xdot + y*ydot + z*zdot)/r
        hx = y*zdot - z*ydot
        hy = z*xdot - x*zdot
        hz = x*ydot - y*xdot
        hsq = hx*hx + hy*hy + hz*hz
        if (hsq==zero) then
            !(rectilinear orbit)
            ei = halfpi
            if (xsqysq==zero) then
                !(axial orbit)
                bom = zero
            else
                !(general rectilinear orbit)
                bom = atan2(y, x)
            end if
            u = atan2(z, sqrt(xsqysq))
            vt = zero
        else
            !(non-degenerate orbit)
            bx = hy*z - hz*y
            by = hz*x - hx*z
            bz = hx*y - hy*x
            hx = y*bz - z*by
            hy = z*bx - x*bz
            hz = x*by - y*bx
            w = hx*hx + hy*hy
            h = sqrt(w + hz*hz)
            ei = atan2(sqrt(w), hz)
            if (w==zero) then
                !(orbit in reference plane)
                bom = zero
                u = atan2(y*sign(one,hz), x)
            else
                !(general orbit)
                bom = atan2(hx, -hy)
                u = atan2(h*z, rsq*bz)
            end if
            vt = h/(r*rsq)
        end if
        call pv2els (gm, r, u, vr, vt, al, q, om, tau)

        e(1) = al
        e(2) = q
        e(3) = ei
        e(4) = bom
        e(5) = om
        e(6) = tau

    end if

    end subroutine pv3els