vincentyDistance static method

double vincentyDistance(
  1. ILatLong latLong1,
  2. ILatLong latLong2
)

Calculates geodetic distance between two LatLongs using Vincenty inverse formula for ellipsoids. This is very accurate but consumes more resources and time than the sphericalDistance method.

Adaptation of Chriss Veness' JavaScript Code on http://www.movable-type.co.uk/scripts/latlong-vincenty.html

Paper: Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975 (http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf)

@param latLong1 first LatLong @param latLong2 second LatLong @return distance in meters between points as a double

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 * (pow(EQUATORIAL_RADIUS, 2) - pow(POLAR_RADIUS, 2)) / pow(POLAR_RADIUS, 2);
  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;
}