1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.models.earth;
18
19 import java.io.BufferedReader;
20 import java.io.File;
21 import java.io.FileInputStream;
22 import java.io.FileNotFoundException;
23 import java.io.InputStream;
24 import java.io.InputStreamReader;
25 import java.util.Collection;
26 import java.util.StringTokenizer;
27
28 import org.hipparchus.geometry.euclidean.threed.Vector3D;
29 import org.hipparchus.util.FastMath;
30 import org.junit.Assert;
31 import org.junit.BeforeClass;
32 import org.junit.Test;
33 import org.orekit.Utils;
34 import org.orekit.data.DataProvidersManager;
35 import org.orekit.errors.OrekitException;
36 import org.orekit.forces.gravity.potential.EGMFormatReader;
37 import org.orekit.forces.gravity.potential.GravityFieldFactory;
38 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
39 import org.orekit.frames.FramesFactory;
40 import org.orekit.models.earth.GeoMagneticFieldFactory.FieldModel;
41 import org.orekit.time.AbsoluteDate;
42 import org.orekit.time.TimeScalesFactory;
43
44 public class GeoMagneticFieldTest {
45
46
47 private static final int maxOrder = 360;
48
49 private static final int maxDegree = 360;
50
51 private static ReferenceEllipsoid WGS84 = new ReferenceEllipsoid(
52 6378137.00, 1 / 298.257223563, FramesFactory.getGCRF(),
53 3.986004418e14, 7292115e-11);
54
55
56
57
58 private static NormalizedSphericalHarmonicsProvider potential;
59
60
61 private final double[][] wmmTestValues = {
62
63
64 {2015.0, 0, FastMath.toRadians(80), FastMath.toRadians(0), 6627.1, -445.9, 54432.3, 6642.1, 54836.0, FastMath.toRadians(83.04), FastMath.toRadians(-3.85)},
65 {2015.0, 0, FastMath.toRadians(0), FastMath.toRadians(120), 39518.2, 392.9, -11252.4, 39520.2, 41090.9, FastMath.toRadians(-15.89), FastMath.toRadians(0.57)},
66 {2015.0, 0, FastMath.toRadians(-80), FastMath.toRadians(240), 5797.3, 15761.1, -52919.1, 16793.5, 55519.8, FastMath.toRadians(-72.39), FastMath.toRadians(69.81)},
67 {2015.0, 100000, FastMath.toRadians(80), FastMath.toRadians(0), 6314.3, -471.6, 52269.8, 6331.9, 52652.0, FastMath.toRadians(83.09), FastMath.toRadians(-4.27)},
68 {2015.0, 100000, FastMath.toRadians(0), FastMath.toRadians(120), 37535.6, 364.4, -10773.4, 37537.3, 39052.7, FastMath.toRadians(-16.01), FastMath.toRadians(0.56)},
69 {2015.0, 100000, FastMath.toRadians(-80), FastMath.toRadians(240), 5613.1, 14791.5, -50378.6, 15820.7, 52804.4, FastMath.toRadians(-72.57), FastMath.toRadians(69.22)},
70 {2017.5, 0, FastMath.toRadians(80), FastMath.toRadians(0), 6599.4, -317.1, 54459.2, 6607.0, 54858.5, FastMath.toRadians(83.08), FastMath.toRadians(-2.75)},
71 {2017.5, 0, FastMath.toRadians(0), FastMath.toRadians(120), 39571.4, 222.5, -11030.1, 39572.0, 41080.5, FastMath.toRadians(-15.57), FastMath.toRadians(0.32)},
72 {2017.5, 0, FastMath.toRadians(-80), FastMath.toRadians(240), 5873.8, 15781.4, -52687.9, 16839.1, 55313.4, FastMath.toRadians(-72.28), FastMath.toRadians(69.58)},
73 {2017.5, 100000, FastMath.toRadians(80), FastMath.toRadians(0), 6290.5, -348.5, 52292.7, 6300.1, 52670.9, FastMath.toRadians( 83.13), FastMath.toRadians(-3.17)},
74 {2017.5, 100000, FastMath.toRadians(0), FastMath.toRadians(120), 37585.5, 209.5, -10564.2, 37586.1, 39042.5, FastMath.toRadians(-15.70), FastMath.toRadians(0.32)},
75 {2017.5, 100000, FastMath.toRadians(-80), FastMath.toRadians(240), 5683.5, 14808.8, -50163.0, 15862.0, 52611.1, FastMath.toRadians(-72.45), FastMath.toRadians(69.00)}
76 };
77
78
79
80
81 private final double[][] igrfTestValues = {
82
83
84 {2015.0, 0, FastMath.toRadians(80), FastMath.toRadians(0), 6630.9, -447.2, 54434.5, 6645.9, 54838.7, FastMath.toRadians(83.039), FastMath.toRadians(-3.858)},
85 {2015.0, 0, FastMath.toRadians(0), FastMath.toRadians(120), 39519.3, 388.6, -11251.7, 39521.3, 41091.7, FastMath.toRadians(-15.891), FastMath.toRadians(0.563)},
86 {2015.0, 0, FastMath.toRadians(-80), FastMath.toRadians(240), 5808.8, 15754.8, -52945.5, 16791.5, 55544.4, FastMath.toRadians(-72.403), FastMath.toRadians(69.761)},
87 {2015.0, 100000, FastMath.toRadians(80), FastMath.toRadians(0), 6317.2, -472.6, 52272.0, 6334.9, 52654.5, FastMath.toRadians(83.090), FastMath.toRadians(-4.278)},
88 {2015.0, 100000, FastMath.toRadians(0), FastMath.toRadians(120), 37536.9, 361.2, -10773.1, 37538.6, 39053.9, FastMath.toRadians(-16.012), FastMath.toRadians(0.551)},
89 {2015.0, 100000, FastMath.toRadians(-80), FastMath.toRadians(240), 5622.8, 14786.8, -50401.4, 15819.8, 52825.8, FastMath.toRadians(-72.574), FastMath.toRadians(69.180)},
90 {2017.5, 0, FastMath.toRadians(80), FastMath.toRadians(0), 6601.0, -316.4, 54455.5, 6608.5, 54855.0, FastMath.toRadians(83.080), FastMath.toRadians(-2.744)},
91 {2017.5, 0, FastMath.toRadians(0), FastMath.toRadians(120), 39568.1, 225.0, -11041.4, 39568.7, 41080.3, FastMath.toRadians(-15.591), FastMath.toRadians(0.325)},
92 {2017.5, 0, FastMath.toRadians(-80), FastMath.toRadians(240), 5894.7, 15768.1, -52696.8, 16833.9, 55320.2, FastMath.toRadians(-72.283), FastMath.toRadians(69.502)},
93 {2017.5, 100000, FastMath.toRadians(80), FastMath.toRadians(0), 6291.6, -347.2, 52289.9, 6301.2, 52668.2, FastMath.toRadians(83.128), FastMath.toRadians(-3.158)},
94 {2017.5, 100000, FastMath.toRadians(0), FastMath.toRadians(120), 37583.0, 212.3, -10575.1, 37583.6, 39043.0, FastMath.toRadians(-15.715), FastMath.toRadians(0.323)},
95 {2017.5, 100000, FastMath.toRadians(-80), FastMath.toRadians(240), 5702.0, 14797.8, -50170.0, 15858.3, 52616.7, FastMath.toRadians(-72.458), FastMath.toRadians(68.927)}
96 };
97
98
99
100
101
102
103 @BeforeClass
104 public static void setUpBefore() throws Exception {
105 Utils.setDataRoot("earth:geoid:regular-data");
106 GravityFieldFactory.clearPotentialCoefficientsReaders();
107 GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("egm96", false));
108 potential = GravityFieldFactory.getConstantNormalizedProvider(maxDegree, maxOrder);
109 }
110
111 @Test
112 public void testInterpolationYYY5() {
113 double decimalYear = GeoMagneticField.getDecimalYear(1, 1, 2005);
114 GeoMagneticField field = GeoMagneticFieldFactory.getIGRF(decimalYear);
115 GeoMagneticElements e = field.calculateField(FastMath.toRadians(1.2), FastMath.toRadians(0.7), -2500);
116 Assert.assertEquals(FastMath.toRadians(-6.0032), e.getDeclination(), 1.0e-4);
117
118 decimalYear = GeoMagneticField.getDecimalYear(2, 1, 2005);
119 field = GeoMagneticFieldFactory.getIGRF(decimalYear);
120 e = field.calculateField(FastMath.toRadians(1.2), FastMath.toRadians(0.7), -2500);
121 Assert.assertEquals(FastMath.toRadians(-6.0029), e.getDeclination(), 1.0e-4);
122 }
123
124 @Test
125 public void testInterpolationAtEndOfEpoch() {
126 double decimalYear = GeoMagneticField.getDecimalYear(31, 12, 2009);
127 GeoMagneticField field1 = GeoMagneticFieldFactory.getIGRF(decimalYear);
128 GeoMagneticField field2 = GeoMagneticFieldFactory.getIGRF(2010.0);
129
130 Assert.assertNotEquals(field1.getEpoch(), field2.getEpoch());
131
132 GeoMagneticElements e1 = field1.calculateField(0, 0, 0);
133 Assert.assertEquals(FastMath.toRadians(-6.1068), e1.getDeclination(), 1.0e-4);
134
135 GeoMagneticElements e2 = field2.calculateField(0, 0, 0);
136 Assert.assertEquals(FastMath.toRadians(-6.1064), e2.getDeclination(), 1.0e-4);
137 }
138
139 @Test
140 public void testInterpolationAtEndOfValidity() {
141 double decimalYear = GeoMagneticField.getDecimalYear(1, 1, 2020);
142 GeoMagneticField field = GeoMagneticFieldFactory.getIGRF(decimalYear);
143
144 GeoMagneticElements e = field.calculateField(0, 0, 0);
145 Assert.assertEquals(FastMath.toRadians(-4.7446), e.getDeclination(), 1.0e-4);
146 }
147
148 @Test
149 public void testContinuityAtPole() {
150 double decimalYear = GeoMagneticField.getDecimalYear(1, 1, 2020);
151 GeoMagneticField field = GeoMagneticFieldFactory.getIGRF(decimalYear);
152
153 GeoMagneticElements eClose = field.calculateField(FastMath.toRadians(89.999999), 0, 0);
154 GeoMagneticElements ePole = field.calculateField(FastMath.toRadians(90.0), 0, 0);
155 Assert.assertEquals("" + (eClose.getDeclination()- ePole.getDeclination()), eClose.getDeclination(), ePole.getDeclination(), 7.0e-7);
156 Assert.assertEquals("" + (eClose.getInclination()- ePole.getInclination()), eClose.getInclination(), ePole.getInclination(), 3.0e-7);
157 Assert.assertEquals("" + (eClose.getTotalIntensity()- ePole.getTotalIntensity()), eClose.getTotalIntensity(), ePole.getTotalIntensity(), 2.0e-4);
158 Assert.assertEquals("" + (eClose.getHorizontalIntensity()- ePole.getHorizontalIntensity()), eClose.getHorizontalIntensity(), ePole.getHorizontalIntensity(), 3.0e-4);
159 }
160
161 @Test(expected=OrekitException.class)
162 public void testTransformationOutsideValidityPeriod() {
163 double decimalYear = GeoMagneticField.getDecimalYear(10, 1, 2020);
164 @SuppressWarnings("unused")
165 GeoMagneticField field = GeoMagneticFieldFactory.getIGRF(decimalYear);
166 }
167
168 @Test
169 public void testWMM() throws Exception {
170
171
172
173
174
175
176 runSampleFile(FieldModel.WMM, "sample_coords.txt", "sample_out_WMM2015.txt");
177
178 final double eps = 1e-1;
179 final double degreeEps = 1e-2;
180 for (int i = 0; i < wmmTestValues.length; i++) {
181 final GeoMagneticField model = GeoMagneticFieldFactory.getWMM(wmmTestValues[i][0]);
182 final GeoMagneticElements result = model.calculateField(wmmTestValues[i][2],
183 wmmTestValues[i][3],
184 wmmTestValues[i][1]);
185
186
187 Assert.assertEquals(wmmTestValues[i][4], result.getFieldVector().getX(), eps);
188
189 Assert.assertEquals(wmmTestValues[i][5], result.getFieldVector().getY(), eps);
190
191 Assert.assertEquals(wmmTestValues[i][6], result.getFieldVector().getZ(), eps);
192
193 Assert.assertEquals(wmmTestValues[i][7], result.getHorizontalIntensity(), eps);
194
195 Assert.assertEquals(wmmTestValues[i][8], result.getTotalIntensity(), eps);
196
197 Assert.assertEquals(wmmTestValues[i][9], result.getInclination(), degreeEps);
198
199 Assert.assertEquals(wmmTestValues[i][10], result.getDeclination(), degreeEps);
200 }
201 }
202
203 @Test
204 public void testWMMWithHeightAboveMSL() throws Exception {
205
206
207
208
209 final double[][] testValues = {
210
211
212 {2015.0, 100000, FastMath.toRadians(80), FastMath.toRadians(0), 6314.2, -471.6, 52269.1, 6331.8, 52651.2, FastMath.toRadians(83.093), FastMath.toRadians(-4.271)},
213 {2015.0, 100000, FastMath.toRadians(0), FastMath.toRadians(120), 37534.4, 364.3, -10773.1, 37536.2, 39051.6, FastMath.toRadians(-16.013), FastMath.toRadians(0.556)},
214 {2015.0, 100000, FastMath.toRadians(-80), FastMath.toRadians(240), 5613.2, 14791.9, -50379.6, 15821.1, 52805.4, FastMath.toRadians(-72.565), FastMath.toRadians(69.219)}
215 };
216
217 final Geoid geoid = new Geoid(potential, WGS84);
218
219 final double eps = 1e-1;
220 final double degreeEps = 1e-2;
221 for (int i = 0; i < testValues.length; i++) {
222 final AbsoluteDate date = new AbsoluteDate(2015, 1, 1, TimeScalesFactory.getUTC());
223 final GeoMagneticField model = GeoMagneticFieldFactory.getWMM(testValues[i][0]);
224 final double undulation = geoid.getUndulation(testValues[i][2],
225 testValues[i][3],
226 date);
227 final GeoMagneticElements result = model.calculateField(testValues[i][2],
228 testValues[i][3],
229 testValues[i][1] + undulation);
230
231
232 Assert.assertEquals(testValues[i][4], result.getFieldVector().getX(), eps);
233
234 Assert.assertEquals(testValues[i][5], result.getFieldVector().getY(), eps);
235
236 Assert.assertEquals(testValues[i][6], result.getFieldVector().getZ(), eps);
237
238 Assert.assertEquals(testValues[i][7], result.getHorizontalIntensity(), eps);
239
240 Assert.assertEquals(testValues[i][8], result.getTotalIntensity(), eps);
241
242 Assert.assertEquals(testValues[i][9], result.getInclination(), degreeEps);
243
244 Assert.assertEquals(testValues[i][10], result.getDeclination(), degreeEps);
245 }
246 }
247
248 @Test
249 public void testIGRF() throws Exception {
250
251
252
253
254
255
256 runSampleFile(FieldModel.IGRF, "sample_coords.txt", "sample_out_IGRF12.txt");
257
258 final double eps = 1e-1;
259 final double degreeEps = 1e-2;
260 for (int i = 0; i < igrfTestValues.length; i++) {
261 final GeoMagneticField model = GeoMagneticFieldFactory.getIGRF(igrfTestValues[i][0]);
262 final GeoMagneticElements result = model.calculateField(igrfTestValues[i][2],
263 igrfTestValues[i][3],
264 igrfTestValues[i][1]);
265
266 final Vector3D b = result.getFieldVector();
267
268
269 Assert.assertEquals(igrfTestValues[i][4], b.getX(), eps);
270
271 Assert.assertEquals(igrfTestValues[i][5], b.getY(), eps);
272
273 Assert.assertEquals(igrfTestValues[i][6], b.getZ(), eps);
274
275 Assert.assertEquals(igrfTestValues[i][7], result.getHorizontalIntensity(), eps);
276
277 Assert.assertEquals(igrfTestValues[i][8], result.getTotalIntensity(), eps);
278
279 Assert.assertEquals(igrfTestValues[i][9], result.getInclination(), degreeEps);
280
281 Assert.assertEquals(igrfTestValues[i][10], result.getDeclination(), degreeEps);
282 }
283 }
284
285 @Test(expected=OrekitException.class)
286 public void testUnsupportedTransform() throws Exception {
287 final GeoMagneticField model = GeoMagneticFieldFactory.getIGRF(1910);
288
289
290 model.transformModel(1950);
291 }
292
293 @Test(expected=OrekitException.class)
294 public void testOutsideValidityTransform() throws Exception {
295 final GeoMagneticField model1 = GeoMagneticFieldFactory.getIGRF(2005);
296 final GeoMagneticField model2 = GeoMagneticFieldFactory.getIGRF(2010);
297
298
299 model1.transformModel(model2, 2012);
300 }
301
302 @Test
303 public void testValidTransform() throws Exception {
304 final GeoMagneticField model = GeoMagneticFieldFactory.getWMM(2015);
305
306 Assert.assertTrue(model.supportsTimeTransform());
307
308 final GeoMagneticField transformedModel = model.transformModel(2017);
309
310 Assert.assertEquals(2015, transformedModel.validFrom(), 1e0);
311 Assert.assertEquals(2020, transformedModel.validTo(), 1e0);
312 Assert.assertEquals(2017, transformedModel.getEpoch(), 1e0);
313 }
314
315 @Test
316 public void testLoadOriginalWMMModel() throws Exception {
317 GeoMagneticModelLoader loader = new GeoMagneticModelLoader();
318
319 InputStream input = getResource("WMM2015.COF");
320 loader.loadData(input, "WMM2015.COF");
321
322 Collection<GeoMagneticField> models = loader.getModels();
323 Assert.assertNotNull(models);
324 Assert.assertEquals(1, models.size());
325
326 GeoMagneticField wmmModel = models.iterator().next();
327 Assert.assertEquals("WMM-2015", wmmModel.getModelName());
328 Assert.assertEquals(2015, wmmModel.getEpoch(), 1e-9);
329
330 final double eps = 1e-1;
331 final double degreeEps = 1e-2;
332 for (int i = 0; i < wmmTestValues.length; i++) {
333 if (wmmTestValues[i][0] != wmmModel.getEpoch()) {
334 continue;
335 }
336 final GeoMagneticElements result = wmmModel.calculateField(wmmTestValues[i][2],
337 wmmTestValues[i][3],
338 wmmTestValues[i][1]);
339
340
341 Assert.assertEquals(wmmTestValues[i][4], result.getFieldVector().getX(), eps);
342
343 Assert.assertEquals(wmmTestValues[i][5], result.getFieldVector().getY(), eps);
344
345 Assert.assertEquals(wmmTestValues[i][6], result.getFieldVector().getZ(), eps);
346
347 Assert.assertEquals(wmmTestValues[i][7], result.getHorizontalIntensity(), eps);
348
349 Assert.assertEquals(wmmTestValues[i][8], result.getTotalIntensity(), eps);
350
351 Assert.assertEquals(wmmTestValues[i][9], result.getInclination(), degreeEps);
352
353 Assert.assertEquals(wmmTestValues[i][10], result.getDeclination(), degreeEps);
354 }
355 }
356
357 public void runSampleFile(final FieldModel type, final String inputFile, final String outputFile)
358 throws Exception {
359
360 final BufferedReader inReader = new BufferedReader(new InputStreamReader(getResource(inputFile)));
361 final BufferedReader outReader = new BufferedReader(new InputStreamReader(getResource(outputFile)));
362
363
364 outReader.readLine();
365
366 String line = null;
367 while ((line = inReader.readLine()) != null) {
368 if (line.trim().length() == 0) {
369 break;
370 }
371
372 final StringTokenizer st = new StringTokenizer(line);
373
374 final double year = getYear(st.nextToken());
375 final String heightType = st.nextToken();
376 final String heightStr = st.nextToken();
377 final double lat = getLatLon(st.nextToken());
378 final double lon = getLatLon(st.nextToken());
379
380 final GeoMagneticField field = GeoMagneticFieldFactory.getField(type, year);
381
382 double height = Double.valueOf(heightStr.substring(1));
383 if (heightStr.startsWith("K")) {
384
385 height *= 1000d;
386 } else if (heightStr.startsWith("F")) {
387
388 height *= 0.3048;
389 }
390
391
392 final GeoMagneticElements ge = field.calculateField(lat, lon, height);
393 final String validateLine = outReader.readLine();
394
395 if (!heightType.startsWith("C")) {
396 validate(ge, validateLine);
397 }
398
399 String geString = ge.toString();
400 Assert.assertNotNull(geString);
401 Assert.assertFalse(geString.isEmpty());
402 }
403 }
404
405 private double getYear(final String yearStr) {
406
407 if (yearStr.contains(",")) {
408 final StringTokenizer st = new StringTokenizer(yearStr, ",");
409 final String y = st.nextToken();
410 final String m = st.nextToken();
411 final String d = st.nextToken();
412 return GeoMagneticField.getDecimalYear(Integer.valueOf(d), Integer.valueOf(m), Integer.valueOf(y));
413 } else {
414 return Double.valueOf(yearStr);
415 }
416 }
417
418 private double getLatLon(final String str) {
419
420 if (str.contains(",")) {
421 final StringTokenizer st = new StringTokenizer(str, ",");
422 final int d = Integer.valueOf(st.nextToken());
423 final int m = Integer.valueOf(st.nextToken());
424 int s = 0;
425 if (st.hasMoreTokens()) {
426 s = Integer.valueOf(st.nextToken());
427 }
428 double deg = FastMath.abs(d) + m / 60d + s / 3600d;
429 if (d < 0) {
430 deg = -deg;
431 }
432 return FastMath.toRadians(deg);
433 } else {
434 return FastMath.toRadians(Double.valueOf(str));
435 }
436 }
437
438 private double getRadians(final String degree, final String minute) {
439 double result = Double.valueOf(degree.substring(0, degree.length() - 1));
440 final double min = Double.valueOf(minute.substring(0, minute.length() - 1)) / 60d;
441 result += (result < 0) ? -min : min;
442 return FastMath.toRadians(result);
443 }
444
445 @SuppressWarnings("unused")
446 private void validate(final GeoMagneticElements ge, final String outputLine)
447 throws Exception {
448
449 final StringTokenizer st = new StringTokenizer(outputLine);
450
451 final double year = getYear(st.nextToken());
452
453 final String coord = st.nextToken();
454 final String heightStr = st.nextToken();
455 final String latStr = st.nextToken();
456 final String lonStr = st.nextToken();
457
458 final double dec = getRadians(st.nextToken(), st.nextToken());
459 final double inc = getRadians(st.nextToken(), st.nextToken());
460
461 final double h = Double.valueOf(st.nextToken());
462 final double x = Double.valueOf(st.nextToken());
463 final double y = Double.valueOf(st.nextToken());
464 final double z = Double.valueOf(st.nextToken());
465 final double f = Double.valueOf(st.nextToken());
466
467 final double eps = 1e-1;
468 Assert.assertEquals(h, ge.getHorizontalIntensity(), eps);
469 Assert.assertEquals(f, ge.getTotalIntensity(), eps);
470 Assert.assertEquals(x, ge.getFieldVector().getX(), eps);
471 Assert.assertEquals(y, ge.getFieldVector().getY(), eps);
472 Assert.assertEquals(z, ge.getFieldVector().getZ(), eps);
473 Assert.assertEquals(dec, ge.getDeclination(), eps);
474 Assert.assertEquals(inc, ge.getInclination(), eps);
475 }
476
477 private InputStream getResource(final String name) throws FileNotFoundException {
478
479 final String separator = System.getProperty("path.separator");
480 final String dataPath = System.getProperty(DataProvidersManager.OREKIT_DATA_PATH).split(separator)[0];
481 return new FileInputStream(new File(dataPath, name));
482 }
483
484 }