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