2
$\begingroup$

I have done an extensive amount of research into this topic, but am really struggling to implement the math necessary to pull this off correctly.

Just for some context: I'm trying to create a semi-realistic 2D simulation of Keplerian orbits in Unity. To achieve this, I have two functions which are able to convert the position and velocity of an object into the necessary orbital elements, and then convert back into position and velocity for a given time around that orbit.

This works perfectly while the orbit is elliptical, but as soon as the eccentricity goes above 1 it all breaks. I've tried a few solutions suggested here and by chatGPT, but no amount of tweaking gets the result I'm looking for. I plan on capping the sphere of influence of each celestial object as a function of its gravitational parameter, so I'm not looking to simulate the entirety of a hyperbolic/parabolic orbit. I only want this feature so I can relatively accurately reproduce complex gravitational manoeuvres like slingshots and orbital captures.

Also, I'm well aware that the specific reason behind the errors is that I'm trying to feed invalid Vector2 coordinates into the Line Renderer. My issue is the fact that those values are invalid in the first place, my code (in an ideal world) should be able to handle hyperbolic orbits and output a valid Vector2 position at a given time.

Here is my code so far as well as some screenshots:

public class Orbit
{
    public info primaryBody;
    public float semiMajorAxis;
    public float semiMinorAxis;
    public float eccentricity;
    public float longitudePeriapsis;
    public float period;
    public float trueAnomaly;

    public void UpdateElements(Vector2 position, Vector2 velocity)
    {
        position = (Vector3)position - primaryBody.transform.position;
        velocity = velocity - primaryBody.rb.velocity;

        float r = position.magnitude;
        float v = velocity.magnitude;
        float mu = Universe.gravitationalConstant * primaryBody.mass;
        float specificEnergy = (v * v) / 2 - mu / r;

        semiMajorAxis = -mu / (2 * specificEnergy);
        semiMinorAxis = semiMajorAxis * Mathf.Sqrt(1 - eccentricity * eccentricity);

        float angularMomentum = r * v;
        eccentricity = Mathf.Sqrt(1 + (2 * specificEnergy * angularMomentum * angularMomentum) / (Mathf.Pow(mu, 2)));

        float dot = Vector2.Dot(position, velocity);
        float cosNu = dot / (r * v);
        float sinNu = Mathf.Sqrt(1 - cosNu * cosNu);
        trueAnomaly = Mathf.Atan2(sinNu, cosNu);
        if (dot < 0)
        {
            if (trueAnomaly >= 0) trueAnomaly -= Mathf.PI;
            else trueAnomaly += Mathf.PI;
        }
        trueAnomaly *= Mathf.Rad2Deg;

        Vector2 eccentricityVector = ((r * v * v - mu) / mu) * position.normalized;
        float w = Mathf.Atan2(eccentricityVector.y, eccentricityVector.x);
        longitudePeriapsis = Mathf.Rad2Deg * w;

        float periodSquared = (4 * Mathf.PI * Mathf.PI * Mathf.Pow(semiMajorAxis, 3)) / mu;
        period = Mathf.Sqrt(periodSquared);
    }

    public Vector2 CalculatePositionOnOrbit(float t)
    {
        float a = semiMajorAxis;
        float e = eccentricity;
        float w = longitudePeriapsis;
        float mA = 2 * Mathf.PI * t; // Calculate the mean anomaly at time t.
        
        //Newton-Raphson method as no direct solution
        float E = mA; //Eccentric anomaly
        float prevE;
        float i = 0;
        do
        {
            prevE = E;
            E = mA + e * Mathf.Sin(prevE);
            i += 1;
        } 
        while (Mathf.Abs(E - prevE) > 1e-6f && i < 100);

        float r = a * (1 - e * Mathf.Cos(E)); //distance to the primary body
        float nu = 2 * Mathf.Atan2(Mathf.Sqrt((1 + e) / (1 - e)) * Mathf.Sin(E / 2), Mathf.Cos(E / 2));  //True anomaly

        float theta = nu + Mathf.Deg2Rad * w;

        float x = r * Mathf.Cos(theta);
        float y = r * Mathf.Sin(theta);

        return new Vector2(x, y) + (Vector2)primaryBody.transform.position; ;
    }
}

Circular-ish orbit

Elliptical orbit with off-centre longitude of periapsis

Broken Hyperbolic orbit

EDIT:

I've since made some changes to the code which I believe to be correct, but the function still doesn't output what is expected. To improve readability, I've copied over a lot of the logic from the elliptical case to the hyperbolic case.

In addition, I've also added the parameter "meanMotion" to the overall class which allows me to adjust the orbit in time increments of seconds, rather than an arbitrary value between 0 and 1. As such, the value "t" now refers to seconds since epoch.

I am unsure about what to try next as I was pretty confident this would work.

public void UpdateElements(Vector2 position, Vector2 velocity)
    {
        position = (Vector3)position - primaryBody.transform.position;
        velocity = velocity - primaryBody.rb.velocity;

        float r = position.magnitude;
        float v = velocity.magnitude;
        float mu = Universe.gravitationalConstant * primaryBody.mass;
        float specificEnergy = (v * v) / 2 - mu / r;

        semiMajorAxis = -mu / (2 * specificEnergy);
        semiMinorAxis = semiMajorAxis * Mathf.Sqrt(1 - eccentricity * eccentricity);

        float angularMomentum = r * v;
        eccentricity = Mathf.Sqrt(1 + (2 * specificEnergy * angularMomentum * angularMomentum) / (Mathf.Pow(mu, 2)));

        Vector2 eccentricityVector = ((r * v * v - mu) / mu) * position.normalized;
        float w = Mathf.Atan2(eccentricityVector.y, eccentricityVector.x);
        longitudePeriapsis = Mathf.Rad2Deg * w;

        period = Mathf.PI * 2 * Mathf.Sqrt(Mathf.Pow(semiMajorAxis, 3) / mu);
        meanMotion = Mathf.Sqrt(mu / Mathf.Pow(Mathf.Abs(semiMajorAxis), 3));
    }

    public Vector2 CalculatePositionOnOrbit(float t)
    {
        float a = semiMajorAxis;
        float e = eccentricity;
        float w = longitudePeriapsis;
        float mA = meanMotion * t;

        if (e == 1) e += Mathf.Epsilon;

        if (e < 1)
        {
            //Newton-Raphson method as no direct solution
            float E = mA; //Eccentric anomaly
            float prevE;
            float i = 0;
            do
            {
                prevE = E;
                E = mA + e * Mathf.Sin(prevE);
                i += 1;
            }
            while (Mathf.Abs(E - prevE) > 1e-6f && i < 100);

            float r = a * (1 - e * Mathf.Cos(E)); //distance to the primary body
            float nu = 2 * Mathf.Atan2(Mathf.Sqrt((1 + e) / (1 - e)) * Mathf.Sin(E / 2), Mathf.Cos(E / 2));  //True anomaly in radians

            float theta = nu + Mathf.Deg2Rad * w;

            float x = r * Mathf.Cos(theta);
            float y = r * Mathf.Sin(theta);

            return new Vector2(x, y) + (Vector2)primaryBody.transform.position;
        }
        else
        {
            //Newton-Raphson method as no direct solution
            float E = mA; //Eccentric anomaly
            float prevE;
            float i = 0;
            do
            {
                prevE = E;
                E = mA + e * Mathf.Sin(prevE);
                i += 1;
            }
            while (Mathf.Abs(E - prevE) > 1e-6f && i < 100);
            float H = e * math.sinh(E) - E;

            float r = (a * (1 - e * e)) / (1 + e * Mathf.Cos(H)); //distance to the primary body
            float nu = 2 * Mathf.Atan2(Mathf.Sqrt(e + 1) * math.sinh(H / 2), Mathf.Sqrt(e - 1) * math.cosh(H / 2)); //True anomaly in radians

            float theta = nu + Mathf.Deg2Rad * w;

            float x = r * math.cosh(theta);
            float y = r * math.sinh(theta);

            return new Vector2(x, y) + (Vector2)primaryBody.transform.position;
        }
    }
```
$\endgroup$

1 Answer 1

4
$\begingroup$

Things wind up changing when you look at hyperbolic orbits instead of elliptical ones when going from cartesian vectors to Keplerian elements, because when you use the elliptical equations some of the values flip into complex domains. So generally, when I write programs to calculate time of flight on hyperbolic orbits, I use the hyperbolic versions of the Keplerian equations.

Once you've got your true anomaly ($f$) calculated, you need to calculated the Hyperbolic anomaly, rather than the eccentric anomaly.

The Eccentric anomaly is the area (in units of semi-major axis -squared) of the circular sector of the osculating circle (a circle with the same center as the orbit, and radius equal to the semimajor axis) of an elliptical orbit that goes from the direction of periapsis, to the center of the orbit, to the vertical projection the object's position parallel to the minor axis on that osculating circle. |Elliptic Orbit Anomalies| |-| |Elliptic orbit anomalies}| |GeoGebra Graph| |The orbit is in Green. The angle of True anomaly is shaded in green, the osculating circle is in blue and the area of the Eccentric Anomaly is shaded blue. The hypothetical mean anomaly object's orbit is in red, with the mean anomaly the shaded area.|

For a hyperbolic orbit, you need the Hyperbolic anomaly, which is analagous, except it's a hyperbolic sector projected onto an equilateral hyperbola with the same semi-major axis as your orbit.

Hyperbolic orbit and Anomalies
enter image description here
GeoGebra Graph
Hyperbolic orbit is in black. Asymptotes and phantom lobe are dotted, True anomaly is marked in black. Equilateral hyperbola is in blue. Sector that is half the area of the region designated by the Hyperbolic anomaly is in blue. The hypothetical mean anomaly object's orbit is in red, with the mean anomaly the shaded area.

Calculating Hyperbolic anomaly $H$ from True Anomaly $\nu$:

$$H = \mathrm{atanh}\left(\sqrt{\frac{e -1}{e + 1}}\right)\cdot \tan\frac{\nu}{2}$$

If going in the other direction, True Anomaly from Hyperbolic Anomaly: $$\nu = \mathrm{atan2}(\sqrt{e + 1}\cdot \sinh\frac{H}{2},\sqrt{e - 1}\cdot \sinh\frac{H}{2} )$$

This will wind up in the range $(-\infty, \infty)$

Once you've got that, you can get Mean Anomaly $M$, which is almost the same. It's still a measurement of time since or to periapsis in terms of the swept out area (in units of semi-major-axis-squared) of a hypothetical body in a circular orbit whose orbital radius is the absolute value of your semi-major axis. The equation, however, flips signs, and uses hyperbolic trigonometric functions. Mean anomaly will have valid values in the range $(-\infty, \infty)$

Calculating Mean anomaly from hyperbolic anomaly $(-\infty, \infty)$: $$M = e * \sinh H - H$$

Going in the other direction, you'll have to use Newton's method on these as well, but using the above equation and its first derivative in the algorithm.

From there, to find position as a function of time, you do the same type of mean motion calculation you would to set your mean anomaly:

$$ M = M_0 + n * t $$

Where $M_0$ is mean anomaly at epoch, $n$ is the mean motion $ n= \sqrt{\frac{\mu}{|a|^3}}$, and $t$ is the time since epoch, and work backwards.

Which brings me to another issue; if in the function Vector2 CalculatePositionOnOrbit(float t), your value of t is coming in as an ordinary measurement of time, you're definitely calculating your mean anomaly incorrectly. As mentioned above, Mean Anomaly is both dependent on the time since periapsis and the semimajor axis of the orbit, you cannot calculate it from time since periapsis alone.

$\endgroup$
10
  • $\begingroup$ Thank you for your very swift reply, just a couple questions and clarifications. In my equations, I always wind up using the true anomaly and the distance between both bodies to calculate the position vector, are you saying to basically swap out the true anomaly for hyperbolic anomaly (when eccentricity > 1), and use cosh and sinh instead? $\endgroup$
    – Ethan
    Oct 12, 2023 at 6:49
  • $\begingroup$ As to your other issue, "t" is a value between 0 and 1, 0 being at periapsis, 1 being at periapsis again, but after a full period. I believe I've mislabelled my variables, the mA calculation is not the actual mean anomaly, it's actually the mean motion. I have quite a basic understanding of the maths and have been working to translate other people's examples into simplified c# code, so forgive me for making silly errors like that. $\endgroup$
    – Ethan
    Oct 12, 2023 at 7:03
  • 1
    $\begingroup$ @Ethan With regards to the hyperbolic equations, sometimes you substitute in the hyperbolic trig functions, and sometimes you also flip signs on terms in the equations, sometimes you use different equations entirely. And then there's the parabolic case, which is directly solvable as a function of time. Ultimately, I really found the Wikipedia articles very helpful when hashing through these equations myself before KSP had a map screen. $\endgroup$
    – notovny
    Oct 13, 2023 at 1:07
  • $\begingroup$ With regards to t, yeah, that's not mean motion, it's just "fraction of orbital period elapsed since the last periapsis," then. Also, it doesn't look like you're doing Newton's Method in the "newton's method" section, and I don't have the math chops to know at a glance if it will still converge if you substitute in the hyperbolic anomaly equation. $\endgroup$
    – notovny
    Oct 13, 2023 at 1:10
  • $\begingroup$ Yeah, as I said I dont have a great understanding of the maths behind it, so I've mainly been Frankensteining other people's code. I'll read through what you've sent as well as the Wikipedia articles and see if I can generate a solution. When I do I'll post it here so others in future can see. Thanks for your help $\endgroup$
    – Ethan
    Oct 13, 2023 at 8:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.