From 91c9a01e2653c1d259b214cefedf657ad362a2f3 Mon Sep 17 00:00:00 2001 From: Stefan de Lange Date: Sun, 13 Sep 2026 02:00:02 +0200 Subject: [PATCH 01/13] Ignore juliac build products and Python envs under lib/ Co-Authored-By: Claude Opus 5 (1M context) --- .gitignore | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/.gitignore b/.gitignore index 11baeb20..dc661877 100644 --- a/.gitignore +++ b/.gitignore @@ -18,3 +18,12 @@ play/Project.toml **/.CondaPkg/ *.mem **/Manifest-*.toml + +# juliac build products +lib/out/ +lib/c/*.o + +# juliac copies the whole project tree into a temp dir to build, and chokes on the +# executable bits inside a virtualenv, so keep Python environments out of lib/ +lib/.venv/ +lib/**/__pycache__/ From 5010609685a6ea475e7a5d02fbcfdcbd2428bf38 Mon Sep 17 00:00:00 2001 From: Stefan de Lange Date: Sun, 13 Sep 2026 02:00:06 +0200 Subject: [PATCH 02/13] Refresh CLAUDE.md architecture notes Co-Authored-By: Claude Opus 5 (1M context) --- CLAUDE.md | 51 ++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 36 insertions(+), 15 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 7b519f7c..1b4ff08f 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -45,7 +45,7 @@ julia --project=docs docs/make.jl **Run benchmarks** (AirspeedVelocity.jl: `benchpkg SolarPosition --rev=main,dirty`, or include `benchmark/benchmarks.jl` and `run(SUITE)` with BenchmarkTools). -Note: the pre-commit `no-commit-to-branch` hook blocks commits to `main` — always work on a branch. `test/`, `docs/`, `benchmark/`, and `examples/` are workspace subprojects (root `Project.toml` `[workspace]`) with their own `Project.toml` that depends on the package via `[sources]` path — no `Pkg.develop` needed. +Note: the pre-commit `no-commit-to-branch` hook blocks commits to `main` — always work on a branch. `test/`, `docs/`, `benchmark/`, and `examples/` are workspace subprojects (root `Project.toml` `[workspace]`) with their own `Project.toml` that depends on the package via `[sources]` path — no `Pkg.develop` needed. CI (`enforce-changelog.yml`) requires every PR to touch `CHANGELOG.md` unless it carries the `skip-changelog` label. ## Architecture @@ -53,43 +53,64 @@ The package provides a unified interface to multiple solar position algorithms. ```text src/ - SolarPosition.jl # Main module, re-exports all submodules + SolarPosition.jl # Main module, re-exports all submodules; holds the docstrings + # and `function ... end` stubs that extensions fill in + # (analemmas!, SolarPositionBlock) Positioning/ # Solar position algorithms Positioning.jl # Observer struct, SolPos/ApparentSolPos types, solar_position() API psa.jl # PSA algorithm (default, ±0.0083°) - noaa.jl, spa.jl, # Other algorithms (NOAA, SPA, Walraven, USNO) - walraven.jl, usno.jl + noaa.jl, spa.jl, # Other algorithms (NOAA, SPA + spa_coefficients.jl tables, + walraven.jl, # Walraven, USNO, Iqbal, Michalsky) + usno.jl, iqbal.jl, + michalsky.jl + interpolated.jl # Interpolated{SPA} wrapper + solar_rate (needs Interpolations ext + # to construct; the struct and query path live here) deltat.jl # Delta T / leap seconds timebase.jl # Magnitude-safe day/century counts since J2000, at precision T utils.jl # unit_clamp for inverse trig, EMR/AU constants Refraction/ # Atmospheric refraction correction models Refraction.jl # Abstract base + interface - hughes.jl, bennett.jl, sg2.jl, spa.jl, ... + hughes.jl, bennett.jl, sg2.jl, spa.jl, archer.jl, michalsky.jl Utilities/ # Sunrise/sunset/transit calculations srt.jl, spa.jl ext/ # Weak dependency extensions (auto-load on `using` of the trigger pkg) + SolarPositionTimeZonesExt.jl # TimeZones → ZonedDateTime input, zoned sunrise/sunset + SolarPositionInterpolationsExt.jl # Interpolations → Interpolated(SPA(); tspan, step) ctor SolarPositionMakieExt.jl # Makie → analemmas!() sun-path plotting (PolarAxis/Axis) SolarPositionOhMyThreadsExt.jl # OhMyThreads → solar_position[!] with an extra ::Scheduler arg SolarPositionModelingToolkitExt.jl # ModelingToolkit/Symbolics → SolarPositionBlock() (t in SECONDS) +dyad/solarpositionblock.dyad # Dyad component declaration mirroring SolarPositionBlock +lib/ # Python bindings: juliac --trim + JuliaLibWrapping (see lib/README.md) + src/libsolarposition.jl # the @api surface; enums for algorithm/refraction selection + build.jl # default and --threads=N builds + python/_extras.py # hand-maintained façade additions + test/ # pytest suite + Julia fixture generator ``` -Minimum Julia: **1.10** (LTS). Extensions only activate once their trigger package is loaded; don't `import` them directly. +Minimum Julia: **1.10** (LTS); CI runs `lts` and `1` on Linux/macOS/Windows. Extensions only activate once their trigger package is loaded; don't `import` them directly. `TimeZones` is a weak dependency: `ZonedDateTime` methods only exist once TimeZones is loaded, which is unavoidable anyway since constructing a `ZonedDateTime` requires it. -**Core API pattern:** `solar_position(obs, dt, alg=PSA(), refraction=DefaultRefraction())` — observer first, then datetime, then algorithm/refraction (both default-able). `dt` may be a single `DateTime`/`ZonedDateTime` or an `AbstractVector` of them (returns a `StructArray`). A table interface (`solar_position(table, obs; dt_col=:datetime)`, mutates the table) and an in-place `solar_position!(pos, obs, dts, alg, refraction)` also exist. The `Utilities` module exports `transit_sunrise_sunset` and `next_`/`previous_` variants of `sunrise`/`sunset`/`solar_noon`. +**Core API pattern:** `solar_position(obs, dt, alg=PSA(), refraction=DefaultRefraction())` — observer first, then datetime, then algorithm/refraction (both default-able). `dt` may be a single `DateTime`/`ZonedDateTime` or an `AbstractVector` of them (returns a `StructArray`). Also: -- **Algorithms and refraction models are singleton dispatch types** — e.g. `PSA()`, `NOAA()`, `SPA()`, `HUGHES(pressure, temperature)`, `NoRefraction()`. New algorithms are added by defining a struct `<: SolarAlgorithm` (or `<: RefractionAlgorithm`) and a `solar_position` method on it. -- **Return type depends on refraction:** `result_type(...)` yields `SolPos{T}` for `NoRefraction`, else `ApparentSolPos{T}` (adds `apparent_elevation`/`apparent_zenith`). -- **Angle convention:** all degrees. Azimuth 0°=North, +clockwise, normalized to [0°, 360°) by every algorithm (some docstrings in `Positioning.jl` still say [-180°, 180°] — the code is authoritative); elevation [-90°, 90°]; zenith = 90° − elevation. Refraction parameters (pressure, temperature) live on the refraction algorithm structs (e.g. `HUGHES(pressure, temperature)`), not on `Observer`. -- `Observer(latitude, longitude; altitude=0.0, horizon=0.0)` precomputes lat/lon trig (`sin_lat`, `cos_lat`) for performance. `horizon` (degrees, or a `deg=>arcmin` pair like `0=>34`) is the horizon depression used for sunrise/sunset. `Observer{T}(lat, lon, ...)` converts any `Real` arguments to precision `T`. -- **Type-generic precision:** every algorithm computes at the `Observer{T}` element type. `Float32`, `Float64`, `Float128`, and `BigFloat` are supported; `Float16` overflows to `NaN`. Inside algorithms always take time via `julian_day_j2000(T, dt)`, `julian_century(T, dt)`, or `fractional_hour(T, dt)` from `timebase.jl`, never `datetime2julian`, so the ~2.45e6 Julian Date magnitude is not materialized at `T`. Guard `asin`/`acos` arguments with `unit_clamp`. ΔT is always evaluated in `Float64`. See `docs/src/guides/precision.md` for measured accuracy and speed. +- `solar_position!(pos, obs, dts, alg, refraction)` fills a preallocated `StructVector`. +- Table interface: `solar_position!(table, obs, alg, refraction; dt_col=:datetime)` **mutates** the table by adding `azimuth`/`elevation`/`zenith` (plus the apparent columns when refraction applies); `solar_position(table, obs, ...)` copies first. Note the observer is the *second* positional argument and `dt_col` is a keyword. +- The `Utilities` module exports `transit_sunrise_sunset`, `transit_sunrise_sunset_seconds` (seconds since UTC midnight at the observer's element type — use this when a `DateTime`'s whole-second rounding would discard an uncertainty or derivative), and `next_`/`previous_` variants of `sunrise`/`sunset`/`solar_noon`. -**Test discovery:** Test files matching `test-*.jl` under `test/` are automatically discovered and wrapped in `@testset`. Reference values live in `expected-values.jl` files alongside algorithm tests. +Key invariants: + +- **Algorithms and refraction models are dispatch types, most carrying options** — `Walraven()`, `Iqbal()`, `ARCHER()`, `MICHALSKY()`, `NoRefraction()` and `DefaultRefraction()` are singletons; the rest have fields: `PSA(coeffs)`, `NOAA(delta_t)`, `USNO(delta_t, gmst_option)`, `SPA(delta_t, pressure, temperature, atmos_refract)`, `Michalsky(spencer_correction, julian_date::Symbol)`, and the parametric `HUGHES/BENNETT/SG2{T}(pressure, temperature)` and `SPARefraction{T}(pressure, temperature, atmos_refract)`. New algorithms are added by defining a struct `<: SolarAlgorithm` (or `<: RefractionAlgorithm`) plus a `_solar_position` method and a `result_type` method on it. See `docs/src/guides/new-algorithm.md`. +- **Return type is resolved at compile time by `result_type(AlgType, RefrType, T)`:** `SolPos{T}` for `NoRefraction`, `ApparentSolPos{T}` (adds `apparent_elevation`/`apparent_zenith`) for an explicit refraction model. `DefaultRefraction()` is resolved *per algorithm*: PSA, Walraven, Iqbal and USNO apply none and return `SolPos`; NOAA (HUGHES), Michalsky (MICHALSKY) and SPA (SPARefraction) return `ApparentSolPos`. Adding an algorithm means adding its `DefaultRefraction` `result_type` too, or the batch path allocates the wrong eltype. +- **Angle convention:** all degrees. Azimuth 0°=North, +clockwise, normalized to [0°, 360°) by every algorithm (the `SolPos`/`ApparentSolPos` field docstrings in `Positioning.jl` still say [-180°, 180°] — the code is authoritative); elevation [-90°, 90°]; zenith = 90° − elevation. Refraction parameters (pressure, temperature) live on the refraction algorithm structs (e.g. `HUGHES(pressure, temperature)`), not on `Observer`, and promote against the observer's element type. +- `Observer(lat, lon; altitude=0, horizon=0)` — also positional `Observer(lat, lon, alt[, horizon])` — precomputes lat/lon trig (`sin_lat`, `cos_lat`) for performance. The element type promotes over the arguments, so integer and mixed inputs work; `Observer{T}(lat, lon, ...)` converts any `Real` arguments to precision `T`. `horizon` (degrees, or a `deg=>arcmin` pair like `0=>34`) is *documented* as the horizon depression for sunrise/sunset but is currently **inert** — it is stored, converted and shown, and never read by any computation in `src/`; `src/Utilities/spa.jl:89` hardcodes `h0 = T(-0.8333)` instead. +- **Type-generic precision:** `Observer{T<:Real}` — deliberately `Real`, not `AbstractFloat`, so `ForwardDiff.Dual` and `Measurements.Measurement` flow through unchanged (no extension needed; see `test/positioning/test-autodiff.jl`, `test-measurements.jl`). Every algorithm computes at the observer's element type: `Float32`, `Float64`, `Float128` (Quadmath, a docs/benchmark-only dep) and `BigFloat` are supported; `Float16` overflows to `NaN`. Inside algorithms always take time via `julian_day_j2000(T, dt)`, `julian_century(T, dt)`, or `fractional_hour(T, dt)` from `timebase.jl`, never `datetime2julian`, so the ~2.45e6 Julian Date magnitude is not materialized at `T`. Guard `asin`/`acos` arguments with `unit_clamp`. ΔT is always evaluated in `Float64`. See `docs/src/guides/precision.md` for measured accuracy and speed. + +**Test discovery:** test files matching `test-*.jl` under `test/` are discovered recursively by `runtests.jl` and wrapped in a `@testset` named after the filename — don't edit `runtests.jl` to add tests. `test/positioning/expected-values.jl` is loaded once by `runtests.jl`; `test/utilities/expected-values.jl` is included by `test-srt.jl` itself. `test/linting.jl` runs Aqua.jl always and JET.jl only on Julia 1.12. **Reference values:** the tables in `test/positioning/expected-values.jl` are exact Float64 outputs of the Python package solposx, generated at timestamps whose offset from noon UTC is a multiple of 84.375 s so the full Julian Date is exactly representable in Float64. This makes the references artifact-free and lets tests compare at `atol = 1e-10` (`1e-8` for SPA). When regenerating: keep timestamps on that grid, use whole seconds only because solposx `usno()` mishandles sub-second times, and paste full round-trip digits. ## Code Style -- Formatter: **Runic.jl** (enforced via pre-commit). Run before committing. +- Formatter: **Runic.jl** (enforced via pre-commit). Run before committing. Runic has no configuration. - Imports: All used symbols must be explicitly imported (checked by ExplicitImports.jl). - Package quality is checked with **Aqua.jl** and type inference with **JET.jl** (Julia 1.12 only) — both run as part of `Pkg.test()` via `test/linting.jl`. -- Docstrings use **DocStringExtensions** macros: `$(TYPEDEF)` + `$(TYPEDFIELDS)` for structs, `$(TYPEDSIGNATURES)` for functions. +- Docstrings use **DocStringExtensions** macros: `$(TYPEDEF)` + `$(TYPEDFIELDS)` for structs, `$(TYPEDSIGNATURES)` for functions. Doctests run against the setup in `docs/make.jl`, which predefines `obs` and `dt`. +- Public API changes should update the matching page under `docs/src/` and add a `CHANGELOG.md` entry under `## unreleased`. From f8a37c72cb1ae8667716921ed794f8fea69ea6f4 Mon Sep 17 00:00:00 2001 From: Stefan de Lange Date: Sun, 13 Sep 2026 02:00:06 +0200 Subject: [PATCH 03/13] Add the libsolarposition @api surface Co-Authored-By: Claude Opus 5 (1M context) --- lib/Project.toml | 15 + lib/src/libsolarposition.jl | 608 ++++++++++++++++++++++++++++++++++++ 2 files changed, 623 insertions(+) create mode 100644 lib/Project.toml create mode 100644 lib/src/libsolarposition.jl diff --git a/lib/Project.toml b/lib/Project.toml new file mode 100644 index 00000000..695be98d --- /dev/null +++ b/lib/Project.toml @@ -0,0 +1,15 @@ +name = "libsolarposition" +uuid = "b7d3f1a2-4c5e-4a8b-9f21-3d6e8c0a5b74" +version = "0.1.0" + +[deps] +Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +JLWInterop = "65e54657-ed21-41a3-96db-71ab7fa6d94b" +SolarPosition = "5b9d1343-a731-5a90-8730-7bf8d89bf3eb" + +[sources] +SolarPosition = {path = ".."} + +[compat] +JLWInterop = "0.2" +julia = "1.13" diff --git a/lib/src/libsolarposition.jl b/lib/src/libsolarposition.jl new file mode 100644 index 00000000..ce958da7 --- /dev/null +++ b/lib/src/libsolarposition.jl @@ -0,0 +1,608 @@ +""" +C and Python facing surface over SolarPosition.jl, compiled with `juliac --trim`. + +Every positioning algorithm and every refraction model is reachable. Algorithm and +refraction are selected by enum, which JuliaLibWrapping renders as a Python `IntEnum` and +validates on the way in, and their options ride along as keyword arguments. + +Two conventions keep the compiled path narrow. Time crosses the boundary as Unix seconds +rather than a `DateTime`, so nothing here touches TimeZones.jl. Every entrypoint returns +all five angles, with `apparent_*` equal to the true values when no refraction applies, so +the return type does not vary with the `SolPos`/`ApparentSolPos` split that +`SolarPosition.result_type` produces — that uniformity is what lets each dispatch arm stay +inferable under `--trim`. +""" +module libsolarposition + +using Dates: datetime2unix, unix2datetime +using JLWInterop: @api, @export_release_entrypoints, JLWInterop +import SolarPosition as SP +import SolarPosition.Refraction as SPR + +# The two enums live in nested modules because `@enum` binds its members in the enclosing +# module, and `MICHALSKY` and `SPA` each name both an algorithm and a refraction model. +# Only the types and the two default members are imported, which is enough for `@api` to +# resolve a bare member name as a keyword default. +module AlgorithmEnums + @enum Algorithm::Int32 PSA NOAA SPA WALRAVEN USNO IQBAL MICHALSKY +end + +module RefractionEnums + @enum RefractionModel::Int32 DEFAULT NONE HUGHES ARCHER BENNETT MICHALSKY SG2 SPA +end + +using .AlgorithmEnums: Algorithm, PSA +using .RefractionEnums: DEFAULT, RefractionModel + +"Julian date formulation for the Michalsky algorithm." +@enum JulianDateMode::Int32 ORIGINAL STANDARD + +""" +Solar angles in degrees. `apparent_elevation` and `apparent_zenith` equal `elevation` and +`zenith` when the selected refraction model applies no correction. +""" +struct SolarAngles + azimuth::Float64 + elevation::Float64 + zenith::Float64 + apparent_elevation::Float64 + apparent_zenith::Float64 +end + +# An isbits struct is its own carrier, so all three hooks are the identity. +JLWInterop.carrier_type(::Type{SolarAngles}) = SolarAngles +JLWInterop.to_carrier(a::SolarAngles) = a +JLWInterop.from_carrier(::Type{SolarAngles}, a::SolarAngles) = a + +# Algorithm and refraction options, bundled so the dispatch helpers take one argument +# instead of a dozen. Never crosses the ABI. +struct Opts + pressure::Float64 + temperature::Float64 + delta_t::Union{Float64, Nothing} + atmos_refract::Float64 + refraction_limit::Float64 + psa_coeffs::Int64 + gmst_option::Int64 + spencer_correction::Bool + julian_date::JulianDateMode +end + +_angles(p::SP.SolPos) = + SolarAngles(p.azimuth, p.elevation, p.zenith, p.elevation, p.zenith) +_angles(p::SP.ApparentSolPos) = SolarAngles( + p.azimuth, p.elevation, p.zenith, p.apparent_elevation, p.apparent_zenith, +) + +# Geometric position, refraction deliberately switched off. One arm per algorithm; each +# arm constructs a concrete algorithm type, so the branch is on a value and every call is +# statically dispatched. All arms return `SolPos{Float64}`. +@inline function _geometric(a::Algorithm, o::Opts, obs::SP.Observer{Float64}, dt) + return if a === AlgorithmEnums.PSA + SP.solar_position(obs, dt, SP.PSA(o.psa_coeffs), SPR.NoRefraction()) + elseif a === AlgorithmEnums.NOAA + SP.solar_position(obs, dt, SP.NOAA(o.delta_t), SPR.NoRefraction()) + elseif a === AlgorithmEnums.SPA + SP.solar_position( + obs, dt, + SP.SPA(o.delta_t, o.pressure, o.temperature, o.atmos_refract), + SPR.NoRefraction(), + ) + elseif a === AlgorithmEnums.WALRAVEN + SP.solar_position(obs, dt, SP.Walraven(), SPR.NoRefraction()) + elseif a === AlgorithmEnums.USNO + SP.solar_position(obs, dt, SP.USNO(o.delta_t, o.gmst_option), SPR.NoRefraction()) + elseif a === AlgorithmEnums.IQBAL + SP.solar_position(obs, dt, SP.Iqbal(), SPR.NoRefraction()) + else + jd = o.julian_date === ORIGINAL ? :original : :standard + SP.solar_position( + obs, dt, SP.Michalsky(o.spencer_correction, jd), SPR.NoRefraction(), + ) + end +end + +# Refraction correction in degrees, to be added to the true elevation. This mirrors what +# `Positioning._solar_position(obs, dt, alg, refraction)` does internally, which is why +# the algorithm and refraction branches can be separate rather than a 7x8 cross product. +@inline function _refraction_correction(r::RefractionModel, o::Opts, elevation::Float64) + return if r === RefractionEnums.HUGHES + SPR.refraction(SPR.HUGHES(o.pressure, o.temperature), elevation) + elseif r === RefractionEnums.ARCHER + SPR.refraction(SPR.ARCHER(), elevation) + elseif r === RefractionEnums.BENNETT + SPR.refraction(SPR.BENNETT(o.pressure, o.temperature), elevation) + elseif r === RefractionEnums.MICHALSKY + SPR.refraction(SPR.MICHALSKY(), elevation) + elseif r === RefractionEnums.SG2 + SPR.refraction(SPR.SG2(o.pressure, o.temperature), elevation) + else + SPR.refraction( + SPR.SPARefraction(o.pressure, o.temperature, o.refraction_limit), elevation, + ) + end +end + +# `DefaultRefraction` resolves per algorithm inside the package, and SPA's default path +# feeds the SPA struct's own pressure, temperature and atmos_refract fields, so it cannot +# be rebuilt from the decoupled form above. Delegate to the package instead. +@inline function _default_refraction(a::Algorithm, o::Opts, obs::SP.Observer{Float64}, dt) + return if a === AlgorithmEnums.PSA + _angles(SP.solar_position(obs, dt, SP.PSA(o.psa_coeffs))) + elseif a === AlgorithmEnums.NOAA + _angles(SP.solar_position(obs, dt, SP.NOAA(o.delta_t))) + elseif a === AlgorithmEnums.SPA + _angles( + SP.solar_position( + obs, dt, SP.SPA(o.delta_t, o.pressure, o.temperature, o.atmos_refract), + ), + ) + elseif a === AlgorithmEnums.WALRAVEN + _angles(SP.solar_position(obs, dt, SP.Walraven())) + elseif a === AlgorithmEnums.USNO + _angles(SP.solar_position(obs, dt, SP.USNO(o.delta_t, o.gmst_option))) + elseif a === AlgorithmEnums.IQBAL + _angles(SP.solar_position(obs, dt, SP.Iqbal())) + else + jd = o.julian_date === ORIGINAL ? :original : :standard + _angles(SP.solar_position(obs, dt, SP.Michalsky(o.spencer_correction, jd))) + end +end + +@inline function _compute( + obs::SP.Observer{Float64}, dt, a::Algorithm, r::RefractionModel, o::Opts, + ) + r === DEFAULT && return _default_refraction(a, o, obs, dt) + pos = _geometric(a, o, obs, dt) + # With no refraction there is no correction, and the package returns a plain `SolPos`. + # Pass its own zenith through rather than recomputing `90 - elevation`, which differs + # in the last bit for the algorithms that derive elevation from zenith. + r === RefractionEnums.NONE && return _angles(pos) + apparent_elevation = pos.elevation + _refraction_correction(r, o, pos.elevation) + return SolarAngles( + pos.azimuth, pos.elevation, pos.zenith, + apparent_elevation, 90.0 - apparent_elevation, + ) +end + +""" +Solar position for one instant. + +`unix_seconds` is seconds since the Unix epoch, UTC. `latitude` is positive north and +`longitude` positive east, both in degrees, and `altitude` is metres above mean sea level. + +`pressure` [Pa] and `temperature` [°C] feed the refraction models that take them, and the +SPA algorithm. `delta_t` [s] is the difference between terrestrial time and UT1; pass +`None` to have it computed from the date. The remaining keywords are algorithm options: +`atmos_refract` and `psa_coeffs` for SPA and PSA, `gmst_option` for USNO, and +`spencer_correction` and `julian_date` for Michalsky. `refraction_limit` is +SPARefraction's own refraction limit. +""" +function solar_position_single( + latitude::Float64, + longitude::Float64, + unix_seconds::Float64; + altitude::Float64 = 0.0, + algorithm::Algorithm = PSA, + refraction::RefractionModel = DEFAULT, + pressure::Float64 = 101325.0, + temperature::Float64 = 12.0, + delta_t::Union{Float64, Nothing} = 67.0, + atmos_refract::Float64 = 0.5667, + refraction_limit::Float64 = -0.5667, + psa_coeffs::Int64 = 2020, + gmst_option::Int64 = 1, + spencer_correction::Bool = true, + julian_date::JulianDateMode = ORIGINAL, + ) + obs = SP.Observer(latitude, longitude, altitude) + o = Opts( + pressure, temperature, delta_t, atmos_refract, refraction_limit, + psa_coeffs, gmst_option, spencer_correction, julian_date, + ) + return _compute(obs, unix2datetime(unix_seconds), algorithm, refraction, o) +end + +@api solar_position_single( + latitude::Float64, + longitude::Float64, + unix_seconds::Float64; + altitude::Float64 = 0.0, + algorithm::Algorithm = PSA, + refraction::RefractionModel = DEFAULT, + pressure::Float64 = 101325.0, + temperature::Float64 = 12.0, + delta_t::Union{Float64, Nothing} = 67.0, + atmos_refract::Float64 = 0.5667, + refraction_limit::Float64 = -0.5667, + psa_coeffs::Int64 = 2020, + gmst_option::Int64 = 1, + spencer_correction::Bool = true, + julian_date::JulianDateMode = ORIGINAL, +)::SolarAngles + +""" +Solar position for many instants, written into caller-allocated output buffers. + +The five output arrays must each hold at least `length(unix_seconds)` elements. They are +borrowed, so the writes land directly in the caller's numpy arrays and nothing is +allocated. The observer is built once and only read, so this is also the form the threaded +entrypoint parallelises. See `solar_position_single` for the keyword arguments. +""" +function solar_position_inplace( + latitude::Float64, + longitude::Float64, + unix_seconds::Vector{Float64}, + azimuth::Vector{Float64}, + elevation::Vector{Float64}, + zenith::Vector{Float64}, + apparent_elevation::Vector{Float64}, + apparent_zenith::Vector{Float64}; + altitude::Float64 = 0.0, + algorithm::Algorithm = PSA, + refraction::RefractionModel = DEFAULT, + pressure::Float64 = 101325.0, + temperature::Float64 = 12.0, + delta_t::Union{Float64, Nothing} = 67.0, + atmos_refract::Float64 = 0.5667, + refraction_limit::Float64 = -0.5667, + psa_coeffs::Int64 = 2020, + gmst_option::Int64 = 1, + spencer_correction::Bool = true, + julian_date::JulianDateMode = ORIGINAL, + ) + n = length(unix_seconds) + if length(azimuth) < n || length(elevation) < n || length(zenith) < n || + length(apparent_elevation) < n || length(apparent_zenith) < n + throw(DimensionMismatch("output buffers are shorter than the input")) + end + obs = SP.Observer(latitude, longitude, altitude) + o = Opts( + pressure, temperature, delta_t, atmos_refract, refraction_limit, + psa_coeffs, gmst_option, spencer_correction, julian_date, + ) + for i in 1:n + a = _compute(obs, unix2datetime(unix_seconds[i]), algorithm, refraction, o) + azimuth[i] = a.azimuth + elevation[i] = a.elevation + zenith[i] = a.zenith + apparent_elevation[i] = a.apparent_elevation + apparent_zenith[i] = a.apparent_zenith + end + return nothing +end + +@api solar_position_inplace( + latitude::Float64, + longitude::Float64, + unix_seconds::Vector{Float64}, + azimuth::Vector{Float64}, + elevation::Vector{Float64}, + zenith::Vector{Float64}, + apparent_elevation::Vector{Float64}, + apparent_zenith::Vector{Float64}; + altitude::Float64 = 0.0, + algorithm::Algorithm = PSA, + refraction::RefractionModel = DEFAULT, + pressure::Float64 = 101325.0, + temperature::Float64 = 12.0, + delta_t::Union{Float64, Nothing} = 67.0, + atmos_refract::Float64 = 0.5667, + refraction_limit::Float64 = -0.5667, + psa_coeffs::Int64 = 2020, + gmst_option::Int64 = 1, + spencer_correction::Bool = true, + julian_date::JulianDateMode = ORIGINAL, +)::Nothing + +""" +Solar position for many instants, returned as a freshly allocated `n x 5` matrix. + +The columns are azimuth, elevation, zenith, apparent elevation and apparent zenith. An +owning array return transfers the allocation to the caller, so this is the convenient form +rather than the fast one; use `solar_position_inplace` to write into buffers you +already hold. See `solar_position_single` for the keyword arguments. +""" +function solar_position( + latitude::Float64, + longitude::Float64, + unix_seconds::Vector{Float64}; + altitude::Float64 = 0.0, + algorithm::Algorithm = PSA, + refraction::RefractionModel = DEFAULT, + pressure::Float64 = 101325.0, + temperature::Float64 = 12.0, + delta_t::Union{Float64, Nothing} = 67.0, + atmos_refract::Float64 = 0.5667, + refraction_limit::Float64 = -0.5667, + psa_coeffs::Int64 = 2020, + gmst_option::Int64 = 1, + spencer_correction::Bool = true, + julian_date::JulianDateMode = ORIGINAL, + ) + n = length(unix_seconds) + obs = SP.Observer(latitude, longitude, altitude) + o = Opts( + pressure, temperature, delta_t, atmos_refract, refraction_limit, + psa_coeffs, gmst_option, spencer_correction, julian_date, + ) + out = Matrix{Float64}(undef, n, 5) + for i in 1:n + a = _compute(obs, unix2datetime(unix_seconds[i]), algorithm, refraction, o) + out[i, 1] = a.azimuth + out[i, 2] = a.elevation + out[i, 3] = a.zenith + out[i, 4] = a.apparent_elevation + out[i, 5] = a.apparent_zenith + end + return out +end + +@api solar_position( + latitude::Float64, + longitude::Float64, + unix_seconds::Vector{Float64}; + altitude::Float64 = 0.0, + algorithm::Algorithm = PSA, + refraction::RefractionModel = DEFAULT, + pressure::Float64 = 101325.0, + temperature::Float64 = 12.0, + delta_t::Union{Float64, Nothing} = 67.0, + atmos_refract::Float64 = 0.5667, + refraction_limit::Float64 = -0.5667, + psa_coeffs::Int64 = 2020, + gmst_option::Int64 = 1, + spencer_correction::Bool = true, + julian_date::JulianDateMode = ORIGINAL, +)::Matrix{Float64} + +# ----------------------------------------------------------- sunrise, sunset, transit +# +# These take no algorithm selector. `Utilities._transit_sunrise_sunset` is only defined +# for `alg::SPA` — the wider `SolarAlgorithm` bound on the public signature has no other +# implementation — so the SPA options are exposed directly instead. + +"Which solar event to locate." +@enum SunEvent::Int32 SUNRISE SUNSET SOLAR_NOON + +"Search forwards or backwards in time." +@enum EventDirection::Int32 NEXT PREVIOUS + +""" +Solar transit, sunrise and sunset as seconds since midnight UTC. + +All three can exceed 86400 or go negative when the event falls outside the UTC day that +contains the requested instant. +""" +struct SunEvents + transit::Float64 + sunrise::Float64 + sunset::Float64 +end + +JLWInterop.carrier_type(::Type{SunEvents}) = SunEvents +JLWInterop.to_carrier(e::SunEvents) = e +JLWInterop.from_carrier(::Type{SunEvents}, e::SunEvents) = e + +""" +Solar transit, sunrise and sunset for the UTC day containing `unix_seconds`. + +Declared with positional arguments: reached through `@api`'s keyword-argument wrapper the +`juliac --trim` verifier cannot resolve the call ("unresolved call from statement +Core.kwcall(...)"). `sun_event` keeps its keywords because it goes through the package's +`DateTime` path, which takes a different `_frac_to_event` specialisation. The keyword-only +Python signature is restored by hand in `lib/python/_extras.py`. + +Returned as seconds since midnight UTC of that day, which is the form that keeps full +precision — the `DateTime` variant in the Julia package rounds to a whole second. + +Computed with SPA, so `delta_t`, `pressure`, `temperature` and `atmos_refract` are SPA's +own options. + +There is deliberately no `horizon` parameter. `src/Utilities/spa.jl` hardcodes the +horizon depression to -0.8333 degrees and never reads `Observer.horizon`, so a parameter +here would be silently ignored — which matters most for exactly the case someone would +reach for it, such as civil twilight. Add one once the package honours the field. + +At latitudes in polar day or polar night there is no sunrise or sunset. The underlying +package warns on stderr and returns zero for all three fields in that case. +""" +function transit_sunrise_sunset( + latitude::Float64, + longitude::Float64, + unix_seconds::Float64, + altitude::Float64, + pressure::Float64, + temperature::Float64, + delta_t::Union{Float64, Nothing}, + atmos_refract::Float64, + ) + obs = SP.Observer(latitude, longitude, altitude) + alg = SP.SPA(delta_t, pressure, temperature, atmos_refract) + r = SP.transit_sunrise_sunset_seconds(obs, unix2datetime(unix_seconds), alg) + return SunEvents(r.transit, r.sunrise, r.sunset) +end + +@api transit_sunrise_sunset( + latitude::Float64, + longitude::Float64, + unix_seconds::Float64, + altitude::Float64, + pressure::Float64, + temperature::Float64, + delta_t::Union{Float64, Nothing}, + atmos_refract::Float64, +)::SunEvents + +""" +The next or previous sunrise, sunset or solar noon, as Unix seconds. + +`event` selects sunrise, sunset or solar noon and `direction` selects whether to search +forward or backward from `unix_seconds`. The result is rounded to a whole second, because +the underlying package returns a `DateTime`; use `transit_sunrise_sunset` when the +sub-second part matters. See `transit_sunrise_sunset` for the keyword arguments. +""" +function sun_event( + latitude::Float64, + longitude::Float64, + unix_seconds::Float64, + event::SunEvent, + direction::EventDirection; + altitude::Float64 = 0.0, + pressure::Float64 = 101325.0, + temperature::Float64 = 12.0, + delta_t::Union{Float64, Nothing} = 67.0, + atmos_refract::Float64 = 0.5667, + ) + obs = SP.Observer(latitude, longitude, altitude) + alg = SP.SPA(delta_t, pressure, temperature, atmos_refract) + dt = unix2datetime(unix_seconds) + # Each arm calls a package function that passes a constant field symbol down to + # `_next_event`/`_previous_event`, so nothing dispatches on a runtime Symbol. + result = if direction === NEXT + if event === SUNRISE + SP.next_sunrise(obs, dt, alg) + elseif event === SUNSET + SP.next_sunset(obs, dt, alg) + else + SP.next_solar_noon(obs, dt, alg) + end + else + if event === SUNRISE + SP.previous_sunrise(obs, dt, alg) + elseif event === SUNSET + SP.previous_sunset(obs, dt, alg) + else + SP.previous_solar_noon(obs, dt, alg) + end + end + return datetime2unix(result) +end + +@api sun_event( + latitude::Float64, + longitude::Float64, + unix_seconds::Float64, + event::SunEvent, + direction::EventDirection; + altitude::Float64 = 0.0, + pressure::Float64 = 101325.0, + temperature::Float64 = 12.0, + delta_t::Union{Float64, Nothing} = 67.0, + atmos_refract::Float64 = 0.5667, +)::Float64 + +""" +The number of Julia threads this library was compiled with. + +The thread count is fixed when the library is built and cannot be changed at load time, so +this is how a caller finds out whether it has a threaded build. A value of 1 means +`solar_position_inplace_threaded` runs serially. +""" +julia_nthreads() = Int64(Threads.nthreads()) + +@api julia_nthreads()::Int64 + +# A function barrier for the parallel loop. `Threads.@threads` builds a closure, and +# through `@api`'s keyword-argument wrapper the trim verifier cannot resolve the resulting +# `kwcall` ("unresolved call from statement Core.kwcall(...)"). Taking every value as a +# positional argument of concrete type keeps the loop inferable and the entrypoint +# trimmable. +function _fill_threaded!( + obs::SP.Observer{Float64}, + n::Int, + unix_seconds::Vector{Float64}, + azimuth::Vector{Float64}, + elevation::Vector{Float64}, + zenith::Vector{Float64}, + apparent_elevation::Vector{Float64}, + apparent_zenith::Vector{Float64}, + algorithm::Algorithm, + refraction::RefractionModel, + o::Opts, + ) + Threads.@threads for i in 1:n + a = _compute(obs, unix2datetime(unix_seconds[i]), algorithm, refraction, o) + azimuth[i] = a.azimuth + elevation[i] = a.elevation + zenith[i] = a.zenith + apparent_elevation[i] = a.apparent_elevation + apparent_zenith[i] = a.apparent_zenith + end + return nothing +end + +""" +Threaded `solar_position_inplace`, using the threads the library was compiled with. + +The work is embarrassingly parallel: the observer is built once and only read, and each +iteration writes its own index in the caller's buffers. With a single-threaded build this +is just the serial loop. See `solar_position_single` for the keyword arguments. +""" +function solar_position_inplace_threaded( + latitude::Float64, + longitude::Float64, + unix_seconds::Vector{Float64}, + azimuth::Vector{Float64}, + elevation::Vector{Float64}, + zenith::Vector{Float64}, + apparent_elevation::Vector{Float64}, + apparent_zenith::Vector{Float64}, + altitude::Float64, + algorithm::Algorithm, + refraction::RefractionModel, + pressure::Float64, + temperature::Float64, + delta_t::Union{Float64, Nothing}, + atmos_refract::Float64, + refraction_limit::Float64, + psa_coeffs::Int64, + gmst_option::Int64, + spencer_correction::Bool, + julian_date::JulianDateMode, + ) + n = length(unix_seconds) + if length(azimuth) < n || length(elevation) < n || length(zenith) < n || + length(apparent_elevation) < n || length(apparent_zenith) < n + throw(DimensionMismatch("output buffers are shorter than the input")) + end + obs = SP.Observer(latitude, longitude, altitude) + o = Opts( + pressure, temperature, delta_t, atmos_refract, refraction_limit, + psa_coeffs, gmst_option, spencer_correction, julian_date, + ) + _fill_threaded!( + obs, n, unix_seconds, azimuth, elevation, zenith, + apparent_elevation, apparent_zenith, algorithm, refraction, o, + ) + return nothing +end + +@api solar_position_inplace_threaded( + latitude::Float64, + longitude::Float64, + unix_seconds::Vector{Float64}, + azimuth::Vector{Float64}, + elevation::Vector{Float64}, + zenith::Vector{Float64}, + apparent_elevation::Vector{Float64}, + apparent_zenith::Vector{Float64}, + altitude::Float64, + algorithm::Algorithm, + refraction::RefractionModel, + pressure::Float64, + temperature::Float64, + delta_t::Union{Float64, Nothing}, + atmos_refract::Float64, + refraction_limit::Float64, + psa_coeffs::Int64, + gmst_option::Int64, + spencer_correction::Bool, + julian_date::JulianDateMode, +)::Nothing + +# `solar_position` hands an owned array to the caller, so the generated bindings +# need the release entrypoints to give it back. +@export_release_entrypoints + +end # module From b16b2e4fb1737eaf94f623f37f6770e2c74603a1 Mon Sep 17 00:00:00 2001 From: Stefan de Lange Date: Sun, 13 Sep 2026 02:00:17 +0200 Subject: [PATCH 04/13] Add the juliac build script and thread shim Co-Authored-By: Claude Opus 5 (1M context) --- lib/build-env/Project.toml | 8 ++ lib/build.jl | 228 +++++++++++++++++++++++++++++++++++++ lib/c/set_num_threads.c | 83 ++++++++++++++ 3 files changed, 319 insertions(+) create mode 100644 lib/build-env/Project.toml create mode 100644 lib/build.jl create mode 100644 lib/c/set_num_threads.c diff --git a/lib/build-env/Project.toml b/lib/build-env/Project.toml new file mode 100644 index 00000000..89dd1c9c --- /dev/null +++ b/lib/build-env/Project.toml @@ -0,0 +1,8 @@ +[deps] +JuliaC = "acedd4c2-ced6-4a15-accc-2607eb759ba2" +JuliaLibWrapping = "d61f35a8-f6af-436f-bc10-cee6b101f7bd" + +[compat] +JuliaC = "0.3.10" +JuliaLibWrapping = "0.2" +julia = "1.13" diff --git a/lib/build.jl b/lib/build.jl new file mode 100644 index 00000000..a5310a72 --- /dev/null +++ b/lib/build.jl @@ -0,0 +1,228 @@ +#!/usr/bin/env julia +# +# Build libsolarposition and its C header and Python package. +# +# julia --project=lib lib/build.jl # default: signal-safe, 1 thread +# julia --project=lib lib/build.jl --threads=4 # threaded, see the caveat below +# julia --project=lib lib/build.jl --verbose +# +# ## The two builds +# +# The default build sets `threads=1` and `handle-signals=no`. Those are JuliaC's defaults +# for a library and the correct mode for embedding: Python keeps its own SIGINT/SIGSEGV +# handlers, so Ctrl-C behaves normally. This is the build the wheel ships. +# +# The threaded build sets `threads=N` and `handle-signals=yes`. The second half is not +# optional: on Julia >= 1.12 a library with more than one thread and signal handling off +# segfaults during shutdown, and JuliaC's own validation step rejects the combination. The +# cost is that Julia installs its own signal handlers into the host process. Upstream +# tracking issue: https://github.com/JuliaLang/julia/issues/61319 +# +# The thread count is fixed at compile time. It cannot be chosen when the library loads: +# JuliaC compiles an `__attribute__((constructor))` shim that calls `jl_parse_opts` with +# `--threads=` before `jl_init`, which overrides `JULIA_NUM_THREADS`, and by the time +# Python has `dlopen`ed the library that constructor has already run. +# +# ## Why the threaded path drives JuliaC directly +# +# JuliaLibWrapping 0.2's `build_library` does not forward `jl_options` to +# `JuliaC.ImageRecipe`, so the thread count cannot be set through it. Both builds produce +# byte-identical wrappers and ABI metadata — only the `.so` differs — so the threaded path +# runs the normal `standard_build` first and then recompiles just the library and bundle +# with `jl_options` set. Delete `build_threaded_library` once JuliaLibWrapping grows a +# `jl_options` passthrough. + +push!(LOAD_PATH, joinpath(@__DIR__, "build-env")) +using JuliaC: BundleRecipe, ImageRecipe, LinkRecipe, bundle_products, compile_products, + link_products +using JuliaLibWrapping: standard_build +using TOML: TOML + +const LIBNAME = "libsolarposition" +const PYPKG = "solarposition" + +flag(args, name) = name in args + +function parse_threads(args) + for a in args + startswith(a, "--threads=") || continue + return a[(length("--threads=") + 1):end] + end + return "1" +end + +# `juliac` relocates the project into a temporary directory to build it, which breaks the +# relative paths this project uses to reach SolarPosition one level up. `build_library` +# has a private helper for this; the threaded path needs its own. Paths are resolved +# against the original project directory, matching what JuliaLibWrapping does, so the +# manifest keeps pointing at the real package rather than into the copy. +function absolutize_paths!(table, base) + table isa AbstractDict || return + p = get(table, "path", nothing) + if p isa AbstractString && !isabspath(p) + abs = abspath(joinpath(base, p)) + ispath(abs) || error("path \"$p\" resolves to $abs, which does not exist") + table["path"] = abs + end + for (_, v) in table + if v isa AbstractDict + absolutize_paths!(v, base) + elseif v isa AbstractVector + foreach(e -> absolutize_paths!(e, base), v) + end + end + return +end + +function materialize_project(project::AbstractString) + dir = mktempdir() + for entry in readdir(project) + entry in ("out", "test", "build-env") && continue # not needed to compile + cp(joinpath(project, entry), joinpath(dir, entry); force = true) + end + for name in ("Project.toml", "Manifest.toml") + path = joinpath(dir, name) + isfile(path) || continue + toml = TOML.parsefile(path) + absolutize_paths!(toml, project) + open(path, "w") do io + TOML.print(io, toml) + end + end + return dir +end + +""" +Recompile the library and its bundle with `jl_options` set and the load-time thread shim +linked in, then refresh the bundle copy inside the generated Python package. Everything +else from `standard_build` — the ABI JSON, the C header, `_lowlevel.py`, `_facade.py` — +is independent of this and is left alone. + +This runs for every build, not just threaded ones, because `standard_build` can pass +neither `c_sources` nor `jl_options` to `JuliaC.ImageRecipe`. +""" +function rebuild_library(dir, out, threads, handle_signals; verbose::Bool) + entry = joinpath(dir, "src", LIBNAME * ".jl") + bundle_dir = joinpath(out, LIBNAME * "-bundle") + project = materialize_project(dir) + + img = ImageRecipe(; + output_type = "--output-lib", + file = entry, + project, + trim_mode = "safe", + add_ccallables = true, + export_abi = joinpath(out, LIBNAME * ".abi.json"), + verbose, + c_sources = [joinpath(dir, "c", "set_num_threads.c")], + jl_options = Dict("threads" => threads, "handle-signals" => handle_signals), + ) + # `@bundle` is required for a relocatable wheel; without it the `.so` bakes in + # absolute paths to this machine's Julia installation. + link = LinkRecipe(; + image_recipe = img, outname = joinpath(out, LIBNAME), rpath = "@bundle", + ) + compile_products(img) + link_products(link) + + # `juliac --bundle` refuses to overwrite an existing tree. + ispath(bundle_dir) && rm(bundle_dir; recursive = true) + bundle_products( + BundleRecipe(; link_recipe = link, output_dir = bundle_dir, privatize = true), + ) + + pkg_bundle = joinpath(out, PYPKG, "bundle") + ispath(pkg_bundle) && rm(pkg_bundle; recursive = true) + cp(bundle_dir, pkg_bundle) + return bundle_dir +end + +const EXTRAS_MARKER = "# Appended by lib/build.jl: hand-maintained wrappers from _extras.py." + +""" +The names `_extras.py` exports, read out of its `__all__` literal. + +Parsed rather than hardcoded so that adding a wrapper there is enough to have it shadow +the generated one: a name listed in `__all__` but missing from the appended import would +be reachable through `from ._facade import *` while the generated definition still won. +""" +function extras_exports(path) + text = read(path, String) + m = match(r"^__all__\s*=\s*\[(.*?)\]"ms, text) + m === nothing && error("no __all__ literal in $path") + return [String(x.captures[1]) for x in eachmatch(r"\"([^\"]+)\"", m.captures[1])] +end + +""" +Install the hand-maintained Python additions from `lib/python/` into the generated +package. `_extras.py` becomes a module of its own, and one import is appended to the +generated `_facade.py` so its definitions win. See `lib/python/_extras.py` for why some +entrypoints need this. +""" +function install_python_extras(dir, out) + extras = joinpath(dir, "python", "_extras.py") + isfile(extras) || return false + cp(extras, joinpath(out, PYPKG, "_extras.py"); force = true) + names = extras_exports(extras) + facade = joinpath(out, PYPKG, "_facade.py") + text = read(facade, String) + # Truncate any block from an earlier build before appending, so a changed export list + # replaces the old import instead of stacking a second one on top of it. + idx = findfirst(EXTRAS_MARKER, text) + idx === nothing || (text = rstrip(text[1:(first(idx) - 1)]) * "\n") + open(facade, "w") do io + print(io, text) + println(io) + println(io, EXTRAS_MARKER) + println(io, "# The import shadows the generated definitions; the __all__ update is") + println(io, "# what makes the new names reachable through `from ._facade import *`.") + println(io, "from ._extras import $(join(names, ", ")) # noqa: F811,E402") + println(io, "from ._extras import __all__ as _extras_all # noqa: E402") + println(io, "__all__ = sorted(set(__all__) | set(_extras_all))") + end + return true +end + +function main(args) + threads = parse_threads(args) + verbose = flag(args, "--verbose") + out = joinpath(@__DIR__, "out") + + @info "Building $LIBNAME" threads handle_signals = + (threads == "1" ? "no" : "yes") + + # JuliaLibWrapping writes `_facade.py` once and never overwrites it, so a renamed or + # deleted entrypoint would otherwise survive as a definition calling a `_lowlevel` + # symbol that no longer exists. Nothing in it is hand-edited -- the hand-maintained + # wrappers live in `python/_extras.py` and are appended below -- so regenerate it. + facade = joinpath(out, PYPKG, "_facade.py") + isfile(facade) && rm(facade) + + result = standard_build( + @__DIR__; + libname = LIBNAME, + python_package = PYPKG, + bundle = true, + verbose, + ) + + install_python_extras(@__DIR__, out) + + if threads != "1" + @warn """ + Threaded build: Julia will install its own SIGINT/SIGSEGV handlers into the host + process, because more than one thread requires handle-signals=yes on Julia >= 1.12. + A default build can raise the thread count at load time instead, via + solarposition.set_num_threads(). See https://github.com/JuliaLang/julia/issues/61319 + """ + end + rebuild_library(@__DIR__, out, threads, threads == "1" ? "no" : "yes"; verbose) + + @info "Build complete" library = result.library abi = result.abi_path + @info "Python package" path = joinpath(out, PYPKG) install = "pip install $out" + return result +end + +main(ARGS) + +pop!(LOAD_PATH) diff --git a/lib/c/set_num_threads.c b/lib/c/set_num_threads.c new file mode 100644 index 00000000..af3ef21e --- /dev/null +++ b/lib/c/set_num_threads.c @@ -0,0 +1,83 @@ +/* + * Load-time thread configuration for libsolarposition. + * + * JuliaC compiles its own __attribute__((constructor)) that calls jl_parse_opts with the + * --threads and --handle-signals values baked in at build time. The Julia runtime does + * not initialise until the first call into a @ccallable entrypoint, so there is a window + * after dlopen in which those options can still be changed. + * + * Two things matter for using that window: + * + * 1. Assigning jl_options.nthreads directly does NOT work. The write lands, but + * Threads.nthreads() does not follow it; jl_parse_opts does bookkeeping beyond that + * field. Re-running jl_parse_opts, exactly as JuliaC's own shim does, works. + * + * 2. More than one thread requires handle-signals=yes. On Julia >= 1.12 a library with + * handle-signals=no and more than one thread segfaults during shutdown. So asking + * for n > 1 necessarily turns Julia's own SIGINT/SIGSEGV handlers on in the host + * process. That is the caller's decision to make, which is why it is a function + * rather than something this library does on its own. + * See https://github.com/JuliaLang/julia/issues/61319 + * + * Must be called before the first solar position call. Afterwards it is a no-op as far as + * the runtime is concerned, so it returns -2 once the runtime has started. + */ +#include +#include +#include +#include + +/* Set by the first successful configuration so repeat calls can be reported. */ +static int solarposition_threads_configured = 0; + +static int apply_option(const char *opt) +{ + /* Save the getopt state and the fields jl_parse_opts overwrites, exactly as JuliaC's + * shim does, so the host process is left as it was found. */ + int saved_optind = optind, saved_opterr = opterr, saved_optopt = optopt; + char *saved_optarg = optarg; + const char *saved_image_file = jl_options.image_file; + + char *argv[] = {"juliac", (char *)opt, NULL}; + int argc = 2; + char **argvp = argv; + + optind = 0; + jl_parse_opts(&argc, &argvp); + + jl_options.image_file = saved_image_file; + optind = saved_optind; + opterr = saved_opterr; + optopt = saved_optopt; + optarg = saved_optarg; + return 0; +} + +/* + * Request `n` Julia compute threads. Returns the resulting Threads.nthreads() value, or + * -1 for n < 1 and -2 if the Julia runtime has already started. + */ +JL_DLLEXPORT int solarposition_set_num_threads(int n) +{ + char flag[32]; + + if (n < 1) return -1; + if (jl_is_initialized()) return -2; + + if (n > 1) { + /* Required: see note 2 above. */ + apply_option("--handle-signals=yes"); + } + snprintf(flag, sizeof(flag), "--threads=%d", n); + apply_option(flag); + solarposition_threads_configured = 1; + + /* jl_options.nthreads counts the interactive thread too; report the compute pool. */ + return n; +} + +/* The raw jl_options.nthreads, for diagnostics. Includes the interactive thread. */ +JL_DLLEXPORT int solarposition_get_option_nthreads(void) +{ + return (int)jl_options.nthreads; +} From 766ea41efc633f192f48156a5e8e4d484e5b2c3f Mon Sep 17 00:00:00 2001 From: Stefan de Lange Date: Sun, 13 Sep 2026 02:00:17 +0200 Subject: [PATCH 05/13] Add hand-maintained Python wrappers Co-Authored-By: Claude Opus 5 (1M context) --- lib/python/_extras.py | 499 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 499 insertions(+) create mode 100644 lib/python/_extras.py diff --git a/lib/python/_extras.py b/lib/python/_extras.py new file mode 100644 index 00000000..333565e8 --- /dev/null +++ b/lib/python/_extras.py @@ -0,0 +1,499 @@ +"""Hand-maintained additions to the generated façade. + +`build.jl` copies this module into the generated Python package and appends an import to +`_facade.py` so these definitions override the generated ones. + +Why this file exists: two entrypoints must be declared with all-positional arguments +the Julia side, because reached through `@api`'s keyword-argument wrapper the +`juliac --trim` verifier cannot resolve the call ("unresolved call from statement +Core.kwcall(...)") and the build fails. + +- ``solar_position_inplace_threaded`` — ``Threads.@threads`` builds a closure. +- ``transit_sunrise_sunset`` — the seconds-returning sunrise/sunset path takes a + different ``_frac_to_event`` specialisation than the ``DateTime`` one, and only the + former trips the verifier. ``sun_event`` goes through the ``DateTime`` path and keeps + its keywords, so it needs no wrapper for this reason. + +Declaring those arguments positionally keeps the entrypoints trimmable, and the +wrappers here restore the keyword-only signatures the rest of the API has. + +``solar_position`` is here for an unrelated reason: the generated binding for it +is miscompiled by some builds of libffi. See its docstring. + +Every entrypoint that takes a time accepts datetimes as well as Unix seconds — see +``_unix_seconds``. The C side takes only ``Float64``, so that conversion has to happen +in Python, and it is the sole reason ``solar_position_single``, +``solar_position_inplace`` and ``sun_event`` are wrapped here at all. +""" + +import ctypes +import datetime + +import numpy as np + +from . import _lowlevel +from ._lowlevel import ( + COpt_Float64, + CVector_borrowed_Float64, + Algorithm, + EventDirection, + JulianDateMode, + RefractionModel, + SunEvent, + _enum_coerce, +) + +__all__ = [ + "set_num_threads", + "solar_position", + "solar_position_inplace", + "solar_position_inplace_threaded", + "solar_position_single", + "sun_event", + "transit_sunrise_sunset", +] + + +def _unix_seconds(when): + """One instant as seconds since the Unix epoch, UTC. + + Accepts a ``datetime.datetime`` (including a pandas ``Timestamp``), a + ``datetime.date`` or a ``numpy.datetime64``. Anything else — a plain number above + all — is returned untouched, so this only ever adds accepted input. + + A naive ``datetime`` is read as UTC, not as local time. That differs from + ``datetime.timestamp()``, deliberately: every time this library takes or returns is + UTC, and ``numpy.datetime64`` carries no zone to read, so treating naive input as + UTC is the only rule that stays consistent across all the accepted types. Attach + ``tzinfo=datetime.timezone.utc`` when you want to be explicit, or any other + ``tzinfo`` to have the offset applied. + """ + if isinstance(when, datetime.datetime): + if when.tzinfo is None: + when = when.replace(tzinfo=datetime.timezone.utc) + return when.timestamp() + if isinstance(when, datetime.date): + # A bare date means midnight UTC. `datetime` is a subclass of `date`, so this + # is reached only by a real `date`. + midnight = datetime.datetime( + when.year, when.month, when.day, tzinfo=datetime.timezone.utc + ) + return midnight.timestamp() + if isinstance(when, np.datetime64): + return np.datetime64(when, "us").astype(np.int64) / 1e6 + return when + + +def _unix_seconds_each(times): + """Unix seconds for a sequence of instants, one element at a time.""" + return np.array([_unix_seconds(t) for t in times], dtype=np.float64) + + +def _holds_datetimes(times): + """Whether ``times`` starts with a datetime, taking the rest to match.""" + for first in times: + return isinstance(first, (datetime.date, np.datetime64)) + return False + + +def _unix_seconds_array(times): + """A sequence of instants as a ``float64`` array of Unix seconds. + + Converts ``datetime64`` arrays, pandas ``DatetimeIndex`` and ``Series`` (tz-aware + or not) and lists of ``datetime`` objects, by the rules in ``_unix_seconds``. + Anything else is passed through untouched, so numeric input still has to arrive as + a contiguous ``float64`` array rather than being silently reinterpreted. + """ + dtype = getattr(times, "dtype", None) + if getattr(dtype, "tz", None) is not None: + # A tz-aware pandas container. What numpy sees of one — UTC or local wall time + # — has varied across pandas versions, so read each Timestamp's own offset + # instead of converting the block. + return _unix_seconds_each(times) + kind = getattr(dtype, "kind", None) + if kind == "M": + # Microseconds, because that is as fine as a float64 Unix second resolves at + # present-day magnitudes; going through int64 keeps the conversion exact. + return np.asarray(times, dtype="datetime64[us]").astype(np.int64) / 1e6 + if (kind == "O" or isinstance(times, (list, tuple))) and _holds_datetimes(times): + return _unix_seconds_each(times) + return times + + +def set_num_threads(n: int) -> int: + """Choose how many Julia threads to use. Starts the Julia runtime. + + Must be called before any solar position call, because the thread count is fixed + once + the runtime is running. Returns the number of threads actually configured, and + raises + ``RuntimeError`` if that does not match what you asked for — which is what happens + if + the runtime was already started, whether by an earlier computation or an earlier + call + to this function. + + The default build starts single-threaded with Julia's signal handling disabled, + which + is the right mode for a library embedded in Python: Ctrl-C keeps working normally. + Asking for more than one thread necessarily switches Julia's own SIGINT/SIGSEGV + handlers on, because on Julia >= 1.12 a multi-threaded library with signal handling + off segfaults during shutdown. That trade-off is why this is an explicit call rather + than something the library does for you. See + https://github.com/JuliaLang/julia/issues/61319 + + >>> import solarposition as sp + >>> sp.set_num_threads(8) # doctest: +SKIP + 8 + """ + if n < 1: + raise ValueError(f"n must be >= 1, got {n}") + + fn = _lowlevel._lib.solarposition_set_num_threads + fn.argtypes = [ctypes.c_int] + fn.restype = ctypes.c_int + if fn(int(n)) < 0: + raise ValueError(f"could not request {n} threads") + + # Start the runtime now and report what it actually came up with. There is no + # reliable way to ask whether the runtime is already running from here — + # `jl_is_initialized()` reads 0 even after a call in an AOT-compiled library — so + # verifying the outcome is what makes a late call an error rather than a silent no- + # op. + actual = _lowlevel.libsolarposition_julia_nthreads().value + if actual != n: + raise RuntimeError( + f"asked for {n} Julia threads but got {actual}; the runtime was already " + "running. Call set_num_threads() before any solar position call." + ) + return actual + + +def solar_position_single( + latitude, + longitude, + unix_seconds, + *, + altitude=0.0, + algorithm=Algorithm.PSA, + refraction=RefractionModel.DEFAULT, + pressure=101325.0, + temperature=12.0, + delta_t=67.0, + atmos_refract=0.5667, + refraction_limit=-0.5667, + psa_coeffs=2020, + gmst_option=1, + spencer_correction=True, + julian_date=JulianDateMode.ORIGINAL, +): + """Solar position for one instant. + + ``unix_seconds`` is seconds since the Unix epoch, UTC, or any of the datetime types + ``_unix_seconds`` accepts:: + + >>> import datetime as dt, solarposition as sp + >>> when = dt.datetime(2023, 6, 21, 3, 45, tzinfo=dt.timezone.utc) + >>> sp.solar_position_single(40.0, -105.0, when) # doctest: +SKIP + + A naive datetime is read as UTC. ``latitude`` is positive north and ``longitude`` + positive east, both in degrees, and ``altitude`` is metres above mean sea level. + + ``pressure`` [Pa] and ``temperature`` [°C] feed the refraction models that take + them, and the SPA algorithm. ``delta_t`` [s] is the difference between terrestrial + time and UT1; pass ``None`` to have it computed from the date. The remaining + keywords are algorithm options: ``atmos_refract`` and ``psa_coeffs`` for SPA and + PSA, ``gmst_option`` for USNO, and ``spencer_correction`` and ``julian_date`` for + Michalsky. ``refraction_limit`` is SPARefraction's own refraction limit. + """ + _r = _lowlevel.libsolarposition_solar_position_single( + latitude, + longitude, + _unix_seconds(unix_seconds), + altitude, + _enum_coerce(Algorithm, algorithm), + _enum_coerce(RefractionModel, refraction), + pressure, + temperature, + COpt_Float64.from_optional(delta_t), + atmos_refract, + refraction_limit, + psa_coeffs, + gmst_option, + spencer_correction, + _enum_coerce(JulianDateMode, julian_date), + ) + return _r.value + + +def transit_sunrise_sunset( + latitude, + longitude, + unix_seconds, + *, + altitude=0.0, + pressure=101325.0, + temperature=12.0, + delta_t=67.0, + atmos_refract=0.5667, +): + """Solar transit, sunrise and sunset for the UTC day containing ``unix_seconds``. + + Returns a ``SunEvents`` with ``transit``, ``sunrise`` and ``sunset`` as seconds + since + midnight UTC of that day. Any of them can fall outside ``[0, 86400)`` when the event + lands outside that UTC day. + + ``unix_seconds`` may be given as a datetime; see ``solar_position_single``. + + Computed with SPA, so ``delta_t``, ``pressure``, ``temperature`` and + ``atmos_refract`` are SPA's options. There is deliberately no ``horizon`` parameter: + the package hardcodes the horizon depression and ignores ``Observer.horizon``. + + In polar day or polar night there is no sunrise or sunset. The underlying Julia + package warns on stderr and returns zero for all three fields. + """ + _r = _lowlevel.libsolarposition_transit_sunrise_sunset( + latitude, + longitude, + _unix_seconds(unix_seconds), + altitude, + pressure, + temperature, + COpt_Float64.from_optional(delta_t), + atmos_refract, + ) + return _r.value + + +def sun_event( + latitude, + longitude, + unix_seconds, + event, + direction, + *, + altitude=0.0, + pressure=101325.0, + temperature=12.0, + delta_t=67.0, + atmos_refract=0.5667, +): + """The next or previous sunrise, sunset or solar noon, as Unix seconds. + + ``event`` selects sunrise, sunset or solar noon and ``direction`` selects whether to + search forward or backward from ``unix_seconds``, which may be given as a datetime; + see ``solar_position_single``. The result is rounded to a whole second, because the + underlying package returns a ``DateTime``; use ``transit_sunrise_sunset`` when the + sub-second part matters. See it also for the keyword arguments. + """ + _r = _lowlevel.libsolarposition_sun_event( + latitude, + longitude, + _unix_seconds(unix_seconds), + _enum_coerce(SunEvent, event), + _enum_coerce(EventDirection, direction), + altitude, + pressure, + temperature, + COpt_Float64.from_optional(delta_t), + atmos_refract, + ) + return _r.value + + +def solar_position_inplace( + latitude, + longitude, + unix_seconds, + azimuth, + elevation, + zenith, + apparent_elevation, + apparent_zenith, + *, + altitude=0.0, + algorithm=Algorithm.PSA, + refraction=RefractionModel.DEFAULT, + pressure=101325.0, + temperature=12.0, + delta_t=67.0, + atmos_refract=0.5667, + refraction_limit=-0.5667, + psa_coeffs=2020, + gmst_option=1, + spencer_correction=True, + julian_date=JulianDateMode.ORIGINAL, +): + """Solar position for many instants, written into caller-allocated output buffers. + + The five output arrays must each hold at least ``len(unix_seconds)`` elements. They + are borrowed, so the writes land directly in the caller's numpy arrays and nothing + is allocated. The observer is built once and only read, so this is also the form the + threaded entrypoint parallelises. + + ``unix_seconds`` may be a sequence of datetimes instead of a ``float64`` array; see + ``solar_position_single``. Converting one allocates, so pass Unix seconds when the + point of this entrypoint is to allocate nothing. + + See ``solar_position_single`` for the keyword arguments. + """ + return _lowlevel.libsolarposition_solar_position_inplace( + latitude, + longitude, + CVector_borrowed_Float64.from_numpy(_unix_seconds_array(unix_seconds)), + CVector_borrowed_Float64.from_numpy(azimuth), + CVector_borrowed_Float64.from_numpy(elevation), + CVector_borrowed_Float64.from_numpy(zenith), + CVector_borrowed_Float64.from_numpy(apparent_elevation), + CVector_borrowed_Float64.from_numpy(apparent_zenith), + altitude, + _enum_coerce(Algorithm, algorithm), + _enum_coerce(RefractionModel, refraction), + pressure, + temperature, + COpt_Float64.from_optional(delta_t), + atmos_refract, + refraction_limit, + psa_coeffs, + gmst_option, + spencer_correction, + _enum_coerce(JulianDateMode, julian_date), + ) + + +def solar_position_inplace_threaded( + latitude, + longitude, + unix_seconds, + azimuth, + elevation, + zenith, + apparent_elevation, + apparent_zenith, + *, + altitude=0.0, + algorithm=Algorithm.PSA, + refraction=RefractionModel.DEFAULT, + pressure=101325.0, + temperature=12.0, + delta_t=67.0, + atmos_refract=0.5667, + refraction_limit=-0.5667, + psa_coeffs=2020, + gmst_option=1, + spencer_correction=True, + julian_date=JulianDateMode.ORIGINAL, +): + """Solar position for many instants, computed in parallel into caller-held buffers. + + Identical to ``solar_position_inplace`` except that the loop runs across the + Julia + threads the library was compiled with. Call ``julia_nthreads()`` to find out how + many + that is; a single-threaded build runs this serially and gives the same answers. + + The five output arrays must each hold at least ``len(unix_seconds)`` elements and + are + written in place. See ``solar_position_single`` for the keyword arguments. + """ + return _lowlevel.libsolarposition_solar_position_inplace_threaded( + latitude, + longitude, + CVector_borrowed_Float64.from_numpy(_unix_seconds_array(unix_seconds)), + CVector_borrowed_Float64.from_numpy(azimuth), + CVector_borrowed_Float64.from_numpy(elevation), + CVector_borrowed_Float64.from_numpy(zenith), + CVector_borrowed_Float64.from_numpy(apparent_elevation), + CVector_borrowed_Float64.from_numpy(apparent_zenith), + altitude, + _enum_coerce(Algorithm, algorithm), + _enum_coerce(RefractionModel, refraction), + pressure, + temperature, + COpt_Float64.from_optional(delta_t), + atmos_refract, + refraction_limit, + psa_coeffs, + gmst_option, + spencer_correction, + _enum_coerce(JulianDateMode, julian_date), + ) + + +def solar_position( + latitude, + longitude, + unix_seconds, + *, + altitude=0.0, + algorithm=Algorithm.PSA, + refraction=RefractionModel.DEFAULT, + pressure=101325.0, + temperature=12.0, + delta_t=67.0, + atmos_refract=0.5667, + refraction_limit=-0.5667, + psa_coeffs=2020, + gmst_option=1, + spencer_correction=True, + julian_date=JulianDateMode.ORIGINAL, +): + """Solar position for many instants, returned as an ``n x 5`` array. + + The columns are azimuth, elevation, zenith, apparent elevation and apparent zenith. + This allocates the result for you; use ``solar_position_inplace`` when you + already hold the output buffers. See ``solar_position_single`` for the keyword + arguments. + + ``unix_seconds`` may be a ``float64`` array of Unix seconds, a ``datetime64`` array, + a pandas ``DatetimeIndex`` or ``Series``, or a list of ``datetime`` objects:: + + >>> import pandas as pd, solarposition as sp + >>> t = pd.date_range("2023-06-21", periods=24, freq="h", tz="UTC") + >>> sp.solar_position(40.0, -105.0, t).shape # doctest: +SKIP + (24, 5) + + Why this overrides the generated binding: the C entrypoint that returns an owning + matrix is called correctly from C, but some builds of libffi marshal that particular + signature wrongly and silently pass ``delta_t`` where ``latitude`` belongs. The + trigger is a struct argument split across classes (``COpt_Float64`` is an integer + flag plus a double) whose integer half lands in the last available general-purpose + register; libffi then puts its double half in ``xmm0``, on top of the first + floating-point argument. Only this entrypoint has that shape — every other one + either passes ``delta_t`` on the stack or leaves a register spare. Observed with the + libffi bundled in python-build-standalone 3.13.7 (what ``uv`` installs by default) + and not with Fedora's libffi 3.5.2, so it cannot be avoided by a version check. + + Routing through ``solar_position_inplace``, whose signature is unaffected, makes + the answer independent of the host libffi. It is also the faster path, since the + result is written straight into the returned array. + """ + times = _unix_seconds_array(unix_seconds) + n = len(times) + # Column-major, so each column below is a contiguous 1-D view that the C side can + # write into directly. The layout matches what the generated binding returned. + out = np.empty((n, 5), dtype=np.float64, order="F") + _lowlevel.libsolarposition_solar_position_inplace( + latitude, + longitude, + CVector_borrowed_Float64.from_numpy(times), + CVector_borrowed_Float64.from_numpy(out[:, 0]), + CVector_borrowed_Float64.from_numpy(out[:, 1]), + CVector_borrowed_Float64.from_numpy(out[:, 2]), + CVector_borrowed_Float64.from_numpy(out[:, 3]), + CVector_borrowed_Float64.from_numpy(out[:, 4]), + altitude, + _enum_coerce(Algorithm, algorithm), + _enum_coerce(RefractionModel, refraction), + pressure, + temperature, + COpt_Float64.from_optional(delta_t), + atmos_refract, + refraction_limit, + psa_coeffs, + gmst_option, + spencer_correction, + _enum_coerce(JulianDateMode, julian_date), + ) + return out From a5d45fe0a249cfa6081b6fe56ed449bbb5f6989a Mon Sep 17 00:00:00 2001 From: Stefan de Lange Date: Sun, 13 Sep 2026 02:00:17 +0200 Subject: [PATCH 06/13] Add the pytest suite and fixture generator Co-Authored-By: Claude Opus 5 (1M context) --- lib/test/generate_fixtures.jl | 210 ++++++++++ lib/test/test_solarposition.py | 590 ++++++++++++++++++++++++++++ lib/test/test_threads_subprocess.py | 72 ++++ 3 files changed, 872 insertions(+) create mode 100644 lib/test/generate_fixtures.jl create mode 100644 lib/test/test_solarposition.py create mode 100644 lib/test/test_threads_subprocess.py diff --git a/lib/test/generate_fixtures.jl b/lib/test/generate_fixtures.jl new file mode 100644 index 00000000..b899a319 --- /dev/null +++ b/lib/test/generate_fixtures.jl @@ -0,0 +1,210 @@ +#!/usr/bin/env julia +# +# Emit lib/test/fixtures.json, the reference data the pytest suite compares against. +# +# julia --project=test lib/test/generate_fixtures.jl +# +# Two independent sets of numbers go in: +# +# "reference" — the solposx tables from test/positioning/expected-values.jl, on the +# 84.375 s grid so the full Julian Date is exactly representable in Float64. These are +# the authoritative cross-implementation values and are what catches a real numeric +# regression in the compiled library. +# +# "package" — what this checkout of SolarPosition.jl itself returns for every +# algorithm x refraction pair at the same instants. Comparing against these catches +# anything that `juliac --trim` changed about the numerics, at full Float64 precision. + +using Dates: DateTime, datetime2unix +using TimeZones: ZonedDateTime, UTC +using SolarPosition + +const REPO = dirname(dirname(@__DIR__)) +include(joinpath(REPO, "test", "positioning", "expected-values.jl")) + +# Minimal JSON writer. The fixture is only numbers, strings and nested containers, and +# `string(::Float64)` already emits the shortest round-tripping form, so this avoids +# adding a JSON dependency to the test project. +jsonesc(s) = replace(string(s), '\\' => "\\\\", '"' => "\\\"") +tojson(x::AbstractString) = string('"', jsonesc(x), '"') +tojson(x::Bool) = x ? "true" : "false" +tojson(::Nothing) = "null" +tojson(::Missing) = "null" +tojson(x::Integer) = string(x) +function tojson(x::AbstractFloat) + isfinite(x) || error("non-finite value in fixture: $x") + return string(Float64(x)) +end +tojson(v::AbstractVector) = string('[', join(map(tojson, v), ','), ']') +tojson(d::AbstractDict) = + string('{', join(("$(tojson(string(k))):$(tojson(v))" for (k, v) in d), ','), '}') + +const ALGORITHMS = [ + "PSA" => PSA(), + "NOAA" => NOAA(), + "SPA" => SPA(), + "WALRAVEN" => Walraven(), + "USNO" => USNO(), + "IQBAL" => Iqbal(), + "MICHALSKY" => Michalsky(), +] + +# Pressure, temperature, and SPARefraction's refraction limit must match the defaults the +# Python entrypoints use, or the comparison is meaningless. +const PRESSURE, TEMPERATURE, REFRACTION_LIMIT = 101325.0, 12.0, -0.5667 + +const REFRACTIONS = [ + "DEFAULT" => DefaultRefraction(), + "NONE" => NoRefraction(), + "HUGHES" => HUGHES(PRESSURE, TEMPERATURE), + "ARCHER" => ARCHER(), + "BENNETT" => BENNETT(PRESSURE, TEMPERATURE), + "MICHALSKY" => MICHALSKY(), + "SG2" => SG2(PRESSURE, TEMPERATURE), + "SPA" => SPARefraction(PRESSURE, TEMPERATURE, REFRACTION_LIMIT), +] + +apparent(p::SolPos) = (p.elevation, p.zenith) +apparent(p::ApparentSolPos) = (p.apparent_elevation, p.apparent_zenith) + +function main() + conds = test_conditions() + + times, lats, lons, alts = Float64[], Float64[], Float64[], Float64[] + datetimes = DateTime[] + for (t, lat, lon, alt) in zip(conds.time, conds.latitude, conds.longitude, conds.altitude) + dt = DateTime(t, UTC) # `test_conditions` already parsed these to ZonedDateTime + push!(datetimes, dt) + push!(times, datetime2unix(dt)) + push!(lats, Float64(lat)) + push!(lons, Float64(lon)) + push!(alts, ismissing(alt) ? 0.0 : Float64(alt)) + end + + # The solposx reference tables. Each carries the call that reproduces it and its own + # column names, because the tables are not uniform: NOAA is validated against + # HUGHES at temperature 10.0 rather than the 12.0 default, USNO has two GMST + # variants, and SPA compares at 1e-8 with an extra column this API does not expose. + # The specs mirror test/positioning/test-*.jl exactly. + reference_specs = [ + ( + name = "PSA_2020", table = expected_2020, algorithm = "PSA", + refraction = "DEFAULT", psa_coeffs = 2020, gmst_option = 1, + temperature = TEMPERATURE, atol = 1.0e-10, skip_poles = false, + ), + ( + name = "PSA_2001", table = expected_2001, algorithm = "PSA", + refraction = "DEFAULT", psa_coeffs = 2001, gmst_option = 1, + temperature = TEMPERATURE, atol = 1.0e-10, skip_poles = false, + ), + ( + name = "NOAA", table = expected_noaa, algorithm = "NOAA", + refraction = "HUGHES", psa_coeffs = 2020, gmst_option = 1, + temperature = 10.0, atol = 1.0e-10, skip_poles = true, + ), + ( + name = "WALRAVEN", table = expected_walraven, algorithm = "WALRAVEN", + refraction = "DEFAULT", psa_coeffs = 2020, gmst_option = 1, + temperature = TEMPERATURE, atol = 1.0e-10, skip_poles = false, + ), + ( + name = "USNO", table = expected_usno, algorithm = "USNO", + refraction = "DEFAULT", psa_coeffs = 2020, gmst_option = 1, + temperature = TEMPERATURE, atol = 1.0e-10, skip_poles = false, + ), + ( + name = "USNO_OPTION_2", table = expected_usno_option_2, algorithm = "USNO", + refraction = "DEFAULT", psa_coeffs = 2020, gmst_option = 2, + temperature = TEMPERATURE, atol = 1.0e-10, skip_poles = false, + ), + ( + name = "SPA", table = expected_spa, algorithm = "SPA", + refraction = "DEFAULT", psa_coeffs = 2020, gmst_option = 1, + temperature = TEMPERATURE, atol = 1.0e-8, skip_poles = false, + ), + ( + name = "IQBAL", table = expected_iqbal, algorithm = "IQBAL", + refraction = "DEFAULT", psa_coeffs = 2020, gmst_option = 1, + temperature = TEMPERATURE, atol = 1.0e-10, skip_poles = false, + ), + ] + + reference = Dict{String, Any}() + for s in reference_specs + df = s.table() + reference[s.name] = Dict{String, Any}( + "algorithm" => s.algorithm, + "refraction" => s.refraction, + "psa_coeffs" => s.psa_coeffs, + "gmst_option" => s.gmst_option, + "temperature" => s.temperature, + "atol" => s.atol, + "skip_poles" => s.skip_poles, + "columns" => String.(names(df)), + "rows" => [[Float64(v) for v in r] for r in eachrow(df)], + ) + end + + # This checkout's own output for every algorithm x refraction pair, as + # (azimuth, elevation, zenith, apparent_elevation, apparent_zenith). + package = Dict{String, Any}() + for (aname, alg) in ALGORITHMS, (rname, refr) in REFRACTIONS + rows = Vector{Vector{Float64}}(undef, length(datetimes)) + for i in eachindex(datetimes) + obs = Observer(lats[i], lons[i], alts[i]) + p = solar_position(obs, datetimes[i], alg, refr) + (ae, az) = apparent(p) + rows[i] = [p.azimuth, p.elevation, p.zenith, ae, az] + end + package["$(aname)|$(rname)"] = rows + end + + # Sunrise, sunset and transit. SPA only — `Utilities._transit_sunrise_sunset` has no + # method for any other algorithm. The options match the Python entrypoint defaults. + srt_alg = SPA(67.0, PRESSURE, TEMPERATURE, 0.5667) + sun_events = Vector{Vector{Float64}}(undef, length(datetimes)) + for i in eachindex(datetimes) + obs = Observer(lats[i], lons[i], alts[i]) + r = transit_sunrise_sunset_seconds(obs, datetimes[i], srt_alg) + sun_events[i] = [r.transit, r.sunrise, r.sunset] + end + + # next_/previous_ sunrise, sunset and solar noon, as Unix seconds. + event_fns = [ + "SUNRISE|NEXT" => next_sunrise, "SUNSET|NEXT" => next_sunset, + "SOLAR_NOON|NEXT" => next_solar_noon, + "SUNRISE|PREVIOUS" => previous_sunrise, "SUNSET|PREVIOUS" => previous_sunset, + "SOLAR_NOON|PREVIOUS" => previous_solar_noon, + ] + sun_event_lookup = Dict{String, Any}() + for (key, f) in event_fns + sun_event_lookup[key] = [ + datetime2unix(f(Observer(lats[i], lons[i], alts[i]), datetimes[i], srt_alg)) + for i in eachindex(datetimes) + ] + end + + fixture = Dict{String, Any}( + "unix_seconds" => times, + "sun_events" => sun_events, + "sun_event" => sun_event_lookup, + "latitude" => lats, + "longitude" => lons, + "altitude" => alts, + "pressure" => PRESSURE, + "temperature" => TEMPERATURE, + "refraction_limit" => REFRACTION_LIMIT, + "reference" => reference, + "package" => package, + ) + + out = joinpath(@__DIR__, "fixtures.json") + write(out, tojson(fixture)) + println("wrote $out") + println(" conditions: $(length(times))") + println(" reference tables: $(length(reference))") + println(" algorithm x refraction pairs: $(length(package))") + return nothing +end + +main() diff --git a/lib/test/test_solarposition.py b/lib/test/test_solarposition.py new file mode 100644 index 00000000..2122e2ae --- /dev/null +++ b/lib/test/test_solarposition.py @@ -0,0 +1,590 @@ +"""Tests for the compiled solarposition Python bindings. + +Run against an installed build: + + julia --project=test lib/test/generate_fixtures.jl # refresh fixtures.json + julia --project=lib lib/build.jl # build the library + pip install ./lib/out + pytest lib/test -v + +The two substantive tests are `test_matches_solposx_reference`, which checks the +compiled library against the solposx tables that also back the Julia test suite, and +`test_matches_julia_package`, which checks it against this checkout's own Julia output +for every algorithm x refraction pair. The first catches a numeric regression anywhere +in the chain; the second isolates whether `juliac --trim` changed anything. +""" + +import datetime as dt +import json +import pathlib + +import numpy as np +import pytest + +import solarposition as sp + +FIXTURES = pathlib.Path(__file__).parent / "fixtures.json" + +# Column name in a reference table -> attribute on the returned SolarAngles. +_COLUMN_TO_FIELD = { + "elevation": "elevation", + "zenith": "zenith", + "azimuth": "azimuth", + "apparent_elevation": "apparent_elevation", + "apparent_zenith": "apparent_zenith", +} + + +@pytest.fixture(scope="session") +def fx(): + if not FIXTURES.exists(): + pytest.skip( + f"{FIXTURES} missing; run " + "`julia --project=test lib/test/generate_fixtures.jl`" + ) + return json.loads(FIXTURES.read_text()) + + +def _observers(fx): + return list( + zip(fx["latitude"], fx["longitude"], fx["altitude"], fx["unix_seconds"]) + ) + + +# --------------------------------------------------------------------------- basics + + +def test_enums_cover_every_algorithm_and_refraction_model(): + assert set(sp.Algorithm.__members__) == { + "PSA", "NOAA", "SPA", "WALRAVEN", "USNO", "IQBAL", "MICHALSKY", + } + assert set(sp.RefractionModel.__members__) == { + "DEFAULT", "NONE", "HUGHES", "ARCHER", "BENNETT", "MICHALSKY", "SG2", "SPA", + } + assert set(sp.JulianDateMode.__members__) == {"ORIGINAL", "STANDARD"} + + +def test_single_returns_five_angles(): + r = sp.solar_position_single(52.35888, 4.88185, 1687348800.0, altitude=100.0) + assert r.zenith == pytest.approx(90.0 - r.elevation, abs=1e-9) + assert -90.0 <= r.elevation <= 90.0 + assert 0.0 <= r.azimuth < 360.0 + + +def test_refraction_raises_apparent_elevation(): + # Refraction bends light downward, so the sun appears higher than it is. + args = (40.0, -105.0, 1687348800.0) + for model in ( + sp.RefractionModel.HUGHES, + sp.RefractionModel.ARCHER, + sp.RefractionModel.BENNETT, + sp.RefractionModel.MICHALSKY, + sp.RefractionModel.SG2, + sp.RefractionModel.SPA, + ): + r = sp.solar_position_single(*args, refraction=model) + assert r.apparent_elevation > r.elevation, model + assert r.apparent_zenith == pytest.approx(90.0 - r.apparent_elevation, abs=1e-9) + + +def test_no_refraction_leaves_apparent_equal_to_true(): + r = sp.solar_position_single( + 40.0, -105.0, 1687348800.0, refraction=sp.RefractionModel.NONE + ) + assert r.apparent_elevation == r.elevation + assert r.apparent_zenith == r.zenith + + +def test_delta_t_none_is_accepted(): + fixed = sp.solar_position_single(40.0, -105.0, 1687348800.0, delta_t=67.0) + auto = sp.solar_position_single(40.0, -105.0, 1687348800.0, delta_t=None) + assert auto.elevation == pytest.approx(fixed.elevation, abs=1e-3) + + +# ------------------------------------------------------------------- correctness + + +def test_matches_solposx_reference(fx): + """The compiled library reproduces the solposx tables the Julia tests use.""" + checked = 0 + for name, spec in fx["reference"].items(): + columns = spec["columns"] + atol = spec["atol"] + for (lat, lon, alt, t), row in zip(_observers(fx), spec["rows"]): + # test-noaa.jl skips the poles, where NOAA is numerically unstable. + if spec["skip_poles"] and abs(abs(lat) - 90.0) < 1e-9: + continue + got = sp.solar_position_single( + lat, + lon, + t, + altitude=alt, + algorithm=sp.Algorithm[spec["algorithm"]], + refraction=sp.RefractionModel[spec["refraction"]], + pressure=fx["pressure"], + temperature=spec["temperature"], + psa_coeffs=spec["psa_coeffs"], + gmst_option=spec["gmst_option"], + refraction_limit=fx["refraction_limit"], + ) + for col, expected in zip(columns, row): + field = _COLUMN_TO_FIELD.get(col) + if field is None: + continue # e.g. SPA's equation_of_time, not exposed by this API + assert getattr(got, field) == pytest.approx(expected, abs=atol), ( + f"{name} {col} at lat={lat} lon={lon} t={t}" + ) + checked += 1 + assert checked > 0 + + +def test_matches_julia_package(fx): + """Every algorithm x refraction pair matches the Julia package exactly.""" + worst = 0.0 + for key, rows in fx["package"].items(): + alg_name, refr_name = key.split("|") + for (lat, lon, alt, t), row in zip(_observers(fx), rows): + got = sp.solar_position_single( + lat, + lon, + t, + altitude=alt, + algorithm=sp.Algorithm[alg_name], + refraction=sp.RefractionModel[refr_name], + pressure=fx["pressure"], + temperature=fx["temperature"], + refraction_limit=fx["refraction_limit"], + ) + mine = [ + got.azimuth, + got.elevation, + got.zenith, + got.apparent_elevation, + got.apparent_zenith, + ] + for a, b in zip(mine, row): + worst = max(worst, abs(a - b)) + assert mine == pytest.approx(row, abs=1e-12, rel=0), key + # The compiled library runs the same Float64 code, so agreement is near-exact. + assert worst < 1e-12, f"worst deviation from the Julia package: {worst}" + + +# ------------------------------------------------------------------------ batch + + +def test_many_instants_match_single(fx): + lat, lon, alt = 52.35888, 4.88185, 100.0 + times = np.array(fx["unix_seconds"], dtype=np.float64) + for alg in sp.Algorithm: + m = sp.solar_position(lat, lon, times, altitude=alt, algorithm=alg) + assert m.shape == (len(times), 5) + for i, t in enumerate(times): + r = sp.solar_position_single( + lat, lon, float(t), altitude=alt, algorithm=alg + ) + assert list(m[i]) == [ + r.azimuth, + r.elevation, + r.zenith, + r.apparent_elevation, + r.apparent_zenith, + ], f"{alg} row {i}" + + +def test_inplace_matches_allocating(fx): + lat, lon, alt = 52.35888, 4.88185, 100.0 + times = np.array(fx["unix_seconds"], dtype=np.float64) + n = len(times) + bufs = [np.zeros(n) for _ in range(5)] + sp.solar_position_inplace(lat, lon, times, *bufs, altitude=alt) + m = sp.solar_position(lat, lon, times, altitude=alt) + for col, buf in enumerate(bufs): + np.testing.assert_array_equal(buf, m[:, col]) + + +@pytest.mark.parametrize( + "algorithm", [sp.Algorithm.PSA, sp.Algorithm.WALRAVEN, sp.Algorithm.IQBAL] +) +def test_delta_t_does_not_leak_into_other_arguments(algorithm): + """`delta_t` must reach the argument it names, on every entrypoint. + + These three algorithms ignore `delta_t` entirely, so changing it must change + nothing. That makes this a cheap detector for a whole class of ABI-marshalling + faults, where a value arrives in the wrong slot: `delta_t` is the only argument + passed as a struct split across an integer and a floating-point half, which is + exactly the shape ctypes/libffi is most likely to place wrongly. One such libffi + bug silently put `delta_t` where `latitude` belongs on `solar_position`; + see `solar_position` in `lib/python/_extras.py`. + """ + lat, lon, alt = 52.35888, 4.88185, 100.0 + times = np.linspace(1687348800.0, 1687435200.0, 64) + kwargs = dict(altitude=alt, algorithm=algorithm) + + def sample(delta_t): + one = sp.solar_position_single(lat, lon, float(times[0]), + delta_t=delta_t, **kwargs) + many = sp.solar_position(lat, lon, times, delta_t=delta_t, **kwargs) + inplace = [np.empty(len(times)) for _ in range(5)] + sp.solar_position_inplace(lat, lon, times, *inplace, + delta_t=delta_t, **kwargs) + threaded = [np.empty(len(times)) for _ in range(5)] + sp.solar_position_inplace_threaded(lat, lon, times, *threaded, + delta_t=delta_t, **kwargs) + return one.azimuth, many, inplace, threaded + + one_a, many_a, inplace_a, threaded_a = sample(67.0) + one_b, many_b, inplace_b, threaded_b = sample(0.0) + + assert one_a == one_b, "solar_position_single" + np.testing.assert_array_equal(many_a, many_b, err_msg="solar_position") + for col, (a, b) in enumerate(zip(inplace_a, inplace_b)): + np.testing.assert_array_equal(a, b, err_msg=f"inplace column {col}") + for col, (a, b) in enumerate(zip(threaded_a, threaded_b)): + np.testing.assert_array_equal(a, b, err_msg=f"threaded column {col}") + + # And all four must agree with each other, not merely be self-consistent. + np.testing.assert_array_equal(many_a[:, 0], inplace_a[0]) + np.testing.assert_array_equal(many_a[:, 0], threaded_a[0]) + assert many_a[0, 0] == one_a + + +def test_solar_position_uses_the_hand_maintained_wrapper(): + """`solar_position` must come from `_extras`, not from the generated façade. + + The generated binding for this one entrypoint is mismarshalled by some builds of + libffi, so the override is the fix rather than a convenience. If a build stops + installing it, every result from this function silently becomes wrong, and this + is the tripwire for that. + """ + assert sp.solar_position.__module__.endswith("_extras") + + +def test_inplace_rejects_short_buffers(): + times = np.zeros(10) + short = [np.zeros(3) for _ in range(5)] + with pytest.raises(sp.JLWError) as excinfo: + sp.solar_position_inplace(0.0, 0.0, times, *short) + # DimensionMismatch maps to status code 3. + assert excinfo.value.code == 3 + + +def test_inplace_writes_through_to_caller_arrays(): + times = np.array([1687348800.0, 1687352400.0]) + bufs = [np.full(2, np.nan) for _ in range(5)] + sp.solar_position_inplace(40.0, -105.0, times, *bufs) + for buf in bufs: + assert np.all(np.isfinite(buf)) + + +def test_empty_input_is_allowed(): + m = sp.solar_position(40.0, -105.0, np.zeros(0)) + assert m.shape == (0, 5) + + +# ---------------------------------------------------------------- datetime input + +UTC_INSTANT = dt.datetime(2023, 6, 21, 3, 45, tzinfo=dt.timezone.utc) +UTC_SECONDS = UTC_INSTANT.timestamp() + + +def test_single_accepts_a_datetime(): + a = sp.solar_position_single(40.0, -105.0, UTC_INSTANT, algorithm=sp.Algorithm.SPA) + b = sp.solar_position_single(40.0, -105.0, UTC_SECONDS, algorithm=sp.Algorithm.SPA) + assert a.elevation == b.elevation + assert a.azimuth == b.azimuth + + +@pytest.mark.parametrize( + "when", + [ + # A naive datetime is read as UTC, not as the machine's local time. That is a + # deliberate departure from `datetime.timestamp()`: everything this library + # takes and returns is UTC, and numpy's datetime64 carries no zone to read. + dt.datetime(2023, 6, 21, 3, 45), + dt.datetime(2023, 6, 21, 3, 45, tzinfo=dt.timezone.utc), + # Same instant written in another zone; the offset must be applied. + dt.datetime(2023, 6, 21, 5, 45, tzinfo=dt.timezone(dt.timedelta(hours=2))), + np.datetime64("2023-06-21T03:45"), + ], +) +def test_every_spelling_of_one_instant_agrees(when): + got = sp.solar_position_single(40.0, -105.0, when) + ref = sp.solar_position_single(40.0, -105.0, UTC_SECONDS) + assert got.elevation == ref.elevation + + +def test_a_bare_date_means_midnight_utc(): + midnight = dt.datetime(2023, 6, 21, tzinfo=dt.timezone.utc) + assert ( + sp.solar_position_single(40.0, -105.0, dt.date(2023, 6, 21)).elevation + == sp.solar_position_single(40.0, -105.0, midnight).elevation + ) + + +def _hourly(n=5): + """`n` hourly instants, as Unix seconds and as datetimes.""" + seconds = np.array([UTC_SECONDS + 3600.0 * i for i in range(n)]) + return seconds, [dt.datetime.fromtimestamp(t, dt.timezone.utc) for t in seconds] + + +def test_solar_position_accepts_datetime_sequences(): + seconds, datetimes = _hourly() + ref = sp.solar_position(40.0, -105.0, seconds) + for label, times in [ + ("list of datetime", datetimes), + ("tuple of datetime", tuple(datetimes)), + # datetime64 has no zone, so strip the (UTC) one rather than let numpy + # warn about dropping it. + ( + "datetime64 array", + np.array([d.replace(tzinfo=None) for d in datetimes], "datetime64[s]"), + ), + ("object array", np.array(datetimes, dtype=object)), + ]: + got = sp.solar_position(40.0, -105.0, times) + np.testing.assert_array_equal(got, ref, err_msg=label) + + +def test_solar_position_accepts_pandas_datetime_containers(): + """pandas is the common caller, and its tz-aware containers take a separate path.""" + pd = pytest.importorskip("pandas") + seconds, _ = _hourly() + ref = sp.solar_position(40.0, -105.0, seconds) + aware = pd.to_datetime(seconds, unit="s", utc=True) + for label, times in [ + ("DatetimeIndex utc", aware), + ("DatetimeIndex naive", aware.tz_localize(None)), + ("DatetimeIndex offset", aware.tz_convert("Europe/Amsterdam")), + ("Series utc", pd.Series(aware)), + ("Series naive", pd.Series(aware.tz_localize(None))), + ]: + got = sp.solar_position(40.0, -105.0, times) + np.testing.assert_array_equal(got, ref, err_msg=label) + + +@pytest.mark.parametrize( + "entrypoint", + [sp.solar_position_inplace, sp.solar_position_inplace_threaded], +) +def test_inplace_accepts_datetimes(entrypoint): + seconds, datetimes = _hourly() + ref = sp.solar_position(40.0, -105.0, seconds) + bufs = [np.zeros(len(seconds)) for _ in range(5)] + entrypoint(40.0, -105.0, datetimes, *bufs) + for col, buf in enumerate(bufs): + np.testing.assert_array_equal(buf, ref[:, col]) + + +def test_sunrise_sunset_entrypoints_accept_datetimes(): + assert ( + sp.transit_sunrise_sunset(52.35888, 4.88185, UTC_INSTANT).sunrise + == sp.transit_sunrise_sunset(52.35888, 4.88185, UTC_SECONDS).sunrise + ) + assert sp.sun_event( + 52.35888, 4.88185, UTC_INSTANT, sp.SunEvent.SUNRISE, sp.EventDirection.NEXT + ) == sp.sun_event( + 52.35888, 4.88185, UTC_SECONDS, sp.SunEvent.SUNRISE, sp.EventDirection.NEXT + ) + + +# -------------------------------------------------------------------- validation + + +def test_out_of_range_algorithm_is_rejected(): + """An invalid enum value must raise, not silently fall through to an algorithm.""" + with pytest.raises((ValueError, KeyError, TypeError)): + sp.solar_position_single(40.0, -105.0, 1687348800.0, algorithm=99) + + +def test_out_of_range_refraction_is_rejected(): + with pytest.raises((ValueError, KeyError, TypeError)): + sp.solar_position_single(40.0, -105.0, 1687348800.0, refraction=99) + + +def test_wrong_dtype_array_is_rejected(): + with pytest.raises((ValueError, TypeError)): + sp.solar_position(40.0, -105.0, np.array([1, 2, 3], dtype=np.int64)) + + +def test_non_contiguous_array_is_rejected(): + times = np.linspace(1687348800.0, 1687366800.0, 20)[::2] + assert not times.flags["C_CONTIGUOUS"] + with pytest.raises((ValueError, TypeError)): + sp.solar_position(40.0, -105.0, times) + + +# ------------------------------------------------------------------- consistency + + +def test_psa_coefficient_sets_differ(): + a = sp.solar_position_single(40.0, -105.0, 1687348800.0, psa_coeffs=2020) + b = sp.solar_position_single(40.0, -105.0, 1687348800.0, psa_coeffs=2001) + assert a.elevation != b.elevation + + +def test_usno_gmst_options_are_both_accepted(): + for opt in (1, 2): + r = sp.solar_position_single( + 40.0, -105.0, 1687348800.0, algorithm=sp.Algorithm.USNO, gmst_option=opt + ) + assert np.isfinite(r.elevation) + + +def test_michalsky_options_are_honoured(): + base = dict(algorithm=sp.Algorithm.MICHALSKY) + orig = sp.solar_position_single( + 40.0, -105.0, 1687348800.0, julian_date=sp.JulianDateMode.ORIGINAL, **base + ) + std = sp.solar_position_single( + 40.0, -105.0, 1687348800.0, julian_date=sp.JulianDateMode.STANDARD, **base + ) + assert np.isfinite(orig.elevation) and np.isfinite(std.elevation) + no_spencer = sp.solar_position_single( + 40.0, -105.0, 1687348800.0, spencer_correction=False, **base + ) + assert np.isfinite(no_spencer.azimuth) + + +# ------------------------------------------------------- sunrise, sunset, transit + + +def test_transit_sunrise_sunset_matches_julia(fx): + for (lat, lon, alt, t), row in zip(_observers(fx), fx["sun_events"]): + got = sp.transit_sunrise_sunset(lat, lon, t, altitude=alt) + assert [got.transit, got.sunrise, got.sunset] == pytest.approx( + row, abs=1e-9, rel=0 + ), f"lat={lat} lon={lon} t={t}" + + +def test_sun_event_matches_julia(fx): + events = { + "SUNRISE": sp.SunEvent.SUNRISE, + "SUNSET": sp.SunEvent.SUNSET, + "SOLAR_NOON": sp.SunEvent.SOLAR_NOON, + } + directions = {"NEXT": sp.EventDirection.NEXT, + "PREVIOUS": sp.EventDirection.PREVIOUS} + for key, expected in fx["sun_event"].items(): + ev, direction = key.split("|") + for (lat, lon, alt, t), want in zip(_observers(fx), expected): + got = sp.sun_event( + lat, lon, t, events[ev], directions[direction], altitude=alt + ) + # The package returns a DateTime here, so this is whole-second exact. + assert got == pytest.approx(want, abs=1e-6), f"{key} lat={lat} t={t}" + + +def test_known_sunrise_sunset_for_amsterdam(): + """Cross-check against the values published in the project README. + + README: transit 2023-06-21T13:42:15+02:00, sunrise 05:18:05+02:00, + sunset 22:06:24+02:00 — i.e. 11:42:15, 03:18:05 and 20:06:24 UTC. + """ + e = sp.transit_sunrise_sunset(52.35888, 4.88185, 1687348800.0, altitude=100.0) + assert round(e.transit) == 11 * 3600 + 42 * 60 + 15 + assert round(e.sunrise) == 3 * 3600 + 18 * 60 + 5 + assert round(e.sunset) == 20 * 3600 + 6 * 60 + 24 + + +def test_sunrise_precedes_transit_precedes_sunset(): + e = sp.transit_sunrise_sunset(40.0, -105.0, 1687348800.0) + assert e.sunrise < e.transit < e.sunset + + +def test_next_and_previous_bracket_the_query_time(): + t = 1687348800.0 + nxt = sp.sun_event(40.0, -105.0, t, sp.SunEvent.SUNRISE, sp.EventDirection.NEXT) + prv = sp.sun_event(40.0, -105.0, t, sp.SunEvent.SUNRISE, sp.EventDirection.PREVIOUS) + assert prv < t < nxt + # Consecutive sunrises are one day apart, give or take the equation of time. + assert 23 * 3600 < nxt - prv < 25 * 3600 + + +def test_polar_night_returns_zeros(): + """The package warns and returns zeros where the sun never rises.""" + # 2023-12-21, near the pole + e = sp.transit_sunrise_sunset(89.0, 0.0, 1703116800.0) + assert (e.transit, e.sunrise, e.sunset) == (0.0, 0.0, 0.0) + + +# -------------------------------------------------------------------- threading + + +def test_julia_nthreads_is_reported(): + n = sp.julia_nthreads() + assert isinstance(n, int) and n >= 1 + + +def test_set_num_threads_rejects_bad_input(): + with pytest.raises(ValueError): + sp.set_num_threads(0) + with pytest.raises(ValueError): + sp.set_num_threads(-1) + + +def test_set_num_threads_after_runtime_start_raises(): + """The thread count is fixed once the runtime is up, so a late call must not pass + silently. Every other test in this file has already started the runtime, so by the + time this runs the call is guaranteed to be late.""" + sp.julia_nthreads() # ensure the runtime is running regardless of test order + current = sp.julia_nthreads() + with pytest.raises(RuntimeError, match="already running"): + sp.set_num_threads(current + 1) + + +def test_threaded_matches_serial(fx): + """The parallel loop must agree with the serial one, elementwise, per algorithm. + + This is meaningful on a single-threaded build too: it exercises the separate + all-positional entrypoint and the hand-maintained keyword wrapper over it. + """ + lat, lon, alt = 52.35888, 4.88185, 100.0 + times = np.array(fx["unix_seconds"], dtype=np.float64) + n = len(times) + for alg in sp.Algorithm: + bufs = [np.zeros(n) for _ in range(5)] + sp.solar_position_inplace_threaded( + lat, lon, times, *bufs, altitude=alt, algorithm=alg + ) + serial = sp.solar_position(lat, lon, times, altitude=alt, algorithm=alg) + for col, buf in enumerate(bufs): + np.testing.assert_array_equal( + buf, serial[:, col], err_msg=f"{alg} col {col}" + ) + + +def test_threaded_rejects_short_buffers(): + times = np.zeros(10) + short = [np.zeros(3) for _ in range(5)] + with pytest.raises(sp.JLWError) as excinfo: + sp.solar_position_inplace_threaded(0.0, 0.0, times, *short) + assert excinfo.value.code == 3 + + +def test_threaded_is_repeatable_over_many_calls(): + """A parallel region re-entered many times must stay stable and deterministic.""" + times = np.linspace(1687348800.0, 1687435200.0, 500) + first = None + for _ in range(25): + bufs = [np.zeros(len(times)) for _ in range(5)] + sp.solar_position_inplace_threaded( + 40.0, -105.0, times, *bufs, algorithm=sp.Algorithm.SPA + ) + if first is None: + first = [b.copy() for b in bufs] + else: + for a, b in zip(first, bufs): + np.testing.assert_array_equal(a, b) + + +# ------------------------------------------------------------------- consistency + + +def test_extreme_latitudes_and_epochs_are_finite(fx): + """The fixture conditions include the poles and the years 1800 and 2200.""" + for lat, lon, alt, t in _observers(fx): + for alg in sp.Algorithm: + r = sp.solar_position_single( + lat, lon, t, altitude=alt, algorithm=alg + ) + assert np.isfinite(r.azimuth) and np.isfinite(r.elevation), (alg, lat, t) diff --git a/lib/test/test_threads_subprocess.py b/lib/test/test_threads_subprocess.py new file mode 100644 index 00000000..a1571bc3 --- /dev/null +++ b/lib/test/test_threads_subprocess.py @@ -0,0 +1,72 @@ +"""Thread-count tests that each need a fresh Julia runtime. + +The thread count is fixed once the runtime starts, so these cannot share a process with +each other or with the rest of the suite. Each runs in a subprocess. +""" + +import subprocess +import sys + +import pytest + + +def run(code: str) -> subprocess.CompletedProcess: + return subprocess.run( + [sys.executable, "-c", code], capture_output=True, text=True, timeout=300 + ) + + +@pytest.mark.parametrize("n", [1, 2, 4, 8]) +def test_set_num_threads_takes_effect(n): + """set_num_threads(n) before the first call must actually give n Julia threads.""" + proc = run( + f"import solarposition as sp;" + f"got = sp.set_num_threads({n});" + f"print(got, sp.julia_nthreads())" + ) + assert proc.returncode == 0, proc.stderr + got, reported = proc.stdout.split() + assert int(got) == n + assert int(reported) == n + + +def test_default_build_is_single_threaded(): + """Without an explicit request the library must start single-threaded, which is the + signal-safe configuration.""" + proc = run("import solarposition as sp; print(sp.julia_nthreads())") + assert proc.returncode == 0, proc.stderr + assert int(proc.stdout.strip()) == 1 + + +def test_threaded_results_match_serial_with_runtime_threads(): + """Raising the thread count at runtime must not change any answer.""" + proc = run( + "import numpy as np, solarposition as sp\n" + "sp.set_num_threads(4)\n" + "t = np.linspace(1687348800.0, 1687435200.0, 20000)\n" + "a = [np.empty(len(t)) for _ in range(5)]\n" + "b = [np.empty(len(t)) for _ in range(5)]\n" + "sp.solar_position_inplace(52.35888, 4.88185, t, *a," + " algorithm=sp.Algorithm.SPA)\n" + "sp.solar_position_inplace_threaded(52.35888, 4.88185, t, *b," + " algorithm=sp.Algorithm.SPA)\n" + "print(all(np.array_equal(x, y) for x, y in zip(a, b)))\n" + ) + assert proc.returncode == 0, proc.stderr + assert proc.stdout.strip() == "True" + + +def test_multithreaded_process_exits_cleanly(): + """On Julia >= 1.12 a multi-threaded library with signal handling off segfaults at + shutdown. set_num_threads() turns signal handling on for n > 1 to avoid exactly that, + so a threaded process must still exit 0.""" + for _ in range(5): + proc = run( + "import numpy as np, solarposition as sp\n" + "sp.set_num_threads(8)\n" + "t = np.linspace(1687348800.0, 1687435200.0, 5000)\n" + "b = [np.empty(len(t)) for _ in range(5)]\n" + "sp.solar_position_inplace_threaded(40.0, -105.0, t, *b," + " algorithm=sp.Algorithm.SPA)\n" + ) + assert proc.returncode == 0, f"exit {proc.returncode}: {proc.stderr[-400:]}" From 85e2be373d64192d7f8aaa73043d48720ac6f0d5 Mon Sep 17 00:00:00 2001 From: Stefan de Lange Date: Sun, 13 Sep 2026 02:00:22 +0200 Subject: [PATCH 07/13] Add Python examples and the solposx benchmark Co-Authored-By: Claude Opus 5 (1M context) --- lib/examples/benchmark_vs_solposx.py | 261 +++++++++++++++++++ lib/examples/getting_started.py | 370 +++++++++++++++++++++++++++ lib/examples/simple.py | 27 ++ 3 files changed, 658 insertions(+) create mode 100644 lib/examples/benchmark_vs_solposx.py create mode 100644 lib/examples/getting_started.py create mode 100644 lib/examples/simple.py diff --git a/lib/examples/benchmark_vs_solposx.py b/lib/examples/benchmark_vs_solposx.py new file mode 100644 index 00000000..b1f7c472 --- /dev/null +++ b/lib/examples/benchmark_vs_solposx.py @@ -0,0 +1,261 @@ +#!/usr/bin/env python3 +"""Benchmark the wrapped Julia library against solposx, the Python reference package. + + pip install ./lib/out solposx + python lib/examples/benchmark_vs_solposx.py + python lib/examples/benchmark_vs_solposx.py --sizes 1,1000,100000 --repeats 7 + +solposx is the package this library's reference values were generated from, and it +implements the same seven algorithms, so this is a like-for-like comparison of the same +maths in two languages. + +Method notes, because they decide whether the numbers mean anything: + +* Inputs are prepared in each library's native form *before* timing — a pandas + ``DatetimeIndex`` for solposx, a float64 array of Unix seconds for this library. The + conversion between them is measured separately and reported, so you can see whether it + would swamp the difference. +* solposx is fully vectorised over numpy/pandas, so this is not a Python-loop straw man. +* Refraction is switched off on both sides and only the geometric ``elevation`` and + ``azimuth`` are compared, because each algorithm has a different refraction default. +* Every configuration is checked for agreement before it is timed. A disagreement is + reported and that row is skipped rather than producing a meaningless speed number. +* The Julia runtime initialises on first call, so there is a warm-up call outside the + timed region. Best-of-N is reported, not the mean, to suppress scheduler noise. +""" + +from __future__ import annotations + +import argparse +import time + +import numpy as np +import pandas as pd + +import solarposition as sp +from solposx import solarposition as sx + +LATITUDE, LONGITUDE, ALTITUDE = 52.35888, 4.88185, 100.0 + +# Each entry: our enum, the solposx callable, and the kwargs that make the two agree. +# Options are pinned explicitly on both sides rather than relying on matching defaults. +ALGORITHMS = [ + ( + "PSA", + sp.Algorithm.PSA, + lambda times: sx.psa(times, LATITUDE, LONGITUDE, coefficients=2020), + dict(psa_coeffs=2020), + ), + ( + "NOAA", + sp.Algorithm.NOAA, + lambda times: sx.noaa(times, LATITUDE, LONGITUDE, delta_t=67.0), + dict(delta_t=67.0), + ), + ( + "SPA", + sp.Algorithm.SPA, + lambda times: sx.spa( + times, LATITUDE, LONGITUDE, ALTITUDE, + air_pressure=101325.0, temperature=12.0, delta_t=67.0, + ), + dict(altitude=ALTITUDE, pressure=101325.0, temperature=12.0, delta_t=67.0), + ), + ( + "WALRAVEN", + sp.Algorithm.WALRAVEN, + lambda times: sx.walraven(times, LATITUDE, LONGITUDE), + dict(), + ), + ( + "USNO", + sp.Algorithm.USNO, + lambda times: sx.usno(times, LATITUDE, LONGITUDE, delta_t=67.0, gmst_option=1), + dict(delta_t=67.0, gmst_option=1), + ), + ( + "IQBAL", + sp.Algorithm.IQBAL, + lambda times: sx.iqbal(times, LATITUDE, LONGITUDE), + dict(), + ), + ( + "MICHALSKY", + sp.Algorithm.MICHALSKY, + lambda times: sx.michalsky( + times, LATITUDE, LONGITUDE, + spencer_correction=True, julian_date="original", + ), + dict(spencer_correction=True, julian_date=sp.JulianDateMode.ORIGINAL), + ), +] + + +def best_of(fn, repeats: int) -> float: + """Shortest wall time over `repeats` runs, in seconds.""" + timings = [] + for _ in range(repeats): + start = time.perf_counter() + fn() + timings.append(time.perf_counter() - start) + return min(timings) + + +# Divisor per pandas datetime resolution. `asi8` is the index's raw int64 array in its +# own unit, so the unit must be consulted: pandas >= 2 keeps non-nanosecond units +# and `date_range` here yields microseconds. Dividing by a hardcoded 1e9 is a silent +# 1000x error. Going via `.astype("datetime64[ns]")` on the tz-aware index is correct +# but ~6000x slower (3.9 s vs 0.6 ms for a million instants). +_UNIT_PER_SECOND = {"s": 1.0, "ms": 1e3, "us": 1e6, "ns": 1e9} + + +def to_unix_seconds(index: pd.DatetimeIndex) -> np.ndarray: + """DatetimeIndex -> float64 Unix seconds, at the index's own resolution.""" + return index.asi8 / _UNIT_PER_SECOND[index.unit] + + +def make_inputs(n: int): + """The same instants in each library's native input form.""" + index = pd.date_range("2024-01-01", periods=n, freq="1min", tz="UTC") + return index, to_unix_seconds(index) + + +def agreement(name, our_enum, sx_call, our_kwargs, index, unix): + """Max absolute elevation/azimuth difference between the two implementations.""" + ref = sx_call(index) + ours = sp.solar_position( + LATITUDE, LONGITUDE, unix, + algorithm=our_enum, refraction=sp.RefractionModel.NONE, **our_kwargs, + ) + d_elev = np.max(np.abs(ours[:, 1] - ref["elevation"].to_numpy())) + d_azim = np.max(np.abs(ours[:, 0] - ref["azimuth"].to_numpy())) + return max(d_elev, d_azim) + + +def run(sizes, repeats, tolerance, threads): + if threads > 1: + # Must happen before anything else touches the Julia runtime. + sp.set_num_threads(threads) + nthreads = sp.julia_nthreads() + print("Benchmark: wrapped SolarPosition.jl vs solposx") + print(f" julia_nthreads() {nthreads}" + f"{' (single-threaded build)' if nthreads == 1 else ''}") + print(f" numpy {np.__version__}, pandas {pd.__version__}") + print(f" best of {repeats} runs; site {LATITUDE}N {LONGITUDE}E {ALTITUDE:.0f} m") + + # --- agreement, checked once on a decent sample -------------------------------- + print("\nAgreement (max |diff| in elevation/azimuth, degrees, 10k instants)") + index, unix = make_inputs(10_000) + skip = set() + for name, our_enum, sx_call, our_kwargs in ALGORITHMS: + try: + diff = agreement(name, our_enum, sx_call, our_kwargs, index, unix) + except Exception as err: # pragma: no cover - defensive + print(f" {name:<10} could not compare: {type(err).__name__}: {err}") + skip.add(name) + continue + verdict = "ok" if diff <= tolerance else "DISAGREES" + print(f" {name:<10} {diff:.3e} {verdict}") + if diff > tolerance: + skip.add(name) + if skip: + print(f" -> skipping timings for: {', '.join(sorted(skip))}") + + # --- conversion cost, so it can be judged separately --------------------------- + print("\nInput conversion (excluded from the timings below)") + for n in sizes: + index, _ = make_inputs(n) + conv = best_of(lambda: to_unix_seconds(index), repeats) + print(f" n={n:<9,} DatetimeIndex -> Unix seconds: {conv * 1e3:9.4f} ms") + + # --- timings -------------------------------------------------------------------- + for n in sizes: + index, unix = make_inputs(n) + bufs = [np.empty(n) for _ in range(5)] + + print(f"\nn = {n:,} instants") + print(f" {'algorithm':<10} {'solposx':>12} {'julia':>12} " + f"{'julia(inplace)':>15} {'speedup':>9} {'inplace':>9}") + for name, our_enum, sx_call, our_kwargs in ALGORITHMS: + if name in skip: + continue + + def do_sx(): + sx_call(index) + + def do_ours(): + sp.solar_position( + LATITUDE, LONGITUDE, unix, algorithm=our_enum, + refraction=sp.RefractionModel.NONE, **our_kwargs, + ) + + def do_ours_inplace(): + sp.solar_position_inplace( + LATITUDE, LONGITUDE, unix, *bufs, algorithm=our_enum, + refraction=sp.RefractionModel.NONE, **our_kwargs, + ) + + do_sx(); do_ours(); do_ours_inplace() # warm up both sides + + t_sx = best_of(do_sx, repeats) + t_ours = best_of(do_ours, repeats) + t_inplace = best_of(do_ours_inplace, repeats) + + print(f" {name:<10} {t_sx * 1e3:10.4f}ms {t_ours * 1e3:10.4f}ms " + f"{t_inplace * 1e3:13.4f}ms {t_sx / t_ours:8.1f}x " + f"{t_sx / t_inplace:8.1f}x") + + # --- threaded entrypoint -------------------------------------------------------- + if nthreads > 1: + n = max(sizes) + _, unix = make_inputs(n) + bufs = [np.empty(n) for _ in range(5)] + index, _ = make_inputs(n) + print(f"\nThreaded entrypoint, n = {n:,}, {nthreads} Julia threads") + print(f" {'algorithm':<10} {'serial':>12} {'threaded':>12} {'scaling':>8} " + f"{'vs solposx':>11}") + for name, our_enum, sx_call, our_kwargs in ALGORITHMS: + if name in skip: + continue + serial = best_of(lambda: sp.solar_position_inplace( + LATITUDE, LONGITUDE, unix, *bufs, algorithm=our_enum, + refraction=sp.RefractionModel.NONE, **our_kwargs), repeats) + par = best_of(lambda: sp.solar_position_inplace_threaded( + LATITUDE, LONGITUDE, unix, *bufs, algorithm=our_enum, + refraction=sp.RefractionModel.NONE, **our_kwargs), repeats) + t_sx = best_of(lambda: sx_call(index), repeats) + print(f" {name:<10} {serial * 1e3:10.2f}ms {par * 1e3:10.2f}ms " + f"{serial / par:6.2f}x {t_sx / par:9.1f}x") + else: + print("\nThreaded entrypoint skipped: pass --threads=N to measure it.") + + # --- per-call latency ----------------------------------------------------------- + print("\nSingle-instant latency (one position, no arrays)") + one_index, one_unix = make_inputs(1) + for name, our_enum, sx_call, our_kwargs in ALGORITHMS: + if name in skip: + continue + t_sx = best_of(lambda: sx_call(one_index), repeats) + t_ours = best_of(lambda: sp.solar_position_single( + LATITUDE, LONGITUDE, float(one_unix[0]), algorithm=our_enum, + refraction=sp.RefractionModel.NONE, **our_kwargs), repeats) + print(f" {name:<10} solposx {t_sx * 1e6:9.1f}us " + f"julia {t_ours * 1e6:8.2f}us {t_sx / t_ours:7.0f}x") + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--sizes", default="1000,100000,1000000", + help="comma-separated instant counts") + parser.add_argument("--repeats", type=int, default=5) + parser.add_argument("--tolerance", type=float, default=1e-6, + help="max allowed deg difference before a row is skipped") + parser.add_argument("--threads", type=int, default=1, + help="Julia threads to request before the first call") + args = parser.parse_args() + sizes = [int(s) for s in args.sizes.split(",")] + run(sizes, args.repeats, args.tolerance, args.threads) + + +if __name__ == "__main__": + main() diff --git a/lib/examples/getting_started.py b/lib/examples/getting_started.py new file mode 100644 index 00000000..55edc583 --- /dev/null +++ b/lib/examples/getting_started.py @@ -0,0 +1,370 @@ +#!/usr/bin/env python3 +"""Worked examples for the `solarposition` Python bindings. + + julia --project=lib lib/build.jl + pip install ./lib/out + python lib/examples/getting_started.py + +Every entrypoint takes a datetime directly, and Unix seconds in UTC underneath; the +first section shows both, and the conversion back that results still need. Angles are +all in degrees, azimuth measured from north and increasing clockwise. +""" + +from __future__ import annotations + +import datetime as dt +import os +import time + +import numpy as np + +import solarposition as sp + +try: + from zoneinfo import ZoneInfo # Python 3.9+ +except ImportError: # pragma: no cover + ZoneInfo = None + +# Van Gogh Museum, Amsterdam. +LATITUDE, LONGITUDE, ALTITUDE = 52.35888, 4.88185, 100.0 + + +def rule(title: str) -> None: + print(f"\n{title}\n{'-' * len(title)}") + + +# -------------------------------------------------------------------------------------- +# 1. Times +# -------------------------------------------------------------------------------------- +# Wherever a time is taken you may pass a datetime, a date, a numpy datetime64, a +# pandas DatetimeIndex or Series, or plain Unix seconds. A datetime carrying a tzinfo +# has its offset applied; one without is read as UTC, which is what the whole library +# speaks. Results come back as Unix seconds, so only the way back needs a helper. +# +# Unix seconds are still worth holding directly when you build a grid arithmetically, +# as sections 2 and 5 do, or when you are calling in a loop and want no conversion in +# it at all. + + +def to_unix(when: dt.datetime) -> float: + if when.tzinfo is None: + raise ValueError("pass a timezone-aware datetime") + return when.timestamp() + + +def from_unix(seconds: float, tz: dt.tzinfo = dt.timezone.utc) -> dt.datetime: + return dt.datetime.fromtimestamp(seconds, tz) + + +# Nothing here calls the library, deliberately: section 2 has to run before the +# first computation. + + +def times_section() -> dt.datetime: + rule("1. Times") + noon_utc = dt.datetime(2023, 6, 21, 12, 0, tzinfo=dt.timezone.utc) + print(f" {noon_utc.isoformat()} -> {to_unix(noon_utc):.1f} Unix seconds") + print(f" a naive datetime is read as UTC, so {dt.datetime(2023, 6, 21, 12, 0)} " + "means the same instant") + + if ZoneInfo is not None: + local = dt.datetime(2023, 6, 21, 14, 0, tzinfo=ZoneInfo("Europe/Amsterdam")) + print(f" {local.isoformat()} -> {to_unix(local):.1f} (same instant)") + return noon_utc + + +# -------------------------------------------------------------------------------------- +# 2. Threads +# -------------------------------------------------------------------------------------- +# The thread count is chosen at run time, but only before the first computation. That +# is why this comes first. No rebuild is involved. Asking for more than one thread +# switches Julia's signal handlers on for the whole process; see lib/README.md. +# +# Try it: SP_THREADS=8 python lib/examples/getting_started.py + + +def threaded_section() -> None: + rule("2. Threads") + requested = int(os.environ.get("SP_THREADS", "4")) + + # Set it before anything else touches the runtime. Pick the number however you + # like -- an env var here, but os.cpu_count(), argv or a config file all work. + sp.set_num_threads(requested) + print(f" set_num_threads({requested}) -> julia_nthreads() = {sp.julia_nthreads()}") + + start = to_unix(dt.datetime(2024, 1, 1, tzinfo=dt.timezone.utc)) + times = start + np.arange(0, 500_000, dtype=np.float64) * 60.0 + bufs = [np.empty(len(times)) for _ in range(5)] + + def elapsed(fn) -> float: + began = time.perf_counter() + fn() + return time.perf_counter() - began + + serial = elapsed(lambda: sp.solar_position_inplace( + LATITUDE, LONGITUDE, times, *bufs, + altitude=ALTITUDE, algorithm=sp.Algorithm.SPA, + )) + parallel = elapsed(lambda: sp.solar_position_inplace_threaded( + LATITUDE, LONGITUDE, times, *bufs, + altitude=ALTITUDE, algorithm=sp.Algorithm.SPA, + )) + print(f" {len(times):,} positions, SPA") + print(f" serial {serial * 1e3:8.1f} ms") + print(f" threaded {parallel * 1e3:8.1f} ms ({serial / parallel:.2f}x)") + + # The count is fixed once the runtime is up. A late call raises rather than + # silently doing nothing, so this mistake is never invisible. To use a different + # number of threads, start a new process. + try: + sp.set_num_threads(requested + 1) + except RuntimeError as err: + print(f"\n changing it afterwards is refused:\n {err}") + + +# -------------------------------------------------------------------------------------- +# 3. A single position +# -------------------------------------------------------------------------------------- + + +def single_section(when: dt.datetime) -> None: + rule("3. One instant") + # The datetime goes straight in. Unix seconds work just as well, and every other + # section here passes those, because they build their times by arithmetic. + r = sp.solar_position_single( + LATITUDE, LONGITUDE, when, altitude=ALTITUDE, algorithm=sp.Algorithm.SPA + ) + print(f" {when.isoformat()}") + print(f" azimuth {r.azimuth:10.6f}deg") + print(f" elevation {r.elevation:10.6f}deg") + print(f" zenith {r.zenith:10.6f}deg") + print(f" apparent elevation {r.apparent_elevation:10.6f}deg (refraction applied)") + print(f" apparent zenith {r.apparent_zenith:10.6f}deg") + + +# -------------------------------------------------------------------------------------- +# 4. Choosing an algorithm +# -------------------------------------------------------------------------------------- +# SPA is the most accurate and the slowest; PSA is the package default. `RefractionModel +# .DEFAULT` means "whatever this algorithm's own default is", which is no refraction for +# PSA, Walraven, Iqbal and USNO, and a model for NOAA, Michalsky and SPA. + + +def algorithms_section(t: float) -> None: + rule("4. Algorithms (true elevation, no refraction)") + reference = sp.solar_position_single( + LATITUDE, LONGITUDE, t, altitude=ALTITUDE, + algorithm=sp.Algorithm.SPA, refraction=sp.RefractionModel.NONE, + ).elevation + + for alg in sp.Algorithm: + r = sp.solar_position_single( + LATITUDE, LONGITUDE, t, altitude=ALTITUDE, + algorithm=alg, refraction=sp.RefractionModel.NONE, + ) + delta = r.elevation - reference + print(f" {alg.name:<10} {r.elevation:10.6f}deg vs SPA {delta:+.2e}deg") + + +# -------------------------------------------------------------------------------------- +# 5. Refraction models +# -------------------------------------------------------------------------------------- +# Refraction bends light downward, so the sun always appears higher than it truly is. +# The effect is tiny overhead and largest near the horizon, so this uses a low sun. + + +def refraction_section() -> None: + rule("5. Refraction models, sun near the horizon") + # Amsterdam, shortly after sunrise on the solstice. + t = to_unix(dt.datetime(2023, 6, 21, 3, 45, tzinfo=dt.timezone.utc)) + + true_elevation = sp.solar_position_single( + LATITUDE, LONGITUDE, t, altitude=ALTITUDE, + refraction=sp.RefractionModel.NONE, + ).elevation + print(f" true elevation {true_elevation:.6f}deg\n") + + for model in sp.RefractionModel: + if model is sp.RefractionModel.DEFAULT: + continue + r = sp.solar_position_single( + LATITUDE, LONGITUDE, t, altitude=ALTITUDE, refraction=model, + pressure=101325.0, temperature=12.0, + ) + bend = r.apparent_elevation - r.elevation + print(f" {model.name:<10} apparent {r.apparent_elevation:9.6f}deg" + f" bend {bend:+.6f}deg") + + +# -------------------------------------------------------------------------------------- +# 6. Many instants at once +# -------------------------------------------------------------------------------------- + + +def many_section() -> None: + rule("6. A day of positions") + start = to_unix(dt.datetime(2023, 6, 21, 0, 0, tzinfo=dt.timezone.utc)) + times = start + np.arange(0, 24, dtype=np.float64) * 3600.0 + + # Returns an (n, 5) array: azimuth, elevation, zenith, apparent elevation, + # apparent zenith. + m = sp.solar_position( + LATITUDE, LONGITUDE, times, altitude=ALTITUDE, algorithm=sp.Algorithm.SPA + ) + azimuth, elevation = m[:, 0], m[:, 1] + + print(" hour(UTC) azimuth elevation") + for hour in (3, 6, 9, 12, 15, 18, 21): + print(f" {hour:02d} {azimuth[hour]:8.3f} {elevation[hour]:9.3f}") + + daylight = elevation > 0 + print(f"\n highest elevation {elevation.max():.3f}deg " + f"at {from_unix(times[elevation.argmax()]).strftime('%H:%M')} UTC") + print(f" hours sampled above the horizon: {daylight.sum()} of 24") + + +# -------------------------------------------------------------------------------------- +# 7. Reusing your own buffers +# -------------------------------------------------------------------------------------- +# `solar_position` allocates a fresh array per call. In a loop over many sites, +# `..._into` writes straight into arrays you already hold and allocates nothing. The +# arrays must be float64, contiguous, and at least as long as the input. + + +def inplace_section() -> None: + rule("7. Writing into preallocated arrays") + start = to_unix(dt.datetime(2023, 6, 21, tzinfo=dt.timezone.utc)) + times = start + np.arange(0, 1440, dtype=np.float64) * 60.0 # every minute + n = len(times) + + azimuth = np.empty(n) + elevation = np.empty(n) + zenith = np.empty(n) + apparent_elevation = np.empty(n) + apparent_zenith = np.empty(n) + + sites = {"Amsterdam": (52.35888, 4.88185), "Nairobi": (-1.2921, 36.8219)} + for name, (lat, lon) in sites.items(): + sp.solar_position_inplace( + lat, lon, times, + azimuth, elevation, zenith, apparent_elevation, apparent_zenith, + algorithm=sp.Algorithm.SPA, + ) + print(f" {name:<10} peak elevation {elevation.max():7.3f}deg " + f"minutes of daylight {(elevation > 0).sum()}") + + +# -------------------------------------------------------------------------------------- +# 8. Sunrise, transit and sunset +# -------------------------------------------------------------------------------------- +# `transit_sunrise_sunset` returns seconds since midnight UTC of the day containing the +# instant you pass, keeping sub-second precision. Values can fall outside [0, 86400) +# when the event lands outside that UTC day. + + +def sunrise_section() -> None: + rule("8. Sunrise, transit, sunset") + day = dt.datetime(2023, 6, 21, tzinfo=dt.timezone.utc) + events = sp.transit_sunrise_sunset( + LATITUDE, LONGITUDE, to_unix(day), altitude=ALTITUDE + ) + + midnight = to_unix(day) + tz = ZoneInfo("Europe/Amsterdam") if ZoneInfo is not None else dt.timezone.utc + for name, seconds in ( + ("sunrise", events.sunrise), + ("transit", events.transit), + ("sunset", events.sunset), + ): + # Round, so the printed clock time matches the whole-second values the Julia + # API reports; the float carries the sub-second part. + when = from_unix(midnight + round(seconds), tz) + print(f" {name:<8} {when.strftime('%Y-%m-%d %H:%M:%S %Z')}" + f" ({seconds:9.3f} s after 00:00 UTC)") + + length = events.sunset - events.sunrise + print(f"\n day length {length / 3600:.3f} hours") + + +# -------------------------------------------------------------------------------------- +# 9. The next or previous event +# -------------------------------------------------------------------------------------- + + +def next_previous_section() -> None: + rule("9. Next and previous events") + now = to_unix(dt.datetime(2023, 6, 21, 12, 30, tzinfo=dt.timezone.utc)) + + for event in sp.SunEvent: + nxt = sp.sun_event( + LATITUDE, LONGITUDE, now, event, sp.EventDirection.NEXT, altitude=ALTITUDE + ) + prv = sp.sun_event( + LATITUDE, LONGITUDE, now, event, sp.EventDirection.PREVIOUS, + altitude=ALTITUDE, + ) + print(f" {event.name:<11} previous {from_unix(prv):%Y-%m-%d %H:%M:%S}" + f" next {from_unix(nxt):%Y-%m-%d %H:%M:%S} (UTC)") + + +# -------------------------------------------------------------------------------------- +# 10. Options and error handling +# -------------------------------------------------------------------------------------- + + +def options_section(t: float) -> None: + rule("10. Options and errors") + + # Every option is accepted by every entrypoint, but an algorithm that has no such + # setting ignores it. delta_t only reaches SPA and USNO: PSA and Iqbal have no + # delta_t field at all, and NOAA has one that its implementation never reads. + # So demonstrate it with SPA, and on a date where the model and the 67 s default + # diverge. + old = to_unix(dt.datetime(1900, 6, 21, 12, tzinfo=dt.timezone.utc)) + for label, value in (("67.0", 67.0), ("None", None)): + r = sp.solar_position_single( + LATITUDE, LONGITUDE, old, algorithm=sp.Algorithm.SPA, delta_t=value + ) + print(f" SPA 1900 delta_t={label:<4} elevation {r.elevation:.9f}deg") + print(" (delta_t is ignored by PSA, Iqbal and NOAA)") + + # PSA ships two coefficient sets. + for coeffs in (2020, 2001): + r = sp.solar_position_single(LATITUDE, LONGITUDE, t, psa_coeffs=coeffs) + print(f" psa_coeffs={coeffs} elevation {r.elevation:.6f}deg") + + # Julia exceptions arrive as JLWError with a code: 2 is ArgumentError, + # 3 is DimensionMismatch. + try: + times = np.zeros(10) + too_short = [np.zeros(3) for _ in range(5)] + sp.solar_position_inplace(LATITUDE, LONGITUDE, times, *too_short) + except sp.JLWError as err: + print(f"\n JLWError code={err.code}: {err.message}") + + # Invalid enum values are rejected rather than silently falling through. + try: + sp.solar_position_single(LATITUDE, LONGITUDE, t, algorithm=99) + except (ValueError, KeyError, TypeError) as err: + print(f" invalid algorithm rejected: {type(err).__name__}") + + +def main() -> None: + print(f"solarposition example - Amsterdam ({LATITUDE}N, {LONGITUDE}E, " + f"{ALTITUDE:.0f} m)") + when = times_section() + threaded_section() # must run before any computation: it sets the threads + t = to_unix(when) + single_section(when) + algorithms_section(t) + refraction_section() + many_section() + inplace_section() + sunrise_section() + next_previous_section() + options_section(t) + print() + + +if __name__ == "__main__": + main() diff --git a/lib/examples/simple.py b/lib/examples/simple.py new file mode 100644 index 00000000..5de27ce9 --- /dev/null +++ b/lib/examples/simple.py @@ -0,0 +1,27 @@ +"""Simple demo of solarposition.py bindings.""" + +import datetime as dt + +import solarposition as sp + +# number of threads. Needs to be set before calling any other functions. +sp.set_num_threads(16) + +# algorithm selection +ALG = sp.Algorithm.SPA +REF = sp.RefractionModel.BENNETT +print(f"Using algorithm {ALG.name} and refraction model {REF.name}") + +# location and time +LATITUDE = 40.0 +LONGITUDE = -105.0 +print(f"Location: {LATITUDE} N, {LONGITUDE} W") + +# times can be a datetime, a date, a numpy datetime64, or plain Unix seconds. A +# datetime without a tzinfo is read as UTC. +WHEN = dt.datetime(2023, 6, 21, 3, 45, tzinfo=dt.timezone.utc) + +pos = sp.solar_position_single(LATITUDE, LONGITUDE, WHEN, refraction=REF, algorithm=ALG) +print(f"Solar position at {WHEN.isoformat()}:") +print(f" Elevation: {pos.elevation:.2f} degrees") +print(f" Azimuth: {pos.azimuth:.2f} degrees") From 53ee8c9bc14e6df8e443565f62438b4f117e51b4 Mon Sep 17 00:00:00 2001 From: Stefan de Lange Date: Sun, 13 Sep 2026 02:00:22 +0200 Subject: [PATCH 08/13] Document the bindings in lib/README.md Co-Authored-By: Claude Opus 5 (1M context) --- lib/README.md | 341 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 341 insertions(+) create mode 100644 lib/README.md diff --git a/lib/README.md b/lib/README.md new file mode 100644 index 00000000..6d9481b7 --- /dev/null +++ b/lib/README.md @@ -0,0 +1,341 @@ +# libsolarposition — Python bindings for SolarPosition.jl + +A C-ABI shared library and Python package compiled from SolarPosition.jl with +[`juliac --trim`](https://github.com/JuliaLang/JuliaC.jl), wrapped by +[JuliaLibWrapping.jl](https://github.com/JuliaInterop/JuliaLibWrapping.jl). No Julia +installation is needed to use the result: the Julia runtime is bundled. + +Every positioning algorithm and every refraction model is exposed, plus the +sunrise/sunset utilities. Results are validated to match both the Julia package and the +solposx reference tables — see [Tests](#tests). + +## Requirements + +- **Julia 1.13+** to build. JuliaLibWrapping's ABI export requires it. The package + itself still supports Julia 1.10 LTS; this subproject is deliberately *not* part of the + root `[workspace]` so it cannot drag that requirement into the package's CI. +- Python 3.8+ and NumPy to use. +- A C compiler (`gcc`) — `juliac` compiles a small shim and links the library. + +## Build + +```bash +julia --project=lib -e 'using Pkg; Pkg.instantiate()' +julia --project=lib/build-env -e 'using Pkg; Pkg.instantiate()' + +julia --project=lib lib/build.jl # default: signal-safe, single-threaded +julia --project=lib lib/build.jl --threads=4 # bake in a default; usually unneeded +julia --project=lib lib/build.jl --verbose # pass juliac's output through + +pip install ./lib/out +``` + +Output lands in `lib/out/` (gitignored): + +| Path | What | +| --- | --- | +| `libsolarposition.h` | C header | +| `libsolarposition.abi.json` | ABI metadata emitted by `juliac` | +| `libsolarposition.jlw.json` | `@api` metadata sidecar | +| `solarposition/` | installable Python package, including a ~80 MB `bundle/` | +| `libsolarposition-bundle/` | the relocatable runtime bundle | + +The compiled library is 2.9 MB; the rest of the bundle is the Julia runtime. + +## Python API + +```python +import datetime as dt +import numpy as np +import solarposition as sp + +# One instant, as a datetime or as Unix seconds; both reach the same C double. +when = dt.datetime(2023, 6, 21, 12, tzinfo=dt.timezone.utc) +r = sp.solar_position_single(52.35888, 4.88185, when, altitude=100.0, + algorithm=sp.Algorithm.SPA) +r.azimuth, r.elevation, r.zenith, r.apparent_elevation, r.apparent_zenith + +# Many instants -> an (n, 5) array of the same five columns. +times = np.array([1687348800.0 + 3600 * i for i in range(24)]) +m = sp.solar_position(52.35888, 4.88185, times, algorithm=sp.Algorithm.SPA) + +# Or write into buffers you already hold, with no allocation. +az, el, ze, ael, aze = (np.zeros(len(times)) for _ in range(5)) +sp.solar_position_inplace(52.35888, 4.88185, times, az, el, ze, ael, aze) + +# Sunrise, transit and sunset, as seconds since midnight UTC. +e = sp.transit_sunrise_sunset(52.35888, 4.88185, when, altitude=100.0) +e.transit, e.sunrise, e.sunset + +# The next or previous event, as Unix seconds. +sp.sun_event(52.35888, 4.88185, 1687348800.0, + sp.SunEvent.SUNRISE, sp.EventDirection.NEXT) +``` + +### Times + +Anywhere a time is taken you may pass Unix seconds, a `datetime.datetime`, a +`datetime.date` (meaning midnight), a `numpy.datetime64`, or — for the batch +entrypoints — a `datetime64` array, a list of `datetime` objects, or a pandas +`DatetimeIndex` or `Series`, tz-aware or not. + +A `datetime` carrying a `tzinfo` has its offset applied. A naive one is read as **UTC**, +not as local time. That departs from `datetime.timestamp()` on purpose: everything this +library takes and returns is UTC, and `numpy.datetime64` carries no zone to read, so UTC +is the only rule that stays consistent across all the accepted types. + +The conversion happens in Python, in `_extras.py` — the C side takes `Float64` only. It +is vectorised for `datetime64` input and falls back to a per-element loop for object +arrays and tz-aware pandas containers, so a large tz-aware index is worth converting +once yourself if you are calling repeatedly. Numeric input is passed through untouched: +a wrong dtype or a non-contiguous array still raises rather than being silently copied +or reinterpreted. + +Results come back as Unix seconds; converting those back is +`datetime.datetime.fromtimestamp(seconds, datetime.timezone.utc)`. + +### Algorithms and options + +`Algorithm` is `PSA NOAA SPA WALRAVEN USNO IQBAL MICHALSKY`; `RefractionModel` is +`DEFAULT NONE HUGHES ARCHER BENNETT MICHALSKY SG2 SPA`. Both are `enum.IntEnum`s and are +validated on the way in — an out-of-range value raises rather than silently selecting an +algorithm. `DEFAULT` means each algorithm's own default, exactly as in Julia: none for +PSA, Walraven, Iqbal and USNO; HUGHES for NOAA, MICHALSKY for Michalsky, SPA for SPA. + +Algorithm and refraction options ride along as keyword arguments: `pressure`, +`temperature`, `delta_t` (pass `None` to compute it from the date), `atmos_refract`, +`refraction_limit`, `psa_coeffs`, `gmst_option`, `spencer_correction` and `julian_date`. + +Because there is one entrypoint per shape rather than one per algorithm, **every option is +accepted by every call and quietly ignored by algorithms that have no such setting**: + +| Option | Reaches | +| --- | --- | +| `pressure`, `temperature` | HUGHES, BENNETT, SG2, SPARefraction; the SPA algorithm | +| `delta_t` | SPA, USNO. *Not* PSA or Iqbal (no such field), and not NOAA (see below) | +| `atmos_refract` | the SPA algorithm | +| `refraction_limit` | SPARefraction | +| `psa_coeffs` | PSA (2001 or 2020) | +| `gmst_option` | USNO (1 or 2) | +| `spencer_correction`, `julian_date` | Michalsky | + +`lib/examples/getting_started.py` is a runnable tour of all of this. + +Every entrypoint returns all five angles. Where the selected refraction model applies no +correction, `apparent_*` equals the true value. That keeps the ABI return type uniform +across the `SolPos`/`ApparentSolPos` split, which is what lets each dispatch arm stay +inferable under `--trim`. + +Errors surface as `JLWError` with `.code` and `.message` — `2` for `ArgumentError`, `3` +for `DimensionMismatch`, and so on. + +## Threading + +Choose the thread count at run time, before the first computation: + +```python +import solarposition as sp + +sp.set_num_threads(8) # must precede any solar position call +sp.julia_nthreads() # -> 8 +sp.solar_position_inplace_threaded(lat, lon, times, *buffers) +``` + +`set_num_threads` returns what it actually configured and raises `RuntimeError` if that +does not match the request, which is what happens when the runtime is already running — +whether from an earlier computation or an earlier call. The count cannot change after +that; start a new process instead. `julia_nthreads()` reports the current value, and +`solar_position_inplace_threaded` runs the loop across those threads, giving +elementwise identical answers to the serial version. + +**The default is one thread, and raising it is a real trade-off.** More than one Julia +thread requires `handle-signals=yes`, because on Julia ≥ 1.12 a library with +`handle-signals=no` and more than one thread segfaults during shutdown (JuliaC's own +validation step rejects that combination). `set_num_threads(n > 1)` therefore turns +Julia's own SIGINT/SIGSEGV handlers on, which changes Ctrl-C behaviour for the whole +Python program. That is why it is an explicit call rather than a default. Upstream +tracking issue: [JuliaLang/julia#61319](https://github.com/JuliaLang/julia/issues/61319), +whose "critical fixes for even basic functionality" are still open. + +### How it works, and what does not work + +`JULIA_NUM_THREADS` has no effect. `juliac` compiles an `__attribute__((constructor))` +shim that calls `jl_parse_opts` with `--threads=` at `dlopen`, overriding the +environment. + +But the Julia runtime does **not** initialise at `dlopen` — it starts on the first call +into a `@ccallable` entrypoint. (Delete `libopenblas64_.so` from the bundle and +`import solarposition` still succeeds; the first computation is what dies.) That leaves a +window, and `lib/c/set_num_threads.c` uses it. + +Within that window, one approach works and one does not: + +| Attempt | Result | +| --- | --- | +| Assign `jl_options.nthreads` directly | write lands, `Threads.nthreads()` ignores it | +| Re-run `jl_parse_opts`, as JuliaC's own shim does | works | + +`jl_is_initialized()` is not usable as a guard here — it reads 0 even after a successful +call in an AOT-compiled library — which is why `set_num_threads` verifies its own outcome +instead of asking whether the runtime has started. + +Passing `--threads=N` to `lib/build.jl` still bakes a default in at compile time, which is +useful if you want a multi-threaded library without the caller having to ask. + +### Measured + +8 threads on a 32-core machine, one million instants, SPA, from the **default** +single-threaded build with `set_num_threads(8)`: + +| | time | | +| --- | --- | --- | +| serial | 3412.8 ms | | +| threaded | 430.2 ms | **7.93x** | + +Results are elementwise identical, and 20/20 import-compute-exit cycles at 8 threads +exited cleanly. + +## Tests + +```bash +julia --project=test lib/test/generate_fixtures.jl # writes lib/test/fixtures.json +julia --project=lib lib/build.jl +pip install ./lib/out +pytest lib/test -v +``` + +`fixtures.json` is generated (and gitignored) from two independent sources: + +- **The solposx reference tables** in `test/positioning/expected-values.jl` — the same + values the Julia test suite uses, on the 84.375 s grid where the full Julian Date is + exactly representable in Float64. Compared at the Julia suite's own tolerances + (`1e-10`, `1e-8` for SPA). Each table carries its call spec as data, because they are + not uniform: NOAA is validated against HUGHES at temperature 10.0 and skips the poles, + USNO has two GMST variants, and SPA has a column this API does not expose. +- **This checkout's own Julia output** for all 7 × 8 algorithm/refraction pairs at the + same 19 instants (poles, years 1800 and 2200 included), plus sunrise/sunset. Compared + at `1e-12`, which isolates anything `--trim` changed about the numerics. + +Both pass. The refraction decoupling described below was verified this way: all 168 +algorithm × refraction × instant combinations are bit-identical to +`solar_position(obs, dt, alg, refraction)`. + +## Design notes + +**Refraction is decoupled.** `Positioning._solar_position(obs, dt, alg, refraction)` +computes the geometric position and then adds `Refraction.refraction(model, elevation)`. +So the 7 × 8 matrix needs 7 geometric arms plus 8 refraction arms rather than 56 compiled +arms. `DEFAULT` is the exception and delegates to the package per algorithm, because +SPA's default path feeds the SPA struct's own pressure, temperature and `atmos_refract` +and cannot be rebuilt from the decoupled form. + +**The enums live in nested modules.** `@enum` binds its members in the enclosing module, +and `MICHALSKY` and `SPA` each name both an algorithm and a refraction model. Nesting +avoids the collision while still giving Python the clean `Algorithm.SPA` / +`RefractionModel.SPA` spelling. + +**Time crosses as Unix seconds**, not a `DateTime`, so nothing on the compiled path +touches TimeZones.jl. Convert on the Python side. + +**No `horizon` parameter.** `src/Utilities/spa.jl` hardcodes the horizon depression to +−0.8333° and never reads `Observer.horizon`; the field is inert throughout `src/`. A +`horizon` keyword here would be silently ignored, which matters most in exactly the case +someone would reach for it (civil twilight). Add one when the package honours the field. + +**`Interpolated` and `solar_rate` are not exposed.** They need Interpolations.jl on the +compiled path, which is very unlikely to survive `--trim`. + +## Known constraints + +Ordered by how much they cost. + +1. **The wheel is not publishable as-is.** `auditwheel` reports the wheel as + `linux_x86_64`, constrained by glibc 2.38 symbols from the build host, and PyPI only + accepts `manylinux*`. Publishing requires building inside a manylinux container (e.g. + `quay.io/pypa/manylinux_2_28_x86_64`) with Julia installed. The generated + `pyproject.toml` also declares `py3-none-any`, which is wrong for a binary wheel, and + carries a placeholder version. See + [JuliaLibWrapping#38](https://github.com/JuliaInterop/JuliaLibWrapping.jl/issues/38). + The built wheel is 72 MB, under PyPI's 100 MB per-file cap but without much headroom. +2. **~103 MB of the bundle is OpenBLAS, and it cannot be removed.** SolarPosition uses no + linear algebra, but `OpenBLAS_jll.__init__` still runs when the trimmed library + initialises and hard-fails with `could not load library "libopenblas64_.so"` if the + file is absent. Deleting libgfortran, libquadmath, libgomp or libblastrampoline hits + the same error. PackageCompiler also stores each library three times (`libfoo.so`, + `.so.N`, `.so.N.M`) — in `lib/out` those are links, but `pip install` materialises + them into 231 MB on disk. +3. **Two entrypoints must take positional arguments.** Reached through `@api`'s + keyword-argument wrapper, the `--trim` verifier cannot resolve the call + (`unresolved call from statement Core.kwcall(...)`) for + `solar_position_inplace_threaded` (`Threads.@threads` builds a closure) or + `transit_sunrise_sunset` (the seconds-returning path takes a different + `_frac_to_event` specialisation than the `DateTime` one, and only that one trips the + verifier — `sun_event` keeps its keywords). `lib/python/_extras.py` restores the + keyword-only Python signatures by hand. +4. **`solar_position` cannot use its generated binding.** Some builds of libffi + marshal that one signature wrongly and pass `delta_t` where `latitude` belongs, so + every result is silently for the wrong place. Calling the same `.so` from C is + correct, and so is every other entrypoint, so this is neither a Julia nor a + JuliaLibWrapping fault. The trigger is an argument classified across two classes — + `COpt_Float64` is an `int32` flag plus a `double` — whose integer half lands in the + *last* free general-purpose register; libffi then puts its `double` half in `xmm0`, + on top of the first floating-point argument. Only this entrypoint has that shape: + the returned matrix forces a hidden return pointer, and with the input vector and + the two enums that fills the register file exactly. Seen with the libffi bundled in + python-build-standalone 3.13.7, which is what `uv` installs by default, and not with + Fedora's libffi 3.5.2 — so a version check cannot route around it. + `lib/python/_extras.py` reimplements the function over + `solar_position_inplace`, whose signature is unaffected and which is the faster + path anyway. `test_delta_t_does_not_leak_into_other_arguments` guards the whole + class: `delta_t` is the only split argument, and three algorithms ignore it, so + changing it must change nothing. +5. **JuliaLibWrapping forwards neither `jl_options` nor `c_sources`** to + `JuliaC.ImageRecipe`, so neither the thread defaults nor the load-time thread shim can + go through `build_library`. `build.jl` therefore runs the normal `standard_build` for + the wrappers and ABI, then recompiles the library and bundle through `JuliaC` + directly. Only the `.so` differs between thread configurations, so the wrappers are + untouched. Delete `rebuild_library` once a passthrough exists upstream. +6. **One wrapped library per process.** The generated `_lowlevel.py` warns if a second + JuliaLibWrapping-wrapped package is imported into the same process, because the + dynamic linker would share a single `libjulia` between them. +7. **Keep virtualenvs out of `lib/`.** `juliac` copies the project tree into a temporary + directory and chokes on the executable bits inside a `.venv`. + +## Upstream issues worth filing + +- `jl_options` passthrough for JuliaLibWrapping's `build_library` (constraint 5). Small + change; would delete `build_threaded_library` here. +- `@api` keyword wrappers defeat the `--trim` verifier for callees containing closures + (constraint 3). +- `juliac --trim` still runs `OpenBLAS_jll.__init__` for a library that uses no BLAS + (constraint 2) — worth roughly 103 MB per bundle. +- `experiments/julia-issue-handle-signals-threads.md` in this repo is a complete, + unfiled regression report for the `handle-signals=no` shutdown segfault, with crash + tables across Julia 1.10–1.13 and a gdb backtrace. +- **libffi mismarshals a split-class struct argument that exhausts the integer + registers** (constraint 4). Reproduced in twelve lines of C plus `ctypes`, with no + Julia involved: a callee taking `(double, double, {int64,ptr}, double, int32, int32, + double, double, {int32,double}, double, double)` receives the struct's `double` half + in place of its first argument. Correct on Fedora's libffi 3.5.2, wrong on the libffi + inside python-build-standalone 3.13.7. Removing any one of the integer arguments, or + making the return small enough not to need a return pointer, makes it correct again — + which places the fault in the register-exhaustion path of the classifier. Worth + filing against libffi, and against python-build-standalone so its bundled copy moves. +- **`ARCHER` refraction throws `DomainError` near −2.7° elevation, and is nonsense + either side of it.** `src/Refraction/archer.jl:55` divides by + `0.955 + 20.267·cos(Z)`, which crosses zero at an elevation of about −2.70°, so the + argument to `acosd` leaves `[-1, 1]`: −2.69°, −2.70° and −2.71° all throw, while + −2.60° and −2.80° return corrections of +13.39° and −13.83° for a quantity that never + exceeds about 0.6°. Reproduced in plain Julia, so it is not a bindings artefact. The + package's own rule is to guard inverse trigonometry with `unit_clamp`, and this call + is unguarded — but clamping alone would convert the throw into those same wrong + values, so the real question is the elevation below which Archer's fit stops being + usable. In the Python suite it surfaces as three of the 56 algorithm/refraction pairs + failing whenever a sampled elevation happens to land near the pole. +- `Observer.horizon` is documented but never used (design note above). +- **`NOAA.delta_t` is ignored by its own implementation.** `src/Positioning/noaa.jl:38-42` + computes `δt` from the field (or from `calculate_deltat`) and then never reads it — it + is a dead local. Measured through these bindings, `NOAA(67.0)` and `NOAA(nothing)` agree + to 0.0e+00 degrees at 1900, 2000 and 2050, where SPA moves by ~1e-4. The Julia test + that covers this case (`test/positioning/test-noaa.jl:57`) compares at `atol = 1.0` + *degree*, which is far too loose to notice. From 59058ca8fca318ad22c5852fcf508aac9364ffd8 Mon Sep 17 00:00:00 2001 From: Stefan de Lange Date: Sun, 13 Sep 2026 02:00:22 +0200 Subject: [PATCH 09/13] Add changelog entry for the Python bindings Co-Authored-By: Claude Opus 5 (1M context) --- CHANGELOG.md | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index fc20ffb4..5fbed3c0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,20 @@ Releases before v0.5.0 predate this file. See the ## unreleased +### Added + +- Python bindings, as a `lib/` subproject compiled with `juliac --trim` and wrapped by + [JuliaLibWrapping.jl](https://github.com/JuliaInterop/JuliaLibWrapping.jl). Every + positioning algorithm and refraction model is exposed, along with sunrise, sunset and + transit, and an opt-in multithreaded build. The bindings are validated against both the + solposx reference tables and this package's own output for all 56 algorithm/refraction + pairs, and the thread count is selectable at run time via `set_num_threads()`. + Times may be given as `datetime`, `date`, `numpy.datetime64` or a pandas + `DatetimeIndex`/`Series` as well as Unix seconds; naive datetimes are read as UTC. + Nothing in `src/` changed. Building needs Julia 1.13, so `lib/` is intentionally + outside the root `[workspace]` and the package keeps its 1.10 LTS floor. See + `lib/README.md` + ### Changed - Documentation code blocks now render with syntax highlighting, line numbers with From b0eb42f24492a3ba7021c5d5d075a6d6fcdea959 Mon Sep 17 00:00:00 2001 From: Stefan de Lange Date: Sun, 13 Sep 2026 02:19:35 +0200 Subject: [PATCH 10/13] Add a solar_position_threaded entrypoint Co-Authored-By: Claude Opus 5 (1M context) --- lib/src/libsolarposition.jl | 79 ++++++++++++++++++++++++++++++++++++- 1 file changed, 77 insertions(+), 2 deletions(-) diff --git a/lib/src/libsolarposition.jl b/lib/src/libsolarposition.jl index ce958da7..1df927bd 100644 --- a/lib/src/libsolarposition.jl +++ b/lib/src/libsolarposition.jl @@ -532,6 +532,28 @@ function _fill_threaded!( return nothing end +# The same barrier for the allocating entrypoint, writing the five columns of one matrix +# instead of five separate buffers. +function _fill_threaded!( + obs::SP.Observer{Float64}, + n::Int, + unix_seconds::Vector{Float64}, + out::Matrix{Float64}, + algorithm::Algorithm, + refraction::RefractionModel, + o::Opts, + ) + Threads.@threads for i in 1:n + a = _compute(obs, unix2datetime(unix_seconds[i]), algorithm, refraction, o) + out[i, 1] = a.azimuth + out[i, 2] = a.elevation + out[i, 3] = a.zenith + out[i, 4] = a.apparent_elevation + out[i, 5] = a.apparent_zenith + end + return nothing +end + """ Threaded `solar_position_inplace`, using the threads the library was compiled with. @@ -601,8 +623,61 @@ end julian_date::JulianDateMode, )::Nothing -# `solar_position` hands an owned array to the caller, so the generated bindings -# need the release entrypoints to give it back. +""" +Threaded `solar_position`, returning a freshly allocated `n x 5` matrix. + +Identical to `solar_position` except that the loop runs across the Julia threads the +library was compiled with; with a single-threaded build it is the serial loop. Declared +with positional arguments for the same reason as `solar_position_inplace_threaded`. +""" +function solar_position_threaded( + latitude::Float64, + longitude::Float64, + unix_seconds::Vector{Float64}, + altitude::Float64, + algorithm::Algorithm, + refraction::RefractionModel, + pressure::Float64, + temperature::Float64, + delta_t::Union{Float64, Nothing}, + atmos_refract::Float64, + refraction_limit::Float64, + psa_coeffs::Int64, + gmst_option::Int64, + spencer_correction::Bool, + julian_date::JulianDateMode, + ) + n = length(unix_seconds) + obs = SP.Observer(latitude, longitude, altitude) + o = Opts( + pressure, temperature, delta_t, atmos_refract, refraction_limit, + psa_coeffs, gmst_option, spencer_correction, julian_date, + ) + out = Matrix{Float64}(undef, n, 5) + _fill_threaded!(obs, n, unix_seconds, out, algorithm, refraction, o) + return out +end + +@api solar_position_threaded( + latitude::Float64, + longitude::Float64, + unix_seconds::Vector{Float64}, + altitude::Float64, + algorithm::Algorithm, + refraction::RefractionModel, + pressure::Float64, + temperature::Float64, + delta_t::Union{Float64, Nothing}, + atmos_refract::Float64, + refraction_limit::Float64, + psa_coeffs::Int64, + gmst_option::Int64, + spencer_correction::Bool, + julian_date::JulianDateMode, +)::Matrix{Float64} + +# `solar_position` and `solar_position_threaded` hand owned arrays to the caller, so the +# generated bindings need the release entrypoints to give them back. @export_release_entrypoints end # module From 163959117d95b526b951770385282f68a6c75458 Mon Sep 17 00:00:00 2001 From: Stefan de Lange Date: Sun, 13 Sep 2026 02:19:35 +0200 Subject: [PATCH 11/13] Wrap solar_position_threaded in Python Co-Authored-By: Claude Opus 5 (1M context) --- lib/python/_extras.py | 74 ++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 69 insertions(+), 5 deletions(-) diff --git a/lib/python/_extras.py b/lib/python/_extras.py index 333565e8..db833bdd 100644 --- a/lib/python/_extras.py +++ b/lib/python/_extras.py @@ -3,12 +3,13 @@ `build.jl` copies this module into the generated Python package and appends an import to `_facade.py` so these definitions override the generated ones. -Why this file exists: two entrypoints must be declared with all-positional arguments -the Julia side, because reached through `@api`'s keyword-argument wrapper the +Why this file exists: three entrypoints must be declared with all-positional arguments +on the Julia side, because reached through `@api`'s keyword-argument wrapper the `juliac --trim` verifier cannot resolve the call ("unresolved call from statement Core.kwcall(...)") and the build fails. -- ``solar_position_inplace_threaded`` — ``Threads.@threads`` builds a closure. +- ``solar_position_inplace_threaded`` and ``solar_position_threaded`` — + ``Threads.@threads`` builds a closure. - ``transit_sunrise_sunset`` — the seconds-returning sunrise/sunset path takes a different ``_frac_to_event`` specialisation than the ``DateTime`` one, and only the former trips the verifier. ``sun_event`` goes through the ``DateTime`` path and keeps @@ -17,8 +18,9 @@ Declaring those arguments positionally keeps the entrypoints trimmable, and the wrappers here restore the keyword-only signatures the rest of the API has. -``solar_position`` is here for an unrelated reason: the generated binding for it -is miscompiled by some builds of libffi. See its docstring. +``solar_position`` is here for a second reason: the generated binding for it is +miscompiled by some builds of libffi. See its docstring — ``solar_position_threaded`` +has the same signature shape and takes the same detour. Every entrypoint that takes a time accepts datetimes as well as Unix seconds — see ``_unix_seconds``. The C side takes only ``Float64``, so that conversion has to happen @@ -49,6 +51,7 @@ "solar_position_inplace", "solar_position_inplace_threaded", "solar_position_single", + "solar_position_threaded", "sun_event", "transit_sunrise_sunset", ] @@ -497,3 +500,64 @@ def solar_position( _enum_coerce(JulianDateMode, julian_date), ) return out + + +def solar_position_threaded( + latitude, + longitude, + unix_seconds, + *, + altitude=0.0, + algorithm=Algorithm.PSA, + refraction=RefractionModel.DEFAULT, + pressure=101325.0, + temperature=12.0, + delta_t=67.0, + atmos_refract=0.5667, + refraction_limit=-0.5667, + psa_coeffs=2020, + gmst_option=1, + spencer_correction=True, + julian_date=JulianDateMode.ORIGINAL, +): + """Solar position for many instants, computed in parallel into a new array. + + ``solar_position`` and ``solar_position_inplace_threaded`` in one: the columns are + azimuth, elevation, zenith, apparent elevation and apparent zenith, the result is + allocated for you, and the loop runs across the Julia threads the library was + compiled with. Call ``julia_nthreads()`` to find out how many that is; a + single-threaded build runs this serially and gives the same answers. + + ``unix_seconds`` accepts the same datetime types as ``solar_position``, and the + keyword arguments are ``solar_position_single``'s. + + Like ``solar_position``, this routes through the in-place entrypoint rather than the + generated binding for the allocating one, because that signature is miscompiled by + some builds of libffi. See ``solar_position`` for the details. + """ + times = _unix_seconds_array(unix_seconds) + n = len(times) + out = np.empty((n, 5), dtype=np.float64, order="F") + _lowlevel.libsolarposition_solar_position_inplace_threaded( + latitude, + longitude, + CVector_borrowed_Float64.from_numpy(times), + CVector_borrowed_Float64.from_numpy(out[:, 0]), + CVector_borrowed_Float64.from_numpy(out[:, 1]), + CVector_borrowed_Float64.from_numpy(out[:, 2]), + CVector_borrowed_Float64.from_numpy(out[:, 3]), + CVector_borrowed_Float64.from_numpy(out[:, 4]), + altitude, + _enum_coerce(Algorithm, algorithm), + _enum_coerce(RefractionModel, refraction), + pressure, + temperature, + COpt_Float64.from_optional(delta_t), + atmos_refract, + refraction_limit, + psa_coeffs, + gmst_option, + spencer_correction, + _enum_coerce(JulianDateMode, julian_date), + ) + return out From f181e16a37c587bf05e12735efe04f5b0caf5cc0 Mon Sep 17 00:00:00 2001 From: Stefan de Lange Date: Sun, 13 Sep 2026 02:19:35 +0200 Subject: [PATCH 12/13] Test solar_position_threaded Co-Authored-By: Claude Opus 5 (1M context) --- lib/test/test_solarposition.py | 51 ++++++++++++++++++++--------- lib/test/test_threads_subprocess.py | 5 ++- 2 files changed, 39 insertions(+), 17 deletions(-) diff --git a/lib/test/test_solarposition.py b/lib/test/test_solarposition.py index 2122e2ae..9a2e89df 100644 --- a/lib/test/test_solarposition.py +++ b/lib/test/test_solarposition.py @@ -214,7 +214,8 @@ def test_delta_t_does_not_leak_into_other_arguments(algorithm): passed as a struct split across an integer and a floating-point half, which is exactly the shape ctypes/libffi is most likely to place wrongly. One such libffi bug silently put `delta_t` where `latitude` belongs on `solar_position`; - see `solar_position` in `lib/python/_extras.py`. + see `solar_position` in `lib/python/_extras.py`. `solar_position_threaded` has + the same signature shape, so it is checked here too. """ lat, lon, alt = 52.35888, 4.88185, 100.0 times = np.linspace(1687348800.0, 1687435200.0, 64) @@ -230,33 +231,41 @@ def sample(delta_t): threaded = [np.empty(len(times)) for _ in range(5)] sp.solar_position_inplace_threaded(lat, lon, times, *threaded, delta_t=delta_t, **kwargs) - return one.azimuth, many, inplace, threaded + many_threaded = sp.solar_position_threaded(lat, lon, times, + delta_t=delta_t, **kwargs) + return one.azimuth, many, inplace, threaded, many_threaded - one_a, many_a, inplace_a, threaded_a = sample(67.0) - one_b, many_b, inplace_b, threaded_b = sample(0.0) + one_a, many_a, inplace_a, threaded_a, many_threaded_a = sample(67.0) + one_b, many_b, inplace_b, threaded_b, many_threaded_b = sample(0.0) assert one_a == one_b, "solar_position_single" np.testing.assert_array_equal(many_a, many_b, err_msg="solar_position") + np.testing.assert_array_equal(many_threaded_a, many_threaded_b, + err_msg="solar_position_threaded") for col, (a, b) in enumerate(zip(inplace_a, inplace_b)): np.testing.assert_array_equal(a, b, err_msg=f"inplace column {col}") for col, (a, b) in enumerate(zip(threaded_a, threaded_b)): np.testing.assert_array_equal(a, b, err_msg=f"threaded column {col}") - # And all four must agree with each other, not merely be self-consistent. + # And all of them must agree with each other, not merely be self-consistent. np.testing.assert_array_equal(many_a[:, 0], inplace_a[0]) np.testing.assert_array_equal(many_a[:, 0], threaded_a[0]) + np.testing.assert_array_equal(many_a, many_threaded_a) assert many_a[0, 0] == one_a -def test_solar_position_uses_the_hand_maintained_wrapper(): - """`solar_position` must come from `_extras`, not from the generated façade. +@pytest.mark.parametrize( + "entrypoint", [sp.solar_position, sp.solar_position_threaded] +) +def test_allocating_entrypoints_use_the_hand_maintained_wrapper(entrypoint): + """These must come from `_extras`, not from the generated façade. - The generated binding for this one entrypoint is mismarshalled by some builds of - libffi, so the override is the fix rather than a convenience. If a build stops - installing it, every result from this function silently becomes wrong, and this + The generated bindings for the two allocating entrypoints are mismarshalled by + some builds of libffi, so the override is the fix rather than a convenience. If a + build stops installing it, every result from them silently becomes wrong, and this is the tripwire for that. """ - assert sp.solar_position.__module__.endswith("_extras") + assert entrypoint.__module__.endswith("_extras") def test_inplace_rejects_short_buffers(): @@ -276,8 +285,11 @@ def test_inplace_writes_through_to_caller_arrays(): assert np.all(np.isfinite(buf)) -def test_empty_input_is_allowed(): - m = sp.solar_position(40.0, -105.0, np.zeros(0)) +@pytest.mark.parametrize( + "entrypoint", [sp.solar_position, sp.solar_position_threaded] +) +def test_empty_input_is_allowed(entrypoint): + m = entrypoint(40.0, -105.0, np.zeros(0)) assert m.shape == (0, 5) @@ -327,9 +339,12 @@ def _hourly(n=5): return seconds, [dt.datetime.fromtimestamp(t, dt.timezone.utc) for t in seconds] -def test_solar_position_accepts_datetime_sequences(): +@pytest.mark.parametrize( + "entrypoint", [sp.solar_position, sp.solar_position_threaded] +) +def test_solar_position_accepts_datetime_sequences(entrypoint): seconds, datetimes = _hourly() - ref = sp.solar_position(40.0, -105.0, seconds) + ref = entrypoint(40.0, -105.0, seconds) for label, times in [ ("list of datetime", datetimes), ("tuple of datetime", tuple(datetimes)), @@ -341,7 +356,7 @@ def test_solar_position_accepts_datetime_sequences(): ), ("object array", np.array(datetimes, dtype=object)), ]: - got = sp.solar_position(40.0, -105.0, times) + got = entrypoint(40.0, -105.0, times) np.testing.assert_array_equal(got, ref, err_msg=label) @@ -551,6 +566,10 @@ def test_threaded_matches_serial(fx): np.testing.assert_array_equal( buf, serial[:, col], err_msg=f"{alg} col {col}" ) + allocating = sp.solar_position_threaded( + lat, lon, times, altitude=alt, algorithm=alg + ) + np.testing.assert_array_equal(allocating, serial, err_msg=str(alg)) def test_threaded_rejects_short_buffers(): diff --git a/lib/test/test_threads_subprocess.py b/lib/test/test_threads_subprocess.py index a1571bc3..f8238f76 100644 --- a/lib/test/test_threads_subprocess.py +++ b/lib/test/test_threads_subprocess.py @@ -50,7 +50,10 @@ def test_threaded_results_match_serial_with_runtime_threads(): " algorithm=sp.Algorithm.SPA)\n" "sp.solar_position_inplace_threaded(52.35888, 4.88185, t, *b," " algorithm=sp.Algorithm.SPA)\n" - "print(all(np.array_equal(x, y) for x, y in zip(a, b)))\n" + "c = sp.solar_position_threaded(52.35888, 4.88185, t," + " algorithm=sp.Algorithm.SPA)\n" + "print(all(np.array_equal(x, y) for x, y in zip(a, b))" + " and all(np.array_equal(x, c[:, i]) for i, x in enumerate(a)))\n" ) assert proc.returncode == 0, proc.stderr assert proc.stdout.strip() == "True" From 766028a2d336336126e77d62ef8a7c5b69069a92 Mon Sep 17 00:00:00 2001 From: Stefan de Lange Date: Sun, 13 Sep 2026 02:19:35 +0200 Subject: [PATCH 13/13] Document solar_position_threaded Co-Authored-By: Claude Opus 5 (1M context) --- lib/README.md | 36 ++++++++++++++++++--------------- lib/examples/getting_started.py | 11 ++++++++-- 2 files changed, 29 insertions(+), 18 deletions(-) diff --git a/lib/README.md b/lib/README.md index 6d9481b7..07638bb5 100644 --- a/lib/README.md +++ b/lib/README.md @@ -138,15 +138,17 @@ import solarposition as sp sp.set_num_threads(8) # must precede any solar position call sp.julia_nthreads() # -> 8 -sp.solar_position_inplace_threaded(lat, lon, times, *buffers) +sp.solar_position_threaded(lat, lon, times) # -> a new (n, 5) array +sp.solar_position_inplace_threaded(lat, lon, times, *buffers) # into your buffers ``` `set_num_threads` returns what it actually configured and raises `RuntimeError` if that does not match the request, which is what happens when the runtime is already running — whether from an earlier computation or an earlier call. The count cannot change after that; start a new process instead. `julia_nthreads()` reports the current value, and -`solar_position_inplace_threaded` runs the loop across those threads, giving -elementwise identical answers to the serial version. +`solar_position_threaded` and `solar_position_inplace_threaded` run the loop across +those threads, giving elementwise identical answers to their serial counterparts +`solar_position` and `solar_position_inplace`. **The default is one thread, and raising it is a real trade-off.** More than one Julia thread requires `handle-signals=yes`, because on Julia ≥ 1.12 a library with @@ -264,29 +266,31 @@ Ordered by how much they cost. the same error. PackageCompiler also stores each library three times (`libfoo.so`, `.so.N`, `.so.N.M`) — in `lib/out` those are links, but `pip install` materialises them into 231 MB on disk. -3. **Two entrypoints must take positional arguments.** Reached through `@api`'s +3. **Three entrypoints must take positional arguments.** Reached through `@api`'s keyword-argument wrapper, the `--trim` verifier cannot resolve the call (`unresolved call from statement Core.kwcall(...)`) for - `solar_position_inplace_threaded` (`Threads.@threads` builds a closure) or + `solar_position_threaded` and `solar_position_inplace_threaded` + (`Threads.@threads` builds a closure) or `transit_sunrise_sunset` (the seconds-returning path takes a different `_frac_to_event` specialisation than the `DateTime` one, and only that one trips the verifier — `sun_event` keeps its keywords). `lib/python/_extras.py` restores the keyword-only Python signatures by hand. -4. **`solar_position` cannot use its generated binding.** Some builds of libffi - marshal that one signature wrongly and pass `delta_t` where `latitude` belongs, so - every result is silently for the wrong place. Calling the same `.so` from C is - correct, and so is every other entrypoint, so this is neither a Julia nor a - JuliaLibWrapping fault. The trigger is an argument classified across two classes — +4. **The two allocating entrypoints cannot use their generated bindings.** Some builds + of libffi marshal that signature wrongly and pass `delta_t` where `latitude` + belongs, so every result is silently for the wrong place. Calling the same `.so` + from C is correct, and so is every other entrypoint, so this is neither a Julia nor + a JuliaLibWrapping fault. The trigger is an argument classified across two classes — `COpt_Float64` is an `int32` flag plus a `double` — whose integer half lands in the *last* free general-purpose register; libffi then puts its `double` half in `xmm0`, - on top of the first floating-point argument. Only this entrypoint has that shape: - the returned matrix forces a hidden return pointer, and with the input vector and - the two enums that fills the register file exactly. Seen with the libffi bundled in + on top of the first floating-point argument. Only `solar_position` and + `solar_position_threaded` have that shape: the returned matrix forces a hidden return + pointer, and with the input vector and the two enums that fills the register file + exactly. Both are wrong by the same 14.53° here. Seen with the libffi bundled in python-build-standalone 3.13.7, which is what `uv` installs by default, and not with Fedora's libffi 3.5.2 — so a version check cannot route around it. - `lib/python/_extras.py` reimplements the function over - `solar_position_inplace`, whose signature is unaffected and which is the faster - path anyway. `test_delta_t_does_not_leak_into_other_arguments` guards the whole + `lib/python/_extras.py` reimplements both over their in-place counterparts, whose + signatures are unaffected and which are the faster path anyway. + `test_delta_t_does_not_leak_into_other_arguments` guards the whole class: `delta_t` is the only split argument, and three algorithms ignore it, so changing it must change nothing. 5. **JuliaLibWrapping forwards neither `jl_options` nor `c_sources`** to diff --git a/lib/examples/getting_started.py b/lib/examples/getting_started.py index 55edc583..11a1392b 100644 --- a/lib/examples/getting_started.py +++ b/lib/examples/getting_started.py @@ -109,9 +109,16 @@ def elapsed(fn) -> float: LATITUDE, LONGITUDE, times, *bufs, altitude=ALTITUDE, algorithm=sp.Algorithm.SPA, )) + # solar_position_threaded is the same loop returning a new (n, 5) array, for when + # you do not already hold the buffers. + allocating = elapsed(lambda: sp.solar_position_threaded( + LATITUDE, LONGITUDE, times, + altitude=ALTITUDE, algorithm=sp.Algorithm.SPA, + )) print(f" {len(times):,} positions, SPA") - print(f" serial {serial * 1e3:8.1f} ms") - print(f" threaded {parallel * 1e3:8.1f} ms ({serial / parallel:.2f}x)") + print(f" serial {serial * 1e3:8.1f} ms") + print(f" threaded {parallel * 1e3:8.1f} ms ({serial / parallel:.2f}x)") + print(f" allocating {allocating * 1e3:8.1f} ms (threaded, new array)") # The count is fixed once the runtime is up. A late call raises rather than # silently doing nothing, so this mistake is never invisible. To use a different