1   /* Copyright 2002-2024 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  
19  import java.io.Serializable;
20  import java.util.List;
21  import java.util.ListIterator;
22  import java.util.function.ToDoubleFunction;
23  
24  import org.hipparchus.util.FastMath;
25  import org.hipparchus.util.MathUtils;
26  import org.orekit.utils.Constants;
27  import org.orekit.utils.SecularAndHarmonic;
28  
29  /** Fitter for one Earth Orientation Parameter.
30   * @see PredictedEOPHistory
31   * @see EOPFitter
32   * @see SecularAndHarmonic
33   * @since 12.0
34   * @author Luc Maisonobe
35   */
36  public class SingleParameterFitter implements Serializable {
37  
38      /** Sun pulsation, one year period. */
39      public static final double SUN_PULSATION = MathUtils.TWO_PI / Constants.JULIAN_YEAR;
40  
41      /** Moon pulsation (one Moon draconic period). */
42      public static final double MOON_DRACONIC_PULSATION = MathUtils.TWO_PI / (27.212221 * Constants.JULIAN_DAY);
43  
44      /** Serializable UID. */
45      private static final long serialVersionUID = 20230309L;
46  
47      /** Time constant of the exponential decay weight. */
48      private final double timeConstant;
49  
50      /** Convergence on fitted parameter. */
51      private final double convergence;
52  
53      /** Degree of the polynomial model. */
54      private final int degree;
55  
56      /** Pulsations of harmonic part (rad/s). */
57      private final double[] pulsations;
58  
59      /** Simple constructor.
60       * @param fittingDuration ignored parameter since 12.0
61       * @param timeConstant time constant \(\tau\) of the exponential decay weight, point weight is \(e^{\frac{t-t_0}{\tau}}\),
62       * i.e. points far in the past before \(t_0\) have smaller weights
63       * @param convergence convergence on fitted parameter
64       * @param degree degree of the polynomial model
65       * @param pulsations pulsations of harmonic part (rad/s)
66       * @see #createDefaultDut1FitterShortTermPrediction()
67       * @see #createDefaultDut1FitterLongTermPrediction()
68       * @see #createDefaultPoleFitterShortTermPrediction()
69       * @see #createDefaultPoleFitterLongTermPrediction()
70       * @see #createDefaultNutationFitterShortTermPrediction()
71       * @see #createDefaultNutationFitterLongTermPrediction()
72       * @see SecularAndHarmonic
73       * @deprecated replaced by {@link #SingleParameterFitter(double, double, int, double...)}
74       */
75      @Deprecated
76      public SingleParameterFitter(final double fittingDuration, final double timeConstant, final double convergence,
77                                   final int degree, final double... pulsations) {
78          this(timeConstant, convergence, degree, pulsations);
79      }
80  
81      /** Simple constructor.
82       * @param timeConstant time constant \(\tau\) of the exponential decay weight, point weight is \(e^{\frac{t-t_0}{\tau}}\),
83       * i.e. points far in the past before \(t_0\) have smaller weights
84       * @param convergence convergence on fitted parameter
85       * @param degree degree of the polynomial model
86       * @param pulsations pulsations of harmonic part (rad/s)
87       * @see #createDefaultDut1FitterShortTermPrediction()
88       * @see #createDefaultDut1FitterLongTermPrediction()
89       * @see #createDefaultPoleFitterShortTermPrediction()
90       * @see #createDefaultPoleFitterLongTermPrediction()
91       * @see #createDefaultNutationFitterShortTermPrediction()
92       * @see #createDefaultNutationFitterLongTermPrediction()
93       * @see SecularAndHarmonic
94       * @since 12.0.1
95       */
96      public SingleParameterFitter(final double timeConstant, final double convergence,
97                                   final int degree, final double... pulsations) {
98          this.timeConstant    = timeConstant;
99          this.convergence     = convergence;
100         this.degree          = degree;
101         this.pulsations      = pulsations.clone();
102     }
103 
104     /** Perform secular and harmonic fitting.
105      * @param rawHistory EOP history to fit
106      * @param extractor extractor for Earth Orientation Parameter
107      * @return configured fitter
108      */
109     public SecularAndHarmonic fit(final EOPHistory rawHistory, final ToDoubleFunction<EOPEntry> extractor) {
110 
111         final List<EOPEntry> rawEntries = rawHistory.getEntries();
112         final EOPEntry       last       = rawEntries.get(rawEntries.size() - 1);
113 
114         // create fitter
115         final SecularAndHarmonic sh = new SecularAndHarmonic(degree, pulsations);
116 
117         // set up convergence
118         sh.setConvergenceRMS(convergence);
119 
120         // set up reference date and initial guess to a constant value
121         final double[] initialGuess = new double[degree + 1 + 2 * pulsations.length];
122         initialGuess[0] = extractor.applyAsDouble(last);
123         sh.resetFitting(last.getDate(), initialGuess);
124 
125         // sample history
126         final ListIterator<EOPEntry> backwardIterator = rawEntries.listIterator(rawEntries.size());
127         while (backwardIterator.hasPrevious()) {
128             final EOPEntry entry = backwardIterator.previous();
129             sh.addWeightedPoint(entry.getDate(), extractor.applyAsDouble(entry),
130                                 FastMath.exp(entry.getDate().durationFrom(last.getDate()) / timeConstant));
131         }
132 
133         // perform fitting
134         sh.fit();
135 
136         return sh;
137 
138     }
139 
140     /** Create fitter with default parameters adapted for fitting orientation parameters dUT1 and LOD
141      * for short term prediction.
142      * <p>
143      * The main difference between these settings and {@link #createDefaultDut1FitterLongTermPrediction()
144      * the settings for long prediction} is the much smaller \(\tau\). This means more
145      * weight is set to the points at the end of the history, hence forcing the fitted prediction
146      * model to be closer to these points, hence the prediction error to be smaller just after
147      * raw history end. On the other hand, this implies that the model will diverge on long term.
148      * These settings are intended when prediction is used for at most 5 days after raw EOP end.
149      * </p>
150      * <ul>
151      *   <li>time constant \(\tau\) of the exponential decay set to 6 {@link Constants#JULIAN_DAY days}</li>
152      *   <li>convergence set to 10⁻¹² s</li>
153      *   <li>polynomial part set to degree 3</li>
154      *   <li>one harmonic term at {@link #SUN_PULSATION}}</li>
155      *   <li>one harmonic term at 2 times {@link #SUN_PULSATION}}</li>
156      *   <li>one harmonic term at 3 times {@link #SUN_PULSATION}}</li>
157      *   <li>one harmonic term at {@link #MOON_DRACONIC_PULSATION}}</li>
158      *   <li>one harmonic term at 2 times {@link #MOON_DRACONIC_PULSATION}}</li>
159      *   <li>one harmonic term at 3 times {@link #MOON_DRACONIC_PULSATION}}</li>
160      * </ul>
161      * @return fitter with default configuration for orientation parameters dUT1 and LOD
162      * @see #createDefaultDut1FitterShortTermPrediction()
163      */
164     public static SingleParameterFitter createDefaultDut1FitterShortTermPrediction() {
165         return new SingleParameterFitter(6 * Constants.JULIAN_DAY, 1.0e-12, 3,
166                                          SUN_PULSATION, 2 * SUN_PULSATION, 3 * SUN_PULSATION,
167                                          MOON_DRACONIC_PULSATION, 2 * MOON_DRACONIC_PULSATION, 3 * MOON_DRACONIC_PULSATION);
168     }
169 
170     /** Create fitter with default parameters adapted for fitting orientation parameters dUT1 and LOD
171      * for long term prediction.
172      * <p>
173      * The main difference between these settings and {@link #createDefaultDut1FitterShortTermPrediction()
174      * the settings for short prediction} is the much larger \(\tau\). This means weight
175      * is spread throughout history, hence forcing the fitted prediction model to be remain very stable
176      * on the long term. On the other hand, this implies that the model will start with already a much
177      * larger error just after raw history end.
178      * These settings are intended when prediction is used for 5 days after raw EOP end or more.
179      * </p>
180      * <ul>
181      *   <li>time constant \(\tau\) of the exponential decay set to 60 {@link Constants#JULIAN_DAY days}</li>
182      *   <li>convergence set to 10⁻¹² s</li>
183      *   <li>polynomial part set to degree 3</li>
184      *   <li>one harmonic term at {@link #SUN_PULSATION}}</li>
185      *   <li>one harmonic term at 2 times {@link #SUN_PULSATION}}</li>
186      *   <li>one harmonic term at 3 times {@link #SUN_PULSATION}}</li>
187      *   <li>one harmonic term at {@link #MOON_DRACONIC_PULSATION}}</li>
188      *   <li>one harmonic term at 2 times {@link #MOON_DRACONIC_PULSATION}}</li>
189      *   <li>one harmonic term at 3 times {@link #MOON_DRACONIC_PULSATION}}</li>
190      * </ul>
191      * @return fitter with default configuration for orientation parameters dUT1 and LOD
192      * @see #createDefaultDut1FitterShortTermPrediction()
193      */
194     public static SingleParameterFitter createDefaultDut1FitterLongTermPrediction() {
195         return new SingleParameterFitter(60 * Constants.JULIAN_DAY, 1.0e-12, 3,
196                                          SUN_PULSATION, 2 * SUN_PULSATION, 3 * SUN_PULSATION,
197                                          MOON_DRACONIC_PULSATION, 2 * MOON_DRACONIC_PULSATION, 3 * MOON_DRACONIC_PULSATION);
198     }
199 
200     /** Create fitter with default parameters adapted for fitting pole parameters Xp and Yp
201      * for long term prediction.
202      * <p>
203      * The main difference between these settings and {@link #createDefaultPoleFitterLongTermPrediction()
204      * the settings for long prediction} is the much smaller \(\tau\). This means more
205      * weight is set to the points at the end of the history, hence forcing the fitted prediction
206      * model to be closer to these points, hence the prediction error to be smaller just after
207      * raw history end. On the other hand, this implies that the model will diverge on long term.
208      * These settings are intended when prediction is used for at most 5 days after raw EOP end.
209      * </p>
210      * <ul>
211      *   <li>time constant \(\tau\) of the exponential decay set to 12 {@link Constants#JULIAN_DAY days}</li>
212      *   <li>convergence set to 10⁻¹² rad</li>
213      *   <li>polynomial part set to degree 3</li>
214      *   <li>one harmonic term at {@link #SUN_PULSATION}}</li>
215      *   <li>one harmonic term at 2 times {@link #SUN_PULSATION}}</li>
216      *   <li>one harmonic term at 3 times {@link #SUN_PULSATION}}</li>
217      *   <li>one harmonic term at {@link #MOON_DRACONIC_PULSATION}}</li>
218      *   <li>one harmonic term at 2 times {@link #MOON_DRACONIC_PULSATION}}</li>
219      *   <li>one harmonic term at 3 times {@link #MOON_DRACONIC_PULSATION}}</li>
220      * </ul>
221      * @return fitter with default configuration for pole parameters Xp and Yp
222      */
223     public static SingleParameterFitter createDefaultPoleFitterShortTermPrediction() {
224         return new SingleParameterFitter(12 * Constants.JULIAN_DAY, 1.0e-12, 3,
225                                          SUN_PULSATION, 2 * SUN_PULSATION, 3 * SUN_PULSATION,
226                                          MOON_DRACONIC_PULSATION, 2 * MOON_DRACONIC_PULSATION, 3 * MOON_DRACONIC_PULSATION);
227     }
228 
229     /** Create fitter with default parameters adapted for fitting pole parameters Xp and Yp
230      * for long term prediction.
231      * <p>
232      * The main difference between these settings and {@link #createDefaultPoleFitterShortTermPrediction()
233      * the settings for short prediction} is the much larger \(\tau\). This means weight
234      * is spread throughout history, hence forcing the fitted prediction model to be remain very stable
235      * on the long term. On the other hand, this implies that the model will start with already a much
236      * larger error just after raw history end.
237      * These settings are intended when prediction is used for 5 days after raw EOP end or more.
238      * </p>
239      * <ul>
240      *   <li>time constant \(\tau\) of the exponential decay set to 60 {@link Constants#JULIAN_DAY days}</li>
241      *   <li>convergence set to 10⁻¹² rad</li>
242      *   <li>polynomial part set to degree 3</li>
243      *   <li>one harmonic term at {@link #SUN_PULSATION}}</li>
244      *   <li>one harmonic term at 2 times {@link #SUN_PULSATION}}</li>
245      *   <li>one harmonic term at 3 times {@link #SUN_PULSATION}}</li>
246      *   <li>one harmonic term at {@link #MOON_DRACONIC_PULSATION}}</li>
247      *   <li>one harmonic term at 2 times {@link #MOON_DRACONIC_PULSATION}}</li>
248      *   <li>one harmonic term at 3 times {@link #MOON_DRACONIC_PULSATION}}</li>
249      * </ul>
250      * @return fitter with default configuration for pole parameters Xp and Yp
251      */
252     public static SingleParameterFitter createDefaultPoleFitterLongTermPrediction() {
253         return new SingleParameterFitter(60 * Constants.JULIAN_DAY, 1.0e-12, 3,
254                                          SUN_PULSATION, 2 * SUN_PULSATION, 3 * SUN_PULSATION,
255                                          MOON_DRACONIC_PULSATION, 2 * MOON_DRACONIC_PULSATION, 3 * MOON_DRACONIC_PULSATION);
256     }
257 
258     /** Create fitter with default parameters adapted for fitting nutation parameters dx and dy
259      * for long term prediction.
260      * <p>
261      * The main difference between these settings and {@link #createDefaultNutationFitterLongTermPrediction()
262      * the settings for long prediction} is the much smaller \(\tau\). This means more
263      * weight is set to the points at the end of the history, hence forcing the fitted prediction
264      * model to be closer to these points, hence the prediction error to be smaller just after
265      * raw history end. On the other hand, this implies that the model will diverge on long term.
266      * These settings are intended when prediction is used for at most 5 days after raw EOP end.
267      * </p>
268      * <ul>
269      *   <li>time constant \(\tau\) of the exponential decay set to 12 {@link Constants#JULIAN_DAY days}</li>
270      *   <li>convergence set to 10⁻¹² s</li>
271      *   <li>polynomial part set to degree 3</li>
272      *   <li>one harmonic term at {@link #SUN_PULSATION}}</li>
273      *   <li>one harmonic term at 2 times {@link #SUN_PULSATION}}</li>
274      *   <li>one harmonic term at 3 times {@link #SUN_PULSATION}}</li>
275      *   <li>one harmonic term at {@link #MOON_DRACONIC_PULSATION}}</li>
276      *   <li>one harmonic term at 2 times {@link #MOON_DRACONIC_PULSATION}}</li>
277      *   <li>one harmonic term at 3 times {@link #MOON_DRACONIC_PULSATION}}</li>
278      * </ul>
279      * @return fitter with default configuration for pole nutation parameters dx and dy
280      */
281     public static SingleParameterFitter createDefaultNutationFitterShortTermPrediction() {
282         return new SingleParameterFitter(12 * Constants.JULIAN_DAY, 1.0e-12, 3,
283                                          SUN_PULSATION, 2 * SUN_PULSATION, 3 * SUN_PULSATION,
284                                          MOON_DRACONIC_PULSATION, 2 * MOON_DRACONIC_PULSATION, 3 * MOON_DRACONIC_PULSATION);
285     }
286 
287     /** Create fitter with default parameters adapted for fitting nutation parameters dx and dy
288      * for long term prediction.
289      * <p>
290      * The main difference between these settings and {@link #createDefaultNutationFitterShortTermPrediction()
291      * the settings for short prediction} is the much larger \(\tau\). This means weight
292      * is spread throughout history, hence forcing the fitted prediction model to be remain very stable
293      * on the long term. On the other hand, this implies that the model will start with already a much
294      * larger error just after raw history end.
295      * These settings are intended when prediction is used for 5 days after raw EOP end or more.
296      * </p>
297      * <ul>
298      *   <li>time constant \(\tau\) of the exponential decay set to 60 {@link Constants#JULIAN_DAY days}</li>
299      *   <li>convergence set to 10⁻¹² s</li>
300      *   <li>polynomial part set to degree 3</li>
301      *   <li>one harmonic term at {@link #SUN_PULSATION}}</li>
302      *   <li>one harmonic term at 2 times {@link #SUN_PULSATION}}</li>
303      *   <li>one harmonic term at 3 times {@link #SUN_PULSATION}}</li>
304      *   <li>one harmonic term at {@link #MOON_DRACONIC_PULSATION}}</li>
305      *   <li>one harmonic term at 2 times {@link #MOON_DRACONIC_PULSATION}}</li>
306      *   <li>one harmonic term at 3 times {@link #MOON_DRACONIC_PULSATION}}</li>
307      * </ul>
308      * @return fitter with default configuration for pole nutation parameters dx and dy
309      */
310     public static SingleParameterFitter createDefaultNutationFitterLongTermPrediction() {
311         return new SingleParameterFitter(60 * Constants.JULIAN_DAY, 1.0e-12, 3,
312                                          SUN_PULSATION, 2 * SUN_PULSATION, 3 * SUN_PULSATION,
313                                          MOON_DRACONIC_PULSATION, 2 * MOON_DRACONIC_PULSATION, 3 * MOON_DRACONIC_PULSATION);
314     }
315 
316 }