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