1   /* Copyright 2002-2024 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.integration;
18  
19  import java.util.ArrayList;
20  import java.util.List;
21  
22  import org.hipparchus.Field;
23  import org.hipparchus.analysis.differentiation.Gradient;
24  import org.hipparchus.analysis.differentiation.GradientField;
25  import org.hipparchus.geometry.euclidean.threed.FieldRotation;
26  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
27  import org.hipparchus.geometry.euclidean.threed.Vector3D;
28  import org.orekit.attitudes.AttitudeProvider;
29  import org.orekit.attitudes.FieldAttitude;
30  import org.orekit.orbits.OrbitType;
31  import org.orekit.orbits.FieldCartesianOrbit;
32  import org.orekit.orbits.FieldEquinoctialOrbit;
33  import org.orekit.orbits.FieldOrbit;
34  import org.orekit.orbits.PositionAngleType;
35  import org.orekit.propagation.FieldSpacecraftState;
36  import org.orekit.propagation.SpacecraftState;
37  import org.orekit.time.AbsoluteDate;
38  import org.orekit.time.FieldAbsoluteDate;
39  import org.orekit.utils.FieldAngularCoordinates;
40  import org.orekit.utils.FieldPVCoordinates;
41  import org.orekit.utils.FieldAbsolutePVCoordinates;
42  import org.orekit.utils.ParameterDriver;
43  import org.orekit.utils.ParameterDriversProvider;
44  import org.orekit.utils.TimeStampedFieldAngularCoordinates;
45  import org.orekit.utils.TimeStampedFieldPVCoordinates;
46  import org.orekit.utils.TimeSpanMap.Span;
47  
48  /** Converter for states and parameters arrays.
49   *  @author Luc Maisonobe
50   *  @author Bryan Cazabonne
51   *  @since 10.2
52   */
53  public abstract class AbstractGradientConverter {
54  
55      /** Dimension of the state. */
56      private final int freeStateParameters;
57  
58      /** States with various number of additional parameters. */
59      private final List<FieldSpacecraftState<Gradient>> gStates;
60  
61      /** Simple constructor.
62       * @param freeStateParameters number of free parameters
63       */
64      protected AbstractGradientConverter(final int freeStateParameters) {
65          this.freeStateParameters = freeStateParameters;
66          this.gStates             = new ArrayList<>();
67      }
68  
69      /** Get the number of free state parameters.
70       * @return number of free state parameters
71       */
72      public int getFreeStateParameters() {
73          return freeStateParameters;
74      }
75  
76      /** Initialize first state with 0 parameters.
77       * @param zeroParametersState state with zero parameters
78       * @since 11.2
79       */
80      protected void initStates(final FieldSpacecraftState<Gradient> zeroParametersState) {
81          gStates.clear();
82          gStates.add(zeroParametersState);
83      }
84  
85      /** Add zero derivatives.
86       * @param original original scalar
87       * @param freeParameters total number of free parameters in the gradient
88       * @return extended scalar
89       */
90      protected Gradient extend(final Gradient original, final int freeParameters) {
91          final double[] originalDerivatives = original.getGradient();
92          final double[] extendedDerivatives = new double[freeParameters];
93          System.arraycopy(originalDerivatives, 0, extendedDerivatives, 0, originalDerivatives.length);
94          return new Gradient(original.getValue(), extendedDerivatives);
95      }
96  
97      /**
98       * Add zero derivatives.
99       *
100      * @param original       original date
101      * @param freeParameters total number of free parameters in the gradient
102      * @return extended date
103      */
104     protected FieldAbsoluteDate<Gradient> extend(
105             final FieldAbsoluteDate<Gradient> original,
106             final int freeParameters) {
107         final AbsoluteDate date = original.toAbsoluteDate();
108         final Gradient gradient = original.durationFrom(date);
109         return new FieldAbsoluteDate<>(date, extend(gradient, freeParameters));
110     }
111 
112     /** Add zero derivatives.
113      * @param original original vector
114      * @param freeParameters total number of free parameters in the gradient
115      * @return extended vector
116      */
117     protected FieldVector3D<Gradient> extend(final FieldVector3D<Gradient> original, final int freeParameters) {
118         return new FieldVector3D<>(extend(original.getX(), freeParameters),
119                                    extend(original.getY(), freeParameters),
120                                    extend(original.getZ(), freeParameters));
121     }
122 
123     /** Add zero derivatives.
124      * @param original original rotation
125      * @param freeParameters total number of free parameters in the gradient
126      * @return extended rotation
127      */
128     protected FieldRotation<Gradient> extend(final FieldRotation<Gradient> original, final int freeParameters) {
129         return new FieldRotation<>(extend(original.getQ0(), freeParameters),
130                                    extend(original.getQ1(), freeParameters),
131                                    extend(original.getQ2(), freeParameters),
132                                    extend(original.getQ3(), freeParameters),
133                                    false);
134     }
135 
136     /** Process a state into a Gradient version without force model parameter.
137      * @param state state
138      * @param freeStateParameters number of free parameters
139      * @param provider attitude provider
140      * @return Gradient version of the state
141      * @since 12.0
142      */
143     protected static FieldSpacecraftState<Gradient> buildBasicGradientSpacecraftState(final SpacecraftState state,
144                                                                                       final int freeStateParameters,
145                                                                                       final AttitudeProvider provider) {
146 
147         // Derivative field
148         final Field<Gradient> field =  GradientField.getField(freeStateParameters);
149 
150         // position always has derivatives
151         final Vector3D pos = state.getPosition();
152         final FieldVector3D<Gradient> posG = new FieldVector3D<>(
153                 Gradient.variable(freeStateParameters, 0, pos.getX()),
154                 Gradient.variable(freeStateParameters, 1, pos.getY()),
155                 Gradient.variable(freeStateParameters, 2, pos.getZ()));
156 
157         // velocity may have derivatives or not
158         final Vector3D vel = state.getPVCoordinates().getVelocity();
159         final FieldVector3D<Gradient> velG;
160         if (freeStateParameters > 3) {
161             velG = new FieldVector3D<>(
162                     Gradient.variable(freeStateParameters, 3, vel.getX()),
163                     Gradient.variable(freeStateParameters, 4, vel.getY()),
164                     Gradient.variable(freeStateParameters, 5, vel.getZ()));
165         } else {
166             velG = new FieldVector3D<>(field, vel);
167         }
168 
169         // acceleration never has derivatives
170         final Vector3D acc = state.getPVCoordinates().getAcceleration();
171         final FieldVector3D<Gradient> accG = new FieldVector3D<>(field, acc);
172 
173         // mass never has derivatives
174         final Gradient gMass = Gradient.constant(freeStateParameters, state.getMass());
175 
176         final TimeStampedFieldPVCoordinates<Gradient> timeStampedFieldPVCoordinates = new TimeStampedFieldPVCoordinates<>(
177                 state.getDate(), posG, velG, accG);
178 
179         final FieldCartesianOrbit<Gradient> gOrbit;
180         final FieldAbsolutePVCoordinates<Gradient> gAbsolutePV;
181         if (state.isOrbitDefined()) {
182             final Gradient gMu = Gradient.constant(freeStateParameters, state.getMu());
183             gOrbit = new FieldCartesianOrbit<>(timeStampedFieldPVCoordinates, state.getFrame(), gMu);
184             gAbsolutePV = null;
185         } else {
186             gOrbit = null;
187             gAbsolutePV = new FieldAbsolutePVCoordinates<>(state.getFrame(), timeStampedFieldPVCoordinates);
188         }
189 
190         final FieldAttitude<Gradient> gAttitude;
191         if (freeStateParameters > 3) {
192             // compute attitude partial derivatives with respect to position/velocity
193             gAttitude = provider.getAttitude((state.isOrbitDefined()) ? gOrbit : gAbsolutePV,
194                     timeStampedFieldPVCoordinates.getDate(), state.getFrame());
195         } else {
196             // force model does not depend on attitude, don't bother recomputing it
197             gAttitude = new FieldAttitude<>(field, state.getAttitude());
198         }
199 
200         if (state.isOrbitDefined()) {
201             return new FieldSpacecraftState<>(gOrbit, gAttitude, gMass);
202         } else {
203             return new FieldSpacecraftState<>(gAbsolutePV, gAttitude, gMass);
204         }
205     }
206 
207     /**
208      * Get the state with the number of parameters consistent with parametric model.
209      * @param parametricModel parametric model
210      * @return state with the number of parameters consistent with parametric model
211      */
212     public FieldSpacecraftState<Gradient> getState(final ParameterDriversProvider parametricModel) {
213 
214         // count the required number of parameters
215         int nbParams = 0;
216         for (final ParameterDriver driver : parametricModel.getParametersDrivers()) {
217             if (driver.isSelected()) {
218                 nbParams += driver.getNbOfValues();
219             }
220         }
221 
222         // fill in intermediate slots
223         while (gStates.size() < nbParams + 1) {
224             gStates.add(null);
225         }
226 
227         if (gStates.get(nbParams) == null) {
228             // it is the first time we need this number of parameters
229             // we need to create the state
230             final int freeParameters = freeStateParameters + nbParams;
231             final FieldSpacecraftState<Gradient> s0 = gStates.get(0);
232             final AbsoluteDate date = s0.getDate().toAbsoluteDate();
233 
234             // attitude
235             final FieldAngularCoordinates<Gradient> ac0 = s0.getAttitude().getOrientation();
236             final FieldAttitude<Gradient> gAttitude =
237                     new FieldAttitude<>(s0.getAttitude().getReferenceFrame(),
238                             new TimeStampedFieldAngularCoordinates<>(date,
239                                     extend(ac0.getRotation(), freeParameters),
240                                     extend(ac0.getRotationRate(), freeParameters),
241                                     extend(ac0.getRotationAcceleration(), freeParameters)));
242 
243             // mass
244             final Gradient gMass = extend(s0.getMass(), freeParameters);
245 
246             // orbit or absolute position-velocity coordinates
247             final FieldPVCoordinates<Gradient> pv0 = s0.getPVCoordinates();
248             final TimeStampedFieldPVCoordinates<Gradient> timeStampedFieldPVCoordinates = new TimeStampedFieldPVCoordinates<>(
249                     date,
250                     extend(pv0.getPosition(),     freeParameters),
251                     extend(pv0.getVelocity(),     freeParameters),
252                     extend(pv0.getAcceleration(), freeParameters));
253             final FieldSpacecraftState<Gradient> spacecraftState;
254             if (s0.isOrbitDefined()) {
255                 final FieldOrbit<Gradient> orbit = s0.getOrbit();
256                 if (orbit.getType().equals(OrbitType.EQUINOCTIAL)) {
257                     // for DSST, which always uses EquinoctialOrbit, not CartesianOrbit
258                     spacecraftState = new FieldSpacecraftState<>(
259                             new FieldEquinoctialOrbit<>(
260                                     extend(orbit.getA(), freeParameters),
261                                     extend(orbit.getEquinoctialEx(), freeParameters),
262                                     extend(orbit.getEquinoctialEy(), freeParameters),
263                                     extend(orbit.getHx(), freeParameters),
264                                     extend(orbit.getHy(), freeParameters),
265                                     extend(orbit.getLM(), freeParameters),
266                                     PositionAngleType.MEAN,
267                                     s0.getFrame(),
268                                     extend(s0.getDate(), freeParameters),
269                                     extend(s0.getMu(), freeParameters)
270                             ),
271                             gAttitude,
272                             gMass);
273                 } else {
274                     spacecraftState = new FieldSpacecraftState<>(new FieldCartesianOrbit<>(timeStampedFieldPVCoordinates,
275                             s0.getFrame(), extend(s0.getMu(), freeParameters)), gAttitude, gMass);
276                 }
277             } else {
278                 spacecraftState = new FieldSpacecraftState<>(new FieldAbsolutePVCoordinates<>(s0.getFrame(),
279                         timeStampedFieldPVCoordinates), gAttitude, gMass);
280             }
281 
282             gStates.set(nbParams, spacecraftState);
283 
284         }
285 
286         return gStates.get(nbParams);
287 
288     }
289 
290     /** Get the parametric model parameters, return gradient values for each span of each driver (several gradient
291      * values for each parameter).
292      * Different from {@link #getParametersAtStateDate(FieldSpacecraftState, ParameterDriversProvider)}
293      * which return a Gradient list containing for each driver the gradient value at state date (only 1 gradient
294      * value for each parameter).
295      * @param state state as returned by {@link #getState(ParameterDriversProvider) getState(parametricModel)}
296      * @param parametricModel parametric model associated with the parameters
297      * @return parametric model parameters (for all span of each driver)
298      */
299     public Gradient[] getParameters(final FieldSpacecraftState<Gradient> state,
300                                     final ParameterDriversProvider parametricModel) {
301         final int freeParameters = state.getMass().getFreeParameters();
302         final List<ParameterDriver> drivers = parametricModel.getParametersDrivers();
303         int sizeDrivers = 0;
304         for ( ParameterDriver driver : drivers) {
305             sizeDrivers += driver.getNbOfValues();
306         }
307         final Gradient[] parameters = new Gradient[sizeDrivers];
308         int index = freeStateParameters;
309         int i = 0;
310         for (ParameterDriver driver : drivers) {
311             // Loop on the spans
312             for (Span<Double> span = driver.getValueSpanMap().getFirstSpan(); span != null; span = span.next()) {
313 
314                 parameters[i++] = driver.isSelected() ?
315                                   Gradient.variable(freeParameters, index++, span.getData()) :
316                                   Gradient.constant(freeParameters, span.getData());
317             }
318         }
319         return parameters;
320     }
321 
322     /** Get the parametric model parameters, return gradient values at state date for each driver (only 1 gradient
323      * value for each parameter).
324      * Different from {@link #getParameters(FieldSpacecraftState, ParameterDriversProvider)}
325      * which return a Gradient list containing for each driver the gradient values for each span value (several gradient
326      * values for each parameter).
327      * @param state state as returned by {@link #getState(ParameterDriversProvider) getState(parametricModel)}
328      * @param parametricModel parametric model associated with the parameters
329      * @return parametric model parameters (for all span of each driver)
330      */
331     public Gradient[] getParametersAtStateDate(final FieldSpacecraftState<Gradient> state,
332             final ParameterDriversProvider parametricModel) {
333         final int freeParameters = state.getMass().getFreeParameters();
334         final List<ParameterDriver> drivers = parametricModel.getParametersDrivers();
335 
336         final Gradient[] parameters = new Gradient[drivers.size()];
337         int index = freeStateParameters;
338         int i = 0;
339         for (ParameterDriver driver : drivers) {
340             for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
341                 if (span.getData().equals(driver.getNameSpan(state.getDate().toAbsoluteDate()))) {
342                     parameters[i++] = driver.isSelected() ?
343                                           Gradient.variable(freeParameters, index, driver.getValue(state.getDate().toAbsoluteDate())) :
344                                           Gradient.constant(freeParameters, driver.getValue(state.getDate().toAbsoluteDate()));
345                 }
346                 index = driver.isSelected() ? index + 1 : index;
347             }
348         }
349         return parameters;
350     }
351 
352 
353 }