Finding the sun is easier than finding sunrise

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 has had a few bugs too, which aren’t always easy to distinguish from limitations of the published algorithm. The replacement described in this post keeps SPA’s position calculation and uses it to search for sunrise, sunset and twilight events. The main design goals are correctness, accuracy, and simplicity. Efficiency is a lesser consideration, as there are plenty of fast algorithms that give decent results for everyday use.

Precise positioning isn’t enough

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 and has become something of a reference algorithm.

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 uses SPA’s position for a sea-level observer (topocentric coordinates). This shift in apparent direction (parallax) moves the crossing times by only 16 seconds in this example, 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 span parts of two solar cycles and contain more than two events. On 16 April 2020, Longyearbyen has a sunrise around 00:20, a sunset around 21:55 and another sunrise around 23:57 UTC.

At the equator, 179.9° east, transits fall just before and after 10 June 2020, leaving that UTC date with a sunrise and sunset but no transit. 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?

We could solve a simplified model analytically by holding the Sun’s declination (its angle north or south of the celestial equator) fixed while Earth rotates. Appendix A.2.4 of the SPA report gives the crossing angle in 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$ the sampled declination, and $h'_0$ the selected elevation (−0.8333° in SPA). $H_0$ is the angle from the local meridian to the crossing, which we can convert into a time.

The Sun’s coordinates and observer-dependent corrections change throughout the day, though. 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 evaluate positions at trial times to solve elevation(time) = horizon numerically. The approximation can give us a starting guess, but near a grazing crossing even small coordinate changes can shift the answer considerably.

Other astronomical software uses numerical searches too. Skyfield refines rise and set times from ephemeris positions, while JPL’s SPICE Geometry Finder separates finding intervals that bracket events from refining the crossings inside them. Finding those intervals, especially around closely spaced crossings, 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.

A nice thing about this approach: search is independent of the positioning algorithm. It needs elevation and hour angle at a given time, whether those come from SPA, Grena3 or another solar model (subject to the smoothness limits below). The choice of model determines the position accuracy and supported dates.

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 = SolarEvents.spa();
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.

The examples here use SPA. SolarEvents.grena3() selects Grena3, and a custom provider can supply elevation and hour angle without changing the event search.

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 it. A smaller fixed step makes this less likely, but still possible.

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, we have a standard root-finding problem, with established solutions such as bisection and Brent’s method.

I use the newer ITP method by Oliveira and Takahashi (interpolate, truncate, project). It uses interpolation for speed, with safeguards that keep the worst-case number of steps close to bisection’s.

The search stops when the time interval is at most one millisecond wide, with the calculated Sun on opposite sides of the chosen horizon. 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 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 ITP steps. The grey curve is 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.

During subdivision, each halving reduces the room needed for bending to a quarter of its previous size (apart from a small allowance for 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.

A larger margin 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 comparisons 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 one-millisecond stopping interval, closer crossings may fall in the same interval and cannot reliably be distinguished.

Does it work?

The comparison here is between the SPA paper’s Appendix A.2 sunrise/sunset approximation (as implemented in solarpositioning 2.1.2) and our new numerical search. The new search uses SPA’s precise position calculation to find horizon crossings, replacing the Appendix A.2 procedure.

Against the USNO tables

I compared both implementations with the US Naval Observatory’s annual sunrise/sunset tables, covering every date in 2020 and 2025 at ten locations. These range from Singapore to Longyearbyen and McMurdo (a research station in Antarctica), with test points near the polar circles and on either side of the date line. I used the same standard sunrise/sunset horizon, grouped events by UTC date and rounded the calculated times to the nearest minute to match the tables.

The new search finds all 12,472 reference events, with no extras. The old implementation misses nine and reports eleven that aren’t in the tables. For events found by both sides, the largest timing differences look like this:

Figure 4: Largest absolute difference from USNO at each location, after rounding to minutes. Missing and extra events are counted separately in the text.

Against JPL positions

To compare times in seconds, I checked the new search against Skyfield using the JPL DE440 ephemeris (the smaller DE440s file). Skyfield calculates and searches its own solar positions. This comparison covers 14 date/location combinations at all four horizons (sunrise/sunset and the three twilights), including both poles, with matching time scales and sea-level observers.

Both searches find the same events. Of the 59 horizon crossings, 57 agree within 0.15 seconds. The two exceptions are sunrises at the exact poles: about 1.3 seconds apart at the north pole and 8.9 seconds at the south pole. There, the Sun’s elevation changes so slowly that the tiniest difference in the calculated position makes a noticeable difference to the crossing time.

Figure 5: Absolute timing differences between the new search with SPA positions and Skyfield/JPL. Each point is one crossing, ordered by difference. The vertical scale is logarithmic (each major step is a factor of ten).

These are comparisons between calculations under the same assumptions. As always, it bears repeating that agreement to a fraction of a second doesn’t mean we can predict an observed sunrise or sunset so closely. Local topography and weather will have a significant influence on that.

The cost of accuracy

Repeatedly evaluating SPA costs much more than interpolating three days of coordinates. A local test covers 42,822 combinations of date and coordinates at all four standard horizons:

Calculation	Cases/s	Time per case
SPA Appendix A.2 approximation	15,600	0.064 ms
Search with SPA positions, bisection	1,400	0.71 ms
Search with SPA positions, ITP	1,610	0.62 ms
Search with Grena3 positions, bisection	17,150	0.058 ms
Search with Grena3 positions, ITP	18,970	0.053 ms

These are medians of warmed-up runs on my machine, not a general performance claim. The old API returns one solar cycle, while the new API collects every event in the requested calendar date, so the workloads are comparable but not identical. ITP improves SPA throughput by about 15% over bisection, but the full SPA search still takes roughly ten times as long as the original approximation.

For this library, that seems a reasonable cost for using the precise positions we already calculate.