CFLib.org – Common Function Library Project

getVincentyDistance(lat1, lon1, lat2, lon2[, units])

Last updated August 26, 2012

author

Stephen Withington

Version: 1 | Requires: CF9 | Library: CFMLLib

Description:
Vincenty Inverse Solution of Geodesics on the Ellipsoid Original formula found at http://www.movable-type.co.uk/scripts/latlong-vincenty.html This work is licensed under a Create Commons Attribution 3.0 Unported License. http://creativecommons.org/licenses/by/3.0/ from: 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 )

Return Values:
Returns a numeric value that is the distance between the two points

Example:

getVincentyDistance( lat1 = 41.4994954, lon1 = -81.6954088, lat2 = 41.8781136, lon2 = -87.6297982, units = 'sm' )

Parameters:

Name Description Required
lat1 Latitude of first point as a decimal Yes
lon1 Longitude of first point as a decimal Yes
lat2 Latitude of second point as a decimal Yes
lon2 Longitude of second point as a decimal Yes
units Units of measure, one of: m (metres, default), km (kilometres),sm (statutory miles), nm (nautical miles), yd (yards), ft (feet), in (inches), cm (centimetres), mm (millimetres). No

Full UDF Source:

/**
 * Calculates geodetic distance between two points specified by latitude/longitude using Vincenty inverse formula for ellipsoids: http://www.movable-type.co.uk/scripts/latlong-vincenty.html
 * version 0.1 by Stephen Withington
 * version 1.0 by Adam Cameron: throwing exceptions for error conditions instead of returning invalid strings
 * 
 * @param lat1 	 Latitude of first point as a decimal (Required)
 * @param lon1 	 Longitude of first point as a decimal (Required)
 * @param lat2 	 Latitude of second point as a decimal (Required)
 * @param lon2 	 Longitude of second point as a decimal (Required)
 * @param units 	 Units of measure, one of: m (metres, default), km (kilometres),sm (statutory miles), nm (nautical miles), yd (yards), ft (feet), in (inches), cm (centimetres), mm (millimetres). (Optional)
 * @return Returns a numeric value that is the distance between the two points 
 * @author Stephen Withington (stephenwithington@gmail.com) 
 * @version 1, August 26, 2012 
 */
numeric function getVincentyDistance(required numeric lat1, required numeric lon1, required numeric lat2, required numeric lon2, string units="m" ) {
	// WGS-84 ellipsoid params
	var a = 6378137;
	var b = 6356752.314245;
	var f = 1 / 298.257223563;

	var L = ( lon2 - lon1 ) * pi() / 180;
	var U1 = atn( ( 1 - f ) * tan( lat1 * pi() / 180 ) );
	var U2 = atn( ( 1 - f ) * tan( lat2 * pi() / 180 ) );
	var sinU1 = sin( U1 );
	var cosU1 = cos( U1 );
	var sinU2 = sin( U2 );
	var cosU2 = cos( U2 );
	var lambda = L;
	var lambdaP = '';
	var iterLimit = 100;
	var sinLambda = '';
	var cosLambda = '';
	var sinSigma = '';
	var cosSigma = '';
	var sigma = '';
	var sinAlpha = '';
	var cosSqAlpha = '';
	var cos2SigmaM = '';
	var C = '';
	var uSq = '';
	var AA = '';
	var BB = '';
	var deltaSigma = '';
	var s = '';
	var Math = createObject( 'java', 'java.lang.Math' );

	if ( lat1 > 90 || lat1 < -90 || lon1 > 180 || lon1 < -180 || lat2 > 90 || lat2 < -90 || lon2 > 180 || lon2 < -180 ) {
		throw(
			type	= "InvalidCoordinatesException",
			message	= "Invalid coordinates",
			detail	= "Latitude must be between 0deg and +/-90deg. (south latitude is negative). Longitude must be between 0deg and +/-180deg. (west longitude is negative)"
		);
	};

	if ( not structKeyExists( arguments, 'units' ) or not listFindNoCase( 'm,km,sm,nm,yd,ft,in,cm,mm', arguments.units ) ) {
		arguments.units = 'sm'; // default to statute miles
	};

	do {
		sinLambda = sin( lambda );
		cosLambda = cos( lambda );
		sinSigma = sqr( ( 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 = Math.ATan2( JavaCast( 'double', sinSigma ), javaCast( 'double', cosSigma ) ); // CFML doesn't have a native ATan2() method avail.
		sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
		cosSqAlpha = 1 - sinAlpha * sinAlpha;
		cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
		if ( not IsNumeric( cos2SigmaM ) ) { cos2SigmaM = 0; }; // equatorial line: cosSqAlpha=0 (ยง6)
		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 ( abs( lambda - lambdaP ) > 0.000000000001 and --iterLimit > 0 );

	if ( iterLimit == 0 ) {
		throw(
			type	= "FormulaFailureException",
			message	= "Formula failed to converge",
			detail	= "Nearly-antipodal points may fail to give a solution"
		);
	};

	uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
	AA = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
	BB = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
	deltaSigma = BB * sinSigma * ( cos2SigmaM + BB / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) - BB / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
	s = b * AA * ( sigma - deltaSigma );
	s = NumberFormat( s, .999 ); // round to 1mm precision
	
	switch( arguments.units ) {
		case 'm' : // meters 
			break;
		case 'km' : s = s * 0.001; // kilometers
			break;
		case 'nm' : s = s * 0.000539956803; // nautical miles
			break;
		case 'yd' : s = s * 1.0936133; // yards
			break;
		case 'ft' : s = s * 3.2808399; // feet
			break;
		case 'in' : s = s * 39.3700787; // inches
			break;
		case 'cm' : s = s * 100; // centimeters
			break;
		case 'mm' : s = s * 1000; // millimeters
			break;
		default : s = s * 0.000621371192; // statute miles
	};

	return s;
};
blog comments powered by Disqus

Search CFLib.org


Latest Additions

Kevin Cotton added
date2ExcelDate
May 5, 2016

Raymond Camden added
CapFirst
April 25, 2016

Chris Wigginton added
loremIpsum
January 18, 2016

Gary Stanton added
calculateArrival...
November 19, 2015

Sebastiaan Naafs - van Dijk added
getDaysInQuarter
November 13, 2015

Created by Raymond Camden / Design by Justin Johnson