From d0a6dad24ad9b7592d078ca1ca957bb1f6fe6969 Mon Sep 17 00:00:00 2001
From: Rahil Doshi <rahil.doshi@fau.de>
Date: Thu, 12 Sep 2024 14:31:06 +0200
Subject: [PATCH] Update the raytracer and add plots for better visualization
---
src/ray_tracer/__init__.py | 187 -----
src/ray_tracer/raytracerMM.py | 421 +++++++----
src/ray_tracer/raytracer_plots.py | 1101 +++++++++++++++++++++++++++++
3 files changed, 1378 insertions(+), 331 deletions(-)
create mode 100644 src/ray_tracer/raytracer_plots.py
diff --git a/src/ray_tracer/__init__.py b/src/ray_tracer/__init__.py
index 57f07dd66..e69de29bb 100644
--- a/src/ray_tracer/__init__.py
+++ b/src/ray_tracer/__init__.py
@@ -1,187 +0,0 @@
-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/raytracerMM.py b/src/ray_tracer/raytracerMM.py
index eced7951d..7dcb79dfd 100644
--- a/src/ray_tracer/raytracerMM.py
+++ b/src/ray_tracer/raytracerMM.py
@@ -1,5 +1,6 @@
import numpy as np
import sympy as sp
+# import matplotlib.pyplot as plt
from typing import Tuple
from src.materials.alloys.Alloy import Alloy
from src.materials.alloys import Ti6Al4V
@@ -22,25 +23,28 @@ print('Ti6Al4V.composition: ', Ti6Al4V.composition)
# 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
+dx: float = 1000e-6 # Cell size in meters
+dy: float = 1000e-6 # Cell size in meters
+dz: float = 1000e-6 # Cell size in meters
+sigma: float = 3000e-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)
+# Define the work distance (height of the electron beam gun)
+w_d = 12e-3 # Work distance in meters (272 mm)
# 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_min, x_max = -5000e-6, 5000e-6 # -500e-6, 500e-6
+y_min, y_max = -5000e-6, 5000e-6 # -250e-6, 250e-6
+z_min, z_max = 0, 4000e-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
+print(f"Domain dimensions: x={x_num}, y={y_num}, z={z_num}")
domain = np.zeros((z_num, y_num, x_num)) # (z, y , x)
print('domain shape: ', np.shape(domain))
domain[0:2, :, :] = 1 # Solid cells
@@ -84,17 +88,19 @@ 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 = np.arange(-2, 3)
+y = np.arange(-2, 3)
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
+print('beam_profile: ', 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_height = 12e-3 # Detector height in meters (272 mm)
+detector_spacing = 5e-3 # Detector spacing in meters (127 mm)
+detector_diameter = 5e-3 # Detector diameter in meters (50 mm)
+
detector_positions = np.array([
[detector_spacing, 0, detector_height], # Right detector
[-detector_spacing, 0, detector_height], # Left detector
@@ -103,6 +109,33 @@ detector_positions = np.array([
])
+
+
+# Define the origin coordinates
+origin_x = (x_num - 1) // 2 # Assuming x_num = 101, origin_x = 50
+print('origin_x: ',origin_x)
+origin_y = (y_num - 1) // 2 # Assuming y_num = 51, origin_y = 25
+print('origin_y: ',origin_y)
+origin_z = 2 # Assuming the origin is at z=2 (interface)
+
+origin_point = np.array([origin_x * dx, origin_y * dy, origin_z * dz])
+# origin_point = np.array([0, 0, 0])
+print('origin_x * dx: ', origin_x * dx)
+print("origin_point: ", origin_point)
+
+# Loop over detector positions
+for i in range(len(detector_positions)):
+ for j in range(i+1, len(detector_positions)):
+ dist = np.linalg.norm(detector_positions[i] - detector_positions[j])
+ print(f"Distance between detectors {i+1} and {j+1}: {dist}")
+
+ # Calculate distance from detector to origin
+ dist_to_origin = np.linalg.norm(detector_positions[i] - origin_point)
+ print(f"Distance from detector {i+1} to origin ({origin_x}, {origin_y}, {origin_z}): {dist_to_origin}")
+
+
+
+
# Define detector vertices for solid angle calculation
def get_detector_vertices(detector_position: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
@@ -115,11 +148,13 @@ def get_detector_vertices(detector_position: np.ndarray) -> Tuple[np.ndarray, np
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
+ side_length = np.sqrt(3) * (detector_diameter / 2)
+ height = np.sqrt(3) * (side_length / 2)
# Determine the orientation based on the detector position
- if detector_position[0] > 0: # Right detector
+ if np.allclose(detector_position[:2], [0, 0]):
+ raise ValueError("Detector position cannot be at the origin (0, 0, 0)")
+ elif 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])
@@ -131,10 +166,20 @@ def get_detector_vertices(detector_position: np.ndarray) -> Tuple[np.ndarray, np
_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
+ elif detector_position[1] < 0: # 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])
+ else:
+ raise ValueError("Invalid detector position. Must be in one of the four expected positions.")
+
+ print(f"Detector {det_idx+1} vertices:")
+ print(f"v1: {_v1}")
+ print(f"v2: {_v2}")
+ print(f"v3: {_v3}")
+ _normal = np.cross(_v2 - _v1, _v3 - _v1)
+ _normal /= np.linalg.norm(_normal)
+ print(f"Detector {det_idx+1} normal vector: {_normal}")
return _v1, _v2, _v3
@@ -251,139 +296,227 @@ def g_BSE(_theta: float, _alpha: float, _beta: float) -> float:
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)
+# Calculate the beam direction vector p_E
+def calculate_p_E(_beam_x, _beam_y):
+ # Calculate the vector from the beam gun to the beam location
+ beam_vector = np.array([_beam_x * dx, _beam_y * dy, -w_d])
+ # Normalize the beam vector
+ _p_E = beam_vector / np.linalg.norm(beam_vector)
+ return _p_E
+
+
+# Calculate phi as the angle between p_E and the negative z-axis
+def calculate_phi(_p_E):
+ z_axis = np.array([0, 0, -1])
+ dot_product = np.dot(_p_E, z_axis)
+ _phi_rad = np.arccos(np.clip(dot_product, -1.0, 1.0))
+ return _phi_rad
-# 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('x_num: ', x_num)
print('y_num: ', y_num)
-beam_loc = ((x_num+1) // 2, (y_num+1) // 2) # (x, y)
+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
+domain_center = ((x_num-1)/2, (y_num-1)/2)
+print(f"Domain center: {domain_center}")
+print(f"Is beam centered: {beam_loc == domain_center}")'''
+
+
+# Initialize arrays to store angles
+phi_angles = []
+theta_local_angles = []
+alpha_angles = []
+beta_angles = []
+
+# Initialize arrays to store electron counts
+electrons_per_detector_per_step = np.zeros((x_num, y_num, len(detector_positions)))
+total_electrons_per_step = np.zeros((x_num, y_num))
+
+# Calculate the center indices
+x_center = (x_num - 1) // 2 # 5
+y_center = (y_num - 1) // 2 # 5
+print('x_center, y_center: ', x_center, y_center)
+
+# Main simulation loop to iterate over the simulation domain
+# for beam_x in range(x_num//4, (x_num*3)//4+2):
+for beam_x in range(x_center, x_center + 1):
+ for beam_y in range(y_center, y_center + 1):
+ beam_loc = (beam_x, beam_y)
+ print(f"Beam location: {beam_loc}")
+
+ # Calculate the beam direction vector p_E
+ p_E = calculate_p_E(beam_x - x_center, beam_y - y_center)
+ print('p_E: ', p_E)
+ # dx = (beam_x - x_center) * dx
+ # dy = (beam_y - y_center) * dy
+
+ # Calculate phi as the angle between p_E and the negative z-axis
+ # phi_rad = np.arctan2(np.sqrt(dx**2 + dy**2), w_d)
+ phi_rad = calculate_phi(p_E)
+ phi_angles.append(phi_rad) # Store phi angles
+ print(f"phi_rad: {phi_rad}")
+
+ # p_E = np.array([np.sin(phi_rad), 0, -np.cos(phi_rad)])
+ # print('p_E: ', p_E)
+
+ # Reset counters for this beam location
+ total_electrons: int = 0
+ electrons_per_detector = np.zeros(len(detector_positions))
+
+ # 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
+ if not np.allclose(np.linalg.norm(n_S_local), 1.0):
+ raise ValueError("n_S_local is not normalized.")
# 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"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 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)
+ )
+ theta_local_angles.append(theta_local) # Store theta_local angles
+
+ scattering_point = np.array([
+ x_min + domain_x * dx,
+ y_min + domain_y * dy,
+ z_min + interface_layer * dz
+ ])
+ print('scattering_point: ', scattering_point)
+
+ for det_idx, det_pos in enumerate(detector_positions):
+ print(f"Detector {det_idx+1} position: {det_pos}")
+ det_idx: int
+ det_pos: np.ndarray[int]
+
+ # Calculate the vector from the scattering point to the detector
+ det_vec = det_pos - scattering_point
+ det_vec /= np.linalg.norm(det_vec)
+ if not np.allclose(np.linalg.norm(det_vec), 1.0):
+ raise ValueError("det_vec is not normalized.")
+
+ # Calculate the distance from the scattering point to the detector
+ distance = np.linalg.norm(det_pos - scattering_point)
+ print(f"Distance to detector {det_idx + 1} ({domain_x}, {domain_y}): {distance:.9f} m")
+
+ # Debugging print statements
+ print(f"Detector {det_idx+1} position: {det_pos}")
+ print(f"Detector Position Type: {type(det_pos)}")
+ print(f"Scattering Point: {scattering_point}")
+ print(f"Detector Vector (normalized): {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))
+ alpha_angles.append(alpha) # Store alpha angles
+
+ 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))
+ beta_angles.append(beta) # Store beta angles
+
+ # 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)
+ normal = np.cross(v2 - v1, v3 - v1)
+ normal /= np.linalg.norm(normal)
+ print(f"Detector {det_idx+1}:")
+ print(f" Vertices: {v1}, {v2}, {v3}")
+ print(f" Normal vector: {normal}")
+ # print("Type of vertices:", type(v1), type(v2), type(v3))
+ dOmega = compute_solid_angle(v1 - scattering_point,
+ v2 - scattering_point,
+ v3 - scattering_point)
+ print(f"Detector {det_idx+1} solid angle: {dOmega}")
+
+ # 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}")'''
+ # 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}")
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
+ # Store results for this beam location
+ electrons_per_detector_per_step[beam_x, beam_y] = electrons_per_detector
+ total_electrons_per_step[beam_x, beam_y] = total_electrons
+
+ print(f"Electrons per detector at this step: {electrons_per_detector}")
+ print(f"Total electrons at this step: {total_electrons}")
+
+
+# Compute final results
+total_electrons_all_steps = np.sum(total_electrons_per_step)
+total_electrons_per_detector = np.sum(electrons_per_detector_per_step, axis=(0, 1))
+
+# print(f"Beam location: {beam_loc}")
-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}")
+print(f"Total electrons hitting all detectors: {total_electrons_all_steps}")
+for i, electrons in enumerate(total_electrons_per_detector):
+ print(f"Total electrons hitting detector {i + 1}: {electrons}")
diff --git a/src/ray_tracer/raytracer_plots.py b/src/ray_tracer/raytracer_plots.py
new file mode 100644
index 000000000..53d639b7d
--- /dev/null
+++ b/src/ray_tracer/raytracer_plots.py
@@ -0,0 +1,1101 @@
+import numpy as np
+import sympy as sp
+import matplotlib.pyplot as plt
+import matplotlib.colors as colors
+from typing import Tuple
+from src.materials.alloys.Alloy import Alloy
+from src.materials.alloys import Ti6Al4V
+
+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)
+print('A: ', A)
+dx: float = 1000e-6 # Cell size in meters
+dy: float = 1000e-6 # Cell size in meters
+dz: float = 1000e-6 # Cell size in meters
+sigma: float = 3000e-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)
+# Define the work distance (height of the electron beam gun)
+# w_d = 12e-3 # Work distance in meters (272 mm)
+
+# Domain setup (Section 2.1.1)
+# Define the 3D domain for the simulation
+x_min, x_max = -40000e-6, 40000e-6 # -40 mm to 40 mm
+y_min, y_max = -40000e-6, 40000e-6 # -40 mm to 40 mm
+z_min, z_max = -2000e-6, 2000e-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
+print(f"Domain dimensions: x={x_num}, y={y_num}, z={z_num}")
+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(-2, 3)
+y = np.arange(-2, 3)
+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
+print('beam_profile: ', beam_profile)
+
+# Detector setup (Section 2.1.1)
+# Define positions of the detectors
+detector_height = 272e-3 # Detector height in meters (272 mm)
+detector_spacing = 127e-3 # Detector spacing in meters (127 mm)
+detector_diameter = 50e-3 # Detector diameter in meters (50 mm)
+
+# Define detector positions relative to the origin
+detector_positions = np.array([
+ [detector_spacing, 0, detector_height],
+ [-detector_spacing, 0, detector_height],
+ [0, detector_spacing, detector_height],
+ [0, -detector_spacing, detector_height]
+])
+'''detector_positions = np.array([
+ [x_max, y_max / 2, detector_height], # Right detector
+ [0, y_max / 2, detector_height], # Left detector
+ [x_max / 2, y_max, detector_height], # Back detector
+ [x_max / 2, 0, detector_height] # Front detector
+])'''
+print('detector_positions:', detector_positions)
+
+# Define the origin coordinates
+'''origin_x = (x_num - 1) // 2 # Assuming x_num = 101, origin_x = 50
+print('origin_x: ', origin_x)
+origin_y = (y_num - 1) // 2 # Assuming y_num = 51, origin_y = 25
+print('origin_y: ', origin_y)
+origin_z = (z_num - 1) // 2 # Assuming the origin is at z=2 (interface)
+
+origin_point = np.array([origin_x * dx, origin_y * dy, origin_z * dz])
+# origin_point = np.array([0, 0, 0])
+print('origin_x * dx: ', origin_x * dx)
+print("origin_point: ", origin_point)
+
+# Loop over detector positions
+for i in range(len(detector_positions)):
+ for j in range(i + 1, len(detector_positions)):
+ dist = np.linalg.norm(detector_positions[i] - detector_positions[j])
+ print(f"Distance between detectors {i + 1} and {j + 1}: {dist}")
+
+ # Calculate distance from detector to origin
+ dist_to_origin = np.linalg.norm(detector_positions[i] - origin_point)
+ print(f"Distance from detector {i + 1} to origin ({origin_x}, {origin_y}, {origin_z}): {dist_to_origin}")'''
+
+
+# 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) * (side_length / 2)
+
+ # Determine the orientation based on the detector position
+ if np.allclose(detector_position[:2], [0, 0]):
+ raise ValueError("Detector position cannot be at the origin (0, 0, 0)")
+ elif 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])
+ elif detector_position[1] < 0: # 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])
+ else:
+ raise ValueError("Invalid detector position. Must be in one of the four expected positions.")
+
+ '''print(f"Detector {det_idx+1} vertices:")
+ print(f"v1: {_v1}")
+ print(f"v2: {_v2}")
+ print(f"v3: {_v3}")
+ _normal = np.cross(_v2 - _v1, _v3 - _v1)
+ _normal /= np.linalg.norm(_normal)
+ print(f"Detector {det_idx+1} normal vector: {_normal}")'''
+
+ 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
+
+
+# Calculate the beam direction vector p_E
+def calculate_p_E(_beam_x, _beam_y):
+ # Calculate the vector from the beam gun to the beam location
+ beam_vector = np.array([_beam_x * dx, _beam_y * dy, -detector_height])
+ # Normalize the beam vector
+ _p_E = beam_vector / np.linalg.norm(beam_vector)
+ return _p_E
+
+
+# Calculate phi as the angle between p_E and the negative z-axis
+def calculate_phi(_p_E):
+ z_axis = np.array([0, 0, -1])
+ dot_product = np.dot(_p_E, z_axis)
+ _phi_rad = np.arccos(np.clip(dot_product, -1.0, 1.0))
+ return _phi_rad
+
+
+def set_beam_movement(_x_center, _y_center, _x_offset, _y_offset, _movement_type):
+ if _movement_type == 'X':
+ if _x_offset >= 0:
+ _beam_x_start = _x_center + _x_offset
+ _beam_x_end = _x_center - _x_offset - 1
+ else:
+ _beam_x_start = _x_center + _x_offset
+ _beam_x_end = _x_center - _x_offset + 1
+ _beam_y_start = _y_center + _y_offset
+ _beam_y_end = _y_center + _y_offset + 1
+ elif _movement_type == 'Y':
+ _beam_x_start = _x_center + _x_offset
+ _beam_x_end = _x_center + _x_offset + 1
+ if _y_offset >= 0:
+ _beam_y_start = _y_center + _y_offset
+ _beam_y_end = _y_center - _y_offset - 1
+ else:
+ _beam_y_start = _y_center + _y_offset
+ _beam_y_end = _y_center - _y_offset + 1
+ elif _movement_type == 'XY':
+ if _x_offset >= 0:
+ _beam_x_start = _x_center + _x_offset
+ _beam_x_end = _x_center - _x_offset - 1
+ else:
+ _beam_x_start = _x_center + _x_offset
+ _beam_x_end = _x_center - _x_offset + 1
+ if _y_offset >= 0:
+ _beam_y_start = _y_center + _y_offset
+ _beam_y_end = _y_center - _y_offset - 1
+ else:
+ _beam_y_start = _y_center + _y_offset
+ _beam_y_end = _y_center - _y_offset + 1
+ else:
+ raise ValueError("movement_type must be 'X' or 'Y' or 'XY'")
+
+ return _beam_x_start, _beam_x_end, _beam_y_start, _beam_y_end
+
+
+
+'''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)
+domain_center = ((x_num-1)/2, (y_num-1)/2)
+print(f"Domain center: {domain_center}")
+print(f"Is beam centered: {beam_loc == domain_center}")'''
+
+
+# Initialize arrays to store electron counts
+electrons_per_detector_per_step = np.zeros((x_num, y_num, len(detector_positions)))
+total_electrons_per_step = np.zeros((x_num, y_num))
+
+# Additional arrays to capture detailed information
+electrons_per_detector_per_beam_profile = np.zeros(
+ (beam_profile.shape[0], beam_profile.shape[1], len(detector_positions)))
+total_electrons_per_beam_profile = np.zeros((beam_profile.shape[0], beam_profile.shape[1]))
+
+# Calculate the center indices
+x_center = (x_num - 1) // 2
+y_center = (y_num - 1) // 2
+
+x_offset: int = -5
+y_offset: int = 20
+movement_type: str = 'X'
+
+# Beam X-Movement
+# beam_x_start, beam_x_end = x_center - x_offset, x_center + x_offset + 1
+# beam_y_start, beam_y_end = y_center + y_offset, y_center + y_offset + 1
+
+# Beam Y-Movement
+# beam_x_start, beam_x_end = x_center + x_offset, x_center + x_offset + 1
+# beam_y_start, beam_y_end = y_center - y_offset, y_center + y_offset + 1
+
+# For X or Y Movement
+beam_x_start, beam_x_end, beam_y_start, beam_y_end = set_beam_movement(
+ x_center, y_center, x_offset, y_offset, movement_type)
+
+# Determine the step size based on the direction of movement
+step_x = 1 if beam_x_start <= beam_x_end else -1
+step_y = 1 if beam_y_start <= beam_y_end else -1
+
+print(beam_x_start, beam_x_end, beam_y_start, beam_y_end, step_x, step_y)
+
+# exit()
+
+# Main simulation loop to iterate over the simulation domain
+for beam_x in range(beam_x_start, beam_x_end, step_x):
+ for beam_y in range(beam_y_start, beam_y_end, step_y):
+ beam_loc = (beam_x, beam_y)
+ print(f"Beam location: {beam_loc}")
+
+ # Calculate the beam direction vector p_E
+ p_E = calculate_p_E(beam_x - x_center, beam_y - y_center)
+ # print('p_E: ', p_E)
+
+ # Calculate phi as the angle between p_E and the negative z-axis
+ # phi_rad = np.arctan2(np.sqrt(dx**2 + dy**2), w_d)
+ phi_rad = calculate_phi(p_E)
+ # print(f"phi_rad: {phi_rad}")
+
+ # Reset counters for this beam location
+ total_electrons: int = 0
+ electrons_per_detector = np.zeros(len(detector_positions))
+
+ # 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
+ if not np.allclose(np.linalg.norm(n_S_local), 1.0):
+ raise ValueError("n_S_local is not normalized.")
+
+ # 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 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)
+ )
+
+ scattering_point = np.array([
+ x_min + domain_x * dx,
+ y_min + domain_y * dy,
+ z_min + interface_layer * dz
+ ])
+ print('scattering_point: ', scattering_point)
+
+ for det_idx, det_pos in enumerate(detector_positions):
+ print(f"Detector {det_idx + 1} position: {det_pos}")
+ det_idx: int
+ det_pos: np.ndarray[float]
+
+ # Calculate the vector from the scattering point to the detector
+ det_vec = det_pos - scattering_point
+ det_vec /= np.linalg.norm(det_vec)
+ if not np.allclose(np.linalg.norm(det_vec), 1.0):
+ raise ValueError("det_vec is not normalized.")
+
+ # Calculate the distance from the scattering point to the detector
+ distance = np.linalg.norm(det_pos - scattering_point)
+ print(f"Distance to detector {det_idx + 1} ({domain_x}, {domain_y}): {distance:.9f} m")
+
+ # Debugging print statements
+ '''print(f"Detector {det_idx+1} position: {det_pos}")
+ print(f"Detector Position Type: {type(det_pos)}")
+ print(f"Scattering Point: {scattering_point}")
+ print(f"Detector Vector (normalized): {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)
+ normal = np.cross(v2 - v1, v3 - v1)
+ normal /= np.linalg.norm(normal)
+ '''print(f"Detector {det_idx+1}:")
+ print(f" Vertices: {v1}, {v2}, {v3}")
+ print(f" Normal vector: {normal}")'''
+
+ dOmega = compute_solid_angle(v1 - scattering_point, v2 - scattering_point,
+ v3 - scattering_point)
+ print(f"Detector {det_idx + 1} solid angle: {dOmega}")
+
+ # 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 for this beam profile point
+ electrons_per_detector_per_beam_profile[profile_x, profile_y, det_idx] += N_1_d
+
+ # Accumulate electrons per detector for this domain point
+ electrons_per_detector_per_step[beam_x, beam_y, det_idx] += N_1_d
+
+ # Add to total electrons for this beam profile point
+ total_electrons_per_beam_profile[profile_x, profile_y] += N_1_d
+
+ # Add to total electrons for this domain point
+ total_electrons_per_step[beam_x, beam_y] += N_1_d
+
+ # Debugging print statements
+ print(f"Before addition: electrons_per_detector[{det_idx + 1}]: {electrons_per_detector[det_idx]}")
+ print(f"N_0: {N_0}, N_1: {N_1}, g_bse_value: {g_bse_value}, dOmega: {dOmega}")
+ print(f"N_1_d: {N_1_d}")
+
+ # Accumulate electrons per detector
+ electrons_per_detector[det_idx] += N_1_d
+
+ # Add to total electrons
+ total_electrons += N_1_d
+
+ print(f"After addition: electrons_per_detector[{det_idx + 1}]: {electrons_per_detector[det_idx]}")
+
+ # Print detailed information for this beam profile point
+ print(f"Electrons per detector at this beam profile point ({profile_x}, {profile_y}): {electrons_per_detector_per_beam_profile[profile_x, profile_y]}")
+ print(f"Total electrons at this beam profile point ({profile_x}, {profile_y}): {total_electrons_per_beam_profile[profile_x, profile_y]}")
+
+ # Print detailed information for this domain point after iterating over the beam profile
+ print(f"Electrons per detector at this domain point ({beam_x}, {beam_y}): {electrons_per_detector_per_step[beam_x, beam_y]}")
+ print(f"Total electrons at this domain point ({beam_x}, {beam_y}): {total_electrons_per_step[beam_x, beam_y]}")
+
+# Compute final results
+total_electrons_all_steps = np.sum(total_electrons_per_step)
+total_electrons_per_detector = np.sum(electrons_per_detector_per_step, axis=(0, 1))
+
+# Print final results
+print(f"Total electrons hitting all detectors: {total_electrons_all_steps}")
+for i, electrons in enumerate(total_electrons_per_detector):
+ print(f"Total electrons hitting detector {i + 1}: {electrons}")
+
+for beam_x in range(beam_x_start, beam_x_end, step_x):
+ for beam_y in range(beam_y_start, beam_y_end, step_y):
+ print(f"Electrons per detector at this domain point ({beam_x}, {beam_y}): {electrons_per_detector_per_step[beam_x, beam_y]}")
+
+
+
+
+# Plot results
+# Define the domain points
+domain_points_x = np.arange(beam_x_start, beam_x_end, step_x)
+domain_points_y = np.arange(beam_y_start, beam_y_end, step_y)
+
+# Extract electron counts for each detector based on movement type
+if movement_type == 'X':
+ detector_1_electrons_x = electrons_per_detector_per_step[beam_x_start:beam_x_end:step_x, y_center+y_offset, 0].flatten()
+ detector_2_electrons_x = electrons_per_detector_per_step[beam_x_start:beam_x_end:step_x, y_center+y_offset, 1].flatten()
+ detector_3_electrons_x = electrons_per_detector_per_step[beam_x_start:beam_x_end:step_x, y_center+y_offset, 2].flatten()
+ detector_4_electrons_x = electrons_per_detector_per_step[beam_x_start:beam_x_end:step_x, y_center+y_offset, 3].flatten()
+
+ # Plot for detectors 1 and 2
+ plt.figure(figsize=(10, 6))
+ plt.plot(domain_points_x, detector_1_electrons_x, label='Detector 1', marker='o')
+ plt.plot(domain_points_x, detector_2_electrons_x, label='Detector 2', marker='s')
+ # Add annotations
+ for x, y in zip(domain_points_x, detector_1_electrons_x):
+ plt.annotate(f'{y:.2e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+ for x, y in zip(domain_points_x, detector_2_electrons_x):
+ plt.annotate(f'{y:.2e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+ plt.title("Electrons Captured by Detectors 1 and 2 (X-axis Movement)")
+ plt.xlabel("Domain Point (X-axis)")
+ plt.ylabel("Electrons Captured")
+ plt.legend()
+ plt.grid(True)
+ plt.tight_layout()
+ plt.savefig("detectors_1_2_x_axis.png", dpi=300)
+ plt.close()
+
+ # Plot for detectors 3 and 4
+ plt.figure(figsize=(10, 6))
+ plt.plot(domain_points_x, detector_3_electrons_x, label='Detector 3', marker='o')
+ plt.plot(domain_points_x, detector_4_electrons_x, label='Detector 4', marker='s')
+ # Add annotations
+ for x, y in zip(domain_points_x, detector_3_electrons_x):
+ plt.annotate(f'{y:.8e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+ for x, y in zip(domain_points_x, detector_4_electrons_x):
+ plt.annotate(f'{y:.8e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+ plt.title("Electrons Captured by Detectors 3 and 4 (X-axis Movement)")
+ plt.xlabel("Domain Point (X-axis)")
+ plt.ylabel("Electrons Captured")
+ plt.legend()
+ plt.grid(True)
+ plt.tight_layout()
+ plt.savefig("detectors_3_4_x_axis.png", dpi=300)
+ plt.close()
+
+ # Normalize the data by subtracting the minimum value
+ detector_3_normalized = detector_3_electrons_x - np.min(detector_3_electrons_x)
+ detector_4_normalized = detector_4_electrons_x - np.min(detector_4_electrons_x)
+
+ plt.figure(figsize=(12, 8)) # Increase figure size for better visibility
+
+ plt.plot(domain_points_x, detector_3_normalized, label='Detector 3', marker='o')
+ plt.plot(domain_points_x, detector_4_normalized, label='Detector 4', marker='s')
+
+ # Add annotations for Detector 3
+ for x, y in zip(domain_points_x, detector_3_normalized):
+ plt.annotate(f'{y:.2f}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+ # Add annotations for Detector 4
+ for x, y in zip(domain_points_x, detector_4_normalized):
+ plt.annotate(f'{y:.2f}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+
+ plt.title("Normalized Electrons Captured by Detectors 3 and 4 (X-axis Movement)")
+ plt.xlabel("Domain Point (X-axis)")
+ plt.ylabel("Electrons Captured (Normalized)")
+ plt.legend()
+ plt.grid(True)
+ plt.tight_layout()
+ plt.savefig("detectors_3_4_x_axis_normalized.png", dpi=300)
+ plt.close()
+
+elif movement_type == 'Y':
+ detector_1_electrons_y = electrons_per_detector_per_step[x_center+x_offset, beam_y_start:beam_y_end:step_y, 0].flatten()
+ detector_2_electrons_y = electrons_per_detector_per_step[x_center+x_offset, beam_y_start:beam_y_end:step_y, 1].flatten()
+ detector_3_electrons_y = electrons_per_detector_per_step[x_center+x_offset, beam_y_start:beam_y_end:step_y, 2].flatten()
+ detector_4_electrons_y = electrons_per_detector_per_step[x_center+x_offset, beam_y_start:beam_y_end:step_y, 3].flatten()
+
+ # Plot for detectors 1 and 2
+ plt.figure(figsize=(10, 6))
+ plt.plot(domain_points_y, detector_1_electrons_y, label='Detector 1', marker='o')
+ plt.plot(domain_points_y, detector_2_electrons_y, label='Detector 2', marker='s')
+ # Add annotations
+ for x, y in zip(domain_points_y, detector_1_electrons_y):
+ plt.annotate(f'{y:.8e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+ for x, y in zip(domain_points_y, detector_2_electrons_y):
+ plt.annotate(f'{y:.8e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+ plt.title("Electrons Captured by Detectors 1 and 2 (Y-axis Movement)")
+ plt.xlabel("Domain Point (Y-axis)")
+ plt.ylabel("Electrons Captured")
+ plt.legend()
+ plt.grid(True)
+ plt.tight_layout()
+ plt.savefig("detectors_1_2_y_axis.png", dpi=300)
+ plt.close()
+
+ # Normalize the data by subtracting the minimum value
+ detector_1_normalized = detector_1_electrons_y - np.min(detector_1_electrons_y)
+ detector_2_normalized = detector_2_electrons_y - np.min(detector_2_electrons_y)
+
+ plt.figure(figsize=(12, 8)) # Increase figure size for better visibility
+
+ plt.plot(domain_points_y, detector_1_normalized, label='Detector 1', marker='o')
+ plt.plot(domain_points_y, detector_2_normalized, label='Detector 2', marker='s')
+
+ # Add annotations for Detector 1
+ for x, y in zip(domain_points_y, detector_1_normalized):
+ plt.annotate(f'{y:.2f}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+ # Add annotations for Detector 2
+ for x, y in zip(domain_points_y, detector_2_normalized):
+ plt.annotate(f'{y:.2f}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+
+ plt.title("Normalized Electrons Captured by Detectors 1 and 2 (Y-axis Movement)")
+ plt.xlabel("Domain Point (Y-axis)")
+ plt.ylabel("Electrons Captured (Normalized)")
+ plt.legend()
+ plt.grid(True)
+ plt.tight_layout()
+ plt.savefig("detectors_1_2_y_axis_normalized.png", dpi=300)
+ plt.close()
+
+ # Plot for detectors 3 and 4
+ plt.figure(figsize=(10, 6))
+ plt.plot(domain_points_y, detector_3_electrons_y, label='Detector 3', marker='o')
+ plt.plot(domain_points_y, detector_4_electrons_y, label='Detector 4', marker='s')
+ # Add annotations
+ for x, y in zip(domain_points_y, detector_3_electrons_y):
+ plt.annotate(f'{y:.8e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+ for x, y in zip(domain_points_y, detector_4_electrons_y):
+ plt.annotate(f'{y:.8e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+ plt.title("Electrons Captured by Detectors 3 and 4 (Y-axis Movement)")
+ plt.xlabel("Domain Point (Y-axis)")
+ plt.ylabel("Electrons Captured")
+ plt.legend()
+ plt.grid(True)
+ plt.tight_layout()
+ plt.savefig("detectors_3_4_y_axis.png", dpi=300)
+ plt.close()
+
+elif movement_type == 'XY':
+ # Extract 2D data for XY movement
+ detector_1_electrons_xy = electrons_per_detector_per_step[beam_x_start:beam_x_end:step_x, beam_y_start:beam_y_end:step_y, 0]
+ detector_2_electrons_xy = electrons_per_detector_per_step[beam_x_start:beam_x_end:step_x, beam_y_start:beam_y_end:step_y, 1]
+ detector_3_electrons_xy = electrons_per_detector_per_step[beam_x_start:beam_x_end:step_x, beam_y_start:beam_y_end:step_y, 2]
+ detector_4_electrons_xy = electrons_per_detector_per_step[beam_x_start:beam_x_end:step_x, beam_y_start:beam_y_end:step_y, 3]
+
+ # Create meshgrid for x and y coordinates
+ X, Y = np.meshgrid(np.arange(beam_x_start, beam_x_end, step_x), np.arange(beam_y_start, beam_y_end, step_y))
+
+ # Plot for detectors 1 and 2 (XY movement)
+ '''plt.figure(figsize=(12, 8))
+ plt.pcolormesh(X, Y, detector_1_electrons_xy.T, shading='auto', cmap='viridis')
+ plt.colorbar(label='Electrons Captured (Detector 1)')
+ plt.title("Electrons Captured by Detector 1 (XY Movement)")
+ plt.xlabel("Domain Point (X-axis)")
+ plt.ylabel("Domain Point (Y-axis)")
+ plt.tight_layout()
+ plt.savefig("detector_1_xy_movement.png", dpi=300)
+ plt.close()
+
+ plt.figure(figsize=(12, 8))
+ plt.pcolormesh(X, Y, detector_2_electrons_xy.T, shading='auto', cmap='viridis')
+ plt.colorbar(label='Electrons Captured (Detector 2)')
+ plt.title("Electrons Captured by Detector 2 (XY Movement)")
+ plt.xlabel("Domain Point (X-axis)")
+ plt.ylabel("Domain Point (Y-axis)")
+ plt.tight_layout()
+ plt.savefig("detector_2_xy_movement.png", dpi=300)
+ plt.close()
+
+ # Plot for detectors 3 and 4 (XY movement)
+ plt.figure(figsize=(12, 8))
+ plt.pcolormesh(X, Y, detector_3_electrons_xy.T, shading='auto', cmap='viridis')
+ plt.colorbar(label='Electrons Captured (Detector 3)')
+ plt.title("Electrons Captured by Detector 3 (XY Movement)")
+ plt.xlabel("Domain Point (X-axis)")
+ plt.ylabel("Domain Point (Y-axis)")
+ plt.tight_layout()
+ plt.savefig("detector_3_xy_movement.png", dpi=300)
+ plt.close()
+
+ plt.figure(figsize=(12, 8))
+ plt.pcolormesh(X, Y, detector_4_electrons_xy.T, shading='auto', cmap='viridis')
+ plt.colorbar(label='Electrons Captured (Detector 4)')
+ plt.title("Electrons Captured by Detector 4 (XY Movement)")
+ plt.xlabel("Domain Point (X-axis)")
+ plt.ylabel("Domain Point (Y-axis)")
+ plt.tight_layout()
+ plt.savefig("detector_4_xy_movement.png", dpi=300)
+ plt.close()'''
+
+
+
+
+
+ # Plot for each detector
+ '''fig, axs = plt.subplots(2, 2, figsize=(12, 10))
+
+ # Detector 1
+ vmin = np.min(detector_1_electrons_xy)
+ vmax = np.max(detector_1_electrons_xy)
+ axs[0, 0].pcolormesh(X, Y, detector_1_electrons_xy, shading='auto', cmap='viridis', vmin=vmin, vmax=vmax, norm=None)
+ axs[0, 0].set_title("Electrons Captured by Detector 1")
+ axs[0, 0].set_xlabel("Domain Point (X-axis)")
+ axs[0, 0].set_ylabel("Domain Point (Y-axis)")
+ # plt.colorbar(label='Electrons Captured')
+ axs[0, 0].set_aspect('equal')
+
+ # Detector 2
+ axs[0, 1].pcolormesh(X, Y, detector_2_electrons_xy, shading='auto', cmap='viridis')
+ axs[0, 1].set_title("Electrons Captured by Detector 2")
+ axs[0, 1].set_xlabel("Domain Point (X-axis)")
+ axs[0, 1].set_ylabel("Domain Point (Y-axis)")
+ axs[0, 1].set_aspect('equal')
+
+ # Detector 3
+ axs[1, 0].pcolormesh(X, Y, detector_3_electrons_xy, shading='auto', cmap='viridis')
+ axs[1, 0].set_title("Electrons Captured by Detector 3")
+ axs[1, 0].set_xlabel("Domain Point (X-axis)")
+ axs[1, 0].set_ylabel("Domain Point (Y-axis)")
+ axs[1, 0].set_aspect('equal')
+
+ # Detector 4
+ axs[1, 1].pcolormesh(X, Y, detector_4_electrons_xy, shading='auto', cmap='viridis')
+ axs[1, 1].set_title("Electrons Captured by Detector 4")
+ axs[1, 1].set_xlabel("Domain Point (X-axis)")
+ axs[1, 1].set_ylabel("Domain Point (Y-axis)")
+ axs[1, 1].set_aspect('equal')
+
+ # Add colorbars
+ for ax in axs.flat:
+ fig.colorbar(plt.cm.ScalarMappable(cmap='viridis'), ax=ax, shrink=0.6)
+
+ # Layout so plots do not overlap
+ fig.tight_layout()
+
+ # Save the figure
+ plt.savefig("electron_heatmaps.png", dpi=300)
+ plt.close()'''
+
+
+
+
+ # Create a figure with multiple subplots
+ fig, axs = plt.subplots(2, 2, figsize=(16, 16))
+
+ # Plot for each detector
+ for i in range(4):
+ row = i // 2
+ col = i % 2
+ im = axs[row, col].pcolormesh(X, Y, electrons_per_detector_per_step[beam_x_start:beam_x_end:step_x, beam_y_start:beam_y_end:step_y, i].T, shading='auto', cmap='viridis', norm=colors.LogNorm())
+ axs[row, col].set_title(f"Electrons Captured by Detector {i+1}")
+ axs[row, col].set_xlabel("Domain Point (X-axis)")
+ axs[row, col].set_ylabel("Domain Point (Y-axis)")
+ fig.colorbar(im, ax=axs[row, col])
+
+ # Adjust layout and save
+ plt.tight_layout()
+ plt.savefig("_all_detectors_xy_movement_.png", dpi=300)
+ plt.close()
+
+ # Now let's create individual plots for each detector
+ for i in range(4):
+ plt.figure(figsize=(10, 8))
+ im = plt.pcolormesh(X, Y, electrons_per_detector_per_step[beam_x_start:beam_x_end:step_x, beam_y_start:beam_y_end:step_y, i].T, shading='auto', cmap='viridis', norm=colors.LogNorm())
+ plt.colorbar(im)
+ plt.title(f"Electrons Captured by Detector {i+1}")
+ plt.xlabel("Domain Point (X-axis)")
+ plt.ylabel("Domain Point (Y-axis)")
+ plt.savefig(f"_detector_{i+1}_xy_movement_.png", dpi=300)
+ plt.close()
+
+
+
+
+ detector_1_electrons_xy = electrons_per_detector_per_step[beam_x_start:beam_x_end, beam_y_start:beam_y_end, 0]
+ print(detector_1_electrons_xy)
+
+ # Create the plot
+ plt.figure(figsize=(10, 8))
+ im = plt.pcolormesh(X, Y, detector_1_electrons_xy, shading='auto', cmap='viridis', norm=colors.LogNorm())
+ plt.colorbar(im, label='Electrons Captured')
+ plt.title("Electrons Captured by Detector 1 (XY Movement)")
+ plt.xlabel("Domain Point (X-axis)")
+ plt.ylabel("Domain Point (Y-axis)")
+
+ # Set the correct extent for the plot
+ # plt.xlim(beam_x_start, beam_x_end - 1)
+ # plt.ylim(beam_y_start, beam_y_end - 1)
+
+ # Adjust color scale to highlight variations
+ '''vmin = np.min(detector_1_electrons_xy)
+ vmax = np.max(detector_1_electrons_xy)
+ plt.clim(float(vmin), float(vmax))'''
+
+ # Add text annotations for each cell
+ for i in range(detector_1_electrons_xy.shape[0]):
+ for j in range(detector_1_electrons_xy.shape[1]):
+ plt.text(beam_x_start + i, beam_y_start + j, f'{detector_1_electrons_xy[i, j]:.2f}',
+ ha='center', va='center', fontsize=6, color='white')
+
+ plt.tight_layout()
+ plt.savefig("detector_1_xy_movement_corrected.png", dpi=300)
+ plt.close()
+
+ print(f"Min: {np.min(detector_1_electrons_xy)}, Max: {np.max(detector_1_electrons_xy)}")
+
+# --------------------------------------------------
+'''
+# Plot for detectors 1 and 2 (X-axis movement)
+plt.figure(figsize=(10, 6))
+plt.plot(domain_points_x, detector_1_electrons_x, label='Detector 1', marker='o')
+plt.plot(domain_points_x, detector_2_electrons_x, label='Detector 2', marker='s')
+# Add annotations for Detector 1
+for x, y in zip(domain_points_x, detector_1_electrons_x):
+ plt.annotate(f'{y:.8e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+# Add annotations for Detector 2
+for x, y in zip(domain_points_x, detector_2_electrons_x):
+ plt.annotate(f'{y:.8e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+plt.title("Electrons Captured by Detectors 1 and 2 (X-axis Movement)")
+plt.xlabel("Domain Point (X-axis)")
+plt.ylabel("Electrons Captured")
+plt.legend()
+plt.grid(True)
+plt.tight_layout()
+plt.savefig("detectors_1_2_x_axis.png", dpi=300)
+plt.close()
+
+# --------------------------------------------------
+
+# Plot for detectors 3 and 4 (X-axis movement)
+plt.figure(figsize=(10, 6))
+plt.plot(domain_points_x, detector_3_electrons_x, label='Detector 3', marker='o')
+plt.plot(domain_points_x, detector_4_electrons_x, label='Detector 4', marker='s')
+
+# Add annotations for Detector 3
+for x, y in zip(domain_points_x, detector_3_electrons_x):
+ plt.annotate(f'{y:.8e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+
+# Add annotations for Detector 4
+for x, y in zip(domain_points_x, detector_4_electrons_x):
+ plt.annotate(f'{y:.8e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+
+# Calculate y-axis limits
+y_min = min(min(detector_3_electrons_x), min(detector_4_electrons_x))
+y_max = max(max(detector_3_electrons_x), max(detector_4_electrons_x))
+y_range = y_max - y_min
+
+# Set y-axis limits to focus on the region of interest
+plt.ylim(y_min - 0.1*y_range, y_max + 0.1*y_range)
+
+# Use scientific notation for y-axis
+plt.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
+
+plt.title("Electrons Captured by Detectors 3 and 4 (X-axis Movement)")
+plt.xlabel("Domain Point (X-axis)")
+plt.ylabel("Electrons Captured")
+plt.legend()
+plt.grid(True)
+plt.tight_layout()
+plt.savefig("detectors_3_4_x_axis.png", dpi=300)
+plt.close()
+
+
+# Normalize the data by subtracting the minimum value
+detector_3_normalized = detector_3_electrons_x - np.min(detector_3_electrons_x)
+detector_4_normalized = detector_4_electrons_x - np.min(detector_4_electrons_x)
+
+plt.figure(figsize=(12, 8)) # Increase figure size for better visibility
+
+plt.plot(domain_points_x, detector_3_normalized, label='Detector 3', marker='o')
+plt.plot(domain_points_x, detector_4_normalized, label='Detector 4', marker='s')
+
+# Add annotations for Detector 3
+for x, y in zip(domain_points_x, detector_3_normalized):
+ plt.annotate(f'{y:.2f}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+
+# Add annotations for Detector 4
+for x, y in zip(domain_points_x, detector_4_normalized):
+ plt.annotate(f'{y:.2f}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+
+plt.title("Normalized Electrons Captured by Detectors 3 and 4 (X-axis Movement)")
+plt.xlabel("Domain Point (X-axis)")
+plt.ylabel("Electrons Captured (Normalized)")
+plt.legend()
+plt.grid(True)
+plt.tight_layout()
+plt.savefig("detectors_3_4_x_axis_normalized.png", dpi=300)
+plt.close()
+
+# --------------------------------------------------
+
+# Plot for detectors 1 and 2 (Y-axis movement)
+plt.figure(figsize=(10, 6))
+plt.plot(domain_points_y, detector_1_electrons_y, label='Detector 1', marker='o')
+plt.plot(domain_points_y, detector_2_electrons_y, label='Detector 2', marker='s')
+
+# Add annotations for Detector 1
+for x, y in zip(domain_points_y, detector_1_electrons_y):
+ plt.annotate(f'{y:.8e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+
+# Add annotations for Detector 2
+for x, y in zip(domain_points_y, detector_2_electrons_y):
+ plt.annotate(f'{y:.8e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+
+plt.title("Electrons Captured by Detectors 1 and 2 (Y-axis Movement)")
+plt.xlabel("Domain Point (Y-axis)")
+plt.ylabel("Electrons Captured")
+plt.legend()
+plt.grid(True)
+plt.tight_layout()
+plt.savefig("detectors_1_2_y_axis.png", dpi=300)
+plt.close()
+
+
+# Normalize the data by subtracting the minimum value
+detector_1_normalized = detector_1_electrons_y - np.min(detector_1_electrons_y)
+detector_2_normalized = detector_2_electrons_y - np.min(detector_2_electrons_y)
+
+plt.figure(figsize=(12, 8)) # Increase figure size for better visibility
+
+plt.plot(domain_points_y, detector_1_normalized, label='Detector 1', marker='o')
+plt.plot(domain_points_y, detector_2_normalized, label='Detector 2', marker='s')
+
+# Add annotations for Detector 1
+for x, y in zip(domain_points_y, detector_1_normalized):
+ plt.annotate(f'{y:.2f}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+
+# Add annotations for Detector 2
+for x, y in zip(domain_points_y, detector_2_normalized):
+ plt.annotate(f'{y:.2f}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+
+plt.title("Normalized Electrons Captured by Detectors 1 and 2 (Y-axis Movement)")
+plt.xlabel("Domain Point (Y-axis)")
+plt.ylabel("Electrons Captured (Normalized)")
+plt.legend()
+plt.grid(True)
+plt.tight_layout()
+plt.savefig("detectors_1_2_y_axis_normalized.png", dpi=300)
+plt.close()
+
+# --------------------------------------------------
+
+# Plot for detectors 3 and 4 (Y-axis movement)
+plt.figure(figsize=(10, 6))
+plt.plot(domain_points_y, detector_3_electrons_y, label='Detector 3', marker='o')
+plt.plot(domain_points_y, detector_4_electrons_y, label='Detector 4', marker='s')
+
+# Add annotations for Detector 3
+for x, y in zip(domain_points_y, detector_3_electrons_y):
+ plt.annotate(f'{y:.8e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+
+# Add annotations for Detector 4
+for x, y in zip(domain_points_y, detector_4_electrons_y):
+ plt.annotate(f'{y:.8e}', (x, y), textcoords="offset points", xytext=(5,-15), ha='center', fontsize=8, rotation=45)
+
+plt.title("Electrons Captured by Detectors 3 and 4 (Y-axis Movement)")
+plt.xlabel("Domain Point (Y-axis)")
+plt.ylabel("Electrons Captured")
+plt.legend()
+plt.grid(True)
+plt.tight_layout()
+plt.savefig("detectors_3_4_y_axis.png", dpi=300)
+plt.close()
+'''
+# --------------------------------------------------
+
+'''
+# Create meshgrid for x and y coordinates
+X, Y = np.meshgrid(np.arange(beam_x_start, beam_x_end), np.arange(beam_y_start, beam_y_end))
+
+# Extract electron counts for each detector
+detector_1_electrons = electrons_per_detector_per_step[beam_x_start:beam_x_end, beam_y_start:beam_y_end, 0]
+detector_2_electrons = electrons_per_detector_per_step[beam_x_start:beam_x_end, beam_y_start:beam_y_end, 1]
+detector_3_electrons = electrons_per_detector_per_step[beam_x_start:beam_x_end, beam_y_start:beam_y_end, 2]
+detector_4_electrons = electrons_per_detector_per_step[beam_x_start:beam_x_end, beam_y_start:beam_y_end, 3]
+
+# Plot for detectors 1 and 2
+plt.figure(figsize=(12, 8))
+plt.pcolormesh(X, Y, detector_1_electrons, shading='auto', cmap='viridis')
+plt.colorbar(label='Electrons Captured (Detector 1)')
+plt.title("Electrons Captured by Detector 1")
+plt.xlabel("Domain Point (X-axis)")
+plt.ylabel("Domain Point (Y-axis)")
+plt.tight_layout()
+plt.savefig("detector_1_heatmap.png", dpi=300)
+plt.close()
+
+plt.figure(figsize=(12, 8))
+plt.pcolormesh(X, Y, detector_2_electrons, shading='auto', cmap='viridis')
+plt.colorbar(label='Electrons Captured (Detector 2)')
+plt.title("Electrons Captured by Detector 2")
+plt.xlabel("Domain Point (X-axis)")
+plt.ylabel("Domain Point (Y-axis)")
+plt.tight_layout()
+plt.savefig("detector_2_heatmap.png", dpi=300)
+plt.close()
+
+# Plot for detectors 3 and 4
+plt.figure(figsize=(12, 8))
+plt.pcolormesh(X, Y, detector_3_electrons, shading='auto', cmap='viridis')
+plt.colorbar(label='Electrons Captured (Detector 3)')
+plt.title("Electrons Captured by Detector 3")
+plt.xlabel("Domain Point (X-axis)")
+plt.ylabel("Domain Point (Y-axis)")
+plt.tight_layout()
+plt.savefig("detector_3_heatmap.png", dpi=300)
+plt.close()
+
+plt.figure(figsize=(12, 8))
+plt.pcolormesh(X, Y, detector_4_electrons, shading='auto', cmap='viridis')
+plt.colorbar(label='Electrons Captured (Detector 4)')
+plt.title("Electrons Captured by Detector 4")
+plt.xlabel("Domain Point (X-axis)")
+plt.ylabel("Domain Point (Y-axis)")
+plt.tight_layout()
+plt.savefig("detector_4_heatmap.png", dpi=300)
+plt.close()
+'''
--
GitLab