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