delete personal scripts

Analyser/FitAnalyser.py

@@ -25,11 +25,16 @@ from scipy.special import erf, erfc
 from scipy.special import gamma as gamfcn
 from scipy.special import wofz
 from scipy.optimize import curve_fit
+from scipy.interpolate import CubicSpline
+from scipy.ndimage import gaussian_filter, rotate
+import scipy.constants as const
 
 import xarray as xr
 
 import copy
 
+import matplotlib.pyplot as plt
+
 
 log2 = log(2)
 s2pi = sqrt(2*pi)
@@ -82,40 +87,100 @@ def two_gaussian2d(x, y=0.0, A_amplitude=1.0, A_centerx=0.0, A_centery=0.0, A_si
     return z
 
 
-def ThomasFermi_2d(x, y=0.0, centerx=0.0, centery=0.0, amplitude=1.0, sigmax=1.0, sigmay=1.0): 
+def polylog_tab(pow, x, order=100):
-    res = (1- ((x-centerx)/(sigmax))**2 - ((y-centery)/(sigmay))**2)**(3 / 2)
+    """Calculationg the polylog sum_i(x^i/i^pow), up to the order-th element of the sum
-    return amplitude * 5 / 2 / np.pi / max(tiny, sigmax * sigmay) * np.where(res > 0, res, 0)
+
+    :param pow: power in denominator of sum
+    :type pow: int (can be float)
+    :param x: argument of Polylogarithm
+    :type x: float or 1D numpy array
+    :return: value of polylog(x)
+    :rtype: same as x, float or 1D numpy array
+    """
+
+    sum = 0
+    for k in range(1,order):
+        sum  += x ** k /k **pow
+    return sum
+
+# Creating array for interpolation
+x_int = np.linspace(0, 1.00001, 100000)
+poly_tab = polylog_tab(2,x_int)
+
+# Creating function interpolating between the Polylog values calculated before
+polylog_int = CubicSpline(x_int, poly_tab)
 
 
-def polylog(power, numerator):
+def thermal(x, x0, amp, sigma):
-    
+    """Calculating thermal density distribution in 1D (scaled such that if amp=1, return = 1)
-    order = 20
-    
-    dataShape = numerator.shape
-    numerator = np.tile(numerator, (order, 1))
-    numerator = np.power(numerator.T, np.arange(1, order+1)).T
 
-    denominator = np.arange(1, order+1)
+    :param x: axis
-    denominator = np.tile(denominator, (dataShape[0], 1))
+    :type x: float or 1d array
-    denominator = denominator.T
+    :param x0: position of peak along axis
-
+    :type x0: float
-    data = numerator/ np.power(denominator, power)
+    :param amp: amplitude of function
-
+    :type amp: float
-    return np.sum(data, axis=0)
+    :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)
+    return amp/1.643 * polylog_int(res)
 
 
-def polylog2_2d(x, y=0.0, centerx=0.0, centery=0.0, amplitude=1.0, sigmax=1.0, sigmay=1.0): 
+def Thomas_Fermi_1d(x, x0, amp, sigma):
+    res = (1- ((x-x0)/sigma)**2)
+    res = np.where(res > 0, res, 0)
+    res = res**(3/2)
+    return amp * res
+
+
+def density_1d(x, x0_bec, x0_th, amp_bec, amp_th, sigma_bec, sigma_th):
+    return thermal(x, x0_th, amp_th, sigma_th) + Thomas_Fermi_1d(x, x0_bec, amp_bec, sigma_bec)
+
+
+def polylog(pow, x):
+    order = 15
+    sum = 0
+    for k in range(1,order):
+        sum  += x ** k /k **pow
+    return sum
+
+
+def ThomasFermi_2d(x, y=0.0, centerx=0.0, centery=0.0, amplitude=1.0, sigmax=1.0, sigmay=1.0):
+
+    res = (1- ((x-centerx)/(sigmax))**2 - ((y-centery)/(sigmay))**2)
+    res = np.where(res > 0, res, 0)
+    res = res**(3/2)
+    return amplitude * res
+
+
+def polylog2_2d(x, y=0.0, centerx=0.0, centery=0.0, amplitude=1.0, sigmax=1.0, sigmay=1.0):
     ## Approximation of the polylog function with 2D gaussian as argument. -> discribes the thermal part of the cloud
-    return amplitude / 2 / np.pi / 1.20206 / max(tiny, sigmax * sigmay) * polylog(2, np.exp( -((x-centerx)**2/(2 * (sigmax)**2))-((y-centery)**2/( 2 * (sigmay)**2)) ))
+    return amplitude/1.643  * polylog_int(np.exp( -((x-centerx)**2/(2 * sigmax**2))-((y-centery)**2/( 2 * sigmay**2)) ))
 
-
+def rotate_coord(x,y, rot_angle):
-def density_profile_BEC_2d(x, y=0.0, BEC_amplitude=1.0, thermal_amplitude=1.0, BEC_centerx=0.0, BEC_centery=0.0, thermal_centerx=0.0, thermal_centery=0.0, 
-                        BEC_sigmax=1.0, BEC_sigmay=1.0, thermal_sigmax=1.0, thermal_sigmay=1.0):
     
-    return ThomasFermi_2d(x=x, y=y, centerx=BEC_centerx, centery=BEC_centery, 
+    rot_angle *= 2*np.pi/360
-                          amplitude=BEC_amplitude, sigmax=BEC_sigmax, sigmay=BEC_sigmay
+    x_og = x
-            ) + polylog2_2d(x=x, y=y, centerx=thermal_centerx, centery=thermal_centery, 
+    x = x_og * np.cos(rot_angle) + y * np.sin(rot_angle)
-                            amplitude=thermal_amplitude, sigmax=thermal_sigmax, sigmay=thermal_sigmay)
+    y = -x_og * np.sin(rot_angle) + y * np.cos(rot_angle)
+    return x, y
+    
+def density_profile_BEC_2d(x, y=0.0, amp_bec=1.0, amp_th=1.0, x0_bec=0.0, y0_bec=0.0, x0_th=0.0, y0_th=0.0,
+                           sigmax_bec=1.0, sigmay_bec=1.0, sigma_th=1.0, rot_angle=None):
+    rot_angle *= -1
+    if rot_angle is not None:
+        x, y = rotate_coord(x,y, rot_angle)
+        x0_bec, y0_bec = rotate_coord(x0_bec, y0_bec, rot_angle)
+        x0_th, y0_th = rotate_coord(x0_th, y0_th, rot_angle)
+
+    
+    return ThomasFermi_2d(x=x, y=y, centerx=x0_bec, centery=y0_bec,
+                          amplitude=amp_bec, sigmax=sigmax_bec, sigmay=sigmay_bec
+                          ) + polylog2_2d(x=x, y=y, centerx=x0_th, centery=y0_th,
+                                          amplitude=amp_th, sigmax=sigma_th,sigmay=sigma_th)
 
 
 class GaussianWithOffsetModel(Model):
@@ -318,115 +383,479 @@ class ThomasFermi2dModel(Model):
         return update_param_vals(pars, self.prefix, **kwargs)
     
 
-class DensityProfileBEC2dModel(Model):
+class DensityProfileBEC2dModel(lmfit.Model):
-
+    """ Fitting class to do 2D bimodal fit on OD of absorption images
-    fwhm_factor = 2*np.sqrt(2*np.log(2))
+    """
-    height_factor = 1./2*np.pi
+    def __init__(self,
-
+                independent_vars=['x', 'y'],
-    def __init__(self, independent_vars=['x', 'y'], prefix='', nan_policy='raise',
+                prefix='',
-                 **kwargs):
+                nan_policy='raise',
+                atom_n_conv=144,
+                pre_check=False,
+                post_check=False,
+                is_debug=False,
+                **kwargs
+        ):
         kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,
                        'independent_vars': independent_vars})
+
+        self.atom_n_conv = atom_n_conv
+        self.pre_check = pre_check
+        self.post_check = post_check
+        self.is_debug=is_debug
         super().__init__(density_profile_BEC_2d, **kwargs)
         self._set_paramhints_prefix()
 
     def _set_paramhints_prefix(self):
         # self.set_param_hint('BEC_sigmax', min=0)
-        self.set_param_hint('deltax', min=0)
+        self.set_param_hint('amp_bec', min=0)
-        self.set_param_hint('BEC_sigmax', expr=f'3 * {self.prefix}thermal_sigmax - {self.prefix}deltax')
+        self.set_param_hint('amp_th', min=0)
-        
+        self.set_param_hint('x0_bec', min=0)
-        self.set_param_hint('BEC_sigmay', min=0)
+        self.set_param_hint('y0_bec', min=0)
-        self.set_param_hint('thermal_sigmax', min=0)
+        self.set_param_hint('x0_th', min=0)
-        # self.set_param_hint('thermal_sigmay', min=0)
+        self.set_param_hint('y0_th', min=0)
-        self.set_param_hint('BEC_amplitude', min=0)
+        self.set_param_hint('sigmax_bec', min=0)
-        self.set_param_hint('thermal_amplitude', min=0)
+        self.set_param_hint('sigmay_bec', min=0)
-        
+        self.set_param_hint('sigma_th', min=0)
-        self.set_param_hint('thermalAspectRatio', min=0.8, max=1.2)
-        self.set_param_hint('thermal_sigmay', expr=f'{self.prefix}thermalAspectRatio * {self.prefix}thermal_sigmax')
-        
-        # self.set_param_hint('betax', value=0)
-        # self.set_param_hint('BEC_centerx', expr=f'{self.prefix}thermal_sigmax - {self.prefix}betax')
-        
-        self.set_param_hint('condensate_fraction', expr=f'{self.prefix}BEC_amplitude / ({self.prefix}BEC_amplitude + {self.prefix}thermal_amplitude)')
 
-    def guess(self, data, x, y, negative=False, pureBECThreshold=0.5, noBECThThreshold=0.0, **kwargs):
+        self.set_param_hint('rot_angle', min=-90, max=90)
-        """Estimate initial model parameter values from data."""
-        fitModel = TwoGaussian2dModel()
-        pars = fitModel.guess(data, x=x, y=y, negative=negative)
-        pars['A_amplitude'].set(min=0)
-        pars['B_amplitude'].set(min=0)
-        pars['A_centerx'].set(min=pars['A_centerx'].value - 3 * pars['A_sigmax'], 
-                              max=pars['A_centerx'].value + 3 * pars['A_sigmax'],)
-        pars['A_centery'].set(min=pars['A_centery'].value - 3 * pars['A_sigmay'], 
-                              max=pars['A_centery'].value + 3 * pars['A_sigmay'],)
-        pars['B_centerx'].set(min=pars['B_centerx'].value - 3 * pars['B_sigmax'], 
-                              max=pars['B_centerx'].value + 3 * pars['B_sigmax'],)
-        pars['B_centery'].set(min=pars['B_centery'].value - 3 * pars['B_sigmay'], 
-                              max=pars['B_centery'].value + 3 * pars['B_sigmay'],)
         
-        fitResult = fitModel.fit(data, x=x, y=y, params=pars, **kwargs)
+        self.set_param_hint('atom_number_bec', expr=f'{self.prefix}amp_bec / 5 * 2 * 3.14159265359 * {self.prefix}sigmax_bec * {self.prefix}sigmay_bec')
-        pars_guess = fitResult.params
+        self.set_param_hint('atom_number_th', expr=f'{self.prefix}amp_th * 2 * 3.14159265359 * 1.20206 / 1.643 * {self.prefix}sigma_th * {self.prefix}sigma_th')
+        self.set_param_hint('condensate_fraction', expr=f'{self.prefix}atom_number_bec / ({self.prefix}atom_number_bec + {self.prefix}atom_number_th)')
+
+    def guess(self, data, x, y, rot_angle=0, vary_rot=False, is_debug=False, 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 rot_angle: angle in degrees, The image is rotated counterclockwise by this angle to match the two axes of the cloud with x and y 
+                            for the guessing procedure. The 2d-fit is done by rotating the fitting function clockwise with this angle
+        :type rot_angle: float
+        :param vary_rot: if True the angle is varied in the 2d-fit
+        :type vary_rot: bool, optional
+        :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)
+        """
+        self.pre_check = pre_check 
+        self.post_check = post_check
+        self.is_debug = is_debug
         
-        BEC_amplitude = pars_guess['A_amplitude'].value
+        # 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
-        thermal_amplitude = pars_guess['B_amplitude'].value
+        x_width = len(np.unique(x))
-                
+        y_width = len(np.unique(y))
-        pars = self.make_params(BEC_amplitude=BEC_amplitude,
+        x_1d = np.linspace(x[0], x[-1], x_width)
-                                thermal_amplitude=thermal_amplitude, 
+        y_1d = np.linspace(y[0], y[-1], y_width)
-                                BEC_centerx=pars_guess['A_centerx'].value, BEC_centery=pars_guess['A_centery'].value,
+
-                                # BEC_sigmax=(pars_guess['A_sigmax'].value / 2.355), 
+        data = np.reshape(data, (y_width, x_width))
-                                deltax = 3 * (pars_guess['B_sigmax'].value * s2) - (pars_guess['A_sigmax'].value / 2.355),
+        data = data.T
-                                BEC_sigmay=(pars_guess['A_sigmay'].value / 2.355), 
+
-                                thermal_centerx=pars_guess['B_centerx'].value, thermal_centery=pars_guess['B_centery'].value,
+        if is_debug:
-                                thermal_sigmax=(pars_guess['B_sigmax'].value * s2), 
+            X, Y = np.meshgrid(x_1d,y_1d)
-                                thermalAspectRatio=(pars_guess['B_sigmax'].value * s2) / (pars_guess['B_sigmay'].value * s2)
+            plt.pcolormesh(X,Y, data.T, cmap='jet')
-                                # thermal_sigmay=(pars_guess['B_sigmay'].value * s2)
+            plt.gca().set_aspect('equal')
-                                )
+            plt.title(f'Input data')
+            plt.xlabel('x_axis')
+            plt.ylabel('y_axis')
+            plt.show()
+
+        # the image is rotated counterclockwise by rot_angle, CAREFUL: The image has the form a_xy (last coordinate y) and therefore the rotation is done counter-clockwise. 
+        # Doing the same with a standard image with a_yx rotates it clockwise!
+        if rot_angle != 0:
+            data = rotate(data, rot_angle, reshape=False)
+
+        shape = np.shape(data)
+        if self.is_debug:
+            print(f'shape: {shape}')
+        max_width = np.max(shape)
+
+        # 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_pix = self.calc_cen_pix(thresh)
+        center = self.center_pix_conv(center_pix, x_1d, y_1d)
+        # guessing BEC width, or better of width of center blob if no BEC is present
+        BEC_width_guess = self.guess_BEC_width(thresh, center_pix)
+
+        # plot binarized image and center position for debugging
+        if self.is_debug:
+            X, Y = np.meshgrid(x_1d,y_1d)
+            plt.pcolormesh(X,Y, 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 self.is_debug:
+                print(f'x smaller y, 1d fit along x')
+            s_width_ind = 0
+            x_fit = x_1d
+            # slice of the image along the short BEC axis with width of BEC width is taken
+            X_guess = np.sum(data[:, round(center_pix[1] - BEC_width_guess[1]/2) : round(center_pix[1] + BEC_width_guess[1]/2)], 1) / len(data[0,round(center_pix[1] - BEC_width_guess[1]/2) : round(center_pix[1] + BEC_width_guess[1]/2)])
+        else:
+            if self.is_debug:
+                print(f'y smaller x, 1d fit along y')
+            s_width_ind = 1
+            x_fit = y_1d
+            X_guess = np.sum(data[round(center_pix[0] - BEC_width_guess[0]/2) : round(center_pix[0] + BEC_width_guess[0]/2), :], 0) / len(data[0,round(center_pix[0] - BEC_width_guess[0]/2) : round(center_pix[0] + BEC_width_guess[0]/2)])
+
+
+        # Creating 1d fit init params + Performing fit 
         
-        nBEC = pars[f'{self.prefix}BEC_amplitude'] / 2 / np.pi / 5.546 / pars[f'{self.prefix}BEC_sigmay'] / pars[f'{self.prefix}BEC_sigmax']
+
-        if (pars[f'{self.prefix}condensate_fraction']>0.95) and (np.max(data) > 1.05 * nBEC):
+        max_val = np.max(X_guess)
-            temp = ((np.max(data) - nBEC) * s2pi * pars[f'{self.prefix}thermal_sigmay'] / pars[f'{self.prefix}thermal_sigmax'])
+
-            if temp > pars[f'{self.prefix}BEC_amplitude']:
+        fitmodel_1d = lmfit.Model(density_1d, independent_vars=['x'])
-                pars[f'{self.prefix}thermal_amplitude'].set(value=pars[f'{self.prefix}BEC_amplitude'] / 2)
+        params_1d = lmfit.Parameters()
+        params_1d.add_many(
+            ('x0_bec', center[s_width_ind], True, center[s_width_ind]-10, center[s_width_ind]+10),
+            ('x0_th',center[s_width_ind], True, center[s_width_ind]-10, center[s_width_ind]+10),
+            ('amp_bec', 0.5 * 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, max_width),
+            # ('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, 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')
+
+
+        res_1d = fitmodel_1d.fit(X_guess, x=x_fit, params=params_1d)
+        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)
+
+            plt.plot(x_fit, X_guess, label='1d int. data')
+            plt.plot(x_fit, density_1d(x_fit,**bval_1d), label='bimodal fit')
+            plt.plot(x_fit, thermal(x_fit,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_')
             else:
-                pars[f'{self.prefix}thermal_amplitude'].set(value=temp * 10)
+                plt.title('1d fit of data along y-axis')
-        
+                plt.xlabel('y_axis')
-        if BEC_amplitude / (thermal_amplitude + BEC_amplitude) > pureBECThreshold:
+            plt.show()
-            pars[f'{self.prefix}thermal_amplitude'].set(value=0)
+
-            pars[f'{self.prefix}BEC_amplitude'].set(value=(thermal_amplitude + BEC_amplitude))
+        # 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'])
-        if BEC_amplitude / (thermal_amplitude + BEC_amplitude) < noBECThThreshold:
+        max_val = np.max(data)
-            pars[f'{self.prefix}BEC_amplitude'].set(value=0)
+
-            pars[f'{self.prefix}thermal_amplitude'].set(value=(thermal_amplitude + BEC_amplitude))
+        params = self.make_params()
-        
+
-        pars[f'{self.prefix}BEC_centerx'].set(
+        # 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
-            min=pars[f'{self.prefix}BEC_centerx'].value - 10 * pars[f'{self.prefix}BEC_sigmax'].value,
+        if bval_1d['amp_th']/bval_1d['amp_bec'] > 7 and self.pre_check:
-            max=pars[f'{self.prefix}BEC_centerx'].value + 10 * pars[f'{self.prefix}BEC_sigmax'].value,  
+            print(f'Image seems to be purely thermal (guessed from 1d fit amplitude)')
-        )
+
-        
+            params[f'{self.prefix}amp_bec'].set(value=0, vary=False)
-        pars[f'{self.prefix}thermal_centerx'].set(
+            params[f'{self.prefix}amp_th'].set(value=amp_conv_1d_2d * bval_1d['amp_th'], max=1.3 * max_val, vary=True)
-            min=pars[f'{self.prefix}thermal_centerx'].value - 3 * pars[f'{self.prefix}thermal_sigmax'].value,
+            params[f'{self.prefix}x0_bec'].set(value=1, vary=False)
-            max=pars[f'{self.prefix}thermal_centerx'].value + 3 * pars[f'{self.prefix}thermal_sigmax'].value,  
+            params[f'{self.prefix}y0_bec'].set(value=1, vary=False)
-        )
+            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)
-        pars[f'{self.prefix}BEC_centery'].set(
+            params[f'{self.prefix}sigmax_bec'].set(value=1, vary=False)
-            min=pars[f'{self.prefix}BEC_centery'].value - 10 * pars[f'{self.prefix}BEC_sigmay'].value,
+            params[f'{self.prefix}sigmay_bec'].set(value=1, vary=False)
-            max=pars[f'{self.prefix}BEC_centery'].value + 10 * pars[f'{self.prefix}BEC_sigmay'].value,  
+            params[f'{self.prefix}sigma_th'].set(value=bval_1d['sigma_th'], max=max_width, vary=True)
-        )
+
-        
+        # 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
-        pars[f'{self.prefix}thermal_centery'].set(
+        elif bval_1d['amp_bec']/bval_1d['amp_th'] > 10 and self.pre_check:
-            min=pars[f'{self.prefix}thermal_centery'].value - 3 * pars[f'{self.prefix}thermal_sigmay'].value,
+            print('Image seems to be pure BEC (guessed from 1d fit amplitude)')
-            max=pars[f'{self.prefix}thermal_centery'].value + 3 * pars[f'{self.prefix}thermal_sigmay'].value,  
+
-        )
+            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=0, vary=False)
-        pars[f'{self.prefix}BEC_sigmay'].set(
+            params[f'{self.prefix}x0_bec'].set(value=center[0], min=center[0] -10, max=center[0] + 10, vary=True)
-            max=5 * pars[f'{self.prefix}BEC_sigmay'].value,  
+            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=1, vary=False)
-        
+            params[f'{self.prefix}y0_th'].set(value=1, vary=False)
-        pars[f'{self.prefix}thermal_sigmax'].set(
+            params[f'{self.prefix}sigma_th'].set(value=1, vary=False)
-            max=5 * pars[f'{self.prefix}thermal_sigmax'].value,   
+
-        )
+            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)
-        return update_param_vals(pars, self.prefix, **kwargs)
+                params[f'{self.prefix}sigmay_bec'].set(value=BEC_width_guess[1]/1.22, max= 2*BEC_width_guess[1], vary=True)
+            elif s_width_ind == 1:
+                params[f'{self.prefix}sigmax_bec'].set(value=BEC_width_guess[0]/1.22, max= 2*BEC_width_guess[0], vary=True)
+                params[f'{self.prefix}sigmay_bec'].set(value=bval_1d['sigma_bec'], max= 2*BEC_width_guess[1], vary=True)
+            else:
+                print('Error in small width BEC recogintion, s_width_ind should be 0 or 1')
+
+        # params for normal 2d bimodal fit are initialized
+        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_th'].set(value=amp_conv_1d_2d * bval_1d['amp_th'], max=1.3 * max_val, vary=True)
+            params[f'{self.prefix}x0_bec'].set(value=center[0], min=center[0] -10, max=center[0] + 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}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=max_width, vary=True)
+
+            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}sigmay_bec'].set(value=BEC_width_guess[1]/1.22, max= 2*BEC_width_guess[1], vary=True)
+            elif s_width_ind == 1:
+                params[f'{self.prefix}sigmax_bec'].set(value=BEC_width_guess[0]/1.22, max= 2*BEC_width_guess[0], vary=True)
+                params[f'{self.prefix}sigmay_bec'].set(value=bval_1d['sigma_bec'], max= 2*BEC_width_guess[1], vary=True)
+            else:
+                print('Error in small width BEC recogintion, s_width_ind should be 0 or 1')
+
+        params[f'{self.prefix}rot_angle'].set(value=rot_angle, min=rot_angle-30, max=rot_angle+30, vary=vary_rot)
+
+        if self.is_debug:
+            print('')
+            print('Init Params')
+            params.pretty_print()
+            print('')
+        return lmfit.models.update_param_vals(params, self.prefix, **kwargs)
+
+
+    def fit(self, data, **kwargs):
+        """fitting function overwrites parent class fitting function of lmfit, in order to check (if post_check is enabled) 
+            if thermal fit completely lies in BEC fit by counting sourrounding number of atoms and comparing it to threshold value
+
+        :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:
+            print('bval first fit')
+            self.print_bval(res)
+
+        bval = res.best_values
+        # Do described post_check if enabled
+
+        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'])
+            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)
+            mask = np.where(tf_fit > 0, np.nan, data)
+            mask = np.where(tf_fit_2 > 0, mask, np.nan)
+
+            N_c = np.nansum(mask)
+            # conversion N_count to Pixels
+            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:
+                print('No thermal part detected, performing fit without thermal function')
+                params = res.params
+                params[f'{self.prefix}amp_th'].set(value=0, vary=False)
+                params[f'{self.prefix}x0_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)
+
+                res = super().fit(data, x=kwargs['x'], y=kwargs['y'], params=params)
+
+                return res
+
+        return res
+
+    def calc_thresh(self, data, thresh_val=0.3, sigma=0.4):
+        """Returns thresholded binary image after blurring to guess BEC size
+
+        :param data: 2d image 
+        :type data: 2d numpy array
+        :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)
+        :return: binary 2d image 
+        :rtype: 2d numpy array
+        """
+        shape = np.shape(data)
+        thresh = np.zeros(shape)
+
+        blurred = gaussian_filter(data, sigma=sigma)
+
+        thresh = np.where(blurred < np.max(blurred)*thresh_val, 0, 1)
+
+        return thresh
+
+    def calc_cen_pix(self, thresh1):
+        """Calculating the center (in pixel) of a blob (atom cloud) in a binarized image by first calculating the probability distribution along both axes and afterwards the expectation value
+
+        :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)
+        (X,Y) = np.shape(thresh1)
+
+        thresh1 = thresh1 /np.sum(thresh1)
+
+        # marginal distributions
+        dx = np.sum(thresh1, 1)
+        dy = np.sum(thresh1, 0)
+
+        # expected values
+        cen[0] = np.sum(dx * np.arange(X))
+        cen[1] = np.sum(dy * np.arange(Y))
+        return cen
+    
+    def center_pix_conv(self, center_pix, x, y):
+        """Converts center in pixel to center in values of x and y 
+
+        :param center_pix: pixel values of center
+        :type center_pix: numpy array (shape=(1,2))
+        :param x: x-axis
+        :type x: 1d numpy array
+        :param y: y-axis
+        :type y: 1d numpy array
+        :return: center coordinates in form [x_center, y_center] with respect to the axes
+        :rtype: numpy array (shap=(1,2))
+        """
+        center = np.empty(2)
+        center[0] = x[round(center_pix[0])]
+        center[1] = y[round(center_pix[1])]
+        return center
+
+    def guess_BEC_width(self, thresh, center):
+        """ returns width of blob in binarized image along both axis through the center 
+
+        :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] in pixel
+        :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)
+
+        if len(shape) == 2:
+            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
+
+        else:
+            print("Shape of data is wrong, output is empty")
+
+        return BEC_width_guess
+
+    def cond_frac(self, results, X, Y):
+        """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
+        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)
+        fit = density_profile_BEC_2d(X,Y, **bval)
+        N_ges = np.sum(fit)
+        return N_bec/N_ges
+
+    def return_atom_number(self, result, X, Y, is_print=True):
+        """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
+        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)
+
+        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 = N_bec + N_th
+        frac = N_bec/N
+        # fit = density_profile_BEC_2d(X,Y, **bval)
+        # N_ges = self.atom_n_conv * np.sum(fit)
+
+        if is_print:
+            print()
+            print('Atom numbers:')
+            print(f'   N_bec: {N_bec :.0f}')
+            print(f'   N_th: {N_th :.0f}')
+            print(f'   N_ges: {N:.0f}')
+            print(f'   Cond. frac: {frac *1e2:.2f} %')
+            print('')
+
+        atom_n = {'N' : N, 'N_bec' : N_bec, 'N_th' : N_th, 'cond_f' : frac}
+        return atom_n
+
+
+    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)
+        # print(R_th)
+        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:
+            print(f'Temperature: {T*1e9:.2f} nK')
+        return T
+
+    def print_bval(self, res_s):
+        """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()
+        bval = res_s.best_values
+        init = res_s.init_params
+
+        for item in keys:
+            print(f'{item}: {bval[item]:.3f}, (init = {init[item].value:.3f}), bounds = [{init[item].min:.2f} : {init[item].max :.2f}] ')
+        print('')
     
 
 class NewFitModel(Model):

Analyser/ImagingAnalyser.py

@@ -167,6 +167,8 @@ class ImageAnalyser():
             center = self._center
         if span is None:
             span = self._span
+            
+        if not x in 
 
         x_start = int(center[0] - span[0] / 2)
         x_end = int(center[0] + span[0] / 2)
@@ -193,7 +195,14 @@ class ImageAnalyser():
                 dataSet[key].attrs['x_span'] = span[0]
                 dataSet[key].attrs['y_span'] = span[1]
 
-        return dataSet.isel(x=slice(x_start, x_end), y=slice(y_start, y_end))
+        res = dataSet.isel(x=slice(x_start, x_end), y=slice(y_start, y_end))
+        res = res.assign_coords(
+            {
+                'x': np.linspace(x_start, x_end - 1, span[0]),
+                'y': np.linspace(y_start, y_end - 1, span[1]),
+            }
+        )
+        return res
     
     def get_OD(self, imageAtom, imageBackground, imageDrak):
         """Calculate the OD image for absorption imaging.