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