Skip to content

Model the sun's gravitational light deflection (C library) - #31

Open
arahlin wants to merge 1 commit into
faster-sidereal-timefrom
light-deflection
Open

arahlin wants to merge 1 commit into
faster-sidereal-timefrom
light-deflection

Conversation

@arahlin

@arahlin arahlin commented Sep 27, 2026 •

Copy link
Copy Markdown
Owner

Replaces #29, which GitHub marked merged when light-deflection was briefly rebased
below its own base.

Based on faster-sidereal-time (#30). C library and qpoint only.

The sun bends the incoming ray, and qpoint did not account for it. It was the whole of
the 14 mas by which azel2radec differed from astropy, and of the 43 mas coming back:

against astropy before after
azel2radec 14.30 mas 0.79 mas
radec2azel 42.56 mas 0.84 mas

test_astropy.py's tolerance drops from 0.1 arcsec to 5 mas, and make_qpoint there
switches the term on beside the IERS rates, astropy applying it too.

How it is wired

rate_defl turns it on and defaults to never, like rate_dut1 and rate_wobble:
~20 mas for about a tenth of azel2bore, so a caller who does not ask gets byte-for-byte
what they got before. Only the sun's position is cached at that rate — the deflection
depends on where the telescope points, so the per-sample part runs every sample, as
aberration does. It sits last in the forward chain and first in the inverse, physics
applying it before aberration.

Aberration and deflection are one function, qp_apply_aaber_defl, replacing
qp_apply_annual_aberration. eraEpv00 costs 7.8 µs a call and serves both — aberration
wants the earth's barycentric velocity, deflection its heliocentric position, and it
returns both. Calling it twice doubled the whole transform with both rates on always
(3202 ms against 1610 for 200k samples). Sharing the pointing vector and summing the two
rotation vectors into one quaternion halved the per-sample cost again, 35.5 → 17.5 ns.

Switched on it costs 1.10x on azel2bore. Switched off the output is bit-identical to the
pre-deflection build over 240k values, including the mean_aber=False and fast_aber=False
paths.

Three implementation choices:

  • eraLdsun does the arithmetic, not the textbook 4.07 mas · cot(elongation/2):
    ERFA is already vendored, and its routine takes the sun-to-observer vector directly
    rather than needing the elongation.
  • The rotation always takes the small-angle form, unlike aberration's fast_aber
    switch: the deflection never exceeds ~0.1″, so the dropped cubic term is ~2e-20 rad and
    the exact path measured nearly twice the cost. fast_aber=False keeps its own quaternion.
  • The inverse is a sign flip. eraLdsun has no closed-form inverse — ERFA iterates in
    eraAticq — but negating is second order, 1e-9 mas.

Also

eraEpv00 was being called with its two outputs aliased to one array, so the heliocentric
position was overwritten by the barycentric one. Harmless while only the velocity was
wanted; deflection needs both.

qp_reset_rates enumerates the rates one at a time and was missing qp_reset_rate_defl,
so reset_rates carried the previous chunk's sun position into the next. TestResetRates
is rewritten to pin the property rather than call it and assert nothing: with a rate set to
'once' the correction freezes at the first sample, so after reset_rates the answer must
match a freshly built QPoint — parametrised over every rate the package reports, so a
rate added later is covered without anyone remembering.

C API: qp_memory_t gains state_defl, state_defl_inv, e_sun and em_sun; the
ctypes mirror tracks them and all 55 field offsets were checked against the compiled struct.

Tests

465 passing, 14 skipped. TestLightDeflection asserts the term is present, that switching
it off reproduces the old 14 mas and 43 mas exactly, and pins size and direction against
erfa.ldsun.

The module's shared TOL_ARCSEC moves with the term rather than being left loose, since
every comparison in the file now has the deflection switched on. make_qpoint lets a
caller override one of the rates it sets, which is how TestLightDeflection gets the same
QPoint with the term off to show what it is worth.

Coverage, from a temporary -O0 --coverage build under gcc: qp_apply_aaber_defl at 95%
of its 41 lines, qpoint.c unmoved at 95%. Worth stating because rate_defl defaults to
never, so a term nobody switched on would read as zero coverage and look identical to one
that does not work.

🤖 Generated with Claude Code

@arahlin
arahlin force-pushed the faster-sidereal-time branch from 9a79b6c to 98d0506 Compare September 27, 2026 07:45
@arahlin
arahlin force-pushed the faster-sidereal-time branch from 98d0506 to 0d79e23 Compare September 27, 2026 13:40
@arahlin
arahlin added this pull request to stack #32 September 27, 2026 13:43
@arahlin arahlin self-assigned this Sep 27, 2026
@arahlin
arahlin removed this pull request from stack #32 September 27, 2026 15:43
@arahlin
arahlin force-pushed the faster-sidereal-time branch from 0d79e23 to 0c6af84 Compare September 27, 2026 15:44
@arahlin
arahlin force-pushed the faster-sidereal-time branch from 0c6af84 to eee5e24 Compare September 27, 2026 16:59
@arahlin
arahlin force-pushed the faster-sidereal-time branch from eee5e24 to e5ff09a Compare September 27, 2026 17:11
@arahlin
arahlin force-pushed the faster-sidereal-time branch from e5ff09a to ad33d11 Compare September 27, 2026 17:19
@arahlin
arahlin force-pushed the faster-sidereal-time branch from ad33d11 to e5017c6 Compare September 27, 2026 18:03
@arahlin
arahlin force-pushed the faster-sidereal-time branch from e5017c6 to 1d14e05 Compare September 27, 2026 18:14
@arahlin
arahlin force-pushed the faster-sidereal-time branch from 1d14e05 to 830ded8 Compare September 27, 2026 18:21
@arahlin
arahlin force-pushed the faster-sidereal-time branch from 830ded8 to 6804334 Compare September 27, 2026 18:29
@arahlin
arahlin force-pushed the faster-sidereal-time branch from 6804334 to 8516cd6 Compare September 27, 2026 18:40
@arahlin
arahlin force-pushed the faster-sidereal-time branch from 8516cd6 to a3762e9 Compare September 27, 2026 19:21
@arahlin
arahlin force-pushed the faster-sidereal-time branch from a3762e9 to 03640e6 Compare September 27, 2026 19:58
@arahlin
arahlin force-pushed the faster-sidereal-time branch from 03640e6 to a42ea2d Compare September 27, 2026 20:24
@arahlin
arahlin force-pushed the faster-sidereal-time branch from a42ea2d to 829c185 Compare September 27, 2026 20:41
@arahlin
arahlin force-pushed the faster-sidereal-time branch from 829c185 to ff48a45 Compare September 27, 2026 20:55
@arahlin
arahlin force-pushed the faster-sidereal-time branch from ff48a45 to b975315 Compare September 27, 2026 21:14
@arahlin
arahlin force-pushed the faster-sidereal-time branch from b975315 to 7f3454e Compare September 28, 2026 01:57
@arahlin
arahlin force-pushed the faster-sidereal-time branch from 7f3454e to b975315 Compare September 28, 2026 02:14
@arahlin
arahlin force-pushed the faster-sidereal-time branch from b975315 to 37b6265 Compare September 28, 2026 05:35
@arahlin
arahlin force-pushed the faster-sidereal-time branch from 37b6265 to 6aa1efb Compare September 28, 2026 05:50
@arahlin
arahlin force-pushed the light-deflection branch 2 times, most recently from c045451 to 8b4c391 Compare September 28, 2026 15:46
The sun bends the incoming ray and neither package accounted for it.
That is the whole of the 14 mas by which azel2radec disagreed with
astropy, and of the 43 mas on the way back -- one missing term rather
than an accumulation of small ones. Adding it takes azel2radec to 0.79
mas and radec2azel to 0.84, so test_astropy's tolerance drops from a
tenth of an arcsecond to five milliarcseconds, resting on two models
agreeing rather than on room for a term one of them lacks.

rate_defl turns it on and defaults to 'never', like rate_dut1 and
rate_wobble: at ~20 mas for about a tenth of azel2bore it is opt-in, and
a caller who does not ask gets exactly what they got before. Only the
sun's position is cached at that rate; the deflection depends on where
the telescope points, so the per-sample part runs every sample, as
aberration does. It goes last in the forward chain and first in the
inverse, physics applying it before aberration and the forward chain
undoing an observation.

Aberration and deflection are one function rather than two, because
eraEpv00 costs 7.8 us a call and serves both -- aberration wants the
earth's barycentric velocity, deflection its heliocentric position, and
eraEpv00 returns both. Computing it twice doubled the whole transform
with both rates on 'always', 3202 ms against 1610 per 200k samples.
Sharing the pointing vector and summing the two rotation vectors into
one quaternion halves the per-sample cost again, 35.5 to 17.5 ns.

Summing drops the commutator of the two rotations, 1e-4 * 5e-7 rad or
0.002 uas. With deflection off the sum is the aberration vector
untouched: with the default rates the output is bit-identical to the
pre-deflection build over 240k values, including fast_aber=False, which
is not a small-angle rotation and keeps its own quaternion.

Switched on it costs 1.10x on azel2bore for at most ~20 mas at the
elongations a telescope observes at, which is the trade and the reason
the default is off.

eraLdsun does the arithmetic rather than the closed form 4.07 mas /
tan(elongation / 2), so a port of this shares the vendored ERFA and
agrees bit for bit by construction rather than by transcription. The
rotation always takes the small-angle form, unlike aberration's
fast_aber switch: the deflection never exceeds ~0.1 arcsec, so the
dropped cubic term is ~2e-20 rad, where the exact path measured nearly
twice the cost. The inverse is a sign flip -- eraLdsun has no
closed-form inverse, ERFA iterating in eraAticq -- which is second order
in the deflection, 1e-9 mas.

qp_reset_rates enumerates the rates one at a time and was missing
qp_reset_rate_defl, so reset_rates carried the previous chunk's sun
position into the next one, having been documented as the call to make
at the start of each chunk. TestResetRates now pins the property rather
than calling it and asserting nothing: with one rate set to 'once' the
correction is computed at the first sample and frozen, so after
reset_rates the answer has to match a freshly built QPoint. It is
parametrised over every rate the package reports rather than a list
written out in the test, so a rate added later is covered, and a
companion case checks the comparison can fail.

TestLightDeflection checks the term against erfa.ldsun rather than the
closed form, because get_sun returns a GCRS position with a distance
attached: compared to an ICRS coordinate, astropy applies parallax,
which for a body one au away swings the direction by 80 degrees.
test_astropy switches the term on in make_qpoint beside the IERS rates,
astropy applying it too, so leaving it off would compare two different
models.

The module's shared TOL_ARCSEC moves with the term rather than being
left loose, since every comparison in the file now has the deflection
switched on. make_qpoint lets a caller override one of the rates it
sets, which is how TestLightDeflection gets the same QPoint with the
term off to show what it is worth.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant