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