1 /* Copyright 2002-2021 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.Transform;
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(DEFAULT_MAXCHECK, DEFAULT_THRESHOLD, DEFAULT_MAX_ITER,
107 new StopOnIncreasing<FootprintOverlapDetector>(),
108 fov, body, zone, samplingStep, sample(body, zone, samplingStep));
109 }
110
111 /** Private constructor with full parameters.
112 * <p>
113 * This constructor is private as users are expected to use the builder
114 * API with the various {@code withXxx()} methods to set up the instance
115 * in a readable manner without using a huge amount of parameters.
116 * </p>
117 * @param maxCheck maximum checking interval (s)
118 * @param threshold convergence threshold (s)
119 * @param maxIter maximum number of iterations in the event time search
120 * @param handler event handler to call at event occurrences
121 * @param body body on which the geographic zone is defined
122 * @param zone geographic zone to consider
123 * @param fov sensor field of view
124 * @param sampledZone sampling of the geographic zone
125 * @param samplingStep linear step used for sampling the geographic zone (in meters)
126 */
127 private FootprintOverlapDetector(final double maxCheck, final double threshold,
128 final int maxIter, final EventHandler<? super FootprintOverlapDetector> handler,
129 final FieldOfView fov,
130 final OneAxisEllipsoid body,
131 final SphericalPolygonsSet zone,
132 final double samplingStep,
133 final List<SamplingPoint> sampledZone) {
134
135 super(maxCheck, threshold, maxIter, handler);
136 this.fov = fov;
137 this.body = body;
138 this.samplingStep = samplingStep;
139 this.zone = zone;
140 this.sampledZone = sampledZone;
141
142 final EnclosingBall<Sphere2D, S2Point> cap = zone.getEnclosingCap();
143 final SinCos sc = FastMath.sinCos(cap.getRadius());
144 this.capCenter = cap.getCenter().getVector();
145 this.capCos = sc.cos();
146 this.capSin = sc.sin();
147
148 }
149
150 /** Sample the region.
151 * @param body body on which the geographic zone is defined
152 * @param zone geographic zone to consider
153 * @param samplingStep linear step used for sampling the geographic zone (in meters)
154 * @return sampling points
155 */
156 private static List<SamplingPoint> sample(final OneAxisEllipsoid body,
157 final SphericalPolygonsSet zone,
158 final double samplingStep) {
159
160 final List<SamplingPoint> sampledZone = new ArrayList<SamplingPoint>();
161
162 // sample the zone boundary
163 final List<Vertex> boundary = zone.getBoundaryLoops();
164 for (final Vertex loopStart : boundary) {
165 int count = 0;
166 for (Vertex v = loopStart; count == 0 || v != loopStart; v = v.getOutgoing().getEnd()) {
167 ++count;
168 final Edge edge = v.getOutgoing();
169 final int n = (int) FastMath.ceil(edge.getLength() * body.getEquatorialRadius() / samplingStep);
170 for (int i = 0; i < n; ++i) {
171 final S2Point intermediate = new S2Point(edge.getPointAt(i * edge.getLength() / n));
172 final GeodeticPoint gp = new GeodeticPoint(0.5 * FastMath.PI - intermediate.getPhi(),
173 intermediate.getTheta(), 0.0);
174 sampledZone.add(new SamplingPoint(body.transform(gp), gp.getZenith()));
175 }
176 }
177 }
178
179 // sample the zone interior
180 final EllipsoidTessellator tessellator =
181 new EllipsoidTessellator(body, new DivertedSingularityAiming(zone), 1);
182 final List<List<GeodeticPoint>> gpSample = tessellator.sample(zone, samplingStep, samplingStep);
183 for (final List<GeodeticPoint> list : gpSample) {
184 for (final GeodeticPoint gp : list) {
185 sampledZone.add(new SamplingPoint(body.transform(gp), gp.getZenith()));
186 }
187 }
188
189 return sampledZone;
190
191 }
192
193 /** {@inheritDoc} */
194 @Override
195 protected FootprintOverlapDetector create(final double newMaxCheck, final double newThreshold,
196 final int newMaxIter,
197 final EventHandler<? super FootprintOverlapDetector> newHandler) {
198 return new FootprintOverlapDetector(newMaxCheck, newThreshold, newMaxIter, newHandler,
199 fov, body, zone, samplingStep, sampledZone);
200 }
201
202 /** Get the geographic zone triggering the events.
203 * <p>
204 * The zone is mapped on the unit sphere
205 * </p>
206 * @return geographic zone triggering the events
207 */
208 public SphericalPolygonsSet getZone() {
209 return zone;
210 }
211
212 /** Get the Field Of View.
213 * @return Field Of View
214 * @since 10.1
215 */
216 public FieldOfView getFOV() {
217 return fov;
218 }
219
220 /** Get the body on which the geographic zone is defined.
221 * @return body on which the geographic zone is defined
222 */
223 public BodyShape getBody() {
224 return body;
225 }
226
227 /** {@inheritDoc}
228 * <p>
229 * The g function value is the minimum offset among the region points
230 * with respect to the Field Of View boundary. It is positive if all region
231 * points are outside of the Field Of View, and negative if at least some
232 * of the region points are inside of the Field Of View. The minimum is
233 * computed by sampling the region, considering only the points for which
234 * the spacecraft is above the horizon. The accuracy of the detection
235 * depends on the linear sampling step set at detector construction. If
236 * the spacecraft is below horizon for all region points, an arbitrary
237 * positive value is returned.
238 * </p>
239 * <p>
240 * As per the previous definition, when the region enters the Field Of
241 * View, a decreasing event is generated, and when the region leaves
242 * the Field Of View, an increasing event is generated.
243 * </p>
244 */
245 public double g(final SpacecraftState s) {
246
247 // initial arbitrary positive value
248 double value = FastMath.PI;
249
250 // get spacecraft position in body frame
251 final Vector3D scBody = s.getPVCoordinates(body.getBodyFrame()).getPosition();
252
253 // map the point to a sphere
254 final GeodeticPoint gp = body.transform(scBody, body.getBodyFrame(), s.getDate());
255 final S2Point s2p = new S2Point(gp.getLongitude(), 0.5 * FastMath.PI - gp.getLatitude());
256
257 // for faster computation, we start using only the surrounding cap, to filter out
258 // far away points (which correspond to most of the points if the zone is small)
259 final Vector3D p = s2p.getVector();
260 final double dot = Vector3D.dotProduct(p, capCenter);
261 if (dot < capCos) {
262 // the spacecraft is outside of the cap, look for the closest cap point
263 final Vector3D t = p.subtract(dot, capCenter).normalize();
264 final Vector3D close = new Vector3D(capCos, capCenter, capSin, t);
265 if (Vector3D.dotProduct(p, close) < -0.01) {
266 // the spacecraft is not visible from the cap edge,
267 // even taking some margin into account for sphere/ellipsoid different shapes
268 // this induces no points in the zone can see the spacecraft,
269 // we can return the arbitrary initial positive value without performing further computation
270 return value;
271 }
272 }
273
274 // the spacecraft may be visible from some points in the zone, check them all
275 final Transform bodyToSc = new Transform(s.getDate(),
276 body.getBodyFrame().getTransformTo(s.getFrame(), s.getDate()),
277 s.toTransform());
278 for (final SamplingPoint point : sampledZone) {
279 final Vector3D lineOfSightBody = point.getPosition().subtract(scBody);
280 if (Vector3D.dotProduct(lineOfSightBody, point.getZenith()) <= 0) {
281 // spacecraft is above this sample point local horizon
282 // get line of sight in spacecraft frame
283 final double offset = fov.offsetFromBoundary(bodyToSc.transformVector(lineOfSightBody),
284 0.0, VisibilityTrigger.VISIBLE_ONLY_WHEN_FULLY_IN_FOV);
285 value = FastMath.min(value, offset);
286 }
287 }
288
289 return value;
290
291 }
292
293 /** Container for sampling points. */
294 private static class SamplingPoint {
295
296 /** Position of the point. */
297 private final Vector3D position;
298
299 /** Zenith vector of the point. */
300 private final Vector3D zenith;
301
302 /** Simple constructor.
303 * @param position position of the point
304 * @param zenith zenith vector of the point
305 */
306 SamplingPoint(final Vector3D position, final Vector3D zenith) {
307 this.position = position;
308 this.zenith = zenith;
309 }
310
311 /** Get the point position.
312 * @return point position
313 */
314 public Vector3D getPosition() {
315 return position;
316 }
317
318 /** Get the point zenith vector.
319 * @return point zenith vector
320 */
321 public Vector3D getZenith() {
322 return zenith;
323 }
324
325 }
326
327 }