in endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/GeodesicsOnEllipsoid.java [564:851]
final void computeDistance() {
canComputeDistance();
/*
* The algorithm in this method requires the following canonical configuration:
*
* Negative latitude of starting point: φ₁ ≤ 0
* Ending point latitude smaller in magnitude: φ₁ ≤ φ₂ ≤ -φ₁
* Positive longitude difference: 0 ≤ ∆λ ≤ π
* (Consequence of above): 0 ≤ α₀ ≤ π/2
*
* If the given points do not met above conditions, then we need to swap start and end points or to
* swap coordinate signs. We apply those changes on local variables only, not on the class fields.
* We will need to apply the converse of those changes in the final results.
*/
double φ1 = this.φ1;
double φ2 = this.φ2;
double Δλ = IEEEremainder(λ2 - λ1, 2*PI); // In [-π … +π] range.
final boolean swapPoints = abs(φ2) > abs(φ1);
if (swapPoints) {
φ1 = this.φ2;
φ2 = this.φ1;
Δλ = -Δλ;
}
final boolean inverseLatitudeSigns = φ1 > 0;
if (inverseLatitudeSigns) {
φ1 = -φ1;
φ2 = -φ2;
}
final boolean inverseLongitudeSigns = Δλ < 0;
if (inverseLongitudeSigns) {
Δλ = -Δλ;
}
/*
* Compute an approximation of the azimuth α₁ at starting point. This estimation will be refined by iteration
* in the loop later, but that iteration will not converge if the first α₁ estimation is not good enough. The
* general formula does not give good α₁ initial value for antipodal points, so we need to check for special
* cases first:
*
* 1) Nearly antipodal points with φ₁ = -φ₂.
* 2) Nearly antipodal points with φ₁ slightly different than -φ₂.
* 3) Equatorial case: φ₁ = φ₂ = 0 but restricted to ∆λ ≤ (1-f)⋅π.
* 4) Meridional case: ∆λ = 0 or ∆λ = π (handled by general case in this method).
*/
if (φ1 > -LATITUDE_THRESHOLD) { // Sufficient test because φ₁ ≤ 0 and |φ₂| ≤ φ₁
/*
* Points on equator but not nearly anti-podal. The geodesic is an arc on equator and the azimuths
* are α₁ = α₂ = ±90°. We need this special case because when φ = 0, the general case get sinβ = 0
* then σ = atan2(sinβ, cosα⋅cosβ) = 0, which result in a distance of 0. I have not yet understood
* how to use the general formulas in such case. This code is a workaround in the meantime.
*
* See https://issues.apache.org/jira/browse/SIS-467
*/
if (Δλ > axisRatio * PI) {
// Karney's special case documented before equation 45.
throw new GeodeticException("Cannot compute geodesics for antipodal points on equator.");
}
super.computeDistance();
return;
}
/*
* Reduced latitudes β (Karney 6). Actually we don't need the β angles
* (except for a special case), but rather their sine and cosine values.
*/
final double tanβ1 = axisRatio * tan(φ1);
final double tanβ2 = axisRatio * tan(φ2);
final double cosβ1 = 1 / sqrt(1 + tanβ1*tanβ1);
final double cosβ2 = 1 / sqrt(1 + tanβ2*tanβ2);
final double sinβ1 = tanβ1 * cosβ1;
final double sinβ2 = tanβ2 * cosβ2;
double α1;
if (Δλ >= PI - NEARLY_ANTIPODAL_Δλ && abs(φ1 + φ2) <= NEARLY_ANTIPODAL_Δλ) {
/*
* Nearly antipodal points. Karney's equations 53 are reproduced on the left side
* (with f = 1 - b/a), which we replace by the equations on the right side below:
*
* Δ = f⋅a⋅π⋅cos²β₁ ┃ Δ₁ = f⋅π⋅cosβ₁
* ∆λ = π + Δ/(a⋅cosβ₁)⋅x ┃ ∆λ = π + Δ₁⋅x
* β₂ = -β₁ + (Δ/a)⋅y ┃ β₂ = -β₁ + (Δ₁⋅cosβ₁)⋅y
*/
final double β1 = atan2(sinβ1, cosβ1);
final double β2 = atan2(sinβ2, cosβ2);
final double Δ1 = (1 - axisRatio) * PI * cosβ1; // Differ from Karney by a⋅cosβ₁ factor.
final double y = (β2 + β1) / (Δ1*cosβ1);
final double x = (PI - Δλ) / Δ1; // Opposite sign of Karney. We have x ≥ 0.
final double x2 = x*x;
final double y2 = y*y;
if (STORE_LOCAL_VARIABLES) { // For comparing values with Karney table 4.
store("x", -x); // Sign used by Karney.
store("y", y);
}
if (y2 < 1E-12) { // Empirical threshold. See μ(…) for more information.
α1 = (x2 > 1) ? PI/2 : atan(x / sqrt(1 - x2)); // (Karney 57) with opposite sign of x. Result in α₁ ≥ 0.
} else {
final double μ = μ(x2, y2);
α1 = atan2(x*μ, (y*(1+μ))); // (Karney 56) with opposite sign of x.
if (STORE_LOCAL_VARIABLES) {
store("μ", μ);
}
}
} else {
/*
* Usual case (non-antipodal). Estimation is based on variation of geodetic longitude λ
* and spherical longitude ω on the auxiliary sphere. Karney makes a special case for
* Δλ = 0 and Δλ = π by defining α₁ = Δλ. However, formulas below produce the same result.
*/
double w = (cosβ1 + cosβ2) / 2;
w = sqrt(1 - eccentricitySquared * (w*w));
final double Δω = Δλ / w;
final double cω = cos(Δω);
final double sω = sin(Δω);
final double αx = cosβ1*sinβ2 - sinβ1*cosβ2*cω;
final double αy = cosβ2*sω;
α1 = atan2(αy, αx); // (Karney 49)
if (STORE_LOCAL_VARIABLES) { // For comparing values with Karney table 3.
store("ωb", w);
store("Δω", Δω);
store("α₂", atan2(cosβ1*sω, cosβ1*cω*sinβ2 - sinβ1*cosβ2)); // (Karney 50)
store("Δσ", atan2(hypot(αy, αx), cosβ1*cω*cosβ2 + sinβ1*sinβ2)); // (Karney 51)
// For small distances we could stop here. But it is easier to let iteration does its work.
}
}
/*
* Stores locale variables for comparison against Karney tables 4, 5 and 6. Values β₁ and β₂ are kept
* constant during all this method. Value α₁ is a first estimation and will be updated during iteration.
* Note that because Karney separate calculation of α₁ and remaining calculation in two separated tables,
* we need to truncate α₁ to the same number of digits than Karney in order to get the same numbers in
* the rest of this method.
*/
if (STORE_LOCAL_VARIABLES) {
store("β₁", atan(tanβ1));
store("β₂", atan(tanβ2));
store("α₁", α1);
α1 = computedToGiven(α1);
}
/*
* Refine α₁ using Newton's method. Karney proposes an hybrid approach: approximate α₁,
* then compute σ₁, σ₂, ω₁ and ω₂ (actually the sine and cosine of those angles in our
* implementation).
*/
double σ1, σ2;
int moreRefinements = Formulas.MAXIMUM_ITERATIONS;
do {
α0(msinα1 = sin(α1), mcosα1 = cos(α1), sinβ1, cosβ1);
final double k2 = computeSeriesExpansionCoefficients();
/*
* Compute α₂ from Karney equation 5: sin(α₀) = sin(α₁)⋅cos(β₁) = sin(α₂)⋅cos(β₂)
* The cos(α₂) term could be computed from sin(α₂), but Karnay recommends instead
* equation 45. An older publication (Karney 2010) went one step further with the
* replacement of (cos²β₂ - cos²β₁) by:
*
* (cosβ₂ - cosβ₁)⋅(cosβ₂ + cosβ₁) if β₁ < -π/4 (Reminder: β₁ is always negative)
* (sinβ₁ - sinβ₂)⋅(sinβ₁ + sinβ₂) otherwise.
*
* Actually we don't need α values directly, but instead cos(α)⋅cos(β).
* Note that cos(α₁) can be negative because α₁ ∈ [0…2π].
*/
final double cosα1_cosβ1 = mcosα1 * cosβ1;
final double cosα2_cosβ2 = sqrt(cosα1_cosβ1*cosα1_cosβ1 +
(cosβ1 <= 1/MathFunctions.SQRT_2 ? (cosβ2 - cosβ1)*(cosβ2 + cosβ1)
: (sinβ1 - sinβ2)*(sinβ1 + sinβ2)));
msinα2 = sinα0;
mcosα2 = cosα2_cosβ2;
/*
* Karney gives the following formulas:
*
* σ = atan2(sinβ, cosα⋅cosβ) — spherical arc length.
* ω = atan2(sin(α₀)⋅sinσ, cosσ) — spherical longitude.
*
* We perform the following substitutions where c is an unknown constant:
* That unknown coefficient will disaspear in atan2(c⋅y, c⋅x) expressions:
*
* sin(σ) = c⋅sinβ
* cos(σ) = c⋅cosα⋅cosβ
*/
final double ω1, ω2;
σ1 = atan2(sinβ1, cosα1_cosβ1); // (Karney 11)
σ2 = atan2(sinβ2, cosα2_cosβ2);
ω1 = atan2(sinβ1*sinα0, cosα1_cosβ1); // (Karney 12) with above-cited substitutions.
ω2 = atan2(sinβ2*sinα0, cosα2_cosβ2);
/*
* Compute the difference in longitude using the current start point and end point.
* If this difference is close enough to the requested accuracy, we are almost done.
* Otherwise refine the results. Note that if we detect that accuracy is good enough,
* we still complete the computation in order to not waste what we have computed so
* far in current iteration.
*/
final double λ2E;
λ1E = sphericalToGeodeticLongitude(ω1, σ1);
λ2E = sphericalToGeodeticLongitude(ω2, σ2);
final double Δλ_error = IEEEremainder(λ2E - λ1E - Δλ, 2*PI);
if (abs(Δλ_error) <= ITERATION_TOLERANCE) {
moreRefinements = 0;
} else if (--moreRefinements == 0) {
throw new GeodeticException(Resources.format(Resources.Keys.NoConvergence));
}
/*
* Special case for α₁ = π/2 and β₂ = ±β₁ (Karney's equation 47). We replace the β₂ = ±β₁
* condition by |β₂| - |β₁| ≈ 0. Assuming tan(θ) ≈ θ for small angles we take the tangent
* of above difference and use tan(β₂ - β₁) = (tanβ₂ - tanβ₁)/(1 + tanβ₂⋅tanβ₁) identity.
* Note that tanβ₁ ≤ 0 and |tanβ₂| ≤ |tanβ₁| in this method.
*/
final double dΔλ_dα1;
if (abs(mcosα1) < LATITUDE_THRESHOLD && (-tanβ1 - abs(tanβ2)) < (1 + abs(tanβ1*tanβ2)) * LATITUDE_THRESHOLD) {
dΔλ_dα1 = -2 * sqrt(1 - eccentricitySquared * (cosβ1*cosβ1)) / sinβ1;
} else {
/*
* Karney's equation 38 combined with equation 46. The substitutions
* for sin(σ) and cos(σ) described above are applied again here.
*/
final double h1 = hypot(sinβ1, cosα1_cosβ1);
final double h2 = hypot(sinβ2, cosα2_cosβ2);
final double sinσ1 = sinβ1 / h1;
final double sinσ2 = sinβ2 / h2;
final double cosσ1 = cosα1_cosβ1 / h1;
final double cosσ2 = cosα2_cosβ2 / h2;
final double J2 = sphericalToEllipsoidalAngle(σ2, true);
final double J1 = sphericalToEllipsoidalAngle(σ1, true);
final double Δm = (sqrt(1 + k2*(sinσ2*sinσ2))*cosσ1*sinσ2
- sqrt(1 + k2*(sinσ1*sinσ1))*sinσ1*cosσ2
- cosσ1*cosσ2*(J2 - J1));
dΔλ_dα1 = Δm * axisRatio / cosα2_cosβ2;
if (STORE_LOCAL_VARIABLES) {
store("J(σ₁)", J1);
store("J(σ₂)", J2);
store("Δm", Δm * semiMinorAxis());
}
}
final double dα1 = Δλ_error / dΔλ_dα1; // Opposite sign of Karney δα₁ term.
/*
* We need to compute α₁ -= dα₁ then iterate again. But sometimes the subtraction has no effect
* because dα₁ ≪ α₁ and iteration continues with unchanged α₁ value until no convergence error.
* If we detect this situation, assume that we have the best accuracy that we can get.
*
* Note: we tried Kahan summation algorithm but it didn't solved the problem.
* No convergence were still happening, but in more indirect ways (after a cycle in iterations).
*/
if (α1 == (α1 -= dα1)) {
moreRefinements = 0;
if (STORE_LOCAL_VARIABLES) {
store("dα₁ ≪ α₁", dα1); // Flag for `iterationReachedPrecisionLimit` in tests.
}
}
if (STORE_LOCAL_VARIABLES) { // For comparing values against Karney table 5 and 6.
final double I1_σ2;
I1_σ1 = sphericalToEllipsoidalAngle(σ1, false); // Required for computation of s₁ in `snapshot()`.
I1_σ2 = sphericalToEllipsoidalAngle(σ2, false);
snapshot();
store("k²", k2);
store("σ₁", σ1);
store("σ₂", σ2);
store("ω₁", ω1);
store("ω₂", ω2);
store("α₂", atan2(msinα2, mcosα2));
store("λ₂", λ2E);
store("Δλ", λ2E - λ1E);
store("δλ", Δλ_error);
store("dλ/dα", dΔλ_dα1);
store("δσ₁", -dα1);
store("α₁", α1);
store("s₂", I1_σ2 * semiMinorAxis());
store("Δs", (I1_σ2 - I1_σ1) * semiMinorAxis());
}
} while (moreRefinements != 0);
final double I1_σ2;
I1_σ1 = sphericalToEllipsoidalAngle(σ1, false);
I1_σ2 = sphericalToEllipsoidalAngle(σ2, false);
geodesicDistance = (I1_σ2 - I1_σ1) * semiMinorAxis();
/*
* Restore the coordinate sign and order to the original configuration.
*/
if (swapPoints) {
double t;
t = msinα1; msinα1 = msinα2; msinα2 = t;
t = mcosα1; mcosα1 = mcosα2; mcosα2 = t;
}
if (inverseLongitudeSigns ^ swapPoints) {
msinα1 = -msinα1;
msinα2 = -msinα2;
}
if (inverseLatitudeSigns ^ swapPoints) {
mcosα1 = -mcosα1;
mcosα2 = -mcosα2;
}
setValid(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE);
if (!(swapPoints | inverseLongitudeSigns | inverseLatitudeSigns)) {
setValid(COEFFICIENTS_FOR_START_POINT);
}
}