physics

Documentation

Overview

Package physics implements the spacecraft propagation layer: state vectors, integrators (Verlet + RK4), and patched-conic SOI handling.

Index

• Constants
• func Accel(r orbital.Vec3, mu float64) orbital.Vec3
• func AirRelativeVelocity(r, v orbital.Vec3, primary bodies.CelestialBody) orbital.Vec3
• func AtmosphereOmega(primary bodies.CelestialBody) orbital.Vec3
• func AtmosphericDensity(primary bodies.CelestialBody, altitude float64) float64
• func DragAccel(r, v orbital.Vec3, primary bodies.CelestialBody, bc float64) orbital.Vec3
• func DynamicPressure(r, v orbital.Vec3, primary bodies.CelestialBody) float64
• func GravityOnly(mu float64) func(r, v orbital.Vec3, t float64) orbital.Vec3
• func SOIRadius(body, primary bodies.CelestialBody) float64
• func SemimajorAxis(s StateVector, mu float64) float64
• func SpecificEnergy(s StateVector, mu float64) float64
• func SurfaceGravity(b bodies.CelestialBody) float64
• type Primary
  • func FindPrimary(system bodies.System, rInertial orbital.Vec3, ...) Primary
• type StateVector
  • func ClampToSurface(s StateVector, primary bodies.CelestialBody) (StateVector, bool)
  • func KeplerStep(s StateVector, mu, dt float64) (StateVector, bool)
  • func Rebase(s StateVector, oldPrimaryInertial, newPrimaryInertial orbital.Vec3, ...) StateVector
  • func StepRCSPulse(s StateVector, dirUnit orbital.Vec3, dvQuantum, mAfter float64) StateVector
  • func StepRK4(s StateVector, dt float64, ...) StateVector
  • func StepVerlet(s StateVector, mu, dt float64) StateVector
  • func StepVerletWithAccel(s StateVector, mu, dt float64, extraAccel func(r, v orbital.Vec3) orbital.Vec3) StateVector

Constants

const MinSpinPeriodSeconds = 1.0

MinSpinPeriodSeconds is the floor on a body's rotation period when computing its spin vector. No physical body rotates faster than once per second; clamping here prevents corrupt catalog data from producing an unbounded ω that would corrupt the drag integrator.

Functions

func Accel

func Accel(r orbital.Vec3, mu float64) orbital.Vec3

Accel returns the two-body gravitational acceleration on the spacecraft due to the primary with gravitational parameter mu. Direction is toward the primary (origin in the relative frame).

func AirRelativeVelocity

func AirRelativeVelocity(r, v orbital.Vec3, primary bodies.CelestialBody) orbital.Vec3

AirRelativeVelocity returns the craft's velocity relative to the co-rotating atmosphere: v_rel = v − ω×r, where ω is the body's spin vector (AtmosphereOmega). Single source of truth for "air-relative velocity" — the quantity the drag model opposes, the parachute auto-deploy gate measures, the surface-arrival chute route tests, and the descent HUD reads. Sharing it keeps those consumers in agreement by construction instead of re-deriving v_rel from independent ω sources that could drift apart. v0.12 Slice 3 (ADR 0008).

func AtmosphereOmega

func AtmosphereOmega(primary bodies.CelestialBody) orbital.Vec3

AtmosphereOmega returns the body's spin angular-velocity vector (rad/s) in the world inertial frame. Tilted along the body's physical spin axis per AxialTilt + AxialAzimuth, matching render.BodyRotationAxisWorld so the drag frame agrees with the integrator's surface-co-rotation frame and the renderer's body- fixed projection (ADR 0003). Returns zero when the body has no rotation period set.

Period selection mirrors render.rotationPeriodSeconds: tidally- locked bodies use SideralOrbit, free bodies use SideralRotation (both stored as hours / days respectively per catalog convention).

Pre-v0.11.2 this was a Z-aligned approximation ("the atmosphere co-rotates about world +Z"). The Z-aligned shortcut diverged from the integrator's tilted-axis surface velocity for any body with non-zero AxialTilt; ADR 0003 unifies the convention.

func AtmosphericDensity

func AtmosphericDensity(primary bodies.CelestialBody, altitude float64) float64

AtmosphericDensity returns ρ at altitude h above primary's surface using the body's exponential-density model: ρ(h) = ρ₀ · exp(-h / H) for h < cutoff, 0 above. Returns 0 when the primary has no atmosphere.

func DragAccel

func DragAccel(r, v orbital.Vec3, primary bodies.CelestialBody, bc float64) orbital.Vec3

DragAccel returns the atmospheric-drag acceleration vector on a craft with state (r, v) relative to primary, given a ballistic coefficient BC = C_D · A / m (m²/kg). Returns zero outside the body's atmosphere (no Atmosphere defined, altitude above cutoff, or BC ≤ 0).

Drag direction opposes the craft's velocity *relative to the atmosphere*: v_rel = v - ω × r where ω is the body's spin vector. That distinction matters for low-altitude flight — Earth's surface rotates ~465 m/s eastward at the equator, so a craft sitting at LEO orbital speed (~7800 m/s) feels drag against ~7335 m/s of relative flow, not the full inertial speed.

Magnitude: a = 0.5 · ρ · |v_rel|² · BC, applied along -v̂_rel. v0.8.4+.

func DynamicPressure

func DynamicPressure(r, v orbital.Vec3, primary bodies.CelestialBody) float64

DynamicPressure returns q = 0.5 · ρ · |v_rel|² (Pa) for a craft at state (r, v) relative to primary, using the same air-relative velocity (v_rel = v − ω × r) and exponential-density model as DragAccel. Returns 0 outside the body's atmosphere (no Atmosphere, altitude above cutoff, below the surface). v0.12 Slice 3 (ADR 0008): the parachute auto-deploy gate is expressed in q, which is body-agnostic — no body-specific deploy-altitude constant.

func GravityOnly

func GravityOnly(mu float64) func(r, v orbital.Vec3, t float64) orbital.Vec3

GravityOnly wraps Accel as an RK4-compatible accelFn closure.

func SOIRadius

func SOIRadius(body, primary bodies.CelestialBody) float64

SOIRadius returns the sphere-of-influence radius of `body` orbiting `primary`, in meters. SOI = a_body * (m_body / m_primary)^(2/5). Returns 0 for bodies without orbital data (e.g. the system primary itself).

func SemimajorAxis

func SemimajorAxis(s StateVector, mu float64) float64

SemimajorAxis returns the orbit's semimajor axis a = -μ/(2ε) for bound orbits; NaN for parabolic (ε≈0) and negative for hyperbolic trajectories.

func SpecificEnergy

func SpecificEnergy(s StateVector, mu float64) float64

SpecificEnergy returns the two-body orbital specific energy ε = v²/2 − μ/r (J/kg). Used as an invariant by energy-conservation tests.

func SurfaceGravity

func SurfaceGravity(b bodies.CelestialBody) float64

SurfaceGravity returns the magnitude of gravitational acceleration at the body's mean radius (m/s²). Diagnostic helper for tests.

Types

type Primary

type Primary struct {
	Body     bodies.CelestialBody
	Inertial orbital.Vec3
}

Primary captures which body the spacecraft is currently bound to and its inertial position at the relevant sim-time. Physics operates relative to the primary; world code looks this up every few ticks.

func FindPrimary

func FindPrimary(
	system bodies.System,
	rInertial orbital.Vec3,
	positions map[string]orbital.Vec3,
) Primary

FindPrimary picks the best primary for a spacecraft at inertial position rInertial. Rule: smallest SOI that contains the spacecraft wins; default to the system primary if none contain it.

Each body's SOI is computed against its actual parent (v0.5.0+): planets against the system primary, moons against their planet (per CelestialBody.ParentID). Pre-v0.5.0 every body's SOI was computed against the system primary, which silently treated moons as planets — a Luna at 384k km from Earth had its SOI sized as if it orbited the Sun at 384k km, an absurd value. Fixing this is what enables the Earth-and-Luna nested-SOI walk.

type StateVector

type StateVector struct {
	R orbital.Vec3 // position relative to primary
	V orbital.Vec3 // velocity relative to primary
	M float64      // current total mass (wet)
}

StateVector is the spacecraft's kinematic state relative to the current primary body (see patched-conic model in plan §Phase 1). All SI: position in m, velocity in m/s, mass in kg.

func ClampToSurface

func ClampToSurface(s StateVector, primary bodies.CelestialBody) (StateVector, bool)

ClampToSurface tests whether the state has crossed below the primary's mean radius. On hit it projects the position back to the surface along r̂ and zeros velocity — the craft is treated as "landed" / sitting at altitude 0. Returns (state, false) when the craft is above the surface.

v0.8.5: minimum-viable surface-impact handling. Without this, a craft aerobraking past altitude 0 keeps falling toward r=0; gravity (1/r²) blows up and the symplectic integrator slingshots the craft back out at huge velocity. A real "crashed" status with destruction / structural damage is deferred to v0.9+ — for now landed and zero-velocity is a valid resting state the HUD can read.

Degenerate r=0 input picks +X for the surface direction so callers never see NaN, though in practice the integrator never reaches r=0 because every sub-step is followed by this clamp.

func KeplerStep

func KeplerStep(s StateVector, mu, dt float64) (StateVector, bool)

KeplerStep advances a state vector by dt using analytic Kepler propagation. Snapshots orbital elements, advances mean anomaly by n·dt, and reconstructs (r, v) from the new true anomaly. Exact (modulo Newton iteration tolerance in SolveKepler) — no numerical drift, regardless of dt size relative to orbital period.

Returns ok=false for hyperbolic / parabolic orbits (e ≥ 1) and degenerate inputs; the caller falls back to numerical integration (Verlet). Used for the v0.4.3 "warp lock": when warp > 1× and no active burn, integrateSpacecraft takes one KeplerStep per tick instead of looping Verlet sub-steps that drift eccentricity over many orbits.

SOI handling is the caller's responsibility: KeplerStep does not detect or rebase across primary boundaries. integrateSpacecraft guards against SOI crossings by checking apo < primary's SOI before taking the analytic path.

func Rebase

func Rebase(
	s StateVector,
	oldPrimaryInertial, newPrimaryInertial orbital.Vec3,
	dvInertial orbital.Vec3,
) StateVector

Rebase converts a state vector expressed relative to oldPrimary into one expressed relative to newPrimary. Both primaries' inertial positions are required. Velocity transforms by vector subtraction only if the primaries have non-zero relative velocity; in v0.1 planets are on fixed Keplerian tracks with velocities computed on demand, so we accept a relative-velocity parameter (dvInertial = v_oldPrimary - v_newPrimary).

func StepRCSPulse

func StepRCSPulse(s StateVector, dirUnit orbital.Vec3, dvQuantum, mAfter float64) StateVector

StepRCSPulse applies one RCS pulse to a state vector: a Δv of magnitude dvQuantum along the given unit direction, with mass updated to mAfter. No sub-stepping needed at the cm/s scale a pulse delivers — the velocity change is instantaneous from the integrator's perspective.

Position is unchanged by the pulse itself (impulsive approximation); position evolves under the next gravity step. v0.8.0+.

dirUnit must be a unit vector (caller is responsible — typically spacecraft.DirectionUnit). Caller passes mAfter = total mass after the monoprop debit so the returned state's M reflects post-pulse mass without StepRCSPulse needing to know about Isp / propellant.

func StepRK4

func StepRK4(
	s StateVector,
	dt float64,
	accelFn func(r, v orbital.Vec3, t float64) orbital.Vec3,
	t float64,
) StateVector

StepRK4 advances a StateVector by dt using classical 4th-order Runge-Kutta. Not symplectic — energy drifts over long horizons — but handles non- conservative forces (thrust during a finite burn, drag) cleanly, which Verlet does not. Parked for Phase 2 burns per plan §Phase 1 rationale; until then, Verlet is the default integrator.

accelFn(r, v, t) lets callers inject an applied thrust on top of two-body gravity. For pure gravity, pass a closure that computes Accel(r, μ).

func StepVerlet

func StepVerlet(s StateVector, mu, dt float64) StateVector

StepVerlet advances a StateVector by dt using the symplectic velocity- Verlet integrator. Single force evaluation per step (cheap under warp), conserves specific orbital energy over long horizons — see plan §Phase 1 rationale. Returns the new state; input is untouched.

The algorithm (for a = a(r) only, true for pure two-body gravity):

r' = r + v·dt + ½·a·dt²
a' = Accel(r')
v' = v + ½·(a + a')·dt

func StepVerletWithAccel

func StepVerletWithAccel(s StateVector, mu, dt float64, extraAccel func(r, v orbital.Vec3) orbital.Vec3) StateVector

StepVerletWithAccel is StepVerlet with an extra non-conservative acceleration term a_extra(r, v) added to the gravitational accel at both half-kicks. Used by v0.8.4+ to fold atmospheric drag into the integrator without splitting it across separate Verlet + drag steps (which would lose the symplectic property's stability margin).

Velocity-dependent forces aren't strictly symplectic in this naive fold-in (the canonical fix is operator splitting), but for drag magnitudes typical of LEO and below the energy drift over a single tick is negligible compared with the actual deorbit physics we want to capture.

extraAccel = nil collapses to plain two-body Verlet — same numerical path as the original StepVerlet.