1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.gravity;
18
19 import java.util.List;
20 import java.util.Map;
21 import java.util.stream.Collectors;
22
23 import org.hipparchus.geometry.euclidean.threed.Vector3D;
24 import org.hipparchus.ode.AbstractIntegrator;
25 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
26 import org.hipparchus.util.Binary64;
27 import org.hipparchus.util.Binary64Field;
28 import org.hipparchus.util.FastMath;
29 import org.junit.jupiter.api.Assertions;
30 import org.junit.jupiter.api.BeforeEach;
31 import org.junit.jupiter.api.Test;
32 import org.orekit.Utils;
33 import org.orekit.bodies.OneAxisEllipsoid;
34 import org.orekit.data.DataContext;
35 import org.orekit.errors.OrekitException;
36 import org.orekit.errors.OrekitMessages;
37 import org.orekit.forces.ForceModel;
38 import org.orekit.forces.gravity.potential.AstronomicalAmplitudeReader;
39 import org.orekit.forces.gravity.potential.FESCHatEpsilonReader;
40 import org.orekit.forces.gravity.potential.GravityFieldFactory;
41 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
42 import org.orekit.forces.gravity.potential.OceanLoadDeformationCoefficients;
43 import org.orekit.frames.Frame;
44 import org.orekit.frames.FramesFactory;
45 import org.orekit.models.earth.ReferenceEllipsoid;
46 import org.orekit.orbits.CartesianOrbit;
47 import org.orekit.orbits.KeplerianOrbit;
48 import org.orekit.orbits.Orbit;
49 import org.orekit.orbits.OrbitType;
50 import org.orekit.orbits.PositionAngleType;
51 import org.orekit.propagation.SpacecraftState;
52 import org.orekit.propagation.ToleranceProvider;
53 import org.orekit.propagation.events.DateDetector;
54 import org.orekit.propagation.events.EventDetector;
55 import org.orekit.propagation.events.EventDetectorsProvider;
56 import org.orekit.propagation.events.FieldDateDetector;
57 import org.orekit.propagation.events.FieldEventDetector;
58 import org.orekit.propagation.numerical.NumericalPropagator;
59 import org.orekit.time.AbsoluteDate;
60 import org.orekit.time.FieldAbsoluteDate;
61 import org.orekit.time.TimeScale;
62 import org.orekit.time.TimeScalesFactory;
63 import org.orekit.time.TimeStamped;
64 import org.orekit.time.UT1Scale;
65 import org.orekit.utils.Constants;
66 import org.orekit.utils.IERSConventions;
67 import org.orekit.utils.PVCoordinates;
68 import org.orekit.utils.ParameterDriver;
69
70 public class OceanTidesTest {
71
72
73
74
75 @Test
76 void testIssue1055() {
77
78
79 AstronomicalAmplitudeReader aaReader =
80 new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
81 DataContext.getDefault().getDataProvidersManager().feed(aaReader.getSupportedNames(), aaReader);
82 FESCHatEpsilonReader otReader = new FESCHatEpsilonReader("fes2004-7x7.dat", 0.01, FastMath.toRadians(1.0),
83 OceanLoadDeformationCoefficients.IERS_2010,
84 aaReader.getAstronomicalAmplitudesMap());
85 GravityFieldFactory.addOceanTidesReader(otReader);
86 NormalizedSphericalHarmonicsProvider gravityProvider = GravityFieldFactory.getNormalizedProvider(5, 5);
87 Frame earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
88 UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
89 AbsoluteDate epoch = new AbsoluteDate("2021-03-29T23:59:42.000", TimeScalesFactory.getUTC());
90 CartesianOrbit orbit = new CartesianOrbit(new PVCoordinates(new Vector3D(531.7545638326519,-718035.0404060499, 6835039.371647774),
91 new Vector3D(1165.678605794716, -7476.55626016528, -789.2919601239232)),
92 FramesFactory.getGCRF(), epoch, gravityProvider.getMu());
93 SpacecraftState state = new SpacecraftState(orbit);
94
95
96 Vector3D acc3 = oceanTidesAcceleration(3, ut1, earthFrame, state, gravityProvider);
97 Assertions.assertEquals(3, otReader.getMaxAvailableDegree());
98 Assertions.assertEquals(3, otReader.getMaxParseDegree());
99 Assertions.assertEquals(3, otReader.getMaxAvailableOrder());
100 Assertions.assertEquals(3, otReader.getMaxParseOrder());
101 Vector3D acc6 = oceanTidesAcceleration(6, ut1, earthFrame, state, gravityProvider);
102 Assertions.assertEquals(6, otReader.getMaxAvailableDegree());
103 Assertions.assertEquals(6, otReader.getMaxParseDegree());
104 Assertions.assertEquals(6, otReader.getMaxAvailableOrder());
105 Assertions.assertEquals(6, otReader.getMaxParseOrder());
106 Vector3D acc2 = oceanTidesAcceleration(2, ut1, earthFrame, state, gravityProvider);
107 Assertions.assertEquals(2, otReader.getMaxAvailableDegree());
108 Assertions.assertEquals(2, otReader.getMaxParseDegree());
109 Assertions.assertEquals(2, otReader.getMaxAvailableOrder());
110 Assertions.assertEquals(2, otReader.getMaxParseOrder());
111
112
113 Assertions.assertEquals(6.004665607951679E-9, acc3.getX());
114 Assertions.assertEquals(2.379362826744579E-8, acc3.getY());
115 Assertions.assertEquals(1.2474166439853716E-9, acc3.getZ());
116
117 Assertions.assertEquals(-4.443030668255491E-9, acc6.getX());
118 Assertions.assertEquals(1.7782200620821885E-9, acc6.getY());
119 Assertions.assertEquals(1.3663400321897177E-9, acc6.getZ());
120
121 Assertions.assertEquals(-1.6754788749741251E-9, acc2.getX());
122 Assertions.assertEquals(1.1119622068766252E-8, acc2.getY());
123 Assertions.assertEquals(8.891976691804308E-9, acc2.getZ());
124
125 }
126
127 @Test
128 void testDefaultInterpolation() {
129
130 IERSConventions conventions = IERSConventions.IERS_2010;
131 Frame eme2000 = FramesFactory.getEME2000();
132 Frame itrf = FramesFactory.getITRF(conventions, true);
133 TimeScale utc = TimeScalesFactory.getUTC();
134 UT1Scale ut1 = TimeScalesFactory.getUT1(conventions, true);
135 AstronomicalAmplitudeReader aaReader =
136 new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
137 DataContext.getDefault().getDataProvidersManager().feed(aaReader.getSupportedNames(), aaReader);
138 Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
139 GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
140 0.01, FastMath.toRadians(1.0),
141 OceanLoadDeformationCoefficients.IERS_2010,
142 map));
143 NormalizedSphericalHarmonicsProvider gravityField =
144 GravityFieldFactory.getNormalizedProvider(5, 5);
145
146
147 AbsoluteDate date = new AbsoluteDate(1970, 07, 01, 13, 59, 27.816, utc);
148 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
149 FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
150 0, PositionAngleType.MEAN, eme2000, date,
151 gravityField.getMu());
152
153 AbsoluteDate target = date.shiftedBy(7 * Constants.JULIAN_DAY);
154 ForceModel hf = new HolmesFeatherstoneAttractionModel(itrf, gravityField);
155 SpacecraftState raw = propagate(orbit, target, hf, new OceanTides(itrf, gravityField.getAe(), gravityField.getMu(),
156 true, Double.NaN, -1,
157 6, 6, conventions, ut1));
158 SpacecraftState interpolated = propagate(orbit, target, hf, new OceanTides(itrf, gravityField.getAe(), gravityField.getMu(),
159 6, 6, IERSConventions.IERS_2010, ut1));
160 Assertions.assertEquals(0.0,
161 Vector3D.distance(raw.getPosition(),
162 interpolated.getPosition()),
163 9.9e-6);
164
165 }
166
167 @Test
168 void testTideEffect1996() {
169 doTestTideEffect(IERSConventions.IERS_1996, 3.66948, 0.00000);
170 }
171
172 @Test
173 void testTideEffect2003() {
174 doTestTideEffect(IERSConventions.IERS_2003, 3.66941, 0.00000);
175 }
176
177 @Test
178 void testTideEffect2010() {
179 doTestTideEffect(IERSConventions.IERS_2010, 3.66939, 0.08981);
180 }
181
182 private void doTestTideEffect(IERSConventions conventions, double delta1, double delta2) {
183
184 Frame eme2000 = FramesFactory.getEME2000();
185 Frame itrf = FramesFactory.getITRF(conventions, true);
186 TimeScale utc = TimeScalesFactory.getUTC();
187 UT1Scale ut1 = TimeScalesFactory.getUT1(conventions, true);
188 AstronomicalAmplitudeReader aaReader =
189 new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
190 DataContext.getDefault().getDataProvidersManager().feed(aaReader.getSupportedNames(), aaReader);
191 Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
192 GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
193 0.01, FastMath.toRadians(1.0),
194 OceanLoadDeformationCoefficients.IERS_2010,
195 map));
196 NormalizedSphericalHarmonicsProvider gravityField =
197 GravityFieldFactory.getNormalizedProvider(5, 5);
198
199
200 AbsoluteDate date = new AbsoluteDate(2003, 07, 01, 13, 59, 27.816, utc);
201 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, FastMath.toRadians(98.7),
202 FastMath.toRadians(93.0), FastMath.toRadians(15.0 * 22.5),
203 0, PositionAngleType.MEAN, eme2000, date,
204 gravityField.getMu());
205
206 AbsoluteDate target = date.shiftedBy(7 * Constants.JULIAN_DAY);
207 ForceModel hf = new HolmesFeatherstoneAttractionModel(itrf, gravityField);
208 SpacecraftState noTides = propagate(orbit, target, hf);
209 SpacecraftState oceanTidesNoPoleTide = propagate(orbit, target, hf, new OceanTides(itrf, gravityField.getAe(), gravityField.getMu(),
210 false, SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS,
211 6, 6, conventions, ut1));
212 SpacecraftState oceanTidesPoleTide = propagate(orbit, target, hf, new OceanTides(itrf, gravityField.getAe(), gravityField.getMu(),
213 true, SolidTides.DEFAULT_STEP, SolidTides.DEFAULT_POINTS,
214 6, 6, conventions, ut1));
215 Assertions.assertEquals(delta1,
216 Vector3D.distance(noTides.getPosition(),
217 oceanTidesNoPoleTide.getPosition()),
218 0.01);
219 Assertions.assertEquals(delta2,
220 Vector3D.distance(oceanTidesNoPoleTide.getPosition(),
221 oceanTidesPoleTide.getPosition()),
222 0.01);
223
224 }
225
226 @Test
227 void testNoGetParameter() {
228 AstronomicalAmplitudeReader aaReader =
229 new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
230 DataContext.getDefault().getDataProvidersManager().feed(aaReader.getSupportedNames(), aaReader);
231 Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
232 GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
233 0.01, FastMath.toRadians(1.0),
234 OceanLoadDeformationCoefficients.IERS_2010,
235 map));
236 ForceModel fm = new OceanTides(FramesFactory.getITRF(IERSConventions.IERS_1996, false),
237 Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
238 Constants.WGS84_EARTH_MU,
239 5, 5, IERSConventions.IERS_1996,
240 TimeScalesFactory.getUT1(IERSConventions.IERS_1996, false));
241 Assertions.assertTrue(fm.dependsOnPositionOnly());
242 Assertions.assertEquals(1, fm.getParametersDrivers().size());
243 try {
244 fm.getParameterDriver("unknown");
245 Assertions.fail("an exception should have been thrown");
246 } catch (OrekitException miae) {
247 Assertions.assertEquals(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, miae.getSpecifier());
248 }
249 }
250
251 @Test
252 void testNoSetParameter() {
253 AstronomicalAmplitudeReader aaReader =
254 new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
255 DataContext.getDefault().getDataProvidersManager().feed(aaReader.getSupportedNames(), aaReader);
256 Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
257 GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
258 0.01, FastMath.toRadians(1.0),
259 OceanLoadDeformationCoefficients.IERS_2010,
260 map));
261 ForceModel fm = new OceanTides(FramesFactory.getITRF(IERSConventions.IERS_1996, false),
262 Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
263 Constants.WGS84_EARTH_MU,
264 5, 5, IERSConventions.IERS_1996,
265 TimeScalesFactory.getUT1(IERSConventions.IERS_1996, false));
266 Assertions.assertEquals(1, fm.getParametersDrivers().size());
267 try {
268 fm.getParameterDriver("unknown").setValue(0.0);
269 Assertions.fail("an exception should have been thrown");
270 } catch (OrekitException miae) {
271 Assertions.assertEquals(OrekitMessages.UNSUPPORTED_PARAMETER_NAME, miae.getSpecifier());
272 }
273 }
274
275
276
277
278 @Test
279 void testGetEventDetectors() {
280
281
282
283
284 final IERSConventions conventions = IERSConventions.IERS_2010;
285 final Frame itrf = FramesFactory.getITRF(conventions, true);
286 final AbsoluteDate t0 = AbsoluteDate.ARBITRARY_EPOCH;
287
288 final OneAxisEllipsoid body = ReferenceEllipsoid.getIers2010(itrf);
289
290 final NormalizedSphericalHarmonicsProvider gravityField =
291 GravityFieldFactory.getNormalizedProvider(5, 5);
292
293
294 final AstronomicalAmplitudeReader aaReader =
295 new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
296 DataContext.getDefault().getDataProvidersManager().feed(aaReader.getSupportedNames(), aaReader);
297 final Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
298 GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
299 0.01, FastMath.toRadians(1.0),
300 OceanLoadDeformationCoefficients.IERS_2010,
301 map));
302
303
304 final ForceModel oceanTidesModel = new OceanTides(body.getBodyFrame(),
305 gravityField.getAe(), gravityField.getMu(),
306 5, 5, conventions,
307 TimeScalesFactory.getUT1(conventions, true));
308
309
310 List<EventDetector> detectors = oceanTidesModel.getEventDetectors().collect(Collectors.toList());
311 List<FieldEventDetector<Binary64>> fieldDetectors = oceanTidesModel.getFieldEventDetectors(Binary64Field.getInstance()).collect(Collectors.toList());
312
313
314 Assertions.assertTrue(detectors.isEmpty());
315 Assertions.assertTrue(fieldDetectors.isEmpty());
316
317
318 final List<ParameterDriver> drivers = oceanTidesModel.getParametersDrivers();
319
320 for (final ParameterDriver driver : drivers) {
321 driver.addSpanAtDate(t0);
322 }
323
324 detectors = oceanTidesModel.getEventDetectors().collect(Collectors.toList());
325 DateDetector dateDetector = (DateDetector) detectors.get(0);
326 List<TimeStamped> dates = dateDetector.getDates();
327
328 fieldDetectors = oceanTidesModel.getFieldEventDetectors(Binary64Field.getInstance()).collect(Collectors.toList());
329 FieldDateDetector<Binary64> fieldDateDetector = (FieldDateDetector<Binary64>) fieldDetectors.get(0);
330 FieldAbsoluteDate<Binary64> fieldDate = fieldDateDetector.getDate();
331
332
333 Assertions.assertFalse(detectors.isEmpty());
334 Assertions.assertEquals(1, detectors.size());
335 Assertions.assertTrue(detectors.get(0) instanceof DateDetector);
336
337 Assertions.assertEquals(1, dates.size());
338 Assertions.assertEquals(0., dates.get(0).durationFrom(t0), 0.);
339
340 Assertions.assertFalse(fieldDetectors.isEmpty());
341 Assertions.assertEquals(1, fieldDetectors.size());
342 Assertions.assertTrue(fieldDetectors.get(0) instanceof FieldDateDetector);
343 Assertions.assertEquals(0., fieldDate.durationFrom(t0).getReal(), 0.);
344 }
345
346 private SpacecraftState propagate(Orbit orbit, AbsoluteDate target, ForceModel... forceModels)
347 {
348 double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(orbit, OrbitType.KEPLERIAN);
349 AbstractIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 300, tolerances[0], tolerances[1]);
350 NumericalPropagator propagator = new NumericalPropagator(integrator);
351 for (ForceModel forceModel : forceModels) {
352 propagator.addForceModel(forceModel);
353 }
354 propagator.setInitialState(new SpacecraftState(orbit));
355 return propagator.propagate(target);
356 }
357
358 private Vector3D oceanTidesAcceleration(int degree, UT1Scale ut1, Frame earthFrame, SpacecraftState state,
359 NormalizedSphericalHarmonicsProvider gravityProvider) {
360 OceanTides force = new OceanTides(earthFrame, gravityProvider.getAe(), gravityProvider.getMu(), true,
361 600.0, 12, degree, degree, IERSConventions.IERS_2010, ut1);
362 return force.acceleration(state, force.getParameters());
363 }
364
365 @BeforeEach
366 public void setUp() {
367 Utils.setDataRoot("regular-data:potential/icgem-format:tides");
368 }
369
370 }