1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.bodies;
18
19 import java.util.ArrayList;
20 import java.util.List;
21 import java.util.stream.DoubleStream;
22 import java.util.stream.Stream;
23
24 import org.hipparchus.CalculusFieldElement;
25 import org.hipparchus.Field;
26 import org.hipparchus.analysis.UnivariateFunction;
27 import org.hipparchus.analysis.differentiation.DSFactory;
28 import org.hipparchus.analysis.differentiation.DerivativeStructure;
29 import org.hipparchus.analysis.differentiation.FiniteDifferencesDifferentiator;
30 import org.hipparchus.analysis.differentiation.UnivariateDifferentiableFunction;
31 import org.hipparchus.analysis.differentiation.UnivariateDerivative2;
32 import org.hipparchus.geometry.euclidean.oned.Vector1D;
33 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
34 import org.hipparchus.geometry.euclidean.threed.Line;
35 import org.hipparchus.geometry.euclidean.threed.Vector3D;
36 import org.hipparchus.geometry.euclidean.twod.Vector2D;
37 import org.hipparchus.random.SobolSequenceGenerator;
38 import org.hipparchus.util.Binary64;
39 import org.hipparchus.util.Binary64Field;
40 import org.hipparchus.util.FastMath;
41 import org.hipparchus.util.MathUtils;
42 import org.junit.jupiter.api.Assertions;
43 import org.junit.jupiter.api.BeforeEach;
44 import org.junit.jupiter.api.Test;
45 import org.orekit.Utils;
46 import org.orekit.errors.OrekitException;
47 import org.orekit.frames.Frame;
48 import org.orekit.frames.FramesFactory;
49 import org.orekit.orbits.CircularOrbit;
50 import org.orekit.orbits.EquinoctialOrbit;
51 import org.orekit.orbits.Orbit;
52 import org.orekit.orbits.PositionAngleType;
53 import org.orekit.time.AbsoluteDate;
54 import org.orekit.time.DateComponents;
55 import org.orekit.time.FieldAbsoluteDate;
56 import org.orekit.time.TimeComponents;
57 import org.orekit.time.TimeInterpolator;
58 import org.orekit.time.TimeScalesFactory;
59 import org.orekit.utils.CartesianDerivativesFilter;
60 import org.orekit.utils.Constants;
61 import org.orekit.utils.IERSConventions;
62 import org.orekit.utils.PVCoordinates;
63 import org.orekit.utils.PVCoordinatesProvider;
64 import org.orekit.utils.TimeStampedPVCoordinates;
65 import org.orekit.utils.TimeStampedPVCoordinatesHermiteInterpolator;
66
67
68 class OneAxisEllipsoidTest {
69
70 @Test
71 void testStandard() {
72 checkCartesianToEllipsoidic(6378137.0, 1.0 / 298.257222101,
73 4637885.347, 121344.608, 4362452.869,
74 0.026157811533131, 0.757987116290729, 260.455572965555);
75 }
76
77 @Test
78 void testLongitudeZero() {
79 checkCartesianToEllipsoidic(6378137.0, 1.0 / 298.257222101,
80 6378400.0, 0, 6379000.0,
81 0.0, 0.787815771252351, 2653416.77864152);
82 }
83
84 @Test
85 void testLongitudePi() {
86 checkCartesianToEllipsoidic(6378137.0, 1.0 / 298.257222101,
87 -6379999.0, 0, 6379000.0,
88 3.14159265358979, 0.787690146758403, 2654544.7767725);
89 }
90
91 @Test
92 void testNorthPole() {
93 checkCartesianToEllipsoidic(6378137.0, 1.0 / 298.257222101,
94 0.0, 0.0, 7000000.0,
95 0.0, 1.57079632679490, 643247.685859644);
96 }
97
98 @Test
99 void testEquator() {
100 checkCartesianToEllipsoidic(6378137.0, 1.0 / 298.257222101,
101 6379888.0, 6377000.0, 0.0,
102 0.785171775899913, 0.0, 2642345.24279301);
103 }
104
105 @Test
106 void testNoFlattening() {
107 final double r = 7000000.0;
108 final double lambda = 2.345;
109 final double phi = -1.23;
110 final double cL = FastMath.cos(lambda);
111 final double sL = FastMath.sin(lambda);
112 final double cH = FastMath.cos(phi);
113 final double sH = FastMath.sin(phi);
114 checkCartesianToEllipsoidic(6378137.0, 0,
115 r * cL * cH, r * sL * cH, r * sH,
116 lambda, phi, r - 6378137.0);
117 }
118
119 @Test
120 void testNoFlatteningPolar() {
121 final double r = 1000.0;
122 final double h = 100;
123 checkCartesianToEllipsoidic(r, 0,
124 0, 0, r + h,
125 0, 0.5 * FastMath.PI, h);
126 checkCartesianToEllipsoidic(r, 0,
127 0, 0, r - h,
128 0, 0.5 * FastMath.PI, -h);
129 checkCartesianToEllipsoidic(r, 0,
130 0, 0, -(r + h),
131 0, -0.5 * FastMath.PI, h);
132 checkCartesianToEllipsoidic(r, 0,
133 0, 0, -(r - h),
134 0, -0.5 * FastMath.PI, -h);
135 }
136
137 @Test
138 void testOnSurface() {
139 Vector3D surfacePoint = new Vector3D(-1092200.775949484,
140 -3944945.7282234835,
141 4874931.946956173);
142 OneAxisEllipsoid earthShape = new OneAxisEllipsoid(6378136.460, 1 / 298.257222101,
143 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
144 GeodeticPoint gp = earthShape.transform(surfacePoint, earthShape.getBodyFrame(),
145 AbsoluteDate.J2000_EPOCH);
146 Vector3D rebuilt = earthShape.transform(gp);
147 Assertions.assertEquals(0, rebuilt.distance(surfacePoint), 3.0e-9);
148 }
149
150 @Test
151 void testInside3Roots() {
152 checkCartesianToEllipsoidic(6378137.0, 1.0 / 298.257,
153 9219.0, -5322.0, 6056743.0,
154 5.75963470503781, 1.56905114598949, -300000.009586231);
155 }
156
157 @Test
158 void testInsideLessThan3Roots() {
159 checkCartesianToEllipsoidic(6378137.0, 1.0 / 298.257,
160 1366863.0, -789159.0, -5848.988,
161 -0.523598928689, -0.00380885831963, -4799808.27951);
162 }
163
164 @Test
165 void testOutside() {
166 checkCartesianToEllipsoidic(6378137.0, 1.0 / 298.257,
167 5722966.0, -3304156.0, -24621187.0,
168 5.75958652642615, -1.3089969725151, 19134410.3342696);
169 }
170
171 @Test
172 void testGeoCar() {
173 OneAxisEllipsoid model =
174 new OneAxisEllipsoid(6378137.0, 1.0 / 298.257222101,
175 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
176 GeodeticPoint nsp =
177 new GeodeticPoint(0.852479154923577, 0.0423149994747243, 111.6);
178 Vector3D p = model.transform(nsp);
179 Assertions.assertEquals(4201866.69291890, p.getX(), 1.0e-6);
180 Assertions.assertEquals(177908.184625686, p.getY(), 1.0e-6);
181 Assertions.assertEquals(4779203.64408617, p.getZ(), 1.0e-6);
182 }
183
184 @Test
185 void testGroundProjectionPosition() {
186 OneAxisEllipsoid model =
187 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
188 Constants.WGS84_EARTH_FLATTENING,
189 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
190
191 TimeStampedPVCoordinates initPV =
192 new TimeStampedPVCoordinates(AbsoluteDate.J2000_EPOCH.shiftedBy(584.),
193 new Vector3D(3220103., 69623., 6449822.),
194 new Vector3D(6414.7, -2006., -3180.),
195 Vector3D.ZERO);
196 Frame eme2000 = FramesFactory.getEME2000();
197 Orbit orbit = new EquinoctialOrbit(initPV, eme2000, Constants.EIGEN5C_EARTH_MU);
198
199 for (double dt = 0; dt < 3600.0; dt += 60.0) {
200
201 AbsoluteDate date = orbit.getDate().shiftedBy(dt);
202 Vector3D pos = orbit.getPosition(date, eme2000);
203 Vector3D groundPV = model.projectToGround(pos, date, eme2000);
204 Vector3D groundP = model.projectToGround(pos, date, eme2000);
205
206
207 Assertions.assertEquals(model.transform(pos, eme2000, date).getLatitude(),
208 model.transform(groundPV, eme2000, date).getLatitude(),
209 1.0e-10);
210 Assertions.assertEquals(model.transform(pos, eme2000, date).getLongitude(),
211 model.transform(groundPV, eme2000, date).getLongitude(),
212 1.0e-10);
213 Assertions.assertEquals(0.0, Vector3D.distance(groundP, groundPV), 1.0e-15 * groundP.getNorm());
214
215 }
216
217 }
218
219 @Test
220 void testGroundProjectionDerivatives() {
221 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
222 Frame eme2000 = FramesFactory.getEME2000();
223 OneAxisEllipsoid model =
224 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
225 Constants.WGS84_EARTH_FLATTENING,
226 itrf);
227
228 TimeStampedPVCoordinates initPV =
229 new TimeStampedPVCoordinates(AbsoluteDate.J2000_EPOCH.shiftedBy(584.),
230 new Vector3D(3220103., 69623., 6449822.),
231 new Vector3D(6414.7, -2006., -3180.),
232 Vector3D.ZERO);
233 Orbit orbit = new EquinoctialOrbit(initPV, eme2000, Constants.EIGEN5C_EARTH_MU);
234
235 double[] errors = derivativesErrors(orbit, orbit.getDate(), eme2000, model);
236 Assertions.assertEquals(0, errors[0], 1.0e-16);
237 Assertions.assertEquals(0, errors[1], 2.0e-12);
238 Assertions.assertEquals(0, errors[2], 2.0e-4);
239
240 }
241
242 @Test
243 void testGroundToGroundIssue181() {
244 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
245 Frame eme2000 = FramesFactory.getEME2000();
246 OneAxisEllipsoid model =
247 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
248 Constants.WGS84_EARTH_FLATTENING,
249 itrf);
250
251 TimeStampedPVCoordinates initPV =
252 new TimeStampedPVCoordinates(AbsoluteDate.J2000_EPOCH.shiftedBy(584.),
253 new Vector3D(3220103., 69623., 6449822.),
254 new Vector3D(6414.7, -2006., -3180.),
255 Vector3D.ZERO);
256 TimeStampedPVCoordinates body = itrf.getTransformTo(eme2000, initPV.getDate()).transformPVCoordinates(initPV);
257 TimeStampedPVCoordinates ground1 = model.projectToGround(body, itrf);
258 TimeStampedPVCoordinates ground2 = model.projectToGround(ground1, itrf);
259 Assertions.assertEquals(0.0, Vector3D.distance(ground1.getPosition(), ground2.getPosition()), 1.0e-12);
260 Assertions.assertEquals(0.0, Vector3D.distance(ground1.getVelocity(), ground2.getVelocity()), 2.0e-12);
261 Assertions.assertEquals(0.0, Vector3D.distance(ground1.getAcceleration(), ground2.getAcceleration()), 1.0e-12);
262
263 }
264
265 private double[] derivativesErrors(PVCoordinatesProvider provider, AbsoluteDate date, Frame frame,
266 OneAxisEllipsoid model) {
267 List<TimeStampedPVCoordinates> pvList = new ArrayList<TimeStampedPVCoordinates>();
268 List<TimeStampedPVCoordinates> groundPVList = new ArrayList<TimeStampedPVCoordinates>();
269 for (double dt = -0.25; dt <= 0.25; dt += 0.125) {
270 TimeStampedPVCoordinates shiftedPV = provider.getPVCoordinates(date.shiftedBy(dt), frame);
271 Vector3D p = model.projectToGround(shiftedPV.getPosition(), shiftedPV.getDate(), frame);
272 pvList.add(shiftedPV);
273 groundPVList.add(new TimeStampedPVCoordinates(shiftedPV.getDate(),
274 p, Vector3D.ZERO, Vector3D.ZERO));
275 }
276
277
278 final TimeInterpolator<TimeStampedPVCoordinates> interpolator =
279 new TimeStampedPVCoordinatesHermiteInterpolator(pvList.size(), CartesianDerivativesFilter.USE_P);
280
281 final TimeInterpolator<TimeStampedPVCoordinates> interpolatorGround =
282 new TimeStampedPVCoordinatesHermiteInterpolator(groundPVList.size(), CartesianDerivativesFilter.USE_P);
283
284 TimeStampedPVCoordinates computed = model.projectToGround(interpolator.interpolate(date, pvList), frame);
285 TimeStampedPVCoordinates reference = interpolatorGround.interpolate(date, groundPVList);
286
287 TimeStampedPVCoordinates pv0 = provider.getPVCoordinates(date, frame);
288 Vector3D p0 = pv0.getPosition();
289 Vector3D v0 = pv0.getVelocity();
290 Vector3D a0 = pv0.getAcceleration();
291
292 return new double[] {
293 Vector3D.distance(computed.getPosition(), reference.getPosition()) / p0.getNorm(),
294 Vector3D.distance(computed.getVelocity(), reference.getVelocity()) / v0.getNorm(),
295 Vector3D.distance(computed.getAcceleration(), reference.getAcceleration()) / a0.getNorm(),
296 };
297
298 }
299
300 @Test
301 void testGroundProjectionTaylor() {
302 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
303 Frame eme2000 = FramesFactory.getEME2000();
304 OneAxisEllipsoid model =
305 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
306 Constants.WGS84_EARTH_FLATTENING,
307 itrf);
308
309 TimeStampedPVCoordinates initPV =
310 new TimeStampedPVCoordinates(AbsoluteDate.J2000_EPOCH.shiftedBy(584.),
311 new Vector3D(3220103., 69623., 6449822.),
312 new Vector3D(6414.7, -2006., -3180.),
313 Vector3D.ZERO);
314 Orbit orbit = new EquinoctialOrbit(initPV, eme2000, Constants.EIGEN5C_EARTH_MU);
315
316 TimeStampedPVCoordinates pv0 = orbit.getPVCoordinates(orbit.getDate(), model.getBodyFrame());
317 PVCoordinatesProvider groundTaylor =
318 model.projectToGround(pv0, model.getBodyFrame()).toTaylorProvider(model.getBodyFrame());
319
320 TimeStampedPVCoordinates g0 = groundTaylor.getPVCoordinates(orbit.getDate(), model.getBodyFrame());
321 Vector3D zenith = pv0.getPosition().subtract(g0.getPosition()).normalize();
322 Vector3D acrossTrack = Vector3D.crossProduct(zenith, g0.getVelocity()).normalize();
323 Vector3D alongTrack = Vector3D.crossProduct(acrossTrack, zenith).normalize();
324 for (double dt = -1; dt < 1; dt += 0.01) {
325 AbsoluteDate date = orbit.getDate().shiftedBy(dt);
326 Vector3D taylorP = groundTaylor.getPosition(date, model.getBodyFrame());
327 Vector3D refP = model.projectToGround(orbit.getPosition(date, model.getBodyFrame()),
328 date, model.getBodyFrame());
329 Vector3D delta = taylorP.subtract(refP);
330 Assertions.assertEquals(0.0, Vector3D.dotProduct(delta, alongTrack), 0.0015);
331 Assertions.assertEquals(0.0, Vector3D.dotProduct(delta, acrossTrack), 0.0007);
332 Assertions.assertEquals(0.0, Vector3D.dotProduct(delta, zenith), 0.00002);
333 }
334
335 }
336
337 @Test
338 void testLineIntersection() {
339 AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
340 Frame frame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
341
342 OneAxisEllipsoid model = new OneAxisEllipsoid(100.0, 0.9, frame);
343 Vector3D point = new Vector3D(0.0, 93.7139699, 3.5930796);
344 Vector3D direction = new Vector3D(0.0, 1.0, 1.0);
345 Line line = new Line(point, point.add(direction), 1.0e-10);
346 GeodeticPoint gp = model.getIntersectionPoint(line, point, frame, date);
347 Assertions.assertEquals(gp.getAltitude(), 0.0, 1.0e-12);
348 Assertions.assertTrue(line.contains(model.transform(gp)));
349
350 model = new OneAxisEllipsoid(100.0, 0.9, frame);
351 point = new Vector3D(0.0, -93.7139699, -3.5930796);
352 direction = new Vector3D(0.0, -1.0, -1.0);
353 line = new Line(point, point.add(direction), 1.0e-10).revert();
354 gp = model.getIntersectionPoint(line, point, frame, date);
355 Assertions.assertTrue(line.contains(model.transform(gp)));
356
357 model = new OneAxisEllipsoid(100.0, 0.9, frame);
358 point = new Vector3D(0.0, -93.7139699, 3.5930796);
359 direction = new Vector3D(0.0, -1.0, 1.0);
360 line = new Line(point, point.add(direction), 1.0e-10);
361 gp = model.getIntersectionPoint(line, point, frame, date);
362 Assertions.assertTrue(line.contains(model.transform(gp)));
363
364 model = new OneAxisEllipsoid(100.0, 0.9, frame);
365 point = new Vector3D(-93.7139699, 0.0, 3.5930796);
366 direction = new Vector3D(-1.0, 0.0, 1.0);
367 line = new Line(point, point.add(direction), 1.0e-10);
368 gp = model.getIntersectionPoint(line, point, frame, date);
369 Assertions.assertTrue(line.contains(model.transform(gp)));
370 Assertions.assertFalse(line.contains(new Vector3D(0, 0, 7000000)));
371
372 point = new Vector3D(0.0, 0.0, 110);
373 direction = new Vector3D(0.0, 0.0, 1.0);
374 line = new Line(point, point.add(direction), 1.0e-10);
375 gp = model.getIntersectionPoint(line, point, frame, date);
376 Assertions.assertEquals(gp.getLatitude(), FastMath.PI/2, 1.0e-12);
377
378 point = new Vector3D(0.0, 110, 0);
379 direction = new Vector3D(0.0, 1.0, 0.0);
380 line = new Line(point, point.add(direction), 1.0e-10);
381 gp = model.getIntersectionPoint(line, point, frame, date);
382 Assertions.assertEquals(gp.getLatitude(), 0, 1.0e-12);
383
384 }
385
386 @Test
387 void testNoLineIntersection() {
388 AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
389 Frame frame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
390 OneAxisEllipsoid model = new OneAxisEllipsoid(100.0, 0.9, frame);
391 Vector3D point = new Vector3D(0.0, 93.7139699, 3.5930796);
392 Vector3D direction = new Vector3D(0.0, 9.0, -2.0);
393 Line line = new Line(point, point.add(direction), 1.0e-10);
394 Assertions.assertNull(model.getIntersectionPoint(line, point, frame, date));
395 }
396
397 @Test
398 void testNegativeZ() {
399 AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
400 Frame frame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
401 OneAxisEllipsoid model = new OneAxisEllipsoid(90.0, 5.0 / 9.0, frame);
402 Vector3D point = new Vector3D(140.0, 0.0, -30.0);
403 GeodeticPoint gp = model.transform(point, frame, date);
404 Vector3D rebuilt = model.transform(gp);
405 Assertions.assertEquals(0.0, rebuilt.distance(point), 1.0e-10);
406 }
407
408 @Test
409 void testNumerousIteration() {
410
411
412
413
414 AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
415 Frame frame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
416 OneAxisEllipsoid model = new OneAxisEllipsoid(100.0, 999.0 / 1000.0, frame);
417 Vector3D point = new Vector3D(100.001, 0.0, 1.0);
418 GeodeticPoint gp = model.transform(point, frame, date);
419 Vector3D rebuilt = model.transform(gp);
420 Assertions.assertEquals(0.0, rebuilt.distance(point), 1.2e-11);
421 }
422
423 @Test
424 void testEquatorialInside() {
425 AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
426 Frame frame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
427 OneAxisEllipsoid model = new OneAxisEllipsoid(90.0, 5.0 / 9.0, frame);
428 for (double rho = 0; rho < model.getEquatorialRadius(); rho += 0.01) {
429 Vector3D point = new Vector3D(rho, 0.0, 0.0);
430 GeodeticPoint gp = model.transform(point, frame, date);
431 Vector3D rebuilt = model.transform(gp);
432 Assertions.assertEquals(0.0, rebuilt.distance(point), 1.0e-10);
433 }
434 }
435
436 @Test
437 void testFarPoint() {
438 AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
439 Frame frame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
440 OneAxisEllipsoid model = new OneAxisEllipsoid(90.0, 5.0 / 9.0, frame);
441 Vector3D point = new Vector3D(1.0e15, 2.0e15, -1.0e12);
442 GeodeticPoint gp = model.transform(point, frame, date);
443 Vector3D rebuilt = model.transform(gp);
444 Assertions.assertEquals(0.0, rebuilt.distance(point), 1.0e-15 * point.getNorm());
445 }
446
447 @Test
448 void testIssue141() {
449 AbsoluteDate date = new AbsoluteDate("2002-03-06T20:50:20.44188731559965033", TimeScalesFactory.getUTC());
450 Frame frame = FramesFactory.getGTOD(IERSConventions.IERS_1996, true);
451 OneAxisEllipsoid model = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
452 Constants.WGS84_EARTH_FLATTENING,
453 frame);
454 Vector3D point = new Vector3D(-6838696.282102453, -2148321.403361013, -0.011907944179711194);
455 GeodeticPoint gp = model.transform(point, frame, date);
456 Vector3D rebuilt = model.transform(gp);
457 Assertions.assertEquals(0.0, rebuilt.distance(point), 1.0e-15 * point.getNorm());
458 }
459
460 @Test
461 void testIntersectionFromPoints() {
462 AbsoluteDate date = new AbsoluteDate(new DateComponents(2008, 03, 21),
463 TimeComponents.H12,
464 TimeScalesFactory.getUTC());
465
466 Frame frame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
467 OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.460, 1 / 298.257222101, frame);
468
469
470
471 final double mu = 3.9860047e14;
472 CircularOrbit circ =
473 new CircularOrbit(7178000.0, 0.5e-4, 0., FastMath.toRadians(90.), FastMath.toRadians(60.),
474 FastMath.toRadians(90.), PositionAngleType.MEAN,
475 FramesFactory.getEME2000(), date, mu);
476
477
478 PVCoordinates pvSatEME2000 = circ.getPVCoordinates();
479 Vector3D pSatItrf = frame.getStaticTransformTo(FramesFactory.getEME2000(), date).transformPosition(pvSatEME2000.getPosition());
480
481
482 GeodeticPoint geoPoint = new GeodeticPoint(FastMath.toRadians(70.), FastMath.toRadians(60.), 0.);
483 Vector3D pointItrf = earth.transform(geoPoint);
484 Line line = new Line(pSatItrf, pointItrf, 1.0e-10);
485 GeodeticPoint geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
486 Assertions.assertEquals(geoPoint.getLongitude(), geoInter.getLongitude(), Utils.epsilonAngle);
487 Assertions.assertEquals(geoPoint.getLatitude(), geoInter.getLatitude(), Utils.epsilonAngle);
488
489
490 geoPoint = new GeodeticPoint(FastMath.toRadians(65.), FastMath.toRadians(-120.), 0.);
491 pointItrf = earth.transform(geoPoint);
492 line = new Line(pSatItrf, pointItrf, 1.0e-10);
493 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
494 Assertions.assertEquals(geoPoint.getLongitude(), geoInter.getLongitude(), Utils.epsilonAngle);
495 Assertions.assertEquals(geoPoint.getLatitude(), geoInter.getLatitude(), Utils.epsilonAngle);
496
497
498 geoPoint = new GeodeticPoint(FastMath.toRadians(30.), FastMath.toRadians(60.), 0.);
499 pointItrf = earth.transform(geoPoint);
500 line = new Line(pSatItrf, pointItrf, 1.0e-10);
501
502 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
503
504
505 Assertions.assertEquals(1.04437199, geoInter.getLongitude(), Utils.epsilonAngle);
506 Assertions.assertEquals(1.36198012, geoInter.getLatitude(), Utils.epsilonAngle);
507
508
509
510 circ =
511 new CircularOrbit(7178000.0, 0.5e-4, 0., FastMath.toRadians(1.e-4), FastMath.toRadians(0.),
512 FastMath.toRadians(0.), PositionAngleType.MEAN,
513 FramesFactory.getEME2000(), date, mu);
514
515
516 pvSatEME2000 = circ.getPVCoordinates();
517 pSatItrf = frame.getStaticTransformTo(FramesFactory.getEME2000(), date).transformPosition(pvSatEME2000.getPosition());
518
519
520 geoPoint = new GeodeticPoint(FastMath.toRadians(5.), FastMath.toRadians(0.), 0.);
521 pointItrf = earth.transform(geoPoint);
522 line = new Line(pSatItrf, pointItrf, 1.0e-10);
523 Assertions.assertTrue(line.toSubSpace(pSatItrf).getX() < 0);
524 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
525 Assertions.assertEquals(geoPoint.getLongitude(), geoInter.getLongitude(), Utils.epsilonAngle);
526 Assertions.assertEquals(geoPoint.getLatitude(), geoInter.getLatitude(), Utils.epsilonAngle);
527
528
529 GeodeticPoint geoInter2 = earth.getIntersectionPoint(line, line.toSpace(new Vector1D(-line.toSubSpace(pSatItrf).getX())), frame, date);
530 Assertions.assertTrue(FastMath.abs(geoInter.getLongitude() - geoInter2.getLongitude()) > FastMath.toRadians(0.1));
531 Assertions.assertTrue(FastMath.abs(geoInter.getLatitude() - geoInter2.getLatitude()) > FastMath.toRadians(0.1));
532
533
534 geoPoint = new GeodeticPoint(FastMath.toRadians(-5.), FastMath.toRadians(0.), 0.);
535 pointItrf = earth.transform(geoPoint);
536 line = new Line(pSatItrf, pointItrf, 1.0e-10);
537 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
538 Assertions.assertEquals(geoPoint.getLongitude(), geoInter.getLongitude(), Utils.epsilonAngle);
539 Assertions.assertEquals(geoPoint.getLatitude(), geoInter.getLatitude(), Utils.epsilonAngle);
540
541
542 geoPoint = new GeodeticPoint(FastMath.toRadians(40.), FastMath.toRadians(0.), 0.);
543 pointItrf = earth.transform(geoPoint);
544 line = new Line(pSatItrf, pointItrf, 1.0e-10);
545 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
546 Assertions.assertEquals(-0.00768481, geoInter.getLongitude(), Utils.epsilonAngle);
547 Assertions.assertEquals( 0.32180410, geoInter.getLatitude(), Utils.epsilonAngle);
548
549
550
551
552 circ =
553 new CircularOrbit(7178000.0, 0.5e-4, 0., FastMath.toRadians(50.), FastMath.toRadians(0.),
554 FastMath.toRadians(90.), PositionAngleType.MEAN,
555 FramesFactory.getEME2000(), date, mu);
556
557
558 pvSatEME2000 = circ.getPVCoordinates();
559 pSatItrf = frame.getStaticTransformTo(FramesFactory.getEME2000(), date).transformPosition(pvSatEME2000.getPosition());
560
561
562 geoPoint = new GeodeticPoint(FastMath.toRadians(40.), FastMath.toRadians(90.), 0.);
563 pointItrf = earth.transform(geoPoint);
564 line = new Line(pSatItrf, pointItrf, 1.0e-10);
565 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
566 Assertions.assertEquals(geoPoint.getLongitude(), geoInter.getLongitude(), Utils.epsilonAngle);
567 Assertions.assertEquals(geoPoint.getLatitude(), geoInter.getLatitude(), Utils.epsilonAngle);
568
569
570 geoPoint = new GeodeticPoint(FastMath.toRadians(60.), FastMath.toRadians(90.), 0.);
571 pointItrf = earth.transform(geoPoint);
572 line = new Line(pSatItrf, pointItrf, 1.0e-10);
573 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
574 Assertions.assertEquals(geoPoint.getLongitude(), geoInter.getLongitude(), Utils.epsilonAngle);
575 Assertions.assertEquals(geoPoint.getLatitude(), geoInter.getLatitude(), Utils.epsilonAngle);
576
577
578 geoPoint = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(90.), 0.);
579 pointItrf = earth.transform(geoPoint);
580 line = new Line(pSatItrf, pointItrf, 1.0e-10);
581 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
582 Assertions.assertEquals(FastMath.toRadians(89.5364061088196), geoInter.getLongitude(), Utils.epsilonAngle);
583 Assertions.assertEquals(FastMath.toRadians(35.555543683351125), geoInter.getLatitude(), Utils.epsilonAngle);
584
585 }
586
587 @Test
588 void testMovingGeodeticPointSymmetry() {
589
590 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
591 Constants.WGS84_EARTH_FLATTENING,
592 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
593 double lat0 = FastMath.toRadians(60.0);
594 double lon0 = FastMath.toRadians(25.0);
595 double alt0 = 100.0;
596 double lat1 = 1.0e-3;
597 double lon1 = -2.0e-3;
598 double alt1 = 1.2;
599 double lat2 = -1.0e-5;
600 double lon2 = -3.0e-5;
601 double alt2 = -0.01;
602 final DSFactory factory = new DSFactory(1, 2);
603 final DerivativeStructure latDS = factory.build(lat0, lat1, lat2);
604 final DerivativeStructure lonDS = factory.build(lon0, lon1, lon2);
605 final DerivativeStructure altDS = factory.build(alt0, alt1, alt2);
606
607
608 PVCoordinates pv = new PVCoordinates(earth.transform(new FieldGeodeticPoint<>(latDS, lonDS, altDS)));
609 FieldGeodeticPoint<UnivariateDerivative2> rebuilt = earth.transform(pv, earth.getBodyFrame(), null);
610 Assertions.assertEquals(lat0, rebuilt.getLatitude().getReal(), 1.0e-16);
611 Assertions.assertEquals(lat1, rebuilt.getLatitude().getPartialDerivative(1), 5.0e-19);
612 Assertions.assertEquals(lat2, rebuilt.getLatitude().getPartialDerivative(2), 5.0e-14);
613 Assertions.assertEquals(lon0, rebuilt.getLongitude().getReal(), 1.0e-16);
614 Assertions.assertEquals(lon1, rebuilt.getLongitude().getPartialDerivative(1), 5.0e-19);
615 Assertions.assertEquals(lon2, rebuilt.getLongitude().getPartialDerivative(2), 1.0e-20);
616 Assertions.assertEquals(alt0, rebuilt.getAltitude().getReal(), 2.0e-11);
617 Assertions.assertEquals(alt1, rebuilt.getAltitude().getPartialDerivative(1), 6.0e-13);
618 Assertions.assertEquals(alt2, rebuilt.getAltitude().getPartialDerivative(2), 2.0e-14);
619
620 }
621
622 @Test
623 void testMovingGeodeticPoint() {
624
625 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
626 Constants.WGS84_EARTH_FLATTENING,
627 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
628 double lat0 = FastMath.toRadians(60.0);
629 double lon0 = FastMath.toRadians(25.0);
630 double alt0 = 100.0;
631 double lat1 = 1.0e-3;
632 double lon1 = -2.0e-3;
633 double alt1 = 1.2;
634 double lat2 = -1.0e-5;
635 double lon2 = -3.0e-5;
636 double alt2 = -0.01;
637 final DSFactory factory = new DSFactory(1, 2);
638 final DerivativeStructure latDS = factory.build(lat0, lat1, lat2);
639 final DerivativeStructure lonDS = factory.build(lon0, lon1, lon2);
640 final DerivativeStructure altDS = factory.build(alt0, alt1, alt2);
641
642
643 PVCoordinates pv = new PVCoordinates(earth.transform(new FieldGeodeticPoint<>(latDS, lonDS, altDS)));
644
645
646 FiniteDifferencesDifferentiator differentiator =
647 new FiniteDifferencesDifferentiator(5, 0.1);
648 UnivariateDifferentiableFunction fx =
649 differentiator.differentiate(new UnivariateFunction() {
650 public double value(double dt) {
651 GeodeticPoint gp =
652 new GeodeticPoint(latDS.taylor(dt), lonDS.taylor(dt), altDS.taylor(dt));
653 return earth.transform(gp).getX();
654 }
655 });
656 UnivariateDifferentiableFunction fy =
657 differentiator.differentiate(new UnivariateFunction() {
658 public double value(double dt) {
659 GeodeticPoint gp =
660 new GeodeticPoint(latDS.taylor(dt), lonDS.taylor(dt), altDS.taylor(dt));
661 return earth.transform(gp).getY();
662 }
663 });
664 UnivariateDifferentiableFunction fz =
665 differentiator.differentiate(new UnivariateFunction() {
666 public double value(double dt) {
667 GeodeticPoint gp =
668 new GeodeticPoint(latDS.taylor(dt), lonDS.taylor(dt), altDS.taylor(dt));
669 return earth.transform(gp).getZ();
670 }
671 });
672 DerivativeStructure dtZero = factory.variable(0, 0.0);
673 DerivativeStructure xDS = fx.value(dtZero);
674 DerivativeStructure yDS = fy.value(dtZero);
675 DerivativeStructure zDS = fz.value(dtZero);
676 Assertions.assertEquals(xDS.getValue(), pv.getPosition().getX(),
677 2.0e-20 * FastMath.abs(xDS.getValue()));
678 Assertions.assertEquals(xDS.getPartialDerivative(1), pv.getVelocity().getX(),
679 2.0e-12 * FastMath.abs(xDS.getPartialDerivative(1)));
680 Assertions.assertEquals(xDS.getPartialDerivative(2), pv.getAcceleration().getX(),
681 2.0e-9 * FastMath.abs(xDS.getPartialDerivative(2)));
682 Assertions.assertEquals(yDS.getValue(), pv.getPosition().getY(),
683 2.0e-20 * FastMath.abs(yDS.getValue()));
684 Assertions.assertEquals(yDS.getPartialDerivative(1), pv.getVelocity().getY(),
685 2.0e-12 * FastMath.abs(yDS.getPartialDerivative(1)));
686 Assertions.assertEquals(yDS.getPartialDerivative(2), pv.getAcceleration().getY(),
687 2.0e-9 * FastMath.abs(yDS.getPartialDerivative(2)));
688 Assertions.assertEquals(zDS.getValue(), pv.getPosition().getZ(),
689 2.0e-20 * FastMath.abs(zDS.getValue()));
690 Assertions.assertEquals(zDS.getPartialDerivative(1), pv.getVelocity().getZ(),
691 2.0e-12 * FastMath.abs(zDS.getPartialDerivative(1)));
692 Assertions.assertEquals(zDS.getPartialDerivative(2), pv.getAcceleration().getZ(),
693 2.0e-9 * FastMath.abs(zDS.getPartialDerivative(2)));
694 }
695
696 private void checkCartesianToEllipsoidic(double ae, double f,
697 double x, double y, double z,
698 double longitude, double latitude,
699 double altitude) {
700
701 AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
702 Frame frame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
703 OneAxisEllipsoid model = new OneAxisEllipsoid(ae, f, frame);
704 GeodeticPoint gp = model.transform(new Vector3D(x, y, z), frame, date);
705 Assertions.assertEquals(longitude, MathUtils.normalizeAngle(gp.getLongitude(), longitude), 1.0e-10);
706 Assertions.assertEquals(latitude, gp.getLatitude(), 1.0e-10);
707 Assertions.assertEquals(altitude, gp.getAltitude(), 1.0e-10 * FastMath.abs(ae));
708 Vector3D rebuiltNadir = Vector3D.crossProduct(gp.getSouth(), gp.getWest());
709 Assertions.assertEquals(0, rebuiltNadir.subtract(gp.getNadir()).getNorm(), 1.0e-15);
710
711 FieldGeodeticPoint<Binary64> gp64 = model.transform(new FieldVector3D<Binary64>(new Binary64(x),
712 new Binary64(y),
713 new Binary64(z)),
714 frame,
715 new FieldAbsoluteDate<>(Binary64Field.getInstance(), date));
716 Assertions.assertEquals(longitude, MathUtils.normalizeAngle(gp64.getLongitude().getReal(), longitude), 1.0e-10);
717 Assertions.assertEquals(latitude, gp64.getLatitude().getReal(), 1.0e-10);
718 Assertions.assertEquals(altitude, gp64.getAltitude().getReal(), 1.0e-10 * FastMath.abs(ae));
719 FieldVector3D<Binary64> rebuiltNadir64 = FieldVector3D.crossProduct(gp64.getSouth(), gp64.getWest());
720 Assertions.assertEquals(0, rebuiltNadir64.subtract(gp64.getNadir()).getNorm().getReal(), 1.0e-15);
721
722
723 gp = model.transform(model.projectToGround(new Vector3D(x, y, z), date, frame), frame, date);
724 Assertions.assertEquals(longitude, MathUtils.normalizeAngle(gp.getLongitude(), longitude), 1.0e-10);
725 Assertions.assertEquals(latitude, gp.getLatitude(), 1.0e-10);
726 Assertions.assertEquals(0.0, gp.getAltitude(), 1.0e-10 * FastMath.abs(ae));
727
728
729 FieldGeodeticPoint<UnivariateDerivative2> gpDs = model.transform(
730 model.projectToGround(
731 new TimeStampedPVCoordinates(
732 date,
733 new Vector3D(x, y, z),
734 new Vector3D(1, 2, 3)),
735 frame),
736 frame,
737 date);
738 Assertions.assertEquals(longitude, MathUtils.normalizeAngle(gpDs.getLongitude().getReal(), longitude), 1.0e-10);
739 Assertions.assertEquals(latitude, gpDs.getLatitude().getReal(), 1.0e-10);
740 Assertions.assertEquals(0.0, gpDs.getAltitude().getReal(), 1.0e-10 * FastMath.abs(ae));
741
742 }
743
744 @Test
745 void testTransformVsOldIterativeSobol() {
746
747 OneAxisEllipsoid model = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
748 Constants.WGS84_EARTH_FLATTENING,
749 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
750 SobolSequenceGenerator sobol = new SobolSequenceGenerator(3);
751 final double rMax = 10 * model.getEquatorialRadius();
752 Stream<Vector3D> points = Stream.generate(() -> {
753 final double[] v = sobol.nextVector();
754 return new Vector3D(rMax * (2 * v[0] - 1), rMax * (2 * v[1] - 1), rMax * (2 * v[2] - 1));
755 }).limit(1000000);
756
757 doTestTransformVsOldIterative(model, points, 2.0e-15, 1.0e-15, 8.0e-14 * rMax);
758 }
759
760 @Test
761 void testTransformVsOldIterativePolarAxis() {
762 OneAxisEllipsoid model = new OneAxisEllipsoid(90, 5.0 / 9.0,
763 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
764 Stream<Vector3D> points = DoubleStream.iterate(0, x -> x + 1.0).limit(150).mapToObj(z -> new Vector3D(0, 0, z));
765 doTestTransformVsOldIterative(model, points, 2.0e-15, 1.0e-15, 1.0e-14 * model.getEquatorialRadius());
766 }
767
768 @Test
769 void testTransformVsOldIterativeEquatorial() {
770 OneAxisEllipsoid model = new OneAxisEllipsoid(90, 5.0 / 9.0,
771 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
772 Stream<Vector3D> points = DoubleStream.iterate(0, x -> x + 1.0).limit(150).mapToObj(x -> new Vector3D(x, 0, 0));
773 doTestTransformVsOldIterative(model, points, 2.0e-15, 1.0e-15, 1.0e-14 * model.getEquatorialRadius());
774 }
775
776 @Test
777 void testIssue373() {
778 final Frame ecef = FramesFactory.getITRF(IERSConventions.IERS_2010,true);
779 final OneAxisEllipsoid earth = new OneAxisEllipsoid(6378137, 1./298.257223563, ecef);
780 final Vector3D sunPos = new Vector3D(-149757851422.23358, 8410610314.781021, 14717269835.161688 );
781 final GeodeticPoint sunGP = earth.transform(sunPos, ecef, null);
782 Assertions.assertEquals(5.603878, FastMath.toDegrees(sunGP.getLatitude()), 1.0e-6);
783 }
784
785 @Test
786 void testIsometricLatitude() {
787 final Frame ecef = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
788 final OneAxisEllipsoid earth = new OneAxisEllipsoid(6378137, 1. / 298.257223563, ecef);
789
790 Assertions.assertEquals(0., earth.geodeticToIsometricLatitude(0), 1.0e-9);
791 Assertions.assertEquals(0., earth.geodeticToIsometricLatitude(2.0e-13), 1.0e-9);
792
793 for (final double lat: new double[] {
794 FastMath.toRadians(10),
795 FastMath.toRadians(-45),
796 FastMath.toRadians(80),
797 FastMath.toRadians(-90)}) {
798 final double eSinLat = earth.getEccentricity() * FastMath.sin(lat);
799 final double term1 = FastMath.log(FastMath.tan(FastMath.PI / 4. + lat / 2.));
800 final double term2 = (earth.getEccentricity() / 2.) * FastMath.log((1. - eSinLat) / (1 + eSinLat));
801
802 Assertions.assertEquals(term1 + term2, earth.geodeticToIsometricLatitude(lat), 1.0e-12);
803 }
804 }
805
806 @Test
807 void testFieldIsometricLatitude() {
808 final Frame ecef = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
809 final OneAxisEllipsoid earth = new OneAxisEllipsoid(6378137, 1. / 298.257223563, ecef);
810
811 final Binary64Field field = Binary64Field.getInstance();
812 Assertions.assertEquals(Binary64Field.getInstance().getZero(), earth.geodeticToIsometricLatitude(field.getOne().newInstance(0.)));
813 Assertions.assertEquals(Binary64Field.getInstance().getZero(), earth.geodeticToIsometricLatitude(field.getOne().newInstance(2.0e-13)));
814
815 final Binary64 ecc = field.getZero().newInstance(earth.getEccentricity());
816 for (final Binary64 lat: new Binary64[] {
817 field.getOne().newInstance(FastMath.toRadians(10)),
818 field.getOne().newInstance(FastMath.toRadians(-45)),
819 field.getOne().newInstance(FastMath.toRadians(80)),
820 field.getOne().newInstance(FastMath.toRadians(-90))}) {
821 final Binary64 eSinLat = ecc.multiply(FastMath.sin(lat));
822 final Binary64 term1 = FastMath.log(FastMath.tan(lat.getPi().divide(4.).add(lat.divide(2.))));
823 final Binary64 term2 = ecc.divide(2.).multiply(FastMath.log(field.getOne().subtract(eSinLat).divide(field.getOne().add(eSinLat))));
824
825 Assertions.assertEquals(term1.add(term2), earth.geodeticToIsometricLatitude(lat));
826 }
827 }
828
829 @Test
830 void testAzimuthBetweenPoints() {
831 final Frame ecef = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
832 final OneAxisEllipsoid earth = new OneAxisEllipsoid(6378137, 1. / 298.257223563, ecef);
833
834
835 final GeodeticPoint newYork = new GeodeticPoint(FastMath.toRadians(40.71427), FastMath.toRadians(-74.00597), 0);
836 final GeodeticPoint chicago = new GeodeticPoint(FastMath.toRadians(41.85003), FastMath.toRadians(-87.65005), 0);
837 final GeodeticPoint london = new GeodeticPoint(FastMath.toRadians(51.5), FastMath.toRadians(-0.16667), 0);
838 final GeodeticPoint berlin = new GeodeticPoint(FastMath.toRadians(52.523403), FastMath.toRadians(13.4114), 0);
839 final GeodeticPoint perth = new GeodeticPoint(FastMath.toRadians(-31.952712), FastMath.toRadians(115.8604796), 0);
840
841
842 Assertions.assertEquals(276.297499, FastMath.toDegrees(earth.azimuthBetweenPoints(newYork, chicago)), 1.0e-6);
843 Assertions.assertEquals(78.0854898, FastMath.toDegrees(earth.azimuthBetweenPoints(newYork, london)), 1.0e-6);
844 Assertions.assertEquals(83.0357553, FastMath.toDegrees(earth.azimuthBetweenPoints(london, berlin)), 1.0e-6);
845 Assertions.assertEquals(132.894864, FastMath.toDegrees(earth.azimuthBetweenPoints(berlin, perth)), 1.0e-6);
846 Assertions.assertEquals(65.3853195, FastMath.toDegrees(earth.azimuthBetweenPoints(perth, newYork)), 1.0e-6);
847 }
848
849 @Test
850 void testAzimuthBetweenFieldPoints() {
851 final Frame ecef = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
852 final OneAxisEllipsoid earth = new OneAxisEllipsoid(6378137, 1. / 298.257223563, ecef);
853
854 final Binary64Field field = Binary64Field.getInstance();
855
856
857 final FieldGeodeticPoint<Binary64> newYork = new FieldGeodeticPoint<>(
858 FastMath.toRadians(field.getZero().add(40.71427)),
859 FastMath.toRadians(field.getZero().add(-74.00597)), field.getZero());
860 final FieldGeodeticPoint<Binary64> chicago = new FieldGeodeticPoint<>(
861 FastMath.toRadians(field.getZero().add(41.85003)),
862 FastMath.toRadians(field.getZero().add(-87.65005)), field.getZero());
863
864 final FieldGeodeticPoint<Binary64> london = new FieldGeodeticPoint<>(
865 FastMath.toRadians(field.getZero().add(51.5)),
866 FastMath.toRadians(field.getZero().add(-0.16667)), field.getZero());
867 final FieldGeodeticPoint<Binary64> berlin = new FieldGeodeticPoint<>(
868 FastMath.toRadians(field.getZero().add(52.523403)),
869 FastMath.toRadians(field.getZero().add(13.4114)), field.getZero());
870 final FieldGeodeticPoint<Binary64> perth = new FieldGeodeticPoint<>(
871 FastMath.toRadians(field.getZero().add(-31.952712)),
872 FastMath.toRadians(field.getZero().add(115.8604796)), field.getZero());
873
874
875 Assertions.assertEquals(276.297499, FastMath.toDegrees(earth.azimuthBetweenPoints(newYork, chicago)).getReal(), 1.0e-6);
876 Assertions.assertEquals(78.0854898, FastMath.toDegrees(earth.azimuthBetweenPoints(newYork, london)).getReal(), 1.0e-6);
877 Assertions.assertEquals(83.0357553, FastMath.toDegrees(earth.azimuthBetweenPoints(london, berlin)).getReal(), 1.0e-6);
878 Assertions.assertEquals(132.894864, FastMath.toDegrees(earth.azimuthBetweenPoints(berlin, perth)).getReal(), 1.0e-6);
879 Assertions.assertEquals(65.3853195, FastMath.toDegrees(earth.azimuthBetweenPoints(perth, newYork)).getReal(), 1.0e-6);
880 }
881
882 @Test
883 void testPointNearCenter1() {
884 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
885 Constants.WGS84_EARTH_FLATTENING,
886 FramesFactory.getITRF(IERSConventions.IERS_2010, false));
887 final Vector3D p1 = new Vector3D( 14605530.402633, 7681001.886684, 24582223.005261);
888 final Vector3D p2 = new Vector3D(-14650836.411867, -7561887.405778, -24575352.170908);
889 final Vector3D pMid = new Vector3D(0.5, p1, 0.5, p2);
890 final GeodeticPoint gp = earth.transform(pMid, earth.getFrame(), null);
891 Vector3D rebuilt = earth.transform(gp);
892 Assertions.assertEquals(0.0, rebuilt.distance(pMid), 1.5e-9);
893 }
894
895 @Test
896 void testPointNearCenter2() {
897 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
898 Constants.WGS84_EARTH_FLATTENING,
899 FramesFactory.getITRF(IERSConventions.IERS_2010, false));
900 final Vector3D pMid = new Vector3D(-20923.23737959098, 56464.586571323685, -7647.317096056417);
901 final GeodeticPoint gp = earth.transform(pMid, earth.getFrame(), null);
902 Vector3D rebuilt = earth.transform(gp);
903
904 Assertions.assertEquals(540.598, rebuilt.distance(pMid), 1.0e-3);
905 }
906
907
908 @Test
909 void testLowestAltitudeIntermediate() {
910 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
911 Constants.WGS84_EARTH_FLATTENING,
912 FramesFactory.getITRF(IERSConventions.IERS_2010, false));
913 final Vector3D close = earth.transform(new GeodeticPoint(FastMath.toRadians(12.0),
914 FastMath.toRadians(34.0),
915 100000.0));
916 final Vector3D far = earth.transform(new GeodeticPoint(FastMath.toRadians(-47.0),
917 FastMath.toRadians( 43.0),
918 10000000.0));
919 final GeodeticPoint lowestA = earth.lowestAltitudeIntermediate(close, far);
920 Assertions.assertEquals(0.0, Vector3D.distance(close, earth.transform(lowestA)), 2.0e-9);
921 final GeodeticPoint lowestB = earth.lowestAltitudeIntermediate(far, close);
922 Assertions.assertEquals(0.0, Vector3D.distance(close, earth.transform(lowestB)), 2.0e-9);
923
924 final double h = 5000000.0;
925 final Vector3D p1 = earth.transform(new GeodeticPoint(FastMath.toRadians(-20.0),
926 FastMath.toRadians(-12.0),
927 h));
928 final Vector3D p2 = earth.transform(new GeodeticPoint(FastMath.toRadians(17.0),
929 FastMath.toRadians(13.0),
930 h));
931 final GeodeticPoint lowest = earth.lowestAltitudeIntermediate(p1, p2);
932 Assertions.assertTrue(lowest.getAltitude() < h);
933 final GeodeticPoint oneCentimeterBefore = earth.transform(new Vector3D( 1.0, earth.transform(lowest),
934 -1.0e-2, p2.subtract(p1).normalize()),
935 earth.getBodyFrame(), null);
936 Assertions.assertTrue(oneCentimeterBefore.getAltitude() > lowest.getAltitude());
937 final GeodeticPoint oneCentimeterAfter = earth.transform(new Vector3D( 1.0, earth.transform(lowest),
938 +1.0e-2, p2.subtract(p1).normalize()),
939 earth.getBodyFrame(), null);
940 Assertions.assertTrue(oneCentimeterAfter.getAltitude() > lowest.getAltitude());
941
942 }
943
944 @Test
945 void testLowestAltitudeIntermediateField() {
946 doTestLowestAltitudeIntermediateField(Binary64Field.getInstance());
947 }
948
949 private <T extends CalculusFieldElement<T>> void doTestLowestAltitudeIntermediateField(final Field<T> field) {
950
951 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
952 Constants.WGS84_EARTH_FLATTENING,
953 FramesFactory.getITRF(IERSConventions.IERS_2010, false));
954 final FieldVector3D<T> close = earth.transform(new FieldGeodeticPoint<>(FastMath.toRadians(field.getZero().newInstance(12.0)),
955 FastMath.toRadians(field.getZero().newInstance(34.0)),
956 field.getZero().newInstance(100000.0)));
957 final FieldVector3D<T> far = earth.transform(new FieldGeodeticPoint<>(FastMath.toRadians(field.getZero().newInstance(-47.0)),
958 FastMath.toRadians(field.getZero().newInstance( 43.0)),
959 field.getZero().newInstance(10000000.0)));
960 final FieldGeodeticPoint<T> lowestA = earth.lowestAltitudeIntermediate(close, far);
961 Assertions.assertEquals(0.0, FieldVector3D.distance(close, earth.transform(lowestA)).getReal(), 2.0e-9);
962 final FieldGeodeticPoint<T> lowestB = earth.lowestAltitudeIntermediate(far, close);
963 Assertions.assertEquals(0.0, FieldVector3D.distance(close, earth.transform(lowestB)).getReal(), 2.0e-9);
964
965 final double h =5000000.0;
966 final FieldVector3D<T> p1 = earth.transform(new FieldGeodeticPoint<>(FastMath.toRadians(field.getZero().newInstance(-20.0)),
967 FastMath.toRadians(field.getZero().newInstance(-12.0)),
968 field.getZero().newInstance(h)));
969 final FieldVector3D<T> p2 = earth.transform(new FieldGeodeticPoint<>(FastMath.toRadians(field.getZero().newInstance(17.0)),
970 FastMath.toRadians(field.getZero().newInstance(13.0)),
971 field.getZero().newInstance(h)));
972 final FieldGeodeticPoint<T> lowest = earth.lowestAltitudeIntermediate(p1, p2);
973 Assertions.assertTrue(lowest.getAltitude().getReal() < h);
974 final FieldGeodeticPoint<T> oneCentimeterBefore = earth.transform(new FieldVector3D<>( 1.0, earth.transform(lowest),
975 -1.0e-2, p2.subtract(p1).normalize()),
976 earth.getBodyFrame(), null);
977 Assertions.assertTrue(oneCentimeterBefore.getAltitude().subtract(lowest.getAltitude()).getReal() > 0);
978 final FieldGeodeticPoint<T> oneCentimeterAfter = earth.transform(new FieldVector3D<>( 1.0, earth.transform(lowest),
979 +1.0e-2, p2.subtract(p1).normalize()),
980 earth.getBodyFrame(), null);
981 Assertions.assertTrue(oneCentimeterAfter.getAltitude().subtract(lowest.getAltitude()).getReal() > 0);
982
983 }
984
985 private void doTestTransformVsOldIterative(OneAxisEllipsoid model,
986 Stream<Vector3D> points,
987 double latitudeTolerance, double longitudeTolerance,
988 double altitudeTolerance) {
989 points.forEach(point -> {
990 try {
991 GeodeticPoint reference = transformOldIterative(model, point, model.getBodyFrame(), null);
992 GeodeticPoint result = model.transform(point, model.getBodyFrame(), null);
993 Assertions.assertEquals(reference.getLatitude(), result.getLatitude(), latitudeTolerance);
994 Assertions.assertEquals(reference.getLongitude(), result.getLongitude(), longitudeTolerance);
995 Assertions.assertEquals(reference.getAltitude(), result.getAltitude(), altitudeTolerance);
996 } catch (OrekitException oe) {
997 Assertions.fail(oe.getLocalizedMessage());
998 }
999 });
1000
1001 }
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015 private GeodeticPoint transformOldIterative(final OneAxisEllipsoid model,
1016 final Vector3D point,
1017 final Frame frame,
1018 final AbsoluteDate date) {
1019
1020
1021 final Vector3D pointInBodyFrame = frame.getTransformTo(model.getBodyFrame(), date).transformPosition(point);
1022 final double r2 = pointInBodyFrame.getX() * pointInBodyFrame.getX() +
1023 pointInBodyFrame.getY() * pointInBodyFrame.getY();
1024 final double r = FastMath.sqrt(r2);
1025 final double z = pointInBodyFrame.getZ();
1026
1027
1028 final Ellipse meridian = new Ellipse(Vector3D.ZERO,
1029 new Vector3D(pointInBodyFrame.getX() / r, pointInBodyFrame.getY() / r, 0),
1030 Vector3D.PLUS_K,
1031 model.getA(), model.getC(), model.getBodyFrame());
1032
1033
1034 final Vector2D ellipsePoint = meridian.projectToEllipse(new Vector2D(r, z));
1035
1036
1037 final double dr = r - ellipsePoint.getX();
1038 final double dz = z - ellipsePoint.getY();
1039 final double ae2 = model.getA() * model.getA();
1040 final double f = model.getFlattening();
1041 final double g = 1.0 - f;
1042 final double g2 = g * g;
1043 final double insideIfNegative = g2 * (r2 - ae2) + z * z;
1044
1045 return new GeodeticPoint(FastMath.atan2(ellipsePoint.getY(), g2 * ellipsePoint.getX()),
1046 FastMath.atan2(pointInBodyFrame.getY(), pointInBodyFrame.getX()),
1047 FastMath.copySign(FastMath.hypot(dr, dz), insideIfNegative));
1048
1049 }
1050
1051 @BeforeEach
1052 public void setUp() {
1053 Utils.setDataRoot("regular-data");
1054 }
1055
1056 }
1057