Author: desruisseaux
Date: Fri Jul 31 10:42:33 2015
New Revision: 1693564
URL: http://svn.apache.org/r1693564
Log:
Partial rollback of the use of double-double arithmetic in map projection
initialization.
Our usage of double-double arithmetic has proven its value in matrix
operations, but has
less value in NormalizedProjection subclasses after the point where we use
transcendental
functions (sine, logarithmic, etc.) because we have no double-double versions
of those functions.
By reducing double-double arithmetic usage in those cases, we keep the code
more readable
and avoid to give a false sense of accuracy.
Modified:
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/TimeDependentBWP.java
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
Modified:
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/TimeDependentBWP.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/TimeDependentBWP.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
---
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/TimeDependentBWP.java
[UTF-8] (original)
+++
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/TimeDependentBWP.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -23,6 +23,7 @@ import org.apache.sis.internal.util.Nume
import org.apache.sis.internal.util.DoubleDouble;
import static org.apache.sis.util.ArgumentChecks.*;
+import static org.apache.sis.internal.util.DoubleDouble.verbatim;
import static org.apache.sis.internal.referencing.Formulas.JULIAN_YEAR_LENGTH;
@@ -159,7 +160,7 @@ public class TimeDependentBWP extends Bu
if (time != null) {
final long millis = time.getTime() - timeReference;
if (millis != 0) { // Returns null for 0 as an optimization.
- final DoubleDouble period = new DoubleDouble(millis, 0);
+ final DoubleDouble period = verbatim(millis);
period.divide(1000 * JULIAN_YEAR_LENGTH, 0);
return period;
}
Modified:
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
---
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
[UTF-8] (original)
+++
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -28,14 +28,20 @@ import org.apache.sis.referencing.operat
import static java.lang.Math.*;
import static org.apache.sis.util.ArgumentChecks.ensureNonNull;
+import static org.apache.sis.internal.util.DoubleDouble.verbatim;
/**
* Helper class for map projection constructions, providing formulas normally
needed only at construction time.
- * Since map projection constructions should not happen very often, we afford
using double-double arithmetic here.
+ * Since map projection constructions should not happen very often, we afford
using some double-double arithmetic.
* The main intend is not to provide more accurate coordinate conversions
(while it may be a nice side-effect),
- * but rather to increase the chances that the concatenations of
(de)normalization matrices with the matrices of
- * other transforms give back identity matrices when such result is expected.
+ * but to improve the result of concatenations of (de)normalization matrices
with the matrices of other transforms,
+ * as found in transformation chains.
+ *
+ * <p>As a general rule, we stop storing result with double-double precision
after the point where we need
+ * transcendental functions (sine, logarithm, <i>etc.</i>), since we do not
have double-double versions of
+ * those functions. Digits after the {@code double} part are usually not
significant in such cases, except
+ * in some relatively rare scenarios like 1 ± (a result much smaller than
1).</p>
*
* @author Martin Desruisseaux (Geomatys)
* @since 0.6
@@ -71,7 +77,13 @@ final class Initializer {
final byte variant;
/**
- * Creates a new initializer.
+ * Creates a new initializer. The parameters are described in
+ * {@link NormalizedProjection#NormalizedProjection(OperationMethod,
Parameters, Map)}.
+ *
+ * @param method Description of the map projection parameters.
+ * @param parameters The parameters of the projection to be created.
+ * @param roles Parameters to look for <cite>central meridian</cite>,
<cite>scale factor</cite>,
+ * <cite>false easting</cite>, <cite>false northing</cite> and
other values.
*/
Initializer(final OperationMethod method, final Parameters parameters,
final Map<ParameterRole, ? extends ParameterDescriptor<Double>>
roles,
@@ -102,25 +114,23 @@ final class Initializer {
- getAndStore(roles.get(ParameterRole.FALSE_SOUTHING));
excentricitySquared = new DoubleDouble();
- final DoubleDouble k = new DoubleDouble(a); // The value by which to
multiply all results of normalized projection.
+ DoubleDouble k = new DoubleDouble(a); // The value by which to
multiply all results of normalized projection.
if (a != b) {
/*
- * Equivalent Java code for the following lines:
- *
- * final double rs = b / a;
- * excentricitySquared = 1 - (rs * rs);
+ * ℯ² = 1 - (b/a)²
*
- * Test show that double-double arithmetic here makes a difference
in the 3 last digits for WGS84 ellipsoid.
- * Those 3 digits are not significant since the parameter are not
so accurate (furthermore the 'b' parameter
- * used below may have been computed from the inverse flattening
factor).
+ * Double-double arithmetic here makes a difference in the 3 last
digits for WGS84 ellipsoid.
*/
- final DoubleDouble rs = new DoubleDouble(b);
- final double eb = rs.error;
- rs.divide(k); // rs = b/a
- rs.multiply(rs);
- excentricitySquared.value = 1;
- excentricitySquared.subtract(rs);
-
+ if (DoubleDouble.DISABLED) {
+ final double rs = b / a;
+ excentricitySquared.value = 1 - (rs * rs);
+ } else {
+ final DoubleDouble rs = new DoubleDouble(b);
+ rs.divide(k); // rs = b/a
+ rs.multiply(rs);
+ excentricitySquared.value = 1;
+ excentricitySquared.subtract(rs);
+ }
final ParameterDescriptor<Double> radius =
roles.get(ParameterRole.LATITUDE_OF_CONFORMAL_SPHERE_RADIUS);
if (radius != null) {
/*
@@ -139,13 +149,8 @@ final class Initializer {
* final double sinφ =
sin(toRadians(parameters.doubleValue(radius)));
* k = b / (1 - excentricitySquared * (sinφ*sinφ));
*/
- final DoubleDouble t = new
DoubleDouble(sin(toRadians(parameters.doubleValue(radius))), 0);
- t.multiply(t);
- t.multiply(excentricitySquared);
- k.clear();
- k.value = 1;
- k.subtract(t);
- k.inverseDivide(b, eb);
+ k = rν2(sin(toRadians(parameters.doubleValue(radius))));
+ k.inverseDivide(b, 0);
}
}
context.normalizeGeographicInputs(λ0);
@@ -198,28 +203,19 @@ final class Initializer {
return value;
}
-
-
-
-
//////////////////////////////////////////////////////////////////////////////////////////
- ////////
////////
- //////// FORMULAS FROM EPSG or SNYDER
////////
- ////////
////////
-
//////////////////////////////////////////////////////////////////////////////////////////
-
/**
- * Computes the reciprocal of the radius of curvature of the ellipsoid
perpendicular to the meridian at latitude φ.
- * That radius of curvature is:
+ * Computes the square of the reciprocal of the radius of curvature of the
ellipsoid
+ * perpendicular to the meridian at latitude φ. That radius of curvature
is:
*
* <blockquote>ν = 1 / √(1 - ℯ²⋅sin²φ)</blockquote>
*
- * This method returns 1/ν.
+ * This method returns 1/ν².
*
* <div class="section">Relationship with Snyder</div>
* This is related to functions (14-15) from Snyder (used for computation
of scale factors
* at the true scale latitude) as below:
*
- * <blockquote>m = cosφ / rν</blockquote>
+ * <blockquote>m = cosφ / sqrt(rν²)</blockquote>
*
* Special cases:
* <ul>
@@ -228,24 +224,46 @@ final class Initializer {
* (otherwise we get {@link Double#NaN}).</li>
* </ul>
*
- * @param sinφ The sine of the φ latitude in radians.
- * @return Reciprocal of the radius of curvature of the ellipsoid
perpendicular to the meridian at latitude φ.
- */
- final DoubleDouble rν(final double sinφ) {
- /*
- * Equivalent Java code:
- *
- * return sqrt(1 - excentricitySquared * (sinφ*sinφ));
- */
- final DoubleDouble t = new DoubleDouble(sinφ, 0);
+ * @param sinφ The sine of the φ latitude.
+ * @return Reciprocal squared of the radius of curvature of the ellipsoid
+ * perpendicular to the meridian at latitude φ.
+ */
+ private DoubleDouble rν2(final double sinφ) {
+ if (DoubleDouble.DISABLED) {
+ return verbatim(1 - excentricitySquared.value * (sinφ*sinφ));
+ }
+ final DoubleDouble t = verbatim(sinφ);
t.multiply(t);
t.multiply(excentricitySquared);
+
+ // Save result of ℯ²⋅sin²φ
final double value = t.value;
final double error = t.error;
+
+ // Compute 1 - above. Since above may be small, this
+ // is where double-double arithmetic has more value.
t.clear();
t.value = 1;
t.subtract(value, error);
- t.sqrt();
return t;
}
+
+ /**
+ * Returns the scale factor at latitude φ. This is computed as:
+ *
+ * <blockquote>cosφ / sqrt(rν2(sinφ))</blockquote>
+ *
+ * The result is returned as a {@code double} because the limited
precision of {@code sinφ} and {@code cosφ}
+ * makes the error term meaningless. We use double-double arithmetic only
for intermediate calculation.
+ *
+ * @param sinφ The sine of the φ latitude.
+ * @param cosφ The cosine of the φ latitude.
+ * @return Scale factor at latitude φ.
+ */
+ final double scaleAtφ(final double sinφ, final double cosφ) {
+ final DoubleDouble s = rν2(sinφ);
+ s.sqrt();
+ s.inverseDivide(cosφ, 0);
+ return s.value;
+ }
}
Modified:
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
---
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
[UTF-8] (original)
+++
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -40,6 +40,7 @@ import org.apache.sis.util.Workaround;
import static java.lang.Math.*;
import static java.lang.Double.*;
import static org.apache.sis.math.MathFunctions.isPositive;
+import static org.apache.sis.internal.util.DoubleDouble.verbatim;
/**
@@ -96,8 +97,20 @@ public class LambertConicConformal exten
/**
* Constant for the Belgium 2SP case. This is 29.2985 seconds, given here
in radians.
+ *
+ * <div class="note"><b>Tip:</b> how to verify the value:
+ * {@preformat java
+ * BigDecimal a = new BigDecimal(BELGE_A.value);
+ * a = a.add (new BigDecimal(BELGE_A.error));
+ * a = a.multiply(new
BigDecimal("57.29577951308232087679815481410517"));
+ * a = a.multiply(new BigDecimal(60 * 60));
+ * System.out.println(a);
+ * }
+ * </div>
*/
- private static final double BELGE_A = 0.00014204313635987700;
+ static Number belgeA() {
+ return new DoubleDouble(-1.420431363598774E-4,
-1.1777378450498224E-20);
+ }
/**
* Internal coefficients for computation, depending only on values of
standards parallels.
@@ -243,62 +256,35 @@ public class LambertConicConformal exten
φ1 = toRadians(φ1);
φ2 = toRadians(φ2);
/*
- * Computes constants. We do not need to use special formulas for the
spherical case below,
+ * Compute constants. We do not need to use special formulas for the
spherical case below,
* since rν(sinφ) = 1 and expOfNorthing(φ) = tan(π/4 + φ/2) when the
excentricity is zero.
* However we need special formulas for φ1 ≈ φ2 in the calculation of
n, otherwise we got
* a 0/0 indetermination.
- *
- * Opportunistically use double-double arithmetic below since this is
what we will store in the
- * (de)normalization matrices. The extra precision that we get is not
necessarily significant,
- * but we do that more in an attempt to reduce rounding errors in
concatenations of a sequence
- * of MathTransforms (through matrix multiplications) than for map
projection precisions.
- * Equivalent Java code for the following double-double arithmetic:
- *
- * final double m1 = cos(φ1) / rν(sinφ1);
*/
final double sinφ1 = sin(φ1);
- final DoubleDouble m1 = initializer.rν(sinφ1);
- m1.inverseDivide(cos(φ1), 0);
+ final double m1 = initializer.scaleAtφ(sinφ1, cos(φ1));
final double t1 = expOfNorthing(φ1, excentricity*sinφ1);
/*
- * Computes n = (ln m₁ – ln m₂) / (ln t₁ – ln t₂), which we rewrite as
ln(m₁/m₂) / ln(t₁/t₂)
+ * Compute n = (ln m₁ – ln m₂) / (ln t₁ – ln t₂), which we rewrite as
ln(m₁/m₂) / ln(t₁/t₂)
* since division is less at risk of precision lost than subtraction.
Note that this equation
* tends toward 0/0 if φ₁ ≈ φ₂, which force us to do a special check
for the SP1 case.
- *
- * Equivalent Java code for the following double-double arithmetic:
- *
- * final double sinφ2 = sin(φ2);
- * final double m2 = cos(φ2) / rν(sinφ2);
- * final double t2 = expOfNorthing(φ2, excentricity*sinφ2);
- * n = log(m1/m2) / log(t1/t2); // Tend toward 0/0 if φ1 ≈ φ2.
*/
- final DoubleDouble F = new DoubleDouble();
if (abs(φ1 - φ2) >= ANGULAR_TOLERANCE) { // Should be 'true' for 2SP
case.
final double sinφ2 = sin(φ2);
- final DoubleDouble m2 = initializer.rν(sinφ2);
- m2.inverseDivide(cos(φ2), 0);
+ final double m2 = initializer.scaleAtφ(sinφ2, cos(φ2));
final double t2 = expOfNorthing(φ2, excentricity*sinφ2);
- m2.inverseDivide(m1);
- F.value = log(m2.value);
- F.divide(log(t1/t2), 0);
+ n = log(m1/m2) / log(t1/t2); // Tend toward 0/0 if φ1 ≈ φ2.
} else {
- F.value = -sinφ1;
+ n = -sinφ1;
}
- n = F.value;
/*
- * Scale factor for longitudes, stored now before we modify the F
value.
- */
- final DoubleDouble sx = new DoubleDouble(F);
- if (!isNorth) {
- sx.negate();
- }
- /*
- * Computes F = m₁/(n⋅t₁ⁿ) from Geomatics Guidance Note number 7.
+ * Compute F = m₁/(n⋅t₁ⁿ) from Geomatics Guidance Note number 7.
* Following constants will be stored in the denormalization matrix,
to be applied
* after the non-linear formulas implemented by this
LambertConicConformal class.
* Opportunistically use double-double arithmetic since the matrix
coefficients will
* be stored in that format anyway. This makes a change in the 2 or 3
last digits.
*/
+ final DoubleDouble F = verbatim(n);
F.multiply(pow(t1, n), 0);
F.inverseDivide(m1);
if (!isNorth) {
@@ -306,16 +292,16 @@ public class LambertConicConformal exten
}
/*
* Compute the radius of the parallel of latitude of the false origin.
- * This is related to the "ρ0" term in Snyder. From EPG guide:
+ * This is related to the "ρ₀" term in Snyder. From EPG guide:
*
* r = a⋅F⋅tⁿ where (in our case) a=1 and t is our
'expOfNorthing' function.
*
* EPSG uses this term in the computation of y = FN + rF – r⋅cos(θ).
*/
- final DoubleDouble rF = new DoubleDouble(); // Initialized to zero.
+ DoubleDouble rF = null;
if (φ0 != copySign(PI/2, -n)) { // For reducing the rounding error
documented in expOfNorthing(+π/2).
- rF.value = pow(expOfNorthing(φ0, excentricity*sin(φ0)), n);
- rF.multiply(F);
+ rF = new DoubleDouble(F);
+ rF.multiply(pow(expOfNorthing(φ0, excentricity*sin(φ0)), n), 0);
}
/*
* At this point, all parameters have been processed. Now store
@@ -334,14 +320,19 @@ public class LambertConicConformal exten
* - Multiply by the scale factor (done by the super-class
constructor).
* - Add false easting and false northing (done by the super-class
constructor).
*/
- final MatrixSIS normalize = context.getMatrix(true);
- normalize.convertAfter(0, sx, (initializer.variant == BELGIUM) ? new
DoubleDouble(-BELGE_A, 0) : null);
+ DoubleDouble sλ = verbatim(n);
+ DoubleDouble sφ = null;
if (isNorth) {
- normalize.convertAfter(1, new DoubleDouble(-1, 0), null);
+ // Reverse the sign of either longitude or latitude values before
map projection.
+ sφ = verbatim(-1);
+ } else {
+ sλ.negate();
}
+ final MatrixSIS normalize = context.getMatrix(true);
final MatrixSIS denormalize = context.getMatrix(false);
- denormalize.convertBefore(0, F, null);
- F.negate();
+ normalize .convertAfter(0, sλ, (initializer.variant == BELGIUM) ?
belgeA() : null);
+ normalize .convertAfter(1, sφ, null);
+ denormalize.convertBefore(0, F, null); F.negate();
denormalize.convertBefore(1, F, rF);
}
Modified:
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
---
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
[UTF-8] (original)
+++
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -38,6 +38,7 @@ import org.apache.sis.util.Workaround;
import static java.lang.Math.*;
import static java.lang.Double.*;
import static org.apache.sis.math.MathFunctions.isPositive;
+import static org.apache.sis.internal.util.DoubleDouble.verbatim;
/**
@@ -228,8 +229,7 @@ public class Mercator extends ConformalP
* if they really want, since we sometime see such CRS definitions.
*/
final double φ1 =
toRadians(initializer.getAndStore(Mercator2SP.STANDARD_PARALLEL));
- final DoubleDouble k0 = new DoubleDouble(cos(φ1), 0);
- k0.divide(initializer.rν(sin(φ1)));
+ final Number k0 = verbatim(initializer.scaleAtφ(sin(φ1), cos(φ1)));
/*
* In principle we should rotate the central meridian (λ0) in the
normalization transform, as below:
*
@@ -251,11 +251,11 @@ public class Mercator extends ConformalP
denormalize.convertBefore(0, null, offset);
}
if (φ0 != 0) {
- denormalize.convertBefore(1, null, new
DoubleDouble(-log(expOfNorthing(φ0, excentricity * sin(φ0)))));
+ denormalize.convertBefore(1, null, verbatim(-log(expOfNorthing(φ0,
excentricity * sin(φ0)))));
}
if (variant == MILLER) {
- normalize.convertBefore(1, new DoubleDouble(0.80), null);
- denormalize.convertBefore(1, new DoubleDouble(1.25), null);
+ normalize .convertBefore(1, 0.80, null);
+ denormalize.convertBefore(1, 1.25, null);
}
/*
* At this point we are done, but we add here a little bit a maniac
precision hunting.
@@ -281,9 +281,9 @@ public class Mercator extends ConformalP
* those remaning lines of code.
*/
if (φ0 == 0 && isPositive(φ1 != 0 ? φ1 : φ0)) {
- final DoubleDouble revert = new DoubleDouble(-1, 0);
- normalize.convertBefore(1, revert, null);
- denormalize.convertBefore(1, revert, null); // Must be before
false easting/northing.
+ final Number reverseSign = verbatim(-1);
+ normalize .convertBefore(1, reverseSign, null);
+ denormalize.convertBefore(1, reverseSign, null); // Must be
before false easting/northing.
}
}
Modified:
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
---
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
[UTF-8] (original)
+++
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -29,7 +29,6 @@ import org.apache.sis.internal.referenci
import org.apache.sis.internal.referencing.provider.PolarStereographicB;
import org.apache.sis.internal.referencing.provider.PolarStereographicC;
import org.apache.sis.internal.referencing.Formulas;
-import org.apache.sis.internal.util.DoubleDouble;
import org.apache.sis.parameter.Parameters;
import org.apache.sis.util.resources.Errors;
import org.apache.sis.util.Workaround;
@@ -37,6 +36,7 @@ import org.apache.sis.measure.Latitude;
import org.apache.sis.math.MathFunctions;
import static java.lang.Math.*;
+import static org.apache.sis.internal.util.DoubleDouble.verbatim;
/**
@@ -185,8 +185,7 @@ public class PolarStereographic extends
φ1 = -φ1;
}
φ1 = toRadians(φ1); // May be anything in [-π/2 … 0] range.
- final DoubleDouble ρ;
- DoubleDouble ρF = null; // Actually -ρF (compared to EPSG guide).
+ final Number ρ, ρF; // This ρF is actually -ρF in EPSG guide.
if (abs(φ1 + PI/2) < ANGULAR_TOLERANCE) {
/*
* Polar Stereographic (variant A)
@@ -199,19 +198,9 @@ public class PolarStereographic extends
* - the 't' factor, because it needs to be computed in the
transform(…) method.
*
* In the spherical case, should give ρ == 2.
- *
- * Opportunistically use double-double arithmetic below since this
is what we will store in the
- * (de)normalization matrices. The extra precision that we get is
not necessarily significant,
- * but we do that more in an attempt to reduce rounding errors in
concatenations of a sequence
- * of MathTransforms (through matrix multiplications) than for map
projection precisions.
- * Equivalent Java code for the following double-double arithmetic:
- *
- * ρ = 2 / sqrt(pow(1+excentricity, 1+excentricity) *
pow(1-excentricity, 1-excentricity));
*/
- ρ = new DoubleDouble(pow(1+excentricity, 1+excentricity), 0);
- ρ.multiply (pow(1-excentricity, 1-excentricity), 0);
- ρ.sqrt();
- ρ.inverseDivide(2, 0);
+ ρ = verbatim(2 / sqrt(pow(1+excentricity, 1+excentricity) *
pow(1-excentricity, 1-excentricity)));
+ ρF = null;
} else {
/*
* Polar Stereographic (variant B or C)
@@ -229,21 +218,11 @@ public class PolarStereographic extends
* k₀ = ρ⋅√[…]/2 but we do not need that value.
*
* In the spherical case, should give ρ = 1 + sinφ1 (Synder 21-7
and 21-11).
- *
- * Equivalent Java code for the following double-double arithmetic:
- *
- * final double mF = cos(φ1) / rν(sinφ1);
- * ρ = mF / expOfNorthing(φ1, excentricity*sinφ1);
- * if (variant == C) ρF = -mF;
*/
final double sinφ1 = sin(φ1);
- ρ = initializer.rν(sinφ1);
- ρ.inverseDivide(cos(φ1), 0);
- if (variant == C) {
- ρF = new DoubleDouble(ρ);
- ρF.negate();
- }
- ρ.divide(expOfNorthing(φ1, excentricity*sinφ1), 0);
+ final double mF = initializer.scaleAtφ(sinφ1, cos(φ1));
+ ρ = verbatim(mF / expOfNorthing(φ1, excentricity*sinφ1));
+ ρF = (variant == C) ? verbatim(-mF) : null;
}
/*
* At this point, all parameters have been processed. Now process to
their
@@ -253,8 +232,10 @@ public class PolarStereographic extends
denormalize.convertBefore(0, ρ, null);
denormalize.convertBefore(1, ρ, ρF);
if (isNorth) {
- context.getMatrix(true).convertAfter(1, -1, null);
- denormalize.convertBefore(1, -1, null);
+ final Number reverseSign = verbatim(-1);
+ final MatrixSIS normalize = context.getMatrix(true);
+ normalize .convertAfter (1, reverseSign, null);
+ denormalize.convertBefore(1, reverseSign, null);
}
}
Modified:
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
---
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
[UTF-8] (original)
+++
sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -31,6 +31,7 @@ import org.apache.sis.util.Workaround;
import static java.lang.Math.*;
import static org.apache.sis.math.MathFunctions.asinh;
import static org.apache.sis.math.MathFunctions.atanh;
+import static org.apache.sis.internal.util.DoubleDouble.verbatim;
/**
@@ -129,7 +130,7 @@ public class TransverseMercator extends
* Compute B = (n4/64 + n2/4 + 1) / (n + 1)
* Opportunistically uses double-double arithmetic since we use it
anyway for denormalization matrix.
*/
- final DoubleDouble B = new DoubleDouble(n);
+ final DoubleDouble B = verbatim(n);
B.add(1);
B.inverseDivide(1, n4/64 + n2/4);
/*
@@ -160,7 +161,7 @@ public class TransverseMercator extends
*/
final double Q = asinh(tan(φ0)) - excentricity * atanh(excentricity *
sin(φ0));
final double β = atan(sinh(Q));
- final DoubleDouble M0 = new DoubleDouble(β, 0);
+ final DoubleDouble M0 = verbatim(β);
M0.add(h1 * sin(2*β), 0);
M0.add(h2 * sin(4*β), 0);
M0.add(h3 * sin(6*β), 0);
Modified:
sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
---
sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
[UTF-8] (original)
+++
sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -17,6 +17,7 @@
package org.apache.sis.referencing.operation.projection;
import java.util.Random;
+import java.math.BigDecimal;
import org.opengis.util.FactoryException;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
@@ -28,6 +29,7 @@ import org.apache.sis.internal.referenci
import org.apache.sis.internal.referencing.provider.LambertConformalMichigan;
import org.apache.sis.referencing.operation.transform.CoordinateDomain;
import org.apache.sis.parameter.Parameters;
+import org.apache.sis.internal.util.DoubleDouble;
import org.apache.sis.test.DependsOnMethod;
import org.apache.sis.test.DependsOn;
import org.apache.sis.test.TestUtilities;
@@ -53,6 +55,22 @@ import static org.junit.Assert.*;
@DependsOn(ConformalProjectionTest.class)
public final strictfp class LambertConicConformalTest extends
MapProjectionTestCase {
/**
+ * Verifies the value of the constant used in <cite>"Lambert Conic
Conformal (2SP Belgium)"</cite> projection.
+ *
+ * @see #testLambertConicConformalBelgium()
+ */
+ @Test
+ public void verifyBelgeConstant() {
+ final DoubleDouble BELGE_A = (DoubleDouble)
LambertConicConformal.belgeA();
+ BigDecimal a = new BigDecimal(BELGE_A.value);
+ a = a.add (new BigDecimal(BELGE_A.error));
+ a = a.multiply(new BigDecimal("57.29577951308232087679815481410517"));
// Conversion from radians to degrees.
+ a = a.multiply(new BigDecimal(60 * 60));
// Conversion from degrees to seconds.
+ a = a.add (new BigDecimal("29.2985"));
// The standard value.
+ assertTrue(Math.abs(a.doubleValue()) < 1E-31);
+ }
+
+ /**
* Creates a new instance of {@link LambertConicConformal}. See the class
javadoc for an explanation
* about why we ask only for the latitude of origin and not the standard
parallels.
*
@@ -199,7 +217,7 @@ public final strictfp class LambertConic
* @see
org.opengis.test.referencing.ParameterizedTransformTest#testLambertConicConformal1SP()
*/
@Test
- @DependsOnMethod("testLambertConicConformal2SP")
+ @DependsOnMethod({"testLambertConicConformal2SP", "verifyBelgeConstant"})
public void testLambertConicConformalBelgium() throws FactoryException,
TransformException {
createGeoApiTest(new
LambertConformalBelgium()).testLambertConicConformalBelgium();
}
Modified:
sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
URL:
http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
---
sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
[UTF-8] (original)
+++
sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -237,6 +237,20 @@ public final class DoubleDouble extends
}
/**
+ * Uses the given value verbatim, without inferring an error term for
double-double arithmetic.
+ * We use this method when the value has been computed using
transcendental functions (cosine,
+ * logarithm, <i>etc.</i>) in which case there is no way we can infer a
meaningful error term.
+ *
+ * <p>We use this method both for readability and for making easier to
search where such thing occur.</p>
+ *
+ * @param value The value to wrap in a {@code DoubleDouble} instance.
+ * @return A {@code DoubleDouble} containing exactly the given value,
without error term.
+ */
+ public static DoubleDouble verbatim(final double value) {
+ return new DoubleDouble(value, 0);
+ }
+
+ /**
* Returns a new {@code DoubleDouble} instance initialized to the
conversion factor
* from radians to angular degrees.
*
