[dba-SQLServer] C# / SQL Server: select all zips within 5 miles of a ZIP list

I wrote an Access application to get a list of all the zips within X miles of another zip, based on 
Lat/Long.  Today my client asked me to perform population counts of all zips within 5 miles of a 
list of zips.

I do not particularly want to do an exhaustive calculation (cartesian product) of each zip in the 
list against every other zip, there are 33K zips in my specific zip table.  The list itself (in this 
instance) contains 400 zips and another that contains 250 zips.  The distance calc has a ton of math 
and doing that math on 12 million lines is probably not going to be particularly speedy.

In thinking about how to narrow down the results, it occurred to me that if i did a preliminary 
calculation on the Lat so that I only pulled other zips within 5 miles of the zip(s) in the list, 
this would narrow things down pretty well.

So I did, and the result is 7K zips within 5 miles of the latitude of the 400 zips in the list. 
Does that make sense.

OK so I have a list of 400 zips, and another list of 7K zips "in the 5 mile latitude band" as those 
400 zips.  That is certainly better.  I could have done the same thing with the longitude.  The 
problem is that the longitude is not a fixed distance apart, but rather starts at zero at the poles 
and becomes greatest at the equator.  Although the max distance between two adjacent longitudes (for 
the US) is going to be found at a specific latitude point in the Hawaiian islands. If I just knew 
what that the distance was at that latitude in Hawaii I could use that factor as well.  But I don't 
know where to find that, though I will Google more.

Anyway...

I ran out and found C# code to perform that calculation.  It would be NICE to do this in TSQL in a 
query right inside of SQL Server.  I found code for both.  Ain't the internet GRAND?

The C# code:

     // Calculate Distance in Milesdouble
     //d = GeoCodeCalc.CalcDistance(47.8131545175277, -122.783203125, 42.0982224111897, -87.890625);
     // Calculate Distance in Kilometersdouble
     //d = GeoCodeCalc.CalcDistance(47.8131545175277, -122.783203125, 42.0982224111897, -87.890625, 
GeoCodeCalcMeasurement.Kilometers);

     //GeoCodeCalc C# Class:

     public static class GeoCodeCalc{
     public const double EarthRadiusInMiles = 3956.0;
     public const double EarthRadiusInKilometers = 6367.0;
     public static double ToRadian(double val) { return val * (Math.PI / 180); }
     public static double DiffRadian(double val1, double val2) { return ToRadian(val2) - 
ToRadian(val1); }
     /// <summary>
     /// Calculate the distance between two geocodes. Defaults to using Miles.
     /// </summary>
     public static double CalcDistance(double lat1, double lng1, double lat2, double lng2) {
     return CalcDistance(lat1, lng1, lat2, lng2, GeoCodeCalcMeasurement.Miles);
     }
      /// <summary>
     /// Calculate the distance between two geocodes.
     /// </summary>
         public static double CalcDistance(double lat1, double lng1, double lat2, double lng2, 
GeoCodeCalcMeasurement m) {
         double radius = GeoCodeCalc.EarthRadiusInMiles;
         if (m == GeoCodeCalcMeasurement.Kilometers) { radius = GeoCodeCalc.EarthRadiusInKilometers; }
         return radius * 2 * Math.Asin( Math.Min(1, Math.Sqrt( ( Math.Pow(Math.Sin((DiffRadian(lat1, 
lat2)) / 2.0), 2.0) + Math.Cos(ToRadian(lat1)) * Math.Cos(ToRadian(lat2)) * 
Math.Pow(Math.Sin((DiffRadian(lng1, lng2)) / 2.0), 2.0) ) ) ) );
         }
     }
     public enum GeoCodeCalcMeasurement : int
     {
     Miles = 0,
     Kilometers = 1
     }


The TSQL code:

Create Function [dbo].[fnGetDistance]
(
       @Lat1 Float(18),
       @Long1 Float(18),
       @Lat2 Float(18),
       @Long2 Float(18),
       @ReturnType VarChar(10)
)

Returns Float(18)

AS

Begin

       Declare @R Float(8);
       Declare @dLat Float(18);
       Declare @dLon Float(18);
       Declare @a Float(18);
       Declare @c Float(18);
       Declare @d Float(18);

       Set @R =
             Case @ReturnType
             When 'Miles' Then 3956.55
             When 'Kilometers' Then 6367.45
             When 'Feet' Then 20890584
             When 'Meters' Then 6367450
             Else 20890584 -- Default feet (Garmin rel elev)
             End

       Set @dLat = Radians(@lat2 - @lat1);

       Set @dLon = Radians(@long2 - @long1);

       Set @a = Sin(@dLat / 2)
                  * Sin(@dLat / 2)
                  + Cos(Radians(@lat1))
                  * Cos(Radians(@lat2))
                  * Sin(@dLon / 2)
                  * Sin(@dLon / 2);
       Set @c = 2 * Asin(Min(Sqrt(@a)));

       Set @d = @R * @c;
       Return @d;

End



-- 
John W. Colby
www.ColbyConsulting.com