Example usage for org.apache.commons.math3.analysis.interpolation HermiteInterpolator HermiteInterpolator

List of usage examples for org.apache.commons.math3.analysis.interpolation HermiteInterpolator HermiteInterpolator

Introduction

In this page you can find the example usage for org.apache.commons.math3.analysis.interpolation HermiteInterpolator HermiteInterpolator.

Prototype

public HermiteInterpolator() 

Source Link

Document

Create an empty interpolator.

Usage

From source file:be.ugent.maf.cellmissy.analysis.singlecell.processing.impl.interpolation.TrackHermiteInterpolator.java

@Override
public InterpolatedTrack interpolateTrack(double[] time, double[] x, double[] y) {
    // create interpolators for X and Y
    HermiteInterpolator xHermite = new HermiteInterpolator();
    HermiteInterpolator yHermite = new HermiteInterpolator();
    int interpolationPoints = PropertiesConfigurationHolder.getInstance().getInt("numberOfInterpolationPoints");

    // create arrays to hold the interpolant time, the interpolated X and the interpolated Y
    double[] interpolantTime = new double[interpolationPoints];
    double[] interpolatedX = new double[interpolationPoints];
    double[] interpolatedY = new double[interpolationPoints];
    // the step used for the interpolation in both direction
    double interpolationStep = (time[time.length - 1] - time[0]) / interpolationPoints;

    // check for monotonicity
    boolean monotonic = MathArrays.isMonotonic(time, MathArrays.OrderDirection.INCREASING, false);
    // in case time is not monotonic, sort in place time, x and y coordinates
    if (!monotonic) {
        MathArrays.sortInPlace(time, x, y);
    }/* w w w.ja  v  a2 s . co m*/

    double[] internalPointsDerivativeX = internalPointsDerivative(time, x);
    double[] internalPointsDerivativeY = internalPointsDerivative(time, y);

    double[] endPointsDerivativeX = endPointsDerivative(time, x);
    double[] endPointsDerivativeY = endPointsDerivative(time, y);

    // call the interpolator and add sample points to it
    // we do add only the values, and not their derivatives
    for (int i = 0; i < time.length; i++) {
        if (i == 0) {
            xHermite.addSamplePoint(time[i], new double[] { x[i] }, new double[] { endPointsDerivativeX[0] });
            yHermite.addSamplePoint(time[i], new double[] { y[i] }, new double[] { endPointsDerivativeY[0] });
        } else if (i == time.length - 1) {
            xHermite.addSamplePoint(time[i], new double[] { x[i] }, new double[] { endPointsDerivativeX[1] });
            yHermite.addSamplePoint(time[i], new double[] { y[i] }, new double[] { endPointsDerivativeY[1] });
        } else {
            xHermite.addSamplePoint(time[i], new double[] { x[i] },
                    new double[] { internalPointsDerivativeX[i - 1] });
            yHermite.addSamplePoint(time[i], new double[] { y[i] },
                    new double[] { internalPointsDerivativeY[i - 1] });
        }
    }

    for (int i = 0; i < interpolationPoints; i++) {
        interpolantTime[i] = time[0] + (i * interpolationStep);
        interpolatedX[i] = xHermite.value(interpolantTime[i])[0];
        interpolatedY[i] = yHermite.value(interpolantTime[i])[0];
    }

    // get the polynomial functions in both directions
    PolynomialFunction polynomialFunctionX = xHermite.getPolynomials()[0];
    PolynomialFunction polynomialFunctionY = yHermite.getPolynomials()[0];

    return new InterpolatedTrack(interpolantTime, interpolatedX, interpolatedY, polynomialFunctionX,
            polynomialFunctionY);
}

From source file:org.orekit.frames.EOPHistory.java

/** Get the UT1-UTC value.
 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
 * @param date date at which the value is desired
 * @return UT1-UTC in seconds (0 if date is outside covered range)
 *///from  ww w  .  ja  va2  s.com
public double getUT1MinusUTC(final AbsoluteDate date) {
    //check if there is data for date
    if (!this.hasDataFor(date)) {
        // no EOP data available for this date, we use a default 0.0 offset
        return (tidalCorrection == null) ? 0.0 : tidalCorrection.value(date)[2];
    }
    //we have EOP data -> interpolate offset
    try {
        final List<EOPEntry> neighbors = getNeighbors(date);
        final HermiteInterpolator interpolator = new HermiteInterpolator();
        final double firstDUT = neighbors.get(0).getUT1MinusUTC();
        boolean beforeLeap = true;
        for (final EOPEntry neighbor : neighbors) {
            final double dut;
            if (neighbor.getUT1MinusUTC() - firstDUT > 0.9) {
                // there was a leap second between the entries
                dut = neighbor.getUT1MinusUTC() - 1.0;
                if (neighbor.getDate().compareTo(date) <= 0) {
                    beforeLeap = false;
                }
            } else {
                dut = neighbor.getUT1MinusUTC();
            }
            interpolator.addSamplePoint(neighbor.getDate().durationFrom(date), new double[] { dut });
        }
        double interpolated = interpolator.value(0)[0];
        if (tidalCorrection != null) {
            interpolated += tidalCorrection.value(date)[2];
        }
        return beforeLeap ? interpolated : interpolated + 1.0;
    } catch (TimeStampedCacheException tce) {
        //this should not happen because of date check above
        throw new OrekitInternalError(tce);
    }
}

From source file:org.orekit.frames.EOPHistory.java

/** Get the LoD (Length of Day) value.
 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
 * @param date date at which the value is desired
 * @return LoD in seconds (0 if date is outside covered range)
 *///www .ja va2 s  .  c  om
public double getLOD(final AbsoluteDate date) {
    //check if there is data for date
    if (!this.hasDataFor(date)) {
        // no EOP data available for this date, we use a default null correction
        return (tidalCorrection == null) ? 0.0 : tidalCorrection.value(date)[3];
    }
    //we have EOP data for date -> interpolate correction
    try {
        final HermiteInterpolator interpolator = new HermiteInterpolator();
        for (final EOPEntry entry : getNeighbors(date)) {
            interpolator.addSamplePoint(entry.getDate().durationFrom(date), new double[] { entry.getLOD() });
        }
        double interpolated = interpolator.value(0)[0];
        if (tidalCorrection != null) {
            interpolated += tidalCorrection.value(date)[3];
        }
        return interpolated;
    } catch (TimeStampedCacheException tce) {
        // this should not happen because of date check above
        throw new OrekitInternalError(tce);
    }
}

From source file:org.orekit.frames.EOPHistory.java

/** Get the pole IERS Reference Pole correction.
 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
 * @param date date at which the correction is desired
 * @return pole correction ({@link PoleCorrection#NULL_CORRECTION
 * PoleCorrection.NULL_CORRECTION} if date is outside covered range)
 *//*from  w w  w .ja  va  2s .  c  om*/
public PoleCorrection getPoleCorrection(final AbsoluteDate date) {
    // check if there is data for date
    if (!this.hasDataFor(date)) {
        // no EOP data available for this date, we use a default null correction
        if (tidalCorrection == null) {
            return PoleCorrection.NULL_CORRECTION;
        } else {
            final double[] correction = tidalCorrection.value(date);
            return new PoleCorrection(correction[0], correction[1]);
        }
    }
    //we have EOP data for date -> interpolate correction
    try {
        final HermiteInterpolator interpolator = new HermiteInterpolator();
        for (final EOPEntry entry : getNeighbors(date)) {
            interpolator.addSamplePoint(entry.getDate().durationFrom(date),
                    new double[] { entry.getX(), entry.getY() });
        }
        final double[] interpolated = interpolator.value(0);
        if (tidalCorrection != null) {
            final double[] correction = tidalCorrection.value(date);
            interpolated[0] += correction[0];
            interpolated[1] += correction[1];
        }
        return new PoleCorrection(interpolated[0], interpolated[1]);
    } catch (TimeStampedCacheException tce) {
        // this should not happen because of date check above
        throw new OrekitInternalError(tce);
    }
}

From source file:org.orekit.frames.EOPHistory.java

/** Get the correction to the nutation parameters for equinox-based paradigm.
 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
 * @param date date at which the correction is desired
 * @return nutation correction in longitude  and in obliquity 
 * (zero if date is outside covered range)
 *//* www  .  j  a  v  a2 s.co m*/
public double[] getEquinoxNutationCorrection(final AbsoluteDate date) {
    // check if there is data for date
    if (!this.hasDataFor(date)) {
        // no EOP data available for this date, we use a default null correction
        return new double[2];
    }
    //we have EOP data for date -> interpolate correction
    try {
        final HermiteInterpolator interpolator = new HermiteInterpolator();
        for (final EOPEntry entry : getNeighbors(date)) {
            interpolator.addSamplePoint(entry.getDate().durationFrom(date),
                    new double[] { entry.getDdPsi(), entry.getDdEps() });
        }
        return interpolator.value(0);
    } catch (TimeStampedCacheException tce) {
        // this should not happen because of date check above
        throw new OrekitInternalError(tce);
    }
}

From source file:org.orekit.frames.EOPHistory.java

/** Get the correction to the nutation parameters for Non-Rotating Origin paradigm.
 * <p>The data provided comes from the IERS files. It is smoothed data.</p>
 * @param date date at which the correction is desired
 * @return nutation correction in Celestial Intermediat Pole coordinates
 * X and Y (zero if date is outside covered range)
 *///from   w  w w .j  a v  a 2s  .c o  m
public double[] getNonRotatinOriginNutationCorrection(final AbsoluteDate date) {
    // check if there is data for date
    if (!this.hasDataFor(date)) {
        // no EOP data available for this date, we use a default null correction
        return new double[2];
    }
    //we have EOP data for date -> interpolate correction
    try {
        final HermiteInterpolator interpolator = new HermiteInterpolator();
        for (final EOPEntry entry : getNeighbors(date)) {
            interpolator.addSamplePoint(entry.getDate().durationFrom(date),
                    new double[] { entry.getDx(), entry.getDy() });
        }
        return interpolator.value(0);
    } catch (TimeStampedCacheException tce) {
        // this should not happen because of date check above
        throw new OrekitInternalError(tce);
    }
}

From source file:org.orekit.models.earth.tessellation.AlongTrackAiming.java

/** {@inheritDoc} */
@Override//from ww w  . ja v  a  2s  .co  m
public Vector3D alongTileDirection(final Vector3D point, final GeodeticPoint gp) throws OrekitException {

    final double lStart = halfTrack.get(0).getFirst().getLatitude();
    final double lEnd = halfTrack.get(halfTrack.size() - 1).getFirst().getLatitude();
    // check the point can be reached
    if (gp.getLatitude() < FastMath.min(lStart, lEnd) || gp.getLatitude() > FastMath.max(lStart, lEnd)) {
        throw new OrekitException(OrekitMessages.OUT_OF_RANGE_LATITUDE, FastMath.toDegrees(gp.getLatitude()),
                FastMath.toDegrees(FastMath.min(lStart, lEnd)), FastMath.toDegrees(FastMath.max(lStart, lEnd)));
    }

    // bracket the point in the half track sample
    int iInf = 0;
    int iSup = halfTrack.size() - 1;
    while (iSup - iInf > 1) {
        final int iMiddle = (iSup + iInf) / 2;
        if ((lStart < lEnd) ^ (halfTrack.get(iMiddle).getFirst().getLatitude() > gp.getLatitude())) {
            // the specified latitude is in the second half
            iInf = iMiddle;
        } else {
            // the specified latitude is in the first half
            iSup = iMiddle;
        }
    }

    // ensure we can get points at iStart, iStart + 1, iStart + 2 and iStart + 3
    final int iStart = FastMath.max(0, FastMath.min(iInf - 1, halfTrack.size() - 4));

    // interpolate ground sliding point at specified latitude
    final HermiteInterpolator interpolator = new HermiteInterpolator();
    for (int i = iStart; i < iStart + 4; ++i) {
        final Vector3D position = halfTrack.get(i).getSecond().getPosition();
        final Vector3D velocity = halfTrack.get(i).getSecond().getVelocity();
        interpolator.addSamplePoint(halfTrack.get(i).getFirst().getLatitude(), new double[] { position.getX(),
                position.getY(), position.getZ(), velocity.getX(), velocity.getY(), velocity.getZ() });
    }
    final DerivativeStructure[] p = interpolator.value(new DerivativeStructure(1, 1, 0, gp.getLatitude()));

    // extract interpolated ground position/velocity
    final Vector3D position = new Vector3D(p[0].getValue(), p[1].getValue(), p[2].getValue());
    final Vector3D velocity = new Vector3D(p[3].getValue(), p[4].getValue(), p[5].getValue());

    // adjust longitude to match the specified one
    final Rotation rotation = new Rotation(Vector3D.PLUS_K, position, Vector3D.PLUS_K, point);
    final Vector3D fixedVelocity = rotation.applyTo(velocity);

    // the tile direction is aligned with sliding point velocity
    return fixedVelocity.normalize();

}

From source file:org.orekit.orbits.CircularOrbit.java

/** {@inheritDoc}
 * <p>//from w w  w .ja  v a2s.  c o  m
 * The interpolated instance is created by polynomial Hermite interpolation
 * on circular elements, without derivatives (which means the interpolation
 * falls back to Lagrange interpolation only).
 * </p>
 * <p>
 * As this implementation of interpolation is polynomial, it should be used only
 * with small samples (about 10-20 points) in order to avoid <a
 * href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
 * and numerical problems (including NaN appearing).
 * </p>
 * <p>
 * If orbit interpolation on large samples is needed, using the {@link
 * org.orekit.propagation.analytical.Ephemeris} class is a better way than using this
 * low-level interpolation. The Ephemeris class automatically handles selection of
 * a neighboring sub-sample with a predefined number of point from a large global sample
 * in a thread-safe way.
 * </p>
 */
public CircularOrbit interpolate(final AbsoluteDate date, final Collection<Orbit> sample) {

    // set up an interpolator
    final HermiteInterpolator interpolator = new HermiteInterpolator();

    // add sample points
    AbsoluteDate previousDate = null;
    double previousRAAN = Double.NaN;
    double previousAlphaM = Double.NaN;
    for (final Orbit orbit : sample) {
        final CircularOrbit circ = (CircularOrbit) OrbitType.CIRCULAR.convertType(orbit);
        final double continuousRAAN;
        final double continuousAlphaM;
        if (previousDate == null) {
            continuousRAAN = circ.getRightAscensionOfAscendingNode();
            continuousAlphaM = circ.getAlphaM();
        } else {
            final double dt = circ.getDate().durationFrom(previousDate);
            final double keplerAM = previousAlphaM + circ.getKeplerianMeanMotion() * dt;
            continuousRAAN = MathUtils.normalizeAngle(circ.getRightAscensionOfAscendingNode(), previousRAAN);
            continuousAlphaM = MathUtils.normalizeAngle(circ.getAlphaM(), keplerAM);
        }
        previousDate = circ.getDate();
        previousRAAN = continuousRAAN;
        previousAlphaM = continuousAlphaM;
        interpolator.addSamplePoint(circ.getDate().durationFrom(date), new double[] { circ.getA(),
                circ.getCircularEx(), circ.getCircularEy(), circ.getI(), continuousRAAN, continuousAlphaM });
    }

    // interpolate
    final double[] interpolated = interpolator.value(0);

    // build a new interpolated instance
    return new CircularOrbit(interpolated[0], interpolated[1], interpolated[2], interpolated[3],
            interpolated[4], interpolated[5], PositionAngle.MEAN, getFrame(), date, getMu());

}

From source file:org.orekit.orbits.EquinoctialOrbit.java

/** {@inheritDoc}
 * <p>/*from ww  w. java 2s  .com*/
 * The interpolated instance is created by polynomial Hermite interpolation
 * on equinoctial elements, without derivatives (which means the interpolation
 * falls back to Lagrange interpolation only).
 * </p>
 * <p>
 * As this implementation of interpolation is polynomial, it should be used only
 * with small samples (about 10-20 points) in order to avoid <a
 * href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
 * and numerical problems (including NaN appearing).
 * </p>
 * <p>
 * If orbit interpolation on large samples is needed, using the {@link
 * org.orekit.propagation.analytical.Ephemeris} class is a better way than using this
 * low-level interpolation. The Ephemeris class automatically handles selection of
 * a neighboring sub-sample with a predefined number of point from a large global sample
 * in a thread-safe way.
 * </p>
 */
public EquinoctialOrbit interpolate(final AbsoluteDate date, final Collection<Orbit> sample) {

    // set up an interpolator
    final HermiteInterpolator interpolator = new HermiteInterpolator();

    // add sample points
    AbsoluteDate previousDate = null;
    double previousLm = Double.NaN;
    for (final Orbit orbit : sample) {
        final EquinoctialOrbit equi = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(orbit);
        final double continuousLm;
        if (previousDate == null) {
            continuousLm = equi.getLM();
        } else {
            final double dt = equi.getDate().durationFrom(previousDate);
            final double keplerLm = previousLm + equi.getKeplerianMeanMotion() * dt;
            continuousLm = MathUtils.normalizeAngle(equi.getLM(), keplerLm);
        }
        previousDate = equi.getDate();
        previousLm = continuousLm;
        interpolator.addSamplePoint(equi.getDate().durationFrom(date), new double[] { equi.getA(),
                equi.getEquinoctialEx(), equi.getEquinoctialEy(), equi.getHx(), equi.getHy(), continuousLm });
    }

    // interpolate
    final double[] interpolated = interpolator.value(0);

    // build a new interpolated instance
    return new EquinoctialOrbit(interpolated[0], interpolated[1], interpolated[2], interpolated[3],
            interpolated[4], interpolated[5], PositionAngle.MEAN, getFrame(), date, getMu());

}

From source file:org.orekit.orbits.KeplerianOrbit.java

/** {@inheritDoc}
 * <p>//from ww  w .  j a v  a2s. c o  m
 * The interpolated instance is created by polynomial Hermite interpolation
 * on Keplerian elements, without derivatives (which means the interpolation
 * falls back to Lagrange interpolation only).
 * </p>
 * <p>
 * As this implementation of interpolation is polynomial, it should be used only
 * with small samples (about 10-20 points) in order to avoid <a
 * href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
 * and numerical problems (including NaN appearing).
 * </p>
 * <p>
 * If orbit interpolation on large samples is needed, using the {@link
 * org.orekit.propagation.analytical.Ephemeris} class is a better way than using this
 * low-level interpolation. The Ephemeris class automatically handles selection of
 * a neighboring sub-sample with a predefined number of point from a large global sample
 * in a thread-safe way.
 * </p>
 */
public KeplerianOrbit interpolate(final AbsoluteDate date, final Collection<Orbit> sample) {

    // set up an interpolator
    final HermiteInterpolator interpolator = new HermiteInterpolator();

    // add sample points
    AbsoluteDate previousDate = null;
    double previousPA = Double.NaN;
    double previousRAAN = Double.NaN;
    double previousM = Double.NaN;
    for (final Orbit orbit : sample) {
        final KeplerianOrbit kep = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(orbit);
        final double continuousPA;
        final double continuousRAAN;
        final double continuousM;
        if (previousDate == null) {
            continuousPA = kep.getPerigeeArgument();
            continuousRAAN = kep.getRightAscensionOfAscendingNode();
            continuousM = kep.getMeanAnomaly();
        } else {
            final double dt = kep.getDate().durationFrom(previousDate);
            final double keplerM = previousM + kep.getKeplerianMeanMotion() * dt;
            continuousPA = MathUtils.normalizeAngle(kep.getPerigeeArgument(), previousPA);
            continuousRAAN = MathUtils.normalizeAngle(kep.getRightAscensionOfAscendingNode(), previousRAAN);
            continuousM = MathUtils.normalizeAngle(kep.getMeanAnomaly(), keplerM);
        }
        previousDate = kep.getDate();
        previousPA = continuousPA;
        previousRAAN = continuousRAAN;
        previousM = continuousM;
        interpolator.addSamplePoint(kep.getDate().durationFrom(date),
                new double[] { kep.getA(), kep.getE(), kep.getI(), continuousPA, continuousRAAN, continuousM });
    }

    // interpolate
    final double[] interpolated = interpolator.value(0);

    // build a new interpolated instance
    return new KeplerianOrbit(interpolated[0], interpolated[1], interpolated[2], interpolated[3],
            interpolated[4], interpolated[5], PositionAngle.MEAN, getFrame(), date, getMu());

}