celestial_core/cio/
locator.rs

1//! CIO locator (s) for the IAU 2006/2000A precession-nutation model.
2//!
3//! The CIO locator `s` positions the Celestial Intermediate Origin on the CIP equator.
4//! It's the arc length from the GCRS x-axis intersection to the CIO, measured along
5//! the CIP equator. This small angle (microarcseconds) completes the transformation
6//! from GCRS to CIRS coordinates.
7//!
8//! # The transformation chain
9//!
10//! GCRS -> CIRS requires three pieces:
11//! 1. CIP coordinates (X, Y) - where the pole is
12//! 2. CIO locator (s) - where the origin is
13//! 3. Earth Rotation Angle - how much Earth has rotated
14//!
15//! This module provides piece #2.
16//!
17//! # When to use this
18//!
19//! You need the CIO locator when:
20//! - Building the GCRS-to-CIRS rotation matrix via `gcrs_to_cirs_matrix(x, y, s)`
21//! - Computing equation of the origins (difference between CEO-based and equinox-based sidereal time)
22//! - Implementing IAU 2000/2006 compliant coordinate transformations
23//!
24//! For most uses, [`CioSolution::calculate`](super::CioSolution::calculate) handles this automatically.
25//!
26//! # Algorithm
27//!
28//! Uses the IAU 2006/2000A series expansion with 66 periodic terms across 5 polynomial
29//! orders, plus a polynomial part. The full expression is:
30//!
31//! ```text
32//! s = series(t) - X*Y/2
33//! ```
34//!
35//! where `t` is TT centuries from J2000.0. The `-X*Y/2` term accounts for the
36//! frame rotation induced by the CIP motion.
37//!
38//! # References
39//!
40//! - Capitaine et al. (2003), A&A 400, 1145-1154
41//! - IERS Conventions (2010), Chapter 5
42//! - SOFA library: `iauS06` function
43
44use crate::constants::{ARCSEC_TO_RAD, TWOPI};
45use crate::errors::{AstroError, AstroResult};
46
47/// Computes the CIO locator angle `s` for a given epoch.
48///
49/// The locator is model-dependent. Currently only IAU 2006A is implemented,
50/// which uses the IAU 2006 precession with IAU 2000A nutation.
51///
52/// # Example
53///
54/// ```
55/// use celestial_core::cio::CioLocator;
56///
57/// // Compute s for J2000.0 + 0.5 centuries (year ~2050)
58/// let locator = CioLocator::iau2006a(0.5);
59///
60/// // X, Y from CIP coordinates (typically from precession-nutation matrix)
61/// let x = 1.0e-7;  // radians
62/// let y = 2.0e-7;  // radians
63///
64/// let s = locator.calculate(x, y).unwrap();
65/// // s is in radians, typically on the order of 10^-8
66/// ```
67#[derive(Debug, Clone)]
68pub struct CioLocator {
69    tt_centuries: f64,
70    model: CioModel,
71}
72
73#[derive(Clone, Copy)]
74struct SeriesTerm {
75    coeffs: [i8; 8],
76    sine: f64,
77    cosine: f64,
78}
79
80#[allow(clippy::excessive_precision)]
81const SP: [f64; 6] = [
82    94.00e-6,
83    3808.65e-6,
84    -122.68e-6,
85    -72574.11e-6,
86    27.98e-6,
87    15.62e-6,
88];
89
90#[allow(clippy::excessive_precision)]
91const S0: [SeriesTerm; 33] = [
92    SeriesTerm {
93        coeffs: [0, 0, 0, 0, 1, 0, 0, 0],
94        sine: -2640.73e-6,
95        cosine: 0.39e-6,
96    },
97    SeriesTerm {
98        coeffs: [0, 0, 0, 0, 2, 0, 0, 0],
99        sine: -63.53e-6,
100        cosine: 0.02e-6,
101    },
102    SeriesTerm {
103        coeffs: [0, 0, 2, -2, 3, 0, 0, 0],
104        sine: -11.75e-6,
105        cosine: -0.01e-6,
106    },
107    SeriesTerm {
108        coeffs: [0, 0, 2, -2, 1, 0, 0, 0],
109        sine: -11.21e-6,
110        cosine: -0.01e-6,
111    },
112    SeriesTerm {
113        coeffs: [0, 0, 2, -2, 2, 0, 0, 0],
114        sine: 4.57e-6,
115        cosine: 0.00e-6,
116    },
117    SeriesTerm {
118        coeffs: [0, 0, 2, 0, 3, 0, 0, 0],
119        sine: -2.02e-6,
120        cosine: 0.00e-6,
121    },
122    SeriesTerm {
123        coeffs: [0, 0, 2, 0, 1, 0, 0, 0],
124        sine: -1.98e-6,
125        cosine: 0.00e-6,
126    },
127    SeriesTerm {
128        coeffs: [0, 0, 0, 0, 3, 0, 0, 0],
129        sine: 1.72e-6,
130        cosine: 0.00e-6,
131    },
132    SeriesTerm {
133        coeffs: [0, 1, 0, 0, 1, 0, 0, 0],
134        sine: 1.41e-6,
135        cosine: 0.01e-6,
136    },
137    SeriesTerm {
138        coeffs: [0, 1, 0, 0, -1, 0, 0, 0],
139        sine: 1.26e-6,
140        cosine: 0.01e-6,
141    },
142    SeriesTerm {
143        coeffs: [1, 0, 0, 0, -1, 0, 0, 0],
144        sine: 0.63e-6,
145        cosine: 0.00e-6,
146    },
147    SeriesTerm {
148        coeffs: [1, 0, 0, 0, 1, 0, 0, 0],
149        sine: 0.63e-6,
150        cosine: 0.00e-6,
151    },
152    SeriesTerm {
153        coeffs: [0, 1, 2, -2, 3, 0, 0, 0],
154        sine: -0.46e-6,
155        cosine: 0.00e-6,
156    },
157    SeriesTerm {
158        coeffs: [0, 1, 2, -2, 1, 0, 0, 0],
159        sine: -0.45e-6,
160        cosine: 0.00e-6,
161    },
162    SeriesTerm {
163        coeffs: [0, 0, 4, -4, 4, 0, 0, 0],
164        sine: -0.36e-6,
165        cosine: 0.00e-6,
166    },
167    SeriesTerm {
168        coeffs: [0, 0, 1, -1, 1, -8, 12, 0],
169        sine: 0.24e-6,
170        cosine: 0.12e-6,
171    },
172    SeriesTerm {
173        coeffs: [0, 0, 2, 0, 0, 0, 0, 0],
174        sine: -0.32e-6,
175        cosine: 0.00e-6,
176    },
177    SeriesTerm {
178        coeffs: [0, 0, 2, 0, 2, 0, 0, 0],
179        sine: -0.28e-6,
180        cosine: 0.00e-6,
181    },
182    SeriesTerm {
183        coeffs: [1, 0, 2, 0, 3, 0, 0, 0],
184        sine: -0.27e-6,
185        cosine: 0.00e-6,
186    },
187    SeriesTerm {
188        coeffs: [1, 0, 2, 0, 1, 0, 0, 0],
189        sine: -0.26e-6,
190        cosine: 0.00e-6,
191    },
192    SeriesTerm {
193        coeffs: [0, 0, 2, -2, 0, 0, 0, 0],
194        sine: 0.21e-6,
195        cosine: 0.00e-6,
196    },
197    SeriesTerm {
198        coeffs: [0, 1, -2, 2, -3, 0, 0, 0],
199        sine: -0.19e-6,
200        cosine: 0.00e-6,
201    },
202    SeriesTerm {
203        coeffs: [0, 1, -2, 2, -1, 0, 0, 0],
204        sine: -0.18e-6,
205        cosine: 0.00e-6,
206    },
207    SeriesTerm {
208        coeffs: [0, 0, 0, 0, 0, 8, -13, -1],
209        sine: 0.10e-6,
210        cosine: -0.05e-6,
211    },
212    SeriesTerm {
213        coeffs: [0, 0, 0, 2, 0, 0, 0, 0],
214        sine: -0.15e-6,
215        cosine: 0.00e-6,
216    },
217    SeriesTerm {
218        coeffs: [2, 0, -2, 0, -1, 0, 0, 0],
219        sine: 0.14e-6,
220        cosine: 0.00e-6,
221    },
222    SeriesTerm {
223        coeffs: [0, 1, 2, -2, 2, 0, 0, 0],
224        sine: 0.14e-6,
225        cosine: 0.00e-6,
226    },
227    SeriesTerm {
228        coeffs: [1, 0, 0, -2, 1, 0, 0, 0],
229        sine: -0.14e-6,
230        cosine: 0.00e-6,
231    },
232    SeriesTerm {
233        coeffs: [1, 0, 0, -2, -1, 0, 0, 0],
234        sine: -0.14e-6,
235        cosine: 0.00e-6,
236    },
237    SeriesTerm {
238        coeffs: [0, 0, 4, -2, 4, 0, 0, 0],
239        sine: -0.13e-6,
240        cosine: 0.00e-6,
241    },
242    SeriesTerm {
243        coeffs: [0, 0, 2, -2, 4, 0, 0, 0],
244        sine: 0.11e-6,
245        cosine: 0.00e-6,
246    },
247    SeriesTerm {
248        coeffs: [1, 0, -2, 0, -3, 0, 0, 0],
249        sine: -0.11e-6,
250        cosine: 0.00e-6,
251    },
252    SeriesTerm {
253        coeffs: [1, 0, -2, 0, -1, 0, 0, 0],
254        sine: -0.11e-6,
255        cosine: 0.00e-6,
256    },
257];
258
259#[allow(clippy::excessive_precision)]
260const S1: [SeriesTerm; 3] = [
261    SeriesTerm {
262        coeffs: [0, 0, 0, 0, 2, 0, 0, 0],
263        sine: -0.07e-6,
264        cosine: 3.57e-6,
265    },
266    SeriesTerm {
267        coeffs: [0, 0, 0, 0, 1, 0, 0, 0],
268        sine: 1.73e-6,
269        cosine: -0.03e-6,
270    },
271    SeriesTerm {
272        coeffs: [0, 0, 2, -2, 3, 0, 0, 0],
273        sine: 0.00e-6,
274        cosine: 0.48e-6,
275    },
276];
277
278#[allow(clippy::excessive_precision)]
279const S2: [SeriesTerm; 25] = [
280    SeriesTerm {
281        coeffs: [0, 0, 0, 0, 1, 0, 0, 0],
282        sine: 743.52e-6,
283        cosine: -0.17e-6,
284    },
285    SeriesTerm {
286        coeffs: [0, 0, 2, -2, 2, 0, 0, 0],
287        sine: 56.91e-6,
288        cosine: 0.06e-6,
289    },
290    SeriesTerm {
291        coeffs: [0, 0, 2, 0, 2, 0, 0, 0],
292        sine: 9.84e-6,
293        cosine: -0.01e-6,
294    },
295    SeriesTerm {
296        coeffs: [0, 0, 0, 0, 2, 0, 0, 0],
297        sine: -8.85e-6,
298        cosine: 0.01e-6,
299    },
300    SeriesTerm {
301        coeffs: [0, 1, 0, 0, 0, 0, 0, 0],
302        sine: -6.38e-6,
303        cosine: -0.05e-6,
304    },
305    SeriesTerm {
306        coeffs: [1, 0, 0, 0, 0, 0, 0, 0],
307        sine: -3.07e-6,
308        cosine: 0.00e-6,
309    },
310    SeriesTerm {
311        coeffs: [0, 1, 2, -2, 2, 0, 0, 0],
312        sine: 2.23e-6,
313        cosine: 0.00e-6,
314    },
315    SeriesTerm {
316        coeffs: [0, 0, 2, 0, 1, 0, 0, 0],
317        sine: 1.67e-6,
318        cosine: 0.00e-6,
319    },
320    SeriesTerm {
321        coeffs: [1, 0, 2, 0, 2, 0, 0, 0],
322        sine: 1.30e-6,
323        cosine: 0.00e-6,
324    },
325    SeriesTerm {
326        coeffs: [0, 1, -2, 2, -2, 0, 0, 0],
327        sine: 0.93e-6,
328        cosine: 0.00e-6,
329    },
330    SeriesTerm {
331        coeffs: [1, 0, 0, -2, 0, 0, 0, 0],
332        sine: 0.68e-6,
333        cosine: 0.00e-6,
334    },
335    SeriesTerm {
336        coeffs: [0, 0, 2, -2, 1, 0, 0, 0],
337        sine: -0.55e-6,
338        cosine: 0.00e-6,
339    },
340    SeriesTerm {
341        coeffs: [1, 0, -2, 0, -2, 0, 0, 0],
342        sine: 0.53e-6,
343        cosine: 0.00e-6,
344    },
345    SeriesTerm {
346        coeffs: [0, 0, 0, 2, 0, 0, 0, 0],
347        sine: -0.27e-6,
348        cosine: 0.00e-6,
349    },
350    SeriesTerm {
351        coeffs: [1, 0, 0, 0, 1, 0, 0, 0],
352        sine: -0.27e-6,
353        cosine: 0.00e-6,
354    },
355    SeriesTerm {
356        coeffs: [1, 0, -2, -2, -2, 0, 0, 0],
357        sine: -0.26e-6,
358        cosine: 0.00e-6,
359    },
360    SeriesTerm {
361        coeffs: [1, 0, 0, 0, -1, 0, 0, 0],
362        sine: -0.25e-6,
363        cosine: 0.00e-6,
364    },
365    SeriesTerm {
366        coeffs: [1, 0, 2, 0, 1, 0, 0, 0],
367        sine: 0.22e-6,
368        cosine: 0.00e-6,
369    },
370    SeriesTerm {
371        coeffs: [2, 0, 0, -2, 0, 0, 0, 0],
372        sine: -0.21e-6,
373        cosine: 0.00e-6,
374    },
375    SeriesTerm {
376        coeffs: [2, 0, -2, 0, -1, 0, 0, 0],
377        sine: 0.20e-6,
378        cosine: 0.00e-6,
379    },
380    SeriesTerm {
381        coeffs: [0, 0, 2, 2, 2, 0, 0, 0],
382        sine: 0.17e-6,
383        cosine: 0.00e-6,
384    },
385    SeriesTerm {
386        coeffs: [2, 0, 2, 0, 2, 0, 0, 0],
387        sine: 0.13e-6,
388        cosine: 0.00e-6,
389    },
390    SeriesTerm {
391        coeffs: [2, 0, 0, 0, 0, 0, 0, 0],
392        sine: -0.13e-6,
393        cosine: 0.00e-6,
394    },
395    SeriesTerm {
396        coeffs: [1, 0, 2, -2, 2, 0, 0, 0],
397        sine: -0.12e-6,
398        cosine: 0.00e-6,
399    },
400    SeriesTerm {
401        coeffs: [0, 0, 2, 0, 0, 0, 0, 0],
402        sine: -0.11e-6,
403        cosine: 0.00e-6,
404    },
405];
406
407#[allow(clippy::excessive_precision)]
408const S3: [SeriesTerm; 4] = [
409    SeriesTerm {
410        coeffs: [0, 0, 0, 0, 1, 0, 0, 0],
411        sine: 0.30e-6,
412        cosine: -23.42e-6,
413    },
414    SeriesTerm {
415        coeffs: [0, 0, 2, -2, 2, 0, 0, 0],
416        sine: -0.03e-6,
417        cosine: -1.46e-6,
418    },
419    SeriesTerm {
420        coeffs: [0, 0, 2, 0, 2, 0, 0, 0],
421        sine: -0.01e-6,
422        cosine: -0.25e-6,
423    },
424    SeriesTerm {
425        coeffs: [0, 0, 0, 0, 2, 0, 0, 0],
426        sine: 0.00e-6,
427        cosine: 0.23e-6,
428    },
429];
430
431#[allow(clippy::excessive_precision)]
432const S4: [SeriesTerm; 1] = [SeriesTerm {
433    coeffs: [0, 0, 0, 0, 1, 0, 0, 0],
434    sine: -0.26e-6,
435    cosine: -0.01e-6,
436}];
437
438#[inline]
439fn sum_terms(w: &mut f64, terms: &[SeriesTerm], fa: &[f64; 8]) {
440    for term in terms.iter().rev() {
441        let mut arg = 0.0;
442        for (i, item) in fa.iter().enumerate() {
443            arg += f64::from(term.coeffs[i]) * item;
444        }
445        *w += term.sine * libm::sin(arg) + term.cosine * libm::cos(arg);
446    }
447}
448
449fn cio_series_s_plus_xy_half(t: f64) -> f64 {
450    let fa = fundamental_arguments(t);
451
452    let mut w0 = SP[0];
453    sum_terms(&mut w0, &S0, &fa);
454
455    let mut w1 = SP[1];
456    sum_terms(&mut w1, &S1, &fa);
457
458    let mut w2 = SP[2];
459    sum_terms(&mut w2, &S2, &fa);
460
461    let mut w3 = SP[3];
462    sum_terms(&mut w3, &S3, &fa);
463
464    let mut w4 = SP[4];
465    sum_terms(&mut w4, &S4, &fa);
466
467    let w5 = SP[5];
468
469    let series_arcsec = w0 + (w1 + (w2 + (w3 + (w4 + w5 * t) * t) * t) * t) * t;
470    series_arcsec * ARCSEC_TO_RAD
471}
472
473fn fundamental_arguments(t: f64) -> [f64; 8] {
474    [
475        mean_anomaly_moon(t),
476        mean_anomaly_sun(t),
477        mean_longitude_moon_minus_node(t),
478        mean_elongation_moon_sun(t),
479        mean_longitude_ascending_node_moon(t),
480        mean_longitude_venus(t),
481        mean_longitude_earth(t),
482        general_precession_longitude(t),
483    ]
484}
485
486#[allow(clippy::excessive_precision)]
487fn mean_anomaly_moon(t: f64) -> f64 {
488    (485868.249036 + t * (1717915923.2178 + t * (31.8792 + t * (0.051635 + t * (-0.00024470)))))
489        % crate::constants::CIRCULAR_ARCSECONDS
490        * ARCSEC_TO_RAD
491}
492
493#[allow(clippy::excessive_precision)]
494fn mean_anomaly_sun(t: f64) -> f64 {
495    (1287104.793048 + t * (129596581.0481 + t * (-0.5532 + t * (0.000136 + t * (-0.00001149)))))
496        % crate::constants::CIRCULAR_ARCSECONDS
497        * ARCSEC_TO_RAD
498}
499
500#[allow(clippy::excessive_precision)]
501fn mean_longitude_moon_minus_node(t: f64) -> f64 {
502    (335779.526232 + t * (1739527262.8478 + t * (-12.7512 + t * (-0.001037 + t * (0.00000417)))))
503        % crate::constants::CIRCULAR_ARCSECONDS
504        * ARCSEC_TO_RAD
505}
506
507#[allow(clippy::excessive_precision)]
508fn mean_elongation_moon_sun(t: f64) -> f64 {
509    (1072260.703692 + t * (1602961601.2090 + t * (-6.3706 + t * (0.006593 + t * (-0.00003169)))))
510        % crate::constants::CIRCULAR_ARCSECONDS
511        * ARCSEC_TO_RAD
512}
513
514#[allow(clippy::excessive_precision)]
515fn mean_longitude_ascending_node_moon(t: f64) -> f64 {
516    (450160.398036 + t * (-6962890.5431 + t * (7.4722 + t * (0.007702 + t * (-0.00005939)))))
517        % crate::constants::CIRCULAR_ARCSECONDS
518        * ARCSEC_TO_RAD
519}
520
521#[allow(clippy::excessive_precision)]
522fn mean_longitude_venus(t: f64) -> f64 {
523    (3.176146697 + 1021.3285546211 * t) % TWOPI
524}
525
526#[allow(clippy::excessive_precision)]
527fn mean_longitude_earth(t: f64) -> f64 {
528    (1.753470314 + 628.3075849991 * t) % TWOPI
529}
530
531#[allow(clippy::excessive_precision)]
532fn general_precession_longitude(t: f64) -> f64 {
533    (0.024381750 + 0.00000538691 * t) * t
534}
535
536/// The precession-nutation model used for CIO locator computation.
537#[derive(Debug, Clone, Copy, PartialEq)]
538enum CioModel {
539    /// IAU 2006 precession with IAU 2000A nutation (66 terms).
540    IAU2006A,
541}
542
543impl CioLocator {
544    /// Creates a CIO locator using the IAU 2006/2000A model.
545    ///
546    /// # Parameters
547    ///
548    /// * `tt_centuries` - TT (Terrestrial Time) as Julian centuries from J2000.0.
549    ///   Computed as `(JD_TT - 2451545.0) / celestial_core::constants::DAYS_PER_JULIAN_CENTURY`.
550    pub fn iau2006a(tt_centuries: f64) -> Self {
551        Self {
552            tt_centuries,
553            model: CioModel::IAU2006A,
554        }
555    }
556
557    /// Computes the CIO locator `s` given the CIP coordinates.
558    ///
559    /// # Parameters
560    ///
561    /// * `x` - CIP X coordinate in radians (from NPB matrix element [2][0])
562    /// * `y` - CIP Y coordinate in radians (from NPB matrix element [2][1])
563    ///
564    /// # Returns
565    ///
566    /// The CIO locator `s` in radians.
567    ///
568    /// # Errors
569    ///
570    /// Returns an error if the epoch is more than 20 centuries from J2000.0,
571    /// where the series expansion becomes unreliable.
572    pub fn calculate(&self, x: f64, y: f64) -> AstroResult<f64> {
573        match self.model {
574            CioModel::IAU2006A => self.calculate_iau2006a(x, y),
575        }
576    }
577
578    fn calculate_iau2006a(&self, x: f64, y: f64) -> AstroResult<f64> {
579        let t = self.tt_centuries;
580
581        if t.abs() > 20.0 {
582            return Err(AstroError::math_error(
583                "CIO locator calculation",
584                crate::errors::MathErrorKind::InvalidInput,
585                &format!(
586                    "Time too far from J2000.0 for CIO locator: {:.1} centuries",
587                    t
588                ),
589            ));
590        }
591
592        let s_series = cio_series_s_plus_xy_half(t);
593        let s = s_series - 0.5 * x * y;
594
595        Ok(s)
596    }
597}
598
599#[cfg(test)]
600mod tests {
601    use super::*;
602
603    #[test]
604    fn test_cio_locator_at_j2000() {
605        let locator = CioLocator::iau2006a(0.0);
606        let s = locator.calculate(0.0, 0.0).unwrap();
607
608        assert!(
609            s.abs() < 1e-6,
610            "CIO locator at J2000.0 should be small: {}",
611            s
612        );
613    }
614
615    #[test]
616    fn test_cio_locator_time_dependency() {
617        let locator_past = CioLocator::iau2006a(-1.0);
618        let locator_future = CioLocator::iau2006a(1.0);
619
620        let s_past = locator_past.calculate(0.0, 0.0).unwrap();
621        let s_future = locator_future.calculate(0.0, 0.0).unwrap();
622
623        assert_ne!(s_past, s_future);
624        assert!(
625            (s_future - s_past).abs() > 1e-8,
626            "CIO locator should show time dependence"
627        );
628    }
629
630    #[test]
631    fn test_cio_locator_cip_dependency() {
632        let locator = CioLocator::iau2006a(0.0);
633
634        let s_zero = locator.calculate(0.0, 0.0).unwrap();
635        let s_offset = locator.calculate(1e-6, 1e-6).unwrap();
636
637        assert_ne!(s_zero, s_offset);
638    }
639
640    #[test]
641    fn test_extreme_time_validation() {
642        let locator = CioLocator::iau2006a(25.0);
643        let result = locator.calculate(0.0, 0.0);
644
645        assert!(result.is_err());
646    }
647}
celestial_core/cio/locator.rs

celestial_core/cio/
locator.rs
