diff --git a/.lychee.toml b/.lychee.toml
index fe0b7f0f..cca3f552 100644
--- a/.lychee.toml
+++ b/.lychee.toml
@@ -14,6 +14,8 @@ exclude = [
   "^https://www\\.timeanddate\\.com/.*$",
   # This site does not send its intermediate certificate, so TLS verification fails
   "^https://pvpmc\\.sandia\\.gov/.*$",
+  # This site resets connections from automated clients
+  "^https://aa\\.usno\\.navy\\.mil/.*$",
 ]
 
 exclude_path = ["docs/build"]
diff --git a/CHANGELOG.md b/CHANGELOG.md
index cafef5fd..f6ce5c20 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -22,6 +22,13 @@ Releases before v0.5.0 predate this file. See the
   linkable permalinks, and hover tooltips on identifiers that link to their docstrings,
   via [DocumenterCodeBlocks.jl](https://github.com/fredrikekre/DocumenterCodeBlocks.jl)
 - Compat widened to DocStringExtensions 0.10
+- A new algorithm without a `DefaultRefraction` method now falls back to `NoRefraction`
+  [#136](https://github.com/JuliaAstro/SolarPosition.jl/pull/136)
+
+### Fixed
+
+- `Walraven` azimuth is now always in [0°, 360°). It could previously return exactly 360°
+  [#136](https://github.com/JuliaAstro/SolarPosition.jl/pull/136)
 
 ## v0.5.0 - 2026-07-31
 
diff --git a/docs/src/guides/new-algorithm.md b/docs/src/guides/new-algorithm.md
index 25761c8f..44cd139d 100644
--- a/docs/src/guides/new-algorithm.md
+++ b/docs/src/guides/new-algorithm.md
@@ -120,22 +120,14 @@ end
 
 ## [Step 3: Handle Default Refraction](@id step-3-handle-refraction)
 
-Each algorithm must specify how it handles `DefaultRefraction`. There are two common
+Each algorithm decides how it handles `DefaultRefraction`. There are two common
 patterns:
 
 ### Pattern A: No Refraction by Default (like [`PSA`](@ref SolarPosition.Positioning.PSA))
 
-If your algorithm should NOT apply refraction by default:
-
-```julia
-function _solar_position(obs, dt, alg::SimpleAlgorithm, ::DefaultRefraction)
-    return _solar_position(obs, dt, alg, NoRefraction())
-end
-
-# Return type for DefaultRefraction
-result_type(::Type{SimpleAlgorithm}, ::Type{DefaultRefraction}, ::Type{T}) where {T} =
-    SolPos{T}
-```
+If your algorithm should NOT apply refraction by default, there is nothing to define.
+`DefaultRefraction` falls back to [`NoRefraction`](@ref SolarPosition.Refraction.NoRefraction)
+for every `SolarAlgorithm` and returns a `SolPos`.
 
 ### Pattern B: Apply Refraction by Default (like [`NOAA`](@ref SolarPosition.Positioning.NOAA))
 
@@ -302,8 +294,7 @@ Before submitting your algorithm for review, ensure you've completed the followi
 | Algorithm struct   | Subtypes [`SolarAlgorithm`](@ref SolarPosition.Positioning.SolarAlgorithm)                  |
 | Docstring          | Includes `TYPEDEF`, `TYPEDFIELDS`, accuracy, and references                                 |
 | `_solar_position`  | Function implemented with correct signature                                                 |
-| Default refraction | Handling defined for [`DefaultRefraction`](@ref SolarPosition.Refraction.DefaultRefraction) |
-| `result_type`      | Function defined for [`DefaultRefraction`](@ref SolarPosition.Refraction.DefaultRefraction) |
+| Default refraction | Only with a default model: `_solar_position` and `result_type` for `DefaultRefraction`      |
 | Time base          | Uses `timebase.jl` helpers, never `datetime2julian`                                         |
 | Inverse trig       | `asin`/`acos` arguments guarded with `unit_clamp`                                           |
 | Export             | Algorithm exported from `Positioning`, which re-exports automatically                       |
diff --git a/ext/SolarPositionInterpolationsExt.jl b/ext/SolarPositionInterpolationsExt.jl
index 6dd6d481..e68d1da4 100644
--- a/ext/SolarPositionInterpolationsExt.jl
+++ b/ext/SolarPositionInterpolationsExt.jl
@@ -1,14 +1,12 @@
 """
     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.
+Extension that builds cubic B-spline interpolants to provide the `Interpolated` solar
+position algorithm.
 """
 module SolarPositionInterpolationsExt
 
-using SolarPosition.Positioning: Positioning, SPA, calculate_deltat
+using SolarPosition.Positioning: Positioning, SPA
 using Interpolations: interpolate, scale, BSpline, Cubic, Line, OnGrid
 using Dates: Dates, DateTime, Millisecond
 using Base.Threads: @threads
@@ -28,11 +26,7 @@ function Positioning._build_interpolants(
     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
+        δt = Positioning._resolve_deltat(Float64, algorithm, dti)
         p = Positioning._compute_spa_srt_parameters(Float64, dti, δt)
         αs[i] = p.α
         δs[i] = p.δ
diff --git a/ext/SolarPositionTimeZonesExt.jl b/ext/SolarPositionTimeZonesExt.jl
index 76dfbbc1..bc2698b8 100644
--- a/ext/SolarPositionTimeZonesExt.jl
+++ b/ext/SolarPositionTimeZonesExt.jl
@@ -1,14 +1,7 @@
 """
 TimeZones support for SolarPosition.jl.
 
-`ZonedDateTime` inputs and timezone aware sunrise/sunset results live here rather than in
-the package proper, so that loading SolarPosition does not pull in TimeZones and, through
-it, TZJData, Downloads and LibCURL. This is invisible in practice: a `ZonedDateTime` or a
-`TimeZone` cannot be constructed without TimeZones being loaded, which is exactly what
-triggers this extension.
-
-Every method here converts to UTC and delegates to the `DateTime` method in the package,
-converting back into the given zone where the result is itself a time.
+This module adds support for `ZonedDateTime` inputs and timezone aware sunrise/sunset.
 """
 module SolarPositionTimeZonesExt
 
@@ -25,7 +18,9 @@ using SolarPosition.Utilities: SPA, TransitSunriseSunset, _transit_sunrise_sunse
     _transit_sunrise_sunset_impl, next_solar_noon, next_sunrise, next_sunset,
     previous_solar_noon, previous_sunrise, previous_sunset, transit_sunrise_sunset
 
-# ------------------------------------------------------------------ solar position
+const ZonedVector = Union{
+    AbstractVector{ZonedDateTime}, AbstractVector{Union{DateTime, ZonedDateTime}},
+}
 
 function Positioning.solar_position(
         obs::AbstractObserver{T},
@@ -39,7 +34,7 @@ end
 function Positioning.solar_position!(
         pos::StructArrays.StructVector{T},
         obs::AbstractObserver,
-        dts::AbstractVector{ZonedDateTime},
+        dts::ZonedVector,
         alg::SolarAlgorithm = PSA(),
         refraction::RefractionAlgorithm = DefaultRefraction(),
     ) where {T <: AbstractSolPos}
@@ -49,36 +44,9 @@ function Positioning.solar_position!(
     return pos
 end
 
-function Positioning.solar_position!(
-        pos::StructArrays.StructVector{T},
-        obs::AbstractObserver,
-        dts::AbstractVector{Union{DateTime, ZonedDateTime}},
-        alg::SolarAlgorithm = PSA(),
-        refraction::RefractionAlgorithm = DefaultRefraction(),
-    ) where {T <: AbstractSolPos}
-    @inbounds for i in eachindex(dts, pos)
-        pos[i] = solar_position(obs, dts[i], alg, refraction)
-    end
-    return pos
-end
-
-function Positioning.solar_position(
-        obs::AbstractObserver{T},
-        dts::AbstractVector{ZonedDateTime},
-        alg::SolarAlgorithm = PSA(),
-        refraction::RefractionAlgorithm = DefaultRefraction(),
-    ) where {T <: Real}
-    RetType = result_type(
-        typeof(alg), typeof(refraction), _result_eltype(T, refraction),
-    )
-    pos = StructArrays.StructVector{RetType}(undef, length(dts))
-    solar_position!(pos, obs, dts, alg, refraction)
-    return pos
-end
-
 function Positioning.solar_position(
         obs::AbstractObserver{T},
-        dts::AbstractVector{Union{DateTime, ZonedDateTime}},
+        dts::ZonedVector,
         alg::SolarAlgorithm = PSA(),
         refraction::RefractionAlgorithm = DefaultRefraction(),
     ) where {T <: Real}
@@ -101,7 +69,6 @@ Positioning.calculate_deltat(datetime::ZonedDateTime) =
 Positioning.calculate_deltat(::Type{T}, datetime::ZonedDateTime) where {T <: Real} =
     T(calculate_deltat(datetime))
 
-# the interpolant span accepts zoned endpoints
 Positioning._as_utc(zdt::ZonedDateTime) = DateTime(zdt, UTC)
 
 # --------------------------------------------------------- sunrise, sunset, transit
diff --git a/src/Positioning/Positioning.jl b/src/Positioning/Positioning.jl
index 40d3ec77..9e8b1986 100644
--- a/src/Positioning/Positioning.jl
+++ b/src/Positioning/Positioning.jl
@@ -332,6 +332,10 @@ function _solar_position(obs, dt, alg::SolarAlgorithm, ::NoRefraction)
     return _solar_position(obs, dt, alg)
 end
 
+function _solar_position(obs, dt, alg::SolarAlgorithm, ::DefaultRefraction)
+    return _solar_position(obs, dt, alg, NoRefraction())
+end
+
 function _solar_position(obs, dt, alg::SolarAlgorithm, refraction::RefractionAlgorithm)
     pos = _solar_position(obs, dt, alg)
 
@@ -449,6 +453,7 @@ end
 
 # helper function to determine return type based on refraction
 result_type(::Type{<:SolarAlgorithm}, ::Type{NoRefraction}, ::Type{T}) where {T} = SolPos{T}
+result_type(::Type{<:SolarAlgorithm}, ::Type{DefaultRefraction}, ::Type{T}) where {T} = SolPos{T}
 result_type(::Type{<:SolarAlgorithm}, ::Type{<:RefractionAlgorithm}, ::Type{T}) where {T} =
     ApparentSolPos{T}
 
diff --git a/src/Positioning/deltat.jl b/src/Positioning/deltat.jl
index 4416248e..86dd3357 100644
--- a/src/Positioning/deltat.jl
+++ b/src/Positioning/deltat.jl
@@ -204,3 +204,9 @@ calculate_deltat(::Type{T}, year::Real, month::Real) where {T <: Real} =
 
 calculate_deltat(::Type{T}, date::Union{DateTime, Date}) where {T <: Real} =
     T(calculate_deltat(date))
+
+# ΔT [seconds] for an algorithm with a `delta_t` field: its fixed value when set, otherwise
+# the ΔT model evaluated at `dt`.
+@inline function _resolve_deltat(::Type{T}, alg, dt) where {T <: Real}
+    return alg.delta_t === nothing ? calculate_deltat(T, dt) : T(alg.delta_t)
+end
diff --git a/src/Positioning/interpolated.jl b/src/Positioning/interpolated.jl
index 666c29d2..2398f9cd 100644
--- a/src/Positioning/interpolated.jl
+++ b/src/Positioning/interpolated.jl
@@ -149,47 +149,33 @@ function _solar_position(
 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)
+    return _solar_position(_spa_observer(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),
-        ),
-    )
+    return _solar_position(obs, dt, alg, _spa_refraction(T, alg.algorithm))
 end
 
 function solar_position!(
         pos::StructArrays.StructVector{S},
         obs::AbstractObserver{T},
         dts::AbstractVector{DateTime},
-        alg::Interpolated,
+        alg::Union{SPA, Interpolated},
         refraction::RefractionAlgorithm = DefaultRefraction(),
     ) where {S <: AbstractSolPos, T <: Real}
-    spa_obs = SPAObserver{T}(obs.latitude, obs.longitude, obs.altitude)
+    spa_obs = _spa_observer(obs)
     @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} =
+result_type(::Type{<:Interpolated}, ::Type{DefaultRefraction}, ::Type{T}) where {T} =
     ApparentSolPos{T}
 
 """
@@ -208,7 +194,7 @@ function solar_rate(
         dt::DateTime,
         alg::Interpolated,
     ) where {T <: Real}
-    spa_obs = SPAObserver{T}(obs.latitude, obs.longitude, obs.altitude)
+    spa_obs = _spa_observer(obs)
     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
diff --git a/src/Positioning/iqbal.jl b/src/Positioning/iqbal.jl
index b1cffe9b..cfbde849 100644
--- a/src/Positioning/iqbal.jl
+++ b/src/Positioning/iqbal.jl
@@ -61,10 +61,3 @@ function _solar_position(obs::Observer{T}, dt::DateTime, ::Iqbal) where {T <: Re
 
     return SolPos{T}(azimuth, 90 - zenith, zenith)
 end
-
-function _solar_position(obs::AbstractObserver, dt, alg::Iqbal, ::DefaultRefraction)
-    return _solar_position(obs, dt, alg, NoRefraction())
-end
-
-# Iqbal applies no refraction by default, so DefaultRefraction returns SolPos
-result_type(::Type{Iqbal}, ::Type{DefaultRefraction}, ::Type{T}) where {T} = SolPos{T}
diff --git a/src/Positioning/michalsky.jl b/src/Positioning/michalsky.jl
index 5e709e54..0ca255cf 100644
--- a/src/Positioning/michalsky.jl
+++ b/src/Positioning/michalsky.jl
@@ -4,8 +4,7 @@
 Michalsky solar position algorithm.
 
 Implements the approximate solar position algorithm of the Astronomical Almanac. The
-default refraction model is [`MICHALSKY`](@ref), matching the original publication, so
-`solar_position` returns an [`ApparentSolPos`](@ref) unless another model is given.
+default refraction model is [`MICHALSKY`](@ref).
 
 # Accuracy
 ±0.01° between 1950 and 2050 with the original Julian date formulation.
@@ -48,30 +47,35 @@ function Michalsky(; spencer_correction::Bool = true, julian_date::Symbol = :ori
     return Michalsky(spencer_correction, julian_date)
 end
 
-# Days since J2000.0 noon, the quantity the almanac series are expressed in.
+# Days since J2000.0 noon. `leapdiv` counts the leap days since 1949:
+# - Michalsky floors (`fld`)
+# - SOLPOS truncates towards zero (`div`).
+# The definitions also differ before 1949.
 #
 # The paper builds a Julian Date from a naive leap year count since 1949, which would
 # materialise the ~2.43e6 Julian Date magnitude at precision T. The constant J2000 offset is
 # folded into the integer day count instead, so nothing larger than the day count itself is
 # ever represented at T. The half day is the paper's midnight epoch against J2000 noon.
-@inline function _michalsky_days(::Type{T}, dt::DateTime, julian_date::Symbol) where {T <: Real}
-    julian_date === :original || return julian_day_j2000(T, dt)
+@inline function _almanac_days(::Type{T}, dt::DateTime, leapdiv::F) where {T <: Real, F}
     Δyear = year(dt) - 1949
-    days = Δyear * 365 + fld(Δyear, 4) + dayofyear(dt) - 18629
+    days = Δyear * 365 + leapdiv(Δyear, 4) + dayofyear(dt) - 18629
     return T(days) + T(0.5) + fractional_hour(T, dt) / 24
 end
 
-function _solar_position(obs::Observer{T}, dt::DateTime, alg::Michalsky) where {T <: Real}
-    n = _michalsky_days(T, dt, alg.julian_date)
+@inline function _michalsky_days(::Type{T}, dt::DateTime, julian_date::Symbol) where {T <: Real}
+    julian_date === :original || return julian_day_j2000(T, dt)
+    return _almanac_days(T, dt, fld)
+end
+
+@inline function _almanac_ecliptic_longitude(q::T, g::T) where {T <: Real}
+    return mod(q + T(1.915) * sind(g) + T(0.02) * sind(2 * g), 360)
+end
 
-    # mean longitude, mean anomaly and ecliptic longitude in degrees
+
+@inline function _almanac_equatorial(n::T, dt::DateTime) where {T <: Real}
     mean_longitude = mod(T(280.46) + T(0.9856474) * n, 360)
     mean_anomaly = mod(T(357.528) + T(0.9856003) * n, 360)
-    ecliptic_longitude = mod(
-        mean_longitude + T(1.915) * sind(mean_anomaly) +
-            T(0.02) * sind(2 * mean_anomaly),
-        360,
-    )
+    ecliptic_longitude = _almanac_ecliptic_longitude(mean_longitude, mean_anomaly)
 
     # obliquity of the ecliptic in degrees
     obliquity = T(23.439) - T(0.0000004) * n
@@ -83,8 +87,17 @@ function _solar_position(obs::Observer{T}, dt::DateTime, alg::Michalsky) where {
     declination = asind(unit_clamp(sin_ε * sin_λ))
     right_ascension = mod(atand(cos_ε * sin_λ, cos_λ), 360)
 
-    # Greenwich and local mean sidereal time in hours
+    # Greenwich mean sidereal time in hours
     gmst = mod(T(6.697375) + T(0.0657098242) * n + fractional_hour(T, dt), 24)
+
+    return (declination, right_ascension, gmst)
+end
+
+function _solar_position(obs::Observer{T}, dt::DateTime, alg::Michalsky) where {T <: Real}
+    n = _michalsky_days(T, dt, alg.julian_date)
+    (declination, right_ascension, gmst) = _almanac_equatorial(n, dt)
+
+    # local mean sidereal time in hours
     lmst = mod(gmst + obs.longitude / 15, 24)
 
     # hour angle in hours, wrapped to [-12, 12)
@@ -97,7 +110,7 @@ function _solar_position(obs::Observer{T}, dt::DateTime, alg::Michalsky) where {
     azimuth = asind(unit_clamp(-cos_δ * sin_ω / cosd(elevation)))
 
     azimuth = if alg.spencer_correction
-        _spencer_quadrant(azimuth, elevation, sin_δ, obs.sin_lat)
+        _spencer_quadrant(azimuth, sin_δ - sind(elevation) * obs.sin_lat)
     else
         _original_quadrant(azimuth, elevation, sin_δ, obs.sin_lat, hour_angle)
     end
@@ -105,10 +118,8 @@ function _solar_position(obs::Observer{T}, dt::DateTime, alg::Michalsky) where {
     return SolPos{T}(azimuth, elevation, 90 - elevation)
 end
 
-# Spencer (1989) quadrant correction. The sign of the azimuth cosine picks the branch, which
-# makes the result valid in both hemispheres.
-@inline function _spencer_quadrant(azimuth::T, elevation, sin_δ, sin_lat) where {T <: Real}
-    cos_azimuth = sin_δ - sind(elevation) * sin_lat
+# Spencer (1989) quadrant correction for an arcsine azimuth in [-90, 90] degrees
+@inline function _spencer_quadrant(azimuth::T, cos_azimuth) where {T <: Real}
     if cos_azimuth < 0
         azimuth = 180 - azimuth
     elseif sind(azimuth) < 0
@@ -139,7 +150,5 @@ function _solar_position(obs::AbstractObserver, dt, alg::Michalsky, ::DefaultRef
     return _solar_position(obs, dt, alg, MICHALSKY())
 end
 
-# Michalsky applies the MICHALSKY refraction model by default, matching the paper, so
-# DefaultRefraction returns ApparentSolPos
 result_type(::Type{Michalsky}, ::Type{DefaultRefraction}, ::Type{T}) where {T} =
     ApparentSolPos{T}
diff --git a/src/Positioning/noaa.jl b/src/Positioning/noaa.jl
index b5abad5e..e86d21ba 100644
--- a/src/Positioning/noaa.jl
+++ b/src/Positioning/noaa.jl
@@ -35,12 +35,6 @@ NOAA() = NOAA(67.0)  # default delta_t value (2020 default from pvlib)
 
 
 function _solar_position(obs::Observer{T}, dt::DateTime, alg::NOAA) where {T}
-    δt = if alg.delta_t === nothing
-        calculate_deltat(T, dt)
-    else
-        T(alg.delta_t)
-    end
-
     # Julian century since J2000.0, at precision T
     jc = julian_century(T, dt)
 
@@ -117,6 +111,5 @@ function _solar_position(obs::Observer{T}, dt, alg::NOAA, ::DefaultRefraction) w
     return _solar_position(obs, dt, alg, HUGHES{T}())
 end
 
-# NOAA with DefaultRefraction returns ApparentSolPos (uses HUGHES refraction)
 result_type(::Type{NOAA}, ::Type{DefaultRefraction}, ::Type{T}) where {T} =
     ApparentSolPos{T}
diff --git a/src/Positioning/psa.jl b/src/Positioning/psa.jl
index ce4d9400..d9851e9f 100644
--- a/src/Positioning/psa.jl
+++ b/src/Positioning/psa.jl
@@ -108,10 +108,3 @@ function _solar_position(obs::Observer{T}, dt::DateTime, alg::PSA) where {T}
 
     return SolPos{T}(mod(rad2deg(γ), 360), rad2deg(T(π) / 2 - θz), rad2deg(θz))
 end
-
-function _solar_position(obs, dt, alg::PSA, ::DefaultRefraction)
-    return _solar_position(obs, dt, alg, NoRefraction())
-end
-
-# PSA with DefaultRefraction returns SolPos (no refraction by default)
-result_type(::Type{PSA}, ::Type{DefaultRefraction}, ::Type{T}) where {T} = SolPos{T}
diff --git a/src/Positioning/solpos.jl b/src/Positioning/solpos.jl
index 174914d0..3c6db9f5 100644
--- a/src/Positioning/solpos.jl
+++ b/src/Positioning/solpos.jl
@@ -47,39 +47,12 @@ solar_position(obs, DateTime(1999, 7, 22, 14, 45, 37), SOLPOS(), refraction)
 """
 struct SOLPOS <: SolarAlgorithm end
 
-# Days since J2000.0 noon from the original Julian date, whose leap year count truncates
-# towards zero. The J2000 epoch is subtracted in integer arithmetic, so only the small day
-# count is converted to T and the large Julian date never is.
-@inline function _solpos_days(::Type{T}, dt::DateTime) where {T <: Real}
-    Δyear = year(dt) - 1949
-    days = Δyear * 365 + div(Δyear, 4) + dayofyear(dt) - 18629
-    return T(days) + T(0.5) + fractional_hour(T, dt) / 24
-end
-
 function _solar_position(obs::Observer{T}, dt::DateTime, ::SOLPOS) where {T <: Real}
-    n = _solpos_days(T, dt)
-
-    # mean longitude, mean anomaly and ecliptic longitude in degrees
-    mean_longitude = mod(T(280.46) + T(0.9856474) * n, 360)
-    mean_anomaly = mod(T(357.528) + T(0.9856003) * n, 360)
-    ecliptic_longitude = mod(
-        mean_longitude + T(1.915) * sind(mean_anomaly) +
-            T(0.02) * sind(2 * mean_anomaly),
-        360,
-    )
-
-    # obliquity of the ecliptic in degrees
-    obliquity = T(23.439) - T(4.0e-7) * n
-
-    (sin_λ, cos_λ) = sincosd(ecliptic_longitude)
-    (sin_ε, cos_ε) = sincosd(obliquity)
 
-    # declination and right ascension in degrees
-    declination = asind(unit_clamp(sin_ε * sin_λ))
-    right_ascension = mod(atand(cos_ε * sin_λ, cos_λ), 360)
+    n = _almanac_days(T, dt, div)
+    (declination, right_ascension, gmst) = _almanac_equatorial(n, dt)
 
-    # Greenwich mean sidereal time in hours, local mean sidereal time in degrees
-    gmst = mod(T(6.697375) + T(0.0657098242) * n + fractional_hour(T, dt), 24)
+    # local mean sidereal time in degrees
     lmst = mod(15 * gmst + obs.longitude, 360)
 
     # hour angle in degrees, wrapped once into [-180, 180] like the original code
@@ -116,6 +89,5 @@ function _solar_position(obs::Observer{T}, dt, alg::SOLPOS, ::DefaultRefraction)
     return _solar_position(obs, dt, alg, HUGHES{T}(T(101325), T(15)))
 end
 
-# SOLPOS applies its Zimmerman refraction correction by default (which is Hughes)
 result_type(::Type{SOLPOS}, ::Type{DefaultRefraction}, ::Type{T}) where {T} =
     ApparentSolPos{T}
diff --git a/src/Positioning/spa.jl b/src/Positioning/spa.jl
index d765b857..d7840401 100644
--- a/src/Positioning/spa.jl
+++ b/src/Positioning/spa.jl
@@ -403,11 +403,7 @@ function _solar_position(
         dt::DateTime,
         alg::SPA,
     ) where {T <: Real}
-    δt::T = if alg.delta_t === nothing
-        calculate_deltat(T, dt)
-    else
-        T(alg.delta_t)
-    end
+    δt = _resolve_deltat(T, alg, dt)
 
     # Compute sidereal time, right ascension, declination, and related parameters
     (; ν, α, δ, R) = _compute_spa_srt_parameters(T, dt, δt)
@@ -415,9 +411,22 @@ function _solar_position(
     return _spa_topocentric(obs, ν, α, δ, R)
 end
 
+# SPA and Interpolated run on an SPAObserver, which caches the parallax terms
+@inline function _spa_observer(obs::AbstractObserver{T}) where {T <: Real}
+    return SPAObserver{T}(obs.latitude, obs.longitude, obs.altitude)
+end
+
+# SPARefraction built at precision T from the atmosphere fields of an SPA
+@inline function _spa_refraction(::Type{T}, alg::SPA) where {T <: Real}
+    return SPARefraction{T}(
+        pressure = T(alg.pressure),
+        temperature = T(alg.temperature),
+        atmos_refract = T(alg.atmos_refract),
+    )
+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)
+    return _solar_position(_spa_observer(obs), dt, alg)
 end
 
 function _solar_position(
@@ -426,33 +435,7 @@ function _solar_position(
         alg::SPA,
         ::DefaultRefraction,
     ) where {T <: Real}
-    return _solar_position(
-        obs,
-        dt,
-        alg,
-        SPARefraction{T}(
-            pressure = T(alg.pressure),
-            temperature = T(alg.temperature),
-            atmos_refract = T(alg.atmos_refract),
-        ),
-    )
-end
-
-function solar_position!(
-        pos::StructArrays.StructVector{S},
-        obs::AbstractObserver{T},
-        dts::AbstractVector{DateTime},
-        alg::SPA,
-        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
+    return _solar_position(obs, dt, alg, _spa_refraction(T, alg))
 end
 
-# SPA returns SolPos with NoRefraction, ApparentSolPos with any refraction
-result_type(::Type{SPA}, ::Type{NoRefraction}, ::Type{T}) where {T} = SolPos{T}
-result_type(::Type{SPA}, ::Type{<:RefractionAlgorithm}, ::Type{T}) where {T} =
-    ApparentSolPos{T}
+result_type(::Type{SPA}, ::Type{DefaultRefraction}, ::Type{T}) where {T} = ApparentSolPos{T}
diff --git a/src/Positioning/usno.jl b/src/Positioning/usno.jl
index e4123412..f0848b39 100644
--- a/src/Positioning/usno.jl
+++ b/src/Positioning/usno.jl
@@ -32,11 +32,7 @@ USNO() = USNO(67.0, 1)  # default delta_t value and gmst_option
 
 
 function _solar_position(obs::Observer{T}, dt::DateTime, alg::USNO) where {T <: Real}
-    δt::T = if alg.delta_t === nothing
-        calculate_deltat(T, dt)
-    else
-        T(alg.delta_t)
-    end
+    δt = _resolve_deltat(T, alg, dt)
 
     # days since J2000.0 (UT), magnitude-safe at precision T
     D = julian_day_j2000(T, dt)
@@ -50,8 +46,7 @@ function _solar_position(obs::Observer{T}, dt::DateTime, alg::USNO) where {T <:
     q = mod(q, 360)
 
     # geocentric apparent ecliptic longitude of the sun (adjusted for aberration) [deg]
-    L = q + T(1.915) * sind(g) + T(0.02) * sind(2 * g)
-    L = mod(L, 360)
+    L = _almanac_ecliptic_longitude(q, g)
 
     # mean obliquity of the ecliptic [deg]
     ϵ = T(23.439) - T(0.00000036) * D
@@ -112,10 +107,3 @@ function _solar_position(obs::Observer{T}, dt::DateTime, alg::USNO) where {T <:
 
     return SolPos{T}(mod(azimuth, 360), elevation, 90 - elevation)
 end
-
-function _solar_position(obs, dt, alg::USNO, ::DefaultRefraction)
-    return _solar_position(obs, dt, alg, NoRefraction())
-end
-
-# USNO with DefaultRefraction returns SolPos (no refraction by default)
-result_type(::Type{USNO}, ::Type{DefaultRefraction}, ::Type{T}) where {T} = SolPos{T}
diff --git a/src/Positioning/walraven.jl b/src/Positioning/walraven.jl
index 40f646d0..2fb5d11b 100644
--- a/src/Positioning/walraven.jl
+++ b/src/Positioning/walraven.jl
@@ -87,23 +87,8 @@ function _solar_position(obs::Observer{T}, dt::DateTime, ::Walraven) where {T}
     az = rad2deg(asin(unit_clamp(cos_d * sin_ha / cos_el)))
 
     # azimuth quadrant assignment - Spencer (1989) correction for all longitudes
-    cos_az = sin_d - sin_el * obs.sin_lat
-
-    if (cos_az >= 0) && (sin(deg2rad(az)) < 0)
-        az = 360 + az
-    end
-
-    if cos_az < 0
-        az = 180 - az
-    end
+    az = _spencer_quadrant(az, sin_d - sin_el * obs.sin_lat)
 
     elevation_deg = rad2deg(el)
     return SolPos{T}(az, elevation_deg, 90 - elevation_deg)
 end
-
-function _solar_position(obs, dt, alg::Walraven, ::DefaultRefraction)
-    return _solar_position(obs, dt, alg, NoRefraction())
-end
-
-# Walraven with DefaultRefraction returns SolPos (no refraction by default)
-result_type(::Type{Walraven}, ::Type{DefaultRefraction}, ::Type{T}) where {T} = SolPos{T}
diff --git a/src/Utilities/Utilities.jl b/src/Utilities/Utilities.jl
index a384c349..03c2c697 100644
--- a/src/Utilities/Utilities.jl
+++ b/src/Utilities/Utilities.jl
@@ -1,7 +1,7 @@
 module Utilities
 
 using DocStringExtensions: TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES
-using ..Positioning: Observer, SPA, SolarAlgorithm, calculate_deltat
+using ..Positioning: Observer, SPA, SolarAlgorithm
 import Dates: DateTime, Date, Day
 
 include("spa.jl")
diff --git a/src/Utilities/spa.jl b/src/Utilities/spa.jl
index 08c994c8..3f914ce8 100644
--- a/src/Utilities/spa.jl
+++ b/src/Utilities/spa.jl
@@ -2,7 +2,7 @@
 
 using Dates: Dates
 
-using ..Positioning: _compute_spa_srt_parameters
+using ..Positioning: _compute_spa_srt_parameters, _resolve_deltat
 
 _frac_to_dt(dt_midnight, frac) =
     dt_midnight + Dates.Second(round(Int, frac * 86_400))
@@ -55,11 +55,7 @@ function _transit_sunrise_sunset_impl(
     # This creates a new DateTime at midnight
     dt_midnight = DateTime(Date(dt))
 
-    δt = if alg.delta_t === nothing
-        calculate_deltat(T, dt_midnight)
-    else
-        T(alg.delta_t)
-    end
+    δt = _resolve_deltat(T, alg, dt_midnight)
 
     lon = obs.longitude