diff --git a/src/ray_tracer/__init__.py b/src/ray_tracer/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..57f07dd6671c593d9a4d36ea57839ff187b0a416
--- /dev/null
+++ b/src/ray_tracer/__init__.py
@@ -0,0 +1,187 @@
+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}")
diff --git a/src/ray_tracer/raytracer.py b/src/ray_tracer/raytracer.py
new file mode 100644
index 0000000000000000000000000000000000000000..014b2832f098808b3b8c939cd30b784de8be7f7e
--- /dev/null
+++ b/src/ray_tracer/raytracer.py
@@ -0,0 +1,575 @@
+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()
diff --git a/src/ray_tracer/raytracerMM.py b/src/ray_tracer/raytracerMM.py
new file mode 100644
index 0000000000000000000000000000000000000000..eced7951d93ffb0b0a2a297e76eb12349fe17570
--- /dev/null
+++ b/src/ray_tracer/raytracerMM.py
@@ -0,0 +1,389 @@
+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}")