Two-body orbit
A light planet orbiting a heavy star under Newtonian gravity, as the Kepler problem of their separation.
- Name:
two_body_orbit - Field:
physics/celestial_mechanics - Formalism: ordinary differential equations
- References: Kepler 1609; Newton 1687
Factory
@register_model(
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 model, with the coordinates `planet.x` and `planet.y` and the momenta `planet.px` and `planet.py`.
Examples:
>>> from skcomplex.models.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,
),
],
)
Components
class KeplerBodyElement(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.
Args:
label:
The label of the planet.
visual:
How the planet is drawn, or `None` for its class visual.
x:
The initial first coordinate.
y:
The initial second coordinate.
px:
The initial first momentum per unit mass.
py:
The initial second momentum per unit mass.
G:
The gravitational constant.
mass:
The mass of the star.
Examples:
>>> from skcomplex.models.physics.celestial_mechanics import KeplerBodyElement
>>> KeplerBodyElement('planet', mass=10.0).mass
10.0
"""
visual = VisualSpec('bodies', {'positions': ('x', 'y'), 'trail': None, 'marks': _star})
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
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
planet.x: diff(x_planet(t), t) == px_planet(t)
planet.y: diff(y_planet(t), t) == py_planet(t)
planet.px: diff(px_planet(t), t) == -G_planet*mass_planet*x_planet(t)/(x_planet(t)^2 + y_planet(t)^2)^(3/2)
planet.py: diff(py_planet(t), t) == -G_planet*mass_planet*y_planet(t)/(x_planet(t)^2 + y_planet(t)^2)^(3/2)