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 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 Assert.assertEquals(0, errors[0], 1.0e-16);
237 Assert.assertEquals(0, errors[1], 2.0e-12);
238 Assert.assertEquals(0, errors[2], 2.0e-4);
239
240 }
241
242 @Test
243 public 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 Assert.assertEquals(0.0, Vector3D.distance(ground1.getPosition(), ground2.getPosition()), 1.0e-12);
260 Assert.assertEquals(0.0, Vector3D.distance(ground1.getVelocity(), ground2.getVelocity()), 2.0e-12);
261 Assert.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 TimeStampedPVCoordinates computed =
277 model.projectToGround(TimeStampedPVCoordinates.interpolate(date,
278 CartesianDerivativesFilter.USE_P,
279 pvList),
280 frame);
281 TimeStampedPVCoordinates reference =
282 TimeStampedPVCoordinates.interpolate(date,
283 CartesianDerivativesFilter.USE_P,
284 groundPVList);
285
286 TimeStampedPVCoordinates pv0 = provider.getPVCoordinates(date, frame);
287 Vector3D p0 = pv0.getPosition();
288 Vector3D v0 = pv0.getVelocity();
289 Vector3D a0 = pv0.getAcceleration();
290
291 return new double[] {
292 Vector3D.distance(computed.getPosition(), reference.getPosition()) / p0.getNorm(),
293 Vector3D.distance(computed.getVelocity(), reference.getVelocity()) / v0.getNorm(),
294 Vector3D.distance(computed.getAcceleration(), reference.getAcceleration()) / a0.getNorm(),
295 };
296
297 }
298
299 @Test
300 public void testGroundProjectionTaylor() {
301 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
302 Frame eme2000 = FramesFactory.getEME2000();
303 OneAxisEllipsoid model =
304 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
305 Constants.WGS84_EARTH_FLATTENING,
306 itrf);
307
308 TimeStampedPVCoordinates initPV =
309 new TimeStampedPVCoordinates(AbsoluteDate.J2000_EPOCH.shiftedBy(584.),
310 new Vector3D(3220103., 69623., 6449822.),
311 new Vector3D(6414.7, -2006., -3180.),
312 Vector3D.ZERO);
313 Orbit orbit = new EquinoctialOrbit(initPV, eme2000, Constants.EIGEN5C_EARTH_MU);
314
315 TimeStampedPVCoordinates pv0 = orbit.getPVCoordinates(orbit.getDate(), model.getBodyFrame());
316 PVCoordinatesProvider groundTaylor =
317 model.projectToGround(pv0, model.getBodyFrame()).toTaylorProvider(model.getBodyFrame());
318
319 TimeStampedPVCoordinates g0 = groundTaylor.getPVCoordinates(orbit.getDate(), model.getBodyFrame());
320 Vector3D zenith = pv0.getPosition().subtract(g0.getPosition()).normalize();
321 Vector3D acrossTrack = Vector3D.crossProduct(zenith, g0.getVelocity()).normalize();
322 Vector3D alongTrack = Vector3D.crossProduct(acrossTrack, zenith).normalize();
323 for (double dt = -1; dt < 1; dt += 0.01) {
324 AbsoluteDate date = orbit.getDate().shiftedBy(dt);
325 Vector3D taylorP = groundTaylor.getPVCoordinates(date, model.getBodyFrame()).getPosition();
326 Vector3D refP = model.projectToGround(orbit.getPVCoordinates(date, model.getBodyFrame()).getPosition(),
327 date, model.getBodyFrame());
328 Vector3D delta = taylorP.subtract(refP);
329 Assert.assertEquals(0.0, Vector3D.dotProduct(delta, alongTrack), 0.0015);
330 Assert.assertEquals(0.0, Vector3D.dotProduct(delta, acrossTrack), 0.0007);
331 Assert.assertEquals(0.0, Vector3D.dotProduct(delta, zenith), 0.00002);
332 }
333
334 }
335
336 @Test
337 public void testLineIntersection() {
338 AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
339 Frame frame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
340
341 OneAxisEllipsoid model = new OneAxisEllipsoid(100.0, 0.9, frame);
342 Vector3D point = new Vector3D(0.0, 93.7139699, 3.5930796);
343 Vector3D direction = new Vector3D(0.0, 1.0, 1.0);
344 Line line = new Line(point, point.add(direction), 1.0e-10);
345 GeodeticPoint gp = model.getIntersectionPoint(line, point, frame, date);
346 Assert.assertEquals(gp.getAltitude(), 0.0, 1.0e-12);
347 Assert.assertTrue(line.contains(model.transform(gp)));
348
349 model = new OneAxisEllipsoid(100.0, 0.9, frame);
350 point = new Vector3D(0.0, -93.7139699, -3.5930796);
351 direction = new Vector3D(0.0, -1.0, -1.0);
352 line = new Line(point, point.add(direction), 1.0e-10).revert();
353 gp = model.getIntersectionPoint(line, point, frame, date);
354 Assert.assertTrue(line.contains(model.transform(gp)));
355
356 model = new OneAxisEllipsoid(100.0, 0.9, frame);
357 point = new Vector3D(0.0, -93.7139699, 3.5930796);
358 direction = new Vector3D(0.0, -1.0, 1.0);
359 line = new Line(point, point.add(direction), 1.0e-10);
360 gp = model.getIntersectionPoint(line, point, frame, date);
361 Assert.assertTrue(line.contains(model.transform(gp)));
362
363 model = new OneAxisEllipsoid(100.0, 0.9, frame);
364 point = new Vector3D(-93.7139699, 0.0, 3.5930796);
365 direction = new Vector3D(-1.0, 0.0, 1.0);
366 line = new Line(point, point.add(direction), 1.0e-10);
367 gp = model.getIntersectionPoint(line, point, frame, date);
368 Assert.assertTrue(line.contains(model.transform(gp)));
369 Assert.assertFalse(line.contains(new Vector3D(0, 0, 7000000)));
370
371 point = new Vector3D(0.0, 0.0, 110);
372 direction = new Vector3D(0.0, 0.0, 1.0);
373 line = new Line(point, point.add(direction), 1.0e-10);
374 gp = model.getIntersectionPoint(line, point, frame, date);
375 Assert.assertEquals(gp.getLatitude(), FastMath.PI/2, 1.0e-12);
376
377 point = new Vector3D(0.0, 110, 0);
378 direction = new Vector3D(0.0, 1.0, 0.0);
379 line = new Line(point, point.add(direction), 1.0e-10);
380 gp = model.getIntersectionPoint(line, point, frame, date);
381 Assert.assertEquals(gp.getLatitude(), 0, 1.0e-12);
382
383 }
384
385 @Test
386 public void testNoLineIntersection() {
387 AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
388 Frame frame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
389 OneAxisEllipsoid model = new OneAxisEllipsoid(100.0, 0.9, frame);
390 Vector3D point = new Vector3D(0.0, 93.7139699, 3.5930796);
391 Vector3D direction = new Vector3D(0.0, 9.0, -2.0);
392 Line line = new Line(point, point.add(direction), 1.0e-10);
393 Assert.assertNull(model.getIntersectionPoint(line, point, frame, date));
394 }
395
396 @Test
397 public void testNegativeZ() {
398 AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
399 Frame frame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
400 OneAxisEllipsoid model = new OneAxisEllipsoid(90.0, 5.0 / 9.0, frame);
401 Vector3D point = new Vector3D(140.0, 0.0, -30.0);
402 GeodeticPoint gp = model.transform(point, frame, date);
403 Vector3D rebuilt = model.transform(gp);
404 Assert.assertEquals(0.0, rebuilt.distance(point), 1.0e-10);
405 }
406
407 @Test
408 public void testNumerousIteration() {
409
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 Assert.assertEquals(0.0, rebuilt.distance(point), 8.0e-12);
421 }
422
423 @Test
424 public 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 Assert.assertEquals(0.0, rebuilt.distance(point), 1.0e-10);
433 }
434 }
435
436 @Test
437 public 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 Assert.assertEquals(0.0, rebuilt.distance(point), 1.0e-15 * point.getNorm());
445 }
446
447 @Test
448 public 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 Assert.assertEquals(0.0, rebuilt.distance(point), 1.0e-15 * point.getNorm());
458 }
459
460 @Test
461 public void testSerialization() throws IOException, ClassNotFoundException {
462 OneAxisEllipsoid original = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
463 Constants.WGS84_EARTH_FLATTENING,
464 FramesFactory.getITRFEquinox(IERSConventions.IERS_1996, true));
465 original.setAngularThreshold(1.0e-3);
466
467 ByteArrayOutputStream bos = new ByteArrayOutputStream();
468 ObjectOutputStream oos = new ObjectOutputStream(bos);
469 oos.writeObject(original);
470 Assert.assertTrue(bos.size() > 250);
471 Assert.assertTrue(bos.size() < 350);
472
473 ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
474 ObjectInputStream ois = new ObjectInputStream(bis);
475 OneAxisEllipsoid deserialized = (OneAxisEllipsoid) ois.readObject();
476 Assert.assertEquals(original.getEquatorialRadius(), deserialized.getEquatorialRadius(), 1.0e-12);
477 Assert.assertEquals(original.getFlattening(), deserialized.getFlattening(), 1.0e-12);
478
479 }
480
481 @Test
482 public void testIntersectionFromPoints() {
483 AbsoluteDate date = new AbsoluteDate(new DateComponents(2008, 03, 21),
484 TimeComponents.H12,
485 TimeScalesFactory.getUTC());
486
487 Frame frame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
488 OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.460, 1 / 298.257222101, frame);
489
490
491
492 final double mu = 3.9860047e14;
493 CircularOrbit circ =
494 new CircularOrbit(7178000.0, 0.5e-4, 0., FastMath.toRadians(90.), FastMath.toRadians(60.),
495 FastMath.toRadians(90.), PositionAngle.MEAN,
496 FramesFactory.getEME2000(), date, mu);
497
498
499 PVCoordinates pvSatEME2000 = circ.getPVCoordinates();
500 PVCoordinates pvSatItrf = frame.getTransformTo(FramesFactory.getEME2000(), date).transformPVCoordinates(pvSatEME2000);
501 Vector3D pSatItrf = pvSatItrf.getPosition();
502
503
504 GeodeticPoint geoPoint = new GeodeticPoint(FastMath.toRadians(70.), FastMath.toRadians(60.), 0.);
505 Vector3D pointItrf = earth.transform(geoPoint);
506 Line line = new Line(pSatItrf, pointItrf, 1.0e-10);
507 GeodeticPoint geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
508 Assert.assertEquals(geoPoint.getLongitude(), geoInter.getLongitude(), Utils.epsilonAngle);
509 Assert.assertEquals(geoPoint.getLatitude(), geoInter.getLatitude(), Utils.epsilonAngle);
510
511
512 geoPoint = new GeodeticPoint(FastMath.toRadians(65.), FastMath.toRadians(-120.), 0.);
513 pointItrf = earth.transform(geoPoint);
514 line = new Line(pSatItrf, pointItrf, 1.0e-10);
515 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
516 Assert.assertEquals(geoPoint.getLongitude(), geoInter.getLongitude(), Utils.epsilonAngle);
517 Assert.assertEquals(geoPoint.getLatitude(), geoInter.getLatitude(), Utils.epsilonAngle);
518
519
520 geoPoint = new GeodeticPoint(FastMath.toRadians(30.), FastMath.toRadians(60.), 0.);
521 pointItrf = earth.transform(geoPoint);
522 line = new Line(pSatItrf, pointItrf, 1.0e-10);
523
524 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
525
526
527 Assert.assertEquals(1.04437199, geoInter.getLongitude(), Utils.epsilonAngle);
528 Assert.assertEquals(1.36198012, geoInter.getLatitude(), Utils.epsilonAngle);
529
530
531
532 circ =
533 new CircularOrbit(7178000.0, 0.5e-4, 0., FastMath.toRadians(1.e-4), FastMath.toRadians(0.),
534 FastMath.toRadians(0.), PositionAngle.MEAN,
535 FramesFactory.getEME2000(), date, mu);
536
537
538 pvSatEME2000 = circ.getPVCoordinates();
539 pvSatItrf = frame.getTransformTo(FramesFactory.getEME2000(), date).transformPVCoordinates(pvSatEME2000);
540 pSatItrf = pvSatItrf.getPosition();
541
542
543 geoPoint = new GeodeticPoint(FastMath.toRadians(5.), FastMath.toRadians(0.), 0.);
544 pointItrf = earth.transform(geoPoint);
545 line = new Line(pSatItrf, pointItrf, 1.0e-10);
546 Assert.assertTrue(line.toSubSpace(pSatItrf).getX() < 0);
547 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
548 Assert.assertEquals(geoPoint.getLongitude(), geoInter.getLongitude(), Utils.epsilonAngle);
549 Assert.assertEquals(geoPoint.getLatitude(), geoInter.getLatitude(), Utils.epsilonAngle);
550
551
552 GeodeticPoint geoInter2 = earth.getIntersectionPoint(line, line.toSpace(new Vector1D(-line.toSubSpace(pSatItrf).getX())), frame, date);
553 Assert.assertTrue(FastMath.abs(geoInter.getLongitude() - geoInter2.getLongitude()) > FastMath.toRadians(0.1));
554 Assert.assertTrue(FastMath.abs(geoInter.getLatitude() - geoInter2.getLatitude()) > FastMath.toRadians(0.1));
555
556
557 geoPoint = new GeodeticPoint(FastMath.toRadians(-5.), FastMath.toRadians(0.), 0.);
558 pointItrf = earth.transform(geoPoint);
559 line = new Line(pSatItrf, pointItrf, 1.0e-10);
560 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
561 Assert.assertEquals(geoPoint.getLongitude(), geoInter.getLongitude(), Utils.epsilonAngle);
562 Assert.assertEquals(geoPoint.getLatitude(), geoInter.getLatitude(), Utils.epsilonAngle);
563
564
565 geoPoint = new GeodeticPoint(FastMath.toRadians(40.), FastMath.toRadians(0.), 0.);
566 pointItrf = earth.transform(geoPoint);
567 line = new Line(pSatItrf, pointItrf, 1.0e-10);
568 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
569 Assert.assertEquals(-0.00768481, geoInter.getLongitude(), Utils.epsilonAngle);
570 Assert.assertEquals( 0.32180410, geoInter.getLatitude(), Utils.epsilonAngle);
571
572
573
574
575 circ =
576 new CircularOrbit(7178000.0, 0.5e-4, 0., FastMath.toRadians(50.), FastMath.toRadians(0.),
577 FastMath.toRadians(90.), PositionAngle.MEAN,
578 FramesFactory.getEME2000(), date, mu);
579
580
581 pvSatEME2000 = circ.getPVCoordinates();
582 pvSatItrf = frame.getTransformTo(FramesFactory.getEME2000(), date).transformPVCoordinates(pvSatEME2000);
583 pSatItrf = pvSatItrf.getPosition();
584
585
586 geoPoint = new GeodeticPoint(FastMath.toRadians(40.), FastMath.toRadians(90.), 0.);
587 pointItrf = earth.transform(geoPoint);
588 line = new Line(pSatItrf, pointItrf, 1.0e-10);
589 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
590 Assert.assertEquals(geoPoint.getLongitude(), geoInter.getLongitude(), Utils.epsilonAngle);
591 Assert.assertEquals(geoPoint.getLatitude(), geoInter.getLatitude(), Utils.epsilonAngle);
592
593
594 geoPoint = new GeodeticPoint(FastMath.toRadians(60.), FastMath.toRadians(90.), 0.);
595 pointItrf = earth.transform(geoPoint);
596 line = new Line(pSatItrf, pointItrf, 1.0e-10);
597 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
598 Assert.assertEquals(geoPoint.getLongitude(), geoInter.getLongitude(), Utils.epsilonAngle);
599 Assert.assertEquals(geoPoint.getLatitude(), geoInter.getLatitude(), Utils.epsilonAngle);
600
601
602 geoPoint = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(90.), 0.);
603 pointItrf = earth.transform(geoPoint);
604 line = new Line(pSatItrf, pointItrf, 1.0e-10);
605 geoInter = earth.getIntersectionPoint(line, pSatItrf, frame, date);
606 Assert.assertEquals(FastMath.toRadians(89.5364061088196), geoInter.getLongitude(), Utils.epsilonAngle);
607 Assert.assertEquals(FastMath.toRadians(35.555543683351125), geoInter.getLatitude(), Utils.epsilonAngle);
608
609 }
610
611 @Test
612 public void testMovingGeodeticPointSymmetry() {
613
614 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
615 Constants.WGS84_EARTH_FLATTENING,
616 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
617 double lat0 = FastMath.toRadians(60.0);
618 double lon0 = FastMath.toRadians(25.0);
619 double alt0 = 100.0;
620 double lat1 = 1.0e-3;
621 double lon1 = -2.0e-3;
622 double alt1 = 1.2;
623 double lat2 = -1.0e-5;
624 double lon2 = -3.0e-5;
625 double alt2 = -0.01;
626 final DSFactory factory = new DSFactory(1, 2);
627 final DerivativeStructure latDS = factory.build(lat0, lat1, lat2);
628 final DerivativeStructure lonDS = factory.build(lon0, lon1, lon2);
629 final DerivativeStructure altDS = factory.build(alt0, alt1, alt2);
630
631
632 PVCoordinates pv = new PVCoordinates(earth.transform(new FieldGeodeticPoint<>(latDS, lonDS, altDS)));
633 FieldGeodeticPoint<DerivativeStructure> rebuilt = earth.transform(pv, earth.getBodyFrame(), null);
634 Assert.assertEquals(lat0, rebuilt.getLatitude().getReal(), 1.0e-16);
635 Assert.assertEquals(lat1, rebuilt.getLatitude().getPartialDerivative(1), 5.0e-19);
636 Assert.assertEquals(lat2, rebuilt.getLatitude().getPartialDerivative(2), 5.0e-14);
637 Assert.assertEquals(lon0, rebuilt.getLongitude().getReal(), 1.0e-16);
638 Assert.assertEquals(lon1, rebuilt.getLongitude().getPartialDerivative(1), 5.0e-19);
639 Assert.assertEquals(lon2, rebuilt.getLongitude().getPartialDerivative(2), 1.0e-20);
640 Assert.assertEquals(alt0, rebuilt.getAltitude().getReal(), 2.0e-11);
641 Assert.assertEquals(alt1, rebuilt.getAltitude().getPartialDerivative(1), 6.0e-13);
642 Assert.assertEquals(alt2, rebuilt.getAltitude().getPartialDerivative(2), 2.0e-14);
643
644 }
645
646 @Test
647 public void testMovingGeodeticPoint() {
648
649 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
650 Constants.WGS84_EARTH_FLATTENING,
651 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
652 double lat0 = FastMath.toRadians(60.0);
653 double lon0 = FastMath.toRadians(25.0);
654 double alt0 = 100.0;
655 double lat1 = 1.0e-3;
656 double lon1 = -2.0e-3;
657 double alt1 = 1.2;
658 double lat2 = -1.0e-5;
659 double lon2 = -3.0e-5;
660 double alt2 = -0.01;
661 final DSFactory factory = new DSFactory(1, 2);
662 final DerivativeStructure latDS = factory.build(lat0, lat1, lat2);
663 final DerivativeStructure lonDS = factory.build(lon0, lon1, lon2);
664 final DerivativeStructure altDS = factory.build(alt0, alt1, alt2);
665
666
667 PVCoordinates pv = new PVCoordinates(earth.transform(new FieldGeodeticPoint<>(latDS, lonDS, altDS)));
668
669
670 FiniteDifferencesDifferentiator differentiator =
671 new FiniteDifferencesDifferentiator(5, 0.1);
672 UnivariateDifferentiableFunction fx =
673 differentiator.differentiate(new UnivariateFunction() {
674 public double value(double dt) {
675 GeodeticPoint gp =
676 new GeodeticPoint(latDS.taylor(dt), lonDS.taylor(dt), altDS.taylor(dt));
677 return earth.transform(gp).getX();
678 }
679 });
680 UnivariateDifferentiableFunction fy =
681 differentiator.differentiate(new UnivariateFunction() {
682 public double value(double dt) {
683 GeodeticPoint gp =
684 new GeodeticPoint(latDS.taylor(dt), lonDS.taylor(dt), altDS.taylor(dt));
685 return earth.transform(gp).getY();
686 }
687 });
688 UnivariateDifferentiableFunction fz =
689 differentiator.differentiate(new UnivariateFunction() {
690 public double value(double dt) {
691 GeodeticPoint gp =
692 new GeodeticPoint(latDS.taylor(dt), lonDS.taylor(dt), altDS.taylor(dt));
693 return earth.transform(gp).getZ();
694 }
695 });
696 DerivativeStructure dtZero = factory.variable(0, 0.0);
697 DerivativeStructure xDS = fx.value(dtZero);
698 DerivativeStructure yDS = fy.value(dtZero);
699 DerivativeStructure zDS = fz.value(dtZero);
700 Assert.assertEquals(xDS.getValue(), pv.getPosition().getX(),
701 2.0e-20 * FastMath.abs(xDS.getValue()));
702 Assert.assertEquals(xDS.getPartialDerivative(1), pv.getVelocity().getX(),
703 2.0e-12 * FastMath.abs(xDS.getPartialDerivative(1)));
704 Assert.assertEquals(xDS.getPartialDerivative(2), pv.getAcceleration().getX(),
705 2.0e-9 * FastMath.abs(xDS.getPartialDerivative(2)));
706 Assert.assertEquals(yDS.getValue(), pv.getPosition().getY(),
707 2.0e-20 * FastMath.abs(yDS.getValue()));
708 Assert.assertEquals(yDS.getPartialDerivative(1), pv.getVelocity().getY(),
709 2.0e-12 * FastMath.abs(yDS.getPartialDerivative(1)));
710 Assert.assertEquals(yDS.getPartialDerivative(2), pv.getAcceleration().getY(),
711 2.0e-9 * FastMath.abs(yDS.getPartialDerivative(2)));
712 Assert.assertEquals(zDS.getValue(), pv.getPosition().getZ(),
713 2.0e-20 * FastMath.abs(zDS.getValue()));
714 Assert.assertEquals(zDS.getPartialDerivative(1), pv.getVelocity().getZ(),
715 2.0e-12 * FastMath.abs(zDS.getPartialDerivative(1)));
716 Assert.assertEquals(zDS.getPartialDerivative(2), pv.getAcceleration().getZ(),
717 2.0e-9 * FastMath.abs(zDS.getPartialDerivative(2)));
718 }
719
720 private void checkCartesianToEllipsoidic(double ae, double f,
721 double x, double y, double z,
722 double longitude, double latitude,
723 double altitude) {
724
725 AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
726 Frame frame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
727 OneAxisEllipsoid model = new OneAxisEllipsoid(ae, f, frame);
728 GeodeticPoint gp = model.transform(new Vector3D(x, y, z), frame, date);
729 Assert.assertEquals(longitude, MathUtils.normalizeAngle(gp.getLongitude(), longitude), 1.0e-10);
730 Assert.assertEquals(latitude, gp.getLatitude(), 1.0e-10);
731 Assert.assertEquals(altitude, gp.getAltitude(), 1.0e-10 * FastMath.abs(ae));
732 Vector3D rebuiltNadir = Vector3D.crossProduct(gp.getSouth(), gp.getWest());
733 Assert.assertEquals(0, rebuiltNadir.subtract(gp.getNadir()).getNorm(), 1.0e-15);
734
735 FieldGeodeticPoint<Decimal64> gp64 = model.transform(new FieldVector3D<Decimal64>(new Decimal64(x),
736 new Decimal64(y),
737 new Decimal64(z)),
738 frame,
739 new FieldAbsoluteDate<>(Decimal64Field.getInstance(), date));
740 Assert.assertEquals(longitude, MathUtils.normalizeAngle(gp64.getLongitude().getReal(), longitude), 1.0e-10);
741 Assert.assertEquals(latitude, gp64.getLatitude().getReal(), 1.0e-10);
742 Assert.assertEquals(altitude, gp64.getAltitude().getReal(), 1.0e-10 * FastMath.abs(ae));
743 FieldVector3D<Decimal64> rebuiltNadir64 = FieldVector3D.crossProduct(gp64.getSouth(), gp64.getWest());
744 Assert.assertEquals(0, rebuiltNadir64.subtract(gp64.getNadir()).getNorm().getReal(), 1.0e-15);
745
746
747 gp = model.transform(model.projectToGround(new Vector3D(x, y, z), date, frame), frame, date);
748 Assert.assertEquals(longitude, MathUtils.normalizeAngle(gp.getLongitude(), longitude), 1.0e-10);
749 Assert.assertEquals(latitude, gp.getLatitude(), 1.0e-10);
750 Assert.assertEquals(0.0, gp.getAltitude(), 1.0e-10 * FastMath.abs(ae));
751
752
753 FieldGeodeticPoint<DerivativeStructure> gpDs = model.transform(
754 model.projectToGround(
755 new TimeStampedPVCoordinates(
756 date,
757 new Vector3D(x, y, z),
758 new Vector3D(1, 2, 3)),
759 frame),
760 frame,
761 date);
762 Assert.assertEquals(longitude, MathUtils.normalizeAngle(gpDs.getLongitude().getReal(), longitude), 1.0e-10);
763 Assert.assertEquals(latitude, gpDs.getLatitude().getReal(), 1.0e-10);
764 Assert.assertEquals(0.0, gpDs.getAltitude().getReal(), 1.0e-10 * FastMath.abs(ae));
765
766 }
767
768 @Test
769 public void testTransformVsOldIterativeSobol() {
770
771 OneAxisEllipsoid model = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
772 Constants.WGS84_EARTH_FLATTENING,
773 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
774 SobolSequenceGenerator sobol = new SobolSequenceGenerator(3);
775 final double rMax = 10 * model.getEquatorialRadius();
776 Stream<Vector3D> points = Stream.generate(() -> {
777 final double[] v = sobol.nextVector();
778 return new Vector3D(rMax * (2 * v[0] - 1), rMax * (2 * v[1] - 1), rMax * (2 * v[2] - 1));
779 }).limit(1000000);
780
781 doTestTransformVsOldIterative(model, points, 2.0e-15, 1.0e-15, 8.0e-14 * rMax);
782 }
783
784 @Test
785 public void testTransformVsOldIterativePolarAxis() {
786 OneAxisEllipsoid model = new OneAxisEllipsoid(90, 5.0 / 9.0,
787 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
788 Stream<Vector3D> points = DoubleStream.iterate(0, x -> x + 1.0).limit(150).mapToObj(z -> new Vector3D(0, 0, z));
789 doTestTransformVsOldIterative(model, points, 2.0e-15, 1.0e-15, 1.0e-14 * model.getEquatorialRadius());
790 }
791
792 @Test
793 public void testTransformVsOldIterativeEquatorial() {
794 OneAxisEllipsoid model = new OneAxisEllipsoid(90, 5.0 / 9.0,
795 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
796 Stream<Vector3D> points = DoubleStream.iterate(0, x -> x + 1.0).limit(150).mapToObj(x -> new Vector3D(x, 0, 0));
797 doTestTransformVsOldIterative(model, points, 2.0e-15, 1.0e-15, 1.0e-14 * model.getEquatorialRadius());
798 }
799
800 @Test
801 public void testIssue373() {
802 final Frame ecef = FramesFactory.getITRF(IERSConventions.IERS_2010,true);
803 final OneAxisEllipsoid earth = new OneAxisEllipsoid(6378137, 1./298.257223563, ecef);
804 final Vector3D sunPos = new Vector3D(-149757851422.23358, 8410610314.781021, 14717269835.161688 );
805 final GeodeticPoint sunGP = earth.transform(sunPos, ecef, null);
806 Assert.assertEquals(5.603878, FastMath.toDegrees(sunGP.getLatitude()), 1.0e-6);
807 }
808
809 private void doTestTransformVsOldIterative(OneAxisEllipsoid model,
810 Stream<Vector3D> points,
811 double latitudeTolerance, double longitudeTolerance,
812 double altitudeTolerance) {
813 points.forEach(point -> {
814 try {
815 GeodeticPoint reference = transformOldIterative(model, point, model.getBodyFrame(), null);
816 GeodeticPoint result = model.transform(point, model.getBodyFrame(), null);
817 Assert.assertEquals(reference.getLatitude(), result.getLatitude(), latitudeTolerance);
818 Assert.assertEquals(reference.getLongitude(), result.getLongitude(), longitudeTolerance);
819 Assert.assertEquals(reference.getAltitude(), result.getAltitude(), altitudeTolerance);
820 } catch (OrekitException oe) {
821 Assert.fail(oe.getLocalizedMessage());
822 }
823 });
824
825 }
826
827
828
829
830
831
832
833
834
835
836
837
838
839 private GeodeticPoint transformOldIterative(final OneAxisEllipsoid model,
840 final Vector3D point,
841 final Frame frame,
842 final AbsoluteDate date) {
843
844
845 final Vector3D pointInBodyFrame = frame.getTransformTo(model.getBodyFrame(), date).transformPosition(point);
846 final double r2 = pointInBodyFrame.getX() * pointInBodyFrame.getX() +
847 pointInBodyFrame.getY() * pointInBodyFrame.getY();
848 final double r = FastMath.sqrt(r2);
849 final double z = pointInBodyFrame.getZ();
850
851
852 final Ellipse meridian = new Ellipse(Vector3D.ZERO,
853 new Vector3D(pointInBodyFrame.getX() / r, pointInBodyFrame.getY() / r, 0),
854 Vector3D.PLUS_K,
855 model.getA(), model.getC(), model.getBodyFrame());
856
857
858 final Vector2D ellipsePoint = meridian.projectToEllipse(new Vector2D(r, z));
859
860
861 final double dr = r - ellipsePoint.getX();
862 final double dz = z - ellipsePoint.getY();
863 final double ae2 = model.getA() * model.getA();
864 final double f = model.getFlattening();
865 final double g = 1.0 - f;
866 final double g2 = g * g;
867 final double insideIfNegative = g2 * (r2 - ae2) + z * z;
868
869 return new GeodeticPoint(FastMath.atan2(ellipsePoint.getY(), g2 * ellipsePoint.getX()),
870 FastMath.atan2(pointInBodyFrame.getY(), pointInBodyFrame.getX()),
871 FastMath.copySign(FastMath.hypot(dr, dz), insideIfNegative));
872
873 }
874
875 @Before
876 public void setUp() {
877 Utils.setDataRoot("regular-data");
878 }
879
880 }
881