MKGeometry+MKGeographicalDistanceBetweenCoordinates.j

if (Number.prototype.toRad === undefined)
Number.prototype.toRad = /* (Number) */ function  () { return this * Math.PI / 180; }

function MKGeographicalDistanceBetweenCoordinates (fromCoords, toCoords) {
	
	if (CLLocationCoordinate2DEqualToCLLocationCoordinate2D(fromCoords, toCoords))
	return 0;
	
	//	Vincenty Inverse Solution.
	//	Google Maps API v3 does not provide distance calculation, so we have to roll our own.
	//	Original: Chris Veness http://www.movable-type.co.uk/scripts/latlong-vincenty.html

	//	Formatting.  Sanitization.
	var	fromLatitude = fromCoords.latitude, fromLongitude = fromCoords.longitude,
		toLatitude = toCoords.latitude, toLongitude = toCoords.longitude;
	var	fromLatitudeRadians = Number(fromLatitude).toRad(), fromLongitudeRadians = Number(fromLongitude).toRad(),
		toLatitudeRadians = Number(toLatitude).toRad(), toLongitudeRadians = Number(toLongitude).toRad();

	//	Ellipsoid Parameters. Using WGS 1984 Data.
	var 	ellipsoidEquatorialAxis = 6378137,
		ellipsoidPolarAxis = 6356752.314245,
		ellipsoidInverseFlattening = 1/298.257223563;		
	var	ellipsoidEquatorialAxisSq = Math.pow(ellipsoidEquatorialAxis, 2),
		ellipsoidPolarAxisSq = Math.pow(ellipsoidPolarAxis, 2);

	//	Difference in Longitude
	var	longitudeDifferenceInRadians = toLongitudeRadians - fromLongitudeRadians;

	//	Reduced Latitude
	var	U1 = Math.atan( (1 - ellipsoidInverseFlattening) * Math.tan(fromLatitudeRadians) ),
		U2 = Math.atan( (1 - ellipsoidInverseFlattening) * Math.tan(toLatitudeRadians) );
		
	var	sinU1 = Math.sin(U1), cosU1 = Math.cos(U1), sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
	var	lambda = longitudeDifferenceInRadians, lambdaP = null, iterationLimit = 100;
	
	
	do {
		
	//	Iterate till lambdaP reached accuracy of 10^-12, approximately 0.06mm

		var sinLambda = Math.sin(lambda), cosLambda = Math.cos(lambda);
		var sinSigma = Math.sqrt(
			
			Math.pow( cosU2 * sinLambda , 2 ) + 
			Math.pow( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda , 2 )

		);
						
		if (sinSigma == 0) return 0;
	//	Distance between co-incident points is zero

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

	//	Equatorial line: cosSqAlpha = 0 (§6)
		if ( isNaN(cos2SigmaM) ) cos2SigmaM = 0;

		var C = ( ellipsoidInverseFlattening / 16 ) * cosSqAlpha * 
			( 4 + ellipsoidInverseFlattening * ( 4 - 3 * cosSqAlpha ) );
		
		lambdaP = lambda;
		
		lambda = longitudeDifferenceInRadians + ( 1 - C ) * ellipsoidInverseFlattening * sinAlpha * 
			( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
				
	} while ( (Math.abs(lambda - lambdaP) > (1e-12)) && ( --iterationLimit > 0 ));
	
	
//	If formula failed to converge, return NaN
	if (iterationLimit == 0) return NaN;

	var uSq = cosSqAlpha * ( ellipsoidEquatorialAxisSq - ellipsoidPolarAxisSq ) / ellipsoidPolarAxisSq;

	var A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
	var B = uSq / 1024 * ( 256 + uSq *( -128 + uSq * ( 74 - 47 * uSq ) ) );

	var deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( 

		cosSigma * ( -1 + 2 * Math.pow(cos2SigmaM, 2) ) -
		B / 6 * cos2SigmaM * ( -3 + 4 * Math.pow(sinSigma, 2) ) * ( -3 + 4 * Math.pow(cos2SigmaM, 2) )

	) );

	var distanceInMeters = ellipsoidPolarAxis * A * ( sigma - deltaSigma );  
	
//	Round to 1mm precision
	return distanceInMeters.toFixed(3); 

}