How do I calculate Earth's axial precession in an idealized Sun-Earth system?
I am experimenting with an idealized two-body system based on the Sun-Earth system. I want to calculate how much Earth's axis precesses and nutates given the simplified setup.
In my system, the orbit of the Earth about the Sun is a perfect circle. The sun is centered on the origin of the reference frame. At time $t = 0$, the Earth is located on the positive $x$-axis. 394
Earth's axis of rotation is tilted 22.5° from vertical about the line of the equinoxes, which at $t = 0$ is parallel to and aligned with the $x$-axis. There are no other bodies in the system; the only thing gravitationally acting on the Earth is the Sun. The Sun is a perfect sphere, while the Earth is a spheroid with an equatorial radius of 6371 km and a polar radius of 6356 km.
Wikipedia gives this equation for the torque caused by a celestial body's gravity acting on the Earth:
$$\vec{T} = \frac{3GM}{r^3} (C − A) \sin(δ) \cos(δ) \begin{pmatrix} \sin(α)\ −\cos(α)\ 0\ \end{pmatrix}$$
Where
* $GM$ is the standard gravitational parameter, the product of the gravitational constant $G$ and the mass $M$ of the perturbing body;
* $r$ is the distance between the center of the Earth and the center of the perturbing body;
* $C$ is the moment of inertia around Earth's axis of rotation;
* $A$ is the moment of inertia around any equatorial diameter of Earth;
* $(C − A)$ is the moment of inertia of Earth's equatorial bulge (C > A);
* $δ$ is the declination of the perturbing body (positive for north of the equator, negative for south of the equator); and
* $α$ is the right ascension of the perturbing body (east from vernal equinox)
I have values for all these variables, making the equation for the sun:
$$\vec{T_s}(t) = 6.916 × 10^{44} \frac{kg \: m^2}{s^2} \sin(δ_{s}(t)) \cos(δ_{s}(t)) \begin{pmatrix} \sin(α_{s}(t))\ −\cos(α_{s}(t))\ 0\ \end{pmatrix}$$
Where - $δ_{s}(t) = \arcsin(-\sin(t) \sin(π/8))$ - $α_{s}(t) = \arctan(\tan(t) \cos(π/8))$, adjusted to cover the range [0, 2π)
Wikipedia says that the $y$ component of the vector averages to zero and can be neglected, leaving only the $x$ component. In my simplified system, the average value of the $x$ component of the vector should be
$$T_x = \left(\frac{3}{2}\right) \left(\frac{GM}{r^3}\right) (C-A) \sin(π/8) \cos(π/8)$$
Which works out to
$$T_x = 1.546 × 10^{22} \frac{kg \: m^2}{s^2} \sin(π/8) \cos(π/8)$$
Then Wikipedia says precession is
$$\frac{dψ}{dt} = \left(\frac{3}{2}\right) \left(\frac{GM}{r^3}\right) \left(\frac{C - A}{C}\right) \left(\frac{\cos(π/8)}{ω}\right)$$
Which works out to
$$\frac{dψ}{dt} = 2.46623 × 10^{-12} rad^{-1} s^{-1}$$
Which somehow converts to arcseconds per year.
That gives an average value, but is there a way to get the precession over time more exactly, rather than as an average? 371