Vincenty Direct

CREATE FUNCTION [dbo].[vdd](@lon1 FLOAT, @lat1 FLOAT, @brg FLOAT, @dist FLOAT) 
RETURNS @results TABLE (lon float,lat float, bearing float)
    -- converted to tsql from http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html Chris Veness' Java code by pdc
    -- copyright Chris Veness, pls check his site for license details
    -- pubdc yahoo com
    --
    -- usage : select * from vdd (69.213 , 34.566 , 300, 4500) 
    --
    -- parameters are lon, lat, true bearing, and dist in meters
    -- where lon and lat are coordinates in decimals (as in : 179.9999,89.9999)
    -- to convert from deg/min/secs first do deg+(min/60)+(secs/3600)
    --
    -- the result is a table :
    -- lon                        lat                    bearing
    -- 69.170522            34.586275        300
    --
    -- alternative use in update statement : 
    --    update PEP_raw set 
    -- offxlon = (select lon from vdd(179,89,359,(select ref1.offx from PEP_raw ref1 where ref1.ID = PEP_raw.id))),
    -- offxlat = (select lat from vdd(179,89,359,(select ref2.offx from PEP_raw ref2 where ref2.ID = PEP_raw.id)))
    -- where offx is not null
    --
AS
BEGIN

DECLARE @lon1_rad FLOAT;
DECLARE @lat1_rad FLOAT;
DECLARE @lon2_rad FLOAT;
DECLARE @lat2_rad FLOAT;
DECLARE @lon2 FLOAT;
DECLARE @lat2 FLOAT;

DECLARE @a FLOAT;
DECLARE @b FLOAT;
DECLARE @f FLOAT;
DECLARE @s FLOAT;

DECLARE @alpha1 FLOAT;
DECLARE @sinAlpha1 FLOAT;
DECLARE @cosAlpha1 FLOAT;

DECLARE @tanU1 FLOAT;
DECLARE @cosU1 FLOAT;
DECLARE @sinU1 FLOAT;

DECLARE @sigma1 FLOAT;
DECLARE @sin_alpha FLOAT;
DECLARE @u_sq FLOAT;
DECLARE @cos_sq_alpha FLOAT;
DECLARE @bigA FLOAT;
DECLARE @bigB FLOAT;

DECLARE @sigma FLOAT;
DECLARE @sigmaP FLOAT;
DECLARE @cos2_sigma_m FLOAT;
DECLARE @sin_sigma FLOAT;
DECLARE @cos_sigma FLOAT;
DECLARE @delta_sigma FLOAT;
DECLARE @temp FLOAT;
DECLARE @lambda FLOAT;
DECLARE @C FLOAT;
DECLARE @L FLOAT;
DECLARE @revAz_rad FLOAT;
DECLARE @revAz FLOAT;

SET @lon1_rad = @lon1*PI()/180;
SET @lat1_rad = @lat1*PI()/180;

--WGS-84 ellipsoid
SET @a = 6378137; --wgs84_major
SET @b = 6356752.3142; -- wgs84_minor
SET @f = 1 / 298.257223563; -- wgs84_flattening
SET @s = @dist;
SET @alpha1 = @brg*PI()/180;
SET @sinAlpha1 = SIN(@alpha1);
SET @cosAlpha1 = COS(@alpha1);

SET @tanU1 = (1-@f) * tan(@lat1_rad);
SET @cosU1 = 1/sqrt(1+@tanU1*@tanU1);
SET @sinU1 = @tanU1 * @cosU1;
SET @sigma1 = ATN2(@tanU1,@cosAlpha1);
SET @sin_alpha = @cosU1 * @sinAlpha1;
SET @cos_sq_alpha = 1- (@sin_alpha*@sin_alpha);
SET @u_sq = @cos_sq_alpha * (@a*@a-@b*@b) / (@b*@b);
SET @bigA = 1 + @u_sq/16384*(4096+@u_sq*(-768+@u_sq*(320-175*@u_sq)));
SET @bigB = @u_sq/1024*(256+@u_sq*(-128+@u_sq*(74-47*@u_sq)));

SET @sigma = @s / (@b*@bigA)
SET @sigmaP = 2*PI();

while ABS(@sigma-@sigmaP) > 1e-12 --select cast(1e-12 as numeric(24,18))
    BEGIN
        SET @cos2_sigma_m = COS(2*@sigma1 + @sigma);
        SET @sin_sigma = SIN(@sigma);
        SET @cos_sigma = COS(@sigma);
        SET @delta_sigma = @bigB * @sin_sigma * (@cos2_sigma_m+@bigB/4*(@cos_sigma*(-1+2*@cos2_sigma_m*@cos2_sigma_m)
            - @bigB/6*@cos2_sigma_m*(-3+4*@sin_sigma*@sin_sigma)*(-3+4*@cos2_sigma_m*@cos2_sigma_m)));
        SET @sigmaP = @sigma;
        SET @sigma = @s / (@b*@bigA) + @delta_sigma;
    END

SET @temp = @sinU1*@sin_sigma - @cosU1*@cos_sigma*@cosAlpha1;
SET @lat2_rad = ATN2(@sinU1*@cos_sigma + @cosU1*@sin_sigma*@cosAlpha1,
 (1-@f)*sqrt(@sin_alpha*@sin_alpha + @temp * @temp));
SET @lambda = ATN2(@sin_sigma*@sinAlpha1,
            @cosU1*@cos_sigma - @sinU1*@sin_sigma*@cosAlpha1);
SET @C = @f/16*@cos_sq_alpha*(4+@f*(4-3*@cos_sq_alpha));
SET @L = @lambda - (1-@C) * @f * @sin_alpha *
    (@sigma + @C * @sin_sigma*(@cos2_sigma_m+@C*@cos_sigma*(-1+2*@cos2_sigma_m*@cos2_sigma_m)));
SET @lon2_rad = @lon1_rad+@L+3*PI();
SET @revAz_rad = ATN2(@sin_alpha,-@temp); -- final bearing, if required

SET @lon2 = @lon2_rad*180/PI();
SET @lon2 = @lon2-floor(@lon2/360)*360-180;
SET @lat2 = @lat2_rad*180/PI();
SET @revAz = @revAz_rad*180/PI();
SET @revAz = @revAz-floor(@revAz/360)*360;

insert into @results values(cast(@lon2 as numeric(9,6)),CAST(@lat2 as numeric(9,6)), cast(@revAz as numeric(3,0)));

RETURN

END