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/README.md b/README.md
index 9711df8d..fdae596f 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.jl/dev/guides/precision/)
+for details, including why `NOAA` and `USNO` are currently broken at `Float128`.
 
 ## Refraction correction algorithms
 
diff --git a/docs/Project.toml b/docs/Project.toml
index 0d02934f..6064c94f 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -4,10 +4,13 @@ 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"
 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/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
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/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/guides/modelingtoolkit.md b/docs/src/guides/modelingtoolkit.md
index cccf29b1..5052791a 100644
--- a/docs/src/guides/modelingtoolkit.md
+++ b/docs/src/guides/modelingtoolkit.md
@@ -201,6 +201,145 @@ lines!(ax3, sol.t ./ 3600, sol[sys_building.sun.elevation])
 fig
 ```
 
+## High Accuracy Forcing with Interpolated
+
+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
+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.
+
+## 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
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
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
+```
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..93d4ce1c
--- /dev/null
+++ b/src/Positioning/interpolated.jl
@@ -0,0 +1,222 @@
+"""
+    $(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. 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.
+
+# 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"))
+    # 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,
+        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 = ", Dates.canonicalize(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/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)
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..544a0a76
--- /dev/null
+++ b/test/extensions/test-interpolations.jl
@@ -0,0 +1,196 @@
+"""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_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))
+        @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
+        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}}
+
+        # 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)
+        )
+
+        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
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