navispatial
commented
5 days ago Implementation ported from https://github.com/chrisveness/geodesy/blob/master/latlon-ellipsoidal-vincenty.js by Chris Veness (MIT License).
Ported Dart code for direct and inverse Vincenty functions based on the ellipsoidal Earth models.
/// Vincenty direct calculation.
GeodeticArcSegment _direct(double distance, double initialBearing) {
if (distance.isNaN) {
throw throw FormatException('invalid distance $distance');
}
if (distance == 0) {
return GeodeticArcSegment(
origin: origin,
bearing: double.nan,
destination: origin,
finalBearing: double.nan,
distance: 0.0,
);
}
if (initialBearing.isNaN) {
throw throw FormatException('invalid bearing $initialBearing');
}
// symbols: α = alfa, σ = sigma, Δ = delta
final lat1 = origin.lat.toRadians();
final lon1 = origin.lon.toRadians();
final alfa1 = initialBearing.toRadians();
final s = distance;
// ellipsoidal parameters
final a = ellipsoid.a;
final b = ellipsoid.b;
final f = ellipsoid.f;
final sinAlfa1 = sin(alfa1);
final cosAlfa1 = cos(alfa1);
final tanU1 = (1.0 - f) * tan(lat1);
final cosU1 = 1.0 / sqrt(1.0 + tanU1 * tanU1);
final sinU1 = tanU1 * cosU1;
// σ1 = angular distance on the sphere from the equator to P1
final sigma1 = atan2(tanU1, cosAlfa1);
// α = azimuth of the geodesic at the equator
final sinAlfa = cosU1 * sinAlfa1;
final cosSqAlfa = 1.0 - sinAlfa * sinAlfa;
final uSq = cosSqAlfa * (a * a - b * b) / (b * b);
final A = 1.0 +
uSq / 16384.0 * (4096.0 + uSq * (-768.0 + uSq * (320.0 - 175.0 * uSq)));
final B =
uSq / 1024.0 * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq)));
// σ = angular distance P₁ P₂ on the sphere
var sigma = s / (b * A);
double? sinSigma;
double? cosSigma;
// σₘ = angular distance on the sphere from the equator to the midpoint of the line
double? cos2SigmaM;
double? sigmaI;
var iterations = 0;
do {
cos2SigmaM = cos(2 * sigma1 + sigma);
sinSigma = sin(sigma);
cosSigma = cos(sigma);
final deltaSigma = B *
sinSigma *
(cos2SigmaM +
B /
4.0 *
(cosSigma * (-1.0 + 2.0 * cos2SigmaM * cos2SigmaM) -
B /
6.0 *
cos2SigmaM *
(-3.0 + 4.0 * sinSigma * sinSigma) *
(-3.0 + 4.0 * cos2SigmaM * cos2SigmaM)));
sigmaI = sigma;
sigma = s / (b * A) + deltaSigma;
// TV: 'iterate until negligible change in λ' (≈0.006mm)
} while ((sigma - sigmaI).abs() > 1e-12 && ++iterations < 100);
if (iterations >= 100) {
// not possible?
throw const FormatException('Vincenty formula failed to converge');
}
final x = sinU1 * sinSigma - cosU1 * cosSigma * cosAlfa1;
final lat2 = atan2(
sinU1 * cosSigma + cosU1 * sinSigma * cosAlfa1,
(1.0 - f) * sqrt(sinAlfa * sinAlfa + x * x),
);
final lon = atan2(
sinSigma * sinAlfa1,
cosU1 * cosSigma - sinU1 * sinSigma * cosAlfa1,
);
final C = f / 16.0 * cosSqAlfa * (4.0 + f * (4.0 - 3.0 * cosSqAlfa));
final L = lon -
(1.0 - C) *
f *
sinAlfa *
(sigma +
C *
sinSigma *
(cos2SigmaM +
C * cosSigma * (-1.0 + 2.0 * cos2SigmaM * cos2SigmaM)));
final lon2 = lon1 + L;
final alfa2 = atan2(sinAlfa, -x);
return GeodeticArcSegment(
origin: origin,
bearing: initialBearing,
distance: distance,
finalBearing: alfa2.toDegrees().wrap360(),
destination: Geographic(lat: lat2.toDegrees(), lon: lon2.toDegrees()),
);
}
/// Vincenty inverse calculation.
GeodeticArcSegment _inverse(Geographic destination) {
final lat1 = origin.lat.toRadians();
final lon1 = origin.lon.toRadians();
final lat2 = destination.lat.toRadians();
final lon2 = destination.lon.toRadians();
// symbols: α = alfa, σ = sigma, Δ = delta
// ellipsoidal parameters
final a = ellipsoid.a;
final b = ellipsoid.b;
final f = ellipsoid.f;
// L = difference in longitude, U = reduced latitude, defined by tan U = (1-f)·tanφ.
final L = lon2 - lon1;
final tanU1 = (1.0 - f) * tan(lat1);
final cosU1 = 1.0 / sqrt(1.0 + tanU1 * tanU1);
final sinU1 = tanU1 * cosU1;
final tanU2 = (1.0 - f) * tan(lat2);
final cosU2 = 1.0 / sqrt(1.0 + tanU2 * tanU2);
final sinU2 = tanU2 * cosU2;
final antipodal = L.abs() > pi / 2.0 || (lat2 - lat1).abs() > pi / 2.0;
// λ = difference in longitude on an auxiliary sphere
var lon = L;
double? sinLon;
double? cosLon;
// σ = angular distance P₁ P₂ on the sphere
var sigma = antipodal ? pi : 0.0;
var sinSigma = 0.0;
var cosSigma = antipodal ? -1.0 : 1.0;
double? sinSqSigma;
// σₘ = angular distance on the sphere from the equator to the midpoint of the line
var cos2SigmaM = 1.0;
// α = azimuth of the geodesic at the equator
var cosSqAlfa = 1.0;
double? lonI;
var iterations = 0;
do {
sinLon = sin(lon);
cosLon = cos(lon);
sinSqSigma = (cosU2 * sinLon) * (cosU2 * sinLon) +
(cosU1 * sinU2 - sinU1 * cosU2 * cosLon) *
(cosU1 * sinU2 - sinU1 * cosU2 * cosLon);
if (sinSqSigma.abs() < 1e-24) {
// co-incident/antipodal points (σ < ≈0.006mm)
break;
}
sinSigma = sqrt(sinSqSigma);
cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLon;
sigma = atan2(sinSigma, cosSigma);
final sinAlfa = cosU1 * cosU2 * sinLon / sinSigma;
cosSqAlfa = 1.0 - sinAlfa * sinAlfa;
// on equatorial line cos²α = 0 (§6)
cos2SigmaM = (cosSqAlfa != 0.0)
? (cosSigma - 2.0 * sinU1 * sinU2 / cosSqAlfa)
: 0.0;
final C = f / 16.0 * cosSqAlfa * (4.0 + f * (4.0 - 3.0 * cosSqAlfa));
lonI = lon;
lon = L +
(1.0 - C) *
f *
sinAlfa *
(sigma +
C *
sinSigma *
(cos2SigmaM +
C *
cosSigma *
(-1.0 + 2.0 * cos2SigmaM * cos2SigmaM)));
final iterationCheck = antipodal ? lon.abs() - pi : lon.abs();
if (iterationCheck > pi) {
throw const FormatException('λ > π');
}
// TV: 'iterate until negligible change in λ' (≈0.006mm)
} while ((lon - lonI).abs() > 1e-12 && ++iterations < 1000);
if (iterations >= 1000) {
throw const FormatException('Vincenty formula failed to converge');
}
final uSq = cosSqAlfa * (a * a - b * b) / (b * b);
final A = 1.0 +
uSq / 16384.0 * (4096.0 + uSq * (-768.0 + uSq * (320.0 - 175.0 * uSq)));
final B =
uSq / 1024.0 * (256.0 + uSq * (-128.0 + uSq * (74.0 - 47.0 * uSq)));
final deltaSigma = B *
sinSigma *
(cos2SigmaM +
B /
4.0 *
(cosSigma * (-1.0 + 2.0 * cos2SigmaM * cos2SigmaM) -
B /
6.0 *
cos2SigmaM *
(-3.0 + 4.0 * sinSigma * sinSigma) *
(-3.0 + 4.0 * cos2SigmaM * cos2SigmaM)));
// s = length of the geodesic
final s = b * A * (sigma - deltaSigma);
// note special handling of exactly antipodal points where sin²σ = 0 (due to
// discontinuity)
//
// atan2(0, 0) = 0 but atan2(ε, 0) = π/2 / 90°) - in which case bearing is
// always meridional, due north (or due south!)
// α = azimuths of the geodesic; α2 the direction P₁ P₂ produced
final epsilon = doublePrecisionEpsilon;
final alfa1 = sinSqSigma.abs() < epsilon
? 0.0
: atan2(cosU2 * sinLon, cosU1 * sinU2 - sinU1 * cosU2 * cosLon);
final alfa2 = sinSqSigma.abs() < epsilon
? pi
: atan2(cosU1 * sinLon, -sinU1 * cosU2 + cosU1 * sinU2 * cosLon);
return GeodeticArcSegment(
origin: origin,
destination: destination,
distance: s,
bearing: s.abs() < epsilon ? double.nan : alfa1.toDegrees().wrap360(),
finalBearing:
s.abs() < epsilon ? double.nan : alfa2.toDegrees().wrap360(),
);
}
Based on #242 and #243 add also ellipsoidal Vincenty calculations (distances and bearings) between geographic points.
Also use same interfaces as previous (already published) spherical geodesy functions based on rhumb line and great circle.