sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 01/02: Apply some algebraic rearrangements and update the documentation about which Snyder formulas are used. Verify the validity of latitude arguments compared to satellite orbit inclination. Tune the tolerance factors used in tests.
Date Fri, 04 Oct 2019 15:07:57 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 130cec114323d2bf0ea412dfb6fb4975e88daca4
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Fri Oct 4 14:18:09 2019 +0200

    Apply some algebraic rearrangements and update the documentation about which Snyder formulas
are used.
    Verify the validity of latitude arguments compared to satellite orbit inclination.
    Tune the tolerance factors used in tests.
---
 .../operation/projection/SatelliteTracking.java    | 304 ++++++++++++---------
 .../projection/SatelliteTrackingTest.java          | 142 ++++++----
 2 files changed, 263 insertions(+), 183 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/SatelliteTracking.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/SatelliteTracking.java
index d678cec..cbd9a5a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/SatelliteTracking.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/SatelliteTracking.java
@@ -20,12 +20,15 @@ import java.util.EnumMap;
 import org.opengis.parameter.ParameterDescriptor;
 import org.opengis.referencing.operation.Matrix;
 import org.opengis.referencing.operation.OperationMethod;
+import org.opengis.parameter.InvalidParameterValueException;
 import org.apache.sis.internal.referencing.Formulas;
 import org.apache.sis.internal.referencing.Resources;
 import org.apache.sis.referencing.operation.matrix.Matrix2;
 import org.apache.sis.referencing.operation.matrix.MatrixSIS;
 import org.apache.sis.referencing.operation.transform.ContextualParameters;
 import org.apache.sis.parameter.Parameters;
+import org.apache.sis.measure.Latitude;
+import org.apache.sis.util.resources.Errors;
 import org.apache.sis.util.Workaround;
 
 import static java.lang.Math.*;
@@ -110,8 +113,9 @@ public class SatelliteTracking extends NormalizedProjection {
     /**
      * Creates a Satellite Tracking projection from the given parameters.
      *
-     * @param method      description of the projection parameters.
-     * @param parameters  the parameter values of the projection to create.
+     * @param  method      description of the projection parameters.
+     * @param  parameters  the parameter values of the projection to create.
+     * @throws InvalidParameterValueException if some parameters have incompatible values.
      */
     public SatelliteTracking(final OperationMethod method, final Parameters parameters) {
         this(initializer(method, parameters));
@@ -126,132 +130,188 @@ public class SatelliteTracking extends NormalizedProjection {
         final double φ0 = toRadians(initializer.getAndStore(LATITUDE_OF_ORIGIN));
         final double φ1 = toRadians(initializer.getAndStore(STANDARD_PARALLEL_1));
         final double φ2 = toRadians(initializer.getAndStore(STANDARD_PARALLEL_2));
-        isConic = Math.abs(φ2 + φ1) > ANGULAR_TOLERANCE;
+        final double i  = toRadians(initializer.getAndStore(SATELLITE_ORBIT_INCLINATION));
+        p2_on_p1 =                  initializer.getAndStore(SATELLITE_ORBITAL_PERIOD) /
+                                    initializer.getAndStore(ASCENDING_NODE_PERIOD);
         /*
-         * Common for both cylindrical and conic sattelite tracking projections.
-         * Symbols from Snyder:
+         * Symbols from Snyder used in the constructor:
          *
          *   i  =  angle of inclination between the plane of Earth Equator and the plane
of satellite orbit.
          *   P₁ =  length of Earth's rotation with respect to precessed ascending node.
          *   P₂ =  time required for revolution of the satellite.
-         *   φ₁ =  standard parallel (North and South in cylindrical case).
-         *   φ₂ =  second standard parallel (conic case only).
+         *   φ₁ =  standard parallel.
+         *   φ₂ =  second standard parallel. Equals to -φ₁ in the cylindrical case.
+         *   F₁ =  angle on the globe between the ground-track and the meridian at latitude
φ₁.
+         *   n  :  used for conic projection only.
+         *   s₀ :  user for conic projection only.
+         *
+         * Next terms below are common to both cylindrical and conic sattelite tracking projections.
          */
-        final double i = toRadians(initializer.getAndStore(SATELLITE_ORBIT_INCLINATION));
-        cos_i    = cos(i);
-        sin_i    = sin(i);
-        cos2_i   = cos_i * cos_i;
-        p2_on_p1 = initializer.getAndStore(SATELLITE_ORBITAL_PERIOD) / initializer.getAndStore(ASCENDING_NODE_PERIOD);
-        final double cos_φ1  = cos(φ1);
-        final double cos2_φ1 = cos_φ1 * cos_φ1;
+        sin_i   = sin(i);
+        cos_i   = cos(i);
+        cos2_i  = cos_i * cos_i;
+        isConic = Math.abs(φ2 + φ1) > ANGULAR_TOLERANCE;
+        final double cosφ1   = cos(φ1);
+        final double cos2_φ1 = cosφ1 * cosφ1;
         /*
-         *
+         * Some terms will be stored as scale factor and offset applied by matrix before
projection
+         * (normalization) and after projection (denormalization).  Those factors will be
different
+         * for cylindric and conic projections.
          */
-        double scale;
-        final Double ρ0;
+        final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
         if (isConic) {
-            final double sin_φ1 = sin(φ1);
+            final double sinφ1 = sin(φ1);
+            final double L0 =  L(sin(φ0), LATITUDE_OF_ORIGIN);
+            final double L1 =  L(sinφ1,   STANDARD_PARALLEL_1);
+            final double F1 =  F(cos2_φ1, STANDARD_PARALLEL_1);
             /*
-             * Conic projection case.
-             * Inclination of the groundtrack to the meridian at latitude φ1
+             * For conic projection with one standard parallel (φ₁ = φ₂), the general
case implemented by
+             * equation 28-10 become indeterminate. Equation 28-17 below resolves that indetermination.
+             * If furthermore φ₁ is equal to the upper tracking limit π-i, that equation
28-17 could be
+             * simplified by equation 28-18:
+             *
+             *     f = p2_on_p1 * cos_i - 1;
+             *     n = sin_i / (f*f);                                       // Snyder equation
28-18.
+             *
+             * However since equation 28-17 still work, we keep it for avoiding discontinuity.
              */
-            final double F1  = computeFn(cos2_φ1);
-            final double dλ0 = -asin(sin(φ0) / sin_i);       // eq.28-2a in Snyder
-            final double dλ1 = -asin(sin_φ1  / sin_i);
-            final double λt0 = computeλtn(dλ0);
-            final double λt1 = computeλtn(dλ1);
-            final double L0  = λt0 - p2_on_p1 * dλ0;
-            final double L1  = λt1 - p2_on_p1 * dλ1;
-
-            //tracking limit computed as 180 - i from Snyder's manual p.238
-            if (φ1 == PI - i) {
-                final double factor = (p2_on_p1 * cos_i - 1);
-                n = sin_i / (factor * factor); //eq. 28-18 in Snyder
-            } else if (φ2 != φ1) {
-                final double cos_φ2 = cos(φ2);
-                final double dλ2 = -asin(sin(φ2) / sin_i);
-                final double λt2 = computeλtn(dλ2);
-                final double L2 = λt2 - p2_on_p1 * dλ2;
-                final double F2 = computeFn(cos_φ2 * cos_φ2);
-                n = (F2 - F1) / (L2 - L1);
+            if (abs(φ2 - φ1) < ANGULAR_TOLERANCE) {
+                n = sinφ1 * (p2_on_p1 * (2*cos2_i - cos2_φ1) - cos_i)
+                          / (p2_on_p1 * (           cos2_φ1) - cos_i);      // Equation
28-17.
             } else {
-                n = sin_φ1 * (p2_on_p1 * (2 * cos2_i - cos2_φ1) - cos_i) / (p2_on_p1 *
cos2_φ1 - cos_i); //eq. 28-17 in Snyder
+                final double cosφ2 = cos(φ2);
+                final double F2 = F(cosφ2*cosφ2, STANDARD_PARALLEL_2);
+                final double L2 = L(sin(φ2),     STANDARD_PARALLEL_2);
+                n = (F2 - F1) / (L2 - L1);                                  // Snyder equation
28-10.
+            }
+            s0 = F1 - n*L1;                                                 // Snyder equation
28-11.
+            final double ρf = cosφ1 * sin(F1) / n;                          // Part of
Snyder 28-12 fraction.
+            final double ρ0 = ρf / sin(n*L0 + s0);                          // Remaining
of Snyder 28-12 without R.
+            if (!Double.isFinite(ρf) || ρf == 0) {
+                throw invalid(STANDARD_PARALLEL_1);
             }
-            // cos(φ₁) × sin(F₁) / n
-            scale = cos_φ1 * sin(F1) / n;
-            s0 = F1 - n * L1;
-            ρ0 = scale / sin(n * L0 + s0); // *R in eq.28-12 in Snyder
-            scale = -scale;
+            final MatrixSIS normalize = context.getMatrix(ContextualParameters.MatrixRole.NORMALIZATION);
+            normalize  .convertAfter (0,  n,  null);
+            denormalize.convertBefore(0, +ρf, null);
+            denormalize.convertBefore(1, -ρf, ρ0);
         } else {
             /*
-             * Cylindrical projection case.
+             * Cylindrical projection case. The equations are (ignoring R and λ₀):
+             *
+             *     x   =  λ⋅cosφ₁
+             *     y   =  L⋅cosφ₁/F₁′  where  L depends on φ but F₁′ is constant.
+             *     F₁′ =  tan(F₁)
+             *
+             * The cosφ₁ (for x at dimension 0) and cosφ₁/F₁′ (for y at dimension
1) factors are computed
+             * in advance and stored below. The remaining factor to compute in transform(…)
method is L.
              */
             n = s0 = Double.NaN;
-            ρ0 = null;
-            scale = sqrt(cos2_φ1 - cos2_i) * cos_φ1 / (p2_on_p1 * cos2_φ1 - cos_i);
+            final double cotF = sqrt(cos2_φ1 - cos2_i) / (p2_on_p1*cos2_φ1 - cos_i);  
 // Cotangente of F₁.
+            denormalize.convertBefore(0, cosφ1,      null);
+            denormalize.convertBefore(1, cosφ1*cotF, null);
+            if (!Double.isFinite(cotF) || cotF == 0) {
+                throw invalid(STANDARD_PARALLEL_1);
+            }
         }
-        /*
-         * At this point, all parameters have been processed. Now process to their
-         * validation and the initialization of (de)normalize affine transforms.
-         */
-        final MatrixSIS normalize = context.getMatrix(ContextualParameters.MatrixRole.NORMALIZATION);
-        normalize.convertAfter(0, isConic ? n : cos_φ1, null);
+    }
 
-        final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
-        denormalize.convertBefore(1, scale, ρ0);
-        if (isConic) {
-            denormalize.convertBefore(0, -scale, null);
+    /**
+     * Computes the F₁ or F₂ coefficient using Snyder equation 28-9. Note that this is
the same equation
+     * than F₁′ in above cylindrical case, but with the addition of arc-tangent. This
value is constant
+     * after construction time.
+     *
+     * @param  cos2_φ  square of cosine of φ₁ or φ₂.
+     * @param  source  description of φ argument. Used for error message only.
+     * @return F coefficient for the given φ latitude.
+     */
+    private double F(final double cos2_φ, final ParameterDescriptor<Double> source)
{
+        final double F = atan((p2_on_p1*cos2_φ - cos_i) / sqrt(cos2_φ - cos2_i));
+        if (Double.isFinite(F)) {
+            return F;
         }
+        throw invalid(source);
+    }
+
+    /**
+     * Computes the L₀, L₁ or L₂ coefficient using Snyder equation 28-2a to 28-4a.
+     * This value is constant after construction time.
+     *
+     * @param  sinφ    sine of φ₀, φ₁ or φ₂.
+     * @param  source  description of φ argument. Used for error message only.
+     * @return L coefficient for the given φ latitude.
+     */
+    private double L(final double sinφ, final ParameterDescriptor<Double> source)
{
+        final double λp = -asin(sinφ / sin_i);                      // λ′ in Snyder
equation 28-2a.
+        final double L = atan(tan(λp) * cos_i) - p2_on_p1 * λp;     // Snyder equation
28-3a and 28-4a.
+        if (Double.isFinite(L)) {
+            return L;
+        }
+        throw invalid(source);
+    }
+
+    /**
+     * Returns an exception for a latitude parameter out of range.
+     * The range is assumed given by satellite orbit inclination.
+     *
+     * @param  source  description of invalid φ argument.
+     */
+    private InvalidParameterValueException invalid(final ParameterDescriptor<Double>
source) {
+        final String name  =     source.getName().getCode();
+        final double value =     context.doubleValue(source);
+        final double limit = abs(context.doubleValue(SATELLITE_ORBIT_INCLINATION));
+        return new InvalidParameterValueException(Errors.format(Errors.Keys.ValueOutOfRange_4,
+                                                  name, -limit, limit, new Latitude(value)),
name, value);
     }
 
     /**
-     * Converts the specified (λ,φ) coordinate (units in radians) and stores the result
in {@code dstPts}
-     * (linear distance on a unit sphere). In addition, opportunistically computes the projection
derivative
-     * if {@code derivate} is {@code true}.
+     * Converts the specified (λ,φ) coordinate (units in radians) and stores the result
in {@code dstPts}.
+     * In addition, opportunistically computes the projection derivative if {@code derivate}
is {@code true}.
+     * The results must be multiplied by the denormalization matrix before to get linear
distances.
      *
-     * <cite> The Yaxis lies along the central meridian λ0, y increasing
-     * northerly, and X axis intersects perpendicularly at LATITUDE_OF_ORIGIN
-     * φ0, x increasing easterly.
-     * </cite>
+     * <p>The <var>y</var> axis lies along the central meridian λ₀,
<var>y</var> increasing northerly, and
+     * <var>x</var> axis intersects perpendicularly at latitude of origin φ₀,
<var>x</var> increasing easterly.</p>
      *
      * @return the matrix of the projection derivative at the given source position,
      *         or {@code null} if the {@code derivate} argument is {@code false}.
      * @throws ProjectionException if the coordinates can not be converted.
      */
     @Override
-    public Matrix transform(double[] srcPts, int srcOff,
-                            double[] dstPts, int dstOff,
-                            boolean derivate) throws ProjectionException
+    public Matrix transform(final double[] srcPts, final int srcOff,
+                            final double[] dstPts, final int dstOff,
+                            final boolean derivate) throws ProjectionException
     {
         final double φ         = srcPts[srcOff + 1];
         final double sinφ_sini = sin(φ) / sin_i;
-        final double dλ        = -asin(sinφ_sini);
-        final double tan_dλ    = tan(dλ);
-        final double λt        = atan(tan_dλ * cos_i);
-        double x     = srcPts[srcOff];          // In cylindrical case, x=λ (assuming Earth
radius of 1).
-        double y     = λt - p2_on_p1 * dλ;      // In cylindrical case, y=L. Otherwise
will be adjusted.
-        double dx_dL = 0;       // Part of dx_dφ
-        double dy_dL = 1;       // Part of dy_dφ
+              double λpm       = -asin(sinφ_sini);          // Initialized to λ′ (Snyder
equation 28-2), to be modified later.
+        final double tanλp     = tan(λpm);                  // tan(λ′), saved because
also used in derivative computation.
+        final double λt        = atan(tanλp * cos_i);       // λₜ in Snyder equation
28-3.
+        /*
+         * The (x,y) coordinates below are all is needed for the cylindrical case.
+         * With R and cosφ₁ omitted because multiplied outside this method:
+         *
+         *     x = λ
+         *     y = L    with L = λₜ − (P₂/P₁)λ′     (Snyder equation 28-6)
+         *
+         * Only in conic case, we need to continue with some additional computation.
+         * Note that λpm is NaN if latitude φ is closer to pole than tracking limit.
+         */
+        double x = srcPts[srcOff];
+        double y = λt - p2_on_p1 * λpm;
         if (isConic) {
-            double nλs0 = n*y + s0;
-            if (nλs0 * n < 0) {
+            λpm = n*y + s0;                     // Use this variable for a new purpose.
Needed for derivative.
+            if ((Double.doubleToRawLongBits(λpm) ^ Double.doubleToRawLongBits(n)) < 0)
{
                 /*
-                 * if nλs0 does not have the sign than n, the (x,y) values computed below
would suddenly
+                 * if λpm does not have the sign than n, the (x,y) values computed below
would suddenly
                  * change their sign. The y values lower or greater (depending of n sign)
than -s0/n can
                  * not be plotted. Snyder suggests to use another projection if cosmetic
output is wanted.
                  * For now, we just set the results to NaN (meaning "no result", which is
not the same than
                  * TransformException which means that a result exists but can not be computed).
                  */
-                nλs0 = Double.NaN;
+                λpm = Double.NaN;
             }
-            final double iρ = sin(nλs0);        // Inverse of ρ.
-            y = cos(x) / iρ;
+            final double iρ = sin(λpm);         // Inverse of ρ and without the ρf =
cos(φ₁)⋅sin(F₁)/n part.
+            y = cos(x) / iρ;                    // x already multiplied by n (was done by
normalization matrix).
             x = sin(x) / iρ;                    // Must be last.
-            if (derivate) {
-                final double dρ_dφ = -n / tan(nλs0);
-                dx_dL = x * dρ_dφ;
-                dy_dL = y * dρ_dφ;
-            }
         }
         if (dstPts != null) {
             dstPts[dstOff    ] = x;
@@ -260,12 +320,21 @@ public class SatelliteTracking extends NormalizedProjection {
         if (!derivate) {
             return null;
         }
-        // first computed term associated with derivative of (-dλ/dφ)
-        final double dL_dφ = ((cos(φ) / sin_i) / sqrt(1 - sinφ_sini*sinφ_sini))
-                           * (p2_on_p1 - ((1 + tan_dλ*tan_dλ)*cos_i/(1 + λt*λt)));
-
-        return new Matrix2(isConic ? +y : 1,    dx_dL * dL_dφ,
-                           isConic ? -x : 0,    dy_dL * dL_dφ);
+        /*
+         * Create a derivative matrix initializez with ∂L/∂φ. For cylindrical case where
y = L (in this method),
+         * this is all we need to do. For the conic case, we need to multiply the last column
∂x/∂L and ∂y/∂L.
+         */
+        final Matrix2 d = new Matrix2();
+        d.m11 = ((cos(φ) / sin_i) / sqrt(1 - sinφ_sini*sinφ_sini))
+              * (p2_on_p1 - ((1 + tanλp*tanλp)*cos_i/(1 + λt*λt)));
+        if (isConic) {
+            final double dρ_dφ = -n / tan(λpm);
+            d.m00  = +y;                            // ∂x/∂λ
+            d.m10  = -x;                            // ∂y/∂λ
+            d.m01  =  x * dρ_dφ * d.m11;            // ∂x/∂φ  =  ∂x/∂L ⋅ ∂L/∂φ
+            d.m11 *=  y * dρ_dφ;                    // ∂y/∂φ  =  ∂y/∂L ⋅ ∂L/∂φ
+        }
+        return d;
     }
 
     /**
@@ -275,12 +344,11 @@ public class SatelliteTracking extends NormalizedProjection {
      * @throws ProjectionException if the coordinates can not be converted.
      */
     @Override
-    protected void inverseTransform(double[] srcPts, int srcOff,
-                                    double[] dstPts, int dstOff)
-                                    throws ProjectionException {
-
+    protected void inverseTransform(final double[] srcPts, final int srcOff,
+                                    final double[] dstPts, final int dstOff) throws ProjectionException
+    {
         double x = srcPts[srcOff];
-        double L = srcPts[srcOff + 1];
+        double L = srcPts[srcOff + 1];                  // Multiplication e.g. by F₁′/(R⋅cosφ₁)
already done.
         if (isConic) {
             final double ρ = copySign(hypot(x,L), n);
             x = atan(x / L);                            // Undefined if x = y = 0.
@@ -290,46 +358,20 @@ public class SatelliteTracking extends NormalizedProjection {
          * Approximation of the latitude associated with L coefficient by applying Newton-Raphson
method.
          * Equations 28-24, 28-25 and then 28-22 in Snyder's book.
          */
-        double Δdλ, dλ = -PI/2;
+        final double ic  = 1 / cos2_i;
+        final double pc = p2_on_p1 * cos_i;
+        double λp = -PI/2, Δλp;
         int iter = Formulas.MAXIMUM_ITERATIONS;
         do {
             if (--iter < 0) {
                 throw new ProjectionException(Resources.format(Resources.Keys.NoConvergence));
             }
-            final double dλn = dλ;
-
-            // Alternative calculation with Snyder's eq.  28-20 and 28-21
-//            λt = l + p2_on_p1 * dλ_n;
-//            dλ = atan(tan(λt) / cos_i);
-
-            final double A   = tan(L + p2_on_p1 * dλn) / cos_i;
+            final double A   = tan(L + p2_on_p1 * λp) / cos_i;          // Snyder equation
28-24.
             final double A2  = A*A;
-            Δdλ = (atan(A) - dλn) / (1 - (A2 + 1/cos2_i) * (p2_on_p1*cos_i/(A2 + 1)));
-            dλ  = dλn + Δdλ;
-        } while (abs(Δdλ) >= ANGULAR_TOLERANCE);
+            Δλp = (atan(A) - λp) / (1 - pc*((A2 + ic) / (A2 + 1)));     // Snyder equation
28-25.
+            λp += Δλp;
+        } while (abs(Δλp) >= ANGULAR_TOLERANCE);
         dstPts[dstOff  ] = x;
-        dstPts[dstOff+1] = -asin(sin(dλ) * sin_i);
-    }
-
-    /**
-     * Method to compute the φn coefficient according to equation 28-9
-     * in Snyder's Map Projections manual.
-     *
-     * @param cos2_φn : square of the φn 's cosine.
-     * @return Fn  coefficient associated with the φn latitude.
-     */
-    private double computeFn(final double cos2_φn) {
-        return atan((p2_on_p1 * cos2_φn - cos_i) / sqrt(cos2_φn - cos2_i)); // eq.28-9
in Snyder
-    }
-
-    /**
-     * Method to compute the φn coefficient according to equation 28-3a
-     * in Snyder's Map Projections manual.
-     *
-     * @param dλn  coefficient associated with the φn latitude.
-     * @return λtn  coefficient associated with the φn latitude.
-     */
-    private double computeλtn(final double dλn) {
-        return atan(tan(dλn) * cos_i); // eq.28-3a in Snyder
+        dstPts[dstOff+1] = -asin(sin(λp) * sin_i);                      // Snyder equation
28-22.
     }
 }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SatelliteTrackingTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SatelliteTrackingTest.java
index 38d550b..75aa143 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SatelliteTrackingTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SatelliteTrackingTest.java
@@ -20,7 +20,6 @@ import java.util.Collections;
 import org.opengis.util.FactoryException;
 import org.opengis.parameter.ParameterValueGroup;
 import org.opengis.referencing.operation.TransformException;
-import org.apache.sis.internal.referencing.Formulas;
 import org.apache.sis.internal.referencing.NilReferencingObject;
 import org.apache.sis.internal.referencing.provider.SatelliteTracking;
 import org.apache.sis.measure.Units;
@@ -68,7 +67,7 @@ public final strictfp class SatelliteTrackingTest extends MapProjectionTestCase
     {
         final SatelliteTracking provider = new SatelliteTracking();
         final ParameterValueGroup values = provider.getParameters().createValue();
-        final DefaultEllipsoid sphere = DefaultEllipsoid.createEllipsoid(
+        final DefaultEllipsoid    sphere = DefaultEllipsoid.createEllipsoid(
                 Collections.singletonMap(DefaultEllipsoid.NAME_KEY, NilReferencingObject.UNNAMED),
1, 1, Units.METRE);
 
         values.parameter("semi_major")                 .setValue(sphere.getSemiMajorAxis());
@@ -82,13 +81,28 @@ public final strictfp class SatelliteTrackingTest extends MapProjectionTestCase
         values.parameter("ascending_node_period")      .setValue(1440.0,   Units.MINUTE);
         transform = new MathTransformFactoryMock(provider).createParameterizedTransform(values);
         validate();
+        /*
+         * Assuming that tolerance has been set to the number of fraction digits published
in Snyder tables,
+         * relax the tolerance during inverse transforms for taking in account the increase
in magnitude of
+         * coordinate values. The transform results are between 0 and 1, while the inverse
transform results
+         * are between -90° and 90°, which is an increase in magnitude close to ×100.
+         */
+        toleranceModifier = (tolerance, coordinate, mode) -> {
+            switch (mode) {
+                case INVERSE_TRANSFORM: {
+                    for (int i=0; i<tolerance.length; i++) {
+                        tolerance[i] *= 50;
+                    }
+                    break;
+                }
+            }
+        };
     }
 
     /**
-     * Tests the projection of a few points on a sphere.
-     *
-     * Test based on the numerical example given by Snyder p. 360 to 363 of
-     * <cite> Map Projections - A working Manual</cite>
+     * Tests the projection of a few points using spherical formulas.
+     * Test based on the numerical example given by Snyder pages 360 to 363
+     * of <cite>Map Projections - A working Manual</cite>.
      *
      * @throws FactoryException if an error occurred while creating the map projection.
      * @throws TransformException if an error occurred while projecting a point.
@@ -96,7 +110,6 @@ public final strictfp class SatelliteTrackingTest extends MapProjectionTestCase
     @Test
     public void testCylindricalTransform() throws FactoryException, TransformException {
         tolerance = 1E-7;   // Number of digits in the output values provided by Snyder.
-        tolerance = 1E-5;   // TODO
         createForLandsat(-90, 0, 30, -30);
         assertTrue(isInverseTransformSupported);
         verifyTransform(
@@ -109,9 +122,9 @@ public final strictfp class SatelliteTrackingTest extends MapProjectionTestCase
     }
 
     /**
-     * Tests the projection of a few points on a sphere.
-     * Test based on the numerical example given by Snyder pages 360 to 363 of
-     * <cite>Map Projections - A working Manual</cite>.
+     * Tests the projection of a few points using conic formulas.
+     * Test based on the numerical example given by Snyder pages 360 to 363
+     * of <cite>Map Projections - A working Manual</cite>.
      *
      * @throws FactoryException if an error occurred while creating the map projection.
      * @throws TransformException if an error occurred while projecting a point.
@@ -119,7 +132,6 @@ public final strictfp class SatelliteTrackingTest extends MapProjectionTestCase
     @Test
     public void testConicTransform() throws FactoryException, TransformException {
         tolerance = 1E-7;   // Number of digits in the output values provided by Snyder.
-        tolerance = 1E-5;   // TODO
         createForLandsat(-90, 30, 45, 70);
         assertTrue(isInverseTransformSupported);
         verifyTransform(
@@ -140,9 +152,9 @@ public final strictfp class SatelliteTrackingTest extends MapProjectionTestCase
 
     /**
      * Compares the projection of a few points against expected values computed from intermediate
values
-     * published by Snyder. As Snyder tables in chapter 28 introduce only the values for
<var>n</var>,
-     * <var>F₁</var> and <var>ρ</var> coefficients, the test was
realized by checking these coefficients
-     * in debugger and by extracting the computed results of the projection.
+     * published by Snyder. Snyder tables in chapter 28 do not give directly the (x,y) values.
Instead
+     * the tables give some intermediate values like <var>F₁</var>, which are
verified in debugger.
+     * This method converts intermediate values to final coordinate values to compare.
      *
      * @param  xScale       scale to apply on <var>x</var> values.
      * @param  coordinates  the points to transform.
@@ -160,9 +172,9 @@ public final strictfp class SatelliteTrackingTest extends MapProjectionTestCase
 
     /**
      * Compares the projection of a few points against expected values computed from intermediate
values
-     * published by Snyder. As Snyder tables in chapter 28 introduce only the values for
<var>n</var>,
-     * <var>F₁</var> and <var>ρ</var> coefficients, the test was
realized by checking these coefficients
-     * in debugger and by extracting the computed results of the projection.
+     * published by Snyder. Snyder tables in chapter 28 do not give directly the (x,y) values.
Instead
+     * the tables give some intermediate values like <var>F₁</var> and <var>n</var>,
which are verified
+     * in debugger. This method converts intermediate values to final coordinate values to
compare.
      *
      * @param  λ0           central meridian.
      * @param  n            cone factor <var>n</var>.
@@ -180,7 +192,7 @@ public final strictfp class SatelliteTrackingTest extends MapProjectionTestCase
     }
 
     /**
-     * Tests the projection of a few points on a sphere.
+     * Tests the projection of a few points using cylindrical formulas.
      * Test based on the sample coordinates for several of the Satellite-Tracking projections
      * shown in table 38 of <cite>Map Projections - A working Manual</cite>.
      *
@@ -189,41 +201,52 @@ public final strictfp class SatelliteTrackingTest extends MapProjectionTestCase
      */
     @Test
     public void testCylindricalInternal() throws FactoryException, TransformException {
-        tolerance = NormalizedProjection.ANGULAR_TOLERANCE;
-        tolerance = Formulas.LINEAR_TOLERANCE;                  // TODO
         /*
-         * φ₁ = 0°
+         * First group of 3 columns in Snyder table 38, for φ₁ = 0°.
+         * Snyder gives the following values, which can be verified in a debugger:
+         *
+         *     F₁  =  13.09724°         can be verified with toDegrees(atan(1/cotF))
+         *     x   =  0.017453⋅λ°       can be verified with toRadians(cosφ1)
+         *
+         * Accuracy is set to the number of fraction digits published by Snyder (5 digits)
+         * with tolerance relaxed on the last digit. The verifyCylindricInternal(…) first
+         * argument is the factor in above x equation, with additional digits obtained by
+         * inspecting the value in a debugging session.
          */
+        tolerance = 4E-5;
         createForLandsat(0, 0, 0, 0);
-        verifyCylindricInternal(0.017453,
-                new double[] {                  // (λ,φ) coordinates in degrees to project.
+        verifyCylindricInternal(0.017453292519943295,   // See x in above comment.
+                new double[] {                          // (λ,φ) coordinates in degrees
to project.
                        0,   0,
                       10,   0,
                      -10,  10,
                       60,  40,
                       80,  70,
-                    -120,  80.908               // Tracking limit.
+                    -120,  80.908                       // Tracking limit.
                 },
-                new double[] {                  // Expected (x,y) results on a unit sphere.
+                new double[] {                          // Expected (x,y) results on a unit
sphere.
                        0,   0,
                       10,   0,
                      -10,   0.17579,
                       60,   0.79741,
                       80,   2.34465,
-                    -120,   7.23571             // Projection of tracking limit.
+                    -120,   7.23571                     // Projection of tracking limit.
                 });
         /*
-         * φ₁ = -30°
+         * Second group of 3 columns for φ₁ = -30°.
+         *
+         *     F₁  =  13.96868°
+         *     x   =  0.015115⋅λ°
          */
         createForLandsat(0, 0, -30, 30);
-        verifyCylindricInternal(0.015115,
+        verifyCylindricInternal(0.015114994701951816,   // See x in above comment.
                 new double[] {
                        0,   0,
                       10,   0,
                      -10,  10,
                       60,  40,
                       80,  70,
-                    -120,  80.908
+                    -120,  80.908                       // Tracking limit.
                 },
                 new double[] {
                        0,  0,
@@ -234,17 +257,20 @@ public final strictfp class SatelliteTrackingTest extends MapProjectionTestCase
                     -120,  5.86095
                 });
         /*
-         * φ₁ = 45°
+         * Third group of 3 columns for φ₁ = 45°
+         *
+         *     F₁  =  15.71115°
+         *     x   =  0.012341⋅λ°
          */
         createForLandsat(0, 0, 45, -45);
-        verifyCylindricInternal(0.012341,
+        verifyCylindricInternal(0.012341341494884351,   // See x in above comment.
                 new double[] {
                        0,   0,
                       10,   0,
                      -10,  10,
                       60,  40,
                       80,  70,
-                    -120,  80.908
+                    -120,  80.908                       // Tracking limit.
                 },
                 new double[] {
                        0, 0,
@@ -257,7 +283,7 @@ public final strictfp class SatelliteTrackingTest extends MapProjectionTestCase
     }
 
     /**
-     * Tests the projection of a few points on a sphere.
+     * Tests the projection of a few points using conic formulas.
      * Test based on the sample coordinates for several of the Satellite-Tracking projections
      * shown in table 39 of <cite>Map Projections - A working Manual</cite>,
page 238.
      *
@@ -266,57 +292,71 @@ public final strictfp class SatelliteTrackingTest extends MapProjectionTestCase
      */
     @Test
     public void testConicInternal() throws FactoryException, TransformException {
-        tolerance = NormalizedProjection.ANGULAR_TOLERANCE;
-        tolerance = Formulas.LINEAR_TOLERANCE;                  // TODO
         /*
-         * φ₁ = 30° ; φ₂ = 60°
+         * First group of 3 columns in Snyder table 38, for φ₁ = 30° and φ₂ = 60°.
+         * Snyder gives the following values, which can be verified in a debugger:
+         *
+         *     F₁  =  13.96868°         can be verified with toDegrees(F1)
+         *     n   =   0.49073
+         *
+         * Accuracy is set to the number of fraction digits published by Snyder (5 digits)
+         * with tolerance relaxed on the last digit. The verifyCylindricInternal(…) first
+         * argument is the factor in above x equation, with additional digits obtained by
+         * inspecting the value in a debugging session.
          */
+        tolerance = 3E-5;
         createForLandsat(-90, 0, 30, 60);
-        verifyConicInternal(-90, 0.49073,
-                new double[] {                  // (λ,φ) coordinates in degrees to project.
+        verifyConicInternal(-90, 0.4907267554554259,    // See n in above comment.
+                new double[] {                          // (λ,φ) coordinates in degrees
to project.
                        0, -10,
                        0,   0,
                        0,  10,
                        0,  70,
-                    -120,  80.908               // Tracking limit.
+                    -120,  80.908                       // Tracking limit.
                 },
-                new double[] {                  // Expected (x,y) results on a unit sphere.
+                new double[] {                          // Expected (x,y) results on a unit
sphere.
                     2.67991,  0.46093,
                     2.38332,  0.67369,
                     2.14662,  0.84348,
                     0.98470,  1.67697,
-                    0.50439,  1.89549           // Projection of tracking limit.
+                    0.50439,  1.89549                   // Projection of tracking limit.
                 });
         /*
-         * φ₁ = 45° ; φ₂ = 70°
+         * Second group of 3 columns for φ₁ = 45° and φ₂ = 70°.
+         *
+         *     F₁  =  15.71115°
+         *     n   =  0.69478
          */
         createForLandsat(-90, 0, 45, 70);
-        verifyConicInternal(-90, 0.69478,
+        verifyConicInternal(-90, 0.6947829166657693,    // See n in above comment.
                 new double[] {
                     0, -10,
                     0,   0,
                     0,  10,
                     0,  70,
-                 -120,  80.908                  // Tracking limit.
+                 -120,  80.908                          // Tracking limit.
                 },
                 new double[] {
                     2.92503,  0.90110,
                     2.25035,  1.21232,
                     1.82978,  1.40632,
                     0.57297,  1.98605,
-                    0.28663,  1.982485          // Projection of tracking limit.
+                    0.28663,  1.982485
                 });
         /*
-         * φ₁ = 45° ; φ₂ = 80.908° (the tracking limit).
+         * Second group of 3 columns for φ₁ = 45° and φ₂ = 80.908° (the tracking
limit).
+         *
+         *     F₁  =  15.71115°
+         *     n   =  0.88475
          */
         createForLandsat(-90, 0, 45, 80.908);
-        verifyConicInternal(-90, 0.88475,
+        verifyConicInternal(-90, 0.8847514352390218,    // See n in above comment.
                 new double[] {
                        0, -10,
                        0,   0,
                        0,  10,
                        0,  70,
-                    -120,  80.908
+                    -120,  80.908                       // Tracking limit.
                 },
                 new double[] {
                     4.79153,  1.80001,
@@ -339,8 +379,7 @@ public final strictfp class SatelliteTrackingTest extends MapProjectionTestCase
         createForLandsat(-90, 0, 30, -30);
         final double delta = (1.0 / 60) / 1852;                 // Approximately 1 metre.
         derivativeDeltas = new double[] {delta, delta};
-        tolerance = NormalizedProjection.ANGULAR_TOLERANCE;
-        tolerance = Formulas.LINEAR_TOLERANCE / 100;            // TODO
+        tolerance = 1E-4;
         verifyDerivative( -75, 40);
         verifyDerivative(-100,  3);
         verifyDerivative( -56, 50);
@@ -359,8 +398,7 @@ public final strictfp class SatelliteTrackingTest extends MapProjectionTestCase
         createForLandsat(-90, 30, 45, 70);
         final double delta = (1.0 / 60) / 1852;                 // Approximately 1 metre.
         derivativeDeltas = new double[] {delta, delta};
-        tolerance = NormalizedProjection.ANGULAR_TOLERANCE;
-        tolerance = Formulas.LINEAR_TOLERANCE/100 ;             // TODO
+        tolerance = 1E-4;
         verifyDerivative( -75, 40);
         verifyDerivative(-100,  3);
         verifyDerivative( -56, 50);


Mime
View raw message