vincentyDistance static method

double vincentyDistance(ILatLong latLong1, ILatLong latLong2)

Calculates the geodetic distance in meters between two points using the Vincenty inverse formula for ellipsoids.

This formula is very accurate but computationally more expensive than the Haversine formula.

Implementation

static double vincentyDistance(ILatLong latLong1, ILatLong latLong2) {
  double f = 1 / INVERSE_FLATTENING;
  double L = degToRadian(latLong2.longitude - latLong1.longitude);
  double U1 = atan((1 - f) * tan(degToRadian(latLong1.latitude)));
  double U2 = atan((1 - f) * tan(degToRadian(latLong2.latitude)));
  double sinU1 = sin(U1), cosU1 = cos(U1);
  double sinU2 = sin(U2), cosU2 = cos(U2);

  double lambda = L, lambdaP, iterLimit = 100;

  double cosSqAlpha = 0, sinSigma = 0, cosSigma = 0, cos2SigmaM = 0, sigma = 0, sinLambda = 0, sinAlpha = 0, cosLambda = 0;
  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) return 0; // co-incident points
    cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
    sigma = atan2(sinSigma, cosSigma);
    sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
    cosSqAlpha = 1 - sinAlpha * sinAlpha;
    if (cosSqAlpha != 0) {
      cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
    } else {
      cos2SigmaM = 0;
    }
    double 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() > 1e-12 && --iterLimit > 0);

  if (iterLimit == 0) return 0; // formula failed to converge

  double uSq = cosSqAlpha * (EQUATORIAL_RADIUS * EQUATORIAL_RADIUS - POLAR_RADIUS * POLAR_RADIUS) / POLAR_RADIUS / POLAR_RADIUS;
  double A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
  double B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
  double deltaSigma =
      B *
      sinSigma *
      (cos2SigmaM +
          B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - B / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)));
  double s = POLAR_RADIUS * A * (sigma - deltaSigma);

  return s;
}