Use the Vincenty formula rather than the Haversine formula for earth_distance.

Vincenty is much better behaved with small distances. This is a translation of
http://www.movable-type.co.uk/scripts/latlong-vincenty.html

While I'm here, put a couple of WGS84 constants into a header.

1 files changed, 2 insertions, 3 deletions

diff --git a/geoid.c b/geoid.c
index fa79104d..5db3ca68 100644
--- a/geoid.c
+++ b/geoid.c
@@ -90,9 +90,8 @@ void ecef_to_wgs84fix(struct gps_data_t *gpsdata,
 /* fill in WGS84 position/velocity fields from ECEF coordinates */
 {
     double lambda,phi,p,theta,n,h,vnorth,veast,heading;
-    const double a = 6378137;			/* equatorial radius */
-    const double f = 1 / 298.257223563;		/* flattening */
-    const double b = a * (1 - f);		/* polar radius */
+    const double a = WGS84A;			/* equatorial radius */
+    const double b = WGS84B;		/* polar radius */
     const double e2 = (a*a - b*b) / (a*a);
     const double e_2 = (a*a - b*b) / (b*b);