import numpy as np
from scipy import special
from scipy.optimize import curve_fit
__all__ = ['knife_edge', 'knife_edge_base']
[docs]def knife_edge(positions, intensities, axisRange=None, p0=None):
"""
Calculates the beam radius at 1/e^2 from a knife-edge scan by
fitting with erfc function: f(a,b,u) = a*erfc(u) + b or
where u = sqrt(2)*(x-x0)/w0 with w0 the beam radius at 1/e^2
and x0 the beam center.
Parameters
----------
positions : np.ndarray
Motor position values, typically 1D
intensities : np.ndarray
Intensity values, could be either 1D or 2D, with the number or rows
equivalent with the position size
axisRange : sequence of two floats or None
Edges of the scanning axis between which to apply the fit.
p0 : list of floats, numpy 1D array
Initial parameters used for the fit: x0, w0, a, b. If None, a beam
radius of 100 um is assumed.
Returns
-------
width : float
The beam radius at 1/e^2
std : float
The standard deviation of the width
"""
popt, pcov = knife_edge_base(positions, intensities,
axisRange=axisRange, p0=p0)
width, std = 0, 0
if popt is not None and pcov is not None:
width, std = np.abs(popt[1]), pcov[1, 1]**0.5
return width, std
[docs]def knife_edge_base(positions, intensities, axisRange=None, p0=None):
"""
The base implementation of the knife-edge scan analysis.
Calculates the beam radius at 1/e^2 from a knife-edge scan by
fitting with erfc function: f(a,b,u) = a*erfc(u) + b or
where u = sqrt(2)*(x-x0)/w0 with w0 the beam radius at 1/e^2
and x0 the beam center.
Parameters
----------
positions : np.ndarray
Motor position values, typically 1D
intensities : np.ndarray
Intensity values, could be either 1D or 2D, with the number or rows
equivalent with the position size
axisRange : sequence of two floats or None
Edges of the scanning axis between which to apply the fit.
p0 : list of floats, numpy 1D array
Initial parameters used for the fit: x0, w0, a, b. If None, a beam
radius of 100 um is assumed.
Returns
-------
popt : sequence of floats or None
The parameters of the resulting fit.
pcov : sequence of floats
The covariance matrix of the resulting fit.
"""
# Prepare arrays
positions, intensities = prepare_arrays(positions, intensities,
xRange=axisRange)
# Estimate initial fitting params
if p0 is None:
p0 = [np.mean(positions), 0.1, np.max(intensities) / 2, 0]
# Fit
popt, pcov = function_fit(erfc, positions, intensities, p0=p0)
return popt, pcov
def function_fit(func, x, y, **kwargs):
"""A wrapper for scipy.optimize curve_fit()
"""
# Fit
try:
popt, pcov = curve_fit(func, x, y, **kwargs)
except (TypeError, RuntimeError) as err:
print("Fit did not converge:", err)
popt, pcov = (None, None)
return popt, pcov
def prepare_arrays(arrX: np.ndarray, arrY: np.ndarray,
xRange=None, yRange=None):
"""
Preprocessing of the input x and y arrays.
This involves the following steps.
1. Converting the arrays to 1D of the same size
2. Select the ranges from the input x- and y-ranges
3. Retrieve finite values.
"""
# Convert both arrays to 1D of the same size
arrX, arrY = arrays_to1d(arrX, arrY)
# Select ranges
if xRange is not None:
low, high = xRange
if low == high:
raise ValueError('The supplied xRange is not a valid range.')
mask_ = range_mask(arrX, low, high)
arrX = arrX[mask_]
arrY = arrY[mask_]
if yRange is not None:
low, high = yRange
if low == high:
raise ValueError('The supplied xRange is not a valid range.')
mask_ = range_mask(arrY, low, high)
arrX = arrX[mask_]
arrY = arrY[mask_]
# Clean both arrays by only getting finite values
finite_idx = np.isfinite(arrX) & np.isfinite(arrY)
arrX = arrX[finite_idx]
arrY = arrY[finite_idx]
return arrX, arrY
def arrays_to1d(arrX: np.ndarray, arrY: np.ndarray):
"""Flatten two arrays and matches their sizes
"""
assert arrX.shape[0] == arrY.shape[0]
arrX, arrY = arrX.flatten(), arrY.flatten()
if len(arrX) > len(arrY):
arrY = np.repeat(arrY, len(arrX) // len(arrY))
else:
arrX = np.repeat(arrX, len(arrY) // len(arrX))
return arrX, arrY
def range_mask(array, minimum=None, maximum=None):
"""Retrieve the resulting array from the given minimum and maximum
"""
default = np.ones(array.shape, dtype=bool)
min_slice, max_slice = default, default
if minimum is not None:
if minimum > np.nanmax(array):
raise ValueError('The range minimum is too large.')
min_slice = array >= minimum
if maximum is not None:
if maximum < np.nanmin(array):
raise ValueError('The range maximum is too small.')
max_slice = array <= maximum
return min_slice & max_slice
def erfc(x, x0, w0, a, b):
return a * special.erfc(np.sqrt(2) * (x - x0) / w0) + b