Declination angle

Question

Currently I'm studying spherical geometry. I have come across the declination angle $\delta$ of the sun and there is an approximation to calculate $\delta$: $$\delta \approx -23.45° \cdot \cos\left(\frac{360}{365} \cdot (d+10)\right)$$ (where d is the day of the year with Jan 1 as d = 1, source here).
The idea behind this function is simple: You take a sinusodial function with amplitude $23.45°$.

Is it possible to derive a more accurate formula based on geometrical reasoning and where can I read about this? (In a German book I found the mystical hint that there exists a function that provides (almost) exact values - but without a source).

If not, can anyone recommend me a book that explains the algorithms for calculating $\delta$ (e.g. the PSA algorithm)?

Answer 1

The general idea behind how to tackle this problem is simple, and is stated in detail in a previous answer of mine here. There, a fixed geocentric system of coordinates is used to derive the declination and right ascension angles for a circular trajectory of the Earth around the Sun. Now, if we take into account the fact that this trajectory is actually better approximated by an ellipse, define the vector starting at the Sun and ending at the Earth

$$\mathbf{r}(t)=R(t)(\cos\phi(t), \sin\phi(t),0)$$

with $r(0)=r_{\max}, \theta(0)=0, \delta(0)=-\alpha\approx-23.45^\circ$. These initial conditions place the Earth at time zero at the winter solstice. Following the procedure outlined in the answer above and making appropriate adjustments in the coordinate system one can readily show

$$\sin\delta(t)=-\sin \alpha\cos\phi(t)$$

This formula is exact, but it needs to be supplemented by a dynamical equation coming from Newton's laws that helps us calculate the angle swept as a function of time. Given that the distance of the Earth from the Sun is given as a function of the angle by $$R(\phi)=\frac{2r_{\min}r_{\max}}{r_{\max}+r_{\min}-(r_{\max}-r_{\min})\cos\phi}$$

one can readily derive from Newton's second law that

$$\frac{d\phi}{dt}=\frac{J}{4mr_{\min}^2r^2_{\max}}\left(r_{\max}+r_{\min}-(r_{\max}-r_{\min})\cos\phi\right)^2$$

where $J^2=Gm^2M\frac{2r_\min r_\max}{r_{\min}+r_{\max}}$ is the conserved angular momentum. The ODE is separable and can in principle be solved exactly. After separation of variables, note that the integral on the LHS admits an elementary antiderivative

$$\int \frac{dx}{(a+b\cos x)^2}=\frac{2a}{(b^2-a^2)^{3/2}}\text{arctanh}\left(\sqrt{\frac{a - b}{a + b}}\tan (x/2)\right)-\frac{b \sin x}{(a^2 - b^2)(a + b \cos x)}+C$$

and hence we can write the solution explicitly

$$\frac{2a}{(b^2-a^2)^{3/2}}\text{arctanh}\left(\sqrt{\frac{a - b}{a + b}}\tan (\phi/2)\right)-\frac{b \sin \phi}{(a^2 - b^2)(a + b \cos \phi)}=t~~,~~ 0\leq\phi\leq \pi\\ \frac{2a}{(b^2-a^2)^{3/2}}\text{arctanh}\left(\sqrt{\frac{a - b}{a + b}}\tan (\phi/2)\right)-\frac{b \sin \phi}{(a^2 - b^2)(a + b \cos \phi)}=t-T~~,~~\pi<\phi< 2\pi$$

for appropriate $a,b$ and $T$ the period of the trajectory but in my opinion finding the angle as a function of time is best done by numerical integration, since the solution above is implicit and will have to be solved numerically anyway.