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 }