Tengo problemas para calcular mis elementos orbitales clásicos y no sé dónde buscar

Estoy usando el motor Unity y C#, y estoy tratando de encontrar una forma de trazar futuras trayectorias orbitales 2D. Actualmente tengo una configuración de dos cuerpos. El planeta tiene una masa de 1000000, el satélite una masa de 1.

En la imagen de arriba, puede ver una línea rosa que indica el camino orbital (PASADO, NO FUTURO) creado por un renderizador de senderos, lo que indica que las matemáticas son incorrectas ya que la excentricidad no está entre 0 y 1.

También he notado que el AP y el PE son incorrectos ya que el PE se muestra como -AP. Sé que mis matemáticas están mal en alguna parte, pero no tengo idea de dónde.

Si alguien mucho más inteligente que yo es capaz de resolverlo/arreglar esto, sería genial si también pudiera indicarme cómo usar esto para trazar una trayectoria futura en tiempo real.

void JustKeplerThings(Vector3 curPos, Vector3 shipCurVel, GameObject orbiting) {        
    //Gravitational Constant
    float G = 6.67408f * Mathf.Pow(10, -11);

    //Mass of ship
    float m1 = thisShip.mass;

    //Mass of the larger body
    float m2 = -orbiting.GetComponent<PointEffector2D>().forceMagnitude;
    float M = m2;

    //Standard Gravitational Parameter
    float µ = G * (m1 + m2);
    //float µ = G * M;

    //Relative Position Vector
    Vector3 r = FindRelativePosition(transform, orbitee.transform.position);

    //Relative Velocity Vector
    Vector3 v = shipCurVel - new Vector3(orbitee.GetComponent<Rigidbody2D>().velocity.x, orbitee.GetComponent<Rigidbody2D>().velocity.y, 0);

    //Specific Angular Momentum
    Vector3 h = Vector3.Cross(r, v);

    //Eccentricity Vector
    Vector3 evec = (Vector3.Cross(v, h) / µ) - (r / Vector3.Magnitude(r));

    //Eccentricity
    float e = Vector3.Magnitude(evec);

    //Mechanical Energy
    float E = (Mathf.Pow(Vector3.Magnitude(v), 2) / 2) - (µ / Vector3.Magnitude(r));

    //Semi-Major Axis
    float a = µ / (2 * E);
    if (Mathf.Abs(e - 1.0f) <= Mathf.Epsilon)
        a = Mathf.Infinity;

    //Semi-Latus Rectum Check
    float p = a * (1 - Mathf.Pow(e, 2));
    if (Mathf.Abs(e - 1.0f) <= Mathf.Epsilon)
        p = Mathf.Pow(Vector3.Magnitude(h), 2) / µ;

    //Inclination (2D)
    float i = 0;

    //Longitude of Ascending Node (2D)
    float Ω = 0;

    //Argument of Periapsis
    float ω = Mathf.Atan2(evec.y, evec.x);
    //ω = (2 * Mathf.PI) - ω;

    //True Anomaly
    float t = Mathf.Acos((Vector3.Dot(evec, r)) / (e * Vector3.Magnitude(r)));
    if (Vector3.Dot(r, v) < 0)
        t = (2 * Mathf.PI) - t;

    //Apoapsis
    float ap = a * (1 + e);

    //Periapsis
    float pe = a * (1 - e);

    //Orbital Period
    float T = (2 * Mathf.PI) * Mathf.Sqrt(Mathf.Pow(a, 3) / µ);
}

-

public static Vector3 FindRelativePosition(Transform origin, Vector3 position) {
    Vector3 distance = position - origin.position;
    Vector3 relativePosition = Vector3.zero;
    relativePosition.x = Vector3.Dot(distance, origin.right.normalized);
    relativePosition.y = Vector3.Dot(distance, origin.up.normalized);
    relativePosition.z = Vector3.Dot(distance, origin.forward.normalized);
    return relativePosition;
}