A stone sundial in a courtyard. Choose a latitude, a date and a time and watch the real shadow of the gnomon move, while the sky shows the Sun's arc for that day, the three key arcs of the year and the figure-8 analemma. Then change Earth's axial tilt or orbit shape to see why the seasons exist.
For each axial tilt the sweep scans every minute of the solstice day, finding the lowest latitude where the Sun never drops below the horizon (-0.833 deg), and compares it with the closed form: latitude = 90 deg - tilt - 0.833 deg.
The science
Why seasons exist. Earth's spin axis is tilted by ε = 23.44° to its orbit and keeps pointing the same way. The Sun's ecliptic longitude λ advances through the year, and the Sun's declination is sin δ = sin ε sin λ, swinging between ±ε at the solstices. The distance from the Sun is a minor effect: Earth is closest in early January.
Orbit. Mean anomaly M = 2π(t - tp)/365.2422; Kepler's equation E - e sin E = M is solved by Newton iteration; true anomaly ν follows from E and λ = ν + ϖ with the longitude of perihelion ϖ = 282.9 deg; distance r/a = 1 - e cos E.
Sun in the sky. With latitude φ and hour angle H (15 deg per hour from solar noon): sin h = sinφ sinδ + cosφ cosδ cos H, azimuth A = atan2(sin H, cos H sinφ - tanδ cosφ) + 180°. Noon altitude is 90° - |φ - δ|.
Day length. Sunrise/sunset at H0 where cos H0 = (sin h0 - sinφ sinδ)/(cosφ cosδ), h0 = -0.833 deg (refraction + solar radius). Day length = 2H0/15 deg per hour. If the right side is above 1 there is polar night, below -1 midnight sun.
Equation of time. The true Sun runs fast and slow because of the eccentric orbit and the tilt: EoT = (M + ϖ) - α, where α is the Sun's right ascension (tan α = cosε tanλ). A sundial reads mean time + EoT: about -14 min in mid-February and +16 min in early November. The analemma is the Sun's position at the same clock time through the year: height = δ, width = EoT.
Gnomon. The shadow of a vertical pointer of height h has length h / tan(altitude) and points opposite the Sun's azimuth. The red curve on the dial is the locus of the shadow tip during the day (a hyperbola on a flat plane).
Sunlight at the top of the atmosphere. Daily mean Q = (S0/π)(r/a)-2(H0 sinφ sinδ + cosφ cosδ sin H0), S0 = 1361 W/m² (typical value; h0 = 0 here). The global annual mean is S0/(4√(1-e²)).
Honest simplifications. No nutation, precession or leap years; declination and EoT are evaluated at local noon of the chosen date; the 3D Sun position uses the date's noon declination for the whole day; the sky is a stylised dome; clouds and the atmosphere's effect on sunlight at the ground are not modelled. The Spencer (1971) series used as a cross-check is itself an approximation good to about a minute. The "what if" presets (tilt 0, 82 deg, e = 0.15) are hypothetical planets, not forecasts.