distanceCoordsLambert function

double distanceCoordsLambert(double latADeg, double longADeg, double latBDeg, double longBDeg,)

Calculates the distance in meters between point A at spherical (latitude, longitude) coordinates (latADeg, longADeg) and B at (latBDeg, longBDeg) - all in degrees, using Lambert's formula for maximum accuracy, taking into account the earth's ellipsoid shape.

The result is very accurate, tens of meters over thousands of kilometers. Based on https://www.calculator.net/distance-calculator.html, https://en.wikipedia.org/wiki/Geographical_distance#Lambert's_formula_for_long_lines and https://python.algorithms-library.com/geodesy/lamberts_ellipsoidal_distance.

Implementation

double distanceCoordsLambert(
    double latADeg, double longADeg, double latBDeg, double longBDeg) {
  // Calculate the reduced latitudes.
  final latARad = degToRad(latADeg);
  final beta1 = atan((1 - earthFlattening) * tan(latARad));
  final latBRad = degToRad(latBDeg);
  final beta2 = atan((1 - earthFlattening) * tan(latBRad));

  // The correct value for the radius to use in the haversine calculation is
  // not clear, but the implementations at
  // https://www.calculator.net/distance-calculator.html and
  // https://python.algorithms-library.com/geodesy/lamberts_ellipsoidal_distance
  // use the radius at the equator rather than the mean radius.
  const haversineEarthRadius = EarthRadiusMeters.equatorial;

  // Calculate the haversine distance and with that the central angle sigma.
  final sigma = distanceCoordsHaversine(
          radToDeg(beta1), longADeg, radToDeg(beta2), longBDeg,
          earthRadiusMeters: haversineEarthRadius) /
      haversineEarthRadius;
  final sinSigma = sin(sigma);
  final sinHalfSigma = sin(sigma / 2);
  final cosHalfSigma = cos(sigma / 2);

  // Prevent division by zero further down.
  if (sinHalfSigma == 0 || cosHalfSigma == 0) {
    return 0;
  }

  // Calculate various components of the final formula.
  final p = (beta1 + beta2) / 2;
  final q = (beta2 - beta1) / 2;
  final sinP = sin(p);
  final cosP = cos(p);
  final sinQ = sin(q);
  final cosQ = cos(q);
  final x = (sigma - sinSigma) *
      (sinP * sinP * cosQ * cosQ / (cosHalfSigma * cosHalfSigma));
  final y = (sigma + sinSigma) *
      (cosP * cosP * sinQ * sinQ / (sinHalfSigma * sinHalfSigma));

  // Calculate and return the distance.
  final distance =
      EarthRadiusMeters.equatorial * (sigma - (earthFlattening / 2) * (x + y));
  return distance;
}