SecularAndHarmonic.java

  1. /* Copyright 2002-2018 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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;

  19. import java.util.ArrayList;
  20. import java.util.Collection;
  21. import java.util.List;

  23. import org.hipparchus.analysis.ParametricUnivariateFunction;
  24. import org.hipparchus.fitting.AbstractCurveFitter;
  25. import org.hipparchus.fitting.PolynomialCurveFitter;
  26. import org.hipparchus.fitting.WeightedObservedPoint;
  27. import org.hipparchus.linear.DiagonalMatrix;
  28. import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresBuilder;
  29. import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem;
  30. import org.hipparchus.util.FastMath;
  31. import org.orekit.time.AbsoluteDate;

  33. /** Class for fitting evolution of osculating orbital parameters.
  34.  * <p>
  35.  * This class allows conversion from osculating parameters to mean parameters.
  36.  * </p>
  37.  *
  38.  * @author Luc Maisonobe
  39.  */
  40. public class SecularAndHarmonic {

  42.     /** Degree of polynomial secular part. */
  43.     private final int secularDegree;

  45.     /** Pulsations of harmonic part. */
  46.     private final double[] pulsations;

  48.     /** Reference date for the model. */
  49.     private AbsoluteDate reference;

  51.     /** Fitted parameters. */
  52.     private double[] fitted;

  54.     /** Observed points. */
  55.     private List<WeightedObservedPoint> observedPoints;

  57.     /** Simple constructor.
  58.      * @param secularDegree degree of polynomial secular part
  59.      * @param pulsations pulsations of harmonic part
  60.      */
  61.     public SecularAndHarmonic(final int secularDegree, final double... pulsations) {
  62.         this.secularDegree  = secularDegree;
  63.         this.pulsations     = pulsations.clone();
  64.         this.observedPoints = new ArrayList<WeightedObservedPoint>();
  65.     }

  67.     /** Reset fitting.
  68.      * @param date reference date
  69.      * @param initialGuess initial guess for the parameters
  70.      * @see #getReferenceDate()
  71.      */
  72.     public void resetFitting(final AbsoluteDate date, final double... initialGuess) {
  73.         reference = date;
  74.         fitted    = initialGuess.clone();
  75.         observedPoints.clear();
  76.     }

  78.     /** Add a fitting point.
  79.      * @param date date of the point
  80.      * @param osculatingValue osculating value
  81.      */
  82.     public void addPoint(final AbsoluteDate date, final double osculatingValue) {
  83.         observedPoints.add(new WeightedObservedPoint(1.0, date.durationFrom(reference), osculatingValue));
  84.     }

  86.     /** Get the reference date.
  87.      * @return reference date
  88.      * @see #resetFitting(AbsoluteDate, double...)
  89.      */
  90.     public AbsoluteDate getReferenceDate() {
  91.         return reference;
  92.     }

  94.     /** Get an upper bound of the fitted harmonic amplitude.
  95.      * @return upper bound of the fitted harmonic amplitude
  96.      */
  97.     public double getHarmonicAmplitude() {
  98.         double amplitude = 0;
  99.         for (int i = 0; i < pulsations.length; ++i) {
  100.             amplitude += FastMath.hypot(fitted[secularDegree + 2 * i + 1],
  101.                                         fitted[secularDegree + 2 * i + 2]);
  102.         }
  103.         return amplitude;
  104.     }

  106.     /** Fit parameters.
  107.      * @see #getFittedParameters()
  108.      */
  109.     public void fit() {

  111.         final AbstractCurveFitter fitter = new AbstractCurveFitter() {
  112.             /** {@inheritDoc} */
  113.             @Override
  114.             protected LeastSquaresProblem getProblem(final Collection<WeightedObservedPoint> observations) {
  115.                 // Prepare least-squares problem.
  116.                 final int len = observations.size();
  117.                 final double[] target  = new double[len];
  118.                 final double[] weights = new double[len];

  120.                 int i = 0;
  121.                 for (final WeightedObservedPoint obs : observations) {
  122.                     target[i]  = obs.getY();
  123.                     weights[i] = obs.getWeight();
  124.                     ++i;
  125.                 }

  127.                 final AbstractCurveFitter.TheoreticalValuesFunction model =
  128.                         new AbstractCurveFitter.TheoreticalValuesFunction(new LocalParametricFunction(), observations);

  130.                 // build a new least squares problem set up to fit a secular and harmonic curve to the observed points
  131.                 return new LeastSquaresBuilder().
  132.                         maxEvaluations(Integer.MAX_VALUE).
  133.                         maxIterations(Integer.MAX_VALUE).
  134.                         start(fitted).
  135.                         target(target).
  136.                         weight(new DiagonalMatrix(weights)).
  137.                         model(model.getModelFunction(), model.getModelFunctionJacobian()).
  138.                         build();

  140.             }
  141.         };

  143.         fitted = fitter.fit(observedPoints);

  145.     }

  147.     /** Local parametric function used for fitting. */
  148.     private class LocalParametricFunction implements ParametricUnivariateFunction {

  150.         /** {@inheritDoc} */
  151.         public double value(final double x, final double... parameters) {
  152.             return truncatedValue(secularDegree, pulsations.length, x, parameters);
  153.         }

  155.         /** {@inheritDoc} */
  156.         public double[] gradient(final double x, final double... parameters) {
  157.             final double[] gradient = new double[secularDegree + 1 + 2 * pulsations.length];

  159.             // secular part
  160.             double xN = 1.0;
  161.             for (int i = 0; i <= secularDegree; ++i) {
  162.                 gradient[i] = xN;
  163.                 xN *= x;
  164.             }

  166.             // harmonic part
  167.             for (int i = 0; i < pulsations.length; ++i) {
  168.                 gradient[secularDegree + 2 * i + 1] = FastMath.cos(pulsations[i] * x);
  169.                 gradient[secularDegree + 2 * i + 2] = FastMath.sin(pulsations[i] * x);
  170.             }

  172.             return gradient;
  173.         }

  175.     }

  177.     /** Get a copy of the last fitted parameters.
  178.      * @return copy of the last fitted parameters.
  179.      * @see #fit()
  180.      */
  181.     public double[] getFittedParameters() {
  182.         return fitted.clone();
  183.     }

  185.     /** Get fitted osculating value.
  186.      * @param date current date
  187.      * @return osculating value at current date
  188.      */
  189.     public double osculatingValue(final AbsoluteDate date) {
  190.         return truncatedValue(secularDegree, pulsations.length,
  191.                               date.durationFrom(reference), fitted);
  192.     }

  194.     /** Get fitted osculating derivative.
  195.      * @param date current date
  196.      * @return osculating derivative at current date
  197.      */
  198.     public double osculatingDerivative(final AbsoluteDate date) {
  199.         return truncatedDerivative(secularDegree, pulsations.length,
  200.                                    date.durationFrom(reference), fitted);
  201.     }

  203.     /** Get fitted osculating second derivative.
  204.      * @param date current date
  205.      * @return osculating second derivative at current date
  206.      */
  207.     public double osculatingSecondDerivative(final AbsoluteDate date) {
  208.         return truncatedSecondDerivative(secularDegree, pulsations.length,
  209.                                          date.durationFrom(reference), fitted);
  210.     }

  212.     /** Get mean value, truncated to first components.
  213.      * @param date current date
  214.      * @param degree degree of polynomial secular part to consider
  215.      * @param harmonics number of harmonics terms to consider
  216.      * @return mean value at current date
  217.      */
  218.     public double meanValue(final AbsoluteDate date, final int degree, final int harmonics) {
  219.         return truncatedValue(degree, harmonics, date.durationFrom(reference), fitted);
  220.     }

  222.     /** Get mean derivative, truncated to first components.
  223.      * @param date current date
  224.      * @param degree degree of polynomial secular part to consider
  225.      * @param harmonics number of harmonics terms to consider
  226.      * @return mean derivative at current date
  227.      */
  228.     public double meanDerivative(final AbsoluteDate date, final int degree, final int harmonics) {
  229.         return truncatedDerivative(degree, harmonics, date.durationFrom(reference), fitted);
  230.     }

  232.     /** Approximate an already fitted model to polynomial only terms.
  233.      * <p>
  234.      * This method is mainly used in order to combine the large amplitude long
  235.      * periods with the secular part as a new approximate polynomial model over
  236.      * some time range. This should be used rather than simply extracting the
  237.      * polynomial coefficients from {@link #getFittedParameters()} when some
  238.      * periodic terms amplitudes are large (for example Sun resonance effects
  239.      * on local solar time in sun synchronous orbits). In theses cases, the pure
  240.      * polynomial secular part in the coefficients may be far from the mean model.
  241.      * </p>
  242.      * @param combinedDegree desired degree for the combined polynomial
  243.      * @param combinedReference desired reference date for the combined polynomial
  244.      * @param meanDegree degree of polynomial secular part to consider
  245.      * @param meanHarmonics number of harmonics terms to consider
  246.      * @param start start date of the approximation time range
  247.      * @param end end date of the approximation time range
  248.      * @param step sampling step
  249.      * @return coefficients of the approximate polynomial (in increasing degree order),
  250.      * using the user provided reference date
  251.      */
  252.     public double[] approximateAsPolynomialOnly(final int combinedDegree, final AbsoluteDate combinedReference,
  253.                                                 final int meanDegree, final int meanHarmonics,
  254.                                                 final AbsoluteDate start, final AbsoluteDate end,
  255.                                                 final double step) {
  256.         final List<WeightedObservedPoint> points = new ArrayList<WeightedObservedPoint>();
  257.         for (AbsoluteDate date = start; date.compareTo(end) < 0; date = date.shiftedBy(step)) {
  258.             points.add(new WeightedObservedPoint(1.0,
  259.                                                  date.durationFrom(combinedReference),
  260.                                                  meanValue(date, meanDegree, meanHarmonics)));
  261.         }
  262.         return PolynomialCurveFitter.create(combinedDegree).fit(points);
  263.     }

  265.     /** Get mean second derivative, truncated to first components.
  266.      * @param date current date
  267.      * @param degree degree of polynomial secular part
  268.      * @param harmonics number of harmonics terms to consider
  269.      * @return mean second derivative at current date
  270.      */
  271.     public double meanSecondDerivative(final AbsoluteDate date, final int degree, final int harmonics) {
  272.         return truncatedSecondDerivative(degree, harmonics, date.durationFrom(reference), fitted);
  273.     }

  275.     /** Get value truncated to first components.
  276.      * @param degree degree of polynomial secular part
  277.      * @param harmonics number of harmonics terms to consider
  278.      * @param time time parameter
  279.      * @param parameters models parameters (must include all parameters,
  280.      * including the ones ignored due to model truncation)
  281.      * @return truncated value
  282.      */
  283.     private double truncatedValue(final int degree, final int harmonics,
  284.                                   final double time, final double... parameters) {

  286.         double value = 0;

  288.         // secular part
  289.         double tN = 1.0;
  290.         for (int i = 0; i <= degree; ++i) {
  291.             value += parameters[i] * tN;
  292.             tN    *= time;
  293.         }

  295.         // harmonic part
  296.         for (int i = 0; i < harmonics; ++i) {
  297.             value += parameters[secularDegree + 2 * i + 1] * FastMath.cos(pulsations[i] * time) +
  298.                      parameters[secularDegree + 2 * i + 2] * FastMath.sin(pulsations[i] * time);
  299.         }

  301.         return value;

  303.     }

  305.     /** Get derivative truncated to first components.
  306.      * @param degree degree of polynomial secular part
  307.      * @param harmonics number of harmonics terms to consider
  308.      * @param time time parameter
  309.      * @param parameters models parameters (must include all parameters,
  310.      * including the ones ignored due to model truncation)
  311.      * @return truncated derivative
  312.      */
  313.     private double truncatedDerivative(final int degree, final int harmonics,
  314.                                        final double time, final double... parameters) {

  316.         double derivative = 0;

  318.         // secular part
  319.         double tN = 1.0;
  320.         for (int i = 1; i <= degree; ++i) {
  321.             derivative += i * parameters[i] * tN;
  322.             tN         *= time;
  323.         }

  325.         // harmonic part
  326.         for (int i = 0; i < harmonics; ++i) {
  327.             derivative += pulsations[i] * (-parameters[secularDegree + 2 * i + 1] * FastMath.sin(pulsations[i] * time) +
  328.                                             parameters[secularDegree + 2 * i + 2] * FastMath.cos(pulsations[i] * time));
  329.         }

  331.         return derivative;

  333.     }

  335.     /** Get second derivative truncated to first components.
  336.      * @param degree degree of polynomial secular part
  337.      * @param harmonics number of harmonics terms to consider
  338.      * @param time time parameter
  339.      * @param parameters models parameters (must include all parameters,
  340.      * including the ones ignored due to model truncation)
  341.      * @return truncated second derivative
  342.      */
  343.     private double truncatedSecondDerivative(final int degree, final int harmonics,
  344.                                              final double time, final double... parameters) {

  346.         double d2 = 0;

  348.         // secular part
  349.         double tN = 1.0;
  350.         for (int i = 2; i <= degree; ++i) {
  351.             d2 += (i - 1) * i * parameters[i] * tN;
  352.             tN *= time;
  353.         }

  355.         // harmonic part
  356.         for (int i = 0; i < harmonics; ++i) {
  357.             d2 += -pulsations[i] * pulsations[i] *
  358.                   (parameters[secularDegree + 2 * i + 1] * FastMath.cos(pulsations[i] * time) +
  359.                    parameters[secularDegree + 2 * i + 2] * FastMath.sin(pulsations[i] * time));
  360.         }

  362.         return d2;

  364.     }

  366. }
Created with JaCoCo 0.7.9.201702052155