Add the Michalsky positioning algorithm - #120
Merged
Merged
Conversation
Implements the Astronomical Almanac series as presented by Michalsky (1988), with the Spencer (1989) azimuth quadrant correction on by default so the result is valid in both hemispheres. The default refraction model is MICHALSKY, matching the paper, so a default call returns ApparentSolPos. The paper's Julian date uses a naive leap year count from 1949 and only holds the stated 0.01 degree accuracy between 1950 and 2050. Both variants are available through julian_date, :original and :standard, and they agree exactly inside that window. Rather than build the ~2.43e6 Julian Date and subtract, which would materialise that magnitude at precision T, the constant J2000 offset is folded into the integer day count. The reference tables were stale and are regenerated -------------------------------------------------- The tables in test-michalsky.jl were generated against a test_conditions() whose base time was 12:30 UTC. That time is 1800 s from noon, which is not a multiple of the 84.375 s grid the reference values rely on for an exactly representable Julian Date, and the conditions were later moved to 13:30, which is. The two tables in this file and in test-sg2.jl were never regenerated, and because both testsets were skipped nobody noticed. Anyone implementing against them would have seen roughly 15 degrees of error on 14 of the 19 rows. Solving for the instant each stored row corresponded to put all 19 rows on the old times to within 4e-9, which is what confirmed the implementation before the tables were replaced. The three tables are now exact solposx output at the current conditions, which lets the tolerance go from 1e-6 to 1e-10, matching the other algorithms. benchmark/CondaPkg.toml is added so that generation is reproducible; it mirrors the tracked docs/CondaPkg.toml exactly. test-interface.jl kept its own copy of the algorithm to result type mapping, which sent Michalsky down the SolPos branch and made the in-place buffer fail to convert. That copy is replaced by a call to result_type, so the next algorithm with a default refraction model does not hit it. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01KhdBYRmHazfXZHWCRWt5fT
Contributor
Benchmark Results (Julia vlts)Time benchmarks
Memory benchmarks
|
Contributor
Benchmark Results (Julia v1)Time benchmarks
Memory benchmarks
|
Contributor
Benchmark Results (Julia vpre)Time benchmarks
Memory benchmarks
|
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #120 +/- ##
=========================================
Coverage 100.00% 100.00%
=========================================
Files 20 21 +1
Lines 766 811 +45
=========================================
+ Hits 766 811 +45 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Codecov flagged one uncovered line, the azimuth += 360 arm of _original_quadrant. That branch needs the sun below the critical elevation after solar noon, and every one of the 19 shared test conditions is either daytime or before noon, so it was unreachable from test_conditions(). A branch tally over those conditions confirmed it: 15 hits on the 180 - azimuth arm, 4 on neither, 1 on the undefined ratio at the equator, and 0 here. test_conditions() is shared with every other algorithm's reference table, so it cannot grow a row. Two evening cases at latitude plus and minus 45 cover the branch instead, and their azimuths above 180 degrees are the full turn taking effect. The timestamp is 29025 s from noon, which is 344 times the 84.375 s grid and a whole number of seconds, so the Julian Date stays exactly representable. Values are solposx output, the same provenance as the tables. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01KhdBYRmHazfXZHWCRWt5fT
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Implements the Astronomical Almanac series as presented by Michalsky (1988), closing the
@test_skipthat has been sitting intest/positioning/test-michalsky.jl.spencer_correction = falsereproduces the original northern hemisphere formulation.MICHALSKY, matching the paper, so a default call returnsApparentSolPos.julian_dateselects the paper's integer based Julian date,:original, or the exact one,:standard. The paper's naive leap year count from 1949 only holds the stated 0.01° accuracy between 1950 and 2050; the two agree exactly inside that window and diverge outside it.Rather than build the ~2.43e6 Julian Date and subtract J2000, which would materialise that magnitude at precision
Tand break the magnitude-safety rule inCLAUDE.md, the constant offset is folded into the integer day count. Every literal is converted atTand everyasinargument goes throughunit_clamp.The stored reference tables were stale
This is the part worth a close look, and probably why this algorithm stalled before.
The tables in
test-michalsky.jlwere generated against atest_conditions()whose base time was 12:30 UTC. That is 1800 s from noon, which is not a multiple of the 84.375 s grid the reference values depend on for an exactly representable Float64 Julian Date. The conditions were later moved to 13:30 UTC, which is exactly 64 × 84.375. The tables here and intest-sg2.jlwere never regenerated, and because both testsets were skipped, nobody noticed. Implementing against them shows roughly 15° of error on 14 of the 19 rows, which looks exactly like a broken algorithm.I verified this rather than assuming it: solving for the instant each stored row actually corresponds to put all 19 rows within 4e-9 of this implementation at the old times, including latitude ±90°, longitude ±180°, and the years 1800 and 2200. That validated the implementation against 285 independent numbers before any table was touched.
test_conditions()12:30+02:00)12:30-02:00)All three tables are now exact
solposxoutput at the current conditions, generated through the existing CondaPkg setup. This implementation matchessolposxbit for bit on row 1 (26.45700256925404), which let the tolerance tighten from 1e-6 to 1e-10, the same as every other algorithm.benchmark/CondaPkg.tomlis added so that generation is reproducible; it mirrors the already trackeddocs/CondaPkg.tomlexactly.test-sg2.jlhas the identical staleness and is left alone here, since SG2 positioning is still unimplemented and its name would collide with theSG2refraction model.A duplicated mapping in the interface tests
test-interface.jlkept its own copy of the algorithm to result type mapping, hardcoded asSPA/NOAA→ApparentSolPosand everything else →SolPos. Michalsky fell into the wrong branch, and the in-place test built aSolPosbuffer that then failed withCannot convert ApparentSolPos{Float64} to SolPos{Float64}. Instead of adding Michalsky to the list, the copy is replaced by a call toresult_type, with the buffer's fields taken fromfieldnames, so the next algorithm with a default refraction model will not hit this.expected-values.jlwas also missing theMichalskyimport;test_algorithms()had been resolving it by luck through another test file's import.Verification
Full
Pkg.test()passes, exit 0, no failures or errors. Michalsky 350/350 where it was previously a single skip, Interface 755/755, Autodiff 111/111, Typestability 274/274, Aqua 11/11, JET 1/1. Adding Michalsky totest_algorithms()means the interface, precision, type-stability and autodiff testsets all cover it.docs/make.jlexits 0 and the fullprek run -ahook set passes.One naming question
solposxspells the optionjulian_date='original'|'pandas'. I used:originaland:standard, since "pandas" names a Python library that has nothing to do with the Julian date in a Julia package; the docstring records the correspondence. Happy to switch to:pandasif you would rather match upstream exactly.🤖 Generated with Claude Code