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