Link to home
Start Free TrialLog in
Avatar of AndyH79
AndyH79Flag for United States of America

asked on

Conversion from Excel VBA to T-SQL not generating the same results.

I'm needing to convert the following code from Excel VBA to T-SQL. The problem that I am running into is that I can't get the system to generate the same result when run in T-SQL as when run in VBA. Both processes should return a result of 193585.85. Currently the T-SQL process is returning a result of -99 which is due to the debugging process identifying that all 20 iterations have occurred prior to a result being generated which makes Abs(@lambda - @lambdaP) > @EPSILON = FALSE

Since the VBA code uses a constant value for Pi I have used the same constant value as a declared variable rather than use the PI() SQL function since I'm unsure of the number of significant digits the SQL function uses.

The first issue that I run into is the calculation of:
@U1
Excel VBA code Line 62         U1 = Atn((1 - f) * Tan(toRad(lat1)))
returns   0.791295464197368
SQL Server code Line 68     SET @U1 = Atan((1 - @f) * Tan(RADIANS(@lat1)))
returns   0.791295464197410

@U2
Excel VBA code Line 62         U2 = Atn((1 - f) * Tan(toRad(lat2)))
returns   0.821498274972993
SQL Server code Line 69 --  SET @U2 = Atan((1 - @f) * Tan(RADIANS(@lat2)))
returns   0.821498274973035

Neither difference seems to be significant,  I'm just not sure what is causing them to happen. Originally I was also manually converting the degrees to radians in SQL Server using the following function. Since the result was the same using the Radian() sql function I removed the references to this udf in the code.

CREATE FUNCTION [dbo].[udfMath_DegreesToRadian]
(
      -- Add the parameters for the function here
      @degrees Decimal(38,30)
)
RETURNS Decimal(38,30)
AS
BEGIN
      DECLARE @res Decimal(38,30)
      DECLARE @PI Decimal(38,30)
      SET @PI = 3.14159265358979
      SELECT @res = (@degrees * (@PI /180))
      
      RETURN @res
END


Complete VBA Code

Option Explicit

Private Const PI = 3.14159265358979
Private Const EPSILON As Double = 0.000000000001

Public Sub test()
    Dim metersDistance As Double
    metersDistance = Geo.GetDistance(45.434094, -122.5791498, 47.1643339, -122.293198)
End Sub

Public Function GetDistance(ByVal lat1 As Double, ByVal lon1 As Double, ByVal lat2 As Double, ByVal lon2 As Double) As Double
'=================================================================================
' Calculate geodesic distance (in m) between two points specified by latitude/longitude (in numeric degrees)
' using Vincenty inverse formula for ellipsoids

    Dim low_a As Double
    Dim low_b As Double
    Dim f As Double
    Dim L As Double
    Dim U1 As Double
    Dim U2 As Double
    Dim sinU1 As Double
    Dim sinU2 As Double
    Dim cosU1 As Double
    Dim cosU2 As Double
    Dim lambda As Double
      Dim lambdaSub1 As Double
      Dim lambdaSub2 As Double
      Dim lambdaSub3 As Double
   Dim lambdaP As Double
    Dim iterLimit As Integer
    Dim sinLambda As Double
    Dim cosLambda As Double
    Dim sinSigma As Double
    Dim cosSigma As Double
    Dim sigma As Double
    
    Dim sinAlpha As Double
    Dim cosSqAlpha As Double
    Dim cos2SigmaM As Double
    Dim C As Double
    Dim uSq As Double
    Dim upper_A As Double
      Dim upper_ASub1 As Double
      Dim upper_ASub2 As Double
      Dim upper_ASub3 As Double
    Dim upper_B As Double
      Dim upper_BSub1 As Double
      Dim upper_BSub2 As Double
      Dim upper_BSub3 As Double
    Dim deltaSigma As Double
      Dim deltaSigmaSub1 As Double
      Dim deltaSigmaSub2 As Double
      Dim deltaSigmaSub3 As Double
    Dim s As Double
    
    low_a = 6378137
    low_b = 6356752.31414
    f = 1 / 298.257222101 'GRS-80 ellipsiod
    
    L = toRad(lon2 - lon1)
    U1 = Atn((1 - f) * Tan(toRad(lat1)))
    U2 = Atn((1 - f) * Tan(toRad(lat2)))
    sinU1 = Sin(U1)
    cosU1 = Cos(U1)
    sinU2 = Sin(U2)
    cosU2 = Cos(U2)
    
    lambda = L
    lambdaP = 2 * PI
    iterLimit = 20
    
     While (Abs(lambda - lambdaP) > EPSILON) And (iterLimit > 0)
        iterLimit = iterLimit - 1
        
        sinLambda = Sin(lambda)
        cosLambda = Cos(lambda)
        sinSigma = Sqr(((cosU2 * sinLambda) ^ 2) + ((cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) ^ 2))
        If sinSigma = 0 Then
            GetDistance = 0  'co-incident points
            Exit Function
        End If
        cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
        sigma = Atan2(cosSigma, sinSigma)
        sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
        cosSqAlpha = 1 - sinAlpha * sinAlpha

        If cosSqAlpha = 0 Then 'check for a divide by zero
            cos2SigmaM = 0 '2 points on the equator
        Else
            cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha
        End If
        
        C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
        lambdaP = lambda
        'Had to split the calc of lambda since excel says is too complex
        'lambda = L + (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * (cos2SigmaM ^ 2))))
        lambdaSub1 = (1 - C) * f * sinAlpha
        lambdaSub2 = (-1 + 2 * (cos2SigmaM ^ 2))
        lambdaSub3 = (cos2SigmaM + C * cosSigma * lambdaSub2)
        lambda = L + lambdaSub1 * (sigma + C * sinSigma * lambdaSub3)
     Wend
    
    If iterLimit < 1 Then
        MsgBox "iteration limit has been reached, something didn't work."
        Exit Function
    End If

    uSq = cosSqAlpha * (low_a ^ 2 - low_b ^ 2) / (low_b ^ 2)
    'Had to split the calc of upper_A since excel says is too complex
    'upper_A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
    upper_ASub1 = 320 - 175 * uSq
    upper_ASub2 = (-768 + uSq * (upper_ASub1))
    upper_ASub3 = (4096 + uSq * upper_ASub2)
    upper_A = 1 + uSq / 16384 * upper_ASub3
    'Had to split the calc of upper_B since excel says is too complex
    'upper_B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
    upper_BSub1 = (74 - 47 * uSq)
    upper_BSub2 = (-128 + uSq * upper_BSub1)
    upper_BSub3 = (256 + uSq * upper_BSub2)
    upper_B = uSq / 1024 * upper_BSub3
    'Had to split the calc of deltaSigma since excel says is too complex
    'deltaSigma = upper_B * sinSigma * (cos2SigmaM + upper_B / 4 * (cosSigma * (-1 + 2 * cos2SigmaM ^ 2) - upper_B / 6 * cos2SigmaM * (-3 + 4 * sinSigma ^ 2) * (-3 + 4 * cos2SigmaM ^ 2))
    deltaSigmaSub1 = (-3 + 4 * sinSigma ^ 2) * (-3 + 4 * cos2SigmaM ^ 2)
    deltaSigmaSub2 = (cosSigma * (-1 + 2 * cos2SigmaM ^ 2) - upper_B / 6 * cos2SigmaM * deltaSigmaSub1)
    deltaSigmaSub3 = (cos2SigmaM + upper_B / 4 * deltaSigmaSub2)
    deltaSigma = upper_B * sinSigma * deltaSigmaSub3
    s = low_b * upper_A * (sigma - deltaSigma)
    GetDistance = Round(s, 3) 'Rounds distance to nearest millimeter.

End Function

Private Function toRad(ByVal degrees As Double) As Double
    toRad = degrees * (PI / 180)
End Function

Private Function Atan2(ByVal X As Double, ByVal y As Double) As Double
   If y > 0 Then
      If X >= y Then
         Atan2 = Atn(y / X)
      ElseIf X <= -y Then
         Atan2 = Atn(y / X) + PI
      Else
         Atan2 = PI / 2 - Atn(X / y)
      End If
   Else
      If X >= -y Then
         Atan2 = Atn(y / X)
      ElseIf X <= y Then
         Atan2 = Atn(y / X) - PI
      Else
         Atan2 = -Atn(X / y) - PI / 2
      End If
   End If
End Function

Open in new window


Complete T-SQL Code:

DECLARE @lat1 FLOAT
DECLARE @lon1 FLOAT
DECLARE @lat2 FLOAT
DECLARE @lon2 FLOAT

-- 0027		
SET @lat1 = 45.4340940
SET @lon1 = -122.5791498
-- 0601
SET @lat2 = 47.1643339
SET @lon2 =-122.2931980

DECLARE @EPSILON FLOAT
DECLARE @PI FLOAT
SET @EPSILON = 0.000000000001
SET @PI = 3.14159265358979

DECLARE @low_a FLOAT,
		@low_b FLOAT,
		@f FLOAT,
		@L FLOAT,
		@U1 FLOAT,
		@U2 FLOAT,
		@sinU1 FLOAT,
		@sinU2 FLOAT,
		@cosU1 FLOAT,
		@cosU2 FLOAT,
		@lambda FLOAT,
			@lambdaSub1 FLOAT,
			@lambdaSub2 FLOAT,
			@lambdaSub3 FLOAT,
		@lambdaP FLOAT,
		@iterLimit TinyInt,
		@sinLambda FLOAT,
		@cosLambda FLOAT,
		@sinSigma FLOAT,
		@cosSigma FLOAT,
		@sigma FLOAT,
		
		@sinAlpha FLOAT,
		@cosSqAlpha FLOAT,
		@cos2SigmaM FLOAT,
		@C FLOAT,
		@uSq FLOAT,
		@upper_A FLOAT,
			@upper_ASub1 FLOAT,
			@upper_ASub2 FLOAT,
			@upper_ASub3 FLOAT,
		@upper_B FLOAT,
			@upper_BSub1 FLOAT,
			@upper_BSub2 FLOAT,
			@upper_BSub3 FLOAT,
		@deltaSigma FLOAT,
			@deltaSigmaSub1 FLOAT,
			@deltaSigmaSub2 FLOAT,
			@deltaSigmaSub3 FLOAT,
		@s FLOAT


DECLARE @res FLOAT


SET @low_a = 6378137
SET @low_b = 6356752.31414
SET @f = 1 / 298.257222101 --GRS-80 ellipsiod
    
SET @L = RADIANS(@lon2 - @lon1)
SET @U1 = Atan((1 - @f) * Tan(RADIANS(@lat1)))
SET @U2 = Atan((1 - @f) * Tan(RADIANS(@lat2)))

--Included for debugging purposes
SELECT @L as L, @U1 as U1, @U2 as U2

SET @sinU1 = Sin(@U1)
SET @cosU1 = Cos(@U1)
SET @sinU2 = Sin(@U2)
SET @cosU2 = Cos(@U2)
    
SET @lambda = @L
SET @lambdaP = 2 * @PI
SET @iterLimit = 20

--Included for debugging purposes
--SELECT @sinU1 as sinU1, @cosU1 as cosU1, @sinU2 as sinU2,@cosU2 as cosU2, @lambda as lambda, @lambdaP as lambdaP, @iterLimit as iterLimit

While (Abs(@lambda - @lambdaP) > @EPSILON) AND (@iterLimit > 0)
	BEGIN
	SET @iterLimit = @iterLimit - 1
    
   SET @sinLambda = Sin(@lambda)
   SET @cosLambda = Cos(@lambda)
   SET @sinSigma = Square((Square((@cosU2 * @sinLambda))) + (Square(@cosU1 * @sinU2 - @sinU1 * @cosU2 * @cosLambda)))
   IF @sinSigma = 0
      BEGIN  
        SELECT @iterLimit as Iterations, -98  AS ErrorCode--co-incident points
        --RETURN @res
		END
	SET @cosSigma = @sinU1 * @sinU2 + @cosU1 * @cosU2 * @cosLambda
   SET @sigma = atn2(@cosSigma, @sinSigma)
   SET @sinAlpha = @cosU1 * @cosU2 * @sinLambda / @sinSigma
   SET @cosSqAlpha = 1 - @sinAlpha * @sinAlpha

   If @cosSqAlpha = 0 --check for a divide by zero
       Begin
       SET @cos2SigmaM = 0 --2 points on the equator
       END
   Else
		Begin
        SET @cos2SigmaM = @cosSigma - 2 * @sinU1 * @sinU2 / @cosSqAlpha
		End
    
   SET @C = @f / 16 * @cosSqAlpha * (4 + @f * (4 - 3 * @cosSqAlpha))
   SET @lambdaP = @lambda
   --Had to split the calc of lambda since excel says is too complex
   --lambda = L + (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * (Power(cos2SigmaM,2)))))
   SET @lambdaSub1 = (1 - @C) * @f * @sinAlpha
   SET @lambdaSub2 = (-1 + 2 * (Square(@cos2SigmaM)))
   SET @lambdaSub3 = (@cos2SigmaM + @C * @cosSigma * @lambdaSub2)
   SET @lambda = @L + @lambdaSub1 * (@sigma + @C * @sinSigma * @lambdaSub3)
   END
   
   If @iterLimit < 1
		BEGIN
			SET @res = -99 AS ErrorCode
			SELECT @res --iteration limit has been reached, something didn't work."
			--RETURN @res
		END

    SET @uSq = @cosSqAlpha * (Square(@low_a) - Square(@low_b) / Square(@low_b))
    --Had to split the calc of upper_A since excel says is too complex
    --upper_A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
    SET @upper_ASub1 = 320 - 175 * @uSq
    SET @upper_ASub2 = (-768 + @uSq * (@upper_ASub1))
    SET @upper_ASub3 = (4096 + @uSq * @upper_ASub2)
    SET @upper_A = 1 + @uSq / 16384 * @upper_ASub3
    --Had to split the calc of upper_B since excel says is too complex
    --upper_B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
    SET @upper_BSub1 = (74 - 47 * @uSq)
    SET @upper_BSub2 = (-128 + @uSq * @upper_BSub1)
    SET @upper_BSub3 = (256 + @uSq * @upper_BSub2)
    SET @upper_B = @uSq / 1024 * @upper_BSub3
    --Had to split the calc of deltaSigma since excel says is too complex
    --deltaSigma = upper_B * sinSigma * (cos2SigmaM + upper_B / 4 * (cosSigma * (-1 + 2 * Power(cos2SigmaM,2)) - upper_B / 6 * cos2SigmaM * (-3 + 4 * Power(sinSigma,2)) * (-3 + 4 * Power(cos2SigmaM,2)))
    SET @deltaSigmaSub1 = (-3 + 4 * Square(@sinSigma)) * (-3 + 4 * Square(@cos2SigmaM))
    SET @deltaSigmaSub2 = (@cosSigma * (-1 + 2 * Square(@cos2SigmaM)) - @upper_B / 6 * @cos2SigmaM * @deltaSigmaSub1)
    SET @deltaSigmaSub3 = (@cos2SigmaM + @upper_B / 4 * @deltaSigmaSub2)
    SET @deltaSigma = @upper_B * @sinSigma * @deltaSigmaSub3
    SET @s = @low_b * @upper_A * (@sigma - @deltaSigma)
    
    If @res <> -98 and @res <> -99
		BEGIN
			SET @res = Round(@s, 3)--Rounds distance to nearest millimeter.
			SELECT @res 
		END

Open in new window

Avatar of David Todd
David Todd
Flag of New Zealand image

Hi,

In SQL PI is
3.1415926535897931

so try that in your VBA code
select convert( decimal( 38, 30 ), pi()) as pie_as_decimal, convert( float( 53 ), pi()) as pie_as_float;

Your differences start in the 13th digit and to me this is the difference in the two values of pi.

HTH
  David
ASKER CERTIFIED SOLUTION
Avatar of David Todd
David Todd
Flag of New Zealand image

Link to home
membership
This solution is only available to members.
To access this solution, you must be a member of Experts Exchange.
Start Free Trial
Avatar of AndyH79

ASKER

Both the SQL and VBA code make use of a locally declared value for PI = 3.14159265358979. For sql it is @PI. As a test when I update the VBA PI constant to 3.1415926535897931 it only results in a change in value of 0.0000000000000001where the values for @U1 and @U2 are 0.000000000000042 different from their VBA counterparts. Could this be due to some rounding error because I'm using float values instead of decimal values? I may look at doing this as a CLR in the future.  
Hi,

I've been working on an iterative process-calculation. What I found is that is division using decimal only returns 6 decimal places, so I'm using float, which equates to double in VBA.

Potentially there maybe a little difference in the returns from tan and arctan functions. Have you checked them individually?

Something to try is in every place you want float values returned specify float constants ie 180.0 instead of 180 and see if that makes a difference. It made a significant difference at one point for my project.

Regards
  David

Avatar of AndyH79

ASKER

I took your advice and converted it into a CLR function.

[Microsoft.SqlServer.Server.SqlFunction]
    public static SqlString udfMath_InverseSolutionOfGeodesics
        (
        Double lat1,
        Double lon1,
        Double lat2,
        Double lon2
        )
    {
        //Calculate geodesic distance (in m) between two points specified by latitude/longitude (in numeric degrees)
        //using Vincenty inverse formula for ellipsoids
        Double result;
        Double a;
        Double b;
        Double f;
        Double L;
        Double U1;
        Double U2;
        Double sinU1;
        Double sinU2;
        Double cosU1;
        Double cosU2;
        Double lambda;
        Double lambdaP;
        int iterLimit;
        Double sinLambda;
        Double cosLambda;
        Double sinSigma = 0;
        Double cosSigma = 0;
        Double sigma = 0;
   
        Double sinAlpha;
        Double cosSqAlpha = 0;
        Double cos2SigmaM = 0;
        Double C;
        Double uSq;
        Double A;
        Double B;
        Double deltaSigma;
        Double s;
       
        Double PI = 3.14159265358979;
        Double EPSILON = 0.000000000001;

    a = 6378137;
    b = 6356752.31414;
    f = 1 / 298.257222101; //GRS-80 ellipsiod
    L = PI * (lon2-lon1) / 180.0;
    U1 = Math.Atan((1 - f) * Math.Tan(PI * (lat1) / 180.0));
    U2 = Math.Atan((1 - f) * Math.Tan(PI * (lat2) / 180.0));
    sinU1 = Math.Sin(U1);
    cosU1 = Math.Cos(U1);
    sinU2 = Math.Sin(U2);
    cosU2 = Math.Cos(U2);
   
    lambda = L;
    lambdaP = 2 * PI;
    iterLimit = 20;
   
     while (Math.Abs(lambda - lambdaP) > EPSILON && iterLimit > 0)
     {
        iterLimit = iterLimit - 1;
       
        sinLambda = Math.Sin(lambda);
        cosLambda = Math.Cos(lambda);
        sinSigma = Math.Sqrt(Math.Pow((cosU2 * sinLambda),2) + Math.Pow((cosU1 * sinU2 - sinU1 * cosU2 * cosLambda),2));
        if (sinSigma == 0)
        {
            result = -98;  //co-incident points
        }
        cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
        sigma = Math.Atan2(sinSigma,cosSigma);
        sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
        cosSqAlpha = 1 - sinAlpha * sinAlpha;

        if (cosSqAlpha == 0) //check for a divide by zero
        {
            cos2SigmaM = 0; //2 points on the equator
        }
        else
        {
            cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
        }
       
        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 * Math.Pow(cos2SigmaM,2))));
    }
   
    if (iterLimit < 1)
       {
        result = -99; //iteration limit has been reached, something didn't work.
       }


    uSq = cosSqAlpha * (Math.Pow(a, 2) - Math.Pow(b, 2)) / Math.Pow(b, 2);
    A = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
    B = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
    deltaSigma = B * sinSigma * (cos2SigmaM + B / 4 * (cosSigma * (-1 + 2 * (Math.Pow(cos2SigmaM,2) - B / 6 * cos2SigmaM * (-3 + 4 * Math.Pow(sinSigma,2) * (-3 + 4 * Math.Pow(cos2SigmaM,2)))))));
    s = b * A * (sigma - deltaSigma);
    result = Math.Round(s, 3); //Rounds distance to nearest millimeter
    return new SqlString(result.ToString());
    }
Avatar of AndyH79

ASKER

This is used for determining the distance between to coordinates in meters.  Credit for the original code goes to http://www.movable-type.co.uk/scripts/latlong-vincenty.html