1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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.orbits.OrbitType;
27 import org.orekit.orbits.PositionAngleType;
28 import org.orekit.propagation.AbstractMatricesHarvester;
29 import org.orekit.propagation.AdditionalDataProvider;
30 import org.orekit.propagation.FieldSpacecraftState;
31 import org.orekit.propagation.SpacecraftState;
32 import org.orekit.time.AbsoluteDate;
33 import org.orekit.time.FieldAbsoluteDate;
34 import org.orekit.utils.DoubleArrayDictionary;
35 import org.orekit.utils.ParameterDriver;
36 import org.orekit.utils.TimeSpanMap;
37 import org.orekit.utils.TimeSpanMap.Span;
38
39
40
41
42
43
44
45
46 public abstract class AbstractAnalyticalMatricesHarvester extends AbstractMatricesHarvester
47 implements AdditionalDataProvider<double[]> {
48
49
50 private List<String> columnsNames;
51
52
53 private AbsoluteDate epoch;
54
55
56 private double[][] analyticalDerivativesStm;
57
58
59 private final DoubleArrayDictionary analyticalDerivativesJacobianColumns;
60
61
62 private final AbstractAnalyticalPropagator propagator;
63
64
65
66
67
68
69
70
71
72 protected AbstractAnalyticalMatricesHarvester(final AbstractAnalyticalPropagator propagator) {
73 this.propagator = propagator;
74 this.epoch = propagator.getInitialState().getDate();
75 this.columnsNames = null;
76 this.analyticalDerivativesJacobianColumns = new DoubleArrayDictionary();
77 }
78
79
80 @Override
81 protected void setInitialStm(final String stmName, final RealMatrix stm) {
82 super.setInitialStm(stmName, stm);
83 this.analyticalDerivativesStm = getInitialStateTransitionMatrix().getData();
84 }
85
86
87 @Override
88 public List<String> getJacobiansColumnsNames() {
89 return columnsNames == null ? propagator.getJacobiansColumnsNames() : columnsNames;
90 }
91
92
93 @Override
94 public void freezeColumnsNames() {
95 columnsNames = getJacobiansColumnsNames();
96 }
97
98
99 @Override
100 public String getName() {
101 return getStmName();
102 }
103
104
105 @Override
106 public double[] getAdditionalData(final SpacecraftState state) {
107
108 updateDerivativesIfNeeded(state);
109
110 return toArray(analyticalDerivativesStm);
111 }
112
113
114 @Override
115 public RealMatrix getStateTransitionMatrix(final SpacecraftState state) {
116
117 if (!state.hasAdditionalData(getName())) {
118 return null;
119 }
120
121 return toSquareMatrix(state.getAdditionalState(getName()));
122 }
123
124
125 @Override
126 public RealMatrix getParametersJacobian(final SpacecraftState state) {
127
128
129 updateDerivativesIfNeeded(state);
130
131
132 final List<String> names = getJacobiansColumnsNames();
133 if (names == null || names.isEmpty()) {
134 return null;
135 }
136
137
138 final RealMatrix dYdP = MatrixUtils.createRealMatrix(getStateDimension(), names.size());
139
140
141 for (int j = 0; j < names.size(); ++j) {
142 final double[] column = analyticalDerivativesJacobianColumns.get(names.get(j));
143 if (column != null) {
144 for (int i = 0; i < getStateDimension(); i++) {
145 dYdP.addToEntry(i, j, column[i]);
146 }
147 }
148 }
149
150
151 return dYdP;
152
153 }
154
155
156 @Override
157 public void setReferenceState(final SpacecraftState reference) {
158
159
160 for (final double[] row : analyticalDerivativesStm) {
161 Arrays.fill(row, 0.0);
162 }
163 analyticalDerivativesJacobianColumns.clear();
164
165 final AbstractAnalyticalGradientConverter converter = getGradientConverter();
166 final FieldAbstractAnalyticalPropagator<Gradient> gPropagator = converter.getPropagator();
167
168
169 final AbsoluteDate target = reference.getDate();
170 final FieldAbsoluteDate<Gradient> start = gPropagator.getInitialState().getDate();
171 final double dt = target.durationFrom(start.toAbsoluteDate());
172 final FieldSpacecraftState<Gradient> state = gPropagator.getInitialState();
173 final Gradient[] parameters = converter.getParameters(state, converter);
174 final FieldOrbit<Gradient> gOrbit = gPropagator.propagateOrbit(start.shiftedBy(dt), parameters);
175 final Gradient[] orbitDerivatives = new Gradient[6];
176 getOrbitType().mapOrbitToArray(gOrbit, getPositionAngleType(), orbitDerivatives, null);
177
178
179 for (int i = 0; i < orbitDerivatives.length; ++i) {
180 addToRow(orbitDerivatives[i].getGradient(), i);
181 }
182
183
184 int paramsIndex = converter.getFreeStateParameters();
185 for (ParameterDriver driver : converter.getParametersDrivers()) {
186 if (driver.isSelected()) {
187
188 final TimeSpanMap<String> driverNameSpanMap = driver.getNamesSpanMap();
189
190 for (Span<String> span = driverNameSpanMap.getFirstSpan(); span != null; span = span.next()) {
191
192 DoubleArrayDictionary.Entry entry = analyticalDerivativesJacobianColumns.getEntry(span.getData());
193 if (entry == null) {
194
195 analyticalDerivativesJacobianColumns.put(span.getData(), new double[getStateDimension()]);
196 entry = analyticalDerivativesJacobianColumns.getEntry(span.getData());
197 }
198
199
200 entry.increment(new double[] {
201 orbitDerivatives[0].getPartialDerivative(paramsIndex),
202 orbitDerivatives[1].getPartialDerivative(paramsIndex),
203 orbitDerivatives[2].getPartialDerivative(paramsIndex),
204 orbitDerivatives[3].getPartialDerivative(paramsIndex),
205 orbitDerivatives[4].getPartialDerivative(paramsIndex),
206 orbitDerivatives[5].getPartialDerivative(paramsIndex)
207 });
208 ++paramsIndex;
209 }
210 }
211 }
212
213
214 epoch = target;
215
216 }
217
218
219
220
221 private void updateDerivativesIfNeeded(final SpacecraftState state) {
222 if (!state.getDate().isEqualTo(epoch)) {
223 setReferenceState(state);
224 }
225 }
226
227
228
229
230
231 private void addToRow(final double[] derivatives, final int index) {
232 for (int i = 0; i < 6; i++) {
233 analyticalDerivativesStm[index][i] += derivatives[i];
234 }
235 }
236
237
238 @Override
239 public OrbitType getOrbitType() {
240
241 return OrbitType.CARTESIAN;
242 }
243
244
245 @Override
246 public PositionAngleType getPositionAngleType() {
247
248 return PositionAngleType.MEAN;
249 }
250
251
252
253
254
255 public abstract AbstractAnalyticalGradientConverter getGradientConverter();
256
257 }