Add missing deg2rad conversion to angular distance calc

thor/observation_source.py

@@ -101,22 +101,34 @@ def _within_radius(
     radius: float,
 ) -> detections.PointSourceDetections:
     """
-    Return the detections within a given radius of a given ra and dec
+    Return the detections within a given radius of a given ra and dec.
+
+    ra, dec, and radius should be in degrees.
     """
-    sdlon = np.sin(detections.ra.to_numpy() - ra)
-    cdlon = np.cos(detections.ra.to_numpy() - ra)
-    slat1 = np.sin(dec)
-    slat2 = np.sin(detections.dec.to_numpy())
-    clat1 = np.cos(dec)
-    clat2 = np.cos(detections.dec.to_numpy())
+    det_ra = np.deg2rad(detections.ra.to_numpy())
+    det_dec = np.deg2rad(detections.dec.to_numpy())
+
+    center_ra = np.deg2rad(ra)
+    center_dec = np.deg2rad(dec)
+
+    dist_lon = det_ra - center_ra
+    sin_dist_lon = np.sin(dist_lon)
+    cos_dist_lon = np.cos(dist_lon)
+
+    sin_center_lat = np.sin(center_dec)
+    sin_det_lat = np.sin(det_dec)
+    cos_center_lat = np.cos(center_dec)
+    cos_det_lat = np.cos(det_dec)
 
-    num1 = clat2 * sdlon
-    num2 = clat1 * slat2 - slat1 * clat2 * cdlon
-    denominator = slat1 * slat2 + clat1 * clat2 * cdlon
+    num1 = cos_det_lat * sin_dist_lon
+    num2 = cos_center_lat * sin_det_lat - sin_center_lat * cos_det_lat * cos_dist_lon
+    denominator = (
+        sin_center_lat * sin_det_lat + cos_center_lat * cos_det_lat * cos_dist_lon
+    )
 
     distances = np.arctan2(np.hypot(num1, num2), denominator)
 
-    mask = distances <= radius
+    mask = distances <= np.deg2rad(radius)
     return detections.apply_mask(mask)
 
 

thor/tests/test_observation_source.py

@@ -17,7 +17,7 @@ def fixed_test_orbit():
         y=[0],
         z=[0],
         vx=[0],
-        vy=[2*np.pi/365.25],
+        vy=[2 * np.pi / 365.25],
         vz=[0],
         time=coordinates.Times.from_astropy(astropy.time.Time("2020-01-01T00:00:00")),
         origin=coordinates.Origin.from_kwargs(code=["SUN"]),
@@ -118,6 +118,6 @@ def test_fixed_radius_observation_source(fixed_test_orbit, fixed_observations):
     have = fos.gather_observations(fixed_test_orbit)
     assert len(have.exposures) == 5
     assert have.exposures == fixed_observations.exposures
-    assert len(have.detections) < len(fixed_observations.detections)
-    assert len(have.detections) > 0.75 * len(fixed_observations.detections)
-
+    # Should be about pi/4 fraction of the detections (0.785
+    assert len(have.detections) < 0.80 * len(fixed_observations.detections)
+    assert len(have.detections) > 0.76 * len(fixed_observations.detections)