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.bodies;
18  
19  
20  import org.hipparchus.CalculusFieldElement;
21  import org.hipparchus.Field;
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.RotationConvention;
26  import org.hipparchus.geometry.euclidean.threed.Vector3D;
27  import org.hipparchus.util.Precision;
28  import org.orekit.frames.FieldStaticTransform;
29  import org.orekit.frames.FieldTransform;
30  import org.orekit.frames.Frame;
31  import org.orekit.frames.StaticTransform;
32  import org.orekit.frames.Transform;
33  import org.orekit.frames.TransformProvider;
34  import org.orekit.time.AbsoluteDate;
35  import org.orekit.time.FieldAbsoluteDate;
36  import org.orekit.utils.FieldPVCoordinates;
37  import org.orekit.utils.PVCoordinates;
38  import org.orekit.utils.TimeStampedFieldPVCoordinates;
39  import org.orekit.utils.TimeStampedPVCoordinates;
40  
41  /** Implementation of the {@link CelestialBody} interface using JPL or INPOP ephemerides.
42   * @author Luc Maisonobe
43   */
44  class JPLCelestialBody implements CelestialBody {
45  
46      /** Name of the body. */
47      private final String name;
48  
49      /** Regular expression for supported files names. */
50      private final String supportedNames;
51  
52      /** Ephemeris type to generate. */
53      private final JPLEphemeridesLoader.EphemerisType generateType;
54  
55      /** Raw position-velocity provider. */
56      private final transient JPLEphemeridesLoader.RawPVProvider rawPVProvider;
57  
58      /** Attraction coefficient of the body (m³/s²). */
59      private final double gm;
60  
61      /** Scaling factor for position-velocity. */
62      private final double scale;
63  
64      /** IAU pole. */
65      private final IAUPole iauPole;
66  
67      /** Inertially oriented, body-centered frame. */
68      private final Frame inertialFrame;
69  
70      /** Body oriented, body-centered frame. */
71      private final Frame bodyFrame;
72  
73      /** Build an instance and the underlying frame.
74       * @param name name of the body
75       * @param supportedNames regular expression for supported files names
76       * @param generateType ephemeris type to generate
77       * @param rawPVProvider raw position-velocity provider
78       * @param gm attraction coefficient (in m³/s²)
79       * @param scale scaling factor for position-velocity
80       * @param iauPole IAU pole implementation
81       * @param definingFrameAlignedWithICRF frame in which celestial body coordinates are defined,
82       * this frame <strong>must</strong> be aligned with ICRF
83       * @param inertialFrameName name to use for inertial frame (if null a default name will be built)
84       * @param bodyOrientedFrameName name to use for body-oriented frame (if null a default name will be built)
85       */
86      JPLCelestialBody(final String name, final String supportedNames,
87                       final JPLEphemeridesLoader.EphemerisType generateType,
88                       final JPLEphemeridesLoader.RawPVProvider rawPVProvider,
89                       final double gm, final double scale,
90                       final IAUPole iauPole, final Frame definingFrameAlignedWithICRF,
91                       final String inertialFrameName, final String bodyOrientedFrameName) {
92          this.name           = name;
93          this.gm             = gm;
94          this.scale          = scale;
95          this.supportedNames = supportedNames;
96          this.generateType   = generateType;
97          this.rawPVProvider  = rawPVProvider;
98          this.iauPole        = iauPole;
99          this.inertialFrame  = new InertiallyOriented(definingFrameAlignedWithICRF, inertialFrameName);
100         this.bodyFrame      = new BodyOriented(bodyOrientedFrameName);
101     }
102 
103     /** {@inheritDoc} */
104     public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate date, final Frame frame) {
105 
106         // apply the scale factor to raw position-velocity
107         final PVCoordinates rawPV    = rawPVProvider.getRawPV(date);
108         final TimeStampedPVCoordinates scaledPV = new TimeStampedPVCoordinates(date, scale, rawPV);
109 
110         // the raw PV are relative to the parent of the body centered inertially oriented frame
111         final Transform transform = getInertiallyOrientedFrame().getParent().getTransformTo(frame, date);
112 
113         // convert to requested frame
114         return transform.transformPVCoordinates(scaledPV);
115 
116     }
117 
118     /** Get the {@link FieldPVCoordinates} of the body in the selected frame.
119      * @param date current date
120      * @param frame the frame where to define the position
121      * @param <T> type of the field elements
122      * @return time-stamped position/velocity of the body (m and m/s)
123      */
124     public <T extends CalculusFieldElement<T>> TimeStampedFieldPVCoordinates<T> getPVCoordinates(final FieldAbsoluteDate<T> date,
125                                                                                                  final Frame frame) {
126 
127         // apply the scale factor to raw position-velocity
128         final FieldPVCoordinates<T>            rawPV    = rawPVProvider.getRawPV(date);
129         final TimeStampedFieldPVCoordinates<T> scaledPV = new TimeStampedFieldPVCoordinates<>(date, scale, rawPV);
130 
131         // the raw PV are relative to the parent of the body centered inertially oriented frame
132         final FieldTransform<T> transform = getInertiallyOrientedFrame().getParent().getTransformTo(frame, date);
133 
134         // convert to requested frame
135         return transform.transformPVCoordinates(scaledPV);
136 
137     }
138 
139     /** {@inheritDoc} */
140     @Override
141     public Vector3D getPosition(final AbsoluteDate date, final Frame frame) {
142 
143         // apply the scale factor to raw position
144         final Vector3D rawPosition    = rawPVProvider.getRawPosition(date);
145         final Vector3D scaledPosition = rawPosition.scalarMultiply(scale);
146 
147         // the raw position is relative to the parent of the body centered inertially oriented frame
148         final StaticTransform transform = getInertiallyOrientedFrame().getParent().getStaticTransformTo(frame, date);
149 
150         // convert to requested frame
151         return transform.transformPosition(scaledPosition);
152     }
153 
154     /** {@inheritDoc} */
155     @Override
156     public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPosition(final FieldAbsoluteDate<T> date, final Frame frame) {
157 
158         // apply the scale factor to raw position
159         final FieldVector3D<T> rawPosition     = rawPVProvider.getRawPosition(date);
160         final FieldVector3D<T> scaledPosition  = rawPosition.scalarMultiply(scale);
161 
162         // the raw position is relative to the parent of the body centered inertially oriented frame
163         final FieldStaticTransform<T> transform = getInertiallyOrientedFrame().getParent().getStaticTransformTo(frame, date);
164 
165         // convert to requested frame
166         return transform.transformPosition(scaledPosition);
167     }
168 
169     /** {@inheritDoc} */
170     public String getName() {
171         return name;
172     }
173 
174     /** {@inheritDoc} */
175     public double getGM() {
176         return gm;
177     }
178 
179     /** {@inheritDoc} */
180     public Frame getInertiallyOrientedFrame() {
181         return inertialFrame;
182     }
183 
184     /** {@inheritDoc} */
185     public Frame getBodyOrientedFrame() {
186         return bodyFrame;
187     }
188 
189     /** Inertially oriented body centered frame. */
190     private class InertiallyOriented extends Frame {
191 
192         /** Suffix for inertial frame name. */
193         private static final String INERTIAL_FRAME_SUFFIX = "/inertial";
194 
195         /** Simple constructor.
196          * @param definingFrame frame in which celestial body coordinates are defined
197          * @param frameName name to use (if null a default name will be built)
198          */
199         InertiallyOriented(final Frame definingFrame, final String frameName) {
200             super(definingFrame, new TransformProvider() {
201 
202                 /** {@inheritDoc} */
203                 public Transform getTransform(final AbsoluteDate date) {
204 
205                     // compute translation from parent frame to self
206                     final PVCoordinates pv = getPVCoordinates(date, definingFrame);
207                     final Transform translation = new Transform(date, pv.negate());
208 
209                     // compute rotation from ICRF frame to self,
210                     // as per the "Report of the IAU/IAG Working Group on Cartographic
211                     // Coordinates and Rotational Elements of the Planets and Satellites"
212                     // These definitions are common for all recent versions of this report
213                     // published every three years, the precise values of pole direction
214                     // and W angle coefficients may vary from publication year as models are
215                     // adjusted. These coefficients are not in this class, they are in the
216                     // specialized classes that do implement the getPole and getPrimeMeridianAngle
217                     // methods
218                     final Vector3D pole  = iauPole.getPole(date);
219                     final Vector3D qNode = iauPole.getNode(date);
220                     final Transform rotation =
221                                     new Transform(date, new Rotation(pole, qNode, Vector3D.PLUS_K, Vector3D.PLUS_I));
222 
223                     // update transform from parent to self
224                     return new Transform(date, translation, rotation);
225 
226                 }
227 
228                 @Override
229                 public StaticTransform getStaticTransform(final AbsoluteDate date) {
230                     // compute translation from parent frame to self
231                     final Vector3D position = getPosition(date, definingFrame);
232 
233                     // compute rotation from ICRF frame to self,
234                     // as per the "Report of the IAU/IAG Working Group on Cartographic
235                     // Coordinates and Rotational Elements of the Planets and Satellites"
236                     // These definitions are common for all recent versions of this report
237                     // published every three years, the precise values of pole direction
238                     // and W angle coefficients may vary from publication year as models are
239                     // adjusted. These coefficients are not in this class, they are in the
240                     // specialized classes that do implement the getPole and getPrimeMeridianAngle
241                     // methods
242                     final Vector3D pole  = iauPole.getPole(date);
243                     final Vector3D qNode = iauPole.getNode(date);
244                     final Rotation rotation =
245                                     new Rotation(pole, qNode, Vector3D.PLUS_K, Vector3D.PLUS_I);
246 
247                     // update transform from parent to self
248                     return StaticTransform.of(date, position.negate(), rotation);
249                 }
250 
251                 /** {@inheritDoc} */
252                 public <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
253 
254                     // compute translation from parent frame to self
255                     final FieldPVCoordinates<T> pv = getPVCoordinates(date, definingFrame);
256                     final FieldTransform<T> translation = new FieldTransform<>(date, pv.negate());
257 
258                     // compute rotation from ICRF frame to self,
259                     // as per the "Report of the IAU/IAG Working Group on Cartographic
260                     // Coordinates and Rotational Elements of the Planets and Satellites"
261                     // These definitions are common for all recent versions of this report
262                     // published every three years, the precise values of pole direction
263                     // and W angle coefficients may vary from publication year as models are
264                     // adjusted. These coefficients are not in this class, they are in the
265                     // specialized classes that do implement the getPole and getPrimeMeridianAngle
266                     // methods
267                     final FieldVector3D<T> pole  = iauPole.getPole(date);
268                     FieldVector3D<T> qNode = FieldVector3D.crossProduct(Vector3D.PLUS_K, pole);
269                     if (qNode.getNormSq().getReal() < Precision.SAFE_MIN) {
270                         qNode = FieldVector3D.getPlusI(date.getField());
271                     }
272                     final FieldTransform<T> rotation =
273                                     new FieldTransform<>(date,
274                                                     new FieldRotation<>(pole,
275                                                                     qNode,
276                                                                     FieldVector3D.getPlusK(date.getField()),
277                                                                     FieldVector3D.getPlusI(date.getField())));
278 
279                     // update transform from parent to self
280                     return new FieldTransform<>(date, translation, rotation);
281 
282                 }
283 
284                 @Override
285                 public <T extends CalculusFieldElement<T>> FieldStaticTransform<T> getStaticTransform(final FieldAbsoluteDate<T> date) {
286                     // field
287                     final Field<T> field = date.getField();
288                     // compute translation from parent frame to self
289                     final FieldVector3D<T> position = getPosition(date, definingFrame);
290 
291                     // compute rotation from ICRF frame to self,
292                     // as per the "Report of the IAU/IAG Working Group on Cartographic
293                     // Coordinates and Rotational Elements of the Planets and Satellites"
294                     // These definitions are common for all recent versions of this report
295                     // published every three years, the precise values of pole direction
296                     // and W angle coefficients may vary from publication year as models are
297                     // adjusted. These coefficients are not in this class, they are in the
298                     // specialized classes that do implement the getPole and getPrimeMeridianAngle
299                     // methods
300                     final FieldVector3D<T> pole  = iauPole.getPole(date);
301                     final FieldVector3D<T> qNode = iauPole.getNode(date);
302                     final FieldRotation<T> rotation =
303                                     new FieldRotation<>(pole, qNode, FieldVector3D.getPlusK(field), FieldVector3D.getPlusI(field));
304 
305                     // update transform from parent to self
306                     return FieldStaticTransform.of(date, position.negate(), rotation);
307                 }
308 
309             }, frameName == null ? name + INERTIAL_FRAME_SUFFIX : frameName, true);
310         }
311 
312     }
313 
314     /** Body oriented body centered frame. */
315     private class BodyOriented extends Frame {
316 
317         /**
318          * Suffix for body frame name.
319          */
320         private static final String BODY_FRAME_SUFFIX = "/rotating";
321 
322         /**
323          * Simple constructor.
324          *
325          * @param frameName name to use (if null a default name will be built)
326          */
327         BodyOriented(final String frameName) {
328             super(inertialFrame, new TransformProvider() {
329 
330                 /** {@inheritDoc} */
331                 public Transform getTransform(final AbsoluteDate date) {
332                     final double dt = 10.0;
333                     final double w0 = iauPole.getPrimeMeridianAngle(date);
334                     final double w1 = iauPole.getPrimeMeridianAngle(date.shiftedBy(dt));
335                     return new Transform(date,
336                             new Rotation(Vector3D.PLUS_K, w0, RotationConvention.FRAME_TRANSFORM),
337                             new Vector3D((w1 - w0) / dt, Vector3D.PLUS_K));
338                 }
339 
340                 /** {@inheritDoc} */
341                 public <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
342                     final double dt = 10.0;
343                     final T w0 = iauPole.getPrimeMeridianAngle(date);
344                     final T w1 = iauPole.getPrimeMeridianAngle(date.shiftedBy(dt));
345                     return new FieldTransform<>(date,
346                             new FieldRotation<>(FieldVector3D.getPlusK(date.getField()), w0,
347                                     RotationConvention.FRAME_TRANSFORM),
348                             new FieldVector3D<>(w1.subtract(w0).divide(dt), Vector3D.PLUS_K));
349                 }
350 
351             }, frameName == null ? name + BODY_FRAME_SUFFIX : frameName, false);
352         }
353     }
354 }