An integral for the Earth's insolation Consider the function $$ [-\pi/2,\pi/2] \ni \theta \mapsto s_\beta(\theta) = \int_0^{2\pi} \sqrt{ 1 - \left(\cos \theta \sin \beta \cos \gamma \- \sin \theta \cos \beta \right)^2} \, d \gamma $$ for $\beta \in [0, \pi/2]$. This function is proportional to the average insolation (amount of solar energy) that reaches the top of the atmosphere of a planet with obliquity (axial tilt) $\beta$, at latitude $\theta$. So $\theta = 0$ at the equator and $\theta = \pi/2$ at the north pole. If $\beta = 0$, the planet's axis of rotation is perpendicular to its ecliptic plane, and if $\beta = \pi/2$, this axis is parallel to the plane. Numerical evidence suggests that $s_\beta$ is an even function of $\theta$. I am looking for a proof.

The integrand is invariant under $\theta \to -\theta$ and $\gamma \to \gamma + \pi$.