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