rfa/prec_nut/
nut80.rs

1use crate::rfam::*;
2use crate::vector_matrix::angle_ops::anpm::*;
3use crate::utils::*;
4///  Nutation, IAU 1980 model.
5///
6///  Given:
7///   * date1,date2  TT as a 2-part Julian Date (Note 1)
8///
9///  Returned:
10///   * dpsi    nutation in longitude (radians)
11///   * deps    nutation in obliquity (radians)
12///
13/// # Notes:
14///
15/// 1) The TT date date1+date2 is a Julian Date, apportioned in any
16///    convenient way between the two arguments.  For example,
17///    JD(TT)=2450123.7 could be expressed in any of these ways,
18///    among others:
19///
20///    |    date1    |      date2   |                      |
21///    |-------------|--------------|----------------------|
22///    |2450123.7    |       0.0    |  (JD method)         |
23///    |2451545.0    |   -1421.3    |  (J2000 method)      |
24///    |2400000.5    |   50123.2    |  (MJD method)        |
25///    |2450123.5    |       0.2    | (date & time method) |
26///
27///    The JD method is the most natural and convenient to use in
28///    cases where the loss of several decimal digits of resolution
29///    is acceptable.  The J2000 method is best matched to the way
30///    the argument is handled internally and will deliver the
31///    optimum resolution.  The MJD method and the date & time methods
32///    are both good compromises between resolution and convenience.
33///
34///  2) The nutation components are with respect to the ecliptic of
35///     date.
36///
37/// # Called:
38///   * anpm      normalize angle into range +/- pi
39///
40/// # Reference:
41///   * Explanatory Supplement to the Astronomical Almanac,
42///     P. Kenneth Seidelmann (ed), University Science Books (1992),
43///     Section 3.222 (p111).
44///
45///  This revision:  2021 May 11
46pub fn nut80(date1: f64, date2: f64, dpsi: &mut f64, deps: &mut f64)
47{
48
49// Units of 0.1 milliarcsecond to radians 
50   const U2R: f64 = URSA_DAS2R / 1e4;
51
52// ------------------------------------------------ 
53// Table of multiples of arguments and coefficients 
54// ------------------------------------------------ 
55
56// The units for the sine and cosine coefficients are 0.1 mas and 
57// the same per Julian century 
58struct Xls {
59   nl: i32, nlp: i32, nf: i32, nd: i32, nom: i32, /* coefficients of l,l',F,D,Om */
60   sp: f64, spt: f64,        /*  longitude sine, 1 and t coefficients */
61   ce: f64, cet: f64        /* obliquity cosine, 1 and t coefficients  */
62}
63impl Xls {
64   pub const fn new(nl: i32,nlp: i32,nf: i32,nd: i32,nom: i32,sp: f64, spt: f64,ce: f64,cet: f64)->Self{
65       Xls { nl: nl, nlp: nlp, nf: nf, nd: nd, nom: nom, sp: sp, spt: spt, ce: ce, cet: cet }
66   }
67}
68const X:[Xls; 106] = [
69   // 1-10 
70      Xls::new(0,  0,  0,  0,  1, -171996.0, -174.2,  92025.0,    8.9 ),
71      Xls::new(  0,  0,  0,  0,  2,    2062.0,    0.2,   -895.0,    0.5 ),
72      Xls::new( -2,  0,  2,  0,  1,      46.0,    0.0,    -24.0,    0.0 ),
73      Xls::new(  2,  0, -2,  0,  0,      11.0,    0.0,      0.0,    0.0 ),
74      Xls::new( -2,  0,  2,  0,  2,      -3.0,    0.0,      1.0,    0.0 ),
75      Xls::new(  1, -1,  0, -1,  0,      -3.0,    0.0,      0.0,    0.0 ),
76      Xls::new(  0, -2,  2, -2,  1,      -2.0,    0.0,      1.0,    0.0 ),
77      Xls::new(  2,  0, -2,  0,  1,       1.0,    0.0,      0.0,    0.0 ),
78      Xls::new(  0,  0,  2, -2,  2,  -13187.0,   -1.6,   5736.0,   -3.1 ),
79      Xls::new(  0,  1,  0,  0,  0,    1426.0,   -3.4,     54.0,   -0.1 ),
80
81   // 11-20 
82      Xls::new(  0,  1,  2, -2,  2,    -517.0,    1.2,    224.0,   -0.6 ),
83      Xls::new(  0, -1,  2, -2,  2,     217.0,   -0.5,    -95.0,    0.3 ),
84      Xls::new(  0,  0,  2, -2,  1,     129.0,    0.1,    -70.0,    0.0 ),
85      Xls::new(  2,  0,  0, -2,  0,      48.0,    0.0,      1.0,    0.0 ),
86      Xls::new(  0,  0,  2, -2,  0,     -22.0,    0.0,      0.0,    0.0 ),
87      Xls::new(  0,  2,  0,  0,  0,      17.0,   -0.1,      0.0,    0.0 ),
88      Xls::new(  0,  1,  0,  0,  1,     -15.0,    0.0,      9.0,    0.0 ),
89      Xls::new(  0,  2,  2, -2,  2,     -16.0,    0.1,      7.0,    0.0 ),
90      Xls::new(  0, -1,  0,  0,  1,     -12.0,    0.0,      6.0,    0.0 ),
91      Xls::new( -2,  0,  0,  2,  1,      -6.0,    0.0,      3.0,    0.0 ),
92
93   // 21-30 
94      Xls::new(  0, -1,  2, -2,  1,      -5.0,    0.0,      3.0,    0.0 ),
95      Xls::new(  2,  0,  0, -2,  1,       4.0,    0.0,     -2.0,    0.0 ),
96      Xls::new(  0,  1,  2, -2,  1,       4.0,    0.0,     -2.0,    0.0 ),
97      Xls::new(  1,  0,  0, -1,  0,      -4.0,    0.0,      0.0,    0.0 ),
98      Xls::new(  2,  1,  0, -2,  0,       1.0,    0.0,      0.0,    0.0 ),
99      Xls::new(  0,  0, -2,  2,  1,       1.0,    0.0,      0.0,    0.0 ),
100      Xls::new(  0,  1, -2,  2,  0,      -1.0,    0.0,      0.0,    0.0 ),
101      Xls::new(  0,  1,  0,  0,  2,       1.0,    0.0,      0.0,    0.0 ),
102      Xls::new( -1,  0,  0,  1,  1,       1.0,    0.0,      0.0,    0.0 ),
103      Xls::new(  0,  1,  2, -2,  0,      -1.0,    0.0,      0.0,    0.0 ),
104
105   // 31-40 
106      Xls::new(  0,  0,  2,  0,  2,   -2274.0,   -0.2,    977.0,   -0.5 ),
107      Xls::new(  1,  0,  0,  0,  0,     712.0,    0.1,     -7.0,    0.0 ),
108      Xls::new(  0,  0,  2,  0,  1,    -386.0,   -0.4,    200.0,    0.0 ),
109      Xls::new(  1,  0,  2,  0,  2,    -301.0,    0.0,    129.0,   -0.1 ),
110      Xls::new(  1,  0,  0, -2,  0,    -158.0,    0.0,     -1.0,    0.0 ),
111      Xls::new( -1,  0,  2,  0,  2,     123.0,    0.0,    -53.0,    0.0 ),
112      Xls::new(  0,  0,  0,  2,  0,      63.0,    0.0,     -2.0,    0.0 ),
113      Xls::new(  1,  0,  0,  0,  1,      63.0,    0.1,    -33.0,    0.0 ),
114      Xls::new( -1,  0,  0,  0,  1,     -58.0,   -0.1,     32.0,    0.0 ),
115      Xls::new( -1,  0,  2,  2,  2,     -59.0,    0.0,     26.0,    0.0 ),
116
117   // 41-50 
118      Xls::new(  1,  0,  2,  0,  1,     -51.0,    0.0,     27.0,    0.0 ),
119      Xls::new(  0,  0,  2,  2,  2,     -38.0,    0.0,     16.0,    0.0 ),
120      Xls::new(  2,  0,  0,  0,  0,      29.0,    0.0,     -1.0,    0.0 ),
121      Xls::new(  1,  0,  2, -2,  2,      29.0,    0.0,    -12.0,    0.0 ),
122      Xls::new(  2,  0,  2,  0,  2,     -31.0,    0.0,     13.0,    0.0 ),
123      Xls::new(  0,  0,  2,  0,  0,      26.0,    0.0,     -1.0,    0.0 ),
124      Xls::new( -1,  0,  2,  0,  1,      21.0,    0.0,    -10.0,    0.0 ),
125      Xls::new( -1,  0,  0,  2,  1,      16.0,    0.0,     -8.0,    0.0 ),
126      Xls::new(  1,  0,  0, -2,  1,     -13.0,    0.0,      7.0,    0.0 ),
127      Xls::new( -1,  0,  2,  2,  1,     -10.0,    0.0,      5.0,    0.0 ),
128
129   // 51-60 
130      Xls::new(  1,  1,  0, -2,  0,      -7.0,    0.0,      0.0,    0.0 ),
131      Xls::new(  0,  1,  2,  0,  2,       7.0,    0.0,     -3.0,    0.0 ),
132      Xls::new(  0, -1,  2,  0,  2,      -7.0,    0.0,      3.0,    0.0 ),
133      Xls::new(  1,  0,  2,  2,  2,      -8.0,    0.0,      3.0,    0.0 ),
134      Xls::new(  1,  0,  0,  2,  0,       6.0,    0.0,      0.0,    0.0 ),
135      Xls::new(  2,  0,  2, -2,  2,       6.0,    0.0,     -3.0,    0.0 ),
136      Xls::new(  0,  0,  0,  2,  1,      -6.0,    0.0,      3.0,    0.0 ),
137      Xls::new(  0,  0,  2,  2,  1,      -7.0,    0.0,      3.0,    0.0 ),
138      Xls::new(  1,  0,  2, -2,  1,       6.0,    0.0,     -3.0,    0.0 ),
139      Xls::new(  0,  0,  0, -2,  1,      -5.0,    0.0,      3.0,    0.0 ),
140
141   // 61-70 
142      Xls::new(  1, -1,  0,  0,  0,       5.0,    0.0,      0.0,    0.0 ),
143      Xls::new(  2,  0,  2,  0,  1,      -5.0,    0.0,      3.0,    0.0 ),
144      Xls::new(  0,  1,  0, -2,  0,      -4.0,    0.0,      0.0,    0.0 ),
145      Xls::new(  1,  0, -2,  0,  0,       4.0,    0.0,      0.0,    0.0 ),
146      Xls::new(  0,  0,  0,  1,  0,      -4.0,    0.0,      0.0,    0.0 ),
147      Xls::new(  1,  1,  0,  0,  0,      -3.0,    0.0,      0.0,    0.0 ),
148      Xls::new(  1,  0,  2,  0,  0,       3.0,    0.0,      0.0,    0.0 ),
149      Xls::new(  1, -1,  2,  0,  2,      -3.0,    0.0,      1.0,    0.0 ),
150      Xls::new( -1, -1,  2,  2,  2,      -3.0,    0.0,      1.0,    0.0 ),
151      Xls::new( -2,  0,  0,  0,  1,      -2.0,    0.0,      1.0,    0.0 ),
152
153   // 71-80 
154      Xls::new(  3,  0,  2,  0,  2,      -3.0,    0.0,      1.0,    0.0 ),
155      Xls::new(  0, -1,  2,  2,  2,      -3.0,    0.0,      1.0,    0.0 ),
156      Xls::new(  1,  1,  2,  0,  2,       2.0,    0.0,     -1.0,    0.0 ),
157      Xls::new( -1,  0,  2, -2,  1,      -2.0,    0.0,      1.0,    0.0 ),
158      Xls::new(  2,  0,  0,  0,  1,       2.0,    0.0,     -1.0,    0.0 ),
159      Xls::new(  1,  0,  0,  0,  2,      -2.0,    0.0,      1.0,    0.0 ),
160      Xls::new(  3,  0,  0,  0,  0,       2.0,    0.0,      0.0,    0.0 ),
161      Xls::new(  0,  0,  2,  1,  2,       2.0,    0.0,     -1.0,    0.0 ),
162      Xls::new( -1,  0,  0,  0,  2,       1.0,    0.0,     -1.0,    0.0 ),
163      Xls::new(  1,  0,  0, -4,  0,      -1.0,    0.0,      0.0,    0.0 ),
164
165   // 81-90 
166      Xls::new( -2,  0,  2,  2,  2,       1.0,    0.0,     -1.0,    0.0 ),
167      Xls::new( -1,  0,  2,  4,  2,      -2.0,    0.0,      1.0,    0.0 ),
168      Xls::new(  2,  0,  0, -4,  0,      -1.0,    0.0,      0.0,    0.0 ),
169      Xls::new(  1,  1,  2, -2,  2,       1.0,    0.0,     -1.0,    0.0 ),
170      Xls::new(  1,  0,  2,  2,  1,      -1.0,    0.0,      1.0,    0.0 ),
171      Xls::new( -2,  0,  2,  4,  2,      -1.0,    0.0,      1.0,    0.0 ),
172      Xls::new( -1,  0,  4,  0,  2,       1.0,    0.0,      0.0,    0.0 ),
173      Xls::new(  1, -1,  0, -2,  0,       1.0,    0.0,      0.0,    0.0 ),
174      Xls::new(  2,  0,  2, -2,  1,       1.0,    0.0,     -1.0,    0.0 ),
175      Xls::new(  2,  0,  2,  2,  2,      -1.0,    0.0,      0.0,    0.0 ),
176
177   // 91-100 
178      Xls::new(  1,  0,  0,  2,  1,      -1.0,    0.0,      0.0,    0.0 ),
179      Xls::new(  0,  0,  4, -2,  2,       1.0,    0.0,      0.0,    0.0 ),
180      Xls::new(  3,  0,  2, -2,  2,       1.0,    0.0,      0.0,    0.0 ),
181      Xls::new(  1,  0,  2, -2,  0,      -1.0,    0.0,      0.0,    0.0 ),
182      Xls::new(  0,  1,  2,  0,  1,       1.0,    0.0,      0.0,    0.0 ),
183      Xls::new( -1, -1,  0,  2,  1,       1.0,    0.0,      0.0,    0.0 ),
184      Xls::new(  0,  0, -2,  0,  1,      -1.0,    0.0,      0.0,    0.0 ),
185      Xls::new(  0,  0,  2, -1,  2,      -1.0,    0.0,      0.0,    0.0 ),
186      Xls::new(  0,  1,  0,  2,  0,      -1.0,    0.0,      0.0,    0.0 ),
187      Xls::new(  1,  0, -2, -2,  0,      -1.0,    0.0,      0.0,    0.0 ),
188
189   // 101-106 
190      Xls::new(  0, -1,  2,  0,  1,      -1.0,    0.0,      0.0,    0.0 ),
191      Xls::new(  1,  1,  0, -2,  1,      -1.0,    0.0,      0.0,    0.0 ),
192      Xls::new(  1,  0, -2,  2,  0,      -1.0,    0.0,      0.0,    0.0 ),
193      Xls::new(  2,  0,  0,  2,  0,       1.0,    0.0,      0.0,    0.0 ),
194      Xls::new(  0,  0,  2,  4,  2,      -1.0,    0.0,      0.0,    0.0 ),
195      Xls::new(  0,  1,  0,  1,  0,       1.0,    0.0,      0.0,    0.0 )
196   ];
197
198// ------------------------------------------------------------------ 
199
200// Interval between fundamental epoch J2000.0 and given date (JC). 
201let t = ((date1 - URSA_DJ00) + date2) / URSA_DJC;
202
203// --------------------- 
204// Fundamental arguments 
205// --------------------- 
206
207// Mean longitude of Moon minus mean longitude of Moon's perigee. 
208let el = anpm(
209        (485866.733 + (715922.633 + (31.310 + 0.064 * t) * t) * t)
210        * URSA_DAS2R + fmod(1325.0 * t, 1.0) * URSA_D2PI);
211
212// Mean longitude of Sun minus mean longitude of Sun's perigee. 
213let elp = anpm(
214         (1287099.804 + (1292581.224 + (-0.577 - 0.012 * t) * t) * t)
215         * URSA_DAS2R + fmod(99.0 * t, 1.0) * URSA_D2PI);
216
217// Mean longitude of Moon minus mean longitude of Moon's node. 
218let f = anpm(
219       (335778.877 + (295263.137 + (-13.257 + 0.011 * t) * t) * t)
220       * URSA_DAS2R + fmod(1342.0 * t, 1.0) * URSA_D2PI);
221
222// Mean elongation of Moon from Sun. 
223let d = anpm(
224       (1072261.307 + (1105601.328 + (-6.891 + 0.019 * t) * t) * t)
225       * URSA_DAS2R + fmod(1236.0 * t, 1.0) * URSA_D2PI);
226
227// Longitude of the mean ascending node of the lunar orbit on the 
228// ecliptic, measured from the mean equinox of date. 
229let om = anpm(
230        (450160.280 + (-482890.539 + (7.455 + 0.008 * t) * t) * t)
231        * URSA_DAS2R + fmod(-5.0 * t, 1.0) * URSA_D2PI);
232
233// --------------- 
234// Nutation series 
235// --------------- 
236
237// Initialize nutation components. 
238let mut dp = 0.0;
239let mut de = 0.0;
240
241// Sum the nutation terms, ending with the biggest. 
242   for x_j in X.iter() {
243
244   // Form argument for current term. 
245      let arg = x_j.nl as f64 * el
246          + x_j.nlp as f64 * elp
247          + x_j.nf as f64  * f
248          + x_j.nd as f64  * d
249          + x_j.nom as f64 * om;
250
251   // Accumulate current nutation term. 
252      let s = x_j.sp + x_j.spt * t;
253      let c = x_j.ce + x_j.cet * t;
254      if s != 0.0 {dp += s * sin(arg)};
255      if c != 0.0 {de += c * cos(arg)};
256   }
257
258// Convert results from 0.1 mas units to radians. 
259   *dpsi = dp * U2R;
260   *deps = de * U2R;
261
262// Finished. 
263
264}
rfa/prec_nut/nut80.rs

rfa/prec_nut/
nut80.rs
