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.geometry.euclidean.threed.FieldVector3D;
24 import org.hipparchus.geometry.euclidean.threed.Vector3D;
25 import org.hipparchus.util.Binary64Field;
26 import org.hipparchus.util.FastMath;
27 import org.hipparchus.util.Precision;
28 import org.junit.jupiter.api.Assertions;
29 import org.junit.jupiter.api.BeforeEach;
30 import org.junit.jupiter.api.Test;
31 import org.orekit.Utils;
32 import org.orekit.attitudes.Attitude;
33 import org.orekit.bodies.GeodeticPoint;
34 import org.orekit.bodies.OneAxisEllipsoid;
35 import org.orekit.estimation.measurements.GroundStation;
36 import org.orekit.frames.Frame;
37 import org.orekit.frames.FramesFactory;
38 import org.orekit.frames.TopocentricFrame;
39 import org.orekit.gnss.PredefinedGnssSignal;
40 import org.orekit.orbits.FieldKeplerianOrbit;
41 import org.orekit.orbits.FieldOrbit;
42 import org.orekit.orbits.KeplerianOrbit;
43 import org.orekit.orbits.Orbit;
44 import org.orekit.orbits.OrbitType;
45 import org.orekit.orbits.PositionAngleType;
46 import org.orekit.propagation.FieldSpacecraftState;
47 import org.orekit.propagation.SpacecraftState;
48 import org.orekit.propagation.ToleranceProvider;
49 import org.orekit.time.AbsoluteDate;
50 import org.orekit.time.FieldAbsoluteDate;
51 import org.orekit.time.TimeScalesFactory;
52 import org.orekit.utils.Constants;
53 import org.orekit.utils.IERSConventions;
54 import org.orekit.utils.ParameterDriver;
55 import org.orekit.utils.ParameterDriversList;
56
57 import java.util.List;
58
59 public class EstimatedIonosphericModelTest {
60
61 @BeforeEach
62 public void setUp() throws Exception {
63 Utils.setDataRoot("regular-data");
64 }
65
66 @Test
67 public void testL1GPS() {
68
69 final IonosphericMappingFunction mapping = new SingleLayerModelMappingFunction(0.0);
70 final EstimatedIonosphericModel model = new EstimatedIonosphericModel(mapping, 1.0);
71
72 final double delay = model.pathDelay(0.5 * FastMath.PI,
73 PredefinedGnssSignal.G01.getFrequency(),
74 model.getParameters(new AbsoluteDate()));
75
76 Assertions.assertEquals(0.162, delay, 0.001);
77 }
78
79 @Test
80 public void testFieldL1GPS() {
81 doTestFieldL1GPS(Binary64Field.getInstance());
82 }
83
84 private <T extends CalculusFieldElement<T>> void doTestFieldL1GPS(final Field<T> field) {
85
86 final T zero = field.getZero();
87
88 final IonosphericMappingFunction mapping = new SingleLayerModelMappingFunction(0.0);
89 final EstimatedIonosphericModel model = new EstimatedIonosphericModel(mapping, 1.0);
90
91 final T delay = model.pathDelay(zero.add(0.5 * FastMath.PI),
92 PredefinedGnssSignal.G01.getFrequency(),
93 model.getParameters(field));
94
95 Assertions.assertEquals(0.162, delay.getReal(), 0.001);
96 }
97
98 @Test
99 public void testDelay() {
100 final double elevation = 70.;
101
102 final IonosphericMappingFunction mapping = new SingleLayerModelMappingFunction();
103 final EstimatedIonosphericModel model = new EstimatedIonosphericModel(mapping, 50.0);
104
105
106
107 double delayMeters = model.pathDelay(FastMath.toRadians(elevation),
108 PredefinedGnssSignal.G01.getFrequency(),
109 model.getParameters());
110
111 Assertions.assertTrue(Precision.compareTo(delayMeters, 12., 1.0e-6) < 0);
112 Assertions.assertTrue(Precision.compareTo(delayMeters, 0., 1.0e-6) > 0);
113 }
114
115 @Test
116 public void testFieldDelay() {
117 doTestFieldDelay(Binary64Field.getInstance());
118 }
119
120 private <T extends CalculusFieldElement<T>> void doTestFieldDelay(final Field<T> field) {
121 final double elevation = 70.;
122
123 final IonosphericMappingFunction mapping = new SingleLayerModelMappingFunction();
124 final EstimatedIonosphericModel model = new EstimatedIonosphericModel(mapping, 50.0);
125 T zero = field.getZero();
126 T delayMeters = model.pathDelay(zero.add(FastMath.toRadians(elevation)),
127 PredefinedGnssSignal.G01.getFrequency(),
128 model.getParameters(field));
129
130 Assertions.assertTrue(Precision.compareTo(delayMeters.getReal(), 12., 1.0e-6) < 0);
131 Assertions.assertTrue(Precision.compareTo(delayMeters.getReal(), 0., 1.0e-6) > 0);
132 }
133
134 @Test
135 public void testZeroDelay() {
136
137 final double frequency = PredefinedGnssSignal.G01.getFrequency();
138
139
140 final double height = 0.0;
141 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.0), FastMath.toRadians(45.0), height);
142
143 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
144 Constants.WGS84_EARTH_FLATTENING,
145 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
146
147 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "topo");
148
149
150 final IonosphericMappingFunction mapping = new SingleLayerModelMappingFunction();
151 final EstimatedIonosphericModel model = new EstimatedIonosphericModel(mapping, 10.0);
152
153
154 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
155 final Frame frame = FramesFactory.getEME2000();
156 final Orbit orbit = new KeplerianOrbit(24464560.0, 0.05, 0.122138,
157 3.10686, 1.00681, 0.048363,
158 PositionAngleType.MEAN, frame, date, Constants.WGS84_EARTH_MU);
159 final SpacecraftState state = new SpacecraftState(orbit);
160
161
162 final double delay = model.pathDelay(state, baseFrame, frequency, model.getParameters());
163 Assertions.assertEquals(0.0, delay, Double.MIN_VALUE);
164 }
165
166 @Test
167 public void testFieldZeroDelay() {
168 doTestFieldZeroDelay(Binary64Field.getInstance());
169 }
170
171 private <T extends CalculusFieldElement<T>> void doTestFieldZeroDelay(final Field<T> field) {
172 final T zero = field.getZero();
173
174 final double frequency = PredefinedGnssSignal.G01.getFrequency();
175
176
177 final double height = 0.0;
178 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.0), FastMath.toRadians(45.0), height);
179
180 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
181 Constants.WGS84_EARTH_FLATTENING,
182 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
183
184 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "topo");
185
186
187 final IonosphericMappingFunction mapping = new SingleLayerModelMappingFunction();
188 final EstimatedIonosphericModel model = new EstimatedIonosphericModel(mapping, 10.0);
189
190
191 final FieldAbsoluteDate<T> date = FieldAbsoluteDate.getJ2000Epoch(field);
192 final Frame frame = FramesFactory.getEME2000();
193 final FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(zero.add(24464560.0), zero.add(0.05), zero.add(0.122138),
194 zero.add(3.10686), zero.add(1.00681), zero.add(0.048363),
195 PositionAngleType.MEAN, frame, date, zero.add(Constants.WGS84_EARTH_MU));
196 final FieldSpacecraftState<T> state = new FieldSpacecraftState<>(orbit);
197
198
199 final T delay = model.pathDelay(state, baseFrame, frequency, model.getParameters(field));
200 Assertions.assertEquals(0.0, delay.getReal(), Double.MIN_VALUE);
201 }
202
203 @Test
204 public void testEquality() {
205 doTestEquality(Binary64Field.getInstance());
206 }
207
208 private <T extends CalculusFieldElement<T>> void doTestEquality(final Field<T> field) {
209 final double elevation = 70.;
210
211 final IonosphericMappingFunction mapping = new SingleLayerModelMappingFunction();
212 final EstimatedIonosphericModel model = new EstimatedIonosphericModel(mapping, 50.0);
213 T zero = field.getZero();
214 T delayMetersF = model.pathDelay(zero.add(FastMath.toRadians(elevation)),
215 PredefinedGnssSignal.G01.getFrequency(),
216 model.getParameters(field));
217
218 double delayMetersR = model.pathDelay(FastMath.toRadians(elevation),
219 PredefinedGnssSignal.G01.getFrequency(),
220 model.getParameters());
221
222 Assertions.assertEquals(delayMetersR, delayMetersF.getReal(), 1.0e-15);
223 }
224
225 @Test
226 public void testDelayStateDerivatives() {
227
228
229 final double frequency = PredefinedGnssSignal.G01.getFrequency();
230
231
232 final double height = 0.0;
233 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(0.0), FastMath.toRadians(0.0), height);
234
235 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
236 Constants.WGS84_EARTH_FLATTENING,
237 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
238
239 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "topo");
240
241
242 final GroundStation station = new GroundStation(baseFrame);
243
244
245 final IonosphericMappingFunction mapping = new SingleLayerModelMappingFunction();
246 final EstimatedIonosphericModel model = new EstimatedIonosphericModel(mapping, 10.0);
247
248
249 final DSFactory factory = new DSFactory(6, 1);
250 final DerivativeStructure a0 = factory.variable(0, 24464560.0);
251 final DerivativeStructure e0 = factory.variable(1, 0.05);
252 final DerivativeStructure i0 = factory.variable(2, 0.122138);
253 final DerivativeStructure pa0 = factory.variable(3, 3.10686);
254 final DerivativeStructure raan0 = factory.variable(4, 1.00681);
255 final DerivativeStructure anomaly0 = factory.variable(5, 0.048363);
256 final Field<DerivativeStructure> field = a0.getField();
257 final DerivativeStructure zero = field.getZero();
258
259
260 final FieldAbsoluteDate<DerivativeStructure> dsDate = new FieldAbsoluteDate<>(field);
261
262 final Frame frame = FramesFactory.getEME2000();
263 final FieldOrbit<DerivativeStructure> dsOrbit = new FieldKeplerianOrbit<>(a0, e0, i0, pa0, raan0, anomaly0,
264 PositionAngleType.MEAN, frame,
265 dsDate, zero.add(3.9860047e14));
266
267 final FieldSpacecraftState<DerivativeStructure> dsState = new FieldSpacecraftState<>(dsOrbit);
268
269
270 final FieldVector3D<DerivativeStructure> position = dsState.getPosition();
271 final DerivativeStructure dsElevation = baseFrame.getTrackingCoordinates(position, frame, dsDate).getElevation();
272
273
274 for (final ParameterDriver driver : model.getParametersDrivers()) {
275 driver.setReferenceDate(dsDate.toAbsoluteDate());
276 }
277
278
279 final double delayR = model.pathDelay(dsState.toSpacecraftState(), baseFrame, frequency, model.getParameters());
280 final DerivativeStructure delayD = model.pathDelay(dsState, baseFrame, frequency, model.getParameters(field));
281 Assertions.assertEquals(delayR, delayD.getValue(), 5e-15);
282
283
284 final DerivativeStructure delay = model.pathDelay(dsElevation, frequency, model.getParameters(field));
285
286 final double[] compDeriv = delay.getAllDerivatives();
287
288
289 final Orbit orbit = dsOrbit.toOrbit();
290 final SpacecraftState state = dsState.toSpacecraftState();
291
292
293 final double[][] refDeriv = new double[1][6];
294 final OrbitType orbitType = OrbitType.KEPLERIAN;
295 final PositionAngleType angleType = PositionAngleType.MEAN;
296 double dP = 0.001;
297 double[] steps = ToleranceProvider.getDefaultToleranceProvider(1000000 * dP).getTolerances(orbit, orbitType)[0];
298 for (int i = 0; i < 6; i++) {
299 SpacecraftState stateM4 = shiftState(state, orbitType, angleType, -4 * steps[i], i);
300 final Vector3D positionM4 = stateM4.getPosition();
301 final double elevationM4 = station.getBaseFrame().
302 getTrackingCoordinates(positionM4, stateM4.getFrame(), stateM4.getDate()).
303 getElevation();
304 double delayM4 = model.pathDelay(elevationM4, frequency, model.getParameters());
305
306 SpacecraftState stateM3 = shiftState(state, orbitType, angleType, -3 * steps[i], i);
307 final Vector3D positionM3 = stateM3.getPosition();
308 final double elevationM3 = station.getBaseFrame().
309 getTrackingCoordinates(positionM3, stateM3.getFrame(), stateM3.getDate()).
310 getElevation();
311 double delayM3 = model.pathDelay(elevationM3, frequency, model.getParameters());
312
313 SpacecraftState stateM2 = shiftState(state, orbitType, angleType, -2 * steps[i], i);
314 final Vector3D positionM2 = stateM2.getPosition();
315 final double elevationM2 = station.getBaseFrame().
316 getTrackingCoordinates(positionM2, stateM2.getFrame(), stateM2.getDate()).
317 getElevation();
318 double delayM2 = model.pathDelay(elevationM2, frequency, model.getParameters());
319
320 SpacecraftState stateM1 = shiftState(state, orbitType, angleType, -1 * steps[i], i);
321 final Vector3D positionM1 = stateM1.getPosition();
322 final double elevationM1 = station.getBaseFrame().
323 getTrackingCoordinates(positionM1, stateM1.getFrame(), stateM1.getDate()).
324 getElevation();
325 double delayM1 = model.pathDelay(elevationM1, frequency, model.getParameters());
326
327 SpacecraftState stateP1 = shiftState(state, orbitType, angleType, 1 * steps[i], i);
328 final Vector3D positionP1 = stateP1.getPosition();
329 final double elevationP1 = station.getBaseFrame().
330 getTrackingCoordinates(positionP1, stateP1.getFrame(), stateP1.getDate()).
331 getElevation();
332 double delayP1 = model.pathDelay(elevationP1, frequency, model.getParameters());
333
334 SpacecraftState stateP2 = shiftState(state, orbitType, angleType, 2 * steps[i], i);
335 final Vector3D positionP2 = stateP2.getPosition();
336 final double elevationP2 = station.getBaseFrame().
337 getTrackingCoordinates(positionP2, stateP2.getFrame(), stateP2.getDate()).
338 getElevation();
339 double delayP2 = model.pathDelay(elevationP2, frequency, model.getParameters());
340
341 SpacecraftState stateP3 = shiftState(state, orbitType, angleType, 3 * steps[i], i);
342 final Vector3D positionP3 = stateP3.getPosition();
343 final double elevationP3 = station.getBaseFrame().
344 getTrackingCoordinates(positionP3, stateP3.getFrame(), stateP3.getDate()).
345 getElevation();
346 double delayP3 = model.pathDelay(elevationP3, frequency, model.getParameters());
347
348 SpacecraftState stateP4 = shiftState(state, orbitType, angleType, 4 * steps[i], i);
349 final Vector3D positionP4 = stateP4.getPosition();
350 final double elevationP4 = station.getBaseFrame().
351 getTrackingCoordinates(positionP4, stateP4.getFrame(), stateP4.getDate()).
352 getElevation();
353 double delayP4 = model.pathDelay(elevationP4, frequency, model.getParameters());
354
355 fillJacobianColumn(refDeriv, i, steps[i],
356 delayM4, delayM3, delayM2, delayM1,
357 delayP1, delayP2, delayP3, delayP4);
358 }
359
360 for (int i = 0; i < 6; i++) {
361 Assertions.assertEquals(compDeriv[i + 1], refDeriv[0][i], 8.3e-11);
362 }
363 }
364
365 private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType, PositionAngleType angleType,
366 double delta, int column) {
367
368 double[][] array = stateToArray(state, orbitType, angleType, true);
369 array[0][column] += delta;
370
371 return arrayToState(array, orbitType, angleType, state.getFrame(), state.getDate(),
372 state.getOrbit().getMu(), state.getAttitude());
373
374 }
375
376 @Test
377 public void testParametersDerivatives() {
378
379
380 final double frequency = PredefinedGnssSignal.G01.getFrequency();
381
382
383 final double latitude = FastMath.toRadians(45.0);
384 final double longitude = FastMath.toRadians(45.0);
385 final double height = 0.0;
386 final GeodeticPoint point = new GeodeticPoint(latitude, longitude, height);
387
388 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
389 Constants.WGS84_EARTH_FLATTENING,
390 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
391
392 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "topo");
393
394
395 final IonosphericMappingFunction mapping = new SingleLayerModelMappingFunction();
396 final EstimatedIonosphericModel model = new EstimatedIonosphericModel(mapping, 50.0);
397
398
399 for (final ParameterDriver driver : model.getParametersDrivers()) {
400 driver.setValue(driver.getReferenceValue());
401 driver.setSelected(driver.getName().equals(EstimatedIonosphericModel.VERTICAL_TOTAL_ELECTRON_CONTENT));
402 }
403
404
405 int nbParams = 0;
406 for (final ParameterDriver driver : model.getParametersDrivers()) {
407 if (driver.isSelected()) {
408 ++nbParams;
409 }
410 }
411
412
413 final DSFactory factory = new DSFactory(6 + nbParams, 1);
414 final DerivativeStructure a0 = factory.variable(0, 24464560.0);
415 final DerivativeStructure e0 = factory.variable(1, 0.05);
416 final DerivativeStructure i0 = factory.variable(2, 0.122138);
417 final DerivativeStructure pa0 = factory.variable(3, 3.10686);
418 final DerivativeStructure raan0 = factory.variable(4, 1.00681);
419 final DerivativeStructure anomaly0 = factory.variable(5, 0.048363);
420 final Field<DerivativeStructure> field = a0.getField();
421 final DerivativeStructure zero = field.getZero();
422
423
424 final FieldAbsoluteDate<DerivativeStructure> dsDate = new FieldAbsoluteDate<>(field, 2018, 11, 19, 18, 0, 0.0,
425 TimeScalesFactory.getUTC());
426
427
428 for (final ParameterDriver driver : model.getParametersDrivers()) {
429 driver.setReferenceDate(dsDate.toAbsoluteDate());
430 }
431
432
433 final Frame frame = FramesFactory.getEME2000();
434 final FieldOrbit<DerivativeStructure> dsOrbit = new FieldKeplerianOrbit<>(a0, e0, i0, pa0, raan0, anomaly0,
435 PositionAngleType.MEAN, frame,
436 dsDate, zero.add(3.9860047e14));
437
438
439 final FieldSpacecraftState<DerivativeStructure> dsState = new FieldSpacecraftState<>(dsOrbit);
440
441
442 final FieldVector3D<DerivativeStructure> position = dsState.getPosition();
443 final DerivativeStructure dsElevation = baseFrame.getTrackingCoordinates(position, frame, dsState.getDate()).getElevation();
444
445
446 final List<ParameterDriver> drivers = model.getParametersDrivers();
447 final DerivativeStructure[] parameters = new DerivativeStructure[drivers.size()];
448 int index = 6;
449 for (int i = 0; i < drivers.size(); ++i) {
450 parameters[i] = drivers.get(i).isSelected() ?
451 factory.variable(index++, drivers.get(i).getValue()) :
452 factory.constant(drivers.get(i).getValue());
453 }
454
455
456 final DerivativeStructure delay = model.pathDelay(dsElevation, frequency, parameters);
457
458 final double[] compDeriv = delay.getAllDerivatives();
459
460
461 final double elevation = dsElevation.getReal();
462
463
464 final double[][] refDeriv = new double[1][1];
465 ParameterDriversList bound = new ParameterDriversList();
466 for (final ParameterDriver driver : model.getParametersDrivers()) {
467 if (driver.getName().equals(EstimatedIonosphericModel.VERTICAL_TOTAL_ELECTRON_CONTENT)) {
468 driver.setSelected(true);
469 bound.add(driver);
470 } else {
471 driver.setSelected(false);
472 }
473 }
474 ParameterDriver selected = bound.getDrivers().get(0);
475 double p0 = selected.getReferenceValue();
476 double h = selected.getScale();
477
478 selected.setValue(p0 - 4 * h);
479 double delayM4 = model.pathDelay(elevation, frequency, model.getParameters());
480
481 selected.setValue(p0 - 3 * h);
482 double delayM3 = model.pathDelay(elevation, frequency, model.getParameters());
483
484 selected.setValue(p0 - 2 * h);
485 double delayM2 = model.pathDelay(elevation, frequency, model.getParameters());
486
487 selected.setValue(p0 - 1 * h);
488 double delayM1 = model.pathDelay(elevation, frequency, model.getParameters());
489
490 selected.setValue(p0 + 1 * h);
491 double delayP1 = model.pathDelay(elevation, frequency, model.getParameters());
492
493 selected.setValue(p0 + 2 * h);
494 double delayP2 = model.pathDelay(elevation, frequency, model.getParameters());
495
496 selected.setValue(p0 + 3 * h);
497 double delayP3 = model.pathDelay(elevation, frequency, model.getParameters());
498
499 selected.setValue(p0 + 4 * h);
500 double delayP4 = model.pathDelay(elevation, frequency, model.getParameters());
501
502 fillJacobianColumn(refDeriv, 0, h,
503 delayM4, delayM3, delayM2, delayM1,
504 delayP1, delayP2, delayP3, delayP4);
505
506 Assertions.assertEquals(compDeriv[7], refDeriv[0][0], 1.0e-15);
507
508 }
509
510 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType, PositionAngleType angleType,
511 boolean withMass) {
512 double[][] array = new double[2][withMass ? 7 : 6];
513 orbitType.mapOrbitToArray(state.getOrbit(), angleType, array[0], array[1]);
514 if (withMass) {
515 array[0][6] = state.getMass();
516 }
517 return array;
518 }
519
520 private SpacecraftState arrayToState(double[][] array, OrbitType orbitType, PositionAngleType angleType,
521 Frame frame, AbsoluteDate date, double mu,
522 Attitude attitude) {
523 Orbit orbit = orbitType.mapArrayToOrbit(array[0], array[1], angleType, date, mu, frame);
524 return (array.length > 6) ?
525 new SpacecraftState(orbit, attitude) :
526 new SpacecraftState(orbit, attitude).withMass(array[0][6]);
527 }
528
529 private void fillJacobianColumn(double[][] jacobian, int column, double h,
530 double sM4h, double sM3h,
531 double sM2h, double sM1h,
532 double sP1h, double sP2h,
533 double sP3h, double sP4h) {
534
535 jacobian[0][column] = ( -3 * (sP4h - sM4h) +
536 32 * (sP3h - sM3h) -
537 168 * (sP2h - sM2h) +
538 672 * (sP1h - sM1h)) / (840 * h);
539 }
540
541 }