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.analysis.differentiation.FieldUnivariateDerivative1;
20  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
21  import org.hipparchus.util.Binary64;
22  import org.hipparchus.util.Binary64Field;
23  import org.hipparchus.util.FastMath;
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.TestUtils;
28  import org.orekit.Utils;
29  import org.orekit.bodies.GeodeticPoint;
30  import org.orekit.bodies.OneAxisEllipsoid;
31  import org.orekit.frames.FieldKinematicTransform;
32  import org.orekit.frames.FramesFactory;
33  import org.orekit.frames.TopocentricFrame;
34  import org.orekit.orbits.FieldEquinoctialOrbit;
35  import org.orekit.orbits.FieldOrbit;
36  import org.orekit.propagation.FieldPropagator;
37  import org.orekit.propagation.FieldSpacecraftState;
38  import org.orekit.propagation.SpacecraftState;
39  import org.orekit.propagation.analytical.FieldEcksteinHechlerPropagator;
40  import org.orekit.propagation.events.FieldEventsLogger.FieldLoggedEvent;
41  import org.orekit.propagation.events.handlers.FieldContinueOnEvent;
42  import org.orekit.time.AbsoluteDate;
43  import org.orekit.time.FieldAbsoluteDate;
44  import org.orekit.time.TimeScale;
45  import org.orekit.time.TimeScalesFactory;
46  import org.orekit.utils.Constants;
47  import org.orekit.utils.FieldPVCoordinates;
48  import org.orekit.utils.IERSConventions;
49  import org.orekit.utils.TimeStampedFieldPVCoordinates;
50  
51  class FieldElevationExtremumDetectorTest {
52  
53  
54      @Test
55      void testG() {
56          // GIVEN
57          final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
58                  Constants.WGS84_EARTH_FLATTENING,
59                  FramesFactory.getITRF(IERSConventions.IERS_2010, true));
60          final GeodeticPoint gp = new GeodeticPoint(FastMath.toRadians(51.0), FastMath.toRadians(66.6), 300.0);
61          final FieldElevationExtremumDetector<Binary64> detector =
62                  new FieldElevationExtremumDetector<>(new Binary64(0.001), new Binary64(1.e-6), new TopocentricFrame(earth, gp, "test")).
63                          withHandler(new FieldContinueOnEvent<>());
64          final SpacecraftState state = new SpacecraftState(TestUtils.getDefaultOrbit(AbsoluteDate.ARBITRARY_EPOCH));
65          final FieldSpacecraftState<Binary64> fieldState = new FieldSpacecraftState<>(Binary64Field.getInstance(), state);
66          // WHEN
67          final Binary64 actualG = detector.g(fieldState);
68          // THEN
69          final FieldKinematicTransform<Binary64> inertToTopo = fieldState.getFrame().getKinematicTransformTo(detector.getTopocentricFrame(),
70                  fieldState.getDate());
71          final TimeStampedFieldPVCoordinates<Binary64> pvTopo = inertToTopo.transformOnlyPV(fieldState.getPVCoordinates());
72          final FieldVector3D<FieldUnivariateDerivative1<Binary64>> pvDS = pvTopo.toUnivariateDerivative1Vector();
73          final FieldUnivariateDerivative1<Binary64> elevation = pvDS.getZ().divide(pvDS.getNorm()).asin();
74          final Binary64 expectedG = elevation.getFirstDerivative();
75          Assertions.assertEquals(expectedG.getReal(), actualG.getReal(), 1e-12);
76      }
77  
78      @Test
79      void testLEO() {
80  
81          final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
82                                                              Constants.WGS84_EARTH_FLATTENING,
83                                                              FramesFactory.getITRF(IERSConventions.IERS_2010, true));
84  
85          final GeodeticPoint gp = new GeodeticPoint(FastMath.toRadians(51.0), FastMath.toRadians(66.6), 300.0);
86          final FieldElevationExtremumDetector<Binary64> raw =
87                  new FieldElevationExtremumDetector<>(Binary64Field.getInstance(), new TopocentricFrame(earth, gp, "test")).
88                  withMaxCheck(60.0).
89                  withThreshold(new Binary64(1.e-6)).
90                  withHandler(new FieldContinueOnEvent<>());
91          final FieldEventSlopeFilter<FieldElevationExtremumDetector<Binary64>, Binary64> maxElevationDetector =
92                  new FieldEventSlopeFilter<>(raw, FilterType.TRIGGER_ONLY_DECREASING_EVENTS);
93  
94          Assertions.assertEquals(60.0, raw.getMaxCheckInterval().currentInterval(null, true), 1.0e-15);
95          Assertions.assertEquals(1.0e-6, raw.getThreshold().getReal(), 1.0e-15);
96          Assertions.assertEquals(AbstractDetector.DEFAULT_MAX_ITER, raw.getMaxIterationCount());
97          Assertions.assertEquals("test", raw.getTopocentricFrame().getName());
98  
99          final TimeScale utc = TimeScalesFactory.getUTC();
100         final FieldVector3D<Binary64> position = new FieldVector3D<>(new Binary64(-6142438.668),
101                                                                      new Binary64(3492467.56),
102                                                                      new Binary64(-25767.257));
103         final FieldVector3D<Binary64> velocity = new FieldVector3D<>(new Binary64(505.848),
104                                                                      new Binary64(942.781),
105                                                                      new Binary64(7435.922));
106         final FieldAbsoluteDate<Binary64> date = new FieldAbsoluteDate<>(Binary64Field.getInstance(),
107                                                                          new AbsoluteDate(2003, 9, 16, utc));
108         final FieldOrbit<Binary64> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position,  velocity),
109                                                                        FramesFactory.getEME2000(), date,
110                                                                        new Binary64(Constants.EIGEN5C_EARTH_MU));
111 
112         FieldPropagator<Binary64> propagator =
113             new FieldEcksteinHechlerPropagator<>(orbit,
114                                                  Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
115                                                  new Binary64(Constants.EIGEN5C_EARTH_MU),
116                                                  Constants.EIGEN5C_EARTH_C20,
117                                                  Constants.EIGEN5C_EARTH_C30,
118                                                  Constants.EIGEN5C_EARTH_C40,
119                                                  Constants.EIGEN5C_EARTH_C50,
120                                                  Constants.EIGEN5C_EARTH_C60);
121 
122         FieldEventsLogger<Binary64> logger = new FieldEventsLogger<>();
123         propagator.addEventDetector(logger.monitorDetector(maxElevationDetector));
124 
125         propagator.propagate(date.shiftedBy(Constants.JULIAN_DAY));
126         int visibleEvents = 0;
127         for (FieldLoggedEvent<Binary64> e : logger.getLoggedEvents()) {
128             final double eMinus = raw.getElevation(e.getState().shiftedBy(-10.0)).getReal();
129             final double e0     = raw.getElevation(e.getState()).getReal();
130             final double ePlus  = raw.getElevation(e.getState().shiftedBy(+10.0)).getReal();
131             if (e0 > FastMath.toRadians(5.0)) {
132                 ++visibleEvents;
133             }
134             Assertions.assertTrue(e0 > eMinus);
135             Assertions.assertTrue(e0 > ePlus);
136         }
137         Assertions.assertEquals(15, logger.getLoggedEvents().size());
138         Assertions.assertEquals( 6, visibleEvents);
139 
140     }
141 
142     @BeforeEach
143     void setUp() {
144         Utils.setDataRoot("regular-data");
145     }
146 
147 }
148