Why does NASA Horizons say that Earth inclination is not 0 at J2000 epoch?

Question

I'm using NASA Horizons system to get Earth orbital elements with following settings:

Ephemeris Type : ELEMENTS
Target Body :    Earth-Moon Barycenter [EMB] [3]
Center :         Solar System Barycenter (SSB) [500@0]
Time Span :      Start=2000-01-01, Stop=2001-01-01, Step=12 h
Table Settings : output units=KM-S

This is what it says in the output:

Output units    : KM-S, deg, Julian Day Number (Tp)
Output type     : GEOMETRIC osculating elements
Output format   : 10
Reference frame : ICRF/J2000.0
Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch

And here are the orbital elements for the J2000 epoch:

2451545.000000000 = A.D. 2000-Jan-01 12:00:00.0000 TDB 
 EC= 1.487370652262564E-02 QR= 1.466664481240442E+08 IN= 1.196622093669550E-02
 OM= 9.348586355904454E+00 W = 6.534678781572920E+01 Tp=  2451519.451250376645
 N = 1.149765722413008E-05 MA= 2.538006615432194E+01 TA= 2.612305949633834E+01
 A = 1.488808583174952E+08 AD= 1.510952685109462E+08 PR= 3.131072643603165E+07

So as you can see the inclination is about 0.012° which is not 0 (by the way, same goes both for Earth and Earth-Moon Barycenter). I thought inclination of the Earth orbit in J2000 reference frame must be 0 at J2000 epoch by definition. Could someone clarify this please?

P.S. Now that I think about it and look at this picture my question becomes more like "Why is the inclination not 23.4°?"

Answer 1

Here's the heart of your problem:

Center :         Solar System Barycenter (SSB) [500@0]

The ecliptic is defined in terms of the Earth's mean (or average) orbit about the Sun at the epoch rather than about the solar system barycenter. Had you instead chosen the Sun as the center you would have found an inclination of 0.000266 degrees (about one arc second) at noon TDB on 1 Jan 20000. Aside: The Earth-Moon barycenter had an inclination of about 0.0001 degrees at that time per the Horizons web site.

There are a number of reasons for that small one arc second discrepancy. One is that JPL's HORIZONS system currently uses the DE431 JPL ephemerides to compute the positions and velocities of solar system bodies. This ephemeris was released in 2014, decades after the J2000.0 frames were defined.

Another reason is the way JPL generates its ephemerides. JPL uses a batch least squares technique to iteratively refine the initial states for an integrated model of the solar system in an attempt to minimize the weighted sum squared error between observations and what this integrated model predicts would be observed. JPL does not release those epoch states. Instead, it releases Chebyshev polynomial coefficients that describe the positions of solar system bodies in something close to a minimum absolute error sense. This is one of the beauties of Chebyshev polynomials; they come very close to the ideal of a minimum absolute error predictor. However, the generator only looks at position; it does not look at velocity. Both position and velocity are needed to compute the osculating orbital elements.

Yet another reason is that those Chebyshev polynomial coefficients are expressed in terms of the ICRF frame rather than an ecliptic frame. A transformation is needed to convert those ICRF positions and velocities to an ecliptic frame. Small numerical errors (e.g., see https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html) will creep in while performing those transformations. Other numerical errors will creep in while calculating the inclination. These numerical errors alone explain a good portion of that one arc second discrepancy.

Finally, the transformation from the J2000 Earth mean equator / Earth mean equinox frame (which is very close to the ICRF frame) to the J2000 ecliptic / Earth mean equinox frame was computed using a much older version of JPL's solar system model, and it was done using the integrated model rather than the Chebyshev polynomial coefficients.

Answer 2

I'll venture an answer.

First, if you look at the actual data output from the ephemeris, you can see that around the year 2000 the Earth+Moon center of mass (barycenter) is moving in a plane parallel to the ecliptic within a few hundred kilometers. Both the Earth and the Sun slowly drift up and down because the origin of the coordinate system is actually the center of mass of the solar system, and the other planets, Jupiter and Uranus in particular, have to be taken into account.

I never really use those generated orbital elements. They are an approximation - a fit - to the motion of the body at that particular time. You can propagate the motion of the body for a short time using them, but they are not meant to be used to describe long-term motion.

The solar system is a bunch of objects all interacting with each other. The sun is not (really) in the center, and the planets do not (exactly) even orbit around the Sun.

So, for example if you look at the same orbital elements around the Sun instead of around the entire solar system barycenter, you'll see that the inclination drops to 0.0001 degrees. If at all possible, stick to the cartesian (x, y, z) data to avoid having to think about these issues.

above: data from JPL Horizons database. These are the z-coordiantes only - perpendicular to the ecliptic. Units are kilometers. Solid smooth line is the Earth-Moon barycenter. Solid wiggly line is the Earth Geocenter, and dashed line is the sun. TOP: raw z-coordinates, BOTTOM: z-coordinates minus z-coordinates of the Sun.

Away from the year 2000, the Earth is moving up and down once per year relative to the sun. It flattens out around 2000, when the Earth's orbit coincides with the ecliptic by definition.

Here is the python script I used to read the data and make the plot. It's not pretty but it does the job.

import numpy as np
import matplotlib.pyplot as plt

fnames = ['J2000 Earth Barycenter horizons_results.txt',
          'J2000 Earth Geocenter horizons_results.txt',
          'J2000 Sun horizons_results.txt']

names  = ['Earth_Barycenter', 'Earth_Geocenter', 'Sun']

JDs, positions = [], []

for fname in fnames:

    with open(fname, 'r') as infile:

        lines = infile.read().splitlines()

        iSOE = [i for i, line in enumerate(lines) if "$$SOE" in line][0]
        iEOE = [i for i, line in enumerate(lines) if "$$EOE" in line][0]

        print iSOE, iEOE, lines[iSOE], lines[iEOE]

        lines = [line.split(',') for line in lines[iSOE+1:iEOE]]
        JD  = np.array([float(line[0]) for line in lines])
        pos = np.array([[float(item) for item in line[2:5]] for line in lines])

        JDs.append(JD)
        positions.append(pos)

JD = JDs[0]  # assume they are identical 

class Obj(object):
    def __init__(self, name):
        self.name = name

things = []

for name, pos, JD in zip(names, positions, JDs):

    thing     = Obj(name)
    thing.pos = pos
    thing.JD  = JD

    things.append(thing)

earth_bary, earth_geo, sun = things

years = 2000. + (JD - JD[1826])/365.2564

if 1 == 1:
    fig = plt.figure()

    ax1 = plt.subplot(2,1,1)
    ax1.plot(years, earth_bary.pos.T[2])
    ax1.plot(years, earth_geo.pos.T[2])
    ax1.plot(years, sun.pos.T[2], '--')
    ax1.set_xlim(1995, 2005)
    ax1.get_xaxis().get_major_formatter().set_useOffset(False)

    ax2 = plt.subplot(2,1,2)
    ax2.plot(years, earth_bary.pos.T[2]-sun.pos.T[2])
    ax2.plot(years, earth_geo.pos.T[2]-sun.pos.T[2])
    ax2.plot(years, sun.pos.T[2]-sun.pos.T[2], '--')
    ax2.set_xlim(1995, 2005)
    ax2.get_xaxis().get_major_formatter().set_useOffset(False)

    plt.show()