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 org.hipparchus.geometry.euclidean.threed.Rotation;
20  import org.hipparchus.geometry.euclidean.threed.RotationConvention;
21  import org.hipparchus.geometry.euclidean.threed.Vector3D;
22  import org.hipparchus.util.FastMath;
23  import org.junit.jupiter.api.Assertions;
24  import org.junit.jupiter.api.BeforeEach;
25  import org.junit.jupiter.api.Test;
26  import org.orekit.Utils;
27  import org.orekit.data.DataContext;
28  import org.orekit.time.AbsoluteDate;
29  import org.orekit.time.DateComponents;
30  import org.orekit.time.TimeComponents;
31  import org.orekit.time.TimeScalesFactory;
32  import org.orekit.utils.AngularDerivativesFilter;
33  import org.orekit.utils.CartesianDerivativesFilter;
34  import org.orekit.utils.Constants;
35  import org.orekit.utils.IERSConventions;
36  import org.orekit.utils.OrekitConfiguration;
37  import org.orekit.utils.PVCoordinates;
38  
39  import java.io.FileNotFoundException;
40  
41  
42  public class TODProviderTest {
43  
44      @Test
45      public void testRotationRate() {
46          TransformProvider provider =
47                  new InterpolatingTransformProvider(new TODProvider(IERSConventions.IERS_1996, null, DataContext.getDefault().getTimeScales()),
48                                                     CartesianDerivativesFilter.USE_PVA,
49                                                     AngularDerivativesFilter.USE_R,
50                                                     3, 1.0, 5, Constants.JULIAN_DAY, 100.0);
51          AbsoluteDate tMin = new AbsoluteDate(2035, 3, 2, 15, 58, 59, TimeScalesFactory.getUTC());
52          double minRate = provider.getTransform(tMin).getRotationRate().getNorm();
53          Assertions.assertEquals(6.4e-14, minRate, 1.0e-15);
54          AbsoluteDate tMax = new AbsoluteDate(2043, 12, 16, 14, 18, 9, TimeScalesFactory.getUTC());
55          double maxRate = provider.getTransform(tMax).getRotationRate().getNorm();
56          Assertions.assertEquals(1.4e-11, maxRate, 1.0e-12);
57      }
58  
59      @Test
60      public void testAASReferenceLEO() {
61  
62          // this reference test has been extracted from the following paper:
63          // Implementation Issues Surrounding the New IAU Reference Systems for Astrodynamics
64          // David A. Vallado, John H. Seago, P. Kenneth Seidelmann
65          // http://www.centerforspace.com/downloads/files/pubs/AAS-06-134.pdf
66          Utils.setLoaders(IERSConventions.IERS_1996,
67                           Utils.buildEOPList(IERSConventions.IERS_1996, ITRFVersion.ITRF_2008, new double[][] {
68                               { 53098, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
69                               { 53099, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
70                               { 53100, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
71                               { 53101, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
72                               { 53102, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
73                               { 53103, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
74                               { 53104, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
75                               { 53105, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN }
76                           }));
77          AbsoluteDate t0 = new AbsoluteDate(new DateComponents(2004, 04, 06),
78                                             new TimeComponents(07, 51, 28.386009),
79                                             TimeScalesFactory.getUTC());
80  
81          Transform tt = FramesFactory.getMOD(IERSConventions.IERS_1996).
82                  getTransformTo(FramesFactory.getTOD(IERSConventions.IERS_1996, true), t0);
83          Transform ff = FramesFactory.getMOD(false).getTransformTo(FramesFactory.getTOD(false), t0);
84  
85          //TOD iau76
86          PVCoordinates pvTODiau76 =
87              new PVCoordinates(new Vector3D(5094514.7804, 6127366.4612, 6380344.5328),
88                                new Vector3D(-4746.088567, 786.077222, 5531.931288));
89          //MOD iau76
90          PVCoordinates pvMODiau76WithoutNutCorr =
91              new PVCoordinates(new Vector3D(5094029.0167, 6127870.9363, 6380247.8885),
92                                new Vector3D(-4746.262495, 786.014149, 5531.791025));
93          //MOD iau76
94          PVCoordinates pvMODiau76 =
95              new PVCoordinates(new Vector3D(5094028.3745, 6127870.8164, 6380248.5164),
96                                new Vector3D(-4746.263052, 786.014045, 5531.790562));
97  
98          // it seems the induced effect of pole nutation correction δΔψ on the equation of the equinoxes
99          // was not taken into account in the reference paper, so we fix it here for the test
100         final double dDeltaPsi =
101                 FramesFactory.getEOPHistory(IERSConventions.IERS_1996, true).getEquinoxNutationCorrection(t0)[0];
102         final double epsilonA = IERSConventions.IERS_1996.getMeanObliquityFunction().value(t0);
103         final Transform fix =
104                 new Transform(t0, new Rotation(Vector3D.PLUS_K,
105                                                dDeltaPsi * FastMath.cos(epsilonA),
106                                                RotationConvention.FRAME_TRANSFORM));
107 
108         checkPV(pvTODiau76, fix.transformPVCoordinates(tt.transformPVCoordinates(pvMODiau76)), 1.13e-3, 5.3e-5);
109         checkPV(pvTODiau76, ff.transformPVCoordinates(pvMODiau76WithoutNutCorr), 1.07e-3, 5.3e-5);
110 
111     }
112 
113     @Test
114     public void testAASReferenceGEO() {
115 
116         // this reference test has been extracted from the following paper:
117         // Implementation Issues Surrounding the New IAU Reference Systems for Astrodynamics
118         // David A. Vallado, John H. Seago, P. Kenneth Seidelmann
119         // http://www.centerforspace.com/downloads/files/pubs/AAS-06-134.pdf
120         Utils.setLoaders(IERSConventions.IERS_1996,
121                          Utils.buildEOPList(IERSConventions.IERS_1996, ITRFVersion.ITRF_2008, new double[][] {
122                              { 53153, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
123                              { 53154, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
124                              { 53155, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
125                              { 53156, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
126                              { 53157, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
127                              { 53158, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
128                              { 53159, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
129                              { 53160, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN }
130                          }));
131         AbsoluteDate t0 = new AbsoluteDate(new DateComponents(2004, 06, 01),
132                                            TimeComponents.H00,
133                                            TimeScalesFactory.getUTC());
134 
135         Transform tt = FramesFactory.getMOD(IERSConventions.IERS_1996).
136                 getTransformTo(FramesFactory.getTOD(IERSConventions.IERS_1996, true), t0);
137         Transform ff = FramesFactory.getMOD(false).getTransformTo(FramesFactory.getTOD(false), t0);
138 
139         // TOD iau76
140         PVCoordinates pvTODiau76 =
141             new PVCoordinates(new Vector3D(-40577427.7501, -11500096.1306, 10293.2583),
142                               new Vector3D(837.552338, -2957.524176, -0.928772));
143         // MOD iau76
144         PVCoordinates pvMODiau76WithoutNutCorr =
145             new PVCoordinates(new Vector3D(-40576822.6385, -11502231.5013, 9738.2304),
146                               new Vector3D(837.708020, -2957.480118, -0.814275));
147 
148         // MOD iau76
149         PVCoordinates pvMODiau76 =
150             new PVCoordinates(new Vector3D(-40576822.6395, -11502231.5015, 9733.7842),
151                               new Vector3D(837.708020, -2957.480117, -0.814253));
152 
153 
154         // it seems the induced effect of pole nutation correction δΔψ on the equation of the equinoxes
155         // was not taken into account in the reference paper, so we fix it here for the test
156         final double dDeltaPsi =
157                 FramesFactory.getEOPHistory(IERSConventions.IERS_1996, true).getEquinoxNutationCorrection(t0)[0];
158         final double epsilonA = IERSConventions.IERS_1996.getMeanObliquityFunction().value(t0);
159         final Transform fix =
160                 new Transform(t0, new Rotation(Vector3D.PLUS_K,
161                                                dDeltaPsi * FastMath.cos(epsilonA),
162                                                RotationConvention.FRAME_TRANSFORM));
163 
164         checkPV(pvTODiau76, fix.transformPVCoordinates(tt.transformPVCoordinates(pvMODiau76)), 4.86e-4, 6.2e-5);
165         checkPV(pvTODiau76, ff.transformPVCoordinates(pvMODiau76WithoutNutCorr), 4.87e-4, 6.31e-5);
166 
167     }
168 
169     @Test
170     public void testInterpolationAccuracyWithEOP() throws FileNotFoundException {
171 
172         // max interpolation error observed on a one month period with 60 seconds step
173         //
174         // number of sample points    time between sample points    max error
175         //        6                          86400s /  8 =  3h       19.56e-12 rad
176         //        6                          86400s / 12 =  2h       13.02e-12 rad
177         //        6                          86400s / 16 =  1h30      9.75e-12 rad
178         //        6                          86400s / 20 =  1h12      7.79e-12 rad
179         //        6                          86400s / 24 =  1h        6.48e-12 rad
180         //        8                          86400s /  8 =  3h       20.91e-12 rad
181         //        8                          86400s / 12 =  2h       13.91e-12 rad
182         //        8                          86400s / 16 =  1h30     10.42e-12 rad
183         //        8                          86400s / 20 =  1h12      8.32e-12 rad
184         //        8                          86400s / 24 =  1h        6.92e-12 rad
185         //       10                          86400s /  8 =  3h       21.65e-12 rad
186         //       10                          86400s / 12 =  2h       14.41e-12 rad
187         //       10                          86400s / 16 =  1h30     10.78e-12 rad
188         //       10                          86400s / 20 =  1h12      8.61e-12 rad
189         //       10                          86400s / 24 =  1h        7.16e-12 rad
190         //       12                          86400s /  8 =  3h       22.12e-12 rad
191         //       12                          86400s / 12 =  2h       14.72e-12 rad
192         //       12                          86400s / 16 =  1h30     11.02e-12 rad
193         //       12                          86400s / 20 =  1h12      8.80e-12 rad
194         //       12                          86400s / 24 =  1h        7.32e-12 rad
195         //
196         // looking at error behavior during along the sample show the max error is
197         // a peak at 00h00 each day for all curves, which matches the EOP samples
198         // points used for correction (eopHistoru is set to non null at construction here).
199         // So looking only at max error does not allow to select an interpolation
200         // setting as they all fall in a similar 6e-12 to 8e-12 range. Looking at
201         // the error behavior between these peaks however shows that there is still
202         // some signal if the time interval is between sample points is too large,
203         // in order to get only numerical noise, we have to go as far as 1h between
204         // the points.
205         // We finally select 6 interpolation points separated by 1 hour each
206         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_1996, false);
207         TransformProvider nonInterpolating = new TODProvider(IERSConventions.IERS_1996, eopHistory, DataContext.getDefault().getTimeScales());
208         final TransformProvider interpolating =
209                 new InterpolatingTransformProvider(nonInterpolating,
210                                                    CartesianDerivativesFilter.USE_PVA,
211                                                    AngularDerivativesFilter.USE_R,
212                                                    6, Constants.JULIAN_DAY / 24,
213                                                    OrekitConfiguration.getCacheSlotsNumber(),
214                                                    Constants.JULIAN_YEAR, 30 * Constants.JULIAN_DAY);
215 
216         // the following time range is located around the maximal observed error
217         AbsoluteDate start = new AbsoluteDate(2002, 11, 11, 0, 0, 0.0, TimeScalesFactory.getTAI());
218         AbsoluteDate end   = new AbsoluteDate(2002, 11, 15, 6, 0, 0.0, TimeScalesFactory.getTAI());
219         double maxError = 0.0;
220         for (AbsoluteDate date = start; date.compareTo(end) < 0; date = date.shiftedBy(60)) {
221             final Transform transform =
222                     new Transform(date,
223                                   interpolating.getTransform(date),
224                                   nonInterpolating.getTransform(date).getInverse());
225             final double error = transform.getRotation().getAngle();
226             maxError = FastMath.max(maxError, error);
227         }
228 
229         Assertions.assertTrue(maxError < 7e-12);
230 
231     }
232 
233     @Test
234     public void testInterpolationAccuracyWithoutEOP() throws FileNotFoundException {
235 
236         // max interpolation error observed on a one month period with 60 seconds step
237         //
238         // number of sample points    time between sample points    max error
239         //        5                          86400s /  3 =  8h     3286.90e-15 rad
240         //        5                          86400s /  6 =  4h      103.90e-15 rad
241         //        5                          86400s /  8 =  3h       24.74e-15 rad
242         //        5                          86400s / 12 =  2h        4.00e-15 rad
243         //        6                          86400s /  3 =  8h      328.91e-15 rad
244         //        6                          86400s /  6 =  4h        5.92e-15 rad
245         //        6                          86400s /  8 =  3h        3.95e-15 rad
246         //        6                          86400s / 12 =  2h        3.94e-15 rad
247         //        8                          86400s /  3 =  8h        5.87e-15 rad
248         //        8                          86400s /  6 =  4h        4.73e-15 rad
249         //        8                          86400s /  8 =  3h        4.45e-15 rad
250         //        8                          86400s / 12 =  2h        3.87e-15 rad
251         //       10                          86400s /  3 =  8h        5.29e-15 rad
252         //       10                          86400s /  6 =  4h        5.36e-15 rad
253         //       10                          86400s /  8 =  3h        5.86e-15 rad
254         //       10                          86400s / 12 =  2h        5.76e-15 rad
255         //
256         //
257         // We don't see anymore the peak at 00h00 so this confirms it is related to EOP
258         // sampling. All values between 3e-15 and 6e-15 are really equivalent: it is
259         // mostly numerical noise. The best settings are 6 or 8 points every 2 or 3 hours.
260         // We finally select 6 interpolation points separated by 3 hours each
261         TransformProvider nonInterpolating = new TODProvider(IERSConventions.IERS_1996, null, DataContext.getDefault().getTimeScales());
262                 final TransformProvider interpolating =
263                         new InterpolatingTransformProvider(nonInterpolating,
264                                                            CartesianDerivativesFilter.USE_PVA,
265                                                            AngularDerivativesFilter.USE_R,
266                                                            6, Constants.JULIAN_DAY / 8,
267                                                            OrekitConfiguration.getCacheSlotsNumber(),
268                                                            Constants.JULIAN_YEAR, 30 * Constants.JULIAN_DAY);
269 
270                 // the following time range is located around the maximal observed error
271                 AbsoluteDate start = new AbsoluteDate(2002, 11, 11, 0, 0, 0.0, TimeScalesFactory.getTAI());
272                 AbsoluteDate end   = new AbsoluteDate(2002, 11, 15, 6, 0, 0.0, TimeScalesFactory.getTAI());
273                 double maxError = 0.0;
274                 for (AbsoluteDate date = start; date.compareTo(end) < 0; date = date.shiftedBy(60)) {
275                     final Transform transform =
276                             new Transform(date,
277                                           interpolating.getTransform(date),
278                                           nonInterpolating.getTransform(date).getInverse());
279                     final double error = transform.getRotation().getAngle();
280                     maxError = FastMath.max(maxError, error);
281                 }
282 
283                 Assertions.assertTrue(maxError < 4.0e-15);
284 
285     }
286 
287     @Test
288     public void testSofaPnm80() {
289 
290         // the reference value has been computed using the March 2012 version of the SOFA library
291         // http://www.iausofa.org/2012_0301_C.html, with the following code
292         //
293         //        double utc1, utc2, tai1, tai2, tt1, tt2, rmatpn[3][3];
294         //
295         //        // 2004-02-14:00:00:00Z, MJD = 53049, UT1-UTC = -0.4093509
296         //        utc1  = DJM0 + 53049.0;
297         //        utc2  = 0.0;
298         //        iauUtctai(utc1, utc2, &tai1, &tai2);
299         //        iauTaitt(tai1, tai2, &tt1, &tt2);
300         //
301         //        iauPnm80(tt1, tt2, rmatpn);
302         //
303         //        printf("iauPnm80(%.20g, %.20g, rmatpn)\n"
304         //               "  --> %.20g %.20g %.20g\n"
305         //               "      %.20g %.20g %.20g\n"
306         //               "      %.20g %.20g %.20g\n",
307         //               tt1, tt2,
308         //               rmatpn[0][0], rmatpn[0][1], rmatpn[0][2],
309         //               rmatpn[1][0], rmatpn[1][1], rmatpn[1][2],
310         //               rmatpn[2][0], rmatpn[2][1], rmatpn[2][2]);
311         //
312         // the output of this test reads:
313         //        iauNutm80(2453049.5, 0.00074287037037037029902, nut)
314         //         --> 0.99999999859236310407 4.8681019508684473249e-05 2.1105264333587349032e-05
315         //            -4.8680343021901595118e-05 0.99999999830143670998 -3.205231683600651138e-05
316         //            -2.1106824637199909505e-05 3.2051289379386727063e-05 0.99999999926360838565
317         //        iauPnm80(2453049.5, 0.00074287037037037029902, rmatpn)
318         //         --> 0.99999954755358466674 -0.00087243169070689370777 -0.00037915111913272635073
319         //            0.0008724195377896877112 0.99999961892302935418 -3.2217171614061089913e-05
320         //            0.00037917908192846747854 3.1886378193416632805e-05 0.99999992760323874741
321 
322         // As the iauNutm80 and iauPnm80 do not allow user to specify EOP corrections,
323         // the test is done with Predefined.TOD_WITHOUT_EOP_CORRECTIONS.
324 
325         AbsoluteDate date = new AbsoluteDate(2004, 2, 14, TimeScalesFactory.getUTC());
326         Frame tod  = FramesFactory.getFrame(Predefined.TOD_WITHOUT_EOP_CORRECTIONS);
327         checkRotation(new double[][] {
328             { 0.99999999859236310407, 4.8681019508684473249e-05, 2.1105264333587349032e-05 },
329             { -4.8680343021901595118e-05, 0.99999999830143670998, -3.205231683600651138e-05 },
330             { -2.1106824637199909505e-05, 3.2051289379386727063e-05, 0.99999999926360838565    }
331 
332         }, tod.getParent().getTransformTo(tod, date), 5.0e-11);
333         checkRotation(new double[][] {
334             { 0.99999954755358466674,   -0.00087243169070689370777, -0.00037915111913272635073 },
335             { 0.0008724195377896877112,  0.99999961892302935418,    -3.2217171614061089913e-05 },
336             { 0.00037917908192846747854, 3.1886378193416632805e-05,  0.99999992760323874741    }
337 
338         }, tod.getParent().getParent().getTransformTo(tod, date), 5.0e-11);
339 
340     }
341 
342     @Test
343     public void testTOD1976vs2006() {
344 
345         final Frame tod1976 = FramesFactory.getTOD(IERSConventions.IERS_1996, true);
346         final Frame tod2006 = FramesFactory.getTOD(IERSConventions.IERS_2010, true);
347         for (double dt = 0; dt < 2 * Constants.JULIAN_YEAR; dt += 100 * Constants.JULIAN_DAY) {
348             AbsoluteDate date = new AbsoluteDate(AbsoluteDate.J2000_EPOCH, dt);
349             double delta = tod1976.getTransformTo(tod2006, date).getRotation().getAngle();
350             // TOD2006 and TOD2000 are similar to about 65 milli-arcseconds
351             // between 2000 and 2002, with EOP corrections taken into account in both cases
352             Assertions.assertEquals(0.0, delta, 3.2e-7);
353         }
354 
355     }
356 
357     @Test
358     public void testTOD2000vs2006() {
359 
360         final Frame tod2000 = FramesFactory.getTOD(IERSConventions.IERS_2003, true);
361         final Frame tod2006 = FramesFactory.getTOD(IERSConventions.IERS_2010, true);
362         for (double dt = 0; dt < 2 * Constants.JULIAN_YEAR; dt += 100 * Constants.JULIAN_DAY) {
363             AbsoluteDate date = new AbsoluteDate(AbsoluteDate.J2000_EPOCH, dt);
364             double delta = tod2000.getTransformTo(tod2006, date).getRotation().getAngle();
365             // TOD2006 and TOD2000 are similar to about 30 micro-arcseconds
366             // between 2000 and 2002, with EOP corrections taken into account in both cases
367             Assertions.assertEquals(0.0, delta, 1.5e-10);
368         }
369 
370     }
371 
372     @BeforeEach
373     public void setUp() {
374         Utils.setDataRoot("compressed-data");
375     }
376 
377     private void checkPV(PVCoordinates reference, PVCoordinates result,
378                          double expectedPositionError, double expectedVelocityError) {
379 
380         Vector3D dP = result.getPosition().subtract(reference.getPosition());
381         Vector3D dV = result.getVelocity().subtract(reference.getVelocity());
382         Assertions.assertEquals(expectedPositionError, dP.getNorm(), 0.01 * expectedPositionError);
383         Assertions.assertEquals(expectedVelocityError, dV.getNorm(), 0.01 * expectedVelocityError);
384     }
385 
386     private void checkRotation(double[][] reference, Transform t, double epsilon) {
387         double[][] mat = t.getRotation().getMatrix();
388         for (int i = 0; i < 3; ++i) {
389             for (int j = 0; j < 3; ++j) {
390                 Assertions.assertEquals(reference[i][j], mat[i][j], epsilon);
391             }
392         }
393     }
394 
395 }