Geo: fix computation of bounding box.

A bug was reported in the context in issue #3631. The root cause of the
bug was that certain neighbor boxes were zeroed after the "inside the
bounding box or not" check, simply because the bounding box computation
function was wrong.

A few debugging infos where enhanced and moved in other parts of the
code. A check to avoid steps=0 was added, but is unrelated to this
issue and I did not verified it was an actual bug in practice.

1 files changed, 9 insertions, 43 deletions

diff --git a/src/geohash_helper.c b/src/geohash_helper.c
index 139bcea11..77d8ab392 100644
--- a/src/geohash_helper.c
+++ b/src/geohash_helper.c
@@ -82,30 +82,18 @@ uint8_t geohashEstimateStepsByRadius(double range_meters, double lat) {
     return step;
 }
 
+/* Return the bounding box of the search area centered at latitude,longitude
+ * having a radius of radius_meter. bounds[0] - bounds[2] is the minimum
+ * and maxium longitude, while bounds[1] - bounds[3] is the minimum and
+ * maximum latitude. */
 int geohashBoundingBox(double longitude, double latitude, double radius_meters,
                        double *bounds) {
     if (!bounds) return 0;
 
-    double lonr, latr;
-    lonr = deg_rad(longitude);
-    latr = deg_rad(latitude);
-
-    if (radius_meters > EARTH_RADIUS_IN_METERS)
-        radius_meters = EARTH_RADIUS_IN_METERS;
-    double distance = radius_meters / EARTH_RADIUS_IN_METERS;
-    double min_latitude = latr - distance;
-    double max_latitude = latr + distance;
-
-    /* Note: we're being lazy and not accounting for coordinates near poles */
-    double min_longitude, max_longitude;
-    double difference_longitude = asin(sin(distance) / cos(latr));
-    min_longitude = lonr - difference_longitude;
-    max_longitude = lonr + difference_longitude;
-
-    bounds[0] = rad_deg(min_longitude);
-    bounds[1] = rad_deg(min_latitude);
-    bounds[2] = rad_deg(max_longitude);
-    bounds[3] = rad_deg(max_latitude);
+    bounds[0] = longitude - rad_deg(radius_meters/EARTH_RADIUS_IN_METERS/cos(deg_rad(latitude)));
+    bounds[2] = longitude + rad_deg(radius_meters/EARTH_RADIUS_IN_METERS/cos(deg_rad(latitude)));
+    bounds[1] = latitude - rad_deg(radius_meters/EARTH_RADIUS_IN_METERS);
+    bounds[3] = latitude + rad_deg(radius_meters/EARTH_RADIUS_IN_METERS);
     return 1;
 }
 
@@ -158,35 +146,13 @@ GeoHashRadius geohashGetAreasByRadius(double longitude, double latitude, double
             < radius_meters) decrease_step = 1;
     }
 
-    if (decrease_step) {
+    if (steps > 1 && decrease_step) {
         steps--;
         geohashEncode(&long_range,&lat_range,longitude,latitude,steps,&hash);
         geohashNeighbors(&hash,&neighbors);
         geohashDecode(long_range,lat_range,hash,&area);
     }
 
-    /* Example debug info. This turns to be very useful every time there is
-     * to investigate radius search potential bugs. So better to leave it
-     * here. */
-    if (0) {
-        GeoHashArea myarea = {{0}};
-        geohashDecode(long_range, lat_range, neighbors.west, &myarea);
-
-        /* Dump West. */
-        D("Neighbors");
-        D("area.longitude.min: %f\n", myarea.longitude.min);
-        D("area.longitude.max: %f\n", myarea.longitude.max);
-        D("area.latitude.min: %f\n", myarea.latitude.min);
-        D("area.latitude.max: %f\n", myarea.latitude.max);
-
-        /* Dump center square. */
-        D("Area");
-        D("area.longitude.min: %f\n", area.longitude.min);
-        D("area.longitude.max: %f\n", area.longitude.max);
-        D("area.latitude.min: %f\n", area.latitude.min);
-        D("area.latitude.max: %f\n", area.latitude.max);
-    }
-
     /* Exclude the search areas that are useless. */
     if (area.latitude.min < min_lat) {
         GZERO(neighbors.south);