Merge branch 'main' into paper

README.md

@@ -1,6 +1,7 @@
 [![GitHub release (with filter)](https://img.shields.io/github/v/release/BlackHolePerturbationToolkit/KerrGeoPy)](https://github.com/BlackHolePerturbationToolkit/KerrGeoPy/releases)
 [![Test Status](https://github.com/BlackHolePerturbationToolkit/KerrGeoPy/actions/workflows/tests.yml/badge.svg)](https://github.com/BlackHolePerturbationToolkit/KerrGeoPy/actions)
 [![PyPI - Version](https://img.shields.io/pypi/v/kerrgeopy)](https://pypi.org/project/kerrgeopy/)
+[![Conda Version](https://img.shields.io/conda/vn/conda-forge/kerrgeopy.svg)](https://anaconda.org/conda-forge/kerrgeopy)
 [![Documentation Status](https://readthedocs.org/projects/kerrgeopy/badge/?version=latest)](https://kerrgeopy.readthedocs.io/en/latest/?badge=latest)
 [![GitHub License](https://img.shields.io/github/license/BlackHolePerturbationToolkit/KerrGeoPy)](https://github.com/BlackHolePerturbationToolkit/KerrGeoPy/blob/main/LICENSE)
 [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.10126824.svg)](https://doi.org/10.5281/zenodo.10126824)
@@ -11,7 +12,7 @@ KerrGeoPy is a python implementation of the [KerrGeodesics](https://bhptoolkit.o
 
 ## Installation
 
-Install using anaconda
+Install using Anaconda
 ```bash
 conda install -c conda-forge kerrgeopy
 ```
@@ -22,11 +23,11 @@ pip install kerrgeopy
 
 > **Note**
 >
-> KerrGeoPy uses functions introduced in scipy 1.8, so it may also be necessary to update scipy by running `pip install scipy -U`, although in most cases this should be done automatically by pip. Certain plotting and animation functions also make use of features introduced in matplotlib 3.7 and rely on [ffmpeg](https://ffmpeg.org/download.html), which can be easily installed using [homebrew](https://formulae.brew.sh/formula/ffmpeg) or [anaconda](https://anaconda.org/conda-forge/ffmpeg).
+> This library uses functions introduced in scipy 1.8, so it may also be necessary to update scipy by running `pip install scipy -U`, although in most cases this should be done automatically by pip. Certain plotting and animation functions also make use of features introduced in matplotlib 3.7 and rely on [ffmpeg](https://ffmpeg.org/download.html), which can be easily installed using [homebrew](https://formulae.brew.sh/formula/ffmpeg) or [anaconda](https://anaconda.org/conda-forge/ffmpeg).
 
 ## Stable Bound Orbits
 
-Kerrgeopy computes orbits in Boyer-Lindquist coordinates $(t,r,\theta,\phi)$. Let $M$ to represent the mass of the primary body and let $J$ represent its angular momentum. Working in geometrized units where $G=c=1$, stable bound orbits are parametrized using the following variables:
+KerrGeoPy computes orbits in Boyer-Lindquist coordinates $(t,r,\theta,\phi)$. Let $M$ to represent the mass of the primary body and let $J$ represent its angular momentum. Working in geometrized units where $G=c=1$, stable bound orbits are parametrized using the following variables:
 
 $a$ - spin of the primary body
 <br>

docs/source/index.rst

@@ -22,7 +22,7 @@ The library also provides a set of methods for computing constants of motion and
 
 Installation
 ------------
-Install using anaconda
+Install using Anaconda
 
 .. code-block:: bash
 

tests/test_constants.py

@@ -50,15 +50,22 @@ def test_invalid_arguments(self):
 
     def test_constants_random(self):
         """
-        Test constants_of_motion method against output from the KerrGeodesics Mathematica library for a random set of orbits.
+        Test constants_of_motion method against output from the KerrGeodesics
+        Mathematica library for a random set of orbits.
         """
-        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)
+        parameters = 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, parameters
+        )
 
-        for i, params in enumerate(values):
+        for i, params in enumerate(parameters):
             with self.subTest(i=i, params=params):
-                self.assertTrue(np.allclose(mathematica_const_output[i], python_const_output[i]))
+                self.assertTrue(
+                    np.allclose(mathematica_const_output[i], python_const_output[i])
+                )
 
 
 class TestSeparatrix(unittest.TestCase):
@@ -78,15 +85,22 @@ def test_invalid_arguments(self):
 
     def test_separatrix_random(self):
         """
-        Test separatrix method against output from the KerrGeodesics Mathematica library for a random set of inputs.
+        Test separatrix method against output from the KerrGeodesics
+        Mathematica library for a random set of inputs.
         """
-        sep_values = np.genfromtxt(DATA_DIR / "separatrix_parameters.txt", delimiter=",")
+        sep_parameters = 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):
+        python_separatrix_output = np.apply_along_axis(
+            lambda x: separatrix(*x), 1, sep_parameters
+        )
+        for i, params in enumerate(sep_parameters):
             with self.subTest(i=i, params=params):
                 self.assertTrue(
-                    np.allclose(python_separatrix_output[i], mathematica_separatrix_output[i])
+                    np.allclose(
+                        python_separatrix_output[i], mathematica_separatrix_output[i]
+                    )
                 )

tests/test_four_velocity.py

@@ -8,27 +8,30 @@
 
 DATA_DIR = THIS_DIR.parent / "tests/data"
 
-stable_orbit_values = np.genfromtxt(DATA_DIR / "stable_orbit_parameters.txt", delimiter=",")
-complex_plunging_orbit_values = np.genfromtxt(
+stable_orbit_parameters = np.genfromtxt(
+    DATA_DIR / "stable_orbit_parameters.txt", delimiter=","
+)
+complex_plunging_orbit_parameters = np.genfromtxt(
     DATA_DIR / "plunging_orbit_parameters_complex.txt", delimiter=","
 )
-real_plunging_orbit_values = np.genfromtxt(
+real_plunging_orbit_parameters = 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
+plunging_orbit_parameters = np.concatenate(
+    (complex_plunging_orbit_parameters, real_plunging_orbit_parameters), axis=0
 )
 times = np.genfromtxt(DATA_DIR / "stable_orbit_times.txt", delimiter=",")
 
 
 class TestFourVelocity(unittest.TestCase):
     def test_stable_orbit_four_velocity(self):
         """
-        Test four_velocity method against output from the KerrGeodesics Mathematica library for a random set of orbits.
+        Test four_velocity method against output from the KerrGeodesics
+        Mathematica library for a random set of orbits.
         """
         components = ["ut", "ur", "utheta", "uphi"]
 
-        for i, orbit in enumerate(stable_orbit_values):
+        for i, orbit in enumerate(stable_orbit_parameters):
             mathematica_trajectory = np.genfromtxt(
                 DATA_DIR / f"four_velocity/trajectory{i}.txt", delimiter=","
             )
@@ -48,11 +51,15 @@ def test_stable_orbit_four_velocity(self):
                     i=i,
                     component=component,
                     params="a = {}, p = {}, e = {}, x = {}".format(*orbit),
-                    diff=np.max(np.abs(mathematica_trajectory[:, j] - python_trajectory[:, j])),
+                    diff=np.max(
+                        np.abs(mathematica_trajectory[:, j] - python_trajectory[:, j])
+                    ),
                 ):
                     self.assertTrue(
                         np.allclose(
-                            mathematica_trajectory[:, j], python_trajectory[:, j], atol=1e-6
+                            mathematica_trajectory[:, j],
+                            python_trajectory[:, j],
+                            atol=1e-6,
                         )
                     )
 
@@ -62,26 +69,30 @@ def test_norm(self):
         """
         times = np.linspace(0, 10, 10)
         # stable orbits
-        for i, orbit in enumerate(stable_orbit_values):
+        for i, orbit in enumerate(stable_orbit_parameters):
             a, p, e, x = orbit
             stable_orbit = StableOrbit(a, p, e, x)
             norm = stable_orbit._four_velocity_norm()
-            for j, time in enumerate(times):
+            for time in times:
                 with self.subTest(
-                    i=i, params="a = {}, p = {}, e = {}, x = {}".format(*orbit), norm=norm(time)
+                    i=i,
+                    params="a = {}, p = {}, e = {}, x = {}".format(*orbit),
+                    norm=norm(time),
                 ):
                     self.assertTrue(abs(norm(time) + 1) < 1e-8)
 
         # plunging orbits
-        for i, orbit in enumerate(plunging_orbit_values):
+        for i, orbit in enumerate(plunging_orbit_parameters):
             a, E, L, Q = orbit
             plunging_orbit = PlungingOrbit(a, E, L, Q)
             norm = plunging_orbit._four_velocity_norm()
             for j, time in enumerate(times):
                 with self.subTest(
-                    i=i, params="a = {}, E = {}, L = {}, Q = {}".format(*orbit), norm=norm(time)
+                    i=i,
+                    params="a = {}, E = {}, L = {}, Q = {}".format(*orbit),
+                    norm=norm(time),
                 ):
-                    self.assertTrue(abs(norm(time) + 1) < 1e-5)
+                    self.assertTrue(abs(norm(time) + 1) < 1e-3)
 
     def test_using_numerical_differentiation(self):
         """
@@ -91,11 +102,16 @@ def test_using_numerical_differentiation(self):
         components = ["ut", "ur", "utheta", "uphi"]
 
         # stable orbits
-        for i, orbit in enumerate(stable_orbit_values):
+        for i, orbit in enumerate(stable_orbit_parameters):
             a, p, e, x = orbit
             stable_orbit = StableOrbit(a, p, e, x)
             u_t, u_r, u_theta, u_phi = stable_orbit.four_velocity()
-            delta_t, delta_r, delta_theta, delta_phi = stable_orbit.numerical_four_velocity()
+            (
+                delta_t,
+                delta_r,
+                delta_theta,
+                delta_phi,
+            ) = stable_orbit.numerical_four_velocity()
 
             analytic_four_velocity = np.transpose(
                 np.apply_along_axis(
@@ -105,7 +121,9 @@ def test_using_numerical_differentiation(self):
 
             numerical_four_velocity = np.transpose(
                 np.apply_along_axis(
-                    lambda x: np.array([delta_t(x), delta_r(x), delta_theta(x), delta_phi(x)]),
+                    lambda x: np.array(
+                        [delta_t(x), delta_r(x), delta_theta(x), delta_phi(x)]
+                    ),
                     0,
                     times,
                 )
@@ -117,51 +135,15 @@ def test_using_numerical_differentiation(self):
                     component=component,
                     params="a = {}, p = {}, e = {}, x = {}".format(*orbit),
                     diff=np.max(
-                        np.abs(analytic_four_velocity[:, j] - numerical_four_velocity[:, j])
+                        np.abs(
+                            analytic_four_velocity[:, j] - numerical_four_velocity[:, j]
+                        )
                     ),
                 ):
                     self.assertTrue(
                         np.allclose(
-                            analytic_four_velocity[:, j], numerical_four_velocity[:, j], atol=1e-3
+                            analytic_four_velocity[:, j],
+                            numerical_four_velocity[:, j],
+                            atol=1e-3,
                         )
-                    )
-
-        # plunging orbits
-
-        # for i, orbit in enumerate(plunging_orbit_values):
-        #     a, E, L, Q = orbit
-        #     plunging_orbit = PlungingOrbit(a, E, L, Q)
-        #     u_t, u_r, u_theta, u_phi = plunging_orbit.four_velocity()
-        #     delta_t, delta_r, delta_theta, delta_phi = plunging_orbit.numerical_four_velocity()
-
-        #     analytic_four_velocity = np.transpose(
-        #         np.apply_along_axis(
-        #             lambda x: np.array([u_t(x), u_r(x), u_theta(x), u_phi(x)]), 0, times
-        #         )
-        #     )
-
-        #     numerical_four_velocity = np.transpose(
-        #         np.apply_along_axis(
-        #             lambda x: np.array([delta_t(x), delta_r(x), delta_theta(x), delta_phi(x)]),
-        #             0,
-        #             times,
-        #         )
-        #     )
-
-        #     for j, component in enumerate(components):
-        #         with self.subTest(
-        #             i=i,
-        #             component=component,
-        #             params="a = {}, E = {}, L = {}, Q = {}".format(*orbit),
-        #             diff=np.max(
-        #                 np.abs(analytic_four_velocity[:, j] - numerical_four_velocity[:, j])
-        #             ),
-        #         ):
-        #             self.assertTrue(
-        #                 np.allclose(
-        #                     analytic_four_velocity[:, j],
-        #                     numerical_four_velocity[:, j],
-        #                     rtol=1e-3,
-        #                     atol=1e-1,
-        #                 )
-        #             )
+                    )

tests/test_frequencies.py

@@ -46,14 +46,20 @@ def test_unstable(self):
 
     def test_random(self):
         """
-        Test mino_frequencies method against output from the KerrGeodesics Mathematica library for a random set of orbits.
+        Test mino_frequencies method against output from the KerrGeodesics 
+        Mathematica library for a random set of orbits.
         """
-        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):
+        parameters = 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, parameters
+        )
+        for i, params in enumerate(parameters):
             with self.subTest(i=i, params=params):
-                E, L, Q = constants_of_motion(*params)
                 self.assertTrue(
-                    np.allclose(abs(mathematica_freq_output[i]), abs(python_freq_output[i]))
+                    np.allclose(
+                        abs(mathematica_freq_output[i]), abs(python_freq_output[i])
+                    )
                 )

tests/test_initial_conditions.py

@@ -7,11 +7,13 @@
 
 DATA_DIR = THIS_DIR.parent / "tests/data"
 
-stable_orbit_values = np.genfromtxt(DATA_DIR / "stable_orbit_parameters.txt", delimiter=",")
-complex_plunging_orbit_values = np.genfromtxt(
+stable_orbit_parameters = np.genfromtxt(
+    DATA_DIR / "stable_orbit_parameters.txt", delimiter=","
+)
+complex_plunging_orbit_parameters = np.genfromtxt(
     DATA_DIR / "plunging_orbit_parameters_complex.txt", delimiter=","
 )
-real_plunging_orbit_values = np.genfromtxt(
+real_plunging_orbit_parameters = np.genfromtxt(
     DATA_DIR / "plunging_orbit_parameters_real.txt", delimiter=","
 )