Add detector architecture + isotropic detectors

src/pg_rad/physics/fluence.py

@@ -2,6 +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
 
@@ -12,6 +15,7 @@ def phi(
         branching_ratio: float,
         mu_mass_air: float,
         air_density: float,
+        eff: float
         ) -> float:
     """Compute the contribution of a single point source to the
     primary photon fluence rate phi at position (x,y,z).
@@ -34,6 +38,7 @@ def phi(
 
     phi_r = (
         activity
+        * eff
         * branching_ratio
         * np.exp(-mu_air * r)
         / (4 * np.pi * r**2)
@@ -42,12 +47,20 @@ def phi(
     return phi_r
 
 
-def calculate_fluence_at(landscape: "Landscape", pos: np.ndarray, scaling=1E6):
+def calculate_fluence_at(
+    landscape: "Landscape",
+    pos: np.ndarray,
+    detector: IsotropicDetector | AngularDetector,
+    tangent_vectors: np.ndarray,
+    scaling=1E6
+):
     """Compute fluence at an arbitrary position in the landscape.
 
     Args:
         landscape (Landscape): The landscape to compute.
         pos (np.ndarray): (N, 3) array of positions.
+        detector (IsotropicDetector | AngularDetector):
+            Detector object, needed to compute correct efficiency.
 
     Returns:
         total_phi (np.ndarray): (N,) array of fluences.
@@ -56,15 +69,30 @@ def calculate_fluence_at(landscape: "Landscape", pos: np.ndarray, scaling=1E6):
     total_phi = np.zeros(pos.shape[0])
 
     for source in landscape.point_sources:
-        r = np.linalg.norm(pos - np.array(source.pos), axis=1)
+        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)
+        else:
+            eff = detector.get_efficiency(energy=source.isotope.E)
+
         phi_source = phi(
             r=r,
             activity=source.activity * scaling,
             branching_ratio=source.isotope.b,
             mu_mass_air=source.isotope.mu_mass_air,
-            air_density=landscape.air_density
+            air_density=landscape.air_density,
+            eff=eff
         )
 
         total_phi += phi_source
@@ -74,6 +102,7 @@ def calculate_fluence_at(landscape: "Landscape", pos: np.ndarray, scaling=1E6):
 
 def calculate_fluence_along_path(
     landscape: "Landscape",
+    detector: "IsotropicDetector | AngularDetector",
     points_per_segment: int = 10
 ) -> Tuple[np.ndarray, np.ndarray]:
     path = landscape.path
@@ -98,6 +127,17 @@ def calculate_fluence_along_path(
     z = np.full(xnew.shape, path.z)
     full_positions = np.c_[xnew, ynew, z]
 
-    phi_result = calculate_fluence_at(landscape, full_positions)
+    # 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)
+
+    phi_result = calculate_fluence_at(
+        landscape,
+        full_positions,
+        detector,
+        tangent_vectors)
 
     return s, phi_result