Fixed bugs + commented out
This commit is contained in:
parent
8a96cb841d
commit
a841f10379
@ -111,6 +111,19 @@ polylog_int = CubicSpline(x_int, poly_tab)
|
|||||||
|
|
||||||
|
|
||||||
def thermal(x, x0, amp, sigma):
|
def thermal(x, x0, amp, sigma):
|
||||||
|
"""Calculating thermal density distribution in 1D (scaled such that if amp=1, return = 1)
|
||||||
|
|
||||||
|
:param x: axis
|
||||||
|
:type x: float or 1d array
|
||||||
|
:param x0: position of peak along axis
|
||||||
|
:type x0: float
|
||||||
|
:param amp: amplitude of function
|
||||||
|
:type amp: float
|
||||||
|
:param sigma: width of function
|
||||||
|
:type sigma: float
|
||||||
|
:return: calculated function value
|
||||||
|
:rtype: float or 1D array
|
||||||
|
"""
|
||||||
res = np.exp(-0.5 * (x-x0)**2 / sigma**2)
|
res = np.exp(-0.5 * (x-x0)**2 / sigma**2)
|
||||||
return amp/1.643 * polylog_int(res)
|
return amp/1.643 * polylog_int(res)
|
||||||
|
|
||||||
@ -363,6 +376,8 @@ class DensityProfileBEC2dModel(lmfit.Model):
|
|||||||
prefix='',
|
prefix='',
|
||||||
nan_policy='raise',
|
nan_policy='raise',
|
||||||
atom_n_conv=144,
|
atom_n_conv=144,
|
||||||
|
pre_check=False,
|
||||||
|
post_check=False,
|
||||||
is_debug=False,
|
is_debug=False,
|
||||||
**kwargs
|
**kwargs
|
||||||
):
|
):
|
||||||
@ -370,6 +385,8 @@ class DensityProfileBEC2dModel(lmfit.Model):
|
|||||||
'independent_vars': independent_vars})
|
'independent_vars': independent_vars})
|
||||||
|
|
||||||
self.atom_n_conv = atom_n_conv
|
self.atom_n_conv = atom_n_conv
|
||||||
|
self.pre_check = pre_check
|
||||||
|
self.post_check = post_check
|
||||||
self.is_debug=is_debug
|
self.is_debug=is_debug
|
||||||
super().__init__(density_profile_BEC_2d, **kwargs)
|
super().__init__(density_profile_BEC_2d, **kwargs)
|
||||||
self._set_paramhints_prefix()
|
self._set_paramhints_prefix()
|
||||||
@ -387,48 +404,71 @@ class DensityProfileBEC2dModel(lmfit.Model):
|
|||||||
self.set_param_hint('sigma_th', min=0)
|
self.set_param_hint('sigma_th', min=0)
|
||||||
|
|
||||||
|
|
||||||
def guess(self, data, x, y, **kwargs):
|
def guess(self, data, x, y, pre_check=False, post_check=False, **kwargs):
|
||||||
|
"""Estimate and create initial model parameters for 2d bimodal fit, by doing a 1d bimodal fit along an integrated slice of the image
|
||||||
|
|
||||||
|
:param data: Flattened 2d array, in form [a_00, a_10, a_20, ..., a_01, a_02, .. ,a_XY] with a_xy, x_dim=X, y_dim=Y
|
||||||
|
:type data: 1d numpy array
|
||||||
|
:param x: flattened X output of np.meshgrid(x_axis,y_axis) in form: [x1, x2, .., xX, x1, x2, .., xX, .. Y times ..]
|
||||||
|
:type x: 1d numpy array
|
||||||
|
:param y: flattened Y output of np.meshgrid(x_axis,y_axis) in form: [y1, y1, .., y1 (X times), y2, y2, .., y2 (X times), .. Y times ..]
|
||||||
|
:type y: 1d numpy array
|
||||||
|
:param pre_check: if True the amplitude of the 1d fit is used to guess if the image is purely BEC or thermal and
|
||||||
|
the corresponding amplitude of the 2d fit is set to zero and not varied to speed up the fitting, defaults to False
|
||||||
|
:type pre_check: bool, optional
|
||||||
|
:param post_check: if True, after doing a 2d bimodal fit the number of atoms surrounding the fitted BEC is counted and if the value is
|
||||||
|
below a certain threshhold the fit is done again with the thermal amplitude set to zero, defaults to False
|
||||||
|
:type post_check: bool, optional
|
||||||
|
:return: initial parameters for 2d fit
|
||||||
|
:rtype: params object (lmfit)
|
||||||
"""
|
"""
|
||||||
Estimate and create initial model parameters for 2d fit
|
self.pre_check = pre_check
|
||||||
:param data: Flattened 2d array, flattened from array with (x,y) -->
|
self.post_check = post_check
|
||||||
:param x:
|
|
||||||
:param y:
|
# reshaping the image to 2D in the form [[a_00, a_01, .., a_0Y], [a_10,.., a_1Y], .., [a_X0, .., a_XY]], with a_xy
|
||||||
:param kwargs:
|
|
||||||
:return:
|
|
||||||
"""
|
|
||||||
#
|
|
||||||
# global X_guess
|
|
||||||
# global bval_1d
|
|
||||||
x_width = len(np.unique(x))
|
x_width = len(np.unique(x))
|
||||||
y_width = len(np.unique(y))
|
y_width = len(np.unique(y))
|
||||||
|
|
||||||
data = np.reshape(data, (y_width, x_width))
|
data = np.reshape(data, (y_width, x_width))
|
||||||
data = data.T
|
data = data.T
|
||||||
|
|
||||||
shape = np.shape(data)
|
shape = np.shape(data)
|
||||||
cut_width = np.max(shape)
|
if self.is_debug:
|
||||||
|
print(f'shape: {shape}')
|
||||||
|
max_width = np.max(shape)
|
||||||
|
|
||||||
thresh = self.calc_thresh(data)
|
# binarizing image to guess BEC width and calculate center
|
||||||
|
thresh = self.calc_thresh(data,thresh_val=0.5)
|
||||||
|
# calculating center of cloud by statistical distribution of binarized image
|
||||||
center = self.calc_cen(thresh)
|
center = self.calc_cen(thresh)
|
||||||
|
# guessing BEC width, or better of width of center blob if no BEC is present
|
||||||
BEC_width_guess = self.guess_BEC_width(thresh, center)
|
BEC_width_guess = self.guess_BEC_width(thresh, center)
|
||||||
|
|
||||||
|
# plot binarized image and center position for debugging
|
||||||
|
if self.is_debug:
|
||||||
|
plt.pcolormesh(thresh.T, cmap='jet')
|
||||||
|
plt.plot(center[0], center[1], marker='x', markersize=25, color='green')
|
||||||
|
plt.gca().set_aspect('equal')
|
||||||
|
plt.title(f'Binarized image for guessing BEC width + center position (BEC_width: x={BEC_width_guess[0]:.0f}, y={BEC_width_guess[1]:.0f} pix)')
|
||||||
|
plt.xlabel('x_axis')
|
||||||
|
plt.ylabel('y_axis')
|
||||||
|
plt.show()
|
||||||
|
|
||||||
|
# The 1d fit is done along the short axis of the BEC (decided via the BEC_width guess)
|
||||||
if BEC_width_guess[0] < BEC_width_guess[1]:
|
if BEC_width_guess[0] < BEC_width_guess[1]:
|
||||||
if self.is_debug:
|
if self.is_debug:
|
||||||
print(f'x smaller y')
|
print(f'x smaller y, 1d fit along x')
|
||||||
s_width_ind = 0
|
s_width_ind = 0
|
||||||
|
# slice of the image along the short BEC axis with width of BEC width is taken
|
||||||
X_guess = np.sum(data[:, round(center[1] - BEC_width_guess[1]/2) : round(center[1] + BEC_width_guess[1]/2)], 1) / len(data[0,round(center[1] - BEC_width_guess[1]/2) : round(center[1] + BEC_width_guess[1]/2)])
|
X_guess = np.sum(data[:, round(center[1] - BEC_width_guess[1]/2) : round(center[1] + BEC_width_guess[1]/2)], 1) / len(data[0,round(center[1] - BEC_width_guess[1]/2) : round(center[1] + BEC_width_guess[1]/2)])
|
||||||
else:
|
else:
|
||||||
if self.is_debug:
|
if self.is_debug:
|
||||||
print(f'y smaller x')
|
print(f'y smaller x, 1d fit along y')
|
||||||
s_width_ind = 1
|
s_width_ind = 1
|
||||||
X_guess = np.sum(data[round(center[0] - BEC_width_guess[0]/2) : round(center[0] + BEC_width_guess[0]/2), :], 0) / len(data[0,round(center[0] - BEC_width_guess[0]/2) : round(center[0] + BEC_width_guess[0]/2)])
|
X_guess = np.sum(data[round(center[0] - BEC_width_guess[0]/2) : round(center[0] + BEC_width_guess[0]/2), :], 0) / len(data[0,round(center[0] - BEC_width_guess[0]/2) : round(center[0] + BEC_width_guess[0]/2)])
|
||||||
|
|
||||||
if self.is_debug:
|
|
||||||
print(f'center = {center}')
|
|
||||||
print(f'BEC widths: {BEC_width_guess}')
|
|
||||||
plt.plot(X_guess)
|
|
||||||
plt.show()
|
|
||||||
|
|
||||||
|
# Creating 1d fit init params + Performing fit
|
||||||
x = np.linspace(0, len(X_guess), len(X_guess))
|
x = np.linspace(0, len(X_guess), len(X_guess))
|
||||||
|
|
||||||
max_val = np.max(X_guess)
|
max_val = np.max(X_guess)
|
||||||
@ -436,32 +476,52 @@ class DensityProfileBEC2dModel(lmfit.Model):
|
|||||||
fitmodel_1d = lmfit.Model(density_1d, independent_vars=['x'])
|
fitmodel_1d = lmfit.Model(density_1d, independent_vars=['x'])
|
||||||
params_1d = lmfit.Parameters()
|
params_1d = lmfit.Parameters()
|
||||||
params_1d.add_many(
|
params_1d.add_many(
|
||||||
('x0_bec', center[s_width_ind], True,0, 200),
|
('x0_bec', center[s_width_ind], True, center[s_width_ind]-10, center[s_width_ind]+10),
|
||||||
('x0_th',center[s_width_ind], True,0, 200),
|
('x0_th',center[s_width_ind], True, center[s_width_ind]-10, center[s_width_ind]+10),
|
||||||
('amp_bec', 0.7 * max_val, True, 0, 1.3 * max_val),
|
('amp_bec', 0.5 * max_val, True, 0, 1.3 * max_val),
|
||||||
('amp_th', 0.3 * max_val, True, 0, 1.3 * max_val),
|
('amp_th', 0.5 * max_val, True, 0, 1.3 * max_val),
|
||||||
('deltax', 3*BEC_width_guess[s_width_ind], True, 0,cut_width),
|
('deltax', 3*BEC_width_guess[s_width_ind], True, 0, max_width),
|
||||||
# ('sigma_bec',BEC_width_guess[i,j,0]/1.22, True, 0, 50)
|
# ('sigma_bec',BEC_width_guess[i,j,0]/1.22, True, 0, 50)
|
||||||
('sigma_bec',BEC_width_guess[s_width_ind]/1.22, True, 0, 50)
|
('sigma_bec',BEC_width_guess[s_width_ind]/1.22, True, 0, BEC_width_guess[s_width_ind]*2)
|
||||||
)
|
)
|
||||||
params_1d.add('sigma_th', 3*BEC_width_guess[0], min=0, expr=f'0.632*sigma_bec + 0.518*deltax')
|
params_1d.add('sigma_th', 3*BEC_width_guess[0], min=0, expr=f'0.632*sigma_bec + 0.518*deltax')
|
||||||
|
|
||||||
|
|
||||||
res_1d = fitmodel_1d.fit(X_guess, x=x, params=params_1d)
|
res_1d = fitmodel_1d.fit(X_guess, x=x, params=params_1d)
|
||||||
|
|
||||||
if self.is_debug:
|
|
||||||
params_1d.pretty_print()
|
|
||||||
self.print_bval(res_1d)
|
|
||||||
|
|
||||||
|
|
||||||
bval_1d = res_1d.best_values
|
bval_1d = res_1d.best_values
|
||||||
|
|
||||||
|
if self.is_debug:
|
||||||
|
print('')
|
||||||
|
print('1d fit initialization')
|
||||||
|
print(f'center = {center}')
|
||||||
|
print(f'BEC widths: {BEC_width_guess}')
|
||||||
|
print('')
|
||||||
|
print('1d init fit values')
|
||||||
|
params_1d.pretty_print()
|
||||||
|
print('1d fitted values')
|
||||||
|
self.print_bval(res_1d)
|
||||||
|
|
||||||
|
x = np.linspace(0, len(X_guess), len(X_guess))
|
||||||
|
plt.plot(x, X_guess, label='1d int. data')
|
||||||
|
plt.plot(x, density_1d(x,**bval_1d), label='bimodal fit')
|
||||||
|
plt.plot(x, thermal(x,x0=bval_1d['x0_th'], amp=bval_1d['amp_th'], sigma=bval_1d['sigma_th']), label='thermal part')
|
||||||
|
plt.legend()
|
||||||
|
if s_width_ind==0:
|
||||||
|
plt.title('1d fit of data along x-axis')
|
||||||
|
plt.xlabel('x_axis (pix)')
|
||||||
|
else:
|
||||||
|
plt.title('1d fit of data along y-axis')
|
||||||
|
plt.xlabel('y_axis (pix)')
|
||||||
|
plt.show()
|
||||||
|
|
||||||
|
# scaling amplitudes of 1d fit with the maximum value of blurred 2d data
|
||||||
amp_conv_1d_2d = np.max(gaussian_filter(data, sigma=1)) / (bval_1d['amp_bec'] + bval_1d['amp_th'])
|
amp_conv_1d_2d = np.max(gaussian_filter(data, sigma=1)) / (bval_1d['amp_bec'] + bval_1d['amp_th'])
|
||||||
max_val = np.max(data)
|
max_val = np.max(data)
|
||||||
|
|
||||||
params = self.make_params()
|
params = self.make_params()
|
||||||
|
|
||||||
if bval_1d['amp_th']/bval_1d['amp_bec'] > 3:
|
# if precheck enabled and amp_th is 7x higher than amp_bec (value might be changed), amplitude of BEC in 2d fit is set to zero
|
||||||
|
if bval_1d['amp_th']/bval_1d['amp_bec'] > 7 and self.pre_check:
|
||||||
print(f'Image seems to be purely thermal (guessed from 1d fit amplitude)')
|
print(f'Image seems to be purely thermal (guessed from 1d fit amplitude)')
|
||||||
|
|
||||||
params[f'{self.prefix}amp_bec'].set(value=0, vary=False)
|
params[f'{self.prefix}amp_bec'].set(value=0, vary=False)
|
||||||
@ -472,9 +532,10 @@ class DensityProfileBEC2dModel(lmfit.Model):
|
|||||||
params[f'{self.prefix}y0_th'].set(value=center[1], min=center[1] -10, max=center[1] + 10, vary=True)
|
params[f'{self.prefix}y0_th'].set(value=center[1], min=center[1] -10, max=center[1] + 10, vary=True)
|
||||||
params[f'{self.prefix}sigmax_bec'].set(value=1, vary=False)
|
params[f'{self.prefix}sigmax_bec'].set(value=1, vary=False)
|
||||||
params[f'{self.prefix}sigmay_bec'].set(value=1, vary=False)
|
params[f'{self.prefix}sigmay_bec'].set(value=1, vary=False)
|
||||||
params[f'{self.prefix}sigma_th'].set(value=bval_1d['sigma_th'], max=cut_width, vary=True)
|
params[f'{self.prefix}sigma_th'].set(value=bval_1d['sigma_th'], max=max_width, vary=True)
|
||||||
|
|
||||||
elif bval_1d['amp_bec']/bval_1d['amp_th'] > 10:
|
# if precheck enabled and amp_bec is 10x higher than amp_th (value might be changed), amplitude of thermal part in 2d fit is set to zero
|
||||||
|
elif bval_1d['amp_bec']/bval_1d['amp_th'] > 10 and self.pre_check:
|
||||||
print('Image seems to be pure BEC (guessed from 1d fit amplitude)')
|
print('Image seems to be pure BEC (guessed from 1d fit amplitude)')
|
||||||
|
|
||||||
params[f'{self.prefix}amp_bec'].set(value=amp_conv_1d_2d * bval_1d['amp_bec'], max=1.3 * max_val, vary=True)
|
params[f'{self.prefix}amp_bec'].set(value=amp_conv_1d_2d * bval_1d['amp_bec'], max=1.3 * max_val, vary=True)
|
||||||
@ -494,6 +555,7 @@ class DensityProfileBEC2dModel(lmfit.Model):
|
|||||||
else:
|
else:
|
||||||
print('Error in small width BEC recogintion, s_width_ind should be 0 or 1')
|
print('Error in small width BEC recogintion, s_width_ind should be 0 or 1')
|
||||||
|
|
||||||
|
# params for normal 2d bimodal fit are initialized
|
||||||
else:
|
else:
|
||||||
params[f'{self.prefix}amp_bec'].set(value=amp_conv_1d_2d * bval_1d['amp_bec'], max=1.3 * max_val, vary=True)
|
params[f'{self.prefix}amp_bec'].set(value=amp_conv_1d_2d * bval_1d['amp_bec'], max=1.3 * max_val, vary=True)
|
||||||
params[f'{self.prefix}amp_th'].set(value=amp_conv_1d_2d * bval_1d['amp_th'], max=1.3 * max_val, vary=True)
|
params[f'{self.prefix}amp_th'].set(value=amp_conv_1d_2d * bval_1d['amp_th'], max=1.3 * max_val, vary=True)
|
||||||
@ -501,7 +563,7 @@ class DensityProfileBEC2dModel(lmfit.Model):
|
|||||||
params[f'{self.prefix}y0_bec'].set(value=center[1], min=center[1] -10, max=center[1] + 10, vary=True)
|
params[f'{self.prefix}y0_bec'].set(value=center[1], min=center[1] -10, max=center[1] + 10, vary=True)
|
||||||
params[f'{self.prefix}x0_th'].set(value=center[0], min=center[0] -10, max=center[0] + 10, vary=True)
|
params[f'{self.prefix}x0_th'].set(value=center[0], min=center[0] -10, max=center[0] + 10, vary=True)
|
||||||
params[f'{self.prefix}y0_th'].set(value=center[1], min=center[1] -10, max=center[1] + 10, vary=True)
|
params[f'{self.prefix}y0_th'].set(value=center[1], min=center[1] -10, max=center[1] + 10, vary=True)
|
||||||
params[f'{self.prefix}sigma_th'].set(value=bval_1d['sigma_th'], max=cut_width, vary=True)
|
params[f'{self.prefix}sigma_th'].set(value=bval_1d['sigma_th'], max=max_width, vary=True)
|
||||||
|
|
||||||
if s_width_ind == 0:
|
if s_width_ind == 0:
|
||||||
params[f'{self.prefix}sigmax_bec'].set(value=bval_1d['sigma_bec'], max= 2*BEC_width_guess[0], vary=True)
|
params[f'{self.prefix}sigmax_bec'].set(value=bval_1d['sigma_bec'], max= 2*BEC_width_guess[0], vary=True)
|
||||||
@ -516,33 +578,43 @@ class DensityProfileBEC2dModel(lmfit.Model):
|
|||||||
print('')
|
print('')
|
||||||
print('Init Params')
|
print('Init Params')
|
||||||
params.pretty_print()
|
params.pretty_print()
|
||||||
|
print('')
|
||||||
return lmfit.models.update_param_vals(params, self.prefix, **kwargs)
|
return lmfit.models.update_param_vals(params, self.prefix, **kwargs)
|
||||||
|
|
||||||
|
|
||||||
def fit(self, data, **kwargs):
|
def fit(self, data, **kwargs):
|
||||||
|
"""fitting function overwrites parent class fitting function of lmfit, in order to check (if post_check is enabled)
|
||||||
data_1d = data
|
if thermal fit completely lies in BEC fit by counting sourrounding number of atoms and comparing it to threshold value
|
||||||
|
|
||||||
res = super().fit(data_1d, **kwargs)
|
:param data: Flattened 2d array, in form [a_00, a_10, a_20, ..., a_01, a_02, .. ,a_XY] with a_xy, x_dim=X, y_dim=Y
|
||||||
|
:type data: 1d numpy array
|
||||||
|
:return: result of 2d fit
|
||||||
|
:rtype: result object (lmfit)
|
||||||
|
"""
|
||||||
|
|
||||||
|
res = super().fit(data, **kwargs)
|
||||||
|
|
||||||
if self.is_debug:
|
if self.is_debug:
|
||||||
print('bval first fit')
|
print('bval first fit')
|
||||||
self.print_bval(res)
|
self.print_bval(res)
|
||||||
|
|
||||||
|
bval = res.best_values
|
||||||
if res.params['amp_bec'].vary and res.params['amp_th'].vary:
|
# Do described post_check if enabled
|
||||||
bval = res.best_values
|
if res.params['amp_bec'].vary and res.params['amp_th'].vary and bval['amp_bec']>0.5*bval['amp_th'] and self.post_check:
|
||||||
|
|
||||||
|
# creating image by cutting out region around BEC and counting number of atoms
|
||||||
sigma_cut = max(bval['sigmay_bec'], bval['sigmax_bec'])
|
sigma_cut = max(bval['sigmay_bec'], bval['sigmax_bec'])
|
||||||
tf_fit = ThomasFermi_2d(kwargs['x'],kwargs['y'],centerx=bval['x0_bec'], centery=bval['y0_bec'], amplitude=bval['amp_bec'], sigmax=bval['sigmax_bec'], sigmay=bval['sigmay_bec'])
|
tf_fit = ThomasFermi_2d(kwargs['x'],kwargs['y'],centerx=bval['x0_bec'], centery=bval['y0_bec'], amplitude=bval['amp_bec'], sigmax=bval['sigmax_bec'], sigmay=bval['sigmay_bec'])
|
||||||
tf_fit_2 = ThomasFermi_2d(kwargs['x'],kwargs['y'],centerx=bval['x0_bec'], centery=bval['y0_bec'], amplitude=bval['amp_bec'], sigmax=1.5 * sigma_cut, sigmay=1.5*sigma_cut)
|
tf_fit_2 = ThomasFermi_2d(kwargs['x'],kwargs['y'],centerx=bval['x0_bec'], centery=bval['y0_bec'], amplitude=bval['amp_bec'], sigmax=1.5 * sigma_cut, sigmay=1.5*sigma_cut)
|
||||||
mask = np.where(tf_fit > 0, np.nan, data_1d)
|
mask = np.where(tf_fit > 0, np.nan, data)
|
||||||
#mask[i,j] = gaussian_filter(mask[i,j], sigma = 0.4)
|
|
||||||
mask = np.where(tf_fit_2 > 0, mask, np.nan)
|
mask = np.where(tf_fit_2 > 0, mask, np.nan)
|
||||||
|
|
||||||
N_c = np.nansum(mask)
|
N_c = np.nansum(mask)
|
||||||
# conversion N_count to Pixels
|
# conversion N_count to Pixels
|
||||||
N_a = self.atom_n_conv * N_c
|
N_a = self.atom_n_conv * N_c
|
||||||
|
|
||||||
|
#TODO change fixed threshhold to variable
|
||||||
|
# If number of atoms around BEC is small the image is guessed to be purely BEC and another 2d fit is performed with setting the thermal amplitude to zero
|
||||||
if N_a < 6615:
|
if N_a < 6615:
|
||||||
print('No thermal part detected, performing fit without thermal function')
|
print('No thermal part detected, performing fit without thermal function')
|
||||||
params = res.params
|
params = res.params
|
||||||
@ -551,7 +623,7 @@ class DensityProfileBEC2dModel(lmfit.Model):
|
|||||||
params[f'{self.prefix}y0_th'].set(value=1, vary=False)
|
params[f'{self.prefix}y0_th'].set(value=1, vary=False)
|
||||||
params[f'{self.prefix}sigma_th'].set(value=1, vary=False)
|
params[f'{self.prefix}sigma_th'].set(value=1, vary=False)
|
||||||
|
|
||||||
res = super().fit(data_1d, x=kwargs['x'], y=kwargs['y'], params=params)
|
res = super().fit(data, x=kwargs['x'], y=kwargs['y'], params=params)
|
||||||
|
|
||||||
return res
|
return res
|
||||||
|
|
||||||
@ -560,49 +632,38 @@ class DensityProfileBEC2dModel(lmfit.Model):
|
|||||||
def calc_thresh(self, data, thresh_val=0.3, sigma=0.4):
|
def calc_thresh(self, data, thresh_val=0.3, sigma=0.4):
|
||||||
"""Returns thresholded binary image after blurring to guess BEC size
|
"""Returns thresholded binary image after blurring to guess BEC size
|
||||||
|
|
||||||
:param data: 2d image or 1D or 2D array containing 2d images
|
:param data: 2d image
|
||||||
:type data: 2d, 3d or 4d numpy array
|
:type data: 2d numpy array
|
||||||
:param thresh_val: relative threshhold value for binarization with respect to maximum of blurred image
|
:param thresh_val: relative threshhold value for binarization with respect to maximum of blurred image
|
||||||
:param sigma: sigma of gaussian blur filter (see scipy.ndimage.gaussian_filter
|
:param sigma: sigma of gaussian blur filter (see scipy.ndimage.gaussian_filter)
|
||||||
:return: binary 2d image or 1D or 2D array containing 2d images
|
:return: binary 2d image
|
||||||
:rtype: 2d, 3d or 4d numpy array
|
:rtype: 2d numpy array
|
||||||
"""
|
"""
|
||||||
shape = np.shape(data)
|
shape = np.shape(data)
|
||||||
thresh = np.zeros(shape)
|
thresh = np.zeros(shape)
|
||||||
|
|
||||||
blurred = gaussian_filter(data, sigma=sigma)
|
blurred = gaussian_filter(data, sigma=sigma)
|
||||||
|
|
||||||
if len(shape) == 4:
|
thresh = np.where(blurred < np.max(blurred)*thresh_val, 0, 1)
|
||||||
for i in range(0,shape[0]):
|
|
||||||
for j in range(0, shape[1]):
|
|
||||||
thresh[i,j] = np.where(blurred[i,j] < np.max(blurred[i,j])*thresh_val, 0, 1)
|
|
||||||
|
|
||||||
elif len(shape) == 3:
|
|
||||||
for i in range(0,shape[0]):
|
|
||||||
thresh[i] = np.where(blurred[i] < np.max(blurred[i])*thresh_val, 0, 1)
|
|
||||||
|
|
||||||
elif len(shape) == 2:
|
|
||||||
thresh = np.where(blurred < np.max(blurred)*thresh_val, 0, 1)
|
|
||||||
|
|
||||||
|
|
||||||
else:
|
|
||||||
print("Shape of data is wrong, output is empty")
|
|
||||||
|
|
||||||
return thresh
|
return thresh
|
||||||
|
|
||||||
def calc_cen(self, thresh1):
|
def calc_cen(self, thresh1):
|
||||||
"""
|
"""Calculating the center of a blob (atom cloud) in a binarized image by first calculating the probability distribution along both axes and afterwards the expectation value
|
||||||
returns array: [X_center,Y_center]
|
|
||||||
|
:param thresh1: Binary 2D image in the form [[a_00, a_01, .., a_0Y], [a_10,.., a_1Y], .., [a_X0, .., a_XY]], with a_xy, x_dim=X, y_dim=Y
|
||||||
|
:type thresh1: 2D numpy array
|
||||||
|
:return: center coordinates of blob in form [x_center, y_center]
|
||||||
|
:rtype: 1d numpy array (shape=(1,2))
|
||||||
"""
|
"""
|
||||||
cen = np.zeros(2)
|
cen = np.zeros(2)
|
||||||
(Y,X) = np.shape(thresh1)
|
(X,Y) = np.shape(thresh1)
|
||||||
|
|
||||||
|
|
||||||
thresh1 = thresh1 /np.sum(thresh1)
|
thresh1 = thresh1 /np.sum(thresh1)
|
||||||
|
|
||||||
# marginal distributions
|
# marginal distributions
|
||||||
dx = np.sum(thresh1, 0)
|
dx = np.sum(thresh1, 1)
|
||||||
dy = np.sum(thresh1, 1)
|
dy = np.sum(thresh1, 0)
|
||||||
|
|
||||||
# expected values
|
# expected values
|
||||||
cen[0] = np.sum(dx * np.arange(X))
|
cen[0] = np.sum(dx * np.arange(X))
|
||||||
@ -610,42 +671,61 @@ class DensityProfileBEC2dModel(lmfit.Model):
|
|||||||
return cen
|
return cen
|
||||||
|
|
||||||
def guess_BEC_width(self, thresh, center):
|
def guess_BEC_width(self, thresh, center):
|
||||||
"""
|
""" returns width of blob in binarized image along both axis through the center
|
||||||
returns width of thresholded area along both axis through the center with shape of thresh and [X_width, Y_width] for each image
|
|
||||||
|
:param thresh: Binary 2D image in the form [[a_00, a_01, .., a_0Y], [a_10,.., a_1Y], .., [a_X0, .., a_XY]], with a_xy, x_dim=X, y_dim=Y
|
||||||
|
:type thresh: 2d numpy array
|
||||||
|
:param center: center of blob in image in form [x_center, y_center]
|
||||||
|
:type center: 1d numpy array (shape=(1,2))
|
||||||
|
:return: width of blob in image as [x_width, y_width]
|
||||||
|
:rtype: 1d numpy array (shape=(1,2))
|
||||||
"""
|
"""
|
||||||
shape = np.shape(thresh)
|
shape = np.shape(thresh)
|
||||||
|
|
||||||
if len(shape) == 2:
|
if len(shape) == 2:
|
||||||
BEC_width_guess = np.array([np.sum(thresh[round(center[1]), :]), np.sum(thresh[:, round(center[0])]) ])
|
BEC_width_guess = np.array([np.sum(thresh[:, round(center[1])]), np.sum(thresh[round(center[0]), :]) ])
|
||||||
|
for i in range(2):
|
||||||
|
if BEC_width_guess[i] <= 0:
|
||||||
|
BEC_width_guess[i] = 1
|
||||||
|
|
||||||
elif len(shape) == 3:
|
|
||||||
BEC_width_guess = np.zeros((shape[0], 2))
|
|
||||||
for i in range(0, shape[0]):
|
|
||||||
BEC_width_guess[i, 0] = np.sum(thresh[i, round(center[i,j,1]), :])
|
|
||||||
BEC_width_guess[i, 1] = np.sum(thresh[i, :, round(center[i,j,0])])
|
|
||||||
|
|
||||||
elif len(shape) == 4:
|
|
||||||
BEC_width_guess = np.zeros((shape[0], shape[1], 2))
|
|
||||||
for i in range(0, shape[0]):
|
|
||||||
for j in range(0, shape[1]):
|
|
||||||
BEC_width_guess[i, j, 0] = np.sum(thresh[i, j, round(center[i,j,1]), :])
|
|
||||||
BEC_width_guess[i, j, 1] = np.sum(thresh[i, j, :, round(center[i,j,0])])
|
|
||||||
else:
|
else:
|
||||||
print("Shape of data is wrong, output is empty")
|
print("Shape of data is wrong, output is empty")
|
||||||
|
|
||||||
return BEC_width_guess
|
return BEC_width_guess
|
||||||
|
|
||||||
def cond_frac(self, results):
|
def cond_frac(self, results, X, Y):
|
||||||
"""Returns condensate fraction"""
|
"""Returns condensate fraction of 2d fitting result
|
||||||
|
|
||||||
|
:param results: result of 2d bimodal fit
|
||||||
|
:type results: result object (lmfit)
|
||||||
|
:param X: X output of np.meshgrid(x_axis,y_axis) in form: [[x1, x2, .., xX], [x1, x2, .., xX] .. Y times ..]
|
||||||
|
:type X: 2d numpy array
|
||||||
|
:param Y: Y output of np.meshgrid(x_axis,y_axis) in form: [[y1, y1, .., y1 (X times)], [y2, y2, .., y2 (X times)], .. Y times ..]
|
||||||
|
:type Y: 2d numpy array
|
||||||
|
:return: condensate fraction
|
||||||
|
:rtype: float between 0 and 1
|
||||||
|
"""
|
||||||
bval = results.best_values
|
bval = results.best_values
|
||||||
tf_fit = ThomasFermi_2d(X, Y, centerx=bval['x0_bec'], centery=bval['y0_bec'], amplitude=bval['amp_bec'], sigmax=bval['sigmax_bec'], sigmay=bval['sigmay_bec'])
|
tf_fit = ThomasFermi_2d(X,Y,centerx=bval['x0_bec'], centery=bval['y0_bec'], amplitude=bval['amp_bec'], sigmax=bval['sigmax_bec'], sigmay=bval['sigmay_bec'])
|
||||||
N_bec = np.sum(tf_fit)
|
N_bec = np.sum(tf_fit)
|
||||||
fit = density_profile_BEC_2d(X, Y, **bval)
|
fit = density_profile_BEC_2d(X,Y, **bval)
|
||||||
N_ges = np.sum(fit)
|
N_ges = np.sum(fit)
|
||||||
return N_bec/N_ges
|
return N_bec/N_ges
|
||||||
|
|
||||||
def return_atom_number(self, result, X, Y, is_print=True):
|
def return_atom_number(self, result, X, Y, is_print=True):
|
||||||
"""Printing fitted atom number in bec + thermal state"""
|
"""Calculating (and printing if enabled) fitted atom number in BEC + thermal state, and condensate fraction
|
||||||
|
|
||||||
|
:param result: result of 2d bimodal fit
|
||||||
|
:type result: result object (lmfit)
|
||||||
|
:param X: X output of np.meshgrid(x_axis,y_axis) in form: [[x1, x2, .., xX], [x1, x2, .., xX] .. Y times ..]
|
||||||
|
:type X: 2d numpy array
|
||||||
|
:param Y: Y output of np.meshgrid(x_axis,y_axis) in form: [[y1, y1, .., y1 (X times)], [y2, y2, .., y2 (X times)], .. Y times ..]
|
||||||
|
:type Y: 2d numpy array
|
||||||
|
:param is_print: if true results are printed, defaults to True
|
||||||
|
:type is_print: bool, optional
|
||||||
|
:return: dictionary with total atom number N, BEC N_bec, thermal N_th and condensate fraction cond_f
|
||||||
|
:rtype: dictionary
|
||||||
|
"""
|
||||||
bval = result.best_values
|
bval = result.best_values
|
||||||
tf_fit = ThomasFermi_2d(X,Y,centerx=bval['x0_bec'], centery=bval['y0_bec'], amplitude=bval['amp_bec'], sigmax=bval['sigmax_bec'], sigmay=bval['sigmay_bec'])
|
tf_fit = ThomasFermi_2d(X,Y,centerx=bval['x0_bec'], centery=bval['y0_bec'], amplitude=bval['amp_bec'], sigmax=bval['sigmax_bec'], sigmay=bval['sigmay_bec'])
|
||||||
N_bec = self.atom_n_conv * np.sum(tf_fit)
|
N_bec = self.atom_n_conv * np.sum(tf_fit)
|
||||||
@ -653,6 +733,8 @@ class DensityProfileBEC2dModel(lmfit.Model):
|
|||||||
th_fit = polylog2_2d(X, Y, centerx=bval['x0_th'], centery=bval['y0_th'], amplitude=bval['amp_th'], sigmax=bval['sigma_th'], sigmay=bval['sigma_th'])
|
th_fit = polylog2_2d(X, Y, centerx=bval['x0_th'], centery=bval['y0_th'], amplitude=bval['amp_th'], sigmax=bval['sigma_th'], sigmay=bval['sigma_th'])
|
||||||
N_th = self.atom_n_conv * np.sum(th_fit)
|
N_th = self.atom_n_conv * np.sum(th_fit)
|
||||||
|
|
||||||
|
N = N_bec + N_th
|
||||||
|
frac = N_bec/N
|
||||||
# fit = density_profile_BEC_2d(X,Y, **bval)
|
# fit = density_profile_BEC_2d(X,Y, **bval)
|
||||||
# N_ges = self.atom_n_conv * np.sum(fit)
|
# N_ges = self.atom_n_conv * np.sum(fit)
|
||||||
|
|
||||||
@ -661,22 +743,47 @@ class DensityProfileBEC2dModel(lmfit.Model):
|
|||||||
print('Atom numbers:')
|
print('Atom numbers:')
|
||||||
print(f' N_bec: {N_bec :.0f}')
|
print(f' N_bec: {N_bec :.0f}')
|
||||||
print(f' N_th: {N_th :.0f}')
|
print(f' N_th: {N_th :.0f}')
|
||||||
print(f' N_ges: {N_bec + N_th :.0f}')
|
print(f' N_ges: {N:.0f}')
|
||||||
print(f' Cond. frac: {N_bec/(N_bec + N_th):.2f}')
|
print(f' Cond. frac: {frac *1e2:.2f} %')
|
||||||
print('')
|
print('')
|
||||||
|
|
||||||
|
atom_n = {'N' : N, 'N_bec' : N_bec, 'N_th' : N_th, 'cond_f' : frac}
|
||||||
|
return atom_n
|
||||||
|
|
||||||
def return_temperature(self, result, omg, tof, is_print=True, eff_pix=2.472e-6):
|
|
||||||
"""Returns temperature of thermal cloud"""
|
def return_temperature(self, result, tof, omg=None, is_print=True, eff_pix=2.472e-6):
|
||||||
|
"""Returns temperature of thermal cloud
|
||||||
|
|
||||||
|
:param result: result of 2d bimodal fit
|
||||||
|
:type result: result object (lmfit)
|
||||||
|
:param tof: time of flight
|
||||||
|
:type tof: float
|
||||||
|
:param omg: geometric average of trapping frequencies, defaults to None
|
||||||
|
:type omg: float, if NONE initial cloud size is neglected optional
|
||||||
|
:param is_print: if True temperature is printed, defaults to True
|
||||||
|
:type is_print: bool, optional
|
||||||
|
:param eff_pix: effective pixel size of imaging setup, defaults to 2.472e-6
|
||||||
|
:type eff_pix: float, optional
|
||||||
|
:return: temperature of atom cloud
|
||||||
|
:rtype: float
|
||||||
|
"""
|
||||||
R_th = result.best_values['sigma_th'] * eff_pix * np.sqrt(2)
|
R_th = result.best_values['sigma_th'] * eff_pix * np.sqrt(2)
|
||||||
print(R_th)
|
# print(R_th)
|
||||||
T = R_th**2 * 164*const.u/const.k * (1/omg**2 + tof**2)**(-1)
|
if omg is None:
|
||||||
|
T = R_th**2 * 164*const.u/const.k * (tof**2)**(-1)
|
||||||
|
else:
|
||||||
|
T = R_th**2 * 164*const.u/const.k * (1/omg**2 + tof**2)**(-1)
|
||||||
|
|
||||||
if is_print:
|
if is_print:
|
||||||
print(f'Temperature: {T*1e9:.2f} nK')
|
print(f'Temperature: {T*1e9:.2f} nK')
|
||||||
return T
|
return T
|
||||||
|
|
||||||
def print_bval(self, res_s):
|
def print_bval(self, res_s):
|
||||||
"""nicely print best fitted values + init values + bounds """
|
"""nicely print best fitted values + init values + bounds
|
||||||
|
|
||||||
|
:param res_s: result of 2d bimodal fit
|
||||||
|
:type res_s: result object (lmfit)
|
||||||
|
"""
|
||||||
keys = res_s.best_values.keys()
|
keys = res_s.best_values.keys()
|
||||||
bval = res_s.best_values
|
bval = res_s.best_values
|
||||||
init = res_s.init_params
|
init = res_s.init_params
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user