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.forces.gravity;
18  
19  import org.hipparchus.Field;
20  import org.hipparchus.analysis.differentiation.DerivativeStructure;
21  import org.hipparchus.analysis.differentiation.Gradient;
22  import org.hipparchus.analysis.differentiation.GradientField;
23  import org.hipparchus.geometry.euclidean.threed.FieldRotation;
24  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
25  import org.hipparchus.geometry.euclidean.threed.Vector3D;
26  import org.hipparchus.ode.AbstractIntegrator;
27  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
28  import org.hipparchus.util.FastMath;
29  import org.junit.Assert;
30  import org.junit.Before;
31  import org.junit.Test;
32  import org.orekit.Utils;
33  import org.orekit.attitudes.AttitudeProvider;
34  import org.orekit.attitudes.FieldAttitude;
35  import org.orekit.attitudes.LofOffset;
36  import org.orekit.bodies.CelestialBodyFactory;
37  import org.orekit.forces.AbstractLegacyForceModelTest;
38  import org.orekit.forces.ForceModel;
39  import org.orekit.forces.gravity.potential.GravityFieldFactory;
40  import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
41  import org.orekit.frames.Frame;
42  import org.orekit.frames.FramesFactory;
43  import org.orekit.frames.LOFType;
44  import org.orekit.orbits.FieldCartesianOrbit;
45  import org.orekit.orbits.KeplerianOrbit;
46  import org.orekit.orbits.Orbit;
47  import org.orekit.orbits.OrbitType;
48  import org.orekit.orbits.PositionAngle;
49  import org.orekit.propagation.FieldSpacecraftState;
50  import org.orekit.propagation.SpacecraftState;
51  import org.orekit.propagation.numerical.NumericalPropagator;
52  import org.orekit.time.AbsoluteDate;
53  import org.orekit.time.FieldAbsoluteDate;
54  import org.orekit.time.TimeScale;
55  import org.orekit.time.TimeScalesFactory;
56  import org.orekit.time.UT1Scale;
57  import org.orekit.utils.Constants;
58  import org.orekit.utils.IERSConventions;
59  import org.orekit.utils.TimeStampedFieldAngularCoordinates;
60  import org.orekit.utils.TimeStampedFieldPVCoordinates;
61  
62  
63  public class SolidTidesTest extends AbstractLegacyForceModelTest {
64  
65      private static final AttitudeProvider DEFAULT_LAW = Utils.defaultLaw();
66  
67      @Override
68      protected FieldVector3D<DerivativeStructure> accelerationDerivatives(final ForceModel forceModel,
69                                                                           final AbsoluteDate date, final  Frame frame,
70                                                                           final FieldVector3D<DerivativeStructure> position,
71                                                                           final FieldVector3D<DerivativeStructure> velocity,
72                                                                           final FieldRotation<DerivativeStructure> rotation,
73                                                                           final DerivativeStructure mass)
74          {
75          try {
76              java.lang.reflect.Field attractionModelField = SolidTides.class.getDeclaredField("attractionModel");
77              attractionModelField.setAccessible(true);
78              ForceModel attractionModel = (ForceModel) attractionModelField.get(forceModel);
79              double mu = GravityFieldFactory.getConstantNormalizedProvider(5, 5).getMu();
80              Field<DerivativeStructure> field = position.getX().getField();
81              FieldAbsoluteDate<DerivativeStructure> dsDate = new FieldAbsoluteDate<>(field, date);
82              FieldVector3D<DerivativeStructure> zero = FieldVector3D.getZero(field);
83              FieldSpacecraftState<DerivativeStructure> dState =
84                              new FieldSpacecraftState<>(new FieldCartesianOrbit<>(new TimeStampedFieldPVCoordinates<>(dsDate,
85                                                                                                                       position,
86                                                                                                                       velocity,
87                                                                                                                       zero),
88                                                                                   frame, field.getZero().add(mu)),
89                                                        new FieldAttitude<>(frame,
90                                                                        new TimeStampedFieldAngularCoordinates<>(dsDate,
91                                                                                                                 rotation,
92                                                                                                                 zero,
93                                                                                                                 zero)),
94                                                        mass);
95              return attractionModel.acceleration(dState, attractionModel.getParameters(field));
96  
97          } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
98              return null;
99          }
100     }
101 
102     @Override
103     protected FieldVector3D<Gradient> accelerationDerivativesGradient(final ForceModel forceModel,
104                                                                       final AbsoluteDate date, final  Frame frame,
105                                                                       final FieldVector3D<Gradient> position,
106                                                                       final FieldVector3D<Gradient> velocity,
107                                                                       final FieldRotation<Gradient> rotation,
108                                                                       final Gradient mass)
109         {
110         try {
111             java.lang.reflect.Field attractionModelField = SolidTides.class.getDeclaredField("attractionModel");
112             attractionModelField.setAccessible(true);
113             ForceModel attractionModel = (ForceModel) attractionModelField.get(forceModel);
114             double mu = GravityFieldFactory.getConstantNormalizedProvider(5, 5).getMu();
115             final int freeParameters = position.getX().getFreeParameters();
116             Field<Gradient> field = GradientField.getField(freeParameters);
117             FieldAbsoluteDate<Gradient> dsDate = new FieldAbsoluteDate<>(field, date);
118             FieldVector3D<Gradient> zero = FieldVector3D.getZero(field);
119             FieldSpacecraftState<Gradient> dState =
120                             new FieldSpacecraftState<>(new FieldCartesianOrbit<>(new TimeStampedFieldPVCoordinates<>(dsDate,
121                                                                                                                      position,
122                                                                                                                      velocity,
123                                                                                                                      zero),
124                                                                                  frame, field.getZero().add(mu)),
125                                                       new FieldAttitude<>(frame,
126                                                                       new TimeStampedFieldAngularCoordinates<>(dsDate,
127                                                                                                                rotation,
128                                                                                                                zero,
129                                                                                                                zero)),
130                                                       mass);
131             return attractionModel.acceleration(dState, attractionModel.getParameters(field));
132 
133         } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
134             return null;
135         }
136     }
137 
138     @Test
139     public void testDefaultInterpolation() {
140 
141         IERSConventions conventions = IERSConventions.IERS_2010;
142         Frame eme2000 = FramesFactory.getEME2000();
143         Frame itrf    = FramesFactory.getITRF(conventions, true);
144         TimeScale utc = TimeScalesFactory.getUTC();
145         UT1Scale  ut1 = TimeScalesFactory.getUT1(conventions, true);
146         NormalizedSphericalHarmonicsProvider gravityField =
147                 GravityFieldFactory.getConstantNormalizedProvider(5, 5);
148 
149         // initialization
150         AbsoluteDate date = new AbsoluteDate(1970, 07, 01, 13, 59, 27.816, utc);
151         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
152                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
153                                          0, PositionAngle.MEAN, eme2000, date,
154                                          gravityField.getMu());
155 
156         AbsoluteDate target = date.shiftedBy(7 * Constants.JULIAN_DAY);
157         ForceModel hf = new HolmesFeatherstoneAttractionModel(itrf, gravityField);
158         SpacecraftState raw = propagate(orbit, target, hf,
159                                         new SolidTides(itrf, gravityField.getAe(), gravityField.getMu(),
160                                                        gravityField.getTideSystem(), true, Double.NaN, -1,
161                                                        conventions, ut1,
162                                                        CelestialBodyFactory.getSun(),
163                                                        CelestialBodyFactory.getMoon()));
164         SpacecraftState interpolated = propagate(orbit, target, hf,
165                                                  new SolidTides(itrf, gravityField.getAe(), gravityField.getMu(),
166                                                                 gravityField.getTideSystem(),
167                                                                 conventions, ut1,
168                                                                 CelestialBodyFactory.getSun(),
169                                                                 CelestialBodyFactory.getMoon()));
170         Assert.assertEquals(0.0,
171                             Vector3D.distance(raw.getPVCoordinates().getPosition(),
172                                               interpolated.getPVCoordinates().getPosition()),
173                             2.0e-5); // threshold would be 1.2e-3 for 30 days propagation
174 
175     }
176 
177     @Test
178     public void testTideEffect1996() {
179         Frame eme2000 = FramesFactory.getEME2000();
180         TimeScale utc = TimeScalesFactory.getUTC();
181         AbsoluteDate date = new AbsoluteDate(2003, 07, 01, 13, 59, 27.816, utc);
182         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
183                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
184                                          0, PositionAngle.MEAN, eme2000, date,
185                                          Constants.EIGEN5C_EARTH_MU);
186         doTestTideEffect(orbit, IERSConventions.IERS_1996, 44.09481, 0.00000);
187     }
188 
189     @Test
190     public void testTideEffect2003WithinAnnualPoleRange() {
191         Frame eme2000 = FramesFactory.getEME2000();
192         TimeScale utc = TimeScalesFactory.getUTC();
193         AbsoluteDate date = new AbsoluteDate(1969, 07, 01, 13, 59, 27.816, utc);
194         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
195                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
196                                          0, PositionAngle.MEAN, eme2000, date,
197                                          Constants.EIGEN5C_EARTH_MU);
198         doTestTideEffect(orbit, IERSConventions.IERS_2003, 73.14011, 0.87360);
199     }
200 
201     @Test
202     public void testTideEffect2003AfterAnnualPoleRange() {
203         Frame eme2000 = FramesFactory.getEME2000();
204         TimeScale utc = TimeScalesFactory.getUTC();
205         AbsoluteDate date = new AbsoluteDate(2003, 07, 01, 13, 59, 27.816, utc);
206         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
207                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
208                                          0, PositionAngle.MEAN, eme2000, date,
209                                          Constants.EIGEN5C_EARTH_MU);
210         doTestTideEffect(orbit, IERSConventions.IERS_2003, 44.24999, 0.61752);
211     }
212 
213     @Test
214     public void testTideEffect2010BeforePoleModelChange() {
215         Frame eme2000 = FramesFactory.getEME2000();
216         TimeScale utc = TimeScalesFactory.getUTC();
217         AbsoluteDate date = new AbsoluteDate(2003, 07, 01, 13, 59, 27.816, utc);
218         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
219                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
220                                          0, PositionAngle.MEAN, eme2000, date,
221                                          Constants.EIGEN5C_EARTH_MU);
222         doTestTideEffect(orbit, IERSConventions.IERS_2010, 44.25001, 0.70710);
223     }
224 
225     @Test
226     public void testTideEffect2010AfterModelChange() {
227         Frame eme2000 = FramesFactory.getEME2000();
228         TimeScale utc = TimeScalesFactory.getUTC();
229         AbsoluteDate date = new AbsoluteDate(2964, 8, 12, 11, 30, 00.000, utc);
230         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
231                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
232                                          0, PositionAngle.MEAN, eme2000, date,
233                                          Constants.EIGEN5C_EARTH_MU);
234         doTestTideEffect(orbit, IERSConventions.IERS_2010, 24.02815, 30.37047);
235     }
236 
237     @Test
238     public void testStateJacobianVs80ImplementationNoPoleTide()
239         {
240         Frame eme2000 = FramesFactory.getEME2000();
241         TimeScale utc = TimeScalesFactory.getUTC();
242         AbsoluteDate date = new AbsoluteDate(2964, 8, 12, 11, 30, 00.000, utc);
243         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
244                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
245                                          0, PositionAngle.MEAN, eme2000, date,
246                                          Constants.EIGEN5C_EARTH_MU);
247         Frame itrf    = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
248         UT1Scale  ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
249         NormalizedSphericalHarmonicsProvider gravityField =
250                         GravityFieldFactory.getConstantNormalizedProvider(5, 5);
251 
252         ForceModel forceModel = new SolidTides(itrf, gravityField.getAe(), gravityField.getMu(),
253                                                gravityField.getTideSystem(), false,
254                                                SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS,
255                                                IERSConventions.IERS_2010, ut1,
256                                                CelestialBodyFactory.getSun(),
257                                                CelestialBodyFactory.getMoon());
258         Assert.assertTrue(forceModel.dependsOnPositionOnly());
259 
260         checkStateJacobianVs80Implementation(new SpacecraftState(orbit), forceModel,
261                                              new LofOffset(orbit.getFrame(), LOFType.LVLH_CCSDS),
262                                              2.0e-15, false);
263 
264     }
265 
266     @Test
267     public void testStateJacobianVs80ImplementationGradientNoPoleTide()
268         {
269         Frame eme2000 = FramesFactory.getEME2000();
270         TimeScale utc = TimeScalesFactory.getUTC();
271         AbsoluteDate date = new AbsoluteDate(2964, 8, 12, 11, 30, 00.000, utc);
272         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
273                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
274                                          0, PositionAngle.MEAN, eme2000, date,
275                                          Constants.EIGEN5C_EARTH_MU);
276         Frame itrf    = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
277         UT1Scale  ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
278         NormalizedSphericalHarmonicsProvider gravityField =
279                         GravityFieldFactory.getConstantNormalizedProvider(5, 5);
280 
281         ForceModel forceModel = new SolidTides(itrf, gravityField.getAe(), gravityField.getMu(),
282                                                gravityField.getTideSystem(), false,
283                                                SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS,
284                                                IERSConventions.IERS_2010, ut1,
285                                                CelestialBodyFactory.getSun(),
286                                                CelestialBodyFactory.getMoon());
287         Assert.assertTrue(forceModel.dependsOnPositionOnly());
288 
289         checkStateJacobianVs80ImplementationGradient(new SpacecraftState(orbit), forceModel,
290                                              new LofOffset(orbit.getFrame(), LOFType.LVLH_CCSDS),
291                                              2.0e-15, false);
292 
293     }
294 
295     @Test
296     public void testStateJacobianVs80ImplementationPoleTide()
297         {
298         Frame eme2000 = FramesFactory.getEME2000();
299         TimeScale utc = TimeScalesFactory.getUTC();
300         AbsoluteDate date = new AbsoluteDate(2964, 8, 12, 11, 30, 00.000, utc);
301         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
302                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
303                                          0, PositionAngle.MEAN, eme2000, date,
304                                          Constants.EIGEN5C_EARTH_MU);
305         Frame itrf    = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
306         UT1Scale  ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
307         NormalizedSphericalHarmonicsProvider gravityField =
308                         GravityFieldFactory.getConstantNormalizedProvider(5, 5);
309 
310         ForceModel forceModel = new SolidTides(itrf, gravityField.getAe(), gravityField.getMu(),
311                                                gravityField.getTideSystem(), true,
312                                                SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS,
313                                                IERSConventions.IERS_2010, ut1,
314                                                CelestialBodyFactory.getSun(),
315                                                CelestialBodyFactory.getMoon());
316 
317         checkStateJacobianVs80Implementation(new SpacecraftState(orbit), forceModel,
318                                              new LofOffset(orbit.getFrame(), LOFType.LVLH_CCSDS),
319                                              2.0e-15, false);
320 
321     }
322 
323     @Test
324     public void testStateJacobianVs80ImplementationGradientPoleTide()
325         {
326         Frame eme2000 = FramesFactory.getEME2000();
327         TimeScale utc = TimeScalesFactory.getUTC();
328         AbsoluteDate date = new AbsoluteDate(2964, 8, 12, 11, 30, 00.000, utc);
329         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
330                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
331                                          0, PositionAngle.MEAN, eme2000, date,
332                                          Constants.EIGEN5C_EARTH_MU);
333         Frame itrf    = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
334         UT1Scale  ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
335         NormalizedSphericalHarmonicsProvider gravityField =
336                         GravityFieldFactory.getConstantNormalizedProvider(5, 5);
337 
338         ForceModel forceModel = new SolidTides(itrf, gravityField.getAe(), gravityField.getMu(),
339                                                gravityField.getTideSystem(), true,
340                                                SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS,
341                                                IERSConventions.IERS_2010, ut1,
342                                                CelestialBodyFactory.getSun(),
343                                                CelestialBodyFactory.getMoon());
344 
345         checkStateJacobianVs80ImplementationGradient(new SpacecraftState(orbit), forceModel,
346                                              new LofOffset(orbit.getFrame(), LOFType.LVLH_CCSDS),
347                                              2.0e-15, false);
348 
349     }
350 
351     @Test
352     public void testStateJacobianVsFiniteDifferencesNoPoleTide()
353         {
354         Frame eme2000 = FramesFactory.getEME2000();
355         TimeScale utc = TimeScalesFactory.getUTC();
356         AbsoluteDate date = new AbsoluteDate(2964, 8, 12, 11, 30, 00.000, utc);
357         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
358                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
359                                          0, PositionAngle.MEAN, eme2000, date,
360                                          Constants.EIGEN5C_EARTH_MU);
361         Frame itrf    = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
362         UT1Scale  ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
363         NormalizedSphericalHarmonicsProvider gravityField =
364                         GravityFieldFactory.getConstantNormalizedProvider(5, 5);
365 
366         ForceModel forceModel = new SolidTides(itrf, gravityField.getAe(), gravityField.getMu(),
367                                                gravityField.getTideSystem(), false,
368                                                SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS,
369                                                IERSConventions.IERS_2010, ut1,
370                                                CelestialBodyFactory.getSun(),
371                                                CelestialBodyFactory.getMoon());
372 
373         checkStateJacobianVsFiniteDifferences(new SpacecraftState(orbit), forceModel, DEFAULT_LAW,
374                                               10.0, 2.0e-10, false);
375 
376     }
377 
378     @Test
379     public void testStateJacobianVsFiniteDifferencesGradientNoPoleTide()
380         {
381         Frame eme2000 = FramesFactory.getEME2000();
382         TimeScale utc = TimeScalesFactory.getUTC();
383         AbsoluteDate date = new AbsoluteDate(2964, 8, 12, 11, 30, 00.000, utc);
384         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
385                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
386                                          0, PositionAngle.MEAN, eme2000, date,
387                                          Constants.EIGEN5C_EARTH_MU);
388         Frame itrf    = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
389         UT1Scale  ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
390         NormalizedSphericalHarmonicsProvider gravityField =
391                         GravityFieldFactory.getConstantNormalizedProvider(5, 5);
392 
393         ForceModel forceModel = new SolidTides(itrf, gravityField.getAe(), gravityField.getMu(),
394                                                gravityField.getTideSystem(), false,
395                                                SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS,
396                                                IERSConventions.IERS_2010, ut1,
397                                                CelestialBodyFactory.getSun(),
398                                                CelestialBodyFactory.getMoon());
399 
400         checkStateJacobianVsFiniteDifferencesGradient(new SpacecraftState(orbit), forceModel, DEFAULT_LAW,
401                                               10.0, 2.0e-10, false);
402 
403     }
404 
405     @Test
406     public void testStateJacobianVsFiniteDifferencesPoleTide()
407         {
408         Frame eme2000 = FramesFactory.getEME2000();
409         TimeScale utc = TimeScalesFactory.getUTC();
410         AbsoluteDate date = new AbsoluteDate(2964, 8, 12, 11, 30, 00.000, utc);
411         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
412                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
413                                          0, PositionAngle.MEAN, eme2000, date,
414                                          Constants.EIGEN5C_EARTH_MU);
415         Frame itrf    = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
416         UT1Scale  ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
417         NormalizedSphericalHarmonicsProvider gravityField =
418                         GravityFieldFactory.getConstantNormalizedProvider(5, 5);
419 
420         ForceModel forceModel = new SolidTides(itrf, gravityField.getAe(), gravityField.getMu(),
421                                                gravityField.getTideSystem(), true,
422                                                SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS,
423                                                IERSConventions.IERS_2010, ut1,
424                                                CelestialBodyFactory.getSun(),
425                                                CelestialBodyFactory.getMoon());
426 
427         checkStateJacobianVsFiniteDifferences(new SpacecraftState(orbit), forceModel, DEFAULT_LAW,
428                                               10.0, 2.0e-10, false);
429 
430     }
431 
432     @Test
433     public void testStateJacobianVsFiniteDifferencesGradientPoleTide()
434         {
435         Frame eme2000 = FramesFactory.getEME2000();
436         TimeScale utc = TimeScalesFactory.getUTC();
437         AbsoluteDate date = new AbsoluteDate(2964, 8, 12, 11, 30, 00.000, utc);
438         Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
439                                          FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
440                                          0, PositionAngle.MEAN, eme2000, date,
441                                          Constants.EIGEN5C_EARTH_MU);
442         Frame itrf    = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
443         UT1Scale  ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
444         NormalizedSphericalHarmonicsProvider gravityField =
445                         GravityFieldFactory.getConstantNormalizedProvider(5, 5);
446 
447         ForceModel forceModel = new SolidTides(itrf, gravityField.getAe(), gravityField.getMu(),
448                                                gravityField.getTideSystem(), true,
449                                                SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS,
450                                                IERSConventions.IERS_2010, ut1,
451                                                CelestialBodyFactory.getSun(),
452                                                CelestialBodyFactory.getMoon());
453 
454         checkStateJacobianVsFiniteDifferencesGradient(new SpacecraftState(orbit), forceModel, DEFAULT_LAW,
455                                               10.0, 2.0e-10, false);
456 
457     }
458 
459     private void doTestTideEffect(Orbit orbit, IERSConventions conventions, double delta1, double delta2) {
460 
461         Frame itrf    = FramesFactory.getITRF(conventions, true);
462         UT1Scale  ut1 = TimeScalesFactory.getUT1(conventions, true);
463         NormalizedSphericalHarmonicsProvider gravityField =
464                 GravityFieldFactory.getConstantNormalizedProvider(5, 5);
465 
466         // initialization
467 
468         AbsoluteDate target = orbit.getDate().shiftedBy(7 * Constants.JULIAN_DAY);
469         ForceModel hf = new HolmesFeatherstoneAttractionModel(itrf, gravityField);
470         SpacecraftState noTides              = propagate(orbit, target, hf);
471         SpacecraftState solidTidesNoPoleTide = propagate(orbit, target, hf,
472                                                          new SolidTides(itrf, gravityField.getAe(), gravityField.getMu(),
473                                                                         gravityField.getTideSystem(), false,
474                                                                         SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS,
475                                                                         conventions, ut1,
476                                                                         CelestialBodyFactory.getSun(),
477                                                                         CelestialBodyFactory.getMoon()));
478         SpacecraftState solidTidesPoleTide = propagate(orbit, target, hf,
479                                                        new SolidTides(itrf, gravityField.getAe(), gravityField.getMu(),
480                                                                       gravityField.getTideSystem(), true,
481                                                                       SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS,
482                                                                       conventions, ut1,
483                                                                       CelestialBodyFactory.getSun(),
484                                                                       CelestialBodyFactory.getMoon()));
485         Assert.assertEquals(delta1,
486                             Vector3D.distance(noTides.getPVCoordinates().getPosition(),
487                                               solidTidesNoPoleTide.getPVCoordinates().getPosition()),
488                             0.01);
489         Assert.assertEquals(delta2,
490                             Vector3D.distance(solidTidesNoPoleTide.getPVCoordinates().getPosition(),
491                                               solidTidesPoleTide.getPVCoordinates().getPosition()),
492                             0.01);
493 
494     }
495 
496     private SpacecraftState propagate(Orbit orbit, AbsoluteDate target, ForceModel... forceModels)
497         {
498         double[][] tolerances = NumericalPropagator.tolerances(10, orbit, OrbitType.KEPLERIAN);
499         AbstractIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 300, tolerances[0], tolerances[1]);
500         NumericalPropagator propagator = new NumericalPropagator(integrator);
501         for (ForceModel forceModel : forceModels) {
502             propagator.addForceModel(forceModel);
503         }
504         propagator.setInitialState(new SpacecraftState(orbit));
505         return propagator.propagate(target);
506     }
507 
508     @Before
509     public void setUp() {
510         Utils.setDataRoot("regular-data:potential/icgem-format");
511     }
512 
513 }