1   /* Copyright 2002-2026 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.analytical;
18  
19  import java.util.Arrays;
20  import java.util.List;
21  
22  import org.hipparchus.analysis.differentiation.Gradient;
23  import org.hipparchus.linear.MatrixUtils;
24  import org.hipparchus.linear.RealMatrix;
25  import org.orekit.orbits.FieldOrbit;
26  import org.orekit.orbits.OrbitType;
27  import org.orekit.orbits.PositionAngleType;
28  import org.orekit.propagation.AbstractMatricesHarvester;
29  import org.orekit.propagation.AdditionalDataProvider;
30  import org.orekit.propagation.FieldSpacecraftState;
31  import org.orekit.propagation.SpacecraftState;
32  import org.orekit.time.AbsoluteDate;
33  import org.orekit.time.FieldAbsoluteDate;
34  import org.orekit.utils.DoubleArrayDictionary;
35  import org.orekit.utils.ParameterDriver;
36  import org.orekit.utils.TimeSpanMap;
37  import org.orekit.utils.TimeSpanMap.Span;
38  
39  /**
40   * Base class harvester between two-dimensional Jacobian
41   * matrices and analytical orbit propagator.
42   * @author Thomas Paulet
43   * @author Bryan Cazabonne
44   * @since 11.1
45   */
46  public abstract class AbstractAnalyticalMatricesHarvester extends AbstractMatricesHarvester
47          implements AdditionalDataProvider<double[]> {
48  
49      /** Columns names for parameters. */
50      private List<String> columnsNames;
51  
52      /** Epoch of the last computed state transition matrix. */
53      private AbsoluteDate epoch;
54  
55      /** Analytical derivatives that apply to State Transition Matrix. */
56      private double[][] analyticalDerivativesStm;
57  
58      /** Analytical derivatives that apply to Jacobians columns. */
59      private final DoubleArrayDictionary analyticalDerivativesJacobianColumns;
60  
61      /** Propagator bound to this harvester. */
62      private final AbstractAnalyticalPropagator propagator;
63  
64      /** Simple constructor.
65       * <p>
66       * The arguments for initial matrices <em>must</em> be compatible with the
67       * {@link org.orekit.orbits.OrbitType orbit type}
68       * and {@link PositionAngleType position angle} that will be used by propagator
69       * </p>
70       * @param propagator propagator bound to this harvester
71       */
72      protected AbstractAnalyticalMatricesHarvester(final AbstractAnalyticalPropagator propagator) {
73          this.propagator                           = propagator;
74          this.epoch                                = propagator.getInitialState().getDate();
75          this.columnsNames                         = null;
76          this.analyticalDerivativesJacobianColumns = new DoubleArrayDictionary();
77      }
78  
79      /** {@inheritDoc} */
80      @Override
81      protected void setInitialStm(final String stmName, final RealMatrix stm) {
82          super.setInitialStm(stmName, stm);
83          this.analyticalDerivativesStm = getInitialStateTransitionMatrix().getData();
84      }
85  
86      /** {@inheritDoc} */
87      @Override
88      public List<String> getJacobiansColumnsNames() {
89          return columnsNames == null ? propagator.getJacobiansColumnsNames() : columnsNames;
90      }
91  
92      /** {@inheritDoc} */
93      @Override
94      public void freezeColumnsNames() {
95          columnsNames = getJacobiansColumnsNames();
96      }
97  
98      /** {@inheritDoc} */
99      @Override
100     public String getName() {
101         return getStmName();
102     }
103 
104     /** {@inheritDoc} */
105     @Override
106     public double[] getAdditionalData(final SpacecraftState state) {
107         // Update the partial derivatives if needed
108         updateDerivativesIfNeeded(state);
109         // Return the state transition matrix in an array
110         return toArray(analyticalDerivativesStm);
111     }
112 
113     /** {@inheritDoc} */
114     @Override
115     public RealMatrix getStateTransitionMatrix(final SpacecraftState state) {
116         // Check if additional state is defined
117         if (!state.hasAdditionalData(getName())) {
118             return null;
119         }
120         // Return the state transition matrix
121         return toSquareMatrix(state.getAdditionalState(getName()));
122     }
123 
124     /** {@inheritDoc} */
125     @Override
126     public RealMatrix getParametersJacobian(final SpacecraftState state) {
127 
128         // Update the partial derivatives if needed
129         updateDerivativesIfNeeded(state);
130 
131         // Estimated parameters
132         final List<String> names = getJacobiansColumnsNames();
133         if (names == null || names.isEmpty()) {
134             return null;
135         }
136 
137         // Initialize Jacobian
138         final RealMatrix dYdP = MatrixUtils.createRealMatrix(getStateDimension(), names.size());
139 
140         // Add the derivatives
141         for (int j = 0; j < names.size(); ++j) {
142             final double[] column = analyticalDerivativesJacobianColumns.get(names.get(j));
143             if (column != null) {
144                 for (int i = 0; i < getStateDimension(); i++) {
145                     dYdP.addToEntry(i, j, column[i]);
146                 }
147             }
148         }
149 
150         // Return
151         return dYdP;
152 
153     }
154 
155     /** {@inheritDoc} */
156     @Override
157     public void setReferenceState(final SpacecraftState reference) {
158 
159         // reset derivatives to zero
160         for (final double[] row : analyticalDerivativesStm) {
161             Arrays.fill(row, 0.0);
162         }
163         analyticalDerivativesJacobianColumns.clear();
164 
165         final AbstractAnalyticalGradientConverter converter           = getGradientConverter();
166         final FieldAbstractAnalyticalPropagator<Gradient> gPropagator = converter.getPropagator();
167 
168         // Compute Jacobian
169         final AbsoluteDate                   target           = reference.getDate();
170         final FieldAbsoluteDate<Gradient>    start            = gPropagator.getInitialState().getDate();
171         final double                         dt               = target.durationFrom(start.toAbsoluteDate());
172         final FieldSpacecraftState<Gradient> state            = gPropagator.getInitialState();
173         final Gradient[]                     parameters       = converter.getParameters(state, converter);
174         final FieldOrbit<Gradient>           gOrbit           = gPropagator.propagateOrbit(start.shiftedBy(dt), parameters);
175         final Gradient[]                     orbitDerivatives = new Gradient[6];
176         getOrbitType().mapOrbitToArray(gOrbit, getPositionAngleType(), orbitDerivatives, null);
177 
178         // Update Jacobian with respect to state
179         for (int i = 0; i < orbitDerivatives.length; ++i) {
180             addToRow(orbitDerivatives[i].getGradient(), i);
181         }
182 
183         // Partial derivatives of the state with respect to propagation parameters
184         int paramsIndex = converter.getFreeStateParameters();
185         for (ParameterDriver driver : converter.getParametersDrivers()) {
186             if (driver.isSelected()) {
187 
188                 final TimeSpanMap<String> driverNameSpanMap = driver.getNamesSpanMap();
189                 // for each span (for each estimated value) corresponding name is added
190                 for (Span<String> span = driverNameSpanMap.getFirstSpan(); span != null; span = span.next()) {
191                     // get the partials derivatives for this driver
192                     DoubleArrayDictionary.Entry entry = analyticalDerivativesJacobianColumns.getEntry(span.getData());
193                     if (entry == null) {
194                         // create an entry filled with zeroes
195                         analyticalDerivativesJacobianColumns.put(span.getData(), new double[getStateDimension()]);
196                         entry = analyticalDerivativesJacobianColumns.getEntry(span.getData());
197                     }
198 
199                     // add the contribution of the current force model
200                     entry.increment(new double[] {
201                         orbitDerivatives[0].getPartialDerivative(paramsIndex),
202                         orbitDerivatives[1].getPartialDerivative(paramsIndex),
203                         orbitDerivatives[2].getPartialDerivative(paramsIndex),
204                         orbitDerivatives[3].getPartialDerivative(paramsIndex),
205                         orbitDerivatives[4].getPartialDerivative(paramsIndex),
206                         orbitDerivatives[5].getPartialDerivative(paramsIndex)
207                     });
208                     ++paramsIndex;
209                 }
210             }
211         }
212 
213         // Update the epoch of the last computed partial derivatives
214         epoch = target;
215 
216     }
217 
218     /** Update the partial derivatives (if needed).
219      * @param state current spacecraft state
220      */
221     private void updateDerivativesIfNeeded(final SpacecraftState state) {
222         if (!state.getDate().isEqualTo(epoch)) {
223             setReferenceState(state);
224         }
225     }
226 
227     /** Fill State Transition Matrix rows.
228      * @param derivatives derivatives of a component
229      * @param index component index
230      */
231     private void addToRow(final double[] derivatives, final int index) {
232         for (int i = 0; i < 6; i++) {
233             analyticalDerivativesStm[index][i] += derivatives[i];
234         }
235     }
236 
237     /** {@inheritDoc} */
238     @Override
239     public OrbitType getOrbitType() {
240         // Set to CARTESIAN because analytical gradient converters uses Cartesian representation
241         return OrbitType.CARTESIAN;
242     }
243 
244     /** {@inheritDoc} */
245     @Override
246     public PositionAngleType getPositionAngleType() {
247         // Irrelevant: set a default value
248         return PositionAngleType.MEAN;
249     }
250 
251     /**
252      * Get the gradient converter related to the analytical orbit propagator.
253      * @return the gradient converter
254      */
255     public abstract AbstractAnalyticalGradientConverter getGradientConverter();
256 
257 }