function compute_celestial_positions contains the line
- 22438.42 * np.cos(l_prime-f_omega_d)
where f_omega_d = 2 * (F + omega - D)
But I think the line should be
- 22438.42 * np.cos(l_prime+f_omega_d)
because IERS table 5.2b says
i b_{s,j})i b{c,j})_i l l' F D Om
5 -17.40 22438.42 0 1 2 -2 2