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