Refactoring XS code

100idPyApp/Db/pyDevCRLsingle.db

@@ -14,7 +14,7 @@ record(ao, "$(P)$(CRL)$(N):min_thick"){
 
 
 record(ao, "$(P)$(CRL)$(N):energy"){
-	field(EGU, "eV")
+	field(EGU, "keV")
 	field(DTYP, "pydev")
 	field(OUT, "@$(OBJ).updateE('VAL')")
 }

100idPyApp/Db/pyDevFilters.db

@@ -72,6 +72,7 @@ record(ai, "$(P)Filter$(N):energy_RBV"){
   field(DTYP, "pydev")
   field(INP,  "@pydev.iointr('new_energy')")
   field(SCAN, "I/O Intr")
+  field(EGU, "keV")
   field(PREC, "4")
 
 }

100idPyApp/python/transfocator_calcs.py

+import xraylib
+
+__all__ = """
+	lookup_diameter 
+	materials_to_deltas 
+	materials_to_linear_attenuation 
+	calc_lookup_table
+""".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}
+
+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 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 calc_lookup_table(num_configs, radii, materials, energy_keV, numlens, lens_loc):
+	lookupTable = np.empty(num_configs)
+	
+	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(material, energy_keV) # mu values for CRL1 stacks
+	L1_Feq      = radii/(2*numlens*L1_delta) + lens_loc                       # CRL1 equivalent f in m for each stack
+
+	#------------------> Needs refactoring <--------------------------------
+
+	
+	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*L1_n*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*L1_n/L1_R))
+		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(L1_HE*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
+	#------------------> End refactoring <--------------------------------
+
+	return FWHM_list

iocBoot/ioc100idPyCRL/examples/pydeviceCRL.cmd

@@ -6,27 +6,45 @@
 epicsEnvSet("MONOE","testMonoE")
 dbLoadRecords("${TOP}/db/energyTestTools.db","P=$(PREFIX), MONOE=$(MONOE), IDENERGY=testIDE")
 
-#CRL 1 (substitutions file NEEDS WORK)
+# CRL DBs and defining substitution file to get stack properties
 dbLoadTemplate("substitutions/pyDeviceTfor.substitutions","P=$(PREFIX)")
 pydev("stack_subFile = 'substitutions/pyDeviceTfor.substitutions'")
 
-# python code NEEDS WORK
+# Import Transfocator class
 pydev("from pyTransfocator_single import singleTF")
-pydev("CRL1 = singleTF()")
 
-# Update the next line for the beamline energy
-epicsEnvSet("BLE","$(PREFIX)$(MONOE)")
-dbLoadRecords(""${TOP}/db/pyDevCRLsingle.db","P=$(PREFIX), N=1, OBJ=CRL1, STACK_THICK=0.050, MIN_THICK=0.003")
+# Set up beam properties -- these may only need to be set up at creation of IOC
+# energy	 - energy in keV (can be updated after boot-time)
+# L_und      - Undulator length in m
+# sigmaH_e   - Sigma electron source size in H direction in m
+# sigmaV_e   - Sigma electron source size in V direction in m
+# sigmaHp_e  - Sigma electron divergence in H direction in rad
+# sigmaVp_e  - Sigma electron divergence in V direction in rad
+pydev("b_properties = {'energy': 15, 'L_und':4.7, 'sigmaH_e':14.8e-6, 'sigmaV_e':3.7e-6, 'sigmaHp_e':2.8e-6, 'sigmaVp_e':1.5e-6})
 
+# Set up beamline properties -- these may only need to be set up at creation of IOC
+#
+#
+pydev("bl_properties = {'':,'':,'':}
+
+# Set up crl properties -- these are properties common to the stacks
+# 
+#
+#
+pydev("crl_properties = {'':,'':,'':}
+
+# Create Transfocator object
+pydev("CRL1 = singleTF(beam = b_properties, bl = bl_properties, crl = crl_properties)")
 
 # Update the next line for the beamline energy
 epicsEnvSet("BLE","$(PREFIX)$(MONOE)")
-dbLoadRecords("${TOP}/db/pyDevFilters.db","P=$(PREFIX), N=1, OBJ=filters, NUM=12, KEV=$(BLE)")
+# DB for Transfocator system
+dbLoadRecords(""${TOP}/db/pyDevCRLsingle.db","P=$(PREFIX), N=1, OBJ=CRL1, STACK_THICK=0.050, MIN_THICK=0.003")
 
 # After iocInit, need to setup attenuation lookup table
-# 	filters_subFile - string holding name/rel. location of substitutions file loaded
+# 	stack_subFile - string holding name/rel. location of substitutions file loaded
 #                     loaded in  dbLoadTemplate
-#   12  - number of filters in this case
+#   12  - number of lenses in this case
 # 	8.5 - beam energy in keV used for the initial lookup table setup -- can be updated
 #          after startup
-doAfterIocInit("pydev('filters.setupLookupTable(filters_subFile, 12, energy = 8.5)')")
+doAfterIocInit("pydev('CRL1.setupLookupTable(stack_subFile, 12, energy = 8.5)')")