analyseScript/Joschka/Parameter guessing_new_data.ipynb

926 lines
2.9 MiB
Plaintext
Raw Normal View History

{
"cells": [
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 16,
"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",
2023-07-20 20:34:19 +02:00
"import lmfit\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",
"\n",
"\n",
"\n",
"#dataSet"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:29:14.809415800Z",
"start_time": "2023-07-20T12:29:14.694793900Z"
}
}
},
{
"cell_type": "markdown",
"source": [
"## Some functions"
],
"metadata": {
"collapsed": false
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 19,
"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",
"\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",
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, order=15) + Thomas_Fermi_1d(x, x0_bec, amp_bec, sigma_bec)\n"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:30:00.165795400Z",
"start_time": "2023-07-20T12:30:00.055461200Z"
}
}
},
{
"cell_type": "markdown",
"source": [
"## Import Data"
],
"metadata": {
"collapsed": false
}
},
{
"cell_type": "code",
"execution_count": 3,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"//DyLabNAS/Data/Evaporative_Cooling/2023/04/17/import/*.h5\n"
]
}
],
"source": [
"\n",
"# import data\n",
"img_dir = '//DyLabNAS/Data/'\n",
"SequenceName = \"Evaporative_Cooling\" + \"/\"\n",
"folderPath = img_dir + SequenceName + '2023/04/17'# get_date()\n",
"\n",
"\n",
"shotNum = \"import\"\n",
"filePath = folderPath + \"/\" + shotNum + \"/*.h5\"\n",
"print(filePath)\n",
"#filePath = folderPath + \"/\" + shotNum + \"/2023-04-24_0009_Evaporative_Cooling_*0.h5\"\n",
"\n",
"# # load the data from HDF5 files\n",
"# dataSetDict = {\n",
"# dskey[groupList[i]]: read_hdf5_file(filePath, groupList[i])\n",
"# for i in [0] # range(len(groupList))\n",
"# }\n",
"\n",
"# selecte the data for centain camera\n",
"dataSet = read_hdf5_file(filePath, \"images/MOT_3D_Camera/in_situ_absorption\")\n",
"# flip the x and y axis"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:28:09.861106300Z",
"start_time": "2023-07-20T12:28:05.953128700Z"
}
}
},
{
"cell_type": "code",
"execution_count": 4,
"outputs": [],
"source": [
"dataSet = swap_xy(dataSet)\n",
"\n",
"# get the scan axis name of the shot\n",
"scanAxis = get_scanAxis(dataSet)\n",
"\n",
"# rechunck the data for parallel computing\n",
"# dataSet = auto_rechunk(dataSet)\n",
"\n",
"# calculate the absorption imaging\n",
"dataSet = imageAnalyser.get_absorption_images(dataSet)\n",
"\n",
"dataSet\n",
"\n",
"OD = dataSet[\"OD\"]\n",
"\n",
"OD_np = OD.to_numpy()\n",
"\n",
"#dataSet"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:28:14.590676200Z",
"start_time": "2023-07-20T12:28:09.859108500Z"
}
}
},
{
"cell_type": "code",
"execution_count": 5,
"outputs": [
{
"data": {
"text/plain": "<Figure size 700x1500 with 11 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAqQAAAXcCAYAAAD0kJyLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOy9f3yWV33//wzc2IY2FFIKUkCCC0q+GlYig2HakX6Ms2zLJljBFjtwxVmmTFu7flXmbjOtdV1b7Jc5qqOfwmqL1Ap+zD6japQgyxBkoKRb6ogjGCgF06RN2qaVwP3947zf13mf674Tghaw7Xk9Hnnkuu7r17nOdc77vM771ynK5XI5IiIiIiIiIiIiIs4Thp3vAkRERERERERERLy2EQlpRERERERERETEeUUkpBEREREREREREecVkZBGREREREREREScV0RCGhERERERERERcV4RCWlERERERERERMR5RSSkERERERERERER5xWRkEZERERERERERJxXREIaERERERERERFxXhEJaURERERERERExHlFJKQRERERERERZw2PPPIIGzZsON/FiPgNR1Fcyz4iIiIiIiLibOHyyy/nqaee4tSpU+e7KBG/wcic7wJEREREREREvLLQ1tbGhg0bePzxx+np6RmUbHZ1dQFw9dVXJ78tW7aMpUuXnvVyRrxyEAlpRERERERExJCxe/duampqePHFFykqKhrydT/4wQ8AyOVy1NTUnKXSRbxSEQlpRERERERExJBxxx138NJLLzFt2jT+4A/+gNGjRw9KTO+66y6ee+45stls8lskpBFpRB/SiIiIiIiIiCHj8ssv5/jx4/z85z/n8ssvH9L50Yc04nSIUfYRERERERERQ8bTTz/NZZddNiQyGhExVERCGhERERERETFk9Pf3c+GFF57vYkS8yhAJaURERERERMSQcSaBTL/K+RGvTcSgpoiIiIiIiIgh43/+53/IZIZOH3bv3s3JkyfPYokiXg2IQU0RERERERERERHnFdFkHxERERERERERcV4RTfYRERERERERQ8awYb+aLiumfYoYDJGQRkRERERERAwZw4YN40y9/aJ3YMTpEAlpRERERERExJDx/e9/f8BjJ0+e5Pjx4zQ3N/Pggw9y8uRJVq9ezbRp085hCSNeiYhBTREREREREREvO44ePco73vEOurq6aGlp4bLLLjvfRYr4DUYMaoqIiIiIiIh42TFhwgTWrl3LL37xC+64447zXZyI33BEDWlERERERETEWcNFF13EpEmT+OlPf3q+ixLxG4yoIY2IiIiIiIg4axgxYgSHDx8+38WI+A1HDGqKiIiIiIiIOCtoa2vjueeeiymfIk6LaLKPiIiIiIiIiIg4r4gm+4iIiIiIiIiIiPOKaLKPiIiIiIiIOGMcOHCARx55hJ/85Cd0d3fT398/4Lm5XI6mpqZzV7iIVxwiIY2IiIiIiIg4I3zmM5/hc5/7HKdOnaKoqOi050fvwIjTIRLSiIiIiIiIiCHj4Ycf5rOf/SwAY8eO5V3veheTJ0/mwgsvPM8li3glIxLSiIiIiIiIiCHjS1/6Erlcjrq6Or72ta9RXFx8vosU8SpAjLKPiIiIiIiIGDJGjRrFc889R0dHBxMnTjzfxYl4lSAS0oiIiIiIiIgh45JLLmHYsGF0d3ef76JEvIoQ0z5FREREREREDBkVFRU8//zzvPjii+e7KBGvIkRCGhERERERETFkLF++nJMnT7Jp06bzXZSIVxGiyT4iIiIiIiLijHDttdfS2NjI5s2b+V//63+d7+JEvAoQCWlERERERETEkFFfX8/JkydZu3YtTz/9NHPnzmXOnDlccsklg16XzWbPUQkjXomIhDQiIiIiIiJiyBg+fDi5XC5JiD9UGnHq1KmzWayIVzhiHtKIiIiIiIiIIaOmpiauvBTxsiNqSCMiIiIiIiIiIs4rYpR9RERERERERETEeUUkpBEREREREREREecVkZBGREREREREREScV8SgpoiIiIiIiIgh4+qrrz7ja3K5HE1NTS9/YSJeNYhBTRERERERERFDxvDhw097jlILmxoqpn2KGAxRQxoRERERERExZJwuwX1vby+7du2iubmZ0aNH85GPfGRIJDbitY2oIY2IOEP88pe/5B/+4R/YuHEj//3f/83JkyeZOnUqCxcu5NZbb6WkpOTXun9vby9333033/jGNzh48CDDhw/nTW96E+973/tYuXIlr3vd616mN4mIeHWiv7+f7du3s3fvXv7jP/6DvXv30tbWRi6XY+nSpaxfv/5leU7sq4PjBz/4Ae9+97v5vd/7Pb75zW+e7+JE/IYjEtKIiDNAd3c373jHO9i3bx8AF1xwAcOHD+eFF14AYMqUKWzfvp0pU6b8Svc/dOgQNTU1tLe3AzBy5EhOnjzJSy+9BMDMmTP53ve+x5gxY379l4mIeJWivb2dqVOnFjz2chHS2FeHhq9+9av86Z/+KV/+8pf54Ac/eL6LE/EbjBhlHxFxBliyZAn79u1j1KhRbNq0iRdeeIHnn3+e73znO0yYMIFDhw5RV1fHyZMnz/je/f391NXV0d7ezoQJE/jud7/L888/zwsvvMDXvvY1SkpK2LdvH+9///vPwptFRLy6UFJSwpVXXslHP/pRNmzYwBVXXPGy3Tv21aHjve99LyNGjGDdunXnuygRv+nIRUREDAmNjY05IAfkNm7cmHf83//935Pj69atO+P7r1u3Lrn+3//93/OOP/zww8nxxsbGX+kdIiJeCzh58mTu1KlTwW/z5s3LAbmlS5f+2vePffXMMHr06NyoUaPOdzEifsMRNaSvYNTU1FBUVMRnPvMZTpw4wd13382sWbMYPXo0RUVFSYqNz3zmMxQVFVFTUzPgvZqamigqKkoiIi3S13/ve9/jD//wD7nsssu48MILqaiooL6+nhdffHHA+3/7299m4cKFTJo0ide97nWMGjWKN77xjfz+7/8+d911F11dXb9OVZwTbNiwAYA3vvGNLF68OO/43Llzkzr653/+51/5/ldffTVz587NO/6+970vMUP+KvePOH+IffXcYtiwYQXr5+VC7KtDx5NPPklPT08SdR8RMRAiIX0V4MUXX6SmpoZbb72Vn/zkJ2dVGP/93/8973znO9m6dSv9/f388pe/5IknnuAzn/kMf/AHf1DQVP23f/u3XHPNNWzZsoUjR44wYsQIcrkcBw8e5Lvf/S5/9Vd/xf79+89KeV9OfPe73wXgmmuuGbB+58+fD8C//du/0dfXN+R7v/DCCzQ3Nwf3SKOoqIhrrrkGgO985ztDvnfEbw5iX33lI/bVoeOXv/wlf/EXf0Eul+Otb33r+S5OxG84YtqnVwG+9KUvAfDAAw+wePFiiouLefrpp1/2ge4nP/kJO3bs4BOf+AS33HILY8eOpaenh7vvvpu//du/Zdu2bWzYsIE/+7M/S645dOgQ9fX1ANxyyy18/OMf5/LLLwfg2WefpaWlhY0bN/7akelnG08//TRPPfUUwKCCVY+dOnWK1tZWqqqqhnT/1tbWJEffUO7/1FNP0dXVRWlp6ZDuH/GbgdhXX/mIfZWknQyE/v5+nnzySb797W/z5JNPAvDhD3/4XBQt4hWMSEhfBXjuuef41re+RV1dXfLbpZde+rI/55lnniGbzfKZz3wm+W3UqFHU19fz+OOPs3nzZjZu3BgMcrt27eLUqVO86U1v4u677w7ud8kll3DllVdy5ZVXnnFZBouiHQrS73E6qFAFmDhx4oDn2WNPPvnkkAnpr3r/V9Mg91pA7KtnjjPtq2cbsa86TfpQTPBFRUUMHz6cv/7rv2bJkiXnoGQRr2REQvoqwFve8pZggDtbuOCCC7j11lsLHvuTP/kTNm/enGfOGz16NODy9T3//PNcdNFFL0tZhg8fzvjx43/l6y+++OIzOr+3tzfZHjly5IDn2WP2mvN9/4jfDMS+euY40756thH7qvOJHoyQjhgxgtLSUmbMmMGiRYv4rd/6rXNYuohXKiIhfRWgurr6nDznLW95y4CDg5r20gEPs2fPZuz
},
"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[1], row=scanAxis[0])\n",
"plt.show()"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:28:17.402410Z",
"start_time": "2023-07-20T12:28:14.604021700Z"
}
}
},
{
"cell_type": "markdown",
"source": [
"## Guess center\n",
"\n",
"ToDo: Crop from center guess"
],
"metadata": {
"collapsed": false
}
},
{
"cell_type": "code",
"execution_count": 6,
"outputs": [],
"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",
"# thresh = gaussian_filter(thresh, sigma=0.1)\n",
"# thresh = np.where(thresh<0.1,0,1)\n",
"\n",
"#M = moments(thresh)"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:28:17.681942100Z",
"start_time": "2023-07-20T12:28:17.406531600Z"
}
}
},
{
"cell_type": "code",
"execution_count": 7,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[[76.98612316 83.06157849]\n",
" [75.11961722 83.70813397]\n",
" [75.59055118 84.61023622]\n",
" [74.91266376 80.6069869 ]\n",
" [76.890625 91.13541667]]\n",
"\n",
" [[77.40971719 84.36548223]\n",
" [75.53406593 86.14505495]\n",
" [76.15422886 86.86069652]\n",
" [75.89119171 88.92746114]\n",
" [77.05294118 80.42352941]]]\n"
]
}
],
"source": [
"center = calc_cen_bulk(thresh)\n",
"print(center)"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:28:17.682940500Z",
"start_time": "2023-07-20T12:28:17.434799900Z"
}
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 8,
"outputs": [
{
"data": {
"text/plain": "<Figure size 1200x800 with 10 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA98AAAKTCAYAAAAAHrsTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABYQ0lEQVR4nO3db4xc9X0/+s+4S9ci9a6FLVMHr3Gos/AgFguNK36bmsBtEW1AxbdQgX6i4JbEfhIhSNRiP8J9YucBFws1reorIXSDdY3bOhdRoUolNOFPNkqxYGOJimDHod5QYyMT7zY0THE494Gzk91h/u6cM3/Oeb2sEczMmfNn9/PZmfd8v3OmlCRJEgAAAEBmlvV6BwAAACDvhG8AAADImPANAAAAGRO+AQAAIGPCNwAAAGRM+AYAAICMCd8AAACQMeEbAAAAMiZ8AwAAQMbaCt8ffPBBbN26NcbHx+Pqq6+Om266KY4fPx4REWfOnIk/+IM/iE9/+tPxmc98Jl588cXK4xrdB4NC/VN0eoAiU/8UnR6AFCRt+PnPf548++yzyUcffZQkSZL89V//dfL5z38+SZIk+bM/+7Pk4YcfTpIkSf7t3/4tueyyy5L/+Z//aXofDAr1T9HpAYpM/VN0egA6V0qSJFlqcD9y5Ejccccd8dZbb8Vv/MZvxPHjx+M3f/M3IyLid37nd2LPnj3x+7//+w3vq1Yul6NcLleuf/TRR/Hee+/FqlWrolQqLXVXYUmSJIn/+q//ik9+8pOxbNniiSLqnyLQAxSZ+qfIGtV/hB4g/5r1wFIMdfLgxx57LG677bY4e/ZsfPjhh5WmiojYsGFDnDx5suF9tezduzf+6q/+qpPdgtTNzMzEunXrFt2m/ikSPUCRqX+KrFb9R+gBiqNeDyzFksP3nj174vjx4/H888/Hz3/+81R2JiJi165d8ZWvfKVyfXZ2NtavXx8RD0bEcGrbgdaUI2JfrFixYtGt6p/i0AMUmfqnyGrXf4QeoCjq98BSLSl8P/LII/HNb34zvvWtb8XFF18cF198cQwNDcU777xTeWfrrbfeivXr18eqVavq3lfL8PBwDA/Xaq7hiFi+lN2Fji2c6qT+KSI9QJGpf4qserq3HqBo0vzIQ9uT1x999NE4ePBgPPfcc7Fy5crK7X/yJ38Sf/d3fxcREa+88kq8/fbb8fnPf77pfTBI1D9FpwcoMvVP0ekB6ExbJ1z7yU9+EmNjY3HFFVdUht+Hh4fj+9//fpw+fTr+9E//NH784x/Hr//6r8fXv/71uPHGGyMiGt7XzNzcXIyOjkbEzvCOF933QUR8LWZnZ2Nubk79U0B6gCJT/xTZr+p/ZGREDqCAFvdAGjo623k3aDp6K/2ma4f6p/f0AEWm/imy3tZ/hB6g19LvgXTOmQ4AAADUJXwDAABAxoRvAAAAyJjwDQAAABkTvgEAACBjwjcAAABkTPgGAACAjAnfAAAAkDHhGwAAADImfAMAAEDGhG8AAADImPANAAAAGRO+AQAAIGPCNwAAAGRM+AYAAICMCd8AAACQMeEbAAAAMiZ8AwAAQMaEbwAAAMiY8A0AAAAZE74BAAAgY8I3AAAAZEz4BgAAgIwJ3wAAAJAx4RsAAAAyJnwDAABAxtoK3/fff39s2LAhSqVSTE9PR0TE2bNnY2JionIZHx+PoaGheO+99yIi4oYbbohPfepTlfv37duX+kFAt+gBikz9U2Tqn6LTA9C5oXYWvuOOO+Iv//Iv43d/93crt61atarSgBERjzzySLzwwgtxySWXVG7bt29fbN26teOdhV7TAxSZ+qfI1D9Fpwegc22F7+uvv77pMo8//njs3bt3yTtULpejXC5Xrs/NzS15XZC2rHtA/dPPPAdQZOqfotMD0LlUP/M9NTUVP/3pT+PWW29ddPvOnTtj06ZNceedd8aJEycarmPv3r0xOjpauYyNjaW5i5CpTntA/TPIPAdQZOqfotMD0Fyq4fvxxx+Pe+65J4aGfjWg/uSTT8Ybb7wRR48ejS1btnysIavt2rUrZmdnK5eZmZk0dxEy1WkPqH8GmecAikz9U3R6AJpLLXz/7Gc/i7//+7+PP//zP190+/w7VqVSKb785S/HiRMn4uzZs3XXMzw8HCMjI4suMAjS6AH1z6DyHECRqX+KTg9Aa1IL34cOHYqrr746rrrqqspt58+fj9OnT1euHz58OC699NJYtWpVWpuFvqEHKDL1T5Gpf4pOD0Br2jrh2o4dO+LZZ5+Nd955J26++eZYsWJFHD9+PCIuTDX50pe+tGj5crkct9xyS5TL5Vi2bFmsXr06nnnmmfT2HrpMD1Bk6p8iU/8UnR6AzpWSJEl6vRONzM3NxejoaETsjIjlvd4dCueDiPhazM7O9mTqk/qn9/QARab+KbLe1n+EHqDX0u+BVE+4BgAAAHyc8A0AAAAZE74BAAAgY8I3AAAAZEz4BgAAgIwJ3wAAAJAx4RsAAAAyJnwDAABAxoRvAAAAyJjwDQAAABkTvgEAACBjwjcAAABkTPgGAACAjAnfAAAAkDHhGwAAADImfAMAAEDGhG8AAADImPANAAAAGRO+AQAAIGPCNwAAAGRM+AYAAICMCd8AAACQMeEbAAAAMiZ8AwAAQMaEbwAAAMiY8A0AAAAZayt833///bFhw4YolUoxPT1duX3Dhg1x5ZVXxsTERExMTMShQ4cq9x07diwmJydjfHw8Nm/eHK+//npqOw/dpgcoMvVP0ekBikz9Q+faCt933HFHvPzyy3H55Zd/7L5Dhw7F9PR0TE9Px5133lm5fceOHbF9+/Z4880346GHHopt27Z1vNPQK3qAIlP/FJ0eoMjUP3SurfB9/fXXx7p161pe/syZM3HkyJG4++67IyLi9ttvj5mZmTh+/Hjdx5TL5Zibm1t0gX6RdQ+of/qZ5wCKznMAReY5ADqX2me+77nnnti0aVPcd9998e6770ZExMzMTKxduzaGhoYiIqJUKsX69evj5MmTddezd+/eGB0drVzGxsbS2kXIVBo9oP4ZVJ4DKDrPARSZ5wBoTSrh+8UXX4yjR4/Gq6++GqtXr4577713yevatWtXzM7OVi4zMzNp7CJkKq0eUP8MIs8BFJ3nAIrMcwC0biiNlaxfvz4iIi666KJ44IEHYnx8PCIixsbG4tSpU3H+/PkYGhqKJEni5MmTleVrGR4ejuHh4TR2C7omrR5Q/wwizwEUnecAisxzALSu45Hv999/P86dO1e5fvDgwbjmmmsiImLNmjVx7bXXxoEDByIi4vDhw7Fu3brYuHFjp5uFvqEHKDL1T9HpAYpM/UN72hr53rFjRzz77LPxzjvvxM033xwrVqyIf/mXf4nbb789fvGLX0SSJHHFFVfEN77xjcpj9u/fH9u2bYs9e/bEyMhIPPHEE6kfBHSLHqDI1D9FpwcoMvUPnSslSZL0eicamZubi9HR0YjYGRHLe707FM4HEfG1mJ2djZGRka5vXf3Te3qAIlP/FFlv6z9CD9Br6fdAamc7BwAAAGoTvgEAACBjwjcAAABkTPgGAACAjAnfAAAAkDHhGwAAADImfAMAAEDGhG8AAADImPANAAAAGRO+AQAAIGPCNwAAAGRM+AYAAICMCd8AAACQMeEbAAAAMiZ8AwAAQMaEbwAAAMiY8A0AAAAZE74BAAAgY8I3AAAAZEz4BgAAgIwJ3wAAAJAx4RsAAAAyJnwDAABAxoRvAAAAyJjwDQAAABkTvgEAACBjbYXv+++/PzZs2BClUimmp6cjIuKDDz6IrVu3xvj4eFx99dVx0003xfHjxyuPueGGG+JTn/pUTExMxMTEROzbty/VAyAdu2N3r3dhIOgBikz9U2Tqn6LTA9C5tsL3HXfcES+//HJcfvnli27fvn17/PCHP4wf/OAHcdttt8UXv/jFRffv27cvpqenY3p6Oh588MHO9xp6RA9QZOqfIlP/FJ0egM61Fb6vv/76WLdu3aLbli9fHl/4wheiVCpFRMR1110Xb7311pJ3qFwux9zc3KIL2TPy3Zqse0D90888B1Bk6p+i0wPQudQ/8/3
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": "<Figure size 1200x800 with 10 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA98AAAKTCAYAAAAAHrsTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9f3Dd53Xfib/YXta6bkBXkAxGhiTAXtABm4BWsAmUDZKAnoVnoe8EbqgmVJdqQ2WpNvTul5lQ3TC2/AeG328td7VZqV11a3VGypjbEbvityUz4aZia3QiJkESIS7iJZIlGiIxYAlRiFhQBaSGXd6E3z/O88Y5z3M/ICVbMsHV58zcucC9n/v5PD/OOc95znmf82y7evXqVWqqqaaaaqqppppqqqmmmmqqqaZ3jP7SjW5ATTXVVFNNNdVUU0011VRTTTX9P53qzXdNNdVUU0011VRTTTXVVFNNNb3DVG++a6qppppqqqmmmmqqqaaaaqrpHaZ6811TTTXVVFNNNdVUU0011VRTTe8w1Zvvmmqqqaaaaqqppppqqqmmmmp6h6nefNdUU0011VRTTTXVVFNNNdVU0ztM9ea7pppqqqmmmmqqqaaaaqqpppreYao33zXVVFNNNdVUU0011VRTTTXV9A5TvfmuqaaaaqqppppqqqmmmmqqqaZ3mN7S5vvrX/86P/ZjP8aHP/xhPvKRj/Cxj32M+fl5AJaXlxkbG2PXrl18z/d8D7/2a7+28btrfVdTTTcL1fxf07udahmo6d1MNf/X9G6nWgZqqultoKtvgdbX16/+yq/8ytW/+Iu/uHr16tWrTz311NWRkZGrV69evfpTP/VTVycmJq5evXr16vT09NXu7u6r//k//+frfldTTTcL1fxf07udahmo6d1MNf/X9G6nWgZqqulbp21Xr169+s1u3L/4xS/y4z/+4ywsLPAd3/EdzM/P853f+Z0ADA0N8dhjjzE6OnrN70r6xje+wTe+8Y2N///iL/6ClZUVbrvtNrZt2/bNNrWmmr4punr1Kmtra3zgAx/gL/2lHChS839N7waqZaCmdzPV/F/Tu5muxf9Qy0BN/8+n68nAN0ONb+XH//gf/2P+xt/4G7z22mtcuXJlQ6gAent7+cpXvnLN76ros5/9LMePH/9WmlVTTW87vfzyy9x5553ZZzX/1/RuoloGano3U83/Nb2bqYr/oZaBmt49tJkMfDP0TW++H3vsMebn5/l3/+7fsb6+/rY0BuBTn/oUjzzyyMb/b7zxBnfffTdwFHhPuPK7gd8HPgZ8obhLB7AG3Ap8IF0H1t0WcA/wHcBX0kvD8DHgN9JvI30gXf9H6feke48D/3tx7fuBP02/+ePQHtJ9fxLoBj5b0ftbgdeB+4D/K70/m37fSveWsvpQak8VvR/4HuA/hDZcj+4Bfg/v34eAj2Nj1wBeKO7fDP+/nvr2fuC7gN8O9ylpHDgL3J3+/yHgZPr7J9N3r6f/PwDsTfd6AXr+37BYjlsTEP81gAOprd3ASnre17HxBPhiRZvEL/fDX/4u+PN/FO75DeBJOjo6sl98+/n/f8H4qqTvw+ZIbXg/0JU+uwe4BfjD4rcNbPybwBXgD8J34tv3h99IbkQ/i/HWC/jYdWDzNwV8DfgoNq8d5Dz4oXT9R9P9/zDdvxObmwbGQx9K159N10om7wjPfZ1q+nDq0wFs3iX/3x3+Fh3C5KhBrkd07d342H4ck8dI5diQ2lveD5xXfxaTkd9ObSX1pYnJ9/vT/y1cTiT395HLYknN1Pa7gdPps3swWdC9Iy+8Pz1HYzmO9Tvqlq0iA1oD7k59uZYOBJuH14FFbFzEhx8CXsX0rT67NV3zEWxMAC4Ad6VnfATowfW2+F50D/Cl6/Tsh1K7X8d4fDndQ3Mf21NFmrv/EfiF8Lnack/6P7ZDshDH6gAmmx9IbfgANqb/Pn3/ZufyQ+m3v3Gd9kbS3H0A+G645fvg61VroegAcGaTNnVg8nCqeKbW0N/G52U7JpN/iumQP7LncysuW3r/KPCr6f4/DZwA/v4W4P+Xgf/NPvzLn4I/17iVNlBcE8F11N0Y35V6STSO6fUubE7FO98syZ75cLrncrpf2T5RB/l83I3Lh+ijmDyD9Uvtux/XdyWVfPjR9P6r2PyeTt83MH5pYLyjNt2RnvNDVPP6D2E2xhfZXCeJ7xP9tU/Bfyz5XuMV6SfTZ7fSbmuCrf/RprkHW3//AOdjMF7XvV/Hxl9646eBf5a+K/UabKb/4UbIwAzwX9BuPzeAH8bsiXKtvC/9/wKVvPeXPwV/TrpnOcdpb/AdPwd/pvX6C9j6sQ78AM4rYPPUg/Hcd+PjSrq+5FPZw/H3kK9Xasf3Fc/6vnTP/xKYTNffn74rZUHPFY+vY7bfq+kzrUHr+LoA8HP4WMe2Rz6XHv0Ncl66HmmsC9nYoMS/3/Ep+LNfxsb6j3CeBp9jrd2/Tb7fgs31WNXauBltB/6/lTLwzdI3tfn+hV/4BU6fPs3k5CTvfe97ee9730uj0eBP/uRPNjxbCwsL3H333dx2222bfldF73nPe3jPe95T8c3tGLOI/hDbWPw68N/hGzjSde/HBv8PgR/FmFNdXgbmgEFgB2b4A/xl4HuB6fS/BGMFY8oGPmQd2KTfEp7bCRwBjqff6LtbsM3ge4HdwHPF70jfL6XPfzW1bTX935W+a4XfLWNKqC/9fz69H8aMka8A+/DNggRrEFNgGgvSfT+OLWpd6f2P05jNhz6I1lI/Hwuf3ZI+1wZK924Ae4BL6fvXgB/Bx/g0JrhX0qsTV44rwGVsQwcsTmGCuh3f/PcCL+GCNpP+HgTOwAfvhC8/kf6fpX2+VoAxbJ7/Nfz5vy76a/2IUKcbw/8PAk/jiqWJ8cDvFe1dS6/vwvitGzM+b8H5fCVdI+X+WWweLoYx6cUcHy+ksdN8kT4TL96CGQjPYkZyN8Z7M9h8NrD51Zzqmf83xlvi+6/jfPGH+MbkFmzxVz9/L91z3b7bNgFXjwNDqf1r6flx86mx+TrmRNNmeTj95gLGo98DjGLj/HXgJzDeejBd/3S4Vxe2GD2d2jMUxui3sAU48lpXGs9T6Tdd5HIANvaSo1tT+86FPjyartd9+7G50kLZjfHH79mrmRbN9QVsA/FfY3L43+D6YjQ945b097/HN4JR120FGXhPaucg3P4J+OpLOD+Jj8DGFowPxjAdP42vAzswg+VsGoud2LzI6Bav9gFfxvjxD9NnGvsr6e8jGJ/N4boafE5i26KRPIfx4SFs7dL6Uq4Lfdg6IHn4UeA3MQOpA1iExqPQOp7uSXGPrwDfn9qgsVpI10g+lzEZVwbajwJNuG8QXjidni/9up5eTeA7U5/+W5xPR3H9IPnV2hPXuBXgr8LXnyTXDyX9Msa3F/H51f2upO9jf5upDZPp+f+R9s3EH+OO3ehMuyV8f0u6Rrp0K/D//86GzvjzaAv8Oqa7tFZfTW0fwHX5ErYp+63wu+H0/lp6//eYzXMCdzp8B7bergP7cce2eE3rRRUlHb0R5CD9f5Vchy2nexzF5lmyo/e7Mf08mdofSfr8X2N8MlV835faK32njQbYxkL9GU7t1XXiySup7bdgvP5z6fuLqT1j2FrXTNf8cWrvMjb+s7jjIdip//EE8GPp+1XcFpXsgI1tdCyVuqEDt4N7MF7dHtr7W8CngCfTPeM8Sf6/H5tn2ctXMHugPepcwr1vjAycxfhHYww+zr+FzaXG6auYLutO//9E+v0t+e/+/EQaG82Pvg9r+p+dxuziyH/vx20r0TomG8O4bQA2L1cxPtX1E9g4RxlbD/cp5/uPi89+D+Pvi/h+4z+Fv/el685gjru7MT6LfC5dKrm8J7X/J9L3S+Qydgsms3O4vTSf2vB+jJe+P30O8BR5cCIGK94D/A+YDi77Chu66s8uA38bkw3ZPwfTvfTbW/C50D5RY3Am3P8QrvO/go+9ghBxHRrD1rVjwP8HaJeBb4XeMnj9iSee4F/8i3/BF77wBf7aX/trG5//xE/8BE8//TQAv/M7v8PS0hIjIyPX/e7N0xo
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": "<Figure size 1200x800 with 10 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA98AAAKTCAYAAAAAHrsTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9f3Sd53XfiX6QHjWCazAhxIChIRqQAzpgYjAyxgHrIC6oFWQVWhOkoZZD3ZCT0AnZCXk7zJSahLXE3p6FmZE0l82Y7bAdsbOkGbFZUka8M2KW2Q7Zml0i6zKuGA/sEduSqdAItMQoZGQwJtzCHqLh/WPv79n7fc8LyrLlEBy9e60DnPP+fH7sZz/72fu799Nx8+bNm9RUU0011VRTTTXVVFNNNdVUU03fNfqe212Ammqqqaaaaqqppppqqqmmmmr6fzrVi++aaqqppppqqqmmmmqqqaaaavouU734rqmmmmqqqaaaaqqppppqqqmm7zLVi++aaqqppppqqqmmmmqqqaaaavouU734rqmmmmqqqaaaaqqppppqqqmm7zLVi++aaqqppppqqqmmmmqqqaaaavouU734rqmmmmqqqaaaaqqppppqqqmm7zLVi++aaqqppppqqqmmmmqqqaaaavouU734rqmmmmqqqaaaaqqppppqqqmm7zK9o8X3N77xDX7u536OD3/4w/zYj/0YP/3TP83MzAwAV69eZWJignXr1vGRj3yEf/7P/3nrvludq6mmO4Vq/q/pvU71GKjpvUw1/9f0Xqd6DNRU07tAN98BLSws3PzH//gf3/zTP/3Tmzdv3rx56NChm2NjYzdv3rx585d/+ZdvNpvNmzdv3rx57ty5m729vTf/7//7/37bczXVdKdQzf81vdepHgM1vZep5v+a3utUj4GaavrOqePmzZs3v92F+xe/+EU++clPMjs7y/vf/35mZmb4wR/8QQBGRkZ44oknGB8fv+W5Mn3zm9/km9/8Zuv3n/7pnzI3N8c999xDR0fHt1vUmmr6tujmzZvMz8/zgQ98gO/5niJQpOb/mt4LVI+Bmt7LVPN/Te9luhX/Qz0Gavp/Pr3dGPh2qPGd3Px3/+7f5a/8lb/CV7/6VW7cuNEaVAD9/f185StfueW5KnryySeZmpr6TopVU03vOr3++uvce++9hWM1/9f0XqJ6DNT0Xqaa/2t6L1MV/0M9Bmp679BSY+DboW978f3EE08wMzPDP/tn/4yFhYV3pTAAjz76KI888kjr99e+9jU++MEPAnuB701XPgC8BGwFnk/HG8BK4I+BHwD6gC/6uR/w4x8D1gC/D/xbohke8mf+cXrWIvBh4C6/fjE9azvwm6Ua6B0fBL6SjuHHd8DdPfCNJytq/yHgD7xOJ/z5f8fvv+bP/AO/9i8C/7JU71y2vwj8C7/vW6Gf9OfpGfcDDwL/CmvPf5je88H0HX/HH3v57wc+m55TJvXXh/33zxJt+KvAi0T7fxjrZ4Bn4BO/AZ8vt9tKoo4N4K/btfQB1/35C8D/CXRh/Zsp88sO6OuBS387lf+bwEG6uroKd/3Z8//fB96kvV0fBD7n38WrKzCe/4t+/N9hbaR7G1hfka4ViQc/jI0NXZ/ee/ej8I1LWD9+APhD+3/3dviGyvLTwNNYm/9B3MvHgEte7q8DXwa6sfF4HOuLDwA/5t//vl877/f/MHDEy6mxW26Tj3mdfh3rd/V5HjO671FsXN/w94t+Ehs/H0rv3gmU+a/q/Vv92FGKIlZ89qjX+wTwo378D/z/v/ZjV/1ajZM/8Gdu8ecu9e6G1/PHsLZTe8z5M9S/oiyb8Od/meh7WD5jQHPA/V5G/e8i+ijTFqyuv4/VU/UWf+S2kEz7MeAj/v0l//1/2v8P9cAfPJmuz4rjTxPjcCl6EOvfa/7ca/6Mea/DGortXiYfk9//KPxJ5kPVf9J/Zz5WXTPv/zomc1X/D2LjSuer2rKKftSfcXyJ8+U2IpXjA8BPwkd/CL5UNReKHqV9zIm6sDkjz8GSR7+KyadJ4LSfW+nl+VEvw1/0Z4j3v+z/J71OXXD3fwHfeBb4L28//3/f6/A1H9M/8ij8G7XLDmzOE+U5MdOHMLn2DyvO6TkvYePqRdp1jHdK6oufBHowXvvyLcq3EusP8cz9mNzKPLSFmK+60/e/julKVVTmwzxOHgUOp/I8ALyf4OmVXo6XgF+iuu0eBL7h12i8lel+gr+Ajz8KXyjztXTHTL/uZf8g7bqm3n0i/dac/y8p6sZ/EZtTwOoqffIr8Ocehf+4lFz7AeANquQ/3I454DX4c92pvKJOrF//tX8g+uuX/Pc/9OtK5ex7FL6Gy1SN/fxc4Ef+OvwbyYzPYry8QDtPfABrQ/HN30nnFjC59A9SOf566ZoP+P/3U5wLOoFxL5vu/UlMhxrH+vkPsTHc8Hdk+muYPtDpdRC//76X9Q/TczVugfc/Cl9XW2dZfD/Bzz/qdT7B0vK6ot1bvCtdq0wPAS+6rDsCbMbGx9F0jdq/E9jkzyvrNEs9X2M167pL0Urgb1SOgW+Xvi3/+W/+5m/y4osvcuLECd73vvdxzz330Gg0+KM/+qPWNbOzs3zwgx+85bkq+t7v/V5WrFjR+nzf932fn7kXuDt9vuD/XwT+VjrewDr5g5gS8a8whr/bf9+NMd7nsEb/CDDmnz8HjKZn/Rf+/yuYcG1gg+JurLMPl8r0QVj5iH+/mo53Ar3AT0FjAL5xwo+/P13z4xjDq06j2ILgbr9Xirau/1dex1/1T8OPP+nv+7+AfcWycTfw8+nY++PT+E9TOe4GZrCB/jmM2XP7/iF0/lVvFymPd3v5/49UFj3/Z7ABfTfwVf+te/+el+39fl9netdXsMn3iJ37/Fm/9iPY4L/f+011awC/6/3/sJXngdXAc8B9mKDJ/fURv2eb/34OLh1M5b8b+AsABajTbeH/7//PCMOHPncDn/fjKvNXMN74ceD7MEE77+d/3D/eh/yhXXvfk8BPeR/8oD9nLfArfmxrsd2+cQJ41c9dx3huDr7xHPA+f/95f+/7UlnFt/PABYxXdP3/CfyC990lTFg+Q/DRF/3znJf/K3buo/+NX/Or6T3/yvv27xFy4m5sosq88gvAOeCfAn/kbfNkase/ho25R/xz0M99wJ//ZHp2HlefxYS53tPwez7lxw5iCsLPYYaRL1r7cZef/3fYOPhVYpw0YM2TFGXGT3g5Mj9/1J/3DHzsSfvwk8SC5BupX+/2Mi3491/xPvnzfl58NgQsgzHA93o5N8DHnyRk9Q3gL6V2+GH/vIQZ3z7q9f8VP/8+rG//EOOB38Dk9QKmtB70z/cAZzF+eg7+4GB6h8v3ricxOfT50rt13Q8QMvUlv+8GJov+HTQe82tuAB3pPn3+MtbP4slfgj855+/wOj/wmJ//nH/y/f8Kk7f/STo24/+/keryIS/DDYxPfgn+O9XtJ4AfwvjrB1K9ftDfl+ff30jXXPXv+v3jmMJ2N8bvPwBfOlhqr/LnoJWlcM0PpX7/e6XrOzGF9kU//6q/a87amxtehh8C/or320X/6Bn/Jt7zjU5gEFgG/P+1v4/JjLvh36xI5X0O62P9lgL9895/v+Ln34cZInTdL/nnfoyXngP+KvB/QOPJ1Feau59M1+oZ91f0mT5z/v+L9sxWGy+ka/5yPL9zH6aQ69xFjIc+gvWpZGuau1rXHibGd/78FLAx/f4A8Lp9Op/ElPVhjOd/BZsvPpfq3EnMIUe9DZ4kdMpfxXStP6Y499yd+uQnvC5pvvjCEb/3I/6uX8UMC7nsmsM+S+hJeS79YX/33d4nH8D0FfXbiy6f7sZkmdpN/HHV2v8/HiRkzFXXB/SOeeT0KsO9b88c8A+9vJKp4qGb3k7vS8dnMP7uIgx1N9P5n7P/l16EPzmILRyz/Pwlv/4m/JszGB8cxOTITUyGZN1YPP9l74+DpXN3Y7rs3Xb/2icxvl1JyDTJqq+U7rtJyCXV4YuYznCGGGvXCfn+twjZfN7f8SuY3M4yby499y/770ft8/U5L59kscbURWzs/9fAj2D6zg/
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fsize = (12,8)\n",
"\n",
"vmax = 1\n",
"fig, ax = plt.subplots(shape[0],shape[1],figsize=fsize)\n",
"\n",
"for i in range(0,shape[0]):\n",
" for j in range(0,shape[1]):\n",
" # ax[i][j].pcolormesh(blurred[i][j], cmap='jet', vmin=0, vmax=1.5, alpha=1)\n",
" ax[i][j].pcolormesh(thresh[i][j], cmap='jet', vmin=0, vmax=1, alpha=1)\n",
" #ax[i][j].pcolormesh(cropOD[i][j], cmap='hot', vmin=0, vmax=1, alpha=1)\n",
" ax[i][j].plot(center[i,j,0],center[i,j,1], marker='x', markersize=12)\n",
"plt.show()\n",
"\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",
" ax[i][j].pcolormesh(blurred[i][j], cmap='jet', vmin=0, vmax=vmax, alpha=1)\n",
" #ax[i][j].pcolormesh(thresh[i][j], cmap='jet', vmin=0, vmax=1.5, alpha=1)\n",
" #ax[i][j].pcolormesh(cropOD[i][j], cmap='hot', vmin=0, vmax=1, alpha=1)\n",
" ax[i][j].plot(center[i,j,0],center[i,j,1], marker='x', markersize=12)\n",
"plt.show()\n",
"\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",
" ax[i][j].pcolormesh(cropOD[i][j], cmap='jet', vmin=0, vmax=vmax)\n",
" #ax[i][j].plot(max[i,j,1],max[i,j,0], marker='x', markersize=12)\n",
" ax[i][j].plot(center[i,j,0],center[i,j,1], marker='x', color='g', markersize=12)\n",
"plt.show()"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:28:21.174757200Z",
"start_time": "2023-07-20T12:28:17.453258900Z"
}
}
},
{
"cell_type": "markdown",
"source": [
"## Guess width"
],
"metadata": {
"collapsed": false
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 9,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[[30. 39.]\n",
" [15. 30.]\n",
" [ 9. 28.]\n",
" [10. 25.]\n",
" [10. 24.]]\n",
"\n",
" [[36. 44.]\n",
" [15. 32.]\n",
" [11. 28.]\n",
" [ 8. 24.]\n",
" [11. 22.]]]\n"
]
}
],
"source": [
"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",
"BEC_width_guess = guess_BEC_width(thresh, center)\n",
"\n",
"print(BEC_width_guess)"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:28:21.189199100Z",
"start_time": "2023-07-20T12:28:21.178971800Z"
}
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 10,
"outputs": [
{
"data": {
"text/plain": "<Figure size 1200x800 with 10 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA98AAAKTCAYAAAAAHrsTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABbvklEQVR4nO3dbYxc5Xkw/mucJWvRetfClqnBawx1Fj7EYqFxg9yawNMi2oCKU1wRPaLglsT+EpmXKMX+Uujzwc6jB9kiTavylxBSsB7jFicVFapUQp9AKIFiUdcSFcGOQ70BYyMT7zY0bDCc/wcyw+x4XnfPmZdzfj9pZM+ZM+dl97r2zDX3fe67lCRJEgAAAEBmFvT6AAAAACDvFN8AAACQMcU3AAAAZEzxDQAAABlTfAMAAEDGFN8AAACQMcU3AAAAZEzxDQAAABlTfAMAAEDGOiq+33vvvdiwYUOMj4/H5ZdfHtddd10cOXIkIiJOnjwZv/d7vxef+tSn4tOf/nQ8++yzlfc1ew0Ghfin6OQARSb+KTo5AClIOvDzn/88efLJJ5MPP/wwSZIk+cu//Mvkc5/7XJIkSfInf/InyX333ZckSZL867/+a3LhhRcmv/jFL1q+BoNC/FN0coAiE/8UnRyA+SslSZLMtXA/cOBAbNy4MV5//fX41V/91Thy5Ej82q/9WkRE/OZv/mbs2LEjfvd3f7fpa7VmZmZiZmam8vzDDz+Md955J5YsWRKlUmmuhwpzkiRJ/Nd//VdccMEFsWDB7I4i4p8ikAMUmfinyJrFf4QcIP9a5cBcDM3nzQ8++GDcdNNNcerUqXj//fcrSRURsWrVqjh27FjT1+rZuXNn/MVf/MV8DgtSNzk5GStWrJi1TPxTJHKAIhP/FFm9+I+QAxRHoxyYizkX3zt27IgjR47E008/HT//+c9TOZiIiO3bt8c999xTeT41NRUrV66MiLsjYji1/UB7ZiJidyxatGjWUvFPccgBikz8U2T14z9CDlAUjXNgruZUfD/wwAPx7W9/O7773e/GueeeG+eee24MDQ3FW2+9Vflm6/XXX4+VK1fGkiVLGr5Wz/DwcAwP10uu4YhYOJfDhXmr7uok/ikiOUCRiX+KrLa7txygaNK85aHjzuu7du2KvXv3xlNPPRWLFy+uLP+jP/qj+Ju/+ZuIiHjppZfijTfeiM997nMtX4NBIv4pOjlAkYl/ik4OwPx0NODaT37ykxgbG4tLLrmk0vw+PDwcL774Ypw4cSL++I//OH784x/HJz/5yfjmN78Z1157bURE09damZ6ejtHR0YjYFr7xovvei4ivx9TUVExPT4t/CkgOUGTinyL7OP5HRkbUARTQ7BxIw7xGO+8GSUdvpZ90nRD/9J4coMjEP0XW2/iPkAP0Wvo5kM6Y6QAAAEBDim8AAADImOIbAAAAMqb4BgAAgIwpvgEAACBjim8AAADImOIbAAAAMqb4BgAAgIwpvgEAACBjim8AAADImOIbAAAAMqb4BgAAgIwpvgEAACBjim8AAADImOIbAAAAMqb4BgAAgIwpvgEAACBjim8AAADImOIbAAAAMqb4BgAAgIwpvgEAACBjim8AAADImOIbAAAAMqb4BgAAgIwpvgEAACBjim8AAADIWEfF99atW2PVqlVRKpXi4MGDERFx6tSpmJiYqDzGx8djaGgo3nnnnYiIuOaaa+Liiy+uvL579+7UTwK6RQ5QZOKfIhP/FJ0cgPkb6mTljRs3xp/92Z/Fb//2b1eWLVmypJKAEREPPPBAPPPMM3HeeedVlu3evTs2bNgw74OFXpMDFJn4p8jEP0UnB2D+Oiq+r7766pbrPPzww7Fz5845H9DMzEzMzMxUnk9PT895W5C2rHNA/NPPXAMoMvFP0ckBmL9U7/l+/vnn46c//WnceOONs5Zv27Yt1qxZE7fcckscPXq06TZ27twZo6OjlcfY2FiahwiZmm8OiH8GmWsARSb+KTo5AK2lWnw//PDDcdttt8XQ0McN6o8++mi8+uqrcejQoVi/fv1ZCVlr+/btMTU1VXlMTk6meYiQqfnmgPhnkLkGUGTin6KTA9BaasX3z372s/jbv/3b+NM//dNZy8vfWJVKpfjKV74SR48ejVOnTjXczvDwcIyMjMx6wCBIIwfEP4PKNYAiE/8UnRyA9qRWfO/bty8uv/zyuOyyyyrLzpw5EydOnKg8379/f5x//vmxZMmStHYLfUMOUGTinyIT/xSdHID2dDTg2pYtW+LJJ5+Mt956K66//vpYtGhRHDlyJCI+6mry5S9/edb6MzMzccMNN8TMzEwsWLAgli5dGk888UR6Rw9dJgcoMvFPkYl/ik4OwPyVkiRJen0QzUxPT8fo6GhEbIuIhb0+HArnvYj4ekxNTfWk65P4p/fkAEUm/imy3sZ/hByg19LPgVQHXAMAAADOpvgGAACAjCm+AQAAIGOKbwAAAMiY4hsAAAAypvgGAACAjCm+AQAAIGOKbwAAAMiY4hsAAAAypvgGAACAjCm+AQAAIGOKbwAAAMiY4hsAAAAypvgGAACAjCm+AQAAIGOKbwAAAMiY4hsAAAAypvgGAACAjCm+AQAAIGOKbwAAAMiY4hsAAAAypvgGAACAjCm+AQAAIGOKbwAAAMiY4hsAAAAypvgGAACAjHVUfG/dujVWrVoVpVIpDh48WFm+atWquPTSS2NiYiImJiZi3759ldcOHz4c69ati/Hx8Vi7dm288sorqR08dJscoMjEP0UnBygy8Q/z11HxvXHjxnjuuefioosuOuu1ffv2xcGDB+PgwYNxyy23VJZv2bIlNm/eHK+99lrce++9sWnTpnkfNPSKHKDIxD9FJwcoMvEP89dR8X311VfHihUr2l7/5MmTceDAgbj11lsjIuLmm2+OycnJOHLkSMP3zMzMxPT09KwH9Iusc0D8089cAyg61wCKzDUA5i+1e75vu+22WLNmTdxxxx3x9ttvR0TE5ORkLF++PIaGhiIiolQqxcqVK+PYsWMNt7Nz584YHR2tPMbGxtI6RMhUGjkg/hlUrgEUnWsAReYaAO1Jpfh+9tln49ChQ/Hyyy/H0qVL4/bbb5/ztrZv3x5TU1OVx+TkZBqHCJlKKwfEP4PINYCicw2gyFwDoH1DaWxk5cqVERFxzjnnxF133RXj4+MRETE2NhbHjx+PM2fOxNDQUCRJEseOHausX8/w8HAMDw+ncVjQNWnlgPhnELkGUHSuARSZawC0b94t3++++26cPn268nzv3r1xxRVXRETEsmXL4sorr4w9e/ZERMT+/ftjxYoVsXr16vnuFvqGHKDIxD9FJwcoMvEPnemo5XvLli3x5JNPxltvvRXXX399LFq0KP7pn/4pbr755vjggw8iSZK45JJL4lvf+lblPQ899FBs2rQpduzYESMjI/HII4+kfhLQLXKAIhP/FJ0coMjEP8xfKUmSpNcH0cz09HSMjo5GxLaIWNjrw6Fw3ouIr8fU1FSMjIx0fe/in96TAxSZ+KfIehv/EXKAXks/B1Ib7RwAAACoT/ENAAAAGVN8AwAAQMYU3wAAAJAxxTcAAABkTPENAAAAGVN8AwAAQMYU3wAAAJAxxTcAAABkTPENAAAAGVN8AwAAQMYU3wAAAJAxxTcAAABkTPENAAAAGVN8AwAAQMYU3wAAAJAxxTcAAABkTPENAAAAGVN8AwAAQMYU3wAAAJAxxTcAAABkTPENAAAAGVN8AwAAQMYU3wAAAJAxxTcAAABkTPENAAAAGeuo+N66dWusWrUqSqVSHDx4MCIi3nvvvdiwYUOMj4/H5ZdfHtddd10cOXKk8p5rrrkmLr744piYmIiJiYnYvXt3qidAOpI37+/1IQwEOUCRiX+KTPxTdHIA5q+j4nvjxo3x3HPPxUUXXTRr+ebNm+OHP/xh/Pu//3vcdNNN8aUvfWnW67t3746DBw/GwYMH4+67757/UUOPyAGKTPxTZOKfopMDMH8dFd9XX311rFixYtayhQsXxuc///kolUoREXHVVVfF66+/PucDmpmZienp6VkPsle64P5eH8JAyDoHxD/
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# %matplotlib notebook\n",
"fig, ax = plt.subplots(shape[0],shape[1],figsize=fsize)\n",
"\n",
"for i in range(0,shape[0]):\n",
" for j in range(0,shape[1]):\n",
" # ax[i][j].pcolormesh(blurred[i][j], cmap='jet', vmin=0, vmax=1.5, alpha=1)\n",
" ax[i][j].pcolormesh(thresh[i][j], cmap='jet', vmin=0, vmax=1.5, alpha=1)\n",
" #ax[i][j].pcolormesh(cropOD[i][j], cmap='hot', vmin=0, vmax=1, alpha=1)\n",
" ax[i][j].plot(center[i,j,0],center[i,j,1], marker='x', markersize=12)\n",
"plt.show()\n"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:28:22.222142300Z",
"start_time": "2023-07-20T12:28:21.241262800Z"
}
}
},
{
"cell_type": "markdown",
"source": [
"## Mask array"
],
"metadata": {
"collapsed": false
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 11,
"outputs": [
{
"data": {
"text/plain": "(2, 5, 200, 150)"
},
2023-07-20 20:34:19 +02:00
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"masked = np.where(thresh==0, cropOD, np.nan)\n",
"np.shape(masked)"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:28:22.372628Z",
"start_time": "2023-07-20T12:28:22.224112500Z"
}
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 12,
"outputs": [
{
"data": {
"text/plain": "<Figure size 1100x1100 with 10 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA5EAAAN2CAYAAABgkw2FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOy9fXycR3nu/7W9kr1ypBCFSDQysQoy2E3lA2pjp3Ubhx63JKUpB0rSnqQt6Y/XnhQOLaWHkraQtoTyGijvUBrSQ9LiFsghtEkbU2LAQAwVaQyVwQLkYAWkYCWRYq2tXUu/P+a5dV/zaGVr7bVXcub6fPTRs8/OzDMv19wzu3tfcy+bmZmZISEhISEhISEhISEhISFhAVje6AokJCQkJCQkJCQkJCQkLB2kD5EJCQkJCQkJCQkJCQkJC0b6EJmQkJCQkJCQkJCQkJCwYKQPkQkJCQkJCQkJCQkJCQkLRvoQmZCQkJCQkJCQkJCQkLBgpA+RCQkJCQkJCQkJCQkJCQtG+hCZkJCQkJCQkJCQkJCQsGAUGl2B42F6epoHH3yQ1tZWli1b1ujqJCxhzMzMMDExwfnnn8/y5afn+5PE34R6oRH8hcThhPoh2eCEpYzE34SljFPB30X/IfLBBx/kyU9+cqOrkXAG4fvf/z5r1qw5Lc9K/E2oN04nfyFxOKH+SDY4YSkj8TdhKaOe/F30HyJbW1uzqw8Bxezva9m93wb+Lru+EChl14elhIcJzVwl9x6y0oEJyf9wdv3gCdb24uz/A8A52fW3gGdJXe5jBUe5hq8DcCsbOcqzgfuyNM+Q53+L0F6yel4MfEfacGF2/U3gAnn2T2fXX8vKAxiVcq1PrVzDM7Ly9V5rVt+vACu4lWdylCfK+w8DV2bXfydlVyTNU7J6Wbnn4H19AdCeXZ8H3J1dPw34ttTrvuz6/Kw8a5OluYDQJyVJp+N4HoEX1wunTj38WbcB9+O8eJjAX4BPAb+YXX/yOCXqOCuMJyWOjwuq5DfonDgfnzcPEI/tedl1AZ9PDxPG5rusYJpruB8wjuu3XiFNXB97hnH6HOCLc2oXyv1eVu5FGRcfy959kPn74WnZ/28Tj4HhQrl+OoGH1l6be+D8NGia4+E8vK9sDmifWp11DLSvngxcd1r5C8rh3wdWZtfnAE/Nrr82J0+AzXuI7cE5OH++jY//g7l0BWL7F9VqnvvHmwdFoCx10iWwkkt7vpRVyezgl4CmzA6uOMZzDE9hBR1cw60A3MqPy1y4GO+7p0u9vo3b7fukrALO3YeI8Szgc1Wf7/wp4vP54SppF4Jq/XUBbmub8D5RvkPozxLwfxpkg20P8U28HUWqz11do/JrSTVcAdwpr/NcsjUrb/dsfjwN+EF2Ha+/4XU1Xmv/Pg3f92jffhMwu/kd4FBmj88F1mZpxpi7Htj8fDjj/X8C53ErncJf22fch8/h7+bKsD6sZHW0tf4r0TOgI7vWfla+P4zPR12LJnDbXsk9X1FtPwAxP/Pzw/Ztn8P5cAS4qUH8fXNWh/PkXav/FcAdcl/rXq2tMHd+KqqtkeA8PA/vq2/n0mhf6z7ufJzbE/gcrOBcGpVnnoPz8qcxW+l7C7PDU/O0QXEe8FjG5e8Bj0hee/a3iOdeAZ/H5+fey6eD2JY8BZ/P+TXiQmxeBjwj+/8NvH/B+/cc5vbxsWC2xPI/JM+4l3rzd9F/iPSf74tAC/FEb5XrFsDS6k/+JeLNCJJnFb5wt+BGWD9w1oLVubpaWXb/CLCKo8DfzW5KLV9RrjWv1aWcS7cql26++3at/abtK8t1C3GfIPX9WcIGwcoylHCDrvXVRTRfbhGfWNpXq4nHs9p1Ee/P/P1VwEyV9sZ1Pp0uIf4sq6tuBrTfzpLrY0HHWWGvZzg+8n2TL0fHSZ9n12V8zHRulbB2zuW4QsfNnmPPqMYFRyj3WXKnJbtr+efrh2pzpVTlfQhjkW+vcW6+fi9zfGi/2xzQ+WJ1ztslzXN6+Rs/b2WuLtqn1aB1V3uQz6vjn/8QOV//Hu/+fPNgFbBC6nSsD5Fq58rCab1/vPnWwlHOFs6Oynt5e2dtydtwg861fJ9Xny/xGMzH/VpQrb+U100sXhtsvFuFtyO/3iFpdY06nl0+K5cmz6Vq46nl6jjF6294XY3X+fy61mh+s5tbCR8YLa+lUzukZYf3fP3vAfbO06b52md9WMnSqB1VLlbjdX5t0TTaV/bsSpV2aF3mW3OQa31evq4+ZxrHX+0HqxvM5V+1fq42xsfqL5hrJ9T+L6Tc/Ppl3C4Tf4jUfWq1eeflzLXDC3HLLAKVLO+zgWHJq9zVuacfIou59/Lp8vm13dX2Iqtyry1/NY7m0x8P+b2mzs/wJXA9+ZsO1klISEhISEhISEhISEhYMBb9L5FzsRvoyq715/sx/JP3MP4rT4HwE3M1lxV1E9kD9GbXvcBdx6lHB/4thX37cFCeb990rAf6s+tRoJdlTLOB7wAPMMATmWEM2CJtsjptyeoFoc275PlbgZ3yDPuGuwcYkXT7sv9FYheQDrm25/UTXFysLncBXVl9VwKfzepbkvwl4Jbs+jJ57gRwILvek5VpddFv4/dnfwBX4y42+o3nPrkuSLvb5b6Nf0leK5qY++3w6YS5L1i9unD+dgPbc+nXZ/+Hibm7X651PK3dXcTc037UMrqkLu34nNot6Ybxca4A43JtdVLXS4BBgIwzPwIuYIDzmeEH8oxBSd+Ra9MeuTYu+PvLeC4buA04lHERoE/q0ibl6Leo2g/Wb8qRPXK/CGwGdkj6u+Ra29uZ/Z/A7Ye2AZyno7l7YzgnO/Bf+rV/BnEzXYtLy6nGGDFX2uW+YZx4zlmaHtwmgo+vjt9IltfyF3LX1ez5ernfg9sJRQ/x+BTxMRyX+ncTj0OP2O2VDNDBDE/E55qVBcG1zObzIMv4Nhv4SQAGmGKGR7L31J6X5HldxP1j/VYgnrMlnONalkLbUCaeE8ardpybHcQ8VZvaI9fD+DoxQDzONs9HiH9tGCGWmpxufJO5vyrleWS2RPt/GNiUXQ/iY7BWrj+V/Z/PBgxVqY/avWp2GsJ8UPsO3r/Kvb04/7RcgNaMu+EZvuewfcsEoS2az8qqiC1/MMtrvz3Y/O8i3nOYHVUe2XOsnTqHm3JtsTE4wNy9gpal9QTnus0XtT+6lo3idmaU2HZZP7cTj5nZ60bC3DzNXpXw/vkUcfts7dJ9qr1v/WtpIP6l7fnE+xFdh8tybXO7B++rdkJfG0eV1zrGSL1aidcSwwS+D9qN8WoZj7KB5cAEAzwp249aXfpyZV2R/d8FlDIu/wDYL1y2uZ7fH1aIuVTNhVXXKH1/mOq2tpswV2yOlCCT/cT7KnCOVnNvz5e7AZ+DpSyN1WuTlPWLwF9RTyzBD5FLH0Wm+CavBWA1r2OywfU5HkJ9rwasvs0NrlHCYkeRCt/kfQCs5m/qxvEih/kmb83KXfxzJ+HMQbCDxum/WDD3wly4Ksv31sTZhNOOImXhbm1reGzL0/qf0FgEPt4IwGpuqsmehryvyfImLtcD6UNkg/AQrbiWa/HjIc4haDoTEhaGh2jB9Wf1LHc1C9N+JiTUFw9F2rNa8j2h7nVJSKgFD0VaydOXNyGh3nho9gyJE8m7tPbeix3pQ2QDMMkqOng/savR4kWo7z24y05CwrExSTMd/BHB1aWe5bbQwZ+zsNNQExLqh2AH/4jYhW4h+Zrp4LPZq13AoXpXLSHhmJhkZcbdE8nbfMJ5ExLqjcDHm1jYYXbV8n6Y+d3HE2rFEvoQOUpYvFXnpYt53t9adQOqgVifS6saGtMx7CH2k9eTqlQzYfqHIWJ9nj6/W+q7BSf+GEHXCMFf2T6gteI+73uy/Pa8VnnOTlyTsk/uD0p9W4F1cm3tVo3oevzkpo4szV3y2upu9TBYjJlOXCuwS8rtk7Trs/fMJz/v123YzVw9DsyvCVQjYpo0y1PIlTVOY/U45p9u46QczGtHYGFGrprP/Lhc58swP3zTQ9p47Mf7rQf/cuMyYt3Exuz/AebqgRV6dLzVUbVEimPpVDW9aja
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots(shape[0],shape[1],figsize=(11,11))\n",
"\n",
"cut_factor = 1\n",
"\n",
"for i in range(0,shape[0]):\n",
" for j in range(0,shape[1]):\n",
" #print(np.nanmax(masked[i,j]))\n",
" # ax[i][j].pcolormesh(blurred[i][j], cmap='jet', vmin=0, vmax=1.5, alpha=1)\n",
" #ax[i][j].pcolormesh(masked[i][j], cmap='jet', vmin=0, vmax=0.5, alpha=1)\n",
" ax[i][j].pcolormesh(cropOD[i][j], cmap='jet', vmin=0, vmax=2.5, alpha=1)\n",
" #ax[i][j].pcolormesh(thresh[i][j], cmap='jet', vmin=0, vmax=1.5, alpha=1)\n",
" #ax[i][j].pcolormesh(cropOD[i][j], cmap='hot', vmin=0, vmax=1, alpha=1)\n",
" ax[i][j].plot(center[i,j,0],center[i,j,1], marker='x', markersize=12)\n",
" alpha=1\n",
" ax[i][j].hlines(center[i,j,1] - cut_factor * BEC_width_guess[i,j,1]/2, 0, 150, color='r',alpha=alpha,linestyles='dotted')\n",
" ax[i][j].hlines(center[i,j,1] + cut_factor * BEC_width_guess[i,j,1]/2, 0, 150, color='r',alpha=alpha,linestyles='dotted')\n",
"\n",
" ax[i][j].vlines(center[i,j,0] - cut_factor * BEC_width_guess[i,j,0]/2, 0, 200, color='r',alpha=alpha,linestyles='dotted')\n",
" ax[i][j].vlines(center[i,j,0] + cut_factor * BEC_width_guess[i,j,0]/2, 0, 200, color='r',alpha=alpha,linestyles='dotted')\n",
"\n",
" ax[i][j].set_xlim(25,125)\n",
" ax[i][j].set_ylim(45,145)\n",
"plt.show()"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:28:23.363914600Z",
"start_time": "2023-07-20T12:28:22.244413900Z"
}
}
},
{
"cell_type": "markdown",
"source": [
"## Try with not masked array\n"
],
"metadata": {
"collapsed": false
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 13,
"outputs": [],
"source": [
"shape = np.shape(masked)\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(masked[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(masked[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)])"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:28:23.534159700Z",
"start_time": "2023-07-20T12:28:23.369082300Z"
}
}
},
{
"cell_type": "markdown",
"source": [
"## Bimodal 1d Fit"
],
"metadata": {
"collapsed": false
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 14,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
2023-07-20 20:34:19 +02:00
"fitting time: 138.66639137268066 ms\n"
]
}
],
"source": [
"#[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] / 1.22 # 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')\n",
" #popt[i,j, 1], pcov = curve_fit(density_1d, y, Y_guess_og[i,j] , p0[i, j, 1], nan_policy='omit')"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:28:23.535125600Z",
"start_time": "2023-07-20T12:28:23.379370400Z"
}
}
},
{
"cell_type": "code",
2023-07-20 20:34:19 +02:00
"execution_count": 15,
"outputs": [
{
"data": {
"text/plain": "<Figure size 1200x800 with 10 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA9wAAAKPCAYAAACfPcLwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdeXxddZ34/9c5d783uTdp071J0y3dKE0LBSwoi4JVNgdcUFELssw4jN8Zf/MFGXS+4jIw21cdvyp1YGAUN2RRRmURK1igQEsplBa60TRpmzRdktx9O+f8/vicm4UmbdLcJbn3/Xw88rC9OffmXbkn97zP+/15fzTLsiyEEEIIIYQQQgiRV3qpAxBCCCGEEEIIIcqRJNxCCCGEEEIIIUQBSMIthBBCCCGEEEIUgCTcQgghhBBCCCFEAUjCLYQQQgghhBBCFIAk3EIIIYQQQgghRAFIwi2EEEIIIYQQQhSAs9QBjIRpmhw8eJDq6mo0TSt1OKICWZZFJBJh+vTp6Hrx71fJOSBKSd7/otLJOSAqmbz/RSUbzft/XCXcBw8epL6+vtRhCEFbWxszZ84s+s+Vc0CMBfL+F5VOzgFRyeT9LyrZqbz/x1XCXV1dDah/aDAYLHE0ohKFw2Hq6+t734vFJueAKCV5/4tKJ+eAqGTy/heVbDTv/3GVcOfaR4LBoJxooqRK1cok54AYC+T9LyqdnAOiksn7X1SyU3n/y9A0IYQQQgghhBCiACThFkIIIYQQQgghCkASbiGEEEIIIcQp27VrF6tWraKpqYmVK1eybdu24465//77aW5u7v2qq6vjqquuKkG0QhSXJNxCCCGEEEKIU3bzzTdz0003sXPnTm677TbWrFlz3DHXXXcdW7Zs6f2aOnUqn/70p4sfrBBFJgm3EEIIIYQQ4pR0dnayadMmrr32WgCuvvpq2tra2L1795DPefnll+ns7OSKK64oVphClIwk3EIIIYQQQohT0tbWxrRp03A61eZHmqbR0NBAa2vrkM+57777+MxnPoPL5RrymFQqRTgcHvAlxHgkCbcQQgghhBCiKGKxGL/4xS/4/Oc/f8Lj7rrrLkKhUO9XfX19kSIUIr8k4RZCCCGEEEKckvr6etrb28lmswBYlkVraysNDQ2DHv+rX/2KJUuWsHjx4hO+7u23305PT0/vV1tbW95jF6IYJOGuYJZllToEIYpG3u9CKJZlyfkghMibyZMns2LFCh588EEAHnnkEWbOnMm8efMGPf6+++47aXUbwOPxEAwGB3yJ/JHPgeKRhLtC3fjjTVz87T+TzBilDkWIgtvUcoxldz7NQxvl7riobBnD5EPfXc91D2wsdShCFNWvNrVxxjf+wGutXaUOpSytXbuWtWvX0tTUxN133839998PwA033MDjjz/ee9yOHTvYsmULn/jEJ0oVqgBe3XeMs/7pj3znmZ2lDqUiOEsdgCi+dNbkD9sPAbC5tYtVc+tKHJEQhfXFn79GOJnl1kfe4OMrZQ2YqFzbD4Z5uyPC2x0RsoaJ0yH33UVl+N8PvwHAR+/ZwJ5/+nCJoyk/CxYsYMOGDcc9fu+99x53XCQSKVZYYhDHYmmu/qH6b/X/1u3mbz/QVOKIyp980lagzkiy98+pjFnCSIQojqwpbVNCAJj9WgijqWwJIxGiNAzTIp2Vax9Ruf7r+b29f55U7SlhJJVDEu4KdCjcl3D3T76FKFeaVuoIhBgbUv0SjUhSEm5ROfonFi/vPVrCSIQonXg6y09e2tf7d1PWcReFJNwVqL2nL8k+FE6VMBIhikNDMm4hAKL9kmxJuEWlME2LY7E0K7Sd/JvrHvbuG3p/aCHK2av7uuhJZHr/HkvJLKdikIS7AnUMSLilwi3Kn1S4hVD6t5FLS7moFMfiadxmgkc9X+Ojjj8zb+9PSh2SECURTWYJEuOhwL/wccefiKWzmLLsruAk4S5jv9rUxnX3v3LcRVX/JFsq3KISSL4thBJJZgb9sxDl7HAkxXWOp3r/XhOWycyiMsXTBp93/p6zjC38i+s/sSyIy45FBScJdxn73w+/wZ92HB4wHAHe3VIuFW5R/jQpcQsBQCSVBVQ1Q1rKRaXojKQ4U9/R+/e6ZEvpghGihOIZg0n09P69jh5i0u1UcJJwV4CueHrA3wdWuCXhFuWvf75tyYAQUcGykaP8wX0r/+26206+hSh/hyMppmt9g9ImZw5AsucEzxCiPCXSWWZqh3v/frq+R5YXFYEk3BWoo1+SfSSaImvI9hiivPVPuNPyfhcV7KI9dzFfP8D5jjdIRw6f/AlClIHOSJLp2pGBD3a8WZpghCiheNqgSd/f+/dl+jtS4S4CSbgrTGL/Vv4legdPuG9jsdaCacGB7kSpwxKioPpPKY/LRE5RqWJHOa37T71/9XS/U8JghCienq6jBDV1rfOisVg9ePjtEkYkRGlY8WNM1bp6/366JhXuYpCEu0wZ/SYO9iYb6RiOn/4F79G3s0hv42HvN5lMF396u7NEUQpRHGl772ENk3TrJsimT/IMIcpQz8CtkAJhSbhFhehRFb24I0SLNRUAKyYdHqLyVId3D/j7ZK1btgYrAkm4y1Si38TB3nbajffhThym2woQ1mvwW3Gudz7BH946VJoghSiSeFoNivpn538y5ZcfgvX/XuqQhCi+cPuAvwZje4c4UIjy4okfBCDhn8YRggBkI1JsEJXHlVSzDAzNBUBIi0lLeRFIwl2m4v1OnqxhgmnCSz8E4FvZT/Pcwv8DwKcdf2TrOwfthESI8pTIGFyqv8zHnc+pB7Y+VNqAhCiF8IEBf52Y2FeiQIQoLn9C3WzKVs+gR6sBIBOWhFtUHkdaDQuM+GcCECImLeVFIAl3mYql+yrc8bQB7VsgcpAEXh43VhFcdinUNlKtJXgPr3MsJi22ojxlDJOMYbFS77der+cAGLIHsagwYVXl22nOAGByuvVERwtRNoIp1clnVk8n4Z4AgBWVlnJReZzpCAAJO+Gu1hIkEjLLqdAk4S5T/dtD4hkDdj4FwJ+NpaRws2haCBZ8GIAP6JtVUi5EGcotr5itdfQ9aKTg0LYSRSREaVh2wv2SqYZGTTY6QLbJExUglFETyrXgDLLeierPsSMneooQZcmdDQOQqprR+1g23l2iaCqHJNxlqn8CHU9lYZdKuJ8xl1PlcTKp2gNNqwG40LGFWCJVkjiFKLSEfS402gm3obvVNw5sKlVIQhSdYVq8tm07ANutWQA4MSAVKWVYQhSF31RJhjs4CcOnEm5nUhJuUXk8WfU73/JPJOmoAsCQhLvgJOEuU7F+a7KzySi0vw7AemMps+sCaJoGs1YRx0edFobDb5UqVCEKKpE2cJFlpqbaB1umqRtNHNpewqiEKK6j0RShjH0OWFNJWmpgDoljJYxKiMJLZQ2qrSgAnuqJ+CdMA8Cd7palRaLieLPqXNC8IdJONUBQS3Sd6CkiDyThLlP9W8pnJt4CyyTimUoHE5ldF1DfcLjY6VoIgPfgK6UIU4iCi6cNZmqHcWgWMctDe2CJ+kZUpvOLyhFOpJmqqeS63ZpAF9UAZKNHSxmWEAUXSWapQSUZ3uBELjljIYaltm+J98jngKgsPlOdCw5/DRm3SrhJdZcuoAox4oR7165drFq1iqamJlauXMm2bcevg9ywYQPNzc00NzezZMkSbr75ZlKpvpbl++67j/nz5zN37lxuvPFGMhm5w5hv8X576tUe3QLAO161bq834Qb2+E4DoOrwq8ULTogiSmSyve3k+6ypdDnUwBwi7Sd4lhDlJRLuJqCpz+Evf+xCui3VShjpkknNoryFExlCWgwAh7+Wc+dPoUdTicar23aVMjQhis7fm3DXknWH1J+T3SWMqDKMOOG++eabuemmm9i5cye33XYba9asOe6YZcuWsXHjRrZs2cLWrVvp7OzkBz/4AQB79+7lq1/9KuvXr2f37t0cOnSIH/3oR6P+h4iB+reUr9DVB8pvjqgBCXMm9SXcrf6lANQefa2I0QlRPIm02Ztwt1hTOKLlEm6pbIjKkYioSnYKNx9aMYeYQyUcUUm4RZkLJzK
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": "<Figure size 1200x800 with 10 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA98AAAKTCAYAAAAAHrsTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9f3Sd53XfiX6QHjWCazAhxIChIRqQAzpgYjAyxgHrIC6oFWQVWhOkoZZD3ZCT0AnZCXk7zJSahLXE3p6FmZE0l82Y7bAdsbOkGbFZUka8M2KW2Q7Zml0i6zKuGA/sEduSqdAItMQoZGQwJtzCHqLh/WPv79n7fc8LyrLlEBy9e60DnPP+fH7sZz/72fu799Nx8+bNm9RUU0011VRTTTXVVFNNNdVUU03fNfqe212Ammqqqaaaaqqppppqqqmmmmr6fzrVi++aaqqppppqqqmmmmqqqaaaavouU734rqmmmmqqqaaaaqqppppqqqmm7zLVi++aaqqppppqqqmmmmqqqaaaavouU734rqmmmmqqqaaaaqqppppqqqmm7zLVi++aaqqppppqqqmmmmqqqaaaavouU734rqmmmmqqqaaaaqqppppqqqmm7zLVi++aaqqppppqqqmmmmqqqaaaavouU734rqmmmmqqqaaaaqqppppqqqmm7zK9o8X3N77xDX7u536OD3/4w/zYj/0YP/3TP83MzAwAV69eZWJignXr1vGRj3yEf/7P/3nrvludq6mmO4Vq/q/pvU71GKjpvUw1/9f0Xqd6DNRU07tAN98BLSws3PzH//gf3/zTP/3Tmzdv3rx56NChm2NjYzdv3rx585d/+ZdvNpvNmzdv3rx57ty5m729vTf/7//7/37bczXVdKdQzf81vdepHgM1vZep5v+a3utUj4GaavrOqePmzZs3v92F+xe/+EU++clPMjs7y/vf/35mZmb4wR/8QQBGRkZ44oknGB8fv+W5Mn3zm9/km9/8Zuv3n/7pnzI3N8c999xDR0fHt1vUmmr6tujmzZvMz8/zgQ98gO/5niJQpOb/mt4LVI+Bmt7LVPN/Te9luhX/Qz0Gavp/Pr3dGPh2qPGd3Px3/+7f5a/8lb/CV7/6VW7cuNEaVAD9/f185StfueW5KnryySeZmpr6TopVU03vOr3++uvce++9hWM1/9f0XqJ6DNT0Xqaa/2t6L1MV/0M9Bmp679BSY+DboW978f3EE08wMzPDP/tn/4yFhYV3pTAAjz76KI888kjr99e+9jU++MEPAnuB701XPgC8BGwFnk/HG8BK4I+BHwD6gC/6uR/w4x8D1gC/D/xbohke8mf+cXrWIvBh4C6/fjE9azvwm6Ua6B0fBL6SjuHHd8DdPfCNJytq/yHgD7xOJ/z5f8fvv+bP/AO/9i8C/7JU71y2vwj8C7/vW6Gf9OfpGfcDDwL/CmvPf5je88H0HX/HH3v57wc+m55TJvXXh/33zxJt+KvAi0T7fxjrZ4Bn4BO/AZ8vt9tKoo4N4K/btfQB1/35C8D/CXRh/Zsp88sO6OuBS387lf+bwEG6uroKd/3Z8//fB96kvV0fBD7n38WrKzCe/4t+/N9hbaR7G1hfka4ViQc/jI0NXZ/ee/ej8I1LWD9+APhD+3/3dviGyvLTwNNYm/9B3MvHgEte7q8DXwa6sfF4HOuLDwA/5t//vl877/f/MHDEy6mxW26Tj3mdfh3rd/V5HjO671FsXN/w94t+Ehs/H0rv3gmU+a/q/Vv92FGKIlZ89qjX+wTwo378D/z/v/ZjV/1ajZM/8Gdu8ecu9e6G1/PHsLZTe8z5M9S/oiyb8Od/meh7WD5jQHPA/V5G/e8i+ijTFqyuv4/VU/UWf+S2kEz7MeAj/v0l//1/2v8P9cAfPJmuz4rjTxPjcCl6EOvfa/7ca/6Mea/DGortXiYfk9//KPxJ5kPVf9J/Zz5WXTPv/zomc1X/D2LjSuer2rKKftSfcXyJ8+U2IpXjA8BPwkd/CL5UNReKHqV9zIm6sDkjz8GSR7+KyadJ4LSfW+nl+VEvw1/0Z4j3v+z/J71OXXD3fwHfeBb4L28//3/f6/A1H9M/8ij8G7XLDmzOE+U5MdOHMLn2DyvO6TkvYePqRdp1jHdK6oufBHowXvvyLcq3EusP8cz9mNzKPLSFmK+60/e/julKVVTmwzxOHgUOp/I8ALyf4OmVXo6XgF+iuu0eBL7h12i8lel+gr+Ajz8KXyjztXTHTL/uZf8g7bqm3n0i/dac/y8p6sZ/EZtTwOoqffIr8Ocehf+4lFz7AeANquQ/3I454DX4c92pvKJOrF//tX8g+uuX/Pc/9OtK5ex7FL6Gy1SN/fxc4Ef+OvwbyYzPYry8QDtPfABrQ/HN30nnFjC59A9SOf566ZoP+P/3U5wLOoFxL5vu/UlMhxrH+vkPsTHc8Hdk+muYPtDpdRC//76X9Q/TczVugfc/Cl9XW2dZfD/Bzz/qdT7B0vK6ot1bvCtdq0wPAS+6rDsCbMbGx9F0jdq/E9jkzyvrNEs9X2M167pL0Urgb1SOgW+Xvi3/+W/+5m/y4osvcuLECd73vvdxzz330Gg0+KM/+qPWNbOzs3zwgx+85bkq+t7v/V5WrFjR+nzf932fn7kXuDt9vuD/XwT+VjrewDr5g5gS8a8whr/bf9+NMd7nsEb/CDDmnz8HjKZn/Rf+/yuYcG1gg+JurLMPl8r0QVj5iH+/mo53Ar3AT0FjAL5xwo+/P13z4xjDq06j2ILgbr9Xirau/1dex1/1T8OPP+nv+7+AfcWycTfw8+nY++PT+E9TOe4GZrCB/jmM2XP7/iF0/lVvFymPd3v5/49UFj3/Z7ABfTfwVf+te/+el+39fl9netdXsMn3iJ37/Fm/9iPY4L/f+011awC/6/3/sJXngdXAc8B9mKDJ/fURv2eb/34OLh1M5b8b+AsABajTbeH/7//PCMOHPncDn/fjKvNXMN74ceD7MEE77+d/3D/eh/yhXXvfk8BPeR/8oD9nLfArfmxrsd2+cQJ41c9dx3huDr7xHPA+f/95f+/7UlnFt/PABYxXdP3/CfyC990lTFg+Q/DRF/3znJf/K3buo/+NX/Or6T3/yvv27xFy4m5sosq88gvAOeCfAn/kbfNkase/ho25R/xz0M99wJ//ZHp2HlefxYS53tPwez7lxw5iCsLPYYaRL1r7cZef/3fYOPhVYpw0YM2TFGXGT3g5Mj9/1J/3DHzsSfvwk8SC5BupX+/2Mi3491/xPvnzfl58NgQsgzHA93o5N8DHnyRk9Q3gL6V2+GH/vIQZ3z7q9f8VP/8+rG//EOOB38Dk9QKmtB70z/cAZzF+eg7+4GB6h8v3ricxOfT50rt13Q8QMvUlv+8GJov+HTQe82tuAB3pPn3+MtbP4slfgj855+/wOj/wmJ//nH/y/f8Kk7f/STo24/+/keryIS/DDYxPfgn+O9XtJ4AfwvjrB1K9ftDfl+ff30jXXPXv+v3jmMJ2N8bvPwBfOlhqr/LnoJWlcM0PpX7/e6XrOzGF9kU//6q/a87amxtehh8C/or320X/6Bn/Jt7zjU5gEFgG/P+1v4/JjLvh36xI5X0O62P9lgL9895/v+Ln34cZInTdL/nnfoyXngP+KvB/QOPJ1Feau59M1+oZ91f0mT5z/v+L9sxWGy+ka/5yPL9zH6aQ69xFjIc+gvWpZGuau1rXHibGd/78FLAx/f4A8Lp9Op/ElPVhjOd/BZsvPpfq3EnMIUe9DZ4kdMpfxXStP6Y499yd+uQnvC5pvvjCEb/3I/6uX8UMC7nsmsM+S+hJeS79YX/33d4nH8D0FfXbiy6f7sZkmdpN/HHV2v8/HiRkzFXXB/SOeeT0KsO9b88c8A+9vJKp4qGb3k7vS8dnMP7uIgx1N9P5n7P/l16EPzmILRyz/Pwlv/4m/JszGB8cxOTITUyGZN1YPP9l74+DpXN3Y7rs3Xb/2icxvl1JyDTJqq+U7rtJyCXV4YuYznCGGGvXCfn+twjZfN7f8SuY3M4yby499y/770ft8/U5L59kscbURWzs/9fAj2D6zg/
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"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)\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",
"\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",
" ax[i][j].pcolormesh(cropOD[i][j], cmap='jet', vmin=0, vmax=vmax)\n",
" #ax[i][j].plot(max[i,j,1],max[i,j,0], marker='x', markersize=12)\n",
" ax[i][j].plot(center[i,j,0],center[i,j,1], marker='x', color='g', markersize=12)\n",
"plt.show()"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
2023-07-20 20:34:19 +02:00
"end_time": "2023-07-20T12:28:25.674215800Z",
"start_time": "2023-07-20T12:28:23.538256400Z"
}
}
},
{
"cell_type": "code",
"execution_count": 291,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[75.69178082 75.69178082 1.28576093 0.5510404 6.55737705 24. ]\n",
"[75.41570362 77.33387879 2.56706777 0.32435402 9.02534088 21.78912987]\n",
"\n",
"[77.85648148 77.85648148 1.26012635 0.54005415 6.55737705 24. ]\n",
"[77.87651572 81.11229756 2.46976235 0.33019909 9.77966756 21.30091787]\n",
"\n",
"[74.20467836 74.20467836 1.24055688 0.53166723 7.37704918 27. ]\n",
"[73.89357612 76.27664601 2.43121779 0.2744056 9.98643603 20.45775088]\n",
"\n",
"[77.33333333 77.33333333 1.34319384 0.5756545 5.73770492 21. ]\n",
"[77.11393649 81.4023359 2.8892113 0.26104137 8.62446159 18.7349952 ]\n",
"\n",
"[74.9516129 74.9516129 1.37694159 0.59011782 5.73770492 21. ]\n",
"[74.59730235 78.58391875 2.78343113 0.22862231 9.03244268 18.55359115]\n",
"\n",
"[78.32758621 78.32758621 1.11153483 0.47637207 6.55737705 24. ]\n",
"[78.50292398 82.4025832 2.35517042 0.20782941 9.17053528 15.99072519]\n",
"\n",
"[76.02054795 76.02054795 0.47105793 0.20188197 6.55737705 24. ]\n",
"[75.5176545 76.75506323 0.87438665 0.14660684 8.96570206 5.96922762]\n",
"\n",
"[76.22058824 76.22058824 0.47981987 0.20563709 7.37704918 27. ]\n",
"[76.17321865 76.98708716 0.61643543 0.23667781 9.40435385 5.34968882]\n",
"\n",
"[75.00787402 75.00787402 0.43698764 0.18728042 6.55737705 24. ]\n",
"[74.88995229 73.93081685 0.90813163 0.09113896 9.42575582 7.98491696]\n",
"\n"
]
}
],
"source": [
"for i in range(0,3):\n",
" for j in range(0,3):\n",
" print(p0[i,j])\n",
" print(popt[i,j])\n",
" print(\"\")"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-07-19T13:29:00.708064200Z",
"start_time": "2023-07-19T13:29:00.599583900Z"
}
}
},
2023-07-20 20:34:19 +02:00
{
"cell_type": "markdown",
"source": [
"## Fit Y"
],
"metadata": {
"collapsed": false
}
},
{
"cell_type": "code",
"execution_count": 17,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[1.29990174 0.83054958 0.81102892 0.78729 0.61898862]\n",
" [1.46785552 0.86056717 0.78554798 0.88165254 0.66979222]]\n"
]
}
],
"source": [
"S = np.zeros((shape[0], shape[1]))\n",
"for i in range(0, shape[0]):\n",
" for j in range(0, shape[1]):\n",
" S[i,j] = np.max(Y_guess_og[i,j])/(popt[i,j,2] + popt[i,j,3])\n",
"print(S)\n"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-07-20T12:29:33.071841100Z",
"start_time": "2023-07-20T12:29:33.020523300Z"
}
}
},
{
"cell_type": "code",
"execution_count": 20,
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.32201295482940323\n"
]
}
],
"source": [
"# without mathematical constraint\n",
"# Fix most\n",
"fitmodel = lmfit.Model(density_1d,independent_vars=['x'])\n",
"\n",
"p0_y = np.ones((shape[0], shape[1], 6))\n",
"popt_y = np.zeros((shape[0], shape[1], 6))\n",
"\n",
"\n",
"print(S[0,0]* popt[0,0,3]+S[0,0]* popt[0,0,2])\n",
"\n",
"result = []\n",
"start = time.time()\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,1], False,0, 200),\n",
" ('x0_th',center[i,j,1], False,0, 200),\n",
" ('amp_bec', S[i,j]* popt[i,j,2], True, 0, 1.3 * np.max(Y_guess_og[i,j])),\n",
" ('amp_th', S[i,j]* popt[i,j,3], False),\n",
" ('sigma_bec',BEC_width_guess[i,j,1], True, 0, 50),\n",
" ('sigma_th',popt[i,j,5], False)\n",
" )\n",
" res = fitmodel.fit(Y_guess_og[i,j], x=y, params=params)\n",
" temp_res.append(res)\n",
" result.append(temp_res)\n",
"stop = time.time()"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-07-20T12:30:04.185312900Z",
"start_time": "2023-07-20T12:30:04.030255800Z"
}
}
},
{
"cell_type": "code",
"execution_count": 22,
"outputs": [
{
"data": {
"text/plain": "<Figure size 1200x800 with 10 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA9wAAAKPCAYAAACfPcLwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOy9eXxcd3nv/z6zz0ia0S4vsizvsR3bshOHEAKEJWBocNqEQCihNZfEBsq97c3tL8ENoYQSkpbbS1kaYurULQ0khSzgUiAszkpMbCex4yWJl1iWZGvfZiTNPuf3x/fMaEaLtc0qPe/Xa16aOfOdM19Jc+acz/d5ns+j6bquIwiCIAiCIAiCIAhCWjHlegKCIAiCIAiCIAiCMBsRwS0IgiAIgiAIgiAIGUAEtyAIgiAIgiAIgiBkABHcgiAIgiAIgiAIgpABRHALgiAIgiAIgiAIQgYQwS0IgiAIgiAIgiAIGUAEtyAIgiAIgiAIgiBkAEuuJzAVYrEYFy5coKSkBE3Tcj0dYQ6i6zo+n48FCxZgMmV/vUqOASGXyOdfmOvIMSDMZeTzL8xlZvL5LyjBfeHCBRYtWpTraQgCzc3N1NbWZv195RgQ8gH5/AtzHTkGhLmMfP6Fucx0Pv8FJbhLSkoA9Yu63e4cz0aYi3i9XhYtWpT4LGYbOQaEXCKff2GuI8eAMJeRz78wl5nJ57+gBHc8fcTtdsuBJuSUXKUyyTEg5APy+RfmOnIMCHMZ+fwLc5npfP7FNE0QBEEQBEEQBEEQMoAIbkEQBEEQhBly6tQprrrqKlauXMnmzZs5fvz4qDHPPPMMTqeThoaGxM3v9+dgtoIgCEK2KKiUckEQBEEQhHxkx44dbN++nW3btvHYY4+xbds2Dh48OGrcqlWrOHz4cPYnKAiCIOQEiXALgiAIgiDMgI6ODg4dOsQtt9wCwI033khzczOnT5+e0X6DwSBerzflJgiCIBQWIrgFQRAEQRBmQHNzM/Pnz8diUYmDmqZRV1dHU1PTqLFnzpxh06ZNbN68mQceeOCi+73vvvvweDyJm7REEgRBKDwkpVwQBEEQBCELbNq0iZaWFjweDy0tLXz4wx+msrKSj33sY2OO37lzJ7fffnvicbwtjSAIglA4SIRbEARBEARhBixatIjW1lYikQgAuq7T1NREXV1dyji3243H4wGgtraWT3ziEzz//PPj7tdutydaIEkrJCGfmYxp4J49e1IMAysrK7nhhhtyMFtByC4iuAVBEIS0IC7NwlylurqaTZs28fDDDwPw+OOPU1tby/Lly1PGtba2EovFAPD5fPz85z9n48aNWZ+vIKSbuGngyZMnufPOO9m2bduoMZ/+9Kc5fPhw4jZv3jw++clPZn+ygpBlRHALgjAr0HU911OY80zmgguGXZrjN6fTmd2JzjLks58f7Nq1i127drFy5Uruv/9+9uzZA8Ctt97K3r17ASXE161bx4YNG7jyyiu59tpr+fSnP53LaRcc8nnPP6ZjGvjSSy/R0dHB1q1bszXNOYUcJ/mF1HAXKgOd8LO/gDO/g5pL4fp/hnmX5npWgpATdj7xGs+d7OKXf/VO3A5rrqczJ4lfcP36178G1AXXF77wBU6fPj0qyjdZgsEgwWAw8VgcmkczEIyw5Z+e421LKvjHj23I9XTmNKtWrWL//v2jtu/evTtx/wtf+AJf+MIXsjmtWYU/FGXLt55jU10Z3/x4Q66nIxhczDRwvO//hx56iE996lNYreOfs+UcMD16B0N86FvP88G1NdxzvWiDfEAi3IVIJAj/8cdw6imIRaD1MPzbh6GvOdczE4Sc8MiBZs73+fmvIxdyPZU5SyZcmsWheWJ+93o7Lb1+Hn+lJddTEYSM8tJb3az+8q841z3Ek6+ez/V0hBkwODjIo48+ymc+85mLjpNzwPT4jz+co80b4N/3n8v1VAQDEdyFyO+/Be3HoKgKPvkYzG+AQD/84v/L9cwEIaeYNS3XUxAmIO7S/Morr/Dkk0/y4IMP8uMf/3jMsTt37qS/vz9xa26WRcWR2C3Dp3F/KJrDmQhCZvn49/+Q6ykI4zBZ08A4P/nJT1i7di1r1qy56H7lHDA9huRckHeI4C40Al548bvq/pb7YcW1cMO/gGaCk7+E86/kdn6CkGWiseE6JatZvtJyRSZcmsWheWK0pEWmTl/wIiMFQRAyw2RNA+M89NBDE0a3Qc4B0yUQFsGdb8jVaaFx+IcQ7IfKVbDWaKVQtRLW3aTuv/Rg7uYmCFmmayDI1u++kHhsMUuEO1eIS3NuSI5qdw4EcjgTQcguYgqVX0zGNBDgzTff5PDhw3z84x/P1VRnNdGYTmP3YK6nIYxATNMKjcM/VD+vuA1MSeslV2yH1/4TTuyFD38DHJ7czE8Qssj/fepNjl8QE5V8YdeuXWzbto2vf/3ruN3ulAuurVu3snXrVh5//HG+973vYbFYiEQi3HTTTeLSPAOSUwclwi3MJaIxXRZZ84jJmAbGx/l8vmxNa87xN08c5Zk3O3M9DWEEUxbcp06d4s///M/p6urC4/Hwb//2b6xduzZlzP79+/nc5z4HQDgc5uqrr+bb3/42drudZ555hg996EOsWrUqZby0hZkEnW9C21EwWeHSG1OfW3iZinp3vQknn4L1H8vNHAUhizT3DqU8DkZiOZqJAOLSnAuGQpHE/Q4R3MIcIhzVsZhzPQtByC/+85DUuecjU04pn0yf1Q0bNnDw4EEOHz7M0aNH6ejoSHGilR6s0+TkU+rnkneBqzz1OU2D1dcZ436V3XkJQo4IhFMFtghuYa4hEW5hrhKKyve9IExEss+NkDumJLgn29je5XIl+uqFQiH8fn+KsctkCQaDeL3elNuc5pTqb8vKD479/Apj++nfQjQy9hhBmEWMNAYJilGIMMcYTIpwt3ulhluYO4RFcAvChMhxkh9MSXBPpc9qY2MjGzZsoLKyEo/Hw+c///nEc9KDdRoE+qHJSNVcce3YY2ovB2e5Gtv8UvbmJgg5YpTglgi3MMdINk073+fP4UwEIbuIkBCEVMY6JuS6KD/ImEt5fX09R44coa2tjWAwyBNPPAFID9Zpc+ZpiEWgYgWULx17jMkMKz6g7ktauTAH6B4MpTz+xlNvsvOJo/QPhXM0I0HILoPBYcHd1DN0kZGCMLsIRyRVVhCS6R0KjdomC1P5wZQE91Qb2wMUFxdz880388MfKndt6cE6Tc78Tv2MC+rxiEe/33o6s/MRhBzjD0XpG0NYP3KgiV8ca83BjAQh+/jDwynlF/oCcnElzBmkhlsQUukeEMGdr0xJcE+2z+rp06cJh9WFcCgU4sknn2T9+vWA9GCdNk1/UD+XvPPi4xZfpX62H4egtF0QZi8X6zM51iqvIMxGkiPc0ZjOBUkrF+YIIiQEIZWewdHXPiFJKc8LppxSPpnG9vv27WPjxo1s2LCBjRs3UlNTw9133w0okb5u3To2bNjAlVdeybXXXis9WCdiqAe6Tqr7i9528bHuBVBaB3oMWg5mfm6CkCPOdA6wSTvJ/zQ/wXrtTMpzAwExDRTmBsk13ADv/sYz4korzDoiY4hrEdyCkMrIMjuQ4yRfmHIf7sn0Wd2+fTvbt28f8/XSg3UaxA3QKleObgc2FnVvh74mFRVf9t7Mzk0QckT4zd/xE9s9mDWd/6U/yWfCf81zsQ0ADARFcAtzg7hL+QKPgwv9yqW80xdknseRy2kJQloZy/hJhIQgpNJltIb8o3XzeelsN10DITFNyxMyZpompJF4OvlE0e048XHx1wnCbCMa4V0nv4ZZU5E8qxblby0/wIISHz6JcAtzhHiE+/4b1ye2SQqhMNsY2ZECIBiO8a3fnuL3p7tyMCNByD/ixpm1ZU7sFjMA4ahkPOUDIrgLgXiEe7KCu+7t6mfLIenHLcxO3vgvKsJtdOslPLf1RQLWMpaZWnmP6TAggluYO8Qj3OVFNkpdVgBCUelHL8wuxorSPfZKC9/87Uk+uVvaoAoCDAvuugoXdouSeLIAmx+I4M53IiE4/4q6X3fl5F5TdQk4PBAehPZjmZubIOQ
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"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(y, Y_guess_og[i,j])\n",
" ax[i,j].plot(y, density_1d(y, **result[i][j].best_values), label = lab)\n",
" #ax[i,j].plot(y, thermal(y, **result[i][j].best_values))\n",
"\n",
"\n",
" #ax[i,j].legend(fontsize=10)\n",
" ax[i,j].set_facecolor('#FFFFFF')\n",
"plt.show()"
],
"metadata": {
"collapsed": false,
"ExecuteTime": {
"end_time": "2023-07-20T12:30:32.349861100Z",
"start_time": "2023-07-20T12:30:30.895968200Z"
}
}
},
{
"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
}