Skip to content

Celestial mechanics

Celestial mechanics: orbits under gravity and relativity, spin-orbit coupling, and charges in magnetic fields.

DipoleChargeElement(label, visual=None, *, x=1.0, x_d1=0.0, y=0.0, y_d1=0.0, z=0.0, z_d1=0.0, alpha=1.0)

Bases: Element

Charge in the magnetic dipole field of a planet, at the position (x, y, z) with coordinates of order 2.

The parameter alpha is the product of the charge and the dipole moment, divided by the mass.

Parameters:

Name Type Description Default
label str

The label of the charge.

 required
visual VisualSpec | None

How the charge is drawn, or None for its class visual.

 None
x Annotated[float, StateSpec(order=2)]

The initial first coordinate.

 1.0
x_d1 float

The initial first velocity.

 0.0
y Annotated[float, StateSpec(order=2)]

The initial second coordinate.

 0.0
y_d1 float

The initial second velocity.

 0.0
z Annotated[float, StateSpec(order=2)]

The initial third coordinate, along the dipole.

 0.0
z_d1 float

The initial third velocity.

 0.0
alpha Annotated[float, ParameterSpec()]

The product of the charge and the dipole moment, divided by the mass.

 1.0

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import DipoleChargeElement
>>> DipoleChargeElement('charge', alpha=2.0).alpha
2.0
Source code in src/skcomplex/systems/physics/celestial_mechanics/_dipoles.py 
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(order=2)] = 1.0,
    x_d1: float = 0.0,
    y: Annotated[float, StateSpec(order=2)] = 0.0,
    y_d1: float = 0.0,
    z: Annotated[float, StateSpec(order=2)] = 0.0,
    z_d1: float = 0.0,
    alpha: Annotated[float, ParameterSpec()] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.x_d1 = x_d1
    self.y = y
    self.y_d1 = y_d1
    self.z = z
    self.z_d1 = z_d1
    self.alpha = alpha

DumbbellElement(label, visual=None, *, u=1.0, u_d1=0.0, theta=0.0, theta_d1=1.0, v=0.0, v_d1=0.0, phi=0.0, phi_d1=1.0, psi=math.pi / 2, psi_d1=0.0, delta=0.5, epsilon=0.1, mu=0.9)

Bases: Element

Rigid dumbbell orbiting a primary: the position of its centre of mass and the orientation of its arm.

The cylindrical coordinates (u, theta, v) place the centre of mass of the dumbbell relative to the primary, and the spherical angles (phi, psi) give the direction of its arm, of length epsilon in units of the mean orbit radius. The mass 1 - delta of the dumbbell sits at delta epsilon along the arm and the mass delta at 1 - delta against it. All five coordinates are states of order 2 that a Lagrangian moves.

Parameters:

Name Type Description Default
label str

The label of the dumbbell.

 required
visual VisualSpec | None

How the dumbbell is drawn, or None for default_visual_.

 None
u Annotated[float, StateSpec(order=2)]

The initial distance from the axis through the primary.

 1.0
u_d1 float

The initial radial velocity.

 0.0
theta Annotated[float, StateSpec(order=2)]

The initial orbital angle.

 0.0
theta_d1 float

The initial orbital angular velocity.

 1.0
v Annotated[float, StateSpec(order=2)]

The initial height above the orbital plane.

 0.0
v_d1 float

The initial vertical velocity.

 0.0
phi Annotated[float, StateSpec(order=2)]

The initial azimuth of the arm.

 0.0
phi_d1 float

The initial spin of the arm about the vertical.

 1.0
psi Annotated[float, StateSpec(order=2)]

The initial polar angle of the arm, pi / 2 in the orbital plane.

 pi / 2
psi_d1 float

The initial rate of the polar angle of the arm.

 0.0
delta Annotated[float, ParameterSpec(domain=open(0, 1))]

The share of the second mass in the mass of the dumbbell, 1 / 2 for equal masses.

 0.5
epsilon Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The length of the arm in units of the mean orbit radius.

 0.1
mu Annotated[float, ParameterSpec(domain=open(0, 1))]

The share of the primary in the mass of the three bodies.

 0.9

Attributes:

Name Type Description
default_visual_

The visual of the dumbbell when it sets none.

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import DumbbellElement
>>> DumbbellElement('dumbbell').delta
0.5
Source code in src/skcomplex/systems/physics/celestial_mechanics/_spin_orbit.py 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    u: Annotated[float, StateSpec(order=2)] = 1.0,
    u_d1: float = 0.0,
    theta: Annotated[float, StateSpec(order=2)] = 0.0,
    theta_d1: float = 1.0,
    v: Annotated[float, StateSpec(order=2)] = 0.0,
    v_d1: float = 0.0,
    phi: Annotated[float, StateSpec(order=2)] = 0.0,
    phi_d1: float = 1.0,
    psi: Annotated[float, StateSpec(order=2)] = math.pi / 2,
    psi_d1: float = 0.0,
    delta: Annotated[float, ParameterSpec(domain=RealSet.open(0, 1))] = 0.5,
    epsilon: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    mu: Annotated[float, ParameterSpec(domain=RealSet.open(0, 1))] = 0.9,
) -> None:
    super().__init__(label, visual=visual)
    self.u = u
    self.u_d1 = u_d1
    self.theta = theta
    self.theta_d1 = theta_d1
    self.v = v
    self.v_d1 = v_d1
    self.phi = phi
    self.phi_d1 = phi_d1
    self.psi = psi
    self.psi_d1 = psi_d1
    self.delta = delta
    self.epsilon = epsilon
    self.mu = mu

GeodesicParticleElement(label, x=0.0, x_d1=0.0, visual=None)

Bases: Element

Test particle at the spacetime position x = (t, r, theta, phi) around a central mass, with its four-velocity.

The position is a vector state of order 2 in the Schwarzschild or Boyer-Lindquist coordinates of the space, and its companion x_d1 is the four-velocity, the derivative of the position by the proper time of a massive particle or by the affine parameter of a light ray. The particle is drawn in the equatorial plane, at (r cos phi, r sin phi).

Parameters:

Name Type Description Default
label str

The label of the particle.

 required
x Annotated[ArrayLike, StateSpec(vector=True, order=2)]

The initial position (t, r, theta, phi).

 0.0
x_d1 ArrayLike

The initial four-velocity.

 0.0
visual VisualSpec | None

How the particle is drawn, or None for its class visual.

 None

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import GeodesicParticleElement
>>> GeodesicParticleElement('particle', x=[0.0, 10.0, 1.5708, 0.0]).x
[0.0, 10.0, 1.5708, 0.0]
Source code in src/skcomplex/systems/physics/celestial_mechanics/_geodesics.py 
65
66
67
68
69
70
71
72
73
74
def __init__(
    self,
    label: str,
    x: Annotated[npt.ArrayLike, StateSpec(vector=True, order=2)] = 0.0,
    x_d1: npt.ArrayLike = 0.0,
    visual: VisualSpec | None = None,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.x_d1 = x_d1

KeplerBodyElement(label, visual=None, *, x=10.0, y=0.0, px=0.0, py=10.0, G=1.0, mass=1000.0)

Bases: Element

Planet at (x, y) relative to a heavy star, with the canonical momenta (px, py) of the Kepler problem.

The states are of order 1, and the star of mass mass sits at the origin.

Parameters:

Name Type Description Default
label str

The label of the planet.

 required
visual VisualSpec | None

How the planet is drawn, or None for its class visual.

 None
x Annotated[float, StateSpec()]

The initial first coordinate.

 10.0
y Annotated[float, StateSpec()]

The initial second coordinate.

 0.0
px Annotated[float, StateSpec()]

The initial first momentum per unit mass.

 0.0
py Annotated[float, StateSpec()]

The initial second momentum per unit mass.

 10.0
G Annotated[float, ParameterSpec(domain=open(0, oo))]

The gravitational constant.

 1.0
mass Annotated[float, ParameterSpec(domain=open(0, oo))]

The mass of the star.

 1000.0

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import KeplerBodyElement
>>> KeplerBodyElement('planet', mass=10.0).mass
10.0
Source code in src/skcomplex/systems/physics/celestial_mechanics/_orbits.py 
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 10.0,
    y: Annotated[float, StateSpec()] = 0.0,
    px: Annotated[float, StateSpec()] = 0.0,
    py: Annotated[float, StateSpec()] = 10.0,
    G: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    mass: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1000.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.px = px
    self.py = py
    self.G = G
    self.mass = mass

KerrSpace(mass=1.0, spin=0.5, *, curvature=False, label='space')

Bases: PseudoRiemannianSpace

Spacetime outside a spinning mass M of angular momentum M a, the Kerr metric in Boyer-Lindquist coordinates.

With Σ = r² + a² cos²θ and Δ = r² - 2 M r + a², the line element in (t, r, theta, phi) is ds² = -(1 - 2Mr/Σ) dt² - (4 M a r sin²θ / Σ) dt dφ + Σ/Δ dr² + Σ dθ² + (r² + a² + 2 M a² r sin²θ / Σ) sin²θ dφ², in units with G = c = 1. The term in dt dφ drags the frames of free particles around the spin axis. The mass and the spin enter the metric as exact rationals, which keeps the symbolic derivation fast.

Parameters:

Name Type Description Default
mass float

The mass M.

 1.0
spin float

The spin parameter a, the angular momentum per unit mass, between -M and M.

 0.5
curvature bool

Whether prepare derives the curvature tensors besides the Christoffel symbols, which takes long.

 False
label str

The label of the space in its system.

 'space'

Raises:

Type Description
ValidationError

The mass is not positive, or the spin exceeds the mass.

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import KerrSpace
>>> KerrSpace(mass=1.0, spin=0.5).prepare().metric_[0, 3]
-4*r*sin(theta)^2/(4*r^2 + cos(theta)^2)
Source code in src/skcomplex/systems/physics/celestial_mechanics/_geodesics.py 
162
163
164
165
166
167
168
169
170
171
172
def __init__(
    self,
    mass: float = 1.0,
    spin: float = 0.5,
    *,
    curvature: bool = False,
    label: str = 'space',
) -> None:
    super().__init__(curvature=curvature, label=label)
    self.mass = mass
    self.spin = spin

RelativisticOrbitElement(label, visual=None, *, u=1.5, u_d1=0.0, alpha=1.0, epsilon=0.05)

Bases: Element

Inverse radius of a planet along its orbit angle under the relativistic correction of Newtonian gravity.

The inverse radius u = 1 / r follows Einstein's orbit equation u'' + u = alpha + epsilon u², whose independent variable is the polar angle of the planet, which the simulation time stands for. For epsilon = 0 it is the Newtonian conic u = alpha (1 + e cos θ), and the quadratic term makes the perihelion advance every revolution.

Parameters:

Name Type Description Default
label str

The label of the planet.

 required
visual VisualSpec | None

How the planet is drawn, or None for default_visual_.

 None
u Annotated[float, StateSpec(order=2)]

The initial inverse radius.

 1.5
u_d1 float

The initial derivative of the inverse radius along the angle.

 0.0
alpha Annotated[float, ParameterSpec(domain=open(0, oo))]

The Newtonian inverse semi-latus rectum G M / h², with h the angular momentum per unit mass.

 1.0
epsilon Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The strength of the relativistic correction, 3 G M / c².

 0.05

Attributes:

Name Type Description
default_visual_

The visual of the planet when it sets none.

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import RelativisticOrbitElement
>>> RelativisticOrbitElement('planet').epsilon
0.05
Source code in src/skcomplex/systems/physics/celestial_mechanics/_spin_orbit.py 
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    u: Annotated[float, StateSpec(order=2)] = 1.5,
    u_d1: float = 0.0,
    alpha: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    epsilon: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
) -> None:
    super().__init__(label, visual=visual)
    self.u = u
    self.u_d1 = u_d1
    self.alpha = alpha
    self.epsilon = epsilon

RotatingFrameBodyElement(label, visual=None, *, x=0.5, y=0.0, px=0.0, py=0.5, mu=0.01215)

Bases: Element

Light body at (x, y) with the canonical momenta (px, py) in the frame that rotates with two primaries.

The states are of order 1, and the momenta are px = x' - y and py = y' + x for the unit angular velocity of the frame. The primaries of masses 1 - mu and mu sit at (-mu, 0) and (1 - mu, 0).

Parameters:

Name Type Description Default
label str

The label of the body.

 required
visual VisualSpec | None

How the body is drawn, or None for its class visual.

 None
x Annotated[float, StateSpec()]

The initial first coordinate.

 0.5
y Annotated[float, StateSpec()]

The initial second coordinate.

 0.0
px Annotated[float, StateSpec()]

The initial first momentum.

 0.0
py Annotated[float, StateSpec()]

The initial second momentum.

 0.5
mu Annotated[float, ParameterSpec(domain=open(0, oo))]

The mass of the lighter primary, as a fraction of the total.

 0.01215

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import RotatingFrameBodyElement
>>> RotatingFrameBodyElement('body', mu=0.1).mu
0.1
Source code in src/skcomplex/systems/physics/celestial_mechanics/_orbits.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.5,
    y: Annotated[float, StateSpec()] = 0.0,
    px: Annotated[float, StateSpec()] = 0.0,
    py: Annotated[float, StateSpec()] = 0.5,
    mu: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.01215,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.px = px
    self.py = py
    self.mu = mu

SchwarzschildSpace(mass=1.0, *, curvature=False, label='space')

Bases: PseudoRiemannianSpace

Spacetime outside a static, spherical mass M, the Schwarzschild metric in units with G = c = 1.

The line element in the coordinates (t, r, theta, phi) is ds² = -(1 - 2M/r) dt² + dr² / (1 - 2M/r) + r² (dθ² + sin²θ dφ²), with the horizon at r = 2M. The mass enters the metric as an exact rational, which keeps the symbolic derivation fast.

Parameters:

Name Type Description Default
mass float

The mass M.

 1.0
curvature bool

Whether prepare derives the curvature tensors besides the Christoffel symbols.

 False
label str

The label of the space in its system.

 'space'

Raises:

Type Description
ValidationError

The mass is not positive.

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import SchwarzschildSpace
>>> SchwarzschildSpace(mass=1.0, curvature=False).prepare().christoffel_[1][0][0]
(r - 2)/r^3
Source code in src/skcomplex/systems/physics/celestial_mechanics/_geodesics.py 
104
105
106
def __init__(self, mass: float = 1.0, *, curvature: bool = False, label: str = 'space') -> None:
    super().__init__(curvature=curvature, label=label)
    self.mass = mass

SitnikovBodyElement(label, visual=None, *, z=0.5, z_d1=0.0, e=0.0)

Bases: Element

Light body at the height z, of order 2, on the axis through the centre of mass of two equal primaries.

The primaries orbit with the eccentricity e.

Parameters:

Name Type Description Default
label str

The label of the body.

 required
visual VisualSpec | None

How the body is drawn, or None for none.

 None
z Annotated[float, StateSpec(order=2)]

The initial height.

 0.5
z_d1 float

The initial vertical velocity.

 0.0
e Annotated[float, ParameterSpec(domain=closed_open(0, oo))]

The eccentricity of the primaries' orbits.

 0.0

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import SitnikovBodyElement
>>> SitnikovBodyElement('body', e=0.1).e
0.1
Source code in src/skcomplex/systems/physics/celestial_mechanics/_orbits.py 
179
180
181
182
183
184
185
186
187
188
189
190
191
def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    z: Annotated[float, StateSpec(order=2)] = 0.5,
    z_d1: float = 0.0,
    e: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.z = z
    self.z_d1 = z_d1
    self.e = e

keplerian_dumbbell(u=1.0, vu=0.0, theta=0.0, omega=1.0, *, phi=0.0, spin=3.0, v=0.0, vv=0.0, psi=1.5707963267948966, vpsi=0.0, delta=0.5, epsilon=0.1, mu=0.9)

Return the Keplerian dumbbell: two masses on a rigid massless arm that orbit a primary under gravity.

The centre of mass of the dumbbell at the cylindrical coordinates (u, theta, v) about the primary and the arm at the spherical angles (phi, psi) move under the dimensionless Lagrangian L = (u'² + u² theta'² + v'²) / 2 + I (psi'² + sin² psi phi'²) / 2 + (1 - delta) / rho12 + delta / rho13, with the moment of inertia I = delta (1 - delta) epsilon² / mu of the arm. The distances of the two masses from the primary are rho12² = u² + v² + delta² epsilon² + 2 delta epsilon P and rho13² = u² + v² + (1 - delta)² epsilon² - 2 (1 - delta) epsilon P, with P = u sin psi cos(theta - phi) + v cos psi. Its Euler-Lagrange equations are the equations of motion of Dilão and Murteira. The gravity gradient across the arm couples the spin to the orbit, the simplest spin-orbit coupling of a rigid body. For epsilon = 0 the orbit is a Kepler orbit, of period 2π at u = 1. The energy and the angular momentum about the normal to the orbital plane are conserved. By default the arm spins three times faster than it orbits, in the orbital plane, and its spin and orbit exchange angular momentum.

Parameters:

Name Type Description Default
u float

The initial distance of the centre of mass from the axis through the primary.

 1.0
vu float

The initial radial velocity.

 0.0
theta float

The initial orbital angle.

 0.0
omega float

The initial orbital angular velocity.

 1.0
phi float

The initial azimuth of the arm.

 0.0
spin float

The initial spin of the arm about the normal to the orbital plane.

 3.0
v float

The initial height above the orbital plane.

 0.0
vv float

The initial vertical velocity.

 0.0
psi float

The initial polar angle of the arm, pi / 2 in the orbital plane.

 1.5707963267948966
vpsi float

The initial rate of the polar angle of the arm.

 0.0
delta float

The share of the second mass in the mass of the dumbbell, 1 / 2 for equal masses.

 0.5
epsilon float

The length of the arm in units of the mean orbit radius.

 0.1
mu float

The share of the primary in the mass of the three bodies.

 0.9

Returns:

Name Type Description
system System

The system, with the coordinates dumbbell.u, dumbbell.theta, dumbbell.v, dumbbell.phi, and dumbbell.psi of order 2, and the observables dumbbell.energy and dumbbell.angular_momentum.

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import keplerian_dumbbell
>>> keplerian_dumbbell(epsilon=0.2).elements[0].epsilon
0.2
Source code in src/skcomplex/systems/physics/celestial_mechanics/_spin_orbit.py 
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
@register_system(
    field='physics/celestial_mechanics',
    tags=('mechanics', 'ode', 'second order', 'gravity', 'spin-orbit coupling', 'lagrangian'),
    references=('Dilão and Murteira 2020', 'Dilão 2023'),
    defaults={
        'until': 60.0,
        'dt': 0.05,
        'solver': ScipyODESolver(method='DOP853', rtol=1e-10, atol=1e-12),
        'invariants': {'dumbbell.energy': 1e-6, 'dumbbell.angular_momentum': 1e-6},
    },
    observables=_get_dumbbell_observables,
)
def keplerian_dumbbell(
    u: float = 1.0,
    vu: float = 0.0,
    theta: float = 0.0,
    omega: float = 1.0,
    *,
    phi: float = 0.0,
    spin: float = 3.0,
    v: float = 0.0,
    vv: float = 0.0,
    psi: float = 1.5707963267948966,
    vpsi: float = 0.0,
    delta: float = 0.5,
    epsilon: float = 0.1,
    mu: float = 0.9,
) -> System:
    """Return the Keplerian dumbbell: two masses on a rigid massless arm that orbit a primary under gravity.

    The centre of mass of the dumbbell at the cylindrical coordinates `(u, theta, v)` about the primary and the arm at
    the spherical angles `(phi, psi)` move under the dimensionless Lagrangian
    `L = (u'² + u² theta'² + v'²) / 2 + I (psi'² + sin² psi phi'²) / 2 + (1 - delta) / rho12 + delta / rho13`, with the
    moment of inertia `I = delta (1 - delta) epsilon² / mu` of the arm. The distances of the two masses from the primary
    are `rho12² = u² + v² + delta² epsilon² + 2 delta epsilon P` and
    `rho13² = u² + v² + (1 - delta)² epsilon² - 2 (1 - delta) epsilon P`, with
    `P = u sin psi cos(theta - phi) + v cos psi`. Its Euler-Lagrange equations are the equations of motion of Dilão
    and Murteira. The gravity gradient across the arm couples the spin to the orbit, the simplest spin-orbit coupling
    of a rigid body. For `epsilon = 0` the orbit is a Kepler orbit, of period `2π` at `u = 1`. The energy and the
    angular momentum about the normal to the orbital plane are conserved. By default the arm spins three times faster
    than it orbits, in the orbital plane, and its spin and orbit exchange angular momentum.

    Args:
        u:
            The initial distance of the centre of mass from the axis through the primary.

        vu:
            The initial radial velocity.

        theta:
            The initial orbital angle.

        omega:
            The initial orbital angular velocity.

        phi:
            The initial azimuth of the arm.

        spin:
            The initial spin of the arm about the normal to the orbital plane.

        v:
            The initial height above the orbital plane.

        vv:
            The initial vertical velocity.

        psi:
            The initial polar angle of the arm, `pi / 2` in the orbital plane.

        vpsi:
            The initial rate of the polar angle of the arm.

        delta:
            The share of the second mass in the mass of the dumbbell, `1 / 2` for equal masses.

        epsilon:
            The length of the arm in units of the mean orbit radius.

        mu:
            The share of the primary in the mass of the three bodies.

    Returns:
        system:
            The system, with the coordinates `dumbbell.u`, `dumbbell.theta`, `dumbbell.v`, `dumbbell.phi`, and
            `dumbbell.psi` of order 2, and the observables `dumbbell.energy` and `dumbbell.angular_momentum`.

    Examples:
        >>> from skcomplex.systems.physics.celestial_mechanics import keplerian_dumbbell
        >>> keplerian_dumbbell(epsilon=0.2).elements[0].epsilon
        0.2
    """
    dumbbell = DumbbellElement(
        'dumbbell',
        u=u,
        u_d1=vu,
        theta=theta,
        theta_d1=omega,
        v=v,
        v_d1=vv,
        phi=phi,
        phi_d1=spin,
        psi=psi,
        psi_d1=vpsi,
        delta=delta,
        epsilon=epsilon,
        mu=mu,
    )
    motion = EulerLagrangeInteraction(
        'motion',
        target='dumbbell',
        states=COORDINATES,
        lagrangian=_get_dumbbell_lagrangian,
    )
    return System(elements=[dumbbell], interactions=[motion])

kerr_frame_dragging(radius=10.0, *, mass=1.0, spin=0.9)

Return a particle without angular momentum that falls into a spinning mass and is dragged around it.

The particle follows a timelike geodesic of the Kerr spacetime KerrSpace in the equatorial plane, in units with G = c = 1, with its proper time as the time of the system. It starts at rest in the radius at r₀, with zero angular momentum L = g_tφ dt/dτ + g_φφ dφ/dτ = 0. A particle without angular momentum around a static mass falls straight in, but the spin drags the frames around it: the particle turns with the angular velocity dφ/dt = -g_tφ / g_φφ = 2 M a / (r³ + a² r + 2 M a²) of the frames at its radius, faster as it falls.

Parameters:

Name Type Description Default
radius float

The initial radius r₀.

 10.0
mass float

The mass M of the central body.

 1.0
spin float

The spin parameter a, between -M and M.

 0.9

Returns:

Name Type Description
system System

The system, with the position particle.x = (t, r, theta, phi) and the four-velocity particle.x_d1.

Raises:

Type Description
ValidationError

The particle starts inside the ergosphere, at a radius of at most 2 M.

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import kerr_frame_dragging
>>> round(kerr_frame_dragging().elements[0].x_d1[3], 5)
0.00199
Source code in src/skcomplex/systems/physics/celestial_mechanics/_geodesics.py 
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
@register_system(
    field='physics/celestial_mechanics',
    tags=('general relativity', 'ode', 'second order', 'gravity', 'geodesics', 'frame dragging'),
    references=('Schwichtenberg 2020', 'Kerr 1963', 'Bardeen, Press, and Teukolsky 1972'),
    defaults={'until': 33.0, 'dt': 0.1, 'solver': GEODESIC_SOLVER},
)
def kerr_frame_dragging(
    radius: float = 10.0,
    *,
    mass: float = 1.0,
    spin: float = 0.9,
) -> System:
    """Return a particle without angular momentum that falls into a spinning mass and is dragged around it.

    The particle follows a timelike geodesic of the Kerr spacetime `KerrSpace` in the equatorial plane, in units with
    `G = c = 1`, with its proper time as the time of the system. It starts at rest in the radius at `r₀`, with zero
    angular momentum `L = g_tφ dt/dτ + g_φφ dφ/dτ = 0`. A particle without angular momentum around a static mass falls
    straight in, but the spin drags the frames around it: the particle turns with the angular velocity
    `dφ/dt = -g_tφ / g_φφ = 2 M a / (r³ + a² r + 2 M a²)` of the frames at its radius, faster as it falls.

    Args:
        radius:
            The initial radius `r₀`.

        mass:
            The mass `M` of the central body.

        spin:
            The spin parameter `a`, between `-M` and `M`.

    Returns:
        system:
            The system, with the position `particle.x = (t, r, theta, phi)` and the four-velocity `particle.x_d1`.

    Raises:
        ValidationError:
            The particle starts inside the ergosphere, at a radius of at most `2 M`.

    Examples:
        >>> from skcomplex.systems.physics.celestial_mechanics import kerr_frame_dragging
        >>> round(kerr_frame_dragging().elements[0].x_d1[3], 5)
        0.00199
    """
    if radius <= 2 * mass:
        error_msg = f'The radius {radius} is inside the ergosphere at 2 M = {2 * mass}. Give a larger radius.'
        raise ValidationError(error_msg)
    space = KerrSpace(mass=mass, spin=spin)
    g_tt = -(1 - 2 * mass / radius)
    g_tphi = -2 * mass * spin / radius
    g_phiphi = radius**2 + spin**2 + 2 * mass * spin**2 / radius
    omega = -g_tphi / g_phiphi
    rate = 1 / math.sqrt(-(g_tt + 2 * g_tphi * omega + g_phiphi * omega**2))
    particle = GeodesicParticleElement(
        'particle',
        x=[0.0, radius, EQUATOR, 0.0],
        x_d1=[rate, 0.0, 0.0, omega * rate],
    )
    geodesic = GeodesicInteraction('geodesic', target='particle', space=space)
    return System(elements=[particle], interactions=[geodesic], space=space)

light_bending(impact_parameter=20.0, distance=2000.0, *, mass=1.0)

Return a light ray that passes a static mass and is deflected towards it.

The ray follows a null geodesic of the Schwarzschild spacetime SchwarzschildSpace, in units with G = c = 1, with the affine parameter as the time of the system. It starts at the distance D along the axis behind the mass, offset by the impact parameter b, aimed parallel to the axis, with the energy E = 1 and the angular momentum L = b. Far from the mass the ray bends by 4 M / b, twice the Newtonian deflection of a particle at the speed of light, as the eclipse of 1919 measured for starlight that grazes the Sun.

Parameters:

Name Type Description Default
impact_parameter float

The impact parameter b, the distance of the unbent ray from the mass.

 20.0
distance float

The distance D along the axis from which the ray starts.

 2000.0
mass float

The mass M of the central body.

 1.0

Returns:

Name Type Description
system System

The system, with the position photon.x = (t, r, theta, phi) and the tangent photon.x_d1.

Raises:

Type Description
ValidationError

The ray would be captured by the mass, its impact parameter being at most 3 √3 M, or it starts closer than its impact parameter.

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import light_bending
>>> round(light_bending().elements[0].x_d1[0], 4)
1.001
Source code in src/skcomplex/systems/physics/celestial_mechanics/_geodesics.py 
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
@register_system(
    field='physics/celestial_mechanics',
    tags=('general relativity', 'ode', 'second order', 'gravity', 'geodesics', 'light'),
    references=('Schwichtenberg 2020', 'Satz 2022', 'Einstein 1915', 'Dyson, Eddington, and Davidson 1920'),
    defaults={'until': 4000.0, 'dt': 4.0, 'solver': GEODESIC_SOLVER},
)
def light_bending(
    impact_parameter: float = 20.0,
    distance: float = 2000.0,
    *,
    mass: float = 1.0,
) -> System:
    """Return a light ray that passes a static mass and is deflected towards it.

    The ray follows a null geodesic of the Schwarzschild spacetime `SchwarzschildSpace`, in units with `G = c = 1`, with
    the affine parameter as the time of the system. It starts at the distance `D` along the axis behind the mass,
    offset by the impact parameter `b`, aimed parallel to the axis, with the energy `E = 1` and the angular momentum
    `L = b`. Far from the mass the ray bends by `4 M / b`, twice the Newtonian deflection of a particle at the speed of
    light, as the eclipse of 1919 measured for starlight that grazes the Sun.

    Args:
        impact_parameter:
            The impact parameter `b`, the distance of the unbent ray from the mass.

        distance:
            The distance `D` along the axis from which the ray starts.

        mass:
            The mass `M` of the central body.

    Returns:
        system:
            The system, with the position `photon.x = (t, r, theta, phi)` and the tangent `photon.x_d1`.

    Raises:
        ValidationError:
            The ray would be captured by the mass, its impact parameter being at most `3 √3 M`, or it starts closer
            than its impact parameter.

    Examples:
        >>> from skcomplex.systems.physics.celestial_mechanics import light_bending
        >>> round(light_bending().elements[0].x_d1[0], 4)
        1.001
    """
    if impact_parameter <= 3 * math.sqrt(3) * mass:
        error_msg = (
            f'The impact parameter {impact_parameter} is at most 3 √3 M = {3 * math.sqrt(3) * mass}, so the mass '
            'captures the ray. Give a larger impact parameter.'
        )
        raise ValidationError(error_msg)
    if distance <= impact_parameter:
        error_msg = f'The distance {distance} is not larger than the impact parameter {impact_parameter}.'
        raise ValidationError(error_msg)
    radius = math.hypot(distance, impact_parameter)
    lapse = 1 - 2 * mass / radius
    radial = -math.sqrt(1 - lapse * impact_parameter**2 / radius**2)
    photon = GeodesicParticleElement(
        'photon',
        x=[0.0, radius, EQUATOR, math.atan2(impact_parameter, -distance)],
        x_d1=[1 / lapse, radial, 0.0, -impact_parameter / radius**2],
    )
    space = SchwarzschildSpace(mass=mass)
    geodesic = GeodesicInteraction('geodesic', target='photon', space=space)
    return System(elements=[photon], interactions=[geodesic], space=space)

n_body(preset='figure_eight', *, masses=None, positions=None, velocities=None, G=None)

Return bodies in the plane that attract each other pairwise by Newtonian gravity.

The preset 'figure_eight' is the periodic orbit of three equal masses found by Chenciner and Montgomery, with period 6.3259 for G = 1. The preset 'inner_solar_system' is the Sun with Mercury, Venus, the Earth, and Mars on circular orbits, in astronomical units and years, with G = 4 π². Masses, positions, velocities, or G given explicitly replace those of the preset.

Parameters:

Name Type Description Default
preset str

'figure_eight' or 'inner_solar_system'.

 'figure_eight'
masses Sequence[float] | None

The masses of the bodies, or None for those of the preset.

 None
positions Sequence[Sequence[float]] | None

The initial positions of the bodies, or None for those of the preset.

 None
velocities Sequence[Sequence[float]] | None

The initial velocities of the bodies, or None for those of the preset.

 None
G float | None

The gravitational constant, or None for that of the preset.

 None

Returns:

Name Type Description
system System

The system, with the particles body_<i> in two dimensions.

Raises:

Type Description
ValidationError

The preset is unknown.

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import n_body
>>> len(n_body('inner_solar_system').interactions)
10
Source code in src/skcomplex/systems/physics/celestial_mechanics/_orbits.py 
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
@register_system(
    field='physics/celestial_mechanics',
    tags=('mechanics', 'ode', 'gravity', 'chaos'),
    references=('Argun et al. 2021', 'Dilão 2023', 'Chenciner and Montgomery 2000'),
    defaults={'until': 6.33, 'dt': 0.01},
)
def n_body(
    preset: str = 'figure_eight',
    *,
    masses: Sequence[float] | None = None,
    positions: Sequence[Sequence[float]] | None = None,
    velocities: Sequence[Sequence[float]] | None = None,
    G: float | None = None,
) -> System:
    """Return bodies in the plane that attract each other pairwise by Newtonian gravity.

    The preset `'figure_eight'` is the periodic orbit of three equal masses found by Chenciner and Montgomery, with
    period 6.3259 for `G = 1`. The preset `'inner_solar_system'` is the Sun with Mercury, Venus, the Earth, and Mars on
    circular orbits, in astronomical units and years, with `G = 4 π²`. Masses, positions, velocities, or `G` given
    explicitly replace those of the preset.

    Args:
        preset:
            `'figure_eight'` or `'inner_solar_system'`.

        masses:
            The masses of the bodies, or `None` for those of the preset.

        positions:
            The initial positions of the bodies, or `None` for those of the preset.

        velocities:
            The initial velocities of the bodies, or `None` for those of the preset.

        G:
            The gravitational constant, or `None` for that of the preset.

    Returns:
        system:
            The system, with the particles `body_<i>` in two dimensions.

    Raises:
        ValidationError:
            The preset is unknown.

    Examples:
        >>> from skcomplex.systems.physics.celestial_mechanics import n_body
        >>> len(n_body('inner_solar_system').interactions)
        10
    """
    preset_G, preset_masses, preset_positions, preset_velocities = _preset(preset)
    masses = preset_masses if masses is None else masses
    positions = preset_positions if positions is None else positions
    velocities = preset_velocities if velocities is None else velocities
    G = preset_G if G is None else G
    bodies = [
        MassPointElement(f'body_{index}', x=list(position), x_d1=list(velocity), m=mass)
        for index, (mass, position, velocity) in enumerate(zip(masses, positions, velocities, strict=True))
    ]
    attractions = [
        GravityInteraction(
            f'gravity_{first}_{second}',
            first=f'body_{first}',
            second=f'body_{second}',
            G=G,
        )
        for first, second in itertools.combinations(range(len(bodies)), 2)
    ]
    return System(elements=bodies, interactions=attractions, space=ContinuousSpace(n_dim=2))

relativistic_orbit(u=1.5, du=0.0, *, alpha=1.0, epsilon=0.05)

Return Einstein's orbit equation u'' + u = alpha + epsilon u² of a planet, whose perihelion precesses.

The inverse radius u = 1 / r of the planet is a function of its polar angle, which the simulation time stands for, so the orbit in the plane is (cos t / u, sin t / u). For epsilon = 0 the orbit is a closed Newtonian ellipse of eccentricity |u(0) / alpha - 1| for u'(0) = 0. The term epsilon u² is the first relativistic correction, and with it the orbit is a rosette whose perihelion advances by about 2π epsilon alpha every revolution. The equation has two fixed points, the roots (1 ∓ sqrt(1 - 4 alpha epsilon)) / (2 epsilon) of epsilon u² - u + alpha. The smaller one, close to alpha, is a centre, the stable circular orbit, and the larger one is a saddle, the unstable circular orbit close to the black hole, beyond which a planet falls in.

Parameters:

Name Type Description Default
u float

The initial inverse radius.

 1.5
du float

The initial derivative of the inverse radius along the angle.

 0.0
alpha float

The Newtonian inverse semi-latus rectum.

 1.0
epsilon float

The strength of the relativistic correction.

 0.05

Returns:

Name Type Description
system System

The system, with the inverse radius planet.u of order 2 along the angle and its derivative planet.u_d1.

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import relativistic_orbit
>>> relativistic_orbit(epsilon=0.0).elements[0].epsilon
0.0
Source code in src/skcomplex/systems/physics/celestial_mechanics/_spin_orbit.py 
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
@register_system(
    field='physics/celestial_mechanics',
    tags=('mechanics', 'ode', 'second order', 'gravity', 'general relativity', 'precession'),
    references=('Einstein 1915', 'Dilão 2023'),
    defaults={'until': 60.0, 'dt': 0.02},
)
def relativistic_orbit(u: float = 1.5, du: float = 0.0, *, alpha: float = 1.0, epsilon: float = 0.05) -> System:
    """Return Einstein's orbit equation `u'' + u = alpha + epsilon u²` of a planet, whose perihelion precesses.

    The inverse radius `u = 1 / r` of the planet is a function of its polar angle, which the simulation time stands
    for, so the orbit in the plane is `(cos t / u, sin t / u)`. For `epsilon = 0` the orbit is a closed Newtonian
    ellipse of eccentricity `|u(0) / alpha - 1|` for `u'(0) = 0`. The term `epsilon u²` is the first relativistic
    correction, and with it the orbit is a rosette whose perihelion advances by about `2π epsilon alpha` every
    revolution. The equation has two fixed points, the roots `(1 ∓ sqrt(1 - 4 alpha epsilon)) / (2 epsilon)` of
    `epsilon u² - u + alpha`. The smaller one, close to `alpha`, is a centre, the stable circular orbit, and the
    larger one is a saddle, the unstable circular orbit close to the black hole, beyond which a planet falls in.

    Args:
        u:
            The initial inverse radius.

        du:
            The initial derivative of the inverse radius along the angle.

        alpha:
            The Newtonian inverse semi-latus rectum.

        epsilon:
            The strength of the relativistic correction.

    Returns:
        system:
            The system, with the inverse radius `planet.u` of order 2 along the angle and its derivative
            `planet.u_d1`.

    Examples:
        >>> from skcomplex.systems.physics.celestial_mechanics import relativistic_orbit
        >>> relativistic_orbit(epsilon=0.0).elements[0].epsilon
        0.0
    """
    return System(elements=[RelativisticOrbitElement('planet', u=u, u_d1=du, alpha=alpha, epsilon=epsilon)])

restricted_three_body(x=0.5, y=0.0, vx=0.0, vy=1.15, *, mu=0.01215)

Return a light body moving in the plane of two primaries on circular orbits, in the frame that rotates with them.

The primaries of masses 1 - μ and μ sit at (-μ, 0) and (1 - μ, 0), and the light body moves under the Hamiltonian H = (px² + py²) / 2 + y px - x py - (1 - μ) / r1 - μ / r2 of the rotating frame, with the momenta px = x' - y and py = y' + x. Its terms y px - x py couple the coordinates and the momenta, so the Hamiltonian is not separable. It equals minus half the Jacobi constant x² + y² + 2 (1 - μ) / r1 + 2 μ / r2 - x'² - y'², which is conserved. The default μ is that of the Earth and the Moon.

The default start has the Jacobi constant 2.835, below that of the neck towards the Moon, and its orbit is chaotic: it passes near the Moon after 3.5 time units and wanders between the two primaries, with a Lyapunov exponent near 0.2. Starts of a larger constant are confined around the Earth, such as vy = 0.8 with the constant 3.52, which traces a regular precessing rosette.

Parameters:

Name Type Description Default
x float

The initial first coordinate.

 0.5
y float

The initial second coordinate.

 0.0
vx float

The initial first velocity.

 0.0
vy float

The initial second velocity.

 1.15
mu float

The mass of the lighter primary, as a fraction of the total.

 0.01215

Returns:

Name Type Description
system System

The system, with the coordinates body.x and body.y and the momenta body.px and body.py.

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import restricted_three_body
>>> restricted_three_body(mu=0.1).elements[0].mu
0.1
Source code in src/skcomplex/systems/physics/celestial_mechanics/_orbits.py 
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
@register_system(
    field='physics/celestial_mechanics',
    tags=('mechanics', 'ode', 'gravity', 'chaos', 'rotating frame', 'hamiltonian'),
    references=('Dilão 2023', 'Szebehely 1967'),
    defaults={'until': 20.0, 'dt': 0.02, 'solver': RK4Solver(), 'invariants': {'body.energy': 1e-6}},
    observables=_get_rotating_frame_observables,
)
def restricted_three_body(
    x: float = 0.5,
    y: float = 0.0,
    vx: float = 0.0,
    vy: float = 1.15,
    *,
    mu: float = 0.01215,
) -> System:
    """Return a light body moving in the plane of two primaries on circular orbits, in the frame that rotates with them.

    The primaries of masses `1 - μ` and `μ` sit at `(-μ, 0)` and `(1 - μ, 0)`, and the light body moves under the
    Hamiltonian `H = (px² + py²) / 2 + y px - x py - (1 - μ) / r1 - μ / r2` of the rotating frame, with the momenta
    `px = x' - y` and `py = y' + x`. Its terms `y px - x py` couple the coordinates and the momenta, so the Hamiltonian
    is not separable. It equals minus half the Jacobi constant `x² + y² + 2 (1 - μ) / r1 + 2 μ / r2 - x'² - y'²`, which
    is conserved. The default `μ` is that of the Earth and the Moon.

    The default start has the Jacobi constant 2.835, below that of the neck towards the Moon, and its orbit is chaotic:
    it passes near the Moon after 3.5 time units and wanders between the two primaries, with a Lyapunov exponent near
    0.2. Starts of a larger constant are confined around the Earth, such as `vy = 0.8` with the constant 3.52, which
    traces a regular precessing rosette.

    Args:
        x:
            The initial first coordinate.

        y:
            The initial second coordinate.

        vx:
            The initial first velocity.

        vy:
            The initial second velocity.

        mu:
            The mass of the lighter primary, as a fraction of the total.

    Returns:
        system:
            The system, with the coordinates `body.x` and `body.y` and the momenta `body.px` and `body.py`.

    Examples:
        >>> from skcomplex.systems.physics.celestial_mechanics import restricted_three_body
        >>> restricted_three_body(mu=0.1).elements[0].mu
        0.1
    """
    return System(
        elements=[RotatingFrameBodyElement('body', x=x, y=y, px=vx - y, py=vy + x, mu=mu)],
        interactions=[
            HamiltonianInteraction(
                'motion',
                target='body',
                coordinates=('x', 'y'),
                momenta=('px', 'py'),
                hamiltonian=_get_rotating_frame_hamiltonian,
            ),
        ],
    )

schwarzschild_precession(semi_major_axis=100.0, eccentricity=0.5, *, mass=1.0)

Return a massive particle on a bound orbit around a static mass, whose perihelion advances every turn.

The particle follows a timelike geodesic of the Schwarzschild spacetime SchwarzschildSpace, in units with G = c = 1 and with its proper time as the time of the system. It starts at the periapsis r₁ = a (1 - e) of an orbit whose apoapsis is r₂ = a (1 + e), with the energy E and the angular momentum L for which both are turning points of E² = (1 - 2M/r) (1 + L²/r²) + (dr/dτ)². Unlike a Kepler ellipse, the orbit does not close: the periapsis advances by 6 π M / (a (1 - e²)) per turn when the orbit is far from the horizon, or 6 π G M / (c² a (1 - e²)) in other units, which accounts for the anomalous 43 arcseconds per century of Mercury.

Parameters:

Name Type Description Default
semi_major_axis float

The semi-major axis a of the orbit.

 100.0
eccentricity float

The eccentricity e of the orbit, in [0, 1).

 0.5
mass float

The mass M of the central body.

 1.0

Returns:

Name Type Description
system System

The system, with the position particle.x = (t, r, theta, phi) and the four-velocity particle.x_d1.

Raises:

Type Description
ValidationError

The eccentricity is outside [0, 1), or the orbit comes closer than the innermost stable circular orbit.

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import schwarzschild_precession
>>> round(schwarzschild_precession().elements[0].x[1], 6)
50.0
Source code in src/skcomplex/systems/physics/celestial_mechanics/_geodesics.py 
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
@register_system(
    field='physics/celestial_mechanics',
    tags=('general relativity', 'ode', 'second order', 'gravity', 'geodesics'),
    references=('Schwichtenberg 2020', 'Satz 2022', 'Schwarzschild 1916', 'Einstein 1915', 'Carroll 2004'),
    defaults={'until': 20000.0, 'dt': 10.0, 'solver': GEODESIC_SOLVER},
)
def schwarzschild_precession(
    semi_major_axis: float = 100.0,
    eccentricity: float = 0.5,
    *,
    mass: float = 1.0,
) -> System:
    """Return a massive particle on a bound orbit around a static mass, whose perihelion advances every turn.

    The particle follows a timelike geodesic of the Schwarzschild spacetime `SchwarzschildSpace`, in units with
    `G = c = 1` and with its proper time as the time of the system. It starts at the periapsis
    `r₁ = a (1 - e)` of an orbit whose apoapsis is `r₂ = a (1 + e)`, with the energy `E` and the angular momentum `L`
    for which both are turning points of `E² = (1 - 2M/r) (1 + L²/r²) + (dr/dτ)²`. Unlike a Kepler ellipse, the orbit
    does not close: the periapsis advances by `6 π M / (a (1 - e²))` per turn when the orbit is far from the
    horizon, or `6 π G M / (c² a (1 - e²))` in other units, which accounts for the anomalous 43 arcseconds per century
    of Mercury.

    Args:
        semi_major_axis:
            The semi-major axis `a` of the orbit.

        eccentricity:
            The eccentricity `e` of the orbit, in `[0, 1)`.

        mass:
            The mass `M` of the central body.

    Returns:
        system:
            The system, with the position `particle.x = (t, r, theta, phi)` and the four-velocity `particle.x_d1`.

    Raises:
        ValidationError:
            The eccentricity is outside `[0, 1)`, or the orbit comes closer than the innermost stable circular orbit.

    Examples:
        >>> from skcomplex.systems.physics.celestial_mechanics import schwarzschild_precession
        >>> round(schwarzschild_precession().elements[0].x[1], 6)
        50.0
    """
    if not 0 <= eccentricity < 1:
        error_msg = f'The eccentricity is {eccentricity!r}. Give a value in [0, 1) for a bound orbit.'
        raise ValidationError(error_msg)
    periapsis, apoapsis = semi_major_axis * (1 - eccentricity), semi_major_axis * (1 + eccentricity)
    if periapsis <= 6 * mass:
        error_msg = (
            f'The periapsis {periapsis} lies inside the innermost stable circular orbit at 6 M = {6 * mass}. Give a '
            'larger semi-major axis or a smaller eccentricity.'
        )
        raise ValidationError(error_msg)
    lapses = [1 - 2 * mass / radius for radius in (periapsis, apoapsis)]
    if eccentricity == 0:
        momentum_squared = mass * periapsis**2 / (periapsis - 3 * mass)
    else:
        momentum_squared = (lapses[1] - lapses[0]) / (lapses[0] / periapsis**2 - lapses[1] / apoapsis**2)
    energy = math.sqrt(lapses[0] * (1 + momentum_squared / periapsis**2))
    particle = GeodesicParticleElement(
        'particle',
        x=[0.0, periapsis, EQUATOR, 0.0],
        x_d1=[energy / lapses[0], 0.0, 0.0, math.sqrt(momentum_squared) / periapsis**2],
    )
    space = SchwarzschildSpace(mass=mass)
    geodesic = GeodesicInteraction('geodesic', target='particle', space=space)
    return System(elements=[particle], interactions=[geodesic], space=space)

sitnikov(z=0.5, vz=0.0, *, e=0.0)

Return a light body moving on the axis through the centre of mass of two equal primaries, normal to their orbit.

The primaries of mass 1/2 orbit their centre of mass with unit mean separation and period 2π, at distance r(t) = (1 - e cos t) / 2 from it to first order in the eccentricity e, and the light body follows z'' = -z / (z² + r(t)²)^(3/2). For circular primaries the energy z'²/2 - 1 / √(z² + 1/4) is conserved and the motion is integrable, and eccentric primaries make bounded orbits chaotic.

Parameters:

Name Type Description Default
z float

The initial height above the plane of the primaries.

 0.5
vz float

The initial vertical velocity.

 0.0
e float

The eccentricity of the primaries' orbits.

 0.0

Returns:

Name Type Description
system System

The system, with the height body.z of order 2.

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import sitnikov
>>> sitnikov(e=0.1).elements[0].e
0.1
Source code in src/skcomplex/systems/physics/celestial_mechanics/_orbits.py 
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
@register_system(
    field='physics/celestial_mechanics',
    tags=('mechanics', 'ode', 'second order', 'gravity', 'chaos'),
    references=('Dilão 2023', 'Sitnikov 1960'),
    defaults={'until': 50.0, 'dt': 0.05},
)
def sitnikov(z: float = 0.5, vz: float = 0.0, *, e: float = 0.0) -> System:
    """Return a light body moving on the axis through the centre of mass of two equal primaries, normal to their orbit.

    The primaries of mass 1/2 orbit their centre of mass with unit mean separation and period `2π`, at distance
    `r(t) = (1 - e cos t) / 2` from it to first order in the eccentricity `e`, and the light body follows
    `z'' = -z / (z² + r(t)²)^(3/2)`. For circular primaries the energy `z'²/2 - 1 / √(z² + 1/4)` is conserved and the
    motion is integrable, and eccentric primaries make bounded orbits chaotic.

    Args:
        z:
            The initial height above the plane of the primaries.

        vz:
            The initial vertical velocity.

        e:
            The eccentricity of the primaries' orbits.

    Returns:
        system:
            The system, with the height `body.z` of order 2.

    Examples:
        >>> from skcomplex.systems.physics.celestial_mechanics import sitnikov
        >>> sitnikov(e=0.1).elements[0].e
        0.1
    """
    return System(
        elements=[SitnikovBodyElement('body', z=z, z_d1=vz, e=e)],
        interactions=[GradientInteraction('gravity', target='body', states='z', potential=_get_sitnikov_potential)],
    )

stormer_problem(x=1.0, y=0.0, z=0.0, *, vx=0.0, vy=0.1, vz=0.1, alpha=1.0)

Return a charged particle moving in the magnetic dipole field of a planet, the Störmer problem.

The charge follows the Euler-Lagrange equations of |x'|² / 2 + alpha (x y' - y x') / r³, the kinetic energy and the coupling of the velocity to the vector potential of the dipole, whose Lorentz force gives x'' = 3 alpha z (y' z - z' y) / r⁵ - alpha y' / r³, y'' = -3 alpha z (x' z - z' x) / r⁵ + alpha x' / r³, and z'' = 3 alpha z (x' y - y' x) / r⁵, with alpha = q M / m for a dipole of moment M along the z axis. The magnetic force does no work, so the speed is conserved, and the symmetry about the axis conserves x y' - y x' + alpha (x² + y²) / r³. Bounded orbits near the planet trace the radiation belts, and many of them are chaotic.

Parameters:

Name Type Description Default
x float

The initial first coordinate, in units of the planet's radius.

 1.0
y float

The initial second coordinate.

 0.0
z float

The initial third coordinate, along the dipole.

 0.0
vx float

The initial first velocity.

 0.0
vy float

The initial second velocity.

 0.1
vz float

The initial third velocity.

 0.1
alpha float

The product of the charge and the dipole moment, divided by the mass.

 1.0

Returns:

Name Type Description
system System

The system, with the coordinates charge.x, charge.y, and charge.z, all of order 2.

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import stormer_problem
>>> stormer_problem(alpha=2.0).elements[0].alpha
2.0
Source code in src/skcomplex/systems/physics/celestial_mechanics/_dipoles.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
@register_system(
    field='physics/celestial_mechanics',
    tags=('mechanics', 'ode', 'second order', 'magnetism', 'chaos', 'hamiltonian'),
    references=('Störmer 1955', 'Dilão 2023'),
    defaults={'until': 20.0, 'dt': 0.02},
)
def stormer_problem(
    x: float = 1.0,
    y: float = 0.0,
    z: float = 0.0,
    *,
    vx: float = 0.0,
    vy: float = 0.1,
    vz: float = 0.1,
    alpha: float = 1.0,
) -> System:
    """Return a charged particle moving in the magnetic dipole field of a planet, the Störmer problem.

    The charge follows the Euler-Lagrange equations of `|x'|² / 2 + alpha (x y' - y x') / r³`, the kinetic energy and
    the coupling of the velocity to the vector potential of the dipole, whose Lorentz force gives
    `x'' = 3 alpha z (y' z - z' y) / r⁵ - alpha y' / r³`, `y'' = -3 alpha z (x' z - z' x) / r⁵ + alpha x' / r³`,
    and `z'' = 3 alpha z (x' y - y' x) / r⁵`, with `alpha = q M / m` for a dipole of moment `M` along the `z` axis.
    The magnetic force does no work, so the speed is conserved, and the symmetry about the axis conserves
    `x y' - y x' + alpha (x² + y²) / r³`. Bounded orbits near the planet trace the radiation belts, and many of them are
    chaotic.

    Args:
        x:
            The initial first coordinate, in units of the planet's radius.

        y:
            The initial second coordinate.

        z:
            The initial third coordinate, along the dipole.

        vx:
            The initial first velocity.

        vy:
            The initial second velocity.

        vz:
            The initial third velocity.

        alpha:
            The product of the charge and the dipole moment, divided by the mass.

    Returns:
        system:
            The system, with the coordinates `charge.x`, `charge.y`, and `charge.z`, all of order 2.

    Examples:
        >>> from skcomplex.systems.physics.celestial_mechanics import stormer_problem
        >>> stormer_problem(alpha=2.0).elements[0].alpha
        2.0
    """
    return System(
        elements=[DipoleChargeElement('charge', x=x, y=y, z=z, x_d1=vx, y_d1=vy, z_d1=vz, alpha=alpha)],
        interactions=[
            EulerLagrangeInteraction(
                'lorentz',
                target='charge',
                states=('x', 'y', 'z'),
                lagrangian=_get_dipole_lagrangian,
            ),
        ],
    )

two_body_orbit(eccentricity=0.0, *, G=1.0, mass=1000.0, radius=10.0)

Return a light planet orbiting a heavy star under Newtonian gravity, as the Kepler problem of their separation.

The planet moves in the coordinate (x, y) relative to the star with the canonical momenta (px, py) per unit mass, under the Hamiltonian H = (px² + py²) / 2 - G M / sqrt(x² + y²), which is separable. The energy and the angular momentum x py - y px are conserved, and a bound orbit is an ellipse with the star at a focus.

Parameters:

Name Type Description Default
eccentricity float

The eccentricity of the orbit, between 0 for a circle and 1 for an escape.

 0.0
G float

The gravitational constant.

 1.0
mass float

The mass of the star.

 1000.0
radius float

The distance of the planet from the star at its closest approach.

 10.0

Returns:

Name Type Description
system System

The system, with the coordinates planet.x and planet.y and the momenta planet.px and planet.py.

Examples:

>>> from skcomplex.systems.physics.celestial_mechanics import two_body_orbit
>>> [element.label for element in two_body_orbit().elements]
['planet']
Source code in src/skcomplex/systems/physics/celestial_mechanics/_orbits.py 
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
@register_system(
    field='physics/celestial_mechanics',
    tags=('mechanics', 'ode', 'gravity', 'hamiltonian'),
    references=('Kepler 1609', 'Newton 1687'),
    defaults={'until': 6.0, 'dt': 0.01, 'solver': Yoshida4Solver(), 'invariants': {'planet.energy': 1e-6}},
    observables=_get_kepler_observables,
)
def two_body_orbit(
    eccentricity: float = 0.0,
    *,
    G: float = 1.0,
    mass: float = 1000.0,
    radius: float = 10.0,
) -> System:
    """Return a light planet orbiting a heavy star under Newtonian gravity, as the Kepler problem of their separation.

    The planet moves in the coordinate `(x, y)` relative to the star with the canonical momenta `(px, py)` per unit
    mass, under the Hamiltonian `H = (px² + py²) / 2 - G M / sqrt(x² + y²)`, which is separable. The energy and the
    angular momentum `x py - y px` are conserved, and a bound orbit is an ellipse with the star at a focus.

    Args:
        eccentricity:
            The eccentricity of the orbit, between 0 for a circle and 1 for an escape.

        G:
            The gravitational constant.

        mass:
            The mass of the star.

        radius:
            The distance of the planet from the star at its closest approach.

    Returns:
        system:
            The system, with the coordinates `planet.x` and `planet.y` and the momenta `planet.px` and `planet.py`.

    Examples:
        >>> from skcomplex.systems.physics.celestial_mechanics import two_body_orbit
        >>> [element.label for element in two_body_orbit().elements]
        ['planet']
    """
    speed = float(np.sqrt(G * mass * (1 + eccentricity) / radius))
    return System(
        elements=[KeplerBodyElement('planet', x=radius, y=0.0, px=0.0, py=speed, G=G, mass=mass)],
        interactions=[
            HamiltonianInteraction(
                'gravity',
                target='planet',
                coordinates=('x', 'y'),
                momenta=('px', 'py'),
                hamiltonian=_get_kepler_hamiltonian,
            ),
        ],
    )