periapsis_apoapsis Subroutine

    pure subroutine periapsis_apoapsis(mu,a,e,rp,ra,vp,va)

    implicit none

    real(wp),intent(in)  :: mu  !! gravitational parameter [\(km^{3}/s^{2}\)]
    real(wp),intent(in)  :: a   !! semimajor axis [km]
    real(wp),intent(in)  :: e   !! eccentricity [--]
    real(wp),intent(out) :: rp  !! periapsis position magnitude [km]
    real(wp),intent(out) :: ra  !! apoapsis position magnitude [km]
    real(wp),intent(out) :: vp  !! periapsis velocity magnitude [km/s]
    real(wp),intent(out) :: va  !! apoapsis velocity magnitude [km/s]

    real(wp) :: rarp   !! \(r_a + r_p \)
    real(wp) :: twomu  !! \( 2 \mu \)

    twomu = two * mu

    rp   = a*(one-e)
    ra   = a*(one+e)
    rarp = ra + rp
    vp   = sqrt(twomu*ra/(rp*rarp))
    va   = sqrt(twomu*rp/(ra*rarp))

    end subroutine periapsis_apoapsis