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.frames;
18  
19  import java.util.ArrayList;
20  import java.util.List;
21  import java.util.SortedSet;
22  import java.util.TreeSet;
23  import java.util.stream.Collectors;
24  
25  import org.hipparchus.CalculusFieldElement;
26  import org.hipparchus.analysis.UnivariateVectorFunction;
27  import org.hipparchus.analysis.differentiation.DSFactory;
28  import org.hipparchus.analysis.differentiation.DerivativeStructure;
29  import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
30  import org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction;
31  import org.hipparchus.geometry.euclidean.threed.Rotation;
32  import org.hipparchus.geometry.euclidean.threed.Vector3D;
33  import org.hipparchus.util.Binary64;
34  import org.hipparchus.util.Binary64Field;
35  import org.hipparchus.util.FastMath;
36  import org.hipparchus.util.MathArrays;
37  import org.hipparchus.util.MathUtils;
38  import org.junit.jupiter.api.Assertions;
39  import org.junit.jupiter.api.BeforeEach;
40  import org.junit.jupiter.api.Test;
41  import org.orekit.Utils;
42  import org.orekit.bodies.CelestialBodyFactory;
43  import org.orekit.bodies.GeodeticPoint;
44  import org.orekit.bodies.OneAxisEllipsoid;
45  import org.orekit.data.DataContext;
46  import org.orekit.data.DataProvidersManager;
47  import org.orekit.errors.OrekitException;
48  import org.orekit.errors.OrekitIllegalStateException;
49  import org.orekit.errors.OrekitMessages;
50  import org.orekit.time.AbsoluteDate;
51  import org.orekit.time.ChronologicalComparator;
52  import org.orekit.time.DateComponents;
53  import org.orekit.time.FieldAbsoluteDate;
54  import org.orekit.time.TimeComponents;
55  import org.orekit.time.TimeScalesFactory;
56  import org.orekit.time.UTCScale;
57  import org.orekit.utils.AngularCoordinates;
58  import org.orekit.utils.AngularDerivativesFilter;
59  import org.orekit.utils.CartesianDerivativesFilter;
60  import org.orekit.utils.Constants;
61  import org.orekit.utils.IERSConventions;
62  import org.orekit.utils.IERSConventions.NutationCorrectionConverter;
63  import org.orekit.utils.PVCoordinates;
64  
65  public class FramesFactoryTest {
66  
67      @Test
68      public void testTreeRoot() {
69          Assertions.assertNull(FramesFactory.getFrame(Predefined.GCRF).getParent());
70      }
71  
72      @Test
73      public void testWrongSupportedFileNames1980() {
74          FramesFactory.addDefaultEOP1980HistoryLoaders("wrong-rapidDataColumns-1980",
75                                                        "wrong-rapidDataXML-1980",
76                                                        "wrong-eopC04-1980",
77                                                        "wrong-bulletinB-1980",
78                                                        "wrong-bulletinA-1980",
79                                                        "wrong-csv");
80          try {
81              FramesFactory.getEOPHistory(IERSConventions.IERS_1996, true).getStartDate();
82              Assertions.fail("an exception should have been thrown");
83          } catch (OrekitIllegalStateException oe) {
84              Assertions.assertEquals(OrekitMessages.NO_CACHED_ENTRIES, oe.getSpecifier());
85          }
86      }
87  
88      @Test
89      public void testWrongSupportedFileNames2000() {
90          FramesFactory.addDefaultEOP2000HistoryLoaders("wrong-rapidDataColumns-2000",
91                                                        "wrong-rapidDataXML-2000",
92                                                        "wrong-eopC04-2000",
93                                                        "wrong-bulletinB-2000",
94                                                        "wrong-bulletinA-2000",
95                                                        "wrong-csv");
96          try {
97              FramesFactory.getEOPHistory(IERSConventions.IERS_2010, true).getStartDate();
98              Assertions.fail("an exception should have been thrown");
99          } catch (OrekitIllegalStateException oe) {
100             Assertions.assertEquals(OrekitMessages.NO_CACHED_ENTRIES, oe.getSpecifier());
101         }
102     }
103 
104     @Test
105     public void testWrongConventions() {
106         // set up only 1980 conventions
107         FramesFactory.addDefaultEOP1980HistoryLoaders(null, null, null, null, null, null);
108         try {
109             // attempt to retrieve 2000 conventions
110             FramesFactory.getEOPHistory(IERSConventions.IERS_2010, true).getStartDate();
111             Assertions.fail("an exception should have been thrown");
112         } catch (OrekitIllegalStateException oe) {
113             Assertions.assertEquals(OrekitMessages.NO_CACHED_ENTRIES, oe.getSpecifier());
114         }
115     }
116 
117     @Test
118     public void testEOPLoaderException() {
119         final boolean[] flags = new boolean[2];
120         try {
121             FramesFactory.addEOPHistoryLoader(IERSConventions.IERS_2010, new EopHistoryLoader() {
122                 @Override
123                 public void fillHistory(NutationCorrectionConverter converter, SortedSet<EOPEntry> history) {
124                     // don't really fill history here
125                     flags[0] = true;
126                 }
127             });
128             FramesFactory.addEOPHistoryLoader(IERSConventions.IERS_2010, new EopHistoryLoader() {
129                 @Override
130                 public void fillHistory(NutationCorrectionConverter converter, SortedSet<EOPEntry> history)
131                     {
132                     // generate exception
133                     flags[1] = true;
134                     throw new OrekitException(OrekitMessages.NO_DATA_GENERATED, AbsoluteDate.J2000_EPOCH);
135                 }
136             });
137             FramesFactory.getEOPHistory(IERSConventions.IERS_2010, true);
138             Assertions.fail("an exception should have been thrown");
139         } catch (OrekitException oe) {
140             Assertions.assertTrue(flags[0]);
141             Assertions.assertTrue(flags[1]);
142             Assertions.assertEquals(OrekitMessages.NO_DATA_GENERATED, oe.getSpecifier());
143         }
144     }
145 
146     @Test
147     public void testUnwrapInterpolatingTransformProvider() {
148         TransformProvider raw = new TransformProvider() {
149             public Transform getTransform(final AbsoluteDate date) {
150                 double dt = date.durationFrom(AbsoluteDate.J2000_EPOCH);
151                 double sin = FastMath.sin(dt * MathUtils.TWO_PI / Constants.JULIAN_DAY);
152                 return new Transform(date,
153                                      new PVCoordinates(new Vector3D(sin, Vector3D.PLUS_I),
154                                                        Vector3D.ZERO));
155             }
156             public <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
157                 throw new UnsupportedOperationException("never called in this test");
158             }
159         };
160         Frame parent = FramesFactory.getGCRF();
161         Assertions.assertEquals("GCRF", parent.getName());
162         Frame frame  = new Frame(parent,
163                                  new InterpolatingTransformProvider(raw,
164                                                                     CartesianDerivativesFilter.USE_P,
165                                                                     AngularDerivativesFilter.USE_R,
166                                                                     4, Constants.JULIAN_DAY, 10,
167                                                                     Constants.JULIAN_YEAR, 2 * Constants.JULIAN_DAY),
168                                  "sine");
169         double maxErrorNonInterpolating = 0;
170         double maxErrorInterpolating    = 0;
171         for (double dt = 0; dt < Constants.JULIAN_DAY; dt += 60.0) {
172             AbsoluteDate date            = AbsoluteDate.J2000_EPOCH.shiftedBy(dt);
173             Transform reference          = raw.getTransform(date);
174             Transform nonInterpolating   = FramesFactory.getNonInterpolatingTransform(parent, frame, date);
175             Transform interpolating      = parent.getTransformTo(frame, date);
176             double errorNonInterpolating = Vector3D.distance(reference.getTranslation(),
177                                                              nonInterpolating.getTranslation());
178             maxErrorNonInterpolating     = FastMath.max(maxErrorNonInterpolating, errorNonInterpolating);
179             double errorInterpolating    = Vector3D.distance(reference.getTranslation(),
180                                                              interpolating.getTranslation());
181             maxErrorInterpolating        = FastMath.max(maxErrorInterpolating, errorInterpolating);
182         }
183         Assertions.assertEquals(0.0, maxErrorNonInterpolating, 1.0e-15);
184         Assertions.assertEquals(1.0, maxErrorInterpolating,    1.0e-15);
185     }
186 
187     @Test
188     public void testUnwrapShiftingTransformProvider() {
189         TransformProvider raw = new TransformProvider() {
190             public Transform getTransform(final AbsoluteDate date) {
191                 double dt = date.durationFrom(AbsoluteDate.J2000_EPOCH);
192                 double sin = FastMath.sin(dt * MathUtils.TWO_PI / Constants.JULIAN_DAY);
193                 return new Transform(date,
194                                      new PVCoordinates(new Vector3D(sin, Vector3D.PLUS_I),
195                                                        Vector3D.ZERO));
196             }
197             public <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
198                 throw new UnsupportedOperationException("never called in this test");
199             }
200         };
201         Frame parent = FramesFactory.getGCRF();
202         Frame frame  = new Frame(parent,
203                                  new ShiftingTransformProvider(raw,
204                                                                CartesianDerivativesFilter.USE_P,
205                                                                AngularDerivativesFilter.USE_R,
206                                                                4, Constants.JULIAN_DAY, 10,
207                                                                Constants.JULIAN_YEAR, 2 * Constants.JULIAN_DAY),
208                                  "sine");
209         double maxErrorNonShifting = 0;
210         double maxErrorShifting    = 0;
211         for (double dt = 0; dt < Constants.JULIAN_DAY; dt += 60.0) {
212             AbsoluteDate date       = AbsoluteDate.J2000_EPOCH.shiftedBy(dt);
213             Transform reference     = raw.getTransform(date);
214             Transform nonShifting   = FramesFactory.getNonInterpolatingTransform(parent, frame, date);
215             Transform shifting      = parent.getTransformTo(frame, date);
216             double errorNonShifting = Vector3D.distance(reference.getTranslation(),
217                                                         nonShifting.getTranslation());
218             maxErrorNonShifting     = FastMath.max(maxErrorNonShifting, errorNonShifting);
219             double errorShifting    = Vector3D.distance(reference.getTranslation(),
220                                                         shifting.getTranslation());
221             maxErrorShifting        = FastMath.max(maxErrorShifting, errorShifting);
222         }
223         Assertions.assertEquals(0.0, maxErrorNonShifting, 1.0e-15);
224         Assertions.assertEquals(1.0, maxErrorShifting,    1.0e-15);
225     }
226 
227     @Test
228     public void testTreeICRF() {
229         Frame icrf = FramesFactory.getFrame(Predefined.ICRF);
230         Assertions.assertEquals("ICRF", icrf.getName());
231         Transform t = icrf.getTransformTo(FramesFactory.getGCRF(),
232                                           new AbsoluteDate(1969, 6, 25, TimeScalesFactory.getTT()));
233         Assertions.assertEquals(0.0, t.getRotation().getAngle(), 1.0e-15);
234         Assertions.assertEquals(CelestialBodyFactory.EARTH_MOON + "/inertial", icrf.getParent().getName());
235         Assertions.assertEquals(Predefined.GCRF.getName(), icrf.getParent().getParent().getName());
236     }
237 
238     @Test
239     public void testTree() {
240         Predefined[][] reference = new Predefined[][] {
241             { Predefined.EME2000,                                Predefined.GCRF },
242             { Predefined.ITRF_CIO_CONV_1996_ACCURATE_EOP,        Predefined.TIRF_CONVENTIONS_1996_ACCURATE_EOP },
243             { Predefined.ITRF_CIO_CONV_1996_SIMPLE_EOP,          Predefined.TIRF_CONVENTIONS_1996_SIMPLE_EOP   },
244             { Predefined.ITRF_CIO_CONV_2003_ACCURATE_EOP,        Predefined.TIRF_CONVENTIONS_2003_ACCURATE_EOP },
245             { Predefined.ITRF_CIO_CONV_2003_SIMPLE_EOP,          Predefined.TIRF_CONVENTIONS_2003_SIMPLE_EOP   },
246             { Predefined.ITRF_CIO_CONV_2010_ACCURATE_EOP,        Predefined.TIRF_CONVENTIONS_2010_ACCURATE_EOP },
247             { Predefined.ITRF_CIO_CONV_2010_SIMPLE_EOP,          Predefined.TIRF_CONVENTIONS_2010_SIMPLE_EOP   },
248             { Predefined.TIRF_CONVENTIONS_1996_ACCURATE_EOP,     Predefined.CIRF_CONVENTIONS_1996_ACCURATE_EOP },
249             { Predefined.TIRF_CONVENTIONS_1996_SIMPLE_EOP,       Predefined.CIRF_CONVENTIONS_1996_SIMPLE_EOP   },
250             { Predefined.TIRF_CONVENTIONS_2003_ACCURATE_EOP,     Predefined.CIRF_CONVENTIONS_2003_ACCURATE_EOP },
251             { Predefined.TIRF_CONVENTIONS_2003_SIMPLE_EOP,       Predefined.CIRF_CONVENTIONS_2003_SIMPLE_EOP   },
252             { Predefined.TIRF_CONVENTIONS_2010_ACCURATE_EOP,     Predefined.CIRF_CONVENTIONS_2010_ACCURATE_EOP },
253             { Predefined.TIRF_CONVENTIONS_2010_SIMPLE_EOP,       Predefined.CIRF_CONVENTIONS_2010_SIMPLE_EOP   },
254             { Predefined.CIRF_CONVENTIONS_1996_ACCURATE_EOP,     Predefined.GCRF },
255             { Predefined.CIRF_CONVENTIONS_1996_SIMPLE_EOP,       Predefined.GCRF },
256             { Predefined.CIRF_CONVENTIONS_2003_ACCURATE_EOP,     Predefined.GCRF },
257             { Predefined.CIRF_CONVENTIONS_2003_SIMPLE_EOP,       Predefined.GCRF },
258             { Predefined.CIRF_CONVENTIONS_2010_ACCURATE_EOP,     Predefined.GCRF },
259             { Predefined.CIRF_CONVENTIONS_2010_SIMPLE_EOP,       Predefined.GCRF },
260             { Predefined.VEIS_1950,                              Predefined.GTOD_WITHOUT_EOP_CORRECTIONS },
261             { Predefined.ITRF_EQUINOX_CONV_1996_ACCURATE_EOP,    Predefined.GTOD_CONVENTIONS_1996_ACCURATE_EOP },
262             { Predefined.ITRF_EQUINOX_CONV_1996_SIMPLE_EOP,      Predefined.GTOD_CONVENTIONS_1996_SIMPLE_EOP   },
263             { Predefined.ITRF_EQUINOX_CONV_2003_ACCURATE_EOP,    Predefined.GTOD_CONVENTIONS_2003_ACCURATE_EOP },
264             { Predefined.ITRF_EQUINOX_CONV_2003_SIMPLE_EOP,      Predefined.GTOD_CONVENTIONS_2003_SIMPLE_EOP   },
265             { Predefined.ITRF_EQUINOX_CONV_2010_ACCURATE_EOP,    Predefined.GTOD_CONVENTIONS_2010_ACCURATE_EOP },
266             { Predefined.ITRF_EQUINOX_CONV_2010_SIMPLE_EOP,      Predefined.GTOD_CONVENTIONS_2010_SIMPLE_EOP   },
267             { Predefined.GTOD_WITHOUT_EOP_CORRECTIONS,           Predefined.TOD_WITHOUT_EOP_CORRECTIONS },
268             { Predefined.GTOD_CONVENTIONS_1996_ACCURATE_EOP,     Predefined.TOD_CONVENTIONS_1996_ACCURATE_EOP },
269             { Predefined.GTOD_CONVENTIONS_1996_SIMPLE_EOP,       Predefined.TOD_CONVENTIONS_1996_SIMPLE_EOP   },
270             { Predefined.GTOD_CONVENTIONS_2003_ACCURATE_EOP,     Predefined.TOD_CONVENTIONS_2003_ACCURATE_EOP },
271             { Predefined.GTOD_CONVENTIONS_2003_SIMPLE_EOP,       Predefined.TOD_CONVENTIONS_2003_SIMPLE_EOP   },
272             { Predefined.GTOD_CONVENTIONS_2010_ACCURATE_EOP,     Predefined.TOD_CONVENTIONS_2010_ACCURATE_EOP },
273             { Predefined.GTOD_CONVENTIONS_2010_SIMPLE_EOP,       Predefined.TOD_CONVENTIONS_2010_SIMPLE_EOP   },
274             { Predefined.TOD_WITHOUT_EOP_CORRECTIONS,            Predefined.MOD_WITHOUT_EOP_CORRECTIONS },
275             { Predefined.TOD_CONVENTIONS_1996_ACCURATE_EOP,      Predefined.MOD_CONVENTIONS_1996 },
276             { Predefined.TOD_CONVENTIONS_1996_SIMPLE_EOP,        Predefined.MOD_CONVENTIONS_1996 },
277             { Predefined.TOD_CONVENTIONS_2003_ACCURATE_EOP,      Predefined.MOD_CONVENTIONS_2003 },
278             { Predefined.TOD_CONVENTIONS_2003_SIMPLE_EOP,        Predefined.MOD_CONVENTIONS_2003 },
279             { Predefined.TOD_CONVENTIONS_2010_ACCURATE_EOP,      Predefined.MOD_CONVENTIONS_2010 },
280             { Predefined.TOD_CONVENTIONS_2010_SIMPLE_EOP,        Predefined.MOD_CONVENTIONS_2010 },
281             { Predefined.MOD_WITHOUT_EOP_CORRECTIONS,            Predefined.EME2000 },
282             { Predefined.MOD_CONVENTIONS_1996,                   Predefined.GCRF    },
283             { Predefined.MOD_CONVENTIONS_2003,                   Predefined.EME2000 },
284             { Predefined.MOD_CONVENTIONS_2010,                   Predefined.EME2000 },
285             { Predefined.TEME,                                   Predefined.TOD_WITHOUT_EOP_CORRECTIONS }
286         };
287         for (final Predefined[] pair : reference) {
288             Frame child  = FramesFactory.getFrame(pair[0]);
289             Frame parent = FramesFactory.getFrame(pair[1]);
290             Assertions.assertEquals(parent.getName(), child.getParent().getName(), "wrong parent for " + child.getName());
291         }
292     }
293 
294     @Test
295     public void testEOPConversionSymetry1980() {
296         Utils.setDataRoot("rapid-data-columns");
297         IERSConventions.NutationCorrectionConverter converter =
298                 IERSConventions.IERS_1996.getNutationCorrectionConverter();
299         SortedSet<EOPEntry> rawEquinox = new TreeSet<EOPEntry>(new ChronologicalComparator());
300         DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
301         UTCScale utc = DataContext.getDefault().getTimeScales().getUTC();
302         new RapidDataAndPredictionColumnsLoader(false, "^finals\\.daily$", manager, () -> utc)
303                 .fillHistory(converter, rawEquinox);
304         Assertions.assertEquals(181, rawEquinox.size());
305         for (final EOPEntry entry : rawEquinox) {
306             final double[] rebuiltEquinox = converter.toEquinox(entry.getDate(),
307                                                                 entry.getDx(), entry.getDy());
308             Assertions.assertEquals(entry.getDdPsi(), rebuiltEquinox[0], 2.0e-22);
309             Assertions.assertEquals(entry.getDdEps(), rebuiltEquinox[1], 2.0e-23);
310         }
311     }
312 
313     @Test
314     public void testEOPConversionSymetry2003() {
315         Utils.setDataRoot("rapid-data-columns");
316         IERSConventions.NutationCorrectionConverter converter =
317                 IERSConventions.IERS_2003.getNutationCorrectionConverter();
318         final SortedSet<EOPEntry> rawNRO = new TreeSet<EOPEntry>(new ChronologicalComparator());
319         DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
320         UTCScale utc = DataContext.getDefault().getTimeScales().getUTC();
321         new RapidDataAndPredictionColumnsLoader(true, "^finals2000A\\.daily$", manager, () -> utc)
322                 .fillHistory(converter, rawNRO);
323         Assertions.assertEquals(181, rawNRO.size());
324         for (final EOPEntry entry : rawNRO) {
325             final double[] rebuiltNRO = converter.toNonRotating(entry.getDate(),
326                                                                 entry.getDdPsi(), entry.getDdEps());
327             Assertions.assertEquals(entry.getDx(), rebuiltNRO[0], 6.0e-23);
328             Assertions.assertEquals(entry.getDy(), rebuiltNRO[1], 2.0e-23);
329         }
330     }
331 
332     @Test
333     public void testCIP2003() {
334         testCIP(IERSConventions.IERS_2003, 1.2e-11);
335     }
336 
337     @Test
338     public void testCIP2010() {
339         testCIP(IERSConventions.IERS_2010, 1.2e-11);
340     }
341 
342     private void testCIP(IERSConventions conventions, double threshold) {
343         Utils.setLoaders(conventions, new ArrayList<EOPEntry>());
344         Frame cirf = FramesFactory.getCIRF(conventions, false);
345         Frame tod  = FramesFactory.getTOD(conventions, false);
346         AbsoluteDate t0 = new AbsoluteDate(new DateComponents(2003, 06, 21), TimeComponents.H00,
347                                            TimeScalesFactory.getUTC());
348         for (double dt = -10 * Constants.JULIAN_DAY; dt < 10 * Constants.JULIAN_DAY; dt += 7200) {
349             // CIRF and TOD should both have the Celestial Intermediate Pole as their Z axis
350             AbsoluteDate date = t0.shiftedBy(dt);
351             Transform t = FramesFactory.getNonInterpolatingTransform(tod, cirf, date);
352             Vector3D z = t.transformVector(Vector3D.PLUS_K);
353             Assertions.assertEquals(0.0, Vector3D.angle(z, Vector3D.PLUS_K), threshold);
354         }
355     }
356 
357     @Test
358     public void testEOPConversion() {
359 
360         // real data from buletinb-298.txt
361 
362         // first use case: don't propagate the dx, dy correction to TOD, set dPsi, dEpsilon to 0.0
363         final List<EOPEntry> forced = Utils.buildEOPList(IERSConventions.IERS_2010, ITRFVersion.ITRF_2008, new double[][] {
364             { 56202, 0.3726886, 0.0008843, 0.168556, 0.332869, 0.0, 0.0, -0.000118, 0.000091 },
365             { 56203, 0.3719108, 0.0007204, 0.168261, 0.331527, 0.0, 0.0, -0.000140, 0.000111 },
366             { 56204, 0.3712561, 0.0006217, 0.168218, 0.330668, 0.0, 0.0, -0.000165, 0.000148 },
367             { 56205, 0.3706736, 0.0005530, 0.167775, 0.329688, 0.0, 0.0, -0.000188, 0.000189 },
368             { 56206, 0.3701593, 0.0005139, 0.166829, 0.328457, 0.0, 0.0, -0.000180, 0.000203 }
369         });
370 
371         Utils.setLoaders(IERSConventions.IERS_2010, forced);
372         Frame cirf            = FramesFactory.getCIRF(IERSConventions.IERS_2010, false);
373         Frame todNoCorrection = FramesFactory.getTOD(IERSConventions.IERS_2010, false);
374 
375         // second use case: convert dx, dy data into dDPsi, dDEpsilon
376         final List<EOPEntry> converted = Utils.buildEOPList(IERSConventions.IERS_2010, ITRFVersion.ITRF_2008, new double[][] {
377             { 56202, 0.3726886, 0.0008843, 0.168556, 0.332869, Double.NaN, Double.NaN, -0.000118, 0.000091 },
378             { 56203, 0.3719108, 0.0007204, 0.168261, 0.331527, Double.NaN, Double.NaN, -0.000140, 0.000111 },
379             { 56204, 0.3712561, 0.0006217, 0.168218, 0.330668, Double.NaN, Double.NaN, -0.000165, 0.000148 },
380             { 56205, 0.3706736, 0.0005530, 0.167775, 0.329688, Double.NaN, Double.NaN, -0.000188, 0.000189 },
381             { 56206, 0.3701593, 0.0005139, 0.166829, 0.328457, Double.NaN, Double.NaN, -0.000180, 0.000203 }
382         });
383         Utils.setLoaders(IERSConventions.IERS_2010, converted);
384         Frame todConvertedCorrection  = FramesFactory.getTOD(IERSConventions.IERS_2010, false);
385 
386         for (AbsoluteDate date = forced.get(0).getDate();
387              date.compareTo(forced.get(forced.size() - 1).getDate()) < 0;
388              date = date.shiftedBy(3600)) {
389             Transform tNoCorrection =
390                     FramesFactory.getNonInterpolatingTransform(todNoCorrection, cirf, date);
391 
392             // when we forget the correction on TOD,
393             // its Z axis is slightly offset from CIRF Z axis
394             Vector3D zNoCorrection  = tNoCorrection.transformVector(Vector3D.PLUS_K);
395             Assertions.assertTrue(Vector3D.angle(zNoCorrection, Vector3D.PLUS_K) > 7.2e-10);
396             Transform tConverted =
397                     FramesFactory.getNonInterpolatingTransform(todConvertedCorrection, cirf, date);
398             // when we convert the correction and apply it to TOD,
399             // its Z axis is much better aligned with CIRF Z axis
400             Vector3D zConverted  = tConverted.transformVector(Vector3D.PLUS_K);
401             Assertions.assertTrue(Vector3D.angle(zConverted, Vector3D.PLUS_K) < 6e-12);
402 
403         }
404 
405     }
406 
407     @Test
408     public void testEOPConversionUAI2000Package() {
409         // the reference value has been computed using the uai2000.package routines
410         // provided by Ch. Bizouard on page http://hpiers.obspm.fr/eop-pc/models/models_fr.html
411         // using the following main program
412         //
413         //        program test_EOP_conversion
414         //        double precision dmjd,dpsi,deps,dx1,dy1,dx2,dy2
415         //C       2004-02-14:00:00:00Z, MJD = 53049,
416         //C                             UT1-UTC=-0.4093475, LOD = 0.4676,
417         //C                             X = -0.076804, Y = 0.204671,
418         //C                             dx  = -0.075, dy  = -0.189
419         //C       values extracted from finals2000A.all file, bulletinA columns
420         //        dmjd = 53049.0
421         //        dx1  = -0.075
422         //        dy1  = -0.189
423         //        call DPSIDEPS2000_DXDY2000(dmjd,dx1,dy1,dpsi,deps)
424         //        write (6, *) 'dPsi = ', dpsi, 'dEpsilon = ', deps
425         //        call DXDY2000_DPSIDEPS2000(dmjd,dpsi,deps,dX2,dY2)
426         //        write (6, *) 'dx = ', dx2, 'dy = ', dy2
427         //      end
428         //
429         // the output of this test reads:
430         //     dPsi =  -0.18810999708158463      dEpsilon =  -0.18906891450729962
431         //     dx =   -7.5000002980232239E-002 dy =  -0.18899999558925629
432 
433         IERSConventions.NutationCorrectionConverter converter =
434                 IERSConventions.IERS_2003.getNutationCorrectionConverter();
435         AbsoluteDate date =
436                 new AbsoluteDate(new DateComponents(DateComponents.MODIFIED_JULIAN_EPOCH, 53049),
437                                  TimeScalesFactory.getUTC());
438         double dx  = Constants.ARC_SECONDS_TO_RADIANS * -0.075;
439         double dy  = Constants.ARC_SECONDS_TO_RADIANS * -0.189;
440         double[] equinox = converter.toEquinox(date, dx, dy);
441 
442         // The code from uai2000.package uses sin(epsilon0), cos(epsilon0) in the
443         // formula, whereas IERS conventions use sin(epsilonA), cos(epsilon0), i.e.
444         // the sine is computed on current date instead of epoch. This explains why
445         // the following threshold had to be raised to 2e-11.
446         // We still decided to stick with sin(epsilonA) in our implementation, in
447         // order to remain consistent with IERS conventions
448         Assertions.assertEquals(Constants.ARC_SECONDS_TO_RADIANS * -0.18810999708158463,
449                             equinox[0], 2.0e-11);
450 
451         Assertions.assertEquals(Constants.ARC_SECONDS_TO_RADIANS * -0.18906891450729962,
452                             equinox[1], 2.2e-14);
453         double[] nro = converter.toNonRotating(date, equinox[0], equinox[1]);
454         Assertions.assertEquals(dx, nro[0], 1.0e-20);
455         Assertions.assertEquals(dy, nro[1], 1.0e-20);
456 
457     }
458 
459     @Test
460     public void testFieldConsistency() {
461         for (final Predefined predefined : Predefined.values()) {
462             final Frame frame = FramesFactory.getFrame(predefined);
463             final Frame parent = frame.getParent();
464             if (parent != null) {
465                 double maxPositionError             = 0;
466                 double maxVelocityError             = 0;
467                 double maxAccelerationError         = 0;
468                 double maxRotationError             = 0;
469                 double maxRotationRateError         = 0;
470                 double maxRotationAccelerationError = 0;
471                 for (double dt = 0; dt < Constants.JULIAN_DAY; dt += 60.0) {
472                     final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(dt);
473                     final Transform transformDouble = parent.getTransformTo(frame, date);
474                     final FieldTransform<Binary64> transformD64 =
475                                     parent.getTransformTo(frame, new FieldAbsoluteDate<>(Binary64Field.getInstance(), date));
476                     maxPositionError             = FastMath.max(maxPositionError,
477                                                                 Vector3D.distance(transformDouble.getTranslation(),
478                                                                                   transformD64.getTranslation().toVector3D()));
479                     maxVelocityError             = FastMath.max(maxVelocityError,
480                                                                 Vector3D.distance(transformDouble.getVelocity(),
481                                                                                   transformD64.getVelocity().toVector3D()));
482                     maxAccelerationError         = FastMath.max(maxAccelerationError,
483                                                                 Vector3D.distance(transformDouble.getAcceleration(),
484                                                                                   transformD64.getAcceleration().toVector3D()));
485                     maxRotationError             = FastMath.max(maxRotationError,
486                                                                Rotation.distance(transformDouble.getRotation(),
487                                                                                  transformD64.getRotation().toRotation()));
488                     maxRotationRateError         = FastMath.max(maxRotationRateError,
489                                                                 Vector3D.distance(transformDouble.getRotationRate(),
490                                                                                   transformD64.getRotationRate().toVector3D()));
491                     maxRotationAccelerationError = FastMath.max(maxRotationAccelerationError,
492                                                                Vector3D.distance(transformDouble.getRotationAcceleration(),
493                                                                                  transformD64.getRotationAcceleration().toVector3D()));
494                 }
495                 Assertions.assertEquals(0, maxPositionError,             1.0e-100);
496                 Assertions.assertEquals(0, maxVelocityError,             1.0e-100);
497                 Assertions.assertEquals(0, maxAccelerationError,         1.0e-100);
498                 Assertions.assertEquals(0, maxRotationError,             2.0e-14);
499                 Assertions.assertEquals(0, maxRotationRateError,         2.0e-18);
500                 Assertions.assertEquals(0, maxRotationAccelerationError, 8.0e-22);
501             }
502         }
503     }
504 
505     @Test
506     public void testDerivatives2000WithInterpolation() {
507         doTestDerivatives(AbsoluteDate.J2000_EPOCH, Constants.JULIAN_DAY, 60.0, false,
508                           8.0e-5, 2.0e-5, 3.0e-8, 4.0e-11, 7.0e-13, 2.0e-14);
509     }
510 
511     @Test
512     public void testDerivatives2000WithoutInterpolation() {
513         // when we forbid interpolation, the test is really slow (almost two hours
514         // runtime on a very old machine for one day time span and one minute rate),
515         // we drastically reduce sampling to circumvent this drawback
516         doTestDerivatives(AbsoluteDate.J2000_EPOCH, Constants.JULIAN_DAY / 4, 3600, true,
517                           8.0e-5, 2.0e-5, 3.0e-8, 4.0e-11, 7.0e-13, 2.0e-14);
518     }
519 
520     @Test
521     public void testDerivatives2003WithInterpolation() {
522         doTestDerivatives(AbsoluteDate.J2000_EPOCH.shiftedBy(3 * Constants.JULIAN_YEAR),
523                           Constants.JULIAN_DAY, 60.0, false,
524                           8.0e-5, 2.0e-5, 4.0e-8, 8.0e-12, 3.0e-13, 4.0e-15);
525     }
526 
527     @Test
528     public void testDerivatives2003WithoutInterpolation() {
529         // when we forbid interpolation, the test is really slow (almost two hours
530         // runtime on a very old machine for one day time span and one minute rate),
531         // we drastically reduce sampling to circumvent this drawback
532         doTestDerivatives(AbsoluteDate.J2000_EPOCH.shiftedBy(3 * Constants.JULIAN_YEAR),
533                           Constants.JULIAN_DAY / 4, 3600, true,
534                           8.0e-5, 2.0e-5, 4.0e-8, 8.0e-12, 3.0e-13, 4.0e-15);
535     }
536 
537     @Test
538     public void testNoEOPHistoryFound() {
539         Assertions.assertNull(FramesFactory.findEOP(null));
540         Assertions.assertNull(FramesFactory.findEOP(FramesFactory.getEcliptic(IERSConventions.IERS_2010)));
541         Assertions.assertNull(FramesFactory.findEOP(FramesFactory.getGCRF()));
542         Assertions.assertNull(FramesFactory.findEOP(FramesFactory.getTEME()));
543         Assertions.assertNull(FramesFactory.findEOP(FramesFactory.getICRF()));
544         Assertions.assertNull(FramesFactory.findEOP(CelestialBodyFactory.getEarth().getBodyOrientedFrame()));
545         Assertions.assertNull(FramesFactory.findEOP(FramesFactory.getGTOD(false)));
546     }
547 
548     @Test
549     public void testEOPHistoryFound() {
550         Assertions.assertNotNull(FramesFactory.findEOP(FramesFactory.getCIRF(IERSConventions.IERS_2010, true)));
551         Assertions.assertNotNull(FramesFactory.findEOP(FramesFactory.getCIRF(IERSConventions.IERS_2010, false)));
552         Assertions.assertNotNull(FramesFactory.findEOP(FramesFactory.getITRF(IERSConventions.IERS_2010, true)));
553         Assertions.assertNotNull(FramesFactory.findEOP(FramesFactory.getITRF(IERSConventions.IERS_2010, false)));
554         Assertions.assertNotNull(FramesFactory.findEOP(FramesFactory.getTOD(IERSConventions.IERS_2010, true)));
555         Assertions.assertNotNull(FramesFactory.findEOP(FramesFactory.getTOD(IERSConventions.IERS_2010, false)));
556         Assertions.assertNotNull(FramesFactory.findEOP(FramesFactory.getGTOD(true)));
557         Assertions.assertNotNull(FramesFactory.findEOP(new TopocentricFrame(new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
558                                                                                              Constants.WGS84_EARTH_FLATTENING,
559                                                                                              FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
560                                                                         new GeodeticPoint(FastMath.toRadians(-22.27113),
561                                                                                           FastMath.toRadians(166.44332),
562                                                                                           0.0), "Nouméa")));
563     }
564 
565     @Test
566     public void testCustomEOP() {
567         final double deltaUT1 = 0.1;
568         final Frame      baseITRF    = FramesFactory.getITRF(IERSConventions.IERS_2010, false);
569         final EOPHistory baseEOP     = FramesFactory.findEOP(baseITRF);
570         final EOPHistory modifiedEOP = new EOPHistory(baseEOP.getConventions(),
571                                                       EOPHistory.DEFAULT_INTERPOLATION_DEGREE,
572                                                       baseEOP.
573                                                       getEntries().
574                                                       stream().
575                                                       map(e -> new EOPEntry(e.getMjd(),
576                                                                             e.getUT1MinusUTC() + deltaUT1, e.getLOD(),
577                                                                             e.getX(), e.getY(), e.getXRate(), e.getYRate(),
578                                                                             e.getDdPsi(), e.getDdEps(),
579                                                                             e.getDx(), e.getDy(),
580                                                                             e.getITRFType(), e.getDate())).
581                                                       collect(Collectors.toList()),
582                                                       baseEOP.isSimpleEop());
583         final Frame modifiedITRF = FramesFactory.buildUncachedITRF(modifiedEOP, TimeScalesFactory.getUTC());
584         final AbsoluteDate t0 = new AbsoluteDate(2003, 2, 14, 13, 59, 43.0, TimeScalesFactory.getUTC());
585         for (double dt = 0; dt < 7 * Constants.JULIAN_DAY; dt += 3600) {
586             final Transform t = baseITRF.getTransformTo(modifiedITRF, t0.shiftedBy(dt));
587             Assertions.assertEquals(Constants.WGS84_EARTH_ANGULAR_VELOCITY * deltaUT1,
588                                     t.getRotation().getAngle(),
589                                     1.6e-13);
590         }
591     }
592 
593     private void doTestDerivatives(AbsoluteDate ref,
594                                    double duration, double step, boolean forbidInterpolation,
595                                    double cartesianTolerance, double cartesianDotTolerance, double cartesianDotDotTolerance,
596                                    double rodriguesTolerance, double rodriguesDotTolerance, double rodriguesDotDotTolerance) {
597 
598         final DSFactory factory = new DSFactory(1, 2);
599         final FieldAbsoluteDate<DerivativeStructure> refDS = new FieldAbsoluteDate<>(factory.getDerivativeField(), ref);
600         FiniteDifferencesDifferentiator differentiator = new FiniteDifferencesDifferentiator(8, 60.0);
601 
602         for (final Predefined predefined : Predefined.values()) {
603             final Frame frame = FramesFactory.getFrame(predefined);
604             final Frame parent = frame.getParent();
605             if (parent != null) {
606 
607                 UnivariateDifferentiableVectorFunction dCartesian = differentiator.differentiate(new UnivariateVectorFunction() {
608                     @Override
609                     public double[] value(double t) {
610                         return forbidInterpolation ?
611                                FramesFactory.getNonInterpolatingTransform(parent, frame, ref.shiftedBy(t)).getTranslation().toArray() :
612                                parent.getTransformTo(frame, ref.shiftedBy(t)).getTranslation().toArray();
613                     }
614                 });
615 
616                 UnivariateDifferentiableVectorFunction dOrientation = differentiator.differentiate(new UnivariateVectorFunction() {
617                     double sign = +1.0;
618                     Rotation previous = Rotation.IDENTITY;
619                     @Override
620                     public double[] value(double t) {
621                         AngularCoordinates ac = forbidInterpolation ?
622                                                 FramesFactory.getNonInterpolatingTransform(parent, frame, ref.shiftedBy(t)).getAngular() :
623                                                 parent.getTransformTo(frame, ref.shiftedBy(t)).getAngular();
624                         final double dot = MathArrays.linearCombination(ac.getRotation().getQ0(), previous.getQ0(),
625                                                                         ac.getRotation().getQ1(), previous.getQ1(),
626                                                                         ac.getRotation().getQ2(), previous.getQ2(),
627                                                                         ac.getRotation().getQ3(), previous.getQ3());
628                         sign = FastMath.copySign(1.0, dot * sign);
629                         previous = ac.getRotation();
630                         return ac.getModifiedRodrigues(sign)[0];
631                     }
632                 });
633 
634                 double maxCartesianError       = 0;
635                 double maxCartesianDotError    = 0;
636                 double maxCartesianDotDotError = 0;
637                 double maxRodriguesError       = 0;
638                 double maxRodriguesDotError    = 0;
639                 double maxRodriguesDotDotError = 0;
640                 for (double dt = 0; dt < duration; dt += step) {
641 
642                     final DerivativeStructure dtDS = factory.variable(0, dt);
643                     final FieldTransform<DerivativeStructure> tDS = forbidInterpolation ?
644                                                                     FramesFactory.getNonInterpolatingTransform(parent, frame, refDS.shiftedBy(dtDS)) :
645                                                                     parent.getTransformTo(frame, refDS.shiftedBy(dtDS));
646 
647                     final DerivativeStructure[] refCart   = dCartesian.value(dtDS);
648                     final DerivativeStructure[] fieldCart = tDS.getTranslation().toArray();
649                     for (int i = 0; i < 3; ++i) {
650                         maxCartesianError       = FastMath.max(maxCartesianError,       FastMath.abs(refCart[i].getValue()              - fieldCart[i].getValue()));
651                         maxCartesianDotError    = FastMath.max(maxCartesianDotError,    FastMath.abs(refCart[i].getPartialDerivative(1) - fieldCart[i].getPartialDerivative(1)));
652                         maxCartesianDotDotError = FastMath.max(maxCartesianDotDotError, FastMath.abs(refCart[i].getPartialDerivative(2) - fieldCart[i].getPartialDerivative(2)));
653                     }
654 
655                     final DerivativeStructure[] refOr   = dOrientation.value(dtDS);
656                     DerivativeStructure[] fieldOr = tDS.getAngular().getModifiedRodrigues(1.0)[0];
657                     final double dot = refOr[0].linearCombination(refOr, fieldOr).getReal();
658                     if (dot < 0 || Double.isNaN(dot)) {
659                         fieldOr = tDS.getAngular().getModifiedRodrigues(-1.0)[0];
660                     }
661                     for (int i = 0; i < 3; ++i) {
662                         maxRodriguesError       = FastMath.max(maxRodriguesError,       FastMath.abs(refOr[i].getValue()              - fieldOr[i].getValue()));
663                         maxRodriguesDotError    = FastMath.max(maxRodriguesDotError,    FastMath.abs(refOr[i].getPartialDerivative(1) - fieldOr[i].getPartialDerivative(1)));
664                         maxRodriguesDotDotError = FastMath.max(maxRodriguesDotDotError, FastMath.abs(refOr[i].getPartialDerivative(2) - fieldOr[i].getPartialDerivative(2)));
665                     }
666 
667                 }
668                 Assertions.assertEquals(0, maxCartesianError,             cartesianTolerance,frame.getName());
669                 Assertions.assertEquals(0, maxCartesianDotError,          cartesianDotTolerance,frame.getName());
670                 Assertions.assertEquals(0, maxCartesianDotDotError,       cartesianDotDotTolerance,frame.getName());
671                 Assertions.assertEquals(0, maxRodriguesError,             rodriguesTolerance,frame.getName());
672                 Assertions.assertEquals(0, maxRodriguesDotError,          rodriguesDotTolerance,frame.getName());
673                 Assertions.assertEquals(0, maxRodriguesDotDotError,       rodriguesDotDotTolerance,frame.getName());
674             }
675         }
676     }
677 
678     @BeforeEach
679     public void setUp() {
680         Utils.setDataRoot("regular-data");
681     }
682 
683 }