1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.models.earth.displacement;
18
19 import org.hipparchus.geometry.euclidean.threed.Vector3D;
20 import org.hipparchus.util.FastMath;
21 import org.junit.jupiter.api.Assertions;
22 import org.junit.jupiter.api.BeforeEach;
23 import org.junit.jupiter.api.Test;
24 import org.orekit.Utils;
25 import org.orekit.bodies.OneAxisEllipsoid;
26 import org.orekit.data.BodiesElements;
27 import org.orekit.data.FundamentalNutationArguments;
28 import org.orekit.frames.FramesFactory;
29 import org.orekit.time.AbsoluteDate;
30 import org.orekit.time.TimeScale;
31 import org.orekit.time.TimeScalesFactory;
32 import org.orekit.utils.Constants;
33 import org.orekit.utils.IERSConventions;
34
35 import java.lang.reflect.Field;
36 import java.util.Collections;
37 import java.util.List;
38 import java.util.Map;
39 import java.util.stream.Collectors;
40
41 public class OceanLoadingTest {
42
43 private OneAxisEllipsoid earth;
44
45 @Test
46 public void testSemiDiurnal() {
47 TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
48 FundamentalNutationArguments fna = IERSConventions.IERS_2010.getNutationArguments(ut1);
49 BodiesElements elements = fna.evaluateAll(new AbsoluteDate(2009, 6, 25, 0, 0, 0.0, ut1));
50 for (Tide tide : getTides()) {
51 if (tide.getDoodsonMultipliers()[0] == 2) {
52 double f = tide.getRate(elements) * Constants.JULIAN_DAY/ (2 * FastMath.PI);
53 Assertions.assertTrue(f > 1.5);
54 Assertions.assertTrue(f <= 2.5);
55 }
56 }
57 }
58
59 @Test
60 public void testDiurnal() {
61 TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
62 FundamentalNutationArguments fna = IERSConventions.IERS_2010.getNutationArguments(ut1);
63 BodiesElements elements = fna.evaluateAll(new AbsoluteDate(2009, 6, 25, 0, 0, 0.0, ut1));
64 for (Tide tide : getTides()) {
65 if (tide.getDoodsonMultipliers()[0] == 1) {
66 double f = tide.getRate(elements) * Constants.JULIAN_DAY/ (2 * FastMath.PI);
67 Assertions.assertTrue(f > 0.5);
68 Assertions.assertTrue(f <= 1.5);
69 }
70 }
71 }
72
73 @Test
74 public void testLongPeriod() {
75 TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
76 FundamentalNutationArguments fna = IERSConventions.IERS_2010.getNutationArguments(ut1);
77 BodiesElements elements = fna.evaluateAll(new AbsoluteDate(2009, 6, 25, 0, 0, 0.0, ut1));
78 for (Tide tide : getTides()) {
79 if (tide.getDoodsonMultipliers()[0] == 0) {
80 double f = tide.getRate(elements) * Constants.JULIAN_DAY/ (2 * FastMath.PI);
81 Assertions.assertTrue(f > 0.0);
82 Assertions.assertTrue(f <= 0.5);
83 }
84 }
85 }
86
87 @Test
88 public void testNoExtra() {
89 for (Tide tide : getTides()) {
90 if (tide.getDoodsonMultipliers()[0] > 2) {
91 Assertions.fail("unexpected tide " + tide.getDoodsonNumber());
92 }
93 }
94 }
95
96 @Test
97 public void testStableRates() {
98
99
100
101
102
103
104
105 final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
106 final FundamentalNutationArguments fna = IERSConventions.IERS_2010.getNutationArguments(ut1);
107
108
109 final BodiesElements el2000 = fna.evaluateAll(AbsoluteDate.J2000_EPOCH);
110 List<Tide> tides = getTides();
111 Collections.sort(tides, (t1, t2) -> Double.compare(t1.getRate(el2000), t2.getRate(el2000)));
112
113 for (double dt = -122000; dt < 54000; dt += 100) {
114 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(dt * Constants.JULIAN_YEAR);
115 final BodiesElements el = fna.evaluateAll(date);
116 for (int i = 1; i < tides.size(); ++i) {
117 final Tide t1 = tides.get(i - 1);
118 final Tide t2 = tides.get(i);
119 Assertions.assertTrue(t1.getRate(el) < t2.getRate(el));
120 }
121 }
122 }
123
124 @Test
125 public void testTidesRatesPastInversion() {
126
127 final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
128 final FundamentalNutationArguments fna = IERSConventions.IERS_2010.getNutationArguments(ut1);
129 final Tide tide1 = new Tide(245556);
130 final Tide tide2 = new Tide(245635);
131
132 final AbsoluteDate t0530 = new AbsoluteDate(-122502, 11, 9, 5, 30, 0.0, TimeScalesFactory.getTAI());
133 final BodiesElements el0530 = fna.evaluateAll(t0530);
134 Assertions.assertTrue(tide1.getRate(el0530) < tide2.getRate(el0530));
135
136 final AbsoluteDate t0430 = t0530.shiftedBy(-3600.0);
137 final BodiesElements el0430 = fna.evaluateAll(t0430);
138 Assertions.assertTrue(tide1.getRate(el0430) > tide2.getRate(el0430));
139
140 }
141
142 @Test
143 public void testTidesRatesFutureInversion() {
144
145 final TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
146 final FundamentalNutationArguments fna = IERSConventions.IERS_2010.getNutationArguments(ut1);
147 final Tide tide1 = new Tide(274554);
148 final Tide tide2 = new Tide(274556);
149
150 final AbsoluteDate t1700 = new AbsoluteDate(56824, 11, 2, 17, 0, 0.0, TimeScalesFactory.getTAI());
151 final BodiesElements el1700 = fna.evaluateAll(t1700);
152 Assertions.assertTrue(tide1.getRate(el1700) < tide2.getRate(el1700));
153
154 final AbsoluteDate t1800 = t1700.shiftedBy(3600.0);
155 final BodiesElements el1800 = fna.evaluateAll(t1800);
156 Assertions.assertTrue(tide1.getRate(el1800) > tide2.getRate(el1800));
157
158 }
159
160 @Test
161 public void testOnsalaOriginalEarthRotation() {
162
163 doTestOnsala(true, 4.7e-6);
164 }
165
166 @Test
167 public void testOnsalaIERSEarthRotation() {
168
169 doTestOnsala(false, 3.4e-6);
170 }
171
172 private void doTestOnsala(boolean patchEarthRotation, double tolerance) {
173 OceanLoadingCoefficientsBLQFactory factory = new OceanLoadingCoefficientsBLQFactory("^hardisp\\.blq$");
174 OceanLoadingCoefficients coefficients = factory.getCoefficients("Onsala");
175 Vector3D refPoint = earth.transform(coefficients.getSiteLocation());
176 OceanLoading loading = new OceanLoading(earth, coefficients);
177 TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
178 FundamentalNutationArguments fna = IERSConventions.IERS_2010.getNutationArguments(ut1);
179 if (patchEarthRotation) {
180 TideTest.patchEarthRotationModel(fna, ut1);
181 }
182 AbsoluteDate t0 = new AbsoluteDate("2009-06-25T01:10:45", ut1 );
183 double[][] ref = new double[][] {
184 { 0.003094, -0.001538, -0.000895 },
185 { 0.001812, -0.000950, -0.000193 },
186 { 0.000218, -0.000248, 0.000421 },
187 { -0.001104, 0.000404, 0.000741 },
188 { -0.001668, 0.000863, 0.000646 },
189 { -0.001209, 0.001042, 0.000137 },
190 { 0.000235, 0.000926, -0.000667 },
191 { 0.002337, 0.000580, -0.001555 },
192 { 0.004554, 0.000125, -0.002278 },
193 { 0.006271, -0.000291, -0.002615 },
194 { 0.006955, -0.000537, -0.002430 },
195 { 0.006299, -0.000526, -0.001706 },
196 { 0.004305, -0.000244, -0.000559 },
197 { 0.001294, 0.000245, 0.000793 },
198 { -0.002163, 0.000819, 0.002075 },
199 { -0.005375, 0.001326, 0.003024 },
200 { -0.007695, 0.001622, 0.003448 },
201 { -0.008669, 0.001610, 0.003272 },
202 { -0.008143, 0.001262, 0.002557 },
203 { -0.006290, 0.000633, 0.001477 },
204 { -0.003566, -0.000155, 0.000282 },
205 { -0.000593, -0.000941, -0.000766 },
206 { 0.001992, -0.001561, -0.001457 },
207 { 0.003689, -0.001889, -0.001680 }
208 };
209 for (int i = 0; i < ref.length; ++i) {
210 BodiesElements elements = fna.evaluateAll(t0.shiftedBy(i * 3600.0));
211 final Vector3D d = loading.displacement(elements, earth.getBodyFrame(), refPoint);
212 Assertions.assertEquals(ref[i][0], Vector3D.dotProduct(d, coefficients.getSiteLocation().getZenith()), tolerance);
213 Assertions.assertEquals(ref[i][1], Vector3D.dotProduct(d, coefficients.getSiteLocation().getSouth()), tolerance);
214 Assertions.assertEquals(ref[i][2], Vector3D.dotProduct(d, coefficients.getSiteLocation().getWest()), tolerance);
215 }
216 }
217
218 @Test
219 public void testReykjavikOriginalEarthRotation() {
220
221 doTestReykjavik(true, 9.3e-6);
222 }
223
224 @Test
225 public void testReykjavikIERSEarthRotation() {
226
227 doTestReykjavik(false, 16.9e-6);
228 }
229
230 private void doTestReykjavik(boolean patchEarthRotation, double tolerance) {
231
232
233
234
235
236 OceanLoadingCoefficientsBLQFactory factory = new OceanLoadingCoefficientsBLQFactory("^hardisp\\.blq$");
237 OceanLoadingCoefficients coefficients = factory.getCoefficients("Reykjavik");
238 Vector3D refPoint = earth.transform(coefficients.getSiteLocation());
239 OceanLoading loading = new OceanLoading(earth, coefficients);
240 TimeScale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
241 FundamentalNutationArguments fna = IERSConventions.IERS_2010.getNutationArguments(ut1);
242 if (patchEarthRotation) {
243 TideTest.patchEarthRotationModel(fna, ut1);
244 }
245 AbsoluteDate t0 = new AbsoluteDate("2009-06-25T01:10:45", ut1 );
246 double[][] ref = new double[][] {
247 { -0.005940, -0.001245, -0.000278 },
248 { 0.013516, -0.001086, 0.003212 },
249 { 0.029599, -0.000353, 0.005483 },
250 { 0.038468, 0.000699, 0.005997 },
251 { 0.038098, 0.001721, 0.004690 },
252 { 0.028780, 0.002363, 0.001974 },
253 { 0.013016, 0.002371, -0.001369 },
254 { -0.005124, 0.001653, -0.004390 },
255 { -0.021047, 0.000310, -0.006225 },
256 { -0.030799, -0.001383, -0.006313 },
257 { -0.032056, -0.003048, -0.004549 },
258 { -0.024698, -0.004288, -0.001314 },
259 { -0.010814, -0.004794, 0.002623 },
260 { 0.005849, -0.004416, 0.006291 },
261 { 0.020857, -0.003208, 0.008766 },
262 { 0.030226, -0.001413, 0.009402 },
263 { 0.031437, 0.000594, 0.007996 },
264 { 0.024079, 0.002389, 0.004844 },
265 { 0.009945, 0.003606, 0.000663 },
266 { -0.007426, 0.004022, -0.003581 },
267 { -0.023652, 0.003601, -0.006911 },
268 { -0.034618, 0.002505, -0.008585 },
269 { -0.037515, 0.001044, -0.008270 },
270 { -0.031544, -0.000402, -0.006125 }
271 };
272 for (int i = 0; i < ref.length; ++i) {
273 BodiesElements elements = fna.evaluateAll(t0.shiftedBy(i * 3600.0));
274 final Vector3D d = loading.displacement(elements, earth.getBodyFrame(), refPoint);
275 Assertions.assertEquals(ref[i][0], Vector3D.dotProduct(d, coefficients.getSiteLocation().getZenith()), tolerance);
276 Assertions.assertEquals(ref[i][1], Vector3D.dotProduct(d, coefficients.getSiteLocation().getSouth()), tolerance);
277 Assertions.assertEquals(ref[i][2], Vector3D.dotProduct(d, coefficients.getSiteLocation().getWest()), tolerance);
278 }
279 }
280
281 private List<Tide> getTides() {
282 try {
283 Field mapField = OceanLoading.class.getDeclaredField("CARTWRIGHT_EDDEN_AMPLITUDE_MAP");
284 mapField.setAccessible(true);
285 @SuppressWarnings("unchecked")
286 Map<Tide, Double> map = (Map<Tide, Double>) mapField.get(null);
287 return map.entrySet().stream().map(e -> e.getKey()).collect(Collectors.toList());
288 } catch (NoSuchFieldException | IllegalArgumentException | IllegalAccessException e) {
289 Assertions.fail(e.getLocalizedMessage());
290 return null;
291 }
292 }
293
294 @BeforeEach
295 public void setUp() throws Exception {
296 Utils.setDataRoot("regular-data:oso-blq");
297 earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
298 Constants.WGS84_EARTH_FLATTENING,
299 FramesFactory.getITRF(IERSConventions.IERS_2010, false));
300 }
301
302 }