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.io.ByteArrayInputStream;
21  import java.io.ByteArrayOutputStream;
22  import java.io.IOException;
23  import java.io.ObjectInputStream;
24  import java.io.ObjectOutputStream;
25  
26  import org.hipparchus.CalculusFieldElement;
27  import org.hipparchus.geometry.euclidean.threed.FieldRotation;
28  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
29  import org.hipparchus.geometry.euclidean.threed.Rotation;
30  import org.hipparchus.geometry.euclidean.threed.RotationConvention;
31  import org.hipparchus.geometry.euclidean.threed.Vector3D;
32  import org.hipparchus.util.Decimal64;
33  import org.hipparchus.util.Decimal64Field;
34  import org.hipparchus.util.FastMath;
35  import org.junit.Assert;
36  import org.junit.Before;
37  import org.junit.Test;
38  import org.orekit.Utils;
39  import org.orekit.data.DataContext;
40  import org.orekit.time.AbsoluteDate;
41  import org.orekit.time.DateComponents;
42  import org.orekit.time.FieldAbsoluteDate;
43  import org.orekit.time.TimeComponents;
44  import org.orekit.time.TimeScalesFactory;
45  import org.orekit.utils.FieldPVCoordinates;
46  import org.orekit.utils.IERSConventions;
47  import org.orekit.utils.PVCoordinates;
48  
49  
50  public class GTODProviderTest {
51  
52      @Test
53      public void testAASReferenceLEO() {
54  
55          // this reference test has been extracted from the following paper:
56          // Implementation Issues Surrounding the New IAU Reference Systems for Astrodynamics
57          // David A. Vallado, John H. Seago, P. Kenneth Seidelmann
58          // http://www.centerforspace.com/downloads/files/pubs/AAS-06-134.pdf
59          Utils.setLoaders(IERSConventions.IERS_1996,
60                           Utils.buildEOPList(IERSConventions.IERS_1996, ITRFVersion.ITRF_2008, new double[][] {
61                               { 53098, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
62                               { 53099, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
63                               { 53100, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
64                               { 53101, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
65                               { 53102, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
66                               { 53103, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
67                               { 53104, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN },
68                               { 53105, -0.4399619, 0.0015563, -0.140682, 0.333309, -0.052195, -0.003875, Double.NaN, Double.NaN }
69                           }));
70          AbsoluteDate t0 = new AbsoluteDate(new DateComponents(2004, 04, 06),
71                                             new TimeComponents(07, 51, 28.386009),
72                                             TimeScalesFactory.getUTC());
73  
74          // PEF iau76
75          PVCoordinates pvPEF =
76             new PVCoordinates(new Vector3D(-1033475.0313, 7901305.5856, 6380344.5328),
77                               new Vector3D(-3225.632747, -2872.442511, 5531.931288));
78  
79          // it seems the induced effect of pole nutation correction δΔψ on the equation of the equinoxes
80          // was not taken into account in the reference paper, so we fix it here for the test
81          final double dDeltaPsi =
82                  FramesFactory.getEOPHistory(IERSConventions.IERS_1996, true).getEquinoxNutationCorrection(t0)[0];
83          final double epsilonA = IERSConventions.IERS_1996.getMeanObliquityFunction().value(t0);
84          final Transform fix =
85                  new Transform(t0, new Rotation(Vector3D.PLUS_K,
86                                                 dDeltaPsi * FastMath.cos(epsilonA),
87                                                 RotationConvention.FRAME_TRANSFORM));
88  
89          // TOD iau76
90          PVCoordinates pvTOD =
91              new PVCoordinates(new Vector3D(5094514.7804, 6127366.4612, 6380344.5328),
92                                new Vector3D(-4746.088567, 786.077222, 5531.931288));
93  
94          Transform t = FramesFactory.getTOD(IERSConventions.IERS_1996, true).
95                  getTransformTo(FramesFactory.getGTOD(IERSConventions.IERS_1996, true), t0);
96          checkPV(fix.transformPVCoordinates(pvPEF), t.transformPVCoordinates(pvTOD), 0.00942, 3.12e-5);
97  
98          // if we forget to apply nutation corrections, results are much worse, which is expected
99          t = FramesFactory.getTOD(false).getTransformTo(FramesFactory.getGTOD(false), t0);
100         checkPV(fix.transformPVCoordinates(pvPEF), t.transformPVCoordinates(pvTOD), 257.49, 0.13955);
101 
102     }
103 
104     @Test
105     public void testAASReferenceGEO() {
106 
107         // this reference test has been extracted from the following paper:
108         // Implementation Issues Surrounding the New IAU Reference Systems for Astrodynamics
109         // David A. Vallado, John H. Seago, P. Kenneth Seidelmann
110         // http://www.centerforspace.com/downloads/files/pubs/AAS-06-134.pdf
111         Utils.setLoaders(IERSConventions.IERS_1996,
112                          Utils.buildEOPList(IERSConventions.IERS_1996, ITRFVersion.ITRF_2008, new double[][] {
113                              { 53153, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
114                              { 53154, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
115                              { 53155, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
116                              { 53156, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
117                              { 53157, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
118                              { 53158, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
119                              { 53159, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
120                              { 53160, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN }
121                          }));
122         AbsoluteDate t0 = new AbsoluteDate(new DateComponents(2004, 06, 01),
123                                            TimeComponents.H00,
124                                            TimeScalesFactory.getUTC());
125 
126         Transform t = FramesFactory.getTOD(IERSConventions.IERS_1996, true).
127                 getTransformTo(FramesFactory.getGTOD(IERSConventions.IERS_1996, true), t0);
128         // TOD iau76
129         PVCoordinates pvTOD =
130             new PVCoordinates(new Vector3D(-40577427.7501, -11500096.1306, 10293.2583),
131                               new Vector3D(837.552338, -2957.524176, -0.928772));
132 
133         // PEF iau76
134         PVCoordinates pvPEF =
135             new PVCoordinates(new Vector3D(24796919.2956, -34115870.9001, 10293.2583),
136                               new Vector3D(-0.979178, -1.476540, -0.928772));
137 
138         // it seems the induced effect of pole nutation correction δΔψ on the equation of the equinoxes
139         // was not taken into account in the reference paper, so we fix it here for the test
140         final double dDeltaPsi =
141                 FramesFactory.getEOPHistory(IERSConventions.IERS_1996, true).getEquinoxNutationCorrection(t0)[0];
142         final double epsilonA = IERSConventions.IERS_1996.getMeanObliquityFunction().value(t0);
143         final Transform fix =
144                 new Transform(t0, new Rotation(Vector3D.PLUS_K,
145                                                dDeltaPsi * FastMath.cos(epsilonA),
146                                                RotationConvention.FRAME_TRANSFORM));
147 
148         checkPV(fix.transformPVCoordinates(pvPEF), t.transformPVCoordinates(pvTOD), 0.0503, 3.59e-4);
149 
150         // if we forget to apply nutation corrections, results are much worse, which is expected
151         t = FramesFactory.getTOD(false).getTransformTo(FramesFactory.getGTOD(false), t0);
152         checkPV(fix.transformPVCoordinates(pvPEF), t.transformPVCoordinates(pvTOD), 1458.27, 3.847e-4);
153 
154     }
155 
156     @Test
157     public void testAASReferenceGEOField() {
158 
159         // this reference test has been extracted from the following paper:
160         // Implementation Issues Surrounding the New IAU Reference Systems for Astrodynamics
161         // David A. Vallado, John H. Seago, P. Kenneth Seidelmann
162         // http://www.centerforspace.com/downloads/files/pubs/AAS-06-134.pdf
163         Utils.setLoaders(IERSConventions.IERS_1996,
164                          Utils.buildEOPList(IERSConventions.IERS_1996, ITRFVersion.ITRF_2008, new double[][] {
165                              { 53153, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
166                              { 53154, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
167                              { 53155, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
168                              { 53156, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
169                              { 53157, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
170                              { 53158, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
171                              { 53159, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN },
172                              { 53160, -0.4709050,  0.0000000, -0.083853,  0.467217, -0.053614, -0.004494, Double.NaN, Double.NaN }
173                          }));
174         FieldAbsoluteDate<Decimal64> t0 = new FieldAbsoluteDate<>(Decimal64Field.getInstance(),
175                                                                   new DateComponents(2004, 06, 01),
176                                                                   TimeComponents.H00,
177                                                                   TimeScalesFactory.getUTC());
178 
179         FieldTransform<Decimal64> t = FramesFactory.getTOD(IERSConventions.IERS_1996, true).
180                 getTransformTo(FramesFactory.getGTOD(IERSConventions.IERS_1996, true), t0);
181         // TOD iau76
182         PVCoordinates pvTOD =
183             new PVCoordinates(new Vector3D(-40577427.7501, -11500096.1306, 10293.2583),
184                               new Vector3D(837.552338, -2957.524176, -0.928772));
185 
186         // PEF iau76
187         PVCoordinates pvPEF =
188             new PVCoordinates(new Vector3D(24796919.2956, -34115870.9001, 10293.2583),
189                               new Vector3D(-0.979178, -1.476540, -0.928772));
190 
191         // it seems the induced effect of pole nutation correction δΔψ on the equation of the equinoxes
192         // was not taken into account in the reference paper, so we fix it here for the test
193         final Decimal64 dDeltaPsi =
194                 FramesFactory.getEOPHistory(IERSConventions.IERS_1996, true).getEquinoxNutationCorrection(t0)[0];
195         final Decimal64 epsilonA = IERSConventions.IERS_1996.getMeanObliquityFunction().value(t0);
196         final FieldTransform<Decimal64> fix =
197                 new FieldTransform<>(t0, new FieldRotation<>(FieldVector3D.getPlusK(Decimal64Field.getInstance()),
198                                                              dDeltaPsi.multiply(epsilonA.cos()),
199                                                              RotationConvention.FRAME_TRANSFORM));
200 
201         checkPV(fix.transformPVCoordinates(pvPEF), t.transformPVCoordinates(pvTOD), 0.0503, 3.59e-4);
202 
203         // if we forget to apply nutation corrections, results are much worse, which is expected
204         t = FramesFactory.getTOD(false).getTransformTo(FramesFactory.getGTOD(false), t0);
205         checkPV(fix.transformPVCoordinates(pvPEF), t.transformPVCoordinates(pvTOD), 1458.27, 3.847e-4);
206 
207     }
208 
209     @Test
210     public void testSerialization() throws IOException, ClassNotFoundException {
211         GTODProvider provider = new GTODProvider(IERSConventions.IERS_2010,
212                                                  FramesFactory.getEOPHistory(IERSConventions.IERS_2010, true),
213                                                  DataContext.getDefault().getTimeScales());
214 
215         ByteArrayOutputStream bos = new ByteArrayOutputStream();
216         ObjectOutputStream    oos = new ObjectOutputStream(bos);
217         oos.writeObject(provider);
218 
219         Assert.assertTrue(bos.size() > 295000);
220         Assert.assertTrue(bos.size() < 300000);
221 
222         ByteArrayInputStream  bis = new ByteArrayInputStream(bos.toByteArray());
223         ObjectInputStream     ois = new ObjectInputStream(bis);
224         GTODProvider deserialized  = (GTODProvider) ois.readObject();
225         for (int i = 0; i < FastMath.min(100, provider.getEOPHistory().getEntries().size()); ++i) {
226             AbsoluteDate date = provider.getEOPHistory().getEntries().get(i).getDate();
227             Transform expectedIdentity = new Transform(date,
228                                                        provider.getTransform(date).getInverse(),
229                                                        deserialized.getTransform(date));
230             Assert.assertEquals(0.0, expectedIdentity.getTranslation().getNorm(), 1.0e-15);
231             Assert.assertEquals(0.0, expectedIdentity.getRotation().getAngle(),   1.0e-15);
232         }
233 
234     }
235 
236     @Before
237     public void setUp() {
238         Utils.setDataRoot("compressed-data");
239     }
240 
241     private void checkPV(PVCoordinates reference, PVCoordinates result,
242                          double expectedPositionError, double expectedVelocityError) {
243 
244         Vector3D dP = result.getPosition().subtract(reference.getPosition());
245         Vector3D dV = result.getVelocity().subtract(reference.getVelocity());
246         Assert.assertEquals(expectedPositionError, dP.getNorm(), 0.01 * expectedPositionError);
247         Assert.assertEquals(expectedVelocityError, dV.getNorm(), 0.01 * expectedVelocityError);
248     }
249 
250     private <T extends CalculusFieldElement<T>> void checkPV(FieldPVCoordinates<T> reference,
251                                                          FieldPVCoordinates<T> result,
252                                                          double expectedPositionError,
253                                                          double expectedVelocityError) {
254 
255         FieldVector3D<T> dP = result.getPosition().subtract(reference.getPosition());
256         FieldVector3D<T> dV = result.getVelocity().subtract(reference.getVelocity());
257         Assert.assertEquals(expectedPositionError, dP.getNorm().getReal(), 0.01 * expectedPositionError);
258         Assert.assertEquals(expectedVelocityError, dV.getNorm().getReal(), 0.01 * expectedVelocityError);
259     }
260 
261 }