Vincenty Inverse

CREATE FUNCTION [dbo].[vd](@lng1 FLOAT, @lat1 FLOAT, @lng2 FLOAT, @lat2 FLOAT, @metric VARCHAR(2)) 
RETURNS FLOAT
    -- converted to tsql from http://forge.mysql.com/tools/tool.php?id=222 mysql code by pdc
   -- 'Vincenty Distance WGS-84 http://code.google.com/p/geopy/'
    -- metric : 'mi', 'nm', anything else result is in 'km'
AS
BEGIN
DECLARE @gcdx FLOAT;
DECLARE @lng_rad1 FLOAT;
DECLARE @lat_rad1 FLOAT;
DECLARE @lng_rad2 FLOAT;
DECLARE @lat_rad2 FLOAT;

DECLARE @wgs84_major FLOAT;
DECLARE @wgs84_minor FLOAT;
DECLARE @wgs84_flattening FLOAT;

DECLARE @delta_lng FLOAT;
DECLARE @reduced_lat1 FLOAT;
DECLARE @reduced_lat2 FLOAT;
DECLARE @sin_reduced1 FLOAT;
DECLARE @cos_reduced1 FLOAT;
DECLARE @sin_reduced2 FLOAT;
DECLARE @cos_reduced2 FLOAT;

DECLARE @lambda_lng FLOAT;
DECLARE @lambda_prime FLOAT;

DECLARE @iter_limit INT;

DECLARE @sin_lambda_lng FLOAT;
DECLARE @cos_lambda_lng FLOAT;
DECLARE @sin_sigma FLOAT;
DECLARE @cos_sigma FLOAT;
DECLARE @sigma FLOAT;
DECLARE @sin_alpha FLOAT;
DECLARE @cos_sq_alpha FLOAT;
DECLARE @cos2_sigma_m FLOAT;
DECLARE @C FLOAT;
DECLARE @u_sq FLOAT;
DECLARE @A FLOAT;
DECLARE @B FLOAT;
DECLARE @delta_sigma FLOAT;


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) > 10e-11 and @iter_limit > 0 
BEGIN
 SET @sin_lambda_lng = sin(@lambda_lng);
 SET @cos_lambda_lng = cos(@lambda_lng);

 SET @sin_sigma = sqrt(power(@cos_reduced2 * @sin_lambda_lng,2) + 
            power(@cos_reduced1 * @sin_reduced2 - @sin_reduced1 * @cos_reduced2 * @cos_lambda_lng,2));

 IF @sin_sigma = 0
        BEGIN
 RETURN 0
        END

 SET @cos_sigma = (@sin_reduced1 * @sin_reduced2 +
 @cos_reduced1 * @cos_reduced2 * @cos_lambda_lng);

 SET @sigma = ATN2(@sin_sigma, @cos_sigma);

 SET @sin_alpha = @cos_reduced1 * @cos_reduced2 * @sin_lambda_lng / @sin_sigma;
 SET @cos_sq_alpha = 1 - power(@sin_alpha,2);

 IF @cos_sq_alpha != 0
        BEGIN
 SET @cos2_sigma_m = @cos_sigma - 2 * (@sin_reduced1 * @sin_reduced2 /
 @cos_sq_alpha);
        END
 ELSE
        BEGIN
 SET @cos2_sigma_m = 0.0;
        END

 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 * power(@cos2_sigma_m,2)))));
 SET @iter_limit = @iter_limit - 1;
END

IF @iter_limit = 0
    BEGIN
 RETURN NULL
    END

SET @u_sq = @cos_sq_alpha * (power(@wgs84_major,2) - power(@wgs84_minor,2)) / power(@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 * power(@cos2_sigma_m,2)) -
 @B / 6. * @cos2_sigma_m * (-3 + 4 * power(@sin_sigma,2)) *
 (-3 + 4 * power(@cos2_sigma_m,2)))));

SET @gcdx = @wgs84_minor * @A * (@sigma - @delta_sigma);

IF @metric = 'mi' 
    BEGIN
    SET @gcdx = @gcdx * 0.621371192;
    END
IF @metric = 'nm' 
    BEGIN
    SET @gcdx = @gcdx / 1.852;
    END

RETURN @gcdx;

END