1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.radiation;
18
19 import java.io.FileNotFoundException;
20 import java.text.ParseException;
21
22 import org.hipparchus.Field;
23 import org.hipparchus.analysis.differentiation.DSFactory;
24 import org.hipparchus.analysis.differentiation.DerivativeStructure;
25 import org.hipparchus.analysis.differentiation.Gradient;
26 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
27 import org.hipparchus.geometry.euclidean.threed.Vector3D;
28 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
29 import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaFieldIntegrator;
30 import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator;
31 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
32 import org.hipparchus.util.FastMath;
33 import org.junit.Assert;
34 import org.junit.Before;
35 import org.junit.Test;
36 import org.orekit.Utils;
37 import org.orekit.attitudes.Attitude;
38 import org.orekit.attitudes.AttitudeProvider;
39 import org.orekit.bodies.CelestialBodyFactory;
40 import org.orekit.bodies.OneAxisEllipsoid;
41 import org.orekit.errors.OrekitException;
42 import org.orekit.forces.AbstractForceModelTest;
43 import org.orekit.forces.ForceModel;
44 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
45 import org.orekit.forces.gravity.potential.GravityFieldFactory;
46 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
47 import org.orekit.frames.Frame;
48 import org.orekit.frames.FramesFactory;
49 import org.orekit.orbits.CartesianOrbit;
50 import org.orekit.orbits.EquinoctialOrbit;
51 import org.orekit.orbits.FieldCartesianOrbit;
52 import org.orekit.orbits.FieldKeplerianOrbit;
53 import org.orekit.orbits.FieldOrbit;
54 import org.orekit.orbits.KeplerianOrbit;
55 import org.orekit.orbits.Orbit;
56 import org.orekit.orbits.OrbitType;
57 import org.orekit.orbits.PositionAngle;
58 import org.orekit.propagation.FieldSpacecraftState;
59 import org.orekit.propagation.SpacecraftState;
60 import org.orekit.propagation.integration.AbstractIntegratedPropagator;
61 import org.orekit.propagation.numerical.FieldNumericalPropagator;
62 import org.orekit.propagation.numerical.NumericalPropagator;
63 import org.orekit.propagation.sampling.OrekitFixedStepHandler;
64 import org.orekit.time.AbsoluteDate;
65 import org.orekit.time.DateComponents;
66 import org.orekit.time.FieldAbsoluteDate;
67 import org.orekit.time.TimeComponents;
68 import org.orekit.time.TimeScalesFactory;
69 import org.orekit.utils.Constants;
70 import org.orekit.utils.ExtendedPVCoordinatesProvider;
71 import org.orekit.utils.FieldPVCoordinates;
72 import org.orekit.utils.IERSConventions;
73 import org.orekit.utils.PVCoordinates;
74
75 public class ECOM2Test extends AbstractForceModelTest {
76
77 private static final AttitudeProvider DEFAULT_LAW = Utils.defaultLaw();
78
79 @Before
80 public void setUp() throws OrekitException {
81 Utils.setDataRoot("potential/shm-format:regular-data");
82 }
83
84 private NumericalPropagator setUpPropagator(Orbit orbit, double dP,
85 OrbitType orbitType, PositionAngle angleType,
86 ForceModel... models)
87 {
88
89 final double minStep = 0.001;
90 final double maxStep = 1000;
91
92 double[][] tol = NumericalPropagator.tolerances(dP, orbit, orbitType);
93 NumericalPropagator propagator =
94 new NumericalPropagator(new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]));
95 propagator.setOrbitType(orbitType);
96 propagator.setPositionAngleType(angleType);
97 for (ForceModel model : models) {
98 propagator.addForceModel(model);
99 }
100 return propagator;
101 }
102
103 private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType, PositionAngle angleType,
104 double delta, int column) {
105
106 double[][] array = stateToArray(state, orbitType, angleType, true);
107 array[0][column] += delta;
108
109 return arrayToState(array, orbitType, angleType, state.getFrame(), state.getDate(),
110 state.getMu(), state.getAttitude());
111
112 }
113
114 private SpacecraftState arrayToState(double[][] array, OrbitType orbitType, PositionAngle angleType,
115 Frame frame, AbsoluteDate date, double mu,
116 Attitude attitude) {
117 Orbit orbit = orbitType.mapArrayToOrbit(array[0], array[1], angleType, date, mu, frame);
118 return (array.length > 6) ?
119 new SpacecraftState(orbit, attitude) :
120 new SpacecraftState(orbit, attitude, array[0][6]);
121 }
122
123 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType, PositionAngle angleType,
124 boolean withMass) {
125 double[][] array = new double[2][withMass ? 7 : 6];
126 orbitType.mapOrbitToArray(state.getOrbit(), angleType, array[0], array[1]);
127 if (withMass) {
128 array[0][6] = state.getMass();
129 }
130 return array;
131 }
132
133 private void fillJacobianModelColumn(double[][] jacobian, int column,
134 OrbitType orbitType, PositionAngle angleType, double h,
135 Vector3D sM4h, Vector3D sM3h,
136 Vector3D sM2h, Vector3D sM1h,
137 Vector3D sP1h, Vector3D sP2h,
138 Vector3D sP3h, Vector3D sP4h) {
139
140 jacobian[0][column] = ( -3 * (sP4h.getX() - sM4h.getX()) +
141 32 * (sP3h.getX() - sM3h.getX()) -
142 168 * (sP2h.getX() - sM2h.getX()) +
143 672 * (sP1h.getX() - sM1h.getX())) / (840 * h);
144 jacobian[1][column] = ( -3 * (sP4h.getY() - sM4h.getY()) +
145 32 * (sP3h.getY() - sM3h.getY()) -
146 168 * (sP2h.getY() - sM2h.getY()) +
147 672 * (sP1h.getY() - sM1h.getY())) / (840 * h);
148 jacobian[2][column] = ( -3 * (sP4h.getZ() - sM4h.getZ()) +
149 32 * (sP3h.getZ() - sM3h.getZ()) -
150 168 * (sP2h.getZ() - sM2h.getZ()) +
151 672 * (sP1h.getZ() - sM1h.getZ())) / (840 * h);
152 }
153
154 @Test
155 public void testJacobianModelMatrix() {
156 final DSFactory factory = new DSFactory(6, 1);
157 NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(2, 0);
158 ForceModel gravityField =
159 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
160
161
162 final DerivativeStructure x = factory.variable(0, -2747600.0);
163 final DerivativeStructure y = factory.variable(1, 22572100.0);
164 final DerivativeStructure z = factory.variable(2, 13522760.0);
165 final DerivativeStructure xDot = factory.variable(3, -2729.5);
166 final DerivativeStructure yDot = factory.variable(4, 1142.7);
167 final DerivativeStructure zDot = factory.variable(5, -2523.9);
168
169
170 final Field<DerivativeStructure> field = x.getField();
171 final DerivativeStructure one = field.getOne();
172 final FieldVector3D<DerivativeStructure> pos = new FieldVector3D<DerivativeStructure>(x, y, z);
173 final FieldVector3D<DerivativeStructure> vel = new FieldVector3D<DerivativeStructure>(xDot, yDot, zDot);
174 final FieldPVCoordinates<DerivativeStructure> dsPV = new FieldPVCoordinates<DerivativeStructure>(pos, vel);
175 final FieldAbsoluteDate<DerivativeStructure> dsDate = new FieldAbsoluteDate<>(field, new AbsoluteDate(2003, 2, 13, 2, 31, 30, TimeScalesFactory.getUTC()));
176 final DerivativeStructure mu = one.multiply(Constants.EGM96_EARTH_MU);
177 final FieldOrbit<DerivativeStructure> dsOrbit = new FieldCartesianOrbit<DerivativeStructure>(dsPV, FramesFactory.getEME2000(), dsDate, mu);
178 final FieldSpacecraftState<DerivativeStructure> dsState = new FieldSpacecraftState<DerivativeStructure>(dsOrbit);
179
180
181 final ECOM2 forceModel = new ECOM2(2, 2, 1e-7, CelestialBodyFactory.getSun(), Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
182
183
184 final FieldVector3D<DerivativeStructure> acc = forceModel.acceleration(dsState, forceModel.getParameters(field));
185 final double[] accX = acc.getX().getAllDerivatives();
186 final double[] accY = acc.getY().getAllDerivatives();
187 final double[] accZ = acc.getZ().getAllDerivatives();
188
189
190 final Orbit orbit = dsOrbit.toOrbit();
191 final SpacecraftState state = dsState.toSpacecraftState();
192 final double[][] refDeriv = new double[3][6];
193 final OrbitType orbitType = orbit.getType();
194 final PositionAngle angleType = PositionAngle.MEAN;
195 double dP = 0.001;
196 double[] steps = NumericalPropagator.tolerances(1000000 * dP, orbit, orbitType)[0];
197 AbstractIntegratedPropagator propagator = setUpPropagator(orbit, dP, orbitType, angleType, gravityField, forceModel);
198
199
200 for(int i = 0; i < 6; i++) {
201 propagator.resetInitialState(shiftState(state, orbitType, angleType, -4 * steps[i], i));
202 SpacecraftState sM4h = propagator.propagate(state.getDate());
203 Vector3D accM4 = forceModel.acceleration(sM4h, forceModel.getParameters());
204
205 propagator.resetInitialState(shiftState(state, orbitType, angleType, -3 * steps[i], i));
206 SpacecraftState sM3h = propagator.propagate(state.getDate());
207 Vector3D accM3 = forceModel.acceleration(sM3h, forceModel.getParameters());
208
209 propagator.resetInitialState(shiftState(state, orbitType, angleType, -2 * steps[i], i));
210 SpacecraftState sM2h = propagator.propagate(state.getDate());
211 Vector3D accM2 = forceModel.acceleration(sM2h, forceModel.getParameters());
212
213 propagator.resetInitialState(shiftState(state, orbitType, angleType, -1 * steps[i] , i));
214 SpacecraftState sM1h = propagator.propagate(state.getDate());
215 Vector3D accM1 = forceModel.acceleration(sM1h, forceModel.getParameters());
216
217 propagator.resetInitialState(shiftState(state, orbitType, angleType, 1 * steps[i], i));
218 SpacecraftState sP1h = propagator.propagate(state.getDate());
219 Vector3D accP1 = forceModel.acceleration(sP1h, forceModel.getParameters());
220
221 propagator.resetInitialState(shiftState(state, orbitType, angleType, 2 * steps[i], i));
222 SpacecraftState sP2h = propagator.propagate(state.getDate());
223 Vector3D accP2 = forceModel.acceleration(sP2h, forceModel.getParameters());
224
225 propagator.resetInitialState(shiftState(state, orbitType, angleType, 3 * steps[i], i));
226 SpacecraftState sP3h = propagator.propagate(state.getDate());
227 Vector3D accP3 = forceModel.acceleration(sP3h, forceModel.getParameters());
228
229 propagator.resetInitialState(shiftState(state, orbitType, angleType, 4 * steps[i], i));
230 SpacecraftState sP4h = propagator.propagate(state.getDate());
231 Vector3D accP4 = forceModel.acceleration(sP4h, forceModel.getParameters());
232 fillJacobianModelColumn(refDeriv, i, orbitType, angleType, steps[i],
233 accM4, accM3, accM2, accM1,
234 accP1, accP2, accP3, accP4);
235 }
236
237
238 for (int i = 0; i < 6; i++) {
239 final double errorX = (accX[i + 1] - refDeriv[0][i]) / refDeriv[0][i];
240 Assert.assertEquals(0, errorX, 1.0e-10);
241 final double errorY = (accY[i + 1] - refDeriv[1][i]) / refDeriv[1][i];
242 Assert.assertEquals(0, errorY, 1.5e-10);
243 final double errorZ = (accZ[i + 1] - refDeriv[2][i]) / refDeriv[2][i];
244 Assert.assertEquals(0, errorZ, 1.0e-10);
245 }
246
247 }
248
249 @Test
250 public void testRealField() {
251
252
253
254 DSFactory factory = new DSFactory(6, 4);
255 DerivativeStructure a_0 = factory.variable(0, 7e6);
256 DerivativeStructure e_0 = factory.variable(1, 0.01);
257 DerivativeStructure i_0 = factory.variable(2, 1.2);
258 DerivativeStructure R_0 = factory.variable(3, 0.7);
259 DerivativeStructure O_0 = factory.variable(4, 0.5);
260 DerivativeStructure n_0 = factory.variable(5, 0.1);
261
262 Field<DerivativeStructure> field = a_0.getField();
263 DerivativeStructure zero = field.getZero();
264
265
266 FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
267
268
269 Frame EME = FramesFactory.getEME2000();
270
271
272 FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
273 PositionAngle.MEAN,
274 EME,
275 J2000,
276 zero.add(Constants.EIGEN5C_EARTH_MU));
277
278
279 FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
280 SpacecraftState iSR = initialState.toSpacecraftState();
281
282
283 ClassicalRungeKuttaFieldIntegrator<DerivativeStructure> integrator =
284 new ClassicalRungeKuttaFieldIntegrator<>(field, zero.add(6));
285 ClassicalRungeKuttaIntegrator RIntegrator =
286 new ClassicalRungeKuttaIntegrator(6);
287 OrbitType type = OrbitType.EQUINOCTIAL;
288
289
290 FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
291 FNP.setOrbitType(type);
292 FNP.setInitialState(initialState);
293
294 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
295 NP.setOrbitType(type);
296 NP.setInitialState(iSR);
297
298
299 final ECOM2 forceModel = new ECOM2(0, 0, 1e-7, CelestialBodyFactory.getSun(), Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
300
301 FNP.addForceModel(forceModel);
302 NP.addForceModel(forceModel);
303
304
305 checkRealFieldPropagation(FKO, PositionAngle.MEAN, 300., NP, FNP,
306 1.0e-30, 1.3e-8, 6.7e-11, 1.4e-10,
307 1, false);
308 }
309
310 @Test
311 public void testRealFieldGradient() {
312
313
314
315 int freeParameters = 6;
316 Gradient a_0 = Gradient.variable(freeParameters, 0, 7e6);
317 Gradient e_0 = Gradient.variable(freeParameters, 1, 0.01);
318 Gradient i_0 = Gradient.variable(freeParameters, 2, 1.2);
319 Gradient R_0 = Gradient.variable(freeParameters, 3, 0.7);
320 Gradient O_0 = Gradient.variable(freeParameters, 4, 0.5);
321 Gradient n_0 = Gradient.variable(freeParameters, 5, 0.1);
322
323 Field<Gradient> field = a_0.getField();
324 Gradient zero = field.getZero();
325
326
327 FieldAbsoluteDate<Gradient> J2000 = new FieldAbsoluteDate<>(field);
328
329
330 Frame EME = FramesFactory.getEME2000();
331
332
333 FieldKeplerianOrbit<Gradient> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
334 PositionAngle.MEAN,
335 EME,
336 J2000,
337 zero.add(Constants.EIGEN5C_EARTH_MU));
338
339
340 FieldSpacecraftState<Gradient> initialState = new FieldSpacecraftState<>(FKO);
341 SpacecraftState iSR = initialState.toSpacecraftState();
342
343
344 ClassicalRungeKuttaFieldIntegrator<Gradient> integrator =
345 new ClassicalRungeKuttaFieldIntegrator<>(field, zero.add(6));
346 ClassicalRungeKuttaIntegrator RIntegrator =
347 new ClassicalRungeKuttaIntegrator(6);
348 OrbitType type = OrbitType.EQUINOCTIAL;
349
350
351 FieldNumericalPropagator<Gradient> FNP = new FieldNumericalPropagator<>(field, integrator);
352 FNP.setOrbitType(type);
353 FNP.setInitialState(initialState);
354
355 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
356 NP.setOrbitType(type);
357 NP.setInitialState(iSR);
358
359
360 final ECOM2 forceModel = new ECOM2(0, 0, 1e-7, CelestialBodyFactory.getSun(), Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
361
362 FNP.addForceModel(forceModel);
363 NP.addForceModel(forceModel);
364
365
366 checkRealFieldPropagationGradient(FKO, PositionAngle.MEAN, 300., NP, FNP,
367 1.0e-30, 1.3e-2, 9.6e-5, 1.4e-4,
368 1, false);
369 }
370
371 @Test
372 public void testParameterDerivative() {
373
374 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
375 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
376 final SpacecraftState state =
377 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
378 FramesFactory.getGCRF(),
379 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
380 Constants.EIGEN5C_EARTH_MU));
381
382
383 final ECOM2 forceModel = new ECOM2(0, 0, 1e-7, CelestialBodyFactory.getSun(), Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
384
385 Assert.assertFalse(forceModel.dependsOnPositionOnly());
386
387 checkParameterDerivative(state, forceModel, ECOM2.ECOM_COEFFICIENT + " B0", 1.0e-4, 3.0e-12);
388 checkParameterDerivative(state, forceModel, ECOM2.ECOM_COEFFICIENT + " D0", 1.0e-4, 3.0e-12);
389 checkParameterDerivative(state, forceModel, ECOM2.ECOM_COEFFICIENT + " Y0", 1.0e-4, 3.0e-12);
390 }
391
392 @Test
393 public void testParameterDerivativeGradient() {
394
395 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
396 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
397 final SpacecraftState state =
398 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
399 FramesFactory.getGCRF(),
400 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
401 Constants.EIGEN5C_EARTH_MU));
402
403
404 final ECOM2 forceModel = new ECOM2(0, 0, 1e-7, CelestialBodyFactory.getSun(), Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
405
406 Assert.assertFalse(forceModel.dependsOnPositionOnly());
407
408 checkParameterDerivativeGradient(state, forceModel, ECOM2.ECOM_COEFFICIENT + " B0", 1.0e-4, 3.0e-12);
409 checkParameterDerivativeGradient(state, forceModel, ECOM2.ECOM_COEFFICIENT + " D0", 1.0e-4, 3.0e-12);
410 checkParameterDerivativeGradient(state, forceModel, ECOM2.ECOM_COEFFICIENT + " Y0", 1.0e-4, 3.0e-12);
411 }
412
413 @Test
414 public void testJacobianVsFiniteDifferences() {
415
416
417 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
418 new TimeComponents(13, 59, 27.816),
419 TimeScalesFactory.getUTC());
420 double i = FastMath.toRadians(98.7);
421 double omega = FastMath.toRadians(93.0);
422 double OMEGA = FastMath.toRadians(15.0 * 22.5);
423 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
424 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
425 Constants.EIGEN5C_EARTH_MU);
426
427 final ECOM2 forceModel = new ECOM2(2, 2, 1e-7, CelestialBodyFactory.getSun(), Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
428
429 SpacecraftState state = new SpacecraftState(orbit,
430 DEFAULT_LAW.getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
431 checkStateJacobianVsFiniteDifferences(state, forceModel, DEFAULT_LAW, 1.0, 5.0e-6, false);
432
433 }
434
435 @Test
436 public void testJacobianVsFiniteDifferencesGradient() {
437
438
439 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
440 new TimeComponents(13, 59, 27.816),
441 TimeScalesFactory.getUTC());
442 double i = FastMath.toRadians(98.7);
443 double omega = FastMath.toRadians(93.0);
444 double OMEGA = FastMath.toRadians(15.0 * 22.5);
445 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
446 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
447 Constants.EIGEN5C_EARTH_MU);
448
449 final ECOM2 forceModel = new ECOM2(2, 2, 1e-7, CelestialBodyFactory.getSun(), Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
450
451 SpacecraftState state = new SpacecraftState(orbit,
452 DEFAULT_LAW.getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
453 checkStateJacobianVsFiniteDifferencesGradient(state, forceModel, DEFAULT_LAW, 1.0, 5.0e-6, false);
454
455 }
456
457 @Test
458 public void testRoughOrbitalModifs() throws ParseException, OrekitException, FileNotFoundException {
459
460
461 AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 7, 1),
462 new TimeComponents(13, 59, 27.816),
463 TimeScalesFactory.getUTC());
464 Orbit orbit = new EquinoctialOrbit(42164000, 10e-3, 10e-3,
465 FastMath.tan(0.001745329)*FastMath.cos(2*FastMath.PI/3),
466 FastMath.tan(0.001745329)*FastMath.sin(2*FastMath.PI/3),
467 0.1, PositionAngle.TRUE, FramesFactory.getEME2000(), date, Constants.WGS84_EARTH_MU);
468 final double period = orbit.getKeplerianPeriod();
469 Assert.assertEquals(86164, period, 1);
470 ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
471
472
473 OneAxisEllipsoid earth =
474 new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765,
475 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
476 final ECOM2 SRP = new ECOM2(2, 2, 1e-7, sun, earth.getEquatorialRadius());
477
478
479 double[] absTolerance = {
480 0.1, 1.0e-9, 1.0e-9, 1.0e-5, 1.0e-5, 1.0e-5, 0.001
481 };
482 double[] relTolerance = {
483 1.0e-4, 1.0e-4, 1.0e-4, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-7
484 };
485 AdaptiveStepsizeIntegrator integrator =
486 new DormandPrince853Integrator(900.0, 60000, absTolerance, relTolerance);
487 integrator.setInitialStepSize(3600);
488 final NumericalPropagator calc = new NumericalPropagator(integrator);
489 calc.addForceModel(SRP);
490
491
492 calc.setStepHandler(FastMath.floor(period), new SolarStepHandler());
493 AbsoluteDate finalDate = date.shiftedBy(10 * period);
494 calc.setInitialState(new SpacecraftState(orbit, 1500.0));
495 calc.propagate(finalDate);
496 Assert.assertTrue(calc.getCalls() < 7100);
497 }
498
499 @Test
500 public void testRealAndFieldComparison() {
501
502
503 final int freeParameters = 6;
504 final Gradient sma = Gradient.variable(freeParameters, 0, 26559614.1);
505 final Gradient ecc = Gradient.variable(freeParameters, 1, 0.00522136);
506 final Gradient inc = Gradient.variable(freeParameters, 2, 0.963785748);
507 final Gradient aop = Gradient.variable(freeParameters, 3, 0.451712027);
508 final Gradient raan = Gradient.variable(freeParameters, 4, -1.159458779);
509 final Gradient lm = Gradient.variable(freeParameters, 4, -2.105941778);
510
511
512 final Field<Gradient> field = sma.getField();
513
514
515 final FieldAbsoluteDate<Gradient> epoch = FieldAbsoluteDate.getJ2000Epoch(field);
516
517
518 FieldKeplerianOrbit<Gradient> orbit = new FieldKeplerianOrbit<>(sma, ecc, inc, aop, raan, lm,
519 PositionAngle.MEAN,
520 FramesFactory.getEME2000(),
521 epoch,
522 field.getZero().add(Constants.EIGEN5C_EARTH_MU));
523
524
525 final ECOM2 forceModel = new ECOM2(2, 2, 1e-7, CelestialBodyFactory.getSun(), Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
526
527
528 final FieldVector3D<Gradient> accField = forceModel.acceleration(new FieldSpacecraftState<>(orbit), forceModel.getParameters(field));
529
530
531 final Vector3D accReal = forceModel.acceleration(new SpacecraftState(orbit.toOrbit()), forceModel.getParameters());
532
533
534 Assert.assertEquals(0.0, accReal.distance(accField.toVector3D()), 1.0e-20);
535
536 }
537
538 private static class SolarStepHandler implements OrekitFixedStepHandler {
539
540 public void handleStep(SpacecraftState currentState) {
541 final double dex = currentState.getEquinoctialEx() - 0.01071166;
542 final double dey = currentState.getEquinoctialEy() - 0.00654848;
543 final double alpha = FastMath.toDegrees(FastMath.atan2(dey, dex));
544 Assert.assertTrue(alpha > 100.0);
545 Assert.assertTrue(alpha < 112.0);
546 checkRadius(FastMath.sqrt(dex * dex + dey * dey), 0.003469, 0.003525);
547 }
548
549 }
550
551 public static void checkRadius(double radius , double min , double max) {
552 Assert.assertTrue(radius >= min);
553 Assert.assertTrue(radius <= max);
554 }
555
556
557 }