Note
Go to the end to download the full example code.
Interactive Plotly visualizers: SO(3) momentum sphere and orbit#
Unlike the static Matplotlib plots elsewhere in
physicskit.classical.visualizers, these are draggable/zoomable in a
live Python session or a browser – useful for directly inspecting how
close a rigid-body trajectory passes to the unstable intermediate-axis
pole, or for zooming into a precessing orbit’s perihelion. The two
systems are governed by the equations already introduced elsewhere in
the gallery:
EulerTop’s torque-free
Euler equations,
\(I_i\dot\omega_i = (I_j - I_k)\,\omega_j\omega_k\) (cyclic in
\(i,j,k\)) for the body-frame angular velocity (see
The intermediate axis theorem (Dzhanibekov effect)), and
KeplerSystem’s
separable Hamiltonian \(H = |p|^2/(2\mu) + V(r)\) with the
post-Newtonian perturbation \(c_\mathrm{pn}/r^3\) in \(V(r)\)
that drives the perihelion precession (see
Kepler orbits and perihelion precession).
from physicskit.classical.core.base_system import SimulationResult
from physicskit.classical.systems.newtonian import KeplerSystem
from physicskit.classical.systems.rotations import EulerTop
from physicskit.classical.visualizers.interactive import interactive_orbit, interactive_so3_momentum_sphere
# The integrations below use a fine dt for numerical accuracy, but an
# interactive plot of every single step (potentially hundreds of
# thousands of points) is both unnecessarily large and sluggish to pan
# and zoom -- downsample with a stride to a few hundred points, plenty
# to see the geometry clearly.
def _stride(result: SimulationResult, stride: int) -> SimulationResult:
idx = list(range(0, len(result.t), stride))
if idx[-1] != len(result.t) - 1:
idx.append(len(result.t) - 1)
q = result.q[idx] if result.q is not None else None
p = result.p[idx] if result.p is not None else None
return SimulationResult(t=result.t[idx], y=result.y[idx], q=q, p=p)
SO(3) momentum sphere for a tumbling (intermediate-axis) Euler top#
top = EulerTop([0.01, 1.0, 0.01], I1=1.0, I2=2.0, I3=3.0)
result_top = _stride(top.integrate((0, 100), dt=1e-1, method="implicit_midpoint"), stride=200)
fig_sphere = interactive_so3_momentum_sphere(result_top.y[:, :3], I1=1.0, I2=2.0, I3=3.0, title="EulerTop: tumbling about the intermediate axis")
A precessing Kepler orbit with its LRL vector rotating#
kepler = KeplerSystem.from_orbital_elements(a=1.0, e=0.3, c_pn=0.02)
result_orbit = _stride(kepler.integrate((0, 300), dt=1e-3, method="yoshida4"), stride=300)
fig_orbit = interactive_orbit(result_orbit, system=kepler, title="Precessing Kepler orbit (c_pn=0.02)")
Total running time of the script: (0 minutes 0.434 seconds)