Problemas con el cálculo de las posiciones X, Y, Z de los planetas en un momento determinado en Unity3D (C#)

Estoy tratando de crear un pequeño programa en Unity3D que posiciona los planetas en el espacio 3D en un momento dado. He seguido este documento: http://stjarnhimlen.se/comp/ppcomp.html .

Usando los datos de la sección 4 y el código de las secciones 6 y 7, debería tener el algoritmo correcto, ¿verdad? Omití los cambios en los atributos orbitales (por ejemplo, + 4.70935E-5 * d) porque solo necesito que se muevan en una órbita establecida. Realmente no necesito que cambien con el tiempo. Por ahora, solo estoy aumentando la anomalía media con el tiempo sin tener en cuenta el período.

El problema es que los planetas no parecen tener el desplazamiento correcto con respecto a la eclíptica, y su movimiento es realmente extraño: aceleran, luego disminuyen la velocidad y retroceden en sus caminos. Aquí hay un video de su movimiento: https://www.youtube.com/watch?v=zHOBChFpan0

Parece que la excentricidad es correcta en sus elipses y las distancias al sol también son correctas, pero el resto es simplemente incorrecto. No soy tan bueno en matemáticas o astronomía como para entender completamente lo que está pasando en este algoritmo, pero espero que alguien aquí pueda detectar el error.

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;
}

}