Celestial Programming : Low Precision Moon Phase - Simpson

$$ \begin{align*} X = &[ 383.0 \sin(8399.685t + 5.381) \\ &+ 31.5 \sin(70.990t + 6.169) \\ &+ 10.6 \sin(16728.377t + 1.453) \\ &+ 6.2 \sin(1185.622t + 0.481) \\ &+ 3.2 \sin(7143.070t + 5.017) \\ &+ 2.3 \sin(15613.745t + 0.857) \\ &+ 0.8 \sin(8467.263t + 1.010)] \times 10^6 m \\ \\ Y = &[351.0 \sin(8399.687t + 3.811) \\ &+ 28.9 \sin(70.997t + 4.596) \\ &+ 13.7 \sin(8433.466t + 4.766) \\ &+ 9.7 \sin(16728.380t + 6.165) \\ &+ 5.7 \sin(1185.667t + 5.164) \\ &+ 2.9 \sin(7143.058t + 0.300) \\ &+ 2.1 \sin(15 613.755t + 5.565)] \times 10^6 m \\ \\ Z = &[153.2 \sin(8399.672t + 3.807) \\ &+ 31.5 \sin(8433.464t + 1.629) \\ &+ 12.5 \sin(70.996t + 4.595) \\ &+ 4.2 \sin(16728.364t + 6.162) \\ &+ 2.5 \sin(1185.645t + 5.167) \\ &+ 3.0 \sin(104.881t + 2.555) \\ &+ 1.8 \sin(8399.116t + 6.248)] \times 10^6 m \\ \end{align*} $$ t is Julian centuries since J2000 ((jd-2451545)/36525).
Arguments to the sin function are already in radians.