Keplerian dumbbell and spin-orbit coupling
The Keplerian dumbbell: two masses on a rigid massless arm that orbit a primary under gravity.
- Name:
keplerian_dumbbell - Field:
physics/celestial_mechanics - Formalism: ordinary differential equations
- References: Dilão and Murteira 2020; Dilão 2023
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('keplerian_dumbbell')
entry = describe('keplerian_dumbbell')
trajectory = create_simulation('keplerian_dumbbell').prepare().run().trajectory_
Factory
@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])
Components
Elements:
DumbbellElement, fromElement
Interactions:
EulerLagrangeInteraction
class DumbbellElement(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.
Args:
label:
The label of the dumbbell.
visual:
How the dumbbell is drawn, or `None` for `default_visual_`.
u:
The initial distance from the axis through the primary.
u_d1:
The initial radial velocity.
theta:
The initial orbital angle.
theta_d1:
The initial orbital angular velocity.
v:
The initial height above the orbital plane.
v_d1:
The initial vertical velocity.
phi:
The initial azimuth of the arm.
phi_d1:
The initial spin of the arm about the vertical.
psi:
The initial polar angle of the arm, `pi / 2` in the orbital plane.
psi_d1:
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.
Attributes:
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
"""
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
class EulerLagrangeInteraction(Interaction):
"""Equations of motion `d/dt dL/dq' - dL/dq = Q - dR/dq'` of the coordinates `q` of a Lagrangian `L(ns)`.
The coordinates are states of order 2, and the Lagrangian, the Rayleigh dissipation function `R`, and the
generalized forcing `Q` read them and their first derivatives `derivative(q, 1)`. The interaction returns one
equation law per coordinate, which the compiler solves jointly for the second derivatives, so a Lagrangian with a
mass matrix, such as a double pendulum's, needs no inversion by hand. The equations of one coordinate from several
interactions add, so their Lagrangians, dissipations, and forcings sum: a coupling, such as a spring between two
pendulums, is an interaction of its own whose Lagrangian holds only the coupling energy. The coordinates are the
states of every role that binds one element or group, `target` by default, and every component of a vector state
is a coordinate.
A Lagrangian, dissipation, or forcing 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. The forcing returns one generalized
force per coordinate, in the order of the roles and of the states, or one expression for a single coordinate.
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 coordinates.
states:
The name of the state of order 2 that is a coordinate, or the names of several states.
lagrangian:
A function of the namespace that returns the Lagrangian.
dissipation:
A function of the namespace that returns the Rayleigh dissipation function, or `None` for none.
forcing:
A function of the namespace that returns the generalized forces, or `None` for none.
Raises:
ValidationError:
A state is not a name, or a function is not callable.
Examples:
>>> from typing import Annotated
>>> from sage.all__sagemath_symbolics import cos
>>> from skcomplex.core.components.elements import Element
>>> from skcomplex.core.components.interactions import EulerLagrangeInteraction
>>> from skcomplex.core.components.laws.utils import derivative
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.specs import StateSpec
>>> class AngleElement(Element):
... def __init__(self, label, q: Annotated[float, StateSpec(order=2)] = 0.1, q_d1: float = 0.0):
... super().__init__(label)
... self.q = q
... self.q_d1 = q_d1
>>> swing = EulerLagrangeInteraction(
... 'swing',
... target='bob',
... states='q',
... lagrangian=lambda s: derivative(s.target.q, 1) ** 2 / 2 + cos(s.target.q),
... )
>>> system = System(elements=[AngleElement('bob')], interactions=[swing])
>>> trajectory = Simulation(system, dt=0.01).prepare().run(until=3.2).trajectory_
>>> round(float(trajectory.values('bob.q')[-1]), 2)
-0.1
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (EquationLaw,)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('states',))],
states: str | Sequence[str] = 'x',
lagrangian: Callable[[SimpleNamespace], object] | None = None,
dissipation: Callable[[SimpleNamespace], object] | None = None,
forcing: Callable[[SimpleNamespace], object] | None = None,
) -> None:
super().__init__(label)
self.target = target
self.states = states
self.lagrangian = lagrangian
self.dissipation = dissipation
self.forcing = forcing
def _get_names(self) -> tuple[str, ...]:
"""Return the names of the states that are coordinates."""
return (self.states,) if isinstance(self.states, str) else tuple(self.states)
def _check(self) -> None:
"""Raise an error for a state that is not a name, or a function that is not callable.
Raises:
ValidationError:
The type declares states, a state is not a name, or a function is not callable.
"""
super()._check()
names = self._get_names()
if not names or not all(isinstance(name, str) and name.isidentifier() for name in names):
error_msg = f'{type(self).__name__} {self.label!r} reads the states {self.states!r}. Give names of states.'
raise ValidationError(error_msg)
for name in ('lagrangian', 'dissipation', 'forcing'):
value = getattr(self, name)
if value is not None and not callable(value):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace.'
raise ValidationError(error_msg)
def _get_lagrangian(self, namespace: SimpleNamespace) -> object:
"""Return the Lagrangian, the hook a subclass overrides to derive it from its parameters."""
if self.lagrangian is None:
error_msg = f'{type(self).__name__} {self.label!r} has no Lagrangian. Give a function of the namespace.'
raise ValidationError(error_msg)
return self.lagrangian(namespace)
def _get_coordinates(self, namespace: SimpleNamespace) -> list[Expression]:
"""Return the coordinates, every component of the states of every role that binds one path."""
coordinates = []
for name, spec in self._get_specs().items():
if not isinstance(spec, RoleSpec) or spec.variadic or spec.noise:
continue
for state in self._get_names():
value = getattr(getattr(namespace, name), state)
coordinates += list(value) if isinstance(value, FreeModuleElement) else [value]
return coordinates
def _get_forces(self, namespace: SimpleNamespace, count: int) -> list[Expression]:
"""Return the generalized force of every coordinate, zero without a forcing."""
if self.forcing is None:
return [SR(0)] * count
forces = self.forcing(namespace)
forces = list(forces) if isinstance(forces, FreeModuleElement | list | tuple) else [forces]
if len(forces) != count:
error_msg = f'The forcing of {self.label!r} returns {len(forces)} forces for {count} coordinates.'
raise ValidationError(error_msg)
return [SR(force) for force in forces]
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the Euler-Lagrange equation of every coordinate.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the roles.
Returns:
laws:
The equation laws of the coordinates.
"""
coordinates = self._get_coordinates(namespace)
velocities = {derivative(q, 1): SR.var(f'_velocity_{index}') for index, q in enumerate(coordinates)}
back = {symbol: item for item, symbol in velocities.items()}
slots = list(velocities.values())
lagrangian = SR(self._get_lagrangian(namespace)).subs(velocities)
dissipation = SR(0 if self.dissipation is None else self.dissipation(namespace)).subs(velocities)
forces = self._get_forces(namespace, len(coordinates))
time = SR(namespace.t)
laws: list[BasePart] = []
for q, slot, force in zip(coordinates, slots, forces, strict=True):
momentum = lagrangian.diff(slot)
change = momentum.diff(time) + sum(
momentum.diff(other) * other_slot + momentum.diff(other_slot) * derivative(other, 2)
for other, other_slot in zip(coordinates, slots, strict=True)
)
lhs = (change - lagrangian.diff(q)).subs(back)
rhs = (force - dissipation.diff(slot)).subs(back)
laws.append(EquationLaw(lhs == rhs, target=q))
return laws
Equations
dumbbell.u: diff(u_dumbbell(t), t, t) == u_dumbbell(t)*diff(theta_dumbbell(t), t)^2 + (delta_dumbbell*epsilon_dumbbell*cos(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t)) + u_dumbbell(t))*(delta_dumbbell - 1)/(delta_dumbbell^2*epsilon_dumbbell^2 + 2*(cos(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t))*u_dumbbell(t) + cos(psi_dumbbell(t))*v_dumbbell(t))*delta_dumbbell*epsilon_dumbbell + u_dumbbell(t)^2 + v_dumbbell(t)^2)^(3/2) - ((delta_dumbbell - 1)*epsilon_dumbbell*cos(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t)) + u_dumbbell(t))*delta_dumbbell/((delta_dumbbell - 1)^2*epsilon_dumbbell^2 + 2*(cos(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t))*u_dumbbell(t) + cos(psi_dumbbell(t))*v_dumbbell(t))*(delta_dumbbell - 1)*epsilon_dumbbell + u_dumbbell(t)^2 + v_dumbbell(t)^2)^(3/2)
dumbbell.theta: diff(theta_dumbbell(t), t, t) == ((delta_dumbbell - 1)*delta_dumbbell*epsilon_dumbbell*sin(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t))*u_dumbbell(t)/((delta_dumbbell - 1)^2*epsilon_dumbbell^2 + 2*(cos(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t))*u_dumbbell(t) + cos(psi_dumbbell(t))*v_dumbbell(t))*(delta_dumbbell - 1)*epsilon_dumbbell + u_dumbbell(t)^2 + v_dumbbell(t)^2)^(3/2) - (delta_dumbbell - 1)*delta_dumbbell*epsilon_dumbbell*sin(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t))*u_dumbbell(t)/(delta_dumbbell^2*epsilon_dumbbell^2 + 2*(cos(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t))*u_dumbbell(t) + cos(psi_dumbbell(t))*v_dumbbell(t))*delta_dumbbell*epsilon_dumbbell + u_dumbbell(t)^2 + v_dumbbell(t)^2)^(3/2) - 2*u_dumbbell(t)*diff(theta_dumbbell(t), t)*diff(u_dumbbell(t), t))/u_dumbbell(t)^2
dumbbell.v: diff(v_dumbbell(t), t, t) == (delta_dumbbell*epsilon_dumbbell*cos(psi_dumbbell(t)) + v_dumbbell(t))*(delta_dumbbell - 1)/(delta_dumbbell^2*epsilon_dumbbell^2 + 2*(cos(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t))*u_dumbbell(t) + cos(psi_dumbbell(t))*v_dumbbell(t))*delta_dumbbell*epsilon_dumbbell + u_dumbbell(t)^2 + v_dumbbell(t)^2)^(3/2) - ((delta_dumbbell - 1)*epsilon_dumbbell*cos(psi_dumbbell(t)) + v_dumbbell(t))*delta_dumbbell/((delta_dumbbell - 1)^2*epsilon_dumbbell^2 + 2*(cos(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t))*u_dumbbell(t) + cos(psi_dumbbell(t))*v_dumbbell(t))*(delta_dumbbell - 1)*epsilon_dumbbell + u_dumbbell(t)^2 + v_dumbbell(t)^2)^(3/2)
dumbbell.phi: diff(phi_dumbbell(t), t, t) == -(2*(delta_dumbbell - 1)*delta_dumbbell*epsilon_dumbbell^2*cos(psi_dumbbell(t))*sin(psi_dumbbell(t))*diff(phi_dumbbell(t), t)*diff(psi_dumbbell(t), t)/mu_dumbbell - (delta_dumbbell - 1)*delta_dumbbell*epsilon_dumbbell*sin(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t))*u_dumbbell(t)/((delta_dumbbell - 1)^2*epsilon_dumbbell^2 + 2*(cos(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t))*u_dumbbell(t) + cos(psi_dumbbell(t))*v_dumbbell(t))*(delta_dumbbell - 1)*epsilon_dumbbell + u_dumbbell(t)^2 + v_dumbbell(t)^2)^(3/2) + (delta_dumbbell - 1)*delta_dumbbell*epsilon_dumbbell*sin(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t))*u_dumbbell(t)/(delta_dumbbell^2*epsilon_dumbbell^2 + 2*(cos(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t))*u_dumbbell(t) + cos(psi_dumbbell(t))*v_dumbbell(t))*delta_dumbbell*epsilon_dumbbell + u_dumbbell(t)^2 + v_dumbbell(t)^2)^(3/2))*mu_dumbbell/((delta_dumbbell - 1)*delta_dumbbell*epsilon_dumbbell^2*sin(psi_dumbbell(t))^2)
dumbbell.psi: diff(psi_dumbbell(t), t, t) == ((delta_dumbbell - 1)*delta_dumbbell*epsilon_dumbbell^2*cos(psi_dumbbell(t))*sin(psi_dumbbell(t))*diff(phi_dumbbell(t), t)^2/mu_dumbbell + (cos(-phi_dumbbell(t) + theta_dumbbell(t))*cos(psi_dumbbell(t))*u_dumbbell(t) - sin(psi_dumbbell(t))*v_dumbbell(t))*(delta_dumbbell - 1)*delta_dumbbell*epsilon_dumbbell/((delta_dumbbell - 1)^2*epsilon_dumbbell^2 + 2*(cos(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t))*u_dumbbell(t) + cos(psi_dumbbell(t))*v_dumbbell(t))*(delta_dumbbell - 1)*epsilon_dumbbell + u_dumbbell(t)^2 + v_dumbbell(t)^2)^(3/2) - (cos(-phi_dumbbell(t) + theta_dumbbell(t))*cos(psi_dumbbell(t))*u_dumbbell(t) - sin(psi_dumbbell(t))*v_dumbbell(t))*(delta_dumbbell - 1)*delta_dumbbell*epsilon_dumbbell/(delta_dumbbell^2*epsilon_dumbbell^2 + 2*(cos(-phi_dumbbell(t) + theta_dumbbell(t))*sin(psi_dumbbell(t))*u_dumbbell(t) + cos(psi_dumbbell(t))*v_dumbbell(t))*delta_dumbbell*epsilon_dumbbell + u_dumbbell(t)^2 + v_dumbbell(t)^2)^(3/2))*mu_dumbbell/((delta_dumbbell - 1)*delta_dumbbell*epsilon_dumbbell^2)