AbsolutePVCoordinatesHermiteInterpolator.java

  1. /* Copyright 2002-2023 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.utils;

  18. import org.hipparchus.analysis.interpolation.HermiteInterpolator;
  19. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  20. import org.orekit.errors.OrekitInternalError;
  21. import org.orekit.frames.Frame;
  22. import org.orekit.time.AbsoluteDate;
  23. import org.orekit.time.AbstractTimeInterpolator;

  24. import java.util.List;

  25. /**
  26.  * Class using a Hermite interpolator to interpolate absolute position-velocity-acceleration coordinates.
  27.  * <p>
  28.  * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation points
  29.  * (about 10-20 points) in order to avoid <a href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
  30.  * and numerical problems (including NaN appearing).
  31.  *
  32.  * @author Luc Maisonobe
  33.  * @author Vincent Cucchietti
  34.  * @see HermiteInterpolator
  35.  * @see AbsolutePVCoordinates
  36.  */
  37. public class AbsolutePVCoordinatesHermiteInterpolator extends AbstractTimeInterpolator<AbsolutePVCoordinates> {

  38.     /** Filter for derivatives from the sample to use in interpolation. */
  39.     private final CartesianDerivativesFilter filter;

  40.     /** Output frame for the interpolated instance. */
  41.     private final Frame outputFrame;

  42.     /**
  43.      * Constructor with :
  44.      * <ul>
  45.      *     <li>Default number of interpolation points of {@code DEFAULT_INTERPOLATION_POINTS}</li>
  46.      *     <li>Default extrapolation threshold value ({@code DEFAULT_EXTRAPOLATION_THRESHOLD_SEC} s)</li>
  47.      *     <li>Use of position and two time derivatives during interpolation</li>
  48.      * </ul>
  49.      * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation
  50.      * points (about 10-20 points) in order to avoid <a href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's
  51.      * phenomenon</a> and numerical problems (including NaN appearing).
  52.      *
  53.      * @param outputFrame frame for the interpolated instance
  54.      */
  55.     public AbsolutePVCoordinatesHermiteInterpolator(final Frame outputFrame) {
  56.         this(DEFAULT_INTERPOLATION_POINTS, DEFAULT_EXTRAPOLATION_THRESHOLD_SEC, outputFrame,
  57.              CartesianDerivativesFilter.USE_PVA);
  58.     }

  59.     /**
  60.      * Constructor with :
  61.      * <ul>
  62.      *     <li>Default extrapolation threshold value ({@code DEFAULT_EXTRAPOLATION_THRESHOLD_SEC} s)</li>
  63.      *     <li>Use of position and two time derivatives during interpolation</li>
  64.      * </ul>
  65.      * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation
  66.      * points (about 10-20 points) in order to avoid <a href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's
  67.      * phenomenon</a> and numerical problems (including NaN appearing).
  68.      *
  69.      * @param interpolationPoints number of interpolation points
  70.      * @param outputFrame frame for the interpolated instance
  71.      */
  72.     public AbsolutePVCoordinatesHermiteInterpolator(final int interpolationPoints, final Frame outputFrame) {
  73.         this(interpolationPoints, DEFAULT_EXTRAPOLATION_THRESHOLD_SEC, outputFrame, CartesianDerivativesFilter.USE_PVA);
  74.     }

  75.     /**
  76.      * Constructor with default extrapolation threshold value ({@code DEFAULT_EXTRAPOLATION_THRESHOLD_SEC} s).
  77.      * <p>
  78.      * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation
  79.      * points (about 10-20 points) in order to avoid <a href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's
  80.      * phenomenon</a> and numerical problems (including NaN appearing).
  81.      *
  82.      * @param interpolationPoints number of interpolation points
  83.      * @param outputFrame frame for the interpolated instance
  84.      * @param filter filter for derivatives from the sample to use in interpolation
  85.      */
  86.     public AbsolutePVCoordinatesHermiteInterpolator(final int interpolationPoints, final Frame outputFrame,
  87.                                                     final CartesianDerivativesFilter filter) {
  88.         this(interpolationPoints, DEFAULT_EXTRAPOLATION_THRESHOLD_SEC, outputFrame, filter);
  89.     }

  90.     /**
  91.      * Constructor.
  92.      * <p>
  93.      * As this implementation of interpolation is polynomial, it should be used only with small number of interpolation
  94.      * points (about 10-20 points) in order to avoid <a href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's
  95.      * phenomenon</a> and numerical problems (including NaN appearing).
  96.      *
  97.      * @param interpolationPoints number of interpolation points
  98.      * @param extrapolationThreshold extrapolation threshold beyond which the propagation will fail
  99.      * @param outputFrame frame for the interpolated instance
  100.      * @param filter filter for derivatives from the sample to use in interpolation
  101.      */
  102.     public AbsolutePVCoordinatesHermiteInterpolator(final int interpolationPoints, final double extrapolationThreshold,
  103.                                                     final Frame outputFrame, final CartesianDerivativesFilter filter) {
  104.         super(interpolationPoints, extrapolationThreshold);
  105.         this.outputFrame = outputFrame;
  106.         this.filter      = filter;
  107.     }

  108.     /** Get the filter for derivatives from the sample to use in interpolation.
  109.      * @return filter for derivatives from the sample to use in interpolation.
  110.      */
  111.     public CartesianDerivativesFilter getFilter() {
  112.         return filter;
  113.     }

  114.     /** Get output frame for the interpolated instance.
  115.      * @return output frame for the interpolated instance
  116.      */
  117.     public Frame getOutputFrame() {
  118.         return outputFrame;
  119.     }

  120.     /**
  121.      * {@inheritDoc}
  122.      * <p>
  123.      * The interpolated instance is created by polynomial Hermite interpolation ensuring velocity remains the exact
  124.      * derivative of position.
  125.      * <p>
  126.      * Note that even if first time derivatives (velocities) from sample can be ignored, the interpolated instance always
  127.      * includes interpolated derivatives. This feature can be used explicitly to compute these derivatives when it would be
  128.      * too complex to compute them from an analytical formula: just compute a few sample points from the explicit formula and
  129.      * set the derivatives to zero in these sample points, then use interpolation to add derivatives consistent with the
  130.      * positions.
  131.      */
  132.     @Override
  133.     protected AbsolutePVCoordinates interpolate(final InterpolationData interpolationData) {

  134.         // Get date
  135.         final AbsoluteDate date = interpolationData.getInterpolationDate();

  136.         // Get sample
  137.         final List<AbsolutePVCoordinates> sample = interpolationData.getNeighborList();

  138.         // Set up an interpolator taking derivatives into account
  139.         final HermiteInterpolator interpolator = new HermiteInterpolator();

  140.         // Add sample points
  141.         switch (filter) {
  142.             case USE_P:
  143.                 // Populate sample with position data, ignoring velocity
  144.                 sample.forEach(pv -> {
  145.                     final Vector3D position = pv.getPosition();
  146.                     interpolator.addSamplePoint(pv.getDate().durationFrom(date),
  147.                                                 position.toArray());
  148.                 });
  149.                 break;
  150.             case USE_PV:
  151.                 // Populate sample with position and velocity data
  152.                 sample.forEach(pv -> {
  153.                     final Vector3D position = pv.getPosition();
  154.                     final Vector3D velocity = pv.getVelocity();
  155.                     interpolator.addSamplePoint(pv.getDate().durationFrom(date),
  156.                                                 position.toArray(), velocity.toArray());
  157.                 });
  158.                 break;
  159.             case USE_PVA:
  160.                 // Populate sample with position, velocity and acceleration data
  161.                 sample.forEach(pv -> {
  162.                     final Vector3D position     = pv.getPosition();
  163.                     final Vector3D velocity     = pv.getVelocity();
  164.                     final Vector3D acceleration = pv.getAcceleration();
  165.                     interpolator.addSamplePoint(pv.getDate().durationFrom(date),
  166.                                                 position.toArray(), velocity.toArray(), acceleration.toArray());
  167.                 });
  168.                 break;
  169.             default:
  170.                 // this should never happen
  171.                 throw new OrekitInternalError(null);
  172.         }

  173.         // interpolate
  174.         final double[][] pva = interpolator.derivatives(0.0, 2);

  175.         // build a new interpolated instance
  176.         return new AbsolutePVCoordinates(outputFrame, date, new Vector3D(pva[0]), new Vector3D(pva[1]),
  177.                                          new Vector3D(pva[2]));
  178.     }
  179. }