Thanks Tom,

My good friend Dave Bernard gave me his distance function which uses teh great-circle distance calculations. Here's his:

> SELECT * FROM MAPWANTS WHERE
> DISTANCE(2,VAL(M.LAT),VAL(M.LONG),VAL(LAT),VAL(LONG),1)<=1
>
> The <= 1 is 1 mile.
>
> Here's his distance function:
> http://www.intellectiongroup.com/
> Dave Bernard
> bernard.dave@comcast.net
>
> *---------------------------------------------------------------------
> -*
> * Given a pair of lattitude and longitudes, return the approximate
> * distance in nautical miles between points A and B.
> *
> * P_FUNC  1 = DD:MM:SS format, 2 = DD.MMMM format
> * P_LT1   Lattitude of point A
> * P_LG1   Longitude of point A
> * P_LT1   Lattitude of point B
> * P_LG2   Longitude of point B
> * P_TYPE  0=Nautical Miles, 1=Statute Miles, 2=Kilometers
> *
> * All lattitude and longitude values are passed as strings in function 
> 1
> * and as numbers in function 2.
> *
> * Assumes N lattitude and W longitude
> *
> * Examples:
> * a = Distance(1, "33:57:00", "118:24:00", "40:38:00", "73:47:00", 2)
> * a = Distance(2, 33.95, 118.4, 40.6333, 73.7833, 2)
> *
> * Test with http://www.indo.com/distance
> *
> *---------------------------------------------------------------------
> -*
> function Distance
> parameters p_func, p_lt1, p_lg1, p_lt2, p_lg2, p_type
>
>    private p_func, p_lt1, p_lg1, p_lt2, p_lg2, p_svdec, p_pi,;
>            p1_degN, p1_minN, p1_secN, p1_degW, p1_minW, p1_secW,;
>            p2_degN, p2_minN, p2_secN, p2_degW, p2_minW, p2_secW,;
>            p_lat1, p_lon1, p_lat2, p_lon2, p_dist, p_type
>
>    p_svdec = set("decimals")
>
>    set decimals to 9
>    
>    if type("p_type") <> "N"
>       p_type = 0
>    endif
>    
>    if p_type > 2
>       p_type = 0
>    endif
>    
>    p_pi = pi()
>
>    do case
>       case p_func = 1
>            p_lt1 = p_lt1 + ":00"
>            p_lg1 = p_lg1 + ":00"
>            p_lt2 = p_lt2 + ":00"
>            p_lg2 = p_lg2 + ":00"
>            
>            p1_degN = val(left(p_lt1, at(":", p_lt1) - 1))
>            p_lt1   = substr(p_lt1, at(":", p_lt1) + 1)
>            p1_minN = val(left(p_lt1, at(":", p_lt1) - 1))
>            p_lt1   = substr(p_lt1, at(":", p_lt1) + 1)
>            p1_secN = val(p_lt1)
>
>            p1_degW = val(left(p_lg1, at(":", p_lg1) - 1))
>            p_lg1   = substr(p_lg1, at(":", p_lg1) + 1)
>            p1_minW = val(left(p_lg1, at(":", p_lg1) - 1))
>            p_lg1   = substr(p_lg1, at(":", p_lg1) + 1)
>            p1_secW = val(p_lg1)
>
>            p2_degN = val(left(p_lt2, at(":", p_lt2) - 1))
>            p_lt2   = substr(p_lt2, at(":", p_lt2) + 1)
>            p2_minN = val(left(p_lt2, at(":", p_lt2) - 1))
>            p_lt2   = substr(p_lt2, at(":", p_lt2) + 1)
>            p2_secN = val(p_lt2)
>
>            p2_degW = val(left(p_lg2, at(":", p_lg2) - 1))
>            p_lg2   = substr(p_lg2, at(":", p_lg2) + 1)
>            p2_minW = val(left(p_lg2, at(":", p_lg2) - 1))
>            p_lg2   = substr(p_lg2, at(":", p_lg2) + 1)
>            p2_secW = val(p_lg2)
>
>            p_lat1 = (p1_degN + ((p1_minN + (p1_secN / 60)) / 60)) * 
> p_pi / 180
> *                    p_pi * (p1_degN + (p1_minN / 60) + (p1_secN / 3600))
/
> 180
>            p_lon1 = (p1_degW + ((p1_minW + (p1_secW / 60)) / 60)) * 
> p_pi / 180
>            p_lat2 = (p2_degN + ((p2_minN + (p2_secN / 60)) / 60)) * 
> p_pi / 180
>            p_lon2 = (p2_degW + ((p2_minW + (p2_secW / 60)) / 60)) * 
> p_pi / 180
>       
>       case p_func = 2
>            p1_degN = p_lt1
>            p1_minN = 0
>            p1_secN = 0
>            p1_degW = p_lg1
>            p1_minW = 0
>            p1_secW = 0
>            p2_degN = p_lt2
>            p2_minN = 0
>            p2_secN = 0
>            p2_degW = p_lg2
>            p2_minW = 0
>            p2_secW = 0
>
>            p_lat1 = p1_degN * p_pi / 180
>            p_lon1 = p1_degW * p_pi / 180
>            p_lat2 = p2_degN * p_pi / 180
>            p_lon2 = p2_degW * p_pi / 180
>            
>       otherwise
>            set decimals to &p_svdec
>            return 0
>    endcase
>
> ***Method 1   
>    p_dist = acos((sin(p_lat1) * sin(p_lat2)) + (cos(p_lat1) * 
> cos(p_lat2)
*
> cos(p_lon1 - p_lon2)))
>    p_dist = int((p_dist * 180 * 60 / p_pi) +.5)		&& Distance
in
> Nautical Miles
>
> ***Method 2, with conversion to UTM (Mercator) coordinates
>    p_er = 6371.315				&& Average radius of the
> Earth
>
> ***Assume N lattitude and W longitude
> ***  For S lattitude: p_radlat1 = (p_pi / 2) + p_lat1
> ***  Others remain the same
>    p_radlat1 = (p_pi / 2) - p_lat1
>    p_radlon1 = (p_pi * 2) - p_lon1
>    p_radlat2 = (p_pi / 2) - p_lat2
>    p_radlon2 = (p_pi * 2) - p_lon2
>
> ***Spherical coordinates: x=r*cos(long)sin(lat), 
> y=r*sin(long)*cos(lat),
> z=r*cos(lat)
>    p_x1 = p_er * cos(p_radlon1) * sin(p_radlat1)
>    p_y1 = p_er * sin(p_radlon1) * sin(p_radlat1)
>    p_z1 = p_er * cos(p_radlat1)
>
>    p_x2 = p_er * cos(p_radlon2) * sin(p_radlat2)
>    p_y2 = p_er * sin(p_radlon2) * sin(p_radlat2)
>    p_z2 = p_er * cos(p_radlat2)
>
>    p_d = sqrt((p_x1 - p_x2)^2 + (p_y1 - p_y2)^2 + (p_z1 - p_z2)^2)
>    
> ***Side, side, side, law of cosines and arccos
>    p_theta = acos(((p_er^2 * 2) - (p_d^2)) / (p_er^2 * 2))
>
> ***These need to be rounded
>    p_dist2 = p_theta * p_er				&& Distance in
> Kilometers
>    p_dist3 = p_dist2 / 1.609344			&& Distance in Miles
>    p_dist4 = p_dist2 / 1.852			&& Distance in Nautical
> Miles
>    
>    set decimals to &p_svdec
>
>    do case
>       case p_type = 0
>            return p_dist
>       case p_type = 1
>            return p_dist3
>       case p_type = 2
>            return p_dist2
>    endcase
>
> return p_dist

>The radius of the earth is assumed to be 3963.189 miles. That makes the circumference of the earth 24,901.4297 miles using the equation of Radius * 2 * PI (3.14159). If you divide the circumference by 360 (number of degrees in a circle) you will get 69.1706 linear miles per latitude or longitude. However, the longitude value changes considerably as one moves from the equator toward the polar caps.
>
>The correct calculation needs to make use of spherical geometry and apply what is popularly known as the great-circle distance ( http://en.wikipedia.org/wiki/Great_circle_distance ) formula. The first step is to convert the latitude / longitude coordinates from degrees to radians via the DTOR() function. The rest of the equation makes uses of arc cosines, sines and cosines as follows:
>
>
>nLat1 = DTOR( m.nLat1 )
>nLon1 = DTOR( m.nLon1 )
>nLat2 = DTOR( m.nLat2 )
>nLon2 = DTOR( m.nLon2 )
>
>nRetVal = 3693.189 * ACOS( SIN( m.nLat1 ) * SIN( m.nLat2 ) + COS( m.nLat1 ) * ;
> COS( m.nLat2 ) * COS( m.nLon2 - m.nLon1 ) )
Thanks Tom. My good friend Dave had a similar solution using the

>
>Here is a code snippet of how I use this in an actual application (notice that .nLat1 and .nLon1 are pre-established just use the degrees to radians function and you are set):
>
>

> FUNCTION Distance( tnLat2, tnLon2 )
>    LOCAL nRetVal AS Number
>    WITH This
>     *-- Convert degrees to radians
>     .nLat2 = DTOR( m.tnLat2 )
>     .nLon2 = DTOR( m.tnLon2 )
>
>      *-- Radius of the earth is 3,963.189 miles
>      nRetVal = 3963.189 * ACOS( SIN( .nLat1 ) * SIN( .nLat2 ) + COS( .nLat1 ) * COS( .nLat2 ) * COS( .nLon2 - .nLon1 ) )
>    ENDWITH
>
>    RETURN m.nRetVal
>  ENDFUNC
>

>