1   /* Copyright 2002-2022 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.frames;
18  
19  
20  
21  import java.io.IOException;
22  
23  import org.hipparchus.Field;
24  import org.hipparchus.CalculusFieldElement;
25  import org.hipparchus.analysis.polynomials.PolynomialFunction;
26  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
27  import org.hipparchus.geometry.euclidean.threed.Vector3D;
28  import org.hipparchus.random.RandomGenerator;
29  import org.hipparchus.random.Well1024a;
30  import org.hipparchus.util.Decimal64Field;
31  import org.hipparchus.util.FastMath;
32  import org.hipparchus.util.MathUtils;
33  import org.junit.After;
34  import org.junit.Assert;
35  import org.junit.Before;
36  import org.junit.Test;
37  import org.orekit.Utils;
38  import org.orekit.bodies.BodyShape;
39  import org.orekit.bodies.FieldGeodeticPoint;
40  import org.orekit.bodies.GeodeticPoint;
41  import org.orekit.bodies.OneAxisEllipsoid;
42  import org.orekit.errors.OrekitException;
43  import org.orekit.orbits.CircularOrbit;
44  import org.orekit.orbits.FieldCircularOrbit;
45  import org.orekit.orbits.PositionAngle;
46  import org.orekit.propagation.FieldSpacecraftState;
47  import org.orekit.propagation.SpacecraftState;
48  import org.orekit.propagation.analytical.FieldKeplerianPropagator;
49  import org.orekit.propagation.analytical.KeplerianPropagator;
50  import org.orekit.time.AbsoluteDate;
51  import org.orekit.time.DateComponents;
52  import org.orekit.time.FieldAbsoluteDate;
53  import org.orekit.time.TimeComponents;
54  import org.orekit.time.TimeScalesFactory;
55  import org.orekit.utils.Constants;
56  import org.orekit.utils.FieldPVCoordinates;
57  import org.orekit.utils.IERSConventions;
58  import org.orekit.utils.PVCoordinates;
59  
60  
61  public class TopocentricFrameTest {
62  
63      // Computation date
64      private AbsoluteDate date;
65  
66      // Reference frame = ITRF
67      private Frame itrf;
68  
69      // Earth shape
70      OneAxisEllipsoid earthSpheric;
71  
72      // Body mu
73      private double mu;
74  
75  
76      @Test
77      public void testZero() {
78  
79          final GeodeticPoint point = new GeodeticPoint(0., 0., 0.);
80          final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "zero");
81  
82          // Check that frame directions are aligned
83          final double xDiff = Vector3D.dotProduct(topoFrame.getEast(), Vector3D.PLUS_J);
84          final double yDiff = Vector3D.dotProduct(topoFrame.getNorth(), Vector3D.PLUS_K);
85          final double zDiff = Vector3D.dotProduct(topoFrame.getZenith(), Vector3D.PLUS_I);
86          Assert.assertEquals(1., xDiff, Utils.epsilonTest);
87          Assert.assertEquals(1., yDiff, Utils.epsilonTest);
88          Assert.assertEquals(1., zDiff, Utils.epsilonTest);
89     }
90  
91      @Test
92      public void testPole() {
93  
94          final GeodeticPoint point = new GeodeticPoint(FastMath.PI/2., 0., 0.);
95          final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "north pole");
96  
97          // Check that frame directions are aligned
98          final double xDiff = Vector3D.dotProduct(topoFrame.getEast(), Vector3D.PLUS_J);
99          final double yDiff = Vector3D.dotProduct(topoFrame.getSouth(), Vector3D.PLUS_I);
100         final double zDiff = Vector3D.dotProduct(topoFrame.getZenith(), Vector3D.PLUS_K);
101         Assert.assertEquals(1., xDiff, Utils.epsilonTest);
102         Assert.assertEquals(1., yDiff, Utils.epsilonTest);
103         Assert.assertEquals(1., zDiff, Utils.epsilonTest);
104    }
105 
106     @Test
107     public void testNormalLatitudes() {
108 
109         // First point at latitude 45°
110         final GeodeticPoint point1 = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(30.), 0.);
111         final TopocentricFrame topoFrame1 = new TopocentricFrame(earthSpheric, point1, "lat 45");
112 
113         // Second point at latitude -45° and same longitude
114         final GeodeticPoint point2 = new GeodeticPoint(FastMath.toRadians(-45.), FastMath.toRadians(30.), 0.);
115         final TopocentricFrame topoFrame2 = new TopocentricFrame(earthSpheric, point2, "lat -45");
116 
117         // Check that frame North and Zenith directions are all normal to each other, and East are the same
118         final double xDiff = Vector3D.dotProduct(topoFrame1.getEast(), topoFrame2.getEast());
119         final double yDiff = Vector3D.dotProduct(topoFrame1.getNorth(), topoFrame2.getNorth());
120         final double zDiff = Vector3D.dotProduct(topoFrame1.getZenith(), topoFrame2.getZenith());
121 
122         Assert.assertEquals(1., xDiff, Utils.epsilonTest);
123         Assert.assertEquals(0., yDiff, Utils.epsilonTest);
124         Assert.assertEquals(0., zDiff, Utils.epsilonTest);
125   }
126 
127     @Test
128     public void testOppositeLongitudes() {
129 
130         // First point at latitude 45°
131         final GeodeticPoint point1 = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(30.), 0.);
132         final TopocentricFrame topoFrame1 = new TopocentricFrame(earthSpheric, point1, "lon 30");
133         final GeodeticPoint p1 = topoFrame1.getPoint();
134         Assert.assertEquals(point1.getLatitude(), p1.getLatitude(), 1.0e-15);
135         Assert.assertEquals(point1.getLongitude(), p1.getLongitude(), 1.0e-15);
136         Assert.assertEquals(point1.getAltitude(), p1.getAltitude(), 1.0e-15);
137 
138         // Second point at latitude -45° and same longitude
139         final GeodeticPoint point2 = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(210.), 0.);
140         final TopocentricFrame topoFrame2 = new TopocentricFrame(earthSpheric, point2, "lon 210");
141 
142         // Check that frame North and Zenith directions are all normal to each other,
143         // and East of the one is West of the other
144         final double xDiff = Vector3D.dotProduct(topoFrame1.getEast(), topoFrame2.getWest());
145         final double yDiff = Vector3D.dotProduct(topoFrame1.getNorth(), topoFrame2.getNorth());
146         final double zDiff = Vector3D.dotProduct(topoFrame1.getZenith(), topoFrame2.getZenith());
147 
148         Assert.assertEquals(1., xDiff, Utils.epsilonTest);
149         Assert.assertEquals(0., yDiff, Utils.epsilonTest);
150         Assert.assertEquals(0., zDiff, Utils.epsilonTest);
151   }
152 
153     @Test
154     public void testAntipodes()
155         {
156 
157         // First point at latitude 45° and longitude 30
158         final GeodeticPoint point1 = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(30.), 0.);
159         final TopocentricFrame topoFrame1 = new TopocentricFrame(earthSpheric, point1, "lon 30");
160 
161         // Second point at latitude -45° and longitude 210
162         final GeodeticPoint point2 = new GeodeticPoint(FastMath.toRadians(-45.), FastMath.toRadians(210.), 0.);
163         final TopocentricFrame topoFrame2 = new TopocentricFrame(earthSpheric, point2, "lon 210");
164 
165         // Check that frame Zenith directions are opposite to each other,
166         // and East and North are the same
167         final double xDiff = Vector3D.dotProduct(topoFrame1.getEast(), topoFrame2.getWest());
168         final double yDiff = Vector3D.dotProduct(topoFrame1.getNorth(), topoFrame2.getNorth());
169         final double zDiff = Vector3D.dotProduct(topoFrame1.getZenith(), topoFrame2.getZenith());
170 
171         Assert.assertEquals(1., xDiff, Utils.epsilonTest);
172         Assert.assertEquals(1., yDiff, Utils.epsilonTest);
173         Assert.assertEquals(-1., zDiff, Utils.epsilonTest);
174 
175         Assert.assertEquals(1, Vector3D.dotProduct(topoFrame1.getNadir(), topoFrame2.getZenith()), Utils.epsilonTest);
176         Assert.assertEquals(1, Vector3D.dotProduct(topoFrame1.getZenith(), topoFrame2.getNadir()), Utils.epsilonTest);
177 
178     }
179 
180     @Test
181     public void testSiteAtZenith()
182         {
183 
184         // Surface point at latitude 45°
185         final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(30.), 0.);
186         final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 30 lat 45");
187 
188         // Point at 800 km over zenith
189         final GeodeticPoint satPoint = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(30.), 800000.);
190 
191         // Zenith point elevation = 90 deg
192         final double site = topoFrame.getElevation(earthSpheric.transform(satPoint), earthSpheric.getBodyFrame(), date);
193         Assert.assertEquals(FastMath.PI/2., site, Utils.epsilonAngle);
194 
195         // Zenith point range = defined altitude
196         final double range = topoFrame.getRange(earthSpheric.transform(satPoint), earthSpheric.getBodyFrame(), date);
197         Assert.assertEquals(800000., range, 1e-8);
198   }
199 
200     @Test
201     public void testFieldSiteAtZenith() {
202         doTestFieldSiteAtZenith(Decimal64Field.getInstance());
203     }
204 
205     private <T extends CalculusFieldElement<T>> void doTestFieldSiteAtZenith(final Field<T> field)
206         {
207 
208         // zero
209         final T zero = field.getZero();
210 
211         // Surface point at latitude 45°
212         final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(30.), 0.);
213         final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 30 lat 45");
214 
215         // Point at 800 km over zenith
216         final FieldGeodeticPoint<T> satPoint = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(45.)),
217                                                                         zero.add(FastMath.toRadians(30.)),
218                                                                         zero.add(800000.));
219 
220         // Field date
221         final FieldAbsoluteDate<T> fieldDate = new FieldAbsoluteDate<>(field, date);
222         // Zenith point elevation = 90 deg
223         final T site = topoFrame.getElevation(earthSpheric.transform(satPoint), earthSpheric.getBodyFrame(), fieldDate);
224         Assert.assertEquals(FastMath.PI/2., site.getReal(), Utils.epsilonAngle);
225 
226         // Zenith point range = defined altitude
227         final T range = topoFrame.getRange(earthSpheric.transform(satPoint), earthSpheric.getBodyFrame(), fieldDate);
228         Assert.assertEquals(800000., range.getReal(), 1e-8);
229   }
230 
231     @Test
232     public void testAzimuthEquatorial()
233         {
234 
235         // Surface point at latitude 0
236         final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(30.), 0.);
237         final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 30 lat 0");
238 
239         // Point at infinite, separated by +20 deg in longitude
240         // *****************************************************
241         GeodeticPoint infPoint = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(50.), 1000000000.);
242 
243         // Azimuth = pi/2
244         double azi = topoFrame.getAzimuth(earthSpheric.transform(infPoint), earthSpheric.getBodyFrame(), date);
245         Assert.assertEquals(FastMath.PI/2., azi, Utils.epsilonAngle);
246 
247         // Site = pi/2 - longitude difference
248         double site = topoFrame.getElevation(earthSpheric.transform(infPoint), earthSpheric.getBodyFrame(), date);
249         Assert.assertEquals(FastMath.PI/2. - FastMath.abs(point.getLongitude() - infPoint.getLongitude()), site, 1.e-2);
250 
251         // Point at infinite, separated by -20 deg in longitude
252         // *****************************************************
253         infPoint = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(10.), 1000000000.);
254 
255         // Azimuth = pi/2
256         azi = topoFrame.getAzimuth(earthSpheric.transform(infPoint), earthSpheric.getBodyFrame(), date);
257         Assert.assertEquals(3*FastMath.PI/2., azi, Utils.epsilonAngle);
258 
259         // Site = pi/2 - longitude difference
260         site = topoFrame.getElevation(earthSpheric.transform(infPoint), earthSpheric.getBodyFrame(), date);
261         Assert.assertEquals(FastMath.PI/2. - FastMath.abs(point.getLongitude() - infPoint.getLongitude()), site, 1.e-2);
262 
263     }
264 
265     @Test
266     public void testFieldAzimuthEquatorial() {
267         doTestFieldAzimuthEquatorial(Decimal64Field.getInstance());
268     }
269 
270     private <T extends CalculusFieldElement<T>> void doTestFieldAzimuthEquatorial(final Field<T> field)
271         {
272 
273         // zero
274         final T zero = field.getZero();
275 
276         // Surface point at latitude 0
277         final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(30.), 0.);
278         final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 30 lat 0");
279 
280         // Point at infinite, separated by +20 deg in longitude
281         // *****************************************************
282         FieldGeodeticPoint<T> infPoint = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(0.)),
283                                                                   zero.add(FastMath.toRadians(50.)),
284                                                                   zero.add(1000000000.));
285 
286         // Field date
287         final FieldAbsoluteDate<T> fieldDate = new FieldAbsoluteDate<>(field, date);
288         // Azimuth = pi/2
289         T azi = topoFrame.getAzimuth(earthSpheric.transform(infPoint), earthSpheric.getBodyFrame(), fieldDate);
290         Assert.assertEquals(FastMath.PI/2., azi.getReal(), Utils.epsilonAngle);
291 
292         // Site = pi/2 - longitude difference
293         T site = topoFrame.getElevation(earthSpheric.transform(infPoint), earthSpheric.getBodyFrame(), fieldDate);
294         Assert.assertEquals(FastMath.abs(infPoint.getLongitude().negate().add(point.getLongitude())).negate().add(FastMath.PI/2.).getReal(), site.getReal(), 1.e-2);
295 
296         // Point at infinite, separated by -20 deg in longitude
297         // *****************************************************
298         infPoint = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(0.)),
299                                             zero.add(FastMath.toRadians(10.)),
300                                             zero.add(1000000000.));
301 
302         // Azimuth = pi/2
303         azi = topoFrame.getAzimuth(earthSpheric.transform(infPoint), earthSpheric.getBodyFrame(), fieldDate);
304         Assert.assertEquals(3*FastMath.PI/2., azi.getReal(), Utils.epsilonAngle);
305 
306         // Site = pi/2 - longitude difference
307         site = topoFrame.getElevation(earthSpheric.transform(infPoint), earthSpheric.getBodyFrame(), fieldDate);
308         Assert.assertEquals(FastMath.abs(infPoint.getLongitude().negate().add(point.getLongitude())).negate().add(FastMath.PI/2.).getReal(), site.getReal(), 1.e-2);
309 
310     }
311 
312     @Test
313     public void testAzimuthPole()
314         {
315 
316         // Surface point at latitude 0
317         final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(89.999), FastMath.toRadians(0.), 0.);
318         final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 0 lat 90");
319 
320         // Point at 30 deg longitude
321         // **************************
322         GeodeticPoint satPoint = new GeodeticPoint(FastMath.toRadians(28.), FastMath.toRadians(30.), 800000.);
323 
324         // Azimuth =
325         double azi = topoFrame.getAzimuth(earthSpheric.transform(satPoint), earthSpheric.getBodyFrame(), date);
326         Assert.assertEquals(FastMath.PI - satPoint.getLongitude(), azi, 1.e-5);
327 
328         // Point at -30 deg longitude
329         // ***************************
330         satPoint = new GeodeticPoint(FastMath.toRadians(28.), FastMath.toRadians(-30.), 800000.);
331 
332         // Azimuth =
333         azi = topoFrame.getAzimuth(earthSpheric.transform(satPoint), earthSpheric.getBodyFrame(), date);
334         Assert.assertEquals(FastMath.PI - satPoint.getLongitude(), azi, 1.e-5);
335 
336     }
337 
338     @Test
339     public void testFieldAzimuthPole() {
340         doTestFieldAzimuthPole(Decimal64Field.getInstance());
341     }
342 
343     private <T extends CalculusFieldElement<T>> void doTestFieldAzimuthPole(final Field<T> field)
344         {
345 
346         // zero 
347         final T zero = field.getZero();
348 
349         // Surface point at latitude 0
350         final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(89.999), FastMath.toRadians(0.), 0.);
351         final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 0 lat 90");
352 
353         // Point at 30 deg longitude
354         // **************************
355         FieldGeodeticPoint<T> satPoint = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(28.)),
356                                                                   zero.add(FastMath.toRadians(30.)), zero.add(800000.));
357 
358         // Field date
359         final FieldAbsoluteDate<T> fieldDate = new FieldAbsoluteDate<>(field, date);
360 
361         // Azimuth =
362         T azi = topoFrame.getAzimuth(earthSpheric.transform(satPoint), earthSpheric.getBodyFrame(), fieldDate);
363         Assert.assertEquals(satPoint.getLongitude().negate().add(FastMath.PI).getReal(), azi.getReal(), 1.e-5);
364 
365         // Point at -30 deg longitude
366         // ***************************
367         satPoint = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(28.)),
368                                             zero.add(FastMath.toRadians(-30.)), zero.add(800000.));
369 
370         // Azimuth =
371         azi = topoFrame.getAzimuth(earthSpheric.transform(satPoint), earthSpheric.getBodyFrame(), fieldDate);
372         Assert.assertEquals(satPoint.getLongitude().negate().add(FastMath.PI).getReal(), azi.getReal(), 1.e-5);
373 
374     }
375 
376     @Test
377     public void testDoppler()
378         {
379 
380         // Surface point at latitude 45, longitude 5
381         final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(5.), 0.);
382         final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 5 lat 45");
383 
384         // Point at 30 deg longitude
385         // ***************************
386         final CircularOrbit orbit =
387             new CircularOrbit(7178000.0, 0.5e-8, -0.5e-8, FastMath.toRadians(50.), FastMath.toRadians(120.),
388                                    FastMath.toRadians(90.), PositionAngle.MEAN,
389                                    FramesFactory.getEME2000(), date, mu);
390 
391         // Transform satellite position to position/velocity parameters in body frame
392         final Transform eme2000ToItrf = FramesFactory.getEME2000().getTransformTo(earthSpheric.getBodyFrame(), date);
393         final PVCoordinates pvSatItrf = eme2000ToItrf.transformPVCoordinates(orbit.getPVCoordinates());
394 
395         // Compute range rate directly
396         //********************************************
397         final double dop = topoFrame.getRangeRate(pvSatItrf, earthSpheric.getBodyFrame(), date);
398 
399         // Compare to finite difference computation (2 points)
400         //*****************************************************
401         final double dt = 0.1;
402         KeplerianPropagator extrapolator = new KeplerianPropagator(orbit);
403 
404         // Extrapolate satellite position a short while after reference date
405         AbsoluteDate dateP = date.shiftedBy(dt);
406         Transform j2000ToItrfP = FramesFactory.getEME2000().getTransformTo(earthSpheric.getBodyFrame(), dateP);
407         SpacecraftState orbitP = extrapolator.propagate(dateP);
408         Vector3D satPointGeoP = j2000ToItrfP.transformPVCoordinates(orbitP.getPVCoordinates()).getPosition();
409 
410         // Retropolate satellite position a short while before reference date
411         AbsoluteDate dateM = date.shiftedBy(-dt);
412         Transform j2000ToItrfM = FramesFactory.getEME2000().getTransformTo(earthSpheric.getBodyFrame(), dateM);
413         SpacecraftState orbitM = extrapolator.propagate(dateM);
414         Vector3D satPointGeoM = j2000ToItrfM.transformPVCoordinates(orbitM.getPVCoordinates()).getPosition();
415 
416         // Compute ranges at both instants
417         double rangeP = topoFrame.getRange(satPointGeoP, earthSpheric.getBodyFrame(), dateP);
418         double rangeM = topoFrame.getRange(satPointGeoM, earthSpheric.getBodyFrame(), dateM);
419         final double dopRef2 = (rangeP - rangeM) / (2. * dt);
420         Assert.assertEquals(dopRef2, dop, 1.e-3);
421 
422     }
423 
424     @Test
425     public void testFieldDoppler() {
426         doTestFieldDoppler(Decimal64Field.getInstance());
427     }
428 
429     private <T extends CalculusFieldElement<T>> void doTestFieldDoppler(final Field<T> field)
430         {
431 
432         // zero
433         final T zero = field.getZero();
434 
435         // Surface point at latitude 45, longitude 5
436         final GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(45.), FastMath.toRadians(5.), 0.);
437         final TopocentricFrame topoFrame = new TopocentricFrame(earthSpheric, point, "lon 5 lat 45");
438 
439         // Field date
440         final FieldAbsoluteDate<T> fieldDate = new FieldAbsoluteDate<>(field, date);
441 
442         // Point at 30 deg longitude
443         // ***************************
444         final FieldCircularOrbit<T> orbit =
445             new FieldCircularOrbit<>(zero.add(7178000.0), zero.add(0.5e-8), zero.add(-0.5e-8),
446                                      zero.add(FastMath.toRadians(50.)), zero.add(FastMath.toRadians(120.)),
447                                      zero.add(FastMath.toRadians(90.)), PositionAngle.MEAN,
448                                      FramesFactory.getEME2000(), fieldDate, zero.add(mu));
449 
450         // Transform satellite position to position/velocity parameters in body frame
451         final FieldTransform<T> eme2000ToItrf = FramesFactory.getEME2000().getTransformTo(earthSpheric.getBodyFrame(), fieldDate);
452         final FieldPVCoordinates<T> pvSatItrf = eme2000ToItrf.transformPVCoordinates(orbit.getPVCoordinates());
453 
454         // Compute range rate directly
455         //********************************************
456         final T dop = topoFrame.getRangeRate(pvSatItrf, earthSpheric.getBodyFrame(), fieldDate);
457 
458         // Compare to finite difference computation (2 points)
459         //*****************************************************
460         final double dt = 0.1;
461         FieldKeplerianPropagator<T> extrapolator = new FieldKeplerianPropagator<>(orbit);
462 
463         // Extrapolate satellite position a short while after reference date
464         FieldAbsoluteDate<T> dateP = fieldDate.shiftedBy(dt);
465         FieldTransform<T> j2000ToItrfP = FramesFactory.getEME2000().getTransformTo(earthSpheric.getBodyFrame(), dateP);
466         FieldSpacecraftState<T> orbitP = extrapolator.propagate(dateP);
467         FieldVector3D<T> satPointGeoP = j2000ToItrfP.transformPVCoordinates(orbitP.getPVCoordinates()).getPosition();
468 
469         // Retropolate satellite position a short while before reference date
470         FieldAbsoluteDate<T> dateM = fieldDate.shiftedBy(-dt);
471         FieldTransform<T> j2000ToItrfM = FramesFactory.getEME2000().getTransformTo(earthSpheric.getBodyFrame(), dateM);
472         FieldSpacecraftState<T> orbitM = extrapolator.propagate(dateM);
473         FieldVector3D<T> satPointGeoM = j2000ToItrfM.transformPVCoordinates(orbitM.getPVCoordinates()).getPosition();
474 
475         // Compute ranges at both instants
476         T rangeP = topoFrame.getRange(satPointGeoP, earthSpheric.getBodyFrame(), dateP);
477         T rangeM = topoFrame.getRange(satPointGeoM, earthSpheric.getBodyFrame(), dateM);
478         final T dopRef2 = (rangeP.subtract(rangeM)).divide(2. * dt);
479         Assert.assertEquals(dopRef2.getReal(), dop.getReal(), 1.e-3);
480 
481     }
482 
483     @Test
484     public void testEllipticEarth()  {
485 
486         // Elliptic earth shape
487         final OneAxisEllipsoid earthElliptic =
488             new OneAxisEllipsoid(6378136.460, 1 / 298.257222101, itrf);
489 
490         // Satellite point
491         // Caution !!! Sat point target shall be the same whatever earth shape chosen !!
492         final GeodeticPoint satPointGeo = new GeodeticPoint(FastMath.toRadians(30.), FastMath.toRadians(15.), 800000.);
493         final Vector3D satPoint = earthElliptic.transform(satPointGeo);
494 
495         // ****************************
496         // Test at equatorial position
497         // ****************************
498         GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(5.), 0.);
499         TopocentricFrame topoElliptic = new TopocentricFrame(earthElliptic, point, "elliptic, equatorial lon 5");
500         TopocentricFrame topoSpheric = new TopocentricFrame(earthSpheric, point, "spheric, equatorial lon 5");
501 
502         // Compare azimuth/elevation/range of satellite point : shall be strictly identical
503         // ***************************************************
504         double aziElli = topoElliptic.getAzimuth(satPoint, earthElliptic.getBodyFrame(), date);
505         double aziSphe = topoSpheric.getAzimuth(satPoint, earthSpheric.getBodyFrame(), date);
506         Assert.assertEquals(aziElli, aziSphe, Utils.epsilonAngle);
507 
508         double eleElli = topoElliptic.getElevation(satPoint, earthElliptic.getBodyFrame(), date);
509         double eleSphe = topoSpheric.getElevation(satPoint, earthSpheric.getBodyFrame(), date);
510         Assert.assertEquals(eleElli, eleSphe, Utils.epsilonAngle);
511 
512         double disElli = topoElliptic.getRange(satPoint, earthElliptic.getBodyFrame(), date);
513         double disSphe = topoSpheric.getRange(satPoint, earthSpheric.getBodyFrame(), date);
514         Assert.assertEquals(disElli, disSphe, Utils.epsilonTest);
515 
516         // Infinite point separated by -20 deg in longitude
517         // *************************************************
518         GeodeticPoint infPointGeo = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(-15.), 1000000000.);
519         Vector3D infPoint = earthElliptic.transform(infPointGeo);
520 
521         // Azimuth = pi/2
522         aziElli = topoElliptic.getAzimuth(infPoint, earthElliptic.getBodyFrame(), date);
523         Assert.assertEquals(3*FastMath.PI/2., aziElli, Utils.epsilonAngle);
524 
525         // Site = pi/2 - longitude difference
526         eleElli = topoElliptic.getElevation(infPoint, earthElliptic.getBodyFrame(), date);
527         Assert.assertEquals(FastMath.PI/2. - FastMath.abs(point.getLongitude() - infPointGeo.getLongitude()), eleElli, 1.e-2);
528 
529         // Infinite point separated by +20 deg in longitude
530         // *************************************************
531         infPointGeo = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(25.), 1000000000.);
532         infPoint = earthElliptic.transform(infPointGeo);
533 
534         // Azimuth = pi/2
535         aziElli = topoElliptic.getAzimuth(infPoint, earthElliptic.getBodyFrame(), date);
536         Assert.assertEquals(FastMath.PI/2., aziElli, Utils.epsilonAngle);
537 
538         // Site = pi/2 - longitude difference
539         eleElli = topoElliptic.getElevation(infPoint, earthElliptic.getBodyFrame(), date);
540         Assert.assertEquals(FastMath.PI/2. - FastMath.abs(point.getLongitude() - infPointGeo.getLongitude()), eleElli, 1.e-2);
541 
542         // ************************
543         // Test at polar position
544         // ************************
545         point = new GeodeticPoint(FastMath.toRadians(89.999), FastMath.toRadians(0.), 0.);
546         topoSpheric  = new TopocentricFrame(earthSpheric, point, "lon 0 lat 90");
547         topoElliptic = new TopocentricFrame(earthElliptic, point, "lon 0 lat 90");
548 
549         // Compare azimuth/elevation/range of satellite point : slight difference due to earth flatness
550         // ***************************************************
551         aziElli = topoElliptic.getAzimuth(satPoint, earthElliptic.getBodyFrame(), date);
552         aziSphe = topoSpheric.getAzimuth(satPoint, earthSpheric.getBodyFrame(), date);
553         Assert.assertEquals(aziElli, aziSphe, 1.e-7);
554 
555         eleElli = topoElliptic.getElevation(satPoint, earthElliptic.getBodyFrame(), date);
556         eleSphe = topoSpheric.getElevation(satPoint, earthSpheric.getBodyFrame(), date);
557         Assert.assertEquals(eleElli, eleSphe, 1.e-2);
558 
559         disElli = topoElliptic.getRange(satPoint, earthElliptic.getBodyFrame(), date);
560         disSphe = topoSpheric.getRange(satPoint, earthSpheric.getBodyFrame(), date);
561         Assert.assertEquals(disElli, disSphe, 20.e+3);
562 
563 
564         // *********************
565         // Test at any position
566         // *********************
567         point = new GeodeticPoint(FastMath.toRadians(60), FastMath.toRadians(30.), 0.);
568         topoSpheric  = new TopocentricFrame(earthSpheric, point, "lon 10 lat 45");
569         topoElliptic = new TopocentricFrame(earthElliptic, point, "lon 10 lat 45");
570 
571         // Compare azimuth/elevation/range of satellite point : slight difference
572         // ***************************************************
573         aziElli = topoElliptic.getAzimuth(satPoint, earthElliptic.getBodyFrame(), date);
574         aziSphe = topoSpheric.getAzimuth(satPoint, earthSpheric.getBodyFrame(), date);
575         Assert.assertEquals(aziElli, aziSphe, 1.e-2);
576 
577         eleElli = topoElliptic.getElevation(satPoint, earthElliptic.getBodyFrame(), date);
578         eleSphe = topoSpheric.getElevation(satPoint, earthSpheric.getBodyFrame(), date);
579         Assert.assertEquals(eleElli, eleSphe, 1.e-2);
580 
581         disElli = topoElliptic.getRange(satPoint, earthElliptic.getBodyFrame(), date);
582         disSphe = topoSpheric.getRange(satPoint, earthSpheric.getBodyFrame(), date);
583         Assert.assertEquals(disElli, disSphe, 20.e+3);
584 
585     }
586 
587     @Test
588     public void testFieldEllipticEarth() {
589         doTestFieldEllipticEarth(Decimal64Field.getInstance());
590     }
591 
592     private <T extends CalculusFieldElement<T>> void doTestFieldEllipticEarth(final Field<T> field)  {
593 
594         // zero
595         final T zero = field.getZero();
596 
597         // Elliptic earth shape
598         final OneAxisEllipsoid earthElliptic =
599             new OneAxisEllipsoid(6378136.460, 1 / 298.257222101, itrf);
600 
601         // Field date
602         final FieldAbsoluteDate<T> fieldDate = new FieldAbsoluteDate<>(field, date);
603 
604         // Satellite point
605         // Caution !!! Sat point target shall be the same whatever earth shape chosen !!
606         final FieldGeodeticPoint<T> satPointGeo = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(30.)),
607                                                                            zero.add(FastMath.toRadians(15.)),
608                                                                            zero.add(800000.));
609         final FieldVector3D<T> satPoint = earthElliptic.transform(satPointGeo);
610 
611         // ****************************
612         // Test at equatorial position
613         // ****************************
614         GeodeticPoint point = new GeodeticPoint(FastMath.toRadians(0.), FastMath.toRadians(5.), 0.);
615         TopocentricFrame topoElliptic = new TopocentricFrame(earthElliptic, point, "elliptic, equatorial lon 5");
616         TopocentricFrame topoSpheric = new TopocentricFrame(earthSpheric, point, "spheric, equatorial lon 5");
617 
618         // Compare azimuth/elevation/range of satellite point : shall be strictly identical
619         // ***************************************************
620         T aziElli = topoElliptic.getAzimuth(satPoint, earthElliptic.getBodyFrame(), fieldDate);
621         T aziSphe = topoSpheric.getAzimuth(satPoint, earthSpheric.getBodyFrame(), fieldDate);
622         Assert.assertEquals(aziElli.getReal(), aziSphe.getReal(), Utils.epsilonAngle);
623 
624         T eleElli = topoElliptic.getElevation(satPoint, earthElliptic.getBodyFrame(), fieldDate);
625         T eleSphe = topoSpheric.getElevation(satPoint, earthSpheric.getBodyFrame(), fieldDate);
626         Assert.assertEquals(eleElli.getReal(), eleSphe.getReal(), Utils.epsilonAngle);
627 
628         T disElli = topoElliptic.getRange(satPoint, earthElliptic.getBodyFrame(), fieldDate);
629         T disSphe = topoSpheric.getRange(satPoint, earthSpheric.getBodyFrame(), fieldDate);
630         Assert.assertEquals(disElli.getReal(), disSphe.getReal(), Utils.epsilonTest);
631 
632         // Infinite point separated by -20 deg in longitude
633         // *************************************************
634         FieldGeodeticPoint<T> infPointGeo = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(0.)),
635                                                                      zero.add(FastMath.toRadians(-15.)),
636                                                                      zero.add(1000000000.));
637         FieldVector3D<T> infPoint = earthElliptic.transform(infPointGeo);
638 
639         // Azimuth = pi/2
640         aziElli = topoElliptic.getAzimuth(infPoint, earthElliptic.getBodyFrame(), fieldDate);
641         Assert.assertEquals(3*FastMath.PI/2., aziElli.getReal(), Utils.epsilonAngle);
642 
643         // Site = pi/2 - longitude difference
644         eleElli = topoElliptic.getElevation(infPoint, earthElliptic.getBodyFrame(), fieldDate);
645         Assert.assertEquals(FastMath.abs(infPointGeo.getLongitude().negate().add(point.getLongitude())).negate().add(FastMath.PI/2.).getReal(), eleElli.getReal(), 1.e-2);
646 
647         // Infinite point separated by +20 deg in longitude
648         // *************************************************
649         infPointGeo = new FieldGeodeticPoint<>(zero.add(FastMath.toRadians(0.)),
650                                                zero.add(FastMath.toRadians(25.)),
651                                                zero.add(1000000000.));
652         infPoint = earthElliptic.transform(infPointGeo);
653 
654         // Azimuth = pi/2
655         aziElli = topoElliptic.getAzimuth(infPoint, earthElliptic.getBodyFrame(), fieldDate);
656         Assert.assertEquals(FastMath.PI/2., aziElli.getReal(), Utils.epsilonAngle);
657 
658         // Site = pi/2 - longitude difference
659         eleElli = topoElliptic.getElevation(infPoint, earthElliptic.getBodyFrame(), fieldDate);
660         Assert.assertEquals(FastMath.abs(infPointGeo.getLongitude().negate().add(point.getLongitude())).negate().add(FastMath.PI/2.).getReal(), eleElli.getReal(), 1.e-2);
661 
662         // ************************
663         // Test at polar position
664         // ************************
665         point = new GeodeticPoint(FastMath.toRadians(89.999), FastMath.toRadians(0.), 0.);
666         topoSpheric  = new TopocentricFrame(earthSpheric, point, "lon 0 lat 90");
667         topoElliptic = new TopocentricFrame(earthElliptic, point, "lon 0 lat 90");
668 
669         // Compare azimuth/elevation/range of satellite point : slight difference due to earth flatness
670         // ***************************************************
671         aziElli = topoElliptic.getAzimuth(satPoint, earthElliptic.getBodyFrame(), fieldDate);
672         aziSphe = topoSpheric.getAzimuth(satPoint, earthSpheric.getBodyFrame(), fieldDate);
673         Assert.assertEquals(aziElli.getReal(), aziSphe.getReal(), 1.e-7);
674 
675         eleElli = topoElliptic.getElevation(satPoint, earthElliptic.getBodyFrame(), fieldDate);
676         eleSphe = topoSpheric.getElevation(satPoint, earthSpheric.getBodyFrame(), fieldDate);
677         Assert.assertEquals(eleElli.getReal(), eleSphe.getReal(), 1.e-2);
678 
679         disElli = topoElliptic.getRange(satPoint, earthElliptic.getBodyFrame(), fieldDate);
680         disSphe = topoSpheric.getRange(satPoint, earthSpheric.getBodyFrame(), fieldDate);
681         Assert.assertEquals(disElli.getReal(), disSphe.getReal(), 20.e+3);
682 
683 
684         // *********************
685         // Test at any position
686         // *********************
687         point = new GeodeticPoint(FastMath.toRadians(60), FastMath.toRadians(30.), 0.);
688         topoSpheric  = new TopocentricFrame(earthSpheric, point, "lon 10 lat 45");
689         topoElliptic = new TopocentricFrame(earthElliptic, point, "lon 10 lat 45");
690 
691         // Compare azimuth/elevation/range of satellite point : slight difference
692         // ***************************************************
693         aziElli = topoElliptic.getAzimuth(satPoint, earthElliptic.getBodyFrame(), fieldDate);
694         aziSphe = topoSpheric.getAzimuth(satPoint, earthSpheric.getBodyFrame(), fieldDate);
695         Assert.assertEquals(aziElli.getReal(), aziSphe.getReal(), 1.e-2);
696 
697         eleElli = topoElliptic.getElevation(satPoint, earthElliptic.getBodyFrame(), fieldDate);
698         eleSphe = topoSpheric.getElevation(satPoint, earthSpheric.getBodyFrame(), fieldDate);
699         Assert.assertEquals(eleElli.getReal(), eleSphe.getReal(), 1.e-2);
700 
701         disElli = topoElliptic.getRange(satPoint, earthElliptic.getBodyFrame(), fieldDate);
702         disSphe = topoSpheric.getRange(satPoint, earthSpheric.getBodyFrame(), fieldDate);
703         Assert.assertEquals(disElli.getReal(), disSphe.getReal(), 20.e+3);
704 
705     }
706 
707     @Test
708     public void testPointAtDistance() {
709 
710         RandomGenerator random = new Well1024a(0xa1e6bd5cd0578779l);
711         final OneAxisEllipsoid earth =
712             new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING,
713                                  itrf);
714         final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
715 
716         for (int i = 0; i < 20; ++i) {
717             // we don't need uniform point on the sphere, just a few different test configurations
718             double latitude  = FastMath.PI * (0.5 - random.nextDouble());
719             double longitude = 2 * FastMath.PI * random.nextDouble();
720             TopocentricFrame topo = new TopocentricFrame(earth,
721                                                          new GeodeticPoint(latitude, longitude, 0.0),
722                                                          "topo");
723             Transform transform = earth.getBodyFrame().getTransformTo(topo, date);
724             for (int j = 0; j < 20; ++j) {
725                 double elevation      = FastMath.PI * (0.5 - random.nextDouble());
726                 double azimuth        = 2 * FastMath.PI * random.nextDouble();
727                 double range          = 500000.0 * (1.0 + random.nextDouble());
728                 Vector3D absolutePoint = earth.transform(topo.pointAtDistance(azimuth, elevation, range));
729                 Vector3D relativePoint = transform.transformPosition(absolutePoint);
730                 double rebuiltElevation = topo.getElevation(relativePoint, topo, AbsoluteDate.J2000_EPOCH);
731                 double rebuiltAzimuth   = topo.getAzimuth(relativePoint, topo, AbsoluteDate.J2000_EPOCH);
732                 double rebuiltRange     = topo.getRange(relativePoint, topo, AbsoluteDate.J2000_EPOCH);
733                 Assert.assertEquals(elevation, rebuiltElevation, 1.0e-12);
734                 Assert.assertEquals(azimuth, MathUtils.normalizeAngle(rebuiltAzimuth, azimuth), 1.0e-12);
735                 Assert.assertEquals(range, rebuiltRange, 1.0e-12 * range);
736             }
737         }
738     }
739 
740     @Test
741     public void testIssue145() {
742         Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
743         BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
744                                                Constants.WGS84_EARTH_FLATTENING,
745                                                itrf);
746         TopocentricFrame staFrame = new TopocentricFrame(earth, new GeodeticPoint(0.0, 0.0, 0.0), "test");
747         GeodeticPoint gp = staFrame.computeLimitVisibilityPoint(Constants.WGS84_EARTH_EQUATORIAL_RADIUS+600000,
748                                                                 0.0, FastMath.toRadians(5.0));
749         Assert.assertEquals(0.0, gp.getLongitude(), 1.0e-15);
750         Assert.assertTrue(gp.getLatitude() > 0);
751         Assert.assertEquals(0.0, staFrame.getNorth().distance(Vector3D.PLUS_K), 1.0e-15);
752 
753     }
754 
755     @Before
756     public void setUp() {
757         try {
758 
759             Utils.setDataRoot("regular-data");
760 
761             // Reference frame = ITRF
762             itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
763 
764             // Elliptic earth shape
765             earthSpheric = new OneAxisEllipsoid(6378136.460, 0., itrf);
766 
767             // Reference date
768             date = new AbsoluteDate(new DateComponents(2008, 04, 07),
769                                     TimeComponents.H00,
770                                     TimeScalesFactory.getUTC());
771 
772             // Body mu
773             mu = 3.9860047e14;
774 
775         } catch (OrekitException oe) {
776             Assert.fail(oe.getMessage());
777         }
778     }
779 
780     @Test
781     public void testVisibilityCircle() throws IOException {
782 
783         // a few random from International Laser Ranging Service
784         final BodyShape earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
785                                                      Constants.WGS84_EARTH_FLATTENING,
786                                                      FramesFactory.getITRF(IERSConventions.IERS_2010, true));
787         final TopocentricFrame[] ilrs = {
788             new TopocentricFrame(earth,
789                                  new GeodeticPoint(FastMath.toRadians(52.3800), FastMath.toRadians(3.0649), 133.745),
790                                  "Potsdam"),
791             new TopocentricFrame(earth,
792                                  new GeodeticPoint(FastMath.toRadians(36.46273), FastMath.toRadians(-6.20619), 64.0),
793                                  "San Fernando"),
794             new TopocentricFrame(earth,
795                                  new GeodeticPoint(FastMath.toRadians(35.5331), FastMath.toRadians(24.0705), 157.0),
796                                  "Chania")
797         };
798 
799         PolynomialFunction distanceModel =
800                 new PolynomialFunction(new double[] { 7.0892e+05, 3.1913, -8.2181e-07, 1.4033e-13 });
801         for (TopocentricFrame station : ilrs) {
802             for (double altitude = 500000; altitude < 2000000; altitude += 100000) {
803                 for (double azimuth = 0; azimuth < 2 * FastMath.PI; azimuth += 0.05) {
804                     GeodeticPoint p = station.computeLimitVisibilityPoint(Constants.WGS84_EARTH_EQUATORIAL_RADIUS + altitude,
805                                                                           azimuth, FastMath.toRadians(5.0));
806                     double d = station.getRange(earth.transform(p), earth.getBodyFrame(), AbsoluteDate.J2000_EPOCH);
807                     Assert.assertEquals(distanceModel.value(altitude), d, 40000.0);
808                 }
809             }
810         }
811 
812     }
813 
814     @After
815     public void tearDown() {
816         date = null;
817         itrf = null;
818         earthSpheric = null;
819     }
820 
821 
822 }