1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces;
18
19
20 import java.util.List;
21 import java.util.concurrent.atomic.AtomicReference;
22
23 import org.hipparchus.Field;
24 import org.hipparchus.analysis.differentiation.DSFactory;
25 import org.hipparchus.analysis.differentiation.DerivativeStructure;
26 import org.hipparchus.analysis.differentiation.Gradient;
27 import org.hipparchus.analysis.differentiation.GradientField;
28 import org.hipparchus.geometry.euclidean.threed.FieldRotation;
29 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
30 import org.hipparchus.geometry.euclidean.threed.Vector3D;
31 import org.hipparchus.linear.RealMatrix;
32 import org.hipparchus.random.GaussianRandomGenerator;
33 import org.hipparchus.random.RandomGenerator;
34 import org.hipparchus.random.UncorrelatedRandomVectorGenerator;
35 import org.hipparchus.random.Well19937a;
36 import org.hipparchus.util.FastMath;
37 import org.hipparchus.util.Precision;
38 import org.junit.Assert;
39 import org.orekit.attitudes.AttitudeProvider;
40 import org.orekit.attitudes.FieldAttitude;
41 import org.orekit.orbits.FieldCartesianOrbit;
42 import org.orekit.orbits.FieldKeplerianOrbit;
43 import org.orekit.orbits.KeplerianOrbit;
44 import org.orekit.orbits.OrbitType;
45 import org.orekit.orbits.PositionAngle;
46 import org.orekit.propagation.FieldSpacecraftState;
47 import org.orekit.propagation.MatricesHarvester;
48 import org.orekit.propagation.SpacecraftState;
49 import org.orekit.propagation.numerical.FieldNumericalPropagator;
50 import org.orekit.propagation.numerical.NumericalPropagator;
51 import org.orekit.propagation.sampling.OrekitStepHandler;
52 import org.orekit.propagation.sampling.OrekitStepInterpolator;
53 import org.orekit.time.AbsoluteDate;
54 import org.orekit.time.FieldAbsoluteDate;
55 import org.orekit.utils.Constants;
56 import org.orekit.utils.Differentiation;
57 import org.orekit.utils.FieldPVCoordinates;
58 import org.orekit.utils.PVCoordinates;
59 import org.orekit.utils.ParameterDriver;
60 import org.orekit.utils.TimeStampedFieldAngularCoordinates;
61 import org.orekit.utils.TimeStampedFieldPVCoordinates;
62
63
64 public abstract class AbstractForceModelTest {
65
66 protected void checkParameterDerivative(SpacecraftState state,
67 ForceModel forceModel, String name,
68 double hFactor, double tol) {
69
70 final DSFactory factory11 = new DSFactory(1, 1);
71 final Field<DerivativeStructure> field = factory11.getDerivativeField();
72 final FieldSpacecraftState<DerivativeStructure> stateF = new FieldSpacecraftState<DerivativeStructure>(field, state);
73 final List<ParameterDriver> drivers = forceModel.getParametersDrivers();
74 final DerivativeStructure[] parametersDS = new DerivativeStructure[drivers.size()];
75 for (int i = 0; i < parametersDS.length; ++i) {
76 if (drivers.get(i).getName().equals(name)) {
77 parametersDS[i] = factory11.variable(0, drivers.get(i).getValue());
78 } else {
79 parametersDS[i] = factory11.constant(drivers.get(i).getValue());
80 }
81 }
82 FieldVector3D<DerivativeStructure> accDer = forceModel.acceleration(stateF, parametersDS);
83 Vector3D derivative = new Vector3D(accDer.getX().getPartialDerivative(1),
84 accDer.getY().getPartialDerivative(1),
85 accDer.getZ().getPartialDerivative(1));
86
87 int selected = -1;
88 final double[] parameters = new double[drivers.size()];
89 for (int i = 0; i < drivers.size(); ++i) {
90 parameters[i] = drivers.get(i).getValue();
91 if (drivers.get(i).getName().equals(name)) {
92 selected = i;
93 }
94 }
95 double p0 = parameters[selected];
96 double hParam = hFactor * p0;
97 drivers.get(selected).setValue(p0 - 1 * hParam);
98 parameters[selected] = drivers.get(selected).getValue();
99 Assert.assertEquals(p0 - 1 * hParam, parameters[selected], 1.0e-10);
100 final Vector3D gammaM1h = forceModel.acceleration(state, parameters);
101 drivers.get(selected).setValue(p0 + 1 * hParam);
102 parameters[selected] = drivers.get(selected).getValue();
103 Assert.assertEquals(p0 + 1 * hParam, parameters[selected], 1.0e-10);
104 final Vector3D gammaP1h = forceModel.acceleration(state, parameters);
105 drivers.get(selected).setValue(p0);
106
107 final Vector3D reference = new Vector3D( 1 / (2 * hParam), gammaP1h.subtract(gammaM1h));
108 final Vector3D delta = derivative.subtract(reference);
109 Assert.assertEquals(0., delta.getNorm(), tol * reference.getNorm());
110
111 }
112
113 protected void checkParameterDerivativeGradient(SpacecraftState state,
114 ForceModel forceModel, String name,
115 double hFactor, double tol) {
116
117 final int freeParameters = 1;
118 final Field<Gradient> field = GradientField.getField(freeParameters);
119 final FieldSpacecraftState<Gradient> stateF = new FieldSpacecraftState<>(field, state);
120 final List<ParameterDriver> drivers = forceModel.getParametersDrivers();
121 final Gradient[] parametersDS = new Gradient[drivers.size()];
122 for (int i = 0; i < parametersDS.length; ++i) {
123 if (drivers.get(i).getName().equals(name)) {
124 parametersDS[i] = Gradient.variable(freeParameters, 0, drivers.get(i).getValue());
125 } else {
126 parametersDS[i] = Gradient.constant(freeParameters, drivers.get(i).getValue());
127 }
128 }
129 FieldVector3D<Gradient> accDer = forceModel.acceleration(stateF, parametersDS);
130 Vector3D derivative = new Vector3D(accDer.getX().getPartialDerivative(0),
131 accDer.getY().getPartialDerivative(0),
132 accDer.getZ().getPartialDerivative(0));
133
134 int selected = -1;
135 final double[] parameters = new double[drivers.size()];
136 for (int i = 0; i < drivers.size(); ++i) {
137 parameters[i] = drivers.get(i).getValue();
138 if (drivers.get(i).getName().equals(name)) {
139 selected = i;
140 }
141 }
142 double p0 = parameters[selected];
143 double hParam = hFactor * p0;
144 drivers.get(selected).setValue(p0 - 1 * hParam);
145 parameters[selected] = drivers.get(selected).getValue();
146 Assert.assertEquals(p0 - 1 * hParam, parameters[selected], 1.0e-10);
147 final Vector3D gammaM1h = forceModel.acceleration(state, parameters);
148 drivers.get(selected).setValue(p0 + 1 * hParam);
149 parameters[selected] = drivers.get(selected).getValue();
150 Assert.assertEquals(p0 + 1 * hParam, parameters[selected], 1.0e-10);
151 final Vector3D gammaP1h = forceModel.acceleration(state, parameters);
152 drivers.get(selected).setValue(p0);
153
154 final Vector3D reference = new Vector3D( 1 / (2 * hParam), gammaP1h.subtract(gammaM1h));
155 final Vector3D delta = derivative.subtract(reference);
156 Assert.assertEquals(0, delta.getNorm(), tol * reference.getNorm());
157
158 }
159
160 protected FieldSpacecraftState<DerivativeStructure> toDS(final SpacecraftState state,
161 final AttitudeProvider attitudeProvider) {
162
163 final Vector3D p = state.getPVCoordinates().getPosition();
164 final Vector3D v = state.getPVCoordinates().getVelocity();
165 final Vector3D a = state.getPVCoordinates().getAcceleration();
166 DSFactory factory = new DSFactory(6, 1);
167 Field<DerivativeStructure> field = factory.getDerivativeField();
168 final FieldAbsoluteDate<DerivativeStructure> fDate = new FieldAbsoluteDate<>(field, state.getDate());
169 final TimeStampedFieldPVCoordinates<DerivativeStructure> fPVA =
170 new TimeStampedFieldPVCoordinates<>(fDate,
171 new FieldVector3D<>(factory.variable(0, p.getX()),
172 factory.variable(1, p.getY()),
173 factory.variable(2, p.getZ())),
174 new FieldVector3D<>(factory.variable(3, v.getX()),
175 factory.variable(4, v.getY()),
176 factory.variable(5, v.getZ())),
177 new FieldVector3D<>(factory.constant(a.getX()),
178 factory.constant(a.getY()),
179 factory.constant(a.getZ())));
180 final FieldCartesianOrbit<DerivativeStructure> orbit =
181 new FieldCartesianOrbit<>(fPVA, state.getFrame(), field.getZero().add(state.getMu()));
182 final FieldAttitude<DerivativeStructure> attitude =
183 attitudeProvider.getAttitude(orbit, orbit.getDate(), orbit.getFrame());
184
185 return new FieldSpacecraftState<>(orbit, attitude, field.getZero().add(state.getMass()));
186
187 }
188
189 protected FieldSpacecraftState<Gradient> toGradient(final SpacecraftState state,
190 final AttitudeProvider attitudeProvider) {
191
192 final Vector3D p = state.getPVCoordinates().getPosition();
193 final Vector3D v = state.getPVCoordinates().getVelocity();
194 final Vector3D a = state.getPVCoordinates().getAcceleration();
195 final int freeParameters = 6;
196 Field<Gradient> field = GradientField.getField(freeParameters);
197 final FieldAbsoluteDate<Gradient> fDate = new FieldAbsoluteDate<>(field, state.getDate());
198 final TimeStampedFieldPVCoordinates<Gradient> fPVA =
199 new TimeStampedFieldPVCoordinates<>(fDate,
200 new FieldVector3D<>(Gradient.variable(freeParameters, 0, p.getX()),
201 Gradient.variable(freeParameters, 1, p.getY()),
202 Gradient.variable(freeParameters, 2, p.getZ())),
203 new FieldVector3D<>(Gradient.variable(freeParameters, 3, v.getX()),
204 Gradient.variable(freeParameters, 4, v.getY()),
205 Gradient.variable(freeParameters, 5, v.getZ())),
206 new FieldVector3D<>(Gradient.constant(freeParameters, a.getX()),
207 Gradient.constant(freeParameters, a.getY()),
208 Gradient.constant(freeParameters, a.getZ())));
209 final FieldCartesianOrbit<Gradient> orbit =
210 new FieldCartesianOrbit<>(fPVA, state.getFrame(), field.getZero().add(state.getMu()));
211 final FieldAttitude<Gradient> attitude =
212 attitudeProvider.getAttitude(orbit, orbit.getDate(), orbit.getFrame());
213
214 return new FieldSpacecraftState<>(orbit, attitude, field.getZero().add(state.getMass()));
215
216 }
217
218 protected void checkStateJacobianVsFiniteDifferences(final SpacecraftState state0, final ForceModel forceModel,
219 final AttitudeProvider provider, final double dP,
220 final double checkTolerance, final boolean print) {
221
222 double[][] finiteDifferencesJacobian =
223 Differentiation.differentiate(state -> forceModel.acceleration(state, forceModel.getParameters()).toArray(),
224 3, provider, OrbitType.CARTESIAN, PositionAngle.MEAN,
225 dP, 5).
226 value(state0);
227
228 DSFactory factory = new DSFactory(6, 1);
229 Field<DerivativeStructure> field = factory.getDerivativeField();
230 final FieldAbsoluteDate<DerivativeStructure> fDate = new FieldAbsoluteDate<>(field, state0.getDate());
231 final Vector3D p = state0.getPVCoordinates().getPosition();
232 final Vector3D v = state0.getPVCoordinates().getVelocity();
233 final Vector3D a = state0.getPVCoordinates().getAcceleration();
234 final TimeStampedFieldPVCoordinates<DerivativeStructure> fPVA =
235 new TimeStampedFieldPVCoordinates<>(fDate,
236 new FieldVector3D<>(factory.variable(0, p.getX()),
237 factory.variable(1, p.getY()),
238 factory.variable(2, p.getZ())),
239 new FieldVector3D<>(factory.variable(3, v.getX()),
240 factory.variable(4, v.getY()),
241 factory.variable(5, v.getZ())),
242 new FieldVector3D<>(factory.constant(a.getX()),
243 factory.constant(a.getY()),
244 factory.constant(a.getZ())));
245 final TimeStampedFieldAngularCoordinates<DerivativeStructure> fAC =
246 new TimeStampedFieldAngularCoordinates<>(fDate,
247 new FieldRotation<>(field, state0.getAttitude().getRotation()),
248 new FieldVector3D<>(field, state0.getAttitude().getSpin()),
249 new FieldVector3D<>(field, state0.getAttitude().getRotationAcceleration()));
250 final FieldSpacecraftState<DerivativeStructure> fState =
251 new FieldSpacecraftState<>(new FieldCartesianOrbit<>(fPVA, state0.getFrame(), field.getZero().add(state0.getMu())),
252 new FieldAttitude<>(state0.getFrame(), fAC),
253 field.getZero().add(state0.getMass()));
254 FieldVector3D<DerivativeStructure> dsJacobian = forceModel.acceleration(fState,
255 forceModel.getParameters(fState.getDate().getField()));
256
257 Vector3D dFdPXRef = new Vector3D(finiteDifferencesJacobian[0][0],
258 finiteDifferencesJacobian[1][0],
259 finiteDifferencesJacobian[2][0]);
260 Vector3D dFdPXRes = new Vector3D(dsJacobian.getX().getPartialDerivative(1, 0, 0, 0, 0, 0),
261 dsJacobian.getY().getPartialDerivative(1, 0, 0, 0, 0, 0),
262 dsJacobian.getZ().getPartialDerivative(1, 0, 0, 0, 0, 0));
263 Vector3D dFdPYRef = new Vector3D(finiteDifferencesJacobian[0][1],
264 finiteDifferencesJacobian[1][1],
265 finiteDifferencesJacobian[2][1]);
266 Vector3D dFdPYRes = new Vector3D(dsJacobian.getX().getPartialDerivative(0, 1, 0, 0, 0, 0),
267 dsJacobian.getY().getPartialDerivative(0, 1, 0, 0, 0, 0),
268 dsJacobian.getZ().getPartialDerivative(0, 1, 0, 0, 0, 0));
269 Vector3D dFdPZRef = new Vector3D(finiteDifferencesJacobian[0][2],
270 finiteDifferencesJacobian[1][2],
271 finiteDifferencesJacobian[2][2]);
272 Vector3D dFdPZRes = new Vector3D(dsJacobian.getX().getPartialDerivative(0, 0, 1, 0, 0, 0),
273 dsJacobian.getY().getPartialDerivative(0, 0, 1, 0, 0, 0),
274 dsJacobian.getZ().getPartialDerivative(0, 0, 1, 0, 0, 0));
275 Vector3D dFdVXRef = new Vector3D(finiteDifferencesJacobian[0][3],
276 finiteDifferencesJacobian[1][3],
277 finiteDifferencesJacobian[2][3]);
278 Vector3D dFdVXRes = new Vector3D(dsJacobian.getX().getPartialDerivative(0, 0, 0, 1, 0, 0),
279 dsJacobian.getY().getPartialDerivative(0, 0, 0, 1, 0, 0),
280 dsJacobian.getZ().getPartialDerivative(0, 0, 0, 1, 0, 0));
281 Vector3D dFdVYRef = new Vector3D(finiteDifferencesJacobian[0][4],
282 finiteDifferencesJacobian[1][4],
283 finiteDifferencesJacobian[2][4]);
284 Vector3D dFdVYRes = new Vector3D(dsJacobian.getX().getPartialDerivative(0, 0, 0, 0, 1, 0),
285 dsJacobian.getY().getPartialDerivative(0, 0, 0, 0, 1, 0),
286 dsJacobian.getZ().getPartialDerivative(0, 0, 0, 0, 1, 0));
287 Vector3D dFdVZRef = new Vector3D(finiteDifferencesJacobian[0][5],
288 finiteDifferencesJacobian[1][5],
289 finiteDifferencesJacobian[2][5]);
290 Vector3D dFdVZRes = new Vector3D(dsJacobian.getX().getPartialDerivative(0, 0, 0, 0, 0, 1),
291 dsJacobian.getY().getPartialDerivative(0, 0, 0, 0, 0, 1),
292 dsJacobian.getZ().getPartialDerivative(0, 0, 0, 0, 0, 1));
293 if (print) {
294 System.out.println("dF/dPX ref: " + dFdPXRef.getX() + " " + dFdPXRef.getY() + " " + dFdPXRef.getZ());
295 System.out.println("dF/dPX res: " + dFdPXRes.getX() + " " + dFdPXRes.getY() + " " + dFdPXRes.getZ());
296 System.out.println("dF/dPY ref: " + dFdPYRef.getX() + " " + dFdPYRef.getY() + " " + dFdPYRef.getZ());
297 System.out.println("dF/dPY res: " + dFdPYRes.getX() + " " + dFdPYRes.getY() + " " + dFdPYRes.getZ());
298 System.out.println("dF/dPZ ref: " + dFdPZRef.getX() + " " + dFdPZRef.getY() + " " + dFdPZRef.getZ());
299 System.out.println("dF/dPZ res: " + dFdPZRes.getX() + " " + dFdPZRes.getY() + " " + dFdPZRes.getZ());
300 System.out.println("dF/dPX ref norm: " + dFdPXRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPXRef, dFdPXRes) + ", rel error: " + (Vector3D.distance(dFdPXRef, dFdPXRes) / dFdPXRef.getNorm()));
301 System.out.println("dF/dPY ref norm: " + dFdPYRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPYRef, dFdPYRes) + ", rel error: " + (Vector3D.distance(dFdPYRef, dFdPYRes) / dFdPYRef.getNorm()));
302 System.out.println("dF/dPZ ref norm: " + dFdPZRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPZRef, dFdPZRes) + ", rel error: " + (Vector3D.distance(dFdPZRef, dFdPZRes) / dFdPZRef.getNorm()));
303 System.out.println("dF/dVX ref norm: " + dFdVXRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVXRef, dFdVXRes) + ", rel error: " + (Vector3D.distance(dFdVXRef, dFdVXRes) / dFdVXRef.getNorm()));
304 System.out.println("dF/dVY ref norm: " + dFdVYRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVYRef, dFdVYRes) + ", rel error: " + (Vector3D.distance(dFdVYRef, dFdVYRes) / dFdVYRef.getNorm()));
305 System.out.println("dF/dVZ ref norm: " + dFdVZRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVZRef, dFdVZRes) + ", rel error: " + (Vector3D.distance(dFdVZRef, dFdVZRes) / dFdVZRef.getNorm()));
306 }
307 checkdFdP(dFdPXRef, dFdPXRes, checkTolerance);
308 checkdFdP(dFdPYRef, dFdPYRes, checkTolerance);
309 checkdFdP(dFdPZRef, dFdPZRes, checkTolerance);
310 checkdFdP(dFdVXRef, dFdVXRes, checkTolerance);
311 checkdFdP(dFdVYRef, dFdVYRes, checkTolerance);
312 checkdFdP(dFdVZRef, dFdVZRes, checkTolerance);
313 }
314
315 protected void checkStateJacobianVsFiniteDifferencesGradient(final SpacecraftState state0, final ForceModel forceModel,
316 final AttitudeProvider provider, final double dP,
317 final double checkTolerance, final boolean print) {
318
319 double[][] finiteDifferencesJacobian =
320 Differentiation.differentiate(state -> forceModel.acceleration(state, forceModel.getParameters()).toArray(),
321 3, provider, OrbitType.CARTESIAN, PositionAngle.MEAN,
322 dP, 5).
323 value(state0);
324
325 final int freePrameters = 6;
326 Field<Gradient> field = GradientField.getField(freePrameters);
327 final FieldAbsoluteDate<Gradient> fDate = new FieldAbsoluteDate<>(field, state0.getDate());
328 final Vector3D p = state0.getPVCoordinates().getPosition();
329 final Vector3D v = state0.getPVCoordinates().getVelocity();
330 final Vector3D a = state0.getPVCoordinates().getAcceleration();
331 final TimeStampedFieldPVCoordinates<Gradient> fPVA =
332 new TimeStampedFieldPVCoordinates<>(fDate,
333 new FieldVector3D<>(Gradient.variable(freePrameters, 0, p.getX()),
334 Gradient.variable(freePrameters, 1, p.getY()),
335 Gradient.variable(freePrameters, 2, p.getZ())),
336 new FieldVector3D<>(Gradient.variable(freePrameters, 3, v.getX()),
337 Gradient.variable(freePrameters, 4, v.getY()),
338 Gradient.variable(freePrameters, 5, v.getZ())),
339 new FieldVector3D<>(Gradient.constant(freePrameters, a.getX()),
340 Gradient.constant(freePrameters, a.getY()),
341 Gradient.constant(freePrameters, a.getZ())));
342 final TimeStampedFieldAngularCoordinates<Gradient> fAC =
343 new TimeStampedFieldAngularCoordinates<>(fDate,
344 new FieldRotation<>(field, state0.getAttitude().getRotation()),
345 new FieldVector3D<>(field, state0.getAttitude().getSpin()),
346 new FieldVector3D<>(field, state0.getAttitude().getRotationAcceleration()));
347 final FieldSpacecraftState<Gradient> fState =
348 new FieldSpacecraftState<>(new FieldCartesianOrbit<>(fPVA, state0.getFrame(), field.getZero().add(state0.getMu())),
349 new FieldAttitude<>(state0.getFrame(), fAC),
350 field.getZero().add(state0.getMass()));
351 FieldVector3D<Gradient> gJacobian = forceModel.acceleration(fState,
352 forceModel.getParameters(fState.getDate().getField()));
353
354 Vector3D dFdPXRef = new Vector3D(finiteDifferencesJacobian[0][0],
355 finiteDifferencesJacobian[1][0],
356 finiteDifferencesJacobian[2][0]);
357 Vector3D dFdPXRes = new Vector3D(gJacobian.getX().getPartialDerivative(0),
358 gJacobian.getY().getPartialDerivative(0),
359 gJacobian.getZ().getPartialDerivative(0));
360 Vector3D dFdPYRef = new Vector3D(finiteDifferencesJacobian[0][1],
361 finiteDifferencesJacobian[1][1],
362 finiteDifferencesJacobian[2][1]);
363 Vector3D dFdPYRes = new Vector3D(gJacobian.getX().getPartialDerivative(1),
364 gJacobian.getY().getPartialDerivative(1),
365 gJacobian.getZ().getPartialDerivative(1));
366 Vector3D dFdPZRef = new Vector3D(finiteDifferencesJacobian[0][2],
367 finiteDifferencesJacobian[1][2],
368 finiteDifferencesJacobian[2][2]);
369 Vector3D dFdPZRes = new Vector3D(gJacobian.getX().getPartialDerivative(2),
370 gJacobian.getY().getPartialDerivative(2),
371 gJacobian.getZ().getPartialDerivative(2));
372 Vector3D dFdVXRef = new Vector3D(finiteDifferencesJacobian[0][3],
373 finiteDifferencesJacobian[1][3],
374 finiteDifferencesJacobian[2][3]);
375 Vector3D dFdVXRes = new Vector3D(gJacobian.getX().getPartialDerivative(3),
376 gJacobian.getY().getPartialDerivative(3),
377 gJacobian.getZ().getPartialDerivative(3));
378 Vector3D dFdVYRef = new Vector3D(finiteDifferencesJacobian[0][4],
379 finiteDifferencesJacobian[1][4],
380 finiteDifferencesJacobian[2][4]);
381 Vector3D dFdVYRes = new Vector3D(gJacobian.getX().getPartialDerivative(4),
382 gJacobian.getY().getPartialDerivative(4),
383 gJacobian.getZ().getPartialDerivative(4));
384 Vector3D dFdVZRef = new Vector3D(finiteDifferencesJacobian[0][5],
385 finiteDifferencesJacobian[1][5],
386 finiteDifferencesJacobian[2][5]);
387 Vector3D dFdVZRes = new Vector3D(gJacobian.getX().getPartialDerivative(5),
388 gJacobian.getY().getPartialDerivative(5),
389 gJacobian.getZ().getPartialDerivative(5));
390 if (print) {
391 System.out.println("dF/dPX ref: " + dFdPXRef.getX() + " " + dFdPXRef.getY() + " " + dFdPXRef.getZ());
392 System.out.println("dF/dPX res: " + dFdPXRes.getX() + " " + dFdPXRes.getY() + " " + dFdPXRes.getZ());
393 System.out.println("dF/dPY ref: " + dFdPYRef.getX() + " " + dFdPYRef.getY() + " " + dFdPYRef.getZ());
394 System.out.println("dF/dPY res: " + dFdPYRes.getX() + " " + dFdPYRes.getY() + " " + dFdPYRes.getZ());
395 System.out.println("dF/dPZ ref: " + dFdPZRef.getX() + " " + dFdPZRef.getY() + " " + dFdPZRef.getZ());
396 System.out.println("dF/dPZ res: " + dFdPZRes.getX() + " " + dFdPZRes.getY() + " " + dFdPZRes.getZ());
397 System.out.println("dF/dPX ref norm: " + dFdPXRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPXRef, dFdPXRes) + ", rel error: " + (Vector3D.distance(dFdPXRef, dFdPXRes) / dFdPXRef.getNorm()));
398 System.out.println("dF/dPY ref norm: " + dFdPYRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPYRef, dFdPYRes) + ", rel error: " + (Vector3D.distance(dFdPYRef, dFdPYRes) / dFdPYRef.getNorm()));
399 System.out.println("dF/dPZ ref norm: " + dFdPZRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPZRef, dFdPZRes) + ", rel error: " + (Vector3D.distance(dFdPZRef, dFdPZRes) / dFdPZRef.getNorm()));
400 System.out.println("dF/dVX ref norm: " + dFdVXRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVXRef, dFdVXRes) + ", rel error: " + (Vector3D.distance(dFdVXRef, dFdVXRes) / dFdVXRef.getNorm()));
401 System.out.println("dF/dVY ref norm: " + dFdVYRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVYRef, dFdVYRes) + ", rel error: " + (Vector3D.distance(dFdVYRef, dFdVYRes) / dFdVYRef.getNorm()));
402 System.out.println("dF/dVZ ref norm: " + dFdVZRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVZRef, dFdVZRes) + ", rel error: " + (Vector3D.distance(dFdVZRef, dFdVZRes) / dFdVZRef.getNorm()));
403 }
404 checkdFdP(dFdPXRef, dFdPXRes, checkTolerance);
405 checkdFdP(dFdPYRef, dFdPYRes, checkTolerance);
406 checkdFdP(dFdPZRef, dFdPZRes, checkTolerance);
407 checkdFdP(dFdVXRef, dFdVXRes, checkTolerance);
408 checkdFdP(dFdVYRef, dFdVYRes, checkTolerance);
409 checkdFdP(dFdVZRef, dFdVZRes, checkTolerance);
410 }
411
412 private void checkdFdP(final Vector3D reference, final Vector3D result, final double checkTolerance) {
413 if (reference.getNorm() == 0) {
414
415
416 Assert.assertEquals(0, result.getNorm(), Precision.SAFE_MIN);
417 } else {
418 Assert.assertEquals(0, Vector3D.distance(reference, result), checkTolerance * reference.getNorm());
419 }
420 }
421
422 protected void checkStateJacobian(NumericalPropagator propagator, SpacecraftState state0,
423 AbsoluteDate targetDate, double hFactor,
424 double[] integratorAbsoluteTolerances, double checkTolerance) {
425
426 propagator.setInitialState(state0);
427 double[][] reference = new double[][] {
428 jacobianColumn(propagator, state0, targetDate, 0, hFactor * integratorAbsoluteTolerances[0]),
429 jacobianColumn(propagator, state0, targetDate, 1, hFactor * integratorAbsoluteTolerances[1]),
430 jacobianColumn(propagator, state0, targetDate, 2, hFactor * integratorAbsoluteTolerances[2]),
431 jacobianColumn(propagator, state0, targetDate, 3, hFactor * integratorAbsoluteTolerances[3]),
432 jacobianColumn(propagator, state0, targetDate, 4, hFactor * integratorAbsoluteTolerances[4]),
433 jacobianColumn(propagator, state0, targetDate, 5, hFactor * integratorAbsoluteTolerances[5])
434 };
435 for (int j = 0; j < 6; ++j) {
436 for (int k = j + 1; k < 6; ++k) {
437 double tmp = reference[j][k];
438 reference[j][k] = reference[k][j];
439 reference[k][j] = tmp;
440 }
441 }
442
443 propagator.setInitialState(state0);
444 MatricesHarvester harvester = propagator.setupMatricesComputation("stm", null, null);
445 final AtomicReference<RealMatrix> dYdY0 = new AtomicReference<>();
446 propagator.setStepHandler(new OrekitStepHandler() {
447 public void handleStep(OrekitStepInterpolator interpolator) {
448 }
449 public void finish(SpacecraftState finalState) {
450
451 dYdY0.set(harvester.getStateTransitionMatrix(finalState));
452 }
453
454 });
455 propagator.propagate(targetDate);
456
457 for (int j = 0; j < 6; ++j) {
458 for (int k = 0; k < 6; ++k) {
459 double scale = integratorAbsoluteTolerances[j] / integratorAbsoluteTolerances[k];
460 Assert.assertEquals(reference[j][k], dYdY0.get().getEntry(j, k), checkTolerance * scale);
461 }
462 }
463
464 }
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482 protected void checkRealFieldPropagation(final FieldKeplerianOrbit<DerivativeStructure> initialOrbit,
483 final PositionAngle positionAngle,
484 final double duration,
485 final NumericalPropagator propagator,
486 final FieldNumericalPropagator<DerivativeStructure> fieldpropagator,
487 final double positionRelativeTolerancePropag,
488 final double positionRelativeToleranceTaylor,
489 final double velocityRelativeToleranceTaylor,
490 final double accelerationRelativeToleranceTaylor,
491 final int nbTests,
492 final boolean print) {
493
494
495
496
497
498
499 FieldAbsoluteDate<DerivativeStructure> target = initialOrbit.getDate().shiftedBy(duration);
500 FieldSpacecraftState<DerivativeStructure> finalState_DS = fieldpropagator.propagate(target);
501 SpacecraftState finalState_R = propagator.propagate(target.toAbsoluteDate());
502 FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
503 PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
504
505
506 final Vector3D finPosition_DS = finPVC_DS.toPVCoordinates().getPosition();
507 final Vector3D finPosition_R = finPVC_R.getPosition();
508 final double finPositionDelta = Vector3D.distance(finPosition_DS, finPosition_R);
509 final double propagPosTol = finPosition_R.getNorm() * positionRelativeTolerancePropag;
510 if (print) {
511 System.out.println("1 - Check ΔP after propagation");
512 System.out.println("\t Pf_DS = " + finPosition_DS.getX() + " " + finPosition_DS.getY() + " " + finPosition_DS.getZ());
513 System.out.println("\t Pf_R = " + finPosition_R.getX() + " " + finPosition_R.getY() + " " + finPosition_R.getZ());
514 System.out.println("\t ΔPf = " + finPositionDelta);
515 System.out.println("\t εP = " + propagPosTol);
516 System.out.println("\tΔPf / ||Pf_R|| = " + finPositionDelta / finPosition_R.getNorm());
517 }
518 Assert.assertEquals(0., finPositionDelta, propagPosTol);
519
520
521
522
523
524
525
526 long number = 23091991;
527 RandomGenerator RG = new Well19937a(number);
528 GaussianRandomGenerator NGG = new GaussianRandomGenerator(RG);
529 UncorrelatedRandomVectorGenerator URVG = new UncorrelatedRandomVectorGenerator(new double[] {0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 },
530 new double[] {1e3, 0.005, 0.005, 0.01, 0.01, 0.01},
531 NGG);
532
533 double a_R = initialOrbit.getA().getReal();
534 double e_R = initialOrbit.getE().getReal();
535 double i_R = initialOrbit.getI().getReal();
536 double R_R = initialOrbit.getPerigeeArgument().getReal();
537 double O_R = initialOrbit.getRightAscensionOfAscendingNode().getReal();
538 double n_R = initialOrbit.getAnomaly(positionAngle).getReal();
539
540
541 double maxP = 0.;
542 double maxV = 0.;
543 double maxA = 0.;
544
545
546 for (int ii = 0; ii < nbTests; ii++){
547
548
549 double[] rand_next = URVG.nextVector();
550 double a_shift = a_R + rand_next[0];
551 double e_shift = e_R + rand_next[1];
552 double i_shift = i_R + rand_next[2];
553 double R_shift = R_R + rand_next[3];
554 double O_shift = O_R + rand_next[4];
555 double n_shift = n_R + rand_next[5];
556
557 KeplerianOrbit shiftedOrb = new KeplerianOrbit(a_shift, e_shift, i_shift, R_shift, O_shift, n_shift,
558 PositionAngle.MEAN,
559 initialOrbit.getFrame(),
560 initialOrbit.getDate().toAbsoluteDate(),
561 Constants.EIGEN5C_EARTH_MU
562 );
563
564 SpacecraftState shift_iSR = new SpacecraftState(shiftedOrb);
565
566
567 propagator.setInitialState(shift_iSR);
568 SpacecraftState finalState_shift = propagator.propagate(target.toAbsoluteDate());
569 PVCoordinates finPVC_shift = finalState_shift.getPVCoordinates();
570
571
572
573
574
575
576 FieldVector3D<DerivativeStructure> pos_DS = finPVC_DS.getPosition();
577 double x_DS = pos_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
578 double y_DS = pos_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
579 double z_DS = pos_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
580
581
582 double x = finPVC_shift.getPosition().getX();
583 double y = finPVC_shift.getPosition().getY();
584 double z = finPVC_shift.getPosition().getZ();
585
586
587 maxP = FastMath.max(maxP, FastMath.abs((x_DS - x) / (x - pos_DS.getX().getReal())));
588 maxP = FastMath.max(maxP, FastMath.abs((y_DS - y) / (y - pos_DS.getY().getReal())));
589 maxP = FastMath.max(maxP, FastMath.abs((z_DS - z) / (z - pos_DS.getZ().getReal())));
590
591
592
593
594
595
596 FieldVector3D<DerivativeStructure> vel_DS = finPVC_DS.getVelocity();
597 double vx_DS = vel_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
598 double vy_DS = vel_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
599 double vz_DS = vel_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
600
601
602 double vx = finPVC_shift.getVelocity().getX();
603 double vy = finPVC_shift.getVelocity().getY();
604 double vz = finPVC_shift.getVelocity().getZ();
605
606
607 maxV = FastMath.max(maxV, FastMath.abs((vx_DS - vx) / vx));
608 maxV = FastMath.max(maxV, FastMath.abs((vy_DS - vy) / vy));
609 maxV = FastMath.max(maxV, FastMath.abs((vz_DS - vz) / vz));
610
611
612
613
614
615
616 FieldVector3D<DerivativeStructure> acc_DS = finPVC_DS.getAcceleration();
617 double ax_DS = acc_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
618 double ay_DS = acc_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
619 double az_DS = acc_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
620
621
622 double ax = finPVC_shift.getAcceleration().getX();
623 double ay = finPVC_shift.getAcceleration().getY();
624 double az = finPVC_shift.getAcceleration().getZ();
625
626
627 maxA = FastMath.max(maxA, FastMath.abs((ax_DS - ax) / ax));
628 maxA = FastMath.max(maxA, FastMath.abs((ay_DS - ay) / ay));
629 maxA = FastMath.max(maxA, FastMath.abs((az_DS - az) / az));
630 }
631
632
633 if (print) {
634 System.out.println("\n2 - Check P, V, A with Taylor expansion");
635 System.out.println("\tMax ΔP = " + maxP);
636 System.out.println("\tMax ΔV = " + maxV);
637 System.out.println("\tMax ΔA = " + maxA);
638 }
639 Assert.assertEquals(0, maxP, positionRelativeToleranceTaylor);
640 Assert.assertEquals(0, maxV, velocityRelativeToleranceTaylor);
641 Assert.assertEquals(0, maxA, accelerationRelativeToleranceTaylor);
642 }
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659 protected void checkRealFieldPropagationGradient(final FieldKeplerianOrbit<Gradient> initialOrbit,
660 final PositionAngle positionAngle,
661 final double duration,
662 final NumericalPropagator propagator,
663 final FieldNumericalPropagator<Gradient> fieldpropagator,
664 final double positionRelativeTolerancePropag,
665 final double positionRelativeToleranceTaylor,
666 final double velocityRelativeToleranceTaylor,
667 final double accelerationRelativeToleranceTaylor,
668 final int nbTests,
669 final boolean print) {
670
671
672
673
674
675
676 FieldAbsoluteDate<Gradient> target = initialOrbit.getDate().shiftedBy(duration);
677 FieldSpacecraftState<Gradient> finalState_G = fieldpropagator.propagate(target);
678 SpacecraftState finalState_R = propagator.propagate(target.toAbsoluteDate());
679 FieldPVCoordinates<Gradient> finPVC_G = finalState_G.getPVCoordinates();
680 PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
681
682
683 final Vector3D finPosition_G = finPVC_G.toPVCoordinates().getPosition();
684 final Vector3D finPosition_R = finPVC_R.getPosition();
685 final double finPositionDelta = Vector3D.distance(finPosition_G, finPosition_R);
686 final double propagPosTol = finPosition_R.getNorm() * positionRelativeTolerancePropag;
687 if (print) {
688 System.out.println("1 - Check ΔP after propagation");
689 System.out.println("\t Pf_DS = " + finPosition_G.getX() + " " + finPosition_G.getY() + " " + finPosition_G.getZ());
690 System.out.println("\t Pf_R = " + finPosition_R.getX() + " " + finPosition_R.getY() + " " + finPosition_R.getZ());
691 System.out.println("\t ΔPf = " + finPositionDelta);
692 System.out.println("\t εP = " + propagPosTol);
693 System.out.println("\tΔPf / ||Pf_R|| = " + finPositionDelta / finPosition_R.getNorm());
694 }
695 Assert.assertEquals(0., finPositionDelta, propagPosTol);
696
697
698
699
700
701
702
703 long number = 23091991;
704 RandomGenerator RG = new Well19937a(number);
705 GaussianRandomGenerator NGG = new GaussianRandomGenerator(RG);
706 UncorrelatedRandomVectorGenerator URVG = new UncorrelatedRandomVectorGenerator(new double[] {0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 },
707 new double[] {1e3, 0.005, 0.005, 0.01, 0.01, 0.01},
708 NGG);
709
710 double a_R = initialOrbit.getA().getReal();
711 double e_R = initialOrbit.getE().getReal();
712 double i_R = initialOrbit.getI().getReal();
713 double R_R = initialOrbit.getPerigeeArgument().getReal();
714 double O_R = initialOrbit.getRightAscensionOfAscendingNode().getReal();
715 double n_R = initialOrbit.getAnomaly(positionAngle).getReal();
716
717
718 double maxP = 0.;
719 double maxV = 0.;
720 double maxA = 0.;
721
722
723 for (int ii = 0; ii < nbTests; ii++){
724
725
726 double[] rand_next = URVG.nextVector();
727 double a_shift = a_R + rand_next[0];
728 double e_shift = e_R + rand_next[1];
729 double i_shift = i_R + rand_next[2];
730 double R_shift = R_R + rand_next[3];
731 double O_shift = O_R + rand_next[4];
732 double n_shift = n_R + rand_next[5];
733
734 KeplerianOrbit shiftedOrb = new KeplerianOrbit(a_shift, e_shift, i_shift, R_shift, O_shift, n_shift,
735 PositionAngle.MEAN,
736 initialOrbit.getFrame(),
737 initialOrbit.getDate().toAbsoluteDate(),
738 Constants.EIGEN5C_EARTH_MU
739 );
740
741 SpacecraftState shift_iSR = new SpacecraftState(shiftedOrb);
742
743
744 propagator.setInitialState(shift_iSR);
745 SpacecraftState finalState_shift = propagator.propagate(target.toAbsoluteDate());
746 PVCoordinates finPVC_shift = finalState_shift.getPVCoordinates();
747
748
749
750
751
752
753 FieldVector3D<Gradient> pos_DS = finPVC_G.getPosition();
754 double x_DS = pos_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
755 double y_DS = pos_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
756 double z_DS = pos_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
757
758
759 double x = finPVC_shift.getPosition().getX();
760 double y = finPVC_shift.getPosition().getY();
761 double z = finPVC_shift.getPosition().getZ();
762
763
764 maxP = FastMath.max(maxP, FastMath.abs((x_DS - x) / (x - pos_DS.getX().getReal())));
765 maxP = FastMath.max(maxP, FastMath.abs((y_DS - y) / (y - pos_DS.getY().getReal())));
766 maxP = FastMath.max(maxP, FastMath.abs((z_DS - z) / (z - pos_DS.getZ().getReal())));
767
768
769
770
771
772
773 FieldVector3D<Gradient> vel_G = finPVC_G.getVelocity();
774 double vx_G = vel_G.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
775 double vy_G = vel_G.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
776 double vz_G = vel_G.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
777
778
779 double vx = finPVC_shift.getVelocity().getX();
780 double vy = finPVC_shift.getVelocity().getY();
781 double vz = finPVC_shift.getVelocity().getZ();
782
783
784 maxV = FastMath.max(maxV, FastMath.abs((vx_G - vx) / vx));
785 maxV = FastMath.max(maxV, FastMath.abs((vy_G - vy) / vy));
786 maxV = FastMath.max(maxV, FastMath.abs((vz_G - vz) / vz));
787
788
789
790
791
792
793 FieldVector3D<Gradient> acc_G = finPVC_G.getAcceleration();
794 double ax_G = acc_G.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
795 double ay_G = acc_G.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
796 double az_G = acc_G.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
797
798
799 double ax = finPVC_shift.getAcceleration().getX();
800 double ay = finPVC_shift.getAcceleration().getY();
801 double az = finPVC_shift.getAcceleration().getZ();
802
803
804 maxA = FastMath.max(maxA, FastMath.abs((ax_G - ax) / ax));
805 maxA = FastMath.max(maxA, FastMath.abs((ay_G - ay) / ay));
806 maxA = FastMath.max(maxA, FastMath.abs((az_G - az) / az));
807 }
808
809
810 if (print) {
811 System.out.println("\n2 - Check P, V, A with Taylor expansion");
812 System.out.println("\tMax ΔP = " + maxP);
813 System.out.println("\tMax ΔV = " + maxV);
814 System.out.println("\tMax ΔA = " + maxA);
815 }
816 Assert.assertEquals(0, maxP, positionRelativeToleranceTaylor);
817 Assert.assertEquals(0, maxV, velocityRelativeToleranceTaylor);
818 Assert.assertEquals(0, maxA, accelerationRelativeToleranceTaylor);
819 }
820
821 private double[] jacobianColumn(final NumericalPropagator propagator, final SpacecraftState state0,
822 final AbsoluteDate targetDate, final int index,
823 final double h)
824 {
825 return differential4(integrateShiftedState(propagator, state0, targetDate, index, -2 * h),
826 integrateShiftedState(propagator, state0, targetDate, index, -1 * h),
827 integrateShiftedState(propagator, state0, targetDate, index, +1 * h),
828 integrateShiftedState(propagator, state0, targetDate, index, +2 * h),
829 h);
830 }
831
832 private double[] integrateShiftedState(final NumericalPropagator propagator,
833 final SpacecraftState state0,
834 final AbsoluteDate targetDate,
835 final int index, final double h) {
836 OrbitType orbitType = propagator.getOrbitType();
837 PositionAngle angleType = propagator.getPositionAngleType();
838 double[] a = new double[6];
839 double[] aDot = new double[6];
840 orbitType.mapOrbitToArray(state0.getOrbit(), angleType, a, aDot);
841 a[index] += h;
842 SpacecraftState shiftedState = new SpacecraftState(orbitType.mapArrayToOrbit(a, aDot, angleType, state0.getDate(),
843 state0.getMu(), state0.getFrame()),
844 state0.getAttitude(),
845 state0.getMass());
846 propagator.setInitialState(shiftedState);
847 SpacecraftState integratedState = propagator.propagate(targetDate);
848 orbitType.mapOrbitToArray(integratedState.getOrbit(), angleType, a, null);
849 return a;
850 }
851
852 protected double[] differential8(final double[] fM4h, final double[] fM3h, final double[] fM2h, final double[] fM1h,
853 final double[] fP1h, final double[] fP2h, final double[] fP3h, final double[] fP4h,
854 final double h) {
855
856 double[] a = new double[fM4h.length];
857 for (int i = 0; i < a.length; ++i) {
858 a[i] = differential8(fM4h[i], fM3h[i], fM2h[i], fM1h[i], fP1h[i], fP2h[i], fP3h[i], fP4h[i], h);
859 }
860 return a;
861 }
862
863 protected double differential8(final double fM4h, final double fM3h, final double fM2h, final double fM1h,
864 final double fP1h, final double fP2h, final double fP3h, final double fP4h,
865 final double h) {
866
867
868
869 return (-3 * (fP4h - fM4h) + 32 * (fP3h - fM3h) - 168 * (fP2h - fM2h) + 672 * (fP1h - fM1h)) / (840 * h);
870
871 }
872
873 protected double[] differential4(final double[] fM2h, final double[] fM1h,
874 final double[] fP1h, final double[] fP2h,
875 final double h) {
876
877 double[] a = new double[fM2h.length];
878 for (int i = 0; i < a.length; ++i) {
879 a[i] = differential4(fM2h[i], fM1h[i], fP1h[i], fP2h[i], h);
880 }
881 return a;
882 }
883
884 protected double differential4(final double fM2h, final double fM1h,
885 final double fP1h, final double fP2h,
886 final double h) {
887
888
889
890 return (-1 * (fP2h - fM2h) + 8 * (fP1h - fM1h)) / (12 * h);
891
892 }
893
894 }
895
896