sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 02/02: Rewrite Bennett equation 1 by merging ΔΨ = Ψ(φ₂) - Ψ(φ₁) in a single step. It reduces the number of calls to Math.log and Math.pow.
Date Wed, 07 Aug 2019 15:52:23 GMT
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 ad6c9d8b656ecacb918b65d5387bea44183255b8
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed Aug 7 07:28:17 2019 +0200

    Rewrite Bennett equation 1 by merging ΔΨ = Ψ(φ₂) - Ψ(φ₁) in a single step.
    It reduces the number of calls to Math.log and Math.pow.
---
 .../sis/referencing/GeodesicsOnEllipsoid.java      | 42 +++++++++++++---------
 .../sis/referencing/GeodesicsOnEllipsoidTest.java  |  4 ---
 2 files changed, 25 insertions(+), 21 deletions(-)

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 9a2fc18..9d4c05c 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
@@ -944,10 +944,32 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
     @Override
     final void computeRhumbLine() {
         canComputeDistance();
+        /*
+         * Bennett (1996) equation 1 computes isometric latitudes Ψ for given geodetic latitudes
φ:
+         *
+         *     Ψ(φ) = log(tan(PI/4 + φ/2) * pow((1 - ℯsinφ) / (1 + ℯsinφ), ℯ/2));
+         *
+         * (ℯ is the eccentricity, not squared). However we need only the isometric latitudes
difference:
+         *
+         *     ΔΨ = Ψ(φ₂) - Ψ(φ₁)
+         *
+         * We rewrite the equation using log(Ψ₁) - log(Ψ₂) = log(Ψ₁/Ψ₂) and other
identities:
+         *
+         *     ΔΨ = log(tan(PI/4 + φ₂/2)) + log(pow((1 - ℯsinφ₂) / (1 + ℯsinφ₂),
ℯ/2))
+         *        - log(tan(PI/4 + φ₁/2)) - log(pow((1 - ℯsinφ₁) / (1 + ℯsinφ₁),
ℯ/2))
+         *
+         *        = log(tan(PI/4 + φ₂/2) /
+         *              tan(PI/4 + φ₁/2) *
+         *              pow(((1 - ℯsinφ₂)*(1 + ℯsinφ₁)) /
+         *                  ((1 + ℯsinφ₂)*(1 - ℯsinφ₁)), ℯ/2))
+         */
+        final double eccentricity = sqrt(eccentricitySquared);      // TODO: avoid computing
on each invocation.
+        final double sinφ1 = sin(φ1);
+        final double sinφ2 = sin(φ2);
+        final double sd = eccentricity * (sinφ1 - sinφ2);
+        final double sm = 1 - eccentricitySquared * (sinφ1 * sinφ2);
+        final double ΔΨ = log(tan(PI/4 + φ2/2) / tan(PI/4 + φ1/2) * pow((sm+sd)/(sm-sd),
eccentricity/2));
         final double Δλ = λ2 - λ1;
-        final double Ψ1 = Ψ(φ1);                    // Bennett equation 1.
-        final double Ψ2 = Ψ(φ2);
-        final double ΔΨ = Ψ2 - Ψ1;
         final double azimuth = atan2(Δλ, ΔΨ);
         final double S;
         if (abs(φ1 - φ2) < LATITUDE_THRESHOLD) {
@@ -967,26 +989,12 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
         rhumblineLength = S * (semiMinor / axisRatio);      // TODO: compute semiMajor only
once.
         if (STORE_LOCAL_VARIABLES) {
             store("Δλ", Δλ);
-            store("Ψ₁", Ψ1);
-            store("Ψ₂", Ψ2);
             store("ΔΨ", ΔΨ);
             store("C",  azimuth);
         }
     }
 
     /**
-     * Returns the isometric latitude Ψ for the given geodetic latitude φ.
-     * This is equation 1 in Bennett (1996).
-     *
-     * @see org.apache.sis.referencing.operation.projection.ConformalProjection#expΨ
-     */
-    private double Ψ(final double φ) {
-        final double eccentricity = sqrt(eccentricitySquared);      // TODO: avoid computing
on each invocation.
-        final double ℯsinφ = eccentricity * sin(φ);
-        return log(tan(PI/4 + φ/2) * pow((1 - ℯsinφ) / (1 + ℯsinφ), eccentricity/2));
-    }
-
-    /**
      * Computes Bennett (1996) equation 2.
      */
     private double m(final double φ) {
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
index 80c27a1..7344e19 100644
--- 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
@@ -524,8 +524,6 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT
         final double distance = testedEarth.getRhumblineLength();
         final double scale = testedEarth.ellipsoid.getSemiMajorAxis() / NAUTICAL_MILE;
         assertValueEquals("Δλ", 0, 75+35.4 / 60,         1E-11, true);
-        assertValueEquals("Ψ₁", 0,  617.64 / (10800/PI), 1E-5, false);
-        assertValueEquals("Ψ₂", 0, 3794.54 / (10800/PI), 1E-5, false);
         assertValueEquals("ΔΨ", 0, 3176.89 / (10800/PI), 1E-5, false);
         assertValueEquals("m₁", 0,  615.43 / scale,      1E-6, false);
         assertValueEquals("m₂", 0, 3201.59 / scale,      1E-6, false);
@@ -545,8 +543,6 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT
         testedEarth.setEndGeographicPoint  (-53-10.8/60, -41-34.6/60);
         final double distance = testedEarth.getRhumblineLength();
         assertValueEquals("Δλ", 0,  55+57.0 / 60,         1E-11, true);
-        assertValueEquals("Ψ₁", 0, -3725.18 / (10800/PI), 1E-5, false);
-        assertValueEquals("Ψ₂", 0, -3763.30 / (10800/PI), 1E-5, false);
         assertValueEquals("ΔΨ", 0,   -38.12 / (10800/PI), 1E-5, false);
         assertValueEquals("C",  0,  90.6505,              1E-4, true);
         assertEquals("distance", 2028.9 * NAUTICAL_MILE, distance, 0.05 * NAUTICAL_MILE);


Mime
View raw message