pyTransfocator_single.py

import numpy as np
from scipy.optimize import root_scalar
import xraylib

"""
pyDevice TO DO: 

WHAT inputs change the focal size arrays? Energy, what else?
WHAT inputs change the search through the arrays? desired focal size, what else?

IOC init functions
-get lens stack parameters (# of lenses in each stack, radius, location, thickness, thickness error) -- from substitution file but put into PVs? Update with autosave?
-get source info
	-energy from from ID IOC
	-hor/vert sizes and divergence (also energy dependent)
-lens diameter table? What is it doing?

-desired focal size is changed --> what needs updating? --> nothing, just need to search focal size array again
	-multiple flags: is focal size achievable? is it achievable at sample?

recalc function -- should probably be same as init function
-energy is updated --> what needs updating?

-what else could user/staff change? sample position?
"""

# Beamline input block
energy = 15000.0            # Energy in eV
energy_keV = energy/1000.0  # Energy in keV
wl = 1239.84 / (energy * 10**9)
d_StoL1 = 51.9              # Source-to-CRL1 distance, in m
d_StoL2 = 62.1              # Source-to-CRL2 distance, in m
d_Stof  = 66.2              # Source-to-focus distance, in m
#slit1_H = 500.0e-6          # H slit size before CRL 1
#slit1_V = 300.0e-6          # V slit size before CRL 1

# CRL input block
d_min   = 3.0e-5            # Minimum thickness at the apex in m
stack_d = 50.0e-3           # Stack thickness in m
L1_n    = np.array([1,      1,      1,      1,      1,      1,      2,      4,      8,      16])                # CRL1 number of lenses in each stack
L1_R    = np.array([2.0e-3, 1.0e-3, 5.0e-4, 3.0e-4, 2.0e-4, 1.0e-4, 1.0e-4, 1.0e-4, 1.0e-4, 1.0e-4])            # CRL1 lens radius in each stack
L1_mater= np.array(["Be",   "Be",   "Be",   "Be",   "Be",   "Be",   "Be",   "Be",   "Be",   "Be"])              # CRL1 lens material in each stack
L1_loc  = np.array([4.5,    3.5,    2.5,    1.5,    0.5,    -0.5,   -1.5,   -2.5,   -3.5,   -4.5])*stack_d      # CRL1 lens stack location relative to center stack, positive means upstream
L1_HE   = np.array([1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6, 1.4e-6, 2.0e-6, 2.8e-6, 4.0e-6])            # CRL1 lens RMS thickness error


# Source size input block
L_und = 4.7                 # undulator length
sigmaH_e = 14.8e-6          # Sigma electron source size in H direction in m
sigmaV_e = 3.7e-6           # Sigma electron source size in V direction in m
sigmaHp_e = 2.8e-6          # Sigma electron divergence in H direction in rad
sigmaVp_e = 1.5e-6          # Sigma electron divergence in V direction in rad
sigmaH = (sigmaH_e**2 + wl*L_und/2/np.pi/np.pi)**0.5
sigmaV = (sigmaV_e**2 + wl*L_und/2/np.pi/np.pi)**0.5
sigmaHp = (sigmaHp_e**2 + wl/L_und/2)**0.5
sigmaVp = (sigmaVp_e**2 + wl/L_und/2)**0.5

# 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 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 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 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 absorptionaperture(x, n1mud, sigma, n1mur):
    numerator = np.exp(-(x**2/(2*sigma**2))) * np.exp(-n1mur*(x**2) - n1mud)
    denominator = np.exp(-n1mud)
    return numerator / denominator - 0.5

def find_levels(array, levels, direction='forward'):
    """
    Find the first indices at which the array crosses specified levels and the corresponding crossed values.

    Parameters:
        array (numpy.ndarray): An array of numbers.
        levels (float or numpy.ndarray): A number or an array of levels to find crossings.
        direction (str, optional): The searching direction. Defaults to 'forward'.
                                   Can be either 'forward' or 'backward'.

    Returns:
        tuple: A tuple containing two arrays:
            - An array of first indices at which the array crosses the specified levels.
            - An array of first crossed values at the corresponding indices.
    """

    # Convert a single level to a numpy array
    if isinstance(levels, (int, float)):
        levels = np.array([levels])

    indices = []
    values = []

    # Compute the max and min of the array ignoring NaNs
    max_val = np.nanmax(array)
    min_val = np.nanmin(array)

    for level in levels:
        # If level is out of bounds
        if level > max_val or level < min_val:
            indices.append(-1)
            values.append(np.nan)
            continue

        crossings = []

        if direction == 'forward':
            for i in range(1, len(array)):
                if np.isnan(array[i - 1]) or np.isnan(array[i]):
                    continue
                if (array[i - 1] < level <= array[i]) or (array[i - 1] > level >= array[i]):
                    crossings.append(i - 1)
                    break

        elif direction == 'backward':
            for i in range(len(array) - 2, -1, -1):
                if np.isnan(array[i + 1]) or np.isnan(array[i]):
                    continue
                if (array[i + 1] < level <= array[i]) or (array[i + 1] > level >= array[i]):
                    crossings.append(i)
                    break

        else:
            raise ValueError("Invalid direction. It should be either 'forward' or 'backward'.")

        if len(crossings) > 0:
            idx = crossings[0]
            indices.append(idx)
            values.append(array[idx])
        else:
            # In case no crossing is found within the range
            indices.append(-1)
            values.append(np.nan)

    return np.array(indices), np.array(values)


def Single_CRL2D_control(fsize):

    L1_D        = np.zeros(L1_R.size)                                   # CRL1 diameters for each stack
    for i in range(L1_R.size):
        L1_D[i] = lookup_diameter(L1_R[i])
    L1_delta    = materials_to_deltas(L1_mater, energy_keV)             # delta values for CRL1 stacks
    L1_mu       = materials_to_linear_attenuation(L1_mater, energy_keV) # mu values for CRL1 stacks
    L1_Feq      = L1_R/(2*L1_n*L1_delta) + L1_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*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

    indices, values = find_levels(FWHM_list, fsize, direction='backward')
    index = indices[0]
    if index == -1:
        print(f"Cannot achieve the focal size {fsize*1.0e6:.2f} μm")
    else:
        print("======== Find size at focus ========================================")
        print(f"Energy: {energy_keV} keV")
        print(f"CRL1 configuration index in sorted list is {index}")
        print(f"CRL1 configuration index is {L1_invF_list_sort_indices[index]} or {index_to_binary_list(L1_invF_list_sort_indices[index], L1_Feq.size)}")
        print(f"CRL1 f is {1/L1_invF_list_sorted[index]:.2f} m, focus at q1 = {q1_list[index]:.2f} m")
        print(f"Focal size is {FWHM1H_list[index]*1.0e6:.2f} μm x {FWHM1V_list[index]*1.0e6:.2f} μm at the focal point ({dq1_list[index]*1e3:.1f} mm from sample)")

    FWHM1H_atsample_list = (FWHM1H_list**2 + (aperL1H_list*dq1_list/q1_list)**2)**0.5
    FWHM1V_atsample_list = (FWHM1V_list**2 + (aperL1V_list*dq1_list/q1_list)**2)**0.5
    FWHM_atsample_list   = (FWHM1H_atsample_list*FWHM1V_atsample_list)**0.5
    indices, values = find_levels(FWHM_atsample_list, fsize, direction='forward')
    index2 = indices[0]
    if index2 == -1:
        print(f"Cannot achieve the bame size {fsize*1.0e6:.2f} μm at sample")
    else:
        print("======== Find size at sample =======================================")
        print(f"CRL1 configuration index in sorted list is {index2}")
        print(f"CRL1 configuration index is {L1_invF_list_sort_indices[index2]} or {index_to_binary_list(L1_invF_list_sort_indices[index2], L1_Feq.size)}")
        print(f"CRL1 f is {1/L1_invF_list_sorted[index2]:.2f} m, focus at q1 = {q1_list[index2]:.2f} m ({dq1_list[index2]*1e3:.1f} mm from sample)")
        print(f"Beam size is {FWHM1H_atsample_list[index2]*1.0e6:.2f} μm x {FWHM1V_atsample_list[index2]*1.0e6:.2f} μm at the sample position)")

    indices, values = find_levels(dq1_list, 0.0, direction='backward')
    index3 = indices[0]
    if index == -1:
        print(f"Cannot find combination to focus close to sample")
    else:
        print("======== Find configuration focus close to the sample ==============")
        print(f"CRL1 configuration index in sorted list is {index3}")
        print(f"CRL1 configuration index is {L1_invF_list_sort_indices[index3]} or {index_to_binary_list(L1_invF_list_sort_indices[index3], L1_Feq.size)}")
        print(f"CRL1 f is {1/L1_invF_list_sorted[index3]:.2f} m, focus at q1 = {q1_list[index3]:.2f} m ({dq1_list[index3]*1e3:.1f} mm from sample)")
        print(f"Beam size is {FWHM1H_atsample_list[index3]*1.0e6:.2f} μm x {FWHM1V_atsample_list[index3]*1.0e6:.2f} μm at the sample position)")

    return


if __name__ == "__main__":

    flag_HE = True

    fsize  = 50.0e-6            # Desired focal size in m (area average of h and v size)
    #Single_CRL2D_control(fsize)  # Find the best configuration for a single transfocator system


'''
Update the following to accommodate XS code
'''

class singleTF():
	
	def __init__(self):
		# Initialize beamline layout variables
		self.d_StoL = 51.9              # Source-to-CRL1 distance, in m
		self.d_Stof  = 66.2              # Source-to-focus distance, in m
			
		# Initialize source variables
		self.L_und = 4.7
		self.energy = 15000.0            # Energy in eV
		self.energy_keV = self.energy/1000.0  # Energy in keV
		self.wl = 1239.84 / (self.energy * 10**9)	#Wavelength in nm(?)
		
		self.sigmaH_e = 14.8e-6          # Sigma electron source size in H direction in m
		self.sigmaV_e = 3.7e-6           # Sigma electron source size in V direction in m
		self.sigmaHp_e = 2.8e-6          # Sigma electron divergence in H direction in rad
		self.sigmaVp_e = 1.5e-6          # Sigma electron divergence in V direction in rad
		
		self.sigmaH =  (self.sigmaH_e**2 +  self.wl*self.L_und/2/np.pi/np.pi)**0.5
		self.sigmaV =  (self.sigmaV_e**2 +  self.wl*self.L_und/2/np.pi/np.pi)**0.5
		self.sigmaHp = (self.sigmaHp_e**2 + self.wl/self.L_und/2)**0.5
		self.sigmaVp = (self.sigmaVp_e**2 + self.wl/self.L_und/2)**0.5
		
		# Initialize lens variables
		self.d_min   = 3.0e-5            # Minimum thickness at the apex in m
		self.stack_d = 50.0e-3           # Stack thickness in m
		self.L1_n    = np.array([1,      1,      1,      1,      1,      1,      2,      4,      8,      16])                # CRL1 number of lenses in each stack
		self.L1_R    = np.array([2.0e-3, 1.0e-3, 5.0e-4, 3.0e-4, 2.0e-4, 1.0e-4, 1.0e-4, 1.0e-4, 1.0e-4, 1.0e-4])            # CRL1 lens radius in each stack
		self.L1_mater= np.array(["Be",   "Be",   "Be",   "Be",   "Be",   "Be",   "Be",   "Be",   "Be",   "Be"])              # CRL1 lens material in each stack
		self.L1_loc  = np.array([4.5,    3.5,    2.5,    1.5,    0.5,    -0.5,   -1.5,   -2.5,   -3.5,   -4.5])*stack_d      # CRL1 lens stack location relative to center stack, positive means upstream
		self.L1_HE   = np.array([1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6, 1.4e-6, 2.0e-6, 2.8e-6, 4.0e-6])            # CRL1 lens RMS thickness error

		self.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		
		self.Lens_diameter_dict = {int(col1): col2 for col1, col2 in Lens_diameter_table}
		
		
		# Initialize pre-CRL slit size
		self.slit1_H = 500.0e-6          # H slit size before CRL 1
		self.slit1_V = 300.0e-6          # V slit size before CRL 1

		#Calc lookup table
		self.lookup_table=calc_lookup_table(self, ...)


		self.energy = 0  # gets value from an ao (incoming beam energy)
		self.focalSize = 0 # get value from an ao (desired focal length)
		self.lenses = 0 # sets integer (2^12) whose binary representation indicates which lenses are in or out
		

		
		
	def setSource(self):
		self.L_und = 4.7
		self.energy = 15000.0            # Energy in eV
		self.energy_keV = self.energy/1000.0  # Energy in keV
		self.wl = 1239.84 / (self.energy * 10**9)	#Wavelength in nm(?)
		
		self.sigmaH_e = 14.8e-6          # Sigma electron source size in H direction in m
		self.sigmaV_e = 3.7e-6           # Sigma electron source size in V direction in m
		self.sigmaHp_e = 2.8e-6          # Sigma electron divergence in H direction in rad
		self.sigmaVp_e = 1.5e-6          # Sigma electron divergence in V direction in rad
		
		self.sigmaH =  (self.sigmaH_e**2 +  self.wl*self.L_und/2/np.pi/np.pi)**0.5
		self.sigmaV =  (self.sigmaV_e**2 +  self.wl*self.L_und/2/np.pi/np.pi)**0.5
		self.sigmaHp = (self.sigmaHp_e**2 + self.wl/self.L_und/2)**0.5
		self.sigmaVp = (self.sigmaVp_e**2 + self.wl/self.L_und/2)**0.5

		#any update calcs to be called?

		
	def setBeamline(self):
		self.d_StoL =               # Source-to-CRL1 distance, in m
		self.d_Stof  =              # Source-to-focus distance, in m

	def setLenses(self, propertyType, propertyVal, lensNum):
		self.d_min   = 3.0e-5            # Minimum thickness at the apex in m
		self.stack_d = 50.0e-3           # Stack thickness in m
		self.L1_n    = np.array([1,      1,      1,      1,      1,      1,      2,      4,      8,      16])                # CRL1 number of lenses in each stack
		self.L1_R    = np.array([2.0e-3, 1.0e-3, 5.0e-4, 3.0e-4, 2.0e-4, 1.0e-4, 1.0e-4, 1.0e-4, 1.0e-4, 1.0e-4])            # CRL1 lens radius in each stack
		self.L1_mater= np.array(["Be",   "Be",   "Be",   "Be",   "Be",   "Be",   "Be",   "Be",   "Be",   "Be"])              # CRL1 lens material in each stack
		self.L1_loc  = np.array([4.5,    3.5,    2.5,    1.5,    0.5,    -0.5,   -1.5,   -2.5,   -3.5,   -4.5])*stack_d      # CRL1 lens stack location relative to center stack, positive means upstream
		self.L1_HE   = np.array([1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-6, 1.4e-6, 2.0e-6, 2.8e-6, 4.0e-6])            # CRL1 lens RMS thickness error

		#any update calcs to be called?
		
	def setLensCount(self, lensCount):
		pass


	def calc_lookup_table(self):
		#------------------> Needs refactoring <--------------------------------
		L1_D        = np.zeros(L1_R.size)                                   # CRL1 diameters for each stack
		for i in range(L1_R.size):
			L1_D[i] = lookup_diameter(L1_R[i])
		L1_delta    = materials_to_deltas(L1_mater, energy_keV)             # delta values for CRL1 stacks
		L1_mu       = materials_to_linear_attenuation(L1_mater, energy_keV) # mu values for CRL1 stacks
		L1_Feq      = L1_R/(2*L1_n*L1_delta) + L1_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*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 <--------------------------------

	def find_config(self):
		# match desired fsize to lookup table
		pass		
				
	def updateSlitSize(self, size, slit):
		if slit = 'hor':
			self.slit1_H = float(size)          # H slit size before CRL 1
		elif slit == 'vert':
			self.slit1_V = float(size)     # V slit size before CRL 1
		else 
			# Need error handling
			break
		self.calc_lookup_table()

	def updateE(self, energy):
		# Energy variable sent from IOC as a string
		self.energy = float(energy)
		self.calc_lookup_table()
	
	def updateFsize(self, focalSize):

		# focalPoint variable sent from IOC as a string
		self.focalSize = float(focalSize)
		self.find_config()
		
	def calc_lenses(self):
		self.lenses = (self.energy * self.focalPoint) % 4096
		pydev.iointr('new_lens_config', self.lenses)