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.geometry.euclidean.threed.Vector3D;
20  import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
21  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
22  import org.hipparchus.util.FastMath;
23  import org.hipparchus.util.MathUtils;
24  import org.junit.jupiter.api.Assertions;
25  import org.junit.jupiter.api.BeforeEach;
26  import org.junit.jupiter.api.Test;
27  import org.orekit.Utils;
28  import org.orekit.errors.OrekitIllegalArgumentException;
29  import org.orekit.errors.OrekitMessages;
30  import org.orekit.forces.ForceModel;
31  import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
32  import org.orekit.forces.gravity.potential.GravityFieldFactory;
33  import org.orekit.forces.gravity.potential.ICGEMFormatReader;
34  import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
35  import org.orekit.frames.Frame;
36  import org.orekit.frames.FramesFactory;
37  import org.orekit.orbits.CartesianOrbit;
38  import org.orekit.orbits.KeplerianOrbit;
39  import org.orekit.orbits.Orbit;
40  import org.orekit.orbits.OrbitType;
41  import org.orekit.orbits.PositionAngleType;
42  import org.orekit.propagation.*;
43  import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
44  import org.orekit.propagation.analytical.KeplerianPropagator;
45  import org.orekit.propagation.events.EventsLogger.LoggedEvent;
46  import org.orekit.propagation.events.handlers.ContinueOnEvent;
47  import org.orekit.propagation.numerical.NumericalPropagator;
48  import org.orekit.time.AbsoluteDate;
49  import org.orekit.time.TimeScale;
50  import org.orekit.time.TimeScalesFactory;
51  import org.orekit.utils.Constants;
52  import org.orekit.utils.IERSConventions;
53  import org.orekit.utils.PVCoordinates;
54  
55  public class PositionAngleDetectorTest {
56  
57      @Test
58      public void testCartesian() {
59          try {
60              new PositionAngleDetector(OrbitType.CARTESIAN, PositionAngleType.TRUE, 0.0).
61              withMaxCheck(600.0).
62              withThreshold(1.0e-6);
63              Assertions.fail("an exception should habe been thrown");
64          } catch (OrekitIllegalArgumentException oiae) {
65              Assertions.assertEquals(OrekitMessages.ORBIT_TYPE_NOT_ALLOWED, oiae.getSpecifier());
66              Assertions.assertEquals(OrbitType.CARTESIAN, oiae.getParts()[0]);
67          }
68      }
69  
70      @Test
71      public void testTrueAnomalyForward() {
72          doTest(OrbitType.KEPLERIAN, PositionAngleType.TRUE, FastMath.toRadians(10.0), Constants.JULIAN_DAY, 15);
73      }
74  
75      @Test
76      public void testTrueAnomalyBackward() {
77          doTest(OrbitType.KEPLERIAN, PositionAngleType.TRUE, FastMath.toRadians(10.0), -Constants.JULIAN_DAY, 14);
78      }
79  
80      @Test
81      public void testMeanAnomalyForward() {
82          doTest(OrbitType.KEPLERIAN, PositionAngleType.MEAN, FastMath.toRadians(10.0), Constants.JULIAN_DAY, 15);
83      }
84  
85      @Test
86      public void testMeanAnomalyBackward() {
87          doTest(OrbitType.KEPLERIAN, PositionAngleType.MEAN, FastMath.toRadians(10.0), -Constants.JULIAN_DAY, 14);
88      }
89  
90      @Test
91      public void testEccentricAnomalyForward() {
92          doTest(OrbitType.KEPLERIAN, PositionAngleType.ECCENTRIC, FastMath.toRadians(10.0), Constants.JULIAN_DAY, 15);
93      }
94  
95      @Test
96      public void testEccentricAnomalyBackward() {
97          doTest(OrbitType.KEPLERIAN, PositionAngleType.ECCENTRIC, FastMath.toRadians(10.0), -Constants.JULIAN_DAY, 14);
98      }
99  
100     @Test
101     public void testTrueLatitudeArgumentForward() {
102         doTest(OrbitType.CIRCULAR, PositionAngleType.TRUE, FastMath.toRadians(730.0), Constants.JULIAN_DAY, 15);
103     }
104 
105     @Test
106     public void testTrueLatitudeArgumentBackward() {
107         doTest(OrbitType.CIRCULAR, PositionAngleType.TRUE, FastMath.toRadians(730.0), -Constants.JULIAN_DAY, 14);
108     }
109 
110     @Test
111     public void testMeanLatitudeArgumentForward() {
112         doTest(OrbitType.CIRCULAR, PositionAngleType.MEAN, FastMath.toRadians(730.0), Constants.JULIAN_DAY, 15);
113     }
114 
115     @Test
116     public void testMeanLatitudeArgumentBackward() {
117         doTest(OrbitType.CIRCULAR, PositionAngleType.MEAN, FastMath.toRadians(730.0), -Constants.JULIAN_DAY, 14);
118     }
119 
120     @Test
121     public void testEccentricLatitudeArgumentForward() {
122         doTest(OrbitType.CIRCULAR, PositionAngleType.ECCENTRIC, FastMath.toRadians(730.0), Constants.JULIAN_DAY, 15);
123     }
124 
125     @Test
126     public void testEccentricLatitudeArgumentBackward() {
127         doTest(OrbitType.CIRCULAR, PositionAngleType.ECCENTRIC, FastMath.toRadians(730.0), -Constants.JULIAN_DAY, 14);
128     }
129 
130     @Test
131     public void testTrueLongitudeArgumentForward() {
132         doTest(OrbitType.EQUINOCTIAL, PositionAngleType.TRUE, FastMath.toRadians(-45.0), Constants.JULIAN_DAY, 15);
133     }
134 
135     @Test
136     public void testTrueLongitudeArgumentBackward() {
137         doTest(OrbitType.EQUINOCTIAL, PositionAngleType.TRUE, FastMath.toRadians(-45.0), -Constants.JULIAN_DAY, 14);
138     }
139 
140     @Test
141     public void testMeanLongitudeArgumentForward() {
142         doTest(OrbitType.EQUINOCTIAL, PositionAngleType.MEAN, FastMath.toRadians(-45.0), Constants.JULIAN_DAY, 15);
143     }
144 
145     @Test
146     public void testMeanLongitudeArgumentBackward() {
147         doTest(OrbitType.EQUINOCTIAL, PositionAngleType.MEAN, FastMath.toRadians(-45.0), -Constants.JULIAN_DAY, 14);
148     }
149 
150     @Test
151     public void testEccentricLongitudeArgumentForward() {
152         doTest(OrbitType.EQUINOCTIAL, PositionAngleType.ECCENTRIC, FastMath.toRadians(-45.0), Constants.JULIAN_DAY, 15);
153     }
154 
155     @Test
156     public void testEccentricLongitudeArgumentBackward() {
157         doTest(OrbitType.EQUINOCTIAL, PositionAngleType.ECCENTRIC, FastMath.toRadians(-45.0), -Constants.JULIAN_DAY, 14);
158     }
159 
160     @Test
161     public void testIssue493() {
162 
163         GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", false));
164         NormalizedSphericalHarmonicsProvider provider =
165                         GravityFieldFactory.getNormalizedProvider(10, 10);
166 
167         Frame inertialFrame = FramesFactory.getEME2000();
168 
169         TimeScale utc = TimeScalesFactory.getUTC();
170         AbsoluteDate initialDate = new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, utc);
171 
172         double mu =  provider.getMu();
173 
174         double a = 24396159;                 // semi major axis in meters
175         double e = 0.72831215;               // eccentricity
176         double i = FastMath.toRadians(7);        // inclination
177         double omega = FastMath.toRadians(180);  // perigee argument
178         double raan = FastMath.toRadians(261);   // right ascension of ascending node
179         double lM = 0;                       // mean anomaly
180 
181         Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngleType.MEAN,
182                                                 inertialFrame, initialDate, mu);
183 
184         // Initial state definition
185         SpacecraftState initialState = new SpacecraftState(initialOrbit);
186 
187         // Adaptive step integrator
188         // with a minimum step of 0.001 and a maximum step of 1000
189         double minStep = 0.001;
190         double maxstep = 1000.0;
191         double positionTolerance = 10.0;
192         OrbitType propagationType = OrbitType.KEPLERIAN;
193         double[][] tolerances =
194                 ToleranceProvider.getDefaultToleranceProvider(positionTolerance).getTolerances(initialOrbit, propagationType);
195         AdaptiveStepsizeIntegrator integrator =
196                         new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
197 
198         // Propagator in Keplerian mode
199         NumericalPropagator propagator = new NumericalPropagator(integrator);
200         propagator.setOrbitType(propagationType);
201 
202         // Simple gravity field force model
203         ForceModel holmesFeatherstone =
204                         new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010,true),
205                                                               provider);
206 
207         propagator.addForceModel(holmesFeatherstone);
208 
209         final double maxCheck  = 600.0;
210         final double threshold = 1.0e-6;
211         PositionAngleDetector detector01 = new PositionAngleDetector(maxCheck,
212                                                                      threshold,
213                                                                      propagationType,
214                                                                      PositionAngleType.TRUE,
215                                                                      FastMath.toRadians(01.0)).
216                                            withHandler(new ContinueOnEvent());
217         PositionAngleDetector detector90 = new PositionAngleDetector(maxCheck,
218                                                                      threshold,
219                                                                      propagationType,
220                                                                      PositionAngleType.TRUE,
221                                                                      FastMath.toRadians(90.0)).
222                                            withHandler(new ContinueOnEvent());
223 
224         // detect events with numerical propagator (and generate ephemeris)
225         final EphemerisGenerator generator = propagator.getEphemerisGenerator();
226         propagator.setInitialState(initialState);
227         EventsLogger logger1 = new EventsLogger();
228         propagator.addEventDetector(logger1.monitorDetector(detector01));
229         propagator.addEventDetector(logger1.monitorDetector(detector90));
230         final AbsoluteDate finalDate = propagator.propagate(new AbsoluteDate(initialDate, Constants.JULIAN_DAY)).getDate();
231         final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
232         Assertions.assertEquals(6, logger1.getLoggedEvents().size());
233 
234         // detect events with generated ephemeris
235         EventsLogger logger2 = new EventsLogger();
236         ephemeris.addEventDetector(logger2.monitorDetector(detector01));
237         ephemeris.addEventDetector(logger2.monitorDetector(detector90));
238         ephemeris.propagate(initialDate, finalDate);
239         Assertions.assertEquals(logger1.getLoggedEvents().size(), logger2.getLoggedEvents().size());
240         for (int k = 0; k < logger1.getLoggedEvents().size(); ++k) {
241             AbsoluteDate date1 = logger1.getLoggedEvents().get(k).getState().getDate();
242             AbsoluteDate date2 = logger2.getLoggedEvents().get(k).getState().getDate();
243             Assertions.assertEquals(0.0, date2.durationFrom(date1), threshold);
244         }
245 
246     }
247 
248     
249     @Test
250     public void testIssue779() {
251 
252         final TimeScale utc = TimeScalesFactory.getUTC();
253         final AbsoluteDate initialDate = new AbsoluteDate(2022, 05, 10, 00, 00, 00.000, utc);
254 
255         // General orbit
256         double a = 24396159;                     // semi major axis in meters
257         double e = 0.72831215;                   // eccentricity
258         double i = FastMath.toRadians(7);        // inclination
259         double omega = FastMath.toRadians(180);  // perigee argument
260         double raan = FastMath.toRadians(261);   // right ascension of ascending node
261         double lM = 0;                           // mean anomaly
262 
263         Orbit orbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngleType.MEAN,
264         									   FramesFactory.getEME2000(), initialDate, 
265         									   Constants.WGS84_EARTH_MU);
266                 
267        
268         
269 
270         // Event detectors configurations
271         final double maxCheck  = 600.0;
272         final double threshold = 1.0e-6;
273         final double angle = 90.0;
274         
275         
276         // First, configure a propagation with the default event handler (expected to stop on event)
277         Propagator np1 = new KeplerianPropagator(orbit);
278         PositionAngleDetector detectorDefault = new PositionAngleDetector(maxCheck,
279                                                                      threshold,
280                                                                      OrbitType.KEPLERIAN,
281                                                                      PositionAngleType.MEAN,
282                                                                      FastMath.toRadians(angle));
283         
284         // Create and add event logger
285         EventsLogger logger1 = new EventsLogger();
286         np1.addEventDetector(logger1.monitorDetector(detectorDefault));
287         
288         // Propagate
289         np1.propagate(initialDate.shiftedBy(Constants.JULIAN_DAY));
290         
291 
292         
293         // Then, repeat the test with a Continue On Event handler
294         Propagator np2 = new KeplerianPropagator(orbit);
295         
296         // Create and configure the Continue On Event handler
297         PositionAngleDetector detectorContinue = new PositionAngleDetector(maxCheck,
298                                                                      threshold,
299                                                                      OrbitType.KEPLERIAN,
300                                                                      PositionAngleType.MEAN,
301                                                                      FastMath.toRadians(angle)).
302                                            				withHandler(new ContinueOnEvent());
303 
304         // Create and add the event logger
305         EventsLogger logger2 = new EventsLogger();
306         np2.addEventDetector(logger2.monitorDetector(detectorContinue));
307         
308         // Propagate
309         np2.propagate(initialDate.shiftedBy(Constants.JULIAN_DAY));
310         
311         // Check test results
312         Assertions.assertEquals(1, logger1.getLoggedEvents().size());
313         Assertions.assertTrue(logger2.getLoggedEvents().size() > 1);
314     }
315     
316     
317     private void doTest(final OrbitType orbitType, final PositionAngleType positionAngleType,
318                         final double angle, final double deltaT, final int expectedCrossings) {
319 
320         PositionAngleDetector d =
321                 new PositionAngleDetector(orbitType, positionAngleType, angle).
322                 withMaxCheck(60).
323                 withThreshold(1.e-10).
324                 withHandler(new ContinueOnEvent());
325 
326         Assertions.assertEquals(60.0, d.getMaxCheckInterval().currentInterval(null, true), 1.0e-15);
327         Assertions.assertEquals(1.0e-10, d.getThreshold(), 1.0e-15);
328         Assertions.assertEquals(orbitType, d.getOrbitType());
329         Assertions.assertEquals(positionAngleType, d.getPositionAngleType());
330         Assertions.assertEquals(angle, d.getAngle(), 1.0e-14);
331         Assertions.assertEquals(AbstractDetector.DEFAULT_MAX_ITER, d.getMaxIterationCount());
332 
333         final TimeScale utc = TimeScalesFactory.getUTC();
334         final Vector3D position = new Vector3D(-6142438.668, 3492467.56, -25767.257);
335         final Vector3D velocity = new Vector3D(506.0, 943.0, 7450);
336         final AbsoluteDate date = new AbsoluteDate(2003, 9, 16, utc);
337         final Orbit orbit = new CartesianOrbit(new PVCoordinates(position,  velocity),
338                                                FramesFactory.getEME2000(), date,
339                                                Constants.EIGEN5C_EARTH_MU);
340 
341         Propagator propagator =
342             new EcksteinHechlerPropagator(orbit,
343                                           Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
344                                           Constants.EIGEN5C_EARTH_MU,
345                                           Constants.EIGEN5C_EARTH_C20,
346                                           Constants.EIGEN5C_EARTH_C30,
347                                           Constants.EIGEN5C_EARTH_C40,
348                                           Constants.EIGEN5C_EARTH_C50,
349                                           Constants.EIGEN5C_EARTH_C60);
350 
351         EventsLogger logger = new EventsLogger();
352         propagator.addEventDetector(logger.monitorDetector(d));
353 
354         propagator.propagate(date.shiftedBy(deltaT));
355 
356         double[] array = new double[6];
357         for (LoggedEvent e : logger.getLoggedEvents()) {
358             SpacecraftState state = e.getState();
359             orbitType.mapOrbitToArray(state.getOrbit(), positionAngleType, array, null);
360             Assertions.assertEquals(angle, MathUtils.normalizeAngle(array[5], angle), 1.0e-10);
361             Assertions.assertEquals(state.getDate(), e.getDate());
362         }
363         Assertions.assertEquals(expectedCrossings, logger.getLoggedEvents().size());
364 
365     }
366 
367     @BeforeEach
368     public void setUp() {
369         Utils.setDataRoot("regular-data:potential");
370     }
371 
372 }
373