1 /* Copyright 2002-2025 CS GROUP 2 * Licensed to CS GROUP (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.orbits; 18 19 import org.hipparchus.geometry.euclidean.threed.Vector3D; 20 import org.hipparchus.linear.DecompositionSolver; 21 import org.hipparchus.linear.MatrixUtils; 22 import org.hipparchus.linear.QRDecomposition; 23 import org.hipparchus.linear.RealMatrix; 24 import org.hipparchus.util.FastMath; 25 import org.hipparchus.util.MathArrays; 26 import org.orekit.errors.OrekitIllegalArgumentException; 27 import org.orekit.errors.OrekitInternalError; 28 import org.orekit.errors.OrekitMessages; 29 import org.orekit.frames.Frame; 30 import org.orekit.frames.StaticTransform; 31 import org.orekit.frames.Transform; 32 import org.orekit.time.AbsoluteDate; 33 import org.orekit.time.TimeOffset; 34 import org.orekit.time.TimeShiftable; 35 import org.orekit.time.TimeStamped; 36 import org.orekit.utils.PVCoordinates; 37 import org.orekit.utils.PVCoordinatesProvider; 38 import org.orekit.utils.TimeStampedPVCoordinates; 39 40 /** 41 * This class handles orbital parameters. 42 43 * <p> 44 * For user convenience, both the Cartesian and the equinoctial elements 45 * are provided by this class, regardless of the canonical representation 46 * implemented in the derived class (which may be classical Keplerian 47 * elements for example). 48 * </p> 49 * <p> 50 * The parameters are defined in a frame specified by the user. It is important 51 * to make sure this frame is consistent: it probably is inertial and centered 52 * on the central body. This information is used for example by some 53 * force models. 54 * </p> 55 * <p> 56 * Instance of this class are guaranteed to be immutable. 57 * </p> 58 * @author Luc Maisonobe 59 * @author Guylaine Prat 60 * @author Fabien Maussion 61 * @author Véronique Pommier-Maurussane 62 */ 63 public abstract class Orbit 64 implements TimeStamped, TimeShiftable<Orbit>, PVCoordinatesProvider { 65 66 /** Absolute tolerance when checking if the rate of the position angle is Keplerian or not. */ 67 protected static final double TOLERANCE_POSITION_ANGLE_RATE = 1e-15; 68 69 /** Frame in which are defined the orbital parameters. */ 70 private final Frame frame; 71 72 /** Date of the orbital parameters. */ 73 private final AbsoluteDate date; 74 75 /** Value of mu used to compute position and velocity (m³/s²). */ 76 private final double mu; 77 78 /** Computed position. 79 * @since 12.0 80 */ 81 private Vector3D position; 82 83 /** Computed PVCoordinates. */ 84 private TimeStampedPVCoordinates pvCoordinates; 85 86 /** Jacobian of the orbital parameters with mean angle with respect to the Cartesian coordinates. */ 87 private double[][] jacobianMeanWrtCartesian; 88 89 /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with mean angle. */ 90 private double[][] jacobianWrtParametersMean; 91 92 /** Jacobian of the orbital parameters with eccentric angle with respect to the Cartesian coordinates. */ 93 private double[][] jacobianEccentricWrtCartesian; 94 95 /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with eccentric angle. */ 96 private double[][] jacobianWrtParametersEccentric; 97 98 /** Jacobian of the orbital parameters with true angle with respect to the Cartesian coordinates. */ 99 private double[][] jacobianTrueWrtCartesian; 100 101 /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with true angle. */ 102 private double[][] jacobianWrtParametersTrue; 103 104 /** Default constructor. 105 * Build a new instance with arbitrary default elements. 106 * @param frame the frame in which the parameters are defined 107 * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame}) 108 * @param date date of the orbital parameters 109 * @param mu central attraction coefficient (m^3/s^2) 110 * @exception IllegalArgumentException if frame is not a {@link 111 * Frame#isPseudoInertial pseudo-inertial frame} 112 */ 113 protected Orbit(final Frame frame, final AbsoluteDate date, final double mu) 114 throws IllegalArgumentException { 115 ensurePseudoInertialFrame(frame); 116 this.date = date; 117 this.mu = mu; 118 this.pvCoordinates = null; 119 this.frame = frame; 120 jacobianMeanWrtCartesian = null; 121 jacobianWrtParametersMean = null; 122 jacobianEccentricWrtCartesian = null; 123 jacobianWrtParametersEccentric = null; 124 jacobianTrueWrtCartesian = null; 125 jacobianWrtParametersTrue = null; 126 } 127 128 /** Set the orbit from Cartesian parameters. 129 * 130 * <p> The acceleration provided in {@code pvCoordinates} is accessible using 131 * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods 132 * use {@code mu} and the position to compute the acceleration, including 133 * {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}. 134 * 135 * @param pvCoordinates the position and velocity in the inertial frame 136 * @param frame the frame in which the {@link TimeStampedPVCoordinates} are defined 137 * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame}) 138 * @param mu central attraction coefficient (m^3/s^2) 139 * @exception IllegalArgumentException if frame is not a {@link 140 * Frame#isPseudoInertial pseudo-inertial frame} 141 */ 142 protected Orbit(final TimeStampedPVCoordinates pvCoordinates, final Frame frame, final double mu) 143 throws IllegalArgumentException { 144 ensurePseudoInertialFrame(frame); 145 this.date = pvCoordinates.getDate(); 146 this.mu = mu; 147 if (pvCoordinates.getAcceleration().getNormSq() == 0) { 148 // the acceleration was not provided, 149 // compute it from Newtonian attraction 150 final double r2 = pvCoordinates.getPosition().getNormSq(); 151 final double r3 = r2 * FastMath.sqrt(r2); 152 this.pvCoordinates = new TimeStampedPVCoordinates(pvCoordinates.getDate(), 153 pvCoordinates.getPosition(), 154 pvCoordinates.getVelocity(), 155 new Vector3D(-mu / r3, pvCoordinates.getPosition())); 156 } else { 157 this.pvCoordinates = pvCoordinates; 158 } 159 this.frame = frame; 160 } 161 162 /** Check if Cartesian coordinates include non-Keplerian acceleration. 163 * @param pva Cartesian coordinates 164 * @param mu central attraction coefficient 165 * @return true if Cartesian coordinates include non-Keplerian acceleration 166 */ 167 protected static boolean hasNonKeplerianAcceleration(final PVCoordinates pva, final double mu) { 168 169 final Vector3D a = pva.getAcceleration(); 170 if (a == null) { 171 return false; 172 } 173 174 final Vector3D p = pva.getPosition(); 175 final double r2 = p.getNormSq(); 176 177 // Check if acceleration is relatively close to 0 compared to the Keplerian acceleration 178 final double tolerance = mu * 1e-9; 179 final Vector3D aTimesR2 = a.scalarMultiply(r2); 180 if (aTimesR2.getNorm() < tolerance) { 181 return false; 182 } 183 184 if ((aTimesR2.add(p.normalize().scalarMultiply(mu))).getNorm() > tolerance) { 185 // we have a relevant acceleration, we can compute derivatives 186 return true; 187 } else { 188 // the provided acceleration is either too small to be reliable (probably even 0), or NaN 189 return false; 190 } 191 } 192 193 /** Returns true if and only if the orbit is elliptical i.e. has a non-negative semi-major axis. 194 * @return true if getA() is strictly greater than 0 195 * @since 12.0 196 */ 197 public boolean isElliptical() { 198 return getA() > 0.; 199 } 200 201 /** Get the orbit type. 202 * @return orbit type 203 */ 204 public abstract OrbitType getType(); 205 206 /** Ensure the defining frame is a pseudo-inertial frame. 207 * @param frame frame to check 208 * @exception IllegalArgumentException if frame is not a {@link 209 * Frame#isPseudoInertial pseudo-inertial frame} 210 */ 211 private static void ensurePseudoInertialFrame(final Frame frame) 212 throws IllegalArgumentException { 213 if (!frame.isPseudoInertial()) { 214 throw new OrekitIllegalArgumentException(OrekitMessages.NON_PSEUDO_INERTIAL_FRAME, 215 frame.getName()); 216 } 217 } 218 219 /** Get the frame in which the orbital parameters are defined. 220 * @return frame in which the orbital parameters are defined 221 */ 222 public Frame getFrame() { 223 return frame; 224 } 225 226 /** Get the semi-major axis. 227 * <p>Note that the semi-major axis is considered negative for hyperbolic orbits.</p> 228 * @return semi-major axis (m) 229 */ 230 public abstract double getA(); 231 232 /** Get the semi-major axis derivative. 233 * <p>Note that the semi-major axis is considered negative for hyperbolic orbits.</p> 234 * <p> 235 * If the orbit was created without derivatives, the value returned is {@link Double#NaN}. 236 * </p> 237 * @return semi-major axis derivative (m/s) 238 * @since 9.0 239 */ 240 public abstract double getADot(); 241 242 /** Get the first component of the equinoctial eccentricity vector. 243 * @return first component of the equinoctial eccentricity vector 244 */ 245 public abstract double getEquinoctialEx(); 246 247 /** Get the first component of the equinoctial eccentricity vector derivative. 248 * <p> 249 * If the orbit was created without derivatives, the value returned is {@link Double#NaN}. 250 * </p> 251 * @return first component of the equinoctial eccentricity vector derivative 252 * @since 9.0 253 */ 254 public abstract double getEquinoctialExDot(); 255 256 /** Get the second component of the equinoctial eccentricity vector. 257 * @return second component of the equinoctial eccentricity vector 258 */ 259 public abstract double getEquinoctialEy(); 260 261 /** Get the second component of the equinoctial eccentricity vector derivative. 262 * <p> 263 * If the orbit was created without derivatives, the value returned is {@link Double#NaN}. 264 * </p> 265 * @return second component of the equinoctial eccentricity vector derivative 266 * @since 9.0 267 */ 268 public abstract double getEquinoctialEyDot(); 269 270 /** Get the first component of the inclination vector. 271 * @return first component of the inclination vector 272 */ 273 public abstract double getHx(); 274 275 /** Get the first component of the inclination vector derivative. 276 * <p> 277 * If the orbit was created without derivatives, the value returned is {@link Double#NaN}. 278 * </p> 279 * @return first component of the inclination vector derivative 280 281 * @since 9.0 282 */ 283 public abstract double getHxDot(); 284 285 /** Get the second component of the inclination vector. 286 * @return second component of the inclination vector 287 */ 288 public abstract double getHy(); 289 290 /** Get the second component of the inclination vector derivative. 291 * <p> 292 * If the orbit was created without derivatives, the value returned is {@link Double#NaN}. 293 * </p> 294 * @return second component of the inclination vector derivative 295 * @since 9.0 296 */ 297 public abstract double getHyDot(); 298 299 /** Get the eccentric longitude argument. 300 * @return E + ω + Ω eccentric longitude argument (rad) 301 */ 302 public abstract double getLE(); 303 304 /** Get the eccentric longitude argument derivative. 305 * <p> 306 * If the orbit was created without derivatives, the value returned is {@link Double#NaN}. 307 * </p> 308 * @return d(E + ω + Ω)/dt eccentric longitude argument derivative (rad/s) 309 * @since 9.0 310 */ 311 public abstract double getLEDot(); 312 313 /** Get the true longitude argument. 314 * @return v + ω + Ω true longitude argument (rad) 315 */ 316 public abstract double getLv(); 317 318 /** Get the true longitude argument derivative. 319 * <p> 320 * If the orbit was created without derivatives, the value returned is {@link Double#NaN}. 321 * </p> 322 * @return d(v + ω + Ω)/dt true longitude argument derivative (rad/s) 323 * @since 9.0 324 */ 325 public abstract double getLvDot(); 326 327 /** Get the mean longitude argument. 328 * @return M + ω + Ω mean longitude argument (rad) 329 */ 330 public abstract double getLM(); 331 332 /** Get the mean longitude argument derivative. 333 * <p> 334 * If the orbit was created without derivatives, the value returned is {@link Double#NaN}. 335 * </p> 336 * @return d(M + ω + Ω)/dt mean longitude argument derivative (rad/s) 337 * @since 9.0 338 */ 339 public abstract double getLMDot(); 340 341 // Additional orbital elements 342 343 /** Get the eccentricity. 344 * @return eccentricity 345 */ 346 public abstract double getE(); 347 348 /** Get the eccentricity derivative. 349 * <p> 350 * If the orbit was created without derivatives, the value returned is {@link Double#NaN}. 351 * </p> 352 * @return eccentricity derivative 353 * @since 9.0 354 */ 355 public abstract double getEDot(); 356 357 /** Get the inclination. 358 * @return inclination (rad) 359 */ 360 public abstract double getI(); 361 362 /** Get the inclination derivative. 363 * <p> 364 * If the orbit was created without derivatives, the value returned is {@link Double#NaN}. 365 * </p> 366 * @return inclination derivative (rad/s) 367 * @since 9.0 368 */ 369 public abstract double getIDot(); 370 371 /** Check if orbit includes non-Keplerian rates. 372 * @return true if orbit includes non-Keplerian derivatives 373 * @see #getADot() 374 * @see #getEquinoctialExDot() 375 * @see #getEquinoctialEyDot() 376 * @see #getHxDot() 377 * @see #getHyDot() 378 * @see #getLEDot() 379 * @see #getLvDot() 380 * @see #getLMDot() 381 * @see #getEDot() 382 * @see #getIDot() 383 * @since 13.0 384 */ 385 public boolean hasNonKeplerianAcceleration() { 386 return hasNonKeplerianAcceleration(getPVCoordinates(), getMu()); 387 } 388 389 /** Get the central acceleration constant. 390 * @return central acceleration constant 391 */ 392 public double getMu() { 393 return mu; 394 } 395 396 /** Get the Keplerian period. 397 * <p>The Keplerian period is computed directly from semi major axis 398 * and central acceleration constant.</p> 399 * @return Keplerian period in seconds, or positive infinity for hyperbolic orbits 400 */ 401 public double getKeplerianPeriod() { 402 final double a = getA(); 403 return isElliptical() ? 2.0 * FastMath.PI * a * FastMath.sqrt(a / mu) : Double.POSITIVE_INFINITY; 404 } 405 406 /** Get the Keplerian mean motion. 407 * <p>The Keplerian mean motion is computed directly from semi major axis 408 * and central acceleration constant.</p> 409 * @return Keplerian mean motion in radians per second 410 */ 411 public double getKeplerianMeanMotion() { 412 final double absA = FastMath.abs(getA()); 413 return FastMath.sqrt(mu / absA) / absA; 414 } 415 416 /** Get the derivative of the mean anomaly with respect to the semi major axis. 417 * @return derivative of the mean anomaly with respect to the semi major axis 418 */ 419 public double getMeanAnomalyDotWrtA() { 420 return -1.5 * getKeplerianMeanMotion() / getA(); 421 } 422 423 /** Get the date of orbital parameters. 424 * @return date of the orbital parameters 425 */ 426 public AbsoluteDate getDate() { 427 return date; 428 } 429 430 /** Get the {@link TimeStampedPVCoordinates} in a specified frame. 431 * @param outputFrame frame in which the position/velocity coordinates shall be computed 432 * @return pvCoordinates in the specified output frame 433 * @see #getPVCoordinates() 434 */ 435 public TimeStampedPVCoordinates getPVCoordinates(final Frame outputFrame) { 436 if (pvCoordinates == null) { 437 pvCoordinates = initPVCoordinates(); 438 } 439 440 // If output frame requested is the same as definition frame, 441 // PV coordinates are returned directly 442 if (outputFrame == frame) { 443 return pvCoordinates; 444 } 445 446 // Else, PV coordinates are transformed to output frame 447 final Transform t = frame.getTransformTo(outputFrame, date); 448 return t.transformPVCoordinates(pvCoordinates); 449 } 450 451 /** {@inheritDoc} */ 452 public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate otherDate, final Frame otherFrame) { 453 return shiftedBy(otherDate.durationFrom(getDate())).getPVCoordinates(otherFrame); 454 } 455 456 /** {@inheritDoc} */ 457 @Override 458 public Vector3D getPosition(final AbsoluteDate otherDate, final Frame otherFrame) { 459 return shiftedBy(otherDate.durationFrom(getDate())).getPosition(otherFrame); 460 } 461 462 /** Get the position in a specified frame. 463 * @param outputFrame frame in which the position coordinates shall be computed 464 * @return position in the specified output frame 465 * @see #getPosition() 466 * @since 12.0 467 */ 468 public Vector3D getPosition(final Frame outputFrame) { 469 if (position == null) { 470 position = initPosition(); 471 } 472 473 // If output frame requested is the same as definition frame, 474 // Position vector is returned directly 475 if (outputFrame == frame) { 476 return position; 477 } 478 479 // Else, position vector is transformed to output frame 480 final StaticTransform t = frame.getStaticTransformTo(outputFrame, date); 481 return t.transformPosition(position); 482 483 } 484 485 /** Get the position in definition frame. 486 * @return position in the definition frame 487 * @see #getPVCoordinates() 488 * @since 12.0 489 */ 490 public Vector3D getPosition() { 491 if (position == null) { 492 position = initPosition(); 493 } 494 return position; 495 } 496 497 /** Get the {@link TimeStampedPVCoordinates} in definition frame. 498 * @return pvCoordinates in the definition frame 499 * @see #getPVCoordinates(Frame) 500 */ 501 public TimeStampedPVCoordinates getPVCoordinates() { 502 if (pvCoordinates == null) { 503 pvCoordinates = initPVCoordinates(); 504 position = pvCoordinates.getPosition(); 505 } 506 return pvCoordinates; 507 } 508 509 /** Compute the position coordinates from the canonical parameters. 510 * @return computed position coordinates 511 * @since 12.0 512 */ 513 protected abstract Vector3D initPosition(); 514 515 /** Compute the position/velocity coordinates from the canonical parameters. 516 * @return computed position/velocity coordinates 517 */ 518 protected abstract TimeStampedPVCoordinates initPVCoordinates(); 519 520 /** 521 * Create a new object representing the same physical orbital state, but attached to a different reference frame. 522 * If the new frame is not inertial, an exception will be thrown. 523 * 524 * @param inertialFrame reference frame of output orbit 525 * @return orbit with different frame 526 * @since 13.0 527 */ 528 public abstract Orbit inFrame(Frame inertialFrame); 529 530 /** Get a time-shifted orbit. 531 * <p> 532 * The orbit can be slightly shifted to close dates. The shifting model is a 533 * Keplerian one if no derivatives are available in the orbit, or Keplerian 534 * plus quadratic effect of the non-Keplerian acceleration if derivatives are 535 * available. Shifting is <em>not</em> intended as a replacement for proper 536 * orbit propagation but should be sufficient for small time shifts or coarse 537 * accuracy. 538 * </p> 539 * @param dt time shift in seconds 540 * @return a new orbit, shifted with respect to the instance (which is immutable) 541 */ 542 @Override 543 public abstract Orbit shiftedBy(double dt); 544 545 /** Get a time-shifted orbit. 546 * <p> 547 * The orbit can be slightly shifted to close dates. The shifting model is a 548 * Keplerian one if no derivatives are available in the orbit, or Keplerian 549 * plus quadratic effect of the non-Keplerian acceleration if derivatives are 550 * available. Shifting is <em>not</em> intended as a replacement for proper 551 * orbit propagation but should be sufficient for small time shifts or coarse 552 * accuracy. 553 * </p> 554 * @param dt time shift 555 * @return a new orbit, shifted with respect to the instance (which is immutable) 556 */ 557 @Override 558 public abstract Orbit shiftedBy(TimeOffset dt); 559 560 /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters. 561 * <p> 562 * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with 563 * respect to Cartesian coordinate j. This means each row corresponds to one orbital parameter 564 * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot. 565 * </p> 566 * @param type type of the position angle to use 567 * @param jacobian placeholder 6x6 (or larger) matrix to be filled with the Jacobian, if matrix 568 * is larger than 6x6, only the 6x6 upper left corner will be modified 569 */ 570 public void getJacobianWrtCartesian(final PositionAngleType type, final double[][] jacobian) { 571 572 final double[][] cachedJacobian; 573 synchronized (this) { 574 switch (type) { 575 case MEAN : 576 if (jacobianMeanWrtCartesian == null) { 577 // first call, we need to compute the Jacobian and cache it 578 jacobianMeanWrtCartesian = computeJacobianMeanWrtCartesian(); 579 } 580 cachedJacobian = jacobianMeanWrtCartesian; 581 break; 582 case ECCENTRIC : 583 if (jacobianEccentricWrtCartesian == null) { 584 // first call, we need to compute the Jacobian and cache it 585 jacobianEccentricWrtCartesian = computeJacobianEccentricWrtCartesian(); 586 } 587 cachedJacobian = jacobianEccentricWrtCartesian; 588 break; 589 case TRUE : 590 if (jacobianTrueWrtCartesian == null) { 591 // first call, we need to compute the Jacobian and cache it 592 jacobianTrueWrtCartesian = computeJacobianTrueWrtCartesian(); 593 } 594 cachedJacobian = jacobianTrueWrtCartesian; 595 break; 596 default : 597 throw new OrekitInternalError(null); 598 } 599 } 600 601 // fill the user provided array 602 for (int i = 0; i < cachedJacobian.length; ++i) { 603 System.arraycopy(cachedJacobian[i], 0, jacobian[i], 0, cachedJacobian[i].length); 604 } 605 606 } 607 608 /** Compute the Jacobian of the Cartesian parameters with respect to the orbital parameters. 609 * <p> 610 * Element {@code jacobian[i][j]} is the derivative of Cartesian coordinate i of the orbit with 611 * respect to orbital parameter j. This means each row corresponds to one Cartesian coordinate 612 * x, y, z, xdot, ydot, zdot whereas columns 0 to 5 correspond to the orbital parameters. 613 * </p> 614 * @param type type of the position angle to use 615 * @param jacobian placeholder 6x6 (or larger) matrix to be filled with the Jacobian, if matrix 616 * is larger than 6x6, only the 6x6 upper left corner will be modified 617 */ 618 public void getJacobianWrtParameters(final PositionAngleType type, final double[][] jacobian) { 619 620 final double[][] cachedJacobian; 621 synchronized (this) { 622 switch (type) { 623 case MEAN : 624 if (jacobianWrtParametersMean == null) { 625 // first call, we need to compute the Jacobian and cache it 626 jacobianWrtParametersMean = createInverseJacobian(type); 627 } 628 cachedJacobian = jacobianWrtParametersMean; 629 break; 630 case ECCENTRIC : 631 if (jacobianWrtParametersEccentric == null) { 632 // first call, we need to compute the Jacobian and cache it 633 jacobianWrtParametersEccentric = createInverseJacobian(type); 634 } 635 cachedJacobian = jacobianWrtParametersEccentric; 636 break; 637 case TRUE : 638 if (jacobianWrtParametersTrue == null) { 639 // first call, we need to compute the Jacobian and cache it 640 jacobianWrtParametersTrue = createInverseJacobian(type); 641 } 642 cachedJacobian = jacobianWrtParametersTrue; 643 break; 644 default : 645 throw new OrekitInternalError(null); 646 } 647 } 648 649 // fill the user-provided array 650 for (int i = 0; i < cachedJacobian.length; ++i) { 651 System.arraycopy(cachedJacobian[i], 0, jacobian[i], 0, cachedJacobian[i].length); 652 } 653 654 } 655 656 /** Create an inverse Jacobian. 657 * @param type type of the position angle to use 658 * @return inverse Jacobian 659 */ 660 private double[][] createInverseJacobian(final PositionAngleType type) { 661 662 // get the direct Jacobian 663 final double[][] directJacobian = new double[6][6]; 664 getJacobianWrtCartesian(type, directJacobian); 665 666 // invert the direct Jacobian 667 final RealMatrix matrix = MatrixUtils.createRealMatrix(directJacobian); 668 final DecompositionSolver solver = new QRDecomposition(matrix).getSolver(); 669 return solver.getInverse().getData(); 670 671 } 672 673 /** Compute the Jacobian of the orbital parameters with mean angle with respect to the Cartesian parameters. 674 * <p> 675 * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with 676 * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter 677 * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot. 678 * </p> 679 * <p> 680 * The array returned by this method will not be modified. 681 * </p> 682 * @return 6x6 Jacobian matrix 683 * @see #computeJacobianEccentricWrtCartesian() 684 * @see #computeJacobianTrueWrtCartesian() 685 */ 686 protected abstract double[][] computeJacobianMeanWrtCartesian(); 687 688 /** Compute the Jacobian of the orbital parameters with eccentric angle with respect to the Cartesian parameters. 689 * <p> 690 * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with 691 * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter 692 * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot. 693 * </p> 694 * <p> 695 * The array returned by this method will not be modified. 696 * </p> 697 * @return 6x6 Jacobian matrix 698 * @see #computeJacobianMeanWrtCartesian() 699 * @see #computeJacobianTrueWrtCartesian() 700 */ 701 protected abstract double[][] computeJacobianEccentricWrtCartesian(); 702 703 /** Compute the Jacobian of the orbital parameters with true angle with respect to the Cartesian parameters. 704 * <p> 705 * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with 706 * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter 707 * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot. 708 * </p> 709 * <p> 710 * The array returned by this method will not be modified. 711 * </p> 712 * @return 6x6 Jacobian matrix 713 * @see #computeJacobianMeanWrtCartesian() 714 * @see #computeJacobianEccentricWrtCartesian() 715 */ 716 protected abstract double[][] computeJacobianTrueWrtCartesian(); 717 718 /** Add the contribution of the Keplerian motion to parameters derivatives 719 * <p> 720 * This method is used by integration-based propagators to evaluate the part of Keplerian 721 * motion to evolution of the orbital state. 722 * </p> 723 * @param type type of the position angle in the state 724 * @param gm attraction coefficient to use 725 * @param pDot array containing orbital state derivatives to update (the Keplerian 726 * part must be <em>added</em> to the array components, as the array may already 727 * contain some non-zero elements corresponding to non-Keplerian parts) 728 */ 729 public abstract void addKeplerContribution(PositionAngleType type, double gm, double[] pDot); 730 731 /** Fill a Jacobian half row with a single vector. 732 * @param a coefficient of the vector 733 * @param v vector 734 * @param row Jacobian matrix row 735 * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set) 736 */ 737 protected static void fillHalfRow(final double a, final Vector3D v, final double[] row, final int j) { 738 row[j] = a * v.getX(); 739 row[j + 1] = a * v.getY(); 740 row[j + 2] = a * v.getZ(); 741 } 742 743 /** Fill a Jacobian half row with a linear combination of vectors. 744 * @param a1 coefficient of the first vector 745 * @param v1 first vector 746 * @param a2 coefficient of the second vector 747 * @param v2 second vector 748 * @param row Jacobian matrix row 749 * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set) 750 */ 751 protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2, 752 final double[] row, final int j) { 753 row[j] = MathArrays.linearCombination(a1, v1.getX(), a2, v2.getX()); 754 row[j + 1] = MathArrays.linearCombination(a1, v1.getY(), a2, v2.getY()); 755 row[j + 2] = MathArrays.linearCombination(a1, v1.getZ(), a2, v2.getZ()); 756 } 757 758 /** Fill a Jacobian half row with a linear combination of vectors. 759 * @param a1 coefficient of the first vector 760 * @param v1 first vector 761 * @param a2 coefficient of the second vector 762 * @param v2 second vector 763 * @param a3 coefficient of the third vector 764 * @param v3 third vector 765 * @param row Jacobian matrix row 766 * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set) 767 */ 768 protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2, 769 final double a3, final Vector3D v3, 770 final double[] row, final int j) { 771 row[j] = MathArrays.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX()); 772 row[j + 1] = MathArrays.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY()); 773 row[j + 2] = MathArrays.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ()); 774 } 775 776 /** Fill a Jacobian half row with a linear combination of vectors. 777 * @param a1 coefficient of the first vector 778 * @param v1 first vector 779 * @param a2 coefficient of the second vector 780 * @param v2 second vector 781 * @param a3 coefficient of the third vector 782 * @param v3 third vector 783 * @param a4 coefficient of the fourth vector 784 * @param v4 fourth vector 785 * @param row Jacobian matrix row 786 * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set) 787 */ 788 protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2, 789 final double a3, final Vector3D v3, final double a4, final Vector3D v4, 790 final double[] row, final int j) { 791 row[j] = MathArrays.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX(), a4, v4.getX()); 792 row[j + 1] = MathArrays.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY(), a4, v4.getY()); 793 row[j + 2] = MathArrays.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ(), a4, v4.getZ()); 794 } 795 796 /** Fill a Jacobian half row with a linear combination of vectors. 797 * @param a1 coefficient of the first vector 798 * @param v1 first vector 799 * @param a2 coefficient of the second vector 800 * @param v2 second vector 801 * @param a3 coefficient of the third vector 802 * @param v3 third vector 803 * @param a4 coefficient of the fourth vector 804 * @param v4 fourth vector 805 * @param a5 coefficient of the fifth vector 806 * @param v5 fifth vector 807 * @param row Jacobian matrix row 808 * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set) 809 */ 810 protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2, 811 final double a3, final Vector3D v3, final double a4, final Vector3D v4, 812 final double a5, final Vector3D v5, 813 final double[] row, final int j) { 814 final double[] a = new double[] { 815 a1, a2, a3, a4, a5 816 }; 817 row[j] = MathArrays.linearCombination(a, new double[] { 818 v1.getX(), v2.getX(), v3.getX(), v4.getX(), v5.getX() 819 }); 820 row[j + 1] = MathArrays.linearCombination(a, new double[] { 821 v1.getY(), v2.getY(), v3.getY(), v4.getY(), v5.getY() 822 }); 823 row[j + 2] = MathArrays.linearCombination(a, new double[] { 824 v1.getZ(), v2.getZ(), v3.getZ(), v4.getZ(), v5.getZ() 825 }); 826 } 827 828 /** Fill a Jacobian half row with a linear combination of vectors. 829 * @param a1 coefficient of the first vector 830 * @param v1 first vector 831 * @param a2 coefficient of the second vector 832 * @param v2 second vector 833 * @param a3 coefficient of the third vector 834 * @param v3 third vector 835 * @param a4 coefficient of the fourth vector 836 * @param v4 fourth vector 837 * @param a5 coefficient of the fifth vector 838 * @param v5 fifth vector 839 * @param a6 coefficient of the sixth vector 840 * @param v6 sixth vector 841 * @param row Jacobian matrix row 842 * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set) 843 */ 844 protected static void fillHalfRow(final double a1, final Vector3D v1, final double a2, final Vector3D v2, 845 final double a3, final Vector3D v3, final double a4, final Vector3D v4, 846 final double a5, final Vector3D v5, final double a6, final Vector3D v6, 847 final double[] row, final int j) { 848 final double[] a = new double[] { 849 a1, a2, a3, a4, a5, a6 850 }; 851 row[j] = MathArrays.linearCombination(a, new double[] { 852 v1.getX(), v2.getX(), v3.getX(), v4.getX(), v5.getX(), v6.getX() 853 }); 854 row[j + 1] = MathArrays.linearCombination(a, new double[] { 855 v1.getY(), v2.getY(), v3.getY(), v4.getY(), v5.getY(), v6.getY() 856 }); 857 row[j + 2] = MathArrays.linearCombination(a, new double[] { 858 v1.getZ(), v2.getZ(), v3.getZ(), v4.getZ(), v5.getZ(), v6.getZ() 859 }); 860 } 861 862 }