calculate_bplane_data Subroutine

public subroutine calculate_bplane_data(mu, state, bdotr, bdott, bmag, theta, istat)

Compute B-plane parameters from position and velocity -- alternate version.

See also

Reference

  • This one is based on the algorithm in GMAT: GmatCalcUtil::CalculateBPlaneData

Arguments

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

central body grav parameter

real(kind=wp), intent(in), dimension(6) :: state

state vector (km,km/s)
real(kind=wp), intent(out) :: bdotr

(km)
real(kind=wp), intent(out) :: bdott

(km)
real(kind=wp), intent(out) :: bmag

magnitude of B vector (km)
real(kind=wp), intent(out) :: theta

aim point orientation [rad]
integer, intent(out) :: istat

status flag:

  • 0 if no errors.
  • -1 if state is not hyperbolic

Calls

proc~~calculate_bplane_data~~CallsGraph proc~calculate_bplane_data bplane_module::calculate_bplane_data proc~cross vector_module::cross proc~calculate_bplane_data->proc~cross proc~unit vector_module::unit proc~calculate_bplane_data->proc~unit

Called by

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

Source Code

    subroutine calculate_bplane_data(mu,state,bdotr,bdott,bmag,theta,istat)

    implicit none

    real(wp),intent(in)              :: mu     !! central body grav parameter \( (km^3/s^2) \)
    real(wp),dimension(6),intent(in) :: state  !! state vector (km,km/s)
    real(wp),intent(out)             :: bdotr  !! \( \mathbf{B} \cdot \mathbf{R} \) (km)
    real(wp),intent(out)             :: bdott  !! \( \mathbf{B} \cdot \mathbf{T} \) (km)
    real(wp),intent(out)             :: bmag   !! magnitude of B vector (km)
    real(wp),intent(out)             :: theta  !! aim point orientation [rad]
    integer,intent(out)              :: istat  !! status flag:
                                               !!
                                               !! * 0 if no errors.
                                               !! * -1 if state is not hyperbolic

    real(wp),dimension(3) :: r        !! position vector
    real(wp),dimension(3) :: v        !! velocity vector
    real(wp),dimension(3) :: evec     !! eccentricity vector
    real(wp),dimension(3) :: hvec     !! angular momentum vector
    real(wp),dimension(3) :: nvec     !! orbit normal vector
    real(wp),dimension(3) :: svec     !! incoming asymptote
    real(wp),dimension(3) :: bvec     !! B vector
    real(wp),dimension(3) :: tvec     !! T vector
    real(wp),dimension(3) :: rvec     !! R vector
    real(wp)              :: rmag     !! magnitude of `r`
    real(wp)              :: vmag     !! magnitude of `v`
    real(wp)              :: e        !! eccentricity
    real(wp)              :: hmag     !! magnitude of `hvec` vector
    real(wp)              :: b        !! semiminor axis
    real(wp)              :: oneovere !! 1/e
    real(wp)              :: temp

    r    = state(1:3)
    v    = state(4:6)
    hvec = cross(r, v)
    rmag = norm2(r)
    vmag = norm2(v)
    evec = cross(v,hvec)/mu - r/rmag
    e    = norm2(evec)

    if (e <= 1.0) then ! not hyperbolic
        write(*,*) 'error: state is not hyperbolic.'
        istat    = -1
        bdotr    = 0.0_wp
        bdott    = 0.0_wp
        bmag     = 0.0_wp
        theta    = 0.0_wp
    else
        istat    = 0
        evec     = unit(evec)
        hmag     = norm2(hvec)
        hvec     = unit(hvec)
        nvec     = cross(hvec, evec)
        b        = (hmag*hmag) / (mu * sqrt(e*e - 1.0_wp))
        oneovere = 1.0_wp/e
        temp     = sqrt(1.0_wp - oneovere*oneovere)
        svec     = (evec/e) + (temp*nvec)
        bvec     = b * (temp * evec - oneovere*nvec)
        tvec     = [svec(2), -svec(1), 0.0_wp] / sqrt(svec(1)*svec(1) + svec(2)*svec(2))
        rvec     = cross(svec, tvec)
        bdott    = dot_product(bvec, tvec)
        bdotr    = dot_product(bvec, rvec)
        bmag     = sqrt(bdott*bdott + bdotr*bdotr)
        theta    = atan(bdotr, bdott)
    end if

    end subroutine calculate_bplane_data