import xraylib
from scipy.optimize import root_scalar
import numpy as np
__all__ = """
lookup_diameter
materials_to_deltas
materials_to_linear_attenuation
calc_lookup_table
get_densities
""".split()
# Lookup table where each entry is a tuple (column1, column2)
Lens_diameter_table = [
(50, 450.0),
(100, 632.0),
(200, 894.0),
(300, 1095.0),
(500, 1414.0),
(1000, 2000.0),
(1500, 2450.0),
]
# Convert the lookup table to a dictionary for faster lookup
Lens_diameter_dict = {int(col1): col2 for col1, col2 in Lens_diameter_table}
#Using local density definitions until package/library found with longer list
#in g/cm^3
DENSITY = {'Si': 2.33, 'TiO2': 4.23, 'InSb': 5.78}
def get_densities(materials):
densities = dict.fromkeys(materials)
for material in list(densities):
try:
matdb = xraylib.CompoundParser(material)
except ValueError as err:
print(f"{material} not found in xraylib.")
else:
if matdb['nAtomsAll'] == 1:
density = xraylib.ElementDensity(matdb['Elements'][0])
else:
if material in list(DENSITY):
density = DENSITY[material]
else:
raise ValueError(f"{material} not found in DENSITY keys.")
densities[material]=density
return densities
def index_to_binary_list(index, length):
"""
Converts an index number to its binary representation as a list of digits,
and pads the list with zeros in front to achieve the desired length.
Parameters:
index (int): The index number to be converted.
length (int): The desired length of the binary list.
Returns:
list: A list of digits representing the binary representation of the index.
"""
# Convert the index to a binary string and remove the '0b' prefix
binary_str = bin(index)[2:]
# Pad the binary string with zeros in front to achieve the desired length
#padded_binary_str = binary_str.zfill(length)
# Reverse the binary string
reversed_binary_str = binary_str[::-1]
# Convert the reversed binary string to a list of integers
binary_list = [int(digit) for digit in reversed_binary_str]
# Pad the list with zeros at the end to achieve the desired length
while len(binary_list) < length:
binary_list.append(0)
return binary_list
def binary_list_to_index(binary_list, length):
"""
Converts a list of binary digits in reverse order to its integer representation,
padding the list with zeros at the end to have a fixed number of elements.
Parameters:
binary_list (list): A list of digits representing the binary number in reverse order.
length (int): The fixed number of elements the list should have.
Returns:
int: The integer representation of the binary number.
"""
# Pad the list with zeros at the end to achieve the desired length
while len(binary_list) < length:
binary_list.append(0)
# Convert the binary list to an integer
index = 0
for i, digit in enumerate(binary_list):
index += digit * 2**i
return index
def lookup_diameter(lens_radius):
# Convert the input float to an integer
input_int = int(round(lens_radius*1.0e6))
return Lens_diameter_dict.get(input_int, (lens_radius*1.0e6)**0.5*63.222+ 0.73)/1.0e6
def materials_to_deltas(material_list, energy):
"""
Convert a list of material names to a list of delta values at a given energy.
Parameters:
material_list (list): A list of material names.
energy (float): The energy in keV.
Returns:
list: A list of delta values for the given materials at the given energy.
"""
# The list to store delta values
delta_list = []
# Iterate through each material in the input list
for material in material_list:
# Compute the delta value for the current material at the given energy
Z = xraylib.SymbolToAtomicNumber(material)
density = xraylib.ElementDensity(Z)
delta = 1.0-xraylib.Refractive_Index_Re(material, energy, density)
# Add the delta value to the delta list
delta_list.append(delta)
return np.array(delta_list)
def materials_to_linear_attenuation(material_list, energy):
"""
Convert a list of material names to a list of linear attenuation coefficients at a given energy.
Parameters:
material_list (list): A list of material names.
energy (float): The energy in keV.
Returns:
list: A list of linear attenuation coefficient values (in m^-1) for the given materials at the given energy.
"""
# The list to store linear attenuation coefficient values
mu_list = []
# Iterate through each material in the input list
for material in material_list:
# Compute the delta value for the current material at the given energy
Z = xraylib.SymbolToAtomicNumber(material)
density = xraylib.ElementDensity(Z)
# Compute the mass attenuation coefficient in cm^2/g
#mass_attenuation = xraylib.CS_Photo(Z, energy)
mass_attenuation = xraylib.CS_Total(Z, energy)
# Convert mass attenuation coefficient to linear attenuation coefficient in m^-1
mu = mass_attenuation * density * 100.0
# Add the linear attenuation coefficient value to the list
mu_list.append(mu)
return mu_list
def absorptionaperture(x, n1mud, sigma, n1mur):
'''
Description:
TODO
Parameters:
TODO
Returns:
TODO
'''
numerator = np.exp(-(x**2/(2*sigma**2))) * np.exp(-n1mur*(x**2) - n1mud)
denominator = np.exp(-n1mud)
return numerator / denominator - 0.5
def calc_lookup_table(num_configs, radii, materials, energy_keV, wl, numlens,
lens_loc, beam, bl, crl, slit1_H, slit1_V, thickerr,
flag_HE = False):
lookupTable = np.empty(num_configs)
sigmaH = beam['sigmaH']
sigmaV = beam['sigmaV']
sigmaHp = beam['sigmaHp']
sigmaVp = beam['sigmaVp']
d_StoL1 = bl['d_StoL']
d_StoF = bl['d_Stof']
d_min = crl['d_min']
L1_D = np.zeros(len(radii)) # CRL1 diameters for each stack
for i in range(len(radii)):
L1_D[i] = lookup_diameter(radii[i])
L1_delta = materials_to_deltas(materials, energy_keV) # delta values for CRL1 stacks
L1_mu = materials_to_linear_attenuation(materials, energy_keV) # mu values for CRL1 stacks
L1_Feq = radii/(2*numlens*L1_delta) + lens_loc # CRL1 equivalent f in m for each stack
L1_index_n = 2**L1_Feq.size # Total number of combinations for CRL1
L1_invF_list= np.zeros(L1_index_n) # List of equivalent 1/f in m^-1 for CRL1
for i in range(L1_index_n):
L1_invF_list[i] = np.sum(index_to_binary_list(i, L1_Feq.size)/L1_Feq)
# Sort the L1_invF list (to avoid zigzagging)
L1_invF_list_sort_indices = np.argsort(L1_invF_list)
L1_invF_list_sorted = L1_invF_list[L1_invF_list_sort_indices]
q1_list = 1/(L1_invF_list_sorted - 1/d_StoL1) # focal position of CRL1 for all configurations (sorted)
dq1_list = q1_list - (d_StoF - d_StoL1)
# Start generating focal size list as a function of CRL1 setting
sigma1H = (sigmaH**2 + (sigmaHp*d_StoL1)**2)**0.5 # sigma beam size before CRL1
sigma1V = (sigmaV**2 + (sigmaVp*d_StoL1)**2)**0.5 # sigma beam size before CRL1
L1_n1mud_list = np.zeros(L1_index_n) # List of n1*mu*d_min for all possible CRL1 configurations
L1_n1muR_list = np.zeros(L1_index_n) # List of n1*mu/R for all possible CRL1 configurations
aperL1H_list = np.zeros(L1_index_n) # absorption H aperture of CRL1 for all configurations
aperL1V_list = np.zeros(L1_index_n) # absorption V aperture of CRL1 for all configurations
diameter1_list = np.zeros(L1_index_n) # CRL1 diameter for all possible configurations
FWHM1H_list = np.zeros(L1_index_n) # H focal size at the focus of CRL1
FWHM1V_list = np.zeros(L1_index_n) # V focal size at the focus of CRL1
Strehl_list = np.zeros(L1_index_n) # Strehl ratio based on lens thickness error
for i in range(L1_index_n):
# absorption aperture is a function of CRL absorption/physical aperture, incident beam size, and physical slits
L1_n1mud_list[i] = np.sum(index_to_binary_list(L1_invF_list_sort_indices[i], L1_Feq.size)*np.array(L1_mu*numlens*d_min))
L1_n1muR_list[i] = np.sum(index_to_binary_list(L1_invF_list_sort_indices[i], L1_Feq.size)*np.array(L1_mu*numlens/radii))
solution = root_scalar(absorptionaperture, args=(L1_n1mud_list[i], sigma1H, L1_n1muR_list[i]), bracket=[0.0, 2*sigma1H], method='bisect')
aperL1H_list[i] = solution.root*2.0
solution = root_scalar(absorptionaperture, args=(L1_n1mud_list[i], sigma1V, L1_n1muR_list[i]), bracket=[0.0, 2*sigma1V], method='bisect')
aperL1V_list[i] = solution.root*2.0
mask = (np.array(index_to_binary_list(L1_invF_list_sort_indices[i], L1_Feq.size)) == 1)
if np.all(mask == False):
diameter1_list[i] = np.inf
else:
diameter1_list[i] = np.min(L1_D[mask])
aperL1H_list[i] = min(aperL1H_list[i], diameter1_list[i], slit1_H)
aperL1V_list[i] = min(aperL1V_list[i], diameter1_list[i], slit1_V)
phase_error_tmp = np.linalg.norm(index_to_binary_list(L1_invF_list_sort_indices[i], L1_Feq.size)*np.array(thickerr*L1_delta)*2*np.pi/wl)
Strehl_list[i] = np.exp(-phase_error_tmp**2)
# FWHMbeam size at CRL1 focus
FWHM1H_list = ((0.88*wl*q1_list/aperL1H_list)**2 + (2.355*sigmaH*q1_list/d_StoL1)**2)**0.5
FWHM1V_list = ((0.88*wl*q1_list/aperL1V_list)**2 + (2.355*sigmaV*q1_list/d_StoL1)**2)**0.5
if flag_HE:
FWHM1H_list *= (Strehl_list)**(-0.5)
FWHM1V_list *= (Strehl_list)**(-0.5)
FWHM_list = (FWHM1H_list*FWHM1V_list)**0.5
return FWHM_list