Elliptical orbits are hard
This page contains math formulas. If it doesn't load, reload the page
Elliptical orbits
My game, i’m working on, uses circular orbits. They are really easy to compute and maintain. But they, well… circular. I want to add stretched orbits for some objects. So i made a little research about how they works. By now i only need time function, so in this article i will present equations and code required to achive it.
Elliptical orbit is an orbit with eccentricity in range 0-1. Elliptical orbits also include circular orbits (eccentricity = 0).
The problem with elliptical orbits, is ecentricity. Eccentricity determines how much orbit deviates from circular orbit. 0 for circular, 1 for parabolic trajectory and >1 for hyperbolic trajectory.
Elliptical can’t been described with \(\displaystyle{ \frac {x^{2}}{a^{2}}}+{\frac {y^{2}}{b^{2}}}=1\) like ellipse - gravity is not linear and affects velocity differently at different points.
Minor, major axis and eccentricity
Semi major axis:
\(\displaystyle a={\frac {r_{\text{max}}+r_{\text{min}}}{2}}\)
public float semiMajorAxis => (aphelion + perihelion) * .5f;Semi minor axis:
\(\displaystyle b={\sqrt {r_{\text{max}}r_{\text{min}}}}\)
public float semiMinorAxis => math.sqrt(aphelion * perihelion);Eccentricity:
\(\displaystyle e={\sqrt {1-{\frac {b^{2}}{a^{2}}}}}\)
public float eccentricity =>
math.sqrt(1 - (aphelion * perihelion) / (semiMajorAxis * semiMajorAxis));https://en.wikipedia.org/wiki/Semi-major_and_semi-minor_axes
Orbital period
Orbital period is the amount of time required to complete one orbit turn. For example, Earth have orbital period equal to 1 year, Mercury - 0.24 years and Neptune - 164.8 years.
Orbital period can be calculated as:
\(\displaystyle T=2\pi\sqrt{a^3\over{\mu}}\)
where:
- \(\mu\) - standard gravitational parameter
- a - semi-major axis
I don’t use gravitational parameters in my orbit class, so my formula will be:
public override float orbitalPeriod =>
2 * math.PI *
math.sqrt(semiMajorAxis * semiMajorAxis * semiMajorAxis) / orbitalSpeed;Polar coordinates
Polar coordinated are an easy way to work with round objects. They require angle and radius from heliocentric center. Heliocentric coordinates using Sun as a center. Angle in heliocentric polar coordinates will be true anomaly.
There are 3 types of anomalies:
- Mean (elapsed angle on circular orbit with same orbital period, if body moving at constant speed)
- Eccentric (elapsed angle on circular orbit with same orbital period)
- True anomaly (angle from sun to body)
To find true anomaly, we need to find mean anomaly and then eccentric anomaly.
Mean anomaly
Mean anomaly is a linear function. We just need to divide time by orbital period
\(\displaystyle M=\frac{2\pi t}{T}\)
public float meanMotion => 2 * math.PI / orbitalPeriod;
public double MeanAnomaly(Time time) => meanMotion * time.seconds;Eccentric anomaly
Eccentric anomaly is pretty simple:
\(\displaystyle M=E-e\sin E\)
But this is a fun part: We need to calculate \(E\) in equation, which doesn’t have closed-form solution.
Simple equation above become not so simple:
\(\displaystyle E=M+2\sum _{n=1}^{\infty }{\frac {J_{n}(ne)}{n}}\sin(nM)\)
\(J_{n}(x)\) is the Bessel function:
\(\displaystyle J_{n}(x)=\sum _{m=0}^{\infty }{\frac {(-1)^{m}}{m!\Gamma (m+\alpha +1)}}{\left({\frac {x}{2}}\right)}^{2m+\alpha }\)
I took C# bessel implementation from here.
public double EccentricAnomaly(Time time) {
double fEccentricity = eccentricity; // e
double fMeanAnomaly = MeanAnomaly(time); // M
fMeanAnomaly %= math.PI * 2; // anomalies are angles, normalize to save precision
double sum = 0;
const int iter = 32;
for (int i = 1; i <= iter; i++) {
double bessel = Bessel.J(i, i * fEccentricity);
sum += bessel / i * math.sin(i * fMeanAnomaly);
}
double eccentricAnomaly = fMeanAnomaly + 2 * sum;
return eccentricAnomaly;
}True anomaly
True anomaly is just a
2 * math.atan(math.sqrt((1 + eccentricity) / (1 - eccentricity)) *
math.tan(eccentricAnomaly * .5));Radius
public double RadiusFromEccentricAnomaly(double eccentricAnomaly) =>
semiMajorAxis * (1 - eccentricity * math.cos(eccentricAnomaly));Position
public override float2 GetPosition(Time time) {
double eccentricAnomaly = ParallelEccentricAnomaly(time);
double angle = TrueAnomalyFromEccentricAnomaly(eccentricAnomaly);
double radius = RadiusFromEccentricAnomaly(eccentricAnomaly);
return new((float)(math.sin(angle) * radius), (float)(math.cos(angle) * radius));
}Precision of coordinates
You can control performance and precision by changing iterations in eccentric anomaly calculation. Orbits with high eccentricity require more iterations.
\(\displaystyle StdDev(1-\frac{M}{E-\epsilon\sin{E}})\), \(e \approx 0.97\)
| Iterations | 2 | 4 | 8 | 16 | 32 | 64 | 128 | 256 | 512 | 1024 | 2048 | 4096 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| StdDev | 167.95% | 89.54% | 42.43% | 17.15% | 4.99% | 0.74% | 0.36% | 0.09% | 0.01% | 0.00% | 0.00% | 0.00% |
I’m using 32 iterations. It’s fine if orbit is not so stretched.
Equation is not linear because of Bessel function:
m = 2 * ((n + (int)math.sqrt(40.0 * n)) / 2);
for (j = m; j > 0; j--) {...}Multithreading
\(\sum\) can be easly parallelized. It’s a great place to use Unity job system. Note that job version have lot’s of allocations, which can be avoided in different scenarios.
| Type | Basic | Jobs |
|---|---|---|
| time | 1539ms | 188ms |
| orbits | 100k | 100k |
| iter | 32 | 32 |
[BurstCompile]
public struct EccentricAnomalyJob : IJobParallelFor {
[WriteOnly] public NativeArray<double> results;
[ReadOnly] public NativeArray<double> eccentricity;
[ReadOnly] public NativeArray<double> meanAnomaly;
public int iterations;
public void Execute(int index) {
int bodyIndex = index / iterations;
int i = (index % iterations) + 1;
double bessel = Bessel.J(i, i * eccentricity[bodyIndex]);
results[index] = bessel / i * math.sin(i * meanAnomaly[bodyIndex]);
}
}