1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.events;
18
19 import org.hipparchus.Field;
20 import org.hipparchus.CalculusFieldElement;
21 import org.hipparchus.exception.MathRuntimeException;
22 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
23 import org.hipparchus.geometry.euclidean.threed.Vector3D;
24 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
25 import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
26 import org.hipparchus.util.Decimal64Field;
27 import org.hipparchus.util.FastMath;
28 import org.junit.Assert;
29 import org.junit.Before;
30 import org.junit.Test;
31 import org.orekit.Utils;
32 import org.orekit.bodies.CelestialBody;
33 import org.orekit.bodies.CelestialBodyFactory;
34 import org.orekit.errors.OrekitException;
35 import org.orekit.frames.FramesFactory;
36 import org.orekit.orbits.FieldCartesianOrbit;
37 import org.orekit.orbits.FieldEquinoctialOrbit;
38 import org.orekit.orbits.FieldOrbit;
39 import org.orekit.orbits.OrbitType;
40 import org.orekit.propagation.FieldSpacecraftState;
41 import org.orekit.propagation.analytical.FieldKeplerianPropagator;
42 import org.orekit.propagation.events.handlers.FieldStopOnDecreasing;
43 import org.orekit.propagation.numerical.FieldNumericalPropagator;
44 import org.orekit.time.AbsoluteDate;
45 import org.orekit.time.FieldAbsoluteDate;
46 import org.orekit.time.TimeScalesFactory;
47 import org.orekit.utils.FieldPVCoordinates;
48 import org.orekit.utils.TimeStampedFieldPVCoordinates;
49
50 public class FieldEclipseDetectorTest {
51
52 private double mu;
53 private CelestialBody sun;
54 private CelestialBody earth;
55 private double sunRadius;
56 private double earthRadius;
57
58
59 @Before
60 public void setUp() {
61 try {
62 Utils.setDataRoot("regular-data");
63 sun = CelestialBodyFactory.getSun();
64 earth = CelestialBodyFactory.getEarth();
65 sunRadius = 696000000.;
66 earthRadius = 6400000.;
67 mu = 3.9860047e14;
68 } catch (OrekitException oe) {
69 Assert.fail(oe.getLocalizedMessage());
70 }
71 }
72
73 @Test
74 public void testEclipse() {
75 doTestEclipse(Decimal64Field.getInstance());
76 }
77 @Test
78 public void testPenumbra() {
79 doTestPenumbra(Decimal64Field.getInstance());
80 }
81 @Test
82 public void testWithMethods() {
83 doTestWithMethods(Decimal64Field.getInstance());
84 }
85
86 @Test
87 public void testInsideOcculting() {
88 doTestInsideOcculting(Decimal64Field.getInstance());
89 }
90 @Test
91 public void testInsideOcculted() {
92 doTestInsideOcculted(Decimal64Field.getInstance());
93 }
94 @Test
95 public void testTooSmallMaxIterationCount() {
96 testTooSmallMaxIterationCount(Decimal64Field.getInstance());
97 }
98
99
100
101 private <T extends CalculusFieldElement<T>> void doTestEclipse(Field<T> field) {
102 T zero = field.getZero();
103 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.560), zero.add(-25767.25680));
104 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(505.8479685), zero.add(942.7809215), zero.add(7435.922231));
105 FieldAbsoluteDate<T> iniDate = new FieldAbsoluteDate<>(field, 1969, 7, 28, 4, 0, 0.0, TimeScalesFactory.getTT());
106 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
107 FramesFactory.getGCRF(), iniDate, zero.add(mu));
108 FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(orbit);
109 FieldKeplerianPropagator<T> propagator = new FieldKeplerianPropagator<>(orbit);
110 propagator.resetInitialState(initialState);
111
112 FieldEclipseDetector<T> e = new FieldEclipseDetector<>(field.getZero().add(60.), field.getZero().add(1e-3),
113 sun, sunRadius,
114 earth, earthRadius).
115 withHandler(new FieldStopOnDecreasing<FieldEclipseDetector<T>, T>()).
116 withUmbra();
117 Assert.assertEquals(60.0, e.getMaxCheckInterval().getReal(), 1.0e-15);
118 Assert.assertEquals(1.0e-3, e.getThreshold().getReal(), 1.0e-15);
119 Assert.assertEquals(AbstractDetector.DEFAULT_MAX_ITER, e.getMaxIterationCount());
120 Assert.assertSame(sun, e.getOcculted());
121 Assert.assertEquals(sunRadius, e.getOccultedRadius(), 1.0);
122 Assert.assertSame(earth, e.getOcculting());
123 Assert.assertEquals(earthRadius, e.getOccultingRadius(), 1.0);
124 Assert.assertTrue(e.getTotalEclipse());
125 propagator.addEventDetector(e);
126 final FieldSpacecraftState<T> finalState = propagator.propagate(iniDate.shiftedBy(6000));
127 Assert.assertEquals(2303.1835, finalState.getDate().durationFrom(iniDate).getReal(), 1.0e-3);
128 }
129
130 private <T extends CalculusFieldElement<T>> void doTestPenumbra(Field<T> field) {
131 T zero = field.getZero();
132 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.560), zero.add(-25767.25680));
133 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(505.8479685), zero.add(942.7809215), zero.add(7435.922231));
134 FieldAbsoluteDate<T> iniDate = new FieldAbsoluteDate<>(field, 1969, 7, 28, 4, 0, 0.0, TimeScalesFactory.getTT());
135 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
136 FramesFactory.getGCRF(), iniDate, zero.add(mu));
137 FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(orbit);
138 double[] absTolerance = {
139 0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001
140 };
141 double[] relTolerance = {
142 1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7
143 };
144 AdaptiveStepsizeFieldIntegrator<T> integrator =
145 new DormandPrince853FieldIntegrator<>(field, 0.001, 1000, absTolerance, relTolerance);
146 integrator.setInitialStepSize(60.);
147 FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
148 propagator.setOrbitType(OrbitType.EQUINOCTIAL);
149 propagator.setInitialState(initialState);
150 sun = CelestialBodyFactory.getSun();
151 earth = CelestialBodyFactory.getEarth();
152 sunRadius = 696000000.;
153 earthRadius = 6400000.;
154
155 FieldEclipseDetector<T> e = new FieldEclipseDetector<>(zero.add(60.), zero.add(1.e-3), sun, sunRadius,
156 earth, earthRadius).
157 withPenumbra();
158 Assert.assertFalse(e.getTotalEclipse());
159 propagator.addEventDetector(e);
160 final FieldSpacecraftState<T> finalState = propagator.propagate(iniDate.shiftedBy(6000));
161 Assert.assertEquals(4388.155852, finalState.getDate().durationFrom(iniDate).getReal(), 2.0e-6);
162 }
163
164 private <T extends CalculusFieldElement<T>> void doTestWithMethods(Field<T> field) {
165 T zero = field.getZero();
166 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.560), zero.add(-25767.25680));
167 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(505.8479685), zero.add(942.7809215), zero.add(7435.922231));
168 FieldAbsoluteDate<T> iniDate = new FieldAbsoluteDate<>(field, 1969, 7, 28, 4, 0, 0.0, TimeScalesFactory.getTT());
169 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
170 FramesFactory.getGCRF(), iniDate, zero.add(mu));
171 FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(orbit);
172 double[] absTolerance = {
173 0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001
174 };
175 double[] relTolerance = {
176 1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7
177 };
178 AdaptiveStepsizeFieldIntegrator<T> integrator =
179 new DormandPrince853FieldIntegrator<>(field, 0.001, 1000, absTolerance, relTolerance);
180 integrator.setInitialStepSize(60);
181 FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
182 propagator.setOrbitType(OrbitType.EQUINOCTIAL);
183 propagator.setInitialState(initialState);
184 sun = CelestialBodyFactory.getSun();
185 earth = CelestialBodyFactory.getEarth();
186 sunRadius = 696000000.;
187 earthRadius = 6400000.;
188
189 FieldEclipseDetector<T> e = new FieldEclipseDetector<>(field.getZero().add(60.), field.getZero().add(1e-3),
190 sun, sunRadius,
191 earth, earthRadius).
192 withHandler(new FieldStopOnDecreasing<FieldEclipseDetector<T>, T>()).
193 withMaxCheck(field.getZero().add(120.0)).
194 withThreshold(field.getZero().add(1.0e-4)).
195 withMaxIter(12);
196 Assert.assertEquals(120.0, e.getMaxCheckInterval().getReal(), 1.0e-15);
197 Assert.assertEquals(1.0e-4, e.getThreshold().getReal(), 1.0e-15);
198 Assert.assertEquals(12, e.getMaxIterationCount());
199 propagator.addEventDetector(e);
200 final FieldSpacecraftState<T> finalState = propagator.propagate(iniDate.shiftedBy(6000));
201 Assert.assertEquals(2303.1835, finalState.getDate().durationFrom(iniDate).getReal(), 1.0e-3);
202
203 }
204
205 private <T extends CalculusFieldElement<T>> void doTestInsideOcculting(Field<T> field) {
206 T zero = field.getZero();
207 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.560), zero.add(-25767.25680));
208 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(505.8479685), zero.add(942.7809215), zero.add(7435.922231));
209 FieldAbsoluteDate<T> iniDate = new FieldAbsoluteDate<>(field, 1969, 7, 28, 4, 0, 0.0, TimeScalesFactory.getTT());
210 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
211 FramesFactory.getGCRF(), iniDate, zero.add(mu));
212 FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(orbit);
213 double[] absTolerance = {
214 0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001
215 };
216 double[] relTolerance = {
217 1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7
218 };
219 AdaptiveStepsizeFieldIntegrator<T> integrator =
220 new DormandPrince853FieldIntegrator<>(field, 0.001, 1000, absTolerance, relTolerance);
221 integrator.setInitialStepSize(60);
222 FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
223 propagator.setOrbitType(OrbitType.EQUINOCTIAL);
224 propagator.setInitialState(initialState);
225 sun = CelestialBodyFactory.getSun();
226 earth = CelestialBodyFactory.getEarth();
227 sunRadius = 696000000.;
228 earthRadius = 6400000.;
229
230 FieldEclipseDetector<T> e = new FieldEclipseDetector<>(field.getZero().add(60.), field.getZero().add(1.e-3),
231 sun, sunRadius,
232 earth, earthRadius);
233 FieldSpacecraftState<T> s = new FieldSpacecraftState<>(new FieldCartesianOrbit<>(new TimeStampedFieldPVCoordinates<>(FieldAbsoluteDate.getJ2000Epoch(field),
234 new FieldPVCoordinates<>(new FieldVector3D<>(field.getZero().add(1e6),
235 field.getZero().add(2e6),
236 field.getZero().add(3e6)),
237 new FieldVector3D<>(field.getZero().add(1000),
238 field.getZero().add(0),
239 field.getZero().add(0)))),
240 FramesFactory.getGCRF(),
241 zero.add(mu)));
242 Assert.assertEquals(-FastMath.PI, e.g(s).getReal(), 1.0e-15);
243 }
244
245 private <T extends CalculusFieldElement<T>> void doTestInsideOcculted(Field<T> field) {
246 T zero = field.getZero();
247 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.560), zero.add(-25767.25680));
248 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(505.8479685), zero.add(942.7809215), zero.add(7435.922231));
249 FieldAbsoluteDate<T> iniDate = new FieldAbsoluteDate<>(field, 1969, 7, 28, 4, 0, 0.0, TimeScalesFactory.getTT());
250 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
251 FramesFactory.getGCRF(), iniDate, zero.add(mu));
252 FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(orbit);
253 double[] absTolerance = {
254 0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001
255 };
256 double[] relTolerance = {
257 1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7
258 };
259 AdaptiveStepsizeFieldIntegrator<T> integrator =
260 new DormandPrince853FieldIntegrator<>(field, 0.001, 1000, absTolerance, relTolerance);
261 integrator.setInitialStepSize(60);
262 FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
263 propagator.setOrbitType(OrbitType.EQUINOCTIAL);
264 propagator.setInitialState(initialState);
265 sun = CelestialBodyFactory.getSun();
266 earth = CelestialBodyFactory.getEarth();
267 sunRadius = 696000000.;
268 earthRadius = 6400000.;
269
270 FieldEclipseDetector<T> e = new FieldEclipseDetector<>(field.getZero().add(60.), field.getZero().add(1.e-3),
271 sun, sunRadius,
272 earth, earthRadius);
273 Vector3D p = sun.getPVCoordinates(AbsoluteDate.J2000_EPOCH,
274 FramesFactory.getGCRF()).getPosition();
275 FieldSpacecraftState<T> s = new FieldSpacecraftState<>(new FieldCartesianOrbit<>(new TimeStampedFieldPVCoordinates<>(FieldAbsoluteDate.getJ2000Epoch(field),
276 new FieldPVCoordinates<>(new FieldVector3D<>(field.getOne(),
277 field.getZero(),
278 field.getZero()).add(p),
279 new FieldVector3D<>(field.getZero(),
280 field.getZero(),
281 field.getOne()))),
282 FramesFactory.getGCRF(),
283 zero.add(mu)));
284 Assert.assertEquals(FastMath.PI, e.g(s).getReal(), 1.0e-15);
285 }
286
287 private <T extends CalculusFieldElement<T>> void testTooSmallMaxIterationCount(Field<T> field) {
288 T zero = field.getZero();
289 final FieldVector3D<T> position = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.560), zero.add(-25767.25680));
290 final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(505.8479685), zero.add(942.7809215), zero.add(7435.922231));
291 FieldAbsoluteDate<T> iniDate = new FieldAbsoluteDate<>(field, 1969, 7, 28, 4, 0, 0.0, TimeScalesFactory.getTT());
292 final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
293 FramesFactory.getGCRF(), iniDate, zero.add(mu));
294 FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(orbit);
295 double[] absTolerance = {
296 0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001
297 };
298 double[] relTolerance = {
299 1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7
300 };
301 AdaptiveStepsizeFieldIntegrator<T> integrator =
302 new DormandPrince853FieldIntegrator<>(field, 0.001, 1000, absTolerance, relTolerance);
303 integrator.setInitialStepSize(60);
304 FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
305 propagator.setOrbitType(OrbitType.EQUINOCTIAL);
306 propagator.setInitialState(initialState);
307 sun = CelestialBodyFactory.getSun();
308 earth = CelestialBodyFactory.getEarth();
309 sunRadius = 696000000.;
310 earthRadius = 6400000.;
311
312 int n = 5;
313 FieldEclipseDetector<T> e = new FieldEclipseDetector<>(field.getZero().add(60.), field.getZero().add(1.e-3),
314 sun, sunRadius,
315 earth, earthRadius).
316 withHandler(new FieldStopOnDecreasing<FieldEclipseDetector<T>, T>()).
317 withMaxCheck(field.getZero().add(120.0)).
318 withThreshold(field.getZero().add(1.0e-4)).
319 withMaxIter(n);
320 propagator.addEventDetector(e);
321 try {
322 propagator.propagate(iniDate.shiftedBy(6000));
323 Assert.fail("an exception should have been thrown");
324 } catch (OrekitException oe) {
325 Assert.assertEquals(n, ((Integer) ((MathRuntimeException) oe.getCause()).getParts()[0]).intValue());
326 }
327 }
328
329 }
330