1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.bodies;
18
19 import java.io.BufferedReader;
20 import java.io.IOException;
21 import java.io.InputStream;
22 import java.io.InputStreamReader;
23 import java.io.UnsupportedEncodingException;
24 import java.nio.charset.StandardCharsets;
25 import java.util.Collections;
26 import java.util.List;
27 import java.util.stream.Stream;
28
29 import org.hipparchus.CalculusFieldElement;
30 import org.hipparchus.Field;
31 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
32 import org.hipparchus.geometry.euclidean.threed.Vector3D;
33 import org.hipparchus.ode.AbstractIntegrator;
34 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
35 import org.hipparchus.util.Binary64;
36 import org.hipparchus.util.Binary64Field;
37 import org.hipparchus.util.FastMath;
38 import org.junit.jupiter.api.Assertions;
39 import org.junit.jupiter.api.Test;
40 import org.orekit.Utils;
41 import org.orekit.errors.OrekitException;
42 import org.orekit.forces.ForceModel;
43 import org.orekit.frames.Frame;
44 import org.orekit.frames.FramesFactory;
45 import org.orekit.frames.Transform;
46 import org.orekit.orbits.CartesianOrbit;
47 import org.orekit.orbits.KeplerianOrbit;
48 import org.orekit.orbits.Orbit;
49 import org.orekit.orbits.OrbitType;
50 import org.orekit.propagation.FieldSpacecraftState;
51 import org.orekit.propagation.SpacecraftState;
52 import org.orekit.propagation.ToleranceProvider;
53 import org.orekit.propagation.analytical.KeplerianPropagator;
54 import org.orekit.propagation.events.EventDetector;
55 import org.orekit.propagation.events.FieldEventDetector;
56 import org.orekit.propagation.numerical.NumericalPropagator;
57 import org.orekit.time.AbsoluteDate;
58 import org.orekit.time.FieldAbsoluteDate;
59 import org.orekit.time.TimeScale;
60 import org.orekit.time.TimeScalesFactory;
61 import org.orekit.utils.Constants;
62 import org.orekit.utils.PVCoordinates;
63 import org.orekit.utils.PVCoordinatesProvider;
64 import org.orekit.utils.ParameterDriver;
65
66 public class SolarBodyTest {
67
68 @Test
69 void testNaif() throws UnsupportedEncodingException, IOException {
70 Utils.setDataRoot("regular-data");
71 final Frame refFrame = FramesFactory.getICRF();
72 final TimeScale tdb = TimeScalesFactory.getTDB();
73 final InputStream inEntry = getClass().getResourceAsStream("/naif/DE431-ephemeris-NAIF.txt");
74 try (BufferedReader reader = new BufferedReader(new InputStreamReader(inEntry, StandardCharsets.UTF_8))) {
75 for (String line = reader.readLine(); line != null; line = reader.readLine()) {
76 line = line.trim();
77 if (!line.isEmpty() && !line.startsWith("#")) {
78
79
80 String[] fields = line.split("\\s+");
81 final AbsoluteDate date1 = new AbsoluteDate(fields[0], tdb);
82 final AbsoluteDate date2 = new AbsoluteDate(AbsoluteDate.J2000_EPOCH,
83 Double.parseDouble(fields[1]),
84 tdb);
85 String name = fields[2];
86 final String barycenter = fields[3];
87 final Vector3D pRef = new Vector3D(Double.parseDouble(fields[4]) * 1000.0,
88 Double.parseDouble(fields[5]) * 1000.0,
89 Double.parseDouble(fields[6]) * 1000.0);
90 final Vector3D vRef = new Vector3D(Double.parseDouble(fields[7]) * 1000.0,
91 Double.parseDouble(fields[8]) * 1000.0,
92 Double.parseDouble(fields[9]) * 1000.0);
93
94
95 Assertions.assertEquals("BARYCENTER", barycenter);
96 if (name.equals("EARTH")) {
97 name = "EARTH-MOON BARYCENTER";
98 }
99 Assertions.assertEquals(0.0, date2.durationFrom(date1), 8.0e-5);
100 final PVCoordinates pv = CelestialBodyFactory.getBody(name).getPVCoordinates(date2,
101 refFrame);
102
103 Assertions.assertEquals(0.0, Vector3D.distance(pRef, pv.getPosition()), 15.0);
104 Assertions.assertEquals(0.0, Vector3D.distance(vRef, pv.getVelocity()), 1.0e-5);
105 }
106 }
107 }
108 }
109
110 @Test
111 void testPO405() {
112
113 Utils.setDataRoot("regular-data");
114 double threshold = 4.0e-11;
115
116
117
118
119
120
121
122
123
124 testPOCoordinate(1969, 6, 1, 6, 8, 2, 28.3804268378833, threshold);
125 testPOCoordinate(1969, 7, 1, 10, 4, 3, 0.2067143944892, threshold);
126 testPOCoordinate(1969, 8, 1, 11, 5, 6, 0.0028060503833, threshold);
127 testPOCoordinate(1969, 9, 1, 11, 8, 4, -0.0026546031502, threshold);
128 testPOCoordinate(1969, 10, 1, 13, 9, 2, 1.3012071081901, threshold);
129
130
131
132
133
134
135 testPOCoordinate(1970, 1, 1, 3, 2, 4, -0.0372587543468, threshold);
136 testPOCoordinate(1970, 2, 1, 12, 11, 4, 0.0000020665428, threshold);
137 testPOCoordinate(1970, 3, 1, 8, 7, 1, 2.7844971346089, threshold);
138 testPOCoordinate(1970, 4, 1, 2, 8, 6, 0.0067657404049, threshold);
139
140
141
142
143 testPOCoordinate(1970, 7, 1, 15, 0, 4, 0.0002194918273, threshold);
144 testPOCoordinate(1970, 8, 1, 7, 12, 5, -0.0037257497269, threshold);
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161 testPOCoordinate(2003, 1, 1, 12, 9, 5, 0.0008047877725, threshold);
162 testPOCoordinate(2003, 2, 1, 14, 0, 1, -0.0000669416537, threshold);
163 testPOCoordinate(2003, 3, 1, 10, 5, 3, -1.4209598221498, threshold);
164 testPOCoordinate(2003, 4, 1, 14, 0, 3, -0.0000007196601, threshold);
165 testPOCoordinate(2003, 5, 1, 7, 12, 3, -4.2775692622201, threshold);
166 testPOCoordinate(2003, 6, 1, 3, 8, 3, 8.5963291192940, threshold);
167 testPOCoordinate(2003, 7, 1, 9, 6, 1, -5.4895120744145, threshold);
168 testPOCoordinate(2003, 8, 1, 4, 5, 2, -3.4742823298269, threshold);
169 testPOCoordinate(2003, 9, 1, 2, 5, 5, -0.0124587966663, threshold);
170 testPOCoordinate(2003, 10, 1, 10, 9, 6, 0.0078155256966, threshold);
171 testPOCoordinate(2003, 11, 1, 10, 7, 3, 4.2838045135189, threshold);
172 testPOCoordinate(2003, 12, 1, 5, 3, 6, -0.0050290663372, threshold);
173 testPOCoordinate(2004, 1, 1, 11, 6, 5, 0.0009776712155, threshold);
174 testPOCoordinate(2004, 2, 1, 2, 7, 4, -0.0175274499718, threshold);
175
176 }
177
178 @Test
179 void testPO406() {
180
181 Utils.setDataRoot("regular-data");
182 double threshold = 2.0e-13;
183
184
185
186
187
188
189
190 testPOCoordinate(2964, 8, 1, 9, 12, 6, -0.0011511788059, threshold);
191 testPOCoordinate(2964, 10, 1, 6, 8, 5, 0.0046432313657, threshold);
192 testPOCoordinate(2964, 12, 1, 10, 11, 5, 0.0090766356095, threshold);
193
194 }
195
196 private void testPOCoordinate(final int year, final int month, final int day,
197 final int targetNumber, final int centerNumber,
198 final int coordinateNumber, final double coordinateValue,
199 final double threshold)
200 {
201 final AbsoluteDate date = new AbsoluteDate(year, month, day, TimeScalesFactory.getTDB());
202 final CelestialBody target = getBody(targetNumber);
203 final CelestialBody center = getBody(centerNumber);
204 if (target != null && center != null) {
205 final PVCoordinates relativePV =
206 new PVCoordinates(center.getPVCoordinates(date, FramesFactory.getICRF()),
207 target.getPVCoordinates(date, FramesFactory.getICRF()));
208 Assertions.assertEquals(coordinateValue, getCoordinate(coordinateNumber, relativePV), threshold);
209 }
210 }
211
212 private CelestialBody getBody(final int number) {
213 switch (number) {
214 case 1 :
215 return CelestialBodyFactory.getMercury();
216 case 2 :
217 return CelestialBodyFactory.getVenus();
218 case 3 :
219 return CelestialBodyFactory.getEarth();
220 case 4 :
221 return CelestialBodyFactory.getMars();
222 case 5 :
223 return CelestialBodyFactory.getJupiter();
224 case 6 :
225 return CelestialBodyFactory.getSaturn();
226 case 7 :
227 return CelestialBodyFactory.getUranus();
228 case 8 :
229 return CelestialBodyFactory.getNeptune();
230 case 9 :
231 return CelestialBodyFactory.getPluto();
232 case 10 :
233 return CelestialBodyFactory.getMoon();
234 case 11 :
235 return CelestialBodyFactory.getSun();
236 case 12 :
237 return CelestialBodyFactory.getSolarSystemBarycenter();
238 case 13 :
239 return CelestialBodyFactory.getEarthMoonBarycenter();
240 default :
241 return null;
242 }
243 }
244
245 public double getCoordinate(final int number, final PVCoordinates pv) {
246 switch (number) {
247 case 1 :
248 return pv.getPosition().getX() / Constants.JPL_SSD_ASTRONOMICAL_UNIT;
249 case 2 :
250 return pv.getPosition().getY() / Constants.JPL_SSD_ASTRONOMICAL_UNIT;
251 case 3 :
252 return pv.getPosition().getZ() / Constants.JPL_SSD_ASTRONOMICAL_UNIT;
253 case 4 :
254 return pv.getVelocity().getX() * Constants.JULIAN_DAY / Constants.JPL_SSD_ASTRONOMICAL_UNIT;
255 case 5 :
256 return pv.getVelocity().getY() * Constants.JULIAN_DAY / Constants.JPL_SSD_ASTRONOMICAL_UNIT;
257 case 6 :
258 return pv.getVelocity().getZ() * Constants.JULIAN_DAY / Constants.JPL_SSD_ASTRONOMICAL_UNIT;
259 default :
260 return Double.NaN;
261 }
262 }
263
264 @Test
265 public void noMercury() {
266 Assertions.assertThrows(OrekitException.class, () -> {
267 Utils.setDataRoot("no-data");
268 CelestialBodyFactory.getMercury();
269 });
270 }
271
272 @Test
273 public void noVenus() {
274 Assertions.assertThrows(OrekitException.class, () -> {
275 Utils.setDataRoot("no-data");
276 CelestialBodyFactory.getVenus();
277 });
278 }
279
280 @Test
281 public void noEarthMoonBarycenter() {
282 Assertions.assertThrows(OrekitException.class, () -> {
283 Utils.setDataRoot("no-data");
284 CelestialBodyFactory.getEarthMoonBarycenter();
285 });
286 }
287
288 @Test
289 public void noMars() {
290 Assertions.assertThrows(OrekitException.class, () -> {
291 Utils.setDataRoot("no-data");
292 CelestialBodyFactory.getMars();
293 });
294 }
295
296 @Test
297 public void noJupiter() {
298 Assertions.assertThrows(OrekitException.class, () -> {
299 Utils.setDataRoot("no-data");
300 CelestialBodyFactory.getJupiter();
301 });
302 }
303
304 @Test
305 public void noSaturn() {
306 Assertions.assertThrows(OrekitException.class, () -> {
307 Utils.setDataRoot("no-data");
308 CelestialBodyFactory.getSaturn();
309 });
310 }
311
312 @Test
313 public void noUranus() {
314 Assertions.assertThrows(OrekitException.class, () -> {
315 Utils.setDataRoot("no-data");
316 CelestialBodyFactory.getUranus();
317 });
318 }
319
320 @Test
321 public void noNeptune() {
322 Assertions.assertThrows(OrekitException.class, () -> {
323 Utils.setDataRoot("no-data");
324 CelestialBodyFactory.getNeptune();
325 });
326 }
327
328 @Test
329 public void noPluto() {
330 Assertions.assertThrows(OrekitException.class, () -> {
331 Utils.setDataRoot("no-data");
332 CelestialBodyFactory.getPluto();
333 });
334 }
335
336 @Test
337 public void noMoon() {
338 Assertions.assertThrows(OrekitException.class, () -> {
339 Utils.setDataRoot("no-data");
340 CelestialBodyFactory.getMoon();
341 });
342 }
343
344 @Test
345 public void noSun() {
346 Assertions.assertThrows(OrekitException.class, () -> {
347 Utils.setDataRoot("no-data");
348 CelestialBodyFactory.getSun();
349 });
350 }
351
352 @Test
353 void testFrameShift() {
354 Utils.setDataRoot("regular-data");
355 final Frame moon = CelestialBodyFactory.getMoon().getBodyOrientedFrame();
356 final Frame earth = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
357 final AbsoluteDate date0 = new AbsoluteDate(1969, 06, 25, TimeScalesFactory.getTDB());
358
359 for (double t = 0; t < Constants.JULIAN_DAY; t += 3600) {
360 final AbsoluteDate date = date0.shiftedBy(t);
361 final Transform transform = earth.getTransformTo(moon, date);
362 for (double dt = -10; dt < 10; dt += 0.125) {
363 final Transform shifted = transform.shiftedBy(dt);
364 final Transform computed = earth.getTransformTo(moon, transform.getDate().shiftedBy(dt));
365 final Transform error = new Transform(computed.getDate(), computed, shifted.getInverse());
366 Assertions.assertEquals(0.0, error.getTranslation().getNorm(), 100.0);
367 Assertions.assertEquals(0.0, error.getVelocity().getNorm(), 20.0);
368 Assertions.assertEquals(0.0, error.getRotation().getAngle(), 4.0e-8);
369 Assertions.assertEquals(0.0, error.getRotationRate().getNorm(), 8.0e-10);
370 }
371 }
372 }
373
374 @Test
375 void testPropagationVsEphemeris() {
376
377 Utils.setDataRoot("regular-data");
378
379
380 final CelestialBody sun = CelestialBodyFactory.getSun();
381 final CelestialBody mercury = CelestialBodyFactory.getMercury();
382 final CelestialBody venus = CelestialBodyFactory.getVenus();
383 final CelestialBody earth = CelestialBodyFactory.getEarth();
384 final CelestialBody mars = CelestialBodyFactory.getMars();
385 final CelestialBody jupiter = CelestialBodyFactory.getJupiter();
386 final CelestialBody saturn = CelestialBodyFactory.getSaturn();
387 final CelestialBody uranus = CelestialBodyFactory.getUranus();
388 final CelestialBody neptune = CelestialBodyFactory.getNeptune();
389 final CelestialBody pluto = CelestialBodyFactory.getPluto();
390
391
392 final AbsoluteDate startingDate = new AbsoluteDate(2000, 1, 2, TimeScalesFactory.getUTC());
393 AbsoluteDate endDate = startingDate.shiftedBy(30 * Constants.JULIAN_DAY);
394
395 final Frame icrf = FramesFactory.getICRF();
396
397
398 double negligibleMu = 1.0e-3;
399 SpacecraftState initialState = new SpacecraftState(new CartesianOrbit(venus.getPVCoordinates(startingDate, icrf),
400 icrf, startingDate, negligibleMu));
401
402
403 final double[][] tol = ToleranceProvider.getDefaultToleranceProvider(1000.).getTolerances(initialState.getOrbit(), OrbitType.CARTESIAN);
404 AbstractIntegrator dop1 = new DormandPrince853Integrator(1.0, 1.0e5, tol[0], tol[1]);
405 NumericalPropagator propag = new NumericalPropagator(dop1);
406 propag.setOrbitType(OrbitType.CARTESIAN);
407 propag.setInitialState(initialState);
408 propag.setMu(negligibleMu);
409
410
411 propag.addForceModel(new BodyAttraction(sun));
412 propag.addForceModel(new BodyAttraction(mercury));
413 propag.addForceModel(new BodyAttraction(earth));
414 propag.addForceModel(new BodyAttraction(mars));
415 propag.addForceModel(new BodyAttraction(jupiter));
416 propag.addForceModel(new BodyAttraction(saturn));
417 propag.addForceModel(new BodyAttraction(uranus));
418 propag.addForceModel(new BodyAttraction(neptune));
419 propag.addForceModel(new BodyAttraction(pluto));
420
421
422 propag.setStepHandler(1000.0, currentState -> {
423
424 Vector3D propagatedP = currentState.getPosition(icrf);
425 Vector3D ephemerisP = venus.getPosition(currentState.getDate(), icrf);
426 Assertions.assertEquals(0, Vector3D.distance(propagatedP, ephemerisP), 1400.0);
427 });
428
429 propag.propagate(startingDate, endDate);
430 }
431
432 private static class BodyAttraction implements ForceModel {
433
434
435 public static final String ATTRACTION_COEFFICIENT_SUFFIX = " attraction coefficient";
436
437
438 private final ParameterDriver parameterDriver;
439
440
441 private final CelestialBody body;
442
443
444
445
446
447
448 public BodyAttraction(final CelestialBody body) {
449 parameterDriver = new ParameterDriver(body.getName() + ATTRACTION_COEFFICIENT_SUFFIX,
450 body.getGM(), 1.0e-5 * body.getGM(),
451 0.0, Double.POSITIVE_INFINITY);
452 this.body = body;
453 }
454
455
456 @Override
457 public boolean dependsOnPositionOnly() {
458 return true;
459 }
460
461
462 @Override
463 public Vector3D acceleration(final SpacecraftState s, final double[] parameters)
464 {
465
466 final double gm = parameters[0];
467
468
469 final Vector3D centralToBody = body.getPosition(s.getDate(), s.getFrame());
470 final Vector3D satToBody = centralToBody.subtract(s.getPosition());
471 final double r2Sat = satToBody.getNormSq();
472
473
474 return new Vector3D(gm / (r2Sat * FastMath.sqrt(r2Sat)), satToBody);
475
476 }
477
478
479 @Override
480 public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
481 final T[] parameters)
482 {
483
484 final T gm = parameters[0];
485
486
487 final FieldVector3D<T> centralToBody = body.getPosition(s.getDate(), s.getFrame());
488 final FieldVector3D<T> satToBody = centralToBody.subtract(s.getPosition());
489 final T r2Sat = satToBody.getNormSq();
490
491
492 return new FieldVector3D<>(r2Sat.multiply(r2Sat.sqrt()).reciprocal().multiply(gm), satToBody);
493
494 }
495
496
497 @Override
498 public Stream<EventDetector> getEventDetectors() {
499 return Stream.empty();
500 }
501
502
503 @Override
504 public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventDetectors(final Field<T> field) {
505 return Stream.empty();
506 }
507
508
509 @Override
510 public List<ParameterDriver> getParametersDrivers() {
511 return Collections.singletonList(parameterDriver);
512 }
513
514 }
515
516 @Test
517 void testKepler() {
518 Utils.setDataRoot("regular-data");
519 AbsoluteDate date = new AbsoluteDate(1969, 06, 28, TimeScalesFactory.getTT());
520 final double au = 149597870691.0;
521 checkKepler(CelestialBodyFactory.getMoon(), CelestialBodyFactory.getEarth(), date, 3.844e8, 0.012);
522 checkKepler(CelestialBodyFactory.getMercury(), CelestialBodyFactory.getSun(), date, 0.387 * au, 4.0e-9);
523 checkKepler(CelestialBodyFactory.getVenus(), CelestialBodyFactory.getSun(), date, 0.723 * au, 8.0e-9);
524 checkKepler(CelestialBodyFactory.getEarth(), CelestialBodyFactory.getSun(), date, 1.000 * au, 2.0e-5);
525 checkKepler(CelestialBodyFactory.getMars(), CelestialBodyFactory.getSun(), date, 1.52 * au, 2.0e-7);
526 checkKepler(CelestialBodyFactory.getJupiter(), CelestialBodyFactory.getSun(), date, 5.20 * au, 2.0e-6);
527 checkKepler(CelestialBodyFactory.getSaturn(), CelestialBodyFactory.getSun(), date, 9.58 * au, 8.0e-7);
528 checkKepler(CelestialBodyFactory.getUranus(), CelestialBodyFactory.getSun(), date, 19.20 * au, 6.0e-7);
529 checkKepler(CelestialBodyFactory.getNeptune(), CelestialBodyFactory.getSun(), date, 30.05 * au, 4.0e-7);
530 checkKepler(CelestialBodyFactory.getPluto(), CelestialBodyFactory.getSun(), date, 39.24 * au, 3.0e-7);
531 }
532
533 private void checkKepler(final PVCoordinatesProvider orbiting, final CelestialBody central,
534 final AbsoluteDate start, final double a, final double epsilon)
535 {
536
537
538 Orbit orbit = new KeplerianOrbit(orbiting.getPVCoordinates(start, central.getInertiallyOrientedFrame()),
539 central.getInertiallyOrientedFrame(), start, central.getGM());
540 KeplerianPropagator propagator = new KeplerianPropagator(orbit);
541 Assertions.assertEquals(a, orbit.getA(), 0.02 * a);
542 double duration = FastMath.min(50 * Constants.JULIAN_DAY, 0.01 * orbit.getKeplerianPeriod());
543
544 double max = 0;
545 for (AbsoluteDate date = start; date.durationFrom(start) < duration; date = date.shiftedBy(duration / 100)) {
546 PVCoordinates ephemPV = orbiting.getPVCoordinates(date, central.getInertiallyOrientedFrame());
547 PVCoordinates keplerPV = propagator.propagate(date).getPVCoordinates();
548 Vector3D error = keplerPV.getPosition().subtract(ephemPV.getPosition());
549 max = FastMath.max(max, error.getNorm());
550 }
551 Assertions.assertTrue(max < epsilon * a);
552 }
553
554
555
556
557
558 @Test
559 void testGetPosition() {
560 Utils.setDataRoot("regular-data");
561
562
563
564
565
566 final Frame j2000 = FramesFactory.getEME2000();
567 final AbsoluteDate date = new AbsoluteDate();
568
569
570 for (int iBody = 1; iBody <= 13; iBody++) {
571
572 final CelestialBody body = getBody(iBody);
573 final double dP = Vector3D.distance(body.getPosition(date, j2000),
574 body.getPosition(date, j2000));
575
576 Assertions.assertEquals(0., dP, 0.);
577 }
578
579
580
581
582
583 final FieldAbsoluteDate<Binary64> fDate = new FieldAbsoluteDate<>(Binary64Field.getInstance());
584
585
586 for (int iBody = 1; iBody <= 13; iBody++) {
587
588 final CelestialBody body = getBody(iBody);
589 final Binary64 dP = FieldVector3D.distance(body.getPVCoordinates(fDate, j2000).getPosition(),
590 body.getPosition(fDate, j2000));
591
592 Assertions.assertEquals(0., dP.getReal(), 0.);
593 }
594 }
595
596 }