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