subroutine bplane(mu,rv,vinfvec,bmag,theta,BdotT,BdotR,status_ok)

    implicit none

    real(wp),intent(in)               :: mu        !! central body grav parameter \( (km^3/s^2) \)
    real(wp),dimension(6),intent(in)  :: rv        !! state vector (km,km/s)
    real(wp),dimension(3),intent(out) :: vinfvec   !! incoming V-infinity vector (km/s)
    real(wp),intent(out)              :: bmag      !! magnitude of B vector (km)
    real(wp),intent(out)              :: theta     !! aim point orientation [rad]
    real(wp),intent(out)              :: BdotT     !! \( \mathbf{B} \cdot \mathbf{T} \) (km)
    real(wp),intent(out)              :: BdotR     !! \( \mathbf{B} \cdot \mathbf{R} \) (km)
    logical,intent(out)               :: status_ok !! false if there were errors (non-hyperbolic or degenerate state)

    !local variables:
    real(wp),dimension(3) :: r,v,h,evec,bvec,tvec,ehat,hhat,hehat,Shat,That,Rhat,Bhat
    real(wp) :: rmag,vmag,vinf2,rdv,vmag2,a,alpha,sd2,cd2,ct,st,vinf,e

    status_ok = .true.

    r         = rv(1:3)
    v         = rv(4:6)
    h         = cross(r,v)
    hhat      = unit(h)

    if (all(hhat==zero)) then
        write(*,*) 'error: degenerate state.'
        status_ok = .false.
        return
    end if

    rdv       = dot_product(r,v)
    rmag      = norm2(r)
    vmag      = norm2(v)
    vmag2     = vmag*vmag
    vinf2     = vmag2 - two*mu/rmag
    vinf      = sqrt(vinf2)                     ! magnitude of v-infinity vector
    evec      = cross(v,h)/mu - r/rmag          ! eccentricity vector
    e         = norm2(evec)                     ! eccentricity

    if (e<=one) then
        write(*,*) 'error: state is not hyperbolic.'
        status_ok = .false.
        return
    end if

    ehat      = evec/e                          ! eccentricity unit vector
    a         = one / (two/rmag - vmag2/mu)     ! semi-major axis
    hehat     = cross(hhat,ehat)                ! h x e unit vector
    sd2       = one/e                           ! sin(delta/2)
    cd2       = sqrt(one-sd2*sd2)               ! cos(delta/2)
    Shat      = cd2*hehat + sd2*ehat            ! incoming vinf unit vector
    !Shat     = cd2*hehat - sd2*ehat            ! outgoing vinf unit vector
    That      = ucross(Shat,[zero, zero, one])  ! here we define Tvec relative to the Z-axis of
                                                ! the frame in which the state is defined

    if (all(That==zero)) then
        write(*,*) 'error: vinf vector is parallel to z-axis.'
        status_ok = .false.
        return
    end if

    Rhat      = cross(Shat,That)             ! Eqn 1 in [1]
    Bhat      = ucross(Shat,h)
    ct        = dot_product(Bhat,That)       ! cos(theta)
    st        = dot_product(Bhat,Rhat)       ! sin(theta)

    !outputs:
    Bmag      = abs(a)*sqrt( e*e - one )     ! magnitude of B vector
    theta     = atan2(st,ct)                 ! aim point orientation
    vinfvec   = vinf * Shat                  ! incoming vinf vector
    BdotT     = bmag*ct                      ! B dot T
    BdotR     = bmag*st                      ! B dot R

    end subroutine bplane