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.MathArrays;
28 import org.hipparchus.util.Precision;
29 import org.junit.Assert;
30 import org.junit.Before;
31 import org.junit.BeforeClass;
32 import org.junit.Test;
33 import org.orekit.Utils;
34 import org.orekit.attitudes.Attitude;
35 import org.orekit.bodies.FieldGeodeticPoint;
36 import org.orekit.bodies.GeodeticPoint;
37 import org.orekit.bodies.OneAxisEllipsoid;
38 import org.orekit.errors.OrekitException;
39 import org.orekit.estimation.measurements.GroundStation;
40 import org.orekit.frames.Frame;
41 import org.orekit.frames.FramesFactory;
42 import org.orekit.frames.TopocentricFrame;
43 import org.orekit.orbits.FieldKeplerianOrbit;
44 import org.orekit.orbits.FieldOrbit;
45 import org.orekit.orbits.Orbit;
46 import org.orekit.orbits.OrbitType;
47 import org.orekit.orbits.PositionAngle;
48 import org.orekit.propagation.FieldSpacecraftState;
49 import org.orekit.propagation.SpacecraftState;
50 import org.orekit.propagation.numerical.NumericalPropagator;
51 import org.orekit.time.AbsoluteDate;
52 import org.orekit.time.FieldAbsoluteDate;
53 import org.orekit.time.TimeScalesFactory;
54 import org.orekit.utils.Constants;
55 import org.orekit.utils.IERSConventions;
56
57 public class FieldGlobalMappingFunctionModelTest {
58
59 @BeforeClass
60 public static void setUpGlobal() {
61 Utils.setDataRoot("atmosphere");
62 }
63
64 @Before
65 public void setUp() throws OrekitException {
66 Utils.setDataRoot("regular-data:potential/shm-format");
67 }
68
69 @Test
70 public void testMappingFactors() {
71 doTestMappingFactors(Decimal64Field.getInstance());
72 }
73
74 private <T extends CalculusFieldElement<T>> void doTestMappingFactors(final Field<T> field) {
75 final T zero = field.getZero();
76
77
78
79
80
81
82
83
84
85
86
87
88 final FieldAbsoluteDate<T> date = FieldAbsoluteDate.createMJDDate(55055, zero, TimeScalesFactory.getUTC());
89
90 final double latitude = 0.6708665767;
91 final double longitude = -1.393397187;
92 final double height = 844.715;
93 final FieldGeodeticPoint<T> point = new FieldGeodeticPoint<T>(zero.add(latitude), zero.add(longitude), zero.add(height));
94
95 final double elevation = 0.5 * FastMath.PI - 1.278564131;
96 final double expectedHydro = 3.425246;
97 final double expectedWet = 3.449589;
98
99 final MappingFunction model = new GlobalMappingFunctionModel();
100
101 final T[] computedMapping = model.mappingFactors(zero.add(elevation), point, date);
102
103 Assert.assertEquals(expectedHydro, computedMapping[0].getReal(), 1.0e-6);
104 Assert.assertEquals(expectedWet, computedMapping[1].getReal(), 1.0e-6);
105 }
106
107 @Test
108 public void testFixedHeight() {
109 doTestFixedHeight(Decimal64Field.getInstance());
110 }
111
112 private <T extends CalculusFieldElement<T>> void doTestFixedHeight(final Field<T> field) {
113 final T zero = field.getZero();
114 final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
115 MappingFunction model = new GlobalMappingFunctionModel();
116 FieldGeodeticPoint<T> point = new FieldGeodeticPoint<T>(zero.add(FastMath.toRadians(45.0)), zero.add(FastMath.toRadians(45.0)), zero.add(350.0));
117 final T[] lastFactors = MathArrays.buildArray(field, 2);
118 lastFactors[0] = zero.add(Double.MAX_VALUE);
119 lastFactors[1] = zero.add(Double.MAX_VALUE);
120
121
122 for (double elev = 10d; elev < 90d; elev += 8d) {
123 final T[] factors = model.mappingFactors(zero.add(FastMath.toRadians(elev)), point, date);
124 Assert.assertTrue(Precision.compareTo(factors[0].getReal(), lastFactors[0].getReal(), 1.0e-6) < 0);
125 Assert.assertTrue(Precision.compareTo(factors[1].getReal(), lastFactors[1].getReal(), 1.0e-6) < 0);
126 lastFactors[0] = factors[0];
127 lastFactors[1] = factors[1];
128 }
129 }
130
131 @Test
132 public void testMFStateDerivatives() {
133
134
135 final double latitude = FastMath.toRadians(45.0);
136 final double longitude = FastMath.toRadians(45.0);
137 final double height = 0.0;
138 final GeodeticPoint point = new GeodeticPoint(latitude, longitude, height);
139
140 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
141 Constants.WGS84_EARTH_FLATTENING,
142 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
143
144 final TopocentricFrame baseFrame = new TopocentricFrame(earth, point, "topo");
145
146
147 final GroundStation station = new GroundStation(baseFrame);
148
149
150 final MappingFunction model = new GlobalMappingFunctionModel();
151
152
153 final DSFactory factory = new DSFactory(6, 1);
154 final DerivativeStructure a0 = factory.variable(0, 24464560.0);
155 final DerivativeStructure e0 = factory.variable(1, 0.05);
156 final DerivativeStructure i0 = factory.variable(2, 0.122138);
157 final DerivativeStructure pa0 = factory.variable(3, 3.10686);
158 final DerivativeStructure raan0 = factory.variable(4, 1.00681);
159 final DerivativeStructure anomaly0 = factory.variable(5, 0.048363);
160 final Field<DerivativeStructure> field = a0.getField();
161 final DerivativeStructure zero = field.getZero();
162
163
164 final FieldAbsoluteDate<DerivativeStructure> dsDate = new FieldAbsoluteDate<>(field);
165
166 final Frame frame = FramesFactory.getEME2000();
167 final FieldOrbit<DerivativeStructure> dsOrbit = new FieldKeplerianOrbit<>(a0, e0, i0, pa0, raan0, anomaly0,
168 PositionAngle.MEAN, frame,
169 dsDate, zero.add(3.9860047e14));
170
171 final FieldSpacecraftState<DerivativeStructure> dsState = new FieldSpacecraftState<>(dsOrbit);
172
173
174 final FieldVector3D<DerivativeStructure> position = dsState.getPVCoordinates().getPosition();
175 final DerivativeStructure dsElevation = baseFrame.getElevation(position, frame, dsDate);
176
177
178 final FieldGeodeticPoint<DerivativeStructure> dsPoint = new FieldGeodeticPoint<>(zero.add(latitude), zero.add(longitude), zero.add(height));
179 final DerivativeStructure[] factors = model.mappingFactors(dsElevation, dsPoint, dsDate);
180
181 final double[] compMFH = factors[0].getAllDerivatives();
182 final double[] compMFW = factors[1].getAllDerivatives();
183
184
185 final Orbit orbit = dsOrbit.toOrbit();
186 final SpacecraftState state = dsState.toSpacecraftState();
187
188
189 final double[][] refMF = new double[2][6];
190 final OrbitType orbitType = OrbitType.KEPLERIAN;
191 final PositionAngle angleType = PositionAngle.MEAN;
192 double dP = 0.001;
193 double[] steps = NumericalPropagator.tolerances(1000000 * dP, orbit, orbitType)[0];
194 for (int i = 0; i < 6; i++) {
195 SpacecraftState stateM4 = shiftState(state, orbitType, angleType, -4 * steps[i], i);
196 final Vector3D positionM4 = stateM4.getPVCoordinates().getPosition();
197 final double elevationM4 = station.getBaseFrame().getElevation(positionM4, stateM4.getFrame(), stateM4.getDate());
198 double[] delayM4 = model.mappingFactors(elevationM4, point, stateM4.getDate());
199
200 SpacecraftState stateM3 = shiftState(state, orbitType, angleType, -3 * steps[i], i);
201 final Vector3D positionM3 = stateM3.getPVCoordinates().getPosition();
202 final double elevationM3 = station.getBaseFrame().getElevation(positionM3, stateM3.getFrame(), stateM3.getDate());
203 double[] delayM3 = model.mappingFactors(elevationM3, point, stateM3.getDate());
204
205 SpacecraftState stateM2 = shiftState(state, orbitType, angleType, -2 * steps[i], i);
206 final Vector3D positionM2 = stateM2.getPVCoordinates().getPosition();
207 final double elevationM2 = station.getBaseFrame().getElevation(positionM2, stateM2.getFrame(), stateM2.getDate());
208 double[] delayM2 = model.mappingFactors(elevationM2, point, stateM2.getDate());
209
210 SpacecraftState stateM1 = shiftState(state, orbitType, angleType, -1 * steps[i], i);
211 final Vector3D positionM1 = stateM1.getPVCoordinates().getPosition();
212 final double elevationM1 = station.getBaseFrame().getElevation(positionM1, stateM1.getFrame(), stateM1.getDate());
213 double[] delayM1 = model.mappingFactors(elevationM1, point, stateM1.getDate());
214
215 SpacecraftState stateP1 = shiftState(state, orbitType, angleType, 1 * steps[i], i);
216 final Vector3D positionP1 = stateP1.getPVCoordinates().getPosition();
217 final double elevationP1 = station.getBaseFrame().getElevation(positionP1, stateP1.getFrame(), stateP1.getDate());
218 double[] delayP1 = model.mappingFactors(elevationP1, point, stateP1.getDate());
219
220 SpacecraftState stateP2 = shiftState(state, orbitType, angleType, 2 * steps[i], i);
221 final Vector3D positionP2 = stateP2.getPVCoordinates().getPosition();
222 final double elevationP2 = station.getBaseFrame().getElevation(positionP2, stateP2.getFrame(), stateP2.getDate());
223 double[] delayP2 = model.mappingFactors(elevationP2, point, stateP2.getDate());
224
225 SpacecraftState stateP3 = shiftState(state, orbitType, angleType, 3 * steps[i], i);
226 final Vector3D positionP3 = stateP3.getPVCoordinates().getPosition();
227 final double elevationP3 = station.getBaseFrame().getElevation(positionP3, stateP3.getFrame(), stateP3.getDate());
228 double[] delayP3 = model.mappingFactors(elevationP3, point, stateP3.getDate());
229
230 SpacecraftState stateP4 = shiftState(state, orbitType, angleType, 4 * steps[i], i);
231 final Vector3D positionP4 = stateP4.getPVCoordinates().getPosition();
232 final double elevationP4 = station.getBaseFrame().getElevation(positionP4, stateP4.getFrame(), stateP4.getDate());
233 double[] delayP4 = model.mappingFactors(elevationP4, point, stateP4.getDate());
234
235 fillJacobianColumn(refMF, i, orbitType, angleType, steps[i],
236 delayM4, delayM3, delayM2, delayM1,
237 delayP1, delayP2, delayP3, delayP4);
238 }
239
240 for (int i = 0; i < 6; i++) {
241 Assert.assertEquals(compMFH[i + 1], refMF[0][i], 2.1e-11);
242 Assert.assertEquals(compMFW[i + 1], refMF[1][i], 1.7e-11);
243 }
244 }
245
246 private void fillJacobianColumn(double[][] jacobian, int column,
247 OrbitType orbitType, PositionAngle angleType, double h,
248 double[] sM4h, double[] sM3h,
249 double[] sM2h, double[] sM1h,
250 double[] sP1h, double[] sP2h,
251 double[] sP3h, double[] sP4h) {
252 for (int i = 0; i < jacobian.length; ++i) {
253 jacobian[i][column] = ( -3 * (sP4h[i] - sM4h[i]) +
254 32 * (sP3h[i] - sM3h[i]) -
255 168 * (sP2h[i] - sM2h[i]) +
256 672 * (sP1h[i] - sM1h[i])) / (840 * h);
257 }
258 }
259
260 private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType, PositionAngle angleType,
261 double delta, int column) {
262
263 double[][] array = stateToArray(state, orbitType, angleType, true);
264 array[0][column] += delta;
265
266 return arrayToState(array, orbitType, angleType, state.getFrame(), state.getDate(),
267 state.getMu(), state.getAttitude());
268
269 }
270
271 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType, PositionAngle angleType,
272 boolean withMass) {
273 double[][] array = new double[2][withMass ? 7 : 6];
274 orbitType.mapOrbitToArray(state.getOrbit(), angleType, array[0], array[1]);
275 if (withMass) {
276 array[0][6] = state.getMass();
277 }
278 return array;
279 }
280
281 private SpacecraftState arrayToState(double[][] array, OrbitType orbitType, PositionAngle angleType,
282 Frame frame, AbsoluteDate date, double mu,
283 Attitude attitude) {
284 Orbit orbit = orbitType.mapArrayToOrbit(array[0], array[1], angleType, date, mu, frame);
285 return (array.length > 6) ?
286 new SpacecraftState(orbit, attitude) :
287 new SpacecraftState(orbit, attitude, array[0][6]);
288 }
289
290 }