Add an Interpolated algorithm for fast repeated evaluation

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"

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
 

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"

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",

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

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

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
+```

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

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.

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