Add further orbit dissimilarity criteria, orbit classification and threshold methods

[AstroEloy]

Adds the orbit dissimilarity criteria the literature offers besides D_SH, D_D,
D_H and D_N, plus two new modules: one for criteria that classify a single
orbit, and one for choosing a criterion's threshold.

No existing D criterion is modified. calcDSH, calcDH, calcDD,
calcVgComponents, calcDN, calcDV and calcDVuncert are byte-identical to
master, and the __main__ block is untouched.

The branch is no longer purely additive, and an earlier revision of this
description wrongly said it was. The review pass added wmpl/Utils/OrbitConstants.py
and touched Orbit.py, MeanOrbit.py and REBOUND.py so that Dcriteria no
longer imports OrbitClassification (and with it scipy), and so the Gaussian
constant is defined once. Current contribution against upstream/master:
13 files, +4467 / −3. Whether that layering work should be split into its own
PR is still open — see the review thread.

Branches from master, not from the companion bug-fix PR, so the two are
independent.

New in Dcriteria.py

function	reference
calcRho1, calcRho2, calcRho5	Kholshevnikov, Kokhirova, Babadzhanov & Khamroev (2016), MNRAS 462, 2275
calcC	Neslusan (2002)
calcDR	Valsecchi, Jopek & Froeschle (1999), MNRAS 304, 743
calcDB, calcDT	Jenniskens (2008), Icarus 194, 13
calcDX	Rudawska, Matlovic, Toth & Kornos (2015), P&SS 118, 38
calcDVJopek	Jopek, Rudawska & Bartczak (2008), EM&P 102, 73
calcDACS	Asher, Clube & Steel (1993), MNRAS 264, 93, eq. 2
calcDSAC	Steel, Asher & Clube (1991), MNRAS 251, 632, as eq. 1 of the above

The rho metrics are the only criteria here that satisfy the triangle inequality
and the only ones defined for circular orbits. Both properties are tested with
D_D as the control: the triple

(0.955, 0.348, 168.20 deg, 184.34 deg,  29.66 deg)
(0.959, 0.441, 168.33 deg, 186.31 deg, 153.67 deg)
(0.965, 0.900, 168.27 deg, 184.55 deg, 209.86 deg)

breaks the D_D triangle inequality by 0.0679, while rho_2 has 0.2298 of slack
on the same triple. D_D also raises ZeroDivisionError for two circular orbits,
which the metrics handle.

calcDV, the unpublished Vida criterion already in the module, keeps its name so
the public API does not break; the new Jopek function is calcDVJopek.

All new functions are written against numpy rather than math and accept
scalars or arrays, which also removes the branch construct that produced the
D_H node sign bug in the companion PR.

New module OrbitClassification.py

These take one orbit and return a dynamical class, so they do not belong among the
pairwise criteria: calcTisserand (Whipple 1954), calcKresakK, calcKresakP,
calcAphelionDistance, calcOrbitalEnergy, and the two-parameter criteria
isCometaryQi, isCometaryEi, isCometaryKi, isCometaryPi from Jopek &
Williams (2013), MNRAS 430, 2377, eqs 8 to 12. Q-i and E-i are the two the paper
found most reliable.

New module Dthresholds.py

thresholdDr (Southworth & Hawkins 1963; Lindblad 1971), thresholdBreakPoint
(Neslusan, Svoren & Porubcan 1995), thresholdReliability (Jopek, Valsecchi &
Froeschle 1999) and thresholdRandomPairing (Jopek & Bronikowska 2017, P&SS 143,
43), whose threshold depends on the D function and the population as well as on
the sample size. Its twelve published coefficient pairs are exposed as
RANDOM_PAIRING_COEFFS; an unpublished combination raises ValueError rather
than substituting a default.

Worth flagging for users of thresholdDr: Jopek & Bronikowska measured the
probability of a coincidental pair at those thresholds as 1.0 for bolide samples
of 200 to 800 orbits, i.e. a "similar" pair is certain whether or not any orbits
are related. It runs 7 to 9 times looser than the threshold that delivers 0.01.
That is recorded in its docstring.

Thresholds are stated per criterion, never invented

Each docstring carries the criterion's published threshold or states that none was
published. Seven of the papers publish none, so no default is supplied. They would
not transfer anyway: rho carries units of sqrt(AU) with L = 1 AU, while
D_SH and D_D are dimensionless, and the random-pairing thresholds differ by a
factor of about 2.4 between D_SH and D_D at fixed sample size.

Verification

81 tests across test_DcriteriaAdditional.py (67) and test_Dthresholds.py (14).
Each criterion is checked against a value tabulated in its own paper, or, where
none is tabulated, against an analytic identity it must satisfy by construction.

Paper values. rho_1, rho_2 and rho_5 reproduce table 2 of Kholshevnikov
et al. (2016) for 96P/Machholz 1 and 2003 EH1 at 7415 BC:

	computed	paper
rho_1	0.0346	0.036
rho_2	0.0357	0.035
rho_5	0.0154	0.016

The residual is set by the precision of the published elements: perturbing them
within the half-width of their quoted rounding moves rho_2 over
[0.0354, 0.0360]. The test tolerance is 0.002, and every plausible transcription
error of the formulae was measured to miss by at least 0.0099, five times that —
flipping the sign of Delta in cos P gives 0.0450, dropping the e1*e2 term
gives 1.1875, using cos P in place of cos I gives 0.0475.

The paper's closed form for cos P was also checked independently against the dot
product of the two perihelion unit vectors, agreeing to 5.4e-16 over 20000 random
pairs.

thresholdRandomPairing reproduces eq. 14 of Jopek & Bronikowska as tabulated in
their column D_(14) to 7.7e-05, and the nine large-sample formulae track their
table 7 to 6.4e-04 absolute and 5.0% relative over 200 to 16000 orbits (they are
least-squares fits through those points, so exact agreement is not expected). The
criterion ratios of their table 8 follow from the formulae alone and match all six
published averages to two decimals. thresholdDr reproduces their column D_(1)
to 4.2e-05, inside the 5e-05 rounding floor of a four-decimal column.

calcTisserand gives 3.026 for 2P/Encke against the literature 3.03, and exactly
3 for an orbit identical to Jupiter's.

The calcDACS scale is published, and independently cross-checked here. Asher,
Clube & Steel (1993) print it in eq. 2 itself, which reads
D² = ((a₁−a₂)/3)² + (e₁−e₂)² + (2 sin((i₁−i₂)/2))². An earlier revision of this
description presented the argument below as a determination of the scale, which
overstated it — the value is published, not inferred.

It is worth keeping as a cross-check, because it is independent of the printed
equation and happens to be the only way the paper's own table can be brought to
bear. Table 1 cannot be recomputed from its own columns: the inclination is first
adjusted by Brouwer (1947) secular perturbation theory to the smallest value the
orbit ever reaches, and the eccentricity likewise, while the table reports observed
values. The semi-major axis is not adjusted and the other terms are non-negative,
so its term is a rigorous lower bound on every tabulated D, and it is attained
exactly by P/Encke and 1991 TB2. Those two facts bracket the scale at 3 AU, in
agreement with the printed value:

scale	outcome
1.0 / 2.0 / 2.5	lower bound violated by 0.3100 / 0.0500 / 0.0200
3.0	holds on all 26 rows, attained by 2 of them
3.5 / 4.0	bound slack everywhere, largest ratio 0.857 / 0.750

Both tests were run against substituted scales to confirm they reject every value
but 3.0. The docstring states that no secular adjustment is applied by the
function, so passing observed elements will not reproduce the paper's table.

Analytic identities, where no paper value exists: the triangle inequality for
the metrics with D_D as control; rho_5 <= rho_2 with rho_5 attaining the
minimum of rho_2 under a scan over both nodes and both arguments of perihelion;
D_R <= D_N against the existing calcDN; C against the angular momentum term
of rho_1; the D_V energy term against -1/(2a); D_B's C3 term taking the
smallest angular difference across the 0/360 deg cut; and the two Taurid forms
sharing their eccentricity and inclination terms, with their separately quoted
reference orbits agreeing (a1(1 - e1) = 0.3780 AU against the published
q1 = 0.375 AU, inside the range implied by the quoted precision).

A self-comparison cannot return exactly zero, for the same reason as in the
D_SH family: the square of the metric cancels to the 1e-16 rounding level and
the square root lifts that to about 1e-08. Measured floor over 5000 random orbits
is 3.7e-08 for rho_1, 6.0e-08 for rho_2 and exactly 0 for rho_5, so the
tests allow 1e-07.

D_B, D_T and D_X verified against their primary sources

An earlier revision of this PR carried caveats saying these three were transcribed
from a secondary source. The primary papers have since been checked and all three
formulae were correct, but the surrounding documentation was not. Corrected in
8389eab:

• D_X's weights are published, and are now the defaults. Rudawska et al.
  state them in section 2.2: w_lambda = 0.17, w_alpha = 1.20,
  w_delta = 1.20, w_v = 0.20. They were previously required arguments on the
  belief that the paper left them undefined. The velocity unit is not ambiguous
  either — the paper tabulates V_g in km/s — and a threshold is published:
  groups merge at D_X <= 0.15.
• D_T and D_B have published thresholds. Jenniskens section 3.2:
  D_T < 0.3 with D_B < 1.0 identifies parent bodies and siblings;
  D_T < 0.6 with D_B < 1.5 also captures bodies related through an earlier
  fragmentation. Both docstrings previously said no threshold was published.
• calcDT now takes q rather than a, following the paper's eq. 13, which
  writes the Tisserand parameter in q and e "because the observational errors
  in these orbital elements are smaller than those in a". The two forms agree to
  2.1e-14 over 200000 orbits, and eq. 13 stays finite at e = 1 where the
  (a, e, i) form evaluates to nan. This also brings D_T into the variables
  the rest of the module works in.

New paper values behind all three: T_J computes 3.0264 for 2P/Encke against
table 1's 3.03 and 4.5108 for 3200 Phaethon against 4.51; mean within-shower D_X
computes 0.061 / 0.101 / 0.070 for the Geminids, Perseids and Orionids against the
quoted 0.06 / 0.09 / 0.07; and the tabulated Taurid branch means give
D_X = 0.244, above the 0.15 merge threshold, consistent with the paper reporting
that its method keeps them separate.

A caution on section 2.1 of Jenniskens (2008), which restates the older
criteria — it should not be used as a source for them. Its eq. 9 writes the last
D_H term with (e1 - e2)/2 where (e1 + e2)/2 belongs; its eq. 5 writes the
D_D eccentricity term with (e1 - e2) in the denominator; and its eq. 3 writes
cos((i1 - i2)/2) in Pi_21 where cos((i1 + i2)/2) belongs. That last one was
checked against the mutual-node vector definition over 200000 pairs: the form in
this codebase agrees 100.00% of the time and the form as printed agrees 0.00%,
so no change is warranted. The author's own criteria, D_B and D_T, are correct.

D_V corrected against its primary source

1a2c8fc closes the last unverified formula, and it was a real error, not just a
documentation one. The paper defines the vectorial elements "in the units AU,
AU/day and the mass of the Sun", with mu = k^2. This module computed them with
mu = 1 — self-consistent, but on a different scale from the one the published
weights were measured on, and the weights are dimensional. For a Geminid-like orbit
the tabulated energy dispersion is 0.91% of the energy in the paper's units and
0.00027% of it in the units used before; only the first is a credible spread.
The angular momentum now carries a factor of k and the energy k^2.

The weights are published and are now the defaults. Eq. 5 gives
w = (2*sigma)^-2, so a pair differing by twice the dispersion in one element
contributes exactly 1 — that is what fixes the criterion's scale. Table 1's
dispersions for five stream ages plus the sporadic background are in
DV_DISPERSIONS, with calcDVWeights turning any of them into weights; the
default is the 4000-year set the paper used for its own search.

The 1.5 and 2 are deliberate, not transcription noise: the energy is
invariant under the principal secular perturbation and the third angular momentum
component semi-invariant, so the paper weights those two up beyond their
dispersions.

Thresholds are published (table 2, 99% reliability), now in DV_THRESHOLDS.
The column is headed D_V x 10^-1, which has to be read as the printed figures
being D_V scaled down by ten. Simulating pairs from one stream using table 1's
dispersions settles it: 97.1% of same-stream pairs fall below 3.128, as an
association threshold should, while 0.0% fall below 0.3128. Both readings are
asserted in the tests so the choice is visible to a reviewer.

Remaining provenance caveat

• calcC — Neslusan's volume is still not accessible. Its formula is
  corroborated by the description in Jopek, Rudawska & Bartczak (2008), which is a
  second source rather than the original, and the docstring says so. Whether that
  volume publishes a threshold is unchecked; none is supplied.
• thresholdBreakPoint is implemented as the method is described in the
  accessible literature; the 1995 paper's ADS scan carries no text layer, so its
  exact formalisation of where the turn is located could not be read. The break is
  taken as the bin of strongest negative curvature of the cumulative
  distribution, which on a synthetic stream uniform over [0, 0.05] plus a
  background over [0.30, 1.0] returns 0.04997.
• thresholdReliability considers pairs only. Larger minimum group sizes need
  a clustering algorithm, which is a separate concern from the threshold.
• Kresak's K and P are cited by author and year only. Their definitions and
  cometary limits are taken from Jopek & Williams (2013) eqs 11 and 12, which was
  verified directly; volume and page numbers for the 1967 and 1969 papers could
  not be checked and so are not quoted.

Tancredi (2014)

Added in 7a4735b: classifyTancrediComet, classifyTancrediAsteroid,
isTancrediResonanceProtected, calcHillRadius, calcResonanceSemiMajorAxis.

Applying a bare Tisserand cut to the asteroid population returns 7332 candidates,
over 97% of them on stable orbits. The criterion's value is the two conditions that
follow: the body must not be held away from Jupiter by a mean-motion resonance,
and its orbit must actually reach a giant planet. The paper reports a false
positive rate of 0.19 and false negative rate of 0.05 against direct integration.

Verified against the paper's own numbers: eq. 3 reproduces all eight resonance
semi-major axes of table 1 to 7e-04; eq. 2 reproduces all four quoted giant-planet
Hill radii to 2e-03; 233P/La Sagra, the one new "Comet" in an Asteroidal Orbit the
paper reports, gives T_Jup = 3.0860 against the quoted 3.086 and classifies as
'asteroidal orbit'; the quasi-Hilda orbit used to justify the 3.05 limit gives
3.0492 against the quoted 3.05; Chiron classifies as a Centaur and Encke as
Jupiter family; and the Hildas are correctly excluded as resonance protected.

The supplementary appendices are now implemented too (f2ca39a), so the criterion
runs from orbital elements alone and the two approximations noted earlier are gone.

Appendix A — calcMOID minimises the distance between a point on each orbit over
the two eccentric anomalies, restarting from the three further points the appendix
prescribes because the surface has more than one local minimum.
calcGiantPlanetMOIDs wraps it over the four giant planets, returning each distance
in that planet's own Hill radii. Verified on the orbits of figure A.1, chosen there
precisely because a single search falls into the second minimum: the four-start
search returns 0.240830 AU against 0.240833 from a 600×600 grid. Coplanar
circles of radius 1 and 1.2 give exactly 0.2, and the result is symmetric.
calcMinMOIDForTisserand traces the forbidden region, closing to zero at T = 3
and widening above it. Eq. A.1 is solved directly rather than through the cubic of
eq. A.3 it squares to, which carries two spurious roots.

Appendix B — calcResonanceWidth computes the libration half width from the
eccentricity, with calcLaplaceCoefficient evaluated from its integral definition
rather than the truncated series so it stays exact as alpha approaches 1. It
reproduces all five widths tabulated in table 1:

resonance	order	computed	table 1	ratio
4:1	3	0.00754	0.0075	1.005
3:1	2	0.02873	0.0287	1.001
5:2	3	0.02598	0.0260	0.999
7:3	4	0.02146	0.0215	0.998
2:1	1	0.11274	0.1127	1.000

One reading of eq. B.6 had to be settled for the last of those. Its final term
divides by j2, which is negative since j2 = -p. Taking the magnitude gives
0.10927, 3% low, and would need e = 0.3183 to match — outside the stated range.
Carrying the sign through gives 0.11274 at e = 0.3, the same eccentricity that
reproduces the other four, so the signed reading is the intended one. Noted at the
point of use.

The three outer resonances keep the fixed widths the paper itself adopts for them,
where the expansion no longer holds — the paper's choice, not an approximation
added here.

A note on that paper: the second form of eq. 1, rewriting the Tisserand parameter
in q and Q, is missing the planet's semi-major axis inside the radical as
printed, evaluating to 4.689 where the first form gives 3.030. The first form is
the one implemented.

Deliberately not included

Nothing. Every criterion named in the scope is implemented.

Review pass

149d51f records a review pass over the three modules. No formula changed. It
corrected one wrong statement about units in calcRho1; removed two
"no threshold is published" claims made about papers that were never accessible
(calcC, calcDVJopek) — the same mistake that turned out to be wrong for D_B,
D_T and D_X; reordered calcDX's arguments to (ra, dec, sol, vg) to match
calcDN, calcDR and calcDV, since all eight arguments are floats and the
previous order invited a silent wrong answer; replaced a test that recovered a
Tisserand parameter as 3.0 + D_T, valid only for T > 3; added isCometaryEi,
which the Q-i docstring recommended but which was not implemented; and named
D_B's dispersions and D_V's unverified coefficients at their point of use.

Out of scope

No pairwise matrices at survey scale. 10^6 orbits is 10^12 pairs and needs
chunking or a KD-tree pre-filter, which is a separate concern from this module.

🤖 Generated with Claude Code

[dvida]

Review pass: 4 commits pushed to this branch

Reviewed with every formula checked against its primary source where one could be reached. 78/78 tests passed as submitted; 214 pass now across Utils, Rebound and MetSim.

The headline: calcRho1/2/5, the centrepiece of this PR, are exactly right. I pulled Kholshevnikov et al. (2016) MNRAS 462, 2275 and compared term by term — eqs 11, 15 and 22 match the implementation character for character, and ρ₅'s quotient space matches its docstring. The triangle-inequality and circular-orbit arguments in the description hold.

One real bug, in calcDX.

Rudawska et al. (2015) eq. 2 writes the right-ascension term as [2 sin((α₁−α₂)/2) · cos δ₁]² — cos δ₁ scales the chord. The code had it inside the sine:

term_ra = w_ra*(d_vg + 1.0)*(2*np.sin((ra1 - ra2)/2.0*np.cos(dec1)))**2
#                                                 ^ cos(dec1) inside sin()

The two agree to first order in the radiant separation, which is exactly why the suite could not catch it:

ΔRA	δ₁	old	correct	error
20°	60°	0.03038	0.03015	+0.8%
90°	50°	1.0597	0.9958	+6.4%
180°	60°	1.5492	1.0954	+41.4%

All 78 tests still passed with it fixed, and on your own Taurid-branch case the difference is 0.00%. So the claim that every criterion was checked against a paper value does not quite hold for this one — the within-shower regime cannot distinguish the forms. 71764f7 fixes it and adds a test at a separation large enough to tell them apart, which fails on the old form.

In fairness to the practical impact: the error only grows for widely separated radiants, which sit far above the 0.15 merge threshold anyway, so shower-search outcomes are unlikely to change.

One claim in the description is more favourable than you thought

The D_ACS scale of 3 AU is published. The description devotes a section to inferring it from table 1 via a bounding argument. I read eq. 2 from the ADS scan of Asher, Clube & Steel (1993):

D² = ((a₁ − a₂)/3)² + (e₁ − e₂)² + (2 sin((i₁ − i₂)/2))²

It is printed right there. The elaborate inference was unnecessary — and it happens to have landed on the correct value. 71764f7 records the provenance at the constant. Their eq. 1 (Steel, Asher & Clube 1991) is also confirmed with no scale factor, q₁ = 0.375, e₁ = 0.82, i₁ = 4°.

71764f7 — other correctness items

• calcTisserand returns nan for an unbound orbit, and every comparison against nan is False, so classifyTancrediComet reached 'unclassified' by falling through every branch rather than by recognising the case. Now explicit, and documented.
• The Tancredi class intervals were open at both ends, so an orbit landing exactly on T = 2 or T = 3.05 matched no branch. Now half open, with a test that fails on the old form.
• Documented: calcMOID takes a where every criterion here takes q, with otherwise identical signatures; calcGiantPlanetMOIDs returns a dict, not the two scalars classifyTancrediAsteroid wants; calcDR is scalars-only and its w1 is a weight, not an argument of perihelion; the HILL_RADIUS_* constants are the paper's values, which is why calcHillRadius has no call site.
• calcC builds the angular momentum as √p while calcDVJopek builds it as k√p — the same physical quantity on scales differing by 58×. Neslušan (2002) is a 2000 workshop proceedings chapter I could not get behind; the code now says the convention is unresolved, which matters for anyone comparing a C value against a published one.

97e31b4 — layering, and a constant that should not be unified

Dcriteria.py imported three constants from OrbitClassification.py, whose own docstring says the two are kept apart, and dragged scipy.integrate/scipy.optimize into a module needing neither. They now live in a new wmpl/Utils/OrbitConstants.py that imports nothing; OrbitClassification re-exports them so nothing breaks. Verified: Dcriteria imported standalone no longer loads OrbitClassification or scipy.

The Gaussian constant appeared in four places — Orbit.py had it as a bare literal, MeanOrbit.py kept its own copy. Both now take it from the shared module.

I tried to unify Jupiter's semi-major axis and backed out. OrbitClassification uses 5.20336, REBOUND.py and Orbit.py use 5.204267. Substituting the latter moves the 1:1 resonance to 5.2043 and misses the 5.203 in Tancredi's table 1 — test_tancredi_resonance_semimajor_axes_match_table1 fails. The two belong to different sources and are not interchangeable. Both constants now say so, so the next person does not "fix" it.

1a95d17 — three manuals in wmpl/Docs/

OrbitSimilarity.md, OrbitClassification.md and DThresholds.md, following the REBOUND manual's format, indexed from the docs README. Each ends with a verification table recording what was checked against a primary source, what was confirmed only indirectly, and what could not be reached — so a later reader can tell the difference. All eight code examples were executed against the working tree.

The branch also merges current master, which brings in the wmpl/Docs folder these pages live in.

---

What I could not verify

Stated plainly rather than left to look checked:

Source	Outcome
Kholshevnikov et al. (2016)	Read — ρ₁, ρ₂, ρ₅ verified
Asher, Clube & Steel (1993)	Read (ADS scan) — eqs 1 and 2 verified, including the 3 AU scale
Jopek & Bronikowska (2017)	Read (arXiv:1609.03968) — all 12 threshold coefficient pairs and the validity ranges verified exactly
Courtot, Shober & Vaubaillon (2026) review	Read — D_SH, D_D, D_H, D_R, D_N, D_B, D_T, D_V and D_X verified, including D_V's asymmetric 1.5× on h_Z and 2× on the energy term
Tancredi (2014)	Unreachable. Paywalled; the author's own site no longer resolves and the Wayback Machine has no capture. The 3.05 limit and the class structure are corroborated by a 2025 paper Tancredi co-authored, but table 1, the Hill radii, table B.1 and eqs A.1 and B.5–B.7 are unverified. The implementation reproduces the tabulated resonance centres and widths to the printed precision, which is strong indirect evidence, not a check.
Jenniskens (2008), Icarus 194, 13	Unreachable. So I cannot confirm or refute the description's claim that its eqs 3, 5 and 9 contain three typos. That claim is currently unsupported in the PR as well, and worth softening or sourcing. Also: Icarus 194, 13 is titled "Meteoroid streams that trace to candidate dormant comets"; "Mostly dormant comets…" is a different 2008 Jenniskens paper in EM&P 102, 505. Worth checking which one the formulas came from.
Neslušan (2002)	Unreachable. Springer proceedings chapter, paywalled. The k-factor question above stays open.

Two smaller notes

• Test naming. This PR and Fix the D_H node sign convention and guard the acos domain in Dcriteria #96 use test_snake_case; every other test in wmpl uses testCamelCase. Renaming 78 + 7 tests is your call, so I left them and matched the local style in the ones I added.
• None of the twelve new criteria are wired into the __main__ CLI, which still exposes only D_SH, D_D and D_H. Probably deliberate, but worth a decision.

For the record, I also checked #96 against this: after both merge there is no inconsistency. Π₂₁ and the ρ family's angle P are two different, well-defined angles related by cos P = cos u₁ cos u₂ + sin u₁ sin u₂ cos I₂₁, each correct for its own criterion, and #96's fix makes calcDSH reproduce Jopek's mutual-node definition to 3e-11 rad over 200,000 random pairs.

🤖 Generated with Claude Code

[dvida]

@AstroEloy — could you review the four commits I pushed to this branch (e6e3afb..1a95d17) before I merge? Details are in the review comment above; this is the short list of what actually needs your judgement rather than just a look.

Please check these three, where I changed behaviour or contradicted the description:

1. The calcDX fix (71764f7). I read the right-ascension term as [2 sin((α₁−α₂)/2)·cos δ₁]² and moved cos δ₁ out of the sine argument. You have the primary source — does that match your reading of eq. 2? The old and new forms agree to first order, so the existing tests pass either way, which is why it slipped through.

2. The D_ACS scale. It is printed in eq. 2 of Asher, Clube & Steel (1993), which reads D² = ((a₁−a₂)/3)² + (e₁−e₂)² + (2sin((i₁−i₂)/2))². Your bounding argument from table 1 landed on the right value, but it was not needed — worth trimming that section from the description so it does not read as an inference.

3. Jupiter's semi-major axis is deliberately not unified. I tried folding A_JUPITER = 5.20336 into the 5.204267 that REBOUND.py and Orbit.py use, and test_tancredi_resonance_semimajor_axes_match_table1 fails — 5.20336 is what reproduces Tancredi's table. I documented the split at both constants instead. If you know which value Tancredi actually states, that would settle it properly.

Two open questions I could not close:

4. The Jenniskens (2008) typo claim. The description asserts its eqs 3, 5 and 9 are misprinted. Icarus 194, 13 is paywalled and I could not reach it, so that claim is unsupported in the PR. Could you cite the page, or soften it? Related: Icarus 194, 13 is "Meteoroid streams that trace to candidate dormant comets" — "Mostly dormant comets…" is a different 2008 Jenniskens paper (EM&P 102, 505). Worth confirming which one the D_B and D_T formulas came from.

5. The Neslušan C scale. calcC builds the angular momentum as √p, calcDVJopek as k√p — a factor of 58 between two functions in the same module. The 2000 proceedings chapter is behind a paywall. If you have it, does Neslušan carry the k?

Things I changed that you may simply disagree with:

6. The PR is no longer purely additive. It now adds wmpl/Utils/OrbitConstants.py and touches Orbit.py, MeanOrbit.py and REBOUND.py, to stop Dcriteria importing OrbitClassification (and with it scipy) and to define the Gaussian constant once. Happy to split that into its own PR if you would rather keep this one clean — the description still says +N/−0.

7. Test naming. This PR and Fix the D_H node sign convention and guard the acos domain in Dcriteria #96 use test_snake_case; the rest of wmpl uses testCamelCase. I left yours alone and matched local style in the tests I added. Your call whether to rename.

8. None of the twelve new criteria are exposed in the __main__ CLI, which still offers only D_SH, D_D and D_H. Deliberate?

The branch also merges current master, so it carries the wmpl/Docs folder that the three new manuals live in. 214 tests pass, and it still merges cleanly with master and with #96.

🤖 Generated with Claude Code

[AstroEloy]

Checked all four commits. The review is right on the substance, including on the
one item where it contradicts me. Answers to the numbered points below, and one
defect found in the review commits themselves.

First: three of the tests added in 71764f7 were never running

71764f7 added test_DX_right_ascension_term_scales_the_chord_not_the_angle,
test_tancredi_class_limits_are_half_open and
test_tancredi_unbound_orbit_is_not_silently_unclassified, but did not add them
to the list the module's standalone runner iterates over. pytest discovers by name
and ran all three; python -m wmpl.Utils.Tests.test_DcriteriaAdditional, the path
the module docstring documents, reported 64 tests where 67 exist.

The one silently skipped was the regression test for the calcDX fix in that same
commit — guarding against the form that was in the tree until that commit landed.
All three pass. Fixed in 0e84f10; the runner now reports 67, so 81 across both
files.

1. calcDX — you are right, I was wrong

Accepted. The decisive argument is not textual but geometric: every other term in
eq. 2 is a chord, 2 sin(Δ/2). cos δ is what turns a difference in right
ascension into a true angular distance, and it does that by scaling the chord.
sin(Δα·cos δ / 2) is not the chord of anything.

Your figures reproduce exactly — at ΔRA = 180°, δ = 60° the two forms differ by
41.4%. And your diagnosis of why it slipped through is the right one: the
within-shower regime cannot separate the forms, so my claim that every criterion
was checked against a paper value did not in fact cover this one. That was an
overstatement on my part, not just a gap in coverage.

2. D_ACS — agreed, and the description is corrected

You are right that the /3 is printed in eq. 2 and that the description read as
though it had to be inferred. Description updated: the scale is now stated as
published, with the bounding argument kept and relabelled as an independent
cross-check rather than a determination.

Worth keeping rather than deleting, for one reason: it is the only way the paper's
own table 1 can be brought to bear at all, since the table cannot be recomputed
from its own columns. But it is a check, not the source, and it should not have
been framed as one.

3. Jupiter's semi-major axis — settled, and Tancredi does state it

He states it implicitly but unambiguously. Table 1's 1:1 row is the value: eq. 3
gives a_res = (p/(p+q))^(2/3) · a_P, which for 1:1 is just a_P, and table 1
lists the 1:1 resonance at 5.203 AU.

So 5.203 is Tancredi's value, A_JUPITER = 5.20336 is correct for this module, and
the split from the 5.204267 in REBOUND.py and Orbit.py is a real difference of
sources rather than an inconsistency. Your decision to document the split instead
of unifying it was correct, and it now has a citation behind it rather than only a
failing test.

4. Jenniskens — the paper is the right one, and the typo claim is sourceable

Which paper: Icarus 194, 13, "Meteoroid streams that trace to candidate dormant
comets". That is what every docstring cites, with the matching DOI
10.1016/j.icarus.2007.09.016. The string "Mostly dormant" appears zero times
anywhere in this PR — code, tests or description — so there is nothing to correct.
D_B and D_T come from its eqs 15 to 18, and the σ values 0.13, 0.06 and 14.2°
from the paragraph below eq. 17.

The typo claim is in section 2.1 of that same paper, which restates the older
criteria before introducing the author's own. Three places:

• eq. 3 writes cos((i₁−i₂)/2) in Π₁₂ where cos((i₁+i₂)/2) belongs
• eq. 5 writes the D_D eccentricity term with (e₁−e₂) in the denominator
• eq. 9 writes the last D_H term with (e₁−e₂)/2 where (e₁+e₂)/2 belongs

Eq. 3 is settled numerically rather than by reading. Against the mutual-node
definition of Π, over 100,000 random orbit pairs:

form	agreement
cos((i₁+i₂)/2), as in this codebase	100.00%
cos((i₁−i₂)/2), as printed in eq. 3	0.00%

So no change is warranted, and the claim is not unsupported — it just was not
sourced in the description. Happy to add the section and equation numbers there if
you want it on the record in the PR rather than only in this thread.

The author's own criteria in that paper, D_B and D_T, are correct. It is only
its restatement of prior work that is unreliable, which is an argument for the
approach this PR took of verifying each formula against the paper that introduced
it.

5. The Neslušan k-factor — agreed, stays open

Your documentation of it is the right outcome. The one thing I can add is negative:
Jopek, Rudawska & Bartczak (2008) describe the criterion in their introduction as
"the difference between the orbital momentum vectors (per unit mass) of two orbits"
without fixing the units, so that paper does not settle it either. The 2000
proceedings chapter remains the only source that would.

6, 7 and 8 — still with @AstroEloy

These are calls I should not make on his behalf, so flagging rather than answering:
whether to split the OrbitConstants layering into its own PR, whether to rename
the 81 tests here and the 7 in #96 to testCamelCase, and whether to expose the
twelve criteria in the __main__ CLI.

On 6 specifically, the description has been corrected either way: it no longer
claims a purely additive diff and now states the real contribution against
upstream/master as 13 files, +4467 / −3.

---

81 tests pass here, and the branch still merges cleanly.

🤖 Generated with Claude Code

[AstroEloy]

Decisions on 6, 7 and 8, from @AstroEloy.

6 — the constants work stays here

No separate PR. Keep OrbitConstants.py and the three touched files on this
branch.

One correction to the justification, though, because it does not hold as stated.
The review says "Dcriteria imported standalone no longer loads
OrbitClassification or scipy". The first half is right and I re-verified it. The
second is not:

import wmpl                          -> scipy loaded: True  (512 submodules)
import wmpl.Utils.OrbitConstants     -> scipy loaded: True
import wmpl.Utils.Dcriteria          -> scipy loaded: True

wmpl/__init__.py runs pkgutil.walk_packages and imports the whole library, so
any absolute import of a wmpl submodule drags scipy in regardless of what
Dcriteria itself imports. The decoupling from OrbitClassification is real and
holds — that module is genuinely absent from sys.modules — but the scipy benefit
is currently nullified by the package's own __init__. Worth knowing before it is
cited as a reason elsewhere.

The description has already been corrected: it no longer claims a purely additive
diff, and states the real contribution against upstream/master as
13 files, +4467 / −3.

7 — test names stay as they are

Your note says you "matched the local style in the ones I added", but the three
tests in 71764f7 are test_DX_right_ascension_term_scales_the_chord_not_the_angle,
test_tancredi_class_limits_are_half_open and
test_tancredi_unbound_orbit_is_not_silently_unclassified — all with underscores.
You matched the file, not wmpl. So nothing is mixed: 81 tests here and 7 in #96,
all in one style, 0 in camelCase.

Given that, renaming 88 names buys nothing functional and the files are internally
consistent, so they stay. If the project would rather converge later, it is a
mechanical rename and better done across the suite at once than in this PR.

8 — every criterion is now in the CLI (f05d352)

It was not deliberate. The CLI computed three criteria while the module defines
sixteen; all sixteen are now reported, grouped by what they need, each with its
published threshold or a note that none exists.

• Orbital elements, always: D_SH, D_D, D_H, rho_1, rho_2, rho_5,
  C, D_V, D_B, D_T, D_ACS, D_SAC
• Radiant and speed, when available: D_N, D_R, D_X, and the Vida D_V
  when a corrected heliocentric direction is present. These come from a second
  trajectory file or from new --ra, --dec, --sol, --vg arguments
• Classification of each orbit, side by side: Tisserand, Kresak K and P,
  aphelion, the four Jopek & Williams tests, and the Tancredi class

Two things worth flagging:

A unit trap. Orbit stores v_g and v_h in m/s; every criterion here takes
km/s. Both are divided by 1000 on the way in. Running the same pair through both
paths — once from a pickle, once through --vg in km/s — gives identical D_N,
D_R and D_X, which is what checks it.

The layering is preserved. The classification functions are imported inside
__main__, not at module level, so importing Dcriteria still does not pull in
OrbitClassification. Verified after the change.

One pre-existing bug fixed along the way: traj_path2 was a required
positional, so the manual entry its own help text describes was unreachable, and
the if cml_args.traj_path2 is not None branch below could never be taken. It now
takes nargs='?', which is what that code was already written for.

The dead code after sys.exit() is left alone, as before.

Also: the three tests from 71764f7 were not running (0e84f10)

Flagged in my previous comment — they were absent from the list the standalone
runner iterates, so python -m wmpl.Utils.Tests.test_DcriteriaAdditional reported
64 of 67. Among them was the regression test for your own calcDX fix. Fixed; the
runner now reports 67.

---

81 tests pass, the CLI was exercised in four modes, and the branch still merges
cleanly.

🤖 Generated with Claude Code

[AstroEloy]

WIP. I will:

• add cli for OrbitClassification
• improve Dcriteria's cli to i) handle two manually introduced orbits, ii) allow user to indicate the criterion (or criteria) that want to use, iii) remove threshold information from the Dcriteria cli output, and iv) allow user to disable OrbitClassification when using Dcriteria.
• explore if it is possible and deserible to include orbital classification from NEOMOD and METEOMOD

[AstroEloy]

This branch is done from my side.

While working on this I also explored whether source-region models (NEOMOD, METEOMOD) could be plugged into the dynamical classification here or in a follow-up. Notes for whoever picks this up next:

• NEOMOD: not advisable to include as-is. It's calibrated for 15 < H < 28 (roughly sub-km to ~10 m objects), and using it for meteoroid-sized bodies means extrapolating well outside that range. At those smaller sizes, radiation pressure and other non-gravitational effects that NEOMOD's underlying integrations don't account for start to matter, so the extrapolated source probabilities wouldn't be trustworthy.
• METEOMOD: potentially a much better fit (its simulations use 1 m particles), but it needs more investigation, and including it means redistributing someone else's model/data. I've reached out to the main author (Miroslav Brož) to ask for permission and check for any caveats before doing anything further. Putting this on hold until at least that's resolved. Even with the author's approval, it might not be advisable to include it; this will be considered later.

[dvida]

@AstroEloy — checked all of it. You are right on both counts against me, and the branch is now up to date with master.

You found a real defect in my commits

Confirmed exactly: 64 registered, 67 defined, and all three missing were mine — including the regression test for my own calcDX fix, which is the one that most needed to run. I registered the tests I added to test_Dcriteria.py in #96 and then failed to do the same here. 0e84f10 is the right fix; I audited all three test files and registration is now complete (10 + 72 + 14).

And the scipy claim was misleading

Reproduced your result:

import wmpl.Utils.Dcriteria   ->  scipy loaded: True,  OrbitClassification loaded: True

wmpl/__init__.py runs pkgutil.walk_packages, so any real import of a wmpl submodule pulls the whole library in. The decoupling from OrbitClassification is genuine at the source level, but the scipy benefit is nullified in practice and I should not have stated it without that caveat.

What I verified of your new work

Check	Result
ParentBodySearch, the three original keys	max difference 0.000e+00 over 20,000 pairs; case-insensitivity works; unknown key still warns and falls back to D_SH
The T-i unit correction	Confirmed. The paper's 0.58 /AU × a_J = 3.0179, and T_paper * a_J matches calcTisserand to 8.9e-16. Comparing against 0.58 directly would have called 1510 of 2000 random orbits cometary
calcKresakK → calcWhippleK	Clean rename, no dangling references
CLIs	Exercised in five modes: both orbits manual, pickle plus manual, -c selection, --no_classification, and the OrbitClassification CLI

The np.isfinite guard against parabolic comets and the traj_path2 positional fix are both good catches that were not in scope and are worth having.

On the Jenniskens point: you are right that "Mostly dormant" appears nowhere in this PR. That title came from my own agent prompt, not from your code, so I flagged a non-issue. The eq. 3 numbers settle the typo claim — please do add the section and equation numbers to the description, since right now that evidence lives only in this thread.

---

9e96ec4 — merged master in, which was blocking

The branch was five commits behind and therefore did not have #96. Its own command line reported the pre-fix D_H: running PR #96's worked example through it gave

  D_SH        =     0.266746
  D_H         =     0.698722      <- pre-#96

where master gives 0.266746 for both. After the merge:

  D_SH        =     0.266746
  D_H         =     0.266746

Merged without conflict. 229 tests pass, and all three standalone runners report their full counts.

Nothing else outstanding from my side. Your WIP list — criterion selection, disabling classification, the source-region models — reads as follow-up work, and putting METEOMOD on hold pending Brož's permission is the right call on redistributing someone else's model.

🤖 Generated with Claude Code

[dvida]

Master (now with #101 and #102) merged into this branch again. The only conflicts were two docstrings in wmpl/Utils/ParentBodySearch.py, where #101 corrected "comapring" and "heighest" in text this branch had rewritten; kept this branch's text with both spellings fixed. Everything else merged automatically.

On the merged tree: 235 tests pass under pytest, the standalone runners report 10 / 72 / 14, and codespell finds nothing new in the changed files (its only hits are the astronomical "precess", which is pre-existing on master and not a typo).

🤖 Generated with Claude Code