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  
20  import java.util.ArrayList;
21  
22  import org.hipparchus.geometry.euclidean.threed.Rotation;
23  import org.hipparchus.geometry.euclidean.threed.RotationConvention;
24  import org.hipparchus.geometry.euclidean.threed.Vector3D;
25  import org.hipparchus.util.FastMath;
26  import org.hipparchus.util.MathUtils;
27  import org.junit.Assert;
28  import org.junit.Before;
29  import org.junit.Test;
30  import org.orekit.Utils;
31  import org.orekit.time.AbsoluteDate;
32  import org.orekit.time.DateComponents;
33  import org.orekit.time.TTScale;
34  import org.orekit.time.TimeComponents;
35  import org.orekit.time.TimeScale;
36  import org.orekit.time.TimeScalesFactory;
37  import org.orekit.time.UT1Scale;
38  import org.orekit.utils.Constants;
39  import org.orekit.utils.IERSConventions;
40  import org.orekit.utils.PVCoordinates;
41  
42  
43  public class ITRFEquinoxProviderTest {
44  
45      @Test
46      public void testEquinoxVersusCIO() {
47          Frame itrfEquinox  = FramesFactory.getITRFEquinox(IERSConventions.IERS_1996, true);
48          Frame itrfCIO      = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
49          AbsoluteDate start = new AbsoluteDate(2011, 4, 10, TimeScalesFactory.getUTC());
50          AbsoluteDate end   = new AbsoluteDate(2011, 7,  4, TimeScalesFactory.getUTC());
51          for (AbsoluteDate date = start; date.compareTo(end) < 0; date = date.shiftedBy(10000)) {
52              double angularOffset =
53                      itrfEquinox.getTransformTo(itrfCIO, date).getRotation().getAngle();
54              Assert.assertEquals(0, angularOffset / Constants.ARC_SECONDS_TO_RADIANS, 0.07);
55          }
56      }
57  
58      @Test
59      public void testAASReferenceLEO() {
60  
61          // this reference test has been extracted from the following paper:
62          // Implementation Issues Surrounding the New IAU Reference Systems for Astrodynamics
63          // David A. Vallado, John H. Seago, P. Kenneth Seidelmann
64          // http://www.centerforspace.com/downloads/files/pubs/AAS-06-134.pdf
65          Utils.setLoaders(IERSConventions.IERS_1996,
66                           Utils.buildEOPList(IERSConventions.IERS_1996, ITRFVersion.ITRF_2008, new double[][] {
67                               { 53098, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
68                               { 53099, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
69                               { 53100, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
70                               { 53101, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
71                               { 53102, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
72                               { 53103, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
73                               { 53104, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
74                               { 53105, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN }
75                           }));
76          AbsoluteDate t0 = new AbsoluteDate(new DateComponents(2004, 04, 06),
77                                             new TimeComponents(07, 51, 28.386009),
78                                             TimeScalesFactory.getUTC());
79  
80          // ITRF
81          PVCoordinates pvITRF =
82             new PVCoordinates(new Vector3D(-1033479.3830, 7901295.2754, 6380356.5958),
83                               new Vector3D(-3225.636520, -2872.451450, 5531.924446));
84  
85          // GTOD
86          PVCoordinates pvGTOD =
87              new PVCoordinates(new Vector3D(-1033475.0313, 7901305.5856, 6380344.5328),
88                                new Vector3D(-3225.632747, -2872.442511, 5531.931288));
89  
90          Transform t = FramesFactory.getGTOD(IERSConventions.IERS_1996, true).
91                  getTransformTo(FramesFactory.getITRFEquinox(IERSConventions.IERS_1996, true), t0);
92          checkPV(pvITRF, t.transformPVCoordinates(pvGTOD), 8.08e-5, 3.78e-7);
93  
94      }
95  
96      @Test
97      public void testAASReferenceGEO() {
98  
99          // this reference test has been extracted from the following paper:
100         // Implementation Issues Surrounding the New IAU Reference Systems for Astrodynamics
101         // David A. Vallado, John H. Seago, P. Kenneth Seidelmann
102         // http://www.centerforspace.com/downloads/files/pubs/AAS-06-134.pdf
103         Utils.setLoaders(IERSConventions.IERS_1996,
104                          Utils.buildEOPList(IERSConventions.IERS_1996, ITRFVersion.ITRF_2008, new double[][] {
105                              { 53153, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
106                              { 53154, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
107                              { 53155, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
108                              { 53156, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
109                              { 53157, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
110                              { 53158, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
111                              { 53159, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
112                              { 53160, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN }
113                          }));
114         AbsoluteDate t0 = new AbsoluteDate(new DateComponents(2004, 06, 01),
115                                            TimeComponents.H00,
116                                            TimeScalesFactory.getUTC());
117 
118         Transform t = FramesFactory.getGTOD(IERSConventions.IERS_1996, true).
119                 getTransformTo(FramesFactory.getITRFEquinox(IERSConventions.IERS_1996, true), t0);
120 
121         // GTOD
122         PVCoordinates pvGTOD =
123             new PVCoordinates(new Vector3D(24796919.2956, -34115870.9001, 10293.2583),
124                               new Vector3D(-0.979178, -1.476540, -0.928772));
125 
126         // ITRF
127         PVCoordinates pvITRF =
128             new PVCoordinates(new Vector3D(24796919.2915, -34115870.9234, 10226.0621),
129                               new Vector3D(-0.979178, -1.476538, -0.928776));
130 
131         checkPV(pvITRF, t.transformPVCoordinates(pvGTOD), 3.954e-4, 4.69e-7);
132 
133     }
134 
135     @Test
136     public void testSofaCookbook() {
137 
138         // SOFA cookbook test case:
139         //     date       2007 April 05, 12h00m00s.0 UTC
140         //     xp         +0′′.0349282
141         //     yp         +0′′.4833163
142         //     UT1 − UTC  -0s.072073685
143         //     dψ 1980    -0′′.0550655
144         //     dε 1980    -0′′.0063580
145         //     dX 2000    +0′′.0001725
146         //     dY 2000    -0′′.0002650
147         //     dX 2006    +0′′.0001750
148         //     dY 2006    -0′′.0002259
149 
150         Utils.setLoaders(IERSConventions.IERS_1996,
151                          Utils.buildEOPList(IERSConventions.IERS_1996, ITRFVersion.ITRF_2008, new double[][] {
152                              { 54192, -0.072073685,  1.4020, 0.0349282, 0.4833163, -0.0550666, -0.0063580, Double.NaN, Double.NaN },
153                              { 54193, -0.072073685,  1.4020, 0.0349282, 0.4833163, -0.0550666, -0.0063580, Double.NaN, Double.NaN },
154                              { 54194, -0.072073685,  1.4020, 0.0349282, 0.4833163, -0.0550666, -0.0063580, Double.NaN, Double.NaN },
155                              { 54195, -0.072073685,  1.4020, 0.0349282, 0.4833163, -0.0550666, -0.0063580, Double.NaN, Double.NaN },
156                              { 54196, -0.072073685,  1.4020, 0.0349282, 0.4833163, -0.0550666, -0.0063580, Double.NaN, Double.NaN },
157                              { 54197, -0.072073685,  1.4020, 0.0349282, 0.4833163, -0.0550666, -0.0063580, Double.NaN, Double.NaN },
158                              { 54198, -0.072073685,  1.4020, 0.0349282, 0.4833163, -0.0550666, -0.0063580, Double.NaN, Double.NaN },
159                              { 54199, -0.072073685,  1.4020, 0.0349282, 0.4833163, -0.0550666, -0.0063580, Double.NaN, Double.NaN }
160                          }));
161 
162         EOPHistory eopHistory = FramesFactory.getEOPHistory(IERSConventions.IERS_1996, true);
163 
164         TimeScale utc = TimeScalesFactory.getUTC();
165         TTScale tt    = TimeScalesFactory.getTT();
166         UT1Scale ut1  = TimeScalesFactory.getUT1(eopHistory);
167         Frame gcrf    = FramesFactory.getGCRF();
168         Frame itrf    = FramesFactory.getITRFEquinox(IERSConventions.IERS_1996, true);
169         Frame gtod    = itrf.getParent();
170         Frame tod     = gtod.getParent();
171 
172         // time scales checks
173         AbsoluteDate date = new AbsoluteDate(new DateComponents(2007, 4, 5), TimeComponents.H12, utc);
174         Assert.assertEquals(0.50075444444444,
175                             date.getComponents(tt).getTime().getSecondsInUTCDay() / Constants.JULIAN_DAY,
176                             5.0e-15);
177         Assert.assertEquals(0.499999165813831,
178                             date.getComponents(ut1).getTime().getSecondsInUTCDay() / Constants.JULIAN_DAY,
179                             1.0e-15);
180 
181         // sidereal time check
182         double gast = IERSConventions.IERS_1996.getGASTFunction(ut1, eopHistory).value(date);
183         Assert.assertEquals(13.412402380740 * 3600 * 1.0e6,
184                             radToMicroAS(MathUtils.normalizeAngle(gast, 0)),
185                             25);
186 
187         // nutation/precession/bias matrix check
188         Rotation refNPB = new Rotation(new double[][] {
189             { +0.999998403176203, -0.001639032970562, -0.000712190961847 },
190             { +0.001639000942243, +0.999998655799521, -0.000045552846624 },
191             { +0.000712264667137, +0.000044385492226, +0.999999745354454 }
192         }, 1.0e-13);
193         Rotation npb = gcrf.getTransformTo(tod, date).getRotation();
194         Assert.assertEquals(0.0, radToMicroAS(Rotation.distance(refNPB, npb)), 27.0);
195 
196         // celestial to terrestrial frames matrix, without polar motion
197         Rotation refWithoutPolarMotion = new Rotation(new double[][] {
198             { +0.973104317592265, +0.230363826166883, -0.000703332813776 },
199             { -0.230363798723533, +0.973104570754697, +0.000120888299841 },
200             { +0.000712264667137, +0.000044385492226, +0.999999745354454 }
201         }, 1.0e-13);
202         Rotation withoutPM = gcrf.getTransformTo(gtod, date).getRotation();
203         Assert.assertEquals(0.0, radToMicroAS(Rotation.distance(refWithoutPolarMotion, withoutPM)), 9);
204 
205         // celestial to terrestrial frames matrix, with polar motion
206         Rotation refWithPolarMotion = new Rotation(new double[][] {
207             { +0.973104317712772, +0.230363826174782, -0.000703163477127 },
208             { -0.230363800391868, +0.973104570648022, +0.000118545116892 },
209             { +0.000711560100206, +0.000046626645796, +0.999999745754058 }
210         }, 1.0e-13);
211         Rotation withPM = gcrf.getTransformTo(itrf, date).getRotation();
212         Assert.assertEquals(0.0, radToMicroAS(Rotation.distance(refWithPolarMotion, withPM)), 10);
213 
214     }
215 
216     @Test
217     public void testNROvsEquinoxRealEOP() {
218         Utils.setDataRoot("regular-data");
219         checkFrames(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
220                     FramesFactory.getITRFEquinox(IERSConventions.IERS_2010, true),
221                     1.7);
222     }
223 
224     @Test
225     public void testNROvsEquinoxNoEOP2010() {
226         Utils.setLoaders(IERSConventions.IERS_2010, new ArrayList<EOPEntry>());
227         checkFrames(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
228                     FramesFactory.getITRFEquinox(IERSConventions.IERS_2010, true),
229                     1.7);
230     }
231 
232     @Test
233     public void testNROvsEquinoxNoEOP2003() {
234         Utils.setLoaders(IERSConventions.IERS_2003, new ArrayList<EOPEntry>());
235         checkFrames(FramesFactory.getITRF(IERSConventions.IERS_2003, true),
236                     FramesFactory.getITRFEquinox(IERSConventions.IERS_2003, true),
237                     1.9);
238     }
239 
240     @Test
241     public void testNROvsEquinoxNoEOP1996() {
242         Utils.setLoaders(IERSConventions.IERS_1996, new ArrayList<EOPEntry>());
243         checkFrames(FramesFactory.getITRF(IERSConventions.IERS_1996, true),
244                     FramesFactory.getITRFEquinox(IERSConventions.IERS_1996, true),
245                     100);
246     }
247 
248     private void checkFrames(Frame frame1, Frame frame2, double toleranceMicroAS)
249         {
250         AbsoluteDate t0 = new AbsoluteDate(2005, 5, 30, TimeScalesFactory.getUTC());
251         for (double dt = 0; dt < Constants.JULIAN_YEAR; dt += Constants.JULIAN_DAY / 4) {
252             AbsoluteDate date = t0.shiftedBy(dt);
253             Transform t = FramesFactory.getNonInterpolatingTransform(frame1, frame2, date);
254             Vector3D a = t.getRotation().getAxis(RotationConvention.VECTOR_OPERATOR);
255             double delta = FastMath.copySign(radToMicroAS(t.getRotation().getAngle()), a.getZ());
256             Assert.assertEquals(0.0, delta, toleranceMicroAS);
257         }
258     }
259 
260     @Before
261     public void setUp() {
262         Utils.setDataRoot("rapid-data-columns");
263     }
264 
265     private void checkPV(PVCoordinates reference, PVCoordinates result,
266                          double expectedPositionError, double expectedVelocityError) {
267 
268         Vector3D dP = result.getPosition().subtract(reference.getPosition());
269         Vector3D dV = result.getVelocity().subtract(reference.getVelocity());
270         Assert.assertEquals(expectedPositionError, dP.getNorm(), 0.01 * expectedPositionError);
271         Assert.assertEquals(expectedVelocityError, dV.getNorm(), 0.01 * expectedVelocityError);
272     }
273 
274     double radToMicroAS(double deltaRad) {
275         return deltaRad * 1.0e6 / Constants.ARC_SECONDS_TO_RADIANS;
276     }
277 
278 }