Converting Keplerian Orbital Elements into 2D State Vectors in Unity C#

Question

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; ;
    }
}

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;
        }
    }
```

Answer 1

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| |-| |}| |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
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.