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.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.propagation.AbstractMatricesHarvester;
27  import org.orekit.propagation.FieldSpacecraftState;
28  import org.orekit.propagation.SpacecraftState;
29  import org.orekit.time.AbsoluteDate;
30  import org.orekit.time.FieldAbsoluteDate;
31  import org.orekit.utils.DoubleArrayDictionary;
32  import org.orekit.utils.FieldPVCoordinates;
33  import org.orekit.utils.ParameterDriver;
34  
35  /**
36   * Base class harvester between two-dimensional Jacobian
37   * matrices and analytical orbit propagator.
38   * @author Thomas Paulet
39   * @author Bryan Cazabonne
40   * @since 11.1
41   */
42  public abstract class AbstractAnalyticalMatricesHarvester extends AbstractMatricesHarvester {
43  
44      /** Columns names for parameters. */
45      private List<String> columnsNames;
46  
47      /** Epoch of the last computed state transition matrix. */
48      private AbsoluteDate epoch;
49  
50      /** Analytical derivatives that apply to State Transition Matrix. */
51      private final double[][] analyticalDerivativesStm;
52  
53      /** Analytical derivatives that apply to Jacobians columns. */
54      private final DoubleArrayDictionary analyticalDerivativesJacobianColumns;
55  
56      /** Propagator bound to this harvester. */
57      private final AbstractAnalyticalPropagator propagator;
58  
59      /** Simple constructor.
60       * <p>
61       * The arguments for initial matrices <em>must</em> be compatible with the
62       * {@link org.orekit.orbits.OrbitType orbit type}
63       * and {@link org.orekit.orbits.PositionAngle position angle} that will be used by propagator
64       * </p>
65       * @param propagator propagator bound to this harvester
66       * @param stmName State Transition Matrix state name
67       * @param initialStm initial State Transition Matrix ∂Y/∂Y₀,
68       * if null (which is the most frequent case), assumed to be 6x6 identity
69       * @param initialJacobianColumns initial columns of the Jacobians matrix with respect to parameters,
70       * if null or if some selected parameters are missing from the dictionary, the corresponding
71       * initial column is assumed to be 0
72       */
73      protected AbstractAnalyticalMatricesHarvester(final AbstractAnalyticalPropagator propagator, final String stmName,
74                                                    final RealMatrix initialStm, final DoubleArrayDictionary initialJacobianColumns) {
75          super(stmName, initialStm, initialJacobianColumns);
76          this.propagator                           = propagator;
77          this.epoch                                = propagator.getInitialState().getDate();
78          this.columnsNames                         = null;
79          this.analyticalDerivativesStm             = getInitialStateTransitionMatrix().getData();
80          this.analyticalDerivativesJacobianColumns = new DoubleArrayDictionary();
81      }
82  
83      /** {@inheritDoc} */
84      @Override
85      public List<String> getJacobiansColumnsNames() {
86          return columnsNames == null ? propagator.getJacobiansColumnsNames() : columnsNames;
87      }
88  
89      /** {@inheritDoc} */
90      @Override
91      public void freezeColumnsNames() {
92          columnsNames = getJacobiansColumnsNames();
93      }
94  
95      /** {@inheritDoc} */
96      @Override
97      public RealMatrix getStateTransitionMatrix(final SpacecraftState state) {
98          // Update the partial derivatives if needed
99          updateDerivativesIfNeeded(state);
100         //return the state transition matrix
101         return MatrixUtils.createRealMatrix(analyticalDerivativesStm);
102     }
103 
104     /** {@inheritDoc} */
105     @Override
106     public RealMatrix getParametersJacobian(final SpacecraftState state) {
107         // Update the partial derivatives if needed
108         updateDerivativesIfNeeded(state);
109 
110         // Estimated parameters
111         final List<String> names = getJacobiansColumnsNames();
112         if (names == null || names.isEmpty()) {
113             return null;
114         }
115 
116         // Initialize Jacobian
117         final RealMatrix dYdP = MatrixUtils.createRealMatrix(STATE_DIMENSION, names.size());
118 
119         // Add the derivatives
120         for (int j = 0; j < names.size(); ++j) {
121             final double[] column = analyticalDerivativesJacobianColumns.get(names.get(j));
122             if (column != null) {
123                 for (int i = 0; i < STATE_DIMENSION; i++) {
124                     dYdP.addToEntry(i, j, column[i]);
125                 }
126             }
127         }
128 
129         // Return
130         return dYdP;
131     }
132 
133     /** {@inheritDoc} */
134     @Override
135     public void setReferenceState(final SpacecraftState reference) {
136 
137         // reset derivatives to zero
138         for (final double[] row : analyticalDerivativesStm) {
139             Arrays.fill(row, 0.0);
140         }
141         analyticalDerivativesJacobianColumns.clear();
142 
143         final AbstractAnalyticalGradientConverter converter           = getGradientConverter();
144         final FieldSpacecraftState<Gradient> gState                   = converter.getState();
145         final Gradient[] gParameters                                  = converter.getParameters(gState);
146         final FieldAbstractAnalyticalPropagator<Gradient> gPropagator = converter.getPropagator(gState, gParameters);
147 
148         // Compute Jacobian
149         final AbsoluteDate target               = reference.getDate();
150         final FieldAbsoluteDate<Gradient> start = gPropagator.getInitialState().getDate();
151         final double dt                         = target.durationFrom(start.toAbsoluteDate());
152         final FieldOrbit<Gradient> gOrbit       = gPropagator.propagateOrbit(start.shiftedBy(dt), gParameters);
153         final FieldPVCoordinates<Gradient> gPv  = gOrbit.getPVCoordinates();
154 
155         final double[] derivativesX   = gPv.getPosition().getX().getGradient();
156         final double[] derivativesY   = gPv.getPosition().getY().getGradient();
157         final double[] derivativesZ   = gPv.getPosition().getZ().getGradient();
158         final double[] derivativesVx  = gPv.getVelocity().getX().getGradient();
159         final double[] derivativesVy  = gPv.getVelocity().getY().getGradient();
160         final double[] derivativesVz  = gPv.getVelocity().getZ().getGradient();
161 
162         // Update Jacobian with respect to state
163         addToRow(derivativesX,  0);
164         addToRow(derivativesY,  1);
165         addToRow(derivativesZ,  2);
166         addToRow(derivativesVx, 3);
167         addToRow(derivativesVy, 4);
168         addToRow(derivativesVz, 5);
169 
170         // Partial derivatives of the state with respect to propagation parameters
171         int paramsIndex = converter.getFreeStateParameters();
172         for (ParameterDriver driver : converter.getParametersDrivers()) {
173             if (driver.isSelected()) {
174 
175                 // get the partials derivatives for this driver
176                 DoubleArrayDictionary.Entry entry = analyticalDerivativesJacobianColumns.getEntry(driver.getName());
177                 if (entry == null) {
178                     // create an entry filled with zeroes
179                     analyticalDerivativesJacobianColumns.put(driver.getName(), new double[STATE_DIMENSION]);
180                     entry = analyticalDerivativesJacobianColumns.getEntry(driver.getName());
181                 }
182 
183                 // add the contribution of the current force model
184                 entry.increment(new double[] {
185                     derivativesX[paramsIndex], derivativesY[paramsIndex], derivativesZ[paramsIndex],
186                     derivativesVx[paramsIndex], derivativesVy[paramsIndex], derivativesVz[paramsIndex]
187                 });
188                 ++paramsIndex;
189 
190             }
191         }
192 
193         // Update the epoch of the last computed partial derivatives
194         epoch = target;
195 
196     }
197 
198     /** Update the partial derivatives (if needed).
199      * @param state current spacecraft state
200      */
201     private void updateDerivativesIfNeeded(final SpacecraftState state) {
202         if (!state.getDate().isEqualTo(epoch)) {
203             setReferenceState(state);
204         }
205     }
206 
207     /** Fill State Transition Matrix rows.
208      * @param derivatives derivatives of a component
209      * @param index component index
210      */
211     private void addToRow(final double[] derivatives, final int index) {
212         for (int i = 0; i < 6; i++) {
213             analyticalDerivativesStm[index][i] += derivatives[i];
214         }
215     }
216 
217     /**
218      * Get the gradient converter related to the analytical orbit propagator.
219      * @return the gradient converter
220      */
221     public abstract AbstractAnalyticalGradientConverter getGradientConverter();
222 
223 }