fixed initial conditions

kerrgeopy/constants.py

@@ -2,8 +2,7 @@
 Module containing functions for computing the constants of motion for bound orbits using the :math:`(a,p,e,x)` parametrization.
 Constants of motion are computed using the method described in Appendix B of `Schmidt <https://doi.org/10.48550/arXiv.gr-qc/0202090>`_.
 """
-from numpy import sign, sqrt, copysign, inf, nan
-from math import pi
+from numpy import sign, sqrt, copysign, pi, nan, inf
 from scipy.optimize import root_scalar
 from .units import *
 from scipy.interpolate import RectBivariateSpline

kerrgeopy/initial_conditions.py

@@ -7,7 +7,7 @@
 import numpy as np
 from numpy import sin, cos, sqrt, pi, arcsin, arccos
 from numpy.polynomial import Polynomial
-from scipy.special import ellipkinc, ellipk
+from scipy.special import ellipkinc
 
 def apex_from_constants(a,E,L,Q):
     r"""

kerrgeopy/orbit.py

@@ -6,10 +6,8 @@
 from .units import *
 from .constants import scale_constants
 from .frequencies import mino_frequencies, fundamental_frequencies
-from numpy import sin, cos, sqrt, pi, arcsin
-from numpy.polynomial import Polynomial
+from numpy import sin, cos, sqrt, pi
 import numpy as np
-from scipy.special import ellipkinc
 import matplotlib.pyplot as plt
 from matplotlib.animation import FuncAnimation
 from matplotlib.animation import FFMpegWriter
@@ -35,23 +33,22 @@ class Orbit:
     :ivar E: dimensionless energy
     :ivar L: dimensionless angular momentum
     :ivar Q: dimensionless carter constant
+    :ivar stable: boolean indicating whether the orbit is stable
+    :ivar upsilon_r: dimensionless radial orbital frequency in Mino time
+    :ivar upsilon_theta: dimensionless polar orbital frequency in Mino time
     """
     def __init__(self,a,initial_position,initial_velocity, M=None, mu=None):
         self.a, self.initial_position, self.initial_velocity, self.M, self.mu = a, initial_position, initial_velocity, M, mu
 
-        t0, r0, theta0, phi0 = initial_position
-        dt0, dr0, dtheta0, dphi0 = initial_velocity
-
         # check if initial four-velocity is valid
         spacetime = KerrSpacetime(a)
-        initial_norm =  spacetime.norm(t0,r0,theta0,phi0,initial_velocity)
+        initial_norm =  spacetime.norm(*initial_position,initial_velocity)
         if initial_norm >= 0: raise ValueError("Initial velocity is not timelike")
         if abs(initial_norm + 1) > 1e-6: raise ValueError("Initial velocity is not normalized")
 
         E, L, Q = constants_from_initial_conditions(a,initial_position,initial_velocity)
         self.E, self.L, self.Q = E, L, Q
 
-
         if is_stable(a,initial_position,initial_velocity):
             self.stable = True
             a, p, e, x = apex_from_constants(a,E,L,Q)
@@ -63,8 +60,19 @@ def __init__(self,a,initial_position,initial_velocity, M=None, mu=None):
             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.
 
-    def trajectory(self,distance_units="natural",time_units="natural"):
+        :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"""
         Computes the components of the trajectory as a function of Mino time
 
@@ -76,14 +84,15 @@ def trajectory(self,distance_units="natural",time_units="natural"):
         :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 self.stable:
             from .stable_orbit import StableOrbit
             orbit = StableOrbit(self.a,self.p,self.e,self.x,self.M,self.mu)
-            return orbit.trajectory(self.initial_phases,distance_units,time_units)
+            return orbit.trajectory(initial_phases,distance_units,time_units)
         else:
             from .plunging_orbit import PlungingOrbit
             orbit = PlungingOrbit(self.a,self.E,self.L,self.Q,self.M,self.mu)
-            return orbit.trajectory(self.initial_phases,distance_units,time_units)
+            return orbit.trajectory(initial_phases,distance_units,time_units)
 
     def constants_of_motion(self, units="natural"):
         """
@@ -112,7 +121,7 @@ def constants_of_motion(self, units="natural"):
 
         raise ValueError("units must be one of 'natural', 'mks', or 'cgs'")
 
-    def four_velocity(self,initial_phases=(0,0,0,0)):
+    def four_velocity(self,initial_phases=None):
         r"""
         Computes the four velocity of the orbit as a function of Mino time
 
@@ -122,13 +131,14 @@ def four_velocity(self,initial_phases=(0,0,0,0)):
         :return: components of the four velocity (i.e. :math:`u^t,u^r,u^\theta,u^\phi`)
         :rtype: tuple(function, function, function, function)
         """
+        if initial_phases is None: initial_phases = self.initial_phases
         t, r, theta, phi = self.trajectory(initial_phases)
         spacetime = KerrSpacetime(self.a)
         constants = self.E, self.L, self.Q
 
         return spacetime.four_velocity(t,r,theta,phi,constants,self.upsilon_r,self.upsilon_theta,initial_phases)
 
-    def four_velocity_norm(self,initial_phases=(0,0,0,0)):
+    def four_velocity_norm(self,initial_phases=None):
         r"""
         Computes the norm of the four velocity of the orbit as a function of Mino time
 
@@ -138,6 +148,7 @@ def four_velocity_norm(self,initial_phases=(0,0,0,0)):
         :return: norm of the four velocity :math:`g_{\mu\nu}u^\mu u^\nu`
         :rtype: function
         """
+        if initial_phases is None: initial_phases = self.initial_phases
         t, r, theta, phi = self.trajectory(initial_phases)
         spacetime = KerrSpacetime(self.a)
         constants = self.E, self.L, self.Q
@@ -178,7 +189,7 @@ def plot(self,lambda0=0, lambda1=20, elevation=30 ,azimuth=-60, initial_phases=(
         num_pts = int(lambda_range*point_density)
         time = np.linspace(lambda0,lambda1,num_pts)
 
-        t, r, theta, phi = self.trajectory(initial_phases)
+        t, r, theta, phi = self.trajectory(initial_phases=initial_phases)
 
         # compute trajectory in cartesian coordinates
         trajectory_x = r(time)*sin(theta(time))*cos(phi(time))

kerrgeopy/plunging_orbit.py

@@ -5,10 +5,9 @@
 from .plunging_solutions import _plunging_radial_roots
 from .stable_solutions import *
 from .orbit import Orbit
-import matplotlib.pyplot as plt
 
 class PlungingOrbit(Orbit):
-    """
+    r"""
     Class representing a plunging orbit in Kerr spacetime.
 
     :param a: dimensionaless spin parameter
@@ -28,16 +27,24 @@ class PlungingOrbit(Orbit):
     :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,M=None,mu=None):
-        self.a, self.E, self.L, self.Q, self.M, self.mu = a, E, L, Q, M, mu
+    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.upsilon_r, self.upsilon_theta = plunging_mino_frequencies(a,E,L,Q)
+        self.stable = False
+
+        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=(0,0,0,0),distance_units="natural",time_units="natural"):
+    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
 
@@ -63,7 +70,7 @@ def trajectory(self,initial_phases=(0,0,0,0),distance_units="natural",time_units
         else:
             return self._real_trajectory(initial_phases,distance_units,time_units)
 
-    def _complex_trajectory(self,initial_phases=(0,0,0,0), distance_units="natural",time_units="natural"):
+    def _complex_trajectory(self,initial_phases=None, distance_units="natural",time_units="natural"):
         """
         Computes the components of the trajectory in the case of two complex and two real radial roots
 
@@ -72,6 +79,7 @@ def _complex_trajectory(self,initial_phases=(0,0,0,0), distance_units="natural",
 
         :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)
@@ -90,7 +98,7 @@ def _complex_trajectory(self,initial_phases=(0,0,0,0), distance_units="natural",
 
 
         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) + a*L*mino_time - C_t + q_t0, self.M)
+            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)
@@ -99,11 +107,11 @@ 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) - a*E*mino_time - C_phi + q_phi0
+            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=(0,0,0,0),distance_units="natural",time_units="natural"):
+    def _real_trajectory(self,initial_phases=None,distance_units="natural",time_units="natural"):
         """
         Computes the components of the trajectory in the case of four real radial roots
 
@@ -112,6 +120,7 @@ def _real_trajectory(self,initial_phases=(0,0,0,0),distance_units="natural",time
 
         :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)

kerrgeopy/plunging_solutions.py

@@ -1,7 +1,7 @@
 """
 Module containing functions to compute intermediate terms in the plunging orbit solutions of `Dyson and van de Meent <https://doi.org/10.48550/arXiv.2302.03704>`_
 """
-from numpy import sort_complex, sqrt, arctan, arctan2, arccos, log, sin, cos, pi
+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
@@ -132,7 +132,7 @@ def I_r(q_r):
         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))) \
+        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):
@@ -224,15 +224,15 @@ def t_r(q_r):
             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)
-            ) - a*L*mino_time
+            )
 
     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)
-            ) + a*E*mino_time
+            )
 
     return r, t_r, phi_r
 
@@ -277,6 +277,6 @@ def phi_theta(q_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 L/z2*_ellippiinc(xi_theta,z1**2,k_theta)
 
     return theta, t_theta, phi_theta

kerrgeopy/spacetime.py

@@ -3,7 +3,6 @@
 """
 from .constants import *
 from numpy import cos, sin
-from numpy.polynomial import Polynomial
 import numpy as np
 
 class KerrSpacetime:

kerrgeopy/stable_orbit.py

@@ -9,7 +9,7 @@
 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)
@@ -34,6 +34,8 @@ class StableOrbit(Orbit):
     :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
@@ -42,19 +44,25 @@ class StableOrbit(Orbit):
     :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,M = None,mu=None):
+    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.M, self.mu = a, p, e, x, M, mu
+        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,M=None,mu=None):
+    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.
 
@@ -73,7 +81,7 @@ def from_constants(cls,a,E,L,Q,M=None,mu=None):
         """
         a, p, e, x = apex_from_constants(a,E,L,Q)
 
-        return cls(a,p,e,x,M,mu)
+        return cls(a,p,e,x,initial_phases,M,mu)
 
     def mino_frequencies(self, units="natural"):
         r"""
@@ -120,7 +128,7 @@ def fundamental_frequencies(self, units="natural"):
 
         raise ValueError("units must be one of 'natural', 'mks', 'cgs', or 'mHz'")
 
-    def trajectory(self,initial_phases=(0,0,0,0),distance_units="natural",time_units="natural"):
+    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.
 
@@ -134,7 +142,7 @@ def trajectory(self,initial_phases=(0,0,0,0),distance_units="natural",time_units
         :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")
 

kerrgeopy/stable_solutions.py

@@ -5,7 +5,7 @@
 from .frequencies_from_constants import _ellippi
 from .frequencies_from_constants import *
 from scipy.special import ellipj, ellipeinc
-from numpy import sin, cos, arcsin, arccos, floor, where
+from numpy import pi, sin, cos, arcsin, arccos, floor, where
 
 def _ellippiinc(phi,n,k):
     r"""