You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sis.apache.org by de...@apache.org on 2019/08/05 15:32:42 UTC
[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
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 <ma...@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);
+ }
}