1   /* Copyright 2002-2025 Joseph Reed
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    * Joseph Reed 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.bodies;
18  
19  import org.hipparchus.geometry.euclidean.threed.Vector3D;
20  import org.hipparchus.util.FastMath;
21  import org.hipparchus.util.MathUtils;
22  import org.hipparchus.util.SinCos;
23  import org.junit.jupiter.api.Assertions;
24  import org.junit.jupiter.api.BeforeEach;
25  import org.junit.jupiter.api.Test;
26  import org.orekit.Utils;
27  import org.orekit.frames.Frame;
28  import org.orekit.frames.FramesFactory;
29  import org.orekit.frames.TopocentricFrame;
30  import org.orekit.time.AbsoluteDate;
31  import org.orekit.utils.IERSConventions;
32  
33  /** Unit tests for {@link Loxodrome} and {@link LoxodromeArc}. */
34  public class LoxodromeTest {
35      private GeodeticPoint newYork;
36      private GeodeticPoint london;
37      private GeodeticPoint berlin;
38      private GeodeticPoint perth;
39      private GeodeticPoint philadelphia;
40      private OneAxisEllipsoid earth;
41  
42      @BeforeEach
43      public void setup() {
44          Utils.setDataRoot("regular-data");
45  
46          newYork = new GeodeticPoint(FastMath.toRadians(40.714268), FastMath.toRadians(-74.005974), 0);
47          london = new GeodeticPoint(FastMath.toRadians(51.5), FastMath.toRadians(-0.16667), 0);
48          berlin = new GeodeticPoint(FastMath.toRadians(52.523403), FastMath.toRadians(13.4114), 0);
49          perth = new GeodeticPoint(FastMath.toRadians(-31.952712), FastMath.toRadians(115.8604796), 0);
50          philadelphia = new GeodeticPoint(FastMath.toRadians(39.952330), FastMath.toRadians(-75.16379), 0);
51  
52          final Frame ecef = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
53          earth  = new OneAxisEllipsoid(6378137, 1. / 298.257223563, ecef);
54      }
55  
56      /** Verify short distance.
57       */
58      @Test
59      public void verifyShortDistance() {
60          // short distance: new york - philadelphia
61          executeTest("new york - philadelphia", newYork, philadelphia, FastMath.toRadians(229.31), 1., 5.);
62      }
63  
64      @Test
65      public void verifyLongDistance() {
66          executeTest("new york - london", newYork, london, FastMath.toRadians(78.09), 30., 6_000.);
67          executeTest("berlin - perth", berlin, perth, FastMath.toRadians(132.89), 90., 35_000.);
68      }
69  
70      @Test
71      public void verifyParallel() {
72          executeTest("perfect parallel",
73                      berlin,
74                      new GeodeticPoint(berlin.getLatitude(), berlin.getLongitude() + 0.125, berlin.getAltitude()),
75                      MathUtils.SEMI_PI, 0.1, 1.0e-20);
76      }
77  
78      @Test
79      public void verifyMeridian() {
80          executeTest("perfect meridian",
81                      berlin,
82                      new GeodeticPoint(berlin.getLatitude() + 0.125, berlin.getLongitude(), berlin.getAltitude()),
83                      0.0, 8.0, 0.003);
84      }
85  
86      void executeTest(final String header, final GeodeticPoint start, final GeodeticPoint stop, final double expectedAzimuth, final double numericalError, final double pointError) {
87          final double az = earth.azimuthBetweenPoints(start, stop);
88          Assertions.assertEquals(expectedAzimuth, az, 1e-4);
89  
90          final Loxodrome lox = new Loxodrome(start, az, earth);
91          final LoxodromeArc arc = new LoxodromeArc(start, stop, earth);
92  
93          final GeodeticPoint expected = highFidelityPointAtDistance(start, az, arc.getDistance());
94          final GeodeticPoint actual = lox.pointAtDistance(arc.getDistance());
95  
96          final double error1 = earth.transform(actual).distance(earth.transform(expected));
97          final double error2 = earth.transform(actual).distance(earth.transform(stop));
98  
99          // over short distances, analytic answer within 1m of numerical computation
100         Assertions.assertTrue(error1 < numericalError,header + " analytic answer not within " + numericalError + "m of numerical answer [error=" + error1 + "m]");
101         // over short distances, analytic answer is within 5m of actual lat/lon
102         Assertions.assertTrue(error2 < pointError,header + " computed destination not within " + pointError + "m of actual answer.[error=" + error2 + "m]");
103         // verify accuracy to 0.001%
104         Assertions.assertEquals(0.0, error1 / arc.getDistance(), 1e-5,header + " accuracy beyond allowable tolerance");
105     }
106 
107     GeodeticPoint highFidelityPointAtDistance(final GeodeticPoint start, final double azimuth, final double distance) {
108         GeodeticPoint point = start;
109         double d = distance;
110         final SinCos scAzimuth = FastMath.sinCos(azimuth);
111         final Vector3D step = new Vector3D(scAzimuth.sin(), scAzimuth.cos(), 0);
112         while (d > 0) {
113             final GeodeticPoint tmp = earth.transform(
114                 new TopocentricFrame(earth, point, "frame").getStaticTransformTo(earth.getBodyFrame(), AbsoluteDate.ARBITRARY_EPOCH).transformPosition(step),
115                 earth.getBodyFrame(),
116                 AbsoluteDate.ARBITRARY_EPOCH);
117             d -= 1.;
118             point = tmp;
119         }
120         return point;
121     }
122 }