Solved

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

Posted on 2011-02-28
6
475 Views
Last Modified: 2013-11-05
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

0
Comment
Question by:AndyH79
  • 3
  • 3
6 Comments
 
LVL 35

Expert Comment

by:David Todd
ID: 35002350
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
0
 
LVL 35

Accepted Solution

by:
David Todd earned 500 total points
ID: 35002354
PS Suggest coding this as a CLR function instead of using TSQL
0
 

Author Comment

by:AndyH79
ID: 35002892
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.  
0
Three Reasons Why Backup is Strategic

Backup is strategic to your business because your data is strategic to your business. Without backup, your business will fail. This white paper explains why it is vital for you to design and immediately execute a backup strategy to protect 100 percent of your data.

 
LVL 35

Expert Comment

by:David Todd
ID: 35009679
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

0
 

Author Closing Comment

by:AndyH79
ID: 35009787
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());
    }
0
 

Author Comment

by:AndyH79
ID: 35009801
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
0

Featured Post

Optimizing Cloud Backup for Low Bandwidth

With cloud storage prices going down a growing number of SMBs start to use it for backup storage. Unfortunately, business data volume rarely fits the average Internet speed. This article provides an overview of main Internet speed challenges and reveals backup best practices.

Question has a verified solution.

If you are experiencing a similar issue, please ask a related question

Long way back, we had to take help from third party tools in order to encrypt and decrypt data.  Gradually Microsoft understood the need for this feature and started to implement it by building functionality into SQL Server. Finally, with SQL 2008, …
After restoring a Microsoft SQL Server database (.bak) from backup or attaching .mdf file, you may run into "Error '15023' User or role already exists in the current database" when you use the "User Mapping" SQL Management Studio functionality to al…

713 members asked questions and received personalized solutions in the past 7 days.

Join the community of 500,000 technology professionals and ask your questions.

Join & Ask a Question