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