From 65935096e739ab0054b06cc755c53c7776509bf7 Mon Sep 17 00:00:00 2001 From: Stefan de Lange Date: Wed, 3 Jun 2026 13:17:00 +0200 Subject: [PATCH] Add Iqbal, Michalsky and SG2 positioning algorithms Implement three additional solar position algorithms: - Iqbal: Fourier-series declination and equation of time (no refraction). - Michalsky: Astronomical Almanac approximate algorithm, with options for the Spencer (1989) azimuth correction and the Julian date formulation; defaults to the MICHALSKY refraction model. - SG2: fast multi-decadal algorithm valid for 1980-2030 (throws outside that range); defaults to the SG2 refraction model. To free the SG2 name for the positioning algorithm, the SG2 refraction model is renamed SG2 -> SG2Refraction, mirroring the existing SPA / SPARefraction split. Enables the previously skipped iqbal/michalsky/sg2 test stubs and updates the package and documentation tables, references and bibliography. Co-Authored-By: Claude Opus 4.8 (1M context) --- README.md | 3 + docs/src/index.md | 5 +- docs/src/positioning.md | 40 +++++++++ docs/src/refraction.md | 8 +- docs/src/refs.bib | 19 +++++ src/Positioning/Positioning.jl | 6 +- src/Positioning/iqbal.jl | 72 ++++++++++++++++ src/Positioning/michalsky.jl | 117 +++++++++++++++++++++++++ src/Positioning/sg2.jl | 132 +++++++++++++++++++++++++++++ src/Refraction/Refraction.jl | 2 +- src/Refraction/sg2.jl | 19 +++-- test/extensions/test-mtk.jl | 2 +- test/positioning/test-interface.jl | 4 +- test/positioning/test-iqbal.jl | 38 ++++----- test/positioning/test-michalsky.jl | 66 +++++++++------ test/positioning/test-psa.jl | 2 +- test/positioning/test-sg2.jl | 96 ++++++++------------- test/refraction/test-refraction.jl | 4 +- 18 files changed, 505 insertions(+), 130 deletions(-) create mode 100644 src/Positioning/iqbal.jl create mode 100644 src/Positioning/michalsky.jl create mode 100644 src/Positioning/sg2.jl diff --git a/README.md b/README.md index 6744942..2d6e48c 100644 --- a/README.md +++ b/README.md @@ -98,6 +98,9 @@ accuracy and implementation status. | Walraven | [Walraven, 1978]() | ±0.0100° | None | ✅ | | USNO | [U.S. Naval Observatory](https://aa.usno.navy.mil/faq/sun_approx) | ±0.0500° | None | ✅ | | SPA | [Reda & Andreas, 2004](https://doi.org/10.1016/j.solener.2003.12.003) | ±0.0003° | Built-in | ✅ | +| Iqbal | [Iqbal, 1983](https://doi.org/10.1016/B978-0-12-373750-2.X5001-0) | ±0.0100° | None | ✅ | +| Michalsky | [Michalsky, 1988]() | ±0.0100° | MICHALSKY | ✅ | +| SG2 | [Blanc & Wald, 2012](https://doi.org/10.1016/j.solener.2012.07.018) | ±0.0030° | SG2Refraction | ✅ | ## Refraction correction algorithms diff --git a/docs/src/index.md b/docs/src/index.md index 3b37486..2371fad 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -85,6 +85,9 @@ accuracy and implementation status. | [`Walraven`](@ref SolarPosition.Positioning.Walraven) | [Walraven, 1978](https://doi.org/10.1016/0038-092X(78)90155-X) | ±0.0100° | None | ✅ | | [`USNO`](@ref SolarPosition.Positioning.USNO) | [U.S. Naval Observatory](https://aa.usno.navy.mil/faq/sun_approx) | ±0.0500° | None | ✅ | | [`SPA`](@ref SolarPosition.Positioning.SPA) | [Reda & Andreas, 2004](https://doi.org/10.1016/j.solener.2003.12.003) | ±0.0003° | Built-in | ✅ | +| [`Iqbal`](@ref SolarPosition.Positioning.Iqbal) | [Iqbal, 1983](https://doi.org/10.1016/B978-0-12-373750-2.X5001-0) | ±0.0100° | None | ✅ | +| [`Michalsky`](@ref SolarPosition.Positioning.Michalsky) | [Michalsky, 1988](https://doi.org/10.1016/0038-092X(88)90045-X) | ±0.0100° | [`MICHALSKY`](@ref SolarPosition.Refraction.MICHALSKY) | ✅ | +| [`SG2`](@ref SolarPosition.Positioning.SG2) | [Blanc & Wald, 2012](https://doi.org/10.1016/j.solener.2012.07.018) | ±0.0030° | [`SG2Refraction`](@ref SolarPosition.Refraction.SG2Refraction) | ✅ | ## Refraction correction algorithms @@ -96,7 +99,7 @@ Atmospheric refraction correction algorithms available in SolarPosition.jl. | [`ARCHER`](@ref SolarPosition.Refraction.ARCHER) | Archer et al., 1980 | None | ✅ | | [`BENNETT`](@ref SolarPosition.Refraction.BENNETT) | [Bennett, 1982](https://doi.org/10.1017/S0373463300022037) | Pressure, Temperature | ✅ | | [`MICHALSKY`](@ref SolarPosition.Refraction.MICHALSKY) | [Michalsky, 1988](https://doi.org/10.1016/0038-092X(88)90045-X) | None | ✅ | -| [`SG2`](@ref SolarPosition.Refraction.SG2) | [Blanc & Wald, 2012](https://doi.org/10.1016/j.solener.2012.07.018) | Pressure, Temperature | ✅ | +| [`SG2Refraction`](@ref SolarPosition.Refraction.SG2Refraction) | [Blanc & Wald, 2012](https://doi.org/10.1016/j.solener.2012.07.018) | Pressure, Temperature | ✅ | | [`SPARefraction`](@ref SolarPosition.Refraction.SPARefraction) | [Reda & Andreas, 2004](https://doi.org/10.1016/j.solener.2003.12.003) | Pressure, Temperature | ✅ | ## Extensions diff --git a/docs/src/positioning.md b/docs/src/positioning.md index e72fc49..e7ad391 100644 --- a/docs/src/positioning.md +++ b/docs/src/positioning.md @@ -63,6 +63,9 @@ The following solar position algorithms are currently implemented in SolarPositi | [`Walraven`](@ref SolarPosition.Positioning.Walraven) | [Wal78](@cite) | ±0.0100° | None | ✅ | | [`USNO`](@ref SolarPosition.Positioning.USNO) | [USNO](@cite) | ±0.0500° | None | ✅ | | [`SPA`](@ref SolarPosition.Positioning.SPA) | [RA04](@cite) | ±0.0003° | Built-in | ✅ | +| [`Iqbal`](@ref SolarPosition.Positioning.Iqbal) | [Iqb83](@cite) | ±0.0100° | None | ✅ | +| [`Michalsky`](@ref SolarPosition.Positioning.Michalsky) | [Mic88](@cite) | ±0.0100° | [`MICHALSKY`](@ref SolarPosition.Refraction.MICHALSKY) | ✅ | +| [`SG2`](@ref SolarPosition.Positioning.SG2) | [BW12](@cite) | ±0.0030° | [`SG2Refraction`](@ref SolarPosition.Refraction.SG2Refraction) | ✅ | ## [PSA](@id psa-algorithm) @@ -121,3 +124,40 @@ position calculation with periodic terms for Earth heliocentric longitude and la ```@docs SolarPosition.Positioning.SPA ``` + +## [Iqbal](@id iqbal-algorithm) + +The Iqbal algorithm is a fast, low-complexity method that derives the solar declination +and equation of time from a truncated Fourier series in the day angle. + +The formulation was compiled by [Iqb83](@cite) and builds on the Fourier-series +representation of [Spe71](@cite). + +```@docs +SolarPosition.Positioning.Iqbal +``` + +## [Michalsky](@id michalsky-algorithm) + +The Michalsky algorithm implements the Astronomical Almanac's approximate solar position +algorithm, with a stated accuracy of ±0.01° between 1950 and 2050. It exposes options for +the azimuth-quadrant correction and the Julian date formulation. + +The algorithm was published by [Mic88](@cite); the azimuth correction that makes it valid +for all latitudes is from [Spe89](@cite). + +```@docs +SolarPosition.Positioning.Michalsky +``` + +## [SG2](@id sg2-algorithm) + +The SG2 (Second Generation) algorithm is optimized for fast and accurate computation over +multi-decadal periods, with a stated accuracy of ±0.003° between 1980 and 2030. Dates +outside this range raise an `ArgumentError`. + +The algorithm was published by [BW12](@cite). + +```@docs +SolarPosition.Positioning.SG2 +``` diff --git a/docs/src/refraction.md b/docs/src/refraction.md index 9503f8d..5fb82f4 100644 --- a/docs/src/refraction.md +++ b/docs/src/refraction.md @@ -29,7 +29,7 @@ of the available algorithms: | [`ARCHER`](@ref SolarPosition.Refraction.ARCHER) | [Arc80](@cite) | None | ✅ | | [`BENNETT`](@ref SolarPosition.Refraction.BENNETT) | [Ben82](@cite) | Pressure, Temperature | ✅ | | [`MICHALSKY`](@ref SolarPosition.Refraction.MICHALSKY) | [Mic88](@cite) | None | ✅ | -| [`SG2`](@ref SolarPosition.Refraction.SG2) | [BW12](@cite) | Pressure, Temperature | ✅ | +| [`SG2Refraction`](@ref SolarPosition.Refraction.SG2Refraction) | [BW12](@cite) | Pressure, Temperature | ✅ | | [`SPARefraction`](@ref SolarPosition.Refraction.SPARefraction) | [RA04](@cite) | Pressure, Temperature | ✅ | To calculate refraction, we can use the [`refraction`](@ref SolarPosition.Refraction.refraction) function: @@ -83,7 +83,7 @@ using CairoMakie # Define models and elevation range models = [("Archer", SolarPosition.Refraction.ARCHER()), ("Bennett", SolarPosition.Refraction.BENNETT()), ("Hughes", SolarPosition.Refraction.HUGHES()), ("Michalsky", SolarPosition.Refraction.MICHALSKY()), - ("SG2", SolarPosition.Refraction.SG2()), ("SPA", SolarPosition.Refraction.SPARefraction())] + ("SG2", SolarPosition.Refraction.SG2Refraction()), ("SPA", SolarPosition.Refraction.SPARefraction())] elevation = -1.5:0.1:90.0 # Create figure with two subplots @@ -161,7 +161,7 @@ solar position calculations. It includes special handling for very low elevation SolarPosition.Refraction.MICHALSKY ``` -## [SG2](@id sg2-refraction) +## [SG2Refraction](@id sg2-refraction) The SG2 (Second Generation) refraction algorithm is optimized for fast computation over multi-decadal periods. @@ -169,7 +169,7 @@ Developed by [BW12](@cite), this algorithm uses a two-regime approach with diffe for elevations above and below a threshold. It accounts for atmospheric pressure and temperature. ```@docs -SolarPosition.Refraction.SG2 +SolarPosition.Refraction.SG2Refraction ``` ## [SPARefraction](@id spa-refraction) diff --git a/docs/src/refs.bib b/docs/src/refs.bib index 4e35be8..34ff1fd 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -127,6 +127,25 @@ @article{Spe89 url = {https://doi.org/10.1016/0038-092X(89)90039-X} } +@book{Iqb83, + author = {Iqbal, M.}, + title = {An Introduction to Solar Radiation}, + publisher = {Academic Press}, + year = {1983}, + doi = {10.1016/B978-0-12-373750-2.X5001-0}, + url = {https://doi.org/10.1016/B978-0-12-373750-2.X5001-0} +} + +@article{Spe71, + author = {Spencer, J. W.}, + title = {Fourier series representation of the position of the Sun}, + journal = {Search}, + volume = {2}, + number = {5}, + pages = {172}, + year = {1971} +} + @article{Ben82, author = {Bennett, G. G.}, title = {The Calculation of Astronomical Refraction in Marine Navigation}, diff --git a/src/Positioning/Positioning.jl b/src/Positioning/Positioning.jl index 1a5f78b..6d31343 100644 --- a/src/Positioning/Positioning.jl +++ b/src/Positioning/Positioning.jl @@ -452,9 +452,13 @@ include("noaa.jl") include("walraven.jl") include("usno.jl") include("spa.jl") +include("iqbal.jl") +include("michalsky.jl") +include("sg2.jl") export Observer, - PSA, NOAA, Walraven, USNO, SPA, solar_position, solar_position!, SolPos, ApparentSolPos + PSA, NOAA, Walraven, USNO, SPA, Iqbal, Michalsky, SG2, + solar_position, solar_position!, SolPos, ApparentSolPos export SolarAlgorithm, AbstractSolPos, AbstractApparentSolPos export calculate_deltat diff --git a/src/Positioning/iqbal.jl b/src/Positioning/iqbal.jl new file mode 100644 index 0000000..1c8f0ef --- /dev/null +++ b/src/Positioning/iqbal.jl @@ -0,0 +1,72 @@ +""" + $(TYPEDEF) + +Iqbal solar position algorithm. + +A lightweight algorithm that obtains the solar declination and equation of time from a +truncated Fourier series in the day angle, then derives the zenith and azimuth from the +standard spherical-trigonometry relations. No atmospheric refraction is applied by default. + +# Accuracy +The truncated Fourier expansion gives a declination accurate to about ±0.01°, which makes +this algorithm a good choice when speed matters more than sub-arcminute precision. + +# Literature +Based on the formulation compiled by [Iqb83](@cite), which builds on the Fourier-series +representation of [Spe71](@cite). + +# Fields +$(TYPEDFIELDS) +""" +struct Iqbal <: SolarAlgorithm end + +function _solar_position(obs::Observer{T}, dt::DateTime, ::Iqbal) where {T} + # day angle [radians] + day_angle = 2π * (dayofyear(dt) - 1) / 365 + + (sin_Γ, cos_Γ) = sincos(day_angle) + (sin_2Γ, cos_2Γ) = sincos(2 * day_angle) + (sin_3Γ, cos_3Γ) = sincos(3 * day_angle) + + # solar declination [degrees] + declination = rad2deg( + 0.006918 - 0.399912 * cos_Γ + 0.070257 * sin_Γ - + 0.006758 * cos_2Γ + 0.000907 * sin_2Γ - + 0.002697 * cos_3Γ + 0.00148 * sin_3Γ, + ) + + # equation of time [minutes] + eot = + ( + 0.0000075 + 0.001868 * cos_Γ - 0.032077 * sin_Γ - + 0.014615 * cos_2Γ - 0.040849 * sin_2Γ + ) * 1440 / (2π) + + # hour angle [degrees] + hour_angle = (fractional_hour(dt) - 12) * 15 + obs.longitude + eot / 4 + + (sin_δ, cos_δ) = sincosd(declination) + (sin_ω, cos_ω) = sincosd(hour_angle) + + # zenith angle [degrees] + zenith = acosd(sin_δ * obs.sin_lat + cos_δ * obs.cos_lat * cos_ω) + + # azimuth [degrees], measured eastward from north. The arctan form selects the correct + # quadrant where the closed-form expression would be ambiguous. + azimuth = mod( + atand( + sin_ω * cos_δ, + cos_ω * obs.sin_lat * cos_δ - obs.cos_lat * sin_δ, + ) + 180, + 360, + ) + + return SolPos{T}(azimuth, 90 - zenith, zenith) +end + +function _solar_position(obs, dt, alg::Iqbal, ::DefaultRefraction) + return _solar_position(obs, dt, alg, NoRefraction()) +end + +# Iqbal with DefaultRefraction returns SolPos (no refraction by default) +result_type(::Type{Iqbal}, ::Type{DefaultRefraction}, ::Type{T}) where {T} = SolPos{T} diff --git a/src/Positioning/michalsky.jl b/src/Positioning/michalsky.jl new file mode 100644 index 0000000..7ad5589 --- /dev/null +++ b/src/Positioning/michalsky.jl @@ -0,0 +1,117 @@ +using ..Refraction: MICHALSKY + +""" + $(TYPEDEF) + +Michalsky solar position algorithm. + +Implements the Astronomical Almanac's approximate solar position algorithm. By default the +[`MICHALSKY`](@ref) atmospheric refraction model is applied, matching the original +publication. + +# Accuracy +Claimed accuracy: ±0.01° between 1950 and 2050 when using the original Julian date +formulation. + +# Options +- `spencer_correction`: when `true` (default), applies the azimuth-quadrant correction of + [Spe89](@cite) so the algorithm is valid at all latitudes. The original formulation + (`false`) is only valid in the northern hemisphere. +- `julian_date`: `:original` (default) uses the integer-based Julian date from the original + paper, which guarantees the stated accuracy only between 1950 and 2050; `:standard` uses + the exact Julian date and remains usable outside that window. + +# Literature +Based on the algorithm of [Mic88](@cite) with the azimuth correction of [Spe89](@cite). + +# Fields +$(TYPEDFIELDS) +""" +struct Michalsky <: SolarAlgorithm + "Apply the Spencer (1989) azimuth-quadrant correction (valid for all latitudes)." + spencer_correction::Bool + "Julian date formulation: `:original` or `:standard`." + julian_date::Symbol +end + +Michalsky(; spencer_correction::Bool = true, julian_date::Symbol = :original) = + Michalsky(spencer_correction, julian_date) + +function _solar_position(obs::Observer{T}, dt::DateTime, alg::Michalsky) where {T} + hour = fractional_hour(dt) + + # Julian date + jd = if alg.julian_date === :original + delta = year(dt) - 1949 + leap = floor(delta / 4) + 2432916.5 + delta * 365 + leap + dayofyear(dt) + hour / 24 + elseif alg.julian_date === :standard + datetime2julian(dt) + else + throw( + ArgumentError( + "`julian_date` must be :original or :standard, got :$(alg.julian_date)", + ), + ) + end + + # days since J2000.0 + n = jd - 2451545.0 + + # mean longitude, mean anomaly and ecliptic longitude [degrees] + L = mod(280.46 + 0.9856474 * n, 360) + g = mod(357.528 + 0.9856003 * n, 360) + l = mod(L + 1.915 * sind(g) + 0.02 * sind(2 * g), 360) + + # obliquity of the ecliptic [degrees] + ep = 23.439 - 0.0000004 * n + + # declination and right ascension [degrees] + dec = asind(sind(ep) * sind(l)) + ra = mod(rad2deg(atan(cosd(ep) * sind(l), cosd(l))), 360) + + # Greenwich and local mean sidereal time [hours] + gmst = mod(6.697375 + 0.0657098242 * n + hour, 24) + lmst = mod(gmst + obs.longitude / 15, 24) + + # hour angle [hours], wrapped to [-12, 12) + ha = mod(lmst - ra / 15 + 12, 24) - 12 + + # elevation and azimuth [degrees] + el = asind(sind(dec) * obs.sin_lat + cosd(dec) * obs.cos_lat * cosd(15 * ha)) + az = asind(-cosd(dec) * sind(15 * ha) / cosd(el)) + + if alg.spencer_correction + # Spencer (1989) quadrant correction + cos_az = sind(dec) - sind(el) * obs.sin_lat + if cos_az >= 0 && sind(az) < 0 + az += 360 + end + if cos_az < 0 + az = 180 - az + end + else + # Original Michalsky quadrant assignment via the critical elevation. The ratio can + # exceed the asin domain (e.g. at the equator); treat that as undefined so neither + # correction applies, matching the reference behaviour. + ratio = sind(dec) / obs.sin_lat + elc = abs(ratio) <= 1 ? asind(ratio) : T(NaN) + if el >= elc + az = 180 - az + end + if el <= elc && ha > 0 + az += 360 + end + az = mod(az, 360) + end + + return SolPos{T}(az, el, 90 - el) +end + +function _solar_position(obs, dt, alg::Michalsky, ::DefaultRefraction) + return _solar_position(obs, dt, alg, MICHALSKY()) +end + +# Michalsky with DefaultRefraction returns ApparentSolPos (uses MICHALSKY refraction) +result_type(::Type{Michalsky}, ::Type{DefaultRefraction}, ::Type{T}) where {T} = + ApparentSolPos{T} diff --git a/src/Positioning/sg2.jl b/src/Positioning/sg2.jl new file mode 100644 index 0000000..1a74041 --- /dev/null +++ b/src/Positioning/sg2.jl @@ -0,0 +1,132 @@ +using ..Refraction: SG2Refraction + +# Earth heliocentric longitude periodic terms (frequency [1/day], amplitude [rad], phase [rad]). +const _SG2_HELIO = ( + (1 / 365.261278, 3.401508e-2, 1.60078), + (1 / 182.632412, 3.48644e-4, 1.662976), + (1 / 29.530634, 3.136227e-5, -1.195905), + (1 / 399.52985, 3.578979e-5, -1.042052), + (1 / 291.956812, 2.676185e-5, 2.012613), + (1 / 583.598201, 2.333925e-5, -2.867714), + (1 / 4652.629372, 1.221214e-5, 1.225038), + (1 / 1450.236684, 1.217941e-5, -0.828601), + (1 / 199.459709, 1.343914e-5, -3.108253), + (1 / 365.355291, 8.499475e-4, -2.353709), +) + +""" + $(TYPEDEF) + +SG2 (Second Generation) solar position algorithm. + +A fast, accurate algorithm tuned for multi-decadal time periods. By default the +[`SG2Refraction`](@ref) atmospheric refraction model is applied. + +# Accuracy +Claimed accuracy: ±0.003° between 1980 and 2030. The algorithm is only defined within this +range; calling it for a date outside it throws an `ArgumentError`. + +# Literature +Based on the algorithm of [BW12](@cite). + +# Fields +$(TYPEDFIELDS) +""" +struct SG2 <: SolarAlgorithm end + +function _solar_position(obs::Observer{T}, dt::DateTime, ::SG2) where {T} + yr = year(dt) + mo = month(dt) + dom = day(dt) + hour = fractional_hour(dt) + + # year in decimal form, used both for the validity range and the ΔT branch + year_dec = yr + (mo - 0.5) / 12 + if year_dec < 1980 || year_dec > 2030 + throw( + ArgumentError( + "SG2 is only valid for years 1980–2030, got $(year_dec)", + ), + ) + end + + # ΔT [seconds] from a piecewise polynomial in the (integer) year + (y_ref, c0, c1, c2, c3, c4, c5) = if year_dec <= 1986 + (1975, 45.45, 1.067, -1 / 260, -1 / 718, 0.0, 0.0) + elseif year_dec <= 2005 + (2000, 63.86, 0.3345, -0.060374, 0.0017275, 6.518e-4, 2.374e-5) + else + (2000, 63.48, 0.204, 0.005576, 0.0, 0.0, 0.0) + end + Δy = yr - y_ref + delta_t = c0 + c1 * Δy + c2 * Δy^2 + c3 * Δy^3 + c4 * Δy^4 + c5 * Δy^5 + + # Julian day (UT and TT) + (year_mod, month_mod) = (mo == 1 || mo == 2) ? (yr - 1, mo + 12) : (yr, mo) + jd_ut = + 1721028.0 + dom + floor((153.0 * month_mod - 2.0) / 5.0) + + 365.0 * year_mod + floor(year_mod / 4.0) + hour / 24 - 0.5 - + floor(year_mod / 100.0) + floor(year_mod / 400.0) + jd_tt = jd_ut + delta_t / 86400 + + jd_ut_mod = jd_ut - 2444239.5 + jd_tt_mod = jd_tt - 2444239.5 + + # Earth heliocentric longitude [rad] + sums = 0.0 + for (f_L, ρ_L, φ_L) in _SG2_HELIO + sums += ρ_L * cos(2π * f_L * jd_tt_mod - φ_L) + end + L = mod(sums + (1 / 58.130101) * jd_tt_mod + 1.742145, 2π) + + # geocentric nutation and true obliquity [rad] + D_t = -9.933735e-5 + ξ = 4.263521e-5 + D_ψ = 8.329092e-5 * cos(2π * (1 / 6791.164405) * jd_tt_mod - (-2.052757)) + ε = + 4.456183e-5 * cos(2π * (1 / 6791.164405) * jd_tt_mod - 2.660352) + + (-6.216374e-9) * jd_tt_mod + 4.091383e-1 + + # apparent sun geocentric longitude [rad] + Θ = mod(L + π + D_ψ + D_t, 2π) + + # geocentric declination and right ascension [rad] + decl_g = asin(sin(Θ) * sin(ε)) + ra = atan(sin(Θ) * cos(ε), cos(Θ)) + + # mean sidereal time [rad] + mst = mod(6.300388099 * jd_ut_mod + 1.742079, 2π) + + # observer geocentric coordinates (oblate Earth) + a = 6378140.0 + f = 1 / 298.257282697 + u = atan((1 - f) * tan(obs.latitude_rad)) + x = cos(u) + obs.altitude / a * obs.cos_lat + y = (1 - f) * sin(u) + obs.altitude / a * obs.sin_lat + + # topocentric correction (parallax) + ν = mst + D_ψ * cos(ε) + ω_g = ν + obs.longitude_rad - ra + D_ra = -x * sin(ω_g) / cos(decl_g) * ξ + declination = decl_g + (x * cos(ω_g) * sin(decl_g) - y * cos(decl_g)) * ξ + ω = mst + D_ψ * cos(ε) - ra + obs.longitude_rad - D_ra + + # topocentric azimuth and elevation [rad] + azimuth_rad = + atan( + sin(ω), + cos(ω) * obs.sin_lat - tan(declination) * obs.cos_lat, + ) + π + elevation_rad = + asin(obs.sin_lat * sin(declination) + obs.cos_lat * cos(declination) * cos(ω)) + + elevation = rad2deg(elevation_rad) + return SolPos{T}(mod(rad2deg(azimuth_rad), 360), elevation, 90 - elevation) +end + +function _solar_position(obs, dt, alg::SG2, ::DefaultRefraction) + return _solar_position(obs, dt, alg, SG2Refraction()) +end + +# SG2 with DefaultRefraction returns ApparentSolPos (uses SG2Refraction) +result_type(::Type{SG2}, ::Type{DefaultRefraction}, ::Type{T}) where {T} = ApparentSolPos{T} diff --git a/src/Refraction/Refraction.jl b/src/Refraction/Refraction.jl index 2c0f5fa..121f9de 100644 --- a/src/Refraction/Refraction.jl +++ b/src/Refraction/Refraction.jl @@ -75,7 +75,7 @@ include("sg2.jl") include("spa.jl") export RefractionAlgorithm, NoRefraction, DefaultRefraction -export HUGHES, ARCHER, BENNETT, MICHALSKY, SG2, SPARefraction +export HUGHES, ARCHER, BENNETT, MICHALSKY, SG2Refraction, SPARefraction export refraction end diff --git a/src/Refraction/sg2.jl b/src/Refraction/sg2.jl index 82a3286..a9c17da 100644 --- a/src/Refraction/sg2.jl +++ b/src/Refraction/sg2.jl @@ -1,9 +1,10 @@ """ $(TYPEDEF) -SG2 refraction model. +SG2 (Second Generation) atmospheric refraction model. -Atmospheric refraction correction based on the algorithm in SG2. +Atmospheric refraction correction based on the SG2 algorithm. The matching solar position +algorithm is [`SG2`](@ref). This function calculates the atmospheric refraction correction of the solar elevation angle using the method developed by Ph. Blanc and L. Wald [1]. @@ -12,8 +13,8 @@ elevation angle using the method developed by Ph. Blanc and L. Wald [1]. $(TYPEDFIELDS) # Constructor -- `SG2()`: Uses default parameters: pressure = 101325 Pa, temperature = 12 °C -- `SG2(pressure, temperature)`: Specify custom pressure [Pa] and temperature [°C] +- `SG2Refraction()`: Uses default parameters: pressure = 101325 Pa, temperature = 12 °C +- `SG2Refraction(pressure, temperature)`: Specify custom pressure [Pa] and temperature [°C] # Notes The equation to calculate the refraction correction is given by: @@ -39,10 +40,10 @@ This method was described by [BW12](@cite). using SolarPosition # Create SG2 refraction model with default parameters -sg2 = SG2() +sg2 = SG2Refraction() # Or specify custom atmospheric conditions -sg2_custom = SG2(101325.0, 25.0) # 25°C temperature +sg2_custom = SG2Refraction(101325.0, 25.0) # 25°C temperature # Apply refraction correction to elevation angle elevation = 30.0 # degrees @@ -50,16 +51,16 @@ refraction_correction = refraction(sg2, elevation) apparent_elevation = elevation + refraction_correction ``` """ -struct SG2{T} <: RefractionAlgorithm where {T <: AbstractFloat} +struct SG2Refraction{T} <: RefractionAlgorithm where {T <: AbstractFloat} "Annual average atmospheric pressure [Pascal]" pressure::T "Annual average temperature [°C]" temperature::T end -SG2() = SG2{Float64}(101325.0, 12.0) +SG2Refraction() = SG2Refraction{Float64}(101325.0, 12.0) -function _refraction(model::SG2{T}, elevation_deg::T) where {T <: AbstractFloat} +function _refraction(model::SG2Refraction{T}, elevation_deg::T) where {T <: AbstractFloat} # Convert pressure from Pascal to hPa (hectopascal) pressure_hPa = model.pressure / T(100.0) diff --git a/test/extensions/test-mtk.jl b/test/extensions/test-mtk.jl index cc25eb7..2d52289 100644 --- a/test/extensions/test-mtk.jl +++ b/test/extensions/test-mtk.jl @@ -10,7 +10,7 @@ using SolarPosition: DefaultRefraction, SolarPositionBlock -using SolarPosition: HUGHES, BENNETT, ARCHER, MICHALSKY, SG2 +using SolarPosition: HUGHES, BENNETT, ARCHER, MICHALSKY, SG2Refraction using ModelingToolkit: @named, @variables, @parameters, unknowns, System, mtkcompile using ModelingToolkit: t_nounits as t, D_nounits as D using Dates: DateTime diff --git a/test/positioning/test-interface.jl b/test/positioning/test-interface.jl index 86073a8..a679c87 100644 --- a/test/positioning/test-interface.jl +++ b/test/positioning/test-interface.jl @@ -234,7 +234,7 @@ end @testset "Refraction Integration" begin using SolarPosition.Positioning: ApparentSolPos - using SolarPosition.Refraction: BENNETT, HUGHES, ARCHER, MICHALSKY, SG2 + using SolarPosition.Refraction: BENNETT, HUGHES, ARCHER, MICHALSKY, SG2Refraction obs = Observer(45.0, 10.0, 100.0) dt = DateTime(2023, 6, 21, 12, 0, 0) @@ -250,7 +250,7 @@ end ("HUGHES", HUGHES()), ("ARCHER", ARCHER()), ("MICHALSKY", MICHALSKY()), - ("SG2", SG2()), + ("SG2", SG2Refraction()), ] res = solar_position(obs, dt, alg, refr) @test res isa ApparentSolPos diff --git a/test/positioning/test-iqbal.jl b/test/positioning/test-iqbal.jl index 7553372..9c6c0ab 100644 --- a/test/positioning/test-iqbal.jl +++ b/test/positioning/test-iqbal.jl @@ -28,30 +28,22 @@ function expected_iqbal() return DataFrame(reduce(hcat, values)', columns) end -@testset "Iqbal (not implemented)" begin - @test_skip begin - df_expected = expected_iqbal() - conds = test_conditions() - @test size(df_expected, 1) == 19 - @test size(df_expected, 2) == 3 - @test size(conds, 1) == 19 - @test size(conds, 2) == 4 +@testset "Iqbal" begin + df_expected = expected_iqbal() + conds = test_conditions() + @test size(df_expected, 1) == 19 + @test size(df_expected, 2) == 3 + @test size(conds, 1) == 19 + @test size(conds, 2) == 4 - # TODO: Implement Iqbal algorithm - # struct Iqbal <: SolarAlgorithm end + for ((dt, lat, lon, alt), (exp_elev, exp_zen, exp_az)) in + zip(eachrow(conds), eachrow(df_expected)) + obs = ismissing(alt) ? Observer(lat, lon) : Observer(lat, lon, altitude = alt) - # for ((dt, lat, lon, alt), (exp_elev, exp_zen, exp_az)) in - # zip(eachrow(conds), eachrow(df_expected)) - # if ismissing(alt) - # obs = Observer(lat, lon) - # else - # obs = Observer(lat, lon, altitude = alt) - # end - # - # res = solar_position(obs, dt, Iqbal()) - # @test isapprox(res.elevation, exp_elev, atol = 1e-8) - # @test isapprox(res.zenith, exp_zen, atol = 1e-8) - # @test isapprox(res.azimuth, exp_az, atol = 1e-8) - # end + res = solar_position(obs, dt, Iqbal()) + @test res isa SolPos + @test isapprox(res.elevation, exp_elev, atol = 1.0e-6) + @test isapprox(res.zenith, exp_zen, atol = 1.0e-6) + @test isapprox(res.azimuth, exp_az, atol = 1.0e-6) end end diff --git a/test/positioning/test-michalsky.jl b/test/positioning/test-michalsky.jl index f9aa75b..c595549 100644 --- a/test/positioning/test-michalsky.jl +++ b/test/positioning/test-michalsky.jl @@ -84,35 +84,49 @@ function expected_michalsky_spencer_false() return DataFrame(reduce(hcat, values)', columns) end -@testset "Michalsky (not implemented)" begin - @test_skip begin - # TODO: Implement Michalsky algorithm with options: - # - spencer_correction: true/false - # - julian_date: "original" or "pandas" - # struct Michalsky <: SolarAlgorithm - # spencer_correction::Bool - # julian_date::String - # end - - @testset "Original Julian Date" begin - df_expected = expected_michalsky_original_julian() - conds = test_conditions() - - # Test implementation would go here - end +# Compare every condition against the expected actual and apparent angles. The default +# refraction for Michalsky is the MICHALSKY model, so the result carries apparent angles. +function check_michalsky(alg, df_expected) + conds = test_conditions() + @test size(df_expected, 1) == 19 + @test size(df_expected, 2) == 5 + for ((dt, lat, lon, alt), row) in zip(eachrow(conds), eachrow(df_expected)) + obs = ismissing(alt) ? Observer(lat, lon) : Observer(lat, lon, altitude = alt) + res = solar_position(obs, dt, alg) + @test res isa ApparentSolPos + @test isapprox(res.elevation, row.elevation, atol = 1.0e-6) + @test isapprox(res.apparent_elevation, row.apparent_elevation, atol = 1.0e-6) + @test isapprox(res.zenith, row.zenith, atol = 1.0e-6) + @test isapprox(res.apparent_zenith, row.apparent_zenith, atol = 1.0e-6) + @test isapprox(res.azimuth, row.azimuth, atol = 1.0e-6) + end + return +end - @testset "Pandas Julian Date" begin - df_expected = expected_michalsky_pandas_julian() - conds = test_conditions() +@testset "Michalsky" begin + @testset "Original Julian Date" begin + check_michalsky(Michalsky(), expected_michalsky_original_julian()) + end - # Test implementation would go here - end + @testset "Standard Julian Date" begin + check_michalsky( + Michalsky(julian_date = :standard), + expected_michalsky_pandas_julian(), + ) + end - @testset "Spencer Correction False" begin - df_expected = expected_michalsky_spencer_false() - conds = test_conditions() + @testset "Spencer Correction False" begin + check_michalsky( + Michalsky(spencer_correction = false), + expected_michalsky_spencer_false(), + ) + end - # Test implementation would go here - end + @testset "Invalid Julian date option" begin + @test_throws ArgumentError solar_position( + Observer(45.0, 10.0), + DateTime(2020, 10, 17, 12, 30), + Michalsky(julian_date = :bogus), + ) end end diff --git a/test/positioning/test-psa.jl b/test/positioning/test-psa.jl index 3263b15..6b962c8 100644 --- a/test/positioning/test-psa.jl +++ b/test/positioning/test-psa.jl @@ -47,7 +47,7 @@ end # test with other refraction algorithms - algs = [ARCHER(), MICHALSKY(), SG2()] + algs = [ARCHER(), MICHALSKY(), SG2Refraction()] for alg in algs res_alg = solar_position(obs, dt, PSA(), alg) @test res_alg isa ApparentSolPos diff --git a/test/positioning/test-sg2.jl b/test/positioning/test-sg2.jl index 663d4bf..715c907 100644 --- a/test/positioning/test-sg2.jl +++ b/test/positioning/test-sg2.jl @@ -30,73 +30,51 @@ function expected_sg2() return DataFrame(reduce(hcat, values)', columns) end -@testset "SG2 (not implemented)" begin - @test_skip begin - df_expected = expected_sg2() - conds = test_conditions() - @test size(df_expected, 1) == 19 - @test size(df_expected, 2) == 5 - @test size(conds, 1) == 19 - @test size(conds, 2) == 4 +@testset "SG2" begin + df_expected = expected_sg2() + conds = test_conditions() + @test size(df_expected, 1) == 19 + @test size(df_expected, 2) == 5 + @test size(conds, 1) == 19 + @test size(conds, 2) == 4 - # TODO: Implement SG2 algorithm - # Note: SG2 is valid only between 1980 and 2030 - # struct SG2 <: SolarAlgorithm end + # The reference values were generated with an independent SG2 implementation, so the + # tolerance matches the algorithm's stated accuracy of 0.003°. + for ((dt, lat, lon, alt), row) in zip(eachrow(conds), eachrow(df_expected)) + obs = ismissing(alt) ? Observer(lat, lon) : Observer(lat, lon, altitude = alt) - # for ((dt, lat, lon, alt), row) in zip(eachrow(conds), eachrow(df_expected)) - # year_val = year(dt) - # - # # Skip years outside valid range (1980-2030) - # if year_val < 1980 || year_val > 2030 - # # These should return NaN or throw error - # @test all(isnan.([row.elevation, row.zenith, row.azimuth])) - # continue - # end - # - # if ismissing(alt) - # obs = Observer(lat, lon) - # else - # obs = Observer(lat, lon, altitude = alt) - # end - # - # # SG2 includes built-in refraction correction - # res = solar_position(obs, dt, SG2()) - # @test isapprox(res.elevation, row.elevation, atol = 1e-3) - # @test isapprox(res.zenith, row.zenith, atol = 1e-3) - # @test isapprox(res.azimuth, row.azimuth, atol = 1e-3) - # @test isapprox(res.apparent_elevation, row.apparent_elevation, atol = 1e-3) - # @test isapprox(res.apparent_zenith, row.apparent_zenith, atol = 1e-3) - # end + # SG2 is only defined for 1980–2030; out-of-range dates throw instead of returning. + if isnan(row.elevation) + @test_throws ArgumentError solar_position(obs, dt, SG2()) + continue + end + + res = solar_position(obs, dt, SG2()) + @test res isa ApparentSolPos + @test isapprox(res.elevation, row.elevation, atol = 3.0e-3) + @test isapprox(res.apparent_elevation, row.apparent_elevation, atol = 3.0e-3) + @test isapprox(res.zenith, row.zenith, atol = 3.0e-3) + @test isapprox(res.apparent_zenith, row.apparent_zenith, atol = 3.0e-3) + @test isapprox(res.azimuth, row.azimuth, atol = 3.0e-3) end @testset "SG2 year validation" begin - @test_skip begin - # Test that SG2 throws error for years outside 1980-2030 range - obs = Observer(50.0, 10.0) - - # Should throw error for 1960 - dt_1960 = ZonedDateTime(1960, 1, 1, 12, 0, 0, tz"UTC") - @test_throws ArgumentError solar_position(obs, dt_1960, SG2()) - - # Should throw error for 2035 - dt_2035 = ZonedDateTime(2035, 1, 1, 12, 0, 0, tz"UTC") - @test_throws ArgumentError solar_position(obs, dt_2035, SG2()) - end + obs = Observer(50.0, 10.0) + @test_throws ArgumentError solar_position( + obs, ZonedDateTime(1960, 1, 1, 12, 0, 0, tz"UTC"), SG2(), + ) + @test_throws ArgumentError solar_position( + obs, ZonedDateTime(2035, 1, 1, 12, 0, 0, tz"UTC"), SG2(), + ) end @testset "SG2 site elevation" begin - @test_skip begin - # Make sure that site elevations are properly used by the SG2 algorithm - times = [ZonedDateTime(2020, 1, 1, 12, 0, 0, tz"UTC+2")] - latitude, longitude = 50.0, 10.0 - - obs_zero = Observer(latitude, longitude, altitude = 0.0) - obs_high = Observer(latitude, longitude, altitude = 4000.0) + # Make sure that site elevations are properly used by the SG2 algorithm + t = ZonedDateTime(2020, 1, 1, 12, 0, 0, tz"UTC+2") + obs_zero = Observer(50.0, 10.0, altitude = 0.0) + obs_high = Observer(50.0, 10.0, altitude = 4000.0) - zero_elevation = solar_position(obs_zero, times[1], SG2()) - high_elevation = solar_position(obs_high, times[1], SG2()) - - @test zero_elevation.elevation ≠ high_elevation.elevation - end + @test solar_position(obs_zero, t, SG2()).elevation ≠ + solar_position(obs_high, t, SG2()).elevation end end diff --git a/test/refraction/test-refraction.jl b/test/refraction/test-refraction.jl index c5a53ae..21a6ab8 100644 --- a/test/refraction/test-refraction.jl +++ b/test/refraction/test-refraction.jl @@ -1,7 +1,7 @@ """Unit tests for atmospheric refraction algorithms""" using SolarPosition.Refraction: - HUGHES, ARCHER, BENNETT, MICHALSKY, SG2, SPARefraction, refraction + HUGHES, ARCHER, BENNETT, MICHALSKY, SG2Refraction, SPARefraction, refraction # test elevation angles in degrees test_elevation_angles() = [-1.0, -0.6, 0.0, 1.0, 4.0, 6.0, 10.0, 90.0] @@ -76,7 +76,7 @@ test_algorithms = [ ("Archer", () -> ARCHER(), expected["Archer"]), ("Bennett", () -> BENNETT(101325.0, 12.0), expected["Bennett"]), ("Michalsky", () -> MICHALSKY(), expected["Michalsky"]), - ("SG2", () -> SG2(101325.0, 12.0), expected["SG2"]), + ("SG2", () -> SG2Refraction(101325.0, 12.0), expected["SG2"]), ("SPARefraction", () -> SPARefraction(101325.0, 12.0), expected["SPA"]), ]