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 java.util.ArrayList;
20  import java.util.List;
21  
22  import org.hipparchus.geometry.enclosing.EnclosingBall;
23  import org.hipparchus.geometry.euclidean.threed.Vector3D;
24  import org.hipparchus.geometry.spherical.twod.Edge;
25  import org.hipparchus.geometry.spherical.twod.S2Point;
26  import org.hipparchus.geometry.spherical.twod.Sphere2D;
27  import org.hipparchus.geometry.spherical.twod.SphericalPolygonsSet;
28  import org.hipparchus.geometry.spherical.twod.Vertex;
29  import org.hipparchus.ode.events.Action;
30  import org.hipparchus.util.FastMath;
31  import org.hipparchus.util.SinCos;
32  import org.orekit.bodies.BodyShape;
33  import org.orekit.bodies.GeodeticPoint;
34  import org.orekit.bodies.OneAxisEllipsoid;
35  import org.orekit.frames.StaticTransform;
36  import org.orekit.geometry.fov.FieldOfView;
37  import org.orekit.models.earth.tessellation.DivertedSingularityAiming;
38  import org.orekit.models.earth.tessellation.EllipsoidTessellator;
39  import org.orekit.propagation.SpacecraftState;
40  import org.orekit.propagation.events.handlers.EventHandler;
41  import org.orekit.propagation.events.handlers.StopOnIncreasing;
42  
43  /** Detector triggered by geographical region entering/leaving a spacecraft sensor
44   * {@link FieldOfView Field Of View}.
45   * <p>
46   * This detector is a mix between to {@link FieldOfViewDetector} and {@link
47   * GeographicZoneDetector}. Similar to the first detector above, it triggers events
48   * related to entry/exit of targets in a Field Of View, taking attitude into account.
49   * Similar to the second detector above, its target is an entire geographic region
50   * (which can even be split in several non-connected patches and can have holes).
51   * </p>
52   * <p>
53   * This detector is typically used for ground observation missions with agile
54   * satellites than can look away from nadir.
55   * </p>
56   * <p>The default implementation behavior is to {@link Action#CONTINUE continue}
57   * propagation at FOV entry and to {@link Action#STOP stop} propagation
58   * at FOV exit. This can be changed by calling
59   * {@link #withHandler(EventHandler)} after construction.</p>
60   * @see org.orekit.propagation.Propagator#addEventDetector(EventDetector)
61   * @see FieldOfViewDetector
62   * @see GeographicZoneDetector
63   * @author Luc Maisonobe
64   * @since 7.1
65   */
66  public class FootprintOverlapDetector extends AbstractDetector<FootprintOverlapDetector> {
67  
68      /** Field of view. */
69      private final FieldOfView fov;
70  
71      /** Body on which the geographic zone is defined. */
72      private final OneAxisEllipsoid body;
73  
74      /** Geographic zone to consider. */
75      private final SphericalPolygonsSet zone;
76  
77      /** Linear step used for sampling the geographic zone. */
78      private final double samplingStep;
79  
80      /** Sampling of the geographic zone. */
81      private final List<SamplingPoint> sampledZone;
82  
83      /** Center of the spherical cap surrounding the zone. */
84      private final Vector3D capCenter;
85  
86      /** Cosine of the radius of the spherical cap surrounding the zone. */
87      private final double capCos;
88  
89      /** Sine of the radius of the spherical cap surrounding the zone. */
90      private final double capSin;
91  
92      /** Build a new instance.
93       * <p>The maximal interval between distance to FOV boundary checks should
94       * be smaller than the half duration of the minimal pass to handle,
95       * otherwise some short passes could be missed.</p>
96       * @param fov sensor field of view
97       * @param body body on which the geographic zone is defined
98       * @param zone geographic zone to consider
99       * @param samplingStep linear step used for sampling the geographic zone (in meters)
100      * @since 10.1
101      */
102     public FootprintOverlapDetector(final FieldOfView fov,
103                                     final OneAxisEllipsoid body,
104                                     final SphericalPolygonsSet zone,
105                                     final double samplingStep) {
106         this(EventDetectionSettings.getDefaultEventDetectionSettings(), new StopOnIncreasing(),
107              fov, body, zone, samplingStep, sample(body, zone, samplingStep));
108     }
109 
110     /** Protected constructor with full parameters.
111      * <p>
112      * This constructor is not public as users are expected to use the builder
113      * API with the various {@code withXxx()} methods to set up the instance
114      * in a readable manner without using a huge amount of parameters.
115      * </p>
116      * @param detectionSettings event detection settings
117      * @param handler event handler to call at event occurrences
118      * @param body body on which the geographic zone is defined
119      * @param zone geographic zone to consider
120      * @param fov sensor field of view
121      * @param sampledZone sampling of the geographic zone
122      * @param samplingStep linear step used for sampling the geographic zone (in meters)
123      * @since 13.0
124      */
125     protected FootprintOverlapDetector(final EventDetectionSettings detectionSettings, final EventHandler handler,
126                                        final FieldOfView fov,
127                                        final OneAxisEllipsoid body,
128                                        final SphericalPolygonsSet zone,
129                                        final double samplingStep,
130                                        final List<SamplingPoint> sampledZone) {
131 
132         super(detectionSettings, handler);
133         this.fov          = fov;
134         this.body         = body;
135         this.samplingStep = samplingStep;
136         this.zone         = zone;
137         this.sampledZone  = sampledZone;
138 
139         final EnclosingBall<Sphere2D, S2Point> cap = zone.getEnclosingCap();
140         final SinCos sc = FastMath.sinCos(cap.getRadius());
141         this.capCenter    = cap.getCenter().getVector();
142         this.capCos       = sc.cos();
143         this.capSin       = sc.sin();
144 
145     }
146 
147     /** Sample the region.
148      * @param body body on which the geographic zone is defined
149      * @param zone geographic zone to consider
150      * @param samplingStep  linear step used for sampling the geographic zone (in meters)
151      * @return sampling points
152      */
153     private static List<SamplingPoint> sample(final OneAxisEllipsoid body,
154                                               final SphericalPolygonsSet zone,
155                                               final double samplingStep) {
156 
157         final List<SamplingPoint> sampledZone = new ArrayList<>();
158 
159         // sample the zone boundary
160         final List<Vertex> boundary = zone.getBoundaryLoops();
161         for (final Vertex loopStart : boundary) {
162             int count = 0;
163             for (Vertex v = loopStart; count == 0 || v != loopStart; v = v.getOutgoing().getEnd()) {
164                 ++count;
165                 final Edge edge = v.getOutgoing();
166                 final int n = (int) FastMath.ceil(edge.getLength() * body.getEquatorialRadius() / samplingStep);
167                 for (int i = 0; i < n; ++i) {
168                     final S2Point intermediate = new S2Point(edge.getPointAt(i * edge.getLength() / n));
169                     final GeodeticPoint gp = new GeodeticPoint(0.5 * FastMath.PI - intermediate.getPhi(),
170                                                                intermediate.getTheta(), 0.0);
171                     sampledZone.add(new SamplingPoint(body.transform(gp), gp.getZenith()));
172                 }
173             }
174         }
175 
176         // sample the zone interior
177         final EllipsoidTessellator tessellator =
178                         new EllipsoidTessellator(body, new DivertedSingularityAiming(zone), 1);
179         final List<List<GeodeticPoint>> gpSample = tessellator.sample(zone, samplingStep, samplingStep);
180         for (final List<GeodeticPoint> list : gpSample) {
181             for (final GeodeticPoint gp : list) {
182                 sampledZone.add(new SamplingPoint(body.transform(gp), gp.getZenith()));
183             }
184         }
185 
186         return sampledZone;
187 
188     }
189 
190     /** {@inheritDoc} */
191     @Override
192     protected FootprintOverlapDetector create(final EventDetectionSettings detectionSettings,
193                                               final EventHandler newHandler) {
194         return new FootprintOverlapDetector(detectionSettings, newHandler,
195                                             fov, body, zone, samplingStep, sampledZone);
196     }
197 
198     /** Get the geographic zone triggering the events.
199      * <p>
200      * The zone is mapped on the unit sphere
201      * </p>
202      * @return geographic zone triggering the events
203      */
204     public SphericalPolygonsSet getZone() {
205         return zone;
206     }
207 
208     /** Get the Field Of View.
209      * @return Field Of View
210      * @since 10.1
211      */
212     public FieldOfView getFOV() {
213         return fov;
214     }
215 
216     /** Get the body on which the geographic zone is defined.
217      * @return body on which the geographic zone is defined
218      */
219     public BodyShape getBody() {
220         return body;
221     }
222 
223     /** {@inheritDoc}
224      * <p>
225      * The g function value is the minimum offset among the region points
226      * with respect to the Field Of View boundary. It is positive if all region
227      * points are outside of the Field Of View, and negative if at least some
228      * of the region points are inside of the Field Of View. The minimum is
229      * computed by sampling the region, considering only the points for which
230      * the spacecraft is above the horizon. The accuracy of the detection
231      * depends on the linear sampling step set at detector construction. If
232      * the spacecraft is below horizon for all region points, an arbitrary
233      * positive value is returned.
234      * </p>
235      * <p>
236      * As per the previous definition, when the region enters the Field Of
237      * View, a decreasing event is generated, and when the region leaves
238      * the Field Of View, an increasing event is generated.
239      * </p>
240      */
241     public double g(final SpacecraftState s) {
242 
243         // initial arbitrary positive value
244         double value = FastMath.PI;
245 
246         // get spacecraft position in body frame
247         final Vector3D      scBody      = s.getPosition(body.getBodyFrame());
248 
249         // map the point to a sphere
250         final GeodeticPoint gp          = body.transform(scBody, body.getBodyFrame(), s.getDate());
251         final S2Point       s2p         = new S2Point(gp.getLongitude(), 0.5 * FastMath.PI - gp.getLatitude());
252 
253         // for faster computation, we start using only the surrounding cap, to filter out
254         // far away points (which correspond to most of the points if the zone is small)
255         final Vector3D p   = s2p.getVector();
256         final double   dot = Vector3D.dotProduct(p, capCenter);
257         if (dot < capCos) {
258             // the spacecraft is outside of the cap, look for the closest cap point
259             final Vector3D t     = p.subtract(dot, capCenter).normalize();
260             final Vector3D close = new Vector3D(capCos, capCenter, capSin, t);
261             if (Vector3D.dotProduct(p, close) < -0.01) {
262                 // the spacecraft is not visible from the cap edge,
263                 // even taking some margin into account for sphere/ellipsoid different shapes
264                 // this induces no points in the zone can see the spacecraft,
265                 // we can return the arbitrary initial positive value without performing further computation
266                 return value;
267             }
268         }
269 
270         // the spacecraft may be visible from some points in the zone, check them all
271         final StaticTransform bodyToSc = StaticTransform.compose(
272                 s.getDate(),
273                 body.getBodyFrame().getStaticTransformTo(s.getFrame(), s.getDate()),
274                 s.toStaticTransform());
275         for (final SamplingPoint point : sampledZone) {
276             final Vector3D lineOfSightBody = point.getPosition().subtract(scBody);
277             if (Vector3D.dotProduct(lineOfSightBody, point.getZenith()) <= 0) {
278                 // spacecraft is above this sample point local horizon
279                 // get line of sight in spacecraft frame
280                 final double offset = fov.offsetFromBoundary(bodyToSc.transformVector(lineOfSightBody),
281                                                              0.0, VisibilityTrigger.VISIBLE_ONLY_WHEN_FULLY_IN_FOV);
282                 value = FastMath.min(value, offset);
283             }
284         }
285 
286         return value;
287 
288     }
289 
290     /** Container for sampling points. */
291     private static class SamplingPoint {
292 
293         /** Position of the point. */
294         private final Vector3D position;
295 
296         /** Zenith vector of the point. */
297         private final Vector3D zenith;
298 
299         /** Simple constructor.
300          * @param position position of the point
301          * @param zenith zenith vector of the point
302          */
303         SamplingPoint(final Vector3D position, final Vector3D zenith) {
304             this.position = position;
305             this.zenith   = zenith;
306         }
307 
308         /** Get the point position.
309          * @return point position
310          */
311         public Vector3D getPosition() {
312             return position;
313         }
314 
315         /** Get the point zenith vector.
316          * @return point zenith vector
317          */
318         public Vector3D getZenith() {
319             return zenith;
320         }
321 
322     }
323 
324 }