Design proposal: GCRF/ICRF Frames

[translunar]

Wrote this with Claude, but I supervised it pretty closely. Would love feedback before I move on to writing an implementation plan.

Add ICRS Orientation and GCRF/ICRF Frames to ANISE

Summary

ANISE currently treats ICRF as an alias for J2000 (id_from_orientation_name("ICRF") returns J2000), matching SPICE behavior. This conflates two distinct reference systems:

• J2000 (EME2000): Earth mean equator and dynamical equinox at J2000.0 epoch
• ICRS: International Celestial Reference System, defined by extragalactic radio sources (quasars) via VLBI

The two are related by a constant rotation of ~23 milliarcseconds (~0.7 meters at Earth's surface), known as the IERS 2006 frame bias (IERS Conventions 2010, Chapter 5; Hilton et al. 2006).

This proposal adds ICRS as a first-class orientation in ANISE, enabling:

• GCRF — Geocentric Celestial Reference Frame (Earth-centered, ICRS axes)
• ICRF — International Celestial Reference Frame (SSB-centered, ICRS axes)

The implementation follows the existing ECLIPJ2000 precedent: a hard-coded constant rotation from J2000 with no kernel data required.

Motivation

ICRS is the IAU's primary celestial reference system since 1998. GCRF is the standard for precise geocentric work (e.g., GNSS, precise orbit determination). SPICE's simplification of ICRF ≈ J2000 is documented as intentional but approximate. ANISE, as a modern reimplementation, has the opportunity to model this correctly.

For most applications the ~0.7m bias is negligible. But for users who need it — precise geodesy, inter-agency data exchange, or simply physical correctness — the current workaround is manual matrix application outside of ANISE's frame system. This proposal makes it a first-class citizen.

API Surface

Orientation Constant

New orientation ID in constants::orientations:

/// ICRS orientation axes (International Celestial Reference System).
/// Related to J2000 by the IERS 2006 frame bias.
pub const ICRS: NaifId = 22;

ID 22 is the next sequential ID after the SPICE built-in inertial frames (1–21). This is an ANISE-specific extension; SPICE does not define a separate ICRS orientation. The final ID choice is at the maintainer's discretion.

Frame Bias Constants

Three constants from SOFA iauBi00 (Chapront et al. 2002; Mathews et al. 2002), stored in arcseconds to match the SOFA source:

/// SOFA iauBi00: longitude bias
pub const FRAME_BIAS_DPSIBI_ARCSEC: f64 = -0.041775;
/// SOFA iauBi00: obliquity bias
pub const FRAME_BIAS_DEPSBI_ARCSEC: f64 = -0.0068192;
/// SOFA iauBi00: right ascension correction
pub const FRAME_BIAS_DRA0_ARCSEC: f64 = -0.0146;

The frame bias matrix is (SOFA iauBp00; IERS 2010 eq. 5.4):

B = R1(-depsbi) * R2(dpsibi * sin(EPS0)) * R3(dra0)

where EPS0 = 84381.448 arcseconds (J2000.0 obliquity), already available in ANISE as J2000_TO_ECLIPJ2000_ANGLE_RAD.

This is the exact formulation (USNO Circular 179, eq. 3.4), not the small-angle approximation (eq. 3.3). The difference is ~6 nanometers at Pluto, so both are equivalent in practice, but the exact form is simpler code (compose three rotation matrices ANISE already has) and matches SOFA.

Frame Constants

/// Geocentric Celestial Reference Frame (Earth-centered, ICRS orientation)
pub const GCRF: Frame = Frame::new(EARTH, ICRS);
/// International Celestial Reference Frame (SSB-centered, ICRS orientation)
pub const ICRF: Frame = Frame::new(SOLAR_SYSTEM_BARYCENTER, ICRS);

Following the EME2000 = Frame::new(EARTH, J2000) precedent — the conventional names encode the origin.

No other body-centered ICRS constants (e.g., MOON_ICRS) are defined initially. Users can construct them with Frame::new(MOON, ICRS) if needed.

Name Mappings

orientation_name_from_id(22) returns "ICRS" (the orientation system name).

id_from_orientation_name maps:

"ICRS" | "GCRF" | "ICRF" => Ok(ICRS),

The existing "ICRF" alias is moved from J2000 to ICRS.

Python Bindings

Frames.GCRF        # Frame(EARTH, ICRS)
Frames.ICRF        # Frame(SSB, ICRS)
Orientations.ICRS  # 22

String parsing (Frame.from_name("Earth", "ICRS")) resolves via the updated name mappings.

Correctness Criteria

These define what "correct" means for this feature. Implementation must satisfy all of these.

1. Round-trip identity

Rotate a state from J2000 to ICRS and back. The result must equal the input to machine precision (~1e-15 relative error).

2. Frame bias magnitude

Transform a position vector at Earth's surface (~6378 km along X in J2000) to ICRS. The position difference magnitude must be approximately 0.7 meters (23 mas × 6378 km). This validates that the bias is neither zero (broken) nor wildly wrong (swapped frames).

3. SOFA cross-check

Compute the 3×3 frame bias matrix independently using the SOFA formulation:

B = R1(-depsbi) * R2(dpsibi * sin(eps0)) * R3(dra0)

with the iauBi00 constants. Compare all 9 matrix elements against what almanac.rotate(EME2000, GCRF, epoch).rot_mat produces. Maximum element-wise difference must be < 1e-15.

4. Chain composition

Transform a GEO orbit from GCRF to ITRF93 (internally: ICRS → J2000 → ITRF93). Compare against the same orbit transformed from EME2000 to ITRF93. Position difference must be ~0.7m (the frame bias), not zero (which would mean ICRS isn't being applied) and not large (which would mean a bug in the chain).

This test requires a high-precision Earth orientation BPC file (e.g., earth_latest_high_prec.bpc).

5. Name resolution

• id_from_orientation_name("ICRS") returns 22
• id_from_orientation_name("GCRF") returns 22
• id_from_orientation_name("ICRF") returns 22
• orientation_name_from_id(22) returns "ICRS"
• Frame::from_name("Earth", "GCRF") returns GCRF
• Frame::from_name("Earth", "ICRF") returns Frame::new(EARTH, ICRS) (same orientation as GCRF, same origin)
• Frame::from_name("SSB", "ICRF") returns Frame::new(SSB, ICRS) (which equals ICRF)

6. Angular velocity

angular_velocity_rad_s(GCRF, EME2000, epoch) must either return zero or error consistently with how ECLIPJ2000 behaves (both are constant rotations with rot_mat_dt: None).

Implementation

anise/src/constants.rs

• Add ICRS orientation ID (22) with doc comment
• Add the three FRAME_BIAS_*_ARCSEC constants with SOFA reference
• Add ICRS => Some("ICRS") to orientation_name_from_id
• Add "ICRS" | "GCRF" | "ICRF" => Ok(ICRS) to id_from_orientation_name; remove "ICRF" from the J2000 match arm
• Add GCRF and ICRF frame constants

anise/src/orientations/rotate_to_parent.rs

Add a new branch after the ECLIPJ2000 case, following the identical pattern:

} else if source.orient_origin_id_match(ICRS) {
    let das2r = std::f64::consts::PI / (180.0 * 3600.0);
    let dra0   = FRAME_BIAS_DRA0_ARCSEC * das2r;
    let dpsibi = FRAME_BIAS_DPSIBI_ARCSEC * das2r;
    let depsbi = FRAME_BIAS_DEPSBI_ARCSEC * das2r;

    // B = R1(-depsbi) * R2(dpsibi * sin(EPS0)) * R3(dra0)
    // B transforms a vector expressed in ICRS to J2000 coordinates;
    // transpose gives J2000 -> ICRS
    let b = r1(-depsbi)
          * r2(dpsibi * J2000_TO_ECLIPJ2000_ANGLE_RAD.sin())
          * r3(dra0);

    return Ok(DCM {
        rot_mat: b.transpose(),
        rot_mat_dt: None,
        from: J2000,
        to: ICRS,
    });
}

• rot_mat_dt: None matches ECLIPJ2000 (constant rotation, no time derivative)
• from: J2000, to: ICRS matches the convention used by ECLIPJ2000 and BPC rotations: the rot_mat transforms vectors from the parent to the child frame
• Reuses J2000_TO_ECLIPJ2000_ANGLE_RAD for the obliquity (same value SOFA uses for EPS0)

anise/src/orientations/paths.rs

Two changes to teach the path finder that ICRS's parent is J2000:

Initial parent resolution (~line 98): Add ICRS before the BPC → planetary → EPA fallback chain:

let mut inertial_frame_id = if source.orient_origin_id_match(ICRS) {
    J2000
} else {
    match self.bpc_summary_at_epoch(source.orientation_id, epoch) {
        // existing chain
    }
};

Intermediate node handling (line 115): Add ICRS alongside ECLIPJ2000:

if inertial_frame_id == ECLIPJ2000 || inertial_frame_id == ICRS {
    inertial_frame_id = J2000;
    of_path[of_path_len] = Some(inertial_frame_id);
    of_path_len += 1;
}

Python bindings

Expose GCRF, ICRF, and Orientations.ICRS in the Python constants module.

Alternatives Considered

Precomputed constant matrix

Store the 9 matrix elements of the frame bias DCM instead of computing r1()*r2()*r3() at each call. Marginally faster (avoids 4 trig calls per rotation), but harder to audit — 9 opaque floats versus 3 readable angle constants traceable to SOFA. This optimization could apply equally to ECLIPJ2000, which also recomputes r1() from a constant angle at every call. Both could be changed together in a future PR. Not proposed here to keep the changeset minimal and consistent with ANISE as currently authored.

EulerParameter dataset

Encode the frame bias as a quaternion in an EPA file loaded by the almanac. Fully data-driven with no changes to the rotation or path-finding logic — the existing EPA fallback handles everything. However, every GCRF rotation would pay for failed BPC and planetary data lookups before reaching the EPA fallback. The overhead is negligible in practice, but the hard-coded approach is simpler and consistent with the ECLIPJ2000 precedent.

Keep ICRF as J2000 alias (status quo)

SPICE treats ICRF = J2000. Sufficient for most applications (~0.7m at Earth surface). But ANISE already provides higher-fidelity Earth orientation (ITRF93 via BPC) and Moon orientation (MOON_PA, MOON_ME) than SPICE's defaults. Adding the ~23 mas frame bias is consistent with ANISE's philosophy of improving on SPICE where the data supports it.

References

• [IERS Conventions (2010), Chapter 5](https://iers-conventions.obspm.fr/content/tn36.pdf) — IERS Technical Note No. 36, Petit & Luzum (eds.)
• [USNO Circular 179](https://aa.usno.navy.mil/downloads/Circular_179.pdf) (Kaplan, 2005) — The IAU Resolutions on Astronomical Reference Systems, Time Scales, and Earth Rotation Models
• [Hilton et al. (2006)](https://link.springer.com/article/10.1007/s10569-006-0001-2) — Report of the International Astronomical Union Division I Working Group on Precession and the Ecliptic, Celest. Mech. Dyn. Astr., 94, 351–367
• [Chapront et al. (2002)](https://www.aanda.org/articles/aa/abs/2002/20/aa2201/aa2201.html) — A New Determination of Lunar Orbital Parameters, Precession Constant and Tidal Acceleration from LLR Measurements, Astron. Astrophys., 387, 700–709
• [Mathews, Herring & Buffett (2002)](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2001JB000390) — Modeling of Nutation and Precession: New Nutation Series for Nonrigid Earth, J. Geophys. Res., 107, B4 (MHB2000)
• SOFA iauBi00 / eraBi00 (ERFA open-source equivalent) — IAU Standards of Fundamental Astronomy
• [NAIF Frames Required Reading](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/frames.html) — Built-in Inertial Reference Frames (Appendix)

[ChristopherRabotin]

This is excellent Juno, thank you for the proposal and the write-up. This is precisely how I would envision implementing it. Would you like to implement this? I'll soon start on version 0.10, which will also implement #675 , and provide support for CK kernels.

[translunar]

Paging @ChristopherRabotin :) mind having a look and letting me know what you think?

[ChristopherRabotin]

@translunar , yes, I've had a look, this is excellent and exactly how I would implement it: #686 (comment) .

[translunar]

Oh that's so weird, I didn't get your first reply. Yes, I'll get started on it.

[c-type]

This looks very good Juno. I like the clear distinction between the reference system (ICRS) and the frames.

[c-type]

Juno, one more recommendation for the tests that check against expected offsets of the new frames. The offset is only expected for vectors perpendicular to the axis of rotation defined by the bias direction cosine matrix. The offset should actually be zero along the axis of rotation. You could determine the direction of the axis of rotation and make the tests more precise that way.

[ChristopherRabotin]

The last commit needed for this feature was the exporting of the constants, and is now in #706 as a cherry-pick from #692 to keep original authorship.

Thanks for your work on this, Juno!