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