Solar Geometry and Insolation

This page provides the mathematical formulations used in Insolation.jl to calculate solar geometry and insolation. The equations are from the textbook draft "Physics of Earth's Climate", Chapter 3.

Overview

Calculating solar position in the sky and insolation requires a sequence of astronomical computations that determine:

1. The planet's position in its orbit (mean and true anomaly)
2. The latitude of the subsolar point (declination)
3. The longitude relative to solar noon (hour angle)
4. The geometric relationships between the sun and a surface location (zenith and azimuth angles)
5. The planet-star distance for radiation calculations

Mean Anomaly

The mean anomaly $M$ represents the angular position of a planet in its orbit, assuming uniform circular motion. For an elliptical orbit, this is a convenient starting point for more accurate calculations.

The mean anomaly at current time $t$ is

\[M = \frac{2\pi (t - t_0)}{Y_a} + M_0,\]

where we have:

• The time at the epoch (J2000), $t_0$, defined as January 1, 2000 at 12:00 Terrestrial Time (TT); the package default is the corresponding UTC instant, 11:58:55.816
• The mean anomaly at the epoch, $M_0$ [radians]
• The length of the anomalistic year, $Y_a$ (period from perihelion to perihelion) [seconds]

The mean anomaly increases linearly with time at a rate of $2\pi/Y_a$ radians per second.

True Anomaly

The true anomaly $A$ is the actual angular position of the planet in its elliptical orbit, measured from perihelion (the point of closest approach to the star). Unlike the mean anomaly, the true anomaly accounts for the planet's varying orbital speed due to Kepler's laws.

The true anomaly is computed from the mean anomaly using a series expansion,

\[A = M + \left( 2e - \frac{1}{4}e^{3} \right) \sin(M) + \frac{5}{4} e^2 \sin(2M) + \frac{13}{12} e^3 \sin(3M) + \mathcal{O}(e^4),\]

where $e$ is the orbital eccentricity. This series approximation is accurate to order $e^3$ and is sufficient for Earth's relatively circular orbit ($e \approx 0.017$).

Solar Longitude

The solar longitude $L_s$ (also called ecliptic longitude or true longitude) is the ecliptic longitude of the Sun, measured from the vernal equinox. It combines the orbital phase (true anomaly $A$) with the longitude of perihelion $\varpi$: the Sun's ecliptic longitude at perihelion, measured from the vernal equinox (≈ 283° for present-day Earth):

\[L_s = A + \varpi.\]

For Earth, $\varpi$ varies slowly due to precession (period $\sim 21{,}000$ years).

Declination

The solar declination $\delta$ is the angle between the sun's rays and the equatorial plane — equivalently, the latitude of the subsolar point, where the sun is directly overhead at solar noon. The declination varies between $\pm\gamma$ (obliquity) over the course of a year.

The sine of the declination angle is

\[\sin \delta = \sin \gamma \sin L_s,\]

where we have:

• The orbital obliquity $\gamma$ (axial tilt) [radians]
• The solar longitude $L_s$ [radians]

For Earth's current obliquity of $\gamma \approx 23.44°$, the declination ranges from $-23.44°$ (northern winter solstice) to $+23.44°$ (northern summer solstice), passing through $0°$ at the equinoxes. Its extremes are the latitudes of the Tropics of Capricorn and Cancer. The package default, obliq_epoch, is the mean obliquity of the ecliptic at J2000 — the IAU 2006 value, $23.43927944° = 84381.406''$.

Equation of Time

The equation of time corrects for the difference between apparent solar time (when the sun is actually at its highest point) and mean solar time (uniform clock time). This discrepancy arises from two effects:

1. The elliptical orbit causes varying orbital speed (eccentricity effect)
2. The tilted axis projects the sun's motion onto the equatorial plane (obliquity effect)

The equation-of-time hour angle correction [radians] is computed as the difference between the mean longitude $L = M + \varpi$ and the right ascension $\alpha$ of the true Sun,

\[\Delta \eta = L - \alpha, \qquad \alpha = \operatorname{atan2}\left(\cos\gamma \, \sin L_s,\; \cos L_s\right),\]

where $L_s = A + \varpi$ is the solar longitude. The right ascension follows from the exact projection of the ecliptic longitude onto the equatorial plane, so the obliquity contribution is valid for any tilt $\gamma$ (including $\gamma > 90°$), while the eccentricity contribution enters through the true anomaly $A$.

For small eccentricity and obliquity this reduces to the familiar perturbative form

\[\Delta \eta \approx -2 e \sin(M) + \tan^2(\gamma/2) \sin(2M+2\varpi),\]

where the first term accounts for orbital eccentricity and the second for axial tilt. The correction can be converted to a time offset through $\Delta t = \Delta\eta T_d/(2\pi)$, where $T_d$ is the length of the solar day. For present-day Earth it runs from about $-14$ minutes in February to $+16$ minutes in November, which is why sundials and clocks disagree.

Hour Angle

The hour angle $\eta$ measures the angular distance of the sun from the local meridian (north-south line). It quantifies how far past (or before) solar noon we are at a given location:

• Local solar noon: $\eta = 0$
• In the afternoon: $\eta > 0$
• In the morning: $\eta < 0$

The hour angle is calculated from the time of day, with equation of time correction, and adjusted for longitude:

\[\eta = \left( \eta_\text{uncorrected} + \Delta\eta \right) + \lambda,\]

where we have:

• The uncorrected hour angle at the prime meridian (0° longitude), $\eta_\text{uncorrected} = 2\pi t_\text{day}$ [radians]
• The fractional part of the Julian date, $t_\text{day}$ [dimensionless]
• The equation of time hour angle correction, $\Delta\eta$ [radians]
• The longitude $\lambda$ [radians]

The Julian day begins at noon, so $t_\text{day}$ is 0 at 12:00 UTC and 0.5 at 00:00 UTC; this is what makes $\eta = 0$ at mean solar noon on the prime meridian, with no further offset. The factor $2\pi$ converts fractional day to angle.

All terms are taken modulo $2\pi$ for proper angle wrapping, so the returned hour angle lies in $[0, 2\pi)$: morning hours appear as $\eta$ slightly below $2\pi$ rather than as negative values. Since $\cos\eta$ is what enters the zenith angle, the wrapping is immaterial to the insolation.

Zenith Angle

The zenith angle $\theta$ is the angle between the sun's rays and the vertical direction (zenith) at a location. It determines how directly sunlight strikes a surface:

• Sun directly overhead (zenith): $\theta = 0°$
• Sun at the horizon (sunrise/sunset): $\theta = 90°$
• Sun below the horizon (night): $\theta > 90°$

The cosine of the zenith angle is

\[\cos \theta = \cos \phi \cos \delta \cos \eta + \sin \phi \sin \delta,\]

where we have:

• The latitude, $\phi$ [radians], positive northward
• The solar declination angle, $\delta$ [radians]
• The hour angle, $\eta$ [radians]

This is a fundamental equation in solar geometry. The incident solar radiation is proportional to $\cos \theta$, which is why high solar zenith angles (low sun) produce less heating than low zenith angles (overhead sun).

Sunrise/Sunset Angle

The sunrise/sunset hour angle $\eta_d$ is the hour angle at which the sun crosses the horizon. It determines the length of day and night:

• Day length: the sun is above the horizon for hour angles between $-\eta_d$ and $+\eta_d$, a fraction $\eta_d/\pi$ of the day, i.e., for a time $2\eta_d \, T_d/(2\pi)$ with $T_d$ the length of the solar day
• Polar day when $\tan \phi \tan \delta \ge 1$ (the sun never sets, $\eta_d = \pi$)
• Polar night when $\tan \phi \tan \delta \le -1$ (the sun never rises, $\eta_d = 0$)

The sunrise/sunset hour angle is given by

\[\cos \eta_d = - \tan \phi \tan \delta,\]

which follows from setting $\theta = 90°$ (sun at horizon) in the zenith angle formula. Sunrise occurs at $-\eta_d$ and sunset at $+\eta_d$, symmetrically about solar noon. In the code, the two polar cases need no branching: clamping the right-hand side to $[-1, 1]$ before taking the arccosine returns $\eta_d = \pi$ for polar day and $\eta_d = 0$ for polar night automatically.

Note that both polar cases require $|\phi| \ge 90° - |\delta|$, so they can only occur poleward of the polar circles at $\pm(90° - \gamma)$ — for present-day Earth, poleward of 66.6°.

Diurnally Averaged Insolation

The diurnally averaged (daily-mean) insolation requires averaging $\cos \theta$ over a full 24-hour period. This is helpful for conceptual models and simpler climate models that do not resolve the full diurnal cycle.

Since insolation is proportional to $\cos \theta$, we need $\overline{\cos \theta}$, the time-averaged cosine of the zenith angle. This can be computed analytically from the sunrise/sunset hour angle:

\[\overline{\cos \theta} = \frac{1}{\pi} \left( \eta_d \sin \phi \sin \delta + \cos \phi \cos \delta \sin \eta_d \right).\]

Physical interpretation:

• The $1/\pi$ prefactor combines the $1/(2\pi)$ average over the full 24-hour day with a factor of 2 from the symmetric sunrise-to-sunset integration limits ($\pm\eta_d$)
• When $\eta_d = 0$ (polar night), $\overline{\cos \theta} = 0$ (no insolation)
• When $\eta_d = \pi$ (polar day), $\overline{\cos \theta} = \sin \phi \sin \delta$: the sun circles the horizon without setting, and only the term that survives a full rotation remains
• At the equator ($\phi = 0$) during equinox ($\delta = 0$), we have $\eta_d = \pi/2$ and the formula reduces to $(1/\pi) \cos \phi \cos \delta \sin \eta_d = 1/\pi \approx 0.318$

This expression is what makes summer insolation maximal at the pole rather than at the subsolar latitude near the solstices: the polar-day term $\eta_d \sin\phi \sin\delta$ grows with latitude fast enough to beat the falling solar elevation. For a planet on a nearly circular orbit, this happens whenever the obliquity exceeds about $20.7°$, which Earth's $23.4°$ does.

Azimuth Angle

The azimuth angle $\zeta$ specifies the horizontal direction to the sun. It is essential for tracking solar panels, understanding shading, and computing radiation on tilted surfaces. Note that this package uses a non-standard convention, measuring $\zeta$ from due East and increasing counter-clockwise (see below), rather than the more common convention of measuring clockwise from due North.

The azimuth angle is

\[\zeta = \frac{3\pi}{2} - \operatorname{atan2}\left( \cos\delta \sin \eta,\; \cos\delta \cos \eta \sin \phi - \sin\delta \cos \phi \right),\]

where the two-argument arctangent ($\operatorname{atan2}$) resolves the correct quadrant. The numerator and denominator are written with a common factor of $\cos\delta \ge 0$ (rather than the more familiar $\sin\eta / (\cos\eta\sin\phi - \tan\delta\cos\phi)$) to avoid the $\tan\delta$ singularity at high declinations.

Convention in this package:

• Sun due East: $\zeta = 0$ (or $2\pi$)
• Sun due North: $\zeta = \pi/2$
• Sun due West: $\zeta = \pi$
• Sun due South: $\zeta = 3\pi/2$

The azimuth increases counter-clockwise when viewed from above. At local solar noon ($\eta = 0$) the sun is exactly due south ($\zeta = 3\pi/2$) or due north ($\zeta = \pi/2$), depending on which side of the subsolar latitude the observer is on: due south when $\phi > \delta$, due north when $\phi < \delta$. This is not the same as the hemisphere — in the tropics in summer the subsolar point can lie poleward of the observer, and the noon sun then stands to the north in the Northern Hemisphere.

Planet-Star Distance

The planet-star distance $d$ varies throughout the year due to the elliptical orbit. This variation affects the solar flux received at the top of atmosphere through the inverse square law ($S \propto 1/d^2$).

The distance is calculated from the equation for the orbital ellipse:

\[d = \frac{1-e^2}{1+e\cos A} d_0,\]

where we have:

• The orbital eccentricity $e$ (0 for circular, $0<e<1$ for elliptical) [unitless]
• The true anomaly $A$ [radians]
• The semi-major axis $d_0$, the mean of the perihelion and aphelion distances, at which the total solar irradiance $S_0$ is defined

Since the insolation depends only on the ratio $d/d_0$ (the inverse-square law gives $S = S_0 (d_0/d)^2$ with $S_0$ the irradiance at the mean distance), Insolation.jl works with the dimensionless distance $d/d_0 = (1-e^2)/(1+e\cos A)$ and never needs the absolute orbit size. The returned distance is therefore expressed in units of the semi-major axis.

For Earth:

• Eccentricity $e \approx 0.0167$ (current value, varies over millennia)
• Perihelion (closest): $d/d_0 \approx 0.983$ (early January)
• Aphelion (farthest): $d/d_0 \approx 1.017$ (early July)
• The $\pm 1.7\%$ variation in distance causes a $\pm 3.4\%$ variation in solar flux
• The semi-major axis is $d_0 \approx 1.496 \times 10^{11}$ m $= 1$ AU; multiply $d/d_0$ by this to recover a physical distance

Earth is closest to the sun during Northern Hemisphere winter, but the obliquity effect dominates over the distance effect in determining seasons: a $\pm 3.4\%$ swing in flux is small next to the factor-of-several swing in $\overline{\cos\theta}$ that the tilt produces at middle and high latitudes.

There is also an exact cancellation worth knowing about. The flux varies as $d^{-2}$, and by Kepler's second law so does the orbital angular velocity, so the planet lingers in a given arc of its orbit for a time proportional to $d^{2}$. The two effects cancel, and the energy received over any fixed interval of solar longitude — an astronomical season, or the whole year — is independent of where in the orbit that interval falls (Herschel's law). Changing the longitude of perihelion therefore redistributes insolation within the seasonal cycle without changing any seasonal total; see Milankovitch Cycles for a numerical demonstration.