Just an essay regarding celestial mechanics in 2D

I spent a couple of days investigating some new approaches regarding celestial mechanics in 2D (i.e. orbits and planetary motion), and I wanted to share some observations and techniques I've developed.

First of all, it is relatively easy to implement a complex system of bodies, including moons, moons of moons, and stellar systems where two stars orbit around a common center of mass, even for a junior. However, if you want a stable simulation (and you do unless you want chaos), each gravity pair should be treated in isolation, so that the whole system is set up as a group of independent two-body interactions.

Second, I'm obviously talking about simplified models, not true astrophysics. Things can be intentionally exaggerated, scales can be non-realistic, distances shrunk etc.

To implement such physics, semi-implicit Euler method seems to be the easiest and a very robust approach. It's superior to Verlet imo, and can be implemented by rookies. It is, however, sensitive to sudden changes in framerate and time steps should be done in FixedUpdate for improved numerical stability.

Semi-implicit Euler differs from the usual Euler integration such that instead of doing (pseudo-code)

Code (csharp):

1. velocity[t+1] = velocity[t] + acceleration * deltaTime;
2. position[t+1] = position[t] + velocity[t] * deltaTime;

you do

Code (csharp):

1. velocity[t+1] = velocity[t] + acceleration * deltaTime;
2. position[t+1] = position[t] + velocity[t+1] * deltaTime;

It's a trivial modification but has the effect of diminishing translational errors, while conserving the total energy of the system, thus it provides physical consistency throughout.

Such a simulation can be neatly time-scaled on the fly (the errors are negligible unless you scale it to values greater than multiple frames, because this will inevitable scale up the vectors to large magnitudes), and velocities, positions and relationships can be neatly configured and even predicted during edit mode. The way I've implemented it has a serialized (and fixed aka "non-reorderable") list of properties for each body, where I can set the mass, the color, etc. and when it comes to picking the parent body, you just set its index. But before the simulation starts, I build a tree out of this data so that it recursively updates from the top down. I.e. Sun -> planet1 -> its moons -> planet2 -> planet 3 -> its moon X.



I simulate each body as if it was relative to the world origin, and once its velocity and position are updated internally, I then position the body relative to where its parent truly lies in space (which is always a correct position thanks to iterating top to bottom). Because each branch is relatively short (1 to 3 steps at most, if we ignore moons of moons which are only theoretically possible according to sciences) I integrate the body's absolute position on the fly, through backward propagation (basically if I need to refresh the moon X's world position, I first find where planet 3 is, but to do that I need to find where the Sun is ... Once I hit "I have no parent" I return everything added up, just like transform.position works).

With all that said, the system is rich enough that you may introduce more than just circular orbits. Elliptical, parabolical, and even hyperbolical orbits (technically all conic sections) all arise naturally, depending on the initial values. It is also easy to come up with a system of naturally arising gravitational resonances (Kepler's 3rd law). Kepler's 1st and 2nd law (accelerating near the Sun/focus when following eccentric orbits) are also obeyed as a by-product of simulating actual gravitational acceleration.



This is the formula which is plugged in as an 'acceleration' just like any other continuous force you might add to the system

Code (csharp):

1. G_constant * parent_mass / distance^3 * 2d_vector_from_body_to_parent

I simply ignore the mass of the pulled object, the masses in my model are strictly one-way, so there is no tidal effect, for performance and simplicity.

G constant is big G in physics, and realistically it's a very small (empirical) value (it's not 9.81, that's small g), but for gaming purposes and simple simulations you may keep it at 1, because you can always tweak the masses and velocities, so it can be totally ignored. Similarly, masses are only relative to each other, you don't have to bother yourself with realistic units. Though, conveniently, Earth's mass is around 6 petagrams (5.927x10^12 kgs), so you might base everything on some coarse convention. It's not pentagrams, goddamnit.

https://www.google.com/search?q=petagram

Finally, if you want to simulate dozens of body in your game, constantly, while also doing other game stuff, it makes sense to think of a system that will simulate the motion up to a certain point (perhaps only during development), and then "bake" that into deterministic paths to simply remove rigid and typically immovable objects (such as planets) from having to compute complex formulas all the time. Maybe do that only for ships and other small bodies subjected to merciless gravity of true celestial objects?

If you do want to "bake" your paths, however, you then need to find a way to create actual parametric orbital representations. For example, a circle has a center and radius, and that's usually more than enough to come up with all its points, it's tangents, all kinds of distances, chords, angles, circumference etc.

This also allows for more complex phenomena, like orbital precession. This is extremely hard to simulate based on a naively plugged Newton's formula for gravitational attraction, especially if the system is designed to isolate all interactions into two-body problems.



And while circles are simple enough, ellipses are a HELL to work with. They're not that much different from the circles -- VISUALLY -- you only need one more parameter, right, and that's the horizontal stretching. Oh but now you also need rotation because stretching breaks symmetry. But the real problem lies in the fact that the simulation can give you only points, so it's a matter of fitting a parametric ellipse to existing data.

Implicit quadratic form of an ellipse has 6 polynomial coefficients (ax^2+2bxy+cy^2+2dx+2fy+g=0) and solving a system of equations which emerges when trying to fit existing data, requires knowledge about Hilbertian matrices, eigenvalues/eigenvectors, null space, and yes it's true, I'm just throwing funny words at you. The only legit code I've been able to scour is either GIS, data visualization, image recognition and similarly ultra-dedicated jaw-dropping pieces of someone's academic work, typically relying on heavy libraries designed for hard core linear algebra solvers and complex plane shenanigans. Not to mention that all of it is either OpenCV (cpp; relatively intelligible but not useful), NumPy (python; S*** I can't read this AT ALL).

So what do we do? Well, we do it ourselves lol. Let's do some shmackademic work!

Truth be told, we have a trick in our sleeve, because what we do isn't fitting to some noisy data sample resembling a curve. We actually have a neat elliptic path defined through ordered points. But that's not all! Remember that the bodies always revolve around the world origin in the simulation? This means that we can easily produce an absolute orbital path for each body, instead of having to analyze all kinds of floral designs and Lissajous patterns of relative prograde/retrograde motion (video). By Kepler's 1st law, it just happens that the parent body (the gravity well) sits exactly in one focus of the ellipse (or it's the center of a circle, but I wouldn't bet on it, as true circles are rare). So we know that one focus at all times, and if we could only find the true center, we would have it all.



What we need are the following 2D parameters:
a) true elliptic center (right in-between the two foci),
b) axis lines,
c) rotation,
d) semi-minor axis length,
e) semi-major / semi-minor axis ratio (aka the stretching factor, analog to eccentricity, but much easier to work with).

First we need a method to analytically discern circles and ellipses from parabolas and hyperbolas. How? Well I don't know. That part is surely solvable but not really that simple or even necessary. Just make sure your paths are closed if you want parametrization. Only circles and ellipses exhibit indefinite periodic motion anyway and that's what we're after.

Ok then, once we know the path is closed, we can proceed by naively finding the average of all points (aka the finite points centroid), by summing them all up and dividing by their number. This will get us the center. Whoa, that was easy.

Now we know that the two foci always lie on the semi-major axis. We can use this to find the line of this axis. Any line can be defined by a point (vector) and a slope (scalar), however to make it more robust, I say let's do this with a point (vector) and a direction (unit vector). We know the "point" part has to be the center. I called these parameters p and k (don't ask me why).

Code (csharp):

1. void estimateAxisLines(Vector2 center, Vector2 focus, out (Vector2 p, Vector2 k) major, out (Vector2 p, Vector2 k) minor) {
2.   major = ( center, center.To(focus).normalized ); // where a.To(b) equals (b - a)
3.   minor = ( center, major.k.Perp(normalize: false) ); // axes are always orthogonal to each other
4. }

where Perp looks like this, and because we know that k was already unit, we can bypass normalization.

Code (csharp):

1. static public Vector2 Perp(this Vector2 v, bool normalize = true)
2.   => normalize? new Vector2(-v.y, v.x).normalized
3.               : new Vector2(-v.y, v.x);

Generally in 2D I rotate stuff like this (we don't need it tho, at least not for the simulation per se, but you might find it useful, it's faster than quaternions)

Code (csharp):

1. static public Vector2 Rotated(this Vector2 point, float radians)
2.   => point.Rotated(Vector2.zero, radians);
3.  
4. static public Vector2 Rotated(this Vector2 point, Vector2 center, float radians) {
5.   Vector2 trig = Polar(radians), c2p = point - center;
6.   return center + new Vector2(c2p.PerpDot(trig), c2p.Dot(trig));
7. }
8.  
9. static public Vector2 Polar(float theta) => new Vector2(MathF.Cos(theta), MathF.Sin(theta));
10.  
11. static public float PerpDot(this Vector2 lhs, Vector2 rhs)
12.   => lhs.x * rhs.y - lhs.y * rhs.x;
13.  
14. static public float Dot(this Vector2 lhs, Vector2 rhs)
15.   => lhs.x * rhs.x + lhs.y * rhs.y;

(Technically multiplying by 2x2 matrix {{cos(a), -sin(a)}, {sin(a), cos(a)}} with some origin shifting.)

This PerpDot thing is a peculiar thing, specific to 2D, similar to how Cross (⨯) is specific to 3D. In fact, Cross IS a component-wise combination of PerpDots per each plane (YZ, ZX, XY in that order). It's like a weird cousin to Dot (⋅), where the original vector is orthogonalized prior to operation. We'll need that again later.

But for now we have the axis lines, and we can easily find rotation by passing semi-major axis's k parameter to Atan2. In my case I had to flip the sign of X, but that's probably my sloppy math.

Code (csharp):

1. void someFunc() {
2.   Vectx.ComputePolar(_estMajor.k.ScaledBy(x: -1f), out _estRotation);
3. }
4.  
5. // where
6. static public Vector2 ScaledBy(this Vector2 vec, float x = 1f, float y = 1f)
7.   => new Vector2(vec.x * x, vec.y * y);
8.  
9. static public bool ComputePolar(this Vector2 v, out float polar) {
10.   polar = Mathx.Atan2(v.y, v.x);
11.   if(polar.IsNaN()) { polar = 0f; return false; } // this is the case only if
12.                      // input was NaN, but I've also made ComputeSpherical
13.                      // and these two share the same general design principle
14.                      // (don't judge me)
15.   return true;
16. }

(The reason why I call this polar is because it's one half of polar coordinates θ, the other half r being the magnitude of the vector, but I don't need that, I just want that Atan2 out of sight. The reason behind this is because I've developed a style where I intentionally push grittiness of the numeric algebra below a certain line, as this makes the code much easier to parse and harder to produce a typo with.)

Now comes the interesting part. We can find the rotated bounding box of the ellipse if we take one axis line (say semi-major) and project all the points left and right onto it. That's easy to do, here's the math (when the line is considered as infinite and defined by a point and a direction).

Code (csharp):

1. static public Vector2 ClosestPointOnLine(Vector2 lp, Vector2 dir, Vector2 point)
2.   => lp.Directed(dir, lp.To(point).Dot(dir));
3.  
4. static public Vector2 Directed(this Vector2 v, Vector2 dir, float length)
5.   => v + length * dir;

I then find the point farthest from the center. I then take this point and find the point farthest from it. For this I use custom utilities I've developed over the years (earlier version of this can be found in this thread). All of this boils down to this kind of code

Code (csharp):

1. int fIndex;
2. List<Vector2> proj = getProjected(_estMajor, _pts);
3. _estCenter.FarthestOf(out fIndex, proj); // gives us the index of the farthest point
4. proj[fIndex].FarthestOf(out _estMajorDist, proj); // gives us the extreme distance

Voila, we get our semi-major length estimation.

Do the same thing for the semi-minor axis, and now you know the sides of a rotated bounding box which is guaranteed to perfectly contain the ellipse.

However, if I did I good job here of visualizing how orbits work, you've probably noticed that there is one problem with this approach -- Kepler's 3rd law basically ruins everything for me because it introduces a velocity gradient along the path. If the body accelerates when near the gravity well focus (periapsis) and slows down when away from it (apoapsis) this means our simulation does not produce evenly distributed points along the path. This messes up our naive approach to finding the center of the ellipse, because averaging favors the side with more densely packed points, and introduces a hefty bias.



Let's do something else then. Instead of averaging the points, let's find a true centroid of a simple polygon. "Simple" means it's not intersecting with itself (that would be a "complex" polygon), and what I've been calling "an ellipse" all this time is honestly just a convex polygon made out of ordered and finite points, which is why this works.

It's a little bit more involved, but not too much. We previously had this (wiki)

Code (csharp):

1. Vector2 estimateCloudCentroid(List<Vector2> pts) {
2.   var sum = Vector2.zero;
3.   for(int i = 0; i < pts.Count; i++) sum += pts[i];
4.   return sum / pts.Count;
5. }

and now we can upgrade to this (wiki)

Code (csharp):

1. Vector2 computePolygonCentroid(List<Vector2> pts) {
2.   var sum = Vector2.zero;
3.  
4.   var area = 0f; // aka area accumulator
5.   for(int i = 0, j = pts.Count - 1; i < pts.Count; i++) {
6.     var t = pts[j].PerpDot(pts[i]);
7.     area += t;
8.     sum += t * (pts[i] + pts[j]);
9.     j = i;
10.   }
11.  
12.   if(area.IsZero(1E-7f)) return estimateCloudCentroid(pts);
13.   return sum / (3f * area);
14. }
15.  
16. static public bool IsZero(this float value, float epsilon = kEPSILON_32)
17.   => Abs(value) < epsilon;

Notice the perpdot? This particular application is called a shoelace formula.
I keep my kEPSILON_32 at 1E-5f btw, but it could be finer. This is to avoid issues with 32-bit floating point precision.

Centroids in general are supposed to be imaginary points where you can support the shape with a pencil and it will stay in balance (assuming its weight is evenly distributed over area). With convex shapes, it is guaranteed that the centroid will be contained by the shape. Our naive approach had a bias, however with polygons we take the area into account, so it doesn't matter how well distributed the points are.

And that's all. You have now successfully parametrized a physically simulated elliptic path. The only thing that's left is to mark up the velocity gradient, so that the body can pretend to be affected by gravity as it orbits. You are now free to accelerate (or decelerate) your simulation as much as you need to and it will also resist any sudden spikes in performance. Of course, we've sacrificed the dynamic part of the simulation, but for many practical purposes, no one really needs to move entire planets off course, right?

Well, sometimes.

[edit]
Oh, and I forgot to show how to actually sample a parametrized ellipse.

First we start with the obvious. It's kind of circular, so there must be some polar coordinates involved. t is a simple interpolant 0..1 that corresponds to θ, but saves us from having to think about radians. We can think about revolutions instead. A full revolution of 1 lands on the same place as 0, so this value is also modular.

Code (csharp):

1. Vector2 getPoint(float radius, float t) => radius * Polar(t.Mod(1f) * TAU);

Right now this is just a circle centered at world origin. To turn it into an ellipse we change radius to semiMinorAxis and introduce a stretching factor (a ratio). The rules of my parametrization is that ratio can never get below 1 (practically semiMinorAxis is always the minor one, or they are both equal).

Code (csharp):

1. // ratio must be >= 1
2. Vector2 getPoint(float semiMinorAxis, float ratio, float t)
3.   => (semiMinorAxis * Polar(t.Mod(1f) * TAU)).ScaledBy(y: majorRatio);

Mod or 'modulo' (not the same thing as

%

) can be defined as

Code (csharp):

1. // m is supposed to be > 0
2. static public float Mod(this float n, float m)
3.   => (m <= 0f)? 0f : (n %= m) < 0f? n + m : n;

(This provides plenty of benefits when working with cyclic structures, which are typical for all kinds of modular geometry (vertices, edges, loops, angles etc.), but also with data arrays. For example 7.Mod(7) would yield 0, allowing you to cross the boundary and start the same sequence again, but there is also -1.Mod(7) which would yield 6, making it easy to refer to a last element without having to invoke an upper bound etc.)

Now we should introduce both rotation and translation so that it becomes independent from the world origin. But if we think about it, we already did scale, and intend to do rotation and translation. That's standard transformation.

So to make this more universal (and conveniently separable if we ever need the straightforward non-transformed geometry for some reason), let's add TRS method (short for translation/rotation/scale). By convention/intuition the three operations are performed in reversed order, so first go the scalar multipliers, followed by matrix rotation, then we offset the whole thing. (A great benefit when doing this on your own is that you can easily extend the process to batch transform a load of points, with precomputed trigonometry.)

Code (csharp):

1. // rotation in radians
2. Vector2 trs(Vector2 pt, Vector2 scale, float rotation, Vector2 translate)
3.   => translate + pt.ScaledBy(scale).Rotated(rotation);
4.  
5. // rotation in radians; ratio must be >= 1
6. Vector2 getPoint(Vector2 center, float semiMinorAxis, float ratio, float rotation, float t)
7.   => trs((semiMinorAxis * Polar(t.Mod(1f) * TAU)), new Vector2(1f, ratio), rotation, center);

Finally, here's how to get the two foci points.

Code (csharp):

1. // rotation in radians; ratio must be >= 1
2. (Vector2 focus1, Vector2 focus2) getFoci(Vector2 center, float semiMinorAxis, float ratio, float rotation) {
3.   if((ratio - 1f).IsZero()) return ( center, center ); // this must be a circle
4.   var f = MathF.Sqrt((ratio * semiMinorAxis).Sqr() - semiMinorAxis.Sqr());
5.   return ( // the scaling is already integrated above for better precision
6.     trs(new Vector2(0f, -f), Vector2.one, rotation, center),
7.     trs(new Vector2(0f, +f), Vector2.one, rotation, center)
8.   );
9. }

[/edit]

Btw here are some more details on the actual semi-implicit Euler implementation.
This is the core update routine, for example.

Code (csharp):

1. void FixedUpdate() {
2.   recurseUpdate(_bt, Time.fixedDeltaTime * _timeScale); // _bt is the root of 'body tree'
3. }
4.  
5. void recurseUpdate(BodyNode body, float deltaTime) {
6.   updateMotionOf(body, deltaTime);
7.   if(body.satellites != null)
8.     for(int i = 0; i < body.satellites.Count; i++)
9.       recurseUpdate(body.satellites[i], deltaTime);
10. }
11.  
12. void updateMotionOf(BodyNode body, float deltaTime) {
13.   if(body.parent != null) { // only localPos is affected
14.     body.vel += g_acc(-body.localPos, body.parent.mass) * deltaTime;
15.     body.localPos += body.vel * deltaTime;
16.   }
17.  
18.   // true global pos will integrate the parents' offsets up the chain
19.   body.gameObject.transform.position = body.pos.ToVector3_XZ();
20. }

This is how my data is organized (ymmv). I start with linearly serializable BodyDef structs

Code (csharp):

1. [SerializeField] [NonReorderable] BodyDef[] _bodies;
2.  
3. // ...
4.  
5. [Serializable]
6. private struct BodyDef {
7.  
8.   public BodyDefType type;
9.   public bool showProjection; // I use this to selectively toggle path gizmos in edit mode
10.   public int parent; // just an index
11.   [Min(1f)] public float mass; // in petagrams?
12.   [ColorUsage(showAlpha: false)] public Color color;
13.   public Vector2 position; // initial position, some coordinate relative to the parent
14.   public Vector2 velocity; // typically tangential velocity, or the body would free-fall
15.  
16.   public BodyDef(BodyDefType type, int parent, float mass, Color color, Vector2 position, Vector2 velocity)
17.     => (this.type, this.showProjection, this.parent, this.mass, this.color, this.position, this.velocity)
18.       = (type, false, parent, mass, color.ToColor(alpha: 1f), position, velocity);
19. }
20.  
21. private enum BodyDefType {
22.   Barycenter,
23.   Stellar,
24.   Substellar
25. }

This data is then reorganized to a tree structure

Code (csharp):

1. private class BodyNode {
2.  
3.   public GameObject gameObject;
4.   public bool isBarycenter; // I don't instantiate a prefab if the gravity well is a mathematical point
5.   public bool isEmissive; // I use this to configure the sun's material
6.   public Color color;
7.   public BodyNode parent; // actual hierarchy is replicated
8.   public List<BodyNode> satellites;
9.   public float mass;
10.   public Vector2 localPos; // the value that is updated by the sim
11.   public Vector2 pos => parent is null? localPos : localPos + parent.pos; // the true position is computed on the fly
12.   public float scale => ComputeLogScale(mass);
13.   public Vector2 vel; // the value that is updated by the sim
14.  
15.   // I use this just to visualize the wildly ranging masses
16.   // Here are some examples of masses: The Sun = 230, The Earth = 12, The Moon = 1
17.   static public float ComputeLogScale(float val) => (1f + Mathx.Log(val)) * .1f;
18.  
19. }

When converting from an array to a tree I also need to make sure I didn't accidentally made a cyclic graph. Trees are acyclic by definition, and besides, running a loop over a cyclic graph would halt Unity. This was just an experiment, so I'm not too strict about this, but here's one piece of code that tries to prevent this from happening. This method behaves recursively, but is written in a non-recursive manner, and attempts to find a world position of an object while in edit mode, because I draw gizmos and whatnot, which lets me inspect some behaviors before I set everything in motion.

Code (csharp):

1. bool getPositionOf(int index, out Vector2 pos) {
2.   pos = Vector2.zero;
3.   int count = 0;
4.   while(true) { // <- usually a tell-tale sign that dragons dwell inside
5.     if(index == -1) return true; // -1 is a special signal used by the top-level body, this means it's finished
6.     if(index < 0 || index >= _bodies.Length) return false; // index was out of bounds, report failure
7.     pos += _bodies[index].position; // accumulate local positions up the hierarchy
8.     index = _bodies[index].parent; // re-index in order to process the parent next
9.     if(++count > 20) break; // this prevents a deadlock; I like to keep safety mechanisms at the bottom
10.                                   // do..while also makes sense
11.   }
12.   throw new InvalidOperationException("Cyclic graph error.");
13. }

This is how I build the tree

Code (csharp):

1. List<BodyNode> recurseTree(BodyDef[] defs, BodyNode parent = null, int target = -1) {
2.   if(target < -1) throw new NotSupportedException($"Invalid index '{target}'.");
3.  
4.   var result = new List<BodyNode>();
5.   for(int i = 0; i < defs.Length; i++) {
6.     if(defs[i].parent == target) {
7.       var node = makeFromDef(defs[i]);
8.       node.parent = parent;
9.       node.satellites = recurseTree(defs, node, i);
10.       result.Add(node);
11.     }
12.   }
13.  
14.   // here we make sure this isn't a so-called forest, but a single tree
15.   if(target == -1 && result.Count > 1)
16.     throw new NotSupportedException("System can't have more than one dominant gravity well.");
17.  
18.   return result;
19. }
20.  
21. BodyNode makeFromDef(BodyDef def)
22.   => new BodyNode() {
23.     isBarycenter = def.type == BodyDefType.Barycenter,
24.     isEmissive = def.type == BodyDefType.Stellar,
25.     color = def.color,
26.     mass = def.mass,
27.     localPos = def.position,
28.     vel = def.velocity
29.   };

This is all put into motion by doing

Code (csharp):

1. void OnEnable() => _bt = recurseTree(_bodies)[0];

I then recursively instantiate my prefabs in Start, similar to how Update works. At this point you've seen all the relevant pieces of the puzzle.

That's all, hopefully this research helps you make a fantastic game.

 

Thank you so much, this was a very interesting read.

 

+1 really great writeup. I hope it bubbles to the top of google searches.

I still remember profound disappointment by the instability in my first such simulation, way back on a Commodore Vic-20 in BASIC... and that was SINGLE-body simulation! The pixel always went flying or else fell into the center. Grr.

 

Thanks!

I feel it's more about the code itself than it is about physics. I'm not a physicist myself, and I've learned math beyond anything I ever imagined only because I wanted to be able to make competent games. I would hate if I'd be heavily scrutinized for my shmackademic approach, mostly because that would miss the point entirely.

I trust I'm not alone in this and that people are more likely to approach the problem from the implementation perspective, instead of chasing some abstract purity. Semi-implicit Euler is childishly simple, it just happens that it also provides a rigorous physical approximation -- which is why I think this forum is perfect: maybe exactly to encourage anyone who wouldn't dare to go to physics forum in the first place.

I also made a mistake once and tried to participate in the game design sub forum. Let's just say I find the crowd here much more pleasant. Maybe it's the daily volume of threads that helps this forum feel less like someone else's turf.