Is this C# code to obtain the coordinates of the planets correct?

Question

This should be a comment to an 2012rcampion's answer to this question: Determining orbital position at a future point in time

I'm using the linked data of that answer, under "Computing Position from Orbital Elements".

As each parameter has two values, I have declared the variables as arrays of two values:

public double[] semimajorAxis = new double[2];
public double[] eccentricity = new double[2];
public double[] inclination = new double[2];
public double[] meanLongitude = new double[2];
public double[] longitudeOfPerihelion = new double[2]; 
public double[] longitudeOfTheAscendingNode = new double[2];

So, for example, in the case of Earth, semimajorAxis[0] would be 1.00000261, and semimajorAxis[1] would be 0.00000562.

Now, the method is:

void ComputePosition(){
    double tMillisFromJ2000 = DateTime.Now.ToUniversalTime().Subtract(new DateTime(2000, 1, 1, 12, 0, 0, DateTimeKind.Utc)).TotalMilliseconds;
    double tCenturiesFromJ2000 = tMillisFromJ2000 / (1000 * 60 * 60 * 24 * 365.25 * 100);
    double a = semimajorAxis [0] + semimajorAxis [1] * tCenturiesFromJ2000;
    double e = eccentricity [0] + eccentricity [1] * tCenturiesFromJ2000;
    double i = inclination [0] + inclination [1] * tCenturiesFromJ2000;
    double L = meanLongitude [0] + meanLongitude [1] * tCenturiesFromJ2000;

    Debug.Log (L);

    double p = longitudeOfPerihelion [0] + longitudeOfPerihelion [1] * tCenturiesFromJ2000;
    double W = longitudeOfTheAscendingNode [0] + longitudeOfTheAscendingNode [1] * tCenturiesFromJ2000;

    double M = L - p;
    double w = p - W;

    double E = M;
    while (true) {
        double dE = (E - e * Mathd.Sin (E) - M) / (1 - e * Mathd.Cos (E));
        E -= dE;
        if (Mathd.Abs (dE) < 1e-6)
            break;
    }

    double P = a * (Mathd.Cos (E) - e);
    double Q = a * Mathd.Sin (E) * Math.Sqrt (1 - Mathd.Pow (e, 2d));

    // rotate by argument of periapsis
    double x = Mathd.Cos(w) * P - Mathd.Sin(w) * Q;
    double y = Mathd.Sin (w) * P + Mathd.Cos (w) * Q;
    // rotate by inclination
    double z = Mathd.Sin(i) * x;
    x = Mathd.Cos (i * x);
    // rotate by longitude of ascending node
    double xTemp = x;
    x = Mathd.Cos (W) * xTemp - Mathd.Sin (W) * y;
    y = Mathd.Sin (W) * xTemp + Mathd.Cos (W) * y;

    position = new Vector3d (x * 149597870700d, y * 149597870700d, z * 149597870700d);
}

I've checked the code and orbital parameters three times but it's giving me wrong positions for at least some planets. Neptune, 29 AU above the Sun, totally normal:

I can't see any sense in these positions. They don't look any similar to what we can see in other simulators.

Is there something wrong or missing in the code? should I check the parameters again?

Answer 1

The document to which you wished to link is Keplerian Elements for Approximate Positions of the Major Planets. You did not read the instructions! Step 1 says "Compute the value of each of that planet's six elements: $a= a_0 + \dot a T$, etc." You didn't do the et cetera part of that step. Be very careful when computing the mean longitude $L*$. You need to bring the result into the range $(-180^\circ,180^\circ]$ using the modulus function (% in the C family of languages).

Another place you didn't read the instructions: The angles in that document are in degrees. I would diverge from that document and convert to radians after computing M (modulo 360 degrees). Then you don't have to do the silly step of computing eccentricity in degrees.

---

Based on requests, and a closed follow-up question, how to compute the modulus of a real number with respect to another real number is language-dependent. With python, it's easy: Just use the % operator. You'll get almost what you expect. For example,

>>> (10.123456789*math.pi % math.pi)/math.pi
0.12345678900000093

That little bit at the end? That's floating point arithmetic for you.

The % operator doesn't work this way in the C family of languages. In C, you'll need to use fmod or remainder, prefaced by std:: in C++. The language in question is C#. Here, you'll need to use math.IEEERemainder.