FieldOrbit.java

  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.orbits;



  21. import org.hipparchus.RealFieldElement;
  22. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  23. import org.hipparchus.linear.FieldDecompositionSolver;
  24. import org.hipparchus.linear.FieldLUDecomposition;
  25. import org.hipparchus.linear.FieldMatrix;
  26. import org.hipparchus.linear.MatrixUtils;
  27. import org.hipparchus.util.FastMath;
  28. import org.hipparchus.util.MathArrays;
  29. import org.orekit.errors.OrekitException;
  30. import org.orekit.errors.OrekitIllegalArgumentException;
  31. import org.orekit.errors.OrekitInternalError;
  32. import org.orekit.errors.OrekitMessages;
  33. import org.orekit.frames.Frame;
  34. import org.orekit.frames.Transform;
  35. import org.orekit.time.FieldAbsoluteDate;
  36. import org.orekit.time.FieldTimeInterpolable;
  37. import org.orekit.time.FieldTimeShiftable;
  38. import org.orekit.time.FieldTimeStamped;
  39. import org.orekit.utils.FieldPVCoordinates;
  40. import org.orekit.utils.FieldPVCoordinatesProvider;
  41. import org.orekit.utils.TimeStampedFieldPVCoordinates;
  42. import org.orekit.utils.TimeStampedPVCoordinates;

  44. /**
  45.  * This class handles orbital parameters.

  47.  * <p>
  48.  * For user convenience, both the Cartesian and the equinoctial elements
  49.  * are provided by this class, regardless of the canonical representation
  50.  * implemented in the derived class (which may be classical Keplerian
  51.  * elements for example).
  52.  * </p>
  53.  * <p>
  54.  * The parameters are defined in a frame specified by the user. It is important
  55.  * to make sure this frame is consistent: it probably is inertial and centered
  56.  * on the central body. This information is used for example by some
  57.  * force models.
  58.  * </p>
  59.  * <p>
  60.  * Instance of this class are guaranteed to be immutable.
  61.  * </p>
  62.  * @author Luc Maisonobe
  63.  * @author Guylaine Prat
  64.  * @author Fabien Maussion
  65.  * @author V&eacute;ronique Pommier-Maurussane
  66.  * @since 9.0
  67.  */
  68. public abstract class FieldOrbit<T extends RealFieldElement<T>>
  69.     implements FieldPVCoordinatesProvider<T>, FieldTimeStamped<T>, FieldTimeShiftable<FieldOrbit<T>, T>, FieldTimeInterpolable<FieldOrbit<T>, T> {

  71.     /** Frame in which are defined the orbital parameters. */
  72.     private final Frame frame;

  74.     /** Date of the orbital parameters. */
  75.     private final FieldAbsoluteDate<T> date;

  77.     /** Value of mu used to compute position and velocity (m³/s²). */
  78.     private final double mu;

  80.     /** Computed PVCoordinates. */
  81.     private transient TimeStampedFieldPVCoordinates<T> pvCoordinates;

  83.     /** Jacobian of the orbital parameters with mean angle with respect to the Cartesian coordinates. */
  84.     private transient T[][] jacobianMeanWrtCartesian;

  86.     /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with mean angle. */
  87.     private transient T[][] jacobianWrtParametersMean;

  89.     /** Jacobian of the orbital parameters with eccentric angle with respect to the Cartesian coordinates. */
  90.     private transient T[][] jacobianEccentricWrtCartesian;

  92.     /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with eccentric angle. */
  93.     private transient T[][] jacobianWrtParametersEccentric;

  95.     /** Jacobian of the orbital parameters with true angle with respect to the Cartesian coordinates. */
  96.     private transient T[][] jacobianTrueWrtCartesian;

  98.     /** Jacobian of the Cartesian coordinates with respect to the orbital parameters with true angle. */
  99.     private transient T[][] jacobianWrtParametersTrue;

  101.     /** Default constructor.
  102.      * Build a new instance with arbitrary default elements.
  103.      * @param frame the frame in which the parameters are defined
  104.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  105.      * @param date date of the orbital parameters
  106.      * @param mu central attraction coefficient (m^3/s^2)
  107.      * @exception IllegalArgumentException if frame is not a {@link
  108.      * Frame#isPseudoInertial pseudo-inertial frame}
  109.      */
  110.     protected FieldOrbit(final Frame frame, final FieldAbsoluteDate<T> date, final double mu)
  111.         throws IllegalArgumentException {
  112.         ensurePseudoInertialFrame(frame);
  113.         this.date                      = date;
  114.         this.mu                        = mu;
  115.         this.pvCoordinates             = null;
  116.         this.frame                     = frame;
  117.         jacobianMeanWrtCartesian       = null;
  118.         jacobianWrtParametersMean      = null;
  119.         jacobianEccentricWrtCartesian  = null;
  120.         jacobianWrtParametersEccentric = null;
  121.         jacobianTrueWrtCartesian       = null;
  122.         jacobianWrtParametersTrue      = null;
  123.     }

  125.     /** Set the orbit from Cartesian parameters.
  126.      *
  127.      * <p> The acceleration provided in {@code FieldPVCoordinates} is accessible using
  128.      * {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
  129.      * use {@code mu} and the position to compute the acceleration, including
  130.      * {@link #shiftedBy(RealFieldElement)} and {@link #getPVCoordinates(FieldAbsoluteDate, Frame)}.
  131.      *
  132.      * @param FieldPVCoordinates the position and velocity in the inertial frame
  133.      * @param frame the frame in which the {@link TimeStampedPVCoordinates} are defined
  134.      * (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
  135.      * @param mu central attraction coefficient (m^3/s^2)
  136.      * @exception IllegalArgumentException if frame is not a {@link
  137.      * Frame#isPseudoInertial pseudo-inertial frame}
  138.      */
  139.     protected FieldOrbit(final TimeStampedFieldPVCoordinates<T> FieldPVCoordinates, final Frame frame, final double mu)
  140.         throws IllegalArgumentException {
  141.         ensurePseudoInertialFrame(frame);
  142.         this.date = FieldPVCoordinates.getDate();
  143.         this.mu = mu;
  144.         if (FieldPVCoordinates.getAcceleration().getNormSq().getReal() == 0.0) {
  145.             // the acceleration was not provided,
  146.             // compute it from Newtonian attraction
  147.             final T r2 = FieldPVCoordinates.getPosition().getNormSq();
  148.             final T r3 = r2.multiply(r2.sqrt());
  149.             this.pvCoordinates = new TimeStampedFieldPVCoordinates<>(FieldPVCoordinates.getDate(),
  150.                                                                      FieldPVCoordinates.getPosition(),
  151.                                                                      FieldPVCoordinates.getVelocity(),
  152.                                                                      new FieldVector3D<>(r3.reciprocal().multiply(-mu), FieldPVCoordinates.getPosition()));
  153.         } else {
  154.             this.pvCoordinates = FieldPVCoordinates;
  155.         }
  156.         this.frame = frame;
  157.     }

  159.     /** Check if Cartesian coordinates include non-Keplerian acceleration.
  160.      * @param pva Cartesian coordinates
  161.      * @param mu central attraction coefficient
  162.      * @param <T> type of the field elements
  163.      * @return true if Cartesian coordinates include non-Keplerian acceleration
  164.      */
  165.     protected static <T extends RealFieldElement<T>> boolean hasNonKeplerianAcceleration(final FieldPVCoordinates<T> pva, final double mu) {

  167.         final FieldVector3D<T> a = pva.getAcceleration();
  168.         if (a == null) {
  169.             return false;
  170.         }

  172.         final FieldVector3D<T> p = pva.getPosition();
  173.         final T r2 = p.getNormSq();
  174.         final T r  = r2.sqrt();
  175.         final FieldVector3D<T> keplerianAcceleration = new FieldVector3D<>(r.multiply(r2).reciprocal().multiply(-mu), p);
  176.         if (a.getNorm().getReal() > 1.0e-9 * keplerianAcceleration.getNorm().getReal()) {
  177.             // we have a relevant acceleration, we can compute derivatives
  178.             return true;
  179.         } else {
  180.             // the provided acceleration is either too small to be reliable (probably even 0), or NaN
  181.             return false;
  182.         }

  184.     }

  186.     /** Get the orbit type.
  187.      * @return orbit type
  188.      */
  189.     public abstract OrbitType getType();

  191.     /** Ensure the defining frame is a pseudo-inertial frame.
  192.      * @param frame frame to check
  193.      * @exception IllegalArgumentException if frame is not a {@link
  194.      * Frame#isPseudoInertial pseudo-inertial frame}
  195.      */
  196.     private static void ensurePseudoInertialFrame(final Frame frame)
  197.         throws IllegalArgumentException {
  198.         if (!frame.isPseudoInertial()) {
  199.             throw new OrekitIllegalArgumentException(OrekitMessages.NON_PSEUDO_INERTIAL_FRAME,
  200.                                                      frame.getName());
  201.         }
  202.     }

  204.     /** Get the frame in which the orbital parameters are defined.
  205.      * @return frame in which the orbital parameters are defined
  206.      */
  207.     public Frame getFrame() {
  208.         return frame;
  209.     }

  211.     /**Transforms the FieldOrbit instance into an Orbit instance.
  212.      * @return Orbit instance with same properties
  213.      */
  214.     public abstract Orbit toOrbit();

  216.     /** Get the semi-major axis.
  217.      * <p>Note that the semi-major axis is considered negative for hyperbolic orbits.</p>
  218.      * @return semi-major axis (m)
  219.      */
  220.     public abstract T getA();

  222.     /** Get the semi-major axis derivative.
  223.      * <p>Note that the semi-major axis is considered negative for hyperbolic orbits.</p>
  224.      * <p>
  225.      * If the orbit was created without derivatives, the value returned is null.
  226.      * </p>
  227.      * @return semi-major axis  derivative (m/s)
  228.      */
  229.     public abstract T getADot();

  231.     /** Get the first component of the equinoctial eccentricity vector.
  232.      * @return first component of the equinoctial eccentricity vector
  233.      */
  234.     public abstract T getEquinoctialEx();

  236.     /** Get the first component of the equinoctial eccentricity vector.
  237.      * <p>
  238.      * If the orbit was created without derivatives, the value returned is null.
  239.      * </p>
  240.      * @return first component of the equinoctial eccentricity vector
  241.      */
  242.     public abstract T getEquinoctialExDot();

  244.     /** Get the second component of the equinoctial eccentricity vector.
  245.      * @return second component of the equinoctial eccentricity vector
  246.      */
  247.     public abstract T getEquinoctialEy();

  249.     /** Get the second component of the equinoctial eccentricity vector.
  250.      * <p>
  251.      * If the orbit was created without derivatives, the value returned is null.
  252.      * </p>
  253.      * @return second component of the equinoctial eccentricity vector
  254.      */
  255.     public abstract T getEquinoctialEyDot();

  257.     /** Get the first component of the inclination vector.
  258.      * @return first component of the inclination vector
  259.      */
  260.     public abstract T getHx();

  262.     /** Get the first component of the inclination vector derivative.
  263.      * <p>
  264.      * If the orbit was created without derivatives, the value returned is null.
  265.      * </p>
  266.      * @return first component of the inclination vector derivative
  267.      */
  268.     public abstract T getHxDot();

  270.     /** Get the second component of the inclination vector.
  271.      * @return second component of the inclination vector
  272.      */
  273.     public abstract T getHy();

  275.     /** Get the second component of the inclination vector derivative.
  276.      * <p>
  277.      * </p>
  278.      * @return second component of the inclination vector derivative
  279.      */
  280.     public abstract T getHyDot();

  282.     /** Get the eccentric longitude argument.
  283.      * @return E + ω + Ω eccentric longitude argument (rad)
  284.      */
  285.     public abstract T getLE();

  287.     /** Get the eccentric longitude argument derivative.
  288.      * <p>
  289.      * If the orbit was created without derivatives, the value returned is null.
  290.      * </p>
  291.      * @return d(E + ω + Ω)/dt eccentric longitude argument derivative (rad/s)
  292.      */
  293.     public abstract T getLEDot();

  295.     /** Get the true longitude argument.
  296.      * @return v + ω + Ω true longitude argument (rad)
  297.      */
  298.     public abstract T getLv();

  300.     /** Get the true longitude argument derivative.
  301.      * <p>
  302.      * If the orbit was created without derivatives, the value returned is null.
  303.      * </p>
  304.      * @return d(v + ω + Ω)/dt true longitude argument derivative (rad/s)
  305.      */
  306.     public abstract T getLvDot();

  308.     /** Get the mean longitude argument.
  309.      * @return M + ω + Ω mean longitude argument (rad)
  310.      */
  311.     public abstract T getLM();

  313.     /** Get the mean longitude argument derivative.
  314.      * <p>
  315.      * If the orbit was created without derivatives, the value returned is null.
  316.      * </p>
  317.      * @return d(M + ω + Ω)/dt mean longitude argument derivative (rad/s)
  318.      */
  319.     public abstract T getLMDot();

  321.     // Additional orbital elements

  323.     /** Get the eccentricity.
  324.      * @return eccentricity
  325.      */
  326.     public abstract T getE();

  328.     /** Get the eccentricity derivative.
  329.      * <p>
  330.      * If the orbit was created without derivatives, the value returned is null.
  331.      * </p>
  332.      * @return eccentricity derivative
  333.      */
  334.     public abstract T getEDot();

  336.     /** Get the inclination.
  337.      * <p>
  338.      * If the orbit was created without derivatives, the value returned is null.
  339.      * </p>
  340.      * @return inclination (rad)
  341.      */
  342.     public abstract T getI();

  344.     /** Get the inclination derivative.
  345.      * @return inclination derivative (rad/s)
  346.      */
  347.     public abstract T getIDot();

  349.     /** Get the central acceleration constant.
  350.      * @return central acceleration constant
  351.      */

  353.     /** Check if orbit includes derivatives.
  354.      * @return true if orbit includes derivatives
  355.      * @see #getADot()
  356.      * @see #getEquinoctialExDot()
  357.      * @see #getEquinoctialEyDot()
  358.      * @see #getHxDot()
  359.      * @see #getHyDot()
  360.      * @see #getLEDot()
  361.      * @see #getLvDot()
  362.      * @see #getLMDot()
  363.      * @see #getEDot()
  364.      * @see #getIDot()
  365.      * @since 9.0
  366.      */
  367.     public abstract boolean hasDerivatives();

  369.     public double getMu() {
  370.         return mu;
  371.     }

  373.     /** Get the Keplerian period.
  374.      * <p>The Keplerian period is computed directly from semi major axis
  375.      * and central acceleration constant.</p>
  376.      * @return Keplerian period in seconds, or positive infinity for hyperbolic orbits
  377.      */
  378.     public T getKeplerianPeriod() {
  379.         final T a = getA();
  380.         return (a.getReal() < 0) ? getA().getField().getZero().add(Double.POSITIVE_INFINITY) : a.multiply(2 * FastMath.PI).multiply(a.divide(mu).sqrt());
  381.     }

  383.     /** Get the Keplerian mean motion.
  384.      * <p>The Keplerian mean motion is computed directly from semi major axis
  385.      * and central acceleration constant.</p>
  386.      * @return Keplerian mean motion in radians per second
  387.      */
  388.     public T getKeplerianMeanMotion() {
  389.         final T absA = getA().abs();
  390.         return absA.reciprocal().multiply(mu).sqrt().divide(absA);
  391.     }

  393.     /** Get the date of orbital parameters.
  394.      * @return date of the orbital parameters
  395.      */
  396.     public FieldAbsoluteDate<T> getDate() {
  397.         return date;
  398.     }

  400.     /** Get the {@link TimeStampedPVCoordinates} in a specified frame.
  401.      * @param outputFrame frame in which the position/velocity coordinates shall be computed
  402.      * @return FieldPVCoordinates in the specified output frame
  403.      * @exception OrekitException if transformation between frames cannot be computed
  404.      * @see #getPVCoordinates()
  405.      */
  406.     public TimeStampedFieldPVCoordinates<T> getPVCoordinates(final Frame outputFrame)
  407.         throws OrekitException {
  408.         if (pvCoordinates == null) {
  409.             pvCoordinates = initPVCoordinates();
  410.         }

  412.         // If output frame requested is the same as definition frame,
  413.         // PV coordinates are returned directly
  414.         if (outputFrame == frame) {
  415.             return pvCoordinates;
  416.         }

  418.         // Else, PV coordinates are transformed to output frame
  419.         final Transform t = frame.getTransformTo(outputFrame, date.toAbsoluteDate()); //TODO CHECK THIS
  420.         return t.transformPVCoordinates(pvCoordinates);
  421.     }

  423.     /** {@inheritDoc} */
  424.     public TimeStampedFieldPVCoordinates<T> getPVCoordinates(final FieldAbsoluteDate<T> otherDate, final Frame otherFrame)
  425.         throws OrekitException {
  426.         return shiftedBy(otherDate.durationFrom(getDate())).getPVCoordinates(otherFrame);
  427.     }


  430.     /** Get the {@link TimeStampedPVCoordinates} in definition frame.
  431.      * @return FieldPVCoordinates in the definition frame
  432.      * @see #getPVCoordinates(Frame)
  433.      */
  434.     public TimeStampedFieldPVCoordinates<T> getPVCoordinates() {
  435.         if (pvCoordinates == null) {
  436.             pvCoordinates = initPVCoordinates();

  438.         }
  439.         return pvCoordinates;
  440.     }

  442.     /** Compute the position/velocity coordinates from the canonical parameters.
  443.      * @return computed position/velocity coordinates
  444.      */
  445.     protected abstract TimeStampedFieldPVCoordinates<T> initPVCoordinates();

  447.     /** Get a time-shifted orbit.
  448.      * <p>
  449.      * The orbit can be slightly shifted to close dates. This shift is based on
  450.      * a simple Keplerian model. It is <em>not</em> intended as a replacement
  451.      * for proper orbit and attitude propagation but should be sufficient for
  452.      * small time shifts or coarse accuracy.
  453.      * </p>
  454.      * @param dt time shift in seconds
  455.      * @return a new orbit, shifted with respect to the instance (which is immutable)
  456.      */
  457.     public abstract FieldOrbit<T> shiftedBy(T dt);

  459.     /** Compute the Jacobian of the orbital parameters with respect to the Cartesian parameters.
  460.      * <p>
  461.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  462.      * respect to Cartesian coordinate j. This means each row corresponds to one orbital parameter
  463.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  464.      * </p>
  465.      * @param type type of the position angle to use
  466.      * @param jacobian placeholder 6x6 (or larger) matrix to be filled with the Jacobian, if matrix
  467.      * is larger than 6x6, only the 6x6 upper left corner will be modified
  468.      */
  469.     public void getJacobianWrtCartesian(final PositionAngle type, final T[][] jacobian) {

  471.         final T[][] cachedJacobian;
  472.         synchronized (this) {
  473.             switch (type) {
  474.                 case MEAN :
  475.                     if (jacobianMeanWrtCartesian == null) {
  476.                         // first call, we need to compute the Jacobian and cache it
  477.                         jacobianMeanWrtCartesian = computeJacobianMeanWrtCartesian();
  478.                     }
  479.                     cachedJacobian = jacobianMeanWrtCartesian;
  480.                     break;
  481.                 case ECCENTRIC :
  482.                     if (jacobianEccentricWrtCartesian == null) {
  483.                         // first call, we need to compute the Jacobian and cache it
  484.                         jacobianEccentricWrtCartesian = computeJacobianEccentricWrtCartesian();
  485.                     }
  486.                     cachedJacobian = jacobianEccentricWrtCartesian;
  487.                     break;
  488.                 case TRUE :
  489.                     if (jacobianTrueWrtCartesian == null) {
  490.                         // first call, we need to compute the Jacobian and cache it
  491.                         jacobianTrueWrtCartesian = computeJacobianTrueWrtCartesian();
  492.                     }
  493.                     cachedJacobian = jacobianTrueWrtCartesian;
  494.                     break;
  495.                 default :
  496.                     throw new OrekitInternalError(null);
  497.             }
  498.         }

  500.         // fill the user provided array
  501.         for (int i = 0; i < cachedJacobian.length; ++i) {
  502.             System.arraycopy(cachedJacobian[i], 0, jacobian[i], 0, cachedJacobian[i].length);
  503.         }

  505.     }

  507.     /** Compute the Jacobian of the Cartesian parameters with respect to the orbital parameters.
  508.      * <p>
  509.      * Element {@code jacobian[i][j]} is the derivative of Cartesian coordinate i of the orbit with
  510.      * respect to orbital parameter j. This means each row corresponds to one Cartesian coordinate
  511.      * x, y, z, xdot, ydot, zdot whereas columns 0 to 5 correspond to the orbital parameters.
  512.      * </p>
  513.      * @param type type of the position angle to use
  514.      * @param jacobian placeholder 6x6 (or larger) matrix to be filled with the Jacobian, if matrix
  515.      * is larger than 6x6, only the 6x6 upper left corner will be modified
  516.      */
  517.     public void getJacobianWrtParameters(final PositionAngle type, final T[][] jacobian) {

  519.         final T[][] cachedJacobian;
  520.         synchronized (this) {
  521.             switch (type) {
  522.                 case MEAN :
  523.                     if (jacobianWrtParametersMean == null) {
  524.                         // first call, we need to compute the Jacobian and cache it
  525.                         jacobianWrtParametersMean = createInverseJacobian(type);
  526.                     }
  527.                     cachedJacobian = jacobianWrtParametersMean;
  528.                     break;
  529.                 case ECCENTRIC :
  530.                     if (jacobianWrtParametersEccentric == null) {
  531.                         // first call, we need to compute the Jacobian and cache it
  532.                         jacobianWrtParametersEccentric = createInverseJacobian(type);
  533.                     }
  534.                     cachedJacobian = jacobianWrtParametersEccentric;
  535.                     break;
  536.                 case TRUE :
  537.                     if (jacobianWrtParametersTrue == null) {
  538.                         // first call, we need to compute the Jacobian and cache it
  539.                         jacobianWrtParametersTrue = createInverseJacobian(type);
  540.                     }
  541.                     cachedJacobian = jacobianWrtParametersTrue;
  542.                     break;
  543.                 default :
  544.                     throw new OrekitInternalError(null);
  545.             }
  546.         }

  548.         // fill the user-provided array
  549.         for (int i = 0; i < cachedJacobian.length; ++i) {
  550.             System.arraycopy(cachedJacobian[i], 0, jacobian[i], 0, cachedJacobian[i].length);
  551.         }

  553.     }

  555.     /** Create an inverse Jacobian.
  556.      * @param type type of the position angle to use
  557.      * @return inverse Jacobian
  558.      */
  559.     private T[][] createInverseJacobian(final PositionAngle type) {

  561.         // get the direct Jacobian
  562.         final T[][] directJacobian = MathArrays.buildArray(getA().getField(), 6, 6);
  563.         getJacobianWrtCartesian(type, directJacobian);

  565.         // invert the direct Jacobian
  566.         final FieldMatrix<T> matrix = MatrixUtils.createFieldMatrix(directJacobian);
  567.         final FieldDecompositionSolver<T> solver = new FieldLUDecomposition<>(matrix).getSolver();
  568.         return solver.getInverse().getData();

  570.     }

  572.     /** Compute the Jacobian of the orbital parameters with mean angle with respect to the Cartesian parameters.
  573.      * <p>
  574.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  575.      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
  576.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  577.      * </p>
  578.      * @return 6x6 Jacobian matrix
  579.      * @see #computeJacobianEccentricWrtCartesian()
  580.      * @see #computeJacobianTrueWrtCartesian()
  581.      */
  582.     protected abstract T[][] computeJacobianMeanWrtCartesian();

  584.     /** Compute the Jacobian of the orbital parameters with eccentric angle with respect to the Cartesian parameters.
  585.      * <p>
  586.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  587.      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
  588.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  589.      * </p>
  590.      * @return 6x6 Jacobian matrix
  591.      * @see #computeJacobianMeanWrtCartesian()
  592.      * @see #computeJacobianTrueWrtCartesian()
  593.      */
  594.     protected abstract T[][] computeJacobianEccentricWrtCartesian();

  596.     /** Compute the Jacobian of the orbital parameters with true angle with respect to the Cartesian parameters.
  597.      * <p>
  598.      * Element {@code jacobian[i][j]} is the derivative of parameter i of the orbit with
  599.      * respect to Cartesian coordinate j. This means each row correspond to one orbital parameter
  600.      * whereas columns 0 to 5 correspond to the Cartesian coordinates x, y, z, xDot, yDot and zDot.
  601.      * </p>
  602.      * @return 6x6 Jacobian matrix
  603.      * @see #computeJacobianMeanWrtCartesian()
  604.      * @see #computeJacobianEccentricWrtCartesian()
  605.      */
  606.     protected abstract T[][] computeJacobianTrueWrtCartesian();

  608.     /** Add the contribution of the Keplerian motion to parameters derivatives
  609.      * <p>
  610.      * This method is used by integration-based propagators to evaluate the part of Keplerian
  611.      * motion to evolution of the orbital state.
  612.      * </p>
  613.      * @param type type of the position angle in the state
  614.      * @param gm attraction coefficient to use
  615.      * @param pDot array containing orbital state derivatives to update (the Keplerian
  616.      * part must be <em>added</em> to the array components, as the array may already
  617.      * contain some non-zero elements corresponding to non-Keplerian parts)
  618.      */
  619.     public abstract void addKeplerContribution(PositionAngle type, double gm, T[] pDot);

  621.         /** Fill a Jacobian half row with a single vector.
  622.      * @param a coefficient of the vector
  623.      * @param v vector
  624.      * @param row Jacobian matrix row
  625.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  626.      */

  628.     protected void fillHalfRow(final T a, final FieldVector3D<T> v, final T[] row, final int j) {
  629.         row[j]     = a.multiply(v.getX());
  630.         row[j + 1] = a.multiply(v.getY());
  631.         row[j + 2] = a.multiply(v.getZ());
  632.     }

  634.     /** Fill a Jacobian half row with a linear combination of vectors.
  635.      * @param a1 coefficient of the first vector
  636.      * @param v1 first vector
  637.      * @param a2 coefficient of the second vector
  638.      * @param v2 second vector
  639.      * @param row Jacobian matrix row
  640.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  641.      */
  642.     protected void fillHalfRow(final T a1, final FieldVector3D<T> v1, final T a2, final FieldVector3D<T> v2,
  643.                                       final T[] row, final int j) {
  644.         row[j]     = a1.linearCombination(a1, v1.getX(), a2, v2.getX());
  645.         row[j + 1] = a1.linearCombination(a1, v1.getY(), a2, v2.getY());
  646.         row[j + 2] = a1.linearCombination(a1, v1.getZ(), a2, v2.getZ());
  647. //        row[j]     = a1.multiply(v1.getX()).add(a2.multiply(v2.getX()));
  648. //        row[j + 1] = a1.multiply(v1.getY()).add(a2.multiply(v2.getY()));
  649. //        row[j + 2] = a1.multiply(v1.getZ()).add(a2.multiply(v2.getZ()));
  650.     }

  652.     /** Fill a Jacobian half row with a linear combination of vectors.
  653.      * @param a1 coefficient of the first vector
  654.      * @param v1 first vector
  655.      * @param a2 coefficient of the second vector
  656.      * @param v2 second vector
  657.      * @param a3 coefficient of the third vector
  658.      * @param v3 third vector
  659.      * @param row Jacobian matrix row
  660.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  661.      */
  662.     protected void fillHalfRow(final T a1, final FieldVector3D<T> v1,
  663.                                final T a2, final FieldVector3D<T> v2,
  664.                                final T a3, final FieldVector3D<T> v3,
  665.                                       final T[] row, final int j) {
  666.         row[j]     = a1.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX());
  667.         row[j + 1] = a1.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY());
  668.         row[j + 2] = a1.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ());



  672. //        row[j]     = a1.multiply(v1.getX()).add(a2.multiply(v2.getX())).add(a3.multiply(v3.getX()));
  673. //        row[j + 1] = a1.multiply(v1.getY()).add(a2.multiply(v2.getY())).add(a3.multiply(v3.getY()));
  674. //        row[j + 2] = a1.multiply(v1.getZ()).add(a2.multiply(v2.getZ())).add(a3.multiply(v3.getZ()));
  675.     }

  677.     /** Fill a Jacobian half row with a linear combination of vectors.
  678.      * @param a1 coefficient of the first vector
  679.      * @param v1 first vector
  680.      * @param a2 coefficient of the second vector
  681.      * @param v2 second vector
  682.      * @param a3 coefficient of the third vector
  683.      * @param v3 third vector
  684.      * @param a4 coefficient of the fourth vector
  685.      * @param v4 fourth vector
  686.      * @param row Jacobian matrix row
  687.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  688.      */
  689.     protected void fillHalfRow(final T a1, final FieldVector3D<T> v1, final T a2, final FieldVector3D<T> v2,
  690.                                       final T a3, final FieldVector3D<T> v3, final T a4, final FieldVector3D<T> v4,
  691.                                       final T[] row, final int j) {
  692.         row[j]     = a1.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX(), a4, v4.getX());
  693.         row[j + 1] = a1.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY(), a4, v4.getY());
  694.         row[j + 2] = a1.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ(), a4, v4.getZ());
  695. //        row[j]     = a1.multiply(v1.getX()).add(a2.multiply(v2.getX())).add(a3.multiply(v3.getX())).add(a4.multiply(v4.getX()));
  696. //        row[j + 1] = a1.multiply(v1.getY()).add(a2.multiply(v2.getY())).add(a3.multiply(v3.getY())).add(a4.multiply(v4.getY()));
  697. //        row[j + 2] = a1.multiply(v1.getZ()).add(a2.multiply(v2.getZ())).add(a3.multiply(v3.getZ())).add(a4.multiply(v4.getZ()));
  698.     }

  700.     /** Fill a Jacobian half row with a linear combination of vectors.
  701.      * @param a1 coefficient of the first vector
  702.      * @param v1 first vector
  703.      * @param a2 coefficient of the second vector
  704.      * @param v2 second vector
  705.      * @param a3 coefficient of the third vector
  706.      * @param v3 third vector
  707.      * @param a4 coefficient of the fourth vector
  708.      * @param v4 fourth vector
  709.      * @param a5 coefficient of the fifth vector
  710.      * @param v5 fifth vector
  711.      * @param row Jacobian matrix row
  712.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  713.      */
  714.     protected void fillHalfRow(final T a1, final FieldVector3D<T> v1, final T a2, final FieldVector3D<T> v2,
  715.                                       final T a3, final FieldVector3D<T> v3, final T a4, final FieldVector3D<T> v4,
  716.                                       final T a5, final FieldVector3D<T> v5,
  717.                                       final T[] row, final int j) {
  718.         row[j]     = a1.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX(), a4, v4.getX()).add(a5.multiply(v5.getX()));
  719.         row[j + 1] = a1.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY(), a4, v4.getY()).add(a5.multiply(v5.getY()));
  720.         row[j + 2] = a1.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ(), a4, v4.getZ()).add(a5.multiply(v5.getZ()));

  722. //        row[j]     = a1.multiply(v1.getX()).add(a2.multiply(v2.getX())).add(a3.multiply(v3.getX())).add(a4.multiply(v4.getX())).add(a5.multiply(v5.getX()));
  723. //        row[j + 1] = a1.multiply(v1.getY()).add(a2.multiply(v2.getY())).add(a3.multiply(v3.getY())).add(a4.multiply(v4.getY())).add(a5.multiply(v5.getY()));
  724. //        row[j + 2] = a1.multiply(v1.getZ()).add(a2.multiply(v2.getZ())).add(a3.multiply(v3.getZ())).add(a4.multiply(v4.getZ())).add(a5.multiply(v5.getZ()));
  725.     }

  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 a3 coefficient of the third vector
  733.      * @param v3 third vector
  734.      * @param a4 coefficient of the fourth vector
  735.      * @param v4 fourth vector
  736.      * @param a5 coefficient of the fifth vector
  737.      * @param v5 fifth vector
  738.      * @param a6 coefficient of the sixth vector
  739.      * @param v6 sixth vector
  740.      * @param row Jacobian matrix row
  741.      * @param j index of the first element to set (row[j], row[j+1] and row[j+2] will all be set)
  742.      */
  743.     protected void fillHalfRow(final T a1, final FieldVector3D<T> v1, final T a2, final FieldVector3D<T> v2,
  744.                                       final T a3, final FieldVector3D<T> v3, final T a4, final FieldVector3D<T> v4,
  745.                                       final T a5, final FieldVector3D<T> v5, final T a6, final FieldVector3D<T> v6,
  746.                                       final T[] row, final int j) {
  747.         row[j]     = a1.linearCombination(a1, v1.getX(), a2, v2.getX(), a3, v3.getX(), a4, v4.getX()).add(a1.linearCombination(a5, v5.getX(), a6, v6.getX()));
  748.         row[j + 1] = a1.linearCombination(a1, v1.getY(), a2, v2.getY(), a3, v3.getY(), a4, v4.getY()).add(a1.linearCombination(a5, v5.getY(), a6, v6.getY()));
  749.         row[j + 2] = a1.linearCombination(a1, v1.getZ(), a2, v2.getZ(), a3, v3.getZ(), a4, v4.getZ()).add(a1.linearCombination(a5, v5.getZ(), a6, v6.getZ()));


  752. //        row[j]     = a1.multiply(v1.getX()).add(a2.multiply(v2.getX())).add(a3.multiply(v3.getX())).add(a4.multiply(v4.getX())).add(a5.multiply(v5.getX())).add(a6.multiply(v6.getX()));
  753. //        row[j + 1] = a1.multiply(v1.getY()).add(a2.multiply(v2.getY())).add(a3.multiply(v3.getY())).add(a4.multiply(v4.getY())).add(a5.multiply(v5.getY())).add(a6.multiply(v6.getY()));
  754. //        row[j + 2] = a1.multiply(v1.getZ()).add(a2.multiply(v2.getZ())).add(a3.multiply(v3.getZ())).add(a4.multiply(v4.getZ())).add(a5.multiply(v5.getZ())).add(a6.multiply(v6.getZ()));
  755.     }


  758. }
Created with JaCoCo 0.7.9.201702052155