LoxodromeArc.java

/* Copyright 2022 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;

/** Loxodrome defined by a start and ending point.
 *
 * @author Joe Reed
 * @since 11.3
 */
public class LoxodromeArc extends Loxodrome {

    /** Threshold for considering latitudes equal, in radians. */
    private static final double LATITUDE_THRESHOLD = 1e-6;

    /** Maximum number of iterations used when computing distance between points. */
    private static final int MAX_ITER = 50;

    /** Ending point of the arc. */
    private final GeodeticPoint endPoint;

    /** Delta longitude, cached, radians. */
    private final double deltaLon;

    /** Cached arc distance, meters. */
    private double distance = -1;

    /** Class constructor where the arc's altitude is the average of the initial and final points.
     * @param point the starting point
     * @param endPoint the ending point
     * @param body the body on which the loxodrome is defined
     */
    public LoxodromeArc(final GeodeticPoint point, final GeodeticPoint endPoint, final OneAxisEllipsoid body) {
        this(point, endPoint, body, (point.getAltitude() + endPoint.getAltitude()) / 2.);
    }

    /** Class constructor.
     * @param point the starting point
     * @param endPoint the ending point
     * @param body the body on which the loxodrome is defined
     * @param altitude the altitude above the reference body (meters)
     */
    public LoxodromeArc(final GeodeticPoint point, final GeodeticPoint endPoint, final OneAxisEllipsoid body,
                        final double altitude) {
        super(point, body.azimuthBetweenPoints(point, endPoint), body, altitude);
        this.endPoint = endPoint;
        this.deltaLon = MathUtils.normalizeAngle(endPoint.getLongitude(), point.getLongitude()) -
                        point.getLongitude();
    }

    /** Get the final point of the arc.
     *
     * @return the ending point of the arc
     */
    public GeodeticPoint getFinalPoint() {
        return this.endPoint;
    }

    /** Compute the distance of the arc along the surface of the ellipsoid.
     *
     * @return the distance (meters)
     */
    public double getDistance() {
        if (distance >= 0) {
            return distance;
        }

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

        // compute intermediate values
        final double t1 = 1. - getBody().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);

        if (FastMath.abs(endPoint.getLatitude() - ptLat) < LATITUDE_THRESHOLD) {
            distance = (semiMajorAxis / t3) * FastMath.abs(FastMath.cos(ptLat) * deltaLon);
        }
        else {
            final double eccSq34 = 0.75 * getBody().getEccentricitySquared();
            final double halfEccSq34 = eccSq34 / 2.;
            final double t6 = 1. / (1. - eccSq34);
            final double t7 = t1 * semiMajorAxis / meridianCurve;

            final double t8 = ptLat + t6 *
                (t7 * (endPoint.getLatitude() - ptLat) + halfEccSq34 * FastMath.sin(ptLat * 2.));
            final double t9 = halfEccSq34 * t6;

            double guess = 0;
            double lat = endPoint.getLatitude();
            for (int i = 0; i < MAX_ITER; i++) {
                guess = lat;
                lat = t8 - t9 * FastMath.sin(2. * guess);

                if (FastMath.abs(lat - guess) < LATITUDE_THRESHOLD) {
                    break;
                }
            }

            final double azimuth = FastMath.atan2(deltaLon,
                getBody().geodeticToIsometricLatitude(lat) - getBody().geodeticToIsometricLatitude(ptLat));
            distance = meridianCurve * FastMath.abs((lat - ptLat) / FastMath.cos(azimuth));
        }

        return distance;
    }

    /** Calculate a point at a specific percentage along the arc.
     *
     * @param fraction the fraction along the arc to compute the point
     * @return the point along the arc
     */
    public GeodeticPoint calculatePointAlongArc(final double fraction) {
        if (fraction == 0.) {
            return getPoint();
        }
        else if (fraction == 1.) {
            return getFinalPoint();
        }
        else {
            final double d = getDistance() * fraction;
            return this.pointAtDistance(d);
        }
    }
}