1   /* Copyright 2002-2025 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.frames;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.geometry.euclidean.threed.FieldRotation;
21  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
22  import org.hipparchus.geometry.euclidean.threed.Rotation;
23  import org.hipparchus.geometry.euclidean.threed.Vector3D;
24  import org.hipparchus.util.FastMath;
25  import org.hipparchus.util.FieldSinCos;
26  import org.hipparchus.util.SinCos;
27  import org.orekit.time.AbsoluteDate;
28  import org.orekit.time.FieldAbsoluteDate;
29  import org.orekit.time.TimeVectorFunction;
30  
31  /** Celestial Intermediate Reference Frame.
32   * <p>This provider includes precession effects according to either the IAU 2006 precession
33   * (also known as Nicole Capitaines's P03 precession theory) and IAU 2000A_R06 nutation
34   * for IERS 2010 conventions or the IAU 2000A precession-nutation model for IERS 2003
35   * conventions. These models replaced the older IAU-76 precession (Lieske) and IAU-80
36   * theory of nutation (Wahr) which were used in the classical equinox-based paradigm.
37   * It <strong>must</strong> be used with the Earth Rotation Angle (REA) defined by
38   * Capitaine's model and <strong>not</strong> IAU-82 sidereal
39   * time which is consistent with the older models only.</p>
40   * <p>Its parent frame is the GCRF frame.
41   */
42  class CIRFProvider implements EOPBasedTransformProvider {
43  
44      /** Function computing CIP/CIO components. */
45      private final TimeVectorFunction xysPxy2Function;
46  
47      /** EOP history. */
48      private final EOPHistory eopHistory;
49  
50      /** Simple constructor.
51       * @param eopHistory EOP history
52       * @see Frame
53       */
54      CIRFProvider(final EOPHistory eopHistory) {
55  
56          // load the nutation model
57          xysPxy2Function = eopHistory.getConventions()
58                  .getXYSpXY2Function(eopHistory.getTimeScales());
59  
60          // store correction to the model
61          this.eopHistory = eopHistory;
62  
63      }
64  
65      /** {@inheritDoc} */
66      @Override
67      public EOPHistory getEOPHistory() {
68          return eopHistory;
69      }
70  
71      /** {@inheritDoc} */
72      @Override
73      public CIRFProvider getNonInterpolatingProvider() {
74          return new CIRFProvider(eopHistory.getEOPHistoryWithoutCachedTidalCorrection());
75      }
76  
77      /** {@inheritDoc} */
78      @Override
79      public Transform getTransform(final AbsoluteDate date) {
80  
81          final double[] xys  = xysPxy2Function.value(date);
82          final double[] dxdy = eopHistory.getNonRotatinOriginNutationCorrection(date);
83  
84          // position of the Celestial Intermediate Pole (CIP)
85          final double xCurrent = xys[0] + dxdy[0];
86          final double yCurrent = xys[1] + dxdy[1];
87  
88          // position of the Celestial Intermediate Origin (CIO)
89          final double sCurrent = xys[2] - xCurrent * yCurrent / 2;
90  
91          // set up the bias, precession and nutation rotation
92          final double x2Py2  = xCurrent * xCurrent + yCurrent * yCurrent;
93          final double zP1    = 1 + FastMath.sqrt(1 - x2Py2);
94          final double r      = FastMath.sqrt(x2Py2);
95          final double sPe2   = 0.5 * (sCurrent + FastMath.atan2(yCurrent, xCurrent));
96          final SinCos sc     = FastMath.sinCos(sPe2);
97          final double xPr    = xCurrent + r;
98          final double xPrCos = xPr * sc.cos();
99          final double xPrSin = xPr * sc.sin();
100         final double yCos   = yCurrent * sc.cos();
101         final double ySin   = yCurrent * sc.sin();
102         final Rotation bpn  = new Rotation(zP1 * (xPrCos + ySin), -r * (yCos + xPrSin),
103                                            r * (xPrCos - ySin), zP1 * (yCos - xPrSin),
104                                            true);
105 
106         return new Transform(date, bpn, Vector3D.ZERO);
107 
108     }
109 
110     /** {@inheritDoc} */
111     @Override
112     public <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
113 
114         final T[] xys  = xysPxy2Function.value(date);
115         final T[] dxdy = eopHistory.getNonRotatinOriginNutationCorrection(date);
116 
117         // position of the Celestial Intermediate Pole (CIP)
118         final T xCurrent = xys[0].add(dxdy[0]);
119         final T yCurrent = xys[1].add(dxdy[1]);
120 
121         // position of the Celestial Intermediate Origin (CIO)
122         final T sCurrent = xys[2].subtract(xCurrent.multiply(yCurrent).multiply(0.5));
123 
124         // set up the bias, precession and nutation rotation
125         final T x2Py2           = xCurrent.multiply(xCurrent).add(yCurrent.multiply(yCurrent));
126         final T zP1             = x2Py2.subtract(1).negate().sqrt().add(1);
127         final T r               = x2Py2.sqrt();
128         final T sPe2            = sCurrent.add(yCurrent.atan2(xCurrent)).multiply(0.5);
129         final FieldSinCos<T> sc = FastMath.sinCos(sPe2);
130         final T xPr             = xCurrent.add(r);
131         final T xPrCos          = xPr.multiply(sc.cos());
132         final T xPrSin          = xPr.multiply(sc.sin());
133         final T yCos            = yCurrent.multiply(sc.cos());
134         final T ySin            = yCurrent.multiply(sc.sin());
135         final FieldRotation<T> bpn  = new FieldRotation<>(zP1.multiply(xPrCos.add(ySin)),
136                                                           r.multiply(yCos.add(xPrSin)).negate(),
137                                                           r.multiply(xPrCos.subtract(ySin)),
138                                                           zP1.multiply(yCos.subtract(xPrSin)),
139                                                           true);
140 
141         return new FieldTransform<>(date, bpn, FieldVector3D.getZero(date.getField()));
142 
143     }
144 
145 }