1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.radiation;
18
19
20 import java.io.FileNotFoundException;
21 import java.lang.reflect.InvocationTargetException;
22 import java.text.ParseException;
23
24 import org.hipparchus.Field;
25 import org.hipparchus.analysis.differentiation.DSFactory;
26 import org.hipparchus.analysis.differentiation.DerivativeStructure;
27 import org.hipparchus.analysis.differentiation.Gradient;
28 import org.hipparchus.geometry.euclidean.threed.FieldRotation;
29 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
30 import org.hipparchus.geometry.euclidean.threed.Vector3D;
31 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
32 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
33 import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
34 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
35 import org.hipparchus.util.Decimal64Field;
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.LofOffset;
42 import org.orekit.bodies.CelestialBodyFactory;
43 import org.orekit.bodies.OneAxisEllipsoid;
44 import org.orekit.errors.OrekitException;
45 import org.orekit.forces.AbstractLegacyForceModelTest;
46 import org.orekit.forces.BoxAndSolarArraySpacecraft;
47 import org.orekit.forces.ForceModel;
48 import org.orekit.frames.Frame;
49 import org.orekit.frames.FramesFactory;
50 import org.orekit.frames.LOFType;
51 import org.orekit.orbits.CartesianOrbit;
52 import org.orekit.orbits.EquinoctialOrbit;
53 import org.orekit.orbits.FieldKeplerianOrbit;
54 import org.orekit.orbits.KeplerianOrbit;
55 import org.orekit.orbits.Orbit;
56 import org.orekit.orbits.OrbitType;
57 import org.orekit.orbits.PositionAngle;
58 import org.orekit.propagation.FieldSpacecraftState;
59 import org.orekit.propagation.Propagator;
60 import org.orekit.propagation.SpacecraftState;
61 import org.orekit.propagation.analytical.KeplerianPropagator;
62 import org.orekit.propagation.numerical.FieldNumericalPropagator;
63 import org.orekit.propagation.numerical.NumericalPropagator;
64 import org.orekit.propagation.sampling.OrekitFixedStepHandler;
65 import org.orekit.time.AbsoluteDate;
66 import org.orekit.time.DateComponents;
67 import org.orekit.time.FieldAbsoluteDate;
68 import org.orekit.time.TimeComponents;
69 import org.orekit.time.TimeScalesFactory;
70 import org.orekit.utils.Constants;
71 import org.orekit.utils.ExtendedPVCoordinatesProvider;
72 import org.orekit.utils.FieldPVCoordinates;
73 import org.orekit.utils.IERSConventions;
74 import org.orekit.utils.PVCoordinates;
75 import org.orekit.utils.PVCoordinatesProvider;
76
77
78 public class SolarRadiationPressureTest extends AbstractLegacyForceModelTest {
79
80 @Override
81 protected FieldVector3D<DerivativeStructure> accelerationDerivatives(final ForceModel forceModel,
82 final AbsoluteDate date, final Frame frame,
83 final FieldVector3D<DerivativeStructure> position,
84 final FieldVector3D<DerivativeStructure> velocity,
85 final FieldRotation<DerivativeStructure> rotation,
86 final DerivativeStructure mass)
87 {
88 try {
89 java.lang.reflect.Field kRefField = SolarRadiationPressure.class.getDeclaredField("kRef");
90 kRefField.setAccessible(true);
91 double kRef = kRefField.getDouble(forceModel);
92 java.lang.reflect.Field sunField = SolarRadiationPressure.class.getDeclaredField("sun");
93 sunField.setAccessible(true);
94 PVCoordinatesProvider sun = (PVCoordinatesProvider) sunField.get(forceModel);
95 java.lang.reflect.Field spacecraftField = SolarRadiationPressure.class.getDeclaredField("spacecraft");
96 spacecraftField.setAccessible(true);
97 RadiationSensitive spacecraft = (RadiationSensitive) spacecraftField.get(forceModel);
98 java.lang.reflect.Method getLightingRatioMethod = SolarRadiationPressure.class.getDeclaredMethod("getLightingRatio",
99 FieldVector3D.class,
100 Frame.class,
101 FieldAbsoluteDate.class);
102 getLightingRatioMethod.setAccessible(true);
103
104 final Field<DerivativeStructure> field = position.getX().getField();
105 final FieldVector3D<DerivativeStructure> sunSatVector = position.subtract(sun.getPVCoordinates(date, frame).getPosition());
106 final DerivativeStructure r2 = sunSatVector.getNormSq();
107
108
109 final DerivativeStructure ratio = (DerivativeStructure) getLightingRatioMethod.invoke(forceModel, position, frame, new FieldAbsoluteDate<>(field, date));
110 final DerivativeStructure rawP = ratio.multiply(kRef).divide(r2);
111 final FieldVector3D<DerivativeStructure> flux = new FieldVector3D<>(rawP.divide(r2.sqrt()), sunSatVector);
112
113
114 return spacecraft.radiationPressureAcceleration(new FieldAbsoluteDate<>(field, date),
115 frame, position, rotation, mass, flux,
116 forceModel.getParameters(field));
117 } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException |
118 SecurityException | NoSuchMethodException | InvocationTargetException e) {
119 return null;
120 }
121 }
122
123 @Override
124 protected FieldVector3D<Gradient> accelerationDerivativesGradient(final ForceModel forceModel,
125 final AbsoluteDate date, final Frame frame,
126 final FieldVector3D<Gradient> position,
127 final FieldVector3D<Gradient> velocity,
128 final FieldRotation<Gradient> rotation,
129 final Gradient mass)
130 {
131 try {
132 java.lang.reflect.Field kRefField = SolarRadiationPressure.class.getDeclaredField("kRef");
133 kRefField.setAccessible(true);
134 double kRef = kRefField.getDouble(forceModel);
135 java.lang.reflect.Field sunField = SolarRadiationPressure.class.getDeclaredField("sun");
136 sunField.setAccessible(true);
137 PVCoordinatesProvider sun = (PVCoordinatesProvider) sunField.get(forceModel);
138 java.lang.reflect.Field spacecraftField = SolarRadiationPressure.class.getDeclaredField("spacecraft");
139 spacecraftField.setAccessible(true);
140 RadiationSensitive spacecraft = (RadiationSensitive) spacecraftField.get(forceModel);
141 java.lang.reflect.Method getLightingRatioMethod = SolarRadiationPressure.class.getDeclaredMethod("getLightingRatio",
142 FieldVector3D.class,
143 Frame.class,
144 FieldAbsoluteDate.class);
145 getLightingRatioMethod.setAccessible(true);
146
147 final Field<Gradient> field = position.getX().getField();
148 final FieldVector3D<Gradient> sunSatVector = position.subtract(sun.getPVCoordinates(date, frame).getPosition());
149 final Gradient r2 = sunSatVector.getNormSq();
150
151
152 final Gradient ratio = (Gradient) getLightingRatioMethod.invoke(forceModel, position, frame, new FieldAbsoluteDate<>(field, date));
153 final Gradient rawP = ratio.multiply(kRef).divide(r2);
154 final FieldVector3D<Gradient> flux = new FieldVector3D<>(rawP.divide(r2.sqrt()), sunSatVector);
155
156
157 return spacecraft.radiationPressureAcceleration(new FieldAbsoluteDate<>(field, date),
158 frame, position, rotation, mass, flux,
159 forceModel.getParameters(field));
160 } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException |
161 SecurityException | NoSuchMethodException | InvocationTargetException e) {
162 return null;
163 }
164 }
165
166 @Test
167 public void testLightingInterplanetary() throws ParseException {
168
169 AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 3, 21),
170 new TimeComponents(13, 59, 27.816),
171 TimeScalesFactory.getUTC());
172 Orbit orbit = new KeplerianOrbit(1.0e11, 0.1, 0.2, 0.3, 0.4, 0.5, PositionAngle.TRUE,
173 CelestialBodyFactory.getSolarSystemBarycenter().getInertiallyOrientedFrame(),
174 date, Constants.JPL_SSD_SUN_GM);
175 ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
176 SolarRadiationPressure srp =
177 new SolarRadiationPressure(sun, Constants.SUN_RADIUS,
178 (RadiationSensitive) new IsotropicRadiationClassicalConvention(50.0, 0.5, 0.5));
179 Assert.assertFalse(srp.dependsOnPositionOnly());
180
181 Vector3D position = orbit.getPVCoordinates().getPosition();
182 Frame frame = orbit.getFrame();
183 Assert.assertEquals(1.0,
184 srp.getLightingRatio(position, frame, date),
185 1.0e-15);
186
187 Assert.assertEquals(1.0,
188 srp.getLightingRatio(new FieldVector3D<>(Decimal64Field.getInstance(), position),
189 frame,
190 new FieldAbsoluteDate<>(Decimal64Field.getInstance(), date)).getReal(),
191 1.0e-15);
192 }
193
194 @Test
195 public void testLighting() throws ParseException {
196
197 AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 3, 21),
198 new TimeComponents(13, 59, 27.816),
199 TimeScalesFactory.getUTC());
200 Orbit orbit = new EquinoctialOrbit(42164000, 10e-3, 10e-3,
201 FastMath.tan(0.001745329)*FastMath.cos(2*FastMath.PI/3), FastMath.tan(0.001745329)*FastMath.sin(2*FastMath.PI/3),
202 0.1, PositionAngle.TRUE, FramesFactory.getEME2000(), date, mu);
203 ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
204 OneAxisEllipsoid earth =
205 new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765,
206 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
207 SolarRadiationPressure SRP =
208 new SolarRadiationPressure(sun, earth.getEquatorialRadius(),
209 (RadiationSensitive) new IsotropicRadiationCNES95Convention(50.0, 0.5, 0.5));
210
211 double period = 2*FastMath.PI*FastMath.sqrt(orbit.getA()*orbit.getA()*orbit.getA()/orbit.getMu());
212 Assert.assertEquals(86164, period, 1);
213
214
215 KeplerianPropagator k = new KeplerianPropagator(orbit);
216
217
218 AbsoluteDate currentDate;
219 double changed = 1;
220 int count=0;
221
222 for (int t=1;t<3*period;t+=1000) {
223 currentDate = date.shiftedBy(t);
224 try {
225
226 double ratio = SRP.getLightingRatio(k.propagate(currentDate).getPVCoordinates().getPosition(),
227 FramesFactory.getEME2000(), currentDate);
228
229 if (FastMath.floor(ratio)!=changed) {
230 changed = FastMath.floor(ratio);
231 if (changed == 0) {
232 count++;
233 }
234 }
235 } catch (OrekitException e) {
236 e.printStackTrace();
237 }
238 }
239 Assert.assertTrue(3==count);
240 }
241
242 @Test
243 public void testGlobalStateJacobianIsotropicSingle()
244 {
245
246
247 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
248 new TimeComponents(13, 59, 27.816),
249 TimeScalesFactory.getUTC());
250 double i = FastMath.toRadians(98.7);
251 double omega = FastMath.toRadians(93.0);
252 double OMEGA = FastMath.toRadians(15.0 * 22.5);
253 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
254 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
255 Constants.EIGEN5C_EARTH_MU);
256 OrbitType integrationType = OrbitType.CARTESIAN;
257 double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
258
259 NumericalPropagator propagator =
260 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
261 tolerances[0], tolerances[1]));
262 propagator.setOrbitType(integrationType);
263 SolarRadiationPressure forceModel =
264 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
265 new IsotropicRadiationSingleCoefficient(2.5, 0.7));
266 propagator.addForceModel(forceModel);
267 SpacecraftState state0 = new SpacecraftState(orbit);
268
269 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
270 1e3, tolerances[0], 2.0e-5);
271
272 }
273
274 @Test
275 public void testLocalJacobianIsotropicClassicalVs80Implementation()
276 {
277
278
279 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
280 new TimeComponents(13, 59, 27.816),
281 TimeScalesFactory.getUTC());
282 double i = FastMath.toRadians(98.7);
283 double omega = FastMath.toRadians(93.0);
284 double OMEGA = FastMath.toRadians(15.0 * 22.5);
285 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
286 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
287 Constants.EIGEN5C_EARTH_MU);
288 final SolarRadiationPressure forceModel =
289 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
290 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
291
292 checkStateJacobianVs80Implementation(new SpacecraftState(orbit), forceModel,
293 new LofOffset(orbit.getFrame(), LOFType.VVLH),
294 1.0e-15, false);
295
296 }
297
298 @Test
299 public void testLocalJacobianIsotropicClassicalVs80ImplementationGradient()
300 {
301
302
303 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
304 new TimeComponents(13, 59, 27.816),
305 TimeScalesFactory.getUTC());
306 double i = FastMath.toRadians(98.7);
307 double omega = FastMath.toRadians(93.0);
308 double OMEGA = FastMath.toRadians(15.0 * 22.5);
309 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
310 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
311 Constants.EIGEN5C_EARTH_MU);
312 final SolarRadiationPressure forceModel =
313 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
314 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
315
316 checkStateJacobianVs80ImplementationGradient(new SpacecraftState(orbit), forceModel,
317 new LofOffset(orbit.getFrame(), LOFType.VVLH),
318 1.0e-15, false);
319
320 }
321
322 @Test
323 public void testLocalJacobianIsotropicClassicalVsFiniteDifferencesFullLight()
324 {
325
326 doTestLocalJacobianIsotropicClassicalVsFiniteDifferences(250.0, 1000.0, 3.0e-8, false);
327 }
328
329 @Test
330 public void testLocalJacobianIsotropicClassicalVsFiniteDifferencesGradientFullLight()
331 {
332
333 doTestLocalJacobianIsotropicClassicalVsFiniteDifferencesGradient(250.0, 1000.0, 3.0e-8, false);
334 }
335
336 @Test
337 public void testLocalJacobianIsotropicClassicalVsFiniteDifferencesPenumbra()
338 {
339
340
341 doTestLocalJacobianIsotropicClassicalVsFiniteDifferences(275.5, 100.0, 8.0e-7, false);
342 }
343
344 @Test
345 public void testLocalJacobianIsotropicClassicalVsFiniteDifferencesGradientPenumbra()
346 {
347
348
349 doTestLocalJacobianIsotropicClassicalVsFiniteDifferencesGradient(275.5, 100.0, 8.0e-7, false);
350 }
351
352 @Test
353 public void testLocalJacobianIsotropicClassicalVsFiniteDifferencesEclipse()
354 {
355
356 doTestLocalJacobianIsotropicClassicalVsFiniteDifferences(300.0, 1000.0, 1.0e-50, false);
357 }
358
359 @Test
360 public void testLocalJacobianIsotropicClassicalVsFiniteDifferencesGradientEclipse()
361 {
362
363 doTestLocalJacobianIsotropicClassicalVsFiniteDifferencesGradient(300.0, 1000.0, 1.0e-50, false);
364 }
365
366 private void doTestLocalJacobianIsotropicClassicalVsFiniteDifferences(double deltaT, double dP,
367 double checkTolerance,
368 boolean print)
369 {
370
371
372 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
373 new TimeComponents(13, 59, 27.816),
374 TimeScalesFactory.getUTC());
375 double i = FastMath.toRadians(98.7);
376 double omega = FastMath.toRadians(93.0);
377 double OMEGA = FastMath.toRadians(15.0 * 22.5);
378 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
379 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
380 Constants.EIGEN5C_EARTH_MU);
381 final SolarRadiationPressure forceModel =
382 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
383 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
384
385 checkStateJacobianVsFiniteDifferences(new SpacecraftState(orbit.shiftedBy(deltaT)), forceModel,
386 Propagator.DEFAULT_LAW, dP, checkTolerance, print);
387
388 }
389
390 private void doTestLocalJacobianIsotropicClassicalVsFiniteDifferencesGradient(double deltaT, double dP,
391 double checkTolerance,
392 boolean print)
393 {
394
395
396 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
397 new TimeComponents(13, 59, 27.816),
398 TimeScalesFactory.getUTC());
399 double i = FastMath.toRadians(98.7);
400 double omega = FastMath.toRadians(93.0);
401 double OMEGA = FastMath.toRadians(15.0 * 22.5);
402 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
403 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
404 Constants.EIGEN5C_EARTH_MU);
405 final SolarRadiationPressure forceModel =
406 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
407 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
408
409 checkStateJacobianVsFiniteDifferencesGradient(new SpacecraftState(orbit.shiftedBy(deltaT)), forceModel,
410 Propagator.DEFAULT_LAW, dP, checkTolerance, print);
411
412 }
413
414 @Test
415 public void testGlobalStateJacobianIsotropicClassical()
416 {
417
418
419 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
420 new TimeComponents(13, 59, 27.816),
421 TimeScalesFactory.getUTC());
422 double i = FastMath.toRadians(98.7);
423 double omega = FastMath.toRadians(93.0);
424 double OMEGA = FastMath.toRadians(15.0 * 22.5);
425 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
426 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
427 Constants.EIGEN5C_EARTH_MU);
428 OrbitType integrationType = OrbitType.CARTESIAN;
429 double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
430
431 NumericalPropagator propagator =
432 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
433 tolerances[0], tolerances[1]));
434 propagator.setOrbitType(integrationType);
435 SolarRadiationPressure forceModel =
436 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
437 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
438 propagator.addForceModel(forceModel);
439 SpacecraftState state0 = new SpacecraftState(orbit);
440
441 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
442 1e6, tolerances[0], 2.0e-5);
443
444 }
445
446 @Test
447 public void testGlobalStateJacobianIsotropicCnes()
448 {
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 OrbitType integrationType = OrbitType.CARTESIAN;
461 double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
462
463 NumericalPropagator propagator =
464 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
465 tolerances[0], tolerances[1]));
466 propagator.setOrbitType(integrationType);
467 SolarRadiationPressure forceModel =
468 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
469 new IsotropicRadiationCNES95Convention(2.5, 0.7, 0.2));
470 propagator.addForceModel(forceModel);
471 SpacecraftState state0 = new SpacecraftState(orbit);
472
473 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
474 1e3, tolerances[0], 3.0e-5);
475
476 }
477
478 @Test
479 public void testParameterDerivativeBox() {
480
481 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
482 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
483 final SpacecraftState state =
484 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
485 FramesFactory.getGCRF(),
486 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
487 Constants.EIGEN5C_EARTH_MU));
488
489 SolarRadiationPressure forceModel =
490 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
491 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
492 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
493
494 checkParameterDerivative(state, forceModel, RadiationSensitive.ABSORPTION_COEFFICIENT, 0.25, 1.9e-15);
495 checkParameterDerivative(state, forceModel, RadiationSensitive.REFLECTION_COEFFICIENT, 0.25, 6.9e-15);
496
497 }
498
499 @Test
500 public void testParameterDerivativeGradientBox() {
501
502 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
503 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
504 final SpacecraftState state =
505 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
506 FramesFactory.getGCRF(),
507 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
508 Constants.EIGEN5C_EARTH_MU));
509
510 SolarRadiationPressure forceModel =
511 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
512 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
513 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
514
515 checkParameterDerivativeGradient(state, forceModel, RadiationSensitive.ABSORPTION_COEFFICIENT, 0.25, 1.9e-15);
516 checkParameterDerivativeGradient(state, forceModel, RadiationSensitive.REFLECTION_COEFFICIENT, 0.25, 6.9e-10);
517
518 }
519
520 @Test
521 public void testGlobalStateJacobianBox()
522 {
523
524
525 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
526 new TimeComponents(13, 59, 27.816),
527 TimeScalesFactory.getUTC());
528 double i = FastMath.toRadians(98.7);
529 double omega = FastMath.toRadians(93.0);
530 double OMEGA = FastMath.toRadians(15.0 * 22.5);
531 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
532 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
533 Constants.EIGEN5C_EARTH_MU);
534 OrbitType integrationType = OrbitType.CARTESIAN;
535 double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
536
537 NumericalPropagator propagator =
538 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
539 tolerances[0], tolerances[1]));
540 propagator.setOrbitType(integrationType);
541 SolarRadiationPressure forceModel =
542 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
543 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
544 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
545 propagator.addForceModel(forceModel);
546 SpacecraftState state0 = new SpacecraftState(orbit);
547
548 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
549 1e3, tolerances[0], 5.0e-4);
550
551 }
552
553 @Test
554 public void testRoughOrbitalModifs() throws ParseException, OrekitException, FileNotFoundException {
555
556
557 AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 7, 1),
558 new TimeComponents(13, 59, 27.816),
559 TimeScalesFactory.getUTC());
560 Orbit orbit = new EquinoctialOrbit(42164000, 10e-3, 10e-3,
561 FastMath.tan(0.001745329)*FastMath.cos(2*FastMath.PI/3),
562 FastMath.tan(0.001745329)*FastMath.sin(2*FastMath.PI/3),
563 0.1, PositionAngle.TRUE, FramesFactory.getEME2000(), date, mu);
564 final double period = orbit.getKeplerianPeriod();
565 Assert.assertEquals(86164, period, 1);
566 ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
567
568
569 OneAxisEllipsoid earth =
570 new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765,
571 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
572 SolarRadiationPressure SRP =
573 new SolarRadiationPressure(sun, earth.getEquatorialRadius(),
574 new IsotropicRadiationCNES95Convention(500.0, 0.7, 0.7));
575
576
577 double[] absTolerance = {
578 0.1, 1.0e-9, 1.0e-9, 1.0e-5, 1.0e-5, 1.0e-5, 0.001
579 };
580 double[] relTolerance = {
581 1.0e-4, 1.0e-4, 1.0e-4, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-7
582 };
583 AdaptiveStepsizeIntegrator integrator =
584 new DormandPrince853Integrator(900.0, 60000, absTolerance, relTolerance);
585 integrator.setInitialStepSize(3600);
586 final NumericalPropagator calc = new NumericalPropagator(integrator);
587 calc.addForceModel(SRP);
588
589
590 calc.setMasterMode(FastMath.floor(period), new SolarStepHandler());
591 AbsoluteDate finalDate = date.shiftedBy(10 * period);
592 calc.setInitialState(new SpacecraftState(orbit, 1500.0));
593 calc.propagate(finalDate);
594 Assert.assertTrue(calc.getCalls() < 7100);
595 }
596
597 public static void checkRadius(double radius , double min , double max) {
598 Assert.assertTrue(radius >= min);
599 Assert.assertTrue(radius <= max);
600 }
601
602 private double mu = 3.98600E14;
603
604 private static class SolarStepHandler implements OrekitFixedStepHandler {
605
606 public void handleStep(SpacecraftState currentState, boolean isLast) {
607 final double dex = currentState.getEquinoctialEx() - 0.01071166;
608 final double dey = currentState.getEquinoctialEy() - 0.00654848;
609 final double alpha = FastMath.toDegrees(FastMath.atan2(dey, dex));
610 Assert.assertTrue(alpha > 100.0);
611 Assert.assertTrue(alpha < 112.0);
612 checkRadius(FastMath.sqrt(dex * dex + dey * dey), 0.003524, 0.003541);
613 }
614
615 }
616
617
618
619
620
621
622 @Test
623 public void RealFieldIsotropicTest() {
624
625
626 DSFactory factory = new DSFactory(6, 5);
627 DerivativeStructure a_0 = factory.variable(0, 7e7);
628 DerivativeStructure e_0 = factory.variable(1, 0.4);
629 DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
630 DerivativeStructure R_0 = factory.variable(3, 0.7);
631 DerivativeStructure O_0 = factory.variable(4, 0.5);
632 DerivativeStructure n_0 = factory.variable(5, 0.1);
633
634 Field<DerivativeStructure> field = a_0.getField();
635 DerivativeStructure zero = field.getZero();
636
637
638 FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
639
640
641 Frame EME = FramesFactory.getEME2000();
642
643
644 FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
645 PositionAngle.MEAN,
646 EME,
647 J2000,
648 zero.add(Constants.EIGEN5C_EARTH_MU));
649
650 FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
651 SpacecraftState iSR = initialState.toSpacecraftState();
652
653
654 final OrbitType type = OrbitType.KEPLERIAN;
655 double[][] tolerance = NumericalPropagator.tolerances(10.0, FKO.toOrbit(), type);
656 AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
657 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
658 integrator.setInitialStepSize(zero.add(60));
659
660 AdaptiveStepsizeIntegrator RIntegrator =
661 new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
662 RIntegrator.setInitialStepSize(60);
663
664
665 FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
666 FNP.setOrbitType(type);
667 FNP.setInitialState(initialState);
668
669 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
670 NP.setOrbitType(type);
671 NP.setInitialState(iSR);
672
673 ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
674
675
676 OneAxisEllipsoid earth =
677 new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765,
678 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
679 SolarRadiationPressure forceModel =
680 new SolarRadiationPressure(sun, earth.getEquatorialRadius(),
681 new IsotropicRadiationCNES95Convention(500.0, 0.7, 0.7));
682
683 FNP.addForceModel(forceModel);
684 NP.addForceModel(forceModel);
685
686
687 checkRealFieldPropagation(FKO, PositionAngle.MEAN, 1000., NP, FNP,
688 1.0e-30, 5.0e-10, 3.0e-11, 3.0e-10,
689 1, false);
690 }
691
692
693
694
695
696 @Test
697 public void RealFieldExpectErrorTest() {
698 DSFactory factory = new DSFactory(6, 0);
699 DerivativeStructure a_0 = factory.variable(0, 7e7);
700 DerivativeStructure e_0 = factory.variable(1, 0.4);
701 DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
702 DerivativeStructure R_0 = factory.variable(3, 0.7);
703 DerivativeStructure O_0 = factory.variable(4, 0.5);
704 DerivativeStructure n_0 = factory.variable(5, 0.1);
705
706 Field<DerivativeStructure> field = a_0.getField();
707 DerivativeStructure zero = field.getZero();
708
709 FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
710
711 Frame EME = FramesFactory.getEME2000();
712
713 FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
714 PositionAngle.MEAN,
715 EME,
716 J2000,
717 zero.add(Constants.EIGEN5C_EARTH_MU));
718
719 FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
720
721 SpacecraftState iSR = initialState.toSpacecraftState();
722
723 final OrbitType type = OrbitType.KEPLERIAN;
724 double[][] tolerance = NumericalPropagator.tolerances(0.001, FKO.toOrbit(), type);
725
726
727 AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
728 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
729 integrator.setInitialStepSize(zero.add(60));
730 AdaptiveStepsizeIntegrator RIntegrator =
731 new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
732 RIntegrator.setInitialStepSize(60);
733
734 FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
735 FNP.setOrbitType(type);
736 FNP.setInitialState(initialState);
737
738 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
739 NP.setInitialState(iSR);
740 ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
741
742
743 OneAxisEllipsoid earth =
744 new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765,
745 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
746 SolarRadiationPressure forceModel =
747 new SolarRadiationPressure(sun, earth.getEquatorialRadius(),
748 new IsotropicRadiationCNES95Convention(500.0, 0.7, 0.7));
749
750 FNP.addForceModel(forceModel);
751
752
753 FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1000.);
754 FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
755 SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
756 FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
757 PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
758
759 Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getX() - finPVC_R.getPosition().getX()) < FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
760 Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getY() - finPVC_R.getPosition().getY()) < FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
761 Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getZ() - finPVC_R.getPosition().getZ()) < FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
762 }
763
764 @Before
765 public void setUp() {
766 Utils.setDataRoot("regular-data");
767 }
768
769 }