Differentiate with respect to refraction model parameters

.gitignore

@@ -15,6 +15,6 @@ node_modules
 .vscode/settings.json
 **.old
 play/Project.toml
-docs/.CondaPkg/
+**/.CondaPkg/
 *.mem
 **/Manifest-*.toml

docs/src/guides/autodiff.md

@@ -27,6 +27,21 @@ 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.
 
+A refraction model's own parameters are differentiable inputs too. The model types and
+the result types promote, so a dual valued pressure needs no change to the observer:
+
+```@example autodiff
+ForwardDiff.derivative(
+    p -> solar_position(Observer(45.0, 10.0), dt, PSA(), HUGHES(p, 10.0)).apparent_elevation,
+    101325.0,
+)  # degrees of apparent elevation per pascal
+```
+
+[`Interpolated`](@ref SolarPosition.Positioning.Interpolated) differentiates as well. Its
+interpolants cover only the geocentric quantities, which depend on time alone, so the
+duals travel through the topocentric half it shares with the wrapped algorithm and the
+derivatives match that algorithm's to roundoff.
+
 ## DifferentiationInterface
 
 Because the differentiability comes from the code being generic rather than from any

docs/src/guides/modelingtoolkit.md

@@ -296,16 +296,40 @@ query uses the state interpolant, whose accuracy is set by the solver tolerances
 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 first is nonsmoothness. Solar forcing models clip at the horizon with `max(0, ...)`,
+and a step spanning sunrise or sunset hits a kink and rejects. `d_discontinuities` fixes
+that, but only with the times the model actually breaks at.
+[`transit_sunrise_sunset`](@ref) returns the almanac event, when the sun's upper limb
+reaches −0.8333° with refraction allowed for, whereas `max(0, sind(elevation))` breaks
+when the geometric elevation crosses zero. In Amsterdam at the solstice the two are
+seven minutes apart, so the declared discontinuity lands where nothing happens. The
+almanac event brackets the geometric one, so bisect inside it:
 
-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
+events = transit_sunrise_sunset(obs, t0)
+almanac = [Dates.value(dt - t0) / 1000 for dt in (events.sunrise, events.sunset)]
+
+function elevation_root(lo, hi)
+    elevation(x) = solar_position(
+        obs, t0 + Millisecond(round(Int, 1000x)), PSA(), NoRefraction()
+    ).elevation
+    for _ in 1:60
+        mid = (lo + hi) / 2
+        elevation(lo) * elevation(mid) <= 0 ? (hi = mid) : (lo = mid)
+    end
+    return (lo + hi) / 2
+end
+
+kinks = [elevation_root(a - 1800, a + 1800) for a in almanac]
+(kinks .- almanac) ./ 60  # minutes from the almanac event to the model's kink
+```
+
+The second failure mode is smooth blindness. A state with a large time constant filters
+the forcing, so the controller sees little state error and under resolves the forcing's
+integral. Adding that integral as a state, here `E_sol`, puts its quadrature into the
+error budget. A vector `abstol` matched to `unknowns(sys)` then scales the tolerance to
+physical units, since the default 1e-6 on a joule count reaching 2.5e7 is far tighter
+than the problem needs:
 
 ```@example mtk
 @parameters C = 5.0e5 k = 25.0
@@ -327,18 +351,152 @@ pmap = [
 ]
 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)]
+abstol = [isequal(u, E_sol) ? 1.0e-1 : 1.0e-6 for u in unknowns(sys)]
+sol = solve(prob; d_discontinuities = kinks, reltol = 1.0e-5, abstol)
+(steps = length(sol.t), rejected = sol.stats.nreject, daily_insolation = sol[sys.E_sol][end])
+```
+
+Where the solver put its steps says more than a step count does. A loose `reltol` of
+1e-5 keeps a day to a countable number of steps, so each can be drawn as a vertical line
+under the forcing. The first three rows vary only the declared kink times; the fourth
+adds the quadrature state to the best of them:
+
+```@example mtk
+using CairoMakie: Figure, Axis, Label, Point2f, RGBf, Relative, colgap!, colsize!,
+    hidespines!, hidexdecorations!, hideydecorations!, linesegments!, lines!,
+    rowgap!, rowsize!, text!, vlines!
+
+# the same house without the quadrature state, for comparison
+@named house_filtered = System(eqs[1:2], t; systems = [sun])
+sys_f = mtkcompile(house_filtered)
+prob_f = ODEProblem(
+    sys_f,
+    [
+        sys_f.sun.observer => obs,
+        sys_f.sun.t0 => t0,
+        sys_f.sun.algorithm => PSA(),
+        sys_f.sun.refraction => NoRefraction(),
+    ],
+    (0.0, 86400.0),
+)
+
+variants = [
+    ("no kinks declared", solve(prob_f; reltol = 1.0e-5), RGBf(0, 0.447, 0.698)),
+    (
+        "almanac kinks",
+        solve(prob_f; d_discontinuities = almanac, reltol = 1.0e-5),
+        RGBf(0.902, 0.624, 0),
+    ),
+    (
+        "model kinks",
+        solve(prob_f; d_discontinuities = kinks, reltol = 1.0e-5),
+        RGBf(0, 0.62, 0.451),
+    ),
+    ("model kinks + quadrature", sol, RGBf(0.8, 0.475, 0.655)),
+]
+
+ts = range(0, 86400; length = 3000)
+Qs = collect(sol(ts; idxs = sys.Q))
+rug!(ax, y, s, c) = linesegments!(
+    ax, [Point2f(x / 3600, y + dy) for x in s.t for dy in (-0.36, 0.36)];
+    color = c, linewidth = 1.5,
+)
+
+fig = Figure(size = (980, 540))
+zoom = (3.22, 3.62)  # a window around sunrise
+
+ax_day = Axis(fig[1, 1]; ylabel = "Solar heat gain (W)")
+ax_rug = Axis(
+    fig[2, 1];
+    xlabel = "Time of day (hours)", xticks = 0:3:24,
+    yticks = (
+        1:4,
+        [
+            "$n\n$(length(s.t)) steps, $(s.stats.nreject) rejected"
+                for (n, s, _) in reverse(variants)
+        ],
+    ),
+)
+ax_zoom = Axis(fig[1, 2]; title = "sunrise, zoomed", titlesize = 12)
+ax_zoomrug = Axis(fig[2, 2]; xlabel = "Time of day (hours)", xticks = 3.3:0.1:3.6)
+
+for ax in (ax_day, ax_zoom)
+    lines!(ax, ts ./ 3600, Qs; color = :grey25, linewidth = 2)
+    hidexdecorations!(ax; grid = false)
+end
+for ax in (ax_day, ax_rug)
+    vlines!(ax, kinks ./ 3600; color = :grey55, linestyle = :dash, linewidth = 1)
+end
+for ax in (ax_zoom, ax_zoomrug)
+    vlines!(ax, almanac[1] / 3600; color = :grey55, linestyle = :dot, linewidth = 1.5)
+    vlines!(ax, kinks[1] / 3600; color = :grey55, linestyle = :dash, linewidth = 1.5)
+    hideydecorations!(ax; grid = false)
+end
+for (y, (_, s, c)) in zip(4:-1:1, variants)
+    rug!(ax_rug, y, s, c)
+    rug!(ax_zoomrug, y, s, c)
+end
 
-sol = solve(prob; d_discontinuities = kinks, reltol = 1.0e-8)
-(steps = length(sol.t), daily_insolation = sol[sys.E_sol][end])
+text!(
+    ax_zoom, almanac[1] / 3600, 30; text = " almanac\n sunrise",
+    align = (:left, :bottom), fontsize = 10, color = :grey35,
+)
+text!(
+    ax_zoom, kinks[1] / 3600, 100; text = " model\n kink",
+    align = (:left, :bottom), fontsize = 10, color = :grey35,
+)
+
+ax_day.limits = ((0, 24), (-40, 830))
+ax_rug.limits = ((0, 24), (0.4, 4.6))
+ax_zoom.limits = (zoom, (-40, 830))
+ax_zoomrug.limits = (zoom, (0.4, 4.6))
+for ax in (ax_day, ax_rug, ax_zoom, ax_zoomrug)
+    hidespines!(ax, :t, :r)
+end
+
+Label(
+    fig[0, 1:2], "Where the solver steps over one day (reltol = 1e-5)";
+    fontsize = 15, font = :bold, padding = (0, 0, 4, 0),
+)
+rowsize!(fig.layout, 1, Relative(0.36))
+colsize!(fig.layout, 2, Relative(0.26))
+rowgap!(fig.layout, 6)
+colgap!(fig.layout, 14)
+fig
 ```
 
-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.
+The top row declares nothing and steps across the day at near constant spacing,
+indifferent to whether the sun is up. The second declares the almanac times, and the
+zoom shows why that barely helps: the dotted line sits where the forcing is still flat,
+so the step is spent on a discontinuity that is not there and the real kink at the dashed
+line is met unprepared. The third lands a step exactly on the break. The fourth changes
+the sampling very little.
+
+Against a `reltol = 1e-13` reference, where the error columns are the largest room
+temperature deviation over the day and the relative error in daily insolation:
+
+| variant | steps | rejected | `f` evals | max ΔT_room | rel. err. insolation |
+| :--- | ---: | ---: | ---: | ---: | ---: |
+| no kinks declared | 21 | 1 | 129 | 1.7e-2 K | — |
+| almanac kinks | 23 | 0 | 137 | 7.7e-3 K | — |
+| model kinks | 21 | 0 | 125 | 1.4e-4 K | — |
+| model kinks + quadrature | 21 | 0 | 125 | 7.0e-5 K | 4.3e-8 |
+
+Almost all of the benefit comes from one change. Bisecting for the model's own kink times
+buys a factor of 125 on the temperature and costs slightly less than declaring nothing,
+since the rejections it avoids more than pay for the two steps it forces. The almanac
+times alone buy a factor of two.
+
+The quadrature state does less than it appears to. `T_room` was never the under resolved
+state, so it gains only a factor of two; what the extra state controls is `E_sol`, and it
+earns its keep only when the integrand depends on a state. Here `Q` depends on time and
+parameters alone, so the dense output evaluates it exactly and a trapezoid rule over the
+filtered solve reaches 1.1e-5 relative error at a hundred query points and 2.4e-7 at a
+thousand, for no solver cost. Add the state when the integral feeds back into the
+dynamics; post-process when it does not.
+
+Whatever you settle on, verify it once against a reference solve at `reltol = 1e-10`.
+That check, not the step count, is what shows the recipe is sufficient.
 
 ## Implementation Details
 

docs/src/guides/precision.md

@@ -25,6 +25,24 @@ A magnitude-safe time base keeps full intra-day resolution at every precision. T
 carried as an exact integer day count since J2000 plus a fraction of a day, so low
 precision types never ride on the ~2.45e6 Julian Date.
 
+## Refraction models carry their own precision
+
+A refraction model's pressure and temperature have an element type of their own, and the
+apparent angles are computed at it. The result type is the promotion of the two, so a
+narrow observer paired with a default model widens:
+
+```@example precision
+(
+    typeof(solar_position(obs32, dt, PSA(), HUGHES())),          # Float64 pressure widens
+    typeof(solar_position(obs32, dt, PSA(), HUGHES{Float32}())), # matched, stays narrow
+    typeof(solar_position(obs32, dt, PSA(), ARCHER())),          # no parameters to widen
+)
+```
+
+Construct the model at the observer's precision to keep the result narrow. Models without
+parameters never widen anything, and `DefaultRefraction` builds its model at the observer's
+precision, so the default path is unaffected.
+
 ## Supported types
 
 - `Float64` is the default and the reference. Every algorithm agrees with a 256-bit

docs/src/index.md

@@ -91,7 +91,9 @@ accuracy and implementation status.
 
 The computation runs at the precision of the
 [`Observer{T}`](@ref SolarPosition.Positioning.Observer) element type. `Float32`,
-`Float64`, `Float128`, and `BigFloat` are supported. See the
+`Float64`, `Float128`, and `BigFloat` are supported. A refraction model's own parameter
+type promotes with the observer's, so build the model at the same precision to keep a
+narrow result narrow. See the
 [Numeric Precision](@ref numeric-precision) guide for measured accuracy and runtime of
 every algorithm at each precision, including multithreaded benchmarks.