Generalize SSB<->LISA coordinate transforms to arbitrary constellations (Taiji/TianQin), add numeric/real orbit support - #5394
Open
WuShichao wants to merge 23 commits into
Conversation
Generalizes the fixed circular-orbit LISA machinery in pycbc.coordinates.space to work with any 3-spacecraft triangular constellation orbit (analytic or numerical), for LISA/Taiji/TianQin, as groundwork for native single-link response and TDI generation. - pycbc/coordinates/space_orbit.py (new): NumericOrbits (scipy-only spline interpolation, duck-type compatible with lisaorbits.Orbits), constellation_frame (derives instantaneous centroid + rotation matrix from 3 spacecraft positions), and t_detector_from_ssb / t_ssb_from_t_detector (generalized light-travel-time root-finding). Validated to reproduce lisa_position_ssb / rotation_matrix_ssb_to_lisa / t_lisa_from_ssb exactly for the LDC manual's analytic circular orbit, and includes an optional (skipped if unavailable) cross-check against a real lisaorbits.Orbits instance. - pycbc/detector/space.py: register Taiji/TianQin in _space_detectors and add a _Generic_detector backend stub whose project_wave raises NotImplementedError, documenting that response/TDI is follow-on work. - pycbc/coordinates/__init__.py: export the new public names. pycbc.coordinates.space itself is untouched; existing test_coordinates_space.py and test_detector.py pass unchanged.
Cross-check constellation_frame/NumericOrbits against constellation geometries genuinely different from the LISA special case already covered: Taiji (heliocentric, leading the Earth by 20 degrees instead of trailing it) and TianQin (a fast rigidly-rotating triangle around a geocentric guiding center, in a plane fixed in inertial space rather than precessing over a year). The Taiji fixture follows the same first-order-in-eccentricity Keplerian expansion as the LDC manual's Eq. 48-52 above (Rubbo, Cornish & Poujade 2004, Phys. Rev. D 69, 082003), with Taiji's own arm length/eccentricity and a +20 degree lead angle. The TianQin fixture follows Hu et al 2018 (Class. Quantum Grav. 35, 095008), simplified to a pure circular heliocentric guiding center coincident with the Earth. Verifies: arm lengths match design values, Taiji's guiding center leads an idealized Earth proxy by exactly 20 degrees, TianQin's centroid coincides with that same Earth proxy, TianQin's constellation-plane normal is time-independent over a year (in contrast to Taiji's, which precesses), rotation matrices stay orthonormal, light-travel-time round trips are self-consistent, and NumericOrbits interpolates both shapes accurately (including TianQin's much shorter ~3.65 day rotation period, which needs a denser sampling grid than LISA/Taiji's yearly one).
to accept an arbitrary constellation orbit These are the functions actually called in real-time during parameter estimation, so the orbit-provider generality added in space_orbit.py needs to be reachable from here, not just from a parallel module nothing calls into. Each of the four gains an optional `orbit` argument (default None): with no orbit given, the code path is byte-for-byte the same as before (verified by the full existing test_coordinates_space.py suite passing unchanged); with an orbit given, the time-delay solve and rotation matrix are computed via space_orbit.t_detector_from_ssb / constellation_frame instead of the hard-coded circular LISA orbit, so any orbit provider (NumericOrbits, a real lisaorbits.Orbits instance, or an analytic Taiji/TianQin orbit) can be used instead. rotation_matrix_ssb_to_lisa and lisa_position_ssb themselves are left untouched -- they remain the LISA-specific default path, not something that needs to be made generic in place. New tests confirm the generalized path reproduces the existing hard-coded LISA path exactly when fed the LISA-equivalent analytic orbit (not just self-consistency), and that full ssb -> detector -> ssb round trips recover the original parameters for Taiji and TianQin.
These two compose lisa_to_ssb/ssb_to_lisa (already generalized) with ssb_to_geo/geo_to_ssb (Earth-specific, nothing to generalize there), but were silently dropping any orbit argument and always using the hard-coded circular LISA orbit on their LISA-side leg -- so Taiji/ TianQin round trips through the geocentric frame were not actually reachable despite ssb_to_lisa/lisa_to_ssb supporting them directly. Default (orbit=None) behavior is unchanged.
No code change. t_lisa_from_ssb, t_ssb_from_t_lisa, ssb_to_lisa, lisa_to_ssb, lisa_to_geo, geo_to_lisa keep their original names for backward compatibility, but now accept any constellation via `orbit` -- clarify in the module docstring that "lisa" in those names/params means "the constellation frame described by orbit" once orbit is given, not literally LISA.
parameter names, add orbit-file support to LISA/GEO transforms Two changes to the same classes, found together while wiring pycbc.coordinates.space's orbit= generalization through to the PE-config layer: 1. Pre-existing bug: SSBToLISA, GEOToLISA and SSBToGEO fully redefine __init__ (rather than extending their parent's) purely to share transform/inverse_transform via class-attribute assignment. Each called `super().__init__()`, but since their MRO is e.g. SSBToLISA -> LISAToSSB -> BaseTransform, that silently invoked LISAToSSB.__init__ (the next class in the MRO, not BaseTransform), re-running its default-parameter substitution with no arguments and clobbering any custom parameter names (`tc_lisa_param=...` etc.) the subclass's own __init__ had just set -- and, since `.outputs`/ `.inputs` were otherwise never set at all, `from_config` was already broken for these three classes (AttributeError). Fixed by calling `BaseTransform.__init__(self)` directly instead of `super().__init__()`, bypassing the MRO entirely. 2. New: LISAToSSB, SSBToLISA, LISAToGEO and GEOToLISA accept an optional `orbit_file` (ini: `orbit-file`), an HDF5 file readable by `space_orbit.NumericOrbits.from_file`, threaded through to coordinates.ssb_to_lisa/lisa_to_ssb/lisa_to_geo/geo_to_lisa's already-generalized `orbit=` argument. Default (no orbit-file) behavior is unchanged. This is what actually makes the orbit generalization reachable from a PE config, not just the Python API.
Reads an HDF5 file (or a group within it) containing 't' and 'positions' datasets (and optionally 'velocities'), matching NumericOrbits.__init__'s existing input contract. This is the file format any of LISA, Taiji, TianQin can be supplied in -- and is what pycbc.transforms' new `orbit-file` option (previous commit) loads to make a PE config's orbit provider selectable from the ini file. Also adds end-to-end tests confirming the transforms' orbit-file option reproduces pycbc.coordinates.space's orbit= functions exactly, both when constructing the transform class directly and via from_config, using the existing Taiji/TianQin analytic fixtures to build a real orbit file.
Neither constellation_frame nor NumericOrbits currently expose the per-arm unit vector and arm length between two spacecraft -- only the whole-constellation centroid and orientation. Single-link response formulas need exactly this: the direction light travels along one laser link, and that link's instantaneous length, for the sinc-type finite-armlength factor. link_vector(t, orbit, sc_emitter, sc_receiver) computes both directly from orbit.compute_position, for any orbit provider (analytic fixture, NumericOrbits, or a real lisaorbits.Orbits instance). It evaluates both spacecraft at the same time t (neglecting the light travel time between emission and reception), which is what the spatial projection factors in single-link response formulas need; the true retarded emission time, if ever needed, can be solved for the same way t_detector_from_ssb already does for the constellation centroid. Tested against all three constellation fixtures: arm length matches each mission's design value, unit vectors are normalized, opposite- direction links are antiparallel with equal length, and a vectorized multi-link call matches individual per-link calls.
LISA previously only had TestConstellationFrame (byte-level agreement with the pre-existing hardcoded pycbc.coordinates.space formulas) -- a stronger check, but not the same kind of test as the standalone arm-length/orthonormality/round-trip sanity checks Taiji and TianQin each get in their own dedicated test classes. Add the same direct checks for LISA for consistency. A "trails Earth by 20 degrees" check, analogous to Taiji's "leads by 20 degrees" one, is intentionally omitted: this fixture's T0 offset is calibrated to match the pre-existing TIME_OFFSET_20_DEGREES/real-epoch convention rather than a simple angular offset from a zero-phase toy reference, so it doesn't give a clean constant angle (verified numerically: ~84 degrees, not 20).
a toy zero-phase circular reference The previous version of this test class claimed a "trails Earth by 20 degrees" check couldn't give a clean result for this fixture and omitted it. That was based on comparing against an arbitrary, self-invented zero-phase circular Earth proxy (the same kind used for Taiji's fixture), which gives a constant but physically meaningless ~84 degrees for LISA -- because T0 is a time shift tuned to match a specific external convention, not an angle offset from that toy reference, so the comparison was never apples to apples to begin with. Comparing instead against the real Earth position (space.earth_position_ssb, real ephemeris) gives ~19-23 degrees, matching the range already documented in space.py's own TIME_OFFSET_20_DEGREES comment -- confirmed stable across several different multi-year time spans (not just near the particular epoch T0 happens to have been tuned at). Replaced the omitted check with this one.
against real_position_ssb directly Following up on the same issue found for LISA: the Taiji/TianQin test fixtures' guiding-center phase (alpha = OMEGA_0 * t) started at an arbitrary zero point with no connection to any real epoch, so comparing them directly against the real Earth position (space.earth_position_ssb) gave physically meaningless numbers (Taiji ~83-87 degrees instead of ~20; TianQin's centroid off by 161% of 1 AU) -- not because the generalized code is wrong, but because these fixtures were never anchored to reality in the first place. Fixed by anchoring both fixtures' phase to PHI0_REAL, real Earth's ecliptic longitude at GPS t=0 -- obtained directly from pycbc's own, already-existing space.earth_position_ssb (astropy-based real ephemeris), with no external orbital-element constants of any kind. This is a pure constant phase shift, so it doesn't change any existing test's pass/fail behavior (arm length, orthonormality, round-trip time delay, and the previous "vs an arbitrarily-phased circular Earth proxy" checks, which were updated to add the same PHI0_REAL shift to their own reference so they still hold exactly as before). With this anchor in place, Taiji now leads the real Earth by ~17.6- 21.6 degrees (verified over a dense grid spanning multiple decades, not just the randomly-sampled test times) -- close to the intended 20 degrees, with the residual explained by the fixture's guiding center still being an idealized circle rather than Earth's real eccentric orbit. TianQin's centroid now stays within ~4.3% of 1 AU of the real Earth position, the expected residual from the same circular approximation. Added test_leads_real_earth_by_approximately_20_degrees and test_centroid_close_to_real_earth_position to check this directly, with tolerance bands verified safe via a dense multi-decade grid, not just the specific randomly-sampled times a given test run happens to draw.
Before this, there was no way to do any Taiji/TianQin analysis in PyCBC without first obtaining or building a numeric orbit file -- orbit=None only has a working default for LISA (the existing lisa_position_ssb/rotation_matrix_ssb_to_lisa), and Taiji/TianQin have no built-in orbit at all despite being registered in pycbc.detector.space._space_detectors and despite space_orbit.py's machinery already being fully orbit-agnostic. These two classes close that gap with idealized analytic reference orbits (same functional form as the existing LISA formula, per Rubbo, Cornish & Poujade 2004 for Taiji; per Hu et al 2018 for TianQin), usable directly as an `orbit=` argument anywhere pycbc.coordinates.space's generalized functions accept one. Both default to anchoring their reference phase (kappa0) to the real Earth's ecliptic longitude at SSB t=0, via pycbc's own space.earth_position_ssb (astropy-based ephemeris, no external constants) -- so TaijiAnalyticOrbit()/TianQinAnalyticOrbit() with no arguments give a roughly-realistic-today constellation out of the box, while kappa0 can still be set explicitly for arbitrary or scenario-specific reference epochs. The lazy (function-local, not module-level) import of earth_position_ssb avoids a circular import with pycbc.coordinates.space, which imports from this module. These are explicitly documented as idealized reference orbits for prototyping/methods development (e.g. single-link response and TDI work) ahead of an official numeric orbit product, not a substitute for real mission ephemeris in science-quality analysis. Tested against the independently hand-written AnalyticTaijiOrbit/ AnalyticTianQinOrbit fixtures (bit-for-bit position agreement given the same kappa0, catching any transcription error in promoting the fixtures' math to a real class), the default kappa0's real-Earth anchoring, that a custom kappa0 actually changes the orbit, arm-length correctness, and usability as an orbit= provider through ssb_to_lisa/lisa_to_ssb. Also fixes a pre-existing omission: link_vector (added in an earlier commit) was missing from pycbc/coordinates/__init__.py's __all__ list.
TianQinAnalyticOrbit for consistency LISA previously had no explicit OrbitProvider class in space_orbit.py -- only the implicit orbit=None default baked into pycbc.coordinates.space's own functions (lisa_position_ssb/ rotation_matrix_ssb_to_lisa), with no importable object a user could pass directly to constellation_frame/link_vector or use as an explicit comparison baseline, unlike the two new Taiji/TianQin classes. LisaAnalyticOrbit() reproduces that existing orbit=None default exactly: same t0 (TIME_OFFSET_20_DEGREES), same formula, verified both against pycbc.coordinates.space's own hardcoded functions directly and against the independent AnalyticEqualArmOrbit test fixture. Unlike Taiji/TianQin, it deliberately does *not* default to a real-Earth- anchored reference epoch: t0 is already a real, pre-existing constant tuned for BBHx compatibility, and silently changing that default here would disagree with pycbc.coordinates.space's own default. Also factored the shared first-order-in-eccentricity Keplerian expansion (Rubbo, Cornish & Poujade 2004) out of TaijiAnalyticOrbit into a private _equal_arm_orbit_position helper, now used by both LisaAnalyticOrbit and TaijiAnalyticOrbit instead of duplicating it.
…iles Adds TestOptionalESAOemOrbitFiles, which fetches the official ESA lisa-orbit-files trailing-constellation dataset (CCSDS OEM format, via lisaorbits.OEMOrbits.from_included, network-cached ~600 kB total) and runs it through the same link_vector/constellation_frame machinery used elsewhere in this file. OEM files are given in the EME2000 (~ICRS) equatorial frame, not the SSB ecliptic frame space_orbit/space assume, so this test also exercises the ICRS -> BarycentricMeanEcliptic conversion (the same astropy transform already used by space.earth_position_ssb) that any caller must apply first. Skips gracefully if lisaorbits isn't installed or the download fails for any reason (no network, hash mismatch, etc.), mirroring the existing TestOptionalLisaorbitsDuckTyping pattern -- never depends on a local, machine-specific data path.
Real spacecraft ephemerides (e.g. ESA's CCSDS OEM files, read by lisaorbits.OEMOrbits) are most often published in an equatorial (ICRS/ EME2000) frame, not the SSB ecliptic frame space_orbit/space assume. Until now, using such a source required the caller to do their own ICRS -> BarycentricMeanEcliptic conversion before treating it as an OrbitProvider -- easy to get wrong (arm length comes out identical either way since it's a rotation invariant, but constellation_frame's centroid/rotation-matrix output does not, and would be silently wrong if the conversion were skipped). ICRSOrbitAdapter wraps any ICRS-frame compute_position/compute_velocity object and rotates its output into pycbc's ecliptic convention on the fly, using the same fixed rotation matrix pattern as space.earth_position_ssb (equinox-pinned, so it's a one-time computation, cached at module scope, rather than a per-call astropy transform). This makes a real lisaorbits.Orbits instance (OEMOrbits included) a true drop-in OrbitProvider with zero manual conversion. TestICRSOrbitAdapter validates the wrapper end to end, including a check against a fresh, independent astropy transform of a non-basis vector -- a round trip through the adapter's own cached matrix alone cannot catch a transpose/sign error in that matrix's construction (composing any orthogonal matrix with its own transpose is the identity regardless of whether the matrix itself is correct), which is exactly the bug this test caught during development. The existing ESA-OEM-orbit-files test is simplified to use the new adapter directly in place of its previous inline conversion.
…data The existing ESA lisa-orbit-files test only exercised space_orbit's lower-level geometric primitives (link_vector, constellation_frame) against real data; the actual PE-facing coordinate transforms (ssb_to_lisa/lisa_to_ssb -- sky localization + polarization angle) had only ever been checked against the analytic toy fixtures in TestSpaceAcceptsOrbitProvider, never against a real numeric orbit. Adds test_ssb_to_lisa_round_trip_with_real_esa_orbit, which feeds ICRSOrbitAdapter(oem_orbit) (the real ESA orbit, wrapped) directly as the `orbit=` argument to ssb_to_lisa/lisa_to_ssb and checks the full round trip recovers the original time/sky-location/polarization parameters. Same optional/skip-gated pattern as the rest of this class.
test_coordinates_space_orbit.py has several lisaorbits-gated tests (TestOptionalLisaorbitsDuckTyping, TestOptionalESAOemOrbitFiles) that have always skipped gracefully when lisaorbits isn't installed -- which, until now, was every leg of this CI matrix, since lisaorbits isn't a pycbc dependency and was never installed here. That made those tests dead weight on GitHub Actions: correct, but never actually exercised. lisaorbits requires Python >= 3.12, so this only runs on the 3.12/3.13 unittest legs (skipped on 3.11, where it can't install at all). It's a plain `pip install` into the already-solved `unittest` pixi environment for that job, not a change to pyproject.toml's dependency graph/lockfile, so it can't affect the other environments (search/inference/docs) or the 3.11 leg's solve.
Single-link response models need spacecraft velocity (e.g. Doppler/ finite-armlength corrections) and acceleration, which nothing in this module previously exposed for the analytic mission orbits, and NumericOrbits only went as far as velocity. - NumericOrbits.compute_acceleration: the analytic derivative of the velocity spline (itself either user-supplied or derived from position), exactly mirroring lisaorbits.InterpolatedOrbits' interp_ax = interp_vx.derivative() construction. - LisaAnalyticOrbit/TaijiAnalyticOrbit.compute_velocity/ compute_acceleration: exact analytic derivatives of the shared first-order-in-eccentricity Keplerian position formula, via two new module-level helpers (_equal_arm_orbit_velocity/_acceleration) that parallel the existing _equal_arm_orbit_position. Independently re-derived by hand from that same formula (chain rule through the constant-angular-frequency phase), not copied from any other implementation. - TianQinAnalyticOrbit.compute_velocity/compute_acceleration: both the Earth-like guiding center and the fast-rotation triangle term are uniform circular motion at their own constant angular frequency, so each differentiates independently (acceleration of a uniform circular term is simply -omega^2 times that same term's position). - ICRSOrbitAdapter.compute_acceleration: rotates through with the same fixed matrix as position/velocity (the ICRS -> ecliptic transform commutes with time differentiation). Precision is verified two ways: a dependency-free finite-difference cross-check (velocity vs. central difference of position, acceleration vs. central difference of velocity) for all three analytic orbits, and -- addressing the "match lisaorbits' precision" requirement directly -- an exact, machine-precision comparison against a real lisaorbits.EqualArmlengthOrbits instance, using its own angular frequency (sqrt(GM_SUN/a**3), not this module's EARTH_ORBIT_ANGULAR_FREQUENCY constant) to isolate a pure formula-for-formula check from that unrelated frequency-convention difference. NumericOrbits.compute_acceleration is checked against LisaAnalyticOrbit's exact acceleration via the same interpolation- accuracy pattern already used for position/velocity.
Benchmarked space_orbit against the real lisaorbits package (both analytic-orbit and numeric-orbit paths) at N=1e3/1e5/1e6-point calls, found and fixed three concrete inefficiencies: 1. _equal_arm_orbit_position recomputed np.sin(alpha)/np.cos(alpha) repeatedly inside its per-spacecraft loop (unlike the velocity/ acceleration versions, which already cached them) -- ~4x more transcendental evaluations than necessary over the full time array. Now caches sin(alpha)/cos(alpha) (and cos(alpha-beta_n)'s expansion via the angle-subtraction identity, avoiding a third transcendental call per spacecraft) once, matching the pattern already used correctly in velocity/acceleration. This was the one case where position was measurably *slower* than lisaorbits.EqualArmlengthOrbits at N=1e6 (1.15x) before this fix. 2. NumericOrbits built one scipy BSpline per (spacecraft, xyz) column (up to 9 independent scalar splines) instead of a single vector- valued spline over all columns at once. scipy shares the B-spline basis-function evaluation across all output columns of one BSpline call, so this was leaving a large, free speedup on the table: measured ~5-8x faster construction and evaluation using one combined spline (reshaping (N, M, 3) <-> (N, M*3) around scipy calls), identical numerical results (BSpline.derivative() applies per-column independently, verified). This was the single biggest win: went from roughly at-parity with lisaorbits.InterpolatedOrbits to ~4-6x faster. 3. ICRSOrbitAdapter rotated (N, M, 3) arrays by broadcasting a (3, 3) matrix directly over the leading axes via `@`, which numpy handles as many small batched matmuls; reshaping to (N*M, 3) first and rotating with a single 2D matmul measured ~35x faster in isolation (confirmed real, not noise, via a dedicated micro-benchmark separating the rotation cost from the wrapped orbit's own compute cost). Net effect, measured directly against the real lisaorbits package: LisaAnalyticOrbit is now consistently ~3-7x faster than lisaorbits.EqualArmlengthOrbits across all three of compute_position/ velocity/acceleration and all tested array sizes (no more N-dependent crossover to slower); NumericOrbits is ~4-6x faster than lisaorbits.InterpolatedOrbits (previously roughly at parity). ICRSOrbitAdapter's overhead over a wrapped lisaorbits object dropped significantly but is not fully eliminated -- it involves reshaping a transposed, non-contiguous array that lisaorbits itself returns, an inherent (if now much smaller) cost of the frame conversion, not something a smarter matmul alone removes. All numerical results are byte-for-byte compatible with before this change (pure algebraic rearrangement / restructuring, not a change in formula); full test suite, including the machine-precision cross-check against a real lisaorbits.EqualArmlengthOrbits instance, still passes.
…orbits Three additions, all requested directly: 1. NumericOrbits.from_oem_files(oem_1, oem_2, oem_3): a small, self-contained CCSDS OEM (v2.0) parser -- no dependency on lisaorbits or the external oem package it uses. Converts ICRS/ EME2000 km, km/s to pycbc's SSB-ecliptic convention in meters/GPS seconds (reusing ICRSOrbitAdapter's cached rotation matrix), and validates REF_FRAME/epoch consistency across the three files. Verified against real ESA lisa-orbit-files data (fetched via lisaorbits' own pooch cache, network-gated, no local path) to 4 mm absolute agreement with lisaorbits.OEMOrbits reading the exact same files -- i.e. matching to floating-point precision once the two are queried at the same physical instants. 2. NumericOrbits.from_triangle_dat_files(orbit_dir): reads Taiji Data Challenge / Triangle-Simulator's SCP*.dat/SCV*.dat format directly (AU / AU-per-day, no absolute timestamps -- caller supplies t0/dt). Both readers are covered by synthetic-fixture unit tests (a known analytic orbit re-encoded in each format, so they never depend on a local, machine-specific data checkout) plus one real-data, network-gated integration test for the OEM path. 2. LisaKeplerianOrbit/TaijiKeplerianOrbit: analytic orbits built from a genuine two-body Kepler ellipse per spacecraft (Newton-Raphson-solved eccentric anomaly, standard orbital-element-to-Cartesian projection), rather than LisaAnalyticOrbit/TaijiAnalyticOrbit's flat, first-order- in-eccentricity expansion -- closing the "no second-order-in- eccentricity option" gap relative to lisaorbits.KeplerianOrbits. The eccentricity/inclination relationship that keeps the three ellipses' mutual distance near the design arm length (a standard "tilted formation" equal-arm construction, delta=5/8 in _kepler_orbit_elements) was independently derived and checked here (attempted from-scratch symbolic/numeric re-derivation via sympy/ mpmath, consistent with though not as cleanly conclusive as intended given the numerical precision needed to isolate a second-order coefficient by direct root-finding) before being implemented with fully independent code (own variable names/structure, no code ported from lisaorbits) and validated: position/velocity/acceleration agree with a real lisaorbits.KeplerianOrbits instance (using its own angular frequency to isolate the comparison from the unrelated EARTH_ORBIT_ANGULAR_FREQUENCY-vs-sqrt(GM_SUN/a**3) difference already documented for the first-order case) to ~1e-2 m / 1e-8 m/s / 1e-14 m/s^2 -- floating-point precision at these scales -- and the derived eccentricity matches lisaorbits' own to 12 decimal places. Important, counterintuitive finding surfaced while testing this: confirmed directly against real lisaorbits output that EqualArmlengthOrbits (first order) keeps LISA's arm length constant to ~1e-5 m over a year, while KeplerianOrbits (second order) has a genuine ~0.2% arm-length variation and undershoots the nominal design arm length by a similar amount -- an inherent property of the true- Kepler-ellipse construction itself, not a numerical error, and not "more accurate" in the arm-length-constancy sense despite the name. Documented explicitly in LisaKeplerianOrbit's docstring so this isn't mistaken for a bug or for de facto superiority over the first-order default.
TianQin's known approximation gap isn't a "second order in eccentricity" issue the way LISA/Taiji's was: the fast-rotation triangle is a rigid rotation by design (not an independent Kepler ellipse per spacecraft, so there's no eccentricity to correct there), and neither lisaorbits (LISA-only) nor Triangle-Simulator's own GeocentricEqualArmAnalyticOrbit (also a circular-Earth approximation, confirmed by reading its source) offer a reference implementation of a more realistic TianQin guiding center to validate a Kepler-ellipse correction against. The real gap is the guiding center itself: it's approximated as a pure circular heliocentric orbit, neglecting Earth's real ~1.7% orbital eccentricity (larger than LISA/Taiji's own ~0.005-0.006 arm-length- driven eccentricity). Rather than build a closed-form eccentric-Kepler approximation of Earth's orbit, this adds a `guiding_center='real_earth'` option that uses the real astropy/JPL ephemeris directly (a new _real_earth_position_velocity() lazy-imported helper, using get_body_barycentric_posvel for position and velocity in one query, rotated into pycbc's ecliptic convention with the same cached matrix ICRSOrbitAdapter uses) -- strictly more accurate than any two-body Kepler ellipse, real or not, since it captures actual perturbations from the Moon and other planets too. compute_acceleration in this mode uses a small central finite difference of the real ephemeris velocity, since astropy does not expose Earth's acceleration directly; falling back to an idealized two-body acceleration would have reintroduced the very approximation this option exists to avoid. Backward compatible: guiding_center='circular' remains the default, identical to the pre-existing behavior (confirmed: all pre-existing tests pass unchanged). The new mode's guiding center matches space.earth_position_ssb directly (to sub-millimeter precision), and its velocity/acceleration match finite differences of its own position/ velocity. Verified the two modes deviate from each other by the expected ~2.5-6.5 million km (consistent with Earth's real eccentricity times 1 AU), not by an arbitrary phase-mismatch artifact or by nothing at all.
Generalizes the existing SSB-hub coordinate transform architecture (space.py) to a single lunar-surface detector, enabling arrival-time, sky-localization, and polarization-angle conversion between Moon, SSB, GEO, and space-based (LISA/Taiji/TianQin) frames -- the prerequisite for coherent LGWA+LISA+ET multiband parameter estimation. Per Tissino et al. 2026 (arXiv:2606.04918), sky position and polarization are effectively invariant to the observer's location within the Solar System, so the Moon frame uses an identity rotation relative to SSB; only arrival time needs the detector's precise position, optionally including real lunar libration via the lunarsky package (lazily imported, degrading cleanly to the Moon's barycenter when absent). - pycbc/coordinates/moon.py: rotation_matrix_ssb_to_moon, moon_site_position_ssb, t_moon_from_ssb/t_ssb_from_t_moon, ssb_to_moon/moon_to_ssb, and composite moon_to_geo/geo_to_moon, moon_to_lisa/lisa_to_moon (the latter transparently accepting orbit= for Taiji/TianQin). - pycbc/coordinates/space_orbit.py: generalize _real_earth_position_velocity into _real_body_position_velocity(t, body), reused by moon.py for the Moon's real barycenter position. - pycbc/transforms.py: MoonToSSB/SSBToMoon, MoonToGEO/GEOToMoon, MoonToLISA/LISAToMoon transform classes for use in inference configs. - pycbc/detector/space.py: register LGWA as a Generic-backend placeholder (response function is future work, matching the existing Taiji/TianQin placeholders). - CI: install lunarsky (no Python-version restriction, unlike lisaorbits) so the optional site-position tests get real coverage.
![github-advanced-security[bot]](https://avatars.githubusercontent.com/in/57789?s=60&v=4){kind=small}
There was a problem hiding this comment.
qlty found more than 20 potential problems in the proposed changes. Check the Files changed tab for more details.
WuShichao
force-pushed
the
feature/orbit-coordinate-generalization
branch
from
August 2, 2026 19:34
bc9b982 to
526265e
Compare
Contributor
|
@WuShichao This seems a huge amount of code addition for what is a relatively small motivation. This is also Claude-written (it could be another LLM, but is definitely LLM-written and looks like Claude). I've no issue with this, but in line with scipy's policy (https://docs.scipy.org/doc/scipy-1.18.0/dev/conduct/ai_policy.html) we need to make sure that users are reviewing code before submitting and understand all of what is being added. So can I ask you (without using Claude/similar) to explain what all this is doing and why such a large amount of code is actually needed. |
Member
Author
Thanks for the comment! My guiding principle is that I lead the development of the core algorithms and code (building upon my previous PR for coordinate system transformation #4289), while relying on Claude to assist in writing the extensive unit test code—though I define the requirements for each test myself, and refine my handwritten code. Currently, I am still manually and repeatedly testing the PR (it is far from being merged), and I welcome reviews and discussions regarding the details. After all, I'm not a professional software engineer, so I do need an LLM to provide some assistance. |
spxiwh
removed their request for review
August 3, 2026 07:23
Member
|
@vikasjadhav-why @acorreia61201 Please look at this PR carefully. Check the (1) interfaces (2) generality (3) unit test coverge. Also check if this is consistent with what you've been playing around with (e.g. do you get consitant results. |
Member
|
@WuShichao I also think this is quite a large amount of code. I would ask if you can explain in more detail the struture here and also consider if this can be split up or reduced in size (e.g. do a human check to see that this actually makes sense for a single PR and that there isn't any code duplication / bloat). |
The existing native readers (from_oem_files/from_triangle_dat_files) had no counterpart writer: PE configs select a numeric orbit via a pycbc-schema HDF5 file (orbit-file, read by from_file), but nothing produced one from an orbit loaded through those readers except writing the t/positions datasets by hand. NumericOrbits.to_file writes exactly the schema from_file reads (t, positions, optional velocities, optional HDF5 group), so any orbit built from an OEM/Triangle-Simulator/ lisaorbits source can be converted to a usable orbit-file in one call. Velocities are only written if the instance was built with them explicitly, so a round trip through from_file reproduces the original construction (spline-derived vs. explicit) rather than adding a second layer of interpolation on re-sampled velocities; acceleration is never stored, consistent with it always being the velocity spline's analytic derivative elsewhere in this class. NumericOrbits.from_lisaorbits_file reads the file format produced by lisaorbits.Orbits.write directly (tcb/x positions in ICRS, uniform TCB sampling from t0/dt/size attributes), without requiring lisaorbits itself to be installed -- the same self-contained-parser approach as from_oem_files, reusing its cached ICRS -> ecliptic rotation. Only position data is used, matching from_oem_files' treatment of velocity/ acceleration. Verified end to end: a lisaorbits-native file loaded via from_lisaorbits_file, written back out via to_file, and pointed to by transforms.SSBToLISA(orbit_file=...) reproduces the same transform output as querying the orbit directly; round-trip position/velocity/ acceleration checks pass for both the with- and without-explicit- velocities cases, and from_lisaorbits_file's output matches a real lisaorbits.EqualArmlengthOrbits.write() fixture to ~1e-3 m with the correct 2.5e6 km LISA design arm length.
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
4 participants
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.
Standard information about the request
This is a: new feature
This change affects: inference
This change: has appropriate unit tests, follows style guidelines (See e.g. PEP8), has been proposed using the contribution guidelines
This change will: require additional dependencies (optional:
lisaorbits, only for validating against real/numeric orbit data -- everything else has no new hard dependency)Motivation
This PR bundles two related but distinct pieces of work, kept together because the second is a direct prerequisite for a companion PR:
pycbc.coordinates.spacewas hardcoded to LISA's own analytic circular orbit. Extending PyCBC to other space-based GW missions (Taiji, TianQin) and to numeric/real orbit data (e.g. ESA's own LISA orbit files, mission-specific simulators) needs these transforms to accept an arbitrary orbit provider instead of assuming LISA specifically.pycbc.coordinates.moon). This is not an instance of (1) -- the Moon is a single detector, not a multi-spacecraft constellation with an orbit provider -- but a separate extension of the same SSB-hub architecture, included here because a companion PR (LGWA response backend + relative-binning support) depends on it.Of the non-test code in this PR, roughly 2400 lines implement (1) and roughly 600 lines implement (2).
Contents
OrbitProvider-style analytic orbit classes:TaijiAnalyticOrbit,TianQinAnalyticOrbit, and (for parity/consistency) a productionLisaAnalyticOrbit.space.ssb_to_lisa/lisa_to_ssb/t_lisa_from_ssb/t_ssb_from_t_lisa, andlisa_to_geo/geo_to_lisa, generalized to accept anorbit=/sc=argument selecting the constellation and spacecraft labels; default behavior (noorbitgiven) is unchanged from before this PR.pycbc.coordinates.space_orbit.NumericOrbits: load an arbitrary numeric orbit (interpolated) from a plain time/position array or from a saved HDF5 file (.from_file).space_orbit.ICRSOrbitAdapter: wraps any ICRS-frameOrbitProvider(including a reallisaorbits.Orbitsinstance) into pycbc's own SSB-ecliptic convention, so numeric/real orbit sources are usable as drop-in providers everywhere inpycbc.coordinates.space/space_orbit.space_orbit.link_vector,compute_velocity/compute_accelerationadded throughout, for single-link response inputs and higher-derivative quantities.pycbc.transforms:orbit-fileini option added to the LISA/GEO transform classes, so an inference config can point directly at a numeric orbit file.SSBToLISA/GEOToLISA/SSBToGEOwere clobbering custom parameter names in some configurations.lisaorbitson Python >= 3.12 legs, enabling several tests that validate against real, independently-produced orbit data (ESA's own OEM orbit files, Taiji's Triangle-Simulator TDC-II products) -- these tests skip gracefully whereverlisaorbitsisn't installed.pycbc.coordinates.moon, the SSB-hub coordinate-system unification for lunar detectors -- arrival-time/sky-localization/polarization conversion between Moon, SSB, GEO, and LISA-family frames, with an identity rotation (per Tissino et al. 2026, arXiv:2606.04918, Eq. 7-9) and real lunar libration via the optionallunarskypackage.Links to any issues or associated PRs
Base branch for a companion PR adding LGWA (lunar GW detector) response-function support and relative-binning compatibility.
Testing performed
test/test_coordinates_space_orbit.py,test/test_transforms.py,test/test_coordinates_space.py,test/test_coordinates_moon.pyall pass.orbit=generalization is a zero-behavior-change default: with noorbitgiven, output is bit-identical to pre-PR behavior.NumericOrbits/ICRSOrbitAdapteragainst three independent, real data sources:lisaorbits's ownInterpolatedOrbits, ESA's officially-publishedlisa-orbit-files(CCSDS OEM format, fetched vialisaorbits.oem.OEMOrbits), and Taiji's Triangle-Simulator Data Challenge II orbit products.Additional notes
This branch is the base for a companion PR adding LGWA (lunar GW detector) response-function support and relative-binning compatibility, which depends on the Moon-frame coordinate unification included here (item (2) above). Roughly 40% of this PR's total diff is test code (
test/test_coordinates_space_orbit.pyalone cross-validates against three independent real orbit-data sources).