3
$\begingroup$

Sorry bad english.

This should be a comment to an 2012rcampion's answer to this question: Determining orbital position at a future point in time, but as I'm new and I don't have enough reputation, I have to publish a new question.

I'm using the data of the PDF linked in his answer, under "Computing Position from Orbital Elements". I don't have enough reputation so I can't post it here... (no more than two links allowed, I suppose to avoid spam, but hey, the first link is a link to this webpage! and the second one I was trying to post is a link to NASA webpage! this limits are a total pain to new users!)

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 semimajorAxis1 would be 0.00000562.

Now, the method (updated after first answer to do some fixes):

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:

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 others simulators.

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

Thanks in advance.

$\endgroup$
3
$\begingroup$

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.

$\endgroup$
  • $\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 '17 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$ – a CVn May 11 '17 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 '17 at 23:05
  • 1
    $\begingroup$ @vistaero -- No! Use L = Math.Remainder(Math.IEEERemainder(Ldot*T,360) + L0, 360). This is basic math. $\endgroup$ – David Hammen May 13 '17 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$ – 2012rcampion May 14 '17 at 0:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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