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 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, .
real(kind=wp), intent(in) :: Lz

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

On output the semi-latus rectum, .
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