# 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[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?

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 [0] + meanLongitude [1] * tCenturiesFromJ2000; I have to put the modulus function between two elements. Where exactly? – vistaero May 11 '17 at 0:13
• @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