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.URISyntaxException;
22  import java.net.URL;
23  import java.nio.file.Paths;
24  
25  import org.hipparchus.CalculusFieldElement;
26  import org.hipparchus.Field;
27  import org.hipparchus.util.Binary64Field;
28  import org.hipparchus.util.FastMath;
29  import org.junit.jupiter.api.Assertions;
30  import org.junit.jupiter.api.BeforeEach;
31  import org.junit.jupiter.api.Test;
32  import org.orekit.Utils;
33  import org.orekit.bodies.GeodeticPoint;
34  import org.orekit.bodies.OneAxisEllipsoid;
35  import org.orekit.data.DataSource;
36  import org.orekit.data.DirectoryCrawlerTest;
37  import org.orekit.errors.OrekitException;
38  import org.orekit.errors.OrekitMessages;
39  import org.orekit.frames.FramesFactory;
40  import org.orekit.frames.TopocentricFrame;
41  import org.orekit.gnss.PredefinedGnssSignal;
42  import org.orekit.orbits.FieldKeplerianOrbit;
43  import org.orekit.orbits.KeplerianOrbit;
44  import org.orekit.orbits.Orbit;
45  import org.orekit.orbits.PositionAngleType;
46  import org.orekit.propagation.FieldSpacecraftState;
47  import org.orekit.propagation.SpacecraftState;
48  import org.orekit.time.AbsoluteDate;
49  import org.orekit.time.FieldAbsoluteDate;
50  import org.orekit.time.TimeScalesFactory;
51  import org.orekit.utils.Constants;
52  import org.orekit.utils.IERSConventions;
53  
54  public class GlobalIonosphereMapModelTest {
55  
56      private static final double epsilonParser = 1.0e-16;
57      private static final double epsilonDelay  = 0.001;
58      private SpacecraftState state;
59      private OneAxisEllipsoid earth;
60  
61      @BeforeEach
62      public void setUp() {
63          Utils.setDataRoot("regular-data:ionex");
64          final Orbit orbit = new KeplerianOrbit(24464560.0, 0.0, 1.122138, 1.10686, 1.00681,
65                                                 0.048363, PositionAngleType.MEAN,
66                                                 FramesFactory.getEME2000(),
67                                                 new AbsoluteDate(2019, 1, 14, 23, 59, 59.0, TimeScalesFactory.getUTC()),
68                                                 Constants.WGS84_EARTH_MU);
69          state = new SpacecraftState(orbit);
70          earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
71                                       Constants.WGS84_EARTH_FLATTENING,
72                                       FramesFactory.getITRF(IERSConventions.IERS_2010, true));
73      }
74  
75      @Test
76      public void testDelayAtIPP() {
77          GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
78          final double latitude  = FastMath.toRadians(30.0);
79          final double longitude = FastMath.toRadians(-130.0);
80          try {
81              final Method pathDelay = GlobalIonosphereMapModel.class.getDeclaredMethod("pathDelayAtIPP",
82                                                                                        AbsoluteDate.class,
83                                                                                        GeodeticPoint.class,
84                                                                                        Double.TYPE, Double.TYPE);
85              pathDelay.setAccessible(true);
86              final double delay = (Double) pathDelay.invoke(model,
87                                                             new AbsoluteDate(2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC()),
88                                                             new GeodeticPoint(latitude, longitude, 0.0),
89                                                             0.5 * FastMath.PI, PredefinedGnssSignal.G01.getFrequency());
90              Assertions.assertEquals(1.557, delay, epsilonDelay);
91          } catch (NoSuchMethodException | SecurityException | IllegalAccessException |
92                   IllegalArgumentException | InvocationTargetException e) {
93              Assertions.fail(e);
94          }
95      }
96  
97      @Test
98      public void testFieldDelayAtIPP() {
99          doTestFieldDelayAtIPP(Binary64Field.getInstance());
100     }
101 
102     private <T extends CalculusFieldElement<T>> void doTestFieldDelayAtIPP(final Field<T> field) {
103         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
104         final T zero = field.getZero();
105         final double latitude  = FastMath.toRadians(30.0);
106         final double longitude = FastMath.toRadians(-130.0);
107         try {
108             final Method pathDelay = GlobalIonosphereMapModel.class.getDeclaredMethod("pathDelayAtIPP",
109                                                                                       FieldAbsoluteDate.class,
110                                                                                       GeodeticPoint.class,
111                                                                                       CalculusFieldElement.class,
112                                                                                       Double.TYPE);
113             pathDelay.setAccessible(true);
114             @SuppressWarnings("unchecked")
115             final T delay = (T) pathDelay.invoke(model,
116                                                  new FieldAbsoluteDate<>(field, 2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC()),
117                                                  new GeodeticPoint(latitude, longitude, 0.0),
118                                                  zero.add(0.5 * FastMath.PI),
119                                                  PredefinedGnssSignal.G01.getFrequency());
120             Assertions.assertEquals(1.557, delay.getReal(), epsilonDelay);
121         } catch (NoSuchMethodException | SecurityException | IllegalAccessException |
122                         IllegalArgumentException | InvocationTargetException e) {
123             Assertions.fail(e);
124         }
125     }
126 
127     @Test
128     public void testSpacecraftState() {
129         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
130         double a = 7187990.1979844316;
131         double e = 0.5e-4;
132         double i = FastMath.toRadians(98.01);
133         double omega = FastMath.toRadians(131.88);
134         double OMEGA = FastMath.toRadians(252.24);
135         double lv = FastMath.toRadians(250.00);
136 
137         AbsoluteDate date = new AbsoluteDate(2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC());
138         final SpacecraftState state =
139                         new SpacecraftState(new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
140                                                                FramesFactory.getEME2000(), date,
141                                                                Constants.EIGEN5C_EARTH_MU));
142         final TopocentricFrame topo = new TopocentricFrame(earth,
143                                                            new GeodeticPoint(FastMath.toRadians(20.5236),
144                                                                              FastMath.toRadians(85.7881),
145                                                                              36.0),
146                                                            "Cuttack");
147         final double delay = model.pathDelay(state, topo, PredefinedGnssSignal.G01.getFrequency(), null);
148         Assertions.assertEquals(2.810, delay, epsilonDelay);
149 
150         // the delay at station longitude is different, due to IPP
151         try {
152             final Method pathDelay = GlobalIonosphereMapModel.class.getDeclaredMethod("pathDelayAtIPP",
153                                                                                       AbsoluteDate.class,
154                                                                                       GeodeticPoint.class,
155                                                                                       Double.TYPE, Double.TYPE);
156             pathDelay.setAccessible(true);
157             final double delayIPP = (Double) pathDelay.invoke(model, date, topo.getPoint(),
158                                                               0.5 * FastMath.PI,
159                                                               PredefinedGnssSignal.G01.getFrequency());
160             Assertions.assertEquals(2.173, delayIPP, epsilonDelay);
161         } catch (NoSuchMethodException | SecurityException | IllegalAccessException |
162                  IllegalArgumentException | InvocationTargetException ex) {
163             Assertions.fail(ex);
164         }
165     }
166 
167     @Test
168     public void testAboveIono() {
169         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
170         double a = 7187990.1979844316;
171         double e = 0.5e-4;
172         double i = FastMath.toRadians(98.01);
173         double omega = FastMath.toRadians(131.88);
174         double OMEGA = FastMath.toRadians(252.24);
175         double lv = FastMath.toRadians(250.00);
176 
177         AbsoluteDate date = new AbsoluteDate(2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC());
178         final SpacecraftState state =
179                         new SpacecraftState(new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
180                                                                FramesFactory.getEME2000(), date,
181                                                                Constants.EIGEN5C_EARTH_MU));
182         final TopocentricFrame topo = new TopocentricFrame(earth,
183                                                            new GeodeticPoint(FastMath.toRadians(20.5236),
184                                                                              FastMath.toRadians(85.7881),
185                                                                              650000.0),
186                                                            "very-high");
187         final double delay = model.pathDelay(state, topo, PredefinedGnssSignal.G01.getFrequency(), null);
188         Assertions.assertEquals(0.0, delay, epsilonDelay);
189 
190     }
191 
192     @Test
193     public void testSpacecraftStateField() {
194         doTestSpacecraftStateField(Binary64Field.getInstance());
195     }
196 
197     private <T extends CalculusFieldElement<T>> void doTestSpacecraftStateField(final Field<T> field) {
198         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
199         T a = field.getZero().newInstance(7187990.1979844316);
200         T e = field.getZero().newInstance(0.5e-4);
201         T i = FastMath.toRadians(field.getZero().newInstance(98.01));
202         T omega = FastMath.toRadians(field.getZero().newInstance(131.88));
203         T OMEGA = FastMath.toRadians(field.getZero().newInstance(252.24));
204         T lv = FastMath.toRadians(field.getZero().newInstance(250.00));
205 
206         FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field,
207                                                             new AbsoluteDate(2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC()));
208         final FieldSpacecraftState<T> state =
209                         new FieldSpacecraftState<>(new FieldKeplerianOrbit<>(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
210                                                                              FramesFactory.getEME2000(), date,
211                                                                              field.getZero().newInstance(Constants.EIGEN5C_EARTH_MU)));
212         final TopocentricFrame topo = new TopocentricFrame(earth,
213                                                            new GeodeticPoint(FastMath.toRadians(20.5236),
214                                                                              FastMath.toRadians(85.7881),
215                                                                              36.0),
216                                                            "Cuttack");
217         final T delay = model.pathDelay(state, topo, PredefinedGnssSignal.G01.getFrequency(), null);
218         Assertions.assertEquals(2.810, delay.getReal(), epsilonDelay);
219 
220         // the delay at station longitude is different, due to IPP
221         try {
222             final Method pathDelay = GlobalIonosphereMapModel.class.getDeclaredMethod("pathDelayAtIPP",
223                                                                                       FieldAbsoluteDate.class,
224                                                                                       GeodeticPoint.class,
225                                                                                       CalculusFieldElement.class,
226                                                                                       Double.TYPE);
227             pathDelay.setAccessible(true);
228             @SuppressWarnings("unchecked")
229             final T delayIPP = (T) pathDelay.invoke(model, date, topo.getPoint(),
230                                                     field.getZero().newInstance(0.5 * FastMath.PI),
231                                                     PredefinedGnssSignal.G01.getFrequency());
232             Assertions.assertEquals(2.173, delayIPP.getReal(), epsilonDelay);
233         } catch (NoSuchMethodException | SecurityException | IllegalAccessException |
234                  IllegalArgumentException | InvocationTargetException ex) {
235             Assertions.fail(ex);
236         }
237     }
238 
239     @Test
240     public void testAboveIonoField() {
241         doTestAboveIonoField(Binary64Field.getInstance());
242     }
243 
244     private <T extends CalculusFieldElement<T>> void doTestAboveIonoField(final Field<T> field) {
245         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
246         T a = field.getZero().newInstance(7187990.1979844316);
247         T e = field.getZero().newInstance(0.5e-4);
248         T i = FastMath.toRadians(field.getZero().newInstance(98.01));
249         T omega = FastMath.toRadians(field.getZero().newInstance(131.88));
250         T OMEGA = FastMath.toRadians(field.getZero().newInstance(252.24));
251         T lv = FastMath.toRadians(field.getZero().newInstance(250.00));
252 
253         FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field,
254                                                             new AbsoluteDate(2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC()));
255         final FieldSpacecraftState<T> state =
256                         new FieldSpacecraftState<>(new FieldKeplerianOrbit<>(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
257                                                                              FramesFactory.getEME2000(), date,
258                                                                              field.getZero().newInstance(Constants.EIGEN5C_EARTH_MU)));
259         final TopocentricFrame topo = new TopocentricFrame(earth,
260                                                            new GeodeticPoint(FastMath.toRadians(20.5236),
261                                                                              FastMath.toRadians(85.7881),
262                                                                              650000.0),
263                                                            "very-high");
264         final T delay = model.pathDelay(state, topo, PredefinedGnssSignal.G01.getFrequency(), null);
265         Assertions.assertEquals(0.0, delay.getReal(), epsilonDelay);
266 
267     }
268 
269     @Test
270     public void testParser() throws URISyntaxException {
271 
272         Utils.setDataRoot("regular-data");
273         URL url1 = DirectoryCrawlerTest.class.getClassLoader().getResource("ionex/split-1.19i");
274         DataSource ds1 = new DataSource(Paths.get(url1.toURI()).toString());
275         URL url2 = DirectoryCrawlerTest.class.getClassLoader().getResource("ionex/split-2.19i");
276         DataSource ds2 = new DataSource(Paths.get(url2.toURI()).toString());
277         URL url3 = DirectoryCrawlerTest.class.getClassLoader().getResource("ionex/split-3.19i");
278         DataSource ds3 = new DataSource(Paths.get(url3.toURI()).toString());
279         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel(TimeScalesFactory.getUTC(), ds1, ds2, ds3);
280 
281         // Commons parameters
282         final double latitude = FastMath.toRadians(45.0);
283 
284         double longitude1;
285         double longitude2;
286 
287         try {
288             final Method pathDelay = GlobalIonosphereMapModel.class.getDeclaredMethod("pathDelayAtIPP",
289                                                                                       AbsoluteDate.class,
290                                                                                       GeodeticPoint.class,
291                                                                                       Double.TYPE, Double.TYPE);
292             pathDelay.setAccessible(true);
293 
294             // Test longitude = 181° and longitude = -179°
295             longitude1 = FastMath.toRadians(181.0);
296             longitude2 = FastMath.toRadians(-179.0);
297             AbsoluteDate date1 = new AbsoluteDate(2019, 1, 15, 1, 0, 0.0, TimeScalesFactory.getUTC());
298             Assertions.assertEquals((Double) pathDelay.invoke(model, date1, new GeodeticPoint(latitude, longitude1, 0.0),
299                                                               0.01, PredefinedGnssSignal.G01.getFrequency()),
300                                     ((Double) pathDelay.invoke(model, date1, new GeodeticPoint(latitude, longitude2, 0.0),
301                                                                0.01, PredefinedGnssSignal.G01.getFrequency())),
302                                     epsilonParser);
303 
304             // Test longitude = 180° and longitude = -180°
305             AbsoluteDate date2 = new AbsoluteDate(2019, 1, 15, 3, 0, 0.0, TimeScalesFactory.getUTC());
306             longitude1 = FastMath.toRadians(180.0);
307             longitude2 = FastMath.toRadians(-180.0);
308 
309             Assertions.assertEquals(((Double) pathDelay.invoke(model, date2, new GeodeticPoint(latitude, longitude1, 0.0),
310                                                                0.01, PredefinedGnssSignal.G01.getFrequency())),
311                                     ((Double) pathDelay.invoke(model, date2, new GeodeticPoint(latitude, longitude2, 0.0),
312                                                                0.01, PredefinedGnssSignal.G01.getFrequency())),
313                                     epsilonParser);
314 
315             // Test longitude = 0° and longitude = 360°
316             AbsoluteDate date3 = new AbsoluteDate(2019, 1, 15, 5, 0, 0.0, TimeScalesFactory.getUTC());
317             longitude1 =  FastMath.toRadians(0.);
318             longitude2 =  FastMath.toRadians(360.0);
319 
320             Assertions.assertEquals(((Double) pathDelay.invoke(model, date3, new GeodeticPoint(latitude, longitude1, 0.0),
321                                                                0.01, PredefinedGnssSignal.G01.getFrequency())),
322                                     ((Double) pathDelay.invoke(model, date3, new GeodeticPoint(latitude, longitude2, 0.0),
323                                                                0.01, PredefinedGnssSignal.G01.getFrequency())),
324                                     epsilonParser);
325 
326         } catch (NoSuchMethodException | SecurityException | IllegalAccessException |
327                         IllegalArgumentException | InvocationTargetException e) {
328             Assertions.fail(e);
329         }
330     }
331 
332     @Test
333     public void testCorruptedFileBadData() {
334         final String fileName = "corrupted-bad-data-gpsg0150.19i";
335 
336         try {
337             new GlobalIonosphereMapModel(fileName);
338             Assertions.fail("An exception should have been thrown");
339 
340         } catch (OrekitException oe) {
341             Assertions.assertEquals(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE, oe.getSpecifier());
342         }
343     }
344 
345     @Test
346     public void testEarlierDate() {
347         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
348         final double latitude  = FastMath.toRadians(60.0);
349         final double longitude = FastMath.toRadians(-130.0);
350         final GeodeticPoint point = new GeodeticPoint(latitude, longitude, 0.0);
351 
352         try {
353             model.pathDelay(state, new TopocentricFrame(earth, point, null),
354                             PredefinedGnssSignal.G01.getFrequency(),
355                             model.getParameters(new AbsoluteDate()));
356             Assertions.fail("An exception should have been thrown");
357         } catch (OrekitException oe) {
358             Assertions.assertEquals(OrekitMessages.NO_TEC_DATA_IN_FILES_FOR_DATE, oe.getSpecifier());
359         }
360     }
361 
362     @Test
363     public void testFieldEarlierDate() {
364         doTestFieldEarlierDate(Binary64Field.getInstance());
365     }
366 
367     private <T extends CalculusFieldElement<T>> void doTestFieldEarlierDate(final Field<T> field) {
368         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
369         final double latitude  = FastMath.toRadians(60.0);
370         final double longitude = FastMath.toRadians(-130.0);
371         final GeodeticPoint point = new GeodeticPoint(latitude, longitude, 0.0);
372 
373         try {
374             model.pathDelay(new FieldSpacecraftState<>(field, state), new TopocentricFrame(earth, point, null),
375                             PredefinedGnssSignal.G01.getFrequency(),
376                             model.getParameters(field, new FieldAbsoluteDate<>(field)));
377             Assertions.fail("An exception should have been thrown");
378         } catch (OrekitException oe) {
379             Assertions.assertEquals(OrekitMessages.NO_TEC_DATA_IN_FILES_FOR_DATE, oe.getSpecifier());
380         }
381     }
382 
383     @Test
384     public void testLaterDate() {
385         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
386         final double latitude  = FastMath.toRadians(60.0);
387         final double longitude = FastMath.toRadians(-130.0);
388         final GeodeticPoint point = new GeodeticPoint(latitude, longitude, 0.0);
389 
390         try {
391             model.pathDelay(state, new TopocentricFrame(earth, point, null),
392                             PredefinedGnssSignal.G01.getFrequency(),
393                             model.getParameters(new AbsoluteDate()));
394             Assertions.fail("An exception should have been thrown");
395         } catch (OrekitException oe) {
396             Assertions.assertEquals(OrekitMessages.NO_TEC_DATA_IN_FILES_FOR_DATE, oe.getSpecifier());
397         }
398 
399     }
400 
401     @Test
402     public void testFieldLaterDate() {
403         doTestFieldLaterDate(Binary64Field.getInstance());
404     }
405 
406     private <T extends CalculusFieldElement<T>> void doTestFieldLaterDate(final Field<T> field) {
407         GlobalIonosphereMapModel model = new GlobalIonosphereMapModel("gpsg0150.19i");
408         final double latitude  = FastMath.toRadians(60.0);
409         final double longitude = FastMath.toRadians(-130.0);
410         final GeodeticPoint point = new GeodeticPoint(latitude, longitude, 0.0);
411         try {
412             model.pathDelay(new FieldSpacecraftState<>(field, state), new TopocentricFrame(earth, point, null),
413                             PredefinedGnssSignal.G01.getFrequency(),
414                             model.getParameters(field, new FieldAbsoluteDate<>(field)));
415             Assertions.fail("An exception should have been thrown");
416         } catch (OrekitException oe) {
417             Assertions.assertEquals(OrekitMessages.NO_TEC_DATA_IN_FILES_FOR_DATE, oe.getSpecifier());
418         }
419 
420     }
421 
422     @Test
423     /**
424      * The goal of this test is to verify if an OrekitException is thrown when latitude or longitude
425      * boundaries are not present in the header section of the Global Ionosphere Map.
426      */
427     public void testIssue621() {
428         final String fileName  = "missing-lat-lon-header-gpsg0150.19i";
429 
430         try {
431             new GlobalIonosphereMapModel(fileName);
432             Assertions.fail("An exception should have been thrown");
433 
434         } catch (OrekitException oe) {
435             Assertions.assertEquals(OrekitMessages.NO_LATITUDE_LONGITUDE_BOUNDARIES_IN_IONEX_HEADER, oe.getSpecifier());
436         }
437     }
438 
439     @Test
440     public void testJammedFields() {
441         Assertions.assertNotNull(new GlobalIonosphereMapModel("jammed-fields.14i"));
442     }
443 
444     @Test
445     public void testMissingEpochInHeader() {
446         final String fileName  = "missing-epoch-header-gpsg0150.19i";
447 
448         try {
449             new GlobalIonosphereMapModel(fileName);
450             Assertions.fail("An exception should have been thrown");
451 
452         } catch (OrekitException oe) {
453             Assertions.assertEquals(OrekitMessages.NO_EPOCH_IN_IONEX_HEADER, oe.getSpecifier());
454         }
455 
456     }
457 
458 }