Celestial Programming : Nutation Matrix

Produces a Nutation Matrix given the nutation quantities \(\Delta\psi, \Delta\epsilon, \) and \(\epsilon_0 \). The values \(\Delta\psi, \Delta\epsilon, \) come from the nutation algorithms, and the value \(\epsilon_0 \) is computed from the functions at the bottom. Even though the matrix is formed the same regardless of the nutation algorithm used, care must still be taken to ensure compatable versions of precession and nutation algorithms are used. For example, IAU1976 precession and IAU1980 nutation are to be used together, while IAU2006 Precession can be used with IAU2000 or IAU2006 nutation algorithms. The obliquity algorithms given below should match the nutation algorithm.

Using rotation matrixes $$ \vec N=R_1(-\epsilon)R_3(-\Delta\psi)R_1(-\epsilon_0) $$ or explaning out the rotations into a single matrix:
$$ \vec N = \begin{bmatrix} \cos \Delta\psi & -\sin \Delta\psi \cos \epsilon_0 & -\sin\Delta\psi \sin \epsilon_0 \\ \sin \Delta\psi \cos \epsilon & \cos \Delta\psi \cos \epsilon \cos \epsilon_0 + \sin \epsilon \sin \epsilon_0 & \cos \Delta\psi \cos \epsilon \sin \epsilon_0 - \sin \epsilon \cos \epsilon_0 \\ \sin \Delta\psi \sin \epsilon & \cos \Delta\psi \sin \epsilon \cos \epsilon_0 - \cos \epsilon \sin \epsilon_0 & \cos \Delta\psi \sin \epsilon \sin \epsilon_0 + \cos \epsilon \cos \epsilon_0 \end{bmatrix} $$