Milanković’s orbital insolation (1920)

Milutin Milanković computed the sunlight each latitude receives through the year and how it changes as the tilt and orbit of the Earth vary. He proposed that ice ages begin when high northern summers receive little sunlight. The daily mean at the top of the atmosphere depends on latitude, the Sun’s declination and the Earth-Sun distance. Whatever the season, the global mean is a quarter of the solar constant.

import matplotlib.pyplot as plt
import numpy as np

from earthsciencekit.climate import daily_insolation

obliquity = np.deg2rad(23.44)
day = np.arange(365)
declination = np.arcsin(np.sin(obliquity) * np.sin(2 * np.pi * (day - 80) / 365.25))
latitude = np.deg2rad(np.linspace(-90, 90, 181))
flux = daily_insolation(1361, latitude[:, None], declination[None, :], 1)
global_mean = np.trapezoid(flux * np.cos(latitude)[:, None], latitude, axis=0) / 2
np.testing.assert_allclose(global_mean, 1361 / 4, rtol=1e-3)
assert np.isclose(daily_insolation(1361, np.pi / 2, obliquity, 1), 1361 * np.sin(obliquity))

tilt = np.deg2rad(np.linspace(22.1, 24.5, 50))
summer_65n = daily_insolation(1361, np.deg2rad(65), tilt, 1)
assert np.all(np.diff(summer_65n) > 0)

fig, (left, right) = plt.subplots(1, 2, figsize=(10, 4))
image = left.contourf(day, np.rad2deg(latitude), flux, levels=20, cmap="magma")
fig.colorbar(image, ax=left, label="daily-mean insolation [W/m²]")
left.set(xlabel="day of year", ylabel="latitude [degree]", title="Sunlight through the year")
right.plot(np.rad2deg(tilt), summer_65n)
right.set(xlabel="obliquity [degree]", ylabel="insolation [W/m²]", title="65°N at the summer solstice")
fig.tight_layout()

plt.show()