1   /* Copyright 2002-2025 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.bodies;
18  
19  import java.io.IOException;
20  
21  import org.hipparchus.CalculusFieldElement;
22  import org.hipparchus.Field;
23  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
24  import org.hipparchus.geometry.euclidean.threed.Vector3D;
25  import org.hipparchus.random.RandomGenerator;
26  import org.hipparchus.random.Well1024a;
27  import org.hipparchus.random.Well19937a;
28  import org.hipparchus.util.Binary64Field;
29  import org.hipparchus.util.FastMath;
30  import org.hipparchus.util.MathArrays;
31  import org.junit.jupiter.api.Assertions;
32  import org.junit.jupiter.api.Test;
33  import org.orekit.errors.OrekitException;
34  import org.orekit.errors.OrekitMessages;
35  import org.orekit.frames.FactoryManagedFrame;
36  import org.orekit.frames.FramesFactory;
37  import org.orekit.frames.Predefined;
38  import org.orekit.models.earth.ReferenceEllipsoid;
39  
40  
41  public class EllipsoidTest {
42  
43      @Test
44      void testGetters() {
45  
46          final Ellipsoid ellipsoid =
47                  new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
48  
49          Assertions.assertEquals(Predefined.EME2000,
50                              ((FactoryManagedFrame) ellipsoid.getFrame()).getFactoryKey());
51          Assertions.assertEquals(1.0, ellipsoid.getA(), 1.0e-15);
52          Assertions.assertEquals(2.0, ellipsoid.getB(), 1.0e-15);
53          Assertions.assertEquals(3.0, ellipsoid.getC(), 1.0e-15);
54  
55      }
56  
57      @Test
58      void testPrincipalPlanesIntersections() {
59  
60          final Ellipsoid ellipsoid =
61                  new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
62  
63          final Ellipse xy = ellipsoid.getPlaneSection(Vector3D.ZERO, Vector3D.PLUS_K);
64          Assertions.assertEquals(0, Vector3D.distance(Vector3D.ZERO, xy.getCenter()), 1.0e-15);
65          checkPrincipalAxes(xy, Vector3D.PLUS_J, Vector3D.MINUS_I);
66          Assertions.assertEquals(2.0, xy.getA(), 1.0e-15);
67          Assertions.assertEquals(1.0, xy.getB(), 1.0e-15);
68          Assertions.assertTrue(xy.getFrame() == ellipsoid.getFrame());
69          Assertions.assertEquals(0.0, errorOnEllipsoid(xy, ellipsoid), 1.0e-15);
70          Assertions.assertEquals(0.0, errorOnPlane(xy, Vector3D.ZERO, Vector3D.PLUS_K), 1.0e-15);
71  
72          final Ellipse yz = ellipsoid.getPlaneSection(Vector3D.ZERO, Vector3D.PLUS_I);
73          Assertions.assertEquals(0, Vector3D.distance(Vector3D.ZERO, yz.getCenter()), 1.0e-15);
74          checkPrincipalAxes(yz, Vector3D.PLUS_K, Vector3D.MINUS_J);
75          Assertions.assertEquals(3.0, yz.getA(), 1.0e-15);
76          Assertions.assertEquals(2.0, yz.getB(), 1.0e-15);
77          Assertions.assertTrue(yz.getFrame() == ellipsoid.getFrame());
78          Assertions.assertEquals(0.0, errorOnEllipsoid(yz, ellipsoid), 1.0e-15);
79          Assertions.assertEquals(0.0, errorOnPlane(yz, Vector3D.ZERO, Vector3D.PLUS_I), 1.0e-15);
80  
81          final Ellipse zx = ellipsoid.getPlaneSection(Vector3D.ZERO, Vector3D.PLUS_J);
82          Assertions.assertEquals(0, Vector3D.distance(Vector3D.ZERO, zx.getCenter()), 1.0e-15);
83          checkPrincipalAxes(zx, Vector3D.PLUS_K, Vector3D.PLUS_I);
84          Assertions.assertEquals(3.0, zx.getA(), 1.0e-15);
85          Assertions.assertEquals(1.0, zx.getB(), 1.0e-15);
86          Assertions.assertTrue(zx.getFrame() == ellipsoid.getFrame());
87          Assertions.assertEquals(0.0, errorOnEllipsoid(zx, ellipsoid), 1.0e-15);
88          Assertions.assertEquals(0.0, errorOnPlane(zx, Vector3D.ZERO, Vector3D.PLUS_J), 1.0e-15);
89  
90      }
91  
92      @Test
93      void testFieldPrincipalPlanesIntersections() {
94          doTestFieldPrincipalPlanesIntersections(Binary64Field.getInstance());
95      }
96  
97      private <T extends CalculusFieldElement<T>> void doTestFieldPrincipalPlanesIntersections(final Field<T> field) {
98  
99          final Ellipsoid ellipsoid =
100                 new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
101 
102         final FieldVector3D<T> zero   = FieldVector3D.getZero(field);
103         final FieldVector3D<T> plusI  = FieldVector3D.getPlusI(field);
104         final FieldVector3D<T> plusJ  = FieldVector3D.getPlusJ(field);
105         final FieldVector3D<T> plusK  = FieldVector3D.getPlusK(field);
106         final FieldVector3D<T> minusI = FieldVector3D.getMinusI(field);
107         final FieldVector3D<T> minusJ = FieldVector3D.getMinusJ(field);
108 
109         final FieldEllipse<T> xy = ellipsoid.getPlaneSection(zero, plusK);
110         Assertions.assertEquals(0, FieldVector3D.distance(zero, xy.getCenter()).getReal(), 1.0e-15);
111         checkPrincipalAxes(xy, plusJ, minusI);
112         Assertions.assertEquals(2.0, xy.getA().getReal(), 1.0e-15);
113         Assertions.assertEquals(1.0, xy.getB().getReal(), 1.0e-15);
114         Assertions.assertTrue(xy.getFrame() == ellipsoid.getFrame());
115         Assertions.assertEquals(0.0, errorOnEllipsoid(xy, ellipsoid).getReal(), 1.0e-15);
116         Assertions.assertEquals(0.0, errorOnPlane(xy, zero, plusK).getReal(), 1.0e-15);
117 
118         final FieldEllipse<T> yz = ellipsoid.getPlaneSection(zero, plusI);
119         Assertions.assertEquals(0, FieldVector3D.distance(zero, yz.getCenter()).getReal(), 1.0e-15);
120         checkPrincipalAxes(yz, plusK, minusJ);
121         Assertions.assertEquals(3.0, yz.getA().getReal(), 1.0e-15);
122         Assertions.assertEquals(2.0, yz.getB().getReal(), 1.0e-15);
123         Assertions.assertTrue(yz.getFrame() == ellipsoid.getFrame());
124         Assertions.assertEquals(0.0, errorOnEllipsoid(yz, ellipsoid).getReal(), 1.0e-15);
125         Assertions.assertEquals(0.0, errorOnPlane(yz, zero, plusI).getReal(), 1.0e-15);
126 
127         final FieldEllipse<T> zx = ellipsoid.getPlaneSection(zero, plusJ);
128         Assertions.assertEquals(0, FieldVector3D.distance(zero, zx.getCenter()).getReal(), 1.0e-15);
129         checkPrincipalAxes(zx, plusK, plusI);
130         Assertions.assertEquals(3.0, zx.getA().getReal(), 1.0e-15);
131         Assertions.assertEquals(1.0, zx.getB().getReal(), 1.0e-15);
132         Assertions.assertTrue(zx.getFrame() == ellipsoid.getFrame());
133         Assertions.assertEquals(0.0, errorOnEllipsoid(zx, ellipsoid).getReal(), 1.0e-15);
134         Assertions.assertEquals(0.0, errorOnPlane(zx, zero, plusJ).getReal(), 1.0e-15);
135 
136     }
137 
138     @Test
139     void testNoIntersections() {
140         final Ellipsoid ellipsoid =
141                 new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
142         final Ellipse ps = ellipsoid.getPlaneSection(new Vector3D(0, 0, 4), Vector3D.PLUS_K);
143         Assertions.assertNull(ps);
144     }
145 
146     @Test
147     void testFieldNoIntersections() {
148         doTestFieldNoIntersections(Binary64Field.getInstance());
149     }
150 
151     private <T extends CalculusFieldElement<T>> void doTestFieldNoIntersections(final Field<T> field) {
152         final T zero = field.getZero();
153         final Ellipsoid ellipsoid =
154                         new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
155         final FieldEllipse<T> ps = ellipsoid.getPlaneSection(new FieldVector3D<>(zero, zero, zero.newInstance(4)),
156                                                              FieldVector3D.getPlusK(field));
157         Assertions.assertNull(ps);
158     }
159 
160     @Test
161     void testSinglePoint() throws IOException {
162         final Ellipsoid ellipsoid =
163                 new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
164         final Ellipse ps = ellipsoid.getPlaneSection(new Vector3D(0, 0, 3), Vector3D.PLUS_K);
165         Assertions.assertEquals(0, Vector3D.distance(new Vector3D(0, 0, 3), ps.getCenter()), 1.0e-15);
166         Assertions.assertEquals(0.0, ps.getA(), 1.0e-15);
167         Assertions.assertEquals(0.0, ps.getB(), 1.0e-15);
168     }
169 
170     @Test
171     void testFieldSinglePoint() throws IOException {
172         doTestFieldSinglePoint(Binary64Field.getInstance());
173     }
174 
175     private <T extends CalculusFieldElement<T>> void doTestFieldSinglePoint(final Field<T> field) {
176         final T zero = field.getZero();
177         final Ellipsoid ellipsoid =
178                         new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
179         final FieldEllipse<T> ps = ellipsoid.getPlaneSection(new FieldVector3D<>(zero, zero, zero.newInstance(3)),
180                                                              FieldVector3D.getPlusK(field));
181         Assertions.assertEquals(0,
182                                 FieldVector3D.distance(new FieldVector3D<>(zero, zero, zero.newInstance(3)), ps.getCenter()).getReal(),
183                                 1.0e-15);
184         Assertions.assertEquals(0.0, ps.getA().getReal(), 1.0e-15);
185         Assertions.assertEquals(0.0, ps.getB().getReal(), 1.0e-15);
186     }
187 
188     @Test
189     void testRandomNormalSections() throws IOException {
190         RandomGenerator random = new Well19937a(0x573c54d152aeafe4l);
191         for (int i = 0; i < 100; ++i) {
192             double a = 10 * random.nextDouble();
193             double b = 10 * random.nextDouble();
194             double c = 10 * random.nextDouble();
195             double size = FastMath.max(FastMath.max(a, b), c);
196             final Ellipsoid ellipsoid = new Ellipsoid(FramesFactory.getEME2000(), a, b, c);
197             for (int j = 0; j < 50; ++j) {
198                 double phi     = FastMath.PI * (random.nextDouble() - 0.5);
199                 double lambda  = 2 * FastMath.PI * random.nextDouble();
200                 double cPhi    = FastMath.cos(phi);
201                 double sPhi    = FastMath.sin(phi);
202                 double cLambda = FastMath.cos(lambda);
203                 double sLambda = FastMath.sin(lambda);
204                 Vector3D surfacePoint = new Vector3D(ellipsoid.getA() * cPhi * cLambda,
205                                                      ellipsoid.getB() * cPhi * sLambda,
206                                                      ellipsoid.getC() * sPhi);
207                 Vector3D t1 = new Vector3D(-ellipsoid.getA() * cPhi * sLambda,
208                                             ellipsoid.getB() * cPhi * cLambda,
209                                             0).normalize();
210                 Vector3D t2 = new Vector3D(-ellipsoid.getA() * sPhi * cLambda,
211                                            -ellipsoid.getB() * sPhi * sLambda,
212                                             ellipsoid.getC() * cPhi).normalize();
213                 final double azimuth = 2 * FastMath.PI * random.nextDouble();
214                 double cAzimuth = FastMath.cos(azimuth);
215                 double sAzimuth = FastMath.sin(azimuth);
216                 Vector3D tAz = new Vector3D(cAzimuth, t1, sAzimuth, t2);
217 
218                 final Ellipse ps = ellipsoid.getPlaneSection(surfacePoint, tAz);
219                 Assertions.assertEquals(0.0, errorOnEllipsoid(ps, ellipsoid), 1.0e-12 * size);
220                 Assertions.assertEquals(0.0, errorOnPlane(ps, surfacePoint, tAz), 1.0e-10 * size);
221                 double cos = Vector3D.dotProduct(surfacePoint.subtract(ps.getCenter()), ps.getU()) / ps.getA();
222                 double sin = Vector3D.dotProduct(surfacePoint.subtract(ps.getCenter()), ps.getV()) / ps.getB();
223                 final Vector3D rebuilt = ps.pointAt(FastMath.atan2(sin, cos));
224                 Assertions.assertEquals(0, Vector3D.distance(surfacePoint, rebuilt), 1.0e-11 * size);
225             }
226         }
227     }
228 
229     @Test
230     void testFieldRandomNormalSections() throws IOException {
231         doTestFieldRandomNormalSections(Binary64Field.getInstance());
232     }
233 
234     private <T extends CalculusFieldElement<T>> void doTestFieldRandomNormalSections(final Field<T> field) {
235         final T zero = field.getZero();
236         RandomGenerator random = new Well19937a(0x573c54d152aeafe4l);
237         for (int i = 0; i < 100; ++i) {
238             double a = 10 * random.nextDouble();
239             double b = 10 * random.nextDouble();
240             double c = 10 * random.nextDouble();
241             double size = FastMath.max(FastMath.max(a, b), c);
242             final Ellipsoid ellipsoid = new Ellipsoid(FramesFactory.getEME2000(), a, b, c);
243             for (int j = 0; j < 50; ++j) {
244                 double phi     = FastMath.PI * (random.nextDouble() - 0.5);
245                 double lambda  = 2 * FastMath.PI * random.nextDouble();
246                 double cPhi    = FastMath.cos(phi);
247                 double sPhi    = FastMath.sin(phi);
248                 double cLambda = FastMath.cos(lambda);
249                 double sLambda = FastMath.sin(lambda);
250                 FieldVector3D<T> surfacePoint = new FieldVector3D<>(zero.newInstance(ellipsoid.getA() * cPhi * cLambda),
251                                                                     zero.newInstance(ellipsoid.getB() * cPhi * sLambda),
252                                                                     zero.newInstance(ellipsoid.getC() * sPhi));
253                 FieldVector3D<T> t1 = new FieldVector3D<>(zero.newInstance(-ellipsoid.getA() * cPhi * sLambda),
254                                                           zero.newInstance(ellipsoid.getB() * cPhi * cLambda),
255                                                            zero).normalize();
256                 FieldVector3D<T> t2 = new FieldVector3D<>(zero.newInstance(-ellipsoid.getA() * sPhi * cLambda),
257                                                           zero.newInstance(-ellipsoid.getB() * sPhi * sLambda),
258                                                           zero.newInstance(ellipsoid.getC() * cPhi)).normalize();
259                 final double azimuth = 2 * FastMath.PI * random.nextDouble();
260                 double cAzimuth = FastMath.cos(azimuth);
261                 double sAzimuth = FastMath.sin(azimuth);
262                 FieldVector3D<T> tAz = new FieldVector3D<>(cAzimuth, t1, sAzimuth, t2);
263 
264                 final FieldEllipse<T> ps = ellipsoid.getPlaneSection(surfacePoint, tAz);
265                 Assertions.assertEquals(0.0, errorOnEllipsoid(ps, ellipsoid).getReal(), 1.0e-12 * size);
266                 Assertions.assertEquals(0.0, errorOnPlane(ps, surfacePoint, tAz).getReal(), 1.0e-10 * size);
267                 T cos = FieldVector3D.dotProduct(surfacePoint.subtract(ps.getCenter()), ps.getU()).divide(ps.getA());
268                 T sin = FieldVector3D.dotProduct(surfacePoint.subtract(ps.getCenter()), ps.getV()).divide(ps.getB());
269                 final FieldVector3D<T> rebuilt = ps.pointAt(FastMath.atan2(sin, cos));
270                 Assertions.assertEquals(0, FieldVector3D.distance(surfacePoint, rebuilt).getReal(), 1.0e-11 * size);
271             }
272         }
273     }
274 
275     @Test
276     void testInside() {
277         final Ellipsoid ellipsoid =
278                         new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
279         for (double f = -2.0; f <= 2.0; f += 1.0 / 128.0) {
280             final boolean inside = FastMath.abs(f) <= 1.0;
281             Assertions.assertEquals(inside, ellipsoid.isInside(new Vector3D(f * ellipsoid.getA(), 0, 0)));
282             Assertions.assertEquals(inside, ellipsoid.isInside(new Vector3D(0, f * ellipsoid.getB(), 0)));
283             Assertions.assertEquals(inside, ellipsoid.isInside(new Vector3D(0, 0, f * ellipsoid.getC())));
284         }
285     }
286 
287     @Test
288     void testFieldInside() {
289         doTestFieldInside(Binary64Field.getInstance());
290     }
291 
292     private <T extends CalculusFieldElement<T>> void doTestFieldInside(final Field<T> field) {
293         final T zero = field.getZero();
294         final Ellipsoid ellipsoid =
295                         new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
296         for (double f = -2.0; f <= 2.0; f += 1.0 / 128.0) {
297             final boolean inside = FastMath.abs(f) <= 1.0;
298             Assertions.assertEquals(inside, ellipsoid.isInside(new FieldVector3D<>(zero.newInstance(f * ellipsoid.getA()), zero, zero)));
299             Assertions.assertEquals(inside, ellipsoid.isInside(new FieldVector3D<>(zero, zero.newInstance(f * ellipsoid.getB()), zero)));
300             Assertions.assertEquals(inside, ellipsoid.isInside(new FieldVector3D<>(zero, zero, zero.newInstance(f * ellipsoid.getC()))));
301         }
302     }
303 
304     @Test
305     void testLimb() {
306         final Ellipsoid ellipsoid =
307                         new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
308         RandomGenerator random = new Well1024a(0xa69c430a67475af7l);
309         for (int i = 0; i < 5000; ++i) {
310             Vector3D observer = new Vector3D((random.nextDouble() - 0.5) * 5,
311                                              (random.nextDouble() - 0.5) * 5,
312                                              (random.nextDouble() - 0.5) * 5);
313             Vector3D outside  = new Vector3D((random.nextDouble() - 0.5) * 5,
314                                              (random.nextDouble() - 0.5) * 5,
315                                              (random.nextDouble() - 0.5) * 5);
316             if (ellipsoid.isInside(observer)) {
317                 try {
318                     ellipsoid.pointOnLimb(observer, outside);
319                     Assertions.fail("an exception should have been thrown");
320                 } catch (OrekitException oe) {
321                     Assertions.assertEquals(OrekitMessages.POINT_INSIDE_ELLIPSOID, oe.getSpecifier());
322                 }
323             } else {
324                 checkOnLimb(ellipsoid, observer, outside);
325             }
326         }
327     }
328 
329     @Test
330     void testFieldLimb() {
331         doTestFieldLimb(Binary64Field.getInstance());
332     }
333 
334     private <T extends CalculusFieldElement<T>> void doTestFieldLimb(final Field<T> field) {
335         final T zero = field.getZero();
336         final Ellipsoid ellipsoid =
337                         new Ellipsoid(FramesFactory.getEME2000(), 1, 2, 3);
338         RandomGenerator random = new Well1024a(0xa69c430a67475af7l);
339         for (int i = 0; i < 5000; ++i) {
340             FieldVector3D<T> observer = new FieldVector3D<>(zero.newInstance((random.nextDouble() - 0.5) * 5),
341                                                             zero.newInstance((random.nextDouble() - 0.5) * 5),
342                                                             zero.newInstance((random.nextDouble() - 0.5) * 5));
343             FieldVector3D<T> outside  = new FieldVector3D<>(zero.newInstance((random.nextDouble() - 0.5) * 5),
344                                                             zero.newInstance((random.nextDouble() - 0.5) * 5),
345                                                             zero.newInstance((random.nextDouble() - 0.5) * 5));
346             if (ellipsoid.isInside(observer)) {
347                 try {
348                     ellipsoid.pointOnLimb(observer, outside);
349                     Assertions.fail("an exception should have been thrown");
350                 } catch (OrekitException oe) {
351                     Assertions.assertEquals(OrekitMessages.POINT_INSIDE_ELLIPSOID, oe.getSpecifier());
352                 }
353             } else {
354                 checkOnLimb(ellipsoid, observer, outside);
355             }
356         }
357     }
358 
359     /** A test case where x >> y and y << 1, which is stressing numerically. */
360     @Test
361     void testIssue639() {
362         // setup
363         Vector3D observer = new Vector3D(
364                 5621586.021199942, -4496118.751975084, 0.000000008);
365         Vector3D outside = new Vector3D(
366                 69159195202.69193, 123014642034.89732, -44866184753.460625);
367         Ellipsoid oae = ReferenceEllipsoid.getWgs84(null);
368 
369         // action
370         Vector3D actual = oae.pointOnLimb(observer, outside);
371 
372         // verify
373         checkOnLimb(oae, observer, outside);
374         Assertions.assertFalse(actual.isNaN(),"" + actual);
375     }
376 
377     @Test
378     void testFieldIssue639() {
379         doTestFieldIssue639(Binary64Field.getInstance());
380     }
381 
382     private <T extends CalculusFieldElement<T>> void doTestFieldIssue639(final Field<T> field) {
383         final T zero = field.getZero();
384         // setup
385         FieldVector3D<T> observer = new FieldVector3D<>(zero.newInstance(5621586.021199942),
386                                                         zero.newInstance(-4496118.751975084),
387                                                         zero.newInstance(0.000000008));
388         FieldVector3D<T>  outside = new FieldVector3D<>(zero.newInstance(69159195202.69193),
389                                                         zero.newInstance(123014642034.89732),
390                                                         zero.newInstance(-44866184753.460625));
391         Ellipsoid oae = ReferenceEllipsoid.getWgs84(null);
392 
393         // action
394         FieldVector3D<T>  actual = oae.pointOnLimb(observer, outside);
395 
396         // verify
397         checkOnLimb(oae, observer, outside);
398         Assertions.assertFalse(actual.isNaN(),"" + actual);
399     }
400 
401     private void checkPrincipalAxes(Ellipse ps, Vector3D expectedU, Vector3D expectedV) {
402         if (Vector3D.dotProduct(expectedU, ps.getU()) >= 0) {
403             Assertions.assertEquals(0, Vector3D.distance(expectedU,  ps.getU()), 1.0e-15);
404             Assertions.assertEquals(0, Vector3D.distance(expectedV,  ps.getV()), 1.0e-15);
405         } else {
406             Assertions.assertEquals(0, Vector3D.distance(expectedU.negate(), ps.getU()), 1.0e-15);
407             Assertions.assertEquals(0, Vector3D.distance(expectedV.negate(), ps.getV()), 1.0e-15);
408         }
409     }
410 
411     private <T extends CalculusFieldElement<T>> void checkPrincipalAxes(FieldEllipse<T> ps, FieldVector3D<T> expectedU, FieldVector3D<T> expectedV) {
412         if (FieldVector3D.dotProduct(expectedU, ps.getU()).getReal() >= 0) {
413             Assertions.assertEquals(0, FieldVector3D.distance(expectedU,  ps.getU()).getReal(), 1.0e-15);
414             Assertions.assertEquals(0, FieldVector3D.distance(expectedV,  ps.getV()).getReal(), 1.0e-15);
415         } else {
416             Assertions.assertEquals(0, FieldVector3D.distance(expectedU.negate(), ps.getU()).getReal(), 1.0e-15);
417             Assertions.assertEquals(0, FieldVector3D.distance(expectedV.negate(), ps.getV()).getReal(), 1.0e-15);
418         }
419     }
420 
421     private double errorOnEllipsoid(Ellipse ps, Ellipsoid ellipsoid) {
422         double max = 0;
423         for (double theta = 0; theta < 2 * FastMath.PI; theta += 0.1) {
424             Vector3D p = ps.pointAt(theta);
425             double xOa = p.getX() / ellipsoid.getA();
426             double yOb = p.getY() / ellipsoid.getB();
427             double zOc = p.getZ() / ellipsoid.getC();
428             max = FastMath.max(max, FastMath.abs(MathArrays.linearCombination(xOa, xOa, yOb, yOb, zOc, zOc, 1, -1)));
429         }
430         return max;
431     }
432 
433     private <T extends CalculusFieldElement<T>> T errorOnEllipsoid(FieldEllipse<T> ps, Ellipsoid ellipsoid) {
434         T one = ps.getA().getField().getOne();
435         T max = ps.getA().getField().getZero();
436         for (double theta = 0; theta < 2 * FastMath.PI; theta += 0.1) {
437             FieldVector3D<T> p = ps.pointAt(max.newInstance(theta));
438             T xOa = p.getX().divide(ellipsoid.getA());
439             T yOb = p.getY().divide(ellipsoid.getB());
440             T zOc = p.getZ().divide(ellipsoid.getC());
441             max = FastMath.max(max, FastMath.abs(max.linearCombination(xOa, xOa, yOb, yOb, zOc, zOc, one, one.negate())));
442         }
443         return max;
444     }
445 
446     private double errorOnPlane(Ellipse ps, Vector3D planePoint, Vector3D planeNormal) {
447         double max = 0;
448         for (double theta = 0; theta < 2 * FastMath.PI; theta += 0.1) {
449             Vector3D p = ps.pointAt(theta);
450             max = FastMath.max(max, FastMath.abs(Vector3D.dotProduct(p.subtract(planePoint).normalize(), planeNormal)));
451         }
452         return max;
453     }
454 
455     private <T extends CalculusFieldElement<T>> T errorOnPlane(FieldEllipse<T> ps, FieldVector3D<T> planePoint, FieldVector3D<T> planeNormal) {
456         T max = ps.getA().getField().getZero();
457         for (double theta = 0; theta < 2 * FastMath.PI; theta += 0.1) {
458             FieldVector3D<T> p = ps.pointAt(max.newInstance(theta));
459             max = FastMath.max(max, FastMath.abs(FieldVector3D.dotProduct(p.subtract(planePoint).normalize(), planeNormal)));
460         }
461         return max;
462     }
463 
464     private void checkOnLimb(Ellipsoid ellipsoid, Vector3D observer, Vector3D outside) {
465         final Vector3D onLimb = ellipsoid.pointOnLimb(observer, outside);
466         Assertions.assertEquals(0,
467                             FastMath.sin(Vector3D.angle(Vector3D.crossProduct(observer, outside),
468                                                         Vector3D.crossProduct(observer, onLimb))),
469                             2e-14);
470         final double scaledX = onLimb.getX() / ellipsoid.getA();
471         final double scaledY = onLimb.getY() / ellipsoid.getB();
472         final double scaledZ = onLimb.getZ() / ellipsoid.getC();
473         Assertions.assertEquals(1.0, scaledX * scaledX + scaledY * scaledY + scaledZ * scaledZ, 9e-11);
474         final Vector3D normal = new Vector3D(scaledX / ellipsoid.getA(),
475                                              scaledY / ellipsoid.getB(),
476                                              scaledZ / ellipsoid.getC()).normalize();
477         final Vector3D lineOfSight = onLimb.subtract(observer).normalize();
478         Assertions.assertEquals(0.0, Vector3D.dotProduct(normal, lineOfSight), 5e-10);
479     }
480 
481     private <T extends CalculusFieldElement<T>> void checkOnLimb(Ellipsoid ellipsoid, FieldVector3D<T> observer, FieldVector3D<T> outside) {
482         final FieldVector3D<T> onLimb = ellipsoid.pointOnLimb(observer, outside);
483         Assertions.assertEquals(0,
484                                 FastMath.sin(FieldVector3D.angle(FieldVector3D.crossProduct(observer, outside),
485                                                                  FieldVector3D.crossProduct(observer, onLimb))).getReal(),
486                                 2e-14);
487         final T scaledX = onLimb.getX().divide(ellipsoid.getA());
488         final T scaledY = onLimb.getY().divide(ellipsoid.getB());
489         final T scaledZ = onLimb.getZ().divide(ellipsoid.getC());
490         Assertions.assertEquals(1.0,
491                                 scaledX.multiply(scaledX).add(scaledY.multiply(scaledY)).add(scaledZ.multiply(scaledZ)).getReal(),
492                                 9e-11);
493         final FieldVector3D<T> normal = new FieldVector3D<>(scaledX.divide(ellipsoid.getA()),
494                                                             scaledY.divide(ellipsoid.getB()),
495                                                             scaledZ.divide(ellipsoid.getC())).normalize();
496         final FieldVector3D<T> lineOfSight = onLimb.subtract(observer).normalize();
497         Assertions.assertEquals(0.0, FieldVector3D.dotProduct(normal, lineOfSight).getReal(), 5e-10);
498     }
499 
500 }
501