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 FieldViennaOneModelTest {
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 final T zero = field.getZero();
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100 final FieldAbsoluteDate<T> date = FieldAbsoluteDate.createMJDDate(55055, zero, TimeScalesFactory.getUTC());
101
102 final double latitude = FastMath.toRadians(38.0);
103 final double longitude = FastMath.toRadians(280.0);
104 final double height = 824.17;
105
106 final double elevation = 0.5 * FastMath.PI - 1.278564131;
107 final double expectedHydro = 3.425088;
108 final double expectedWet = 3.448300;
109
110 final double[] a = { 0.00127683, 0.00060955 };
111 final double[] z = {2.0966, 0.2140};
112
113 final FieldGeodeticPoint<T> point = new FieldGeodeticPoint<>(zero.add(latitude), zero.add(longitude), zero.add(height));
114 final ViennaOneModel model = new ViennaOneModel(a, z);
115
116 final T[] computedMapping = model.mappingFactors(zero.add(elevation), point, date);
117
118 Assert.assertEquals(expectedHydro, computedMapping[0].getReal(), 4.1e-6);
119 Assert.assertEquals(expectedWet, computedMapping[1].getReal(), 1.0e-6);
120 }
121
122 @Test
123 public void testDelay() {
124 doTestDelay(Decimal64Field.getInstance());
125 }
126
127 private <T extends CalculusFieldElement<T>> void doTestDelay(final Field<T> field) {
128 final T zero = field.getZero();
129 final double elevation = 10d;
130 final double height = 100d;
131 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
132 final double[] a = { 0.00127683, 0.00060955 };
133 final double[] z = {2.0966, 0.2140};
134 final FieldGeodeticPoint<T> point = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(45.0)), zero.add(FastMath.toRadians(45.0)), zero.add(height));
135 ViennaOneModel model = new ViennaOneModel(a, z);
136 final T path = model.pathDelay(zero.add(FastMath.toRadians(elevation)), point, model.getParameters(field), date);
137 Assert.assertTrue(Precision.compareTo(path.getReal(), 20d, epsilon) < 0);
138 Assert.assertTrue(Precision.compareTo(path.getReal(), 0d, epsilon) > 0);
139 }
140
141 @Test
142 public void testFixedHeight() {
143 doTestFixedHeight(Decimal64Field.getInstance());
144 }
145
146 private <T extends CalculusFieldElement<T>> void doTestFixedHeight(final Field<T> field) {
147 final T zero = field.getZero();
148 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
149 final double[] a = { 0.00127683, 0.00060955 };
150 final double[] z = {2.0966, 0.2140};
151 final FieldGeodeticPoint<T> point = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(45.0)), zero.add(FastMath.toRadians(45.0)), zero.add(350.0));
152 ViennaOneModel model = new ViennaOneModel(a, z);
153 T lastDelay = zero.add(Double.MAX_VALUE);
154
155 for (double elev = 10d; elev < 90d; elev += 8d) {
156 final T delay = model.pathDelay(zero.add(FastMath.toRadians(elev)), point, model.getParameters(field), date);
157 Assert.assertTrue(Precision.compareTo(delay.getReal(), lastDelay.getReal(), epsilon) < 0);
158 lastDelay = delay;
159 }
160 }
161
162 @Test
163 public void testDelayStateDerivatives() {
164
165
166 final double latitude = FastMath.toRadians(45.0);
167 final double longitude = FastMath.toRadians(45.0);
168 final double height = 0.0;
169 final GeodeticPoint point = new GeodeticPoint(latitude, longitude, height);
170
171 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
172 Constants.WGS84_EARTH_FLATTENING,
173 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
174
175 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "topo");
176
177
178 final GroundStation station = new GroundStation(baseFrame);
179
180
181 final double[] a = { 0.00127683, 0.00060955 };
182 final double[] z = {2.0966, 0.2140};
183 final DiscreteTroposphericModel model = new ViennaOneModel(a, z);
184
185
186 final DSFactory factory = new DSFactory(6, 1);
187 final DerivativeStructure a0 = factory.variable(0, 24464560.0);
188 final DerivativeStructure e0 = factory.variable(1, 0.05);
189 final DerivativeStructure i0 = factory.variable(2, 0.122138);
190 final DerivativeStructure pa0 = factory.variable(3, 3.10686);
191 final DerivativeStructure raan0 = factory.variable(4, 1.00681);
192 final DerivativeStructure anomaly0 = factory.variable(5, 0.048363);
193 final Field<DerivativeStructure> field = a0.getField();
194 final DerivativeStructure zero = field.getZero();
195
196
197 final FieldAbsoluteDate<DerivativeStructure> dsDate = new FieldAbsoluteDate<>(field, 2018, 11, 19, 18, 0, 0.0,
198 TimeScalesFactory.getUTC());
199
200 final Frame frame = FramesFactory.getEME2000();
201 final FieldOrbit<DerivativeStructure> dsOrbit = new FieldKeplerianOrbit<>(a0, e0, i0, pa0, raan0, anomaly0,
202 PositionAngle.MEAN, frame,
203 dsDate, zero.add(3.9860047e14));
204
205 final FieldSpacecraftState<DerivativeStructure> dsState = new FieldSpacecraftState<>(dsOrbit);
206
207
208 final FieldVector3D<DerivativeStructure> position = dsState.getPVCoordinates().getPosition();
209 final DerivativeStructure dsElevation = baseFrame.getElevation(position, frame, dsDate);
210
211
212 final FieldGeodeticPoint<DerivativeStructure> dsPoint = new FieldGeodeticPoint<>(zero.add(latitude), zero.add(longitude), zero.add(height));
213 final DerivativeStructure delay = model.pathDelay(dsElevation, dsPoint, model.getParameters(field), dsDate);
214
215 final double[] compDelay = delay.getAllDerivatives();
216
217
218 final Orbit orbit = dsOrbit.toOrbit();
219 final SpacecraftState state = dsState.toSpacecraftState();
220
221
222 final double[][] refDeriv = new double[1][6];
223 final OrbitType orbitType = OrbitType.KEPLERIAN;
224 final PositionAngle angleType = PositionAngle.MEAN;
225 double dP = 0.001;
226 double[] steps = NumericalPropagator.tolerances(1000000 * dP, orbit, orbitType)[0];
227 for (int i = 0; i < 6; i++) {
228 SpacecraftState stateM4 = shiftState(state, orbitType, angleType, -4 * steps[i], i);
229 final Vector3D positionM4 = stateM4.getPVCoordinates().getPosition();
230 final double elevationM4 = station.getBaseFrame().getElevation(positionM4, stateM4.getFrame(), stateM4.getDate());
231 double delayM4 = model.pathDelay(elevationM4, point, model.getParameters(), stateM4.getDate());
232
233 SpacecraftState stateM3 = shiftState(state, orbitType, angleType, -3 * steps[i], i);
234 final Vector3D positionM3 = stateM3.getPVCoordinates().getPosition();
235 final double elevationM3 = station.getBaseFrame().getElevation(positionM3, stateM3.getFrame(), stateM3.getDate());
236 double delayM3 = model.pathDelay(elevationM3, point, model.getParameters(), stateM3.getDate());
237
238 SpacecraftState stateM2 = shiftState(state, orbitType, angleType, -2 * steps[i], i);
239 final Vector3D positionM2 = stateM2.getPVCoordinates().getPosition();
240 final double elevationM2 = station.getBaseFrame().getElevation(positionM2, stateM2.getFrame(), stateM2.getDate());
241 double delayM2 = model.pathDelay(elevationM2, point, model.getParameters(), stateM2.getDate());
242
243 SpacecraftState stateM1 = shiftState(state, orbitType, angleType, -1 * steps[i], i);
244 final Vector3D positionM1 = stateM1.getPVCoordinates().getPosition();
245 final double elevationM1 = station.getBaseFrame().getElevation(positionM1, stateM1.getFrame(), stateM1.getDate());
246 double delayM1 = model.pathDelay(elevationM1, point, model.getParameters(), stateM1.getDate());
247
248 SpacecraftState stateP1 = shiftState(state, orbitType, angleType, 1 * steps[i], i);
249 final Vector3D positionP1 = stateP1.getPVCoordinates().getPosition();
250 final double elevationP1 = station.getBaseFrame().getElevation(positionP1, stateP1.getFrame(), stateP1.getDate());
251 double delayP1 = model.pathDelay(elevationP1, point, model.getParameters(), stateP1.getDate());
252
253 SpacecraftState stateP2 = shiftState(state, orbitType, angleType, 2 * steps[i], i);
254 final Vector3D positionP2 = stateP2.getPVCoordinates().getPosition();
255 final double elevationP2 = station.getBaseFrame().getElevation(positionP2, stateP2.getFrame(), stateP2.getDate());
256 double delayP2 = model.pathDelay(elevationP2, point, model.getParameters(), stateP2.getDate());
257
258 SpacecraftState stateP3 = shiftState(state, orbitType, angleType, 3 * steps[i], i);
259 final Vector3D positionP3 = stateP3.getPVCoordinates().getPosition();
260 final double elevationP3 = station.getBaseFrame().getElevation(positionP3, stateP3.getFrame(), stateP3.getDate());
261 double delayP3 = model.pathDelay(elevationP3, point, model.getParameters(), stateP3.getDate());
262
263 SpacecraftState stateP4 = shiftState(state, orbitType, angleType, 4 * steps[i], i);
264 final Vector3D positionP4 = stateP4.getPVCoordinates().getPosition();
265 final double elevationP4 = station.getBaseFrame().getElevation(positionP4, stateP4.getFrame(), stateP4.getDate());
266 double delayP4 = model.pathDelay(elevationP4, point, model.getParameters(), stateP4.getDate());
267
268 fillJacobianColumn(refDeriv, i, orbitType, angleType, steps[i],
269 delayM4, delayM3, delayM2, delayM1,
270 delayP1, delayP2, delayP3, delayP4);
271 }
272
273 for (int i = 0; i < 6; i++) {
274 Assert.assertEquals(compDelay[i + 1], refDeriv[0][i], 3.0e-11);
275 }
276 }
277
278 private void fillJacobianColumn(double[][] jacobian, int column,
279 OrbitType orbitType, PositionAngle angleType, double h,
280 double sM4h, double sM3h,
281 double sM2h, double sM1h,
282 double sP1h, double sP2h,
283 double sP3h, double sP4h) {
284 for (int i = 0; i < jacobian.length; ++i) {
285 jacobian[i][column] = ( -3 * (sP4h - sM4h) +
286 32 * (sP3h - sM3h) -
287 168 * (sP2h - sM2h) +
288 672 * (sP1h - sM1h)) / (840 * h);
289 }
290 }
291
292 private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType, PositionAngle angleType,
293 double delta, int column) {
294
295 double[][] array = stateToArray(state, orbitType, angleType, true);
296 array[0][column] += delta;
297
298 return arrayToState(array, orbitType, angleType, state.getFrame(), state.getDate(),
299 state.getMu(), state.getAttitude());
300
301 }
302
303 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType, PositionAngle angleType,
304 boolean withMass) {
305 double[][] array = new double[2][withMass ? 7 : 6];
306 orbitType.mapOrbitToArray(state.getOrbit(), angleType, array[0], array[1]);
307 if (withMass) {
308 array[0][6] = state.getMass();
309 }
310 return array;
311 }
312
313 private SpacecraftState arrayToState(double[][] array, OrbitType orbitType, PositionAngle angleType,
314 Frame frame, AbsoluteDate date, double mu,
315 Attitude attitude) {
316 Orbit orbit = orbitType.mapArrayToOrbit(array[0], array[1], angleType, date, mu, frame);
317 return (array.length > 6) ?
318 new SpacecraftState(orbit, attitude) :
319 new SpacecraftState(orbit, attitude, array[0][6]);
320 }
321
322 }