erfa/time/
s06.rs

1// This Source Code Form is subject to the terms of the Mozilla Public
2// License, v. 2.0. If a copy of the MPL was not distributed with this
3// file, You can obtain one at http://mozilla.org/MPL/2.0/.
4
5use crate::{constants::*, fundamental_argument::*};
6
7/// The CIO locator s, positioning the Celestial Intermediate Origin on the
8/// equator of the Celestial Intermediate Pole, given the CIP's X,Y coordinates.
9/// Compatible with IAU 2006/2000A precession-nutation. (`eraS06`)
10///
11/// Given:
12///  * `date1`,`date2`: TT as a 2-part Julian Date (Note 1)
13///  * `x`,`y`: CIP coordinates (Note 3)
14///
15/// Returned:
16///  * the CIO locator s in radians (Note 2)
17///
18/// # Notes:
19///
20/// 1) The TT date `date1+date2` is a Julian Date, apportioned in any convenient
21///    way between the two arguments.  For example, `JD(TT)=2450123.7` could be
22///    expressed in any of these ways, among others:
23///
24///    | `date1`   | `date2` |                    |
25///    |-----------|---------|--------------------|
26///    | 2450123.7 |     0.0 | JD method          |
27///    | 2451545.0 | -1421.3 | J2000 method       |
28///    | 2400000.5 | 50123.2 | MJD method         |
29///    | 2450123.5 |     0.2 | date & time method |
30///
31///    The JD method is the most natural and convenient to use in cases where
32///    the loss of several decimal digits of resolution is acceptable.  The
33///    J2000 method is best matched to the way the argument is handled
34///    internally and will deliver the optimum resolution.  The MJD method and
35///    the date & time methods are both good compromises between resolution and
36///    convenience.
37///
38/// 2) The CIO locator s is the difference between the right ascensions of the
39///    same point in two systems:  the two systems are the GCRS and the CIP,CIO,
40///    and the point is the ascending node of the CIP equator.  The quantity s
41///    remains below 0.1 arcsecond throughout 1900-2100.
42///
43/// 3) The series used to compute s is in fact for s+XY/2, where X and Y are the
44///    x and y components of the CIP unit vector;  this series is more compact
45///    than a direct series for s would be.  This function requires X,Y to be
46///    supplied by the caller, who is responsible for providing values that are
47///    consistent with the supplied date.
48///
49/// 4) The model is consistent with the "P03" precession (Capitaine et al.
50///    2003), adopted by IAU 2006 Resolution 1, 2006, and the IAU 2000A nutation
51///    (with P03 adjustments).
52///
53/// # References:
54///
55/// * Capitaine, N., Wallace, P.T. & Chapront, J., 2003, Astron. Astrophys. 432,
56///   355
57///
58/// * McCarthy, D.D., Petit, G. (eds.) 2004, IERS Conventions (2003), IERS
59///   Technical Note No. 32, BKG
60///
61#[allow(non_snake_case)]
62pub fn S06(date1: f64, date2: f64, x: f64, y: f64) -> f64 {
63    /* Interval between fundamental epoch J2000.0 and current date (JC). */
64    let t = (date1 - ERFA_DJ00 + date2) / ERFA_DJC;
65
66    /* Fundamental Arguments (from IERS Conventions 2003) */
67    let fa: [f64; 8] = [
68        /* Mean anomaly of the Moon. */
69        l03(t),
70        /* Mean anomaly of the Sun. */
71        lp03(t),
72        /* Mean longitude of the Moon minus that of the ascending node. */
73        f03(t),
74        /* Mean elongation of the Moon from the Sun. */
75        d03(t),
76        /* Mean longitude of the ascending node of the Moon. */
77        om03(t),
78        /* Mean longitude of Venus. */
79        ve03(t),
80        /* Mean longitude of Earth. */
81        e03(t),
82        /* General precession in longitude. */
83        pa03(t),
84    ];
85
86    /* Evaluate s. */
87    let mut w0 = SP[0];
88    let mut w1 = SP[1];
89    let mut w2 = SP[2];
90    let mut w3 = SP[3];
91    let mut w4 = SP[4];
92    let w5 = SP[5];
93
94    for s0 in S0.iter().rev() {
95        let a = s0
96            .nfa
97            .iter()
98            .copied()
99            .zip(fa.iter().copied())
100            .fold(0.0, |acc, (nfa, fa)| acc + f64::from(nfa) * fa);
101        w0 += s0.s * a.sin() + s0.c * a.cos();
102    }
103    for s1 in S1.iter().rev() {
104        let a = s1
105            .nfa
106            .iter()
107            .copied()
108            .zip(fa.iter().copied())
109            .fold(0.0, |acc, (nfa, fa)| acc + f64::from(nfa) * fa);
110        w1 += s1.s * a.sin() + s1.c * a.cos();
111    }
112    for s2 in S2.iter().rev() {
113        let a = s2
114            .nfa
115            .iter()
116            .copied()
117            .zip(fa.iter().copied())
118            // .fold(0.0, |acc, (nfa, fa)| acc + f64::from(nfa) + fa);
119            .fold(0.0, |acc, (nfa, fa)| acc + f64::from(nfa) * fa);
120        w2 += s2.s * a.sin() + s2.c * a.cos();
121    }
122    for s3 in S3.iter().rev() {
123        let a = s3
124            .nfa
125            .iter()
126            .copied()
127            .zip(fa.iter().copied())
128            .fold(0.0, |acc, (nfa, fa)| acc + f64::from(nfa) * fa);
129        w3 += s3.s * a.sin() + s3.c * a.cos();
130    }
131    for s4 in S4.iter().rev() {
132        let a = s4
133            .nfa
134            .iter()
135            .copied()
136            .zip(fa.iter().copied())
137            .fold(0.0, |acc, (nfa, fa)| acc + f64::from(nfa) * fa);
138        w4 += s4.s * a.sin() + s4.c * a.cos();
139    }
140
141    (w0 + (w1 + (w2 + (w3 + (w4 + w5 * t) * t) * t) * t) * t) * ERFA_DAS2R - x * y / 2.0
142}
143
144/* Polynomial coefficients */
145const SP: [f64; 6] = [
146    94.00e-6,
147    3808.65e-6,
148    -122.68e-6,
149    -72574.11e-6,
150    27.98e-6,
151    15.62e-6,
152];
153
154/* --------------------- */
155/* The series for s+XY/2 */
156/* --------------------- */
157struct Term {
158    /// coefficients of l,l',F,D,Om,LVe,LE,pA
159    nfa: [i32; 8],
160    /// sine coefficients
161    s: f64,
162    /// cosine coefficients
163    c: f64,
164}
165
166/* Terms of order t^0 */
167const S0: [Term; 33] = [
168    Term {
169        nfa: [0, 0, 0, 0, 1, 0, 0, 0],
170        s: -2640.73e-6,
171        c: 0.39e-6,
172    },
173    Term {
174        nfa: [0, 0, 0, 0, 2, 0, 0, 0],
175        s: -63.53e-6,
176        c: 0.02e-6,
177    },
178    Term {
179        nfa: [0, 0, 2, -2, 3, 0, 0, 0],
180        s: -11.75e-6,
181        c: -0.01e-6,
182    },
183    Term {
184        nfa: [0, 0, 2, -2, 1, 0, 0, 0],
185        s: -11.21e-6,
186        c: -0.01e-6,
187    },
188    Term {
189        nfa: [0, 0, 2, -2, 2, 0, 0, 0],
190        s: 4.57e-6,
191        c: 0.00e-6,
192    },
193    Term {
194        nfa: [0, 0, 2, 0, 3, 0, 0, 0],
195        s: -2.02e-6,
196        c: 0.00e-6,
197    },
198    Term {
199        nfa: [0, 0, 2, 0, 1, 0, 0, 0],
200        s: -1.98e-6,
201        c: 0.00e-6,
202    },
203    Term {
204        nfa: [0, 0, 0, 0, 3, 0, 0, 0],
205        s: 1.72e-6,
206        c: 0.00e-6,
207    },
208    Term {
209        nfa: [0, 1, 0, 0, 1, 0, 0, 0],
210        s: 1.41e-6,
211        c: 0.01e-6,
212    },
213    Term {
214        nfa: [0, 1, 0, 0, -1, 0, 0, 0],
215        s: 1.26e-6,
216        c: 0.01e-6,
217    },
218    Term {
219        nfa: [1, 0, 0, 0, -1, 0, 0, 0],
220        s: 0.63e-6,
221        c: 0.00e-6,
222    },
223    Term {
224        nfa: [1, 0, 0, 0, 1, 0, 0, 0],
225        s: 0.63e-6,
226        c: 0.00e-6,
227    },
228    Term {
229        nfa: [0, 1, 2, -2, 3, 0, 0, 0],
230        s: -0.46e-6,
231        c: 0.00e-6,
232    },
233    Term {
234        nfa: [0, 1, 2, -2, 1, 0, 0, 0],
235        s: -0.45e-6,
236        c: 0.00e-6,
237    },
238    Term {
239        nfa: [0, 0, 4, -4, 4, 0, 0, 0],
240        s: -0.36e-6,
241        c: 0.00e-6,
242    },
243    Term {
244        nfa: [0, 0, 1, -1, 1, -8, 12, 0],
245        s: 0.24e-6,
246        c: 0.12e-6,
247    },
248    Term {
249        nfa: [0, 0, 2, 0, 0, 0, 0, 0],
250        s: -0.32e-6,
251        c: 0.00e-6,
252    },
253    Term {
254        nfa: [0, 0, 2, 0, 2, 0, 0, 0],
255        s: -0.28e-6,
256        c: 0.00e-6,
257    },
258    Term {
259        nfa: [1, 0, 2, 0, 3, 0, 0, 0],
260        s: -0.27e-6,
261        c: 0.00e-6,
262    },
263    Term {
264        nfa: [1, 0, 2, 0, 1, 0, 0, 0],
265        s: -0.26e-6,
266        c: 0.00e-6,
267    },
268    Term {
269        nfa: [0, 0, 2, -2, 0, 0, 0, 0],
270        s: 0.21e-6,
271        c: 0.00e-6,
272    },
273    Term {
274        nfa: [0, 1, -2, 2, -3, 0, 0, 0],
275        s: -0.19e-6,
276        c: 0.00e-6,
277    },
278    Term {
279        nfa: [0, 1, -2, 2, -1, 0, 0, 0],
280        s: -0.18e-6,
281        c: 0.00e-6,
282    },
283    Term {
284        nfa: [0, 0, 0, 0, 0, 8, -13, -1],
285        s: 0.10e-6,
286        c: -0.05e-6,
287    },
288    Term {
289        nfa: [0, 0, 0, 2, 0, 0, 0, 0],
290        s: -0.15e-6,
291        c: 0.00e-6,
292    },
293    Term {
294        nfa: [2, 0, -2, 0, -1, 0, 0, 0],
295        s: 0.14e-6,
296        c: 0.00e-6,
297    },
298    Term {
299        nfa: [0, 1, 2, -2, 2, 0, 0, 0],
300        s: 0.14e-6,
301        c: 0.00e-6,
302    },
303    Term {
304        nfa: [1, 0, 0, -2, 1, 0, 0, 0],
305        s: -0.14e-6,
306        c: 0.00e-6,
307    },
308    Term {
309        nfa: [1, 0, 0, -2, -1, 0, 0, 0],
310        s: -0.14e-6,
311        c: 0.00e-6,
312    },
313    Term {
314        nfa: [0, 0, 4, -2, 4, 0, 0, 0],
315        s: -0.13e-6,
316        c: 0.00e-6,
317    },
318    Term {
319        nfa: [0, 0, 2, -2, 4, 0, 0, 0],
320        s: 0.11e-6,
321        c: 0.00e-6,
322    },
323    Term {
324        nfa: [1, 0, -2, 0, -3, 0, 0, 0],
325        s: -0.11e-6,
326        c: 0.00e-6,
327    },
328    Term {
329        nfa: [1, 0, -2, 0, -1, 0, 0, 0],
330        s: -0.11e-6,
331        c: 0.00e-6,
332    },
333];
334
335/* Terms of order t^1 */
336const S1: [Term; 3] = [
337    Term {
338        nfa: [0, 0, 0, 0, 2, 0, 0, 0],
339        s: -0.07e-6,
340        c: 3.57e-6,
341    },
342    Term {
343        nfa: [0, 0, 0, 0, 1, 0, 0, 0],
344        s: 1.73e-6,
345        c: -0.03e-6,
346    },
347    Term {
348        nfa: [0, 0, 2, -2, 3, 0, 0, 0],
349        s: 0.00e-6,
350        c: 0.48e-6,
351    },
352];
353
354/* Terms of order t^2 */
355const S2: [Term; 25] = [
356    Term {
357        nfa: [0, 0, 0, 0, 1, 0, 0, 0],
358        s: 743.52e-6,
359        c: -0.17e-6,
360    },
361    Term {
362        nfa: [0, 0, 2, -2, 2, 0, 0, 0],
363        s: 56.91e-6,
364        c: 0.06e-6,
365    },
366    Term {
367        nfa: [0, 0, 2, 0, 2, 0, 0, 0],
368        s: 9.84e-6,
369        c: -0.01e-6,
370    },
371    Term {
372        nfa: [0, 0, 0, 0, 2, 0, 0, 0],
373        s: -8.85e-6,
374        c: 0.01e-6,
375    },
376    Term {
377        nfa: [0, 1, 0, 0, 0, 0, 0, 0],
378        s: -6.38e-6,
379        c: -0.05e-6,
380    },
381    Term {
382        nfa: [1, 0, 0, 0, 0, 0, 0, 0],
383        s: -3.07e-6,
384        c: 0.00e-6,
385    },
386    Term {
387        nfa: [0, 1, 2, -2, 2, 0, 0, 0],
388        s: 2.23e-6,
389        c: 0.00e-6,
390    },
391    Term {
392        nfa: [0, 0, 2, 0, 1, 0, 0, 0],
393        s: 1.67e-6,
394        c: 0.00e-6,
395    },
396    Term {
397        nfa: [1, 0, 2, 0, 2, 0, 0, 0],
398        s: 1.30e-6,
399        c: 0.00e-6,
400    },
401    Term {
402        nfa: [0, 1, -2, 2, -2, 0, 0, 0],
403        s: 0.93e-6,
404        c: 0.00e-6,
405    },
406    Term {
407        nfa: [1, 0, 0, -2, 0, 0, 0, 0],
408        s: 0.68e-6,
409        c: 0.00e-6,
410    },
411    Term {
412        nfa: [0, 0, 2, -2, 1, 0, 0, 0],
413        s: -0.55e-6,
414        c: 0.00e-6,
415    },
416    Term {
417        nfa: [1, 0, -2, 0, -2, 0, 0, 0],
418        s: 0.53e-6,
419        c: 0.00e-6,
420    },
421    Term {
422        nfa: [0, 0, 0, 2, 0, 0, 0, 0],
423        s: -0.27e-6,
424        c: 0.00e-6,
425    },
426    Term {
427        nfa: [1, 0, 0, 0, 1, 0, 0, 0],
428        s: -0.27e-6,
429        c: 0.00e-6,
430    },
431    Term {
432        nfa: [1, 0, -2, -2, -2, 0, 0, 0],
433        s: -0.26e-6,
434        c: 0.00e-6,
435    },
436    Term {
437        nfa: [1, 0, 0, 0, -1, 0, 0, 0],
438        s: -0.25e-6,
439        c: 0.00e-6,
440    },
441    Term {
442        nfa: [1, 0, 2, 0, 1, 0, 0, 0],
443        s: 0.22e-6,
444        c: 0.00e-6,
445    },
446    Term {
447        nfa: [2, 0, 0, -2, 0, 0, 0, 0],
448        s: -0.21e-6,
449        c: 0.00e-6,
450    },
451    Term {
452        nfa: [2, 0, -2, 0, -1, 0, 0, 0],
453        s: 0.20e-6,
454        c: 0.00e-6,
455    },
456    Term {
457        nfa: [0, 0, 2, 2, 2, 0, 0, 0],
458        s: 0.17e-6,
459        c: 0.00e-6,
460    },
461    Term {
462        nfa: [2, 0, 2, 0, 2, 0, 0, 0],
463        s: 0.13e-6,
464        c: 0.00e-6,
465    },
466    Term {
467        nfa: [2, 0, 0, 0, 0, 0, 0, 0],
468        s: -0.13e-6,
469        c: 0.00e-6,
470    },
471    Term {
472        nfa: [1, 0, 2, -2, 2, 0, 0, 0],
473        s: -0.12e-6,
474        c: 0.00e-6,
475    },
476    Term {
477        nfa: [0, 0, 2, 0, 0, 0, 0, 0],
478        s: -0.11e-6,
479        c: 0.00e-6,
480    },
481];
482
483/* Terms of order t^3 */
484const S3: [Term; 4] = [
485    Term {
486        nfa: [0, 0, 0, 0, 1, 0, 0, 0],
487        s: 0.30e-6,
488        c: -23.42e-6,
489    },
490    Term {
491        nfa: [0, 0, 2, -2, 2, 0, 0, 0],
492        s: -0.03e-6,
493        c: -1.46e-6,
494    },
495    Term {
496        nfa: [0, 0, 2, 0, 2, 0, 0, 0],
497        s: -0.01e-6,
498        c: -0.25e-6,
499    },
500    Term {
501        nfa: [0, 0, 0, 0, 2, 0, 0, 0],
502        s: 0.00e-6,
503        c: 0.23e-6,
504    },
505];
506
507/* Terms of order t^4 */
508const S4: [Term; 1] = [Term {
509    nfa: [0, 0, 0, 0, 1, 0, 0, 0],
510    s: -0.26e-6,
511    c: -0.01e-6,
512}];
erfa/time/s06.rs

erfa/time/
s06.rs
