Initial raytracing model implementation

src/ray_tracer/__init__.py

+import numpy as np
+
+# Constants and parameters
+pi = np.pi
+e = 1.60217662e-19  # Elementary charge
+atomic_number = 21.5
+
+# User-defined parameters
+beam_current = 3e-3  # Beam current in Amperes
+dwell_time = 0.4e-6  # Dwell time in seconds
+transimpedance = 1e5  # Transimpedance in Volts per Ampere
+
+# Geometry
+build_surface_width = 80e-3  # Build surface width in meters
+num_scan_lines = 1500  # Number of scan lines
+pixel_size = 53.3e-6  # Pixel size in meters
+detector_height = 272e-3  # Detector height in meters
+detector_spacing = 127e-3  # Detector spacing in meters
+
+# Scattering function parameters
+theta = np.deg2rad(0)
+k = np.deg2rad(theta / 13)  # Scattering lobe width parameter
+
+# Define domain size
+dx = 5e-6  # Cell size in meters
+dy = dx
+dz = dx
+x_min, x_max = -250e-6, 250e-6
+y_min, y_max = -250e-6, 250e-6
+z_min, z_max = 0, 20e-6
+x_num = int((x_max - x_min) / dx) + 1
+y_num = int((y_max - y_min) / dy) + 1
+z_num = int((z_max - z_min) / dz) + 1
+
+# Initialize the domain array
+domain = np.zeros((z_num, y_num, x_num))
+domain[:, :, 0:2] = 1  # Solid cells
+domain[:, :, 2] = 0.5  # Interface layer
+
+# Calculate gradients to determine surface normals
+fill_level = domain
+grad_x = np.gradient(fill_level, axis=2)
+grad_y = np.gradient(fill_level, axis=1)
+grad_z = np.gradient(fill_level, axis=0)
+n_s = np.stack((grad_x, grad_y, grad_z), axis=-1)
+norm = np.linalg.norm(n_s, axis=-1)
+
+interface_mask = (fill_level > 0) & (fill_level < 1)
+nonzero_mask = norm > 0
+valid_mask = interface_mask & nonzero_mask
+n_s[valid_mask] /= norm[valid_mask][..., np.newaxis]
+n_s[fill_level == 1] = np.array([0, 0, -1])
+n_s[fill_level == 0] = np.array([0, 0, 1])
+
+# Beam profile calculation
+sigma = 300e-6 / (2 * dx)
+x = np.arange(-30, 31)
+y = np.arange(-30, 31)
+x, y = np.meshgrid(x, y)
+beam_profile = np.exp(-(x ** 2 + y ** 2) / (2 * sigma ** 2))
+beam_profile /= np.sum(beam_profile)
+
+# Extend beam_profile to 3D to match the domain shape
+beam_profile_3d = np.tile(beam_profile[:, :, np.newaxis], (1, 1, z_num))
+
+# Extend beam profile to 3D to match domain shape
+'''beam_profile_3d = np.zeros((z_num, y_num, x_num))
+beam_profile_3d[z_num // 2 - beam_profile.shape[0] // 2:z_num // 2 + beam_profile.shape[0] // 2,
+y_num // 2 - beam_profile.shape[1] // 2:y_num // 2 + beam_profile.shape[1] // 2,
+x_num // 2] = beam_profile'''
+
+# Beam direction
+beam_dir = np.array([np.sin(theta), 0, -np.cos(theta)])
+
+# Define detector
+detector_distance = int(120e-3 / dx)
+detector_height_cells = int(300e-3 / dx)
+detector_width = 10
+detector_vertices = np.array([
+    [detector_distance, -detector_width / 2, detector_height_cells],
+    [detector_distance, detector_width / 2, detector_height_cells],
+    [detector_distance, 0, 0]
+])
+
+# Initialize arrays
+beam_interactions = np.zeros_like(domain)
+scattered_electrons = np.zeros((num_scan_lines, int(build_surface_width / pixel_size)))
+
+
+def compute_reflection_vector(beam_dir, surface_normal):
+    beam_dir = beam_dir / np.linalg.norm(beam_dir)
+    surface_normal = surface_normal / np.linalg.norm(surface_normal)
+    p_reflected = beam_dir - 2 * np.dot(beam_dir, surface_normal) * surface_normal
+    return p_reflected
+
+
+def bse_coefficient_0(A):
+    return -0.0254 + 0.016 * A - 1.86e-4 * A ** 2 + 8.3e-7 * A ** 3
+
+
+def bse_coefficient(theta, A):
+    B = 0.89
+    if theta == 0:
+        return bse_coefficient_0(A)
+    return B * (bse_coefficient_0(A) / B) ** np.cos(theta)
+
+
+def scattering_function_correction(theta):
+    return 1.01 - 4.06 * theta + 1.36e-5 * theta ** 2 - 1.71 * theta ** 3
+
+
+def scattering_function(theta, alpha, beta, A):
+    g_BSE = (bse_coefficient_0(A) / pi) * np.cos(alpha) + \
+            (((bse_coefficient(theta, A) - bse_coefficient_0(A)) / scattering_function_correction(theta)) * \
+             ((k + 1) / (2 * pi))) * np.cos(beta) ** k
+    return g_BSE
+
+
+# Iterate over the beam profile
+for profile_x in range(beam_profile.shape[0]):
+    for profile_y in range(beam_profile.shape[1]):
+        for profile_z in range(z_num):
+            domain_x = x_num // 2 - beam_profile.shape[0] // 2 + profile_x
+            domain_y = y_num // 2 - beam_profile.shape[1] // 2 + profile_y
+            if (0 <= domain_x < x_num) and (0 <= domain_y < y_num):
+                if beam_profile_3d[profile_x, profile_y, profile_z] > 0.001:
+                    beam_interactions[profile_z, domain_y, domain_x] += beam_profile_3d[profile_x, profile_y, profile_z]
+
+                    cell_center = np.array([domain_x, domain_y, profile_z]) * dx
+                    v_1 = detector_vertices[0] - cell_center
+                    v_2 = detector_vertices[1] - cell_center
+                    v_3 = detector_vertices[2] - cell_center
+
+                    omega = np.linalg.norm(np.cross(v_1, v_2)) + \
+                            np.linalg.norm(np.cross(v_2, v_3)) + \
+                            np.linalg.norm(np.cross(v_3, v_1))
+
+                    norm_surface_normal = n_s[profile_z, domain_y, domain_x]
+                    if norm_surface_normal.shape[0] == 5:
+                        norm_surface_normal = norm_surface_normal[0]
+
+                    alpha = np.arccos(np.dot(norm_surface_normal, beam_dir) /
+                                      (np.linalg.norm(norm_surface_normal) * np.linalg.norm(beam_dir)))
+                    p_reflected = compute_reflection_vector(beam_dir, norm_surface_normal)
+                    beta = np.arccos(np.dot(p_reflected, norm_surface_normal) /
+                                     (np.linalg.norm(p_reflected) * np.linalg.norm(norm_surface_normal)))
+
+                    g_value = scattering_function(theta, alpha, beta, atomic_number)
+                    N1_d = beam_interactions[profile_z, domain_y, domain_x] * g_value * omega
+                    scattered_electrons[:, domain_x] += N1_d
+
+
+# Ray tracing function
+def ray_tracing(build_surface, detectors, alpha, beta):
+    electrons_per_detector = np.zeros_like(detectors)
+    for i in range(build_surface.shape[0]):
+        for j in range(build_surface.shape[1]):
+            scattering_value = scattering_function(theta, alpha, beta, atomic_number)
+            bse_value = bse_coefficient(theta, atomic_number)
+            correction_value = scattering_function_correction(theta)
+            interaction = scattering_value * bse_value * correction_value * np.random.random()
+            for d in range(detectors.shape[0]):
+                electrons_per_detector[d, i, j] += interaction
+    return electrons_per_detector, scattering_value, bse_value, correction_value
+
+
+# Define parameters based on the paper
+alpha = 0.5  # Example value for alpha (related to Lambert reflection)
+beta = 0.5  # Example value for beta (related to the scattering lobe)
+
+# Initialize build surface and detectors
+build_surface = np.zeros((num_scan_lines, int(build_surface_width / pixel_size)))
+detectors = np.zeros((4, num_scan_lines, int(build_surface_width / pixel_size)))
+
+# Perform ray tracing
+electrons_per_detector, scattering_value, bse_value, correction_value = ray_tracing(build_surface, detectors, alpha,
+                                                                                    beta)
+
+# Convert detected electrons to voltage
+voltage_per_detector = electrons_per_detector * e * transimpedance
+
+# Log the results
+print(f"Scattering Value: {scattering_value}")
+print(f"BSE Coefficient: {bse_value}")
+print(f"Correction Value: {correction_value}")
+print(f"Electrons per detector:\n{electrons_per_detector}")
+print(f"Voltage per detector:\n{voltage_per_detector}")

src/ray_tracer/raytracer.py

+import numpy as np
+import sympy as sp
+from typing import Tuple, List
+from src.materials.alloys.Alloy import Alloy
+from src.materials.alloys import Ti6Al4V
+# from src.materials.alloys.SS316L import SS316L
+
+T = sp.Symbol('T')
+
+# Constants
+pi: float = np.pi  # Mathematical constant π
+e: float = 1.60217662e-19  # Elementary charge in Coulombs
+beam_current: float = 3e-3  # Beam current in Amperes
+dwell_time: float = 0.4e-6  # Dwell time in seconds
+# A: float = 21.5  # Approximation for Ti6Al4V
+Ti6Al4V = Ti6Al4V.create_Ti6Al4V(T)
+A: float = Ti6Al4V.atomic_number
+print('Ti6Al4V.elements: ', Ti6Al4V.elements)
+print('Ti6Al4V.composition: ', Ti6Al4V.composition)
+# SS316L = SS316L.create_SS316L(T)
+# A: float = SS316L.atomic_number
+# print('SS316L.elements: ', SS316L.elements)
+# print('SS316L.composition: ', SS316L.composition)
+print('A: ', A)
+dx: float = 10e-6  # Cell size in meters
+dy: float = 10e-6  # Cell size in meters
+dz: float = 5e-6  # Cell size in meters
+sigma: float = 300e-6 / (2 * dx)  # Beam profile sigma
+threshold: float = 0.0001  # Beam profile threshold
+# The voltage at detector d can be computed using the transimpedance T (in units of V /A ) of a virtual amplifier in the context of the model
+transimpedance: float = 1e5  # Transimpedance in Volts per Ampere
+eV_threshold = 50  # Threshold energy in eV
+threshold_energy = eV_threshold * e  # Threshold energy in joules
+print('threshold_energy: ', threshold_energy)
+
+# Domain setup (Section 2.1.1)
+# Define the 3D domain for the simulation
+x_min, x_max = -500e-6, 500e-6
+y_min, y_max = -250e-6, 250e-6
+z_min, z_max = 0, 20e-6
+x_num = int((x_max - x_min) / dx) + 1
+y_num = int((y_max - y_min) / dy) + 1
+z_num = int((z_max - z_min) / dz) + 1
+domain = np.zeros((z_num, y_num, x_num))  # (z, y , x)
+print('domain shape: ', np.shape(domain))
+domain[0:2, :, :] = 1  # Solid cells
+domain[2, :, :] = 0.5  # Interface layer
+
+
+def compute_surface_normals(_domain: np.ndarray) -> Tuple[np.ndarray, Tuple[np.ndarray, np.ndarray, np.ndarray], int]:
+    """
+    Compute surface normals from the domain.
+
+    Args:
+        _domain (ndarray): The 3D domain array.
+
+    Returns:
+        tuple: A tuple containing the surface normals array and the gradients (grad_x, grad_y, grad_z).
+    """
+    _grad_x = np.gradient(_domain, axis=2)  # x-direction
+    _grad_y = np.gradient(_domain, axis=1)  # y-direction
+    _grad_z = np.gradient(_domain, axis=0)  # z-direction
+    _n_S = np.stack((_grad_x, _grad_y, _grad_z), axis=-1)  # (z, y, x, 3)
+    norm = np.linalg.norm(_n_S, axis=-1)  # (z, y, x)
+    interface_mask = (_domain > 0) & (_domain < 1)
+    nonzero_mask = norm > 0
+    valid_mask = interface_mask & nonzero_mask
+
+    _n_S[valid_mask] /= norm[valid_mask][..., np.newaxis]
+    _n_S[valid_mask] *= -1  # Invert the surface normals to point towards the detectors
+
+    if not np.allclose(np.linalg.norm(_n_S[valid_mask], axis=-1), 1.0, atol=1e-8):
+        raise ValueError("n_S is not normalized.")
+
+    _interface_layer = np.argmax(np.any(interface_mask, axis=(1, 2)))
+
+    return _n_S, (_grad_x, _grad_y, _grad_z), _interface_layer
+
+
+# Compute surface normals (Section 2.1.1)
+n_S, (grad_x, grad_y, grad_z), interface_layer = compute_surface_normals(domain)
+print('np.shape(n_S): ', np.shape(n_S))  # (z, y, x, 3)
+print('interface_layer: ', interface_layer)
+
+# Beam profile (Section 2.1)
+# Define the Gaussian beam profile
+x = np.arange(-30, 31)
+y = np.arange(-30, 31)
+X, Y = np.meshgrid(x, y)
+beam_profile = np.exp(-(X ** 2 + Y ** 2) / (2 * sigma ** 2))
+beam_profile /= np.sum(beam_profile)  # Normalize beam profile
+
+# Detector setup (Section 2.1.1)
+# Define positions of the detectors
+detector_height = 272e-3  # Detector height in meters
+detector_spacing = 127e-3  # Detector spacing in meters
+detector_diameter = 50e-3  # Detector diameter in meters
+detector_positions = np.array([
+    [detector_spacing, 0, detector_height],  # Right detector
+    [-detector_spacing, 0, detector_height],  # Left detector
+    [0, detector_spacing, detector_height],  # Back detector
+    [0, -detector_spacing, detector_height]  # Front detector
+])
+
+
+# Define detector vertices for solid angle calculation
+def get_detector_vertices(detector_position: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """
+    Get the vertices of an equilateral triangle detector surface.
+
+    Parameters:
+        detector_position (array-like): The position of the center of the detector.
+
+    Returns:
+        tuple: The vertices of the detector surface.
+    """
+    print('detector_positions: ', detector_position)
+    print('detector_position type: ', type(detector_position))
+    # Calculate the side length of the detector
+    side_length = np.sqrt(3) * detector_diameter / 2
+    height = (np.sqrt(3) / 2) * side_length
+
+    _v1 = detector_position + np.array([0, 2 / 3 * height, 0])
+    _v2 = detector_position + np.array([-0.5 * side_length, -1 / 3 * height, 0])
+    _v3 = detector_position + np.array([0.5 * side_length, -1 / 3 * height, 0])
+
+    # print("Vertices:", v1, v2, v3)
+    # print("Type of vertices:", type(v1), type(v2), type(v3))
+
+    return _v1, _v2, _v3
+
+
+# Compute solid angle using Equation A1 (Appendix A)
+def compute_solid_angle(vector1: np.ndarray, vector2: np.ndarray, vector3: np.ndarray) -> float:
+    """
+    Compute the solid angle of a 3D triangle.
+
+    Args:
+        vector1: The first vector of the triangle.
+        vector2: The second vector of the triangle.
+        vector3: The third vector of the triangle.
+
+    Returns:
+        The solid angle of the triangle.
+    """
+    # Input validation
+    if not isinstance(vector1, np.ndarray) or not isinstance(vector2, np.ndarray) or not isinstance(vector3, np.ndarray):
+        raise ValueError("Inputs must be numpy arrays")
+    if vector1.shape != (3,) or vector2.shape != (3,) or vector3.shape != (3,):
+        raise ValueError("Inputs must be 3D vectors")
+    if not np.isfinite(vector1).all() or not np.isfinite(vector2).all() or not np.isfinite(vector3).all():
+        raise ValueError("Inputs must be finite")
+
+    cross_product = np.cross(vector2, vector3)
+    # Check if the cross product is nearly zero
+    if np.linalg.norm(cross_product) < 1e-8:
+        raise ValueError("Vectors are nearly coplanar")
+
+    # Compute solid angle
+    try:
+        # Compute the numerator of the solid angle formula
+        numerator: float = np.dot(vector1, cross_product)
+        # Check if the numerator is nearly zero
+        if np.abs(numerator) < 1e-8:
+            raise ValueError("Vectors are nearly parallel")
+
+        # Compute the denominator of the solid angle formula
+        denominator: float = (
+                np.linalg.norm(vector1) * np.linalg.norm(vector2) * np.linalg.norm(vector3)
+                + np.dot(vector1, vector2) * np.linalg.norm(vector3)
+                + np.dot(vector1, vector3) * np.linalg.norm(vector2)
+                + np.dot(vector2, vector3) * np.linalg.norm(vector1)
+        )
+
+        # Compute the solid angle in radians
+        solid_angle = 2 * np.arctan2(numerator, denominator)
+        # print("Solid angle:", solid_angle)
+        return solid_angle
+    except ZeroDivisionError:
+        raise ValueError("Error computing solid angle: division by zero")
+    except np.linalg.LinAlgError:
+        raise ValueError("Error computing solid angle: linear algebra error")
+
+
+# BSE coefficient functions (Equations 11 and 10, Section 2.1.2)
+def eta_0(_A: float) -> float:
+    """Equation 11 from the paper"""
+    _eta_0: float = -0.0254 + 0.016 * _A - 1.86e-4 * _A ** 2 + 8.3e-7 * _A ** 3
+    # print('eta_0: ', _eta_0)
+    return _eta_0
+
+
+# BSE Coefficient for single atomic targets, found to hold in the range of 10–100 keV
+def eta(theta: float, _A: float) -> float:
+    """Equation 10 from the paper"""
+    B: float = 0.89
+    _eta: float = B * (eta_0(_A) / B) ** np.cos(theta)
+    # print('eta: ', _eta)
+    return _eta
+
+
+# Compute the weighted average BSE coefficient for the alloy
+def compute_weighted_eta(theta: float, alloy: Alloy) -> float:
+    """Equation 12 from the paper"""
+    weighted_eta: float = 0
+    for element, weight_fraction in zip(alloy.elements, alloy.composition):
+        # print('element: ', element)
+        # print('weight_fraction: ', weight_fraction)
+        _A: float = element.atomic_number
+        # print('A: ', _A)
+        eta_i: float = eta(theta, _A)
+        # print('eta_i: ', eta_i)
+        weighted_eta += weight_fraction * eta_i
+        # print('weighted_eta: ', weighted_eta)
+    return weighted_eta
+
+
+# Correction function C(theta) (Equation 14, Section 2.1.3)
+def C(theta: float) -> float:
+    """Equation 14 from the paper"""
+    # theta = np.rad2deg(theta) # Ensure angle is in degrees for C function (?)
+    _C: float = 1.01 - 4.06 * theta + 1.36e-5 * theta ** 2 - 1.71e-6 * theta ** 3
+    # print('C: ', _C)
+    # Check if C is negative and raise an exception
+    if _C < 0:
+        raise ValueError(f"C function returned a negative value for theta: {theta}. "
+                         f"Please check the input and the C function implementation.")
+    return _C
+
+
+# Scattering function g_BSE (Equation 13, Section 2.1.3)
+def g_BSE(_theta: float, _alpha: float, _beta: float) -> float:
+    """Equation 13 from the paper"""
+    # k: float = np.rad2deg(_theta) / 13  # Correct implementation as per the paper
+    diffusive_part: float = (eta_0(A) / pi) * np.cos(_alpha)
+    # reflective_part: float = ((eta(_theta, A) - eta_0(A)) / C(_theta)) * ((k + 1) / (2 * pi)) * (np.cos(_beta) ** k)
+
+    # Compute a weighted average BSE coefficient for the alloy
+    # reflective_part: float = ((compute_weighted_eta(_theta, Ti6Al4V) - eta_0(A)) / C(_theta)) * ((k + 1) / (2 * pi)) * (np.cos(_beta) ** k)
+
+    # print('g_BSE: ', diffusive_part + reflective_part)
+    return diffusive_part #+ reflective_part
+
+
+# Beam direction
+phi: float = 0.0  # Angle in degrees (for phi > ~14.25°, C(theta) becomes negative)
+phi_rad: float = np.deg2rad(phi)  # Convert to radians
+p_E = np.array([np.sin(phi_rad), 0, -np.cos(phi_rad)])  # Beam direction vector
+p_E /= np.linalg.norm(p_E)
+
+# Initialize total electrons counter and electrons per detector array
+total_electrons: int = 0
+electrons_per_detector = np.zeros(len(detector_positions))
+beam_loc = (x_num // 2, y_num // 2)  # (x, y)
+print('beam_loc: ', beam_loc)
+
+# Iterate over the beam profile
+for profile_x in range(beam_profile.shape[0]):
+    for profile_y in range(beam_profile.shape[1]):
+        # print(f"beam_profile.shape ({profile_x}, {profile_y}):")
+        domain_x = beam_loc[0] - beam_profile.shape[0] // 2 + profile_x
+        domain_y = beam_loc[1] - beam_profile.shape[1] // 2 + profile_y
+        # print(f"Domain ({domain_x}, {domain_y}):")
+        if (0 <= domain_x < x_num) and (0 <= domain_y < y_num):
+            if beam_profile[profile_x, profile_y] > threshold:
+                n_S_local = n_S[interface_layer, domain_y, domain_x]  # Interface layer
+
+                # Debugging print statements
+                print(f"grad_x: {grad_x[2, domain_y, domain_x]}")
+                print(f"grad_y: {grad_y[2, domain_y, domain_x]}")
+                print(f"grad_z: {grad_z[2, domain_y, domain_x]}")
+                print(f"n_S_local: {n_S_local}")
+
+                if not np.allclose(np.linalg.norm(n_S_local), 1.0):
+                    raise ValueError("n_S_local is not normalized.")
+
+                if np.allclose(n_S_local, np.array([0, 0, 1])) or np.allclose(n_S_local, np.array([0, 0, -1])):
+                    print('hellooo')
+                    theta_local: float = phi_rad
+                else:
+                    theta_local: float = np.arccos(
+                        np.clip(np.dot(p_E, n_S_local) / (np.linalg.norm(p_E) * np.linalg.norm(n_S_local)),
+                                -1.0, 1.0)
+                    )
+
+                for det_idx, det_pos in enumerate(detector_positions):
+                    det_idx: int
+                    det_pos: np.ndarray[int]
+
+                    det_vec = det_pos - np.array([domain_x * dx, domain_y * dx, 0])
+                    det_vec /= np.linalg.norm(det_vec)
+
+                    # Debugging print statements
+                    print(f"Detector Position: {det_pos}")
+                    print(f"Detector Position Type: {type(det_pos)}")
+                    print(f"Domain Point: {np.array([domain_x * dx, domain_y * dx, 0])}")
+                    print(f"Detector Vector: {det_vec}")
+                    print(f"Dot product: {np.dot(det_vec, n_S_local)}")
+
+                    alpha = np.arccos(np.clip(np.dot(det_vec, n_S_local), -1.0, 1.0))
+
+                    p_E_reflected = p_E - 2 * np.dot(p_E, n_S_local) * n_S_local
+                    p_E_reflected /= np.linalg.norm(p_E_reflected)
+
+                    beta = np.arccos(np.clip(np.dot(p_E_reflected, det_vec), -1.0, 1.0))
+
+                    # Apply energy threshold
+                    # electron_energy = N_0 * e * g_bse_value
+                    # print('electron_energy: ', electron_energy)
+
+                    # Calculate the initial energy of the primary electrons (E_0)
+                    V_acc = 60e3  # Accelerating voltage in Volts
+                    E_0 = V_acc * e  # Initial energy of primary electrons in joules
+
+                    # Calculate the energy of the backscatter electrons
+                    # E_BSE = E_0 * eta(theta_local, A)
+                    # Compute a weighted average BSE coefficient for the alloy
+                    E_BSE = E_0 * compute_weighted_eta(theta_local, Ti6Al4V)
+                    print('E_BSE: ', E_BSE)
+
+                    # Bifurcation logic for PE and SE
+                    if E_BSE > threshold_energy:
+                        # Backscatter electrons (BSE)
+                        # Compute scattering function g_bse_value
+                        g_bse_value = g_BSE(theta_local, alpha, beta)
+
+                        # Ensure g_bse_value is non-negative
+                        if g_bse_value < 0:
+                            raise ValueError("g_bse_value is negative, check calculations.")
+
+                        # Compute solid angle (Equation A1)
+                        v1, v2, v3 = get_detector_vertices(det_pos)
+                        print("Vertices:", v1, v2, v3)
+                        print("Type of vertices:", type(v1), type(v2), type(v3))
+                        dOmega = compute_solid_angle(v1, v2, v3)
+
+                        # Compute number of electrons collected by the detector (Equation A2)
+                        # The beam ray stays on each grid point ij for the dwell time t_d and distributes N_0 PE into it
+                        N_0 = beam_current * dwell_time / e  # Equation 3, Section 2.1.1
+
+                        # There are N1 electrons in the chamber after the 1st scattering event
+                        # N_1 = N_0 * eta(theta_local, A)  # Equation 4, Section 2.1.1
+                        # Compute a weighted average BSE coefficient for the alloy
+                        N_1 = N_0 * compute_weighted_eta(theta_local, Ti6Al4V)  # Equation 4, Section 2.1.1
+
+                        # The number of electrons N_1_d collected by a general detector d
+                        N_1_d = N_1 * np.sum(g_bse_value * dOmega)  # Equation 5, Section 2.1.1
+
+                        # Accumulate electrons per detector
+                        electrons_per_detector[det_idx] += N_1_d
+
+                        # Add to total electrons
+                        total_electrons += N_1_d
+
+                        # Debugging print statements
+                        print(f"Domain ({domain_x}, {domain_y}):")
+                        print(f"n_S_local: {n_S_local}")
+                        print(f"theta_local: {np.rad2deg(theta_local)}, "
+                              f"alpha: {np.rad2deg(alpha)}, "
+                              f"beta: {np.rad2deg(beta)}")
+                        print(f"det_vec: {det_vec}")
+                        print(f"p_E_reflected: {p_E_reflected}")
+                        print(f"g_bse_value: {g_bse_value}, N_0: {N_0}, N_1: {N_1}, N_1_d: {N_1_d}")
+                        print(f"total_electrons: {total_electrons}")
+                    else:
+                        # Secondary electrons (SE)
+                        # Note: This part is currently neglected in the model
+                        # but can be implemented if necessary
+                        continue  # Skip to the next iteration if energy is below threshold
+
+# Initialize build surface and detectors (Section 2.2)
+num_scan_lines = 1500  # y (rows)
+build_surface_width = 80e-3  # x (columns)
+pixel_size = 53.3e-6
+build_surface = np.zeros((num_scan_lines, int(build_surface_width / pixel_size)))  # (y, x) -> (rows, columns)
+print('build_surface: ', np.shape(build_surface))  # (1500, 1500)
+detectors = np.zeros((4, num_scan_lines, int(build_surface_width / pixel_size)))  # (4, y, x)
+
+
+# Ray tracing function
+def ray_tracing(_build_surface, _detectors, _interface_layer):
+    _electrons_per_detector_ray_tracing = np.zeros_like(_detectors)
+    # Iterate over the scan lines and simulate ray interactions
+    for _j in range(_build_surface.shape[0]):  # j -> y
+        for _i in range(_build_surface.shape[1]):  # i -> x
+            # Calculate the ray direction using detector vectors
+            for _det_idx, _det_pos in enumerate(detector_positions):
+                if _j < n_S.shape[1] and _i < n_S.shape[2]:  # Ensure indices are within bounds
+                    _n_S_local = n_S[_interface_layer, _j, _i, :]  # Ensure correct surface normal is used
+
+                    if not np.allclose(np.linalg.norm(_n_S_local), 1.0):
+                        raise ValueError("n_S_local is not normalized.")
+
+                    if np.allclose(_n_S_local, np.array([-0, -0, 1])):
+                        _theta_local = phi_rad
+                    else:
+                        _theta_local = np.arccos(
+                            np.clip(np.dot(_n_S_local, p_E) / (np.linalg.norm(_n_S_local) * np.linalg.norm(p_E)),
+                                    -1.0, 1.0))
+
+                    _det_vec = _det_pos - np.array([_i * pixel_size, _j * pixel_size, 0])
+                    _det_vec /= np.linalg.norm(_det_vec)
+
+                    _alpha = np.arccos(np.clip(np.dot(_det_vec, _n_S_local), -1.0, 1.0))
+
+                    _p_E_reflected = p_E - 2 * np.dot(p_E, _n_S_local) * _n_S_local
+                    _p_E_reflected /= np.linalg.norm(_p_E_reflected)
+
+                    _beta = np.arccos(np.clip(np.dot(_p_E_reflected, _det_vec), -1.0, 1.0))
+
+                    # Compute scattering function
+                    scattering_value = g_BSE(_theta_local, _alpha, _beta)
+                    # Apply the BSE coefficient
+                    # bse_value = eta(_theta_local, A)
+                    # Compute a weighted average BSE coefficient for the alloy
+                    bse_value = compute_weighted_eta(_theta_local, Ti6Al4V)
+                    # Apply scattering function correction
+                    correction_value = C(_theta_local)
+                    interaction = scattering_value * bse_value * correction_value
+                    _electrons_per_detector_ray_tracing[_det_idx, _j, _i] += interaction
+
+                    # print(electrons_per_detector)
+                    # print(electrons_per_detector.shape)
+
+                    # Debugging print statements
+                    print(f"Ray Tracing - Scan line ({_i}, {_j}):")
+                    print(f"n_S_local: {_n_S_local}")
+                    print(f"theta_local: {np.rad2deg(_theta_local)}, "
+                          f"alpha: {np.rad2deg(_alpha)}, "
+                          f"beta: {np.rad2deg(_beta)}")
+                    print(f"det_vec: {_det_vec}")
+                    print(f"p_E_reflected: {_p_E_reflected}")
+                    print(f"scattering_value: {scattering_value}, "
+                          f"bse_value: {bse_value}, "
+                          f"correction_value: {correction_value}")
+                    print(f"interaction: {interaction}")
+
+    return _electrons_per_detector_ray_tracing, np.sum(_electrons_per_detector_ray_tracing)
+
+
+# Perform ray tracing
+electrons_per_detector_ray_tracing, total_electrons_ray_tracing = ray_tracing(build_surface, detectors, interface_layer)
+
+# Convert detected electrons to voltage using the transimpedance amplifier
+voltage_per_detector = electrons_per_detector * e * transimpedance
+
+# Print results
+print(f"Electrons per detector (ray_tracing):\n {electrons_per_detector_ray_tracing}")
+for i, electrons in enumerate(electrons_per_detector_ray_tracing):
+    print(f"Electrons hitting detector {i + 1} (ray_tracing): {np.sum(electrons)}")
+print(f"Total electrons hitting all detectors (ray_tracing): {total_electrons_ray_tracing}")
+# print(f"Voltage per detector (ray tracing):\n {voltage_per_detector}")
+print(f"Electrons per detector: {electrons_per_detector}")
+for i, electrons in enumerate(electrons_per_detector):
+    print(f"Electrons hitting detector {i + 1}: {np.sum(electrons)}")
+print(f"Total electrons hitting all detectors: {total_electrons}")
+
+
+# Test functions
+def test_eta_0() -> None:
+    assert np.isclose(eta_0(21.5), 0.24087)
+
+
+def test_eta() -> None:
+    theta: float = np.deg2rad(10)
+    assert np.isclose(eta(theta, A), 0.24570)
+
+
+def test_C() -> None:
+    assert np.isclose(C(0), 1.01)
+    theta: float = np.deg2rad(10)
+    assert np.isclose(C(theta), 0.30139672887864666)
+
+
+def test_g_BSE() -> None:
+    _theta: float = np.deg2rad(10)
+    _alpha: float = np.deg2rad(25)
+    _beta: float = np.deg2rad(25)
+    assert np.isclose(g_BSE(_theta, _alpha, _beta), 0.07367186921200412)
+
+
+def test_compute_solid_angle() -> None:
+    _v1: np.ndarray = np.array([1, 0, 0])
+    _v2: np.ndarray = np.array([0, 1, 0])
+    _v3: np.ndarray = np.array([0, 0, 1])
+    solid_angle: float = compute_solid_angle(_v1, _v2, _v3)
+    assert np.isclose(solid_angle, np.pi / 2)
+
+
+def test_incident_angles() -> None:
+    incident_angles: List[float] = [0, 2, 4, 6, 8, 10]  # in degrees
+    for _theta in incident_angles:
+        _theta_rad: float = np.deg2rad(_theta)
+        _p_E: np.ndarray = np.array([np.sin(_theta_rad), 0, -np.cos(_theta_rad)])
+        _p_E /= np.linalg.norm(_p_E)
+
+        _eta_val: float = eta(_theta_rad, A)
+        _C_val: float = C(_theta_rad)
+        # print('C_val: ', C_val)
+
+        for _alpha in np.linspace(0, np.pi / 4, 5):
+            for _beta in np.linspace(0, np.pi / 4, 5):
+                _g_BSE_val: float = g_BSE(_theta_rad, _alpha, _beta)
+                assert 0 < _eta_val < 1
+                assert 0 < _C_val <= 1.01
+                assert _g_BSE_val >= 0
+
+
+def test_vector_calculations() -> None:
+    # Test reflection
+    _p_E = np.array([0, 0, -1])
+    _n_S = np.array([0, 1, 0])
+    dot_product = np.dot(_p_E, _n_S)
+    reflection_vector = _p_E - 2 * dot_product * _n_S
+    _p_E_reflected = reflection_vector / np.linalg.norm(reflection_vector)
+
+    # Check if the reflected vector is orthogonal to the surface normal
+    assert np.isclose(np.dot(_p_E_reflected, _n_S), 0, atol=1e-8)
+    # Check if the magnitude of the reflected vector is correct
+    assert np.isclose(np.linalg.norm(_p_E_reflected), 1, atol=1e-8)
+
+    # Test detector vectors
+    detector_pos = np.array([1, 1, 1], dtype=np.float64)
+    domain_point = np.array([0, 0, 0], dtype=np.float64)
+    detector_vector = detector_pos - domain_point
+    detector_vector = detector_vector.astype(np.float64)  # Ensure floating-point type
+    detector_vector /= np.linalg.norm(detector_vector)
+    assert np.allclose(detector_vector, [0.57735027, 0.57735027, 0.57735027])
+
+    # Test angles between vectors
+    vec1 = np.array([1, 0, 0], dtype=np.float64)
+    vec2 = np.array([1, 1, 0], dtype=np.float64)  # Ensure floating-point type
+    vec2 /= np.linalg.norm(vec2)
+    angle = np.arccos(np.dot(vec1, vec2))
+    assert np.isclose(np.rad2deg(angle), 45)
+
+
+def test_compute_surface_normals() -> None:
+    _domain = np.zeros((5, 10, 10))
+    _domain[0:2, :, :] = 1
+    _domain[2, :, :] = 0.5
+    _n_S, _, _ = compute_surface_normals(_domain)
+
+    # Check normalization for valid mask
+    '''valid_mask = (domain > 0) & (domain < 1)
+    nonzero_mask = np.linalg.norm(n_S, axis=-1) > 0
+    valid_mask = valid_mask & nonzero_mask'''
+
+    # Check if the surface normals are normalized
+    for _i in range(_n_S.shape[1]):
+        for _j in range(_n_S.shape[2]):
+            # if valid_mask[2, _i, _j]:
+            assert np.isclose(np.linalg.norm(_n_S[2, _i, _j]), 1.0, atol=1e-8)
+
+            # Check if the normal is close to [0, 0, 1] or [0, 0, -1]
+            assert (np.allclose(_n_S[2, _i, _j], [0, 0, 1], atol=1e-8) or
+                    np.allclose(_n_S[2, _i, _j], [0, 0, -1], atol=1e-8))
+
+
+def main() -> None:
+    test_eta_0()
+    test_eta()
+    test_C()
+    test_g_BSE()
+    test_compute_solid_angle()
+    test_incident_angles()
+    test_vector_calculations()
+    test_compute_surface_normals()
+    print("All tests passed!")
+
+
+if __name__ == "__main__":
+    weighted_eta_value = compute_weighted_eta(np.deg2rad(10), Ti6Al4V)
+    print("Weighted eta value:", weighted_eta_value)
+    # main()

src/ray_tracer/raytracerMM.py

+import numpy as np
+import sympy as sp
+from typing import Tuple
+from src.materials.alloys.Alloy import Alloy
+from src.materials.alloys import Ti6Al4V
+# from src.materials.alloys.SS316L import SS316L
+
+T = sp.Symbol('T')
+
+# Constants
+pi: float = np.pi  # Mathematical constant π
+e: float = 1.60217662e-19  # Elementary charge in Coulombs
+beam_current: float = 3e-3  # Beam current in Amperes
+dwell_time: float = 0.4e-6  # Dwell time in seconds
+# A: float = 21.5  # Approximation for Ti6Al4V
+Ti6Al4V = Ti6Al4V.create_Ti6Al4V(T)
+A: float = Ti6Al4V.atomic_number
+print('Ti6Al4V.elements: ', Ti6Al4V.elements)
+print('Ti6Al4V.composition: ', Ti6Al4V.composition)
+# SS316L = SS316L.create_SS316L(T)
+# A: float = SS316L.atomic_number
+# print('SS316L.elements: ', SS316L.elements)
+# print('SS316L.composition: ', SS316L.composition)
+print('A: ', A)
+dx: float = 10e-6  # Cell size in meters
+dy: float = 10e-6  # Cell size in meters
+dz: float = 5e-6  # Cell size in meters
+sigma: float = 300e-6 / (2 * dx)  # Beam profile sigma
+threshold: float = 0.0001  # Beam profile threshold
+# The voltage at detector d can be computed using the transimpedance T (in units of V /A ) of a virtual amplifier in the context of the model
+transimpedance: float = 1e5  # Transimpedance in Volts per Ampere
+eV_threshold = 50  # Threshold energy in eV
+threshold_energy = eV_threshold * e  # Threshold energy in joules
+print('threshold_energy: ', threshold_energy)
+
+# Domain setup (Section 2.1.1)
+# Define the 3D domain for the simulation
+x_min, x_max = -250e-6, 250e-6
+y_min, y_max = -250e-6, 250e-6
+z_min, z_max = 0, 20e-6
+x_num = int((x_max - x_min) / dx) + 1
+y_num = int((y_max - y_min) / dy) + 1
+z_num = int((z_max - z_min) / dz) + 1
+domain = np.zeros((z_num, y_num, x_num))  # (z, y , x)
+print('domain shape: ', np.shape(domain))
+domain[0:2, :, :] = 1  # Solid cells
+domain[2, :, :] = 0.5  # Interface layer
+
+
+def compute_surface_normals(_domain: np.ndarray) -> Tuple[np.ndarray, Tuple[np.ndarray, np.ndarray, np.ndarray], int]:
+    """
+    Compute surface normals from the domain.
+
+    Args:
+        _domain (ndarray): The 3D domain array.
+
+    Returns:
+        tuple: A tuple containing the surface normals array and the gradients (grad_x, grad_y, grad_z).
+    """
+    _grad_x = np.gradient(_domain, axis=2)  # x-direction
+    _grad_y = np.gradient(_domain, axis=1)  # y-direction
+    _grad_z = np.gradient(_domain, axis=0)  # z-direction
+    _n_S = np.stack((_grad_x, _grad_y, _grad_z), axis=-1)  # (z, y, x, 3)
+    norm = np.linalg.norm(_n_S, axis=-1)  # (z, y, x)
+    interface_mask = (_domain > 0) & (_domain < 1)
+    nonzero_mask = norm > 0
+    valid_mask = interface_mask & nonzero_mask
+
+    _n_S[valid_mask] /= norm[valid_mask][..., np.newaxis]
+    _n_S[valid_mask] *= -1  # Invert the surface normals to point towards the detectors
+
+    if not np.allclose(np.linalg.norm(_n_S[valid_mask], axis=-1), 1.0, atol=1e-8):
+        raise ValueError("n_S is not normalized.")
+
+    _interface_layer = np.argmax(np.any(interface_mask, axis=(1, 2)))
+
+    return _n_S, (_grad_x, _grad_y, _grad_z), _interface_layer
+
+
+# Compute surface normals (Section 2.1.1)
+n_S, (grad_x, grad_y, grad_z), interface_layer = compute_surface_normals(domain)
+print('np.shape(n_S): ', np.shape(n_S))  # (z, y, x, 3)
+print('interface_layer: ', interface_layer)
+
+# Beam profile (Section 2.1)
+# Define the Gaussian beam profile
+x = np.arange(-30, 31)
+y = np.arange(-30, 31)
+X, Y = np.meshgrid(x, y)
+beam_profile = np.exp(-(X ** 2 + Y ** 2) / (2 * sigma ** 2))
+beam_profile /= np.sum(beam_profile)  # Normalize beam profile
+
+# Detector setup (Section 2.1.1)
+# Define positions of the detectors
+detector_height = 272e-3  # Detector height in meters
+detector_spacing = 127e-3  # Detector spacing in meters
+detector_diameter = 50e-3  # Detector diameter in meters
+detector_positions = np.array([
+    [detector_spacing, 0, detector_height],  # Right detector
+    [-detector_spacing, 0, detector_height],  # Left detector
+    [0, detector_spacing, detector_height],  # Back detector
+    [0, -detector_spacing, detector_height]  # Front detector
+])
+
+
+# Define detector vertices for solid angle calculation
+def get_detector_vertices(detector_position: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """
+    Get the vertices of an equilateral triangle detector surface oriented symmetrically.
+
+    Parameters:
+    detector_position (array-like): The position of the center of the detector.
+
+    Returns:
+    tuple: The vertices of the detector surface.
+    """
+    # Calculate the side length of the detector
+    side_length = np.sqrt(3) * detector_diameter / 2
+    height = (np.sqrt(3) / 2) * side_length
+
+    # Determine the orientation based on the detector position
+    if detector_position[0] > 0:  # Right detector
+        _v1 = detector_position + np.array([2/3*height, 0, 0])
+        _v2 = detector_position + np.array([-height/3, side_length/2, 0])
+        _v3 = detector_position + np.array([-height/3, -side_length/2, 0])
+    elif detector_position[0] < 0:  # Left detector
+        _v1 = detector_position + np.array([-2/3*height, 0, 0])
+        _v2 = detector_position + np.array([height/3, -side_length/2, 0])
+        _v3 = detector_position + np.array([height/3, side_length/2, 0])
+    elif detector_position[1] > 0:  # Back detector
+        _v1 = detector_position + np.array([0, 2/3*height, 0])
+        _v2 = detector_position + np.array([-side_length/2, -height/3, 0])
+        _v3 = detector_position + np.array([side_length/2, -height/3, 0])
+    else:  # Front detector
+        _v1 = detector_position + np.array([0, -2/3*height, 0])
+        _v2 = detector_position + np.array([side_length/2, height/3, 0])
+        _v3 = detector_position + np.array([-side_length/2, height/3, 0])
+
+    return _v1, _v2, _v3
+
+
+# Compute solid angle using Equation A1 (Appendix A)
+def compute_solid_angle(vector1: np.ndarray, vector2: np.ndarray, vector3: np.ndarray) -> float:
+    """
+    Compute the solid angle of a 3D triangle.
+
+    Args:
+        vector1: The first vector of the triangle.
+        vector2: The second vector of the triangle.
+        vector3: The third vector of the triangle.
+
+    Returns:
+        The solid angle of the triangle.
+    """
+    # Input validation
+    if not isinstance(vector1, np.ndarray) or not isinstance(vector2, np.ndarray) or not isinstance(vector3, np.ndarray):
+        raise ValueError("Inputs must be numpy arrays")
+    if vector1.shape != (3,) or vector2.shape != (3,) or vector3.shape != (3,):
+        raise ValueError("Inputs must be 3D vectors")
+    if not np.isfinite(vector1).all() or not np.isfinite(vector2).all() or not np.isfinite(vector3).all():
+        raise ValueError("Inputs must be finite")
+
+    cross_product = np.cross(vector2, vector3)
+    # Check if the cross product is nearly zero
+    if np.linalg.norm(cross_product) < 1e-8:
+        raise ValueError("Vectors are nearly coplanar")
+
+    # Compute solid angle
+    try:
+        # Compute the numerator of the solid angle formula
+        numerator: float = np.dot(vector1, cross_product)
+        # Check if the numerator is nearly zero
+        if np.abs(numerator) < 1e-8:
+            raise ValueError("Vectors are nearly parallel")
+
+        # Compute the denominator of the solid angle formula
+        denominator: float = (
+                np.linalg.norm(vector1) * np.linalg.norm(vector2) * np.linalg.norm(vector3)
+                + np.dot(vector1, vector2) * np.linalg.norm(vector3)
+                + np.dot(vector1, vector3) * np.linalg.norm(vector2)
+                + np.dot(vector2, vector3) * np.linalg.norm(vector1)
+        )
+
+        # Compute the solid angle in radians
+        solid_angle = 2 * np.arctan2(numerator, denominator)
+        # print("Solid angle:", solid_angle)
+        return solid_angle
+    except ZeroDivisionError:
+        raise ValueError("Error computing solid angle: division by zero")
+    except np.linalg.LinAlgError:
+        raise ValueError("Error computing solid angle: linear algebra error")
+
+
+# BSE coefficient functions (Equations 11 and 10, Section 2.1.2)
+def eta_0(_A: float) -> float:
+    """Equation 11 from the paper"""
+    _eta_0: float = -0.0254 + 0.016 * _A - 1.86e-4 * _A ** 2 + 8.3e-7 * _A ** 3
+    # print('eta_0: ', _eta_0)
+    return _eta_0
+
+
+# BSE Coefficient for single atomic targets, found to hold in the range of 10–100 keV
+def eta(theta: float, _A: float) -> float:
+    """Equation 10 from the paper"""
+    B: float = 0.89
+    _eta: float = B * (eta_0(_A) / B) ** np.cos(theta)
+    # print('eta: ', _eta)
+    return _eta
+
+
+# Compute the weighted average BSE coefficient for the alloy
+def compute_weighted_eta(theta: float, alloy: Alloy) -> float:
+    """Equation 12 from the paper"""
+    weighted_eta: float = 0
+    for element, weight_fraction in zip(alloy.elements, alloy.composition):
+        # print('element: ', element)
+        # print('weight_fraction: ', weight_fraction)
+        _A: float = element.atomic_number
+        # print('A: ', _A)
+        eta_i: float = eta(theta, _A)
+        # print('eta_i: ', eta_i)
+        weighted_eta += weight_fraction * eta_i
+        # print('weighted_eta: ', weighted_eta)
+    return weighted_eta
+
+
+# Correction function C(theta) (Equation 14, Section 2.1.3)
+def C(theta: float) -> float:
+    """Equation 14 from the paper"""
+    theta = np.rad2deg(theta) # Ensure angle is in degrees for C function
+    _C: float = 1.01 - 4.06 * theta + 1.36e-5 * theta ** 2 - 1.71e-6 * theta ** 3
+    # print('C: ', _C)
+    # Check if C is negative and raise an exception
+    if _C < 0:
+        raise ValueError(f"C function returned a negative value for theta: {theta}. "
+                         f"Please check the input and the C function implementation.")
+    return _C
+
+
+# Scattering function g_BSE (Equation 13, Section 2.1.3)
+def g_BSE(_theta: float, _alpha: float, _beta: float) -> float:
+    """Equation 13 from the paper"""
+    # k: float = np.rad2deg(_theta) / 13  # Correct implementation as per the paper
+    diffusive_part: float = (eta_0(A) / pi) * np.cos(_alpha)
+    # reflective_part: float = ((eta(_theta, A) - eta_0(A)) / C(_theta)) * ((k + 1) / (2 * pi)) * (np.cos(_beta) ** k)
+
+    # Compute a weighted average BSE coefficient for the alloy
+    # reflective_part: float = ((compute_weighted_eta(_theta, Ti6Al4V) - eta_0(A)) / C(_theta)) * ((k + 1) / (2 * pi)) * (np.cos(_beta) ** k)
+
+    # print('g_BSE: ', diffusive_part + reflective_part)
+    return diffusive_part #+ reflective_part
+
+
+# Beam direction
+phi: float = 0.0  # Angle in degrees (for phi > ~14.25°, C(theta) becomes negative)
+phi_rad: float = np.deg2rad(phi)  # Convert to radians
+p_E = np.array([np.sin(phi_rad), 0, -np.cos(phi_rad)])  # Beam direction vector
+p_E /= np.linalg.norm(p_E)
+
+# Initialize total electrons counter and electrons per detector array
+total_electrons: int = 0
+electrons_per_detector = np.zeros(len(detector_positions))
+print('x_num: ', x_num)
+print('y_num: ', y_num)
+beam_loc = ((x_num+1) // 2, (y_num+1) // 2)  # (x, y)
+print('beam_loc: ', beam_loc)
+
+# loop over beam_loc from x_num // 4 to (x_num*3)//4
+#    compute phi depending on beam origin and beam_loc
+
+# Iterate over the beam profile
+for profile_x in range(beam_profile.shape[0]):
+    for profile_y in range(beam_profile.shape[1]):
+        # print(f"beam_profile.shape ({profile_x}, {profile_y}):")
+        domain_x = beam_loc[0] - beam_profile.shape[0] // 2 + profile_x
+        domain_y = beam_loc[1] - beam_profile.shape[1] // 2 + profile_y
+        # print(f"Domain ({domain_x}, {domain_y}):")
+        if (0 <= domain_x < x_num) and (0 <= domain_y < y_num):
+            if beam_profile[profile_x, profile_y] > threshold:
+                n_S_local = n_S[interface_layer, domain_y, domain_x]  # Interface layer
+
+                # Debugging print statements
+                print(f"grad_x: {grad_x[2, domain_y, domain_x]}")
+                print(f"grad_y: {grad_y[2, domain_y, domain_x]}")
+                print(f"grad_z: {grad_z[2, domain_y, domain_x]}")
+                print(f"n_S_local: {n_S_local}")
+
+                if not np.allclose(np.linalg.norm(n_S_local), 1.0):
+                    raise ValueError("n_S_local is not normalized.")
+
+                if np.allclose(n_S_local, np.array([0, 0, 1])) or np.allclose(n_S_local, np.array([0, 0, -1])):
+                    print('hellooo')
+                    theta_local: float = phi_rad
+                else:
+                    theta_local: float = np.arccos(
+                        np.clip(np.dot(p_E, n_S_local) / (np.linalg.norm(p_E) * np.linalg.norm(n_S_local)),
+                                -1.0, 1.0)
+                    )
+
+                for det_idx, det_pos in enumerate(detector_positions):
+                    det_idx: int
+                    det_pos: np.ndarray[int]
+
+                    det_vec = det_pos - np.array([domain_x * dx, domain_y * dx, 0])
+                    det_vec /= np.linalg.norm(det_vec)
+
+                    # Debugging print statements
+                    print(f"Detector Position: {det_pos}")
+                    print(f"Detector Position Type: {type(det_pos)}")
+                    print(f"Domain Point: {np.array([domain_x * dx, domain_y * dx, 0])}")
+                    print(f"Detector Vector: {det_vec}")
+                    print(f"Dot product: {np.dot(det_vec, n_S_local)}")
+
+                    alpha = np.arccos(np.clip(np.dot(det_vec, n_S_local), -1.0, 1.0))
+
+                    p_E_reflected = p_E - 2 * np.dot(p_E, n_S_local) * n_S_local
+                    p_E_reflected /= np.linalg.norm(p_E_reflected)
+
+                    beta = np.arccos(np.clip(np.dot(p_E_reflected, det_vec), -1.0, 1.0))
+
+                    # Apply energy threshold
+                    # electron_energy = N_0 * e * g_bse_value
+                    # print('electron_energy: ', electron_energy)
+
+                    # Calculate the initial energy of the primary electrons (E_0)
+                    V_acc = 60e3  # Accelerating voltage in Volts
+                    E_0 = V_acc * e  # Initial energy of primary electrons in joules
+
+                    # Calculate the energy of the backscatter electrons
+                    # E_BSE = E_0 * eta(theta_local, A)
+                    # Compute a weighted average BSE coefficient for the alloy
+                    E_BSE = E_0 * compute_weighted_eta(theta_local, Ti6Al4V)
+                    print('E_BSE: ', E_BSE)
+
+                    # Bifurcation logic for PE and SE
+                    if E_BSE > threshold_energy:
+                        # Backscatter electrons (BSE)
+                        # Compute scattering function g_bse_value
+                        g_bse_value = g_BSE(theta_local, alpha, beta)
+
+                        # Ensure g_bse_value is non-negative
+                        if g_bse_value < 0:
+                            raise ValueError("g_bse_value is negative, check calculations.")
+
+                        # Compute solid angle (Equation A1)
+                        v1, v2, v3 = get_detector_vertices(det_pos)
+                        print("Vertices:", v1, v2, v3)
+                        # print("Type of vertices:", type(v1), type(v2), type(v3))
+                        dOmega = compute_solid_angle(v1, v2, v3)
+
+                        # Compute number of electrons collected by the detector (Equation A2)
+                        # The beam ray stays on each grid point ij for the dwell time t_d and distributes N_0 PE into it
+                        N_0 = beam_current * dwell_time / e  # Equation 3, Section 2.1.1
+
+                        # There are N1 electrons in the chamber after the 1st scattering event
+                        # N_1 = N_0 * eta(theta_local, A)  # Equation 4, Section 2.1.1
+                        # Compute a weighted average BSE coefficient for the alloy
+                        N_1 = N_0 * compute_weighted_eta(theta_local, Ti6Al4V)  # Equation 4, Section 2.1.1
+
+                        # The number of electrons N_1_d collected by a general detector d
+                        N_1_d = N_1 * np.sum(g_bse_value * dOmega)  # Equation 5, Section 2.1.1
+
+                        # Accumulate electrons per detector
+                        electrons_per_detector[det_idx] += N_1_d
+
+                        # Add to total electrons
+                        total_electrons += N_1_d
+
+                        # Debugging print statements
+                        print(f"Domain ({domain_x}, {domain_y}):")
+                        print(f"n_S_local: {n_S_local}")
+                        print(f"theta_local: {np.rad2deg(theta_local)}, "
+                              f"alpha: {np.rad2deg(alpha)}, "
+                              f"beta: {np.rad2deg(beta)}")
+                        print(f"det_vec: {det_vec}")
+                        print(f"p_E_reflected: {p_E_reflected}")
+                        print(f"g_bse_value: {g_bse_value}, N_0: {N_0}, N_1: {N_1}, N_1_d: {N_1_d}")
+                        print(f"total_electrons: {total_electrons}")
+                    else:
+                        # Secondary electrons (SE)
+                        # Note: This part is currently neglected in the model
+                        # but can be implemented if necessary
+                        continue  # Skip to the next iteration if energy is below threshold
+
+
+print(f"Electrons per detector: {electrons_per_detector}")
+for i, electrons in enumerate(electrons_per_detector):
+    print(f"Electrons hitting detector {i + 1}: {np.sum(electrons)}")
+print(f"Total electrons hitting all detectors: {total_electrons}")