Calculate/measure distance between two poins using Vincenty Distance WGS-84 algorithm.

Code translate from geopy http://code.google.com/p/geopy/.

DELIMITER ;;
DROP FUNCTION IF EXISTS `vd`;;
 
CREATE FUNCTION `vd`(lng1 DOUBLE, lat1 DOUBLE, lng2 DOUBLE, lat2 DOUBLE, metric VARCHAR(2)) RETURNS DOUBLE
    DETERMINISTIC
    COMMENT 'Vincenty Distance WGS-84 http://code.google.com/p/geopy/'
BEGIN
DECLARE gcdx DOUBLE;
DECLARE lng_rad1 DOUBLE;
DECLARE lat_rad1 DOUBLE;
DECLARE lng_rad2 DOUBLE;
DECLARE lat_rad2 DOUBLE;
 
DECLARE wgs84_major DOUBLE;
DECLARE wgs84_minor DOUBLE;
DECLARE wgs84_flattening DOUBLE;
 
DECLARE delta_lng DOUBLE;
DECLARE reduced_lat1 DOUBLE;
DECLARE reduced_lat2 DOUBLE;
DECLARE sin_reduced1 DOUBLE;
DECLARE cos_reduced1 DOUBLE;
DECLARE sin_reduced2 DOUBLE;
DECLARE cos_reduced2 DOUBLE;
 
DECLARE lambda_lng DOUBLE;
DECLARE lambda_prime DOUBLE;
 
DECLARE iter_limit INT;
 
DECLARE sin_lambda_lng DOUBLE;
DECLARE cos_lambda_lng DOUBLE;
DECLARE sin_sigma DOUBLE;
DECLARE cos_sigma DOUBLE;
DECLARE sigma DOUBLE;
DECLARE sin_alpha DOUBLE;
DECLARE cos_sq_alpha DOUBLE;
DECLARE cos2_sigma_m DOUBLE;
DECLARE C DOUBLE;
DECLARE u_sq DOUBLE;
DECLARE A DOUBLE;
DECLARE B DOUBLE;
DECLARE delta_sigma DOUBLE;
 
SET lng_rad1 := RADIANS(lng1);
SET lat_rad1 := RADIANS(lat1);
SET lng_rad2 := RADIANS(lng2);
SET lat_rad2 := RADIANS(lat2);
 
SET wgs84_major := 6378.137;
SET wgs84_minor := 6356.7523142;
SET wgs84_flattening := 1 / 298.257223563;
 
SET delta_lng := lng_rad2 - lng_rad1;
 
SET reduced_lat1 := atan((1 - wgs84_flattening) * tan(lat_rad1));
SET reduced_lat2 := atan((1 - wgs84_flattening) * tan(lat_rad2));
 
SET sin_reduced1 := sin(reduced_lat1);
SET cos_reduced1 := cos(reduced_lat1);
SET sin_reduced2 := sin(reduced_lat2);
SET cos_reduced2 := cos(reduced_lat2);
 
SET lambda_lng := delta_lng;
SET lambda_prime := 2 * pi();
 
SET iter_limit = 20;
 
WHILE abs(lambda_lng - lambda_prime) > pow(10, -11) AND iter_limit > 0 DO
     SET sin_lambda_lng := sin(lambda_lng);
     SET cos_lambda_lng := cos(lambda_lng);
 
     SET sin_sigma := sqrt(pow((cos_reduced2 * sin_lambda_lng), 2) +
                      pow((cos_reduced1 * sin_reduced2 - sin_reduced1 *
                       cos_reduced2 * cos_lambda_lng), 2));
 
     IF sin_sigma = 0 THEN
        RETURN 0;
     END IF;
 
     SET cos_sigma := (sin_reduced1 * sin_reduced2 +
                         cos_reduced1 * cos_reduced2 * cos_lambda_lng);
 
     SET sigma := atan2(sin_sigma, cos_sigma);
 
     SET sin_alpha := cos_reduced1 * cos_reduced2 * sin_lambda_lng / sin_sigma;
     SET cos_sq_alpha := 1 - pow(sin_alpha, 2);
 
     IF cos_sq_alpha != 0 THEN
         SET cos2_sigma_m := cos_sigma - 2 * (sin_reduced1 * sin_reduced2 /
                                         cos_sq_alpha);
     ELSE
         SET cos2_sigma_m := 0.0;
     END IF;
 
     SET C := wgs84_flattening / 16.0 * cos_sq_alpha * (4 + wgs84_flattening * (4 - 3 * cos_sq_alpha));
 
     SET lambda_prime := lambda_lng;
     SET lambda_lng := (delta_lng + (1 - C) * wgs84_flattening * sin_alpha *
                   (sigma + C * sin_sigma *
                    (cos2_sigma_m + C * cos_sigma *
                     (-1 + 2 * pow(cos2_sigma_m, 2)))));
     SET iter_limit := iter_limit - 1;
END WHILE;
 
IF iter_limit = 0 THEN
    RETURN NULL;
END IF;
 
SET u_sq := cos_sq_alpha * (pow(wgs84_major, 2) - pow(wgs84_minor, 2)) / pow(wgs84_minor, 2);
 
SET A := 1 + u_sq / 16384.0 * (4096 + u_sq * (-768 + u_sq *
                                        (320 - 175 * u_sq)));
 
SET B := u_sq / 1024.0 * (256 + u_sq * (-128 + u_sq * (74 - 47 * u_sq)));
 
SET delta_sigma := (B * sin_sigma *
               (cos2_sigma_m + B / 4. *
                (cos_sigma * (-1 + 2 * pow(cos2_sigma_m, 2)) -
                 B / 6. * cos2_sigma_m * (-3 + 4 * pow(sin_sigma, 2)) *
                 (-3 + 4 * pow(cos2_sigma_m, 2)))));
 
SET gcdx := wgs84_minor * A * (sigma - delta_sigma);
 
IF metric = 'km' THEN
 RETURN gcdx;
ELSEIF metric = 'mi' THEN
 RETURN gcdx * 0.621371192;
ELSEIF metric = 'nm' THEN
 RETURN gcdx / 1.852;
ELSE
 RETURN gcdx;
END IF;
END;;
 
DELIMITER ;

Current Tags

You must be logged in to tag this tool

No Comments yet