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