1   /* Copyright 2002-2025 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;
18  
19  import java.util.List;
20  
21  import org.hipparchus.linear.MatrixUtils;
22  import org.hipparchus.linear.RealMatrix;
23  import org.orekit.orbits.PositionAngleType;
24  import org.orekit.utils.DoubleArrayDictionary;
25  
26  /** Base harvester between two-dimensional Jacobian matrices and one-dimensional {@link
27   * SpacecraftState#getAdditionalState(String) additional state arrays}.
28   * @author Luc Maisonobe
29   * @since 11.1
30   */
31  public abstract class AbstractMatricesHarvester implements MatricesHarvester {
32  
33      /** State dimension, fixed to 6.
34       * @deprecated as of 13.1, use DEFAULT_STATE_DIMENSION */
35      @Deprecated
36      public static final int STATE_DIMENSION = 6;
37  
38      /** Default state dimension, equivalent to position and velocity vectors. */
39      public static final int DEFAULT_STATE_DIMENSION = 6;
40  
41      /** Identity conversion matrix for Cartesian-like coordinates. */
42      private static final double[][] IDENTITY6 = {
43          { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
44          { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
45          { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
46          { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
47          { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
48          { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 }
49      };
50  
51      /** Initial State Transition Matrix. */
52      private final RealMatrix initialStm;
53  
54      /** Initial columns of the Jacobians matrix with respect to parameters. */
55      private final DoubleArrayDictionary initialJacobianColumns;
56  
57      /** State Transition Matrix state name. */
58      private final String stmName;
59  
60      /** Simple constructor.
61       * <p>
62       * The arguments for initial matrices <em>must</em> be compatible with the {@link org.orekit.orbits.OrbitType orbit type}
63       * and {@link PositionAngleType position angle} that will be used by propagator
64       * </p>
65       * @param stmName State Transition Matrix state name
66       * @param initialStm initial State Transition Matrix ∂Y/∂Y₀,
67       * if null (which is the most frequent case), assumed to be 6x6 identity
68       * @param initialJacobianColumns initial columns of the Jacobians matrix with respect to parameters,
69       * if null or if some selected parameters are missing from the dictionary, the corresponding
70       * initial column is assumed to be 0
71       */
72      protected AbstractMatricesHarvester(final String stmName, final RealMatrix initialStm, final DoubleArrayDictionary initialJacobianColumns) {
73          this.stmName                = stmName;
74          this.initialStm             = initialStm == null ? MatrixUtils.createRealIdentityMatrix(DEFAULT_STATE_DIMENSION) : initialStm;
75          this.initialJacobianColumns = initialJacobianColumns == null ? new DoubleArrayDictionary() : initialJacobianColumns;
76      }
77  
78      /**
79       * Getter for the state dimension.
80       * @return state dimension
81       * @since 13.1
82       */
83      public int getStateDimension() {
84          return initialStm.getColumnDimension();
85      }
86  
87      /** Get the State Transition Matrix state name.
88       * @return State Transition Matrix state name
89       */
90      public String getStmName() {
91          return stmName;
92      }
93  
94      /** Get the initial State Transition Matrix.
95       * @return initial State Transition Matrix
96       */
97      public RealMatrix getInitialStateTransitionMatrix() {
98          return initialStm;
99      }
100 
101     /** Get the initial column of Jacobian matrix with respect to named parameter.
102      * @param columnName name of the column
103      * @return initial column of the Jacobian matrix
104      */
105     public double[] getInitialJacobianColumn(final String columnName) {
106         final DoubleArrayDictionary.Entry entry = initialJacobianColumns.getEntry(columnName);
107         return entry == null ? new double[getStateDimension()] : entry.getValue();
108     }
109 
110     /** Convert a flattened array to a square matrix.
111      * @param array input array
112      * @return the corresponding matrix
113      * @since 13.1
114      */
115     public RealMatrix toSquareMatrix(final double[] array) {
116         final int stateDimension = getStateDimension();
117         final RealMatrix matrix = MatrixUtils.createRealMatrix(stateDimension, stateDimension);
118         int index = 0;
119         for (int i = 0; i < stateDimension; ++i) {
120             for (int j = 0; j < stateDimension; ++j) {
121                 matrix.setEntry(i, j, array[index++]);
122             }
123         }
124         return matrix;
125     }
126 
127     /** Set the STM data into an array.
128      * @param matrix STM matrix
129      * @return an array containing the STM data
130      * @since 13.1
131      */
132     public double[] toArray(final double[][] matrix) {
133         final int stateDimension = matrix.length;
134         final double[] array = new double[stateDimension * stateDimension];
135         int index = 0;
136         for (final double[] row : matrix) {
137             for (int j = 0; j < stateDimension; ++j) {
138                 array[index++] = row[j];
139             }
140         }
141         return array;
142     }
143 
144     /** Get the conversion Jacobian between state parameters and parameters used for derivatives.
145      * <p>
146      * The base implementation returns identity, which is suitable for DSST and TLE propagators,
147      * as state parameters and parameters used for derivatives are the same.
148      * </p>
149      * <p>
150      * For Numerical propagator, parameters used for derivatives are Cartesian
151      * and they can be different from state parameters because the numerical propagator can accept different type
152      * of orbits, so the method is overridden in derived classes.
153      * </p>
154      * @param state spacecraft state
155      * @return conversion Jacobian
156      */
157     protected double[][] getConversionJacobian(final SpacecraftState state) {
158         return IDENTITY6;
159     }
160 
161     /** {@inheritDoc} */
162     @Override
163     public void setReferenceState(final SpacecraftState reference) {
164         // nothing to do
165     }
166 
167     /** {@inheritDoc} */
168     @Override
169     public RealMatrix getStateTransitionMatrix(final SpacecraftState state) {
170 
171         if (!state.hasAdditionalData(stmName)) {
172             return null;
173         }
174 
175         // get the conversion Jacobian
176         final double[][] dYdC = getConversionJacobian(state);
177 
178         // extract the additional state
179         final double[] p = state.getAdditionalState(stmName);
180 
181         // compute dYdY0 = dYdC * dCdY0
182         final int stateDimension = getStateDimension();
183         final RealMatrix  dYdY0 = MatrixUtils.createRealMatrix(stateDimension, stateDimension);
184         for (int i = 0; i < stateDimension; i++) {
185             final double[] rowC = dYdC[i];
186             for (int j = 0; j < stateDimension; ++j) {
187                 double sum = 0;
188                 int pIndex = j;
189                 for (int k = 0; k < stateDimension; ++k) {
190                     sum += rowC[k] * p[pIndex];
191                     pIndex += stateDimension;
192                 }
193                 dYdY0.setEntry(i, j, sum);
194             }
195         }
196 
197         return dYdY0;
198 
199     }
200 
201     /** {@inheritDoc} */
202     @Override
203     public RealMatrix getParametersJacobian(final SpacecraftState state) {
204 
205         final List<String> columnsNames = getJacobiansColumnsNames();
206 
207         if (columnsNames == null || columnsNames.isEmpty()) {
208             return null;
209         }
210 
211         // get the conversion Jacobian
212         final RealMatrix dYdC = MatrixUtils.createRealIdentityMatrix(getStateDimension());
213         dYdC.setSubMatrix(getConversionJacobian(state), 0, 0);
214 
215         // compute dYdP = dYdC * dCdP
216         final RealMatrix dYdP = MatrixUtils.createRealMatrix(getStateDimension(), columnsNames.size());
217         for (int j = 0; j < columnsNames.size(); j++) {
218             final double[] p = state.getAdditionalState(columnsNames.get(j));
219             for (int i = 0; i < getStateDimension(); ++i) {
220                 final double[] dYdCi = dYdC.getRow(i);
221                 double sum = 0;
222                 for (int k = 0; k < getStateDimension(); ++k) {
223                     sum += dYdCi[k] * p[k];
224                 }
225                 dYdP.setEntry(i, j, sum);
226             }
227         }
228 
229         return dYdP;
230 
231     }
232 
233     /** Freeze the names of the Jacobian columns.
234      * <p>
235      * This method is called when propagation starts, i.e. when configuration is completed
236      * </p>
237      */
238     public abstract void freezeColumnsNames();
239 
240 }