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.function.Function;
20
21 import org.hipparchus.analysis.UnivariateFunction;
22 import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
23 import org.hipparchus.util.FastMath;
24 import org.hipparchus.util.MathUtils;
25 import org.orekit.errors.OrekitIllegalArgumentException;
26 import org.orekit.errors.OrekitMessages;
27 import org.orekit.orbits.CircularOrbit;
28 import org.orekit.orbits.EquinoctialOrbit;
29 import org.orekit.orbits.KeplerianOrbit;
30 import org.orekit.orbits.Orbit;
31 import org.orekit.orbits.OrbitType;
32 import org.orekit.orbits.PositionAngleType;
33 import org.orekit.propagation.SpacecraftState;
34 import org.orekit.propagation.events.handlers.EventHandler;
35 import org.orekit.propagation.events.handlers.StopOnEvent;
36 import org.orekit.time.AbsoluteDate;
37 import org.orekit.utils.TimeSpanMap;
38
39 /** Detector for in-orbit position angle.
40 * <p>
41 * The detector is based on anomaly for {@link OrbitType#KEPLERIAN Keplerian}
42 * orbits, latitude argument for {@link OrbitType#CIRCULAR circular} orbits,
43 * or longitude argument for {@link OrbitType#EQUINOCTIAL equinoctial} orbits.
44 * It does not support {@link OrbitType#CARTESIAN Cartesian} orbits. The
45 * angles can be either {@link PositionAngleType#TRUE true}, {@link PositionAngleType#MEAN
46 * mean} or {@link PositionAngleType#ECCENTRIC eccentric} angles.
47 * </p>
48 * @author Luc Maisonobe
49 * @since 7.1
50 */
51 public class PositionAngleDetector extends AbstractDetector<PositionAngleDetector> {
52
53 /** Orbit type defining the angle type. */
54 private final OrbitType orbitType;
55
56 /** Type of position angle. */
57 private final PositionAngleType positionAngleType;
58
59 /** Fixed angle to be crossed. */
60 private final double angle;
61
62 /** Position angle extraction function. */
63 private final Function<Orbit, Double> positionAngleExtractor;
64
65 /** Estimators for the offset angle, taking care of 2π wrapping and g function continuity. */
66 private TimeSpanMap<OffsetEstimator> offsetEstimators;
67
68 /** Build a new detector.
69 * <p>The new instance uses default values for maximal checking interval
70 * ({@link #DEFAULT_MAX_CHECK}) and convergence threshold ({@link
71 * #DEFAULT_THRESHOLD}).</p>
72 * @param orbitType orbit type defining the angle type
73 * @param positionAngleType type of position angle
74 * @param angle fixed angle to be crossed
75 * @exception OrekitIllegalArgumentException if orbit type is {@link OrbitType#CARTESIAN}
76 */
77 public PositionAngleDetector(final OrbitType orbitType, final PositionAngleType positionAngleType,
78 final double angle)
79 throws OrekitIllegalArgumentException {
80 this(DEFAULT_MAX_CHECK, DEFAULT_THRESHOLD, orbitType, positionAngleType, angle);
81 }
82
83 /** Build a detector.
84 * <p> This instance uses by default the {@link StopOnEvent} handler </p>
85 * @param maxCheck maximal checking interval (s)
86 * @param threshold convergence threshold (s)
87 * @param orbitType orbit type defining the angle type
88 * @param positionAngleType type of position angle
89 * @param angle fixed angle to be crossed
90 * @exception OrekitIllegalArgumentException if orbit type is {@link OrbitType#CARTESIAN}
91 */
92 public PositionAngleDetector(final double maxCheck, final double threshold,
93 final OrbitType orbitType, final PositionAngleType positionAngleType,
94 final double angle)
95 throws OrekitIllegalArgumentException {
96 this(new EventDetectionSettings(maxCheck, threshold, DEFAULT_MAX_ITER), new StopOnEvent(),
97 orbitType, positionAngleType, angle);
98 }
99
100 /** Protected constructor with full parameters.
101 * <p>
102 * This constructor is not public as users are expected to use the builder
103 * API with the various {@code withXxx()} methods to set up the instance
104 * in a readable manner without using a huge amount of parameters.
105 * </p>
106 * @param detectionSettings event detection settings
107 * @param handler event handler to call at event occurrences
108 * @param orbitType orbit type defining the angle type
109 * @param positionAngleType type of position angle
110 * @param angle fixed angle to be crossed
111 * @exception OrekitIllegalArgumentException if orbit type is {@link OrbitType#CARTESIAN}
112 * @since 13.0
113 */
114 protected PositionAngleDetector(final EventDetectionSettings detectionSettings, final EventHandler handler,
115 final OrbitType orbitType, final PositionAngleType positionAngleType,
116 final double angle)
117 throws OrekitIllegalArgumentException {
118
119 super(detectionSettings, handler);
120
121 this.orbitType = orbitType;
122 this.positionAngleType = positionAngleType;
123 this.angle = angle;
124 this.offsetEstimators = null;
125
126 switch (orbitType) {
127 case KEPLERIAN:
128 positionAngleExtractor = o -> ((KeplerianOrbit) orbitType.convertType(o)).getAnomaly(positionAngleType);
129 break;
130 case CIRCULAR:
131 positionAngleExtractor = o -> ((CircularOrbit) orbitType.convertType(o)).getAlpha(positionAngleType);
132 break;
133 case EQUINOCTIAL:
134 positionAngleExtractor = o -> ((EquinoctialOrbit) orbitType.convertType(o)).getL(positionAngleType);
135 break;
136 default:
137 final String sep = ", ";
138 throw new OrekitIllegalArgumentException(OrekitMessages.ORBIT_TYPE_NOT_ALLOWED,
139 orbitType,
140 OrbitType.KEPLERIAN + sep +
141 OrbitType.CIRCULAR + sep +
142 OrbitType.EQUINOCTIAL);
143 }
144
145 }
146
147 /** {@inheritDoc} */
148 @Override
149 protected PositionAngleDetector create(final EventDetectionSettings detectionSettings,
150 final EventHandler newHandler) {
151 return new PositionAngleDetector(detectionSettings, newHandler, orbitType, positionAngleType, angle);
152 }
153
154 /** Get the orbit type defining the angle type.
155 * @return orbit type defining the angle type
156 */
157 public OrbitType getOrbitType() {
158 return orbitType;
159 }
160
161 /** Get the type of position angle.
162 * @return type of position angle
163 */
164 public PositionAngleType getPositionAngleType() {
165 return positionAngleType;
166 }
167
168 /** Get the fixed angle to be crossed (radians).
169 * @return fixed angle to be crossed (radians)
170 */
171 public double getAngle() {
172 return angle;
173 }
174
175 /** {@inheritDoc} */
176 @Override
177 public void init(final SpacecraftState s0, final AbsoluteDate t) {
178 super.init(s0, t);
179 offsetEstimators = new TimeSpanMap<>(new OffsetEstimator(s0.getOrbit(), +1.0));
180 }
181
182 /** Compute the value of the detection function.
183 * <p>
184 * The value is the angle difference between the spacecraft and the fixed
185 * angle to be crossed, with some sign tweaks to ensure continuity.
186 * These tweaks imply the {@code increasing} flag in events detection becomes
187 * irrelevant here! As an example, the angle always increase in a Keplerian
188 * orbit, but this g function will increase and decrease so it
189 * will cross the zero value once per orbit, in increasing and decreasing
190 * directions on alternate orbits..
191 * </p>
192 * @param s the current state information: date, kinematics, attitude
193 * @return angle difference between the spacecraft and the fixed
194 * angle, with some sign tweaks to ensure continuity
195 */
196 public double g(final SpacecraftState s) {
197
198 final Orbit orbit = s.getOrbit();
199
200 // angle difference
201 OffsetEstimator estimator = offsetEstimators.get(s.getDate());
202 double delta = estimator.delta(orbit);
203
204 // we use a value greater than π for handover in order to avoid
205 // several switches to be estimated as the calling propagator
206 // and Orbit.shiftedBy have different accuracy. It is sufficient
207 // to have a handover roughly opposite to the detected position angle
208 while (FastMath.abs(delta) >= 3.5) {
209 // we are too far away from the current estimator, we need to set up a new one
210 // ensuring that we do have a crossing event in the current orbit
211 // and we ensure sign continuity with the current estimator
212
213 // find when the previous estimator becomes invalid
214 final AbsoluteDate handover = estimator.dateForOffset(FastMath.copySign(FastMath.PI, delta), orbit);
215
216 // perform handover to a new estimator at this date
217 estimator = new OffsetEstimator(orbit, delta);
218 delta = estimator.delta(orbit);
219 if (isForward()) {
220 offsetEstimators.addValidAfter(estimator, handover.getDate(), false);
221 } else {
222 offsetEstimators.addValidBefore(estimator, handover.getDate(), false);
223 }
224
225 }
226
227 return delta;
228
229 }
230
231 /** Local class for estimating offset angle, handling 2π wrap-up and sign continuity. */
232 private class OffsetEstimator {
233
234 /** Target angle. */
235 private final double target;
236
237 /** Sign correction to offset. */
238 private final double sign;
239
240 /** Reference angle. */
241 private final double r0;
242
243 /** Slope of the linearized model. */
244 private final double r1;
245
246 /** Reference date. */
247 private final AbsoluteDate t0;
248
249 /** Simple constructor.
250 * @param orbit current orbit
251 * @param currentSign desired sign of the offset at current orbit time (magnitude is ignored)
252 */
253 OffsetEstimator(final Orbit orbit, final double currentSign) {
254 r0 = positionAngleExtractor.apply(orbit);
255 target = MathUtils.normalizeAngle(angle, r0);
256 sign = FastMath.copySign(1.0, (r0 - target) * currentSign);
257 r1 = orbit.getKeplerianMeanMotion();
258 t0 = orbit.getDate();
259 }
260
261 /** Compute offset from reference angle.
262 * @param orbit current orbit
263 * @return offset between current angle and reference angle
264 */
265 public double delta(final Orbit orbit) {
266 final double rawAngle = positionAngleExtractor.apply(orbit);
267 final double linearReference = r0 + r1 * orbit.getDate().durationFrom(t0);
268 final double linearizedAngle = MathUtils.normalizeAngle(rawAngle, linearReference);
269 return sign * (linearizedAngle - target);
270 }
271
272 /** Find date at which offset reaches specified value.
273 * <p>
274 * This computation is an approximation because it relies on
275 * {@link Orbit#shiftedBy(double)} only.
276 * </p>
277 * @param offset target value for offset angle
278 * @param orbit current orbit
279 * @return approximate date at which offset reached specified value
280 */
281 public AbsoluteDate dateForOffset(final double offset, final Orbit orbit) {
282
283 // bracket the search
284 final double period = orbit.getKeplerianPeriod();
285 final double delta0 = delta(orbit);
286 final double searchInf;
287 final double searchSup;
288 if ((delta0 - offset) * sign >= 0) {
289 // the date is before current orbit
290 searchInf = -period;
291 searchSup = 0;
292 } else {
293 // the date is after current orbit
294 searchInf = 0;
295 searchSup = +period;
296 }
297
298 // find the date as an offset from current orbit
299 final BracketingNthOrderBrentSolver solver = new BracketingNthOrderBrentSolver(getThreshold(), 5);
300 final UnivariateFunction f = dt -> delta(orbit.shiftedBy(dt)) - offset;
301 final double root = solver.solve(getMaxIterationCount(), f, searchInf, searchSup);
302
303 return orbit.getDate().shiftedBy(root);
304
305 }
306
307 }
308
309 }