ciudades y distancia por latitud-longitud

Aquí está el nuestro. Es posible que deba modificarlo para la estructura de su tabla. El nuestro busca ubicaciones minoristas (y servicios), no ciudades, pero la parte difícil es el "más cercano por distancia" que funciona en esta declaración.

CREATE PROCEDURE [dbo].[GetNearbyLocations] @CenterLatitude FLOAT, @CenterLongitude FLOAT
AS

DECLARE @CntXAxis FLOAT 
DECLARE @CntYAxis FLOAT 
DECLARE @CntZAxis FLOAT 

SET @CntXAxis = COS(RADIANS(@CenterLatitude)) * COS(RADIANS(@CenterLongitude)) 
SET @CntYAxis = COS(RADIANS(@CenterLatitude)) * SIN(RADIANS(@CenterLongitude)) 
SET @CntZAxis = SIN(RADIANS(@CenterLatitude)) 

SELECT LocationId, LocationName, Address, City, State, Zip, Phone, Latitude, Longitude,
               hasATM, hasCarWash, hasDiesel, hasE85, is24hr, hasTrendar, hasWiFi, isTravelCenter, isMiniTravelCenter, isTruckerFriendly, hasScale, hasHotFood,
               ProxDistance = 3961 * ACOS( dbo.XAxis(latitude, longitude)*@CntXAxis + dbo.YAxis(latitude, longitude)*@CntYAxis + dbo.ZAxis(latitude)*@CntZAxis) 
FROM Locations 
WHERE latitude IS NOT NULL 
ORDER BY ProxDistance ASC
GO

Editar:agregado (lo siento, me los perdí originalmente)

-- USER-DEFINED FUNCTIONS
-- XAxis #########################################
CREATE FUNCTION [dbo].[XAxis] (@lat float, @lon float)  
RETURNS float
AS  
BEGIN 
   RETURN COS(4 * (4 * atn2(1, 5) - atn2(1, 239)) / 180 * @lat) * COS(4 * (4 * atn2(1, 5) - atn2(1, 239)) / 180 * @lon) 
END

CREATE FUNCTION [dbo].[YAxis] (@lat float, @lon float)  
RETURNS float AS  
BEGIN 
RETURN COS(4 * (4 * atn2(1,5) - atn2(1,239)) / 180 * @lat) * SIN(4 * (4 * atn2(1,5) - atn2(1,239)) / 180 * @lon)
END

CREATE FUNCTION [dbo].[ZAxis] (@lat float)  
RETURNS float AS  
BEGIN 
RETURN SIN(4 * (4 * atn2(1,5) - atn2(1,239)) / 180 * @lat)
END