Loading [MathJax]/jax/output/HTML-CSS/jax.js

invert_pe Subroutine

  1. submodule_geod_schw_implementation.f90
  2. geodesic_schwarzschild_implementation
  3. invert_pe

subroutine invert_pe(En, Lz, p, e)

Uses

  • proc~~invert_pe~~UsesGraph proc~invert_pe invert_pe module~parameters parameters proc~invert_pe->module~parameters module~kinds kinds module~parameters->module~kinds

Routine to convert from Energy and Angular momentum per unit mass to semi-latus rectum and eccentricity.

The equations in this routine was obtained using Mathematica by inverting the relationships E2=(p−2−2e)(p−2+2e)p(p−3−e2), L2z=p2M2p−3−e2. Note that this routine only works outside the separatrix and will give wrong results inside.

Arguments

Type IntentOptional AttributesName
real(kind=wp), intent(in) :: En

The energy per unit mass, E.
real(kind=wp), intent(in) :: Lz

The angular momentum per unit mass, Lz.
real(kind=wp), intent(out) :: p

On output the semi-latus rectum, p.
real(kind=wp), intent(out) :: e

On output the eccentricity.

Called by

proc~~invert_pe~~CalledByGraph proc~invert_pe invert_pe proc~geod_schw_save_globals_1 geod_schw_save_globals_1 proc~geod_schw_save_globals_1->proc~invert_pe interface~geod_schw_save_globals_1 geod_schw_save_globals_1 interface~geod_schw_save_globals_1->proc~geod_schw_save_globals_1

Contents

Source Code

invert_pe

Source Code

  subroutine invert_pe ( En, Lz, p, e )
  !! Routine to convert from Energy and Angular momentum per unit mass to
  !! semi-latus rectum and eccentricity.
  !!
  !! The equations in this routine was obtained using Mathematica by
  !! inverting the relationships
  !! \[
  !!    E^2 = \frac{(p-2-2e)(p-2+2e)}{p(p-3-e^2)},
  !! \]
  !! \[
  !!    L_z^2 = \frac{p^2 M^2}{p-3-e^2}.
  !! \]
  !! Note that this routine only works outside the separatrix and will give
  !! wrong results inside.
    use parameters, only : mass

    implicit none

    real(wp), intent(in) :: En
    !! The energy per unit mass, \(E\).
    real(wp), intent(in) :: Lz
    !! The angular momentum per unit mass, \(L_z\).
    real(wp), intent(out) :: p
    !! On output the semi-latus rectum, \(p\).
    real(wp), intent(out) :: e
    !! On output the eccentricity.
    complex(wp), parameter :: zi = (0.0_wp ,1.0_wp)
    complex(wp), parameter :: onepisqrt3 = 1.0_wp+zi*sqrt(3.0_wp)
    complex(wp), parameter :: onemisqrt3 = 1.0_wp-zi*sqrt(3.0_wp)
    real(wp), parameter :: threesqrt3 = 3.0_wp*sqrt(3.0_wp)
    real(wp), parameter :: third = 1.0_wp/3.0_wp
    complex(wp) :: pc, ec

    real(wp) :: E2, Lz2, M2, M4, M6, E22, Lz22, Lz23, onem3E2, E22Lz2, Lz2M2
    real(wp) :: onemE2, fac, twoLz2p4M2, sixE2M2inv, thirtysixE22M2inv
    complex(wp) :: sqrtfac, sqrtfac2, cuberootfac, cuberootfacinv, mainfac

    M2 = mass**2
    M4 = M2**2
    M6 = M2*M4

    E2 = En**2
    E22 = E2**2
    Lz2 = Lz**2
    Lz22 = Lz2**2
    Lz23 = Lz2*Lz22
    onem3E2 = 1.0_wp - 3.0_wp*E2
    E22Lz2 = E22*Lz2
    Lz2M2 = Lz2*M2
    onemE2 = 1.0_wp - E2

    sqrtfac = sqrt(cmplx(E22Lz2*M6*((8.0_wp-36.0*E2+27.0*E22)*Lz2M2 &
                                +16.0_wp*M4+Lz22*onemE2),0.0_wp,wp))
    sqrtfac2 = 8.0_wp*(-8.0_wp*M6+threesqrt3*sqrtfac)
    cuberootfac = (-Lz23-24.0_wp*(2.0_wp-6.0_wp*E2+9.0_wp*E22)*Lz2*M4 &
                   -12.0_wp*Lz22*M2*onem3E2+sqrtfac2)**third
    cuberootfacinv = 1.0_wp/cuberootfac
    fac = Lz22+16.0_wp*M4+8.0_wp*Lz2M2*onem3E2
    twoLz2p4M2 = 2.0_wp*(Lz2+4.0_wp*M2)
    mainfac = onemisqrt3*cuberootfac+onepisqrt3*cuberootfacinv*fac+twoLz2p4M2
    sixE2M2inv = 1.0_wp/(6.0_wp*E2*M2)
    thirtysixE22M2inv = sixE2M2inv/(6.0_wp*E2)

    pc = mainfac*sixe2M2inv
    ec = sqrt(-3.0_wp+mainfac*sixE2M2inv-mainfac**2*thirtysixE22M2inv/Lz2)

    p = real(pc,wp)
    e = real(ec,wp)
  end subroutine invert_pe