<p><strong>MAJOR EDIT!</strong></p>
<p>Upon further inspection, it appears the actual issue is that the axes of the rotation matrix are incorrect. Swapping r_z and r_y in the <code>GetPositionInReferenceFrame</code> function fixes the issue entirely. I believe this isn't an issue with the equation itself, rather that in the paper the z direction points up, while in Unity it points to the side. I'll leave the original solution up in case anyone has this issue.</p>
<p>As per user uhoh,</p>
<blockquote>
<p>Therefore, the variable i for inclination may be 0 at the north pole
and π at the south and π/2 at the equator. I'll bet if you use 90° or π/2
for your &quot;inclination&quot; you'll get an equatorial orbit.</p>
</blockquote>
<p>Adding this line of code at the beginning of the main function solved the problem:<br />
<code>float i = Mathf.Repeat(inclination + pi/2, tau);</code></p>
<p>Then you can use i for the rotation matrix instead of inclination.</p>

CC BY-SA 4.0
2024-06-16T14:30:27.797Z
2024-06-16T16:21:38.717Z
2024-06-16T16:21:38.717Z
Eduardo Dieppa
<p>Recently I've become interested in simulating orbital mechanics using Kepler's equations, and in Unity I created my own equation solver by following <a href="https://downloads.rene-schwarz.com/download/M001-Keplerian_Orbit_Elements_to_Cartesian_State_Vectors.pdf" rel="nofollow noreferrer">this paper I found.</a> However, inputting Earth's Keplerian elements found <a href="https://ssd.jpl.nasa.gov/planets/approx_pos.html" rel="nofollow noreferrer">here,</a> it gave me a polar orbit instead of what I was expecting. I'm aware that some of the elements were given in degrees, and I converted them to radians for my solver.</p>
<p>Here's the code for the solver, and yes, I've triple-checked that the huge rotation matrix in <code>GetPositionInReferenceFrame</code> is correct:</p>
<pre><code>public class OrbitingScript : MonoBehaviour
{
    GameObject orbitingObject;

    [Header(&quot;Keplerian Elements&quot;)]
    public float semiMajorAxis = 10;
    [Range(0, 1)]
    public float eccentricity = 0;
    [Range(0, Mathf.PI * 2)]
    public float argOfPeriapsis = 0;
    [Range(0, Mathf.PI * 2)]
    public float ascendingNode = 0;
    [Range(0, Mathf.PI * 2)]
    public float inclination = 0;
    [Range(0, Mathf.PI * 2)]
    public float meanAnomaly0 = 0;

    [Header(&quot;Simulation Parameters&quot;)]
    public float centralBodyMass = 10;
    public float consideredEpoch = 0;
    public float gravityConst = 1;
    public float scale;

    [HideInInspector]
    public float standardGrav;

    float pi = Mathf.PI;
    float tau = Mathf.PI * 2;

    private void OnValidate()
    {
        gameObject.transform.position =  GetPositionAndVelocityAtEpoch().Item1 / scale;
    }

    // Start is called before the first frame update
    void Start()
    {

    }

    // Update is called once per frame
    void Update()
    {
        
    }

    public (Vector3, Vector3) GetPositionAndVelocityAtEpoch()
    {
        standardGrav = gravityConst * centralBodyMass;
        float meanAnomaly = CalculateMeanAnomaly();
        float eccAnomaly = CalculateEccentricAnomaly(meanAnomaly, 5);
        float trueAnomaly = CalculateTrueAnomaly(eccAnomaly);
        float distance = GetDistanceToCenter(eccAnomaly);
        Vector3 orbitalPosition = GetPositionInOrbitalFrame(trueAnomaly, distance);
        Vector3 referencePosition = GetPositionInReferenceFrame(orbitalPosition);
        return (referencePosition, Vector3.zero);
    }

    public float CalculateMeanAnomaly()
    {
        if (consideredEpoch == 0)
        {
            return meanAnomaly0;
        } else
        {
            float d_t = 86400 * consideredEpoch;
            float unclampedMeanAnomaly = meanAnomaly0 + d_t * Mathf.Sqrt(standardGrav / Mathf.Pow(semiMajorAxis, 3));
            return Mathf.Repeat(unclampedMeanAnomaly, tau);
        }
    }

    public float CalculateEccentricAnomaly(float meanAnomaly, int decimalPlaces)
    {
        float eccAnomaly;
        if (eccentricity &lt; 0.8) eccAnomaly = meanAnomaly; else eccAnomaly = pi;
        float delta = Mathf.Pow(10, -decimalPlaces);
        float F = eccAnomaly - eccentricity * Mathf.Sin(meanAnomaly) - meanAnomaly;
        int maxIterations = 30;
        int i = 0;
        while (Mathf.Abs(F) &gt; delta &amp;&amp; i &lt; maxIterations) 
        {
            eccAnomaly = eccAnomaly - F / (1 - eccentricity * Mathf.Cos(eccAnomaly));
            F = eccAnomaly - eccentricity * Mathf.Sin(eccAnomaly) - meanAnomaly;
            i = i + 1;
        }

        return Mathf.Round(eccAnomaly * Mathf.Pow(10, decimalPlaces)) / Mathf.Pow(10, decimalPlaces);
    }

    public float CalculateTrueAnomaly(float eccAnomaly)
    {
        float y = Mathf.Sqrt(1 + eccentricity) * Mathf.Sin(eccAnomaly / 2);
        float x = Mathf.Sqrt(1 - eccentricity) * Mathf.Cos(eccAnomaly / 2);

        return 2 * Mathf.Atan2(y, x);
    }

    public float GetDistanceToCenter(float eccAnomaly)
    {
        return semiMajorAxis * (1 - eccentricity * Mathf.Cos(eccAnomaly));
    }

    public Vector3 GetPositionInOrbitalFrame(float trueAnomaly, float distance)
    {
        float o_x = distance * Mathf.Cos(trueAnomaly);
        float o_y = distance * Mathf.Sin(trueAnomaly);
        float o_z = 0;
        return new Vector3(o_x, o_y, o_z);
    }

    public Vector3 GetPositionInReferenceFrame(Vector3 orbitalPosition)
    {

        float r_x = orbitalPosition.x * (Mathf.Cos(argOfPeriapsis) * Mathf.Cos(ascendingNode) - Mathf.Sin(argOfPeriapsis) * Mathf.Cos(inclination) * Mathf.Sin(ascendingNode)) - orbitalPosition.y * (Mathf.Sin(argOfPeriapsis) * Mathf.Cos(ascendingNode) + Mathf.Cos(argOfPeriapsis) * Mathf.Cos(inclination) * Mathf.Sin(ascendingNode));
        float r_y = orbitalPosition.x * (Mathf.Cos(argOfPeriapsis) * Mathf.Sin(ascendingNode) + Mathf.Sin(argOfPeriapsis) * Mathf.Cos(inclination) * Mathf.Cos(ascendingNode)) + orbitalPosition.y * (Mathf.Cos(argOfPeriapsis) * Mathf.Cos(inclination) * Mathf.Cos(ascendingNode) - Mathf.Sin(argOfPeriapsis) * Mathf.Sin(ascendingNode));
        float r_z = orbitalPosition.x * (Mathf.Sin(argOfPeriapsis) * Mathf.Sin(inclination)) + orbitalPosition.y * (Mathf.Cos(argOfPeriapsis) * Mathf.Sin(inclination));

        return new Vector3(r_x, r_y, r_z);
    }
}
</code></pre>
<p>Here's an image of the solver at t=60.
<a href="https://i.sstatic.net/37yN9KlD.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/37yN9KlD.png" alt="unity solver screen" /></a></p>
<p>I'm not an astronomer and I'm very new to these equations and elements, so I know there's a good chance I'm misunderstanding something. Any ideas how to fix this?</p>

CC BY-SA 4.0
2024-06-16T04:32:37.807Z
I think you are using equations 9 and 10, which according to the footnotes are from Standish, E. Myles; Williams, James G.: "Orbital Ephemerides of the Sun, Moon and Planets" which looks like it's written in [spherical coordinates](https://en.wikipedia.org/wiki/Spherical_coordinate_system). Therefore, the variable *i* for inclination may be 0 at the north pole and $\pi$ at the south and $\pi/2$ at the equator. I'll bet if you use 90° or $\pi/2$ for your "inclination" you'll get an equatorial orbit.
uhoh
CC BY-SA 4.0
2024-06-16T04:33:56.683Z
and if you add another 180° or $\pi$ to it you'll get the same orbit but going backwards. If that works, please feel free to post it as an answer. *Welcome to Stack Exchange!*
uhoh
CC BY-SA 4.0
2024-06-16T14:26:18.737Z
Thank you so much! I assumed i=0 would make a flat orbit because that's how everyone describes it, but I guess it's a quirk of the equations that I chose.
Eduardo Dieppa
CC BY-SA 4.0
2024-06-16T14:28:20.477Z
I'm looking forward to hearing if that works or not, please let us know!
uhoh
<orbit><mathematics><orbital-elements><keplers-laws>
Orbit equation solver giving wrong orbits
Eduardo Dieppa
# Orbit equation solver giving wrong orbits

(Asked by: username_0 on 2024-06-16)

Recently I've become interested in simulating orbital mechanics using Kepler's equations, and in Unity I created my own equation solver by following [this paper I found.](https://downloads.rene-schwarz.com/download/M001-Keplerian_Orbit_Elements_to_Cartesian_State_Vectors.pdf) However, inputting Earth's Keplerian elements found [here,](https://ssd.jpl.nasa.gov/planets/approx_pos.html) it gave me a polar orbit instead of what I was expecting. I'm aware that some of the elements were given in degrees, and I converted them to radians for my solver.

Here's the code for the solver, and yes, I've triple-checked that the huge rotation matrix in `GetPositionInReferenceFrame` is correct:

```
public class OrbitingScript : MonoBehaviour
{
    GameObject orbitingObject;

    [Header("Keplerian Elements")]
    public float semiMajorAxis = 10;
    [Range(0, 1)]
    public float eccentricity = 0;
    [Range(0, Mathf.PI * 2)]
    public float argOfPeriapsis = 0;
    [Range(0, Mathf.PI * 2)]
    public float ascendingNode = 0;
    [Range(0, Mathf.PI * 2)]
    public float inclination = 0;
    [Range(0, Mathf.PI * 2)]
    public float meanAnomaly0 = 0;

    [Header("Simulation Parameters")]
    public float centralBodyMass = 10;
    public float consideredEpoch = 0;
    public float gravityConst = 1;
    public float scale;

    [HideInInspector]
    public float standardGrav;

    float pi = Mathf.PI;
    float tau = Mathf.PI * 2;

    private void OnValidate()
    {
        gameObject.transform.position =  GetPositionAndVelocityAtEpoch().Item1 / scale;
    }

    // Start is called before the first frame update
    void Start()
    {

    }

    // Update is called once per frame
    void Update()
    {
        
    }

    public (Vector3, Vector3) GetPositionAndVelocityAtEpoch()
    {
        standardGrav = gravityConst * centralBodyMass;
        float meanAnomaly = CalculateMeanAnomaly();
        float eccAnomaly = CalculateEccentricAnomaly(meanAnomaly, 5);
        float trueAnomaly = CalculateTrueAnomaly(eccAnomaly);
        float distance = GetDistanceToCenter(eccAnomaly);
        Vector3 orbitalPosition = GetPositionInOrbitalFrame(trueAnomaly, distance);
        Vector3 referencePosition = GetPositionInReferenceFrame(orbitalPosition);
        return (referencePosition, Vector3.zero);
    }

    public float CalculateMeanAnomaly()
    {
        if (consideredEpoch == 0)
        {
            return meanAnomaly0;
        } else
        {
            float d_t = 86400 * consideredEpoch;
            float unclampedMeanAnomaly = meanAnomaly0 + d_t * Mathf.Sqrt(standardGrav / Mathf.Pow(semiMajorAxis, 3));
            return Mathf.Repeat(unclampedMeanAnomaly, tau);
        }
    }

    public float CalculateEccentricAnomaly(float meanAnomaly, int decimalPlaces)
    {
        float eccAnomaly;
        if (eccentricity < 0.8) eccAnomaly = meanAnomaly; else eccAnomaly = pi;
        float delta = Mathf.Pow(10, -decimalPlaces);
        float F = eccAnomaly - eccentricity * Mathf.Sin(meanAnomaly) - meanAnomaly;
        int maxIterations = 30;
        int i = 0;
        while (Mathf.Abs(F) > delta && i < maxIterations) 
        {
            eccAnomaly = eccAnomaly - F / (1 - eccentricity * Mathf.Cos(eccAnomaly));
            F = eccAnomaly - eccentricity * Mathf.Sin(eccAnomaly) - meanAnomaly;
            i = i + 1;
        }

        return Mathf.Round(eccAnomaly * Mathf.Pow(10, decimalPlaces)) / Mathf.Pow(10, decimalPlaces);
    }

    public float CalculateTrueAnomaly(float eccAnomaly)
    {
        float y = Mathf.Sqrt(1 + eccentricity) * Mathf.Sin(eccAnomaly / 2);
        float x = Mathf.Sqrt(1 - eccentricity) * Mathf.Cos(eccAnomaly / 2);

        return 2 * Mathf.Atan2(y, x);
    }

    public float GetDistanceToCenter(float eccAnomaly)
    {
        return semiMajorAxis * (1 - eccentricity * Mathf.Cos(eccAnomaly));
    }

    public Vector3 GetPositionInOrbitalFrame(float trueAnomaly, float distance)
    {
        float o_x = distance * Mathf.Cos(trueAnomaly);
        float o_y = distance * Mathf.Sin(trueAnomaly);
        float o_z = 0;
        return new Vector3(o_x, o_y, o_z);
    }

    public Vector3 GetPositionInReferenceFrame(Vector3 orbitalPosition)
    {

        float r_x = orbitalPosition.x * (Mathf.Cos(argOfPeriapsis) * Mathf.Cos(ascendingNode) - Mathf.Sin(argOfPeriapsis) * Mathf.Cos(inclination) * Mathf.Sin(ascendingNode)) - orbitalPosition.y * (Mathf.Sin(argOfPeriapsis) * Mathf.Cos(ascendingNode) + Mathf.Cos(argOfPeriapsis) * Mathf.Cos(inclination) * Mathf.Sin(ascendingNode));
        float r_y = orbitalPosition.x * (Mathf.Cos(argOfPeriapsis) * Mathf.Sin(ascendingNode) + Mathf.Sin(argOfPeriapsis) * Mathf.Cos(inclination) * Mathf.Cos(ascendingNode)) + orbitalPosition.y * (Mathf.Cos(argOfPeriapsis) * Mathf.Cos(inclination) * Mathf.Cos(ascendingNode) - Mathf.Sin(argOfPeriapsis) * Mathf.Sin(ascendingNode));
        float r_z = orbitalPosition.x * (Mathf.Sin(argOfPeriapsis) * Mathf.Sin(inclination)) + orbitalPosition.y * (Mathf.Cos(argOfPeriapsis) * Mathf.Sin(inclination));

        return new Vector3(r_x, r_y, r_z);
    }
}

```

Here's an image of the solver at t=60.
[![unity solver screen](https://i.sstatic.net/37yN9KlD.png)](https://i.sstatic.net/37yN9KlD.png)

I'm not an astronomer and I'm very new to these equations and elements, so I know there's a good chance I'm misunderstanding something. Any ideas how to fix this?

**Tags:** orbit, mathematics, orbital-elements, keplers-laws

> **Comments:**
> **username_1**: I think you are using equations 9 and 10, which according to the footnotes are from Standish, E. Myles; Williams, James G.: "Orbital Ephemerides of the Sun, Moon and Planets" which looks like it's written in [spherical coordinates](https://en.wikipedia.org/wiki/Spherical_coordinate_system). Therefore, the variable *i* for inclination may be 0 at the north pole and $\pi$ at the south and $\pi/2$ at the equator. I'll bet if you use 90° or $\pi/2$ for your "inclination" you'll get an equatorial orbit.  
> **username_1**: and if you add another 180° or $\pi$ to it you'll get the same orbit but going backwards. If that works, please feel free to post it as an answer. *Welcome to Stack Exchange!*  
> **username_0**: Thank you so much! I assumed i=0 would make a flat orbit because that's how everyone describes it, but I guess it's a quirk of the equations that I chose.  
> **username_1**: I'm looking forward to hearing if that works or not, please let us know!  

---

## 1 answers

### Answer by username_0 on 2024-06-16. Score: 3 | Accepted answer

**MAJOR EDIT!**

Upon further inspection, it appears the actual issue is that the axes of the rotation matrix are incorrect. Swapping r\_z and r\_y in the `GetPositionInReferenceFrame` function fixes the issue entirely. I believe this isn't an issue with the equation itself, rather that in the paper the z direction points up, while in Unity it points to the side. I'll leave the original solution up in case anyone has this issue.

As per user username_1,

> Therefore, the variable i for inclination may be 0 at the north pole
> and π at the south and π/2 at the equator. I'll bet if you use 90° or π/2
> for your "inclination" you'll get an equatorial orbit.

Adding this line of code at the beginning of the main function solved the problem:  
`float i = Mathf.Repeat(inclination + pi/2, tau);`

Then you can use i for the rotation matrix instead of inclination.


---

