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.nequick;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.Field;
21  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
22  import org.hipparchus.geometry.euclidean.threed.Vector3D;
23  import org.hipparchus.util.Binary64Field;
24  import org.hipparchus.util.FastMath;
25  import org.junit.jupiter.api.Assertions;
26  import org.junit.jupiter.api.BeforeEach;
27  import org.junit.jupiter.api.Test;
28  import org.orekit.Utils;
29  import org.orekit.bodies.FieldGeodeticPoint;
30  import org.orekit.bodies.GeodeticPoint;
31  import org.orekit.bodies.OneAxisEllipsoid;
32  import org.orekit.frames.FramesFactory;
33  import org.orekit.frames.TopocentricFrame;
34  import org.orekit.gnss.PredefinedGnssSignal;
35  import org.orekit.orbits.CartesianOrbit;
36  import org.orekit.orbits.FieldCartesianOrbit;
37  import org.orekit.orbits.FieldOrbit;
38  import org.orekit.orbits.Orbit;
39  import org.orekit.propagation.FieldSpacecraftState;
40  import org.orekit.propagation.SpacecraftState;
41  import org.orekit.time.AbsoluteDate;
42  import org.orekit.time.FieldAbsoluteDate;
43  import org.orekit.time.TimeScalesFactory;
44  import org.orekit.utils.Constants;
45  import org.orekit.utils.FieldPVCoordinates;
46  import org.orekit.utils.IERSConventions;
47  import org.orekit.utils.PVCoordinates;
48  
49  /**
50   * Reference values for the tests are from : "European Union (2016). European GNSS (Galileo)
51   * Open Service-Ionospheric Correction Algorithm for Galileo Single Frequency Users. 1.2."
52   */
53  public class NeQuickGalileoTest {
54  
55      private double[] medium;
56      private double[] high;
57  
58      @BeforeEach
59      public void setUp() {
60          Utils.setDataRoot("regular-data");
61          high = new double[] {
62              236.831641, -0.39362878, 0.00402826613
63          };
64          medium = new double[] {
65              121.129893, 0.351254133, 0.0134635348
66          };
67  
68      }
69  
70      @Test
71      public void testHighSolarActivityGalileo() {
72  
73          // Model
74          final NeQuickGalileo model = new NeQuickGalileo(high);
75  
76          // Getters
77          Assertions.assertEquals(236.831641,    model.getAlpha()[0], 1.0e-6);
78          Assertions.assertEquals(-0.39362878,   model.getAlpha()[1], 1.0e-8);
79          Assertions.assertEquals(0.00402826613, model.getAlpha()[2], 1.0e-11);
80          Assertions.assertEquals("UTC", model.getUtc().getName());
81  
82          // Geodetic points
83          final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(82.49), FastMath.toRadians(297.66), 78.11);
84          final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(54.29), FastMath.toRadians(8.23), 20281546.18);
85  
86          // Date
87          final AbsoluteDate date = new AbsoluteDate(2018, 4, 2, 0, 0, 0, TimeScalesFactory.getUTC());
88  
89          // STEC
90          final double stec = model.stec(date, recP, satP);
91          Assertions.assertEquals(20.319, stec, 1.0e-3);
92      }
93  
94      @Test
95      public void testFieldHighSolarActivityGalileo() {
96          doTestFieldHighSolarActivityGalileo(Binary64Field.getInstance());
97      }
98  
99      private <T extends CalculusFieldElement<T>> void doTestFieldHighSolarActivityGalileo(final Field<T> field) {
100 
101         // Zero
102         final T zero = field.getZero();
103 
104         // Model
105         final NeQuickGalileo model = new NeQuickGalileo(high);
106 
107         // Geodetic points
108         final FieldGeodeticPoint<T> recP = new FieldGeodeticPoint<>(zero.newInstance(FastMath.toRadians(82.49)),
109                                                                     zero.newInstance(FastMath.toRadians(297.66)),
110                                                                     zero.newInstance(78.11));
111         final FieldGeodeticPoint<T> satP = new FieldGeodeticPoint<>(zero.newInstance(FastMath.toRadians(54.29)),
112                                                                     zero.newInstance(FastMath.toRadians(8.23)),
113                                                                     zero.newInstance(20281546.18));
114 
115         // Date
116         final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2018, 4, 2, 0, 0, 0, TimeScalesFactory.getUTC());
117 
118         // STEC
119         final T stec = model.stec(date, recP, satP);
120         Assertions.assertEquals(20.319, stec.getReal(), 1.0e-3);
121     }
122 
123     @Test
124     public void testMediumSolarActivityGalileo() {
125 
126         // Model
127         final NeQuickGalileo model = new NeQuickGalileo(medium);
128 
129         // Geodetic points
130         final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(-31.80), FastMath.toRadians(115.89), 12.78);
131         final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(-14.31), FastMath.toRadians(124.09), 20100697.90);
132 
133         // Date
134         final AbsoluteDate date = new AbsoluteDate(2018, 1, 2, 16, 0, 0, TimeScalesFactory.getUTC());
135 
136         // STEC
137         final double stec = model.stec(date, recP, satP);
138         Assertions.assertEquals(20.481, stec, 1.0e-3);
139     }
140 
141     @Test
142     public void testFieldMediumSolarActivityGalileo() {
143         doTestFieldMediumSolarActivityGalileo(Binary64Field.getInstance());
144     }
145 
146     private <T extends CalculusFieldElement<T>> void doTestFieldMediumSolarActivityGalileo(final Field<T> field) {
147 
148         // Zero
149         final T zero = field.getZero();
150 
151         // Model
152         final NeQuickGalileo model = new NeQuickGalileo(medium);
153 
154         // Geodetic points
155         final FieldGeodeticPoint<T> recP = new FieldGeodeticPoint<>(zero.newInstance(FastMath.toRadians(-31.80)),
156                                                                     zero.newInstance(FastMath.toRadians(115.89)),
157                                                                     zero.newInstance(12.78));
158         final FieldGeodeticPoint<T> satP = new FieldGeodeticPoint<>(zero.newInstance(FastMath.toRadians(-14.31)),
159                                                                     zero.newInstance(FastMath.toRadians(124.09)),
160                                                                     zero.newInstance(20100697.90));
161 
162         // Date
163         final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2018, 1, 2, 16, 0, 0, TimeScalesFactory.getUTC());
164 
165         // STEC
166         final T stec = model.stec(date, recP, satP);
167         Assertions.assertEquals(20.481, stec.getReal(), 1.0e-3);
168     }
169 
170     @Test
171     public void testDelay() {
172 
173         // Model
174         final NeQuickGalileo model = new NeQuickGalileo(medium);
175 
176         // Geodetic points
177         final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(-31.80), FastMath.toRadians(115.89), 12.78);
178         final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(-14.31), FastMath.toRadians(124.09), 20100697.90);
179 
180         // Date
181         final AbsoluteDate date = new AbsoluteDate(2018, 4, 2, 16, 0, 0, TimeScalesFactory.getUTC());
182 
183         // Earth
184         final OneAxisEllipsoid ellipsoid = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
185                                                                 Constants.WGS84_EARTH_FLATTENING,
186                                                                 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
187         // Satellite position
188         final Vector3D satPosInITRF    = ellipsoid.transform(satP);
189         final Vector3D satPosInEME2000 = ellipsoid.getBodyFrame().getStaticTransformTo(FramesFactory.getEME2000(), date).transformPosition(satPosInITRF);
190 
191         // Spacecraft state
192         final PVCoordinates   pv      = new PVCoordinates(satPosInEME2000, new Vector3D(1.0, 1.0, 1.0));
193         final Orbit           orbit   = new CartesianOrbit(pv, FramesFactory.getEME2000(), date, Constants.WGS84_EARTH_MU);
194         final SpacecraftState state   = new SpacecraftState(orbit);
195 
196         final double delay = model.pathDelay(state, new TopocentricFrame(ellipsoid, recP, null),
197                                              PredefinedGnssSignal.G01.getFrequency(), model.getParameters());
198        
199         // Verify
200         Assertions.assertEquals(1.13, delay, 0.01);
201     }
202 
203     @Test
204     public void testFieldDelay() {
205         doTestFieldDelay(Binary64Field.getInstance());
206     }
207 
208     private <T extends CalculusFieldElement<T>> void doTestFieldDelay(final Field<T> field) {
209 
210         // Zero and One
211         final T zero = field.getZero();
212         final T one  = field.getOne();
213 
214         // Model
215         final NeQuickGalileo model = new NeQuickGalileo(medium);
216 
217         // Geodetic points
218         final double recLat = FastMath.toRadians(-31.80);
219         final double recLon = FastMath.toRadians(115.89);
220         final double recAlt = 12.78;
221         final GeodeticPoint         recP = new GeodeticPoint(recLat, recLon, recAlt);
222         final FieldGeodeticPoint<T> satP = new FieldGeodeticPoint<>(zero.newInstance(FastMath.toRadians(-14.31)),
223                         zero.newInstance(FastMath.toRadians(124.09)), zero.newInstance(20100697.90));
224 
225         // Date
226         final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2018, 4, 2, 16, 0, 0, TimeScalesFactory.getUTC());
227 
228         // Earth
229         final OneAxisEllipsoid ellipsoid = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
230                                                                 Constants.WGS84_EARTH_FLATTENING,
231                                                                 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
232         // Satellite position
233         final FieldVector3D<T> satPosInITRF    = ellipsoid.transform(satP);
234         final FieldVector3D<T> satPosInEME2000 = ellipsoid.getBodyFrame().getStaticTransformTo(FramesFactory.getEME2000(), date).transformPosition(satPosInITRF);
235 
236         // Spacecraft state
237         final FieldPVCoordinates<T>   pv      = new FieldPVCoordinates<>(satPosInEME2000, new FieldVector3D<>(one, one, one));
238         final FieldOrbit<T>           orbit   = new FieldCartesianOrbit<>(pv, FramesFactory.getEME2000(), date, zero.newInstance(Constants.WGS84_EARTH_MU));
239         final FieldSpacecraftState<T> state   = new FieldSpacecraftState<>(orbit);
240 
241         final T delay = model.pathDelay(state, new TopocentricFrame(ellipsoid, recP, null),
242                                         PredefinedGnssSignal.G01.getFrequency(), model.getParameters(field));
243        
244         // Verify
245         Assertions.assertEquals(1.13, delay.getReal(), 0.01);
246     }
247 
248     @Test
249     public void testAntiMeridian() {
250 
251         // Model
252         final NeQuickGalileo model = new NeQuickGalileo(medium);
253 
254         // Date
255         final AbsoluteDate date = new AbsoluteDate(2018,  11,  2, 16, 0, 0, TimeScalesFactory.getUTC());
256 
257         // Geodetic points
258         final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(-31.80), FastMath.toRadians(-179.99), 12.78);
259         final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(-14.31), FastMath.toRadians(-177.43), 20100697.90);
260         double stec = model.stec(date, recP, satP);
261         Assertions.assertEquals(20.471, stec, 0.001);
262 
263     }
264 
265     @Test
266     public void testFieldAntiMeridian() {
267         doTestFieldAntiMeridian(Binary64Field.getInstance());
268     }
269 
270     private <T extends CalculusFieldElement<T>> void doTestFieldAntiMeridian(final Field<T> field) {
271 
272         final T zero = field.getZero();
273 
274         // Model
275         final NeQuickGalileo model = new NeQuickGalileo(medium);
276 
277         // Date
278         final FieldAbsoluteDate<T> date =
279                         new FieldAbsoluteDate<>(field,
280                                                 new AbsoluteDate(2018,  11,  2, 16, 0, 0, TimeScalesFactory.getUTC()));
281 
282         // Geodetic points
283         final FieldGeodeticPoint<T> recP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(-31.80)),
284                                                                     FastMath.toRadians(zero.newInstance(-179.99)),
285                                                                     zero.newInstance(12.78));
286         final FieldGeodeticPoint<T> satP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(-14.31)),
287                                                                     FastMath.toRadians(zero.newInstance(-177.43)),
288                                                                     zero.newInstance(20100697.90));
289         T stec = model.stec(date, recP, satP);
290         Assertions.assertEquals(20.471, stec.getReal(), 0.001);
291 
292     }
293 
294     @Test
295     public void testZenith() {
296 
297         // Model
298         final NeQuickGalileo model = new NeQuickGalileo(medium);
299 
300         // Date
301         final AbsoluteDate date = new AbsoluteDate(2018,  4,  2, 16, 0, 0, TimeScalesFactory.getUTC());
302 
303         // Geodetic points
304         final GeodeticPoint recP = new GeodeticPoint(FastMath.toRadians(51.678), FastMath.toRadians(-9.448), 0.0);
305         final GeodeticPoint satP = new GeodeticPoint(FastMath.toRadians(51.678), FastMath.toRadians(-9.448), 600000.0);
306         double stec = model.stec(date, recP, satP);
307         Assertions.assertEquals(26.346, stec, 0.001);
308 
309     }
310 
311     @Test
312     public void testFieldZenith() {
313         doTestFieldZenith(Binary64Field.getInstance());
314     }
315 
316     private <T extends CalculusFieldElement<T>> void doTestFieldZenith(final Field<T> field) {
317 
318         final T zero = field.getZero();
319 
320         // Model
321         final NeQuickGalileo model = new NeQuickGalileo(medium);
322 
323         // Date
324         final FieldAbsoluteDate<T> date =
325                 new FieldAbsoluteDate<>(field,
326                                         new AbsoluteDate(2018,  4,  2, 16, 0, 0, TimeScalesFactory.getUTC()));
327 
328         // Geodetic points
329         final FieldGeodeticPoint<T> recP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(51.678)),
330                                                                     FastMath.toRadians(zero.newInstance(-9.448)),
331                                                                     zero.newInstance(0.0));
332         final FieldGeodeticPoint<T> satP = new FieldGeodeticPoint<>(FastMath.toRadians(zero.newInstance(51.678)),
333                                                                     FastMath.toRadians(zero.newInstance(-9.448)),
334                                                                     zero.newInstance(600000.0));
335         T stec = model.stec(date, recP, satP);
336         Assertions.assertEquals(26.346, stec.getReal(), 0.001);
337 
338     }
339 
340 }