FieldShortPeriodicsInterpolatedCoefficient

/* Copyright 2002-2025 CS GROUP
 * Licensed to CS GROUP (CS) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * CS licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.orekit.propagation.semianalytical.dsst.utilities;

import java.util.ArrayList;

import org.hipparchus.Field;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.analysis.interpolation.FieldHermiteInterpolator;
import org.hipparchus.util.FastMath;
import org.orekit.time.FieldAbsoluteDate;

/** Interpolated short periodics coefficients.
 * <p>
 * Representation of a coefficient that need to be interpolated over time.
 * </p><p>
 * The short periodics coefficients can be interpolated for faster computation.
 * This class stores computed values of the coefficients through the method
 * {@link #addGridPoint} and gives an interpolated result through the method
 * {@link #value}.
 * </p>
 * @author Nicolas Bernard
 * @param <T> type of the field elements
 */
public class FieldShortPeriodicsInterpolatedCoefficient <T extends CalculusFieldElement<T>> {

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

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

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

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

    /**Simple constructor.
     * @param interpolationPoints number of points used in the interpolation
     */
    public FieldShortPeriodicsInterpolatedCoefficient(final int interpolationPoints) {
        this.interpolationPoints = interpolationPoints;
        this.abscissae = new ArrayList<>();
        this.values = new ArrayList<>();
        this.latestClosestNeighbor = 0;
    }

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

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

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

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

        //interpolation
        return interpolator.value(zero);

    }

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

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

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

                    if (lowerNeighborDistance.getReal() <= upperNeighborDistance.getReal()) {
                        inf--;
                    } else {
                        sup++;
                    }
                }
            }

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

        }

        return neighborsIndices;
    }

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

        //case where the date is before the available points
        if (date.compareTo(abscissae.get(0)) <= 0) {
            closestNeighbor = 0;
        }
        //case where the date is after the available points
        else if (date.compareTo(abscissae.get(abscissae.size() - 1)) >= 0) {
            closestNeighbor = abscissae.size() - 1;
        }
        //general case: one is looking for the two consecutives entries that surround the input date
        //then one choose the closest one
        else {
            int lowerBorder = latestClosestNeighbor;
            int upperBorder = latestClosestNeighbor;

            final int searchDirection = date.compareTo(abscissae.get(latestClosestNeighbor));
            if (searchDirection > 0) {
                upperBorder++;
                while (date.compareTo(abscissae.get(upperBorder)) > 0) {
                    upperBorder++;
                    lowerBorder++;
                }
            }
            else {
                lowerBorder--;
                while (date.compareTo(abscissae.get(lowerBorder)) < 0) {
                    upperBorder--;
                    lowerBorder--;
                }
            }

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

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

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

    /** Clear the recorded values from the interpolation grid.
     */
    public void clearHistory() {
        abscissae.clear();
        values.clear();
    }

    /** Add a point to the interpolation grid.
     * @param date abscissa of the point
     * @param value value of the element
     */
    public void addGridPoint(final FieldAbsoluteDate<T> date, final T[] value) {
        //If the grid is empty, the value is directly added to both arrays
        if (abscissae.isEmpty()) {
            abscissae.add(date);
            values.add(value);
        }
        //If the grid already contains this point, only its value is changed
        else if (abscissae.contains(date)) {
            values.set(abscissae.indexOf(date), value);
        }
        //If the grid does not contain this point, the position of the point
        //in the grid is computed first
        else {
            final int closestNeighbor = getClosestNeighbor(date);
            final int index = (date.compareTo(abscissae.get(closestNeighbor)) < 0) ? closestNeighbor : closestNeighbor + 1;
            abscissae.add(index, date);
            values.add(index, value);
        }
    }

}