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  import org.hipparchus.geometry.euclidean.threed.Rotation;
20  import org.hipparchus.geometry.euclidean.threed.Vector3D;
21  import org.junit.jupiter.api.Assertions;
22  import org.junit.jupiter.api.Test;
23  import org.orekit.Utils;
24  import org.orekit.errors.OrekitException;
25  import org.orekit.frames.Frame;
26  import org.orekit.frames.FramesFactory;
27  import org.orekit.time.AbsoluteDate;
28  import org.orekit.time.TimeScalesFactory;
29  import org.orekit.utils.Constants;
30  import org.orekit.utils.IERSConventions;
31  import org.orekit.utils.PVCoordinates;
32  import org.orekit.utils.TimeStampedPVCoordinates;
33  
34  import java.util.ArrayList;
35  import java.util.Arrays;
36  import java.util.List;
37  import java.util.concurrent.ExecutorService;
38  import java.util.concurrent.Executors;
39  import java.util.concurrent.Future;
40  import java.util.concurrent.TimeUnit;
41  import java.util.concurrent.atomic.AtomicReference;
42  
43  public class CelestialBodyFactoryTest {
44  
45      @Test
46      public void getSun() {
47          Utils.setDataRoot("regular-data");
48  
49          CelestialBody sun = CelestialBodyFactory.getSun();
50          Assertions.assertNotNull(sun);
51      }
52  
53      @Test
54      public void clearCache() {
55          Utils.setDataRoot("regular-data");
56  
57          CelestialBody sun = CelestialBodyFactory.getSun();
58          Assertions.assertNotNull(sun);
59          CelestialBodyFactory.clearCelestialBodyCache();
60          CelestialBody sun2 = CelestialBodyFactory.getSun();
61          Assertions.assertNotNull(sun2);
62          Assertions.assertNotSame(sun, sun2);
63      }
64  
65      @Test
66      public void clearLoaders() {
67          Utils.setDataRoot("regular-data");
68  
69          CelestialBody sun = CelestialBodyFactory.getSun();
70          Assertions.assertNotNull(sun);
71          CelestialBodyFactory.clearCelestialBodyLoaders();
72          CelestialBody sun2 = CelestialBodyFactory.getSun();
73          Assertions.assertNotNull(sun2);
74          Assertions.assertNotSame(sun, sun2);
75          CelestialBodyFactory.clearCelestialBodyLoaders(CelestialBodyFactory.SUN);
76          CelestialBodyFactory.clearCelestialBodyCache(CelestialBodyFactory.SUN);
77          CelestialBodyFactory.addDefaultCelestialBodyLoader(JPLEphemeridesLoader.DEFAULT_DE_SUPPORTED_NAMES);
78          CelestialBody sun3 = CelestialBodyFactory.getSun();
79          Assertions.assertNotNull(sun3);
80          Assertions.assertNotSame(sun,  sun3);
81          Assertions.assertNotSame(sun2, sun3);
82      }
83  
84      @Test
85      void testHorizon() {
86  
87          // The following data are an excerpt from a telnet session with JPL Horizon system
88          // note that in Horizon we selected Jupiter barycenter rather than Jupiter body center
89          // this seems to match better the content of the DE-431 ephemeris
90  
91          //  *******************************************************************************
92          //  Ephemeris / PORT_LOGIN Mon Oct 26 04:53:43 2015 Pasadena, USA    / Horizons
93          //  *******************************************************************************
94          //  Target body name: Jupiter Barycenter (5)          {source: DE-0431LE-0431}
95          //  Center body name: Solar System Barycenter (0)     {source: DE-0431LE-0431}
96          //  Center-site name: BODY CENTER
97          //  *******************************************************************************
98          //  Start time      : A.D. 2000-Jan-01 00:00:00.0000 CT
99          //  Stop  time      : A.D. 2003-Dec-31 23:59:00.0000 CT
100         //  Step-size       : 1440 minutes
101         //  *******************************************************************************
102         //  Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)}
103         //  Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)}
104         //  Center radii    : (undefined)
105         //  Output units    : KM-S
106         //  Output format   : 02
107         //  Reference frame : ICRF/J2000.0
108         //  Output type     : GEOMETRIC cartesian states
109         //  Coordinate systm: Earth Mean Equator and Equinox of Reference Epoch
110         //  *******************************************************************************
111         //  JDCT
112         //     X     Y     Z
113         //     VX    VY    VZ
114         //  *******************************************************************************
115         //  $$SOE
116         //  2451544.500000000 = A.D. 2000-Jan-01 00:00:00.0000 (TDB)
117         //   X = 5.978411018921824E+08 Y = 4.085508359611598E+08 Z = 1.605595308103096E+08
118         //   VX=-7.892151874487445E+00 VY= 1.017751699703826E+01 VZ= 4.554715748011852E+00
119         //  2451545.500000000 = A.D. 2000-Jan-02 00:00:00.0000 (TDB)
120         //   X = 5.971584965869523E+08 Y = 4.094296790808872E+08 Z = 1.609528639632485E+08
121         //   VX=-7.908893450088906E+00 VY= 1.016606978596496E+01 VZ= 4.550216570971850E+00
122         //  2451546.500000000 = A.D. 2000-Jan-03 00:00:00.0000 (TDB)
123         //   X = 5.964744456934582E+08 Y = 4.103075321378759E+08 Z = 1.613458079269412E+08
124         //   VX=-7.925614558638352E+00 VY= 1.015459888397081E+01 VZ= 4.545706740033853E+00
125         //  2451547.500000000 = A.D. 2000-Jan-04 00:00:00.0000 (TDB)
126         //   X = 5.957889509819047E+08 Y = 4.111843930867567E+08 Z = 1.617383617815004E+08
127         //   VX=-7.942315157290042E+00 VY= 1.014310432640078E+01 VZ= 4.541186269306714E+00
128         //  2451548.500000000 = A.D. 2000-Jan-05 00:00:00.0000 (TDB)
129         //   X = 5.951020142261952E+08 Y = 4.120602598852173E+08 Z = 1.621305246082588E+08
130         //   VX=-7.958995203281466E+00 VY= 1.013158614867071E+01 VZ= 4.536655172931710E+00
131         //  2451549.500000000 = A.D. 2000-Jan-06 00:00:00.0000 (TDB)
132         //   X = 5.944136372039236E+08 Y = 4.129351304940084E+08 Z = 1.625222954897725E+08
133         //   VX=-7.975654653933506E+00 VY= 1.012004438626713E+01 VZ= 4.532113465082445E+00
134         final TimeStampedPVCoordinates[] refPV = new TimeStampedPVCoordinates[] {
135             createPV(2000, 1, 1,
136                      5.978411018921824E+08, 4.085508359611598E+08, 1.605595308103096E+08,
137                      -7.892151874487445E+00, 1.017751699703826E+01, 4.554715748011852E+00),
138             createPV(2000, 1, 2,
139                      5.971584965869523E+08, 4.094296790808872E+08, 1.609528639632485E+08,
140              -7.908893450088906E+00, 1.016606978596496E+01, 4.550216570971850E+00),
141             createPV(2000, 1, 3,
142                      5.964744456934582E+08, 4.103075321378759E+08, 1.613458079269412E+08,
143              -7.925614558638352E+00, 1.015459888397081E+01, 4.545706740033853E+00),
144             createPV(2000, 1, 4,
145                      5.957889509819047E+08, 4.111843930867567E+08, 1.617383617815004E+08,
146                      -7.942315157290042E+00, 1.014310432640078E+01, 4.541186269306714E+00),
147             createPV(2000, 1, 5,
148                      5.951020142261952E+08, 4.120602598852173E+08, 1.621305246082588E+08,
149                      -7.958995203281466E+00, 1.013158614867071E+01, 4.536655172931710E+00),
150             createPV(2000, 1, 6,
151                      5.944136372039236E+08, 4.129351304940084E+08, 1.625222954897725E+08,
152                      -7.975654653933506E+00, 1.012004438626713E+01, 4.532113465082445E+00)
153         };
154 
155         Utils.setDataRoot("regular-data");
156         final CelestialBody jupiter = CelestialBodyFactory.getJupiter();
157         final Frame icrf = FramesFactory.getICRF();
158         for (final TimeStampedPVCoordinates ref : refPV) {
159             TimeStampedPVCoordinates testPV = jupiter.getPVCoordinates(ref.getDate(),
160                     icrf);
161             Assertions.assertEquals(0.0,
162                                 Vector3D.distance(ref.getPosition(), testPV.getPosition()),
163                                 4.0e-4);
164             Assertions.assertEquals(0.0,
165                                 Vector3D.distance(ref.getVelocity(), testPV.getVelocity()),
166                                 1.0e-11);
167             Vector3D testP = jupiter.getInertiallyOrientedFrame()
168                     .getStaticTransformTo(icrf, ref.getDate())
169                     .transformPosition(Vector3D.ZERO);
170             Assertions.assertEquals(
171                     0.0,
172                     Vector3D.distance(ref.getPosition(), testP),
173                     8.0e-4);
174             testP = jupiter.getBodyOrientedFrame()
175                     .getStaticTransformTo(icrf, ref.getDate())
176                     .transformPosition(Vector3D.ZERO);
177             Assertions.assertEquals(
178                     0.0,
179                     Vector3D.distance(ref.getPosition(), testP),
180                     8.0e-4);
181         }
182 
183     }
184 
185     private TimeStampedPVCoordinates createPV(int year, int month, int day,
186                                               double xKm, double yKm, double zKM,
187                                               double vxKmS, double vyKms, double vzKms) {
188         return new TimeStampedPVCoordinates(new AbsoluteDate(year, month, day, TimeScalesFactory.getTDB()),
189                                             new Vector3D(  xKm * 1000,   yKm * 1000,   zKM * 1000),
190                                             new Vector3D(vxKmS * 1000, vyKms * 1000, vzKms * 1000));
191     }
192 
193     @Test
194     public void multithreadTest() {
195         Utils.setDataRoot("regular-data");
196         checkMultiThread(10, 100);
197     }
198 
199     private void checkMultiThread(final int threads, final int runs) {
200 
201         final AtomicReference<OrekitException> caught = new AtomicReference<OrekitException>();
202         ExecutorService executorService = Executors.newFixedThreadPool(threads);
203 
204         List<Future<?>> results = new ArrayList<Future<?>>();
205         for (int i = 0; i < threads; i++) {
206             Future<?> result = executorService.submit(new Runnable() {
207                 public void run() {
208                     try {
209                         for (int run = 0; run < runs; run++) {
210                             CelestialBody mars = CelestialBodyFactory.getBody(CelestialBodyFactory.MARS);
211                             Assertions.assertNotNull(mars);
212                             CelestialBodyFactory.clearCelestialBodyLoaders();
213                         }
214                     } catch (OrekitException oe) {
215                         caught.set(oe);
216                     }
217                 }
218             });
219             results.add(result);
220         }
221 
222         try {
223             executorService.shutdown();
224             executorService.awaitTermination(5, TimeUnit.SECONDS);
225         } catch (InterruptedException ie) {
226             Assertions.fail(ie.getLocalizedMessage());
227         }
228 
229         for (Future<?> result : results) {
230             Assertions.assertTrue(result.isDone(),"Not all threads finished -> possible deadlock");
231         }
232 
233         if (caught.get() != null) {
234             throw caught.get();
235         }
236     }
237 
238     @Test
239     void testEarthMoonBarycenter() {
240         Utils.setDataRoot("regular-data/de405-ephemerides");
241         CelestialBody sun = CelestialBodyFactory.getSun();
242         CelestialBody mars = CelestialBodyFactory.getMars();
243         CelestialBody earth = CelestialBodyFactory.getEarth();
244         CelestialBody earthMoonBarycenter = CelestialBodyFactory.getEarthMoonBarycenter();
245         List<Frame> frames = Arrays.asList(FramesFactory.getEME2000(),
246                                            FramesFactory.getGCRF(),
247                                            sun.getInertiallyOrientedFrame(),
248                                            mars.getInertiallyOrientedFrame(),
249                                            earth.getInertiallyOrientedFrame());
250 
251         AbsoluteDate date = new AbsoluteDate(1969, 7, 23, TimeScalesFactory.getTT());
252         final double refDistance = bodyDistance(sun, earthMoonBarycenter, date, frames.get(0));
253         for (Frame frame : frames) {
254             Assertions.assertEquals(refDistance,
255                                 bodyDistance(sun, earthMoonBarycenter, date, frame),
256                                 1.0e-14 * refDistance,frame.toString());
257         }
258     }
259 
260     @Test
261     void testICRFAndGCRFAlignment() {
262         Utils.setDataRoot("regular-data");
263         final CelestialBody earthMoonBarycenter   = CelestialBodyFactory.getEarthMoonBarycenter();
264         final CelestialBody solarSystemBarycenter = CelestialBodyFactory.getSolarSystemBarycenter();
265         final List<Frame> frames = Arrays.asList(earthMoonBarycenter.getInertiallyOrientedFrame(),
266                                                  earthMoonBarycenter.getBodyOrientedFrame(),
267                                                  solarSystemBarycenter.getInertiallyOrientedFrame(),
268                                                  solarSystemBarycenter.getBodyOrientedFrame());
269         final Frame icrf = FramesFactory.getICRF();
270         final Frame gcrf = FramesFactory.getGCRF();
271         for (double dt = 0; dt < Constants.JULIAN_DAY; dt += 60) {
272             final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(dt);
273             for (final Frame frame : frames) {
274                 Assertions.assertEquals(0.0, frame.getTransformTo(icrf, date).getRotation().getAngle(), 1.0e-15);
275                 Assertions.assertEquals(0.0, frame.getTransformTo(gcrf, date).getRotation().getAngle(), 1.0e-15);
276                 Assertions.assertEquals(0.0, frame.getStaticTransformTo(icrf, date).getRotation().getAngle(), 1.0e-15);
277                 Assertions.assertEquals(0.0, frame.getStaticTransformTo(gcrf, date).getRotation().getAngle(), 1.0e-15);
278             }
279         }
280     }
281 
282     @Test
283     void testEarthInertialFrameAroundJ2000() {
284         Utils.setDataRoot("regular-data");
285         final Frame earthFrame = CelestialBodyFactory.getEarth().getInertiallyOrientedFrame();
286         final Frame base       = FramesFactory.getGCRF();
287         final Rotation reference = new Rotation(Vector3D.PLUS_K, Vector3D.PLUS_J,
288                                                 Vector3D.PLUS_K, Vector3D.PLUS_I);
289          for (double dt = -60; dt <= 60; dt += 1.0) {
290              final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(dt);
291              Rotation rotation = base.getTransformTo(earthFrame, date).getRotation();
292              Assertions.assertEquals(0.0, Rotation.distance(reference, rotation), 3.0e-10);
293          }
294     }
295 
296     @Test
297     void testEarthBodyOrientedFrameAroundJ2000() {
298         Utils.setDataRoot("regular-data");
299         final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
300         final Frame base       = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
301          for (double dt = -60; dt <= 60; dt += 1.0) {
302              final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(dt);
303              Rotation rotation = base.getTransformTo(earthFrame, date).getRotation();
304              Assertions.assertEquals(7.9426e-4, Rotation.distance(Rotation.IDENTITY, rotation), 1.0e-8);
305          }
306     }
307 
308     private double bodyDistance(CelestialBody body1, CelestialBody body2, AbsoluteDate date, Frame frame)
309         {
310         Vector3D body1Position = body1.getPosition(date, frame);
311         Vector3D body2Position = body2.getPosition(date, frame);
312         Vector3D bodyPositionDifference = body1Position.subtract(body2Position);
313         return bodyPositionDifference.getNorm();
314     }
315 
316     @Test
317     void testGetVelocity() {
318         // GIVEN
319         Utils.setDataRoot("regular-data");
320         final CelestialBody moon = CelestialBodyFactory.getSun();
321         final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
322         final Frame frame = FramesFactory.getGCRF();
323         // WHEN
324         final Vector3D velocity = moon.getVelocity(date, frame);
325         // THEN
326         final PVCoordinates expectedPV = moon.getPVCoordinates(date, frame);
327         Assertions.assertArrayEquals(expectedPV.getVelocity().toArray(), velocity.toArray(), 1.0e-10);
328     }
329 }