>*circumference is 360 degrees times 60 minutes/degree times 1.852 km/minute = 40003.2 km. >*The implied radius is the circumference divided by 2 pi: R = 6367 km = 3956 mi >lon1= 1698101.25 &&1713881.3750 > lat1= 813831.69 && 833601.8125 >lon2 = 1698215.13 &&1697432.8750 >lat2 = 813458.25 && 810161.0625 > R = 3956 && 3963 && 3956 > dlon = lon2 - lon1 > dlat = lat2 - lat1 > a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2 > ? a > c = 2 * atn2(sqrt(a), sqrt(1-a)) > ? c > d = R * c > ?[Example1: ] > ??d > > >*--is point within a radius >radius = 500 &&(what would 500 ft be?) >centerx = 1698101.25 &&1713881.3750 >centery = 813831.69 && 833601.8125 >x = 1698215.13 &&1697432.8750 >y = 813458.25 && 810161.0625 > >?[Example2: ] >?? (x-centerx)^2 + (y - centery)^2 >?[ Inside radius: ] >?? (x-centerx)^2 + (y - centery)^2 < radius^2 > > > >*--Get distance between two points >#DEFINE KMRAD 6377 >#DEFINE MIRAD 3956 && 3963 >#DEFINE RADIANS DTOR > > >tlKM = .F. > >lon1= 1698101.25 &&1713881.3750 >lat1= 813831.69 && 833601.8125 >lon2 = 1698215.13 &&1697432.8750 >lat2 = 813458.25 && 810161.0625 > > lnDistance=IIF(tlKM, KMRAD, MIRAD)*; > ACOS(COS(DTOR(90-lat1))*COS(DTOR(90-lat2))+; > SIN(DTOR(90-lat1))*SIN(DTOR(90-lat2))*; > COS(DTOR(lon1 - lon2))) > >?[Example3: ] >??lnDistance && *5280 && 5280 feet in a mile > > > >*--Select all within the radius &&Courtesy Carl Karsten >#DEFINE cdbl >#DEFINE _PI 3.1415926535893 >#DEFINE _StatuteMile 5280 &&' 5280 feet = 1 statute mile >#DEFINE _NauticalMile 6076.11549 &&' 6076.11549 feet = 1 nautical mile >#DEFINE _Seconds 60 && ' 60 seconds = 1 nautical mile >lon1= 1698101.25 &&1713881.3750 >lat1= 813831.69 && 833601.8125 >lon2 = 1698215.13 &&1697432.8750 >lat2 = 813458.25 && 810161.0625 > >? InRadius(lat1,lon1,500,lat2,lon2) >return > >FUNCTION InRadius( tlat, tlong, tnRadius, tlat2, tlong2 ) > LOCAL lnResults, lnSeconds, ; > lnLatRange, lnLonRange,; > lnLowLat, lnLowLon, ; > lnHighLat, lnHighLon > > *' 1 degree of latitude = 60 nautical miles or 1 minute = 1 nautical mile > lnLatRange = tnRadius / ((_NauticalMile / _StatuteMile) * _Seconds) > > * Longitude is a bit more complicated; > * 1 degree of longitude at the equator = 60 nautical miles, > * At the poles 1 degree = 0 nautical miles. > lnLonRange = tnRadius / (((COS(cdbl(tlat * _PI / 180)) * ; > _NauticalMile) / _StatuteMile) * _Seconds) > > lnLowLat = tlat - lnLatRange > lnHighLat = tlat + lnLatRange > lnLowLon = tlong - lnLonRange > lnHighLon = tlong+ lnLonRange > >?[Example4: ] >?? tlat2 <= lnHighLat; > AND tlat2 >= lnLowLat; > AND tlong2 >= lnLowLon; > AND tlong2 <= lnHighLon >?[Example5: ] >?? getdistance2(lon1,lat1,lon2,lat2) > >RETURN > > >FUNCTION getdistance2 &&Courtesy Carl Karsten >LPARAMETERS lon1,lat1,lon2,lat2 > >*--Get distance between two points >#DEFINE KMRAD 6377 >#DEFINE MIRAD 3956 && 3963 >#DEFINE RADIANS DTOR > > >tlKM = .F. > >lon1= 1698101.25 &&1713881.3750 >lat1= 813831.69 && 833601.8125 >lon2 = 1698215.13 &&1697432.8750 >lat2 = 813458.25 && 810161.0625 > > lnDistance=IIF(tlKM, KMRAD, MIRAD)*; > ACOS(COS(DTOR(90-lat1))*COS(DTOR(90-lat2))+; > SIN(DTOR(90-lat1))*SIN(DTOR(90-lat2))*; > COS(DTOR(lon1 - lon2))) > >?[Example6: ] >??lnDistance > >return lndistance > >>
* FUNCTION GreatCircleDistance LPARAMETERS ; tnLat1 ; , tnLon1 ; , tnLat2 ; , tnLon2 LOCAL ; lnDeltaLat ; , lnDeltaLon ; , lnA * All 4 parameters are in degrees * Return value is in kilometers * Java implementation from Wikipedia: *!* double haversine(double lat1, double lng1, double lat2, double lng2) { *!* double dLat = Math.toRadians(lat2 - lat1); *!* double dLon = Math.toRadians(lng2 - lng1); *!* double a = Math.sin(dLat / 2) * Math.sin(dLat / 2) + Math.cos(Math.toRadians(lat1)) * *!* Math.cos(Math.toRadians(lat2)) * Math.sin(dLon / 2) * Math.sin(dLon / 2); *!* return 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1 - a)) * 6371; // 6371 = average earth radius * Calculate in radians and kilometers * There are 2 * PI() radians in a circle, one radian is 180 / PI() degrees lnDeltaLat = ( tnLat2 - tnLat1 ) * ( PI( ) / 180 ) lnDeltaLon = ( tnLon2 - tnLon1 ) * ( PI( ) / 180 ) lnA = ( SIN( lnDeltaLat / 2 ) ) ^ 2 ; + COS( tnLat1 * PI( ) / 180 ) * COS( tnLat2 * PI( ) / 180 ) * ( SIN( lnDeltaLon / 2 ) ) ^ 2 RETURN 2 * Atan2( SQRT ( lnA ), SQRT( 1 - lnA ) ) * 6371 ***************************** FUNCTION Atan2 LPARAMETERS ; tnY ; , tnX LOCAL ; lnRetVal * Parameters must be passed in radians * From Wikipedia definition of ATan2() function: IF tnX > 0 lnRetVal = ATAN( tnY / tnX ) ELSE IF tnX < 0 IF tnY < 0 lnRetVal = ATAN( tnY / tnX ) - PI( ) ELSE lnRetVal = ATAN( tnY / tnX ) + PI( ) ENDIF ELSE IF tnY > 0 lnRetVal = PI( ) / 2 ELSE IF tnY < 0 lnRetVal = - ( PI( ) / 2 ) ELSE * Return value where both X and Y are zero is undefined, we'll return 0: lnRetVal = 0 ENDIF ENDIF ENDIF ENDIF RETURN lnRetVal *****************************I tested this exactly once, with a known distance near my home. I got exact lat/lons from Google Maps using a cool site http://www.getlatlon.com/ . I moved the crosshairs to two locations that I know are 350m apart on a curved track, and copy/pasted the exact lat/lons the site provided. The functions above returned 0.29km ( 290m ) distance between them, which is quite believable.