1   /* Copyright 2002-2017 CS Systèmes d'Information
2    * Licensed to CS Systèmes d'Information (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.forces.gravity;
18  
19  import java.util.Map;
20  
21  import org.hipparchus.geometry.euclidean.threed.Vector3D;
22  import org.hipparchus.ode.AbstractIntegrator;
23  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
24  import org.hipparchus.util.FastMath;
25  import org.junit.Assert;
26  import org.junit.Before;
27  import org.junit.Test;
28  import org.orekit.Utils;
29  import org.orekit.data.DataProvidersManager;
30  import org.orekit.errors.OrekitException;
31  import org.orekit.errors.OrekitMessages;
32  import org.orekit.forces.ForceModel;
33  import org.orekit.forces.gravity.potential.AstronomicalAmplitudeReader;
34  import org.orekit.forces.gravity.potential.FESCHatEpsilonReader;
35  import org.orekit.forces.gravity.potential.GravityFieldFactory;
36  import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
37  import org.orekit.forces.gravity.potential.OceanLoadDeformationCoefficients;
38  import org.orekit.frames.Frame;
39  import org.orekit.frames.FramesFactory;
40  import org.orekit.orbits.KeplerianOrbit;
41  import org.orekit.orbits.Orbit;
42  import org.orekit.orbits.OrbitType;
43  import org.orekit.orbits.PositionAngle;
44  import org.orekit.propagation.SpacecraftState;
45  import org.orekit.propagation.numerical.NumericalPropagator;
46  import org.orekit.time.AbsoluteDate;
47  import org.orekit.time.TimeScale;
48  import org.orekit.time.TimeScalesFactory;
49  import org.orekit.time.UT1Scale;
50  import org.orekit.utils.Constants;
51  import org.orekit.utils.IERSConventions;
52  
53  public class OceanTidesTest {
54  
55      @Test
56      public void testDefaultInterpolation() throws OrekitException {
57  
58          IERSConventions conventions = IERSConventions.IERS_2010;
59          Frame eme2000 = FramesFactory.getEME2000();
60          Frame itrf    = FramesFactory.getITRF(conventions, true);
61          TimeScale utc = TimeScalesFactory.getUTC();
62          UT1Scale  ut1 = TimeScalesFactory.getUT1(conventions, true);
63          AstronomicalAmplitudeReader aaReader =
64                  new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
65          DataProvidersManager.getInstance().feed(aaReader.getSupportedNames(), aaReader);
66          Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
67          GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
68                                                                           0.01, FastMath.toRadians(1.0),
69                                                                           OceanLoadDeformationCoefficients.IERS_2010,
70                                                                           map));
71          NormalizedSphericalHarmonicsProvider gravityField =
72                  GravityFieldFactory.getConstantNormalizedProvider(5, 5);
73  
74          // initialization
75          AbsoluteDate date = new AbsoluteDate(1970, 07, 01, 13, 59, 27.816, utc);
76          Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
77                                           FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
78                                           0, PositionAngle.MEAN, eme2000, date,
79                                           gravityField.getMu());
80  
81          AbsoluteDate target = date.shiftedBy(7 * Constants.JULIAN_DAY);
82          ForceModel hf = new HolmesFeatherstoneAttractionModel(itrf, gravityField);
83          SpacecraftState raw = propagate(orbit, target, hf, new OceanTides(itrf, gravityField.getAe(), gravityField.getMu(),
84                         true, Double.NaN, -1,
85                         6, 6, conventions, ut1));
86          SpacecraftState interpolated = propagate(orbit, target, hf, new OceanTides(itrf, gravityField.getAe(), gravityField.getMu(),
87                          6, 6, IERSConventions.IERS_2010, ut1));
88          Assert.assertEquals(0.0,
89                              Vector3D.distance(raw.getPVCoordinates().getPosition(),
90                                                interpolated.getPVCoordinates().getPosition()),
91                              9.9e-6); // threshold would be 3.4e-5 for 30 days propagation
92  
93      }
94  
95      @Test
96      public void testTideEffect1996() throws OrekitException {
97          doTestTideEffect(IERSConventions.IERS_1996, 3.66948, 0.00000);
98      }
99  
100     @Test
101     public void testTideEffect2003() throws OrekitException {
102         doTestTideEffect(IERSConventions.IERS_2003, 3.66941, 0.00000);
103     }
104 
105     @Test
106     public void testTideEffect2010() throws OrekitException {
107         doTestTideEffect(IERSConventions.IERS_2010, 3.66939, 0.08981);
108     }
109 
110     private void doTestTideEffect(IERSConventions conventions, double delta1, double delta2) throws OrekitException {
111 
112         Frame eme2000 = FramesFactory.getEME2000();
113         Frame itrf    = FramesFactory.getITRF(conventions, true);
114         TimeScale utc = TimeScalesFactory.getUTC();
115         UT1Scale  ut1 = TimeScalesFactory.getUT1(conventions, true);
116         AstronomicalAmplitudeReader aaReader =
117                 new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
118         DataProvidersManager.getInstance().feed(aaReader.getSupportedNames(), aaReader);
119         Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
120         GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
121                                                                          0.01, FastMath.toRadians(1.0),
122                                                                          OceanLoadDeformationCoefficients.IERS_2010,
123                                                                          map));
124         NormalizedSphericalHarmonicsProvider gravityField =
125                 GravityFieldFactory.getConstantNormalizedProvider(5, 5);
126 
127         // initialization
128         AbsoluteDate date = new AbsoluteDate(2003, 07, 01, 13, 59, 27.816, utc);
129         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
130                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
131                                          0, PositionAngle.MEAN, eme2000, date,
132                                          gravityField.getMu());
133 
134         AbsoluteDate target = date.shiftedBy(7 * Constants.JULIAN_DAY);
135         ForceModel hf = new HolmesFeatherstoneAttractionModel(itrf, gravityField);
136         SpacecraftState noTides              = propagate(orbit, target, hf);
137         SpacecraftState oceanTidesNoPoleTide = propagate(orbit, target, hf, new OceanTides(itrf, gravityField.getAe(), gravityField.getMu(),
138                         false, SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS,
139                         6, 6, conventions, ut1));
140         SpacecraftState oceanTidesPoleTide = propagate(orbit, target, hf, new OceanTides(itrf, gravityField.getAe(), gravityField.getMu(),
141                           true, SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS,
142                           6, 6, conventions, ut1));
143         Assert.assertEquals(delta1,
144                             Vector3D.distance(noTides.getPVCoordinates().getPosition(),
145                                               oceanTidesNoPoleTide.getPVCoordinates().getPosition()),
146                             0.01);
147         Assert.assertEquals(delta2,
148                             Vector3D.distance(oceanTidesNoPoleTide.getPVCoordinates().getPosition(),
149                                               oceanTidesPoleTide.getPVCoordinates().getPosition()),
150                             0.01);
151 
152     }
153 
154     @Test
155     public void testNoGetParameter() throws OrekitException {
156         AstronomicalAmplitudeReader aaReader =
157                 new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
158         DataProvidersManager.getInstance().feed(aaReader.getSupportedNames(), aaReader);
159         Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
160         GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
161                                                                          0.01, FastMath.toRadians(1.0),
162                                                                          OceanLoadDeformationCoefficients.IERS_2010,
163                                                                          map));
164         ForceModel fm = new OceanTides(FramesFactory.getITRF(IERSConventions.IERS_1996, false),
165                                        Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
166                                        Constants.WGS84_EARTH_MU,
167                                        5, 5, IERSConventions.IERS_1996,
168                                        TimeScalesFactory.getUT1(IERSConventions.IERS_1996, false));
169         Assert.assertTrue(fm.dependsOnPositionOnly());
170         Assert.assertEquals(1, fm.getParametersDrivers().length);
171         try {
172             fm.getParameterDriver("unknown");
173             Assert.fail("an exception should have been thrown");
174         } catch (OrekitException miae) {
175             Assert.assertEquals(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, miae.getSpecifier());
176         }
177     }
178 
179     @Test
180     public void testNoSetParameter() throws OrekitException {
181         AstronomicalAmplitudeReader aaReader =
182                 new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
183         DataProvidersManager.getInstance().feed(aaReader.getSupportedNames(), aaReader);
184         Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
185         GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
186                                                                          0.01, FastMath.toRadians(1.0),
187                                                                          OceanLoadDeformationCoefficients.IERS_2010,
188                                                                          map));
189         ForceModel fm = new OceanTides(FramesFactory.getITRF(IERSConventions.IERS_1996, false),
190                                        Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
191                                        Constants.WGS84_EARTH_MU,
192                                        5, 5, IERSConventions.IERS_1996,
193                                        TimeScalesFactory.getUT1(IERSConventions.IERS_1996, false));
194         Assert.assertEquals(1, fm.getParametersDrivers().length);
195         try {
196             fm.getParameterDriver("unknown").setValue(0.0);
197             Assert.fail("an exception should have been thrown");
198         } catch (OrekitException miae) {
199             Assert.assertEquals(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, miae.getSpecifier());
200         }
201     }
202 
203     private SpacecraftState propagate(Orbit orbit, AbsoluteDate target, ForceModel... forceModels)
204         throws OrekitException {
205         double[][] tolerances = NumericalPropagator.tolerances(10, orbit, OrbitType.KEPLERIAN);
206         AbstractIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 300, tolerances[0], tolerances[1]);
207         NumericalPropagator propagator = new NumericalPropagator(integrator);
208         for (ForceModel forceModel : forceModels) {
209             propagator.addForceModel(forceModel);
210         }
211         propagator.setInitialState(new SpacecraftState(orbit));
212         return propagator.propagate(target);
213     }
214 
215     @Before
216     public void setUp() {
217         Utils.setDataRoot("regular-data:potential/icgem-format:tides");
218     }
219 
220 }