analyseScript/Joschka/Fitting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Import supporting package"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "outputs": [],
   "source": [
    "import lmfit\n",
    "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,
    "ExecuteTime": {
     "end_time": "2023-07-25T07:24:41.519442Z",
     "start_time": "2023-07-25T07:24:32.815578600Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "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):\n",
    "    res = np.exp(-0.5 * (x-x0)**2 / sigma**2)\n",
    "    return amp/1.643 * polylog(2, res)\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):\n",
    "    return thermal(x, x0_th, amp_th, sigma_th) + Thomas_Fermi_1d(x, x0_bec, amp_bec, sigma_bec)\n",
    "\n",
    "\n",
    "def polylog(pow, x):\n",
    "    order = 15\n",
    "    sum = 0\n",
    "    for k in range(1,order):\n",
    "        sum  += x ** k /k **pow\n",
    "    return sum\n",
    "\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",
    "\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, amp_bec=1.0, amp_th=1.0, x0_bec=0.0, y0_bec=0.0, x0_th=0.0, y0_th=0.0,\n",
    "                           sigmax_bec=1.0, sigmay_bec=1.0, sigmax_th=1.0, sigmay_th=1.0):\n",
    "    return ThomasFermi_2d(x=x, y=y, centerx=x0_bec, centery=y0_bec,\n",
    "                          amplitude=amp_bec, sigmax=sigmax_bec, sigmay=sigmay_bec\n",
    "                          ) + polylog2_2d(x=x, y=y, centerx=x0_th, centery=y0_th,\n",
    "                                          amplitude=amp_th, sigmax=sigmax_th, sigmay=sigmay_th)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-26T07:27:08.849310200Z",
     "start_time": "2023-07-26T07:27:08.756787700Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "outputs": [],
   "source": [
    "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/1.643  * polylog_int(np.exp( -((x-centerx)**2/(2 * (sigmax)**2))-((y-centery)**2/( 2 * (sigmay)**2)) ))"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-26T07:27:07.250614700Z",
     "start_time": "2023-07-26T07:27:07.165178200Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "outputs": [],
   "source": [
    "# Set up table for polylog\n",
    "\n",
    "def polylog_tab(pow, x):\n",
    "    order = 100\n",
    "    sum = 0\n",
    "    for k in range(1,order):\n",
    "        sum  += x ** k /k **pow\n",
    "    return sum\n",
    "\n",
    "x_int = np.linspace(0, 1.00001, 100000)\n",
    "\n",
    "poly_tab = polylog_tab(2,x_int)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T14:50:01.270231500Z",
     "start_time": "2023-07-25T14:50:00.931242800Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "outputs": [],
   "source": [
    "from scipy.interpolate import CubicSpline\n",
    "\n",
    "polylog_int = CubicSpline(x_int, poly_tab)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T14:50:02.455416900Z",
     "start_time": "2023-07-25T14:50:02.364988500Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.2525003487835937\n",
      "1.2525003487835928\n"
     ]
    }
   ],
   "source": [
    "x = 0.881\n",
    "print(polylog_int(x))\n",
    "print(polylog_tab(2, x))"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T14:50:03.008503Z",
     "start_time": "2023-07-25T14:50:02.935021500Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "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-25T13:18:51.814713700Z",
     "start_time": "2023-07-25T13:18:50.659769100Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 1000x900 with 10 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4gAAANiCAYAAAAjb+cqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9e3heVZ02jt9pk6YpJKUpJNC0TYQU2i+2QpWCFi0wdTjMixzkoCiCDlB8dVR0PCCOyOgPZEQcRkZlHE+visIlhxFR1CJUrQytFKRoqw2SQgMkQKAJNpSk7N8fe9/Putdnr/0kaXMqrPu6cuV5nr332ut8+Hzuda+KJEkSRERERERERERERERERLziMWm8IxARERERERERERERERExMRAXiBEREREREREREREREREA4gIxIiIiIiIiIiIiIiIiIkNcIEZEREREREREREREREQAiAvEiIiIiIiIiIiIiIiIiAxxgRgREREREREREREREREBIC4QIyIiIiIiIiIiIiIiIjLEBWJEREREREREREREREQEgLhAjIiIiIiIiIiIiIiIiMiw2y4Qv/3tb6OiogItLS0v63cS7e3tqKioQEVFBdrb28f8+Z0F33n33XeP2TsjRgexzY3t8zuL2OYihorYpsf2+Z1FbNMRERFjjd12gVgOt956Kz7zmc/g1ltvHe+ovGzx7//+7/jMZz6DBx54YLyjMuFx4YUXlgb48ZgUjQQefvhhfPvb3wYAbN68Gfvssw+OPfZY3HTTTQB2vc2tW7cO73znOzF79mxUV1djv/32wymnnII//elPQ3r+rrvuwimnnIL99tsP1dXVmD17Nt75zndi3bp1g773q1/9Ks4//3wsXrwY1dXVE7acYpsbOl4ubW7FihV41atehalTp+ba3K5i3bp1uO666wCkbZpt7le/+hWAcJvu7+/H17/+daxYsQKHH3445syZg5qaGkybNg0HHHAAzjrrLKxcuXJI79+2bRuOPfZYNDQ0YOrUqXjVq16FFStWoK2tbdhp+d3vfofJkycPupDavn176fOJJ56IvfbaC1VVVdhnn31w9NFH48tf/jK2bds27PfvLGKbHjpimy6P7du348tf/jKOPPJIzJgxA1OnTkVLSwvOO++8suNokiS455578KlPfQpHHXUUGhsbUVVVhenTp+O1r30tLr74YnR0dJR9d0tLS6lsiv6OPPLIwnjfdttteP/734/Xve51494mIwTJbopvfetbCYCkubk5d+2cc85JACTnnHPOmL1ztPHII48kABIAySOPPDLmz1s0NzcnAJJvfetbZe876KCDkoMOOii59957d/mduyN+9atfJRUVFaW8H4+6s6u4/fbbk2nTppXSUFFRkUyaNKn0/d3vfvcutbmvf/3rSWVlZSm86dOne3k2WL5deumlXtymT59e+l5ZWZl8/etfL3yW9dj+NTc3xza3m+Ll2Obq6upybe6ll17a6fBtm6uoqPDy7NJLLw226aeeesprJxUVFcmMGTOSyZMne7+fc845SX9/f/Dd3/zmN717J02alNTV1ZW+T5s2Lbn99tuHnJa+vr5k/vz5Xph33XVX8N6lS5d691VWVnr9BYCktbU12bRpU/D52KbHB7FNl8cTTzyRHHrooaWwqqqqkhkzZpS+V1dXJ9///veDz37uc5/Ltem99trLy++6urrkf/7nfwrfz3pcV1eXNDY2Bv/e8pa3BJ9dvnz5LrXJiNHDy9KDGDFxsHHjRmzcuBFLliwZ76iMObZt24bzzz8flZWVeN3rXjfe0dkpPPLIIzjjjDOwbds2tLa2AgDmzp2LrVu34tOf/jQA4Fvf+hYeeuihnQr/nnvuwYUXXoiBgQGcfPLJeOyxx/Dcc8/hqaeewooVK0r3FVkQb7zxRlx22WUAgBUrVuCpp57Cc889h8ceewwnn3wyBgYGcOGFF+Kee+4JPj9lyhQccsgheM973oNrr70WZ5999k6lYyIhtrmXT5tbunQp/vznP2Pr1q25NveFL3xhp8LXNnfooYcCSNu0trnLLrssSKGsrq7GP/3TP+GGG25Ae3s7tm/fju7ubrz44otYv3493va2twEAvvOd7+Cqq64Kvv9nP/tZ6fOll15aStvGjRvxhje8Adu2bcMZZ5yBRx55ZEjpueyyy0rPDoaBgYHS55/85CfYvn07nnvuOXR3d+PKK6/EtGnT0NbWhuOPPx4vvPDCkN4/FohtOrbpIiRJgre+9a24//77UVNTg69//evo6elBd3c3Hn/8cbzrXe/C9u3bce655+K+++7LPd/f34+6ujpceOGF+NWvfoW//e1vePbZZ/G3v/0NN910E+bOnYuenh6cfvrp2LBhQ9m4XHPNNXjyySeDf//zP/8TfKa/vx/Nzc34zGc+g3Xr1u02bfIVgfFeoe4sogdxbJ+3GKrl85WMD33oQwmA5JJLLinVyd3N8vnOd74zAZDsu+++yX/+53/m0nDBBReULJY70+aOPPLIBECycOHC5MUXX8xdf/WrX50ASCZPnpwMDAx41wYGBkr18Ljjjss9u3379tLzRx55ZPD9Nkx6I3dnD+IrGS+3Nvfss8/mrrPN1dXVJd3d3cMOX9vcf//3f+fy6Nhjj00AJHvsscew2/RLL72UvOENbyhZ/S26u7uTqVOnJgCSPffcM3h93333TQAk73znOwd933333ZdUVlYmra2tyc9+9rNBPYg33nhj2Tb5gx/8oHT9u9/9bu56bNNjj9imy+O2224r1cl///d/D95zxBFHJACSY445Jnft/vvvL/vOv/71r0lNTU0CIPnHf/zH4D27Uo9/85vfFLINkmTwNhkxehj1BeJdd92VnHbaacmsWbOSKVOmJDNnzkyOOeaY5Jvf/GZucmZxzz33JCeddFIyc+bMZOrUqcmBBx6YfPKTn0x6e3uDi7W77rorSBfTv6KBYyiw7/z973+fnH766cm+++6bTJkyJXnVq16VXHTRRYM28La2tuTCCy9MWltbk6lTpya1tbXJoYcemlx22WXJ1q1bg88MZWDasmVLcsEFFySzZ89OpkyZkjQ1NSXnnntusmnTphEb2JTOV/SnKMp3G5/29vbkvPPOS+bMmZNUV1cn+++/f3LJJZckzz//fOmZ9evXJ+94xzuS2bNnJ9XV1Ulra2vy2c9+NriwsO/64Ac/mPx//9//l+yxxx5JTU1NctBBByUf+MAHks2bN+90XpTDPffck0yaNCk58MADk76+vjEd2EayzZFm8slPfjLY5rTzHm6be/jhh0v3fOc73wne8/GPf7x0z69+9SuvzXFRCiD5yU9+Enz+29/+dumes846a9A2N9wFYmxzYbzS2txI4fnnny9Nxi677LLS77ZNsxzf/e53D6tNV1dXl5697rrrgm36S1/60i6NoxdddFECIJk6dWru2je+8Y1SGE1NTUmS5MfRvfbaq/S81kWLDRs2JPX19QmAZMqUKR59r6g/GKxN7tixo7Qw3nPPPWObNu+KbXr4KGrTCi3Hb37zm8MK/8ILL0yA1KBTVIY33XRTKfydKasTTjghAZAsWrQoeH00DR07duxIamtrEwDJBz7wgREPP6IYo7pA5EABOF6z7lU45phjkp6enuCz3/jGNzx+9vTp00sD4/z585Orr74611GsXr06aWxsLFkop06dmuNBr169eqfTo4Pp97///dIEdfr06V5cDz744KS3tzcYxg033OAN0rW1td73OXPmJH/6059yzw02MN13330e57ympibZc889EyC1St1www0jMrB94QtfSBobG0vpDXHOFUMZ2G666abSpKCurs6rI2984xuTF198MfnJT35SmgDYPWpnnnlmYXy/973veflbXV1d6qyZ/z//+c93Oj9CeOGFF5IFCxYkFRUVyd13350kSTJmA9tItjlOlAAkLS0thW1Owx9Om/va175Weq6zs7MwTrznxBNP9Nqc1oGiNtfZ2ZmbeJVrc8NZIMY2F8Yrrc2NJO64445SPq1ZsyZJknCb1vo8nDat5TBv3rxgm/71r3/tlfVw2vSOHTuSww8/vNQmLd7
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "imageAnalyser.center = (960, 890)\n",
    "imageAnalyser.span = (150, 150)\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-25T13:18:54.107015500Z",
     "start_time": "2023-07-25T13:18:51.857684700Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fitting time: 137.02940940856934 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,149, 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-25T13:18:54.284965700Z",
     "start_time": "2023-07-25T13:18:54.116202500Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 1500x1000 with 9 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABMAAAAMpCAYAAADipk9SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3hUVfrA8e+dnjaTTgJJCF2aFAFRFAEL2FCxrw0URVdXXXfVVde6P/uq664NV8W2VmzYewEVaYbeIaT3ZGYyk0y9vz/uJBJpSZjJpLyf55lHZubOve+MyZyc95zzHkVVVRUhhBBCCCGEEEIIIbopXbQDEEIIIYQQQgghhBAikiQBJoQQQgghhBBCCCG6NUmACSGEEEIIIYQQQohuTRJgQgghhBBCCCGEEKJbkwSYEEIIIYQQQgghhOjWJAEmhBBCCCGEEEIIIbo1SYAJIYQQQgghhBBCiG7NEO0A2iIYDFJSUkJCQgKKokQ7HCGE6PJUVcXpdNK7d290OhkTAWlrhBAi3KStaUnaGSGECK/WtjNdKgFWUlJCdnZ2tMMQQohup7CwkKysrGiH0SlIWyOEEJEhbY1G2hkhhIiMA7UzXSoBlpCQAGhvymq1RjkaIYTo+hwOB9nZ2c3fr0LaGiGECDdpa1qSdkYIIcKrte1Ml0qANU0Rtlqt0lgIIUQYyRKM30hbI4QQkSFtjUbaGSGEiIwDtTOyCF8IIYQQQgghhBBCdGuSABNCCCGEEEIIIYQQ3ZokwIQQQgghhBBCCCFEtyYJMCGEEEIIIYQQQgjRrUkCTAghhBBCCCGEEEJ0a5IAE0IIIYQQQgghhBDdmiTAhIgwjz/Au6uKqHN7ox2KEEKIbu77LZWsKqiNdhhCCCG6uJ+2VbFsZ020wxAirCQBJkSE3fPhBm54azWPf7012qEIIYToxqrqPVz64nIueX4Z/kAw2uEI0Wlde+215ObmoigKeXl5ez1mwYIFjB49uvmWmprKrFmzAMjPz0ev17d4fvv27R34DoSIrDq3l0sWLGP2gmV4/IFohyNE2BiiHYAQ3dn2ynreWF4IwNoie5SjEUII0Z3lV7kIBFWcHj/bK10MyUiIdkhCdEpnnXUWN910E0cdddQ+j5kzZw5z5sxpvj9ixAguuOCC5vsJCQn7TJ4J0dWtL3HgC6j4AgGKahsYkBYf7ZCECAtJgAkRQf/8fDOBoArA1op6VFVFUZQoRyWEEKI7Kq5raP73hlK7JMCE2IfJkye36fhffvmFiooKZs6c2a7reTwePB5P832Hw9Gu8wjRUTaW/vYzuqvaJQkw0W3IEkghIuTXglo+XVeGTgFFAXuDj6p6qQMmxP7IshQh2q+odrcEWIl0sIUIl+eff56LLroIo9HY/JjL5WL8+PGMHTuWe+65h0Bg38vE7r//fmw2W/MtOzu7I8IWot02tEiAuaMYiRDhJQkwISLkvV+LATh9dB/6JscCsLXCGc2QhOj0zjrrLJYsWULfvn33ecycOXPIy8trvmVkZOx1WUrTbcCAAR0RuhBRt/sMsPWSABMiLFwuF2+88QaXXXZZ82OZmZkUFxezfPlyvvrqKxYvXswjjzyyz3Pccsst2O325lthYWFHhC5Eu20s/a3PIgkw0Z1IAkyICCkJdUTG5SYzMF1bhrKtoj6aIQnR6U2ePJmsrKxWH3+wy1JAW5ricDha3IToiop3nwFW6kBV1ShGI0T38PbbbzN8+HCGDRvW/JjZbCY9PR2A5ORkLr30UhYvXrzPc5jNZqxWa4ubEJ2V1x9k226D9vnVrihGI0R4tSkBJktThGi9krpGADJtFgama+vmt5ZLAkyIcDrYZSkgS1NE97H7DLA6t49Se2MUoxGie3j++edbzP4CqKiowOfzAdogyrvvvsuYMWOiEZ4QYbe9sh5f4LcBlAKZASa6kTYlwGRpihCtV+YIJcASLQwKJcBkBpgQ4ROOZSkgS1NE96CqavMMsASztseR1AETYu/mzZtHVlYWRUVFTJ8+nYEDBwIwd+5cFi1a1Hzc5s2bycvL49xzz23x+iVLljBmzBhGjRrF2LFjycjI4LbbbuvQ9yBEpDQVwM9KigGgsNaNPxCMZkhChE2bdoHs6B1ThOiqGn0Bal2NZClVZMaAt1doBpgkwIQIm9YuS3nttde46aab9nkes9mM2WyOeLxCRFKt20eDT5vtOHlIGh+vKWV9iYPjhvWKcmRCdD7z58/f6+PPPfdci/tDhgzB6dyzfuusWbOaV7gI0d00JcCmDknnzRWFeP1BSu2NZIdqGgvRlUW0BtjBLk2RuiyiqyqzN3Kp/jOWmK/H+ng/hn80k7HKFqrqPdS6ZCdIIcJBlqUI8Zum2V9pCWbGZCdyrv5b0je9FOWohBBCdDVNBfCH97aSE0p6SR0w0V1ELAEWjqUpUpdFdFWl9kZO0S8FQAn60Zet5rqYTwHYVimzwITYF1mWIkT7FNdpNVr6JMYwUb+RB43/5fyq/8DWL6McmRBCiK5CVdXmGWBDM63kpmgJMNkJUnQXEUuAhWPHFKnLIrqqyqoKDlVCGzyc+m8AxrAZUKUQvhD7MX/+fIqKivD7/ZSXl7Nt2zZAW5ay+3L6pmUpCQkJLV4/a9Ys1q1bx+rVq1m/fj3/+c9/ZHmj6BGKQjPAchJNDP313ubH6xfdBAFftMISQgjRhVTWe6h2edEpMCQjgb4pcQDskhlgopuIWAIsHEtTZMtg0VXpCn9Gr6hUmvrAqPNAb8YarKO/UsrWij1rSQghhBAHo2kHyJN8n6OvXI9bn0CVaiXeuYOdn/4rusEJIYToEkpDu9j3slqwGPX0TWlaAikzwET30KYEmCxNEaJ1ksp/BqAkaQIYzJA1DoBxus3kV8kIihBCiPAqqm3AjJcpxf8FwHL87XyWPheAlOWPYa+tiWZ4QgghuoAat1arODnOBNA8A6xAEmCim2jTLpCyY4oQrZNjXwGAM/OI0AMTYdePjFc286LTE8XIhBBCdEfFtQ2MUHZi8dVBXDq68Zdx9tggZfe/QgaV5K35ntHHnBHtMIUQQnRiNfUtE2DNNcBqXASDKjqdErXYhAiHiO4CKUSP5Koi27tD+3fuZO2/OVoibLxuM+UOSYAJIYQIr+K6Bkbqdmp3+owFvQGzyURp/HAA3PmrohidEEKIrqA2NAMsKVZLgPVOjEGnQKMvSFW99GFE1ycJMCHCbecPAGwMZpPSq4/2WPYEVBRydeXoXGX4A8EoBiiEEKI7qff4sTf4GKHL1x7IHN38nDf9UABMlWs6PjAhhBBdSo2r5Qwwo17X/O+q0OwwIboySYAJEWb+HdrOpj8Hh9PbFqM9aLFBL20U/jBlizQgQgghwqYkVAD/UH2+9kDmqObn4nMPAyDDtbmjwxJCCNHF1P6uBtju/25KjgnRlUkCTIgw85VvAmCTbiDWmN/K7Cl9jwRggm4T5Y7GqMQmhBCi+6lyejDjZQBF2gO7JcCyhk3U/quWYq+rjkZ4QgghuoimJFdSrAE+uxW+vofkWKP2nFsSYKLrkwSYEGGmq8sHoCE+C0XZrVBk77EADFaKJAEmhBAibGrcXoYqBegJQlwaWHs3P2dLyaBUSQOgaP3SaIUohBCiC2hKgI0ufx+WPgmLH2Eay7XnpAaY6AYkASZEOAV8mFxlAARtOS2fS8oFIEuppEJ2ghRCCBEmNS4vI5oK4GeOAqXlLl1lsUMAcOav6OjQhBBCdCE1Li/p1DJ0/T+bHzu7+mnMeGUJpOgWJAEmRDjZC1EI0qgasST2aflcopYQ661UU2l3RSE4IYQQ3VGNy8sIpSkBNnqP5xtTRwBgKF/bgVEJIYToamrdPu4yvoTBVw+9x4C1D0neUq7Qf0S1JMBENyAJMCHCqTYfgEI1nd5JMS2fS8ggoBgwKgEaqos7PjYhhBDdkjYDLF+7s1v9ryYxfbVC+On1GzswKiGEEF1JIKiS1JDPSfplqIoeZj4BJ/wDgD8aFuF21kY5QiEOniTAhAin2l0AFKppZNgsLZ/T6XHHZGr/thd0cGBCCCG6K7uznsFKoXZnLwmwPkNDhfADxTT
     },
     "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-25T13:18:56.970354800Z",
     "start_time": "2023-07-25T13:18:55.981784100Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.000997304916381836\n",
      "0.015292882919311523\n",
      "0.0\n",
      "\n",
      "0.0009982585906982422\n",
      "0.028937816619873047\n",
      "0.0\n",
      "\n",
      "0.0\n",
      "0.0226593017578125\n",
      "0.0\n",
      "\n",
      "0.0009982585906982422\n",
      "0.015099287033081055\n",
      "0.0\n",
      "\n",
      "0.0011451244354248047\n",
      "0.01217198371887207\n",
      "0.0\n",
      "\n",
      "0.0\n",
      "0.01015615463256836\n",
      "0.0\n",
      "\n",
      "0.0\n",
      "0.020394086837768555\n",
      "0.0\n",
      "\n",
      "0.0\n",
      "0.022302627563476562\n",
      "0.0\n",
      "\n",
      "0.0009968280792236328\n",
      "0.008119821548461914\n",
      "0.0\n",
      "\n",
      "total time: 175.58670043945312 ms\n"
     ]
    }
   ],
   "source": [
    "# from opencv import moments\n",
    "start = time.time()\n",
    "shape = np.shape(cropOD)\n",
    "sigma = 0.4\n",
    "blurred = gaussian_filter(cropOD, sigma=sigma)\n",
    "\n",
    "thresh = np.zeros(shape)\n",
    "for i in range(0,shape[0]):\n",
    "    for j in range(0, shape[1]):\n",
    "        thresh[i,j] = np.where(blurred[i,j] < np.max(blurred[i,j])*0.5,0,1)\n",
    "\n",
    "center = calc_cen_bulk(thresh)\n",
    "\n",
    "BEC_width_guess = guess_BEC_width(thresh, center)\n",
    "\n",
    "X_guess_og = np.zeros((shape[0], shape[1], shape[3]))\n",
    "# Y_guess_og = np.zeros((shape[0], shape[1], shape[2]))\n",
    "\n",
    "for i in range(0, shape[0]):\n",
    "    for j in range(0, shape[1]):\n",
    "        X_guess_og[i,j] = np.sum(cropOD[i,j,round(center[i,j,1] - BEC_width_guess[i,j,1]/2) : round(center[i,j,1] + BEC_width_guess[i,j,1]/2) , :], 0) / len(cropOD[i,j,round(center[i,j,1] - BEC_width_guess[i,j,1]/2) : round(center[i,j,1] + BEC_width_guess[i,j,1]/2),0])\n",
    "\n",
    "        # Y_guess_og[i,j] = np.sum(cropOD[i,j, :, round(center[i,j,0] - BEC_width_guess[i,j,0]/2) : round(center[i,j,0] + BEC_width_guess[i,j,0]/2)], 1) / len(cropOD[i,j,0,round(center[i,j,0] - BEC_width_guess[i,j,0]/2) : round(center[i,j,0] + BEC_width_guess[i,j,0]/2)])\n",
    "\n",
    "#[nr images x, nr images y, X / Y, center / sigma]\n",
    "x = np.linspace(0,shape[3],150)\n",
    "y = np.linspace(0,shape[2], 150)\n",
    "\n",
    "popt = np.zeros((shape[0], shape[1], 6))\n",
    "\n",
    "p0 = np.ones((shape[0], shape[1], 6))\n",
    "\n",
    "max = np.zeros((shape[0], shape[1]))\n",
    "\n",
    "for i in range(0, shape[0]):\n",
    "    max[i] = np.ndarray.max(X_guess_og[i],axis=1)\n",
    "\n",
    "# Fitting x without math constr\n",
    "fitmodel = lmfit.Model(density_1d,independent_vars=['x'])\n",
    "\n",
    "result_x = []\n",
    "\n",
    "for i in range(0, shape[0]):\n",
    "    temp_res = []\n",
    "    for j in range(0, shape[1]):\n",
    "        t1 = time.time()\n",
    "        params = lmfit.Parameters()\n",
    "        params.add_many(\n",
    "            ('x0_bec', center[i,j,0], True,0, 200),\n",
    "            ('x0_th',center[i,j,0], True,0, 200),\n",
    "            ('amp_bec', 0.7 * max[i,j], True, 0, 1.3 * np.max(X_guess_og[i,j])),\n",
    "            ('amp_th', 0.3 * max[i,j], True, 0, 1.3 * np.max(X_guess_og[i,j])),\n",
    "            ('sigma_bec',BEC_width_guess[i,j,0], True, 0, 50),\n",
    "            ('sigma_th', popt[i,j,5], True, 0, 50)\n",
    "        )\n",
    "\n",
    "        t2 = time.time()\n",
    "        res = fitmodel.fit(X_guess_og[i,j], x=x, params=params)\n",
    "        t3 = time.time()\n",
    "        temp_res.append(res)\n",
    "        t4 = time.time()\n",
    "        print(t2 - t1)\n",
    "        print(t3 - t2)\n",
    "        print(t4 - t3)\n",
    "        print(\"\")\n",
    "    result_x.append(temp_res)\n",
    "stop = time.time()\n",
    "\n",
    "print(f'total time: {(stop-start)*1e3} ms')"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T13:19:37.207813100Z",
     "start_time": "2023-07-25T13:19:37.021761400Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "0.017198801040649414\n",
      "0.0\n",
      "\n",
      "0.0\n",
      "0.012012481689453125\n",
      "0.0\n",
      "\n",
      "0.0010266304016113281\n",
      "0.018287181854248047\n",
      "0.0\n",
      "\n",
      "0.0002391338348388672\n",
      "0.011229991912841797\n",
      "0.0\n",
      "\n",
      "0.00096893310546875\n",
      "0.016185998916625977\n",
      "0.0\n",
      "\n",
      "0.0011882781982421875\n",
      "0.011206388473510742\n",
      "0.0\n",
      "\n",
      "0.0\n",
      "0.026421308517456055\n",
      "0.0\n",
      "\n",
      "0.0\n",
      "0.007258415222167969\n",
      "0.0\n",
      "\n",
      "0.0020270347595214844\n",
      "0.07025146484375\n",
      "0.0\n",
      "\n",
      "total time: 206.9249153137207 ms\n"
     ]
    }
   ],
   "source": [
    "# from opencv import moments\n",
    "start = time.time()\n",
    "shape = np.shape(cropOD)\n",
    "sigma = 0.4\n",
    "blurred = gaussian_filter(cropOD, sigma=sigma)\n",
    "\n",
    "thresh = np.zeros(shape)\n",
    "for i in range(0,shape[0]):\n",
    "    for j in range(0, shape[1]):\n",
    "        thresh[i,j] = np.where(blurred[i,j] < np.max(blurred[i,j])*0.5,0,1)\n",
    "\n",
    "center = calc_cen_bulk(thresh)\n",
    "\n",
    "BEC_width_guess = guess_BEC_width(thresh, center)\n",
    "\n",
    "X_guess_og = np.zeros((shape[0], shape[1], shape[3]))\n",
    "# Y_guess_og = np.zeros((shape[0], shape[1], shape[2]))\n",
    "\n",
    "for i in range(0, shape[0]):\n",
    "    for j in range(0, shape[1]):\n",
    "        X_guess_og[i,j] = np.sum(cropOD[i,j,round(center[i,j,1] - BEC_width_guess[i,j,1]/2) : round(center[i,j,1] + BEC_width_guess[i,j,1]/2) , :], 0) / len(cropOD[i,j,round(center[i,j,1] - BEC_width_guess[i,j,1]/2) : round(center[i,j,1] + BEC_width_guess[i,j,1]/2),0])\n",
    "\n",
    "        # Y_guess_og[i,j] = np.sum(cropOD[i,j, :, round(center[i,j,0] - BEC_width_guess[i,j,0]/2) : round(center[i,j,0] + BEC_width_guess[i,j,0]/2)], 1) / len(cropOD[i,j,0,round(center[i,j,0] - BEC_width_guess[i,j,0]/2) : round(center[i,j,0] + BEC_width_guess[i,j,0]/2)])\n",
    "\n",
    "#[nr images x, nr images y, X / Y, center / sigma]\n",
    "x = np.linspace(0,shape[3],150)\n",
    "y = np.linspace(0,shape[2], 150)\n",
    "\n",
    "popt = np.zeros((shape[0], shape[1], 6))\n",
    "\n",
    "p0 = np.ones((shape[0], shape[1], 6))\n",
    "\n",
    "max = np.zeros((shape[0], shape[1]))\n",
    "\n",
    "for i in range(0, shape[0]):\n",
    "    max[i] = np.ndarray.max(X_guess_og[i],axis=1)\n",
    "\n",
    "# Fitting x without math constr\n",
    "fitmodel = lmfit.Model(density_1d,independent_vars=['x'])\n",
    "\n",
    "result_x = []\n",
    "\n",
    "for i in range(0, shape[0]):\n",
    "    temp_res = []\n",
    "    for j in range(0, shape[1]):\n",
    "        t1 = time.time()\n",
    "        params = lmfit.Parameters()\n",
    "        params.add_many(\n",
    "            ('x0_bec', center[i,j,0], True,0, 200),\n",
    "            ('x0_th',center[i,j,0], True,0, 200),\n",
    "            ('amp_bec', 0.7 * max[i,j], True, 0, 1.3 * np.max(X_guess_og[i,j])),\n",
    "            ('amp_th', 0.3 * max[i,j], True, 0, 1.3 * np.max(X_guess_og[i,j])),\n",
    "            ('deltax', 0, True, 0,50),\n",
    "            ('sigma_bec',BEC_width_guess[i,j,0], True, 0, 50)\n",
    "        )\n",
    "        params.add('sigma_th', popt[i,j,5], min=0, expr=f'sigma_bec + deltax')\n",
    "\n",
    "        t2 = time.time()\n",
    "        res = fitmodel.fit(X_guess_og[i,j], x=x, params=params)\n",
    "        t3 = time.time()\n",
    "        temp_res.append(res)\n",
    "        t4 = time.time()\n",
    "        print(t2 - t1)\n",
    "        print(t3 - t2)\n",
    "        print(t4 - t3)\n",
    "        print(\"\")\n",
    "    result_x.append(temp_res)\n",
    "stop = time.time()\n",
    "\n",
    "print(f'total time: {(stop-start)*1e3} ms')"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T13:19:43.156076600Z",
     "start_time": "2023-07-25T13:19:42.942009500Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 1500x1000 with 9 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABMAAAAMpCAYAAADipk9SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3hb1fnA8e/VsCTbkrxHYjvO3mSQsAk0EBJWAgHKLmEGCqWUtuxSSlugA1p+LZRQ0rDKTBhhlL2SEDJJQvZ0vLctybasde/vjys7MUmI7UiWx/t5Hj2JpKt7Xwmio/Oec96jaJqmIYQQQgghhBBCCCFEL2WIdQBCCCGEEEIIIYQQQkSTJMCEEEIIIYQQQgghRK8mCTAhhBBCCCGEEEII0atJAkwIIYQQQgghhBBC9GqSABNCCCGEEEIIIYQQvZokwIQQQgghhBBCCCFEryYJMCGEEEIIIYQQQgjRq5liHUBHqKpKaWkpdrsdRVFiHY4QQvR4mqbh8Xjo168fBoOMiYC0NUIIEWnS1rQl7YwQQkRWe9uZHpUAKy0tJTc3N9ZhCCFEr1NUVEROTk6sw+gWpK0RQojokLZGJ+2MEEJEx+HamR6VALPb7YD+phwOR4yjEUKIns/tdpObm9v6/SqkrRFCiEiTtqYtaWeEECKy2tvO9KgEWMsUYYfDIY2FEEJEkCzB2EfaGiGEiA5pa3TSzgghRHQcrp2RRfhCCCGEEEIIIYQQoleTBJgQQgghhBBCCCGE6NUkASaEEEIIIYQQQgghejVJgAkhhBBCCCGEEEKIXk0SYEIIIYQQQgghhBCiV5MEmBBCCCGEEEIIIYTo1SQBJkSU+YIh3lhbTH2TP9ahCCGE6OW+3F7F2sK6WIchhBCih/t6ZzUr99TGOgwhIkoSYEJE2YPvbOb219bz+Kc7Yh2KEEKIXqy6wcc1z67iqvkrCYbUWIcjRLd16623kp+fj6IorFu37qDHLFiwgPHjx7fe0tLSmD17NgAFBQUYjcY2z+/atasL34EQ0VXf5OeqBSuZs2AlvmAo1uEIETGmWAcgRG+2q6qBV1YVAfBdsSvG0QghhOjNCqobCakaHl+QXVWNDM+yxzokIbqlCy+8kDvuuIOTTjrpkMdcffXVXH311a33x4wZw+WXX9563263HzJ5JkRPt6nUTSCkEQiFKK7zMjg9MdYhCRERkgATIor++uE2QqoGwI7KBjRNQ1GUGEclhBCiNyqp97b+fXOZSxJgQhzClClTOnT8ihUrqKysZObMmZ26ns/nw+fztd53u92dOo8QXWVL2b7/R/fWNEoCTPQasgRSiCj5trCO/20sx6CAooDLG6C6QeqACfFDZFmKEJ1XXLdfAqxUOthCRMr8+fO58sorMZvNrY81NjYyefJkJk6cyIMPPkgodOhlYg8//DBOp7P1lpub2xVhC9Fpm9skwJpiGIkQkSUJMCGi5M1vSwA4b3x/BqTEA7Cj0hPLkITo9i688EKWLl3KgAEDDnnM1Vdfzbp161pvWVlZB12W0nIbPHhwV4QuRMztPwNskyTAhIiIxsZGXnnlFa699trWx7KzsykpKWHVqlV88sknLFmyhEcfffSQ57j77rtxuVytt6Kioq4IXYhO21K2r88iCTDRm0gCTIgoKQ13RCblpzAkQ1+GsrOyIZYhCdHtTZkyhZycnHYff6TLUkBfmuJ2u9vchOiJSvafAVbmRtO0GEYjRO/w+uuvM3r0aEaNGtX6mMViISMjA4CUlBSuueYalixZcshzWCwWHA5Hm5sQ3ZU/qLJzv0H7gprGGEYjRGR1KAEmS1OEaL/S+mYAsp1WhmTo6+Z3VEgCTIhIOtJlKSBLU0Tvsf8MsPqmAGWu5hhGI0TvMH/+/DazvwAqKysJBAKAPojyxhtvMGHChFiEJ0TE7apqIBDaN4BSKDPARC/SoQSYLE0Rov3K3eEEWJKVoeEEmMwAEyJyIrEsBWRpiugdNE1rnQFmt+h7HEkdMCEObu7cueTk5FBcXMz06dMZMmQIANdddx2LFy9uPW7btm2sW7eOiy++uM3rly5dyoQJExg3bhwTJ04kKyuLe++9t0vfgxDR0lIAPyfZBkBRXRPBkBrLkISImA7tAtnVO6YI0VM1B0LUNvrIpI7sBCP+zPAMMEmACREx7V2W8tJLL3HHHXcc8jwWiwWLxRL1eIWIprqmAN6APttxyvB03ttQxqZSN6ePyoxxZEJ0P/PmzTvo488880yb+8OHD8fjObB+6+zZs1tXuAjR27QkwH40PINXVxfhD6qUuZrJDdc0FqIni2oNsCNdmiJ1WURPVe5qZo7xQ1ZYb8Hx94GMfv8Cxik7qW7wUdcoO0EKEQmyLEWIfVpmf6XbLUzITeLHxs9J3/p8jKMSQgjR07QUwB/dz0FeOOkldcBEbxG1BFgklqZIXRbRU5W5mjnXuBwAJeTDWLqaX9jeB2BnlcwCE+JQZFmKEJ1TUq/XaOmfZOM44xb+bP43l1X/H+z8JMaRCSGE6Ck0TWudATYy20F+qp4Ak50gRW8RtQRYJHZMkbosoqeqqqlmnBLe4OFsPck7ga2AJoXwhfgB8+bNo7i4mGAwSEVFBTt37gT0ZSn7L6dvWZZit9vbvH727Nls3LiR9evXs2nTJv7xj3/I8kbRJxSHZ4DlJcUxcu3vWx9vePsOCAVjFZYQQogepKrBR02jH4MCw7PsDEhNAGCvzAATvUTUEmCRWJoiWwaLnspQ+DUmRaXa3A8mXAlGC061nnylnB2VB9aSEEIIIY5Eyw6Q5/rfw1i9hUajkxrNTqJnF7s/+L8YRyeEEKInKAvvYp/psGI1GxmQ2rIEUmaAid6hQwkwWZoiRPskVejLH0tTjgGTBfpPBGCyYRsF1TKCIoQQIrKK67zYaOaUEr2It23GA3yYoQ9Epq78K676uliGJ4QQogeobdJrFackxAG0zgArlASY6CU6tAuk7JgiRPvkuVYD4M46PvzA8VC4nMnKNp7z+GIYmRBCiN6opM7LWGUPcUEPJGZhOPoqLhgXpOLhF8ikhnXrv2T8KefFOkwhhBDdWG1D2wRYaw2w2kZUVcNgUGIWmxCRENVdIIXokxpryPOH638NnKL/macnwiYZtlHhlgSYEEKIyCqp9zLWsEe/0/9oMBixxFkoTRwDQNPetTGMTgghRE9QF54BlhyvJ8D6JdkwKNAcUKlukD6M6PkkASZEpBXoGztsVXNJzczRH8s9Bg2FQYZylMYKgiE1hgEKIYToTRp8QVzeAGNaEmD9xrc+508/CoC4yg0xiEwIIURPUtvYdgaY2Who/Xt1eHaYED2ZJMCEiLDg7q8AWK6Oop/Tpj9oS4JMfUfUo5Xt0oAIIYSImNJwAfxxxgL9gezxrc8l5h8NQGbjti6OSgghRE9T970aYPv/vSU5JkRPJgkwISIsUL4FgK2GIThs+8rsKeFlkMcYtlLhbo5JbEIIIXqfao+PeJrJp1R/IHtc63M5o44DIFcrxVVfE4vwhBBC9BCtM8BsBnj7ZnjvV6TEm/XnmiQBJno+SYAJEWFK/V4AmhJyUZT9CkX210fhhytFkgATQggRMbVNfkYpBRjQwJ4N9szW55xp2ZQr6QAUb14RqxCFEEL0AC0JsIllr8G3L8KqfzNdXao/JzXARC8gCTAhIikUwNJYBoCalNf2ueR8APor1VTKTpBCCCEipLbRv68A/n7LH1uUxQ8DwLNnVRdGJYQQoqepbfTTnyqGb3689bELaucRT7MsgRS9giTAhIgkVxEKKs2aGWtS/7bPhRNi/ZRqqlyNMQhOCCFEb1Tb6D9oAfwWzWljATBVfNeFUQkhhOhp6hr9PGh+FmPIC7nHQXI+jkA1t5jeokYSYKIXkASYEJFUpy9/LNbS6Zdsa/ucPZuQYiROCeGtLY1BcEIIIXqj2kY/Y5QC/c5BZoDFD5gIQHrD1q4LSgghRI8SUjXSmvdwmvFbNIMZZv4fTH8YgGuN79PsljqSoueTBJgQkVRXAEChlkGW09r2OYMRry1b/3t9YdfGJYQQotfyeFwMUUr0O/sVwG/Rb6S+CUtuqBh
     },
     "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",
    "        bval = result_x[i][j].best_values\n",
    "        ax[i,j].plot(x, X_guess_og[i,j])\n",
    "        ax[i,j].plot(x, density_1d(x, **result_x[i][j].best_values), label = lab)\n",
    "        ax[i,j].plot(x, thermal(x, bval['x0_th'], bval['amp_th'], bval['sigma_th']))\n",
    "\n",
    "\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T13:19:48.324787200Z",
     "start_time": "2023-07-25T13:19:46.834362100Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1. 1. 1. 1. 1. 1.]\n",
      "[0. 0. 0. 0. 0. 0.]\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGZCAYAAAC0UVoBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAADXjElEQVR4nO39f5SmZXXnje5HCptCCqqgBBrp5m7txh8D2qgoEZUykREzo6Imk8moGTLjUvO+jiczWWdGzpmzIuvNwLhWxvdNzPtOnGQOyXkxEw0Sl47BGE0KjPEHJLSgwYQmfUOp1WkLq7oLKDoUPOePour57O9de/uoLVUF+7tWrdrPc1/3dV/39eu5ru/e1969fr/ft0KhUCgUCoUNwFM2ugCFQqFQKBSevKiFSKFQKBQKhQ1DLUQKhUKhUChsGGohUigUCoVCYcNQC5FCoVAoFAobhlqIFAqFQqFQ2DDUQqRQKBQKhcKGoRYihUKhUCgUNgwjG12A74VHH33Uvv3tb9vY2Jj1er2NLk6hUCgUCoUh0O/3bXFx0c466yx7ylNi3mPTL0S+/e1v244dOza6GIVCoVAoFH4AzMzM2Nlnnx1e3/QLkbGxscekGTM72czuxtVFSf2PIC9DnpZ0r4W8BHlU0j0M+fiBOC7JFoa4xz1Hn3UP5DFJx/tuh3ySpHsF5O/G+Z2AMj3EC1+DfJbkfSrkP4P8chsO2k5sm9Mhf0vSHW/r4/Tge/NF/W6Yylcfq/iRhyXhIcgs93cG4gkv9rc8FLXnqRYDhThD+uHTIP8d20nrhy/1d5CfCVn7OD/zXbUfRm12d5Lu2QNRX/27rPR7IbOsfyM34X2PQ96PyNh6Pt7pdrbn3/l0dgJk9nmtV5bj6ZBZR6w7M//CfL+nS7rRQNb5gp2ZA/dZFuNGyNzMzUs69l8891RpNDee2MfP8elYrTqciEeiAcry7ZRrbIvzIGt9RfO6puPcxP7xKsgyL53wjIH80D5c+Ec+ncv725DPk3QfhjwF+ZOS7l9DjubGL8lnjid9d4Jj+tthKt+4rFedL77+2P8HzOwn8Tu+PnqbPdbMkSNH7JRTTjGzg7ayELlvcHFSVlhsm1n8aD/3+T7dnV/GB7lGTKCiOTZ0+bbc4gMbmwPthXKT/iCs4rPy+TlButPkM5+LRu9Jh+1HnfHO4HszX/ZvDsRdUv8HWRzOQDKQx5qBvIj8LF4xm93FzCE/V9LhPXagbWdk0hvBJLtsCXgfy8pJ+htyD+uL9XqxpAv64aj0jaXog/ahL6xfhqgfm5nZRyBfinvkR2iek+LrNJP1wd++Gb24tL6ctku0yG8lHX4o7A7kLWNwmW2LeSXthxz8LMOXJR2vcZLXcYu2oea5M06j+ULSbUe62ai+FKiHSdT/3Dd9Ms63c9EGyczPqajXntRrfxofODayDVwwz3V2Hey/2btHdcT21I74DFsfWlb2FS5K9Lfgr4LyNJKO8yjfPeuvRyC3SRn47pxr90g6viPfT+toNd0RM9tlhw8ftpNPPjksZRmrFgqFQqFQ2DDUQqRQKBQKhcKGYQupZuZsRTWzH1dV73TGQKRK4kxJtgDZsWl3+XSOlvp7yEoD7h6IE3juPOk9paVUpbBugaRMpDylDKTSHf2uStoPQv53A3ESXysl62g35qeUIMreQ3nSHob32C7qAEctf3ogTkA1MK9tEenc/17Soa+4Om8lHdspoiw174CCHBF6fXlI6nwC8jzapkN1a/9dBcqamUEts/0kYT9SCWX1yvI0ko7jmPew7yktHOn69aVUVfAYeqKGdfUF2ntkyqdbRplGUOfLHN9aVqgkRnAtVQNSjdf4S1TXLbHPZzY/qCPe72zDzKwf5af9ifMr5L6qdTlmsj7FzkfVBdUJjeTNMlHNpUZI7UCcRB5zX5B0/A3huH04+N7Mq0ioptExTJXLbouB8rnn6twWzVkZmI5567iNTBSytmVdvsQn2/FYXTx6xOxbp5RqplAoFAqFwubFpj81s4YTjl9hOdxCUHfkWL2NR4yDSR7ZDgwJL8Bq9LYjko4sCFexlFu5h5+5mtcdTmQseY9PNs/3ZX6yA5gACzKPMszhnUZldbyEdG43oPWPsvczYy+uzKcG4kFJNol6Zfmc4bDshJzhcAP5DJ/Ool2g9hvsHHroHzQeX9K8g13Ncsac3Bd8b2bz3GmhLvu6Y0p29RFc07D9MnawTTJkf2PbSB8YRT07w2b2PWE2xnBtkW2muzaeYIERZIeZ4259aiAvT0u6l+La0vrfd4Dd8DJ37sricjyxHtRgOblmwaWHyaK0uNDE97Mfjsk8QNb0APtetttHGXS6mIEhsTOezE6XsR+yr6khOJiKOZZVxyDHO8cgx5yOb/SB83HtDs371ZDZftpf+ZnPUoaFfY9lPTF4jpmNIr8l3p+8U2qsSqDdJ6WsM6uHLh5I7h+gGJFCoVAoFAobhlqIFAqFQqFQ2DBsIWPVA2Z2shhBZvQx6bjM4IZ0mtKmNKAi/aiqhjbJIwLpWhogChXWjwyMEm9dI6BUH5FrodEhDatEY8f8eMuiGEzRGG3oXpVRluRySdOjLToGmyzTNGRVuTRDlc7XCylj3t/KPZFDudhYy/tySNSKZGT7+ly2G/sX0u2QvGeisgrNzHrukxLX/o50NKA+TpKRsV1Cu29H+WZbuamBnHmrI9in5N3HUL7FLD92ejYA6O1RoaYjW0KNUtEf0mDZtQf60ZioaRYT/z1raJLnRM808/2fahWpVxplj+P7ORnfTj3H+ueYy/wvsXxa4ZG6b1rSBf5L2K11nnOqC757ZlxKtXrmA4r9S/ODKtC58HAW53JPpG5SBIawHQP7oB+GhvhHzGyyjFULhUKhUChsXmwhRuSvLT2ua2Z+ZwpmobML4eqUq041zOHqrwnukXRulc802XGrxPjMvRMMNkfEkIxHDB2Tk+yu3e6AOwPdST0juKblboN7dAfAeuXSXq3ZAjbI7c61XlkP3B1oWdF3JnFtTo87Mw8+iwbLavPN90iOZIb1oDsX9nMaAatRH/tvZHCm78cyZB4aAXdCVMcCx2ib5I2d2iiuLWX11QbXMmaU13RHNhpcU2N0XmMfQP1rOCznSRZ11JM2c6wW21n660hwSaMnOOaQdY7n8ji4mdk8+wDLkPln54P/Sq5dsv4tE9Kv5+GCfgJhN+YTlsh5n+VYaCXdzyAd+1Q21xKYvzrlZvk+Cvktkgf7CuqYnqXNxLt0BrZn5JG6lXv4vpy7lWkiUF8j8vvhCBeUYVTG99Jq2zxgZj9RjEihUCgUCoXNi61zfNdOtZWdieotCe6EsFrrK3OC405cyXUc/US7AykD9doL+J6r95NUl8vVMle3usMPjk8pOeScmMGGImvh46fW/35JV65cIYOJWdadKFka2hG8WtKxLujY6iuSjk5ysAPrfzJIo+DOVndC2JHfx9W8bjG52+N7JLZFbtcWOHcyM7/TZvl0R85nvTT4Xp57AvqNe6zokB27xPw0b7Tt0vRAnpjy6eYjtlF11+iY7mgqxlzqAI7XlG1hp0e5O8wod8rYLfZkt+jqBcwCm31GjslPgmGcQ1/rK4vFto5imJjZMpiOxcyxVTRnYawu6JzSQs52zcyPtgd6fJcMAup4XnbXdj6uZY68ADcns77kKLXj+TNWmLG9UB721058JtjHjIIF6ZDerCOUdVHjErGeOQdqhujnzrFeFP9FkTkiDGLcLGesGNq58+7L8j9HMSKFQqFQKBQ2DLUQKRQKhUKhsGHYQqqZ0cf+SNs1kobXjg9kM8eppsZxgWGaYpaGSIzBgPsXs6qmOkDLEMTrWNZ3IgWG/Mal3O7YJOm5RI1BWn45Me4lQ7ufvPW0Tzc2hQ+orw5lScqdbcEHyZFAlpX0bOdkdjuQT+L3mpDqGJaBdK8Ym9EYejveYVaydnQo2qkTk4PtxHtEhdPHZ6deY18RqtUZ+ybeJPvBWOjE+mEemcFydHQcNLXSwhNo2/lsbMIIchRGkEsZJQ4VYcdjLQ098X70rnu/ePakN08eye8cuWY/4hjW8c32TOZAqrNc/eH
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(p0[0,0])\n",
    "print(popt[0,0])\n",
    "\n",
    "data = cropOD[0,0]\n",
    "plt.pcolormesh(data,cmap='jet')\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T07:25:05.902906400Z",
     "start_time": "2023-07-25T07:25:05.680454200Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T13:10:47.207582500Z",
     "start_time": "2023-07-25T13:10:47.169500700Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. 0. 0. 0.]\n",
      "[ 0.  0. nan nan  0.  0.]\n",
      "[1. 1. 1. 1. 1. 1.]\n",
      "2.657634270170008\n"
     ]
    }
   ],
   "source": [
    "popt_0 = popt[0,0]\n",
    "print(popt_0)\n",
    "S = np.max(blurred[0,0])/(popt_0[2] + popt_0[3])\n",
    "popt_0[2] *= S\n",
    "popt_0[3] *= S\n",
    "print(popt_0)\n",
    "print(p0[0,0])\n",
    "print(np.max(blurred[0,0]))"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T13:12:07.279465300Z",
     "start_time": "2023-07-25T13:12:07.160965700Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "outputs": [],
   "source": [
    "import time"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T12:27:35.839077Z",
     "start_time": "2023-07-25T12:27:35.780727700Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Fitting without mathematical constraint"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.333087682723999\n",
      "0.3219437599182129\n",
      "0.3473355770111084\n",
      "0.3120737075805664\n",
      "0.3252856731414795\n",
      "0.33151793479919434\n",
      "1.4408035278320312\n",
      "0.42580461502075195\n",
      "0.37368178367614746\n",
      "fitting time = 0.4679482513003879 +- 0.3664738155004792\n"
     ]
    }
   ],
   "source": [
    "result = []\n",
    "times = []\n",
    "x = np.linspace(0,shape[3],150)\n",
    "y = np.linspace(0,shape[2], 150)\n",
    "\n",
    "for i in range(0,shape[0]):\n",
    "    temp_res_arr = []\n",
    "    for j in range(0,shape[1]):\n",
    "        data = cropOD[i,j]\n",
    "        fitModel = lmfit.Model(density_profile_BEC_2d, independent_vars=['x','y'])\n",
    "        #fitModel.set_param_hint('deltax', value=5)\n",
    "\n",
    "        bval = result_x[i][j].best_values\n",
    "        S = np.max(blurred[i,j])/(bval['amp_bec'] + bval['amp_th'])\n",
    "\n",
    "        params = lmfit.Parameters()\n",
    "        #print(bval['sigma_th'])\n",
    "        params.add_many(\n",
    "            ('amp_bec', S * bval['amp_bec'], True, 0, 1.3 * np.max(data)),\n",
    "            ('amp_th',S * bval['amp_th'], True, 0, 1.3 * np.max(data)),\n",
    "            ('x0_bec',center[i,j,0], True, 0, 150),\n",
    "            ('y0_bec',center[i,j,1], True, 0, 150),\n",
    "            ('x0_th',center[i,j,0], True, 0, 150),\n",
    "            ('y0_th',center[i,j,1], True, 0, 150),\n",
    "            ('sigmax_bec',bval['sigma_bec'], True, 0, 50),\n",
    "            ('sigmay_bec',BEC_width_guess[i,j,1], True, 0, 50),\n",
    "            ('sigmax_th',bval['sigma_th'], True, 0, 50),\n",
    "            ('sigmay_th',bval['sigma_th'],True, 0, 50)\n",
    "        )\n",
    "\n",
    "        X,Y = np.meshgrid(x, y)\n",
    "        X_1d = X.flatten()\n",
    "        Y_1d = Y.flatten()\n",
    "\n",
    "        data1d = data.flatten()\n",
    "        start = time.time()\n",
    "        res = fitModel.fit(data1d, x=X_1d, y=Y_1d, params=params)\n",
    "        stop = time.time()\n",
    "        print(stop-start)\n",
    "        times.append(stop-start)\n",
    "        temp_res_arr.append(res)\n",
    "    result.append(temp_res_arr)\n",
    "times = np.array(times)\n",
    "print(f\"fitting time = {np.mean(times)} +- {np.std(times, ddof=1)}\")\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-26T07:27:47.356522500Z",
     "start_time": "2023-07-26T07:27:43.110892300Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Fitting with mathematical constraint"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.6134850978851318\n",
      "1.2167115211486816\n",
      "0.908787727355957\n",
      "1.138059377670288\n",
      "1.0163965225219727\n",
      "1.137791395187378\n",
      "3.49874210357666\n",
      "1.6650748252868652\n",
      "1.1551790237426758\n",
      "fitting time = 1.48335862159729 +- 0.7966618144547525\n"
     ]
    }
   ],
   "source": [
    "result = []\n",
    "times = []\n",
    "x = np.linspace(0,shape[3],150)\n",
    "y = np.linspace(0,shape[2], 150)\n",
    "\n",
    "for i in range(0,shape[0]):\n",
    "    temp_res_arr = []\n",
    "    for j in range(0,shape[1]):\n",
    "        data = cropOD[i,j]\n",
    "        fitModel = lmfit.Model(density_profile_BEC_2d, independent_vars=['x','y'])\n",
    "        #fitModel.set_param_hint('deltax', value=5)\n",
    "\n",
    "        bval = result_x[i][j].best_values\n",
    "        S = np.max(blurred[i,j])/(bval['amp_bec'] + bval['amp_th'])\n",
    "\n",
    "        params = lmfit.Parameters()\n",
    "        #print(bval['sigma_th'])\n",
    "\n",
    "        params.add_many(\n",
    "            ('amp_bec', S * bval['amp_bec'], True, 0, 1.3 * np.max(data)),\n",
    "            ('amp_th',S * bval['amp_th'], True, 0, 1.3 * np.max(data)),\n",
    "            ('x0_bec',center[i,j,0], True, 0, 150),\n",
    "            ('y0_bec',center[i,j,1], True, 0, 150),\n",
    "            ('x0_th',center[i,j,0], True, 0, 150),\n",
    "            ('y0_th',center[i,j,1], True, 0, 150),\n",
    "            ('sigmax_bec',bval['sigma_bec'], True, 0, 50),\n",
    "            ('sigmay_bec',BEC_width_guess[i,j,1], True, 0, 50),\n",
    "            ('D_sigX', 1.93*bval['sigma_th'] - 1.22*bval['sigma_bec'], True, 0.1, 50),\n",
    "            ('D_sig_th', 0, True, -10, 10)\n",
    "        )\n",
    "\n",
    "        params.add('sigmax_th',bval['sigma_th'], min=0, expr=f'0.632*sigmax_bec + 0.518*D_sigX')\n",
    "        params.add('sigmay_th',bval['sigma_th'], min=0, expr=f'sigmax_th + D_sig_th')\n",
    "\n",
    "        X,Y = np.meshgrid(x, y)\n",
    "        X_1d = X.flatten()\n",
    "        Y_1d = Y.flatten()\n",
    "\n",
    "        data1d = data.flatten()\n",
    "        start = time.time()\n",
    "        res = fitModel.fit(data1d, x=X_1d, y=Y_1d, params=params)\n",
    "        stop = time.time()\n",
    "        print(stop-start)\n",
    "        times.append(stop-start)\n",
    "        temp_res_arr.append(res)\n",
    "    result.append(temp_res_arr)\n",
    "times = np.array(times)\n",
    "print(f\"fitting time = {np.mean(times)} +- {np.std(times, ddof=1)}\")\n"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T14:21:05.284099300Z",
     "start_time": "2023-07-25T14:20:51.892911Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Fitting class"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 800x800 with 6 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqgAAAKvCAYAAACmvT3kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOyde5xVU//H34upDE2TJqWbookG48lIYWJKoWJQulBSlNLjCeFHCfscHsLjHhKlSCjdFCVCpVColIznqTTpQuliSkY1rN8fa62z197nnC4kU9anV685Z+21122fvdfan/X9fr5CSilxcHBwcHBwcHBwKCU45K9ugIODg4ODg4ODg4MNt0B1cHBwcHBwcHAoVXALVAcHBwcHBwcHh1IFt0B1cHBwcHBwcHAoVXALVAcHBwcHBwcHh1IFt0B1cHBwcHBwcHAoVXALVAcHBwcHBwcHh1IFt0B1cHBwcHBwcHAoVXALVAcHBwcHBwcHh1IFt0B1cHBw+AsghEAIQdOmTf/qpjg4OPzFiEQisWfCjBkz/urmlAqk/NUNcPAhhEiYXrZsWSpUqEB6ejq1a9cmJyeHxo0bc+GFF5Kamrpf2zhixAgKCwsBdUM5OPydkeyeTYYbb7yRxx9/fI/yTpw4kYULFwJw0003UbFixb1rnIPDHsLNPfsGkUiEaDS6V+ds3rx5j+7tH3/8MfbsaNCgAZdeeuneN/AAg1ugHgDYsWMHGzZsYMOGDSxfvpz3338fgIoVK9K1a1ei0Sjp6en7pS0jRoxg5syZQOl9SDg4HAyYOHEiL774IgDdunVzC1SH/Q4395Qe/Pjjj7HFb9euXd0C1eGvw4QJE2KfpZQUFRWxefNmFi5cyKxZsygsLOTHH3/kiSeeYNy4cbz66qs0adLkL2yxg8PfG/Y9mwx169aNfZZS/pnNcXD4XXBzz75Bx44dufzyy3eb74gjjgDUovvvuPDeFdwCtZRiV29HUkqmTp3KTTfdxNKlS1m9ejUXXXQRc+bM4aSTTtp/jXRwcIjh78BoOBz8cHPPvkH9+vXdM+EPwjlJHYAQQtC6dWs+++yz2JtrUVER7du357fffvuLW+fg4ODgcDDCzT0O+xNugXoAo0KFCowZMyZmm1ZQUMDo0aMT5i0uLmbChAlcf/31NG7cmIyMDMqUKUN6ejonnXQSvXv35osvvkhaV9OmTRFCxGyAwPdCtv+HtyhKSkqYNm0at9xyC02aNKFKlSqULVuWtLQ0jj/+eLp168asWbP+8Fg4OBxoSObF361bN4QQMftTgGOPPTbuXuvWrdv+bbCDg4abe/Y9knnxFxYWIoTg2GOPjaW9+OKLCcfAOJEdLHBb/Ac4qlWrRs+ePXnooYcAeOGFF7jiiivi8p144okJf7xbtmzhq6++4quvvuLZZ5+lf//+3H///fusfeedd15CyYydO3eydOlSli5dyosvvkjXrl157rnnKFu27D6r28HBwcHhz4Gbexz+bLgF6kGATp06xR4SH330ETt37qRMmTKBPMXFxVSqVInzzjuPU089lRo1alCmTBnWrFnD/PnzGTNmDDt37mTgwIFUqVKFm266KXD+v//9bzZs2MCdd97JkiVLgMROIfXr14+rt3z58jRv3pzTTjuNOnXqcNhhh/Hdd9+xZMkSRo0axbZt23jxxRepWLHiHkvwODgcrLjhhhu49NJLefLJJ/nggw8AGDJkCFWqVAnkO+aYY/6K5jk4xODmnj8fVapUYcKECaxfv55evXoB0KxZM2644YaEeQ8qSIdSAyD2f29QUlIijzjiiNi5CxcujMszdepUuXPnzqRlFBYWyvr160tApqWlyS1btiTMl5eXt1dtnD59uvz555+THt+wYYNs0qSJBOQhhxwiv/nmmz0q18GhNOD33rP2uXl5eQmPd+3aNZZnxYoVf6yhDg67gJt79s3c43lerI2e5/3ucz/44IO44ytWrIgd79q16x9u64EAZ4N6EODQQw+lZs2ase8//PBDXJ6WLVuSkpKcMK9duzbPPPMMAFu3buWNN97YJ21r3rz5LgWdMzIyYrZ2v/32G6NGjdon9To47G8ksglzNqMOBzPc3JMc0Wh0t8+Eg81mdF/DLVAPEhx55JGxzxs3bvxdZZx11lmxz3Pnzv3DbdpTHHfccRx99NH7vV4HBwcHhz8GN/c4/FlwNqgHCWyJj2Rh69avX89LL73EO++8w1dffcXmzZv5+eefE+ZdvXr1Pmvbli1bGDVqFFOmTGHx4sVs2LCBbdu2/en1OjjsT+xOqN/ZjDocjHBzT2LsiVD/QWczuo/hFqgHCX788cfY50qVKsUdHz16NL169aKoqGiPytuyZcs+adcHH3xAp06d+P777/drvQ4O+xtOlNvh74iDee555513ki6kYdf3vBPq/+NwC9SDAL/++mvg7e+oo44KHJ81axadOnWKvenm5OTQokUL6tatS3p6OuXKlYvlbdOmTazMP4qlS5dy4YUXUlxcDMAJJ5xAq1atqFevHpUqVeKwww6L5e3Zsyc//PDDPqnXwcHBweHPx8E+9/Ts2ZOVK1cmrUe6cMV/KtwC9SDA4sWLY295RxxxRFzIuUgkEntAPPfcc1x77bUJy0m29fF7MXDgwNgDYsCAAdx7771Jt4CStcnBwcHBoXTCzT0OfybcAvUgwCuvvBL7fNZZZwU8Jnfs2MGHH34IQMOGDXd5M+7qTfH3YPr06YCys7nnnnuSPiC2bt3Kpk2b9mndDg4ODg5/Lg72ucd52f+1cF78Bzi+++47nn/++dj37t27B45v3LiRkpISAOrWrbvLsqZNm7bb+g45xP/J7G57Y926dYAK02ifF8b06dNdHGcHhxD25l5zcNjfcHPP/sXf8XngFqgHMLZu3UqHDh1iRupZWVm0b98+kOfwww+PfV6+fPkuy3rsscd2W2f58uVjn3e3LWPq/uabb5LeUL/++us+DW/n4HCwYG/uNQeH/Qk39+x//B2fB26BegBCSsnUqVNp2LAhs2fPBqBChQq8/vrrcW+L6enp1KtXD4DPPvssoRTOTz/9RPv27Vm1atVu6z722GNjn+fPn7/LvKeffjqgxJsThZHbuXMn1157LZ999tlu63Vw+Lthb+41B4f9ATf3/HWoVKkS6enpACxcuPBvwaI6G9RSiokTJ8Y+SyljtjILFy5k1qxZrFixIna8Zs2avPrqq3EG6gZ9+vSJxe1t164dnTt3pkmTJqSlpfHll18yYsQI1q5dy1VXXcVLL720y3Y1b96cJ598ElBbOn379qV27doceuihAGRmZpKZmRmr99133wXg5ptvZsaMGVxwwQVkZGSwdOlSXnrpJZYuXUqzZs1YunSp00B1cLDQvHnz2OfbbruNH374gRNOOCFm51ejRg2ys7P/quY5HKRwc0/pxbnnnsuECRNYvnw5HTt2pG3btlSsWDF2PC8vb5fRsw44/DURVh0SASse8p78r1ixorzhhhvk5s2bd1nub7/9Jjt37rzLsi655BL5888/7zY+eElJSSx+caL/4fjD/fv332W9ubm5cv369bJ27doSkLVr194nY+ngsD9g/5Z/77nJ7jUppbziiiuS3jt/l3jcDn8+3Nyzb+Yez/OStmdvzv3ggw8S5lmwYIFMTU1N2qcVK1b84T6UJjgG9QBAmTJlqFChAhUqVKBOnTrk5OTQuHFjLrrooj16WxJC8PLLL3PhhRfy/PPPs2DBAn7++WeqVKlCgwYN6NKlCx06dNijthx66KG8++67PPHEE7zxxht8/fXXbNmyJal23f33388555zDU089xdy5cykqKqJy5cpkZWVxxRVX0K1bt13GaXZw+Dtj5MiRnHPOOYwePZovv/ySH3/8MeZ44uDwZ8PNPaULDRo04PPPP+fRRx/lww8/ZNWqVbsMJHCgQ0j5NzBkcHBwcHBwcHBwOGDgnKQcHBwcHBwcHBxKFdwC1cHBwcHBwcHBoVTBLVAdHBwcHBwcHBxKFdwC1cHBwcHBwcHBoVTBLVAdHBwcHBwcHBxKFdwC1cHBwcHBwcHBoVSh1IuA/fbbb6xdu5a
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 800x800 with 6 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqgAAAKvCAYAAACmvT3kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOydf5xN1frH36sGTZkZje5IaKaMolJIJApRoSglhVu4un7kKtJXfqR9TkK5SVKhm4uKSkpxS4kQpVzJTVEhIyRTZEY1YWp9/1hrnb32nnMGJQ2tz7zmdc5ee+31a5+999qf9TyfR0gpJQ4ODg4ODg4ODg4lBMf80Q1wcHBwcHBwcHBwsOEmqA4ODg4ODg4ODiUKboLq4ODg4ODg4OBQouAmqA4ODg4ODg4ODiUKboLq4ODg4ODg4OBQouAmqA4ODg4ODg4ODiUKboLq4ODg4ODg4OBQouAmqA4ODg4ODg4ODiUKboLq4ODg4ODg4OBQouAmqA4ODg5/AIQQCCFo0qTJH90UBweHPxiRSCR2T1i0aNEf3ZwSgaQ/ugEOPoQQcdNLly5NamoqaWlpZGZmUqdOHerXr8+VV15JcnLyYW3jlClTyMnJAdQF5eDwZ0aiazYRbr/9dh5++OEDyvvyyy+zatUqAPr27Uu5cuUOrnEODgcI9+w5NIhEIkSj0YM65rvvvjuga3vXrl2xe0etWrW45pprDr6BRxjcBPUIwN69e/n222/59ttv2bBhA2+99RYA5cqVo3PnzkSjUdLS0g5LW6ZMmcLixYuBknuTcHA4GvDyyy8zdepUALp06eImqA6HHe7ZU3Kwa9eu2OS3c+fOboLq8Mdh1qxZse9SSvLy8vjuu+9YtWoVb7/9Njk5OezatYuxY8fy4osv8uyzz9KoUaM/sMUODn9u2NdsIlStWjX2XUr5ezbHweFXwT17Dg1uuOEGbrzxxv3mO+GEEwA16f4zTryLg5ugllAU93YkpWTu3Ln07duXdevWsWXLFq666ireeecdzj777MPXSAcHhxj+DIyGw9EP9+w5NKhevbq7J/xGOCepIxBCCFq1asWKFStib655eXlcf/31/PLLL39w6xwcHBwcjka4Z4/D4YSboB7BSE1NZcaMGTHbtLVr1/L888/HzVtQUMCsWbPo3bs39evXp3z58pQqVYq0tDTOPvtsevXqxf/+97+EdTVp0gQhRMwGCHwvZPs/vERRWFjIG2+8Qf/+/WnUqBEZGRmULl2alJQUzjjjDLp06cLbb7/9m8fCweFIQyIv/i5duiCEiNmfApx22mlFrrUuXboc3gY7OGi4Z8+hRyIv/pycHIQQnHbaabG0qVOnxh0D40R2tMAt8R/hqFixIt27d2fUqFEA/Pvf/6ZDhw5F8p111llxf7z5+fmsWbOGNWvWMGHCBAYNGsSIESMOWfsuu+yyuJIZ+/btY926daxbt46pU6fSuXNnnnjiCUqXLn3I6nZwcHBw+H3gnj0OvzfcBPUoQMeOHWM3iXfffZd9+/ZRqlSpQJ6CggLS09O57LLLqF27NpUqVaJUqVJs3bqVlStXMmPGDPbt28fIkSPJyMigb9++gePvu+8+vv32W+6++24++eQTIL5TSPXq1YvUW7ZsWZo1a8b5559PVlYWxx13HNu2beOTTz5h2rRp/PDDD0ydOpVy5codsASPg8PRittuu41rrrmGRx55hIULFwIwceJEMjIyAvlOPfXUP6J5Dg4xuGfP74+MjAxmzZpFbm4uPXr0AKBp06bcdtttcfMeVZAOJQZA7P9gUFhYKE844YTYsatWrSqSZ+7cuXLfvn0Jy8jJyZHVq1eXgExJSZH5+flx8zVu3Pig2jh//nz5448/Jtz/7bffykaNGklAHnPMMfKLL744oHIdHEoCfu01ax/buHHjuPs7d+4cy7Nx48bf1lAHh2Lgnj2H5tnjeV6sjZ7n/epjFy5cWGT/xo0bY/s7d+78m9t6JMDZoB4FOPbYY6lcuXJs+5tvvimSp0WLFiQlJSbMMzMzefzxxwHYvXs3r7zyyiFpW7NmzYoVdC5fvnzM1u6XX35h2rRph6ReB4fDjXg2Yc5m1OFohnv2JEY0Gt3vPeFosxk91HAT1KMEJ554Yuz7jh07flUZF110Uez7+++//5vbdKA4/fTTOfnkkw97vQ4ODg4Ovw3u2ePwe8HZoB4lsCU+EoWty83N5amnnmLevHmsWbOG7777jh9//DFu3i1bthyytuXn5zNt2jRee+01Vq9ezbfffssPP/zwu9fr4HA4sT+hfmcz6nA0wj174uNAhPqPOpvRQww3QT1KsGvXrtj39PT0Ivuff/55evToQV5e3gGVl5+ff0jatXDhQjp27MjXX399WOt1cDjccKLcDn9GHM3Pnnnz5iWcSEPx17wT6v/tcBPUowA///xz4O3vL3/5S2D/22+/TceOHWNvunXq1KF58+ZUrVqVtLQ0ypQpE8vbtm3bWJm/FevWrePKK6+koKAAgDPPPJOWLVtSrVo10tPTOe6442J5u3fvzjfffHNI6nVwcHBw+P1xtD97unfvzqZNmxLWI1244t8VboJ6FGD16tWxt7wTTjihSMi5SCQSu0E88cQT/P3vf49bTqKlj1+LkSNHxm4QQ4YMYdiwYQmXgBK1ycHBwcGhZMI9exx+T7gJ6lGA6dOnx75fdNFFAY/JvXv3smTJEgDq1q1b7MVY3Jvir8H8+fMBZWdz7733JrxB7N69m507dx7Suh0cHBwcfl8c7c8e52X/x8J58R/h2LZtG//6179i2926dQvs37FjB4WFhQBUrVq12LLeeOON/dZ3zDH+T2Z/yxvbt28HVJhG+7gw5s+f7+I4OziEcDDXmoPD4YZ79hxe/BnvB26CegRj9+7dtG/fPmakXqNGDa6//vpAnuOPPz72fcOGDcWWNWbMmP3WWbZs2dj3/S3LmLq/+OKLhBfUzz//fEjD2zk4HC04mGvNweFwwj17Dj/+jPcDN0E9AiGlZO7cudStW5elS5cCkJqaygsvvFDkbTEtLY1q1aoBsGLFirhSON9//z3XX389mzdv3m/dp512Wuz7ypUri817wQUXAEq8OV4YuX379vH3v/+dFStW7LdeB4c/Gw7mWnNwOBxwz54/Dunp6aSlpQGwatWqPwWL6mxQSyhefvnl2HcpZcxWZtWqVbz99tts3Lgxtr9y5co8++yzRQzUDfr06ROL29uuXTs6depEo0aNSElJ4eOPP2bKlCl89dVX3HzzzTz11FPFtqtZs2Y88sgjgFrS6devH5mZmRx77LEAZGdnk52dHav3zTffBOCOO+5g0aJFXHHFFZQvX55169bx1FNPsW7dOpo2bcq6deucBqqDg4VmzZrFvg8YMIBvvvmGM888M2bnV6lSJWrWrPlHNc/hKIV79pRcXHrppcyaNYsNGzZwww03cO2111KuXLnY/saNGxcbPeuIwx8TYdUhHrDiIR/If7ly5eRtt90mv/vuu2LL/eWXX2SnTp2KLevqq6+WP/74437jgxcWFsbiF8f7D8cfHjRoULH1NmzYUObm5srMzEwJyMzMzEMylg4OhwP2b/nXHpvoWpNSyg4dOiS8dv4s8bgdfn+4Z8+hefZ4npewPQdz7MKFC+Pm+fDDD2VycnLCPm3cuPE396EkwTGoRwBKlSpFamoqqampZGVlUadOHerXr89VV111QG9LQgieeeYZrrzySv71r3/x4Ycf8uOPP5KRkUGtWrW46aabaN++/QG15dhjj+XNN99k7NixvPLKK3z66afk5+cn1K4bMWIEl1xyCY8++ijvv/8+eXl5nHTSSdSoUYMOHTrQpUuXYuM0Ozj8mfH0009zySWX8Pzzz/Pxxx+za9eumOOJg8PvDffsKVmoVasWH3zwAQ899BBLlixh8+bNxQYSONIhpPwTGDI4ODg4ODg4ODgcMXBOUg4ODg4ODg4ODiUKboLq4ODg4ODg4OBQouAmqA4ODg4ODg4ODiUKboLq4ODg4ODg4OBQouAmqA4ODg4ODg4ODiUKboLq4ODg4ODg4OBQolDiRcB++eUXvvrqK1J
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 800x800 with 6 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqgAAAKvCAYAAACmvT3kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOyde5xX0/7/n5tJBjOTSkm3waTGEUlKpqPSUCHkrqScOhJyJD/CYX8+HOo4cr/l6LgktyNFXxJDhZBLRY6iYlLpomKmw6DR+v2x1vuz1t6fz0xFJ1PWq0eP+XzWXnvd9mfvvdZrvd+vd6CUUnh4eHh4eHh4eHjUEOz0WzfAw8PDw8PDw8PDw4WfoHp4eHh4eHh4eNQo+Amqh4eHh4eHh4dHjYKfoHp4eHh4eHh4eNQo+Amqh4eHh4eHh4dHjYKfoHp4eHh4eHh4eNQo+Amqh4eHh4eHh4dHjYKfoHp4eHh4eHh4eNQo+Amqh4eHh4eHh4dHjYKfoHp4eHj8BgiCgCAI6NKly2/dFA8Pj98YiUQi9UyYPn36b92cGoGs37oBHhZBEGRM32WXXcjNzSUvL4/mzZvTtm1bOnTowPHHH092dvY2bePDDz9MaWkpoG8oD4/fM6q6Z6vCX/7yF26//fbNyjtp0iTmzp0LwKWXXkqdOnW2rHEeHpsJ/+7ZOkgkEiSTyS0655tvvtmse/vbb79NPTvatGnDySefvOUN3M7gJ6jbAX766SfWrFnDmjVrWLx4Ma+99hoAderUoX///iSTSfLy8rZJWx5++GFmzJgB1NyHhIfHjoBJkybxyCOPADBgwAA/QfXY5vDvnpqDb7/9NjX57d+/v5+gevx2mDhxYuqzUoqysjK++eYb5s6dy+uvv05paSnffvstd9xxBxMmTOCJJ56gU6dOv2GLPTx+33Dv2aqw//77pz4rpf6XzfHw+EXw756tgzPPPJOzzjprk/l23313QE+6f48T7+rgJ6g1FNWtjpRSTJkyhUsvvZSFCxeybNkyTjjhBGbOnMkf/vCHbddIDw+PFH4PjIbHjg//7tk6aNWqlX8m/Ep4J6ntEEEQcNxxx/H++++nVq5lZWWcfvrpbNy48TdunYeHh4fHjgj/7vHYlvAT1O0Yubm5PP300ynbtPnz5/PUU09lzFtRUcHEiRO56KKL6NChA/Xq1aNWrVrk5eXxhz/8gSFDhvDhhx9WWVeXLl0IgiBlAwTWC9n9H9+iqKysZOrUqQwfPpxOnTrRoEEDdtllF3JycjjggAMYMGAAr7/++q8eCw+P7Q1VefEPGDCAIAhS9qcA++67b9q9NmDAgG3bYA8PA//u2fqoyou/tLSUIAjYd999U2mPPPJIxjEQJ7IdBX6LfztHo0aNOP/887n55psB+Ne//sXZZ5+dlu/AAw/M+OMtLy/nk08+4ZNPPuH+++/nqquu4qabbtpq7TvmmGMySmZs2LCBhQsXsnDhQh555BH69+/PAw88wC677LLV6vbw8PDw+N/Av3s8/tfwE9QdAH369Ek9JN566y02bNhArVq1InkqKiqoW7cuxxxzDIceeiiNGzemVq1aLF++nNmzZ/P000+zYcMGRo4cSYMGDbj00ksj5//tb39jzZo1/PWvf+U///kPkNkppFWrVmn17rHHHnTr1o3DDjuM/Px8dt11V1asWMF//vMfxo8fz3fffccjjzxCnTp1NluCx8NjR8Ull1zCySefzJ133sm0adMAGDNmDA0aNIjka9as2W/RPA+PFPy753+PBg0aMHHiRFavXs3gwYMB6Nq1K5dccknGvDsUlEeNAZD6vyWorKxUu+++e+rcuXPnpuWZMmWK2rBhQ5VllJaWqlatWilA5eTkqPLy8oz5OnfuvEVtLCkpUd9//32Vx9esWaM6deqkALXTTjupzz//fLPK9fCoCfil96x7bufOnTMe79+/fyrPF1988esa6uFRDfy7Z+u8e8IwTLUxDMNffO60adPSjn/xxRep4/379//Vbd0e4G1QdwDsvPPONGnSJPX966+/TsvTo0cPsrKqJsybN2/OvffeC8D69et57rnntkrbunXrVq2gc7169VK2dhs3bmT8+PFbpV4Pj22NTDZh3mbUY0eGf/dUjWQyuclnwo5mM7q14SeoOwj23HPP1Oe1a9f+ojKOPPLI1OdZs2b96jZtLvbbbz/23nvvbV6vh4eHh8evg3/3ePyv4G1QdxC4Eh9Vha1bvXo1jz76KC+//DKffPIJ33zzDd9//33GvMuWLdtqbSsvL2f8+PG8+OKLzJs3jzVr1vDdd9/9z+v18NiW2JRQv7cZ9dgR4d89mbE5Qv07nM3oVoafoO4g+Pbbb1Of69atm3b8qaeeYvDgwZSVlW1WeeXl5VulXdOmTaNPnz6sXLlym9br4bGt4UW5PX6P2JHfPS+//HKVE2mo/p73Qv2/Hn6CugPg559/jqz+9tprr8jx119/nT59+qRWum3btqW4uJj999+fvLw8ateuncrbu3fvVJm/FgsXLuT444+noqICgJYtW9KzZ09atGhB3bp12XXXXVN5zz//fL7++uutUq+Hh4eHx/8eO/q75/zzz2fJkiVV1qN8uOL/KfwEdQfAvHnzUqu83XffPS3kXCKRSD0gHnjgAf785z9nLKeqrY9fipEjR6YeENdccw033HBDlVtAVbXJw8PDw6Nmwr97PP6X8BPUHQCPP/546vORRx4Z8Zj86aefeOONNwBo165dtTdjdSvFX4KSkhJA29lcf/31VT4g1q9fz7p167Zq3R4eHh4e/1vs6O8e72X/28J78W/nWLFiBf/85z9T3wcOHBg5vnbtWiorKwHYf//9qy1r6tSpm6xvp53sT2ZT2xurVq0CdJhG97w4SkpKfBxnD48YtuRe8/DY1vDvnm2L3+PzwE9Qt2OsX7+eM844I2WkXlhYyOmnnx7Js9tuu6U+L168uNqybrvttk3Wuccee6Q+b2pbRur+/PPPq7yhfv75560a3s7DY0fBltxrHh7bEv7ds+3xe3we+AnqdgilFFOmTKFdu3a8+eabAOTm5vLvf/87bbWYl5dHixYtAHj//fczSuH897//5fTTT2fp0qWbrHvfffdNfZ49e3a1eQ8//HBAizdnCiO3YcMG/vznP/P+++9vsl4Pj98btuRe8/DYFvDvnt8OdevWJS8vD4C5c+f+LlhUb4NaQzFp0qTUZ6VUylZm7ty5vP7663zxxRep402aNOGJJ55IM1AXDB06NBW397TTTqNv37506tSJnJwcPv74Yx5++GG++uorzj33XB599NFq29WtWzfuvPNOQG/pDBs2jObNm7PzzjsDUFBQQEFBQareV155BYDLLruM6dOn0717d+rVq8fChQt59NFHWbhwIV27dmXhwoVeA9XDw0G3bt1Sn6+44gq+/vprWrZsmbLza9y4Ma1bt/6tmuexg8K/e2oujj76aCZOnMjixYs588wzOeWUU6hTp07qeOfOnauNnrXd4beJsOqRCTjxkDfnf506ddQll1yivvnmm2rL3bhxo+rbt2+1ZZ100knq+++/32R88MrKylT84kz/4/GHr7rqqmrrLSoqUqtXr1bNmzdXgGrevPlWGUsPj20B97f8S8+t6l5TSqmzzz67ynvn9xKP2+N/D//u2TrvnjAMq2zPlpw7bdq0jHnmzJmjsrOzq+zTF1988av7UJPgGdTtALVq1SI3N5fc3Fzy8/Np27YtHTp04IQTTtis1VIQBDz22GMcf/zx/POf/2TOnDl8//33NGjQgDZt2tCvXz/OOOOMzWrLzjvvzCuvvMIdd9zBc889x4IFCygvL69Su+6mm27iqKOO4u6772bWrFmUlZVRv359CgsLOfvssxkwYEC1cZo9PH7PGDduHEcddRRPPfUUH3/8Md9++23K8cTD438N/+6pWWjTpg0ffPABt956K2+88QZLly6tNpDA9o5Aqd+BIYOHh4eHh4eHh8d2A+8k5eHh4eHh4eHhUaPgJ6geHh4eHh4eHh41Cn6C6uHh4eHh4eHhUaPgJ6geHh4eHh4eHh41Cn6C6uHh4eHh4eH
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 800x800 with 6 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqgAAAKvCAYAAACmvT3kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOydfZxV0/7H3ys9GEyTylQqU5ooTLcHChOFQUVIKdVNkVtCkVzEZZ/Dj3ARQg+3VJLoSiWU1FUUihKlUKkUPehpCiOl9ftjrXX22vvsMxVhGuszr3mds9deez3ts/de+7O+389XSCklDg4ODg4ODg4ODkUEJf7sBjg4ODg4ODg4ODjYcBNUBwcHBwcHBweHIgU3QXVwcHBwcHBwcChScBNUBwcHBwcHBweHIgU3QXVwcHBwcHBwcChScBNUBwcHBwcHBweHIgU3QXVwcHBwcHBwcChScBNUBwcHBwcHBweHIgU3QXVwcHBwcHBwcChScBNUBwcHhz8BQgiEEDRv3vzPboqDg8OfjFgslrgnzJo1689uTpFAyT+7AQ4+hBCR6aVLl6Zs2bJkZGSQlZVFw4YNadKkCRdddBFpaWl/aBtHjRrF6tWrAXVBOTj8lZHqmk2Fm266iccff3y/8k6aNIlFixYBcPPNN1OuXLkDa5yDw37CPXsODmKxGPF4/ICO2bZt235d29u3b0/cO+rXr89ll1124A08xOAmqIcAfv75ZzZv3szmzZtZuXIl//vf/wAoV64cXbt2JR6Pk5GR8Ye0ZdSoUcyePRsoujcJB4figEmTJjF69GgAunXr5iaoDn843LOn6GD79u2JyW/Xrl3dBNXhz8PEiRMT36WU5Ofns23bNhYtWsQ777zD6tWr2b59O0888QQTJkxg3LhxNG3a9E9ssYPDXxv2NZsKtWrVSnyXUv6ezXFw+FVwz56Dgw4dOnDllVfuM9+RRx4JqEn3X3HiXRjcBLWIorC3IyklU6dO5eabb2b58uWsW7eOiy++mLlz53LyySf/cY10cHBI4K/AaDgUf7hnz8FBnTp13D3hN8I5SR2CEELQqlUrPvroo8Sba35+PldccQV79+79k1vn4ODg4FAc4Z49Dn8k3AT1EEbZsmUZP358wjZt2bJlvPTSS5F5CwoKmDhxIjfccANNmjShQoUKlCpVioyMDE4++WR69erFJ598krKu5s2bI4RI2ACB74Vs/4eXKPbs2cObb75Jv379aNq0KZmZmZQuXZr09HROOOEEunXrxjvvvPObx8LB4VBDKi/+bt26IYRI2J8C1KxZM+la69at2x/bYAcHDffsOfhI5cW/evVqhBDUrFkzkTZ69OjIMTBOZMUFbon/EEeVKlXo0aMHDz/8MADPPvssHTt2TMp30kknRf54d+zYwdKlS1m6dClDhgyhf//+PPDAAwetfeeff36kZMbu3btZvnw5y5cvZ/To0XTt2pVhw4ZRunTpg1a3g4ODg8PvA/fscfi94SaoxQCdOnVK3CTee+89du/eTalSpQJ5CgoKKF++POeffz4NGjSgatWqlCpVim+++YaFCxcyfvx4du/ezYABA8jMzOTmm28OHP9///d/bN68mX/961989tlnQLRTSJ06dZLqPeqoozjvvPNo1KgRNWrU4PDDD2f9+vV89tlnjB07lh9++IHRo0dTrly5/ZbgcXAorujTpw+XXXYZTz75JG+//TYAQ4cOJTMzM5DvuOOO+zOa5+CQgHv2/P7IzMxk4sSJbNq0iZ49ewJwzjnn0KdPn8i8xQrSocgASPwfCPbs2SOPPPLIxLGLFi1KyjN16lS5e/fulGWsXr1a1qlTRwIyPT1d7tixIzJfs2bNDqiNM2bMkD/++GPK/Zs3b5ZNmzaVgCxRooT86quv9qtcB4eigF97zdrHNmvWLHJ/165dE3lWrVr12xrq4FAI3LPn4Dx7PM9LtNHzvF997Ntvv520f9WqVYn9Xbt2/c1tPRTgbFCLAQ477DCqVauW2P7uu++S8rRo0YKSJVMT5llZWTzzzDMA7Ny5k8mTJx+Utp133nmFCjpXqFAhYWu3d+9exo4de1DqdXD4oxFlE+ZsRh2KM9yzJzXi8fg+7wnFzWb0YMNNUIsJjj766MT3LVu2/KoyzjzzzMT3efPm/eY27S+OP/54Kleu/IfX6+Dg4ODw2+CePQ6/F5wNajGBLfGRKmzdpk2beO6555g+fTpLly5l27Zt/Pjjj5F5161bd9DatmPHDsaOHcsbb7zB4sWL2bx5Mz/88MPvXq+Dwx+JfQn1O5tRh+II9+yJxv4I9Rc7m9GDDDdBLSbYvn174nv58uWT9r/00kv07NmT/Pz8/Spvx44dB6Vdb7/9Np06dWLDhg1/aL0ODn80nCi3w18RxfnZM3369JQTaSj8mndC/b8dboJaDPDLL78E3v6OOeaYwP533nmHTp06Jd50GzZsSF5eHrVq1SIjI4MyZcok8rZp0yZR5m/F8uXLueiiiygoKADgxBNPpGXLltSuXZvy5ctz+OGHJ/L26NGD77777qDU6+Dg4ODw+6O4P3t69OjBmjVrUtYjXbji3xVugloMsHjx4sRb3pFHHpkUci4WiyVuEMOGDeMf//hHZDmplj5+LQYMGJC4Qdx1113cd999KZeAUrXJwcHBwaFowj17HH5PuAlqMcALL7yQ+H7mmWcGPCZ//vln3n33XQBOPfXUQi/Gwt4Ufw1mzJgBKDube++9N+UNYufOnWzduvWg1u3g4ODg8PuiuD97nJf9nwvnxX+IY/369fznP/9JbHfv3j2wf8uWLezZsweAWrVqFVrWm2++uc/6SpTwfzL7Wt7YuHEjoMI02seFMWPGDBfH2cEhhAO51hwc/mi4Z88fi7/i/cBNUA9h7Ny5k/bt2yeM1OvWrcsVV1wRyHPEEUckvq9cubLQsgYOHLjPOo866qjE930ty5i6v/rqq5QX1C+//HJQw9s5OBQXHMi15uDwR8I9e/54/BXvB26CeghCSsnUqVM59dRTmTNnDgBly5blv//9b9LbYkZGBrVr1wbgo48+ipTC+f7777niiitYu3btPuuuWbNm4vvChQsLzXvaaacBSrw5Kozc7t27+cc//sFHH320z3odHP5qOJBrzcHhj4B79vx5KF++PBkZGQAsWrToL8GiOhvUIopJkyYlvkspE7YyixYt4p133mHVqlWJ/dWqVWPcuHFJBuoGvXv3TsTtbdeuHZ07d6Zp06akp6ezZMkSRo0axbfffstVV13Fc889V2i7zjvvPJ588klALen07duXrKwsDjvsMACys7PJzs5O1PvWW28BcMsttzBr1iwuvPBCKlSowPLly3nuuedYvnw555xzDsuXL3caqA4OFs4777zE99tuu43vvvuOE088MWHnV7VqVXJycv6s5jkUU7hnT9HFueeey8SJE1m5ciUdOnTg8ssvp1y5con9zZo1KzR61iGHPyfCqkMUsOIh789/uXLlZJ8+feS2bdsKLXfv3r2yc+fOhZZ16aWXyh9//HGf8cH37NmTiF8c9R+OP9y/f/9C683NzZWbNm2SWVlZEpBZWVkHZSwdHP4I2L/lX3tsqmtNSik7duyY8tr5q8Tjdvj94Z49B+fZ43leyvYcyLFvv/12ZJ6PP/5YpqWlpezTqlWrfnMfihIcg3oIoFSpUpQtW5ayZctSo0YNGjZsSJMmTbj44ov3621JCMHzzz/PRRddxH/+8x8+/vhjfvzxRzIzM6lfvz5dunShffv2+9WWww47jLfeeosnnniCyZMn8/nnn7Njx46U2nUPPPAAZ599Nk899RTz5s0jPz+fihUrUrduXTp27Ei3bt0KjdPs4PBXxpgxYzj77LN56aWXWLJkCdu3b084njg4/N5wz56ihfr167NgwQIee+wx3n33XdauXVtoIIFDHULKv4Ahg4ODg4ODg4ODwyED5yTl4ODg4ODg4OBQpOAmqA4ODg4ODg4ODkUKboLq4ODg4ODg4OBQpOAmqA4ODg4ODg4ODkUKboLq4ODg4ODg4OBQpOAmqA4ODg4ODg4ODkUKRV4EbO/evXz77be
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 800x800 with 6 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqgAAAKvCAYAAACmvT3kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOydfZxN1f7H30ukiTGiRp4aMmTUSJP0QEl0UamUVEqIS90ukW55uNnn1C/UTSo9kUpPKjdRKil6pCJJFF0PGXkWalBTkfX7Y6119tr77DMoMab18fKac9Zee+211z777HU+6/P9fIWUUuLg4ODg4ODg4OBQTFDqQHfAwcHBwcHBwcHBwYaboDo4ODg4ODg4OBQruAmqg4ODg4ODg4NDsYKboDo4ODg4ODg4OBQruAmqg4ODg4ODg4NDsYKboDo4ODg4ODg4OBQruAmqg4ODg4ODg4NDsYKboDo4ODg4ODg4OBQruAmqg4ODg4ODg4NDsYKboDo4ODgcAAghEEJw9tlnH+iuODg4HGDEYrHEd8J77713oLtTLFD6QHfAwYcQIrL80EMPpUKFCmRkZJCVlUVeXh6nnnoq559/Pmlpafu1j+PGjSM/Px9QN5SDw18Zqe7ZVLjxxhu577779qju5MmTmT9/PgB9+/alYsWKe9c5B4c9hHv27BvEYjHi8fhe7fP999/v0b39ww8/JL47GjVqxMUXX7z3HTzI4CaoBwF+/fVXNm3axKZNm1i+fDnvvPMOABUrVqRLly7E43EyMjL2S1/GjRvH+++/DxTfLwkHh5KAyZMn89RTTwHQtWtXN0F12O9wz57igx9++CEx+e3SpYuboDocOEyaNCnxWkpJQUEB33//PfPnz+eDDz4gPz+fH374gfvvv5+JEyfy/PPP06xZswPYYweHvzbsezYV6tSpk3gtpfwzu+Pg8Lvgnj37BpdffjlXXHHFbuuVK1cOUJPuv+LEuyi4CWoxRVG/jqSUTJ06lb59+7J06VJWr17NBRdcwKxZszj++OP3XycdHBwS+CswGg4lH+7Zs29Qv359953wB+GCpA5CCCE477zzmDt3buKXa0FBAZdddhm7du06wL1zcHBwcCiJcM8eh/0JN0E9iFGhQgUmTJiQ0KYtXryYF198MbJuYWEhkyZN4oYbbuDUU0+lcuXKlClThoyMDI4//niuv/56vvjii5THOvvssxFCJDRA4Ech2//DSxQ7d+5k2rRp9O/fn2bNmpGZmcmhhx5Keno69erVo2vXrnzwwQd/eCwcHA42pIri79q1K0KIhP4UoHbt2kn3WteuXfdvhx0cNNyzZ98jVRR/fn4+Qghq166dKHvqqacix8AEkZUUuCX+gxxVq1alZ8+e3H333QA88cQTXHnllUn1GjRoEPnh3bp1K4sWLWLRokU8+uijDBw4kKFDh+6z/p177rmRlhk7duxg6dKlLF26lKeeeoouXbowZswYDj300H12bAcHBweHPwfu2ePwZ8NNUEsAOnXqlPiS+Oijj9ixYwdlypQJ1CksLKRSpUqce+65nHTSSVSvXp0yZcqwZs0a5s2bx4QJE9ixYwfDhg0jMzOTvn37Bvb/v//7PzZt2sS///1vvvrqKyA6KKR+/fpJxy1fvjwtW7bk5JNPplatWhx22GGsW7eOr776iueee44ff/yRp556iooVK+6xBY+DQ0lFnz59uPjii3nggQd49913ARg9ejSZmZmBesccc8yB6J6DQwLu2fPnIzMzk0mTJrFx40Z69eoFQIsWLejTp09k3RIF6VBsACT+7w127twpy5Url9h3/vz5SXWmTp0qd+zYkbKN/Px8Wb9+fQnI9PR0uXXr1sh6zZs336s+Tp8+Xf70008pt2/atEk2a9ZMArJUqVLym2++2aN2HRyKA37vPWvv27x588jtXbp0SdRZsWLFH+uog0MRcM+effPs8Twv0UfP8373vu+++27S9hUrViS2d+nS5Q/39WCA06CWABxyyCHUqFEj8f67775LqtOmTRtKl05NmGdlZfHwww8DsG3bNl555ZV90reWLVsWaehcuXLlhNZu165dPPfcc/vkuA4O+xtRmjCnGXUoyXDPntSIx+O7/U4oaZrRfQ03QS0hOOKIIxKvN2/e/LvaOOOMMxKvZ8+e/Yf7tKc49thjOfroo/f7cR0cHBwc/hjcs8fhz4LToJYQ2BYfqdLWbdy4kaeffpq33nqLRYsW8f333/PTTz9F1l29evU+69vWrVt57rnneOONN1i4cCGbNm3ixx9//NOP6+CwP7E7o36nGXUoiXDPnmjsiVF/idOM7mO4CWoJwQ8//JB4XalSpaTtL774Ir169aKgoGCP2tu6des+6de7775Lp06dWL9+/X49roPD/oYz5Xb4K6IkP3veeuutlBNpKPqed0b9fxxugloC8NtvvwV+/R111FGB7R988AGdOnVK/NLNy8ujVatW1KlTh4yMDMqWLZuo2759+0SbfxRLly7l/PPPp7CwEIDjjjuOtm3bUrduXSpVqsRhhx2WqNuzZ0++++67fXJcBwcHB4c/HyX92dOzZ09WrlyZ8jjSpSv+U+EmqCUACxcuTPzKK1euXFLKuVgslviCGDNmDH//+98j20m19PF7MWzYsMQXxODBg7njjjtSLgGl6pODg4ODQ/GEe/Y4/JlwE9QSgPHjxyden3HGGYGIyV9//ZUPP/wQgMaNGxd5Mxb1S/H3YPr06YDS2dx+++0pvyC2bdvGli1b9umxHRwcHBz+XJT0Z4+Lsj+wcFH8BznWrVvHY489lnjfvXv3wPbNmzezc+dOAOrUqVNkW9OmTdvt8UqV8j8yu1ve2LBhA6DSNNr7hTF9+nSXx9nBIYS9udccHPY33LNn/+Kv+H3gJqgHMbZt20bHjh0TIvWcnBwuu+yyQJ3DDz888Xr58uVFtjVy5MjdHrN8+fKJ17tbljHH/uabb1LeUL/99ts+TW/n4FBSsDf3moPD/oR79ux//BW/D9wE9SCElJKpU6fSuHFjZs6cCUCFChX473//m/RrMSMjg7p16wIwd+7cSCuc7du3c9lll7Fq1ardHrt27dqJ1/PmzSuy7imnnAIo8+aoNHI7duzg73//O3Pnzt3tcR0c/mrYm3vNwWF/wD17DhwqVapERkYGAPPnz/9LsKhOg1pMMXny5MRrKWVCKzN//nw++OADVqxYkdheo0YNnn/++SSBukHv3r0TeXs7dOjAVVddRbNmzUhPT+fLL79k3LhxrF27lmuuuYann366yH61bNmSBx54AFBLOv369SMrK4tDDjkEgOzsbLKzsxPHffvttwG46aabeO+992jdujWVK1dm6dKlPP300yxdupQWLVqwdOlS54Hq4GChZcuWide33HIL3333Hccdd1xC51e9enVyc3MPVPccSijcs6f44pxzzmHSpEksX76cyy+/nEsuuYSKFSsmtjdv3rzI7FkHHQ5MhlWHKGDlQ96T/xUrVpR9+vSR33//fZHt7tq1S1511VVFtnXRRRfJn376abf5wXfu3JnIXxz1P5x/eODAgUUet2nTpnLjxo0yKytLAjIrK2ufjKWDw/6A/Vn+vfumuteklPLKK69Mee/8VfJxO/z5cM+effPs8TwvZX/2Zt933303ss7nn38u09LSUp7TihUr/vA5FCc4BvUgQJkyZahQoQIVKlSgVq1a5OXlceqpp3LBBRfs0a8lIQTPPvss559/Po899hiff/45P/30E5mZmTRq1IjOnTvTsWPHPerLIYccwttvv83999/PK6+8wtdff83WrVtTetcNHTqUs846iwcffJDZs2dTUFDAkUceSU5ODldeeSVdu3YtMk+zg8NfGc888wxnnXUWL774Il9++SU//PBDIvDEweHPhnv2FC80atSIzz77jHvvvZcPP/yQVatWFZlI4GCHkPIvIGRwcHBwcHBwcHA4aOCCpBwcHBwcHBwcHIoV3ATVwcHBwcHBwcGhWMFNUB0cHBwcHBwcHIoV3ATVwcHBwcHBwcGhWMFNUB0cHBwcHBwcHIoV3ATVwcHBwcHBwcGhWKHYm4Dt2rW
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 800x800 with 6 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArUAAAKvCAYAAAB07Te8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOydf7wN1f7/n6v86FQcl+7xu6McoZsSIj8iOflRKSKFitL1o0JuvkK3Zu+6oT6VpG50E/pBdRNdlUghSgkJpUKO/FbkUJ1E1vePtdaeNbNnHyrpONbTw2OfPXvNmjVr9uxZ85r3er2FlFLicDgcDofD4XAcxRz3ZzfA4XA4HA6Hw+H4vbhBrcPhcDgcDofjqMcNah0Oh8PhcDgcRz1uUOtwOBwOh8PhOOpxg1qHw+FwOBwOx1GPG9Q6HA6Hw+FwOI563KDW4XA4HA6Hw3HU4wa1DofD4XA4HI6jHjeodTgcDofD4XAc9bhBrcPhcBwlCCEQQnDhhRf+2U1xOBx/MrFYLPGbMHfu3D+7OQWCIn92Axy/DyFE5PJixYpRsmRJ0tPTyczMpE6dOjRo0IBLL72UtLS0I9rGCRMmkJOTA6iT0OE4lkl1zqaif//+PPLII4dUdtq0aSxbtgyA2267jVKlSv26xjkch4i79hweYrEY8Xj8V63z3XffHdK5vWvXrsRvR+3atWnXrt2vb+BRhhvUFlJ+/vlnvv32W7799lvWrl3LO++8A0CpUqXo1q0b8Xic9PT0I9KWCRMmMG/ePKDg/rA4HIWBadOmMXHiRAC6d+/uBrWOI4679hQcdu3alRgwd+vWzQ1qHUcXU6dOTfwtpSQ3N5fvvvuOZcuW8e6775KTk8OuXbsYNWoUU6ZMYfLkyTRp0uRPbLHDcWxjn7OpqFq1auJvKeUf2RyH4zfhrj2Hh6uvvpprrrnmoOVOOukkQA3Uj8XBen64QW0hIr+7MCklM2bM4LbbbmP16tVs3LiRyy67jPfee4+//e1vR66RDocjwbGgnDgKP+7ac3ioUaOG+034nbiJYscIQgguueQSFi9enLhDzs3N5aqrruLAgQN/cuscDofDURhx1x7HkcQNao8xSpYsyUsvvZSItVu1ahUvvvhiZNm8vDymTp3KLbfcQoMGDShTpgxFixYlPT2dv/3tb/Tp04dPPvkk5bYuvPBChBCJmCbwZ2/b/8OPT/bv38/MmTO5/fbbadKkCRkZGRQrVowSJUpwxhln0L17d959993f3RcOx9FGKveD7t27I4RIxNMCnHbaaUnnWvfu3Y9sgx0Ojbv2HH5SuR/k5OQghOC0005LLJs4cWJkH5iJdIUFF35wDFK+fHl69uzJAw88AMDTTz9N586dk8qdeeaZkV/43bt389lnn/HZZ58xZswYhgwZwrBhww5b+y6++OJIe5J9+/axevVqVq9ezcSJE+nWrRtPPvkkxYoVO2zbdjgcDscfg7v2OP5o3KD2GKVLly6JH5b333+fffv2UbRo0UCZvLw8SpcuzcUXX8y5555LxYoVKVq0KJs2bWLp0qW89NJL7Nu3j+HDh5ORkcFtt90WWP9f//oX3377Lf/85z/59NNPgeiJMTVq1Eja7sknn0yLFi2oW7cuVapU4YQTTmDLli18+umnPP/88/zwww9MnDiRUqVKHbLdkcNRWOnXrx/t2rXj0UcfZc6cOQCMHTuWjIyMQLlTTz31z2iew5HAXXv+eDIyMpg6dSrbt2+nV69eADRv3px+/fpFli1USMdRDZD4/2vYv3+/POmkkxLrLlu2LKnMjBkz5L59+1LWkZOTI2vUqCEBWaJECbl79+7Ics2aNftVbZw9e7b88ccfU37+7bffyiZNmkhAHnfccfKrr746pHodjoLAbz1n7XWbNWsW+Xm3bt0SZdatW/f7Gupw5IO79hyea4/neYk2ep73m9edM2dO0ufr1q1LfN6tW7ff3dajARdTe4xy/PHHU6lSpcT7b775JqlM69atKVIktZifmZnJv//9bwD27NnDq6++elja1qJFi3xNusuUKZOIHTxw4ADPP//8Ydmuw3GkiYpxczGwjsKMu/akJh6PH/Q3obDFwB5u3KD2GOYvf/lL4u8dO3b8pjoaNWqU+PvDDz/83W06VE4//XTKlSt3xLfrcDgcjt+Hu/Y4/ihcTO0xjG2nkirl4fbt23nmmWeYNWsWn332Gd999x0//vhjZNmNGzcetrbt3r2b559/njfeeIMVK1bw7bff8sMPP/zh23U4jiQHS77gYmAdhRF37YnmUJIvFLoY2MOMG9Qew+zatSvxd+nSpZM+f/HFF+nVqxe5ubmHVN/u3bsPS7vmzJlDly5d2Lp16xHdrsNxpHFG645jkcJ87Zk1a1bKwTfkf8675Au/HzeoPUb55ZdfAneZf/3rXwOfv/vuu3Tp0iVxR12nTh2ys7OpWrUq6enpFC9ePFG2ffv2iTp/L6tXr+bSSy8lLy8PgOrVq9OmTRuqVatG6dKlOeGEExJle/bsyTfffHNYtutwOByOP57Cfu3p2bMn69evT7kd6VJd/6G4Qe0xyooVKxJ3kyeddFJSusJYLJb4UXnyySf5+9//HllPqscyv5Xhw4cnflTuvPNO7r333pSPp1K1yeFwOBwFE3ftcfyRuEHtMcqkSZMSfzdq1Cgw0/Tnn39m/vz5ANSrVy/fEzi/O9LfwuzZswEVN3TPPfek/FHZs2cPO3fuPKzbdjgcDscfS2G/9jh3gj8X535wDLJlyxb+85//JN736NEj8PmOHTvYv38/AFWrVs23rpkzZx50e8cd53/NDvboZdu2bYBK8WmvF2b27Nkub7jDEeLXnGsOx5HGXXuOLMfi74Eb1B5j7Nmzh06dOiUC9WvWrMlVV10VKHPiiScm/l67dm2+dY0cOfKg2zz55JMTfx/skZHZ9ldffZXyJPzll18Oa2pEh6Ow8GvONYfjSOKuPUeeY/H3wA1qjxGklMyYMYN69eqxYMECAEqWLMl///vfpLvS9PR0qlWrBsDixYsjbYe+//57rrrqKjZs2HDQbZ922mmJv5cuXZpv2fPOOw9QhtxRKQj37dvH3//+dxYvXnzQ7Tocxxq/5lxzOI4E7trz51G6dGnS09MBWLZs2TGh1rqY2kLEtGnTEn9LKROxP8uWLePdd99l3bp1ic8rVarE5MmTk4L0DX379k3kie7YsSNdu3alSZMmlChRgpUrVzJhwgQ2b97M9ddfzzPPPJNvu1q0aMGjjz4KqMdNAwYMIDMzk+OPPx6ArKwssrKyEtt96623APjHP/7B3LlzadWqFWXKlGH16tU888wzrF69mubNm7N69WrnUetwWLRo0SLx96BBg/jmm2+oXr16Im6xYsWK1KpV689qnqOQ4q49BZeLLrqIqVOnsnbtWq6++mquvPJKSpUqlfi8WbNm+WZRO+r4c7LzOg4XWPm3D+V/qVKlZL9+/eR3332Xb70HDhyQXbt2zbeuK664Qv74448HzUe/f//+RL7sqP/hfNdDhgzJd7uNGzeW27dvl5mZmRKQmZmZh6UvHY4jgf1d/q3rpjrXpJSyc+fOKc+dYyX/u+OPx117Ds+1x/O8lO35NevOmTMnsszHH38s09LSUu7TunXrfvc+FCScUltIKVq0KCVLlqRkyZJUqVKFOnXq0KBBAy677LJDuisTQvDcc89x6aWX8p///IePP/6YH3/8kYyMDGrXrs11111Hp06dDqktxx9/PG+99RajRo3i1Vdf5fPPP2f37t0pvQWHDRtG06ZNeeyxx/jwww/Jzc3llFNOoWbNmnTu3Jnu3bvnmxfc4TiWefbZZ2natCkvvvgiK1euZNeuXYnJNw7HH4279hQsateuzZIlS3j44YeZP38+GzZsyDc5xNGOkPIYCLJwOBwOh8PhcBRq3EQxh8PhcDgcDsdRjxvUOhwOh8PhcDiOetyg1uFwOBwOh8Nx1OMGtQ6Hw+FwOByOox43qHU4HA6Hw+FwHPW4Qa3D4XA4HA6H46jn6DJcO0QOHDjA5s2bKVGiBEKIP7s5DsdhQepMPRUqVEhKL/ln4c41R2G
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 800x800 with 6 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqgAAAKvCAYAAACmvT3kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOyde3wVxfnGv6tBjQrhVhABCUrQqEGICFVApMRKtFGDigLFYEkVtKgUfyp4mXO0gm2leEcsFBARpchFqnjBAipakZui0ApIEJBLuQVUVND5/TEzu7N79gRU1JDOw4dPzpndnZ3dPbs788zzPq8npZQ4ODg4ODg4ODg4VBIc8lM3wMHBwcHBwcHBwcGG66A6ODg4ODg4ODhUKrgOqoODg4ODg4ODQ6WC66A6ODg4ODg4ODhUKrgOqoODg4ODg4ODQ6WC66A6ODg4ODg4ODhUKrgOqoODg4ODg4ODQ6WC66A6ODg4ODg4ODhUKrgOqoODg4ODg4ODQ6WC66A6ODg4/ATwPA/P8zjnnHN+6qY4ODj8xEgkEv4zYc6cOT91cyoFMn7qBjgE8Dwvtvywww6jRo0aZGVl0aRJE/Lz82nbti0XXHABmZmZP2obx44dS1lZGaBuKAeH/2Wku2fT4YYbbuD+++/fr3WnTZvGkiVLALjxxhupWbPmt2ucg8N+wr17DgwSiQTJZPJbbbN9+/b9urd37NjhPztatmzJxRdf/O0beJDBdVAPAnz11Vds2bKFLVu2sGrVKv75z38CULNmTUpKSkgmk2RlZf0obRk7dixz584FKu9DwsGhKmDatGmMGzcOgN69e7sOqsOPDvfuqTzYsWOH3/ktKSlxHVSHnw5Tp071P0spKS8vZ/v27SxZsoTXXnuNsrIyduzYwQMPPMCzzz7LxIkTad++/U/YYgeH/23Y92w6nHDCCf5nKeUP2RwHh+8E9+45MLj88su54oor9rneUUcdBahO9/9ix7siuA5qJUVFoyMpJTNnzuTGG29kxYoVrFu3jl/96lfMmzePU0455cdrpIODg4//BUbDoerDvXsODE466ST3TPiecEFSByE8z+P8889nwYIF/si1vLycyy67jG+++eYnbp2Dg4ODQ1WEe/c4/JhwHdSDGDVq1GDSpEm+Nm358uU888wzsevu3r2bqVOnct1119G2bVvq1KlDtWrVyMrK4pRTTqFfv368++67afd1zjnn4HmerwGCIArZ/h+doti7dy8vvfQSAwcOpH379tSrV4/DDjuM6tWr07x5c3r37s1rr732vc+Fg8PBhnRR/L1798bzPF9/CtC0adOUe613794/boMdHDTcu+fAI10Uf1lZGZ7n0bRpU79s3LhxsefABJFVFbgp/oMcDRo04Oqrr+ZPf/oTAH/729/o3r17ynonn3xy7I93586dLFu2jGXLlvHYY48xaNAghgwZcsDad+6558ZaZuzZs4cVK1awYsUKxo0bR0lJCY8//jiHHXbYAdu3g4ODg8MPA/fucfih4TqoVQA9evTwHxJvvvkme/bsoVq1aqF1du/eTe3atTn33HNp1aoVDRs2pFq1aqxfv55FixYxadIk9uzZw9ChQ6lXrx433nhjaPs//OEPbNmyhdtvv50PPvgAiA8KOemkk1L2e/TRR9O5c2dOP/10srOzOeKII9iwYQMffPABEyZM4LPPPmPcuHHUrFlzvy14HByqKq6//nouvvhiHnzwQWbPng3AyJEjqVevXmi944477qdonoODD/fu+eFRr149pk6dyubNm7nmmmsA6NSpE9dff33sulUK0qHSAPD/fxvs3btXHnXUUf62S5YsSVln5syZcs+ePWnrKCsrkyeddJIEZPXq1eXOnTtj1+vYseO3auOsWbPk559/nnb5li1bZPv27SUgDznkEPnRRx/tV70ODpUB3/Wetbft2LFj7PKSkhJ/ndWrV3+/hjo4VAD37jkw7x4hhN9GIcR33nb27Nkpy1evXu0vLykp+d5tPRjgNKhVAIceeiiNGjXyv//3v/9NWadLly5kZKQnzJs0acKjjz4KwK5du5g+ffoBaVvnzp0rNHSuU6eOr7X75ptvmDBhwgHZr4PDj404TZjTjDpUZbh3T3okk8l9PhOqmmb0QMN1UKsIatWq5X/eunXrd6rjrLPO8j+//fbb37tN+4vjjz+eY4455kffr4ODg4PD94N79zj8UHAa1CoC2+IjXdq6zZs388QTT/Dyyy+zbNkytm/fzueffx677rp16w5Y23bu3MmECRN44YUXWLp0KVu2bOGzzz77wffr4PBjYl9G/U4z6lAV4d498dgfo/4qpxk9wHAd1CqCHTt2+J9r166dsvyZZ57hmmuuoby8fL/q27lz5wFp1+zZs+nRowcbN278Uffr4PBjw5lyO/wvoiq/e15++eW0HWmo+J53Rv3fH66DWgXw9ddfh0Z/P/vZz0LLX3vtNXr06OGPdPPz8ykoKOCEE04gKyuLww8/3F+3uLjYr/P7YsWKFVxwwQXs3r0bgBNPPJHCwkJycnKoXbs2RxxxhL/u1VdfzX//+98Dsl8HBwcHhx8eVf3dc/XVV7NmzZq0+5EuXfEPCtdBrQJYunSpP8o76qijUlLOJRIJ/wHx+OOP89vf/ja2nnRTH98VQ4cO9R8Qt912G3fffXfaKaB0bXJwcHBwqJxw7x6HHxKug1oF8NRTT/mfzzrrrFDE5FdffcXrr78OQOvWrSu8GSsaKX4XzJo1C1A6m7vuuivtA2LXrl1s27btgO7bwcHBweGHRVV/97go+58WLor/IMeGDRv461//6n/v06dPaPnWrVvZu3cvACeccEKFdb300kv73N8hhwQ/mX1Nb2zatAlQaRrt7aKYNWuWy+Ps4BDBt7nXHBx+bLh3z4+L/8XngeugHsTYtWsX3bp180Xqubm5XHbZZaF1jjzySP/zqlWrKqxr+PDh+9zn0Ucf7X/e17SM2fdHH32U9ob6+uuvD2h6OweHqoJvc685OPyYcO+eHx//i88D10E9CCGlZObMmbRu3Zo33ngDgBo1avD3v/89ZbSYlZVFTk4OAAsWLIi1wvn000+57LLLWLt27T733bRpU//zokWLKlz3jDPOAJR5c1wauT179vDb3/6WBQsW7HO/Dg7/a/g295qDw48B9+756VC7dm2ysrIAWLJkyf8Ei+o0qJUU06ZN8z9LKX2tzJIlS3jttddYvXq1v7xRo0ZMnDgxRaBu0L9/fz9v76WXXkrPnj1p37491atX5/3332fs2LF88sknXHnllTzxxBMVtqtz5848+OCDgJrSGTBgAE2aNOHQQw8FoFmzZjRr1szf7yuvvALA73//e+bMmcN5551HnTp1WLFiBU888QQrVqygU6dOrFixwnmgOjhY6Ny5s//55ptv5r///S8nnniir/Nr2LAheXl5P1XzHKoo3Lun8uIXv/gFU6dOZdWqVVx++eV07dqVmjVr+ss7duxYYfasgw4/TYZVhzhg5UPen/81a9aU119/vdy+fXuF9X7zzTeyZ8+eFdZ10UUXyc8//3yf+cH37t3r5y+O+x/NPzxo0KAK99uuXTu5efNm2aRJEwnIJk2aHJBz6eDwY8D+LX/XbdPda1JK2b1797T3zv9KPm6HHx7u3XNg3j1CiLTt+Tbbzp49O3adxYsXy8zMzLTHtHr16u99DJUJjkE9CFCtWjVq1KhBjRo1yM7OJj8/n7Zt2/KrX/1qv0ZLnufx5JNPcsEFF/DXv/6VxYsX8/nnn1OvXj1atmxJr1696Nat23615dBDD+WVV17hgQceYPr06fz73/9m586dab3rhgwZwtlnn83DDz/M22+/TXl5OXXr1iU3N5fu3bvTu3fvCvM0Ozj8L2P8+PGcffbZPPPMM7z//vvs2LHDDzxxcPih4d49lQstW7Zk4cKF/OUvf+H1119n7dq1FSYSONjhSfk/IGRwcHBwcHBwcHA4aOCCpBwcHBwcHBwcHCoVXAfVwcHBwcHBwcGhUsF1UB0cHBwcHBwcHCoVXAfVwcHBwcHBwcGhUsF1UB0cHBwcHBwcHCoVXAfVwcH
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 800x800 with 6 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAq4AAAKvCAYAAACro02jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOydf7xVU/7/n0u3cumHfij9covCNWpyNYUijWsqXCNSZNI1NRJq5Gc/fGefY6QYv/NbPgr50aAIiRAVIpUyhaKSJPrhViPpZn3/eK+199r7nBszk9xu69Wjx7lnn7XXXnuds/da+7Ve79dbaa01Hh4eHh4eHh4eHuUce/3aDfDw8PDw8PDw8PD4OfATVw8PDw8PDw8Pj90CfuLq4eHh4eHh4eGxW8BPXD08PDw8PDw8PHYL+Imrh4eHh4eHh4fHbgE/cfXw8PDw8PDw8Ngt4CeuHh4eHh4eHh4euwX8xNXDw8PDw8PDw2O3gJ+4enh4eHh4eHh47BbwE1cPDw+PcgSlFEopTjjhhF+7KR4eHr8yUqlUeE+YPn36r92ccoGcX7sBHj8NpVTW7VWqVKFGjRrUrFmTvLw8CgoKaNeuHaeccgq5ubm7tI1jx45l+fLlgFxoHh57Msq6ZsvCX//6V2677bafVXbSpEnMnz8fgEsvvZT99tvvP2uch8fPhB97dg5SqRTpdPo/2mfDhg0/69r+9ttvw3tH69atOf300//zBu5m8BPX3Rg//PADa9euZe3atXz66ae89tprAOy333706dOHdDpNzZo1d0lbxo4dyxtvvAGU35uHh0dFwKRJkxg3bhwAxcXFfuLqscvhx57yg2+//TacFPfp08dPXD3KHyZOnBj+rbWmpKSEDRs2MH/+fN58802WL1/Ot99+y+23387TTz/N448/TocOHX7FFnt47Nlwr9mycPDBB4d/a61/yeZ4ePxX8GPPzkHPnj05++yzf7LcvvvuC8hkfE+ckO8IfuK6m2FHT1Naa6ZMmcKll17KkiVL+OKLLzj11FOZNWsWv/nNb3ZdIz08PELsCQyIR8WHH3t2Dg477DB/T/gf4YOzKhCUUpx88snMmTMnfNItKSnhrLPO4scff/yVW+fh4eHhURHhxx6PXQk/ca2AqFGjBhMmTAi1b4sXL+bJJ5/MWnbLli1MnDiRiy++mHbt2lGnTh0qV65MzZo1+c1vfsOAAQP44IMPyjzWCSecgFIq1BhBFBXt/k8udZSWljJ16lQuv/xyOnToQL169ahSpQrVq1fnkEMOobi4mDfffPN/7gsPj90NZbkKFBcXo5QK9a0AzZo1y7jWiouLd22DPTwM/Niz81GWq8Dy5ctRStGsWbNw27hx47L2gQ1eqyjwUoEKigYNGnDBBRdw4403AvB///d/nHPOORnlDj/88Kw/6o0bN7Jo0SIWLVrEvffey9ChQ7n++ut3WvtOOumkrNYe27ZtY8mSJSxZsoRx48bRp08f7r//fqpUqbLTju3h4eHh8cvAjz0evzT8xLUCo1evXuHN46233mLbtm1Urlw5VmbLli3Url2bk046iSOPPJJGjRpRuXJlVq1axdy5c5kwYQLbtm1j5MiR1KtXj0svvTS2/3XXXcfatWu55ppr+Ne//gVkD0Y57LDDMo5brVo1TjzxRI466iiaNm3K3nvvzerVq/nXv/7F+PHj+fe//824cePYb7/9frZVkIdHRcWgQYM4/fTTueOOO3j99dcBuO+++6hXr16s3IEHHvhrNM/DI4Qfe3551KtXj4kTJ/L111/Tv39/ADp16sSgQYOylq1Q0B7lHkD4/z9BaWmp3nfffcN958+fn1FmypQpetu2bWXWsXz5cn3YYYdpQFevXl1v3Lgxa7mOHTv+R22cNm2a/u6778r8fO3atbpDhw4a0HvttZf+7LPPfla9Hh7lAf/tNevu27Fjx6yf9+nTJyyzbNmy/62hHh47gB97ds7YEwRB2MYgCP7rfV9//fWMz5ctWxZ+3qdPn/+5rbsDvMa1AqNSpUo0btw4fP/NN99klOnSpQs5OWUT73l5edx9990AbNq0iWeffXantO3EE0/coVF1nTp1Qi3fjz/+yPjx43fKcT08djWyac68JtWjIsOPPWUjnU7/5D2homlSdzb8xLWCo1atWuHf69at+6/qOPbYY8O/Z8+e/T+36efioIMO4oADDtjlx/Xw8PDw+N/gxx6PXwpe41rB4VqRlJW+7+uvv+bhhx/m5ZdfZtGiRWzYsIHvvvsua9kvvvhip7Vt48aNjB8/nhdffJGFCxeydu1a/v3vf//ix/Xw2JX4qQQEXpPqURHhx57s+DkJCCqcJnUnw09cKzi+/fbb8O/atWtnfP7kk0/Sv39/SkpKflZ9Gzdu3Cntev311+nVqxdfffXVLj2uh8euhjcb99gTUZHHnpdffrnMCTbs+Jr3CQj+d/iJawXG9u3bY0+L+++/f+zzN998k169eoVPxgUFBRQWFnLwwQdTs2ZNqlatGpbt1q1bWOf/iiVLlnDKKaewZcsWAA499FC6du1KixYtqF27NnvvvXdY9oILLuCbb77ZKcf18PDw8PjlUdHHngsuuIAVK1aUeRzt0zb/ovAT1wqMhQsXhk+F++67b0bqvVQqFd447r//fv7yl79kraesJZT/FiNHjgxvHMOHD+fvf/97mUtJZbXJw8PDw6N8wo89Hr8k/MS1AuOxxx4L/z722GNjEZw//PADM2bMAKBNmzY7vEh39GT532DatGmA6HiuvfbaMm8cmzZtYv369Tv12B4eHh4evywq+tjjo/5/XXhXgQqK1atX88ADD4Tv+/btG/t83bp1lJaWAnDwwQfvsK6pU6f+5PH22iv6Kf3UMsmaNWsASVfp7pfEtGnTfJ5rD48E/pNrzcNjV8OPPbsWe+L9wE9cKyA2bdpEjx49QnF8fn4+Z511VqzMPvvsE/796aef7rCuW2+99SePWa1atfDvn1rescf+7LPPyrzQtm/fvlPT/Hl4VBT8J9eah8euhB97dj32xPuBn7hWIGitmTJlCm3atGHmzJkA1KhRg3/+858ZT5c1a9akRYsWAMyZMyerZc/mzZs566yzWLly5U8eu1mzZuHfc+fO3WHZ3/3ud4CYUmdLp7dt2zb+8pe/MGfOnJ88rofHnob/5Frz8NgV8GPPr4fatWtTs2ZNAObPn79HsK5e47qbYdKkSeHfWutQizN//nzefPNNli1bFn7euHFjHn/88QxhvMXAgQPDvMbdu3fn3HPPpUOHDlSvXp0PP/yQsWPH8uWXX3Leeefx8MMP77BdJ554InfccQcgS0ODBw8mLy+PSpUqAdC8eXOaN28eHveVV14B4LLLLmP69Ol07tyZOnXqsGTJEh5++GGWLFlCp06dWLJkifdw9fBwcOKJJ4Z/X3XVVXzzzTcceuihoY6wUaNGtGzZ8tdqnkcFhR97yi9+//vfM3HiRD799FN69uzJGWecwX777Rd+3rFjxx1mC9vt8OtkmvX4T4CTL/rn/N9vv/30oEGD9IYNG3ZY748//qjPPffcHdb1xz/+UX/33Xc/mT+9tLQ0zO+c7X8yP/PQoUN3eNz27dvrr7/+Wufl5WlA5+Xl7ZS+9PDYFXB/y//tvmVda1prfc4555R57ewp+co9fnn4sWfnjD1BEJTZnv9k39dffz1rmXnz5unc3Nwyz2nZsmX/8zmUJ3jGdTdG5cqVqVGjBjVq1KBp06YUFBTQrl07Tj311J/1dKWU4tFHH+WUU07hgQceYN68eXz33XfUq1eP1q1b07t3b3r06PGz2lKpUiVeeeUVbr/9dp599lk++ugjNm7cWKb33vXXX8/xxx/PnXfeyezZsykpKaFu3brk5+dzzjnnUFxcvMM81h4eezIeeeQRjj/+eJ588kk+/PBDvv322zDgxcPjl4Yfe8oXWrduzfvvv88tt9zCjBkzWLly5Q4TJOzuUFrvAYIIDw8PDw8PDw+P3R4+OMvDw8PDw8PDw2O3gJ+4enh4eHh4eHh47BbwE1cPDw8PDw8PD4/dAn7i6uHh4eHh4eHhsVvAT1w9PDw
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 800x800 with 6 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAq4AAAKvCAYAAACro02jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOydeXxVxf3+36cGMRUIm2jYEiCIUbFsghgQLHHBijWKqFgkFASXokT9qsDPzr1aQFstKu5iCSKiVARFxQUqqytbRUkRlCCbIAgBlSLg/P6YmXPmnHtvQIsawjy8eOXeOducc+45M/PM83k+npRS4uDg4ODg4ODg4FDB8atfugIODg4ODg4ODg4OBwLXcXVwcHBwcHBwcDgk4DquDg4ODg4ODg4OhwRcx9XBwcHBwcHBweGQgOu4Ojg4ODg4ODg4HBJwHVcHBwcHBwcHB4dDAq7j6uDg4ODg4ODgcEjAdVwdHBwcHBwcHBwOCbiOq4ODg4ODg4ODwyEB13F1cHBwqEDwPA/P8+jatesvXRUHB4dfGLFYzH8nzJ49+5euToVA2i9dAYf9w/O8pOVHHnkkNWrUICMjg6ysLNq0aUOHDh343e9+R3p6+s9ax+LiYkpLSwH1oDk4HM5I9cymwg033MB99913QOtOmzaNpUuXAjBkyBBq1qz5wyrn4HCAcG3PwUEsFiMej/+gbbZt23ZAz/b27dv9d0erVq248MILf3gFDzG4jushjO+++44tW7awZcsWPv30U/71r38BULNmTfr27Us8HicjI+NnqUtxcTFz5swBKu7Lw8GhMmDatGmMHz8egMLCQtdxdfjZ4dqeioPt27f7neK+ffu6jqtDxcPUqVP9z1JKysrK2LZtG0uXLmXu3LmUlpayfft27r//fqZMmcKkSZPo1KnTL1hjB4fDG/YzmwrNmjXzP0spf8rqODj8KLi25+Dg0ksv5bLLLtvvekcffTSgOuOHY4e8PLiO6yGG8kZTUkpmzJjBkCFDWLlyJevWreP8889nwYIFnHTSST9fJR0cHHwcDgyIQ+WHa3sODk444QT3Tvgf4YKzKhE8z+O8885j4cKF/ki3rKyMSy65hO+///4Xrp2Dg4ODQ2WEa3scfk64jmslRI0aNZg8ebKvfSspKeG5555Luu6uXbuYOnUq1113HR06dKBOnTpUqVKFjIwMTjrpJK655hr+/e9/pzxW165d8TzP1xhBEBVt/49Odezdu5fXX3+dm266iU6dOlGvXj2OPPJIqlevzvHHH09hYSFz5879n6+Fg8OhhlSuAoWFhXie5+tbAZo0aZLwrBUWFv68FXZw0HBtz8FHKleB0tJSPM+jSZMmftn48eOTXgMTvFZZ4KQClRSZmZkMHDiQv/71rwD84x//4PLLL09Y78QTT0z6o96xYwfLly9n+fLlPProowwdOpSRI0cetPqdddZZSa099uzZw8qVK1m5ciXjx4+nb9++PP744xx55JEH7dgODg4ODj8NXNvj8FPDdVwrMXr37u2/PN5++2327NlDlSpVQuvs2rWL2rVrc9ZZZ9G6dWsaNGhAlSpVWL9+PYsXL2by5Mns2bOHUaNGUa9ePYYMGRLa/i9/+Qtbtmzh//2//8fHH38MJA9GOeGEExKOW61aNbp160bbtm3Jzs7mqKOOYuPGjXz88cdMnDiRb775hvHjx1OzZs0DtgpycKisuP7667nwwgt54IEHeOuttwB47LHHqFevXmi9xo0b/xLVc3Dw4dqenx716tVj6tSpbN68mUGDBgFw5plncv311yddt1JBOlR4AP7/H4K9e/fKo48+2t926dKlCevMmDFD7tmzJ+U+SktL5QknnCABWb16dbljx46k63Xp0uUH1XHmzJny22+/Tbl8y5YtslOnThKQv/rVr+Rnn312QPt1cKgI+LHPrL1tly5dki7v27evv87q1av/t4o6OJQD1/YcnLZHCOHXUQjxo7d96623EpavXr3aX963b9//ua6HApzGtRLjiCOOoGHDhv73L7/8MmGdc889l7S01MR7VlYWDz/8MAA7d+7kxRdfPCh169atW7lG1XXq1PG1fN9//z0TJ048KMd1cPi5kUxz5jSpDpUZru1JjXg8vt93QmXTpB5suI5rJUetWrX8z1u3bv1R+zj99NP9z++9997/XKcDRdOmTTnuuON+9uM6ODg4OPxvcG2Pw08Fp3Gt5LCtSFKl79u8eTNPPfUUb7zxBsuXL2fbtm18++23Sdddt27dQavbjh07mDhxIq+++irLli1jy5YtfPPNNz/5cR0cfk7sLwGB06Q6VEa4tic5DiQBQaXTpB5kuI5rJcf27dv9z7Vr105Y/txzzzFo0CDKysoOaH87duw4KPV666236N27N1988cXPelwHh58bzmzc4XBEZW573njjjZQdbCj/mXcJCP53uI5rJca+fftCo8VjjjkmtHzu3Ln07t3bHxm3adOG/Px8mjVrRkZGBlWrVvXXLSgo8Pf5v2LlypX87ne/Y9euXQC0aNGC7t2707x5c2rXrs1RRx3lrztw4EC+/PLLg3JcBwcHB4efHpW97Rk4cCBr1qxJeRzp0jb/pHAd10qMZcuW+aPCo48+OiH1XiwW818cjz/+OFdddVXS/aSaQvmxGDVqlP/iGD58OHfeeWfKqaRUdXJwcHBwqJhwbY/DTwnXca3EeOaZZ/zPp59+eiiC87vvvmPevHkAtGvXrtyHtLyR5Y/BzJkzAaXjueOOO1K+OHbu3MlXX311UI/t4ODg4PDTorK3PS7q/5eFcxWopNi4cSNPPPGE/71///6h5Vu3bmXv3r0ANGvWrNx9vf766/s93q9+FfyU9jdNsmnTJkClq7S3i2LmzJkuz7WDQwQ/5FlzcPi54dqenxeH4/vAdVwrIXbu3EmvXr18cXxubi6XXHJJaJ1f//rX/udPP/203H2NHj16v8esVq2a/3l/0zvm2J999lnKB23fvn0HNc2fg0NlwQ951hwcfk64tufnx+H4PnAd10oEKSUzZsygXbt2zJ8/H4AaNWrwz3/+M2F0mZGRQfPmzQFYuHBhUsuer7/+mksuuYS1a9fu99hNmjTxPy9evLjcdU899VRAmVInS6e3Z88errrqKhYuXLjf4zo4HG74Ic+ag8PPAdf2/HKoXbs2GRkZACxduvSwYF2dxvUQw7Rp0/zPUkpfi7N06VLmzp3L6tWr/eUNGzZk0qRJCcJ4g8GDB/t5jXv27MkVV1xBp06dqF69Oh999BHFxcVs2LCBK6+8kqeeeqrcenXr1o0HHngAUFNDRUVFZGVlccQRRwCQk5NDTk6Of9w333wTgBtvvJHZs2dzzjnnUKdOHVauXMlTTz3FypUrOfPMM1m5cqXzcHVwsNCtWzf/8y233MKXX35JixYtfB1hgwYNaNmy5S9VPYdKCtf2VFz89re/ZerUqXz66adceumlXHTRRdSsWdNf3qVLl3KzhR1y+GUyzTr8EGDliz6Q/zVr1pTXX3+93LZtW7n7/f777+UVV1xR7r5+//vfy2+//Xa/+dP37t3r53dO9j+an3no0KHlHjcvL09u3rxZZmVlSUBmZWUdlGvp4PBzwP4t/9htUz1rUkp5+eWXp3x2Dpd85Q4/PVzbc3DaHiFEyvr8kG3feuutpOssWbJEpqenpzyn1atX/8/nUJHgGNdDGFWqVKFGjRrUqFGD7Oxs2rRpQ4cOHTj//PMPaHTleR5PP/00v/vd73jiiSdYsmQJ3377LfXq1aNVq1b06dOHXr16HVBdjjjiCN58803uv/9+XnzxRf7zn/+wY8eOlN57I0eO5IwzzuDBBx/kvffeo6ysjLp165Kbm8vll19OYWFhuXmsHRwOZ0yYMIEzzjiD5557jo8++ojt27f7AS8ODj81XNtTsdCqVSsWLVrE3//+d+bNm8fatWvLTZBwqMOT8jAQRDg4ODg4ODg4OBzycMFZDg4ODg4ODg4OhwRcx9XBwcHBwcHBweGQgOu4Ojg4ODg4ODg4HBJwHVcHBwcHBwcHB4dDAq7j6uDg4ODg4ODgcEjAdVwdHBwcHBwcHBw
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(0,shape[0]):\n",
    "    for j in range(0,shape[1]):\n",
    "        bval = result[i][j].best_values\n",
    "        fit = density_profile_BEC_2d(X,Y, **bval)\n",
    "        vmax = np.max(fit)\n",
    "\n",
    "        fig, axs = plt.subplots(2, 3, figsize=(8, 8))\n",
    "\n",
    "        ax = axs[ 0,0]\n",
    "        art = ax.pcolormesh(X, Y, cropOD[i,j], vmin=0, vmax=vmax, cmap='jet', shading='auto')\n",
    "        #plt.colorbar(art, ax=ax, label='z')\n",
    "        ax.set_title('Data')\n",
    "\n",
    "        # Plot gaussian 2d Fit + legend including Width parameters\n",
    "        ax = axs[0,1]\n",
    "\n",
    "        art = ax.pcolormesh(X, Y, fit, vmin=0, vmax=vmax, cmap='jet', shading='auto')\n",
    "        #plt.colorbar(art, ax=ax, label='z')\n",
    "        ax.set_title('Fit')\n",
    "\n",
    "        ax = axs[0,2]\n",
    "\n",
    "        art = ax.pcolormesh(X, Y, fit-cropOD[i,j], vmin=0, vmax=0.2, cmap='jet', shading='auto')\n",
    "        ax.set_title('Data-Fit')\n",
    "\n",
    "        ax = axs[1,0]\n",
    "\n",
    "        art = ax.plot(x, cropOD[i,j, round(center[i,j,1]), :])\n",
    "        ax.plot(x, fit[round(center[i,j,1]), :])\n",
    "        ax.plot(x, thermal(x, bval['x0_th'], bval['amp_th'], bval['sigmax_th']))\n",
    "\n",
    "        ax.set_title('cut along x')\n",
    "\n",
    "        ax = axs[1,1]\n",
    "\n",
    "        art = ax.plot(y, cropOD[i,j, :, round(center[i,j,0])])\n",
    "        ax.plot(y, fit[:, round(center[i,j,0])])\n",
    "        ax.plot(x, thermal(y, bval['y0_th'], bval['amp_th'], bval['sigmay_th']))\n",
    "        ax.set_title('cut along y')\n",
    "\n",
    "\n",
    "\n",
    "        plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-26T07:28:03.122325300Z",
     "start_time": "2023-07-26T07:27:55.722296100Z"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FWHM_x BEC: 8.88, FWHM_x thermal: 38.18\n",
      "FWHM_y BEC: 31.53, FWHM_y thermal: 35.28\n",
      "\n",
      "FWHM_x BEC: 9.75, FWHM_x thermal: 39.38\n",
      "FWHM_y BEC: 32.46, FWHM_y thermal: 34.81\n",
      "\n",
      "FWHM_x BEC: 10.09, FWHM_x thermal: 37.38\n",
      "FWHM_y BEC: 31.75, FWHM_y thermal: 33.83\n",
      "\n",
      "FWHM_x BEC: 8.57, FWHM_x thermal: 32.57\n",
      "FWHM_y BEC: 31.64, FWHM_y thermal: 28.49\n",
      "\n",
      "FWHM_x BEC: 9.10, FWHM_x thermal: 33.90\n",
      "FWHM_y BEC: 30.85, FWHM_y thermal: 29.66\n",
      "\n",
      "FWHM_x BEC: 9.19, FWHM_x thermal: 27.86\n",
      "FWHM_y BEC: 30.24, FWHM_y thermal: 26.64\n",
      "\n",
      "FWHM_x BEC: 9.40, FWHM_x thermal: 18.84\n",
      "FWHM_y BEC: 22.97, FWHM_y thermal: 22.27\n",
      "\n",
      "FWHM_x BEC: 9.60, FWHM_x thermal: 10.31\n",
      "FWHM_y BEC: 25.20, FWHM_y thermal: 16.63\n",
      "\n",
      "FWHM_x BEC: 9.39, FWHM_x thermal: 13.24\n",
      "FWHM_y BEC: 25.24, FWHM_y thermal: 19.14\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for i in range(0,shape[0]):\n",
    "    for j in range(0,shape[1]):\n",
    "        sx = result[i][j].best_values['sigmax_bec']\n",
    "        sy = result[i][j].best_values['sigmay_bec']\n",
    "        sx_th = result[i][j].best_values['sigmax_th']\n",
    "        sy_th = result[i][j].best_values['sigmay_th']\n",
    "        print(f'FWHM_x BEC: { sx*1.22:.2f}, FWHM_x thermal: { sx_th*1.93:.2f}')\n",
    "        print(f'FWHM_y BEC: { sy*1.22:.2f}, FWHM_y thermal: { sy_th*1.93:.2f}')\n",
    "        print('')"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-26T07:28:08.962590400Z",
     "start_time": "2023-07-26T07:28:08.903953900Z"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "outputs": [
    {
     "data": {
      "text/plain": "<lmfit.model.ModelResult at 0x1f08f1e7d60>",
      "text/html": "<h2> Model</h2> Model(density_profile_BEC_2d) <h2>Fit Statistics</h2><table><tr><td>fitting method</td><td>leastsq</td><td></td></tr><tr><td># function evals</td><td>660</td><td></td></tr><tr><td># data points</td><td>22500</td><td></td></tr><tr><td># variables</td><td>10</td><td></td></tr><tr><td>chi-square</td><td> 115.043851</td><td></td></tr><tr><td>reduced chi-square</td><td> 0.00511533</td><td></td></tr><tr><td>Akaike info crit.</td><td>-118689.037</td><td></td></tr><tr><td>Bayesian info crit.</td><td>-118608.825</td><td></td></tr><tr><td>R-squared</td><td> 0.36517195</td><td></td></tr></table><h2>Variables</h2><table><tr><th> name </th><th> value </th><th> standard error </th><th> relative error </th><th> initial value </th><th> min </th><th> max </th><th> vary </th><th> expression </th></tr><tr><td> amp_bec </td><td>  0.43797350 </td><td>  0.04403209 </td><td> (10.05%) </td><td> 0.6412865073706878 </td><td>  0.00000000 </td><td>  1.20759571 </td><td> True </td><td>  </td></tr><tr><td> amp_th </td><td>  0.33848251 </td><td>  0.05507434 </td><td> (16.27%) </td><td> 0.2748370745874376 </td><td>  0.00000000 </td><td>  1.20759571 </td><td> True </td><td>  </td></tr><tr><td> x0_bec </td><td>  76.9710005 </td><td>  0.14574325 </td><td> (0.19%) </td><td> 76.22058823529412 </td><td>  0.00000000 </td><td>  150.000000 </td><td> True </td><td>  </td></tr><tr><td> y0_bec </td><td>  58.6853106 </td><td>  0.32241760 </td><td> (0.55%) </td><td> 58.63970588235293 </td><td>  0.00000000 </td><td>  150.000000 </td><td> True </td><td>  </td></tr><tr><td> x0_th </td><td>  77.2286423 </td><td>  0.27255981 </td><td> (0.35%) </td><td> 76.22058823529412 </td><td>  0.00000000 </td><td>  150.000000 </td><td> True </td><td>  </td></tr><tr><td> y0_th </td><td>  58.8262204 </td><td>  0.55406928 </td><td> (0.94%) </td><td> 58.63970588235293 </td><td>  0.00000000 </td><td>  150.000000 </td><td> True </td><td>  </td></tr><tr><td> sigmax_bec </td><td>  8.19753077 </td><td>  0.18267507 </td><td> (2.23%) </td><td> 9.0 </td><td>  0.00000000 </td><td>  50.0000000 </td><td> True </td><td>  </td></tr><tr><td> sigmay_bec </td><td>  17.5376093 </td><td>  0.41011798 </td><td> (2.34%) </td><td> 18.0 </td><td>  0.00000000 </td><td>  50.0000000 </td><td> True </td><td>  </td></tr><tr><td> D_sigX </td><td>  0.86059710 </td><td>  0.72666984 </td><td> (84.44%) </td><td> 6.390000000000001 </td><td>  0.00000000 </td><td>  50.0000000 </td><td> True </td><td>  </td></tr><tr><td> D_sigY </td><td>  5.4950e-11 </td><td>  0.00362848 </td><td> (6603174582.12%) </td><td> 6.390000000000001 </td><td>  0.00000000 </td><td>  50.0000000 </td><td> True </td><td>  </td></tr><tr><td> sigmax_th </td><td>  5.62662875 </td><td>  0.32567851 </td><td> (5.79%) </td><td> 8.99802 </td><td>  0.00000000 </td><td>         inf </td><td> False </td><td> 0.632*sigmax_bec + 0.518*D_sigX </td></tr><tr><td> sigmay_th </td><td>  11.0837691 </td><td>  0.25919456 </td><td> (2.34%) </td><td> 14.68602 </td><td>  0.00000000 </td><td>         inf </td><td> False </td><td> 0.632*sigmay_bec + 0.518*D_sigY </td></tr></table><h2>Correlations (unreported correlations are < 0.100)</h2><table><tr><td>amp_bec</td><td>amp_th</td><td>-0.9816</td></tr><tr><td>y0_bec</td><td>y0_th</td><td>-0.8732</td></tr><tr><td>x0_bec</td><td>x0_th</td><td>-0.8519</td></tr><tr><td>sigmay_bec</td><td>D_sigY</td><td>+0.7342</td></tr><tr><td>amp_th</td><td>D_sigX</td><td>-0.6601</td></tr><tr><td>amp_bec</td><td>D_sigX</td><td>+0.6570</td></tr><tr><td>amp_th</td><td>D_sigY</td><td>+0.6209</td></tr><tr><td>amp_bec</td><td>D_sigY</td><td>-0.6140</td></tr><tr><td>sigmax_bec</td><td>D_sigX</td><td>-0.5632</td></tr><tr><td>sigmax_bec</td><td>D_sigY</td><td>-0.2750</td></tr><tr><td>sigmax_bec</td><td>sigmay_bec</td><td>-0.2644</td></tr><tr><td>D_sigX</td><td>D_sigY</td><td>-0.2165</td></tr><tr><td>sigmay_bec</td><td>D_sigX</td><td>+0.1629</td></tr><tr><td>amp_th</td><td>sigmax_bec</td><td>-0.1152</td></tr><tr><td>amp_bec</td><td>sigmay_bec</td><td>-0.1117</td></tr></table>"
     },
     "execution_count": 152,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for i in range(0, shape[0]):\n",
    "    for j in range(0, shape[1]):\n"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2023-07-25T14:17:00.810233Z",
     "start_time": "2023-07-25T14:17:00.767238Z"
    }
   }
  },
  {
   "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
}