Frame dragging around a spinning mass
A particle without angular momentum that falls into a spinning mass and is dragged around it.
- Name:
kerr_frame_dragging - Field:
physics/celestial_mechanics - Formalism: ordinary differential equations
- References: Schwichtenberg 2020; Kerr 1963; Bardeen, Press, and Teukolsky 1972
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('kerr_frame_dragging')
entry = describe('kerr_frame_dragging')
trajectory = create_simulation('kerr_frame_dragging').prepare().run().trajectory_
Factory
@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)
Components
Elements:
GeodesicParticleElement, fromElement
Interactions:
GeodesicInteraction
class GeodesicParticleElement(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)`.
Args:
label:
The label of the particle.
x:
The initial position `(t, r, theta, phi)`.
x_d1:
The initial four-velocity.
visual:
How the particle is drawn, or `None` for its class visual.
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]
"""
visual = VisualSpec('bodies', {'positions': _equatorial_position, 'trail': None, 'marks': _central_mass})
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
class GeodesicInteraction(Interaction):
"""Geodesic equation `x''^mu = -Γ^mu_(alpha beta)(x) x'^alpha x'^beta` of a point of a pseudo-Riemannian space.
The state is a vector of order 2 with one component per coordinate of the space, the position `x^mu` of a test
particle, and its companion holds the tangent `x'^mu`. The time of the system is the affine parameter of the
geodesic, the proper time of a timelike geodesic when the tangent starts with `g_(mu nu) x'^mu x'^nu = -1`. A
tangent that starts null, with `g_(mu nu) x'^mu x'^nu = 0`, follows a light ray. Either norm is conserved along the
geodesic, so it measures the error of a run. The Christoffel symbols are those the space derives from its metric.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state is the position.
state:
The name of the vector state of order 2 that is the position.
space:
The pseudo-Riemannian space whose geodesics the point follows, usually the space of the system too.
Attributes:
space_:
The prepared space, whose Christoffel symbols the rate reads.
Raises:
ValidationError:
The state is not a name.
SpaceError:
The space is not a pseudo-Riemannian space, raised by `prepare`.
Examples:
>>> from sage.all__sagemath_symbolics import SR, matrix
>>> from skcomplex.core.components.elements import PointElement
>>> from skcomplex.core.components.interactions import GeodesicInteraction
>>> from skcomplex.core.components.spaces import PseudoRiemannianSpace
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> r, phi = SR.var('r phi')
>>> plane = PseudoRiemannianSpace(matrix(SR, [[1, 0], [0, r**2]]), (r, phi))
>>> line = GeodesicInteraction('line', target='particle', space=plane)
>>> particle = PointElement('particle', x=[1.0, 0.0], x_d1=[0.0, 1.0])
>>> system = System(elements=[particle], interactions=[line], space=plane)
>>> trajectory = Simulation(system, dt=0.001).prepare().run(until=1.0).trajectory_
>>> [round(float(value), 3) for value in trajectory.values('particle.x')[-1]]
[1.414, 0.785]
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (RateLaw,)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('state',))],
state: str = 'x',
space: PseudoRiemannianSpace | None = None,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.space = space
def _check(self) -> None:
"""Raise an error for a state that is not a name.
Raises:
ValidationError:
The type declares states, or the state is not a name.
"""
super()._check()
if not isinstance(self.state, str) or not self.state.isidentifier():
error_msg = (
f'{type(self).__name__} {self.label!r} moves the state {self.state!r}. Give the name of a state.'
)
raise ValidationError(error_msg)
def _prepare(self, context: BaseSpace | None = None) -> None:
"""Set `space_`, the prepared pseudo-Riemannian space.
Args:
context:
The space the interaction inherits from its system.
Raises:
SpaceError:
The space is not a pseudo-Riemannian space.
"""
super()._prepare(context)
if not isinstance(self.space, PseudoRiemannianSpace):
error_msg = (
f'{type(self).__name__} {self.label!r} moves a point along the geodesics of {self.space!r}. Give a '
'`PseudoRiemannianSpace`.'
)
raise SpaceError(error_msg)
self.space_ = self.space.prepare()
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate law of the second derivative of the position, minus the connection on the tangent.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespace of the target.
Returns:
laws:
The rate law of the position.
Raises:
ValidationError:
The position has not one component per coordinate of the space.
"""
position = getattr(self._get_member(namespace, 'target'), self.state)
components = list(position) if isinstance(position, FreeModuleElement) else [position]
coordinates = self.space_.coordinates_
if len(components) != len(coordinates):
error_msg = (
f'The position {self.state!r} of {self.label!r} has {len(components)} components for the '
f'{len(coordinates)} coordinates of the space. Give a vector state of the dimension of the space.'
)
raise ValidationError(error_msg)
at = dict(zip(coordinates, components, strict=True))
tangent = [derivative(component, 1) for component in components]
indices = range(len(components))
rate = [
-sum(
(
SR(self.space_.christoffel_[mu][alpha][beta]).subs(at) * tangent[alpha] * tangent[beta]
for alpha in indices
for beta in indices
if not SR(self.space_.christoffel_[mu][alpha][beta]).is_zero()
),
SR(0),
)
for mu in indices
]
return [RateLaw(position, vector(SR, rate) if isinstance(position, FreeModuleElement) else rate[0])]
Equations
particle.x[0]: diff(x_particle_0(t), t, t) == -200*((8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2)*(200*x_particle_1(t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 - 1/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))/(8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) - 20000*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2 - 16200*x_particle_1(t)) - 16200*(200*sin(x_particle_2(t))^2*x_particle_1(t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 - sin(x_particle_2(t))^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))*x_particle_1(t)/(8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) - 20000*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2 - 16200*x_particle_1(t)))*diff(x_particle_0(t), t)*diff(x_particle_1(t), t) - 32400*(200*(81*cos(x_particle_2(t))*sin(x_particle_2(t))^3*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 + cos(x_particle_2(t))*sin(x_particle_2(t))*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))*x_particle_1(t)/(8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) - 20000*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2 - 16200*x_particle_1(t)) - (8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2)*cos(x_particle_2(t))*sin(x_particle_2(t))*x_particle_1(t)/((8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) - 20000*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2 - 16200*x_particle_1(t))*(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2))*diff(x_particle_0(t), t)*diff(x_particle_2(t), t) - 180*(200*(16200*sin(x_particle_2(t))^2*x_particle_1(t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 - 81*sin(x_particle_2(t))^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2) - x_particle_1(t))*sin(x_particle_2(t))^2*x_particle_1(t)/(8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) - 20000*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2 - 16200*x_particle_1(t)) - (8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2)*(200*sin(x_particle_2(t))^2*x_particle_1(t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 - sin(x_particle_2(t))^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))/(8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) - 20000*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2 - 16200*x_particle_1(t)))*diff(x_particle_1(t), t)*diff(x_particle_3(t), t) - 360*((8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2)*(81*cos(x_particle_2(t))*sin(x_particle_2(t))^3*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 + cos(x_particle_2(t))*sin(x_particle_2(t))*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))/(8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) - 20000*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2 - 16200*x_particle_1(t)) - ((16200*sin(x_particle_2(t))^2*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2) + 100*x_particle_1(t)^2 + 81)*cos(x_particle_2(t))*sin(x_particle_2(t)) + 16200*(81*cos(x_particle_2(t))*sin(x_particle_2(t))^3*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 + cos(x_particle_2(t))*sin(x_particle_2(t))*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))*sin(x_particle_2(t))^2)*x_particle_1(t)/(8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) - 20000*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2 - 16200*x_particle_1(t)))*diff(x_particle_2(t), t)*diff(x_particle_3(t), t)
particle.x[1]: diff(x_particle_1(t), t, t) == -(16200*sin(x_particle_2(t))^2*x_particle_1(t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 - 81*sin(x_particle_2(t))^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2) - x_particle_1(t))*(100*x_particle_1(t)^2 - 200*x_particle_1(t) + 81)*sin(x_particle_2(t))^2*diff(x_particle_3(t), t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2) - 100*(100*x_particle_1(t)^2 - 200*x_particle_1(t) + 81)*(200*x_particle_1(t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 - 1/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))*diff(x_particle_0(t), t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2) + 100*(100*x_particle_1(t)^2 - 200*x_particle_1(t) + 81)*((81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)*(x_particle_1(t) - 1)/(100*x_particle_1(t)^2 - 200*x_particle_1(t) + 81)^2 - x_particle_1(t)/(100*x_particle_1(t)^2 - 200*x_particle_1(t) + 81))*diff(x_particle_1(t), t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2) + 162*cos(x_particle_2(t))*sin(x_particle_2(t))*diff(x_particle_1(t), t)*diff(x_particle_2(t), t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2) + (100*x_particle_1(t)^2 - 200*x_particle_1(t) + 81)*x_particle_1(t)*diff(x_particle_2(t), t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2) + 180*(200*sin(x_particle_2(t))^2*x_particle_1(t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 - sin(x_particle_2(t))^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))*(100*x_particle_1(t)^2 - 200*x_particle_1(t) + 81)*diff(x_particle_0(t), t)*diff(x_particle_3(t), t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)
particle.x[2]: diff(x_particle_2(t), t, t) == 81*cos(x_particle_2(t))*sin(x_particle_2(t))*diff(x_particle_2(t), t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2) + 1620000*cos(x_particle_2(t))*sin(x_particle_2(t))*x_particle_1(t)*diff(x_particle_0(t), t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^3 - 8100*cos(x_particle_2(t))*sin(x_particle_2(t))*diff(x_particle_1(t), t)^2/((81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)*(100*x_particle_1(t)^2 - 200*x_particle_1(t) + 81)) - 200*x_particle_1(t)*diff(x_particle_1(t), t)*diff(x_particle_2(t), t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2) - 36000*(81*cos(x_particle_2(t))*sin(x_particle_2(t))^3*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 + cos(x_particle_2(t))*sin(x_particle_2(t))*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))*diff(x_particle_0(t), t)*diff(x_particle_3(t), t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2) + ((16200*sin(x_particle_2(t))^2*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2) + 100*x_particle_1(t)^2 + 81)*cos(x_particle_2(t))*sin(x_particle_2(t)) + 16200*(81*cos(x_particle_2(t))*sin(x_particle_2(t))^3*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 + cos(x_particle_2(t))*sin(x_particle_2(t))*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))*sin(x_particle_2(t))^2)*diff(x_particle_3(t), t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)
particle.x[3]: diff(x_particle_3(t), t, t) == -18000*((81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2 - 200*x_particle_1(t))*(200*sin(x_particle_2(t))^2*x_particle_1(t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 - sin(x_particle_2(t))^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))/(8100*cos(x_particle_2(t))^2*sin(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*sin(x_particle_2(t))^2*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^4*x_particle_1(t) - 20000*sin(x_particle_2(t))^2*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2*sin(x_particle_2(t))^2 + 8100*sin(x_particle_2(t))^2*x_particle_1(t)^2 - 16200*sin(x_particle_2(t))^2*x_particle_1(t)) + 200*(200*x_particle_1(t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 - 1/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))*x_particle_1(t)/(8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) - 20000*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2 - 16200*x_particle_1(t)))*diff(x_particle_0(t), t)*diff(x_particle_1(t), t) + 36000*((81*cos(x_particle_2(t))*sin(x_particle_2(t))^3*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 + cos(x_particle_2(t))*sin(x_particle_2(t))*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))*(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2 - 200*x_particle_1(t))/(8100*cos(x_particle_2(t))^2*sin(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*sin(x_particle_2(t))^2*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^4*x_particle_1(t) - 20000*sin(x_particle_2(t))^2*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2*sin(x_particle_2(t))^2 + 8100*sin(x_particle_2(t))^2*x_particle_1(t)^2 - 16200*sin(x_particle_2(t))^2*x_particle_1(t)) + 16200*cos(x_particle_2(t))*sin(x_particle_2(t))*x_particle_1(t)^2/((8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) - 20000*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2 - 16200*x_particle_1(t))*(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2))*diff(x_particle_0(t), t)*diff(x_particle_2(t), t) + 200*((81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2 - 200*x_particle_1(t))*(16200*sin(x_particle_2(t))^2*x_particle_1(t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 - 81*sin(x_particle_2(t))^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2) - x_particle_1(t))*sin(x_particle_2(t))^2/(8100*cos(x_particle_2(t))^2*sin(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*sin(x_particle_2(t))^2*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^4*x_particle_1(t) - 20000*sin(x_particle_2(t))^2*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2*sin(x_particle_2(t))^2 + 8100*sin(x_particle_2(t))^2*x_particle_1(t)^2 - 16200*sin(x_particle_2(t))^2*x_particle_1(t)) + 16200*(200*sin(x_particle_2(t))^2*x_particle_1(t)^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 - sin(x_particle_2(t))^2/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))*x_particle_1(t)/(8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) - 20000*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2 - 16200*x_particle_1(t)))*diff(x_particle_1(t), t)*diff(x_particle_3(t), t) - 2*(((16200*sin(x_particle_2(t))^2*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2) + 100*x_particle_1(t)^2 + 81)*cos(x_particle_2(t))*sin(x_particle_2(t)) + 16200*(81*cos(x_particle_2(t))*sin(x_particle_2(t))^3*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 + cos(x_particle_2(t))*sin(x_particle_2(t))*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))*sin(x_particle_2(t))^2)*(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2 - 200*x_particle_1(t))/(8100*cos(x_particle_2(t))^2*sin(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*sin(x_particle_2(t))^2*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^4*x_particle_1(t) - 20000*sin(x_particle_2(t))^2*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2*sin(x_particle_2(t))^2 + 8100*sin(x_particle_2(t))^2*x_particle_1(t)^2 - 16200*sin(x_particle_2(t))^2*x_particle_1(t)) + 3240000*(81*cos(x_particle_2(t))*sin(x_particle_2(t))^3*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2)^2 + cos(x_particle_2(t))*sin(x_particle_2(t))*x_particle_1(t)/(81*cos(x_particle_2(t))^2 + 100*x_particle_1(t)^2))*x_particle_1(t)/(8100*cos(x_particle_2(t))^2*x_particle_1(t)^2 + 10000*x_particle_1(t)^4 + 16200*sin(x_particle_2(t))^2*x_particle_1(t) - 20000*x_particle_1(t)^3 + 6561*cos(x_particle_2(t))^2 + 8100*x_particle_1(t)^2 - 16200*x_particle_1(t)))*diff(x_particle_2(t), t)*diff(x_particle_3(t), t)