sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 03/05: Replace authalic radius by semi-major axis length in GeodeticCalculator. The intent is to allow GeodesicsOnEllipsoid to delegate to its parent class in the equatorial case.
Date Fri, 09 Aug 2019 10:37:07 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 dba1bd6b65acf7a05ff8ce4acd2f745a69bab25f
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Fri Aug 9 11:41:55 2019 +0200

    Replace authalic radius by semi-major axis length in GeodeticCalculator.
    The intent is to allow GeodesicsOnEllipsoid to delegate to its parent class in the equatorial
case.
---
 .../sis/referencing/GeodesicsOnEllipsoid.java      | 25 +++++++++++-----------
 .../apache/sis/referencing/GeodeticCalculator.java | 18 ++++++----------
 .../sis/referencing/GeodesicsOnEllipsoidTest.java  |  7 +++---
 .../sis/referencing/GeodeticCalculatorTest.java    |  8 +++----
 4 files changed, 27 insertions(+), 31 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 0194d91..650b3f0 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
@@ -136,7 +136,9 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
      *
      * @see org.opengis.referencing.datum.Ellipsoid#getSemiMinorAxis()
      */
-    private final double semiMinor;
+    private double semiMinorAxis() {
+        return semiMajorAxis * axisRatio;
+    }
 
     /**
      * The α value computed from the starting point and starting azimuth.
@@ -222,14 +224,13 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
      */
     GeodesicsOnEllipsoid(final CoordinateReferenceSystem crs, final Ellipsoid ellipsoid)
{
         super(crs, ellipsoid);
-        final double a = ellipsoid.getSemiMajorAxis();
+        final double a = semiMajorAxis;
         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;
         /*
          * For rhumb-line calculation from Bennett (1996) equation 2, modified with Clenshaw
summation.
          */
@@ -509,7 +510,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
         /*
          * 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 s2b = I1_σ1 + geodesicDistance / semiMinorAxis();      // (Karney 18)
+ ∆s/b
         final double σ2  = ellipsoidalToSphericalAngle(s2b / A1);           // (Karney 21)
         final double sinσ2 = sin(σ2);
         final double cosσ2 = cos(σ2);
@@ -541,7 +542,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
         φ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("s₂", s2b * semiMinorAxis());
             store("τ₂", s2b / A1);
             store("σ₂", σ2);
             store("α₂", atan2(msinα2, mcosα2));
@@ -637,7 +638,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
                 throw new GeodesicException("Can not compute geodesics for antipodal points
on equator.");
             }
             final double Δφ = φ2 - φ1;
-            geodesicDistance = hypot(Δλ, Δφ) * ellipsoid.getSemiMajorAxis();
+            geodesicDistance = hypot(Δλ, Δφ) * semiMajorAxis;
             msinα1 = msinα2 = (inverseLongitudeSigns ^ swapPoints) ? -Δλ : Δλ;
             mcosα1 = mcosα2 = (inverseLatitudeSigns  ^ swapPoints) ? -Δφ : Δφ;
             setValid(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE);
@@ -797,7 +798,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
                 if (STORE_LOCAL_VARIABLES) {
                     store("J(σ₁)", J1);
                     store("J(σ₂)", J2);
-                    store("Δm",    Δm * semiMinor);
+                    store("Δm",    Δm * semiMinorAxis());
                 }
             }
             final double dα1 = Δλ_error / dΔλ_dα1;                  // Opposite sign
of Karney δα₁ term.
@@ -819,14 +820,14 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
                 store("dλ/dα", dΔλ_dα1);
                 store("δσ₁",  -dα1);
                 store("α₁",    α1);
-                store("s₂",    I1_σ2 * semiMinor);
-                store("Δs",    (I1_σ2 - I1_σ1) * semiMinor);
+                store("s₂",    I1_σ2 * semiMinorAxis());
+                store("Δs",    (I1_σ2 - I1_σ1) * semiMinorAxis());
             }
         } while (moreRefinements != 0);
         final double I1_σ2;
         I1_σ1 = sphericalToEllipsoidalAngle(σ1, false);
         I1_σ2 = sphericalToEllipsoidalAngle(σ2, false);
-        geodesicDistance = (I1_σ2 - I1_σ1) * semiMinor;
+        geodesicDistance = (I1_σ2 - I1_σ1) * semiMinorAxis();
         /*
          * Restore the coordinate sign and order to the original configuration.
          */
@@ -913,7 +914,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
         store("A₃",     A3);
         store("α₀",     atan2(sinα0(), cosα0));
         store("I₁(σ₁)", I1_σ1);
-        store("s₁",     I1_σ1 * semiMinor);
+        store("s₁",     I1_σ1 * semiMinorAxis());
         store("λ₁",     λ1E);
     }
 
@@ -991,7 +992,7 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
                 store("Δm", m2 - m1);
             }
         }
-        rhumblineLength = S * h * (semiMinor / axisRatio);          // TODO: compute semiMajor
only once.
+        rhumblineLength = S * h * semiMajorAxis;
         if (STORE_LOCAL_VARIABLES) {
             store("Δλ", Δλ);
             store("ΔΨ", ΔΨ);
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 48ed79b..e860382 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
@@ -147,13 +147,9 @@ public class GeodeticCalculator {
     final Ellipsoid ellipsoid;
 
     /**
-     * The radius of a hypothetical sphere having the same surface than the {@linkplain #ellipsoid}.
-     * Used for the approximation using spherical formulas.
-     * Subclasses using ellipsoidal formulas will ignore this field.
-     *
-     * @see org.apache.sis.referencing.datum.DefaultEllipsoid#getAuthalicRadius()
+     * Length of the semi-major axis. For a sphere, this is the radius of the sphere.
      */
-    final double authalicRadius;
+    final double semiMajorAxis;
 
     /**
      * The (<var>latitude</var>, <var>longitude</var>) coordinates
of the start point <strong>in radians</strong>.
@@ -255,7 +251,7 @@ public class GeodeticCalculator {
             throw new IllegalArgumentException(Errors.format(Errors.Keys.IllegalCRSType_1,
ReferencingUtilities.getInterface(crs)));
         }
         this.ellipsoid = ellipsoid;
-        authalicRadius = Formulas.getAuthalicRadius(ellipsoid);
+        semiMajorAxis  = ellipsoid.getSemiMajorAxis();
         userToGeodetic = new PositionTransformer(crs, geographic, null);
     }
 
@@ -631,7 +627,7 @@ public class GeodeticCalculator {
             final double ΔΨ = log(tan(PI/4 + φ2/2) / tan(PI/4 + φ1/2));
             factor = Δφ / ΔΨ * hypot(Δλ, ΔΨ);
         }
-        rhumblineLength = authalicRadius * abs(factor);
+        rhumblineLength = semiMajorAxis * abs(factor);
         setValid(RHUMBLINE_LENGTH);
     }
 
@@ -672,7 +668,7 @@ public class GeodeticCalculator {
          */
         double Δσ = sinφ1*sinφ2 + cosφ1*cosφ2*cosΔλ;        // Actually Δσ = acos(…).
         Δσ = atan2(hypot(msinα1, mcosα1), Δσ);              // Δσ ∈ [0…π].
-        geodesicDistance = authalicRadius * Δσ;
+        geodesicDistance = semiMajorAxis * Δσ;
         setValid(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE);
     }
 
@@ -704,7 +700,7 @@ public class GeodeticCalculator {
     void computeEndPoint() throws GeodesicException {
         canComputeEndPoint();                   // May throw IllegalStateException.
         final double m     = hypot(msinα1, mcosα1);
-        final double Δσ    = geodesicDistance / authalicRadius;
+        final double Δσ    = geodesicDistance / semiMajorAxis;
         final double sinΔσ = sin(Δσ);
         final double cosΔσ = cos(Δσ);
         final double sinφ1 = sin(φ1);
@@ -873,7 +869,7 @@ public class GeodeticCalculator {
             super(ReferencingUtilities.getDimension(userToGeodetic.defaultCRS));
             msinαf = msinα2;  λf = λ2;
             mcosαf = mcosα2;  φf = φ2;
-            tolerance = toDegrees(εx / authalicRadius);
+            tolerance = toDegrees(εx / semiMajorAxis);
             distance  = geodesicDistance;
             length    = rhumblineLength;
             flags     = validity;
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 f52f0b0..7b28707 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
@@ -199,10 +199,9 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT
      * <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("a",    6378137,             testedEarth.semiMajorAxis,            
       STRICT);
         assertEquals("b",    6356752.314245,      testedEarth.ellipsoid.getSemiMinorAxis(),
    1E-6);
-        assertEquals("c",    6371007.180918,      testedEarth.authalicRadius,           
       1E-6);
+        assertEquals("1/f",  298.257223563,       testedEarth.ellipsoid.getInverseFlattening(),
1E-9);
         assertEquals("ℯ²",   0.00669437999014132, testedEarth.eccentricitySquared,   
          1E-17);
         assertEquals("ℯ′²",  0.00673949674227643, testedEarth.secondEccentricitySquared,
       1E-17);
         assertEquals("n",    0.00167922038638370, testedEarth.thirdFlattening,          
       1E-16);
@@ -523,7 +522,7 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT
         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;
+        final double scale = testedEarth.semiMajorAxis / NAUTICAL_MILE;
         assertValueEquals("Δλ", 0, 75+35.4 / 60,         1E-11, true);
         assertValueEquals("ΔΨ", 0, 3176.89 / (10800/PI), 1E-5, false);
         assertValueEquals("m₁", 0,  615.43 / scale,      1E-6, false);
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 11413a2..cd49677 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
@@ -116,7 +116,7 @@ public strictfp class GeodeticCalculatorTest extends TestCase {
      * 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());
+        return new Geodesic(c.semiMajorAxis, 1/c.ellipsoid.getInverseFlattening());
     }
 
     /**
@@ -408,7 +408,7 @@ public strictfp class GeodeticCalculatorTest extends TestCase {
         final Statistics   yError   = new Statistics("∆y/r");
         final Statistics   aErrors  = new Statistics("∆α (°)");
         final double       azimuth  = c.getStartingAzimuth();
-        final double       toMetres = (PI/180) * c.authalicRadius;
+        final double       toMetres = (PI/180) * c.semiMajorAxis;
         final double[]     buffer   = new double[2];
         while (!iterator.isDone()) {
             switch (iterator.currentSegment(buffer)) {
@@ -478,7 +478,7 @@ public strictfp class GeodeticCalculatorTest extends TestCase {
                      * especially near poles and between antipodal points. Following are
empirical thresholds.
                      */
                     linearTolerance    = expected[COLUMN_Δs] * 0.01;
-                    latitudeTolerance  = toDegrees(linearTolerance / c.authalicRadius);
+                    latitudeTolerance  = toDegrees(linearTolerance / c.semiMajorAxis);
                     longitudeTolerance = expected[COLUMN_φ2] > 89.5 ? 180 : latitudeTolerance
/ cosφ2;
                     azimuthTolerance   = 0.5;                                   // About
8.8 metres at distance of 1 km.
                     if (isTestingInverse) {
@@ -492,7 +492,7 @@ public strictfp class GeodeticCalculatorTest extends TestCase {
                 } else {
                     /*
                      * When ellipsoidal formulas are used, we aim for an 1 mm accuracy in
coordinate values.
-                     * We also aim for azimuthd such as the error is less than 1 cm after
the first 10 km.
+                     * We also aim for azimuths 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.
                      */
                     linearTolerance    = 2 * GeodesicsOnEllipsoid.ITERATION_TOLERANCE * AUTHALIC_RADIUS;


Mime
View raw message