# orbit

`import "github.com/matjam/bunyip/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.

## Constants

<a id="G"></a>

```go
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

<a id="AppendPredict"></a>

### AppendPredict

```go
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.

<a id="AppendPredictRelative"></a>

### AppendPredictRelative

```go
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.

<a id="CircularVelocity"></a>

### CircularVelocity

```go
func CircularVelocity(mu, r float64) float64
```

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

<a id="Energy"></a>

### Energy

```go
func Energy(s State, mu float64) float64
```

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

<a id="EscapeVelocity"></a>

### EscapeVelocity

```go
func EscapeVelocity(mu, r float64) float64
```

EscapeVelocity is the speed needed to leave from radius r.

<a id="Hohmann"></a>

### Hohmann

```go
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:

```go
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
```

<a id="MeanAnomalyFromTrue"></a>

### MeanAnomalyFromTrue

```go
func MeanAnomalyFromTrue(trueAnomaly, e float64) float64
```

MeanAnomalyFromTrue converts a true anomaly on a closed orbit.

<a id="Predict"></a>

### Predict

```go
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.

<a id="PredictRelative"></a>

### PredictRelative

```go
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.

<a id="RK4"></a>

### RK4

```go
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.

<a id="SolveKepler"></a>

### SolveKepler

```go
func SolveKepler(meanAnomaly, e float64) float64
```

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

<a id="SphereOfInfluence"></a>

### SphereOfInfluence

```go
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.

<a id="System"></a>

### System

```go
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.

<a id="TrueAnomalyFromMean"></a>

### TrueAnomalyFromMean

```go
func TrueAnomalyFromMean(meanAnomaly, e float64) float64
```

TrueAnomalyFromMean converts a mean anomaly on a closed orbit.

<a id="VisViva"></a>

### VisViva

```go
func VisViva(mu, r, a float64) float64
```

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

## Types

<a id="Body"></a>

<a id="Body.Pos"></a>

<a id="Body.Vel"></a>

<a id="Body.Mass"></a>

### Body

```go
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.

<a id="Elements"></a>

<a id="Elements.SemiMajorAxis"></a>

<a id="Elements.Eccentricity"></a>

<a id="Elements.Inclination"></a>

<a id="Elements.Node"></a>

<a id="Elements.ArgPeriapsis"></a>

<a id="Elements.TrueAnomaly"></a>

### Elements

```go
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.

<a id="Around"></a>

#### Around

```go
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.

<a id="Circular"></a>

#### Circular

```go
func Circular(radius float64) Elements
```

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

<a id="ElementsOf"></a>

#### ElementsOf

```go
func ElementsOf(s State, mu float64) Elements
```

ElementsOf converts a state to classical elements for mu.

<a id="Elements.Apoapsis"></a>

#### Elements.Apoapsis

```go
func (e Elements) Apoapsis() float64
```

Apoapsis is the farthest distance; infinite for open orbits.

<a id="Elements.AtTime"></a>

#### Elements.AtTime

```go
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.

<a id="Elements.MeanMotion"></a>

#### Elements.MeanMotion

```go
func (e Elements) MeanMotion(mu float64) float64
```

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

<a id="Elements.Periapsis"></a>

#### Elements.Periapsis

```go
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.

<a id="Elements.Period"></a>

#### Elements.Period

```go
func (e Elements) Period(mu float64) float64
```

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

<a id="Elements.State"></a>

#### Elements.State

```go
func (e Elements) State(mu float64) State
```

State converts elements to a position and velocity for mu.

<a id="Kepler"></a>

<a id="Kepler.Primary"></a>

<a id="Kepler.Elements"></a>

<a id="Kepler.Mu"></a>

### Kepler

```go
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.

<a id="Kepler.MarshalJSON"></a>

#### Kepler.MarshalJSON

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

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

<a id="Kepler.UnmarshalJSON"></a>

#### Kepler.UnmarshalJSON

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

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

<a id="Settings"></a>

<a id="Settings.G"></a>

<a id="Settings.TimeScale"></a>

<a id="Settings.Substeps"></a>

<a id="Settings.Softening"></a>

<a id="Settings.Scale"></a>

<a id="Settings.Origin"></a>

### Settings

```go
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.

<a id="Ship"></a>

### Ship

```go
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.

<a id="Simulation"></a>

<a id="Simulation.Bodies"></a>

<a id="Simulation.G"></a>

<a id="Simulation.Softening"></a>

<a id="Simulation.Time"></a>

### Simulation

```go
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:

```go
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
```

<a id="Simulation.Accelerations"></a>

#### Simulation.Accelerations

```go
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.

<a id="Simulation.Barycenter"></a>

#### Simulation.Barycenter

```go
func (s *Simulation) Barycenter() Vec3
```

Barycenter is the centre of mass.

<a id="Simulation.Energy"></a>

#### Simulation.Energy

```go
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.

<a id="Simulation.FieldAt"></a>

#### Simulation.FieldAt

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

FieldAt is the gravitational acceleration the bodies produce at p.

<a id="Simulation.Step"></a>

#### Simulation.Step

```go
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.

<a id="State"></a>

<a id="State.Pos"></a>

<a id="State.Vel"></a>

### State

```go
type State struct{ Pos, Vel Vec3 }
```

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

<a id="Propagate"></a>

#### Propagate

```go
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:

```go
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
```

<a id="Thrust"></a>

<a id="Thrust.Accel"></a>

### Thrust

```go
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.

<a id="Vec3"></a>

<a id="Vec3.X"></a>

<a id="Vec3.Y"></a>

<a id="Vec3.Z"></a>

### Vec3

```go
type Vec3 struct{ X, Y, Z float64 }
```

Vec3 is a double-precision vector.

<a id="FromLin"></a>

#### FromLin

```go
func FromLin(v lin.Vec3) Vec3
```

FromLin converts an engine vector.

<a id="V3"></a>

#### V3

```go
func V3(x, y, z float64) Vec3
```

V3 makes a Vec3.

<a id="Vec3.Add"></a>

#### Vec3.Add

```go
func (a Vec3) Add(b Vec3) Vec3
```

Add returns a + b.

<a id="Vec3.Cross"></a>

#### Vec3.Cross

```go
func (a Vec3) Cross(b Vec3) Vec3
```

Cross is the cross product.

<a id="Vec3.Dot"></a>

#### Vec3.Dot

```go
func (a Vec3) Dot(b Vec3) float64
```

Dot is the dot product.

<a id="Vec3.Len"></a>

#### Vec3.Len

```go
func (a Vec3) Len() float64
```

Len is the length.

<a id="Vec3.Lin"></a>

#### Vec3.Lin

```go
func (a Vec3) Lin() lin.Vec3
```

Lin converts to an engine vector, losing precision.

<a id="Vec3.Mul"></a>

#### Vec3.Mul

```go
func (a Vec3) Mul(s float64) Vec3
```

Mul scales the vector.

<a id="Vec3.Norm"></a>

#### Vec3.Norm

```go
func (a Vec3) Norm() Vec3
```

Norm returns the unit vector; zero stays zero.

<a id="Vec3.Sub"></a>

#### Vec3.Sub

```go
func (a Vec3) Sub(b Vec3) Vec3
```

Sub returns a - b.