Use Haversine formula in cone_filter (#624) (#636)

* Use astropy.coordinates.angular_separation in cone_filter (#624)

Replace manual spherical cosine formula with astropy's angular_separation,
converting all coordinates to radians and comparing the separation in radians.
Add an ASV benchmark for the cone_filter angular separation computation.

https://claude.ai/code/session_01D2WEdPXU3F8cfdmewEYXX8

* Rename benchmark method and fix uniform sphere sampling

Append 'benchmark' to the cone filter benchmark name, and use
arcsin(uniform(-1,1)) for declination to produce a uniform distribution
on the sphere instead of biasing toward the poles.

https://claude.ai/code/session_01D2WEdPXU3F8cfdmewEYXX8

* Rename AngularSeparationSuite to ConeFilterSuite

https://claude.ai/code/session_01D2WEdPXU3F8cfdmewEYXX8

* Use haversine formula instead of astropy angular_separation for speed

Replace the slow astropy.coordinates.angular_separation call with an
inline haversine formula using only numpy operations. All math remains
in radians, and the result is compared in radians.

https://claude.ai/code/session_01D2WEdPXU3F8cfdmewEYXX8

* Optimize haversine: compare sin² half-angles, avoid arcsin/sqrt

- Compare haversine value a <= sin²(radius/2) directly, eliminating
  the expensive arcsin and sqrt operations
- Use np.square instead of **2
- Multiply by 0.5 instead of dividing by 2
- Precompute cos(dec0) scalar

Valid for all separations up to 180 degrees since sin²(sep/2) is
monotonically increasing on [0, pi].

https://claude.ai/code/session_01D2WEdPXU3F8cfdmewEYXX8

* Fix pylint: initialize attributes in ConeFilterSuite.__init__

https://claude.ai/code/session_01D2WEdPXU3F8cfdmewEYXX8

* Reformat haversine equation for readability

https://claude.ai/code/session_01D2WEdPXU3F8cfdmewEYXX8

* Split haversine terms for readability, add comment

Break the haversine equation into named intermediate variables
(hav_dec, hav_ra) so the formatter doesn't mangle the layout.

https://claude.ai/code/session_01D2WEdPXU3F8cfdmewEYXX8

---------

Co-authored-by: Claude <noreply@anthropic.com>

Author: hombit

benchmarks/benchmarks.py

@@ -5,6 +5,7 @@
 
 from pathlib import Path
 
+import nested_pandas as npd
 import numpy as np
 
 import hats.pixel_math as hist
@@ -14,6 +15,7 @@
 from hats.pixel_math import HealpixPixel
 from hats.pixel_tree import align_trees
 from hats.pixel_tree.pixel_tree import PixelTree
+from hats.search.region_search import cone_filter
 
 BENCH_DATA_DIR = Path(__file__).parent / "data"
 
@@ -80,3 +82,34 @@ def time_small_cone_large_catalog():
     original_catalog = read_hats(BENCH_DATA_DIR / "large_catalog")
 
     original_catalog.filter_by_cone(315, -66.443, 1)
+
+
+class ConeFilterSuite:
+    """Benchmark the cone_filter angular separation computation."""
+
+    def __init__(self):
+        self.metadata = None
+        self.data_frame = None
+        self.ra = None
+        self.dec = None
+        self.radius_arcsec = None
+
+    def setup(self):
+        rng = np.random.default_rng(42)
+        n = 100_000
+        ra = rng.uniform(0, 360, n)
+        dec = np.degrees(np.arcsin(rng.uniform(-1, 1, n)))
+        self.metadata = TableProperties(
+            catalog_name="bench",
+            catalog_type="object",
+            total_rows=n,
+            ra_column="ra",
+            dec_column="dec",
+        )
+        self.data_frame = npd.NestedFrame({"ra": ra, "dec": dec})
+        self.ra = 180.0
+        self.dec = 0.0
+        self.radius_arcsec = 10.0 * 3600  # 10 degrees
+
+    def time_cone_filter_benchmark(self):
+        cone_filter(self.data_frame, self.ra, self.dec, self.radius_arcsec, self.metadata)

src/hats/search/region_search.py

@@ -92,13 +92,14 @@ def cone_filter(data_frame: npd.NestedFrame, ra, dec, radius_arcsec, metadata: T
     ra0 = np.radians(ra)
     dec0 = np.radians(dec)
 
-    cos_angle = np.sin(dec_rad) * np.sin(dec0) + np.cos(dec_rad) * np.cos(dec0) * np.cos(ra_rad - ra0)
-
-    # Clamp to valid range to avoid numerical issues
-    cos_separation = np.clip(cos_angle, -1.0, 1.0)
-
-    cos_radius = np.cos(np.radians(radius_arcsec / 3600))
-    data_frame = data_frame[cos_separation >= cos_radius]
+    cos_dec0 = np.cos(dec0)
+    # Haversine formula
+    hav_dec = np.square(np.sin((dec_rad - dec0) * 0.5))
+    hav_ra = np.square(np.sin((ra_rad - ra0) * 0.5))
+    a = hav_dec + cos_dec0 * np.cos(dec_rad) * hav_ra
+    radius_rad = np.radians(radius_arcsec / 3600)
+    threshold = np.square(np.sin(radius_rad * 0.5))
+    data_frame = data_frame[a <= threshold]
     return data_frame