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