1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.models.earth.ionosphere;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.Field;
21 import org.hipparchus.analysis.differentiation.DSFactory;
22 import org.hipparchus.analysis.differentiation.DerivativeStructure;
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.attitudes.Attitude;
30 import org.orekit.bodies.GeodeticPoint;
31 import org.orekit.bodies.OneAxisEllipsoid;
32 import org.orekit.frames.Frame;
33 import org.orekit.frames.FramesFactory;
34 import org.orekit.frames.TopocentricFrame;
35 import org.orekit.gnss.PredefinedGnssSignal;
36 import org.orekit.gnss.metric.messages.ssr.subtype.SsrIm201;
37 import org.orekit.gnss.metric.messages.ssr.subtype.SsrIm201Data;
38 import org.orekit.gnss.metric.messages.ssr.subtype.SsrIm201Header;
39 import org.orekit.orbits.FieldKeplerianOrbit;
40 import org.orekit.orbits.FieldOrbit;
41 import org.orekit.orbits.KeplerianOrbit;
42 import org.orekit.orbits.Orbit;
43 import org.orekit.orbits.OrbitType;
44 import org.orekit.orbits.PositionAngleType;
45 import org.orekit.propagation.FieldSpacecraftState;
46 import org.orekit.propagation.SpacecraftState;
47 import org.orekit.propagation.ToleranceProvider;
48 import org.orekit.time.AbsoluteDate;
49 import org.orekit.time.FieldAbsoluteDate;
50 import org.orekit.utils.Constants;
51 import org.orekit.utils.IERSConventions;
52 import org.orekit.utils.ParameterDriver;
53
54 import java.util.Collections;
55
56 public class SsrVtecIonosphericModelTest {
57
58 private SsrIm201 vtecMessage;
59
60 @BeforeEach
61 public void setUp() {
62 Utils.setDataRoot("regular-data");
63
64
65 final SsrIm201Header header = new SsrIm201Header();
66 header.setSsrEpoch1s(302400.0);
67 header.setNumberOfIonosphericLayers(1);
68
69
70 final SsrIm201Data data = new SsrIm201Data();
71 data.setHeightIonosphericLayer(650000.0);
72 data.setSphericalHarmonicsDegree(2);
73 data.setSphericalHarmonicsOrder(1);
74 data.setCnm(new double[][] { {18.2, 0.0},
75 {13.4, 26.2},
76 {14.7, 6.8} });
77 data.setSnm(new double[][] { {0.0, 0.0},
78 {0.0, 6.2},
79 {0.0, 17.6} });
80
81
82
83 vtecMessage = new SsrIm201(201, header, Collections.singletonList(data));
84
85 }
86
87 @Test
88 public void testDelay() {
89
90
91 final double frequency = PredefinedGnssSignal.G01.getFrequency();
92
93
94 final double height = 0.0;
95 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.0), FastMath.toRadians(25.0), height);
96
97 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
98 Constants.WGS84_EARTH_FLATTENING,
99 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
100
101 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "topo");
102
103
104 final SsrVtecIonosphericModel model = new SsrVtecIonosphericModel(vtecMessage);
105
106
107 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
108 final Frame frame = FramesFactory.getEME2000();
109 final Orbit orbit = new KeplerianOrbit(24464560.0, 0.05, 0.122138,
110 3.10686, 1.00681, 0.048363,
111 PositionAngleType.MEAN, frame, date, Constants.WGS84_EARTH_MU);
112 final SpacecraftState state = new SpacecraftState(orbit);
113
114
115 final double delay = model.pathDelay(state, baseFrame, frequency, model.getParameters());
116 Assertions.assertEquals(13.488, delay, 0.001);
117
118 }
119
120 @Test
121 public void testFieldDelay() {
122 doTestFieldDelay(Binary64Field.getInstance());
123 }
124
125 private <T extends CalculusFieldElement<T>> void doTestFieldDelay(final Field<T> field) {
126 final T zero = field.getZero();
127
128 final double frequency = PredefinedGnssSignal.G01.getFrequency();
129
130
131 final double height = 0.0;
132 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.0), FastMath.toRadians(25.0), height);
133
134 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
135 Constants.WGS84_EARTH_FLATTENING,
136 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
137
138 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "topo");
139
140
141 final SsrVtecIonosphericModel model = new SsrVtecIonosphericModel(vtecMessage);
142
143
144 final FieldAbsoluteDate<T> date = FieldAbsoluteDate.getJ2000Epoch(field);
145 final Frame frame = FramesFactory.getEME2000();
146 final FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(24464560.0), zero.add(0.05), zero.add(0.122138),
147 zero.add(3.10686), zero.add(1.00681), zero.add(0.048363),
148 PositionAngleType.MEAN, frame, date, zero.add(Constants.WGS84_EARTH_MU));
149 final FieldSpacecraftState<T> state = new FieldSpacecraftState<>(orbit);
150
151
152 final T delay = model.pathDelay(state, baseFrame, frequency, model.getParameters(field));
153 Assertions.assertEquals(13.488, delay.getReal(), 0.001);
154 }
155
156 @Test
157 public void testZeroDelay() {
158
159 final double frequency = PredefinedGnssSignal.G01.getFrequency();
160
161
162 final double height = 0.0;
163 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.0), FastMath.toRadians(45.0), height);
164
165 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
166 Constants.WGS84_EARTH_FLATTENING,
167 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
168
169 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "topo");
170
171
172 final SsrVtecIonosphericModel model = new SsrVtecIonosphericModel(vtecMessage);
173
174
175 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
176 final Frame frame = FramesFactory.getEME2000();
177 final Orbit orbit = new KeplerianOrbit(24464560.0, 0.05, 0.122138,
178 3.10686, 1.00681, 0.048363,
179 PositionAngleType.MEAN, frame, date, Constants.WGS84_EARTH_MU);
180 final SpacecraftState state = new SpacecraftState(orbit);
181
182
183 final double delay = model.pathDelay(state, baseFrame, frequency, model.getParameters());
184 Assertions.assertEquals(0.0, delay, Double.MIN_VALUE);
185 }
186
187 @Test
188 public void testFieldZeroDelay() {
189 doTestFieldZeroDelay(Binary64Field.getInstance());
190 }
191
192 private <T extends CalculusFieldElement<T>> void doTestFieldZeroDelay(final Field<T> field) {
193 final T zero = field.getZero();
194
195 final double frequency = PredefinedGnssSignal.G01.getFrequency();
196
197
198 final double height = 0.0;
199 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.0), FastMath.toRadians(45.0), height);
200
201 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
202 Constants.WGS84_EARTH_FLATTENING,
203 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
204
205 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "topo");
206
207
208 final SsrVtecIonosphericModel model = new SsrVtecIonosphericModel(vtecMessage);
209
210
211 final FieldAbsoluteDate<T> date = FieldAbsoluteDate.getJ2000Epoch(field);
212 final Frame frame = FramesFactory.getEME2000();
213 final FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(24464560.0), zero.add(0.05), zero.add(0.122138),
214 zero.add(3.10686), zero.add(1.00681), zero.add(0.048363),
215 PositionAngleType.MEAN, frame, date, zero.add(Constants.WGS84_EARTH_MU));
216 final FieldSpacecraftState<T> state = new FieldSpacecraftState<>(orbit);
217
218
219 final T delay = model.pathDelay(state, baseFrame, frequency, model.getParameters(field));
220 Assertions.assertEquals(0.0, delay.getReal(), Double.MIN_VALUE);
221 }
222
223 @Test
224 public void testDelayStateDerivatives() {
225
226
227 final double frequency = PredefinedGnssSignal.G01.getFrequency();
228
229
230 final double height = 0.0;
231 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(0.0), FastMath.toRadians(0.0), height);
232
233 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
234 Constants.WGS84_EARTH_FLATTENING,
235 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
236
237 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "topo");
238
239
240 final SsrVtecIonosphericModel model = new SsrVtecIonosphericModel(vtecMessage);
241
242
243 final DSFactory factory = new DSFactory(6, 1);
244 final DerivativeStructure a0 = factory.variable(0, 24464560.0);
245 final DerivativeStructure e0 = factory.variable(1, 0.05);
246 final DerivativeStructure i0 = factory.variable(2, 0.122138);
247 final DerivativeStructure pa0 = factory.variable(3, 3.10686);
248 final DerivativeStructure raan0 = factory.variable(4, 1.00681);
249 final DerivativeStructure anomaly0 = factory.variable(5, 0.048363);
250 final Field<DerivativeStructure> field = a0.getField();
251 final DerivativeStructure zero = field.getZero();
252
253
254 final FieldAbsoluteDate<DerivativeStructure> dsDate = new FieldAbsoluteDate<>(field);
255
256 final Frame frame = FramesFactory.getEME2000();
257 final FieldOrbit<DerivativeStructure> dsOrbit = new FieldKeplerianOrbit<>(a0, e0, i0, pa0, raan0, anomaly0,
258 PositionAngleType.MEAN, frame,
259 dsDate, zero.add(3.9860047e14));
260
261 final FieldSpacecraftState<DerivativeStructure> dsState = new FieldSpacecraftState<>(dsOrbit);
262
263
264 for (final ParameterDriver driver : model.getParametersDrivers()) {
265 driver.setReferenceDate(dsDate.toAbsoluteDate());
266 }
267
268
269 final DerivativeStructure delay = model.pathDelay(dsState, baseFrame, frequency, model.getParameters(field));
270
271 final double[] compDeriv = delay.getAllDerivatives();
272
273
274 final Orbit orbit = dsOrbit.toOrbit();
275 final SpacecraftState state = dsState.toSpacecraftState();
276
277
278 final double[][] refDeriv = new double[1][6];
279 final OrbitType orbitType = OrbitType.KEPLERIAN;
280 final PositionAngleType angleType = PositionAngleType.MEAN;
281 double dP = 0.001;
282 double[] steps = ToleranceProvider.getDefaultToleranceProvider(1000000 * dP).getTolerances(orbit, orbitType)[0];
283 for (int i = 0; i < 6; i++) {
284 SpacecraftState stateM4 = shiftState(state, orbitType, angleType, -4 * steps[i], i);
285 double delayM4 = model.pathDelay(stateM4, baseFrame, frequency, model.getParameters());
286
287 SpacecraftState stateM3 = shiftState(state, orbitType, angleType, -3 * steps[i], i);
288 double delayM3 = model.pathDelay(stateM3, baseFrame, frequency, model.getParameters());
289
290 SpacecraftState stateM2 = shiftState(state, orbitType, angleType, -2 * steps[i], i);
291 double delayM2 = model.pathDelay(stateM2, baseFrame, frequency, model.getParameters());
292
293 SpacecraftState stateM1 = shiftState(state, orbitType, angleType, -1 * steps[i], i);
294 double delayM1 = model.pathDelay(stateM1, baseFrame, frequency, model.getParameters());
295
296 SpacecraftState stateP1 = shiftState(state, orbitType, angleType, 1 * steps[i], i);
297 double delayP1 = model.pathDelay(stateP1, baseFrame, frequency, model.getParameters());
298
299 SpacecraftState stateP2 = shiftState(state, orbitType, angleType, 2 * steps[i], i);
300 double delayP2 = model.pathDelay(stateP2, baseFrame, frequency, model.getParameters());
301
302 SpacecraftState stateP3 = shiftState(state, orbitType, angleType, 3 * steps[i], i);
303 double delayP3 = model.pathDelay(stateP3, baseFrame, frequency, model.getParameters());
304
305 SpacecraftState stateP4 = shiftState(state, orbitType, angleType, 4 * steps[i], i);
306 double delayP4 = model.pathDelay(stateP4, baseFrame, frequency, model.getParameters());
307
308 fillJacobianColumn(refDeriv, i, steps[i],
309 delayM4, delayM3, delayM2, delayM1,
310 delayP1, delayP2, delayP3, delayP4);
311 }
312
313 for (int i = 0; i < 6; i++) {
314 Assertions.assertEquals(compDeriv[i + 1], refDeriv[0][i], 2.3e-11);
315 }
316 }
317
318 @Test
319 public void testDelayRange() {
320
321
322 final double frequency = PredefinedGnssSignal.G01.getFrequency();
323
324
325 final double height = 0.0;
326
327
328 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
329 Constants.WGS84_EARTH_FLATTENING,
330 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
331
332
333 final SsrVtecIonosphericModel model = new SsrVtecIonosphericModel(vtecMessage);
334
335
336 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
337 final Frame frame = FramesFactory.getEME2000();
338 final Orbit orbit = new KeplerianOrbit(24464560.0, 0.05, 0.122138,
339 3.10686, 1.00681, 0.048363,
340 PositionAngleType.MEAN, frame, date, Constants.WGS84_EARTH_MU);
341 final SpacecraftState state = new SpacecraftState(orbit);
342
343
344 for (double latitude = -90.0; latitude <= 90.0; latitude = latitude + 5.0) {
345 for (double longitude = -180.0; longitude <= 180.0; longitude += 10.0) {
346 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(latitude), FastMath.toRadians(longitude), height);
347 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "topo");
348 final double delay = model.pathDelay(state, baseFrame, frequency, model.getParameters());
349 Assertions.assertTrue(delay >= 0 && delay < 20.0);
350 }
351 }
352
353 }
354
355 @Test
356 public void testFieldDelayRange() {
357 doTestFieldDelayRange(Binary64Field.getInstance());
358 }
359
360 private <T extends CalculusFieldElement<T>> void doTestFieldDelayRange(final Field<T> field) {
361 final T zero = field.getZero();
362
363 final double frequency = PredefinedGnssSignal.G01.getFrequency();
364
365
366 final double height = 0.0;
367
368
369 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
370 Constants.WGS84_EARTH_FLATTENING,
371 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
372
373
374 final SsrVtecIonosphericModel model = new SsrVtecIonosphericModel(vtecMessage);
375
376
377 final FieldAbsoluteDate<T> date = FieldAbsoluteDate.getJ2000Epoch(field);
378 final Frame frame = FramesFactory.getEME2000();
379 final FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(24464560.0), zero.add(0.05), zero.add(0.122138),
380 zero.add(3.10686), zero.add(1.00681), zero.add(0.048363),
381 PositionAngleType.MEAN, frame, date, zero.add(Constants.WGS84_EARTH_MU));
382 final FieldSpacecraftState<T> state = new FieldSpacecraftState<>(orbit);
383
384
385 for (double latitude = -90.0; latitude <= 90.0; latitude = latitude + 5.0) {
386 for (double longitude = -180.0; longitude <= 180.0; longitude += 10.0) {
387 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(latitude), FastMath.toRadians(longitude), height);
388 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "topo");
389 final T delay = model.pathDelay(state, baseFrame, frequency, model.getParameters(field));
390 Assertions.assertTrue(delay.getReal() >= 0 && delay.getReal() < 20.0);
391 }
392 }
393
394 }
395
396 private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType, PositionAngleType angleType,
397 double delta, int column) {
398
399 double[][] array = stateToArray(state, orbitType, angleType, true);
400 array[0][column] += delta;
401
402 return arrayToState(array, orbitType, angleType, state.getFrame(), state.getDate(),
403 state.getOrbit().getMu(), state.getAttitude());
404
405 }
406
407 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType, PositionAngleType angleType,
408 boolean withMass) {
409 double[][] array = new double[2][withMass ? 7 : 6];
410 orbitType.mapOrbitToArray(state.getOrbit(), angleType, array[0], array[1]);
411 if (withMass) {
412 array[0][6] = state.getMass();
413 }
414 return array;
415 }
416
417 private SpacecraftState arrayToState(double[][] array, OrbitType orbitType, PositionAngleType angleType,
418 Frame frame, AbsoluteDate date, double mu,
419 Attitude attitude) {
420 Orbit orbit = orbitType.mapArrayToOrbit(array[0], array[1], angleType, date, mu, frame);
421 return (array.length > 6) ?
422 new SpacecraftState(orbit, attitude) :
423 new SpacecraftState(orbit, attitude, array[0][6]);
424 }
425
426 private void fillJacobianColumn(double[][] jacobian, int column, double h,
427 double sM4h, double sM3h,
428 double sM2h, double sM1h,
429 double sP1h, double sP2h,
430 double sP3h, double sP4h) {
431
432 jacobian[0][column] = ( -3 * (sP4h - sM4h) +
433 32 * (sP3h - sM3h) -
434 168 * (sP2h - sM2h) +
435 672 * (sP1h - sM1h)) / (840 * h);
436 }
437
438 }