Celestial mechanics

Celestial mechanics: orbits of bodies under gravity and charges in planetary 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.models.physics.celestial_mechanics import DipoleChargeElement
>>> DipoleChargeElement('charge', alpha=2.0).alpha
2.0

Source code in src/skcomplex/models/physics/celestial_mechanics/_dipoles.py

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

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.models.physics.celestial_mechanics import KeplerBodyElement
>>> KeplerBodyElement('planet', mass=10.0).mass
10.0

Source code in src/skcomplex/models/physics/celestial_mechanics/_kepler_body.py

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

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.models.physics.celestial_mechanics import RotatingFrameBodyElement
>>> RotatingFrameBodyElement('body', mu=0.1).mu
0.1

Source code in src/skcomplex/models/physics/celestial_mechanics/_rotating_frame_body.py

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

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.models.physics.celestial_mechanics import SitnikovBodyElement
>>> SitnikovBodyElement('body', e=0.1).e
0.1

Source code in src/skcomplex/models/physics/celestial_mechanics/_sitnikov_body.py

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

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 model, with the particles body_<i> in two dimensions.

Raises:

Type	Description
ValidationError	The preset is unknown.

Examples:

>>> from skcomplex.models.physics.celestial_mechanics import n_body
>>> len(n_body('inner_solar_system').interactions)
10

Source code in src/skcomplex/models/physics/celestial_mechanics/_orbits.py

@register_model(
    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 model, with the particles `body_<i>` in two dimensions.

    Raises:
        ValidationError:
            The preset is unknown.

    Examples:
        >>> from skcomplex.models.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))

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 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

Source code in src/skcomplex/models/physics/celestial_mechanics/_orbits.py

@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,
            ),
        ],
    )

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 model, with the height body.z of order 2.

Examples:

>>> from skcomplex.models.physics.celestial_mechanics import sitnikov
>>> sitnikov(e=0.1).elements[0].e
0.1

Source code in src/skcomplex/models/physics/celestial_mechanics/_orbits.py

@register_model(
    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 model, with the height `body.z` of order 2.

    Examples:
        >>> from skcomplex.models.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 model, with the coordinates charge.x, charge.y, and charge.z, all of order 2.

Examples:

>>> from skcomplex.models.physics.celestial_mechanics import stormer_problem
>>> stormer_problem(alpha=2.0).elements[0].alpha
2.0

Source code in src/skcomplex/models/physics/celestial_mechanics/_dipoles.py

@register_model(
    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 model, with the coordinates `charge.x`, `charge.y`, and `charge.z`, all of order 2.

    Examples:
        >>> from skcomplex.models.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 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']

Source code in src/skcomplex/models/physics/celestial_mechanics/_orbits.py

@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,
            ),
        ],
    )