1   /* Copyright 2002-2022 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.utils.DoubleArrayDictionary;
24  
25  /** Base harvester between two-dimensional Jacobian matrices and one-dimensional {@link
26   * SpacecraftState#getAdditionalState(String) additional state arrays}.
27   * @author Luc Maisonobe
28   * @since 11.1
29   */
30  public abstract class AbstractMatricesHarvester implements MatricesHarvester {
31  
32      /** State dimension, fixed to 6. */
33      public static final int STATE_DIMENSION = 6;
34  
35      /** Identity conversion matrix. */
36      private static final double[][] IDENTITY = {
37          { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 },
38          { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
39          { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 },
40          { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 },
41          { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 },
42          { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 }
43      };
44  
45      /** Initial State Transition Matrix. */
46      private final RealMatrix initialStm;
47  
48      /** Initial columns of the Jacobians matrix with respect to parameters. */
49      private final DoubleArrayDictionary initialJacobianColumns;
50  
51      /** State Transition Matrix state name. */
52      private final String stmName;
53  
54      /** Simple constructor.
55       * <p>
56       * The arguments for initial matrices <em>must</em> be compatible with the {@link org.orekit.orbits.OrbitType orbit type}
57       * and {@link org.orekit.orbits.PositionAngle position angle} that will be used by propagator
58       * </p>
59       * @param stmName State Transition Matrix state name
60       * @param initialStm initial State Transition Matrix ∂Y/∂Y₀,
61       * if null (which is the most frequent case), assumed to be 6x6 identity
62       * @param initialJacobianColumns initial columns of the Jacobians matrix with respect to parameters,
63       * if null or if some selected parameters are missing from the dictionary, the corresponding
64       * initial column is assumed to be 0
65       */
66      protected AbstractMatricesHarvester(final String stmName, final RealMatrix initialStm, final DoubleArrayDictionary initialJacobianColumns) {
67          this.stmName                = stmName;
68          this.initialStm             = initialStm == null ? MatrixUtils.createRealIdentityMatrix(STATE_DIMENSION) : initialStm;
69          this.initialJacobianColumns = initialJacobianColumns == null ? new DoubleArrayDictionary() : initialJacobianColumns;
70      }
71  
72      /** Get the State Transition Matrix state name.
73       * @return State Transition Matrix state name
74       */
75      public String getStmName() {
76          return stmName;
77      }
78  
79      /** Get the initial State Transition Matrix.
80       * @return initial State Transition Matrix
81       */
82      public RealMatrix getInitialStateTransitionMatrix() {
83          return initialStm;
84      }
85  
86      /** Get the initial column of Jacobian matrix with respect to named parameter.
87       * @param columnName name of the column
88       * @return initial column of the Jacobian matrix
89       */
90      public double[] getInitialJacobianColumn(final String columnName) {
91          final DoubleArrayDictionary.Entry entry = initialJacobianColumns.getEntry(columnName);
92          return entry == null ? new double[STATE_DIMENSION] : entry.getValue();
93      }
94  
95      /** Get the conversion Jacobian between state parameters and parameters used for derivatives.
96       * <p>
97       * The base implementation returns identity, which is suitable for DSST and TLE propagators,
98       * as state parameters and parameters used for derivatives are the same.
99       * </p>
100      * <p>
101      * For Numerical propagator, parameters used for derivatives are Cartesian
102      * and they can be different from state parameters because the numerical propagator can accept different type
103      * of orbits, so the method is overridden in derived classes.
104      * </p>
105      * @param state spacecraft state
106      * @return conversion Jacobian
107      */
108     protected double[][] getConversionJacobian(final SpacecraftState state) {
109         return IDENTITY;
110     }
111 
112     /** {@inheritDoc} */
113     @Override
114     public void setReferenceState(final SpacecraftState reference) {
115         // nothing to do
116     }
117 
118     /** {@inheritDoc} */
119     @Override
120     public RealMatrix getStateTransitionMatrix(final SpacecraftState state) {
121 
122         if (!state.hasAdditionalState(stmName)) {
123             return null;
124         }
125 
126         // get the conversion Jacobian
127         final double[][] dYdC = getConversionJacobian(state);
128 
129         // extract the additional state
130         final double[] p = state.getAdditionalState(stmName);
131 
132         // compute dYdY0 = dYdC * dCdY0
133         final RealMatrix  dYdY0 = MatrixUtils.createRealMatrix(STATE_DIMENSION, STATE_DIMENSION);
134         for (int i = 0; i < STATE_DIMENSION; i++) {
135             final double[] rowC = dYdC[i];
136             for (int j = 0; j < STATE_DIMENSION; ++j) {
137                 double sum = 0;
138                 int pIndex = j;
139                 for (int k = 0; k < STATE_DIMENSION; ++k) {
140                     sum += rowC[k] * p[pIndex];
141                     pIndex += STATE_DIMENSION;
142                 }
143                 dYdY0.setEntry(i, j, sum);
144             }
145         }
146 
147         return dYdY0;
148 
149     }
150 
151     /** {@inheritDoc} */
152     @Override
153     public RealMatrix getParametersJacobian(final SpacecraftState state) {
154 
155         final List<String> columnsNames = getJacobiansColumnsNames();
156 
157         if (columnsNames == null || columnsNames.isEmpty()) {
158             return null;
159         }
160 
161         // get the conversion Jacobian
162         final double[][] dYdC = getConversionJacobian(state);
163 
164         // compute dYdP = dYdC * dCdP
165         final RealMatrix dYdP = MatrixUtils.createRealMatrix(STATE_DIMENSION, columnsNames.size());
166         for (int j = 0; j < columnsNames.size(); j++) {
167             final double[] p = state.getAdditionalState(columnsNames.get(j));
168             for (int i = 0; i < STATE_DIMENSION; ++i) {
169                 final double[] dYdCi = dYdC[i];
170                 double sum = 0;
171                 for (int k = 0; k < STATE_DIMENSION; ++k) {
172                     sum += dYdCi[k] * p[k];
173                 }
174                 dYdP.setEntry(i, j, sum);
175             }
176         }
177 
178         return dYdP;
179 
180     }
181 
182     /** Freeze the names of the Jacobian columns.
183      * <p>
184      * This method is called when propagation starts, i.e. when configuration is completed
185      * </p>
186      */
187     public abstract void freezeColumnsNames();
188 
189 }