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.CalculusFieldElement;
25 import org.hipparchus.Field;
26 import org.hipparchus.analysis.differentiation.DSFactory;
27 import org.hipparchus.analysis.differentiation.DerivativeStructure;
28 import org.hipparchus.analysis.differentiation.Gradient;
29 import org.hipparchus.geometry.euclidean.threed.FieldRotation;
30 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
31 import org.hipparchus.geometry.euclidean.threed.Vector3D;
32 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
33 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
34 import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
35 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
36 import org.hipparchus.util.Decimal64;
37 import org.hipparchus.util.Decimal64Field;
38 import org.hipparchus.util.FastMath;
39 import org.junit.Assert;
40 import org.junit.Before;
41 import org.junit.Test;
42 import org.orekit.Utils;
43 import org.orekit.attitudes.LofOffset;
44 import org.orekit.bodies.CelestialBodyFactory;
45 import org.orekit.bodies.OneAxisEllipsoid;
46 import org.orekit.errors.OrekitException;
47 import org.orekit.forces.AbstractLegacyForceModelTest;
48 import org.orekit.forces.BoxAndSolarArraySpacecraft;
49 import org.orekit.forces.ForceModel;
50 import org.orekit.frames.Frame;
51 import org.orekit.frames.FramesFactory;
52 import org.orekit.frames.LOFType;
53 import org.orekit.orbits.CartesianOrbit;
54 import org.orekit.orbits.EquinoctialOrbit;
55 import org.orekit.orbits.FieldCartesianOrbit;
56 import org.orekit.orbits.FieldKeplerianOrbit;
57 import org.orekit.orbits.FieldOrbit;
58 import org.orekit.orbits.KeplerianOrbit;
59 import org.orekit.orbits.Orbit;
60 import org.orekit.orbits.OrbitType;
61 import org.orekit.orbits.PositionAngle;
62 import org.orekit.propagation.FieldSpacecraftState;
63 import org.orekit.propagation.SpacecraftState;
64 import org.orekit.propagation.analytical.KeplerianPropagator;
65 import org.orekit.propagation.numerical.FieldNumericalPropagator;
66 import org.orekit.propagation.numerical.NumericalPropagator;
67 import org.orekit.propagation.sampling.FieldOrekitFixedStepHandler;
68 import org.orekit.propagation.sampling.OrekitFixedStepHandler;
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.ExtendedPVCoordinatesProvider;
76 import org.orekit.utils.FieldPVCoordinates;
77 import org.orekit.utils.IERSConventions;
78 import org.orekit.utils.PVCoordinates;
79 import org.orekit.utils.PVCoordinatesProvider;
80 import org.orekit.utils.TimeStampedFieldPVCoordinates;
81 import org.orekit.utils.TimeStampedPVCoordinates;
82
83
84 public class SolarRadiationPressureTest extends AbstractLegacyForceModelTest {
85
86 @Override
87 protected FieldVector3D<DerivativeStructure> accelerationDerivatives(final ForceModel forceModel,
88 final AbsoluteDate date, final Frame frame,
89 final FieldVector3D<DerivativeStructure> position,
90 final FieldVector3D<DerivativeStructure> velocity,
91 final FieldRotation<DerivativeStructure> rotation,
92 final DerivativeStructure mass)
93 {
94 try {
95 java.lang.reflect.Field kRefField = SolarRadiationPressure.class.getDeclaredField("kRef");
96 kRefField.setAccessible(true);
97 double kRef = kRefField.getDouble(forceModel);
98 java.lang.reflect.Field sunField = SolarRadiationPressure.class.getDeclaredField("sun");
99 sunField.setAccessible(true);
100 PVCoordinatesProvider sun = (PVCoordinatesProvider) sunField.get(forceModel);
101 java.lang.reflect.Field spacecraftField = SolarRadiationPressure.class.getDeclaredField("spacecraft");
102 spacecraftField.setAccessible(true);
103 RadiationSensitive spacecraft = (RadiationSensitive) spacecraftField.get(forceModel);
104 java.lang.reflect.Method getLightingRatioMethod = SolarRadiationPressure.class.getDeclaredMethod("getLightingRatio",
105 FieldVector3D.class,
106 Frame.class,
107 FieldAbsoluteDate.class);
108 getLightingRatioMethod.setAccessible(true);
109
110 final Field<DerivativeStructure> field = position.getX().getField();
111 final FieldVector3D<DerivativeStructure> sunSatVector = position.subtract(sun.getPVCoordinates(date, frame).getPosition());
112 final DerivativeStructure r2 = sunSatVector.getNormSq();
113
114
115 final DerivativeStructure ratio = (DerivativeStructure) getLightingRatioMethod.invoke(forceModel, position, frame, new FieldAbsoluteDate<>(field, date));
116 final DerivativeStructure rawP = ratio.multiply(kRef).divide(r2);
117 final FieldVector3D<DerivativeStructure> flux = new FieldVector3D<>(rawP.divide(r2.sqrt()), sunSatVector);
118
119
120 return spacecraft.radiationPressureAcceleration(new FieldAbsoluteDate<>(field, date),
121 frame, position, rotation, mass, flux,
122 forceModel.getParameters(field));
123 } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException |
124 SecurityException | NoSuchMethodException | InvocationTargetException e) {
125 return null;
126 }
127 }
128
129 @Override
130 protected FieldVector3D<Gradient> accelerationDerivativesGradient(final ForceModel forceModel,
131 final AbsoluteDate date, final Frame frame,
132 final FieldVector3D<Gradient> position,
133 final FieldVector3D<Gradient> velocity,
134 final FieldRotation<Gradient> rotation,
135 final Gradient mass)
136 {
137 try {
138 java.lang.reflect.Field kRefField = SolarRadiationPressure.class.getDeclaredField("kRef");
139 kRefField.setAccessible(true);
140 double kRef = kRefField.getDouble(forceModel);
141 java.lang.reflect.Field sunField = SolarRadiationPressure.class.getDeclaredField("sun");
142 sunField.setAccessible(true);
143 PVCoordinatesProvider sun = (PVCoordinatesProvider) sunField.get(forceModel);
144 java.lang.reflect.Field spacecraftField = SolarRadiationPressure.class.getDeclaredField("spacecraft");
145 spacecraftField.setAccessible(true);
146 RadiationSensitive spacecraft = (RadiationSensitive) spacecraftField.get(forceModel);
147 java.lang.reflect.Method getLightingRatioMethod = SolarRadiationPressure.class.getDeclaredMethod("getLightingRatio",
148 FieldVector3D.class,
149 Frame.class,
150 FieldAbsoluteDate.class);
151 getLightingRatioMethod.setAccessible(true);
152
153 final Field<Gradient> field = position.getX().getField();
154 final FieldVector3D<Gradient> sunSatVector = position.subtract(sun.getPVCoordinates(date, frame).getPosition());
155 final Gradient r2 = sunSatVector.getNormSq();
156
157
158 final Gradient ratio = (Gradient) getLightingRatioMethod.invoke(forceModel, position, frame, new FieldAbsoluteDate<>(field, date));
159 final Gradient rawP = ratio.multiply(kRef).divide(r2);
160 final FieldVector3D<Gradient> flux = new FieldVector3D<>(rawP.divide(r2.sqrt()), sunSatVector);
161
162
163 return spacecraft.radiationPressureAcceleration(new FieldAbsoluteDate<>(field, date),
164 frame, position, rotation, mass, flux,
165 forceModel.getParameters(field));
166 } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException |
167 SecurityException | NoSuchMethodException | InvocationTargetException e) {
168 return null;
169 }
170 }
171
172 @Test
173 public void testLightingInterplanetary() throws ParseException {
174
175 AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 3, 21),
176 new TimeComponents(13, 59, 27.816),
177 TimeScalesFactory.getUTC());
178 Orbit orbit = new KeplerianOrbit(1.0e11, 0.1, 0.2, 0.3, 0.4, 0.5, PositionAngle.TRUE,
179 CelestialBodyFactory.getSolarSystemBarycenter().getInertiallyOrientedFrame(),
180 date, Constants.JPL_SSD_SUN_GM);
181 ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
182 SolarRadiationPressure srp =
183 new SolarRadiationPressure(sun, Constants.SUN_RADIUS,
184 new IsotropicRadiationClassicalConvention(50.0, 0.5, 0.5));
185 Assert.assertFalse(srp.dependsOnPositionOnly());
186
187 Vector3D position = orbit.getPVCoordinates().getPosition();
188 Frame frame = orbit.getFrame();
189 Assert.assertEquals(1.0,
190 srp.getLightingRatio(position, frame, date),
191 1.0e-15);
192
193 Assert.assertEquals(1.0,
194 srp.getLightingRatio(new FieldVector3D<>(Decimal64Field.getInstance(), position),
195 frame,
196 new FieldAbsoluteDate<>(Decimal64Field.getInstance(), date)).getReal(),
197 1.0e-15);
198 }
199
200 @Test
201 public void testLighting() throws ParseException {
202
203 AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 3, 21),
204 new TimeComponents(13, 59, 27.816),
205 TimeScalesFactory.getUTC());
206 Orbit orbit = new EquinoctialOrbit(42164000, 10e-3, 10e-3,
207 FastMath.tan(0.001745329)*FastMath.cos(2*FastMath.PI/3), FastMath.tan(0.001745329)*FastMath.sin(2*FastMath.PI/3),
208 0.1, PositionAngle.TRUE, FramesFactory.getEME2000(), date, mu);
209 ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
210 OneAxisEllipsoid earth =
211 new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765,
212 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
213 SolarRadiationPressure SRP =
214 new SolarRadiationPressure(sun, earth.getEquatorialRadius(),
215 new IsotropicRadiationCNES95Convention(50.0, 0.5, 0.5));
216
217 double period = 2*FastMath.PI*FastMath.sqrt(orbit.getA()*orbit.getA()*orbit.getA()/orbit.getMu());
218 Assert.assertEquals(86164, period, 1);
219
220
221 KeplerianPropagator k = new KeplerianPropagator(orbit);
222
223
224 AbsoluteDate currentDate;
225 double changed = 1;
226 int count=0;
227
228 for (int t=1;t<3*period;t+=1000) {
229 currentDate = date.shiftedBy(t);
230 try {
231
232 double ratio = SRP.getLightingRatio(k.propagate(currentDate).getPVCoordinates().getPosition(),
233 FramesFactory.getEME2000(), currentDate);
234
235 if (FastMath.floor(ratio)!=changed) {
236 changed = FastMath.floor(ratio);
237 if (changed == 0) {
238 count++;
239 }
240 }
241 } catch (OrekitException e) {
242 e.printStackTrace();
243 }
244 }
245 Assert.assertTrue(3==count);
246 }
247
248 @Test
249 public void testGlobalStateJacobianIsotropicSingle()
250 {
251
252
253 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
254 new TimeComponents(13, 59, 27.816),
255 TimeScalesFactory.getUTC());
256 double i = FastMath.toRadians(98.7);
257 double omega = FastMath.toRadians(93.0);
258 double OMEGA = FastMath.toRadians(15.0 * 22.5);
259 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
260 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
261 Constants.EIGEN5C_EARTH_MU);
262 OrbitType integrationType = OrbitType.CARTESIAN;
263 double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
264
265 NumericalPropagator propagator =
266 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
267 tolerances[0], tolerances[1]));
268 propagator.setOrbitType(integrationType);
269 SolarRadiationPressure forceModel =
270 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
271 new IsotropicRadiationSingleCoefficient(2.5, 0.7));
272 propagator.addForceModel(forceModel);
273 SpacecraftState state0 = new SpacecraftState(orbit);
274
275 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
276 1e3, tolerances[0], 2.0e-5);
277
278 }
279
280 @Test
281 public void testLocalJacobianIsotropicClassicalVs80Implementation()
282 {
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 final SolarRadiationPressure forceModel =
295 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
296 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
297
298 checkStateJacobianVs80Implementation(new SpacecraftState(orbit), forceModel,
299 new LofOffset(orbit.getFrame(), LOFType.LVLH_CCSDS),
300 1.0e-15, false);
301
302 }
303
304 @Test
305 public void testLocalJacobianIsotropicClassicalVs80ImplementationGradient()
306 {
307
308
309 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
310 new TimeComponents(13, 59, 27.816),
311 TimeScalesFactory.getUTC());
312 double i = FastMath.toRadians(98.7);
313 double omega = FastMath.toRadians(93.0);
314 double OMEGA = FastMath.toRadians(15.0 * 22.5);
315 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
316 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
317 Constants.EIGEN5C_EARTH_MU);
318 final SolarRadiationPressure forceModel =
319 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
320 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
321
322 checkStateJacobianVs80ImplementationGradient(new SpacecraftState(orbit), forceModel,
323 new LofOffset(orbit.getFrame(), LOFType.LVLH_CCSDS),
324 1.0e-15, false);
325
326 }
327
328 @Test
329 public void testLocalJacobianIsotropicClassicalVsFiniteDifferencesFullLight()
330 {
331
332 doTestLocalJacobianIsotropicClassicalVsFiniteDifferences(250.0, 1000.0, 3.0e-8, false);
333 }
334
335 @Test
336 public void testLocalJacobianIsotropicClassicalVsFiniteDifferencesGradientFullLight()
337 {
338
339 doTestLocalJacobianIsotropicClassicalVsFiniteDifferencesGradient(250.0, 1000.0, 3.0e-8, false);
340 }
341
342 @Test
343 public void testLocalJacobianIsotropicClassicalVsFiniteDifferencesPenumbra()
344 {
345
346
347 doTestLocalJacobianIsotropicClassicalVsFiniteDifferences(275.5, 100.0, 8.0e-7, false);
348 }
349
350 @Test
351 public void testLocalJacobianIsotropicClassicalVsFiniteDifferencesGradientPenumbra()
352 {
353
354
355 doTestLocalJacobianIsotropicClassicalVsFiniteDifferencesGradient(275.5, 100.0, 8.0e-7, false);
356 }
357
358 @Test
359 public void testLocalJacobianIsotropicClassicalVsFiniteDifferencesEclipse()
360 {
361
362 doTestLocalJacobianIsotropicClassicalVsFiniteDifferences(300.0, 1000.0, 1.0e-50, false);
363 }
364
365 @Test
366 public void testLocalJacobianIsotropicClassicalVsFiniteDifferencesGradientEclipse()
367 {
368
369 doTestLocalJacobianIsotropicClassicalVsFiniteDifferencesGradient(300.0, 1000.0, 1.0e-50, false);
370 }
371
372 private void doTestLocalJacobianIsotropicClassicalVsFiniteDifferences(double deltaT, double dP,
373 double checkTolerance,
374 boolean print)
375 {
376
377
378 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
379 new TimeComponents(13, 59, 27.816),
380 TimeScalesFactory.getUTC());
381 double i = FastMath.toRadians(98.7);
382 double omega = FastMath.toRadians(93.0);
383 double OMEGA = FastMath.toRadians(15.0 * 22.5);
384 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
385 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
386 Constants.EIGEN5C_EARTH_MU);
387 final SolarRadiationPressure forceModel =
388 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
389 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
390
391 checkStateJacobianVsFiniteDifferences(new SpacecraftState(orbit.shiftedBy(deltaT)), forceModel,
392 Utils.defaultLaw(), dP, checkTolerance, print);
393
394 }
395
396 private void doTestLocalJacobianIsotropicClassicalVsFiniteDifferencesGradient(double deltaT, double dP,
397 double checkTolerance,
398 boolean print)
399 {
400
401
402 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
403 new TimeComponents(13, 59, 27.816),
404 TimeScalesFactory.getUTC());
405 double i = FastMath.toRadians(98.7);
406 double omega = FastMath.toRadians(93.0);
407 double OMEGA = FastMath.toRadians(15.0 * 22.5);
408 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
409 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
410 Constants.EIGEN5C_EARTH_MU);
411 final SolarRadiationPressure forceModel =
412 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
413 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
414
415 checkStateJacobianVsFiniteDifferencesGradient(new SpacecraftState(orbit.shiftedBy(deltaT)), forceModel,
416 Utils.defaultLaw(), dP, checkTolerance, print);
417
418 }
419
420 @Test
421 public void testGlobalStateJacobianIsotropicClassical()
422 {
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 OrbitType integrationType = OrbitType.CARTESIAN;
435 double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
436
437 NumericalPropagator propagator =
438 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
439 tolerances[0], tolerances[1]));
440 propagator.setOrbitType(integrationType);
441 SolarRadiationPressure forceModel =
442 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
443 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
444 propagator.addForceModel(forceModel);
445 SpacecraftState state0 = new SpacecraftState(orbit);
446
447 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
448 1e6, tolerances[0], 2.0e-5);
449
450 }
451
452 @Test
453 public void testGlobalStateJacobianIsotropicCnes()
454 {
455
456
457 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
458 new TimeComponents(13, 59, 27.816),
459 TimeScalesFactory.getUTC());
460 double i = FastMath.toRadians(98.7);
461 double omega = FastMath.toRadians(93.0);
462 double OMEGA = FastMath.toRadians(15.0 * 22.5);
463 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
464 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
465 Constants.EIGEN5C_EARTH_MU);
466 OrbitType integrationType = OrbitType.CARTESIAN;
467 double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
468
469 NumericalPropagator propagator =
470 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
471 tolerances[0], tolerances[1]));
472 propagator.setOrbitType(integrationType);
473 SolarRadiationPressure forceModel =
474 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
475 new IsotropicRadiationCNES95Convention(2.5, 0.7, 0.2));
476 propagator.addForceModel(forceModel);
477 SpacecraftState state0 = new SpacecraftState(orbit);
478
479 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
480 1e3, tolerances[0], 3.0e-5);
481
482 }
483
484 @Test
485 public void testParameterDerivativeBox() {
486
487 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
488 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
489 final SpacecraftState state =
490 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
491 FramesFactory.getGCRF(),
492 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
493 Constants.EIGEN5C_EARTH_MU));
494
495 SolarRadiationPressure forceModel =
496 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
497 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
498 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
499
500 checkParameterDerivative(state, forceModel, RadiationSensitive.ABSORPTION_COEFFICIENT, 0.25, 1.9e-15);
501 checkParameterDerivative(state, forceModel, RadiationSensitive.REFLECTION_COEFFICIENT, 0.25, 6.9e-15);
502
503 }
504
505 @Test
506 public void testParameterDerivativeGradientBox() {
507
508 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
509 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
510 final SpacecraftState state =
511 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
512 FramesFactory.getGCRF(),
513 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
514 Constants.EIGEN5C_EARTH_MU));
515
516 SolarRadiationPressure forceModel =
517 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
518 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
519 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
520
521 checkParameterDerivativeGradient(state, forceModel, RadiationSensitive.ABSORPTION_COEFFICIENT, 0.25, 1.9e-15);
522 checkParameterDerivativeGradient(state, forceModel, RadiationSensitive.REFLECTION_COEFFICIENT, 0.25, 6.9e-10);
523
524 }
525
526 @Test
527 public void testGlobalStateJacobianBox()
528 {
529
530
531 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
532 new TimeComponents(13, 59, 27.816),
533 TimeScalesFactory.getUTC());
534 double i = FastMath.toRadians(98.7);
535 double omega = FastMath.toRadians(93.0);
536 double OMEGA = FastMath.toRadians(15.0 * 22.5);
537 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
538 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date,
539 Constants.EIGEN5C_EARTH_MU);
540 OrbitType integrationType = OrbitType.CARTESIAN;
541 double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType);
542
543 NumericalPropagator propagator =
544 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
545 tolerances[0], tolerances[1]));
546 propagator.setOrbitType(integrationType);
547 SolarRadiationPressure forceModel =
548 new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
549 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
550 Vector3D.PLUS_J, 1.2, 0.7, 0.2));
551 propagator.addForceModel(forceModel);
552 SpacecraftState state0 = new SpacecraftState(orbit);
553
554 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
555 1e3, tolerances[0], 5.0e-4);
556
557 }
558
559 @Test
560 public void testRoughOrbitalModifs() throws ParseException, OrekitException, FileNotFoundException {
561
562
563 AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 7, 1),
564 new TimeComponents(13, 59, 27.816),
565 TimeScalesFactory.getUTC());
566 Orbit orbit = new EquinoctialOrbit(42164000, 10e-3, 10e-3,
567 FastMath.tan(0.001745329)*FastMath.cos(2*FastMath.PI/3),
568 FastMath.tan(0.001745329)*FastMath.sin(2*FastMath.PI/3),
569 0.1, PositionAngle.TRUE, FramesFactory.getEME2000(), date, mu);
570 final double period = orbit.getKeplerianPeriod();
571 Assert.assertEquals(86164, period, 1);
572 ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
573
574
575 OneAxisEllipsoid earth =
576 new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765,
577 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
578 SolarRadiationPressure SRP =
579 new SolarRadiationPressure(sun, earth.getEquatorialRadius(),
580 new IsotropicRadiationCNES95Convention(500.0, 0.7, 0.7));
581
582
583 double[] absTolerance = {
584 0.1, 1.0e-9, 1.0e-9, 1.0e-5, 1.0e-5, 1.0e-5, 0.001
585 };
586 double[] relTolerance = {
587 1.0e-4, 1.0e-4, 1.0e-4, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-7
588 };
589 AdaptiveStepsizeIntegrator integrator =
590 new DormandPrince853Integrator(900.0, 60000, absTolerance, relTolerance);
591 integrator.setInitialStepSize(3600);
592 final NumericalPropagator calc = new NumericalPropagator(integrator);
593 calc.addForceModel(SRP);
594
595
596 calc.setStepHandler(FastMath.floor(period), new SolarStepHandler());
597 AbsoluteDate finalDate = date.shiftedBy(10 * period);
598 calc.setInitialState(new SpacecraftState(orbit, 1500.0));
599 calc.propagate(finalDate);
600 Assert.assertTrue(calc.getCalls() < 7100);
601 }
602
603 public static void checkRadius(double radius , double min , double max) {
604 Assert.assertTrue(radius >= min);
605 Assert.assertTrue(radius <= max);
606 }
607
608 private double mu = 3.98600E14;
609
610 private static class SolarStepHandler implements OrekitFixedStepHandler {
611
612 @Override
613 public void handleStep(SpacecraftState currentState) {
614 final double dex = currentState.getEquinoctialEx() - 0.01071166;
615 final double dey = currentState.getEquinoctialEy() - 0.00654848;
616 final double alpha = FastMath.toDegrees(FastMath.atan2(dey, dex));
617 Assert.assertTrue(alpha > 100.0);
618 Assert.assertTrue(alpha < 112.0);
619 checkRadius(FastMath.sqrt(dex * dex + dey * dey), 0.003524, 0.003541);
620 }
621
622 }
623
624
625
626
627
628
629 @Test
630 public void RealFieldIsotropicTest() {
631
632
633 DSFactory factory = new DSFactory(6, 5);
634 DerivativeStructure a_0 = factory.variable(0, 7e7);
635 DerivativeStructure e_0 = factory.variable(1, 0.4);
636 DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
637 DerivativeStructure R_0 = factory.variable(3, 0.7);
638 DerivativeStructure O_0 = factory.variable(4, 0.5);
639 DerivativeStructure n_0 = factory.variable(5, 0.1);
640
641 Field<DerivativeStructure> field = a_0.getField();
642 DerivativeStructure zero = field.getZero();
643
644
645 FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
646
647
648 Frame EME = FramesFactory.getEME2000();
649
650
651 FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
652 PositionAngle.MEAN,
653 EME,
654 J2000,
655 zero.add(Constants.EIGEN5C_EARTH_MU));
656
657 FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
658 SpacecraftState iSR = initialState.toSpacecraftState();
659
660
661 final OrbitType type = OrbitType.KEPLERIAN;
662 double[][] tolerance = NumericalPropagator.tolerances(10.0, FKO.toOrbit(), type);
663 AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
664 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
665 integrator.setInitialStepSize(60);
666
667 AdaptiveStepsizeIntegrator RIntegrator =
668 new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
669 RIntegrator.setInitialStepSize(60);
670
671
672 FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
673 FNP.setOrbitType(type);
674 FNP.setInitialState(initialState);
675
676 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
677 NP.setOrbitType(type);
678 NP.setInitialState(iSR);
679
680 ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
681
682
683 OneAxisEllipsoid earth =
684 new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765,
685 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
686 SolarRadiationPressure forceModel =
687 new SolarRadiationPressure(sun, earth.getEquatorialRadius(),
688 new IsotropicRadiationCNES95Convention(500.0, 0.7, 0.7));
689
690 FNP.addForceModel(forceModel);
691 NP.addForceModel(forceModel);
692
693
694 checkRealFieldPropagation(FKO, PositionAngle.MEAN, 1000., NP, FNP,
695 1.0e-30, 5.0e-10, 3.0e-11, 3.0e-10,
696 1, false);
697 }
698
699
700
701
702
703 @Test
704 public void RealFieldExpectErrorTest() {
705 DSFactory factory = new DSFactory(6, 0);
706 DerivativeStructure a_0 = factory.variable(0, 7e7);
707 DerivativeStructure e_0 = factory.variable(1, 0.4);
708 DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
709 DerivativeStructure R_0 = factory.variable(3, 0.7);
710 DerivativeStructure O_0 = factory.variable(4, 0.5);
711 DerivativeStructure n_0 = factory.variable(5, 0.1);
712
713 Field<DerivativeStructure> field = a_0.getField();
714 DerivativeStructure zero = field.getZero();
715
716 FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
717
718 Frame EME = FramesFactory.getEME2000();
719
720 FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
721 PositionAngle.MEAN,
722 EME,
723 J2000,
724 zero.add(Constants.EIGEN5C_EARTH_MU));
725
726 FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
727
728 SpacecraftState iSR = initialState.toSpacecraftState();
729
730 final OrbitType type = OrbitType.KEPLERIAN;
731 double[][] tolerance = NumericalPropagator.tolerances(0.001, FKO.toOrbit(), type);
732
733
734 AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
735 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
736 integrator.setInitialStepSize(60);
737 AdaptiveStepsizeIntegrator RIntegrator =
738 new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
739 RIntegrator.setInitialStepSize(60);
740
741 FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
742 FNP.setOrbitType(type);
743 FNP.setInitialState(initialState);
744
745 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
746 NP.setInitialState(iSR);
747 ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
748
749
750 OneAxisEllipsoid earth =
751 new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765,
752 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
753 SolarRadiationPressure forceModel =
754 new SolarRadiationPressure(sun, earth.getEquatorialRadius(),
755 new IsotropicRadiationCNES95Convention(500.0, 0.7, 0.7));
756
757 FNP.addForceModel(forceModel);
758
759
760 FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1000.);
761 FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
762 SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
763 FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
764 PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
765
766 Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getX() - finPVC_R.getPosition().getX()) < FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
767 Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getY() - finPVC_R.getPosition().getY()) < FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
768 Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getZ() - finPVC_R.getPosition().getZ()) < FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
769 }
770
771
772
773
774
775
776
777 @Test
778 public void testMoonPenumbra() {
779 final AbsoluteDate date = new AbsoluteDate(2007, 1, 19, 5, 0, 0, TimeScalesFactory.getGPS());
780 final Vector3D p = new Vector3D(12538484.957505366, 15515522.98001655, -17884023.51839292);
781 final Vector3D v = new Vector3D(-3366.9009055533616, 769.5389825219049, -1679.3840677789601);
782 final Orbit orbit = new CartesianOrbit(new TimeStampedPVCoordinates(date, p, v),
783 FramesFactory.getGCRF(), Constants.EIGEN5C_EARTH_MU);
784 doTestMoonEarth(orbit, 7200.0, 1.0, 0, 0, 0, 1513);
785 }
786
787 @Test
788 public void testEarthPenumbraOnly() {
789 final AbsoluteDate date = new AbsoluteDate(2007, 3, 13, 17, 14, 0, TimeScalesFactory.getGPS());
790 final Vector3D p = new Vector3D(-26168647.4977583, -1516554.3304749255, -3206794.210706205);
791 final Vector3D v = new Vector3D(-213.65557094060222, -2377.3633988328584, 3079.4740070013495);
792 final Orbit orbit = new CartesianOrbit(new TimeStampedPVCoordinates(date, p, v),
793 FramesFactory.getGCRF(), Constants.EIGEN5C_EARTH_MU);
794 doTestMoonEarth(orbit, 720.0, 1.0, 0, 531, 0, 0);
795 }
796
797 @Test
798 public void testEarthPenumbraAndUmbra() {
799 final AbsoluteDate date = new AbsoluteDate(2007, 3, 14, 5, 8, 0, TimeScalesFactory.getGPS());
800 final Vector3D p = new Vector3D(-26101379.998276696, -947280.678355501, -3940992.754483608);
801 final Vector3D v = new Vector3D(-348.8911736753223, -2383.738528546711, 3060.9815784341567);
802 final Orbit orbit = new CartesianOrbit(new TimeStampedPVCoordinates(date, p, v),
803 FramesFactory.getGCRF(), Constants.EIGEN5C_EARTH_MU);
804 doTestMoonEarth(orbit, 1200.0, 1.0, 559, 1004, 0, 0);
805 }
806
807 private void doTestMoonEarth(Orbit orbit, double duration, double step,
808 int expectedEarthUmbraSteps, int expectedEarthPenumbraSteps,
809 int expectedMoonUmbraSteps, int expectedMoonPenumbraSteps) {
810
811 Utils.setDataRoot("2007");
812 final ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
813 final ExtendedPVCoordinatesProvider moon = CelestialBodyFactory.getMoon();
814 final SpacecraftState initialState = new SpacecraftState(orbit);
815
816
817 SolarRadiationPressure srp =
818 new SolarRadiationPressure(sun, Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
819 new IsotropicRadiationClassicalConvention(50.0, 0.5, 0.5));
820 srp.addOccultingBody(moon, Constants.MOON_EQUATORIAL_RADIUS);
821
822
823 final OrbitType type = OrbitType.CARTESIAN;
824 double[][] tol = NumericalPropagator.tolerances(0.1, orbit, type);
825 AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(600, 1000000, tol[0], tol[1]);
826 final NumericalPropagator propagator = new NumericalPropagator(integrator);
827 propagator.setOrbitType(type);
828 propagator.setInitialState(initialState);
829 propagator.addForceModel(srp);
830 MoonEclipseStepHandler handler = new MoonEclipseStepHandler(moon, sun, srp);
831 propagator.setStepHandler(step, handler);
832 propagator.propagate(orbit.getDate().shiftedBy(duration));
833 Assert.assertEquals(expectedEarthUmbraSteps, handler.earthUmbraSteps);
834 Assert.assertEquals(expectedEarthPenumbraSteps, handler.earthPenumbraSteps);
835 Assert.assertEquals(expectedMoonUmbraSteps, handler.moonUmbraSteps);
836 Assert.assertEquals(expectedMoonPenumbraSteps, handler.moonPenumbraSteps);
837 }
838
839
840
841
842
843 private static class MoonEclipseStepHandler implements OrekitFixedStepHandler {
844
845 final ExtendedPVCoordinatesProvider moon;
846 final ExtendedPVCoordinatesProvider sun;
847 final SolarRadiationPressure srp;
848 int earthUmbraSteps;
849 int earthPenumbraSteps;
850 int moonUmbraSteps;
851 int moonPenumbraSteps;
852
853
854
855 MoonEclipseStepHandler(final ExtendedPVCoordinatesProvider moon, final ExtendedPVCoordinatesProvider sun,
856 final SolarRadiationPressure srp) {
857 this.moon = moon;
858 this.sun = sun;
859 this.srp = srp;
860
861 }
862
863
864 @Override
865 public void init(final SpacecraftState s0, final AbsoluteDate t, final double step) {
866 earthUmbraSteps = 0;
867 earthPenumbraSteps = 0;
868 moonUmbraSteps = 0;
869 moonPenumbraSteps = 0;
870 }
871
872
873 @Override
874 public void handleStep(final SpacecraftState currentState) {
875 final AbsoluteDate date = currentState.getDate();
876 final Frame frame = currentState.getFrame();
877
878 final Vector3D moonPos = moon.getPVCoordinates(date, frame).getPosition();
879 final Vector3D sunPos = sun.getPVCoordinates(date, frame).getPosition();
880 final Vector3D statePos = currentState.getPVCoordinates().getPosition();
881
882
883 final double[] moonAngles = srp.getGeneralEclipseAngles(statePos, moonPos, Constants.MOON_EQUATORIAL_RADIUS,
884 sunPos, Constants.SUN_RADIUS);
885 final double moonUmbra = moonAngles[0] - moonAngles[1] + moonAngles[2];
886 final boolean isInMoonUmbra = (moonUmbra < 1.0e-10);
887 if (isInMoonUmbra) {
888 ++moonUmbraSteps;
889 }
890
891 final double moonPenumbra = moonAngles[0] - moonAngles[1] - moonAngles[2];
892 final boolean isInMoonPenumbra = (moonPenumbra < -1.0e-10);
893 if (isInMoonPenumbra) {
894 ++moonPenumbraSteps;
895 }
896
897
898 final double[] earthAngles = srp.getEclipseAngles(sunPos, statePos);
899
900 final double earthUmbra = earthAngles[0] - earthAngles[1] + earthAngles[2];
901 final boolean isInEarthUmbra = (earthUmbra < 1.0e-10);
902 if (isInEarthUmbra) {
903 ++earthUmbraSteps;
904 }
905
906 final double earthPenumbra = earthAngles[0] - earthAngles[1] - earthAngles[2];
907 final boolean isInEarthPenumbra = (earthPenumbra < -1.0e-10);
908 if (isInEarthPenumbra) {
909 ++earthPenumbraSteps;
910 }
911
912
913
914 final double lightingRatio = srp.getTotalLightingRatio(statePos, frame, date);
915
916
917 if (isInMoonUmbra || isInEarthUmbra) {
918 Assert.assertEquals(0.0, lightingRatio, 1e-8);
919 }
920 else if (isInMoonPenumbra || isInEarthPenumbra) {
921 Assert.assertTrue(lightingRatio < 1.0);
922 Assert.assertTrue(lightingRatio > 0.0);
923 }
924 else {
925 Assert.assertEquals(1.0, lightingRatio, 1e-8);
926 }
927 }
928 }
929
930
931
932
933
934
935
936 @Test
937 public void testFieldMoonPenumbra() {
938 Field<Decimal64> field = Decimal64Field.getInstance();
939 final FieldAbsoluteDate<Decimal64> date = new FieldAbsoluteDate<>(field, 2007, 1, 19, 5, 0, 0, TimeScalesFactory.getGPS());
940 final FieldVector3D<Decimal64> p = new FieldVector3D<>(field, new Vector3D(12538484.957505366, 15515522.98001655, -17884023.51839292));
941 final FieldVector3D<Decimal64> v = new FieldVector3D<>(field, new Vector3D(-3366.9009055533616, 769.5389825219049, -1679.3840677789601));
942 final FieldVector3D<Decimal64> a = FieldVector3D.getZero(field);
943 final FieldOrbit<Decimal64> orbit = new FieldCartesianOrbit<>(new TimeStampedFieldPVCoordinates<>(date, p, v, a),
944 FramesFactory.getGCRF(),
945 field.getZero().newInstance(Constants.EIGEN5C_EARTH_MU));
946 doTestFieldMoonEarth(orbit, field.getZero().newInstance(7200.0), field.getZero().newInstance(1.0), 0, 0, 0, 1513);
947 }
948
949 @Test
950 public void testFieldEarthPenumbraOnly() {
951 Field<Decimal64> field = Decimal64Field.getInstance();
952 final FieldAbsoluteDate<Decimal64> date = new FieldAbsoluteDate<>(field, 2007, 3, 13, 17, 14, 0, TimeScalesFactory.getGPS());
953 final FieldVector3D<Decimal64> p = new FieldVector3D<>(field, new Vector3D(-26168647.4977583, -1516554.3304749255, -3206794.210706205));
954 final FieldVector3D<Decimal64> v = new FieldVector3D<>(field, new Vector3D(-213.65557094060222, -2377.3633988328584, 3079.4740070013495));
955 final FieldVector3D<Decimal64> a = FieldVector3D.getZero(field);
956 final FieldOrbit<Decimal64> orbit = new FieldCartesianOrbit<>(new TimeStampedFieldPVCoordinates<>(date, p, v, a),
957 FramesFactory.getGCRF(),
958 field.getZero().newInstance(Constants.EIGEN5C_EARTH_MU));
959 doTestFieldMoonEarth(orbit, field.getZero().newInstance(720.0), field.getZero().newInstance(1.0), 0, 531, 0, 0);
960 }
961
962 @Test
963 public void testFieldEarthPenumbraAndUmbra() {
964 Field<Decimal64> field = Decimal64Field.getInstance();
965 final FieldAbsoluteDate<Decimal64> date = new FieldAbsoluteDate<>(field, 2007, 3, 14, 5, 8, 0, TimeScalesFactory.getGPS());
966 final FieldVector3D<Decimal64> p = new FieldVector3D<>(field, new Vector3D(-26101379.998276696, -947280.678355501, -3940992.754483608));
967 final FieldVector3D<Decimal64> v = new FieldVector3D<>(field, new Vector3D(-348.8911736753223, -2383.738528546711, 3060.9815784341567));
968 final FieldVector3D<Decimal64> a = FieldVector3D.getZero(field);
969 final FieldOrbit<Decimal64> orbit = new FieldCartesianOrbit<>(new TimeStampedFieldPVCoordinates<>(date, p, v, a),
970 FramesFactory.getGCRF(),
971 field.getZero().newInstance(Constants.EIGEN5C_EARTH_MU));
972 doTestFieldMoonEarth(orbit, field.getZero().newInstance(1200.0), field.getZero().newInstance(1.0), 559, 1004, 0, 0);
973 }
974
975 private <T extends CalculusFieldElement<T>> void doTestFieldMoonEarth(FieldOrbit<T> orbit, T duration, T step,
976 int expectedEarthUmbraSteps, int expectedEarthPenumbraSteps,
977 int expectedMoonUmbraSteps, int expectedMoonPenumbraSteps) {
978
979 Utils.setDataRoot("2007");
980 final Field<T> field = orbit.getDate().getField();
981 final ExtendedPVCoordinatesProvider sun = CelestialBodyFactory.getSun();
982 final ExtendedPVCoordinatesProvider moon = CelestialBodyFactory.getMoon();
983 final FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(orbit);
984
985
986 SolarRadiationPressure srp =
987 new SolarRadiationPressure(sun, Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
988 new IsotropicRadiationClassicalConvention(50.0, 0.5, 0.5));
989 srp.addOccultingBody(moon, Constants.MOON_EQUATORIAL_RADIUS);
990
991
992 final OrbitType type = OrbitType.CARTESIAN;
993 double[][] tol = FieldNumericalPropagator.tolerances(field.getZero().newInstance(0.1), orbit, type);
994 AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 600, 1000000, tol[0], tol[1]);
995 final FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
996 propagator.setOrbitType(type);
997 propagator.setInitialState(initialState);
998 propagator.addForceModel(srp);
999 FieldMoonEclipseStepHandler<T> handler = new FieldMoonEclipseStepHandler<>(moon, sun, srp);
1000 propagator.setStepHandler(step, handler);
1001 propagator.propagate(orbit.getDate().shiftedBy(duration));
1002 Assert.assertEquals(expectedEarthUmbraSteps, handler.earthUmbraSteps);
1003 Assert.assertEquals(expectedEarthPenumbraSteps, handler.earthPenumbraSteps);
1004 Assert.assertEquals(expectedMoonUmbraSteps, handler.moonUmbraSteps);
1005 Assert.assertEquals(expectedMoonPenumbraSteps, handler.moonPenumbraSteps);
1006 }
1007
1008
1009
1010
1011
1012 private static class FieldMoonEclipseStepHandler<T extends CalculusFieldElement<T>> implements FieldOrekitFixedStepHandler<T> {
1013
1014 final ExtendedPVCoordinatesProvider moon;
1015 final ExtendedPVCoordinatesProvider sun;
1016 final SolarRadiationPressure srp;
1017 int earthUmbraSteps;
1018 int earthPenumbraSteps;
1019 int moonUmbraSteps;
1020 int moonPenumbraSteps;
1021
1022
1023
1024 FieldMoonEclipseStepHandler(final ExtendedPVCoordinatesProvider moon, final ExtendedPVCoordinatesProvider sun,
1025 final SolarRadiationPressure srp) {
1026 this.moon = moon;
1027 this.sun = sun;
1028 this.srp = srp;
1029
1030 }
1031
1032
1033 @Override
1034 public void init(final FieldSpacecraftState<T> s0, final FieldAbsoluteDate<T> t, final T step) {
1035 earthUmbraSteps = 0;
1036 earthPenumbraSteps = 0;
1037 moonUmbraSteps = 0;
1038 moonPenumbraSteps = 0;
1039 }
1040
1041
1042 @Override
1043 public void handleStep(final FieldSpacecraftState<T> currentState) {
1044 final FieldAbsoluteDate<T> date = currentState.getDate();
1045 final Frame frame = currentState.getFrame();
1046
1047 final FieldVector3D<T> moonPos = moon.getPVCoordinates(date, frame).getPosition();
1048 final FieldVector3D<T> sunPos = sun.getPVCoordinates(date, frame).getPosition();
1049 final FieldVector3D<T> statePos = currentState.getPVCoordinates().getPosition();
1050
1051
1052 final T[] moonAngles = srp.getGeneralEclipseAngles(statePos, moonPos,
1053 date.getField().getZero().newInstance(Constants.MOON_EQUATORIAL_RADIUS),
1054 sunPos,
1055 date.getField().getZero().newInstance(Constants.SUN_RADIUS));
1056 final T moonUmbra = moonAngles[0].subtract(moonAngles[1]).add(moonAngles[2]);
1057 final boolean isInMoonUmbra = (moonUmbra.getReal() < 1.0e-10);
1058 if (isInMoonUmbra) {
1059 ++moonUmbraSteps;
1060 }
1061
1062 final T moonPenumbra = moonAngles[0].subtract(moonAngles[1]).subtract(moonAngles[2]);
1063 final boolean isInMoonPenumbra = (moonPenumbra.getReal() < -1.0e-10);
1064 if (isInMoonPenumbra) {
1065 ++moonPenumbraSteps;
1066 }
1067
1068
1069 final T[] earthAngles = srp.getEclipseAngles(sunPos, statePos);
1070
1071 final T earthUmbra = earthAngles[0].subtract(earthAngles[1]).add(earthAngles[2]);
1072 final boolean isInEarthUmbra = (earthUmbra.getReal() < 1.0e-10);
1073 if (isInEarthUmbra) {
1074 ++earthUmbraSteps;
1075 }
1076
1077 final T earthPenumbra = earthAngles[0].subtract(earthAngles[1]).subtract(earthAngles[2]);
1078 final boolean isInEarthPenumbra = (earthPenumbra.getReal() < -1.0e-10);
1079 if (isInEarthPenumbra) {
1080 ++earthPenumbraSteps;
1081 }
1082
1083
1084 final T lightingRatio = srp.getTotalLightingRatio(statePos, frame, date);
1085
1086
1087 if (isInMoonUmbra || isInEarthUmbra) {
1088 Assert.assertEquals(0.0, lightingRatio.getReal(), 1e-8);
1089 }
1090 else if (isInMoonPenumbra || isInEarthPenumbra) {
1091 Assert.assertTrue(lightingRatio.getReal() < 1.0);
1092 Assert.assertTrue(lightingRatio.getReal() > 0.0);
1093 }
1094 else {
1095 Assert.assertEquals(1.0, lightingRatio.getReal(), 1e-8);
1096 }
1097 }
1098 }
1099 @Before
1100 public void setUp() {
1101 Utils.setDataRoot("regular-data");
1102 }
1103
1104 }