source: trunk/FACT++/erfa/src/nut00b.c@ 19496

Last change on this file since 19496 was 18921, checked in by tbretz, 7 years ago
Updated to ERFA 1.4.0
File size: 17.2 KB
Line 
1#include "erfa.h"
2
3void eraNut00b(double date1, double date2, double *dpsi, double *deps)
4/*
5** - - - - - - - - - -
6** e r a N u t 0 0 b
7** - - - - - - - - - -
8**
9** Nutation, IAU 2000B model.
10**
11** Given:
12** date1,date2 double TT as a 2-part Julian Date (Note 1)
13**
14** Returned:
15** dpsi,deps double nutation, luni-solar + planetary (Note 2)
16**
17** Notes:
18**
19** 1) The TT date date1+date2 is a Julian Date, apportioned in any
20** convenient way between the two arguments. For example,
21** JD(TT)=2450123.7 could be expressed in any of these ways,
22** 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
32** cases where the loss of several decimal digits of resolution
33** is acceptable. The J2000 method is best matched to the way
34** the argument is handled internally and will deliver the
35** optimum resolution. The MJD method and the date & time methods
36** are both good compromises between resolution and convenience.
37**
38** 2) The nutation components in longitude and obliquity are in radians
39** and with respect to the equinox and ecliptic of date. The
40** obliquity at J2000.0 is assumed to be the Lieske et al. (1977)
41** value of 84381.448 arcsec. (The errors that result from using
42** this function with the IAU 2006 value of 84381.406 arcsec can be
43** neglected.)
44**
45** The nutation model consists only of luni-solar terms, but
46** includes also a fixed offset which compensates for certain long-
47** period planetary terms (Note 7).
48**
49** 3) This function is an implementation of the IAU 2000B abridged
50** nutation model formally adopted by the IAU General Assembly in
51** 2000. The function computes the MHB_2000_SHORT luni-solar
52** nutation series (Luzum 2001), but without the associated
53** corrections for the precession rate adjustments and the offset
54** between the GCRS and J2000.0 mean poles.
55**
56** 4) The full IAU 2000A (MHB2000) nutation model contains nearly 1400
57** terms. The IAU 2000B model (McCarthy & Luzum 2003) contains only
58** 77 terms, plus additional simplifications, yet still delivers
59** results of 1 mas accuracy at present epochs. This combination of
60** accuracy and size makes the IAU 2000B abridged nutation model
61** suitable for most practical applications.
62**
63** The function delivers a pole accurate to 1 mas from 1900 to 2100
64** (usually better than 1 mas, very occasionally just outside
65** 1 mas). The full IAU 2000A model, which is implemented in the
66** function eraNut00a (q.v.), delivers considerably greater accuracy
67** at current dates; however, to realize this improved accuracy,
68** corrections for the essentially unpredictable free-core-nutation
69** (FCN) must also be included.
70**
71** 5) The present function provides classical nutation. The
72** MHB_2000_SHORT algorithm, from which it is adapted, deals also
73** with (i) the offsets between the GCRS and mean poles and (ii) the
74** adjustments in longitude and obliquity due to the changed
75** precession rates. These additional functions, namely frame bias
76** and precession adjustments, are supported by the ERFA functions
77** eraBi00 and eraPr00.
78**
79** 6) The MHB_2000_SHORT algorithm also provides "total" nutations,
80** comprising the arithmetic sum of the frame bias, precession
81** adjustments, and nutation (luni-solar + planetary). These total
82** nutations can be used in combination with an existing IAU 1976
83** precession implementation, such as eraPmat76, to deliver GCRS-
84** to-true predictions of mas accuracy at current epochs. However,
85** for symmetry with the eraNut00a function (q.v. for the reasons),
86** the ERFA functions do not generate the "total nutations"
87** directly. Should they be required, they could of course easily
88** be generated by calling eraBi00, eraPr00 and the present function
89** and adding the results.
90**
91** 7) The IAU 2000B model includes "planetary bias" terms that are
92** fixed in size but compensate for long-period nutations. The
93** amplitudes quoted in McCarthy & Luzum (2003), namely
94** Dpsi = -1.5835 mas and Depsilon = +1.6339 mas, are optimized for
95** the "total nutations" method described in Note 6. The Luzum
96** (2001) values used in this ERFA implementation, namely -0.135 mas
97** and +0.388 mas, are optimized for the "rigorous" method, where
98** frame bias, precession and nutation are applied separately and in
99** that order. During the interval 1995-2050, the ERFA
100** implementation delivers a maximum error of 1.001 mas (not
101** including FCN).
102**
103** References:
104**
105** Lieske, J.H., Lederle, T., Fricke, W., Morando, B., "Expressions
106** for the precession quantities based upon the IAU /1976/ system of
107** astronomical constants", Astron.Astrophys. 58, 1-2, 1-16. (1977)
108**
109** Luzum, B., private communication, 2001 (Fortran code
110** MHB_2000_SHORT)
111**
112** McCarthy, D.D. & Luzum, B.J., "An abridged model of the
113** precession-nutation of the celestial pole", Cel.Mech.Dyn.Astron.
114** 85, 37-49 (2003)
115**
116** Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
117** Francou, G., Laskar, J., Astron.Astrophys. 282, 663-683 (1994)
118**
119** Copyright (C) 2013-2017, NumFOCUS Foundation.
120** Derived, with permission, from the SOFA library. See notes at end of file.
121*/
122{
123 double t, el, elp, f, d, om, arg, dp, de, sarg, carg,
124 dpsils, depsls, dpsipl, depspl;
125 int i;
126
127/* Units of 0.1 microarcsecond to radians */
128 static const double U2R = ERFA_DAS2R / 1e7;
129
130/* ---------------------------------------- */
131/* Fixed offsets in lieu of planetary terms */
132/* ---------------------------------------- */
133
134 static const double DPPLAN = -0.135 * ERFA_DMAS2R;
135 static const double DEPLAN = 0.388 * ERFA_DMAS2R;
136
137/* --------------------------------------------------- */
138/* Luni-solar nutation: argument and term coefficients */
139/* --------------------------------------------------- */
140
141/* The units for the sine and cosine coefficients are */
142/* 0.1 microarcsec and the same per Julian century */
143
144 static const struct {
145 int nl,nlp,nf,nd,nom; /* coefficients of l,l',F,D,Om */
146 double ps,pst,pc; /* longitude sin, t*sin, cos coefficients */
147 double ec,ect,es; /* obliquity cos, t*cos, sin coefficients */
148
149 } x[] = {
150
151 /* 1-10 */
152 { 0, 0, 0, 0,1,
153 -172064161.0, -174666.0, 33386.0, 92052331.0, 9086.0, 15377.0},
154 { 0, 0, 2,-2,2,
155 -13170906.0, -1675.0, -13696.0, 5730336.0, -3015.0, -4587.0},
156 { 0, 0, 2, 0,2,-2276413.0,-234.0, 2796.0, 978459.0,-485.0,1374.0},
157 { 0, 0, 0, 0,2,2074554.0, 207.0, -698.0,-897492.0, 470.0,-291.0},
158 { 0, 1, 0, 0,0,1475877.0,-3633.0,11817.0, 73871.0,-184.0,-1924.0},
159 { 0, 1, 2,-2,2,-516821.0, 1226.0, -524.0, 224386.0,-677.0,-174.0},
160 { 1, 0, 0, 0,0, 711159.0, 73.0, -872.0, -6750.0, 0.0, 358.0},
161 { 0, 0, 2, 0,1,-387298.0, -367.0, 380.0, 200728.0, 18.0, 318.0},
162 { 1, 0, 2, 0,2,-301461.0, -36.0, 816.0, 129025.0, -63.0, 367.0},
163 { 0,-1, 2,-2,2, 215829.0, -494.0, 111.0, -95929.0, 299.0, 132.0},
164
165 /* 11-20 */
166 { 0, 0, 2,-2,1, 128227.0, 137.0, 181.0, -68982.0, -9.0, 39.0},
167 {-1, 0, 2, 0,2, 123457.0, 11.0, 19.0, -53311.0, 32.0, -4.0},
168 {-1, 0, 0, 2,0, 156994.0, 10.0, -168.0, -1235.0, 0.0, 82.0},
169 { 1, 0, 0, 0,1, 63110.0, 63.0, 27.0, -33228.0, 0.0, -9.0},
170 {-1, 0, 0, 0,1, -57976.0, -63.0, -189.0, 31429.0, 0.0, -75.0},
171 {-1, 0, 2, 2,2, -59641.0, -11.0, 149.0, 25543.0, -11.0, 66.0},
172 { 1, 0, 2, 0,1, -51613.0, -42.0, 129.0, 26366.0, 0.0, 78.0},
173 {-2, 0, 2, 0,1, 45893.0, 50.0, 31.0, -24236.0, -10.0, 20.0},
174 { 0, 0, 0, 2,0, 63384.0, 11.0, -150.0, -1220.0, 0.0, 29.0},
175 { 0, 0, 2, 2,2, -38571.0, -1.0, 158.0, 16452.0, -11.0, 68.0},
176
177 /* 21-30 */
178 { 0,-2, 2,-2,2, 32481.0, 0.0, 0.0, -13870.0, 0.0, 0.0},
179 {-2, 0, 0, 2,0, -47722.0, 0.0, -18.0, 477.0, 0.0, -25.0},
180 { 2, 0, 2, 0,2, -31046.0, -1.0, 131.0, 13238.0, -11.0, 59.0},
181 { 1, 0, 2,-2,2, 28593.0, 0.0, -1.0, -12338.0, 10.0, -3.0},
182 {-1, 0, 2, 0,1, 20441.0, 21.0, 10.0, -10758.0, 0.0, -3.0},
183 { 2, 0, 0, 0,0, 29243.0, 0.0, -74.0, -609.0, 0.0, 13.0},
184 { 0, 0, 2, 0,0, 25887.0, 0.0, -66.0, -550.0, 0.0, 11.0},
185 { 0, 1, 0, 0,1, -14053.0, -25.0, 79.0, 8551.0, -2.0, -45.0},
186 {-1, 0, 0, 2,1, 15164.0, 10.0, 11.0, -8001.0, 0.0, -1.0},
187 { 0, 2, 2,-2,2, -15794.0, 72.0, -16.0, 6850.0, -42.0, -5.0},
188
189 /* 31-40 */
190 { 0, 0,-2, 2,0, 21783.0, 0.0, 13.0, -167.0, 0.0, 13.0},
191 { 1, 0, 0,-2,1, -12873.0, -10.0, -37.0, 6953.0, 0.0, -14.0},
192 { 0,-1, 0, 0,1, -12654.0, 11.0, 63.0, 6415.0, 0.0, 26.0},
193 {-1, 0, 2, 2,1, -10204.0, 0.0, 25.0, 5222.0, 0.0, 15.0},
194 { 0, 2, 0, 0,0, 16707.0, -85.0, -10.0, 168.0, -1.0, 10.0},
195 { 1, 0, 2, 2,2, -7691.0, 0.0, 44.0, 3268.0, 0.0, 19.0},
196 {-2, 0, 2, 0,0, -11024.0, 0.0, -14.0, 104.0, 0.0, 2.0},
197 { 0, 1, 2, 0,2, 7566.0, -21.0, -11.0, -3250.0, 0.0, -5.0},
198 { 0, 0, 2, 2,1, -6637.0, -11.0, 25.0, 3353.0, 0.0, 14.0},
199 { 0,-1, 2, 0,2, -7141.0, 21.0, 8.0, 3070.0, 0.0, 4.0},
200
201 /* 41-50 */
202 { 0, 0, 0, 2,1, -6302.0, -11.0, 2.0, 3272.0, 0.0, 4.0},
203 { 1, 0, 2,-2,1, 5800.0, 10.0, 2.0, -3045.0, 0.0, -1.0},
204 { 2, 0, 2,-2,2, 6443.0, 0.0, -7.0, -2768.0, 0.0, -4.0},
205 {-2, 0, 0, 2,1, -5774.0, -11.0, -15.0, 3041.0, 0.0, -5.0},
206 { 2, 0, 2, 0,1, -5350.0, 0.0, 21.0, 2695.0, 0.0, 12.0},
207 { 0,-1, 2,-2,1, -4752.0, -11.0, -3.0, 2719.0, 0.0, -3.0},
208 { 0, 0, 0,-2,1, -4940.0, -11.0, -21.0, 2720.0, 0.0, -9.0},
209 {-1,-1, 0, 2,0, 7350.0, 0.0, -8.0, -51.0, 0.0, 4.0},
210 { 2, 0, 0,-2,1, 4065.0, 0.0, 6.0, -2206.0, 0.0, 1.0},
211 { 1, 0, 0, 2,0, 6579.0, 0.0, -24.0, -199.0, 0.0, 2.0},
212
213 /* 51-60 */
214 { 0, 1, 2,-2,1, 3579.0, 0.0, 5.0, -1900.0, 0.0, 1.0},
215 { 1,-1, 0, 0,0, 4725.0, 0.0, -6.0, -41.0, 0.0, 3.0},
216 {-2, 0, 2, 0,2, -3075.0, 0.0, -2.0, 1313.0, 0.0, -1.0},
217 { 3, 0, 2, 0,2, -2904.0, 0.0, 15.0, 1233.0, 0.0, 7.0},
218 { 0,-1, 0, 2,0, 4348.0, 0.0, -10.0, -81.0, 0.0, 2.0},
219 { 1,-1, 2, 0,2, -2878.0, 0.0, 8.0, 1232.0, 0.0, 4.0},
220 { 0, 0, 0, 1,0, -4230.0, 0.0, 5.0, -20.0, 0.0, -2.0},
221 {-1,-1, 2, 2,2, -2819.0, 0.0, 7.0, 1207.0, 0.0, 3.0},
222 {-1, 0, 2, 0,0, -4056.0, 0.0, 5.0, 40.0, 0.0, -2.0},
223 { 0,-1, 2, 2,2, -2647.0, 0.0, 11.0, 1129.0, 0.0, 5.0},
224
225 /* 61-70 */
226 {-2, 0, 0, 0,1, -2294.0, 0.0, -10.0, 1266.0, 0.0, -4.0},
227 { 1, 1, 2, 0,2, 2481.0, 0.0, -7.0, -1062.0, 0.0, -3.0},
228 { 2, 0, 0, 0,1, 2179.0, 0.0, -2.0, -1129.0, 0.0, -2.0},
229 {-1, 1, 0, 1,0, 3276.0, 0.0, 1.0, -9.0, 0.0, 0.0},
230 { 1, 1, 0, 0,0, -3389.0, 0.0, 5.0, 35.0, 0.0, -2.0},
231 { 1, 0, 2, 0,0, 3339.0, 0.0, -13.0, -107.0, 0.0, 1.0},
232 {-1, 0, 2,-2,1, -1987.0, 0.0, -6.0, 1073.0, 0.0, -2.0},
233 { 1, 0, 0, 0,2, -1981.0, 0.0, 0.0, 854.0, 0.0, 0.0},
234 {-1, 0, 0, 1,0, 4026.0, 0.0, -353.0, -553.0, 0.0,-139.0},
235 { 0, 0, 2, 1,2, 1660.0, 0.0, -5.0, -710.0, 0.0, -2.0},
236
237 /* 71-77 */
238 {-1, 0, 2, 4,2, -1521.0, 0.0, 9.0, 647.0, 0.0, 4.0},
239 {-1, 1, 0, 1,1, 1314.0, 0.0, 0.0, -700.0, 0.0, 0.0},
240 { 0,-2, 2,-2,1, -1283.0, 0.0, 0.0, 672.0, 0.0, 0.0},
241 { 1, 0, 2, 2,1, -1331.0, 0.0, 8.0, 663.0, 0.0, 4.0},
242 {-2, 0, 2, 2,2, 1383.0, 0.0, -2.0, -594.0, 0.0, -2.0},
243 {-1, 0, 0, 0,2, 1405.0, 0.0, 4.0, -610.0, 0.0, 2.0},
244 { 1, 1, 2,-2,2, 1290.0, 0.0, 0.0, -556.0, 0.0, 0.0}
245 };
246
247/* Number of terms in the series */
248 const int NLS = (int) (sizeof x / sizeof x[0]);
249
250/*--------------------------------------------------------------------*/
251
252/* Interval between fundamental epoch J2000.0 and given date (JC). */
253 t = ((date1 - ERFA_DJ00) + date2) / ERFA_DJC;
254
255/* --------------------*/
256/* LUNI-SOLAR NUTATION */
257/* --------------------*/
258
259/* Fundamental (Delaunay) arguments from Simon et al. (1994) */
260
261/* Mean anomaly of the Moon. */
262 el = fmod(485868.249036 + (1717915923.2178) * t, ERFA_TURNAS) * ERFA_DAS2R;
263
264/* Mean anomaly of the Sun. */
265 elp = fmod(1287104.79305 + (129596581.0481) * t, ERFA_TURNAS) * ERFA_DAS2R;
266
267/* Mean argument of the latitude of the Moon. */
268 f = fmod(335779.526232 + (1739527262.8478) * t, ERFA_TURNAS) * ERFA_DAS2R;
269
270/* Mean elongation of the Moon from the Sun. */
271 d = fmod(1072260.70369 + (1602961601.2090) * t, ERFA_TURNAS) * ERFA_DAS2R;
272
273/* Mean longitude of the ascending node of the Moon. */
274 om = fmod(450160.398036 + (-6962890.5431) * t, ERFA_TURNAS) * ERFA_DAS2R;
275
276/* Initialize the nutation values. */
277 dp = 0.0;
278 de = 0.0;
279
280/* Summation of luni-solar nutation series (smallest terms first). */
281 for (i = NLS-1; i >= 0; i--) {
282
283 /* Argument and functions. */
284 arg = fmod( (double)x[i].nl * el +
285 (double)x[i].nlp * elp +
286 (double)x[i].nf * f +
287 (double)x[i].nd * d +
288 (double)x[i].nom * om, ERFA_D2PI );
289 sarg = sin(arg);
290 carg = cos(arg);
291
292 /* Term. */
293 dp += (x[i].ps + x[i].pst * t) * sarg + x[i].pc * carg;
294 de += (x[i].ec + x[i].ect * t) * carg + x[i].es * sarg;
295 }
296
297/* Convert from 0.1 microarcsec units to radians. */
298 dpsils = dp * U2R;
299 depsls = de * U2R;
300
301/* ------------------------------*/
302/* IN LIEU OF PLANETARY NUTATION */
303/* ------------------------------*/
304
305/* Fixed offset to correct for missing terms in truncated series. */
306 dpsipl = DPPLAN;
307 depspl = DEPLAN;
308
309/* --------*/
310/* RESULTS */
311/* --------*/
312
313/* Add luni-solar and planetary components. */
314 *dpsi = dpsils + dpsipl;
315 *deps = depsls + depspl;
316
317 return;
318
319}
320/*----------------------------------------------------------------------
321**
322**
323** Copyright (C) 2013-2017, NumFOCUS Foundation.
324** All rights reserved.
325**
326** This library is derived, with permission, from the International
327** Astronomical Union's "Standards of Fundamental Astronomy" library,
328** available from http://www.iausofa.org.
329**
330** The ERFA version is intended to retain identical functionality to
331** the SOFA library, but made distinct through different function and
332** file names, as set out in the SOFA license conditions. The SOFA
333** original has a role as a reference standard for the IAU and IERS,
334** and consequently redistribution is permitted only in its unaltered
335** state. The ERFA version is not subject to this restriction and
336** therefore can be included in distributions which do not support the
337** concept of "read only" software.
338**
339** Although the intent is to replicate the SOFA API (other than
340** replacement of prefix names) and results (with the exception of
341** bugs; any that are discovered will be fixed), SOFA is not
342** responsible for any errors found in this version of the library.
343**
344** If you wish to acknowledge the SOFA heritage, please acknowledge
345** that you are using a library derived from SOFA, rather than SOFA
346** itself.
347**
348**
349** TERMS AND CONDITIONS
350**
351** Redistribution and use in source and binary forms, with or without
352** modification, are permitted provided that the following conditions
353** are met:
354**
355** 1 Redistributions of source code must retain the above copyright
356** notice, this list of conditions and the following disclaimer.
357**
358** 2 Redistributions in binary form must reproduce the above copyright
359** notice, this list of conditions and the following disclaimer in
360** the documentation and/or other materials provided with the
361** distribution.
362**
363** 3 Neither the name of the Standards Of Fundamental Astronomy Board,
364** the International Astronomical Union nor the names of its
365** contributors may be used to endorse or promote products derived
366** from this software without specific prior written permission.
367**
368** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
369** "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
370** LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
371** FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
372** COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
373** INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
374** BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
375** LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
376** CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
377** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
378** ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
379** POSSIBILITY OF SUCH DAMAGE.
380**
381*/
Note: See TracBrowser for help on using the repository browser.