1   /* Copyright 2002-2019 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 only measurement.
30   * <p>
31   * For position-velocity measurement see {@link PV}.
32   * </p>
33   * @see PV
34   * @author Luc Maisonobe
35   * @since 9.3
36   */
37  public class Position extends AbstractMeasurement<Position> {
38  
39      /** Identity matrix, for states derivatives. */
40      private static final double[][] IDENTITY = new double[][] {
41          {
42              1, 0, 0, 0, 0, 0
43          }, {
44              0, 1, 0, 0, 0, 0
45          }, {
46              0, 0, 1, 0, 0, 0
47          }
48      };
49  
50      /** Covariance matrix of the PV measurement (size 3x3). */
51      private final double[][] covarianceMatrix;
52  
53      /** Constructor with one double for the standard deviation.
54       * <p>The double is the position's standard deviation, common to the 3 position's components.</p>
55       * <p>
56       * The measurement must be in the orbit propagation frame.
57       * </p>
58       * <p>This constructor uses 0 as the index of the propagator related
59       * to this measurement, thus being well suited for mono-satellite
60       * orbit determination.</p>
61       * @param date date of the measurement
62       * @param position position
63       * @param sigmaPosition theoretical standard deviation on position components
64       * @param baseWeight base weight
65       * @deprecated as of 9.3, replaced by {@link #Position(AbsoluteDate, Vector3D,
66       * double, double, ObservableSatellite)}
67       */
68      @Deprecated
69      public Position(final AbsoluteDate date, final Vector3D position,
70                      final double sigmaPosition, final double baseWeight) {
71          this(date, position, sigmaPosition, baseWeight, new ObservableSatellite(0));
72      }
73  
74      /** Constructor with one double for the standard deviation.
75       * <p>The double is the position's standard deviation, common to the 3 position's components.</p>
76       * <p>
77       * The measurement must be in the orbit propagation frame.
78       * </p>
79       * @param date date of the measurement
80       * @param position position
81       * @param sigmaPosition theoretical standard deviation on position components
82       * @param baseWeight base weight
83       * @param propagatorIndex index of the propagator related to this measurement
84       * @deprecated as of 9.3, replaced by {@link #Position(AbsoluteDate, Vector3D,
85       * double, double, ObservableSatellite)}
86       */
87      @Deprecated
88      public Position(final AbsoluteDate date, final Vector3D position,
89                      final double sigmaPosition, final double baseWeight,
90                      final int propagatorIndex) {
91          this(date, position, sigmaPosition, baseWeight, new ObservableSatellite(propagatorIndex));
92      }
93  
94      /** Constructor with one double for the standard deviation.
95       * <p>The double is the position's standard deviation, common to the 3 position's components.</p>
96       * <p>
97       * The measurement must be in the orbit propagation frame.
98       * </p>
99       * @param date date of the measurement
100      * @param position position
101      * @param sigmaPosition theoretical standard deviation on position components
102      * @param baseWeight base weight
103      * @param satellite satellite related to this measurement
104      * @since 9.3
105      */
106     public Position(final AbsoluteDate date, final Vector3D position,
107                     final double sigmaPosition, final double baseWeight,
108                     final ObservableSatellite satellite) {
109         this(date, position,
110              new double[] {
111                  sigmaPosition,
112                  sigmaPosition,
113                  sigmaPosition
114              }, baseWeight, satellite);
115     }
116 
117     /** Constructor with one vector for the standard deviation and default value for propagator index..
118      * <p>The 3-sized vector represents the square root of the diagonal elements of the covariance matrix.</p>
119      * <p>The measurement must be in the orbit propagation frame.</p>
120      * <p>This constructor uses 0 as the index of the propagator related
121      * to this measurement, thus being well suited for mono-satellite
122      * orbit determination.</p>
123      * @param date date of the measurement
124      * @param position position
125      * @param sigmaPosition 3-sized vector of the standard deviations of the position
126      * @param baseWeight base weight
127      * @deprecated as of 9.3, replaced by {@link #Position(AbsoluteDate, Vector3D,
128      * double[], double, ObservableSatellite)}
129      */
130     @Deprecated
131     public Position(final AbsoluteDate date, final Vector3D position,
132                     final double[] sigmaPosition, final double baseWeight) {
133         this(date, position, sigmaPosition, baseWeight, new ObservableSatellite(0));
134     }
135 
136     /** Constructor with one vector for the standard deviation.
137      * <p>The 3-sized vector represents the square root of the diagonal elements of the covariance matrix.</p>
138      * <p>The measurement must be in the orbit propagation frame.</p>
139      * @param date date of the measurement
140      * @param position position
141      * @param sigmaPosition 3-sized vector of the standard deviations of the position
142      * @param baseWeight base weight
143      * @param propagatorIndex index of the propagator related to this measurement
144      * @deprecated as of 9.3, replaced by {@link #Position(AbsoluteDate, Vector3D,
145      * double[], double, ObservableSatellite)}
146      */
147     @Deprecated
148     public Position(final AbsoluteDate date, final Vector3D position,
149                     final double[] sigmaPosition, final double baseWeight, final int propagatorIndex) {
150         this(date, position, sigmaPosition, baseWeight, new ObservableSatellite(propagatorIndex));
151     }
152 
153     /** Constructor with one vector for the standard deviation.
154      * <p>The 3-sized vector represents the square root of the diagonal elements of the covariance matrix.</p>
155      * <p>The measurement must be in the orbit propagation frame.</p>
156      * @param date date of the measurement
157      * @param position position
158      * @param sigmaPosition 3-sized vector of the standard deviations of the position
159      * @param baseWeight base weight
160      * @param satellite satellite related to this measurement
161      * @since 9.3
162      */
163     public Position(final AbsoluteDate date, final Vector3D position,
164                     final double[] sigmaPosition, final double baseWeight, final ObservableSatellite satellite) {
165         this(date, position, buildPvCovarianceMatrix(sigmaPosition), baseWeight, satellite);
166     }
167 
168     /**
169      * Constructor with a covariance matrix and default value for propagator index.
170      * <p>The fact that the covariance matrices are symmetric and positive definite is not checked.</p>
171      * <p>The measurement must be in the orbit propagation frame.</p>
172      * <p>This constructor uses 0 as the index of the propagator related
173      * to this measurement, thus being well suited for mono-satellite
174      * orbit determination.</p>
175      * @param date date of the measurement
176      * @param position position
177      * @param positionCovarianceMatrix 3x3 covariance matrix of the position
178      * @param baseWeight base weight
179      * @deprecated as of 9.3, replaced by {@link #Position(AbsoluteDate, Vector3D,
180      * double[][], double, ObservableSatellite)}
181      */
182     @Deprecated
183     public Position(final AbsoluteDate date, final Vector3D position,
184                     final double[][] positionCovarianceMatrix, final double baseWeight) {
185         this(date, position, positionCovarianceMatrix, baseWeight, new ObservableSatellite(0));
186     }
187 
188     /** Constructor with full covariance matrix and all inputs.
189      * <p>The fact that the covariance matrix is symmetric and positive definite is not checked.</p>
190      * <p>The measurement must be in the orbit propagation frame.</p>
191      * @param date date of the measurement
192      * @param position position
193      * @param covarianceMatrix 6x6 covariance matrix of the PV measurement
194      * @param baseWeight base weight
195      * @param propagatorIndex index of the propagator related to this measurement
196      * @deprecated as of 9.3, replaced by {@link #Position(AbsoluteDate, Vector3D,
197      * double[][], double, ObservableSatellite)}
198      */
199     @Deprecated
200     public Position(final AbsoluteDate date, final Vector3D position,
201                     final double[][] covarianceMatrix, final double baseWeight,
202                     final int propagatorIndex) {
203         this(date, position, covarianceMatrix, baseWeight, new ObservableSatellite(propagatorIndex));
204     }
205 
206     /** Constructor with full covariance matrix and all inputs.
207      * <p>The fact that the covariance matrix is symmetric and positive definite is not checked.</p>
208      * <p>The measurement must be in the orbit propagation frame.</p>
209      * @param date date of the measurement
210      * @param position position
211      * @param covarianceMatrix 6x6 covariance matrix of the PV measurement
212      * @param baseWeight base weight
213      * @param satellite satellite related to this measurement
214      * @since 9.3
215      */
216     public Position(final AbsoluteDate date, final Vector3D position,
217                     final double[][] covarianceMatrix, final double baseWeight,
218                     final ObservableSatellite satellite) {
219         super(date,
220               new double[] {
221                   position.getX(), position.getY(), position.getZ()
222               }, extractSigmas(covarianceMatrix),
223               new double[] {
224                   baseWeight, baseWeight, baseWeight
225               }, Arrays.asList(satellite));
226         this.covarianceMatrix = covarianceMatrix;
227     }
228 
229     /** Get the position.
230      * @return position
231      */
232     public Vector3D getPosition() {
233         final double[] pv = getObservedValue();
234         return new Vector3D(pv[0], pv[1], pv[2]);
235     }
236 
237     /** Get the covariance matrix.
238      * @return the covariance matrix
239      */
240     public double[][] getCovarianceMatrix() {
241         return covarianceMatrix;
242     }
243 
244     /** Get the correlation coefficients matrix.
245      * <br>This is the 3x3 matrix M such that:</br>
246      * <br>Mij = Pij/(σi.σj)</br>
247      * <br>Where: <ul>
248      * <li> P is the covariance matrix
249      * <li> σi is the i-th standard deviation (σi² = Pii)
250      * </ul>
251      * @return the correlation coefficient matrix (3x3)
252      */
253     public double[][] getCorrelationCoefficientsMatrix() {
254 
255         // Get the standard deviations
256         final double[] sigmas = getTheoreticalStandardDeviation();
257 
258         // Initialize the correlation coefficients matric to the covariance matrix
259         final double[][] corrCoefMatrix = new double[sigmas.length][sigmas.length];
260 
261         // Divide by the standard deviations
262         for (int i = 0; i < sigmas.length; i++) {
263             for (int j = 0; j < sigmas.length; j++) {
264                 corrCoefMatrix[i][j] = covarianceMatrix[i][j] / (sigmas[i] * sigmas[j]);
265             }
266         }
267         return corrCoefMatrix;
268     }
269 
270     /** {@inheritDoc} */
271     @Override
272     protected EstimatedMeasurement<Position> theoreticalEvaluation(final int iteration, final int evaluation,
273                                                                    final SpacecraftState[] states) {
274 
275         // PV value
276         final ObservableSatellite      satellite = getSatellites().get(0);
277         final SpacecraftState          state     = states[satellite.getPropagatorIndex()];
278         final TimeStampedPVCoordinates pv        = state.getPVCoordinates();
279 
280         // prepare the evaluation
281         final EstimatedMeasurement<Position> estimated =
282                         new EstimatedMeasurement<>(this, iteration, evaluation, states,
283                                                    new TimeStampedPVCoordinates[] {
284                                                        pv
285                                                    });
286 
287         estimated.setEstimatedValue(new double[] {
288             pv.getPosition().getX(), pv.getPosition().getY(), pv.getPosition().getZ()
289         });
290 
291         // partial derivatives with respect to state
292         estimated.setStateDerivatives(0, IDENTITY);
293 
294         return estimated;
295     }
296 
297     /** Extract standard deviations from a 3x3 position covariance matrix.
298      * Check the size of the position covariance matrix first.
299      * @param pCovarianceMatrix the 3x" possition covariance matrix
300      * @return the standard deviations (3-sized vector), they are
301      * the square roots of the diagonal elements of the covariance matrix in input.
302      */
303     private static double[] extractSigmas(final double[][] pCovarianceMatrix) {
304 
305         // Check the size of the covariance matrix, should be 3x3
306         if (pCovarianceMatrix.length != 3 || pCovarianceMatrix[0].length != 3) {
307             throw new OrekitException(LocalizedCoreFormats.DIMENSIONS_MISMATCH_2x2,
308                                       pCovarianceMatrix.length, pCovarianceMatrix[0],
309                                       3, 3);
310         }
311 
312         // Extract the standard deviations (square roots of the diagonal elements)
313         final double[] sigmas = new double[3];
314         for (int i = 0; i < sigmas.length; i++) {
315             sigmas[i] = FastMath.sqrt(pCovarianceMatrix[i][i]);
316         }
317         return sigmas;
318     }
319 
320     /** Build a 3x3 position covariance matrix from a 3-sized vector (position standard deviations).
321      * Check the size of the vector first.
322      * @param sigmaP 3-sized vector with position standard deviations
323      * @return the 3x3 position covariance matrix
324      */
325     private static double[][] buildPvCovarianceMatrix(final double[] sigmaP) {
326         // Check the size of the vector first
327         if (sigmaP.length != 3) {
328             throw new OrekitException(LocalizedCoreFormats.DIMENSIONS_MISMATCH, sigmaP.length, 3);
329 
330         }
331 
332         // Build the 3x3 position covariance matrix
333         final double[][] pvCovarianceMatrix = new double[3][3];
334         for (int i = 0; i < sigmaP.length; i++) {
335             pvCovarianceMatrix[i][i] =  sigmaP[i] * sigmaP[i];
336         }
337         return pvCovarianceMatrix;
338     }
339 
340 }