pascal/lazarus-ccr
1
0
Fork 0
Code Issues Releases Activity

lazmapviewer: Improved numerical stability of geo distance calculation.

Browse Source 
git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@6867 8e941d3f-bd1b-0410-a28a-d453659cc2b4
This commit is contained in:
wp_xxyyzz
2019-04-25 17:43:20 +00:00
parent 23a253de0b
commit fb276886e4
1 changed files with 13 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

19
components/lazmapviewer/source/mvengine.pas
View File

@@ -1131,16 +1131,17 @@ end;
function CalcGeoDistance(Lat1, Lon1, Lat2, Lon2: double; function CalcGeoDistance(Lat1, Lon1, Lat2, Lon2: double;
AUnits: TDistanceUnits = duKilometers): double; AUnits: TDistanceUnits = duKilometers): double;
const const
raduis_km: double = 6371.0; // generalized radius of earth EPS = 1E-12;
var var
d_radians: double; // distance in radians d_radians: double; // distance in radians
lat1r, lon1r, lat2r, lon2r: double; lat1r, lon1r, lat2r, lon2r: double;
arg: Double;
begin begin
// Validate // Validate
if (Lat1 < -90.0) or (Lat1 > 90.0) then exit(-1); if (Lat1 < -90.0) or (Lat1 > 90.0) then exit(NaN);
if (Lon1 < -180.0) or (Lon1 > 180.0) then exit(-1); // if (Lon1 < -180.0) or (Lon1 > 180.0) then exit(NaN);
if (Lat2 < -90.0) or (Lat2 > 90.0) then exit(-1); if (Lat2 < -90.0) or (Lat2 > 90.0) then exit(NaN);
if (Lon2 < -180.0) or (Lon2 > 180.0) then exit(-1); // if (Lon2 < -180.0) or (Lon2 > 180.0) then exit(NaN);
// Turn lat and lon into radian measures // Turn lat and lon into radian measures
lat1r := (PI / 180.0) * Lat1; lat1r := (PI / 180.0) * Lat1;
@@ -1149,7 +1150,13 @@ begin
lon2r := (PI / 180.0) * Lon2; lon2r := (PI / 180.0) * Lon2;
// calc // calc
d_radians := arccos(Sin(lat1r) * sin(lat2r) + cos(lat1r) * cos(lat2r) * cos(lon1r - lon2r)); arg := sin(lat1r) * sin(lat2r) + cos(lat1r) * cos(lat2r) * cos(lon1r - lon2r);
if (arg < -1) or (arg > +1) then
exit(NaN);
if SameValue(abs(Lon1-Lon2), 360, EPS) and SameValue(abs(arg), 1.0, EPS) then
d_radians := PI * 2.0
else
d_radians := arccos(arg);
Result := EARTH_RADIUS * d_radians; Result := EARTH_RADIUS * d_radians;
case AUnits of case AUnits of