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.random.RandomGenerator;
21 import org.hipparchus.random.Well1024a;
22 import org.hipparchus.random.Well19937a;
23 import org.hipparchus.util.FastMath;
24 import org.hipparchus.util.MathArrays;
25 import org.junit.jupiter.api.Assertions;
26 import org.junit.jupiter.api.Test;
27 import org.orekit.errors.OrekitException;
28 import org.orekit.errors.OrekitMessages;
29 import org.orekit.frames.FactoryManagedFrame;
30 import org.orekit.frames.FramesFactory;
31 import org.orekit.frames.Predefined;
32 import org.orekit.models.earth.ReferenceEllipsoid;
33
34 import java.io.IOException;
35
36
37 public class EllipsoidTest {
38
39 @Test
40 public void testGetters() {
41
42 final Ellipsoid ellipsoid =
43 new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
44
45 Assertions.assertEquals(Predefined.EME2000,
46 ((FactoryManagedFrame) ellipsoid.getFrame()).getFactoryKey());
47 Assertions.assertEquals(1.0, ellipsoid.getA(), 1.0e-15);
48 Assertions.assertEquals(2.0, ellipsoid.getB(), 1.0e-15);
49 Assertions.assertEquals(3.0, ellipsoid.getC(), 1.0e-15);
50
51 }
52
53 @Test
54 public void testPrincipalPlanesIntersections() {
55
56 final Ellipsoid ellipsoid =
57 new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
58
59 final Ellipse xy = ellipsoid.getPlaneSection(Vector3D.ZERO, Vector3D.PLUS_K);
60 Assertions.assertEquals(0, Vector3D.distance(Vector3D.ZERO, xy.getCenter()), 1.0e-15);
61 checkPrincipalAxes(xy, Vector3D.PLUS_J, Vector3D.MINUS_I);
62 Assertions.assertEquals(2.0, xy.getA(), 1.0e-15);
63 Assertions.assertEquals(1.0, xy.getB(), 1.0e-15);
64 Assertions.assertTrue(xy.getFrame() == ellipsoid.getFrame());
65 Assertions.assertEquals(0.0, errorOnEllipsoid(xy, ellipsoid), 1.0e-15);
66 Assertions.assertEquals(0.0, errorOnPlane(xy, Vector3D.ZERO, Vector3D.PLUS_K), 1.0e-15);
67
68 final Ellipse yz = ellipsoid.getPlaneSection(Vector3D.ZERO, Vector3D.PLUS_I);
69 Assertions.assertEquals(0, Vector3D.distance(Vector3D.ZERO, yz.getCenter()), 1.0e-15);
70 checkPrincipalAxes(yz, Vector3D.PLUS_K, Vector3D.MINUS_J);
71 Assertions.assertEquals(3.0, yz.getA(), 1.0e-15);
72 Assertions.assertEquals(2.0, yz.getB(), 1.0e-15);
73 Assertions.assertTrue(yz.getFrame() == ellipsoid.getFrame());
74 Assertions.assertEquals(0.0, errorOnEllipsoid(yz, ellipsoid), 1.0e-15);
75 Assertions.assertEquals(0.0, errorOnPlane(yz, Vector3D.ZERO, Vector3D.PLUS_I), 1.0e-15);
76
77 final Ellipse zx = ellipsoid.getPlaneSection(Vector3D.ZERO, Vector3D.PLUS_J);
78 Assertions.assertEquals(0, Vector3D.distance(Vector3D.ZERO, zx.getCenter()), 1.0e-15);
79 checkPrincipalAxes(zx, Vector3D.PLUS_K, Vector3D.PLUS_I);
80 Assertions.assertEquals(3.0, zx.getA(), 1.0e-15);
81 Assertions.assertEquals(1.0, zx.getB(), 1.0e-15);
82 Assertions.assertTrue(zx.getFrame() == ellipsoid.getFrame());
83 Assertions.assertEquals(0.0, errorOnEllipsoid(zx, ellipsoid), 1.0e-15);
84 Assertions.assertEquals(0.0, errorOnPlane(zx, Vector3D.ZERO, Vector3D.PLUS_J), 1.0e-15);
85
86 }
87
88 @Test
89 public void testNoIntersections() {
90 final Ellipsoid ellipsoid =
91 new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
92 final Ellipse ps = ellipsoid.getPlaneSection(new Vector3D(0, 0, 4), Vector3D.PLUS_K);
93 Assertions.assertNull(ps);
94 }
95
96 @Test
97 public void testSinglePoint() throws IOException {
98 final Ellipsoid ellipsoid =
99 new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
100 final Ellipse ps = ellipsoid.getPlaneSection(new Vector3D(0, 0, 3), Vector3D.PLUS_K);
101 Assertions.assertEquals(0, Vector3D.distance(new Vector3D(0, 0, 3), ps.getCenter()), 1.0e-15);
102 Assertions.assertEquals(0.0, ps.getA(), 1.0e-15);
103 Assertions.assertEquals(0.0, ps.getB(), 1.0e-15);
104 }
105
106 @Test
107 public void testRandomNormalSections() throws IOException {
108 RandomGenerator random = new Well19937a(0x573c54d152aeafe4l);
109 for (int i = 0; i < 100; ++i) {
110 double a = 10 * random.nextDouble();
111 double b = 10 * random.nextDouble();
112 double c = 10 * random.nextDouble();
113 double size = FastMath.max(FastMath.max(a, b), c);
114 final Ellipsoid ellipsoid = new Ellipsoid(FramesFactory.getEME2000(), a, b, c);
115 for (int j = 0; j < 50; ++j) {
116 double phi = FastMath.PI * (random.nextDouble() - 0.5);
117 double lambda = 2 * FastMath.PI * random.nextDouble();
118 double cPhi = FastMath.cos(phi);
119 double sPhi = FastMath.sin(phi);
120 double cLambda = FastMath.cos(lambda);
121 double sLambda = FastMath.sin(lambda);
122 Vector3D surfacePoint = new Vector3D(ellipsoid.getA() * cPhi * cLambda,
123 ellipsoid.getB() * cPhi * sLambda,
124 ellipsoid.getC() * sPhi);
125 Vector3D t1 = new Vector3D(-ellipsoid.getA() * cPhi * sLambda,
126 ellipsoid.getB() * cPhi * cLambda,
127 0).normalize();
128 Vector3D t2 = new Vector3D(-ellipsoid.getA() * sPhi * cLambda,
129 -ellipsoid.getB() * sPhi * sLambda,
130 ellipsoid.getC() * cPhi).normalize();
131 final double azimuth = 2 * FastMath.PI * random.nextDouble();
132 double cAzimuth = FastMath.cos(azimuth);
133 double sAzimuth = FastMath.sin(azimuth);
134 Vector3D tAz = new Vector3D(cAzimuth, t1, sAzimuth, t2);
135
136 final Ellipse ps = ellipsoid.getPlaneSection(surfacePoint, tAz);
137 Assertions.assertEquals(0.0, errorOnEllipsoid(ps, ellipsoid), 1.0e-12 * size);
138 Assertions.assertEquals(0.0, errorOnPlane(ps, surfacePoint, tAz), 1.0e-10 * size);
139 double cos = Vector3D.dotProduct(surfacePoint.subtract(ps.getCenter()), ps.getU()) / ps.getA();
140 double sin = Vector3D.dotProduct(surfacePoint.subtract(ps.getCenter()), ps.getV()) / ps.getB();
141 final Vector3D rebuilt = ps.pointAt(FastMath.atan2(sin, cos));
142 Assertions.assertEquals(0, Vector3D.distance(surfacePoint, rebuilt), 1.0e-11 * size);
143 }
144 }
145 }
146
147 @Test
148 public void testInside() {
149 final Ellipsoid ellipsoid =
150 new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
151 for (double f = -2.0; f <= 2.0; f += 1.0 / 128.0) {
152 final boolean inside = FastMath.abs(f) <= 1.0;
153 Assertions.assertEquals(inside, ellipsoid.isInside(new Vector3D(f * ellipsoid.getA(), 0, 0)));
154 Assertions.assertEquals(inside, ellipsoid.isInside(new Vector3D(0, f * ellipsoid.getB(), 0)));
155 Assertions.assertEquals(inside, ellipsoid.isInside(new Vector3D(0, 0, f * ellipsoid.getC())));
156 }
157 }
158
159 @Test
160 public void testLimb() {
161 final Ellipsoid ellipsoid =
162 new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
163 RandomGenerator random = new Well1024a(0xa69c430a67475af7l);
164 for (int i = 0; i < 5000; ++i) {
165 Vector3D observer = new Vector3D((random.nextDouble() - 0.5) * 5,
166 (random.nextDouble() - 0.5) * 5,
167 (random.nextDouble() - 0.5) * 5);
168 Vector3D outside = new Vector3D((random.nextDouble() - 0.5) * 5,
169 (random.nextDouble() - 0.5) * 5,
170 (random.nextDouble() - 0.5) * 5);
171 if (ellipsoid.isInside(observer)) {
172 try {
173 ellipsoid.pointOnLimb(observer, outside);
174 Assertions.fail("an exception should have been thrown");
175 } catch (OrekitException oe) {
176 Assertions.assertEquals(OrekitMessages.POINT_INSIDE_ELLIPSOID, oe.getSpecifier());
177 }
178 } else {
179 checkOnLimb(ellipsoid, observer, outside);
180 }
181 }
182 }
183
184 private void checkOnLimb(Ellipsoid ellipsoid, Vector3D observer, Vector3D outside) {
185 final Vector3D onLimb = ellipsoid.pointOnLimb(observer, outside);
186 Assertions.assertEquals(0,
187 FastMath.sin(Vector3D.angle(Vector3D.crossProduct(observer, outside),
188 Vector3D.crossProduct(observer, onLimb))),
189 2e-14);
190 final double scaledX = onLimb.getX() / ellipsoid.getA();
191 final double scaledY = onLimb.getY() / ellipsoid.getB();
192 final double scaledZ = onLimb.getZ() / ellipsoid.getC();
193 Assertions.assertEquals(1.0, scaledX * scaledX + scaledY * scaledY + scaledZ * scaledZ, 9e-11);
194 final Vector3D normal = new Vector3D(scaledX / ellipsoid.getA(),
195 scaledY / ellipsoid.getB(),
196 scaledZ / ellipsoid.getC()).normalize();
197 final Vector3D lineOfSight = onLimb.subtract(observer).normalize();
198 Assertions.assertEquals(0.0, Vector3D.dotProduct(normal, lineOfSight), 5e-10);
199 }
200
201
202 @Test
203 public void testIssue639() {
204
205 Vector3D observer = new Vector3D(
206 5621586.021199942, -4496118.751975084, 0.000000008);
207 Vector3D outside = new Vector3D(
208 69159195202.69193, 123014642034.89732, -44866184753.460625);
209 Ellipsoid oae = ReferenceEllipsoid.getWgs84(null);
210
211
212 Vector3D actual = oae.pointOnLimb(observer, outside);
213
214
215 checkOnLimb(oae, observer, outside);
216 Assertions.assertFalse(actual.isNaN(),"" + actual);
217 }
218
219 private void checkPrincipalAxes(Ellipse ps, Vector3D expectedU, Vector3D expectedV) {
220 if (Vector3D.dotProduct(expectedU, ps.getU()) >= 0) {
221 Assertions.assertEquals(0, Vector3D.distance(expectedU, ps.getU()), 1.0e-15);
222 Assertions.assertEquals(0, Vector3D.distance(expectedV, ps.getV()), 1.0e-15);
223 } else {
224 Assertions.assertEquals(0, Vector3D.distance(expectedU.negate(), ps.getU()), 1.0e-15);
225 Assertions.assertEquals(0, Vector3D.distance(expectedV.negate(), ps.getV()), 1.0e-15);
226 }
227 }
228
229 private double errorOnEllipsoid(Ellipse ps, Ellipsoid ellipsoid) {
230 double max = 0;
231 for (double theta = 0; theta < 2 * FastMath.PI; theta += 0.1) {
232 Vector3D p = ps.pointAt(theta);
233 double xOa = p.getX() / ellipsoid.getA();
234 double yOb = p.getY() / ellipsoid.getB();
235 double zOc = p.getZ() / ellipsoid.getC();
236 max = FastMath.max(max, FastMath.abs(MathArrays.linearCombination(xOa, xOa, yOb, yOb, zOc, zOc, 1, -1)));
237 }
238 return max;
239 }
240
241 private double errorOnPlane(Ellipse ps, Vector3D planePoint, Vector3D planeNormal) {
242 double max = 0;
243 for (double theta = 0; theta < 2 * FastMath.PI; theta += 0.1) {
244 Vector3D p = ps.pointAt(theta);
245 max = FastMath.max(max, FastMath.abs(Vector3D.dotProduct(p.subtract(planePoint).normalize(), planeNormal)));
246 }
247 return max;
248 }
249
250 }
251