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.Line;
20  import org.hipparchus.geometry.euclidean.threed.RotationOrder;
21  import org.hipparchus.geometry.euclidean.threed.Vector3D;
22  import org.hipparchus.geometry.partitioning.RegionFactory;
23  import org.hipparchus.geometry.spherical.twod.Circle;
24  import org.hipparchus.geometry.spherical.twod.S2Point;
25  import org.hipparchus.geometry.spherical.twod.Sphere2D;
26  import org.hipparchus.geometry.spherical.twod.SphericalPolygonsSet;
27  import org.hipparchus.geometry.spherical.twod.SubCircle;
28  import org.hipparchus.util.FastMath;
29  import org.junit.jupiter.api.Assertions;
30  import org.junit.jupiter.api.BeforeEach;
31  import org.junit.jupiter.api.Test;
32  import org.orekit.Utils;
33  import org.orekit.attitudes.LofOffset;
34  import org.orekit.bodies.GeodeticPoint;
35  import org.orekit.bodies.OneAxisEllipsoid;
36  import org.orekit.errors.OrekitException;
37  import org.orekit.frames.Frame;
38  import org.orekit.frames.FramesFactory;
39  import org.orekit.frames.LOFType;
40  import org.orekit.frames.Transform;
41  import org.orekit.geometry.fov.DoubleDihedraFieldOfView;
42  import org.orekit.geometry.fov.FieldOfView;
43  import org.orekit.geometry.fov.PolygonalFieldOfView;
44  import org.orekit.geometry.fov.PolygonalFieldOfView.DefiningConeType;
45  import org.orekit.orbits.EquinoctialOrbit;
46  import org.orekit.orbits.Orbit;
47  import org.orekit.propagation.Propagator;
48  import org.orekit.propagation.SpacecraftState;
49  import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
50  import org.orekit.propagation.events.EventsLogger.LoggedEvent;
51  import org.orekit.propagation.events.handlers.ContinueOnEvent;
52  import org.orekit.time.AbsoluteDate;
53  import org.orekit.time.TimeScale;
54  import org.orekit.time.TimeScalesFactory;
55  import org.orekit.utils.Constants;
56  import org.orekit.utils.IERSConventions;
57  import org.orekit.utils.PVCoordinates;
58  
59  import java.lang.reflect.Field;
60  import java.util.List;
61  
62  public class FootprintOverlapDetectorTest {
63  
64      private Propagator propagator;
65  
66      private Orbit initialOrbit;
67  
68      private OneAxisEllipsoid earth;
69  
70      @Test
71      public void testSampleAroundPole() throws NoSuchFieldException, IllegalAccessException {
72          S2Point[] polygon = new S2Point[] {
73              new S2Point(FastMath.toRadians(-120.0), FastMath.toRadians(5.0)),
74              new S2Point(FastMath.toRadians(   0.0), FastMath.toRadians(5.0)),
75              new S2Point(FastMath.toRadians( 120.0), FastMath.toRadians(5.0))
76          };
77          SphericalPolygonsSet aoi = new SphericalPolygonsSet(1.e-9, polygon);
78          FieldOfView fov = new DoubleDihedraFieldOfView(Vector3D.PLUS_J,
79                                                         Vector3D.PLUS_I, FastMath.toRadians(5.),
80                                                         Vector3D.PLUS_K, FastMath.toRadians(5.),
81                                                         0.);
82          Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
83          OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf);
84          FootprintOverlapDetector detector = new FootprintOverlapDetector(fov, earth, aoi, 20000.);
85          Assertions.assertEquals(9.91475183e-3, detector.getZone().getSize(), 1.0e-10);
86          Field sampledZoneField = FootprintOverlapDetector.class.getDeclaredField("sampledZone");
87          sampledZoneField.setAccessible(true);
88          List<?> sampledZone = (List<?>) sampledZoneField.get(detector);
89          Assertions.assertEquals(1140, sampledZone.size());
90      }
91  
92      @Test
93      public void testRightForwardView() {
94  
95          propagator.setAttitudeProvider(new LofOffset(initialOrbit.getFrame(), LOFType.LVLH_CCSDS,
96                                                        RotationOrder.XYZ,
97                                                        FastMath.toRadians(-20.0),
98                                                        FastMath.toRadians(+20.0),
99                                                        0.0));
100 
101         // observe continental France plus Corsica
102         final SphericalPolygonsSet france = buildFrance();
103 
104         // square field of view along Z axis (which is pointing sideways), aperture 5°, 0° margin
105         final FieldOfView fov = new PolygonalFieldOfView(Vector3D.PLUS_K,
106                                                          DefiningConeType.INSIDE_CONE_TOUCHING_POLYGON_AT_EDGES_MIDDLE,
107                                                          Vector3D.PLUS_I,
108                                                          FastMath.toRadians(2.5), 4, 0.0);
109         final FootprintOverlapDetector detector =
110                 new FootprintOverlapDetector(fov, earth, france, 50000.0).
111                 withMaxCheck(1.0).
112                 withThreshold(1.0e-6).
113                 withHandler(new ContinueOnEvent());
114         final EventsLogger logger = new EventsLogger();
115         propagator.addEventDetector(logger.monitorDetector(detector));
116 
117         // Extrapolate from the initial to the final date
118         propagator.propagate(initialOrbit.getDate().shiftedBy(635000),
119                              initialOrbit.getDate().shiftedBy(735000));
120 
121         List<LoggedEvent> events = logger.getLoggedEvents();
122         Assertions.assertEquals(8, events.size());
123 
124         // the first two consecutive close events occur during the same ascending orbit
125         // we first see Corsica, then lose visibility over the sea, then see continental France
126 
127         // above Mediterranean sea, between Illes Balears and Sardigna,
128         // pointing to Corsica towards North-East
129         checkEventPair(events.get(0),  events.get(1),
130                        639010.0775,  34.1551, 39.2231,  6.5960, 42.0734,  9.0526);
131 
132         // above Saint-Chamond (Loire), pointing near Saint-Dié-des-Vosges (Vosges) towards North-East
133         checkEventPair(events.get(2),  events.get(3),
134                        639113.5532,  38.3899, 45.5356,  4.4813, 48.4211,  7.1499);
135 
136         // event is on a descending orbit, so the pointing direction,
137         // taking roll and pitch offsets, is towards South-West with respect to spacecraft
138         // above English Channel, pointing near Hanvec (Finistère) towards South-West
139         checkEventPair(events.get(4),  events.get(5),
140                        687772.4531,  27.0852, 50.2693,  0.0493, 48.3243, -4.1510);
141 
142         // event on an ascending orbit
143         // above Atlantic ocean, pointing near to île d'Oléron (Charente-Maritime) towards North-East
144         checkEventPair(events.get(6),  events.get(7),
145                        727696.1034, 112.8867, 42.9486, -4.0325, 45.8192, -1.4565);
146 
147     }
148 
149     private void checkEventPair(final LoggedEvent start, final LoggedEvent end,
150                                 final double expectedStart, final double expectedDuration,
151                                 final double spacecraftLatitude, final double spacecraftLongitude,
152                                 final double fovCenterLatitude, final double fovCenterLongitude) {
153 
154         Assertions.assertFalse(start.isIncreasing());
155         Assertions.assertTrue(end.isIncreasing());
156         Assertions.assertEquals(expectedStart,
157                             start.getState().getDate().durationFrom(initialOrbit.getDate()),
158                             0.001);
159         Assertions.assertEquals(expectedDuration,
160                             end.getState().getDate().durationFrom(start.getState().getDate()),
161                             0.001);
162 
163         SpacecraftState middle = start.getState().shiftedBy(0.5 * expectedDuration);
164 
165         // sub-satellite point
166         Vector3D p = middle.getPosition();
167         GeodeticPoint gpSat = earth.transform(p, middle.getFrame(), middle.getDate());
168         Assertions.assertEquals(spacecraftLatitude,  FastMath.toDegrees(gpSat.getLatitude()),  0.001);
169         Assertions.assertEquals(spacecraftLongitude, FastMath.toDegrees(gpSat.getLongitude()), 0.001);
170 
171         // point at center of Field Of View
172         final Transform scToInert = middle.toTransform().getInverse();
173         GeodeticPoint gpFOV =
174                 earth.getIntersectionPoint(new Line(p, scToInert.transformPosition(Vector3D.PLUS_K), 1.0e-6),
175                                            middle.getPosition(),
176                                            middle.getFrame(), middle.getDate());
177         Assertions.assertEquals(fovCenterLatitude,  FastMath.toDegrees(gpFOV.getLatitude()),  0.001);
178         Assertions.assertEquals(fovCenterLongitude, FastMath.toDegrees(gpFOV.getLongitude()), 0.001);
179 
180     }
181 
182     @BeforeEach
183     public void setUp() {
184         try {
185 
186             Utils.setDataRoot("regular-data");
187 
188             final TimeScale utc = TimeScalesFactory.getUTC();
189             final Vector3D position = new Vector3D(-6142438.668, 3492467.56, -25767.257);
190             final Vector3D velocity = new Vector3D(505.848, 942.781, 7435.922);
191             final AbsoluteDate date = new AbsoluteDate(2003, 9, 16, utc);
192             initialOrbit = new EquinoctialOrbit(new PVCoordinates(position,  velocity),
193                                                      FramesFactory.getEME2000(), date,
194                                                      Constants.EIGEN5C_EARTH_MU);
195 
196             propagator =
197                 new EcksteinHechlerPropagator(initialOrbit,
198                                               Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
199                                               Constants.EIGEN5C_EARTH_MU,
200                                               Constants.EIGEN5C_EARTH_C20,
201                                               Constants.EIGEN5C_EARTH_C30,
202                                               Constants.EIGEN5C_EARTH_C40,
203                                               Constants.EIGEN5C_EARTH_C50,
204                                               Constants.EIGEN5C_EARTH_C60);
205 
206             earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
207                                          Constants.WGS84_EARTH_FLATTENING,
208                                          FramesFactory.getITRF(IERSConventions.IERS_2010, true));
209 
210         } catch (OrekitException oe) {
211             Assertions.fail(oe.getMessage());
212         }
213 
214     }
215 
216     private SphericalPolygonsSet buildFrance() {
217 
218         final SphericalPolygonsSet continental = buildSimpleZone(new double[][] {
219             { 51.14850,  2.51357 }, { 50.94660,  1.63900 }, { 50.12717,  1.33876 }, { 49.34737, -0.98946 },
220             { 49.77634, -1.93349 }, { 48.64442, -1.61651 }, { 48.90169, -3.29581 }, { 48.68416, -4.59234 },
221             { 47.95495, -4.49155 }, { 47.57032, -2.96327 }, { 46.01491, -1.19379 }, { 44.02261, -1.38422 },
222             { 43.42280, -1.90135 }, { 43.03401, -1.50277 }, { 42.34338,  1.82679 }, { 42.47301,  2.98599 },
223             { 43.07520,  3.10041 }, { 43.39965,  4.55696 }, { 43.12889,  6.52924 }, { 43.69384,  7.43518 },
224             { 44.12790,  7.54959 }, { 45.02851,  6.74995 }, { 45.33309,  7.09665 }, { 46.42967,  6.50009 },
225             { 46.27298,  6.02260 }, { 46.72577,  6.03738 }, { 47.62058,  7.46675 }, { 49.01778,  8.09927 },
226             { 49.20195,  6.65822 }, { 49.44266,  5.89775 }, { 49.98537,  4.79922 }
227           });
228 
229         final SphericalPolygonsSet corsica = buildSimpleZone(new double[][] {
230             { 42.15249,  9.56001 }, { 43.00998,  9.39000 }, { 42.62812,  8.74600 }, { 42.25651,  8.54421 },
231             { 41.58361,  8.77572 }, { 41.38000,  9.22975 }
232           });
233 
234           return (SphericalPolygonsSet) new RegionFactory<Sphere2D, S2Point, Circle, SubCircle>().
235                  union(continental, corsica);
236 
237     }
238 
239     private SphericalPolygonsSet buildSimpleZone(double[][] points) {
240         final S2Point[] vertices = new S2Point[points.length];
241         for (int i = 0; i < points.length; ++i) {
242             vertices[i] = new S2Point(FastMath.toRadians(points[i][1]),         // points[i][1] is longitude
243                                       FastMath.toRadians(90.0 - points[i][0])); // points[i][0] is latitude
244         }
245         return new SphericalPolygonsSet(1.0e-10, vertices);
246     }
247 
248 }