1   /* Copyright 2002-2024 Joseph Reed
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    * Joseph Reed 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.geometry.euclidean.threed.Vector3D;
20  import org.hipparchus.ode.events.Action;
21  import org.hipparchus.util.MathUtils;
22  import org.orekit.annotation.DefaultDataContext;
23  import org.orekit.bodies.CelestialBodyFactory;
24  import org.orekit.frames.Frame;
25  import org.orekit.frames.FramesFactory;
26  import org.orekit.propagation.SpacecraftState;
27  import org.orekit.propagation.events.handlers.EventHandler;
28  import org.orekit.propagation.events.handlers.StopOnEvent;
29  import org.orekit.utils.PVCoordinatesProvider;
30  import org.orekit.utils.TimeStampedPVCoordinates;
31  
32  /** Finder for beta angle crossing events.
33   * <p>Locate events when the beta angle (the angle between the orbit plane and the celestial body)
34   * crosses a threshold. The {@link #g(SpacecraftState)} function is negative when the beta angle
35   * is above the threshold and positive when the beta angle is below the threshold.</p>
36   * <p>The inertial frame provided must have it's origin centered at the satellite's orbit plane. The
37   * beta angle is computed as the angle between the celestial body's position in this frame with the
38   * satellite's orbital momentum vector.</p>
39   * <p>The default implementation behavior is to {@link Action#STOP stop}
40   * propagation at the first event date occurrence. This can be changed by calling
41   * {@link #withHandler(EventHandler)} after construction.</p>
42   * @see org.orekit.propagation.Propagator#addEventDetector(EventDetector)
43   * @author Joe Reed
44   * @since 12.1
45   */
46  public class BetaAngleDetector extends AbstractDetector<BetaAngleDetector> {
47      /** Beta angle crossing threshold. */
48      private final double betaAngleThreshold;
49      /** Coordinate provider for the celestial body. */
50      private final PVCoordinatesProvider celestialBodyProvider;
51      /** Inertial frame in which beta angle is calculated. */
52      private final Frame inertialFrame;
53  
54      /**Solar beta angle constructor.
55       * <p>This method uses the default data context, assigns the sun as the celestial
56       * body and uses GCRF as the inertial frame.</p>
57       * @param betaAngleThreshold beta angle threshold (radians)
58       */
59      @DefaultDataContext
60      public BetaAngleDetector(final double betaAngleThreshold) {
61          this(betaAngleThreshold, CelestialBodyFactory.getSun(), FramesFactory.getGCRF());
62      }
63  
64      /** Class constructor.
65       * @param betaAngleThreshold beta angle threshold (radians)
66       * @param celestialBodyProvider coordinate provider for the celestial provider
67       * @param inertialFrame inertial frame in which to compute the beta angle
68       */
69      public BetaAngleDetector(final double betaAngleThreshold, final PVCoordinatesProvider celestialBodyProvider,
70              final Frame inertialFrame) {
71          this(AdaptableInterval.of(DEFAULT_MAXCHECK), DEFAULT_THRESHOLD, DEFAULT_MAX_ITER, new StopOnEvent(),
72                  betaAngleThreshold, celestialBodyProvider, inertialFrame);
73      }
74  
75      /** Protected constructor with full parameters.
76       * <p>This constructor is not public as users are expected to use the builder
77       * API with the various {@code withXxx()} methods to set up the instance
78       * in a readable manner without using a huge amount of parameters.</p>
79       * @param maxCheck maximum checking interval
80       * @param threshold convergence threshold (s)
81       * @param maxIter maximum number of iterations in the event time search
82       * @param handler event handler to call at event occurrences
83       * @param betaAngleThreshold beta angle threshold (radians)
84       * @param celestialBodyProvider coordinate provider for the celestial provider
85       * @param inertialFrame inertial frame in which to compute the beta angle
86       */
87      protected BetaAngleDetector(final AdaptableInterval maxCheck, final double threshold,
88                               final int maxIter, final EventHandler handler,
89                               final double betaAngleThreshold, final PVCoordinatesProvider celestialBodyProvider,
90                               final Frame inertialFrame) {
91          super(maxCheck, threshold, maxIter, handler);
92          this.betaAngleThreshold = betaAngleThreshold;
93          this.celestialBodyProvider = celestialBodyProvider;
94          this.inertialFrame = inertialFrame;
95      }
96  
97      /** Coordinate provider for the celestial body.
98       * @return celestial body's coordinate provider
99       */
100     public PVCoordinatesProvider getCelestialBodyProvider() {
101         return this.celestialBodyProvider;
102     }
103 
104     /** The inertial frame in which beta angle is computed.
105      * @return the inertial frame
106      */
107     public Frame getInertialFrame() {
108         return this.inertialFrame;
109     }
110 
111     /** The beta angle threshold (radians).
112      * @return the beta angle threshold (radians)
113      */
114     public double getBetaAngleThreshold() {
115         return this.betaAngleThreshold;
116     }
117 
118     /** Create a new instance with the provided coordinate provider.
119      * <p>This method does not change the current instance.</p>
120      * @param newProvider the new coordinate provider
121      * @return the new detector instance
122      */
123     public BetaAngleDetector withCelestialProvider(final PVCoordinatesProvider newProvider) {
124         return new BetaAngleDetector(getMaxCheckInterval(), getThreshold(), getMaxIterationCount(),
125                 getHandler(), getBetaAngleThreshold(), newProvider, getInertialFrame());
126     }
127 
128     /** Create a new instance with the provided beta angle threshold.
129      * <p>This method does not change the current instance.</p>
130      * @param newBetaAngleThreshold the beta angle threshold (radians)
131      * @return the new detector instance
132      */
133     public BetaAngleDetector withBetaThreshold(final double newBetaAngleThreshold) {
134         return new BetaAngleDetector(getMaxCheckInterval(), getThreshold(), getMaxIterationCount(),
135                 getHandler(), newBetaAngleThreshold, getCelestialBodyProvider(), getInertialFrame());
136     }
137 
138     /** Create a new instance with the provided inertial frame.
139      * <p>This method does not change the current instance.</p>
140      * @param newFrame the inertial frame
141      * @return the new detector instance
142      */
143     public BetaAngleDetector withInertialFrame(final Frame newFrame) {
144         return new BetaAngleDetector(getMaxCheckInterval(), getThreshold(), getMaxIterationCount(),
145                 getHandler(), getBetaAngleThreshold(), getCelestialBodyProvider(), newFrame);
146     }
147 
148     /** {@inheritDoc} */
149     @Override
150     public double g(final SpacecraftState s) {
151         final double beta = calculateBetaAngle(s, celestialBodyProvider, inertialFrame);
152         return betaAngleThreshold - beta;
153     }
154 
155     /**Calculate the beta angle between the orbit plane and the celestial body.
156      * <p>This method computes the beta angle using the frame from the spacecraft state.</p>
157      * @param state spacecraft state
158      * @param celestialBodyProvider celestial body coordinate provider
159      * @return the beta angle (radians)
160      */
161     public static double calculateBetaAngle(final SpacecraftState state,
162             final PVCoordinatesProvider celestialBodyProvider) {
163         return calculateBetaAngle(state, celestialBodyProvider, state.getFrame());
164     }
165 
166     /**Calculate the beta angle between the orbit plane and the celestial body.
167      * @param state spacecraft state
168      * @param celestialBodyProvider celestial body coordinate provider
169      * @param frame inertial frame in which beta angle will be computed
170      * @return the beta angle (radians)
171      */
172     public static double calculateBetaAngle(final SpacecraftState state,
173             final PVCoordinatesProvider celestialBodyProvider, final Frame frame) {
174         final Vector3D celestialP = celestialBodyProvider.getPosition(state.getDate(), frame);
175         final TimeStampedPVCoordinates pv = state.getPVCoordinates(frame);
176         return MathUtils.SEMI_PI - Vector3D.angle(celestialP, pv.getMomentum());
177     }
178 
179     /** {@inheritDoc} */
180     @Override
181     protected BetaAngleDetector create(final AdaptableInterval newMaxCheck, final double newThreshold,
182             final int newMaxIter, final EventHandler newHandler) {
183         return new BetaAngleDetector(newMaxCheck, newThreshold, newMaxIter, newHandler,
184                 getBetaAngleThreshold(), getCelestialBodyProvider(), getInertialFrame());
185     }
186 }