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