sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 02/02: Initial draft of Rhumb Line distance calculated using formulas published by: G.G. Bennett, 1996. Practical Rhumb Line Calculations on the Spheroid. J. Navigation 49(1), 112--119. https://doi.org/10.1017/S0373463300013151
Date Mon, 05 Aug 2019 15:32:42 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 e3d94efa36c9f2a68158753ea87828dea3fa383f
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Mon Aug 5 17:30:55 2019 +0200

    Initial draft of Rhumb Line distance calculated using formulas published by:
    G.G. Bennett, 1996. Practical Rhumb Line Calculations on the Spheroid. J. Navigation 49(1),
112--119.
    https://doi.org/10.1017/S0373463300013151
---
 .../sis/referencing/GeodesicsOnEllipsoid.java      | 67 ++++++++++++++++++++++
 .../sis/referencing/GeodesicsOnEllipsoidTest.java  | 38 +++++++++---
 2 files changed, 97 insertions(+), 8 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 013cfb6..838f05d 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
@@ -210,6 +210,11 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
     private double A1, A2, A3, C31, C32, C33, C34, C35;
 
     /**
+     * Coefficients for Rhumb line calculation from Bennett (1996) equation 2.
+     */
+    private final double R0, R2, R4, R6;
+
+    /**
      * Constructs a new geodetic calculator expecting coordinates in the supplied CRS.
      *
      * @param  crs         the referencing system for the {@link Position} arguments and
return values.
@@ -225,6 +230,14 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
         thirdFlattening = (a - b) / (a + b);
         axisRatio = b / a;
         semiMinor = b;
+        /*
+         * For rhumb-line calculation.
+         */
+        double fe;
+        R0 = 1  +  (eccentricitySquared)*(-1./4   + eccentricitySquared*(-3./64 + eccentricitySquared*(-5./256)));
+        R2 = (fe  = eccentricitySquared)*( 3./8   + eccentricitySquared*( 3./32 + eccentricitySquared*(45./1024)));
+        R4 = (fe *= eccentricitySquared)*(15./256 + eccentricitySquared*(45./1024));
+        R6 = (fe *  eccentricitySquared)*(35./3072);
     }
 
     /**
@@ -921,4 +934,58 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
     double computedToGiven(double α1) {
         return α1;
     }
+
+    /**
+     * Computes rhumb line using series expansion.
+     *
+     * <p><b>Source:</b> G.G. Bennett, 1996. <a href="https://doi.org/10.1017/S0373463300013151">Practical
+     * Rhumb Line Calculations on the Spheroid</a>. J. Navigation 49(1), 112--119.</p>
+     */
+    @Override
+    void computeRhumbLine() {
+        canComputeDistance();
+        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 m1 = m(φ1);
+        final double m2 = m(φ2);
+        final double φm = (φ1 + φ2)/2;
+        final double sinφ = sin(φm);
+        double S = (m2 - m1) / cos(azimuth);
+        if (φ1 == φ2) {
+            S = Δλ * cos(φm) / (sin(azimuth) * (1 - eccentricitySquared*(sinφ*sinφ)));
    // Bennett equation 4.
+        }
+        rhumblineLength = S * (semiMinor / axisRatio);      // TODO: compute semiMajor only
once.
+        if (STORE_LOCAL_VARIABLES) {
+            store("Δλ", Δλ);
+            store("Ψ₁", Ψ1);
+            store("Ψ₂", Ψ2);
+            store("ΔΨ", ΔΨ);
+            store("m₁", m1);
+            store("m₂", m2);
+            store("Δm", m2 - m1);
+            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 φ) {
+        return R0*φ - R2*sin(2*φ) + R4*sin(4*φ) - R6*sin(6*φ);
+    }
 }
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 3d17a62..aa4ad94 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
@@ -34,6 +34,7 @@ import org.junit.Test;
 
 import static java.lang.StrictMath.*;
 import static org.junit.Assert.*;
+import static org.apache.sis.internal.metadata.ReferencingServices.NAUTICAL_MILE;
 
 
 /**
@@ -152,13 +153,6 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT
                 }
             }
         }
-
-        /** Temporary implementation for allowing tests to pass. */
-        @Override @Deprecated final void computeRhumbLine() {
-            super.computeRhumbLine();
-            rhumblineLength *= ellipsoid.getSemiMajorAxis() / authalicRadius;
-            if (φ1 != 0 || φ2 != 0) rhumblineLength *= 1.1;     // DELETE ME!
-        }
     }
 
     /**
@@ -275,7 +269,7 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT
         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("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);
@@ -511,4 +505,32 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT
         }
         return super.isFailure(expected);
     }
+
+    /**
+     * Tests {@link GeodesicsOnEllipsoid#getRhumblineLength()} using the example given in
Bennett (1996) appendix.
+     */
+    @Test
+    public void testRhumblineLength() {
+        createTracked();
+        verifyParametersForWGS84();
+        /*
+         * Bennett (1996) example 1. For comparing with values given by Bennett:
+         *
+         *   M = Ψ * 10800/PI.
+         *   m = m * semiMajor / 1852
+         */
+        testedEarth.setStartGeographicPoint(10+18.4/60,  37+41.7/60);
+        testedEarth.setEndGeographicPoint  (53+29.5/60, 113+17.1/60);
+        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);
+        assertValueEquals("Δm", 0, 2586.16 / scale,      1E-6, false);
+        assertValueEquals("C",  0, 54.9900,              1E-4, true);
+        assertEquals("distance", 4507.7 * NAUTICAL_MILE, distance, 0.05 * NAUTICAL_MILE);
+    }
 }


Mime
View raw message