Fix the D_H node sign convention and guard the acos domain in Dcriteria

[AstroEloy]

Two correctness fixes in wmpl/Utils/Dcriteria.py, kept separate from any new
functionality so they can be reviewed on their own. One commit per bug, plus the
tests. calcDSH, calcDD, calcDN, calcDV and calcDVuncert keep their
current behaviour except for the acos guard described below.

1. calcDH was missing the node sign convention

calcDH computed the difference of the longitudes of perihelion as

pi21 = w2 - w1 + 2*asin(...)

without the rho = ±1 factor that calcDSH already applies. Southworth &
Hawkins (1963) specify the negative branch when |O2 - O1| > 180 deg, and Jopek
(1993) adopts that definition of Pi unchanged when replacing the perihelion
distance term. Without it, D_H depends on how the ascending node happens to be
written.

Worked example, q = 0.5 AU, e = 0.7, i = 15 deg, O1 = 10 deg,
w1 = 150 deg, w2 = 190 deg, with the second node written the two equivalent
ways:

	O2 = 350 deg	O2 = -10 deg
before	0.69872	0.26675
after	0.26675	0.26675

The same pair of orbits returned two different dissimilarities depending only on
the representation of the node.

Size of the error. rho enters only the last term, which is bounded by
((e1 + e2)/2)^2 * 4, so the discrepancy is bounded by 2*e. Maximising
|D_H(before) - D_H(after)| over the physical parameter box subject to
|O2 - O1| > 180 deg attains 1.9800, exactly 2*e_max for e_max = 0.99. A
plain Monte Carlo sample of 50214 pairs on that branch reaches 1.1478.

Consequence for shower association. Drawing close pairs whose nodes straddle
0/360 deg (nodes within 20 deg either side of 0, partner perturbed by 5% in q,
0.03 in e, 3 deg in i, 5 deg in w), 97379 of 2000000 pairs come out
associated at D_H < 0.10 under at least one node representation. Of those,
36.6% are placed on opposite sides of the threshold by the two representations.

The fix copies the rho block from calcDSH verbatim and multiplies the asin
term by it, so the two functions now differ only in the perihelion distance term.

2. Unguarded acos domain

calcDSH, calcDD and calcDH clip the argument of asin but not the argument
of acos. For the inclination term that argument reduces to cos^2 i + sin^2 i
for two identical orbits, which rounds above 1.0 in double precision, and
math.acos then raises ValueError: math domain error.

Measured over an 18001-point inclination sweep from 0 to 180 deg, comparing one
orbit with itself:

criterion	raises
D_SH	768/18001 = 4.27%
D_H	768/18001 = 4.27%
D_D	1197/18001 = 6.65%

D_D fails more often because calcDD also passes the angular separation of the
two perihelion directions (theta21) through an unguarded acos.

This is not reachable only by a synthetic test. It fires on any self-comparison,
on Monte Carlo clones of one orbit, and on every diagonal element of a distance
matrix, which is the first thing a shower-association or clustering run computes.

Fixed by applying the clipping idiom already used for asin in this file to every
acos argument, including theta21. The guard changes no value that previously
computed successfully: it engages only when the argument exceeds 1 by about one
ulp, where acos would otherwise have raised.

A self-comparison still cannot return exactly 0.0, because acos has an infinite
derivative at 1 and turns a 1 ulp argument error into about 2e-08 in the returned
angle. Measured worst case over the same sweep: 2.107e-08 for D_SH and D_H,
8.188e-09 for D_D. The tests therefore assert zero to within 1e-07.

Tests

New file wmpl/Utils/Tests/test_Dcriteria.py, 7 tests. Runs under pytest or
standalone via python -m wmpl.Utils.Tests.test_Dcriteria, matching the
convention in test_AlphaBeta.py.

Six of the seven fail on master; the seventh is a deliberate control.

• test_DH_pi21_matches_the_vector_definition — the strongest check. Jopek &
  Bronikowska (2017), P&SS 143, 43, eqs 4, 5 and 7 define Pi_21 through the
  mutual node of the two orbital planes, N = h_1 x h_2, with no arcsine and so
  no branch to choose. Comparing sin(Pi_21/2)^2, the quantity that actually
  enters D_H, over 200000 random pairs: with rho 200000/200000 = 100.00%
  agreement, max difference 1.5e-11; without rho 149453/200000 = 74.73%.
  Restricted to the |O2 - O1| > 180 deg branch: with rho 50547/50547 =
  100.00%, without rho 0/50547 = 0.00%. Not one straddling pair agrees
  without the factor. The test works through the public return value, not an
  internal quantity. It fails on master by 0.99933 against a maximum possible
  discrepancy of 1.
• test_DH_coplanar_limit_recovers_longitude_of_perihelion — for i1 = i2 = 0,
  Pi_21 must equal (w2 + O2) - (w1 + O1). Off by 1.3991 on master.
• test_DH_equals_DSH_when_perihelia_sum_to_unity — the two differ only in
  normalising the perihelion term by (q1 + q2), so they must agree exactly when
  q1 + q2 = 1. Disagree by up to 0.83942 on master.
• test_DH_invariant_under_node_representation — rewriting the second node
  across the 0/360 deg cut must not change D_H.
• test_DSH_invariant_under_node_representation — the control. Passes before
  and after, confirming calcDSH was already correct and that the suite is not
  merely asserting the new behaviour.
• test_acos_domain_holds_on_self_comparison and
  test_acos_domain_holds_for_near_identical_orbits — both fail on master with
  math domain error.

One caveat on the invariance test, worth knowing before reviewing it: the rho
test |O2 - O1| > pi is not 2pi periodic, so the convention holds only when both
nodes are given in [0, 360 deg). Adding a full turn to a node can push
|dOmega| past 3pi, where the rule breaks for calcDSH too. That is a
pre-existing property of the Southworth & Hawkins convention and is not addressed
here; the test reduces the node difference to its principal value instead, which
is the case that matters in practice.

Notes

• Reported figures were measured against these functions; the sampling schemes are
  stated above so they can be reproduced.
• No reformatting, renaming or reordering of existing code. The __main__ block
  is untouched.
• Companion PR adds further criteria and branches from master, not from this one,
  so the two are independent and this fix is not blocked by a scope discussion.

🤖 Generated with Claude Code

[AstroEloy]

A third independent primary source for the sign convention, found while checking a
formula for the companion PR.

Jenniskens (2008), Icarus 194, 13, eq. 3 states the rule explicitly: "the sign of
arcsin should be opposite when |Ω₁ − Ω₂| > 180°". So the factor is specified in
Southworth & Hawkins (1963), implied by the branch-cut-free vector construction of
Jopek & Bronikowska (2017) eq. 7, and stated in words by Jenniskens (2008). It is
not an optional refinement.

While there, I checked a discrepancy in that same equation, since it would affect
calcDSH and calcDH if it were real. Jenniskens eq. 3 prints
cos((i₁ − i₂)/2) in Π₁₂, where this codebase and Southworth & Hawkins have
cos((i₁ + i₂)/2). Tested against the mutual-node vector definition over 200000
random orbit pairs, comparing sin²(Π/2):

form	agreement	max difference
cos((i₁ + i₂)/2), as in this codebase	200000/200000 = 100.00%	1.5e-11
cos((i₁ − i₂)/2), as printed in eq. 3	0/200000 = 0.00%	1.0

No change is warranted — the codebase is right and eq. 3 as printed is a
typo. Recording it here so the next reader who finds that equation does not
"fix" the code to match it.

For completeness, section 2.1 of that paper carries two further transcription
errors in its restatement of the older criteria: eq. 9 writes the last D_H term
with (e₁ − e₂)/2 where (e₁ + e₂)/2 belongs, and eq. 5 writes the D_D
eccentricity term with (e₁ − e₂) in the denominator. The author's own criteria
in that paper are correct; it is only the review of prior work that is unreliable.

No code in this PR changed as a result of any of this.

🤖 Generated with Claude Code

[dvida]

Review pass: 1 commit pushed to this branch

Both bugs are real, the fix is right, and the test suite is the best in the four PRs I have reviewed this round. test_DH_pi21_matches_the_vector_definition checks against an independent vector construction rather than against the implementation, which is exactly the right way to settle a sign convention. 7/7 pass on the branch; on master exactly 6 fail and the seventh — the calcDSH control — passes, so the suite is not merely asserting the new behaviour. I reproduced that split before changing anything.

Independent confirmation

Measured myself rather than taken from the description:

D_H values that change, 200,000 random pairs	25.0 %, median change 0.069, max 1.28
Why 25 %	exactly the fraction with abs(dOmega) > 180 deg, as it must be
Associated pairs that flip across a threshold	3–6 % at D_H < 0.05 … 0.20, on shower-like pairs straddling the 0/360 deg cut
master raising on a self-comparison, 18,001 inclinations	D_SH 3.85 %, D_H 3.85 %, D_D 9.14 %; zero after the fix

Note your headline 36.6 % measures something different — two node representations of the same pair, not master against the fix. Both are honest; worth keeping them distinct so a reader does not compare them.

calcDSH and calcDD are unchanged to the last bit, against both master and this PR, over 100,000 pairs. The promise in the description holds.

The literature agrees. During the #97 review an agent verified over 200,000 pairs that your fixed calcDSH/calcDH reproduce the mutual-node construction of Jopek & Bronikowska (2017) eqs 4, 5 and 7 to 3e-11 rad. Courtot, Shober & Vaubaillon (2026) also state that D_H's Pi_1,2 is "defined the same way as D_SH", which is precisely what this restores.

682d926 — what I changed

The fix itself is untouched. This removes the duplication it left and writes down the convention.

• _clippedAcos() replaces the clipping idiom, which was written out at four sites.
• _mutualNodeAngles() returns I21 and pi21, which calcDSH and calcDH now share. Since the two criteria differ only in the perihelion normalisation — which is what this PR restored — the shared part belongs in one place.
• The two pre-existing # Name sure typos went with the blocks they sat in. (Fix spelling errors across the repo (codespell-assisted review) #101 would otherwise have caught them.)
• Both docstrings now state the branch rule and its limit.

One thing to check, because it is not free: calcDH moves by up to 1.9e-13 under the refactor. master wrote its arcsine argument as *1/math.cos(...) where calcDSH wrote *(1/math.cos(...)), and those associate differently in floating point. Folding them onto one expression costs one ulp. Numerically irrelevant, but it is a change, and it is the price of the two functions being genuinely the same code. Say the word if you would rather keep them separate and bit-identical.

Three tests added, 7 to 10: your worked example pinned at 0.26675 for both node representations; the 2π limitation asserted as documented behaviour; and _clippedAcos checked to pass through everything math.acos could already compute.

I also pushed a one-line caveat to the Pitfalls section of wmpl/Docs/OrbitSimilarity.md on the #97 branch, since that page does not exist on master yet.

---

The limitation you flagged — I tried to fix it and made it worse

You note that the rho test is not 2π-periodic and scope it out. I went and measured it, because it looked patchable:

	Agreement with Jopek's vector Pi_21, nodes carrying random extra turns
Your rho rule	78.0 %
rho from the principal value of dOmega	50.3 % — worse
Nodes in [0, 360 deg), your rule	100 %

So the rule cannot be repaired by reducing the angle. You were right to leave it, and I have pinned the behaviour in a test rather than pretending otherwise.

But the robust construction is already in this PR — _piFromVectorDefinition at test_Dcriteria.py:186. It builds Pi_21 from the mutual node N = h1 x h2, has no arcsine and therefore no branch to choose, and is invariant under any node representation. For nodes in [0, 360 deg) it agrees with the fixed code to 3e-11, i.e. numerically invisible.

Worth considering promoting it from the test file into calcDSH and calcDH, which would remove the limitation outright. I did not do it here, for two reasons: it changes calcDSH, which this PR deliberately leaves alone and uses as its control, and D_SH is the most-used criterion in the field — a change there deserves its own decision rather than riding along in a D_H fix. Your call, and the maintainer's.

Smaller notes

• Test naming: test_snake_case here, testCamelCase everywhere else in wmpl. Already raised on Add further orbit dissimilarity criteria, orbit classification and threshold methods #97 covering both PRs; I left yours alone and matched local style in what I added.
• The standalone runner (python -m wmpl.Utils.Tests.test_Dcriteria) still works and now reports 10.

Merges cleanly with master and with #97; 54 tests pass here, 208 on the #97 branch.

🤖 Generated with Claude Code