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

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;
public double[] eccentricity = new double;
public double[] inclination = new double;
public double[] meanLongitude = new double;
public double[] longitudeOfPerihelion = new double;
public double[] longitudeOfTheAscendingNode = new double;


So, for example, in the case of Earth, semimajorAxis would be 1.00000261, and semimajorAxis 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  + semimajorAxis  * tCenturiesFromJ2000;
double e = eccentricity  + eccentricity  * tCenturiesFromJ2000;
double i = inclination  + inclination  * tCenturiesFromJ2000;
double L = meanLongitude  + meanLongitude  * tCenturiesFromJ2000;

Debug.Log (L);

double p = longitudeOfPerihelion  + longitudeOfPerihelion  * tCenturiesFromJ2000;
double W = longitudeOfTheAscendingNode  + longitudeOfTheAscendingNode  * 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?

## 1 Answer

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.

• 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  + meanLongitude  * tCenturiesFromJ2000; I have to put the modulus function between two elements. Where exactly? – vistaero May 11 '17 at 0:13
• @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. – user May 11 '17 at 13:28
• 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? – vistaero May 12 '17 at 23:05
• @vistaero -- No! Use L = Math.Remainder(Math.IEEERemainder(Ldot*T,360) + L0, 360). This is basic math. – David Hammen May 13 '17 at 0:49
• 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$. – 2012rcampion May 14 '17 at 0:33