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 org.hamcrest.MatcherAssert;
20 import org.hamcrest.Matchers;
21 import org.hipparchus.CalculusFieldElement;
22 import org.hipparchus.stat.descriptive.StreamingStatistics;
23 import org.hipparchus.util.FastMath;
24 import org.junit.jupiter.api.Assertions;
25 import org.junit.jupiter.api.BeforeEach;
26 import org.junit.jupiter.api.Test;
27 import org.orekit.Utils;
28 import org.orekit.bodies.CelestialBodyFactory;
29 import org.orekit.data.BodiesElements;
30 import org.orekit.data.FundamentalNutationArguments;
31 import org.orekit.data.PoissonSeries;
32 import org.orekit.data.PoissonSeriesParser;
33 import org.orekit.errors.OrekitInternalError;
34 import org.orekit.forces.gravity.potential.CachedNormalizedSphericalHarmonicsProvider;
35 import org.orekit.forces.gravity.potential.GravityFieldFactory;
36 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
37 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics;
38 import org.orekit.forces.gravity.potential.TideSystem;
39 import org.orekit.frames.Frame;
40 import org.orekit.frames.FramesFactory;
41 import org.orekit.time.AbsoluteDate;
42 import org.orekit.time.FieldAbsoluteDate;
43 import org.orekit.time.TimeScalarFunction;
44 import org.orekit.time.TimeScale;
45 import org.orekit.time.TimeScalesFactory;
46 import org.orekit.time.TimeVectorFunction;
47 import org.orekit.time.UT1Scale;
48 import org.orekit.utils.Constants;
49 import org.orekit.utils.IERSConventions;
50 import org.orekit.utils.LoveNumbers;
51 import org.orekit.utils.OrekitConfiguration;
52
53 import java.lang.reflect.Field;
54 import java.lang.reflect.InvocationTargetException;
55 import java.lang.reflect.Method;
56 import java.util.Arrays;
57
58
59 public class SolidTidesFieldTest {
60
61 @Test
62 void testConventions2003() throws NoSuchFieldException, IllegalAccessException {
63
64 UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, false);
65 SolidTidesField tidesField =
66 new SolidTidesField(IERSConventions.IERS_2003.getLoveNumbers(),
67 IERSConventions.IERS_2003.getTideFrequencyDependenceFunction(ut1),
68 IERSConventions.IERS_2003.getPermanentTide(),
69 IERSConventions.IERS_2003.getSolidPoleTide(ut1.getEOPHistory()),
70 FramesFactory.getITRF(IERSConventions.IERS_2003, false),
71 Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_MU,
72 TideSystem.ZERO_TIDE, CelestialBodyFactory.getSun(), CelestialBodyFactory.getMoon());
73
74 Field fieldReal = tidesField.getClass().getDeclaredField("love");
75 fieldReal.setAccessible(true);
76 LoveNumbers love = (LoveNumbers) fieldReal.get(tidesField);
77
78 Assertions.assertEquals(0.30190, love.getReal(2, 0), 1.0e-10);
79 Assertions.assertEquals(0.29830, love.getReal(2, 1), 1.0e-10);
80 Assertions.assertEquals(0.30102, love.getReal(2, 2), 1.0e-10);
81 Assertions.assertEquals(0.093, love.getReal(3, 0), 1.0e-10);
82 Assertions.assertEquals(0.093, love.getReal(3, 1), 1.0e-10);
83 Assertions.assertEquals(0.093, love.getReal(3, 2), 1.0e-10);
84 Assertions.assertEquals(0.094, love.getReal(3, 3), 1.0e-10);
85
86 Assertions.assertEquals(-0.00000, love.getImaginary(2, 0), 1.0e-10);
87 Assertions.assertEquals(-0.00144, love.getImaginary(2, 1), 1.0e-10);
88 Assertions.assertEquals(-0.00130, love.getImaginary(2, 2), 1.0e-10);
89 Assertions.assertEquals(0.0, love.getImaginary(3, 0), 1.0e-10);
90 Assertions.assertEquals(0.0, love.getImaginary(3, 1), 1.0e-10);
91 Assertions.assertEquals(0.0, love.getImaginary(3, 2), 1.0e-10);
92 Assertions.assertEquals(0.0, love.getImaginary(3, 3), 1.0e-10);
93
94 Assertions.assertEquals(-0.00089, love.getPlus(2, 0), 1.0e-10);
95 Assertions.assertEquals(-0.00080, love.getPlus(2, 1), 1.0e-10);
96 Assertions.assertEquals(-0.00057, love.getPlus(2, 2), 1.0e-10);
97 Assertions.assertEquals(0.0, love.getPlus(3, 0), 1.0e-10);
98 Assertions.assertEquals(0.0, love.getPlus(3, 1), 1.0e-10);
99 Assertions.assertEquals(0.0, love.getPlus(3, 2), 1.0e-10);
100 Assertions.assertEquals(0.0, love.getPlus(3, 3), 1.0e-10);
101
102 }
103
104 @Test
105 void testConventions2010() throws NoSuchFieldException, IllegalAccessException {
106
107 UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
108 SolidTidesField tidesField =
109 new SolidTidesField(IERSConventions.IERS_2010.getLoveNumbers(),
110 IERSConventions.IERS_2010.getTideFrequencyDependenceFunction(ut1),
111 IERSConventions.IERS_2010.getPermanentTide(),
112 IERSConventions.IERS_2010.getSolidPoleTide(ut1.getEOPHistory()),
113 FramesFactory.getITRF(IERSConventions.IERS_2010, false),
114 Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_MU,
115 TideSystem.ZERO_TIDE, CelestialBodyFactory.getSun(), CelestialBodyFactory.getMoon());
116
117 Field fieldReal = tidesField.getClass().getDeclaredField("love");
118 fieldReal.setAccessible(true);
119 LoveNumbers love = (LoveNumbers) fieldReal.get(tidesField);
120
121 Assertions.assertEquals(-0.00000, love.getImaginary(2, 0), 1.0e-10);
122 Assertions.assertEquals(-0.00144, love.getImaginary(2, 1), 1.0e-10);
123 Assertions.assertEquals(-0.00130, love.getImaginary(2, 2), 1.0e-10);
124 Assertions.assertEquals(0.0, love.getImaginary(3, 0), 1.0e-10);
125 Assertions.assertEquals(0.0, love.getImaginary(3, 1), 1.0e-10);
126 Assertions.assertEquals(0.0, love.getImaginary(3, 2), 1.0e-10);
127 Assertions.assertEquals(0.0, love.getImaginary(3, 3), 1.0e-10);
128
129
130 Assertions.assertEquals(-0.00089, love.getPlus(2, 0), 1.0e-10);
131 Assertions.assertEquals(-0.00080, love.getPlus(2, 1), 1.0e-10);
132 Assertions.assertEquals(-0.00057, love.getPlus(2, 2), 1.0e-10);
133 Assertions.assertEquals(0.0, love.getPlus(3, 0), 1.0e-10);
134 Assertions.assertEquals(0.0, love.getPlus(3, 1), 1.0e-10);
135 Assertions.assertEquals(0.0, love.getPlus(3, 2), 1.0e-10);
136 Assertions.assertEquals(0.0, love.getPlus(3, 3), 1.0e-10);
137
138 }
139
140 @Test
141 void testK1Example()
142 throws NoSuchFieldException, IllegalAccessException,
143 NoSuchMethodException, InvocationTargetException {
144
145 final PoissonSeriesParser k21Parser =
146 new PoissonSeriesParser(18).
147 withOptionalColumn(1).
148 withDoodson(4, 3).
149 withFirstDelaunay(10);
150 final String name = "/tides/tab6.5a-only-K1.txt";
151 final double pico = 1.0e-12;
152 final PoissonSeries c21Series =
153 k21Parser.
154 withSinCos(0, 17, pico, 18, pico).
155 parse(getClass().getResourceAsStream(name), name);
156 final PoissonSeries s21Series =
157 k21Parser.
158 withSinCos(0, 18, -pico, 17, pico).
159 parse(getClass().getResourceAsStream(name), name);
160 final UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, false);
161 final TimeScalarFunction gmstFunction = IERSConventions.IERS_2010.getGMSTFunction(ut1);
162 Method getNA = IERSConventions.class.getDeclaredMethod("getNutationArguments", TimeScale.class);
163 getNA.setAccessible(true);
164 final FundamentalNutationArguments arguments =
165 (FundamentalNutationArguments) getNA.invoke(IERSConventions.IERS_2010, ut1);
166 TimeVectorFunction deltaCSFunction = new TimeVectorFunction() {
167 public double[] value(final AbsoluteDate date) {
168 final BodiesElements elements = arguments.evaluateAll(date);
169 return new double[] {
170 0.0, c21Series.value(elements), s21Series.value(elements), 0.0, 0.0
171 };
172 }
173 public <T extends CalculusFieldElement<T>> T[] value(final FieldAbsoluteDate<T> date) {
174
175 throw new OrekitInternalError(null);
176 }
177 };
178
179 SolidTidesField tf = new SolidTidesField(IERSConventions.IERS_2010.getLoveNumbers(),
180 deltaCSFunction,
181 IERSConventions.IERS_2010.getPermanentTide(),
182 IERSConventions.IERS_2010.getSolidPoleTide(ut1.getEOPHistory()),
183 FramesFactory.getITRF(IERSConventions.IERS_2010, false),
184 Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
185 Constants.EIGEN5C_EARTH_MU,
186 TideSystem.ZERO_TIDE,
187 CelestialBodyFactory.getSun(),
188 CelestialBodyFactory.getMoon());
189 Method frequencyDependentPart = SolidTidesField.class.getDeclaredMethod("frequencyDependentPart", AbsoluteDate.class, double[][].class, double[][].class);
190 frequencyDependentPart.setAccessible(true);
191 double[][] cachedCNM = new double[5][5];
192 double[][] cachedSNM = new double[5][5];
193
194 AbsoluteDate t0 = new AbsoluteDate(2003, 5, 6, 13, 43, 32.125, TimeScalesFactory.getUTC());
195 for (double dt = 0; dt < Constants.JULIAN_DAY; dt += 300) {
196 AbsoluteDate date = t0.shiftedBy(dt);
197 for (int i = 0; i < cachedCNM.length; ++i) {
198 Arrays.fill(cachedCNM[i], 0.0);
199 Arrays.fill(cachedSNM[i], 0.0);
200 }
201 frequencyDependentPart.invoke(tf, date, cachedCNM, cachedSNM);
202 double thetaPlusPi = gmstFunction.value(date) + FastMath.PI;
203 Assertions.assertEquals(470.9e-12 * FastMath.sin(thetaPlusPi) - 30.2e-12 * FastMath.cos(thetaPlusPi),
204 cachedCNM[2][1], 2.0e-25);
205 Assertions.assertEquals(470.9e-12 * FastMath.cos(thetaPlusPi) + 30.2e-12 * FastMath.sin(thetaPlusPi),
206 cachedSNM[2][1], 2.0e-25);
207 }
208
209 }
210
211 @Test
212 void testDeltaCnmSnm() {
213 NormalizedSphericalHarmonicsProvider gravityField =
214 GravityFieldFactory.getNormalizedProvider(8, 8);
215 UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
216 TimeScale utc = TimeScalesFactory.getUTC();
217
218 AbsoluteDate date = new AbsoluteDate(2003, 5, 6, 13, 43, 32.125, utc);
219 SolidTidesField tidesField =
220 new SolidTidesField(IERSConventions.IERS_2010.getLoveNumbers(),
221 IERSConventions.IERS_2010.getTideFrequencyDependenceFunction(ut1),
222 IERSConventions.IERS_2010.getPermanentTide(),
223 null,
224 FramesFactory.getITRF(IERSConventions.IERS_2010, true),
225 gravityField.getAe(), gravityField.getMu(), TideSystem.TIDE_FREE,
226 CelestialBodyFactory.getSun(), CelestialBodyFactory.getMoon());
227 NormalizedSphericalHarmonics harmonics = tidesField.onDate(date);
228 double[][] refDeltaCnm = new double[][] {
229 { 0.0, 0.0, 0.0, 0.0, 0.0 },
230 { 0.0, 0.0, 0.0, 0.0, 0.0 },
231 { -2.6732289327355114E-9, 4.9078992447259636E-9, 3.5894110538262888E-9, 0.0 , 0.0 },
232
233
234 { -1.290639603871307E-11, -9.287425756410472E-14, 8.356574033404024E-12, -2.2644465207860626E-12, 0.0 },
235 { 7.888138856951149E-12, -1.4422209452877158E-11, -6.815519349970944E-12, 0.0, 0.0 }
236 };
237 double[][] refDeltaSnm = new double[][] {
238 { 0.0, 0.0, 0.0, 0.0, 0.0 },
239 { 0.0, 0.0, 0.0, 0.0, 0.0 },
240 { 0.0, 1.599927449004677E-9, 2.1815888169727694E-9, 0.0 , 0.0 },
241 { 0.0, -4.6129961143785774E-14, 1.8097527720906976E-11, 1.633889224766215E-11, 0.0 },
242 { 0.0, -4.897228975221076E-12, -4.1034042689652575E-12, 0.0, 0.0 }
243 };
244 for (int n = 0; n < refDeltaCnm.length; ++n) {
245 double threshold = (n == 2) ? 1.3e-17 : 2.3e-20;
246 for (int m = 0; m <= n; ++m) {
247 final String reason = "n" + n + "m" + m;
248 MatcherAssert.assertThat(
249 reason,
250 harmonics.getNormalizedCnm(n, m),
251 Matchers.closeTo(refDeltaCnm[n][m], threshold));
252 MatcherAssert.assertThat(
253 reason,
254 harmonics.getNormalizedSnm(n, m),
255 Matchers.closeTo(refDeltaSnm[n][m], threshold));
256 }
257 }
258 }
259
260 @Test
261 void testInterpolationAccuracy() {
262
263
264
265
266
267
268
269
270 final IERSConventions conventions = IERSConventions.IERS_2010;
271 Frame itrf = FramesFactory.getITRF(conventions, true);
272 TimeScale utc = TimeScalesFactory.getUTC();
273 UT1Scale ut1 = TimeScalesFactory.getUT1(conventions, true);
274 NormalizedSphericalHarmonicsProvider gravityField =
275 GravityFieldFactory.getNormalizedProvider(5, 5);
276
277 SolidTidesField raw = new SolidTidesField(conventions.getLoveNumbers(),
278 conventions.getTideFrequencyDependenceFunction(ut1),
279 conventions.getPermanentTide(),
280 conventions.getSolidPoleTide(ut1.getEOPHistory()),
281 itrf, gravityField.getAe(), gravityField.getMu(),
282 gravityField.getTideSystem(),
283 CelestialBodyFactory.getSun(),
284 CelestialBodyFactory.getMoon());
285 int step = 600;
286 int nbPoints = 12;
287 CachedNormalizedSphericalHarmonicsProvider interpolated =
288 new CachedNormalizedSphericalHarmonicsProvider(raw, step, nbPoints,
289 OrekitConfiguration.getCacheSlotsNumber(),
290 7 * Constants.JULIAN_DAY,
291 0.5 * Constants.JULIAN_DAY);
292
293
294 AbsoluteDate start = new AbsoluteDate(2003, 6, 12, utc);
295 AbsoluteDate end = start.shiftedBy(3 * Constants.JULIAN_DAY);
296 StreamingStatistics stat = new StreamingStatistics();
297 for (AbsoluteDate date = start; date.compareTo(end) < 0; date = date.shiftedBy(60)) {
298 NormalizedSphericalHarmonics rawHarmonics = raw.onDate(date);
299 NormalizedSphericalHarmonics interpolatedHarmonics = interpolated.onDate(date);
300
301 for (int n = 2; n < 5; ++n) {
302 for (int m = 0; m <= n; ++m) {
303
304 if (n < 4 || m < 3) {
305 double cnmRaw = rawHarmonics.getNormalizedCnm(n, m);
306 double cnmInterp = interpolatedHarmonics.getNormalizedCnm(n, m);
307 double errorC = (cnmInterp - cnmRaw) / FastMath.abs(cnmRaw);
308 stat.addValue(errorC);
309
310 if (m > 0) {
311 double snmRaw = rawHarmonics.getNormalizedSnm(n, m);
312 double snmInterp = interpolatedHarmonics.getNormalizedSnm(n, m);
313 double errorS = (snmInterp - snmRaw) / FastMath.abs(snmRaw);
314 stat.addValue(errorS);
315 }
316 }
317 }
318 }
319 }
320 Assertions.assertEquals(0.0, stat.getMean(), 2.7e-12);
321 Assertions.assertTrue(stat.getStandardDeviation() < 2.0e-9);
322 Assertions.assertTrue(stat.getMin() > -9.0e-8);
323 Assertions.assertTrue(stat.getMax() < 2.2e-7);
324
325 }
326
327 @BeforeEach
328 public void setUp() {
329 Utils.setDataRoot("regular-data:potential/icgem-format");
330 }
331
332 }