From a9c7b94d7bcef9532c4da15ad363f8c73435bb1c Mon Sep 17 00:00:00 2001 From: Se Yong Park Date: Thu, 5 Oct 2023 23:14:12 -0400 Subject: [PATCH] refactored code and updated readme --- README.md | 4 + ...rrgeopy.constants.apex_from_constants.rst} | 2 +- ...rgeopy.constants.plunging_radial_roots.rst | 6 + .../_autosummary/kerrgeopy.constants.rst | 4 + ...kerrgeopy.constants.stable_polar_roots.rst | 6 + ...errgeopy.constants.stable_radial_roots.rst | 6 + ...frequencies.plunging_mino_frequencies.rst} | 2 +- .../_autosummary/kerrgeopy.frequencies.rst | 1 + .../kerrgeopy.frequencies_from_constants.rst | 23 - .../kerrgeopy.initial_conditions.rst | 1 - .../_autosummary/kerrgeopy.orbit.Orbit.rst | 3 +- ...rst => kerrgeopy.plunge.PlungingOrbit.rst} | 5 +- ...geopy.plunge.plunging_polar_solutions.rst} | 2 +- ...eopy.plunge.plunging_radial_integrals.rst} | 2 +- ...nge.plunging_radial_solutions_complex.rst} | 2 +- .../kerrgeopy.plunge.plunging_trajectory.rst | 6 + ...ing_solutions.rst => kerrgeopy.plunge.rst} | 16 +- .../_autosummary/kerrgeopy.plunging_orbit.rst | 31 -- docs/source/_autosummary/kerrgeopy.rst | 7 +- .../kerrgeopy.spacetime.KerrSpacetime.rst | 2 - ...t.rst => kerrgeopy.stable.StableOrbit.rst} | 5 +- ...t => kerrgeopy.stable.polar_solutions.rst} | 2 +- ... => kerrgeopy.stable.radial_solutions.rst} | 2 +- ...ble_solutions.rst => kerrgeopy.stable.rst} | 15 +- .../kerrgeopy.stable.stable_trajectory.rst | 6 + .../_autosummary/kerrgeopy.stable_orbit.rst | 31 -- docs/source/index.rst | 8 +- docs/source/notebooks/Getting Started.ipynb | 38 +- docs/source/notebooks/Graphics.ipynb | 6 +- .../source/notebooks/Orbital Properties.ipynb | 51 ++- docs/source/notebooks/Trajectory.ipynb | 28 +- kerrgeopy/__init__.py | 8 +- kerrgeopy/constants.py | 145 +++++- kerrgeopy/frequencies.py | 365 ++++++++++++--- kerrgeopy/frequencies_from_constants.py | 231 ---------- kerrgeopy/initial_conditions.py | 47 +- kerrgeopy/orbit.py | 24 +- kerrgeopy/plunge.py | 411 +++++++++++++++++ kerrgeopy/plunging_orbit.py | 160 ------- kerrgeopy/plunging_solutions.py | 292 ------------ kerrgeopy/spacetime.py | 34 +- kerrgeopy/stable.py | 319 ++++++++++++++ kerrgeopy/stable_orbit.py | 182 -------- kerrgeopy/stable_solutions.py | 144 ------ setup.py | 4 +- tests/{data => }/README.txt | 21 +- tests/data/Unit Test Generation.ipynb | 417 ++++++++++++++++++ ...{const_values.txt => const_parameters.txt} | 0 .../{freq_values.txt => freq_parameters.txt} | 0 ...x_values.txt => separatrix_parameters.txt} | 0 ...values.txt => stable_orbit_parameters.txt} | 0 tests/test_constants.py | 4 +- tests/test_four_velocity.py | 6 +- tests/test_frequencies.py | 17 +- tests/test_initial_conditions.py | 2 +- tests/test_plunging_solutions.py | 2 +- tests/test_stable_solutions.py | 10 +- 57 files changed, 1796 insertions(+), 1372 deletions(-) rename docs/source/_autosummary/{kerrgeopy.initial_conditions.apex_from_constants.rst => kerrgeopy.constants.apex_from_constants.rst} (63%) create mode 100644 docs/source/_autosummary/kerrgeopy.constants.plunging_radial_roots.rst create mode 100644 docs/source/_autosummary/kerrgeopy.constants.stable_polar_roots.rst create mode 100644 docs/source/_autosummary/kerrgeopy.constants.stable_radial_roots.rst rename docs/source/_autosummary/{kerrgeopy.plunging_solutions.plunging_mino_frequencies.rst => kerrgeopy.frequencies.plunging_mino_frequencies.rst} (67%) delete mode 100644 docs/source/_autosummary/kerrgeopy.frequencies_from_constants.rst rename docs/source/_autosummary/{kerrgeopy.plunging_orbit.PlungingOrbit.rst => kerrgeopy.plunge.PlungingOrbit.rst} (78%) rename docs/source/_autosummary/{kerrgeopy.plunging_solutions.plunging_polar_solutions.rst => kerrgeopy.plunge.plunging_polar_solutions.rst} (67%) rename docs/source/_autosummary/{kerrgeopy.plunging_solutions.plunging_radial_integrals.rst => kerrgeopy.plunge.plunging_radial_integrals.rst} (67%) rename docs/source/_autosummary/{kerrgeopy.plunging_solutions.plunging_radial_solutions_complex.rst => kerrgeopy.plunge.plunging_radial_solutions_complex.rst} (72%) create mode 100644 docs/source/_autosummary/kerrgeopy.plunge.plunging_trajectory.rst rename docs/source/_autosummary/{kerrgeopy.plunging_solutions.rst => kerrgeopy.plunge.rst} (57%) delete mode 100644 docs/source/_autosummary/kerrgeopy.plunging_orbit.rst rename docs/source/_autosummary/{kerrgeopy.stable_orbit.StableOrbit.rst => kerrgeopy.stable.StableOrbit.rst} (82%) rename docs/source/_autosummary/{kerrgeopy.stable_solutions.polar_solutions.rst => kerrgeopy.stable.polar_solutions.rst} (60%) rename docs/source/_autosummary/{kerrgeopy.stable_solutions.radial_solutions.rst => kerrgeopy.stable.radial_solutions.rst} (61%) rename docs/source/_autosummary/{kerrgeopy.stable_solutions.rst => kerrgeopy.stable.rst} (51%) create mode 100644 docs/source/_autosummary/kerrgeopy.stable.stable_trajectory.rst delete mode 100644 docs/source/_autosummary/kerrgeopy.stable_orbit.rst delete mode 100644 kerrgeopy/frequencies_from_constants.py create mode 100644 kerrgeopy/plunge.py delete mode 100644 kerrgeopy/plunging_orbit.py delete mode 100644 kerrgeopy/plunging_solutions.py create mode 100644 kerrgeopy/stable.py delete mode 100644 kerrgeopy/stable_orbit.py delete mode 100644 kerrgeopy/stable_solutions.py rename tests/{data => }/README.txt (68%) create mode 100644 tests/data/Unit Test Generation.ipynb rename tests/data/{const_values.txt => const_parameters.txt} (100%) rename tests/data/{freq_values.txt => freq_parameters.txt} (100%) rename tests/data/{separatrix_values.txt => separatrix_parameters.txt} (100%) rename tests/data/{stable_orbit_values.txt => stable_orbit_parameters.txt} (100%) diff --git a/README.md b/README.md index d22063c..8b0af1a 100644 --- a/README.md +++ b/README.md @@ -272,3 +272,7 @@ plt.ylabel(r"$\phi(\lambda)$") ![png](README_files/Getting%20Started_20_1.png) +## Credit + +* Seyong Park +* Zach Nasipak \ No newline at end of file diff --git a/docs/source/_autosummary/kerrgeopy.initial_conditions.apex_from_constants.rst b/docs/source/_autosummary/kerrgeopy.constants.apex_from_constants.rst similarity index 63% rename from docs/source/_autosummary/kerrgeopy.initial_conditions.apex_from_constants.rst rename to docs/source/_autosummary/kerrgeopy.constants.apex_from_constants.rst index 7063a03..37b8f20 100644 --- a/docs/source/_autosummary/kerrgeopy.initial_conditions.apex_from_constants.rst +++ b/docs/source/_autosummary/kerrgeopy.constants.apex_from_constants.rst @@ -1,6 +1,6 @@ apex\_from\_constants ===================== -.. currentmodule:: kerrgeopy.initial_conditions +.. currentmodule:: kerrgeopy.constants .. autofunction:: apex_from_constants \ No newline at end of file diff --git a/docs/source/_autosummary/kerrgeopy.constants.plunging_radial_roots.rst b/docs/source/_autosummary/kerrgeopy.constants.plunging_radial_roots.rst new file mode 100644 index 0000000..331c1c3 --- /dev/null +++ b/docs/source/_autosummary/kerrgeopy.constants.plunging_radial_roots.rst @@ -0,0 +1,6 @@ +plunging\_radial\_roots +======================= + +.. currentmodule:: kerrgeopy.constants + +.. autofunction:: plunging_radial_roots \ No newline at end of file diff --git a/docs/source/_autosummary/kerrgeopy.constants.rst b/docs/source/_autosummary/kerrgeopy.constants.rst index 9faa5c9..31ee9ea 100644 --- a/docs/source/_autosummary/kerrgeopy.constants.rst +++ b/docs/source/_autosummary/kerrgeopy.constants.rst @@ -16,13 +16,17 @@ constants :template: custom-base-template.rst angular_momentum + apex_from_constants carter_constant constants_of_motion energy fast_separatrix is_stable + plunging_radial_roots scale_constants separatrix + stable_polar_roots + stable_radial_roots valid_params diff --git a/docs/source/_autosummary/kerrgeopy.constants.stable_polar_roots.rst b/docs/source/_autosummary/kerrgeopy.constants.stable_polar_roots.rst new file mode 100644 index 0000000..4e153a5 --- /dev/null +++ b/docs/source/_autosummary/kerrgeopy.constants.stable_polar_roots.rst @@ -0,0 +1,6 @@ +stable\_polar\_roots +==================== + +.. currentmodule:: kerrgeopy.constants + +.. autofunction:: stable_polar_roots \ No newline at end of file diff --git a/docs/source/_autosummary/kerrgeopy.constants.stable_radial_roots.rst b/docs/source/_autosummary/kerrgeopy.constants.stable_radial_roots.rst new file mode 100644 index 0000000..9ad7d4a --- /dev/null +++ b/docs/source/_autosummary/kerrgeopy.constants.stable_radial_roots.rst @@ -0,0 +1,6 @@ +stable\_radial\_roots +===================== + +.. currentmodule:: kerrgeopy.constants + +.. autofunction:: stable_radial_roots \ No newline at end of file diff --git a/docs/source/_autosummary/kerrgeopy.plunging_solutions.plunging_mino_frequencies.rst b/docs/source/_autosummary/kerrgeopy.frequencies.plunging_mino_frequencies.rst similarity index 67% rename from docs/source/_autosummary/kerrgeopy.plunging_solutions.plunging_mino_frequencies.rst rename to docs/source/_autosummary/kerrgeopy.frequencies.plunging_mino_frequencies.rst index 1d42fe1..7f9e7ff 100644 --- a/docs/source/_autosummary/kerrgeopy.plunging_solutions.plunging_mino_frequencies.rst +++ b/docs/source/_autosummary/kerrgeopy.frequencies.plunging_mino_frequencies.rst @@ -1,6 +1,6 @@ plunging\_mino\_frequencies =========================== -.. currentmodule:: kerrgeopy.plunging_solutions +.. currentmodule:: kerrgeopy.frequencies .. autofunction:: plunging_mino_frequencies \ No newline at end of file diff --git a/docs/source/_autosummary/kerrgeopy.frequencies.rst b/docs/source/_autosummary/kerrgeopy.frequencies.rst index 1041be7..7ecd9f9 100644 --- a/docs/source/_autosummary/kerrgeopy.frequencies.rst +++ b/docs/source/_autosummary/kerrgeopy.frequencies.rst @@ -19,6 +19,7 @@ frequencies gamma mino_frequencies phi_frequency + plunging_mino_frequencies r_frequency theta_frequency diff --git a/docs/source/_autosummary/kerrgeopy.frequencies_from_constants.rst b/docs/source/_autosummary/kerrgeopy.frequencies_from_constants.rst deleted file mode 100644 index d7c3480..0000000 --- a/docs/source/_autosummary/kerrgeopy.frequencies_from_constants.rst +++ /dev/null @@ -1,23 +0,0 @@ -frequencies\_from\_constants -============================ - -.. automodule:: kerrgeopy.frequencies_from_constants - - - - - - - - - - - - - - - - - - - diff --git a/docs/source/_autosummary/kerrgeopy.initial_conditions.rst b/docs/source/_autosummary/kerrgeopy.initial_conditions.rst index 06e8df5..7daaa2e 100644 --- a/docs/source/_autosummary/kerrgeopy.initial_conditions.rst +++ b/docs/source/_autosummary/kerrgeopy.initial_conditions.rst @@ -15,7 +15,6 @@ initial\_conditions :toctree: :template: custom-base-template.rst - apex_from_constants constants_from_initial_conditions is_stable plunging_orbit_initial_phases diff --git a/docs/source/_autosummary/kerrgeopy.orbit.Orbit.rst b/docs/source/_autosummary/kerrgeopy.orbit.Orbit.rst index f90be49..c4d9e50 100644 --- a/docs/source/_autosummary/kerrgeopy.orbit.Orbit.rst +++ b/docs/source/_autosummary/kerrgeopy.orbit.Orbit.rst @@ -20,9 +20,8 @@ ~Orbit.animate ~Orbit.constants_of_motion ~Orbit.four_velocity - ~Orbit.four_velocity_norm ~Orbit.is_visible - ~Orbit.orbital_parameters + ~Orbit.numerical_four_velocity ~Orbit.plot ~Orbit.trajectory diff --git a/docs/source/_autosummary/kerrgeopy.plunging_orbit.PlungingOrbit.rst b/docs/source/_autosummary/kerrgeopy.plunge.PlungingOrbit.rst similarity index 78% rename from docs/source/_autosummary/kerrgeopy.plunging_orbit.PlungingOrbit.rst rename to docs/source/_autosummary/kerrgeopy.plunge.PlungingOrbit.rst index 72b3497..0663e8a 100644 --- a/docs/source/_autosummary/kerrgeopy.plunging_orbit.PlungingOrbit.rst +++ b/docs/source/_autosummary/kerrgeopy.plunge.PlungingOrbit.rst @@ -1,7 +1,7 @@ PlungingOrbit ============= -.. currentmodule:: kerrgeopy.plunging_orbit +.. currentmodule:: kerrgeopy.plunge .. autoclass:: PlungingOrbit :members: @@ -20,9 +20,8 @@ ~PlungingOrbit.animate ~PlungingOrbit.constants_of_motion ~PlungingOrbit.four_velocity - ~PlungingOrbit.four_velocity_norm ~PlungingOrbit.is_visible - ~PlungingOrbit.orbital_parameters + ~PlungingOrbit.numerical_four_velocity ~PlungingOrbit.plot ~PlungingOrbit.trajectory diff --git a/docs/source/_autosummary/kerrgeopy.plunging_solutions.plunging_polar_solutions.rst b/docs/source/_autosummary/kerrgeopy.plunge.plunging_polar_solutions.rst similarity index 67% rename from docs/source/_autosummary/kerrgeopy.plunging_solutions.plunging_polar_solutions.rst rename to docs/source/_autosummary/kerrgeopy.plunge.plunging_polar_solutions.rst index f44784a..89ade07 100644 --- a/docs/source/_autosummary/kerrgeopy.plunging_solutions.plunging_polar_solutions.rst +++ b/docs/source/_autosummary/kerrgeopy.plunge.plunging_polar_solutions.rst @@ -1,6 +1,6 @@ plunging\_polar\_solutions ========================== -.. currentmodule:: kerrgeopy.plunging_solutions +.. currentmodule:: kerrgeopy.plunge .. autofunction:: plunging_polar_solutions \ No newline at end of file diff --git a/docs/source/_autosummary/kerrgeopy.plunging_solutions.plunging_radial_integrals.rst b/docs/source/_autosummary/kerrgeopy.plunge.plunging_radial_integrals.rst similarity index 67% rename from docs/source/_autosummary/kerrgeopy.plunging_solutions.plunging_radial_integrals.rst rename to docs/source/_autosummary/kerrgeopy.plunge.plunging_radial_integrals.rst index b2362e3..a5916c7 100644 --- a/docs/source/_autosummary/kerrgeopy.plunging_solutions.plunging_radial_integrals.rst +++ b/docs/source/_autosummary/kerrgeopy.plunge.plunging_radial_integrals.rst @@ -1,6 +1,6 @@ plunging\_radial\_integrals =========================== -.. currentmodule:: kerrgeopy.plunging_solutions +.. currentmodule:: kerrgeopy.plunge .. autofunction:: plunging_radial_integrals \ No newline at end of file diff --git a/docs/source/_autosummary/kerrgeopy.plunging_solutions.plunging_radial_solutions_complex.rst b/docs/source/_autosummary/kerrgeopy.plunge.plunging_radial_solutions_complex.rst similarity index 72% rename from docs/source/_autosummary/kerrgeopy.plunging_solutions.plunging_radial_solutions_complex.rst rename to docs/source/_autosummary/kerrgeopy.plunge.plunging_radial_solutions_complex.rst index e83c295..ca6ed3d 100644 --- a/docs/source/_autosummary/kerrgeopy.plunging_solutions.plunging_radial_solutions_complex.rst +++ b/docs/source/_autosummary/kerrgeopy.plunge.plunging_radial_solutions_complex.rst @@ -1,6 +1,6 @@ plunging\_radial\_solutions\_complex ==================================== -.. currentmodule:: kerrgeopy.plunging_solutions +.. currentmodule:: kerrgeopy.plunge .. autofunction:: plunging_radial_solutions_complex \ No newline at end of file diff --git a/docs/source/_autosummary/kerrgeopy.plunge.plunging_trajectory.rst b/docs/source/_autosummary/kerrgeopy.plunge.plunging_trajectory.rst new file mode 100644 index 0000000..baa5566 --- /dev/null +++ b/docs/source/_autosummary/kerrgeopy.plunge.plunging_trajectory.rst @@ -0,0 +1,6 @@ +plunging\_trajectory +==================== + +.. currentmodule:: kerrgeopy.plunge + +.. autofunction:: plunging_trajectory \ No newline at end of file diff --git a/docs/source/_autosummary/kerrgeopy.plunging_solutions.rst b/docs/source/_autosummary/kerrgeopy.plunge.rst similarity index 57% rename from docs/source/_autosummary/kerrgeopy.plunging_solutions.rst rename to docs/source/_autosummary/kerrgeopy.plunge.rst index c7fb6bb..d0a8e8c 100644 --- a/docs/source/_autosummary/kerrgeopy.plunging_solutions.rst +++ b/docs/source/_autosummary/kerrgeopy.plunge.rst @@ -1,7 +1,7 @@ -plunging\_solutions -=================== +plunge +====== -.. automodule:: kerrgeopy.plunging_solutions +.. automodule:: kerrgeopy.plunge @@ -15,15 +15,23 @@ plunging\_solutions :toctree: :template: custom-base-template.rst - plunging_mino_frequencies plunging_polar_solutions plunging_radial_integrals plunging_radial_solutions_complex + plunging_trajectory + .. rubric:: Classes + + .. autosummary:: + :toctree: + :template: custom-class-template.rst + + PlungingOrbit + diff --git a/docs/source/_autosummary/kerrgeopy.plunging_orbit.rst b/docs/source/_autosummary/kerrgeopy.plunging_orbit.rst deleted file mode 100644 index b84902a..0000000 --- a/docs/source/_autosummary/kerrgeopy.plunging_orbit.rst +++ /dev/null @@ -1,31 +0,0 @@ -plunging\_orbit -=============== - -.. automodule:: kerrgeopy.plunging_orbit - - - - - - - - - - - - .. rubric:: Classes - - .. autosummary:: - :toctree: - :template: custom-class-template.rst - - PlungingOrbit - - - - - - - - - diff --git a/docs/source/_autosummary/kerrgeopy.rst b/docs/source/_autosummary/kerrgeopy.rst index 8652a76..534b268 100644 --- a/docs/source/_autosummary/kerrgeopy.rst +++ b/docs/source/_autosummary/kerrgeopy.rst @@ -30,13 +30,10 @@ ~kerrgeopy.constants ~kerrgeopy.frequencies - ~kerrgeopy.frequencies_from_constants ~kerrgeopy.initial_conditions ~kerrgeopy.orbit - ~kerrgeopy.plunging_orbit - ~kerrgeopy.plunging_solutions + ~kerrgeopy.plunge ~kerrgeopy.spacetime - ~kerrgeopy.stable_orbit - ~kerrgeopy.stable_solutions + ~kerrgeopy.stable ~kerrgeopy.units diff --git a/docs/source/_autosummary/kerrgeopy.spacetime.KerrSpacetime.rst b/docs/source/_autosummary/kerrgeopy.spacetime.KerrSpacetime.rst index fbc66c5..586d41a 100644 --- a/docs/source/_autosummary/kerrgeopy.spacetime.KerrSpacetime.rst +++ b/docs/source/_autosummary/kerrgeopy.spacetime.KerrSpacetime.rst @@ -18,11 +18,9 @@ ~KerrSpacetime.__init__ ~KerrSpacetime.four_velocity - ~KerrSpacetime.inner_horizon ~KerrSpacetime.is_stable ~KerrSpacetime.metric ~KerrSpacetime.norm - ~KerrSpacetime.outer_horizon ~KerrSpacetime.separatrix diff --git a/docs/source/_autosummary/kerrgeopy.stable_orbit.StableOrbit.rst b/docs/source/_autosummary/kerrgeopy.stable.StableOrbit.rst similarity index 82% rename from docs/source/_autosummary/kerrgeopy.stable_orbit.StableOrbit.rst rename to docs/source/_autosummary/kerrgeopy.stable.StableOrbit.rst index 7b9c704..85250d3 100644 --- a/docs/source/_autosummary/kerrgeopy.stable_orbit.StableOrbit.rst +++ b/docs/source/_autosummary/kerrgeopy.stable.StableOrbit.rst @@ -1,7 +1,7 @@ StableOrbit =========== -.. currentmodule:: kerrgeopy.stable_orbit +.. currentmodule:: kerrgeopy.stable .. autoclass:: StableOrbit :members: @@ -20,12 +20,11 @@ ~StableOrbit.animate ~StableOrbit.constants_of_motion ~StableOrbit.four_velocity - ~StableOrbit.four_velocity_norm ~StableOrbit.from_constants ~StableOrbit.fundamental_frequencies ~StableOrbit.is_visible ~StableOrbit.mino_frequencies - ~StableOrbit.orbital_parameters + ~StableOrbit.numerical_four_velocity ~StableOrbit.plot ~StableOrbit.trajectory diff --git a/docs/source/_autosummary/kerrgeopy.stable_solutions.polar_solutions.rst b/docs/source/_autosummary/kerrgeopy.stable.polar_solutions.rst similarity index 60% rename from docs/source/_autosummary/kerrgeopy.stable_solutions.polar_solutions.rst rename to docs/source/_autosummary/kerrgeopy.stable.polar_solutions.rst index ecf287d..4bebe1d 100644 --- a/docs/source/_autosummary/kerrgeopy.stable_solutions.polar_solutions.rst +++ b/docs/source/_autosummary/kerrgeopy.stable.polar_solutions.rst @@ -1,6 +1,6 @@ polar\_solutions ================ -.. currentmodule:: kerrgeopy.stable_solutions +.. currentmodule:: kerrgeopy.stable .. autofunction:: polar_solutions \ No newline at end of file diff --git a/docs/source/_autosummary/kerrgeopy.stable_solutions.radial_solutions.rst b/docs/source/_autosummary/kerrgeopy.stable.radial_solutions.rst similarity index 61% rename from docs/source/_autosummary/kerrgeopy.stable_solutions.radial_solutions.rst rename to docs/source/_autosummary/kerrgeopy.stable.radial_solutions.rst index bcc336c..6f1f8da 100644 --- a/docs/source/_autosummary/kerrgeopy.stable_solutions.radial_solutions.rst +++ b/docs/source/_autosummary/kerrgeopy.stable.radial_solutions.rst @@ -1,6 +1,6 @@ radial\_solutions ================= -.. currentmodule:: kerrgeopy.stable_solutions +.. currentmodule:: kerrgeopy.stable .. autofunction:: radial_solutions \ No newline at end of file diff --git a/docs/source/_autosummary/kerrgeopy.stable_solutions.rst b/docs/source/_autosummary/kerrgeopy.stable.rst similarity index 51% rename from docs/source/_autosummary/kerrgeopy.stable_solutions.rst rename to docs/source/_autosummary/kerrgeopy.stable.rst index 9c336ce..c981d9e 100644 --- a/docs/source/_autosummary/kerrgeopy.stable_solutions.rst +++ b/docs/source/_autosummary/kerrgeopy.stable.rst @@ -1,7 +1,7 @@ -stable\_solutions -================= +stable +====== -.. automodule:: kerrgeopy.stable_solutions +.. automodule:: kerrgeopy.stable @@ -17,11 +17,20 @@ stable\_solutions polar_solutions radial_solutions + stable_trajectory + .. rubric:: Classes + + .. autosummary:: + :toctree: + :template: custom-class-template.rst + + StableOrbit + diff --git a/docs/source/_autosummary/kerrgeopy.stable.stable_trajectory.rst b/docs/source/_autosummary/kerrgeopy.stable.stable_trajectory.rst new file mode 100644 index 0000000..726c4c6 --- /dev/null +++ b/docs/source/_autosummary/kerrgeopy.stable.stable_trajectory.rst @@ -0,0 +1,6 @@ +stable\_trajectory +================== + +.. currentmodule:: kerrgeopy.stable + +.. autofunction:: stable_trajectory \ No newline at end of file diff --git a/docs/source/_autosummary/kerrgeopy.stable_orbit.rst b/docs/source/_autosummary/kerrgeopy.stable_orbit.rst deleted file mode 100644 index d4acd9c..0000000 --- a/docs/source/_autosummary/kerrgeopy.stable_orbit.rst +++ /dev/null @@ -1,31 +0,0 @@ -stable\_orbit -============= - -.. automodule:: kerrgeopy.stable_orbit - - - - - - - - - - - - .. rubric:: Classes - - .. autosummary:: - :toctree: - :template: custom-class-template.rst - - StableOrbit - - - - - - - - - diff --git a/docs/source/index.rst b/docs/source/index.rst index 422b92f..1626797 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -18,10 +18,6 @@ The library also provides a set of methods for computing constants of motion and :align: right :width: 45% -.. note:: - - This project is under active development. - .. _Installation: Installation @@ -61,8 +57,8 @@ API Reference :template: custom-module-template.rst :recursive: - ~kerrgeopy.stable_orbit.StableOrbit - ~kerrgeopy.plunging_orbit.PlungingOrbit + ~kerrgeopy.stable.StableOrbit + ~kerrgeopy.plunge.PlungingOrbit ~kerrgeopy.orbit.Orbit ~kerrgeopy.spacetime.KerrSpacetime ~kerrgeopy.constants diff --git a/docs/source/notebooks/Getting Started.ipynb b/docs/source/notebooks/Getting Started.ipynb index 2a25055..f3cc49a 100644 --- a/docs/source/notebooks/Getting Started.ipynb +++ b/docs/source/notebooks/Getting Started.ipynb @@ -33,12 +33,12 @@ "\n", "Note that $a$ and $x$ are restricted to values between -1 and 1, while $e$ is restricted to values between 0 and 1. Retrograde orbits are represented using a negative value for $a$ or for $x$. Polar orbits, marginally bound orbits, and orbits around an extreme Kerr black hole are not supported. \n", "\n", - "First, construct a [`StableOrbit`](stable_orbit.StableOrbit) using the four parameters described above." + "First, construct a [`StableOrbit`](stable.StableOrbit) using the four parameters described above." ] }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, "id": "b68a902b", "metadata": {}, "outputs": [], @@ -54,12 +54,12 @@ "id": "bb1022f0", "metadata": {}, "source": [ - "Plot the orbit from $\\lambda = 0$ to $\\lambda = 10$ using the [`plot()`](stable_orbit.StableOrbit.plot) method" + "Plot the orbit from $\\lambda = 0$ to $\\lambda = 10$ using the [`plot()`](stable.StableOrbit.plot) method" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 2, "id": "61a81092", "metadata": {}, "outputs": [ @@ -83,12 +83,12 @@ "id": "268961e0", "metadata": {}, "source": [ - "Next, compute the time, radial, polar and azimuthal components of the trajectory as a function of Mino time using the [`trajectory()`](stable_orbit.StableOrbit.trajectory) method. By default, the time and radial components of the trajectory are given in geometrized units and are normalized using $M$ so that they are dimensionless." + "Next, compute the time, radial, polar and azimuthal components of the trajectory as a function of Mino time using the [`trajectory()`](stable.StableOrbit.trajectory) method. By default, the time and radial components of the trajectory are given in geometrized units and are normalized using $M$ so that they are dimensionless." ] }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "id": "d682d1a5", "metadata": {}, "outputs": [], @@ -98,7 +98,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 4, "id": "c569de69", "metadata": {}, "outputs": [ @@ -108,7 +108,7 @@ "Text(0, 0.5, '$\\\\phi(\\\\lambda)$')" ] }, - "execution_count": 5, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, @@ -166,7 +166,7 @@ "id": "6fc7c3ec", "metadata": {}, "source": [ - "Use the [`constants_of_motion()`](stable_orbit.StableOrbit.constants_of_motion) method to compute the dimensionless energy, angular momentum and Carter constant. By default, constants of motion are given in geometrized units where $G=c=1$ and are scale-invariant, meaning that they are normalized according to the masses of the two bodies as follows:\n", + "Use the [`constants_of_motion()`](stable.StableOrbit.constants_of_motion) method to compute the dimensionless energy, angular momentum and Carter constant. By default, constants of motion are given in geometrized units where $G=c=1$ and are scale-invariant, meaning that they are normalized according to the masses of the two bodies as follows:\n", "\n", "\\begin{equation}\n", "\\mathcal{E} = \\frac{E}{\\mu}, \\quad \\mathcal{L} = \\frac{L}{\\mu M}, \\quad \\mathcal{Q} = \\frac{Q}{\\mu^2 M^2}\n", @@ -174,12 +174,12 @@ "\n", "Here, $M$ is the mass of the primary body and $\\mu$ is the mass of the secondary body. \n", "\n", - "Frequencies of motion can be computed in Mino time using the [`mino_frequencies()`](stable_orbit.StableOrbit.mino_frequencies) method and in Boyer-Lindquist time using the [`fundamental_frequencies()`](stable_orbit.StableOrbit.fundamental_frequencies) method. As with constants of motion, the frequencies returned by both methods are given in geometrized units and are normalized by $M$ so that they are dimensionless." + "Frequencies of motion can be computed in Mino time using the [`mino_frequencies()`](stable.StableOrbit.mino_frequencies) method and in Boyer-Lindquist time using the [`fundamental_frequencies()`](stable.StableOrbit.fundamental_frequencies) method. As with constants of motion, the frequencies returned by both methods are given in geometrized units and are normalized by $M$ so that they are dimensionless." ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 5, "id": "2b8cb67a", "metadata": {}, "outputs": [ @@ -292,7 +292,7 @@ "a = \\frac{J}{M^2}, \\quad \\mathcal{E} = \\frac{E}{\\mu}, \\quad \\mathcal{L} = \\frac{L}{\\mu M}, \\quad \\mathcal{Q} = \\frac{Q}{\\mu^2 M^2}\n", "\\end{equation}\n", "\n", - "Construct a [`PlungingOrbit`](plunging_orbit.PlungingOrbit) by passing in these four parameters." + "Construct a [`PlungingOrbit`](plunge.PlungingOrbit) by passing in these four parameters." ] }, { @@ -310,7 +310,7 @@ "id": "31281448", "metadata": {}, "source": [ - "As with stable orbits, the components of the trajectory can be computed using the [`trajectory()`](plunging_orbit.PlungingOrbit.trajectory) method" + "As with stable orbits, the components of the trajectory can be computed using the [`trajectory()`](plunge.PlungingOrbit.trajectory) method" ] }, { @@ -386,12 +386,12 @@ "source": [ "## Alternative Parametrizations\n", "\n", - "Use the [`from_constants()`](stable_orbit.StableOrbit.from_constants) class method to construct a [`StableOrbit`](stable_orbit.StableOrbit) from the spin parameter and constants of motion $(a,E,L,Q)$" + "Use the [`from_constants()`](stable.StableOrbit.from_constants) class method to construct a [`StableOrbit`](stable.StableOrbit) from the spin parameter and constants of motion $(a,E,L,Q)$" ] }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 9, "id": "d545e6a3", "metadata": {}, "outputs": [], @@ -409,7 +409,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 12, "id": "14142e83", "metadata": {}, "outputs": [], @@ -424,7 +424,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 13, "id": "9768d7dd", "metadata": {}, "outputs": [ @@ -434,7 +434,7 @@ "Text(0, 0.5, '$\\\\phi(\\\\lambda)$')" ] }, - "execution_count": 4, + "execution_count": 13, "metadata": {}, "output_type": "execute_result" }, @@ -494,7 +494,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.16" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/docs/source/notebooks/Graphics.ipynb b/docs/source/notebooks/Graphics.ipynb index b202c7c..199a94a 100644 --- a/docs/source/notebooks/Graphics.ipynb +++ b/docs/source/notebooks/Graphics.ipynb @@ -21,7 +21,7 @@ "id": "d714c093", "metadata": {}, "source": [ - "To plot a [`StableOrbit`](stable_orbit.StableOrbit), simply use the [`plot()`](stable_orbit.StableOrbit.plot) method and pass in the starting and ending Mino time." + "To plot a [`StableOrbit`](stable.StableOrbit), simply use the [`plot()`](stable.StableOrbit.plot) method and pass in the starting and ending Mino time." ] }, { @@ -121,7 +121,7 @@ "id": "100bca66", "metadata": {}, "source": [ - "To create an animation of a [`StableOrbit`](stable_orbit.StableOrbit) and save it as an mp4 file, use the [`animate()`](stable_orbit.StableOrbit.animate) method and pass in the starting and ending Mino time. Note that this method requires ffmpeg to be installed and may take several minutes to run depending on the length of the animation." + "To create an animation of a [`StableOrbit`](stable.StableOrbit) and save it as an mp4 file, use the [`animate()`](stable.StableOrbit.animate) method and pass in the starting and ending Mino time. Note that this method requires ffmpeg to be installed and may take several minutes to run depending on the length of the animation." ] }, { @@ -144,7 +144,7 @@ " \n", "\n", "\n", - "All of the options available when plotting are also available within the [`animate()`](stable_orbit.StableOrbit.animate) method. In addition, the `tail_length` option can be set to `\"short\"` or `\"none\"` to change the length of the orbital tail." + "All of the options available when plotting are also available within the [`animate()`](stable.StableOrbit.animate) method. In addition, the `tail_length` option can be set to `\"short\"` or `\"none\"` to change the length of the orbital tail." ] }, { diff --git a/docs/source/notebooks/Orbital Properties.ipynb b/docs/source/notebooks/Orbital Properties.ipynb index 2fea4ab..9897578 100644 --- a/docs/source/notebooks/Orbital Properties.ipynb +++ b/docs/source/notebooks/Orbital Properties.ipynb @@ -25,7 +25,7 @@ "\n", "Energy and angular momentum are generated by time translation symmetry and azimuthal symmetry respectively. The Carter constant is generated by a higher order symmetry of the Kerr metric which comes from a second order Killing tensor known as the Carter tensor. \n", "\n", - "These constants of motion can be computed for a [`StableOrbit`](stable_orbit.StableOrbit) using the [`constants_of_motion()`](stable_orbit.StableOrbit.constants_of_motion) method. By default, this method returns dimensionless quantities that are normalized as described in [Getting Started](orbital-properties)." + "These constants of motion can be computed for a [`StableOrbit`](stable.StableOrbit) using the [`constants_of_motion()`](stable.StableOrbit.constants_of_motion) method. By default, this method returns dimensionless quantities that are normalized as described in [Getting Started](orbital-properties)." ] }, { @@ -70,7 +70,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 2, "id": "7c3dfaa7", "metadata": {}, "outputs": [ @@ -111,7 +111,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 3, "id": "ff10d962", "metadata": {}, "outputs": [], @@ -121,7 +121,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 4, "id": "4356062f", "metadata": {}, "outputs": [], @@ -138,12 +138,12 @@ "source": [ "## Frequencies of Motion\n", "\n", - "To compute the frequencies of motion in Mino time for a [`StableOrbit`](stable_orbit.StableOrbit), use the [`mino_frequencies()`](stable_orbit.StableOrbit.mino_frequencies) method. This method returns the frequencies of motion for each of the three spatial coordinates $(\\Upsilon_r,\\Upsilon_\\theta,\\Upsilon_\\phi)$ along with $\\Gamma$, which measures the average rate at which observer time accumulates in Mino time. The frequencies of motion can be computed in observer time by simply dividing each of the Mino frequencies by $\\Gamma$, or by using the [`fundamental_frequencies()`](stable_orbit.StableOrbit.fundamental_frequencies) method. This method returns the Boyer-Lindquist frequencies $(\\Omega_r,\\Omega_\\theta,\\Omega_\\phi)$. By default, the frequencies returned by both methods are given in geometrized units where $G=c=1$ and are normalized by $M$ so that they are dimensionless." + "To compute the frequencies of motion in Mino time for a [`StableOrbit`](stable.StableOrbit), use the [`mino_frequencies()`](stable.StableOrbit.mino_frequencies) method. This method returns the frequencies of motion for each of the three spatial coordinates $(\\Upsilon_r,\\Upsilon_\\theta,\\Upsilon_\\phi)$ along with $\\Gamma$, which measures the average rate at which observer time accumulates in Mino time. The frequencies of motion can be computed in observer time by simply dividing each of the Mino frequencies by $\\Gamma$, or by using the [`fundamental_frequencies()`](stable.StableOrbit.fundamental_frequencies) method. This method returns the Boyer-Lindquist frequencies $(\\Omega_r,\\Omega_\\theta,\\Omega_\\phi)$. By default, the frequencies returned by both methods are given in geometrized units where $G=c=1$ and are normalized by $M$ so that they are dimensionless." ] }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 5, "id": "a60ab129", "metadata": {}, "outputs": [ @@ -206,7 +206,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 6, "id": "b79de692", "metadata": {}, "outputs": [ @@ -268,7 +268,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 7, "id": "2c53e56b", "metadata": {}, "outputs": [], @@ -279,7 +279,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 8, "id": "31f47fc2", "metadata": {}, "outputs": [], @@ -302,7 +302,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "id": "6573b1d7", "metadata": {}, "outputs": [ @@ -312,8 +312,9 @@ "4.462323556641868" ] }, + "execution_count": 9, "metadata": {}, - "output_type": "display_data" + "output_type": "execute_result" } ], "source": [ @@ -333,7 +334,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 10, "id": "e1a1e8da", "metadata": {}, "outputs": [ @@ -343,8 +344,9 @@ "True" ] }, + "execution_count": 10, "metadata": {}, - "output_type": "display_data" + "output_type": "execute_result" } ], "source": [ @@ -353,7 +355,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "id": "5d49a72f", "metadata": {}, "outputs": [ @@ -363,8 +365,9 @@ "False" ] }, + "execution_count": 11, "metadata": {}, - "output_type": "display_data" + "output_type": "execute_result" } ], "source": [ @@ -381,7 +384,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 12, "id": "d6062257", "metadata": {}, "outputs": [ @@ -391,7 +394,7 @@ "array([[4.46232356]])" ] }, - "execution_count": 2, + "execution_count": 12, "metadata": {}, "output_type": "execute_result" } @@ -449,7 +452,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "id": "8fb54304", "metadata": {}, "outputs": [], @@ -469,7 +472,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 14, "id": "c7efe933", "metadata": {}, "outputs": [ @@ -482,8 +485,9 @@ " [-0.396 , 0. , 0. , 26.37214 ]])" ] }, + "execution_count": 14, "metadata": {}, - "output_type": "display_data" + "output_type": "execute_result" } ], "source": [ @@ -500,7 +504,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 15, "id": "689e1f28", "metadata": {}, "outputs": [ @@ -510,8 +514,9 @@ "-1.0000000000000004" ] }, + "execution_count": 15, "metadata": {}, - "output_type": "display_data" + "output_type": "execute_result" } ], "source": [ @@ -548,7 +553,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.9" + "version": "3.10.12" } }, "nbformat": 4, diff --git a/docs/source/notebooks/Trajectory.ipynb b/docs/source/notebooks/Trajectory.ipynb index ad6867c..a4371b3 100644 --- a/docs/source/notebooks/Trajectory.ipynb +++ b/docs/source/notebooks/Trajectory.ipynb @@ -32,7 +32,7 @@ "source": [ "### Trajectory\n", "\n", - "To compute a trajectory, first construct a [`StableOrbit`](stable_orbit.StableOrbit) by specifying a set of orbital parameters $(a,p,e,x)$. Then use the [`trajectory()`](stable_orbit.StableOrbit.trajectory) method to compute the time, radial, polar and azimuthal components of the trajectory as a function of Mino time. Use the `initial_phases` option to set the initial phases $(q_{t_0},q_{r_0},q_{\\theta_0},q_{\\phi_0})$ of the orbit. Phases are defined as follows:\n", + "To compute a trajectory, first construct a [`StableOrbit`](stable.StableOrbit) by specifying a set of orbital parameters $(a,p,e,x)$. Then use the [`trajectory()`](stable.StableOrbit.trajectory) method to compute the time, radial, polar and azimuthal components of the trajectory as a function of Mino time. Use the `initial_phases` option to set the initial phases $(q_{t_0},q_{r_0},q_{\\theta_0},q_{\\phi_0})$ of the orbit. Phases are defined as follows:\n", "\n", "\\begin{equation}\n", "\\begin{aligned}\n", @@ -48,7 +48,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 36, "id": "3e870d54", "metadata": {}, "outputs": [], @@ -63,7 +63,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 37, "id": "220582c2", "metadata": {}, "outputs": [ @@ -73,7 +73,7 @@ "Text(0, 0.5, '$\\\\phi(\\\\lambda)$')" ] }, - "execution_count": 2, + "execution_count": 37, "metadata": {}, "output_type": "execute_result" }, @@ -126,12 +126,12 @@ "\n", "Plunging orbits are parametrized using the $(a,\\mathcal{E},\\mathcal{L},\\mathcal{Q})$ parametrization described in [Getting Started](plunging-orbits).\n", "\n", - "Construct a [`PlungingOrbit`](plunging_orbit.PlungingOrbit) by passing in these four parameters and use the [`trajectory()`](plunging_orbit.PlungingOrbit.trajectory) method to compute the trajectory. As with stable orbits, the `initial_phases` option sets the initial phases $(q_{t_0},q_{r_0},q_{\\theta_0},q_{\\phi_0})$ and the `distance_units` and `time_units` options can be used to specify units if `M` and `mu` are given." + "Construct a [`PlungingOrbit`](plunge.PlungingOrbit) by passing in these four parameters and use the [`trajectory()`](plunge.PlungingOrbit.trajectory) method to compute the trajectory. As with stable orbits, the `initial_phases` option sets the initial phases $(q_{t_0},q_{r_0},q_{\\theta_0},q_{\\phi_0})$ and the `distance_units` and `time_units` options can be used to specify units if `M` and `mu` are given." ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 38, "id": "59954e65", "metadata": {}, "outputs": [], @@ -143,7 +143,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 39, "id": "a2ca073a", "metadata": {}, "outputs": [ @@ -153,7 +153,7 @@ "Text(0, 0.5, '$\\\\phi(\\\\lambda)$')" ] }, - "execution_count": 4, + "execution_count": 39, "metadata": {}, "output_type": "execute_result" }, @@ -209,7 +209,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 40, "id": "345c2851", "metadata": {}, "outputs": [], @@ -221,7 +221,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 41, "id": "2ee806be", "metadata": {}, "outputs": [ @@ -231,7 +231,7 @@ "Text(0, 0.5, '$u^\\\\phi(\\\\lambda)$')" ] }, - "execution_count": 11, + "execution_count": 41, "metadata": {}, "output_type": "execute_result" }, @@ -285,7 +285,7 @@ }, { "cell_type": "code", - "execution_count": 34, + "execution_count": 42, "id": "d4b40f17", "metadata": {}, "outputs": [], @@ -297,7 +297,7 @@ }, { "cell_type": "code", - "execution_count": 35, + "execution_count": 43, "id": "3944a23f", "metadata": {}, "outputs": [ @@ -307,7 +307,7 @@ "Text(0, 0.5, '$u^\\\\phi(\\\\lambda)$')" ] }, - "execution_count": 35, + "execution_count": 43, "metadata": {}, "output_type": "execute_result" }, diff --git a/kerrgeopy/__init__.py b/kerrgeopy/__init__.py index c5e8fd8..c75fe2d 100644 --- a/kerrgeopy/__init__.py +++ b/kerrgeopy/__init__.py @@ -1,12 +1,12 @@ """ Python package for computing plunging and non-plunging geodesics in Kerr spacetime. """ -__all__ = ["units","constants", "frequencies", "stable_orbit", "plunging_orbit","initial_conditions"] +__all__ = ["units","constants", "frequencies","initial_conditions"] from kerrgeopy import * from kerrgeopy.frequencies import * -from kerrgeopy.constants import * from kerrgeopy.initial_conditions import * -from kerrgeopy.stable_orbit import StableOrbit -from kerrgeopy.plunging_orbit import PlungingOrbit +from kerrgeopy.constants import * +from kerrgeopy.stable import StableOrbit +from kerrgeopy.plunge import PlungingOrbit from kerrgeopy.orbit import Orbit from kerrgeopy.spacetime import KerrSpacetime \ No newline at end of file diff --git a/kerrgeopy/constants.py b/kerrgeopy/constants.py index 1260eac..b9ad26b 100644 --- a/kerrgeopy/constants.py +++ b/kerrgeopy/constants.py @@ -1,5 +1,5 @@ """ -Module containing functions for computing the constants of motion for bound orbits using the :math:`(a,p,e,x)` parametrization. +Module containing functions for computing the constants of motion for an orbit in Kerr Spacetime. Constants of motion are computed using the method described in Appendix B of `Schmidt `_. """ from numpy import sign, sqrt, copysign, pi, nan, inf @@ -7,6 +7,115 @@ from .units import * from scipy.interpolate import RectBivariateSpline import numpy as np +from numpy.polynomial import Polynomial + +def stable_radial_roots(a,p,e,x,constants=None): + """ + Computes the radial roots for a stable bound orbit as defined in equation 10 of `Fujita and Hikida `_. + Roots are given in decreasing order. + + :param a: dimensionless spin of the black hole + :type a: double + :param p: orbital semi-latus rectum + :type p: double + :param e: orbital eccentricity + :type e: double + :param x: cosine of the orbital inclination + :type x: tuple(double, double, double, double) + :param constants: dimensionless constants of motion for the orbit + :type constants: tuple(double, double, double) + + :return: tuple containing the four roots of the radial equation + :rtype: tuple(double, double, double, double) + """ + if constants is None: constants = constants_of_motion(a,p,e,x) + E, L, Q = constants + + r1 = p/(1-e) + r2 = p/(1+e) + + A_plus_B = 2/(1-E**2)-r1-r2 + AB = a**2*Q/(r1*r2*(1-E**2)) + + r3 = (A_plus_B+sqrt(A_plus_B**2-4*AB))/2 + r4 = AB/r3 + + return r1, r2, r3, r4 + + +def plunging_radial_roots(a,E,L,Q): + """ + Computes the radial roots for a plunging orbit. + If all roots are real, roots are sorted such that the motion is between r1 and r2 and roots are otherwise in decreasing order. + If there are two complex roots, r1 < r2 are real and r3/r4 are complex conjugates. + + :param a: dimensionless spin parameter + :type a: double + :param E: dimensionless energy + :type E: double + :param L: dimensionless angular momentum + :type L: double + :param Q: dimensionless Carter constant + :type Q: double + + :return: tuple of radial roots + :rtype: tuple(double,double,double,double) + """ + # standard form of the radial polynomial R(r) + R = Polynomial([-a**2*Q, 2*L**2+2*Q+2*a**2*E**2-4*a*E*L, a**2*E**2-L**2-Q-a**2, 2, E**2-1]) + roots = R.roots() + # get the real roots and the complex roots + real_roots = np.sort(np.real(roots[np.isreal(roots)])) + complex_roots = roots[np.iscomplex(roots)] + + r_minus = 1-sqrt(1-a**2) + + # if there are 4 real roots, by convention r4 < r3 < r2 < r1 (consistent with stable orbits) + if len(real_roots) == 4: + # if there are three roots outside the event horizon swap r1/r3 and r2/r4 + if real_roots[1] > r_minus: + r1 = real_roots[1] + r2 = real_roots[0] + r3 = real_roots[3] + r4 = real_roots[2] + else: + r4, r3, r2, r1 = real_roots + + # in the case of two complex roots, r1 < r2 are real and r3/r4 are complex conjugates + elif len(real_roots) == 2: + r1, r2 = real_roots + r3, r4 = complex_roots + + return r1, r2, r3, r4 + +def stable_polar_roots(a,p,e,x,constants=None): + r""" + Computes the polar roots for a stable bound orbit as defined in equation 10 of `Fujita and Hikida `_. + Roots are given in increasing order. + + :param a: dimensionless spin of the black hole + :type a: double + :param p: orbital semi-latus rectum + :type p: double + :param e: orbital eccentricity + :type e: double + :param x: cosine of the orbital inclination + :type x: tuple(double, double, double) + :param constants: dimensionless constants of motion for the orbit + :type constants: tuple(double, double, double) + + :return: tuple of roots in the form :math:`(z_-, z_+)` + :rtype: tuple(double, double, double, double) + """ + if constants is None: constants = constants_of_motion(a,p,e,x) + E, L, Q = constants + epsilon0 = a**2*(1-E**2)/L**2 + z_minus = 1-x**2 + #z_plus = a**2*(1-E**2)/(L**2*epsilon0)+1/(epsilon0*(1-z_minus)) + # simplified using definition of carter constant + z_plus = nan if a == 0 else 1+1/(epsilon0*(1-z_minus)) + + return z_minus, z_plus def _coefficients(r,a,x): """ @@ -225,6 +334,40 @@ def constants_of_motion(a,p,e,x): Q = carter_constant(a,p,e,x,E,L) return E, L, Q +def apex_from_constants(a,E,L,Q): + r""" + Computes the orbital parameters :math:`(a,p,e,x)` for a stable bound orbit with the given constants of motion + + :param a: spin parameter + :type a: double + :param E: dimensionless energy + :type E: double + :param L: dimensionless angular momentum + :type L: double + :param Q: dimensionless Carter constant + :type Q: double + + :return: tuple of orbital parameters :math:`(a,p,e,x)` + :rtype: tuple(double,double,double,double) + """ + # radial polynomial written in terms of z = cos^2(theta) + R = Polynomial([-a**2*Q, 2*L**2+2*Q+2*a**2*E**2-4*a*E*L, a**2*E**2-L**2-Q-a**2, 2, E**2-1]) + radial_roots = R.roots() + # numpy returns roots in increasing order + r4, r3, r2, r1 = radial_roots + + # polar polynomial written in terms of z = cos^2(theta) + Z = Polynomial([Q,-(Q+a**2*(1-E**2)+L**2),a**2*(1-E**2)]) + polar_roots = Z.roots() + if a == 0: polar_roots = [polar_roots[0], polar_roots[0]] + z_minus, z_plus = polar_roots + + p = 2*r1*r2/(r1+r2) + e = (r1-r2)/(r1+r2) + x = np.sign(L)*sqrt(1-z_minus) + + return a, p, e, x + def _S_polar(p,a,e): """ Separatrix polynomial for a polar orbit from equation 37 in `Stein and Warburton `_ diff --git a/kerrgeopy/frequencies.py b/kerrgeopy/frequencies.py index 153c58a..8079a7d 100644 --- a/kerrgeopy/frequencies.py +++ b/kerrgeopy/frequencies.py @@ -1,64 +1,46 @@ """ -Module containing functions for computing frequencies of motion for bound orbits using the :math:`(a,p,e,x)` parametrization. +Module containing functions for computing frequencies of motion for orbits in Kerr spacetime. Frequencies are computed using the method derived in `Fujita and Hikida `_ """ from .constants import _standardize_params -from .frequencies_from_constants import * -from .frequencies_from_constants import _ellippi +from .constants import * +from scipy.special import ellipk, ellipe, elliprj, elliprf +from numpy import sin, cos, sqrt, pi, arcsin, floor, where -def _radial_roots(a,p,e,constants): - """ - Computes r1, r2, r3 and r4 as defined in equation 10 of `Fujita and Hikida `_ - - :param a: dimensionless spin of the black hole - :type a: double - :param p: orbital semi-latus rectum - :type p: double - :param e: orbital eccentricity - :type e: double - :param x: cosine of the orbital inclination - :type x: tuple(double, double, double, double) - :param constants: dimensionless constants of motion for the orbit - :type constants: tuple(double, double, double) +def _ellippi(n,m): + r""" + Complete elliptic integral of the third kind defined as :math:`\Pi(n,m) = \int_0^{\frac{\pi}{2}} \frac{d\theta}{(1-n\sin^2{\theta})\sqrt{1-m\sin^2{\theta}}}` + + :type n: double + :type m: double - :return: tuple containing the four roots of the radial equation - :rtype: tuple(double, double, double, double) + :rtype: double """ - E, L, Q = constants - - r1 = p/(1-e) - r2 = p/(1+e) - - A_plus_B = 2/(1-E**2)-r1-r2 - AB = a**2*Q/(r1*r2*(1-E**2)) - - r3 = (A_plus_B+sqrt(A_plus_B**2-4*AB))/2 - r4 = AB/r3 - - return r1, r2, r3, r4 + # Note: sign of n is reversed from the definition in Fujita and Hikida -def _polar_roots(a,x,constants): + # formula from https://en.wikipedia.org/wiki/Carlson_symmetric_form + return elliprf(0,1-m,1)+1/3*n*elliprj(0,1-m,1,1-n) + +def _ellippiinc(phi,n,m): r""" - Computes z_minus and z_plus as defined in equation 10 of `Fujita and Hikida `_ + Incomplete elliptic integral of the third kind defined as :math:`\Pi(\phi,n,m) = \int_0^{\phi} \frac{1}{1-n\sin^2\theta}\frac{1}{\sqrt{1-m\sin^2\theta}}d\theta`. - :param a: dimensionless spin of the black hole - :type a: double - :param x: cosine of the orbital inclination - :type x: tuple(double, double, double) - :param constants: dimensionless constants of motion for the orbit - :type constants: tuple(double, double, double) + :type phi: double + :type n: double + :type m: double - :return: tuple of roots in the form :math:`(z_-, z_+)` - :rtype: tuple(double, double, double, double) + :rtype: double """ - E, L, Q = constants - epsilon0 = a**2*(1-E**2)/L**2 - z_minus = 1-x**2 - #z_plus = a**2*(1-E**2)/(L**2*epsilon0)+1/(epsilon0*(1-z_minus)) - # simplified using definition of carter constant - z_plus = nan if a == 0 else 1+1/(epsilon0*(1-z_minus)) - - return z_minus, z_plus + # Note: sign of n is reversed from the definition in Fujita and Hikida + + # count the number of half periods + num_cycles = floor(phi/(pi/2)) + # map phi to [0,pi/2] + phi = abs(arcsin(sin(phi))) + # formula from https://en.wikipedia.org/wiki/Carlson_symmetric_form + integral = sin(phi)*elliprf(cos(phi)**2,1-m*sin(phi)**2,1)+1/3*n*sin(phi)**3*elliprj(cos(phi)**2,1-m*sin(phi)**2,1,1-n*sin(phi)**2) + + return where(num_cycles % 2 == 0, num_cycles*_ellippi(n,m)+integral, (num_cycles+1)*_ellippi(n,m)-integral) def r_frequency(a,p,e,x,constants=None): """ @@ -89,7 +71,7 @@ def r_frequency(a,p,e,x,constants=None): if constants is None: constants = constants_of_motion(a,p,e,x) E, L, Q = constants - r1,r2,r3,r4 = _radial_roots(a,p,e,constants) + r1,r2,r3,r4 = stable_radial_roots(a,p,e,x,constants) # equation 13 k_r = sqrt((r1-r2)*(r3-r4)/((r1-r3)*(r2-r4))) # equation 15 @@ -128,7 +110,7 @@ def theta_frequency(a,p,e,x,constants=None): # compute constants if not provided if constants is None: constants = constants_of_motion(a,p,e,x) E, L, Q = constants - z_minus, z_plus = _polar_roots(a,x,constants) + z_minus, z_plus = stable_polar_roots(a,p,e,x,constants) # equation 13 k_theta = sqrt(z_minus/z_plus) @@ -174,8 +156,8 @@ def phi_frequency(a,p,e,x,constants=None,upsilon_r=None,upsilon_theta=None): # compute constants if they are not passed in if constants is None: constants = constants_of_motion(a,p,e,x) E, L, Q = constants - r1,r2,r3,r4 = _radial_roots(a,p,e,constants) - z_minus, z_plus = _polar_roots(a,x,constants) + r1,r2,r3,r4 = stable_radial_roots(a,p,e,x,constants) + z_minus, z_plus = stable_polar_roots(a,p,e,x,constants) # compute frequencies if they are not passed in if upsilon_r is None: upsilon_r = r_frequency(a,p,e,x,constants) @@ -195,10 +177,10 @@ def phi_frequency(a,p,e,x,constants=None,upsilon_r=None,upsilon_theta=None): k_r = sqrt((r1-r2)*(r3-r4)/((r1-r3)*(r2-r4))) # equation 21 - return 2*upsilon_theta/(pi*sqrt(e0zp))*_ellippi(z_minus,k_theta) \ + return 2*upsilon_theta/(pi*sqrt(e0zp))*_ellippi(z_minus,k_theta**2) \ + 2*a*upsilon_r/(pi*(r_plus-r_minus)*sqrt((1-E**2)*(r1-r3)*(r2-r4))) \ - * ((2*E*r_plus-a*L)/(r3-r_plus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_plus)*_ellippi(h_plus,k_r)) \ - - (2*E*r_minus-a*L)/(r3-r_minus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_minus)*_ellippi(h_minus,k_r)) + * ((2*E*r_plus-a*L)/(r3-r_plus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_plus)*_ellippi(h_plus,k_r**2)) \ + - (2*E*r_minus-a*L)/(r3-r_minus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_minus)*_ellippi(h_minus,k_r**2)) ) def gamma(a,p,e,x,constants=None,upsilon_r=None,upsilon_theta=None): @@ -237,8 +219,8 @@ def gamma(a,p,e,x,constants=None,upsilon_r=None,upsilon_theta=None): # compute constants if they are not passed in if constants is None: constants = constants_of_motion(a,p,e,x) E, L, Q = constants - r1,r2,r3,r4 = _radial_roots(a,p,e,constants) - z_minus, z_plus = _polar_roots(a,x,constants) + r1,r2,r3,r4 = stable_radial_roots(a,p,e,x,constants) + z_minus, z_plus = stable_polar_roots(a,p,e,x,constants) epsilon0 = a**2*(1-E**2)/L**2 # simplified form of a**2*sqrt(z_plus/epsilon0) a2sqrt_zp_over_e0 = L**2/((1-E**2)*sqrt(1-z_minus)) if a == 0 else a**2*z_plus/sqrt(epsilon0*z_plus) @@ -261,10 +243,10 @@ def gamma(a,p,e,x,constants=None,upsilon_r=None,upsilon_theta=None): # equation 21 return 4*E + 2*a2sqrt_zp_over_e0*E*upsilon_theta*(ellipk(k_theta**2)-ellipe(k_theta**2))/(pi*L) \ + 2*upsilon_r/(pi*sqrt((1-E**2)*(r1-r3)*(r2-r4))) \ - * (E/2*((r3*(r1+r2+r3)-r1*r2)*ellipk(k_r**2) + (r2-r3)*(r1+r2+r3+r4)*_ellippi(h_r,k_r)+ (r1-r3)*(r2-r4)*ellipe(k_r**2)) \ - + 2*E*(r3*ellipk(k_r**2)+(r2-r3)*_ellippi(h_r,k_r))\ - +2/(r_plus-r_minus)*(((4*E-a*L)*r_plus-2*a**2*E)/(r3-r_plus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_plus)*_ellippi(h_plus,k_r)) \ - -((4*E-a*L)*r_minus-2*a**2*E)/(r3-r_minus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_minus)*_ellippi(h_minus,k_r)) + * (E/2*((r3*(r1+r2+r3)-r1*r2)*ellipk(k_r**2) + (r2-r3)*(r1+r2+r3+r4)*_ellippi(h_r,k_r**2)+ (r1-r3)*(r2-r4)*ellipe(k_r**2)) \ + + 2*E*(r3*ellipk(k_r**2)+(r2-r3)*_ellippi(h_r,k_r**2))\ + +2/(r_plus-r_minus)*(((4*E-a*L)*r_plus-2*a**2*E)/(r3-r_plus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_plus)*_ellippi(h_plus,k_r**2)) \ + -((4*E-a*L)*r_minus-2*a**2*E)/(r3-r_minus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_minus)*_ellippi(h_minus,k_r**2)) ) ) @@ -314,4 +296,263 @@ def fundamental_frequencies(a,p,e,x): upsilon_phi = phi_frequency(a,p,e,x,constants,upsilon_r,upsilon_theta) Gamma = gamma(a,p,e,x,constants,upsilon_r,upsilon_theta) + return upsilon_r/Gamma, abs(upsilon_theta)/Gamma, upsilon_phi/Gamma + +def plunging_mino_frequencies(a,E,L,Q): + r""" + Computes the radial and polar mino frequencies for a plunging orbit. + + :param a: dimensionless spin parameter + :type a: double + :param E: dimensionless energy + :type E: double + :param L: dimensionless angular momentum + :type L: double + :param Q: dimensionless Carter constant + :type Q: double + + :return: radial and polar mino frequencies :math:`(\Upsilon_r,\Upsilon_\theta)` + :rtype: tuple(double,double) + """ + + radial_roots = plunging_radial_roots(a,E,L,Q) + r1, r2, r3, r4 = radial_roots + rho_r = np.real(r3) + rho_i = np.imag(r4) + if np.iscomplex(radial_roots[3]): + # equation 42 + A = sqrt((r2-rho_r)**2+rho_i**2) + B = sqrt((r1-rho_r)**2+rho_i**2) + k_r = sqrt(abs(((r2-r1)**2-(A-B)**2)/(4*A*B))) + + # equation 52 + upsilon_r = pi*sqrt(A*B*(1-E**2))/(2*ellipk(k_r**2)) + + # equation 14 + z1 = sqrt(1/2*(1+(L**2+Q)/(a**2*(1-E**2))-sqrt((1+(L**2+Q)/(a**2*(1-E**2)))**2-4*Q/(a**2*(1-E**2))))) + z2 = sqrt(a**2*(1-E**2)/2*(1+(L**2+Q)/(a**2*(1-E**2))+sqrt((1+(L**2+Q)/(a**2*(1-E**2)))**2-4*Q/(a**2*(1-E**2))))) + k_theta = a*sqrt(1-E**2)*z1/z2 + + # equation 53 + upsilon_theta = pi/2*z2/ellipk(k_theta**2) + else: + # polar polynomial written in terms of z = cos^2(theta) + Z = Polynomial([Q,-(Q+a**2*(1-E**2)+L**2),a**2*(1-E**2)]) + polar_roots = Z.roots() + if a == 0: polar_roots = [polar_roots[0], polar_roots[0]] + constants = (E,L,Q) + upsilon_r, upsilon_theta, upsilon_phi, gamma = _mino_frequencies_from_constants(a,constants,radial_roots,polar_roots) + + + return upsilon_r, upsilon_theta + +def _r_frequency_from_constants(constants,radial_roots): + """ + Computes the frequency of motion in r in Mino time. + + :param constants: dimensionless constants of motion for the orbit in the form :math:`(E,L,Q)` + :type constants: tuple(double, double, double) + :param radial_roots: tuple containing the four roots of the radial equation :math:`(r_1, r_2, r_3, r_4)`. + :type radial_roots: tuple(double, double, double, double) + + :rtype: double + """ + E, L, Q = constants + + r1,r2,r3,r4 = radial_roots + # equation 13 + k_r = sqrt((r1-r2)*(r3-r4)/((r1-r3)*(r2-r4))) + # equation 15 + return pi*sqrt((1-E**2)*(r1-r3)*(r2-r4))/(2*ellipk(k_r**2)) + +def _theta_frequency_from_constants(a,constants,radial_roots,polar_roots): + """ + Computes the frequency of motion in theta in Mino time. + + :param a: dimensionless spin of the black hole + :type a: double + :param constants: dimensionless constants of motion for the orbit in the form :math:`(E,L,Q)` + :type constants: tuple(double, double, double) + :param radial_roots: tuple containing the four roots of the radial polynomial :math:`(r_1, r_2, r_3, r_4)`. + :type radial_roots: tuple(double, double, double, double) + :param polar_roots: tuple containing the roots of the polar equation :math:`(z_-, z_+)` + :type polar_roots: tuple(double, double) + + :rtype: double + """ + r1, r2, r3, r4 = radial_roots + z_minus, z_plus = polar_roots + E, L, Q = constants + + # Schwarzschild case + if a == 0: + p = 2*r1*r2/(r1+r2) + e = (r1-r2)/(r1+r2) + return p/sqrt(-3-e**2+p)*(sign(L) if e == 0 else 1) + + + # equation 13 + k_theta = sqrt(z_minus/z_plus) + # simplified form of epsilon0*z_plus + e0zp = (a**2*(1-E**2)*(1-z_minus)+L**2)/(L**2*(1-z_minus)) + + # equation 15 + return pi*L*sqrt(e0zp)/(2*ellipk(k_theta**2)) + +def _phi_frequency_from_constants(a,constants,radial_roots,polar_roots,upsilon_r=None,upsilon_theta=None): + """ + Computes the frequency of motion in phi in Mino time. + + :param a: dimensionless spin of the black hole + :type a: double + :param constants: dimensionless constants of motion for the orbit in the form :math:`(E,L,Q)` + :type constants: tuple(double, double, double) + :param radial_roots: tuple containing the four roots of the radial polynomial :math:`(r_1, r_2, r_3, r_4)`. + :type radial_roots: tuple(double, double, double, double) + :param polar_roots: tuple containing the roots of the polar equation :math:`(z_-, z_+)` + :type polar_roots: tuple(double, double) + :param upsilon_r: Mino frequency of motion in r can be passed in to speed computation if it is already known + :type upsilon_r: double, optional + :param upsilon_theta: Mino frequency of motion in theta can be passed in to speed computation if it is already known + :type upsilon_theta: double, optional + + :rtype: double + """ + + E, L, Q = constants + r1,r2,r3,r4 = radial_roots + z_minus, z_plus = polar_roots + + # Schwarzschild case + if a == 0: + p = 2*r1*r2/(r1+r2) + e = (r1-r2)/(r1+r2) + return sign(L)*p/sqrt(-3-e**2+p) + + # compute frequencies if they are not passed in + if upsilon_r is None: upsilon_r = _r_frequency_from_constants(a,constants,radial_roots) + if upsilon_theta is None: upsilon_theta = _theta_frequency_from_constants(a,constants,radial_roots,polar_roots) + + # simplified form of epsilon0*z_plus + e0zp = (a**2*(1-E**2)*(1-z_minus)+L**2)/(L**2*(1-z_minus)) + + r_plus = 1+sqrt(1-a**2) + r_minus = 1-sqrt(1-a**2) + + h_plus = (r1-r2)*(r3-r_plus)/((r1-r3)*(r2-r_plus)) + h_minus = (r1-r2)*(r3-r_minus)/((r1-r3)*(r2-r_minus)) + + # equation 13 + k_theta = sqrt(z_minus/z_plus) + k_r = sqrt((r1-r2)*(r3-r4)/((r1-r3)*(r2-r4))) + + # equation 21 + return 2*upsilon_theta/(pi*sqrt(e0zp))*_ellippi(z_minus,k_theta**2) \ + + 2*a*upsilon_r/(pi*(r_plus-r_minus)*sqrt((1-E**2)*(r1-r3)*(r2-r4))) \ + * ((2*E*r_plus-a*L)/(r3-r_plus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_plus)*_ellippi(h_plus,k_r**2)) \ + - (2*E*r_minus-a*L)/(r3-r_minus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_minus)*_ellippi(h_minus,k_r**2)) + ) + +def _gamma_from_constants(a,constants,radial_roots,polar_roots,upsilon_r=None,upsilon_theta=None): + """ + Computes the average rate at which observer time accumulates in Mino time. + + :param a: dimensionless spin of the black hole + :type a: double + :param constants: dimensionless constants of motion for the orbit in the form :math:`(E,L,Q)` + :type constants: tuple(double, double, double) + :param radial_roots: tuple containing the four roots of the radial polynomial :math:`(r_1, r_2, r_3, r_4)`. + :type radial_roots: tuple(double, double, double, double) + :param polar_roots: tuple containing the roots of the polar equation :math:`(z_-, z_+)` + :type polar_roots: tuple(double, double) + :param upsilon_r: Mino frequency of motion in r can be passed in to speed computation if it is already known + :type upsilon_r: double, optional + :param upsilon_theta: Mino frequency of motion in theta can be passed in to speed computation if it is already known + :type upsilon_theta: double, optional + + :rtype: double + """ + r1,r2,r3,r4 = radial_roots + z_minus, z_plus = polar_roots + + e = (r1-r2)/(r1+r2) + # marginally bound case + if e == 1: + return inf + + E, L, Q = constants + + epsilon0 = a**2*(1-E**2)/L**2 + # simplified form of a**2*sqrt(z_plus/epsilon0) + a2sqrt_zp_over_e0 = L**2/((1-E**2)*sqrt(1-z_minus)) if a == 0 else a**2*z_plus/sqrt(epsilon0*z_plus) + + # compute frequencies if they are not passed in + if upsilon_r is None: upsilon_r = _r_frequency_from_constants(a,constants,radial_roots) + if upsilon_theta is None: upsilon_theta = _theta_frequency_from_constants(a,constants,radial_roots,polar_roots) + + r_plus = 1+sqrt(1-a**2) + r_minus = 1-sqrt(1-a**2) + + h_r = (r1-r2)/(r1-r3) + h_plus = (r1-r2)*(r3-r_plus)/((r1-r3)*(r2-r_plus)) + h_minus = (r1-r2)*(r3-r_minus)/((r1-r3)*(r2-r_minus)) + + # equation 13 + k_theta = 0 if a == 0 else sqrt(z_minus/z_plus) + k_r = sqrt((r1-r2)*(r3-r4)/((r1-r3)*(r2-r4))) + + # equation 21 + return 4*E + 2*a2sqrt_zp_over_e0*E*upsilon_theta*(ellipk(k_theta**2)-ellipe(k_theta**2))/(pi*L) \ + + 2*upsilon_r/(pi*sqrt((1-E**2)*(r1-r3)*(r2-r4))) \ + * (E/2*((r3*(r1+r2+r3)-r1*r2)*ellipk(k_r**2) + (r2-r3)*(r1+r2+r3+r4)*_ellippi(h_r,k_r**2)+ (r1-r3)*(r2-r4)*ellipe(k_r**2)) \ + + 2*E*(r3*ellipk(k_r**2)+(r2-r3)*_ellippi(h_r,k_r**2))\ + +2/(r_plus-r_minus)*(((4*E-a*L)*r_plus-2*a**2*E)/(r3-r_plus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_plus)*_ellippi(h_plus,k_r**2)) \ + -((4*E-a*L)*r_minus-2*a**2*E)/(r3-r_minus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_minus)*_ellippi(h_minus,k_r**2)) + ) + ) + +def _mino_frequencies_from_constants(a,constants,radial_roots,polar_roots): + r""" + Computes frequencies of orbital motion in Mino time. + + :param a: dimensionless spin of the black hole (must satisfy -1 < a < 1) + :type a: double + :param constants: dimensionless constants of motion for the orbit in the form :math:`(E,L,Q)` + :type constants: tuple(double, double, double) + :param radial_roots: tuple containing the four roots of the radial polynomial :math:`(r_1, r_2, r_3, r_4)`. + :type radial_roots: tuple(double, double, double, double) + :param polar_roots: tuple containing the roots of the polar equation :math:`(z_-, z_+)` + :type polar_roots: tuple(double, double) + + :return: tuple of frequencies in the form :math:`(\Upsilon_r, \Upsilon_\theta, \Upsilon_\phi, \Gamma)` + :rtype: tuple(double, double, double, double) + """ + upsilon_r = _r_frequency_from_constants(constants,radial_roots) + upsilon_theta = _theta_frequency_from_constants(a,constants,radial_roots,polar_roots) + upsilon_phi = _phi_frequency_from_constants(a,constants,radial_roots,polar_roots,upsilon_r,upsilon_theta) + Gamma = _gamma_from_constants(a,constants,radial_roots,polar_roots,upsilon_r,upsilon_theta) + + return upsilon_r, abs(upsilon_theta), upsilon_phi, Gamma + +def _fundamental_frequencies_from_constants(a,constants,radial_roots,polar_roots): + r""" + Computes frequencies of orbital motion in Boyer-Lindquist time. + + :param a: dimensionless spin of the black hole (must satisfy -1 < a < 1) + :type a: double + :param constants: dimensionless constants of motion for the orbit in the form :math:`(E,L,Q)` + :type constants: tuple(double, double, double) + :param radial_roots: tuple containing the four roots of the radial polynomial :math:`(r_1, r_2, r_3, r_4)`. + :type radial_roots: tuple(double, double, double, double) + :param polar_roots: tuple containing the roots of the polar equation :math:`(z_-, z_+)` + :type polar_roots: tuple(double, double) + + :return: tuple of frequencies in the form :math:`(\Omega_r, \Omega_\theta, \Omega_\phi)` + :rtype: tuple(double, double, double) + """ + upsilon_r = _r_frequency_from_constants(constants,radial_roots) + upsilon_theta = _theta_frequency_from_constants(a,constants,radial_roots,polar_roots) + upsilon_phi = _phi_frequency_from_constants(a,constants,radial_roots,polar_roots,upsilon_r,upsilon_theta) + Gamma = _gamma_from_constants(a,constants,radial_roots,polar_roots,upsilon_r,upsilon_theta) + return upsilon_r/Gamma, abs(upsilon_theta)/Gamma, upsilon_phi/Gamma \ No newline at end of file diff --git a/kerrgeopy/frequencies_from_constants.py b/kerrgeopy/frequencies_from_constants.py deleted file mode 100644 index 3a9a693..0000000 --- a/kerrgeopy/frequencies_from_constants.py +++ /dev/null @@ -1,231 +0,0 @@ -""" -Module containing functions for computing frequencies of motion from the spin parameter, constants of motion, and radial and polar roots. -Frequencies are computed using the method derived in `Fujita and Hikida `_ -""" -from scipy.special import ellipk, ellipe, elliprj, elliprf -from .constants import * - -def _ellippi(n,k): - r""" - Complete elliptic integral of the third kind defined as :math:`\Pi(n,k) = \int_0^{\frac{\pi}{2}} \frac{d\theta}{(1-n\sin^2{\theta})\sqrt{1-k^2\sin^2{\theta}}}` - - :type n: double - :type k: double - - :rtype: double - """ - # Note: sign of n is reversed from the definition in Fujita and Hikida - - # formula from https://en.wikipedia.org/wiki/Carlson_symmetric_form - return elliprf(0,1-k**2,1)+1/3*n*elliprj(0,1-k**2,1,1-n) - -def r_frequency_from_constants(constants,radial_roots): - """ - Computes the frequency of motion in r in Mino time. - - :param constants: dimensionless constants of motion for the orbit in the form :math:`(E,L,Q)` - :type constants: tuple(double, double, double) - :param radial_roots: tuple containing the four roots of the radial equation :math:`(r_1, r_2, r_3, r_4)`. - :type radial_roots: tuple(double, double, double, double) - - :rtype: double - """ - E, L, Q = constants - - r1,r2,r3,r4 = radial_roots - # equation 13 - k_r = sqrt((r1-r2)*(r3-r4)/((r1-r3)*(r2-r4))) - # equation 15 - return pi*sqrt((1-E**2)*(r1-r3)*(r2-r4))/(2*ellipk(k_r**2)) - -def theta_frequency_from_constants(a,constants,radial_roots,polar_roots): - """ - Computes the frequency of motion in theta in Mino time. - - :param a: dimensionless spin of the black hole - :type a: double - :param constants: dimensionless constants of motion for the orbit in the form :math:`(E,L,Q)` - :type constants: tuple(double, double, double) - :param radial_roots: tuple containing the four roots of the radial polynomial :math:`(r_1, r_2, r_3, r_4)`. - :type radial_roots: tuple(double, double, double, double) - :param polar_roots: tuple containing the roots of the polar equation :math:`(z_-, z_+)` - :type polar_roots: tuple(double, double) - - :rtype: double - """ - r1, r2, r3, r4 = radial_roots - z_minus, z_plus = polar_roots - E, L, Q = constants - - # Schwarzschild case - if a == 0: - p = 2*r1*r2/(r1+r2) - e = (r1-r2)/(r1+r2) - return p/sqrt(-3-e**2+p)*(sign(L) if e == 0 else 1) - - - # equation 13 - k_theta = sqrt(z_minus/z_plus) - # simplified form of epsilon0*z_plus - e0zp = (a**2*(1-E**2)*(1-z_minus)+L**2)/(L**2*(1-z_minus)) - - # equation 15 - return pi*L*sqrt(e0zp)/(2*ellipk(k_theta**2)) - -def phi_frequency_from_constants(a,constants,radial_roots,polar_roots,upsilon_r=None,upsilon_theta=None): - """ - Computes the frequency of motion in phi in Mino time. - - :param a: dimensionless spin of the black hole - :type a: double - :param constants: dimensionless constants of motion for the orbit in the form :math:`(E,L,Q)` - :type constants: tuple(double, double, double) - :param radial_roots: tuple containing the four roots of the radial polynomial :math:`(r_1, r_2, r_3, r_4)`. - :type radial_roots: tuple(double, double, double, double) - :param polar_roots: tuple containing the roots of the polar equation :math:`(z_-, z_+)` - :type polar_roots: tuple(double, double) - :param upsilon_r: Mino frequency of motion in r can be passed in to speed computation if it is already known - :type upsilon_r: double, optional - :param upsilon_theta: Mino frequency of motion in theta can be passed in to speed computation if it is already known - :type upsilon_theta: double, optional - - :rtype: double - """ - - E, L, Q = constants - r1,r2,r3,r4 = radial_roots - z_minus, z_plus = polar_roots - - # Schwarzschild case - if a == 0: - p = 2*r1*r2/(r1+r2) - e = (r1-r2)/(r1+r2) - return sign(L)*p/sqrt(-3-e**2+p) - - # compute frequencies if they are not passed in - if upsilon_r is None: upsilon_r = r_frequency_from_constants(a,constants,radial_roots) - if upsilon_theta is None: upsilon_theta = theta_frequency_from_constants(a,constants,radial_roots,polar_roots) - - # simplified form of epsilon0*z_plus - e0zp = (a**2*(1-E**2)*(1-z_minus)+L**2)/(L**2*(1-z_minus)) - - r_plus = 1+sqrt(1-a**2) - r_minus = 1-sqrt(1-a**2) - - h_plus = (r1-r2)*(r3-r_plus)/((r1-r3)*(r2-r_plus)) - h_minus = (r1-r2)*(r3-r_minus)/((r1-r3)*(r2-r_minus)) - - # equation 13 - k_theta = sqrt(z_minus/z_plus) - k_r = sqrt((r1-r2)*(r3-r4)/((r1-r3)*(r2-r4))) - - # equation 21 - return 2*upsilon_theta/(pi*sqrt(e0zp))*_ellippi(z_minus,k_theta) \ - + 2*a*upsilon_r/(pi*(r_plus-r_minus)*sqrt((1-E**2)*(r1-r3)*(r2-r4))) \ - * ((2*E*r_plus-a*L)/(r3-r_plus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_plus)*_ellippi(h_plus,k_r)) \ - - (2*E*r_minus-a*L)/(r3-r_minus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_minus)*_ellippi(h_minus,k_r)) - ) - -def gamma_from_constants(a,constants,radial_roots,polar_roots,upsilon_r=None,upsilon_theta=None): - """ - Computes the average rate at which observer time accumulates in Mino time. - - :param a: dimensionless spin of the black hole - :type a: double - :param constants: dimensionless constants of motion for the orbit in the form :math:`(E,L,Q)` - :type constants: tuple(double, double, double) - :param radial_roots: tuple containing the four roots of the radial polynomial :math:`(r_1, r_2, r_3, r_4)`. - :type radial_roots: tuple(double, double, double, double) - :param polar_roots: tuple containing the roots of the polar equation :math:`(z_-, z_+)` - :type polar_roots: tuple(double, double) - :param upsilon_r: Mino frequency of motion in r can be passed in to speed computation if it is already known - :type upsilon_r: double, optional - :param upsilon_theta: Mino frequency of motion in theta can be passed in to speed computation if it is already known - :type upsilon_theta: double, optional - - :rtype: double - """ - r1,r2,r3,r4 = radial_roots - z_minus, z_plus = polar_roots - - e = (r1-r2)/(r1+r2) - # marginally bound case - if e == 1: - return inf - - E, L, Q = constants - - epsilon0 = a**2*(1-E**2)/L**2 - # simplified form of a**2*sqrt(z_plus/epsilon0) - a2sqrt_zp_over_e0 = L**2/((1-E**2)*sqrt(1-z_minus)) if a == 0 else a**2*z_plus/sqrt(epsilon0*z_plus) - - # compute frequencies if they are not passed in - if upsilon_r is None: upsilon_r = r_frequency_from_constants(a,constants,radial_roots) - if upsilon_theta is None: upsilon_theta = theta_frequency_from_constants(a,constants,radial_roots,polar_roots) - - r_plus = 1+sqrt(1-a**2) - r_minus = 1-sqrt(1-a**2) - - h_r = (r1-r2)/(r1-r3) - h_plus = (r1-r2)*(r3-r_plus)/((r1-r3)*(r2-r_plus)) - h_minus = (r1-r2)*(r3-r_minus)/((r1-r3)*(r2-r_minus)) - - # equation 13 - k_theta = 0 if a == 0 else sqrt(z_minus/z_plus) - k_r = sqrt((r1-r2)*(r3-r4)/((r1-r3)*(r2-r4))) - - # equation 21 - return 4*E + 2*a2sqrt_zp_over_e0*E*upsilon_theta*(ellipk(k_theta**2)-ellipe(k_theta**2))/(pi*L) \ - + 2*upsilon_r/(pi*sqrt((1-E**2)*(r1-r3)*(r2-r4))) \ - * (E/2*((r3*(r1+r2+r3)-r1*r2)*ellipk(k_r**2) + (r2-r3)*(r1+r2+r3+r4)*_ellippi(h_r,k_r)+ (r1-r3)*(r2-r4)*ellipe(k_r**2)) \ - + 2*E*(r3*ellipk(k_r**2)+(r2-r3)*_ellippi(h_r,k_r))\ - +2/(r_plus-r_minus)*(((4*E-a*L)*r_plus-2*a**2*E)/(r3-r_plus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_plus)*_ellippi(h_plus,k_r)) \ - -((4*E-a*L)*r_minus-2*a**2*E)/(r3-r_minus)*(ellipk(k_r**2)-(r2-r3)/(r2-r_minus)*_ellippi(h_minus,k_r)) - ) - ) - -def mino_frequencies_from_constants(a,constants,radial_roots,polar_roots): - r""" - Computes frequencies of orbital motion in Mino time. - - :param a: dimensionless spin of the black hole (must satisfy -1 < a < 1) - :type a: double - :param constants: dimensionless constants of motion for the orbit in the form :math:`(E,L,Q)` - :type constants: tuple(double, double, double) - :param radial_roots: tuple containing the four roots of the radial polynomial :math:`(r_1, r_2, r_3, r_4)`. - :type radial_roots: tuple(double, double, double, double) - :param polar_roots: tuple containing the roots of the polar equation :math:`(z_-, z_+)` - :type polar_roots: tuple(double, double) - - :return: tuple of frequencies in the form :math:`(\Upsilon_r, \Upsilon_\theta, \Upsilon_\phi, \Gamma)` - :rtype: tuple(double, double, double, double) - """ - upsilon_r = r_frequency_from_constants(constants,radial_roots) - upsilon_theta = theta_frequency_from_constants(a,constants,radial_roots,polar_roots) - upsilon_phi = phi_frequency_from_constants(a,constants,radial_roots,polar_roots,upsilon_r,upsilon_theta) - Gamma = gamma_from_constants(a,constants,radial_roots,polar_roots,upsilon_r,upsilon_theta) - - return upsilon_r, abs(upsilon_theta), upsilon_phi, Gamma - -def frequencies_from_constants(a,constants,radial_roots,polar_roots): - r""" - Computes frequencies of orbital motion in Boyer-Lindquist time. - - :param a: dimensionless spin of the black hole (must satisfy -1 < a < 1) - :type a: double - :param constants: dimensionless constants of motion for the orbit in the form :math:`(E,L,Q)` - :type constants: tuple(double, double, double) - :param radial_roots: tuple containing the four roots of the radial polynomial :math:`(r_1, r_2, r_3, r_4)`. - :type radial_roots: tuple(double, double, double, double) - :param polar_roots: tuple containing the roots of the polar equation :math:`(z_-, z_+)` - :type polar_roots: tuple(double, double) - - :return: tuple of frequencies in the form :math:`(\Omega_r, \Omega_\theta, \Omega_\phi)` - :rtype: tuple(double, double, double) - """ - upsilon_r = r_frequency_from_constants(constants,radial_roots) - upsilon_theta = theta_frequency_from_constants(a,constants,radial_roots,polar_roots) - upsilon_phi = phi_frequency_from_constants(a,constants,radial_roots,polar_roots,upsilon_r,upsilon_theta) - Gamma = gamma_from_constants(a,constants,radial_roots,polar_roots,upsilon_r,upsilon_theta) - - return upsilon_r/Gamma, abs(upsilon_theta)/Gamma, upsilon_phi/Gamma \ No newline at end of file diff --git a/kerrgeopy/initial_conditions.py b/kerrgeopy/initial_conditions.py index f55b6a3..cb6c655 100644 --- a/kerrgeopy/initial_conditions.py +++ b/kerrgeopy/initial_conditions.py @@ -1,47 +1,13 @@ """ Module containing functions to compute orbital properties from initial conditions """ - -from .frequencies_from_constants import r_frequency_from_constants, theta_frequency_from_constants -from .plunging_solutions import plunging_radial_roots, plunging_mino_frequencies +from .frequencies import _r_frequency_from_constants, _theta_frequency_from_constants, plunging_mino_frequencies +from .constants import plunging_radial_roots import numpy as np from numpy import sin, cos, sqrt, pi, arcsin, arccos from numpy.polynomial import Polynomial from scipy.special import ellipkinc -def apex_from_constants(a,E,L,Q): - r""" - Computes the orbital parameters :math:`(a,p,e,x)` for a stable bound orbit with the given constants of motion - - :param a: spin parameter - :type a: double - :param E: dimensionless energy - :type E: double - :param L: dimensionless angular momentum - :type L: double - :param Q: dimensionless Carter constant - :type Q: double - - :return: tuple of orbital parameters :math:`(a,p,e,x)` - :rtype: tuple(double,double,double,double) - """ - # radial polynomial written in terms of z = cos^2(theta) - R = Polynomial([-a**2*Q, 2*L**2+2*Q+2*a**2*E**2-4*a*E*L, a**2*E**2-L**2-Q-a**2, 2, E**2-1]) - radial_roots = R.roots() - # numpy returns roots in increasing order - r4, r3, r2, r1 = radial_roots - - # polar polynomial written in terms of z = cos^2(theta) - Z = Polynomial([Q,-(Q+a**2*(1-E**2)+L**2),a**2*(1-E**2)]) - polar_roots = Z.roots() - z_minus, z_plus = polar_roots - - p = 2*r1*r2/(r1+r2) - e = (r1-r2)/(r1+r2) - x = np.sign(L)*sqrt(1-z_minus) - - return a, p, e, x - def constants_from_initial_conditions(a,initial_position,initial_velocity): r""" Computes the constants of motion :math:`(E,L,Q)` of the orbit with the given initial conditions @@ -146,6 +112,10 @@ def stable_orbit_initial_phases(a,initial_position,initial_velocity,constants=No :type constants: tuple(double,double,double), optional :param radial_roots: radial roots :math:`(r_1,r_2,r_3,r_4)` can be passed in to avoid recomputing them :type radial_roots: tuple(double,double,double,double), optional + :param upsilon_r: radial frequency :math:`\Upsilon_r` can be passed in to avoid recomputing it + :type upsilon_r: double, optional + :param upsilon_theta: polar frequency :math:`\Upsilon_\theta` can be passed in to avoid recomputing it + :type upsilon_theta: double, optional :param tol: numerical tolerance used when checking for turning points, defaults to 1e-4 :type tol: double, optional @@ -162,6 +132,7 @@ def stable_orbit_initial_phases(a,initial_position,initial_velocity,constants=No # polar polynomial written in terms of z = cos^2(theta) Z = Polynomial([Q,-(Q+a**2*(1-E**2)+L**2),a**2*(1-E**2)]) polar_roots = Z.roots() + if a == 0: polar_roots = [polar_roots[0], polar_roots[0]] z_minus, z_plus = polar_roots # recompute radial roots if they are not passed in @@ -180,7 +151,7 @@ def stable_orbit_initial_phases(a,initial_position,initial_velocity,constants=No q_r0 = pi else: # compute r frequency - if upsilon_r is None: upsilon_r = r_frequency_from_constants(constants,radial_roots) + if upsilon_r is None: upsilon_r = _r_frequency_from_constants(constants,radial_roots) # Fujita and Hikida equation 13 k_r = sqrt((r1-r2)*(r3-r4)/((r1-r3)*(r2-r4))) @@ -202,7 +173,7 @@ def stable_orbit_initial_phases(a,initial_position,initial_velocity,constants=No q_theta0 = pi else: # compute theta frequency - if upsilon_theta is None: upsilon_theta = theta_frequency_from_constants(a,constants,radial_roots,polar_roots) + if upsilon_theta is None: upsilon_theta = _theta_frequency_from_constants(a,constants,radial_roots,polar_roots) k_theta = sqrt(z_minus/z_plus) y_theta = cos(theta0)/sqrt(z_minus) diff --git a/kerrgeopy/orbit.py b/kerrgeopy/orbit.py index 7eb1b84..44c7b3c 100644 --- a/kerrgeopy/orbit.py +++ b/kerrgeopy/orbit.py @@ -4,7 +4,7 @@ from .spacetime import KerrSpacetime from .initial_conditions import * from .units import * -from .constants import scale_constants +from .constants import scale_constants, apex_from_constants from .frequencies import mino_frequencies, fundamental_frequencies from numpy import sin, cos, sqrt, pi import numpy as np @@ -57,20 +57,10 @@ def __init__(self,a,initial_position,initial_velocity, M=None, mu=None): self.omega_r, self.omega_theta, self.omega_phi = fundamental_frequencies(a,p,e,x) self.initial_phases = stable_orbit_initial_phases(a,initial_position,initial_velocity) else: + if a == 0: raise ValueError("Schwarzschild plunges are not currently supported") self.stable = False self.upsilon_r, self.upsilon_theta = plunging_mino_frequencies(a,E,L,Q) self.initial_phases = plunging_orbit_initial_phases(a,initial_position,initial_velocity) - - def orbital_parameters(self): - """ - Returns the orbital parameters :math:`(a,p,e,x)` of the orbit. Raises a ValueError if the orbit is not stable. - - :return: tuple of orbital parameters :math:`(a,p,e,x)` - :rtype: tuple(double,double,double,double) - """ - if not self.stable: raise ValueError("Orbit is not stable") - - return self.a, self.p, self.e, self.x def trajectory(self,initial_phases=None,distance_units="natural",time_units="natural"): r""" @@ -87,13 +77,11 @@ def trajectory(self,initial_phases=None,distance_units="natural",time_units="nat """ if initial_phases is None: initial_phases = self.initial_phases if self.stable: - from .stable_orbit import StableOrbit - orbit = StableOrbit(self.a,self.p,self.e,self.x,M=self.M,mu=self.mu) - return orbit.trajectory(initial_phases,distance_units,time_units) + from .stable import stable_trajectory + return stable_trajectory(self.a,self.p,self.e,self.x,initial_phases,self.M,distance_units,time_units) else: - from .plunging_orbit import PlungingOrbit - orbit = PlungingOrbit(self.a,self.E,self.L,self.Q,M=self.M,mu=self.mu) - return orbit.trajectory(initial_phases,distance_units,time_units) + from .plunge import plunging_trajectory + return plunging_trajectory(self.a,self.E,self.L,self.Q,initial_phases,self.M,distance_units,time_units) def constants_of_motion(self, units="natural"): """ diff --git a/kerrgeopy/plunge.py b/kerrgeopy/plunge.py new file mode 100644 index 0000000..ce9e12d --- /dev/null +++ b/kerrgeopy/plunge.py @@ -0,0 +1,411 @@ +""" +Module implementing the plunging orbit solutions of `Dyson and van de Meent `_ +""" +from numpy import sqrt, arctan, arctan2, arccos, log, pi +from numpy.polynomial import Polynomial +import numpy as np +from scipy.special import ellipj, ellipeinc, ellipk +from .frequencies import _mino_frequencies_from_constants, _ellippiinc +from .stable import * +from .orbit import Orbit +from .units import * + +class PlungingOrbit(Orbit): + r""" + Class representing a plunging orbit in Kerr spacetime. + + :param a: dimensionaless spin parameter + :type a: double + :param E: dimensionless energy + :type E: double + :param L: dimensionless angular momentum + :type L: double + :param Q: dimensionless Carter constant + :type Q: double + :param initial_phases: tuple of initial phases :math:`(q^t_0, q^r_0, q^\theta_0, q^\phi_0)`, defaults to (0,0,0,0) + :type initial_phases: tuple, optional + :param M: mass of the primary in solar masses, optional + :type M: double + :param mu: mass of the smaller body in solar masses, optional + :type mu: double + + :ivar a: dimensionaless spin parameter + :ivar E: dimensionless energy + :ivar L: dimensionless angular momentum + :ivar Q: dimensionless Carter constant + :ivar initial_position: tuple of initial position coordinates :math:`(t_0, r_0, \theta_0, \phi_0)` + :ivar initial_velocity: tuple of initial four-velocity components :math:`(u^t_0, u^r_0, u^\theta_0, u^\phi_0)` + :ivar M: mass of the primary in solar masses + :ivar mu: mass of the smaller body in solar masses + :ivar upsilon_r: radial Mino frequency + :ivar upsilon_theta: polar Mino frequency + """ + def __init__(self,a,E,L,Q,initial_phases=(0,0,0,0),M=None,mu=None): + self.a, self.E, self.L, self.Q, self.initial_phases, self.M, self.mu = a, E, L, Q, initial_phases, M, mu + + if a == 0: raise ValueError("Schwarzschild plunges are not currently supported") + + self.stable = False + self.upsilon_r, self.upsilon_theta = plunging_mino_frequencies(a,E,L,Q) + + u_t, u_r, u_theta, u_phi = self.four_velocity() + t, r, theta, phi = self.trajectory() + self.initial_position = t(0), r(0), theta(0), phi(0) + self.initial_velocity = u_t(0), u_r(0), u_theta(0), u_phi(0) + + def trajectory(self,initial_phases=None,distance_units="natural",time_units="natural"): + r""" + Computes the components of the trajectory as a function of mino time + + :param initial_phases: tuple of initial phases :math:`(q^t_0, q^r_0, q^\theta_0, q^\phi_0)` + :type initial_phases: tuple, optional + :param distance_units: units to compute the radial component of the trajectory in (options are "natural", "mks", "cgs", "au" and "km"), defaults to "natural" + :type distance_units: str, optional + :param time_units: units to compute the time component of the trajectory in (options are "natural", "mks", "cgs", and "days"), defaults to "natural" + :type time_units: str, optional + + :return: tuple of functions :math:`t(\lambda), r(\lambda), \theta(\lambda), \phi(\lambda)` + :rtype: tuple(function,function,function,function) + """ + if initial_phases is None: initial_phases = self.initial_phases + return plunging_trajectory(self.a,self.E,self.L,self.Q,initial_phases,self.M,distance_units,time_units) + +def plunging_radial_integrals(a,E,L,Q): + r""" + Computes the radial integrals :math:`I_r`, :math:`I_{r^2}` and :math:`I_{r_\pm}` defined in equation 39 of `Dyson and van de Meent `_ as a function of the radial phase. + Used to compute the radial solutions for the case of two complex roots. + + :param a: dimensionless spin parameter + :type a: double + :param E: dimensionless energy + :type E: double + :param L: dimensionless angular momentum + :type L: double + :param Q: dimensionless Carter constant + :type Q: double + + :return: radial integrals :math:`(I_r,I_{r^2},I_{r_+},I_{r_-})` + :rtype: tuple(function,function,function,function) + """ + r1, r2, r3, r4 = plunging_radial_roots(a,E,L,Q) + rho_r = np.real(r3) + rho_i = np.imag(r4) + + # inner and outer horizons + r_plus = 1+sqrt(1-a**2) + r_minus = 1-sqrt(1-a**2) + + # equation 42 + A = sqrt((r2-rho_r)**2+rho_i**2) + B = sqrt((r1-rho_r)**2+rho_i**2) + f = 4*A*B/(A-B)**2 + k_r = sqrt(((r2-r1)**2-(A-B)**2)/(4*A*B)) + + upsilon_r, upsilon_theta = plunging_mino_frequencies(a,E,L,Q) + + # equation 48 + D_plus = sqrt(4*A*B*(r2-r_plus)*(r_plus-r1))/(A*(r_plus-r1)+B*(r2-r_plus)) + D_minus = sqrt(4*A*B*(r2-r_minus)*(r_minus-r1))/(A*(r_minus-r1)+B*(r2-r_minus)) + + def I_r(q_r): + # equation 43 + sn, cn, dn, xi_r = ellipj(2*ellipk(k_r**2)*q_r/pi,k_r**2) + mino_time = q_r/upsilon_r + # equation 46 + return (A*r1-B*r2)/(A-B)*mino_time - 1/sqrt(1-E**2)*arctan((r2-r1)*sn/(2*sqrt(A*B)*dn)) \ + + (A+B)*(r2-r1)/(2*(A-B)*sqrt(A*B*(1-E**2)))*_ellippiinc(xi_r,-1/f,k_r**2) + + def I_r2(q_r): + # equation 43 + sn, cn, dn, xi_r = ellipj(2*ellipk(k_r**2)*q_r/pi,k_r**2) + mino_time = q_r/upsilon_r + # equation 47 + return (A*r1**2-B*r2**2)/(A-B)*mino_time + sqrt(A*B)/sqrt(1-E**2)*ellipeinc(xi_r,k_r**2) \ + - (A+B)*(A**2+2*r1**2-B**2-2*r2**2)/(4*(A-B)*sqrt((1-E**2)*A*B))*_ellippiinc(xi_r,-1/f,k_r**2) \ + - sqrt(A*B)*(A+B-(A-B)*cn)*sn*dn/((A-B)*sqrt(1-E**2)*(f+sn**2)) \ + + (A**2+2*r1**2-B**2-2*r2**2)/(4*(r2-r1)*sqrt(1-E**2))*arctan2(2*sn*dn*sqrt(f*(1+f*k_r**2)),f-(1+2*f*k_r**2)*sn**2) + + def I_r_plus(q_r): + # equation 43 + sn, cn, dn, xi_r = ellipj(2*ellipk(k_r**2)*q_r/pi,k_r**2) + mino_time = q_r/upsilon_r + # equation 48 + return (A-B)*mino_time/(A*(r1-r_plus)-B*(r2-r_plus)) \ + + (r2-r1)*(A*(r1-r_plus)+B*(r2-r_plus))/(2*sqrt(A*B*(1-E**2))*(r_plus-r1)*(r2-r_plus)*(A*(r1-r_plus)-B*(r2-r_plus)))*_ellippiinc(xi_r,1/D_plus**2,k_r**2) \ + - (sqrt(r2-r1)*log( + ((D_plus*sqrt(1-D_plus**2*k_r**2)+dn*sn)**2 + (k_r*(D_plus**2-sn**2))**2) / + ((D_plus*sqrt(1-D_plus**2*k_r**2)-dn*sn)**2 + (k_r*(D_plus**2-sn**2))**2) + ) / + (4*sqrt((1-E**2)*(r2-r_plus)*(r_plus-r1))*sqrt(A**2*(r_plus-r1)-(r2-r_plus)*(r1**2-B**2+r2*r_plus-r1*(r2+r_plus)))) + ) + + def I_r_minus(q_r): + # equation 43 + sn, cn, dn, xi_r = ellipj(2*ellipk(k_r**2)*q_r/pi,k_r**2) + mino_time = q_r/upsilon_r + # equation 48 + return (A-B)*mino_time/(A*(r1-r_minus)-B*(r2-r_minus)) \ + + (r2-r1)*(A*(r1-r_minus)+B*(r2-r_minus))/(2*sqrt(A*B*(1-E**2))*(r_minus-r1)*(r2-r_minus)*(A*(r1-r_minus)-B*(r2-r_minus)))*_ellippiinc(xi_r,1/D_minus**2,k_r**2) \ + - (sqrt(r2-r1)*log( + ((D_minus*sqrt(1-D_minus**2*k_r**2)+dn*sn)**2 + (k_r*(D_minus**2-sn**2))**2) / + ((D_minus*sqrt(1-D_minus**2*k_r**2)-dn*sn)**2 + (k_r*(D_minus**2-sn**2))**2) + ) / + (4*sqrt((1-E**2)*(r2-r_minus)*(r_minus-r1))*sqrt(A**2*(r_minus-r1)-(r2-r_minus)*(r1**2-B**2+r2*r_minus-r1*(r2+r_minus)))) + ) + + return I_r, I_r2, I_r_plus, I_r_minus + +def plunging_radial_solutions_complex(a,E,L,Q): + r""" + Computes the radial solutions :math:`r(q_r), t_r(q_r), \phi_r(q_r)` from equation 50 and 51 of `Dyson and van de Meent `_ for the case of two complex radial roots. + + :param a: dimensionless spin parameter + :type a: double + :param E: dimensionless energy + :type E: double + :param L: dimensionless angular momentum + :type L: double + :param Q: dimensionless Carter constant + :type Q: double + + :return: tuple of radial solutions :math:`(r(q_r), t_r(q_r), \phi_r(q_r))` + :rtype: tuple(function, function, function) + """ + roots = plunging_radial_roots(a,E,L,Q) + if np.iscomplex(roots).sum() != 2: raise ValueError("There should be two complex roots") + + r1, r2, r3, r4 = roots # r1 < r2 are real and r3/r4 are complex conjugates + rho_r = np.real(r3) + rho_i = np.imag(r4) + + # inner and outer horizons + r_plus = 1+sqrt(1-a**2) + r_minus = 1-sqrt(1-a**2) + + # equation 42 + A = sqrt((r2-rho_r)**2+rho_i**2) + B = sqrt((r1-rho_r)**2+rho_i**2) + k_r = sqrt(((r2-r1)**2-(A-B)**2)/(4*A*B)) + + upsilon_r, upsilon_theta = plunging_mino_frequencies(a,E,L,Q) + + I_r, I_r2, I_r_plus, I_r_minus = plunging_radial_integrals(a,E,L,Q) + + def r(q_r): + # equation 43 + sn, cn, dn, xi_r = ellipj(2*ellipk(k_r**2)*q_r/pi,k_r**2) + # equation 49 + return ((A-B)*(A*r1-B*r2)*sn**2+2*A*B*(r1+r2)-2*A*B*(r2-r1)*cn)/(4*A*B+(A-B)**2*sn**2) + + def t_r(q_r): + mino_time = q_r/upsilon_r + # equation 41 + return E*(r_plus**2+r_minus**2+r_plus*r_minus+2*a**2)*mino_time + \ + E*(I_r2(q_r)+I_r(q_r)*(r_minus+r_plus)) + \ + ((r_minus**2+a**2)*(E*(r_minus**2+a**2)-a*L)/(r_minus-r_plus)*I_r_minus(q_r) + + (r_plus**2+a**2)*(E*(r_plus**2+a**2)-a*L)/(r_plus-r_minus)*I_r_plus(q_r) + ) + + def phi_r(q_r): + mino_time = q_r/upsilon_r + # equation 40 + return a*( + (E*(r_minus**2+a**2)-a*L)/(r_minus-r_plus)*I_r_minus(q_r) + + (E*(r_plus**2+a**2)-a*L)/(r_plus-r_minus)*I_r_plus(q_r) + ) + + return r, t_r, phi_r + +def plunging_polar_solutions(a,E,L,Q): + r""" + Computes the polar solutions :math:`\theta(q_\theta), t_\theta(q_\theta), \phi_\theta(q_\theta)` from equation 33 and 37 of `Dyson and van de Meent `_ + + :param a: dimensionless spin parameter + :type a: double + :param E: dimensionless energy + :type E: double + :param L: dimensionless angular momentum + :type L: double + :param Q: dimensionless Carter constant + :type Q: double + + :return: tuple of polar solutions :math:`(\theta(q_\theta), t_\theta(q_\theta), \phi_\theta(q_\theta))` + :rtype: tuple(function, function, function) + """ + z1 = sqrt(1/2*(1+(L**2+Q)/(a**2*(1-E**2))-sqrt((1+(L**2+Q)/(a**2*(1-E**2)))**2-4*Q/(a**2*(1-E**2))))) + z2 = sqrt(a**2*(1-E**2)/2*(1+(L**2+Q)/(a**2*(1-E**2))+sqrt((1+(L**2+Q)/(a**2*(1-E**2)))**2-4*Q/(a**2*(1-E**2))))) + k_theta = a*sqrt(1-E**2)*z1/z2 + + upsilon_r, upsilon_theta = plunging_mino_frequencies(a,E,L,Q) + + def theta(q_theta): + mino_time = q_theta/upsilon_theta + # equation 28 + sn, cn, dn, xi_theta = ellipj(z2*mino_time,k_theta**2) + # equation 27 + return arccos(z1*sn) + + def t_theta(q_theta): + mino_time = q_theta/upsilon_theta + # equation 28 + sn, cn, dn, xi_theta = ellipj(z2*mino_time,k_theta**2) + # equation 35 + return E/(1-E**2)*((z2**2-a**2*(1-E**2))*mino_time-z2*ellipeinc(xi_theta,k_theta**2)) + + def phi_theta(q_theta): + mino_time = q_theta/upsilon_theta + # equation 28 + sn, cn, dn, xi_theta = ellipj(z2*mino_time,k_theta**2) + # equation 31 + return L/z2*_ellippiinc(xi_theta,z1**2,k_theta**2) + + return theta, t_theta, phi_theta + +def plunging_trajectory(a,E,L,Q,initial_phases=(0,0,0,0),M=None,distance_units="natural",time_units="natural"): + r""" + Computes the components of the trajectory as a function of mino time for a plunging orbit with the given parameters. + + :param a: dimensionless spin parameter + :type a: double + :param E: dimensionless energy + :type E: double + :param L: dimensionless angular momentum + :type L: double + :param Q: dimensionless Carter constant + :type Q: double + :param initial_phases: tuple of initial phases :math:`(q^t_0, q^r_0, q^\theta_0, q^\phi_0)` + :type initial_phases: tuple, optional + :param M: mass of the primary in solar masses + :type M: double, optional + :param distance_units: units to compute the radial component of the trajectory in (options are "natural", "mks", "cgs", "au" and "km"), defaults to "natural" + :type distance_units: str, optional + :param time_units: units to compute the time component of the trajectory in (options are "natural", "mks", "cgs", and "days"), defaults to "natural" + :type time_units: str, optional + + :return: components of the trajectory :math:`t(\lambda), r(\lambda), \theta(\lambda), \phi(\lambda)` + :rtype: tuple(function,function,function,function) + """ + radial_roots = plunging_radial_roots(a,E,L,Q) + + if np.iscomplex(radial_roots[3]): + return _complex_plunge_trajectory(a,E,L,Q,initial_phases,M,distance_units,time_units) + else: + return _real_plunge_trajectory(a,E,L,Q,initial_phases,M,distance_units,time_units) + +def _complex_plunge_trajectory(a,E,L,Q,initial_phases=(0,0,0,0), M=None, distance_units="natural",time_units="natural"): + r""" + Computes the components of the trajectory for a plunging orbit with the given parameters assuming that the radial polynomial has two complex roots and two real roots. + + :param a: dimensionless spin parameter + :type a: double + :param E: dimensionless energy + :type E: double + :param L: dimensionless angular momentum + :type L: double + :param Q: dimensionless Carter constant + :type Q: double + :param initial_phases: tuple of initial phases :math:`(q^t_0, q^r_0, q^\theta_0, q^\phi_0)` + :type initial_phases: tuple, optional + :param M: mass of the primary in solar masses + :type M: double, optional + :param distance_units: units to compute the radial component of the trajectory in (options are "natural", "mks", "cgs", "au" and "km"), defaults to "natural" + :type distance_units: str, optional + :param time_units: units to compute the time component of the trajectory in (options are "natural", "mks", "cgs", and "days"), defaults to "natural" + :type time_units: str, optional + + :return: components of the trajectory :math:`t(\lambda), r(\lambda), \theta(\lambda), \phi(\lambda)` + :rtype: tuple(function,function,function,function) + """ + if ((distance_units != "natural") or (time_units != "natural")) and M is None: raise ValueError("M must be specified to convert to physical units") + + upsilon_r, upsilon_theta = plunging_mino_frequencies(a,E,L,Q) + r_phases, t_r, phi_r = plunging_radial_solutions_complex(a,E,L,Q) + theta_phases, t_theta, phi_theta = plunging_polar_solutions(a,E,L,Q) + q_t0, q_r0, q_theta0, q_phi0 = initial_phases + + distance_conversion_func = {"natural": lambda x,M: x, "mks": distance_in_meters, "cgs": distance_in_cm, "au": distance_in_au,"km": distance_in_km} + time_conversion_func = {"natural": lambda x,M: x, "mks": time_in_seconds, "cgs": time_in_seconds, "days": time_in_days} + + # adjust q_theta0 so that initial conditions are consistent with stable orbits + q_theta0 = q_theta0 + pi/2 + + # Calculate normalization constants so that t = 0 and phi = 0 at lambda = 0 when q_t0 = 0 and q_phi0 = 0 + C_t = t_r(q_r0)+t_theta(q_theta0) + C_phi= phi_r(q_r0)+phi_theta(q_theta0) + + def t(mino_time): + return time_conversion_func[time_units](t_r(upsilon_r*mino_time+q_r0) + t_theta(upsilon_theta*mino_time+q_theta0) - C_t + q_t0, M) + + def r(mino_time): + return distance_conversion_func[distance_units](r_phases(upsilon_r*mino_time+q_r0), M) + + def theta(mino_time): + return theta_phases(upsilon_theta*mino_time+q_theta0) + + def phi(mino_time): + return phi_r(upsilon_r*mino_time+q_r0) + phi_theta(upsilon_theta*mino_time+q_theta0) - C_phi + q_phi0 + + return t, r, theta, phi + + +def _real_plunge_trajectory(a,E,L,Q,initial_phases=(0,0,0,0),M=None,distance_units="natural",time_units="natural"): + r""" + Computes the components of the trajectory for a plunging orbit with the given parameters assuming that there are four real roots of the radial polynomial. + + :param a: dimensionless spin parameter + :type a: double + :param E: dimensionless energy + :type E: double + :param L: dimensionless angular momentum + :type L: double + :param Q: dimensionless Carter constant + :type Q: double + :param initial_phases: tuple of initial phases :math:`(q^t_0, q^r_0, q^\theta_0, q^\phi_0)` + :type initial_phases: tuple, optional + :param M: mass of the primary in solar masses + :type M: double, optional + :param distance_units: units to compute the radial component of the trajectory in (options are "natural", "mks", "cgs", "au" and "km"), defaults to "natural" + :type distance_units: str, optional + :param time_units: units to compute the time component of the trajectory in (options are "natural", "mks", "cgs", and "days"), defaults to "natural" + :type time_units: str, optional + + :return: components of the trajectory :math:`t(\lambda), r(\lambda), \theta(\lambda), \phi(\lambda)` + :rtype: tuple(function,function,function,function) + """ + if ((distance_units != "natural") or (time_units != "natural")) and M is None: raise ValueError("M must be specified to convert to physical units") + + constants = (E,L,Q) + radial_roots = plunging_radial_roots(a,E,L,Q) + + # polar polynomial written in terms of z = cos^2(theta) + Z = Polynomial([Q,-(Q+a**2*(1-E**2)+L**2),a**2*(1-E**2)]) + polar_roots = Z.roots() + if a == 0: polar_roots = [polar_roots[0], polar_roots[0]] + + upsilon_r, upsilon_theta, upsilon_phi, gamma = _mino_frequencies_from_constants(a,constants,radial_roots,polar_roots) + r_phases, t_r, phi_r = radial_solutions(a,constants,radial_roots) + theta_phases, t_theta, phi_theta = polar_solutions(a,constants,polar_roots) + q_t0, q_r0, q_theta0, q_phi0 = initial_phases + + distance_conversion_func = {"natural": lambda x,M: x, "mks": distance_in_meters, "cgs": distance_in_cm, "au": distance_in_au,"km": distance_in_km} + time_conversion_func = {"natural": lambda x,M: x, "mks": time_in_seconds, "cgs": time_in_seconds, "days": time_in_days} + + # Calculate normalization constants so that t = 0 and phi = 0 at lambda = 0 when q_t0 = 0 and q_phi0 = 0 + C_t = t_r(q_r0)+t_theta(q_theta0) + C_phi= phi_r(q_r0)+phi_theta(q_theta0) + + def t(mino_time): + # equation 6 + return time_conversion_func[time_units](q_t0 + gamma*mino_time + t_r(upsilon_r*mino_time+q_r0) + t_theta(upsilon_theta*mino_time+q_theta0) - C_t, M) + + def r(mino_time): + return distance_conversion_func[distance_units](r_phases(upsilon_r*mino_time+q_r0), M) + + def theta(mino_time): + return theta_phases(upsilon_theta*mino_time+q_theta0) + + def phi(mino_time): + # equation 6 + return q_phi0 + upsilon_phi*mino_time + phi_r(upsilon_r*mino_time+q_r0) + phi_theta(upsilon_theta*mino_time+q_theta0) - C_phi + + return t, r, theta, phi \ No newline at end of file diff --git a/kerrgeopy/plunging_orbit.py b/kerrgeopy/plunging_orbit.py deleted file mode 100644 index ad699d4..0000000 --- a/kerrgeopy/plunging_orbit.py +++ /dev/null @@ -1,160 +0,0 @@ -""" -Module containing the PlungingOrbit class -""" -from .plunging_solutions import * -from .plunging_solutions import plunging_radial_roots -from .stable_solutions import * -from .orbit import Orbit - -class PlungingOrbit(Orbit): - r""" - Class representing a plunging orbit in Kerr spacetime. - - :param a: dimensionaless spin parameter - :type a: double - :param E: dimensionless energy - :type E: double - :param L: dimensionless angular momentum - :type L: double - :param Q: dimensionless Carter constant - :type Q: double - :param initial_phases: tuple of initial phases :math:`(q^t_0, q^r_0, q^\theta_0, q^\phi_0)`, defaults to (0,0,0,0) - :type initial_phases: tuple, optional - :param M: mass of the primary in solar masses, optional - :type M: double - :param mu: mass of the smaller body in solar masses, optional - :type mu: double - - :ivar a: dimensionaless spin parameter - :ivar E: dimensionless energy - :ivar L: dimensionless angular momentum - :ivar Q: dimensionless Carter constant - :ivar initial_position: tuple of initial position coordinates :math:`(t_0, r_0, \theta_0, \phi_0)` - :ivar initial_velocity: tuple of initial four-velocity components :math:`(u^t_0, u^r_0, u^\theta_0, u^\phi_0)` - :ivar M: mass of the primary in solar masses - :ivar mu: mass of the smaller body in solar masses - :ivar upsilon_r: radial Mino frequency - :ivar upsilon_theta: polar Mino frequency - """ - def __init__(self,a,E,L,Q,initial_phases=(0,0,0,0),M=None,mu=None): - self.a, self.E, self.L, self.Q, self.initial_phases, self.M, self.mu = a, E, L, Q, initial_phases, M, mu - self.stable = False - self.upsilon_r, self.upsilon_theta = plunging_mino_frequencies(a,E,L,Q) - - u_t, u_r, u_theta, u_phi = self.four_velocity() - t, r, theta, phi = self.trajectory() - self.initial_position = t(0), r(0), theta(0), phi(0) - self.initial_velocity = u_t(0), u_r(0), u_theta(0), u_phi(0) - - def trajectory(self,initial_phases=None,distance_units="natural",time_units="natural"): - r""" - Computes the components of the trajectory as a function of mino time - - :param initial_phases: tuple of initial phases :math:`(q^t_0, q^r_0, q^\theta_0, q^\phi_0)` - :type initial_phases: tuple, optional - :param distance_units: units to compute the radial component of the trajectory in (options are "natural", "mks", "cgs", "au" and "km"), defaults to "natural" - :type distance_units: str, optional - :param time_units: units to compute the time component of the trajectory in (options are "natural", "mks", "cgs", and "days"), defaults to "natural" - :type time_units: str, optional - - :return: tuple of functions :math:`t(\lambda), r(\lambda), \theta(\lambda), \phi(\lambda)` - :rtype: tuple(function,function,function,function) - """ - - if distance_units != "natural" and (self.M is None or self.mu is None): raise ValueError("M and mu must be specified to convert r to physical units") - if time_units != "natural" and (self.M is None or self.mu is None): raise ValueError("M and mu must be specified to convert t to physical units") - - a, E, L, Q = self.a, self.E, self.L, self.Q - radial_roots = plunging_radial_roots(a,E,L,Q) - - if np.iscomplex(radial_roots[3]): - return self._complex_trajectory(initial_phases,distance_units,time_units) - else: - return self._real_trajectory(initial_phases,distance_units,time_units) - - def _complex_trajectory(self,initial_phases=None, distance_units="natural",time_units="natural"): - r""" - Computes the components of the trajectory in the case of two complex and two real radial roots - - :param initial_phases: tuple of initial phases :math:`(q^t_0, q^r_0, q^\theta_0, q^\phi_0)`, defaults to (0,0,0,0) - :type initial_phases: tuple, optional - - :rtype: tuple(function,function,function,function) - """ - if initial_phases is None: initial_phases = self.initial_phases - a, E, L, Q = self.a, self.E, self.L, self.Q - upsilon_r, upsilon_theta = self.upsilon_r, self.upsilon_theta - r_phases, t_r, phi_r = plunging_radial_solutions_complex(a,E,L,Q) - theta_phases, t_theta, phi_theta = plunging_polar_solutions(a,E,L,Q) - q_t0, q_r0, q_theta0, q_phi0 = initial_phases - - distance_conversion_func = {"natural": lambda x,M: x, "mks": distance_in_meters, "cgs": distance_in_cm, "au": distance_in_au,"km": distance_in_km} - time_conversion_func = {"natural": lambda x,M: x, "mks": time_in_seconds, "cgs": time_in_seconds, "days": time_in_days} - - # adjust q_theta0 so that initial conditions are consistent with stable case - q_theta0 = q_theta0 + pi/2 - - # Calculate normalization constants so that t = 0 and phi = 0 at lambda = 0 when q_t0 = 0 and q_phi0 = 0 - C_t = t_r(q_r0)+t_theta(q_theta0) - C_phi= phi_r(q_r0)+phi_theta(q_theta0) - - - def t(mino_time): - return time_conversion_func[time_units](t_r(upsilon_r*mino_time+q_r0) + t_theta(upsilon_theta*mino_time+q_theta0) - C_t + q_t0, self.M) - - def r(mino_time): - return distance_conversion_func[distance_units](r_phases(upsilon_r*mino_time+q_r0), self.M) - - def theta(mino_time): - return theta_phases(upsilon_theta*mino_time+q_theta0) - - def phi(mino_time): - return phi_r(upsilon_r*mino_time+q_r0) + phi_theta(upsilon_theta*mino_time+q_theta0) - C_phi + q_phi0 - - return t, r, theta, phi - - def _real_trajectory(self,initial_phases=None,distance_units="natural",time_units="natural"): - r""" - Computes the components of the trajectory in the case of four real radial roots - - :param initial_phases: tuple of initial phases :math:`(q^t_0, q^r_0, q^\theta_0, q^\phi_0)` - :type initial_phases: tuple, optional - - :rtype: tuple(function,function,function,function) - """ - if initial_phases is None: initial_phases = self.initial_phases - a, E, L, Q = self.a, self.E, self.L, self.Q - constants = (E,L,Q) - radial_roots = plunging_radial_roots(a,E,L,Q) - - # polar polynomial written in terms of z = cos^2(theta) - Z = Polynomial([Q,-(Q+a**2*(1-E**2)+L**2),a**2*(1-E**2)]) - polar_roots = Z.roots() - - upsilon_r, upsilon_theta, upsilon_phi, gamma = mino_frequencies_from_constants(a,constants,radial_roots,polar_roots) - r_phases, t_r, phi_r = radial_solutions(a,constants,radial_roots) - theta_phases, t_theta, phi_theta = polar_solutions(a,constants,polar_roots) - q_t0, q_r0, q_theta0, q_phi0 = initial_phases - - distance_conversion_func = {"natural": lambda x,M: x, "mks": distance_in_meters, "cgs": distance_in_cm, "au": distance_in_au,"km": distance_in_km} - time_conversion_func = {"natural": lambda x,M: x, "mks": time_in_seconds, "cgs": time_in_seconds, "days": time_in_days} - - # Calculate normalization constants so that t = 0 and phi = 0 at lambda = 0 when q_t0 = 0 and q_phi0 = 0 - C_t = t_r(q_r0)+t_theta(q_theta0) - C_phi= phi_r(q_r0)+phi_theta(q_theta0) - - def t(mino_time): - # equation 6 - return time_conversion_func[time_units](q_t0 + gamma*mino_time + t_r(upsilon_r*mino_time+q_r0) + t_theta(upsilon_theta*mino_time+q_theta0) - C_t, self.M) - - def r(mino_time): - return distance_conversion_func[distance_units](r_phases(upsilon_r*mino_time+q_r0), self.M) - - def theta(mino_time): - return theta_phases(upsilon_theta*mino_time+q_theta0) - - def phi(mino_time): - # equation 6 - return q_phi0 + upsilon_phi*mino_time + phi_r(upsilon_r*mino_time+q_r0) + phi_theta(upsilon_theta*mino_time+q_theta0) - C_phi - - return t, r, theta, phi \ No newline at end of file diff --git a/kerrgeopy/plunging_solutions.py b/kerrgeopy/plunging_solutions.py deleted file mode 100644 index a7483dd..0000000 --- a/kerrgeopy/plunging_solutions.py +++ /dev/null @@ -1,292 +0,0 @@ -""" -Module containing functions to compute intermediate terms in the plunging orbit solutions of `Dyson and van de Meent `_ -""" -from numpy import sqrt, arctan, arctan2, arccos, log, pi -from numpy.polynomial import Polynomial -import numpy as np -from scipy.special import ellipj, ellipeinc, ellipk -from .stable_solutions import _ellippiinc -from .frequencies_from_constants import mino_frequencies_from_constants - -def plunging_radial_roots(a,E,L,Q): - """ - Computes the radial roots for a plunging orbit. - If all roots are real, roots are sorted such that the motion is between r1 and r2 and roots are otherwise in decreasing order. - If there are two complex roots, r1 < r2 are real and r3/r4 are complex conjugates. - - :param a: dimensionless spin parameter - :type a: double - :param E: dimensionless energy - :type E: double - :param L: dimensionless angular momentum - :type L: double - :param Q: dimensionless Carter constant - :type Q: double - - :return: tuple of radial roots - :rtype: tuple(double,double,double,double) - """ - # standard form of the radial polynomial R(r) - R = Polynomial([-a**2*Q, 2*L**2+2*Q+2*a**2*E**2-4*a*E*L, a**2*E**2-L**2-Q-a**2, 2, E**2-1]) - roots = R.roots() - # get the real roots and the complex roots - real_roots = np.sort(np.real(roots[np.isreal(roots)])) - complex_roots = roots[np.iscomplex(roots)] - - r_minus = 1-sqrt(1-a**2) - - # if there are 4 real roots, by convention r4 < r3 < r2 < r1 (consistent with stable orbits) - if len(real_roots) == 4: - # if there are three roots outside the event horizon swap r1/r3 and r2/r4 - if real_roots[1] > r_minus: - r1 = real_roots[1] - r2 = real_roots[0] - r3 = real_roots[3] - r4 = real_roots[2] - else: - r4, r3, r2, r1 = real_roots - - # in the case of two complex roots, r1 < r2 are real and r3/r4 are complex conjugates - elif len(real_roots) == 2: - r1, r2 = real_roots - r3, r4 = complex_roots - - return r1, r2, r3, r4 - -def plunging_mino_frequencies(a,E,L,Q): - r""" - Computes the radial and polar mino frequencies for a plunging orbit. - - :param a: dimensionless spin parameter - :type a: double - :param E: dimensionless energy - :type E: double - :param L: dimensionless angular momentum - :type L: double - :param Q: dimensionless Carter constant - :type Q: double - - :return: radial and polar mino frequencies :math:`(\Upsilon_r,\Upsilon_\theta)` - :rtype: tuple(double,double) - """ - - radial_roots = plunging_radial_roots(a,E,L,Q) - r1, r2, r3, r4 = radial_roots - rho_r = np.real(r3) - rho_i = np.imag(r4) - if np.iscomplex(radial_roots[3]): - # equation 42 - A = sqrt((r2-rho_r)**2+rho_i**2) - B = sqrt((r1-rho_r)**2+rho_i**2) - k_r = sqrt(abs(((r2-r1)**2-(A-B)**2)/(4*A*B))) - - # equation 52 - upsilon_r = pi*sqrt(A*B*(1-E**2))/(2*ellipk(k_r**2)) - - # equation 14 - z1 = sqrt(1/2*(1+(L**2+Q)/(a**2*(1-E**2))-sqrt((1+(L**2+Q)/(a**2*(1-E**2)))**2-4*Q/(a**2*(1-E**2))))) - z2 = sqrt(a**2*(1-E**2)/2*(1+(L**2+Q)/(a**2*(1-E**2))+sqrt((1+(L**2+Q)/(a**2*(1-E**2)))**2-4*Q/(a**2*(1-E**2))))) - k_theta = a*sqrt(1-E**2)*z1/z2 - - # equation 53 - upsilon_theta = pi/2*z2/ellipk(k_theta**2) - else: - # polar polynomial written in terms of z = cos^2(theta) - Z = Polynomial([Q,-(Q+a**2*(1-E**2)+L**2),a**2*(1-E**2)]) - polar_roots = Z.roots() - constants = (E,L,Q) - upsilon_r, upsilon_theta, upsilon_phi, gamma = mino_frequencies_from_constants(a,constants,radial_roots,polar_roots) - - - return upsilon_r, upsilon_theta - -def plunging_radial_integrals(a,E,L,Q): - r""" - Computes the radial integrals :math:`I_r`, :math:`I_{r^2}` and :math:`I_{r_\pm}` defined in equation 39 of `Dyson and van de Meent `_ as a function of the radial phase. - Used to compute the radial solutions for the case of two complex roots. - - :param a: dimensionless spin parameter - :type a: double - :param E: dimensionless energy - :type E: double - :param L: dimensionless angular momentum - :type L: double - :param Q: dimensionless Carter constant - :type Q: double - - :return: radial integrals :math:`(I_r,I_{r^2},I_{r_+},I_{r_-})` - :rtype: tuple(function,function,function,function) - """ - r1, r2, r3, r4 = plunging_radial_roots(a,E,L,Q) - rho_r = np.real(r3) - rho_i = np.imag(r4) - - # inner and outer horizons - r_plus = 1+sqrt(1-a**2) - r_minus = 1-sqrt(1-a**2) - - # equation 42 - A = sqrt((r2-rho_r)**2+rho_i**2) - B = sqrt((r1-rho_r)**2+rho_i**2) - f = 4*A*B/(A-B)**2 - k_r = sqrt(((r2-r1)**2-(A-B)**2)/(4*A*B)) - - upsilon_r, upsilon_theta = plunging_mino_frequencies(a,E,L,Q) - - # equation 48 - D_plus = sqrt(4*A*B*(r2-r_plus)*(r_plus-r1))/(A*(r_plus-r1)+B*(r2-r_plus)) - D_minus = sqrt(4*A*B*(r2-r_minus)*(r_minus-r1))/(A*(r_minus-r1)+B*(r2-r_minus)) - - def I_r(q_r): - # equation 43 - sn, cn, dn, xi_r = ellipj(2*ellipk(k_r**2)*q_r/pi,k_r**2) - mino_time = q_r/upsilon_r - # equation 46 - return (A*r1-B*r2)/(A-B)*mino_time - 1/sqrt(1-E**2)*arctan((r2-r1)*sn/(2*sqrt(A*B)*dn)) \ - + (A+B)*(r2-r1)/(2*(A-B)*sqrt(A*B*(1-E**2)))*_ellippiinc(xi_r,-1/f,k_r) - - def I_r2(q_r): - # equation 43 - sn, cn, dn, xi_r = ellipj(2*ellipk(k_r**2)*q_r/pi,k_r**2) - mino_time = q_r/upsilon_r - # equation 47 - return (A*r1**2-B*r2**2)/(A-B)*mino_time + sqrt(A*B)/sqrt(1-E**2)*ellipeinc(xi_r,k_r**2) \ - - (A+B)*(A**2+2*r1**2-B**2-2*r2**2)/(4*(A-B)*sqrt((1-E**2)*A*B))*_ellippiinc(xi_r,-1/f,k_r) \ - - sqrt(A*B)*(A+B-(A-B)*cn)*sn*dn/((A-B)*sqrt(1-E**2)*(f+sn**2)) \ - + (A**2+2*r1**2-B**2-2*r2**2)/(4*(r2-r1)*sqrt(1-E**2))*arctan2(2*sn*dn*sqrt(f*(1+f*k_r**2)),f-(1+2*f*k_r**2)*sn**2) - - def I_r_plus(q_r): - # equation 43 - sn, cn, dn, xi_r = ellipj(2*ellipk(k_r**2)*q_r/pi,k_r**2) - mino_time = q_r/upsilon_r - # equation 48 - return (A-B)*mino_time/(A*(r1-r_plus)-B*(r2-r_plus)) \ - + (r2-r1)*(A*(r1-r_plus)+B*(r2-r_plus))/(2*sqrt(A*B*(1-E**2))*(r_plus-r1)*(r2-r_plus)*(A*(r1-r_plus)-B*(r2-r_plus)))*_ellippiinc(xi_r,1/D_plus**2,k_r) \ - - (sqrt(r2-r1)*log( - ((D_plus*sqrt(1-D_plus**2*k_r**2)+dn*sn)**2 + (k_r*(D_plus**2-sn**2))**2) / - ((D_plus*sqrt(1-D_plus**2*k_r**2)-dn*sn)**2 + (k_r*(D_plus**2-sn**2))**2) - ) / - (4*sqrt((1-E**2)*(r2-r_plus)*(r_plus-r1))*sqrt(A**2*(r_plus-r1)-(r2-r_plus)*(r1**2-B**2+r2*r_plus-r1*(r2+r_plus)))) - ) - - def I_r_minus(q_r): - # equation 43 - sn, cn, dn, xi_r = ellipj(2*ellipk(k_r**2)*q_r/pi,k_r**2) - mino_time = q_r/upsilon_r - # equation 48 - return (A-B)*mino_time/(A*(r1-r_minus)-B*(r2-r_minus)) \ - + (r2-r1)*(A*(r1-r_minus)+B*(r2-r_minus))/(2*sqrt(A*B*(1-E**2))*(r_minus-r1)*(r2-r_minus)*(A*(r1-r_minus)-B*(r2-r_minus)))*_ellippiinc(xi_r,1/D_minus**2,k_r) \ - - (sqrt(r2-r1)*log( - ((D_minus*sqrt(1-D_minus**2*k_r**2)+dn*sn)**2 + (k_r*(D_minus**2-sn**2))**2) / - ((D_minus*sqrt(1-D_minus**2*k_r**2)-dn*sn)**2 + (k_r*(D_minus**2-sn**2))**2) - ) / - (4*sqrt((1-E**2)*(r2-r_minus)*(r_minus-r1))*sqrt(A**2*(r_minus-r1)-(r2-r_minus)*(r1**2-B**2+r2*r_minus-r1*(r2+r_minus)))) - ) - - return I_r, I_r2, I_r_plus, I_r_minus - -def plunging_radial_solutions_complex(a,E,L,Q): - r""" - Computes the radial solutions :math:`r(q_r), t_r(q_r), \phi_r(q_r)` from equation 50 and 51 of `Dyson and van de Meent `_ for the case of two complex radial roots. - - :param a: dimensionless spin parameter - :type a: double - :param E: dimensionless energy - :type E: double - :param L: dimensionless angular momentum - :type L: double - :param Q: dimensionless Carter constant - :type Q: double - - :return: tuple of radial solutions :math:`(r(q_r), t_r(q_r), \phi_r(q_r))` - :rtype: tuple(function, function, function) - """ - roots = plunging_radial_roots(a,E,L,Q) - if np.iscomplex(roots).sum() != 2: raise ValueError("There should be two complex roots") - - r1, r2, r3, r4 = roots # r1 < r2 are real and r3/r4 are complex conjugates - rho_r = np.real(r3) - rho_i = np.imag(r4) - - # inner and outer horizons - r_plus = 1+sqrt(1-a**2) - r_minus = 1-sqrt(1-a**2) - - # equation 42 - A = sqrt((r2-rho_r)**2+rho_i**2) - B = sqrt((r1-rho_r)**2+rho_i**2) - k_r = sqrt(((r2-r1)**2-(A-B)**2)/(4*A*B)) - - upsilon_r, upsilon_theta = plunging_mino_frequencies(a,E,L,Q) - - I_r, I_r2, I_r_plus, I_r_minus = plunging_radial_integrals(a,E,L,Q) - - def r(q_r): - # equation 43 - sn, cn, dn, xi_r = ellipj(2*ellipk(k_r**2)*q_r/pi,k_r**2) - # equation 49 - return ((A-B)*(A*r1-B*r2)*sn**2+2*A*B*(r1+r2)-2*A*B*(r2-r1)*cn)/(4*A*B+(A-B)**2*sn**2) - - def t_r(q_r): - mino_time = q_r/upsilon_r - # equation 41 - return E*(r_plus**2+r_minus**2+r_plus*r_minus+2*a**2)*mino_time + \ - E*(I_r2(q_r)+I_r(q_r)*(r_minus+r_plus)) + \ - ((r_minus**2+a**2)*(E*(r_minus**2+a**2)-a*L)/(r_minus-r_plus)*I_r_minus(q_r) + - (r_plus**2+a**2)*(E*(r_plus**2+a**2)-a*L)/(r_plus-r_minus)*I_r_plus(q_r) - ) - - def phi_r(q_r): - mino_time = q_r/upsilon_r - # equation 40 - return a*( - (E*(r_minus**2+a**2)-a*L)/(r_minus-r_plus)*I_r_minus(q_r) + - (E*(r_plus**2+a**2)-a*L)/(r_plus-r_minus)*I_r_plus(q_r) - ) - - return r, t_r, phi_r - -def plunging_polar_solutions(a,E,L,Q): - r""" - Computes the polar solutions :math:`\theta(q_\theta), t_\theta(q_\theta), \phi_\theta(q_\theta)` from equation 33 and 37 of `Dyson and van de Meent `_ - - :param a: dimensionless spin parameter - :type a: double - :param E: dimensionless energy - :type E: double - :param L: dimensionless angular momentum - :type L: double - :param Q: dimensionless Carter constant - :type Q: double - - :return: tuple of polar solutions :math:`(\theta(q_\theta), t_\theta(q_\theta), \phi_\theta(q_\theta))` - :rtype: tuple(function, function, function) - """ - z1 = sqrt(1/2*(1+(L**2+Q)/(a**2*(1-E**2))-sqrt((1+(L**2+Q)/(a**2*(1-E**2)))**2-4*Q/(a**2*(1-E**2))))) - z2 = sqrt(a**2*(1-E**2)/2*(1+(L**2+Q)/(a**2*(1-E**2))+sqrt((1+(L**2+Q)/(a**2*(1-E**2)))**2-4*Q/(a**2*(1-E**2))))) - k_theta = a*sqrt(1-E**2)*z1/z2 - - upsilon_r, upsilon_theta = plunging_mino_frequencies(a,E,L,Q) - - def theta(q_theta): - mino_time = q_theta/upsilon_theta - # equation 28 - sn, cn, dn, xi_theta = ellipj(z2*mino_time,k_theta**2) - # equation 27 - return arccos(z1*sn) - - def t_theta(q_theta): - mino_time = q_theta/upsilon_theta - # equation 28 - sn, cn, dn, xi_theta = ellipj(z2*mino_time,k_theta**2) - # equation 35 - return E/(1-E**2)*((z2**2-a**2*(1-E**2))*mino_time-z2*ellipeinc(xi_theta,k_theta**2)) - - def phi_theta(q_theta): - mino_time = q_theta/upsilon_theta - # equation 28 - sn, cn, dn, xi_theta = ellipj(z2*mino_time,k_theta**2) - # equation 31 - return L/z2*_ellippiinc(xi_theta,z1**2,k_theta) - - return theta, t_theta, phi_theta diff --git a/kerrgeopy/spacetime.py b/kerrgeopy/spacetime.py index d963406..40034f1 100644 --- a/kerrgeopy/spacetime.py +++ b/kerrgeopy/spacetime.py @@ -16,9 +16,13 @@ class KerrSpacetime: :ivar a: dimensionless angular momentum :ivar M: mass of the black hole + :ivar inner_horizon: radius of the inner horizon + :ivar outer_horizon: radius of the outer horizon """ def __init__(self,a,M=None): self.a, self.M = a, M + self.inner_horizon = 1-sqrt(1-self.a**2) + self.outer_horizon = 1+sqrt(1-self.a**2) def separatrix(self,e,x): """ @@ -35,7 +39,7 @@ def separatrix(self,e,x): def is_stable(self,p,e,x): """ - Determines whether a given orbit is stable or not. + Determines whether a given orbit is stable or not :param p: dimensionless semi-latus rectum :type p: double @@ -48,24 +52,6 @@ def is_stable(self,p,e,x): """ return is_stable(self.a,p,e,x) - def inner_horizon(self): - """ - Computes the radius of the inner event horizon - - :return: dimensionless radius of the inner event horizon - :rtype: double - """ - return 1-sqrt(1-self.a**2) - - def outer_horizon(self): - """ - Computes the radius of the outer event horizon - - :return: dimensionless radius of the outer event horizon - :rtype: double - """ - return 1+sqrt(1-self.a**2) - def metric(self,t,r,theta,phi): """ Returns the matrix representation of the metric at a given point expressed in Boyer-Lindquist coordinates. @@ -149,24 +135,24 @@ def four_velocity(self,t,r,theta,phi,constants,upsilon_r,upsilon_theta,initial_p # Z = lambda z: Q-z**2*(a**2*(1-E**2)*(1-z**2)+L**2+Q) # Z = Polynomial([Q,-(Q+a**2*(1-E**2)+L**2),a**2*(1-E**2)]) - def t_prime(time): + def u_t(time): delta = r(time)**2-2*r(time)+a**2 sigma = r(time)**2+a**2*cos(theta(time))**2 return 1/sigma*((r(time)**2+a**2)/delta*(E*(r(time)**2+a**2)-a*L)-a**2*E*(1-cos(theta(time))**2)+a*L) - def r_prime(time): + def u_r(time): q_r = upsilon_r*time + q_r0 sigma = r(time)**2+a**2*cos(theta(time))**2 return np.copysign(1,sin(q_r))*sqrt(abs(R(r(time))))/sigma - def theta_prime(time): + def u_theta(time): q_theta = upsilon_theta*time + q_theta0 sigma = r(time)**2+a**2*cos(theta(time))**2 return np.copysign(1,sin(q_theta))*sqrt(abs(Z(cos(theta(time)))))/(sigma*sin(theta(time))) - def phi_prime(time): + def u_phi(time): sigma = r(time)**2+a**2*cos(theta(time))**2 delta = r(time)**2-2*r(time)+a**2 return 1/sigma*(a/delta*(E*(r(time)**2+a**2)-a*L)+L/(1-cos(theta(time))**2)-a*E) - return t_prime, r_prime, theta_prime, phi_prime \ No newline at end of file + return u_t, u_r, u_theta, u_phi \ No newline at end of file diff --git a/kerrgeopy/stable.py b/kerrgeopy/stable.py new file mode 100644 index 0000000..10741b3 --- /dev/null +++ b/kerrgeopy/stable.py @@ -0,0 +1,319 @@ +""" +Module implementing the stable bound orbit solutions of `Fujita and Hikida `_ +""" +from .constants import * +from .frequencies import _ellippi, _ellippiinc +from .frequencies import * +from scipy.special import ellipj, ellipeinc +from .orbit import Orbit +from numpy import pi, arccos + +class StableOrbit(Orbit): + r""" + Class representing a stable bound orbit in Kerr spacetime. + + :param a: dimensionless angular momentum (must satisfy 0 <= a < 1) + :type a: double + :param p: semi-latus rectum + :type p: double + :param e: orbital eccentricity (must satisfy 0 <= e < 1) + :type e: double + :param x: cosine of the orbital inclination (must satisfy 0 < x^2 <= 1) + :type x: double + :param initial_phases: tuple of initial phases :math:`(q_{t_0},q_{r_0},q_{\theta_0},q_{\phi_0})`, defaults to (0,0,0,0) + :type initial_phases: tuple, optional + :param M: mass of the primary in solar masses, optional + :type M: double + :param mu: mass of the smaller body in solar masses, optional + :type mu: double + + :ivar a: dimensionless angular momentum + :ivar p: semi-latus rectum + :ivar e: orbital eccentricity + :ivar x: cosine of the orbital inclination + :ivar M: mass of the primary in solar masses + :ivar mu: mass of the smaller body in solar masses + :ivar E: dimensionless energy + :ivar L: dimensionless angular momentum + :ivar Q: dimensionless carter constant + :ivar initial_position: tuple of initial position coordinates :math:`(t_0, r_0, \theta_0, \phi_0)` + :ivar initial_velocity: tuple of initial four-velocity components :math:`(u^t_0, u^r_0, u^\theta_0, u^\phi_0)` + :ivar upsilon_r: dimensionless radial orbital frequency in Mino time + :ivar upsilon_theta: dimensionless polar orbital frequency in Mino time + :ivar upsilon_phi: dimensionless azimuthal orbital frequency in Mino time + :ivar gamma: dimensionless time dilation factor + :ivar omega_r: dimensionless radial orbital frequency in Boyer-Lindquist time + :ivar omega_theta: dimensionless polar orbital frequency in Boyer-Lindquist time + :ivar omega_phi: dimensionless azimuthal orbital frequency in Boyer-Lindquist time + """ + def __init__(self,a,p,e,x,initial_phases=(0,0,0,0),M=None,mu=None): + self.a, self.p, self.e, self.x, self.initial_phases, self.M, self.mu = a, p, e, x, initial_phases, M, mu + constants = constants_of_motion(a,p,e,x) + + self.E, self.L, self.Q = constants + self.upsilon_r, self.upsilon_theta, self.upsilon_phi, self.gamma = mino_frequencies(a,p,e,x) + self.omega_r, self.omega_theta, self.omega_phi = fundamental_frequencies(a,p,e,x) + self.stable = True + + u_t, u_r, u_theta, u_phi = self.four_velocity() + t, r, theta, phi = self.trajectory() + self.initial_position = t(0), r(0), theta(0), phi(0) + self.initial_velocity = u_t(0), u_r(0), u_theta(0), u_phi(0) + + @classmethod + def from_constants(cls,a,E,L,Q,initial_phases=(0,0,0,0),M=None,mu=None): + """ + Alternative constructor method that takes the constants of motion as arguments. + + :param a: spin parameter + :type a: double + :param E: dimensionless energy + :type E: double + :param L: dimensionless angular momentum + :type L: double + :param Q: dimensionless Carter constant + :type Q: double + :param M: mass of the primary in solar masses + :type M: double, optional + :param mu: mass of the smaller body in solar masses + :type mu: double, optional + """ + a, p, e, x = apex_from_constants(a,E,L,Q) + + return cls(a,p,e,x,initial_phases,M,mu) + + def mino_frequencies(self, units="natural"): + r""" + Computes orbital frequencies in Mino time. Returns dimensionless frequencies in geometrized units by default. + M and mu must be defined in order to convert to physical units. + + :param units: units to return the frequencies in (options are "natural", "mks" and "cgs"), defaults to "natural" + :type units: str, optional + + :return: tuple of orbital frequencies in the form :math:`(\Upsilon_r, \Upsilon_\theta, \Upsilon_\phi, \Gamma)` + :rtype: tuple(double, double, double, double) + """ + upsilon_r, upsilon_theta, upsilon_phi, gamma = self.upsilon_r, self.upsilon_theta, self.upsilon_phi, self.gamma + if units == "natural": + return upsilon_r, upsilon_theta, upsilon_phi, gamma + + if self.M is None: raise ValueError("M must be specified to convert frequencies to physical units") + + if units == "mks" or units == "cgs": + return time_in_seconds(upsilon_r,self.M), time_in_seconds(upsilon_theta,self.M), time_in_seconds(upsilon_phi,self.M), time2_in_seconds2(gamma,self.M) + + raise ValueError("units must be one of 'natural', 'mks', or 'cgs'") + + def fundamental_frequencies(self, units="natural"): + r""" + Computes orbital frequencies in Boyer-Lindquist time. Returns dimensionless frequencies in geometrized units by default. + M and mu must be defined in order to convert to physical units. + + :param units: units to return the frequencies in (options are "natural", "mks", "cgs" and "mHz"), defaults to "natural" + :type units: str, optional + :return: tuple of orbital frequencies in the form :math:`(\Omega_r, \Omega_\theta, \Omega_\phi)` + :rtype: tuple(double, double, double) + """ + upsilon_r, upsilon_theta, upsilon_phi, gamma = self.upsilon_r, self.upsilon_theta, self.upsilon_phi, self.gamma + if units == "natural": + return upsilon_r/gamma, upsilon_theta/gamma, upsilon_phi/gamma + + if self.M is None: raise ValueError("M must be specified to convert frequencies to physical units") + + if units == "mks" or units == "cgs": + return frequency_in_Hz(upsilon_r/gamma,self.M), frequency_in_Hz(upsilon_theta/gamma,self.M), frequency_in_Hz(upsilon_phi/gamma,self.M) + if units == "mHz": + return frequency_in_mHz(upsilon_r/gamma,self.M), frequency_in_mHz(upsilon_theta/gamma,self.M), frequency_in_mHz(upsilon_phi/gamma,self.M) + + raise ValueError("units must be one of 'natural', 'mks', 'cgs', or 'mHz'") + + def trajectory(self,initial_phases=None,distance_units="natural",time_units="natural"): + r""" + Computes the time, radial, polar, and azimuthal coordinates of the orbit as a function of mino time. + + :param initial_phases: tuple of initial phases :math:`(q_{t_0},q_{r_0},q_{\theta_0},q_{\phi_0})` + :type initial_phases: tuple, optional + :param distance_units: units to compute the radial component of the trajectory in (options are "natural", "mks", "cgs", "au" and "km"), defaults to "natural" + :type distance_units: str, optional + :param time_units: units to compute the time component of the trajectory in (options are "natural", "mks", "cgs", and "days"), defaults to "natural" + :type time_units: str, optional + + :return: tuple of functions in the form :math:`(t(\lambda), r(\lambda), \theta(\lambda), \phi(\lambda))` + :rtype: tuple(function, function, function, function) + """ + if initial_phases is None: initial_phases = self.initial_phases + + return stable_trajectory(self.a,self.p,self.e,self.x,initial_phases,self.M,distance_units,time_units) + + +def radial_solutions(a,constants,radial_roots): + r""" + Computes the radial solutions :math:`r(q_r), t^{(r)}(q_r), \phi^{(r)}(q_r)` from equation 6 of `Fujita and Hikida `_. + :math:`q_r` is defined as :math:`q_r = \Upsilon_r \lambda = 2\pi \frac{\lambda}{\Lambda_r}`. + Assumes the initial conditions :math:`r(0) = r_{\text{min}}` and :math:`\theta(0) = \theta_{\text{min}}`. + + :param a: dimensionless spin parameter + :type a: double + :param constants: tuple of constants :math:`(E,L,Q)` + :type constants: tuple(double,double,double) + :param radial_roots: tuple of roots :math:`(r_1,r_2,r_3,r_4)`. Assumes that motion is between :math:`r_1` and :math:`r_2` and that roots are otherwise in decreasing order. + :type radial_roots: tuple(double,double,double,double) + + :return: tuple of functions in the form :math:`(r, t^{(r)}, \phi^{(r)})` + :rtype: tuple(function, function, function) + """ + E, L, Q = constants + r1, r2, r3, r4 = radial_roots + + r_plus = 1+sqrt(1-a**2) + r_minus = 1-sqrt(1-a**2) + + h_r = (r1-r2)/(r1-r3) + h_plus = (r1-r2)*(r3-r_plus)/((r1-r3)*(r2-r_plus)) + h_minus = (r1-r2)*(r3-r_minus)/((r1-r3)*(r2-r_minus)) + + # equation 13 + k_r = sqrt((r1-r2)*(r3-r4)/((r1-r3)*(r2-r4))) + + def r(q_r): + # equation 27 + u_r = ellipk(k_r**2)*q_r/pi + + sn, cn, dn, psi_r = ellipj(u_r,k_r**2) + return (r3*(r1-r2)*sn**2-r2*(r1-r3))/((r1-r2)*sn**2-(r1-r3)) + + def t_r(q_r): + # equation 27 + u_r = ellipk(k_r**2)*q_r/pi + sn, cn, dn, psi_r = ellipj(u_r,k_r**2) + # equation 28 + return 2/sqrt((1-E**2)*(r1-r3)*(r2-r4))* \ + ( + E/2*( + (r2-r3)*(r1+r2+r3+r4)*(_ellippiinc(psi_r,h_r,k_r**2)-q_r/pi*_ellippi(h_r,k_r**2)) \ + + (r1-r3)*(r2-r4)*(ellipeinc(psi_r,k_r**2)+h_r*sn*cn*sqrt(1-k_r**2*sn**2)/(h_r*sn**2-1) - q_r/pi*ellipe(k_r**2)) + ) + + 2*E*(r2-r3)*(_ellippiinc(psi_r,h_r,k_r**2)-q_r/pi*_ellippi(h_r,k_r**2)) + - 2/(r_plus-r_minus) * \ + ( + ((4*E-a*L)*r_plus-2*a**2*E)*(r2-r3)/((r3-r_plus)*(r2-r_plus))*(_ellippiinc(psi_r,h_plus,k_r**2)-q_r/pi*_ellippi(h_plus,k_r**2)) + - ((4*E-a*L)*r_minus-2*a**2*E)*(r2-r3)/((r3-r_minus)*(r2-r_minus))*(_ellippiinc(psi_r,h_minus,k_r**2)-q_r/pi*_ellippi(h_minus,k_r**2)) + ) + ) + + def phi_r(q_r): + # equation 27 + u_r = ellipk(k_r**2)*q_r/pi + sn, cn, dn, psi_r = ellipj(u_r,k_r**2) + # equation 28 + return -2*a/((r_plus-r_minus)*sqrt((1-E**2)*(r1-r3)*(r2-r4))) * \ + ( + (2*E*r_plus-a*L)*(r2-r3)/((r3-r_plus)*(r2-r_plus))*(_ellippiinc(psi_r,h_plus,k_r**2)-q_r/pi*_ellippi(h_plus,k_r**2)) + - (2*E*r_minus-a*L)*(r2-r3)/((r3-r_minus)*(r2-r_minus))*(_ellippiinc(psi_r,h_minus,k_r**2)-q_r/pi*_ellippi(h_minus,k_r**2)) + ) + + return r, t_r, phi_r + +def polar_solutions(a,constants,polar_roots): + r""" + Computes the polar solutions :math:`\theta(q_\theta), t^{(\theta)}(q_\theta), \phi^{(\theta)}(q_\theta)` from equation 6 of `Fujita and Hikida `_. + :math:`q_\theta` is defined as :math:`q_\theta = \Upsilon_\theta \lambda = 2\pi \frac{\lambda}{\Lambda_\theta}`. + Assumes the initial conditions :math:`r(0) = r_{\text{min}}` and :math:`\theta(0) = \theta_{\text{min}}`. + + :param a: dimensionless spin parameter + :type a: double + :param constants: tuple of constants :math:`(E,L,Q)` + :type constants: tuple(double,double,double) + :param polar_roots: tuple of roots :math:`(z_-,z_+)` + :type polar_roots: tuple(double,double) + + :return: tuple of functions in the form :math:`(\theta, t^{(\theta)}, \phi^{(\theta)})` + :rtype: tuple(function, function, function) + """ + E, L, Q = constants + z_minus, z_plus = polar_roots + epsilon0 = a**2*(1-E**2)/L**2 + # simplified form of epsilon0*z_plus + e0zp = (a**2*(1-E**2)*(1-z_minus)+L**2)/(L**2*(1-z_minus)) + # simplified form of a**2*sqrt(z_plus/epsilon0) + a2sqrt_zp_over_e0 = L**2/((1-E**2)*sqrt(1-z_minus)) if a == 0 else a**2*z_plus/sqrt(epsilon0*z_plus) + + # equation 13 + k_theta = 0 if a == 0 else sqrt(z_minus/z_plus) + + def theta(q_theta): + u_theta = 2/pi*ellipk(k_theta**2)*(q_theta+pi/2) + sn, cn, dn, ph = ellipj(u_theta,k_theta**2) + # equation 38 + return arccos(sqrt(z_minus)*sn) + + def t_theta(q_theta): + u_theta = 2/pi*ellipk(k_theta**2)*(q_theta+pi/2) + sn, cn, dn, psi_theta = ellipj(u_theta,k_theta**2) + # equation 39 + return sign(L)*a2sqrt_zp_over_e0*E/L*(2/pi*ellipe(k_theta**2)*(q_theta+pi/2)-ellipeinc(psi_theta,k_theta**2)) + + def phi_theta(q_theta): + sn, cn, dn, psi_theta = ellipj(2/pi*ellipk(k_theta**2)*(q_theta+pi/2),k_theta**2) + # equation 39 + return sign(L)*1/sqrt(e0zp)*(_ellippiinc(psi_theta,z_minus,k_theta**2)-2/pi*_ellippi(z_minus,k_theta**2)*(q_theta+pi/2)) + + return theta, t_theta, phi_theta + +def stable_trajectory(a,p,e,x,initial_phases=(0,0,0,0),M=None,distance_units="natural",time_units="natural"): + r""" + Computes the time, radial, polar, and azimuthal coordinates of the orbit with the given parameters as a function of mino time. + + :param a: dimensionless spin parameter + :type a: double + :param p: semi-latus rectum + :type p: double + :param e: orbital eccentricity + :type e: double + :param x: cosine of the orbital inclination + :type x: double + :param initial_phases: tuple of initial phases :math:`(q_{t_0},q_{r_0},q_{\theta_0},q_{\phi_0})`, defaults to (0,0,0,0) + :type initial_phases: tuple, optional + :param M: mass of the primary in solar masses + :type M: double, optional + :param distance_units: units to compute the radial component of the trajectory in (options are "natural", "mks", "cgs", "au" and "km"), defaults to "natural" + :type distance_units: str, optional + :param time_units: units to compute the time component of the trajectory in (options are "natural", "mks", "cgs", and "days"), defaults to "natural" + :type time_units: str, optional + + :return: tuple of functions in the form :math:`(t(\lambda), r(\lambda), \theta(\lambda), \phi(\lambda))` + :rtype: tuple(function, function, function, function) + """ + if ((distance_units != "natural") or (time_units != "natural")) and M is None: raise ValueError("M must be specified to convert to physical units") + + upsilon_r, upsilon_theta, upsilon_phi, gamma = mino_frequencies(a,p,e,x) + constants = constants_of_motion(a,p,e,x) + radial_roots = stable_radial_roots(a,p,e,x,constants) + polar_roots = stable_polar_roots(a,p,e,x,constants) + + r_phases, t_r, phi_r = radial_solutions(a,constants,radial_roots) + theta_phases, t_theta, phi_theta = polar_solutions(a,constants,polar_roots) + q_t0, q_r0, q_theta0, q_phi0 = initial_phases + + distance_conversion_func = {"natural": lambda x,M: x, "mks": distance_in_meters, "cgs": distance_in_cm, "au": distance_in_au,"km": distance_in_km} + time_conversion_func = {"natural": lambda x,M: x, "mks": time_in_seconds, "cgs": time_in_seconds, "days": time_in_days} + + # Calculate normalization constants so that t = 0 and phi = 0 at lambda = 0 when q_t0 = 0 and q_phi0 = 0 + C_t = t_r(q_r0)+t_theta(q_theta0) + C_phi= phi_r(q_r0)+phi_theta(q_theta0) + + def t(mino_time): + # equation 6 + return time_conversion_func[time_units](q_t0 + gamma*mino_time + t_r(upsilon_r*mino_time+q_r0) + t_theta(upsilon_theta*mino_time+q_theta0) - C_t, M) + + def r(mino_time): + return distance_conversion_func[distance_units](r_phases(upsilon_r*mino_time+q_r0),M) + + def theta(mino_time): + return theta_phases(upsilon_theta*mino_time+q_theta0) + + def phi(mino_time): + # equation 6 + return q_phi0 + upsilon_phi*mino_time + phi_r(upsilon_r*mino_time+q_r0) + phi_theta(upsilon_theta*mino_time+q_theta0) - C_phi + + return t, r, theta, phi \ No newline at end of file diff --git a/kerrgeopy/stable_orbit.py b/kerrgeopy/stable_orbit.py deleted file mode 100644 index 9086138..0000000 --- a/kerrgeopy/stable_orbit.py +++ /dev/null @@ -1,182 +0,0 @@ -""" -Module containing the BoundOrbit class -""" -from .constants import * -from .frequencies import * -from .stable_solutions import * -from .initial_conditions import apex_from_constants -from .frequencies import _radial_roots, _polar_roots -from .orbit import Orbit - -class StableOrbit(Orbit): - r""" - Class representing a stable bound orbit in Kerr spacetime. - - :param a: dimensionless angular momentum (must satisfy 0 <= a < 1) - :type a: double - :param p: semi-latus rectum - :type p: double - :param e: orbital eccentricity (must satisfy 0 <= e < 1) - :type e: double - :param x: cosine of the orbital inclination (must satisfy 0 < x^2 <= 1) - :type x: double - :param M: mass of the primary in solar masses, optional - :type M: double - :param mu: mass of the smaller body in solar masses, optional - :type mu: double - - :ivar a: dimensionless angular momentum - :ivar p: semi-latus rectum - :ivar e: orbital eccentricity - :ivar x: cosine of the orbital inclination - :ivar M: mass of the primary in solar masses - :ivar mu: mass of the smaller body in solar masses - :ivar E: dimensionless energy - :ivar L: dimensionless angular momentum - :ivar Q: dimensionless carter constant - :ivar initial_position: tuple of initial position coordinates :math:`(t_0, r_0, \theta_0, \phi_0)` - :ivar initial_velocity: tuple of initial four-velocity components :math:`(u^t_0, u^r_0, u^\theta_0, u^\phi_0)` - :ivar upsilon_r: dimensionless radial orbital frequency in Mino time - :ivar upsilon_theta: dimensionless polar orbital frequency in Mino time - :ivar upsilon_phi: dimensionless azimuthal orbital frequency in Mino time - :ivar gamma: dimensionless time dilation factor - :ivar omega_r: dimensionless radial orbital frequency in Boyer-Lindquist time - :ivar omega_theta: dimensionless polar orbital frequency in Boyer-Lindquist time - :ivar omega_phi: dimensionless azimuthal orbital frequency in Boyer-Lindquist time - """ - def __init__(self,a,p,e,x,initial_phases=(0,0,0,0),M = None,mu=None): - """ - Constructor method - """ - self.a, self.p, self.e, self.x, self.initial_phases, self.M, self.mu = a, p, e, x, initial_phases, M, mu - constants = constants_of_motion(a,p,e,x) - - self.E, self.L, self.Q = constants - self.upsilon_r, self.upsilon_theta, self.upsilon_phi, self.gamma = mino_frequencies(a,p,e,x) - self.omega_r, self.omega_theta, self.omega_phi = fundamental_frequencies(a,p,e,x) - self.stable = True - - u_t, u_r, u_theta, u_phi = self.four_velocity() - t, r, theta, phi = self.trajectory() - self.initial_position = t(0), r(0), theta(0), phi(0) - self.initial_velocity = u_t(0), u_r(0), u_theta(0), u_phi(0) - - @classmethod - def from_constants(cls,a,E,L,Q,initial_phases=(0,0,0,0),M=None,mu=None): - """ - Alternative constructor method that takes the constants of motion as arguments. - - :param a: spin parameter - :type a: double - :param E: dimensionless energy - :type E: double - :param L: dimensionless angular momentum - :type L: double - :param Q: dimensionless Carter constant - :type Q: double - :param M: mass of the primary in solar masses - :type M: double, optional - :param mu: mass of the smaller body in solar masses - :type mu: double, optional - """ - a, p, e, x = apex_from_constants(a,E,L,Q) - - return cls(a,p,e,x,initial_phases,M,mu) - - def mino_frequencies(self, units="natural"): - r""" - Computes orbital frequencies in Mino time. Returns dimensionless frequencies in geometrized units by default. - M and mu must be defined in order to convert to physical units. - - :param units: units to return the frequencies in (options are "natural", "mks" and "cgs"), defaults to "natural" - :type units: str, optional - - :return: tuple of orbital frequencies in the form :math:`(\Upsilon_r, \Upsilon_\theta, \Upsilon_\phi, \Gamma)` - :rtype: tuple(double, double, double, double) - """ - upsilon_r, upsilon_theta, upsilon_phi, gamma = self.upsilon_r, self.upsilon_theta, self.upsilon_phi, self.gamma - if units == "natural": - return upsilon_r, upsilon_theta, upsilon_phi, gamma - - if self.M is None: raise ValueError("M must be specified to convert frequencies to physical units") - - if units == "mks" or units == "cgs": - return time_in_seconds(upsilon_r,self.M), time_in_seconds(upsilon_theta,self.M), time_in_seconds(upsilon_phi,self.M), time2_in_seconds2(gamma,self.M) - - raise ValueError("units must be one of 'natural', 'mks', or 'cgs'") - - def fundamental_frequencies(self, units="natural"): - r""" - Computes orbital frequencies in Boyer-Lindquist time. Returns dimensionless frequencies in geometrized units by default. - M and mu must be defined in order to convert to physical units. - - :param units: units to return the frequencies in (options are "natural", "mks", "cgs" and "mHz"), defaults to "natural" - :type units: str, optional - :return: tuple of orbital frequencies in the form :math:`(\Omega_r, \Omega_\theta, \Omega_\phi)` - :rtype: tuple(double, double, double) - """ - upsilon_r, upsilon_theta, upsilon_phi, gamma = self.upsilon_r, self.upsilon_theta, self.upsilon_phi, self.gamma - if units == "natural": - return upsilon_r/gamma, upsilon_theta/gamma, upsilon_phi/gamma - - if self.M is None: raise ValueError("M must be specified to convert frequencies to physical units") - - if units == "mks" or units == "cgs": - return frequency_in_Hz(upsilon_r/gamma,self.M), frequency_in_Hz(upsilon_theta/gamma,self.M), frequency_in_Hz(upsilon_phi/gamma,self.M) - if units == "mHz": - return frequency_in_mHz(upsilon_r/gamma,self.M), frequency_in_mHz(upsilon_theta/gamma,self.M), frequency_in_mHz(upsilon_phi/gamma,self.M) - - raise ValueError("units must be one of 'natural', 'mks', 'cgs', or 'mHz'") - - def trajectory(self,initial_phases=None,distance_units="natural",time_units="natural"): - r""" - Computes the time, radial, polar, and azimuthal coordinates of the orbit as a function of mino time. - - :param initial_phases: tuple of initial phases :math:`(q_{t_0},q_{r_0},q_{\theta_0},q_{\phi_0})`, defaults to (0,0,0,0) - :type initial_phases: tuple, optional - :param distance_units: units to compute the radial component of the trajectory in (options are "natural", "mks", "cgs", "au" and "km"), defaults to "natural" - :type distance_units: str, optional - :param time_units: units to compute the time component of the trajectory in (options are "natural", "mks", "cgs", and "days"), defaults to "natural" - :type time_units: str, optional - - :return: tuple of functions in the form :math:`(t(\lambda), r(\lambda), \theta(\lambda), \phi(\lambda))` - :rtype: tuple(function, function, function, function) - """ - if initial_phases is None: initial_phases = self.initial_phases - if distance_units != "natural" and (self.M is None or self.mu is None): raise ValueError("M and mu must be specified to convert r to physical units") - if time_units != "natural" and (self.M is None or self.mu is None): raise ValueError("M and mu must be specified to convert t to physical units") - - a, p, e, x = self.a, self.p, self.e, self.x - upsilon_r, upsilon_theta, upsilon_phi, gamma = self.upsilon_r, self.upsilon_theta, self.upsilon_phi, self.gamma - constants = self.constants_of_motion() - radial_roots = _radial_roots(a,p,e,constants) - polar_roots = _polar_roots(a,x,constants) - - r_phases, t_r, phi_r = radial_solutions(a,constants,radial_roots) - theta_phases, t_theta, phi_theta = polar_solutions(a,constants,polar_roots) - q_t0, q_r0, q_theta0, q_phi0 = initial_phases - - distance_conversion_func = {"natural": lambda x,M: x, "mks": distance_in_meters, "cgs": distance_in_cm, "au": distance_in_au,"km": distance_in_km} - time_conversion_func = {"natural": lambda x,M: x, "mks": time_in_seconds, "cgs": time_in_seconds, "days": time_in_days} - - # Calculate normalization constants so that t = 0 and phi = 0 at lambda = 0 when q_t0 = 0 and q_phi0 = 0 - C_t = t_r(q_r0)+t_theta(q_theta0) - C_phi= phi_r(q_r0)+phi_theta(q_theta0) - - def t(mino_time): - # equation 6 - return time_conversion_func[time_units](q_t0 + gamma*mino_time + t_r(upsilon_r*mino_time+q_r0) + t_theta(upsilon_theta*mino_time+q_theta0) - C_t, self.M) - - def r(mino_time): - return distance_conversion_func[distance_units](r_phases(upsilon_r*mino_time+q_r0),self.M) - - def theta(mino_time): - return theta_phases(upsilon_theta*mino_time+q_theta0) - - def phi(mino_time): - # equation 6 - return q_phi0 + upsilon_phi*mino_time + phi_r(upsilon_r*mino_time+q_r0) + phi_theta(upsilon_theta*mino_time+q_theta0) - C_phi - - return t, r, theta, phi - - \ No newline at end of file diff --git a/kerrgeopy/stable_solutions.py b/kerrgeopy/stable_solutions.py deleted file mode 100644 index 83feb22..0000000 --- a/kerrgeopy/stable_solutions.py +++ /dev/null @@ -1,144 +0,0 @@ -""" -Module containing functions to compute intermediate terms in the stable bound orbit solutions of `Fujita and Hikida `_ -""" -from .constants import * -from .frequencies_from_constants import _ellippi -from .frequencies_from_constants import * -from scipy.special import ellipj, ellipeinc -from numpy import pi, sin, cos, arcsin, arccos, floor, where - -def _ellippiinc(phi,n,k): - r""" - Incomplete elliptic integral of the third kind defined as :math:`\Pi(\phi,n,k) = \int_0^{\phi} \frac{1}{1-n\sin^2\theta}\frac{1}{\sqrt{1-k^2\sin^2\theta}}d\theta`. - - :type phi: double - :type n: double - :type k: double - - :rtype: double - """ - # Note: sign of n is reversed from the definition in Fujita and Hikida - - # count the number of half periods - num_cycles = floor(phi/(pi/2)) - # map phi to [0,pi/2] - phi = abs(arcsin(sin(phi))) - # formula from https://en.wikipedia.org/wiki/Carlson_symmetric_form - integral = sin(phi)*elliprf(cos(phi)**2,1-k**2*sin(phi)**2,1)+1/3*n*sin(phi)**3*elliprj(cos(phi)**2,1-k**2*sin(phi)**2,1,1-n*sin(phi)**2) - - return where(num_cycles % 2 == 0, num_cycles*_ellippi(n,k)+integral, (num_cycles+1)*_ellippi(n,k)-integral) - -def radial_solutions(a,constants,radial_roots): - r""" - Computes the radial solutions :math:`r(q_r), t^{(r)}(q_r), \phi^{(r)}(q_r)` from equation 6 of `Fujita and Hikida `_. - :math:`q_r` is defined as :math:`q_r = \Upsilon_r \lambda = 2\pi \frac{\lambda}{\Lambda_r}`. - Assumes the initial conditions :math:`r(0) = r_{\text{min}}` and :math:`\theta(0) = \theta_{\text{min}}`. - - :param a: dimensionless spin parameter - :type a: double - :param constants: tuple of constants :math:`(E,L,Q)` - :type constants: tuple(double,double,double) - :param radial_roots: tuple of roots :math:`(r_1,r_2,r_3,r_4)`. Assumes that motion is between :math:`r_1` and :math:`r_2` and that roots are otherwise in decreasing order. - :type radial_roots: tuple(double,double,double,double) - - :return: tuple of functions in the form :math:`(r, t^{(r)}, \phi^{(r)})` - :rtype: tuple(function, function, function) - """ - E, L, Q = constants - r1, r2, r3, r4 = radial_roots - - r_plus = 1+sqrt(1-a**2) - r_minus = 1-sqrt(1-a**2) - - h_r = (r1-r2)/(r1-r3) - h_plus = (r1-r2)*(r3-r_plus)/((r1-r3)*(r2-r_plus)) - h_minus = (r1-r2)*(r3-r_minus)/((r1-r3)*(r2-r_minus)) - - # equation 13 - k_r = sqrt((r1-r2)*(r3-r4)/((r1-r3)*(r2-r4))) - - def r(q_r): - # equation 27 - u_r = ellipk(k_r**2)*q_r/pi - - sn, cn, dn, psi_r = ellipj(u_r,k_r**2) - return (r3*(r1-r2)*sn**2-r2*(r1-r3))/((r1-r2)*sn**2-(r1-r3)) - - def t_r(q_r): - # equation 27 - u_r = ellipk(k_r**2)*q_r/pi - sn, cn, dn, psi_r = ellipj(u_r,k_r**2) - # equation 28 - return 2/sqrt((1-E**2)*(r1-r3)*(r2-r4))* \ - ( - E/2*( - (r2-r3)*(r1+r2+r3+r4)*(_ellippiinc(psi_r,h_r,k_r)-q_r/pi*_ellippi(h_r,k_r)) \ - + (r1-r3)*(r2-r4)*(ellipeinc(psi_r,k_r**2)+h_r*sn*cn*sqrt(1-k_r**2*sn**2)/(h_r*sn**2-1) - q_r/pi*ellipe(k_r**2)) - ) - + 2*E*(r2-r3)*(_ellippiinc(psi_r,h_r,k_r)-q_r/pi*_ellippi(h_r,k_r)) - - 2/(r_plus-r_minus) * \ - ( - ((4*E-a*L)*r_plus-2*a**2*E)*(r2-r3)/((r3-r_plus)*(r2-r_plus))*(_ellippiinc(psi_r,h_plus,k_r)-q_r/pi*_ellippi(h_plus,k_r)) - - ((4*E-a*L)*r_minus-2*a**2*E)*(r2-r3)/((r3-r_minus)*(r2-r_minus))*(_ellippiinc(psi_r,h_minus,k_r)-q_r/pi*_ellippi(h_minus,k_r)) - ) - ) - - def phi_r(q_r): - # equation 27 - u_r = ellipk(k_r**2)*q_r/pi - sn, cn, dn, psi_r = ellipj(u_r,k_r**2) - # equation 28 - return -2*a/((r_plus-r_minus)*sqrt((1-E**2)*(r1-r3)*(r2-r4))) * \ - ( - (2*E*r_plus-a*L)*(r2-r3)/((r3-r_plus)*(r2-r_plus))*(_ellippiinc(psi_r,h_plus,k_r)-q_r/pi*_ellippi(h_plus,k_r)) - - (2*E*r_minus-a*L)*(r2-r3)/((r3-r_minus)*(r2-r_minus))*(_ellippiinc(psi_r,h_minus,k_r)-q_r/pi*_ellippi(h_minus,k_r)) - ) - - return r, t_r, phi_r - -def polar_solutions(a,constants,polar_roots): - r""" - Computes the polar solutions :math:`\theta(q_\theta), t^{(\theta)}(q_\theta), \phi^{(\theta)}(q_\theta)` from equation 6 of `Fujita and Hikida `_. - :math:`q_\theta` is defined as :math:`q_\theta = \Upsilon_\theta \lambda = 2\pi \frac{\lambda}{\Lambda_\theta}`. - Assumes the initial conditions :math:`r(0) = r_{\text{min}}` and :math:`\theta(0) = \theta_{\text{min}}`. - - :param a: dimensionless spin parameter - :type a: double - :param constants: tuple of constants :math:`(E,L,Q)` - :type constants: tuple(double,double,double) - :param polar_roots: tuple of roots :math:`(z_-,z_+)` - :type polar_roots: tuple(double,double) - - :return: tuple of functions in the form :math:`(\theta, t^{(\theta)}, \phi^{(\theta)})` - :rtype: tuple(function, function, function) - """ - E, L, Q = constants - z_minus, z_plus = polar_roots - epsilon0 = a**2*(1-E**2)/L**2 - # simplified form of epsilon0*z_plus - e0zp = (a**2*(1-E**2)*(1-z_minus)+L**2)/(L**2*(1-z_minus)) - # simplified form of a**2*sqrt(z_plus/epsilon0) - a2sqrt_zp_over_e0 = L**2/((1-E**2)*sqrt(1-z_minus)) if a == 0 else a**2*z_plus/sqrt(epsilon0*z_plus) - - # equation 13 - k_theta = 0 if a == 0 else sqrt(z_minus/z_plus) - - def theta(q_theta): - u_theta = 2/pi*ellipk(k_theta**2)*(q_theta+pi/2) - sn, cn, dn, ph = ellipj(u_theta,k_theta**2) - # equation 38 - return arccos(sqrt(z_minus)*sn) - - def t_theta(q_theta): - u_theta = 2/pi*ellipk(k_theta**2)*(q_theta+pi/2) - sn, cn, dn, psi_theta = ellipj(u_theta,k_theta**2) - # equation 39 - return sign(L)*a2sqrt_zp_over_e0*E/L*(2/pi*ellipe(k_theta**2)*(q_theta+pi/2)-ellipeinc(psi_theta,k_theta**2)) - - def phi_theta(q_theta): - sn, cn, dn, psi_theta = ellipj(2/pi*ellipk(k_theta**2)*(q_theta+pi/2),k_theta**2) - # equation 39 - return sign(L)*1/sqrt(e0zp)*(_ellippiinc(psi_theta,z_minus,k_theta)-2/pi*_ellippi(z_minus,k_theta)*(q_theta+pi/2)) - - return theta, t_theta, phi_theta - diff --git a/setup.py b/setup.py index 1fdde74..da782ff 100644 --- a/setup.py +++ b/setup.py @@ -3,9 +3,9 @@ setup( name="kerrgeopy", - version="0.2.0", + version="0.9.0", author="Seyong Park", - description="Library for computing bound and plunging geodesics in Kerr spacetime", + description="Library for computing stable and plunging geodesics in Kerr spacetime", url="https://github.com/syp2001/KerrGeoPy", packages=setuptools.find_packages(exclude=["tests"]), classifiers=( diff --git a/tests/data/README.txt b/tests/README.txt similarity index 68% rename from tests/data/README.txt rename to tests/README.txt index 24335d1..e38e925 100644 --- a/tests/data/README.txt +++ b/tests/README.txt @@ -1,21 +1,22 @@ -Many tests use data files containing a list of orbital parameters to test along with files containing output from the KerrGeodesics mathematica library. +Many tests use data files containing a list of orbital parameters to test along with files containing output from the KerrGeodesics mathematica library. +The jupyter notebook data/Unit Test Generation.ipynb contains the code that was used to generate these files. Below is a list of data files used by each test suite: test_constants.py: -const_values.txt - orbital parameters (a,p,e,x) +const_parameters.txt - orbital parameters (a,p,e,x) mathematica_const_output.txt - (E,L,Q) computed using KerrGeoConstantsOfMotion[a,p,e,x] for each orbit in const_values.txt -separatrix_values.txt - orbital parameters (a,p,e,x) +separatrix_parameters.txt - orbital parameters (a,p,e,x) mathematica_separatrix_output.txt - output of KerrGeoSeparatrix[a,p,e,x] for each orbit in separatrix_values.txt test_frequencies.py: -freq_values.txt - orbital parameters (a,p,e,x) +freq_parameters.txt - orbital parameters (a,p,e,x) mathematica_freq_output.txt - (upsilon_r, upsilon_theta, upsilon_phi, gamma) computed using KerrGeoFrequencies[a,p,e,x] for each orbit in freq_values.txt test_stable_solutions.py: -stable_orbit_values.txt - orbital parameters (a,p,e,x) +stable_orbit_parameters.txt - orbital parameters (a,p,e,x) stable_orbit_times.txt - list of mino time values to test stable_orbits/trajectory{i}.txt - (t, r, theta, phi) evaluated at each time from stable_orbit_times.txt for the i-th orbit defined in stable_orbit_values.txt stable_solutions/trajectory{i}.txt - (t_r, t_theta, phi_r, phi_theta) evaluated at each time from stable_orbit_times.txt for the i-th orbit defined in stable_orbit_values.txt @@ -28,4 +29,12 @@ plunging_orbit_times.txt - list of mino time values to test plunging_integrals/trajectory{i}.txt - (I_r, I_r2, I_r_plus, I_r_minus) evaluated at each time from plunging_orbit_times.txt for the i-th orbit defined in plunging_orbit_parameters_complex.txt plunging_solutions/trajectory{i}.txt - (t_r, phi_r, t_theta, phi_theta) evaluated at each time from plunging_orbit_times.txt for the i-th orbit defined in plunging_orbit_parameters_complex.txt plunging_orbits_real/trajectory{i}.txt - (t, r, theta, phi) evaluated at each time from plunging_orbit_times.txt for the i-th orbit defined in plunging_orbit_parameters_real.txt -plunging_orbits_complex/trajectory{i}.txt - (t, r, theta, phi) evaluated at each time from plunging_orbit_times.txt for the i-th orbit defined in plunging_orbit_parameters_complex.txt \ No newline at end of file +plunging_orbits_complex/trajectory{i}.txt - (t, r, theta, phi) evaluated at each time from plunging_orbit_times.txt for the i-th orbit defined in plunging_orbit_parameters_complex.txt + +test_four_velocity.py: + +stable_orbit_parameters.txt - orbital parameters (a,p,e,x) +stable_orbit_times.txt - list of mino time values to test +plunging_orbit_parameters_real.txt - list of orbital parameters (a,E,L,Q) for which the radial polynomial has all real roots +plunging_orbit_parameters_complex.txt - list of orbital parameters (a,E,L,Q) for which the radial polynomial has complex roots +four_velocity/trajectory{i}.txt - (u_t, u_r, u_theta, u_phi) evaluated at each time from stable_orbit_times.txt for the i-th orbit in stable_orbit_values.txt \ No newline at end of file diff --git a/tests/data/Unit Test Generation.ipynb b/tests/data/Unit Test Generation.ipynb new file mode 100644 index 0000000..8b735e8 --- /dev/null +++ b/tests/data/Unit Test Generation.ipynb @@ -0,0 +1,417 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 4, + "id": "8727d99d", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from itertools import product\n", + "import matplotlib.pyplot as plt\n", + "from re import sub\n", + "from wolframclient.evaluation import WolframLanguageSession\n", + "from wolframclient.language import wl, wlexpr\n", + "from kerrgeopy.constants import *\n", + "import kerrgeopy\n", + "from numpy.polynomial import Polynomial\n", + "\n", + "session = WolframLanguageSession()\n", + "session.evaluate(wlexpr(\"PacletDirectoryLoad[\\\"~/Documents/Wolfram Mathematica/BHPToolkit/KerrGeodesics\\\"]\"))\n", + "session.evaluate(wl.Needs('KerrGeodesics`'))" + ] + }, + { + "cell_type": "markdown", + "id": "9b0d048f", + "metadata": {}, + "source": [ + "## Separatrix" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "56234f75", + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "num_values = int(1000)\n", + "sep_values = np.zeros((num_values,3))\n", + "i = 0\n", + "while i < num_values:\n", + " parameters = np.round([np.random.rand(), np.random.rand(), np.random.rand()*2-1], \\\n", + " decimals = 1)\n", + " if parameters[0] != 1:\n", + " sep_values[i] = parameters\n", + " i += 1\n", + " \n", + "def sub_zeros(txt):\n", + " txt1 = sub(\"0\\.0*,\",\"0,\",txt)\n", + " txt1 = sub(\"0\\.0*]\",\"0]\",txt1)\n", + " txt1 = sub(\"1\\.0*,\",\"1,\",txt1)\n", + " txt1 = sub(\"1\\.0*]\",\"1]\",txt1)\n", + " return txt1\n", + "\n", + "mathematica_sep_output = np.apply_along_axis(lambda x: \\\n", + " session.evaluate(wlexpr(\n", + " sub_zeros(\"KerrGeoSeparatrix[{:},{:},{:}]\".format(*x)))) \\\n", + " ,1,sep_values)\n", + "\n", + "np.savetxt(\"sep_values.txt\",sep_values,fmt=\"%.3f, %.3f, %.3f\")\n", + "np.savetxt(\"mathematica_sep_output.txt\",mathematica_sep_output)" + ] + }, + { + "cell_type": "markdown", + "id": "4e6434c0", + "metadata": {}, + "source": [ + "## Constants" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "544b5fbe", + "metadata": {}, + "outputs": [], + "source": [ + "a_values = [0,0.5]\n", + "p_values = [12]\n", + "e_values = [0, 0.5, 1]\n", + "x_values = [-1,-0.5,0,0.5,1]\n", + "\n", + "const_values = np.array(list(product(a_values,p_values,e_values,x_values)))\n", + "\n", + "num_random_values = int(100)\n", + "random_values = np.zeros((num_random_values,4))\n", + "i = 0\n", + "while i < num_random_values:\n", + " parameters = np.round([np.random.rand(), np.random.rand()*20, np.random.rand(), np.random.rand()*2-1], \\\n", + " decimals = 3)\n", + " if parameters[0] != 1 and parameters[1] > separatrix(parameters[0],parameters[2],parameters[3]):\n", + " random_values[i] = parameters\n", + " i += 1\n", + " \n", + "const_values = np.concatenate([const_values,random_values])\n", + "\n", + "mathematica_const_output = np.apply_along_axis(lambda x: \\\n", + " session.evaluate(\"Values[KerrGeoConstantsOfMotion[{:},{:},{:},{:}]]\".format(*x)) \\\n", + " ,1,const_values)\n", + "\n", + "np.savetxt(\"const_parameters.txt\",const_values,fmt=\"%.3f, %.3f, %.3f, %.3f\")\n", + "np.savetxt(\"mathematica_const_output.txt\",mathematica_const_output)" + ] + }, + { + "cell_type": "markdown", + "id": "53dac606", + "metadata": {}, + "source": [ + "## Frequencies" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "id": "34830563", + "metadata": {}, + "outputs": [], + "source": [ + "a_values = [0,0.5]\n", + "p_values = [12]\n", + "e_values = [0, 0.5]\n", + "x_values = [-1,-0.5,0.5,1]\n", + "\n", + "freq_values = np.array(list(product(a_values,p_values,e_values,x_values)))\n", + "\n", + "num_random_values = int(100)\n", + "random_values = np.zeros((num_random_values,4))\n", + "i = 0\n", + "while i < num_random_values:\n", + " parameters = np.round([np.random.rand(), np.random.rand()*20, np.random.rand(), np.random.rand()*2-1], \\\n", + " decimals = 3)\n", + " if parameters[0] != 1 and parameters[1] > separatrix(parameters[0],parameters[2],parameters[3]):\n", + " random_values[i] = parameters\n", + " i += 1\n", + " \n", + "freq_values = np.concatenate([freq_values,random_values])\n", + "\n", + "mathematica_freq_output = np.apply_along_axis(lambda x: \\\n", + " session.evaluate(\n", + " wlexpr(\"Values[KerrGeoFrequencies[{:},{:},{:},{:},\\\"Time\\\" -> \\\"Mino\\\"]]\".format(*x))) \\\n", + " ,1,freq_values)\n", + "\n", + "\n", + "np.savetxt(\"freq_parameters.txt\",freq_values,fmt=\"%.3f, %.3f, %.3f, %.3f\")\n", + "np.savetxt(\"mathematica_freq_output.txt\",mathematica_freq_output)" + ] + }, + { + "cell_type": "markdown", + "id": "fca4e126", + "metadata": {}, + "source": [ + "## Stable Solutions" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "efc1959f", + "metadata": {}, + "outputs": [], + "source": [ + "a_values = [0,0.5]\n", + "p_values = [12]\n", + "e_values = [0, 0.5]\n", + "x_values = [-1,-0.5,0.5,1]\n", + "\n", + "orbit_values = np.array(list(product(a_values,p_values,e_values,x_values)))\n", + "\n", + "num_random_values = int(50)\n", + "random_values = np.zeros((num_random_values,4))\n", + "i = 0\n", + "while i < num_random_values:\n", + " parameters = np.round([np.random.rand(), np.random.rand()*20, np.random.rand(), np.random.rand()*2-1], \\\n", + " decimals = 3)\n", + " if parameters[0] != 1 and parameters[2] != 1 and parameters[3] != 0 and parameters[1] > separatrix(parameters[0],parameters[2],parameters[3]):\n", + " random_values[i] = parameters\n", + " i += 1\n", + " \n", + "orbit_values = np.concatenate([orbit_values,random_values])\n", + "np.savetxt(\"stable_orbit_parameters.txt\",orbit_values,fmt=\"%.3f, %.3f, %.3f, %.3f\")\n", + "times = np.array(range(10))\n", + "np.savetxt(\"stable_orbit_times.txt\",times)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b73af1e3", + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "for i, orbit in enumerate(orbit_values):\n", + " session.evaluate(wlexpr(\"orbit = KerrGeoOrbitPhases[{},{},{},{}]\".format(*orbit)))\n", + " session.evaluate(wlexpr(\"{tr, ttheta, phir, phitheta} = Values[orbit[\\\"TrajectoryDeltas\\\"]];\"))\n", + " trajectory_deltas = np.zeros((10,4))\n", + " for j,t in enumerate(times):\n", + " trajectory_deltas[j] = np.array(session.evaluate(wlexpr(f'N[{{tr[{t}],ttheta[{t}],phir[{t}],phitheta[{t}]}}]')))\n", + " np.savetxt(\"stable_solutions/trajectory\"+str(i)+\".txt\",trajectory_deltas,delimiter=\",\")" + ] + }, + { + "cell_type": "markdown", + "id": "ea7cbfe2", + "metadata": {}, + "source": [ + "## Stable Orbit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04a14df1", + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "for i, orbit in enumerate(orbit_values):\n", + " session.evaluate(wlexpr(\"orbit = KerrGeoOrbit[{},{},{},{}]\".format(*orbit)))\n", + " session.evaluate(wlexpr(\"{t, r, theta, phi} = orbit[\\\"Trajectory\\\"]\"))\n", + " trajectory = np.zeros((10,4))\n", + " for j,t in enumerate(times):\n", + " trajectory[j] = np.array(session.evaluate(wlexpr(f'N[{{t[{t}],r[{t}],theta[{t}],phi[{t}]}}]')))\n", + " np.savetxt(\"stable_orbits/trajectory\"+str(i)+\".txt\",trajectory,delimiter=\",\")" + ] + }, + { + "cell_type": "markdown", + "id": "5345f2f0", + "metadata": {}, + "source": [ + "## Plunging Solutions and Integrals (Complex Roots)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "1db59f20", + "metadata": {}, + "outputs": [], + "source": [ + "num_random_values = int(100)\n", + "orbit_values = np.zeros((num_random_values,4))\n", + "i = 0\n", + "\n", + "while i < num_random_values:\n", + " parameters = np.round([np.random.rand(), np.random.rand(), (np.random.rand()-0.5)*10, np.random.rand()*10], \\\n", + " decimals = 3)\n", + " # Filter only parameters with complex roots\n", + " a, E, L, Q = parameters\n", + " R = Polynomial([-a**2*Q, 2*L**2+2*Q+2*a**2*E**2-4*a*E*L, a**2*E**2-L**2-Q-a**2, 2, E**2-1])\n", + " radial_roots = R.roots()\n", + " complex_roots = radial_roots[np.iscomplex(radial_roots)]\n", + " \n", + " if parameters[0] != 1 and parameters[1] != 1 and len(complex_roots) > 0:\n", + " orbit_values[i] = parameters\n", + " i += 1\n", + "\n", + "np.savetxt(\"plunging_orbit_parameters_complex.txt\",orbit_values,fmt=\"%.3f, %.3f, %.3f, %.3f\")\n", + "times = np.array(range(10))\n", + "np.savetxt(\"plunging_orbit_times.txt\",times)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "12bd08bd", + "metadata": {}, + "outputs": [], + "source": [ + "for i, orbit in enumerate(orbit_values):\n", + " session.evaluate(wlexpr(\"orbit = KerrGeoPlunge[{},{{{},{},{}}}]\".format(*orbit)))\n", + " session.evaluate(wlexpr(\"{Ir, Ir2, Irp, Irm} = Values[orbit[\\\"RadialIntegrals\\\"]]\"))\n", + " session.evaluate(wlexpr(\"{tr,phir,tz,phiz} = Values[orbit[\\\"TrajectoryDeltas\\\"]]\"))\n", + " trajectory_deltas = np.zeros((10,4))\n", + " integrals = np.zeros((10,4))\n", + " for j,t in enumerate(times):\n", + " trajectory_deltas[j] = np.array(session.evaluate(wlexpr(f'N[{{tr[{t}],phir[{t}],tz[{t}],phiz[{t}]}}]')))\n", + " integrals[j] = np.array(session.evaluate(wlexpr(f'N[{{Ir[{t}],Ir2[{t}],Irp[{t}],Irm[{t}]}}]')))\n", + " np.savetxt(\"plunging_solutions/trajectory\"+str(i)+\".txt\",trajectory_deltas,delimiter=\",\")\n", + " np.savetxt(\"plunging_integrals/trajectory\"+str(i)+\".txt\",integrals,delimiter=\",\")" + ] + }, + { + "cell_type": "markdown", + "id": "0c32215e", + "metadata": {}, + "source": [ + "## Plunging Orbit (Complex Roots)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "0ad8225a", + "metadata": {}, + "outputs": [], + "source": [ + "for i, orbit in enumerate(orbit_values):\n", + " session.evaluate(wlexpr(\"orbit = KerrGeoPlunge[{},{{{},{},{}}}]\".format(*orbit)))\n", + " session.evaluate(wlexpr(\"{t, r, theta, phi} = orbit[\\\"Trajectory\\\"]\"))\n", + " trajectory = np.zeros((len(times),4))\n", + " for j,t in enumerate(times):\n", + " trajectory[j] = np.array(session.evaluate(wlexpr(f'N[{{t[{t}],r[{t}],theta[{t}],phi[{t}]}}]')))\n", + " np.savetxt(\"plunging_orbits_complex/trajectory\"+str(i)+\".txt\",trajectory,delimiter=\",\")" + ] + }, + { + "cell_type": "markdown", + "id": "e971bbb7", + "metadata": {}, + "source": [ + "## Plunging Orbit (Real)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "170bad0e", + "metadata": {}, + "outputs": [], + "source": [ + "num_random_values = int(100)\n", + "orbit_values = np.zeros((num_random_values,4))\n", + "i = 0\n", + "\n", + "while i < num_random_values:\n", + " parameters = np.round([np.random.rand(), np.random.rand(), (np.random.rand()-0.5)*10, np.random.rand()*10], \\\n", + " decimals = 3)\n", + " a, E, L, Q = parameters\n", + " R = Polynomial([-a**2*Q, 2*L**2+2*Q+2*a**2*E**2-4*a*E*L, a**2*E**2-L**2-Q-a**2, 2, E**2-1])\n", + " radial_roots = R.roots()\n", + " complex_roots = radial_roots[np.iscomplex(radial_roots)]\n", + " \n", + " if parameters[0] != 1 and parameters[1] != 1 and len(complex_roots) == 0:\n", + " orbit_values[i] = parameters\n", + " i += 1\n", + "\n", + "np.savetxt(\"plunging_orbit_parameters_real.txt\",orbit_values,fmt=\"%.3f, %.3f, %.3f, %.3f\")" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "888af26c", + "metadata": {}, + "outputs": [], + "source": [ + "for i, orbit in enumerate(orbit_values):\n", + " session.evaluate(wlexpr(\"orbit = KerrGeoPlunge[{},{{{},{},{}}}]\".format(*orbit)))\n", + " session.evaluate(wlexpr(\"{t, r, theta, phi} = orbit[\\\"Trajectory\\\"]\"))\n", + " trajectory = np.zeros((len(times),4))\n", + " for j,t in enumerate(times):\n", + " trajectory[j] = np.array(session.evaluate(wlexpr(f'N[{{t[{t}],r[{t}],theta[{t}],phi[{t}]}}]')))\n", + " np.savetxt(\"plunging_orbits_real/trajectory\"+str(i)+\".txt\",trajectory,delimiter=\",\")" + ] + }, + { + "cell_type": "markdown", + "id": "74cf295a", + "metadata": {}, + "source": [ + "## Four Velocity" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "f0f6d0f8", + "metadata": {}, + "outputs": [], + "source": [ + "orbit_values = np.genfromtxt(\"/users/spark59/Documents/KerrGeoPy/tests/data/stable_orbit_values.txt\",\n", + " delimiter=\",\")\n", + "for i, orbit in enumerate(orbit_values):\n", + " session.evaluate(wlexpr(\"u = KerrGeoFourVelocity[{:},{:},{:},{:}];\".format(*orbit)))\n", + " session.evaluate(wlexpr(\"{ut,ur,utheta,uphi} = Values[u]\"))\n", + " trajectory = np.zeros((10,4))\n", + " for j,t in enumerate(times):\n", + " trajectory[j] = np.array(\n", + " session.evaluate(wlexpr(f'N[{{ut[{t}],ur[{t}],utheta[{t}],uphi[{t}]}}]'))\n", + " )\n", + " np.savetxt(\"four_velocity/trajectory\"+str(i)+\".txt\",trajectory,delimiter=\",\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tests/data/const_values.txt b/tests/data/const_parameters.txt similarity index 100% rename from tests/data/const_values.txt rename to tests/data/const_parameters.txt diff --git a/tests/data/freq_values.txt b/tests/data/freq_parameters.txt similarity index 100% rename from tests/data/freq_values.txt rename to tests/data/freq_parameters.txt diff --git a/tests/data/separatrix_values.txt b/tests/data/separatrix_parameters.txt similarity index 100% rename from tests/data/separatrix_values.txt rename to tests/data/separatrix_parameters.txt diff --git a/tests/data/stable_orbit_values.txt b/tests/data/stable_orbit_parameters.txt similarity index 100% rename from tests/data/stable_orbit_values.txt rename to tests/data/stable_orbit_parameters.txt diff --git a/tests/test_constants.py b/tests/test_constants.py index 2048997..915678a 100644 --- a/tests/test_constants.py +++ b/tests/test_constants.py @@ -27,7 +27,7 @@ def test_invalid_arguments(self): with self.assertRaises(ValueError): carter_constant(-0.5,5,-0.5,0.5) def test_constants_random(self): - values = np.genfromtxt(DATA_DIR / "const_values.txt",delimiter=",") + values = np.genfromtxt(DATA_DIR / "const_parameters.txt",delimiter=",") mathematica_const_output = np.genfromtxt(DATA_DIR / "mathematica_const_output.txt") python_const_output = np.apply_along_axis(lambda x: constants_of_motion(*x),1,values) @@ -44,7 +44,7 @@ def test_invalid_arguments(self): with self.assertRaises(ValueError): separatrix(2,-0.5,-0.5) def test_separatrix_random(self): - sep_values = np.genfromtxt(DATA_DIR / "separatrix_values.txt",delimiter=",") + sep_values = np.genfromtxt(DATA_DIR / "separatrix_parameters.txt",delimiter=",") mathematica_separatrix_output = np.genfromtxt(DATA_DIR / "mathematica_separatrix_output.txt") python_separatrix_output = np.apply_along_axis(lambda x: separatrix(*x),1,sep_values) for i, params in enumerate(sep_values): diff --git a/tests/test_four_velocity.py b/tests/test_four_velocity.py index a1c136c..821fe08 100644 --- a/tests/test_four_velocity.py +++ b/tests/test_four_velocity.py @@ -1,14 +1,14 @@ import unittest import numpy as np -from kerrgeopy.stable_orbit import * -from kerrgeopy.plunging_orbit import * +from kerrgeopy.stable import * +from kerrgeopy.plunge import * from pathlib import Path THIS_DIR = Path(__file__).parent DATA_DIR = THIS_DIR.parent / "tests/data" -stable_orbit_values = np.genfromtxt(DATA_DIR / "stable_orbit_values.txt", delimiter=",") +stable_orbit_values = np.genfromtxt(DATA_DIR / "stable_orbit_parameters.txt", delimiter=",") complex_plunging_orbit_values = np.genfromtxt(DATA_DIR / "plunging_orbit_parameters_complex.txt", delimiter=",") real_plunging_orbit_values = np.genfromtxt(DATA_DIR / "plunging_orbit_parameters_real.txt", delimiter=",") plunging_orbit_values = np.concatenate((complex_plunging_orbit_values,real_plunging_orbit_values),axis=0) diff --git a/tests/test_frequencies.py b/tests/test_frequencies.py index 282e267..d73f9f7 100644 --- a/tests/test_frequencies.py +++ b/tests/test_frequencies.py @@ -10,36 +10,21 @@ class TestFrequencies(unittest.TestCase): def test_extreme_kerr(self): with self.assertRaises(ValueError): mino_frequencies(1,12,0.5,0.5) - with self.assertRaises(ValueError): r_frequency(1,0.5,12,0.5) - with self.assertRaises(ValueError): theta_frequency(1,0.5,12,0.5) - with self.assertRaises(ValueError): phi_frequency(1,0.5,12,0.5) def test_polar(self): with self.assertRaises(ValueError): mino_frequencies(0.5,12,0.5,0) - with self.assertRaises(ValueError): r_frequency(0.5,0.5,12,0) - with self.assertRaises(ValueError): theta_frequency(0.5,0.5,12,0) - with self.assertRaises(ValueError): phi_frequency(0.5,0.5,12,0) def test_marginally_bound(self): with self.assertRaises(ValueError): mino_frequencies(0.5,12,1,0.5) - with self.assertRaises(ValueError): r_frequency(0.5,12,1,0.5) - with self.assertRaises(ValueError): theta_frequency(0.5,12,1,0.5) - with self.assertRaises(ValueError): phi_frequency(0.5,12,1,0.5) def test_invalid_arguments(self): with self.assertRaises(ValueError): mino_frequencies(2,5,0.5,0.5) - with self.assertRaises(ValueError): r_frequency(0.5,5,-0.5,0.5) - with self.assertRaises(ValueError): theta_frequency(0.5,5,0.5,-2) - with self.assertRaises(ValueError): phi_frequency(0.5,5,2,0.5) def test_unstable(self): with self.assertRaises(ValueError): mino_frequencies(0.5,5,0.5,0.5) - with self.assertRaises(ValueError): r_frequency(0.5,5,0.5,0.5) - with self.assertRaises(ValueError): theta_frequency(0.5,5,0.5,0.5) - with self.assertRaises(ValueError): phi_frequency(0.5,5,0.5,0.5) def test_random(self): - values = np.genfromtxt(DATA_DIR / "freq_values.txt",delimiter=",") + values = np.genfromtxt(DATA_DIR / "freq_parameters.txt",delimiter=",") mathematica_freq_output = np.genfromtxt(DATA_DIR / "mathematica_freq_output.txt") python_freq_output = np.apply_along_axis(lambda x: mino_frequencies(*x),1,values) for i, params in enumerate(values): diff --git a/tests/test_initial_conditions.py b/tests/test_initial_conditions.py index beec823..eeabe96 100644 --- a/tests/test_initial_conditions.py +++ b/tests/test_initial_conditions.py @@ -7,7 +7,7 @@ DATA_DIR = THIS_DIR.parent / "tests/data" -stable_orbit_values = np.genfromtxt(DATA_DIR / "stable_orbit_values.txt", delimiter=",") +stable_orbit_values = np.genfromtxt(DATA_DIR / "stable_orbit_parameters.txt", delimiter=",") complex_plunging_orbit_values = np.genfromtxt(DATA_DIR / "plunging_orbit_parameters_complex.txt", delimiter=",") real_plunging_orbit_values = np.genfromtxt(DATA_DIR / "plunging_orbit_parameters_real.txt", delimiter=",") diff --git a/tests/test_plunging_solutions.py b/tests/test_plunging_solutions.py index cd8927e..a225b29 100644 --- a/tests/test_plunging_solutions.py +++ b/tests/test_plunging_solutions.py @@ -1,6 +1,6 @@ import unittest import numpy as np -from kerrgeopy.plunging_orbit import * +from kerrgeopy.plunge import * from pathlib import Path THIS_DIR = Path(__file__).parent diff --git a/tests/test_stable_solutions.py b/tests/test_stable_solutions.py index 78dc992..fdddb39 100644 --- a/tests/test_stable_solutions.py +++ b/tests/test_stable_solutions.py @@ -1,7 +1,7 @@ import unittest import numpy as np -from kerrgeopy.stable_orbit import * -from kerrgeopy.frequencies import _radial_roots, _polar_roots +from kerrgeopy.stable import * +from kerrgeopy.constants import stable_radial_roots, stable_polar_roots from pathlib import Path THIS_DIR = Path(__file__).parent @@ -9,7 +9,7 @@ DATA_DIR = THIS_DIR.parent / "tests/data" -orbit_values = np.genfromtxt(DATA_DIR / "stable_orbit_values.txt", delimiter=",") +orbit_values = np.genfromtxt(DATA_DIR / "stable_orbit_parameters.txt", delimiter=",") orbit_times = np.genfromtxt(DATA_DIR / "stable_orbit_times.txt", delimiter=",") class TestStableSolutions(unittest.TestCase): @@ -21,8 +21,8 @@ def test_random(self): a,p,e,x = orbit constants = constants_of_motion(*orbit) - radial_roots = _radial_roots(a,p,e,constants) - polar_roots = _polar_roots(a,x,constants) + radial_roots = stable_radial_roots(a,p,e,x,constants) + polar_roots = stable_polar_roots(a,p,e,x,constants) r, t_r, phi_r = radial_solutions(a,constants,radial_roots) theta, t_theta, phi_theta = polar_solutions(a,constants,polar_roots)