Loxodrome.java

/* Copyright 2022-2025 Joseph Reed
 * Licensed to CS GROUP (CS) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * Joseph Reed licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.orekit.bodies;

import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathUtils;

/** Perform calculations on a loxodrome (commonly, a rhumb line) on an ellipsoid.
 * <p>
 * A <a href="https://en.wikipedia.org/wiki/Rhumb_line">loxodrome or rhumb line</a>
 * is an arc on an ellipsoid's surface that intersects every meridian at the same angle.
 *
 * @author Joe Reed
 * @since 11.3
 */
public class Loxodrome {

    /** Threshold for cos angles being equal. */
    private static final double COS_ANGLE_THRESHOLD = 1e-6;

    /** Threshold for when distances are close enough to zero. */
    private static final double DISTANCE_THRESHOLD = 1e-9;

    /** Reference point for the loxodrome. */
    private final GeodeticPoint point;

    /** Azimuth-off-north angle of the loxodrome. */
    private final double azimuth;

    /** Reference body. */
    private final OneAxisEllipsoid body;

    /** Altitude above the body. */
    private final double altitude;

    /** Constructor building a loxodrome from an initial point and an azimuth-off-local-north heading.
     *
     * This method is an equivalent to {@code new Loxodrome(point, azimuth, body, point.getAltitude())}
     *
     * @param point the initial loxodrome point
     * @param azimuth the heading, clockwise angle from north (radians, {@code [0,2π]})
     * @param body ellipsoid body on which the loxodrome is defined
     */
    public Loxodrome(final GeodeticPoint point, final double azimuth, final OneAxisEllipsoid body) {
        this(point, azimuth, body, point.getAltitude());
    }

    /** Constructor building a loxodrome from an initial point and an azimuth-off-local-north heading.
     *
     * @param point the initial loxodrome point
     * @param azimuth the heading, clockwise angle from north (radians, {@code [0,2π]})
     * @param body ellipsoid body on which the loxodrome is defined
     * @param altitude altitude above the reference body
     */
    public Loxodrome(final GeodeticPoint point, final double azimuth, final OneAxisEllipsoid body,
                     final double altitude) {
        this.point    = point;
        this.azimuth  = azimuth;
        this.body     = body;
        this.altitude = altitude;
    }

    /** Get the geodetic point defining the loxodrome.
     * @return the geodetic point defining the loxodrome
     */
    public GeodeticPoint getPoint() {
        return this.point;
    }

    /** Get the azimuth.
     * @return the azimuth
     */
    public double getAzimuth() {
        return this.azimuth;
    }

    /** Get the body on which the loxodrome is defined.
     * @return the body on which the loxodrome is defined
     */
    public OneAxisEllipsoid getBody() {
        return this.body;
    }

    /** Get the altitude above the reference body.
     * @return the altitude above the reference body
     */
    public double getAltitude() {
        return this.altitude;
    }

    /** Calculate the point at the specified distance from the origin point along the loxodrome.
     *
     * A positive distance follows the line in the azimuth direction (i.e. northward for arcs with azimuth
     * angles {@code [3π/2, 2π]} or {@code [0, π/2]}). Negative distances travel in the opposite direction along
     * the rhumb line.
     *
     * Distance is computed at the altitude of the origin point.
     *
     * @param distance the distance to travel (meters)
     * @return the point at the specified distance from the origin
     */
    public GeodeticPoint pointAtDistance(final double distance) {
        if (FastMath.abs(distance) < DISTANCE_THRESHOLD) {
            return this.point;
        }

        // compute the e sin(lat)^2
        final double sinLat = FastMath.sin(point.getLatitude());
        final double eccSinLatSq = body.getEccentricitySquared() * sinLat * sinLat;

        // compute intermediate values
        final double t1 = 1. - body.getEccentricitySquared();
        final double t2 = 1. - eccSinLatSq;
        final double t3 = FastMath.sqrt(t2);

        final double semiMajorAxis = getBody().getEquatorialRadius() + getAltitude();

        final double meridianCurve = (semiMajorAxis * t1) / (t2 * t3);

        final double cosAzimuth = FastMath.cos(azimuth);

        final double lat;
        final double lon;
        if (FastMath.abs(cosAzimuth) < COS_ANGLE_THRESHOLD) {
            lat = point.getLatitude();
            lon = point.getLongitude() + ((distance * FastMath.sin(azimuth) * t3) / (semiMajorAxis * FastMath.cos(point.getLatitude())));
        }
        else {
            final double eccSq34 = 0.75 * body.getEccentricitySquared();
            final double halfEccSq34 = eccSq34 / 2.;
            final double t4 = meridianCurve / (t1 * semiMajorAxis);

            final double latPrime = point.getLatitude() + distance * cosAzimuth / meridianCurve;
            final double latOffset = t4 * (
                ((1. - eccSq34) * (latPrime - point.getLatitude())) +
                        (halfEccSq34 * (FastMath.sin(2. * latPrime) - FastMath.sin(2. * point.getLatitude()))));

            lat = fixLatitude(point.getLatitude() + latOffset);

            final double lonOffset = FastMath.tan(azimuth) * (body.geodeticToIsometricLatitude(lat) - body.geodeticToIsometricLatitude(point.getLatitude()));
            lon = point.getLongitude() + lonOffset;
        }

        return new GeodeticPoint(lat, lon, getAltitude());
    }

    /** Adjust the latitude if necessary, ensuring it's always between -π/2 and +π/2.
     *
     * @param lat the latitude value
     * @return the latitude, within {@code [-π/2,+π/2]}
     */
    static double fixLatitude(final double lat) {
        if (lat < -MathUtils.SEMI_PI) {
            return -MathUtils.SEMI_PI;
        }
        else if (lat > MathUtils.SEMI_PI) {
            return MathUtils.SEMI_PI;
        }
        else {
            return lat;
        }
    }
}