Support type generic precision

src/Positioning/Positioning.jl

@@ -8,7 +8,7 @@ refraction corrections.
 """
 module Positioning
 
-using Dates: Dates, datetime2julian, DateTime, Date, daysinmonth, dayofyear
+using Dates: Dates, DateTime, Date, daysinmonth, dayofyear
 using Dates: year, month, day
 using TimeZones: ZonedDateTime, UTC
 using StructArrays: StructArrays
@@ -259,6 +259,17 @@ pos_noaa = solar_position(obs, dt, NOAA())
 - `ZonedDateTime` inputs are automatically converted to UTC
 - For local solar time calculations, use appropriate time zones
 
+# Floating-Point Precision
+The result element type follows the `Observer{T}` element type `T`, and the computation runs
+at that precision:
+- **`Float64`** (default): reference accuracy for every algorithm.
+- **`BigFloat`** (and other wide types, e.g. `Float128`): genuine extended precision for every
+  algorithm — use a higher `setprecision` for more correct digits.
+- **`Float32`**: accurate and faster for every algorithm (a magnitude-safe time base keeps full
+  intra-day resolution instead of riding on the ~2.45e6 Julian Date).
+- **`Float16`**: experimental — its range/precision is too small for these algorithms
+  (overflow and inverse-trig domain errors are likely). Use `Float32` or wider.
+
 See also: [`solar_position!`](@ref), [`Observer`](@ref), [`PSA`](@ref), [`NOAA`](@ref)
 """
 function solar_position end
@@ -446,6 +457,7 @@ result_type(::Type{<:SolarAlgorithm}, ::Type{<:RefractionAlgorithm}, ::Type{T})
     ApparentSolPos{T}
 
 include("utils.jl")
+include("timebase.jl")
 include("deltat.jl")
 include("psa.jl")
 include("noaa.jl")

src/Positioning/deltat.jl

@@ -200,3 +200,14 @@ end
 function calculate_deltat(datetime::ZonedDateTime)
     return calculate_deltat(DateTime(datetime, UTC))
 end
+
+# Type-generic entry points: keep the Float64 polynomial value (ΔT is uncertain to ~seconds),
+# convert to T so the result stays type-stable through the pipeline.
+calculate_deltat(::Type{T}, year::Real, month::Real) where {T <: AbstractFloat} =
+    T(calculate_deltat(year, month))
+
+calculate_deltat(::Type{T}, date::Union{DateTime, Date}) where {T <: AbstractFloat} =
+    T(calculate_deltat(date))
+
+calculate_deltat(::Type{T}, datetime::ZonedDateTime) where {T <: AbstractFloat} =
+    T(calculate_deltat(datetime))

src/Positioning/noaa.jl

@@ -36,83 +36,85 @@ 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(dt)
+        calculate_deltat(T, dt)
     else
-        alg.delta_t
+        T(alg.delta_t)
     end
 
-    # convert to Julian date and Julian century
-    jd = datetime2julian(dt)
-    jc = (jd - 2451545.0) / 36525.0
+    # Julian century since J2000.0, at precision T
+    jc = julian_century(T, dt)
 
     # mean longitude of the sun [degrees]
-    mean_long = mod(280.46646 + jc * (36000.76983 + jc * 0.0003032), 360.0)
+    mean_long = mod(T(280.46646) + jc * (T(36000.76983) + jc * T(0.0003032)), 360)
 
     # mean anomaly [degrees]
-    mean_anom = 357.52911 + jc * (35999.05029 - 0.0001537 * jc)
+    mean_anom = T(357.52911) + jc * (T(35999.05029) - T(0.0001537) * jc)
 
     # cccentricity of Earth's orbit
-    eccent = 0.016708634 - jc * (0.000042037 + 0.0000001267 * jc)
+    eccent = T(0.016708634) - jc * (T(0.000042037) + T(0.0000001267) * jc)
 
     # sun equation of center [degrees]
     sun_eq_ctr = (
-        sind(mean_anom) * (1.914602 - jc * (0.004817 + 0.000014 * jc)) +
-            sind(2 * mean_anom) * (0.019993 - 0.000101 * jc) +
-            sind(3 * mean_anom) * 0.000289
+        sind(mean_anom) * (T(1.914602) - jc * (T(0.004817) + T(0.000014) * jc)) +
+            sind(2 * mean_anom) * (T(0.019993) - T(0.000101) * jc) +
+            sind(3 * mean_anom) * T(0.000289)
     )
 
     # sun true/apparent longitude [degrees]
     sun_true_long = mean_long + sun_eq_ctr
-    sun_app_long = sun_true_long - 0.00569 - 0.00478 * sind(125.04 - 1934.136 * jc)
+    sun_app_long = sun_true_long - T(0.00569) - T(0.00478) * sind(T(125.04) - T(1934.136) * jc)
 
     # mean obliquity of ecliptic [degrees]
     mean_obliq =
-        23.0 +
-        (26.0 + (21.448 - jc * (46.815 + jc * (0.00059 - jc * 0.001813))) / 60.0) / 60.0
+        T(23) +
+        (T(26) + (T(21.448) - jc * (T(46.815) + jc * (T(0.00059) - jc * T(0.001813)))) / 60) / 60
 
     # obliquity correction [degrees]
-    obliq_corr = mean_obliq + 0.00256 * cosd(125.04 - 1934.136 * jc)
-    sun_declin = asind(sind(obliq_corr) * sind(sun_app_long))
+    obliq_corr = mean_obliq + T(0.00256) * cosd(T(125.04) - T(1934.136) * jc)
+    sun_declin = asind(unit_clamp(sind(obliq_corr) * sind(sun_app_long)))
 
     # equation of time [minutes]
-    var_y = tand(obliq_corr / 2.0)^2
+    var_y = tand(obliq_corr / 2)^2
     eot =
-        4.0 * rad2deg(
-        var_y * sind(2.0 * mean_long) - 2.0 * eccent * sind(mean_anom) +
-            4.0 * eccent * var_y * sind(mean_anom) * cosd(2.0 * mean_long) -
-            0.5 * var_y^2 * sind(4.0 * mean_long) - 1.25 * eccent^2 * sind(2.0 * mean_anom),
+        4 * rad2deg(
+        var_y * sind(2 * mean_long) - 2 * eccent * sind(mean_anom) +
+            4 * eccent * var_y * sind(mean_anom) * cosd(2 * mean_long) -
+            T(0.5) * var_y^2 * sind(4 * mean_long) - T(1.25) * eccent^2 * sind(2 * mean_anom),
     )
 
     # true solar time [minutes]
-    hour_frac = fractional_hour(dt)
-    minutes = hour_frac * 60.0
-    true_solar_time = mod(minutes + eot + 4.0 * obs.longitude, 1440.0)
+    hour_frac = fractional_hour(T, dt)
+    minutes = hour_frac * 60
+    true_solar_time = mod(minutes + eot + 4 * obs.longitude, 1440)
 
     # hour angle [degrees]
     # true_solar_time is in [0, 1440) minutes, so true_solar_time/4 is in [0, 360) degrees
     # Convert to standard hour angle range (-180, 180] where 0 is solar noon
-    hour_angle = true_solar_time / 4.0 - 180.0
+    hour_angle = true_solar_time / 4 - 180
 
     # zenith angle [degrees]
     zenith = acosd(
-        obs.sin_lat * sind(sun_declin) + obs.cos_lat * cosd(sun_declin) * cosd(hour_angle),
+        unit_clamp(
+            obs.sin_lat * sind(sun_declin) +
+                obs.cos_lat * cosd(sun_declin) * cosd(hour_angle),
+        ),
     )
 
     # azimuth angle [degrees]
     azimuth_numerator = obs.sin_lat * cosd(zenith) - sind(sun_declin)
     azimuth_denominator = obs.cos_lat * sind(zenith)
 
-    azimuth = if hour_angle > 0.0
-        mod(acosd(azimuth_numerator / azimuth_denominator) + 180.0, 360.0)
+    azimuth = if hour_angle > 0
+        mod(acosd(unit_clamp(azimuth_numerator / azimuth_denominator)) + 180, 360)
     else
-        mod(540.0 - acosd(azimuth_numerator / azimuth_denominator), 360.0)
+        mod(540 - acosd(unit_clamp(azimuth_numerator / azimuth_denominator)), 360)
     end
 
-    return SolPos{T}(azimuth, 90.0 - zenith, zenith)
+    return SolPos{T}(azimuth, 90 - zenith, zenith)
 end
 
-function _solar_position(obs, dt, alg::NOAA, ::DefaultRefraction)
-    return _solar_position(obs, dt, alg, HUGHES())
+function _solar_position(obs::Observer{T}, dt, alg::NOAA, ::DefaultRefraction) where {T}
+    return _solar_position(obs, dt, alg, HUGHES{T}())
 end
 
 # NOAA with DefaultRefraction returns ApparentSolPos (uses HUGHES refraction)

src/Positioning/psa.jl

@@ -65,13 +65,12 @@ PSA() = PSA(2020)
 end
 
 function _solar_position(obs::Observer{T}, dt::DateTime, alg::PSA) where {T}
-    # Get parameters as tuple (allocation-free)
+    # Get parameters as tuple (allocation-free), at the observer's precision T
     p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13, p14, p15 =
-        get_psa_params(alg.coeffs)
+        map(T, get_psa_params(alg.coeffs))
 
-    # elapsed julian days (n) since J2000.0
-    jd = datetime2julian(dt)
-    n = jd - 2451545.0                                                  # Eq. 2
+    # elapsed julian days (n) since J2000.0, at precision T
+    n = julian_day_j2000(T, dt)                                         # Eq. 2
 
     # ecliptic coordinates of the sun
     # ecliptic longitude (λₑ), and obliquity of the ecliptic (ϵ)
@@ -87,26 +86,26 @@ function _solar_position(obs::Observer{T}, dt::DateTime, alg::PSA) where {T}
     (sin_λₑ, cos_λₑ) = sincos(λₑ)
     ra = atan(cos_ϵ * sin_λₑ, cos_λₑ)                                   # Eq. 8
     ra = mod2pi(ra)
-    δ = asin(sin_ϵ * sin_λₑ)                                            # Eq. 9
+    δ = asin(unit_clamp(sin_ϵ * sin_λₑ))                                # Eq. 9
 
     # computes the local coordinates: azimuth (γ) and zenith angle (θz)
     λt = rad2deg(obs.longitude_rad)
     cos_lat = obs.cos_lat
     sin_lat = obs.sin_lat
 
-    hour = fractional_hour(dt)
+    hour = fractional_hour(T, dt)
     gmst = p14 + p15 * n + hour                                         # Eq. 10
     lmst = deg2rad(gmst * 15 + λt)                                      # Eq. 11
     ω = lmst - ra                                                       # Eq. 12
     (sin_δ, cos_δ) = sincos(δ)
     (sin_ω, cos_ω) = sincos(ω)
-    θz = acos(cos_lat * cos_ω * cos_δ + sin_δ * sin_lat)                # Eq. 13
+    θz = acos(unit_clamp(cos_lat * cos_ω * cos_δ + sin_δ * sin_lat))    # Eq. 13
     γ = atan(-sin_ω, (tan(δ) * cos_lat - sin_lat * cos_ω))              # Eq. 14
 
     # parallax correction
-    θz = θz + (EMR / AU) * sin(θz)                                      # Eq. 15,16
+    θz = θz + T(EMR / AU) * sin(θz)                                     # Eq. 15,16
 
-    return SolPos{T}(mod(rad2deg(γ), 360.0), rad2deg(π / 2 - θz), rad2deg(θz))
+    return SolPos{T}(mod(rad2deg(γ), 360), rad2deg(T(π) / 2 - θz), rad2deg(θz))
 end
 
 function _solar_position(obs, dt, alg::PSA, ::DefaultRefraction)