1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.frames;
18
19 import java.io.IOException;
20 import java.util.stream.Stream;
21
22 import org.hipparchus.CalculusFieldElement;
23 import org.hipparchus.Field;
24 import org.hipparchus.analysis.polynomials.PolynomialFunction;
25 import org.hipparchus.complex.Complex;
26 import org.hipparchus.complex.ComplexField;
27 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
28 import org.hipparchus.geometry.euclidean.threed.Vector3D;
29 import org.hipparchus.random.RandomGenerator;
30 import org.hipparchus.random.Well1024a;
31 import org.hipparchus.util.Binary64Field;
32 import org.hipparchus.util.Binary64;
33 import org.hipparchus.util.FastMath;
34 import org.hipparchus.util.MathUtils;
35 import org.hipparchus.util.SinCos;
36 import org.hipparchus.util.FieldSinCos;
37 import org.junit.jupiter.api.AfterEach;
38 import org.junit.jupiter.api.Assertions;
39 import org.junit.jupiter.api.BeforeEach;
40 import org.junit.jupiter.api.Test;
41 import org.junit.jupiter.params.ParameterizedTest;
42 import org.junit.jupiter.params.provider.Arguments;
43 import org.junit.jupiter.params.provider.MethodSource;
44 import org.orekit.Utils;
45 import org.orekit.bodies.BodyShape;
46 import org.orekit.bodies.FieldGeodeticPoint;
47 import org.orekit.bodies.GeodeticPoint;
48 import org.orekit.bodies.OneAxisEllipsoid;
49 import org.orekit.orbits.CircularOrbit;
50 import org.orekit.orbits.FieldCircularOrbit;
51 import org.orekit.orbits.PositionAngleType;
52 import org.orekit.propagation.FieldSpacecraftState;
53 import org.orekit.propagation.SpacecraftState;
54 import org.orekit.propagation.analytical.FieldKeplerianPropagator;
55 import org.orekit.propagation.analytical.KeplerianPropagator;
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.Constants;
62 import org.orekit.utils.FieldPVCoordinates;
63 import org.orekit.utils.FieldTrackingCoordinates;
64 import org.orekit.utils.IERSConventions;
65 import org.orekit.utils.PVCoordinates;
66 import org.orekit.utils.TrackingCoordinates;
67
68
69 class TopocentricFrameTest {
70
71
72 private AbsoluteDate date;
73
74
75 private Frame itrf;
76
77
78 OneAxisEllipsoid earthSpheric;
79
80
81 private double mu;
82
83
84 @Test
85 void testZero() {
86
87 final GeodeticPoint point = new GeodeticPoint(0., 0., 0.);
88 final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "zero");
89
90
91 final double xDiff = Vector3D.dotProduct(topoFrame.getEast(), Vector3D.PLUS_J);
92 final double yDiff = Vector3D.dotProduct(topoFrame.getNorth(), Vector3D.PLUS_K);
93 final double zDiff = Vector3D.dotProduct(topoFrame.getZenith(), Vector3D.PLUS_I);
94 Assertions.assertEquals(1., xDiff, Utils.epsilonTest);
95 Assertions.assertEquals(1., yDiff, Utils.epsilonTest);
96 Assertions.assertEquals(1., zDiff, Utils.epsilonTest);
97 }
98
99 @Test
100 void testPole() {
101
102 final GeodeticPoint point = new GeodeticPoint(FastMath.PI/2., 0., 0.);
103 final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "north pole");
104
105
106 final double xDiff = Vector3D.dotProduct(topoFrame.getEast(), Vector3D.PLUS_J);
107 final double yDiff = Vector3D.dotProduct(topoFrame.getSouth(), Vector3D.PLUS_I);
108 final double zDiff = Vector3D.dotProduct(topoFrame.getZenith(), Vector3D.PLUS_K);
109 Assertions.assertEquals(1., xDiff, Utils.epsilonTest);
110 Assertions.assertEquals(1., yDiff, Utils.epsilonTest);
111 Assertions.assertEquals(1., zDiff, Utils.epsilonTest);
112 }
113
114 @Test
115 void testNormalLatitudes() {
116
117
118 final GeodeticPoint point1 = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(30.), 0.);
119 final TopocentricFrame topoFrame1 = new TopocentricFrame(earthSpheric, point1, "lat 45");
120 Assertions.assertEquals(0.0,
121 Vector3D.distance(topoFrame1.getCartesianPoint(), earthSpheric.transform(point1)),
122 7.0e-9);
123
124
125 final GeodeticPoint point2 = new GeodeticPoint(FastMath.toRadians(-45.), FastMath.toRadians(30.), 0.);
126 final TopocentricFrame topoFrame2 = new TopocentricFrame(earthSpheric, point2, "lat -45");
127 Assertions.assertEquals(0.0,
128 Vector3D.distance(topoFrame2.getCartesianPoint(), earthSpheric.transform(point2)),
129 2.0e-9);
130
131
132 final double xDiff = Vector3D.dotProduct(topoFrame1.getEast(), topoFrame2.getEast());
133 final double yDiff = Vector3D.dotProduct(topoFrame1.getNorth(), topoFrame2.getNorth());
134 final double zDiff = Vector3D.dotProduct(topoFrame1.getZenith(), topoFrame2.getZenith());
135
136 Assertions.assertEquals(1., xDiff, Utils.epsilonTest);
137 Assertions.assertEquals(0., yDiff, Utils.epsilonTest);
138 Assertions.assertEquals(0., zDiff, Utils.epsilonTest);
139 }
140
141 @Test
142 void testOppositeLongitudes() {
143
144
145 final GeodeticPoint point1 = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(30.), 0.);
146 final TopocentricFrame topoFrame1 = new TopocentricFrame(earthSpheric, point1, "lon 30");
147 final GeodeticPoint p1 = topoFrame1.getPoint();
148 Assertions.assertEquals(point1.getLatitude(), p1.getLatitude(), 1.0e-15);
149 Assertions.assertEquals(point1.getLongitude(), p1.getLongitude(), 1.0e-15);
150 Assertions.assertEquals(point1.getAltitude(), p1.getAltitude(), 1.0e-15);
151
152
153 final GeodeticPoint point2 = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(210.), 0.);
154 final TopocentricFrame topoFrame2 = new TopocentricFrame(earthSpheric, point2, "lon 210");
155
156
157
158 final double xDiff = Vector3D.dotProduct(topoFrame1.getEast(), topoFrame2.getWest());
159 final double yDiff = Vector3D.dotProduct(topoFrame1.getNorth(), topoFrame2.getNorth());
160 final double zDiff = Vector3D.dotProduct(topoFrame1.getZenith(), topoFrame2.getZenith());
161
162 Assertions.assertEquals(1., xDiff, Utils.epsilonTest);
163 Assertions.assertEquals(0., yDiff, Utils.epsilonTest);
164 Assertions.assertEquals(0., zDiff, Utils.epsilonTest);
165 }
166
167 @Test
168 void testAntipodes() {
169
170
171 final GeodeticPoint point1 = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(30.), 0.);
172 final TopocentricFrame topoFrame1 = new TopocentricFrame(earthSpheric, point1, "lon 30");
173
174
175 final GeodeticPoint point2 = new GeodeticPoint(FastMath.toRadians(-45.), FastMath.toRadians(210.), 0.);
176 final TopocentricFrame topoFrame2 = new TopocentricFrame(earthSpheric, point2, "lon 210");
177
178
179
180 final double xDiff = Vector3D.dotProduct(topoFrame1.getEast(), topoFrame2.getWest());
181 final double yDiff = Vector3D.dotProduct(topoFrame1.getNorth(), topoFrame2.getNorth());
182 final double zDiff = Vector3D.dotProduct(topoFrame1.getZenith(), topoFrame2.getZenith());
183
184 Assertions.assertEquals(1., xDiff, Utils.epsilonTest);
185 Assertions.assertEquals(1., yDiff, Utils.epsilonTest);
186 Assertions.assertEquals(-1., zDiff, Utils.epsilonTest);
187
188 Assertions.assertEquals(1, Vector3D.dotProduct(topoFrame1.getNadir(), topoFrame2.getZenith()), Utils.epsilonTest);
189 Assertions.assertEquals(1, Vector3D.dotProduct(topoFrame1.getZenith(), topoFrame2.getNadir()), Utils.epsilonTest);
190
191 }
192
193 @Test
194 void testSiteAtZenith() {
195
196
197 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(30.), 0.);
198 final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 30 lat 45");
199
200
201 final GeodeticPoint satPoint = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(30.), 800000.);
202 final TrackingCoordinates tc = topoFrame.getTrackingCoordinates(earthSpheric.transform(satPoint),
203 earthSpheric.getBodyFrame(),
204 date);
205
206
207 Assertions.assertEquals(FastMath.PI/2., tc.getElevation(), Utils.epsilonAngle);
208
209
210 Assertions.assertEquals(800000., tc.getRange(), 1e-8);
211 }
212
213 @Test
214 void testFieldSiteAtZenith() {
215 doTestFieldSiteAtZenith(Binary64Field.getInstance());
216 }
217
218 private <T extends CalculusFieldElement<T>> void doTestFieldSiteAtZenith(final Field<T> field) {
219
220
221 final T zero = field.getZero();
222
223
224 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(30.), 0.);
225 final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 30 lat 45");
226
227
228 final FieldGeodeticPoint<T> satPoint = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(45.)),
229 zero.add(FastMath.toRadians(30.)),
230 zero.add(800000.));
231 final FieldAbsoluteDate<T> fieldDate = new FieldAbsoluteDate<>(field, date);
232 final FieldTrackingCoordinates<T> tc = topoFrame.getTrackingCoordinates(earthSpheric.transform(satPoint),
233 earthSpheric.getBodyFrame(),
234 fieldDate);
235
236
237 Assertions.assertEquals(FastMath.PI/2., tc.getElevation().getReal(), Utils.epsilonAngle);
238
239
240 Assertions.assertEquals(800000., tc.getRange().getReal(), 1e-8);
241
242 }
243
244 @Test
245 void testAzimuthEquatorial() {
246
247
248 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(30.), 0.);
249 final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 30 lat 0");
250
251
252
253 GeodeticPoint infPoint = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(50.), 1000000000.);
254 TrackingCoordinates tc = topoFrame.getTrackingCoordinates(earthSpheric.transform(infPoint), earthSpheric.getBodyFrame(), date);
255
256
257 Assertions.assertEquals(FastMath.PI/2., tc.getAzimuth(), Utils.epsilonAngle);
258
259
260 Assertions.assertEquals(FastMath.PI/2. - FastMath.abs(point.getLongitude() - infPoint.getLongitude()),
261 tc.getElevation(),
262 1.e-2);
263
264
265
266 infPoint = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(10.), 1000000000.);
267 tc = topoFrame.getTrackingCoordinates(earthSpheric.transform(infPoint), earthSpheric.getBodyFrame(), date);
268
269
270 Assertions.assertEquals(3*FastMath.PI/2., tc.getAzimuth(), Utils.epsilonAngle);
271
272
273 Assertions.assertEquals(FastMath.PI/2. - FastMath.abs(point.getLongitude() - infPoint.getLongitude()),
274 tc.getElevation(),
275 1.e-2);
276
277 }
278
279 @Test
280 void testFieldAzimuthEquatorial() {
281 doTestFieldAzimuthEquatorial(Binary64Field.getInstance());
282 }
283
284 private <T extends CalculusFieldElement<T>> void doTestFieldAzimuthEquatorial(final Field<T> field) {
285
286
287 final T zero = field.getZero();
288
289
290 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(30.), 0.);
291 final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 30 lat 0");
292
293
294
295 FieldGeodeticPoint<T> infPoint = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(0.)),
296 zero.add(FastMath.toRadians(50.)),
297 zero.add(1000000000.));
298
299
300 final FieldAbsoluteDate<T> fieldDate = new FieldAbsoluteDate<>(field, date);
301 FieldTrackingCoordinates<T> tc = topoFrame.getTrackingCoordinates(earthSpheric.transform(infPoint),
302 earthSpheric.getBodyFrame(),
303 fieldDate);
304
305 Assertions.assertEquals(FastMath.PI/2., tc.getAzimuth().getReal(), Utils.epsilonAngle);
306
307
308 Assertions.assertEquals(FastMath.abs(infPoint.getLongitude().negate().add(point.getLongitude())).negate().add(FastMath.PI/2.).getReal(),
309 tc.getElevation().getReal(),
310 1.e-2);
311
312
313
314 infPoint = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(0.)),
315 zero.add(FastMath.toRadians(10.)),
316 zero.add(1000000000.));
317
318
319 tc = topoFrame.getTrackingCoordinates(earthSpheric.transform(infPoint), earthSpheric.getBodyFrame(), fieldDate);
320 Assertions.assertEquals(3*FastMath.PI/2., tc.getAzimuth().getReal(), Utils.epsilonAngle);
321
322
323 Assertions.assertEquals(FastMath.abs(infPoint.getLongitude().negate().add(point.getLongitude())).negate().add(FastMath.PI/2.).getReal(),
324 tc.getElevation().getReal(),
325 1.e-2);
326
327 }
328
329 @Test
330 void testAzimuthPole() {
331
332
333 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(89.999), FastMath.toRadians(0.), 0.);
334 final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 0 lat 90");
335
336
337
338 GeodeticPoint satPoint = new GeodeticPoint(FastMath.toRadians(28.), FastMath.toRadians(30.), 800000.);
339 TrackingCoordinates tc = topoFrame.getTrackingCoordinates(earthSpheric.transform(satPoint), earthSpheric.getBodyFrame(), date);
340
341 Assertions.assertEquals(FastMath.PI - satPoint.getLongitude(), tc.getAzimuth(), 1.e-5);
342
343
344
345 satPoint = new GeodeticPoint(FastMath.toRadians(28.), FastMath.toRadians(-30.), 800000.);
346 tc = topoFrame.getTrackingCoordinates(earthSpheric.transform(satPoint), earthSpheric.getBodyFrame(), date);
347
348 Assertions.assertEquals(FastMath.PI - satPoint.getLongitude(), tc.getAzimuth(), 1.e-5);
349
350 }
351
352 @Test
353 void testFieldAzimuthPole() {
354 doTestFieldAzimuthPole(Binary64Field.getInstance());
355 }
356
357 private <T extends CalculusFieldElement<T>> void doTestFieldAzimuthPole(final Field<T> field) {
358
359
360 final T zero = field.getZero();
361
362
363 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(89.999), FastMath.toRadians(0.), 0.);
364 final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 0 lat 90");
365
366
367
368 FieldGeodeticPoint<T> satPoint = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(28.)),
369 zero.add(FastMath.toRadians(30.)), zero.add(800000.));
370 final FieldAbsoluteDate<T> fieldDate = new FieldAbsoluteDate<>(field, date);
371 FieldTrackingCoordinates<T> tc = topoFrame.getTrackingCoordinates(earthSpheric.transform(satPoint), earthSpheric.getBodyFrame(), fieldDate);
372 Assertions.assertEquals(satPoint.getLongitude().negate().add(FastMath.PI).getReal(),
373 tc.getAzimuth().getReal(), 1.e-5);
374
375
376
377 satPoint = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(28.)),
378 zero.add(FastMath.toRadians(-30.)), zero.add(800000.));
379 tc = topoFrame.getTrackingCoordinates(earthSpheric.transform(satPoint), earthSpheric.getBodyFrame(), fieldDate);
380
381 Assertions.assertEquals(satPoint.getLongitude().negate().add(FastMath.PI).getReal(),
382 tc.getAzimuth().getReal(), 1.e-5);
383
384 }
385
386 @Test
387 void testDoppler() {
388
389
390 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(5.), 0.);
391 final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 5 lat 45");
392
393
394
395 final CircularOrbit orbit =
396 new CircularOrbit(7178000.0, 0.5e-8, -0.5e-8, FastMath.toRadians(50.), FastMath.toRadians(120.),
397 FastMath.toRadians(90.), PositionAngleType.MEAN,
398 FramesFactory.getEME2000(), date, mu);
399
400
401 final Transform eme2000ToItrf = FramesFactory.getEME2000().getTransformTo(earthSpheric.getBodyFrame(), date);
402 final PVCoordinates pvSatItrf = eme2000ToItrf.transformPVCoordinates(orbit.getPVCoordinates());
403
404
405
406 final double dop = topoFrame.getRangeRate(pvSatItrf, earthSpheric.getBodyFrame(), date);
407
408
409
410 final double dt = 0.1;
411 KeplerianPropagator extrapolator = new KeplerianPropagator(orbit);
412
413
414 AbsoluteDate dateP = date.shiftedBy(dt);
415 Transform j2000ToItrfP = FramesFactory.getEME2000().getTransformTo(earthSpheric.getBodyFrame(), dateP);
416 SpacecraftState orbitP = extrapolator.propagate(dateP);
417 Vector3D satPointGeoP = j2000ToItrfP.transformPosition(orbitP.getPosition());
418
419
420 AbsoluteDate dateM = date.shiftedBy(-dt);
421 Transform j2000ToItrfM = FramesFactory.getEME2000().getTransformTo(earthSpheric.getBodyFrame(), dateM);
422 SpacecraftState orbitM = extrapolator.propagate(dateM);
423 Vector3D satPointGeoM = j2000ToItrfM.transformPosition(orbitM.getPosition());
424
425
426 double rangeP = topoFrame.getTrackingCoordinates(satPointGeoP, earthSpheric.getBodyFrame(), dateP).getRange();
427 double rangeM = topoFrame.getTrackingCoordinates(satPointGeoM, earthSpheric.getBodyFrame(), dateM).getRange();
428 final double dopRef2 = (rangeP - rangeM) / (2. * dt);
429 Assertions.assertEquals(dopRef2, dop, 1.e-3);
430
431 }
432
433 @Test
434 void testFieldDoppler() {
435 doTestFieldDoppler(Binary64Field.getInstance());
436 }
437
438 private <T extends CalculusFieldElement<T>> void doTestFieldDoppler(final Field<T> field) {
439
440
441 final T zero = field.getZero();
442
443
444 final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(5.), 0.);
445 final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 5 lat 45");
446
447
448 final FieldAbsoluteDate<T> fieldDate = new FieldAbsoluteDate<>(field, date);
449
450
451
452 final FieldCircularOrbit<T> orbit =
453 new FieldCircularOrbit<>(zero.add(7178000.0), zero.add(0.5e-8), zero.add(-0.5e-8),
454 zero.add(FastMath.toRadians(50.)), zero.add(FastMath.toRadians(120.)),
455 zero.add(FastMath.toRadians(90.)), PositionAngleType.MEAN,
456 FramesFactory.getEME2000(), fieldDate, zero.add(mu));
457
458
459 final FieldTransform<T> eme2000ToItrf = FramesFactory.getEME2000().getTransformTo(earthSpheric.getBodyFrame(), fieldDate);
460 final FieldPVCoordinates<T> pvSatItrf = eme2000ToItrf.transformPVCoordinates(orbit.getPVCoordinates());
461
462
463
464 final T dop = topoFrame.getRangeRate(pvSatItrf, earthSpheric.getBodyFrame(), fieldDate);
465
466
467
468 final double dt = 0.1;
469 FieldKeplerianPropagator<T> extrapolator = new FieldKeplerianPropagator<>(orbit);
470
471
472 FieldAbsoluteDate<T> dateP = fieldDate.shiftedBy(dt);
473 FieldTransform<T> j2000ToItrfP = FramesFactory.getEME2000().getTransformTo(earthSpheric.getBodyFrame(), dateP);
474 FieldSpacecraftState<T> orbitP = extrapolator.propagate(dateP);
475 FieldVector3D<T> satPointGeoP = j2000ToItrfP.transformPosition(orbitP.getPosition());
476
477
478 FieldAbsoluteDate<T> dateM = fieldDate.shiftedBy(-dt);
479 FieldTransform<T> j2000ToItrfM = FramesFactory.getEME2000().getTransformTo(earthSpheric.getBodyFrame(), dateM);
480 FieldSpacecraftState<T> orbitM = extrapolator.propagate(dateM);
481 FieldVector3D<T> satPointGeoM = j2000ToItrfM.transformPosition(orbitM.getPosition());
482
483
484 T rangeP = topoFrame.getTrackingCoordinates(satPointGeoP, earthSpheric.getBodyFrame(), dateP).getRange();
485 T rangeM = topoFrame.getTrackingCoordinates(satPointGeoM, earthSpheric.getBodyFrame(), dateM).getRange();
486 final T dopRef2 = (rangeP.subtract(rangeM)).divide(2. * dt);
487 Assertions.assertEquals(dopRef2.getReal(), dop.getReal(), 1.e-3);
488
489 }
490
491 @Test
492 void testEllipticEarth() {
493
494
495 final OneAxisEllipsoid earthElliptic =
496 new OneAxisEllipsoid(6378136.460, 1 / 298.257222101, itrf);
497
498
499
500 final GeodeticPoint satPointGeo = new GeodeticPoint(FastMath.toRadians(30.), FastMath.toRadians(15.), 800000.);
501 final Vector3D satPoint = earthElliptic.transform(satPointGeo);
502
503
504
505
506 GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(5.), 0.);
507 TopocentricFrame topoElliptic = new TopocentricFrame(earthElliptic, point, "elliptic, equatorial lon 5");
508 TopocentricFrame topoSpheric = new TopocentricFrame(earthSpheric, point, "spheric, equatorial lon 5");
509 TrackingCoordinates tcElliptic = topoElliptic.getTrackingCoordinates(satPoint, earthElliptic.getBodyFrame(), date);
510 TrackingCoordinates tcSpheric = topoSpheric.getTrackingCoordinates(satPoint, earthElliptic.getBodyFrame(), date);
511
512
513
514 Assertions.assertEquals(tcElliptic.getAzimuth(), tcSpheric.getAzimuth(), Utils.epsilonAngle);
515 Assertions.assertEquals(tcElliptic.getElevation(), tcSpheric.getElevation(), Utils.epsilonAngle);
516 Assertions.assertEquals(tcElliptic.getRange(), tcSpheric.getRange(), Utils.epsilonTest);
517
518
519
520 GeodeticPoint infPointGeo = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(-15.), 1000000000.);
521 Vector3D infPoint = earthElliptic.transform(infPointGeo);
522 tcElliptic = topoElliptic.getTrackingCoordinates(infPoint, earthElliptic.getBodyFrame(), date);
523
524
525 Assertions.assertEquals(3*FastMath.PI/2., tcElliptic.getAzimuth(), Utils.epsilonAngle);
526
527
528 Assertions.assertEquals(FastMath.PI/2. - FastMath.abs(point.getLongitude() - infPointGeo.getLongitude()),
529 tcElliptic.getElevation(), 1.e-2);
530
531
532
533 infPointGeo = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(25.), 1000000000.);
534 infPoint = earthElliptic.transform(infPointGeo);
535 tcElliptic = topoElliptic.getTrackingCoordinates(infPoint, earthElliptic.getBodyFrame(), date);
536
537
538 Assertions.assertEquals(FastMath.PI/2., tcElliptic.getAzimuth(), Utils.epsilonAngle);
539
540
541 Assertions.assertEquals(FastMath.PI/2. - FastMath.abs(point.getLongitude() - infPointGeo.getLongitude()),
542 tcElliptic.getElevation(), 1.e-2);
543
544
545
546
547 point = new GeodeticPoint(FastMath.toRadians(89.999), FastMath.toRadians(0.), 0.);
548 topoSpheric = new TopocentricFrame(earthSpheric, point, "lon 0 lat 90");
549 topoElliptic = new TopocentricFrame(earthElliptic, point, "lon 0 lat 90");
550
551
552
553 tcElliptic = topoElliptic.getTrackingCoordinates(satPoint, earthElliptic.getBodyFrame(), date);
554 tcSpheric = topoSpheric.getTrackingCoordinates(satPoint, earthSpheric.getBodyFrame(), date);
555 Assertions.assertEquals(tcElliptic.getAzimuth(), tcSpheric.getAzimuth(), 1.e-7);
556 Assertions.assertEquals(tcElliptic.getElevation(), tcSpheric.getElevation(), 1.e-2);
557 Assertions.assertEquals(tcElliptic.getRange(), tcSpheric.getRange(), 20.e+3);
558
559
560
561
562
563 point = new GeodeticPoint(FastMath.toRadians(60), FastMath.toRadians(30.), 0.);
564 topoSpheric = new TopocentricFrame(earthSpheric, point, "lon 10 lat 45");
565 topoElliptic = new TopocentricFrame(earthElliptic, point, "lon 10 lat 45");
566
567
568
569 tcElliptic = topoElliptic.getTrackingCoordinates(satPoint, earthElliptic.getBodyFrame(), date);
570 tcSpheric = topoSpheric.getTrackingCoordinates(satPoint, earthSpheric.getBodyFrame(), date);
571 Assertions.assertEquals(tcElliptic.getAzimuth(), tcSpheric.getAzimuth(), 1.e-2);
572 Assertions.assertEquals(tcElliptic.getElevation(), tcSpheric.getElevation(), 1.e-2);
573 Assertions.assertEquals(tcElliptic.getRange(), tcSpheric.getRange(), 20.e+3);
574
575 }
576
577 @Test
578 void testFieldEllipticEarth() {
579 doTestFieldEllipticEarth(Binary64Field.getInstance());
580 }
581
582 private <T extends CalculusFieldElement<T>> void doTestFieldEllipticEarth(final Field<T> field) {
583
584
585 final T zero = field.getZero();
586
587
588 final OneAxisEllipsoid earthElliptic =
589 new OneAxisEllipsoid(6378136.460, 1 / 298.257222101, itrf);
590
591
592 final FieldAbsoluteDate<T> fieldDate = new FieldAbsoluteDate<>(field, date);
593
594
595
596 final FieldGeodeticPoint<T> satPointGeo = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(30.)),
597 zero.add(FastMath.toRadians(15.)),
598 zero.add(800000.));
599 final FieldVector3D<T> satPoint = earthElliptic.transform(satPointGeo);
600
601
602
603
604 GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(5.), 0.);
605 TopocentricFrame topoElliptic = new TopocentricFrame(earthElliptic, point, "elliptic, equatorial lon 5");
606 TopocentricFrame topoSpheric = new TopocentricFrame(earthSpheric, point, "spheric, equatorial lon 5");
607
608
609
610 FieldTrackingCoordinates<T> tcElliptic =
611 topoElliptic.getTrackingCoordinates(satPoint, earthElliptic.getBodyFrame(), fieldDate);
612 FieldTrackingCoordinates<T> tcSpheric =
613 topoSpheric.getTrackingCoordinates(satPoint, earthSpheric.getBodyFrame(), fieldDate);
614 Assertions.assertEquals(tcElliptic.getAzimuth().getReal(), tcSpheric.getAzimuth().getReal(), Utils.epsilonAngle);
615 Assertions.assertEquals(tcElliptic.getElevation().getReal(), tcSpheric.getElevation().getReal(), Utils.epsilonAngle);
616 Assertions.assertEquals(tcElliptic.getRange().getReal(), tcSpheric.getRange().getReal(), Utils.epsilonTest);
617
618
619
620 FieldGeodeticPoint<T> infPointGeo = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(0.)),
621 zero.add(FastMath.toRadians(-15.)),
622 zero.add(1000000000.));
623 FieldVector3D<T> infPoint = earthElliptic.transform(infPointGeo);
624 tcElliptic = topoElliptic.getTrackingCoordinates(infPoint, earthElliptic.getBodyFrame(), fieldDate);
625
626
627 Assertions.assertEquals(3*FastMath.PI/2., tcElliptic.getAzimuth().getReal(), Utils.epsilonAngle);
628
629
630 Assertions.assertEquals(FastMath.abs(infPointGeo.getLongitude().negate().add(point.getLongitude())).negate().add(FastMath.PI/2.).getReal(),
631 tcElliptic.getElevation().getReal(),
632 1.e-2);
633
634
635
636 infPointGeo = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(0.)),
637 zero.add(FastMath.toRadians(25.)),
638 zero.add(1000000000.));
639 infPoint = earthElliptic.transform(infPointGeo);
640 tcElliptic = topoElliptic.getTrackingCoordinates(infPoint, earthElliptic.getBodyFrame(), fieldDate);
641
642
643 Assertions.assertEquals(FastMath.PI/2., tcElliptic.getAzimuth().getReal(), Utils.epsilonAngle);
644
645
646 Assertions.assertEquals(FastMath.abs(infPointGeo.getLongitude().negate().add(point.getLongitude())).negate().add(FastMath.PI/2.).getReal(),
647 tcElliptic.getElevation().getReal(),
648 1.e-2);
649
650
651
652
653 point = new GeodeticPoint(FastMath.toRadians(89.999), FastMath.toRadians(0.), 0.);
654 topoSpheric = new TopocentricFrame(earthSpheric, point, "lon 0 lat 90");
655 topoElliptic = new TopocentricFrame(earthElliptic, point, "lon 0 lat 90");
656 tcElliptic = topoElliptic.getTrackingCoordinates(satPoint, earthElliptic.getBodyFrame(), fieldDate);
657 tcSpheric = topoSpheric.getTrackingCoordinates(satPoint, earthSpheric.getBodyFrame(), fieldDate);
658
659
660
661 Assertions.assertEquals(tcElliptic.getAzimuth().getReal(), tcSpheric.getAzimuth().getReal(), 1.e-7);
662 Assertions.assertEquals(tcElliptic.getElevation().getReal(), tcSpheric.getElevation().getReal(), 1.e-2);
663 Assertions.assertEquals(tcElliptic.getRange().getReal(), tcSpheric.getRange().getReal(), 20.e+3);
664
665
666
667
668
669 point = new GeodeticPoint(FastMath.toRadians(60), FastMath.toRadians(30.), 0.);
670 topoSpheric = new TopocentricFrame(earthSpheric, point, "lon 10 lat 45");
671 topoElliptic = new TopocentricFrame(earthElliptic, point, "lon 10 lat 45");
672 tcElliptic = topoElliptic.getTrackingCoordinates(satPoint, earthElliptic.getBodyFrame(), fieldDate);
673 tcSpheric = topoSpheric.getTrackingCoordinates(satPoint, earthSpheric.getBodyFrame(), fieldDate);
674
675
676
677 Assertions.assertEquals(tcElliptic.getAzimuth().getReal(), tcSpheric.getAzimuth().getReal(), 1.e-2);
678 Assertions.assertEquals(tcElliptic.getElevation().getReal(), tcSpheric.getElevation().getReal(), 1.e-2);
679 Assertions.assertEquals(tcElliptic.getRange().getReal(), tcSpheric.getRange().getReal(), 20.e+3);
680
681 }
682
683 @Test
684 void testPointAtDistance() {
685
686 RandomGenerator random = new Well1024a(0xa1e6bd5cd0578779l);
687 final OneAxisEllipsoid earth =
688 new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING,
689 itrf);
690 final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
691
692 for (int i = 0; i < 20; ++i) {
693
694 double latitude = FastMath.PI * (0.5 - random.nextDouble());
695 double longitude = 2 * FastMath.PI * random.nextDouble();
696 TopocentricFrame topo = new TopocentricFrame(earth,
697 new GeodeticPoint(latitude, longitude, 0.0),
698 "topo");
699 Transform transform = earth.getBodyFrame().getTransformTo(topo, date);
700 for (int j = 0; j < 20; ++j) {
701 double elevation = FastMath.PI * (0.5 - random.nextDouble());
702 double azimuth = 2 * FastMath.PI * random.nextDouble();
703 double range = 500000.0 * (1.0 + random.nextDouble());
704 Vector3D absolutePoint = earth.transform(topo.pointAtDistance(azimuth, elevation, range));
705 Vector3D relativePoint = transform.transformPosition(absolutePoint);
706 TrackingCoordinates rebuiltTC = topo.getTrackingCoordinates(relativePoint, topo, AbsoluteDate.J2000_EPOCH);
707 Assertions.assertEquals(elevation, rebuiltTC.getElevation(), 1.0e-12);
708 Assertions.assertEquals(azimuth, MathUtils.normalizeAngle(rebuiltTC.getAzimuth(), azimuth), 1.0e-12);
709 Assertions.assertEquals(range, rebuiltTC.getRange(), 1.0e-12 * range);
710 }
711 }
712 }
713
714 @Test
715 void testIssue145() {
716 Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
717 BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
718 Constants.WGS84_EARTH_FLATTENING,
719 itrf);
720 TopocentricFrame staFrame = new TopocentricFrame(earth, new GeodeticPoint(0.0, 0.0, 0.0), "test");
721 GeodeticPoint gp = staFrame.computeLimitVisibilityPoint(Constants.WGS84_EARTH_EQUATORIAL_RADIUS+600000,
722 0.0, FastMath.toRadians(5.0));
723 Assertions.assertEquals(0.0, gp.getLongitude(), 1.0e-15);
724 Assertions.assertTrue(gp.getLatitude() > 0);
725 Assertions.assertEquals(0.0, staFrame.getNorth().distance(Vector3D.PLUS_K), 1.0e-15);
726
727 }
728
729 @Test
730 void testVisibilityCircle() throws IOException {
731
732
733 final BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
734 Constants.WGS84_EARTH_FLATTENING,
735 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
736 final TopocentricFrame[] ilrs = {
737 new TopocentricFrame(earth,
738 new GeodeticPoint(FastMath.toRadians(52.3800), FastMath.toRadians(3.0649), 133.745),
739 "Potsdam"),
740 new TopocentricFrame(earth,
741 new GeodeticPoint(FastMath.toRadians(36.46273), FastMath.toRadians(-6.20619), 64.0),
742 "San Fernando"),
743 new TopocentricFrame(earth,
744 new GeodeticPoint(FastMath.toRadians(35.5331), FastMath.toRadians(24.0705), 157.0),
745 "Chania")
746 };
747
748 PolynomialFunction distanceModel =
749 new PolynomialFunction(new double[] { 7.0892e+05, 3.1913, -8.2181e-07, 1.4033e-13 });
750 for (TopocentricFrame station : ilrs) {
751 for (double altitude = 500000; altitude < 2000000; altitude += 100000) {
752 for (double azimuth = 0; azimuth < 2 * FastMath.PI; azimuth += 0.05) {
753 GeodeticPoint p = station.computeLimitVisibilityPoint(Constants.WGS84_EARTH_EQUATORIAL_RADIUS + altitude,
754 azimuth, FastMath.toRadians(5.0));
755 double d = station.getTrackingCoordinates(earth.transform(p), earth.getBodyFrame(), AbsoluteDate.J2000_EPOCH).getRange();
756 Assertions.assertEquals(distanceModel.value(altitude), d, 40000.0);
757 }
758 }
759 }
760
761 }
762
763 private static Stream<Arguments> testGetTopocentricCoordinatesValues() {
764 return Stream.of(
765 Arguments.of(0,0,1),
766 Arguments.of(0, FastMath.PI / 2, 1),
767 Arguments.of(FastMath.PI, FastMath.PI / 3, 10),
768 Arguments.of(3 * FastMath.PI / 2, FastMath.PI / 4, 1000),
769 Arguments.of(FastMath.PI / 2, -FastMath.PI / 6, 500),
770 Arguments.of(FastMath.PI / 7, -FastMath.PI / 5, 100)
771 );
772 }
773
774 @ParameterizedTest
775 @MethodSource("testGetTopocentricCoordinatesValues")
776 void testGetTopocentricCoordinates(final double az, final double el, final double r) {
777 final TrackingCoordinates coords = new TrackingCoordinates(az,el,r);
778 final Vector3D topoPos = TopocentricFrame.getTopocentricPosition(coords);
779
780 final SinCos sinCosAz = FastMath.sinCos(az);
781 final SinCos sinCosEl = FastMath.sinCos(el);
782 final double expectedX, expectedY, expectedZ;
783 expectedX = sinCosAz.sin() * sinCosEl.cos() * r;
784 expectedY = sinCosAz.cos() * sinCosEl.cos() * r;
785 expectedZ = sinCosEl.sin() * r;
786
787 final double distanceTolerance = 1e-7;
788 Assertions.assertEquals(expectedX,topoPos.getX(), distanceTolerance);
789 Assertions.assertEquals(expectedY,topoPos.getY(), distanceTolerance);
790 Assertions.assertEquals(expectedZ,topoPos.getZ(), distanceTolerance);
791 }
792
793 @ParameterizedTest
794 @MethodSource("testGetTopocentricCoordinatesValues")
795 void testGetFieldTopocentricCoordinates(final double azimuth, final double elevation, final double range) {
796 final Binary64 az, el, r;
797 az = new Binary64(azimuth);
798 el = new Binary64(elevation);
799 r = new Binary64(range);
800
801 final FieldTrackingCoordinates<Binary64> coords = new FieldTrackingCoordinates<>(az, el, r);
802 final FieldVector3D<Binary64> topoPos = TopocentricFrame.getTopocentricPosition(coords);
803
804 final FieldSinCos<Binary64> sinCosAz = FastMath.sinCos(az);
805 final FieldSinCos<Binary64> sinCosEl = FastMath.sinCos(el);
806 final Binary64 expectedX, expectedY, expectedZ;
807 expectedX = r.multiply(sinCosAz.sin().multiply(sinCosEl.cos()));
808 expectedY = r.multiply(sinCosEl.cos().multiply(sinCosAz.cos()));
809 expectedZ = r.multiply(sinCosEl.sin());
810
811 final double distanceTolerance = 1e-7;
812 Assertions.assertEquals(expectedX.getReal(), topoPos.getX().getReal(), distanceTolerance);
813 Assertions.assertEquals(expectedY.getReal(), topoPos.getY().getReal(), distanceTolerance);
814 Assertions.assertEquals(expectedZ.getReal(), topoPos.getZ().getReal(), distanceTolerance);
815 }
816
817 @ParameterizedTest
818 @MethodSource("testGetTopocentricCoordinatesValues")
819 void testInverseGetTopocentricCoordinates(double az, double el, double r) {
820 final TrackingCoordinates expectedCoords = new TrackingCoordinates(az, el, r);
821
822 final Vector3D point = TopocentricFrame.getTopocentricPosition(expectedCoords);
823 final GeodeticPoint geodeticPoint = new GeodeticPoint(0., 0., 0.);
824 final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, geodeticPoint, "geodeticPoint");
825 final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
826 final TrackingCoordinates coords = topoFrame.getTrackingCoordinates(point, topoFrame, date);
827
828 final double angularTolerance = 1e-12;
829 Assertions.assertEquals(expectedCoords.getAzimuth(), coords.getAzimuth(), angularTolerance);
830 Assertions.assertEquals(expectedCoords.getElevation(), coords.getElevation(), angularTolerance);
831 Assertions.assertEquals(expectedCoords.getRange(), coords.getRange());
832 }
833
834 @Test
835 void testGetTrackingCoordinates() {
836
837 final GeodeticPoint geodeticPoint = new GeodeticPoint(0., 0., 0.);
838 final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, geodeticPoint, "geodeticPoint");
839 final AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
840 final Frame frame = FramesFactory.getGCRF();
841 final Vector3D point = new Vector3D(1., 2., -3.);
842
843 final TrackingCoordinates trackingCoordinates = topoFrame.getTrackingCoordinates(point, frame, date);
844
845 final double expectedElevation = topoFrame.getElevation(point, frame, date);
846 final double expectedAzimuth = topoFrame.getAzimuth(point, frame, date);
847 final double expectedRange = topoFrame.getRange(point, frame, date);
848 Assertions.assertEquals(expectedElevation, trackingCoordinates.getElevation());
849 Assertions.assertEquals(expectedAzimuth, trackingCoordinates.getAzimuth());
850 Assertions.assertEquals(expectedRange, trackingCoordinates.getRange());
851 }
852
853 @Test
854 void testGetFieldTrackingCoordinates() {
855
856 final GeodeticPoint geodeticPoint = new GeodeticPoint(0., 0., 0.);
857 final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, geodeticPoint, "geodeticPoint");
858 final ComplexField field = ComplexField.getInstance();
859 final FieldAbsoluteDate<Complex> date = FieldAbsoluteDate.getArbitraryEpoch(field);
860 final Frame frame = FramesFactory.getGCRF();
861 final FieldVector3D<Complex> point = new FieldVector3D<>(field, new Vector3D(1., 2., -3.));
862
863 final FieldTrackingCoordinates<Complex> trackingCoordinates = topoFrame.getTrackingCoordinates(point, frame, date);
864
865 final Complex expectedElevation = topoFrame.getElevation(point, frame, date);
866 final Complex expectedAzimuth = topoFrame.getAzimuth(point, frame, date);
867 final Complex expectedRange = topoFrame.getRange(point, frame, date);
868 Assertions.assertEquals(expectedElevation, trackingCoordinates.getElevation());
869 Assertions.assertEquals(expectedAzimuth, trackingCoordinates.getAzimuth());
870 Assertions.assertEquals(expectedRange, trackingCoordinates.getRange());
871 }
872
873 @BeforeEach
874 public void setUp() {
875
876 Utils.setDataRoot("regular-data");
877
878
879 itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
880
881
882 earthSpheric = new OneAxisEllipsoid(6378136.460, 0., itrf);
883
884
885 date = new AbsoluteDate(new DateComponents(2008, 04, 07),
886 TimeComponents.H00,
887 TimeScalesFactory.getUTC());
888
889
890 mu = 3.9860047e14;
891
892 }
893
894 @AfterEach
895 public void tearDown() {
896 date = null;
897 itrf = null;
898 earthSpheric = null;
899 }
900
901
902 }