in endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/projection/TransverseMercator.java [451:642]
public Matrix transform(final double[] srcPts, final int srcOff,
final double[] dstPts, final int dstOff,
final boolean derivate) throws ProjectionException
{
final double λ = srcPts[srcOff];
final double φ = srcPts[srcOff+1];
if (abs(λ) > 90 * (PI/180)) {
/*
* The Transverse Mercator projection is conceptually a Mercator projection rotated by 90°.
* In Mercator projection, y values tend toward infinity for latitudes close to ±90°.
* In Transverse Mercator, x values tend toward infinity for longitudes close to ±90°
* at equator and after subtraction of central meridian. After we pass the 90° limit,
* the Transverse Mercator formulas produce the same values at (90° + Δ) than at (90° - Δ).
*
* Problem is that 90° is an ordinary longitude value, not even close to the limit of longitude
* values range (±180°). So having f(π/2+Δ, φ) = f(π/2-Δ, φ) results in wrong behavior in some
* algorithms like the one used by `Envelopes.transform(CoordinateOperation, Envelope)`.
* Since a distance of 90° from central meridian is far outside the Transverse Mercator
* domain of validity anyway, we do not let the user go further.
*
* Historical note: in a previous version, we used a limit of 70° instead of 90° because results
* became chaotic after 85°. That limit has been removed in later version because this method now
* behaves like a monotonic function x = f(λ) for fixed φ values. We need to project coordinates
* even in the area where accuracy is bad because projecting those coordinates may happen during
* envelope projections. An envelope may have corners located in invalid projection area even if
* all features inside the envelope have valid coordinates. For "contains" and "intersects" tests
* between envelopes, we do not need accurate coordinates; a monotonic behavior can be sufficient.
*
* Reminder: difference between returning NaN or throwing an exception is as below:
*
* - NaN means "value does not exist or is not a real number".
* - ProjectionException means "value should exist but cannot be computed".
*
* So it is okay to return NaN for values located at Δλ > 90°, but we should throw an exception
* for values at Δλ ≤ 90° if we cannot compute them. Previous version of this method was throwing
* an exception for all Δλ > 70°. Now that we accept all longitudes up to 90°, we return NaN instead.
*/
if (Math.abs(IEEEremainder(λ, 2*PI)) > 90 * (PI/180)) { // More costly check.
return outsideDomainOfValidity(dstPts, dstOff, derivate);
}
}
final double sinλ = sin(λ);
final double ℯsinφ = sin(φ) * eccentricity;
final double Q = asinh(tan(φ)) - atanh(ℯsinφ) * eccentricity;
final double coshQ = cosh(Q); // Cannot be smaller than 1.
final double η0 = atanh(sinλ / coshQ); // Tends toward ±∞ if λ → ±90°.
/*
* Original formula: η0 = atanh(sin(λ) * cos(β)) where
* cos(β) = cos(atan(sinh(Q)))
* = 1 / sqrt(1 + sinh²(Q))
* = 1 / (sqrt(cosh²(Q) - sinh²(Q) + sinh²(Q)))
* = 1 / sqrt(cosh²(Q))
* = 1 / cosh(Q)
*
* So η0 = atanh(sin(λ) / cosh(Q))
*/
final double coshη0 = cosh(η0);
final double ξ0 = asin(tanh(Q) * coshη0);
/*
* Compute sin(2⋅ξ₀), sin(4⋅ξ₀), sin(6⋅ξ₀), sin(8⋅ξ₀) and same for cos, but using the following
* trigonometric identities in order to reduce the number of calls to Math.sin and cos methods.
*/
final double sin_2ξ0 = sin(2*ξ0);
final double cos_2ξ0 = cos(2*ξ0);
final double sin_4ξ0, sin_6ξ0, sin_8ξ0,
cos_4ξ0, cos_6ξ0, cos_8ξ0;
if (!ALLOW_TRIGONOMETRIC_IDENTITIES) {
sin_4ξ0 = sin(4*ξ0);
cos_4ξ0 = cos(4*ξ0);
sin_6ξ0 = sin(6*ξ0);
cos_6ξ0 = cos(6*ξ0);
sin_8ξ0 = sin(8*ξ0);
cos_8ξ0 = cos(8*ξ0);
} else {
final double sin2 = sin_2ξ0 * sin_2ξ0;
final double cos2 = cos_2ξ0 * cos_2ξ0;
sin_4ξ0 = sin_2ξ0 * cos_2ξ0; assert identityEquals(sin_4ξ0, sin(4*ξ0) / 2) : ξ0;
cos_4ξ0 = (cos2 - sin2) * 0.5; assert identityEquals(cos_4ξ0, cos(4*ξ0) / 2) : ξ0;
sin_6ξ0 = (0.75 - sin2) * sin_2ξ0; assert identityEquals(sin_6ξ0, sin(6*ξ0) / 4) : ξ0;
cos_6ξ0 = (cos2 - 0.75) * cos_2ξ0; assert identityEquals(cos_6ξ0, cos(6*ξ0) / 4) : ξ0;
sin_8ξ0 = sin_4ξ0 * cos_4ξ0; assert identityEquals(sin_8ξ0, sin(8*ξ0) / 8) : ξ0;
cos_8ξ0 = 0.125 - sin_4ξ0 * sin_4ξ0; assert identityEquals(cos_8ξ0, cos(8*ξ0) / 8) : ξ0;
}
/*
* Compute sinh(2⋅ξ₀), sinh(4⋅ξ₀), sinh(6⋅ξ₀), sinh(8⋅ξ₀) and same for cosh, but using the following
* hyperbolic identities in order to reduce the number of calls to Math.sinh and cosh methods.
* Note that the formulas are very similar to the above ones, with only some signs reversed.
*/
final double sinh_2η0 = sinh(2*η0);
final double cosh_2η0 = cosh(2*η0);
final double sinh_4η0, sinh_6η0, sinh_8η0,
cosh_4η0, cosh_6η0, cosh_8η0;
if (!ALLOW_TRIGONOMETRIC_IDENTITIES) {
sinh_4η0 = sinh(4*η0);
cosh_4η0 = cosh(4*η0);
sinh_6η0 = sinh(6*η0);
cosh_6η0 = cosh(6*η0);
sinh_8η0 = sinh(8*η0);
cosh_8η0 = cosh(8*η0);
} else {
final double sinh2 = sinh_2η0 * sinh_2η0;
final double cosh2 = cosh_2η0 * cosh_2η0;
cosh_4η0 = (cosh2 + sinh2) * 0.5; assert identityEquals(cosh_4η0, cosh(4*η0) / 2) : η0;
sinh_4η0 = cosh_2η0 * sinh_2η0; assert identityEquals(sinh_4η0, sinh(4*η0) / 2) : η0;
cosh_6η0 = cosh_2η0 * (cosh2 - 0.75); assert identityEquals(cosh_6η0, cosh(6*η0) / 4) : η0;
sinh_6η0 = sinh_2η0 * (sinh2 + 0.75); assert identityEquals(sinh_6η0, sinh(6*η0) / 4) : η0;
cosh_8η0 = sinh_4η0 * sinh_4η0 + 0.125; assert identityEquals(cosh_8η0, cosh(8*η0) / 8) : η0;
sinh_8η0 = sinh_4η0 * cosh_4η0; assert identityEquals(sinh_8η0, sinh(8*η0) / 8) : η0;
}
/*
* The projection of (λ,φ) is given by (η⋅B, ξ⋅B+M₀) — ignoring scale factors and false easting/northing.
* But the B and M₀ parameters have been merged by the constructor with other linear operations in the
* "denormalization" matrix. Consequently, we only need to compute (η,ξ) below.
*/
if (dstPts != null) {
// η(λ,φ)
dstPts[dstOff ] = cf8 * cos_8ξ0 * sinh_8η0
+ cf6 * cos_6ξ0 * sinh_6η0
+ cf4 * cos_4ξ0 * sinh_4η0
+ cf2 * cos_2ξ0 * sinh_2η0
+ η0;
// ξ(λ,φ)
dstPts[dstOff+1] = cf8 * sin_8ξ0 * cosh_8η0
+ cf6 * sin_6ξ0 * cosh_6η0
+ cf4 * sin_4ξ0 * cosh_4η0
+ cf2 * sin_2ξ0 * cosh_2η0
+ ξ0;
}
if (!derivate) {
return null;
}
/*
* Now compute the derivative, if the user asked for it.
*/
final double cosλ = cos(λ); // λ
final double cosφ = cos(φ); // φ
final double cosh2Q = coshQ * coshQ; // Q
final double sinhQ = sinh(Q);
final double tanhQ = tanh(Q);
final double cosh2Q_sin2λ = cosh2Q - (sinλ * sinλ); // Qλ
final double sinhη0 = sinh(η0); // η0
final double sqrt1_thQchη0 = sqrt(1 - (tanhQ * tanhQ) * (coshη0 * coshη0)); // Qη0
// dQ_dλ = 0;
final double dQ_dφ = 1 / cosφ - eccentricitySquared * cosφ / (1 - ℯsinφ * ℯsinφ);
final double dη0_dλ = cosλ * coshQ / cosh2Q_sin2λ;
final double dη0_dφ = -dQ_dφ * sinλ * sinhQ / cosh2Q_sin2λ;
final double dξ0_dλ = sinhQ * sinhη0 * cosλ / (cosh2Q_sin2λ * sqrt1_thQchη0);
final double dξ0_dφ = (dQ_dφ * coshη0 / cosh2Q + dη0_dφ * sinhη0 * tanhQ) / sqrt1_thQchη0;
/*
* Jac(Proj(λ,φ)) is the Jacobian matrix of Proj(λ,φ) function.
* So the derivative of Proj(λ,φ) is defined by:
*
* ┌ ┐
* │ ∂η(λ,φ)/∂λ, ∂η(λ,φ)/∂φ │
* Jac = │ │
* (Proj(λ,φ)) │ ∂ξ(λ,φ)/∂λ, ∂ξ(λ,φ)/∂φ │
* └ ┘
*/
// dξ(λ, φ) / dλ
final double dξ_dλ = dξ0_dλ
+ 2 * (cf2 * (dξ0_dλ * cos_2ξ0 * cosh_2η0 + dη0_dλ * sinh_2η0 * sin_2ξ0)
+ 3 * cf6 * (dξ0_dλ * cos_6ξ0 * cosh_6η0 + dη0_dλ * sinh_6η0 * sin_6ξ0)
+ 2 * (cf4 * (dξ0_dλ * cos_4ξ0 * cosh_4η0 + dη0_dλ * sinh_4η0 * sin_4ξ0)
+ 2 * cf8 * (dξ0_dλ * cos_8ξ0 * cosh_8η0 + dη0_dλ * sinh_8η0 * sin_8ξ0)));
// dξ(λ, φ) / dφ
final double dξ_dφ = dξ0_dφ
+ 2 * (cf2 * (dξ0_dφ * cos_2ξ0 * cosh_2η0 + dη0_dφ * sinh_2η0 * sin_2ξ0)
+ 3 * cf6 * (dξ0_dφ * cos_6ξ0 * cosh_6η0 + dη0_dφ * sinh_6η0 * sin_6ξ0)
+ 2 * (cf4 * (dξ0_dφ * cos_4ξ0 * cosh_4η0 + dη0_dφ * sinh_4η0 * sin_4ξ0)
+ 2 * cf8 * (dξ0_dφ * cos_8ξ0 * cosh_8η0 + dη0_dφ * sinh_8η0 * sin_8ξ0)));
// dη(λ, φ) / dλ
final double dη_dλ = dη0_dλ
+ 2 * (cf2 * (dη0_dλ * cosh_2η0 * cos_2ξ0 - dξ0_dλ * sin_2ξ0 * sinh_2η0)
+ 3 * cf6 * (dη0_dλ * cosh_6η0 * cos_6ξ0 - dξ0_dλ * sin_6ξ0 * sinh_6η0)
+ 2 * (cf4 * (dη0_dλ * cosh_4η0 * cos_4ξ0 - dξ0_dλ * sin_4ξ0 * sinh_4η0)
+ 2 * cf8 * (dη0_dλ * cosh_8η0 * cos_8ξ0 - dξ0_dλ * sin_8ξ0 * sinh_8η0)));
// dη(λ, φ) / dφ
final double dη_dφ = dη0_dφ
+ 2 * (cf2 * (dη0_dφ * cosh_2η0 * cos_2ξ0 - dξ0_dφ * sin_2ξ0 * sinh_2η0)
+ 3 * cf6 * (dη0_dφ * cosh_6η0 * cos_6ξ0 - dξ0_dφ * sin_6ξ0 * sinh_6η0)
+ 2 * (cf4 * (dη0_dφ * cosh_4η0 * cos_4ξ0 - dξ0_dφ * sin_4ξ0 * sinh_4η0)
+ 2 * cf8 * (dη0_dφ * cosh_8η0 * cos_8ξ0 - dξ0_dφ * sin_8ξ0 * sinh_8η0)));
return new Matrix2(dη_dλ, dη_dφ,
dξ_dλ, dξ_dφ);
}