1   /* Copyright 2002-2024 Luc Maisonobe
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  
18  package org.orekit.forces.maneuvers.propulsion;
19  
20  import java.lang.reflect.Array;
21  import java.util.Collections;
22  import java.util.List;
23  import java.util.stream.Stream;
24  
25  import org.hipparchus.CalculusFieldElement;
26  import org.hipparchus.Field;
27  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
28  import org.hipparchus.geometry.euclidean.threed.Vector3D;
29  import org.hipparchus.ode.events.Action;
30  import org.hipparchus.util.FastMath;
31  import org.orekit.forces.maneuvers.Control3DVectorCostType;
32  import org.orekit.propagation.FieldSpacecraftState;
33  import org.orekit.propagation.SpacecraftState;
34  import org.orekit.propagation.events.DateDetector;
35  import org.orekit.propagation.events.EventDetector;
36  import org.orekit.propagation.events.FieldDateDetector;
37  import org.orekit.propagation.events.FieldEventDetector;
38  import org.orekit.time.FieldAbsoluteDate;
39  import org.orekit.time.FieldTimeStamped;
40  import org.orekit.utils.Constants;
41  import org.orekit.utils.ParameterDriver;
42  import org.orekit.utils.TimeSpanMap;
43  
44  /** Thrust propulsion model based on segmented profile.
45   * @author Luc Maisonobe
46   * @since 12.0
47   */
48  public class ProfileThrustPropulsionModel implements ThrustPropulsionModel {
49  
50      /** Accuracy of switching events datation (s). */
51      private static final double DATATION_ACCURACY = 1.0e-10;
52  
53      /** Thrust profile. */
54      private final TimeSpanMap<PolynomialThrustSegment> profile;
55  
56      /** Specific impulse. */
57      private final double isp;
58  
59      /** Name of the maneuver. */
60      private final String name;
61  
62      /** Type of norm linking thrust vector to mass flow rate. */
63      private final Control3DVectorCostType control3DVectorCostType;
64  
65      /** Generic constructor.
66       * @param profile thrust profile (N)
67       * @param isp specific impulse (s)
68       * @param control3DVectorCostType control vector's cost type
69       * @param name name of the maneuver
70       */
71      public ProfileThrustPropulsionModel(final TimeSpanMap<PolynomialThrustSegment> profile, final double isp,
72                                          final Control3DVectorCostType control3DVectorCostType, final String name) {
73          this.name    = name;
74          this.isp     = isp;
75          this.profile = profile;
76          this.control3DVectorCostType = control3DVectorCostType;
77      }
78  
79      /** {@inheritDoc} */
80      @Override
81      public String getName() {
82          return name;
83      }
84  
85      /** {@inheritDoc} */
86      @Override
87      public Control3DVectorCostType getControl3DVectorCostType() {
88          return control3DVectorCostType;
89      }
90  
91      /** {@inheritDoc} */
92      @Override
93      public Vector3D getThrustVector(final SpacecraftState s) {
94          final PolynomialThrustSegment active = profile.get(s.getDate());
95          return active == null ? Vector3D.ZERO : active.getThrustVector(s.getDate());
96      }
97  
98      /** {@inheritDoc} */
99      @Override
100     public double getFlowRate(final SpacecraftState s) {
101         return -control3DVectorCostType.evaluate(getThrustVector(s)) / (Constants.G0_STANDARD_GRAVITY * isp);
102     }
103 
104     /** {@inheritDoc}
105      * <p>
106      * Here the thrust vector does not depend on parameters
107      * </p>
108      */
109     @Override
110     public Vector3D getThrustVector(final SpacecraftState s, final double[] parameters) {
111         return getThrustVector(s);
112     }
113 
114     /** {@inheritDoc}
115      * <p>
116      * Here the flow rate does not depend on parameters
117      * </p>
118      */
119     public double getFlowRate(final SpacecraftState s, final double[] parameters) {
120         return getFlowRate(s);
121     }
122 
123     /** {@inheritDoc}
124      * <p>
125      * Here the thrust vector does not depend on parameters
126      * </p>
127      */
128     public <T extends CalculusFieldElement<T>> FieldVector3D<T> getThrustVector(final FieldSpacecraftState<T> s,
129                                                                                 final T[] parameters) {
130         final PolynomialThrustSegment active = profile.get(s.getDate().toAbsoluteDate());
131         return active == null ? FieldVector3D.getZero(s.getDate().getField()) : active.getThrustVector(s.getDate());
132     }
133 
134     /** {@inheritDoc}
135      * <p>
136      * Here the flow rate does not depend on parameters
137      * </p>
138      */
139     public <T extends CalculusFieldElement<T>> T getFlowRate(final FieldSpacecraftState<T> s, final T[] parameters) {
140         return control3DVectorCostType.evaluate(getThrustVector(s, parameters)).divide(-Constants.G0_STANDARD_GRAVITY * isp);
141     }
142 
143     /** {@inheritDoc}.
144      * <p>
145      * The single detector returned triggers {@link org.hipparchus.ode.events.Action#RESET_DERIVATIVES} events
146      * at every {@link PolynomialThrustSegment thrust segments} boundaries.
147      * </p>
148      */
149     @Override
150     public Stream<EventDetector> getEventDetectors() {
151 
152         final double shortest = shortestSegmentDuration();
153         final DateDetector detector = new DateDetector().
154                                       withMaxCheck(0.5 * shortest).
155                                       withMinGap(0.5 * shortest).
156                                       withThreshold(DATATION_ACCURACY).
157                                       withHandler((state, det, increasing) -> Action.RESET_DERIVATIVES);
158         for (TimeSpanMap.Transition<PolynomialThrustSegment> transition = profile.getFirstTransition();
159              transition != null;
160              transition = transition.next()) {
161             detector.addEventDate(transition.getDate());
162         }
163         return Stream.of(detector);
164     }
165 
166     /** {@inheritDoc}.
167      * <p>
168      * The single detector returned triggers {@link org.hipparchus.ode.events.Action#RESET_DERIVATIVES} events
169      * at every {@link PolynomialThrustSegment thrust segments} boundaries.
170      * </p>
171      */
172     @Override
173     public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventDetectors(final Field<T> field) {
174         final double shortest = shortestSegmentDuration();
175         @SuppressWarnings("unchecked")
176         final FieldDateDetector<T> detector = new FieldDateDetector<>(field,
177                                                                       (FieldTimeStamped<T>[]) Array.newInstance(FieldTimeStamped.class, 0)).
178                                               withMaxCheck(0.5 * shortest).
179                                               withMinGap(0.5 * shortest).
180                                               withThreshold(field.getZero().newInstance(DATATION_ACCURACY)).
181                                               withHandler((state, det, increasing) -> Action.RESET_DERIVATIVES);
182         for (TimeSpanMap.Transition<PolynomialThrustSegment> transition = profile.getFirstTransition();
183              transition != null;
184              transition = transition.next()) {
185             detector.addEventDate(new FieldAbsoluteDate<>(field, transition.getDate()));
186         }
187         return Stream.of(detector);
188     }
189 
190     /** Compute the duration of the shortest segment.
191      * @return duration of the shortest segment
192      */
193     private double shortestSegmentDuration() {
194         double shortest = Double.POSITIVE_INFINITY;
195         for (TimeSpanMap.Span<PolynomialThrustSegment> span = profile.getFirstSpan();
196              span != null;
197              span = span.next()) {
198             shortest = FastMath.min(shortest, span.getEnd().durationFrom(span.getStart()));
199         }
200         return shortest;
201     }
202 
203     @Override
204     public List<ParameterDriver> getParametersDrivers() {
205         return Collections.emptyList();
206     }
207 }