gcrs_to_cirs_matrix

fn gcrs_to_cirs_matrix(x: super::Rad, y: super::Rad, s: super::Rad) -> nalgebra::Matrix3<f64>

Assemble the 3×3 celestial-to-intermediate rotation matrix C from the CIP coordinates (x, y) and CIO locator s.

The returned matrix transforms a column vector in the GCRS to the same column vector in the CIRS:

v_cirs = C · v_gcrs

Algorithm

Equivalent to SOFA iauC2ixys / ERFA c2ixys. Following IERS Conventions 2010 TN36 Eq. (5.6) / (5.10), the CIP position (x, y) defines spherical angles (E, d) via

x = sin(d) · cos(E)
y = sin(d) · sin(E)

from which we recover

E = atan2(y, x)          (when x² + y² > 0, else 0)
d = atan(sqrt(r²) / sqrt(1 − r²))        with r² = x² + y²

and the celestial-to-intermediate matrix is

C = R_z(−(E + s)) · R_y(d) · R_z(E)

where R_y, R_z are right-handed rotations about the y- and z-axes with the SOFA sign convention (positive angle rotates a vector clockwise when viewed from the positive axis).

Numerical behaviour

The r² < 1 guard is implicit: for any realistic CIP offset |x|, |y| < 10⁻⁴ rad so 1 − r² > 0.999…. The atan(sqrt(r²/(1−r²))) form matches SOFA and avoids the half-angle truncation that a naive d = asin(sqrt(r²)) would introduce at sub-µas scales.