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 org.hipparchus.Field;
20 import org.hipparchus.CalculusFieldElement;
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.Decimal64Field;
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
56 public class FieldViennaThreeModelTest {
57
58 private static double epsilon = 1e-6;
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 testMappingFactors() {
72 doTestMappingFactors(Decimal64Field.getInstance());
73 }
74
75 private <T extends CalculusFieldElement<T>> void doTestMappingFactors(final Field<T> field) {
76
77 final T zero = field.getZero();
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2018, 11, 25, TimeScalesFactory.getUTC());
100
101 final double latitude = FastMath.toRadians(37.5);
102 final double longitude = FastMath.toRadians(277.5);
103 final double height = 824.0;
104
105 final double elevation = FastMath.toRadians(38.0);
106 final double expectedHydro = 1.621024;
107 final double expectedWet = 1.623023;
108
109 final double[] a = {0.00123462, 0.00047101};
110 final double[] z = {2.1993, 0.0690};
111
112 final FieldGeodeticPoint<T> point = new FieldGeodeticPoint<>(zero.add(latitude), zero.add(longitude), zero.add(height));
113
114 final ViennaThreeModel model = new ViennaThreeModel(a, z);
115
116 final T[] computedMapping = model.mappingFactors(zero.add(elevation), point,
117 date);
118
119 Assert.assertEquals(expectedHydro, computedMapping[0].getReal(), epsilon);
120 Assert.assertEquals(expectedWet, computedMapping[1].getReal(), epsilon);
121 }
122
123 @Test
124 public void testLowElevation() {
125 doTestLowElevation(Decimal64Field.getInstance());
126 }
127
128 private <T extends CalculusFieldElement<T>> void doTestLowElevation(final Field<T> field) {
129
130 final T zero = field.getZero();
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2018, 11, 25, TimeScalesFactory.getUTC());
153
154 final double latitude = FastMath.toRadians(37.5);
155 final double longitude = FastMath.toRadians(277.5);
156 final double height = 824.0;
157
158 final double elevation = FastMath.toRadians(5.0);
159 final double expectedHydro = 10.132802;
160 final double expectedWet = 10.879154;
161
162 final double[] a = {0.00123462, 0.00047101};
163 final double[] z = {2.1993, 0.0690};
164
165 final FieldGeodeticPoint<T> point = new FieldGeodeticPoint<>(zero.add(latitude), zero.add(longitude), zero.add(height));
166
167 final ViennaThreeModel model = new ViennaThreeModel(a, z);
168
169 final T[] computedMapping = model.mappingFactors(zero.add(elevation), point,
170 date);
171
172 Assert.assertEquals(expectedHydro, computedMapping[0].getReal(), epsilon);
173 Assert.assertEquals(expectedWet, computedMapping[1].getReal(), epsilon);
174 }
175
176 @Test
177 public void testHightElevation() {
178 doTestHightElevation(Decimal64Field.getInstance());
179 }
180
181 private <T extends CalculusFieldElement<T>> void doTestHightElevation(final Field<T> field) {
182
183 final T zero = field.getZero();
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2018, 11, 25, TimeScalesFactory.getUTC());
206
207 final double latitude = FastMath.toRadians(37.5);
208 final double longitude = FastMath.toRadians(277.5);
209 final double height = 824.0;
210
211 final double elevation = FastMath.toRadians(85.0);
212 final double expectedHydro = 1.003810;
213 final double expectedWet = 1.003816;
214
215 final double[] a = {0.00123462, 0.00047101};
216 final double[] z = {2.1993, 0.0690};
217
218 final FieldGeodeticPoint<T> point = new FieldGeodeticPoint<>(zero.add(latitude), zero.add(longitude), zero.add(height));
219
220 final ViennaThreeModel model = new ViennaThreeModel(a, z);
221
222 final T[] computedMapping = model.mappingFactors(zero.add(elevation), point,
223 date);
224
225 Assert.assertEquals(expectedHydro, computedMapping[0].getReal(), epsilon);
226 Assert.assertEquals(expectedWet, computedMapping[1].getReal(), epsilon);
227 }
228
229 @Test
230 public void testDelay() {
231 doTestDelay(Decimal64Field.getInstance());
232 }
233
234 private <T extends CalculusFieldElement<T>> void doTestDelay(final Field<T> field) {
235 final T zero = field.getZero();
236 final double elevation = 10d;
237 final double height = 100d;
238 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
239 final double[] a = { 0.00123462, 0.00047101};
240 final double[] z = {2.1993, 0.0690};
241 final FieldGeodeticPoint<T> point = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(37.5)), zero.add(FastMath.toRadians(277.5)), zero.add(height));
242 ViennaThreeModel model = new ViennaThreeModel(a, z);
243 final T path = model.pathDelay(zero.add(FastMath.toRadians(elevation)), point,
244 model.getParameters(field), date);
245 Assert.assertTrue(Precision.compareTo(path.getReal(), 20d, epsilon) < 0);
246 Assert.assertTrue(Precision.compareTo(path.getReal(), 0d, epsilon) > 0);
247 }
248
249 @Test
250 public void testFixedHeight() {
251 doTestFixedHeight(Decimal64Field.getInstance());
252 }
253
254 private <T extends CalculusFieldElement<T>> void doTestFixedHeight(final Field<T> field) {
255 final T zero = field.getZero();
256 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
257 final double[] a = { 0.00123462, 0.00047101};
258 final double[] z = {2.1993, 0.0690};
259 final FieldGeodeticPoint<T> point = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(37.5)), zero.add(FastMath.toRadians(277.5)), zero.add(350.0));
260 ViennaThreeModel model = new ViennaThreeModel(a, z);
261 T lastDelay = zero.add(Double.MAX_VALUE);
262
263 for (double elev = 10d; elev < 90d; elev += 8d) {
264 final T delay = model.pathDelay(zero.add(FastMath.toRadians(elev)), point,
265 model.getParameters(field), date);
266 Assert.assertTrue(Precision.compareTo(delay.getReal(), lastDelay.getReal(), epsilon) < 0);
267 lastDelay = delay;
268 }
269 }
270
271 @Test
272 public void testDelayStateDerivatives() {
273
274
275 final double latitude = FastMath.toRadians(45.0);
276 final double longitude = FastMath.toRadians(45.0);
277 final double height = 0.0;
278 final GeodeticPoint point = new GeodeticPoint(latitude, longitude, height);
279
280 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
281 Constants.WGS84_EARTH_FLATTENING,
282 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
283
284 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "topo");
285
286
287 final GroundStation station = new GroundStation(baseFrame);
288
289
290 final double[] a = { 0.00127683, 0.00060955 };
291 final double[] z = {2.0966, 0.2140};
292 final DiscreteTroposphericModel model = new ViennaThreeModel(a, z);
293
294
295 final DSFactory factory = new DSFactory(6, 1);
296 final DerivativeStructure a0 = factory.variable(0, 24464560.0);
297 final DerivativeStructure e0 = factory.variable(1, 0.05);
298 final DerivativeStructure i0 = factory.variable(2, 0.122138);
299 final DerivativeStructure pa0 = factory.variable(3, 3.10686);
300 final DerivativeStructure raan0 = factory.variable(4, 1.00681);
301 final DerivativeStructure anomaly0 = factory.variable(5, 0.048363);
302 final Field<DerivativeStructure> field = a0.getField();
303 final DerivativeStructure zero = field.getZero();
304
305
306 final FieldAbsoluteDate<DerivativeStructure> dsDate = new FieldAbsoluteDate<>(field, 2018, 11, 19, 18, 0, 0.0,
307 TimeScalesFactory.getUTC());
308
309 final Frame frame = FramesFactory.getEME2000();
310 final FieldOrbit<DerivativeStructure> dsOrbit = new FieldKeplerianOrbit<>(a0, e0, i0, pa0, raan0, anomaly0,
311 PositionAngle.MEAN, frame,
312 dsDate, zero.add(3.9860047e14));
313
314 final FieldSpacecraftState<DerivativeStructure> dsState = new FieldSpacecraftState<>(dsOrbit);
315
316
317 final FieldVector3D<DerivativeStructure> position = dsState.getPVCoordinates().getPosition();
318 final DerivativeStructure dsElevation = baseFrame.getElevation(position, frame, dsDate);
319
320
321 final FieldGeodeticPoint<DerivativeStructure> dsPoint = new FieldGeodeticPoint<>(zero.add(latitude), zero.add(longitude), zero.add(height));
322 final DerivativeStructure delay = model.pathDelay(dsElevation, dsPoint, model.getParameters(field), dsDate);
323
324 final double[] compDelay = delay.getAllDerivatives();
325
326
327 final Orbit orbit = dsOrbit.toOrbit();
328 final SpacecraftState state = dsState.toSpacecraftState();
329
330
331 final double[][] refDeriv = new double[1][6];
332 final OrbitType orbitType = OrbitType.KEPLERIAN;
333 final PositionAngle angleType = PositionAngle.MEAN;
334 double dP = 0.001;
335 double[] steps = NumericalPropagator.tolerances(1000000 * dP, orbit, orbitType)[0];
336 for (int i = 0; i < 6; i++) {
337 SpacecraftState stateM4 = shiftState(state, orbitType, angleType, -4 * steps[i], i);
338 final Vector3D positionM4 = stateM4.getPVCoordinates().getPosition();
339 final double elevationM4 = station.getBaseFrame().getElevation(positionM4, stateM4.getFrame(), stateM4.getDate());
340 double delayM4 = model.pathDelay(elevationM4, point, model.getParameters(), stateM4.getDate());
341
342 SpacecraftState stateM3 = shiftState(state, orbitType, angleType, -3 * steps[i], i);
343 final Vector3D positionM3 = stateM3.getPVCoordinates().getPosition();
344 final double elevationM3 = station.getBaseFrame().getElevation(positionM3, stateM3.getFrame(), stateM3.getDate());
345 double delayM3 = model.pathDelay(elevationM3, point, model.getParameters(), stateM3.getDate());
346
347 SpacecraftState stateM2 = shiftState(state, orbitType, angleType, -2 * steps[i], i);
348 final Vector3D positionM2 = stateM2.getPVCoordinates().getPosition();
349 final double elevationM2 = station.getBaseFrame().getElevation(positionM2, stateM2.getFrame(), stateM2.getDate());
350 double delayM2 = model.pathDelay(elevationM2, point, model.getParameters(), stateM2.getDate());
351
352 SpacecraftState stateM1 = shiftState(state, orbitType, angleType, -1 * steps[i], i);
353 final Vector3D positionM1 = stateM1.getPVCoordinates().getPosition();
354 final double elevationM1 = station.getBaseFrame().getElevation(positionM1, stateM1.getFrame(), stateM1.getDate());
355 double delayM1 = model.pathDelay(elevationM1, point, model.getParameters(), stateM1.getDate());
356
357 SpacecraftState stateP1 = shiftState(state, orbitType, angleType, 1 * steps[i], i);
358 final Vector3D positionP1 = stateP1.getPVCoordinates().getPosition();
359 final double elevationP1 = station.getBaseFrame().getElevation(positionP1, stateP1.getFrame(), stateP1.getDate());
360 double delayP1 = model.pathDelay(elevationP1, point, model.getParameters(), stateP1.getDate());
361
362 SpacecraftState stateP2 = shiftState(state, orbitType, angleType, 2 * steps[i], i);
363 final Vector3D positionP2 = stateP2.getPVCoordinates().getPosition();
364 final double elevationP2 = station.getBaseFrame().getElevation(positionP2, stateP2.getFrame(), stateP2.getDate());
365 double delayP2 = model.pathDelay(elevationP2, point, model.getParameters(), stateP2.getDate());
366
367 SpacecraftState stateP3 = shiftState(state, orbitType, angleType, 3 * steps[i], i);
368 final Vector3D positionP3 = stateP3.getPVCoordinates().getPosition();
369 final double elevationP3 = station.getBaseFrame().getElevation(positionP3, stateP3.getFrame(), stateP3.getDate());
370 double delayP3 = model.pathDelay(elevationP3, point, model.getParameters(), stateP3.getDate());
371
372 SpacecraftState stateP4 = shiftState(state, orbitType, angleType, 4 * steps[i], i);
373 final Vector3D positionP4 = stateP4.getPVCoordinates().getPosition();
374 final double elevationP4 = station.getBaseFrame().getElevation(positionP4, stateP4.getFrame(), stateP4.getDate());
375 double delayP4 = model.pathDelay(elevationP4, point, model.getParameters(), stateP4.getDate());
376
377 fillJacobianColumn(refDeriv, i, orbitType, angleType, steps[i],
378 delayM4, delayM3, delayM2, delayM1,
379 delayP1, delayP2, delayP3, delayP4);
380 }
381
382 for (int i = 0; i < 6; i++) {
383 Assert.assertEquals(compDelay[i + 1], refDeriv[0][i], 6.2e-12);
384 }
385 }
386
387 private void fillJacobianColumn(double[][] jacobian, int column,
388 OrbitType orbitType, PositionAngle angleType, double h,
389 double sM4h, double sM3h,
390 double sM2h, double sM1h,
391 double sP1h, double sP2h,
392 double sP3h, double sP4h) {
393 for (int i = 0; i < jacobian.length; ++i) {
394 jacobian[i][column] = ( -3 * (sP4h - sM4h) +
395 32 * (sP3h - sM3h) -
396 168 * (sP2h - sM2h) +
397 672 * (sP1h - sM1h)) / (840 * h);
398 }
399 }
400
401 private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType, PositionAngle angleType,
402 double delta, int column) {
403
404 double[][] array = stateToArray(state, orbitType, angleType, true);
405 array[0][column] += delta;
406
407 return arrayToState(array, orbitType, angleType, state.getFrame(), state.getDate(),
408 state.getMu(), state.getAttitude());
409
410 }
411
412 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType, PositionAngle angleType,
413 boolean withMass) {
414 double[][] array = new double[2][withMass ? 7 : 6];
415 orbitType.mapOrbitToArray(state.getOrbit(), angleType, array[0], array[1]);
416 if (withMass) {
417 array[0][6] = state.getMass();
418 }
419 return array;
420 }
421
422 private SpacecraftState arrayToState(double[][] array, OrbitType orbitType, PositionAngle angleType,
423 Frame frame, AbsoluteDate date, double mu,
424 Attitude attitude) {
425 Orbit orbit = orbitType.mapArrayToOrbit(array[0], array[1], angleType, date, mu, frame);
426 return (array.length > 6) ?
427 new SpacecraftState(orbit, attitude) :
428 new SpacecraftState(orbit, attitude, array[0][6]);
429 }
430
431 }