2023-07-20 10:19:32 +02:00
{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Import supporting package"
],
"metadata": {
"collapsed": false
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 11,
2023-07-20 10:19:32 +02:00
"outputs": [],
"source": [
2023-07-20 20:34:19 +02:00
"import lmfit\n",
2023-07-20 10:19:32 +02:00
"import xarray as xr\n",
"import pandas as pd\n",
"import numpy as np\n",
"import copy\n",
"\n",
"import glob\n",
"\n",
"import xrft\n",
"import finufft\n",
"\n",
"from uncertainties import ufloat\n",
"from uncertainties import unumpy as unp\n",
"from uncertainties import umath\n",
"\n",
"from datetime import datetime\n",
"\n",
"import matplotlib.pyplot as plt\n",
"\n",
"\n",
"plt.rcParams['font.size'] = 18\n",
"\n",
"from scipy.ndimage import gaussian_filter\n",
"import matplotlib as mpl\n",
"from scipy.optimize import curve_fit\n",
"mpl.rc('xtick', labelsize=8)\n",
"mpl.rc('ytick', labelsize=8)\n",
"\n",
"from DataContainer.ReadData import read_hdf5_file, read_hdf5_global, read_hdf5_run_time, read_csv_file\n",
"from Analyser.ImagingAnalyser import ImageAnalyser\n",
"from Analyser.FitAnalyser import FitAnalyser\n",
"from Analyser.FitAnalyser import ThomasFermi2dModel, DensityProfileBEC2dModel, Polylog22dModel\n",
"from Analyser.FFTAnalyser import fft, ifft, fft_nutou\n",
"from ToolFunction.ToolFunction import *\n",
"\n",
"import time\n",
"\n",
"from ToolFunction.HomeMadeXarrayFunction import errorbar, dataarray_plot_errorbar\n",
"xr.plot.dataarray_plot.errorbar = errorbar\n",
"xr.plot.accessor.DataArrayPlotAccessor.errorbar = dataarray_plot_errorbar\n",
"\n",
"imageAnalyser = ImageAnalyser()\n",
"\n"
],
"metadata": {
2023-07-20 20:34:19 +02:00
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-07-20T08:25:33.998415800Z",
"start_time": "2023-07-20T08:25:33.964471900Z"
}
2023-07-20 10:19:32 +02:00
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 69,
2023-07-20 10:19:32 +02:00
"outputs": [],
"source": [
"# get center of thresholded image\n",
"\n",
"def calc_cen(thresh1):\n",
" \"\"\"\n",
" returns array: [Y_center,X_center]\n",
" \"\"\"\n",
" cen = np.zeros(2)\n",
" (Y,X) = np.shape(thresh1)\n",
"\n",
"\n",
" thresh1 = thresh1 /np.sum(thresh1)\n",
"\n",
" # marginal distributions\n",
" dx = np.sum(thresh1, 0)\n",
" dy = np.sum(thresh1, 1)\n",
"\n",
" # expected values\n",
" cen[0] = np.sum(dx * np.arange(X))\n",
" cen[1] = np.sum(dy * np.arange(Y))\n",
" return cen\n",
"\n",
"def calc_cen_bulk(thresh):\n",
" \"\"\"\n",
" returns array in shape of input, containing array with [Y_center,X_center] for each image\n",
" \"\"\"\n",
" shape = np.shape(thresh)\n",
" cen = np.zeros((shape[0], shape[1], 2))\n",
" for i in range(0, shape[0]):\n",
" for j in range(0, shape[1]):\n",
" cen[i,j] = calc_cen(thresh[i,j])\n",
" return cen\n",
"\n",
"def guess_BEC_width(thresh, center):\n",
" \"\"\"\n",
" returns width of thresholded area along both axis through the center with shape of thresh and [X_width, Y_width] for each image\n",
" \"\"\"\n",
" shape = np.shape(thresh)\n",
" BEC_width_guess = np.zeros((shape[0], shape[1], 2))\n",
"\n",
" for i in range(0, shape[0]):\n",
" for j in range(0, shape[1]):\n",
" BEC_width_guess[i, j, 0] = np.sum(thresh[i, j, round(center[i,j,1]), :])\n",
" BEC_width_guess[i, j, 1] = np.sum(thresh[i, j, :, round(center[i,j,0])])\n",
"\n",
" return BEC_width_guess\n",
"\n",
"\n",
"\n",
"def gaussian(x, x0, sigma, A):\n",
" return A * np.exp(-0.5 * (x-x0)**2 / sigma**2)\n",
"\n",
"def polylog(power, numerator, order = 15):\n",
"\n",
" dataShape = numerator.shape\n",
" numerator = np.tile(numerator, (order, 1))\n",
" numerator = np.power(numerator.T, np.arange(1, order+1)).T\n",
"\n",
" denominator = np.arange(1, order+1)\n",
" denominator = np.tile(denominator, (dataShape[0], 1))\n",
" denominator = denominator.T\n",
"\n",
" data = numerator/ np.power(denominator, power)\n",
"\n",
" return np.sum(data, axis=0)\n",
"\n",
2023-07-20 20:34:19 +02:00
"def thermal(x, x0, amp, sigma):\n",
" order = 15\n",
2023-07-20 10:19:32 +02:00
" res = np.exp(-0.5 * (x-x0)**2 / sigma**2)\n",
" return amp * polylog(2, res, order)\n",
"\n",
"def Thomas_Fermi_1d(x, x0, amp, sigma):\n",
" res = (1-(( x - x0 ) / sigma) **2) **3/2\n",
" return amp * np.where(res > 0, res, 0)\n",
"\n",
2023-07-20 20:34:19 +02:00
"def density_1d(x, x0_bec, x0_th, amp_bec, amp_th, sigma_bec, sigma_th):\n",
" return thermal(x, x0_th, amp_th, sigma_th) + Thomas_Fermi_1d(x, x0_bec, amp_bec, sigma_bec)"
2023-07-20 10:19:32 +02:00
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T14:17:21.289308300Z",
"start_time": "2023-07-20T14:17:21.240509100Z"
2023-07-20 10:19:32 +02:00
}
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 3,
2023-07-20 10:19:32 +02:00
"outputs": [],
"source": [
"# import data\n",
"img_dir = '//DyLabNAS/Data/'\n",
"SequenceName = \"Evaporative_Cooling\" + \"/\"\n",
"folderPath = img_dir + SequenceName + '2023/04/24'# get_date()\n",
"\n",
"\n",
"shotNum = \"0009\"\n",
"filePath = folderPath + \"/\" + shotNum + \"/2023-04-24_0009_Evaporative_Cooling_*0.h5\"\n",
"\n",
"\n",
"dataSet = read_hdf5_file(filePath, \"images/MOT_3D_Camera/in_situ_absorption\")\n",
"# flip the x and y axis\n",
"dataSet = swap_xy(dataSet)\n",
"\n",
"# get the scan axis name of the shot\n",
"scanAxis = get_scanAxis(dataSet)\n",
"\n",
"# calculate the absorption imaging\n",
"dataSet = imageAnalyser.get_absorption_images(dataSet)\n",
"\n",
"OD = dataSet[\"OD\"]\n",
"\n",
"OD_np = OD.to_numpy()"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T08:21:26.434937300Z",
"start_time": "2023-07-20T08:21:19.016336100Z"
2023-07-20 10:19:32 +02:00
}
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 65,
2023-07-20 10:19:32 +02:00
"outputs": [
{
"data": {
"text/plain": "<Figure size 1000x900 with 10 Axes>",
2023-07-20 20:34:19 +02:00
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA4gAAANiCAYAAAAjb+cqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9e3heVZ02jt9pk6YpJKUpJNC0TYQU2i+2QpWCFi0wdTjMixzkoCiCDlB8dVR0PCCOyOgPZEQcRkZlHE+visIlhxFR1CJUrQytFKRoqw2SQgMkQKAJNpSk7N8fe9/Putdnr/0kaXMqrPu6cuV5nr332ut8+Hzuda+KJEkSRERERERERERERERERLziMWm8IxARERERERERERERERExMRAXiBEREREREREREREREREA4gIxIiIiIiIiIiIiIiIiIkNcIEZEREREREREREREREQAiAvEiIiIiIiIiIiIiIiIiAxxgRgREREREREREREREREBIC4QIyIiIiIiIiIiIiIiIjLEBWJEREREREREREREREQEgLhAjIiIiIiIiIiIiIiIiMiw2y4Qv/3tb6OiogItLS0v63cS7e3tqKioQEVFBdrb28f8+Z0F33n33XeP2TsjRgexzY3t8zuL2OYihorYpsf2+Z1FbNMRERFjjd12gVgOt956Kz7zmc/g1ltvHe+ovGzx7//+7/jMZz6DBx54YLyjMuFx4YUXlgb48ZgUjQQefvhhfPvb3wYAbN68Gfvssw+OPfZY3HTTTQB2vc2tW7cO73znOzF79mxUV1djv/32wymnnII//elPQ3r+rrvuwimnnIL99tsP1dXVmD17Nt75zndi3bp1g773q1/9Ks4//3wsXrwY1dXVE7acYpsbOl4ubW7FihV41atehalTp+ba3K5i3bp1uO666wCkbZpt7le/+hWAcJvu7+/H17/+daxYsQKHH3445syZg5qaGkybNg0HHHAAzjrrLKxcuXJI79+2bRuOPfZYNDQ0YOrUqXjVq16FFStWoK2tbdhp+d3vfofJkycPupDavn176fOJJ56IvfbaC1VVVdhnn31w9NFH48tf/jK2bds27PfvLGKbHjpimy6P7du348tf/jKOPPJIzJgxA1OnTkVLSwvOO++8suNokiS455578KlPfQpHHXUUGhsbUVVVhenTp+O1r30tLr74YnR0dJR9d0tLS6lsiv6OPPLIwnjfdttteP/734/Xve51494mIwTJbopvfetbCYCkubk5d+2cc85JACTnnHPOmL1ztPHII48kABIAySOPPDLmz1s0NzcnAJJvfetbZe876KCDkoMOOii59957d/mduyN+9atfJRUVFaW8H4+6s6u4/fbbk2nTppXSUFFRkUyaNKn0/d3vfvcutbmvf/3rSWVlZSm86dOne3k2WL5deumlXtymT59e+l5ZWZl8/etfL3yW9dj+NTc3xza3m+Ll2Obq6upybe6ll17a6fBtm6uoqPDy7NJLLw226aeeesprJxUVFcmMGTOSyZMne7+fc845SX9/f/Dd3/zmN717J02alNTV1ZW+T5s2Lbn99tuHnJa+vr5k/vz5Xph33XVX8N6lS5d691VWVnr9BYCktbU12bRpU/D52KbHB7FNl8cTTzyRHHrooaWwqqqqkhkzZpS+V1dXJ9///veDz37uc5/Ltem99trLy++6urrkf/7nfwrfz3pcV1eXNDY2Bv/e8pa3BJ9dvnz5LrXJiNHDy9KDGDFxsHHjRmzcuBFLliwZ76iMObZt24bzzz8flZWVeN3rXjfe0dkpPPLIIzjjjDOwbds2tLa2AgDmzp2LrVu34tOf/jQA4Fvf+hYeeuihnQr/nnvuwYUXXoiBgQGcfPLJeOyxx/Dcc8/hqaeewooVK0r3FVkQb7zxRlx22WUAgBUrVuCpp57Cc889h8ceewwnn3wyBgYGcOGFF+Kee+4JPj9lyhQccsgheM973oNrr70WZ5999k6lYyIhtrmXT5tbunQp/vznP2Pr1q25NveFL3xhp8LXNnfooYcCSNu0trnLLrssSKGsrq7GP/3TP+GGG25Ae3s7tm/fju7ubrz44otYv3493va2twEAvvOd7+Cqq64Kvv9nP/tZ6fOll15aStvGjRvxhje8Adu2bcMZZ5yBRx55ZEjpueyyy0rPDoaBgYHS55/85CfYvn07nnvuOXR3d+PKK6/EtGnT0NbWhuOPPx4vvPDCkN4/FohtOrbpIiRJgre+9a24//77UVNTg69//evo6elBd3c3Hn/8cbzrXe/C9u3bce655+K+++7LPd/f34+6ujpceOGF+NWvfoW//e1vePbZZ/G3v/0NN910E+bOnYuenh6cfvrp2LBhQ9m4XHPNNXjyySeDf//zP/8TfKa/vx/Nzc34zGc+g3Xr1u02bfIVgfFeoe4sogdxbJ+3GKrl85WMD33oQwmA5JJLLinVyd3N8vnOd74zAZDsu+++yX/+53/m0nDBBReULJY70+aOPPLIBECycOHC5MUXX8xdf/WrX50ASCZPnpwMDAx41wYGBkr18Ljjjss9u3379tLzRx55ZPD9Nkx6I3dnD+IrGS+3Nvfss8/mrrPN1dXVJd3d3cMOX9vcf//3f+fy6Nhjj00AJHvsscew2/RLL72UvOENbyhZ/S26u7uTqVOnJgCSPffcM3h93333TQAk73znOwd933333ZdUVlYmra2tyc9+9rNBPYg33nhj2Tb5gx/8oHT9u9/9bu56bNNjj9imy+O2224r1cl///d/D95zxBFHJACSY445Jnft/vvvL/vOv/71r0lNTU0CIPnHf/zH4D27Uo9/85vfFLINkmTwNhkxehj1BeJdd92VnHbaacmsWbOSKVOmJDNnzkyOOeaY5Jvf/GZucmZxzz33JCeddFIyc+bMZOrUqcmBBx6YfPKTn0x6e3uDi7W77rorSBfTv6KBYyiw7/z973+fnH766cm+++6bTJkyJXnVq16VXHTRRYM28La2tuTCCy9MWltbk6lTpya1tbXJoYcemlx22WXJ1q1bg88MZWDasmVLcsEFFySzZ89OpkyZkjQ1NSXnnntusmnTphEb2JTOV/SnKMp3G5/29vbkvPPOS+bMmZNUV1cn+++/f3LJJZckzz//fOmZ9evXJ+94xzuS2bNnJ9XV1Ulra2vy2c9+NriwsO/64Ac/mPx//9//l+yxxx5JTU1NctBBByUf+MAHks2bN+90XpTDPffck0yaNCk58MADk76+vjEd2EayzZFm8slPfjLY5rTzHm6be/jhh0v3fOc73wne8/GPf7x0z69+9SuvzXFRCiD5yU9+Enz+29/+dumes846a9A2N9wFYmxzYbzS2txI4fnnny9Nxi677LLS77ZNsxzf/e53D6tNV1dXl5697rrrgm36S1/60i6NoxdddFECIJk6dWru2je+8Y1SGE1NTUmS5MfRvfbaq/S81kWLDRs2JPX19QmAZMqUKR59r6g/GKxN7tixo7Qw3nPPPWObNu+KbXr4KGrTCi3Hb37zm8MK/8ILL0yA1KBTVIY33XRTKfydKasTTjghAZAsWrQoeH00DR07duxIamtrEwDJBz7wgREPP6IYo7pA5EABOF6z7lU45phjkp6enuCz3/jGNzx+9vTp00sD4/z585Orr74611GsXr06aWxsLFkop06dmuNBr169eqfTo4Pp97///dIEdfr06V5cDz744KS3tzcYxg033OAN0rW1td73OXPmJH/6059yzw02MN13330e57ympibZc889EyC1St1www0jMrB94QtfSBobG0vpDXHOFUMZ2G666abSpKCurs6rI2984xuTF198MfnJT35SmgDYPWpnnnlmYXy/973veflbXV1d6qyZ/z//+c93Oj9CeOGFF5IFCxYkFRUVyd13350kSTJmA9tItjlOlAAkLS0thW1Owx9Om/va175Weq6zs7MwTrznxBNP9Nqc1oGiNtfZ2ZmbeJVrc8NZIMY2F8Yrrc2NJO64445SPq1ZsyZJknCb1vo8nDat5TBv3rxgm/71r3/tlfVw2vSOHTuSww8/vNQmLd7
2023-07-20 10:19:32 +02:00
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"imageAnalyser.center = (960, 890)\n",
2023-07-20 20:34:19 +02:00
"imageAnalyser.span = (150, 150)\n",
2023-07-20 10:19:32 +02:00
"imageAnalyser.fraction = (0.1, 0.1)\n",
"\n",
"dataSet_cropOD = imageAnalyser.crop_image(dataSet.OD)\n",
"dataSet_cropOD = imageAnalyser.substract_offset(dataSet_cropOD).load()\n",
"cropOD = dataSet_cropOD.to_numpy()\n",
"dataSet_cropOD.plot.pcolormesh(cmap='jet', vmin=0, col=scanAxis[0], row=scanAxis[1])\n",
"plt.show()"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T14:15:36.973932600Z",
"start_time": "2023-07-20T14:15:31.476880400Z"
2023-07-20 10:19:32 +02:00
}
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 47,
2023-07-20 10:19:32 +02:00
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
2023-07-20 20:34:19 +02:00
"fitting time: 112.75005340576172 ms\n"
2023-07-20 10:19:32 +02:00
]
}
],
"source": [
"# from opencv import moments\n",
"shape = np.shape(cropOD)\n",
"sigma = 0.4\n",
"blurred = gaussian_filter(cropOD, sigma=sigma)\n",
"\n",
"thresh = np.zeros(shape)\n",
"for i in range(0,shape[0]):\n",
" for j in range(0, shape[1]):\n",
" thresh[i,j] = np.where(blurred[i,j] < np.max(blurred[i,j])*0.5,0,1)\n",
"\n",
"center = calc_cen_bulk(thresh)\n",
"\n",
"BEC_width_guess = guess_BEC_width(thresh, center)\n",
"\n",
"X_guess_og = np.zeros((shape[0], shape[1], shape[3]))\n",
"# Y_guess_og = np.zeros((shape[0], shape[1], shape[2]))\n",
"\n",
"for i in range(0, shape[0]):\n",
" for j in range(0, shape[1]):\n",
" X_guess_og[i,j] = np.sum(cropOD[i,j,round(center[i,j,1] - BEC_width_guess[i,j,1]/2) : round(center[i,j,1] + BEC_width_guess[i,j,1]/2) , :], 0) / len(cropOD[i,j,round(center[i,j,1] - BEC_width_guess[i,j,1]/2) : round(center[i,j,1] + BEC_width_guess[i,j,1]/2),0])\n",
"\n",
" # Y_guess_og[i,j] = np.sum(cropOD[i,j, :, round(center[i,j,0] - BEC_width_guess[i,j,0]/2) : round(center[i,j,0] + BEC_width_guess[i,j,0]/2)], 1) / len(cropOD[i,j,0,round(center[i,j,0] - BEC_width_guess[i,j,0]/2) : round(center[i,j,0] + BEC_width_guess[i,j,0]/2)])\n",
"\n",
"#[nr images x, nr images y, X / Y, center / sigma]\n",
"x = np.linspace(0,149,150)\n",
2023-07-20 20:34:19 +02:00
"y = np.linspace(0,149, 200)\n",
2023-07-20 10:19:32 +02:00
"\n",
"popt = np.zeros((shape[0], shape[1], 6))\n",
"\n",
"p0 = np.ones((shape[0], shape[1], 6))\n",
"\n",
"max = np.zeros((shape[0], shape[1]))\n",
"\n",
"for i in range(0, shape[0]):\n",
" max[i] = np.ndarray.max(X_guess_og[i],axis=1)\n",
"\n",
"\n",
"p0[:, :, 0] = center[:, :, 0] # center BEC\n",
"p0[:, :, 1] = center[:, :, 0] # center th\n",
"p0[:, :, 2] = 0.7 * max # amp BEC\n",
"p0[:, :, 3] = 0.3 * max # amp th\n",
"p0[:, :, 4] = BEC_width_guess[:, :, 0] # sigma BEC\n",
"p0[:, :, 5] = BEC_width_guess[:, :, 0] * 3 # sigma th\n",
"\n",
"start = time.time()\n",
"for i in range(0, shape[0]):\n",
" for j in range(0, shape[1]):\n",
" popt[i,j], pcov = curve_fit(density_1d, x, X_guess_og[i,j] , p0[i, j], nan_policy='omit')\n",
"stop = time.time()\n",
"\n",
"print(f'fitting time: {(stop-start)*1e3} ms')"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T13:44:57.754446Z",
"start_time": "2023-07-20T13:44:57.598621Z"
2023-07-20 10:19:32 +02:00
}
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 43,
2023-07-20 10:19:32 +02:00
"outputs": [
{
"data": {
"text/plain": "<Figure size 1500x1000 with 9 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAABMAAAAMpCAYAAADipk9SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3hUVfrA8e+dnjIz6QWSEHqXIqCCAmJviGJdGygsrrq23dVVV3f156pr74qK7LquvWJvWMBGE5AiHZKQ3mYyk0y/vz/uJBJpKTOZlPfzPPPIzNy5950xmZPznnPeo6iqqiKEEEIIIYQQQgghRDeli3UAQgghhBBCCCGEEEJEkyTAhBBCCCGEEEIIIUS3JgkwIYQQQgghhBBCCNGtSQJMCCGEEEIIIYQQQnRrkgATQgghhBBCCCGEEN2aJMCEEEIIIYQQQgghRLcmCTAhhBBCCCGEEEII0a0ZYh1Aa4RCIYqLi7FarSiKEutwhBCiy1NVlbq6Onr16oVOJ2MiIG2NEEJEmrQ1zUk7I4QQkdXSdqZLJcCKi4vJzc2NdRhCCNHtFBYWkpOTE+swOgVpa4QQIjqkrdFIOyOEENFxsHamSyXArFYroL0pm80W42iEEKLrczqd5ObmNn2/CmlrhBAi0qStaU7aGSGEiKyWtjNdKgHWOEXYZrNJYyGEEBEkSzB+JW2NEEJEh7Q1GmlnhBAiOg7WzsgifCGEEEIIIYQQQgjRrUkCTAghhBBCCCGEEEJ0a5IAE0IIIYQQQgghhBDdmiTAhBBCCCGEEEIIIUS3JgkwIYQQQgghhBBCCNGtSQJMCCGEEEIIIYQQQnRrkgATIsq8gSBvrSqitt4X61CEEEJ0c19vrmBVQU2swxBCCNHFfbe1kmU7qmMdhhARJQkwIaLsjvc2cP1ra3jkiy2xDkUIIUQ3Vunycum/l3PJgmUEgqFYhyNEp3X11VeTn5+PoiisXr16n8csXLiQ0aNHN93S0tI488wzAdi5cyd6vb7Z89u2bevAdyBEdNXW+7hk4TJmLVyGNxCMdThCRIwh1gEI0Z1tq3DxyvJCAH4ucsQ4GiGEEN3Zzko3wZBKnTfAtgo3g7OssQ5JiE7prLPO4oYbbuDII4/c7zGzZ89m9uzZTfdHjBjBBRdc0HTfarXuN3kmRFe3vtiJP6jiDwYpqmmgf3pirEMSIiIkASZEFN3/ySaCIRWALeUuVFVFUZQYRyWEEKI72l3b0PTvDSUOSYAJsR+TJ09u1fE//vgj5eXlTJ8+vU3X83q9eL3epvtOp7NN5xGio2ws+fVndFeVWxJgotuQJZBCRMlPBTV8tK4UnQKKAo4GP5UuqQMmxIHIshQh2q6oZo8EWLF0sIWIlAULFnDRRRdhNBqbHnO73YwfP56xY8dyxx13EAzuf5nY3Xffjd1ub7rl5uZ2RNhCtNmGZgmw+hhGIkRkSQJMiCh5+6fdAMwY3Zs+KfEAbCmvi2VIQnR6Z511FkuXLqVPnz77PWb27NmsXr266ZaVlbXPZSmNt/79+3dE6ELE3J4zwNZLAkyIiHC73bzyyitcdtllTY9lZ2eze/duli9fzueff86SJUt44IEH9nuOm266CYfD0XQrLCzsiNCFaLONJb/2WSQBJroTSYAJESXV1ZWcoVvCydZtjEjTftW2lrtiHJUQndvkyZPJyclp8fHtXZYC2tIUp9PZ7CZEV7R7jxlgu4pLURtkN0gh2uv1119n+PDhDBs2rOkxs9lMRkYGACkpKVx66aUsWbJkv+cwm83YbLZmNyE6K18gxNY9Bu13VrljGI0QkSUJMCGiZErpCzxkeopjl13Ko7vOYIZuKVvKJAEmRCS1d1kKyNIU0X00zgBLoIFXQ38i9NgEaKiNbVBCdHELFixoNvsLoLy8HL/fD2iDKG+99RZjxoyJRXhCRNy2Chf+oNp0v0BmgIlupFUJMKnNIkTLDfGuBSBoTEBHkOn672QGmBARFIllKSBLU0T3oKpq0wywuebPyVEq0deXw6oXYhyZEJ3PvHnzyMnJoaioiBNOOIEBAwYAMGfOHBYtWtR03KZNm1i9ejXnnntus9cvXbqUMWPGMGrUKMaOHUtWVha33HJLh74HIaKlsQB+TnIcAIU19QSCoViGJETEtGoXSNkyWIiW8XgaGKTuAAW8x99P/Ad/YIRuJ1skASZExLR0WcpLL73EDTfcsN/zmM1mzGZz1OMVIppq6v00+IPE42GO/gNo7Kv8OB8OvwL0svG3EI3mz5+/z8efe+65ZvcHDx5MXd3e9VvPPPPMpgF+IbqbxgTY0YMzeHVFIb5AiBKHh9xwTWMhurJWzQDr6NosUpdFdFXV21djVgLUqonEHTIdVdGRodSiuEqpcctOkEJEgixLEeJXjbO//hC/mMSQkx2hTJz6JHAWwcZ3YxucEEKILqOxAP7wXjbywkkvqQMmuouo1gCTLYNFT+XZuQyAzYaBKOZElLTBAIzQ7WRrhcwCE2J/ZFmKEG2zu7YeE34uVrXfk8cCZ/C67kTtye8eB1U9wKuFEEIIbTl94wywodk28lO1BJjsBCm6i6jNh2+szfLDDz80PdZYmyUjI4Pq6mrOPfdcHnjggf0uTbnpppu4/vrrm+47nU5JgokuQV+yCoCi+KFMAMgeBRUbGaHsYEuZi/H5KTGNT4jOSpalCNE2RTUNDFSKsKtOQpZk3vVM4us6F7Pj30RXvAqqtkHagFiHKYQQohOrcHmpcvvQKTA4y0qf1AQAdskMMNFNRG0GmGwZLHqyxKqfAahJPkR7IHsUACN1O9hSvnenXQghhGiP3bUNDFR2A6DLHMYFR/SjCjvrA+GBw4qNMYxOCCFEV1BS6wEg02bBYtTTJ7VxCaTMABPdQ9QSYFKbRfRY3jqS67cD4MvQEl+NCbDhup3srJQRFCGEEJFVVNPAQF2Rdid9CH8/bTgnjchis9oLAE+xJMCEEEIcWHW9Vqs4JcEE0DQDrEASYKKbaFUCTGqzCNECxavRoVKkpmFL7609ljUSgN5KFR5HeQyDE0II0R3trmlgUHgGGOlD0OsUHj5vNJWWfABqC36OXXBCCCG6hGrXHgmw5c8xvPw9AHZVuwmFpJak6PpaVQNMarMI0QLFWv2vNaF+9LLHaY9ZbHjtfTE7dpDi/AU4NXbxCSGE6HZ21zYwQGmcAaZtvGI26NFnDIbdoKvaHMPohBBCdAU14Rlgh6k/wwd/Ig0YrLuXTf4cKl1eMmyW2AYoRDtFdRdIIXqkkrUA/BzqR5b910ZCzdKWQeb5thAIhmISmhBCiO7H5Q3gaXCTp4RnGKcPaXrOmjMcALt7J4Sk7RFCCLF/1W4foDKjekHTY1eYPwagMjw7TIiuTBJgQkRYqLYQgF1q5q8zwABTjlbvbpiyUxoQIYQQEVNc20A/pQS9okJcMiRmND2X2384PlWPWfWAozCGUQohhOjsaup9nKBbQU79BtAZAThJ/YZ0asLJMSG6NkmACRFhIYe2BKVan44t7tdVxrrUfgD0Viopc3piEpsQQojup7LO27QDJOlDQFGanhuWk8IONRsA9+4NsQhPCCFEF1HjauBPhte0O5OugdzDMBHgEsOnTQXyhejKJAEmRCSFguhcZQCotl4oe3RCsGk7cWUqNZIAE0IIETHV9T4G6JrX/2qUFG+i2JALQMWOtR0dmhBCiC7E5tjEIN1u/IYEmPhH7QZcqP8ch8MZ4+iEaD9JgAkRSe4KdGqAgKrDlJTV/LnGBBg1VDhlK2EhhBCRUe32NZ8B9htuu7Zrt7dEZoAJIYTYv1T3VgDqU0dAXBIMPhm33k6S4kap3BTb4ISIAEmACRFJTq0DUk4SmfbE5s8lZhJEj0EJ4aoqjkFwQgghuqNqt49BTTtA7p0AM2Roj5lqtnZkWEIIIbqYbO8OAELpQ7UHdHqqE/oDYK6R3YRF1ycJMCEiyakltkrVFHol/WabYJ2eBlMqAL7q3R0dmRBCiG7KUeeij6Itv99XAiwlfwQAqQ07QVU7MDIhhBBdRTCkkhfcBYA+a3jT426bNovYWrctJnEJEUmSABMikhxaYqtETSH
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fsize=(15,10)\n",
"vmax = 1.5\n",
"fig, ax = plt.subplots(shape[0],shape[1],figsize=fsize)\n",
"for i in range(0, shape[0]):\n",
" for j in range(0, shape[1]):\n",
" lab = f\"A$_{{BEC}}$ = {popt[i,j,0]:.1f} \\n A$_{{th}}$ = {popt[i,j,1]:.1f} \"\n",
" ax[i,j].plot(x, X_guess_og[i,j])\n",
" ax[i,j].plot(x, density_1d(x, *popt[i,j]), label = lab, zorder=3)\n",
" ax[i,j].plot(x, thermal(x, popt[i,j,1], popt[i,j, 3], popt[i,j, 5]))\n",
"\n",
"\n",
" #ax[i,j].legend(fontsize=10)\n",
" ax[i,j].set_facecolor('#FFFFFF')\n",
"plt.show()\n"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T13:41:56.067625400Z",
"start_time": "2023-07-20T13:41:54.555950300Z"
2023-07-20 10:19:32 +02:00
}
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 74,
2023-07-20 10:19:32 +02:00
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
2023-07-20 20:34:19 +02:00
"total time: 308.0174922943115 ms\n"
2023-07-20 10:19:32 +02:00
]
}
],
"source": [
2023-07-20 20:34:19 +02:00
"# from opencv import moments\n",
"start = time.time()\n",
"shape = np.shape(cropOD)\n",
"sigma = 0.4\n",
"blurred = gaussian_filter(cropOD, sigma=sigma)\n",
"\n",
"thresh = np.zeros(shape)\n",
"for i in range(0,shape[0]):\n",
" for j in range(0, shape[1]):\n",
" thresh[i,j] = np.where(blurred[i,j] < np.max(blurred[i,j])*0.5,0,1)\n",
"\n",
"center = calc_cen_bulk(thresh)\n",
"\n",
"BEC_width_guess = guess_BEC_width(thresh, center)\n",
"\n",
"X_guess_og = np.zeros((shape[0], shape[1], shape[3]))\n",
"# Y_guess_og = np.zeros((shape[0], shape[1], shape[2]))\n",
"\n",
"for i in range(0, shape[0]):\n",
" for j in range(0, shape[1]):\n",
" X_guess_og[i,j] = np.sum(cropOD[i,j,round(center[i,j,1] - BEC_width_guess[i,j,1]/2) : round(center[i,j,1] + BEC_width_guess[i,j,1]/2) , :], 0) / len(cropOD[i,j,round(center[i,j,1] - BEC_width_guess[i,j,1]/2) : round(center[i,j,1] + BEC_width_guess[i,j,1]/2),0])\n",
"\n",
" # Y_guess_og[i,j] = np.sum(cropOD[i,j, :, round(center[i,j,0] - BEC_width_guess[i,j,0]/2) : round(center[i,j,0] + BEC_width_guess[i,j,0]/2)], 1) / len(cropOD[i,j,0,round(center[i,j,0] - BEC_width_guess[i,j,0]/2) : round(center[i,j,0] + BEC_width_guess[i,j,0]/2)])\n",
"\n",
"#[nr images x, nr images y, X / Y, center / sigma]\n",
"x = np.linspace(0,shape[3],150)\n",
"y = np.linspace(0,shape[2], 150)\n",
"\n",
"popt = np.zeros((shape[0], shape[1], 6))\n",
"\n",
"p0 = np.ones((shape[0], shape[1], 6))\n",
"\n",
"max = np.zeros((shape[0], shape[1]))\n",
"\n",
"for i in range(0, shape[0]):\n",
" max[i] = np.ndarray.max(X_guess_og[i],axis=1)\n",
"\n",
"# Fitting x without math constr\n",
"fitmodel = lmfit.Model(density_1d,independent_vars=['x'])\n",
"\n",
"result_x = []\n",
"\n",
"for i in range(0, shape[0]):\n",
" temp_res = []\n",
" for j in range(0, shape[1]):\n",
" params = lmfit.Parameters()\n",
" params.add_many(\n",
" ('x0_bec', center[i,j,0], True,0, 200),\n",
" ('x0_th',center[i,j,0], True,0, 200),\n",
" ('amp_bec', 0.7 * max[i,j], True, 0, 1.3 * np.max(X_guess_og[i,j])),\n",
" ('amp_th', 0.3 * max[i,j], True, 0, 1.3*np.max(X_guess_og[i,j])),\n",
" ('sigma_bec',BEC_width_guess[i,j,0], True, 0, 50),\n",
" ('sigma_th', 3*BEC_width_guess[i,j,0], True, 0, 50)\n",
" )\n",
" res = fitmodel.fit(X_guess_og[i,j], x=x, params=params)\n",
" temp_res.append(res)\n",
" result_x.append(temp_res)\n",
"stop = time.time()\n",
"\n",
"print(f'total time: {(stop-start)*1e3} ms')"
2023-07-20 10:19:32 +02:00
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T14:19:04.018601600Z",
"start_time": "2023-07-20T14:19:03.703507500Z"
2023-07-20 10:19:32 +02:00
}
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 72,
2023-07-20 10:19:32 +02:00
"outputs": [
{
2023-07-20 20:34:19 +02:00
"data": {
"text/plain": "<Figure size 1500x1000 with 9 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAABMAAAAMpCAYAAADipk9SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd5hTVfrA8e9NmWRaMr3AzDD0XhUbil2woYJ1baAoupZd3VXXdddVd227q6ur64qK+LNiARV7L4AiiHSkDDBM7zNJJplJu/f3x82MzgIywySTKe/nefJIkpt730TIyXnPOe9RNE3TEEIIIYQQQgghhBCilzJEOwAhhBBCCCGEEEIIISJJEmBCCCGEEEIIIYQQoleTBJgQQgghhBBCCCGE6NUkASaEEEIIIYQQQgghejVJgAkhhBBCCCGEEEKIXk0SYEIIIYQQQgghhBCiV5MEmBBCCCGEEEIIIYTo1UzRDqAjVFWlrKyMxMREFEWJdjhCCNHjaZqGy+WiX79+GAwyJgLS1gghRLhJW9OWtDNCCBFe7W1nelQCrKysjNzc3GiHIYQQvU5xcTE5OTnRDqNbkLZGCCEiQ9oanbQzQggRGQdqZ3pUAiwxMRHQ35TNZotyNEII0fM5nU5yc3Nbv1+FtDVCCBFu0ta0Je2MEEKEV3vbmR6VAGuZImyz2aSxEEKIMJIlGD+RtkYIISJD2hqdtDNCCBEZB2pnZBG+EEIIIYQQQgghhOjVJAEmhBBCCCGEEEIIIXo1SYAJIYQQQgghhBBCiF5NEmBCCCGEEEIIIYQQoleTBJgQQgghhBBCCCGE6NUkASaEEEIIIYQQQgghejVJgAkRYd5AkCU/lNDg8UU7FCGEEL3cV9ur+aGoPtphCCGE6OG+Kahh1e66aIchRFhJAkyICLvnnS3c/Np6Hv1sR7RDEUII0YvVNHq54rnVXL5gFYGgGu1whOi2brzxRvLz81EUhXXr1u3zmIULFzJhwoTWW1paGjNnzgSgsLAQo9HY5vmdO3d24TsQIrIaPD4uX7iK2QtX4Q0Eox2OEGFjinYAQvRmO6sbWbS6GICNJY4oRyOEEKI3K6xxE1Q1XN4AO6vdDM9KjHZIQnRL5557LrfeeitHH330fo+ZM2cOc+bMab0/ZswYLr744tb7iYmJ+02eCdHTbS5z4g9q+INBSuqbGJyeEO2QhAgLSYAJEUH//GgbQVUDYEdVI5qmoShKlKMSQgjRG5U2NLX+eUu5QxJgQuzH1KlTO3T8d999R1VVFTNmzDio63m9Xrxeb+t9p9N5UOcRoqv8WP7T39E9tW5JgIleQ5ZAChEha4vq+WBTBQYFFAUcTX5qGqUOmBC/RJalCHHwSup/lgArkw62EOGyYMECLr30Usxmc+tjbrebyZMnM2nSJO655x6Cwf0vE7v//vux2+2tt9zc3K4IW4iDtqVNAswTxUiECC9JgAkRIW+uLQXg7An9GZASB8COKlc0QxKi2zv33HNZvnw5AwYM2O8xc+bMYd26da23rKysfS5LabkNHjy4K0IXIup+PgNssyTAhAgLt9vNokWLuPLKK1sfy87OprS0lNWrV/Ppp5+ybNkyHnroof2e4/bbb8fhcLTeiouLuyJ0IQ7aj+U/9VkkASZ6E0mACREhZaGOyKH5KQzJ0JehFFQ1RjMkIbq9qVOnkpOT0+7jO7ssBfSlKU6ns81NiJ6o9OczwMqdaJoWxWiE6B1ef/11Ro8ezahRo1ofs1gsZGRkAJCSksIVV1zBsmXL9nsOi8WCzWZrcxOiu/IFVAp+NmhfWOuOYjRChFeHEmCyNEWI9itraAYg225lSIa+bn5HpSTAhAinzi5LAVmaInqPn88Aa/D4KXc0RzEaIXqHBQsWtJn9BVBVVYXf7wf0QZQlS5YwceLEaIQnRNjtrG7EH/xpAKVIZoCJXqRDCTBZmiJE+1U4QwmwJCtDQwkwmQEmRPiEY1kKyNIU0TtomtY6AyzRou9xJHXAhNi3efPmkZOTQ0lJCdOmTWPIkCEAzJ07l6VLl7Yet23bNtatW8cFF1zQ5vXLly9n4sSJjB8/nkmTJpGVlcUdd9zRpe9BiEhpKYCfkxwLQHG9h0BQjWZIQoRNh3aB7OodU4ToqZr9QercesH7bFssvoDeaOyQBJgQYdPeZSkvv/wyt956637PY7FYsFgsEY9XiEiq9/hp8uuzHacOT+e9DeVsLnNy0qjMKEcmRPczf/78fT7+zDPPtLk/fPhwXK6967fOnDmzdYWLEL1NSwLs+OEZvPp9Mb6ASrmjmdxQTWMherKI1gDr7NIUqcsieqqK0LKTWLMRW6ypdevgmkYv9W7ZCVKIcJBlKUL8pGX2V3qihYm5SUxStuPbtf+aREIIIcS+tBTAH93PRl4o6SV1wERvEbEEWDiWpkhdFtFT6XVXNMbZXCgBL/EWE/2T9GnEBdUyC0yI/ZFlKUIcnNIGvUZL/6RYJiY2sCjmr9xU9juolVqrQggh2kfTtNYZYCOzbeSn6gkw2QlS9BYRS4CFY8cUqcsieqpyRxOXGz/mVfdcuK8fPHEUp9l2A1IIX4hfMn/+fEpKSggEAlRWVlJQUADoy1J+vpy+ZVlKYmJim9fPnDmTTZs2sX79ejZv3sxjjz0myxtFn1ASmgGWkxzL2IIniVGCmFApeveBKEcmhBCip6hu9FLr9mFQYHhWIgNS4wHYIzPARC8RsQRYOJamyJbBoqcqdzQzzbBav6MFoWozM3367JUdVXvXkhBCCCE6o2UHyLHWSmI2v976eOauN1m98cdohSWEEKIHKQ/tYp9ps2I1GxmQ2rIEUmaAid6hQwkwWZoiRPtUNrgYbwgtOzn5HgDymrcCUFgjIyhCCCHCq2UG2LSqhaCpaMNOZZdlFBbFz7o3HsDh8Uc5QiGEEN1dnUevVZwSHwPQOgOsSBJgopfo0C6QsmOKEO1jrN5KvOLFZ0og5pDZ8MlfiG+uIA0HVS6ZySiEECK8SuubSMNBfsVHACgn/Inc2l3w+qVcwMes3V3KsaPzoxukEEKIbq2usW0CrLUGWJ0bVdUwGJSoxSZEOER0F0gh+qq0hg0ANKaOA6sd0ocDMM6wk0qnN5qhCSGE6IVKG5oYYSjS76QOhawxmEeeQYMpDZvioXbH6ugGKIQQoturD80AS47TE2D9kmIxKNDsV6lplD6M6PkkASZEBOQ1bQYg2H+y/kC/SQCMN+yk1u0lEFSjFZoQQoheptEbwNHkZ7gS2iwoY6T+X4MBp10fgGkq3Ryl6IQQQvQUde62M8DMRkPrn2tCs8OE6MkkASZEmDX7g4wObgMgbtCR+oP99QTYBMMuNE0aECGEEOFTFiqAP8Zcqj+Q8dMO3MZQMsxSv73L4xJCCNGz1P+8BtjOL6BwRWsCrCU5JkRPJgkwIcKsqqKUQYYKAOIGHaE/2DoDbBegUelsjlJ0Qgghepsal74sZYShJQE2svW55PxxAGR7C3E0SSF8IYQQ+9eS5BoU3AkvnAPPn8VQS4P+nEcSYKLnkwSYEGHm3rUSgD2GHJS4ZP3BrDFgMJOEixylWhJgQgghwqbO40NBZaAWqgH2sxlgcTljABhmKGFLmTMa4QkhhOgh9ASYxhE7HgE0UP3M8r2tPyc1wEQvIAkwIcLMUKIXGi60/jQCj8kCmaMBGK/sosolDYgQQojwqHP76K/UYNWawRgDKYN+ejJNrwGWrjgoKCyMToBCCCF6hDq3j6mGDaRVfwvoOz4e7XyPJFyyBFL0CpIAEyLMYur0+l+1ttFtn+h/CKDvBFklM8CEEEKESZ3b91MB/LRhYDT99KQlAYelHwCOPRujEJ0QQoieosHt5Q+mRfqdI34NWWOJUZu5zPgJtZIAE72AJMCECDOrp0z/Q/KAtk9kjwdghFJMpVNmgAkhhAiPOrePYcre9b9a+FOHAaBW/diVYQkhhOhBgqpGdnMBowx7UM1xMPX3MOW3AMw2fYizsTG6AQoRBpIAEyLMEr16AfyY1Ly2TyTlApCt1FLlkhlgQgghwqPW7WOYITQDbB8JsNh+eh2wZPdOmnzBrgxNCCFED+F
},
"metadata": {},
"output_type": "display_data"
2023-07-20 10:19:32 +02:00
}
],
"source": [
2023-07-20 20:34:19 +02:00
"fig, ax = plt.subplots(shape[0],shape[1],figsize=fsize)\n",
"for i in range(0, shape[0]):\n",
" for j in range(0, shape[1]):\n",
" #lab = f\"A$_{{BEC}}$ = {popt[i,j,0]:.1f} \\n A$_{{th}}$ = {popt[i,j,1]:.1f} \"\n",
" bval = result_x[i][j].best_values\n",
" ax[i,j].plot(x, X_guess_og[i,j])\n",
" ax[i,j].plot(x, density_1d(x, **result_x[i][j].best_values), label = lab)\n",
" ax[i,j].plot(x, thermal(x, bval['x0_th'], bval['amp_th'], bval['sigma_th']))\n",
2023-07-20 10:19:32 +02:00
"\n",
2023-07-20 20:34:19 +02:00
"\n",
"plt.show()"
2023-07-20 10:19:32 +02:00
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T14:18:20.304971600Z",
"start_time": "2023-07-20T14:18:18.697965900Z"
2023-07-20 10:19:32 +02:00
}
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 84,
2023-07-20 10:19:32 +02:00
"outputs": [
2023-07-20 20:34:19 +02:00
{
"name": "stdout",
"output_type": "stream",
"text": [
"[1. 1. 1. 1. 1. 1.]\n",
"[0. 0. 0. 0. 0. 0.]\n"
]
},
2023-07-20 10:19:32 +02:00
{
"data": {
2023-07-20 20:34:19 +02:00
"text/plain": "<Figure size 640x480 with 1 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGZCAYAAAC0UVoBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAADXjElEQVR4nO39f5SmZXXnje5HCptCCqqgBBrp5m7txh8D2qgoEZUykREzo6Imk8moGTLjUvO+jiczWWdGzpmzIuvNwLhWxvdNzPtOnGQOyXkxEw0Sl47BGE0KjPEHJLSgwYQmfUOp1WkLq7oLKDoUPOePour57O9de/uoLVUF+7tWrdrPc1/3dV/39eu5ru/e1969fr/ft0KhUCgUCoUNwFM2ugCFQqFQKBSevKiFSKFQKBQKhQ1DLUQKhUKhUChsGGohUigUCoVCYcNQC5FCoVAoFAobhlqIFAqFQqFQ2DDUQqRQKBQKhcKGoRYihUKhUCgUNgwjG12A74VHH33Uvv3tb9vY2Jj1er2NLk6hUCgUCoUh0O/3bXFx0c466yx7ylNi3mPTL0S+/e1v244dOza6GIVCoVAoFH4AzMzM2Nlnnx1e3/QLkbGxscekGTM72czuxtVFSf2PIC9DnpZ0r4W8BHlU0j0M+fiBOC7JFoa4xz1Hn3UP5DFJx/tuh3ySpHsF5O/G+Z2AMj3EC1+DfJbkfSrkP4P8chsO2k5sm9Mhf0vSHW/r4/Tge/NF/W6Yylcfq/iRhyXhIcgs93cG4gkv9rc8FLXnqRYDhThD+uHTIP8d20nrhy/1d5CfCVn7OD/zXbUfRm12d5Lu2QNRX/27rPR7IbOsfyM34X2PQ96PyNh6Pt7pdrbn3/l0dgJk9nmtV5bj6ZBZR6w7M//CfL+nS7rRQNb5gp2ZA/dZFuNGyNzMzUs69l8891RpNDee2MfP8elYrTqciEeiAcry7ZRrbIvzIGt9RfO6puPcxP7xKsgyL53wjIH80D5c+Ec+ncv725DPk3QfhjwF+ZOS7l9DjubGL8lnjid9d4Jj+tthKt+4rFedL77+2P8HzOwn8Tu+PnqbPdbMkSNH7JRTTjGzg7ayELlvcHFSVlhsm1n8aD/3+T7dnV/GB7lGTKCiOTZ0+bbc4gMbmwPthXKT/iCs4rPy+TlButPkM5+LRu9Jh+1HnfHO4HszX/ZvDsRdUv8HWRzOQDKQx5qBvIj8LF4xm93FzCE/V9LhPXagbWdk0hvBJLtsCXgfy8pJ+htyD+uL9XqxpAv64aj0jaXog/ahL6xfhqgfm5nZRyBfinvkR2iek+LrNJP1wd++Gb24tL6ctku0yG8lHX4o7A7kLWNwmW2LeSXthxz8LMOXJR2vcZLXcYu2oea5M06j+ULSbUe62ai+FKiHSdT/3Dd9Ms63c9EGyczPqajXntRrfxofODayDVwwz3V2Hey/2btHdcT21I74DFsfWlb2FS5K9Lfgr4LyNJKO8yjfPeuvRyC3SRn47pxr90g6viPfT+toNd0RM9tlhw8ftpNPPjksZRmrFgqFQqFQ2DDUQqRQKBQKhcKGYQupZuZsRTWzH1dV73TGQKRK4kxJtgDZsWl3+XSOlvp7yEoD7h6IE3juPOk9paVUpbBugaRMpDylDKTSHf2uStoPQv53A3ESXysl62g35qeUIMreQ3nSHob32C7qAEctf3ogTkA1MK9tEenc/17Soa+4Om8lHdspoiw174CCHBF6fXlI6nwC8jzapkN1a/9dBcqamUEts/0kYT9SCWX1yvI0ko7jmPew7yktHOn69aVUVfAYeqKGdfUF2ntkyqdbRplGUOfLHN9aVqgkRnAtVQNSjdf4S1TXLbHPZzY/qCPe72zDzKwf5af9ifMr5L6qdTlmsj7FzkfVBdUJjeTNMlHNpUZI7UCcRB5zX5B0/A3huH04+N7Mq0ioptExTJXLbouB8rnn6twWzVkZmI5567iNTBSytmVdvsQn2/FYXTx6xOxbp5RqplAoFAqFwubFpj81s4YTjl9hOdxCUHfkWL2NR4yDSR7ZDgwJL8Bq9LYjko4sCFexlFu5h5+5mtcdTmQseY9PNs/3ZX6yA5gACzKPMszhnUZldbyEdG43oPWPsvczYy+uzKcG4kFJNol6Zfmc4bDshJzhcAP5DJ/Ool2g9hvsHHroHzQeX9K8g13Ncsac3Bd8b2bz3GmhLvu6Y0p29RFc07D9MnawTTJkf2PbSB8YRT07w2b2PWE2xnBtkW2muzaeYIERZIeZ4259aiAvT0u6l+La0vrfd4Dd8DJ37sricjyxHtRgOblmwaWHyaK0uNDE97Mfjsk8QNb0APtetttHGXS6mIEhsTOezE6XsR+yr6khOJiKOZZVxyDHO8cgx5yOb/SB83HtDs371ZDZftpf+ZnPUoaFfY9lPTF4jpmNIr8l3p+8U2qsSqDdJ6WsM6uHLh5I7h+gGJFCoVAoFAobhlqIFAqFQqFQ2DBsIWPVA2Z2shhBZvQx6bjM4IZ0mtKmNKAi/aiqhjbJIwLpWhogChXWjwyMEm9dI6BUH5FrodEhDatEY8f8eMuiGEzRGG3oXpVRluRySdOjLToGmyzTNGRVuTRDlc7XCylj3t/KPZFDudhYy/tySNSKZGT7+ly2G/sX0u2QvGeisgrNzHrukxLX/o50NKA+TpKRsV1Cu29H+WZbuamBnHmrI9in5N3HUL7FLD92ejYA6O1RoaYjW0KNUtEf0mDZtQf60ZioaRYT/z1raJLnRM808/2fahWpVxplj+P7ORnfTj3H+ueYy/wvsXxa4ZG6b1rSBf5L2K11nnOqC757ZlxKtXrmA4r9S/ODKtC58HAW53JPpG5SBIawHQP7oB+GhvhHzGyyjFULhUKhUChsXmwhRuSvLT2ua2Z+ZwpmobML4eqUq041zOHqrwnukXRulc802XGrxPjMvRMMNkfEkIxHDB2Tk+yu3e6AOwPdST0juKblboN7dAfAeuXSXq3ZAjbI7c61XlkP3B1oWdF3JnFtTo87Mw8+iwbLavPN90iOZIb1oDsX9nMaAatRH/tvZHCm78cyZB4aAXdCVMcCx2ib5I2d2iiuLWX11QbXMmaU13RHNhpcU2N0XmMfQP1rOCznSRZ11JM2c6wW21n660hwSaMnOOaQdY7n8ji4mdk8+wDLkPln54P/Sq5dsv4tE9Kv5+GCfgJhN+YTlsh5n+VYaCXdzyAd+1Q21xKYvzrlZvk+Cvktkgf7CuqYnqXNxLt0BrZn5JG6lXv4vpy7lWkiUF8j8vvhCBeUYVTG99Jq2zxgZj9RjEihUCgUCoXNi61zfNdOtZWdieotCe6EsFrrK3OC405cyXUc/US7AykD9doL+J6r95NUl8vVMle3usMPjk8pOeScmMGGImvh46fW/35JV65cIYOJWdadKFka2hG8WtKxLujY6iuSjk5ysAPrfzJIo+DOVndC2JHfx9W8bjG52+N7JLZFbtcWOHcyM7/TZvl0R85nvTT4Xp57AvqNe6zokB27xPw0b7Tt0vRAnpjy6eYjtlF11+iY7mgqxlzqAI7XlG1hp0e5O8wod8rYLfZkt+jqBcwCm31GjslPgmGcQ1/rK4vFto5imJjZMpiOxcyxVTRnYawu6JzSQs52zcyPtgd6fJcMAup4XnbXdj6uZY68ADcns77kKLXj+TNWmLG9UB721058JtjHjIIF6ZDerCOUdVHjErGeOQdqhujnzrFeFP9FkTkiDGLcLGesGNq58+7L8j9HMSKFQqFQKBQ2DLUQKRQKhUKhsGHYQqqZ0cf+SNs1kobXjg9kM8eppsZxgWGaYpaGSIzBgPsXs6qmOkDLEMTrWNZ3IgWG/Mal3O7YJOm5RI1BWn45Me4lQ7ufvPW0Tzc2hQ+orw5lScqdbcEHyZFAlpX0bOdkdjuQT+L3mpDqGJaBdK8Ym9EYejveYVaydnQo2qkTk4PtxHtEhdPHZ6deY18RqtUZ+ybeJPvBWOjE+mEemcFydHQcNLXSwhNo2/lsbMIIchRGkEsZJQ4VYcdjLQ098X70rnu/ePakN08eye8cuWY/4hjW8c32TOZAqrNc/eH
2023-07-20 10:19:32 +02:00
},
"metadata": {},
2023-07-20 20:34:19 +02:00
"output_type": "display_data"
2023-07-20 10:19:32 +02:00
}
],
"source": [
2023-07-20 20:34:19 +02:00
"print(p0[0,0])\n",
"print(popt[0,0])\n",
"\n",
"data = cropOD[0,0]\n",
"plt.pcolormesh(data,cmap='jet')\n",
"plt.show()"
2023-07-20 10:19:32 +02:00
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T14:24:36.509831700Z",
"start_time": "2023-07-20T14:24:36.334599Z"
2023-07-20 10:19:32 +02:00
}
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 90,
2023-07-20 10:19:32 +02:00
"outputs": [],
"source": [
"\n",
"def polylog(power, numerator):\n",
"\n",
2023-07-20 20:34:19 +02:00
" order = 10\n",
2023-07-20 10:19:32 +02:00
"\n",
" dataShape = numerator.shape\n",
" numerator = np.tile(numerator, (order, 1))\n",
" numerator = np.power(numerator.T, np.arange(1, order+1)).T\n",
"\n",
" denominator = np.arange(1, order+1)\n",
" denominator = np.tile(denominator, (dataShape[0], 1))\n",
" denominator = denominator.T\n",
"\n",
" data = numerator/ np.power(denominator, power)\n",
"\n",
" return np.sum(data, axis=0)\n",
"\n",
"def ThomasFermi_2d(x, y=0.0, centerx=0.0, centery=0.0, amplitude=1.0, sigmax=1.0, sigmay=1.0):\n",
" res = (1- ((x-centerx)/(sigmax))**2 - ((y-centery)/(sigmay))**2)**(3 / 2)\n",
" return amplitude * np.where(res > 0, res, 0)\n",
" # return amplitude * 5 / 2 / np.pi / max(tiny, sigmax * sigmay) * np.where(res > 0, res, 0)\n",
"\n",
"def polylog2_2d(x, y=0.0, centerx=0.0, centery=0.0, amplitude=1.0, sigmax=1.0, sigmay=1.0):\n",
" ## Approximation of the polylog function with 2D gaussian as argument. -> discribes the thermal part of the cloud\n",
" return amplitude * polylog(2, np.exp( -((x-centerx)**2/(2 * (sigmax)**2))-((y-centery)**2/( 2 * (sigmay)**2)) ))\n",
" # 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)) ))\n",
"\n",
"\n",
2023-07-20 20:34:19 +02:00
"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,\n",
" sigmax_bec=1.0, sigmay_bec=1.0, sigmax_th=1.0, sigmay_th=1.0):\n",
" return ThomasFermi_2d(x=x, y=y, centerx=x0_bec, centery=y0_bec,\n",
" amplitude=amp_bec, sigmax=sigmax_bec, sigmay=sigmay_bec\n",
" ) + polylog2_2d(x=x, y=y, centerx=x0_th, centery=y0_th,\n",
" amplitude=amp_th, sigmax=sigmax_th, sigmay=sigmay_th)\n",
2023-07-20 10:19:32 +02:00
"\n",
2023-07-20 20:34:19 +02:00
"\n"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-07-20T14:33:17.937309200Z",
"start_time": "2023-07-20T14:33:17.885650600Z"
}
}
},
{
"cell_type": "code",
"execution_count": 37,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[75.41570333 77.33388115 2.36427288 0.29873033 9.02534473 21.78913789]\n",
"[75.41570333 77.33388115 2.36427288 0.29873033 9.02534473 21.78913789]\n",
"[75.69178082 75.69178082 1.28576093 0.5510404 8. 24. ]\n",
"2.66300321392503\n"
]
}
],
"source": [
"popt_0 = popt[0,0]\n",
"print(popt_0)\n",
"S = np.max(blurred[0,0])/(popt_0[2] + popt_0[3])\n",
"popt_0[2] *= S\n",
"popt_0[3] *= S\n",
"print(popt_0)\n",
"print(p0[0,0])\n",
"print(np.max(blurred[0,0]))"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-07-20T12:16:22.598036600Z",
"start_time": "2023-07-20T12:16:22.566024700Z"
}
}
},
{
"cell_type": "code",
"execution_count": 91,
"outputs": [],
"source": [
"fitModel = lmfit.Model(density_profile_BEC_2d, independent_vars=['x','y'])\n",
"#fitModel.set_param_hint('deltax', value=5)\n",
2023-07-20 10:19:32 +02:00
"\n",
2023-07-20 20:34:19 +02:00
"bval = result_x[0][0].best_values\n",
"S = np.max(blurred[0,0])/(bval['amp_bec'] + bval['amp_th'])\n",
2023-07-20 10:19:32 +02:00
"\n",
2023-07-20 20:34:19 +02:00
"params = lmfit.Parameters()\n",
"#print(bval['sigma_th'])\n",
"params.add_many(\n",
" ('amp_bec', S * bval['amp_bec'], True, 0, 1.3 * np.max(data)),\n",
" ('amp_th',S * bval['amp_th'], True, 0, 1.3 * np.max(data)),\n",
" ('x0_bec',center[0,0,0], True, 0, 150),\n",
" ('y0_bec',center[0,0,1], True, 0, 200),\n",
" ('x0_th',center[0,0,0], True, 0, 150),\n",
" ('y0_th',center[0,0,1], True, 0, 200),\n",
" ('sigmax_bec',bval['sigma_bec'], True, 0, 50),\n",
" ('sigmay_bec',BEC_width_guess[0,0,1], True, 0, 50),\n",
" ('sigmax_th',bval['sigma_th'], True, 0, 50),\n",
" ('sigmay_th',bval['sigma_th'], True, 0, 50)\n",
")\n",
2023-07-20 10:19:32 +02:00
"\n",
2023-07-20 20:34:19 +02:00
"X,Y = np.meshgrid(x, y)\n",
"X_1d = X.flatten()\n",
"Y_1d = Y.flatten()\n",
2023-07-20 10:19:32 +02:00
"\n",
2023-07-20 20:34:19 +02:00
"data1d = data.flatten()\n",
"\n",
"result = fitModel.fit(data1d, x=X_1d, y=Y_1d, params=params)\n",
2023-07-20 10:19:32 +02:00
"\n",
2023-07-20 20:34:19 +02:00
"\n"
2023-07-20 10:19:32 +02:00
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T14:33:22.499517400Z",
"start_time": "2023-07-20T14:33:20.802857800Z"
}
}
},
{
"cell_type": "code",
"execution_count": 92,
"outputs": [
{
"data": {
"text/plain": "<lmfit.model.ModelResult at 0x20d87a58e50>",
"text/html": "<h2> Model</h2> Model(density_profile_BEC_2d) <h2>Fit Statistics</h2><table><tr><td>fitting method</td><td>leastsq</td><td></td></tr><tr><td># function evals</td><td>100</td><td></td></tr><tr><td># data points</td><td>22500</td><td></td></tr><tr><td># variables</td><td>10</td><td></td></tr><tr><td>chi-square</td><td> 142.279826</td><td></td></tr><tr><td>reduced chi-square</td><td> 0.00632636</td><td></td></tr><tr><td>Akaike info crit.</td><td>-113908.185</td><td></td></tr><tr><td>Bayesian info crit.</td><td>-113827.972</td><td></td></tr><tr><td>R-squared</td><td> 0.84849949</td><td></td></tr></table><h2>Variables</h2><table><tr><th> name </th><th> value </th><th> standard error </th><th> relative error </th><th> initial value </th><th> min </th><th> max </th><th> vary </th></tr><tr><td> amp_bec </td><td> 1.36015474 </td><td> 0.00982217 </td><td> (0.72%) </td><td> 2.3177242655279193 </td><td> 0.00000000 </td><td> 3.66861092 </td><td> True </td></tr><tr><td> amp_th </td><td> 0.39877493 </td><td> 0.00413246 </td><td> (1.04%) </td><td> 0.33991000464208837 </td><td> 0.00000000 </td><td> 3.66861092 </td><td> True </td></tr><tr><td> x0_bec </td><td> 75.8839398 </td><td> 0.02102227 </td><td> (0.03%) </td><td> 75.6689655172414 </td><td> 0.00000000 </td><td> 150.000000 </td><td> True </td></tr><tr><td> y0_bec </td><td> 64.2566209 </td><td> 0.08568536 </td><td> (0.13%) </td><td> 63.53103448275861 </td><td> 0.00000000 </td><td> 200.000000 </td><td> True </td></tr><tr><td> x0_th </td><td> 77.5639125 </td><td> 0.11996776 </td><td> (0.15%) </td><td> 75.6689655172414 </td><td> 0.00000000 </td><td> 150.000000 </td><td> True </td></tr><tr><td> y0_th </td><td> 63.2611524 </td><td> 0.12707327 </td><td> (0.20%) </td><td> 63.53103448275861 </td><td> 0.00000000 </td><td> 200.000000 </td><td> True </td></tr><tr><td> sigmax_bec </td><td> 7.29956469 </td><td> 0.04255796 </td><td> (0.58%) </td><td> 9.046333505508725 </td><td> 0.00000000 </td><td> 50.0000000 </td><td> True </td></tr><tr><td> sigmay_bec </td><td> 25.6881035 </td><td> 0.14446069 </td><td> (0.56%) </td><td> 25.0 </td><td> 0.00000000 </td><td> 50.0000000 </td><td> True </td></tr><tr><td> sigmax_th </td><td> 19.7527225 </td><td> 0.16378744 </td><td> (0.83%) </td><td> 20.287472864718893 </td><td> 0.00000000 </td><td> 50.0000000 </td><td> True </td></tr><tr><td> sigmay_th </td><td> 18.2433266 </td><td> 0.14078251 </td><td> (0.77%) </td><td> 20.287472864718893 </td><td> 0.00000000 </td><td> 50.0000000 </td><td> True </td></tr></table><h2>Correlations (unreported correlations are < 0.100)</h2><table><tr><td>amp_th</td><td>sigmax_th</td><td>-0.6954</td></tr><tr><td>amp_bec</td><td>amp_th</td><td>-0.5655</td></tr><tr><td>amp_th</td><td>sigmay_th</td><td>-0.5404</td></tr><tr><td>y0_bec</td><td>y0_th</td><td>-0.5396</td></tr><tr><td>amp_th</td><td>sigmax_bec</td><td>-0.4387</td></tr><tr><td>amp_bec</td><td>sigmax_th</td><td>+0.3938</td></tr><tr><td>sigmay_bec</td><td>sigmay_th</td><td>-0.3383</td></tr><tr><td>amp_bec</td><td>sigmay_bec</td><td>-0.3311</td></tr><tr><td>amp_bec</td><td>sigmay_th</td><td>+0.3007</td></tr><tr><td>sigmax_bec</td><td>sigmax_th</td><td>+0.2814</td></tr><tr><td>x0_bec</td><td>x0_th</td><td>-0.2249</td></tr><tr><td>sigmax_bec</td><td>sigmay_th</td><td>+0.2018</td></tr><tr><td>sigmax_th</td><td>sigmay_th</td><td>+0.1620</td></tr><tr><td>sigmax_bec</td><td>sigmay_bec</td><td>-0.1376</td></tr><tr><td>amp_bec</td><td>x0_th</td><td>+0.1086</td></tr></table>"
},
"execution_count": 92,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"result"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-07-20T14:33:25.858457Z",
"start_time": "2023-07-20T14:33:25.823674300Z"
}
}
},
{
"cell_type": "markdown",
"source": [
"## Fitting class"
],
"metadata": {
"collapsed": false
}
},
{
"cell_type": "code",
"execution_count": 93,
"outputs": [
{
"ename": "type",
"evalue": "operands could not be broadcast together with shapes (150,1500) (10,) ",
"output_type": "error",
"traceback": [
"\u001B[1;31m---------------------------------------------------------------------------\u001B[0m",
"\u001B[1;31mValueError\u001B[0m Traceback (most recent call last)",
"Cell \u001B[1;32mIn[93], line 1\u001B[0m\n\u001B[1;32m----> 1\u001B[0m fit \u001B[38;5;241m=\u001B[39m fitModel\u001B[38;5;241m.\u001B[39mfunc(X,Y, \u001B[38;5;241m*\u001B[39m\u001B[38;5;241m*\u001B[39mresult\u001B[38;5;241m.\u001B[39mbest_values)\n\u001B[0;32m 3\u001B[0m fig, axs \u001B[38;5;241m=\u001B[39m plt\u001B[38;5;241m.\u001B[39msubplots(\u001B[38;5;241m2\u001B[39m, \u001B[38;5;241m2\u001B[39m, figsize\u001B[38;5;241m=\u001B[39m(\u001B[38;5;241m12\u001B[39m, \u001B[38;5;241m10\u001B[39m))\n\u001B[0;32m 5\u001B[0m ax \u001B[38;5;241m=\u001B[39m axs[\u001B[38;5;241m0\u001B[39m, \u001B[38;5;241m0\u001B[39m]\n",
"Cell \u001B[1;32mIn[90], line 32\u001B[0m, in \u001B[0;36mdensity_profile_BEC_2d\u001B[1;34m(x, y, amp_bec, amp_th, x0_bec, y0_bec, x0_th, y0_th, sigmax_bec, sigmay_bec, sigmax_th, sigmay_th)\u001B[0m\n\u001B[0;32m 28\u001B[0m \u001B[38;5;28;01mdef\u001B[39;00m \u001B[38;5;21mdensity_profile_BEC_2d\u001B[39m(x, y\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m0.0\u001B[39m, amp_bec\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m1.0\u001B[39m, amp_th\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m1.0\u001B[39m, x0_bec\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m0.0\u001B[39m, y0_bec\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m0.0\u001B[39m, x0_th\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m0.0\u001B[39m, y0_th\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m0.0\u001B[39m,\n\u001B[0;32m 29\u001B[0m sigmax_bec\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m1.0\u001B[39m, sigmay_bec\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m1.0\u001B[39m, sigmax_th\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m1.0\u001B[39m, sigmay_th\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m1.0\u001B[39m):\n\u001B[0;32m 30\u001B[0m \u001B[38;5;28;01mreturn\u001B[39;00m ThomasFermi_2d(x\u001B[38;5;241m=\u001B[39mx, y\u001B[38;5;241m=\u001B[39my, centerx\u001B[38;5;241m=\u001B[39mx0_bec, centery\u001B[38;5;241m=\u001B[39my0_bec,\n\u001B[0;32m 31\u001B[0m amplitude\u001B[38;5;241m=\u001B[39mamp_bec, sigmax\u001B[38;5;241m=\u001B[39msigmax_bec, sigmay\u001B[38;5;241m=\u001B[39msigmay_bec\n\u001B[1;32m---> 32\u001B[0m ) \u001B[38;5;241m+\u001B[39m \u001B[43mpolylog2_2d\u001B[49m\u001B[43m(\u001B[49m\u001B[43mx\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43mx\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43my\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43my\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mcenterx\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43mx0_th\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mcentery\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43my0_th\u001B[49m\u001B[43m,\u001B[49m\n\u001B[0;32m 33\u001B[0m \u001B[43m \u001B[49m\u001B[43mamplitude\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43mamp_th\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43msigmax\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43msigmax_th\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43msigmay\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43msigmay_th\u001B[49m\u001B[43m)\u001B[49m\n",
"Cell \u001B[1;32mIn[90], line 24\u001B[0m, in \u001B[0;36mpolylog2_2d\u001B[1;34m(x, y, centerx, centery, amplitude, sigmax, sigmay)\u001B[0m\n\u001B[0;32m 22\u001B[0m \u001B[38;5;28;01mdef\u001B[39;00m \u001B[38;5;21mpolylog2_2d\u001B[39m(x, y\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m0.0\u001B[39m, centerx\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m0.0\u001B[39m, centery\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m0.0\u001B[39m, amplitude\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m1.0\u001B[39m, sigmax\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m1.0\u001B[39m, sigmay\u001B[38;5;241m=\u001B[39m\u001B[38;5;241m1.0\u001B[39m):\n\u001B[0;32m 23\u001B[0m \u001B[38;5;66;03m## Approximation of the polylog function with 2D gaussian as argument. -> discribes the thermal part of the cloud\u001B[39;00m\n\u001B[1;32m---> 24\u001B[0m \u001B[38;5;28;01mreturn\u001B[39;00m amplitude \u001B[38;5;241m*\u001B[39m \u001B[43mpolylog\u001B[49m\u001B[43m(\u001B[49m\u001B[38;5;241;43m2\u001B[39;49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mnp\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mexp\u001B[49m\u001B[43m(\u001B[49m\u001B[43m \u001B[49m\u001B[38;5;241;43m-\u001B[39;49m\u001B[43m(\u001B[49m\u001B[43m(\u001B[49m\u001B[43mx\u001B[49m\u001B[38;5;241;43m-\u001B[39;49m\u001B[43mcenterx\u001B[49m\u001B[43m)\u001B[49m\u001B[38;5;241;43m*\u001B[39;49m\u001B[38;5;241;43m*\u001B[39;49m\u001B[38;5;241;43m2\u001B[39;49m\u001B[38;5;241;43m/\u001B[39;49m\u001B[43m(\u001B[49m\u001B[38;5;241;43m2\u001B[39;49m\u001B[43m \u001B[49m\u001B[38;5;241;43m*\u001B[39;49m\u001B[43m \u001B[49m\u001B[43m(\u001B[49m\u001B[43msigmax\u001B[49m\u001B[43m)\u001B[49m\u001B[38;5;241;43m*\u001B[39;49m\u001B[38;5;241;43m*\u001B[39;49m\u001B[38;5;241;43m2\u001B[39;49m\u001B[43m)\u001B[49m\u001B[43m)\u001B[49m\u001B[38;5;241;43m-\u001B[39;49m\u001B[43m(\u001B[49m\u001B[43m(\u001B[49m\u001B[43my\u001B[49m\u001B[38;5;241;43m-\u001B[39;49m\u001B[43mcentery\u001B[49m\u001B[43m)\u001B[49m\u001B[38;5;241;43m*\u001B[39;49m\u001B[38;5;241;43m*\u001B[39;49m\u001B[38;5;241;43m2\u001B[39;49m\u001B[38;5;241;43m/\u001B[39;49m\u001B[43m(\u001B[49m\u001B[43m \u001B[49m\u001B[38;5;241;43m2\u001B[39;49m\u001B[43m \u001B[49m\u001B[38;5;241;43m*\u001B[39;49m\u001B[43m \u001B[49m\u001B[43m(\u001B[49m\u001B[43msigmay\u001B[49m\u001B[43m)\u001B[49m\u001B[38;5;241;43m*\u001B[39;49m\u001B[38;5;241;43m*\u001B[39;49m\u001B[38;5;241;43m2\u001B[39;49m\u001B[43m)\u001B[49m\u001B[43m)\u001B[49m\u001B[43m \u001B[49m\u001B[43m)\u001B[49m\u001B[43m)\u001B[49m\n",
"Cell \u001B[1;32mIn[90], line 7\u001B[0m, in \u001B[0;36mpolylog\u001B[1;34m(power, numerator)\u001B[0m\n\u001B[0;32m 5\u001B[0m dataShape \u001B[38;5;241m=\u001B[39m numerator\u001B[38;5;241m.\u001B[39mshape\n\u001B[0;32m 6\u001B[0m numerator \u001B[38;5;241m=\u001B[39m np\u001B[38;5;241m.\u001B[39mtile(numerator, (order, \u001B[38;5;241m1\u001B[39m))\n\u001B[1;32m----> 7\u001B[0m numerator \u001B[38;5;241m=\u001B[39m \u001B[43mnp\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mpower\u001B[49m\u001B[43m(\u001B[49m\u001B[43mnumerator\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mT\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mnp\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43marange\u001B[49m\u001B[43m(\u001B[49m\u001B[38;5;241;43m1\u001B[39;49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43morder\u001B[49m\u001B[38;5;241;43m+\u001B[39;49m\u001B[38;5;241;43m1\u001B[39;49m\u001B[43m)\u001B[49m\u001B[43m)\u001B[49m\u001B[38;5;241m.\u001B[39mT\n\u001B[0;32m 9\u001B[0m denominator \u001B[38;5;241m=\u001B[39m np\u001B[38;5;241m.\u001B[39marange(\u001B[38;5;241m1\u001B[39m, order\u001B[38;5;241m+\u001B[39m\u001B[38;5;241m1\u001B[39m)\n\u001B[0;32m 10\u001B[0m denominator \u001B[38;5;241m=\u001B[39m np\u001B[38;5;241m.\u001B[39mtile(denominator, (dataShape[\u001B[38;5;241m0\u001B[39m], \u001B[38;5;241m1\u001B[39m))\n",
"\u001B[1;31mValueError\u001B[0m: operands could not be broadcast together with shapes (150,1500) (10,) "
]
}
],
"source": [
"fit = fitModel.func(X,Y, **result.best_values)\n",
"\n",
"fig, axs = plt.subplots(2, 2, figsize=(12, 10))\n",
"\n",
"ax = axs[0, 0]\n",
"art = ax.pcolormesh(X, Y, data, vmin=0, vmax=vmax, shading='auto')\n",
"plt.colorbar(art, ax=ax, label='z')\n",
"ax.set_title('Data')\n",
"\n",
"# Plot gaussian 2d Fit + legend including Width parameters\n",
"ax = axs[0, 1]\n",
"\n",
"art = ax.pcolormesh(X, Y, fit, vmin=0, vmax=vmax, shading='auto')\n",
"plt.colorbar(art, ax=ax, label='z')"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-07-20T14:33:36.731593500Z",
"start_time": "2023-07-20T14:33:36.607896800Z"
}
}
},
{
"cell_type": "code",
"execution_count": 18,
"outputs": [],
"source": [],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-07-20T08:27:55.535706600Z",
"start_time": "2023-07-20T08:27:55.523609400Z"
2023-07-20 10:19:32 +02:00
}
}
},
{
"cell_type": "code",
"execution_count": null,
"outputs": [],
"source": [],
"metadata": {
"collapsed": false
}
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
