2
$\begingroup$

I'm writing a program to compute the heliocentric ecliptic position of planets and asteroids. I want to output rectangular coordinates oriented according to the ecliptic: positive X in the direction of vernal equinox, positive Z toward the north ecliptic pole.

My solution roughly matches other publicly available Kepler simulations (just eyeballing), but I think it is slightly off. What is wrong and how do I verify it?

Here is my process for calculating position at a given JD:

  1. Pull in orbital elements from JPL SSD website. Specifically, the EM-Barycenter elements from this table.

  2. Compute period. For example, in the case of Earth:

METERS_IN_AU = 149597870700
SECONDS_IN_DAY = 86400

GM_SUN = 1.3271244004193938E+20              // m^3/s^2
SEMIMAJOR_AXIS = 1.0000026 * METERS_IN_AU    // meters

period = 2 * PI * sqrt(SEMIMAJOR_AXIS^3 / GM_SUN)
periodInDays = period / SECONDS_IN_DAY

GM is from NASA NAIF. The above formula returns periodInDays=365.2583228291575, as opposed to the commonly cited 365.2422 days in a year.

  1. Compute meanMotion = 2 * PI / periodInDays (radians per day)

  2. Compute mean anomaly at JD=2458563.415278 (spring equinox):

epoch = time of orbital elements (jd)
meanLon = mean longitude (radians)
lonPeri = longitude of perihelion (radians)

meanAnomalyAtEpoch = meanLon - lonPeri;
M = meanAnomalyAtEpoch + meanMotion * (2458563.415278 - epoch);

This can produce an unusually large angle. In this case, M=120.6878 radians!

  1. Estimate eccentric anomaly E using iterative approximation:
    let E = M;
    let lastdiff;
    do {
      const E1 = M + e * sin(E);
      lastdiff = abs(E1 - E);
      E = E1;
    } while (lastdiff > 1e-6);
  1. Compute true anomaly from eccentric anomaly E and eccentricity e:
v = 2 * atan(sqrt((1 + e) / (1 - e)) * tan(E / 2))
  1. Compute radius in AU:
r = a * (1 - e * e) / (1 + e * cos(v)); 
  1. Compute heliocentric coordinates:
o = longitude of ascending node (radians)
p = longitude of perihelion (radians)
i = inclination (radians)

X = r * (cos(o) * cos(v + p - o) - sin(o) * sin(v + p - o) * cos(i))
Y = r * (sin(o) * cos(v + p - o) + cos(o) * sin(v + p - o) * cos(i))
Z = r * (sin(v + p - o) * sin(i))

In the case of Earth with these elements, I get the following outputs for JD=2443568.0 (code here):

Period (days): 365.25832830801806 
Mean motion (rad/day): 0.01720203160400233 
Mean anomaly at epoch (rad): -0.043163916976385774 
Mean anomaly at JD 2458563.415278 (rad): 120.68783750518942 
Eccentric anomaly (rad): 120.70404042143021 
True anomaly (rad): 1.3397564871245002 
Radius (AU): 0.9959122843645666 
X coord (AU): -0.9958986393991808 
Y coord (AU): 0.005213270672998027 
Z coord (AU): -1.3930375794150342e-9 

Given that JD=2458563.415278 is the vernal equinox, I would perhaps expect Y be closer to zero.

$\endgroup$

Your Answer

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

Browse other questions tagged or ask your own question.