1 /* Copyright 2002-2018 CS Systèmes d'Information
2 * Licensed to CS Systèmes d'Information (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.estimation.measurements;
18
19 import java.util.Arrays;
20
21 import org.hipparchus.exception.LocalizedCoreFormats;
22 import org.hipparchus.geometry.euclidean.threed.Vector3D;
23 import org.hipparchus.util.FastMath;
24 import org.orekit.errors.OrekitException;
25 import org.orekit.propagation.SpacecraftState;
26 import org.orekit.time.AbsoluteDate;
27 import org.orekit.utils.TimeStampedPVCoordinates;
28
29 /** Class modeling a position-velocity state.
30 * @author Luc Maisonobe
31 * @since 8.0
32 */
33 public class PV extends AbstractMeasurement<PV> {
34
35 /** Identity matrix, for states derivatives. */
36 private static final double[][] IDENTITY = new double[][] {
37 {
38 1, 0, 0, 0, 0, 0
39 }, {
40 0, 1, 0, 0, 0, 0
41 }, {
42 0, 0, 1, 0, 0, 0
43 }, {
44 0, 0, 0, 1, 0, 0
45 }, {
46 0, 0, 0, 0, 1, 0
47 }, {
48 0, 0, 0, 0, 0, 1
49 }
50 };
51
52 /** Covariance matrix of the PV measurement (size 6x6). */
53 private final double[][] covarianceMatrix;
54
55 /** Constructor with two double for the standard deviations.
56 * The first double is the position's standard deviation, common to the 3 position's components.
57 * The second double is the position's standard deviation, common to the 3 position's components.
58 * <p>
59 * The measurement must be in the orbit propagation frame.
60 * </p>
61 * This constructor uses 0 as the index of the propagator related
62 * to this measurement, thus being well suited for mono-satellite
63 * orbit determination.
64 * @param date date of the measurement
65 * @param position position
66 * @param velocity velocity
67 * @param sigmaPosition theoretical standard deviation on position components
68 * @param sigmaVelocity theoretical standard deviation on velocity components
69 * @param baseWeight base weight
70 * @throws OrekitException if the built inside covariance matrix does not have the proper size
71 */
72 public PV(final AbsoluteDate date, final Vector3D position, final Vector3D velocity,
73 final double sigmaPosition, final double sigmaVelocity, final double baseWeight) throws OrekitException {
74 this(date, position, velocity, sigmaPosition, sigmaVelocity, baseWeight, 0);
75 }
76
77 /** Constructor with two double for the standard deviations.
78 * The first double is the position's standard deviation, common to the 3 position's components.
79 * The second double is the position's standard deviation, common to the 3 position's components.
80 * <p>
81 * The measurement must be in the orbit propagation frame.
82 * </p>
83 * @param date date of the measurement
84 * @param position position
85 * @param velocity velocity
86 * @param sigmaPosition theoretical standard deviation on position components
87 * @param sigmaVelocity theoretical standard deviation on velocity components
88 * @param baseWeight base weight
89 * @param propagatorIndex index of the propagator related to this measurement
90 * @throws OrekitException if the built inside covariance matrix does not have the proper size
91 * @since 9.0
92 */
93 public PV(final AbsoluteDate date, final Vector3D position, final Vector3D velocity,
94 final double sigmaPosition, final double sigmaVelocity, final double baseWeight,
95 final int propagatorIndex) throws OrekitException {
96 this(date, position, velocity,
97 new double[] {
98 sigmaPosition,
99 sigmaPosition,
100 sigmaPosition,
101 sigmaVelocity,
102 sigmaVelocity,
103 sigmaVelocity
104 }, baseWeight, propagatorIndex);
105 }
106
107 /** Constructor with two vectors for the standard deviations and default value for propagator index..
108 * One 3-sized vectors for position standard deviations.
109 * One 3-sized vectors for velocity standard deviations.
110 * The 3-sized vectors are the square root of the diagonal elements of the covariance matrix.
111 * <p>The measurement must be in the orbit propagation frame.</p>
112 * This constructor uses 0 as the index of the propagator related
113 * to this measurement, thus being well suited for mono-satellite
114 * orbit determination.
115 * @param date date of the measurement
116 * @param position position
117 * @param velocity velocity
118 * @param sigmaPosition 3-sized vector of the standard deviations of the position
119 * @param sigmaVelocity 3-sized vector of the standard deviations of the velocity
120 * @param baseWeight base weight
121 * @throws OrekitException if input standard deviations vectors do not have the proper sizes
122 * @since 9.2
123 */
124 public PV(final AbsoluteDate date, final Vector3D position, final Vector3D velocity,
125 final double[] sigmaPosition, final double[] sigmaVelocity, final double baseWeight) throws OrekitException {
126 this(date, position, velocity, sigmaPosition, sigmaVelocity, baseWeight, 0);
127 }
128
129 /** Constructor with two vectors for the standard deviations.
130 * One 3-sized vectors for position standard deviations.
131 * One 3-sized vectors for velocity standard deviations.
132 * The 3-sized vectors are the square root of the diagonal elements of the covariance matrix.
133 * <p>The measurement must be in the orbit propagation frame.</p>
134 * @param date date of the measurement
135 * @param position position
136 * @param velocity velocity
137 * @param sigmaPosition 3-sized vector of the standard deviations of the position
138 * @param sigmaVelocity 3-sized vector of the standard deviations of the velocity
139 * @param baseWeight base weight
140 * @param propagatorIndex index of the propagator related to this measurement
141 * @throws OrekitException if input standard deviations vectors do not have the proper sizes
142 * @since 9.2
143 */
144 public PV(final AbsoluteDate date, final Vector3D position, final Vector3D velocity,
145 final double[] sigmaPosition, final double[] sigmaVelocity,
146 final double baseWeight, final int propagatorIndex) throws OrekitException {
147 this(date, position, velocity,
148 buildPvCovarianceMatrix(sigmaPosition, sigmaVelocity),
149 baseWeight, propagatorIndex);
150 }
151
152 /** Constructor with one vector for the standard deviations and default value for propagator index.
153 * The 6-sized vector is the square root of the diagonal elements of the covariance matrix.
154 * <p>The measurement must be in the orbit propagation frame.</p>
155 * This constructor uses 0 as the index of the propagator related
156 * to this measurement, thus being well suited for mono-satellite
157 * orbit determination.
158 * @param date date of the measurement
159 * @param position position
160 * @param velocity velocity
161 * @param sigmaPV 6-sized vector of the standard deviations
162 * @param baseWeight base weight
163 * @throws OrekitException if input standard deviations vector does not have the proper size
164 * @since 9.2
165 */
166 public PV(final AbsoluteDate date, final Vector3D position, final Vector3D velocity,
167 final double[] sigmaPV, final double baseWeight) throws OrekitException {
168 this(date, position, velocity, sigmaPV, baseWeight, 0);
169 }
170
171 /** Constructor with one vector for the standard deviations.
172 * The 6-sized vector is the square root of the diagonal elements of the covariance matrix.
173 * <p>The measurement must be in the orbit propagation frame.</p>
174 * @param date date of the measurement
175 * @param position position
176 * @param velocity velocity
177 * @param sigmaPV 6-sized vector of the standard deviations
178 * @param baseWeight base weight
179 * @param propagatorIndex index of the propagator related to this measurement
180 * @throws OrekitException if input standard deviations vector does not have the proper size
181 * @since 9.2
182 */
183 public PV(final AbsoluteDate date, final Vector3D position, final Vector3D velocity,
184 final double[] sigmaPV, final double baseWeight, final int propagatorIndex) throws OrekitException {
185 this(date, position, velocity,
186 buildPvCovarianceMatrix(sigmaPV),
187 baseWeight, propagatorIndex);
188 }
189
190 /**
191 * Constructor with 2 smaller covariance matrices and default value for propagator index.
192 * One 3x3 covariance matrix for position and one 3x3 covariance matrix for velocity.
193 * The fact that the covariance matrices are symmetric and positive definite is not checked.
194 * <p>The measurement must be in the orbit propagation frame.</p>
195 * This constructor uses 0 as the index of the propagator related
196 * to this measurement, thus being well suited for mono-satellite
197 * orbit determination.
198 * @param date date of the measurement
199 * @param position position
200 * @param velocity velocity
201 * @param positionCovarianceMatrix 3x3 covariance matrix of the position
202 * @param velocityCovarianceMatrix 3x3 covariance matrix of the velocity
203 * @param baseWeight base weight
204 * @throws OrekitException if input covariance matrices do not have the proper sizes
205 * @since 9.2
206 */
207 public PV(final AbsoluteDate date, final Vector3D position, final Vector3D velocity,
208 final double[][] positionCovarianceMatrix, final double[][] velocityCovarianceMatrix,
209 final double baseWeight) throws OrekitException {
210 this(date, position, velocity, positionCovarianceMatrix, velocityCovarianceMatrix, baseWeight, 0);
211 }
212
213 /**
214 * Constructor with 2 smaller covariance matrices.
215 * One 3x3 covariance matrix for position and one 3x3 covariance matrix for velocity.
216 * The fact that the covariance matrices are symmetric and positive definite is not checked.
217 * <p>The measurement must be in the orbit propagation frame.</p>
218 * @param date date of the measurement
219 * @param position position
220 * @param velocity velocity
221 * @param positionCovarianceMatrix 3x3 covariance matrix of the position
222 * @param velocityCovarianceMatrix 3x3 covariance matrix of the velocity
223 * @param baseWeight base weight
224 * @param propagatorIndex index of the propagator related to this measurement
225 * @throws OrekitException if input covariance matrices do not have the proper sizes
226 * @since 9.2
227 */
228 public PV(final AbsoluteDate date, final Vector3D position, final Vector3D velocity,
229 final double[][] positionCovarianceMatrix, final double[][] velocityCovarianceMatrix,
230 final double baseWeight, final int propagatorIndex) throws OrekitException {
231 this(date, position, velocity,
232 buildPvCovarianceMatrix(positionCovarianceMatrix, velocityCovarianceMatrix),
233 baseWeight, propagatorIndex);
234 }
235
236 /**
237 * Constructor with full covariance matrix but default index for propagator.
238 * The fact that the covariance matrix is symmetric and positive definite is not checked.
239 * <p>The measurement must be in the orbit propagation frame.</p>
240 * This constructor uses 0 as the index of the propagator related
241 * to this measurement, thus being well suited for mono-satellite
242 * orbit determination.
243 * @param date date of the measurement
244 * @param position position
245 * @param velocity velocity
246 * @param covarianceMatrix 6x6 covariance matrix of the PV measurement
247 * @param baseWeight base weight
248 * @throws OrekitException if input covariance matrix does not have the proper size
249 * @since 9.2
250 */
251 public PV(final AbsoluteDate date, final Vector3D position, final Vector3D velocity,
252 final double[][] covarianceMatrix, final double baseWeight) throws OrekitException {
253 this(date, position, velocity, covarianceMatrix, baseWeight, 0);
254 }
255
256 /** Constructor with full covariance matrix and all inputs.
257 * The fact that the covariance matrix is symmetric and positive definite is not checked.
258 * <p>The measurement must be in the orbit propagation frame.</p>
259 * @param date date of the measurement
260 * @param position position
261 * @param velocity velocity
262 * @param covarianceMatrix 6x6 covariance matrix of the PV measurement
263 * @param baseWeight base weight
264 * @param propagatorIndex index of the propagator related to this measurement
265 * @throws OrekitException if input covariance matrix does not have the proper size
266 * @since 9.2
267 */
268 public PV(final AbsoluteDate date, final Vector3D position, final Vector3D velocity,
269 final double[][] covarianceMatrix, final double baseWeight, final int propagatorIndex)
270 throws OrekitException {
271 super(date,
272 new double[] {
273 position.getX(), position.getY(), position.getZ(),
274 velocity.getX(), velocity.getY(), velocity.getZ()
275 }, extractSigmas(covarianceMatrix),
276 new double[] {
277 baseWeight, baseWeight, baseWeight,
278 baseWeight, baseWeight, baseWeight
279 }, Arrays.asList(propagatorIndex));
280 this.covarianceMatrix = covarianceMatrix;
281 }
282
283 /** Get the position.
284 * @return position
285 */
286 public Vector3D getPosition() {
287 final double[] pv = getObservedValue();
288 return new Vector3D(pv[0], pv[1], pv[2]);
289 }
290
291 /** Get the velocity.
292 * @return velocity
293 */
294 public Vector3D getVelocity() {
295 final double[] pv = getObservedValue();
296 return new Vector3D(pv[3], pv[4], pv[5]);
297 }
298
299 /** Get the covariance matrix.
300 * @return the covariance matrix
301 */
302 public double[][] getCovarianceMatrix() {
303 return covarianceMatrix;
304 }
305
306 /** Get the correlation coefficients matrix.
307 * <br>This is the 6x6 matrix M such that:</br>
308 * <br>Mij = Pij/(σi.σj)</br>
309 * <br>Where: <ul>
310 * <li> P is the covariance matrix
311 * <li> σi is the i-th standard deviation (σi² = Pii)
312 * </ul>
313 * @return the correlation coefficient matrix (6x6)
314 */
315 public double[][] getCorrelationCoefficientsMatrix() {
316
317 // Get the standard deviations
318 final double[] sigmas = getTheoreticalStandardDeviation();
319
320 // Initialize the correlation coefficients matric to the covariance matrix
321 final double[][] corrCoefMatrix = new double[sigmas.length][sigmas.length];
322
323 // Divide by the standard deviations
324 for (int i = 0; i < sigmas.length; i++) {
325 for (int j = 0; j < sigmas.length; j++) {
326 corrCoefMatrix[i][j] = covarianceMatrix[i][j] / (sigmas[i] * sigmas[j]);
327 }
328 }
329 return corrCoefMatrix;
330 }
331
332 /** {@inheritDoc} */
333 @Override
334 protected EstimatedMeasurement<PV> theoreticalEvaluation(final int iteration, final int evaluation,
335 final SpacecraftState[] states)
336 throws OrekitException {
337
338 // PV value
339 final TimeStampedPVCoordinates pv = states[getPropagatorsIndices().get(0)].getPVCoordinates();
340
341 // prepare the evaluation
342 final EstimatedMeasurement<PV> estimated =
343 new EstimatedMeasurement<>(this, iteration, evaluation, states,
344 new TimeStampedPVCoordinates[] {
345 pv
346 });
347
348 estimated.setEstimatedValue(new double[] {
349 pv.getPosition().getX(), pv.getPosition().getY(), pv.getPosition().getZ(),
350 pv.getVelocity().getX(), pv.getVelocity().getY(), pv.getVelocity().getZ()
351 });
352
353 // partial derivatives with respect to state
354 estimated.setStateDerivatives(0, IDENTITY);
355
356 return estimated;
357 }
358
359 /** Extract standard deviations from a 6x6 PV covariance matrix.
360 * Check the size of the PV covariance matrix first.
361 * @param pvCovarianceMatrix the 6x6 PV covariance matrix
362 * @return the standard deviations (6-sized vector), they are
363 * the square roots of the diagonal elements of the covariance matrix in input.
364 * @throws OrekitException if the PV covariance matrix is not a 6x6 matrix
365 */
366 private static double[] extractSigmas(final double[][] pvCovarianceMatrix) throws OrekitException {
367
368 // Check the size of the covariance matrix, should be 6x6
369 if (pvCovarianceMatrix[0].length != 6 || pvCovarianceMatrix[1].length != 6) {
370 throw new OrekitException(LocalizedCoreFormats.DIMENSIONS_MISMATCH_2x2,
371 pvCovarianceMatrix[0].length, pvCovarianceMatrix[1],
372 6, 6);
373 }
374
375 // Extract the standard deviations (square roots of the diagonal elements)
376 final double[] sigmas = new double[6];
377 for (int i = 0; i < sigmas.length; i++) {
378 sigmas[i] = FastMath.sqrt(pvCovarianceMatrix[i][i]);
379 }
380 return sigmas;
381 }
382
383 /** Build a 6x6 PV covariance matrix from two 3x3 matrices (covariances in position and velocity).
384 * Check the size of the matrices first.
385 * @param positionCovarianceMatrix the 3x3 covariance matrix in position
386 * @param velocityCovarianceMatrix the 3x3 covariance matrix in velocity
387 * @return the 6x6 PV covariance matrix
388 * @throws OrekitException if the matrices do not have the proper size
389 */
390 private static double[][] buildPvCovarianceMatrix(final double[][] positionCovarianceMatrix,
391 final double[][] velocityCovarianceMatrix)
392 throws OrekitException {
393 // Check the sizes of the matrices first
394 if (positionCovarianceMatrix[0].length != 3 || positionCovarianceMatrix[1].length != 3) {
395 throw new OrekitException(LocalizedCoreFormats.DIMENSIONS_MISMATCH_2x2,
396 positionCovarianceMatrix[0].length, positionCovarianceMatrix[1],
397 3, 3);
398 }
399 if (velocityCovarianceMatrix[0].length != 3 || velocityCovarianceMatrix[1].length != 3) {
400 throw new OrekitException(LocalizedCoreFormats.DIMENSIONS_MISMATCH_2x2,
401 velocityCovarianceMatrix[0].length, velocityCovarianceMatrix[1],
402 3, 3);
403 }
404
405 // Build the PV 6x6 covariance matrix
406 final double[][] pvCovarianceMatrix = new double[6][6];
407 for (int i = 0; i < 3; i++) {
408 for (int j = 0; j < 3; j++) {
409 pvCovarianceMatrix[i][j] = positionCovarianceMatrix[i][j];
410 pvCovarianceMatrix[i + 3][j + 3] = velocityCovarianceMatrix[i][j];
411 }
412 }
413 return pvCovarianceMatrix;
414 }
415
416 /** Build a 6x6 PV covariance matrix from a 6-sized vector (position and velocity standard deviations).
417 * Check the size of the vector first.
418 * @param sigmaPV 6-sized vector with position standard deviations on the first 3 elements
419 * and velocity standard deviations on the last 3 elements
420 * @return the 6x6 PV covariance matrix
421 * @throws OrekitException if the size of the vector is different from 6
422 */
423 private static double[][] buildPvCovarianceMatrix(final double[] sigmaPV) throws OrekitException {
424 // Check the size of the vector first
425 if (sigmaPV.length != 6) {
426 throw new OrekitException(LocalizedCoreFormats.DIMENSIONS_MISMATCH, sigmaPV.length, 6);
427
428 }
429
430 // Build the PV 6x6 covariance matrix
431 final double[][] pvCovarianceMatrix = new double[6][6];
432 for (int i = 0; i < sigmaPV.length; i++) {
433 pvCovarianceMatrix[i][i] = sigmaPV[i] * sigmaPV[i];
434 }
435 return pvCovarianceMatrix;
436 }
437
438 /** Build a 6x6 PV covariance matrix from two 3-sized vectors (position and velocity standard deviations).
439 * Check the sizes of the vectors first.
440 * @param sigmaPosition standard deviations of the position (3-size vector)
441 * @param sigmaVelocity standard deviations of the velocity (3-size vector)
442 * @return the 6x6 PV covariance matrix
443 * @throws OrekitException if the vectors do not have the proper sizes
444 */
445 private static double[][] buildPvCovarianceMatrix(final double[] sigmaPosition,
446 final double[] sigmaVelocity)
447 throws OrekitException {
448
449 // Check the sizes of the vectors first
450 if (sigmaPosition.length != 3) {
451 throw new OrekitException(LocalizedCoreFormats.DIMENSIONS_MISMATCH, sigmaPosition.length, 3);
452
453 }
454 if (sigmaVelocity.length != 3) {
455 throw new OrekitException(LocalizedCoreFormats.DIMENSIONS_MISMATCH, sigmaVelocity.length, 3);
456
457 }
458
459 // Build the PV 6x6 covariance matrix
460 final double[][] pvCovarianceMatrix = new double[6][6];
461 for (int i = 0; i < sigmaPosition.length; i++) {
462 pvCovarianceMatrix[i][i] = sigmaPosition[i] * sigmaPosition[i];
463 pvCovarianceMatrix[i + 3][i + 3] = sigmaVelocity[i] * sigmaVelocity[i];
464 }
465 return pvCovarianceMatrix;
466 }
467 }