Vincenty Inverse

,

This is a TSQL port of a mysql script from "bramsi"

http://forge.mysql.com/tools/tool.php?id=222

It uses the very accurate Vincenty Inverse formula. I have also ported the Vincenty Direct formula, which calculates the coordinates for a given bearing and distance from a pair of reference coordinates, which I will post here soon..

Enjoy,

Peter De Ceulaer

peter_deceulaer_be

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
-- metric : 'mi', 'nm', anything else result is in 'km'
AS
BEGIN
DECLARE @gcdx 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 @wgs84_major = 6378.137;
SET @wgs84_minor = 6356.7523142;
SET @wgs84_flattening = 1 / 298.257223563;

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

5 (1)

5 (1)