SPA’s sunrise problem

Solarpositioning and its Rust counterpart use the NREL Solar Position Algorithm, or SPA. It gives very accurate solar positions, but its sunrise/sunset approximation can go wrong when the Sun barely crosses the horizon.

My implementation of it has had a few bugs too, which aren’t always easy to distinguish from limitations of the published algorithm. The replacement described here keeps SPA’s position calculation and uses it to search for sunrise, sunset and twilight events.

A precise position can still give a bad sunrise

SPA, by Ibrahim Reda and Afshin Andreas at NREL, calculates solar positions with a stated uncertainty of about 0.0003° (about one arcsecond). It draws on Jean Meeus’s Astronomical Algorithms.

Sunrise and sunset come from a separate procedure in the paper’s Appendix A.2, which samples solar coordinates on three neighbouring days, estimates the event times, interpolates the coordinates and applies a single correction. The correction usually works well, but it can overshoot badly when the Sun crosses the horizon at a shallow angle.

Consider one of the world’s northernmost permanent settlements: Longyearbyen, Svalbard, on 16 February 2020. Solarpositioning 2.1.2, using the SPA rise/set procedure, gives sunrise at about 10:09 UTC and sunset at 12:17. The US Naval Observatory’s annual table gives 10:32 and 11:54. That’s roughly 46 extra minutes of daylight (an attractive feature in February, but difficult to justify mathematically).

Evaluating SPA’s own position calculation at those returned times puts the solar centre roughly a quarter of a degree below the requested horizon.

Figure 1: Longyearbyen, 16 February 2020. The curve uses full SPA topocentric positions. The horizontal line is the selected −0.8333° horizon. The old event estimates fall well below it.

Appendix A.2 uses the Sun’s position relative to Earth’s centre (geocentric coordinates), while the new search described here uses SPA’s position for an observer at sea level (topocentric coordinates). This small shift in apparent direction (parallax) is actually not very significant: in this example it moves the crossing times by only 16 seconds, far too little to explain the old 23-minute errors.

Sunrise and sunset don’t always come in pairs

The initial hour-angle calculation in Appendix A.2 uses the Sun’s declination at the start of the day to decide whether a crossing is possible. Near a seasonal transition, the declination can change enough during the day to invalidate that decision.

At 66°34′ south, 0° longitude, on 5 January 2025, the old implementation reported continuous daylight. SPA’s position calculation shows a sunset late that evening, with no sunrise on that UTC date. The old API couldn’t represent this: it required either a rise/set pair or neither.

A calendar date can also contain parts of two solar cycles, with more than two events as a result. Longyearbyen on 16 April 2020 has a sunrise around 00:20, a sunset around 21:55 and another sunrise around 23:57 UTC.

At the equator, 179.9° east, the transits fall just before and just after 10 June 2020, leaving that UTC date without a transit even though it still contains a sunrise and a sunset. An API that requires a transit on the same date will lose those events, quite apart from any errors in SPA’s approximation.

Why search at all?

An analytical solution would be best, and it is available if we simplify the model. Hold the Sun’s declination (its angle north or south of the celestial equator) fixed while Earth rotates, and spherical geometry gives the rotation angle at which it reaches the chosen horizon. Converting that angle into a time gives the familiar sunrise estimate. Appendix A.2.4 of the SPA report introduces equation A4:

Calculate the local hour angle corresponding to the sun elevation equals −0.8333°, $H_0$,

$$H_0 = \arccos\left(\frac{\sin h'_0 - \sin\varphi\,\sin\delta_0}{\cos\varphi\,\cos\delta_0}\right).$$

Here, $\varphi$ is the observer’s latitude, $\delta_0$ is the Sun’s sampled declination, and $h'_0$ is the selected elevation (−0.8333° in SPA). The result $H_0$ is the angle from the local meridian to a crossing, with $\delta_0$ held fixed for this estimate.

The Sun’s coordinates change throughout the day, though, and the small observer-dependent corrections change too. Their values at sunrise depend on the time we’re trying to calculate. SPA’s full position calculation has no comparably simple inverse, so we solve elevation(time) = horizon numerically, evaluating positions at successive trial times. An analytical approximation can still supply a useful starting guess, but near a grazing crossing even a small change in the coordinates can shift the answer considerably.

This seems to be established practice in astronomical software. Skyfield numerically refines rise and set times using ephemeris positions. JPL’s SPICE Geometry Finder explicitly separates finding intervals that bracket events from refining the events within them. Finding all those intervals, including ones around closely spaced crossings, is the part that needs the sampling described below.

Asking for crossings

The basic operation is to find the next crossing of a given elevation within a time interval, with separate queries for rising and setting crossings.

For sunrise, the selected level is 50 arcminutes below the geometric horizon, the conventional allowance for the solar disc and mean refraction. Twilight uses a lower level; custom elevations use exactly the supplied angle. Since the horizon already accounts for refraction where required, the positions are calculated without an additional atmospheric correction. USNO describes the conventions here.

Transit has its own search for passage across the upper meridian. It isn’t needed to find sunrise or sunset, and it isn’t necessarily the instant of maximum elevation.

The low-level Java API looks like this:

var calculator = new SolarEvents(); // SPA by default
var sunrise = calculator.nextRise(
    start, end, latitude, longitude, deltaT,
    SolarEvents.Horizon.SUNRISE_SUNSET);

This returns an Optional<Instant>, empty if there is no crossing. The start is excluded and the end included, so a caller can continue from a returned event to find the next one. For a calendar date, forDate collects all events in the local day:

var events = calculator.forDate(
    LocalDate.of(2020, 4, 16), ZoneId.of("Arctic/Longyearbyen"),
    78.216667, 15.633333, 69.184);

The result contains rises, sets and transits lists (possibly empty or with several events), plus the initial state relative to the horizon. alwaysAbove() and alwaysBelow() indicate continuous daylight or night. The same call supports twilight, custom elevations and multiple horizons, with the time zone’s boundaries accounting for daylight-saving changes and skipped dates.

SPA is the default (and the model used throughout this post), but SolarEvents.grena3() selects Grena3 and a custom provider can supply elevation and hour angle. Accuracy and supported dates depend on the chosen model.

For the interface, I took ideas from Astronomy Engine’s bounded event searches and commons-suncalc’s explicit time windows and independent rise/set results. Time4J also treats absent sunrise and sunset as ordinary optional results.

Fixed steps can miss events

A straightforward search would sample positions every few minutes, look for a change of sign, then refine the crossing with bisection.

Unfortunately, a brief appearance above the horizon can fit entirely between two samples, leaving both below the horizon with a rise and set hidden in the middle. Choosing a smaller fixed step makes this less likely, but doesn’t remove the possibility.

Here’s the same Longyearbyen day with a custom horizon just below the daily maximum. The Sun stays above it for nearly two minutes, entirely between the five-minute samples.

Figure 2: Full SPA positions for Longyearbyen, 16 February 2020, with a custom horizon 0.0001° below the daily maximum.

How much could happen between two samples?

To avoid skipping such a pair, we need to account for how far the curve could bend between two samples. Draw a straight line between them, then leave room above and below it for that bend. If even the top of this band stays below the horizon, there is no need to look inside. If the band reaches the horizon, take another sample halfway through and repeat, searching the earlier half first.

In this approach, called adaptive subdivision, the band’s width comes from the standard error estimate for straight-line interpolation. Once an interval contains a rising crossing and the curve keeps rising throughout it, bisection narrows it to at most one millisecond. The same search finds sunset by reversing the signs.

The direct inspiration for this search was Don Cross’s Astronomy Engine. Its MaxAltitudeSlope supplies conservative experimental limits on how quickly the Sun’s elevation can change, and FindAscent subdivides intervals where a crossing could be hidden. I’ve adapted that approach to SPA, using a limit on how sharply the curve bends rather than how fast it moves.

The panels below show an interval being split, its left quarter discarded, then a rising crossing bracketed and refined through four halvings. The grey curve is included for illustration (the search only evaluates the points it needs).

Figure 3: Selected steps in the same rising-event search, with intermediate subdivisions omitted. Dotted lines join the endpoint samples. Blue bands show the room allowed for bending. Each panel has its own scale.

Each halving reduces the room needed for bending to a quarter of its previous size, apart from a small allowance for numerical rounding.

Choosing the margin

Since the Sun’s daily path is mainly the result of Earth’s steady rotation, we can use that rotation to estimate how sharply the curve bends, then round up generously. The Sun’s motion relative to the stars is much slower, but still needs an allowance of its own, particularly at the poles where the daily rise and fall becomes very small.

The choice of margin affects both the cost and the reliability of the search: a larger one requires extra position calculations, while a smaller one risks skipping events. I’ve chosen a generous margin and checked it against sampled positions and the tests below, though this remains an educated estimate rather than a formally proved bound.

The short numerical version

The code searches sin(elevation) − sin(horizon). Zero still means a crossing, and using the sine avoids awkward behaviour at the zenith. With time measured in hours, a once-per-day sine wave has a maximum curvature of about (2π/24)² = 0.069. The limit used here is 0.1 * abs(cos(latitude)) + 0.0001: round the daily contribution up to 0.1, reduce it towards the poles, and add 0.0001 for slower changes. Transit uses 0.1 without the latitude adjustment.

In a check of 100,000 positions across the supported years and latitudes, the largest sampled curvature was about 69% of the chosen allowance.

With a stopping interval of one millisecond, crossings closer together than that may fall within a single interval and cannot reliably be distinguished.

Does it work?

I compared the new search with the US Naval Observatory’s annual sunrise/sunset tables, covering every date in 2020 and 2025 at ten locations (7,310 location/date combinations). These include Singapore, Longyearbyen and McMurdo, plus test points near the polar circles and on either side of the date line.

Both calculations use a solar-centre elevation of −0.8333°, with no additional refraction correction or horizon dip. Events are collected by UTC date, combining neighbouring solar cycles for the old API, and times are rounded to the nearest minute to match the tables.

“Missing” and “extra” count events reported by only one side. Timing differences are calculated for matched events of the same kind.

Calculation	Missing events	Extra events	Largest timing difference
Released solarpositioning 2.1.2	9	11	23 minutes
New Java search	0	0	2 minutes

The new search matches all 12,472 reference events, including the two-sunrise date and the sunset previously discarded as continuous daylight.

A smaller comparison uses Skyfield with the JPL DE440s ephemeris, covering 56 combinations of date, location and horizon, including both poles. Skyfield calculates and searches its own solar positions, using the same time scales, sea-level observer and horizon, with refraction disabled. Its search functions also helped me work out how to separate event searches from calendar dates.

At the south pole, one crossing differs from JPL by about nine seconds because the Sun’s elevation changes so slowly that even a tiny angular difference becomes noticeable in the event time. To account for this, the test converts SPA’s stated angular uncertainty into a time tolerance using the local rate of elevation change, while still requiring event counts to match exactly.

Tests with mathematically known curves cover crossings about a second apart, a curve that just touches the horizon, and one that stays below it. They also check that splitting a search preserves its results (including a split at a returned event). Calendar tests cover daylight-saving changes, repeated and skipped dates, and midnight crossings.

For a wider check, a separate audit uses Skyfield’s maxima and minima search on SPA’s position curve, then finds crossings between those turning points. Across 48 date/location cases, all 1,816 rise/set searches found the same events as the new search, with times agreeing within two milliseconds. Custom horizons placed close to the extrema produced pairs as little as 31 milliseconds apart. This checks whether the search finds the crossings in SPA’s curve. The JPL comparison checks the underlying positions against a different model.

Cost and accuracy

Repeatedly evaluating SPA costs more than interpolating three days of coordinates. A local smoke test covers 42,822 combinations of date and coordinates with four horizons. In one run on my machine, the old calculation managed roughly 16,000 cases per second; the new first-event search managed about 1,300. The full daily API, collecting every event, took about 1.2 ms per location and date. The old API describes a solar cycle and the new one uses explicit date boundaries, so the workloads are comparable but not identical.

For this library, that seems a reasonable cost for using the precise positions we already calculate. Simpler, faster algorithms are available when that precision isn’t needed.

Refining the bracket to one millisecond doesn’t account for the effects of weather, the observer’s elevation or the local skyline on observed sunrise. Those uncertainties remain, along with the finite precision of the position model and the API’s approximation of UT1 with UTC.

Implementation status

Work on the Java and Rust libraries is in progress, currently in experimental branches.