1   /* Copyright 2022-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.util.FastMath;
20  import org.hipparchus.util.MathUtils;
21  
22  /** Perform calculations on a loxodrome (commonly, a rhumb line) on an ellipsoid.
23   * <p>
24   * A <a href="https://en.wikipedia.org/wiki/Rhumb_line">loxodrome or rhumb line</a>
25   * is an arc on an ellipsoid's surface that intersects every meridian at the same angle.
26   *
27   * @author Joe Reed
28   * @since 11.3
29   */
30  public class Loxodrome {
31  
32      /** Threshold for cos angles being equal. */
33      private static final double COS_ANGLE_THRESHOLD = 1e-6;
34  
35      /** Threshold for when distances are close enough to zero. */
36      private static final double DISTANCE_THRESHOLD = 1e-9;
37  
38      /** Reference point for the loxodrome. */
39      private final GeodeticPoint point;
40  
41      /** Azimuth-off-north angle of the loxodrome. */
42      private final double azimuth;
43  
44      /** Reference body. */
45      private final OneAxisEllipsoid body;
46  
47      /** Altitude above the body. */
48      private final double altitude;
49  
50      /** Constructor building a loxodrome from an initial point and an azimuth-off-local-north heading.
51       *
52       * This method is an equivalent to {@code new Loxodrome(point, azimuth, body, point.getAltitude())}
53       *
54       * @param point the initial loxodrome point
55       * @param azimuth the heading, clockwise angle from north (radians, {@code [0,2π]})
56       * @param body ellipsoid body on which the loxodrome is defined
57       */
58      public Loxodrome(final GeodeticPoint point, final double azimuth, final OneAxisEllipsoid body) {
59          this(point, azimuth, body, point.getAltitude());
60      }
61  
62      /** Constructor building a loxodrome from an initial point and an azimuth-off-local-north heading.
63       *
64       * @param point the initial loxodrome point
65       * @param azimuth the heading, clockwise angle from north (radians, {@code [0,2π]})
66       * @param body ellipsoid body on which the loxodrome is defined
67       * @param altitude altitude above the reference body
68       */
69      public Loxodrome(final GeodeticPoint point, final double azimuth, final OneAxisEllipsoid body,
70                       final double altitude) {
71          this.point    = point;
72          this.azimuth  = azimuth;
73          this.body     = body;
74          this.altitude = altitude;
75      }
76  
77      /** Get the geodetic point defining the loxodrome.
78       * @return the geodetic point defining the loxodrome
79       */
80      public GeodeticPoint getPoint() {
81          return this.point;
82      }
83  
84      /** Get the azimuth.
85       * @return the azimuth
86       */
87      public double getAzimuth() {
88          return this.azimuth;
89      }
90  
91      /** Get the body on which the loxodrome is defined.
92       * @return the body on which the loxodrome is defined
93       */
94      public OneAxisEllipsoid getBody() {
95          return this.body;
96      }
97  
98      /** Get the altitude above the reference body.
99       * @return the altitude above the reference body
100      */
101     public double getAltitude() {
102         return this.altitude;
103     }
104 
105     /** Calculate the point at the specified distance from the origin point along the loxodrome.
106      *
107      * A positive distance follows the line in the azimuth direction (i.e. northward for arcs with azimuth
108      * angles {@code [3π/2, 2π]} or {@code [0, π/2]}). Negative distances travel in the opposite direction along
109      * the rhumb line.
110      *
111      * Distance is computed at the altitude of the origin point.
112      *
113      * @param distance the distance to travel (meters)
114      * @return the point at the specified distance from the origin
115      */
116     public GeodeticPoint pointAtDistance(final double distance) {
117         if (FastMath.abs(distance) < DISTANCE_THRESHOLD) {
118             return this.point;
119         }
120 
121         // compute the e sin(lat)^2
122         final double sinLat = FastMath.sin(point.getLatitude());
123         final double eccSinLatSq = body.getEccentricitySquared() * sinLat * sinLat;
124 
125         // compute intermediate values
126         final double t1 = 1. - body.getEccentricitySquared();
127         final double t2 = 1. - eccSinLatSq;
128         final double t3 = FastMath.sqrt(t2);
129 
130         final double semiMajorAxis = getBody().getEquatorialRadius() + getAltitude();
131 
132         final double meridianCurve = (semiMajorAxis * t1) / (t2 * t3);
133 
134         final double cosAzimuth = FastMath.cos(azimuth);
135 
136         final double lat;
137         final double lon;
138         if (FastMath.abs(cosAzimuth) < COS_ANGLE_THRESHOLD) {
139             lat = point.getLatitude();
140             lon = point.getLongitude() + ((distance * FastMath.sin(azimuth) * t3) / (semiMajorAxis * FastMath.cos(point.getLatitude())));
141         }
142         else {
143             final double eccSq34 = 0.75 * body.getEccentricitySquared();
144             final double halfEccSq34 = eccSq34 / 2.;
145             final double t4 = meridianCurve / (t1 * semiMajorAxis);
146 
147             final double latPrime = point.getLatitude() + distance * cosAzimuth / meridianCurve;
148             final double latOffset = t4 * (
149                 ((1. - eccSq34) * (latPrime - point.getLatitude())) +
150                         (halfEccSq34 * (FastMath.sin(2. * latPrime) - FastMath.sin(2. * point.getLatitude()))));
151 
152             lat = fixLatitude(point.getLatitude() + latOffset);
153 
154             final double lonOffset = FastMath.tan(azimuth) * (body.geodeticToIsometricLatitude(lat) - body.geodeticToIsometricLatitude(point.getLatitude()));
155             lon = point.getLongitude() + lonOffset;
156         }
157 
158         return new GeodeticPoint(lat, lon, getAltitude());
159     }
160 
161     /** Adjust the latitude if necessary, ensuring it's always between -π/2 and +π/2.
162      *
163      * @param lat the latitude value
164      * @return the latitude, within {@code [-π/2,+π/2]}
165      */
166     static double fixLatitude(final double lat) {
167         if (lat < -MathUtils.SEMI_PI) {
168             return -MathUtils.SEMI_PI;
169         }
170         else if (lat > MathUtils.SEMI_PI) {
171             return MathUtils.SEMI_PI;
172         }
173         else {
174             return lat;
175         }
176     }
177 }