1   /* Copyright 2002-2026 CS GROUP
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.bodies;
18  
19  import java.io.Serial;
20  import java.util.ArrayList;
21  import java.util.List;
22  
23  import org.hipparchus.CalculusFieldElement;
24  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
25  import org.hipparchus.geometry.euclidean.threed.Vector3D;
26  import org.hipparchus.util.FastMath;
27  import org.hipparchus.util.FieldSinCos;
28  import org.hipparchus.util.SinCos;
29  import org.orekit.bodies.JPLEphemeridesLoader.EphemerisType;
30  import org.orekit.time.AbsoluteDate;
31  import org.orekit.time.FieldAbsoluteDate;
32  import org.orekit.time.TimeScales;
33  import org.orekit.utils.Constants;
34  
35  /** Abstract class for predefined IAU poles.
36   * <p>The pole models provided here come from the
37   * <a href="https://doi.org/10.1007/s10569-017-9805-5"> Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015 </a>
38   * (WGCCRE 2015).</p>
39   * <p>WGCCRE 2015 provides recommendations that define and relate the coordinate systems of solar system bodies to their rotational elements to support making
40   * cartographic products (i.e., “mapping”) of such bodies, including planets and their satellites, dwarf/minor planets and their satellites, and comets.
41   * Here predefined models for the Sun, Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, Neptune, Pluto, and the Moon, are provided.
42   * </p>
43   * @author Luc Maisonobe
44   * @since 9.0
45   */
46  public abstract class PredefinedIAUPoles implements IAUPole {
47  
48      /** Serializable UID. */
49      @Serial
50      private static final long serialVersionUID = 20200130L;
51  
52      /** Time scales. */
53      private final TimeScales timeScales;
54      /**
55       * {@code true} iff this is GCRF aligned.
56       *
57       * @since 14.0
58       */
59      private final boolean isGcrf;
60  
61      /**
62       * Simple constructor.
63       *
64       * @param timeScales to use when computing the pole, including TDB and J2000.0.
65       */
66      PredefinedIAUPoles(final TimeScales timeScales) {
67          this(timeScales, false);
68      }
69  
70      /**
71       * Simple constructor.
72       *
73       * @param timeScales to use when computing the pole, including TDB and
74       *                   J2000.0.
75       * @param isGcrf     {@code true} iff this is GCRF aligned.
76       * @since 14.0
77       */
78      PredefinedIAUPoles(final TimeScales timeScales,
79                         final boolean isGcrf) {
80          this.timeScales = timeScales;
81          this.isGcrf = isGcrf;
82      }
83  
84      /** IAU pole and prime meridian model for Sun. */
85      private static class Sun extends PredefinedIAUPoles {
86  
87          /** Serializable UID. */
88          @Serial
89          private static final long serialVersionUID = 20200130L;
90  
91          /** Constant term of the prime meridian. */
92          private static final double W0 = 84.176;
93  
94          /** Rate term of the prime meridian. */
95          private static final double W_DOT = 14.1844000;
96  
97          /** Fixed pole. */
98          private final Vector3D pole = new Vector3D(FastMath.toRadians(286.13),
99                                                     FastMath.toRadians(63.87));
100 
101         /**
102          * Simple constructor.
103          *
104          * @param timeScales to use when computing the pole, including TDB and J2000.0.
105          */
106         Sun(final TimeScales timeScales) {
107             super(timeScales);
108         }
109 
110         /** {@inheritDoc} */
111         public Vector3D getPole(final AbsoluteDate date) {
112             return pole;
113         }
114 
115         /** {@inheritDoc} */
116         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPole(final FieldAbsoluteDate<T> date) {
117             return new FieldVector3D<>(date.getField(), pole);
118         }
119 
120         /** {@inheritDoc} */
121         public double getPrimeMeridianAngle(final AbsoluteDate date) {
122             return FastMath.toRadians(d(date) * W_DOT + W0);
123         }
124 
125         /** {@inheritDoc} */
126         public <T extends CalculusFieldElement<T>> T getPrimeMeridianAngle(final FieldAbsoluteDate<T> date) {
127             return FastMath.toRadians(d(date).multiply(W_DOT).add(W0));
128         }
129 
130     }
131 
132     /** IAU pole and prime meridian model for Mercury.
133      * <p>The 20◦ meridian of Mercury is defined by the crater Hun Kal.</p>
134      */
135     private static class Mercury extends PredefinedIAUPoles {
136 
137         /** Serializable UID. */
138         @Serial
139         private static final long serialVersionUID = 20200130L;
140 
141         /** Constant term of the right ascension of the pole. */
142         private static final double ALPHA_0 = 281.0103;
143 
144         /** Rate term of the right ascension of the pole. */
145         private static final double ALPHA_DOT = -0.0328;
146 
147         /** Constant term of the declination of the pole. */
148         private static final double DELTA_0 = 61.4155;
149 
150         /** Rate term of the declination of the pole. */
151         private static final double DELTA_DOT = -0.0049;
152 
153         /** Constant term of the prime meridian with an uncertainty of +-0.0037. */
154         private static final double W_0 = 329.5988;
155 
156         /** Rate term of the prime meridian. */
157         private static final double W_DOT = 6.1385108;
158 
159         /** M1 coefficient of the prime meridian. */
160         private static final double M1_COEFF = 0.01067257;
161 
162         /** M2 coefficient of the prime meridian. */
163         private static final double M2_COEFF = -0.00112309;
164 
165         /** M3 coefficient of the prime meridian. */
166         private static final double M3_COEFF = -0.00011040;
167 
168         /** M4 coefficient of the prime meridian. */
169         private static final double M4_COEFF = -0.00002539;
170 
171         /** M5 coefficient of the prime meridian. */
172         private static final double M5_COEFF = -0.00000571;
173 
174         /** Constant term of the M1 angle. */
175         private static final double M1_0   = 174.7910857;
176 
177         /** Rate term of the M1 angle. */
178         private static final double M1_DOT = 4.092335;
179 
180         /** Constant term of the M2 angle. */
181         private static final double M2_0   = 349.5821714;
182 
183         /** Rate term of the M1 angle. */
184         private static final double M2_DOT = 8.184670;
185 
186         /** Constant term of the M3 angle. */
187         private static final double M3_0   = 164.3732571;
188 
189         /** Rate term of the M1 angle. */
190         private static final double M3_DOT = 12.277005;
191 
192         /** Constant term of the M4 angle. */
193         private static final double M4_0   = 339.1643429;
194 
195         /** Rate term of the M1 angle. */
196         private static final double M4_DOT = 16.369340;
197 
198         /** Constant term of the M5 angle. */
199         private static final double M5_0   = 153.9554286;
200 
201         /** Rate term of the M1 angle. */
202         private static final double M5_DOT = 20.461675;
203 
204         /**
205          * Simple constructor.
206          *
207          * @param timeScales to use when computing the pole, including TDB and J2000.0.
208          */
209         Mercury(final TimeScales timeScales) {
210             super(timeScales);
211         }
212 
213         /** {@inheritDoc} */
214         public Vector3D getPole(final AbsoluteDate date) {
215             final double t = t(date);
216             return new Vector3D(FastMath.toRadians(t * ALPHA_DOT + ALPHA_0),
217                                 FastMath.toRadians(t * DELTA_DOT + DELTA_0));
218         }
219 
220         /** {@inheritDoc} */
221         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPole(final FieldAbsoluteDate<T> date) {
222             final T t = t(date);
223             return new FieldVector3D<>(FastMath.toRadians(t.multiply(ALPHA_DOT).add(ALPHA_0)),
224                                        FastMath.toRadians(t.multiply(DELTA_DOT).add(DELTA_0)));
225         }
226 
227         /** {@inheritDoc} */
228         public double getPrimeMeridianAngle(final AbsoluteDate date) {
229             final double d = d(date);
230             return FastMath.toRadians(d(date) * W_DOT + W_0 +
231                                       FastMath.sin(FastMath.toRadians(d * M1_DOT + M1_0)) * M1_COEFF +
232                                       FastMath.sin(FastMath.toRadians(d * M2_DOT + M2_0)) * M2_COEFF +
233                                       FastMath.sin(FastMath.toRadians(d * M3_DOT + M3_0)) * M3_COEFF +
234                                       FastMath.sin(FastMath.toRadians(d * M4_DOT + M4_0)) * M4_COEFF +
235                                       FastMath.sin(FastMath.toRadians(d * M5_DOT + M5_0)) * M5_COEFF);
236         }
237 
238         /** {@inheritDoc} */
239         public <T extends CalculusFieldElement<T>> T getPrimeMeridianAngle(final FieldAbsoluteDate<T> date) {
240             final T d = d(date);
241             return FastMath.toRadians(d(date).multiply(W_DOT).add(W_0).
242                              add(FastMath.toRadians(d.multiply(M1_DOT).add(M1_0)).sin().multiply(M1_COEFF)).
243                              add(FastMath.toRadians(d.multiply(M2_DOT).add(M2_0)).sin().multiply(M2_COEFF)).
244                              add(FastMath.toRadians(d.multiply(M3_DOT).add(M3_0)).sin().multiply(M3_COEFF)).
245                              add(FastMath.toRadians(d.multiply(M4_DOT).add(M4_0)).sin().multiply(M4_COEFF)).
246                              add(FastMath.toRadians(d.multiply(M5_DOT).add(M5_0)).sin().multiply(M5_COEFF)));
247         }
248 
249     }
250 
251     /** IAU pole and prime meridian model for Venus.
252      * The 0◦ meridian of Venus is defined by the central peak in the crater Ariadne.
253      */
254     private static class Venus extends PredefinedIAUPoles {
255 
256         /** Serializable UID. */
257         @Serial
258         private static final long serialVersionUID = 20200130L;
259 
260         /** Constant term of the prime meridian. */
261         private static final double W_0 = 160.20;
262 
263         /** Rate term of the prime meridian. */
264         private static final double W_DOT = -1.4813688;
265 
266         /** Fixed pole. */
267         private final Vector3D pole = new Vector3D(FastMath.toRadians(272.76),
268                                                    FastMath.toRadians(67.16));
269 
270         /**
271          * Simple constructor.
272          *
273          * @param timeScales to use when computing the pole, including TDB and J2000.0.
274          */
275         Venus(final TimeScales timeScales) {
276             super(timeScales);
277         }
278 
279         /** {@inheritDoc} */
280         public Vector3D getPole(final AbsoluteDate date) {
281             return pole;
282         }
283 
284         /** {@inheritDoc} */
285         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPole(final FieldAbsoluteDate<T> date) {
286             return new FieldVector3D<>(date.getField(), pole);
287         }
288 
289         /** {@inheritDoc} */
290         public double getPrimeMeridianAngle(final AbsoluteDate date) {
291             return FastMath.toRadians(d(date) * W_DOT + W_0);
292         }
293 
294         /** {@inheritDoc} */
295         public <T extends CalculusFieldElement<T>> T getPrimeMeridianAngle(final FieldAbsoluteDate<T> date) {
296             return FastMath.toRadians(d(date).multiply(W_DOT).add(W_0));
297         }
298 
299     }
300 
301     /** Approximate expressions for the Earth
302      * have been removed in order to avoid confusion in WGCCRE 2015.
303      * <p>IAU pole and prime meridian model for Earth.</p>
304      */
305     private static class Earth extends PredefinedIAUPoles {
306 
307         /** Serializable UID. */
308         @Serial
309         private static final long serialVersionUID = 20200130L;
310 
311         /** Constant term of the right ascension of the pole. */
312         private static final double ALPHA_0 =  0.00;
313 
314         /** Rate term of the right ascension of the pole. */
315         private static final double ALPHA_DOT = -0.641;
316 
317         /** Constant term of the declination of the pole. */
318         private static final double DELTA_0 = 90.00;
319 
320         /** Rate term of the declination of the pole. */
321         private static final double DELTA_DOT = -0.557;
322 
323         /** Constant term of the prime meridian. */
324         private static final double W_0 = 190.147;
325 
326         /** Rate term of the prime meridian. */
327         private static final double W_DOT = 360.9856235;
328 
329         /**
330          * Simple constructor.
331          *
332          * @param timeScales to use when computing the pole, including TDB and J2000.0.
333          */
334         Earth(final TimeScales timeScales) {
335             super(timeScales);
336         }
337 
338         /** {@inheritDoc} */
339         public Vector3D getPole(final AbsoluteDate date) {
340             final double t = t(date);
341             return new Vector3D(FastMath.toRadians(t * ALPHA_DOT + ALPHA_0),
342                                 FastMath.toRadians(t * DELTA_DOT + DELTA_0));
343         }
344 
345         /** {@inheritDoc} */
346         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPole(final FieldAbsoluteDate<T> date) {
347             final T t = t(date);
348             return new FieldVector3D<>(FastMath.toRadians(t.multiply(ALPHA_DOT).add(ALPHA_0)),
349                                        FastMath.toRadians(t.multiply(DELTA_DOT).add(DELTA_0)));
350         }
351 
352         /** {@inheritDoc} */
353         @Override
354         public Vector3D getNode(final AbsoluteDate date) {
355             final double t = t(date);
356             return new Vector3D(FastMath.toRadians(t * ALPHA_DOT + ALPHA_0 + 90.0),
357                                 0.0);
358         }
359 
360         /** {@inheritDoc} */
361         @Override
362         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getNode(final FieldAbsoluteDate<T> date) {
363             final T t = t(date);
364             return new FieldVector3D<>(FastMath.toRadians(t.multiply(ALPHA_DOT).add(ALPHA_0 + 90.0)),
365                                        date.getField().getZero());
366         }
367 
368         /** {@inheritDoc} */
369         public double getPrimeMeridianAngle(final AbsoluteDate date) {
370             return FastMath.toRadians(d(date) * W_DOT + W_0);
371         }
372 
373         /** {@inheritDoc} */
374         public <T extends CalculusFieldElement<T>> T getPrimeMeridianAngle(final FieldAbsoluteDate<T> date) {
375             return FastMath.toRadians(d(date).multiply(W_DOT).add(W_0));
376         }
377 
378     }
379 
380     /** Due to low precision of following series expression,
381      * WGCCRE 2015 has removed this model for the Moon’s orientation.
382      * For high accuracy work (e.g., spacecraft operations, high-resolution mapping, and gravity
383      * field determination), it is recommended that a lunar ephemeris be used to obtain the libration
384      * angles for the Moon, from which the pole position and rotation can be derived.
385      * <p>IAU pole and prime meridian model for the Moon.</p>
386      */
387     private static class Moon extends PredefinedIAUPoles {
388 
389         /** Serializable UID. */
390         @Serial
391         private static final long serialVersionUID = 20200130L;
392 
393         /** Constant term of the right ascension of the pole. */
394         private static final double ALPHA_0 = 269.9949;
395 
396         /** Rate term of the right ascension of the pole. */
397         private static final double ALPHA_DOT = 0.0031;
398 
399         /** Constant term of the declination of the pole. */
400         private static final double DELTA_0 = 66.5392;
401 
402         /** Rate term of the declination of the pole. */
403         private static final double DELTA_DOT =  0.0130;
404 
405         /** Constant term of the prime meridian. */
406         private static final double W_0 = 38.3213;
407 
408         /** Rate term of the prime meridian. */
409         private static final double W_DOT = 13.17635815;
410 
411         /** Rate term of the prime meridian. */
412         private static final double W_DOT_DOT = -1.4e-12;
413 
414         /** Constant term of the E1 angle. */
415         private static final double E01_0    = 125.045;
416 
417         /** Rate term of the E1 angle. */
418         private static final double E01_DOT  =  -0.0529921;
419 
420         /** Sine coefficient of the E1 angle. */
421         private static final double E01_SIN  = -3.8787;
422 
423         /** Cosine coefficient of the E1 angle. */
424         private static final double E01_COS  =  1.5419;
425 
426         /** Sine coefficient of the E1 angle, for the prime meridian. */
427         private static final double E01_W_SIN =  3.5610;
428 
429         /** Constant term of the E2 angle. */
430         private static final double E02_0    = 250.089;
431 
432         /** Rate term of the E2 angle. */
433         private static final double E02_DOT  =  -0.1059842;
434 
435         /** Sine coefficient of the E2 angle. */
436         private static final double E02_SIN  = -0.1204;
437 
438         /** Cosine coefficient of the E2 angle. */
439         private static final double E02_COS  =  0.0239;
440 
441         /** Sine coefficient of the E2 angle, for the prime meridian. */
442         private static final double E02_W_SIN =  0.1208;
443 
444         /** Constant term of the E3 angle. */
445         private static final double E03_0    = 260.008;
446 
447         /** Rate term of the E3 angle. */
448         private static final double E03_DOT  =  13.0120009;
449 
450         /** Sine coefficient of the E3 angle. */
451         private static final double E03_SIN  =  0.0700;
452 
453         /** Cosine coefficient of the E3 angle. */
454         private static final double E03_COS  = -0.0278;
455 
456         /** Sine coefficient of the E3 angle, for the prime meridian. */
457         private static final double E03_W_SIN = -0.0642;
458 
459         /** Constant term of the E4 angle. */
460         private static final double E04_0    = 176.625;
461 
462         /** Rate term of the E4 angle. */
463         private static final double E04_DOT  =  13.3407154;
464 
465         /** Sine coefficient of the E4 angle. */
466         private static final double E04_SIN  = -0.0172;
467 
468         /** Cosine coefficient of the E4 angle. */
469         private static final double E04_COS  =  0.0068;
470 
471         /** Sine coefficient of the E4 angle, for the prime meridian. */
472         private static final double E04_W_SIN =  0.0158;
473 
474         /** Constant term of the E5 angle. */
475         private static final double E05_0    = 357.529;
476 
477         /** Rate term of the E5 angle. */
478         private static final double E05_DOT  =   0.9856003;
479 
480         /** Sine coefficient of the E5 angle, for the prime meridian. */
481         private static final double E05_W_SIN =  0.0252;
482 
483         /** Constant term of the E6 angle. */
484         private static final double E06_0    = 311.589;
485 
486         /** Rate term of the E6 angle. */
487         private static final double E06_DOT  =  26.4057084;
488 
489         /** Sine coefficient of the E6 angle. */
490         private static final double E06_SIN  = 0.0072;
491 
492         /** Cosine coefficient of the E6 angle. */
493         private static final double E06_COS  = -0.0029;
494 
495         /** Sine coefficient of the E6 angle, for the prime meridian. */
496         private static final double E06_W_SIN = -0.0066;
497 
498         /** Constant term of the E7 angle. */
499         private static final double E07_0    = 134.963;
500 
501         /** Rate term of the E7 angle. */
502         private static final double E07_DOT  =  13.0649930;
503 
504         /** Cosine coefficient of the E7 angle. */
505         private static final double E07_COS  =  0.0009;
506 
507         /** Sine coefficient of the E7 angle, for the prime meridian. */
508         private static final double E07_W_SIN = -0.0047;
509 
510         /** Constant term of the E8 angle. */
511         private static final double E08_0    = 276.617;
512 
513         /** Rate term of the E8 angle. */
514         private static final double E08_DOT  =   0.3287146;
515 
516         /** Sine coefficient of the E8 angle, for the prime meridian. */
517         private static final double E08_W_SIN = -0.0046;
518 
519         /** Constant term of the E9 angle. */
520         private static final double E09_0    =  34.226;
521 
522         /** Rate term of the E9 angle. */
523         private static final double E09_DOT  =   1.7484877;
524 
525         /** Sine coefficient of the E9 angle, for the prime meridian. */
526         private static final double E09_W_SIN =  0.0028;
527 
528         /** Constant term of the E10 angle. */
529         private static final double E10_0    =  15.134;
530 
531         /** Rate term of the E10 angle. */
532         private static final double E10_DOT  =  -0.1589763;
533 
534         /** Sine coefficient of the E10 angle. */
535         private static final double E10_SIN  = -0.0052;
536 
537         /** Cosine coefficient of the E10 angle. */
538         private static final double E10_COS  = 0.0008;
539 
540         /** Sine coefficient of the E10 angle, for the prime meridian. */
541         private static final double E10_W_SIN =  0.0052;
542 
543         /** Constant term of the E11 angle. */
544         private static final double E11_0    = 119.743;
545 
546         /** Rate term of the E11 angle. */
547         private static final double E11_DOT  =   0.0036096;
548 
549         /** Sine coefficient of the E11 angle, for the prime meridian. */
550         private static final double E11_W_SIN =  0.0040;
551 
552         /** Constant term of the E12 angle. */
553         private static final double E12_0    = 239.961;
554 
555         /** Rate term of the E12 angle. */
556         private static final double E12_DOT  =   0.1643573;
557 
558         /** Sine coefficient of the E12 angle, for the prime meridian. */
559         private static final double E12_W_SIN =  0.0019;
560 
561         /** Constant term of the E13 angle. */
562         private static final double E13_0    =  25.053;
563 
564         /** Rate term of the E13 angle. */
565         private static final double E13_DOT  =  12.9590088;
566 
567         /** Sine coefficient of the E13 angle. */
568         private static final double E13_SIN  = 0.0043;
569 
570         /** Cosine coefficient of the E13 angle. */
571         private static final double E13_COS  = -0.0009;
572 
573         /** Sine coefficient of the E13 angle, for the prime meridian. */
574         private static final double E13_W_SIN = -0.0044;
575 
576         /**
577          * Simple constructor.
578          *
579          * @param timeScales to use when computing the pole, including TDB and J2000.0.
580          */
581         Moon(final TimeScales timeScales) {
582             super(timeScales);
583         }
584 
585         /** {@inheritDoc} */
586         public Vector3D getPole(final AbsoluteDate date) {
587             final double d = d(date);
588             final double t = t(date);
589 
590             final SinCos scE01 = FastMath.sinCos(FastMath.toRadians(d * E01_DOT + E01_0));
591             final SinCos scE02 = FastMath.sinCos(FastMath.toRadians(d * E02_DOT + E02_0));
592             final SinCos scE03 = FastMath.sinCos(FastMath.toRadians(d * E03_DOT + E03_0));
593             final SinCos scE04 = FastMath.sinCos(FastMath.toRadians(d * E04_DOT + E04_0));
594             final SinCos scE06 = FastMath.sinCos(FastMath.toRadians(d * E06_DOT + E06_0));
595             final SinCos scE10 = FastMath.sinCos(FastMath.toRadians(d * E10_DOT + E10_0));
596             final SinCos scE13 = FastMath.sinCos(FastMath.toRadians(d * E13_DOT + E13_0));
597 
598             return new Vector3D(FastMath.toRadians(t * ALPHA_DOT + ALPHA_0 +
599                                                    scE01.sin() * E01_SIN +
600                                                    scE02.sin() * E02_SIN +
601                                                    scE03.sin() * E03_SIN +
602                                                    scE04.sin() * E04_SIN +
603                                                    scE06.sin() * E06_SIN +
604                                                    scE10.sin() * E10_SIN +
605                                                    scE13.sin() * E13_SIN),
606                                 FastMath.toRadians(t * DELTA_DOT + DELTA_0 +
607                                                    scE01.cos() * E01_COS +
608                                                    scE02.cos() * E02_COS +
609                                                    scE03.cos() * E03_COS +
610                                                    scE04.cos() * E04_COS +
611                                                    scE06.cos() * E06_COS +
612                                                    FastMath.cos(FastMath.toRadians(d * E07_DOT + E07_0)) * E07_COS +  // only the cosine is needed
613                                                    scE10.cos() * E10_COS +
614                                                    scE13.cos() * E13_COS));
615         }
616 
617         /** {@inheritDoc} */
618         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPole(final FieldAbsoluteDate<T> date) {
619             final T d = d(date);
620             final T t = t(date);
621 
622             final FieldSinCos<T> scE01 = FastMath.sinCos(FastMath.toRadians(d.multiply(E01_DOT).add(E01_0)));
623             final FieldSinCos<T> scE02 = FastMath.sinCos(FastMath.toRadians(d.multiply(E02_DOT).add(E02_0)));
624             final FieldSinCos<T> scE03 = FastMath.sinCos(FastMath.toRadians(d.multiply(E03_DOT).add(E03_0)));
625             final FieldSinCos<T> scE04 = FastMath.sinCos(FastMath.toRadians(d.multiply(E04_DOT).add(E04_0)));
626             final FieldSinCos<T> scE06 = FastMath.sinCos(FastMath.toRadians(d.multiply(E06_DOT).add(E06_0)));
627             final FieldSinCos<T> scE10 = FastMath.sinCos(FastMath.toRadians(d.multiply(E10_DOT).add(E10_0)));
628             final FieldSinCos<T> scE13 = FastMath.sinCos(FastMath.toRadians(d.multiply(E13_DOT).add(E13_0)));
629 
630             return new FieldVector3D<>(FastMath.toRadians(t.multiply(ALPHA_DOT).add(ALPHA_0).
631                                                  add(scE01.sin().multiply(E01_SIN)).
632                                                  add(scE02.sin().multiply(E02_SIN)).
633                                                  add(scE03.sin().multiply(E03_SIN)).
634                                                  add(scE04.sin().multiply(E04_SIN)).
635                                                  add(scE06.sin().multiply(E06_SIN)).
636                                                  add(scE10.sin().multiply(E10_SIN)).
637                                                  add(scE13.sin().multiply(E13_SIN))),
638                                        FastMath.toRadians(t.multiply(DELTA_DOT).add(DELTA_0).
639                                                  add(scE01.cos().multiply(E01_COS)).
640                                                  add(scE02.cos().multiply(E02_COS)).
641                                                  add(scE03.cos().multiply(E03_COS)).
642                                                  add(scE04.cos().multiply(E04_COS)).
643                                                  add(scE06.cos().multiply(E06_COS)).
644                                                  add(FastMath.toRadians(d.multiply(E07_DOT).add(E07_0)).cos().multiply(E07_COS)).// only the cosine is needed
645                                                  add(scE10.cos().multiply(E10_COS)).
646                                                  add(scE13.cos().multiply(E13_COS))));
647         }
648 
649         /** {@inheritDoc} */
650         public double getPrimeMeridianAngle(final AbsoluteDate date) {
651             final double d = d(date);
652 
653             return FastMath.toRadians(d * (d * W_DOT_DOT + W_DOT) + W_0 +
654                                       FastMath.sin(FastMath.toRadians(d * E01_DOT + E01_0)) * E01_W_SIN +
655                                       FastMath.sin(FastMath.toRadians(d * E02_DOT + E02_0)) * E02_W_SIN +
656                                       FastMath.sin(FastMath.toRadians(d * E03_DOT + E03_0)) * E03_W_SIN +
657                                       FastMath.sin(FastMath.toRadians(d * E04_DOT + E04_0)) * E04_W_SIN +
658                                       FastMath.sin(FastMath.toRadians(d * E05_DOT + E05_0)) * E05_W_SIN +
659                                       FastMath.sin(FastMath.toRadians(d * E06_DOT + E06_0)) * E06_W_SIN +
660                                       FastMath.sin(FastMath.toRadians(d * E07_DOT + E07_0)) * E07_W_SIN +
661                                       FastMath.sin(FastMath.toRadians(d * E08_DOT + E08_0)) * E08_W_SIN +
662                                       FastMath.sin(FastMath.toRadians(d * E09_DOT + E09_0)) * E09_W_SIN +
663                                       FastMath.sin(FastMath.toRadians(d * E10_DOT + E10_0)) * E10_W_SIN +
664                                       FastMath.sin(FastMath.toRadians(d * E11_DOT + E11_0)) * E11_W_SIN +
665                                       FastMath.sin(FastMath.toRadians(d * E12_DOT + E12_0)) * E12_W_SIN +
666                                       FastMath.sin(FastMath.toRadians(d * E13_DOT + E13_0)) * E13_W_SIN);
667         }
668 
669         /** {@inheritDoc} */
670         public <T extends CalculusFieldElement<T>> T getPrimeMeridianAngle(final FieldAbsoluteDate<T> date) {
671             final T d = d(date);
672             return FastMath.toRadians(d.multiply(d.multiply(W_DOT_DOT).add(W_DOT)).add(W_0).
673                                       add(FastMath.toRadians(d.multiply(E01_DOT).add(E01_0)).sin().multiply(E01_W_SIN)).
674                                       add(FastMath.toRadians(d.multiply(E02_DOT).add(E02_0)).sin().multiply(E02_W_SIN)).
675                                       add(FastMath.toRadians(d.multiply(E03_DOT).add(E03_0)).sin().multiply(E03_W_SIN)).
676                                       add(FastMath.toRadians(d.multiply(E04_DOT).add(E04_0)).sin().multiply(E04_W_SIN)).
677                                       add(FastMath.toRadians(d.multiply(E05_DOT).add(E05_0)).sin().multiply(E05_W_SIN)).
678                                       add(FastMath.toRadians(d.multiply(E06_DOT).add(E06_0)).sin().multiply(E06_W_SIN)).
679                                       add(FastMath.toRadians(d.multiply(E07_DOT).add(E07_0)).sin().multiply(E07_W_SIN)).
680                                       add(FastMath.toRadians(d.multiply(E08_DOT).add(E08_0)).sin().multiply(E08_W_SIN)).
681                                       add(FastMath.toRadians(d.multiply(E09_DOT).add(E09_0)).sin().multiply(E09_W_SIN)).
682                                       add(FastMath.toRadians(d.multiply(E10_DOT).add(E10_0)).sin().multiply(E10_W_SIN)).
683                                       add(FastMath.toRadians(d.multiply(E11_DOT).add(E11_0)).sin().multiply(E11_W_SIN)).
684                                       add(FastMath.toRadians(d.multiply(E12_DOT).add(E12_0)).sin().multiply(E12_W_SIN)).
685                                       add(FastMath.toRadians(d.multiply(E13_DOT).add(E13_0)).sin().multiply(E13_W_SIN)));
686         }
687 
688     }
689 
690     /** IAU pole and prime meridian model for Mars.
691      * <p>The longitude of the Viking 1 lander on Mars is defined to be 47◦.95137 west,
692      * maintaining the 0◦ meridian through the crater Airy-0.</p>
693      */
694     private static class Mars extends PredefinedIAUPoles {
695 
696         /** Serializable UID. */
697         @Serial
698         private static final long serialVersionUID = 20200130L;
699 
700         /** Constant term of the right ascension of the pole. */
701         private static final double ALPHA_0 = 317.269202;
702 
703         /** Rate term of the right ascension of the pole. */
704         private static final double ALPHA_DOT = -0.10927547;
705 
706         /** Alpha1 constant term. */
707         private static final double ALPHA1_0 = 198.991226;
708 
709         /** Alpha1 rate term. */
710         private static final double ALPHA1_DOT = 19139.4819985;
711 
712         /** Alpha1 sine coefficient. */
713         private static final double ALPHA1_SIN = 0.000068;
714 
715         /** Alpha2 constant term. */
716         private static final double ALPHA2_0 = 226.292679;
717 
718         /** Alpha2 rate term. */
719         private static final double ALPHA2_DOT = 38280.8511281;
720 
721         /** Alpha2 sine coefficient. */
722         private static final double ALPHA2_SIN = 0.000238;
723 
724         /** Alpha3 constant term. */
725         private static final double ALPHA3_0 = 249.663391;
726 
727         /** Alpha3 rate term. */
728         private static final double ALPHA3_DOT = 57420.7251593;
729 
730         /** Alpha3 sine coefficient. */
731         private static final double ALPHA3_SIN = 0.000052;
732 
733         /** Alpha4 constant term. */
734         private static final double ALPHA4_0 = 266.183510;
735 
736         /** Alpha4 rate term. */
737         private static final double ALPHA4_DOT = 76560.6367950;
738 
739         /** Alpha4 sine coefficient. */
740         private static final double ALPHA4_SIN = 0.000009;
741 
742         /** Alpha5 constant term. */
743         private static final double ALPHA5_0 = 79.398797;
744 
745         /** Alpha5 rate term. */
746         private static final double ALPHA5_DOT = 0.5042615;
747 
748         /** Alpha5 sine coefficient. */
749         private static final double ALPHA5_SIN = 0.419057;
750 
751         /** Constant term of the declination of the pole. */
752         private static final double DELTA_0 = 54.432516;
753 
754         /** Rate term of the declination of the pole. */
755         private static final double DELTA_DOT = -0.05827105;
756 
757         /** Delta1 constant term. */
758         private static final double DELTA1_0 = 122.433576;
759 
760         /** Delta1 rate term. */
761         private static final double DELTA1_DOT = 19139.9407476;
762 
763         /** Delta1 cosine coefficient. */
764         private static final double DELTA1_COS = 0.000051;
765 
766         /** Delta2 constant term. */
767         private static final double DELTA2_0 = 43.058401;
768 
769         /** Delta2 rate term. */
770         private static final double DELTA2_DOT = 38280.8753272;
771 
772         /** Delta2 cosine coefficient. */
773         private static final double DELTA2_COS = 0.000141;
774 
775         /** Delta3 constant term. */
776         private static final double DELTA3_0 = 57.663379;
777 
778         /** Delta3 rate term. */
779         private static final double DELTA3_DOT = 57420.7517205;
780 
781         /** Delta3 cosine coefficient. */
782         private static final double DELTA3_COS = 0.000031;
783 
784         /** Delta4 constant term. */
785         private static final double DELTA4_0 = 79.476401;
786 
787         /** Delta4 rate term. */
788         private static final double DELTA4_DOT = 76560.6495004;
789 
790         /** Delta4 cosine coefficient. */
791         private static final double DELTA4_COS = 0.000005;
792 
793         /** Delta5 constant term. */
794         private static final double DELTA5_0 = 166.325722;
795 
796         /** Delta5 rate term. */
797         private static final double DELTA5_DOT = 0.5042615;
798 
799         /** Delta5 cosine coefficient. */
800         private static final double DELTA5_COS = 1.591274;
801 
802         /** Constant term of the prime meridian. */
803         private static final double W_0 = 176.049863;
804 
805         /** Rate term of the prime meridian. */
806         private static final double W_DOT = 350.891982443297;
807 
808         /** W1 constant term. */
809         private static final double W1_0 = 129.071773;
810 
811         /** W1 rate term. */
812         private static final double W1_DOT = 19140.0328244;
813 
814         /** W1 sine coefficient. */
815         private static final double W1_SIN = 0.000145;
816 
817         /** W2 constant term. */
818         private static final double W2_0 = 36.352167;
819 
820         /** W2 rate term. */
821         private static final double W2_DOT = 38281.0473591;
822 
823         /** W2 sine coefficient. */
824         private static final double W2_SIN = 0.000157;
825 
826         /** W3 constant term. */
827         private static final double W3_0 = 56.668646;
828 
829         /** W3 rate term. */
830         private static final double W3_DOT = 57420.9295360;
831 
832         /** W3 sine coefficient. */
833         private static final double W3_SIN = 0.000040;
834 
835         /** W4 constant term. */
836         private static final double W4_0 = 67.364003;
837 
838         /** W4 rate term. */
839         private static final double W4_DOT = 76560.2552215;
840 
841         /** W4 sine coefficient. */
842         private static final double W4_SIN = 0.000001;
843 
844         /** W5 constant term. */
845         private static final double W5_0 = 104.792680;
846 
847         /** W5 rate term. */
848         private static final double W5_DOT = 95700.4387578;
849 
850         /** W5 sine coefficient. */
851         private static final double W5_SIN = 0.000001;
852 
853         /** W6 constant term. */
854         private static final double W6_0 = 95.391654;
855 
856         /** W6 rate term. */
857         private static final double W6_DOT = 0.5042615;
858 
859         /** W6 sine coefficient. */
860         private static final double W6_SIN = 0.584542;
861 
862         /**
863          * Simple constructor.
864          *
865          * @param timeScales to use when computing the pole, including TDB and J2000.0.
866          */
867         Mars(final TimeScales timeScales) {
868             super(timeScales);
869         }
870 
871         /** {@inheritDoc} */
872         public Vector3D getPole(final AbsoluteDate date) {
873             final double t = t(date);
874 
875             final double alpha = t * ALPHA_DOT + ALPHA_0 +
876                 FastMath.sin(FastMath.toRadians(t * ALPHA1_DOT + ALPHA1_0)) * ALPHA1_SIN +
877                 FastMath.sin(FastMath.toRadians(t * ALPHA2_DOT + ALPHA2_0)) * ALPHA2_SIN +
878                 FastMath.sin(FastMath.toRadians(t * ALPHA3_DOT + ALPHA3_0)) * ALPHA3_SIN +
879                 FastMath.sin(FastMath.toRadians(t * ALPHA4_DOT + ALPHA4_0)) * ALPHA4_SIN +
880                 FastMath.sin(FastMath.toRadians(t * ALPHA5_DOT + ALPHA5_0)) * ALPHA5_SIN;
881 
882             final double delta = t * DELTA_DOT + DELTA_0 +
883                 FastMath.cos(FastMath.toRadians(t * DELTA1_DOT + DELTA1_0)) * DELTA1_COS +
884                 FastMath.cos(FastMath.toRadians(t * DELTA2_DOT + DELTA2_0)) * DELTA2_COS +
885                 FastMath.cos(FastMath.toRadians(t * DELTA3_DOT + DELTA3_0)) * DELTA3_COS +
886                 FastMath.cos(FastMath.toRadians(t * DELTA4_DOT + DELTA4_0)) * DELTA4_COS +
887                 FastMath.cos(FastMath.toRadians(t * DELTA5_DOT + DELTA5_0)) * DELTA5_COS;
888 
889             return new Vector3D(FastMath.toRadians(alpha), FastMath.toRadians(delta));
890         }
891 
892         /** {@inheritDoc} */
893         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPole(final FieldAbsoluteDate<T> date) {
894             final T t = t(date);
895 
896             final T alpha = t.multiply(ALPHA_DOT).add(ALPHA_0)
897                 .add(FastMath.toRadians(t.multiply(ALPHA1_DOT).add(ALPHA1_0)).sin().multiply(ALPHA1_SIN))
898                 .add(FastMath.toRadians(t.multiply(ALPHA2_DOT).add(ALPHA2_0)).sin().multiply(ALPHA2_SIN))
899                 .add(FastMath.toRadians(t.multiply(ALPHA3_DOT).add(ALPHA3_0)).sin().multiply(ALPHA3_SIN))
900                 .add(FastMath.toRadians(t.multiply(ALPHA4_DOT).add(ALPHA4_0)).sin().multiply(ALPHA4_SIN))
901                 .add(FastMath.toRadians(t.multiply(ALPHA5_DOT).add(ALPHA5_0)).sin().multiply(ALPHA5_SIN));
902 
903             final T delta = t.multiply(DELTA_DOT).add(DELTA_0)
904                 .add(FastMath.toRadians(t.multiply(DELTA1_DOT).add(DELTA1_0)).cos().multiply(DELTA1_COS))
905                 .add(FastMath.toRadians(t.multiply(DELTA2_DOT).add(DELTA2_0)).cos().multiply(DELTA2_COS))
906                 .add(FastMath.toRadians(t.multiply(DELTA3_DOT).add(DELTA3_0)).cos().multiply(DELTA3_COS))
907                 .add(FastMath.toRadians(t.multiply(DELTA4_DOT).add(DELTA4_0)).cos().multiply(DELTA4_COS))
908                 .add(FastMath.toRadians(t.multiply(DELTA5_DOT).add(DELTA5_0)).cos().multiply(DELTA5_COS));
909 
910             return new FieldVector3D<>(FastMath.toRadians(alpha), FastMath.toRadians(delta));
911         }
912 
913         /** {@inheritDoc} */
914         public double getPrimeMeridianAngle(final AbsoluteDate date) {
915             final double d = d(date);
916             final double t = t(date);
917 
918             final double w = d * W_DOT + W_0 +
919                 FastMath.sin(FastMath.toRadians(t * W1_DOT + W1_0)) * W1_SIN +
920                 FastMath.sin(FastMath.toRadians(t * W2_DOT + W2_0)) * W2_SIN +
921                 FastMath.sin(FastMath.toRadians(t * W3_DOT + W3_0)) * W3_SIN +
922                 FastMath.sin(FastMath.toRadians(t * W4_DOT + W4_0)) * W4_SIN +
923                 FastMath.sin(FastMath.toRadians(t * W5_DOT + W5_0)) * W5_SIN +
924                 FastMath.sin(FastMath.toRadians(t * W6_DOT + W6_0)) * W6_SIN;
925 
926             return FastMath.toRadians(w);
927         }
928 
929         /** {@inheritDoc} */
930         public <T extends CalculusFieldElement<T>> T getPrimeMeridianAngle(final FieldAbsoluteDate<T> date) {
931             final T d = d(date);
932             final T t = t(date);
933 
934             final T w = d.multiply(W_DOT).add(W_0)
935                 .add(FastMath.toRadians(t.multiply(W1_DOT).add(W1_0)).sin().multiply(W1_SIN))
936                 .add(FastMath.toRadians(t.multiply(W2_DOT).add(W2_0)).sin().multiply(W2_SIN))
937                 .add(FastMath.toRadians(t.multiply(W3_DOT).add(W3_0)).sin().multiply(W3_SIN))
938                 .add(FastMath.toRadians(t.multiply(W4_DOT).add(W4_0)).sin().multiply(W4_SIN))
939                 .add(FastMath.toRadians(t.multiply(W5_DOT).add(W5_0)).sin().multiply(W5_SIN))
940                 .add(FastMath.toRadians(t.multiply(W6_DOT).add(W6_0)).sin().multiply(W6_SIN));
941 
942             return FastMath.toRadians(w);
943         }
944 
945     }
946 
947     /** IAU pole and prime meridian model for Jupiter. */
948     private static class Jupiter extends PredefinedIAUPoles {
949 
950         /** Serializable UID. */
951         @Serial
952         private static final long serialVersionUID = 20200130L;
953 
954         /** Constant term of the right ascension of the pole. */
955         private static final double ALPHA_0 = 268.056595;
956 
957         /** Rate term of the right ascension of the pole. */
958         private static final double ALPHA_DOT = -0.006499;
959 
960         /** Constant term of the declination of the pole. */
961         private static final double DELTA_0 = 64.495303;
962 
963         /** Rate term of the declination of the pole. */
964         private static final double DELTA_DOT = 0.002413;
965 
966         /** Constant term of the ja angle. */
967         private static final double JA_0 =  99.360714;
968 
969         /** Rate term of the ja angle. */
970         private static final double JA_DOT = 4850.4046;
971 
972         /** Sine coefficient of the ja angle. */
973         private static final double JA_SIN = 0.000117;
974 
975         /** Cosine coefficient of the ja angle. */
976         private static final double JA_COS = 0.000050;
977 
978         /** Constant term of the jb angle. */
979         private static final double JB_0 = 175.895369;
980 
981         /** Rate term of the jb angle. */
982         private static final double JB_DOT = 1191.9605;
983 
984         /** Sine coefficient of the jb angle. */
985         private static final double JB_SIN = 0.000938;
986 
987         /** Cosine coefficient of the jb angle. */
988         private static final double JB_COS = 0.000404;
989 
990         /** Constant term of the jc angle. */
991         private static final double JC_0 = 300.323162;
992 
993         /** Rate term of the jc angle. */
994         private static final double JC_DOT = 262.5475;
995 
996         /** Sine coefficient of the jc angle. */
997         private static final double JC_SIN = 0.001432;
998 
999         /** Cosine coefficient of the jc angle. */
1000         private static final double JC_COS = 0.000617;
1001 
1002         /** Constant term of the jd angle. */
1003         private static final double JD_0 = 114.012305;
1004 
1005         /** Rate term of the jd angle. */
1006         private static final double JD_DOT = 6070.2476;
1007 
1008         /** Sine coefficient of the jd angle. */
1009         private static final double JD_SIN = 0.000030;
1010 
1011         /** Cosine coefficient of the jd angle. */
1012         private static final double JD_COS = -0.000013;
1013 
1014         /** Constant term of the je angle. */
1015         private static final double JE_0 = 49.511251;
1016 
1017         /** Rate term of the je angle. */
1018         private static final double JE_DOT = 64.3000;
1019 
1020         /** Sine coefficient of the je angle. */
1021         private static final double JE_SIN = 0.002150;
1022 
1023         /** Cosine coefficient of the je angle. */
1024         private static final double JE_COS = 0.000926;
1025 
1026         /** Constant term of the prime meridian. */
1027         private static final double W_0 = 284.95;
1028 
1029         /** Rate term of the prime meridian. */
1030         private static final double W_DOT = 870.5360000;
1031 
1032         /**
1033          * Simple constructor.
1034          *
1035          * @param timeScales to use when computing the pole, including TDB and J2000.0.
1036          */
1037         Jupiter(final TimeScales timeScales) {
1038             super(timeScales);
1039         }
1040 
1041         /** {@inheritDoc} */
1042         public Vector3D getPole(final AbsoluteDate date) {
1043 
1044             final double t = t(date);
1045             final double ja = FastMath.toRadians(t * JA_DOT + JA_0);
1046             final double jb = FastMath.toRadians(t * JB_DOT + JB_0);
1047             final double jc = FastMath.toRadians(t * JC_DOT + JC_0);
1048             final double jd = FastMath.toRadians(t * JD_DOT + JD_0);
1049             final double je = FastMath.toRadians(t * JE_DOT + JE_0);
1050 
1051             final SinCos scJa = FastMath.sinCos(ja);
1052             final SinCos scJb = FastMath.sinCos(jb);
1053             final SinCos scJc = FastMath.sinCos(jc);
1054             final SinCos scJd = FastMath.sinCos(jd);
1055             final SinCos scJe = FastMath.sinCos(je);
1056 
1057             return new Vector3D(FastMath.toRadians(t * ALPHA_DOT + ALPHA_0 +
1058                                                    scJa.sin() * JA_SIN +
1059                                                    scJb.sin() * JB_SIN +
1060                                                    scJc.sin() * JC_SIN +
1061                                                    scJd.sin() * JD_SIN +
1062                                                    scJe.sin() * JE_SIN),
1063                                 FastMath.toRadians(t * DELTA_DOT + DELTA_0 +
1064                                                    scJa.cos() * JA_COS +
1065                                                    scJb.cos() * JB_COS +
1066                                                    scJc.cos() * JC_COS +
1067                                                    scJd.cos() * JD_COS +
1068                                                    scJe.cos() * JE_COS));
1069         }
1070 
1071         /** {@inheritDoc} */
1072         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPole(final FieldAbsoluteDate<T> date) {
1073 
1074             final T t = t(date);
1075             final T ja = FastMath.toRadians(t.multiply(JA_DOT).add(JA_0));
1076             final T jb = FastMath.toRadians(t.multiply(JB_DOT).add(JB_0));
1077             final T jc = FastMath.toRadians(t.multiply(JC_DOT).add(JC_0));
1078             final T jd = FastMath.toRadians(t.multiply(JD_DOT).add(JD_0));
1079             final T je = FastMath.toRadians(t.multiply(JE_DOT).add(JE_0));
1080 
1081             final FieldSinCos<T> scJa = FastMath.sinCos(ja);
1082             final FieldSinCos<T> scJb = FastMath.sinCos(jb);
1083             final FieldSinCos<T> scJc = FastMath.sinCos(jc);
1084             final FieldSinCos<T> scJd = FastMath.sinCos(jd);
1085             final FieldSinCos<T> scJe = FastMath.sinCos(je);
1086 
1087             return new FieldVector3D<>(FastMath.toRadians(t.multiply(ALPHA_DOT).add(ALPHA_0).
1088                                                  add(scJa.sin().multiply(JA_SIN)).
1089                                                  add(scJb.sin().multiply(JB_SIN)).
1090                                                  add(scJc.sin().multiply(JC_SIN)).
1091                                                  add(scJd.sin().multiply(JD_SIN)).
1092                                                  add(scJe.sin().multiply(JE_SIN))),
1093                                        FastMath.toRadians(t.multiply(DELTA_DOT).add(DELTA_0).
1094                                                  add(scJa.cos().multiply(JA_COS)).
1095                                                  add(scJb.cos().multiply(JB_COS)).
1096                                                  add(scJc.cos().multiply(JC_COS)).
1097                                                  add(scJd.cos().multiply(JD_COS)).
1098                                                  add(scJe.cos().multiply(JE_COS))));
1099 
1100         }
1101 
1102         /** {@inheritDoc} */
1103         public double getPrimeMeridianAngle(final AbsoluteDate date) {
1104             return FastMath.toRadians(d(date) * W_DOT + W_0);
1105         }
1106 
1107         /** {@inheritDoc} */
1108         public <T extends CalculusFieldElement<T>> T getPrimeMeridianAngle(final FieldAbsoluteDate<T> date) {
1109             return FastMath.toRadians(d(date).multiply(W_DOT).add(W_0));
1110         }
1111 
1112     }
1113 
1114     /** IAU pole and prime meridian model for Saturn. */
1115     private static class Saturn extends PredefinedIAUPoles {
1116 
1117         /** Serializable UID. */
1118         @Serial
1119         private static final long serialVersionUID = 20200130L;
1120 
1121         /** Constant term of the right ascension of the pole. */
1122         private static final double ALPHA_0 = 40.589;
1123 
1124         /** Rate term of the right ascension of the pole. */
1125         private static final double ALPHA_DOT = -0.036;
1126 
1127         /** Constant term of the declination of the pole. */
1128         private static final double DELTA_0 = 83.537;
1129 
1130         /** Rate term of the declination of the pole. */
1131         private static final double DELTA_DOT = -0.004;
1132 
1133         /** Constant term of the prime meridian. */
1134         private static final double W_0 = 38.90;
1135 
1136         /** Rate term of the prime meridian. */
1137         private static final double W_DOT = 810.7939024;
1138 
1139         /**
1140          * Simple constructor.
1141          *
1142          * @param timeScales to use when computing the pole, including TDB and J2000.0.
1143          */
1144         Saturn(final TimeScales timeScales) {
1145             super(timeScales);
1146         }
1147 
1148         /** {@inheritDoc} */
1149         public Vector3D getPole(final AbsoluteDate date) {
1150             final double t = t(date);
1151             return new Vector3D(FastMath.toRadians(t * ALPHA_DOT + ALPHA_0),
1152                                 FastMath.toRadians(t * DELTA_DOT + DELTA_0));
1153         }
1154 
1155         /** {@inheritDoc} */
1156         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPole(final FieldAbsoluteDate<T> date) {
1157             final T t = t(date);
1158             return new FieldVector3D<>(FastMath.toRadians(t.multiply(ALPHA_DOT).add(ALPHA_0)),
1159                                        FastMath.toRadians(t.multiply(DELTA_DOT).add(DELTA_0)));
1160         }
1161 
1162         /** {@inheritDoc} */
1163         public double getPrimeMeridianAngle(final AbsoluteDate date) {
1164             return FastMath.toRadians(d(date) * W_DOT + W_0);
1165         }
1166 
1167         /** {@inheritDoc} */
1168         public <T extends CalculusFieldElement<T>> T getPrimeMeridianAngle(final FieldAbsoluteDate<T> date) {
1169             return FastMath.toRadians(d(date).multiply(W_DOT).add(W_0));
1170         }
1171 
1172     }
1173 
1174     /** IAU pole and prime meridian model for Uranus. */
1175     private static class Uranus extends PredefinedIAUPoles {
1176 
1177         /** Serializable UID. */
1178         @Serial
1179         private static final long serialVersionUID = 20200130L;
1180 
1181         /** Constant term of the prime meridian. */
1182         private static final double W_0 = 203.81;
1183 
1184         /** Rate term of the prime meridian. */
1185         private static final double W_DOT = -501.1600928;
1186 
1187         /** Fixed pole. */
1188         private final Vector3D pole = new Vector3D(FastMath.toRadians(257.311),
1189                                                    FastMath.toRadians(-15.175));
1190 
1191         /**
1192          * Simple constructor.
1193          *
1194          * @param timeScales to use when computing the pole, including TDB and J2000.0.
1195          */
1196         Uranus(final TimeScales timeScales) {
1197             super(timeScales);
1198         }
1199 
1200         /** {@inheritDoc} */
1201         public Vector3D getPole(final AbsoluteDate date) {
1202             return pole;
1203         }
1204 
1205         /** {@inheritDoc} */
1206         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPole(final FieldAbsoluteDate<T> date) {
1207             return new FieldVector3D<>(date.getField(), pole);
1208         }
1209 
1210         /** {@inheritDoc} */
1211         public double getPrimeMeridianAngle(final AbsoluteDate date) {
1212             return FastMath.toRadians(d(date) * W_DOT + W_0);
1213         }
1214 
1215         /** {@inheritDoc} */
1216         public <T extends CalculusFieldElement<T>> T getPrimeMeridianAngle(final FieldAbsoluteDate<T> date) {
1217             return FastMath.toRadians(d(date).multiply(W_DOT).add(W_0));
1218         }
1219 
1220     }
1221 
1222     /** IAU pole and prime meridian model for Neptune. */
1223     private static class Neptune extends PredefinedIAUPoles {
1224 
1225         /** Serializable UID. */
1226         @Serial
1227         private static final long serialVersionUID = 20200130L;
1228 
1229         /** Constant term of the right ascension of the pole. */
1230         private static final double ALPHA_0 = 299.36;
1231 
1232         /** Sine term of the right ascension of the pole. */
1233         private static final double ALPHA_SIN = 0.70;
1234 
1235         /** Constant term of the declination of the pole. */
1236         private static final double DELTA_0 = 43.46;
1237 
1238         /** Cosine term of the declination of the pole. */
1239         private static final double DELTA_COS = -0.51;
1240 
1241         /** Constant term of the prime meridian. */
1242         private static final double W_0 = 249.978;
1243 
1244         /** Rate term of the prime meridian. */
1245         private static final double W_DOT = 541.1397757;
1246 
1247         /** Sine term of the prime meridian. */
1248         private static final double W_SIN = -0.48;
1249 
1250         /** Constant term of the N angle. */
1251         private static final double N_0   = 357.85;
1252 
1253         /** Rate term of the M1 angle. */
1254         private static final double N_DOT = 52.316;
1255 
1256         /**
1257          * Simple constructor.
1258          *
1259          * @param timeScales to use when computing the pole, including TDB and J2000.0.
1260          */
1261         Neptune(final TimeScales timeScales) {
1262             super(timeScales);
1263         }
1264 
1265         /** {@inheritDoc} */
1266         public Vector3D getPole(final AbsoluteDate date) {
1267             final double n  = FastMath.toRadians(t(date) * N_DOT + N_0);
1268             final SinCos sc = FastMath.sinCos(n);
1269             return new Vector3D(FastMath.toRadians(sc.sin() * ALPHA_SIN + ALPHA_0),
1270                                 FastMath.toRadians(sc.cos() * DELTA_COS + DELTA_0));
1271         }
1272 
1273         /** {@inheritDoc} */
1274         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPole(final FieldAbsoluteDate<T> date) {
1275             final T n = FastMath.toRadians(t(date).multiply(N_DOT).add(N_0));
1276             final FieldSinCos<T> sc = FastMath.sinCos(n);
1277             return new FieldVector3D<>(FastMath.toRadians(sc.sin().multiply(ALPHA_SIN).add(ALPHA_0)),
1278                                        FastMath.toRadians(sc.cos().multiply(DELTA_COS).add(DELTA_0)));
1279         }
1280 
1281         /** {@inheritDoc} */
1282         public double getPrimeMeridianAngle(final AbsoluteDate date) {
1283             final double n = FastMath.toRadians(t(date) * N_DOT + N_0);
1284             return FastMath.toRadians(d(date) * W_DOT + FastMath.sin(n) * W_SIN + W_0);
1285         }
1286 
1287         /** {@inheritDoc} */
1288         public <T extends CalculusFieldElement<T>> T getPrimeMeridianAngle(final FieldAbsoluteDate<T> date) {
1289             final T n = FastMath.toRadians(t(date).multiply(N_DOT).add(N_0));
1290             return FastMath.toRadians(d(date).multiply(W_DOT).add(n.sin().multiply(W_SIN)).add(W_0));
1291         }
1292 
1293     }
1294 
1295     /** IAU pole and prime meridian model for Pluto.
1296      * <p>The 0 meridian for Pluto is defined as the mean sub-Charon meridian.</p>
1297      */
1298     private static class Pluto extends PredefinedIAUPoles {
1299 
1300         /** Serializable UID. */
1301         @Serial
1302         private static final long serialVersionUID = 20200130L;
1303 
1304         /** Constant term of the prime meridian. */
1305         private static final double W_0 = 302.695;
1306 
1307         /** Rate term of the prime meridian. */
1308         private static final double W_DOT = 56.3625225;
1309 
1310         /** Fixed pole. */
1311         private final Vector3D pole = new Vector3D(FastMath.toRadians(132.993),
1312                                                    FastMath.toRadians(-6.163));
1313 
1314         /**
1315          * Simple constructor.
1316          *
1317          * @param timeScales to use when computing the pole, including TDB and J2000.0.
1318          */
1319         Pluto(final TimeScales timeScales) {
1320             super(timeScales);
1321         }
1322 
1323         /** {@inheritDoc} */
1324         public Vector3D getPole(final AbsoluteDate date) {
1325             return pole;
1326         }
1327 
1328         /** {@inheritDoc} */
1329         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPole(final FieldAbsoluteDate<T> date) {
1330             return new FieldVector3D<>(date.getField(), pole);
1331         }
1332 
1333         /** {@inheritDoc} */
1334         public double getPrimeMeridianAngle(final AbsoluteDate date) {
1335             return FastMath.toRadians(d(date) * W_DOT + W_0);
1336         }
1337 
1338         /** {@inheritDoc} */
1339         public <T extends CalculusFieldElement<T>> T getPrimeMeridianAngle(final FieldAbsoluteDate<T> date) {
1340             return FastMath.toRadians(d(date).multiply(W_DOT).add(W_0));
1341         }
1342 
1343     }
1344 
1345     /** Default IAUPole implementation for barycenters.
1346      * <p>This implementation defines directions such that the inertially oriented and body
1347      * oriented frames are identical and aligned with GCRF. It is used for example
1348      * to define the ICRF.</p>
1349      */
1350     private static class GcrfAligned extends PredefinedIAUPoles {
1351 
1352         /** Serializable UID. */
1353         @Serial
1354         private static final long serialVersionUID = 20200130L;
1355 
1356         /**
1357          * Simple constructor.
1358          *
1359          * @param timeScales to use when computing the pole, including TDB and J2000.0.
1360          */
1361         GcrfAligned(final TimeScales timeScales) {
1362             super(timeScales, true);
1363         }
1364 
1365         /** {@inheritDoc} */
1366         public Vector3D getPole(final AbsoluteDate date) {
1367             return Vector3D.PLUS_K;
1368         }
1369 
1370         /** {@inheritDoc} */
1371         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPole(final FieldAbsoluteDate<T> date) {
1372             return FieldVector3D.getPlusK(date.getField());
1373         }
1374 
1375         /** {@inheritDoc} */
1376         @Override
1377         public Vector3D getNode(final AbsoluteDate date) {
1378             return Vector3D.PLUS_I;
1379         }
1380 
1381         /** {@inheritDoc} */
1382         @Override
1383         public <T extends CalculusFieldElement<T>> FieldVector3D<T> getNode(final FieldAbsoluteDate<T> date) {
1384             return FieldVector3D.getPlusI(date.getField());
1385         }
1386 
1387         /** {@inheritDoc} */
1388         public double getPrimeMeridianAngle(final AbsoluteDate date) {
1389             return 0;
1390         }
1391 
1392         /** {@inheritDoc} */
1393         public <T extends CalculusFieldElement<T>> T getPrimeMeridianAngle(final FieldAbsoluteDate<T> date) {
1394             return date.getField().getZero();
1395         }
1396 
1397     }
1398 
1399 
1400     /** Get a predefined IAU pole.
1401      * @param body body identifier
1402      * @param timeScales to use when computing the pole, including TDB and J2000.0.
1403      * @return predefined IAU pole
1404      */
1405     public static PredefinedIAUPoles getIAUPole(final EphemerisType body,
1406                                                 final TimeScales timeScales) {
1407 
1408         return switch (body) {
1409             case SUN  -> new Sun(timeScales);
1410             case MERCURY  -> new Mercury(timeScales);
1411             case VENUS  -> new Venus(timeScales);
1412             case EARTH  -> new Earth(timeScales);
1413             case MOON  -> new Moon(timeScales);
1414             case MARS  -> new Mars(timeScales);
1415             case JUPITER  -> new Jupiter(timeScales);
1416             case SATURN  -> new Saturn(timeScales);
1417             case URANUS  -> new Uranus(timeScales);
1418             case NEPTUNE  -> new Neptune(timeScales);
1419             case PLUTO  -> new Pluto(timeScales);
1420             default  -> new GcrfAligned(timeScales);
1421         };
1422     }
1423 
1424     /**
1425      * List of predefined IAU poles.
1426      *
1427      * @param timeScales to use when computing the pole, including TDB and J2000.0.
1428      * @return the poles.
1429      */
1430     static List<PredefinedIAUPoles> values(final TimeScales timeScales) {
1431         final List<PredefinedIAUPoles> values = new ArrayList<>(12);
1432         values.add(new Sun(timeScales));
1433         values.add(new Mercury(timeScales));
1434         values.add(new Venus(timeScales));
1435         values.add(new Earth(timeScales));
1436         values.add(new Moon(timeScales));
1437         values.add(new Mars(timeScales));
1438         values.add(new Jupiter(timeScales));
1439         values.add(new Saturn(timeScales));
1440         values.add(new Uranus(timeScales));
1441         values.add(new Neptune(timeScales));
1442         values.add(new Pluto(timeScales));
1443         values.add(new GcrfAligned(timeScales));
1444         return values;
1445     }
1446 
1447     /** Compute the interval in julian centuries from standard epoch.
1448      * @param date date
1449      * @return interval between date and standard epoch in julian centuries
1450      */
1451     protected double t(final AbsoluteDate date) {
1452         return date.offsetFrom(timeScales.getJ2000Epoch(), timeScales.getTDB()) /
1453                 Constants.JULIAN_CENTURY;
1454     }
1455 
1456     /** Compute the interval in julian centuries from standard epoch.
1457      * @param date date
1458      * @param <T> type of the filed elements
1459      * @return interval between date and standard epoch in julian centuries
1460      */
1461     protected <T extends CalculusFieldElement<T>> T t(final FieldAbsoluteDate<T> date) {
1462         final FieldAbsoluteDate<T> j2000Epoch =
1463                 new FieldAbsoluteDate<>(date.getField(), timeScales.getJ2000Epoch());
1464         return date.offsetFrom(j2000Epoch, timeScales.getTDB()).divide(Constants.JULIAN_CENTURY);
1465     }
1466 
1467     /** Compute the interval in julian days from standard epoch.
1468      * @param date date
1469      * @return interval between date and standard epoch in julian days
1470      */
1471     protected double d(final AbsoluteDate date) {
1472         return date.offsetFrom(timeScales.getJ2000Epoch(), timeScales.getTDB()) /
1473                 Constants.JULIAN_DAY;
1474     }
1475 
1476     /** Compute the interval in julian days from standard epoch.
1477      * @param date date
1478      * @param <T> type of the filed elements
1479      * @return interval between date and standard epoch in julian days
1480      */
1481     protected <T extends CalculusFieldElement<T>> T d(final FieldAbsoluteDate<T> date) {
1482         final FieldAbsoluteDate<T> j2000Epoch =
1483                 new FieldAbsoluteDate<>(date.getField(), timeScales.getJ2000Epoch());
1484         return date.offsetFrom(j2000Epoch, timeScales.getTDB()).divide(Constants.JULIAN_DAY);
1485     }
1486 
1487     @Override
1488     public boolean isGcrfAligned() {
1489         return isGcrf;
1490     }
1491 
1492 }