FieldShortPeriodicsInterpolatedCoefficient.java

  1. /* Copyright 2002-2022 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.propagation.semianalytical.dsst.utilities;

  18. import java.util.ArrayList;

  19. import org.hipparchus.Field;
  20. import org.hipparchus.CalculusFieldElement;
  21. import org.hipparchus.analysis.interpolation.FieldHermiteInterpolator;
  22. import org.hipparchus.util.FastMath;
  23. import org.orekit.time.FieldAbsoluteDate;

  24. /** Interpolated short periodics coefficients.
  25.  * <p>
  26.  * Representation of a coefficient that need to be interpolated over time.
  27.  * </p><p>
  28.  * The short periodics coefficients can be interpolated for faster computation.
  29.  * This class stores computed values of the coefficients through the method
  30.  * {@link #addGridPoint} and gives an interpolated result through the method
  31.  * {@link #value}.
  32.  * </p>
  33.  * @author Nicolas Bernard
  34.  *
  35.  */
  36. public class FieldShortPeriodicsInterpolatedCoefficient <T extends CalculusFieldElement<T>> {

  37.     /**Values of the already computed coefficients.*/
  38.     private ArrayList<T[]> values;

  39.     /**Grid points.*/
  40.     private ArrayList<FieldAbsoluteDate<T>> abscissae;

  41.     /**Number of points used in the interpolation.*/
  42.     private int interpolationPoints;

  43.     /**Index of the latest closest neighbor.*/
  44.     private int latestClosestNeighbor;

  45.     /**Simple constructor.
  46.      * @param interpolationPoints number of points used in the interpolation
  47.      */
  48.     public FieldShortPeriodicsInterpolatedCoefficient(final int interpolationPoints) {
  49.         this.interpolationPoints = interpolationPoints;
  50.         this.abscissae = new ArrayList<FieldAbsoluteDate<T>>();
  51.         this.values = new ArrayList<T[]>();
  52.         this.latestClosestNeighbor = 0;
  53.     }

  54.     /**Compute the value of the coefficient.
  55.      * @param date date at which the coefficient should be computed
  56.      * @return value of the coefficient
  57.      */
  58.     public T[] value(final FieldAbsoluteDate<T> date) {

  59.         //Field used by default
  60.         final Field<T> field = date.getField();
  61.         final T zero = field.getZero();

  62.         //Get the closest points from the input date
  63.         final int[] neighbors = getNeighborsIndices(date);

  64.         //Creation and set up of the interpolator
  65.         final FieldHermiteInterpolator<T> interpolator = new FieldHermiteInterpolator<>();
  66.         for (int i : neighbors) {
  67.             interpolator.addSamplePoint(abscissae.get(i).durationFrom(date), values.get(i));
  68.         }

  69.         //interpolation
  70.         return interpolator.value(zero);

  71.     }

  72.     /**Find the closest available points from the specified date.
  73.      * @param date date of interest
  74.      * @return indices corresponding to the closest points on the time scale
  75.      */
  76.     private int[] getNeighborsIndices(final FieldAbsoluteDate<T> date) {
  77.         final int sizeofNeighborhood = FastMath.min(interpolationPoints, abscissae.size());
  78.         final int[] neighborsIndices = new int[sizeofNeighborhood];

  79.         //If the size of the complete sample is less than
  80.         //the desired number of interpolation points,
  81.         //then the entire sample is considered as the neighborhood
  82.         if (interpolationPoints >= abscissae.size()) {
  83.             for (int i = 0; i < sizeofNeighborhood; i++) {
  84.                 neighborsIndices[i] = i;
  85.             }
  86.         } else {
  87.             // get indices around closest neighbor
  88.             int inf = getClosestNeighbor(date);
  89.             int sup = inf + 1;

  90.             while (sup - inf < interpolationPoints) {
  91.                 if (inf == 0) { //This means that we have reached the earliest date
  92.                     sup++;
  93.                 } else if (sup >= abscissae.size()) { //This means that we have reached the latest date
  94.                     inf--;
  95.                 } else { //the choice is made between the two next neighbors
  96.                     final T lowerNeighborDistance = FastMath.abs(abscissae.get(inf - 1).durationFrom(date));
  97.                     final T upperNeighborDistance = FastMath.abs(abscissae.get(sup).durationFrom(date));

  98.                     if (lowerNeighborDistance.getReal() <= upperNeighborDistance.getReal()) {
  99.                         inf--;
  100.                     } else {
  101.                         sup++;
  102.                     }
  103.                 }
  104.             }

  105.             for (int i = 0; i < interpolationPoints; ++i) {
  106.                 neighborsIndices[i] = inf + i;
  107.             }

  108.         }

  109.         return neighborsIndices;
  110.     }

  111.     /**Find the closest point from a specific date amongst the available points.
  112.      * @param date date of interest
  113.      * @return index of the closest abscissa from the date of interest
  114.      */
  115.     private int getClosestNeighbor(final FieldAbsoluteDate<T> date) {
  116.         //the starting point is the latest result of a call to this method.
  117.         //Indeed, as this class is meant to be called during an integration process
  118.         //with an input date evolving often continuously in time, there is a high
  119.         //probability that the result will be the same as for last call of
  120.         //this method.
  121.         int closestNeighbor = latestClosestNeighbor;

  122.         //case where the date is before the available points
  123.         if (date.compareTo(abscissae.get(0)) <= 0) {
  124.             closestNeighbor = 0;
  125.         }
  126.         //case where the date is after the available points
  127.         else if (date.compareTo(abscissae.get(abscissae.size() - 1)) >= 0) {
  128.             closestNeighbor = abscissae.size() - 1;
  129.         }
  130.         //general case: one is looking for the two consecutives entries that surround the input date
  131.         //then one choose the closest one
  132.         else {
  133.             int lowerBorder = latestClosestNeighbor;
  134.             int upperBorder = latestClosestNeighbor;

  135.             final int searchDirection = date.compareTo(abscissae.get(latestClosestNeighbor));
  136.             if (searchDirection > 0) {
  137.                 upperBorder++;
  138.                 while (date.compareTo(abscissae.get(upperBorder)) > 0) {
  139.                     upperBorder++;
  140.                     lowerBorder++;
  141.                 }
  142.             }
  143.             else {
  144.                 lowerBorder--;
  145.                 while (date.compareTo(abscissae.get(lowerBorder)) < 0) {
  146.                     upperBorder--;
  147.                     lowerBorder--;
  148.                 }
  149.             }

  150.             final T lowerDistance = FastMath.abs(date.durationFrom(abscissae.get(lowerBorder)));
  151.             final T upperDistance = FastMath.abs(date.durationFrom(abscissae.get(upperBorder)));

  152.             closestNeighbor = (lowerDistance.getReal() < upperDistance.getReal()) ? lowerBorder : upperBorder;
  153.         }

  154.         //The result is stored in order to speed up the next call to the function
  155.         //Indeed, it is highly likely that the requested result will be the same
  156.         this.latestClosestNeighbor = closestNeighbor;
  157.         return closestNeighbor;
  158.     }

  159.     /** Clear the recorded values from the interpolation grid.
  160.      */
  161.     public void clearHistory() {
  162.         abscissae.clear();
  163.         values.clear();
  164.     }

  165.     /** Add a point to the interpolation grid.
  166.      * @param date abscissa of the point
  167.      * @param value value of the element
  168.      */
  169.     public void addGridPoint(final FieldAbsoluteDate<T> date, final T[] value) {
  170.         //If the grid is empty, the value is directly added to both arrays
  171.         if (abscissae.isEmpty()) {
  172.             abscissae.add(date);
  173.             values.add(value);
  174.         }
  175.         //If the grid already contains this point, only its value is changed
  176.         else if (abscissae.contains(date)) {
  177.             values.set(abscissae.indexOf(date), value);
  178.         }
  179.         //If the grid does not contain this point, the position of the point
  180.         //in the grid is computed first
  181.         else {
  182.             final int closestNeighbor = getClosestNeighbor(date);
  183.             final int index = (date.compareTo(abscissae.get(closestNeighbor)) < 0) ? closestNeighbor : closestNeighbor + 1;
  184.             abscissae.add(index, date);
  185.             values.add(index, value);
  186.         }
  187.     }

  188. }