This is an automated email from the ASF dual-hosted git repository.
desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git
The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
new 8e51394 Creation of GeodesicsOnEllipsoid class for computing geodescis on ellipsoid of revolution. Based on the algorithm published in Karney C.F.F. (2013), Algorithms for geodesics, SRI International, https://doi.org/10.1007/s00190-012-0578-z
8e51394 is described below
commit 8e51394b9a013751d48b62a804e675114bda10fc
Author: Matthieu_Bastianelli <matthieu.bastianelli@geomatys.com>
AuthorDate: Tue Jun 11 19:04:21 2019 +0200
Creation of GeodesicsOnEllipsoid class for computing geodescis on ellipsoid of revolution.
Based on the algorithm published in Karney C.F.F. (2013), Algorithms for geodesics, SRI International, https://doi.org/10.1007/s00190-012-0578-z
---
.../apache/sis/internal/referencing/Formulas.java | 5 +-
.../sis/referencing/GeodesicsOnEllipsoid.java | 922 +++++++++++++++++++++
.../apache/sis/referencing/GeodeticCalculator.java | 155 ++--
.../referencing/operation/GeodesicException.java | 56 ++
.../operation/projection/Mollweide.java | 4 +-
.../sis/referencing/GeodesicsOnEllipsoidTest.java | 506 +++++++++++
.../sis/referencing/GeodeticCalculatorTest.java | 314 ++++---
.../sis/test/suite/ReferencingTestSuite.java | 1 +
.../java/org/apache/sis/math/MathFunctions.java | 249 +++++-
.../org/apache/sis/math/MathFunctionsTest.java | 43 +
.../test/java/org/apache/sis/test/TestCase.java | 3 +
.../java/org/apache/sis/index/tree/QuadTree.java | 12 +-
12 files changed, 2089 insertions(+), 181 deletions(-)
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
index fe86928..b3b6afe 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
@@ -88,8 +88,11 @@ public final class Formulas extends Static {
* Maximum number of iterations for iterative computations. Defined in this {@code Formulas} class as a default value,
* but some classes may use a derived value (for example twice this amount). This constant is mostly useful for identifying
* places where iterations occur.
+ *
+ * <p>In <cite>Algorithms for geodesics</cite>, Karney (2013) reports that in a tiny fraction of cases up to 16 iterations
+ * is required for the Newton method published in his article. We set {@code MAXIMUM_ITERATIONS} to that maximal value.</p>
*/
- public static final int MAXIMUM_ITERATIONS = 15;
+ public static final int MAXIMUM_ITERATIONS = 16;
/**
* Do not allow instantiation of this class.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
new file mode 100644
index 0000000..8ffb3d8
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
@@ -0,0 +1,922 @@
+/*
+ * 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.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>
+ * Charles F. F. Karney, 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 Matthieu Bastianelli (Geomatys)
+ * @author Martin Desruisseaux (Geomatys)
+ * @version 1.0
+ * @since 1.0
+ * @module
+ */
+class GeodesicsOnEllipsoid extends GeodeticCalculator {
+ /**
+ * 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;
+
+ /**
+ * 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);
+
+ /**
+ * Desired accuracy in radians. 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.
+ *
+ * <p><b>Note:</b> when the iteration loop detects that it reached this accuracy, the loop completes the
+ * iteration step which was in progress. Consequently the final accuracy is one iteration better than this
+ * accuracy.</p>
+ */
+ private static final double ACCURACY = Formulas.ANGULAR_TOLERANCE * (PI/180) / 20;
+
+ /**
+ * 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;
+
+ /**
+ * The square of the second eccentricity: ℯ′² = (a²-b²)/b².
+ */
+ final double secondEccentricitySquared;
+
+ /**
+ * Third flattening of the ellipsoid: n = (a-b)/(a+b).
+ * Not to be confused with the third eccentricity squared ℯ″² = (a²-b²)/(a²+b²) which is not used in this class.
+ */
+ final double thirdFlattening;
+
+ /**
+ * Ration between the semi-minor and semi-major axis b/a.
+ * This is related to flattening <var>f</var> = 1 - b/a.
+ */
+ final double axisRatio;
+
+ /**
+ * Length of semi-minor axis.
+ *
+ * @see org.opengis.referencing.datum.Ellipsoid#getSemiMinorAxis()
+ */
+ private final double semiMinor;
+
+ /**
+ * The α value computed from the starting point and starting azimuth.
+ * We use the sine and cosine instead than the angles because those components are more frequently used than angles.
+ * Those values can be kept constant when computing many end points and end azimuths at different geodesic distances.
+ * The {@link #COEFFICIENTS_FOR_START_POINT} flag specifies whether those fields need to be recomputed.
+ */
+ private double sinα0, cosα0;
+
+ /**
+ * Longitude angle from the equatorial point <var>E</var> to starting point P₁ on the ellipsoid.
+ * This longitude is computed from the ω₁ longitude on auxiliary sphere, which is itself computed
+ * from α₀, α₁, β₁ and ω₁ values computed from the starting point and starting azimuth.
+ *
+ * @see #sphericalToGeodeticLongitude(double, double)
+ */
+ private double λ1E;
+
+ /**
+ * Ellipsoidal arc length <var>s₁</var>/<var>b</var> computed from the spherical arc length <var>σ₁</var>.
+ * The <var>σ₁</var> value is an arc length on the auxiliary sphere between equatorial point <var>E</var>
+ * (the point on equator with forward direction toward azimuth α₀) and starting point P₁. This is computed
+ * by I₁(σ₁) in Karney equation 15 and is saved for reuse when the starting point and azimuth do not change.
+ *
+ * @see #sphericalToEllipsoidalAngle(double, boolean)
+ */
+ private double I1_σ1;
+
+ /**
+ * The term to be raised to powers (ε⁰, ε¹, ε², ε³, …) in series expansions. Defined in Karney equation 16 as
+ * ε = (√[k²+1] - 1) / (√[k²+1] + 1) where k² = {@linkplain #secondEccentricitySquared ℯ′²}⋅cos²α₀ (Karney 9).
+ */
+ private double ε;
+
+ /**
+ * Coefficients in series expansion of the form A × (σ + ∑Cₓ⋅sin(…⋅σ)).
+ * There is 3 series expansions used in this class:
+ *
+ * <ul>
+ * <li><var>A₁</var> and <var>C₁ₓ</var> (Karney 15) are used for arc length conversions from auxiliary sphere to ellipsoid.</li>
+ * <li><var>A₂</var> and <var>C₂ₓ</var> (Karney 41) are used in calculation of reduced length.</li>
+ * <li><var>A₃</var> and <var>C₃ₓ</var> (Karney 23) are used for calculation of geodetic longitude from spherical angles.</li>
+ * </ul>
+ *
+ * The <var>Cₓₓ</var> coefficients are hard-coded, except the <var>C₃ₓ</var> coefficients
+ * (used together with <var>A₃</var>) which depend on {@link #sinα0} and {@link #cosα0}.
+ *
+ * <p>All those coefficients must be recomputed when the starting point or starting azimuth change.
+ * The {@link #COEFFICIENTS_FOR_START_POINT} flag specifies whether those fields need to be recomputed.</p>
+ *
+ * @see #computeSeriesExpansionCoefficients()
+ */
+ private double A1, A2, A3, C31, C32, C33, C34, C35;
+
+ /**
+ * Constructs a new geodetic calculator expecting coordinates in the supplied CRS.
+ *
+ * @param crs the referencing system for the {@link Position} arguments and return values.
+ */
+ GeodesicsOnEllipsoid(final CoordinateReferenceSystem crs) {
+ super(crs);
+ final double a = ellipsoid.getSemiMajorAxis();
+ final double b = ellipsoid.getSemiMinorAxis();
+ final double Δ2 = a*a - b*b;
+ eccentricitySquared = Δ2 / (a*a);
+ secondEccentricitySquared = Δ2 / (b*b);
+ thirdFlattening = (a - b) / (a + b);
+ axisRatio = b / a;
+ semiMinor = b;
+ }
+
+ /**
+ * Computes series expansions coefficients.
+ *
+ * <p><b>Preconditions:</b> The {@link #sinα0} and {@link #cosα0} fields shall be set before to invoke this method.
+ * It is caller's responsibility to ensure that sin(α₀)² + cos(α₀)² ≈ 1 (this is verified in assertion).</p>
+ *
+ * <p><b>Post-conditions:</b> this method sets the {@link #ε}, {@link #A1}, {@link #A2}, {@link #A3},
+ * {@link #C31}, {@link #C32}, {@link #C33}, {@link #C34}, {@link #C35} fields.</p>
+ *
+ * It is caller's responsibility to invoke {@code setValid(COEFFICIENTS_FOR_START_POINT)} after it has updated
+ * {@link #λ1E} and {@link #I1_σ1}.
+ *
+ * @return k² = ℯ′²⋅cos²α₀ term (Karney 9).
+ */
+ private double computeSeriesExpansionCoefficients() {
+ assert abs(sinα0*sinα0 + cosα0*cosα0 - 1) <= 1E-10; // Arbitrary threshold.
+ final double k2 = secondEccentricitySquared * (cosα0 * cosα0); // (Karney 9)
+ final double h = sqrt(1 + k2);
+ ε = (h - 1) / (h + 1); // (Karney 16)
+ final double ε2 = ε*ε;
+ A1 = (((( 1./256)*ε2 + 1./64)*ε2 + 1./4)*ε2 + 1) / (1 - ε); // (Karney 17)
+ A2 = ((((25./256)*ε2 + 9./64)*ε2 + 1./4)*ε2 + 1) * (1 - ε); // (Karney 42)
+ /*
+ * Compute A₃ from Karney equation 24 with signs redistributed in a different but equivalent way.
+ * We use + or - operations in a way where multiplication factors of n are positive. By increasing
+ * the number of times where the same value is used (for example a single 1/8 instead of two values
+ * +1/8 and -1/8), our hope is that the compiler can reuse more often values that are already loaded
+ * in registers. Maybe some compilers could even detect when the same multiplication is done two times.
+ *
+ * Note: for each line below, we could compute at construction time the parts at the right of ε.
+ * However since they are only a few additions and multiplications, it may not be worth to vastly
+ * increase the number of fields in `GeodesicOnEllipsoid` for that purpose, especially if compiler
+ * can optimize itself those duplicated multiplications.
+ */
+ final double n = thirdFlattening;
+ A3 = 1 - ε*(1./2 - n*(1./2)
+ + ε*(1./4 + n*(1./8 - n*(3./8 ))
+ + ε*(1./16 + n*(3./16 + n*(1./16))
+ + ε*(3./64 + n*(1./32)
+ + ε*(3./128)))));
+ /*
+ * Karney equation 25 with fε = εⁿ where n = 1, 2, 3, …
+ *
+ * TODO: combine using the technic described in sphericalToEllipsoidalAngle,
+ * then make constant positive.
+ */
+ double fε;
+ C31 = (fε = ε) * ( 1./4 + n*(-1./4)
+ + ε * ( 1./8 + n*( n*(-1./8 ))
+ + ε * ( 3./64 + n*(3./64 + n*(-1./64))
+ + ε * ( 5./128 + n*(1./64)
+ + ε * ( 3./128)))));
+ C32 = (fε *= ε) * ( 1./16 + n*(-3./32 + n*( 1./32))
+ + ε * ( 3./64 + n*(-1./32 + n*(-3./64))
+ + ε * ( 3./128 + n*( 1./128)
+ + ε * ( 5./256))));
+ C33 = (fε *= ε) * ( 5./192 + n*(-3./64 + n*( 5./192))
+ + ε * ( 3./128 + n*(-5./192)
+ + ε * ( 7./512)));
+ C34 = (fε *= ε) * ( 7./512 + n*(-7./256)
+ + ε * ( 7./512));
+ C35 = (fε * ε) * (21./2560);
+ return k2;
+ }
+
+ /**
+ * Computes a term in calculation of ellipsoidal arc length <var>s</var> from the given spherical arc length <var>σ</var>.
+ * The <var>σ</var> value is an arc length on the auxiliary sphere between equatorial point <var>E</var>
+ * (the point on equator with forward direction toward azimuth α₀) and another point.
+ *
+ * <p>If the {@code minusI2} argument is {@code false}, then this method computes the ellipsoidal arc length <var>s</var>
+ * from the given spherical arc length <var>σ</var>. This value is computed by Karney's equation 7:
+ *
+ * <blockquote>s/b = I₁(σ)</blockquote>
+ *
+ * If the {@code minusI2} argument is {@code true}, then this method is modified for computing
+ * Karney's equation 40 instead (a term in calculation of reduced length <var>m</var>):
+ *
+ * <blockquote>J(σ) = I₁(σ) - I₂(σ)</blockquote>
+ *
+ * <p><b>Precondition:</b> {@link #computeSeriesExpansionCoefficients()} must have been invoked before this method.</p>
+ *
+ * @param σ arc length on the auxiliary sphere since equatorial point <var>E</var>.
+ * @param minusI2 whether to subtract I₂(σ) after calculation of I₁(σ).
+ * @return I₁(σ) if {@code minusI2} is false, or J(σ) if {@code minusI2} is true.
+ */
+ private double sphericalToEllipsoidalAngle(final double σ, final boolean minusI2) {
+ /*
+ * Derived from Karney (2013) equation 18 with θ = 2σ:
+ *
+ * ∑Cₓ⋅sin(x⋅θ) ≈ (-1/2 ⋅ε + 3/16 ⋅ε³ + -1/32 ⋅ε⁵)⋅sin(1θ) + // C₁₁
+ * (-1/16⋅ε² + 1/32 ⋅ε⁴ + -9/2048⋅ε⁶)⋅sin(2θ) + // C₁₂
+ * ( -1/48 ⋅ε³ + 3/256 ⋅ε⁵)⋅sin(3θ) + // C₁₃
+ * ( -5/512⋅ε⁴ + 3/512 ⋅ε⁶)⋅sin(4θ) + // C₁₄
+ * ( -7/1280⋅ε⁵)⋅sin(5θ) + // C₁₅
+ * ( -7/2048⋅ε⁶)⋅sin(6θ) // C₁₆
+ *
+ * For performance and simplicity, we rewrite the equations using trigonometric identities
+ * as described in Snyder 3-34 and 3-35, expanded with two more terms. Given:
+ *
+ * s = A⋅sin(θ) + B⋅sin(2θ) + C⋅sin(3θ) + D⋅sin(4θ) + E⋅sin(5θ) + F⋅sin(6θ)
+ *
+ * We rewrite as:
+ *
+ * s = sinθ⋅(A′ + cosθ⋅(B′ + cosθ⋅(C′ + cosθ⋅(D′ + cosθ⋅(E′ + cosθ⋅F′)))))
+ * A′ = A - C + E
+ * B′ = 2B - 4D + 6F
+ * C′ = 4C - 12E
+ * D′ = 8D - 32F
+ * E′ = 16E
+ * F′ = 32F
+ *
+ * See https://svn.apache.org/repos/asf/sis/analysis/Map%20projection%20formulas.ods for calculations.
+ * Additions are done with smallest values first for better precision.
+ *
+ * See Karney 2011, Geodesics on an ellipsoid of revolution, https://arxiv.org/pdf/1102.1215.pdf
+ * if more coefficients (up to ε⁸) are desired. Those coefficients may be needed in a future SIS
+ * version if we want to support celestial bodies with higher flattening factor.
+ * Issue tracker: https://issues.apache.org/jira/browse/SIS-465
+ */
+ final double ε2 = ε*ε;
+ final double ε3 = ε*ε2;
+ final double ε4 = ε2*ε2;
+ final double ε5 = ε2*ε3;
+ final double ε6 = ε3*ε3;
+ final double θ = σ * 2;
+ final double sinθ = sin(θ);
+ final double cosθ = cos(θ);
+ double m = A1 * (σ + sinθ * (-31./640 * ε5 + 5./24 * ε3 + -1./2 * ε
+ + cosθ * (-27./512 * ε6 + 13./128 * ε4 + -1./8 * ε2
+ + cosθ * ( 9./80 * ε5 + -1./12 * ε3
+ + cosθ * ( 5./32 * ε6 + -5./64 * ε4
+ + cosθ * ( -7./80 * ε5
+ + cosθ * ( -7./64 * ε6)))))));
+ if (minusI2) {
+ m -= A2 * (σ + sinθ * ( 39./640 * ε5 + -1./24 * ε3 + 1./2 * ε
+ + cosθ * (105./512 * ε6 + -27./128 * ε4 + 3./8 * ε2
+ + cosθ * ( -41./80 * ε5 + 5./12 * ε3
+ + cosθ * ( -35./32 * ε6 + 35./64 * ε4
+ + cosθ * ( 63./80 * ε5
+ + cosθ * ( 77./64 * ε6)))))));
+ }
+ return m;
+ }
+
+ /**
+ * Computes the spherical arc length <var>σ</var> from the an ellipsoidal arc length <var>s</var>.
+ * This method is the converse of {@link #sphericalToEllipsoidalAngle(double, boolean)}.
+ * Input is τ = s/(b⋅A₁).
+ *
+ * <p><b>Precondition:</b> {@link #computeSeriesExpansionCoefficients()} must have been invoked before this method.</p>
+ *
+ * @param τ arc length on the ellipsoid since equatorial point <var>E</var> divided by {@link #A1}.
+ * @return arc length <var>σ</var> on the auxiliary sphere.
+ *
+ * @see #sphericalToEllipsoidalAngle(double, boolean)
+ */
+ private double ellipsoidalToSphericalAngle(final double τ) {
+ assert !isInvalid(COEFFICIENTS_FOR_START_POINT);
+ /*
+ * Derived from Karney (2013) equation 21 with θ = 2τ:
+ *
+ * ∑C′ₓ⋅sin(x⋅θ) ≈ (1/2 ⋅ε + -9/32 ⋅ε³ + 205/1536 ⋅ε⁵)⋅sin(1θ) + // C′₁₁
+ * (5/16⋅ε² + -37/96 ⋅ε⁴ + 1335/4096 ⋅ε⁶)⋅sin(2θ) + // C′₁₂
+ * ( 29/96 ⋅ε³ + -75/128 ⋅ε⁵)⋅sin(3θ) + // C′₁₃
+ * ( 539/1536⋅ε⁴ + -2391/2560 ⋅ε⁶)⋅sin(4θ) + // C′₁₄
+ * ( 3467/7680 ⋅ε⁵)⋅sin(5θ) + // C′₁₅
+ * ( 38081/61440⋅ε⁶)⋅sin(6θ) // C′₁₆
+ *
+ * For performance and simplicity, we rewrite the equations using trigonometric identities
+ * using the same technic than the one documented above in sphericalToEllipsoidalAngle(σ).
+ */
+ final double ε2 = ε*ε;
+ final double ε3 = ε*ε2;
+ final double ε4 = ε2*ε2;
+ final double ε5 = ε2*ε3;
+ final double ε6 = ε3*ε3;
+ final double θ = τ * 2;
+ final double cosθ = cos(θ);
+ return τ + sin(θ) * ( 281./240 * ε5 + -7./12 * ε3 + 1./2 * ε
+ + cosθ * ( 20753./2560 * ε6 + -835./384 * ε4 + 5./8 * ε2
+ + cosθ * ( -4967./640 * ε5 + 29./24 * ε3
+ + cosθ * (-52427./1920 * ε6 + 539./192 * ε4
+ + cosθ * ( 3467./480 * ε5
+ + cosθ * ( 38081./1920 * ε6))))));
+ }
+
+ /**
+ * Computes the longitude angle λ from the equatorial point <var>E</var> to a point on the ellipsoid
+ * specified by the ω longitude on auxiliary sphere.
+ *
+ * <p><b>Precondition:</b> {@link #computeSeriesExpansionCoefficients()} must have been invoked before this method.</p>
+ *
+ * @param ω longitude on the auxiliary sphere.
+ * @param σ spherical arc length from equatorial point E to point on auxiliary sphere
+ * along a geodesic with azimuth α₀ at equator.
+ * @return geodetic longitude λ from the equatorial point E to the point.
+ *
+ * @todo combine using the technic described in sphericalToEllipsoidalAngle,
+ */
+ private double sphericalToGeodeticLongitude(final double ω, final double σ) {
+ /*
+ * Derived from Karney (2013) equation 23.
+ *
+ * I₃(σ) = A₃⋅(σ + ∑C₃ₓ⋅sin(x⋅θ))
+ */
+ final double I3 = A3*(C35 * sin(2 * 5 * σ)
+ + C34 * sin(2 * 4 * σ)
+ + C33 * sin(2 * 3 * σ)
+ + C32 * sin(2 * 2 * σ)
+ + C31 * sin(2 * 1 * σ) + σ);
+ if (STORE_LOCAL_VARIABLES) store("I₃(σ)", I3);
+ return ω - sinα0 * I3 * (1 - axisRatio);
+ }
+
+ /**
+ * Computes the end point from the start point, the azimuth and the geodesic distance.
+ * This method should be invoked if the end point or ending azimuth is requested while
+ * {@link #END_POINT} validity flag is not set.
+ *
+ * <p>This implementation computes {@link #φ2}, {@link #λ2} and azimuths using the method
+ * described in Karney (2013) <cite>Algorithms for geodesics</cite>. The coefficients of
+ * Fourier and Taylor series expansions are given by equations 15 and 17.</p>
+ *
+ * @throws IllegalStateException if the start point, azimuth or distance has not been set.
+ */
+ @Override
+ final void computeEndPoint() {
+ canComputeEndPoint(); // May throw IllegalStateException.
+ if (isInvalid(COEFFICIENTS_FOR_START_POINT)) {
+ final double vm = hypot(dφ1, dλ1);
+ final double sinα1 = dλ1 / vm; // α is a geographic (not arithmetic) angle.
+ final double cosα1 = dφ1 / vm;
+ final double tanβ1 = axisRatio * tan(φ1); // β₁ is reduced latitude (Karney 6).
+ final double cosβ1 = 1 / sqrt(1 + tanβ1*tanβ1);
+ final double sinβ1 = tanβ1 * cosβ1;
+ sinα0 = sinα1 * cosβ1; // Azimuth at equator (Karney 5) as geographic angle.
+ cosα0 = hypot(cosα1, sinα1*sinβ1); // = sqrt(1 - sinα0*sinα0) with less rounding errors.
+ /*
+ * Note: Karney said that for equatorial geodesics (φ₁=0 and α₁=π/2), calculation of σ₁ is indeterminate
+ * and σ₁=0 should be taken. The indetermination appears as atan2(0,0). However this expression evaluates
+ * to 0 in Java, as required by Math.atan2(y,x) specification, which is what we want. So there is no need
+ * for a special case.
+ */
+ final double cosα1_cosβ1 = cosα1*cosβ1;
+ final double σ1 = atan2(sinβ1, cosα1_cosβ1); // Arc length on great circle (Karney 11).
+ final double ω1 = atan2(sinβ1*sinα0, cosα1_cosβ1);
+ final double k2 = computeSeriesExpansionCoefficients();
+ λ1E = sphericalToGeodeticLongitude(ω1, σ1);
+ I1_σ1 = sphericalToEllipsoidalAngle(σ1, false);
+ setValid(COEFFICIENTS_FOR_START_POINT);
+ if (STORE_LOCAL_VARIABLES) { // For comparing values with Karney table 2.
+ store("β₁", atan2(sinβ1, cosβ1));
+ store("σ₁", σ1);
+ store("ω₁", ω1);
+ store("k²", k2);
+ snapshot();
+ }
+ }
+ /*
+ * Distance from equatorial point E to ending point P₂ along the geodesic: s₂ = s₁ + ∆s.
+ */
+ final double s2b = I1_σ1 + geodesicDistance / semiMinor; // (Karney 18) + ∆s/b
+ final double σ2 = ellipsoidalToSphericalAngle(s2b / A1); // (Karney 21)
+ final double sinσ2 = sin(σ2);
+ final double cosσ2 = cos(σ2);
+ /*
+ * Azimuth at ending point α₂ = atan2(sinα₀, cosα₀⋅cosσ₂) (Karney 14 adapted for ending point).
+ * We do not need atan2(y,x) because we keep x and y components separated. It is not necessary to
+ * normalize to a vector of length 1.
+ */
+ dλ2 = sinα0; // cos(π/2 − θ) = sin(θ)
+ dφ2 = cosα0 * cosσ2; // sin(π/2 − θ) = cos(θ)
+ /*
+ * Ending point coordinates on auxiliary sphere: Latitude β is given by Karney equation 13:
+ *
+ * β₂ = atan2(cos(α₀)⋅sin(σ₂), hypot(cos(α₀)⋅cos(σ₂), sin(α₀))
+ *
+ * We replace cos(α₀)⋅cos(σ₂) by dφ2 since we computed it above. Then we avoid the call to
+ * atan2(y,x) by storing directly the y and x values. Note that `sinβ2` and `cosβ2` are not
+ * really sine and cosine since we do not normalize them to sin² + cos² = 1. We do not need
+ * to normalize because we use either a ratio of those 2 quantities or give them to atan2(…).
+ */
+ final double sinβ2 = cosα0 * sinσ2;
+ final double cosβ2 = hypot(dφ2, sinα0);
+ final double ω2 = atan2(sinα0*sinσ2, cosσ2); // (Karney 12).
+ /*
+ * Convert reduced longitude ω and latitude β on auxiliary sphere
+ * to geodetic longitude λ and latitude φ.
+ */
+ final double λ2E = sphericalToGeodeticLongitude(ω2, σ2);
+ λ2 = IEEEremainder(λ2E - λ1E + λ1, 2*PI);
+ φ2 = atan(sinβ2 / (cosβ2 * axisRatio)); // (Karney 6).
+ setValid(END_POINT | ENDING_AZIMUTH);
+ if (STORE_LOCAL_VARIABLES) { // For comparing values with Karney table 2.
+ store("s₂", s2b * semiMinor);
+ store("τ₂", s2b / A1);
+ store("σ₂", σ2);
+ store("α₂", atan2(sinα0, cosα0*cosσ2));
+ store("β₂", atan2(sinβ2, cosβ2));
+ store("ω₂", ω2);
+ store("λ₂", λ2E);
+ store("Δλ", λ2E - λ1E);
+ }
+ }
+
+ /**
+ * Computes the geodesic distance and azimuths from the start point and end point.
+ * This method should be invoked if the distance or an azimuth is requested while
+ * {@link #STARTING_AZIMUTH}, {@link #ENDING_AZIMUTH} or {@link #GEODESIC_DISTANCE}
+ * validity flag is not set.
+ *
+ * <p>Reminder: given <var>P₁</var> the starting point and <var>E</var> the intersection of the geodesic with equator:</p>
+ * <ul>
+ * <li><var>α₁</var> is the azimuth (0° oriented North) of the geodesic from <var>E</var> to <var>P₁</var>.</li>
+ * <li><var>σ₁</var> is the spherical arc length (in radians) between <var>E</var> and <var>P₁</var>.</li>
+ * <li><var>ω₁</var> is the spherical longitude (in radians) on the auxiliary sphere at <var>P₁</var>.
+ * Spherical longitude is the angle formed by the meridian of <var>E</var> and the meridian of <var>P₁</var>.</li>
+ * </ul>
+ *
+ * @throws IllegalStateException if the distance or azimuth has not been set.
+ * @throws GeodesicException if an azimuth or the distance can not be computed.
+ */
+ @Override
+ final void computeDistance() throws GeodesicException {
+ canComputeDistance();
+ /*
+ * The algorithm in this method requires the following canonical configuration:
+ *
+ * Negative latitude of starting point: φ₁ ≦ 0
+ * Ending point latitude smaller in magnitude: φ₁ ≦ φ₂ ≦ -φ₁
+ * Positive longitude difference: 0 ≦ ∆λ ≦ π
+ * (Consequence of above): 0 ≦ α₀ ≦ π/2
+ *
+ * If the given points do not met above conditions, then we need to swap start and end points or to
+ * swap coordinate signs. We apply those changes on local variables only, not on the class fields.
+ * We will need to apply the converse of those changes in the final results.
+ */
+ double φ1 = this.φ1;
+ double φ2 = this.φ2;
+ double Δλ = IEEEremainder(λ2 - λ1, 2*PI); // In [-π … +π] range.
+ final boolean swapPoints = abs(φ2) > abs(φ1);
+ if (swapPoints) {
+ φ1 = this.φ2;
+ φ2 = this.φ1;
+ Δλ = -Δλ;
+ }
+ final boolean inverseLatitudeSigns = φ1 > 0;
+ if (inverseLatitudeSigns) {
+ φ1 = -φ1;
+ φ2 = -φ2;
+ }
+ final boolean inverseLongitudeSigns = Δλ < 0;
+ if (inverseLongitudeSigns) {
+ Δλ = -Δλ;
+ }
+ /*
+ * Compute reduced latitudes β (Karney 6). Actually we don't need the β angles
+ * (except for a special case), but rather their sine and cosine values.
+ */
+ final double tanβ1 = axisRatio * tan(φ1);
+ final double tanβ2 = axisRatio * tan(φ2);
+ final double cosβ1 = 1 / sqrt(1 + tanβ1*tanβ1);
+ final double cosβ2 = 1 / sqrt(1 + tanβ2*tanβ2);
+ final double sinβ1 = tanβ1 * cosβ1;
+ final double sinβ2 = tanβ2 * cosβ2;
+ /*
+ * Compute an approximation of the azimuth α₁ at starting point. This estimation will be refined by iteration
+ * in the loop later, but that iteration will not converge if the first α₁ estimation is not good enough. The
+ * general formula does not give good α₁ initial value for antipodal points, so we need to check for special
+ * cases first:
+ *
+ * 1) Nearly antipodal points with φ₁ = -φ₂.
+ * 2) Nearly antipodal points with φ₁ slightly different than -φ₂.
+ * 3) Equatorial case: φ₁ = φ₂ = 0 but restricted to ∆λ ≤ (1-f)⋅π.
+ * 4) Meridional case: ∆λ = 0 or ∆λ = π (handled by general case in this method).
+ */
+ if (φ1 > -LATITUDE_THRESHOLD) { // Sufficient test because φ₁ ≦ 0 and |φ₂| ≦ φ₁
+ /*
+ * Points on equator but not nearly anti-podal. The geodesic is an arc on equator and the azimuths
+ * are α₁ = α₂ = ±90°. We need this special case because when φ = 0, the general case get sinβ = 0
+ * then σ = atan2(sinβ, cosα⋅cosβ) = 0, which result in a distance of 0. I have not yet understood
+ * how to use the general formulas in such case. This code is a workaround in the meantime.
+ *
+ * See https://issues.apache.org/jira/browse/SIS-467
+ */
+ if (Δλ > axisRatio * PI) {
+ // Karney's special case documented before equation 45.
+ throw new GeodesicException("Can not compute geodesics for antipodal points on equator.");
+ }
+ final double Δφ = φ2 - φ1;
+ geodesicDistance = hypot(Δλ, Δφ) * ellipsoid.getSemiMajorAxis();
+ dλ1 = dλ2 = (inverseLongitudeSigns ^ swapPoints) ? -Δλ : Δλ;
+ dφ1 = dφ2 = (inverseLatitudeSigns ^ swapPoints) ? -Δφ : Δφ;
+ setValid(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE);
+ return;
+ }
+ double α1;
+ if (Δλ >= PI - NEARLY_ANTIPODAL_Δλ && abs(φ1 + φ2) <= NEARLY_ANTIPODAL_Δλ) {
+ /*
+ * Nearly antipodal points. Karney's equations 53 are reproduced on the left side
+ * (with f = 1 - b/a), which we replace by the equations on the right side below:
+ *
+ * Δ = f⋅a⋅π⋅cos²β₁ ┃ Δ₁ = f⋅π⋅cosβ₁
+ * ∆λ = π + Δ/(a⋅cosβ₁)⋅x ┃ ∆λ = π + Δ₁⋅x
+ * β₂ = -β₁ + (Δ/a)⋅y ┃ β₂ = -β₁ + (Δ₁⋅cosβ₁)⋅y
+ */
+ final double β1 = atan2(sinβ1, cosβ1);
+ final double β2 = atan2(sinβ2, cosβ2);
+ final double Δ1 = (1 - axisRatio) * PI * cosβ1; // Differ from Karney by a⋅cosβ₁ factor.
+ final double y = (β2 + β1) / (Δ1*cosβ1);
+ final double x = (PI - Δλ) / Δ1; // Opposite sign of Karney. We have x ≧ 0.
+ final double x2 = x*x;
+ final double y2 = y*y;
+ if (STORE_LOCAL_VARIABLES) { // For comparing values with Karney table 4.
+ store("x", -x); // Sign used by Karney.
+ store("y", y);
+ }
+ if (y2 < 1E-12) { // Empirical threshold. See μ(…) for more information.
+ α1 = (x2 > 1) ? PI/2 : atan(x / sqrt(1 - x2)); // (Karney 57) with opposite sign of x. Result in α₁ ≧ 0.
+ } else {
+ final double μ = μ(x2, y2);
+ α1 = atan2(x*μ, (y*(1+μ))); // (Karney 56) with opposite sign of x.
+ if (STORE_LOCAL_VARIABLES) {
+ store("μ", μ);
+ }
+ }
+ } else {
+ /*
+ * Usual case (non-antipodal). Estimation is based on variation of geodetic longitude λ
+ * and spherical longitude ω on the auxiliary sphere. Karney makes a special case for
+ * Δλ = 0 and Δλ = π by defining α₁ = Δλ. However formulas below produce the same result.
+ */
+ double w = (cosβ1 + cosβ2) / 2;
+ w = sqrt(1 - eccentricitySquared * (w*w));
+ final double Δω = Δλ / w;
+ final double cω = cos(Δω);
+ final double sω = sin(Δω);
+ final double αx = cosβ1*sinβ2 - sinβ1*cosβ2*cω;
+ final double αy = cosβ2*sω;
+ α1 = atan2(αy, αx); // (Karney 49)
+ if (STORE_LOCAL_VARIABLES) { // For comparing values with Karney table 3.
+ store("ωb", w);
+ store("Δω", Δω);
+ store("α₂", atan2(cosβ1*sω, cosβ1*cω*sinβ2 - sinβ1*cosβ2)); // (Karney 50)
+ store("Δσ", atan2(hypot(αy, αx), cosβ1*cω*cosβ2 + sinβ1*sinβ2)); // (Karney 51)
+ // For small distances we could stop here. But it is easier to let iteration does its work.
+ }
+ }
+ /*
+ * Stores locale variables for comparison against Karney tables 4, 5 and 6. Values β₁ and β₂ are keep
+ * constant during all this method. Value α₁ is a first estimation and will be updated during iteration.
+ * Not that because Karney separate calculation of α₁ and remaining calculation in two separated tables,
+ * we need to truncate α₁ to the same number of digits than Karney in order to get the same numbers in
+ * the rest of this method.
+ */
+ if (STORE_LOCAL_VARIABLES) {
+ store("β₁", atan(tanβ1));
+ store("β₂", atan(tanβ2));
+ store("α₁", α1);
+ α1 = computedToGiven(α1);
+ }
+ /*
+ * Refine α₁ using Newton's method. Karney proposes an hybrid approach: approximate α₁,
+ * then compute σ₁, σ₂, ω₁ and ω₂ (actually the sine and cosine of those angles in our
+ * implementation).
+ */
+ double σ1, σ2;
+ for (int nbIter = Formulas.MAXIMUM_ITERATIONS;;) {
+ final double sinα1 = sin(α1);
+ final double cosα1 = cos(α1);
+ sinα0 = sinα1 * cosβ1;
+ cosα0 = hypot(cosα1, sinα1*sinβ1); // = sqrt(1 - sinα0*sinα0) with less rounding errors.
+ final double k2 = computeSeriesExpansionCoefficients();
+ /*
+ * Compute α₂ from Karney equation 5: sin(α₀) = sin(α₁)⋅cos(β₁) = sin(α₂)⋅cos(β₂)
+ * The cos(α₂) term could be computed from sin(α₂), but Karnay recommends instead
+ * equation 45. An older publication (Karney 2010) went one step further with the
+ * replacement of (cos²β₂ - cos²β₁) by:
+ *
+ * (cosβ₂ - cosβ₁)⋅(cosβ₂ + cosβ₁) if β₁ < -π/4 (Reminder: β₁ is always negative)
+ * (sinβ₁ - sinβ₂)⋅(sinβ₁ + sinβ₂) otherwise.
+ *
+ * Actually we don't need α values directly, but instead cos(α)⋅cos(β).
+ * Note that cos(α₁) can be negative because α₁ ∈ [0…2π].
+ */
+ final double cosα1_cosβ1 = cosα1 * cosβ1;
+ final double cosα2_cosβ2 = sqrt(cosα1_cosβ1*cosα1_cosβ1 +
+ (cosβ1 <= 1/MathFunctions.SQRT_2 ? (cosβ2 - cosβ1)*(cosβ2 + cosβ1)
+ : (sinβ1 - sinβ2)*(sinβ1 + sinβ2)));
+ /*
+ * Karney gives the following formulas:
+ *
+ * σ = atan2(sinβ, cosα⋅cosβ) — spherical arc length.
+ * ω = atan2(sin(α₀)⋅sinσ, cosσ) — spherical longitude.
+ *
+ * We perform the following substitutions where c is an unknown constant:
+ * That unknown coefficient will disaspear in atan2(c⋅y, c⋅x) expressions:
+ *
+ * sin(σ) = c⋅sinβ
+ * cos(σ) = c⋅cosα⋅cosβ
+ */
+ final double ω1, ω2;
+ σ1 = atan2(sinβ1, cosα1_cosβ1); // (Karney 11)
+ σ2 = atan2(sinβ2, cosα2_cosβ2);
+ ω1 = atan2(sinβ1*sinα0, cosα1_cosβ1); // (Karney 12) with above-cited substitutions.
+ ω2 = atan2(sinβ2*sinα0, cosα2_cosβ2);
+ /*
+ * Compute the difference in longitude using the current start point and end point.
+ * If this difference is close enough to the requested accuracy, we are almost done.
+ * Otherwise refine the results. Note that if we detect that accuracy is good enough,
+ * we still complete the computation in order to not waste what we have computed so
+ * far in current iteration.
+ */
+ final double λ2E;
+ λ1E = sphericalToGeodeticLongitude(ω1, σ1);
+ λ2E = sphericalToGeodeticLongitude(ω2, σ2);
+ final double Δλ_error = IEEEremainder(λ2E - λ1E - Δλ, 2*PI);
+ final boolean done = (abs(Δλ_error) <= ACCURACY);
+ if (--nbIter < 0 && !done) {
+ throw new GeodesicException(Resources.format(Resources.Keys.NoConvergence));
+ }
+ /*
+ * Special case for α₁ = π/2 and β₂ = ±β₁ (Karney's equation 47). We replace the β₂ = ±β₁
+ * condition by |β₂| - |β₁| ≈ 0. Assuming tan(θ) ≈ θ for small angles we take the tangent
+ * of above difference and use tan(β₂ - β₁) = (tanβ₂ - tanβ₁)/(1 + tanβ₂⋅tanβ₁) identity.
+ * Note that tanβ₁ ≦ 0 and |tanβ₂| ≦ |tanβ₁| in this method.
+ */
+ final double dΔλ_dα1;
+ if (abs(cosα1) < LATITUDE_THRESHOLD && (-tanβ1 - abs(tanβ2)) < (1 + abs(tanβ1*tanβ2)) * LATITUDE_THRESHOLD) {
+ dΔλ_dα1 = -2 * sqrt(1 - eccentricitySquared * (cosβ1*cosβ1)) / sinβ1;
+ } else {
+ /*
+ * Karney's equation 38 combined with equation 46. The substitutions
+ * for sin(σ) and cos(σ) described above are applied again here.
+ */
+ final double h1 = hypot(sinβ1, cosα1_cosβ1);
+ final double h2 = hypot(sinβ2, cosα2_cosβ2);
+ final double sinσ1 = sinβ1 / h1;
+ final double sinσ2 = sinβ2 / h2;
+ final double cosσ1 = cosα1_cosβ1 / h1;
+ final double cosσ2 = cosα2_cosβ2 / h2;
+ final double J2 = sphericalToEllipsoidalAngle(σ2, true);
+ final double J1 = sphericalToEllipsoidalAngle(σ1, true);
+ final double Δm = (sqrt(1 + k2*(sinσ2*sinσ2))*cosσ1*sinσ2
+ - sqrt(1 + k2*(sinσ1*sinσ1))*sinσ1*cosσ2
+ - cosσ1*cosσ2*(J2 - J1));
+ dΔλ_dα1 = Δm * axisRatio / cosα2_cosβ2;
+ if (STORE_LOCAL_VARIABLES) {
+ store("J(σ₁)", J1);
+ store("J(σ₂)", J2);
+ store("Δm", Δm * semiMinor);
+ }
+ }
+ final double dα1 = Δλ_error / dΔλ_dα1; // Opposite sign of Karney δα₁ term.
+ α1 -= dα1;
+ if (STORE_LOCAL_VARIABLES) { // For comparing values against Karney table 5 and 6.
+ final double I1_σ2;
+ I1_σ1 = sphericalToEllipsoidalAngle(σ1, false); // Required for computation of s₁ in `snapshot()`.
+ I1_σ2 = sphericalToEllipsoidalAngle(σ2, false);
+ snapshot();
+ store("k²", k2);
+ store("σ₁", σ1);
+ store("σ₂", σ2);
+ store("ω₁", ω1);
+ store("ω₂", ω2);
+ store("α₂", atan2(sinα0, cosα2_cosβ2));
+ store("λ₂", λ2E);
+ store("Δλ", λ2E - λ1E);
+ store("δλ", Δλ_error);
+ store("dλ/dα", dΔλ_dα1);
+ store("δσ₁", -dα1);
+ store("α₁", α1);
+ store("s₂", I1_σ2 * semiMinor);
+ store("Δs", (I1_σ2 - I1_σ1) * semiMinor);
+ }
+ if (done) {
+ dλ1 = sinα1;
+ dφ1 = cosα1;
+ dλ2 = sinα0;
+ dφ2 = cosα2_cosβ2;
+ break;
+ }
+ }
+ final double I1_σ2;
+ I1_σ1 = sphericalToEllipsoidalAngle(σ1, false);
+ I1_σ2 = sphericalToEllipsoidalAngle(σ2, false);
+ geodesicDistance = (I1_σ2 - I1_σ1) * semiMinor;
+ /*
+ * Restore the coordinate sign and order to the original configuration.
+ */
+ if (swapPoints) {
+ double t;
+ t = dφ1; dφ1 = dφ2; dφ2 = t;
+ t = dλ1; dλ1 = dλ2; dλ2 = t;
+ }
+ if (inverseLongitudeSigns ^ swapPoints) {
+ dλ1 = -dλ1;
+ dλ2 = -dλ2;
+ }
+ if (inverseLatitudeSigns ^ swapPoints) {
+ dφ1 = -dφ1;
+ dφ2 = -dφ2;
+ }
+ setValid(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE);
+ if (!(swapPoints | inverseLongitudeSigns | inverseLatitudeSigns)) {
+ setValid(COEFFICIENTS_FOR_START_POINT);
+ }
+ }
+
+ /**
+ * Computes the positive root of quartic equation for estimation of α₁ in nearly antipodal case.
+ * Formula is given in appendix B of <a href="https://arxiv.org/pdf/1102.1215.pdf">C.F.F Karney (2011)</a>
+ * given <var>x</var> and <var>y</var> the coordinates on a plane coordinate system centered on the antipodal point:
+ *
+ * {@preformat math
+ * μ⁴ + 2μ³ + (1−x²-y²)μ² − 2y²μ - y² = 0
+ * }
+ *
+ * The results should have only one positive root {@literal (μ > 0)}.
+ *
+ * <div class="section">Condition on <var>y</var> value</div>
+ * This method is indeterminate when <var>y</var> → 0 (it returns {@link Double#NaN}). For values too close to zero,
+ * the result may be non-significative because of rounding errors. For choosing a threshold value for <var>y</var>,
+ * {@code GeodesicsOnEllipsoidTest.Calculator} compares the value computed by this method against the value computed
+ * by {@link org.apache.sis.math.MathFunctions#polynomialRoots(double...)}. If the values differ too much, we presume
+ * that they are mostly noise caused by a <var>y</var> value too low.
+ *
+ * @param x2 the square of <var>x</var>.
+ * @param y2 the square of <var>y</var>.
+ */
+ private static double μ(final double x2, final double y2) {
+ final double r = (x2 + y2 - 1)/6;
+ final double r3 = r*r*r;
+ final double S = x2*y2/4;
+ final double d = S*(S + 2*r3);
+ final double u;
+ if (d < 0) {
+ u = r*(1 + 2*cos(atan2(sqrt(-d), -(S + r3))/3));
+ } else {
+ final double T = cbrt(S + r3 + copySign(sqrt(d), S + r3));
+ u = (T == 0) ? 0 : r + T + r*r/T;
+ }
+ final double v = sqrt(u*u + y2);
+ final double w = (v + u - y2) / (2*v);
+ return (v + u) / (sqrt(v + u + w*w) + w);
+ }
+
+ /**
+ * TODO: not yet implemented.
+ */
+ @Override
+ final void computeRhumbLine() {
+ // Temporary implementation for allowing tests to pass.
+ super.computeRhumbLine();
+ rhumblineLength *= ellipsoid.getSemiMajorAxis() / authalicRadius;
+ if (φ1 != 0 || φ2 != 0) rhumblineLength *= 1.1; // DELETE ME!
+ }
+
+ /**
+ * Takes a snapshot of the current fields in this class. This is used for JUnit tests only. During development phases,
+ * {@link #STORE_LOCAL_VARIABLES} should be {@code true} for allowing {@code GeodesicsOnEllipsoidTest} to verify the
+ * values of a large range of local variables. But when the storage of locale variables is disabled, this method allows
+ * {@code GeodesicsOnEllipsoidTest} to still verify a few variables.
+ */
+ final void snapshot() {
+ store("ε", ε);
+ store("A₁", A1);
+ store("A₂", A2);
+ store("A₃", A3);
+ store("α₀", atan2(sinα0, cosα0));
+ store("I₁(σ₁)", I1_σ1);
+ store("s₁", I1_σ1 * semiMinor);
+ store("λ₁", λ1E);
+ }
+
+ /**
+ * Stores the value of a local variable or a field. This is used in JUnit tests only.
+ *
+ * @param name name of the local variable.
+ * @param value value of the local variable.
+ */
+ void store(String name, double value) {
+ }
+
+ /**
+ * Replaces a computed value by the value given in Karney table. This is used when the result published in Karney
+ * table 4 is used as input for Karney table 5. We need to truncate the value to the same numbers of digits than
+ * Karney, otherwise computation results will differ.
+ */
+ double computedToGiven(double α1) {
+ return α1;
+ }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
index 7687d69..1ebd2aa 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
@@ -37,6 +37,7 @@ import org.apache.sis.measure.AngleFormat;
import org.apache.sis.measure.Latitude;
import org.apache.sis.measure.Units;
import org.apache.sis.geometry.CoordinateFormat;
+import org.apache.sis.referencing.operation.GeodesicException;
import org.apache.sis.internal.referencing.PositionTransformer;
import org.apache.sis.internal.referencing.ReferencingUtilities;
import org.apache.sis.internal.referencing.j2d.ShapeUtilities;
@@ -74,6 +75,7 @@ import static java.lang.Math.*;
* This class is not thread-safe. If geodetic calculations are needed in a multi-threads environment,
* then a distinct instance of {@code GeodeticCalculator} needs to be created for each thread.
*
+ * @author Martin Desruisseaux (Geomatys)
* @version 1.0
*
* @see <a href="https://en.wikipedia.org/wiki/Great-circle_navigation">Great-circle navigation on Wikipedia</a>
@@ -108,7 +110,7 @@ public class GeodeticCalculator {
*
* @see org.apache.sis.referencing.datum.DefaultEllipsoid#getAuthalicRadius()
*/
- private final double radius;
+ final double authalicRadius;
/**
* The (<var>latitude</var>, <var>longitude</var>) coordinates of the start point <strong>in radians</strong>.
@@ -116,7 +118,7 @@ public class GeodeticCalculator {
*
* @see #START_POINT
*/
- private double φ1, λ1;
+ double φ1, λ1;
/**
* The (<var>latitude</var>, <var>longitude</var>) coordinates of the end point <strong>in radians</strong>.
@@ -124,7 +126,7 @@ public class GeodeticCalculator {
*
* @see #END_POINT
*/
- private double φ2, λ2;
+ double φ2, λ2;
/**
* The azimuth at start point and end point as vector components.
@@ -138,11 +140,17 @@ public class GeodeticCalculator {
*
* Those vectors may not be normalized to unitary vectors. This is often not needed because most formulas are
* written in a way that cancel the magnitude. If nevertheless needed, normalization will be applied in formulas.
+ * Those quantities are related to α as below, with α a <em>geographic</em> (not arithmetic) angle:
+ *
+ * <ul>
+ * <li>{@code sinα = dλ / hypot(dφ, dλ)}</li>
+ * <li>{@code cosα = dφ / hypot(dφ, dλ)}</li>
+ * </ul>
*
* @see #STARTING_AZIMUTH
* @see #ENDING_AZIMUTH
*/
- private double dφ1, dλ1, dφ2, dλ2;
+ double dφ1, dλ1, dφ2, dλ2;
/**
* The shortest distance from the starting point ({@link #φ1},{@link #λ1}) to the end point ({@link #φ2},{@link #λ2}).
@@ -151,7 +159,7 @@ public class GeodeticCalculator {
* @see #GEODESIC_DISTANCE
* @see #getDistanceUnit()
*/
- private double geodesicDistance;
+ double geodesicDistance;
/**
* Length of the rhumb line from the starting point ({@link #φ1},{@link #λ1}) to the end point ({@link #φ2},{@link #λ2}).
@@ -160,21 +168,24 @@ public class GeodeticCalculator {
* @see #RHUMBLINE_LENGTH
* @see #getDistanceUnit()
*/
- private double rhumblineLength;
+ double rhumblineLength;
/**
* A bitmask specifying which information are valid. For example if the {@link #END_POINT} bit is not set,
* then {@link #φ2} and {@link #λ2} need to be computed, which implies the computation of ∂φ/∂λ as well.
* If the {@link #GEODESIC_DISTANCE} bit is not set, then {@link #geodesicDistance} needs to be computed,
* which implies recomputation of ∂φ/∂λ as well.
+ *
+ * @see #isInvalid(int)
+ * @see #setValid(int)
*/
private int validity;
/**
* Bitmask specifying which information are valid.
*/
- private static final int START_POINT = 1, END_POINT = 2, STARTING_AZIMUTH = 4, ENDING_AZIMUTH = 8,
- GEODESIC_DISTANCE = 16, RHUMBLINE_LENGTH = 32;
+ static final int START_POINT = 1, END_POINT = 2, STARTING_AZIMUTH = 4, ENDING_AZIMUTH = 8,
+ GEODESIC_DISTANCE = 16, RHUMBLINE_LENGTH = 32, COEFFICIENTS_FOR_START_POINT = 64;
/**
* Constructs a new geodetic calculator expecting coordinates in the supplied CRS.
@@ -195,7 +206,7 @@ public class GeodeticCalculator {
throw new IllegalArgumentException(Errors.format(Errors.Keys.IllegalCRSType_1, ReferencingUtilities.getInterface(crs)));
}
ellipsoid = ReferencingUtilities.getEllipsoid(crs);
- radius = Formulas.getAuthalicRadius(ellipsoid);
+ authalicRadius = Formulas.getAuthalicRadius(ellipsoid);
userToGeodetic = new PositionTransformer(crs, geographic, null);
}
@@ -219,11 +230,18 @@ public class GeodeticCalculator {
/**
* Returns {@code true} if at least one of the properties identified by the given mask is invalid.
*/
- private boolean isInvalid(final int mask) {
+ final boolean isInvalid(final int mask) {
return (validity & mask) != mask;
}
/**
+ * Sets the properties specified by the given bitmask as valid.
+ */
+ final void setValid(final int mask) {
+ validity |= mask;
+ }
+
+ /**
* Returns the Coordinate Reference System (CRS) in which {@code Position}s are represented, unless otherwise specified.
* This is the CRS of all {@link Position} instances returned by methods in this class. This is also the default CRS
* assumed by methods receiving a {@link Position} argument when the given position does not specify its CRS.
@@ -328,7 +346,7 @@ public class GeodeticCalculator {
* and the current azimuth and distance.
*
* @return the destination (end point) represented in the CRS specified at construction time.
- * @throws TransformException if the coordinates can not be transformed to {@linkplain #getPositionCRS() position CRS}.
+ * @throws TransformException if the coordinates can not be computed.
* @throws IllegalStateException if the destination point, azimuth or distance have not been set.
*
* @see #getStartPoint()
@@ -372,8 +390,9 @@ public class GeodeticCalculator {
ArgumentChecks.ensureFinite("longitude", longitude);
φ2 = toRadians(max(Latitude.MIN_VALUE, min(Latitude.MAX_VALUE, latitude)));
λ2 = toRadians(longitude);
- validity |= END_POINT;
- validity &= ~(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE | RHUMBLINE_LENGTH);
+ setValid(END_POINT);
+ validity &= ~(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE | RHUMBLINE_LENGTH | COEFFICIENTS_FOR_START_POINT);
+ // Coefficients for starting point are invalidated because the starting azimuth changed.
}
/**
@@ -385,8 +404,9 @@ public class GeodeticCalculator {
*
* @return the azimuth in degrees from -180° to +180°. 0° is toward North and values are increasing clockwise.
* @throws IllegalStateException if the end point, azimuth or distance have not been set.
+ * @throws GeodesicException if the azimuth can not be computed.
*/
- public double getStartingAzimuth() {
+ public double getStartingAzimuth() throws GeodesicException {
if (isInvalid(STARTING_AZIMUTH)) {
computeDistance();
}
@@ -409,8 +429,8 @@ public class GeodeticCalculator {
azimuth = toRadians(azimuth);
dφ1 = cos(azimuth); // sin(π/2 − θ) = cos(θ)
dλ1 = sin(azimuth); // cos(π/2 − θ) = sin(θ)
- validity |= STARTING_AZIMUTH;
- validity &= ~(END_POINT | ENDING_AZIMUTH | RHUMBLINE_LENGTH);
+ setValid(STARTING_AZIMUTH);
+ validity &= ~(END_POINT | ENDING_AZIMUTH | RHUMBLINE_LENGTH | COEFFICIENTS_FOR_START_POINT);
}
/**
@@ -421,13 +441,14 @@ public class GeodeticCalculator {
*
* @return the azimuth in degrees from -180° to +180°. 0° is toward North and values are increasing clockwise.
* @throws IllegalStateException if the destination point, azimuth or distance have not been set.
+ * @throws GeodesicException if the azimuth can not be computed.
*/
- public double getEndingAzimuth() {
+ public double getEndingAzimuth() throws GeodesicException {
if (isInvalid(ENDING_AZIMUTH)) {
if (isInvalid(END_POINT)) {
computeEndPoint(); // Compute also ending azimuth from start point and distance.
} else {
- computeDistance(); // Compute also ending azimuth from start point and end point.
+ computeDistance(); // Compute also ending azimuth from start point and end point.
}
}
return toDegrees(atan2(dλ2, dφ2));
@@ -441,10 +462,11 @@ public class GeodeticCalculator {
*
* @return the shortest distance in the unit of measurement given by {@link #getDistanceUnit()}.
* @throws IllegalStateException if a point has not been set.
+ * @throws GeodesicException if the distance can not be computed.
*
* @see #getDistanceUnit()
*/
- public double getGeodesicDistance() {
+ public double getGeodesicDistance() throws GeodesicException {
if (isInvalid(GEODESIC_DISTANCE)) {
computeDistance();
}
@@ -463,7 +485,7 @@ public class GeodeticCalculator {
public void setGeodesicDistance(final double distance) {
ArgumentChecks.ensurePositive("distance", distance);
geodesicDistance = distance;
- validity |= GEODESIC_DISTANCE;
+ setValid(GEODESIC_DISTANCE);
validity &= ~(END_POINT | ENDING_AZIMUTH | RHUMBLINE_LENGTH);
}
@@ -509,17 +531,27 @@ public class GeodeticCalculator {
}
/**
- * Computes the length of rhumb line from start point to end point.
+ * Ensures that the start point and end point are set.
+ * This method should be invoked at the beginning of {@link #computeDistance()}.
*
- * @see <a href="https://en.wikipedia.org/wiki/Rhumb_line">Rhumb line on Wikipedia</a>
- *
- * @todo if {@literal Δλ > 180}, must split in two segments.
+ * @throws IllegalStateException if the start point or end point has not been set.
*/
- private void computeRhumbLine() {
+ final void canComputeDistance() {
if (isInvalid(START_POINT | END_POINT)) {
throw new IllegalStateException(Resources.format(
Resources.Keys.StartOrEndPointNotSet_1, Integer.signum(validity & START_POINT)));
}
+ }
+
+ /**
+ * Computes the length of rhumb line from start point to end point.
+ *
+ * @see <a href="https://en.wikipedia.org/wiki/Rhumb_line">Rhumb line on Wikipedia</a>
+ *
+ * @todo if {@literal Δλ > 180}, must split in two segments.
+ */
+ void computeRhumbLine() {
+ canComputeDistance();
final double Δλ = λ2 - λ1;
final double Δφ = φ2 - φ1;
final double factor;
@@ -539,12 +571,12 @@ public class GeodeticCalculator {
final double α = atan(Δλ / ΔG);
factor = Δφ / cos(α);
}
- rhumblineLength = radius * abs(factor);
- validity |= RHUMBLINE_LENGTH;
+ rhumblineLength = authalicRadius * abs(factor);
+ setValid(RHUMBLINE_LENGTH);
}
/**
- * Computes the geodetic distance and azimuths from the start point and end point.
+ * Computes the geodesic distance and azimuths from the start point and end point.
* This method should be invoked if the distance or an azimuth is requested while
* {@link #STARTING_AZIMUTH}, {@link #ENDING_AZIMUTH} or {@link #GEODESIC_DISTANCE}
* validity flag is not set.
@@ -556,12 +588,10 @@ public class GeodeticCalculator {
* </ul>
*
* @throws IllegalStateException if the distance or azimuth has not been set.
+ * @throws GeodesicException if an azimuth or the distance can not be computed.
*/
- private void computeDistance() {
- if (isInvalid(START_POINT | END_POINT)) {
- throw new IllegalStateException(Resources.format(
- Resources.Keys.StartOrEndPointNotSet_1, Integer.signum(validity & START_POINT)));
- }
+ void computeDistance() throws GeodesicException {
+ canComputeDistance();
final double Δλ = λ2 - λ1; // No need to reduce to −π … +π range.
final double sinΔλ = sin(Δλ);
final double cosΔλ = cos(Δλ);
@@ -581,9 +611,23 @@ public class GeodeticCalculator {
* Δσ = atan2(…) computes the same value but with better accuracy.
*/
double Δσ = sinφ1*sinφ2 + cosφ1*cosφ2*cosΔλ; // Actually Δσ = acos(…).
- Δσ = atan2(hypot(dλ1, dφ1), Δσ);
- geodesicDistance = radius * Δσ;
- validity |= (STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE);
+ Δσ = atan2(hypot(dλ1, dφ1), Δσ); // Δσ ∈ [0…π].
+ geodesicDistance = authalicRadius * Δσ;
+ setValid(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE);
+ }
+
+ /**
+ * Ensures that the start point, starting azimuth and geodesic distance are set.
+ * This method should be invoked at the beginning of {@link #computeEndPoint()}.
+ *
+ * @throws IllegalStateException if the start point, azimuth or distance has not been set.
+ */
+ final void canComputeEndPoint() {
+ if (isInvalid(START_POINT | STARTING_AZIMUTH | GEODESIC_DISTANCE)) {
+ throw new IllegalStateException(isInvalid(START_POINT)
+ ? Resources.format(Resources.Keys.StartOrEndPointNotSet_1, 0)
+ : Resources.format(Resources.Keys.AzimuthAndDistanceNotSet));
+ }
}
/**
@@ -594,21 +638,18 @@ public class GeodeticCalculator {
* <p>The default implementation computes {@link #φ2}, {@link #λ2} and ∂φ/∂λ derivatives using
* spherical formulas. Subclasses should override if they can provide ellipsoidal formulas.</p>
*
- * @throws IllegalStateException if the azimuth and the distance have not been set.
+ * @throws IllegalStateException if the start point, azimuth or distance has not been set.
+ * @throws GeodesicException if the end point or azimuth can not be computed.
*/
- void computeEndPoint() {
- if (isInvalid(START_POINT | STARTING_AZIMUTH | GEODESIC_DISTANCE)) {
- throw new IllegalStateException(isInvalid(START_POINT)
- ? Resources.format(Resources.Keys.StartOrEndPointNotSet_1, 0)
- : Resources.format(Resources.Keys.AzimuthAndDistanceNotSet));
- }
+ void computeEndPoint() throws GeodesicException {
+ canComputeEndPoint(); // May throw IllegalStateException.
final double vm = hypot(dφ1, dλ1);
- final double Δσ = geodesicDistance / radius;
+ final double Δσ = geodesicDistance / authalicRadius;
final double sinΔσ = sin(Δσ);
final double cosΔσ = cos(Δσ);
final double sinφ1 = sin(φ1);
final double cosφ1 = cos(φ1);
- final double sinα1 = dλ1 / vm; // α₁ is the azimuth at starting point.
+ final double sinα1 = dλ1 / vm; // α₁ is the azimuth at starting point as a geographic (not arithmetic) angle.
final double cosα1 = dφ1 / vm;
final double sinΔσ_cosα1 = sinΔσ * cosα1;
final double Δλy = sinΔσ * sinα1;
@@ -620,7 +661,7 @@ public class GeodeticCalculator {
λ2 = IEEEremainder(λ1 + Δλ, 2*PI);
dφ2 = cosΔσ*cosα1 - sinφ1/cosφ1 * sinΔσ;
dλ2 = sinα1;
- validity |= END_POINT;
+ setValid(END_POINT | ENDING_AZIMUTH);
}
/**
@@ -630,14 +671,19 @@ public class GeodeticCalculator {
* some of them will need to be specified again.
*
* @see #setStartGeographicPoint(double, double)
+ * @throws GeodesicException if the end point or azimuth can not be computed.
*/
- public void moveToEndPoint() {
+ public void moveToEndPoint() throws GeodesicException {
if (isInvalid(END_POINT)) {
computeEndPoint();
}
- φ1 = φ2; dφ1 = dφ2;
- λ1 = λ2; dλ1 = dλ2;
- validity = START_POINT | STARTING_AZIMUTH;
+ φ1 = φ2; dφ1 = dφ2;
+ λ1 = λ2; dλ1 = dλ2;
+ /*
+ * Following assumes that ENDING_AZIMUTH >>> 1 == STARTING_AZIMUTH.
+ * This is verified by GeodeticCalculatorTest.testMoveToEndPoint().
+ */
+ validity = ((validity & ENDING_AZIMUTH) >>> 1) | START_POINT;
}
/**
@@ -756,7 +802,7 @@ public class GeodeticCalculator {
PathBuilder(final double εx) {
super(ReferencingUtilities.getDimension(userToGeodetic.defaultCRS));
dφf = dφ2; dλf = dλ2; φf = φ2; λf = λ2;
- tolerance = toDegrees(εx / radius);
+ tolerance = toDegrees(εx / authalicRadius);
distance = geodesicDistance;
length = rhumblineLength;
flags = validity;
@@ -780,7 +826,7 @@ public class GeodeticCalculator {
λ2 = λf; dλ2 = dλf;
} else {
geodesicDistance = distance * t;
- validity |= GEODESIC_DISTANCE;
+ setValid(GEODESIC_DISTANCE);
computeEndPoint();
}
evaluateAtEndPoint();
@@ -863,7 +909,8 @@ public class GeodeticCalculator {
final double α1 = t * (2*PI);
dλ1 = sin(α1);
dφ1 = cos(α1);
- validity |= STARTING_AZIMUTH;
+ setValid(STARTING_AZIMUTH);
+ validity &= ~COEFFICIENTS_FOR_START_POINT;
computeEndPoint();
evaluateAtEndPoint();
if (depth <= 1) {
@@ -964,7 +1011,7 @@ public class GeodeticCalculator {
resources.appendLabel(Vocabulary.Keys.GeodesicDistance, buffer);
buffer.append(' ').append(nf.format(distance))
.append(' ').append(unit).append(lineSeparator);
- } catch (IllegalStateException e) {
+ } catch (IllegalStateException | GeodesicException e) {
// Ignore.
}
} catch (IOException e) {
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/GeodesicException.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/GeodesicException.java
new file mode 100644
index 0000000..ca6ce6a
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/GeodesicException.java
@@ -0,0 +1,56 @@
+/*
+ * 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.operation;
+
+import org.opengis.referencing.operation.TransformException;
+
+
+/**
+ * Thrown when an error occurred while computing geodesic between two points.
+ *
+ * <div class="note"><b>API note:</b>
+ * defined as a sub-type of {@link TransformException} because some kind of coordinate operations are involved in
+ * geodesics calculation (e.g. transformation of latitudes and longitudes to coordinates on an auxiliary sphere).
+ * The starting and ending points can also be given in a CRS that require additional coordinate operations.</div>
+ *
+ * @author Martin Desruisseaux (Geomatys)
+ * @version 1.0
+ * @since 1.0
+ * @module
+ */
+public class GeodesicException extends TransformException {
+ /**
+ * For cross-version compatibility.
+ */
+ private static final long serialVersionUID = -7844524593329236175L;
+
+ /**
+ * Constructs a new exception with no message.
+ */
+ public GeodesicException() {
+ super();
+ }
+
+ /**
+ * Constructs a new exception with the specified detail message.
+ *
+ * @param message the detail message, or {@code null} if none.
+ */
+ public GeodesicException(final String message) {
+ super(message);
+ }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java
index 5f96038..e48903b 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java
@@ -114,10 +114,10 @@ public class Mollweide extends NormalizedProjection {
*/
if (abs(sinφ) != 1) {
final double πsinφ = PI * sinφ;
- int nbIte = MAXIMUM_ITERATIONS;
+ int nbIter = MAXIMUM_ITERATIONS;
double Δθ;
do {
- if (--nbIte < 0) {
+ if (--nbIter < 0) {
throw new ProjectionException(Resources.format(Resources.Keys.NoConvergence));
}
Δθ = (θp + sin(θp) - πsinφ) / (1 + cos(θp));
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java
new file mode 100644
index 0000000..24f4145
--- /dev/null
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java
@@ -0,0 +1,506 @@
+/*
+ * 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 java.util.Map;
+import java.util.HashMap;
+import java.util.Arrays;
+import org.opengis.geometry.DirectPosition;
+import org.opengis.referencing.crs.GeographicCRS;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.referencing.operation.GeodesicException;
+import org.apache.sis.internal.referencing.Formulas;
+import org.apache.sis.referencing.crs.HardCodedCRS;
+import org.apache.sis.math.MathFunctions;
+import org.apache.sis.measure.Units;
+import org.apache.sis.test.DependsOn;
+import org.apache.sis.test.DependsOnMethod;
+import org.junit.Test;
+
+import static java.lang.StrictMath.*;
+import static org.junit.Assert.*;
+
+
+/**
+ * Tests {@link GeodesicsOnEllipsoid}. Some of the tests in this class use data published by
+ * <a href="https://link.springer.com/content/pdf/10.1007%2Fs00190-012-0578-z.pdf">Algorithms
+ * for geodesics from Charles F. F. Karney (SRI International)</a>.
+ *
+ * @author Matthieu Bastianelli (Geomatys)
+ * @author Martin Desruisseaux (Geomatys)
+ * @version 1.0
+ * @since 1.0
+ * @module
+ */
+@DependsOn(GeodeticCalculatorTest.class)
+public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorTest {
+ /**
+ * The instance to be tested.
+ */
+ private GeodesicsOnEllipsoid testedEarth;
+
+ /**
+ * Values of local variables in {@link GeodesicsOnEllipsoid} methods. If values for the same key are added
+ * many times, the new values are appended after the current ones. This happen for example during iterations.
+ *
+ * @see GeodesicsOnEllipsoid#STORE_LOCAL_VARIABLES
+ */
+ private Map<String,double[]> localVariables;
+
+ /**
+ * {@code true} if {@link GeodesicsOnEllipsoid#store(String, double)} shall verify consistency instead than
+ * storing local variables.
+ */
+ private boolean verifyConsistency;
+
+ /**
+ * Creates a new test case.
+ */
+ public GeodesicsOnEllipsoidTest() {
+ }
+
+ /**
+ * Returns the calculator to use for testing purpose.
+ *
+ * @param normalized whether to force (longitude, latitude) axis order.
+ */
+ @Override
+ GeodeticCalculator create(final boolean normalized) {
+ testedEarth = new Calculator(normalized ? HardCodedCRS.WGS84 : HardCodedCRS.WGS84_φλ);
+ return testedEarth;
+ }
+
+ /**
+ * Creates a calculator which will store locale variables in the {@link #localVariables} map.
+ * This is used for the test cases that compare intermediate calculations with values provided
+ * in tables of Karney (2013) <cite>Algorithms for geodesics</cite> publication.
+ */
+ private void createTracked() {
+ localVariables = new HashMap<>();
+ testedEarth = new Calculator(HardCodedCRS.WGS84) {
+ @Override double computedToGiven(final double α1) {
+ return (abs(TRUNCATED_α1 - toDegrees(α1)) < 1E-3) ? toRadians(TRUNCATED_α1) : α1;
+ }
+
+ @Override void store(final String name, final double value) {
+ super.store(name, value);
+ if (verifyConsistency) {
+ assertValueEquals(name, 0, value, STRICT, false);
+ } else {
+ double[] values = localVariables.putIfAbsent(name, new double[] {value});
+ if (values != null) {
+ final int i = values.length;
+ values = Arrays.copyOf(values, i+1);
+ values[i] = value;
+ localVariables.put(name, values);
+ }
+ }
+ }
+ };
+ }
+
+ /**
+ * The {@link GeodesicsOnEllipsoid} implementation to use for tests. This implementation compares the
+ * quartic root computed by {@link GeodesicsOnEllipsoid#μ(double, double)} against the values computed
+ * by {@link MathFunctions#polynomialRoots(double...)}.
+ */
+ private static class Calculator extends GeodesicsOnEllipsoid {
+ /** Values needed for computation of μ. */
+ private double x, y;
+
+ /** Creates a new calculator for the given coordinate reference system. */
+ Calculator(final GeographicCRS crs) {
+ super(crs);
+ }
+
+ /** Invoked when {@link GeodesicsOnEllipsoid} computed an intermediate value. */
+ @Override void store(final String name, final double value) {
+ if (name.length() == 1) {
+ switch (name.charAt(0)) {
+ case 'x': x = value; break;
+ case 'y': y = value; break;
+ case 'μ': {
+ double μ = 0;
+ final double y2 = y*y;
+ for (final double r : MathFunctions.polynomialRoots(-y2, -2*y2, (1 - x*x)-y2, 2, 1)) {
+ if (r > μ) μ = r;
+ }
+ if (μ > 0) {
+ /*
+ * This assertion fails if y² is too close to zero. We that as a criterion for choosing
+ * the threshold value that determine when to use Karney (2013) equation 57 instead.
+ * That threshold is in `GeodesicOnEllipsoid.computeDistance()` method.
+ */
+ assertEquals("μ(x²,y²)", μ, value, abs(μ) * 1E-8);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /**
+ * Asserts that variable of the given name is equal to the given value. This is used for comparing an
+ * intermediate value computed by Apache SIS against the value published in a table of Karney (2013)
+ * <cite>Algorithms for geodesics</cite>.
+ *
+ * @param name name of the variable to verify.
+ * @param index if the same variable has been set many times, its index with numbering starting at 1. Otherwise 0.
+ * @param expected the expected value.
+ * @param tolerance tolerance threshold for comparison.
+ * @param angular whether the stored value needs to be converted from radians to degrees.
+ */
+ private void assertValueEquals(final String name, int index, final double expected, final double tolerance, final boolean angular) {
+ final double[] values = localVariables.get(name);
+ if (values != null) {
+ if (index > 0) {
+ index--;
+ } else if (values.length != 1) {
+ fail("Expected exactly one value for " + name + " but got " + Arrays.toString(values));
+ }
+ double value = values[index];
+ if (angular) value = toDegrees(value);
+ assertEquals(name, expected, value, tolerance);
+ } else if (GeodesicsOnEllipsoid.STORE_LOCAL_VARIABLES) {
+ fail("Missing value: " + name);
+ }
+ }
+
+ /**
+ * Takes a snapshot of some intermediate calculations in {@link GeodesicsOnEllipsoid}.
+ * If {@link GeodesicsOnEllipsoid#STORE_LOCAL_VARIABLES} is {@code true}, then we already
+ * have snapshots of those variables and this method will instead verify consistency.
+ */
+ private void snapshot() {
+ verifyConsistency = true;
+ testedEarth.snapshot();
+ }
+
+ /**
+ * Verifies that the constant parameters set by constructor are equal to the ones used by Karney
+ * for his example. Those values do not depend on the point being tested.
+ *
+ * <p><b>Source:</b>Karney (2013) table 1.</p>
+ */
+ private void verifyParametersForWGS84() {
+ assertEquals("1/f", 298.257223563, testedEarth.ellipsoid.getInverseFlattening(), 1E-9);
+ assertEquals("a", 6378137, testedEarth.ellipsoid.getSemiMajorAxis(), STRICT);
+ assertEquals("b", 6356752.314245, testedEarth.ellipsoid.getSemiMinorAxis(), 1E-6);
+ assertEquals("c", 6371007.180918, testedEarth.authalicRadius, 1E-6);
+ assertEquals("ℯ²", 0.00669437999014132, testedEarth.eccentricitySquared, 1E-17);
+ assertEquals("ℯ′²", 0.00673949674227643, testedEarth.secondEccentricitySquared, 1E-17);
+ assertEquals("n", 0.00167922038638370, testedEarth.thirdFlattening, 1E-16);
+ assertEquals("unit", Units.METRE, testedEarth.getDistanceUnit());
+ }
+
+ /**
+ * Tests solving the direct geodesic problem by a call to {@link GeodesicsOnEllipsoid#computeEndPoint()}.
+ * The data used by this test are provided by the example published in Karney (2013) table 2.
+ * Input values are:
+ *
+ * <ul>
+ * <li>Starting longitude: λ₁ = any</li>
+ * <li>Starting latitude: φ₁ = 40°</li>
+ * <li>Starting azimuth: α₁ = 30°</li>
+ * <li>Geodesic distance: s₁₂ = 10 000 000 m</li>
+ * </ul>
+ *
+ * About 20 intermediate values are published in Karney table 2 and copied in this method body.
+ * Those values are verified if {@link GeodesicsOnEllipsoid#STORE_LOCAL_VARIABLES} is {@code true}.
+ *
+ * <p><b>Source:</b> Karney (2013), <u>Algorithms for geodesics</u> table 2.</p>
+ *
+ * @throws TransformException should never happen since this method does not apply map projections.
+ */
+ @Test
+ public void testComputeEndPoint() throws TransformException {
+ createTracked();
+ verifyParametersForWGS84();
+ testedEarth.setStartGeographicPoint(40, 10);
+ testedEarth.setStartingAzimuth(30);
+ testedEarth.setGeodesicDistance(10000000);
+ assertEquals("λ₁", toDegrees(testedEarth.λ1), 10.0, Formulas.ANGULAR_TOLERANCE);
+ assertEquals("φ₁", toDegrees(testedEarth.φ1), 40.0, Formulas.ANGULAR_TOLERANCE);
+ assertEquals("∆s", testedEarth.geodesicDistance, 10000000, STRICT);
+ /*
+ * The following method invocation causes calculation of all intermediate values.
+ * Some intermediate values are verified by `assertValueEquals` statements below.
+ * If some values are wrong, it is easier to check correctness with step-by-step
+ * debugging of `computeDistance()` method while keeping table 3 of Karney aside.
+ */
+ final DirectPosition endPoint = testedEarth.getEndPoint();
+ snapshot();
+ /*
+ * Start point on auxiliary sphere:
+ *
+ * β₁ — reduced latitude of starting point.
+ * α₀ — azimuth at equator in the direction given by α₁ azimuth at the β₁ reduced latitude.
+ * σ₁ — arc length from equatorial point E to starting point on auxiliary sphere.
+ * ω₁ — spherical longitude of the starting point (from equatorial point E) on the auxiliary sphere.
+ * s₁ — distance from equatorial point E to starting point along the geodesic.
+ */
+ assertValueEquals("β₁", 0, 39.90527714601, 1E-11, true);
+ assertValueEquals("α₀", 0, 22.55394020262, 1E-11, true);
+ assertValueEquals("σ₁", 0, 43.99915364500, 1E-11, true);
+ assertValueEquals("ω₁", 0, 20.32371827837, 1E-11, true);
+ /*
+ * End point on auxiliary sphere.
+ *
+ * σ₂ — arc lentgth from equator to the ending point σ₂ on auxiliary sphere.
+ * β₂ — reduced latitude of ending point.
+ * ω₂ — spherical longitude of the ending point (from equatorial point E) on the auxiliary sphere.
+ */
+ assertValueEquals("k²", 0, 0.00574802962857, 1E-14, false);
+ assertValueEquals("ε", 0, 0.00143289220416, 1E-14, false);
+ assertValueEquals("A₁", 0, 1.00143546236207, 1E-14, false);
+ assertValueEquals("I₁(σ₁)", 0, 0.76831538886412, 1E-14, false);
+ assertValueEquals("s₁", 0, 4883990.626232, 1E-6, false);
+ assertValueEquals("s₂", 0, 14883990.626232, 1E-6, false);
+ assertValueEquals("τ₂", 0, 133.96266050208, 1E-11, true);
+ assertValueEquals("σ₂", 0, 133.92164083038, 1E-11, true);
+ assertValueEquals("α₂", 0, 149.09016931807, 1E-11, true);
+ assertValueEquals("β₂", 0, 41.69771809250, 1E-11, true);
+ assertValueEquals("ω₂", 0, 158.28412147112, 1E-11, true);
+ /*
+ * Conversion of end point coordinates to geodetic longitude.
+ *
+ * λ₁ — longitudes of the starting point, measured from equatorial point E.
+ * λ₂ — longitudes of the ending point, measured from Equatorial point E.
+ */
+ assertValueEquals("A₃", 0, 0.99928424306, 1E-11, false);
+ assertValueEquals("I₃(σ)", 1, 0.76773786069, 1E-11, false);
+ assertValueEquals("I₃(σ)", 2, 2.33534322170, 1E-11, false);
+ assertValueEquals("λ₁", 0, 20.26715038016, 1E-11, true);
+ assertValueEquals("λ₂", 0, 158.11205042393, 1E-11, true);
+ assertValueEquals("Δλ", 0, 137.84490004377, 1E-11, true);
+ /*
+ * Final result:
+ *
+ * φ₂ — end point latitude.
+ * λ₂ — end point longitude. Shifted by 10° compared to Karney because of λ₁ = 10°.
+ * α₂ — azimuth at end point.
+ */
+ assertEquals("φ₂", 41.79331020506, endPoint.getOrdinate(1), 1E-11);
+ assertEquals("λ₂", 147.84490004377, endPoint.getOrdinate(0), 1E-11);
+ assertEquals("α₂", 149.09016931807, testedEarth.getEndingAzimuth(), 1E-11);
+ }
+
+ /**
+ * Tests solving the inverse geodesic problem by a call to {@link GeodesicsOnEllipsoid#computeDistance()}
+ * for a short distance. The data used by this test are provided by the example published in Karney (2013)
+ * table 3. Input values are:
+ *
+ * <ul>
+ * <li>Starting longitude: λ₁ = any</li>
+ * <li>Starting latitude: φ₁ = -30.12345°</li>
+ * <li>Ending longitude: λ₂ = λ₁ + 0.00005°</li>
+ * <li>Ending latitude: φ₂ = -30.12344°</li>
+ * </ul>
+ *
+ * About 5 intermediate values are published in Karney table 3 and copied in this method body.
+ * Those values are verified if {@link GeodesicsOnEllipsoid#STORE_LOCAL_VARIABLES} is {@code true}.
+ *
+ * <p><b>Source:</b> Karney (2013), <u>Algorithms for geodesics</u> table 3.</p>
+ *
+ * @throws GeodesicException if a geodesic can not be computed.
+ */
+ @Test
+ @DependsOnMethod("testComputeEndPoint")
+ public void testComputeShortDistance() throws GeodesicException {
+ createTracked();
+ verifyParametersForWGS84();
+ testedEarth.setStartGeographicPoint(-30.12345, 2);
+ testedEarth.setEndGeographicPoint (-30.12344, 2.00005);
+ /*
+ * The following method invocation causes calculation of all intermediate values.
+ * Some intermediate values are verified by `assertValueEquals` statements below.
+ * If some values are wrong, it is easier to check correctness with step-by-step
+ * debugging of `computeDistance()` method while keeping table 3 of Karney aside.
+ *
+ * β₁ — reduced latitude of starting point.
+ * β₂ — reduced latitude of ending point.
+ */
+ final double distance = testedEarth.getGeodesicDistance();
+ assertValueEquals("β₁", 0, -30.03999083821, 1E-11, true);
+ assertValueEquals("β₂", 0, -30.03998085491, 1E-11, true);
+ assertValueEquals("ωb", 0, 0.99748847744, 1E-11, false);
+ assertValueEquals("Δω", 0, 0.00005012589, 1E-11, true);
+ assertValueEquals("Δσ", 0, 0.00004452641, 1E-11, true);
+ assertValueEquals("α₁", 1, 77.04353354237, 1E-8, true);
+ assertValueEquals("α₂", 1, 77.04350844913, 1E-8, true);
+ /*
+ * Final result. Values are slightly different than the values published
+ * in Karney table 3 because GeodeticCalculator has done an iteration.
+ */
+ assertEquals("α₁", 77.04353354237, testedEarth.getStartingAzimuth(), 1E-8); // Last 3 digits differ.
+ assertEquals("α₂", 77.04350844913, testedEarth.getEndingAzimuth(), 1E-8);
+ assertEquals("Δs", 4.944208, distance, 1E-6);
+ }
+
+ /**
+ * Result of Karney table 4, used as input in Karney table 5. We need to truncated that intermediate result
+ * to the same number of digits than Karney in order to get the numbers published in table 5 and 6.
+ * This value is used in {@link #testComputeNearlyAntipodal()} only.
+ */
+ private static final double TRUNCATED_α1 = 161.914;
+
+ /**
+ * Tests solving the inverse geodesic problem by a call to {@link GeodesicsOnEllipsoid#computeDistance()} for nearly
+ * antipodal points. The data used by this test are provided by the example published in Karney (2013) tables 4 to 6.
+ * Input values are:
+ *
+ * <ul>
+ * <li>Starting longitude: λ₁ = any</li>
+ * <li>Starting latitude: φ₁ = -30°</li>
+ * <li>Ending longitude: λ₂ = λ₁ + 179.8°</li>
+ * <li>Ending latitude: φ₂ = 29.9°</li>
+ * </ul>
+ *
+ * Some intermediate values are published in Karney tables 4 to 6 and copied in this method body.
+ * Those values are verified if {@link GeodesicsOnEllipsoid#STORE_LOCAL_VARIABLES} is {@code true}.
+ *
+ * <p><b>Source:</b> Karney (2013), <u>Algorithms for geodesics</u> tables 4, 5 and 6.</p>
+ *
+ * @throws GeodesicException if a geodesic can not be computed.
+ */
+ @Test
+ @DependsOnMethod("testComputeShortDistance")
+ public void testComputeNearlyAntipodal() throws GeodesicException {
+ createTracked();
+ verifyParametersForWGS84();
+ testedEarth.setStartGeographicPoint(-30, 0);
+ testedEarth.setEndGeographicPoint(29.9, 179.8);
+ /*
+ * The following method invocation causes calculation of all intermediate values.
+ * Values β₁ and β₂ are keep constant during all iterations.
+ * Other values are given in Karney table 4.
+ */
+ final double distance = testedEarth.getGeodesicDistance();
+ assertValueEquals("β₁", 0, -29.91674771324, 1E-11, true);
+ assertValueEquals("β₂", 0, 29.81691642189, 1E-11, true);
+ assertValueEquals("x", 0, -0.382344, 1E-6, false);
+ assertValueEquals("y", 0, -0.220189, 1E-6, false);
+ assertValueEquals("μ", 0, 0.231633, 1E-6, false);
+ assertValueEquals("α₁", 1, TRUNCATED_α1, 1E-3, true); // Initial value before iteration.
+ /*
+ * Following values are updated during iterations. Note that in order to get the same values than the ones
+ * published in Karney table 5, we need to truncate the α₁ initial value to the same number of digits than
+ * Karney. This is done automatically if GeodesicsOnEllipsoid.STORE_LOCAL_VARIABLES is true.
+ */
+ assertValueEquals("α₀", 1, 15.60939746414, 1E-11, true);
+ assertValueEquals("σ₁", 1, -148.81253566596, 1E-11, true);
+ assertValueEquals("ω₁", 1, -170.74896696128, 1E-11, true);
+ assertValueEquals("α₂", 1, 18.06728796231, 1E-11, true);
+ assertValueEquals("σ₂", 1, 31.08244976895, 1E-11, true);
+ assertValueEquals("ω₂", 1, 9.21345761110, 1E-11, true);
+ assertValueEquals("k²", 1, 0.00625153791662, 1E-14, false);
+ assertValueEquals("ε", 1, 0.00155801826780, 1E-14, false);
+ assertValueEquals("λ₁", 1, -170.61483552458, 1E-11, true);
+ assertValueEquals("λ₂", 1, 9.18542009839, 1E-11, true);
+ assertValueEquals("Δλ", 1, 179.80025562297, 1E-11, true);
+ assertValueEquals("δλ", 1, 0.00025562297, 1E-11, true);
+ assertValueEquals("J(σ₁)", 1, -0.00948040927640, 1E-14, false);
+ assertValueEquals("J(σ₂)", 1, 0.00031349128630, 1E-14, false);
+ assertValueEquals("Δm", 1, 57288.000110, 1E-6, false);
+ assertValueEquals("dλ/dα", 1, 0.01088931716115, 1E-14, false);
+ assertValueEquals("δσ₁", 1, -0.02347465519, 1E-11, true);
+ assertValueEquals("α₁", 2, 161.89052534481, 1E-11, true);
+ /*
+ * After second iteration.
+ */
+ assertValueEquals("δλ", 2, 0.00000000663, 1E-11, true);
+ assertValueEquals("α₁", 3, 161.89052473633, 1E-11, true);
+ /*
+ * After third iteration.
+ */
+ assertValueEquals("α₀", 3, 15.62947966537, 1E-11, true);
+ assertValueEquals("σ₁", 3, -148.80913691776, 1E-11, true);
+ assertValueEquals("ω₁", 3, -170.73634378066, 1E-11, true);
+ assertValueEquals("α₂", 3, 18.09073724574, 1E-11, true);
+ assertValueEquals("σ₂", 3, 31.08583447040, 1E-11, true);
+ assertValueEquals("ω₂", 3, 9.22602862110, 1E-11, true);
+ assertValueEquals("s₁", 3, -16539979.064227, 1E-6, false);
+ assertValueEquals("s₂", 3, 3449853.763383, 1E-6, false);
+ assertValueEquals("Δs", 3, 19989832.827610, 1E-6, false);
+ assertValueEquals("λ₁", 3, -170.60204712148, 1E-11, true);
+ assertValueEquals("λ₂", 3, 9.19795287852, 1E-11, true);
+ assertValueEquals("Δλ", 3, 179.80000000000, 1E-11, true);
+ /*
+ * Final result:
+ */
+ assertEquals("α₁", 161.89052473633, testedEarth.getStartingAzimuth(), 1E-11);
+ assertEquals("α₂", 18.09073724574, testedEarth.getEndingAzimuth(), 1E-11);
+ assertEquals("Δs", 19989832.827610, distance, 1E-6);
+ }
+
+ /**
+ * Same test than the one defined in parent class, but with expected results modified for ellipsoidal formulas.
+ * Input points are from <a href="https://en.wikipedia.org/wiki/Great-circle_navigation#Example">Wikipedia</a>.
+ * Outputs were computed with GeographicLib.
+ *
+ * @throws TransformException if an error occurred while transforming coordinates.
+ */
+ @Test
+ @Override
+ @DependsOnMethod({"testComputeEndPoint", "testComputeNearlyAntipodal"})
+ public void testGeodesicDistanceAndAzimuths() throws TransformException {
+ final GeodeticCalculator c = create(false);
+ c.setStartGeographicPoint(-33.0, -71.6); // Valparaíso
+ c.setEndGeographicPoint ( 31.4, 121.8); // Shanghai
+ /*
+ * α₁ α₂ distance
+ * Spherical: -94.41° -78.42° 18743 km
+ * Ellipsoidal: -94.82° -78.29° 18752 km
+ */
+ assertEquals(Units.METRE, c.getDistanceUnit());
+ assertEquals("α₁", -94.82071749, c.getStartingAzimuth(), 1E-8);
+ assertEquals("α₂", -78.28609385, c.getEndingAzimuth(), 1E-8);
+ assertEquals("distance", 18752.493521, c.getGeodesicDistance() / 1000, 1E-6);
+ assertPositionEquals(31.4, 121.8, c.getEndPoint(), 1E-12); // Should be the specified value.
+ /*
+ * Keep start point unchanged, but set above azimuth and distance.
+ * Verify that we get the Shanghai coordinates.
+ */
+ c.setStartingAzimuth(-94.82);
+ c.setGeodesicDistance(18752494);
+ assertEquals("α₁", -94.82, c.getStartingAzimuth(), 1E-12); // Should be the specified value.
+ assertEquals("α₂", -78.28678389, c.getEndingAzimuth(), 1E-8);
+ assertEquals("Δs", 18752.494, c.getGeodesicDistance() / 1000, STRICT); // Should be the specified value.
+ assertPositionEquals(31.4, 121.8, c.getEndPoint(), 0.0002);
+ }
+
+ /**
+ * Tells whether failure to compute geodesic for the given data should cause the test case to fail.
+ * This is invoked by {@link #compareAgainstDataset()}.
+ *
+ * @param expected a row from the {@code $SIS_DATA/Tests/GeodTest.dat} file.
+ * @return whether the JUnit test should fail.
+ */
+ @Override
+ boolean isFailure(final double[] expected) {
+ final double φ1 = expected[COLUMN_φ1];
+ final double φ2 = expected[COLUMN_φ2];
+ final double Δλ = expected[COLUMN_λ2] - expected[COLUMN_λ1];
+ if (Δλ > 90 && max(abs(φ1), abs(φ2)) < 2E-4) {
+ return false; // Ignore equatorial case.
+ }
+ if (Δλ > 179 && abs(φ1 + φ2) < 0.002) {
+ return false; // Ignore antipodal case.
+ }
+ return super.isFailure(expected);
+ }
+}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
index b969d95..5b58c76 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
@@ -28,6 +28,7 @@ import java.io.LineNumberReader;
import org.opengis.geometry.DirectPosition;
import org.opengis.referencing.cs.AxisDirection;
import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.referencing.operation.GeodesicException;
import org.apache.sis.internal.referencing.j2d.ShapeUtilitiesExt;
import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.referencing.crs.HardCodedCRS;
@@ -58,11 +59,26 @@ import static org.opengis.test.Assert.*;
* <li>Charles Karney's <a href="https://geographiclib.sourceforge.io/">GeographicLib</a> implementation.</li>
* </ul>
*
+ * This base class tests calculator using spherical formulas.
+ * Subclass executes the same test but using ellipsoidal formulas.
+ *
* @version 1.0
- * @since 1.0
+ * @since 1.0
* @module
*/
-public final strictfp class GeodeticCalculatorTest extends TestCase {
+public strictfp class GeodeticCalculatorTest extends TestCase {
+ /**
+ * A threshold in radians for considering two points as antipodal.
+ * A value of 1E-4 radians correspond to approximately 637 metres at equator.
+ */
+ private static final double ANTIPODAL_THRESHOLD = 1E-4;
+
+ /**
+ * Creates a new test case.
+ */
+ public GeodeticCalculatorTest() {
+ }
+
/**
* Verifies that the given point is equals to the given latitude and longitude.
*
@@ -71,7 +87,7 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
* @param p the actual position to verify.
* @param ε the tolerance threshold.
*/
- private static void assertPositionEquals(final double φ, final double λ, final DirectPosition p, final double ε) {
+ static void assertPositionEquals(final double φ, final double λ, final DirectPosition p, final double ε) {
assertEquals("φ", φ, p.getOrdinate(0), ε);
assertEquals("λ", λ, p.getOrdinate(1), ε);
}
@@ -90,21 +106,31 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
}
/**
- * Returns the calculator to use for testing purpose. This classes uses a calculator for a sphere.
+ * Returns the calculator to use for testing purpose. Default implementation uses a calculator
+ * for a sphere. Subclasses should override for creating instances of the class to be tested.
*
* @param normalized whether to force (longitude, latitude) axis order.
*/
- private static GeodeticCalculator create(final boolean normalized) {
+ GeodeticCalculator create(final boolean normalized) {
return new GeodeticCalculator(normalized ? HardCodedCRS.SPHERE : HardCodedCRS.SPHERE_φλ);
}
/**
+ * Returns a reference implementation for the given geodetic calculator.
+ */
+ private static Geodesic createReferenceImplementation(final GeodeticCalculator c) {
+ return new Geodesic(c.ellipsoid.getSemiMajorAxis(), 1/c.ellipsoid.getInverseFlattening());
+ }
+
+ /**
* Tests some simple azimuth directions. The expected directions are approximately North, East,
* South and West, but not exactly because of Earth curvature. The test verify merely that the
* azimuths are approximately correct.
+ *
+ * @throws GeodesicException if a geodesic can not be computed.
*/
@Test
- public void testCardinalAzimuths() {
+ public void testCardinalAzimuths() throws GeodesicException {
final GeodeticCalculator c = create(false);
final double tolerance = 0.2;
c.setStartGeographicPoint(20, 12);
@@ -116,9 +142,12 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
/**
* Tests azimuths at poles.
+ *
+ * @throws GeodesicException if a geodesic can not be computed.
*/
@Test
- public void testAzimuthAtPoles() {
+ @DependsOnMethod("testCardinalAzimuths")
+ public void testAzimuthAtPoles() throws GeodesicException {
final GeodeticCalculator c = create(false);
final double tolerance = 0.2;
c.setStartGeographicPoint( 90, 30);
@@ -137,15 +166,20 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
/**
* Tests geodesic distances and rhumb line length on the equator.
+ *
+ * @throws GeodesicException if a geodesic can not be computed.
*/
@Test
- public void testDistanceAtEquator() {
- final Random random = TestUtilities.createRandomNumberGenerator();
+ public void testDistanceAtEquator() throws GeodesicException {
+ final Random random = TestUtilities.createRandomNumberGenerator(86909845084528L);
final GeodeticCalculator c = create(false);
- final double r = c.ellipsoid.getSemiMajorAxis() * (PI / 180);
+ final double a = c.ellipsoid.getSemiMajorAxis();
+ final double b = c.ellipsoid.getSemiMinorAxis();
+ final double r = a * (PI / 180);
+ final double λmax = nextDown((b/a) * 180);
c.setStartGeographicPoint(0, 0);
for (int i=0; i<100; i++) {
- final double x = 360 * random.nextDouble() - 180;
+ final double x = (2*λmax) * random.nextDouble() - λmax;
c.setEndGeographicPoint(0, x);
final double expected = abs(x) * r;
assertEquals("Geodesic", expected, c.getGeodesicDistance(), Formulas.LINEAR_TOLERANCE);
@@ -162,6 +196,7 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
* @see <a href="https://en.wikipedia.org/wiki/Great-circle_navigation#Example">Great-circle navigation on Wikipedia</a>
*/
@Test
+ @DependsOnMethod({"testCardinalAzimuths", "testDistanceAtEquator"})
public void testGeodesicDistanceAndAzimuths() throws TransformException {
final GeodeticCalculator c = create(false);
c.setStartGeographicPoint(-33.0, -71.6); // Valparaíso
@@ -208,10 +243,30 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
final double λ = -71.6;
c.setStartPoint(new DirectPosition2D(λ, φ));
assertPositionEquals(λ, φ, c.getStartPoint(), Formulas.ANGULAR_TOLERANCE);
-
+ /*
+ * Test the "Valparaíso to Shanghai" geodesic given by Wikipedia, but with a larger tolerance threshold for allowing
+ * the test to pass with both spherical and ellipsoidal formula. It is not the purpose of this test to check accuracy;
+ * that check is done by `testGeodesicDistanceAndAzimuths()`.
+ */
c.setStartingAzimuth(-94.41);
c.setGeodesicDistance(18743000);
- assertPositionEquals(121.8, 31.4, c.getEndPoint(), 0.01);
+ assertPositionEquals(121.8, 31.4, c.getEndPoint(), 0.2);
+ }
+
+ /**
+ * Tests {@link GeodeticCalculator#moveToEndPoint()}.
+ *
+ * @throws TransformException if an error occurred while transforming coordinates.
+ */
+ @Test
+ public void testMoveToEndPoint() throws TransformException {
+ // Following relationship is required by GeodeticCalculator.moveToEndPoint() implementation.
+ assertEquals(GeodeticCalculator.ENDING_AZIMUTH >>> 1, GeodeticCalculator.STARTING_AZIMUTH);
+ final GeodeticCalculator c = create(false);
+ c.setStartGeographicPoint(-33.0, -71.6); // Valparaíso
+ c.setEndGeographicPoint ( 31.4, 121.8); // Shanghai
+ c.moveToEndPoint();
+ assertPositionEquals(31.4, 121.8, c.getStartPoint(), 1E-12);
}
/**
@@ -248,17 +303,18 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
@DependsOnMethod("testUsingTransform")
public void testGeodesicPath2D() throws TransformException {
final GeodeticCalculator c = create(true);
- final double tolerance = 0.05;
c.setStartGeographicPoint(-33.0, -71.6); // Valparaíso
c.setEndGeographicPoint ( 31.4, 121.8); // Shanghai
final Shape singleCurve = c.createGeodesicPath2D(Double.POSITIVE_INFINITY);
final Shape multiCurves = c.createGeodesicPath2D(10000); // 10 km tolerance.
/*
* The approximation done by a single curve is not very good, but is easier to test.
+ * The coordinate at t=0.5 has a larger tolerance threshold because it varies slightly
+ * depending on whether we are testing with spherical or ellipsoidal formulas.
*/
- assertPointEquals( -71.6, -33.0, ShapeUtilitiesExt.pointOnBezier(singleCurve, 0), tolerance);
- assertPointEquals(-238.2, 31.4, ShapeUtilitiesExt.pointOnBezier(singleCurve, 1), tolerance); // λ₂ = 121.8° - 360°
- assertPointEquals(-159.2, -6.8, ShapeUtilitiesExt.pointOnBezier(singleCurve, 0.5), tolerance);
+ assertPointEquals( -71.6, -33.0, ShapeUtilitiesExt.pointOnBezier(singleCurve, 0), 0.05);
+ assertPointEquals(-238.2, 31.4, ShapeUtilitiesExt.pointOnBezier(singleCurve, 1), 0.05); // λ₂ = 121.8° - 360°
+ assertPointEquals(-159.2, -6.9, ShapeUtilitiesExt.pointOnBezier(singleCurve, 0.5), 0.25);
/*
* The more accurate curve can not be simplified to a Java2D primitive.
*/
@@ -300,8 +356,8 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
public void testBetweenRandomPoints() throws TransformException {
final Random random = TestUtilities.createRandomNumberGenerator();
final GeodeticCalculator c = create(false);
- final Geodesic reference = new Geodesic(c.ellipsoid.getSemiMajorAxis(), 1/c.ellipsoid.getInverseFlattening());
- final StatisticsFormat sf = VERBOSE ? StatisticsFormat.getInstance() : null;
+ final Geodesic reference = createReferenceImplementation(c);
+ Statistics[] errors = null;
for (int i=0; i<100; i++) {
final double φ1 = random.nextDouble() * 180 - 90;
final double λ1 = random.nextDouble() * 360 - 180;
@@ -317,32 +373,45 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
assertEquals("Starting azimuth", expected.azi1, c.getStartingAzimuth(), Formulas.ANGULAR_TOLERANCE);
assertEquals("Ending azimuth", expected.azi2, c.getEndingAzimuth(), Formulas.ANGULAR_TOLERANCE);
assertTrue ("Rhumb ≧ geodesic", rhumbLine >= geodesic);
- if (sf != null) {
+ if (VERBOSE) {
// Checks the geodesic path on only 10% of test data, because this computation is expensive.
if ((i % 10) == 0) {
- out.println(c);
- out.println(sf.format(geodesicPathFitness(c, 1000)));
+ final Statistics[] stats = geodesicPathFitness(c, 1000);
+ if (errors == null) {
+ errors = stats;
+ } else for (int j=0; j<errors.length; j++) {
+ errors[j].combine(stats[j]);
+ }
}
}
}
+ if (errors != null) {
+ out.println("Distance between points on Bézier curve and points on geodesic.");
+ out.println("∆x/r and ∆y/r should be smaller than 1:");
+ out.println(StatisticsFormat.getInstance().format(errors));
+ }
}
/**
* Estimates the differences between the points on the Bézier curves and the points computed by geodetic calculator.
- * This method estimates the length of the path returned by {@link GeodeticCalculator#createGeodesicPath2D(double)}
- * and compares with the expected distance and azimuth at each point, by iterating over line segments and computing
- * the geodesic distance of each segment. The state of the given calculator is modified by this method.
+ * This method approximates the Bézier curve by line segments. Then for each point of the approximated Bézier curve,
+ * this method computes the location of a close point on the geodesic (more specifically a point at the same geodesic
+ * distance from the start point). The distance in metres between the two points is measured and accumulated as a
+ * fraction of the expected resolution <var>r</var>. Consequently the values in ∆x/r and ∆y/r columns should be less
+ * than 1.
+ *
+ * <div class="note"><b>Note:</b> the state of the given calculator is modified by this method.</div>
*
* @param resolution tolerance threshold for the curve approximation, in metres.
* @return statistics about errors relative to the resolution.
*/
private static Statistics[] geodesicPathFitness(final GeodeticCalculator c, final double resolution) throws TransformException {
final PathIterator iterator = c.createGeodesicPath2D(resolution).getPathIterator(null, Formulas.ANGULAR_TOLERANCE);
- final Statistics xError = new Statistics("Δx/r");
- final Statistics yError = new Statistics("Δy/r");
- final Statistics aErrors = new Statistics("Δα (°)");
+ final Statistics xError = new Statistics("∆x/r");
+ final Statistics yError = new Statistics("∆y/r");
+ final Statistics aErrors = new Statistics("∆α (°)");
final double azimuth = c.getStartingAzimuth();
- final double toMetres = (PI/180) * Formulas.getAuthalicRadius(c.ellipsoid);
+ final double toMetres = (PI/180) * c.authalicRadius;
final double[] buffer = new double[2];
while (!iterator.isDone()) {
switch (iterator.currentSegment(buffer)) {
@@ -367,6 +436,14 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
}
/**
+ * Column index for stating/ending latitude/longitude and geodesic distance.
+ * Used by {@link #compareAgainstDataset()} only.
+ */
+ static final int COLUMN_φ1 = 0, COLUMN_λ1 = 1, COLUMN_α1 = 2,
+ COLUMN_φ2 = 3, COLUMN_λ2 = 4, COLUMN_α2 = 5,
+ COLUMN_Δs = 6;
+
+ /**
* Compares computations against values provided in <cite>Karney (2010) Test set for geodesics</cite>.
* This is an optional test executed only if the {@code $SIS_DATA/Tests/GeodTest.dat} file is found.
*
@@ -375,100 +452,117 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
*/
@Test
public void compareAgainstDataset() throws IOException, TransformException {
+ int noConvergenceCount = 0;
try (LineNumberReader reader = OptionalTestData.GEODESIC.reader()) {
- final GeodeticCalculator c = new GeodeticCalculator(HardCodedCRS.WGS84_φλ);
- final Geodesic reference = new Geodesic(Formulas.getAuthalicRadius(c.ellipsoid), 0);
final Random random = TestUtilities.createRandomNumberGenerator();
- final double[] data = new double[7];
- try {
- String line;
- while ((line = reader.readLine()) != null) {
- Arrays.fill(data, Double.NaN);
- final CharSequence[] split = CharSequences.split(line, ' ');
- for (int i=min(split.length, data.length); --i >= 0;) {
- data[i] = Double.parseDouble(split[i].toString());
+ final GeodeticCalculator c = create(false);
+ final boolean isSphere = c.ellipsoid.isSphere();
+ final double[] expected = new double[7];
+ String line;
+ while ((line = reader.readLine()) != null) {
+ Arrays.fill(expected, Double.NaN);
+ final CharSequence[] split = CharSequences.split(line, ' ');
+ for (int i=min(split.length, expected.length); --i >= 0;) {
+ expected[i] = Double.parseDouble(split[i].toString());
+ }
+ /*
+ * Choose randomly whether we will test the direct or inverse geodesic problem.
+ * We execute only one test for each row instead than executing both tests,
+ * for making sure that `GeodeticCalculator` never see the expected values.
+ */
+ final boolean isTestingInverse = random.nextBoolean();
+ final double cosφ1 = abs(cos(toRadians(expected[COLUMN_φ1]))); // For adjusting longitude tolerance.
+ final double cosφ2 = abs(cos(toRadians(expected[COLUMN_φ2])));
+ double linearTolerance, latitudeTolerance, longitudeTolerance, azimuthTolerance;
+ if (isSphere) {
+ /*
+ * When spherical formulas are used instead than ellipsoidal formulas, an error up to 1% is expected
+ * in distance calculations (source: Wikipedia). A yet larger error is observed for azimuth values,
+ * especially near poles and between antipodal points. Following are empirical thresholds.
+ */
+ linearTolerance = expected[COLUMN_Δs] * 0.01;
+ latitudeTolerance = toDegrees(linearTolerance / c.authalicRadius);
+ longitudeTolerance = expected[COLUMN_φ2] > 89.5 ? 180 : latitudeTolerance / cosφ2;
+ azimuthTolerance = 0.5; // About 8.8 metres at distance of 1 km.
+ if (isTestingInverse) {
+ final double Δλ = abs(expected[COLUMN_λ2] - expected[COLUMN_λ1]);
+ if (Δλ > 179) azimuthTolerance = 100;
+ else if (Δλ > 178) azimuthTolerance = 20;
+ else if (Δλ > 175) azimuthTolerance = 10;
+ else if (Δλ > 168) azimuthTolerance = 3;
+ else if (Δλ > 158) azimuthTolerance = 1;
}
+ } else {
/*
- * We aim for an 1 cm accuracy. However when spherical formulas are used instead
- * than ellipsoidal formulas, an error up to 1% is expected (Wikipedia).
+ * When ellipsoidal formulas are used, we aim for an 1 cm accuracy in coordinate values.
+ * We also aim for azimuthd such as the error is less than 1 cm after the first 10 km.
+ * If points are nearly antipodal, we relax the azimuth tolerance threshold to 1 meter.
*/
- final double tolerance = data[6] * 0.01; // 1% of distance.
- final double cosφ = abs(cos(toRadians(data[3]))); // For adjusting longitude tolerance.
- c.setStartGeographicPoint(data[0], data[1]); // (φ₁, λ₁)
- if (random.nextBoolean()) {
- /*
- * Computes the end point from a distance and azimuth. The angular tolerance
- * is derived from the linear tolerance, except at pole where we disable the
- * check of longitude and azimuth values.
- */
- c.setStartingAzimuth (data[2]);
- c.setGeodesicDistance(data[6]);
- final double latitudeTolerance = tolerance * (1d/6371007 * (180/PI));
- final double longitudeTolerance;
- if (data[3] > 89.5) {
- longitudeTolerance = 180; // TODO: remove after we use spherical formulas.
- } else {
- longitudeTolerance = latitudeTolerance / cosφ;
+ linearTolerance = Formulas.LINEAR_TOLERANCE;
+ latitudeTolerance = Formulas.ANGULAR_TOLERANCE;
+ longitudeTolerance = Formulas.ANGULAR_TOLERANCE / cosφ2;
+ azimuthTolerance = Formulas.LINEAR_TOLERANCE * (180/PI) / 10000;
+ if (isTestingInverse) {
+ final double Δ = max(abs(180 - abs(expected[COLUMN_λ2] - expected[COLUMN_λ1])),
+ abs(expected[COLUMN_φ1] + expected[COLUMN_φ2]));
+ if (Δ < 1) {
+ azimuthTolerance = 1 * (180/PI) / 10000; // 1 meter for 10 km.
}
- final double azimuthTolerance = 0.5 / cosφ;
- compareGeodeticData(data, c, latitudeTolerance, longitudeTolerance,
- azimuthTolerance, Formulas.LINEAR_TOLERANCE);
- /*
- * Replace the distance and azimuth values by values computed using spherical formulas,
- * then compare again with values computed by GeodeticCalculator but with low tolerance.
- */
- final GeodesicData gd = reference.Direct(data[0], data[1], data[2], data[6]);
- data[3] = gd.lat2;
- data[4] = gd.lon2;
- data[5] = gd.azi2;
- } else {
- /*
- * Compute the distance and azimuth values between two points. We perform
- * this test or the above test randomly instead of always executing both
- * of them for making sure that GeodeticCalculator never see the expected
- * values.
- */
- c.setEndGeographicPoint(data[3], data[4]); // (φ₂, λ₂)
- compareGeodeticData(data, c,
- Formulas.ANGULAR_TOLERANCE, // Latitude tolerance
- Formulas.ANGULAR_TOLERANCE, // Longitude tolerance
- 100 / cosφ, tolerance); // Azimuth is inaccurate for reason not yet identified.
- /*
- * Replace the distance and azimuth values by values computed using spherical formulas,
- * then compare again with values computed by GeodeticCalculator but with low tolerance.
- */
- final GeodesicData gd = reference.Inverse(data[0], data[1], data[3], data[4]);
- data[2] = gd.azi1;
- data[5] = gd.azi2;
- data[6] = gd.s12;
}
- compareGeodeticData(data, c, Formulas.ANGULAR_TOLERANCE,
- Formulas.ANGULAR_TOLERANCE, 1E-4, Formulas.LINEAR_TOLERANCE);
}
- } catch (AssertionError e) {
- out.printf("Test failure at line %d%nGeodetic calculator is:%n", reader.getLineNumber());
- out.println(c);
- throw e;
+ /*
+ * Set input values, compute then verify results. The azimuth tolerance is divided by cos(φ).
+ * This is consistent with the ∂y/∂φ = 1/cos(φ) term in the Jacobian of Mercator projection.
+ * An error of ε in the calculation of φ causes an error of ε/cos(φ) in calculation of tan(α).
+ * Given that tan(α)+ε ≈ tan(α+ε) for small ε values, we assume that the error on α is also
+ * proportional to 1/cos(φ).
+ */
+ c.setStartGeographicPoint(expected[COLUMN_φ1], expected[COLUMN_λ1]);
+ if (isTestingInverse) {
+ c.setEndGeographicPoint(expected[COLUMN_φ2], expected[COLUMN_λ2]);
+ } else {
+ c.setStartingAzimuth (expected[COLUMN_α1]);
+ c.setGeodesicDistance(expected[COLUMN_Δs]);
+ }
+ final DirectPosition start = c.getStartPoint();
+ final DirectPosition end = c.getEndPoint();
+ try {
+ assertEquals("φ₁", expected[COLUMN_φ1], start.getOrdinate(0), Formulas.ANGULAR_TOLERANCE);
+ assertEquals("λ₁", expected[COLUMN_λ1], start.getOrdinate(1), Formulas.ANGULAR_TOLERANCE);
+ assertEquals("α₁", expected[COLUMN_α1], c.getStartingAzimuth(), azimuthTolerance / cosφ1);
+ assertEquals("φ₂", expected[COLUMN_φ2], end.getOrdinate(0), latitudeTolerance);
+ assertEquals("λ₂", expected[COLUMN_λ2], end.getOrdinate(1), longitudeTolerance);
+ assertEquals("α₂", expected[COLUMN_α2], c.getEndingAzimuth(), azimuthTolerance / cosφ2);
+ assertEquals("∆s", expected[COLUMN_Δs], c.getGeodesicDistance(), linearTolerance);
+ } catch (GeodesicException | AssertionError e) {
+ if (!isTestingInverse || e instanceof AssertionError || isFailure(expected)) {
+ out.printf("Test failure at line %d: %s%n"
+ + "The values provided in the test file are:%n"
+ + "(φ₁,λ₁) = %16.12f %16.12f%n"
+ + "(φ₂,λ₂) = %16.12f %16.12f%n"
+ + "The values computed by the geodesic calculator are:%n",
+ reader.getLineNumber(), e.getLocalizedMessage(),
+ expected[0], expected[1], expected[3], expected[4]);
+ out.println(c);
+ throw e;
+ }
+ noConvergenceCount++;
+ }
}
}
+ assertTrue("Unexpected amount of no convergence errors.", noConvergenceCount <= 8);
}
/**
- * Verifies that geodetic calculator results are equal to the given values.
- * Order in the {@code data} array is as documented in {@link OptionalTestData#GEODESIC}.
+ * Tells whether failure to compute geodesic for the given data should cause the test case to fail.
+ * The default implementation always return {@code true}. Subclass can override if some points are
+ * known to fail.
+ *
+ * @param expected a row from the {@code $SIS_DATA/Tests/GeodTest.dat} file.
+ * Use {@code COLUMN_*} constant for accessing values by column indices.
+ * @return whether the JUnit test should fail.
*/
- private static void compareGeodeticData(final double[] expected, final GeodeticCalculator c,
- final double latitudeTolerance, final double longitudeTolerance, final double azimuthTolerance,
- final double linearTolerance) throws TransformException
- {
- final DirectPosition start = c.getStartPoint();
- final DirectPosition end = c.getEndPoint();
- assertEquals("φ₁", expected[0], start.getOrdinate(0), Formulas.ANGULAR_TOLERANCE);
- assertEquals("λ₁", expected[1], start.getOrdinate(1), Formulas.ANGULAR_TOLERANCE);
- assertEquals("α₁", expected[2], c.getStartingAzimuth(), azimuthTolerance);
- assertEquals("φ₂", expected[3], end.getOrdinate(0), latitudeTolerance);
- assertEquals("λ₂", expected[4], end.getOrdinate(1), longitudeTolerance);
- assertEquals("α₂", expected[5], c.getEndingAzimuth(), azimuthTolerance);
- assertEquals("s₁₂", expected[6], c.getGeodesicDistance(), linearTolerance);
+ boolean isFailure(final double[] expected) {
+ return true;
}
}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
index 0c5c4e8..f46a2db 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
@@ -246,6 +246,7 @@ import org.junit.BeforeClass;
org.apache.sis.referencing.cs.CodesTest.class,
org.apache.sis.referencing.CRSTest.class,
org.apache.sis.referencing.GeodeticCalculatorTest.class,
+ org.apache.sis.referencing.GeodesicsOnEllipsoidTest.class,
org.apache.sis.internal.referencing.DefinitionVerifierTest.class,
org.apache.sis.internal.referencing.CoordinateOperationsTest.class,
diff --git a/core/sis-utility/src/main/java/org/apache/sis/math/MathFunctions.java b/core/sis-utility/src/main/java/org/apache/sis/math/MathFunctions.java
index f0ee80c..48be82e 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/math/MathFunctions.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/math/MathFunctions.java
@@ -23,9 +23,13 @@ import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.resources.Errors;
import org.apache.sis.internal.util.DoubleDouble;
+import static java.lang.Math.abs;
+import static java.lang.Math.sqrt;
import static java.lang.Float.intBitsToFloat;
+import static java.lang.Float.floatToIntBits;
import static java.lang.Float.floatToRawIntBits;
import static java.lang.Double.longBitsToDouble;
+import static java.lang.Double.doubleToLongBits;
import static java.lang.Double.doubleToRawLongBits;
import static org.apache.sis.internal.util.Numerics.SIGN_BIT_MASK;
import static org.apache.sis.internal.util.Numerics.SIGNIFICAND_SIZE;
@@ -223,7 +227,7 @@ public final class MathFunctions extends Static {
// We have found a non-zero element. If it is the only one, returns its absolute value.
double v2;
- do if (i == 0) return Math.abs(v1);
+ do if (i == 0) return abs(v1);
while ((v2 = vector[--i]) == 0);
// If there is exactly 2 elements, use Math.hypot which is more robust than our algorithm.
@@ -406,7 +410,7 @@ public final class MathFunctions extends Static {
* @since 0.6
*/
public static double asinh(final double x) {
- return Math.log(x + Math.sqrt(x*x + 1));
+ return Math.log(x + sqrt(x*x + 1));
}
/**
@@ -421,7 +425,7 @@ public final class MathFunctions extends Static {
* @since 0.6
*/
public static double acosh(final double x) {
- return Math.log(x + Math.sqrt(x*x - 1));
+ return Math.log(x + sqrt(x*x - 1));
}
/**
@@ -601,7 +605,7 @@ public final class MathFunctions extends Static {
* @return {@code true} if both values are equal given the tolerance threshold.
*/
public static boolean epsilonEqual(final float v1, final float v2, final float ε) {
- return (Math.abs(v1 - v2) <= ε) || Float.floatToIntBits(v1) == Float.floatToIntBits(v2);
+ return (abs(v1 - v2) <= ε) || floatToIntBits(v1) == floatToIntBits(v2);
}
/**
@@ -623,7 +627,7 @@ public final class MathFunctions extends Static {
* @return {@code true} if both values are equal given the tolerance threshold.
*/
public static boolean epsilonEqual(final double v1, final double v2, final double ε) {
- return (Math.abs(v1 - v2) <= ε) || Double.doubleToLongBits(v1) == Double.doubleToLongBits(v2);
+ return (abs(v1 - v2) <= ε) || doubleToLongBits(v1) == doubleToLongBits(v2);
}
/**
@@ -777,7 +781,7 @@ public final class MathFunctions extends Static {
primes = Arrays.copyOf(primes, Math.min((index | 0xF) + 1, PRIMES_LENGTH_16_BITS));
do {
testNextNumber: while (true) { // Simulate a "goto" statement (usually not recommanded...)
- final int stopAt = (int) Math.sqrt(n += 2);
+ final int stopAt = (int) sqrt(n += 2);
int prime;
int j = 0;
do {
@@ -842,7 +846,7 @@ testNextNumber: while (true) { // Simulate a "goto" statement (usua
if (number == 0) {
return ArraysExt.EMPTY_INT;
}
- number = Math.abs(number);
+ number = abs(number);
int[] divisors = new int[16];
divisors[0] = 1;
int count = 1;
@@ -852,7 +856,7 @@ testNextNumber: while (true) { // Simulate a "goto" statement (usua
* values before that point, i.e. if n=p1*p2 and p2 is greater than 'sqrt', than p1
* most be lower than 'sqrt'.
*/
- final int sqrt = (int) Math.sqrt(number); // Really want rounding toward 0.
+ final int sqrt = (int) sqrt(number); // Really want rounding toward 0.
for (int p,i=0; (p=primeNumberAt(i)) <= sqrt; i++) {
if (number % p == 0) {
if (count == divisors.length) {
@@ -925,7 +929,7 @@ testNextNumber: while (true) { // Simulate a "goto" statement (usua
*/
int minValue = Integer.MAX_VALUE;
for (int i=0; i<numbers.length; i++) {
- final int n = Math.abs(numbers[i]);
+ final int n = abs(numbers[i]);
if (n <= minValue) {
minValue = n;
}
@@ -937,7 +941,7 @@ testNextNumber: while (true) { // Simulate a "goto" statement (usua
*/
int count = divisors.length;
for (int i=0; i<numbers.length; i++) {
- final int n = Math.abs(numbers[i]);
+ final int n = abs(numbers[i]);
if (n != minValue) {
for (int j=count; --j>0;) { // Do not test j==0, since divisors[0] == 1.
if (n % divisors[j] != 0) {
@@ -948,4 +952,229 @@ testNextNumber: while (true) { // Simulate a "goto" statement (usua
}
return ArraysExt.resize(divisors, count);
}
+
+ /**
+ * Returns the real (non-complex) roots of a polynomial equation having the given coefficients.
+ * This method returns the <var>x</var> values for which <var>y</var>=0 in the following equation:
+ *
+ * <blockquote><var>y</var> =
+ * <var>c<sub>0</sub></var> +
+ * <var>c<sub>1</sub></var>⋅<var>x</var> +
+ * <var>c<sub>2</sub></var>⋅<var>x</var><sup>2</sup> +
+ * <var>c<sub>3</sub></var>⋅<var>x</var><sup>3</sup> + … +
+ * <var>c<sub>n</sub></var>⋅<var>x</var><sup>n</sup>
+ * </blockquote>
+ *
+ * Current implementation can resolve polynomials described by a maximum of 5 coefficients, ignoring
+ * leading and trailing zeros. They correspond to linear, quadratic, cubic and quartic polynomials.
+ *
+ * @param coefficients the <var>c<sub>0</sub></var>, <var>c<sub>1</sub></var>, <var>c<sub>2</sub></var>, …
+ * <var>c<sub>n</sub></var> coefficients, in that order.
+ * @return the non-complex roots, or an empty array if none.
+ * @throws UnsupportedOperationException if given arguments contain more non-zero coefficients than this method can handle.
+ *
+ * @see java.awt.geom.QuadCurve2D#solveQuadratic(double[])
+ * @see java.awt.geom.CubicCurve2D#solveCubic(double[])
+ *
+ * @since 1.0
+ */
+ public static double[] polynomialRoots(final double... coefficients) {
+ int upper = coefficients.length;
+ while (upper > 0) {
+ double a = coefficients[--upper]; // Coefficient of the term with highest exponent.
+ if (a == 0) {
+ continue; // Search the last non-zero coefficient.
+ }
+ double c; // Coefficient of the term with lowest exponent.
+ int lower = 0;
+ while ((c = coefficients[lower]) == 0) lower++;
+ switch (upper - lower) {
+ /*
+ * c = 0
+ *
+ * Can not compute x. We could return an arbitrary value if c = 0, but we rather return
+ * an empty array for keeping the number of root equals to the highest exponent.
+ */
+ case 0: {
+ break;
+ }
+ /*
+ * ax + c = 0 → x = -c/a
+ */
+ case 1: {
+ final double x = -c / a;
+ if (Double.isNaN(x)) break;
+ return new double[] {x};
+ }
+ /*
+ * ax² + bx + c = 0 → x = (-b ± √[b² - 4ac]) / (2a)
+ *
+ * Above equation is numerically unstable. More stable algorithm is given
+ * by Numerical Recipes §5.6 (quadratic equation), which is applied below.
+ */
+ case 2: {
+ final double b = coefficients[lower + 1];
+ final double q = -0.5 * (b + Math.copySign(sqrt(b*b - 4*a*c), b));
+ final double x1 = q/a;
+ final double x2 = c/q;
+ if (Double.isNaN(x1) && Double.isNaN(x2)) break;
+ return (x1 != x2) ? new double[] {x1, x2} : new double[] {x1};
+ }
+ /*
+ * x³ + ax² + bx + c = 0
+ *
+ * Numerical Recipes §5.6 (cubic equation) is applied below.
+ * Solution usually have either 1 or 3 roots.
+ */
+ case 3: {
+ return refineRoots(coefficients, solveCubic(
+ coefficients[lower + 2] / a,
+ coefficients[lower + 1] / a,
+ c / a, false));
+ }
+ /*
+ * x⁴ + ax³ + bx² + cx + d = 0
+ *
+ * https://dlmf.nist.gov/1.11 in "Quartic equations" section.
+ * This algorithm reduces the equation to a cubic equation.
+ */
+ case 4: {
+ double b,d;
+ d = c / a;
+ c = coefficients[lower + 1] / a;
+ b = coefficients[lower + 2] / a;
+ a = coefficients[lower + 3] / a;
+ final double a2 = a*a;
+ final double p = -3./8 * (a2) + b;
+ final double q = 1./8 * (a2*a) - 1./2 * (a*b) + c;
+ final double r = -3./256 * (a2*a2) + 1./16 * (a2*b) - 1./4 * (a*c) + d;
+ final double[] roots = solveCubic(-2*p, p*p-4*r, q*q, true);
+ if (roots.length != 4) break;
+ for (int i=0; i<3; i++) {
+ roots[i] = sqrt(-roots[i]);
+ }
+ if (isPositive(q)) {
+ for (int i=0; i<3; i++) {
+ roots[i] = -roots[i];
+ }
+ }
+ final double α = roots[0], β = roots[1], γ = roots[2];
+ final double s = α + β + γ;
+ if (Double.isNaN(s)) break;
+ roots[0] = s/2 - (a /= 4);
+ roots[1] = (+α - β - γ)/2 - a;
+ roots[2] = (-α + β - γ)/2 - a;
+ roots[3] = (-α - β + γ)/2 - a;
+ return refineRoots(coefficients, removeDuplicated(roots));
+ }
+ default: {
+ throw new UnsupportedOperationException();
+ }
+ }
+ break;
+ }
+ return ArraysExt.EMPTY_DOUBLE;
+ }
+
+ /**
+ * Solves cubic equation x³ + ax² + bx + c = 0. The solution before simplification has either 1 or 3 roots.
+ * The {@code quartic} argument specifies whether this cubic equation is used as a step for solving a quartic equation:
+ *
+ * <ul>
+ * <li>If {@code true}, then we are interested only in the 3 roots solution and we do not check for duplicated values.
+ * The length of returned array is 4 for allowing the caller to reuse the same array.</li>
+ * <li>If {@code false}, then this method may simplify the 3 roots to 2 roots if two of them are equal,
+ * or may return the 1 root solution. The length of returned array is the number of roots.</li>
+ * </ul>
+ */
+ private static double[] solveCubic(double a, double b, double c, final boolean quartic) {
+ final double Q = (a*a - 3*b) / 9; // Q from Numerical Recipes 5.6.10.
+ final double R = (a*(a*a - 4.5*b) + 13.5*c) / 27; // R from Numerical Recipes 5.6.10.
+ final double Q3 = Q*Q*Q;
+ final double R2 = R*R;
+ a /= 3; // Last term of Numerical Recipes 5.6.12, 17 and 18.
+ if (R2 < Q3) {
+ /*
+ * Numerical Recipes 5.6.11 and 5.6.12 uses acos(R/sqrt(Q³)). It is possible to rewrite as
+ * atan2(sqrt(Q3 - R2), R) using the cos(θ) = 1/√[1 + tan²θ] trigonometric identity, but
+ * this substitution seems to decrease accuracy instead of increasing it in our tests.
+ */
+ b = Math.acos(R/sqrt(Q3)) / 3; // θ from Numerical recipes 5.6.11, then b = θ/3.
+ c = -2 * sqrt(Q); // First part of Numerical Recipes 5.6.12.
+ double[] roots = new double[quartic ? 4 : 3];
+ roots[2] = c*Math.cos(b - 2*Math.PI/3) - a; // TODO: try Math.fma with JDK9.
+ roots[1] = c*Math.cos(b + 2*Math.PI/3) - a;
+ roots[0] = c*Math.cos(b) - a;
+ if (!quartic) {
+ roots = removeDuplicated(roots);
+ }
+ return roots;
+ }
+ if (!quartic) {
+ b = -Math.copySign(Math.cbrt(abs(R) + sqrt(R2 - Q3)), R); // A from Numerical Recipes 5.6.15.
+ final double x = (b == 0 ? 0 : b + Q/b) - a;
+ if (!Double.isNaN(x)) {
+ return new double[] {x};
+ }
+ }
+ return ArraysExt.EMPTY_DOUBLE;
+ }
+
+ /**
+ * Remove duplicated values in the given array. This method is okay only for very small arrays (3 or 4 elements).
+ * Duplicated values should be very rare and occur mostly as a consequence of rounding errors while computing the
+ * roots of polynomial equations. Because if the algebraic solution has less roots than what we would expect from
+ * the largest exponent (for example ax² + bx = 0 has only one root instead of two), then {@link #polynomialRoots}
+ * should have reduced the equation to a lower degrees (ax + b = 0 in above example), in which case there is no
+ * duplicated roots to remove.
+ */
+ private static double[] removeDuplicated(double[] roots) {
+ int i = 1;
+next: while (i < roots.length) {
+ for (int j=i; --j >= 0;) {
+ if (roots[j] == roots[i]) {
+ roots = ArraysExt.remove(roots, i, 1);
+ continue next;
+ }
+ }
+ i++;
+ }
+ return roots;
+ }
+
+ /**
+ * Tries to improves accuracy of polynomial roots by applying small displacements
+ * to the <var>x</var> values using ∂y/∂x derivative around those values.
+ *
+ * <div class="note"><b>Purpose:</b>
+ * this refinement is significant in a {@link org.apache.sis.referencing.GeodesicsOnEllipsoid}
+ * test checking the value of an μ(x²,y²) function.</div>
+ *
+ * @param coefficients the user-specified coefficients.
+ * @param roots the roots. This array will be modified in place.
+ * @return {@code roots}.
+ */
+ private static double[] refineRoots(final double[] coefficients, final double[] roots) {
+ for (int i=0; i < roots.length; i++) {
+ double ymin = Double.POSITIVE_INFINITY;
+ double x = roots[i];
+ double dx;
+ do {
+ double px = 1; // Power of x: 1, x¹, x², x³, …
+ double dy = 0; // First derivative of polynomial at x.
+ double y = coefficients[0]; // Value of polynomial at x.
+ double ey = 0, edy = 0; // Error terms for Kahan summation algorithm.
+ for (int j=1; j<coefficients.length; j++) {
+ final double c = coefficients[j];
+ double s;
+ s = c * (px * i) + edy; edy = s + (dy - (dy += s)); // Kahan summation of dy.
+ s = c * (px *= x) + ey; ey = s + ( y - ( y += s)); // Kahan summation of y.
+ }
+ if (!(ymin > (ymin = abs(y)))) break; // If result not better than previous result, stop.
+ roots[i] = x;
+ dx = y/dy;
+ } while (x != (x -= dx) && Double.isFinite(x));
+ }
+ return roots;
+ }
}
diff --git a/core/sis-utility/src/test/java/org/apache/sis/math/MathFunctionsTest.java b/core/sis-utility/src/test/java/org/apache/sis/math/MathFunctionsTest.java
index 210e6f0..b899114 100644
--- a/core/sis-utility/src/test/java/org/apache/sis/math/MathFunctionsTest.java
+++ b/core/sis-utility/src/test/java/org/apache/sis/math/MathFunctionsTest.java
@@ -438,4 +438,47 @@ public final strictfp class MathFunctionsTest extends TestCase {
assertArrayEquals(new int[] {1, 2}, commonDivisors(-4, 2));
assertArrayEquals(new int[] {}, commonDivisors(0));
}
+
+ /**
+ * Tests {@link MathFunctions#polynomialRoots(double...)} with values given in NIST example.
+ * Source: <a href="https://dlmf.nist.gov/1.11">NIST Digital Library of Mathematical Functions</a>.
+ */
+ @Test
+ public void testQuarticRoots() {
+ /*
+ * Given polynomial: f(z) = z⁴ - 4z³ + 5z + 2
+ * Reduced: g(w) = w⁴ - 6w² - 3w + 4
+ * Resolvent cubic: z³ + 12z² + 20z + 9 = 0
+ */
+ assertArrayEquals(new double[] {
+ 3.5615528128088303, // ½(3 + √17)
+ 1.618033988749895, // ½(1 + √5)
+ -0.6180339887498949, // ½(1 - √5)
+ -0.5615528128088303 // ½(3 - √17)
+ }, polynomialRoots(2, 5, 0, -4, 1), 1E-14);
+ }
+
+ /**
+ * Tests {@link MathFunctions#polynomialRoots(double...)} with random values.
+ */
+ @Test
+ public void testPolynomialRoots() {
+ final Random r = TestUtilities.createRandomNumberGenerator(62308330810132L);
+ final double[] coefficients = new double[5];
+ for (int n=2; n < coefficients.length; n++) { // Number of coefficients to fill.
+ for (int t=n*10; --t >= 0;) { // Number of tests to execute.
+ for (int i=0; i<n; i++) {
+ coefficients[i] = (r.nextInt(8) != 0) ? r.nextDouble() * 20 - 10 : 0;
+ }
+ for (final double x : polynomialRoots(coefficients)) {
+ double p = 1;
+ double y = coefficients[0];
+ for (int i=1; i<coefficients.length; i++) {
+ y += coefficients[i] * (p *= x);
+ }
+ assertEquals(0, y, 1E-12);
+ }
+ }
+ }
+ }
}
diff --git a/core/sis-utility/src/test/java/org/apache/sis/test/TestCase.java b/core/sis-utility/src/test/java/org/apache/sis/test/TestCase.java
index da96ed4..b2ae8f1 100644
--- a/core/sis-utility/src/test/java/org/apache/sis/test/TestCase.java
+++ b/core/sis-utility/src/test/java/org/apache/sis/test/TestCase.java
@@ -126,6 +126,9 @@ public abstract strictfp class TestCase {
out = new PrintWriter(buffer = new StringWriter());
VERBOSE = Boolean.getBoolean(TestConfiguration.VERBOSE_OUTPUT_KEY);
RUN_EXTENSIVE_TESTS = Boolean.getBoolean(TestConfiguration.EXTENSIVE_TESTS_KEY);
+ if (VERBOSE) {
+ System.setErr(System.out); // For avoiding log records to be interleaved with block of text.
+ }
}
/**
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/index/tree/QuadTree.java b/storage/sis-storage/src/main/java/org/apache/sis/index/tree/QuadTree.java
index 92fcdea..9bbce43 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/index/tree/QuadTree.java
+++ b/storage/sis-storage/src/main/java/org/apache/sis/index/tree/QuadTree.java
@@ -25,6 +25,7 @@ import org.apache.sis.geometry.Envelope2D;
import org.apache.sis.distance.LatLonPointRadius;
import org.apache.sis.referencing.CommonCRS;
import org.apache.sis.referencing.GeodeticCalculator;
+import org.apache.sis.referencing.operation.GeodesicException;
/**
* Implementation of Quad Tree Index. Insertion algorithm implemented based on
@@ -342,9 +343,9 @@ public class QuadTree {
if (node == null) {
return matches;
} else if (node.getNodeType() != NodeType.GRAY) {
- if (node.getNodeType() == NodeType.WHITE)
+ if (node.getNodeType() == NodeType.WHITE) {
return matches;
- else {
+ } else {
QuadTreeData[] data = node.getData();
for (int i = 0; i < node.getCount(); i++) {
DirectPosition2D latLon = data[i].getLatLon();
@@ -352,7 +353,11 @@ public class QuadTree {
synchronized (calculator) {
calculator.setStartGeographicPoint(latLon.y, latLon.x);
calculator.setEndGeographicPoint(point.y, point.x);
- distance = calculator.getGeodesicDistance();
+ try {
+ distance = calculator.getGeodesicDistance();
+ } catch (GeodesicException e) {
+ throw new RuntimeException(e);
+ }
}
if (distance <= radiusKM) {
matches.add(data[i]);
@@ -360,7 +365,6 @@ public class QuadTree {
}
return matches;
}
-
} else {
Rectangle2D swRectangle = new Rectangle2D.Double(nodeRegion.getX(), nodeRegion.getY(),
nodeRegion.getWidth() / 2, nodeRegion.getHeight() / 2);
|