sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 02/03: Apply to EqualAreaProjection the same Snyder 3-35 formulas (trigonometric identities) than we did for MeridianArcBased. Include coefficients and other fields in serialization for the same reason than ConformalProjection (simplicity, stability). During this work, we identified a bug in the series expansion of φ(β) where the last term was multiplied by e⁸ instead of e⁶. Fixing that bug bring the series expansion to a precision where we no longer need iterations at least with WGS84 ellipsoid.
Date Mon, 29 Apr 2019 16:34: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 6760894f16bca4bfcfbd9e12c8d0daf914b9a089
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Mon Apr 29 14:10:24 2019 +0200

    Apply to EqualAreaProjection the same Snyder 3-35 formulas (trigonometric identities)
than we did for MeridianArcBased.
    Include coefficients and other fields in serialization for the same reason than ConformalProjection
(simplicity, stability).
    During this work, we identified a bug in the series expansion of φ(β) where the last
term was multiplied by e⁸ instead of e⁶.
    Fixing that bug bring the series expansion to a precision where we no longer need iterations
at least with WGS84 ellipsoid.
---
 .../operation/projection/AlbersEqualArea.java      |   1 -
 .../operation/projection/ConformalProjection.java  |  14 +-
 .../operation/projection/CylindricalEqualArea.java |   1 -
 .../operation/projection/EqualAreaProjection.java  | 202 ++++++++++-----------
 .../projection/EqualAreaProjectionTest.java        | 127 +++++++++++++
 .../sis/test/suite/ReferencingTestSuite.java       |   1 +
 6 files changed, 227 insertions(+), 119 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
index c2a137a..c7de81a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
@@ -161,7 +161,6 @@ public class AlbersEqualArea extends EqualAreaProjection {
         denormalize.convertBefore(0, rn, null); rn.negate();
         denormalize.convertBefore(1, rn, ρ0);   rn.inverseDivide(-1);
         normalize.convertAfter(0, rn, null);
-        super.computeCoefficients();
     }
 
     /**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
index 35f0151..c2a9045 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
@@ -23,10 +23,9 @@ import static java.lang.Math.*;
 
 /**
  * Base class of {@link LambertConicConformal}, {@link Mercator} and {@link PolarStereographic}
projections.
- * All those projections have in common the property of being <cite>conformal</cite>,
i.e. they preserve
- * angles locally. However we do not put this base class in public API because not all conformal
projections
- * will extend this base class. For example {@link TransverseMercator}, despite being conformal,
uses very
- * different formulas.
+ * Those projections have in common the property of being <cite>conformal</cite>,
i.e. they preserve angles locally.
+ * However we do not put this base class in public API because not all conformal projections
extend this base class.
+ * For example the {@link TransverseMercator} projection, despite being conformal, uses very
different formulas.
  *
  * <p>Note that no projection can be both conformal and equal-area. This restriction
is implemented in class
  * hierarchy with {@link ConformalProjection} and {@link EqualAreaProjection} being two distinct
classes.</p>
@@ -85,8 +84,9 @@ abstract class ConformalProjection extends NormalizedProjection {
     private final boolean useIterations;
 
     /**
-     * Coefficients in the series expansion of the inverse projection, depending only on
{@linkplain #eccentricity
-     * eccentricity} value. The series expansion is published under the following form, where
χ is the conformal latitude:
+     * Coefficients of the first terms in the series expansion of the inverse projection.
+     * Values of those coefficients depend only on {@linkplain #eccentricity eccentricity}
value.
+     * The series expansion is published under the following form, where χ is the <cite>conformal
latitude</cite>:
      *
      *     <blockquote>c₂⋅sin(2χ) + c₄⋅sin(4χ) + c₆⋅sin(6χ) + c₈⋅sin(8χ)</blockquote>
      *
@@ -200,7 +200,7 @@ abstract class ConformalProjection extends NormalizedProjection {
          * Note that the φ value computed by the line below is called χ in EPSG guide.
          * We name it φ in our code because we will modify that value in-place in order
to get φ.
          */
-        double φ = (PI/2) - 2*atan(expOfSouthing);          // at this point == χ
+        double φ = (PI/2) - 2*atan(expOfSouthing);          // at this point == χ (conformal
latitude).
         /*
          * Add a correction for the flattened shape of the Earth. The correction can be represented
by an
          * infinite series. Here, we apply only the first 4 terms. Those terms are given
by §1.3.3 in the
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
index 6cc9fb6..7fda325 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
@@ -168,7 +168,6 @@ public class CylindricalEqualArea extends EqualAreaProjection {
         ik.divide(k0);
         denormalize.convertAfter(0, k0, null);
         denormalize.convertAfter(1, ik, null);
-        super.computeCoefficients();
     }
 
     /**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java
index 26c6579..3999e54 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java
@@ -16,8 +16,6 @@
  */
 package org.apache.sis.referencing.operation.projection;
 
-import java.io.IOException;
-import java.io.ObjectInputStream;
 import org.apache.sis.internal.referencing.Resources;
 
 import static java.lang.Math.*;
@@ -25,10 +23,16 @@ import static org.apache.sis.math.MathFunctions.atanh;
 
 
 /**
- * Provides formulas common to Equal Area projections.
+ * Base class of {@link AlbersEqualArea} and {@link CylindricalEqualArea} projections.
+ * Those projections have in common the property of being <cite>equal-area</cite>.
+ * However we do not put this base class in public API because not all equal-area projections
extend this base class.
+ * For example the {@link Sinusoidal} projection, despite being equal-area, uses different
formulas.
+ *
+ * <p>Note that no projection can be both conformal and equal-area. This restriction
is implemented in class
+ * hierarchy with {@link ConformalProjection} and {@link EqualAreaProjection} being two distinct
classes.</p>
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 0.8
+ * @version 1.0
  * @since   0.8
  * @module
  */
@@ -36,56 +40,48 @@ abstract class EqualAreaProjection extends NormalizedProjection {
     /**
      * For cross-version compatibility.
      */
-    private static final long serialVersionUID = -6175270149094989517L;
+    private static final long serialVersionUID = 2221537810038553082L;
 
     /**
-     * {@code false} for using the original formulas as published by Snyder, or {@code true}
for using formulas
-     * modified using trigonometric identities. The use of trigonometric identities is for
reducing the amount
-     * of calls to the {@link Math#sin(double)} and similar methods. Some identities used
are:
-     *
-     * <ul>
-     *   <li>sin(2β) = 2⋅sinβ⋅cosβ</li>
-     *   <li>sin(4β) = (2 - 4⋅sin²β)⋅sin(2β)</li>
-     *   <li>sin(8β) = 4⋅sin(2β)⋅(cos²β - sin²β)⋅(8⋅cos⁴β - 8⋅cos²β
+ 1)</li>
-     * </ul>
-     *
-     * Note that since this boolean is static final, the compiler should exclude the code
in the branch that is never
-     * executed (no need to comment-out that code).
+     * The threshold value of {@link #eccentricity} at which we consider the accuracy of
the
+     * series expansion insufficient. This threshold is determined empirically with the help
+     * of the {@code EqualAreaProjectionTest.searchThreshold()} method in the test directory.
      *
-     *
-     * <p><b>BENCHMARK AND ANALYSIS:</b>
-     * as of July 2016, benchmarking shows small benefit in using trigonometric identities
for {@code EqualAreaProjection}
-     * (contrarily to {@link ConformalProjection} where we measured a greater benefit). This
may be because in this class,
-     * the series expansion is unconditionally followed by iterative method in order to reach
the centimetric precision.
-     * We observe that the original series expansion allows convergence in only one iteration,
while the formulas using
-     * trigonometric identifies often requires two iterations.</p>
-     *
-     * @todo Redo the benchmark with JDK9, since {@code sin} seems much faster.
+     * <p>Note: WGS84 eccentricity is about 0.082.</p>
      */
-    private static final boolean ALLOW_TRIGONOMETRIC_IDENTITIES = true;
+    private static final double ECCENTRICITY_THRESHOLD = 0.1;
 
     /**
-     * Coefficients in the series expansion of the inverse projection,
-     * depending only on {@linkplain #eccentricity eccentricity} value.
-     * The series expansion is of the following form:
+     * Coefficients of the first terms in the series expansion of the inverse projection.
+     * Values of those coefficients depend only on {@linkplain #eccentricity eccentricity}
value.
+     * The series expansion is published under the following form, where β is the <cite>authalic
latitude</cite>:
      *
-     *     <blockquote>φ = ci₂⋅sin(2β) + ci₄⋅sin(4β) + ci₈⋅sin(8β)</blockquote>
+     *     <blockquote>φ = c₂⋅sin(2β) + c₄⋅sin(4β) + c₈⋅sin(6β)</blockquote>
      *
-     * This {@code EqualAreaProjection} class uses those coefficients in {@link #φ(double)}.
+     * However we rewrite above series expansion for taking advantage of trigonometric identities.
+     * The equation become (with different <var>c</var> values):
      *
-     * <p><strong>Consider those fields as final!</strong> They are not
final only for sub-class
-     * constructors convenience and for the purpose of {@link #readObject(ObjectInputStream)}.</p>
+     *     <blockquote>sin(2β)⋅(c₂ + cos(2β)⋅(c₄ + cos(2β)⋅c₆))</blockquote>
      *
-     * @see #computeCoefficients()
+     * <div class="note"><b>Serialization note:</b>
+     * we do not strictly need to serialize those fields since they could be computed after
deserialization.
+     * Bu we serialize them anyway in order to simplify a little bit this class (it allows
us to keep those
+     * fields final) and because values computed after deserialization could be slightly
different than the
+     * ones computed after construction if a future version of the constructor uses the double-double
values
+     * provided by {@link Initializer}.</div>
      */
-    private transient double ci2, ci4, ci8;
+    private final double c2β, c4β, c6β;
 
     /**
      * Value of {@link #qm(double)} function (part of Snyder equation (3-12)) at pole (sinφ
= 1).
-     *
-     * @see #computeCoefficients()
      */
-    private transient double qmPolar;
+    private final double qmPolar;
+
+    /**
+     * {@code true} if the {@link #eccentricity} value is greater than or equals to the eccentricity
threshold,
+     * in which case the {@link #φ(double)} method will need to use an iterative method.
+     */
+    private final boolean useIterations;
 
     /**
      * Creates a new normalized projection from the parameters computed by the given initializer.
@@ -94,31 +90,37 @@ abstract class EqualAreaProjection extends NormalizedProjection {
      */
     EqualAreaProjection(final Initializer initializer) {
         super(initializer);
-    }
-
-    /**
-     * Computes the coefficients in the series expansions from the {@link #eccentricitySquared}
value.
-     * This method shall be invoked after {@code EqualAreaProjection} construction or deserialization.
-     */
-    void computeCoefficients() {
         final double e2 = eccentricitySquared;
-        final double e4  = e2 * e2;
-        final double e6  = e2 * e4;
-        ci2  =  517/5040.  * e6  +  31/180. * e4  +  1/3. * e2;
-        ci4  =  251/3780.  * e6  +  23/360. * e4;
-        ci8  =  761/45360. * e6;
+        final double e4 = e2 * e2;
+        final double e6 = e2 * e4;
         /*
-         * When rewriting equations using trigonometric identities, some constants appear.
-         * For example sin(2β) = 2⋅sinβ⋅cosβ, so we can factor out the 2 constant
into the
-         * into the corresponding 'c' field.
+         * Published equation is of the following form:
+         *
+         *     φ  =  β + A⋅sin(2β) + B⋅sin(4β) + C⋅sin(6β)                  (part
of Snyder 3-34)
+         *
+         * We rewrite as:
+         *
+         *     φ  =  β + sin(2β)⋅(A′ + cos(2β)⋅(B′ + cos(2β)⋅C′))      
    (part of Snyder 3-35)
+         *     A′ =  A - C
+         *     B′ =  2B - 4D
+         *     C′ =  4C
+         *
+         * Published coefficients are:
+         *
+         *     A  =  517/5040  ⋅ℯ⁶  +  31/180 ⋅ℯ⁴  +  ℯ²/3
+         *     B  =  251/3780  ⋅ℯ⁶  +  23/360 ⋅ℯ⁴
+         *     C  =  761/45360 ⋅ℯ⁶
+         *
+         * We replace A B C D by A′ B′ C′ D′ using Snyder 3-35, then replace sin2β
by 2sinβ⋅cosβ
+         * and cos2β by 1 – 2sin²β. The result are coefficients below. For each line,
we add the
+         * smallest values first in order to reduce rounding errors. The algebraic changes
are
+         * in https://svn.apache.org/repos/asf/sis/analysis/Map%20projection%20formulas.ods.
          */
-        if (ALLOW_TRIGONOMETRIC_IDENTITIES) {
-            // Multiplication by powers of 2 does not bring any additional rounding error.
-            ci2 *=  2;
-            ci4 *=  8;
-            ci8 *= 64;
-        }
+        c2β     =      4./7    * e6  +    3./5  * e4  +  2./3 * e2;
+        c4β     =  -3028./2835 * e6  +  -23./45 * e4;
+        c6β     =   1522./2835 * e6;
         qmPolar = qm(1);
+        useIterations = (eccentricity >= ECCENTRICITY_THRESHOLD);
     }
 
     /**
@@ -129,10 +131,11 @@ abstract class EqualAreaProjection extends NormalizedProjection {
      */
     EqualAreaProjection(final EqualAreaProjection other) {
         super(other);
-        ci2 = other.ci2;
-        ci4 = other.ci4;
-        ci8 = other.ci8;
+        c2β     = other.c2β;
+        c4β     = other.c4β;
+        c6β     = other.c6β;
         qmPolar = other.qmPolar;
+        useIterations = other.useIterations;
     }
 
     /**
@@ -204,49 +207,36 @@ abstract class EqualAreaProjection extends NormalizedProjection {
      */
     final double φ(final double y) throws ProjectionException {
         final double sinβ = y / qmPolar;
+        final double sinβ2 = sinβ * sinβ;
         final double β = asin(sinβ);
-        double φ;
-        if (!ALLOW_TRIGONOMETRIC_IDENTITIES) {
-            φ = ci8 * sin(8*β)
-              + ci4 * sin(4*β)
-              + ci2 * sin(2*β)
-              + β;                                                                  // Snyder
3-18
-        } else {
-            /*
-             * Same formula than above, but rewriten using trigonometric identities in order
to avoid
-             * multiple calls to sin(double) method. The cost is only one sqrt(double) method
call.
-             */
-            final double sin2_β = sinβ*sinβ;                                        //
= sin²β
-            final double cos2_β = 1 - sin2_β;                                       //
= cos²β
-            final double t2β = sinβ * sqrt(cos2_β);                                 //
= sin(2β) /   2
-            final double t4β = 0.5 - sin2_β;                                        //
= sin(4β) / ( 4⋅sin(2β))
-            final double t8β = (cos2_β - sin2_β)*(cos2_β*cos2_β - cos2_β + 1./8); 
 // = sin(8β) / (32⋅sin(2β))
-
-            φ = (ci8*t8β  +  ci4*t4β  +  ci2) * t2β  +  β;
-        }
         /*
-         * At this point φ is close to the desired value, but may have an error of a few
centimetres.
-         * Use the iterative method for reaching the last part of missing accuracy. Usually
this loop
-         * will perform exactly one iteration, no more, because φ is already quite close
to the result.
-         *
-         * Mathematical note: Snyder 3-16 gives q/(1-ℯ²) instead of y in the calculation
of Δφ below.
-         * For Cylindrical Equal Area projection, Snyder 10-17 gives  q = (qPolar⋅sinβ),
which simplifies
-         * as y.
-         *
-         * For Albers Equal Area projection, Snyder 14-19 gives  q = (C - ρ²n²/a²)/n,
 which we rewrite
-         * as  q = (C - ρ²)/n  (see comment in AlbersEqualArea.inverseTransform(…) for
the mathematic).
-         * The y value given to this method is y = (C - ρ²) / (n⋅(1-ℯ²)) = q/(1-ℯ²),
the desired value.
+         * Snyder 3-18, but rewriten using trigonometric identities in order to avoid
+         * multiple calls to sin(double) method.
          */
-        for (int i=0; i<MAXIMUM_ITERATIONS; i++) {
-            final double sinφ  = sin(φ);
-            final double cosφ  = cos(φ);
-            final double ℯsinφ = eccentricity * sinφ;
-            final double ome   = 1 - ℯsinφ*ℯsinφ;
-            final double Δφ    = ome*ome/(2*cosφ) * (y - sinφ/ome - atanh(ℯsinφ)/eccentricity);
-            φ += Δφ;
-            if (abs(Δφ) <= ITERATION_TOLERANCE) {
-                return φ;
+        double φ = β + cos(β)*sinβ*(c2β + sinβ2*(c4β + sinβ2*c6β));
+        if (useIterations) {
+            /*
+             * Mathematical note: Snyder 3-16 gives q/(1-ℯ²) instead of y in the calculation
of Δφ below.
+             * For Cylindrical Equal Area projection, Snyder 10-17 gives  q = (qPolar⋅sinβ),
which simplifies
+             * as y.
+             *
+             * For Albers Equal Area projection, Snyder 14-19 gives  q = (C - ρ²n²/a²)/n,
 which we rewrite
+             * as  q = (C - ρ²)/n  (see comment in AlbersEqualArea.inverseTransform(…)
for the mathematic).
+             * The y value given to this method is y = (C - ρ²) / (n⋅(1-ℯ²)) = q/(1-ℯ²),
the desired value.
+             */
+            for (int i=0; i<MAXIMUM_ITERATIONS; i++) {
+                final double sinφ  = sin(φ);
+                final double cosφ  = cos(φ);
+                final double ℯsinφ = eccentricity * sinφ;
+                final double ome   = 1 - ℯsinφ*ℯsinφ;
+                final double Δφ    = ome*ome/(2*cosφ) * (y - sinφ/ome - atanh(ℯsinφ)/eccentricity);
+                φ += Δφ;
+                if (abs(Δφ) <= ITERATION_TOLERANCE) {
+                    return φ;
+                }
             }
+        } else if (!Double.isNaN(φ)) {
+            return φ;
         }
         /*
          * In the Albers Equal Area discussion, Snyder said that above algorithm does not
converge if
@@ -275,12 +265,4 @@ abstract class EqualAreaProjection extends NormalizedProjection {
         // Value should have converged but did not.
         throw new ProjectionException(Resources.format(Resources.Keys.NoConvergence));
     }
-
-    /**
-     * Restores transient fields after deserialization.
-     */
-    private void readObject(final ObjectInputStream in) throws IOException, ClassNotFoundException
{
-        in.defaultReadObject();
-        computeCoefficients();
-    }
 }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/EqualAreaProjectionTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/EqualAreaProjectionTest.java
new file mode 100644
index 0000000..40953a4
--- /dev/null
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/EqualAreaProjectionTest.java
@@ -0,0 +1,127 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.referencing.operation.projection;
+
+import java.util.Random;
+import org.apache.sis.parameter.Parameters;
+import org.apache.sis.internal.referencing.provider.LambertCylindricalEqualArea;
+import org.apache.sis.referencing.operation.DefaultOperationMethod;
+import org.apache.sis.internal.util.Constants;
+import org.apache.sis.test.TestUtilities;
+import org.apache.sis.test.DependsOn;
+import org.junit.Test;
+
+import static java.lang.StrictMath.*;
+import static org.apache.sis.math.MathFunctions.atanh;
+import static org.junit.Assert.assertEquals;
+
+
+/**
+ * Tests {@link EqualAreaProjection}.
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.0
+ * @since   1.0
+ * @module
+ */
+@DependsOn(NormalizedProjectionTest.class)
+public final strictfp class EqualAreaProjectionTest extends MapProjectionTestCase {
+    /**
+     * Creates the projection to be tested.
+     *
+     * @param  ellipsoidal   {@code false} for a sphere, or {@code true} for WGS84 ellipsoid.
+     * @return a test instance of the projection.
+     */
+    private EqualAreaProjection create(final boolean ellipsoidal) {
+        final DefaultOperationMethod provider = new LambertCylindricalEqualArea();
+        final CylindricalEqualArea projection = new CylindricalEqualArea(provider, parameters(provider,
ellipsoidal));
+        tolerance = NormalizedProjection.ANGULAR_TOLERANCE;     // = linear tolerance on
a sphere of radius 1.
+        return projection;
+    }
+
+    /**
+     * Computes φ using equation given in EPSG guidance notes, which is also from Snyder
book.
+     * We use this equation as a reference for testing validity of other forms.
+     *
+     * @param  y  in the cylindrical case, this is northing on the normalized ellipsoid.
+     * @return the latitude in radians.
+     */
+    private static double reference(final EqualAreaProjection projection, final double y)
{
+        final double e    = projection.eccentricity;
+        final double e2   = projection.eccentricitySquared;
+        final double e4   = e2 * e2;
+        final double e6   = e2 * e4;
+        final double c2β  = 517./5040  * e6  +  31./180 * e4  +  1./3 * e2;
+        final double c4β  = 251./3780  * e6  +  23./360 * e4;
+        final double c6β  = 761./45360 * e6;
+        final double qmp  = (1/(1 - e*e) + atanh(e)/e);
+        final double sinβ = y / qmp;
+        final double β    = asin(sinβ);
+        return c6β * sin(6*β)
+             + c4β * sin(4*β)
+             + c2β * sin(2*β)
+             + β;                           // Snyder 3-18
+    }
+
+    /**
+     * Compares {@link EqualAreaProjection#φ(double)} with formula taken as references.
+     *
+     * @throws ProjectionException if the function does not converge.
+     */
+    @Test
+    public void compareWithReference() throws ProjectionException {
+        final EqualAreaProjection projection = create(true);
+        final Random random = TestUtilities.createRandomNumberGenerator();
+        for (int i=0; i<100; i++) {
+            final double y = random.nextDouble() * 3 - 1.5;
+            final double reference = reference(projection, y);
+            final double actual    = projection.φ(y);
+            assertEquals(reference, actual, NormalizedProjection.ITERATION_TOLERANCE);
+        }
+    }
+
+    /**
+     * Searches a value for {@link EqualAreaProjection#ECCENTRICITY_THRESHOLD}.
+     * This method is not part of test suite. Steps to enable:
+     *
+     * <ol>
+     *   <li>In {@link EqualAreaProjection#φ(double)} method, for {@code useIterations}
to {@code false}.</li>
+     *   <li>Add a {@link Test} annotation on this method.
+     * </ol>
+     *
+     * @throws ProjectionException if {@link EqualAreaProjection#φ(double)} did not converge.
+     */
+    @SuppressWarnings("UseOfSystemOutOrSystemErr")
+    public void searchThreshold() throws ProjectionException {
+        tolerance = NormalizedProjection.ANGULAR_TOLERANCE;
+        final DefaultOperationMethod provider = new LambertCylindricalEqualArea();
+        final Parameters parameters = parameters(provider, true);
+        for (double e = 0.05; e <= 0.2; e += 0.001) {
+            final double a = parameters.parameter(Constants.SEMI_MAJOR).doubleValue();
+            parameters.parameter(Constants.SEMI_MINOR).setValue(a * sqrt(1 - e*e));
+            final CylindricalEqualArea projection = new CylindricalEqualArea(provider, parameters);
+            for (double y = -1.25; y <= 1.25; y += 0.01) {
+                final double reference = reference(projection, y);
+                final double actual = projection.φ(y);
+                if (abs(actual - reference) > NormalizedProjection.ANGULAR_TOLERANCE)
{
+                    System.out.println("Error exceeds tolerance threshold at eccentricity
" + e);
+                    return;
+                }
+            }
+        }
+    }
+}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
index 43e533b..85d9c85 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
@@ -173,6 +173,7 @@ import org.junit.BeforeClass;
     org.apache.sis.referencing.operation.projection.PolarStereographicTest.class,
     org.apache.sis.referencing.operation.projection.ObliqueStereographicTest.class,
     org.apache.sis.referencing.operation.projection.ObliqueMercatorTest.class,
+    org.apache.sis.referencing.operation.projection.EqualAreaProjectionTest.class,
     org.apache.sis.referencing.operation.projection.CylindricalEqualAreaTest.class,
     org.apache.sis.referencing.operation.projection.AlbersEqualAreaTest.class,
     org.apache.sis.referencing.operation.projection.MeridianArcTest.class,


Mime
View raw message