langestefan · Oct 25, 2026 · Oct 24, 2026 · Oct 24, 2026 · Oct 24, 2026 · Oct 24, 2026
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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,18 +110,16 @@ 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
 
     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}