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 }