Restricted three-body problem
A light body moving in the plane of two primaries on circular orbits, in the frame that rotates with them.
- Name:
restricted_three_body - Field:
physics/celestial_mechanics - Formalism: ordinary differential equations
- References: Dilão 2023; Szebehely 1967
Factory
@register_model(
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 model, with the coordinates `body.x` and `body.y` and the momenta `body.px` and `body.py`.
Examples:
>>> from skcomplex.models.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,
),
],
)
Components
class RotatingFrameBodyElement(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)`.
Args:
label:
The label of the body.
visual:
How the body is drawn, or `None` for its class visual.
x:
The initial first coordinate.
y:
The initial second coordinate.
px:
The initial first momentum.
py:
The initial second momentum.
mu:
The mass of the lighter primary, as a fraction of the total.
Examples:
>>> from skcomplex.models.physics.celestial_mechanics import RotatingFrameBodyElement
>>> RotatingFrameBodyElement('body', mu=0.1).mu
0.1
"""
visual = VisualSpec('bodies', {'positions': ('x', 'y'), 'trail': None, 'marks': _primaries})
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
class HamiltonianInteraction(Interaction):
"""Hamilton's equations `q' = dH/dp` and `p' = -dH/dq` of the canonical pairs `(q, p)` of a Hamiltonian `H(ns)`.
The coordinates and the momenta are states of order 1, which may be angles, of every role that binds one element or
group, `target` by default, and the coordinates pair with the momenta in order. Every component of a vector state
is a coordinate or a momentum. The interaction returns one rate law per coordinate and per momentum, and the rates
of one state add, so the Hamiltonians of separate interactions sum: a coupling between two oscillators is an
interaction of its own whose Hamiltonian holds only the coupling energy. Dissipation or forcing is a rate
interaction on the momenta.
A Hamiltonian function takes the namespace of the interaction, with the namespaces of the roles, the time `t`, and
the symbols of the parameters a subclass declares.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose states are the canonical pairs.
coordinates:
The name of the state of order 1 that is a coordinate, or the names of several states.
momenta:
The name of the state of order 1 that is the momentum of the coordinate, or the names of several states in
the order of the coordinates.
hamiltonian:
A function of the namespace that returns the Hamiltonian.
Raises:
ValidationError:
A state is not a name, the coordinates and the momenta differ in count, the Hamiltonian is not callable,
or a coordinate and its momentum differ in size or are not of order 1.
Examples:
>>> from typing import Annotated
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.interactions import HamiltonianInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class OscillatorElement(Element):
... def __init__(self, label, q: Annotated[float, StateSpec()] = 1.0, p: Annotated[float, StateSpec()] = 0):
... super().__init__(label)
... self.q = q
... self.p = p
>>> energy = lambda s: (s.target.p**2 + s.target.q**2) / 2
>>> spring = HamiltonianInteraction('spring', target='mass', hamiltonian=energy)
>>> system = System(elements=[OscillatorElement('mass')], interactions=[spring])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=3.14159).trajectory_
>>> round(float(trajectory.values('mass.q')[-1]), 3)
-1.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
coordinates: str | Sequence[str] = 'q',
momenta: str | Sequence[str] = 'p',
hamiltonian: Callable[[SimpleNamespace], object] | None = None,
) -> None:
super().__init__(label)
self.target = target
self.coordinates = coordinates
self.momenta = momenta
self.hamiltonian = hamiltonian
def _get_names(self, names: str | Sequence[str]) -> tuple[str, ...]:
"""Return the names of the states of one side of the canonical pairs."""
return (names,) if isinstance(names, str) else tuple(names)
def _check(self) -> None:
"""Raise an error for a state that is not a name, an unpaired state, or a Hamiltonian that is not callable.
Raises:
ValidationError:
The type declares states, a state is not a name, the coordinates and the momenta differ in count, or
the Hamiltonian is not callable.
"""
super()._check()
for side in ('coordinates', 'momenta'):
names = self._get_names(getattr(self, side))
if not names or not all(isinstance(name, str) and name.isidentifier() for name in names):
error_msg = f'The {side} of {self.label!r} are {getattr(self, side)!r}. Give names of states.'
raise ValidationError(error_msg)
coordinates, momenta = self._get_names(self.coordinates), self._get_names(self.momenta)
if len(coordinates) != len(momenta):
unpaired = coordinates[len(momenta) :] or momenta[len(coordinates) :]
error_msg = (
f'{type(self).__name__} {self.label!r} pairs the coordinates {coordinates} with the momenta {momenta}, '
f'so {", ".join(unpaired)} has no partner. Give one momentum per coordinate.'
)
raise ValidationError(error_msg)
if self.hamiltonian is not None and not callable(self.hamiltonian):
error_msg = f'The hamiltonian of {self.label!r} is {self.hamiltonian!r}. Give a function of the namespace.'
raise ValidationError(error_msg)
def _get_hamiltonian(self, namespace: SimpleNamespace) -> object:
"""Return the Hamiltonian, the hook a subclass overrides to derive it from its parameters."""
if self.hamiltonian is None:
error_msg = f'{type(self).__name__} {self.label!r} has no Hamiltonian. Give a function of the namespace.'
raise ValidationError(error_msg)
return self.hamiltonian(namespace)
def _get_members(self, namespace: SimpleNamespace) -> list[SimpleNamespace]:
"""Return the namespaces of the members whose states are canonical pairs, those of every role of one path."""
return [
getattr(namespace, name)
for name, spec in self._get_specs().items()
if isinstance(spec, RoleSpec) and not spec.variadic and not spec.noise
]
def _get_state(self, member: SimpleNamespace, name: str) -> Target:
"""Return a canonical state of a member, after checking that it has order 1."""
if hasattr(member, f'{name}_d1'):
error_msg = (
f'{type(self).__name__} {self.label!r} reads the state {name} of order above 1, but canonical '
'coordinates and momenta have order 1. Declare it with `StateSpec(order=1)`, or derive the equations '
'of a coordinate of order 2 with `EulerLagrangeInteraction`.'
)
raise ValidationError(error_msg)
state: Target = getattr(member, name)
return state
def _get_pairs(self, namespace: SimpleNamespace) -> list[tuple[Target, Target]]:
"""Return the canonical pairs of every member, after checking that each coordinate matches its momentum."""
pairs = []
names = list(zip(self._get_names(self.coordinates), self._get_names(self.momenta), strict=True))
for member in self._get_members(namespace):
for coordinate, momentum in names:
q, p = self._get_state(member, coordinate), self._get_state(member, momentum)
if len(self._get_components(q)) != len(self._get_components(p)):
error_msg = (
f'The coordinate {coordinate} and the momentum {momentum} of {self.label!r} differ in size. '
'Give a momentum of the size of its coordinate.'
)
raise ValidationError(error_msg)
pairs.append((q, p))
return pairs
def _get_components(self, state: Target) -> list[Expression]:
"""Return the components of a state, the state itself for a scalar."""
return list(state) if isinstance(state, FreeModuleElement) else [state]
def _get_derivative(self, hamiltonian: Expression, state: Target) -> Target:
"""Return the derivative of the Hamiltonian with respect to a state, every component of a vector state."""
if isinstance(state, FreeModuleElement):
return vector(SR, [hamiltonian.diff(component) for component in state])
return hamiltonian.diff(state)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate of every coordinate, `dH/dp`, and of every momentum, `-dH/dq`.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the roles.
Returns:
laws:
The rate laws of the coordinates and the momenta.
"""
pairs = self._get_pairs(namespace)
hamiltonian = SR(self._get_hamiltonian(namespace))
laws: list[BasePart] = []
for q, p in pairs:
laws.append(RateLaw(q, self._get_derivative(hamiltonian, p)))
laws.append(RateLaw(p, -self._get_derivative(hamiltonian, q)))
return laws
Equations
body.x: diff(x_body(t), t) == px_body(t) + y_body(t)
body.y: diff(y_body(t), t) == py_body(t) - x_body(t)
body.px: diff(px_body(t), t) == (mu_body + x_body(t))*(mu_body - 1)/((mu_body + x_body(t))^2 + y_body(t)^2)^(3/2) - (mu_body + x_body(t) - 1)*mu_body/((mu_body + x_body(t) - 1)^2 + y_body(t)^2)^(3/2) + py_body(t)
body.py: diff(py_body(t), t) == (mu_body - 1)*y_body(t)/((mu_body + x_body(t))^2 + y_body(t)^2)^(3/2) - mu_body*y_body(t)/((mu_body + x_body(t) - 1)^2 + y_body(t)^2)^(3/2) - px_body(t)