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.bodies;
18  
19  import java.io.Serializable;
20  
21  import org.hipparchus.CalculusFieldElement;
22  import org.hipparchus.Field;
23  import org.hipparchus.geometry.euclidean.threed.FieldRotation;
24  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
25  import org.hipparchus.geometry.euclidean.threed.Rotation;
26  import org.hipparchus.geometry.euclidean.threed.RotationConvention;
27  import org.hipparchus.geometry.euclidean.threed.Vector3D;
28  import org.hipparchus.util.Precision;
29  import org.orekit.annotation.DefaultDataContext;
30  import org.orekit.bodies.JPLEphemeridesLoader.EphemerisType;
31  import org.orekit.data.DataContext;
32  import org.orekit.errors.OrekitException;
33  import org.orekit.errors.OrekitInternalError;
34  import org.orekit.frames.FieldStaticTransform;
35  import org.orekit.frames.FieldTransform;
36  import org.orekit.frames.Frame;
37  import org.orekit.frames.StaticTransform;
38  import org.orekit.frames.Transform;
39  import org.orekit.frames.TransformProvider;
40  import org.orekit.time.AbsoluteDate;
41  import org.orekit.time.FieldAbsoluteDate;
42  import org.orekit.utils.FieldPVCoordinates;
43  import org.orekit.utils.PVCoordinates;
44  import org.orekit.utils.TimeStampedFieldPVCoordinates;
45  import org.orekit.utils.TimeStampedPVCoordinates;
46  
47  /** Implementation of the {@link CelestialBody} interface using JPL or INPOP ephemerides.
48   * @author Luc Maisonobe
49   */
50  class JPLCelestialBody implements CelestialBody {
51  
52      /** Serializable UID. */
53      private static final long serialVersionUID = 3809787672779740923L;
54  
55      /** Name of the body. */
56      private final String name;
57  
58      /** Regular expression for supported files names. */
59      private final String supportedNames;
60  
61      /** Ephemeris type to generate. */
62      private final JPLEphemeridesLoader.EphemerisType generateType;
63  
64      /** Raw position-velocity provider. */
65      private final transient JPLEphemeridesLoader.RawPVProvider rawPVProvider;
66  
67      /** Attraction coefficient of the body (m³/s²). */
68      private final double gm;
69  
70      /** Scaling factor for position-velocity. */
71      private final double scale;
72  
73      /** IAU pole. */
74      private final IAUPole iauPole;
75  
76      /** Inertially oriented, body-centered frame. */
77      private final Frame inertialFrame;
78  
79      /** Body oriented, body-centered frame. */
80      private final Frame bodyFrame;
81  
82      /** Build an instance and the underlying frame.
83       * @param name name of the body
84       * @param supportedNames regular expression for supported files names
85       * @param generateType ephemeris type to generate
86       * @param rawPVProvider raw position-velocity provider
87       * @param gm attraction coefficient (in m³/s²)
88       * @param scale scaling factor for position-velocity
89       * @param iauPole IAU pole implementation
90       * @param definingFrameAlignedWithICRF frame in which celestial body coordinates are defined,
91       * this frame <strong>must</strong> be aligned with ICRF
92       * @param inertialFrameName name to use for inertial frame (if null a default name will be built)
93       * @param bodyOrientedFrameName name to use for body-oriented frame (if null a default name will be built)
94       */
95      JPLCelestialBody(final String name, final String supportedNames,
96                       final JPLEphemeridesLoader.EphemerisType generateType,
97                       final JPLEphemeridesLoader.RawPVProvider rawPVProvider,
98                       final double gm, final double scale,
99                       final IAUPole iauPole, final Frame definingFrameAlignedWithICRF,
100                      final String inertialFrameName, final String bodyOrientedFrameName) {
101         this.name           = name;
102         this.gm             = gm;
103         this.scale          = scale;
104         this.supportedNames = supportedNames;
105         this.generateType   = generateType;
106         this.rawPVProvider  = rawPVProvider;
107         this.iauPole        = iauPole;
108         this.inertialFrame  = new InertiallyOriented(definingFrameAlignedWithICRF, inertialFrameName);
109         this.bodyFrame      = new BodyOriented(bodyOrientedFrameName);
110     }
111 
112     /** {@inheritDoc} */
113     public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate date, final Frame frame) {
114 
115         // apply the scale factor to raw position-velocity
116         final PVCoordinates rawPV    = rawPVProvider.getRawPV(date);
117         final TimeStampedPVCoordinates scaledPV = new TimeStampedPVCoordinates(date, scale, rawPV);
118 
119         // the raw PV are relative to the parent of the body centered inertially oriented frame
120         final Transform transform = getInertiallyOrientedFrame().getParent().getTransformTo(frame, date);
121 
122         // convert to requested frame
123         return transform.transformPVCoordinates(scaledPV);
124 
125     }
126 
127     /** Get the {@link FieldPVCoordinates} of the body in the selected frame.
128      * @param date current date
129      * @param frame the frame where to define the position
130      * @param <T> type of the field elements
131      * @return time-stamped position/velocity of the body (m and m/s)
132      */
133     public <T extends CalculusFieldElement<T>> TimeStampedFieldPVCoordinates<T> getPVCoordinates(final FieldAbsoluteDate<T> date,
134                                                                                                  final Frame frame) {
135 
136         // apply the scale factor to raw position-velocity
137         final FieldPVCoordinates<T>            rawPV    = rawPVProvider.getRawPV(date);
138         final TimeStampedFieldPVCoordinates<T> scaledPV = new TimeStampedFieldPVCoordinates<>(date, scale, rawPV);
139 
140         // the raw PV are relative to the parent of the body centered inertially oriented frame
141         final FieldTransform<T> transform = getInertiallyOrientedFrame().getParent().getTransformTo(frame, date);
142 
143         // convert to requested frame
144         return transform.transformPVCoordinates(scaledPV);
145 
146     }
147 
148     /** {@inheritDoc} */
149     @Override
150     public Vector3D getPosition(final AbsoluteDate date, final Frame frame) {
151 
152         // apply the scale factor to raw position
153         final Vector3D rawPosition    = rawPVProvider.getRawPosition(date);
154         final Vector3D scaledPosition = rawPosition.scalarMultiply(scale);
155 
156         // the raw position is relative to the parent of the body centered inertially oriented frame
157         final StaticTransform transform = getInertiallyOrientedFrame().getParent().getStaticTransformTo(frame, date);
158 
159         // convert to requested frame
160         return transform.transformPosition(scaledPosition);
161     }
162 
163     /** {@inheritDoc} */
164     @Override
165     public <T extends CalculusFieldElement<T>> FieldVector3D<T> getPosition(final FieldAbsoluteDate<T> date, final Frame frame) {
166 
167         // apply the scale factor to raw position
168         final FieldVector3D<T> rawPosition     = rawPVProvider.getRawPosition(date);
169         final FieldVector3D<T> scaledPosition  = rawPosition.scalarMultiply(scale);
170 
171         // the raw position is relative to the parent of the body centered inertially oriented frame
172         final FieldStaticTransform<T> transform = getInertiallyOrientedFrame().getParent().getStaticTransformTo(frame, date);
173 
174         // convert to requested frame
175         return transform.transformPosition(scaledPosition);
176     }
177 
178 
179     /** Replace the instance with a data transfer object for serialization.
180      * <p>
181      * This intermediate class serializes the files supported names, the ephemeris type
182      * and the body name.
183      * </p>
184      * @return data transfer object that will be serialized
185      */
186     @DefaultDataContext
187     private Object writeReplace() {
188         return new DTOCelestialBody(supportedNames, generateType, name);
189     }
190 
191     /** {@inheritDoc} */
192     public String getName() {
193         return name;
194     }
195 
196     /** {@inheritDoc} */
197     public double getGM() {
198         return gm;
199     }
200 
201     /** {@inheritDoc} */
202     public Frame getInertiallyOrientedFrame() {
203         return inertialFrame;
204     }
205 
206     /** {@inheritDoc} */
207     public Frame getBodyOrientedFrame() {
208         return bodyFrame;
209     }
210 
211     /** Inertially oriented body centered frame. */
212     private class InertiallyOriented extends Frame {
213 
214         /** Serializable UID. */
215         private static final long serialVersionUID = -8849993808761896559L;
216 
217         /** Suffix for inertial frame name. */
218         private static final String INERTIAL_FRAME_SUFFIX = "/inertial";
219 
220         /** Simple constructor.
221          * @param definingFrame frame in which celestial body coordinates are defined
222          * @param frameName name to use (if null a default name will be built)
223          */
224         InertiallyOriented(final Frame definingFrame, final String frameName) {
225             super(definingFrame, new TransformProvider() {
226 
227                 /** Serializable UID. */
228                 private static final long serialVersionUID = -8610328386110652400L;
229 
230                 /** {@inheritDoc} */
231                 public Transform getTransform(final AbsoluteDate date) {
232 
233                     // compute translation from parent frame to self
234                     final PVCoordinates pv = getPVCoordinates(date, definingFrame);
235                     final Transform translation = new Transform(date, pv.negate());
236 
237                     // compute rotation from ICRF frame to self,
238                     // as per the "Report of the IAU/IAG Working Group on Cartographic
239                     // Coordinates and Rotational Elements of the Planets and Satellites"
240                     // These definitions are common for all recent versions of this report
241                     // published every three years, the precise values of pole direction
242                     // and W angle coefficients may vary from publication year as models are
243                     // adjusted. These coefficients are not in this class, they are in the
244                     // specialized classes that do implement the getPole and getPrimeMeridianAngle
245                     // methods
246                     final Vector3D pole  = iauPole.getPole(date);
247                     final Vector3D qNode = iauPole.getNode(date);
248                     final Transform rotation =
249                                     new Transform(date, new Rotation(pole, qNode, Vector3D.PLUS_K, Vector3D.PLUS_I));
250 
251                     // update transform from parent to self
252                     return new Transform(date, translation, rotation);
253 
254                 }
255 
256                 @Override
257                 public StaticTransform getStaticTransform(final AbsoluteDate date) {
258                     // compute translation from parent frame to self
259                     final Vector3D position = getPosition(date, definingFrame);
260 
261                     // compute rotation from ICRF frame to self,
262                     // as per the "Report of the IAU/IAG Working Group on Cartographic
263                     // Coordinates and Rotational Elements of the Planets and Satellites"
264                     // These definitions are common for all recent versions of this report
265                     // published every three years, the precise values of pole direction
266                     // and W angle coefficients may vary from publication year as models are
267                     // adjusted. These coefficients are not in this class, they are in the
268                     // specialized classes that do implement the getPole and getPrimeMeridianAngle
269                     // methods
270                     final Vector3D pole  = iauPole.getPole(date);
271                     final Vector3D qNode = iauPole.getNode(date);
272                     final Rotation rotation =
273                                     new Rotation(pole, qNode, Vector3D.PLUS_K, Vector3D.PLUS_I);
274 
275                     // update transform from parent to self
276                     return StaticTransform.of(date, position.negate(), rotation);
277                 }
278 
279                 /** {@inheritDoc} */
280                 public <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
281 
282                     // compute translation from parent frame to self
283                     final FieldPVCoordinates<T> pv = getPVCoordinates(date, definingFrame);
284                     final FieldTransform<T> translation = new FieldTransform<>(date, pv.negate());
285 
286                     // compute rotation from ICRF frame to self,
287                     // as per the "Report of the IAU/IAG Working Group on Cartographic
288                     // Coordinates and Rotational Elements of the Planets and Satellites"
289                     // These definitions are common for all recent versions of this report
290                     // published every three years, the precise values of pole direction
291                     // and W angle coefficients may vary from publication year as models are
292                     // adjusted. These coefficients are not in this class, they are in the
293                     // specialized classes that do implement the getPole and getPrimeMeridianAngle
294                     // methods
295                     final FieldVector3D<T> pole  = iauPole.getPole(date);
296                     FieldVector3D<T> qNode = FieldVector3D.crossProduct(Vector3D.PLUS_K, pole);
297                     if (qNode.getNormSq().getReal() < Precision.SAFE_MIN) {
298                         qNode = FieldVector3D.getPlusI(date.getField());
299                     }
300                     final FieldTransform<T> rotation =
301                                     new FieldTransform<>(date,
302                                                     new FieldRotation<>(pole,
303                                                                     qNode,
304                                                                     FieldVector3D.getPlusK(date.getField()),
305                                                                     FieldVector3D.getPlusI(date.getField())));
306 
307                     // update transform from parent to self
308                     return new FieldTransform<>(date, translation, rotation);
309 
310                 }
311 
312                 @Override
313                 public <T extends CalculusFieldElement<T>> FieldStaticTransform<T> getStaticTransform(final FieldAbsoluteDate<T> date) {
314                     // field
315                     final Field<T> field = date.getField();
316                     // compute translation from parent frame to self
317                     final FieldVector3D<T> position = getPosition(date, definingFrame);
318 
319                     // compute rotation from ICRF frame to self,
320                     // as per the "Report of the IAU/IAG Working Group on Cartographic
321                     // Coordinates and Rotational Elements of the Planets and Satellites"
322                     // These definitions are common for all recent versions of this report
323                     // published every three years, the precise values of pole direction
324                     // and W angle coefficients may vary from publication year as models are
325                     // adjusted. These coefficients are not in this class, they are in the
326                     // specialized classes that do implement the getPole and getPrimeMeridianAngle
327                     // methods
328                     final FieldVector3D<T> pole  = iauPole.getPole(date);
329                     final FieldVector3D<T> qNode = iauPole.getNode(date);
330                     final FieldRotation<T> rotation =
331                                     new FieldRotation<>(pole, qNode, FieldVector3D.getPlusK(field), FieldVector3D.getPlusI(field));
332 
333                     // update transform from parent to self
334                     return FieldStaticTransform.of(date, position.negate(), rotation);
335                 }
336 
337             }, frameName == null ? name + INERTIAL_FRAME_SUFFIX : frameName, true);
338         }
339 
340         /** Replace the instance with a data transfer object for serialization.
341          * <p>
342          * This intermediate class serializes the files supported names, the ephemeris type
343          * and the body name.
344          * </p>
345          * @return data transfer object that will be serialized
346          */
347         @DefaultDataContext
348         private Object writeReplace() {
349             return new DTOInertialFrame(supportedNames, generateType, name);
350         }
351 
352     }
353 
354     /** Body oriented body centered frame. */
355     private class BodyOriented extends Frame {
356 
357         /** Serializable UID. */
358         private static final long serialVersionUID = 20170109L;
359 
360         /** Suffix for body frame name. */
361         private static final String BODY_FRAME_SUFFIX = "/rotating";
362 
363         /** Simple constructor.
364          * @param frameName name to use (if null a default name will be built)
365          */
366         BodyOriented(final String frameName) {
367             super(inertialFrame, new TransformProvider() {
368 
369                 /** Serializable UID. */
370                 private static final long serialVersionUID = 20170109L;
371 
372                 /** {@inheritDoc} */
373                 public Transform getTransform(final AbsoluteDate date) {
374                     final double dt = 10.0;
375                     final double w0 = iauPole.getPrimeMeridianAngle(date);
376                     final double w1 = iauPole.getPrimeMeridianAngle(date.shiftedBy(dt));
377                     return new Transform(date,
378                                          new Rotation(Vector3D.PLUS_K, w0, RotationConvention.FRAME_TRANSFORM),
379                                          new Vector3D((w1 - w0) / dt, Vector3D.PLUS_K));
380                 }
381 
382                 /** {@inheritDoc} */
383                 public <T extends CalculusFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
384                     final double dt = 10.0;
385                     final T w0 = iauPole.getPrimeMeridianAngle(date);
386                     final T w1 = iauPole.getPrimeMeridianAngle(date.shiftedBy(dt));
387                     return new FieldTransform<>(date,
388                                     new FieldRotation<>(FieldVector3D.getPlusK(date.getField()), w0,
389                                                     RotationConvention.FRAME_TRANSFORM),
390                                     new FieldVector3D<>(w1.subtract(w0).divide(dt), Vector3D.PLUS_K));
391                 }
392 
393             }, frameName == null ? name + BODY_FRAME_SUFFIX : frameName, false);
394         }
395 
396         /** Replace the instance with a data transfer object for serialization.
397          * <p>
398          * This intermediate class serializes the files supported names, the ephemeris type
399          * and the body name.
400          * </p>
401          * @return data transfer object that will be serialized
402          */
403         @DefaultDataContext
404         private Object writeReplace() {
405             return new DTOBodyFrame(supportedNames, generateType, name);
406         }
407 
408     }
409 
410     /** Internal class used only for serialization. */
411     @DefaultDataContext
412     private abstract static class DataTransferObject implements Serializable {
413 
414         /** Serializable UID. */
415         private static final long serialVersionUID = 674742836536072422L;
416 
417         /** Regular expression for supported files names. */
418         private final String supportedNames;
419 
420         /** Ephemeris type to generate. */
421         private final EphemerisType generateType;
422 
423         /** Name of the body. */
424         private final String name;
425 
426         /** Simple constructor.
427          * @param supportedNames regular expression for supported files names
428          * @param generateType ephemeris type to generate
429          * @param name name of the body
430          */
431         DataTransferObject(final String supportedNames, final EphemerisType generateType, final String name) {
432             this.supportedNames = supportedNames;
433             this.generateType   = generateType;
434             this.name           = name;
435         }
436 
437         /** Get the body associated with the serialized data.
438          * @return body associated with the serialized data
439          */
440         protected JPLCelestialBody getBody() {
441 
442             try {
443                 // first try to use the factory, in order to avoid building a new instance
444                 // each time we deserialize and have the object properly cached
445                 final CelestialBody factoryProvided =
446                                 DataContext.getDefault().getCelestialBodies().getBody(name);
447                 if (factoryProvided instanceof JPLCelestialBody) {
448                     final JPLCelestialBody jplBody = (JPLCelestialBody) factoryProvided;
449                     if (supportedNames.equals(jplBody.supportedNames) && generateType == jplBody.generateType) {
450                         // the factory created exactly the object we needed, just return it
451                         return jplBody;
452                     }
453                 }
454 
455                 // the factory does not return the object we want
456                 // we create a new one from scratch and don't cache it
457                 return (JPLCelestialBody) new JPLEphemeridesLoader(supportedNames, generateType).loadCelestialBody(name);
458 
459             } catch (OrekitException oe) {
460                 throw new OrekitInternalError(oe);
461             }
462 
463         }
464 
465     }
466 
467     /** Specialization of the data transfer object for complete celestial body serialization. */
468     @DefaultDataContext
469     private static class DTOCelestialBody extends DataTransferObject {
470 
471         /** Serializable UID. */
472         private static final long serialVersionUID = -8287341529741045958L;
473 
474         /** Simple constructor.
475          * @param supportedNames regular expression for supported files names
476          * @param generateType ephemeris type to generate
477          * @param name name of the body
478          */
479         DTOCelestialBody(final String supportedNames, final EphemerisType generateType, final String name) {
480             super(supportedNames, generateType, name);
481         }
482 
483         /** Replace the deserialized data transfer object with a {@link JPLCelestialBody}.
484          * @return replacement {@link JPLCelestialBody}
485          */
486         private Object readResolve() {
487             return getBody();
488         }
489 
490     }
491 
492     /** Specialization of the data transfer object for inertially oriented frame serialization. */
493     @DefaultDataContext
494     private static class DTOInertialFrame extends DataTransferObject {
495 
496         /** Serializable UID. */
497         private static final long serialVersionUID = 7915071664444154948L;
498 
499         /** Simple constructor.
500          * @param supportedNames regular expression for supported files names
501          * @param generateType ephemeris type to generate
502          * @param name name of the body
503          */
504         DTOInertialFrame(final String supportedNames, final EphemerisType generateType, final String name) {
505             super(supportedNames, generateType, name);
506         }
507 
508         /** Replace the deserialized data transfer object with a {@link Frame}.
509          * @return replacement {@link Frame}
510          */
511         private Object readResolve() {
512             return getBody().inertialFrame;
513         }
514 
515     }
516 
517     /** Specialization of the data transfer object for body oriented frame serialization. */
518     @DefaultDataContext
519     private static class DTOBodyFrame extends DataTransferObject {
520 
521         /** Serializable UID. */
522         private static final long serialVersionUID = -3194195019557081000L;
523 
524         /** Simple constructor.
525          * @param supportedNames regular expression for supported files names
526          * @param generateType ephemeris type to generate
527          * @param name name of the body
528          */
529         DTOBodyFrame(final String supportedNames, final EphemerisType generateType, final String name) {
530             super(supportedNames, generateType, name);
531         }
532 
533         /** Replace the deserialized data transfer object with a {@link Frame}.
534          * @return replacement {@link Frame}
535          */
536         private Object readResolve() {
537             return getBody().bodyFrame;
538         }
539 
540     }
541 
542 }