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)

    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:

Neptune at 29 AU above the Sun?!

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?


1 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

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.

  • $\begingroup$ Thanks! Now the code makes more sense (I've updated it in the question), but I don't get how do I use the modulus function. This is the line of code at which I compute the variable L: double L = meanLongitude [0] + meanLongitude [1] * tCenturiesFromJ2000; I have to put the modulus function between two elements. Where exactly? $\endgroup$
    – vistaero
    May 11, 2017 at 0:13
  • $\begingroup$ @vistaero Either this answer does adequately answer your question, and should then be accepted; or it doesn't, and shouldn't be accepted. Don't edit your question after the fact to ask follow-up questions, especially combined with accepting an answer (as this indicates to the community that you feel your question has been adequately answered). It's okay to edit your question if the original question wasn't clear enough for a useful answer to be given, but if you require more help, you should ask a new question instead. $\endgroup$
    – user
    May 11, 2017 at 13:28
  • $\begingroup$ Please David Hammen, telling me just to use the modulus function isn't helpful. I can't find anything on the internet about what has to do the remainder of a division with bringing a value into a range. Can you add some example to your answer? $\endgroup$
    – vistaero
    May 12, 2017 at 23:05
  • 1
    $\begingroup$ @vistaero -- No! Use L = Math.Remainder(Math.IEEERemainder(Ldot*T,360) + L0, 360). This is basic math. $\endgroup$ May 13, 2017 at 0:49
  • 1
    $\begingroup$ I don't think you need to shift the mean longitude at all; Kepler's equation works for any value of $M$, not only values in the range $-\pi<M<+\pi$. $\endgroup$ May 14, 2017 at 0:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.