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