<commit_before>/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF 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.apache.sis.referencing;

import org.opengis.geometry.coordinate.Position;
import org.opengis.referencing.datum.Ellipsoid;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.apache.sis.referencing.operation.GeodesicException;
import org.apache.sis.internal.referencing.Resources;
import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.math.MathFunctions;

import static java.lang.Math.*;
import static org.apache.sis.internal.metadata.ReferencingServices.NAUTICAL_MILE;


/**
 * Performs geodetic calculations on an ellipsoid. This class overrides the spherical
 * formulas implemented in the parent class, replacing them by ellipsoidal formulas.
 * The methods for direct and inverse geodesic problem use the formulas described in
 * the following publication:
 *
 * <blockquote>
 * <NAME>, 2013.
 * <a href="https://doi.org/10.1007/s00190-012-0578-z">Algorithms for geodesics</a>, SRI International.
 * </blockquote>
 *
 * The following symbols are used with the same meaning than in Karney's article,
 * except λ₁₂ which is represented by ∆λ:
 *
 * <ul>
 *   <li><var>a</var>:  equatorial radius of the ellipsoid of revolution.</li>
 *   <li><var>b</var>:  polar semi-axis of the ellipsoid of revolution.</li>
 *   <li><var>ℯ</var>:  first  eccentricity: ℯ  = √[(a²-b²)/a²].</li>
 *   <li><var>ℯ′</var>: second eccentricity: ℯ′ = √[(a²-b²)/b²].</li>
 *   <li><var>n</var>:  third flattening:    n  = (a-b)/(a+b).</li>
 *   <li><var>E</var>:  the point at which the geodesic crosses the equator in the northward direction.</li>
 *   <li><var>P</var>:  the start point (P₁) or end point (P₂).</li>
 *   <li><var>∆λ</var>: longitude difference between start point and end point (λ₁₂ in Karney).</li>
 *   <li><var>β</var>:  reduced latitude, related to φ geodetic latitude on the ellipsoid.</li>
 *   <li><var>ω</var>:  spherical longitude, related to λ geodetic longitude on the ellipsoid.</li>
 *   <li><var>σ</var>:  spherical arc length, related to distance <var>s</var> on the ellipsoid.</li>
 * </ul>
 *
 * Suffix 1 in variable names denotes values computed using starting point (P₁) and starting azimuth (α₁)
 * while suffix 2 denotes values computed using ending point (P₂) and ending azimuth (α₂).
 * All angular values stored in this class are in radians.
 *
 * <p><b>Limitations:</b>
 * Current implementation is still unable to compute the geodesics in some cases.
 * See <a href="https://issues.apache.org/jira/browse/SIS-467">SIS-467</a>.</p>
 *
 * <p>If the following cases where more than one geodesics exist, current implementation returns an arbitrary one:</p>
 * <ul>
 *   <li>Coincident points (distance is zero but azimuths can be anything).</li>
 *   <li>Starting point and ending points are at opposite poles (there are infinitely many geodesics).</li>
 *   <li>φ₁ = -φ₂ and ∆λ is close to 180° (two geodesics may exist).</li>
 *   <li>∆λ = ±180° (two geodesics may exist).</li>
 * </ul>
 *
 * @author  <NAME> (Geomatys)
 * @author  <NAME> (Geomatys)
 * @version 1.0
 * @since   1.0
 * @module
 */
class GeodesicsOnEllipsoid extends GeodeticCalculator {
    /**
     * Factor by which the accuracy is improved compared to {@link Formulas#ANGULAR_TOLERANCE} value.
     * For example the linear accuracy that {@code GeodesicsOnEllipsoid} aims to achieve is:
     *
     * <blockquote><code>
     * accuracy = {@linkplain Formulas#LINEAR_TOLERANCE} / ACCURACY_IMPROVEMENT
     * </code></blockquote>
     *
     * We take a finer accuracy than default SIS configuration in order to met the accuracy of numbers
     * published in Karney (2013). If this value is modified, the effect can be verified by executing
     * the {@code GeodesicsOnEllipsoidTest} methods that compare computed values against Karney's tables.
     * Remember to update {@link GeodeticCalculator} class javadoc if this value is changed.
     *
     * <p><b>Note:</b> when the iteration loop detects that it reached this requested accuracy, the loop
     * completes the iteration step which was in progress. Consequently the final accuracy is one iteration
     * better than the accuracy computed from this value.</p>
     */
    static final double ACCURACY_IMPROVEMENT = 20;

    /**
     * Difference between ending point and antipode of starting point for considering them as nearly antipodal.
     * This is used only for finding an approximate position before to start iteration using Newton's method,
     * so it is okay if the antipodal approximation has some inaccuracy.
     */
    private static final double NEARLY_ANTIPODAL_Δλ = 0.25 * (PI/180);

    /**
     * Maximal difference (in radians) between two latitudes for enabling the use of simplified formulas.
     * This is used in two contexts:
     *
     * <ul>
     *   <li>Maximal difference between latitude φ₁ and equator for using the equatorial approximation.</li>
     *   <li>Maximal difference between |β₁| and |β₂| for enabling the use of Karney's equation 47.</li>
     * </ul>
     *
     * Those special cases are needed when general formulas produce indeterminations like 0/0.
     * Current angular value corresponds to a distance of 1 millimetre on a planet of the size
     * of Earth, which is about 1.57E-10 radians. This value is chosen empirically by trying to
     * minimize the amount of "No convergence errors" reported by {@code GeodesicsOnEllipsoidTest}
     * in the {@code compareAgainstDataset()} method.
     *
     * <p><b>Note:</b> this is an angular tolerance threshold, but is also used with sine and cosine values
     * because sin(θ) ≈ θ for small angles.</p>
     */
    private static final double LATITUDE_THRESHOLD = 0.001 / (NAUTICAL_MILE*60) * (PI/180);

    /**
     * Whether to include code used for JUnit tests only.
     * This value should be {@code false} in releases.
     *
     * @see #snapshot()
     */
    static final boolean STORE_LOCAL_VARIABLES = true;

    /**
     * The square of eccentricity: ℯ² = (a²-b²)/a² where
     * <var>ℯ</var> is the eccentricity,
     * <var>a</var> is the <cite>semi-major</cite> axis length and
     * <var>b</var> is the <cite>semi-minor</cite> axis length.
     */
    final double eccentricitySquared;
<commit_msg>Now that the development of geodesic distance is completed, set STORE_LOCAL_VARIABLES debug flag to false.<commit_after>/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF 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.apache.sis.referencing;

import org.opengis.geometry.coordinate.Position;
import org.opengis.referencing.datum.Ellipsoid;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.apache.sis.referencing.operation.GeodesicException;
import org.apache.sis.internal.referencing.Resources;
import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.math.MathFunctions;

import static java.lang.Math.*;
import static org.apache.sis.internal.metadata.ReferencingServices.NAUTICAL_MILE;


/**
 * Performs geodetic calculations on an ellipsoid. This class overrides the spherical
 * formulas implemented in the parent class, replacing them by ellipsoidal formulas.
 * The methods for direct and inverse geodesic problem use the formulas described in
 * the following publication:
 *
 * <blockquote>
 * <NAME>, 2013.
 * <a href="https://doi.org/10.1007/s00190-012-0578-z">Algorithms for geodesics</a>, SRI International.
 * </blockquote>
 *
 * The following symbols are used with the same meaning than in Karney's article,
 * except λ₁₂ which is represented by ∆λ:
 *
 * <ul>
 *   <li><var>a</var>:  equatorial radius of the ellipsoid of revolution.</li>
 *   <li><var>b</var>:  polar semi-axis of the ellipsoid of revolution.</li>
 *   <li><var>ℯ</var>:  first  eccentricity: ℯ  = √[(a²-b²)/a²].</li>
 *   <li><var>ℯ′</var>: second eccentricity: ℯ′ = √[(a²-b²)/b²].</li>
 *   <li><var>n</var>:  third flattening:    n  = (a-b)/(a+b).</li>
 *   <li><var>E</var>:  the point at which the geodesic crosses the equator in the northward direction.</li>
 *   <li><var>P</var>:  the start point (P₁) or end point (P₂).</li>
 *   <li><var>∆λ</var>: longitude difference between start point and end point (λ₁₂ in Karney).</li>
 *   <li><var>β</var>:  reduced latitude, related to φ geodetic latitude on the ellipsoid.</li>
 *   <li><var>ω</var>:  spherical longitude, related to λ geodetic longitude on the ellipsoid.</li>
 *   <li><var>σ</var>:  spherical arc length, related to distance <var>s</var> on the ellipsoid.</li>
 * </ul>
 *
 * Suffix 1 in variable names denotes values computed using starting point (P₁) and starting azimuth (α₁)
 * while suffix 2 denotes values computed using ending point (P₂) and ending azimuth (α₂).
 * All angular values stored in this class are in radians.
 *
 * <p><b>Limitations:</b>
 * Current implementation is still unable to compute the geodesics in some cases.
 * See <a href="https://issues.apache.org/jira/browse/SIS-467">SIS-467</a>.</p>
 *
 * <p>If the following cases where more than one geodesics exist, current implementation returns an arbitrary one:</p>
 * <ul>
 *   <li>Coincident points (distance is zero but azimuths can be anything).</li>
 *   <li>Starting point and ending points are at opposite poles (there are infinitely many geodesics).</li>
 *   <li>φ₁ = -φ₂ and ∆λ is close to 180° (two geodesics may exist).</li>
 *   <li>∆λ = ±180° (two geodesics may exist).</li>
 * </ul>
 *
 * @author  <NAME> (Geomatys)
 * @author  <NAME> (Geomatys)
 * @version 1.0
 * @since   1.0
 * @module
 */
class GeodesicsOnEllipsoid extends GeodeticCalculator {
    /**
     * Whether to include code used for JUnit tests only. This field should be
     * set to {@code true} during development and to {@code false} in releases.
     *
     * @see #snapshot()
     */
    static final boolean STORE_LOCAL_VARIABLES = false;

    /**
     * Factor by which the accuracy is improved compared to {@link Formulas#ANGULAR_TOLERANCE} value.
     * For example the linear accuracy that {@code GeodesicsOnEllipsoid} aims to achieve is:
     *
     * <blockquote><code>
     * accuracy = {@linkplain Formulas#LINEAR_TOLERANCE} / ACCURACY_IMPROVEMENT
     * </code></blockquote>
     *
     * We take a finer accuracy than default SIS configuration in order to met the accuracy of numbers
     * published in Karney (2013). If this value is modified, the effect can be verified by executing
     * the {@code GeodesicsOnEllipsoidTest} methods that compare computed values against Karney's tables.
     * Remember to update {@link GeodeticCalculator} class javadoc if this value is changed.
     *
     * <p><b>Note:</b> when the iteration loop detects that it reached this requested accuracy, the loop
     * completes the iteration step which was in progress. Consequently the final accuracy is one iteration
     * better than the accuracy computed from this value.</p>
     */
    static final double ACCURACY_IMPROVEMENT = 20;

    /**
     * Difference between ending point and antipode of starting point for considering them as nearly antipodal.
     * This is used only for finding an approximate position before to start iteration using Newton's method,
     * so it is okay if the antipodal approximation has some inaccuracy.
     */
    private static final double NEARLY_ANTIPODAL_Δλ = 0.25 * (PI/180);

    /**
     * Maximal difference (in radians) between two latitudes for enabling the use of simplified formulas.
     * This is used in two contexts:
     *
     * <ul>
     *   <li>Maximal difference between latitude φ₁ and equator for using the equatorial approximation.</li>
     *   <li>Maximal difference between |β₁| and |β₂| for enabling the use of Karney's equation 47.</li>
     * </ul>
     *
     * Those special cases are needed when general formulas produce indeterminations like 0/0.
     * Current angular value corresponds to a distance of 1 millimetre on a planet of the size
     * of Earth, which is about 1.57E-10 radians. This value is chosen empirically by trying to
     * minimize the amount of "No convergence errors" reported by {@code GeodesicsOnEllipsoidTest}
     * in the {@code compareAgainstDataset()} method.
     *
     * <p><b>Note:</b> this is an angular tolerance threshold, but is also used with sine and cosine values
     * because sin(θ) ≈ θ for small angles.</p>
     */
    private static final double LATITUDE_THRESHOLD = 0.001 / (NAUTICAL_MILE*60) * (PI/180);

    /**
     * The square of eccentricity: ℯ² = (a²-b²)/a² where
     * <var>ℯ</var> is the eccentricity,
     * <var>a</var> is the <cite>semi-major</cite> axis length and
     * <var>b</var> is the <cite>semi-minor</cite> axis length.
     */
    final double eccentricitySquared;

{'fraction_non_alphanumeric': 0.1033731553056922, 'fraction_numerical': 0.010541110330288124, 'mean_word_length': 3.578828828828829, 'pattern_counts': {'":': 0, '<': 205, '<?xml version=': 0, '>': 205, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3232753', 'n_tokens_mistral': 4471, 'n_tokens_neox': 4193, 'n_words': 1624}