compute_vinf_vectors Subroutine

private subroutine compute_vinf_vectors(mu, rv, vinfin, vinfout)

Compute the incoming and/or outgoing v-infinity vectors, given the position and velocity of a hyperbola.

Note

This is for testing the other routines.

Arguments

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

central body gravitational parameter
real(kind=wp), intent(in), dimension(6) :: rv

position,velocity vector
real(kind=wp), intent(out), optional, dimension(3) :: vinfin

incoming v-infinity vector
real(kind=wp), intent(out), optional, dimension(3) :: vinfout

outgoing v-infinity vector

Calls

proc~~compute_vinf_vectors~~CallsGraph proc~compute_vinf_vectors bplane_module::compute_vinf_vectors proc~cross vector_module::cross proc~compute_vinf_vectors->proc~cross

Called by

proc~~compute_vinf_vectors~~CalledByGraph proc~compute_vinf_vectors bplane_module::compute_vinf_vectors proc~bplane_test bplane_module::bplane_test proc~bplane_test->proc~compute_vinf_vectors

Source Code

    subroutine compute_vinf_vectors(mu,rv,vinfin,vinfout)

    implicit none

    real(wp),intent(in)                        :: mu      !! central body gravitational parameter
    real(wp),dimension(6),intent(in)           :: rv      !! position,velocity vector
    real(wp),dimension(3),intent(out),optional :: vinfin  !! incoming v-infinity vector
    real(wp),dimension(3),intent(out),optional :: vinfout !! outgoing v-infinity vector

    real(wp),dimension(3) :: h,e,p,q,r,v
    real(wp) :: rmag,vmag,vinfmag,emag,qmag,cbeta,sbeta

    if (present(vinfin) .or. present(vinfout)) then
        r       = rv(1:3)
        v       = rv(4:6)
        rmag    = norm2(r)
        vmag    = norm2(v)
        h       = cross(r,v)
        q       = cross(v,h)
        vinfmag = sqrt(vmag*vmag-two*mu/rmag)
        e       = q/mu-r/rmag
        emag    = norm2(e)
        q       = cross(h,e)
        qmag    = norm2(q)
        cbeta   = one/emag
        sbeta   = sqrt(one-cbeta*cbeta)
        p       = e/emag
        q       = q/qmag
        if (present(vinfin))  vinfin  = vinfmag*( cbeta*p+sbeta*q)
        if (present(vinfout)) vinfout = vinfmag*(-cbeta*p+sbeta*q)
    end if

    end subroutine compute_vinf_vectors