/* 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;
  
  /** 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);
         }
     }
 }