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