1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.forces.radiation;
18
19 import java.text.ParseException;
20 import java.util.Arrays;
21 import java.util.List;
22
23 import org.hipparchus.CalculusFieldElement;
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.FieldVector3D;
29 import org.hipparchus.geometry.euclidean.threed.Vector3D;
30 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
31 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
32 import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
33 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
34 import org.hipparchus.util.Binary64;
35 import org.hipparchus.util.Binary64Field;
36 import org.hipparchus.util.FastMath;
37 import org.junit.jupiter.api.Assertions;
38 import org.junit.jupiter.api.BeforeEach;
39 import org.junit.jupiter.api.Test;
40 import org.mockito.Mockito;
41 import org.orekit.Utils;
42 import org.orekit.attitudes.LofOffset;
43 import org.orekit.bodies.CelestialBody;
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.forces.gravity.HolmesFeatherstoneAttractionModel;
51 import org.orekit.forces.gravity.ThirdBodyAttraction;
52 import org.orekit.forces.gravity.potential.GravityFieldFactory;
53 import org.orekit.forces.gravity.potential.ICGEMFormatReader;
54 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
55 import org.orekit.frames.Frame;
56 import org.orekit.frames.FramesFactory;
57 import org.orekit.frames.LOFType;
58 import org.orekit.orbits.CartesianOrbit;
59 import org.orekit.orbits.CircularOrbit;
60 import org.orekit.orbits.EquinoctialOrbit;
61 import org.orekit.orbits.FieldCartesianOrbit;
62 import org.orekit.orbits.FieldKeplerianOrbit;
63 import org.orekit.orbits.FieldOrbit;
64 import org.orekit.orbits.KeplerianOrbit;
65 import org.orekit.orbits.Orbit;
66 import org.orekit.orbits.OrbitType;
67 import org.orekit.orbits.PositionAngleType;
68 import org.orekit.propagation.FieldSpacecraftState;
69 import org.orekit.propagation.PropagatorsParallelizer;
70 import org.orekit.propagation.SpacecraftState;
71 import org.orekit.propagation.ToleranceProvider;
72 import org.orekit.propagation.analytical.KeplerianPropagator;
73 import org.orekit.propagation.numerical.FieldNumericalPropagator;
74 import org.orekit.propagation.numerical.NumericalPropagator;
75 import org.orekit.propagation.sampling.FieldOrekitFixedStepHandler;
76 import org.orekit.propagation.sampling.MultiSatStepHandler;
77 import org.orekit.propagation.sampling.OrekitFixedStepHandler;
78 import org.orekit.propagation.sampling.OrekitStepInterpolator;
79 import org.orekit.time.AbsoluteDate;
80 import org.orekit.time.DateComponents;
81 import org.orekit.time.FieldAbsoluteDate;
82 import org.orekit.time.TimeComponents;
83 import org.orekit.time.TimeScalesFactory;
84 import org.orekit.utils.*;
85
86
87 public class SolarRadiationPressureTest extends AbstractLegacyForceModelTest {
88
89 @Override
90 protected FieldVector3D<DerivativeStructure> accelerationDerivatives(final ForceModel forceModel,
91 final FieldSpacecraftState<DerivativeStructure> state) {
92 try {
93 final AbsoluteDate date = state.getDate().toAbsoluteDate();
94 final FieldVector3D<DerivativeStructure> position = state.getPVCoordinates().getPosition();
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
105 final Field<DerivativeStructure> field = position.getX().getField();
106 final FieldVector3D<DerivativeStructure> sunSatVector = position.subtract(sun.getPosition(date, state.getFrame()));
107 final DerivativeStructure r2 = sunSatVector.getNormSq();
108
109
110 final DerivativeStructure ratio = ((SolarRadiationPressure) forceModel).getLightingRatio(state);
111 final DerivativeStructure rawP = ratio.multiply(kRef).divide(r2);
112 final FieldVector3D<DerivativeStructure> flux = new FieldVector3D<>(rawP.divide(r2.sqrt()), sunSatVector);
113
114
115 return spacecraft.radiationPressureAcceleration(state, flux, forceModel.getParameters(field));
116 } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
117 return null;
118 }
119 }
120
121 @Override
122 protected FieldVector3D<Gradient> accelerationDerivativesGradient(final ForceModel forceModel,
123 final FieldSpacecraftState<Gradient> state) {
124 try {
125 final AbsoluteDate date = state.getDate().toAbsoluteDate();
126 final FieldVector3D<Gradient> position = state.getPVCoordinates().getPosition();
127 java.lang.reflect.Field kRefField = SolarRadiationPressure.class.getDeclaredField("kRef");
128 kRefField.setAccessible(true);
129 double kRef = kRefField.getDouble(forceModel);
130 java.lang.reflect.Field sunField = SolarRadiationPressure.class.getDeclaredField("sun");
131 sunField.setAccessible(true);
132 PVCoordinatesProvider sun = (PVCoordinatesProvider) sunField.get(forceModel);
133 java.lang.reflect.Field spacecraftField = SolarRadiationPressure.class.getDeclaredField("spacecraft");
134 spacecraftField.setAccessible(true);
135 RadiationSensitive spacecraft = (RadiationSensitive) spacecraftField.get(forceModel);
136
137 final Field<Gradient> field = position.getX().getField();
138 final FieldVector3D<Gradient> sunSatVector = position.subtract(sun.getPosition(date, state.getFrame()));
139 final Gradient r2 = sunSatVector.getNormSq();
140
141
142 final Gradient ratio = ((SolarRadiationPressure) forceModel).getLightingRatio(state);
143 final Gradient rawP = ratio.multiply(kRef).divide(r2);
144 final FieldVector3D<Gradient> flux = new FieldVector3D<>(rawP.divide(r2.sqrt()), sunSatVector);
145
146
147 return spacecraft.radiationPressureAcceleration(state, flux, forceModel.getParameters(field));
148 } catch (IllegalArgumentException | IllegalAccessException | NoSuchFieldException | SecurityException e) {
149 return null;
150 }
151 }
152
153 @Test
154 void testLightingInterplanetary() throws ParseException {
155
156 AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 3, 21),
157 new TimeComponents(13, 59, 27.816),
158 TimeScalesFactory.getUTC());
159 Orbit orbit = new KeplerianOrbit(1.0e11, 0.1, 0.2, 0.3, 0.4, 0.5, PositionAngleType.TRUE,
160 FramesFactory.getICRF(), date, Constants.JPL_SSD_SUN_GM);
161 ExtendedPositionProvider sun = CelestialBodyFactory.getSun();
162 SolarRadiationPressure srp =
163 new SolarRadiationPressure(sun, new OneAxisEllipsoid(1.0e-10, 0.0, FramesFactory.getICRF()),
164 new IsotropicRadiationClassicalConvention(50.0, 0.5, 0.5));
165
166
167 Assertions.assertEquals(1.0,
168 srp.getLightingRatio(new SpacecraftState(orbit)),
169 1.0e-15);
170
171
172 final FieldSpacecraftState<Binary64> fieldSc = new FieldSpacecraftState<>(Binary64Field.getInstance(),
173 new SpacecraftState(orbit));
174 Assertions.assertEquals(1.0,
175 srp.getLightingRatio(fieldSc).getReal(),
176 1.0e-15);
177 }
178
179 @Test
180 void testLighting() throws ParseException {
181
182
183 AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 3, 21),
184 new TimeComponents(13, 59, 27.816),
185 TimeScalesFactory.getUTC());
186 Orbit orbit = new EquinoctialOrbit(42164000, 10e-3, 10e-3,
187 FastMath.tan(0.001745329)*FastMath.cos(2*FastMath.PI/3), FastMath.tan(0.001745329)*FastMath.sin(2*FastMath.PI/3),
188 0.1, PositionAngleType.TRUE, FramesFactory.getEME2000(), date, mu);
189 CelestialBody sun = CelestialBodyFactory.getSun();
190 OneAxisEllipsoid earth =
191 new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765,
192 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
193 SolarRadiationPressure SRP =
194 new SolarRadiationPressure(sun, earth,
195 new IsotropicRadiationCNES95Convention(50.0, 0.5, 0.5));
196
197 double period = 2*FastMath.PI*FastMath.sqrt(orbit.getA()*orbit.getA()*orbit.getA()/orbit.getMu());
198 Assertions.assertEquals(86164, period, 1);
199
200
201 KeplerianPropagator k = new KeplerianPropagator(orbit);
202
203
204 AbsoluteDate currentDate;
205 double changed = 1;
206 int count=0;
207
208
209 for (int t=1;t<3*period;t+=1000) {
210 currentDate = date.shiftedBy(t);
211 try {
212
213 double ratio = SRP.getLightingRatio(k.propagate(currentDate));
214
215 if (FastMath.floor(ratio)!=changed) {
216 changed = FastMath.floor(ratio);
217 if (changed == 0) {
218 count++;
219 }
220 }
221 } catch (OrekitException e) {
222 e.printStackTrace();
223 }
224 }
225 Assertions.assertEquals(3, count);
226 }
227
228 @Test
229 void testGlobalStateJacobianIsotropicSingle() {
230
231
232 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
233 new TimeComponents(13, 59, 27.816),
234 TimeScalesFactory.getUTC());
235 double i = FastMath.toRadians(98.7);
236 double omega = FastMath.toRadians(93.0);
237 double OMEGA = FastMath.toRadians(15.0 * 22.5);
238 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
239 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
240 Constants.EIGEN5C_EARTH_MU);
241 OrbitType integrationType = OrbitType.CARTESIAN;
242 double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(0.01).getTolerances(orbit, integrationType);
243
244 NumericalPropagator propagator =
245 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
246 tolerances[0], tolerances[1]));
247 propagator.setOrbitType(integrationType);
248 SolarRadiationPressure forceModel =
249 new SolarRadiationPressure(CelestialBodyFactory.getSun(),
250 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
251 Constants.WGS84_EARTH_FLATTENING,
252 FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
253 new IsotropicRadiationSingleCoefficient(2.5, 0.7));
254 propagator.addForceModel(forceModel);
255 SpacecraftState state0 = new SpacecraftState(orbit);
256
257 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
258 1e3, tolerances[0], 3.2e-5);
259
260 }
261
262 @Test
263 void testLocalJacobianIsotropicClassicalVs80Implementation() {
264
265
266 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
267 new TimeComponents(13, 59, 27.816),
268 TimeScalesFactory.getUTC());
269 double i = FastMath.toRadians(98.7);
270 double omega = FastMath.toRadians(93.0);
271 double OMEGA = FastMath.toRadians(15.0 * 22.5);
272 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
273 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
274 Constants.EIGEN5C_EARTH_MU);
275 final SolarRadiationPressure forceModel =
276 new SolarRadiationPressure(CelestialBodyFactory.getSun(),
277 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
278 Constants.WGS84_EARTH_FLATTENING,
279 FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
280 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
281
282 checkStateJacobianVs80Implementation(new SpacecraftState(orbit), forceModel,
283 new LofOffset(orbit.getFrame(), LOFType.LVLH_CCSDS),
284 1.0e-15, false);
285
286 }
287
288 @Test
289 void testLocalJacobianIsotropicClassicalVs80ImplementationGradient() {
290
291
292 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
293 new TimeComponents(13, 59, 27.816),
294 TimeScalesFactory.getUTC());
295 double i = FastMath.toRadians(98.7);
296 double omega = FastMath.toRadians(93.0);
297 double OMEGA = FastMath.toRadians(15.0 * 22.5);
298 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
299 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
300 Constants.EIGEN5C_EARTH_MU);
301 final SolarRadiationPressure forceModel =
302 new SolarRadiationPressure(CelestialBodyFactory.getSun(),
303 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
304 Constants.WGS84_EARTH_FLATTENING,
305 FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
306 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
307
308 checkStateJacobianVs80ImplementationGradient(new SpacecraftState(orbit), forceModel,
309 new LofOffset(orbit.getFrame(), LOFType.LVLH_CCSDS),
310 1.0e-15, false);
311
312 }
313
314 @Test
315 void testLocalJacobianIsotropicClassicalVsFiniteDifferencesFullLight() {
316
317 doTestLocalJacobianIsotropicClassicalVsFiniteDifferences(250.0, 1000.0, 3.0e-8, false);
318 }
319
320 @Test
321 void testLocalJacobianIsotropicClassicalVsFiniteDifferencesGradientFullLight() {
322
323 doTestLocalJacobianIsotropicClassicalVsFiniteDifferencesGradient(250.0, 1000.0, 3.0e-8, false);
324 }
325
326 @Test
327 void testLocalJacobianIsotropicClassicalVsFiniteDifferencesPenumbra() {
328
329
330 doTestLocalJacobianIsotropicClassicalVsFiniteDifferences(275.5, 100.0, 8.0e-7, false);
331 }
332
333 @Test
334 void testLocalJacobianIsotropicClassicalVsFiniteDifferencesGradientPenumbra() {
335
336
337 doTestLocalJacobianIsotropicClassicalVsFiniteDifferencesGradient(275.5, 100.0, 8.0e-7, false);
338 }
339
340 @Test
341 void testLocalJacobianIsotropicClassicalVsFiniteDifferencesEclipse() {
342
343 doTestLocalJacobianIsotropicClassicalVsFiniteDifferences(300.0, 1000.0, 1.0e-50, false);
344 }
345
346 @Test
347 void testLocalJacobianIsotropicClassicalVsFiniteDifferencesGradientEclipse() {
348
349 doTestLocalJacobianIsotropicClassicalVsFiniteDifferencesGradient(300.0, 1000.0, 1.0e-50, false);
350 }
351
352 private void doTestLocalJacobianIsotropicClassicalVsFiniteDifferences(double deltaT, double dP,
353 double checkTolerance,
354 boolean print) {
355
356
357 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
358 new TimeComponents(13, 59, 27.816),
359 TimeScalesFactory.getUTC());
360 double i = FastMath.toRadians(98.7);
361 double omega = FastMath.toRadians(93.0);
362 double OMEGA = FastMath.toRadians(15.0 * 22.5);
363 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
364 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
365 Constants.EIGEN5C_EARTH_MU);
366 final SolarRadiationPressure forceModel =
367 new SolarRadiationPressure(CelestialBodyFactory.getSun(),
368 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
369 Constants.WGS84_EARTH_FLATTENING,
370 FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
371 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
372
373 checkStateJacobianVsFiniteDifferences(new SpacecraftState(orbit.shiftedBy(deltaT)), forceModel,
374 Utils.defaultLaw(), dP, checkTolerance, print);
375
376 }
377
378 private void doTestLocalJacobianIsotropicClassicalVsFiniteDifferencesGradient(double deltaT, double dP,
379 double checkTolerance,
380 boolean print) {
381
382
383 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
384 new TimeComponents(13, 59, 27.816),
385 TimeScalesFactory.getUTC());
386 double i = FastMath.toRadians(98.7);
387 double omega = FastMath.toRadians(93.0);
388 double OMEGA = FastMath.toRadians(15.0 * 22.5);
389 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
390 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
391 Constants.EIGEN5C_EARTH_MU);
392 final SolarRadiationPressure forceModel =
393 new SolarRadiationPressure(CelestialBodyFactory.getSun(),
394 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
395 Constants.WGS84_EARTH_FLATTENING,
396 FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
397 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
398
399 checkStateJacobianVsFiniteDifferencesGradient(new SpacecraftState(orbit.shiftedBy(deltaT)), forceModel,
400 Utils.defaultLaw(), dP, checkTolerance, print);
401
402 }
403
404 @Test
405 void testGlobalStateJacobianIsotropicClassical() {
406
407
408 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
409 new TimeComponents(13, 59, 27.816),
410 TimeScalesFactory.getUTC());
411 double i = FastMath.toRadians(98.7);
412 double omega = FastMath.toRadians(93.0);
413 double OMEGA = FastMath.toRadians(15.0 * 22.5);
414 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
415 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
416 Constants.EIGEN5C_EARTH_MU);
417 OrbitType integrationType = OrbitType.CARTESIAN;
418 double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(0.01).getTolerances(orbit, integrationType);
419
420 NumericalPropagator propagator =
421 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
422 tolerances[0], tolerances[1]));
423 propagator.setOrbitType(integrationType);
424 SolarRadiationPressure forceModel =
425 new SolarRadiationPressure(CelestialBodyFactory.getSun(),
426 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
427 Constants.WGS84_EARTH_FLATTENING,
428 FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
429 new IsotropicRadiationClassicalConvention(2.5, 0.7, 0.2));
430 propagator.addForceModel(forceModel);
431 SpacecraftState state0 = new SpacecraftState(orbit);
432
433 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
434 1e6, tolerances[0], 4.2e-5);
435
436 }
437
438 @Test
439 void testGlobalStateJacobianIsotropicCnes() {
440
441
442 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
443 new TimeComponents(13, 59, 27.816),
444 TimeScalesFactory.getUTC());
445 double i = FastMath.toRadians(98.7);
446 double omega = FastMath.toRadians(93.0);
447 double OMEGA = FastMath.toRadians(15.0 * 22.5);
448 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
449 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
450 Constants.EIGEN5C_EARTH_MU);
451 OrbitType integrationType = OrbitType.CARTESIAN;
452 double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(0.01).getTolerances(orbit, integrationType);
453
454 NumericalPropagator propagator =
455 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
456 tolerances[0], tolerances[1]));
457 propagator.setOrbitType(integrationType);
458 SolarRadiationPressure forceModel =
459 new SolarRadiationPressure(CelestialBodyFactory.getSun(),
460 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
461 Constants.WGS84_EARTH_FLATTENING,
462 FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
463 new IsotropicRadiationCNES95Convention(2.5, 0.7, 0.2));
464 propagator.addForceModel(forceModel);
465 SpacecraftState state0 = new SpacecraftState(orbit);
466
467 checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0),
468 1e3, tolerances[0], 3.0e-5);
469
470 }
471
472 @Test
473 void testParameterDerivativeBox() {
474
475 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
476 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
477 final SpacecraftState state =
478 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
479 FramesFactory.getGCRF(),
480 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
481 Constants.EIGEN5C_EARTH_MU));
482
483 SolarRadiationPressure forceModel =
484 new SolarRadiationPressure(CelestialBodyFactory.getSun(),
485 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
486 Constants.WGS84_EARTH_FLATTENING,
487 FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
488 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
489 Vector3D.PLUS_J, 1.2, 0.0, 0.7, 0.2));
490
491 checkParameterDerivative(state, forceModel, RadiationSensitive.GLOBAL_RADIATION_FACTOR, 0.25, 8.8e-16);
492
493 }
494
495 @Test
496 void testParameterDerivativeGradientBox() {
497
498 final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
499 final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
500 final SpacecraftState state =
501 new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel),
502 FramesFactory.getGCRF(),
503 new AbsoluteDate(2003, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()),
504 Constants.EIGEN5C_EARTH_MU));
505
506 SolarRadiationPressure forceModel =
507 new SolarRadiationPressure(CelestialBodyFactory.getSun(),
508 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
509 Constants.WGS84_EARTH_FLATTENING,
510 FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
511 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
512 Vector3D.PLUS_J, 1.2, 0.0, 0.7, 0.2));
513
514 checkParameterDerivativeGradient(state, forceModel, RadiationSensitive.GLOBAL_RADIATION_FACTOR, 0.25, 8.8e-16);
515
516 }
517
518 @Test
519 void testGlobalStateJacobianBox() {
520
521
522 AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01),
523 new TimeComponents(13, 59, 27.816),
524 TimeScalesFactory.getUTC());
525 double i = FastMath.toRadians(98.7);
526 double omega = FastMath.toRadians(93.0);
527 double OMEGA = FastMath.toRadians(15.0 * 22.5);
528 Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
529 0, PositionAngleType.MEAN, FramesFactory.getEME2000(), date,
530 Constants.EIGEN5C_EARTH_MU);
531 OrbitType integrationType = OrbitType.CARTESIAN;
532 double[][] tolerances = ToleranceProvider.getDefaultToleranceProvider(0.01).getTolerances(orbit, integrationType);
533
534 NumericalPropagator propagator =
535 new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120,
536 tolerances[0], tolerances[1]));
537 propagator.setOrbitType(integrationType);
538 SolarRadiationPressure forceModel =
539 new SolarRadiationPressure(CelestialBodyFactory.getSun(),
540 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
541 Constants.WGS84_EARTH_FLATTENING,
542 FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
543 new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0,
544 Vector3D.PLUS_J, 1.2, 0.0, 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 void testRoughOrbitalModifs() {
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, PositionAngleType.TRUE, FramesFactory.getEME2000(), date, mu);
564 final double period = orbit.getKeplerianPeriod();
565 Assertions.assertEquals(86164, period, 1);
566 ExtendedPositionProvider 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, new IsotropicRadiationCNES95Convention(500.0, 0.7, 0.7));
574
575
576 double[] absTolerance = {
577 0.1, 1.0e-9, 1.0e-9, 1.0e-5, 1.0e-5, 1.0e-5, 0.001
578 };
579 double[] relTolerance = {
580 1.0e-4, 1.0e-4, 1.0e-4, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-7
581 };
582 AdaptiveStepsizeIntegrator integrator =
583 new DormandPrince853Integrator(900.0, 60000, absTolerance, relTolerance);
584 integrator.setInitialStepSize(3600);
585 final NumericalPropagator calc = new NumericalPropagator(integrator);
586 calc.addForceModel(SRP);
587
588
589 calc.setStepHandler(FastMath.floor(period), new SolarStepHandler());
590 AbsoluteDate finalDate = date.shiftedBy(10 * period);
591 calc.setInitialState(new SpacecraftState(orbit, 1500.0));
592 calc.propagate(finalDate);
593 Assertions.assertTrue(calc.getCalls() < 7100);
594 }
595
596 public static void checkRadius(double radius , double min , double max) {
597 Assertions.assertTrue(radius >= min);
598 Assertions.assertTrue(radius <= max);
599 }
600
601 private double mu = 3.98600E14;
602
603 private static class SolarStepHandler implements OrekitFixedStepHandler {
604
605 @Override
606 public void handleStep(SpacecraftState currentState) {
607 final double dex = currentState.getOrbit().getEquinoctialEx() - 0.01071166;
608 final double dey = currentState.getOrbit().getEquinoctialEy() - 0.00654848;
609 final double alpha = FastMath.toDegrees(FastMath.atan2(dey, dex));
610 Assertions.assertTrue(alpha > 100.0);
611 Assertions.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 PositionAngleType.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 = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(FKO.toOrbit(), type);
656 AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
657 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
658 integrator.setInitialStepSize(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 ExtendedPositionProvider 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, new IsotropicRadiationCNES95Convention(500.0, 0.7, 0.7));
681
682 FNP.addForceModel(forceModel);
683 NP.addForceModel(forceModel);
684
685
686 checkRealFieldPropagation(FKO, PositionAngleType.MEAN, 1000., NP, FNP,
687 1.0e-30, 5.0e-10, 3.0e-11, 3.0e-10,
688 1, false);
689 }
690
691
692
693
694
695 @Test
696 public void RealFieldExpectErrorTest() {
697 DSFactory factory = new DSFactory(6, 0);
698 DerivativeStructure a_0 = factory.variable(0, 7e7);
699 DerivativeStructure e_0 = factory.variable(1, 0.4);
700 DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
701 DerivativeStructure R_0 = factory.variable(3, 0.7);
702 DerivativeStructure O_0 = factory.variable(4, 0.5);
703 DerivativeStructure n_0 = factory.variable(5, 0.1);
704
705 Field<DerivativeStructure> field = a_0.getField();
706 DerivativeStructure zero = field.getZero();
707
708 FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
709
710 Frame EME = FramesFactory.getEME2000();
711
712 FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0,
713 PositionAngleType.MEAN,
714 EME,
715 J2000,
716 zero.add(Constants.EIGEN5C_EARTH_MU));
717
718 FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
719
720 SpacecraftState iSR = initialState.toSpacecraftState();
721
722 final OrbitType type = OrbitType.KEPLERIAN;
723 double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(0.001).getTolerances(FKO.toOrbit(), type);
724
725
726 AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator =
727 new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
728 integrator.setInitialStepSize(60);
729 AdaptiveStepsizeIntegrator RIntegrator =
730 new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
731 RIntegrator.setInitialStepSize(60);
732
733 FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
734 FNP.setOrbitType(type);
735 FNP.setInitialState(initialState);
736
737 NumericalPropagator NP = new NumericalPropagator(RIntegrator);
738 NP.setInitialState(iSR);
739 ExtendedPositionProvider sun = CelestialBodyFactory.getSun();
740
741
742 OneAxisEllipsoid earth =
743 new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765,
744 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
745 SolarRadiationPressure forceModel =
746 new SolarRadiationPressure(sun, earth, new IsotropicRadiationCNES95Convention(500.0, 0.7, 0.7));
747
748 FNP.addForceModel(forceModel);
749
750
751 FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1000.);
752 FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
753 SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
754 FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
755 PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
756
757 Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getX() - finPVC_R.getPosition().getX()) < FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
758 Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getY() - finPVC_R.getPosition().getY()) < FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
759 Assertions.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getZ() - finPVC_R.getPosition().getZ()) < FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
760 }
761
762 @Test
763 void testFlatteningLEO() {
764
765 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", false));
766 NormalizedSphericalHarmonicsProvider gravityField = GravityFieldFactory.getNormalizedProvider(20, 20);
767
768 final AbsoluteDate initialDate = new AbsoluteDate("2000-02-14T00:00:00.000Z", TimeScalesFactory.getUTC());
769 final Orbit initialOrbit =
770 new CircularOrbit(6992986.636088, -5.3589198032e-04, 1.2502577102e-03,
771 FastMath.toRadians(97.827521740), FastMath.toRadians(124.587089845),
772 FastMath.toRadians(0.0), PositionAngleType.MEAN,
773 FramesFactory.getMOD(IERSConventions.IERS_2010),
774 initialDate, gravityField.getMu());
775 final SpacecraftState initialState = new SpacecraftState(initialOrbit, 30.0);
776
777 final CelestialBody sun = CelestialBodyFactory.getSun();
778 final CelestialBody moon = CelestialBodyFactory.getMoon();
779 final RadiationSensitive spacecraft = new IsotropicRadiationClassicalConvention(0.12, 0.5, 0.5);
780 final OneAxisEllipsoid oblateEarth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
781 Constants.WGS84_EARTH_FLATTENING,
782 FramesFactory.getITRF(IERSConventions.IERS_2010, false));
783 final OneAxisEllipsoid sphericalEarth = new OneAxisEllipsoid(oblateEarth.getEquatorialRadius(),
784 0.0,
785 oblateEarth.getBodyFrame());
786
787 PropagatorsParallelizer parallelizer =
788 new PropagatorsParallelizer(Arrays.asList(createLeoPropagator(sun, moon, oblateEarth,
789 gravityField, spacecraft, initialState),
790 createLeoPropagator(sun, moon, sphericalEarth,
791 gravityField, spacecraft, initialState)),
792 new DistanceCheckerHandler(10.0, 0.816805));
793 parallelizer.propagate(initialDate, initialDate.shiftedBy(7 * Constants.JULIAN_DAY));
794
795 }
796
797 private NumericalPropagator createLeoPropagator(final CelestialBody sun, final CelestialBody moon, final OneAxisEllipsoid earth,
798 final NormalizedSphericalHarmonicsProvider gravityField,
799 final RadiationSensitive spacecraft, SpacecraftState initialState) {
800 final double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1e-6).getTolerances(initialState.getOrbit(), OrbitType.CIRCULAR);
801 final NumericalPropagator propagator =
802 new NumericalPropagator(new DormandPrince853Integrator(1.0e-9, 60.0, tol[0], tol[1]));
803 propagator.setOrbitType(OrbitType.CIRCULAR);
804 propagator.setPositionAngleType(PositionAngleType.TRUE);
805 propagator.addForceModel(new HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityField));
806 propagator.addForceModel(new ThirdBodyAttraction(sun));
807 propagator.addForceModel(new ThirdBodyAttraction(moon));
808 propagator.addForceModel(new SolarRadiationPressure(sun, earth, spacecraft));
809 propagator.setInitialState(initialState);
810 return propagator;
811 }
812
813 private static class DistanceCheckerHandler implements MultiSatStepHandler {
814
815 private final double samplingStep;
816 private final double expectedMax;
817 private double currentMax;
818 private AbsoluteDate nextSample;
819
820 DistanceCheckerHandler(final double samplingStep, final double expectedMax) {
821 this.samplingStep = samplingStep;
822 this.expectedMax = expectedMax;
823 }
824
825
826 @Override
827 public void init(final List<SpacecraftState> states0, final AbsoluteDate t) {
828 currentMax = 0.0;
829 nextSample = states0.get(0).getDate().shiftedBy(samplingStep);
830 }
831
832 public void handleStep(List<OrekitStepInterpolator> interpolators) {
833 while (interpolators.get(0).getPreviousState().getDate().isBefore(nextSample) &&
834 interpolators.get(0).getCurrentState().getDate().isAfterOrEqualTo(nextSample)) {
835 final SpacecraftState state0 = interpolators.get(0).getInterpolatedState(nextSample);
836 final SpacecraftState state1 = interpolators.get(1).getInterpolatedState(nextSample);
837 final double distance = Vector3D.distance(state0.getPosition(),
838 state1.getPosition());
839 currentMax = FastMath.max(currentMax, distance);
840 nextSample = nextSample.shiftedBy(samplingStep);
841 }
842 }
843
844 public void finish(final List<SpacecraftState> finalStates) {
845 Assertions.assertEquals(expectedMax, currentMax, 1.0e-3 * expectedMax);
846 }
847
848 }
849
850 @Test
851 void testDependsOnlyOnPosition() {
852
853 final IsotropicRadiationSingleCoefficient mockedSpacecraft = Mockito.mock(IsotropicRadiationSingleCoefficient.class);
854 final SolarRadiationPressure radiationPressure = new SolarRadiationPressure(null, null, mockedSpacecraft);
855
856 final boolean actualValue = radiationPressure.dependsOnPositionOnly();
857
858 Assertions.assertEquals(mockedSpacecraft, radiationPressure.getRadiationSensitiveSpacecraft());
859 Assertions.assertTrue(actualValue);
860 final BoxAndSolarArraySpacecraft mockedBoxSpacecraft = Mockito.mock(BoxAndSolarArraySpacecraft.class);
861 final SolarRadiationPressure boxRadiationPressure = new SolarRadiationPressure(null, null, mockedBoxSpacecraft);
862 Assertions.assertFalse(boxRadiationPressure.dependsOnPositionOnly());
863 }
864
865
866
867
868
869
870 @Test
871 void testMoonPenumbra() {
872 final AbsoluteDate date = new AbsoluteDate(2007, 1, 19, 5, 0, 0, TimeScalesFactory.getGPS());
873 final Vector3D p = new Vector3D(12538484.957505366, 15515522.98001655, -17884023.51839292);
874 final Vector3D v = new Vector3D(-3366.9009055533616, 769.5389825219049, -1679.3840677789601);
875 final Orbit orbit = new CartesianOrbit(new TimeStampedPVCoordinates(date, p, v),
876 FramesFactory.getGCRF(), Constants.EIGEN5C_EARTH_MU);
877 doTestMoonEarth(orbit, 7200.0, 1.0, 0, 0, 0, 1513, 0.0);
878 }
879
880 @Test
881 void testEarthPenumbraOnly() {
882 final AbsoluteDate date = new AbsoluteDate(2007, 3, 13, 17, 14, 0, TimeScalesFactory.getGPS());
883 final Vector3D p = new Vector3D(-26168647.4977583, -1516554.3304749255, -3206794.210706205);
884 final Vector3D v = new Vector3D(-213.65557094060222, -2377.3633988328584, 3079.4740070013495);
885 final Orbit orbit = new CartesianOrbit(new TimeStampedPVCoordinates(date, p, v),
886 FramesFactory.getGCRF(), Constants.EIGEN5C_EARTH_MU);
887 doTestMoonEarth(orbit, 720.0, 1.0, 0, 525, 0, 0, 2.192e-3);
888 }
889
890 @Test
891 void testEarthPenumbraAndUmbra() {
892 final AbsoluteDate date = new AbsoluteDate(2007, 3, 14, 5, 8, 0, TimeScalesFactory.getGPS());
893 final Vector3D p = new Vector3D(-26101379.998276696, -947280.678355501, -3940992.754483608);
894 final Vector3D v = new Vector3D(-348.8911736753223, -2383.738528546711, 3060.9815784341567);
895 final Orbit orbit = new CartesianOrbit(new TimeStampedPVCoordinates(date, p, v),
896 FramesFactory.getGCRF(), Constants.EIGEN5C_EARTH_MU);
897 doTestMoonEarth(orbit, 3600.0, 1.0, 534, 1003, 0, 0, 11.689e-3);
898 }
899
900 private void doTestMoonEarth(Orbit orbit, double duration, double step,
901 int expectedEarthUmbraSteps, int expectedEarthPenumbraSteps,
902 int expectedMoonUmbraSteps, int expectedMoonPenumbraSteps,
903 double expectedDistance) {
904
905 Utils.setDataRoot("2007");
906 final ExtendedPositionProvider sun = CelestialBodyFactory.getSun();
907 final ExtendedPositionProvider moon = CelestialBodyFactory.getMoon();
908 final SpacecraftState initialState = new SpacecraftState(orbit);
909
910
911 SolarRadiationPressure srpWithFlattening =
912 new SolarRadiationPressure(sun,
913 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
914 Constants.WGS84_EARTH_FLATTENING,
915 FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
916 new IsotropicRadiationClassicalConvention(50.0, 0.5, 0.5));
917 srpWithFlattening.addOccultingBody(moon, Constants.MOON_EQUATORIAL_RADIUS);
918 SolarRadiationPressure srpWithoutFlattening =
919 new SolarRadiationPressure(sun,
920 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
921 0.0,
922 FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
923 new IsotropicRadiationClassicalConvention(50.0, 0.5, 0.5));
924 srpWithoutFlattening.addOccultingBody(moon, Constants.MOON_EQUATORIAL_RADIUS);
925
926
927 final OrbitType type = OrbitType.CARTESIAN;
928 double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1e-6).getTolerances(orbit, type);
929
930 final NumericalPropagator propagatorWithFlattening =
931 new NumericalPropagator(new DormandPrince853Integrator(1.0e-9, 300, tol[0], tol[1]));
932 propagatorWithFlattening.setOrbitType(type);
933 propagatorWithFlattening.setInitialState(initialState);
934 propagatorWithFlattening.addForceModel(srpWithFlattening);
935 MoonEclipseStepHandler handler = new MoonEclipseStepHandler(moon, sun, srpWithFlattening);
936 propagatorWithFlattening.setStepHandler(step, handler);
937 final SpacecraftState withFlattening = propagatorWithFlattening.propagate(orbit.getDate().shiftedBy(duration));
938 Assertions.assertEquals(expectedEarthUmbraSteps, handler.earthUmbraSteps);
939 Assertions.assertEquals(expectedEarthPenumbraSteps, handler.earthPenumbraSteps);
940 Assertions.assertEquals(expectedMoonUmbraSteps, handler.moonUmbraSteps);
941 Assertions.assertEquals(expectedMoonPenumbraSteps, handler.moonPenumbraSteps);
942
943 final NumericalPropagator propagatorWithoutFlattening =
944 new NumericalPropagator(new DormandPrince853Integrator(1.0e-9, 300, tol[0], tol[1]));
945 propagatorWithoutFlattening.setOrbitType(type);
946 propagatorWithoutFlattening.setInitialState(initialState);
947 propagatorWithoutFlattening.addForceModel(srpWithoutFlattening);
948 final SpacecraftState withoutFlattening = propagatorWithoutFlattening.propagate(orbit.getDate().shiftedBy(duration));
949
950 Assertions.assertEquals(expectedDistance,
951 Vector3D.distance(withFlattening.getPosition(),
952 withoutFlattening.getPosition()),
953 1.0e-6);
954
955 }
956
957
958
959
960
961 private static class MoonEclipseStepHandler implements OrekitFixedStepHandler {
962
963 final ExtendedPositionProvider moon;
964 final ExtendedPositionProvider sun;
965 final SolarRadiationPressure srp;
966 int earthUmbraSteps;
967 int earthPenumbraSteps;
968 int moonUmbraSteps;
969 int moonPenumbraSteps;
970
971
972
973 MoonEclipseStepHandler(final ExtendedPositionProvider moon, final ExtendedPositionProvider sun,
974 final SolarRadiationPressure srp) {
975 this.moon = moon;
976 this.sun = sun;
977 this.srp = srp;
978
979 }
980
981
982 @Override
983 public void init(final SpacecraftState s0, final AbsoluteDate t, final double step) {
984 earthUmbraSteps = 0;
985 earthPenumbraSteps = 0;
986 moonUmbraSteps = 0;
987 moonPenumbraSteps = 0;
988 }
989
990
991 @Override
992 public void handleStep(final SpacecraftState currentState) {
993 final AbsoluteDate date = currentState.getDate();
994 final Frame frame = currentState.getFrame();
995
996 final Vector3D moonPos = moon.getPosition(date, frame);
997 final Vector3D sunPos = sun.getPosition(date, frame);
998 final Vector3D statePos = currentState.getPosition();
999
1000
1001 final double[] moonAngles = srp.getGeneralEclipseAngles(statePos, moonPos, Constants.MOON_EQUATORIAL_RADIUS,
1002 sunPos, Constants.SUN_RADIUS);
1003 final double moonUmbra = moonAngles[0] - moonAngles[1] + moonAngles[2];
1004 final boolean isInMoonUmbra = (moonUmbra < 1.0e-10);
1005 if (isInMoonUmbra) {
1006 ++moonUmbraSteps;
1007 }
1008
1009 final double moonPenumbra = moonAngles[0] - moonAngles[1] - moonAngles[2];
1010 final boolean isInMoonPenumbra = (moonPenumbra < -1.0e-10);
1011 if (isInMoonPenumbra) {
1012 ++moonPenumbraSteps;
1013 }
1014
1015
1016 final OccultationEngine.OccultationAngles angles = srp.getOccultingBodies().get(0).angles(currentState);
1017
1018 final double earthUmbra = angles.getSeparation() - angles.getLimbRadius() + angles.getOccultedApparentRadius();
1019 final boolean isInEarthUmbra = (earthUmbra < 1.0e-10);
1020 if (isInEarthUmbra) {
1021 ++earthUmbraSteps;
1022 }
1023
1024 final double earthPenumbra = angles.getSeparation() - angles.getLimbRadius() - angles.getOccultedApparentRadius();
1025 final boolean isInEarthPenumbra = (earthPenumbra < -1.0e-10);
1026 if (isInEarthPenumbra) {
1027 ++earthPenumbraSteps;
1028 }
1029
1030
1031
1032 final double lightingRatio = srp.getLightingRatio(currentState);
1033
1034
1035 if (isInMoonUmbra || isInEarthUmbra) {
1036 Assertions.assertEquals(0.0, lightingRatio, 1e-8);
1037 }
1038 else if (isInMoonPenumbra || isInEarthPenumbra) {
1039 Assertions.assertTrue(lightingRatio < 1.0);
1040 Assertions.assertTrue(lightingRatio > 0.0);
1041 }
1042 else {
1043 Assertions.assertEquals(1.0, lightingRatio, 1e-8);
1044 }
1045 }
1046 }
1047
1048
1049
1050
1051
1052
1053 @Test
1054 void testFieldMoonPenumbra() {
1055 Field<Binary64> field = Binary64Field.getInstance();
1056 final FieldAbsoluteDate<Binary64> date = new FieldAbsoluteDate<>(field, 2007, 1, 19, 5, 0, 0, TimeScalesFactory.getGPS());
1057 final FieldVector3D<Binary64> p = new FieldVector3D<>(field, new Vector3D(12538484.957505366, 15515522.98001655, -17884023.51839292));
1058 final FieldVector3D<Binary64> v = new FieldVector3D<>(field, new Vector3D(-3366.9009055533616, 769.5389825219049, -1679.3840677789601));
1059 final FieldVector3D<Binary64> a = FieldVector3D.getZero(field);
1060 final FieldOrbit<Binary64> orbit = new FieldCartesianOrbit<>(new TimeStampedFieldPVCoordinates<>(date, p, v, a),
1061 FramesFactory.getGCRF(),
1062 field.getZero().newInstance(Constants.EIGEN5C_EARTH_MU));
1063 doTestFieldMoonEarth(orbit, field.getZero().newInstance(7200.0), field.getZero().newInstance(1.0), 0, 0, 0, 1513, 0.0);
1064 }
1065
1066 @Test
1067 void testFieldEarthPenumbraOnly() {
1068 Field<Binary64> field = Binary64Field.getInstance();
1069 final FieldAbsoluteDate<Binary64> date = new FieldAbsoluteDate<>(field, 2007, 3, 13, 17, 14, 0, TimeScalesFactory.getGPS());
1070 final FieldVector3D<Binary64> p = new FieldVector3D<>(field, new Vector3D(-26168647.4977583, -1516554.3304749255, -3206794.210706205));
1071 final FieldVector3D<Binary64> v = new FieldVector3D<>(field, new Vector3D(-213.65557094060222, -2377.3633988328584, 3079.4740070013495));
1072 final FieldVector3D<Binary64> a = FieldVector3D.getZero(field);
1073 final FieldOrbit<Binary64> orbit = new FieldCartesianOrbit<>(new TimeStampedFieldPVCoordinates<>(date, p, v, a),
1074 FramesFactory.getGCRF(),
1075 field.getZero().newInstance(Constants.EIGEN5C_EARTH_MU));
1076 doTestFieldMoonEarth(orbit, field.getZero().newInstance(720.0), field.getZero().newInstance(1.0), 0, 525, 0, 0, 2.192e-3);
1077 }
1078
1079 @Test
1080 void testFieldEarthPenumbraAndUmbra() {
1081 Field<Binary64> field = Binary64Field.getInstance();
1082 final FieldAbsoluteDate<Binary64> date = new FieldAbsoluteDate<>(field, 2007, 3, 14, 5, 8, 0, TimeScalesFactory.getGPS());
1083 final FieldVector3D<Binary64> p = new FieldVector3D<>(field, new Vector3D(-26101379.998276696, -947280.678355501, -3940992.754483608));
1084 final FieldVector3D<Binary64> v = new FieldVector3D<>(field, new Vector3D(-348.8911736753223, -2383.738528546711, 3060.9815784341567));
1085 final FieldVector3D<Binary64> a = FieldVector3D.getZero(field);
1086 final FieldOrbit<Binary64> orbit = new FieldCartesianOrbit<>(new TimeStampedFieldPVCoordinates<>(date, p, v, a),
1087 FramesFactory.getGCRF(),
1088 field.getZero().newInstance(Constants.EIGEN5C_EARTH_MU));
1089 doTestFieldMoonEarth(orbit, field.getZero().newInstance(3600.0), field.getZero().newInstance(1.0), 534, 1003, 0, 0, 11.611e-3);
1090 }
1091
1092 private <T extends CalculusFieldElement<T>> void doTestFieldMoonEarth(FieldOrbit<T> orbit, T duration, T step,
1093 int expectedEarthUmbraSteps, int expectedEarthPenumbraSteps,
1094 int expectedMoonUmbraSteps, int expectedMoonPenumbraSteps,
1095 double expectedDistance) {
1096
1097 Utils.setDataRoot("2007");
1098 final Field<T> field = orbit.getDate().getField();
1099 final ExtendedPositionProvider sun = CelestialBodyFactory.getSun();
1100 final ExtendedPositionProvider moon = CelestialBodyFactory.getMoon();
1101 final FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(orbit);
1102
1103
1104 SolarRadiationPressure srpWithFlattening =
1105 new SolarRadiationPressure(sun,
1106 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
1107 Constants.WGS84_EARTH_FLATTENING,
1108 FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
1109 new IsotropicRadiationClassicalConvention(50.0, 0.5, 0.5));
1110 srpWithFlattening.addOccultingBody(moon, Constants.MOON_EQUATORIAL_RADIUS);
1111 SolarRadiationPressure srpWithoutFlattening =
1112 new SolarRadiationPressure(sun,
1113 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
1114 0.0,
1115 FramesFactory.getITRF(IERSConventions.IERS_2010, true)),
1116 new IsotropicRadiationClassicalConvention(50.0, 0.5, 0.5));
1117 srpWithoutFlattening.addOccultingBody(moon, Constants.MOON_EQUATORIAL_RADIUS);
1118
1119
1120 final OrbitType type = OrbitType.CARTESIAN;
1121 double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1e-6).getTolerances(orbit, type);
1122 final FieldNumericalPropagator<T> propagatorWithFlattening =
1123 new FieldNumericalPropagator<>(field,
1124 new DormandPrince853FieldIntegrator<>(field, 1.0e-9, 300, tol[0], tol[1]));
1125 propagatorWithFlattening.setOrbitType(type);
1126 propagatorWithFlattening.setInitialState(initialState);
1127 propagatorWithFlattening.addForceModel(srpWithFlattening);
1128 FieldMoonEclipseStepHandler<T> handler = new FieldMoonEclipseStepHandler<>(moon, sun, srpWithFlattening);
1129 propagatorWithFlattening.setStepHandler(step, handler);
1130 final FieldSpacecraftState<T> withFlattening = propagatorWithFlattening.propagate(orbit.getDate().shiftedBy(duration));
1131 Assertions.assertEquals(expectedEarthUmbraSteps, handler.earthUmbraSteps);
1132 Assertions.assertEquals(expectedEarthPenumbraSteps, handler.earthPenumbraSteps);
1133 Assertions.assertEquals(expectedMoonUmbraSteps, handler.moonUmbraSteps);
1134 Assertions.assertEquals(expectedMoonPenumbraSteps, handler.moonPenumbraSteps);
1135
1136 final FieldNumericalPropagator<T> propagatorWithoutFlattening =
1137 new FieldNumericalPropagator<>(field,
1138 new DormandPrince853FieldIntegrator<>(field, 1.0e-9, 300, tol[0], tol[1]));
1139 propagatorWithoutFlattening.setOrbitType(type);
1140 propagatorWithoutFlattening.setInitialState(initialState);
1141 propagatorWithoutFlattening.addForceModel(srpWithoutFlattening);
1142 final FieldSpacecraftState<T> withoutFlattening = propagatorWithoutFlattening.propagate(orbit.getDate().shiftedBy(duration));
1143
1144 Assertions.assertEquals(expectedDistance,
1145 FieldVector3D.distance(withFlattening.getPosition(),
1146 withoutFlattening.getPosition()).getReal(),
1147 1.0e-6);
1148
1149 }
1150
1151
1152
1153
1154
1155 private static class FieldMoonEclipseStepHandler<T extends CalculusFieldElement<T>> implements FieldOrekitFixedStepHandler<T> {
1156
1157 final ExtendedPositionProvider moon;
1158 final ExtendedPositionProvider sun;
1159 final SolarRadiationPressure srp;
1160 int earthUmbraSteps;
1161 int earthPenumbraSteps;
1162 int moonUmbraSteps;
1163 int moonPenumbraSteps;
1164
1165
1166
1167 FieldMoonEclipseStepHandler(final ExtendedPositionProvider moon, final ExtendedPositionProvider sun,
1168 final SolarRadiationPressure srp) {
1169 this.moon = moon;
1170 this.sun = sun;
1171 this.srp = srp;
1172
1173 }
1174
1175
1176 @Override
1177 public void init(final FieldSpacecraftState<T> s0, final FieldAbsoluteDate<T> t, final T step) {
1178 earthUmbraSteps = 0;
1179 earthPenumbraSteps = 0;
1180 moonUmbraSteps = 0;
1181 moonPenumbraSteps = 0;
1182 }
1183
1184
1185 @Override
1186 public void handleStep(final FieldSpacecraftState<T> currentState) {
1187 final FieldAbsoluteDate<T> date = currentState.getDate();
1188 final Frame frame = currentState.getFrame();
1189
1190 final FieldVector3D<T> moonPos = moon.getPosition(date, frame);
1191 final FieldVector3D<T> sunPos = sun.getPosition(date, frame);
1192 final FieldVector3D<T> statePos = currentState.getPosition();
1193
1194
1195 final T[] moonAngles = srp.getGeneralEclipseAngles(statePos, moonPos,
1196 date.getField().getZero().newInstance(Constants.MOON_EQUATORIAL_RADIUS),
1197 sunPos,
1198 date.getField().getZero().newInstance(Constants.SUN_RADIUS));
1199 final T moonUmbra = moonAngles[0].subtract(moonAngles[1]).add(moonAngles[2]);
1200 final boolean isInMoonUmbra = (moonUmbra.getReal() < 1.0e-10);
1201 if (isInMoonUmbra) {
1202 ++moonUmbraSteps;
1203 }
1204
1205 final T moonPenumbra = moonAngles[0].subtract(moonAngles[1]).subtract(moonAngles[2]);
1206 final boolean isInMoonPenumbra = (moonPenumbra.getReal() < -1.0e-10);
1207 if (isInMoonPenumbra) {
1208 ++moonPenumbraSteps;
1209 }
1210
1211
1212 final OccultationEngine.FieldOccultationAngles<T> angles = srp.getOccultingBodies().get(0).angles(currentState);
1213
1214 final T earthUmbra = angles.getSeparation().subtract(angles.getLimbRadius()).add(angles.getOccultedApparentRadius());
1215 final boolean isInEarthUmbra = (earthUmbra.getReal() < 1.0e-10);
1216 if (isInEarthUmbra) {
1217 ++earthUmbraSteps;
1218 }
1219
1220 final T earthPenumbra = angles.getSeparation().subtract(angles.getLimbRadius()).subtract(angles.getOccultedApparentRadius());
1221 final boolean isInEarthPenumbra = (earthPenumbra.getReal() < -1.0e-10);
1222 if (isInEarthPenumbra) {
1223 ++earthPenumbraSteps;
1224 }
1225
1226
1227 final T lightingRatio = srp.getLightingRatio(currentState);
1228
1229
1230 if (isInMoonUmbra || isInEarthUmbra) {
1231 Assertions.assertEquals(0.0, lightingRatio.getReal(), 1e-8);
1232 }
1233 else if (isInMoonPenumbra || isInEarthPenumbra) {
1234 Assertions.assertTrue(lightingRatio.getReal() < 1.0);
1235 Assertions.assertTrue(lightingRatio.getReal() > 0.0);
1236 }
1237 else {
1238 Assertions.assertEquals(1.0, lightingRatio.getReal(), 1e-8);
1239 }
1240 }
1241 }
1242
1243 @BeforeEach
1244 public void setUp() {
1245 Utils.setDataRoot("regular-data:potential");
1246 }
1247 }