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