SingleParameterFitter.java

  1. /* Copyright 2023 Luc Maisonobe
  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.frames;

  18. import java.io.Serializable;
  19. import java.util.List;
  20. import java.util.ListIterator;
  21. import java.util.function.ToDoubleFunction;

  22. import org.hipparchus.util.FastMath;
  23. import org.hipparchus.util.MathUtils;
  24. import org.orekit.time.AbsoluteDate;
  25. import org.orekit.utils.Constants;
  26. import org.orekit.utils.SecularAndHarmonic;

  27. /** Fitter for one Earth Orientation Parameter.
  28.  * @see PredictedEOPHistory
  29.  * @see EOPFitter
  30.  * @see SecularAndHarmonic
  31.  * @since 12.0
  32.  * @author Luc Maisonobe
  33.  */
  34. public class SingleParameterFitter implements Serializable {

  35.     /** Sun pulsation, one year period. */
  36.     public static final double SUN_PULSATION = MathUtils.TWO_PI / Constants.JULIAN_YEAR;

  37.     /** Moon pulsation (one Moon draconic period). */
  38.     public static final double MOON_DRACONIC_PULSATION = MathUtils.TWO_PI / (27.212221 * Constants.JULIAN_DAY);

  39.     /** Serializable UID. */
  40.     private static final long serialVersionUID = 20230309L;

  41.     /** Duration of the fitting window at the end of the raw history (s). */
  42.     private final double fittingDuration;

  43.     /** Time constant of the exponential decay weight. */
  44.     private final double timeConstant;

  45.     /** Convergence on fitted parameter. */
  46.     private final double convergence;

  47.     /** Degree of the polynomial model. */
  48.     private final int degree;

  49.     /** Pulsations of harmonic part (rad/s). */
  50.     private final double[] pulsations;

  51.     /** Simple constructor.
  52.      * @param fittingDuration duration of the fitting window at the end of the raw history (s)
  53.      * @param timeConstant time constant \(\tau\) of the exponential decay weight, point weight is \(e^{\frac{t-t_0}{\tau}}\),
  54.      * i.e. points far in the past before \(t_0\) have smaller weights
  55.      * @param convergence convergence on fitted parameter
  56.      * @param degree degree of the polynomial model
  57.      * @param pulsations pulsations of harmonic part (rad/s)
  58.      * @see #createDefaultDut1FitterShortTermPrediction()
  59.      * @see #createDefaultDut1FitterLongTermPrediction()
  60.      * @see #createDefaultPoleFitterShortTermPrediction()
  61.      * @see #createDefaultPoleFitterLongTermPrediction()
  62.      * @see #createDefaultNutationFitterShortTermPrediction()
  63.      * @see #createDefaultNutationFitterLongTermPrediction()
  64.      * @see SecularAndHarmonic
  65.      */
  66.     public SingleParameterFitter(final double fittingDuration, final double timeConstant, final double convergence,
  67.                                  final int degree, final double... pulsations) {
  68.         this.fittingDuration = fittingDuration;
  69.         this.timeConstant    = timeConstant;
  70.         this.convergence     = convergence;
  71.         this.degree          = degree;
  72.         this.pulsations      = pulsations.clone();
  73.     }

  74.     /** Perform secular and harmonic fitting.
  75.      * @param rawHistory EOP history to fit
  76.      * @param extractor extractor for Earth Orientation Parameter
  77.      * @return configured fitter
  78.      */
  79.     public SecularAndHarmonic fit(final EOPHistory rawHistory, final ToDoubleFunction<EOPEntry> extractor) {

  80.         final List<EOPEntry> rawEntries = rawHistory.getEntries();
  81.         final EOPEntry       last       = rawEntries.get(rawEntries.size() - 1);

  82.         // create fitter
  83.         final SecularAndHarmonic sh = new SecularAndHarmonic(degree, pulsations);

  84.         // set up convergence
  85.         sh.setConvergenceRMS(convergence);

  86.         // set up reference date and initial guess to a constant value
  87.         final double[] initialGuess = new double[degree + 1 + 2 * pulsations.length];
  88.         initialGuess[0] = extractor.applyAsDouble(last);
  89.         sh.resetFitting(last.getDate(), initialGuess);

  90.         // sample history
  91.         final AbsoluteDate           fitStart         = last.getDate().shiftedBy(-fittingDuration);
  92.         final ListIterator<EOPEntry> backwardIterator = rawEntries.listIterator(rawEntries.size());
  93.         while (backwardIterator.hasPrevious()) {
  94.             final EOPEntry entry = backwardIterator.previous();
  95.             if (entry.getDate().isAfterOrEqualTo(fitStart)) {
  96.                 // the entry belongs to the fitting interval
  97.                 sh.addWeightedPoint(entry.getDate(), extractor.applyAsDouble(entry),
  98.                                     FastMath.exp(entry.getDate().durationFrom(fitStart) / timeConstant));
  99.             } else {
  100.                 // we have processed all entries from the fitting interval
  101.                 break;
  102.             }
  103.         }

  104.         // perform fitting
  105.         sh.fit();

  106.         return sh;

  107.     }

  108.     /** Create fitter with default parameters adapted for fitting orientation parameters dUT1 and LOD
  109.      * for short term prediction.
  110.      * <p>
  111.      * The main difference between these settings and {@link #createDefaultDut1FitterLongTermPrediction()
  112.      * the settings for long prediction} is the much smaller \(\tau\). This means more
  113.      * weight is set to the points at the end of the history, hence forcing the fitted prediction
  114.      * model to be closer to these points, hence the prediction error to be smaller just after
  115.      * raw history end. On the other hand, this implies that the model will diverge on long term.
  116.      * These settings are intended when prediction is used for at most 5 days after raw EOP end.
  117.      * </p>
  118.      * <ul>
  119.      *   <li>fitting duration set to one {@link Constants#JULIAN_YEAR year}</li>
  120.      *   <li>time constant \(\tau\) of the exponential decay set to 6 {@link Constants#JULIAN_DAY days}</li>
  121.      *   <li>convergence set to 10⁻¹² s</li>
  122.      *   <li>polynomial part set to degree 3</li>
  123.      *   <li>one harmonic term at {@link #SUN_PULSATION}}</li>
  124.      *   <li>one harmonic term at 2 times {@link #SUN_PULSATION}}</li>
  125.      *   <li>one harmonic term at 3 times {@link #SUN_PULSATION}}</li>
  126.      *   <li>one harmonic term at {@link #MOON_DRACONIC_PULSATION}}</li>
  127.      *   <li>one harmonic term at 2 times {@link #MOON_DRACONIC_PULSATION}}</li>
  128.      *   <li>one harmonic term at 3 times {@link #MOON_DRACONIC_PULSATION}}</li>
  129.      * </ul>
  130.      * @return fitter with default configuration for orientation parameters dUT1 and LOD
  131.      * @see #createDefaultDut1FitterShortTermPrediction()
  132.      */
  133.     public static SingleParameterFitter createDefaultDut1FitterShortTermPrediction() {
  134.         return new SingleParameterFitter(Constants.JULIAN_YEAR, 6 * Constants.JULIAN_DAY, 1.0e-12, 3,
  135.                              SUN_PULSATION, 2 * SUN_PULSATION, 3 * SUN_PULSATION,
  136.                              MOON_DRACONIC_PULSATION, 2 * MOON_DRACONIC_PULSATION, 3 * MOON_DRACONIC_PULSATION);
  137.     }

  138.     /** Create fitter with default parameters adapted for fitting orientation parameters dUT1 and LOD
  139.      * for long term prediction.
  140.      * <p>
  141.      * The main difference between these settings and {@link #createDefaultDut1FitterShortTermPrediction()
  142.      * the settings for short prediction} is the much larger \(\tau\). This means weight
  143.      * is spread throughout history, hence forcing the fitted prediction model to be remain very stable
  144.      * on the long term. On the other hand, this implies that the model will start with already a much
  145.      * larger error just after raw history end.
  146.      * These settings are intended when prediction is used for 5 days after raw EOP end or more.
  147.      * </p>
  148.      * <ul>
  149.      *   <li>fitting duration set to three {@link Constants#JULIAN_YEAR years}</li>
  150.      *   <li>time constant \(\tau\) of the exponential decay set to 60 {@link Constants#JULIAN_DAY days}</li>
  151.      *   <li>convergence set to 10⁻¹² s</li>
  152.      *   <li>polynomial part set to degree 3</li>
  153.      *   <li>one harmonic term at {@link #SUN_PULSATION}}</li>
  154.      *   <li>one harmonic term at 2 times {@link #SUN_PULSATION}}</li>
  155.      *   <li>one harmonic term at 3 times {@link #SUN_PULSATION}}</li>
  156.      *   <li>one harmonic term at {@link #MOON_DRACONIC_PULSATION}}</li>
  157.      *   <li>one harmonic term at 2 times {@link #MOON_DRACONIC_PULSATION}}</li>
  158.      *   <li>one harmonic term at 3 times {@link #MOON_DRACONIC_PULSATION}}</li>
  159.      * </ul>
  160.      * @return fitter with default configuration for orientation parameters dUT1 and LOD
  161.      * @see #createDefaultDut1FitterShortTermPrediction()
  162.      */
  163.     public static SingleParameterFitter createDefaultDut1FitterLongTermPrediction() {
  164.         return new SingleParameterFitter(3 * Constants.JULIAN_YEAR, 60 * Constants.JULIAN_DAY, 1.0e-12, 3,
  165.                              SUN_PULSATION, 2 * SUN_PULSATION, 3 * SUN_PULSATION,
  166.                              MOON_DRACONIC_PULSATION, 2 * MOON_DRACONIC_PULSATION, 3 * MOON_DRACONIC_PULSATION);
  167.     }

  168.     /** Create fitter with default parameters adapted for fitting pole parameters Xp and Yp
  169.      * for long term prediction.
  170.      * <p>
  171.      * The main difference between these settings and {@link #createDefaultPoleFitterLongTermPrediction()
  172.      * the settings for long prediction} is the much smaller \(\tau\). This means more
  173.      * weight is set to the points at the end of the history, hence forcing the fitted prediction
  174.      * model to be closer to these points, hence the prediction error to be smaller just after
  175.      * raw history end. On the other hand, this implies that the model will diverge on long term.
  176.      * These settings are intended when prediction is used for at most 5 days after raw EOP end.
  177.      * </p>
  178.      * <ul>
  179.      *   <li>fitting duration set to one {@link Constants#JULIAN_YEAR year}</li>
  180.      *   <li>time constant \(\tau\) of the exponential decay set to 12 {@link Constants#JULIAN_DAY days}</li>
  181.      *   <li>convergence set to 10⁻¹² rad</li>
  182.      *   <li>polynomial part set to degree 3</li>
  183.      *   <li>one harmonic term at {@link #SUN_PULSATION}}</li>
  184.      *   <li>one harmonic term at 2 times {@link #SUN_PULSATION}}</li>
  185.      *   <li>one harmonic term at 3 times {@link #SUN_PULSATION}}</li>
  186.      *   <li>one harmonic term at {@link #MOON_DRACONIC_PULSATION}}</li>
  187.      *   <li>one harmonic term at 2 times {@link #MOON_DRACONIC_PULSATION}}</li>
  188.      *   <li>one harmonic term at 3 times {@link #MOON_DRACONIC_PULSATION}}</li>
  189.      * </ul>
  190.      * @return fitter with default configuration for pole parameters Xp and Yp
  191.      */
  192.     public static SingleParameterFitter createDefaultPoleFitterShortTermPrediction() {
  193.         return new SingleParameterFitter(Constants.JULIAN_YEAR, 12 * Constants.JULIAN_DAY, 1.0e-12, 3,
  194.                              SUN_PULSATION, 2 * SUN_PULSATION, 3 * SUN_PULSATION,
  195.                              MOON_DRACONIC_PULSATION, 2 * MOON_DRACONIC_PULSATION, 3 * MOON_DRACONIC_PULSATION);
  196.     }

  197.     /** Create fitter with default parameters adapted for fitting pole parameters Xp and Yp
  198.      * for long term prediction.
  199.      * <p>
  200.      * The main difference between these settings and {@link #createDefaultPoleFitterShortTermPrediction()
  201.      * the settings for short prediction} is the much larger \(\tau\). This means weight
  202.      * is spread throughout history, hence forcing the fitted prediction model to be remain very stable
  203.      * on the long term. On the other hand, this implies that the model will start with already a much
  204.      * larger error just after raw history end.
  205.      * These settings are intended when prediction is used for 5 days after raw EOP end or more.
  206.      * </p>
  207.      * <ul>
  208.      *   <li>fitting duration set to three {@link Constants#JULIAN_YEAR years}</li>
  209.      *   <li>time constant \(\tau\) of the exponential decay set to 60 {@link Constants#JULIAN_DAY days}</li>
  210.      *   <li>convergence set to 10⁻¹² rad</li>
  211.      *   <li>polynomial part set to degree 3</li>
  212.      *   <li>one harmonic term at {@link #SUN_PULSATION}}</li>
  213.      *   <li>one harmonic term at 2 times {@link #SUN_PULSATION}}</li>
  214.      *   <li>one harmonic term at 3 times {@link #SUN_PULSATION}}</li>
  215.      *   <li>one harmonic term at {@link #MOON_DRACONIC_PULSATION}}</li>
  216.      *   <li>one harmonic term at 2 times {@link #MOON_DRACONIC_PULSATION}}</li>
  217.      *   <li>one harmonic term at 3 times {@link #MOON_DRACONIC_PULSATION}}</li>
  218.      * </ul>
  219.      * @return fitter with default configuration for pole parameters Xp and Yp
  220.      */
  221.     public static SingleParameterFitter createDefaultPoleFitterLongTermPrediction() {
  222.         return new SingleParameterFitter(3 * Constants.JULIAN_YEAR, 60 * Constants.JULIAN_DAY, 1.0e-12, 3,
  223.                              SUN_PULSATION, 2 * SUN_PULSATION, 3 * SUN_PULSATION,
  224.                              MOON_DRACONIC_PULSATION, 2 * MOON_DRACONIC_PULSATION, 3 * MOON_DRACONIC_PULSATION);
  225.     }

  226.     /** Create fitter with default parameters adapted for fitting nutation parameters dx and dy
  227.      * for long term prediction.
  228.      * <p>
  229.      * The main difference between these settings and {@link #createDefaultNutationFitterLongTermPrediction()
  230.      * the settings for long prediction} is the much smaller \(\tau\). This means more
  231.      * weight is set to the points at the end of the history, hence forcing the fitted prediction
  232.      * model to be closer to these points, hence the prediction error to be smaller just after
  233.      * raw history end. On the other hand, this implies that the model will diverge on long term.
  234.      * These settings are intended when prediction is used for at most 5 days after raw EOP end.
  235.      * </p>
  236.      * <ul>
  237.      *   <li>fitting duration set to one {@link Constants#JULIAN_YEAR year}</li>
  238.      *   <li>time constant \(\tau\) of the exponential decay set to 12 {@link Constants#JULIAN_DAY days}</li>
  239.      *   <li>convergence set to 10⁻¹² s</li>
  240.      *   <li>polynomial part set to degree 3</li>
  241.      *   <li>one harmonic term at {@link #SUN_PULSATION}}</li>
  242.      *   <li>one harmonic term at 2 times {@link #SUN_PULSATION}}</li>
  243.      *   <li>one harmonic term at 3 times {@link #SUN_PULSATION}}</li>
  244.      *   <li>one harmonic term at {@link #MOON_DRACONIC_PULSATION}}</li>
  245.      *   <li>one harmonic term at 2 times {@link #MOON_DRACONIC_PULSATION}}</li>
  246.      *   <li>one harmonic term at 3 times {@link #MOON_DRACONIC_PULSATION}}</li>
  247.      * </ul>
  248.      * @return fitter with default configuration for pole nutation parameters dx and dy
  249.      */
  250.     public static SingleParameterFitter createDefaultNutationFitterShortTermPrediction() {
  251.         return new SingleParameterFitter(Constants.JULIAN_YEAR, 12 * Constants.JULIAN_DAY, 1.0e-12, 3,
  252.                              SUN_PULSATION, 2 * SUN_PULSATION, 3 * SUN_PULSATION,
  253.                              MOON_DRACONIC_PULSATION, 2 * MOON_DRACONIC_PULSATION, 3 * MOON_DRACONIC_PULSATION);
  254.     }

  255.     /** Create fitter with default parameters adapted for fitting nutation parameters dx and dy
  256.      * for long term prediction.
  257.      * <p>
  258.      * The main difference between these settings and {@link #createDefaultNutationFitterShortTermPrediction()
  259.      * the settings for short prediction} is the much larger \(\tau\). This means weight
  260.      * is spread throughout history, hence forcing the fitted prediction model to be remain very stable
  261.      * on the long term. On the other hand, this implies that the model will start with already a much
  262.      * larger error just after raw history end.
  263.      * These settings are intended when prediction is used for 5 days after raw EOP end or more.
  264.      * </p>
  265.      * <ul>
  266.      *   <li>fitting duration set to three {@link Constants#JULIAN_YEAR years}</li>
  267.      *   <li>time constant \(\tau\) of the exponential decay set to 60 {@link Constants#JULIAN_DAY days}</li>
  268.      *   <li>convergence set to 10⁻¹² s</li>
  269.      *   <li>polynomial part set to degree 3</li>
  270.      *   <li>one harmonic term at {@link #SUN_PULSATION}}</li>
  271.      *   <li>one harmonic term at 2 times {@link #SUN_PULSATION}}</li>
  272.      *   <li>one harmonic term at 3 times {@link #SUN_PULSATION}}</li>
  273.      *   <li>one harmonic term at {@link #MOON_DRACONIC_PULSATION}}</li>
  274.      *   <li>one harmonic term at 2 times {@link #MOON_DRACONIC_PULSATION}}</li>
  275.      *   <li>one harmonic term at 3 times {@link #MOON_DRACONIC_PULSATION}}</li>
  276.      * </ul>
  277.      * @return fitter with default configuration for pole nutation parameters dx and dy
  278.      */
  279.     public static SingleParameterFitter createDefaultNutationFitterLongTermPrediction() {
  280.         return new SingleParameterFitter(3 * Constants.JULIAN_YEAR, 60 * Constants.JULIAN_DAY, 1.0e-12, 3,
  281.                              SUN_PULSATION, 2 * SUN_PULSATION, 3 * SUN_PULSATION,
  282.                              MOON_DRACONIC_PULSATION, 2 * MOON_DRACONIC_PULSATION, 3 * MOON_DRACONIC_PULSATION);
  283.     }

  284. }