1   /* Copyright 2002-2025 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
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.junit.jupiter.params.ParameterizedTest;
34  import org.junit.jupiter.params.provider.ValueSource;
35  import org.mockito.Mockito;
36  import org.orekit.Utils;
37  import org.orekit.bodies.AnalyticalSolarPositionProvider;
38  import org.orekit.bodies.CelestialBody;
39  import org.orekit.bodies.CelestialBodyFactory;
40  import org.orekit.bodies.OneAxisEllipsoid;
41  import org.orekit.errors.OrekitException;
42  import org.orekit.errors.OrekitMessages;
43  import org.orekit.frames.FramesFactory;
44  import org.orekit.models.earth.ReferenceEllipsoid;
45  import org.orekit.orbits.FieldCartesianOrbit;
46  import org.orekit.orbits.FieldEquinoctialOrbit;
47  import org.orekit.orbits.FieldOrbit;
48  import org.orekit.orbits.OrbitType;
49  import org.orekit.propagation.FieldSpacecraftState;
50  import org.orekit.propagation.analytical.FieldKeplerianPropagator;
51  import org.orekit.propagation.events.handlers.FieldStopOnDecreasing;
52  import org.orekit.propagation.numerical.FieldNumericalPropagator;
53  import org.orekit.time.AbsoluteDate;
54  import org.orekit.time.FieldAbsoluteDate;
55  import org.orekit.time.TimeScalesFactory;
56  import org.orekit.utils.Constants;
57  import org.orekit.utils.FieldPVCoordinates;
58  import org.orekit.utils.IERSConventions;
59  import org.orekit.utils.OccultationEngine;
60  import org.orekit.utils.TimeStampedFieldPVCoordinates;
61  
62  public class FieldEclipseDetectorTest {
63  
64      private double               mu;
65      private CelestialBody        sun;
66      private OneAxisEllipsoid     earth;
67      private double               sunRadius;
68  
69  
70      @BeforeEach
71      public void setUp() {
72          try {
73              Utils.setDataRoot("regular-data");
74              sun = CelestialBodyFactory.getSun();
75              earth = new OneAxisEllipsoid(6400000., 0.0, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
76              sunRadius = 696000000.;
77              mu  = 3.9860047e14;
78          } catch (OrekitException oe) {
79              Assertions.fail(oe.getLocalizedMessage());
80          }
81      }
82  
83      @Test
84      public void testEclipse() {
85          doTestEclipse(Binary64Field.getInstance());
86      }
87      @Test
88      public void testPenumbra() {
89          doTestPenumbra(Binary64Field.getInstance());
90      }
91      @Test
92      public void testWithMethods() {
93          doTestWithMethods(Binary64Field.getInstance());
94      }
95  
96      @Test
97      public void testInsideOcculting() {
98          doTestInsideOcculting(Binary64Field.getInstance());
99      }
100     @Test
101     public void testInsideOcculted() {
102         doTestInsideOcculted(Binary64Field.getInstance());
103     }
104     @Test
105     public void testTooSmallMaxIterationCount() {
106         testTooSmallMaxIterationCount(Binary64Field.getInstance());
107     }
108 
109     @ParameterizedTest
110     @ValueSource(booleans = {true, false})
111     void testVersusNonField(final boolean totalEclipse) {
112         // GIVEN
113         final Binary64Field field = Binary64Field.getInstance();
114         final Binary64 zero = field.getZero();
115         final FieldVector3D<Binary64> position = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.560), zero.add(-25767.25680));
116         final FieldVector3D<Binary64> velocity = new FieldVector3D<>(zero.add(505.8479685), zero.add(942.7809215), zero.add(7435.922231));
117         FieldAbsoluteDate<Binary64> date = new FieldAbsoluteDate<>(field, 1969, 7, 28, 4, 0, 0.0, TimeScalesFactory.getTT());
118         final FieldOrbit<Binary64> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
119                 FramesFactory.getGCRF(), date, zero.add(mu));
120         FieldSpacecraftState<Binary64> state = new FieldSpacecraftState<>(orbit);
121         final OccultationEngine engine = new OccultationEngine(new AnalyticalSolarPositionProvider(),
122                 Constants.SUN_RADIUS, ReferenceEllipsoid.getWgs84(FramesFactory.getGTOD(true)));
123         final Binary64 margin = new Binary64(1e-8);
124         FieldEclipseDetector<Binary64> fieldDetector = new FieldEclipseDetector<>(field, engine).withMargin(margin);
125         fieldDetector = totalEclipse ? fieldDetector.withUmbra() : fieldDetector.withPenumbra();
126         // WHEN
127         final Binary64 fieldG = fieldDetector.g(state);
128         // THEN
129         EclipseDetector detector = new EclipseDetector(engine).withMargin(margin.getReal());
130         detector = totalEclipse ? detector.withUmbra() : detector.withPenumbra();
131         final double g = detector.g(state.toSpacecraftState());
132         Assertions.assertEquals(g, fieldG.getReal());
133     }
134 
135     @Test
136     void testGetDetectionSettings() {
137         // GIVEN
138         final Binary64Field field = Binary64Field.getInstance();
139         final OccultationEngine mockedEngine = Mockito.mock(OccultationEngine.class);
140         // WHEN
141         final FieldEclipseDetector<Binary64> fieldEclipseDetector = new FieldEclipseDetector<>(field,
142                 Mockito.mock(OccultationEngine.class));
143         // THEN
144         final EventDetectionSettings expectedSettings = new EclipseDetector(mockedEngine).getDetectionSettings();
145         Assertions.assertEquals(expectedSettings.getMaxIterationCount(), fieldEclipseDetector.getMaxIterationCount());
146         Assertions.assertEquals(expectedSettings.getThreshold(), fieldEclipseDetector.getThreshold().getReal());
147     }
148 
149     private <T extends CalculusFieldElement<T>> void doTestEclipse(Field<T> field) {
150         T zero = field.getZero();
151         final FieldVector3D<T> position  = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.560), zero.add(-25767.25680));
152         final FieldVector3D<T> velocity  = new FieldVector3D<>(zero.add(505.8479685), zero.add(942.7809215), zero.add(7435.922231));
153         FieldAbsoluteDate<T> iniDate = new FieldAbsoluteDate<>(field, 1969, 7, 28, 4, 0, 0.0, TimeScalesFactory.getTT());
154         final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position,  velocity),
155                                                                 FramesFactory.getGCRF(), iniDate, zero.add(mu));
156         FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(orbit);
157         FieldKeplerianPropagator<T> propagator = new FieldKeplerianPropagator<>(orbit);
158         propagator.resetInitialState(initialState);
159 
160         FieldEclipseDetector<T> e = new FieldEclipseDetector<>(field, sun, sunRadius, earth).
161                                     withMaxCheck(60.0).
162                                     withThreshold(zero.newInstance(1e-3)).
163                                     withHandler(new FieldStopOnDecreasing<T>()).
164                                     withUmbra();
165         Assertions.assertEquals(60.0, e.getMaxCheckInterval().currentInterval(null, true), 1.0e-15);
166         Assertions.assertEquals(1.0e-3, e.getThreshold().getReal(), 1.0e-15);
167         Assertions.assertEquals(AbstractDetector.DEFAULT_MAX_ITER, e.getMaxIterationCount());
168         Assertions.assertEquals(0.0, e.getMargin().getReal(), 1.0e-15);
169         Assertions.assertSame(sun, e.getOccultationEngine().getOcculted());
170         Assertions.assertEquals(sunRadius, e.getOccultationEngine().getOccultedRadius(), 1.0);
171         Assertions.assertSame(earth, e.getOccultationEngine().getOcculting());
172         Assertions.assertTrue(e.getTotalEclipse());
173         propagator.addEventDetector(e);
174         final FieldSpacecraftState<T> finalState = propagator.propagate(iniDate.shiftedBy(6000));
175         Assertions.assertEquals(2303.1835, finalState.getDate().durationFrom(iniDate).getReal(), 1.0e-3);
176     }
177 
178     private <T extends CalculusFieldElement<T>> void doTestPenumbra(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(60.0).
201                                     withThreshold(zero.newInstance(1e-3)).
202                                     withPenumbra();
203         Assertions.assertFalse(e.getTotalEclipse());
204         propagator.addEventDetector(e);
205         final FieldSpacecraftState<T> finalState = propagator.propagate(iniDate.shiftedBy(6000));
206         Assertions.assertEquals(4388.155852, finalState.getDate().durationFrom(iniDate).getReal(), 2.0e-6);
207     }
208 
209     private <T extends CalculusFieldElement<T>> void doTestWithMethods(Field<T> field) {
210         T zero = field.getZero();
211         final FieldVector3D<T> position  = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.560), zero.add(-25767.25680));
212         final FieldVector3D<T> velocity  = new FieldVector3D<>(zero.add(505.8479685), zero.add(942.7809215), zero.add(7435.922231));
213         FieldAbsoluteDate<T> iniDate = new FieldAbsoluteDate<>(field, 1969, 7, 28, 4, 0, 0.0, TimeScalesFactory.getTT());
214         final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position,  velocity),
215                                                                 FramesFactory.getGCRF(), iniDate, zero.add(mu));
216         FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(orbit);
217         double[] absTolerance = {
218             0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001
219         };
220         double[] relTolerance = {
221             1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7
222         };
223         AdaptiveStepsizeFieldIntegrator<T> integrator =
224             new DormandPrince853FieldIntegrator<>(field, 0.001, 1000, absTolerance, relTolerance);
225         integrator.setInitialStepSize(60);
226         FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
227         propagator.setOrbitType(OrbitType.EQUINOCTIAL);
228         propagator.setInitialState(initialState);
229 
230         FieldEclipseDetector<T> e = new FieldEclipseDetector<>(field, sun, sunRadius, earth).
231                                     withMaxCheck(120.0).
232                                     withThreshold(zero.newInstance(1e-4)).
233                                     withHandler(new FieldStopOnDecreasing<>()).
234                                     withMaxIter(12).
235                                     withMargin(zero.newInstance(0.001));
236         Assertions.assertEquals(120.0, e.getMaxCheckInterval().currentInterval(null, true), 1.0e-15);
237         Assertions.assertEquals(1.0e-4, e.getThreshold().getReal(), 1.0e-15);
238         Assertions.assertEquals(12, e.getMaxIterationCount());
239         propagator.addEventDetector(e);
240         final FieldSpacecraftState<T> finalState = propagator.propagate(iniDate.shiftedBy(6000));
241         Assertions.assertEquals(2304.188978, finalState.getDate().durationFrom(iniDate).getReal(), 1.0e-4);
242 
243     }
244 
245     private <T extends CalculusFieldElement<T>> void doTestInsideOcculting(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 
266         FieldEclipseDetector<T> e = new FieldEclipseDetector<>(field, sun, sunRadius, earth).
267                                     withMaxCheck(60.0).
268                                     withThreshold(zero.newInstance(1e-3));
269         FieldSpacecraftState<T> s = new FieldSpacecraftState<>(new FieldCartesianOrbit<>(new TimeStampedFieldPVCoordinates<>(FieldAbsoluteDate.getJ2000Epoch(field),
270                                                                                                                              new FieldPVCoordinates<>(new FieldVector3D<>(zero.newInstance(1e6),
271                                                                                                                                                                           zero.newInstance(2e6),
272                                                                                                                                                                           zero.newInstance(3e6)),
273                                                                                                                                                       new FieldVector3D<>(zero.newInstance(1000),
274                                                                                                                                                                           zero.newInstance(0),
275                                                                                                                                                                           zero.newInstance(0)))),
276                                                                                          FramesFactory.getGCRF(),
277                                                                                          zero.add(mu)));
278         try {
279             e.g(s);
280             Assertions.fail("an exception should have been thrown");
281         } catch (OrekitException oe) {
282             Assertions.assertEquals(OrekitMessages.POINT_INSIDE_ELLIPSOID, oe.getSpecifier());
283         }
284     }
285 
286     private <T extends CalculusFieldElement<T>> void doTestInsideOcculted(Field<T> field) {
287         T zero = field.getZero();
288         T one  = field.getOne();
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 
308         FieldEclipseDetector<T> e = new FieldEclipseDetector<>(field, sun, sunRadius, earth).
309                                     withMaxCheck(60.0).
310                                     withThreshold(zero.newInstance(1e-3));
311         Vector3D p = sun.getPVCoordinates(AbsoluteDate.J2000_EPOCH,
312                                           FramesFactory.getGCRF()).getPosition();
313         FieldSpacecraftState<T> s = new FieldSpacecraftState<>(new FieldCartesianOrbit<>(new TimeStampedFieldPVCoordinates<>(FieldAbsoluteDate.getJ2000Epoch(field),
314                                                                                                                              new FieldPVCoordinates<>(new FieldVector3D<>(one,
315                                                                                                                                                                           zero,
316                                                                                                                                                                           zero).add(p),
317                                                                                                                                                       new FieldVector3D<>(zero,
318                                                                                                                                                                           zero,
319                                                                                                                                                                           one))),
320                                                                                          FramesFactory.getGCRF(),
321                                                                                          zero.add(mu)));
322         Assertions.assertEquals(FastMath.PI, e.g(s).getReal(), 1.0e-15);
323     }
324 
325     private <T extends CalculusFieldElement<T>> void testTooSmallMaxIterationCount(Field<T> field) {
326         T zero = field.getZero();
327         final FieldVector3D<T> position  = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.560), zero.add(-25767.25680));
328         final FieldVector3D<T> velocity  = new FieldVector3D<>(zero.add(505.8479685), zero.add(942.7809215), zero.add(7435.922231));
329         FieldAbsoluteDate<T> iniDate = new FieldAbsoluteDate<>(field, 1969, 7, 28, 4, 0, 0.0, TimeScalesFactory.getTT());
330         final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position,  velocity),
331                                                                 FramesFactory.getGCRF(), iniDate, zero.add(mu));
332         FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(orbit);
333         double[] absTolerance = {
334             0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001
335         };
336         double[] relTolerance = {
337             1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7
338         };
339         AdaptiveStepsizeFieldIntegrator<T> integrator =
340             new DormandPrince853FieldIntegrator<>(field, 0.001, 1000, absTolerance, relTolerance);
341         integrator.setInitialStepSize(60);
342         FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
343         propagator.setOrbitType(OrbitType.EQUINOCTIAL);
344         propagator.setInitialState(initialState);
345 
346         int n = 5;
347         FieldEclipseDetector<T> e = new FieldEclipseDetector<>(field, sun, sunRadius, earth).
348                                     withMaxCheck(120.0).
349                                     withThreshold(zero.newInstance(1e-4)).
350                                     withHandler(new FieldStopOnDecreasing<T>()).
351                                     withMaxIter(n);
352        propagator.addEventDetector(e);
353         try {
354             propagator.propagate(iniDate.shiftedBy(6000));
355             Assertions.fail("an exception should have been thrown");
356         } catch (OrekitException oe) {
357             Assertions.assertEquals(LocalizedODEFormats.FIND_ROOT,
358                                     ((MathRuntimeException) oe.getCause()).getSpecifier());
359         }
360     }
361 
362 }
363 //