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.data;
18  
19  
20  import java.io.BufferedReader;
21  import java.io.ByteArrayInputStream;
22  import java.io.ByteArrayOutputStream;
23  import java.io.IOException;
24  import java.io.InputStream;
25  import java.io.InputStreamReader;
26  import java.io.ObjectInputStream;
27  import java.io.ObjectOutputStream;
28  import java.lang.reflect.InvocationTargetException;
29  import java.lang.reflect.Method;
30  import java.nio.charset.StandardCharsets;
31  import java.util.Arrays;
32  import java.util.function.Function;
33  
34  import org.hipparchus.CalculusFieldElement;
35  import org.hipparchus.analysis.UnivariateFunction;
36  import org.hipparchus.analysis.differentiation.DSFactory;
37  import org.hipparchus.analysis.differentiation.DerivativeStructure;
38  import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
39  import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
40  import org.hipparchus.util.Decimal64;
41  import org.hipparchus.util.Decimal64Field;
42  import org.hipparchus.util.FastMath;
43  import org.hipparchus.util.MathUtils;
44  import org.junit.Assert;
45  import org.junit.Before;
46  import org.junit.Test;
47  import org.orekit.Utils;
48  import org.orekit.errors.OrekitException;
49  import org.orekit.errors.OrekitMessages;
50  import org.orekit.time.AbsoluteDate;
51  import org.orekit.time.FieldAbsoluteDate;
52  import org.orekit.time.TimeScale;
53  import org.orekit.time.TimeScalesFactory;
54  import org.orekit.utils.Constants;
55  import org.orekit.utils.IERSConventions;
56  
57  public class FundamentalNutationArgumentsTest {
58  
59      @Test
60      public void testNoStream() {
61          try {
62              new FundamentalNutationArguments(IERSConventions.IERS_2010, TimeScalesFactory.getTT(), null, "dummy");
63              Assert.fail("an exception should have been thrown");
64          } catch (OrekitException oe) {
65              Assert.assertEquals(OrekitMessages.UNABLE_TO_FIND_FILE, oe.getSpecifier());
66              Assert.assertEquals("dummy", oe.getParts()[0]);
67          }
68      }
69  
70      @Test
71      public void testModifiedData() throws IOException {
72  
73          String directory = "/assets/org/orekit/IERS-conventions/";
74          InputStream is = getClass().getResourceAsStream(directory + "2010/nutation-arguments.txt");
75          BufferedReader reader = new BufferedReader(new InputStreamReader(is, StandardCharsets.UTF_8));
76          StringBuilder builder = new StringBuilder();
77          for (String line = reader.readLine(); line != null; line = reader.readLine()) {
78              builder.append(line);
79              builder.append('\n');
80          }
81          reader.close();
82          try {
83              // remove the line:
84              // F5 ≡ Ω = 125.04455501◦ − 6962890.5431″t + 7.4722″t² + 0.007702″t³ − 0.00005939″t⁴
85              String modified = builder.toString().replaceAll("F5[^\\n]+", "");
86              new FundamentalNutationArguments(IERSConventions.IERS_2010, null,
87                                               new ByteArrayInputStream(modified.getBytes()),
88                                               "modified-data");
89              Assert.fail("an exception should have been thrown");
90          } catch (OrekitException oe) {
91              Assert.assertEquals(OrekitMessages.NOT_A_SUPPORTED_IERS_DATA_FILE, oe.getSpecifier());
92          }
93      }
94  
95      @Test
96      public void testEnum() throws NoSuchMethodException, SecurityException,
97                                    IllegalAccessException, IllegalArgumentException,
98                                    InvocationTargetException {
99          Class<?> e = null;
100         for (final Class<?> c : FundamentalNutationArguments.class.getDeclaredClasses()) {
101             if (c.getName().endsWith("FundamentalName")) {
102                 e = c;
103             }
104         }
105         Method m = e.getDeclaredMethod("valueOf", String.class);
106         m.setAccessible(true);
107         for (String n : Arrays.asList("L", "L_PRIME", "F", "D", "OMEGA",
108                                       "L_ME", "L_VE", "L_E", "L_MA", "L_J", "L_SA", "L_U", "L_NE", "PA")) {
109             Assert.assertEquals(n, m.invoke(null, n).toString());
110         }
111         try {
112             m.invoke(null, "inexistent");
113             Assert.fail("an exception should have been thrown");
114         } catch (InvocationTargetException ite) {
115             Assert.assertTrue(ite.getCause() instanceof IllegalArgumentException);
116         }
117     }
118 
119     @Test
120     public void testDotDouble() {
121         final IERSConventions conventions = IERSConventions.IERS_2010;
122         final TimeScale ut1 = TimeScalesFactory.getUT1(conventions, false);
123         final FundamentalNutationArguments fna = conventions.getNutationArguments(ut1);
124         final AbsoluteDate t0 = new AbsoluteDate(2002, 4, 7, 12, 34, 22.5, TimeScalesFactory.getUTC());
125         final UnivariateDifferentiableFunction gamma  = differentiate(fna, t0, b -> b.getGamma());
126         final UnivariateDifferentiableFunction l      = differentiate(fna, t0, b -> b.getL());
127         final UnivariateDifferentiableFunction lPrime = differentiate(fna, t0, b -> b.getLPrime());
128         final UnivariateDifferentiableFunction f      = differentiate(fna, t0, b -> b.getF());
129         final UnivariateDifferentiableFunction d      = differentiate(fna, t0, b -> b.getD());
130         final UnivariateDifferentiableFunction lMe    = differentiate(fna, t0, b -> b.getLMe());
131         final UnivariateDifferentiableFunction lVe    = differentiate(fna, t0, b -> b.getLVe());
132         final UnivariateDifferentiableFunction lE     = differentiate(fna, t0, b -> b.getLE());
133         final UnivariateDifferentiableFunction lMa    = differentiate(fna, t0, b -> b.getLMa());
134         final UnivariateDifferentiableFunction lJu    = differentiate(fna, t0, b -> b.getLJu());
135         final UnivariateDifferentiableFunction lSa    = differentiate(fna, t0, b -> b.getLSa());
136         final UnivariateDifferentiableFunction lUr    = differentiate(fna, t0, b -> b.getLUr());
137         final UnivariateDifferentiableFunction lNe    = differentiate(fna, t0, b -> b.getLNe());
138         final UnivariateDifferentiableFunction  pa    = differentiate(fna, t0, b -> b.getPa());
139         final DSFactory factory = new DSFactory(1, 1);
140         double maxErrorGamma  = 0;
141         double maxErrorL      = 0;
142         double maxErrorLPrime = 0;
143         double maxErrorF      = 0;
144         double maxErrorD      = 0;
145         double maxErrorLMe    = 0;
146         double maxErrorLVe    = 0;
147         double maxErrorLE     = 0;
148         double maxErrorLMa    = 0;
149         double maxErrorLJu    = 0;
150         double maxErrorLSa    = 0;
151         double maxErrorLUr    = 0;
152         double maxErrorLNe    = 0;
153         double maxErrorPa     = 0;
154         for (double dt = 0; dt < Constants.JULIAN_DAY; dt += 60.0) {
155             BodiesElements be = fna.evaluateAll(t0.shiftedBy(dt));
156             DerivativeStructure dtDS = factory.variable(0, dt);
157             maxErrorGamma  = FastMath.max(maxErrorGamma,  FastMath.abs(gamma .value(dtDS).getPartialDerivative(1) - be.getGammaDot()));
158             maxErrorL      = FastMath.max(maxErrorL,      FastMath.abs(l     .value(dtDS).getPartialDerivative(1) - be.getLDot()));
159             maxErrorLPrime = FastMath.max(maxErrorLPrime, FastMath.abs(lPrime.value(dtDS).getPartialDerivative(1) - be.getLPrimeDot()));
160             maxErrorF      = FastMath.max(maxErrorF,      FastMath.abs(f     .value(dtDS).getPartialDerivative(1) - be.getFDot()));
161             maxErrorD      = FastMath.max(maxErrorD,      FastMath.abs(d     .value(dtDS).getPartialDerivative(1) - be.getDDot()));
162             maxErrorLMe    = FastMath.max(maxErrorLMe,    FastMath.abs(lMe   .value(dtDS).getPartialDerivative(1) - be.getLMeDot()));
163             maxErrorLVe    = FastMath.max(maxErrorLVe,    FastMath.abs(lVe   .value(dtDS).getPartialDerivative(1) - be.getLVeDot()));
164             maxErrorLE     = FastMath.max(maxErrorLE,     FastMath.abs(lE    .value(dtDS).getPartialDerivative(1) - be.getLEDot()));
165             maxErrorLMa    = FastMath.max(maxErrorLMa,    FastMath.abs(lMa   .value(dtDS).getPartialDerivative(1) - be.getLMaDot()));
166             maxErrorLJu    = FastMath.max(maxErrorLJu,    FastMath.abs(lJu   .value(dtDS).getPartialDerivative(1) - be.getLJuDot()));
167             maxErrorLSa    = FastMath.max(maxErrorLSa,    FastMath.abs(lSa   .value(dtDS).getPartialDerivative(1) - be.getLSaDot()));
168             maxErrorLUr    = FastMath.max(maxErrorLUr,    FastMath.abs(lUr   .value(dtDS).getPartialDerivative(1) - be.getLUrDot()));
169             maxErrorLNe    = FastMath.max(maxErrorLNe,    FastMath.abs(lNe   .value(dtDS).getPartialDerivative(1) - be.getLNeDot()));
170             maxErrorPa     = FastMath.max(maxErrorPa,     FastMath.abs(pa    .value(dtDS).getPartialDerivative(1) - be.getPaDot()));
171         }
172         Assert.assertEquals(0, maxErrorGamma,  8.0e-13);
173         Assert.assertEquals(0, maxErrorL,      1.0e-14);
174         Assert.assertEquals(0, maxErrorLPrime, 6.0e-16);
175         Assert.assertEquals(0, maxErrorF,      6.0e-15);
176         Assert.assertEquals(0, maxErrorD,      6.0e-15);
177         Assert.assertEquals(0, maxErrorLMe,    2.0e-15);
178         Assert.assertEquals(0, maxErrorLVe,    5.0e-16);
179         Assert.assertEquals(0, maxErrorLE,     3.0e-16);
180         Assert.assertEquals(0, maxErrorLMa,    4.0e-16);
181         Assert.assertEquals(0, maxErrorLJu,    3.0e-17);
182         Assert.assertEquals(0, maxErrorLSa,    4.0e-17);
183         Assert.assertEquals(0, maxErrorLUr,    1.0e-16);
184         Assert.assertEquals(0, maxErrorLNe,    8.0e-17);
185         Assert.assertEquals(0, maxErrorPa,     3.0e-20);
186     }
187 
188     private UnivariateDifferentiableFunction differentiate(final FundamentalNutationArguments fna, final AbsoluteDate t0,
189                                                            final Function<BodiesElements, Double> f) {
190         return new FiniteDifferencesDifferentiator(8, 10.0).differentiate(new UnivariateFunction() {
191             double angle = 0;
192             @Override
193             public double value(double t) {
194                 double raw = f.apply(fna.evaluateAll(t0.shiftedBy(t)));
195                 angle = MathUtils.normalizeAngle(raw, angle);
196                 return angle;
197             }
198         });
199     }
200 
201     @Test
202     public void testDotField() {
203         final IERSConventions conventions = IERSConventions.IERS_2010;
204         final TimeScale ut1 = TimeScalesFactory.getUT1(conventions, false);
205         final FundamentalNutationArguments fna = conventions.getNutationArguments(ut1);
206         final FieldAbsoluteDate<Decimal64> t0 = new FieldAbsoluteDate<>(Decimal64Field.getInstance(),
207                                                                         2002, 4, 7, 12, 34, 22.5, TimeScalesFactory.getUTC());
208         final UnivariateDifferentiableFunction gamma  = differentiate(fna, t0, b -> b.getGamma());
209         final UnivariateDifferentiableFunction l      = differentiate(fna, t0, b -> b.getL());
210         final UnivariateDifferentiableFunction lPrime = differentiate(fna, t0, b -> b.getLPrime());
211         final UnivariateDifferentiableFunction f      = differentiate(fna, t0, b -> b.getF());
212         final UnivariateDifferentiableFunction d      = differentiate(fna, t0, b -> b.getD());
213         final UnivariateDifferentiableFunction lMe    = differentiate(fna, t0, b -> b.getLMe());
214         final UnivariateDifferentiableFunction lVe    = differentiate(fna, t0, b -> b.getLVe());
215         final UnivariateDifferentiableFunction lE     = differentiate(fna, t0, b -> b.getLE());
216         final UnivariateDifferentiableFunction lMa    = differentiate(fna, t0, b -> b.getLMa());
217         final UnivariateDifferentiableFunction lJu    = differentiate(fna, t0, b -> b.getLJu());
218         final UnivariateDifferentiableFunction lSa    = differentiate(fna, t0, b -> b.getLSa());
219         final UnivariateDifferentiableFunction lUr    = differentiate(fna, t0, b -> b.getLUr());
220         final UnivariateDifferentiableFunction lNe    = differentiate(fna, t0, b -> b.getLNe());
221         final UnivariateDifferentiableFunction  pa    = differentiate(fna, t0, b -> b.getPa());
222         final DSFactory factory = new DSFactory(1, 1);
223         double maxErrorGamma  = 0;
224         double maxErrorL      = 0;
225         double maxErrorLPrime = 0;
226         double maxErrorF      = 0;
227         double maxErrorD      = 0;
228         double maxErrorLMe    = 0;
229         double maxErrorLVe    = 0;
230         double maxErrorLE     = 0;
231         double maxErrorLMa    = 0;
232         double maxErrorLJu    = 0;
233         double maxErrorLSa    = 0;
234         double maxErrorLUr    = 0;
235         double maxErrorLNe    = 0;
236         double maxErrorPa     = 0;
237         for (double dt = 0; dt < Constants.JULIAN_DAY; dt += 60.0) {
238             FieldBodiesElements<Decimal64> be = fna.evaluateAll(t0.shiftedBy(dt));
239             DerivativeStructure dtDS = factory.variable(0, dt);
240             maxErrorGamma  = FastMath.max(maxErrorGamma,  FastMath.abs(gamma .value(dtDS).getPartialDerivative(1) - be.getGammaDot().getReal()));
241             maxErrorL      = FastMath.max(maxErrorL,      FastMath.abs(l     .value(dtDS).getPartialDerivative(1) - be.getLDot().getReal()));
242             maxErrorLPrime = FastMath.max(maxErrorLPrime, FastMath.abs(lPrime.value(dtDS).getPartialDerivative(1) - be.getLPrimeDot().getReal()));
243             maxErrorF      = FastMath.max(maxErrorF,      FastMath.abs(f     .value(dtDS).getPartialDerivative(1) - be.getFDot().getReal()));
244             maxErrorD      = FastMath.max(maxErrorD,      FastMath.abs(d     .value(dtDS).getPartialDerivative(1) - be.getDDot().getReal()));
245             maxErrorLMe    = FastMath.max(maxErrorLMe,    FastMath.abs(lMe   .value(dtDS).getPartialDerivative(1) - be.getLMeDot().getReal()));
246             maxErrorLVe    = FastMath.max(maxErrorLVe,    FastMath.abs(lVe   .value(dtDS).getPartialDerivative(1) - be.getLVeDot().getReal()));
247             maxErrorLE     = FastMath.max(maxErrorLE,     FastMath.abs(lE    .value(dtDS).getPartialDerivative(1) - be.getLEDot().getReal()));
248             maxErrorLMa    = FastMath.max(maxErrorLMa,    FastMath.abs(lMa   .value(dtDS).getPartialDerivative(1) - be.getLMaDot().getReal()));
249             maxErrorLJu    = FastMath.max(maxErrorLJu,    FastMath.abs(lJu   .value(dtDS).getPartialDerivative(1) - be.getLJuDot().getReal()));
250             maxErrorLSa    = FastMath.max(maxErrorLSa,    FastMath.abs(lSa   .value(dtDS).getPartialDerivative(1) - be.getLSaDot().getReal()));
251             maxErrorLUr    = FastMath.max(maxErrorLUr,    FastMath.abs(lUr   .value(dtDS).getPartialDerivative(1) - be.getLUrDot().getReal()));
252             maxErrorLNe    = FastMath.max(maxErrorLNe,    FastMath.abs(lNe   .value(dtDS).getPartialDerivative(1) - be.getLNeDot().getReal()));
253             maxErrorPa     = FastMath.max(maxErrorPa,     FastMath.abs(pa    .value(dtDS).getPartialDerivative(1) - be.getPaDot().getReal()));
254         }
255         Assert.assertEquals(0, maxErrorGamma,  8.0e-13);
256         Assert.assertEquals(0, maxErrorL,      1.0e-14);
257         Assert.assertEquals(0, maxErrorLPrime, 6.0e-16);
258         Assert.assertEquals(0, maxErrorF,      6.0e-15);
259         Assert.assertEquals(0, maxErrorD,      6.0e-15);
260         Assert.assertEquals(0, maxErrorLMe,    2.0e-15);
261         Assert.assertEquals(0, maxErrorLVe,    5.0e-16);
262         Assert.assertEquals(0, maxErrorLE,     3.0e-16);
263         Assert.assertEquals(0, maxErrorLMa,    4.0e-16);
264         Assert.assertEquals(0, maxErrorLJu,    3.0e-17);
265         Assert.assertEquals(0, maxErrorLSa,    4.0e-17);
266         Assert.assertEquals(0, maxErrorLUr,    1.0e-16);
267         Assert.assertEquals(0, maxErrorLNe,    8.0e-17);
268         Assert.assertEquals(0, maxErrorPa,     3.0e-20);
269     }
270 
271     private <T extends CalculusFieldElement<T>> UnivariateDifferentiableFunction differentiate(final FundamentalNutationArguments fna, final FieldAbsoluteDate<T> t0,
272                                                                                            final Function<FieldBodiesElements<T>, T> f) {
273         return new FiniteDifferencesDifferentiator(8, 10.0).differentiate(new UnivariateFunction() {
274             double angle = 0;
275             @Override
276             public double value(double t) {
277                 double raw = f.apply(fna.evaluateAll(t0.shiftedBy(t))).getReal();
278                 angle = MathUtils.normalizeAngle(raw, angle);
279                 return angle;
280             }
281         });
282     }
283 
284     @Test
285     public void testSerializationNoTidalCorrection() throws IOException, ClassNotFoundException {
286         IERSConventions conventions = IERSConventions.IERS_2010;
287         TimeScale ut1 = TimeScalesFactory.getUT1(conventions, true);
288         checkSerialization(295000, 300000, conventions.getNutationArguments(ut1));
289     }
290 
291     @Test
292     public void testSerializationTidalCorrection() throws IOException, ClassNotFoundException {
293         IERSConventions conventions = IERSConventions.IERS_2010;
294         TimeScale ut1 = TimeScalesFactory.getUT1(conventions, false);
295         checkSerialization(295000, 300000, conventions.getNutationArguments(ut1));
296     }
297 
298     @Test
299     public void testSerializationNoUT1Correction() throws IOException, ClassNotFoundException {
300         IERSConventions conventions = IERSConventions.IERS_2010;
301         checkSerialization(850, 950, conventions.getNutationArguments(null));
302     }
303 
304     private void checkSerialization(int low, int high, FundamentalNutationArguments nutation)
305         throws IOException, ClassNotFoundException {
306 
307         ByteArrayOutputStream bos = new ByteArrayOutputStream();
308         ObjectOutputStream    oos = new ObjectOutputStream(bos);
309         oos.writeObject(nutation);
310 
311         Assert.assertTrue(bos.size() > low);
312         Assert.assertTrue(bos.size() < high);
313 
314         ByteArrayInputStream  bis = new ByteArrayInputStream(bos.toByteArray());
315         ObjectInputStream     ois = new ObjectInputStream(bis);
316         FundamentalNutationArguments deserialized  = (FundamentalNutationArguments) ois.readObject();
317         for (double dt = 0; dt < Constants.JULIAN_DAY; dt += 3600) {
318             AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(dt);
319             BodiesElements be1 = nutation.evaluateAll(date);
320             BodiesElements be2 = deserialized.evaluateAll(date);
321             Assert.assertEquals(be1.getGamma(),  be2.getGamma(),  1.0e-15);
322             Assert.assertEquals(be1.getL(),      be2.getL(),      1.0e-15);
323             Assert.assertEquals(be1.getLPrime(), be2.getLPrime(), 1.0e-15);
324             Assert.assertEquals(be1.getF(),      be2.getF(),      1.0e-15);
325             Assert.assertEquals(be1.getD(),      be2.getD(),      1.0e-15);
326             Assert.assertEquals(be1.getOmega(),  be2.getOmega(),  1.0e-15);
327             Assert.assertEquals(be1.getLMe(),    be2.getLMe(),    1.0e-15);
328             Assert.assertEquals(be1.getLVe(),    be2.getLVe(),    1.0e-15);
329             Assert.assertEquals(be1.getLE(),     be2.getLE(),     1.0e-15);
330             Assert.assertEquals(be1.getLMa(),    be2.getLMa(),    1.0e-15);
331             Assert.assertEquals(be1.getLJu(),    be2.getLJu(),    1.0e-15);
332             Assert.assertEquals(be1.getLSa(),    be2.getLSa(),    1.0e-15);
333             Assert.assertEquals(be1.getLUr(),    be2.getLUr(),    1.0e-15);
334             Assert.assertEquals(be1.getLNe(),    be2.getLNe(),    1.0e-15);
335             Assert.assertEquals(be1.getPa(),     be2.getPa(),     1.0e-15);
336         }
337 
338     }
339 
340     @Before
341     public void setUp() {
342         Utils.setDataRoot("compressed-data");
343     }
344 
345 }