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