diff --git a/kerrgeopy/constants.py b/kerrgeopy/constants.py index 8e0c146..1260eac 100644 --- a/kerrgeopy/constants.py +++ b/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 diff --git a/kerrgeopy/initial_conditions.py b/kerrgeopy/initial_conditions.py index 5923291..8e81e67 100644 --- a/kerrgeopy/initial_conditions.py +++ b/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""" diff --git a/kerrgeopy/orbit.py b/kerrgeopy/orbit.py index b8c4c21..af9e4fc 100644 --- a/kerrgeopy/orbit.py +++ b/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)) diff --git a/kerrgeopy/plunging_orbit.py b/kerrgeopy/plunging_orbit.py index c8ef4d1..50bd8de 100644 --- a/kerrgeopy/plunging_orbit.py +++ b/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) diff --git a/kerrgeopy/plunging_solutions.py b/kerrgeopy/plunging_solutions.py index aaa9a6e..a3a242b 100644 --- a/kerrgeopy/plunging_solutions.py +++ b/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,7 +224,7 @@ 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 @@ -232,7 +232,7 @@ def phi_r(q_r): 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 diff --git a/kerrgeopy/spacetime.py b/kerrgeopy/spacetime.py index d931030..d963406 100644 --- a/kerrgeopy/spacetime.py +++ b/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: diff --git a/kerrgeopy/stable_orbit.py b/kerrgeopy/stable_orbit.py index 8d0b6d7..6e1b2d8 100644 --- a/kerrgeopy/stable_orbit.py +++ b/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") diff --git a/kerrgeopy/stable_solutions.py b/kerrgeopy/stable_solutions.py index d674977..83feb22 100644 --- a/kerrgeopy/stable_solutions.py +++ b/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"""