use constant EARTH_RADIUS => 3956; use constant PI => 4 * atan2 1, 1; sub make_quick_reject { # Args: lat (radians), long (radians), dist my(\$lat,\$long,\$dist) = @_; my(\$distcos) = cos(\$dist / EARTH_RADIUS); my(\$maxlatdiff) = atan2(sqrt(1 - (\$distcos * \$distcos)), \$distcos); if (abs(\$lat) + \$maxlatdiff > PI/2) { # special case - near poles return sub { abs(\$lat - \$_[0]) > \$maxlatdiff; }; } my(\$longmult); \$longmult = 1 / cos(abs(\$lat) + \$maxlatdiff); my(\$maxlongdiff) = \$maxlatdiff * \$longmult; return sub { abs(\$lat - \$_[0]) > \$maxlatdiff or abs(\$long - \$_[1]) > \$maxlongdiff; } } sub make_quick_accept { # Args: lat (radians), long (radians), dist my(\$lat,\$long,\$dist) = @_; my(\$distcos) = cos(\$dist / EARTH_RADIUS); my(\$maxlatdiff) = atan2(sqrt(1 - (\$distcos * \$distcos)), \$distcos) / sqrt(2); my(\$longmult); if (abs(\$lat) < \$maxlatdiff ) { # special case - near equator \$longmult = 1; } else { \$longmult = 1 / cos(abs(\$lat) - \$maxlatdiff); } my(\$maxlongdiff) = \$maxlatdiff * \$longmult; return sub { abs(\$lat - \$_[0]) < \$maxlatdiff and abs(\$long - \$_[1]) < \$maxlongdiff; } }