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