SecularAndHarmonic.java
/* Copyright 2002-2024 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.utils;
import java.util.ArrayList;
import java.util.Collection;
import java.util.List;
import org.hipparchus.analysis.ParametricUnivariateFunction;
import org.hipparchus.fitting.AbstractCurveFitter;
import org.hipparchus.fitting.PolynomialCurveFitter;
import org.hipparchus.fitting.WeightedObservedPoint;
import org.hipparchus.linear.DiagonalMatrix;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresBuilder;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.SinCos;
import org.orekit.time.AbsoluteDate;
/** Class for fitting evolution of osculating orbital parameters.
* <p>
* This class allows conversion from osculating parameters to mean parameters.
* </p>
*
* @author Luc Maisonobe
*/
public class SecularAndHarmonic {
/** Degree of polynomial secular part. */
private final int secularDegree;
/** Pulsations of harmonic part. */
private final double[] pulsations;
/** Reference date for the model. */
private AbsoluteDate reference;
/** Fitted parameters. */
private double[] fitted;
/** Observed points. */
private List<WeightedObservedPoint> observedPoints;
/** RMS for convergence.
* @since 10.3
*/
private double convergenceRMS;
/** Maximum number of iterations.
* @since 10.3
*/
private int maxIter;
/** Simple constructor.
* @param secularDegree degree of polynomial secular part
* @param pulsations pulsations of harmonic part
*/
public SecularAndHarmonic(final int secularDegree, final double... pulsations) {
this.secularDegree = secularDegree;
this.pulsations = pulsations.clone();
this.observedPoints = new ArrayList<WeightedObservedPoint>();
this.convergenceRMS = 0.0;
this.maxIter = Integer.MAX_VALUE;
}
/** Reset fitting.
* @param date reference date
* @param initialGuess initial guess for the parameters
* @see #getReferenceDate()
*/
public void resetFitting(final AbsoluteDate date, final double... initialGuess) {
reference = date;
fitted = initialGuess.clone();
observedPoints.clear();
}
/** Set RMS for convergence.
* <p>
* The RMS is the square-root of the sum of squared of
* the residuals, divided by the number of measurements.
* </p>
* @param convergenceRMS RMS below which convergence is considered to have been reached
* @since 10.3
*/
public void setConvergenceRMS(final double convergenceRMS) {
this.convergenceRMS = convergenceRMS;
}
/** Set maximum number of iterations.
* @param maxIter maximum number of iterations
* @since 10.3
*/
public void setMaxIter(final int maxIter) {
this.maxIter = maxIter;
}
/** Add a fitting point.
* <p>
* The point weight is set to 1.0
* </p>
* @param date date of the point
* @param osculatingValue osculating value
* @see #addWeightedPoint(AbsoluteDate, double, double)
*/
public void addPoint(final AbsoluteDate date, final double osculatingValue) {
addWeightedPoint(date, osculatingValue, 1.0);
}
/** Add a weighted fitting point.
* @param date date of the point
* @param osculatingValue osculating value
* @param weight weight of the points
* @since 12.0
*/
public void addWeightedPoint(final AbsoluteDate date, final double osculatingValue, final double weight) {
observedPoints.add(new WeightedObservedPoint(weight, date.durationFrom(reference), osculatingValue));
}
/** Get the reference date.
* @return reference date
* @see #resetFitting(AbsoluteDate, double...)
*/
public AbsoluteDate getReferenceDate() {
return reference;
}
/** Get degree of polynomial secular part.
* @return degree of polynomial secular part
* @since 12.0
*/
public int getSecularDegree() {
return secularDegree;
}
/** Get the pulsations of harmonic part.
* @return pulsations of harmonic part
* @since 12.0
*/
public double[] getPulsations() {
return pulsations.clone();
}
/** Get an upper bound of the fitted harmonic amplitude.
* @return upper bound of the fitted harmonic amplitude
*/
public double getHarmonicAmplitude() {
double amplitude = 0;
for (int i = 0; i < pulsations.length; ++i) {
amplitude += FastMath.hypot(fitted[secularDegree + 2 * i + 1],
fitted[secularDegree + 2 * i + 2]);
}
return amplitude;
}
/** Fit parameters.
* @see #getFittedParameters()
*/
public void fit() {
final AbstractCurveFitter fitter = new AbstractCurveFitter() {
/** {@inheritDoc} */
@Override
protected LeastSquaresProblem getProblem(final Collection<WeightedObservedPoint> observations) {
// Prepare least-squares problem.
final int len = observations.size();
final double[] target = new double[len];
final double[] weights = new double[len];
int i = 0;
for (final WeightedObservedPoint obs : observations) {
target[i] = obs.getY();
weights[i] = obs.getWeight();
++i;
}
final AbstractCurveFitter.TheoreticalValuesFunction model =
new AbstractCurveFitter.TheoreticalValuesFunction(new LocalParametricFunction(), observations);
// build a new least squares problem set up to fit a secular and harmonic curve to the observed points
return new LeastSquaresBuilder().
maxEvaluations(Integer.MAX_VALUE).
maxIterations(maxIter).
checker((iteration, previous, current) -> current.getRMS() <= convergenceRMS).
start(fitted).
target(target).
weight(new DiagonalMatrix(weights)).
model(model.getModelFunction(), model.getModelFunctionJacobian()).
build();
}
};
fitted = fitter.fit(observedPoints);
}
/** Local parametric function used for fitting. */
private class LocalParametricFunction implements ParametricUnivariateFunction {
/** {@inheritDoc} */
public double value(final double x, final double... parameters) {
return truncatedValue(secularDegree, pulsations.length, x, parameters);
}
/** {@inheritDoc} */
public double[] gradient(final double x, final double... parameters) {
final double[] gradient = new double[secularDegree + 1 + 2 * pulsations.length];
// secular part
double xN = 1.0;
for (int i = 0; i <= secularDegree; ++i) {
gradient[i] = xN;
xN *= x;
}
// harmonic part
for (int i = 0; i < pulsations.length; ++i) {
final SinCos sc = FastMath.sinCos(pulsations[i] * x);
gradient[secularDegree + 2 * i + 1] = sc.cos();
gradient[secularDegree + 2 * i + 2] = sc.sin();
}
return gradient;
}
}
/** Get a copy of the last fitted parameters.
* @return copy of the last fitted parameters.
* @see #fit()
*/
public double[] getFittedParameters() {
return fitted.clone();
}
/** Get fitted osculating value.
* @param date current date
* @return osculating value at current date
*/
public double osculatingValue(final AbsoluteDate date) {
return truncatedValue(secularDegree, pulsations.length,
date.durationFrom(reference), fitted);
}
/** Get fitted osculating derivative.
* @param date current date
* @return osculating derivative at current date
*/
public double osculatingDerivative(final AbsoluteDate date) {
return truncatedDerivative(secularDegree, pulsations.length,
date.durationFrom(reference), fitted);
}
/** Get fitted osculating second derivative.
* @param date current date
* @return osculating second derivative at current date
*/
public double osculatingSecondDerivative(final AbsoluteDate date) {
return truncatedSecondDerivative(secularDegree, pulsations.length,
date.durationFrom(reference), fitted);
}
/** Get mean value, truncated to first components.
* @param date current date
* @param degree degree of polynomial secular part to consider
* @param harmonics number of harmonics terms to consider
* @return mean value at current date
*/
public double meanValue(final AbsoluteDate date, final int degree, final int harmonics) {
return truncatedValue(degree, harmonics, date.durationFrom(reference), fitted);
}
/** Get mean derivative, truncated to first components.
* @param date current date
* @param degree degree of polynomial secular part to consider
* @param harmonics number of harmonics terms to consider
* @return mean derivative at current date
*/
public double meanDerivative(final AbsoluteDate date, final int degree, final int harmonics) {
return truncatedDerivative(degree, harmonics, date.durationFrom(reference), fitted);
}
/** Approximate an already fitted model to polynomial only terms.
* <p>
* This method is mainly used in order to combine the large amplitude long
* periods with the secular part as a new approximate polynomial model over
* some time range. This should be used rather than simply extracting the
* polynomial coefficients from {@link #getFittedParameters()} when some
* periodic terms amplitudes are large (for example Sun resonance effects
* on local solar time in sun synchronous orbits). In theses cases, the pure
* polynomial secular part in the coefficients may be far from the mean model.
* </p>
* @param combinedDegree desired degree for the combined polynomial
* @param combinedReference desired reference date for the combined polynomial
* @param meanDegree degree of polynomial secular part to consider
* @param meanHarmonics number of harmonics terms to consider
* @param start start date of the approximation time range
* @param end end date of the approximation time range
* @param step sampling step
* @return coefficients of the approximate polynomial (in increasing degree order),
* using the user provided reference date
*/
public double[] approximateAsPolynomialOnly(final int combinedDegree, final AbsoluteDate combinedReference,
final int meanDegree, final int meanHarmonics,
final AbsoluteDate start, final AbsoluteDate end,
final double step) {
final List<WeightedObservedPoint> points = new ArrayList<WeightedObservedPoint>();
for (AbsoluteDate date = start; date.compareTo(end) < 0; date = date.shiftedBy(step)) {
points.add(new WeightedObservedPoint(1.0,
date.durationFrom(combinedReference),
meanValue(date, meanDegree, meanHarmonics)));
}
return PolynomialCurveFitter.create(combinedDegree).fit(points);
}
/** Get mean second derivative, truncated to first components.
* @param date current date
* @param degree degree of polynomial secular part
* @param harmonics number of harmonics terms to consider
* @return mean second derivative at current date
*/
public double meanSecondDerivative(final AbsoluteDate date, final int degree, final int harmonics) {
return truncatedSecondDerivative(degree, harmonics, date.durationFrom(reference), fitted);
}
/** Get value truncated to first components.
* @param degree degree of polynomial secular part
* @param harmonics number of harmonics terms to consider
* @param time time parameter
* @param parameters models parameters (must include all parameters,
* including the ones ignored due to model truncation)
* @return truncated value
*/
private double truncatedValue(final int degree, final int harmonics,
final double time, final double... parameters) {
double value = 0;
// secular part
double tN = 1.0;
for (int i = 0; i <= degree; ++i) {
value += parameters[i] * tN;
tN *= time;
}
// harmonic part
for (int i = 0; i < harmonics; ++i) {
final SinCos sc = FastMath.sinCos(pulsations[i] * time);
value += parameters[secularDegree + 2 * i + 1] * sc.cos() +
parameters[secularDegree + 2 * i + 2] * sc.sin();
}
return value;
}
/** Get derivative truncated to first components.
* @param degree degree of polynomial secular part
* @param harmonics number of harmonics terms to consider
* @param time time parameter
* @param parameters models parameters (must include all parameters,
* including the ones ignored due to model truncation)
* @return truncated derivative
*/
private double truncatedDerivative(final int degree, final int harmonics,
final double time, final double... parameters) {
double derivative = 0;
// secular part
double tN = 1.0;
for (int i = 1; i <= degree; ++i) {
derivative += i * parameters[i] * tN;
tN *= time;
}
// harmonic part
for (int i = 0; i < harmonics; ++i) {
final SinCos sc = FastMath.sinCos(pulsations[i] * time);
derivative += pulsations[i] * (-parameters[secularDegree + 2 * i + 1] * sc.sin() +
parameters[secularDegree + 2 * i + 2] * sc.cos());
}
return derivative;
}
/** Get second derivative truncated to first components.
* @param degree degree of polynomial secular part
* @param harmonics number of harmonics terms to consider
* @param time time parameter
* @param parameters models parameters (must include all parameters,
* including the ones ignored due to model truncation)
* @return truncated second derivative
*/
private double truncatedSecondDerivative(final int degree, final int harmonics,
final double time, final double... parameters) {
double d2 = 0;
// secular part
double tN = 1.0;
for (int i = 2; i <= degree; ++i) {
d2 += (i - 1) * i * parameters[i] * tN;
tN *= time;
}
// harmonic part
for (int i = 0; i < harmonics; ++i) {
final SinCos sc = FastMath.sinCos(pulsations[i] * time);
d2 += -pulsations[i] * pulsations[i] *
(parameters[secularDegree + 2 * i + 1] * sc.cos() +
parameters[secularDegree + 2 * i + 2] * sc.sin());
}
return d2;
}
}