Orbit equation solver giving wrong orbits

Question

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;

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

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?

Answer 1

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 uhoh,

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.