Avoid violation of acos' domain in map_projection

* There are pathological cases (e.g., setpoint very close to reference
  for certain reference latitudes), where numerical errors lead to a sum
  larger than 1.0 passed to acos, resulting in NaN values.

This should fix issue #2813
This commit is contained in:
Peter Duerr
2015-12-04 14:14:31 +09:00
parent 2a1b1fe11d
commit c070d326e0
+7 -1
View File
@@ -148,7 +148,13 @@ __EXPORT int map_projection_project(const struct map_projection_reference_s *ref
double cos_lat = cos(lat_rad);
double cos_d_lon = cos(lon_rad - ref->lon_rad);
double c = acos(ref->sin_lat * sin_lat + ref->cos_lat * cos_lat * cos_d_lon);
double arg = ref->sin_lat * sin_lat + ref->cos_lat * cos_lat * cos_d_lon;
if (arg > 1.0) {
arg = 1.0;
} else if (arg < -1.0) {
arg = -1.0;
}
double c = acos(arg);
double k = (fabs(c) < DBL_EPSILON) ? 1.0 : (c / sin(c));
*x = k * (ref->cos_lat * sin_lat - ref->sin_lat * cos_lat * cos_d_lon) * CONSTANTS_RADIUS_OF_EARTH;