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