1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.bodies;
18
19 import org.hipparchus.geometry.euclidean.threed.Vector3D;
20 import org.hipparchus.util.FastMath;
21 import org.hipparchus.util.MathUtils;
22 import org.hipparchus.util.SinCos;
23 import org.junit.jupiter.api.Assertions;
24 import org.junit.jupiter.api.BeforeEach;
25 import org.junit.jupiter.api.Test;
26 import org.orekit.Utils;
27 import org.orekit.frames.Frame;
28 import org.orekit.frames.FramesFactory;
29 import org.orekit.frames.TopocentricFrame;
30 import org.orekit.time.AbsoluteDate;
31 import org.orekit.utils.IERSConventions;
32
33
34 public class LoxodromeTest {
35 private GeodeticPoint newYork;
36 private GeodeticPoint london;
37 private GeodeticPoint berlin;
38 private GeodeticPoint perth;
39 private GeodeticPoint philadelphia;
40 private OneAxisEllipsoid earth;
41
42 @BeforeEach
43 public void setup() {
44 Utils.setDataRoot("regular-data");
45
46 newYork = new GeodeticPoint(FastMath.toRadians(40.714268), FastMath.toRadians(-74.005974), 0);
47 london = new GeodeticPoint(FastMath.toRadians(51.5), FastMath.toRadians(-0.16667), 0);
48 berlin = new GeodeticPoint(FastMath.toRadians(52.523403), FastMath.toRadians(13.4114), 0);
49 perth = new GeodeticPoint(FastMath.toRadians(-31.952712), FastMath.toRadians(115.8604796), 0);
50 philadelphia = new GeodeticPoint(FastMath.toRadians(39.952330), FastMath.toRadians(-75.16379), 0);
51
52 final Frame ecef = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
53 earth = new OneAxisEllipsoid(6378137, 1. / 298.257223563, ecef);
54 }
55
56
57
58 @Test
59 public void verifyShortDistance() {
60
61 executeTest("new york - philadelphia", newYork, philadelphia, FastMath.toRadians(229.31), 1., 5.);
62 }
63
64 @Test
65 public void verifyLongDistance() {
66 executeTest("new york - london", newYork, london, FastMath.toRadians(78.09), 30., 6_000.);
67 executeTest("berlin - perth", berlin, perth, FastMath.toRadians(132.89), 90., 35_000.);
68 }
69
70 @Test
71 public void verifyParallel() {
72 executeTest("perfect parallel",
73 berlin,
74 new GeodeticPoint(berlin.getLatitude(), berlin.getLongitude() + 0.125, berlin.getAltitude()),
75 MathUtils.SEMI_PI, 0.1, 1.0e-20);
76 }
77
78 @Test
79 public void verifyMeridian() {
80 executeTest("perfect meridian",
81 berlin,
82 new GeodeticPoint(berlin.getLatitude() + 0.125, berlin.getLongitude(), berlin.getAltitude()),
83 0.0, 8.0, 0.003);
84 }
85
86 void executeTest(final String header, final GeodeticPoint start, final GeodeticPoint stop, final double expectedAzimuth, final double numericalError, final double pointError) {
87 final double az = earth.azimuthBetweenPoints(start, stop);
88 Assertions.assertEquals(expectedAzimuth, az, 1e-4);
89
90 final Loxodrome lox = new Loxodrome(start, az, earth);
91 final LoxodromeArc arc = new LoxodromeArc(start, stop, earth);
92
93 final GeodeticPoint expected = highFidelityPointAtDistance(start, az, arc.getDistance());
94 final GeodeticPoint actual = lox.pointAtDistance(arc.getDistance());
95
96 final double error1 = earth.transform(actual).distance(earth.transform(expected));
97 final double error2 = earth.transform(actual).distance(earth.transform(stop));
98
99
100 Assertions.assertTrue(error1 < numericalError,header + " analytic answer not within " + numericalError + "m of numerical answer [error=" + error1 + "m]");
101
102 Assertions.assertTrue(error2 < pointError,header + " computed destination not within " + pointError + "m of actual answer.[error=" + error2 + "m]");
103
104 Assertions.assertEquals(0.0, error1 / arc.getDistance(), 1e-5,header + " accuracy beyond allowable tolerance");
105 }
106
107 GeodeticPoint highFidelityPointAtDistance(final GeodeticPoint start, final double azimuth, final double distance) {
108 GeodeticPoint point = start;
109 double d = distance;
110 final SinCos scAzimuth = FastMath.sinCos(azimuth);
111 final Vector3D step = new Vector3D(scAzimuth.sin(), scAzimuth.cos(), 0);
112 while (d > 0) {
113 final GeodeticPoint tmp = earth.transform(
114 new TopocentricFrame(earth, point, "frame").getStaticTransformTo(earth.getBodyFrame(), AbsoluteDate.ARBITRARY_EPOCH).transformPosition(step),
115 earth.getBodyFrame(),
116 AbsoluteDate.ARBITRARY_EPOCH);
117 d -= 1.;
118 point = tmp;
119 }
120 return point;
121 }
122 }