FieldTransform.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.frames;

  18. import java.util.ArrayList;
  19. import java.util.Arrays;
  20. import java.util.Collection;
  21. import java.util.List;
  22. import java.util.stream.Stream;

  23. import org.hipparchus.Field;
  24. import org.hipparchus.RealFieldElement;
  25. import org.hipparchus.geometry.euclidean.threed.FieldLine;
  26. import org.hipparchus.geometry.euclidean.threed.FieldRotation;
  27. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  28. import org.hipparchus.geometry.euclidean.threed.Line;
  29. import org.hipparchus.geometry.euclidean.threed.RotationConvention;
  30. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  31. import org.orekit.errors.OrekitException;
  32. import org.orekit.time.AbsoluteDate;
  33. import org.orekit.time.FieldAbsoluteDate;
  34. import org.orekit.time.TimeShiftable;
  35. import org.orekit.time.TimeStamped;
  36. import org.orekit.utils.AngularDerivativesFilter;
  37. import org.orekit.utils.CartesianDerivativesFilter;
  38. import org.orekit.utils.FieldAngularCoordinates;
  39. import org.orekit.utils.FieldPVCoordinates;
  40. import org.orekit.utils.PVCoordinates;
  41. import org.orekit.utils.TimeStampedFieldAngularCoordinates;
  42. import org.orekit.utils.TimeStampedFieldPVCoordinates;
  43. import org.orekit.utils.TimeStampedPVCoordinates;


  44. /** Transformation class in three dimensional space.
  45.  *
  46.  * <p>This class represents the transformation engine between {@link Frame frames}.
  47.  * It is used both to define the relationship between each frame and its
  48.  * parent frame and to gather all individual transforms into one
  49.  * operation when converting between frames far away from each other.</p>
  50.  * <p>The convention used in OREKIT is vectorial transformation. It means
  51.  * that a transformation is defined as a transform to apply to the
  52.  * coordinates of a vector expressed in the old frame to obtain the
  53.  * same vector expressed in the new frame.
  54.  *
  55.  * <p>Instances of this class are guaranteed to be immutable.</p>
  56.  *
  57.  * <h1> Examples </h1>
  58.  *
  59.  * <h2> Example of translation from R<sub>A</sub> to R<sub>B</sub> </h2>
  60.  *
  61.  * <p> We want to transform the {@link FieldPVCoordinates} PV<sub>A</sub> to
  62.  * PV<sub>B</sub> with :
  63.  * <p> PV<sub>A</sub> = ({1, 0, 0}, {2, 0, 0}, {3, 0, 0}); <br>
  64.  *     PV<sub>B</sub> = ({0, 0, 0}, {0, 0, 0}, {0, 0, 0});
  65.  *
  66.  * <p> The transform to apply then is defined as follows :
  67.  *
  68.  * <pre><code>
  69.  * Vector3D translation  = new Vector3D(-1, 0, 0);
  70.  * Vector3D velocity     = new Vector3D(-2, 0, 0);
  71.  * Vector3D acceleration = new Vector3D(-3, 0, 0);
  72.  *
  73.  * Transform R1toR2 = new Transform(date, translation, velocity, acceleration);
  74.  *
  75.  * PVB = R1toR2.transformPVCoordinates<T>(PVA);
  76.  * </code></pre>
  77.  *
  78.  * <h2> Example of rotation from R<sub>A</sub> to R<sub>B</sub> </h2>
  79.  * <p> We want to transform the {@link FieldPVCoordinates} PV<sub>A</sub> to
  80.  * PV<sub>B</sub> with
  81.  *
  82.  * <p> PV<sub>A</sub> = ({1, 0, 0}, { 1, 0, 0}); <br>
  83.  *     PV<sub>B</sub> = ({0, 1, 0}, {-2, 1, 0});
  84.  *
  85.  * <p> The transform to apply then is defined as follows :
  86.  *
  87.  * <pre><code>
  88.  * Rotation rotation = new Rotation(Vector3D.PLUS_K, FastMath.PI / 2);
  89.  * Vector3D rotationRate = new Vector3D(0, 0, -2);
  90.  *
  91.  * Transform R1toR2 = new Transform(rotation, rotationRate);
  92.  *
  93.  * PVB = R1toR2.transformPVCoordinates<T>(PVA);
  94.  * </code></pre>
  95.  *
  96.  * @author Luc Maisonobe
  97.  * @author Fabien Maussion
  98.  * @param <T> the type of the field elements
  99.  * @since 9.0
  100.  */
  101. public class FieldTransform<T extends RealFieldElement<T>>
  102.     implements TimeStamped, TimeShiftable<FieldTransform<T>> {

  103.     /** Date of the transform. */
  104.     private final FieldAbsoluteDate<T> date;

  105.     /** Date of the transform. */
  106.     private final AbsoluteDate aDate;

  107.     /** Cartesian coordinates of the target frame with respect to the original frame. */
  108.     private final FieldPVCoordinates<T> cartesian;

  109.     /** Angular coordinates of the target frame with respect to the original frame. */
  110.     private final FieldAngularCoordinates<T> angular;

  111.     /** Build a transform from its primitive operations.
  112.      * @param date date of the transform
  113.      * @param aDate date of the transform
  114.      * @param cartesian Cartesian coordinates of the target frame with respect to the original frame
  115.      * @param angular angular coordinates of the target frame with respect to the original frame
  116.      */
  117.     private FieldTransform(final FieldAbsoluteDate<T> date, final AbsoluteDate aDate,
  118.                            final FieldPVCoordinates<T> cartesian,
  119.                            final FieldAngularCoordinates<T> angular) {
  120.         this.date      = date;
  121.         this.aDate     = aDate;
  122.         this.cartesian = cartesian;
  123.         this.angular   = angular;
  124.     }

  125.     /** Build a transform from a regular transform.
  126.      * @param field field of the elements
  127.      * @param transform regular transform to convert
  128.      */
  129.     public FieldTransform(final Field<T> field, final Transform transform) {
  130.         this(new FieldAbsoluteDate<>(field, transform.getDate()), transform.getDate(),
  131.              new FieldPVCoordinates<>(field, transform.getCartesian()),
  132.              new FieldAngularCoordinates<>(field, transform.getAngular()));
  133.     }

  134.     /** Build a translation transform.
  135.      * @param date date of the transform
  136.      * @param translation translation to apply (i.e. coordinates of
  137.      * the transformed origin, or coordinates of the origin of the
  138.      * old frame in the new frame)
  139.      */
  140.     public FieldTransform(final FieldAbsoluteDate<T> date, final FieldVector3D<T> translation) {
  141.         this(date, date.toAbsoluteDate(),
  142.              new FieldPVCoordinates<>(translation,
  143.                                       FieldVector3D.getZero(date.getField()),
  144.                                       FieldVector3D.getZero(date.getField())),
  145.              FieldAngularCoordinates.getIdentity(date.getField()));
  146.     }

  147.     /** Build a rotation transform.
  148.      * @param date date of the transform
  149.      * @param rotation rotation to apply ( i.e. rotation to apply to the
  150.      * coordinates of a vector expressed in the old frame to obtain the
  151.      * same vector expressed in the new frame )
  152.      */
  153.     public FieldTransform(final FieldAbsoluteDate<T> date, final FieldRotation<T> rotation) {
  154.         this(date, date.toAbsoluteDate(),
  155.              FieldPVCoordinates.getZero(date.getField()),
  156.              new FieldAngularCoordinates<>(rotation,
  157.                                            FieldVector3D.getZero(date.getField()),
  158.                                            FieldVector3D.getZero(date.getField())));
  159.     }

  160.     /** Build a translation transform, with its first time derivative.
  161.      * @param date date of the transform
  162.      * @param translation translation to apply (i.e. coordinates of
  163.      * the transformed origin, or coordinates of the origin of the
  164.      * old frame in the new frame)
  165.      * @param velocity the velocity of the translation (i.e. origin
  166.      * of the old frame velocity in the new frame)
  167.      */
  168.     public FieldTransform(final FieldAbsoluteDate<T> date,
  169.                           final FieldVector3D<T> translation,
  170.                           final FieldVector3D<T> velocity) {
  171.         this(date,
  172.              new FieldPVCoordinates<>(translation,
  173.                                       velocity,
  174.                                       FieldVector3D.getZero(date.getField())));
  175.     }

  176.     /** Build a translation transform, with its first and second time derivatives.
  177.      * @param date date of the transform
  178.      * @param translation translation to apply (i.e. coordinates of
  179.      * the transformed origin, or coordinates of the origin of the
  180.      * old frame in the new frame)
  181.      * @param velocity the velocity of the translation (i.e. origin
  182.      * of the old frame velocity in the new frame)
  183.      * @param acceleration the acceleration of the translation (i.e. origin
  184.      * of the old frame acceleration in the new frame)
  185.      */
  186.     public FieldTransform(final FieldAbsoluteDate<T> date, final FieldVector3D<T> translation,
  187.                           final FieldVector3D<T> velocity, final FieldVector3D<T> acceleration) {
  188.         this(date,
  189.              new FieldPVCoordinates<>(translation, velocity, acceleration));
  190.     }

  191.     /** Build a translation transform, with its first time derivative.
  192.      * @param date date of the transform
  193.      * @param cartesian Cartesian part of the transformation to apply (i.e. coordinates of
  194.      * the transformed origin, or coordinates of the origin of the
  195.      * old frame in the new frame, with their derivatives)
  196.      */
  197.     public FieldTransform(final FieldAbsoluteDate<T> date, final FieldPVCoordinates<T> cartesian) {
  198.         this(date, date.toAbsoluteDate(),
  199.              cartesian,
  200.              FieldAngularCoordinates.getIdentity(date.getField()));
  201.     }

  202.     /** Build a rotation transform.
  203.      * @param date date of the transform
  204.      * @param rotation rotation to apply ( i.e. rotation to apply to the
  205.      * coordinates of a vector expressed in the old frame to obtain the
  206.      * same vector expressed in the new frame )
  207.      * @param rotationRate the axis of the instant rotation
  208.      * expressed in the new frame. (norm representing angular rate)
  209.      */
  210.     public FieldTransform(final FieldAbsoluteDate<T> date,
  211.                           final FieldRotation<T> rotation,
  212.                           final FieldVector3D<T> rotationRate) {
  213.         this(date,
  214.              new FieldAngularCoordinates<>(rotation,
  215.                                            rotationRate,
  216.                                            FieldVector3D.getZero(date.getField())));
  217.     }

  218.     /** Build a rotation transform.
  219.      * @param date date of the transform
  220.      * @param rotation rotation to apply ( i.e. rotation to apply to the
  221.      * coordinates of a vector expressed in the old frame to obtain the
  222.      * same vector expressed in the new frame )
  223.      * @param rotationRate the axis of the instant rotation
  224.      * @param rotationAcceleration the axis of the instant rotation
  225.      * expressed in the new frame. (norm representing angular rate)
  226.      */
  227.     public FieldTransform(final FieldAbsoluteDate<T> date,
  228.                           final FieldRotation<T> rotation,
  229.                           final FieldVector3D<T> rotationRate,
  230.                           final FieldVector3D<T> rotationAcceleration) {
  231.         this(date, new FieldAngularCoordinates<>(rotation, rotationRate, rotationAcceleration));
  232.     }

  233.     /** Build a rotation transform.
  234.      * @param date date of the transform
  235.      * @param angular angular part of the transformation to apply (i.e. rotation to
  236.      * apply to the coordinates of a vector expressed in the old frame to obtain the
  237.      * same vector expressed in the new frame, with its rotation rate)
  238.      */
  239.     public FieldTransform(final FieldAbsoluteDate<T> date, final FieldAngularCoordinates<T> angular) {
  240.         this(date, date.toAbsoluteDate(),
  241.              FieldPVCoordinates.getZero(date.getField()),
  242.              angular);
  243.     }

  244.     /** Build a transform by combining two existing ones.
  245.      * <p>
  246.      * Note that the dates of the two existing transformed are <em>ignored</em>,
  247.      * and the combined transform date is set to the date supplied in this constructor
  248.      * without any attempt to shift the raw transforms. This is a design choice allowing
  249.      * user full control of the combination.
  250.      * </p>
  251.      * @param date date of the transform
  252.      * @param first first transform applied
  253.      * @param second second transform applied
  254.      */
  255.     public FieldTransform(final FieldAbsoluteDate<T> date,
  256.                           final FieldTransform<T> first,
  257.                           final FieldTransform<T> second) {
  258.         this(date, date.toAbsoluteDate(),
  259.              new FieldPVCoordinates<>(compositeTranslation(first, second),
  260.                                       compositeVelocity(first, second),
  261.                                       compositeAcceleration(first, second)),
  262.              new FieldAngularCoordinates<>(compositeRotation(first, second),
  263.                                            compositeRotationRate(first, second),
  264.                                            compositeRotationAcceleration(first, second)));
  265.     }

  266.     /** Get the identity transform.
  267.      * @param field field for the components
  268.      * @param <T> the type of the field elements
  269.      * @return identity transform
  270.      */
  271.     public static <T extends RealFieldElement<T>> FieldTransform<T> getIdentity(final Field<T> field) {
  272.         return new FieldIdentityTransform<>(field);
  273.     }

  274.     /** Compute a composite translation.
  275.      * @param first first applied transform
  276.      * @param second second applied transform
  277.      * @param <T> the type of the field elements
  278.      * @return translation part of the composite transform
  279.      */
  280.     private static <T extends RealFieldElement<T>> FieldVector3D<T> compositeTranslation(final FieldTransform<T> first, final FieldTransform<T> second) {

  281.         final FieldVector3D<T> p1 = first.cartesian.getPosition();
  282.         final FieldRotation<T> r1 = first.angular.getRotation();
  283.         final FieldVector3D<T> p2 = second.cartesian.getPosition();

  284.         return p1.add(r1.applyInverseTo(p2));

  285.     }

  286.     /** Compute a composite velocity.
  287.      * @param first first applied transform
  288.      * @param second second applied transform
  289.      * @param <T> the type of the field elements
  290.      * @return velocity part of the composite transform
  291.      */
  292.     private static <T extends RealFieldElement<T>> FieldVector3D<T> compositeVelocity(final FieldTransform<T> first, final FieldTransform<T> second) {

  293.         final FieldVector3D<T> v1 = first.cartesian.getVelocity();
  294.         final FieldRotation<T> r1 = first.angular.getRotation();
  295.         final FieldVector3D<T> o1 = first.angular.getRotationRate();
  296.         final FieldVector3D<T> p2 = second.cartesian.getPosition();
  297.         final FieldVector3D<T> v2 = second.cartesian.getVelocity();

  298.         final FieldVector3D<T> crossP = FieldVector3D.crossProduct(o1, p2);

  299.         return v1.add(r1.applyInverseTo(v2.add(crossP)));

  300.     }

  301.     /** Compute a composite acceleration.
  302.      * @param first first applied transform
  303.      * @param second second applied transform
  304.      * @param <T> the type of the field elements
  305.      * @return acceleration part of the composite transform
  306.      */
  307.     private static <T extends RealFieldElement<T>> FieldVector3D<T> compositeAcceleration(final FieldTransform<T> first, final FieldTransform<T> second) {

  308.         final FieldVector3D<T> a1    = first.cartesian.getAcceleration();
  309.         final FieldRotation<T> r1    = first.angular.getRotation();
  310.         final FieldVector3D<T> o1    = first.angular.getRotationRate();
  311.         final FieldVector3D<T> oDot1 = first.angular.getRotationAcceleration();
  312.         final FieldVector3D<T> p2    = second.cartesian.getPosition();
  313.         final FieldVector3D<T> v2    = second.cartesian.getVelocity();
  314.         final FieldVector3D<T> a2    = second.cartesian.getAcceleration();

  315.         final FieldVector3D<T> crossCrossP = FieldVector3D.crossProduct(o1,    FieldVector3D.crossProduct(o1, p2));
  316.         final FieldVector3D<T> crossV      = FieldVector3D.crossProduct(o1,    v2);
  317.         final FieldVector3D<T> crossDotP   = FieldVector3D.crossProduct(oDot1, p2);

  318.         return a1.add(r1.applyInverseTo(new FieldVector3D<>(1, a2, 2, crossV, 1, crossCrossP, 1, crossDotP)));

  319.     }

  320.     /** Compute a composite rotation.
  321.      * @param first first applied transform
  322.      * @param second second applied transform
  323.      * @param <T> the type of the field elements
  324.      * @return rotation part of the composite transform
  325.      */
  326.     private static <T extends RealFieldElement<T>> FieldRotation<T> compositeRotation(final FieldTransform<T> first, final FieldTransform<T> second) {

  327.         final FieldRotation<T> r1 = first.angular.getRotation();
  328.         final FieldRotation<T> r2 = second.angular.getRotation();

  329.         return r1.compose(r2, RotationConvention.FRAME_TRANSFORM);

  330.     }

  331.     /** Compute a composite rotation rate.
  332.      * @param first first applied transform
  333.      * @param second second applied transform
  334.      * @param <T> the type of the field elements
  335.      * @return rotation rate part of the composite transform
  336.      */
  337.     private static <T extends RealFieldElement<T>> FieldVector3D<T> compositeRotationRate(final FieldTransform<T> first, final FieldTransform<T> second) {

  338.         final FieldVector3D<T> o1 = first.angular.getRotationRate();
  339.         final FieldRotation<T> r2 = second.angular.getRotation();
  340.         final FieldVector3D<T> o2 = second.angular.getRotationRate();

  341.         return o2.add(r2.applyTo(o1));

  342.     }

  343.     /** Compute a composite rotation acceleration.
  344.      * @param first first applied transform
  345.      * @param second second applied transform
  346.      * @param <T> the type of the field elements
  347.      * @return rotation acceleration part of the composite transform
  348.      */
  349.     private static <T extends RealFieldElement<T>> FieldVector3D<T> compositeRotationAcceleration(final FieldTransform<T> first, final FieldTransform<T> second) {

  350.         final FieldVector3D<T> o1    = first.angular.getRotationRate();
  351.         final FieldVector3D<T> oDot1 = first.angular.getRotationAcceleration();
  352.         final FieldRotation<T> r2    = second.angular.getRotation();
  353.         final FieldVector3D<T> o2    = second.angular.getRotationRate();
  354.         final FieldVector3D<T> oDot2 = second.angular.getRotationAcceleration();

  355.         return new FieldVector3D<>( 1, oDot2,
  356.                                     1, r2.applyTo(oDot1),
  357.                                    -1, FieldVector3D.crossProduct(o2, r2.applyTo(o1)));

  358.     }

  359.     /** {@inheritDoc} */
  360.     public AbsoluteDate getDate() {
  361.         return aDate;
  362.     }

  363.     /** Get the date.
  364.      * @return date attached to the object
  365.      */
  366.     public FieldAbsoluteDate<T> getFieldDate() {
  367.         return date;
  368.     }

  369.     /** {@inheritDoc} */
  370.     public FieldTransform<T> shiftedBy(final double dt) {
  371.         return new FieldTransform<>(date.shiftedBy(dt), aDate.shiftedBy(dt),
  372.                                     cartesian.shiftedBy(dt), angular.shiftedBy(dt));
  373.     }

  374.     /** Get a time-shifted instance.
  375.      * @param dt time shift in seconds
  376.      * @return a new instance, shifted with respect to instance (which is not changed)
  377.      */
  378.     FieldTransform<T> shiftedBy(final T dt) {
  379.         return new FieldTransform<>(date.shiftedBy(dt), aDate.shiftedBy(dt.getReal()),
  380.                                     cartesian.shiftedBy(dt), angular.shiftedBy(dt));
  381.     }

  382.     /** Interpolate a transform from a sample set of existing transforms.
  383.      * <p>
  384.      * Calling this method is equivalent to call {@link #interpolate(FieldAbsoluteDate,
  385.      * CartesianDerivativesFilter, AngularDerivativesFilter, Collection)} with {@code cFilter}
  386.      * set to {@link CartesianDerivativesFilter#USE_PVA} and {@code aFilter} set to
  387.      * {@link AngularDerivativesFilter#USE_RRA}
  388.      * set to true.
  389.      * </p>
  390.      * @param interpolationDate interpolation date
  391.      * @param sample sample points on which interpolation should be done
  392.      * @param <T> the type of the field elements
  393.      * @return a new instance, interpolated at specified date
  394.      * @exception OrekitException if the number of point is too small for interpolating
  395.      */
  396.     public static <T extends RealFieldElement<T>> FieldTransform<T> interpolate(final FieldAbsoluteDate<T> interpolationDate,
  397.                                                                                 final Collection<FieldTransform<T>> sample)
  398.         throws OrekitException {
  399.         return interpolate(interpolationDate,
  400.                            CartesianDerivativesFilter.USE_PVA, AngularDerivativesFilter.USE_RRA,
  401.                            sample);
  402.     }

  403.     /** Interpolate a transform from a sample set of existing transforms.
  404.      * <p>
  405.      * Note that even if first time derivatives (velocities and rotation rates)
  406.      * from sample can be ignored, the interpolated instance always includes
  407.      * interpolated derivatives. This feature can be used explicitly to
  408.      * compute these derivatives when it would be too complex to compute them
  409.      * from an analytical formula: just compute a few sample points from the
  410.      * explicit formula and set the derivatives to zero in these sample points,
  411.      * then use interpolation to add derivatives consistent with the positions
  412.      * and rotations.
  413.      * </p>
  414.      * <p>
  415.      * As this implementation of interpolation is polynomial, it should be used only
  416.      * with small samples (about 10-20 points) in order to avoid <a
  417.      * href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
  418.      * and numerical problems (including NaN appearing).
  419.      * </p>
  420.      * @param date interpolation date
  421.      * @param cFilter filter for derivatives from the sample to use in interpolation
  422.      * @param aFilter filter for derivatives from the sample to use in interpolation
  423.      * @param sample sample points on which interpolation should be done
  424.      * @return a new instance, interpolated at specified date
  425.      * @exception OrekitException if the number of point is too small for interpolating
  426.      * @param <T> the type of the field elements
  427.      */
  428.     public static <T extends RealFieldElement<T>> FieldTransform<T> interpolate(final FieldAbsoluteDate<T> date,
  429.                                                                                 final CartesianDerivativesFilter cFilter,
  430.                                                                                 final AngularDerivativesFilter aFilter,
  431.                                                                                 final Collection<FieldTransform<T>> sample)
  432.         throws OrekitException {
  433.         return interpolate(date, cFilter, aFilter, sample.stream());
  434.     }

  435.     /** Interpolate a transform from a sample set of existing transforms.
  436.      * <p>
  437.      * Note that even if first time derivatives (velocities and rotation rates)
  438.      * from sample can be ignored, the interpolated instance always includes
  439.      * interpolated derivatives. This feature can be used explicitly to
  440.      * compute these derivatives when it would be too complex to compute them
  441.      * from an analytical formula: just compute a few sample points from the
  442.      * explicit formula and set the derivatives to zero in these sample points,
  443.      * then use interpolation to add derivatives consistent with the positions
  444.      * and rotations.
  445.      * </p>
  446.      * <p>
  447.      * As this implementation of interpolation is polynomial, it should be used only
  448.      * with small samples (about 10-20 points) in order to avoid <a
  449.      * href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
  450.      * and numerical problems (including NaN appearing).
  451.      * </p>
  452.      * @param date interpolation date
  453.      * @param cFilter filter for derivatives from the sample to use in interpolation
  454.      * @param aFilter filter for derivatives from the sample to use in interpolation
  455.      * @param sample sample points on which interpolation should be done
  456.      * @return a new instance, interpolated at specified date
  457.      * @exception OrekitException if the number of point is too small for interpolating
  458.      * @param <T> the type of the field elements
  459.      */
  460.     public static <T extends RealFieldElement<T>> FieldTransform<T> interpolate(final FieldAbsoluteDate<T> date,
  461.                                                                                 final CartesianDerivativesFilter cFilter,
  462.                                                                                 final AngularDerivativesFilter aFilter,
  463.                                                                                 final Stream<FieldTransform<T>> sample)
  464.         throws OrekitException {
  465.         final List<TimeStampedFieldPVCoordinates<T>>      datedPV = new ArrayList<>();
  466.         final List<TimeStampedFieldAngularCoordinates<T>> datedAC = new ArrayList<>();
  467.         sample.forEach(t -> {
  468.             datedPV.add(new TimeStampedFieldPVCoordinates<>(t.getDate(), t.getTranslation(), t.getVelocity(), t.getAcceleration()));
  469.             datedAC.add(new TimeStampedFieldAngularCoordinates<>(t.getDate(), t.getRotation(), t.getRotationRate(), t.getRotationAcceleration()));
  470.         });
  471.         final TimeStampedFieldPVCoordinates<T>      interpolatedPV = TimeStampedFieldPVCoordinates.interpolate(date, cFilter, datedPV);
  472.         final TimeStampedFieldAngularCoordinates<T> interpolatedAC = TimeStampedFieldAngularCoordinates.interpolate(date, aFilter, datedAC);
  473.         return new FieldTransform<>(date, date.toAbsoluteDate(), interpolatedPV, interpolatedAC);
  474.     }

  475.     /** Get the inverse transform of the instance.
  476.      * @return inverse transform of the instance
  477.      */
  478.     public FieldTransform<T> getInverse() {

  479.         final FieldRotation<T> r    = angular.getRotation();
  480.         final FieldVector3D<T> o    = angular.getRotationRate();
  481.         final FieldVector3D<T> oDot = angular.getRotationAcceleration();
  482.         final FieldVector3D<T> rp   = r.applyTo(cartesian.getPosition());
  483.         final FieldVector3D<T> rv   = r.applyTo(cartesian.getVelocity());
  484.         final FieldVector3D<T> ra   = r.applyTo(cartesian.getAcceleration());

  485.         final FieldVector3D<T> pInv        = rp.negate();
  486.         final FieldVector3D<T> crossP      = FieldVector3D.crossProduct(o, rp);
  487.         final FieldVector3D<T> vInv        = crossP.subtract(rv);
  488.         final FieldVector3D<T> crossV      = FieldVector3D.crossProduct(o, rv);
  489.         final FieldVector3D<T> crossDotP   = FieldVector3D.crossProduct(oDot, rp);
  490.         final FieldVector3D<T> crossCrossP = FieldVector3D.crossProduct(o, crossP);
  491.         final FieldVector3D<T> aInv        = new FieldVector3D<>(-1, ra,
  492.                                                                   2, crossV,
  493.                                                                   1, crossDotP,
  494.                                                                  -1, crossCrossP);

  495.         return new FieldTransform<>(date, aDate, new FieldPVCoordinates<>(pInv, vInv, aInv), angular.revert());

  496.     }

  497.     /** Get a frozen transform.
  498.      * <p>
  499.      * This method creates a copy of the instance but frozen in time,
  500.      * i.e. with velocity, acceleration and rotation rate forced to zero.
  501.      * </p>
  502.      * @return a new transform, without any time-dependent parts
  503.      */
  504.     public FieldTransform<T> freeze() {
  505.         return new FieldTransform<>(date, aDate,
  506.                                     new FieldPVCoordinates<>(cartesian.getPosition(),
  507.                                                              FieldVector3D.getZero(date.getField()),
  508.                                                              FieldVector3D.getZero(date.getField())),
  509.                                     new FieldAngularCoordinates<>(angular.getRotation(),
  510.                                                                   FieldVector3D.getZero(date.getField()),
  511.                                                                   FieldVector3D.getZero(date.getField())));
  512.     }

  513.     /** Transform a position vector (including translation effects).
  514.      * @param position vector to transform
  515.      * @return transformed position
  516.      */
  517.     public FieldVector3D<T> transformPosition(final Vector3D position) {
  518.         return angular.getRotation().applyTo(cartesian.getPosition().add(position));
  519.     }

  520.     /** Transform a position vector (including translation effects).
  521.      * @param position vector to transform
  522.      * @return transformed position
  523.      */
  524.     public FieldVector3D<T> transformPosition(final FieldVector3D<T> position) {
  525.         return angular.getRotation().applyTo(position.add(cartesian.getPosition()));
  526.     }

  527.     /** Transform a vector (ignoring translation effects).
  528.      * @param vector vector to transform
  529.      * @return transformed vector
  530.      */
  531.     public FieldVector3D<T> transformVector(final Vector3D vector) {
  532.         return angular.getRotation().applyTo(vector);
  533.     }

  534.     /** Transform a vector (ignoring translation effects).
  535.      * @param vector vector to transform
  536.      * @return transformed vector
  537.      */
  538.     public FieldVector3D<T> transformVector(final FieldVector3D<T> vector) {
  539.         return angular.getRotation().applyTo(vector);
  540.     }

  541.     /** Transform a line.
  542.      * @param line to transform
  543.      * @return transformed line
  544.      */
  545.     public FieldLine<T> transformLine(final Line line) {
  546.         final FieldVector3D<T> transformedP0 = transformPosition(line.getOrigin());
  547.         final FieldVector3D<T> transformedP1 = transformPosition(line.pointAt(1.0e6));
  548.         return new FieldLine<>(transformedP0, transformedP1, 1.0e-10);
  549.     }

  550.     /** Transform a line.
  551.      * @param line to transform
  552.      * @return transformed line
  553.      */
  554.     public FieldLine<T> transformLine(final FieldLine<T> line) {
  555.         final FieldVector3D<T> transformedP0 = transformPosition(line.getOrigin());
  556.         final FieldVector3D<T> transformedP1 = transformPosition(line.pointAt(1.0e6));
  557.         return new FieldLine<>(transformedP0, transformedP1, 1.0e-10);
  558.     }

  559.     /** Transform {@link TimeStampedPVCoordinates} including kinematic effects.
  560.      * <p>
  561.      * In order to allow the user more flexibility, this method does <em>not</em> check for
  562.      * consistency between the transform {@link #getDate() date} and the time-stamped
  563.      * position-velocity {@link TimeStampedPVCoordinates#getDate() date}. The returned
  564.      * value will always have the same {@link TimeStampedPVCoordinates#getDate() date} as
  565.      * the input argument, regardless of the instance {@link #getDate() date}.
  566.      * </p>
  567.      * @param pv time-stamped  position-velocity to transform.
  568.      * @return transformed time-stamped position-velocity
  569.      */
  570.     public FieldPVCoordinates<T> transformPVCoordinates(final PVCoordinates pv) {
  571.         return angular.applyTo(new FieldPVCoordinates<>(cartesian.getPosition().add(pv.getPosition()),
  572.                                                         cartesian.getVelocity().add(pv.getVelocity()),
  573.                                                         cartesian.getAcceleration().add(pv.getAcceleration())));
  574.     }

  575.     /** Transform {@link TimeStampedPVCoordinates} including kinematic effects.
  576.      * <p>
  577.      * In order to allow the user more flexibility, this method does <em>not</em> check for
  578.      * consistency between the transform {@link #getDate() date} and the time-stamped
  579.      * position-velocity {@link TimeStampedPVCoordinates#getDate() date}. The returned
  580.      * value will always have the same {@link TimeStampedPVCoordinates#getDate() date} as
  581.      * the input argument, regardless of the instance {@link #getDate() date}.
  582.      * </p>
  583.      * @param pv time-stamped  position-velocity to transform.
  584.      * @return transformed time-stamped position-velocity
  585.      */
  586.     public TimeStampedFieldPVCoordinates<T> transformPVCoordinates(final TimeStampedPVCoordinates pv) {
  587.         return angular.applyTo(new TimeStampedFieldPVCoordinates<>(pv.getDate(),
  588.                                                                    cartesian.getPosition().add(pv.getPosition()),
  589.                                                                    cartesian.getVelocity().add(pv.getVelocity()),
  590.                                                                    cartesian.getAcceleration().add(pv.getAcceleration())));
  591.     }

  592.     /** Transform {@link TimeStampedFieldPVCoordinates} including kinematic effects.
  593.      * <p>
  594.      * BEWARE! This method does explicit computation of velocity and acceleration by combining
  595.      * the transform velocity, acceleration, rotation rate and rotation acceleration with the
  596.      * velocity and acceleration from the argument. This implies that this method should
  597.      * <em>not</em> be used when derivatives are contained in the {@link RealFieldElement field
  598.      * elements} (typically when using {@link org.hipparchus.analysis.differentiation.DerivativeStructure
  599.      * DerivativeStructure} elements where time is one of the differentiation parameter), otherwise
  600.      * the time derivatives would be computed twice, once explicitly in this method and once implicitly
  601.      * in the field operations. If time derivatives are contained in the field elements themselves,
  602.      * then rather than this method the {@link #transformPosition(FieldVector3D) transformPosition}
  603.      * method should be used, so the derivatives are computed once, as part of the field. This
  604.      * method is rather expected to be used when the field elements are {@link
  605.      * org.hipparchus.analysis.differentiation.DerivativeStructure DerivativeStructure} instances
  606.      * where the differentiation parameters are not time (they can typically be initial state
  607.      * for computing state transition matrices or force models parameters, or ground stations
  608.      * positions, ...).
  609.      * </p>
  610.      * <p>
  611.      * In order to allow the user more flexibility, this method does <em>not</em> check for
  612.      * consistency between the transform {@link #getDate() date} and the time-stamped
  613.      * position-velocity {@link TimeStampedFieldPVCoordinates#getDate() date}. The returned
  614.      * value will always have the same {@link TimeStampedFieldPVCoordinates#getDate() date} as
  615.      * the input argument, regardless of the instance {@link #getDate() date}.
  616.      * </p>
  617.      * @param pv time-stamped position-velocity to transform.
  618.      * @return transformed time-stamped position-velocity
  619.      */
  620.     public FieldPVCoordinates<T> transformPVCoordinates(final FieldPVCoordinates<T> pv) {
  621.         return angular.applyTo(new FieldPVCoordinates<>(pv.getPosition().add(cartesian.getPosition()),
  622.                                                         pv.getVelocity().add(cartesian.getVelocity()),
  623.                                                         pv.getAcceleration().add(cartesian.getAcceleration())));
  624.     }

  625.     /** Transform {@link TimeStampedFieldPVCoordinates} including kinematic effects.
  626.      * <p>
  627.      * BEWARE! This method does explicit computation of velocity and acceleration by combining
  628.      * the transform velocity, acceleration, rotation rate and rotation acceleration with the
  629.      * velocity and acceleration from the argument. This implies that this method should
  630.      * <em>not</em> be used when derivatives are contained in the {@link RealFieldElement field
  631.      * elements} (typically when using {@link org.hipparchus.analysis.differentiation.DerivativeStructure
  632.      * DerivativeStructure} elements where time is one of the differentiation parameter), otherwise
  633.      * the time derivatives would be computed twice, once explicitly in this method and once implicitly
  634.      * in the field operations. If time derivatives are contained in the field elements themselves,
  635.      * then rather than this method the {@link #transformPosition(FieldVector3D) transformPosition}
  636.      * method should be used, so the derivatives are computed once, as part of the field. This
  637.      * method is rather expected to be used when the field elements are {@link
  638.      * org.hipparchus.analysis.differentiation.DerivativeStructure DerivativeStructure} instances
  639.      * where the differentiation parameters are not time (they can typically be initial state
  640.      * for computing state transition matrices or force models parameters, or ground stations
  641.      * positions, ...).
  642.      * </p>
  643.      * <p>
  644.      * In order to allow the user more flexibility, this method does <em>not</em> check for
  645.      * consistency between the transform {@link #getDate() date} and the time-stamped
  646.      * position-velocity {@link TimeStampedFieldPVCoordinates#getDate() date}. The returned
  647.      * value will always have the same {@link TimeStampedFieldPVCoordinates#getDate() date} as
  648.      * the input argument, regardless of the instance {@link #getDate() date}.
  649.      * </p>
  650.      * @param pv time-stamped position-velocity to transform.
  651.      * @return transformed time-stamped position-velocity
  652.      */
  653.     public TimeStampedFieldPVCoordinates<T> transformPVCoordinates(final TimeStampedFieldPVCoordinates<T> pv) {
  654.         return angular.applyTo(new TimeStampedFieldPVCoordinates<>(pv.getDate(),
  655.                                                                    pv.getPosition().add(cartesian.getPosition()),
  656.                                                                    pv.getVelocity().add(cartesian.getVelocity()),
  657.                                                                    pv.getAcceleration().add(cartesian.getAcceleration())));
  658.     }

  659.     /** Compute the Jacobian of the {@link #transformPVCoordinates(FieldPVCoordinates)}
  660.      * method of the transform.
  661.      * <p>
  662.      * Element {@code jacobian[i][j]} is the derivative of Cartesian coordinate i
  663.      * of the transformed {@link FieldPVCoordinates} with respect to Cartesian coordinate j
  664.      * of the input {@link FieldPVCoordinates} in method {@link #transformPVCoordinates(FieldPVCoordinates)}.
  665.      * </p>
  666.      * <p>
  667.      * This definition implies that if we define position-velocity coordinates
  668.      * <pre>
  669.      * PV₁ = transform.transformPVCoordinates<T>(PV₀), then
  670.      * </pre>
  671.      * <p> their differentials dPV₁ and dPV₀ will obey the following relation
  672.      * where J is the matrix computed by this method:
  673.      * <pre>
  674.      * dPV₁ = J &times; dPV₀
  675.      * </pre>
  676.      *
  677.      * @param selector selector specifying the size of the upper left corner that must be filled
  678.      * (either 3x3 for positions only, 6x6 for positions and velocities, 9x9 for positions,
  679.      * velocities and accelerations)
  680.      * @param jacobian placeholder matrix whose upper-left corner is to be filled with
  681.      * the Jacobian, the rest of the matrix remaining untouched
  682.      */
  683.     public void getJacobian(final CartesianDerivativesFilter selector, final T[][] jacobian) {

  684.         final T zero = date.getField().getZero();

  685.         // elementary matrix for rotation
  686.         final T[][] mData = angular.getRotation().getMatrix();

  687.         // dP1/dP0
  688.         System.arraycopy(mData[0], 0, jacobian[0], 0, 3);
  689.         System.arraycopy(mData[1], 0, jacobian[1], 0, 3);
  690.         System.arraycopy(mData[2], 0, jacobian[2], 0, 3);

  691.         if (selector.getMaxOrder() >= 1) {

  692.             // dP1/dV0
  693.             Arrays.fill(jacobian[0], 3, 6, zero);
  694.             Arrays.fill(jacobian[1], 3, 6, zero);
  695.             Arrays.fill(jacobian[2], 3, 6, zero);

  696.             // dV1/dP0
  697.             final FieldVector3D<T> o = angular.getRotationRate();
  698.             final T ox = o.getX();
  699.             final T oy = o.getY();
  700.             final T oz = o.getZ();
  701.             for (int i = 0; i < 3; ++i) {
  702.                 jacobian[3][i] = oz.multiply(mData[1][i]).subtract(oy.multiply(mData[2][i]));
  703.                 jacobian[4][i] = ox.multiply(mData[2][i]).subtract(oz.multiply(mData[0][i]));
  704.                 jacobian[5][i] = oy.multiply(mData[0][i]).subtract(ox.multiply(mData[1][i]));
  705.             }

  706.             // dV1/dV0
  707.             System.arraycopy(mData[0], 0, jacobian[3], 3, 3);
  708.             System.arraycopy(mData[1], 0, jacobian[4], 3, 3);
  709.             System.arraycopy(mData[2], 0, jacobian[5], 3, 3);

  710.             if (selector.getMaxOrder() >= 2) {

  711.                 // dP1/dA0
  712.                 Arrays.fill(jacobian[0], 6, 9, zero);
  713.                 Arrays.fill(jacobian[1], 6, 9, zero);
  714.                 Arrays.fill(jacobian[2], 6, 9, zero);

  715.                 // dV1/dA0
  716.                 Arrays.fill(jacobian[3], 6, 9, zero);
  717.                 Arrays.fill(jacobian[4], 6, 9, zero);
  718.                 Arrays.fill(jacobian[5], 6, 9, zero);

  719.                 // dA1/dP0
  720.                 final FieldVector3D<T> oDot = angular.getRotationAcceleration();
  721.                 final T oDotx  = oDot.getX();
  722.                 final T oDoty  = oDot.getY();
  723.                 final T oDotz  = oDot.getZ();
  724.                 for (int i = 0; i < 3; ++i) {
  725.                     jacobian[6][i] = oDotz.multiply(mData[1][i]).subtract(oDoty.multiply(mData[2][i])).add(oz.multiply(jacobian[4][i]).subtract(oy.multiply(jacobian[5][i])));
  726.                     jacobian[7][i] = oDotx.multiply(mData[2][i]).subtract(oDotz.multiply(mData[0][i])).add(ox.multiply(jacobian[5][i]).subtract(oz.multiply(jacobian[3][i])));
  727.                     jacobian[8][i] = oDoty.multiply(mData[0][i]).subtract(oDotx.multiply(mData[1][i])).add(oy.multiply(jacobian[3][i]).subtract(ox.multiply(jacobian[4][i])));
  728.                 }

  729.                 // dA1/dV0
  730.                 for (int i = 0; i < 3; ++i) {
  731.                     jacobian[6][i + 3] = oz.multiply(mData[1][i]).subtract(oy.multiply(mData[2][i])).multiply(2);
  732.                     jacobian[7][i + 3] = ox.multiply(mData[2][i]).subtract(oz.multiply(mData[0][i])).multiply(2);
  733.                     jacobian[8][i + 3] = oy.multiply(mData[0][i]).subtract(ox.multiply(mData[1][i])).multiply(2);
  734.                 }

  735.                 // dA1/dA0
  736.                 System.arraycopy(mData[0], 0, jacobian[6], 6, 3);
  737.                 System.arraycopy(mData[1], 0, jacobian[7], 6, 3);
  738.                 System.arraycopy(mData[2], 0, jacobian[8], 6, 3);

  739.             }

  740.         }

  741.     }

  742.     /** Get the underlying elementary Cartesian part.
  743.      * <p>A transform can be uniquely represented as an elementary
  744.      * translation followed by an elementary rotation. This method
  745.      * returns this unique elementary translation with its derivative.</p>
  746.      * @return underlying elementary Cartesian part
  747.      * @see #getTranslation()
  748.      * @see #getVelocity()
  749.      */
  750.     public FieldPVCoordinates<T> getCartesian() {
  751.         return cartesian;
  752.     }

  753.     /** Get the underlying elementary translation.
  754.      * <p>A transform can be uniquely represented as an elementary
  755.      * translation followed by an elementary rotation. This method
  756.      * returns this unique elementary translation.</p>
  757.      * @return underlying elementary translation
  758.      * @see #getCartesian()
  759.      * @see #getVelocity()
  760.      * @see #getAcceleration()
  761.      */
  762.     public FieldVector3D<T> getTranslation() {
  763.         return cartesian.getPosition();
  764.     }

  765.     /** Get the first time derivative of the translation.
  766.      * @return first time derivative of the translation
  767.      * @see #getCartesian()
  768.      * @see #getTranslation()
  769.      * @see #getAcceleration()
  770.      */
  771.     public FieldVector3D<T> getVelocity() {
  772.         return cartesian.getVelocity();
  773.     }

  774.     /** Get the second time derivative of the translation.
  775.      * @return second time derivative of the translation
  776.      * @see #getCartesian()
  777.      * @see #getTranslation()
  778.      * @see #getVelocity()
  779.      */
  780.     public FieldVector3D<T> getAcceleration() {
  781.         return cartesian.getAcceleration();
  782.     }

  783.     /** Get the underlying elementary angular part.
  784.      * <p>A transform can be uniquely represented as an elementary
  785.      * translation followed by an elementary rotation. This method
  786.      * returns this unique elementary rotation with its derivative.</p>
  787.      * @return underlying elementary angular part
  788.      * @see #getRotation()
  789.      * @see #getRotationRate()
  790.      * @see #getRotationAcceleration()
  791.      */
  792.     public FieldAngularCoordinates<T> getAngular() {
  793.         return angular;
  794.     }

  795.     /** Get the underlying elementary rotation.
  796.      * <p>A transform can be uniquely represented as an elementary
  797.      * translation followed by an elementary rotation. This method
  798.      * returns this unique elementary rotation.</p>
  799.      * @return underlying elementary rotation
  800.      * @see #getAngular()
  801.      * @see #getRotationRate()
  802.      * @see #getRotationAcceleration()
  803.      */
  804.     public FieldRotation<T> getRotation() {
  805.         return angular.getRotation();
  806.     }

  807.     /** Get the first time derivative of the rotation.
  808.      * <p>The norm represents the angular rate.</p>
  809.      * @return First time derivative of the rotation
  810.      * @see #getAngular()
  811.      * @see #getRotation()
  812.      * @see #getRotationAcceleration()
  813.      */
  814.     public FieldVector3D<T> getRotationRate() {
  815.         return angular.getRotationRate();
  816.     }

  817.     /** Get the second time derivative of the rotation.
  818.      * @return Second time derivative of the rotation
  819.      * @see #getAngular()
  820.      * @see #getRotation()
  821.      * @see #getRotationRate()
  822.      */
  823.     public FieldVector3D<T> getRotationAcceleration() {
  824.         return angular.getRotationAcceleration();
  825.     }

  826.     /** Specialized class for identity transform. */
  827.     private static class FieldIdentityTransform<T extends RealFieldElement<T>> extends FieldTransform<T> {

  828.         /** Simple constructor.
  829.          * @param field field for the components
  830.          */
  831.         FieldIdentityTransform(final Field<T> field) {
  832.             super(FieldAbsoluteDate.getJ2000Epoch(field), AbsoluteDate.J2000_EPOCH,
  833.                   FieldPVCoordinates.getZero(field),
  834.                   FieldAngularCoordinates.getIdentity(field));
  835.         }

  836.         /** {@inheritDoc} */
  837.         @Override
  838.         public FieldTransform<T> shiftedBy(final double dt) {
  839.             return this;
  840.         }

  841.         /** {@inheritDoc} */
  842.         @Override
  843.         public FieldTransform<T> getInverse() {
  844.             return this;
  845.         }

  846.         /** {@inheritDoc} */
  847.         @Override
  848.         public FieldVector3D<T> transformPosition(final FieldVector3D<T> position) {
  849.             return position;
  850.         }

  851.         /** {@inheritDoc} */
  852.         @Override
  853.         public FieldVector3D<T> transformVector(final FieldVector3D<T> vector) {
  854.             return vector;
  855.         }

  856.         /** {@inheritDoc} */
  857.         @Override
  858.         public FieldLine<T> transformLine(final FieldLine<T> line) {
  859.             return line;
  860.         }

  861.         /** {@inheritDoc} */
  862.         @Override
  863.         public FieldPVCoordinates<T> transformPVCoordinates(final FieldPVCoordinates<T> pv) {
  864.             return pv;
  865.         }

  866.         /** {@inheritDoc} */
  867.         @Override
  868.         public void getJacobian(final CartesianDerivativesFilter selector, final T[][] jacobian) {
  869.             final int n = 3 * (selector.getMaxOrder() + 1);
  870.             for (int i = 0; i < n; ++i) {
  871.                 Arrays.fill(jacobian[i], 0, n, getFieldDate().getField().getZero());
  872.                 jacobian[i][i] = getFieldDate().getField().getOne();
  873.             }
  874.         }

  875.     }

  876. }