Sunrise/Sunset Algorithm Fix: Handle RA Wrap and UTC Boundary Events by Normalizing Local Mean Time T

Learn why common sunrise/sunset code fails at year boundaries and equinox RA wraps, and how to fix it: normalize local mean time T and stop force-ranging UT.

Sunrise/sunset calculations based on the well-known Stack Overflow snippet work surprisingly well for many locations, but they break in two places: around the turn of the year and on days where two events can occur in a single UTC day. The root symptoms are familiar: if you force-range UT to [0, 24), you miss events; if you don’t, the first or last day of the year goes wrong. Let’s unpack what actually fails and how to correct it without changing the overall algorithm.

Reproducing the issue

The following compact class performs the usual approximation for sunrise and sunset times. It mirrors the common implementation: UT is rounded and not force-ranged; T (local mean time) is left as-is. That is enough to reproduce the discrepancies described above.

import math
from datetime import datetime, timedelta, time, timezone

RAD = math.pi / 180.0

class SolarError(Exception):
    def __init__(self, msg):
        super(SolarError, self).__init__(msg)

class SolarClock:
    def __init__(self, lat_deg, lon_deg):
        self._lat = lat_deg
        self._lon = lon_deg
        self.lonHours = self._lon / 15.0

    def rise_time(self, on_date=datetime.now(), tz=timezone.utc):
        dt_offset = self.compute_UT(on_date, tz=tz, as_rise=True)
        if dt_offset is None:
            raise SolarError('The sun never rises on this location (on the specified date)')
        elif isinstance(dt_offset, str):
            return dt_offset
        else:
            return datetime.combine(on_date, time(tzinfo=tz)) + dt_offset

    def set_time(self, on_date=datetime.now(), tz=timezone.utc):
        dt_offset = self.compute_UT(on_date, tz=tz, as_rise=False)
        if dt_offset is None:
            raise SolarError('The sun never rises on this location (on the specified date)')
        elif isinstance(dt_offset, str):
            return dt_offset
        else:
            return datetime.combine(on_date, time(tzinfo=tz)) + dt_offset

    def compute_UT(self, on_date, tz, as_rise=True, zenith=90.8):
        if tz is None:
            tz = datetime.now().tzinfo

        day_index = on_date.timetuple().tm_yday

        if as_rise:
            approxT = day_index + ((6 - self.lonHours) / 24.0)
        else:
            approxT = day_index + ((18 - self.lonHours) / 24.0)

        meanAnom = (0.9856 * approxT) - 3.289

        trueLong = meanAnom + (1.916 * math.sin(RAD * meanAnom)) + (0.020 * math.sin(RAD * 2 * meanAnom)) + 282.634
        trueLong = self.wrap_range(trueLong, 360)

        sinDec = 0.39782 * math.sin(RAD * trueLong)
        cosDec = math.cos(math.asin(sinDec))

        cosH = (math.cos(RAD * zenith) - (sinDec * math.sin(RAD * self._lat))) / (cosDec * math.cos(RAD * self._lat))
        if cosH > 1:
            return 'constant_eclipse'
        if cosH < -1:
            return 'constant_sun'

        if as_rise:
            hourAngle = 360 - (1 / RAD) * math.acos(cosH)
        else:
            hourAngle = (1 / RAD) * math.acos(cosH)
        hourAngle /= 15.0

        rightAsc = (1 / RAD) * math.atan(0.91764 * math.tan(RAD * trueLong))
        rightAsc = self.wrap_range(rightAsc, 360)
        Lquad = (math.floor(trueLong / 90)) * 90
        RAquad = (math.floor(rightAsc / 90)) * 90
        rightAsc = (rightAsc + (Lquad - RAquad)) / 15.0

        localMean = hourAngle + rightAsc - (0.06571 * approxT) - 6.622

        utHours = localMean - self.lonHours
        if tz:
            utHours += tz.utcoffset(on_date).total_seconds() / 3600.0

        utHours = round(utHours, 2)

        return timedelta(hours=utHours)

    @staticmethod
    def wrap_range(val, maxv):
        if val < 0:
            return val + maxv
        elif val >= maxv:
            return val - maxv
        return val

Toggle range forcing on UT in this variant and you will observe exactly the trade-off reported earlier: forcing UT loses events within a day, while not forcing UT causes boundary-day errors.

What really goes wrong

The approximation computes a single right ascension and declination for a 24-hour window tied to the input date, changing only at 00:00 UTC. That simplification is fragile where sunrise/sunset shifts quickly, for example near the first or last events at high latitudes. More immediately, there is a discrete jump in the right ascension around the March equinox when it wraps from about 23.9h to 0h. Because T is computed as H + RA − 0.06571·t − 6.622, that RA wrap introduces a spurious −24h jump in T. If UT is then range-forced to [0, 24), you destroy information that legitimately says “this event belongs to the next day,” which is exactly how two events can occur within a UTC day when days are short and the time zone definition slices a local day differently. Conversely, leaving UT unforced exposes the false jump originating in T, which shows up at the year boundary and other sensitive dates.

The fix

Range-force T to [0, 24) and do not force-range UT. By normalizing the local mean time before converting to UT, you remove the artificial −24h step caused by RA wrapping, while preserving whether an event lands on the previous/next UTC day. Here is the corrected core with that one critical change.

    def compute_UT(self, on_date, tz, as_rise=True, zenith=90.8):
        if tz is None:
            tz = datetime.now().tzinfo

        day_index = on_date.timetuple().tm_yday
        approxT = day_index + ((6 - self.lonHours) / 24.0) if as_rise else day_index + ((18 - self.lonHours) / 24.0)

        meanAnom = (0.9856 * approxT) - 3.289
        trueLong = meanAnom + (1.916 * math.sin(RAD * meanAnom)) + (0.020 * math.sin(RAD * 2 * meanAnom)) + 282.634
        trueLong = self.wrap_range(trueLong, 360)

        sinDec = 0.39782 * math.sin(RAD * trueLong)
        cosDec = math.cos(math.asin(sinDec))
        cosH = (math.cos(RAD * zenith) - (sinDec * math.sin(RAD * self._lat))) / (cosDec * math.cos(RAD * self._lat))
        if cosH > 1:
            return 'constant_eclipse'
        if cosH < -1:
            return 'constant_sun'

        hourAngle = (360 - (1 / RAD) * math.acos(cosH)) if as_rise else ((1 / RAD) * math.acos(cosH))
        hourAngle /= 15.0

        rightAsc = (1 / RAD) * math.atan(0.91764 * math.tan(RAD * trueLong))
        rightAsc = self.wrap_range(rightAsc, 360)
        Lquad = (math.floor(trueLong / 90)) * 90
        RAquad = (math.floor(rightAsc / 90)) * 90
        rightAsc = (rightAsc + (Lquad - RAquad)) / 15.0

        localMean = hourAngle + rightAsc - (0.06571 * approxT) - 6.622
        localMean = self.wrap_range(localMean, 24)

        utHours = localMean - self.lonHours
        if tz:
            utHours += tz.utcoffset(on_date).total_seconds() / 3600.0

        utHours = round(utHours, 2)
        return timedelta(hours=utHours)

Once T is normalized and UT is not, you retain next-day or previous-day events and eliminate the artificial discontinuity.

Working with date ranges safely

When building a table over a long interval (such as a full year), compute events for one extra day before the start and one extra day after the end. Then store results keyed by the date of the result, not the input date, and discard anything outside your target range. This sidesteps the “event falls on the next day” situation and keeps the sequence clean. The following minimal driver shows the pattern.

from datetime import datetime, timedelta, timezone

sol = SolarClock(-30, -100)
utz = timezone(timedelta(hours=0))

rows = []
start0 = datetime(2024, 1, 1, tzinfo=timezone.utc)

# Pre-allocate by target date
ptr = start0
count = 0
while ptr.year == start0.year:
    rows.append([ptr.date().isoformat(), [], []])  # date, rises, sets
    ptr += timedelta(days=1)
    count += 1

# Iterate one day before and one day after the target interval
d = start0 - timedelta(days=1)
ix_r = 0
ix_s = 0
for _ in range(count + 2):
    r = sol.rise_time(d, utz)
    if isinstance(r, datetime) and r.year == start0.year:
        key = r.date().isoformat()
        while rows[ix_r][0] != key:
            ix_r += 1
        rows[ix_r][1].append(r.time().isoformat())

    s = sol.set_time(d, utz)
    if isinstance(s, datetime) and s.year == start0.year:
        key = s.date().isoformat()
        while rows[ix_s][0] != key:
            ix_s += 1
        rows[ix_s][2].append(s.time().isoformat())

    d += timedelta(days=1)

print(*rows, sep='\n')

This approach also reflects why two sunsets (or sunrises) can appear on a single UTC day. Event times vary and sometimes fall less than 24 hours apart; the cutoff at 00:00 UTC can slice them into the same calendar date when days and nights are short.

Single-day queries that don’t miss events

To reliably get all events for a specific date, evaluate the previous, the current, and the next day. Then keep only those results whose resulting datetime falls on the target date. This mirrors the range strategy in a minimal form.

def collect_sets_for_date(base_date, tz):
    out = []
    d = base_date - timedelta(days=1)
    for _ in range(3):
        val = sol.set_time(d, tz)
        if isinstance(val, datetime) and val.date() == base_date.date():
            out.append(val.time().isoformat())
        d += timedelta(days=1)
    return out

Why this detail matters

The algorithm intentionally simplifies the Sun’s apparent motion by holding right ascension and declination constant over 24-hour blocks keyed to midnight UTC. That’s convenient, but it amplifies edge effects where daily change is large. The RA wrap at the March equinox injects a false −24h step into T; range-forcing UT then hides which day an event actually belongs to and drops valid events when two occur in one UTC day. Range-forcing T instead removes the false step without discarding real information.

Takeaways

If you see missed sunrises or sunsets when forcing UT to [0, 24), stop forcing UT and normalize T to [0, 24) before converting to UT. When producing daily tables or answering single-day queries, compute across adjacent days and classify by the resulting datetime, not the input date. This preserves correct event attribution across day boundaries while avoiding artificial discontinuities introduced by the RA wrap.