1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.drag;
18
19
20 import org.hipparchus.Field;
21 import org.hipparchus.analysis.differentiation.DSFactory;
22 import org.hipparchus.analysis.differentiation.DerivativeStructure;
23 import org.hipparchus.analysis.differentiation.Gradient;
24 import org.hipparchus.analysis.differentiation.GradientField;
25 import org.hipparchus.geometry.euclidean.threed.FieldRotation;
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.Assert;
38 import org.junit.Before;
39 import org.junit.Test;
40 import org.orekit.Utils;
41 import org.orekit.attitudes.AttitudeProvider;
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.PositionAngle;
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.TimeStampedPVCoordinates;
79
80 public class DragForceTest extends AbstractLegacyForceModelTest {
81
82 private static final AttitudeProvider DEFAULT_LAW = Utils.defaultLaw();
83
84 @Override
85 protected FieldVector3D<DerivativeStructure> accelerationDerivatives(final ForceModel forceModel,
86 final AbsoluteDate date, final Frame frame,
87 final FieldVector3D<DerivativeStructure> position,
88 final FieldVector3D<DerivativeStructure> velocity,
89 final FieldRotation<DerivativeStructure> rotation,
90 final DerivativeStructure mass) {
91 try {
92
93 java.lang.reflect.Field atmosphereField = DragForce.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 = mass.getFactory();
102
103
104 final Frame atmFrame = atmosphere.getFrame();
105 final Transform toBody = frame.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(new FieldAbsoluteDate<>(factory.getDerivativeField(), date),
153 frame, position, rotation, mass, rho, relativeVelocity,
154 forceModel.getParameters(factory.getDerivativeField()));
155
156 } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
157 return null;
158 }
159 }
160
161 @Override
162 protected FieldVector3D<Gradient> accelerationDerivativesGradient(final ForceModel forceModel,
163 final AbsoluteDate date, final Frame frame,
164 final FieldVector3D<Gradient> position,
165 final FieldVector3D<Gradient> velocity,
166 final FieldRotation<Gradient> rotation,
167 final Gradient mass) {
168 try {
169
170 java.lang.reflect.Field atmosphereField = DragForce.class.getDeclaredField("atmosphere");
171 atmosphereField.setAccessible(true);
172 Atmosphere atmosphere = (Atmosphere) atmosphereField.get(forceModel);
173 java.lang.reflect.Field spacecraftField = DragForce.class.getDeclaredField("spacecraft");
174 spacecraftField.setAccessible(true);
175 DragSensitive spacecraft = (DragSensitive) spacecraftField.get(forceModel);
176
177 final int freeParameters = mass.getFreeParameters();
178
179
180 final Frame atmFrame = atmosphere.getFrame();
181 final Transform toBody = frame.getTransformTo(atmFrame, date);
182 final FieldVector3D<Gradient> posBodyG = toBody.transformPosition(position);
183 final Vector3D posBody = posBodyG.toVector3D();
184 final Vector3D vAtmBody = atmosphere.getVelocity(date, posBody, atmFrame);
185
186
187
188
189
190
191 final double delta = 1.0;
192 final double x = posBody.getX();
193 final double y = posBody.getY();
194 final double z = posBody.getZ();
195 final double rho0 = atmosphere.getDensity(date, posBody, atmFrame);
196 final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y, z), atmFrame) - rho0) / delta;
197 final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x, y + delta, z), atmFrame) - rho0) / delta;
198 final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x, y, z + delta), atmFrame) - rho0) / delta;
199 final double[] dXdQ = posBodyG.getX().getGradient();
200 final double[] dYdQ = posBodyG.getY().getGradient();
201 final double[] dZdQ = posBodyG.getZ().getGradient();
202 final double[] rhoAll = new double[dXdQ.length];
203 for (int i = 0; i < rhoAll.length; ++i) {
204 rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
205 }
206 final Gradient rho = new Gradient(rho0, rhoAll);
207
208
209
210
211
212
213 final FieldVector3D<Gradient> vAtmBodyG =
214 new FieldVector3D<>(Gradient.constant(freeParameters, vAtmBody.getX()),
215 Gradient.constant(freeParameters, vAtmBody.getY()),
216 Gradient.constant(freeParameters, vAtmBody.getZ()));
217 final FieldPVCoordinates<Gradient> pvAtmBody = new FieldPVCoordinates<>(posBodyG, vAtmBodyG);
218 final FieldPVCoordinates<Gradient> pvAtm = toBody.getInverse().transformPVCoordinates(pvAtmBody);
219
220
221 final FieldVector3D<Gradient> relativeVelocity = pvAtm.getVelocity().subtract(velocity);
222
223
224 return spacecraft.dragAcceleration(new FieldAbsoluteDate<>(GradientField.getField(freeParameters), date),
225 frame, position, rotation, mass, rho, relativeVelocity,
226 forceModel.getParameters(GradientField.getField(freeParameters)));
227
228 } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
229 return null;
230 }
231 }
232
233 @Test
234 public void testParameterDerivativeSphere() {
235
236 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
237 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
238 final SpacecraftState state =
239 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
240 FramesFactory.getGCRF(),
241 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
242 Constants.EIGEN5C_EARTH_MU));
243
244 final DragForce forceModel =
245 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
246 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
247 Constants.WGS84_EARTH_FLATTENING,
248 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
249 new IsotropicDrag(2.5, 1.2));
250
251 Assert.assertFalse(forceModel.dependsOnPositionOnly());
252
253 checkParameterDerivative(state, forceModel, DragSensitive.DRAG_COEFFICIENT, 1.0e-4, 2.0e-12);
254
255 }
256
257 @Test
258 public void testParameterDerivativeSphereGradient() {
259
260 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
261 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
262 final SpacecraftState state =
263 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
264 FramesFactory.getGCRF(),
265 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
266 Constants.EIGEN5C_EARTH_MU));
267
268 final DragForce forceModel =
269 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
270 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
271 Constants.WGS84_EARTH_FLATTENING,
272 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
273 new IsotropicDrag(2.5, 1.2));
274
275 Assert.assertFalse(forceModel.dependsOnPositionOnly());
276
277 checkParameterDerivativeGradient(state, forceModel, DragSensitive.DRAG_COEFFICIENT, 1.0e-4, 2.0e-12);
278
279 }
280
281 @Test
282 public void testStateJacobianSphere() {
283
284
285 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
286 new TimeComponents(13, 59, 27.816),
287 TimeScalesFactory.getUTC());
288 double i = FastMath.toRadians(98.7);
289 double omega = FastMath.toRadians(93.0);
290 double OMEGA = FastMath.toRadians(15.0 * 22.5);
291 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
292 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
293 Constants.EIGEN5C_EARTH_MU);
294 OrbitType integrationType = OrbitType.CARTESIAN;
295 double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
296
297 NumericalPropagator propagator =
298 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
299 tolerances[0], tolerances[1]));
300 propagator.setOrbitType(integrationType);
301 final DragForce forceModel =
302 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
303 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
304 Constants.WGS84_EARTH_FLATTENING,
305 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
306 new IsotropicDrag(2.5, 1.2));
307 propagator.addForceModel(forceModel);
308 SpacecraftState state0 = new SpacecraftState(orbit);
309
310 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
311 1e3, tolerances[0], 2.0e-8);
312
313 }
314
315 @Test
316 public void testParametersDerivativesBox() {
317
318 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
319 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
320 final SpacecraftState state =
321 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
322 FramesFactory.getGCRF(),
323 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
324 Constants.EIGEN5C_EARTH_MU));
325
326 final DragForce forceModel =
327 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
328 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
329 Constants.WGS84_EARTH_FLATTENING,
330 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
331 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
332 CelestialBodyFactory.getSun(), 20.0, Vector3D.PLUS_J,
333 1.2, 0.1, 0.7, 0.2));
334
335 checkParameterDerivative(state, forceModel, DragSensitive.DRAG_COEFFICIENT, 1.0e-4, 5.0e-13);
336 checkParameterDerivative(state, forceModel, DragSensitive.LIFT_RATIO, 1.0e-4, 2.0e-11);
337
338 }
339
340 @Test
341 public void testParametersDerivativesBoxGradient() {
342
343 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
344 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
345 final SpacecraftState state =
346 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
347 FramesFactory.getGCRF(),
348 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
349 Constants.EIGEN5C_EARTH_MU));
350
351 final DragForce forceModel =
352 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
353 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
354 Constants.WGS84_EARTH_FLATTENING,
355 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
356 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8,
357 CelestialBodyFactory.getSun(), 20.0, Vector3D.PLUS_J,
358 1.2, 0.1, 0.7, 0.2));
359
360 checkParameterDerivativeGradient(state, forceModel, DragSensitive.DRAG_COEFFICIENT, 1.0e-4, 5.0e-13);
361 checkParameterDerivativeGradient(state, forceModel, DragSensitive.LIFT_RATIO, 1.0e-4, 2.0e-11);
362
363 }
364
365 @Test
366 public void testJacobianBoxVs80Implementation() {
367
368
369 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
370 new TimeComponents(13, 59, 27.816),
371 TimeScalesFactory.getUTC());
372 double i = FastMath.toRadians(98.7);
373 double omega = FastMath.toRadians(93.0);
374 double OMEGA = FastMath.toRadians(15.0 * 22.5);
375 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
376 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
377 Constants.EIGEN5C_EARTH_MU);
378 final DragForce forceModel =
379 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
380 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
381 Constants.WGS84_EARTH_FLATTENING,
382 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
383 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
384 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
385 SpacecraftState state = new SpacecraftState(orbit,
386 Utils.defaultLaw().getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
387 checkStateJacobianVs80Implementation(state, forceModel,
388 new LofOffset(state.getFrame(), LOFType.LVLH_CCSDS),
389 5e-6, false);
390
391 }
392
393 @Test
394 public void testJacobianBoxVs80ImplementationGradient() {
395
396
397 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
398 new TimeComponents(13, 59, 27.816),
399 TimeScalesFactory.getUTC());
400 double i = FastMath.toRadians(98.7);
401 double omega = FastMath.toRadians(93.0);
402 double OMEGA = FastMath.toRadians(15.0 * 22.5);
403 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
404 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
405 Constants.EIGEN5C_EARTH_MU);
406 final DragForce forceModel =
407 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
408 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
409 Constants.WGS84_EARTH_FLATTENING,
410 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
411 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
412 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
413 SpacecraftState state = new SpacecraftState(orbit,
414 Utils.defaultLaw().getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
415 checkStateJacobianVs80ImplementationGradient(state, forceModel,
416 new LofOffset(state.getFrame(), LOFType.LVLH_CCSDS),
417 5e-6, false);
418
419 }
420
421 @Test
422 public void testJacobianBoxVsFiniteDifferences() {
423
424
425 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
426 new TimeComponents(13, 59, 27.816),
427 TimeScalesFactory.getUTC());
428 double i = FastMath.toRadians(98.7);
429 double omega = FastMath.toRadians(93.0);
430 double OMEGA = FastMath.toRadians(15.0 * 22.5);
431 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
432 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
433 Constants.EIGEN5C_EARTH_MU);
434 final DragForce forceModel =
435 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
436 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
437 Constants.WGS84_EARTH_FLATTENING,
438 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
439 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
440 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
441 SpacecraftState state = new SpacecraftState(orbit,
442 Utils.defaultLaw().getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
443 checkStateJacobianVsFiniteDifferences(state, forceModel, Utils.defaultLaw(), 1.0, 5.0e-6, false);
444
445 }
446
447 @Test
448 public void testJacobianBoxVsFiniteDifferencesGradient() {
449
450
451 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
452 new TimeComponents(13, 59, 27.816),
453 TimeScalesFactory.getUTC());
454 double i = FastMath.toRadians(98.7);
455 double omega = FastMath.toRadians(93.0);
456 double OMEGA = FastMath.toRadians(15.0 * 22.5);
457 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
458 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
459 Constants.EIGEN5C_EARTH_MU);
460 final DragForce forceModel =
461 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
462 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
463 Constants.WGS84_EARTH_FLATTENING,
464 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
465 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
466 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
467 SpacecraftState state = new SpacecraftState(orbit,
468 DEFAULT_LAW.getAttitude(orbit, orbit.getDate(), orbit.getFrame()));
469 checkStateJacobianVsFiniteDifferencesGradient(state, forceModel, Utils.defaultLaw(), 1.0, 5.0e-6, false);
470
471 }
472
473 @Test
474 public void testGlobalStateJacobianBox() {
475
476
477 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
478 new TimeComponents(13, 59, 27.816),
479 TimeScalesFactory.getUTC());
480 double i = FastMath.toRadians(98.7);
481 double omega = FastMath.toRadians(93.0);
482 double OMEGA = FastMath.toRadians(15.0 * 22.5);
483 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
484 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
485 Constants.EIGEN5C_EARTH_MU);
486 OrbitType integrationType = OrbitType.CARTESIAN;
487 double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
488
489 NumericalPropagator propagator =
490 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
491 tolerances[0], tolerances[1]));
492 propagator.setOrbitType(integrationType);
493 final DragForce forceModel =
494 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
495 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
496 Constants.WGS84_EARTH_FLATTENING,
497 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
498 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
499 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
500 propagator.addForceModel(forceModel);
501 SpacecraftState state0 = new SpacecraftState(orbit);
502
503 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
504 1e3, tolerances[0], 3.0e-8);
505
506 }
507
508 @Test
509 public void testIssue229() {
510 AbsoluteDate initialDate = new AbsoluteDate(2004, 1, 1, 0, 0, 0., TimeScalesFactory.getUTC());
511 Frame frame = FramesFactory.getEME2000();
512 double rpe = 160.e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
513 double rap = 2000.e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
514 double inc = FastMath.toRadians(0.);
515 double aop = FastMath.toRadians(0.);
516 double raan = FastMath.toRadians(0.);
517 double mean = FastMath.toRadians(180.);
518 double mass = 100.;
519 KeplerianOrbit orbit = new KeplerianOrbit(0.5 * (rpe + rap), (rap - rpe) / (rpe + rap),
520 inc, aop, raan, mean, PositionAngle.MEAN,
521 frame, initialDate, Constants.EIGEN5C_EARTH_MU);
522
523 IsotropicDrag shape = new IsotropicDrag(10., 2.2);
524
525 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
526 BodyShape earthShape = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf);
527 Atmosphere atmosphere = new SimpleExponentialAtmosphere(earthShape, 2.6e-10, 200000, 26000);
528
529 double[][] tolerance = NumericalPropagator.tolerances(0.1, orbit, OrbitType.CARTESIAN);
530 AbstractIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 300, tolerance[0], tolerance[1]);
531 NumericalPropagator propagator = new NumericalPropagator(integrator);
532 propagator.setOrbitType(OrbitType.CARTESIAN);
533 propagator.setMu(orbit.getMu());
534 propagator.addForceModel(new DragForce(atmosphere, shape));
535 MatricesHarvester harvester = propagator.setupMatricesComputation("partials", null, null);
536 propagator.setInitialState(new SpacecraftState(orbit, mass));
537
538 SpacecraftState state = propagator.propagate(new AbsoluteDate(2004, 1, 1, 1, 30, 0., TimeScalesFactory.getUTC()));
539 RealMatrix dYdY0 = harvester.getStateTransitionMatrix(state);
540
541 double delta = 0.1;
542 Orbit shifted = new CartesianOrbit(new TimeStampedPVCoordinates(orbit.getDate(),
543 orbit.getPVCoordinates().getPosition().add(new Vector3D(delta, 0, 0)),
544 orbit.getPVCoordinates().getVelocity()),
545 orbit.getFrame(), orbit.getMu());
546 propagator.setInitialState(new SpacecraftState(shifted, mass));
547 SpacecraftState newState = propagator.propagate(new AbsoluteDate(2004, 1, 1, 1, 30, 0., TimeScalesFactory.getUTC()));
548 double[] dPVdX = new double[] {
549 (newState.getPVCoordinates().getPosition().getX() - state.getPVCoordinates().getPosition().getX()) / delta,
550 (newState.getPVCoordinates().getPosition().getY() - state.getPVCoordinates().getPosition().getY()) / delta,
551 (newState.getPVCoordinates().getPosition().getZ() - state.getPVCoordinates().getPosition().getZ()) / delta,
552 (newState.getPVCoordinates().getVelocity().getX() - state.getPVCoordinates().getVelocity().getX()) / delta,
553 (newState.getPVCoordinates().getVelocity().getY() - state.getPVCoordinates().getVelocity().getY()) / delta,
554 (newState.getPVCoordinates().getVelocity().getZ() - state.getPVCoordinates().getVelocity().getZ()) / delta,
555 };
556
557 for (int i = 0; i < 6; ++i) {
558 Assert.assertEquals(dPVdX[i], dYdY0.getEntry(i, 0), 6.2e-6 * FastMath.abs(dPVdX[i]));
559 }
560
561 }
562
563
564
565
566
567
568 @Test
569 public void RealFieldTest() {
570
571
572
573 DSFactory factory = new DSFactory(6, 4);
574 DerivativeStructure a_0 = factory.variable(0, 7e6);
575 DerivativeStructure e_0 = factory.variable(1, 0.01);
576 DerivativeStructure i_0 = factory.variable(2, 1.2);
577 DerivativeStructure R_0 = factory.variable(3, 0.7);
578 DerivativeStructure O_0 = factory.variable(4, 0.5);
579 DerivativeStructure n_0 = factory.variable(5, 0.1);
580
581 Field<DerivativeStructure> field = a_0.getField();
582 DerivativeStructure zero = field.getZero();
583
584
585 FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
586
587
588 Frame EME = FramesFactory.getEME2000();
589
590
591 FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
592 PositionAngle.MEAN,
593 EME,
594 J2000,
595 zero.add(Constants.EIGEN5C_EARTH_MU));
596
597
598 FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
599 SpacecraftState iSR = initialState.toSpacecraftState();
600
601
602 ClassicalRungeKuttaFieldIntegrator<DerivativeStructure> integrator =
603 new ClassicalRungeKuttaFieldIntegrator<>(field, zero.add(6));
604 ClassicalRungeKuttaIntegrator RIntegrator =
605 new ClassicalRungeKuttaIntegrator(6);
606 OrbitType type = OrbitType.EQUINOCTIAL;
607
608
609 FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
610 FNP.setOrbitType(type);
611 FNP.setInitialState(initialState);
612
613 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
614 NP.setOrbitType(type);
615 NP.setInitialState(iSR);
616
617
618 final DragForce forceModel =
619 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
620 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
621 Constants.WGS84_EARTH_FLATTENING,
622 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
623 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
624 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
625 FNP.addForceModel(forceModel);
626 NP.addForceModel(forceModel);
627
628
629 checkRealFieldPropagation(FKO, PositionAngle.MEAN, 1000., NP, FNP,
630 1.0e-30, 9.0e-9, 9.0e-11, 9.0e-11,
631 1, false);
632 }
633
634
635
636
637
638
639 @Test
640 public void RealFieldGradientTest() {
641
642
643
644 final int freeParameters = 6;
645 Gradient a_0 = Gradient.variable(freeParameters, 0, 7e6);
646 Gradient e_0 = Gradient.variable(freeParameters, 1, 0.01);
647 Gradient i_0 = Gradient.variable(freeParameters, 2, 1.2);
648 Gradient R_0 = Gradient.variable(freeParameters, 3, 0.7);
649 Gradient O_0 = Gradient.variable(freeParameters, 4, 0.5);
650 Gradient n_0 = Gradient.variable(freeParameters, 5, 0.1);
651
652 Field<Gradient> field = a_0.getField();
653 Gradient zero = field.getZero();
654
655
656 FieldAbsoluteDate<Gradient> J2000 = new FieldAbsoluteDate<>(field);
657
658
659 Frame EME = FramesFactory.getEME2000();
660
661
662 FieldKeplerianOrbit<Gradient> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
663 PositionAngle.MEAN,
664 EME,
665 J2000,
666 zero.add(Constants.EIGEN5C_EARTH_MU));
667
668
669 FieldSpacecraftState<Gradient> initialState = new FieldSpacecraftState<>(FKO);
670 SpacecraftState iSR = initialState.toSpacecraftState();
671
672
673 ClassicalRungeKuttaFieldIntegrator<Gradient> integrator =
674 new ClassicalRungeKuttaFieldIntegrator<>(field, zero.add(6));
675 ClassicalRungeKuttaIntegrator RIntegrator =
676 new ClassicalRungeKuttaIntegrator(6);
677 OrbitType type = OrbitType.EQUINOCTIAL;
678
679
680 FieldNumericalPropagator<Gradient> FNP = new FieldNumericalPropagator<>(field, integrator);
681 FNP.setOrbitType(type);
682 FNP.setInitialState(initialState);
683
684 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
685 NP.setOrbitType(type);
686 NP.setInitialState(iSR);
687
688
689 final DragForce forceModel =
690 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
691 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
692 Constants.WGS84_EARTH_FLATTENING,
693 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
694 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
695 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
696 FNP.addForceModel(forceModel);
697 NP.addForceModel(forceModel);
698
699
700 checkRealFieldPropagationGradient(FKO, PositionAngle.MEAN, 1000., NP, FNP,
701 1.0e-30, 3.2e-2, 7.7e-5, 2.8e-4,
702 1, false);
703 }
704
705
706
707
708
709
710 @Test
711 public void RealFieldExpectErrorTest() {
712 DSFactory factory = new DSFactory(6, 5);
713 DerivativeStructure a_0 = factory.variable(0, 7e6);
714 DerivativeStructure e_0 = factory.variable(1, 0.01);
715 DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
716 DerivativeStructure R_0 = factory.variable(3, 0.7);
717 DerivativeStructure O_0 = factory.variable(4, 0.5);
718 DerivativeStructure n_0 = factory.variable(5, 0.1);
719
720 Field<DerivativeStructure> field = a_0.getField();
721 DerivativeStructure zero = field.getZero();
722
723 FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
724
725 Frame EME = FramesFactory.getEME2000();
726
727 FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
728 PositionAngle.MEAN,
729 EME,
730 J2000,
731 zero.add(Constants.EIGEN5C_EARTH_MU));
732
733 FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
734
735 SpacecraftState iSR = initialState.toSpacecraftState();
736 OrbitType type = OrbitType.KEPLERIAN;
737 double[][] tolerance = NumericalPropagator.tolerances(10.0, FKO.toOrbit(), type);
738
739
740 AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
741 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
742 integrator.setInitialStepSize(60);
743 AdaptiveStepsizeIntegrator RIntegrator =
744 new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
745 RIntegrator.setInitialStepSize(60);
746
747 FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
748 FNP.setOrbitType(type);
749 FNP.setInitialState(initialState);
750
751 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
752 NP.setOrbitType(type);
753 NP.setInitialState(iSR);
754
755 final DragForce forceModel =
756 new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(),
757 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
758 Constants.WGS84_EARTH_FLATTENING,
759 FramesFactory.getITRF(IERSConventions.IERS_2010, true))),
760 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
761 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
762 FNP.addForceModel(forceModel);
763
764
765 FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1000.);
766 FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
767 SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
768 FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
769 PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
770
771 Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getX() - finPVC_R.getPosition().getX()) < FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
772 Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getY() - finPVC_R.getPosition().getY()) < FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
773 Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getZ() - finPVC_R.getPosition().getZ()) < FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
774 }
775
776 @Before
777 public void setUp() {
778 Utils.setDataRoot("regular-data");
779 }
780
781 }
782
783