1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.drag;
18
19 import java.util.List;
20
21 import org.hipparchus.Field;
22 import org.hipparchus.analysis.differentiation.DSFactory;
23 import org.hipparchus.analysis.differentiation.DerivativeStructure;
24 import org.hipparchus.analysis.differentiation.Gradient;
25 import org.hipparchus.analysis.differentiation.GradientField;
26 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
27 import org.hipparchus.geometry.euclidean.threed.Vector3D;
28 import org.hipparchus.linear.RealMatrix;
29 import org.hipparchus.ode.AbstractIntegrator;
30 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
31 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
32 import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaFieldIntegrator;
33 import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator;
34 import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
35 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
36 import org.hipparchus.util.FastMath;
37 import org.junit.jupiter.api.Assertions;
38 import org.junit.jupiter.api.BeforeEach;
39 import org.junit.jupiter.api.Test;
40 import org.mockito.Mockito;
41 import org.orekit.Utils;
42 import org.orekit.attitudes.LofOffset;
43 import org.orekit.bodies.BodyShape;
44 import org.orekit.bodies.CelestialBodyFactory;
45 import org.orekit.bodies.OneAxisEllipsoid;
46 import org.orekit.errors.OrekitException;
47 import org.orekit.errors.OrekitMessages;
48 import org.orekit.forces.AbstractLegacyForceModelTest;
49 import org.orekit.forces.BoxAndSolarArraySpacecraft;
50 import org.orekit.forces.ForceModel;
51 import org.orekit.frames.Frame;
52 import org.orekit.frames.FramesFactory;
53 import org.orekit.frames.LOFType;
54 import org.orekit.frames.Transform;
55 import org.orekit.models.earth.atmosphere.Atmosphere;
56 import org.orekit.models.earth.atmosphere.HarrisPriester;
57 import org.orekit.models.earth.atmosphere.SimpleExponentialAtmosphere;
58 import org.orekit.orbits.CartesianOrbit;
59 import org.orekit.orbits.FieldKeplerianOrbit;
60 import org.orekit.orbits.KeplerianOrbit;
61 import org.orekit.orbits.Orbit;
62 import org.orekit.orbits.OrbitType;
63 import org.orekit.orbits.PositionAngleType;
64 import org.orekit.propagation.FieldSpacecraftState;
65 import org.orekit.propagation.MatricesHarvester;
66 import org.orekit.propagation.SpacecraftState;
67 import org.orekit.propagation.ToleranceProvider;
68 import org.orekit.propagation.numerical.FieldNumericalPropagator;
69 import org.orekit.propagation.numerical.NumericalPropagator;
70 import org.orekit.time.AbsoluteDate;
71 import org.orekit.time.DateComponents;
72 import org.orekit.time.FieldAbsoluteDate;
73 import org.orekit.time.TimeComponents;
74 import org.orekit.time.TimeScalesFactory;
75 import org.orekit.utils.Constants;
76 import org.orekit.utils.FieldPVCoordinates;
77 import org.orekit.utils.IERSConventions;
78 import org.orekit.utils.PVCoordinates;
79 import org.orekit.utils.ParameterDriver;
80 import org.orekit.utils.TimeStampedPVCoordinates;
81 import org.orekit.utils.TimeSpanMap.Span;
82
83 class DragForceTest extends AbstractLegacyForceModelTest {
84
85 @Override
86 protected FieldVector3D<DerivativeStructure> accelerationDerivatives(final ForceModel forceModel,
87 final FieldSpacecraftState<DerivativeStructure> state) {
88 try {
89
90 final AbsoluteDate date = state.getDate().toAbsoluteDate();
91 final FieldVector3D<DerivativeStructure> position = state.getPVCoordinates().getPosition();
92 final FieldVector3D<DerivativeStructure> velocity = state.getPVCoordinates().getVelocity();
93 java.lang.reflect.Field atmosphereField = AbstractDragForceModel.class.getDeclaredField("atmosphere");
94 atmosphereField.setAccessible(true);
95 Atmosphere atmosphere = (Atmosphere) atmosphereField.get(forceModel);
96 java.lang.reflect.Field spacecraftField = DragForce.class.getDeclaredField("spacecraft");
97 spacecraftField.setAccessible(true);
98 DragSensitive spacecraft = (DragSensitive) spacecraftField.get(forceModel);
99
100
101 final DSFactory factory = state.getMass().getFactory();
102
103
104 final Frame atmFrame = atmosphere.getFrame();
105 final Transform toBody = state.getFrame().getTransformTo(atmFrame, date);
106 final FieldVector3D<DerivativeStructure> posBodyDS = toBody.transformPosition(position);
107 final Vector3D posBody = posBodyDS.toVector3D();
108 final Vector3D vAtmBody = atmosphere.getVelocity(date, posBody, atmFrame);
109
110
111
112
113
114
115 if (factory.getCompiler().getOrder() > 1) {
116 throw new OrekitException(OrekitMessages.OUT_OF_RANGE_DERIVATION_ORDER, factory.getCompiler().getOrder());
117 }
118 final double delta = 1.0;
119 final double x = posBody.getX();
120 final double y = posBody.getY();
121 final double z = posBody.getZ();
122 final double rho0 = atmosphere.getDensity(date, posBody, atmFrame);
123 final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y, z), atmFrame) - rho0) / delta;
124 final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x, y + delta, z), atmFrame) - rho0) / delta;
125 final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x, y, z + delta), atmFrame) - rho0) / delta;
126 final double[] dXdQ = posBodyDS.getX().getAllDerivatives();
127 final double[] dYdQ = posBodyDS.getY().getAllDerivatives();
128 final double[] dZdQ = posBodyDS.getZ().getAllDerivatives();
129 final double[] rhoAll = new double[dXdQ.length];
130 rhoAll[0] = rho0;
131 for (int i = 1; i < rhoAll.length; ++i) {
132 rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
133 }
134 final DerivativeStructure rho = factory.build(rhoAll);
135
136
137
138
139
140
141 final FieldVector3D<DerivativeStructure> vAtmBodyDS =
142 new FieldVector3D<>(factory.constant(vAtmBody.getX()),
143 factory.constant(vAtmBody.getY()),
144 factory.constant(vAtmBody.getZ()));
145 final FieldPVCoordinates<DerivativeStructure> pvAtmBody = new FieldPVCoordinates<>(posBodyDS, vAtmBodyDS);
146 final FieldPVCoordinates<DerivativeStructure> pvAtm = toBody.getInverse().transformPVCoordinates(pvAtmBody);
147
148
149 final FieldVector3D<DerivativeStructure> relativeVelocity = pvAtm.getVelocity().subtract(velocity);
150
151
152 return spacecraft.dragAcceleration(state, rho, relativeVelocity,
153 forceModel.getParameters(factory.getDerivativeField()));
154
155 } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
156 return null;
157 }
158 }
159
160 @Override
161 protected FieldVector3D<Gradient> accelerationDerivativesGradient(final ForceModel forceModel,
162 final FieldSpacecraftState<Gradient> state) {
163 try {
164
165 final AbsoluteDate date = state.getDate().toAbsoluteDate();
166 final FieldVector3D<Gradient> position = state.getPVCoordinates().getPosition();
167 final FieldVector3D<Gradient> velocity = state.getPVCoordinates().getVelocity();
168 java.lang.reflect.Field atmosphereField = AbstractDragForceModel.class.getDeclaredField("atmosphere");
169 atmosphereField.setAccessible(true);
170 Atmosphere atmosphere = (Atmosphere) atmosphereField.get(forceModel);
171 java.lang.reflect.Field spacecraftField = DragForce.class.getDeclaredField("spacecraft");
172 spacecraftField.setAccessible(true);
173 DragSensitive spacecraft = (DragSensitive) spacecraftField.get(forceModel);
174
175 final int freeParameters = state.getMass().getFreeParameters();
176
177
178 final Frame atmFrame = atmosphere.getFrame();
179 final Transform toBody = state.getFrame().getTransformTo(atmFrame, date);
180 final FieldVector3D<Gradient> posBodyG = toBody.transformPosition(position);
181 final Vector3D posBody = posBodyG.toVector3D();
182 final Vector3D vAtmBody = atmosphere.getVelocity(date, posBody, atmFrame);
183
184
185
186
187
188
189 final double delta = 1.0;
190 final double x = posBody.getX();
191 final double y = posBody.getY();
192 final double z = posBody.getZ();
193 final double rho0 = atmosphere.getDensity(date, posBody, atmFrame);
194 final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y, z), atmFrame) - rho0) / delta;
195 final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x, y + delta, z), atmFrame) - rho0) / delta;
196 final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x, y, z + delta), atmFrame) - rho0) / delta;
197 final double[] dXdQ = posBodyG.getX().getGradient();
198 final double[] dYdQ = posBodyG.getY().getGradient();
199 final double[] dZdQ = posBodyG.getZ().getGradient();
200 final double[] rhoAll = new double[dXdQ.length];
201 for (int i = 0; i < rhoAll.length; ++i) {
202 rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
203 }
204 final Gradient rho = new Gradient(rho0, rhoAll);
205
206
207
208
209
210
211 final FieldVector3D<Gradient> vAtmBodyG =
212 new FieldVector3D<>(Gradient.constant(freeParameters, vAtmBody.getX()),
213 Gradient.constant(freeParameters, vAtmBody.getY()),
214 Gradient.constant(freeParameters, vAtmBody.getZ()));
215 final FieldPVCoordinates<Gradient> pvAtmBody = new FieldPVCoordinates<>(posBodyG, vAtmBodyG);
216 final FieldPVCoordinates<Gradient> pvAtm = toBody.getInverse().transformPVCoordinates(pvAtmBody);
217
218
219 final FieldVector3D<Gradient> relativeVelocity = pvAtm.getVelocity().subtract(velocity);
220
221
222 return spacecraft.dragAcceleration(state, rho, relativeVelocity,
223 forceModel.getParameters(GradientField.getField(freeParameters),
224 state.getDate()));
225
226 } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
227 return null;
228 }
229 }
230
231 @Test
232 void testParameterDerivativeSphere() {
233
234 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
235 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
236 final SpacecraftState state =
237 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
238 FramesFactory.getGCRF(),
239 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
240 Constants.EIGEN5C_EARTH_MU));
241
242 final DragForce forceModel =
243 new DragForce(getAtmosphere(),
244 new IsotropicDrag(2.5, 1.2));
245
246 Assertions.assertFalse(forceModel.dependsOnPositionOnly());
247
248 checkParameterDerivative(state, forceModel, DragSensitive.DRAG_COEFFICIENT, 1.0e-4, 2.0e-12);
249
250 }
251
252 @Test
253 void testParameterDerivativeSphereGradient() {
254
255 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
256 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
257 final SpacecraftState state =
258 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
259 FramesFactory.getGCRF(),
260 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
261 Constants.EIGEN5C_EARTH_MU));
262
263 final DragForce forceModel =
264 new DragForce(getAtmosphere(),
265 new IsotropicDrag(2.5, 1.2), true);
266
267 Assertions.assertFalse(forceModel.dependsOnPositionOnly());
268
269 checkParameterDerivativeGradient(state, forceModel, DragSensitive.DRAG_COEFFICIENT, 1.0e-4, 2.0e-12);
270
271 }
272
273 @Test
274 void testStateJacobianSphere() {
275
276
277 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
278 new TimeComponents(13, 59, 27.816),
279 TimeScalesFactory.getUTC());
280 double i = FastMath.toRadians(98.7);
281 double omega = FastMath.toRadians(93.0);
282 double OMEGA = FastMath.toRadians(15.0 * 22.5);
283 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
284 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
285 Constants.EIGEN5C_EARTH_MU);
286 OrbitType integrationType = OrbitType.CARTESIAN;
287 double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(0.01).getTolerances(orbit, integrationType);
288
289 NumericalPropagator propagator =
290 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
291 tolerances[0], tolerances[1]));
292 propagator.setOrbitType(integrationType);
293 final DragForce forceModel =
294 new DragForce(getAtmosphere(),
295 new IsotropicDrag(2.5, 1.2), false);
296 propagator.addForceModel(forceModel);
297 SpacecraftState state0 = new SpacecraftState(orbit);
298
299 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
300 1e3, tolerances[0], 2.0e-8);
301
302 }
303
304 @Test
305 void testParametersDerivativesBox() {
306
307 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
308 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
309 final SpacecraftState state =
310 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
311 FramesFactory.getGCRF(),
312 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
313 Constants.EIGEN5C_EARTH_MU));
314
315 final DragForce forceModel =
316 new DragForce(getAtmosphere(),
317 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
318 CelestialBodyFactory.getSun(), 20.0, Vector3D.PLUS_J,
319 1.2, 0.1, 0.7, 0.2));
320
321 checkParameterDerivative(state, forceModel, DragSensitive.GLOBAL_DRAG_FACTOR, 1.0e-4, 2.0e-11);
322
323 }
324
325 @Test
326 void testParametersDerivativesBoxGradient() {
327
328 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
329 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
330 final SpacecraftState state =
331 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
332 FramesFactory.getGCRF(),
333 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
334 Constants.EIGEN5C_EARTH_MU));
335
336 final DragForce forceModel =
337 new DragForce(getAtmosphere(),
338 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
339 CelestialBodyFactory.getSun(), 20.0, Vector3D.PLUS_J,
340 1.2, 0.1, 0.7, 0.2));
341
342 checkParameterDerivativeGradient(state, forceModel, DragSensitive.GLOBAL_DRAG_FACTOR, 1.0e-4, 2.0e-11);
343
344 }
345
346 @Test
347 void testJacobianBoxVs80Implementation() {
348
349
350 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
351 new TimeComponents(13, 59, 27.816),
352 TimeScalesFactory.getUTC());
353 double i = FastMath.toRadians(98.7);
354 double omega = FastMath.toRadians(93.0);
355 double OMEGA = FastMath.toRadians(15.0 * 22.5);
356 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
357 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
358 Constants.EIGEN5C_EARTH_MU);
359 final DragForce forceModel =
360 new DragForce(getAtmosphere(),
361 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
362 Vector3D.PLUS_J, 1.2, 0.0, 0.7, 0.2));
363 SpacecraftState state = new SpacecraftState(orbit,
364 Utils.defaultLaw().getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
365 checkStateJacobianVs80Implementation(state, forceModel,
366 new LofOffset(state.getFrame(), LOFType.LVLH_CCSDS),
367 5e-6, false);
368
369 }
370
371 @Test
372 void testJacobianBoxVs80ImplementationGradient() {
373
374
375 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
376 new TimeComponents(13, 59, 27.816),
377 TimeScalesFactory.getUTC());
378 double i = FastMath.toRadians(98.7);
379 double omega = FastMath.toRadians(93.0);
380 double OMEGA = FastMath.toRadians(15.0 * 22.5);
381 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
382 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
383 Constants.EIGEN5C_EARTH_MU);
384 final DragForce forceModel =
385 new DragForce(getAtmosphere(),
386 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
387 Vector3D.PLUS_J, 1.2, 0.0, 0.7, 0.2));
388 SpacecraftState state = new SpacecraftState(orbit,
389 Utils.defaultLaw().getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
390 checkStateJacobianVs80ImplementationGradient(state, forceModel,
391 new LofOffset(state.getFrame(), LOFType.LVLH_CCSDS),
392 5e-6, false);
393
394 }
395
396 @Test
397 void testJacobianBoxVsFiniteDifferences() {
398
399
400 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
401 new TimeComponents(13, 59, 27.816),
402 TimeScalesFactory.getUTC());
403 double i = FastMath.toRadians(98.7);
404 double omega = FastMath.toRadians(93.0);
405 double OMEGA = FastMath.toRadians(15.0 * 22.5);
406 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
407 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
408 Constants.EIGEN5C_EARTH_MU);
409 final DragForce forceModel =
410 new DragForce(getAtmosphere(),
411 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
412 Vector3D.PLUS_J, 1.2, 0.0, 0.7, 0.2));
413 SpacecraftState state = new SpacecraftState(orbit,
414 Utils.defaultLaw().getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
415 checkStateJacobianVsFiniteDifferences(state, forceModel, Utils.defaultLaw(), 1.0, 5.0e-6, false);
416
417 }
418
419 @Test
420 void testJacobianBoxVsFiniteDifferencesGradient() {
421
422
423 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
424 new TimeComponents(13, 59, 27.816),
425 TimeScalesFactory.getUTC());
426 double i = FastMath.toRadians(98.7);
427 double omega = FastMath.toRadians(93.0);
428 double OMEGA = FastMath.toRadians(15.0 * 22.5);
429 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
430 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
431 Constants.EIGEN5C_EARTH_MU);
432 final DragForce forceModel =
433 new DragForce(getAtmosphere(),
434 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
435 Vector3D.PLUS_J, 1.2, 0.0, 0.7, 0.2));
436 SpacecraftState state = new SpacecraftState(orbit,
437 Utils.defaultLaw().getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
438 checkStateJacobianVsFiniteDifferencesGradient(state, forceModel, Utils.defaultLaw(), 1.0, 5.0e-6, false);
439
440 }
441
442 @Test
443 void testGlobalStateJacobianBox() {
444
445
446 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
447 new TimeComponents(13, 59, 27.816),
448 TimeScalesFactory.getUTC());
449 double i = FastMath.toRadians(98.7);
450 double omega = FastMath.toRadians(93.0);
451 double OMEGA = FastMath.toRadians(15.0 * 22.5);
452 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
453 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
454 Constants.EIGEN5C_EARTH_MU);
455 OrbitType integrationType = OrbitType.CARTESIAN;
456 double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(0.01).getTolerances(orbit, integrationType);
457
458 NumericalPropagator propagator =
459 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
460 tolerances[0], tolerances[1]));
461 propagator.setOrbitType(integrationType);
462 final DragForce forceModel =
463 new DragForce(getAtmosphere(),
464 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
465 Vector3D.PLUS_J, 1.2, 0.0, 0.7, 0.2));
466 propagator.addForceModel(forceModel);
467 SpacecraftState state0 = new SpacecraftState(orbit);
468
469 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
470 1e3, tolerances[0], 3.0e-8);
471
472 }
473
474 @Test
475 void testIssue229() {
476 AbsoluteDate initialDate = new AbsoluteDate(2004, 1, 1, 0, 0, 0., TimeScalesFactory.getUTC());
477 Frame frame = FramesFactory.getEME2000();
478 double rpe = 160.e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
479 double rap = 2000.e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
480 double inc = FastMath.toRadians(0.);
481 double aop = FastMath.toRadians(0.);
482 double raan = FastMath.toRadians(0.);
483 double mean = FastMath.toRadians(180.);
484 double mass = 100.;
485 KeplerianOrbit orbit = new KeplerianOrbit(0.5 * (rpe + rap), (rap - rpe) / (rpe + rap),
486 inc, aop, raan, mean, PositionAngleType.MEAN,
487 frame, initialDate, Constants.EIGEN5C_EARTH_MU);
488
489 IsotropicDrag shape = new IsotropicDrag(10., 2.2);
490
491 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
492 BodyShape earthShape = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf);
493 Atmosphere atmosphere = new SimpleExponentialAtmosphere(earthShape, 2.6e-10, 200000, 26000);
494
495 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(0.1).getTolerances(orbit, OrbitType.CARTESIAN);
496 AbstractIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 300, tolerance[0], tolerance[1]);
497 NumericalPropagator propagator = new NumericalPropagator(integrator);
498 propagator.setOrbitType(OrbitType.CARTESIAN);
499 propagator.setMu(orbit.getMu());
500 propagator.addForceModel(new DragForce(atmosphere, shape));
501 MatricesHarvester harvester = propagator.setupMatricesComputation("partials", null, null);
502 propagator.setInitialState(new SpacecraftState(orbit, mass));
503
504 SpacecraftState state = propagator.propagate(new AbsoluteDate(2004, 1, 1, 1, 30, 0., TimeScalesFactory.getUTC()));
505 RealMatrix dYdY0 = harvester.getStateTransitionMatrix(state);
506
507 double delta = 0.1;
508 Orbit shifted = new CartesianOrbit(new TimeStampedPVCoordinates(orbit.getDate(),
509 orbit.getPosition().add(new Vector3D(delta, 0, 0)),
510 orbit.getPVCoordinates().getVelocity()),
511 orbit.getFrame(), orbit.getMu());
512 propagator.setInitialState(new SpacecraftState(shifted, mass));
513 SpacecraftState newState = propagator.propagate(new AbsoluteDate(2004, 1, 1, 1, 30, 0., TimeScalesFactory.getUTC()));
514 double[] dPVdX = new double[] {
515 (newState.getPosition().getX() - state.getPVCoordinates().getPosition().getX()) / delta,
516 (newState.getPosition().getY() - state.getPVCoordinates().getPosition().getY()) / delta,
517 (newState.getPosition().getZ() - state.getPVCoordinates().getPosition().getZ()) / delta,
518 (newState.getPVCoordinates().getVelocity().getX() - state.getPVCoordinates().getVelocity().getX()) / delta,
519 (newState.getPVCoordinates().getVelocity().getY() - state.getPVCoordinates().getVelocity().getY()) / delta,
520 (newState.getPVCoordinates().getVelocity().getZ() - state.getPVCoordinates().getVelocity().getZ()) / delta,
521 };
522
523 for (int i = 0; i < 6; ++i) {
524 Assertions.assertEquals(dPVdX[i], dYdY0.getEntry(i, 0), 6.2e-6 * FastMath.abs(dPVdX[i]));
525 }
526
527 }
528
529
530
531
532
533
534 @Test
535 public void RealFieldTest() {
536
537
538
539 DSFactory factory = new DSFactory(6, 4);
540 DerivativeStructure a_0 = factory.variable(0, 7e6);
541 DerivativeStructure e_0 = factory.variable(1, 0.01);
542 DerivativeStructure i_0 = factory.variable(2, 1.2);
543 DerivativeStructure R_0 = factory.variable(3, 0.7);
544 DerivativeStructure O_0 = factory.variable(4, 0.5);
545 DerivativeStructure n_0 = factory.variable(5, 0.1);
546
547 Field<DerivativeStructure> field = a_0.getField();
548 DerivativeStructure zero = field.getZero();
549
550
551 FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
552
553
554 Frame EME = FramesFactory.getEME2000();
555
556
557 FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
558 PositionAngleType.MEAN,
559 EME,
560 J2000,
561 zero.add(Constants.EIGEN5C_EARTH_MU));
562
563
564 FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
565 SpacecraftState iSR = initialState.toSpacecraftState();
566
567
568 ClassicalRungeKuttaFieldIntegrator<DerivativeStructure> integrator =
569 new ClassicalRungeKuttaFieldIntegrator<>(field, zero.add(6));
570 ClassicalRungeKuttaIntegrator RIntegrator =
571 new ClassicalRungeKuttaIntegrator(6);
572 OrbitType type = OrbitType.EQUINOCTIAL;
573
574
575 FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
576 FNP.setOrbitType(type);
577 FNP.setInitialState(initialState);
578
579 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
580 NP.setOrbitType(type);
581 NP.setInitialState(iSR);
582
583
584 final DragForce forceModel =
585 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
586 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
587 Constants.WGS84_EARTH_FLATTENING,
588 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
589 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
590 Vector3D.PLUS_J, 1.2, 0.0, 0.7, 0.2));
591 FNP.addForceModel(forceModel);
592 NP.addForceModel(forceModel);
593
594
595 checkRealFieldPropagation(FKO, PositionAngleType.MEAN, 1000., NP, FNP,
596 1.0e-30, 9.0e-9, 9.0e-11, 9.0e-11,
597 1, false);
598 }
599
600
601
602
603
604
605 @Test
606 public void RealFieldGradientTest() {
607
608
609
610 final int freeParameters = 6;
611 Gradient a_0 = Gradient.variable(freeParameters, 0, 7e6);
612 Gradient e_0 = Gradient.variable(freeParameters, 1, 0.01);
613 Gradient i_0 = Gradient.variable(freeParameters, 2, 1.2);
614 Gradient R_0 = Gradient.variable(freeParameters, 3, 0.7);
615 Gradient O_0 = Gradient.variable(freeParameters, 4, 0.5);
616 Gradient n_0 = Gradient.variable(freeParameters, 5, 0.1);
617
618 Field<Gradient> field = a_0.getField();
619 Gradient zero = field.getZero();
620
621
622 FieldAbsoluteDate<Gradient> J2000 = new FieldAbsoluteDate<>(field);
623
624
625 Frame EME = FramesFactory.getEME2000();
626
627
628 FieldKeplerianOrbit<Gradient> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
629 PositionAngleType.MEAN,
630 EME,
631 J2000,
632 zero.add(Constants.EIGEN5C_EARTH_MU));
633
634
635 FieldSpacecraftState<Gradient> initialState = new FieldSpacecraftState<>(FKO);
636 SpacecraftState iSR = initialState.toSpacecraftState();
637
638
639 ClassicalRungeKuttaFieldIntegrator<Gradient> integrator =
640 new ClassicalRungeKuttaFieldIntegrator<>(field, zero.add(6));
641 ClassicalRungeKuttaIntegrator RIntegrator =
642 new ClassicalRungeKuttaIntegrator(6);
643 OrbitType type = OrbitType.EQUINOCTIAL;
644
645
646 FieldNumericalPropagator<Gradient> FNP = new FieldNumericalPropagator<>(field, integrator);
647 FNP.setOrbitType(type);
648 FNP.setInitialState(initialState);
649
650 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
651 NP.setOrbitType(type);
652 NP.setInitialState(iSR);
653
654
655 final DragForce forceModel =
656 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
657 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
658 Constants.WGS84_EARTH_FLATTENING,
659 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
660 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
661 Vector3D.PLUS_J, 1.2, 0.0, 0.7, 0.2));
662 FNP.addForceModel(forceModel);
663 NP.addForceModel(forceModel);
664
665
666 checkRealFieldPropagationGradient(FKO, PositionAngleType.MEAN, 1000., NP, FNP,
667 1.0e-30, 3.2e-2, 7.7e-5, 2.8e-4,
668 1, false);
669 }
670
671
672
673
674
675
676 @Test
677 public void RealFieldExpectErrorTest() {
678 DSFactory factory = new DSFactory(6, 5);
679 DerivativeStructure a_0 = factory.variable(0, 7e6);
680 DerivativeStructure e_0 = factory.variable(1, 0.01);
681 DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
682 DerivativeStructure R_0 = factory.variable(3, 0.7);
683 DerivativeStructure O_0 = factory.variable(4, 0.5);
684 DerivativeStructure n_0 = factory.variable(5, 0.1);
685
686 Field<DerivativeStructure> field = a_0.getField();
687 DerivativeStructure zero = field.getZero();
688
689 FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
690
691 Frame EME = FramesFactory.getEME2000();
692
693 FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
694 PositionAngleType.MEAN,
695 EME,
696 J2000,
697 zero.add(Constants.EIGEN5C_EARTH_MU));
698
699 FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
700
701 SpacecraftState iSR = initialState.toSpacecraftState();
702 OrbitType type = OrbitType.KEPLERIAN;
703 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(FKO.toOrbit(), type);
704
705
706 AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
707 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
708 integrator.setInitialStepSize(60);
709 AdaptiveStepsizeIntegrator RIntegrator =
710 new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
711 RIntegrator.setInitialStepSize(60);
712
713 FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
714 FNP.setOrbitType(type);
715 FNP.setInitialState(initialState);
716
717 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
718 NP.setOrbitType(type);
719 NP.setInitialState(iSR);
720
721 final DragForce forceModel =
722 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
723 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
724 Constants.WGS84_EARTH_FLATTENING,
725 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
726 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
727 Vector3D.PLUS_J, 1.2, 0.0, 0.7, 0.2));
728 FNP.addForceModel(forceModel);
729
730
731 FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1000.);
732 FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
733 SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
734 FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
735 PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
736
737 Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getX() - finPVC_R.getPosition().getX()) < FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
738 Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getY() - finPVC_R.getPosition().getY()) < FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
739 Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getZ() - finPVC_R.getPosition().getZ()) < FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
740 }
741
742 @Test
743 void testDependsOnAttitudeRateTrue() {
744
745 final DragSensitive mockedSensitive = Mockito.mock(DragSensitive.class);
746 Mockito.when(mockedSensitive.dependsOnAttitudeRate()).thenReturn(true);
747 final DragForce dragForce = new DragForce(Mockito.mock(Atmosphere.class), mockedSensitive);
748
749 final boolean value = dragForce.dependsOnAttitudeRate();
750
751 Assertions.assertTrue(value);
752 }
753
754 @Test
755 void testDependsOnAttitudeRateFalse() {
756
757 final DragSensitive mockedSensitive = Mockito.mock(DragSensitive.class);
758 Mockito.when(mockedSensitive.dependsOnAttitudeRate()).thenReturn(false);
759 final DragForce dragForce = new DragForce(Mockito.mock(Atmosphere.class), mockedSensitive);
760
761 final boolean value = dragForce.dependsOnAttitudeRate();
762
763 Assertions.assertFalse(value);
764 }
765
766
767
768
769
770 @Test
771 void testGetParameterDriversSphereForParameterWithSeveralValues() {
772
773
774 final Atmosphere atmosphere = new HarrisPriester(CelestialBodyFactory.getSun(),
775 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
776 Constants.WGS84_EARTH_FLATTENING,
777 FramesFactory.getITRF(IERSConventions.IERS_2010, true)));
778
779
780 AbsoluteDate date = new AbsoluteDate("2000-01-01T00:00:00.000", TimeScalesFactory.getUTC());
781
782
783
784
785 double dragArea = 2.;
786 double dragCd0 = 0.;
787 DragForce forceModel = new DragForce(atmosphere, new IsotropicDrag(dragArea, dragCd0));
788 Assertions.assertFalse(forceModel.dependsOnPositionOnly());
789 List<ParameterDriver> drivers = forceModel.getParametersDrivers();
790 Assertions.assertEquals(2, drivers.size());
791 Assertions.assertEquals(1.0, drivers.get(0).getValue(), 0.);
792 Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR, drivers.get(0).getName());
793 Assertions.assertEquals(dragCd0, drivers.get(1).getValue(), 0.);
794 Assertions.assertEquals(DragSensitive.DRAG_COEFFICIENT, drivers.get(1).getName());
795
796
797 IsotropicDrag isoDrag = (IsotropicDrag) forceModel.getSpacecraft();
798 drivers = isoDrag.getDragParametersDrivers();
799 Assertions.assertEquals(2, drivers.size());
800 Assertions.assertEquals(1.0, drivers.get(0).getValue(new AbsoluteDate()), 0.);
801 Assertions.assertEquals(DragSensitive.GLOBAL_DRAG_FACTOR, drivers.get(0).getName());
802 Assertions.assertEquals(dragCd0, drivers.get(1).getValue(new AbsoluteDate()), 0.);
803 Assertions.assertEquals(DragSensitive.DRAG_COEFFICIENT, drivers.get(1).getName());
804
805
806
807 double dragCd1 = 1.;
808 double dragCd2 = 2.;
809 double dt = 120.;
810
811 isoDrag = new IsotropicDrag(dragArea, dragCd0);
812 isoDrag.getDragParametersDrivers().get(0).addSpans(date.shiftedBy(-3*dt), date.shiftedBy(2.0*dt), 2*dt);
813 isoDrag.getDragParametersDrivers().get(0).setValue(dragCd2, date.shiftedBy(-2*dt));
814 isoDrag.getDragParametersDrivers().get(0).setValue(dragCd0, date.shiftedBy(-dt));
815 isoDrag.getDragParametersDrivers().get(0).setValue(dragCd1, date.shiftedBy(dt));
816
817 forceModel = new DragForce(atmosphere, isoDrag);
818
819 drivers = forceModel.getParametersDrivers();
820 int nnb = 0;
821 Assertions.assertEquals(3, drivers.get(0).getNbOfValues());
822 for (Span<String> span = isoDrag.getDragParametersDrivers().get(0).getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
823 Assertions.assertEquals("Span" + drivers.get(0).getName() + Integer.toString(nnb++),
824 span.getData());
825 }
826
827
828
829 double eps = 1.e-14;
830
831 Assertions.assertEquals(dragCd2, drivers.get(0).getValue(date.shiftedBy(-2 * dt)), 0.);
832 Assertions.assertEquals(dragCd2, drivers.get(0).getValue(date.shiftedBy(-dt - eps)), 0.);
833 Assertions.assertEquals(dragCd2, drivers.get(0).getValue(date.shiftedBy(-dt - 86400.)), 0.);
834
835 Assertions.assertEquals(dragCd0, drivers.get(0).getValue(date), 0.);
836 Assertions.assertEquals(dragCd0, drivers.get(0).getValue(date.shiftedBy(dt - eps)), 0.);
837 Assertions.assertEquals(dragCd0, drivers.get(0).getValue(date.shiftedBy(-dt)), 0.);
838
839 Assertions.assertEquals(dragCd1, drivers.get(0).getValue(date.shiftedBy(2 * dt)), 0.);
840 Assertions.assertEquals(dragCd1, drivers.get(0).getValue(date.shiftedBy(dt + eps)), 0.);
841 Assertions.assertEquals(dragCd1, drivers.get(0).getValue(date.shiftedBy(dt + 86400.)), 0.);
842
843 }
844
845
846
847
848
849
850 @Test
851 void testParameterDerivativeSphereForParameterWithSeveralValues() {
852
853
854 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
855 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
856 final AbsoluteDate date = new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI());
857 final SpacecraftState state =
858 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
859 FramesFactory.getGCRF(),
860 date,
861 Constants.EIGEN5C_EARTH_MU));
862
863
864 final Atmosphere atmosphere = new HarrisPriester(CelestialBodyFactory.getSun(),
865 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
866 Constants.WGS84_EARTH_FLATTENING,
867 FramesFactory.getITRF(IERSConventions.IERS_2010, true)));
868
869
870 final double dragArea = 2.5;
871
872
873 final double dragCd = 1.2;
874 final IsotropicDrag isotropicDrag = new IsotropicDrag(dragArea, dragCd);
875 isotropicDrag.getDragParametersDrivers().get(0).setName("Cd");
876
877
878 final double dragCd2 = 3.;
879 final double dt2 = 4 * 3600.;
880 final AbsoluteDate date2 = date.shiftedBy(dt2);
881 isotropicDrag.getDragParametersDrivers().get(0).getValueSpanMap().addValidAfter(dragCd2, date2, false);
882 isotropicDrag.getDragParametersDrivers().get(0).getNamesSpanMap().addValidAfter("Cd2", date2, false);
883 isotropicDrag.getDragParametersDrivers().get(0).getNamesSpanMap().addValidBefore("Cd", date2, false);
884
885
886 final double dragCd3 = 3.;
887 final double dt3 = -86400.;
888 final AbsoluteDate date3 = date.shiftedBy(dt3);
889 isotropicDrag.getDragParametersDrivers().get(0).getValueSpanMap().addValidAfter(dragCd3, date3, false);
890 isotropicDrag.getDragParametersDrivers().get(0).getNamesSpanMap().addValidAfter("Cd3", date3, false);
891
892
893 final DragForce forceModel = new DragForce(atmosphere, isotropicDrag);
894
895 Assertions.assertFalse(forceModel.dependsOnPositionOnly());
896
897
898 checkParameterDerivative(state, forceModel, "Cd" , 1.0e-4, 2.0e-12);
899
900
901 checkParameterDerivative(state.shiftedBy(dt2 * 1.1), forceModel, "Cd", 1.0e-4, 2.0e-12);
902
903
904 checkParameterDerivative(state.shiftedBy(dt3 * 1.1), forceModel, "Cd", 1.0e-4, 2.0e-12);
905 }
906
907 private Atmosphere getAtmosphere() {
908 return new HarrisPriester(CelestialBodyFactory.getSun(),
909 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
910 Constants.WGS84_EARTH_FLATTENING,
911 FramesFactory.getITRF(IERSConventions.IERS_2010, true)));
912 }
913
914 @BeforeEach
915 public void setUp() {
916 Utils.setDataRoot("regular-data");
917 }
918
919 }