1   /* Copyright 2002-2025 CS GROUP
2    * Licensed to CS GROUP (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.ionosphere;
18  
19  import java.lang.reflect.InvocationTargetException;
20  import java.lang.reflect.Method;
21  import java.net.URI;
22  import java.net.URISyntaxException;
23  import java.net.URL;
24  import java.nio.file.Paths;
25  
26  import org.hipparchus.CalculusFieldElement;
27  import org.hipparchus.Field;
28  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
29  import org.hipparchus.geometry.euclidean.threed.Vector3D;
30  import org.hipparchus.util.Binary64;
31  import org.hipparchus.util.Binary64Field;
32  import org.hipparchus.util.FastMath;
33  import org.hipparchus.util.MathUtils;
34  import org.junit.jupiter.api.Assertions;
35  import org.junit.jupiter.api.BeforeEach;
36  import org.junit.jupiter.api.Test;
37  import org.orekit.Utils;
38  import org.orekit.bodies.FieldGeodeticPoint;
39  import org.orekit.bodies.GeodeticPoint;
40  import org.orekit.bodies.OneAxisEllipsoid;
41  import org.orekit.data.DataContext;
42  import org.orekit.data.DataSource;
43  import org.orekit.data.DirectoryCrawlerTest;
44  import org.orekit.errors.OrekitException;
45  import org.orekit.errors.OrekitMessages;
46  import org.orekit.frames.Frame;
47  import org.orekit.frames.FramesFactory;
48  import org.orekit.frames.TopocentricFrame;
49  import org.orekit.gnss.PredefinedGnssSignal;
50  import org.orekit.models.earth.Geoid;
51  import org.orekit.models.earth.ReferenceEllipsoid;
52  import org.orekit.orbits.CartesianOrbit;
53  import org.orekit.orbits.FieldCartesianOrbit;
54  import org.orekit.orbits.FieldKeplerianOrbit;
55  import org.orekit.orbits.FieldOrbit;
56  import org.orekit.orbits.KeplerianOrbit;
57  import org.orekit.orbits.Orbit;
58  import org.orekit.orbits.PositionAngleType;
59  import org.orekit.propagation.FieldSpacecraftState;
60  import org.orekit.propagation.SpacecraftState;
61  import org.orekit.time.AbsoluteDate;
62  import org.orekit.time.FieldAbsoluteDate;
63  import org.orekit.time.TimeScale;
64  import org.orekit.time.TimeScalesFactory;
65  import org.orekit.utils.Constants;
66  import org.orekit.utils.FieldPVCoordinates;
67  import org.orekit.utils.IERSConventions;
68  import org.orekit.utils.PVCoordinates;
69  
70  public class GlobalIonosphereMapModelTest {
71  
72      private static final double epsilonParser = 1.0e-16;
73      private static final double epsilonDelay  = 0.001;
74      private SpacecraftState state;
75      private OneAxisEllipsoid earth;
76  
77      @BeforeEach
78      public void setUp() {
79          Utils.setDataRoot("regular-data:ionex");
80          final Orbit orbit = new KeplerianOrbit(24464560.0, 0.0, 1.122138, 1.10686, 1.00681,
81                                                 0.048363, PositionAngleType.MEAN,
82                                                 FramesFactory.getEME2000(),
83                                                 new AbsoluteDate(2019, 1, 14, 23, 59, 59.0, TimeScalesFactory.getUTC()),
84                                                 Constants.WGS84_EARTH_MU);
85          state = new SpacecraftState(orbit);
86          earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
87                                       Constants.WGS84_EARTH_FLATTENING,
88                                       FramesFactory.getITRF(IERSConventions.IERS_2010, true));
89      }
90  
91      @Test
92      public void testDelayAtIPP() {
93          GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
94          final double latitude  = FastMath.toRadians(30.0);
95          final double longitude = FastMath.toRadians(-130.0);
96          try {
97              final Method pathDelay = GlobalIonosphereMapModel.class.getDeclaredMethod("pathDelayAtIPP",
98                                                                                        AbsoluteDate.class,
99                                                                                        GeodeticPoint.class,
100                                                                                       Double.TYPE, Double.TYPE);
101             pathDelay.setAccessible(true);
102             final double delay = (Double) pathDelay.invoke(model,
103                                                            new AbsoluteDate(2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC()),
104                                                            new GeodeticPoint(latitude, longitude, 0.0),
105                                                            0.5 * FastMath.PI, PredefinedGnssSignal.G01.getFrequency());
106             Assertions.assertEquals(1.557, delay, epsilonDelay);
107         } catch (NoSuchMethodException | SecurityException | IllegalAccessException |
108                  IllegalArgumentException | InvocationTargetException e) {
109             Assertions.fail(e);
110         }
111     }
112 
113     @Test
114     public void testFieldDelayAtIPP() {
115         doTestFieldDelayAtIPP(Binary64Field.getInstance());
116     }
117 
118     private <T extends CalculusFieldElement<T>> void doTestFieldDelayAtIPP(final Field<T> field) {
119         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
120         final T zero = field.getZero();
121         final double latitude  = FastMath.toRadians(30.0);
122         final double longitude = FastMath.toRadians(-130.0);
123         try {
124             final Method pathDelay = GlobalIonosphereMapModel.class.getDeclaredMethod("pathDelayAtIPP",
125                                                                                       FieldAbsoluteDate.class,
126                                                                                       FieldGeodeticPoint.class,
127                                                                                       CalculusFieldElement.class,
128                                                                                       Double.TYPE);
129             pathDelay.setAccessible(true);
130             @SuppressWarnings("unchecked")
131             final T delay = (T) pathDelay.invoke(model,
132                                                  new FieldAbsoluteDate<>(field, 2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC()),
133                                                  new FieldGeodeticPoint<>(zero.newInstance(latitude),
134                                                                           zero.newInstance(longitude),
135                                                                           zero),
136                                                  zero.add(0.5 * FastMath.PI),
137                                                  PredefinedGnssSignal.G01.getFrequency());
138             Assertions.assertEquals(1.557, delay.getReal(), epsilonDelay);
139         } catch (NoSuchMethodException | SecurityException | IllegalAccessException |
140                         IllegalArgumentException | InvocationTargetException e) {
141             Assertions.fail(e);
142         }
143     }
144 
145     @Test
146     public void testSpacecraftState() {
147         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
148         double a = 7187990.1979844316;
149         double e = 0.5e-4;
150         double i = FastMath.toRadians(98.01);
151         double omega = FastMath.toRadians(131.88);
152         double OMEGA = FastMath.toRadians(252.24);
153         double lv = FastMath.toRadians(250.00);
154 
155         AbsoluteDate date = new AbsoluteDate(2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC());
156         final SpacecraftState state =
157                         new SpacecraftState(new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
158                                                                FramesFactory.getEME2000(), date,
159                                                                Constants.EIGEN5C_EARTH_MU));
160         final TopocentricFrame topo = new TopocentricFrame(earth,
161                                                            new GeodeticPoint(FastMath.toRadians(20.5236),
162                                                                              FastMath.toRadians(85.7881),
163                                                                              36.0),
164                                                            "Cuttack");
165         final double delay = model.pathDelay(state, topo, PredefinedGnssSignal.G01.getFrequency(), null);
166         Assertions.assertEquals(2.811, delay, epsilonDelay);
167 
168         // the delay at station longitude is different, due to IPP
169         try {
170             final Method pathDelay = GlobalIonosphereMapModel.class.getDeclaredMethod("pathDelayAtIPP",
171                                                                                       AbsoluteDate.class,
172                                                                                       GeodeticPoint.class,
173                                                                                       Double.TYPE, Double.TYPE);
174             pathDelay.setAccessible(true);
175             final double delayIPP = (Double) pathDelay.invoke(model, date, topo.getPoint(),
176                                                               0.5 * FastMath.PI,
177                                                               PredefinedGnssSignal.G01.getFrequency());
178             Assertions.assertEquals(2.173, delayIPP, epsilonDelay);
179         } catch (NoSuchMethodException | SecurityException | IllegalAccessException |
180                  IllegalArgumentException | InvocationTargetException ex) {
181             Assertions.fail(ex);
182         }
183     }
184 
185     @Test
186     public void testAboveIono() {
187         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
188         double a = 7187990.1979844316;
189         double e = 0.5e-4;
190         double i = FastMath.toRadians(98.01);
191         double omega = FastMath.toRadians(131.88);
192         double OMEGA = FastMath.toRadians(252.24);
193         double lv = FastMath.toRadians(250.00);
194 
195         AbsoluteDate date = new AbsoluteDate(2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC());
196         final SpacecraftState state =
197                         new SpacecraftState(new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
198                                                                FramesFactory.getEME2000(), date,
199                                                                Constants.EIGEN5C_EARTH_MU));
200         final TopocentricFrame topo = new TopocentricFrame(earth,
201                                                            new GeodeticPoint(FastMath.toRadians(20.5236),
202                                                                              FastMath.toRadians(85.7881),
203                                                                              650000.0),
204                                                            "very-high");
205         final double delay = model.pathDelay(state, topo, PredefinedGnssSignal.G01.getFrequency(), null);
206         Assertions.assertEquals(0.0, delay, epsilonDelay);
207 
208     }
209 
210     @Test
211     public void testSpacecraftStateField() {
212         doTestSpacecraftStateField(Binary64Field.getInstance());
213     }
214 
215     private <T extends CalculusFieldElement<T>> void doTestSpacecraftStateField(final Field<T> field) {
216         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
217         final T zero = field.getZero();
218         T a = zero.newInstance(7187990.1979844316);
219         T e = zero.newInstance(0.5e-4);
220         T i = FastMath.toRadians(zero.newInstance(98.01));
221         T omega = FastMath.toRadians(zero.newInstance(131.88));
222         T OMEGA = FastMath.toRadians(zero.newInstance(252.24));
223         T lv = FastMath.toRadians(zero.newInstance(250.00));
224 
225         FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field,
226                                                             new AbsoluteDate(2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC()));
227         final FieldSpacecraftState<T> state =
228                         new FieldSpacecraftState<>(new FieldKeplerianOrbit<>(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
229                                                                              FramesFactory.getEME2000(), date,
230                                                                              zero.newInstance(Constants.EIGEN5C_EARTH_MU)));
231         final TopocentricFrame topo = new TopocentricFrame(earth,
232                                                            new GeodeticPoint(FastMath.toRadians(20.5236),
233                                                                              FastMath.toRadians(85.7881),
234                                                                              36.0),
235                                                            "Cuttack");
236         final T delay = model.pathDelay(state, topo, PredefinedGnssSignal.G01.getFrequency(), null);
237         Assertions.assertEquals(2.811, delay.getReal(), epsilonDelay);
238 
239         // the delay at station longitude is different, due to IPP
240         try {
241             final Method pathDelay = GlobalIonosphereMapModel.class.getDeclaredMethod("pathDelayAtIPP",
242                                                                                       FieldAbsoluteDate.class,
243                                                                                       FieldGeodeticPoint.class,
244                                                                                       CalculusFieldElement.class,
245                                                                                       Double.TYPE);
246             pathDelay.setAccessible(true);
247             @SuppressWarnings("unchecked")
248             final T delayIPP = (T) pathDelay.invoke(model, date,
249                                                     new FieldGeodeticPoint<>(field, topo.getPoint()),
250                                                     zero.newInstance(0.5 * FastMath.PI),
251                                                     PredefinedGnssSignal.G01.getFrequency());
252             Assertions.assertEquals(2.173, delayIPP.getReal(), epsilonDelay);
253         } catch (NoSuchMethodException | SecurityException | IllegalAccessException |
254                  IllegalArgumentException | InvocationTargetException ex) {
255             Assertions.fail(ex);
256         }
257     }
258 
259     @Test
260     public void testAboveIonoField() {
261         doTestAboveIonoField(Binary64Field.getInstance());
262     }
263 
264     private <T extends CalculusFieldElement<T>> void doTestAboveIonoField(final Field<T> field) {
265         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
266         final T zero = field.getZero();
267         T a = zero.newInstance(7187990.1979844316);
268         T e = zero.newInstance(0.5e-4);
269         T i = FastMath.toRadians(zero.newInstance(98.01));
270         T omega = FastMath.toRadians(zero.newInstance(131.88));
271         T OMEGA = FastMath.toRadians(zero.newInstance(252.24));
272         T lv = FastMath.toRadians(zero.newInstance(250.00));
273 
274         FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field,
275                                                             new AbsoluteDate(2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC()));
276         final FieldSpacecraftState<T> state =
277                         new FieldSpacecraftState<>(new FieldKeplerianOrbit<>(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
278                                                                              FramesFactory.getEME2000(), date,
279                                                                              zero.newInstance(Constants.EIGEN5C_EARTH_MU)));
280         final TopocentricFrame topo = new TopocentricFrame(earth,
281                                                            new GeodeticPoint(FastMath.toRadians(20.5236),
282                                                                              FastMath.toRadians(85.7881),
283                                                                              650000.0),
284                                                            "very-high");
285         final T delay = model.pathDelay(state, topo, PredefinedGnssSignal.G01.getFrequency(), null);
286         Assertions.assertEquals(0.0, delay.getReal(), epsilonDelay);
287 
288     }
289 
290     @Test
291     public void testParser() throws URISyntaxException {
292 
293         Utils.setDataRoot("regular-data");
294         URL url1 = DirectoryCrawlerTest.class.getClassLoader().getResource("ionex/split-1.19i");
295         DataSource ds1 = new DataSource(Paths.get(url1.toURI()).toString());
296         URL url2 = DirectoryCrawlerTest.class.getClassLoader().getResource("ionex/split-2.19i");
297         DataSource ds2 = new DataSource(Paths.get(url2.toURI()).toString());
298         URL url3 = DirectoryCrawlerTest.class.getClassLoader().getResource("ionex/split-3.19i");
299         DataSource ds3 = new DataSource(Paths.get(url3.toURI()).toString());
300         GlobalIonosphereMapModel model =
301             new GlobalIonosphereMapModel(TimeScalesFactory.getUTC(),
302                                          GlobalIonosphereMapModel.TimeInterpolator.SIMPLE_LINEAR,
303                                          ds1, ds2, ds3);
304 
305         // Commons parameters
306         final double latitude = FastMath.toRadians(45.0);
307 
308         double longitude1;
309         double longitude2;
310 
311         try {
312             final Method pathDelay = GlobalIonosphereMapModel.class.getDeclaredMethod("pathDelayAtIPP",
313                                                                                       AbsoluteDate.class,
314                                                                                       GeodeticPoint.class,
315                                                                                       Double.TYPE, Double.TYPE);
316             pathDelay.setAccessible(true);
317 
318             // Test longitude = 181° and longitude = -179°
319             longitude1 = FastMath.toRadians(181.0);
320             longitude2 = FastMath.toRadians(-179.0);
321             AbsoluteDate date1 = new AbsoluteDate(2019, 1, 15, 1, 0, 0.0, TimeScalesFactory.getUTC());
322             Assertions.assertEquals((Double) pathDelay.invoke(model, date1, new GeodeticPoint(latitude, longitude1, 0.0),
323                                                               0.01, PredefinedGnssSignal.G01.getFrequency()),
324                                     ((Double) pathDelay.invoke(model, date1, new GeodeticPoint(latitude, longitude2, 0.0),
325                                                                0.01, PredefinedGnssSignal.G01.getFrequency())),
326                                     epsilonParser);
327 
328             // Test longitude = 180° and longitude = -180°
329             AbsoluteDate date2 = new AbsoluteDate(2019, 1, 15, 3, 0, 0.0, TimeScalesFactory.getUTC());
330             longitude1 = FastMath.toRadians(180.0);
331             longitude2 = FastMath.toRadians(-180.0);
332 
333             Assertions.assertEquals(((Double) pathDelay.invoke(model, date2, new GeodeticPoint(latitude, longitude1, 0.0),
334                                                                0.01, PredefinedGnssSignal.G01.getFrequency())),
335                                     ((Double) pathDelay.invoke(model, date2, new GeodeticPoint(latitude, longitude2, 0.0),
336                                                                0.01, PredefinedGnssSignal.G01.getFrequency())),
337                                     epsilonParser);
338 
339             // Test longitude = 0° and longitude = 360°
340             AbsoluteDate date3 = new AbsoluteDate(2019, 1, 15, 5, 0, 0.0, TimeScalesFactory.getUTC());
341             longitude1 =  FastMath.toRadians(0.);
342             longitude2 =  FastMath.toRadians(360.0);
343 
344             Assertions.assertEquals(((Double) pathDelay.invoke(model, date3, new GeodeticPoint(latitude, longitude1, 0.0),
345                                                                0.01, PredefinedGnssSignal.G01.getFrequency())),
346                                     ((Double) pathDelay.invoke(model, date3, new GeodeticPoint(latitude, longitude2, 0.0),
347                                                                0.01, PredefinedGnssSignal.G01.getFrequency())),
348                                     epsilonParser);
349 
350         } catch (NoSuchMethodException | SecurityException | IllegalAccessException |
351                         IllegalArgumentException | InvocationTargetException e) {
352             Assertions.fail(e);
353         }
354     }
355 
356     @Test
357     public void testCorruptedFileBadData() {
358         final String fileName = "corrupted-bad-data-gpsg0150.19i";
359 
360         try {
361             new GlobalIonosphereMapModel(fileName);
362             Assertions.fail("An exception should have been thrown");
363 
364         } catch (OrekitException oe) {
365             Assertions.assertEquals(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE, oe.getSpecifier());
366         }
367     }
368 
369     @Test
370     public void testEarlierDate() {
371         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
372         final double latitude  = FastMath.toRadians(60.0);
373         final double longitude = FastMath.toRadians(-130.0);
374         final GeodeticPoint point = new GeodeticPoint(latitude, longitude, 0.0);
375 
376         try {
377             model.pathDelay(state, new TopocentricFrame(earth, point, null),
378                             PredefinedGnssSignal.G01.getFrequency(),
379                             model.getParameters(new AbsoluteDate()));
380             Assertions.fail("An exception should have been thrown");
381         } catch (OrekitException oe) {
382             Assertions.assertEquals(OrekitMessages.NO_TEC_DATA_IN_FILES_FOR_DATE, oe.getSpecifier());
383         }
384     }
385 
386     @Test
387     public void testFieldEarlierDate() {
388         doTestFieldEarlierDate(Binary64Field.getInstance());
389     }
390 
391     private <T extends CalculusFieldElement<T>> void doTestFieldEarlierDate(final Field<T> field) {
392         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
393         final double latitude  = FastMath.toRadians(60.0);
394         final double longitude = FastMath.toRadians(-130.0);
395         final GeodeticPoint point = new GeodeticPoint(latitude, longitude, 0.0);
396 
397         try {
398             model.pathDelay(new FieldSpacecraftState<>(field, state), new TopocentricFrame(earth, point, null),
399                             PredefinedGnssSignal.G01.getFrequency(),
400                             model.getParameters(field, new FieldAbsoluteDate<>(field)));
401             Assertions.fail("An exception should have been thrown");
402         } catch (OrekitException oe) {
403             Assertions.assertEquals(OrekitMessages.NO_TEC_DATA_IN_FILES_FOR_DATE, oe.getSpecifier());
404         }
405     }
406 
407     @Deprecated
408     @Test
409     public void testDeprecated() throws URISyntaxException {
410         final TimeScale utc = DataContext.getDefault().getTimeScales().getUTC();
411         Assertions.assertEquals(GlobalIonosphereMapModel.TimeInterpolator.SIMPLE_LINEAR,
412                                 new GlobalIonosphereMapModel("gpsg0150.19i",
413                                                              DataContext.getDefault().getDataProvidersManager(),
414                                                              utc).getInterpolator());
415         final URL        url = DirectoryCrawlerTest.class.getClassLoader().getResource("ionex/gpsg0150.19i");
416         final DataSource ds  = new DataSource(url.toURI());
417         Assertions.assertEquals(GlobalIonosphereMapModel.TimeInterpolator.SIMPLE_LINEAR,
418                                 new GlobalIonosphereMapModel(utc, ds).getInterpolator());
419     }
420 
421     @Test
422     public void testTimeInterpolation() throws URISyntaxException {
423         final TimeScale  utc = DataContext.getDefault().getTimeScales().getUTC();
424         final URI        uri = DirectoryCrawlerTest.class.getClassLoader().getResource("ionex/gpsg0150.19i").toURI();
425         final GlobalIonosphereMapModel nearest =
426             new GlobalIonosphereMapModel(utc,
427                                          GlobalIonosphereMapModel.TimeInterpolator.NEAREST_MAP,
428                                          new DataSource(uri));
429         final GlobalIonosphereMapModel simple =
430             new GlobalIonosphereMapModel(utc,
431                                          GlobalIonosphereMapModel.TimeInterpolator.SIMPLE_LINEAR,
432                                          new DataSource(uri));
433         final GlobalIonosphereMapModel rotated =
434             new GlobalIonosphereMapModel(utc,
435                                          GlobalIonosphereMapModel.TimeInterpolator.ROTATED_LINEAR,
436                                          new DataSource(uri));
437         final OneAxisEllipsoid sphericalEarth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
438                                                                      0.0,
439                                                                      FramesFactory.getITRF(IERSConventions.IERS_2010,
440                                                                                            true));
441         final TopocentricFrame topo = new TopocentricFrame(sphericalEarth,
442                                                            new GeodeticPoint(FastMath.toRadians(55.0),
443                                                                              FastMath.toRadians(15.0),
444                                                                              0.0),
445                                                            null);
446         final double        frequency = PredefinedGnssSignal.G01.getFrequency();
447         final Frame         inertial  = FramesFactory.getEME2000();
448         final PVCoordinates above     = new PVCoordinates(new Vector3D(0, 0, 1.5e6), new Vector3D(6.4e3, 0, 0));
449         final double        scale     = 0.1;
450         final double        factor    = scale * 40.3e16 / (frequency * frequency);
451 
452         // check at one exact sampling date and for a vertical line of sight
453         final AbsoluteDate    tSampling     = new AbsoluteDate(2019, 1, 15, 14, 0, 0.0, utc);
454         final Orbit           oSampling     = new CartesianOrbit(topo.getTransformTo(inertial, tSampling).
455                                                                  transformPVCoordinates(above),
456                                                                  inertial, tSampling, Constants.WGS84_EARTH_MU);
457         final SpacecraftState sSampling     = new SpacecraftState(oSampling);
458         final double          gridValue     = 103.0;
459 
460         // at exact sampling dates, all interpolation methods should return the same value
461         Assertions.assertEquals(gridValue * factor,
462                                 nearest.pathDelay(sSampling, topo, frequency, nearest.getParameters(tSampling)),
463                                 1.0e-15);
464         Assertions.assertEquals(gridValue * factor,
465                                 simple.pathDelay(sSampling, topo, frequency, simple.getParameters(tSampling)),
466                                 1.0e-15);
467         Assertions.assertEquals(gridValue * factor,
468                                 rotated.pathDelay(sSampling, topo, frequency, rotated.getParameters(tSampling)),
469                                 1.0e-15);
470 
471         // check at midtime between sampling dates and for a vertical line of sight
472         // between sampling dates, interpolation methods should be different
473         final double halfStep = 1800.0;
474         final AbsoluteDate    tInterp   = tSampling.shiftedBy(halfStep);
475 
476         final Orbit           oInterp   = new CartesianOrbit(topo.getTransformTo(inertial, tInterp).
477                                                              transformPVCoordinates(above),
478                                                              inertial, tInterp, Constants.WGS84_EARTH_MU);
479         final SpacecraftState sInterp   = new SpacecraftState(oInterp);
480         final double          gridValue1 = 103.0;
481         final double          gridValue2 = 101.0;
482         Assertions.assertEquals(gridValue1 * factor,
483                                 nearest.pathDelay(sInterp, topo, frequency, nearest.getParameters(tInterp)),
484                                 1.0e-15);
485         Assertions.assertEquals((0.5 * gridValue1 + 0.5 * gridValue2) * factor,
486                                 simple.pathDelay(sInterp, topo, frequency, simple.getParameters(tInterp)),
487                                 1.0e-15);
488 
489         // the half hour offsets is about 7.52 degrees, so maps rotates about 1.504 cell size
490         // it is also exactly at mid-time between the two TEC maps
491         final double siderealDay   = MathUtils.TWO_PI / Constants.WGS84_EARTH_ANGULAR_VELOCITY;
492         final double cellDriftRate = 72.0 / siderealDay;
493         final double cellRatio     = halfStep * cellDriftRate - 1.0; // very slightly above 0.504
494         final double interpMapA    = ((1 - cellRatio) * 104.0 + cellRatio * 105.0) * factor;
495         final double interpMapB    = ((1 - cellRatio) * 101.0 + cellRatio * 100.0) * factor;
496         Assertions.assertEquals(0.5 * (interpMapA + interpMapB),
497                                 rotated.pathDelay(sInterp, topo, frequency, rotated.getParameters(tInterp)),
498                                 1.0e-15);
499 
500     }
501 
502     @Test
503     public void testFieldTimeInterpolation() throws URISyntaxException {
504         doTestFieldTimeInterpolation(Binary64Field.getInstance());
505     }
506 
507     private <T extends CalculusFieldElement<T>> void doTestFieldTimeInterpolation(Field<T> field) throws URISyntaxException {
508         final T zero = field.getZero();
509         final TimeScale  utc = DataContext.getDefault().getTimeScales().getUTC();
510         final URI        uri = DirectoryCrawlerTest.class.getClassLoader().getResource("ionex/gpsg0150.19i").toURI();
511         final GlobalIonosphereMapModel nearest =
512             new GlobalIonosphereMapModel(utc,
513                                          GlobalIonosphereMapModel.TimeInterpolator.NEAREST_MAP,
514                                          new DataSource(uri));
515         final GlobalIonosphereMapModel simple =
516             new GlobalIonosphereMapModel(utc,
517                                          GlobalIonosphereMapModel.TimeInterpolator.SIMPLE_LINEAR,
518                                          new DataSource(uri));
519         final GlobalIonosphereMapModel rotated =
520             new GlobalIonosphereMapModel(utc,
521                                          GlobalIonosphereMapModel.TimeInterpolator.ROTATED_LINEAR,
522                                          new DataSource(uri));
523         final OneAxisEllipsoid sphericalEarth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
524                                                                      0.0,
525                                                                      FramesFactory.getITRF(IERSConventions.IERS_2010,
526                                                                                            true));
527         final TopocentricFrame topo = new TopocentricFrame(sphericalEarth,
528                                                            new GeodeticPoint(FastMath.toRadians(55.0),
529                                                                              FastMath.toRadians(15.0),
530                                                                              0.0),
531                                                            null);
532         final double        frequency = PredefinedGnssSignal.G01.getFrequency();
533         final Frame         inertial  = FramesFactory.getEME2000();
534         final FieldPVCoordinates<T> above     = new FieldPVCoordinates<>(new FieldVector3D<>(zero,
535                                                                                              zero,
536                                                                                              zero.newInstance(1.5e6)),
537                                                                          new FieldVector3D<>(zero.newInstance(6.4e3),
538                                                                                              zero,
539                                                                                              zero));
540         final double        scale     = 0.1;
541         final double        factor    = scale * 40.3e16 / (frequency * frequency);
542 
543         // check at one exact sampling date and for a vertical line of sight
544         final FieldAbsoluteDate<T>    tSampling     = new FieldAbsoluteDate<>(field,
545                                                                               new AbsoluteDate(2019, 1, 15, 14, 0, 0.0, utc));
546         final FieldOrbit<T>           oSampling     = new FieldCartesianOrbit<>(topo.getTransformTo(inertial, tSampling).
547                                                                                 transformPVCoordinates(above),
548                                                                                 inertial, tSampling,
549                                                                                 zero.newInstance(Constants.WGS84_EARTH_MU));
550         final FieldSpacecraftState<T> sSampling     = new FieldSpacecraftState<>(oSampling);
551         final double          gridValue     = 103.0;
552 
553         // at exact sampling dates, all interpolation methods should return the same value
554         Assertions.assertEquals(gridValue * factor,
555                                 nearest.pathDelay(sSampling, topo, frequency, nearest.getParameters(field, tSampling)).getReal(),
556                                 1.0e-15);
557         Assertions.assertEquals(gridValue * factor,
558                                 simple.pathDelay(sSampling, topo, frequency, simple.getParameters(field, tSampling)).getReal(),
559                                 1.0e-15);
560         Assertions.assertEquals(gridValue * factor,
561                                 rotated.pathDelay(sSampling, topo, frequency, rotated.getParameters(field, tSampling)).getReal(),
562                                 1.0e-15);
563 
564         // check at midtime between sampling dates and for a vertical line of sight
565         // between sampling dates, interpolation methods should be different
566         final double halfStep = 1800.0;
567         final FieldAbsoluteDate<T>    tInterp   = tSampling.shiftedBy(halfStep);
568 
569         final FieldOrbit<T>           oInterp   = new FieldCartesianOrbit<>(topo.getTransformTo(inertial, tInterp).
570                                                                             transformPVCoordinates(above),
571                                                                             inertial, tInterp,
572                                                                             zero.newInstance(Constants.WGS84_EARTH_MU));
573         final FieldSpacecraftState<T> sInterp   = new FieldSpacecraftState<>(oInterp);
574         final double          gridValue1 = 103.0;
575         final double          gridValue2 = 101.0;
576         Assertions.assertEquals(gridValue1 * factor,
577                                 nearest.pathDelay(sInterp, topo, frequency, nearest.getParameters(field, tInterp)).getReal(),
578                                 1.0e-15);
579         Assertions.assertEquals((0.5 * gridValue1 + 0.5 * gridValue2) * factor,
580                                 simple.pathDelay(sInterp, topo, frequency, simple.getParameters(field, tInterp)).getReal(),
581                                 1.0e-15);
582 
583         // the half hour offsets is about 7.52 degrees, so maps rotates about 1.504 cell size
584         // it is also exactly at mid-time between the two TEC maps
585         final double siderealDay   = MathUtils.TWO_PI / Constants.WGS84_EARTH_ANGULAR_VELOCITY;
586         final double cellDriftRate = 72.0 / siderealDay;
587         final double cellRatio     = halfStep * cellDriftRate - 1.0; // very slightly above 0.504
588         final double interpMapA    = ((1 - cellRatio) * 104.0 + cellRatio * 105.0) * factor;
589         final double interpMapB    = ((1 - cellRatio) * 101.0 + cellRatio * 100.0) * factor;
590         Assertions.assertEquals(0.5 * (interpMapA + interpMapB),
591                                 rotated.pathDelay(sInterp, topo, frequency, rotated.getParameters(field, tInterp)).getReal(),
592                                 1.0e-15);
593 
594     }
595 
596     @Test
597     public void testNotEllipsoid() throws URISyntaxException {
598         Utils.setDataRoot("regular-data:ionex:potential");
599         final DataContext dataContext = DataContext.getDefault();
600         final TimeScale   utc         = dataContext.getTimeScales().getUTC();
601         final URI         uri         = DirectoryCrawlerTest.class.getClassLoader().getResource("ionex/gpsg0150.19i").toURI();
602         final GlobalIonosphereMapModel gim =
603             new GlobalIonosphereMapModel(utc,
604                                          GlobalIonosphereMapModel.TimeInterpolator.ROTATED_LINEAR,
605                                          new DataSource(uri));
606         try {
607             final SpacecraftState shifted = state.shiftedBy(2000.0);
608             final Frame       itrf        = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
609             final Geoid       geoid       = new Geoid(dataContext.getGravityFields().getNormalizedProvider(2, 0),
610                                                       ReferenceEllipsoid.getWgs84(itrf));
611             gim.pathDelay(shifted,
612                           new TopocentricFrame(geoid, new GeodeticPoint(0, 0, 0.0), "dummy"),
613                           PredefinedGnssSignal.G01.getFrequency(), gim.getParameters(shifted.getDate()));
614             Assertions.fail("an exception should have been thrown");
615         } catch (OrekitException oe) {
616             Assertions.assertEquals(OrekitMessages.BODY_SHAPE_MUST_BE_A_ONE_AXIS_ELLIPSOID, oe.getSpecifier());
617         }
618     }
619 
620    @Test
621     public void testFieldNotEllipsoid() throws URISyntaxException {
622         doTestFieldNotEllipsoid(Binary64Field.getInstance());
623    }
624 
625    private <T extends CalculusFieldElement<T>> void doTestFieldNotEllipsoid(final Field<T> field)
626        throws URISyntaxException {
627         Utils.setDataRoot("regular-data:ionex:potential");
628         final DataContext dataContext = DataContext.getDefault();
629         final TimeScale   utc         = dataContext.getTimeScales().getUTC();
630         final URI         uri         = DirectoryCrawlerTest.class.getClassLoader().getResource("ionex/gpsg0150.19i").toURI();
631         final GlobalIonosphereMapModel gim =
632             new GlobalIonosphereMapModel(utc,
633                                          GlobalIonosphereMapModel.TimeInterpolator.ROTATED_LINEAR,
634                                          new DataSource(uri));
635         try {
636             final FieldSpacecraftState<T> shifted = new FieldSpacecraftState<>(field, state.shiftedBy(2000.0));
637             final Frame       itrf        = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
638             final Geoid       geoid       = new Geoid(dataContext.getGravityFields().getNormalizedProvider(2, 0),
639                                                       ReferenceEllipsoid.getWgs84(itrf));
640             gim.pathDelay(shifted,
641                           new TopocentricFrame(geoid, new GeodeticPoint(0, 0, 0.0), "dummy"),
642                           PredefinedGnssSignal.G01.getFrequency(), gim.getParameters(field, shifted.getDate()));
643             Assertions.fail("an exception should have been thrown");
644         } catch (OrekitException oe) {
645             Assertions.assertEquals(OrekitMessages.BODY_SHAPE_MUST_BE_A_ONE_AXIS_ELLIPSOID, oe.getSpecifier());
646         }
647     }
648 
649     @Test
650     public void testLaterDate() {
651         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
652         final double latitude  = FastMath.toRadians(60.0);
653         final double longitude = FastMath.toRadians(-130.0);
654         final GeodeticPoint point = new GeodeticPoint(latitude, longitude, 0.0);
655 
656         try {
657             model.pathDelay(state, new TopocentricFrame(earth, point, null),
658                             PredefinedGnssSignal.G01.getFrequency(),
659                             model.getParameters(new AbsoluteDate()));
660             Assertions.fail("An exception should have been thrown");
661         } catch (OrekitException oe) {
662             Assertions.assertEquals(OrekitMessages.NO_TEC_DATA_IN_FILES_FOR_DATE, oe.getSpecifier());
663         }
664 
665     }
666 
667     @Test
668     public void testFieldLaterDate() {
669         doTestFieldLaterDate(Binary64Field.getInstance());
670     }
671 
672     private <T extends CalculusFieldElement<T>> void doTestFieldLaterDate(final Field<T> field) {
673         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
674         final double latitude  = FastMath.toRadians(60.0);
675         final double longitude = FastMath.toRadians(-130.0);
676         final GeodeticPoint point = new GeodeticPoint(latitude, longitude, 0.0);
677         try {
678             model.pathDelay(new FieldSpacecraftState<>(field, state), new TopocentricFrame(earth, point, null),
679                             PredefinedGnssSignal.G01.getFrequency(),
680                             model.getParameters(field, new FieldAbsoluteDate<>(field)));
681             Assertions.fail("An exception should have been thrown");
682         } catch (OrekitException oe) {
683             Assertions.assertEquals(OrekitMessages.NO_TEC_DATA_IN_FILES_FOR_DATE, oe.getSpecifier());
684         }
685 
686     }
687 
688     @Test
689     public void testLastDate() {
690         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
691         final double latitude  = FastMath.toRadians(60.0);
692         final double longitude = FastMath.toRadians(-130.0);
693         final GeodeticPoint point = new GeodeticPoint(latitude, longitude, 0.0);
694         final Orbit orbit = new KeplerianOrbit(24464560.0, 0.0, 1.122138, 1.10686, 1.00681,
695                                                0.048363, PositionAngleType.MEAN,
696                                                FramesFactory.getEME2000(),
697                                                new AbsoluteDate(2019, 1, 16, TimeScalesFactory.getUTC()),
698                                                Constants.WGS84_EARTH_MU);
699         final double delay = model.pathDelay(new SpacecraftState(orbit), new TopocentricFrame(earth, point, null),
700                                              PredefinedGnssSignal.G01.getFrequency(),
701                                              model.getParameters(orbit.getDate()));
702         Assertions.assertEquals(3.156, delay, 1.0e-3);
703     }
704 
705     @Test
706     public void testLastDateField() {
707         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
708         final Field<Binary64> field = Binary64Field.getInstance();
709         final double latitude  = FastMath.toRadians(60.0);
710         final double longitude = FastMath.toRadians(-130.0);
711         final GeodeticPoint point = new GeodeticPoint(latitude, longitude, 0.0);
712         final FieldOrbit<Binary64> orbit =
713                 new FieldKeplerianOrbit<>(field,
714                                           new KeplerianOrbit(24464560.0, 0.0, 1.122138, 1.10686, 1.00681,
715                                                               0.048363, PositionAngleType.MEAN,
716                                                               FramesFactory.getEME2000(),
717                                                               new AbsoluteDate(2019, 1, 16, TimeScalesFactory.getUTC()),
718                                                               Constants.WGS84_EARTH_MU));
719         final Binary64 delay = model.pathDelay(new FieldSpacecraftState<>(orbit),
720                                                new TopocentricFrame(earth, point, null),
721                                                PredefinedGnssSignal.G01.getFrequency(),
722                                                model.getParameters(field, orbit.getDate()));
723         Assertions.assertEquals(3.156, delay.getReal(), 1.0e-3);
724     }
725 
726     @Test
727     /**
728      * The goal of this test is to verify if an OrekitException is thrown when latitude or longitude
729      * boundaries are not present in the header section of the Global Ionosphere Map.
730      */
731     public void testIssue621() {
732         final String fileName  = "missing-lat-lon-header-gpsg0150.19i";
733 
734         try {
735             new GlobalIonosphereMapModel(fileName);
736             Assertions.fail("An exception should have been thrown");
737 
738         } catch (OrekitException oe) {
739             Assertions.assertEquals(OrekitMessages.NO_LATITUDE_LONGITUDE_BOUNDARIES_IN_IONEX_HEADER, oe.getSpecifier());
740         }
741     }
742 
743     @Test
744     public void testJammedFields() {
745         Assertions.assertNotNull(new GlobalIonosphereMapModel("jammed-fields.14i"));
746     }
747 
748     @Test
749     public void testMissingEpochInHeader() {
750         final String fileName  = "missing-epoch-header-gpsg0150.19i";
751 
752         try {
753             new GlobalIonosphereMapModel(fileName);
754             Assertions.fail("An exception should have been thrown");
755 
756         } catch (OrekitException oe) {
757             Assertions.assertEquals(OrekitMessages.NO_EPOCH_IN_IONEX_HEADER, oe.getSpecifier());
758         }
759 
760     }
761 
762 }