Keplerian Elements for Approximate Positions of the Major Planets

Kepler.m

// 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;
}

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