# Orbit equation solver giving wrong orbits

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. However, inputting Earth's Keplerian elements found here, 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;

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;

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.

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?

• 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. 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
Commented Jun 16 at 4:32
• 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
Commented Jun 16 at 4:33
• 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. Commented Jun 16 at 14:26
• I'm looking forward to hearing if that works or not, please let us know!
– uhoh
Commented Jun 16 at 14:28

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.
float i = Mathf.Repeat(inclination + pi/2, tau);