summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorGravatar Doug Junkins <junkins@foghead.com>2019-04-30 06:53:31 -0700
committerGravatar bstoeger <32835590+bstoeger@users.noreply.github.com>2019-04-30 23:32:50 +0200
commit14a763a6a0148a42fad06717bd2532181a54c626 (patch)
treed8b71f0c9a2d33d40d60aeccb44e12111abcb2e0
parent154792ffd1c217e5c1c9774b5593a753ae553601 (diff)
downloadsubsurface-14a763a6a0148a42fad06717bd2532181a54c626.tar.gz
Bugfix for algorithm in get_distance()
Fixed bug in the Haversine function in get_distance() based on algorithm at https://www.movable-type.co.uk/scripts/latlong.html and added bounds to the 'a' term to avoid floating point errors for antipodal points. Signed-off-by: Doug Junkins <junkins@foghead.com>
-rw-r--r--CHANGELOG.md1
-rw-r--r--core/divesite.c7
2 files changed, 6 insertions, 2 deletions
diff --git a/CHANGELOG.md b/CHANGELOG.md
index cce954eb1..6b3930ac7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,4 @@
+- Core: fix bug in get_distance() to correctly compute spherical distance
- Desktop: For videos, add save data export as subtitle file
- Desktop: make dive sites 1st class citizens with their own dive site table
- Desktop: only show dives at the dive sites selected in dive site table
diff --git a/core/divesite.c b/core/divesite.c
index b5ce99a8d..45601537c 100644
--- a/core/divesite.c
+++ b/core/divesite.c
@@ -77,15 +77,18 @@ struct dive_site *get_dive_site_by_gps_and_name(char *name, const location_t *lo
// Calculate the distance in meters between two coordinates.
unsigned int get_distance(const location_t *loc1, const location_t *loc2)
{
+ double lat1_r = udeg_to_radians(loc1->lat.udeg);
double lat2_r = udeg_to_radians(loc2->lat.udeg);
double lat_d_r = udeg_to_radians(loc2->lat.udeg - loc1->lat.udeg);
double lon_d_r = udeg_to_radians(loc2->lon.udeg - loc1->lon.udeg);
double a = sin(lat_d_r/2) * sin(lat_d_r/2) +
- cos(lat2_r) * cos(lat2_r) * sin(lon_d_r/2) * sin(lon_d_r/2);
+ cos(lat1_r) * cos(lat2_r) * sin(lon_d_r/2) * sin(lon_d_r/2);
+ if (a < 0.0) a = 0.0;
+ if (a > 1.0) a = 1.0;
double c = 2 * atan2(sqrt(a), sqrt(1.0 - a));
- // Earth radious in metres
+ // Earth radius in metres
return lrint(6371000 * c);
}