Bunyip a game engine in Go
GitHub

Package github.com/matjam/bunyip/orbit

orbit

Package orbit computes celestial mechanics for space games. It provides orbital elements and state vectors, two-body propagation that handles circles, ellipses, parabolas and hyperbolas alike, an N-body integrator for planets and moons pulling on each other, and helpers for the numbers a player sees (periods, apsides, escape velocity, transfer burns). Everything is double precision, because a planetary system does not fit in float32. Nothing assumes our own solar system. Bodies are masses at positions, the gravitational constant is yours to set, and units are whatever you choose (metres and seconds with the real G, or a game's own units with G = 1).

System connects the package to the entity component system. Body components hold positions and velocities at astronomical scale. The system moves them, analytically for entities on a Kepler orbit and numerically for free bodies under thrust and gravity, and writes scaled positions into gfx.Transform relative to a floating origin, so rendering keeps its precision at any camera position. Kepler components preserve their orbital phase through ECS saves and prefab JSON. Register the component and resource types with ecs.Register before saving; omit the orbit system's private query cache with ecs.SaveOptions.SkipUnregistered.

Index

Constants

const G = 6.67430e-11 // m³ kg⁻¹ s⁻²

G is the gravitational constant in SI units. Nothing here depends on SI: set Settings.G or Simulation.G to 1 and choose your own units of mass, distance and time, and every formula still holds. Real-world bodies live in the sol subpackage.

Functions

AppendPredict source

func AppendPredict(dst []lin.Vec3, w *ecs.World, ship ecs.Entity, seconds float64, samples int) []lin.Vec3

AppendPredict is Predict appending the path to dst, for a caller that predicts every frame and keeps one slice for it: pass the last frame's path sliced to zero length and nothing is allocated once it has grown. When Predict would return nil, it returns dst unchanged.

AppendPredictRelative source

func AppendPredictRelative(dst []lin.Vec3, w *ecs.World, ship, ref ecs.Entity, seconds float64, samples int) []lin.Vec3

AppendPredictRelative is PredictRelative appending the path to dst, for a caller that predicts every frame and keeps one slice for it: pass the last frame's path sliced to zero length and nothing is allocated once it has grown. When PredictRelative would return nil, it returns dst unchanged.

CircularVelocity source

func CircularVelocity(mu, r float64) float64

CircularVelocity is the speed of a circular orbit at radius r.

Energy source

func Energy(s State, mu float64) float64

Energy is the specific orbital energy; negative for bound orbits.

EscapeVelocity source

func EscapeVelocity(mu, r float64) float64

EscapeVelocity is the speed needed to leave from radius r.

Hohmann source

func Hohmann(mu, r1, r2 float64) (dv1, dv2, coast float64)

Hohmann plans the two-burn transfer between circular orbits of radii r1 and r2: the nonnegative delta-v magnitude of each burn and the coast time between them. Both radii and mu must be positive and use consistent units.

Example 
package main

import (
	"fmt"

	"github.com/matjam/bunyip/orbit"
	"github.com/matjam/bunyip/orbit/sol"
)

func main() {
	// Low Earth orbit to geostationary: two burns and a coast.
	dv1, dv2, coast := orbit.Hohmann(sol.MuEarth, 6678e3, 42164e3)
	fmt.Printf("burn %.2f km/s, coast %.1f h, burn %.2f km/s\n", dv1/1000, coast/3600, dv2/1000)
}
Output
burn 2.43 km/s, coast 5.3 h, burn 1.47 km/s

MeanAnomalyFromTrue source

func MeanAnomalyFromTrue(trueAnomaly, e float64) float64

MeanAnomalyFromTrue converts a true anomaly on a closed orbit.

Predict source

func Predict(w *ecs.World, ship ecs.Entity, seconds float64, samples int) []lin.Vec3

Predict integrates a copy of a ship forward and returns its path in world units (scaled, relative to the origin), for drawing where it is heading. seconds is a nonnegative duration in simulation time units, independent of Settings.TimeScale. The result contains exactly samples points, excluding the starting point and including the end of the duration. Kepler bodies move during prediction; other Ships are excluded from gravity, and massive bodies with neither Ship nor Kepler stay fixed. Settings and the ship's Body must exist and samples must be positive, otherwise it returns nil. Current Thrust is held constant throughout. To reuse a slice from frame to frame, call AppendPredict.

PredictRelative source

func PredictRelative(w *ecs.World, ship, ref ecs.Entity, seconds float64, samples int) []lin.Vec3

PredictRelative is Predict in the moving frame of another body: each point is shifted by how far that body will have moved by then, so the path of a ship orbiting a planet draws as a loop around the planet where it is now rather than a streak across the sky. With ref None it is the inertial path. Only a reference with Kepler is advanced; a free reference body supplies no moving-frame correction. Requirements, units and sample placement are the same as Predict. To reuse a slice from frame to frame, call AppendPredictRelative.

RK4 source

func RK4(pos, vel Vec3, dt float64, accel func(pos, vel Vec3) Vec3) (Vec3, Vec3)

RK4 advances a particle by dt in an acceleration field, a fourth-order step suited to spacecraft under thrust plus gravity.

SolveKepler source

func SolveKepler(meanAnomaly, e float64) float64

SolveKepler returns the eccentric anomaly for a mean anomaly and eccentricity below one.

SphereOfInfluence source

func SphereOfInfluence(a, m, M float64) float64

SphereOfInfluence returns the Laplace sphere-of-influence approximation a*(m/M)^0.4 for a body's orbital semi-major axis a and masses m and M. It is a patched-conic boundary estimate, not an exact gravity boundary.

System source

func System(w *ecs.World, dt float64)

System advances orbits by dt real seconds (times Settings.TimeScale) and writes scaled positions into existing gfx.Transform components. It creates default Settings when absent. Use a nonnegative scaled timestep: the ship integrator does not integrate backwards. Kepler primary chains should be acyclic and no more than sixteen links deep.

TrueAnomalyFromMean source

func TrueAnomalyFromMean(meanAnomaly, e float64) float64

TrueAnomalyFromMean converts a mean anomaly on a closed orbit.

VisViva source

func VisViva(mu, r, a float64) float64

VisViva is the speed at radius r on an orbit of semi-major axis a.

Types

type Body source

type Body struct {
	Pos, Vel Vec3
	Mass     float64
}

Body is a mass with a position and velocity, used both as an ECS component and inside a Simulation. In a Simulation, Mass zero makes a test particle that feels gravity without exerting it. In the ECS, Ship or Kepler selects motion; a Body with neither stays fixed.

type Elements source

type Elements struct {
	SemiMajorAxis float64 // negative for hyperbolic orbits; periapsis distance when eccentricity is exactly 1
	Eccentricity  float64
	Inclination   float64
	Node          float64 // longitude of the ascending node
	ArgPeriapsis  float64 // argument of periapsis
	TrueAnomaly   float64
}

Elements are classical orbital elements; angles are radians and the reference plane is XY with +Z as its normal. Distances and mu must use consistent units; mu is G times primary mass, with units distance³/time². A zero Elements has no valid orbital radius.

Around source

func Around(w *ecs.World, ship ecs.Entity) (primary ecs.Entity, el Elements, mu float64, ok bool)

Around describes a ship's orbit relative to the massive body whose gravity dominates it, in classical elements, for a readout.

Circular source

func Circular(radius float64) Elements

Circular makes a circular orbit of the given radius in the reference plane.

ElementsOf source

func ElementsOf(s State, mu float64) Elements

ElementsOf converts a state to classical elements for mu.

Apoapsis source

func (e Elements) Apoapsis() float64

Apoapsis is the farthest distance; infinite for open orbits.

AtTime source

func (e Elements) AtTime(mu, dt float64) Elements

AtTime returns the elements after dt simulation time units. Closed orbits advance mean anomaly; open orbits use Propagate and ElementsOf.

MeanMotion source

func (e Elements) MeanMotion(mu float64) float64

MeanMotion is the average angular rate in radians per simulation time unit.

Periapsis source

func (e Elements) Periapsis() float64

Periapsis returns SemiMajorAxis*(1-Eccentricity), the closest distance for elliptic and hyperbolic orbits. For a parabola, read SemiMajorAxis directly: this formula returns zero at eccentricity 1.

Period source

func (e Elements) Period(mu float64) float64

Period is the orbital period for mu; infinite for open orbits.

State source

func (e Elements) State(mu float64) State

State converts elements to a position and velocity for mu.

type Kepler source

type Kepler struct {
	Primary  ecs.Entity
	Elements Elements
	Mu       float64
	// contains filtered or unexported fields
}

Kepler puts an entity on an analytic two-body orbit around Primary: planets around a star, moons around a planet. Its Body is written each update from the elements; Mu zero means G times the primary's mass. JSON encoding preserves the elapsed orbital phase for saves and prefabs; older JSON without elapsed time starts at the elements' epoch.

MarshalJSON source

func (k Kepler) MarshalJSON() ([]byte, error)

MarshalJSON preserves the orbital phase in ECS saves and prefab JSON.

UnmarshalJSON source

func (k *Kepler) UnmarshalJSON(data []byte) error

UnmarshalJSON restores the orbital phase. Older JSON without Elapsed starts at the epoch specified by Elements.

type Settings source

type Settings struct {
	G         float64 // zero means the real constant
	TimeScale float64 // simulation time units per real second; zero means 1
	Substeps  int     // minimum integration steps per update for ships; nonpositive means 8
	Softening float64 // distance softening: squared and added to squared separations; zero disables it
	// Rendering: scene units per simulation distance unit, and the floating origin that is
	// subtracted from every position before scaling so the scene stays
	// near zero where float32 is precise. Move Origin with the camera.
	Scale  float64 // zero means 1
	Origin Vec3
}

Settings is the world resource for the orbit system.

type Ship source

type Ship struct{}

Ship is a marker for free bodies: they are integrated numerically under massive non-Ship bodies' gravity plus their own Thrust. Ships do not exert gravity on one another, regardless of Mass. Without it, a Body with no Kepler stays put (a star at the origin). Use either Ship or Kepler on an entity, not both.

type Simulation source

type Simulation struct {
	Bodies    []Body
	G         float64 // zero means the real constant
	Softening float64 // simulation distance: its square is added to squared separations
	Time      float64 // elapsed simulation time, advanced by Step
	// contains filtered or unexported fields
}

Simulation integrates bodies under their mutual gravity with the leapfrog (kick-drift-kick) scheme. Use a sufficiently small fixed step to limit integration error; energy is not conserved exactly. The zero value is an empty simulation with the SI gravitational constant.

Example 
package main

import (
	"fmt"
	"math"

	"github.com/matjam/bunyip/orbit"
)

func main() {
	// Two stars of equal mass circling each other, in game units.
	sim := orbit.Simulation{G: 1, Bodies: []orbit.Body{
		{Pos: orbit.V3(-5, 0, 0), Vel: orbit.V3(0, -5, 0), Mass: 500},
		{Pos: orbit.V3(5, 0, 0), Vel: orbit.V3(0, 5, 0), Mass: 500},
	}}
	e0 := sim.Energy()
	for range 10000 {
		sim.Step(0.001)
	}
	drift := math.Abs((sim.Energy() - e0) / e0)
	fmt.Printf("energy conserved to 1e-9: %v, separation %.2f\n", drift < 1e-9, sim.Bodies[1].Pos.Sub(sim.Bodies[0].Pos).Len())
}
Output
energy conserved to 1e-9: true, separation 10.00

Accelerations source

func (s *Simulation) Accelerations(out []Vec3)

Accelerations fills out with the gravitational acceleration on each body. out must have length at least len(Bodies). Positive Softening avoids the singularity when distinct bodies occupy the same position. From 256 bodies it shares the bodies out across goroutines, one range per processor; each body's sum still runs over the others in order, so the result is the same as on one goroutine. Do not change Bodies while it runs.

Barycenter source

func (s *Simulation) Barycenter() Vec3

Barycenter is the centre of mass.

Energy source

func (s *Simulation) Energy() float64

Energy returns total kinetic plus unsoftened Newtonian potential energy. It ignores Softening and omits potential terms for coincident bodies, so it is an integration diagnostic for separated, unsoftened systems.

FieldAt source

func (s *Simulation) FieldAt(p Vec3) Vec3

FieldAt is the gravitational acceleration the bodies produce at p.

Step source

func (s *Simulation) Step(dt float64)

Step advances every body and Time by dt simulation time units. The accelerations it ends with are the ones the next Step starts with, so it computes them once per step, not twice. Changing a body's position or mass, adding or removing bodies, or changing G or Softening between steps is noticed, and the next Step computes them afresh.

type State source

type State struct{ Pos, Vel Vec3 }

State is a position and velocity relative to the body being orbited.

Propagate source

func Propagate(s State, mu, dt float64) State

Propagate advances a state by dt simulation time units under two-body gravity using the universal variable formulation for any conic. It solves iteratively in float64 and does not report failure to converge. mu must be positive and the initial distance from the primary nonzero. A zero dt returns s unchanged.

Example 
package main

import (
	"fmt"

	"github.com/matjam/bunyip/orbit"
)

func main() {
	// A fictional world in game units: G is 1 and the star has mass 1000.
	const mu = 1000
	el := orbit.Elements{SemiMajorAxis: 20, Eccentricity: 0.3}
	s := el.State(mu)
	half := orbit.Propagate(s, mu, el.Period(mu)/2)
	fmt.Printf("periapsis %.0f, apoapsis %.0f, at half a period r = %.1f\n", el.Periapsis(), el.Apoapsis(), half.Pos.Len())
}
Output
periapsis 14, apoapsis 26, at half a period r = 26.0

type Thrust source

type Thrust struct{ Accel Vec3 }

Thrust is the constant acceleration applied to a Ship, in simulation distance units per simulation time unit squared, in world axes.

type Vec3 source

type Vec3 struct{ X, Y, Z float64 }

Vec3 is a double-precision vector.

FromLin source

func FromLin(v lin.Vec3) Vec3

FromLin converts an engine vector.

V3 source

func V3(x, y, z float64) Vec3

V3 makes a Vec3.

Add source

func (a Vec3) Add(b Vec3) Vec3

Add returns a + b.

Cross source

func (a Vec3) Cross(b Vec3) Vec3

Cross is the cross product.

Dot source

func (a Vec3) Dot(b Vec3) float64

Dot is the dot product.

Len source

func (a Vec3) Len() float64

Len is the length.

Lin source

func (a Vec3) Lin() lin.Vec3

Lin converts to an engine vector, losing precision.

Mul source

func (a Vec3) Mul(s float64) Vec3

Mul scales the vector.

Norm source

func (a Vec3) Norm() Vec3

Norm returns the unit vector; zero stays zero.

Sub source

func (a Vec3) Sub(b Vec3) Vec3

Sub returns a - b.

Source files

example_test.go lockstep_test.go nbody.go nbody_reuse_test.go orbit.go orbit_test.go perf_bench_test.go persistence.go persistence_test.go predict_test.go system.go