Return sunrise and sunset at the observer's element type

transit_sunrise_sunset rounds the event times to a whole second in order to
build a DateTime. That discards the sub-second part, and with it anything the
element type carries beyond a value. The solve already runs at the observer's
element type, so the information exists right up to that final conversion and
is then thrown away.

transit_sunrise_sunset_seconds returns the same events as seconds since
midnight UTC, keeping the element type. A Measurement observer therefore yields
event times with an uncertainty, and a dual-valued observer yields
differentiable ones. Neither was possible before, and the DateTime variant could
not even be finite-differenced, since a step small enough to be a derivative is
swamped by the rounding.

The conversion now dispatches on the requested type, so the DateTime and
ZonedDateTime paths are untouched, including the polar day and night branch.
TransitSunriseSunset needed no change, having been made generic when TimeZones
moved to an extension.

Covered from both directions: the propagated uncertainty is checked against a
finite-difference gradient, the ForwardDiff derivative against the same, and the
two variants are asserted to describe the same instant for a Float64 observer.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01NQ9zmXUHQtA5suiMB8UdA7

Author: langestefan

docs/src/index.md

@@ -172,10 +172,21 @@ The geometric angles come back exact, since they do not depend on the atmosphere
 the apparent ones carry the pressure and temperature uncertainty. Making the observer
 uncertain as well widens both, and the two contributions combine in the apparent angles.
 
-!!! note
-    [`transit_sunrise_sunset`](@ref) and the `next_`/`previous_` helpers return a
-    `DateTime`, which cannot carry an uncertainty. They compute from the nominal
-    coordinates and drop it.
+### Sunrise and sunset
+
+Event times are a special case, because a `DateTime` cannot carry an uncertainty.
+[`transit_sunrise_sunset`](@ref) rounds to a whole second and so returns the nominal
+times. [`transit_sunrise_sunset_seconds`](@ref) keeps the observer's element type
+instead, giving seconds since midnight UTC:
+
+```@example srt
+transit_sunrise_sunset_seconds(Observer(52.35888 ± 0.01, 4.88185 ± 0.01), Date(2023, 6, 21))
+```
+
+Rounding to a whole second also costs the sub-second part, so this variant is the more
+precise one even with ordinary numbers, and it is the one to use with ForwardDiff. The
+`next_`/`previous_` helpers return a `DateTime` and drop the uncertainty in the same way
+[`transit_sunrise_sunset`](@ref) does.
 
 ## Numeric precision
 

docs/src/utilities.md

@@ -289,5 +289,6 @@ SolarPosition.Utilities.previous_solar_noon
 
 ```@docs
 SolarPosition.Utilities.transit_sunrise_sunset
+SolarPosition.Utilities.transit_sunrise_sunset_seconds
 SolarPosition.Utilities.TransitSunriseSunset
 ```

src/Utilities/Utilities.jl

@@ -9,6 +9,7 @@ include("srt.jl")
 
 export TransitSunriseSunset,
     transit_sunrise_sunset,
+    transit_sunrise_sunset_seconds,
     next_sunrise,
     next_sunset,
     next_solar_noon,

src/Utilities/spa.jl

@@ -7,6 +7,12 @@ using ..Positioning: _compute_spa_srt_parameters
 _frac_to_dt(dt_midnight, frac) =
     dt_midnight + Dates.Second(round(Int, frac * 86_400))
 
+# Render a fraction of the day as the requested event type. Rounding to a whole second
+# discards the sub-second part and anything the element type carries beyond a value, such
+# as an uncertainty or a derivative, so a real-valued request keeps the number instead.
+_frac_to_event(::Type{R}, dt_midnight, frac) where {R} = _frac_to_dt(dt_midnight, frac)
+_frac_to_event(::Type{R}, _, frac) where {R <: Real} = frac * 86_400
+
 # Returns (ν, α, δ): apparent sidereal time, geocentric right ascension, declination (degrees).
 function _compute_srt_parameters(::Type{T}, dt::DateTime, δt) where {T <: Real}
     srt = _compute_spa_srt_parameters(T, dt, δt)
@@ -22,6 +28,15 @@ function _transit_sunrise_sunset(
     return _transit_sunrise_sunset_impl(R, obs, dt, alg, nothing)
 end
 
+function _transit_sunrise_sunset(
+        ::Type{R},
+        obs::Observer{T},
+        dt::DateTime,
+        alg::SPA,
+    ) where {T <: Real, R <: Real}
+    return _transit_sunrise_sunset_impl(R, obs, dt, alg, nothing)
+end
+
 function _transit_sunrise_sunset_impl(
         ::Type{R},
         obs::Observer{T},
@@ -86,7 +101,8 @@ function _transit_sunrise_sunset_impl(
             "polar day (sun above horizon)"
         @warn "Sun does not rise or set on this date at the given location: $polar_condition. Returning midnight UTC for all events." _group =
             :polar_day_night maxlog = 1
-        return TransitSunriseSunset{R}(dt_midnight, dt_midnight, dt_midnight, tz)
+        midnight = _frac_to_event(R, dt_midnight, zero(obs.latitude))
+        return TransitSunriseSunset{R}(midnight, midnight, midnight, tz)
     end
 
     H0 = acosd(cos_H0_arg)
@@ -168,12 +184,11 @@ function _transit_sunrise_sunset_impl(
         S_frac += 1
     end
 
-    # Convert fractions of day to DateTime
-    # Each fraction represents seconds into the day from midnight UTC
+    # Each fraction represents the offset into the day from midnight UTC
     return TransitSunriseSunset{R}(
-        _frac_to_dt(dt_midnight, T_frac),
-        _frac_to_dt(dt_midnight, R_frac),
-        _frac_to_dt(dt_midnight, S_frac),
+        _frac_to_event(R, dt_midnight, T_frac),
+        _frac_to_event(R, dt_midnight, R_frac),
+        _frac_to_event(R, dt_midnight, S_frac),
         tz,
     )
 end

src/Utilities/srt.jl

@@ -43,6 +43,16 @@ function TransitSunriseSunset{DateTime}(
     return TransitSunriseSunset{DateTime}(transit, sunrise, sunset)
 end
 
+# Constructor for a real element type, where the events are seconds since midnight UTC
+function TransitSunriseSunset{T}(
+        transit::T,
+        sunrise::T,
+        sunset::T,
+        ::Nothing,
+    ) where {T <: Real}
+    return TransitSunriseSunset{T}(transit, sunrise, sunset)
+end
+
 """Calculate the sun transit, sunrise, and sunset
 for a given date at an Observer location.
 """
@@ -62,6 +72,41 @@ function transit_sunrise_sunset(
     return transit_sunrise_sunset(obs, DateTime(dt), alg)
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Calculate the sun transit, sunrise, and sunset as seconds since midnight UTC.
+
+[`transit_sunrise_sunset`](@ref) rounds to a whole second in order to return a `DateTime`.
+This variant keeps the event times at the observer's element type instead, so sub-second
+precision survives, and so does anything the element type carries beyond a value. A
+`Measurement` observer therefore yields event times with an uncertainty, and a dual-valued
+observer yields differentiable ones.
+
+# Arguments
+- `obs`: Observer location
+- `dt`: date of interest, any time of day is ignored
+- `alg`: Solar algorithm to use (default: SPA())
+
+# Returns
+- `TransitSunriseSunset{T}` whose fields are seconds since midnight UTC
+"""
+function transit_sunrise_sunset_seconds(
+        obs::Observer{T},
+        dt::DateTime,
+        alg::SolarAlgorithm = SPA(),
+    )::TransitSunriseSunset{T} where {T <: Real}
+    return _transit_sunrise_sunset(T, obs, dt, alg)
+end
+
+function transit_sunrise_sunset_seconds(
+        obs::Observer{T},
+        dt::Date,
+        alg::SolarAlgorithm = SPA(),
+    )::TransitSunriseSunset{T} where {T <: Real}
+    return transit_sunrise_sunset_seconds(obs, DateTime(dt), alg)
+end
+
 # Helper function for next_* functions
 function _next_event(obs::Observer, dt::DateTime, alg::SolarAlgorithm, event_field::Symbol)
     date_only = Date(dt)

test/positioning/test-autodiff.jl

@@ -50,6 +50,23 @@ using StructArrays: StructArrays
         @test ad ≈ fd atol = 1.0e-8
     end
 
+    # only the seconds variant is differentiable. The DateTime one rounds to a whole
+    # second, which flattens the derivative to zero almost everywhere and makes it a step
+    # function elsewhere, so the finite-difference reference uses the seconds variant too
+    @testset "Sunrise and sunset derivatives" begin
+        day = DateTime(2024, 6, 21)
+        h = 0.05
+        @testset "$field" for field in (:transit, :sunrise, :sunset)
+            f = lat -> getproperty(
+                transit_sunrise_sunset_seconds(Observer(lat, 10.0), day), field,
+            )
+            ad = ForwardDiff.derivative(f, 45.0)
+            fd = (f(45.0 + h) - f(45.0 - h)) / 2h
+            @test isfinite(ad)
+            @test ad ≈ fd rtol = 1.0e-3
+        end
+    end
+
     @testset "Second derivatives via nested duals" begin
         f = x -> solar_position(Observer(x[1], x[2]), dt, PSA(), NoRefraction()).elevation
         hess = ForwardDiff.hessian(f, [45.0, 10.0])

test/positioning/test-measurements.jl

@@ -87,17 +87,50 @@ using StructArrays: StructArrays
         @test df.elevation == pos.elevation
     end
 
-    # a DateTime cannot carry an uncertainty, so the sunrise and sunset utilities compute
-    # from the nominal coordinates and drop it. This test pins that down as deliberate
-    # behaviour rather than leaving it to be rediscovered.
-    @testset "Sunrise and sunset drop the uncertainty" begin
+    @testset "Sunrise and sunset" begin
+        day = DateTime(2024, 6, 21)
         nominal = Observer(45.0, 10.0)
-        loose = Observer(45.0 ± 5.0, 10.0 ± 5.0)
-        for f in (next_sunrise, next_sunset, next_solar_noon)
-            @test f(loose, dt) == f(nominal, dt)
+        loose = Observer(45.0 ± σ, 10.0 ± σ)
+
+        # a DateTime cannot carry an uncertainty, so this API computes from the nominal
+        # coordinates and drops it
+        @testset "the DateTime variant drops it" begin
+            wide = Observer(45.0 ± 5.0, 10.0 ± 5.0)
+            for f in (next_sunrise, next_sunset, next_solar_noon)
+                @test f(wide, dt) == f(nominal, dt)
+            end
+            a = transit_sunrise_sunset(wide, Date(2024, 6, 21))
+            b = transit_sunrise_sunset(nominal, Date(2024, 6, 21))
+            @test (a.transit, a.sunrise, a.sunset) == (b.transit, b.sunrise, b.sunset)
+        end
+
+        @testset "the seconds variant keeps it" begin
+            ev = transit_sunrise_sunset_seconds(loose, day)
+            @test ev isa TransitSunriseSunset{Measurement{Float64}}
+            for field in (:transit, :sunrise, :sunset)
+                @test uncertainty(getfield(ev, field)) > 0
+            end
+
+            # the reference gradient is built from the seconds variant as well, since a
+            # step small enough for the DateTime variant is swamped by its rounding
+            h = 0.05
+            sr(lat, lon) = transit_sunrise_sunset_seconds(Observer(lat, lon), day).sunrise
+            ∂lat = (sr(45.0 + h, 10.0) - sr(45.0 - h, 10.0)) / 2h
+            ∂lon = (sr(45.0, 10.0 + h) - sr(45.0, 10.0 - h)) / 2h
+            @test uncertainty(ev.sunrise) ≈ hypot(∂lat * σ, ∂lon * σ) rtol = 1.0e-3
+        end
+
+        # the two variants must describe the same instant, the seconds one simply keeps
+        # the sub-second part that rounding to a DateTime throws away
+        @testset "the variants agree for a Float64 observer" begin
+            ev = transit_sunrise_sunset_seconds(nominal, day)
+            ref = transit_sunrise_sunset(nominal, day)
+            midnight = DateTime(Date(day))
+            for field in (:transit, :sunrise, :sunset)
+                secs = getfield(ev, field)
+                @test midnight + Second(round(Int, secs)) == getfield(ref, field)
+                @test secs != round(secs)
+            end
         end
-        a = transit_sunrise_sunset(loose, Date(2024, 6, 21))
-        b = transit_sunrise_sunset(nominal, Date(2024, 6, 21))
-        @test (a.transit, a.sunrise, a.sunset) == (b.transit, b.sunrise, b.sunset)
     end
 end