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