1   /* Copyright 2002-2022 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.Assert;
25  import org.junit.Before;
26  import org.junit.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.PositionAngle;
42  import org.orekit.propagation.BoundedPropagator;
43  import org.orekit.propagation.EphemerisGenerator;
44  import org.orekit.propagation.Propagator;
45  import org.orekit.propagation.SpacecraftState;
46  import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
47  import org.orekit.propagation.events.EventsLogger.LoggedEvent;
48  import org.orekit.propagation.events.handlers.ContinueOnEvent;
49  import org.orekit.propagation.numerical.NumericalPropagator;
50  import org.orekit.time.AbsoluteDate;
51  import org.orekit.time.TimeScale;
52  import org.orekit.time.TimeScalesFactory;
53  import org.orekit.utils.Constants;
54  import org.orekit.utils.IERSConventions;
55  import org.orekit.utils.PVCoordinates;
56  
57  public class PositionAngleDetectorTest {
58  
59      @Test
60      public void testCartesian() {
61          try {
62              new PositionAngleDetector(OrbitType.CARTESIAN, PositionAngle.TRUE, 0.0).
63              withMaxCheck(600.0).
64              withThreshold(1.0e-6);
65              Assert.fail("an exception should habe been thrown");
66          } catch (OrekitIllegalArgumentException oiae) {
67              Assert.assertEquals(OrekitMessages.ORBIT_TYPE_NOT_ALLOWED, oiae.getSpecifier());
68              Assert.assertEquals(OrbitType.CARTESIAN, oiae.getParts()[0]);
69          }
70      }
71  
72      @Test
73      public void testTrueAnomalyForward() {
74          doTest(OrbitType.KEPLERIAN, PositionAngle.TRUE, FastMath.toRadians(10.0), Constants.JULIAN_DAY, 15);
75      }
76  
77      @Test
78      public void testTrueAnomalyBackward() {
79          doTest(OrbitType.KEPLERIAN, PositionAngle.TRUE, FastMath.toRadians(10.0), -Constants.JULIAN_DAY, 14);
80      }
81  
82      @Test
83      public void testMeanAnomalyForward() {
84          doTest(OrbitType.KEPLERIAN, PositionAngle.MEAN, FastMath.toRadians(10.0), Constants.JULIAN_DAY, 15);
85      }
86  
87      @Test
88      public void testMeanAnomalyBackward() {
89          doTest(OrbitType.KEPLERIAN, PositionAngle.MEAN, FastMath.toRadians(10.0), -Constants.JULIAN_DAY, 14);
90      }
91  
92      @Test
93      public void testEccentricAnomalyForward() {
94          doTest(OrbitType.KEPLERIAN, PositionAngle.ECCENTRIC, FastMath.toRadians(10.0), Constants.JULIAN_DAY, 15);
95      }
96  
97      @Test
98      public void testEccentricAnomalyBackward() {
99          doTest(OrbitType.KEPLERIAN, PositionAngle.ECCENTRIC, FastMath.toRadians(10.0), -Constants.JULIAN_DAY, 14);
100     }
101 
102     @Test
103     public void testTrueLatitudeArgumentForward() {
104         doTest(OrbitType.CIRCULAR, PositionAngle.TRUE, FastMath.toRadians(730.0), Constants.JULIAN_DAY, 15);
105     }
106 
107     @Test
108     public void testTrueLatitudeArgumentBackward() {
109         doTest(OrbitType.CIRCULAR, PositionAngle.TRUE, FastMath.toRadians(730.0), -Constants.JULIAN_DAY, 14);
110     }
111 
112     @Test
113     public void testMeanLatitudeArgumentForward() {
114         doTest(OrbitType.CIRCULAR, PositionAngle.MEAN, FastMath.toRadians(730.0), Constants.JULIAN_DAY, 15);
115     }
116 
117     @Test
118     public void testMeanLatitudeArgumentBackward() {
119         doTest(OrbitType.CIRCULAR, PositionAngle.MEAN, FastMath.toRadians(730.0), -Constants.JULIAN_DAY, 14);
120     }
121 
122     @Test
123     public void testEccentricLatitudeArgumentForward() {
124         doTest(OrbitType.CIRCULAR, PositionAngle.ECCENTRIC, FastMath.toRadians(730.0), Constants.JULIAN_DAY, 15);
125     }
126 
127     @Test
128     public void testEccentricLatitudeArgumentBackward() {
129         doTest(OrbitType.CIRCULAR, PositionAngle.ECCENTRIC, FastMath.toRadians(730.0), -Constants.JULIAN_DAY, 14);
130     }
131 
132     @Test
133     public void testTrueLongitudeArgumentForward() {
134         doTest(OrbitType.EQUINOCTIAL, PositionAngle.TRUE, FastMath.toRadians(-45.0), Constants.JULIAN_DAY, 15);
135     }
136 
137     @Test
138     public void testTrueLongitudeArgumentBackward() {
139         doTest(OrbitType.EQUINOCTIAL, PositionAngle.TRUE, FastMath.toRadians(-45.0), -Constants.JULIAN_DAY, 14);
140     }
141 
142     @Test
143     public void testMeanLongitudeArgumentForward() {
144         doTest(OrbitType.EQUINOCTIAL, PositionAngle.MEAN, FastMath.toRadians(-45.0), Constants.JULIAN_DAY, 15);
145     }
146 
147     @Test
148     public void testMeanLongitudeArgumentBackward() {
149         doTest(OrbitType.EQUINOCTIAL, PositionAngle.MEAN, FastMath.toRadians(-45.0), -Constants.JULIAN_DAY, 14);
150     }
151 
152     @Test
153     public void testEccentricLongitudeArgumentForward() {
154         doTest(OrbitType.EQUINOCTIAL, PositionAngle.ECCENTRIC, FastMath.toRadians(-45.0), Constants.JULIAN_DAY, 15);
155     }
156 
157     @Test
158     public void testEccentricLongitudeArgumentBackward() {
159         doTest(OrbitType.EQUINOCTIAL, PositionAngle.ECCENTRIC, FastMath.toRadians(-45.0), -Constants.JULIAN_DAY, 14);
160     }
161 
162     @Test
163     public void testIssue493() {
164 
165         GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", false));
166         NormalizedSphericalHarmonicsProvider provider =
167                         GravityFieldFactory.getNormalizedProvider(10, 10);
168 
169         Frame inertialFrame = FramesFactory.getEME2000();
170 
171         TimeScale utc = TimeScalesFactory.getUTC();
172         AbsoluteDate initialDate = new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, utc);
173 
174         double mu =  provider.getMu();
175 
176         double a = 24396159;                 // semi major axis in meters
177         double e = 0.72831215;               // eccentricity
178         double i = FastMath.toRadians(7);        // inclination
179         double omega = FastMath.toRadians(180);  // perigee argument
180         double raan = FastMath.toRadians(261);   // right ascension of ascending node
181         double lM = 0;                       // mean anomaly
182 
183         Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.MEAN,
184                                                 inertialFrame, initialDate, mu);
185 
186         // Initial state definition
187         SpacecraftState initialState = new SpacecraftState(initialOrbit);
188 
189         // Adaptive step integrator
190         // with a minimum step of 0.001 and a maximum step of 1000
191         double minStep = 0.001;
192         double maxstep = 1000.0;
193         double positionTolerance = 10.0;
194         OrbitType propagationType = OrbitType.KEPLERIAN;
195         double[][] tolerances =
196                         NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType);
197         AdaptiveStepsizeIntegrator integrator =
198                         new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
199 
200         // Propagator in Keplerian mode
201         NumericalPropagator propagator = new NumericalPropagator(integrator);
202         propagator.setOrbitType(propagationType);
203 
204         // Simple gravity field force model
205         ForceModel holmesFeatherstone =
206                         new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010,true),
207                                                               provider);
208 
209         propagator.addForceModel(holmesFeatherstone);
210 
211         final double maxCheck  = 600.0;
212         final double threshold = 1.0e-6;
213         PositionAngleDetector detector01 = new PositionAngleDetector(maxCheck,
214                                                                      threshold,
215                                                                      propagationType,
216                                                                      PositionAngle.TRUE,
217                                                                      FastMath.toRadians(01.0)).
218                                            withHandler(new ContinueOnEvent<>());
219         PositionAngleDetector detector90 = new PositionAngleDetector(maxCheck,
220                                                                      threshold,
221                                                                      propagationType,
222                                                                      PositionAngle.TRUE,
223                                                                      FastMath.toRadians(90.0)).
224                                            withHandler(new ContinueOnEvent<>());
225 
226         // detect events with numerical propagator (and generate ephemeris)
227         final EphemerisGenerator generator = propagator.getEphemerisGenerator();
228         propagator.setInitialState(initialState);
229         EventsLogger logger1 = new EventsLogger();
230         propagator.addEventDetector(logger1.monitorDetector(detector01));
231         propagator.addEventDetector(logger1.monitorDetector(detector90));
232         final AbsoluteDate finalDate = propagator.propagate(new AbsoluteDate(initialDate, Constants.JULIAN_DAY)).getDate();
233         final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
234         Assert.assertEquals(6, logger1.getLoggedEvents().size());
235 
236         // detect events with generated ephemeris
237         EventsLogger logger2 = new EventsLogger();
238         ephemeris.addEventDetector(logger2.monitorDetector(detector01));
239         ephemeris.addEventDetector(logger2.monitorDetector(detector90));
240         ephemeris.propagate(initialDate, finalDate);
241         Assert.assertEquals(logger1.getLoggedEvents().size(), logger2.getLoggedEvents().size());
242         for (int k = 0; k < logger1.getLoggedEvents().size(); ++k) {
243             AbsoluteDate date1 = logger1.getLoggedEvents().get(k).getState().getDate();
244             AbsoluteDate date2 = logger2.getLoggedEvents().get(k).getState().getDate();
245             Assert.assertEquals(0.0, date2.durationFrom(date1), threshold);
246         }
247 
248     }
249 
250     private void doTest(final OrbitType orbitType, final PositionAngle positionAngle,
251                         final double angle, final double deltaT, final int expectedCrossings) {
252 
253         PositionAngleDetector d =
254                 new PositionAngleDetector(orbitType, positionAngle, angle).
255                 withMaxCheck(60).
256                 withThreshold(1.e-10).
257                 withHandler(new ContinueOnEvent<PositionAngleDetector>());
258 
259         Assert.assertEquals(60.0, d.getMaxCheckInterval(), 1.0e-15);
260         Assert.assertEquals(1.0e-10, d.getThreshold(), 1.0e-15);
261         Assert.assertEquals(orbitType, d.getOrbitType());
262         Assert.assertEquals(positionAngle, d.getPositionAngle());
263         Assert.assertEquals(angle, d.getAngle(), 1.0e-14);
264         Assert.assertEquals(AbstractDetector.DEFAULT_MAX_ITER, d.getMaxIterationCount());
265 
266         final TimeScale utc = TimeScalesFactory.getUTC();
267         final Vector3D position = new Vector3D(-6142438.668, 3492467.56, -25767.257);
268         final Vector3D velocity = new Vector3D(506.0, 943.0, 7450);
269         final AbsoluteDate date = new AbsoluteDate(2003, 9, 16, utc);
270         final Orbit orbit = new CartesianOrbit(new PVCoordinates(position,  velocity),
271                                                FramesFactory.getEME2000(), date,
272                                                Constants.EIGEN5C_EARTH_MU);
273 
274         Propagator propagator =
275             new EcksteinHechlerPropagator(orbit,
276                                           Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
277                                           Constants.EIGEN5C_EARTH_MU,
278                                           Constants.EIGEN5C_EARTH_C20,
279                                           Constants.EIGEN5C_EARTH_C30,
280                                           Constants.EIGEN5C_EARTH_C40,
281                                           Constants.EIGEN5C_EARTH_C50,
282                                           Constants.EIGEN5C_EARTH_C60);
283 
284         EventsLogger logger = new EventsLogger();
285         propagator.addEventDetector(logger.monitorDetector(d));
286 
287         propagator.propagate(date.shiftedBy(deltaT));
288 
289         double[] array = new double[6];
290         for (LoggedEvent e : logger.getLoggedEvents()) {
291             SpacecraftState state = e.getState();
292             orbitType.mapOrbitToArray(state.getOrbit(), positionAngle, array, null);
293             Assert.assertEquals(angle, MathUtils.normalizeAngle(array[5], angle), 1.0e-10);
294             Assert.assertEquals(state.getDate(), e.getDate());
295         }
296         Assert.assertEquals(expectedCrossings, logger.getLoggedEvents().size());
297 
298     }
299 
300     @Before
301     public void setUp() {
302         Utils.setDataRoot("regular-data:potential");
303     }
304 
305 }
306