From 3234d1cf077dedf277e873f2e089b9c27340cee7 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <s.e.d.lange@tue.nl>
Date: Tue, 28 Jul 2026 12:05:16 +0200
Subject: [PATCH 01/10] Extract the topocentric half of SPA into a shared
 function

Pure code motion: _spa_topocentric now holds the hour angle, parallax,
elevation, and azimuth steps so another algorithm can reconstruct
positions from a given geocentric state through the exact same code.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Wg3j4XXUooMvDgUpRajxTq
---
 src/Positioning/spa.jl | 34 ++++++++++++++++++++--------------
 1 file changed, 20 insertions(+), 14 deletions(-)

diff --git a/src/Positioning/spa.jl b/src/Positioning/spa.jl
index d1abc02d..d765b857 100644
--- a/src/Positioning/spa.jl
+++ b/src/Positioning/spa.jl
@@ -361,20 +361,9 @@ function _compute_spa_srt_parameters(::Type{T}, dt::DateTime, δt) where {T <: R
     return (; ν, α, δ, R, ε, δψ, jme)
 end
 
-function _solar_position(
-        obs::SPAObserver{T},
-        dt::DateTime,
-        alg::SPA,
-    ) where {T <: Real}
-    δt::T = if alg.delta_t === nothing
-        calculate_deltat(T, dt)
-    else
-        T(alg.delta_t)
-    end
-
-    # Compute sidereal time, right ascension, declination, and related parameters
-    ν, α, δ, R, ε, δψ, jme = _compute_spa_srt_parameters(T, dt, δt)
-
+# Topocentric half of SPA, shared by SPA and Interpolated so both paths run the exact
+# same code from a given geocentric state (ν, α, δ, R).
+function _spa_topocentric(obs::SPAObserver{T}, ν::T, α::T, δ::T, R::T) where {T <: Real}
     # observer local hour angle
     H = local_hour_angle(ν, obs.longitude, α)
     H_rad = deg2rad(H)
@@ -409,6 +398,23 @@ function _solar_position(
     return SolPos{T}(az, e0, θz0)
 end
 
+function _solar_position(
+        obs::SPAObserver{T},
+        dt::DateTime,
+        alg::SPA,
+    ) where {T <: Real}
+    δt::T = if alg.delta_t === nothing
+        calculate_deltat(T, dt)
+    else
+        T(alg.delta_t)
+    end
+
+    # Compute sidereal time, right ascension, declination, and related parameters
+    (; ν, α, δ, R) = _compute_spa_srt_parameters(T, dt, δt)
+
+    return _spa_topocentric(obs, ν, α, δ, R)
+end
+
 function _solar_position(obs::Observer{T}, dt::DateTime, alg::SPA) where {T <: Real}
     spa_obs = SPAObserver{T}(obs.latitude, obs.longitude, obs.altitude)
     return _solar_position(spa_obs, dt, alg)

From 3f52b2f1b768fee98d282c684d25603dc8dfe2f8 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <s.e.d.lange@tue.nl>
Date: Tue, 28 Jul 2026 12:05:26 +0200
Subject: [PATCH 02/10] Add the Interpolated algorithm and solar_rate

Interpolated wraps SPA with cubic B-splines of its geocentric right
ascension, declination, radius vector, and equation of the equinoxes on
a uniform time grid, then reconstructs topocentric positions through the
shared _spa_topocentric path. Sidereal time stays closed form, so the
only error is the spline error of four smooth series, about 1e-10
degrees at the default one hour step, and one interpolant serves every
observer. Scalar evaluation is allocation free and about 10x faster
than direct SPA.

Construction lives in a new Interpolations.jl package extension; the
struct, evaluation, out_of_range handling (:error or :fallback), and
the finite difference solar_rate live in the core package.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Wg3j4XXUooMvDgUpRajxTq
---
 Project.toml                           |   3 +
 ext/SolarPositionInterpolationsExt.jl  |  57 +++++++
 src/Positioning/Positioning.jl         |   2 +
 src/Positioning/interpolated.jl        | 211 +++++++++++++++++++++++++
 test/Project.toml                      |   2 +
 test/extensions/test-interpolations.jl | 181 +++++++++++++++++++++
 6 files changed, 456 insertions(+)
 create mode 100644 ext/SolarPositionInterpolationsExt.jl
 create mode 100644 src/Positioning/interpolated.jl
 create mode 100644 test/extensions/test-interpolations.jl

diff --git a/Project.toml b/Project.toml
index 9b5787f4..92a78b58 100644
--- a/Project.toml
+++ b/Project.toml
@@ -15,12 +15,14 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 TimeZones = "f269a46b-ccf7-5d73-abea-4c690281aa53"
 
 [weakdeps]
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [extensions]
+SolarPositionInterpolationsExt = "Interpolations"
 SolarPositionMakieExt = "Makie"
 SolarPositionModelingToolkitExt = ["ModelingToolkit", "Symbolics"]
 SolarPositionOhMyThreadsExt = "OhMyThreads"
@@ -29,6 +31,7 @@ SolarPositionOhMyThreadsExt = "OhMyThreads"
 Aqua = "0.8"
 Dates = "1"
 DocStringExtensions = "0.8, 0.9"
+Interpolations = "0.15, 0.16"
 Makie = "0.24"
 ModelingToolkit = "11"
 OhMyThreads = "0.8"
diff --git a/ext/SolarPositionInterpolationsExt.jl b/ext/SolarPositionInterpolationsExt.jl
new file mode 100644
index 00000000..6dd6d481
--- /dev/null
+++ b/ext/SolarPositionInterpolationsExt.jl
@@ -0,0 +1,57 @@
+"""
+    SolarPositionInterpolationsExt
+
+Extension that builds the cubic B-spline interpolants behind the `Interpolated` solar
+position algorithm. It activates when Interpolations.jl is loaded and implements the
+sampling of the wrapped algorithm's geocentric solar coordinates on a uniform time
+grid. Evaluation lives in the main package and needs no extension.
+"""
+module SolarPositionInterpolationsExt
+
+using SolarPosition.Positioning: Positioning, SPA, calculate_deltat
+using Interpolations: interpolate, scale, BSpline, Cubic, Line, OnGrid
+using Dates: Dates, DateTime, Millisecond
+using Base.Threads: @threads
+
+function Positioning._build_interpolants(
+        algorithm::SPA, t0::DateTime, t1::DateTime, step::Millisecond,
+    )
+    # pad the grid two steps beyond each end so queries at the exact span endpoints
+    # stay away from the spline's less accurate boundary cells
+    t0p = t0 - 2 * step
+    t1p = t1 + 2 * step
+    n = Int(cld(Dates.value(t1p) - Dates.value(t0p), Dates.value(step))) + 1
+
+    αs = Vector{Float64}(undef, n)
+    δs = Vector{Float64}(undef, n)
+    Rs = Vector{Float64}(undef, n)
+    eqs = Vector{Float64}(undef, n)
+    @threads for i in 1:n
+        dti = t0p + (i - 1) * step
+        δt = if algorithm.delta_t === nothing
+            calculate_deltat(Float64, dti)
+        else
+            Float64(algorithm.delta_t)
+        end
+        p = Positioning._compute_spa_srt_parameters(Float64, dti, δt)
+        αs[i] = p.α
+        δs[i] = p.δ
+        Rs[i] = p.R
+        eqs[i] = p.δψ * cosd(p.ε)
+    end
+
+    # unwrap right ascension so the fitted series is continuous across 0/360. This is
+    # sequential and relies on the constructor's step cap keeping the per step advance
+    # far below half a turn.
+    for i in 2:n
+        αs[i] -= 360.0 * round((αs[i] - αs[i - 1]) / 360.0)
+    end
+
+    x0 = Positioning.julian_day_j2000(Float64, t0p)
+    dx = Dates.value(step) / 86_400_000
+    xs = range(x0; step = dx, length = n)
+    mk(v) = scale(interpolate(v, BSpline(Cubic(Line(OnGrid())))), xs)
+    return (mk(αs), mk(δs), mk(Rs), mk(eqs))
+end
+
+end
diff --git a/src/Positioning/Positioning.jl b/src/Positioning/Positioning.jl
index bd6e1f84..e225cc34 100644
--- a/src/Positioning/Positioning.jl
+++ b/src/Positioning/Positioning.jl
@@ -488,9 +488,11 @@ include("noaa.jl")
 include("walraven.jl")
 include("usno.jl")
 include("spa.jl")
+include("interpolated.jl")
 
 export Observer,
     PSA, NOAA, Walraven, USNO, SPA, solar_position, solar_position!, SolPos, ApparentSolPos
+export Interpolated, solar_rate
 export SolarAlgorithm, AbstractSolPos, AbstractApparentSolPos
 export calculate_deltat
 
diff --git a/src/Positioning/interpolated.jl b/src/Positioning/interpolated.jl
new file mode 100644
index 00000000..c635bf21
--- /dev/null
+++ b/src/Positioning/interpolated.jl
@@ -0,0 +1,211 @@
+"""
+    $(TYPEDEF)
+
+Fast solar position from cubic B-spline interpolation of the geocentric solar
+coordinates of a wrapped exact algorithm, currently [`SPA`](@ref).
+
+The interpolant samples four observer independent quantities on a uniform time grid:
+geocentric right ascension, geocentric declination, the Earth-Sun radius vector, and
+the equation of the equinoxes. Sidereal time and the topocentric reconstruction stay
+closed form and run through the exact same code as the wrapped algorithm, so the only
+error is the spline interpolation error of the four smooth geocentric series. At the
+default one hour grid spacing this error is far below the wrapped algorithm's own
+accuracy.
+
+One interpolant serves every observer, so a single instance can be shared across a
+whole grid of sites and across threads.
+
+Construction requires Interpolations.jl to be loaded:
+
+    using Interpolations
+    alg = Interpolated(SPA(); tspan = (DateTime(2024, 1, 1), DateTime(2025, 1, 1)))
+    solar_position(obs, dt, alg)
+
+# Arguments
+
+- `algorithm::SPA = SPA()`: the exact algorithm to sample. Its `delta_t` setting is
+  used during sampling. The default constant ΔT keeps the sampled series smooth; with
+  `delta_t = nothing` the piecewise ΔT model introduces slope kinks that are far below
+  the interpolation tolerance.
+- `tspan::Tuple`: the valid query span as a pair of `DateTime` or `ZonedDateTime`.
+  Zoned times are converted to UTC.
+- `step::Period = Hour(1)`: the sampling grid spacing. Must be positive and at most 30
+  days so that right ascension advances much less than half a turn per step.
+- `out_of_range::Symbol = :error`: behaviour for queries outside `tspan`. `:error`
+  throws an `ArgumentError`, `:fallback` silently calls the wrapped exact algorithm.
+
+# Fields
+$(TYPEDFIELDS)
+"""
+struct Interpolated{A <: SolarAlgorithm, ITP} <: SolarAlgorithm
+    "Wrapped exact algorithm, used for sampling and the `:fallback` path"
+    algorithm::A
+    "Cubic B-spline of unwrapped geocentric right ascension in degrees vs days since J2000"
+    itp_α::ITP
+    "Cubic B-spline of geocentric declination in degrees"
+    itp_δ::ITP
+    "Cubic B-spline of the Earth-Sun radius vector in astronomical units"
+    itp_R::ITP
+    "Cubic B-spline of the equation of the equinoxes in degrees"
+    itp_eqeq::ITP
+    "Valid query span as UTC datetimes"
+    tspan::Tuple{DateTime, DateTime}
+    "Sampling grid spacing"
+    step::Dates.Millisecond
+    "Span start in days since J2000"
+    t_min::Float64
+    "Span end in days since J2000"
+    t_max::Float64
+    "Out of range behaviour, `:error` or `:fallback`"
+    out_of_range::Symbol
+end
+
+function Interpolated(
+        algorithm::SPA = SPA();
+        tspan::Tuple{<:Union{DateTime, ZonedDateTime}, <:Union{DateTime, ZonedDateTime}},
+        step::Dates.Period = Dates.Hour(1),
+        out_of_range::Symbol = :error,
+    )
+    out_of_range in (:error, :fallback) || throw(
+        ArgumentError("out_of_range must be :error or :fallback, got :$out_of_range"),
+    )
+    t0 = _as_utc(tspan[1])
+    t1 = _as_utc(tspan[2])
+    t0 < t1 || throw(ArgumentError("tspan must be increasing, got $t0 to $t1"))
+    stepms = Dates.Millisecond(step)
+    # the cap keeps the per step advance of right ascension far below half a turn,
+    # which the unwrap during construction relies on
+    Dates.Millisecond(0) < stepms <= Dates.Millisecond(Dates.Day(30)) ||
+        throw(ArgumentError("step must be positive and at most 30 days, got $step"))
+    (itp_α, itp_δ, itp_R, itp_eqeq) = _build_interpolants(algorithm, t0, t1, stepms)
+    return Interpolated(
+        algorithm, itp_α, itp_δ, itp_R, itp_eqeq, (t0, t1), stepms,
+        julian_day_j2000(Float64, t0), julian_day_j2000(Float64, t1), out_of_range,
+    )
+end
+
+_as_utc(dt::DateTime) = dt
+_as_utc(zdt::ZonedDateTime) = DateTime(zdt, UTC)
+
+# Extension hook. SolarPositionInterpolationsExt defines the working method for SPA;
+# this fallback exists so construction without the extension fails with a clear message.
+function _build_interpolants(
+        ::SolarAlgorithm, ::DateTime, ::DateTime, ::Dates.Millisecond,
+    )
+    throw(
+        ArgumentError(
+            "constructing Interpolated requires Interpolations.jl. " *
+                "Run `using Interpolations` and try again.",
+        ),
+    )
+end
+
+function Base.show(io::IO, alg::Interpolated)
+    print(
+        io, "Interpolated(", alg.algorithm, "; tspan = ", alg.tspan,
+        ", step = ", alg.step, ", out_of_range = :", alg.out_of_range, ")",
+    )
+    return nothing
+end
+
+function _solar_position(
+        obs::SPAObserver{T},
+        dt::DateTime,
+        alg::Interpolated,
+    ) where {T <: Real}
+    t = julian_day_j2000(Float64, dt)
+    if !(alg.t_min <= t <= alg.t_max)
+        alg.out_of_range === :fallback && return _solar_position(obs, dt, alg.algorithm)
+        throw(
+            ArgumentError(
+                "$dt is outside the interpolated span $(alg.tspan[1]) to " *
+                    "$(alg.tspan[2]). Widen tspan or use out_of_range = :fallback.",
+            ),
+        )
+    end
+
+    # interpolated geocentric state, right ascension rewrapped to [0, 360)
+    α = mod(alg.itp_α(t), 360.0)
+    δ = alg.itp_δ(t)
+    R = alg.itp_R(t)
+    eqeq = alg.itp_eqeq(t)
+
+    # apparent sidereal time stays closed form and magnitude safe
+    (n_int, n_frac) = julian_day_j2000_split(Float64, dt)
+    jc = (n_int + n_frac) / 36525.0
+    ν = mean_sidereal_time(n_int, n_frac, jc) + eqeq
+
+    return _spa_topocentric(obs, T(ν), T(α), T(δ), T(R))
+end
+
+function _solar_position(obs::Observer{T}, dt::DateTime, alg::Interpolated) where {T <: Real}
+    spa_obs = SPAObserver{T}(obs.latitude, obs.longitude, obs.altitude)
+    return _solar_position(spa_obs, dt, alg)
+end
+
+# DefaultRefraction resolves exactly as it does for the wrapped SPA
+function _solar_position(
+        obs::AbstractObserver{T},
+        dt,
+        alg::Interpolated,
+        ::DefaultRefraction,
+    ) where {T <: Real}
+    spa = alg.algorithm
+    return _solar_position(
+        obs,
+        dt,
+        alg,
+        SPARefraction{T}(
+            pressure = T(spa.pressure),
+            temperature = T(spa.temperature),
+            atmos_refract = T(spa.atmos_refract),
+        ),
+    )
+end
+
+function solar_position!(
+        pos::StructArrays.StructVector{S},
+        obs::AbstractObserver{T},
+        dts::AbstractVector{DateTime},
+        alg::Interpolated,
+        refraction::RefractionAlgorithm = DefaultRefraction(),
+    ) where {S <: AbstractSolPos, T <: Real}
+    spa_obs = SPAObserver{T}(obs.latitude, obs.longitude, obs.altitude)
+    @inbounds for i in eachindex(dts, pos)
+        pos[i] = solar_position(spa_obs, dts[i], alg, refraction)
+    end
+    return pos
+end
+
+# Interpolated mirrors SPA: SolPos with NoRefraction, ApparentSolPos with any refraction
+result_type(::Type{<:Interpolated}, ::Type{NoRefraction}, ::Type{T}) where {T} = SolPos{T}
+result_type(::Type{<:Interpolated}, ::Type{<:RefractionAlgorithm}, ::Type{T}) where {T} =
+    ApparentSolPos{T}
+
+"""
+    $(TYPEDSIGNATURES)
+
+Rate of change of solar azimuth and elevation in degrees per hour at `dt`, computed
+with a central finite difference of one second half width over the interpolated
+position. Returns a named tuple `(dazimuth_dt, delevation_dt)`. The azimuth difference
+is wrap aware, so the rate is continuous across the 0/360 degree crossing.
+
+`dt` must be at least one second inside the interpolated span, unless the algorithm
+was constructed with `out_of_range = :fallback`.
+"""
+function solar_rate(
+        obs::AbstractObserver{T},
+        dt::DateTime,
+        alg::Interpolated,
+    ) where {T <: Real}
+    spa_obs = SPAObserver{T}(obs.latitude, obs.longitude, obs.altitude)
+    p1 = _solar_position(spa_obs, dt - Dates.Second(1), alg)
+    p2 = _solar_position(spa_obs, dt + Dates.Second(1), alg)
+    # 1800 converts the difference over a two second baseline to degrees per hour
+    daz = mod(p2.azimuth - p1.azimuth + 180, 360) - 180
+    return (dazimuth_dt = daz * 1800, delevation_dt = (p2.elevation - p1.elevation) * 1800)
+end
+
+function solar_rate(obs::AbstractObserver, dt::ZonedDateTime, alg::Interpolated)
+    return solar_rate(obs, DateTime(dt, UTC), alg)
+end
diff --git a/test/Project.toml b/test/Project.toml
index bdd2f087..5044dd1d 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,6 +4,7 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
@@ -21,5 +22,6 @@ SolarPosition = {path = ".."}
 Aqua = "0.8"
 CairoMakie = "0.15.13"
 DataFrames = "1.8.2"
+Interpolations = "0.15, 0.16"
 JET = "0.9, 0.10, 0.11"
 ModelingToolkit = "11.34.1"
diff --git a/test/extensions/test-interpolations.jl b/test/extensions/test-interpolations.jl
new file mode 100644
index 00000000..6ba902be
--- /dev/null
+++ b/test/extensions/test-interpolations.jl
@@ -0,0 +1,181 @@
+"""Interpolated solar position algorithm backed by the Interpolations.jl extension"""
+
+using Interpolations
+using OhMyThreads: DynamicScheduler
+using StructArrays: StructArrays
+using TimeZones: ZonedDateTime, @tz_str
+
+# a prime minute step gives a deterministic sample that never aligns with the
+# interpolation grid
+sample_times(t0, t1; step = Minute(9973)) = collect(t0:step:t1)
+
+wrapdiff(a, b) = abs(mod(a - b + 180, 360) - 180)
+
+const INTERP_OBSERVERS = [
+    Observer(lat, lon; altitude = alt) for (lat, lon, alt) in [
+            (-85.0, 170.0, 0.0),
+            (-66.5, -120.0, 500.0),
+            (-45.0, 30.0, 0.0),
+            (-23.4, -60.0, 100.0),
+            (0.0, 5.0, 0.0),
+            (23.4, 100.0, 2000.0),
+            (45.0, -90.0, 0.0),
+            (52.35888, 4.88185, 100.0),
+            (66.5, 25.0, 300.0),
+            (85.0, -45.0, 0.0),
+        ]
+]
+
+@testset "Interpolated" begin
+    span1y = (DateTime(2024, 1, 1), DateTime(2025, 1, 1))
+    alg = Interpolated(SPA(); tspan = span1y)
+
+    @testset "Constructor validation" begin
+        @test_throws ArgumentError Interpolated(SPA(); tspan = span1y, out_of_range = :clamp)
+        @test_throws ArgumentError Interpolated(SPA(); tspan = (span1y[2], span1y[1]))
+        @test_throws ArgumentError Interpolated(SPA(); tspan = span1y, step = Day(60))
+        @test_throws ArgumentError Interpolated(SPA(); tspan = span1y, step = Hour(0))
+        @test alg isa Interpolated{SPA}
+        @test occursin("out_of_range = :error", string(alg))
+    end
+
+    @testset "Accuracy against direct SPA" begin
+        for (tspan, times) in [
+                (span1y, sample_times(span1y...)),
+                (
+                    (DateTime(2020, 1, 1), DateTime(2030, 1, 1)),
+                    sample_times(DateTime(2020, 1, 1), DateTime(2030, 1, 1); step = Hour(1289)),
+                ),
+            ]
+            interp = tspan == span1y ? alg : Interpolated(SPA(); tspan)
+            maxaz = 0.0
+            maxel = 0.0
+            for obs in INTERP_OBSERVERS, dt in times
+                p1 = solar_position(obs, dt, interp, NoRefraction())
+                p2 = solar_position(obs, dt, SPA(), NoRefraction())
+                maxaz = max(maxaz, wrapdiff(p1.azimuth, p2.azimuth))
+                maxel = max(maxel, abs(p1.elevation - p2.elevation))
+            end
+            @info "Interpolated vs SPA over $(tspan[1]) to $(tspan[2])" maxaz maxel
+            @test maxaz < 1.0e-6
+            @test maxel < 1.0e-6
+        end
+
+        # queries at the exact span endpoints stay in range
+        for dt in span1y
+            @test solar_position(INTERP_OBSERVERS[1], dt, alg, NoRefraction()) isa SolPos{Float64}
+        end
+    end
+
+    @testset "Refraction parity with SPA" begin
+        obs = INTERP_OBSERVERS[8]
+        dt = DateTime(2024, 6, 21, 12, 30)
+        p1 = solar_position(obs, dt, alg)
+        p2 = solar_position(obs, dt, SPA())
+        @test p1 isa ApparentSolPos{Float64}
+        for field in propertynames(p1)
+            @test getproperty(p1, field) ≈ getproperty(p2, field) atol = 1.0e-6
+        end
+        @test solar_position(obs, dt, alg, NoRefraction()) isa SolPos{Float64}
+    end
+
+    @testset "Out of range behaviour" begin
+        obs = INTERP_OBSERVERS[8]
+        outside = DateTime(2025, 6, 1)
+        err = try
+            solar_position(obs, outside, alg, NoRefraction())
+            nothing
+        catch e
+            e
+        end
+        @test err isa ArgumentError
+        @test occursin("2024-01-01", err.msg)
+        @test occursin("2025-01-01", err.msg)
+        @test_throws ArgumentError solar_position(obs, span1y[1] - Second(1), alg, NoRefraction())
+
+        alg_fb = Interpolated(SPA(); tspan = span1y, out_of_range = :fallback)
+        pfb = solar_position(obs, outside, alg_fb, NoRefraction())
+        pspa = solar_position(obs, outside, SPA(), NoRefraction())
+        # the fallback runs the identical SPA code path, so results are bit equal
+        @test pfb.azimuth === pspa.azimuth
+        @test pfb.elevation === pspa.elevation
+    end
+
+    @testset "Batch, in place, and table paths" begin
+        obs = INTERP_OBSERVERS[8]
+        dts = collect(DateTime(2024, 6, 21):Minute(5):DateTime(2024, 6, 22))
+
+        batch = solar_position(obs, dts, alg, NoRefraction())
+        @test batch isa StructArrays.StructVector{SolPos{Float64}}
+        @test all(
+            batch[i].elevation === solar_position(obs, dts[i], alg, NoRefraction()).elevation
+                for i in eachindex(dts)
+        )
+
+        pos = StructArrays.StructVector{SolPos{Float64}}(undef, length(dts))
+        solar_position!(pos, obs, dts, alg, NoRefraction())
+        @test pos.azimuth == batch.azimuth
+
+        df = DataFrame(datetime = dts)
+        solar_position!(df, obs, alg, NoRefraction())
+        @test df.elevation == batch.elevation
+
+        threaded = solar_position(obs, dts, alg, NoRefraction(), DynamicScheduler())
+        @test threaded.azimuth == batch.azimuth
+    end
+
+    @testset "ZonedDateTime support" begin
+        obs = INTERP_OBSERVERS[8]
+        offset = tz"UTC+02"
+        algz = Interpolated(
+            SPA();
+            tspan = (
+                ZonedDateTime(DateTime(2024, 1, 1, 2), offset),
+                ZonedDateTime(DateTime(2024, 2, 1, 2), offset),
+            ),
+        )
+        @test algz.tspan == (DateTime(2024, 1, 1), DateTime(2024, 2, 1))
+        zdt = ZonedDateTime(DateTime(2024, 1, 15, 14, 30), offset)
+        pz = solar_position(obs, zdt, algz, NoRefraction())
+        pu = solar_position(obs, DateTime(2024, 1, 15, 12, 30), algz, NoRefraction())
+        @test pz.azimuth === pu.azimuth
+    end
+
+    @testset "solar_rate" begin
+        fd_rate = function (obs, dt)
+            p1 = solar_position(obs, dt - Second(1), SPA(), NoRefraction())
+            p2 = solar_position(obs, dt + Second(1), SPA(), NoRefraction())
+            daz = mod(p2.azimuth - p1.azimuth + 180, 360) - 180
+            return (daz * 1800, (p2.elevation - p1.elevation) * 1800)
+        end
+        for obs in INTERP_OBSERVERS
+            for dt in (DateTime(2024, 3, 21, 9), DateTime(2024, 6, 21, 12, 30), DateTime(2024, 12, 21, 15))
+                r = solar_rate(obs, dt, alg)
+                (fd_az, fd_el) = fd_rate(obs, dt)
+                @test r.dazimuth_dt ≈ fd_az atol = 0.01
+                @test r.delevation_dt ≈ fd_el atol = 0.01
+            end
+        end
+
+        # around noon in summer at mid northern latitude the azimuth advances faster
+        # than the mean 15 degrees per hour
+        r = solar_rate(INTERP_OBSERVERS[8], DateTime(2024, 6, 21, 12, 30), alg)
+        @test 15 < r.dazimuth_dt < 40
+
+        zdt = ZonedDateTime(DateTime(2024, 6, 21, 14, 30), tz"UTC+02")
+        rz = solar_rate(INTERP_OBSERVERS[8], zdt, alg)
+        @test rz == solar_rate(INTERP_OBSERVERS[8], DateTime(2024, 6, 21, 12, 30), alg)
+    end
+
+    @testset "Allocations and precision" begin
+        obs = INTERP_OBSERVERS[8]
+        dt = DateTime(2024, 6, 21, 12, 30)
+        measure = (obs, dt, alg) -> @allocated solar_position(obs, dt, alg, NoRefraction())
+        measure(obs, dt, alg)
+        @test measure(obs, dt, alg) == 0
+        @test @inferred(solar_position(obs, dt, alg, NoRefraction())) isa SolPos{Float64}
+
+        p32 = solar_position(Observer{Float32}(52.0, 4.9), dt, alg, NoRefraction())
+        @test p32 isa SolPos{Float32}
+    end
+end

From a201075509a971539bf3969065bf084dcc650aa5 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <s.e.d.lange@tue.nl>
Date: Tue, 28 Jul 2026 12:05:34 +0200
Subject: [PATCH 03/10] Document the Interpolated algorithm

Adds a guide with executed construction, accuracy, benchmark, and
solar_rate examples, registers it in the docs pages, and lists the
wrapper on the positioning page with docstring entries.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Wg3j4XXUooMvDgUpRajxTq
---
 docs/Project.toml                |   1 +
 docs/make.jl                     |   1 +
 docs/src/guides/interpolation.md | 107 +++++++++++++++++++++++++++++++
 docs/src/positioning.md          |  18 ++++++
 4 files changed, 127 insertions(+)
 create mode 100644 docs/src/guides/interpolation.md

diff --git a/docs/Project.toml b/docs/Project.toml
index 0d02934f..5e6a43ba 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -8,6 +8,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OhMyThreads = "67456a42-1dca-4109-a031-0a68de7e3ad5"
diff --git a/docs/make.jl b/docs/make.jl
index 17531e53..84cd9d15 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -43,6 +43,7 @@ makedocs(;
             "guides/autodiff.md",
             "guides/plotting.md",
             "guides/parallel.md",
+            "guides/interpolation.md",
             "guides/modelingtoolkit.md",
             "guides/benchmarking.md",
             "guides/new-algorithm.md",
diff --git a/docs/src/guides/interpolation.md b/docs/src/guides/interpolation.md
new file mode 100644
index 00000000..d8309260
--- /dev/null
+++ b/docs/src/guides/interpolation.md
@@ -0,0 +1,107 @@
+# [Interpolated Solar Position](@id interpolated-position)
+
+Dense time series are the common case in solar energy work: a year of positions at
+one minute resolution is over half a million queries, and [`SPA`](@ref) spends about
+2 µs on each. The [`Interpolated`](@ref SolarPosition.Positioning.Interpolated)
+algorithm removes almost all of that cost by precomputing the slow part once.
+
+The sun's geocentric coordinates change smoothly on annual and monthly timescales, so
+`Interpolated` samples them on a uniform grid and fits cubic B-splines. Everything
+fast or observer dependent, sidereal time, the hour angle, parallax, and the
+conversion to azimuth and elevation, stays closed form and runs through the exact
+same code as `SPA` itself. The interpolation error is around 1e-10 degrees at the
+default one hour grid, about seven orders of magnitude below the accuracy of SPA.
+
+Construction requires [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl),
+which is a weak dependency, so load it first:
+
+```@example interp
+using SolarPosition
+using Dates
+using Interpolations
+
+alg = Interpolated(SPA(); tspan = (DateTime(2024, 1, 1), DateTime(2025, 1, 1)))
+```
+
+The result is a drop in replacement for any other algorithm:
+
+```@example interp
+obs = Observer(52.35888, 4.88185; altitude = 100.0)
+dt = DateTime(2024, 6, 21, 12, 30)
+
+solar_position(obs, dt, alg)
+```
+
+## Accuracy
+
+The interpolant reproduces the wrapped algorithm to well below its own accuracy. One
+day of positions at minute resolution against direct `SPA`:
+
+```@example interp
+times = collect(DateTime(2024, 6, 21):Minute(1):DateTime(2024, 6, 22))
+exact = solar_position(obs, times, SPA(), NoRefraction())
+fast = solar_position(obs, times, alg, NoRefraction())
+
+maximum(abs.(fast.elevation .- exact.elevation))
+```
+
+## Speed
+
+Construction samples the geocentric state of `SPA` about 8800 times for a one year
+span at the default `step = Hour(1)`, a few milliseconds of work that is threaded
+over the available Julia threads. Each query afterwards costs a few spline
+evaluations plus the closed form reconstruction:
+
+```@example interp
+using BenchmarkTools
+
+@btime solar_position($obs, $dt, $(SPA()), $(NoRefraction()));
+@btime solar_position($obs, $dt, $alg, $(NoRefraction()));
+nothing # hide
+```
+
+On the machine that built these docs this is roughly a 10x speedup per query, so the
+construction cost is repaid after a few thousand queries. Below that, direct `SPA` is
+the better tool. Two properties make the interpolant attractive beyond raw speed:
+
+- It is observer independent. The splines capture the sun as seen from the Earth's
+  center, so one instance serves every site in a simulation grid.
+- It is immutable, so evaluation is thread safe by construction and composes with the
+  [OhMyThreads extension](@ref parallel-computing).
+
+## Rate of change
+
+Because evaluation is cheap, derivatives come almost for free.
+[`solar_rate`](@ref SolarPosition.Positioning.solar_rate) returns the rate of change
+of azimuth and elevation in degrees per hour, which is what tracker control loops and
+slew rate limits need:
+
+```@example interp
+solar_rate(obs, dt, alg)
+```
+
+## Out of range queries
+
+Queries outside `tspan` throw by default, because silently falling back to the exact
+algorithm would be a hard to notice 10x slowdown:
+
+```@example interp
+try
+    solar_position(obs, DateTime(2026, 1, 1), alg)
+catch err
+    println(err.msg)
+end
+```
+
+Pass `out_of_range = :fallback` to get the wrapped algorithm outside the span
+instead. This is convenient when a handful of stray timestamps should not fail a
+whole pipeline:
+
+```@example interp
+alg_fb = Interpolated(
+    SPA();
+    tspan = (DateTime(2024, 1, 1), DateTime(2025, 1, 1)),
+    out_of_range = :fallback,
+)
+solar_position(obs, DateTime(2026, 1, 1), alg_fb).elevation
+```
diff --git a/docs/src/positioning.md b/docs/src/positioning.md
index e72fc498..98f147f9 100644
--- a/docs/src/positioning.md
+++ b/docs/src/positioning.md
@@ -64,6 +64,11 @@ The following solar position algorithms are currently implemented in SolarPositi
 | [`USNO`](@ref SolarPosition.Positioning.USNO)         | [USNO](@cite)   | ±0.0500° | None                                             | ✅     |
 | [`SPA`](@ref SolarPosition.Positioning.SPA)           | [RA04](@cite)   | ±0.0003° | Built-in                                         | ✅     |
 
+In addition, [`Interpolated`](@ref SolarPosition.Positioning.Interpolated) wraps `SPA`
+with a precomputed spline of its geocentric coordinates for roughly 10x faster
+repeated evaluation at matching accuracy. It is a wrapper rather than a standalone
+algorithm, see the [Interpolated Solar Position](@ref interpolated-position) guide.
+
 ## [PSA](@id psa-algorithm)
 
 The PSA (Plataforma Solar de Almería) algorithm is the default high-accuracy solar
@@ -121,3 +126,16 @@ position calculation with periodic terms for Earth heliocentric longitude and la
 ```@docs
 SolarPosition.Positioning.SPA
 ```
+
+## [Interpolated](@id interpolated-algorithm)
+
+The Interpolated wrapper precomputes cubic B-splines of the geocentric solar
+coordinates of a wrapped exact algorithm over a fixed time span and reconstructs
+topocentric positions analytically at query time. See the
+[Interpolated Solar Position](@ref interpolated-position) guide for usage, accuracy,
+and benchmarks.
+
+```@docs
+SolarPosition.Positioning.Interpolated
+SolarPosition.Positioning.solar_rate
+```

From 7a0dccbff9796e3f0e9c444c122d417736f82bfd Mon Sep 17 00:00:00 2001
From: Stefan de Lange <s.e.d.lange@tue.nl>
Date: Tue, 28 Jul 2026 12:33:51 +0200
Subject: [PATCH 04/10] Cover Interpolated in the ModelingToolkit tests and
 guide

The block test asserts the interpolant is a drop in replacement for the
wrapped SPA with both NoRefraction and DefaultRefraction, and the guide
gains a section on using Interpolated as a high accuracy forcing,
including span sizing and the :fallback recommendation inside solvers.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Wg3j4XXUooMvDgUpRajxTq
---
 docs/src/guides/modelingtoolkit.md | 47 ++++++++++++++++++++++++++++++
 test/extensions/test-mtk.jl        | 34 ++++++++++++++++++++-
 2 files changed, 80 insertions(+), 1 deletion(-)

diff --git a/docs/src/guides/modelingtoolkit.md b/docs/src/guides/modelingtoolkit.md
index cccf29b1..3dcab97f 100644
--- a/docs/src/guides/modelingtoolkit.md
+++ b/docs/src/guides/modelingtoolkit.md
@@ -201,6 +201,53 @@ lines!(ax3, sol.t ./ 3600, sol[sys_building.sun.elevation])
 fig
 ```
 
+## High Accuracy Forcing with Interpolated
+
+The default `PSA()` is fast but carries its ±0.0083° accuracy. Passing `SPA()` gives
+±0.0003° at about 2 µs per evaluation, which the solver pays at every stage of every
+step. The [`Interpolated`](@ref SolarPosition.Positioning.Interpolated) wrapper keeps
+SPA accuracy at close to PSA cost, which makes it the right choice when a model needs
+the best available forcing. Load Interpolations.jl, build the interpolant to cover the
+simulation window with some margin, and pass it like any other algorithm:
+
+```@example mtk
+using Interpolations
+
+t0 = DateTime(2024, 6, 21, 0, 0, 0)
+interp = Interpolated(
+    SPA();
+    tspan = (t0 - Day(1), t0 + Day(2)),
+    out_of_range = :fallback,
+)
+
+@named sun = SolarPositionBlock()
+sys = mtkcompile(sun)
+
+pmap = [
+    sys.observer => Observer(37.7749, -122.4194, 100.0),
+    sys.t0 => t0,
+    sys.algorithm => interp,
+    sys.refraction => NoRefraction(),
+]
+
+prob = ODEProblem(sys, pmap, (0.0, 86400.0))
+sol = solve(prob; saveat = 3600.0)
+
+# identical model with direct SPA for comparison
+pmap_spa = [pmap[1], pmap[2], sys.algorithm => SPA(), pmap[4]]
+sol_spa = solve(ODEProblem(sys, pmap_spa, (0.0, 86400.0)); saveat = 3600.0)
+
+maximum(abs.(sol[sys.elevation] .- sol_spa[sys.elevation]))
+```
+
+Two practical notes. First, size `tspan` to cover the whole solve measured from `t0`
+and pad it generously, since construction costs milliseconds and a few hundred
+kilobytes per year. Second, `out_of_range = :fallback` is a good idea inside a solver,
+because a stray evaluation outside the span then degrades to exact SPA instead of
+aborting the integration. The interpolant is observer independent and immutable, so
+one instance can be shared by every `SolarPositionBlock` in a model and across
+threads.
+
 ## Implementation Details
 
 The extension works by registering the [`solar_position`](@ref) function and helper functions as
diff --git a/test/extensions/test-mtk.jl b/test/extensions/test-mtk.jl
index 71527537..0dfb1de0 100644
--- a/test/extensions/test-mtk.jl
+++ b/test/extensions/test-mtk.jl
@@ -11,9 +11,11 @@ using SolarPosition:
     SolarPositionBlock
 
 using SolarPosition: HUGHES, BENNETT, ARCHER, MICHALSKY, SG2
+using SolarPosition: Interpolated
+using Interpolations
 using ModelingToolkit: @named, @variables, @parameters, unknowns, System, mtkcompile
 using ModelingToolkit: t_nounits as t, D_nounits as D
-using Dates: DateTime
+using Dates: Dates, DateTime
 using OrdinaryDiffEq
 using CairoMakie
 
@@ -213,4 +215,34 @@ using CairoMakie
             end
         end
     end
+
+    @testset "Interpolated algorithm in the block" begin
+        # the interpolant is a drop in replacement for the wrapped SPA, so the block
+        # must produce the same outputs with either within interpolation error
+        t0_interp = DateTime(2024, 6, 21)
+        interp = Interpolated(
+            SPA();
+            tspan = (t0_interp - Dates.Day(1), t0_interp + Dates.Day(2)),
+        )
+        obs_interp = Observer(52.35888, 4.88185, 100.0)
+
+        @named sun = SolarPositionBlock()
+        sys = mtkcompile(sun)
+
+        for refr in (NoRefraction(), DefaultRefraction())
+            sols = map((interp, SPA())) do alg
+                pmap = [
+                    sys.observer => obs_interp,
+                    sys.t0 => t0_interp,
+                    sys.algorithm => alg,
+                    sys.refraction => refr,
+                ]
+                prob = ODEProblem(sys, pmap, (0.0, 86400.0))
+                solve(prob; saveat = 3600.0)
+            end
+            for output in (sys.azimuth, sys.elevation, sys.zenith)
+                @test maximum(abs.(sols[1][output] .- sols[2][output])) < 1.0e-8
+            end
+        end
+    end
 end

From 1b987f21bbc5a31b22844d37147f0a19ca52b4ed Mon Sep 17 00:00:00 2001
From: Stefan de Lange <s.e.d.lange@tue.nl>
Date: Tue, 28 Jul 2026 12:55:58 +0200
Subject: [PATCH 05/10] Warn when the sampling step costs more accuracy than
 SPA offers

Measured against direct SPA the interpolation error crosses SPA's own
3e-4 degree accuracy between 3 and 5 day steps, so steps above 3 days
now warn at construction. The docstring lists the measured error at
several steps and the 30 day hard cap stays as the unwrap safety limit.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Wg3j4XXUooMvDgUpRajxTq
---
 src/Positioning/interpolated.jl        | 12 +++++++++++-
 test/extensions/test-interpolations.jl |  4 ++++
 2 files changed, 15 insertions(+), 1 deletion(-)

diff --git a/src/Positioning/interpolated.jl b/src/Positioning/interpolated.jl
index c635bf21..e865be71 100644
--- a/src/Positioning/interpolated.jl
+++ b/src/Positioning/interpolated.jl
@@ -30,7 +30,11 @@ Construction requires Interpolations.jl to be loaded:
 - `tspan::Tuple`: the valid query span as a pair of `DateTime` or `ZonedDateTime`.
   Zoned times are converted to UTC.
 - `step::Period = Hour(1)`: the sampling grid spacing. Must be positive and at most 30
-  days so that right ascension advances much less than half a turn per step.
+  days so that right ascension advances much less than half a turn per step. Steps
+  above 3 days warn, because the interpolation error then exceeds the wrapped
+  algorithm's own accuracy: measured against direct SPA the maximum error is about
+  1e-10 degrees at `Hour(1)`, 4e-7 at `Day(1)`, 2e-4 at `Day(3)`, and 3e-2 at
+  `Day(30)`.
 - `out_of_range::Symbol = :error`: behaviour for queries outside `tspan`. `:error`
   throws an `ArgumentError`, `:fallback` silently calls the wrapped exact algorithm.
 
@@ -77,6 +81,12 @@ function Interpolated(
     # which the unwrap during construction relies on
     Dates.Millisecond(0) < stepms <= Dates.Millisecond(Dates.Day(30)) ||
         throw(ArgumentError("step must be positive and at most 30 days, got $step"))
+    # beyond 3 days the 13.7 day nutation ripple is sampled too coarsely and the
+    # interpolation error exceeds the wrapped algorithm's own accuracy
+    stepms > Dates.Millisecond(Dates.Day(3)) && @warn(
+        "Interpolated with step = $step has an interpolation error above the wrapped " *
+            "algorithm's own accuracy. Use a step of 3 days or less to stay below it.",
+    )
     (itp_α, itp_δ, itp_R, itp_eqeq) = _build_interpolants(algorithm, t0, t1, stepms)
     return Interpolated(
         algorithm, itp_α, itp_δ, itp_R, itp_eqeq, (t0, t1), stepms,
diff --git a/test/extensions/test-interpolations.jl b/test/extensions/test-interpolations.jl
index 6ba902be..ff32f3c7 100644
--- a/test/extensions/test-interpolations.jl
+++ b/test/extensions/test-interpolations.jl
@@ -35,6 +35,10 @@ const INTERP_OBSERVERS = [
         @test_throws ArgumentError Interpolated(SPA(); tspan = (span1y[2], span1y[1]))
         @test_throws ArgumentError Interpolated(SPA(); tspan = span1y, step = Day(60))
         @test_throws ArgumentError Interpolated(SPA(); tspan = span1y, step = Hour(0))
+        @test_logs (:warn, r"interpolation error above") Interpolated(
+            SPA(); tspan = span1y, step = Day(10),
+        )
+        @test_logs Interpolated(SPA(); tspan = span1y, step = Day(3))
         @test alg isa Interpolated{SPA}
         @test occursin("out_of_range = :error", string(alg))
     end

From bcb95c4a47ee85c2543dd688d5c3ccffd1e491e1 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <s.e.d.lange@tue.nl>
Date: Tue, 28 Jul 2026 13:06:24 +0200
Subject: [PATCH 06/10] Canonicalize the step in show and close the coverage
 gaps

The step now prints as 1 hour instead of 3600000 milliseconds. New
tests hit the extension missing stub through a non SPA algorithm and
the ApparentSolPos result type through the default refraction batch
path, which were the three uncovered lines.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Wg3j4XXUooMvDgUpRajxTq
---
 src/Positioning/interpolated.jl        |  3 ++-
 test/extensions/test-interpolations.jl | 11 +++++++++++
 2 files changed, 13 insertions(+), 1 deletion(-)

diff --git a/src/Positioning/interpolated.jl b/src/Positioning/interpolated.jl
index e865be71..93d4ce1c 100644
--- a/src/Positioning/interpolated.jl
+++ b/src/Positioning/interpolated.jl
@@ -113,7 +113,8 @@ end
 function Base.show(io::IO, alg::Interpolated)
     print(
         io, "Interpolated(", alg.algorithm, "; tspan = ", alg.tspan,
-        ", step = ", alg.step, ", out_of_range = :", alg.out_of_range, ")",
+        ", step = ", Dates.canonicalize(alg.step),
+        ", out_of_range = :", alg.out_of_range, ")",
     )
     return nothing
 end
diff --git a/test/extensions/test-interpolations.jl b/test/extensions/test-interpolations.jl
index ff32f3c7..544a0a76 100644
--- a/test/extensions/test-interpolations.jl
+++ b/test/extensions/test-interpolations.jl
@@ -41,6 +41,13 @@ const INTERP_OBSERVERS = [
         @test_logs Interpolated(SPA(); tspan = span1y, step = Day(3))
         @test alg isa Interpolated{SPA}
         @test occursin("out_of_range = :error", string(alg))
+        @test occursin("step = 1 hour", string(alg))
+
+        # the extension only implements sampling for SPA, so any other algorithm
+        # reaches the stub that asks for Interpolations.jl
+        @test_throws ArgumentError SolarPosition.Positioning._build_interpolants(
+            PSA(), span1y[1], span1y[2], Millisecond(Hour(1)),
+        )
     end
 
     @testset "Accuracy against direct SPA" begin
@@ -111,6 +118,10 @@ const INTERP_OBSERVERS = [
 
         batch = solar_position(obs, dts, alg, NoRefraction())
         @test batch isa StructArrays.StructVector{SolPos{Float64}}
+
+        # the default refraction batch path resolves to ApparentSolPos like SPA
+        batch_ref = solar_position(obs, dts, alg)
+        @test batch_ref isa StructArrays.StructVector{ApparentSolPos{Float64}}
         @test all(
             batch[i].elevation === solar_position(obs, dts[i], alg, NoRefraction()).elevation
                 for i in eachindex(dts)

From 7af6382327b4a5c19b1de99b365033ac197b2d81 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <s.e.d.lange@tue.nl>
Date: Tue, 28 Jul 2026 13:18:03 +0200
Subject: [PATCH 07/10] Link first mentions with @ref and slim the README
 footnote

Adds cross references at the first prose mention of Observer, the five
algorithms, DefaultRefraction, and PSA and SPA in the new MTK section.
The README Quadmath footnote becomes a one line pointer to the
precision guide, which keeps the full explanation.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Wg3j4XXUooMvDgUpRajxTq
---
 README.md                          |  9 ++++-----
 docs/src/guides/getting-started.md |  2 +-
 docs/src/guides/modelingtoolkit.md |  3 ++-
 docs/src/guides/new-algorithm.md   |  2 +-
 docs/src/guides/precision.md       | 10 +++++++---
 5 files changed, 15 insertions(+), 11 deletions(-)

diff --git a/README.md b/README.md
index 9711df8d..0bd12fdd 100644
--- a/README.md
+++ b/README.md
@@ -117,14 +117,13 @@ latitudes from 70°N to 60°S:
 | Algorithm | `Float32` error | `Float32` runtime | `Float64` error | `Float128` error | `Float128` runtime | `BigFloat` runtime |
 | --------- | --------------- | ----------------- | --------------- | ---------------- | ------------------ | ------------------ |
 | PSA       | 0.011°          | 1.3× faster       | 1.3e-11°        | 7.1e-30°         | 93× slower         | 500× slower        |
-| NOAA      | 0.0019°         | 1.2× faster       | 3.5e-12°        | broken[^1]       | n/a                | 330× slower        |
+| NOAA      | 0.0019°         | 1.2× faster       | 3.5e-12°        | broken           | n/a                | 330× slower        |
 | Walraven  | 0.0040°         | 1.3× faster       | 6.4e-12°        | 8.6e-30°         | 74× slower         | 400× slower        |
-| USNO      | 0.0036°         | 1.2× faster       | 9.5e-12°        | broken[^1]       | n/a                | 300× slower        |
+| USNO      | 0.0036°         | 1.2× faster       | 9.5e-12°        | broken           | n/a                | 300× slower        |
 | SPA       | 0.012°          | 1.5× faster       | 1.8e-11°        | 1.1e-29°         | 118× slower        | 330× slower        |
 
-[^1]: Quadmath.jl v1.0.1 implements `rem` with round-to-nearest instead of truncated
-    semantics, which breaks the degree reduction in Base's `sind` and `cosd`, so `NOAA`
-    and `USNO` give wrong results at `Float128` until that is fixed upstream.
+See the [precision guide](https://juliaastro.org/SolarPosition/stable/guides/precision/)
+for details, including why `NOAA` and `USNO` are currently broken at `Float128`.
 
 ## Refraction correction algorithms
 
diff --git a/docs/src/guides/getting-started.md b/docs/src/guides/getting-started.md
index dd9f26e9..595eb905 100644
--- a/docs/src/guides/getting-started.md
+++ b/docs/src/guides/getting-started.md
@@ -28,7 +28,7 @@ using DataFrames
 ## Defining a location
 
 We can observe the sun from anywhere on earth. To define an observer location, we use
-the `Observer` struct, which takes latitude, longitude, and optionally altitude
+the [`Observer`](@ref SolarPosition.Positioning.Observer) struct, which takes latitude, longitude, and optionally altitude
 (in meters) as arguments.
 
 ```@example getting-started
diff --git a/docs/src/guides/modelingtoolkit.md b/docs/src/guides/modelingtoolkit.md
index 3dcab97f..31e12363 100644
--- a/docs/src/guides/modelingtoolkit.md
+++ b/docs/src/guides/modelingtoolkit.md
@@ -203,7 +203,8 @@ fig
 
 ## High Accuracy Forcing with Interpolated
 
-The default `PSA()` is fast but carries its ±0.0083° accuracy. Passing `SPA()` gives
+The default [`PSA`](@ref SolarPosition.Positioning.PSA) is fast but carries its
+±0.0083° accuracy. Passing [`SPA`](@ref SolarPosition.Positioning.SPA) gives
 ±0.0003° at about 2 µs per evaluation, which the solver pays at every stage of every
 step. The [`Interpolated`](@ref SolarPosition.Positioning.Interpolated) wrapper keeps
 SPA accuracy at close to PSA cost, which makes it the right choice when a model needs
diff --git a/docs/src/guides/new-algorithm.md b/docs/src/guides/new-algorithm.md
index 60bf5da4..667b7f07 100644
--- a/docs/src/guides/new-algorithm.md
+++ b/docs/src/guides/new-algorithm.md
@@ -11,7 +11,7 @@ Adding a new algorithm involves these steps:
 
 1. [**Create the algorithm struct**](@ref step-1-create-struct) - Define a type that subtypes [`SolarAlgorithm`](@ref SolarPosition.Positioning.SolarAlgorithm).
 2. [**Implement the core function**](@ref step-2-implement-core) - Write `_solar_position` for your algorithm.
-3. [**Handle refraction**](@ref step-3-handle-refraction) - Define how your algorithm interacts with `DefaultRefraction`.
+3. [**Handle refraction**](@ref step-3-handle-refraction) - Define how your algorithm interacts with [`DefaultRefraction`](@ref SolarPosition.Refraction.DefaultRefraction).
 4. [**Export the algorithm**](@ref step-4-export) - Make it available to users.
 5. [**Write tests**](@ref step-5-write-tests) - Validate correctness against reference values.
 6. [**Document**](@ref step-6-document) - Add docstrings and update documentation.
diff --git a/docs/src/guides/precision.md b/docs/src/guides/precision.md
index 010d8da2..5ce43825 100644
--- a/docs/src/guides/precision.md
+++ b/docs/src/guides/precision.md
@@ -32,7 +32,9 @@ precision types never ride on the ~2.45e6 Julian Date.
 - `Float32` trades a little accuracy for a modest speedup. The error stays within each
   algorithm's own claimed accuracy.
 - `Float128` from [Quadmath.jl](https://github.com/JuliaMath/Quadmath.jl) gives quad
-  precision of roughly 1e-30 degrees for `PSA`, `SPA`, and `Walraven` at a large
+  precision of roughly 1e-30 degrees for [`PSA`](@ref SolarPosition.Positioning.PSA),
+  [`SPA`](@ref SolarPosition.Positioning.SPA), and
+  [`Walraven`](@ref SolarPosition.Positioning.Walraven) at a large
   runtime cost.
 - `BigFloat` gives arbitrary precision. Raise it with `setprecision(BigFloat, bits)`.
   This is the right tool for generating reference values.
@@ -58,8 +60,10 @@ over a grid of 125 combinations of dates from 2015 to 2035 and latitudes from 70
 | SPA       | 0.012°          | 1.5× faster       | 1.8e-11°        | 1.1e-29°         | 118× slower        | 330× slower        |
 
 [^1]: Quadmath.jl v1.0.1 implements `rem` with round-to-nearest instead of truncated
-    semantics, which breaks the degree reduction in Base's `sind` and `cosd`, so `NOAA`
-    and `USNO` give wrong results at `Float128` until that is fixed upstream.
+    semantics, which breaks the degree reduction in Base's `sind` and `cosd`, so
+    [`NOAA`](@ref SolarPosition.Positioning.NOAA) and
+    [`USNO`](@ref SolarPosition.Positioning.USNO) give wrong results at `Float128`
+    until that is fixed upstream.
 
 `Float128` is a fixed 113-bit significand type backed by libquadmath, is allocation
 free, and costs roughly 100× `Float64`. `BigFloat` is arbitrary precision backed by

From 3f96d2f1b61d7dfda7f6eddc3fc920d601c45dab Mon Sep 17 00:00:00 2001
From: Stefan de Lange <s.e.d.lange@tue.nl>
Date: Tue, 28 Jul 2026 13:29:06 +0200
Subject: [PATCH 08/10] Document solver interaction in the ModelingToolkit
 guide

Explains that saveat controls recording rather than stepping, shows the
exact post hoc query of observed sun angles through the solution
object, and covers the two error controller failure modes with a
runnable example: sunrise and sunset from transit_sunrise_sunset as
d_discontinuities, and an insolation quadrature state that puts the
forcing in the error budget. Closes with the reference solve check.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Wg3j4XXUooMvDgUpRajxTq
---
 docs/src/guides/modelingtoolkit.md | 91 ++++++++++++++++++++++++++++++
 1 file changed, 91 insertions(+)

diff --git a/docs/src/guides/modelingtoolkit.md b/docs/src/guides/modelingtoolkit.md
index 31e12363..5052791a 100644
--- a/docs/src/guides/modelingtoolkit.md
+++ b/docs/src/guides/modelingtoolkit.md
@@ -249,6 +249,97 @@ aborting the integration. The interpolant is observer independent and immutable,
 one instance can be shared by every `SolarPositionBlock` in a model and across
 threads.
 
+## Working with the Solver
+
+Two properties of solar forcing surprise people the first time: the solver seems to
+skip straight past the day unless `saveat` is given, and the adaptive error control
+seems unaware of the sun. Both have clean solutions.
+
+### Sampling outputs without saveat
+
+`saveat` does not make the solver take more steps. Save points are filled in from the
+solution's dense interpolant, so it only controls what gets recorded. For the bare
+`SolarPositionBlock` the compiled system has no differential states at all, so the
+solver correctly jumps from start to end in one step, and without `saveat` the
+solution object holds just the two endpoints.
+
+```@example mtk
+obs = Observer(52.35888, 4.88185, 100.0)
+t0 = DateTime(2024, 6, 21, 0, 0, 0)
+
+@named sun = SolarPositionBlock()
+sys = mtkcompile(sun)
+pmap = [
+    sys.observer => obs,
+    sys.t0 => t0,
+    sys.algorithm => PSA(),
+    sys.refraction => NoRefraction(),
+]
+sol = solve(ODEProblem(sys, pmap, (0.0, 86400.0)))
+length(sol.t)
+```
+
+The better tool is the solution object itself. The solar outputs are observed
+variables that depend only on parameters and time, so querying the solution
+re-evaluates the exact solar position at any requested time, at any resolution,
+independent of how coarsely the solver stepped:
+
+```@example mtk
+sol(0.0:600.0:86400.0; idxs = sys.elevation)
+```
+
+This is exact for the sun angles. For observed variables that also involve states the
+query uses the state interpolant, whose accuracy is set by the solver tolerances.
+
+### Making the error controller see the forcing
+
+The embedded error estimator only controls the error of integrating the states it is
+given. Two distinct failure modes follow, each with its own fix.
+
+The first is nonsmoothness. Solar forcing models clip at the horizon, typically with
+`max(0, ...)`, and a step that spans sunrise or sunset sees a kink, rejects, and
+thrashes. The fix is to tell the solver where the kinks are.
+[`transit_sunrise_sunset`](@ref) computes them, and `d_discontinuities` passes them
+in, converted to simulation seconds.
+
+The second is smooth blindness. A state with a large time constant filters the
+forcing, so the controller sees little state error and takes steps that under resolve
+the forcing's integral. The fix is to add the integral as a state, here `E_sol`, so
+the quadrature of the forcing enters the error budget directly:
+
+```@example mtk
+@parameters C = 5.0e5 k = 25.0
+@variables T_room(t) = 18.0 E_sol(t) = 0.0 Q(t)
+
+eqs = [
+    Q ~ 800 * max(0, sind(sun.elevation)),
+    D(T_room) ~ (Q - k * (T_room - 15.0)) / C,
+    D(E_sol) ~ Q,
+]
+@named house = System(eqs, t; systems = [sun])
+sys = mtkcompile(house)
+
+pmap = [
+    sys.sun.observer => obs,
+    sys.sun.t0 => t0,
+    sys.sun.algorithm => PSA(),
+    sys.sun.refraction => NoRefraction(),
+]
+prob = ODEProblem(sys, pmap, (0.0, 86400.0))
+
+events = transit_sunrise_sunset(obs, t0)
+kinks = [Dates.value(dt - t0) / 1000 for dt in (events.sunrise, events.sunset)]
+
+sol = solve(prob; d_discontinuities = kinks, reltol = 1.0e-8)
+(steps = length(sol.t), daily_insolation = sol[sys.E_sol][end])
+```
+
+A vector `abstol` matched to `unknowns(sys)` gives the quadrature state a tolerance
+in its own physical units when it should not share the default. Whatever combination
+you settle on, verify it once against a reference solve at `reltol = 1e-10` and
+compare the quantities you care about. That check, not the step count, is what shows
+the recipe is sufficient.
+
 ## Implementation Details
 
 The extension works by registering the [`solar_position`](@ref) function and helper functions as

From 407b1ff9393c0ab3c019bf4ac51751e2dcdcb8d8 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <s.e.d.lange@tue.nl>
Date: Tue, 28 Jul 2026 13:36:39 +0200
Subject: [PATCH 09/10] Point the README precision link at the dev docs

The guide is not in a stable release yet, so the stable URL 404s in the
link checker. The dev URL resolves.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Wg3j4XXUooMvDgUpRajxTq
---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 0bd12fdd..fdae596f 100644
--- a/README.md
+++ b/README.md
@@ -122,7 +122,7 @@ latitudes from 70°N to 60°S:
 | USNO      | 0.0036°         | 1.2× faster       | 9.5e-12°        | broken           | n/a                | 300× slower        |
 | SPA       | 0.012°          | 1.5× faster       | 1.8e-11°        | 1.1e-29°         | 118× slower        | 330× slower        |
 
-See the [precision guide](https://juliaastro.org/SolarPosition/stable/guides/precision/)
+See the [precision guide](https://juliaastro.org/SolarPosition.jl/dev/guides/precision/)
 for details, including why `NOAA` and `USNO` are currently broken at `Float128`.
 
 ## Refraction correction algorithms

From e5eb8381e77eccd81f5832de000bb5655def2374 Mon Sep 17 00:00:00 2001
From: Stefan de Lange <s.e.d.lange@tue.nl>
Date: Tue, 28 Jul 2026 13:48:29 +0200
Subject: [PATCH 10/10] Show DifferentiationInterface composing with the
 generic code

A small autodiff guide section runs the same gradient through the
ForwardDiff and finite difference backends and prints their agreement,
demonstrating that backend agnostic tooling works without any AD hooks
in the package.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Wg3j4XXUooMvDgUpRajxTq
---
 docs/Project.toml           |  2 ++
 docs/src/guides/autodiff.md | 18 ++++++++++++++++++
 2 files changed, 20 insertions(+)

diff --git a/docs/Project.toml b/docs/Project.toml
index 5e6a43ba..6064c94f 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -4,9 +4,11 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
+FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
diff --git a/docs/src/guides/autodiff.md b/docs/src/guides/autodiff.md
index befbd605..3cf8d6d8 100644
--- a/docs/src/guides/autodiff.md
+++ b/docs/src/guides/autodiff.md
@@ -27,6 +27,24 @@ derivatives. Time is a `DateTime`, not a number, so derivatives with respect to
 are not available through this route. Differentiate through a wrapper that maps a
 number to a `DateTime` if you need them.
 
+## DifferentiationInterface
+
+Because the differentiability comes from the code being generic rather than from any
+AD specific hooks, backend agnostic tooling composes with it out of the box.
+[DifferentiationInterface.jl](https://github.com/JuliaDiff/DifferentiationInterface.jl)
+lets you write the gradient call once and swap the backend freely. Here the
+ForwardDiff and finite difference backends agree to the finite difference accuracy:
+
+```@example autodiff
+import DifferentiationInterface as DI
+import FiniteDiff
+
+f(x) = solar_position(Observer(x[1], x[2]), dt, SPA(), NoRefraction()).elevation
+g_ad = DI.gradient(f, DI.AutoForwardDiff(), [45.0, 10.0])
+g_num = DI.gradient(f, DI.AutoFiniteDiff(), [45.0, 10.0])
+(g_ad, maximum(abs.(g_ad .- g_num)))
+```
+
 ## Example: optimizing a solar panel orientation
 
 Gradient ascent on a plane-of-array irradiance proxy finds the best fixed tilt and