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