A Python framework for simulating spacecraft dynamics near small bodies in the photogravitational circular restricted three-body problem (PGCR3BP), built to test station-keeping controllers at SRP-balanced artificial equilibrium points.
PyAsteroid generates randomized asteroid populations, computes their irregular gravity fields from polyhedral meshes, models solar radiation pressure including self-shadowing, locates artificial equilibrium points, and provides an interactive 3D visualization layer for inspecting the result.
- Motivation
- The Physics
- Installation
- Quick Start
- Module Reference
- Monte Carlo Campaigns
- Validation
- Regenerating the Figures
- Known Limitations
- References
Close-proximity operations at small bodies are hard for a specific reason: the asteroid's gravity is weak enough that non-gravitational forces are not perturbations, they are comparable to the primary attraction. Solar radiation pressure on a spacecraft near a kilometer-scale asteroid can rival the body's own gravitational pull.
That flips from a nuisance into a design opportunity. If SRP is deliberately balanced against asteroid gravity, a spacecraft can be placed at an artificial equilibrium point — a location where, in the rotating frame, the net force vanishes. The spacecraft co-orbits the Sun alongside the asteroid rather than orbiting the asteroid itself, holding proximity for science return without expending propellant on orbital maintenance.
The catch is that these equilibria are saddles. They are unstable, and the environment that creates them is not steady: the gravity field is irregular and time-varying as the body rotates, and the solar flux driving SRP fluctuates with solar activity. Holding position requires active disturbance rejection.
PyAsteroid exists to build that environment so controllers can be tested against it.
The framework uses the nondimensional circular restricted three-body problem with the Sun as primary and the asteroid as secondary, modified so that solar radiation pressure reduces the effective solar attraction.
Working in the rotating frame with the origin at the Sun–asteroid barycenter, distances normalized by the heliocentric separation
where the effective potential is
with
Here
The single modification relative to the classical CR3BP is the factor
Coriolis is not a disturbance. It depends only on measured velocity and the known frame rotation rate, making it exactly computable. In the controller comparison it belongs in the feedforward term, not in the disturbance to be rejected.
Both forces fall off as
This makes
An equilibrium requires
Let
where
The relative size of the first two terms decides the regime. Their ratio is
For kilometer-scale bodies at ~1 AU this is of order
PyAsteroid uses this as the initial guess, then refines it by Newton iteration against the full potential — including the polyhedral gravity field rather than the point-mass approximation, which shifts the solution because a real asteroid is not spherically symmetric.
Worked example. A 1.88 km body at
$\beta = 2.4 \times 10^{-3}$ and$R = 1$ AU yields an AEP 5.64 km anti-sunward, confirming the SRP-dominated regime.
The closed-form solution places the AEP directly anti-sunward — which is inside the asteroid's shadow. A spacecraft there receives no sunlight, meaning no power and, self-defeatingly, no SRP, so the equilibrium it was placed at does not exist.
Illumination state on the equatorial plane, from vectorized Möller–Trumbore ray casting. The collinear AEP falls squarely inside the umbra.
Real mission designs resolve this by displacing the equilibrium toward the terminator plane, where the spacecraft stays illuminated at grazing incidence. The tradeoff is that the SRP force is no longer antiparallel to gravity, so the balance must be closed by the component along the Sun line while a transverse residual is absorbed by control authority.
PyAsteroid supports both configurations. Which one is appropriate depends on the mission profile, and it changes the control problem — see Known Limitations.
Linearizing about the equilibrium yields a Jacobian with a real positive eigenvalue: the point is a saddle. Any displacement grows exponentially with time constant
For representative kilometer-scale cases
That coincidence is not incidental. In most control problems the plant dynamics and the disturbance occupy separable timescales, so stabilization and disturbance rejection can be designed independently. Here they do not separate: the gravity field varies as the body rotates, on the same timescale over which the instability grows. The controller cannot treat the rotating gravity field as quasi-static, and it cannot treat the instability as fast relative to the disturbance. Both must be handled at once.
Every modeled force for a 1.88 km body. Gravity and SRP cross at the AEP; the tidal, albedo, and thermal terms sit three or more orders of magnitude below.
| Force | Role |
|---|---|
| Asteroid gravity | Primary inward attraction; irregular and time-varying with body rotation |
| Solar radiation pressure | Primary outward force; balances gravity at the AEP |
| Coriolis | Exactly known from velocity and frame rate — feedforward, not disturbance |
| Solar tidal | Negligible in the SRP-dominated regime, retained for completeness |
| Albedo | Reflected sunlight from the body surface; small correction |
| Thermal IR | Re-radiated heat from the body surface; small correction |
git clone https://github.com/JCarroll-OU/PyAsteroid.git
cd PyAsteroid
python -m pip install -r requirements.txtDependencies
| Package | Purpose |
|---|---|
numpy |
Array operations throughout |
scipy |
Integration, root finding, eigenvalue analysis |
polyhedral_gravity |
Exact polyhedral gravity via Tsoulis' method |
matplotlib |
3D visualization and interactive plotting |
from pyasteroid import AsteroidGenerator, GravityField, PGCR3BP, Plotter
# 1. Generate a randomized asteroid
gen = AsteroidGenerator(seed=42)
asteroid = gen.generate()
# 2. Build its gravity field from the mesh
field = GravityField(asteroid, density=2000.0) # kg/m^3
# 3. Set up the dynamics model
system = PGCR3BP(
asteroid=asteroid,
gravity=field,
beta=2.4e-3, # spacecraft lightness number
heliocentric_distance=1.0, # AU
)
# 4. Locate the artificial equilibrium point
aep = system.find_aep()
print(f"AEP at {aep} m, |d| = {np.linalg.norm(aep)/1e3:.2f} km")
# 5. Check its stability
eigenvalues = system.linearize(aep)
tau = 1.0 / max(eigenvalues.real)
print(f"Instability time constant: {tau/3600:.2f} hours")
# 6. Visualize
plot = Plotter(asteroid)
plot.add_marker(aep, label="AEP")
plot.show()Represents a single small body: its shape, its mesh, and its rotation state.
The mesh is built by subdividing an icosphere, which gives near-uniform triangle areas across the surface — important because the polyhedral gravity computation sums per-face contributions, and wildly uneven faces degrade accuracy in ways that are hard to diagnose. The base icosphere is then scaled into a tri-axial ellipsoid and perturbed with surface features.
The same body at two subdivision levels. Face count grows as $20 \cdot 4^n$, so cost rises steeply with resolution.
Parameters
| Parameter | Type | Description |
|---|---|---|
semi_axes |
(a, b, c) |
Tri-axial ellipsoid semi-axis lengths in meters. Setting all three equal gives a sphere, useful for validation. |
subdivisions |
int |
Icosphere subdivision level. Face count grows as |
n_features |
int |
Number of surface features (craters and hills) to place. |
feature_params |
— | Controls the Gaussian shape of each feature: amplitude (positive for hills, negative for craters), width, and surface position. |
spin_period |
float |
Rotation period in seconds. |
spin_axis |
(x, y, z) |
Rotation axis direction; need not align with a principal axis. |
Surface features are Gaussian-parameterized: each is a radial displacement of the mesh vertices falling off as
Rotation uses quaternions rather than Euler angles, avoiding gimbal lock and drift over long propagations where a rotation matrix accumulates orthogonality error.
Key methods
| Method | Description |
|---|---|
rotate(t) |
Advance the body's attitude to time t |
vertices() / faces() |
Return the mesh arrays for gravity computation and rendering |
is_illuminated(point, sun_dir) |
Ray-cast shadow test — see GravityField |
Produces randomized asteroid populations for Monte Carlo campaigns. Every geometric parameter is drawn from a Range, so a whole population is specified declaratively and reproduced from a single seed.
Eight bodies from a single generator configuration, rendered at common scale.
Samples from a truncated Gaussian. The truncation matters: a plain Gaussian on a physical quantity will eventually produce a negative semi-axis or an implausibly extreme aspect ratio, and a single invalid body can silently poison a long campaign.
Range(mean=1000.0, std=300.0, min=400.0, max=2500.0)| Parameter | Description |
|---|---|
mean |
Distribution center |
std |
Standard deviation before truncation |
min, max |
Hard bounds; samples outside are rejected and redrawn |
gen = AsteroidGenerator(
size=Range(mean=1000.0, std=300.0, min=400.0, max=2500.0),
elongation=Range(mean=1.3, std=0.2, min=1.0, max=2.0),
flattening=Range(mean=1.2, std=0.15, min=1.0, max=1.8),
n_features=Range(mean=12, std=4, min=0, max=30),
spin_period=Range(mean=6*3600, std=2*3600, min=2*3600, max=20*3600),
seed=42,
)Shape is specified through ratios rather than three independent axes. Drawing
Computes gravitational acceleration from the asteroid's polyhedral mesh using Tsoulis' method, wrapping the polyhedral_gravity library.
This is an exact evaluation for a constant-density polyhedron — not a series truncation. Spherical harmonics diverge inside the circumscribing sphere (the Brillouin sphere), which is precisely the region of interest for close proximity operations. For a strongly elongated body, the Brillouin sphere may enclose the entire operational volume, making the harmonic expansion unusable exactly where accuracy is needed. The polyhedral formulation has no such restriction and remains valid arbitrarily close to the surface.
Acceleration magnitude on the equatorial plane, with the body outline overlaid and direction shown as unit vectors. The departure from spherical symmetry is what the controller has to reject.
Parameters
| Parameter | Type | Description |
|---|---|---|
asteroid |
Asteroid |
Source of the mesh geometry |
density |
float |
Bulk density in kg/m³, assumed constant. Typical asteroid values run 1200–3000. |
Key methods
| Method | Description |
|---|---|
acceleration(points) |
Gravitational acceleration at one or many field points, in the body-fixed frame |
potential(points) |
Gravitational potential |
gm() |
Total |
SRP requires knowing whether a point is illuminated, which means testing whether the ray from the point toward the Sun intersects the asteroid mesh. PyAsteroid uses vectorized Möller–Trumbore ray–triangle intersection, testing all faces simultaneously with NumPy rather than looping. This matters at scale: a Monte Carlo campaign performs this test at every integration step for every run, and a Python-level loop over tens of thousands of faces would dominate total runtime.
The result of that test is shown in the umbra figure above.
The dynamics model. Assembles the gravity field, SRP, and rotating-frame terms into equations of motion, locates equilibria, and analyzes stability.
Parameters
| Parameter | Type | Description |
|---|---|---|
asteroid |
Asteroid |
The secondary body |
gravity |
GravityField |
Gravity model for that body |
beta |
float |
Spacecraft lightness number. Sets the AEP location — see above. |
heliocentric_distance |
float |
Sun–asteroid separation in AU. Sets the mean motion and the solar flux. |
include_tidal |
bool |
Include the solar tidal term. Negligible in the SRP-dominated regime; useful for verifying that claim. |
include_albedo |
bool |
Include reflected sunlight from the body surface. |
include_thermal |
bool |
Include thermal IR re-radiation from the body surface. |
aep_mode |
str |
"collinear" for the pure anti-sunward point, "terminator" for the illuminated displaced solution. |
Key methods
| Method | Description |
|---|---|
dynamics(t, state) |
State derivative for a 6-element scipy.integrate.solve_ivp
|
find_aep() |
Closed-form initial guess refined by Newton iteration against the full field |
linearize(point) |
Jacobian eigenvalues at a point; the largest real part gives |
srp_acceleration(point, t) |
SRP acceleration including the shadow test |
Newton refinement from the closed-form guess, across bodies from 1.0 to 2.6 km and aspect ratios up to 1.84. Convergence is quadratic and reaches the floating-point floor in three to four iterations — the analytic estimate is a good enough starting point that refinement is cheap.
Interactive 3D visualization for asteroid geometry, gravity fields, trajectories, and equilibrium points.
Plotter is deliberately import-free of the domain classes — it never imports Asteroid, GravityField, or PGCR3BP. Interaction is wired through a factory-pattern callback contract: the caller supplies functions producing the data to display, and Plotter invokes them without knowing what produced them.
The payoff is that visualization can be developed, tested, and modified without touching physics code, and the physics classes carry no rendering dependencies. It also means Plotter can display anything conforming to the contract — including data the framework does not currently produce.
The interactive window, with toggles wired through the callback contract.
Key methods
| Method | Description |
|---|---|
add_marker(point, label) |
Place a labeled point, e.g. the AEP |
add_trajectory(states) |
Draw a propagated path |
add_toggle(label, callback) |
Register an interactive button against a data-producing callback |
show() |
Open the interactive window |
The framework's purpose is generating many randomized environments and evaluating a controller across all of them.
gen = AsteroidGenerator(seed=42, **population_config)
results = []
for i in range(N_RUNS):
asteroid = gen.generate()
field = GravityField(asteroid, density=2000.0)
system = PGCR3BP(asteroid=asteroid, gravity=field, beta=beta)
aep = system.find_aep()
results.append(run_controller(system, aep, controller))Because every asteroid derives from a single seeded generator, an entire campaign is reproducible from that seed and the population configuration — and an individual anomalous run can be regenerated in isolation for debugging without re-running the campaign.
SRP disturbance. Beyond geometric randomization, the framework injects heavy noise into the solar flux to represent solar flares and coronal mass ejections. This is the physically motivated worst case for this problem: a flux transient does not merely push the spacecraft, it shifts the equilibrium the spacecraft is trying to hold, because the AEP location depends on the SRP magnitude. The controller must reject a disturbance whose effect is to move the setpoint.
Analytic sphere check. With all semi-axes equal and no surface features, the polyhedral gravity field must converge to the closed-form solution for a uniform sphere,
Mean relative error in $|\mathbf{g}|$ on a shell at $r = 3R$, against the analytic sphere. Convergence is second order in mesh edge length — the expected rate for a linear triangulation of a smooth surface, which indicates the wrapper is behaving correctly rather than merely producing plausible numbers.
Tidal term check. Running with include_tidal=True and False and comparing AEP locations verifies the SRP-dominated regime claim numerically rather than relying on the order-of-magnitude argument alone.
Every figure in this README except Plotter_UI.png is produced by generate_figures.py:
python generate_figures.py # all figures -> Images/
python generate_figures.py --list # available figure names
python generate_figures.py force_budget # regenerate just onePlotter_UI.png is a screenshot of the live interactive window and has no scripted equivalent.
-
Constant density. The polyhedral formulation assumes uniform bulk density. Real asteroids, particularly rubble piles, are heterogeneous. This is a modeling choice, not an implementation gap — heterogeneity would require a mass-concentration or layered model.
-
AEP configuration is unresolved. Whether the target should be the displaced-collinear point or a terminator-plane AEP remains an open design question. The collinear point is analytically cleaner but sits inside the umbra; the terminator solution stays illuminated but does not admit the simple closed-form balance. This choice changes the control problem, not just the setpoint.
-
Timescale coupling. With
$\tau$ comparable to the spin period, the usual separation between stabilization and disturbance rejection does not hold. Controllers designed assuming quasi-static gravity should not be expected to transfer. -
Attitude is not modeled. The spacecraft is treated as a point mass with a fixed effective area-to-mass ratio. Real attitude motion modulates
$\beta$ , and therefore modulates the equilibrium location. -
SRP incidence angle. Work on solar-sail sliding mode control reports failure beyond roughly 30° solar incidence, arising from geometric constraints on the achievable force direction rather than from controller tuning. Since SRP direction is constrained by illumination geometry, this bounds the control authority available in the transverse direction — directly relevant to any terminator-displaced configuration.
- Tsoulis, D. — Analytical computation of the gravity field of a polyhedral source with linearly varying density.
- Möller, T. and Trumbore, B. — Fast, minimum storage ray/triangle intersection.
- McInnes, C. R. — Solar Sailing: Technology, Dynamics and Mission Applications (artificial equilibrium points and the lightness number).
- Scheeres, D. J. — Orbital Motion in Strongly Perturbed Environments (small-body dynamics).