Merge pull request #3 from BlackHolePerturbationToolkit/deltas

Fix #2 and expose deltas

README.md

@@ -1,6 +1,6 @@
 # KerrGeoPy
 
-KerrGeoPy is a python implementation of the [KerrGeodesics](https://bhptoolkit.org/KerrGeodesics/) Mathematica library. It is intended for use in computing orbital trajectories for extreme-mass-ratio inspirals (EMRIs). It implements the analytical solutions for plunging orbits from [Dyson and van de Meent](https://arxiv.org/abs/2302.03704), as well as solutions for stable orbits from [Fujita and Hikida](https://arxiv.org/abs/0906.1420). The library also provides a set of methods for computing constants of motion and orbital frequencies.
+KerrGeoPy is a python implementation of the [KerrGeodesics](https://bhptoolkit.org/KerrGeodesics/) Mathematica library. It is intended for use in computing orbital trajectories for extreme-mass-ratio inspirals (EMRIs). It implements the analytical solutions for plunging orbits from [Dyson and van de Meent](https://arxiv.org/abs/2302.03704), as well as solutions for stable orbits from [Fujita and Hikida](https://arxiv.org/abs/0906.1420). The library also provides a set of methods for computing constants of motion and orbital frequencies. See the [documentation](https://kerrgeopy.readthedocs.io/en/latest/) for more information.
 
 ## Installation
 
@@ -279,4 +279,4 @@ plt.ylabel(r"$\phi(\lambda)$")
 ## Authors
 
 * Seyong Park
-* Zach Nasipak
+* Zach Nasipak

docs/source/_autosummary/kerrgeopy.orbit.Orbit.rst

@@ -24,6 +24,7 @@
       ~Orbit.numerical_four_velocity
       ~Orbit.plot
       ~Orbit.trajectory
+      ~Orbit.trajectory_deltas
 
 
 

docs/source/_autosummary/kerrgeopy.plunge.PlungingOrbit.rst

@@ -24,6 +24,7 @@
       ~PlungingOrbit.numerical_four_velocity
       ~PlungingOrbit.plot
       ~PlungingOrbit.trajectory
+      ~PlungingOrbit.trajectory_deltas
 
 
 

docs/source/_autosummary/kerrgeopy.stable.StableOrbit.rst

@@ -27,6 +27,7 @@
       ~StableOrbit.numerical_four_velocity
       ~StableOrbit.plot
       ~StableOrbit.trajectory
+      ~StableOrbit.trajectory_deltas
 
 
 

docs/source/index.rst

@@ -19,7 +19,7 @@ The library also provides a set of methods for computing constants of motion and
 .. raw:: html
 
    <video width="45%" autoplay loop>
-      <source src="https://raw.githubusercontent.com/BlackHolePerturbationToolkit/KerrGeoPy/main/docs/source/notebooks/animation4.mp4" type="video/mp4">
+      <source src="https://cdn.jsdelivr.net/gh/BlackHolePerturbationToolkit/KerrGeoPy@main/docs/source/notebooks/animation4.mp4" type="video/mp4">
    </video>
 
 .. _Installation:

docs/source/notebooks/Graphics.ipynb

@@ -26,7 +26,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 1,
    "id": "7ae99533",
    "metadata": {},
    "outputs": [
@@ -171,7 +171,7 @@
    "metadata": {},
    "source": [
     "<video width=\"75%\" controls autoplay loop>\n",
-    "  <source src=\"https://raw.githubusercontent.com/BlackHolePerturbationToolkit/KerrGeoPy/main/docs/source/notebooks/animation1.mp4\" type=\"video/mp4\">\n",
+    "  <source src=\"https://cdn.jsdelivr.net/gh/BlackHolePerturbationToolkit/KerrGeoPy@main/docs/source/notebooks/animation1.mp4\" type=\"video/mp4\">\n",
     "  <p>\n",
     "\n",
     "  [](animation1.mp4)\n",
@@ -203,7 +203,7 @@
    "metadata": {},
    "source": [
     "<video width=\"75%\" controls autoplay loop>\n",
-    "  <source src=\"https://raw.githubusercontent.com/BlackHolePerturbationToolkit/KerrGeoPy/main/docs/source/notebooks/animation2.mp4\" type=\"video/mp4\">\n",
+    "  <source src=\"https://cdn.jsdelivr.net/gh/BlackHolePerturbationToolkit/KerrGeoPy@main/docs/source/notebooks/animation2.mp4\" type=\"video/mp4\">\n",
     "  <p>\n",
     "\n",
     "  [](animation2.mp4)\n",
@@ -234,7 +234,7 @@
    "metadata": {},
    "source": [
     "<video width=\"75%\" controls autoplay loop>\n",
-    "  <source src=\"https://raw.githubusercontent.com/BlackHolePerturbationToolkit/KerrGeoPy/main/docs/source/notebooks/animation3.mp4\" type=\"video/mp4\">\n",
+    "  <source src=\"https://cdn.jsdelivr.net/gh/BlackHolePerturbationToolkit/KerrGeoPy@main/docs/source/notebooks/animation3.mp4\" type=\"video/mp4\">\n",
     "  <p>\n",
     "\n",
     "  [](animation3.mp4)\n",
@@ -267,7 +267,7 @@
    "metadata": {},
    "source": [
     "<video width=\"75%\" controls autoplay loop>\n",
-    "  <source src=\"https://raw.githubusercontent.com/BlackHolePerturbationToolkit/KerrGeoPy/main/docs/source/notebooks/animation4.mp4\" type=\"video/mp4\">\n",
+    "  <source src=\"https://cdn.jsdelivr.net/gh/BlackHolePerturbationToolkit/KerrGeoPy@main/docs/source/notebooks/animation4.mp4\" type=\"video/mp4\">\n",
     "  <p>\n",
     "\n",
     "  [](animation4.mp4)\n",
@@ -305,7 +305,7 @@
    "metadata": {},
    "source": [
     "<video width=\"100%\" controls autoplay loop>\n",
-    "  <source src=\"https://raw.githubusercontent.com/BlackHolePerturbationToolkit/KerrGeoPy/main/docs/source/notebooks/animation5.mp4\" type=\"video/mp4\">\n",
+    "  <source src=\"https://cdn.jsdelivr.net/gh/BlackHolePerturbationToolkit/KerrGeoPy@main/docs/source/notebooks/animation5.mp4\" type=\"video/mp4\">\n",
     "  <p>\n",
     "\n",
     "  [](animation5.mp4)\n",

kerrgeopy/frequencies.py

@@ -4,9 +4,26 @@
 """
 from .constants import _standardize_params
 from .constants import *
-from scipy.special import ellipk, ellipe, elliprj, elliprf
+from scipy.special import ellipk, ellipe, elliprj, elliprf, elliprd
 from numpy import sin, cos, sqrt, pi, arcsin, floor, where
 
+def _ellipeinc(phi,m):
+    r"""
+    Incomplete elliptic integral of the second kind defined as :math:`E(\phi,m) = \int_0^{\phi} \sqrt{1-m\sin^2\theta}d\theta`.
+    """
+    # 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*m*sin(phi)**3*elliprd(cos(phi)**2,1-m*sin(phi)**2,1)
+    result = where(num_cycles % 2 == 0, num_cycles*ellipe(m)+integral, (num_cycles+1)*ellipe(m)-integral)
+
+    # return scalar for scalar input
+    return result.item() if np.isscalar(phi) else result
+
 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}}}`
@@ -40,7 +57,10 @@ def _ellippiinc(phi,n,m):
     # 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)
+    result =  where(num_cycles % 2 == 0, num_cycles*_ellippi(n,m)+integral, (num_cycles+1)*_ellippi(n,m)-integral)
+
+    # return scalar for scalar input
+    return result.item() if np.isscalar(phi) else result
 
 def r_frequency(a,p,e,x,constants=None):
     """

kerrgeopy/orbit.py

@@ -4,7 +4,7 @@
 from .spacetime import KerrSpacetime
 from .initial_conditions import *
 from .units import *
-from .constants import scale_constants, apex_from_constants
+from .constants import scale_constants, apex_from_constants, stable_polar_roots, stable_radial_roots
 from .frequencies import mino_frequencies, fundamental_frequencies
 from numpy import sin, cos, sqrt, pi
 import numpy as np
@@ -64,6 +64,46 @@ def __init__(self,a,initial_position,initial_velocity, M=None, mu=None):
             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 trajectory_deltas(self,initial_phases=None):
+        r"""
+        Computes the trajectory deltas :math:`t_r(q_r)`, :math:`t_\theta(q_\theta)`, :math:`\phi_r(q_r)` and :math:`\phi_\theta(q_\theta)`
+
+        :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
+
+        :return: tuple of trajectory deltas :math:`(t_r(q_r), t_\theta(q_\theta), \phi_r(q_r),\phi_\theta(q_\theta))`
+        :rtype: tuple(function, function, function, function)
+        """
+        if initial_phases is None: initial_phases = self.initial_phases
+        q_t0, q_r0, q_theta0, q_phi0 = initial_phases
+
+        if self.stable:
+            from .stable import radial_solutions, polar_solutions
+
+            constants = (self.E,self.L,self.Q)
+            radial_roots = stable_radial_roots(self.a,self.p,self.e,self.x,constants)
+            polar_roots = stable_polar_roots(self.a,self.p,self.e,self.x,constants)
+            r, t_r, phi_r = radial_solutions(self.a,constants,radial_roots)
+            theta, t_theta, phi_theta = polar_solutions(self.a,constants,polar_roots)
+        else:
+            radial_roots = plunging_radial_roots(self.a,self.E,self.L,self.Q)
+            if np.iscomplex(radial_roots[3]):
+                from .plunge import plunging_radial_solutions_complex, plunging_polar_solutions
+
+                # adjust q_theta0 so that initial conditions are consistent with stable orbits
+                q_theta0 = q_theta0 + pi/2
+                r, t_r, phi_r = plunging_radial_solutions_complex(self.a,self.E,self.L,self.Q)
+                theta, t_theta, phi_theta = plunging_polar_solutions(self.a,self.E,self.L,self.Q)
+            else:
+                from .stable import radial_solutions, polar_solutions
+
+                constants = (self.E,self.L,self.Q)
+                r, t_r, phi_r = radial_solutions(self.a,constants,radial_roots)
+                theta, t_theta, phi_theta = polar_solutions(self.a,constants,radial_roots)
+
+        return (lambda q_r: t_r(q_r+q_r0), lambda q_theta: t_theta(q_theta+q_theta0),
+                lambda q_r: phi_r(q_r+q_r0), lambda q_theta: phi_theta(q_theta+q_theta0))
+
     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
@@ -192,16 +232,16 @@ def numerical_four_velocity(self,dx=1e-6,initial_phases=None):
 
         def u_t(mino_time):
             sigma = r(mino_time)**2 + self.a**2*cos(theta(mino_time))**2
-            return (t(mino_time+dx)-t(mino_time-dx))/(2*dx*sigma)
+            return (-t(mino_time+2*dx)+8*t(mino_time+dx)-8*t(mino_time-dx)+t(mino_time-2*dx))/(12*dx*sigma)
         def u_r(mino_time):
             sigma = r(mino_time)**2 + self.a**2*cos(theta(mino_time))**2
-            return (r(mino_time+dx)-r(mino_time-dx))/(2*dx*sigma)
+            return (-r(mino_time+2*dx)+8*r(mino_time+dx)-8*r(mino_time-dx)+r(mino_time-2*dx))/(12*dx*sigma)
         def u_theta(mino_time):
             sigma = r(mino_time)**2 + self.a**2*cos(theta(mino_time))**2
-            return (theta(mino_time+dx)-theta(mino_time-dx))/(2*dx*sigma)
+            return (-theta(mino_time+2*dx)+8*theta(mino_time+dx)-8*theta(mino_time-dx)+theta(mino_time-2*dx))/(12*dx*sigma)
         def u_phi(mino_time):
             sigma = r(mino_time)**2 + self.a**2*cos(theta(mino_time))**2
-            return (phi(mino_time+dx)-phi(mino_time-dx))/(2*dx*sigma)
+            return (-phi(mino_time+2*dx)+8*phi(mino_time+dx)-8*phi(mino_time-dx)+phi(mino_time-2*dx))/(12*dx*sigma)
         return u_t, u_r, u_theta, u_phi
 
     def plot(self,lambda0=0, lambda1=10, elevation=30 ,azimuth=-60, initial_phases=None, grid=True, axes=True, lw=1,color="red",tau=np.inf,point_density=200):
@@ -397,6 +437,7 @@ def animate(self,filename,lambda0=0, lambda1=10, elevation=30 ,azimuth=-60, init
             per_subplot_kw={"O":{"projection":"3d"},"T":{"facecolor":"none"},"R":{"facecolor":"none"},"Θ":{"facecolor":"none"},"Φ":{"facecolor":"none"}}
             )
             ax = ax_dict["O"]
+
             ax_dict["T"].set_ylabel("$t$")
             ax_dict["R"].set_ylabel("$r$")
             ax_dict["Θ"].set_ylabel(r"$\theta$")

kerrgeopy/plunge.py

@@ -54,6 +54,33 @@ def __init__(self,a,E,L,Q,initial_phases=(0,0,0,0),M=None,mu=None):
         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_deltas(self,initial_phases=None):
+        r"""
+        Computes the trajectory deltas :math:`t_r(q_r)`, :math:`t_\theta(q_\theta)`, :math:`\phi_r(q_r)` and :math:`\phi_\theta(q_\theta)`
+
+        :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
+
+        :return: tuple of trajectory deltas :math:`(t_r(q_r), t_\theta(q_\theta), \phi_r(q_r),\phi_\theta(q_\theta))`
+        :rtype: tuple(function, function, function, function)
+        """
+        if initial_phases is None: initial_phases = self.initial_phases
+        q_t0, q_r0, q_theta0, q_phi0 = initial_phases
+
+        radial_roots = plunging_radial_roots(self.a,self.E,self.L,self.Q)
+        if np.iscomplex(radial_roots[3]):
+            # adjust q_theta0 so that initial conditions are consistent with stable orbits
+            q_theta0 = q_theta0 + pi/2
+            r, t_r, phi_r = plunging_radial_solutions_complex(self.a,self.E,self.L,self.Q)
+            theta, t_theta, phi_theta = plunging_polar_solutions(self.a,self.E,self.L,self.Q)
+        else:
+            constants = (self.E,self.L,self.Q)
+            r, t_r, phi_r = radial_solutions(self.a,constants,radial_roots)
+            theta, t_theta, phi_theta = polar_solutions(self.a,constants,radial_roots)
+
+        return (lambda q_r: t_r(q_r+q_r0), lambda q_theta: t_theta(q_theta+q_theta0),
+                lambda q_r: phi_r(q_r+q_r0), lambda q_theta: phi_theta(q_theta+q_theta0))
 
     def trajectory(self,initial_phases=None,distance_units="natural",time_units="natural"):
         r"""

kerrgeopy/stable.py

@@ -3,9 +3,9 @@
 """
 from .constants import *
 from .constants import _standardize_params
-from .frequencies import _ellippi, _ellippiinc
+from .frequencies import _ellippi, _ellippiinc, _ellipeinc
 from .frequencies import *
-from scipy.special import ellipj, ellipeinc
+from scipy.special import ellipj
 from .orbit import Orbit
 from numpy import pi, arccos
 
@@ -37,6 +37,7 @@ class StableOrbit(Orbit):
     :ivar E: dimensionless energy
     :ivar L: dimensionless angular momentum
     :ivar Q: dimensionless carter constant
+    :ivar initial_phases: tuple of initial phases :math:`(q_{t_0},q_{r_0},q_{\theta_0},q_{\phi_0})`
     :ivar stable: boolean indicating whether the orbit is stable
     :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)`
@@ -129,6 +130,28 @@ def fundamental_frequencies(self, units="natural"):
             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_deltas(self,initial_phases=None):
+        r"""
+        Computes the trajectory deltas :math:`t_r(q_r)`, :math:`t_\theta(q_\theta)`, :math:`\phi_r(q_r)` and :math:`\phi_\theta(q_\theta)`
+
+        :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
+
+        :return: tuple of trajectory deltas :math:`(t_r(q_r), t_\theta(q_\theta), \phi_r(q_r),\phi_\theta(q_\theta))`
+        :rtype: tuple(function, function, function, function)
+        """
+        if initial_phases is None: initial_phases = self.initial_phases
+        q_t0, q_r0, q_theta0, q_phi0 = initial_phases
+
+        constants = (self.E,self.L,self.Q)
+        radial_roots = stable_radial_roots(self.a,self.p,self.e,self.x,constants)
+        polar_roots = stable_polar_roots(self.a,self.p,self.e,self.x,constants)
+        r, t_r, phi_r = radial_solutions(self.a,constants,radial_roots)
+        theta, t_theta, phi_theta = polar_solutions(self.a,constants,polar_roots)
+
+        return (lambda q_r: t_r(q_r+q_r0), lambda q_theta: t_theta(q_theta+q_theta0),
+                lambda q_r: phi_r(q_r+q_r0), lambda q_theta: phi_theta(q_theta+q_theta0))
 
     def trajectory(self,initial_phases=None,distance_units="natural",time_units="natural"):
         r"""
@@ -189,12 +212,14 @@ 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)
+        # adding 1e-14 to avoid strange discontinuity in ellipeinc when q_r is set to specific ratios of pi
+        #psi_r = psi_r+1e-14
         # 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))
+            + (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) * \
@@ -254,7 +279,7 @@ 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))
+        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)

setup.py

@@ -1,12 +1,18 @@
 from setuptools import setup
 import setuptools
+# read the contents of your README file
+from pathlib import Path
+this_directory = Path(__file__).parent
+long_description = (this_directory / "README.md").read_text()
 
 setup(
     name="kerrgeopy",
     version="0.9.1",
     author="Seyong Park",
     description="Library for computing stable and plunging geodesics in Kerr spacetime",
-    url="https://github.com/syp2001/KerrGeoPy",
+    long_description=long_description,
+    long_description_content_type='text/markdown',
+    url="https://github.com/BlackHolePerturbationToolkit/KerrGeoPy",
     packages=setuptools.find_packages(exclude=["tests"]),
     classifiers=(
         "Programming Language :: Python :: 3",