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
commit 4b6bb5fd7844d977e14957a5d618b8f91873d4f7 Author: Martin Desruisseaux <[email protected]> AuthorDate: Fri Aug 2 09:48:59 2019 +0200 Omit redundant sinα0 field. --- .../apache/sis/internal/referencing/Formulas.java | 5 +- .../sis/referencing/GeodesicsOnEllipsoid.java | 86 +++++++++++++--------- 2 files changed, 54 insertions(+), 37 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 b3b6afe..6cce1df 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 @@ -89,10 +89,9 @@ public final class Formulas extends Static { * 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> + * <p>Current value has been determined empirically for allowing {@code GeodesicsOnEllipsoidTest} to pass.</p> */ - public static final int MAXIMUM_ITERATIONS = 16; + public static final int MAXIMUM_ITERATIONS = 18; /** * 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 index ae60e24..5da7fd3 100644 --- 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 @@ -170,13 +170,27 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { * 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. + * + * <p>{@code sinα0} field is omitted because it is equal to {@link #msinα2} in this class. + * The {@link #sinα0()} method is used as a synonymous where sin(α₀) value was expected.</p> + */ + private double cosα0; + + /** + * The sin(α₀) value, which is identical to m⋅sin(α₂) in this class. We use the inherited field because we noticed + * that all calculations of sin(α₀) were followed by a {@code msinα2 = sinα0} statement, and using {@code msinα2} + * directly makes some code clearer. We define this method for the cases where we want to make clear that sin(α₀) + * value was intended. */ - private double sinα0, cosα0; + private double sinα0() { + return msinα2; + } /** * 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. + * The {@link #COEFFICIENTS_FOR_START_POINT} flag specifies whether this field needs to be recomputed. * * @see #sphericalToGeodeticLongitude(double, double) */ @@ -187,6 +201,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { * 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. + * The {@link #COEFFICIENTS_FOR_START_POINT} flag specifies whether this field needs to be recomputed. * * @see #sphericalToEllipsoidalAngle(double, boolean) */ @@ -195,6 +210,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { /** * 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). + * The {@link #COEFFICIENTS_FOR_START_POINT} flag specifies whether this field needs to be recomputed. */ private double ε; @@ -251,7 +267,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { * @return k² = ℯ′²⋅cos²α₀ term (Karney 9). */ private double computeSeriesExpansionCoefficients() { - assert abs(sinα0*sinα0 + cosα0*cosα0 - 1) <= 1E-10; // Arbitrary threshold. + 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) @@ -448,7 +464,18 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { + C32 * sin(2 * 2 * σ) + C31 * sin(2 * 1 * σ) + σ); if (STORE_LOCAL_VARIABLES) store("I₃(σ)", I3); - return ω - sinα0 * I3 * (1 - axisRatio); + return ω - sinα0() * I3 * (1 - axisRatio); + } + + /** + * Sets the {@link #sinα0} and {@link #cosα0} terms. Note that the {@code sinα0} value is actually + * stored in the {@link #msinα2} field, because their values are always equal in this class. + * + * @see #sinα0() + */ + private void α0(final double sinα1, final double cosα1, final double sinβ1, final double cosβ1) { + msinα2 = 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. } /** @@ -472,8 +499,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { 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. + α0(sinα1, cosα1, sinβ1, cosβ1); /* * 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 @@ -481,8 +507,8 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { * 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 σ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); @@ -507,7 +533,6 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { * 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. */ - msinα2 = sinα0; mcosα2 = cosα0 * cosσ2; /* * Ending point coordinates on auxiliary sphere: Latitude β is given by Karney equation 13: @@ -520,8 +545,8 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { * 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(sinα0, mcosα2); - final double ω2 = atan2(sinα0*sinσ2, cosσ2); // (Karney 12). + final double cosβ2 = hypot(msinα2, mcosα2); // m⋅sin(α₂) = sin(α₀) in this class. + 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 φ. @@ -534,7 +559,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { store("s₂", s2b * semiMinor); store("τ₂", s2b / A1); store("σ₂", σ2); - store("α₂", atan2(sinα0, cosα0*cosσ2)); + store("α₂", atan2(msinα2, mcosα2)); store("β₂", atan2(sinβ2, cosβ2)); store("ω₂", ω2); store("λ₂", λ2E); @@ -704,11 +729,9 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { * 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. + int moreRefinements = Formulas.MAXIMUM_ITERATIONS; + do { + α0(msinα1 = sin(α1), mcosα1 = cos(α1), sinβ1, cosβ1); final double k2 = computeSeriesExpansionCoefficients(); /* * Compute α₂ from Karney equation 5: sin(α₀) = sin(α₁)⋅cos(β₁) = sin(α₂)⋅cos(β₂) @@ -722,10 +745,11 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { * 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α1_cosβ1 = mcosα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))); + mcosα2 = cosα2_cosβ2; /* * Karney gives the following formulas: * @@ -739,10 +763,10 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { * 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); + σ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. @@ -754,8 +778,9 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { λ1E = sphericalToGeodeticLongitude(ω1, σ1); λ2E = sphericalToGeodeticLongitude(ω2, σ2); final double Δλ_error = IEEEremainder(λ2E - λ1E - Δλ, 2*PI); - final boolean done = (abs(Δλ_error) <= Formulas.ANGULAR_TOLERANCE * (PI/180) / ACCURACY_IMPROVEMENT); - if (--nbIter < 0 && !done) { + if ((abs(Δλ_error) <= Formulas.ANGULAR_TOLERANCE * (PI/180) / ACCURACY_IMPROVEMENT)) { + moreRefinements = 0; + } else if (--moreRefinements == 0) { throw new GeodesicException(Resources.format(Resources.Keys.NoConvergence)); } /* @@ -765,7 +790,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { * 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) { + if (abs(mcosα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 { /* @@ -802,7 +827,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { store("σ₂", σ2); store("ω₁", ω1); store("ω₂", ω2); - store("α₂", atan2(sinα0, cosα2_cosβ2)); + store("α₂", atan2(msinα2, mcosα2)); store("λ₂", λ2E); store("Δλ", λ2E - λ1E); store("δλ", Δλ_error); @@ -812,14 +837,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { store("s₂", I1_σ2 * semiMinor); store("Δs", (I1_σ2 - I1_σ1) * semiMinor); } - if (done) { - msinα1 = sinα1; - mcosα1 = cosα1; - msinα2 = sinα0; - mcosα2 = cosα2_cosβ2; - break; - } - } + } while (moreRefinements != 0); final double I1_σ2; I1_σ1 = sphericalToEllipsoidalAngle(σ1, false); I1_σ2 = sphericalToEllipsoidalAngle(σ2, false); @@ -908,7 +926,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator { store("A₁", A1); store("A₂", A2); store("A₃", A3); - store("α₀", atan2(sinα0, cosα0)); + store("α₀", atan2(sinα0(), cosα0)); store("I₁(σ₁)", I1_σ1); store("s₁", I1_σ1 * semiMinor); store("λ₁", λ1E);
