Numeric Precision

SolarPosition.jl computes at the precision of the Observer{T} element type. Construct the observer with the float type you want and every algorithm runs at that precision:

using SolarPosition
using Dates

dt = DateTime(2024, 6, 21, 12, 0, 0)

obs32 = Observer{Float32}(45.0, 10.0)          # explicit precision
obs64 = Observer(45.0, 10.0)                   # Float64, the default
obsbig = Observer{BigFloat}(45.0, 10.0)

solar_position(obs32, dt).elevation            # Float32 result
67.08369f0

The element type is taken from the arguments, so Observer(45.0f0, 10.0f0) also constructs an Observer{Float32}. The explicit Observer{T} form converts its arguments and avoids literal suffixes.

A magnitude-safe time base keeps full intra-day resolution at every precision. Time is 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.

Supported types

  • Float64 is the default and the reference. Every algorithm agrees with a 256-bit BigFloat computation to about 1e-11 degrees.
  • Float32 trades a little accuracy for a modest speedup. The error stays within each algorithm's own claimed accuracy.
  • Float128 from Quadmath.jl gives quad precision of roughly 1e-30 degrees for PSA, SPA, Walraven, and Iqbal at a large runtime cost.
  • BigFloat gives arbitrary precision. Raise it with setprecision(BigFloat, bits). This is the right tool for generating reference values.
  • Float16 is not supported. The algorithm coefficients overflow its range and results are NaN.

Note that ΔT is only known to about a second and is therefore always evaluated in Float64, which bounds the achievable absolute accuracy regardless of T.

Accuracy and runtime

The table below lists the maximum error against a 256-bit BigFloat reference and the single-threaded runtime relative to Float64 for each algorithm. Errors were measured over a grid of 125 combinations of dates from 2015 to 2035 and latitudes from 70°N to 60°S, covering both elevation and azimuth.

AlgorithmFloat32 errorFloat32 runtimeFloat64 errorFloat128 errorFloat128 runtimeBigFloat runtime
PSA0.011°1.3× faster1.3e-11°7.1e-30°93× slower500× slower
NOAA0.0019°1.2× faster3.5e-12°broken[1]n/a330× slower
Walraven0.0040°1.3× faster6.4e-12°8.6e-30°74× slower400× slower
USNO0.0036°1.2× faster9.5e-12°broken[1]n/a300× slower
SPA0.012°1.5× faster1.8e-11°1.1e-29°118× slower330× slower
Iqbal2.2e-5°1.1× faster2.1e-14°3.1e-32°40× slower180× slower

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 MPFR, allocates every intermediate result, and costs another 3× to 5× on top of Float128 at 256 bits.

Automatic differentiation

The element type is only constrained to Real, so dual numbers flow through every algorithm and solar positions are differentiable with ForwardDiff. See the Automatic Differentiation guide for gradients, panel orientation optimization, and a single axis tracker example.

Combining precision with multithreading

Precision and the OhMyThreads extension compose. Pass a scheduler as the last argument and the vectorized API parallelizes over timestamps at any precision:

using OhMyThreads
using Quadmath

obs = Observer{Float128}(45.0, 10.0)
times = collect(DateTime(2024, 6, 21):Minute(1):DateTime(2024, 6, 22))

positions = solar_position(obs, times, SPA(), NoRefraction(), DynamicScheduler())
positions[1]
SolPos(azimuth=9.3811243509767160759440916938436216°, elevation=-21.011598932136595209954880892522944°, zenith=111.01159893213659520995488089252295°)

Measured on 16 Julia threads on an 8-core AMD Ryzen 7 8845HS with DynamicScheduler():

AlgorithmTSerialThreadedSpeedup
PSAFloat320.084 µs0.024 µs3.6×
PSAFloat640.114 µs0.045 µs2.5×
PSAFloat1289.71 µs1.07 µs9.1×
SPAFloat321.35 µs0.21 µs6.5×
SPAFloat641.94 µs0.30 µs6.4×
SPAFloat128226 µs29.7 µs7.6×

Times are per timestamp. The scaling follows arithmetic intensity. PSA at machine precision is so cheap that writing the output array dominates, and memory bandwidth limits the speedup to about 3×. SPA does enough trigonometry per timestamp to keep all cores busy and reaches the ideal scaling for 8 physical cores. Float128 is pure software arithmetic and scales best of all, so threading recovers most of its cost: threaded Float128 SPA is only about 15× slower than serial Float64 SPA instead of 118×. Threading and reduced precision also compose in the other direction, since threaded Float32 SPA is 9.4× faster than serial Float64.

  • 1Quadmath.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.