Skip to content

Instantly share code, notes, and snippets.

@mayakraft
Last active February 12, 2023 22:16
Show Gist options
  • Save mayakraft/7578514 to your computer and use it in GitHub Desktop.
Save mayakraft/7578514 to your computer and use it in GitHub Desktop.
Keplerian Elements for Approximate Positions of the Major Planets
// http://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf
typedef enum{
Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto
} Planet;
// a e I L long.peri. long.node.
// AU rad deg deg deg deg
static double elements[] = {0.38709927, 0.20563593, 7.00497902, 252.25032350, 77.45779628, 48.33076593, //mercury
0.72333566, 0.00677672, 3.39467605, 181.97909950, 131.60246718, 76.67984255, //venus
1.00000261, 0.01671123, -0.00001531, 100.46457166, 102.93768193, 0.0, //earth moon barycenter
1.52371034, 0.09339410, 1.84969142, -4.55343205, -23.94362959, 49.55953891, //mars
5.20288700, 0.04838624, 1.30439695, 34.39644051, 14.72847983, 100.47390909, //jupiter
9.53667594, 0.05386179, 2.48599187, 49.95424423, 92.59887831, 113.66242448, //saturn
19.18916464, 0.04725744, 0.77263783, 313.23810451, 170.95427630, 74.01692503, //uranus
30.06992276, 0.00859048, 1.77004347, -55.12002969, 44.96476227, 131.78422574, //neptune
39.48211675, 0.24882730, 17.14001206, 238.92903833, 224.06891629, 110.30393684 };//pluto
// AU/Cy rad/Cy deg/Cy deg/Cy deg/Cy deg/Cy
static double rates[] = {0.00000037, 0.00001906, -0.00594749, 149472.67411175, 0.16047689, -0.12534081, //mercury
0.00000390, -0.00004107, -0.00078890, 58517.81538729, 0.00268329, -0.27769418, //venus
0.00000562, -0.00004392, -0.01294668, 35999.37244981, 0.32327364, 0.0, //earth moon barycenter
0.00001847, 0.00007882, -0.00813131, 19140.30268499, 0.44441088, -0.29257343, //mars
-0.00011607, -0.00013253, -0.00183714, 3034.74612775, 0.21252668, 0.20469106, //jupiter
-0.00125060, -0.00050991, 0.00193609, 1222.49362201, -0.41897216, -0.28867794, //saturn
-0.00196176, -0.00004397, -0.00242939, 428.48202785, 0.40805281, 0.04240589, //uranus
0.00026291, 0.00005105, 0.00035372, 218.45945325, -0.32241464, -0.00508664, //neptune
-0.00031596, 0.00005170, 0.00004818, 145.20780515, -0.04062942, -0.01183482 };//pluto
// location of planet in the J2000 ecliptic plane, with the X-axis aligned toward the equinox
-(double*)calculateLocationOfPlanet:(Planet)planet AtTime:(float)time{
static double ecliptic[3];
double *planet0 = &elements[6*planet];
double *per_century = &rates[6*planet];
// step 1
// compute the value of each of that planet's six elements
double a = planet0[0] + per_century[0]*time; // (au) semi_major_axis
double e = planet0[1] + per_century[1]*time; // ( ) eccentricity
double I = planet0[2] + per_century[2]*time; // (°) inclination
double L = planet0[3] + per_century[3]*time; // (°) mean_longitude
double ϖ = planet0[4] + per_century[4]*time; // (°) longitude_of_periapsis
double Ω = planet0[5] + per_century[5]*time; // (°) longitude_of_the_ascending_node
// step 2
// compute the argument of perihelion, ω, and the mean anomaly, M
double ω = ϖ - Ω;
double M = L - ϖ;
// step 3a
// modulus the mean anomaly so that -180° ≤ M ≤ +180°
while(M > 180) M-=360; // in degrees
// step 3b
// obtain the eccentric anomaly, E, from the solution of Kepler's equation
// M = E - e*sinE
// where e* = 180/πe = 57.29578e
double E = M + (e*180./π) * sin(M*π/180.); // E0
for(int i = 0; i < 5; i++){ // iterate for precision, 10^(-6) degrees is sufficient
E = [self KeplersEquation:E M:M e:e];
}
// step 4
// compute the planet's heliocentric coordinates in its orbital plane, r', with the x'-axis aligned from the focus to the perihelion
ω = ω * π/180.;
E = E * π/180.;
I = I * π/180.;
Ω = Ω * π/180.;
double x0 = a*(cos(E)-e);
double y0 = a*sqrt(1-e*e)*sin(E);
// step 5
// compute the coordinates in the J2000 ecliptic plane, with the x-axis aligned toward the equinox:
ecliptic[x] = ( cos(ω)*cos(Ω) - sin(ω)*sin(Ω)*cos(I) )*x0 + ( -sin(ω)*cos(Ω) - cos(ω)*sin(Ω)*cos(I) )*y0;
ecliptic[y] = ( cos(ω)*sin(Ω) + sin(ω)*cos(Ω)*cos(I) )*x0 + ( -sin(ω)*sin(Ω) + cos(ω)*cos(Ω)*cos(I) )*y0;
ecliptic[z] = ( sin(ω)*sin(I) )*x0 + ( cos(ω)*sin(I) )*y0;
return ecliptic;
}
-(double)KeplersEquation:(double)E M:(double)M e:(double)e{
double ΔM = M - ( E - (e*180./π) * sin(E*π/180.) );
double ΔE = ΔM / (1 - e*cos(E*π/180.));
return E + ΔE;
}
@alexgenaud
Copy link

This was a helpful interpretation of JPL's napkin quality summary (I understood 180/(PI * e) whereas you correctly had e*180/PI). Thanks for posting.

I believe there is a sign mistake on line 53. Perhaps it should be instead:

double E = M - (e*180./π) * sin(M*π/180.);  // E0

image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment