1   /* Copyright 2002-2022 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.io.ByteArrayInputStream;
20  import java.io.ByteArrayOutputStream;
21  import java.io.IOException;
22  import java.io.ObjectInputStream;
23  import java.io.ObjectOutputStream;
24  
25  import org.hipparchus.Field;
26  import org.hipparchus.CalculusFieldElement;
27  import org.hipparchus.util.Decimal64Field;
28  import org.hipparchus.util.FastMath;
29  import org.junit.Assert;
30  import org.junit.Before;
31  import org.junit.Test;
32  import org.orekit.Utils;
33  import org.orekit.bodies.GeodeticPoint;
34  import org.orekit.bodies.OneAxisEllipsoid;
35  import org.orekit.errors.OrekitException;
36  import org.orekit.errors.OrekitMessages;
37  import org.orekit.frames.FramesFactory;
38  import org.orekit.frames.TopocentricFrame;
39  import org.orekit.gnss.Frequency;
40  import org.orekit.orbits.KeplerianOrbit;
41  import org.orekit.orbits.Orbit;
42  import org.orekit.orbits.PositionAngle;
43  import org.orekit.propagation.FieldSpacecraftState;
44  import org.orekit.propagation.SpacecraftState;
45  import org.orekit.time.AbsoluteDate;
46  import org.orekit.time.FieldAbsoluteDate;
47  import org.orekit.time.TimeScalesFactory;
48  import org.orekit.utils.Constants;
49  import org.orekit.utils.IERSConventions;
50  
51  public class GlobalIonosphereMapModelTest {
52  
53      private static double epsilonParser = 1.0e-16;
54      private static double epsilonDelay  = 0.001;
55      private SpacecraftState state;
56      private OneAxisEllipsoid earth;
57      private GlobalIonosphereMapModel model;
58  
59      @Before
60      public void setUp() {
61          Utils.setDataRoot("regular-data:ionex");
62          model = new GlobalIonosphereMapModel("gpsg0150.19i");
63          final Orbit orbit = new KeplerianOrbit(24464560.0, 0.0, 1.122138, 1.10686, 1.00681,
64                                                 0.048363, PositionAngle.MEAN,
65                                                 FramesFactory.getEME2000(),
66                                                 new AbsoluteDate(2019, 1, 14, 23, 59, 59.0, TimeScalesFactory.getUTC()),
67                                                 Constants.WGS84_EARTH_MU);
68          state = new SpacecraftState(orbit);
69          earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
70                                       Constants.WGS84_EARTH_FLATTENING,
71                                       FramesFactory.getITRF(IERSConventions.IERS_2010, true));
72      }
73  
74      @Test
75      public void testTEC() {
76          final double latitude  = FastMath.toRadians(30.0);
77          final double longitude = FastMath.toRadians(-130.0); 
78          final double tec = model.getTEC(new AbsoluteDate(2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC()),
79                                          new GeodeticPoint(latitude, longitude, 0.0));
80          Assert.assertEquals(9.592, tec, epsilonDelay);
81      }
82  
83      @Test
84      public void testFieldTEC() {
85          doTestFieldTEC(Decimal64Field.getInstance());
86      }
87  
88      private <T extends CalculusFieldElement<T>> void doTestFieldTEC(final Field<T> field) {
89          final double latitude  = FastMath.toRadians(30.0);
90          final double longitude = FastMath.toRadians(-130.0);
91          final T tec = model.getTEC(new FieldAbsoluteDate<>(field, 2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC()),
92                                     new GeodeticPoint(latitude, longitude, 0.0));
93          Assert.assertEquals(9.592, tec.getReal(), epsilonDelay);
94      }
95  
96      @Test
97      public void testDelay() {
98          final double latitude  = FastMath.toRadians(30.0);
99          final double longitude = FastMath.toRadians(-130.0);
100         final double delay = model.pathDelay(new AbsoluteDate(2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC()),
101                                              new GeodeticPoint(latitude, longitude, 0.0),
102                                              0.5 * FastMath.PI, Frequency.G01.getMHzFrequency() * 1.0e6);
103         Assert.assertEquals(1.557, delay, epsilonDelay);
104     }
105 
106     @Test
107     public void testFieldDelay() {
108         doTestFieldDelay(Decimal64Field.getInstance());
109     }
110 
111     private <T extends CalculusFieldElement<T>> void doTestFieldDelay(final Field<T> field) {
112         final T zero = field.getZero();
113         final double latitude  = FastMath.toRadians(30.0);
114         final double longitude = FastMath.toRadians(-130.0);
115         final T delay = model.pathDelay(new FieldAbsoluteDate<>(field, 2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC()),
116                                         new GeodeticPoint(latitude, longitude, 0.0),
117                                         zero.add(0.5 * FastMath.PI),
118                                         Frequency.G01.getMHzFrequency() * 1.0e6);
119         Assert.assertEquals(1.557, delay.getReal(), epsilonDelay);
120     }
121 
122     @Test
123     public void testParser() {
124 
125         // Commons parameters
126         AbsoluteDate date = new AbsoluteDate(2019, 1, 15, 0, 0, 0.0, TimeScalesFactory.getUTC());
127         final double latitude = FastMath.toRadians(45.0);
128         
129         double longitude1;
130         double longitude2;
131 
132         // Test longitude = 181° and longitude = -179°
133         longitude1 = FastMath.toRadians(181.0);
134         longitude2 = FastMath.toRadians(-179.0);
135 
136         Assert.assertEquals(model.getTEC(date, new GeodeticPoint(latitude, longitude1, 0.0)), model.getTEC(date, new GeodeticPoint(latitude, longitude2, 0.0)), epsilonParser);
137 
138         // Test longitude = 180° and longitude = -180°
139         longitude1 = FastMath.toRadians(180.0);
140         longitude2 = FastMath.toRadians(-180.0);
141 
142         Assert.assertEquals(model.getTEC(date, new GeodeticPoint(latitude, longitude1, 0.0)), model.getTEC(date, new GeodeticPoint(latitude, longitude2, 0.0)), epsilonParser);
143 
144         // Test longitude = 0° and longitude = 360°
145         longitude1 =  FastMath.toRadians(0.);
146         longitude2 =  FastMath.toRadians(360.0);
147 
148         Assert.assertEquals(model.getTEC(date, new GeodeticPoint(latitude, longitude1, 0.0)), model.getTEC(date, new GeodeticPoint(latitude, longitude2, 0.0)), epsilonParser);
149 
150     }
151 
152     @Test
153     public void testCorruptedFileBadData() {
154         final String fileName = "corrupted-bad-data-gpsg0150.19i";
155         final double latitude  = FastMath.toRadians(30.0);
156         final double longitude = FastMath.toRadians(-130.0);
157           
158         try {
159             GlobalIonosphereMapModel corruptedModel = new GlobalIonosphereMapModel(fileName);
160             corruptedModel.getTEC(new AbsoluteDate(2019, 1, 15, 0, 0, 0.0, TimeScalesFactory.getUTC()),
161                          new GeodeticPoint(latitude, longitude, 0.0));
162             Assert.fail("An exception should have been thrown");
163             
164         } catch (OrekitException oe) {
165             Assert.assertEquals(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE, oe.getSpecifier());
166         }
167     }
168 
169     @Test
170     public void testEarlierDate() {
171         final double latitude  = FastMath.toRadians(60.0);
172         final double longitude = FastMath.toRadians(-130.0);
173         final GeodeticPoint point = new GeodeticPoint(latitude, longitude, 0.0);
174 
175         try {
176             model.pathDelay(state, new TopocentricFrame(earth, point, null),
177                             Frequency.G01.getMHzFrequency() * 1.0e6,
178                             model.getParameters());
179             Assert.fail("An exception should have been thrown");
180         } catch (OrekitException oe) {
181             Assert.assertEquals(OrekitMessages.NO_TEC_DATA_IN_FILE_FOR_DATE, oe.getSpecifier());
182         }
183     }
184 
185     @Test
186     public void testFieldEarlierDate() {
187         doTestFieldEarlierDate(Decimal64Field.getInstance());
188     }
189 
190     private <T extends CalculusFieldElement<T>> void doTestFieldEarlierDate(final Field<T> field) {
191         final double latitude  = FastMath.toRadians(60.0);
192         final double longitude = FastMath.toRadians(-130.0);
193         final GeodeticPoint point = new GeodeticPoint(latitude, longitude, 0.0);
194 
195         try {
196             model.pathDelay(new FieldSpacecraftState<>(field, state), new TopocentricFrame(earth, point, null),
197                             Frequency.G01.getMHzFrequency() * 1.0e6,
198                             model.getParameters(field));
199             Assert.fail("An exception should have been thrown");
200         } catch (OrekitException oe) {
201             Assert.assertEquals(OrekitMessages.NO_TEC_DATA_IN_FILE_FOR_DATE, oe.getSpecifier());
202         }
203     }
204 
205     @Test
206     public void testLaterDate() {
207         final double latitude  = FastMath.toRadians(60.0);
208         final double longitude = FastMath.toRadians(-130.0);
209         final GeodeticPoint point = new GeodeticPoint(latitude, longitude, 0.0);
210 
211         try {
212             model.pathDelay(state, new TopocentricFrame(earth, point, null),
213                             Frequency.G01.getMHzFrequency() * 1.0e6,
214                             model.getParameters());
215             Assert.fail("An exception should have been thrown");
216         } catch (OrekitException oe) {
217             Assert.assertEquals(OrekitMessages.NO_TEC_DATA_IN_FILE_FOR_DATE, oe.getSpecifier());
218         }
219 
220     }
221 
222     @Test
223     public void testFieldLaterDate() {
224         doTestFieldLaterDate(Decimal64Field.getInstance());
225     }
226 
227     private <T extends CalculusFieldElement<T>> void doTestFieldLaterDate(final Field<T> field) {
228         final double latitude  = FastMath.toRadians(60.0);
229         final double longitude = FastMath.toRadians(-130.0);
230         final GeodeticPoint point = new GeodeticPoint(latitude, longitude, 0.0);
231         try {
232             model.pathDelay(new FieldSpacecraftState<>(field, state), new TopocentricFrame(earth, point, null),
233                             Frequency.G01.getMHzFrequency() * 1.0e6,
234                             model.getParameters(field));
235             Assert.fail("An exception should have been thrown");
236         } catch (OrekitException oe) {
237             Assert.assertEquals(OrekitMessages.NO_TEC_DATA_IN_FILE_FOR_DATE, oe.getSpecifier());
238         }
239 
240     }
241 
242     @Test
243     /**
244      * The goal of this test is to verify if an OrekitException is thrown when latitude or longitude
245      * bondaries are not present in the header section of the Global Ionosphere Map.
246      */
247     public void testIssue621() {
248         final String fileName  = "missing-lat-lon-header-gpsg0150.19i";
249         final double latitude  = FastMath.toRadians(30.0);
250         final double longitude = FastMath.toRadians(-130.0);
251 
252         try {
253             GlobalIonosphereMapModel corruptedModel = new GlobalIonosphereMapModel(fileName);
254             corruptedModel.getTEC(new AbsoluteDate(2019, 1, 15, 0, 0, 0.0, TimeScalesFactory.getUTC()),
255                          new GeodeticPoint(latitude, longitude, 0.0));
256             Assert.fail("An exception should have been thrown");
257             
258         } catch (OrekitException oe) {
259             Assert.assertEquals(OrekitMessages.NO_LATITUDE_LONGITUDE_BONDARIES_IN_IONEX_HEADER, oe.getSpecifier());
260         }
261     }
262 
263     @Test
264     public void testMissingEpochInHeader() {
265         final String fileName  = "missing-epoch-header-gpsg0150.19i";
266         final double latitude  = FastMath.toRadians(30.0);
267         final double longitude = FastMath.toRadians(-130.0);
268 
269         try {
270             GlobalIonosphereMapModel corruptedModel = new GlobalIonosphereMapModel(fileName);
271             corruptedModel.getTEC(new AbsoluteDate(2019, 1, 15, 0, 0, 0.0, TimeScalesFactory.getUTC()),
272                          new GeodeticPoint(latitude, longitude, 0.0));
273             Assert.fail("An exception should have been thrown");
274             
275         } catch (OrekitException oe) {
276             Assert.assertEquals(OrekitMessages.NO_EPOCH_IN_IONEX_HEADER, oe.getSpecifier());
277         }
278 
279     }
280 
281     @Test
282     public void testSerialization() throws IOException, ClassNotFoundException {
283         GlobalIonosphereMapModel original = new GlobalIonosphereMapModel("gpsg0150.19i");
284 
285         ByteArrayOutputStream bos = new ByteArrayOutputStream();
286         ObjectOutputStream    oos = new ObjectOutputStream(bos);
287         oos.writeObject(original);
288 
289         Assert.assertTrue(bos.size() > 150);
290         Assert.assertTrue(bos.size() < 200);
291 
292         ByteArrayInputStream  bis = new ByteArrayInputStream(bos.toByteArray());
293         ObjectInputStream     ois = new ObjectInputStream(bis);
294         GlobalIonosphereMapModel deserialized = (GlobalIonosphereMapModel) ois.readObject();
295 
296         AbsoluteDate date = new AbsoluteDate(2019, 1, 15, 3, 43, 12.0, TimeScalesFactory.getUTC());
297         GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(30.0), FastMath.toRadians(-130.0), 0.0);
298         Assert.assertEquals(original.getTEC(date, point), deserialized.getTEC(date, point), 1.0e-12);
299     }
300 
301 }