Computing Sunlight: The NOAA Solar Algorithm

There is a seductive simplicity to the premise: give me a latitude, a longitude, and a date, and I will tell you how long the sun will be above the horizon. No network call. No API key. Just math. The DayLength app does exactly this, computing 365 sunrise and sunset times in a single synchronous pass on the CPU. The algorithm fits in roughly 50 lines of Swift.

Why Astronomy Uses a Fractional Year

The NOAA surface radiation budget algorithm begins by converting a calendar date into a single floating-point angle: the fractional year γ (gamma). This collapses the calendar into a unit the trigonometric functions can operate on directly.

// Pseudocode — illustrates the pattern
func fractionalYear(dayNumber: Int, totalDays: Int) -> Double {
    return (2 * π / Double(totalDays)) * Double(dayNumber - 1)
}

For January 1st, γ is 0. For December 31st in a non-leap year, γ approaches 2π. The division by the total number of days means the calculation self-corrects for leap years — a detail worth noting because getting it wrong by one day accumulates into a visible error at year boundaries.

Everything downstream — the equation of time, the solar declination — is a function of this single angle.

Two Things You Need to Know About the Sun

Given γ, the algorithm derives two quantities that together determine when sunrise and sunset occur on any given day.

The first is the equation of time (eqTime): the difference, in minutes, between solar noon (when the sun is highest in the sky) and clock noon (12:00:00 UTC offset by longitude). This difference exists because Earth’s orbit is elliptical, not circular, and because the axial tilt means the sun’s apparent east-west speed varies through the year. The equation of time oscillates between roughly −14 and +16 minutes over the course of a year.

The second is solar declination (δ): the angle between the sun’s rays and Earth’s equatorial plane. This is what causes seasons. At the summer solstice, δ reaches approximately +23.45°. At the winter solstice, −23.45°. At the equinoxes, δ is 0° and every location on Earth gets close to exactly 12 hours of daylight.

Both quantities are computed using Fourier series — sums of sine and cosine terms at the fundamental frequency γ and its harmonics:

// Pseudocode — illustrates the Fourier approximation pattern
func solarParameters(gamma: Double) -> (eqTime: Double, declination: Double) {
    let eqTime = scale * (
        constant
        + a1 * cos(gamma)   - b1 * sin(gamma)
        + a2 * cos(2*gamma) - b2 * sin(2*gamma)
    )
    let declination = (
        c0
        - c1 * cos(gamma)  + c2 * sin(gamma)
        - c3 * cos(2*gamma) + c4 * sin(2*gamma)
        - c5 * cos(3*gamma) + c6 * sin(3*gamma)
    )
    return (eqTime, declination)
}

The coefficients are not derived on the fly — they are empirical constants published by NOAA that encode decades of astronomical observation. The declination uses three harmonics (six terms); the equation of time uses two (four terms). This is not laziness. A three-harmonic approximation is accurate to within a fraction of a degree across the entire year. More terms would add complexity without improving the accuracy enough to matter for a consumer app where the target precision is ±2 minutes against verified reference data.

The Hour Angle and the Polar Cases

With declination in hand, the algorithm computes the solar hour angle at sunrise (HA). This is the angle the Earth needs to rotate for the sun to reach the horizon. The formula comes from spherical trigonometry:

cos(HA) = -tan(latitude) × tan(declination)

Here is where the geometry becomes directly legible. If you are at the equator (latitude = 0°), tan(0) = 0, so cos(HA) = 0, giving HA = 90° regardless of the season. The sun rises and sets exactly 90° from solar noon — always 12 hours of daylight. Latitude matters; the equator is immune to seasons.

flowchart TD
    A["cos(HA) = −tan(lat) × tan(decl)"] --> B{Value range}
    B -->|"cos(HA) > 1"| C["HA = NaN\nPolar night — sun never rises"]
    B -->|"cos(HA) < −1"| D["HA = ∞\nMidnight sun — sun never sets"]
    B -->|"−1 ≤ cos(HA) ≤ 1"| E["HA = acos(cos(HA))\nNormal day"]
    E --> F["sunrise at −HA\nsunset at +HA"]

The polar cases are the most interesting design decision in the whole calculator. When cos(HA) > 1, the sun geometrically cannot rise — polar night. When cos(HA) < −1, the sun cannot set — midnight sun. The implementation encodes these as IEEE 754 sentinel values: .nan for polar night, .infinity for midnight sun. The calling site checks isNaN and isInfinite before converting to a Date.

This is a good use of floating-point semantics. The alternatives — throwing an error, returning an Optional<Double>, adding a separate enum — all require the caller to handle an extra case. The sentinel approach keeps calculateHourAngle returning a single Double, and the mathematical reason for the sentinel (the argument to acos is out of domain) maps cleanly onto the IEEE 754 special values that represent that condition.

Converting an Angle into a UTC Timestamp

Once HA is known, sunrise and sunset times follow from a single formula. The sun rises when the Earth has rotated HA degrees west of solar noon, and sets when it has rotated HA degrees east:

// Pseudocode — illustrates the UTC time derivation
func utcMinutes(longitude: Double, eqTime: Double, hourAngle: Double) -> Double {
    // 720 = 12 hours in minutes (solar noon at UTC midnight offset)
    // 4 minutes per degree of longitude
    return 720.0 - (4.0 * longitude) - eqTime + hourAngle * (4.0 * 180.0 / π)
}

The 4.0 * longitude term converts longitude to minutes: Earth rotates 360° in 1440 minutes, so 1° = 4 minutes. The 4.0 * 180.0 / π converts radians to minutes via the same rate. The result is minutes after UTC midnight on the given date.

This result is then added to the UTC midnight of that day to produce a Date. Importantly, timeUTC can be slightly negative (sunrise before UTC midnight, e.g., in New Zealand in summer) or greater than 1440 (sunset after UTC midnight). Swift’s addingTimeInterval handles both gracefully — the date arithmetic wraps correctly without any explicit clamping.

The Architectural Bet: UTC All the Way Down

The calculator has one job: produce Date values. It knows nothing about the user’s timezone, the target location’s timezone, or daylight saving time. Every output is an absolute UTC timestamp.

This is the right call, and it is non-trivial to get right. Day length — sunset minus sunrise — is a duration, and durations are timezone-agnostic. A day that is 14 hours 23 minutes long is 14 hours 23 minutes long in Tokyo and in Toronto. The dayLength computation never touches the timezone at all; it is pure arithmetic on TimeInterval.

---

Deep Dive: Why Three Harmonics?

The declination equation uses six sine/cosine terms (three harmonic pairs). The equation of time uses four (two pairs). Why not one? Why not ten?

Earth’s axial tilt drives the fundamental annual cycle — the sin(γ) and cos(γ) terms. The first harmonic captures the eccentricity of the orbit and the interaction between tilt and eccentricity, which produces the “figure eight” analemma shape. The second and third harmonics in the declination account for higher-order perturbations in Earth’s motion that are small but measurable.

At one harmonic, the error in computed declination reaches ±0.25° near the solstices. At three harmonics, the error is below ±0.02°. A 0.25° error in declination translates to roughly 1–2 minutes of error in sunrise/sunset time at mid-latitudes — at the edge of the ±2-minute acceptance threshold the tests enforce. Adding a fourth harmonic would improve accuracy by perhaps another 0.005°, which maps to seconds. The three-harmonic choice is an explicit precision budget: good enough for a consumer app, minimal enough to stay readable.

---

Apply This

IEEE 754 sentinels for domain boundaries. When a function returns a floating-point value that has meaningful out-of-domain states (the argument to acos is outside [−1, 1]), encoding those states as .nan and .infinity keeps the return type simple and makes the caller’s intent explicit: if angle.isNaN { return polarNight }. Reserve this for cases where the sentinel values are mathematically motivated — not as a general substitute for Optional or Result.

Fourier approximation as a lookup table replacement. The NOAA coefficients are essentially a compressed lookup table: instead of storing 365 × 2 precomputed values (730 floats), five to seven constants reconstruct any value in the table to within measurement precision. This pattern applies wherever a periodic or quasi-periodic signal needs to be evaluated at arbitrary points: audio synthesis, animation curves, sensor calibration. The cost is upfront derivation of the coefficients; the benefit is a calculation that fits in a register and runs in nanoseconds.