Celestial Programming : Precession IAU 1976

This is an implementation of IAU 1976 Precession from Astronomical Algorithms eq 21.2

$$ \begin{align*} t &=\frac{JD-JD0}{36525} \\ \\ T &= \frac{JD0 - 2451545.0}{36525} \\ \\ \zeta &= (2306.2181 + 1.39656T - 0.000139T^2)t + (0.30188 - 0.000344T)t^2 + 0.017998t^3 \\ z &= (2306.2181 + 1.39656T - 0.000139T^2)t + (1.09468 + 0.000066T)t^2 + 0.018203t^3 \\ \theta &= (2004.3109 -0.85330T - 0.000217T^2)t -(0.42665 + 0.000217T)t^2 - 0.041833t^3 \\ \\ \tan(\alpha - z) &= \frac{\cos(\delta_0) \sin(\alpha_0 + z)}{\cos(\theta) \cos(\delta_0) \cos(\alpha_0 + z) - \sin(\theta) \sin(\delta_0)} \\ \\ \sin(\delta) &= \sin(\theta) \cos(\delta_0) \cos(\alpha_0 + z) + \cos(\theta) \sin(\delta_0) \\ \end{align*} $$ JD - JD to Precess to.
JD0 - Starting epoch (2451545 for J2000)
\(\alpha_0\) - Starting RA
\(\delta_0\) - Starting Dec
\(\alpha\) - Ending RA
\(\delta\) - Ending RA


Alternate matrix method

$$ \vec r=R_3(-z)R_2(-\theta)R_3(-\zeta)\vec r_0 $$


Simplified, if starting from J2000

$$ \begin{align*} \zeta &=2306.2181"t + 0.30188"t^2 + 0.017998"t^3 \\ z &=2306.2181"t + 1.09468"t^2 + 0.018203"t^3 \\ \theta &=2004.3109"t - 0.42665"t^2 - 0.041833"t^3 \\ \end{align*} $$