1   /* Copyright 2002-2024 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 java.io.Serializable;
20  
21  import org.hipparchus.CalculusFieldElement;
22  import org.hipparchus.geometry.euclidean.threed.FieldRotation;
23  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
24  import org.hipparchus.geometry.euclidean.threed.Rotation;
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.errors.OrekitException;
30  import org.orekit.errors.OrekitInternalError;
31  import org.orekit.time.AbsoluteDate;
32  import org.orekit.time.FieldAbsoluteDate;
33  import org.orekit.time.TimeVectorFunction;
34  
35  /** Celestial Intermediate Reference Frame.
36   * <p>This provider includes precession effects according to either the IAU 2006 precession
37   * (also known as Nicole Capitaines's P03 precession theory) and IAU 2000A_R06 nutation
38   * for IERS 2010 conventions or the IAU 2000A precession-nutation model for IERS 2003
39   * conventions. These models replaced the older IAU-76 precession (Lieske) and IAU-80
40   * theory of nutation (Wahr) which were used in the classical equinox-based paradigm.
41   * It <strong>must</strong> be used with the Earth Rotation Angle (REA) defined by
42   * Capitaine's model and <strong>not</strong> IAU-82 sidereal
43   * time which is consistent with the older models only.</p>
44   * <p>Its parent frame is the GCRF frame.
45   */
46  class CIRFProvider implements EOPBasedTransformProvider {
47  
48      /** Serializable UID. */
49      private static final long serialVersionUID = 20130806L;
50  
51      /** Function computing CIP/CIO components. */
52      private final transient TimeVectorFunction xysPxy2Function;
53  
54      /** EOP history. */
55      private final EOPHistory eopHistory;
56  
57      /** Simple constructor.
58       * @param eopHistory EOP history
59       * @see Frame
60       */
61      CIRFProvider(final EOPHistory eopHistory) {
62  
63          // load the nutation model
64          xysPxy2Function = eopHistory.getConventions()
65                  .getXYSpXY2Function(eopHistory.getTimeScales());
66  
67          // store correction to the model
68          this.eopHistory = eopHistory;
69  
70      }
71  
72      /** {@inheritDoc} */
73      @Override
74      public EOPHistory getEOPHistory() {
75          return eopHistory;
76      }
77  
78      /** {@inheritDoc} */
79      @Override
80      public CIRFProvider getNonInterpolatingProvider() {
81          return new CIRFProvider(eopHistory.getEOPHistoryWithoutCachedTidalCorrection());
82      }
83  
84      /** {@inheritDoc} */
85      @Override
86      public Transform getTransform(final AbsoluteDate date) {
87  
88          final double[] xys  = xysPxy2Function.value(date);
89          final double[] dxdy = eopHistory.getNonRotatinOriginNutationCorrection(date);
90  
91          // position of the Celestial Intermediate Pole (CIP)
92          final double xCurrent = xys[0] + dxdy[0];
93          final double yCurrent = xys[1] + dxdy[1];
94  
95          // position of the Celestial Intermediate Origin (CIO)
96          final double sCurrent = xys[2] - xCurrent * yCurrent / 2;
97  
98          // set up the bias, precession and nutation rotation
99          final double x2Py2  = xCurrent * xCurrent + yCurrent * yCurrent;
100         final double zP1    = 1 + FastMath.sqrt(1 - x2Py2);
101         final double r      = FastMath.sqrt(x2Py2);
102         final double sPe2   = 0.5 * (sCurrent + FastMath.atan2(yCurrent, xCurrent));
103         final SinCos sc     = FastMath.sinCos(sPe2);
104         final double xPr    = xCurrent + r;
105         final double xPrCos = xPr * sc.cos();
106         final double xPrSin = xPr * sc.sin();
107         final double yCos   = yCurrent * sc.cos();
108         final double ySin   = yCurrent * sc.sin();
109         final Rotation bpn  = new Rotation(zP1 * (xPrCos + ySin), -r * (yCos + xPrSin),
110                                            r * (xPrCos - ySin), zP1 * (yCos - xPrSin),
111                                            true);
112 
113         return new Transform(date, bpn, Vector3D.ZERO);
114 
115     }
116 
117     /** {@inheritDoc} */
118     @Override
119     public <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
120 
121         final T[] xys  = xysPxy2Function.value(date);
122         final T[] dxdy = eopHistory.getNonRotatinOriginNutationCorrection(date);
123 
124         // position of the Celestial Intermediate Pole (CIP)
125         final T xCurrent = xys[0].add(dxdy[0]);
126         final T yCurrent = xys[1].add(dxdy[1]);
127 
128         // position of the Celestial Intermediate Origin (CIO)
129         final T sCurrent = xys[2].subtract(xCurrent.multiply(yCurrent).multiply(0.5));
130 
131         // set up the bias, precession and nutation rotation
132         final T x2Py2           = xCurrent.multiply(xCurrent).add(yCurrent.multiply(yCurrent));
133         final T zP1             = x2Py2.subtract(1).negate().sqrt().add(1);
134         final T r               = x2Py2.sqrt();
135         final T sPe2            = sCurrent.add(yCurrent.atan2(xCurrent)).multiply(0.5);
136         final FieldSinCos<T> sc = FastMath.sinCos(sPe2);
137         final T xPr             = xCurrent.add(r);
138         final T xPrCos          = xPr.multiply(sc.cos());
139         final T xPrSin          = xPr.multiply(sc.sin());
140         final T yCos            = yCurrent.multiply(sc.cos());
141         final T ySin            = yCurrent.multiply(sc.sin());
142         final FieldRotation<T> bpn  = new FieldRotation<>(zP1.multiply(xPrCos.add(ySin)),
143                                                           r.multiply(yCos.add(xPrSin)).negate(),
144                                                           r.multiply(xPrCos.subtract(ySin)),
145                                                           zP1.multiply(yCos.subtract(xPrSin)),
146                                                           true);
147 
148         return new FieldTransform<>(date, bpn, FieldVector3D.getZero(date.getField()));
149 
150     }
151 
152     /** Replace the instance with a data transfer object for serialization.
153      * <p>
154      * This intermediate class serializes only the frame key.
155      * </p>
156      * @return data transfer object that will be serialized
157      */
158     private Object writeReplace() {
159         return new DataTransferObject(eopHistory);
160     }
161 
162     /** Internal class used only for serialization. */
163     private static class DataTransferObject implements Serializable {
164 
165         /** Serializable UID. */
166         private static final long serialVersionUID = 20131209L;
167 
168         /** EOP history. */
169         private final EOPHistory eopHistory;
170 
171         /** Simple constructor.
172          * @param eopHistory EOP history
173          */
174         DataTransferObject(final EOPHistory eopHistory) {
175             this.eopHistory = eopHistory;
176         }
177 
178         /** Replace the deserialized data transfer object with a {@link CIRFProvider}.
179          * @return replacement {@link CIRFProvider}
180          */
181         private Object readResolve() {
182             try {
183                 // retrieve a managed frame
184                 return new CIRFProvider(eopHistory);
185             } catch (OrekitException oe) {
186                 throw new OrekitInternalError(oe);
187             }
188         }
189 
190     }
191 
192 }