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