1   /* Copyright 2011-2012 Space Applications Services
2    * Licensed to CS Communication & Systèmes (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.models.earth.displacement;
18  
19  import org.hipparchus.geometry.euclidean.threed.Vector3D;
20  import org.hipparchus.util.FastMath;
21  import org.junit.jupiter.api.Assertions;
22  import org.junit.jupiter.api.BeforeEach;
23  import org.junit.jupiter.api.Test;
24  import org.orekit.Utils;
25  import org.orekit.data.BodiesElements;
26  import org.orekit.data.FundamentalNutationArguments;
27  import org.orekit.frames.Frame;
28  import org.orekit.frames.FramesFactory;
29  import org.orekit.time.AbsoluteDate;
30  import org.orekit.time.TimeScale;
31  import org.orekit.time.TimeScalesFactory;
32  import org.orekit.utils.Constants;
33  import org.orekit.utils.IERSConventions;
34  import org.orekit.utils.PVCoordinatesProvider;
35  import org.orekit.utils.TimeStampedPVCoordinates;
36  
37  import java.lang.reflect.Field;
38  import java.lang.reflect.InvocationTargetException;
39  import java.lang.reflect.Method;
40  
41  public class TidalDisplacementTest {
42  
43      @Test
44      public void testIERSDisplacementNumbers() {
45          for (final IERSConventions conventions : IERSConventions.values()) {
46              // as of Orekit 9.0, supported conventions are
47              // IERS conventions 1996, IERS conventions 2003 and IERS conventions 2010
48              // and they all share the same values for anelastic Earth model
49              double[] hl = conventions.getNominalTidalDisplacement();
50              Assertions.assertEquals(13, hl.length);
51              Assertions.assertEquals( 0.6078,  hl[ 0], 1.0e-15); // h⁽⁰⁾
52              Assertions.assertEquals(-0.0006,  hl[ 1], 1.0e-15); // h⁽²⁾
53              Assertions.assertEquals( 0.292,   hl[ 2], 1.0e-15); // h₃
54              Assertions.assertEquals(-0.0025,  hl[ 3], 1.0e-15); // hI diurnal
55              Assertions.assertEquals(-0.0022,  hl[ 4], 1.0e-15); // hI semi-diurnal
56              Assertions.assertEquals( 0.0847,  hl[ 5], 1.0e-15); // l⁽⁰⁾
57              Assertions.assertEquals( 0.0012,  hl[ 6], 1.0e-15); // l⁽¹⁾ diurnal
58              Assertions.assertEquals( 0.0024,  hl[ 7], 1.0e-15); // l⁽¹⁾ semi-diurnal
59              Assertions.assertEquals( 0.0002,  hl[ 8], 1.0e-15); // l⁽²⁾
60              Assertions.assertEquals( 0.015,   hl[ 9], 1.0e-15); // l₃
61              Assertions.assertEquals(-0.0007,  hl[10], 1.0e-15); // lI diurnal
62              Assertions.assertEquals(-0.0007,  hl[11], 1.0e-15); // lI semi-diurnal
63              Assertions.assertEquals(-0.31460, hl[12], 1.0e-15); // H₀ permanent deformation amplitude
64          }
65      }
66  
67      @Test
68      public void testDehantOriginalArgumentsNoRemove() {
69          // this test intends to reproduce as much as possible the DEHANTTIDEINEL.F test case
70          // it does so by replacing the fundamental nutation arguments and frequency correction
71          // models used by Orekit (which come from IERS conventions) by the hard-coded models
72          // in the aforementioned program
73          // The results should be really close to the reference values from the original test case
74          doTestDehant(IERSConventions.IERS_2010, false, true,
75                       0.07700420357108125891, 0.06304056321824967613, 0.05516568152597246810,
76                       1.0e-14);
77      }
78  
79      @Test
80      public void testDehantOriginalArgumentsRemovePermanentTide() {
81          // this test intends to reproduce as much as possible the DEHANTTIDEINEL.F test case
82          // with step 3 activated
83          // it does so by replacing the fundamental nutation arguments and frequency correction
84          // models used by Orekit (which come from IERS conventions) by the hard-coded models
85          // in the aforementioned program
86          // The results should be really close to the reference values from the original test case
87          doTestDehant(IERSConventions.IERS_2010, true, true,
88                       0.085888815560994619, 0.065071968475918243, 0.10369393753580475,
89                       1.0e-14);
90      }
91  
92      @Test
93      public void testDehantIERS1996() {
94          // this test intends to replay the DEHANTTIDEINEL.F test case but using the
95          // fundamental nutation arguments and frequency correction models from IERS
96          // conventions.
97          // The results should be very slightly different from the reference values
98          doTestDehant(IERSConventions.IERS_1996, false, false,
99                       0.07700420357108125891, 0.06304056321824967613, 0.05516568152597246810,
100                      5.8e-4);
101     }
102 
103     @Test
104     public void testDehantIERS2003() {
105         // this test intends to replay the DEHANTTIDEINEL.F test case but using the
106         // fundamental nutation arguments and frequency correction models from IERS
107         // conventions.
108         // The results should be very slightly different from the reference values
109         doTestDehant(IERSConventions.IERS_2003, false, false,
110                      0.07700420357108125891, 0.06304056321824967613, 0.05516568152597246810,
111                      2.3e-5);
112     }
113 
114     @Test
115     public void testDehantIERS2010() {
116         // this test intends to replay the DEHANTTIDEINEL.F test case but using the
117         // fundamental nutation arguments and frequency correction models from IERS
118         // conventions.
119         // The results should be very slightly different from the reference values
120         // the differences are more important than with conventions 2003
121         // because after discussion with Dr. Hana Krásná from TU Wien,
122         // we have fixed a typo in the IERS 2010 table 7.3a (∆Rf(op) for
123         // tide P₁ is +0.07mm but the conventions list it as -0.07mm)
124         doTestDehant(IERSConventions.IERS_2010, false, false,
125                      0.07700420357108125891, 0.06304056321824967613, 0.05516568152597246810,
126                      1.3e-4);
127     }
128 
129     private void doTestDehant(final IERSConventions conventions, final boolean removePermanentDeformation, final boolean replaceModels,
130                               final double expectedDx, final double expectedDy, final double expectedDz,
131                               final double tolerance)
132         {
133 
134         Frame     itrf = FramesFactory.getITRF(conventions, false);
135         TimeScale ut1  = TimeScalesFactory.getUT1(conventions, false);
136 
137         final double re;
138         final double sunEarthSystemMassRatio ;
139         final double earthMoonMassRatio;
140         if (replaceModels) {
141             // constants consistent with DEHANTTIDEINEL.F reference program
142             // available at <ftp://tai.bipm.org/iers/conv2010/chapter7/dehanttideinel/>
143             // and Copyright (C) 2008 IERS Conventions Center
144             re                         = 6378136.6;
145             final double massRatioSun  = 332946.0482;
146             final double massRatioMoon = 0.0123000371;
147             sunEarthSystemMassRatio    = massRatioSun * (1.0 / (1.0 + massRatioMoon));
148             earthMoonMassRatio         = 1.0 / massRatioMoon;
149         } else {
150             // constants consistent with IERS and JPL
151             re                      = Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS;
152             sunEarthSystemMassRatio = Constants.JPL_SSD_SUN_EARTH_PLUS_MOON_MASS_RATIO;
153             earthMoonMassRatio      = Constants.JPL_SSD_EARTH_MOON_MASS_RATIO;
154         }
155 
156         // fake providers generating only the positions from the reference program test
157         PVCoordinatesProvider fakeSun  = (date, frame) -> new TimeStampedPVCoordinates(date,
158                                                                                        new Vector3D(137859926952.015,
159                                                                                                     54228127881.435,
160                                                                                                     23509422341.6960),
161                                                                                        Vector3D.ZERO,
162                                                                                        Vector3D.ZERO);
163         PVCoordinatesProvider fakeMoon = (date, frame) -> new TimeStampedPVCoordinates(date,
164                                                                                        new Vector3D(-179996231.920342,
165                                                                                                     -312468450.131567,
166                                                                                                     -169288918.592160),
167                                                                                        Vector3D.ZERO,
168                                                                                        Vector3D.ZERO);
169 
170         TidalDisplacement td = new TidalDisplacement(re, sunEarthSystemMassRatio, earthMoonMassRatio,
171                                                      fakeSun, fakeMoon, conventions, removePermanentDeformation);
172 
173         FundamentalNutationArguments arguments = null;
174         if (replaceModels) {
175             try {
176                 // we override the official IERS conventions 2010 arguments with fake arguments matching DEHANTTIDEINEL.F code
177                 String regularArguments = "/assets/org/orekit/IERS-conventions/2010/nutation-arguments.txt";
178                 arguments = new FundamentalNutationArguments(conventions, ut1,
179                                                              IERSConventions.class.getResourceAsStream(regularArguments),
180                                                              regularArguments) {
181 
182                     @Override
183                     public BodiesElements evaluateAll(final AbsoluteDate date) {
184                         BodiesElements base = super.evaluateAll(date);
185                         double fhr = date.getComponents(ut1).getTime().getSecondsInUTCDay() / 3600.0;
186                         double t = base.getTC();
187 
188                         // Doodson fundamental arguments as per DEHANTTIDEINEL.F code
189                         double s   = 218.31664563 + (481267.88194 + (-0.0014663889 + (0.00000185139) * t) * t) * t;
190                         double tau = fhr * 15 + 280.4606184 + (36000.7700536 + (0.00038793 + (-0.0000000258) * t) * t) * t - s;
191                         double pr  = (1.396971278 + (0.000308889 + (0.000000021 + (0.000000007) * t) * t) * t) * t ;
192                         double h   = 280.46645 + (36000.7697489 + (0.00030322222  + (0.000000020 + (-0.00000000654) * t) * t) * t) * t;
193                         double p   = 83.35324312 + (4069.01363525 + (-0.01032172222 + (-0.0000124991 + (0.00000005263) * t) * t) * t) * t;
194                         double zns = 234.95544499 + (1934.13626197 + (-0.00207561111 + (-0.00000213944 + (0.00000001650) * t) * t) * t) * t;
195                         double ps  = 282.93734098 + (1.71945766667 + (0.00045688889 + (-0.00000001778 + (-0.00000000334) * t) * t) * t) * t;
196                         s += pr;
197 
198                         // rebuild Delaunay arguments from Doodson arguments, ignoring derivatives
199                         return new BodiesElements(date, base.getTC(),
200                                                   FastMath.toRadians(s + tau), 0.0, FastMath.toRadians(s - p), 0.0, FastMath.toRadians(h - ps), 0.0,
201                                                   FastMath.toRadians(s + zns), 0.0, FastMath.toRadians(s - h), 0.0, FastMath.toRadians(-zns),   0.0,
202                                                   base.getLMe(),               0.0, base.getLVe(),             0.0, base.getLE(),               0.0,
203                                                   base.getLMa(),               0.0, base.getLJu(),             0.0, base.getLSa(),              0.0,
204                                                   base.getLUr(),               0.0, base.getLNe(),             0.0, base.getPa(),               0.0);
205 
206                     }
207                 };
208 
209                 // we override the official IERS conventions 2010 tides displacements with tides displacements matching DEHANTTIDEINEL.F code
210                 String table73a = "/tides/tab7.3a-Dehant.txt";
211                 Field diurnalCorrectionField = td.getClass().getDeclaredField("frequencyCorrectionDiurnal");
212                 diurnalCorrectionField.setAccessible(true);
213                 Method diurnalCorrectionGetter = IERSConventions.class.
214                                                  getDeclaredMethod("getTidalDisplacementFrequencyCorrectionDiurnal",
215                                                                    String.class, Integer.TYPE, Integer.TYPE, Integer.TYPE, Integer.TYPE, Integer.TYPE);
216                 diurnalCorrectionGetter.setAccessible(true);
217                 diurnalCorrectionField.set(td, diurnalCorrectionGetter.invoke(null, table73a, 18, 15, 16, 17, 18));
218 
219             } catch (SecurityException | NoSuchMethodException | NoSuchFieldException |
220                      InvocationTargetException | IllegalArgumentException | IllegalAccessException e) {
221                 Assertions.fail(e.getLocalizedMessage());
222             }
223         } else {
224             arguments = conventions.getNutationArguments(ut1);
225         }
226 
227 
228         Vector3D fundamentalStationWettzell = new Vector3D(4075578.385, 931852.890, 4801570.154);
229         AbsoluteDate date = new AbsoluteDate(2009, 4, 13, 0, 0, 0.0, ut1);
230         Vector3D displacement = td.displacement(arguments.evaluateAll(date), itrf, fundamentalStationWettzell);
231         Assertions.assertEquals(expectedDx, displacement.getX(), tolerance);
232         Assertions.assertEquals(expectedDy, displacement.getY(), tolerance);
233         Assertions.assertEquals(expectedDz, displacement.getZ(), tolerance);
234 
235 
236     }
237 
238     @BeforeEach
239     public void setUp() throws Exception {
240         Utils.setDataRoot("regular-data");
241     }
242 
243 }