update fluence and simulation to work with detectors and correct count rates

src/pg_rad/physics/fluence.py

@@ -2,11 +2,9 @@ from typing import Tuple, TYPE_CHECKING
 
 import numpy as np
 
-from pg_rad.detector.detectors import IsotropicDetector, AngularDetector
-
-
 if TYPE_CHECKING:
     from pg_rad.landscape.landscape import Landscape
+    from pg_rad.detector.detector import Detector
 
 
 def phi(
@@ -33,16 +31,14 @@ def phi(
     """
 
     # Linear photon attenuation coefficient in m^-1.
-    mu_mass_air *= 0.1
-    mu_air = mu_mass_air * air_density
+    mu_air = 0.1 * mu_mass_air * air_density
 
     phi_r = (
         activity
         * eff
         * branching_ratio
         * np.exp(-mu_air * r)
-        / (4 * np.pi * r**2)
-    )
+    ) / (4 * np.pi * r**2)
 
     return phi_r
 
@@ -50,8 +46,7 @@ def phi(
 def calculate_count_rate_per_second(
     landscape: "Landscape",
     pos: np.ndarray,
-    detector: IsotropicDetector | AngularDetector,
-    tangent_vectors: np.ndarray,
+    detector: "Detector",
     scaling=1E6
 ):
     """Compute count rate in s^-1 m^-2 at a position in the landscape.
@@ -59,7 +54,7 @@ def calculate_count_rate_per_second(
     Args:
         landscape (Landscape): The landscape to compute.
         pos (np.ndarray): (N, 3) array of positions.
-        detector (IsotropicDetector | AngularDetector):
+        detector (Detector):
             Detector object, needed to compute correct efficiency.
 
     Returns:
@@ -69,22 +64,29 @@ def calculate_count_rate_per_second(
     total_phi = np.zeros(pos.shape[0])
 
     for source in landscape.point_sources:
+        # See Bukartas (2021) page 25 for incidence angle math
         source_to_detector = pos - np.array(source.pos)
         r = np.linalg.norm(source_to_detector, axis=1)
         r = np.maximum(r, 1E-3)  # enforce minimum distance of 1cm
 
-        if isinstance(detector, AngularDetector):
-            cos_theta = (
-                np.sum(tangent_vectors * source_to_detector, axis=1) / (
-                    np.linalg.norm(source_to_detector, axis=1) *
-                    np.linalg.norm(tangent_vectors, axis=1)
-                )
-            )
-            cos_theta = np.clip(cos_theta, -1, 1)
-            theta = np.arccos(cos_theta)
-            eff = detector.get_efficiency(theta, energy=source.isotope.E)
+        if not detector.is_isotropic:
+            v = np.zeros_like(pos)
+            v[1:] = pos[1:] - pos[:-1]
+            v[0] = v[1]  # handle first point
+            vx, vy = v[:, 0], v[:, 1]
+
+            r_vec = pos - np.array(source.pos)
+            rx, ry = r_vec[:, 0], r_vec[:, 1]
+
+            theta = np.arctan2(vy, vx) - np.arctan2(ry, rx)
+
+            # normalise to [-pi, pi] and convert to degrees
+            theta = (theta + np.pi) % (2 * np.pi) - np.pi
+            theta_deg = np.degrees(theta)
+
+            eff = detector.get_efficiency(source.isotope.E, theta_deg)
         else:
-            eff = detector.get_efficiency(energy=source.isotope.E)
+            eff = detector.get_efficiency(source.isotope.E)
 
         phi_source = phi(
             r=r,
@@ -102,15 +104,15 @@ def calculate_count_rate_per_second(
 
 def calculate_counts_along_path(
     landscape: "Landscape",
-    detector: "IsotropicDetector | AngularDetector",
-    acquisition_time: int,
+    detector: "Detector",
+    velocity: float,
     points_per_segment: int = 10,
 ) -> Tuple[np.ndarray, np.ndarray]:
     """Compute the counts recorded in each acquisition period in the landscape.
 
     Args:
         landscape (Landscape): _description_
-        detector (IsotropicDetector | AngularDetector): _description_
+        detector (Detector): _description_
         points_per_segment (int, optional): _description_. Defaults to 100.
 
     Returns:
@@ -123,7 +125,6 @@ def calculate_counts_along_path(
     num_segments = len(path.segments)
 
     segment_lengths = np.array([seg.length for seg in path.segments])
-    ds = segment_lengths[0]
     original_distances = np.zeros(num_points)
     original_distances[1:] = np.cumsum(segment_lengths)
 
@@ -140,26 +141,17 @@ def calculate_counts_along_path(
     if path.opposite_direction:
         full_positions = np.flip(full_positions, axis=0)
 
-    # to compute the angle between sources and the direction of travel, we
-    # compute tangent vectors along the path.
-    dx_ds = np.gradient(xnew, s)
-    dy_ds = np.gradient(ynew, s)
-    tangent_vectors = np.c_[dx_ds, dy_ds, np.zeros_like(dx_ds)]
-    tangent_vectors /= np.linalg.norm(tangent_vectors, axis=1, keepdims=True)
-
-    count_rate = calculate_count_rate_per_second(
-        landscape, full_positions, detector, tangent_vectors
+    # [counts/s]
+    cps = calculate_count_rate_per_second(
+        landscape, full_positions, detector
     )
 
-    count_rate *= (acquisition_time / points_per_segment)
+    # reshape so each segment is on a row
+    cps_per_seg = cps.reshape(num_segments, points_per_segment)
 
-    count_rate_segs = count_rate.reshape(num_segments, points_per_segment)
-    integrated = np.trapezoid(
-        count_rate_segs,
-        dx=ds/points_per_segment,
-        axis=1
-    )
+    du = s[1] - s[0]
+    integrated_counts = np.trapezoid(cps_per_seg, dx=du, axis=1) / velocity
+    int_counts_result = np.zeros(num_points)
+    int_counts_result[1:] = integrated_counts
 
-    result = np.zeros(num_points)
-    result[1:] = integrated
-    return original_distances, result
+    return original_distances, s, cps, int_counts_result