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

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!

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.

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.

• The SOFA code, while extremely good, is perhaps a bit too good ("'good enough' is the enemy of the best") for most uses by spacecraft. The key driver of the SOFA Earth orientation code as of late has been milliarcsecond astronomy. Getting Earth orientation to that high a level of accuracy takes a lot of CPU time. Those easy-to-use functions are big CPU hogs. The old 1976 precession / 1980 nutation model are a lot faster, but of course not as precise. Jul 1, 2022 at 20:26
• @DavidHammen I agree! incorporating change... Jul 1, 2022 at 20:33
• @RyanC Thank you for the help and showing me SOFA! After fighting with it for a bit, I realized I'm more rusty with C/C++ than I thought I was. I will continue fighting with it, but I'm just looking for the date we cross of the next 30 deg interval. My current estimate is between 2100 and 2400, but I was hoping to get a more precise estimate verify with a numerical tool. Jul 1, 2022 at 23:51
• I guess you could get a reasonable estimate of the axial precession from Horizons using those two RAs... if you neglect nutation. Meuss gives a simple equationfor approximating nutation, but it's only intended for short-term use, not equinox precession cycles. :) Jul 2, 2022 at 4:28
• These guys have easily accessible code for obliquity, and two IAU models of nutation: neoprogrammics.com/obliquity_of_the_ecliptic Jul 2, 2022 at 4:29