How to calculate earth's precession using JPL Horizons System?

Question

I'm new to JPL Horizons and understanding how to use the tool efficiently.

I'm looking to extract earth's precession data and understand where we are in the current cycle. I understand we progress 1 degree every 72 years, but I'm curious to know which degree point we are currently at.

Thank you to anyone who can aid!

Answer 1

There are a few components that make up precession, all of which are combined into a precession matrix, which is then applied to the position vector of the object in question.

Below I have provided code for the currently accepted Precession algorithm. It is in JavaScript, but should be easily ported to any language you're comfortable with.

My suggestion would be to start with some arbitrary vector (e.g. $ \vec{v}=[1,0,0] $ ). Generate two matrixes $P_1$, $P_2$. Compute the Matrix/Vector multiplications:

$$ \vec{v_1}=P_1\vec{v} \\ \vec{v_2}=P_2\vec{v} $$

Then Compute The Angle Between the Two Vectors. If the angle is lesser/greater than what you're looking for, regenerate $ P_2 $ for a time further/closer to $ P_1 $ until you've found the result to the accuracy you desire.

Since the Earth's precession period is about 25,772 years, a good starting point to find one degree change is $ p1 = 2451545.5 $ (Jan 1, 2000) and $ p2 = \frac{25,772years}{360}+p1 = 2477693.34 $.

function getPrecessionMatrix(jd_tdb){
    //Fukushima-Williams IAU 2006
    const t=(jd_tdb-2451545.5)/36525.0;

    const gamma = (-0.052928 + 10.556378*t + 0.4932044*t*t - 0.00031238*t*t*t - 0.000002788*t*t*t*t + 0.0000000260*t*t*t*t*t) /60/60*Math.PI/180;
    const phi = (+84381.412819 - 46.811016*t + 0.0511268*t*t + 0.00053289*t*t*t - 0.000000440*t*t*t*t - 0.0000000176*t*t*t*t*t) /60/60*Math.PI/180;
    const psi = (-0.041775 + 5038.481484*t + 1.5584175*t*t - 0.00018522*t*t*t - 0.000026452*t*t*t*t - 0.0000000148*t*t*t*t*t) /60/60*Math.PI/180;
    const eps = (+84381.406 - 46.836769*t - 0.0001831*t*t + 0.00200340*t*t*t - 0.000000576*t*t*t*t - 0.0000000434*t*t*t*t*t) /60/60*Math.PI/180;

    const a=Vec.getZRotationMatrix(gamma);
    const b=Vec.dot(Vec.getXRotationMatrix(phi),a);
    const c=Vec.dot(Vec.getZRotationMatrix(-psi),b);
    const d=Vec.dot(Vec.getXRotationMatrix(-eps),c);

    return d;
}

The input is a Julian Day, you can use an on-line Julian Day Converter. The code uses a Vector/Matrix library, which you should be able to find an existing one in your language of choice.

Answer 2

I don't think you can extract that data through the Horizons API. The software uses it, but not in a way which lets you separate only the precession from everything else. For example, observation coordinates of right ascension (RA) and declination (DEC) as seen from Earth are available in two forms: either "astrometric" or "apparent". Quoting the Horizons manual, in the section titled "Definition of Observer Table Quantities",

1. Astrometric RA & DEC: Adjusted for light-time aberration only. With respect to the reference plane and equinox of the chosen system (ICRF or FK4/B1950).
2. Apparent RA & DEC: with respect to the true-equator and Earth equinox of-date coordinate system (EOP-corrected IAU76/80 precession and nutation of the spin-pole) and adjusted to model light-time, the gravitational deflection of light, and stellar aberration, with an optional (approximate) correction for atmospheric yellow-light refraction.

To obtain just the precession and nothing else, I would use the Standards of Fundamental Astronomy (SOFA), in particular their earth attitude "cookbook" at https://www.iausofa.org/sofa_pn_c.pdf . That will give you not only the older (1976 precession and 1980 nutation) IAU models Horizon uses, but also the most recent (2006 precession and 2000 nutation) IAU models, a couple of others in between, and more possible choices about how to represent differences among coordinate systems.

Be careful with SOFA for operational use, since the newest stuff is slightly more accurate but much more expensive to compute, while the older ways may be good enough, and perhaps more useful depending on your purpose. If you need milli-arc-second precision, or are only going to use it once a day, then sure, go for 2006 and 2000A. If instead you need to run it many times in inner loops, you might be better off not taking quite so much time to polish just that one number.