diff options
| author | Eric S. Raymond <esr@thyrsus.com> | 2012-04-20 05:14:52 -0400 |
|---|---|---|
| committer | Eric S. Raymond <esr@thyrsus.com> | 2012-04-20 05:14:52 -0400 |
| commit | 4dc7b5468e669dce490951f523d9678e8616b45f (patch) | |
| tree | ee3661fe984fffa6d864ceabe71fdd686e49caa9 /geoid.c | |
| parent | a02226a6461a4bdb245f6c002c0f0e90e6d0b1fb (diff) | |
| download | gpsd-4dc7b5468e669dce490951f523d9678e8616b45f.tar.gz | |
Avoid a possible nunerical instability in interpolating geoid separation.
All regression tesrs pass.
Diffstat (limited to 'geoid.c')
| -rw-r--r-- | geoid.c | 12 |
1 files changed, 7 insertions, 5 deletions
@@ -19,12 +19,14 @@ static double bilinear(double x1, double y1, double x2, double y2, double x, { double delta; - if (y1 == y2 && x1 == x2) +#define EQ(a, b) (fabs((a) - (b)) < 0.001) + if (EQ(y1, y2) && EQ(x1, x2)) return (z11); - if (y1 == y2 && x1 != x2) + if (EQ(y1, y2) && !EQ(x1, x2)) return (z22 * (x - x1) + z11 * (x2 - x)) / (x2 - x1); - if (x1 == x2 && y1 != y2) + if (EQ(x1, x2) && !EQ(y1, y2)) return (z22 * (y - y1) + z11 * (y2 - y)) / (y2 - y1); +#undef EQ delta = (y2 - y1) * (x2 - x1); @@ -82,8 +84,8 @@ double wgs84_separation(double lat, double lon) ilat2 = (ilat < GEOID_ROW - 1) ? ilat + 1 : ilat; ilon2 = (ilon < GEOID_COL - 1) ? ilon + 1 : ilon; - return bilinear(ilon1 * 10. - 180., ilat1 * 10. - 90., - ilon2 * 10. - 180., ilat2 * 10. - 90., + return bilinear(ilon1 * 10.0 - 180.0, ilat1 * 10.0 - 90.0, + ilon2 * 10.0 - 180.0, ilat2 * 10.0 - 90.0, lon, lat, (double)geoid_delta[ilon1 + ilat1 * GEOID_COL], (double)geoid_delta[ilon2 + ilat1 * GEOID_COL], |