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      public static final int STATE_DIMENSION = 6;
35  
36      /** Identity conversion matrix. */
37      private static final double[][] IDENTITY = {
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(STATE_DIMENSION) : initialStm;
70          this.initialJacobianColumns = initialJacobianColumns == null ? new DoubleArrayDictionary() : initialJacobianColumns;
71      }
72  
73      /** Get the State Transition Matrix state name.
74       * @return State Transition Matrix state name
75       */
76      public String getStmName() {
77          return stmName;
78      }
79  
80      /** Get the initial State Transition Matrix.
81       * @return initial State Transition Matrix
82       */
83      public RealMatrix getInitialStateTransitionMatrix() {
84          return initialStm;
85      }
86  
87      /** Get the initial column of Jacobian matrix with respect to named parameter.
88       * @param columnName name of the column
89       * @return initial column of the Jacobian matrix
90       */
91      public double[] getInitialJacobianColumn(final String columnName) {
92          final DoubleArrayDictionary.Entry entry = initialJacobianColumns.getEntry(columnName);
93          return entry == null ? new double[STATE_DIMENSION] : entry.getValue();
94      }
95  
96      /** Get the conversion Jacobian between state parameters and parameters used for derivatives.
97       * <p>
98       * The base implementation returns identity, which is suitable for DSST and TLE propagators,
99       * as state parameters and parameters used for derivatives are the same.
100      * </p>
101      * <p>
102      * For Numerical propagator, parameters used for derivatives are Cartesian
103      * and they can be different from state parameters because the numerical propagator can accept different type
104      * of orbits, so the method is overridden in derived classes.
105      * </p>
106      * @param state spacecraft state
107      * @return conversion Jacobian
108      */
109     protected double[][] getConversionJacobian(final SpacecraftState state) {
110         return IDENTITY;
111     }
112 
113     /** {@inheritDoc} */
114     @Override
115     public void setReferenceState(final SpacecraftState reference) {
116         // nothing to do
117     }
118 
119     /** {@inheritDoc} */
120     @Override
121     public RealMatrix getStateTransitionMatrix(final SpacecraftState state) {
122 
123         if (!state.hasAdditionalData(stmName)) {
124             return null;
125         }
126 
127         // get the conversion Jacobian
128         final double[][] dYdC = getConversionJacobian(state);
129 
130         // extract the additional state
131         final double[] p = state.getAdditionalState(stmName);
132 
133         // compute dYdY0 = dYdC * dCdY0
134         final RealMatrix  dYdY0 = MatrixUtils.createRealMatrix(STATE_DIMENSION, STATE_DIMENSION);
135         for (int i = 0; i < STATE_DIMENSION; i++) {
136             final double[] rowC = dYdC[i];
137             for (int j = 0; j < STATE_DIMENSION; ++j) {
138                 double sum = 0;
139                 int pIndex = j;
140                 for (int k = 0; k < STATE_DIMENSION; ++k) {
141                     sum += rowC[k] * p[pIndex];
142                     pIndex += STATE_DIMENSION;
143                 }
144                 dYdY0.setEntry(i, j, sum);
145             }
146         }
147 
148         return dYdY0;
149 
150     }
151 
152     /** {@inheritDoc} */
153     @Override
154     public RealMatrix getParametersJacobian(final SpacecraftState state) {
155 
156         final List<String> columnsNames = getJacobiansColumnsNames();
157 
158         if (columnsNames == null || columnsNames.isEmpty()) {
159             return null;
160         }
161 
162         // get the conversion Jacobian
163         final double[][] dYdC = getConversionJacobian(state);
164 
165         // compute dYdP = dYdC * dCdP
166         final RealMatrix dYdP = MatrixUtils.createRealMatrix(STATE_DIMENSION, columnsNames.size());
167         for (int j = 0; j < columnsNames.size(); j++) {
168             final double[] p = state.getAdditionalState(columnsNames.get(j));
169             for (int i = 0; i < STATE_DIMENSION; ++i) {
170                 final double[] dYdCi = dYdC[i];
171                 double sum = 0;
172                 for (int k = 0; k < STATE_DIMENSION; ++k) {
173                     sum += dYdCi[k] * p[k];
174                 }
175                 dYdP.setEntry(i, j, sum);
176             }
177         }
178 
179         return dYdP;
180 
181     }
182 
183     /** Freeze the names of the Jacobian columns.
184      * <p>
185      * This method is called when propagation starts, i.e. when configuration is completed
186      * </p>
187      */
188     public abstract void freezeColumnsNames();
189 
190 }