641 lines
663 KiB
Plaintext
641 lines
663 KiB
Plaintext
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"source": [
|
||
|
|
"# Import supporting package"
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": false
|
||
|
|
}
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 2,
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"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": {
|
||
|
|
"collapsed": false
|
||
|
|
}
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 20,
|
||
|
|
"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",
|
||
|
|
"def thermal(x, x0, amp, sigma, order = 15):\n",
|
||
|
|
" 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",
|
||
|
|
"def density_1d(x, x0_bec, x0_th, amp_bec, amp_th, sigma_bec, sigma_th, polyorder=15):\n",
|
||
|
|
" return thermal(x, x0_th, amp_th, sigma_th, polyorder) + Thomas_Fermi_1d(x, x0_bec, amp_bec, sigma_bec)"
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": false,
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2023-07-20T07:56:43.191095100Z",
|
||
|
|
"start_time": "2023-07-20T07:56:43.145698800Z"
|
||
|
|
}
|
||
|
|
}
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 21,
|
||
|
|
"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": {
|
||
|
|
"end_time": "2023-07-20T07:56:51.845589300Z",
|
||
|
|
"start_time": "2023-07-20T07:56:45.635528300Z"
|
||
|
|
}
|
||
|
|
}
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 22,
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": "<Figure size 1000x900 with 10 Axes>",
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA4gAAANiCAYAAAAjb+cqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOx9fXhdVZ3umzahTSEpTSGBpm0iBGgHihC1dSxOAfGCHygIgqIj6ADlXhgd7owfiCN29BG9ogyjXuU6qIzjBzwDOCAjo0VaxsJQICBFWmmAFBogQVKaQEPJKfv+sc971rt+e62dpE3SVPb7PHlyztlrr732+vit3/eqSpIkQYECBQoUKFCgQIECBQoUeM1jyu5uQIECBQoUKFCgQIECBQoUmBwoBMQCBQoUKFCgQIECBQoUKACgEBALFChQoECBAgUKFChQoEAZhYBYoECBAgUKFChQoECBAgUAFAJigQIFChQoUKBAgQIFChQooxAQCxQoUKBAgQIFChQoUKAAgEJALFCgQIECBQoUKFCgQIECZRQCYoECBQoUKFCgQIECBQoUAFAIiAUKFChQoECBAgUKFChQoIw9VkD84Q9/iKqqKrS2tv5JP5Po6upCVVUVqqqq0NXVNeH37yz4zFWrVk3YMwuMD4o1N7H37yyKNVdgpCjW9MTev7Mo1nSBAgUmGnusgJiHn//85/jCF76An//857u7KX+y+Md//Ed84QtfwIMPPri7mzLpccEFF1Q2+N3BFI0FHnvsMfzwhz8EAGzatAn7778/TjzxRNxwww0Adn3NdXR04MMf/jDmzp2LadOm4cADD8Spp56KRx55ZET333HHHTj11FNx4IEHYtq0aZg7dy4+/OEPo6OjY9jnfuc738F5552H9vZ2TJs2bdKOU7HmRo4/lTW3fPlyvO51r8P06dMza25X0dHRgauvvhpAuqa55n7zm98ACK/poaEhfO9738Py5cuxZMkSzJs3D7W1tZgxYwYOPvhgnHXWWVi5cuWInr9t2zaceOKJaGxsxPTp0/G6170Oy5cvR2dn56jf5a677sLUqVOHFaS2b99e+XzyySdj3333RU1NDfbff38cd9xx+OY3v4lt27aN+vk7i2JNjxzFms7H9u3b8c1vfhPHHHMMZs2ahenTp6O1tRXnnntu7j6aJAnuvvtufO5zn8Oxxx6LpqYm1NTUYObMmXjDG96ASy65BN3d3bnPbm1trYxN7O+YY46JtvuWW27BRRddhDe+8Y27fU0WECR7KH7wgx8kAJKWlpbMtbPPPjsBkJx99tkT9szxxhNPPJEASAAkTzzxxITfb9HS0pIASH7wgx/kljvssMOSww47LLnnnnt2+Zl7In7zm98kVVVVlb7fHXNnV3HrrbcmM2bMqLxDVVVVMmXKlMr3j370o7u05r73ve8l1dXVlfpmzpzp9dlw/XbZZZd5bZs5c2ble3V1dfK9730vei/nsf1raWkp1tweij/FNVdfX59Zc6+++upO12/XXFVVlddnl112WXBNP/fcc946qaqqSmbNmpVMnTrV+/3ss89OhoaGgs/+/ve/75WdMmVKUl9fX/k+Y8aM5NZbbx3xuwwODiYLFizw6rzjjjuCZZcuXeqVq66u9ugFgKStrS3ZuHFj8P5iTe8eFGs6H88880xy9NFHV+qqqalJZs2aVfk+bdq05Mc//nHw3i996UuZNb3vvvt6/V1fX5/8+7//e/T5nMf19fVJU1NT8O8973lP8N4TTjhhl9ZkgfHDn6QFscDkwYYNG7BhwwYsXrx4dzdlwrFt2zacd955qK6uxhvf+Mbd3ZydwhNPPIEzzjgD27ZtQ1tbGwBg/vz52Lp1Kz7/+c8DAH7wgx/g4Ycf3qn67777blxwwQUolUo45ZRT8NRTT+GFF17Ac889h+XLl1fKxTSI119/PVasWAEAWL58OZ577jm88MILeOqpp3DKKaegVCrhggsuwN133x28f6+99sJRRx2Fj33sY/jWt76Fv/zLv9yp95hMKNbcn86aW7p0Kf7whz9g69atmTX3ta99bafq1zV39NFHA0jXtK65FStWBF0op02bhr/+67/Gddddh66uLmzfvh19fX145ZVXsG7dOnzgAx8AAFx77bW44oorgs//5S9/Wfl82WWXVd5tw4YNeMtb3oJt27bhjDPOwBNPPDGi91mxYkXl3uFQKpUqn3/xi19g+/bteOGFF9DX14evfvWrmDFjBjo7O/GOd7wDL7/88oiePxEo1nSxpmNIkgSnnXYaHnjgAdTW1uJ73/se+vv70dfXh6effhof+chHsH37dpxzzjm4//77M/cPDQ2hvr4eF1xwAX7zm9/gpZdewpYtW/DSSy/hhhtuwPz589Hf34/3v//9WL9+fW5brrrqKjz77LPBv3//938P3jM0NISWlhZ84QtfQEdHxx6zJl8T2N0S6s6isCBO7P0WI9V8vpbxN3/zNwmA5NJLL63MyT1N8/nhD384AZAccMABybe//e3MO5x//vkVjeXOrLljjjkmAZAsWrQoeeWVVzLXjzjiiARAMnXq1KRUKnnXSqVSZR6edNJJmXu3b99euf+YY44JPt/WSWvknmxBfC3jT23NbdmyJXOda66+vj7p6+sbdf265v75n/8500cnnnhiAiDZe++9R72mX3311eQtb3lLRetv0dfXl0yfPj0BkOyzzz7B6wcccEACIPnwhz887PPuv//+pLq6Omlra0t++ctfDmtBvP7663PX5E9/+tPK9R/96EeZ68WanngUazoft9xyS2VO/uM//mOwzJvf/OYEQHL88cdnrj3wwAO5z3z88ceT2traBEDyV3/1V8EyuzKP/+u//ivqbZAkw6/JAuOHcRcQ77jjjuT0009P5syZk+y1117J7Nmzk+OPPz75/ve/n2HOLO6+++7kve99bzJ79uxk+vTpyaGHHpp89rOfTQYGBoLC2h133BF0F9O/2MYxEthn3nfffcn73//+5IADDkj22muv5HWve11y8cUXD7vAOzs7kwsuuCBpa2tLpk+fntTV1SVHH310smLFimTr1q3Be0ayMW3evDk5//zzk7lz5yZ77bVX0tzcnJxzzjnJxo0bx2xjU3e+2J8i1u+2PV1dXcm5556bzJs3L5k2bVpy0EEHJZdeemny4osvVu5Zt25d8qEPfSiZO3duMm3atKStrS354he/GBQs7LM+8YlPJH/2Z3+W7L333kltbW1y2GGHJR//+MeTTZs27XRf5OHuu+9OpkyZkhx66KHJ4ODghG5sY7nm6Gby2c9+NrjmlHiPds099thjlTLXXnttsMynP/3pSpnf/OY33pqjUAog+cUvfhG8/4c//GGlzFlnnTXsmhutgFisuTBea2turPDiiy9WmLEVK1ZUfrdrmuP40Y9+dFRretq0aZV7r7766uCavvLKK3dpH7344osTAMn06dMz16655ppKHc3NzUmSZPfRfffdt3K/zkWL9evXJw0NDQmAZK+99vLc92L0YLg1uWPHjopgvM8++xRr2jyrWNOjR2xNK3Qcv//974+q/gsuuCABUoVObAxvuOGGSv07M1bvfOc7EwDJkUceGbw+noqOHTt2JHV1dQmA5OMf//iY118gjnEVELlRAM6vWWMVjj/++KS/vz947zXXXOP5Z8+cObOyMS5YsCD5xje+kSEUa9asSZqamioayunTp2f8oNesWbPT76Ob6Y9//OMKgzpz5kyvrYcffngyMDAQrOO6667zNum6ujrv+7x585JHHnkkc99wG9P999/v+ZzX1tYm++yzTwKkWqnrrrtuTDa2r33ta0lTU1PlfUM+54qRbGw33HBDhSmor6/35shb3/rW5JVXXkl+8YtfVBgAG6N25plnRtv7r//6r17/Tps2rUKs2f//+Z//udP9EcLLL7+cLFy4MKmqqkpWrVqVJEkyYRvbWK45MkoAktbW1uia0/pHs+a++93vVu7r6emJtollTj75ZG/N6RyIrbmenp4M45W35kYjIBZrLozX2pobS9x2222Vflq7dm2SJOE1rfN5NGtax+GQQw4Jruk777zTG+v
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"imageAnalyser.center = (960, 890)\n",
|
||
|
|
"imageAnalyser.span = (150, 200)\n",
|
||
|
|
"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": {
|
||
|
|
"end_time": "2023-07-20T07:56:57.392382600Z",
|
||
|
|
"start_time": "2023-07-20T07:56:53.716507300Z"
|
||
|
|
}
|
||
|
|
}
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 14,
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"fitting time: 118.55363845825195 ms\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"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",
|
||
|
|
"y = np.linspace(0,199, 200)\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",
|
||
|
|
"\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": {
|
||
|
|
"end_time": "2023-07-19T14:22:38.020578900Z",
|
||
|
|
"start_time": "2023-07-19T14:22:37.850172600Z"
|
||
|
|
}
|
||
|
|
}
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 17,
|
||
|
|
"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": {
|
||
|
|
"end_time": "2023-07-19T14:29:07.056739700Z",
|
||
|
|
"start_time": "2023-07-19T14:29:05.074878Z"
|
||
|
|
}
|
||
|
|
}
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 16,
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"[75.69178082 75.69178082 1.28576093 0.5510404 8. 24. ]\n",
|
||
|
|
"[75.41570333 77.33388115 2.56706786 0.32435386 9.02534473 21.78913789]\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"print(p0[0,0])\n",
|
||
|
|
"print(popt[0,0])"
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": false,
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2023-07-19T14:23:03.879599500Z",
|
||
|
|
"start_time": "2023-07-19T14:23:03.844141600Z"
|
||
|
|
}
|
||
|
|
}
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 18,
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"ename": "type",
|
||
|
|
"evalue": "name 'data' is not defined",
|
||
|
|
"output_type": "error",
|
||
|
|
"traceback": [
|
||
|
|
"\u001B[1;31m---------------------------------------------------------------------------\u001B[0m",
|
||
|
|
"\u001B[1;31mNameError\u001B[0m Traceback (most recent call last)",
|
||
|
|
"Cell \u001B[1;32mIn[18], line 6\u001B[0m\n\u001B[0;32m 2\u001B[0m \u001B[38;5;66;03m#fitModel.set_param_hint('deltax', value=5)\u001B[39;00m\n\u001B[0;32m 4\u001B[0m params \u001B[38;5;241m=\u001B[39m fitModel\u001B[38;5;241m.\u001B[39mmake_params()\n\u001B[0;32m 5\u001B[0m params\u001B[38;5;241m.\u001B[39madd_many(\n\u001B[1;32m----> 6\u001B[0m (\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mBEC_amplitude\u001B[39m\u001B[38;5;124m'\u001B[39m, np\u001B[38;5;241m.\u001B[39mmax(\u001B[43mdata\u001B[49m)\u001B[38;5;241m/\u001B[39m\u001B[38;5;241m2\u001B[39m, \u001B[38;5;28;01mTrue\u001B[39;00m, \u001B[38;5;241m0\u001B[39m, np\u001B[38;5;241m.\u001B[39mmax(data)),\n\u001B[0;32m 7\u001B[0m (\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mthermal_amplitude\u001B[39m\u001B[38;5;124m'\u001B[39m,np\u001B[38;5;241m.\u001B[39mmax(data)\u001B[38;5;241m/\u001B[39m\u001B[38;5;241m2\u001B[39m, \u001B[38;5;28;01mTrue\u001B[39;00m, \u001B[38;5;241m0\u001B[39m, np\u001B[38;5;241m.\u001B[39mmax(data)),\n\u001B[0;32m 8\u001B[0m (\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mBEC_centerx\u001B[39m\u001B[38;5;124m'\u001B[39m,\u001B[38;5;241m1\u001B[39m),\n\u001B[0;32m 9\u001B[0m (\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mBEC_centery\u001B[39m\u001B[38;5;124m'\u001B[39m,\u001B[38;5;241m1\u001B[39m),\n\u001B[0;32m 10\u001B[0m (\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mthermal_centerx\u001B[39m\u001B[38;5;124m'\u001B[39m,\u001B[38;5;241m1\u001B[39m),\n\u001B[0;32m 11\u001B[0m (\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mthermal_centery\u001B[39m\u001B[38;5;124m'\u001B[39m,\u001B[38;5;241m1\u001B[39m),\n\u001B[0;32m 12\u001B[0m (\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mBEC_sigmay\u001B[39m\u001B[38;5;124m'\u001B[39m,\u001B[38;5;241m1\u001B[39m),\n\u001B[0;32m 13\u001B[0m (\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mthermal_sigmax\u001B[39m\u001B[38;5;124m'\u001B[39m,\u001B[38;5;241m1\u001B[39m),\n\u001B[0;32m 14\u001B[0m (\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mdeltax\u001B[39m\u001B[38;5;124m'\u001B[39m,\u001B[38;5;241m1\u001B[39m),\n\u001B[0;32m 15\u001B[0m (\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mthermalAspectRatio\u001B[39m\u001B[38;5;124m'\u001B[39m,\u001B[38;5;241m1\u001B[39m)\n\u001B[0;32m 16\u001B[0m )\n\u001B[0;32m 17\u001B[0m params\n",
|
||
|
|
"\u001B[1;31mNameError\u001B[0m: name 'data' is not defined"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"fitModel = DensityProfileBEC2dModel()\n",
|
||
|
|
"#fitModel.set_param_hint('deltax', value=5)\n",
|
||
|
|
"\n",
|
||
|
|
"params = fitModel.make_params()\n",
|
||
|
|
"params.add_many(\n",
|
||
|
|
" ('BEC_amplitude', np.max(data)/2, True, 0, np.max(data)),\n",
|
||
|
|
" ('thermal_amplitude',np.max(data)/2, True, 0, np.max(data)),\n",
|
||
|
|
" ('BEC_centerx',1),\n",
|
||
|
|
" ('BEC_centery',1),\n",
|
||
|
|
" ('thermal_centerx',1),\n",
|
||
|
|
" ('thermal_centery',1),\n",
|
||
|
|
" ('BEC_sigmay',1),\n",
|
||
|
|
" ('thermal_sigmax',1),\n",
|
||
|
|
" ('deltax',1),\n",
|
||
|
|
" ('thermalAspectRatio',1)\n",
|
||
|
|
")\n",
|
||
|
|
"params"
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": false,
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2023-07-19T14:37:09.050199900Z",
|
||
|
|
"start_time": "2023-07-19T14:37:08.977867500Z"
|
||
|
|
}
|
||
|
|
}
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 12,
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": "Parameters([('BEC_amplitude', <Parameter 'BEC_amplitude', value=1.0, bounds=[0:inf]>), ('thermal_amplitude', <Parameter 'thermal_amplitude', value=1.0, bounds=[0:inf]>), ('BEC_centerx', <Parameter 'BEC_centerx', value=0.0, bounds=[-inf:inf]>), ('BEC_centery', <Parameter 'BEC_centery', value=0.0, bounds=[-inf:inf]>), ('thermal_centerx', <Parameter 'thermal_centerx', value=0.0, bounds=[-inf:inf]>), ('thermal_centery', <Parameter 'thermal_centery', value=0.0, bounds=[-inf:inf]>), ('BEC_sigmax', <Parameter 'BEC_sigmax', value=inf, bounds=[-inf:inf], expr='3 * thermal_sigmax - deltax'>), ('BEC_sigmay', <Parameter 'BEC_sigmay', value=1.0, bounds=[0:inf]>), ('thermal_sigmax', <Parameter 'thermal_sigmax', value=1.0, bounds=[0:inf]>), ('thermal_sigmay', <Parameter 'thermal_sigmay', value=-inf, bounds=[-inf:inf], expr='thermalAspectRatio * thermal_sigmax'>), ('deltax', <Parameter 'deltax', value=-inf, bounds=[0:inf]>), ('thermalAspectRatio', <Parameter 'thermalAspectRatio', value=-inf, bounds=[0.8:1.2]>), ('condensate_fraction', <Parameter 'condensate_fraction', value=0.5, bounds=[-inf:inf], expr='BEC_amplitude / (BEC_amplitude + thermal_amplitude)'>)])",
|
||
|
|
"text/html": "<table><tr><th> name </th><th> value </th><th> initial value </th><th> min </th><th> max </th><th> vary </th><th> expression </th></tr><tr><td> BEC_amplitude </td><td> 1.00000000 </td><td> None </td><td> 0.00000000 </td><td> inf </td><td> True </td><td> </td></tr><tr><td> thermal_amplitude </td><td> 1.00000000 </td><td> None </td><td> 0.00000000 </td><td> inf </td><td> True </td><td> </td></tr><tr><td> BEC_centerx </td><td> 0.00000000 </td><td> None </td><td> -inf </td><td> inf </td><td> True </td><td> </td></tr><tr><td> BEC_centery </td><td> 0.00000000 </td><td> None </td><td> -inf </td><td> inf </td><td> True </td><td> </td></tr><tr><td> thermal_centerx </td><td> 0.00000000 </td><td> None </td><td> -inf </td><td> inf </td><td> True </td><td> </td></tr><tr><td> thermal_centery </td><td> 0.00000000 </td><td> None </td><td> -inf </td><td> inf </td><td> True </td><td> </td></tr><tr><td> BEC_sigmax </td><td> inf </td><td> None </td><td> -inf </td><td> inf </td><td> False </td><td> 3 * thermal_sigmax - deltax </td></tr><tr><td> BEC_sigmay </td><td> 1.00000000 </td><td> None </td><td> 0.00000000 </td><td> inf </td><td> True </td><td> </td></tr><tr><td> thermal_sigmax </td><td> 1.00000000 </td><td> None </td><td> 0.00000000 </td><td> inf </td><td> True </td><td> </td></tr><tr><td> thermal_sigmay </td><td> -inf </td><td> None </td><td> -inf </td><td> inf </td><td> False </td><td> thermalAspectRatio * thermal_sigmax </td></tr><tr><td> deltax </td><td> -inf </td><td> None </td><td> 0.00000000 </td><td> inf </td><td> True </td><td> </td></tr><tr><td> thermalAspectRatio </td><td> -inf </td><td> None </td><td> 0.80000000 </td><td> 1.20000000 </td><td> True </td><td> </td></tr><tr><td> condensate_fraction </td><td> 0.50000000 </td><td> None </td><td> -inf </td><td> inf </td><td> False </td><td> BEC_amplitude / (BEC_amplitude + thermal_amplitude) </td></tr></table>"
|
||
|
|
},
|
||
|
|
"execution_count": 12,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"params"
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": false,
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2023-07-12T15:04:30.916669900Z",
|
||
|
|
"start_time": "2023-07-12T15:04:30.868798100Z"
|
||
|
|
}
|
||
|
|
}
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"source": [
|
||
|
|
"## Fitting class"
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": false
|
||
|
|
}
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 27,
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"from lmfit import Model"
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": false,
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2023-07-12T16:18:33.140123400Z",
|
||
|
|
"start_time": "2023-07-12T16:18:33.119652Z"
|
||
|
|
}
|
||
|
|
}
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 28,
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"\n",
|
||
|
|
"def polylog(power, numerator):\n",
|
||
|
|
"\n",
|
||
|
|
" order = 2\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",
|
||
|
|
"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",
|
||
|
|
"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,\n",
|
||
|
|
" BEC_sigmax=1.0, BEC_sigmay=1.0, thermal_sigmax=1.0, thermal_sigmay=1.0):\n",
|
||
|
|
"\n",
|
||
|
|
" return ThomasFermi_2d(x=x, y=y, centerx=BEC_centerx, centery=BEC_centery,\n",
|
||
|
|
" amplitude=BEC_amplitude, sigmax=BEC_sigmax, sigmay=BEC_sigmay\n",
|
||
|
|
" ) + polylog2_2d(x=x, y=y, centerx=thermal_centerx, centery=thermal_centery,\n",
|
||
|
|
" amplitude=thermal_amplitude, sigmax=thermal_sigmax, sigmay=thermal_sigmay)\n",
|
||
|
|
"\n",
|
||
|
|
"\n",
|
||
|
|
"class DensityProfileBEC2dModel(Model):\n",
|
||
|
|
"\n",
|
||
|
|
" fwhm_factor = 2*np.sqrt(2*np.log(2))\n",
|
||
|
|
" height_factor = 1./2*np.pi\n",
|
||
|
|
"\n",
|
||
|
|
" def __init__(self, independent_vars=['x', 'y'], prefix='', nan_policy='raise',\n",
|
||
|
|
" **kwargs):\n",
|
||
|
|
" kwargs.update({'prefix': prefix, 'nan_policy': nan_policy,\n",
|
||
|
|
" 'independent_vars': independent_vars})\n",
|
||
|
|
" super().__init__(density_profile_BEC_2d, **kwargs)\n",
|
||
|
|
" self._set_paramhints_prefix()\n",
|
||
|
|
"\n",
|
||
|
|
" def _set_paramhints_prefix(self):\n",
|
||
|
|
" # self.set_param_hint('BEC_sigmax', min=0)\n",
|
||
|
|
" self.set_param_hint('deltax', min=0)\n",
|
||
|
|
" self.set_param_hint('BEC_sigmax', expr=f'3 * {self.prefix}thermal_sigmax - {self.prefix}deltax')\n",
|
||
|
|
"\n",
|
||
|
|
" self.set_param_hint('BEC_sigmay', min=0)\n",
|
||
|
|
" self.set_param_hint('thermal_sigmax', min=0)\n",
|
||
|
|
" # self.set_param_hint('thermal_sigmay', min=0)\n",
|
||
|
|
" self.set_param_hint('BEC_amplitude', min=0)\n",
|
||
|
|
" self.set_param_hint('thermal_amplitude', min=0)\n",
|
||
|
|
"\n",
|
||
|
|
" self.set_param_hint('thermalAspectRatio', min=0.8, max=1.2)\n",
|
||
|
|
" self.set_param_hint('thermal_sigmay', expr=f'{self.prefix}thermalAspectRatio * {self.prefix}thermal_sigmax')\n",
|
||
|
|
"\n",
|
||
|
|
" # self.set_param_hint('betax', value=0)\n",
|
||
|
|
" # self.set_param_hint('BEC_centerx', expr=f'{self.prefix}thermal_sigmax - {self.prefix}betax')\n",
|
||
|
|
"\n",
|
||
|
|
" self.set_param_hint('condensate_fraction', expr=f'{self.prefix}BEC_amplitude / ({self.prefix}BEC_amplitude + {self.prefix}thermal_amplitude)')\n",
|
||
|
|
"\n",
|
||
|
|
" def guess(self, data, x, y, negative=False, pureBECThreshold=0.5, noBECThThreshold=0.0, **kwargs):\n",
|
||
|
|
" \"\"\"Estimate initial model parameter values from data.\"\"\"\n",
|
||
|
|
" fitModel = TwoGaussian2dModel()\n",
|
||
|
|
" pars = fitModel.guess(data, x=x, y=y, negative=negative)\n",
|
||
|
|
" pars['A_amplitude'].set(min=0)\n",
|
||
|
|
" pars['B_amplitude'].set(min=0)\n",
|
||
|
|
" pars['A_centerx'].set(min=pars['A_centerx'].value - 3 * pars['A_sigmax'],\n",
|
||
|
|
" max=pars['A_centerx'].value + 3 * pars['A_sigmax'],)\n",
|
||
|
|
" pars['A_centery'].set(min=pars['A_centery'].value - 3 * pars['A_sigmay'],\n",
|
||
|
|
" max=pars['A_centery'].value + 3 * pars['A_sigmay'],)\n",
|
||
|
|
" pars['B_centerx'].set(min=pars['B_centerx'].value - 3 * pars['B_sigmax'],\n",
|
||
|
|
" max=pars['B_centerx'].value + 3 * pars['B_sigmax'],)\n",
|
||
|
|
" pars['B_centery'].set(min=pars['B_centery'].value - 3 * pars['B_sigmay'],\n",
|
||
|
|
" max=pars['B_centery'].value + 3 * pars['B_sigmay'],)\n",
|
||
|
|
"\n",
|
||
|
|
" fitResult = fitModel.fit(data, x=x, y=y, params=pars, **kwargs)\n",
|
||
|
|
" pars_guess = fitResult.params\n",
|
||
|
|
"\n",
|
||
|
|
" BEC_amplitude = pars_guess['A_amplitude'].value\n",
|
||
|
|
" thermal_amplitude = pars_guess['B_amplitude'].value\n",
|
||
|
|
"\n",
|
||
|
|
" pars = self.make_params(BEC_amplitude=BEC_amplitude,\n",
|
||
|
|
" thermal_amplitude=thermal_amplitude,\n",
|
||
|
|
" BEC_centerx=pars_guess['A_centerx'].value, BEC_centery=pars_guess['A_centery'].value,\n",
|
||
|
|
" # BEC_sigmax=(pars_guess['A_sigmax'].value / 2.355),\n",
|
||
|
|
" deltax = 3 * (pars_guess['B_sigmax'].value * s2) - (pars_guess['A_sigmax'].value / 2.355),\n",
|
||
|
|
" BEC_sigmay=(pars_guess['A_sigmay'].value / 2.355),\n",
|
||
|
|
" thermal_centerx=pars_guess['B_centerx'].value, thermal_centery=pars_guess['B_centery'].value,\n",
|
||
|
|
" thermal_sigmax=(pars_guess['B_sigmax'].value * s2),\n",
|
||
|
|
" thermalAspectRatio=(pars_guess['B_sigmax'].value * s2) / (pars_guess['B_sigmay'].value * s2)\n",
|
||
|
|
" # thermal_sigmay=(pars_guess['B_sigmay'].value * s2)\n",
|
||
|
|
" )\n",
|
||
|
|
"\n",
|
||
|
|
" nBEC = pars[f'{self.prefix}BEC_amplitude'] / 2 / np.pi / 5.546 / pars[f'{self.prefix}BEC_sigmay'] / pars[f'{self.prefix}BEC_sigmax']\n",
|
||
|
|
" if (pars[f'{self.prefix}condensate_fraction']>0.95) and (np.max(data) > 1.05 * nBEC):\n",
|
||
|
|
" temp = ((np.max(data) - nBEC) * s2pi * pars[f'{self.prefix}thermal_sigmay'] / pars[f'{self.prefix}thermal_sigmax'])\n",
|
||
|
|
" if temp > pars[f'{self.prefix}BEC_amplitude']:\n",
|
||
|
|
" pars[f'{self.prefix}thermal_amplitude'].set(value=pars[f'{self.prefix}BEC_amplitude'] / 2)\n",
|
||
|
|
" else:\n",
|
||
|
|
" pars[f'{self.prefix}thermal_amplitude'].set(value=temp * 10)\n",
|
||
|
|
"\n",
|
||
|
|
" if BEC_amplitude / (thermal_amplitude + BEC_amplitude) > pureBECThreshold:\n",
|
||
|
|
" pars[f'{self.prefix}thermal_amplitude'].set(value=0)\n",
|
||
|
|
" pars[f'{self.prefix}BEC_amplitude'].set(value=(thermal_amplitude + BEC_amplitude))\n",
|
||
|
|
"\n",
|
||
|
|
" if BEC_amplitude / (thermal_amplitude + BEC_amplitude) < noBECThThreshold:\n",
|
||
|
|
" pars[f'{self.prefix}BEC_amplitude'].set(value=0)\n",
|
||
|
|
" pars[f'{self.prefix}thermal_amplitude'].set(value=(thermal_amplitude + BEC_amplitude))\n",
|
||
|
|
"\n",
|
||
|
|
" pars[f'{self.prefix}BEC_centerx'].set(\n",
|
||
|
|
" min=pars[f'{self.prefix}BEC_centerx'].value - 10 * pars[f'{self.prefix}BEC_sigmax'].value,\n",
|
||
|
|
" max=pars[f'{self.prefix}BEC_centerx'].value + 10 * pars[f'{self.prefix}BEC_sigmax'].value,\n",
|
||
|
|
" )\n",
|
||
|
|
"\n",
|
||
|
|
" pars[f'{self.prefix}thermal_centerx'].set(\n",
|
||
|
|
" min=pars[f'{self.prefix}thermal_centerx'].value - 3 * pars[f'{self.prefix}thermal_sigmax'].value,\n",
|
||
|
|
" max=pars[f'{self.prefix}thermal_centerx'].value + 3 * pars[f'{self.prefix}thermal_sigmax'].value,\n",
|
||
|
|
" )\n",
|
||
|
|
"\n",
|
||
|
|
" pars[f'{self.prefix}BEC_centery'].set(\n",
|
||
|
|
" min=pars[f'{self.prefix}BEC_centery'].value - 10 * pars[f'{self.prefix}BEC_sigmay'].value,\n",
|
||
|
|
" max=pars[f'{self.prefix}BEC_centery'].value + 10 * pars[f'{self.prefix}BEC_sigmay'].value,\n",
|
||
|
|
" )\n",
|
||
|
|
"\n",
|
||
|
|
" pars[f'{self.prefix}thermal_centery'].set(\n",
|
||
|
|
" min=pars[f'{self.prefix}thermal_centery'].value - 3 * pars[f'{self.prefix}thermal_sigmay'].value,\n",
|
||
|
|
" max=pars[f'{self.prefix}thermal_centery'].value + 3 * pars[f'{self.prefix}thermal_sigmay'].value,\n",
|
||
|
|
" )\n",
|
||
|
|
"\n",
|
||
|
|
" pars[f'{self.prefix}BEC_sigmay'].set(\n",
|
||
|
|
" max=5 * pars[f'{self.prefix}BEC_sigmay'].value,\n",
|
||
|
|
" )\n",
|
||
|
|
"\n",
|
||
|
|
" pars[f'{self.prefix}thermal_sigmax'].set(\n",
|
||
|
|
" max=5 * pars[f'{self.prefix}thermal_sigmax'].value,\n",
|
||
|
|
" )\n",
|
||
|
|
"\n",
|
||
|
|
" return update_param_vals(pars, self.prefix, **kwargs)"
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": false,
|
||
|
|
"ExecuteTime": {
|
||
|
|
"end_time": "2023-07-12T16:18:33.888182200Z",
|
||
|
|
"start_time": "2023-07-12T16:18:33.872222800Z"
|
||
|
|
}
|
||
|
|
}
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"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
|
||
|
|
}
|