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.CalculusFieldElement;
20 import org.hipparchus.ode.events.Action;
21 import org.hipparchus.util.FastMath;
22 import org.hipparchus.util.MathUtils;
23 import org.orekit.frames.Frame;
24 import org.orekit.orbits.FieldOrbit;
25 import org.orekit.orbits.KeplerianOrbit;
26 import org.orekit.orbits.Orbit;
27 import org.orekit.orbits.OrbitType;
28 import org.orekit.orbits.PositionAngleType;
29 import org.orekit.propagation.FieldSpacecraftState;
30 import org.orekit.propagation.events.handlers.FieldEventHandler;
31 import org.orekit.propagation.events.handlers.FieldStopOnIncreasing;
32 import org.orekit.propagation.events.intervals.FieldAdaptableInterval;
33
34 /** Finder for node crossing events.
35 * <p>This class finds equator crossing events (i.e. ascending
36 * or descending node crossing).</p>
37 * <p>The default implementation behavior is to {@link Action#CONTINUE continue}
38 * propagation at descending node crossing and to {@link Action#STOP stop} propagation
39 * at ascending node crossing. This can be changed by calling
40 * {@link #withHandler(FieldEventHandler)} after construction.</p>
41 * <p>Beware that node detection will fail for almost equatorial orbits. If
42 * for example a node detector is used to trigger an {@link
43 * org.orekit.forces.maneuvers.ImpulseManeuver ImpulseManeuver} and the maneuver
44 * turn the orbit plane to equator, then the detector may completely fail just
45 * after the maneuver has been performed! This is a real case that has been
46 * encountered during validation ...</p>
47 * @see org.orekit.propagation.FieldPropagator#addEventDetector(FieldEventDetector)
48 * @author Luc Maisonobe
49 * @param <T> type of the field elements
50 */
51 public class FieldNodeDetector<T extends CalculusFieldElement<T>> extends FieldAbstractDetector<FieldNodeDetector<T>, T> {
52
53 /** Frame in which the equator is defined. */
54 private final Frame frame;
55
56 /** Build a new instance.
57 * <p>The orbit is used only to set an upper bound for the max check interval
58 * to period/3 and to set the convergence threshold according to orbit size.</p>
59 * @param orbit initial orbit
60 * @param frame frame in which the equator is defined (typical
61 * values are {@link org.orekit.frames.FramesFactory#getEME2000() EME<sub>2000</sub>} or
62 * {@link org.orekit.frames.FramesFactory#getITRF(org.orekit.utils.IERSConventions, boolean) ITRF})
63 */
64 public FieldNodeDetector(final FieldOrbit<T> orbit, final Frame frame) {
65 this(orbit.getKeplerianPeriod().multiply(1.0e-13), orbit, frame);
66 }
67
68 /** Build a new instance.
69 * <p>The orbit is used only to set an upper bound for the max check interval
70 * to period/3.</p>
71 * @param threshold convergence threshold (s)
72 * @param orbit initial orbit
73 * @param frame frame in which the equator is defined (typical
74 * values are {@link org.orekit.frames.FramesFactory#getEME2000() EME<sub>2000</sub>} or
75 * {@link org.orekit.frames.FramesFactory#getITRF(org.orekit.utils.IERSConventions, boolean) ITRF})
76 */
77 public FieldNodeDetector(final T threshold, final FieldOrbit<T> orbit, final Frame frame) {
78 this(new FieldEventDetectionSettings<>(FieldAdaptableInterval.of(orbit.getA().getField().getZero().newInstance(2 * estimateNodesTimeSeparation(orbit.toOrbit()) / 3).getReal()),
79 threshold, DEFAULT_MAX_ITER), new FieldStopOnIncreasing<>(), frame);
80 }
81
82 /** Protected constructor with full parameters.
83 * <p>
84 * This constructor is not public as users are expected to use the builder
85 * API with the various {@code withXxx()} methods to set up the instance
86 * in a readable manner without using a huge amount of parameters.
87 * </p>
88 * @param detectionSettings event detection settings
89 * @param handler event handler to call at event occurrences
90 * @param frame frame in which the equator is defined (typical
91 * values are {@link org.orekit.frames.FramesFactory#getEME2000() EME<sub>2000</sub>} or
92 * {@link org.orekit.frames.FramesFactory#getITRF(org.orekit.utils.IERSConventions, boolean) ITRF})
93 * @since 13.0
94 */
95 protected FieldNodeDetector(final FieldEventDetectionSettings<T> detectionSettings,
96 final FieldEventHandler<T> handler, final Frame frame) {
97 super(detectionSettings, handler);
98 this.frame = frame;
99 }
100
101 /** {@inheritDoc} */
102 @Override
103 protected FieldNodeDetector<T> create(final FieldEventDetectionSettings<T> detectionSettings,
104 final FieldEventHandler<T> newHandler) {
105 return new FieldNodeDetector<>(detectionSettings, newHandler, frame);
106 }
107
108 /** Find time separation between nodes.
109 * <p>
110 * The estimation of time separation is based on Keplerian motion, it is only
111 * used as a rough guess for a safe setting of default max check interval for
112 * event detection.
113 * </p>
114 * @param orbit initial orbit
115 * @return minimum time separation between nodes
116 */
117 private static double estimateNodesTimeSeparation(final Orbit orbit) {
118
119 final KeplerianOrbit keplerian = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(orbit);
120
121 // mean anomaly of ascending node
122 final double ascendingM = new KeplerianOrbit(keplerian.getA(), keplerian.getE(),
123 keplerian.getI(),
124 keplerian.getPerigeeArgument(),
125 keplerian.getRightAscensionOfAscendingNode(),
126 -keplerian.getPerigeeArgument(), PositionAngleType.TRUE,
127 keplerian.getFrame(), keplerian.getDate(),
128 keplerian.getMu()).getMeanAnomaly();
129
130 // mean anomaly of descending node
131 final double descendingM = new KeplerianOrbit(keplerian.getA(), keplerian.getE(),
132 keplerian.getI(),
133 keplerian.getPerigeeArgument(),
134 keplerian.getRightAscensionOfAscendingNode(),
135 FastMath.PI - keplerian.getPerigeeArgument(), PositionAngleType.TRUE,
136 keplerian.getFrame(), keplerian.getDate(),
137 keplerian.getMu()).getMeanAnomaly();
138
139 // differences between mean anomalies
140 final double delta1 = MathUtils.normalizeAngle(ascendingM, descendingM + FastMath.PI) - descendingM;
141 final double delta2 = 2 * FastMath.PI - delta1;
142
143 // minimum time separation between the two nodes
144 return FastMath.min(delta1, delta2) / keplerian.getKeplerianMeanMotion();
145
146 }
147
148 /** Get the frame in which the equator is defined.
149 * @return the frame in which the equator is defined
150 */
151 public Frame getFrame() {
152 return frame;
153 }
154
155 /** Compute the value of the switching function.
156 * This function computes the Z position in the defined frame.
157 * @param s the current state information: date, kinematics, attitude
158 * @return value of the switching function
159 */
160 public T g(final FieldSpacecraftState<T> s) {
161 return s.getPosition(frame).getZ();
162 }
163
164 // public NodeDetector toNoField() {
165 // return new NodeDetector(getThreshold().getReal(), orbit.toOrbit(), frame);
166 // }
167
168 }