Probleme beim Berechnen der X-, Y- und Z-Positionen von Planeten zu einem bestimmten Zeitpunkt in Unity3D (C#)

Ich versuche, ein kleines Programm in Unity3D zu erstellen, das die Planeten zu einem bestimmten Zeitpunkt im 3D-Raum positioniert. Ich habe dieses Dokument befolgt: http://stjarnhimlen.se/comp/ppcomp.html .

Mit den Daten aus Abschnitt 4 und dem Code aus Abschnitt 6 und 7 sollte ich den richtigen Algorithmus haben, oder? Ich habe die Änderungen in den Umlaufbahnattributen (z. B. + 4,70935E-5 * d) übersprungen, weil ich sie nur in einer festgelegten Umlaufbahn bewegen muss. Ich brauche nicht wirklich, dass sie sich im Laufe der Zeit ändern. Im Moment erhöhe ich nur die mittlere Anomalie im Laufe der Zeit, ohne den Zeitraum zu berücksichtigen.

Das Problem ist, dass die Planeten nicht den richtigen Versatz zur Ekliptik zu bekommen scheinen und ihre Bewegung wirklich seltsam ist: Sie werden schneller, dann langsamer und kehren ihre Bahnen um. Hier ist ein Video ihrer Bewegung: https://www.youtube.com/watch?v=zHOBChFpan0

Es sieht so aus, als ob die Exzentrizität auf ihren Ellipsen korrekt ist und die Entfernungen zur Sonne auch korrekt sind, aber der Rest ist einfach falsch. Ich bin nicht so gut in Mathematik oder Astronomie, um vollständig zu verstehen, was in diesem Algorithmus vor sich geht, aber ich hoffe, dass jemand hier den Fehler entdecken könnte.

using UnityEngine;
using System;

public class Orbit : MonoBehaviour 
{
public float longitudeOfTheAscendingNode;
public float inclination;
public float argumentOfPerihelion;
public float semimajorAxis;
public float eccentricity;

public float meanAnomalySpeed = 0.1f;

float M = 0; //mean anomaly

void Update()
{
    SetPosition();
}

void SetPosition()
{
    float N = longitudeOfTheAscendingNode;
    float i = inclination;
    float w = argumentOfPerihelion;
    float a = semimajorAxis;
    float e = eccentricity;

    float E = M + e * (180 / Mathf.PI) * Mathf.Sin(M) * (1 + e * Mathf.Cos(M));

    float xv = a * (Mathf.Cos(E) - e);
    float yv = a * (Mathf.Sqrt(1 - e * e) * Mathf.Sin(E));

    float v = Mathf.Atan2(yv, xv);
    float r = Mathf.Sqrt(xv * xv + yv * yv);

    float xx = r * (Mathf.Cos(N) * Mathf.Cos(v + w) - Mathf.Sin(N) * Mathf.Sin(v + w) * Mathf.Cos(i));
    float yy = r * (Mathf.Sin(N) * Mathf.Cos(v + w) + Mathf.Cos(N) * Mathf.Sin(v + w) * Mathf.Cos(i));
    float zz = r * (Mathf.Sin(v + w) * Mathf.Sin(i));

    transform.position = new Vector3(xx, yy, zz);

    M += meanAnomalySpeed * Time.deltaTime;
    if (M > 360)
        M -= 360;
}

}