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