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