1   /* Copyright 2002-2020 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.utils;
18  
19  
20  import org.hipparchus.analysis.UnivariateFunction;
21  import org.hipparchus.analysis.UnivariateVectorFunction;
22  import org.hipparchus.analysis.differentiation.DSFactory;
23  import org.hipparchus.analysis.differentiation.DerivativeStructure;
24  import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
25  import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
26  import org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction;
27  import org.hipparchus.util.Decimal64;
28  import org.hipparchus.util.FastMath;
29  import org.hipparchus.util.MathUtils;
30  import org.junit.Assert;
31  import org.junit.Before;
32  import org.junit.Test;
33  import org.orekit.Utils;
34  import org.orekit.frames.EOPHistory;
35  import org.orekit.frames.FramesFactory;
36  import org.orekit.frames.ITRFVersion;
37  import org.orekit.time.AbsoluteDate;
38  import org.orekit.time.TimeVectorFunction;
39  import org.orekit.time.DateComponents;
40  import org.orekit.time.FieldAbsoluteDate;
41  import org.orekit.time.TimeScalarFunction;
42  import org.orekit.time.TimeScale;
43  import org.orekit.time.TimeScalesFactory;
44  import org.orekit.time.UT1Scale;
45  
46  
47  public class IERSConventionsTest {
48  
49      @Test
50      public void testConventionsNumber() {
51          Assert.assertEquals(3, IERSConventions.values().length);
52      }
53  
54      @Test
55      public void testIERS1996NutationAngles() {
56  
57          // the reference value has been computed using the March 2012 version of the SOFA library
58          // http://www.iausofa.org/2012_0301_C.html, with the following code
59          //
60          //        double utc1, utc2, tai1, tai2, tt1, tt2, dpsi, deps, epsa;
61          //
62          //        // 2004-02-14:00:00:00Z, MJD = 53049, UT1-UTC = -0.4093509
63          //        utc1  = DJM0 + 53049.0;
64          //        utc2  = 0.0;
65          //        iauUtctai(utc1, utc2, &tai1, &tai2);
66          //        iauTaitt(tai1, tai2, &tt1, &tt2);
67          //
68          //        iauNut80(tt1, tt2, &dpsi, &deps);
69          //        epsa = iauObl80(tt1, tt2);
70          //
71          //        printf("iauNut80(%.20g, %.20g, dpsi, deps)\n  --> %.20g %.20g\n",
72          //               tt1, tt2, dpsi, deps);
73          //        printf("iau0bl80(%.20g, %.20g)\n  --> %.20g\n",
74          //               tt1, tt2, epsa);
75          //
76          // the output of this test reads:
77          //      iauNut80(2453049.5, 0.00074287037037037029902, dpsi, deps)
78          //        --> -5.3059154211478291722e-05 3.2051803135750973851e-05
79          //      iau0bl80(2453049.5, 0.00074287037037037029902)
80          //        --> 0.40908345528446415917
81  
82          // the thresholds below have been set up large enough to have the test pass with
83          // the default Orekit setting and the reference SOFA setting. There are implementation
84          // differences between the two libraries, as the Delaunay parameters are not the
85          // same (in Orekit, we are compliant with page 23 in chapter 5 of IERS conventions 1996.
86          // If the content of the IERS-conventions/1996/nutation-arguments.txt file by Orekit is
87          // changed to match SOFA setting instead of IERS conventions as follows:
88          //
89          //        # Mean Anomaly of the Moon
90          //        F1 ≡ l = 485866.733″ + 1717915922.633″t + 31.310″t² + 0.064″t³
91          //
92          //       # Mean Anomaly of the Sun
93          //        F2 ≡ l' = 1287099.804″ + 129596581.224″t − 0.577″t² − 0.012t³
94          //
95          //       # L − Ω
96          //        F3 ≡ F = 335778.877″ + 1739527263.137″t − 13.257″t² + 0.011″t³
97          //
98          //       # Mean Elongation of the Moon from the Sun
99          //        F4 ≡ D = 1072261.307″ + 1602961601.328″t − 6.891″t² + 0.019″t³
100         //
101         //       # Mean Longitude of the Ascending Node of the Moon
102         //        F5 ≡ Ω = 450160.280″ − 6962890.539″t + 7.455″t² + 0.008″t³
103         //
104         // then the thresholds for the test can be reduced to 3e-13 for ∆ψ and 7e-14 for ∆ε.
105         // We decided to stick with IERS published reference values for the default Orekit setting.
106         // The differences are nevertheless quite small (4.8e-11 radians is sub-millimeter level
107         // in low Earth orbit).
108         AbsoluteDate date = new AbsoluteDate(2004, 2, 14, TimeScalesFactory.getUTC());
109         double[] angles= IERSConventions.IERS_1996.getNutationFunction().value(date);
110         Assert.assertEquals(-5.3059154211478291722e-05, angles[0], 4.8e-11); // 3e-13 with SOFA values
111         Assert.assertEquals(3.2051803135750973851e-05,  angles[1], 1.3e-11); // 7e-14 with SOFA values
112 
113         double epsilonA = IERSConventions.IERS_1996.getMeanObliquityFunction().value(date);
114         Assert.assertEquals(0.40908345528446415917,     epsilonA, 1.0e-16);
115 
116     }
117 
118     @Test
119     public void testGMST82FieldConsistency() {
120         final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_1996, true);
121         checkScalarFunctionConsistency(IERSConventions.IERS_1996.getGMSTFunction(ut1),
122                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
123                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 9.0e-16, 8.0e-13);
124     }
125 
126     @Test
127     public void testGMST82Sofa() {
128 
129         // the reference value has been computed using the March 2012 version of the SOFA library
130         // http://www.iausofa.org/2012_0301_C.html, with the following code
131         //
132         //        double utc1, utc2, ut11, ut12, gmst;
133         //
134         //        // 2004-02-14:00:00:00Z, MJD = 53049, UT1-UTC = -0.4093509
135         //        utc1  = DJM0 + 53049.0;
136         //        utc2  = 0.0;
137         //        iauUtcut1(utc1, utc2, -0.4093509, &ut11, &ut12);
138         //        gmst = iauGmst82(ut11, ut12);
139         //        printf("iaugmst82(%.20g, %.20g)\n  --> %.20g\n", ut11, ut12, gmst);
140         //
141         //        // 2004-02-29:00:00:00Z, MJD = 53064, UT1-UTC = -0.4175723
142         //        utc1 = DJM0 + 53064.0;
143         //        utc2 = 0.0;
144         //        iauUtcut1(utc1, utc2, -0.4175723, &ut11, &ut12);
145         //        gmst = iauGmst82(ut11, ut12);
146         //        printf("iaugmst82(%.20g, %.20g)\n  --> %.20g\n", ut11, ut12, gmst);
147         //
148         // the output of this test reads:
149         //      iaugmst82(2453049.5, -4.7378576388888813016e-06)
150         //        --> 2.5021977627453466653
151         //      iaugmst82(2453064.5, -4.8330127314815448519e-06)
152         //        --> 2.7602390405411441066
153 
154         final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_1996, true);
155         final TimeScalarFunction gmst82 = IERSConventions.IERS_1996.getGMSTFunction(ut1);
156         AbsoluteDate date = new AbsoluteDate(2004, 2, 14, TimeScalesFactory.getUTC());
157         double gmst = MathUtils.normalizeAngle(gmst82.value(date), 0.0);
158         Assert.assertEquals(2.5021977627453466653, gmst, 2.0e-13);
159         date = new AbsoluteDate(2004, 2, 29, TimeScalesFactory.getUTC());
160         gmst = MathUtils.normalizeAngle(gmst82.value(date), 0.0);
161         Assert.assertEquals(2.7602390405411441066, gmst, 4.0e-13);
162 
163     }
164 
165     @Test
166     public void testGMST00FieldConsistency() {
167         final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2003, true);
168         checkScalarFunctionConsistency(IERSConventions.IERS_2003.getGMSTFunction(ut1),
169                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
170                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 2.0e-15, 8.0e-13);
171     }
172 
173     @Test
174     public void testGMST00Sofa() {
175 
176         // the reference value has been computed using the March 2012 version of the SOFA library
177         // http://www.iausofa.org/2012_0301_C.html, with the following code
178         //
179         //        double utc1, utc2, tai1, tai2, tt1, tt2, ut11, ut12, gmst;
180         //
181         //        // 2004-02-14:00:00:00Z, MJD = 53049, UT1-UTC = -0.4093509
182         //        utc1  = DJM0 + 53049.0;
183         //        utc2  = 0.0;
184         //        iauUtctai(utc1, utc2, &tai1, &tai2);
185         //        iauTaitt(tai1, tai2, &tt1, &tt2);
186         //        iauUtcut1(utc1, utc2, -0.4093509, &ut11, &ut12);
187         //        gmst = iauGmst00(ut11, ut12, tt1, tt2);
188         //        printf("iaugmst00(%.20g, %.20g, %.20g, %.20g)\n  --> %.20g\n",
189         //               ut11, ut12, tt1, tt2, gmst);
190         //
191         //        // 2004-02-29:00:00:00Z, MJD = 53064, UT1-UTC = -0.4175723
192         //        utc1 = DJM0 + 53064.0;
193         //        utc2 = 0.0;
194         //        iauUtctai(utc1, utc2, &tai1, &tai2);
195         //        iauTaitt(tai1, tai2, &tt1, &tt2);
196         //        iauUtcut1(utc1, utc2, -0.4175723, &ut11, &ut12);
197         //        gmst = iauGmst00(ut11, ut12, tt1, tt2);
198         //        printf("iaugmst00(%.20g, %.20g, %.20g, %.20g)\n  --> %.20g\n",
199         //               ut11, ut12, tt1, tt2, gmst);
200         //
201         // the output of this test reads:
202         //      iaugmst00(2453049.5, -4.7378576388888813016e-06, 2453049.5, 0.00074287037037037029902)
203         //        --> 2.5021977786243789765
204         //      iaugmst00(2453064.5, -4.8330127314815448519e-06, 2453064.5, 0.00074287037037037029902)
205         //        --> 2.7602390558728311376
206 
207         final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2003, true);
208         final TimeScalarFunction gmst00 = IERSConventions.IERS_2003.getGMSTFunction(ut1);
209         AbsoluteDate date = new AbsoluteDate(2004, 2, 14, TimeScalesFactory.getUTC());
210         double gmst = MathUtils.normalizeAngle(gmst00.value(date), 0.0);
211         Assert.assertEquals(2.5021977786243789765, gmst, 1.0e-15);
212         date = new AbsoluteDate(2004, 2, 29, TimeScalesFactory.getUTC());
213         gmst = MathUtils.normalizeAngle(gmst00.value(date), 0.0);
214         Assert.assertEquals(2.7602390558728311376, gmst, 1.0e-15);
215 
216     }
217 
218     @Test
219     public void testGMST06FieldConsistency() {
220         final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
221         checkScalarFunctionConsistency(IERSConventions.IERS_2010.getGMSTFunction(ut1),
222                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
223                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 2.0e-15, 8.0e-13);
224     }
225 
226     @Test
227     public void testGMST06Sofa() {
228 
229         // the reference value has been computed using the March 2012 version of the SOFA library
230         // http://www.iausofa.org/2012_0301_C.html, with the following code
231         //
232         //        double utc1, utc2, tai1, tai2, tt1, tt2, ut11, ut12, gmst;
233         //
234         //        // 2004-02-14:00:00:00Z, MJD = 53049, UT1-UTC = -0.4093509
235         //        utc1  = DJM0 + 53049.0;
236         //        utc2  = 0.0;
237         //        iauUtctai(utc1, utc2, &tai1, &tai2);
238         //        iauTaitt(tai1, tai2, &tt1, &tt2);
239         //        iauUtcut1(utc1, utc2, -0.4093509, &ut11, &ut12);
240         //        gmst = iauGmst06(ut11, ut12, tt1, tt2);
241         //        printf("iaugmst06(%.20g, %.20g, %.20g, %.20g)\n  --> %.20g\n",
242         //               ut11, ut12, tt1, tt2, gmst);
243         //
244         //        // 2004-02-29:00:00:00Z, MJD = 53064, UT1-UTC = -0.4175723
245         //        utc1 = DJM0 + 53064.0;
246         //        utc2 = 0.0;
247         //        iauUtctai(utc1, utc2, &tai1, &tai2);
248         //        iauTaitt(tai1, tai2, &tt1, &tt2);
249         //        iauUtcut1(utc1, utc2, -0.4175723, &ut11, &ut12);
250         //        gmst = iauGmst06(ut11, ut12, tt1, tt2);
251         //        printf("iaugmst06(%.20g, %.20g, %.20g, %.20g)\n  --> %.20g\n",
252         //               ut11, ut12, tt1, tt2, gmst);
253         //
254         // the output of this test reads:
255         //      iaugmst06(2453049.5, -4.7378576388888813016e-06, 2453049.5, 0.00074287037037037029902)
256         //        --> 2.5021977784096232078
257         //      iaugmst06(2453064.5, -4.8330127314815448519e-06, 2453064.5, 0.00074287037037037029902)
258         //        --> 2.7602390556555129741
259 
260         final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
261         final TimeScalarFunction gmst06 = IERSConventions.IERS_2010.getGMSTFunction(ut1);
262         AbsoluteDate date = new AbsoluteDate(2004, 2, 14, TimeScalesFactory.getUTC());
263         double gmst = MathUtils.normalizeAngle(gmst06.value(date), 0.0);
264         Assert.assertEquals(2.5021977784096232078, gmst, 1.0e-15);
265         date = new AbsoluteDate(2004, 2, 29, TimeScalesFactory.getUTC());
266         gmst = MathUtils.normalizeAngle(gmst06.value(date), 0.0);
267         Assert.assertEquals(2.7602390556555129741, gmst, 3.0e-15);
268 
269     }
270 
271     @Test
272     public void testERA1996FieldConsistency() {
273         final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_1996, true);
274         checkScalarFunctionConsistency(IERSConventions.IERS_1996.getEarthOrientationAngleFunction(ut1),
275                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
276                                        0.8 * Constants.JULIAN_DAY, 1800.0, 10.0, 2.0e-15, 8.0e-13);
277     }
278 
279     @Test
280     public void testERA2003FieldConsistency() {
281         final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2003, true);
282         checkScalarFunctionConsistency(IERSConventions.IERS_2003.getEarthOrientationAngleFunction(ut1),
283                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
284                                        0.8 * Constants.JULIAN_DAY, 1800.0, 10.0, 2.0e-15, 8.0e-13);
285     }
286 
287     @Test
288     public void testERA2010FieldConsistency() {
289         final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
290         checkScalarFunctionConsistency(IERSConventions.IERS_2010.getEarthOrientationAngleFunction(ut1),
291                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
292                                        0.8 * Constants.JULIAN_DAY, 1800.0, 10.0, 2.0e-15, 8.0e-13);
293     }
294 
295     @Test
296     public void testERASofa() {
297 
298         // the reference value has been computed using the March 2012 version of the SOFA library
299         // http://www.iausofa.org/2012_0301_C.html, with the following code
300         //
301         //        double tai1, tai2, ut11, ut12, era, taiutc;
302         //        taiutc = 32.0;
303         //
304         //        // 2004-02-14:00:00:00Z, MJD = 53049, UT1-UTC = -0.4093509
305         //        tai1 = DJM0 + 53049.0;
306         //        tai2 = taiutc / DAYSEC;
307         //        iauTaiut1(tai1, tai2, -0.4093509 - taiutc, &ut11, &ut12);
308         //        era = iauEra00(ut11, ut12);
309         //        printf("iauera00(%.20g, %.20g)\n  --> %.20g\n", ut11, ut12, era);
310         //
311         //        // 2004-02-29:00:00:00Z, MJD = 53064, UT1-UTC = -0.4175723
312         //        tai1 = DJM0 + 53064.0;
313         //        tai2 = taiutc / DAYSEC;
314         //        iauTaiut1(tai1, tai2, -0.4175723 - taiutc, &ut11, &ut12);
315         //        era = iauEra00(ut11, ut12);
316         //        printf("iauera00(%.20g, %.20g)\n  --> %.20g\n", ut11, ut12, era);
317         //
318         // the output of this test reads:
319         //      iauera00(2453049.5, -4.7378576388888813016e-06)
320         //        --> 2.5012766511308228701
321         //      iauera00(2453064.5, -4.8330127314815448519e-06)
322         //        --> 2.7593087452455264952
323 
324         final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2003, true);
325         final TimeScalarFunction era00 =
326                 IERSConventions.IERS_2003.getEarthOrientationAngleFunction(ut1);
327         AbsoluteDate date = new AbsoluteDate(2004, 2, 14, TimeScalesFactory.getUTC());
328         double era  = MathUtils.normalizeAngle(era00.value(date), 0.0);
329         Assert.assertEquals(2.5012766511308228701, era, 1.0e-15);
330         date = new AbsoluteDate(2004, 2, 29, TimeScalesFactory.getUTC());
331         era  = MathUtils.normalizeAngle(era00.value(date), 0.0);
332         Assert.assertEquals(2.7593087452455264952, era, 1.0e-15);
333 
334     }
335 
336     @Test
337     public void testGST94FieldConsistency() {
338         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_1996, true);
339         checkScalarFunctionConsistency(IERSConventions.IERS_1996.getGASTFunction(TimeScalesFactory.getUT1(eopHistory),
340                                                                                  eopHistory),
341                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
342                                        0.8 * Constants.JULIAN_DAY, 1800.0, 10.0, 9.0e-16, 8.0e-13);
343     }
344 
345     @Test
346     public void testGST94FieldConsistencyNoEOP() {
347         checkScalarFunctionConsistency(IERSConventions.IERS_1996.getGASTFunction(TimeScalesFactory.getUT1(null), null),
348                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
349                                        0.8 * Constants.JULIAN_DAY, 1800.0, 10.0, 9.0e-16, 2.0e-16);
350     }
351 
352     @Test
353     public void testGST94Sofa() {
354 
355         // the reference value has been computed using the March 2012 version of the SOFA library
356         // http://www.iausofa.org/2012_0301_C.html, with the following code
357         //
358         //        double utc1, utc2, ut11, ut12, tai1, tai2, tt1, tt2, gmst82, eqeq94;
359         //
360         //        // 2004-02-14:00:00:00Z, MJD = 53049, UT1-UTC = -0.4093509
361         //        utc1  = DJM0 + 53049.0;
362         //        utc2  = 0.0;
363         //        iauUtcut1(utc1, utc2, -0.4093509, &ut11, &ut12);
364         //        iauUtctai(utc1, utc2, &tai1, &tai2);
365         //        iauTaitt(tai1, tai2, &tt1, &tt2);
366         //        gmst82 = iauGmst82(ut11, ut12);
367         //        eqeq94 = iauEqeq94(tt1, tt2);
368         //        printf("iauGmst82(%.20g, %.20g)\n  --> %.20g\n", ut11, ut12, gmst82);
369         //        printf("iauEqeq94(%.20g, %.20g)\n  --> %.20g\n", tt1, tt2, eqeq94);
370         //        printf(" gmst82 + eqeq94  --> %.20g\n", gmst82 + eqeq94);
371         //
372         //        // 2004-02-29:00:00:00Z, MJD = 53064, UT1-UTC = -0.4175723
373         //        utc1 = DJM0 + 53064.0;
374         //        utc2 = 0.0;
375         //        iauUtcut1(utc1, utc2, -0.4175723, &ut11, &ut12);
376         //        iauUtctai(utc1, utc2, &tai1, &tai2);
377         //        iauTaitt(tai1, tai2, &tt1, &tt2);
378         //        gmst82 = iauGmst82(ut11, ut12);
379         //        eqeq94 = iauEqeq94(tt1, tt2);
380         //        printf("iauGmst82(%.20g, %.20g)\n  --> %.20g\n", ut11, ut12, gmst82);
381         //        printf("iauEqeq94(%.20g, %.20g)\n  --> %.20g\n", tt1, tt2, eqeq94);
382         //        printf(" gmst82 + eqeq94  --> %.20g\n", gmst82 + eqeq94);
383         //
384         // the output of this test reads:
385         //      iauGmst82(2453049.5, -4.7378576388888813016e-06)
386         //        --> 2.5021977627453466653
387         //      iauEqeq94(2453049.5, 0.00074287037037037029902)
388         //        --> -4.8671604682267536886e-05
389         //       gmst82 + eqeq94  --> 2.5021490911406645274
390         //      iauGmst82(2453064.5, -4.8330127314815448519e-06)
391         //        --> 2.7602390405411441066
392         //      iauEqeq94(2453064.5, 0.00074287037037037029902)
393         //        --> -4.8893054690771762302e-05
394         //       gmst82 + eqeq94  --> 2.7601901474864534158
395 
396         // As can be seen in the code above, we didn't call iauGst94, because as
397         // stated in the function header comment in SOFA source files:
398         //        "accuracy has been compromised for the sake of
399         //         convenience in that UT is used instead of TDB (or TT) to compute
400         //         the equation of the equinoxes."
401         // So we rather performed the date conversion and then called ourselves iauGmst82
402         // with a date in UTC and eqe94 with a date in TT, restoring full accuracy in the
403         // SOFA computation for this test.
404 
405         // The thresholds below have been set up large enough to have the test pass with
406         // the default Orekit setting and the reference SOFA setting. There are implementation
407         // differences between the two libraries, as the Delaunay parameters are not the
408         // same (in Orekit, we are compliant with page 23 in chapter 5 of IERS conventions 1996.
409         // If the content of the IERS-conventions/1996/nutation-arguments.txt file by Orekit is
410         // changed to match SOFA setting instead of IERS conventions as follows:
411         //
412         //        # Mean Anomaly of the Moon
413         //        F1 ≡ l = 485866.733″ + 1717915922.633″t + 31.310″t² + 0.064″t³
414         //
415         //       # Mean Anomaly of the Sun
416         //        F2 ≡ l' = 1287099.804″ + 129596581.224″t − 0.577″t² − 0.012t³
417         //
418         //       # L − Ω
419         //        F3 ≡ F = 335778.877″ + 1739527263.137″t − 13.257″t² + 0.011″t³
420         //
421         //       # Mean Elongation of the Moon from the Sun
422         //        F4 ≡ D = 1072261.307″ + 1602961601.328″t − 6.891″t² + 0.019″t³
423         //
424         //       # Mean Longitude of the Ascending Node of the Moon
425         //        F5 ≡ Ω = 450160.280″ − 6962890.539″t + 7.455″t² + 0.008″t³
426         //
427         // then the thresholds for the test can be reduced to 4.2e-13 for the first test,
428         // and 4.7e-13 for the second test.
429         // We decided to stick with IERS published reference values for the default Orekit setting.
430         // The differences are nevertheless quite small (9e-11 radians is sub-millimeter level
431         // in low Earth orbit).
432         Utils.setLoaders(IERSConventions.IERS_1996,
433                          Utils.buildEOPList(IERSConventions.IERS_1996, ITRFVersion.ITRF_2008, new double[][] {
434                              { 53047, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
435                              { 53048, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
436                              { 53049, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
437                              { 53050, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
438                              { 53051, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
439                              { 53052, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
440                              { 53053, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
441                              { 53054, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
442                              // ...
443                              { 53059, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
444                              { 53060, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
445                              { 53061, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
446                              { 53062, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
447                              { 53063, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
448                              { 53064, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
449                              { 53065, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
450                              { 53066, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 }
451                          }));
452         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_1996, true);
453         final TimeScalarFunction gst94 =
454                 IERSConventions.IERS_1996.getGASTFunction(TimeScalesFactory.getUT1(eopHistory),
455                                                           eopHistory);
456         AbsoluteDate date = new AbsoluteDate(2004, 2, 14, TimeScalesFactory.getUTC());
457         double gst = MathUtils.normalizeAngle(gst94.value(date), 0.0);
458         Assert.assertEquals(2.5021490911406645274, gst, 4.4e-11); // 4.2e-13 with SOFA values
459         date = new AbsoluteDate(2004, 2, 29, TimeScalesFactory.getUTC());
460         gst = MathUtils.normalizeAngle(gst94.value(date), 0.0);
461         Assert.assertEquals(2.7601901474864534158, gst, 9.0e-11); // 4.7e-13 with SOFA values
462 
463     }
464 
465     @Test
466     public void testGST00AFieldConsistency() {
467         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_2003, true);
468         checkScalarFunctionConsistency(IERSConventions.IERS_2003.getGASTFunction(TimeScalesFactory.getUT1(eopHistory), eopHistory),
469                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
470                                        0.8 * Constants.JULIAN_DAY, 1800.0, 10.0, 2.0e-15, 8.0e-13);
471     }
472 
473     @Test
474     public void testGST00AFieldConsistencyNoEOP() {
475         checkScalarFunctionConsistency(IERSConventions.IERS_2003.getGASTFunction(TimeScalesFactory.getUT1(null), null),
476                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
477                                        0.8 * Constants.JULIAN_DAY, 1800.0, 10.0, 2.0e-15, 2.0e-16);
478     }
479 
480     @Test
481     public void testGST00ASofa() {
482 
483         // the reference value has been computed using the March 2012 version of the SOFA library
484         // http://www.iausofa.org/2012_0301_C.html, with the following code
485         //
486         //        double utc1, utc2, tai1, tai2, tt1, tt2, ut11, ut12, gst;
487         //
488         //        // 2004-02-14:00:00:00Z, MJD = 53049, UT1-UTC = -0.4093509
489         //        utc1  = DJM0 + 53049.0;
490         //        utc2  = 0.0;
491         //        iauUtctai(utc1, utc2, &tai1, &tai2);
492         //        iauTaitt(tai1, tai2, &tt1, &tt2);
493         //        iauUtcut1(utc1, utc2, -0.4093509, &ut11, &ut12);
494         //        gst = iauGst00a(ut11, ut12, tt1, tt2);
495         //        printf("iaugst00a(%.20g, %.20g, %.20g, %.20g)\n  --> %.20g\n",
496         //               ut11, ut12, tt1, tt2, gst);
497         //
498         //        // 2004-02-29:00:00:00Z, MJD = 53064, UT1-UTC = -0.4175723
499         //        utc1 = DJM0 + 53064.0;
500         //        utc2 = 0.0;
501         //        iauUtctai(utc1, utc2, &tai1, &tai2);
502         //        iauTaitt(tai1, tai2, &tt1, &tt2);
503         //        iauUtcut1(utc1, utc2, -0.4175723, &ut11, &ut12);
504         //        gst = iauGst00a(ut11, ut12, tt1, tt2);
505         //        printf("iaugst00a(%.20g, %.20g, %.20g, %.20g)\n  --> %.20g\n",
506         //               ut11, ut12, tt1, tt2, gst);
507         //
508         // the output of this test reads:
509         //      iaugst00a(2453049.5, -4.7378576388888813016e-06, 2453049.5, 0.00074287037037037029902)
510         //        --> 2.5021491024360624778
511         //      iaugst00a(2453064.5, -4.8330127314815448519e-06, 2453064.5, 0.00074287037037037029902)
512         //        --> 2.7601901615614221619
513 
514         Utils.setLoaders(IERSConventions.IERS_2003,
515                          Utils.buildEOPList(IERSConventions.IERS_2003, ITRFVersion.ITRF_2008, new double[][] {
516                              { 53047, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
517                              { 53048, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
518                              { 53049, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
519                              { 53050, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
520                              { 53051, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
521                              { 53052, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
522                              { 53053, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
523                              { 53054, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
524                              // ...
525                              { 53059, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
526                              { 53060, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
527                              { 53061, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
528                              { 53062, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
529                              { 53063, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
530                              { 53064, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
531                              { 53065, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
532                              { 53066, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 }
533                          }));
534         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_2003, true);
535         final TimeScalarFunction gst00a =
536                 IERSConventions.IERS_2003.getGASTFunction(TimeScalesFactory.getUT1(eopHistory), eopHistory);
537         AbsoluteDate date = new AbsoluteDate(2004, 2, 14, TimeScalesFactory.getUTC());
538         double gst = MathUtils.normalizeAngle(gst00a.value(date), 0.0);
539         Assert.assertEquals(2.5021491024360624778, gst, 2.0e-13);
540         date = new AbsoluteDate(2004, 2, 29, TimeScalesFactory.getUTC());
541         gst = MathUtils.normalizeAngle(gst00a.value(date), 0.0);
542         Assert.assertEquals(2.7601901615614221619, gst, 3.0e-13);
543 
544     }
545 
546     @Test
547     public void testGST06FieldConsistency() {
548         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_2010, true);
549         checkScalarFunctionConsistency(IERSConventions.IERS_2010.getGASTFunction(TimeScalesFactory.getUT1(eopHistory), eopHistory),
550                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
551                                        0.8 * Constants.JULIAN_DAY, 1800.0, 10.0, 2.0e-15, 8.0e-13);
552     }
553 
554     @Test
555     public void testGST06FieldConsistencyNoEOP() {
556         checkScalarFunctionConsistency(IERSConventions.IERS_2010.getGASTFunction(TimeScalesFactory.getUT1(null), null),
557                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
558                                        0.8 * Constants.JULIAN_DAY, 1800.0, 10.0, 2.0e-15, 2.0e-16);
559     }
560 
561     @Test
562     public void testGST06Sofa() {
563 
564         // the reference value has been computed using the March 2012 version of the SOFA library
565         // http://www.iausofa.org/2012_0301_C.html, with the following code
566         //
567         //        double utc1, utc2, tai1, tai2, tt1, tt2, ut11, ut12, gst;
568         //
569         //        // 2004-02-14:00:00:00Z, MJD = 53049, UT1-UTC = -0.4093509
570         //        utc1  = DJM0 + 53049.0;
571         //        utc2  = 0.0;
572         //        iauUtctai(utc1, utc2, &tai1, &tai2);
573         //        iauTaitt(tai1, tai2, &tt1, &tt2);
574         //        iauUtcut1(utc1, utc2, -0.4093509, &ut11, &ut12);
575         //        gst = iauGst06a(ut11, ut12, tt1, tt2);
576         //        printf("iaugst06a(%.20g, %.20g, %.20g, %.20g)\n  --> %.20g\n",
577         //               ut11, ut12, tt1, tt2, gst);
578         //
579         //        // 2004-02-29:00:00:00Z, MJD = 53064, UT1-UTC = -0.4175723
580         //        utc1 = DJM0 + 53064.0;
581         //        utc2 = 0.0;
582         //        iauUtctai(utc1, utc2, &tai1, &tai2);
583         //        iauTaitt(tai1, tai2, &tt1, &tt2);
584         //        iauUtcut1(utc1, utc2, -0.4175723, &ut11, &ut12);
585         //        gst = iauGst06a(ut11, ut12, tt1, tt2);
586         //        printf("iaugst06a(%.20g, %.20g, %.20g, %.20g)\n  --> %.20g\n",
587         //               ut11, ut12, tt1, tt2, gst);
588         //
589         // the output of this test reads:
590         //      iaugst06a(2453049.5, -4.7378576388888813016e-06, 2453049.5, 0.00074287037037037029902)
591         //        --> 2.5021491022006503435
592         //      iaugst06a(2453064.5, -4.8330127314815448519e-06, 2453064.5, 0.00074287037037037029902)
593         //        --> 2.7601901613234058885
594 
595         Utils.setLoaders(IERSConventions.IERS_2010,
596                          Utils.buildEOPList(IERSConventions.IERS_2010, ITRFVersion.ITRF_2008, new double[][] {
597                              { 53047, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
598                              { 53048, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
599                              { 53049, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
600                              { 53050, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
601                              { 53051, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
602                              { 53052, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
603                              { 53053, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
604                              { 53054, -0.4093509, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
605                              // ...
606                              { 53059, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
607                              { 53060, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
608                              { 53061, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
609                              { 53062, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
610                              { 53063, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
611                              { 53064, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
612                              { 53065, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 },
613                              { 53066, -0.4175723, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 }
614                          }));
615         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_2010, true);
616         final TimeScalarFunction gst06 =
617                 IERSConventions.IERS_2010.getGASTFunction(TimeScalesFactory.getUT1(eopHistory), eopHistory);
618         AbsoluteDate date = new AbsoluteDate(2004, 2, 14, TimeScalesFactory.getUTC());
619         double gst = MathUtils.normalizeAngle(gst06.value(date), 0.0);
620         Assert.assertEquals(2.5021491022006503435, gst, 1.3e-12);
621         date = new AbsoluteDate(2004, 2, 29, TimeScalesFactory.getUTC());
622         gst = MathUtils.normalizeAngle(gst06.value(date), 0.0);
623         Assert.assertEquals(2.7601901613234058885, gst, 1.2e-12);
624 
625     }
626 
627     @Test
628     public void testGMST2000vs82() {
629 
630         final TimeScalarFunction gmst82 =
631                 IERSConventions.IERS_1996.getGMSTFunction(TimeScalesFactory.getUT1(IERSConventions.IERS_1996, true));
632         final TimeScalarFunction gmst00 =
633                 IERSConventions.IERS_2003.getGMSTFunction(TimeScalesFactory.getUT1(IERSConventions.IERS_2003, true));
634         final DSFactory factory = new DSFactory(1, 1);
635         final FieldAbsoluteDate<DerivativeStructure> ds2000 =
636                         FieldAbsoluteDate.getJ2000Epoch(factory.getDerivativeField());
637         for (double dt = 0; dt < 10 * Constants.JULIAN_YEAR; dt += 10 * Constants.JULIAN_DAY) {
638             FieldAbsoluteDate<DerivativeStructure> date = ds2000.shiftedBy(factory.variable(0, dt));
639             DerivativeStructure delta = gmst00.value(date).subtract(gmst82.value(date));
640             // GMST82 and GMST2000 are similar to about 15 milli-arcseconds between 2000 and 2010
641             Assert.assertEquals(0.0, MathUtils.normalizeAngle(delta.getValue(), 0.0), 7.1e-8);
642             Assert.assertEquals(0.0, delta.getPartialDerivative(1), 1.0e-15);
643         }
644     }
645 
646     @Test
647     public void testGMST2000vs2006() {
648 
649         final UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
650         final TimeScalarFunction gmst00 = IERSConventions.IERS_2003.getGMSTFunction(ut1);
651         final TimeScalarFunction gmst06 = IERSConventions.IERS_2010.getGMSTFunction(ut1);
652         final DSFactory factory = new DSFactory(1, 1);
653         final FieldAbsoluteDate<DerivativeStructure> ds2000 =
654                         FieldAbsoluteDate.getJ2000Epoch(factory.getDerivativeField());
655         for (double dt = 0; dt < 10 * Constants.JULIAN_YEAR; dt += 10 * Constants.JULIAN_DAY) {
656             FieldAbsoluteDate<DerivativeStructure> date = ds2000.shiftedBy(factory.variable(0, dt));
657             DerivativeStructure delta = gmst06.value(date).subtract(gmst00.value(date));
658             // GMST2006 and GMST2000 are similar to about 0.15 milli-arcseconds between 2000 and 2010
659             Assert.assertEquals(0.0, MathUtils.normalizeAngle(delta.getValue(), 0.0), 7e-10);
660             Assert.assertEquals(0.0, delta.getPartialDerivative(1), 3.0e-18);
661         }
662     }
663 
664     @Test
665     public void testIAU1994ResolutionC7Discontinuity() {
666         TimeVectorFunction nutation = IERSConventions.IERS_1996.getNutationFunction();
667         AbsoluteDate switchDate = new AbsoluteDate(1997, 2, 27, TimeScalesFactory.getUTC());
668         double h = 0.01;
669         for (double dt = -1.0 - h / 2; dt <= 1.0 + h /2; dt += h) {
670             AbsoluteDate d = switchDate.shiftedBy(dt);
671             final double currentCorr = nutation.value(d)[2];
672                 if (dt < 0) {
673                     Assert.assertEquals(0.0, currentCorr, 1.0e-20);
674                 } else {
675                     Assert.assertEquals(-7.87098e-12, currentCorr, 1.0e-15);
676                 }
677         }
678     }
679 
680     @Test
681     public void testTidalCorrection1996() {
682 
683         // the reference value has been computed using interp_old.f from
684         // ftp://hpiers.obspm.fr/iers/models/ with the following driver program:
685         //
686         //        PROGRAM OR_TEST
687         //        IMPLICIT NONE
688         //        INTEGER J, N
689         //        PARAMETER (N = 4)
690         //        INTEGER MJD(N)
691         //        DOUBLE PRECISION TTTAI, TAIUTC
692         //        PARAMETER (TAIUTC = 32.000D0)
693         //        PARAMETER (TTTAI  = 32.184D0)
694         //        DOUBLE PRECISION RJD(N), X(N), Y(N), T(N), H(5)
695         //        DOUBLE PRECISION RJDINT, XINT, YINT, TINT, CORX, CORY, CORT
696         //        DATA(MJD(J),    T(J),        X(J),      Y(J),
697         //       &     J=1,4)/
698         //       &    52653, -0.2979055D0, -0.120344D0, 0.217095D0,
699         //       &    52654, -0.2984238D0, -0.121680D0, 0.219400D0,
700         //       &    52655, -0.2987682D0, -0.122915D0, 0.221760D0,
701         //       &    52656, -0.2989957D0, -0.124248D0, 0.224294D0/
702         //  C
703         //        DATA(H(J),J=1,5)/0.0D0, 3600.0D0, 7200.0D0, 43200.0D0, 86400.0D0/
704         //  C
705         //        DO 10 J = 1, N
706         //            RJD(J) = MJD(J) + (TTTAI + TAIUTC) / 86400.0D0
707         //   10   CONTINUE
708         //  C
709         //        DO 20 J = 1, 5
710         //            RJDINT = RJD(2) + H(J) / 86400.0D0
711         //            CALL INTERP(RJD,X,Y,T,N,RJDINT,XINT,YINT,TINT)
712         //            WRITE(6, 30) H(J), TINT, XINT, YINT
713         //   20   CONTINUE
714         //   30   FORMAT(F7.1,3(1X, F20.17))
715         //  C
716         //        END PROGRAM
717         //
718         // the output of this test reads:
719         //           0.0 -0.29839889612705234 -0.12191552977388567  0.21921143351558756
720         //        3600.0 -0.29841696994736727 -0.12208276003864491  0.21925416583607277
721         //        7200.0 -0.29843402412052122 -0.12218102086455683  0.21930263880545320
722         //       43200.0 -0.29866785146390035 -0.12250027826538630  0.22103779809979979
723         //       86400.0 -0.29874248853173840 -0.12308592577174847  0.22161565557764881
724 
725         Utils.setLoaders(IERSConventions.IERS_1996,
726                          Utils.buildEOPList(IERSConventions.IERS_1996, ITRFVersion.ITRF_2008, new double[][] {
727                              {  52653,  -0.2979055,   0.0005744,  -0.120344,   0.217095, 0.0, 0.0, 0.0, 0.0 },
728                              {  52654,  -0.2984238,   0.0004224,  -0.121680,   0.219400, 0.0, 0.0, 0.0, 0.0 },
729                              {  52655,  -0.2987682,   0.0002878,  -0.122915,   0.221760, 0.0, 0.0, 0.0, 0.0 },
730                              {  52656,  -0.2989957,   0.0001778,  -0.124248,   0.224294, 0.0, 0.0, 0.0, 0.0 }
731                          }));
732         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_1996, false);
733 
734         final AbsoluteDate t0 = new AbsoluteDate(new DateComponents(DateComponents.MODIFIED_JULIAN_EPOCH, 52654),
735                                            TimeScalesFactory.getUTC());
736         for (double[] row : new double[][] {
737             {     0.0, -0.29839889612705234, -0.12191552977388567,  0.21921143351558756 },
738             {  3600.0, -0.29841696994736727, -0.12208276003864491,  0.21925416583607277 },
739             {  7200.0, -0.29843402412052122, -0.12218102086455683,  0.21930263880545320 },
740             { 43200.0, -0.29866785146390035, -0.12250027826538630,  0.22103779809979979 },
741             { 86400.0, -0.29874248853173840, -0.12308592577174847,  0.22161565557764881 }
742         }) {
743             AbsoluteDate date = t0.shiftedBy(row[0]);
744             Assert.assertEquals(row[1], eopHistory.getUT1MinusUTC(date), 8.8e-11);
745             Assert.assertEquals(row[2] * Constants.ARC_SECONDS_TO_RADIANS,
746                                 eopHistory.getPoleCorrection(date).getXp(),
747                                 3.2e-14);
748             Assert.assertEquals(row[3] * Constants.ARC_SECONDS_TO_RADIANS,
749                                 eopHistory.getPoleCorrection(date).getYp(),
750                                 8.2e-15);
751         }
752 
753     }
754 
755     @Test
756     public void testTidalCorrection2003() {
757 
758         // the reference value has been computed using the September 2007 version
759         // of interp.f from ftp://hpiers.obspm.fr/iers/models/ adding input/output
760         // parameters for LOD (which was already computed in the underlying routines),
761         // with the following driver program:
762         //
763         //        PROGRAM OR_TEST
764         //        IMPLICIT NONE
765         //        INTEGER J, N
766         //        PARAMETER (N = 4)
767         //        INTEGER MJD(N)
768         //        DOUBLE PRECISION TTTAI, TAIUTC
769         //        PARAMETER (TAIUTC = 32.000D0)
770         //        PARAMETER (TTTAI  = 32.184D0)
771         //        DOUBLE PRECISION RJD(N), UT1(N), LOD(N), X(N), Y(N), H(5)
772         //        DOUBLE PRECISION RJD_INT, UT1_INT, LOD_INT, X_INT, Y_INT
773         //        DATA(MJD(J),   UT1(J),     LOD(J),      X(J),      Y(J),
774         //       &     J=1,4)/
775         //       &    52653, -0.2979055D0, 0.0005744D0, -0.120344D0, 0.217095D0,
776         //       &    52654, -0.2984238D0, 0.0004224D0, -0.121680D0, 0.219400D0,
777         //       &    52655, -0.2987682D0, 0.0002878D0, -0.122915D0, 0.221760D0,
778         //       &    52656, -0.2989957D0, 0.0001778D0, -0.124248D0, 0.224294D0/
779         //  C
780         //        DATA(H(J),J=1,5)/0.0D0, 3600.0D0, 7200.0D0, 43200.0D0, 86400.0D0/
781         //  C
782         //        DO 10 J = 1, N
783         //            RJD(J) = MJD(J) + (TTTAI + TAIUTC) / 86400.0D0
784         //   10   CONTINUE
785         //  C
786         //        DO 20 J = 1, 5
787         //            RJD_INT = RJD(2) + H(J) / 86400.0D0
788         //            CALL INTERP(RJD,X,Y,UT1,LOD,N,
789         //       &                RJD_INT,X_INT,Y_INT,UT1_INT,LOD_INT)
790         //            WRITE(6, 30) H(J),UT1_INT,LOD_INT,X_INT,Y_INT
791         //   20   CONTINUE
792         //   30   FORMAT(F7.1,4(1X, F20.17))
793         //  C
794         //        END PROGRAM
795         //
796         // the output of this test reads:
797         //           0.0 -0.29840026968370659  0.00045312852893139 -0.12196223480123573  0.21922730818562719
798         //        3600.0 -0.29841834564816189  0.00041710864863793 -0.12213345007640604  0.21927433626001305
799         //        7200.0 -0.29843503870494986  0.00039207574457087 -0.12222881007999241  0.21932415788122142
800         //       43200.0 -0.29866930257052676  0.00042895046506082 -0.12247697694276605  0.22105450666130921
801         //       86400.0 -0.29874235341010519  0.00035460263868306 -0.12312252389660779  0.22161364352515728
802 
803         Utils.setLoaders(IERSConventions.IERS_2003,
804                          Utils.buildEOPList(IERSConventions.IERS_2003, ITRFVersion.ITRF_2008, new double[][] {
805                              {  52653,  -0.2979055,   0.0005744,  -0.120344,   0.217095, 0.0, 0.0, 0.0, 0.0 },
806                              {  52654,  -0.2984238,   0.0004224,  -0.121680,   0.219400, 0.0, 0.0, 0.0, 0.0 },
807                              {  52655,  -0.2987682,   0.0002878,  -0.122915,   0.221760, 0.0, 0.0, 0.0, 0.0 },
808                              {  52656,  -0.2989957,   0.0001778,  -0.124248,   0.224294, 0.0, 0.0, 0.0, 0.0 }
809                          }));
810         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_2003, false);
811 
812         final AbsoluteDate t0 = new AbsoluteDate(new DateComponents(DateComponents.MODIFIED_JULIAN_EPOCH, 52654),
813                                                  TimeScalesFactory.getUTC());
814 
815         for (double[] row : new double[][] {
816             {     0.0, -0.29840026968370659,  0.00045312852893139, -0.12196223480123573,  0.21922730818562719 },
817             {  3600.0, -0.29841834564816189,  0.00041710864863793, -0.12213345007640604,  0.21927433626001305 },
818             {  7200.0, -0.29843503870494986,  0.00039207574457087, -0.12222881007999241,  0.21932415788122142 },
819             { 43200.0, -0.29866930257052676,  0.00042895046506082, -0.12247697694276605,  0.22105450666130921 },
820             { 86400.0, -0.29874235341010519,  0.00035460263868306, -0.12312252389660779,  0.22161364352515728 }
821         }) {
822             AbsoluteDate date = t0.shiftedBy(row[0]);
823             Assert.assertEquals(row[1], eopHistory.getUT1MinusUTC(date), 4.3e-8);
824             Assert.assertEquals(row[2], eopHistory.getLOD(date), 1.4e-7);
825             Assert.assertEquals(row[3] * Constants.ARC_SECONDS_TO_RADIANS,
826                                 eopHistory.getPoleCorrection(date).getXp(),
827                                 1.6e-10);
828             Assert.assertEquals(row[4] * Constants.ARC_SECONDS_TO_RADIANS,
829                                 eopHistory.getPoleCorrection(date).getYp(),
830                                 0.7e-10);
831         }
832 
833     }
834 
835     @Test
836     public void testTidalCorrection2010() {
837 
838         // the reference value has been computed using the September 2007 version
839         // of interp.f from ftp://hpiers.obspm.fr/iers/models/ adding input/output
840         // parameters for LOD (which was already computed in the underlying routines),
841         // with the following driver program:
842         //
843         //        PROGRAM OR_TEST
844         //        IMPLICIT NONE
845         //        INTEGER J, N
846         //        PARAMETER (N = 4)
847         //        INTEGER MJD(N)
848         //        DOUBLE PRECISION TTTAI, TAIUTC
849         //        PARAMETER (TAIUTC = 32.000D0)
850         //        PARAMETER (TTTAI  = 32.184D0)
851         //        DOUBLE PRECISION RJD(N), UT1(N), LOD(N), X(N), Y(N), H(5)
852         //        DOUBLE PRECISION RJD_INT, UT1_INT, LOD_INT, X_INT, Y_INT
853         //        DATA(MJD(J),   UT1(J),     LOD(J),      X(J),      Y(J),
854         //       &     J=1,4)/
855         //       &    52653, -0.2979055D0, 0.0005744D0, -0.120344D0, 0.217095D0,
856         //       &    52654, -0.2984238D0, 0.0004224D0, -0.121680D0, 0.219400D0,
857         //       &    52655, -0.2987682D0, 0.0002878D0, -0.122915D0, 0.221760D0,
858         //       &    52656, -0.2989957D0, 0.0001778D0, -0.124248D0, 0.224294D0/
859         //  C
860         //        DATA(H(J),J=1,5)/0.0D0, 3600.0D0, 7200.0D0, 43200.0D0, 86400.0D0/
861         //  C
862         //        DO 10 J = 1, N
863         //            RJD(J) = MJD(J) + (TTTAI + TAIUTC) / 86400.0D0
864         //   10   CONTINUE
865         //  C
866         //        DO 20 J = 1, 5
867         //            RJD_INT = RJD(2) + H(J) / 86400.0D0
868         //            CALL INTERP(RJD,X,Y,UT1,LOD,N,
869         //       &                RJD_INT,X_INT,Y_INT,UT1_INT,LOD_INT)
870         //            WRITE(6, 30) H(J),UT1_INT,LOD_INT,X_INT,Y_INT
871         //   20   CONTINUE
872         //   30   FORMAT(F7.1,4(1X, F20.17))
873         //  C
874         //        END PROGRAM
875         //
876         // the output of this test reads:
877         //           0.0 -0.29840026968370659  0.00045312852893139 -0.12196223480123573  0.21922730818562719
878         //        3600.0 -0.29841834564816189  0.00041710864863793 -0.12213345007640604  0.21927433626001305
879         //        7200.0 -0.29843503870494986  0.00039207574457087 -0.12222881007999241  0.21932415788122142
880         //       43200.0 -0.29866930257052676  0.00042895046506082 -0.12247697694276605  0.22105450666130921
881         //       86400.0 -0.29874235341010519  0.00035460263868306 -0.12312252389660779  0.22161364352515728
882 
883         Utils.setLoaders(IERSConventions.IERS_2010,
884                          Utils.buildEOPList(IERSConventions.IERS_2010, ITRFVersion.ITRF_2008, new double[][] {
885                              {  52653,  -0.2979055,   0.0005744,  -0.120344,   0.217095, 0.0, 0.0, 0.0, 0.0 },
886                              {  52654,  -0.2984238,   0.0004224,  -0.121680,   0.219400, 0.0, 0.0, 0.0, 0.0 },
887                              {  52655,  -0.2987682,   0.0002878,  -0.122915,   0.221760, 0.0, 0.0, 0.0, 0.0 },
888                              {  52656,  -0.2989957,   0.0001778,  -0.124248,   0.224294, 0.0, 0.0, 0.0, 0.0 }
889                          }));
890         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_2010, false);
891 
892         final AbsoluteDate t0 = new AbsoluteDate(new DateComponents(DateComponents.MODIFIED_JULIAN_EPOCH, 52654),
893                                                  TimeScalesFactory.getUTC());
894 
895         for (double[] row : new double[][] {
896             {     0.0, -0.29840026968370659,  0.00045312852893139, -0.12196223480123573,  0.21922730818562719 },
897             {  3600.0, -0.29841834564816189,  0.00041710864863793, -0.12213345007640604,  0.21927433626001305 },
898             {  7200.0, -0.29843503870494986,  0.00039207574457087, -0.12222881007999241,  0.21932415788122142 },
899             { 43200.0, -0.29866930257052676,  0.00042895046506082, -0.12247697694276605,  0.22105450666130921 },
900             { 86400.0, -0.29874235341010519,  0.00035460263868306, -0.12312252389660779,  0.22161364352515728 }
901         }) {
902             AbsoluteDate date = t0.shiftedBy(row[0]);
903             Assert.assertEquals(row[1], eopHistory.getUT1MinusUTC(date), 2.5e-12);
904             Assert.assertEquals(row[2], eopHistory.getLOD(date), 4.3e-11);
905             Assert.assertEquals(row[3] * Constants.ARC_SECONDS_TO_RADIANS,
906                                 eopHistory.getPoleCorrection(date).getXp(),
907                                 1.6e-10);
908             Assert.assertEquals(row[4] * Constants.ARC_SECONDS_TO_RADIANS,
909                                 eopHistory.getPoleCorrection(date).getYp(),
910                                 0.7e-10);
911         }
912 
913     }
914 
915     @Test
916     public void testGMSTRate1996() {
917         final UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_1996, true);
918         checkGMSTRate(IERSConventions.IERS_1996.getGMSTFunction(ut1),
919                       IERSConventions.IERS_1996.getGMSTRateFunction(ut1),
920                       AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
921                       0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 8.0e-13);
922     }
923 
924     @Test
925     public void testGMSTRate2003() {
926         final UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2003, true);
927         checkGMSTRate(IERSConventions.IERS_2003.getGMSTFunction(ut1),
928                       IERSConventions.IERS_2003.getGMSTRateFunction(ut1),
929                       AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
930                       0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 8.0e-13);
931     }
932 
933     @Test
934     public void testGMSTRate2010() {
935         final UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
936         checkGMSTRate(IERSConventions.IERS_2010.getGMSTFunction(ut1),
937                       IERSConventions.IERS_2010.getGMSTRateFunction(ut1),
938                       AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
939                       0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 8.0e-13);
940     }
941 
942     private void checkGMSTRate(final TimeScalarFunction gmst, final TimeScalarFunction gmstRate,
943                                final AbsoluteDate date, final double span, final double sampleStep,
944                                final double h, final double tolerance)
945         {
946         UnivariateDifferentiableFunction differentiated =
947                         new FiniteDifferencesDifferentiator(4, h).differentiate(new UnivariateFunction() {
948                             @Override
949                             public double value(final double dt) {
950                                 return gmst.value(date.shiftedBy(dt));
951                             }
952                         });
953         DSFactory factory = new DSFactory(1, 1);
954         double maxRateError = 0;
955         for (double dt = 0; dt < span; dt += sampleStep) {
956             double rateRef = differentiated.value(factory.variable(0, dt)).getPartialDerivative(1);
957             double rate    = gmstRate.value(date.shiftedBy(dt));
958             maxRateError   = FastMath.max(maxRateError, FastMath.abs(rateRef - rate));
959         }
960         Assert.assertEquals(0, maxRateError, tolerance);
961     }
962 
963     @Test
964     public void testFieldGMSTRate1996() {
965         final UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_1996, true);
966         checkFieldGMSTRate(IERSConventions.IERS_1996.getGMSTFunction(ut1),
967                            IERSConventions.IERS_1996.getGMSTRateFunction(ut1),
968                            AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
969                            0.8 * Constants.JULIAN_DAY, 600.0, 10.0, Double.MIN_VALUE);
970     }
971 
972     @Test
973     public void testFieldGMSTRate2003() {
974         final UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2003, true);
975         checkFieldGMSTRate(IERSConventions.IERS_2003.getGMSTFunction(ut1),
976                            IERSConventions.IERS_2003.getGMSTRateFunction(ut1),
977                            AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
978                            0.8 * Constants.JULIAN_DAY, 600.0, 10.0, Double.MIN_VALUE);
979     }
980 
981     @Test
982     public void testFieldGMSTRate2010() {
983         final UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
984         checkFieldGMSTRate(IERSConventions.IERS_2010.getGMSTFunction(ut1),
985                            IERSConventions.IERS_2010.getGMSTRateFunction(ut1),
986                            AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
987                            0.8 * Constants.JULIAN_DAY, 600.0, 10.0, Double.MIN_VALUE);
988     }
989 
990     private void checkFieldGMSTRate(final TimeScalarFunction gmst, final TimeScalarFunction gmstRate,
991                                     final AbsoluteDate date, final double span, final double sampleStep,
992                                     final double h, final double tolerance)
993         {
994         double maxError = 0;
995         for (double dt = 0; dt < span; dt += sampleStep) {
996             double rateDouble    = gmstRate.value(date.shiftedBy(dt));
997             double rateDecimal64 = gmstRate.value(new FieldAbsoluteDate<>(date, new Decimal64(dt))).doubleValue();
998             maxError = FastMath.max(maxError, FastMath.abs(rateDouble - rateDecimal64));
999         }
1000         Assert.assertEquals(0, maxError, tolerance);
1001     }
1002 
1003     @Test
1004     public void testPrecessionFunction1996FieldConsistency() {
1005         checkVectorFunctionConsistency(IERSConventions.IERS_1996.getPrecessionFunction(), 3,
1006                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1007                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 2.0e-22, 3.0e-19);
1008     }
1009 
1010     @Test
1011     public void testPrecessionFunction2003FieldConsistency() {
1012         checkVectorFunctionConsistency(IERSConventions.IERS_2003.getPrecessionFunction(), 3,
1013                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1014                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 6.0e-17, 6.0e-18);
1015     }
1016 
1017     @Test
1018     public void testPrecessionFunction2010FieldConsistency() {
1019         checkVectorFunctionConsistency(IERSConventions.IERS_2010.getPrecessionFunction(), 3,
1020                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1021                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 6.0e-17, 6.0e-18);
1022     }
1023 
1024     @Test
1025     public void testNutationFunction1996FieldConsistency() {
1026         checkVectorFunctionConsistency(IERSConventions.IERS_1996.getNutationFunction(), 3,
1027                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1028                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 2.0e-19, 4.0e-21);
1029     }
1030 
1031     @Test
1032     public void testNutationFunction2003FieldConsistency() {
1033         checkVectorFunctionConsistency(IERSConventions.IERS_2003.getNutationFunction(), 3,
1034                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1035                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 3.0e-19, 6.0e-21);
1036     }
1037 
1038     @Test
1039     public void testNutationFunction2010FieldConsistency() {
1040         checkVectorFunctionConsistency(IERSConventions.IERS_2010.getNutationFunction(), 3,
1041                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1042                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 5.0e-19, 6.0e-21);
1043     }
1044 
1045     @Test
1046     public void testXYSpXY2Function1996FieldConsistency() {
1047         checkVectorFunctionConsistency(IERSConventions.IERS_1996.getXYSpXY2Function(), 3,
1048                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1049                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 6.0e-20, 2.0e-21);
1050     }
1051 
1052     @Test
1053     public void testXYSpXY2Function2003FieldConsistency() {
1054         checkVectorFunctionConsistency(IERSConventions.IERS_2003.getXYSpXY2Function(), 3,
1055                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1056                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 2.0e-19, 3.0e-21);
1057     }
1058 
1059     @Test
1060     public void testXYSpXY2Function2010FieldConsistency() {
1061         checkVectorFunctionConsistency(IERSConventions.IERS_2010.getXYSpXY2Function(), 3,
1062                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1063                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 2.0e-19, 3.0e-21);
1064     }
1065 
1066     @Test
1067     public void testMeanObliquityFunction1996FieldConsistency() {
1068         checkScalarFunctionConsistency(IERSConventions.IERS_1996.getMeanObliquityFunction(),
1069                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1070                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 6.0e-17, 6.0e-18);
1071     }
1072 
1073     @Test
1074     public void testMeanObliquityFunction2003FieldConsistency() {
1075         checkScalarFunctionConsistency(IERSConventions.IERS_2003.getMeanObliquityFunction(),
1076                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1077                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 6.0e-17, 6.0e-18);
1078     }
1079 
1080     @Test
1081     public void testMeanObliquityFunction2010FieldConsistency() {
1082         checkScalarFunctionConsistency(IERSConventions.IERS_2010.getMeanObliquityFunction(),
1083                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1084                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 6.0e-17, 6.0e-18);
1085     }
1086 
1087     @Test
1088     public void testEOPTidalCorrection1996FieldConsistency() {
1089         checkVectorFunctionConsistency(IERSConventions.IERS_1996.getEOPTidalCorrection(), 4,
1090                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1091                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 2.0e-17, 6.0e-20);
1092     }
1093 
1094     @Test
1095     public void testEOPTidalCorrection2003FieldConsistency() {
1096         checkVectorFunctionConsistency(IERSConventions.IERS_2003.getEOPTidalCorrection(), 4,
1097                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1098                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 2.0e-17, 2.0e-19);
1099     }
1100 
1101     @Test
1102     public void testEOPTidalCorrection2010FieldConsistency() {
1103         checkVectorFunctionConsistency(IERSConventions.IERS_2010.getEOPTidalCorrection(), 4,
1104                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1105                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 2.0e-17, 2.0e-19);
1106     }
1107 
1108     @Test
1109     public void testTideFrequencyDependenceFunction1996FieldConsistency() {
1110         final UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_1996, true);
1111         checkVectorFunctionConsistency(IERSConventions.IERS_1996.getTideFrequencyDependenceFunction(ut1), 5,
1112                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1113                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 4.0e-24, 4.0e-22);
1114     }
1115 
1116     @Test
1117     public void testTideFrequencyDependenceFunction2003FieldConsistency() {
1118         final UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2003, true);
1119         checkVectorFunctionConsistency(IERSConventions.IERS_2003.getTideFrequencyDependenceFunction(ut1), 5,
1120                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1121                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 4.0e-24, 4.0e-22);
1122     }
1123 
1124     @Test
1125     public void testTideFrequencyDependenceFunction2010FieldConsistency() {
1126         final UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
1127         checkVectorFunctionConsistency(IERSConventions.IERS_2010.getTideFrequencyDependenceFunction(ut1), 5,
1128                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1129                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 4.0e-24, 4.0e-22);
1130     }
1131 
1132     @Test
1133     public void testSolidPoleTide1996FieldConsistency() {
1134         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_1996, false);
1135         checkVectorFunctionConsistency(IERSConventions.IERS_1996.getSolidPoleTide(eopHistory), 2,
1136                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1137                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 3.0e-25, 4.0e-26);
1138     }
1139 
1140     @Test
1141     public void testSolidPoleTide2003FieldConsistency() {
1142         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_2003, false);
1143         checkVectorFunctionConsistency(IERSConventions.IERS_2003.getSolidPoleTide(eopHistory), 2,
1144                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1145                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 2.0e-25, 3.0e-26);
1146     }
1147 
1148     @Test
1149     public void testSolidPoleTide2010FieldConsistency() {
1150         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_2010, false);
1151         checkVectorFunctionConsistency(IERSConventions.IERS_2010.getSolidPoleTide(eopHistory), 2,
1152                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1153                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 2.0e-25, 2.0e-26);
1154     }
1155 
1156     @Test
1157     public void testOceanPoleTide1996FieldConsistency() {
1158         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_1996, false);
1159         checkVectorFunctionConsistency(IERSConventions.IERS_1996.getOceanPoleTide(eopHistory), 2,
1160                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1161                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, Double.MIN_VALUE, Double.MIN_VALUE);
1162     }
1163 
1164     @Test
1165     public void testOceanPoleTide2003FieldConsistency() {
1166         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_2003, false);
1167         checkVectorFunctionConsistency(IERSConventions.IERS_2003.getOceanPoleTide(eopHistory), 2,
1168                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1169                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, Double.MIN_VALUE, Double.MIN_VALUE);
1170     }
1171 
1172     @Test
1173     public void testOceanPoleTide2010FieldConsistency() {
1174         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_2010, false);
1175         checkVectorFunctionConsistency(IERSConventions.IERS_2010.getOceanPoleTide(eopHistory), 2,
1176                                        AbsoluteDate.J2000_EPOCH.shiftedBy(-0.4 * Constants.JULIAN_DAY),
1177                                        0.8 * Constants.JULIAN_DAY, 600.0, 10.0, 2.0e-26, 2.0e-27);
1178     }
1179 
1180     private void checkScalarFunctionConsistency(final TimeScalarFunction function,
1181                                                 final AbsoluteDate date, final double span, final double sampleStep,
1182                                                 final double h, final double valueTolerance, final double derivativeTolerance) {
1183 
1184         UnivariateDifferentiableFunction differentiated =
1185                 new FiniteDifferencesDifferentiator(4, h).differentiate(new UnivariateFunction() {
1186                     @Override
1187                     public double value(final double dt) {
1188                         return function.value(date.shiftedBy(dt));
1189                     }
1190                 });
1191 
1192         DSFactory factory = new DSFactory(1, 1);
1193         FieldAbsoluteDate<DerivativeStructure> dsDate = new FieldAbsoluteDate<>(date, factory.constant(0.0));
1194         double maxValueError = 0;
1195         double maxDerivativeError = 0;
1196         for (double dt = 0; dt < span; dt += sampleStep) {
1197             DerivativeStructure dsdt = factory.variable(0, dt);
1198             DerivativeStructure yRef = differentiated.value(dsdt);
1199             DerivativeStructure y    = function.value(dsDate.shiftedBy(dsdt));
1200             maxValueError      = FastMath.max(maxValueError,
1201                                               FastMath.abs(yRef.getValue() - y.getValue()));
1202             maxDerivativeError = FastMath.max(maxDerivativeError,
1203                                               FastMath.abs(yRef.getPartialDerivative(1) - y.getPartialDerivative(1)));
1204         }
1205         Assert.assertEquals(0, maxValueError,      valueTolerance);
1206         Assert.assertEquals(0, maxDerivativeError, derivativeTolerance);
1207 
1208     }
1209 
1210     private void checkVectorFunctionConsistency(final TimeVectorFunction function, final int dimension,
1211                                                 final AbsoluteDate date, final double span, final double sampleStep,
1212                                                 final double h, final double valueTolerance, final double derivativeTolerance) {
1213 
1214         UnivariateDifferentiableVectorFunction differentiated =
1215                 new FiniteDifferencesDifferentiator(4, h).differentiate(new UnivariateVectorFunction() {
1216                     @Override
1217                     public double[] value(final double dt) {
1218                         return function.value(date.shiftedBy(dt));
1219                     }
1220                 });
1221 
1222         DSFactory factory = new DSFactory(1, 1);
1223         FieldAbsoluteDate<DerivativeStructure> dsDate = new FieldAbsoluteDate<>(date, factory.constant(0.0));
1224         double maxValueError = 0;
1225         double maxDerivativeError = 0;
1226         for (double dt = 0; dt < span; dt += sampleStep) {
1227             DerivativeStructure dsdt = factory.variable(0, dt);
1228             DerivativeStructure[] yRef = differentiated.value(dsdt);
1229             DerivativeStructure[] y    = function.value(dsDate.shiftedBy(dsdt));
1230             Assert.assertEquals(dimension, yRef.length);
1231             Assert.assertEquals(dimension, y.length);
1232             for (int i = 0; i < dimension; ++i) {
1233                 maxValueError      = FastMath.max(maxValueError,
1234                                                   FastMath.abs(yRef[i].getValue() - y[i].getValue()));
1235                 maxDerivativeError = FastMath.max(maxDerivativeError,
1236                                                   FastMath.abs(yRef[i].getPartialDerivative(1) - y[i].getPartialDerivative(1)));
1237             }
1238         }
1239         Assert.assertEquals(0, maxValueError,      valueTolerance);
1240         Assert.assertEquals(0, maxDerivativeError, derivativeTolerance);
1241 
1242     }
1243 
1244     @Before
1245     public void setUp() {
1246         Utils.setDataRoot("compressed-data");
1247     }
1248 
1249 }