Central Forces & Orbital Mechanics

The central forces module provides specialized tools for analyzing motion in central force fields, including orbital mechanics and the Kepler problem.

Overview

A central force depends only on the distance from a fixed point:

\[\mathbf{F} = F(r) \hat{r}\]

The module implements:

  • Effective Potential: Reduce 2D motion to 1D radial problem

  • Orbit Classification: Bounded, unbounded, circular, elliptical

  • Turning Points: Find perihelion and aphelion

  • Kepler Problem: Specialized solver for gravitational orbits

  • Precession: Calculate apsidal precession for non-Keplerian potentials

Theory

Effective Potential

For a particle with angular momentum \(L\) in a central potential \(V(r)\):

\[V_{\text{eff}}(r) = V(r) + \frac{L^2}{2mr^2}\]

The radial motion is equivalent to 1D motion in this effective potential.

Orbit Classification

Given energy \(E\) and effective potential \(V_{\text{eff}}\):

  • Bounded: \(E < V_{\text{eff}}(\infty)\), motion confined between turning points

  • Unbounded: \(E > V_{\text{eff}}(\infty)\), particle escapes to infinity

  • Circular: \(E = V_{\text{eff,min}}\), particle stays at fixed radius

Kepler Orbits

For gravitational potential \(V = -k/r\):

\[r(\phi) = \frac{p}{1 + e \cos(\phi - \phi_0)}\]

where:

  • \(p = L^2/(mk)\) is the semi-latus rectum

  • \(e = \sqrt{1 + 2EL^2/(mk^2)}\) is the eccentricity

Usage Examples

Effective Potential Analysis

from mechanics_dsl.domains.classical import EffectivePotential
import sympy as sp

# Create effective potential analyzer
eff = EffectivePotential()

# Gravitational potential
r = sp.Symbol('r', positive=True)
k = sp.Symbol('k', positive=True)
m = sp.Symbol('m', positive=True)
L = sp.Symbol('L', positive=True)

V = -k / r

V_eff = eff.compute(V, L, m, 'r')
# V_eff = -k/r + L²/(2mr²)

# Find minimum (circular orbit radius)
r_circular = eff.find_minimum(V_eff, 'r')
# r_circular = L² / (mk)

Finding Turning Points

from mechanics_dsl.domains.classical import CentralForceAnalyzer
import numpy as np

analyzer = CentralForceAnalyzer()

# Define potential and parameters
def V(r):
    return -1.0 / r  # Gravitational (k=1)

# Find turning points for given E and L
E = -0.5  # Negative = bound orbit
L = 1.0
m = 1.0

turning_points = analyzer.find_turning_points(V, E, L, m)

print(f"Perihelion: r_min = {turning_points.r_min:.4f}")
print(f"Aphelion: r_max = {turning_points.r_max:.4f}")

Orbit Classification

from mechanics_dsl.domains.classical import CentralForceAnalyzer, OrbitType

analyzer = CentralForceAnalyzer()

def V(r):
    return -1.0 / r

# Bound orbit (E < 0)
orbit_type = analyzer.classify_orbit(V, E=-0.5, L=1.0, m=1.0)
print(orbit_type)  # OrbitType.BOUNDED

# Unbound orbit (E > 0)
orbit_type = analyzer.classify_orbit(V, E=0.5, L=1.0, m=1.0)
print(orbit_type)  # OrbitType.UNBOUNDED

# Parabolic (E = 0)
orbit_type = analyzer.classify_orbit(V, E=0.0, L=1.0, m=1.0)
print(orbit_type)  # OrbitType.PARABOLIC

Kepler Problem Solver

from mechanics_dsl.domains.classical import KeplerProblem

# Create Kepler problem
kepler = KeplerProblem(G=1.0, M=1.0, m=1.0)

# From initial conditions
r0 = 1.0
v0 = 1.0

elements = kepler.compute_orbital_elements(r0, v0, phi0=0)

print(f"Semi-major axis: a = {elements.semi_major_axis:.4f}")
print(f"Eccentricity: e = {elements.eccentricity:.4f}")
print(f"Period: T = {elements.period:.4f}")

# Compute position at time t
t = 1.0
r, phi = kepler.position_at_time(t, elements)

Computing Orbital Period

from mechanics_dsl.domains.classical import KeplerProblem
import numpy as np

kepler = KeplerProblem(G=6.674e-11, M=1.989e30, m=5.972e24)  # Sun-Earth

# Kepler's third law: T² = (4π²/GM) a³
a = 1.496e11  # 1 AU in meters

T = kepler.compute_period(a)
print(f"Orbital period: {T / (365.25 * 24 * 3600):.4f} years")
# Output: ~1.0 year

API Reference

Classes

class EffectivePotential

Compute and analyze effective potentials.

compute(V, L, m, r_var)

Compute \(V_{\text{eff}} = V + L^2/(2mr^2)\).

find_minimum(V_eff, r_var)

Find radius of minimum (circular orbit).

class CentralForceAnalyzer

Analyzer for central force motion.

find_turning_points(V, E, L, m)

Find radii where \(E = V_{\text{eff}}(r)\).

Returns:

TurningPoints namedtuple with r_min, r_max

classify_orbit(V, E, L, m)

Classify orbit based on energy and effective potential.

Returns:

OrbitType enum

class KeplerProblem

Specialized solver for Keplerian (gravitational) orbits.

compute_orbital_elements(r, v, phi)

Compute orbital elements from position and velocity.

Returns:

OrbitalElements namedtuple

position_at_time(t, elements)

Compute position using Kepler’s equation.

compute_period(a)

Compute orbital period from semi-major axis.

Enums

class OrbitType

  • BOUNDED: Elliptical orbit (E < 0)

  • UNBOUNDED: Hyperbolic orbit (E > 0)

  • PARABOLIC: Escape trajectory (E = 0)

  • CIRCULAR: Circular orbit (E = V_eff_min)

  • COLLISION: Orbit intersects origin

See Also