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 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.hipparchus.util.Precision;
26  import org.junit.jupiter.api.AfterEach;
27  import org.junit.jupiter.api.Assertions;
28  import org.junit.jupiter.api.BeforeEach;
29  import org.junit.jupiter.api.Test;
30  import org.orekit.Utils;
31  import org.orekit.bodies.BodyShape;
32  import org.orekit.bodies.FieldGeodeticPoint;
33  import org.orekit.bodies.GeodeticPoint;
34  import org.orekit.bodies.OneAxisEllipsoid;
35  import org.orekit.errors.OrekitException;
36  import org.orekit.frames.Frame;
37  import org.orekit.frames.FramesFactory;
38  import org.orekit.frames.TopocentricFrame;
39  import org.orekit.orbits.CartesianOrbit;
40  import org.orekit.orbits.FieldCartesianOrbit;
41  import org.orekit.orbits.FieldOrbit;
42  import org.orekit.orbits.Orbit;
43  import org.orekit.propagation.FieldSpacecraftState;
44  import org.orekit.propagation.SpacecraftState;
45  import org.orekit.time.AbsoluteDate;
46  import org.orekit.time.DateTimeComponents;
47  import org.orekit.time.FieldAbsoluteDate;
48  import org.orekit.time.TimeScalesFactory;
49  import org.orekit.time.UTCScale;
50  import org.orekit.utils.Constants;
51  import org.orekit.utils.IERSConventions;
52  import org.orekit.utils.TimeStampedFieldPVCoordinates;
53  import org.orekit.utils.TimeStampedPVCoordinates;
54  
55  public class KlobucharModelTest {
56  
57      private static double epsilon = 1e-6;
58  
59      /** ionospheric model. */
60      private KlobucharIonoModel model;
61  
62      private UTCScale utc;
63  
64      @BeforeEach
65      public void setUp() throws Exception {
66          // Navigation message data
67          // .3820D-07   .1490D-07  -.1790D-06   .0000D-00          ION ALPHA
68          // .1430D+06   .0000D+00  -.3280D+06   .1130D+06          ION BETA
69          model = new KlobucharIonoModel(new double[]{.3820e-07, .1490e-07, -.1790e-06, 0},
70                                         new double[]{.1430e+06, 0, -.3280e+06, .1130e+06});
71  
72          Utils.setDataRoot("regular-data");
73          utc = TimeScalesFactory.getUTC();
74      }
75  
76      @AfterEach
77      public void tearDown() {
78          utc = null;
79      }
80  
81      @Test
82      public void testDelay() {
83          final double latitude = FastMath.toRadians(45);
84          final double longitude = FastMath.toRadians(2);
85          final double altitude = 500;
86          final double elevation = 70.;
87          final double azimuth = 10.;
88  
89          final AbsoluteDate date = new AbsoluteDate();
90  
91          final GeodeticPoint geo = new GeodeticPoint(latitude, longitude, altitude);
92  
93          double delayMeters = model.pathDelay(date, geo,
94                                               FastMath.toRadians(elevation),
95                                               FastMath.toRadians(azimuth),
96                                               1575.42e6, model.getParameters(date));
97  
98          Assertions.assertTrue(Precision.compareTo(delayMeters, 12., epsilon) < 0);
99          Assertions.assertTrue(Precision.compareTo(delayMeters, 0., epsilon) > 0);
100     }
101 
102     @Test
103     public void testFieldDelay() {
104         doTestFieldDelay(Binary64Field.getInstance());
105     }
106 
107     private <T extends CalculusFieldElement<T>> void doTestFieldDelay(final Field<T> field) {
108         final T zero = field.getZero();
109 
110         final T latitude  = zero.add(FastMath.toRadians(45));
111         final T longitude = zero.add(FastMath.toRadians(2));
112         final T altitude  = zero.add(500);
113         final T elevation = zero.add(FastMath.toRadians(70.));
114         final T azimuth   = zero.add(FastMath.toRadians(10.));
115 
116         final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
117 
118         final FieldGeodeticPoint<T> geo = new FieldGeodeticPoint<>(latitude, longitude, altitude);
119 
120         T delayMeters = model.pathDelay(date, geo,
121                                         elevation, azimuth,
122                                         1575.42e6, model.getParameters(field, date));
123 
124         Assertions.assertTrue(Precision.compareTo(delayMeters.getReal(), 12., epsilon) < 0);
125         Assertions.assertTrue(Precision.compareTo(delayMeters.getReal(), 0., epsilon) > 0);
126     }
127 
128     @Test
129     public void compareExpectedValue() throws IllegalArgumentException, OrekitException {
130         final double latitude = FastMath.toRadians(40);
131         final double longitude = FastMath.toRadians(-100);
132         final double altitude = 0.;
133         final double elevation = 20.;
134         final double azimuth = 210.;
135 
136         final AbsoluteDate date = new AbsoluteDate(new DateTimeComponents(2000, 1, 1,
137                                                                           20, 45, 0),
138                                                                           utc);
139 
140         final GeodeticPoint geo = new GeodeticPoint(latitude, longitude, altitude);
141 
142         final double delayMeters = model.pathDelay(date, geo,
143                                                    FastMath.toRadians(elevation),
144                                                    FastMath.toRadians(azimuth),
145                                                    1575.42e6, model.getParameters(date));
146 
147         Assertions.assertEquals(23.784, delayMeters, 0.001);
148     }
149 
150     @Test
151     public <T extends CalculusFieldElement<T>> void compareFieldExpectedValue() {
152         doCompareFieldExpectedValue(Binary64Field.getInstance());
153     }
154 
155     private <T extends CalculusFieldElement<T>> void doCompareFieldExpectedValue(final Field<T> field)
156         throws IllegalArgumentException, OrekitException {
157         final T zero = field.getZero();
158 
159         final T latitude  = zero.add(FastMath.toRadians(40));
160         final T longitude = zero.add(FastMath.toRadians(-100));
161         final T altitude  = zero;
162         final T elevation = zero.add(FastMath.toRadians(20.));
163         final T azimuth   = zero.add(FastMath.toRadians(210.));
164 
165         final AbsoluteDate date = new AbsoluteDate(new DateTimeComponents(2000, 1, 1,
166                                                                           20, 45, 0),
167                                                                           utc);
168         final FieldAbsoluteDate<T> fieldDate = new FieldAbsoluteDate<>(field, date);
169 
170         final FieldGeodeticPoint<T> geo = new FieldGeodeticPoint<>(latitude, longitude, altitude);
171 
172         final T delayMeters = model.pathDelay(fieldDate, geo,
173                                               elevation, azimuth,
174                                               1575.42e6, model.getParameters(field, fieldDate));
175 
176         Assertions.assertEquals(23.784, delayMeters.getReal(), 0.001);
177     }
178 
179     @Test
180     public void testEquality() {
181         // Common parameters
182         final AbsoluteDate date = new AbsoluteDate(2000, 1, 1, 12, 35, 04.245, TimeScalesFactory.getUTC());
183 
184         GeodeticPoint point = new GeodeticPoint(0.389, -2.962, 0);
185 
186         // Delay using azimuth/elevation angles
187         // Reference values are from EcksteinHechlerPropagatorTest class
188         final double azimuth   = 1.70;
189         final double elevation = 0.09;
190 
191         double delayAzEl = model.pathDelay(date, point, elevation, azimuth, 1575.42e6, model.getParameters(date));
192 
193         // Delay using SpacecraftState
194         final Frame ecef = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
195         final Frame eci  = FramesFactory.getEME2000();
196 
197         // Reference values are from EcksteinHechlerPropagatorTest class
198         final TimeStampedPVCoordinates pvaCoordinates = new TimeStampedPVCoordinates(date,
199                                                                                      new Vector3D(-6809069.7032156205, 3568730.540549229, 2042703.5290696376),
200                                                                                      new Vector3D(-3720.8926276445, -5588.6043568734, -2008.3138559626),
201                                                                                      new Vector3D(5.3966003438, -2.8284419922, -1.6223614396));
202         final Orbit orbit = new CartesianOrbit(pvaCoordinates, eci, 3.9860047E14);
203         final SpacecraftState state = new SpacecraftState(orbit);
204         BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
205                                                Constants.WGS84_EARTH_FLATTENING, ecef);
206         TopocentricFrame topo = new TopocentricFrame(earth, point, "Gstation");
207 
208         double delayState = model.pathDelay(state, topo, 1575.42e6, model.getParameters(date));
209 
210         // Verify
211         Assertions.assertEquals(delayAzEl, delayState, 1.0e-6);
212     }
213 
214     @Test
215     public <T extends CalculusFieldElement<T>> void testFieldEquality() {
216         doTestFieldEquality(Binary64Field.getInstance());
217     }
218 
219     private <T extends CalculusFieldElement<T>> void doTestFieldEquality(final Field<T> field) {
220         // Common parameters
221         final T zero = field.getZero();
222         final AbsoluteDate date = new AbsoluteDate(2000, 1, 1, 12, 35, 04.245, TimeScalesFactory.getUTC());
223 
224         final FieldAbsoluteDate<T> dateF = new FieldAbsoluteDate<>(field, date);
225 
226         FieldGeodeticPoint<T> point = new FieldGeodeticPoint<>(zero.add(0.389),
227                                                                zero.add(-2.962),
228                                                                zero);
229 
230         // Delay using azimuth/elevation angles
231         // Reference values are from EcksteinHechlerPropagatorTest class
232         final T azimuth   = zero.add(1.70);
233         final T elevation = zero.add(0.09);
234 
235         T delayAzEl = model.pathDelay(dateF, point, elevation, azimuth, 1575.42e6, model.getParameters(field, dateF));
236 
237         // Delay using SpacecraftState
238         final Frame ecef = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
239         final Frame eci  = FramesFactory.getEME2000();
240 
241         // Reference values are from EcksteinHechlerPropagatorTest class
242         final TimeStampedFieldPVCoordinates<T> pvaCoordinates = new TimeStampedFieldPVCoordinates<>(dateF,
243                                                                                      new FieldVector3D<>(zero.add(-6809069.7032156205), zero.add(3568730.540549229), zero.add(2042703.5290696376)),
244                                                                                      new FieldVector3D<>(zero.add(-3720.8926276445), zero.add(-5588.6043568734), zero.add(-2008.3138559626)),
245                                                                                      new FieldVector3D<>(zero.add(5.3966003438), zero.add(-2.8284419922), zero.add(-1.6223614396)));
246         final FieldOrbit<T> orbit = new FieldCartesianOrbit<>(pvaCoordinates, eci, zero.add(3.9860047E14));
247         final FieldSpacecraftState<T> state = new FieldSpacecraftState<>(orbit);
248         BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
249                                                Constants.WGS84_EARTH_FLATTENING, ecef);
250         TopocentricFrame topo = new TopocentricFrame(earth, new GeodeticPoint(0.389, -2.962, 0), "Gstation");
251 
252         T delayState = model.pathDelay(state, topo, 1575.42e6, model.getParameters(field, dateF));
253 
254         // Verify
255         Assertions.assertEquals(delayAzEl.getReal(), delayState.getReal(), 1.0e-6);
256     }
257 
258 }
259 
260