1   /* Copyright 2011-2012 Space Applications Services
2    * Licensed to CS Communication & Systèmes (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.models.earth;
18  
19  import java.io.BufferedReader;
20  import java.io.File;
21  import java.io.IOException;
22  import java.io.InputStream;
23  import java.io.InputStreamReader;
24  import java.nio.file.Files;
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.jupiter.api.Assertions;
31  import org.junit.jupiter.api.BeforeAll;
32  import org.junit.jupiter.api.Test;
33  import org.orekit.Utils;
34  import org.orekit.data.DataProvidersManager;
35  import org.orekit.data.DataSource;
36  import org.orekit.errors.OrekitException;
37  import org.orekit.forces.gravity.potential.EGMFormatReader;
38  import org.orekit.forces.gravity.potential.GravityFieldFactory;
39  import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
40  import org.orekit.frames.FramesFactory;
41  import org.orekit.models.earth.GeoMagneticFieldFactory.FieldModel;
42  import org.orekit.time.AbsoluteDate;
43  import org.orekit.time.TimeScalesFactory;
44  import org.orekit.utils.units.UnitsConverter;
45  
46  public class GeoMagneticFieldTest {
47  
48      /** maximum degree used in testing {@link Geoid}. */
49      private static final int maxOrder = 360;
50      /** maximum order used in testing {@link Geoid}. */
51      private static final int maxDegree = 360;
52      /** The WGS84 reference ellipsoid. */
53      private static final ReferenceEllipsoid WGS84 = new ReferenceEllipsoid(
54              6378137.00, 1 / 298.257223563, FramesFactory.getGCRF(),
55              3.986004418e14, 7292115e-11);
56  
57      // test results for test values provided as part of the WMM2015 Report
58      private final double[][] wmmTestValues = {
59          // Date    Alt             Lat                       Lon                X        Y         Z        H        F              I                            D
60          //          m              rad                       rad               nT       nT        nT       nT       nT             rad                          rad
61          {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)},
62          {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)},
63          {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)},
64          {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)},
65          {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)},
66          {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)},
67          {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)},
68          {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)},
69          {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)},
70          {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)},
71          {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)},
72          {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)}
73      };
74  
75      // test results for test values provided as part of the WMM2015 Report
76      // the results for the IGRF12 model have been obtained from the NOAA
77      // online calculator: http://www.ngdc.noaa.gov/geomag-web/#igrfwmm
78      private final double[][] igrfTestValues = {
79          // Date     Alt              Lat                       Lon             X        Y         Z        H        F                 I                           D
80          //            m              rad                       rad            nT       nT        nT       nT       nT                rad                         rad
81          {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)},
82          {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)},
83          {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)},
84          {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)},
85          {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)},
86          {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)},
87          {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)},
88          {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)},
89          {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)},
90          {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)},
91          {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)},
92          {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)}
93      };
94  
95      /**
96       * load orekit data and gravity field.
97       *
98       * @throws Exception on error.
99       */
100     @BeforeAll
101     public static void setUpBefore() throws Exception {
102         Utils.setDataRoot("earth:geoid:regular-data");
103     }
104 
105     @Test
106     public void testInterpolationYYY5() {
107         double decimalYear = GeoMagneticField.getDecimalYear(1, 1, 2005);
108         GeoMagneticField field = GeoMagneticFieldFactory.getIGRF(decimalYear);
109         GeoMagneticElements e = field.calculateField(FastMath.toRadians(1.2), FastMath.toRadians(0.7), -2500);
110         Assertions.assertEquals(FastMath.toRadians(-6.0032), e.getDeclination(), 1.0e-4);
111 
112         decimalYear = GeoMagneticField.getDecimalYear(2, 1, 2005);
113         field = GeoMagneticFieldFactory.getIGRF(decimalYear);
114         e = field.calculateField(FastMath.toRadians(1.2), FastMath.toRadians(0.7), -2500);
115         Assertions.assertEquals(FastMath.toRadians(-6.0029), e.getDeclination(), 1.0e-4);
116     }
117 
118     @Test
119     public void testInterpolationAtEndOfEpoch() {
120         double decimalYear = GeoMagneticField.getDecimalYear(31, 12, 2009);
121         GeoMagneticField field1 = GeoMagneticFieldFactory.getIGRF(decimalYear);
122         GeoMagneticField field2 = GeoMagneticFieldFactory.getIGRF(2010.0);
123 
124         Assertions.assertNotEquals(field1.getEpoch(), field2.getEpoch());
125 
126         GeoMagneticElements e1 = field1.calculateField(0, 0, 0);
127         Assertions.assertEquals(FastMath.toRadians(-6.1068), e1.getDeclination(), 1.0e-4);
128 
129         GeoMagneticElements e2 = field2.calculateField(0, 0, 0);
130         Assertions.assertEquals(FastMath.toRadians(-6.1064), e2.getDeclination(), 1.0e-4);
131     }
132 
133     @Test
134     public void testInterpolationAtEndOfValidity() {
135         double decimalYear = GeoMagneticField.getDecimalYear(1, 1, 2020);
136         GeoMagneticField field = GeoMagneticFieldFactory.getIGRF(decimalYear);
137 
138         GeoMagneticElements e = field.calculateField(0, 0, 0);
139         Assertions.assertEquals(FastMath.toRadians(-4.7446), e.getDeclination(), 1.0e-4);
140     }
141 
142     @Test
143     public void testContinuityAtPole() {
144         double decimalYear = GeoMagneticField.getDecimalYear(1, 1, 2020);
145         GeoMagneticField field = GeoMagneticFieldFactory.getIGRF(decimalYear);
146 
147         GeoMagneticElements eClose = field.calculateField(FastMath.toRadians(89.999999), 0, 0);
148         GeoMagneticElements ePole  = field.calculateField(FastMath.toRadians(90.0),      0, 0);
149         Assertions.assertEquals(eClose.getDeclination(),         ePole.getDeclination(),         7.0e-7, "" + (eClose.getDeclination()-         ePole.getDeclination()));
150         Assertions.assertEquals(eClose.getInclination(),         ePole.getInclination(),         3.0e-7, "" + (eClose.getInclination()-         ePole.getInclination()));
151         Assertions.assertEquals(eClose.getTotalIntensity(),      ePole.getTotalIntensity(),      2.0e-13, "" + (eClose.getTotalIntensity()-      ePole.getTotalIntensity()));
152         Assertions.assertEquals(eClose.getHorizontalIntensity(), ePole.getHorizontalIntensity(), 3.0e-13, "" + (eClose.getHorizontalIntensity()- ePole.getHorizontalIntensity()));
153     }
154 
155     @Test
156     public void testTransformationOutsideValidityPeriod() {
157         Assertions.assertThrows(OrekitException.class, () -> {
158             double decimalYear = GeoMagneticField.getDecimalYear(10, 1, 2020);
159             @SuppressWarnings("unused")
160             GeoMagneticField field = GeoMagneticFieldFactory.getIGRF(decimalYear);
161         });
162     }
163 
164     @Test
165     public void testWMM() throws Exception {
166         // test values from sample coordinate file
167         // provided as part of the geomag 7.0 distribution available at
168         // http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html
169         // modification: the julian day calculation of geomag is slightly different
170         // to the one from the WMM code, we use the WMM convention thus the outputs
171         // have been adapted.
172         runSampleFile(FieldModel.WMM, "sample_coords.txt", "sample_out_WMM2015.txt");
173 
174         final double eps = 1e-10;
175         final double degreeEps = 1e-2;
176         for (final double[] wmmTestValue : wmmTestValues) {
177             final GeoMagneticField model = GeoMagneticFieldFactory.getWMM(wmmTestValue[0]);
178             final GeoMagneticElements result = model.calculateField(wmmTestValue[2], wmmTestValue[3], wmmTestValue[1]);
179 
180             // X
181             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(wmmTestValue[4]),
182                                     result.getFieldVector().getX(), eps);
183             // Y
184             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(wmmTestValue[5]),
185                                     result.getFieldVector().getY(), eps);
186             // Z
187             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(wmmTestValue[6]),
188                                     result.getFieldVector().getZ(), eps);
189             // H
190             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(wmmTestValue[7]),
191                                     result.getHorizontalIntensity(), eps);
192             // F
193             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(wmmTestValue[8]),
194                                     result.getTotalIntensity(), eps);
195             // inclination
196             Assertions.assertEquals(wmmTestValue[9], result.getInclination(), degreeEps);
197             // declination
198             Assertions.assertEquals(wmmTestValue[10], result.getDeclination(), degreeEps);
199         }
200     }
201 
202     @Test
203     public void testWMMWithHeightAboveMSL() {
204         // test results for test values provided as part of the WMM2015 Report
205         // using height above MSL instead of height above ellipsoid
206         // the results have been obtained from the NOAA online calculator:
207         // http://www.ngdc.noaa.gov/geomag-web/#igrfwmm
208         final double[][] testValues = {
209             // Date     Alt              Lat                       Lon             X        Y         Z        H        F                  I                            D
210             //            m              rad                       rad            nT       nT        nT       nT       nT                 rad                         rad
211             {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)},
212             {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)},
213             {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)}
214         };
215 
216         GravityFieldFactory.clearPotentialCoefficientsReaders();
217         GravityFieldFactory.addPotentialCoefficientsReader(new EGMFormatReader("egm96", false));
218         final NormalizedSphericalHarmonicsProvider potential = GravityFieldFactory.getNormalizedProvider(maxDegree, maxOrder);
219         final Geoid geoid = new Geoid(potential, WGS84);
220 
221         final double eps = 1e-10;
222         final double degreeEps = 1e-2;
223         for (final double[] testValue : testValues) {
224             final AbsoluteDate date = new AbsoluteDate(2015, 1, 1, TimeScalesFactory.getUTC());
225             final GeoMagneticField model = GeoMagneticFieldFactory.getWMM(testValue[0]);
226             final double undulation = geoid.getUndulation(testValue[2], testValue[3], date);
227             final GeoMagneticElements result = model.calculateField(testValue[2], testValue[3], testValue[1] + undulation);
228 
229             // X
230             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(testValue[4]),
231                                     result.getFieldVector().getX(), eps);
232             // Y
233             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(testValue[5]),
234                                     result.getFieldVector().getY(), eps);
235             // Z
236             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(testValue[6]),
237                                     result.getFieldVector().getZ(), eps);
238             // H
239             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(testValue[7]),
240                                     result.getHorizontalIntensity(), eps);
241             // F
242             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(testValue[8]),
243                                     result.getTotalIntensity(), eps);
244             // inclination
245             Assertions.assertEquals(testValue[9], result.getInclination(), degreeEps);
246             // declination
247             Assertions.assertEquals(testValue[10], result.getDeclination(), degreeEps);
248         }
249     }
250 
251     @Test
252     public void testIGRF() throws Exception {
253         // test values from sample coordinate file
254         // provided as part of the geomag 7.0 distribution available at
255         // http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html
256         // modification: the julian day calculation of geomag is slightly different
257         // to the one from the WMM code, we use the WMM convention thus the outputs
258         // have been adapted.
259         runSampleFile(FieldModel.IGRF, "sample_coords.txt", "sample_out_IGRF12.txt");
260 
261         final double eps = 1e-10;
262         final double degreeEps = 1e-2;
263         for (final double[] igrfTestValue : igrfTestValues) {
264             final GeoMagneticField model = GeoMagneticFieldFactory.getIGRF(igrfTestValue[0]);
265             final GeoMagneticElements result = model.calculateField(igrfTestValue[2], igrfTestValue[3], igrfTestValue[1]);
266 
267             final Vector3D b = result.getFieldVector();
268 
269             // X
270             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(igrfTestValue[4]), b.getX(), eps);
271             // Y
272             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(igrfTestValue[5]), b.getY(), eps);
273             // Z
274             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(igrfTestValue[6]), b.getZ(), eps);
275             // H
276             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(igrfTestValue[7]),
277                                     result.getHorizontalIntensity(), eps);
278             // F
279             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(igrfTestValue[8]),
280                                     result.getTotalIntensity(), eps);
281             // inclination
282             Assertions.assertEquals(igrfTestValue[9], result.getInclination(), degreeEps);
283             // declination
284             Assertions.assertEquals(igrfTestValue[10], result.getDeclination(), degreeEps);
285         }
286     }
287 
288     @Test
289     public void testUnsupportedTransform() {
290         Assertions.assertThrows(OrekitException.class, () -> {
291             final GeoMagneticField model = GeoMagneticFieldFactory.getIGRF(1910);
292 
293             // the IGRF model of 1910 does not have secular variation, thus time transformation is not supported
294             model.transformModel(1950);
295         });
296     }
297 
298     @Test
299     public void testOutsideValidityTransform() {
300         Assertions.assertThrows(OrekitException.class, () -> {
301             final GeoMagneticField model1 = GeoMagneticFieldFactory.getIGRF(2005);
302             final GeoMagneticField model2 = GeoMagneticFieldFactory.getIGRF(2010);
303 
304             // the interpolation transformation is only allowed between 2005 and 2010
305             model1.transformModel(model2, 2012);
306         });
307     }
308 
309     @Test
310     public void testValidTransform() {
311         final GeoMagneticField model = GeoMagneticFieldFactory.getWMM(2015);
312 
313         Assertions.assertTrue(model.supportsTimeTransform());
314 
315         final GeoMagneticField transformedModel = model.transformModel(2017);
316 
317         Assertions.assertEquals(2015, transformedModel.validFrom(), 1e0);
318         Assertions.assertEquals(2020, transformedModel.validTo(), 1e0);
319         Assertions.assertEquals(2017, transformedModel.getEpoch(), 1e0);
320     }
321 
322     @Test
323     public void testParseOriginalWMMModel() throws Exception {
324         final String name = "WMM2015.COF";
325         Collection<GeoMagneticField> models = new GeoMagneticModelParser().
326                                               parse(new DataSource(name, () -> getResource(name)));
327         Assertions.assertNotNull(models);
328         Assertions.assertEquals(1, models.size());
329 
330         GeoMagneticField wmmModel = models.iterator().next();
331         Assertions.assertEquals("WMM-2015", wmmModel.getModelName());
332         Assertions.assertEquals(2015, wmmModel.getEpoch(), 1e-9);
333 
334         final double eps = 1e-10;
335         final double degreeEps = 1e-2;
336         for (final double[] wmmTestValue : wmmTestValues) {
337             if (wmmTestValue[0] != wmmModel.getEpoch()) {
338                 continue;
339             }
340             final GeoMagneticElements result = wmmModel.calculateField(wmmTestValue[2], wmmTestValue[3], wmmTestValue[1]);
341 
342             // X
343             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(wmmTestValue[4]),
344                                     result.getFieldVector().getX(), eps);
345             // Y
346             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(wmmTestValue[5]),
347                                     result.getFieldVector().getY(), eps);
348             // Z
349             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(wmmTestValue[6]),
350                                     result.getFieldVector().getZ(), eps);
351             // H
352             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(wmmTestValue[7]),
353                                     result.getHorizontalIntensity(), eps);
354             // F
355             Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(wmmTestValue[8]),
356                                     result.getTotalIntensity(), eps);
357             // inclination
358             Assertions.assertEquals(wmmTestValue[9], result.getInclination(), degreeEps);
359             // declination
360             Assertions.assertEquals(wmmTestValue[10], result.getDeclination(), degreeEps);
361         }
362     }
363 
364     public void runSampleFile(final FieldModel type, final String inputFile, final String outputFile)
365         throws Exception {
366 
367         final BufferedReader inReader = new BufferedReader(new InputStreamReader(getResource(inputFile)));
368         final BufferedReader outReader = new BufferedReader(new InputStreamReader(getResource(outputFile)));
369 
370         // read header line
371         outReader.readLine();
372 
373         String line;
374         while ((line = inReader.readLine()) != null) {
375             if (line.trim().isEmpty()) {
376                 break;
377             }
378 
379             final StringTokenizer st = new StringTokenizer(line);
380 
381             final double year = getYear(st.nextToken());
382             final String heightType = st.nextToken();
383             final String heightStr = st.nextToken();
384             final double lat = getLatLon(st.nextToken());
385             final double lon = getLatLon(st.nextToken());
386 
387             final GeoMagneticField field = GeoMagneticFieldFactory.getField(type, year);
388 
389             double height = Double.parseDouble(heightStr.substring(1));
390             if (heightStr.startsWith("K")) {
391                 // convert from km to m
392                 height *= 1000d;
393             } else if (heightStr.startsWith("F")) {
394                 // convert from feet to m
395                 height *= 0.3048;
396             }
397 
398             // Calculate the magnetic field with SI base unit inputs
399             final GeoMagneticElements ge = field.calculateField(lat, lon, height);
400             final String validateLine = outReader.readLine();
401             // geocentric altitude is not yet supported, ignore by now
402             if (!heightType.startsWith("C")) {
403                 validate(ge, validateLine);
404             }
405 
406             String geString = ge.toString();
407             Assertions.assertNotNull(geString);
408             Assertions.assertFalse(geString.isEmpty());
409         }
410     }
411 
412     private double getYear(final String yearStr) {
413 
414         if (yearStr.contains(",")) {
415             final StringTokenizer st = new StringTokenizer(yearStr, ",");
416             final String y = st.nextToken();
417             final String m = st.nextToken();
418             final String d = st.nextToken();
419             return GeoMagneticField.getDecimalYear(Integer.parseInt(d), Integer.parseInt(m), Integer.parseInt(y));
420         } else {
421             return Double.parseDouble(yearStr);
422         }
423     }
424 
425     private double getLatLon(final String str) {
426 
427         if (str.contains(",")) {
428             final StringTokenizer st = new StringTokenizer(str, ",");
429             final int d = Integer.parseInt(st.nextToken());
430             final int m = Integer.parseInt(st.nextToken());
431             int s = 0;
432             if (st.hasMoreTokens()) {
433                 s = Integer.parseInt(st.nextToken());
434             }
435             double deg = FastMath.abs(d) + m / 60d + s / 3600d;
436             if (d < 0) {
437                 deg = -deg;
438             }
439             return FastMath.toRadians(deg);
440         } else {
441             return FastMath.toRadians(Double.parseDouble(str));
442         }
443     }
444 
445     private double getRadians(final String degree, final String minute) {
446         double result = Double.parseDouble(degree.substring(0, degree.length() - 1));
447         final double min = Double.parseDouble(minute.substring(0, minute.length() - 1)) / 60d;
448         result += (result < 0) ? -min : min;
449         return FastMath.toRadians(result);
450     }
451 
452     @SuppressWarnings("unused")
453     private void validate(final GeoMagneticElements ge, final String outputLine) {
454 
455         final StringTokenizer st = new StringTokenizer(outputLine);
456 
457         final double year = getYear(st.nextToken());
458 
459         final String coord = st.nextToken();
460         final String heightStr = st.nextToken();
461         final String latStr = st.nextToken();
462         final String lonStr = st.nextToken();
463 
464         final double dec = getRadians(st.nextToken(), st.nextToken());
465         final double inc = getRadians(st.nextToken(), st.nextToken());
466 
467         final double h = Double.parseDouble(st.nextToken());
468         final double x = Double.parseDouble(st.nextToken());
469         final double y = Double.parseDouble(st.nextToken());
470         final double z = Double.parseDouble(st.nextToken());
471         final double f = Double.parseDouble(st.nextToken());
472 
473         final double eps = 1e-1;
474         Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(h), ge.getHorizontalIntensity(), eps);
475         Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(f), ge.getTotalIntensity(), eps);
476         Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(x), ge.getFieldVector().getX(), eps);
477         Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(y), ge.getFieldVector().getY(), eps);
478         Assertions.assertEquals(UnitsConverter.NANO_TESLAS_TO_TESLAS.convert(z), ge.getFieldVector().getZ(), eps);
479         Assertions.assertEquals(dec, ge.getDeclination(), eps);
480         Assertions.assertEquals(inc, ge.getInclination(), eps);
481     }
482 
483     private InputStream getResource(final String name) throws IOException {
484         // the data path has multiple components, the resources are in the first one
485         final String dataPath = System.getProperty(DataProvidersManager.OREKIT_DATA_PATH).
486                                 split(File.pathSeparator)[0];
487         return Files.newInputStream(new File(dataPath, name).toPath());
488     }
489 
490 }