imporve the BEC fitting

Analyser/FitAnalyser.py

@@ -86,33 +86,6 @@ 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): 
-    res = (1- ((x-centerx)/(sigmax))**2 - ((y-centery)/(sigmay))**2)**(3 / 2)
-    return amplitude * 5 / 2 / np.pi / max(tiny, sigmax * sigmay) * np.where(res > 0, res, 0)
-
-
-def polylog(power, numerator):
-    
-    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)
-    denominator = np.tile(denominator, (dataShape[0], 1))
-    denominator = denominator.T
-
-    data = numerator/ np.power(denominator, power)
-
-    return np.sum(data, axis=0)
-
-
-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)) ))
-
-
 def polylog_tab(pow, x, order=100):
     """Calculationg the polylog sum_i(x^i/i^pow), up to the order-th element of the sum
 
@@ -137,13 +110,49 @@ poly_tab = polylog_tab(2,x_int)
 polylog_int = CubicSpline(x_int, poly_tab)
 
 
-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, 
-                          amplitude=BEC_amplitude, sigmax=BEC_sigmax, sigmay=BEC_sigmay
-            ) + polylog2_2d(x=x, y=y, centerx=thermal_centerx, centery=thermal_centery, 
-                            amplitude=thermal_amplitude, sigmax=thermal_sigmax, sigmay=thermal_sigmay)
+def thermal(x, x0, amp, sigma):
+    res = np.exp(-0.5 * (x-x0)**2 / sigma**2)
+    return amp/1.643 * polylog_int(res)
+
+
+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/1.643  * polylog_int(np.exp( -((x-centerx)**2/(2 * sigmax**2))-((y-centery)**2/( 2 * sigmay**2)) ))
+
+
+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):
+    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):
@@ -510,7 +519,9 @@ class DensityProfileBEC2dModel(lmfit.Model):
         return lmfit.models.update_param_vals(params, self.prefix, **kwargs)
 
 
-    def fit(self, data_1d, **kwargs):
+    def fit(self, data, **kwargs):
+        
+        data_1d = data
 
         res = super().fit(data_1d, **kwargs)