distance method

Implementation

@override
double distance(
  final LatLng p1,
  final LatLng p2, {
  final LongitudeDirection lngDir = LongitudeDirection.lazy,
}) {
  const a = equatorRadius,
      b = polarRadius,
      f = flattening; // WGS-84 ellipsoid params

  final effectiveDLng = lngDir.effectiveLongitudinalDelta(p1, p2);

  // Vincenty's solver only converges for |l| <= pi (short arc).
  // Long-arc cases fold l back into (−pi, pi] here and compensate at the end.
  final isLongArc = effectiveDLng.abs() > pi;
  final l = isLongArc
      ? effectiveDLng - (effectiveDLng > 0 ? tau : -tau)
      : effectiveDLng;

  final u1 = atan((1 - f) * tan(p1.latitudeInRad));
  final u2 = atan((1 - f) * tan(p2.latitudeInRad));
  final sinU1 = sin(u1), cosU1 = cos(u1);
  final sinU2 = sin(u2), cosU2 = cos(u2);

  double sinLambda,
      cosLambda,
      sinSigma,
      cosSigma,
      sigma,
      sinAlpha,
      cosSqAlpha,
      cos2SigmaM;
  double lambda = l, lambdaP;
  var maxIterations = this.maxIterations;

  do {
    sinLambda = sin(lambda);
    cosLambda = cos(lambda);
    sinSigma = sqrt((cosU2 * sinLambda) * (cosU2 * sinLambda) +
        (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) *
            (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda));

    if (sinSigma == 0) {
      // Co-incident points (l folded to 0). For long-arc directions this
      // means the caller asked for a full geodesic loop; return the
      // equatorial circumference as the best available scalar approximation.
      return isLongArc ? tau * equatorRadius : 0.0;
    }

    cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
    sigma = atan2(sinSigma, cosSigma);
    sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
    cosSqAlpha = 1 - sinAlpha * sinAlpha;
    cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;

    if (cos2SigmaM.isNaN) {
      cos2SigmaM = 0.0; // equatorial line: cosSqAlpha=0 (§6)
    }

    final C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
    lambdaP = lambda;
    lambda = l +
        (1 - C) *
            f *
            sinAlpha *
            (sigma +
                C *
                    sinSigma *
                    (cos2SigmaM +
                        C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
  } while ((lambda - lambdaP).abs() > accuracy && --maxIterations > 0);

  if (maxIterations == 0) {
    throw StateError('Distance calculation failed to converge!');
  }

  final uSq = cosSqAlpha * (a * a - b * b) / (b * b);
  final A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
  final B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
  final deltaSigma = B *
      sinSigma *
      (cos2SigmaM +
          B /
              4 *
              (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
                  B /
                      6 *
                      cos2SigmaM *
                      (-3 + 4 * sinSigma * sinSigma) *
                      (-3 + 4 * cos2SigmaM * cos2SigmaM)));

  final shortDist = b * A * (sigma - deltaSigma);

  if (!isLongArc) return shortDist;

  // The short and long arcs of a geodesic share the same A (it depends on
  // cosSqAlpha, a property of the geodesic, not the arc direction).
  // Total geodesic length = b*A*2pi, so long arc = total − short arc.
  // deltaSigma for the complementary arc would differ slightly, but the
  // discrepancy is below 1 mm on WGS-84 for any pair of points.
  return b * A * tau - shortDist;
}