{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8b6417c4",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# SCI with transverse motion with precomputed Bragg table"
   ]
  },
  {
   "cell_type": "raw",
   "id": "e4f8d615-db55-4c5d-9367-c0844af3734b",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/restructuredtext",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "We can use the precompute table we made in `this example`_ and the :py:meth:`mwave.simulation_utils.cloud_to_scrbi_ellipse_xy` method to simulate an ellipse with realistic beamsplitters.\n",
    "\n",
    ".. _`this example`: gen_precompute.ipynb\n",
    "\n",
    "First we import the needed libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "642c2d5f-1907-49d0-aeb9-e2257a025d64",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from mwave.precompute import load_fast_bragg_evaluator\n",
    "from mwave.simulation_utils import cloud_to_scrbi_ellipse_xy\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "raw",
   "id": "bfd115d4-0e62-4b4f-a9b5-7af2568a84a3",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/restructuredtext",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Now define the initial conditions, load the precompute table we made previously using :py:meth:`mwave.precompute.load_fast_bragg_evaluator`, and define the laser phase and effective Rabi frequency profiles in space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a285565b-975f-4716-a384-5f6a993bd8dc",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "np.random.seed(13703599)\n",
    "\n",
    "npoints = 100000\n",
    "x0 = np.random.randn(npoints)*0.75e-3\n",
    "y0 = np.random.randn(npoints)*0.75e-3\n",
    "z0 = np.random.randn(npoints)*0.75e-3\n",
    "v0 = np.random.randn(npoints)*0.35e-3\n",
    "vx = np.random.randn(npoints)*3.5e-3*1.8\n",
    "vy = np.random.randn(npoints)*3.5e-3*1.8\n",
    "T = 100e-3\n",
    "Tp = 10e-3\n",
    "phi_c = np.linspace(0,2*np.pi)\n",
    "phi_d = np.pi/4\n",
    "Omega0 = 32\n",
    "w0 = 10e-3\n",
    "\n",
    "n_init = 0\n",
    "n_bragg = 5\n",
    "N_bloch = 10\n",
    "\n",
    "bragglookup = load_fast_bragg_evaluator('sig0.247.h5', n_init, n_bragg, N_bloch)\n",
    "\n",
    "def deltalookup(v, n_bragg):\n",
    "    return 4*n_bragg + 4*(v/3.5e-3) # The modification to delta is 4 times the velocity defined in units of recoil velocities\n",
    "\n",
    "def omegalookup(x, y, z):\n",
    "    return Omega0*np.exp(-2*(x**2 + y**2)/(w0**2))\n",
    "\n",
    "# (l)ocal (p)hase (lookup)\n",
    "def lplookup(x, y, z):\n",
    "    wavelen = 852e-9\n",
    "    zR = np.pi*w0**2/wavelen\n",
    "    kk = 2*np.pi/wavelen\n",
    "    return kk*(x**2 + y**2)/(2*zR)\n",
    "\n",
    "x, y = cloud_to_scrbi_ellipse_xy(x0, y0, z0, vx, vy, v0, T, Tp, n_bragg, N_bloch, np.exp(1j*phi_c), np.exp(1j*phi_d), bragglookup, omegalookup, lplookup, deltalookup)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9462b5a2-a6a6-40f3-ae2c-1fb721071c98",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Now plot and fit the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6d39d0fd-68cd-4dc1-954c-0b6c336a360d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAckAAAGiCAYAAACMDD3oAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAO8dJREFUeJzt3Xt8E2W+P/BPWmgKYlKwtGm1cj9Al0IRbCgHxZVKK6zCWc4uRZTLwXJEQbEo0D3SCrhbbodlVVaU5eZPEcUXFy9YrUV01dhioSpQeNlulYukQGsTKFJo8vz+6Ens0EybtJk0yXzer1deu5k8M3lmLPnOPN/nohFCCBAREVETIe1dASIiIn/FIElERCSDQZKIiEgGgyQREZEMBkkiIiIZDJJEREQyGCSJiIhkMEgSERHJYJAkIiKSwSBJREQkQ9Eg+dlnn+G+++5DbGwsNBoN9uzZ0+I+Bw4cwG233QatVou+ffti69atTcqsX78ePXv2RHh4OIxGI4qKirxfeSIiUj1Fg2RtbS2GDBmC9evXu1W+oqIC48ePx29/+1uUlJRg/vz5ePjhh/Hhhx86y7z55pvIzMxETk4ODh06hCFDhiA1NRXnzp1T6jSIiEilNL6a4Fyj0WD37t2YOHGibJlFixbh/fffx5EjR5zb0tPTUVNTg7y8PACA0WjE7bffjhdffBEAYLfbERcXh3nz5mHx4sWKngMREalLh/auQGMmkwkpKSmSbampqZg/fz4A4OrVqyguLkZWVpbz85CQEKSkpMBkMsket66uDnV1dc73drsd1dXVuOmmm6DRaLx7EkREpDghBC5evIjY2FiEhCjXKOpXQdJsNiM6OlqyLTo6GlarFb/88gt+/vln2Gw2l2WOHz8ue9zc3FwsXbpUkToTEVH7OXXqFG655RbFju9XQVIpWVlZyMzMdL63WCy49dZbcerUKeh0unasGRERtYbVakVcXBxuvPFGRb/Hr4KkwWBAZWWlZFtlZSV0Oh06deqE0NBQhIaGuixjMBhkj6vVaqHVapts1+l0DJJERAFM6ZSZX42TTE5ORkFBgWRbfn4+kpOTAQBhYWEYNmyYpIzdbkdBQYGzDBERkbcoGiQvXbqEkpISlJSUAGgY4lFSUoKTJ08CaGgGnTZtmrP8I488gn/9619YuHAhjh8/jr///e9466238OSTTzrLZGZmYuPGjdi2bRtKS0sxZ84c1NbWYubMmUqeChERqZCiza1ff/01fvvb3zrfO/KC06dPx9atW3H27FlnwASAXr164f3338eTTz6Jv/3tb7jlllvwj3/8A6mpqc4ykydPxvnz55GdnQ2z2YzExETk5eU16cxDRETUVj4bJ+lPrFYr9Ho9LBYLc5JERAHIV7/jfpWTJCIi8icMkkRERDIYJImIiGQwSBIREclgkCQiIpLBIElERCSDQZKIiEgGgyQREZEMBkkiIiIZfrUKCJG/sNkFiiqqce7iFUTdGI6kXt0QGtL21QaUOi4RKYNBkug6eUfOYum7x3DWcsW5LUYfjpz74pE2KMbvjktEymFzKwUFm13AVF6FvSVnYCqvgs3euimJ846cxZzXDkkCGQCYLVcw57VDyDty1q+OS0TK4pMkBTxvPaHZ7AJL3z0GV+FVANAAWPruMdwTb/CoiVSp4zY+PptwiZTBIEkBzfGEdn0AcjyhvfTgbW4HyqKK6iZPeo0JAGctV1BUUY3kPje5XUeljguwCZdIaWxupYDV0hMa0PCE5m7T67mL8oGsNeWUPi6bcImUxyBJAcuTJzR3RN0Y7tVySh7X2zcIROQagyS1C290tPH2E1pSr26I0YdDLpunQUNTZlKvbu5VUMHjevsGgYhcY06SfM5beTRvP6GFhmiQc1885rx2CBpA8pTmCHA598V73ClGieMq1YTrwM5ARA34JEk+5c08mhJPaGmDYvDSg7fBoJcGVoM+3KNOQEofV6mmYaDhv9GolfsxZeNXeGJHCaZs/AqjVu5njpNUSSOEUF3Swmq1Qq/Xw2KxQKfTtXd1VMNmFxi1cr9sM6EGDUHj80V3u/3U4gi6gOsntNYGNn+fccdxLc2WKy7zkq25loB8b+G2Xk8ib/PV7zifJMlnlMijKfXkFxqiQXKfmzAh8WYk97nJa02N3jquowkXQJMn6dY24bIzEFFTzEmSzyiVR0sbFIN74g2qy6E5bhCuz+8aWjlOUsnxnESBikGSfEbJPJrjCU1tvHmDoHRnIKJAxCBJbmtrPs3R0aalPJqnQyzUzls3CErexBAFKgZJcos3hm0oNcSCvIM3MURNseMOtcibwzaU6mhDbadEZ6DGvLVSC5EvcQgIh4A0S4lhG47jqq2jTaBQYtJ0TsRO3uar33EGSQbJZpnKqzBl41ctlnsjY4QqO84EK2/exHDsJSnBV7/jzElSs9jjUZ281RlI6bU0iZTGnCQ1iz0eqS04ETsFOgZJapZSK2OQOrAlggIdgyQ1S+kejxTc2BJBgY5BUiXa0v2ewzaotdgSQYGOHXdUwBvd79U6Pyq1DSeQoEDHISBBPgSE3e/JH3h7nCTH2VJQLZW1fv169OzZE+Hh4TAajSgqKpIte9ddd0Gj0TR5jR8/3llmxowZTT5PS0vzxakEFC59RP4ibVAMPl90N97IGIG/pSfijYwR+HzR3a0KkFwUmnxJ8ebWN998E5mZmdiwYQOMRiPWrVuH1NRUnDhxAlFRUU3K79q1C1evXnW+r6qqwpAhQ/CHP/xBUi4tLQ1btmxxvtdqtcqdRIDi0kfkT7wx9lKuZcQxRSJbRsjbFH+SXLt2LTIyMjBz5kzEx8djw4YN6Ny5MzZv3uyyfLdu3WAwGJyv/Px8dO7cuUmQ1Gq1knJdu3ZV+lQCDrvfUzBhywi1B0WD5NWrV1FcXIyUlJRfvzAkBCkpKTCZTG4dY9OmTUhPT8cNN9wg2X7gwAFERUWhf//+mDNnDqqqqmSPUVdXB6vVKnmpAbvfUzDhxATUHhQNkhcuXIDNZkN0dLRke3R0NMxmc4v7FxUV4ciRI3j44Ycl29PS0vDqq6+ioKAAK1euxKeffop7770XNpvN5XFyc3Oh1+udr7i4uNafVABh93sKJmwZofbg1+MkN23ahISEBCQlJUm2p6en4/7770dCQgImTpyI9957DwcPHsSBAwdcHicrKwsWi8X5OnXqlA9q3/44EQAFE7aMUHtQNEhGRkYiNDQUlZWVku2VlZUwGAzN7ltbW4sdO3Zg1qxZLX5P7969ERkZibKyMpefa7Va6HQ6yUstOBEABQu2jFB7ULR3a1hYGIYNG4aCggJMnDgRAGC321FQUIC5c+c2u+/OnTtRV1eHBx98sMXvOX36NKqqqhATwx98VzgRAAUDTkxA7UHx5tbMzExs3LgR27ZtQ2lpKebMmYPa2lrMnDkTADBt2jRkZWU12W/Tpk2YOHEibrpJ2mX80qVLePrpp/HVV1/hhx9+QEFBASZMmIC+ffsiNTVV6dMJWI7u9xMSb0Zyn5v4Q0IBiS0j5GuKj5OcPHkyzp8/j+zsbJjNZiQmJiIvL8/ZmefkyZMICZHG6hMnTuDzzz/HRx991OR4oaGh+Pbbb7Ft2zbU1NQgNjYWY8eOxfLly4N+rCRnGSFiywj5FqelC5D8pLen9SIiCmRBNS0dtY1jlpHrx4g5ZhnhdFxErdeWFXIo+HEVED/X0iwjGjTMMnJPvIHNTUQeYgsNtYRPkn6Os4wQKYMtNOQOBkk/x1lGiLyP88CSuxgk/RxnGSHyPrbQkLsYJP0cZxkh8j620JC7GCT9HOdfJfI+ttCQuxgkAwBnGSHyLrbQkLs4BCRAcJYRIu/hPLDkLs64EyAz7hCR93GcZODy1e84nySJSLXYQkMtYZBsJ5ysnMg/OFbIIXKFQbIdsImHiCgwsHerj3EqLCKiwMEnSR/iZOVEwelqvR3/z/QDfqy+jB7dOuOh5J4I68BnkGDAIOlDnkyFxRwJUWDI3XcMG/9ZgcbTvP55Xyky7uiFrHHx7Vcx8goGSR/iVFhEwSV33zG8/FlFk+12Aed2BsrAxvYAH+JUWETB42q9HRv/2TRANrbxnxW4Wm/3UY1ICQySPsSpsIiCx/8z/YCWVtKyi4ZyFLgYJH2Ik5UTBY8fqy97tRz5JwZJH+Nk5UTBoUe3zl4tR/6Jc7e209ytnHGHKLBdrbdjwJIPmm1yDdEAx5ffy+EgCuDcrUGOU2ERBbawDiHIuKOXy96tDhl39GKADHAMkkREreQY3nH9OMkQDThOMkiwuZVLZRFRG3HGHd9jcysRUYAI6xCCWXf0bu9qkAIYJL2EHXGIiIIPg6QXcOkrIqLgxEbzNuLSV0REwYtBsg1aWvoKaFj6ytbS3FVEpFo2u4CpvAp7S87AVF7F3ws/w+bWNuDSV0TUFkzV+D8+SbYBl74iotZiqiYwMEi2AZe+IqLWYKomcDBItgGXviKi1vAkVUPti0GyDbj0FRG1BlM1gcMnQXL9+vXo2bMnwsPDYTQaUVRUJFt269at0Gg0kld4uLS5UgiB7OxsxMTEoFOnTkhJScH333+v9Gm4xKWviMhTTNUEDsV7t7755pvIzMzEhg0bYDQasW7dOqSmpuLEiROIiopyuY9Op8OJEyec7zUa6ZPYqlWr8Pzzz2Pbtm3o1asXlixZgtTUVBw7dqxJQPWFtEExuCfewBl3iMgtjlSN2XLFZV5Sg4YbbaZq2p/iT5Jr165FRkYGZs6cifj4eGzYsAGdO3fG5s2bZffRaDQwGAzOV3R0tPMzIQTWrVuHZ555BhMmTMDgwYPx6quv4qeffsKePXuUPh1ZjqWvJiTejOQ+NzFAEpEspmoCh6JB8urVqyguLkZKSsqvXxgSgpSUFJhMJtn9Ll26hB49eiAuLg4TJkzA0aNHnZ9VVFTAbDZLjqnX62E0GmWPWVdXB6vVKnkREbUnpmoCg6LNrRcuXIDNZpM8CQJAdHQ0jh8/7nKf/v37Y/PmzRg8eDAsFgvWrFmDkSNH4ujRo7jllltgNpudx7j+mI7Prpebm4ulS5d64YyIiLyHqRr/53cz7iQnJyM5Odn5fuTIkRg4cCBefvllLF++vFXHzMrKQmZmpvO91WpFXFxcm+tKRNRWjlQN+SdFm1sjIyMRGhqKyspKyfbKykoYDAa3jtGxY0cMHToUZWVlAODcz5NjarVa6HQ6yYuIiKgligbJsLAwDBs2DAUFBc5tdrsdBQUFkqfF5thsNnz33XeIiWlon+/VqxcMBoPkmFarFYWFhW4fk4iIyB2KN7dmZmZi+vTpGD58OJKSkrBu3TrU1tZi5syZAIBp06bh5ptvRm5uLgBg2bJlGDFiBPr27YuamhqsXr0aP/74Ix5++GEADT1f58+fj+eeew79+vVzDgGJjY3FxIkTFTsPLqpMRKQ+igfJyZMn4/z588jOzobZbEZiYiLy8vKcHW9OnjyJkJBfH2h//vlnZGRkwGw2o2vXrhg2bBi+/PJLxMfHO8ssXLgQtbW1mD17NmpqajBq1Cjk5eUpNkaSM/UTEamTRgihuhl0rVYr9Ho9LBZLi/lJx0z9118kxzMku2oTEfmeJ7/jbcG5W5vBmfqJyF9xsWbf8LshIP6EiyoTkT9iCsh3+CTZDM7UT0T+hos1+xaDZDM4Uz8R+ROmgHyPQbIZXFSZiPwJF2v2PQbJZnCmfiLyJ0wB+R6DZAs4Uz8R+QumgHyPvVvdwJn6icgfcLFm32OQdBNn6iei9uZIAc157RA0gCRQMgWkDDa3EhEFEKaAfItPkkREAYYpIN9hkCQiCkBMAfkGm1uJiIhkMEgSERHJYJAkIiKSwZzkdWx2wWQ4EREBYJCU4PIzRETUGJtb/w+XnyEiousxSILLzxBRcLPZBUzlVdhbcgam8ir+lnmAza3wbPkZjksiokDCNFLb8EkSXH6GiIIT00htxyAJLj9DRMGHaSTvYJDEr8vPyA300KCheYLLzxBRoPAkjUTyGCTx6/IzAJoESi4/Q0SBiGkk72CQ/D9cfoaIggnTSN7B3q2NcPkZIgoWjjSS2XLFZV5Sg4aHAKaRmscgeR0uP0NEwcCRRprz2iFoAEmgZBrJfWxuJSIKUkwjtR2fJImIghjTSG3DIElEFOSYRmo9NrcSERHJYJAkIiKSwSBJREQkg0GSiIhIhuo77tjsgr2+iIjIJVUHyfxjZqz55CDXWSMiIpd80ty6fv169OzZE+Hh4TAajSgqKpItu3HjRtxxxx3o2rUrunbtipSUlCblZ8yYAY1GI3mlpaV5XK/MN7/hOmtERCRL8SD55ptvIjMzEzk5OTh06BCGDBmC1NRUnDt3zmX5AwcOYMqUKfjkk09gMpkQFxeHsWPH4syZM5JyaWlpOHv2rPP1xhtveFw3rrNGRNSQdjKVV2FvyRmYyqv429eIRgih6NUwGo24/fbb8eKLLwIA7HY74uLiMG/ePCxevLjF/W02G7p27YoXX3wR06ZNA9DwJFlTU4M9e/a4VYe6ujrU1dU531utVsTFxSFu/lsI0XaW3e+NjBEcgEtEQS3vyFksffdYwKWdrFYr9Ho9LBYLdDqdYt+j6JPk1atXUVxcjJSUlF+/MCQEKSkpMJlMbh3j8uXLuHbtGrp1k85Uf+DAAURFRaF///6YM2cOqqqqZI+Rm5sLvV7vfMXFxbn13VxnjYiCWd6Rs5jz2iGmnZqhaJC8cOECbDYboqOjJdujo6NhNpvdOsaiRYsQGxsrCbRpaWl49dVXUVBQgJUrV+LTTz/FvffeC5vN5vIYWVlZsFgsztepU6fc+m6us0ZEwcpmF1j67jGmnVrg171bV6xYgR07duDAgQMID/81YKWnpzv/f0JCAgYPHow+ffrgwIEDGDNmTJPjaLVaaLXaJtvlBnpwnTUiCnZFFdVNniAbEwDOWq6gqKJa1WknRZ8kIyMjERoaisrKSsn2yspKGAyGZvdds2YNVqxYgY8++giDBw9utmzv3r0RGRmJsrIyj+t4faDkOmtEpAbuppPUnnZSNEiGhYVh2LBhKCgocG6z2+0oKChAcnKy7H6rVq3C8uXLkZeXh+HDh7f4PadPn0ZVVRViYjxLMq+dPITrrBGRKrmbTlJ72knx5tbMzExMnz4dw4cPR1JSEtatW4fa2lrMnDkTADBt2jTcfPPNyM3NBQCsXLkS2dnZ2L59O3r27OnMXXbp0gVdunTBpUuXsHTpUkyaNAkGgwHl5eVYuHAh+vbti9TUVI/qdk+8AROT+nHGHSJSnaRe3RCjD4fZcsVlXpJppwaKB8nJkyfj/PnzyM7OhtlsRmJiIvLy8pydeU6ePImQkF8faF966SVcvXoV//mf/yk5Tk5ODp599lmEhobi22+/xbZt21BTU4PY2FiMHTsWy5cvd5l3bAnXWSMiNQoN0SDnvnjMee0QNJCOG2fa6VeKj5P0R74aX0NE5O84TrJ5ft27lYiIlJU2KAb3xBuYdpLBIElEpHJMO8njepJEREQyGCSJiIhkMEgSERHJYJAkIiKSoeqOO4411Niji4iIXFF1kBz7109xvi7U+T4QxgYREZHvqLq5tdJaJ3nPNdSIiKgxVQfJ63ENNSIiKUdaam/JGZjKq1T326jq5lZXuIYaEVGDQJ2yzpv4JClD7WuoEZG65R05izmvHWqyMLPa0lIMkjLUvoYaEamXzS6w9N1jLpfQUltaikHyOho0NCeofQ01IlKvoorqJk+QjTVOSwU7VQfJ60dEcg01IiL3001qSEupOkhG6aSLNBv04XjpwdtUk5AmInLF3XSTGtJSqu7d+tGTo3G8qp4z7hARNZLUqxti9OEwW664zEtq0PBQoYa0lKqDJNdQIyJqKjREg5z74jHntUPQAJJAqba0lKqbW4mIyLW0QTF46cHbYNBLm1TVlpZS9ZMkERHJSxsUg3viDSiqqFZtWopBkoiIZKk9LcXmViIiIhkMkkRERDIYJImIiGQwSBIREclgkCQiIpLB3q1omPFezV2ciYjINdUHSS4qSkREclTd3Jp/zMxFRYmISJaqg+SKD45zUVEiIjfZ7AKm8irsLTkDU3mVKn4fVd3cWmmtQ4i2s8vPGi8qqubZJoiIAPWmplT9JOkONSwqSkTUnLwjZ1WbmmKQbIEaFhUlIpJjswssffeYalNTqg6S0Tot5AZ6aNDQlKCGRUWJiOQUVVQ3eYJsrHFqKhipOkguvncAADQJlGpbVJSISI67KadgTU35JEiuX78ePXv2RHh4OIxGI4qKipotv3PnTgwYMADh4eFISEjAvn37JJ8LIZCdnY2YmBh06tQJKSkp+P777z2u1z3xBi4qSkTUDHdTTsGamlI8SL755pvIzMxETk4ODh06hCFDhiA1NRXnzp1zWf7LL7/ElClTMGvWLBw+fBgTJ07ExIkTceTIEWeZVatW4fnnn8eGDRtQWFiIG264AampqbhyxfM7mbRBMfh80d14I2ME/paeiDcyRuDzRXczQBIRAUjq1Q0x+nDVpqY0QghFs61GoxG33347XnzxRQCA3W5HXFwc5s2bh8WLFzcpP3nyZNTW1uK9995zbhsxYgQSExOxYcMGCCEQGxuLBQsW4KmnngIAWCwWREdHY+vWrUhPT29yzLq6OtTV1TnfW61WxMXFwWKxQKfTefuUiYiCiqN3KwBJBx5H4GyPljer1Qq9Xq/477iiT5JXr15FcXExUlJSfv3CkBCkpKTAZDK53MdkMknKA0BqaqqzfEVFBcxms6SMXq+H0WiUPWZubi70er3zFRcX19ZTIyJSjbRBMapNTSk6mcCFCxdgs9kQHR0t2R4dHY3jx4+73MdsNrssbzabnZ87tsmVuV5WVhYyMzOd7x1PkkRE5J60QTG4J96gusUgVDHjjlarhVarbe9qEBEFtNAQjepmIFO0uTUyMhKhoaGorKyUbK+srITBYHC5j8FgaLa84389OSYREVFrKBokw8LCMGzYMBQUFDi32e12FBQUIDk52eU+ycnJkvIAkJ+f7yzfq1cvGAwGSRmr1YrCwkLZYxIREbWG4s2tmZmZmD59OoYPH46kpCSsW7cOtbW1mDlzJgBg2rRpuPnmm5GbmwsAeOKJJzB69Gj87//+L8aPH48dO3bg66+/xiuvvAIA0Gg0mD9/Pp577jn069cPvXr1wpIlSxAbG4uJEycqfTpERKQiigfJyZMn4/z588jOzobZbEZiYiLy8vKcHW9OnjyJkJBfH2hHjhyJ7du345lnnsGf/vQn9OvXD3v27MGgQYOcZRYuXIja2lrMnj0bNTU1GDVqFPLy8hAeHpyDWYmIqH0oPk7SH7kzvsZmF6rrxUVEFCh8NU5SFb1bPaXWddOIiEhK1ROcu6LmddOIiEiKQbIRta+bRkTUEptdwFRehb0lZ2Aqrwr630M2tzbiybppahtQS0SkxlQUnyQbUfu6aUREctSaimKQbETt66YREbmi5lQUg2Qjal83jYjIFU9SUcGGQbKR0BANcu6LB4AmgdLxPue+eI6XJCJVUXMqikHyOmpeN42IyBU1p6LYu9UFta6bRkTkiiMVZbZccZmX1KDhQSIYU1EMkjLUuG4aEZErjlTUnNcOQQNIAmWwp6LY3EpERC1SayqKT5JEROQWNaaiGCSJiMhtaktFsbmViIhIBoMkERGRDAZJIiIiGQySREREMthxpxVsdqGq3l1ERGrFIOkhNa6nRkSkVmxu9YBa11MjIlIrBkk3qXk9NSIim13AVF6FvSVnYCqvUs1vHZtb3eTJempqGmhLRMFPzWkmPkm6Sc3rqRGReqk9zcQg6SY1r6dGROrENBODpNsc66nJDfTQoKH5IRjXUyMidfIkzRSsGCTd5FhPDUCTQBns66kRkToxzcQg6RG1rqdGROrENBN7t3pMjeupEZE6OdJMZssVl3lJDRoeEoI5zcQg2QpqW0+NiNTJkWaa89ohaABJoFRLmonNrUREJEvtaSY+SRIRUbPUnGZikCQiohapNc3E5lYiIiIZDJJEREQyFA2S1dXVmDp1KnQ6HSIiIjBr1ixcunSp2fLz5s1D//790alTJ9x66614/PHHYbFYJOU0Gk2T144dO5Q8FSIiUiFFc5JTp07F2bNnkZ+fj2vXrmHmzJmYPXs2tm/f7rL8Tz/9hJ9++glr1qxBfHw8fvzxRzzyyCP46aef8Pbbb0vKbtmyBWlpac73ERERSp4KERGpkEYIocjMtKWlpYiPj8fBgwcxfPhwAEBeXh7GjRuH06dPIzY21q3j7Ny5Ew8++CBqa2vRoUNDTNdoNNi9ezcmTpzYqrpZrVbo9XpYLBbodLpWHcMdNrtQZW8wIiKl+ep3XLEnSZPJhIiICGeABICUlBSEhISgsLAQ//Ef/+HWcRwXwBEgHR577DE8/PDD6N27Nx555BHMnDkTGo3rAFRXV4e6ujrne6vV2ooz8oya118jIgoWiuUkzWYzoqKiJNs6dOiAbt26wWw2u3WMCxcuYPny5Zg9e7Zk+7Jly/DWW28hPz8fkyZNwqOPPooXXnhB9ji5ubnQ6/XOV1xcnOcn5AG1r79GRIHNZhcwlVdhb8kZmMqrgnoprJZ4/CS5ePFirFy5stkypaWlra6Qg9Vqxfjx4xEfH49nn31W8tmSJUuc/3/o0KGora3F6tWr8fjjj7s8VlZWFjIzMyXHVipQtrT+mgYN66/dE29g0ysR+R22gkl5HCQXLFiAGTNmNFumd+/eMBgMOHfunGR7fX09qqurYTAYmt3/4sWLSEtLw4033ojdu3ejY8eOzZY3Go1Yvnw56urqoNVqm3yu1WpdbleCJ+uvqXFgLhH5L0cr2PU3+Y5WMDVMQ3c9j4Nk9+7d0b179xbLJScno6amBsXFxRg2bBgAYP/+/bDb7TAajbL7Wa1WpKamQqvV4p133kF4eMtLsJSUlKBr164+C4TN4fprRBSI2ArmmmI5yYEDByItLQ0ZGRkoKirCF198gblz5yI9Pd3Zs/XMmTMYMGAAioqKADQEyLFjx6K2thabNm2C1WqF2WyG2WyGzWYDALz77rv4xz/+gSNHjqCsrAwvvfQS/vKXv2DevHlKnYpHuP4aEQUiT1rB1ETRcZKvv/465s6dizFjxiAkJASTJk3C888/7/z82rVrOHHiBC5fvgwAOHToEAoLCwEAffv2lRyroqICPXv2RMeOHbF+/Xo8+eSTEEKgb9++WLt2LTIyMpQ8Fbdx/TUiCkRsBXNNsXGS/kzp8TWOdn3A9fpramzXJyL/ZiqvwpSNX7VY7o2MEX7Rn8JX4yQ5d6sC1L7+GhEFHkcrmFy2UYOGXq5qawXjUlkKUfP6a0QUeEJDNMi5Lx5zXjsEDVy3guXcF6+63zA2tyr4mE5EFGgCZZxkwE9LR0REgYetYFIMkkREJBEaovGLzjn+gB13iIiIZDBIEhERyWCQJCIiksGcpB/hIs1ERP6FQdJPBEq3ayIiNWFzqx/gIs1ERP6JQbKdtbQ8DdCwPI2aVwYnIu+z2QVM5VXYW3IGpvIq/sbIYHNrO+MizUTka0zvuI9Pku2My9MQkS8xveMZBsl2xkWaichXmN7xHINkO+PyNETkK56kd6gBg2Q7cyxPA6BJoFTz8jRE5H1M73iOQdIPcJFmIvIFpnc8x96tfoLL0xCR0hzpHbPlisu8pAYNN+dM7/yKQdKPcHkaIlKSI70z57VD0ACSQMn0jmtsbiUiUhGmdzzDJ0kiIpVhesd9DJJERCrE9I572NxKREQkg0+SQYJrURIReR+DZBDgZMVERMpgc2uA42TFRETKYZAMYJysmIiux3UivYvNrQGMa1ESUWNMvXgfnyQDGCcrJiIHpl6UwSAZwDhZMREBTL0oiUEygHEtSiICuE6kkhgkAxjXoiQigKkXJTFIBjhOVkxETL0oh71bgwAnKyZSN64TqRxFnySrq6sxdepU6HQ6REREYNasWbh06VKz+9x1113QaDSS1yOPPCIpc/LkSYwfPx6dO3dGVFQUnn76adTX1yt5Kn7PMVnxhMSbkdznJgZIIhVh6kU5igbJqVOn4ujRo8jPz8d7772Hzz77DLNnz25xv4yMDJw9e9b5WrVqlfMzm82G8ePH4+rVq/jyyy+xbds2bN26FdnZ2UqeChGRX2PqRRkaIYQifYJLS0sRHx+PgwcPYvjw4QCAvLw8jBs3DqdPn0ZsbKzL/e666y4kJiZi3bp1Lj//4IMP8Lvf/Q4//fQToqOjAQAbNmzAokWLcP78eYSFhbVYN6vVCr1eD4vFAp1O17oTJCLyQ2pZ7MBXv+OKPUmaTCZEREQ4AyQApKSkICQkBIWFhc3u+/rrryMyMhKDBg1CVlYWLl++LDluQkKCM0ACQGpqKqxWK44ePeryeHV1dbBarZIXEVEwYurFuxTruGM2mxEVFSX9sg4d0K1bN5jNZtn9HnjgAfTo0QOxsbH49ttvsWjRIpw4cQK7du1yHrdxgATgfC933NzcXCxdurQtpxP01HL3SUTkCY+D5OLFi7Fy5cpmy5SWlra6Qo1zlgkJCYiJicGYMWNQXl6OPn36tOqYWVlZyMzMdL63Wq2Ii4trdR2DDed7JCJyzeMguWDBAsyYMaPZMr1794bBYMC5c+ck2+vr61FdXQ2DweD29xmNRgBAWVkZ+vTpA4PBgKKiIkmZyspKAJA9rlarhVardfs71cQx3+P1iWnHfI9M+BO1L7bytC+Pg2T37t3RvXv3FsslJyejpqYGxcXFGDZsGABg//79sNvtzsDnjpKSEgBATEyM87h//vOfce7cOWdzbn5+PnQ6HeLj4z08G3Vrab5HDRrme7wn3sB/lETtgK087U+xjjsDBw5EWloaMjIyUFRUhC+++AJz585Fenq6s2frmTNnMGDAAOeTYXl5OZYvX47i4mL88MMPeOeddzBt2jTceeedGDx4MABg7NixiI+Px0MPPYRvvvkGH374IZ555hk89thjfFr0EOd7JPJfXNXDPyg6TvL111/HgAEDMGbMGIwbNw6jRo3CK6+84vz82rVrOHHihLP3alhYGD7++GOMHTsWAwYMwIIFCzBp0iS8++67zn1CQ0Px3nvvITQ0FMnJyXjwwQcxbdo0LFu2TMlTCUqc75HIP3FVD/+h6LR03bp1w/bt22U/79mzJxoP04yLi8Onn37a4nF79OiBffv2eaWOasb5Hon8ExdU9x+c4FzFuNQWkX9iK4//YJBUMc73SOSf2MrjPxgkVY7zPRL5H7by+A8ulUVcaovIzzhaeea8dggaQNKBh608vqXYBOf+jBOcE1Eg4DhJeb76HeeTJBGRn2IrT/tjkCSv4hRaRN7lWNWD2geDJHkNm4aIXOPNY+BikCSv4ETpRK7x5jGwcQgItRmn0CJyjfOvBj4GSWozTpRO1BRvHoMDgyS1GafQImqKN4/BgUGS2oxTaBE1xZvH4MAgSW3GKbSImuLNY3BgkKQ240TpRE3x5jE4MEiSV3CidCIp3jwGB87dyrlbvYqDpomkOE5SGb76HWeQZJD0Wwy4FCz4t+x9nOCcVI133+RP2hrkOP9q4GKQJL/DKe7In/CGTd3YcYf8CmcpIX/CaeWIQZL8CmcpIX/BGzYCGCTJz3CWEvIXvGEjgEGS/AxnKSF/wRs2Ahgkyc9wlhLyF7xhI4BBkvyM0rOU2OwCpvIq7C05A1N5FfNJJIs3bARwCAj5IccUd9d3uze0sds9u/KrT1vGNzpu2Oa8dggaQNKBh9PKqQdn3OGMO37Lm7OUyI29dByNYy+Dj7duinhz5Z84LZ2CGCTVxWYXGLVyv2xPRQ0anlI/X3Q3nwqChLdvijitnP/htHREXuJJV35OHRb4WhrfqEHD+MZ74g0eNb3yb0Od2HGHgh678qsLxzeSN/FJkoKe0l352RTnX3hTRN7EIElBz9GV32y54rIJzpGTbE1Xfnbq8L623nRwfCN5E4MkBT2luvJztRLv88ZNh5I3RaQ+zEmSKjjGXhr00qcHgz68VcGMk197n7dW3FB6QgpSF0WDZHV1NaZOnQqdToeIiAjMmjULly5dki3/ww8/QKPRuHzt3LnTWc7V5zt27FDyVCgIpA2KweeL7sYbGSPwt/REvJExAp8vurtVT3tKdw5R28xA3r7p8PZNEamXos2tU6dOxdmzZ5Gfn49r165h5syZmD17NrZv3+6yfFxcHM6eld4tvvLKK1i9ejXuvfdeyfYtW7YgLS3N+T4iIsLr9afg462u/Ep2DgmkPKe3Oi0pMUwnbVAM7ok3sFMVtYliQbK0tBR5eXk4ePAghg8fDgB44YUXMG7cOKxZswaxsbFN9gkNDYXBYJBs2717N/74xz+iS5cuku0RERFNyhL5ilKdQ5TMc3q7F643g7lSNx0c30htpVhzq8lkQkREhDNAAkBKSgpCQkJQWFjo1jGKi4tRUlKCWbNmNfnsscceQ2RkJJKSkrB582Y0N3FQXV0drFar5EXUFkpMfq1knjPvyFmMWrkfUzZ+hSd2lGDKxq8wauV+t/N8ro7njfyhA3ukkr9SLEiazWZERUVJtnXo0AHdunWD2Wx26xibNm3CwIEDMXLkSMn2ZcuW4a233kJ+fj4mTZqERx99FC+88ILscXJzc6HX652vuLg4z0+IqBElOocolef0dkBTIphzxQ3yVx4HycWLF8t2rnG8jh8/3uaK/fLLL9i+fbvLp8glS5bg3//93zF06FAsWrQICxcuxOrVq2WPlZWVBYvF4nydOnWqzfUj8nbnECWaHJUIaEoEc/ZIJX/lcU5ywYIFmDFjRrNlevfuDYPBgHPnzkm219fXo7q62q1c4ttvv43Lly9j2rRpLZY1Go1Yvnw56urqoNVqm3yu1WpdbidqK292DlGiyVGJDjFK5Q+VWiKNqC08DpLdu3dH9+7dWyyXnJyMmpoaFBcXY9iwYQCA/fv3w263w2g0trj/pk2bcP/997v1XSUlJejatSsDIbULb3UOUWIQvBIBTcn8IXukkr9RrHfrwIEDkZaWhoyMDGzYsAHXrl3D3LlzkZ6e7uzZeubMGYwZMwavvvoqkpKSnPuWlZXhs88+w759+5oc991330VlZSVGjBiB8PBw5Ofn4y9/+QueeuoppU6FyCeUmBlIiYCm9Iw27JFK/kTRyQRef/11DBgwAGPGjMG4ceMwatQovPLKK87Pr127hhMnTuDy5cuS/TZv3oxbbrkFY8eObXLMjh07Yv369UhOTkZiYiJefvllrF27Fjk5OUqeCpFPeDvPqUSHGOYPSU246DIXXSY/5M0xjY7erYDrp9PWjr0MpEkPKPj46necQZJBklRAqYDGZcKovTBIKohBktSIAY2Cia9+x7lUFpFKsEMMkee4VBYREZEMBkkiIiIZDJJEREQyGCSJiIhkMEgSERHJYJAkIiKSwSBJREQkg0GSiIhIBoMkERGRDAZJIiIiGQySREREMhgkiYiIZDBIEhERyWCQJCIiksEgSUREJINBkoiISAaDJBERkQwGSSIiIhkMkkRERDIYJImIiGQwSBIREclgkCQiIpLBIElERCSDQZKIiEgGgyQREZEMBkkiIiIZDJJEREQyGCSJiIhkMEgSERHJYJAkIiKSwSBJREQkg0GSiIhIBoMkERGRDMWC5J///GeMHDkSnTt3RkREhFv7CCGQnZ2NmJgYdOrUCSkpKfj+++8lZaqrqzF16lTodDpERERg1qxZuHTpkgJnQEREaqdYkLx69Sr+8Ic/YM6cOW7vs2rVKjz//PPYsGEDCgsLccMNNyA1NRVXrlxxlpk6dSqOHj2K/Px8vPfee/jss88we/ZsJU6BiIhUTiOEEEp+wdatWzF//nzU1NQ0W04IgdjYWCxYsABPPfUUAMBisSA6Ohpbt25Feno6SktLER8fj4MHD2L48OEAgLy8PIwbNw6nT59GbGysy2PX1dWhrq7O+d5iseDWW2/FqVOnoNPpvHOiRETkM1arFXFxcaipqYFer1fui4TCtmzZIvR6fYvlysvLBQBx+PBhyfY777xTPP7440IIITZt2iQiIiIkn1+7dk2EhoaKXbt2yR47JydHAOCLL7744ivIXuXl5R7HJU90gJ8wm80AgOjoaMn26Oho52dmsxlRUVGSzzt06IBu3bo5y7iSlZWFzMxM5/uamhr06NEDJ0+eVPYOxMscd06B9gQcqPUGArfurLdvsd6+52gR7Natm6Lf41GQXLx4MVauXNlsmdLSUgwYMKBNlfI2rVYLrVbbZLterw+4PwwA0Ol0rLePBWrdWW/fYr19LyRE2UEaHgXJBQsWYMaMGc2W6d27d6sqYjAYAACVlZWIiYlxbq+srERiYqKzzLlz5yT71dfXo7q62rk/ERGRt3gUJLt3747u3bsrUpFevXrBYDCgoKDAGRStVisKCwudPWSTk5NRU1OD4uJiDBs2DACwf/9+2O12GI1GRepFRETqpdhz6smTJ1FSUoKTJ0/CZrOhpKQEJSUlkjGNAwYMwO7duwEAGo0G8+fPx3PPPYd33nkH3333HaZNm4bY2FhMnDgRADBw4ECkpaUhIyMDRUVF+OKLLzB37lykp6fL9mx1RavVIicnx2UTrD9jvX0vUOvOevsW6+17Pqu7Uj2Cpk+f7rIn0ieffOIsA0Bs2bLF+d5ut4slS5aI6OhoodVqxZgxY8SJEyckx62qqhJTpkwRXbp0ETqdTsycOVNcvHhRqdMgIiIVU3ycJBERUaDi3K1EREQyGCSJiIhkMEgSERHJYJAkIiKSEZRBMpCX6fL0O3744QdoNBqXr507dzrLufp8x44d7VZvALjrrrua1OmRRx6RlDl58iTGjx+Pzp07IyoqCk8//TTq6+vbrd7V1dWYN28e+vfvj06dOuHWW2/F448/DovFIinn7eu9fv169OzZE+Hh4TAajSgqKmq2/M6dOzFgwACEh4cjISEB+/btk3zuzt+7t3hS940bN+KOO+5A165d0bVrV6SkpDQpP2PGjCbXNi0trV3rvXXr1iZ1Cg8Pl5Tx1TX3pN6u/g1qNBqMHz/eWcYX1/uzzz7Dfffdh9jYWGg0GuzZs6fFfQ4cOIDbbrsNWq0Wffv2xdatW5uU8fTfjUvt2rdWIdnZ2WLt2rUiMzPTrcnVhRBixYoVQq/Xiz179ohvvvlG3H///aJXr17il19+cZZJS0sTQ4YMEV999ZX45z//Kfr27SumTJni1bp7+h319fXi7NmzktfSpUtFly5dJENj8H/DbRqXa3xuvq63EEKMHj1aZGRkSOpksVgk5zZo0CCRkpIiDh8+LPbt2yciIyNFVlZWu9X7u+++E7///e/FO++8I8rKykRBQYHo16+fmDRpkqScN6/3jh07RFhYmNi8ebM4evSoyMjIEBEREaKystJl+S+++EKEhoaKVatWiWPHjolnnnlGdOzYUXz33XfOMu78vXuDp3V/4IEHxPr168Xhw4dFaWmpmDFjhtDr9eL06dPOMtOnTxdpaWmSa1tdXd2u9d6yZYvQ6XSSOpnNZkkZX1xzT+tdVVUlqfORI0dEaGioZGieL673vn37xP/8z/+IXbt2CQBi9+7dzZb/17/+JTp37iwyMzPFsWPHxAsvvCBCQ0NFXl6es4yn10JOUAZJB3dXILHb7cJgMIjVq1c7t9XU1AitViveeOMNIYQQx44dEwDEwYMHnWU++OADodFoxJkzZ7xSX299R2Jioviv//ovyTZ3/vBaq7X1Hj16tHjiiSdkP9+3b58ICQmR/Ni89NJLQqfTibq6unar9/XeeustERYWJq5du+bc5s3rnZSUJB577DHne5vNJmJjY0Vubq7L8n/84x/F+PHjJduMRqP47//+byGEe3/v3uJp3a9XX18vbrzxRrFt2zbntunTp4sJEyZ4tZ7X87TeLf3W+Oqat/V6//WvfxU33nijuHTpknObL653Y+7821m4cKH4zW9+I9k2efJkkZqa6nzf1mvhEJTNrZ6qqKiA2WxGSkqKc5ter4fRaITJZAIAmEwmREREONexBICUlBSEhISgsLDQK/XwxncUFxejpKQEs2bNavLZY489hsjISCQlJWHz5s0QXhoi25Z6v/7664iMjMSgQYOQlZWFy5cvS46bkJAgWRkmNTUVVqsVR48ebdd6N2axWKDT6dChg3SWR29c76tXr6K4uFjytxkSEoKUlBTn3+b1TCaTpDzQcN0c5d35e/eG1tT9epcvX8a1a9earPRw4MABREVFoX///pgzZw6qqqravd6XLl1Cjx49EBcXhwkTJkj+Rn1xzb1xvTdt2oT09HTccMMNku1KXu/WaOlv3BvXwsFvlspqT0ou0+VpPdr6HZs2bcLAgQMxcuRIyfZly5bh7rvvRufOnfHRRx/h0UcfxaVLl/D444+3W70feOAB9OjRA7Gxsfj222+xaNEinDhxArt27XIe19V/E8dn7VXvxi5cuIDly5dj9uzZku3eut4XLlyAzWZzeR2OHz/uch+569b4b9mxTa6MN7Sm7tdbtGgRYmNjJT92aWlp+P3vf49evXqhvLwcf/rTn3DvvffCZDIhNDS0Xerdv39/bN68GYMHD4bFYsGaNWswcuRIHD16FLfccotPrnlbr3dRURGOHDmCTZs2SbYrfb1bQ+5v3Gq14pdffsHPP//c5r89h4AJkoG6TBfgft3b6pdffsH27duxZMmSJp813jZ06FDU1tZi9erVzf5oK13vxoElISEBMTExGDNmDMrLy9GnT59WH9dX19tqtWL8+PGIj4/Hs88+K/msNdebpFasWIEdO3bgwIEDkk4w6enpzv+fkJCAwYMHo0+fPjhw4ADGjBnTHlVFcnIykpOTne9HjhyJgQMH4uWXX8by5cvbpU6e2rRpExISEpCUlCTZ7o/X25cCJkgG8jJd7ta9rUuBvf3227h8+TKmTZvWYlmj0Yjly5ejrq5OdoJgX9W7cZ0AoKysDH369IHBYGjSG62yshIAmj2uL+p98eJFpKWl4cYbb8Tu3bvRsWPHZsu7c71diYyMRGhoqPO8HSorK2XraDAYmi3vzt+7N7Sm7g5r1qzBihUr8PHHH2Pw4MHNlu3duzciIyNRVlbmlR/tttTboWPHjhg6dCjKysoA+Oaat6XetbW12LFjB5YtW9bi93j7ereG3N+4TqdDp06dEBoa2ub/hk4eZTADjKcdd9asWePcZrFYXHbc+frrr51lPvzwQ0U67rT2O0aPHt2kl6Wc5557TnTt2rXVdW3MW9fm888/FwDEN998I4T4teNO495oL7/8stDpdOLKlSvtVm+LxSJGjBghRo8eLWpra936rrZc76SkJDF37lzne5vNJm6++eZmO+787ne/k2xLTk5u0nGnub93b/G07kIIsXLlSqHT6YTJZHLrO06dOiU0Go3Yu3dvm+vr0Jp6N1ZfXy/69+8vnnzySSGE7655a+u9ZcsWodVqxYULF1r8DiWud2Nws+POoEGDJNumTJnSpONOW/4bOuvjUekA8eOPP4rDhw87h0IcPnxYHD58WDIkon///mLXrl3O9ytWrBARERFi79694ttvvxUTJkxwOQRk6NChorCwUHz++eeiX79+igwBae47Tp8+Lfr37y8KCwsl+33//fdCo9GIDz74oMkx33nnHbFx40bx3Xffie+//178/e9/F507dxbZ2dntVu+ysjKxbNky8fXXX4uKigqxd+9e0bt3b3HnnXc693EMARk7dqwoKSkReXl5onv37l4fAuJJvS0WizAajSIhIUGUlZVJusXX19cLIbx/vXfs2CG0Wq3YunWrOHbsmJg9e7aIiIhw9vp96KGHxOLFi53lv/jiC9GhQwexZs0aUVpaKnJyclwOAWnp790bPK37ihUrRFhYmHj77bcl19bxb/fixYviqaeeEiaTSVRUVIiPP/5Y3HbbbaJfv35euXFqbb2XLl0qPvzwQ1FeXi6Ki4tFenq6CA8PF0ePHpWcm9LX3NN6O4waNUpMnjy5yXZfXe+LFy86f6cBiLVr14rDhw+LH3/8UQghxOLFi8VDDz3kLO8YAvL000+L0tJSsX79epdDQJq7Fu4KyiAZyMt0tfQdFRUVTc5FCCGysrJEXFycsNlsTY75wQcfiMTERNGlSxdxww03iCFDhogNGza4LOurep88eVLceeedolu3bkKr1Yq+ffuKp59+WjJOUgghfvjhB3HvvfeKTp06icjISLFgwQLJUAtf1/uTTz5x+bcFQFRUVAghlLneL7zwgrj11ltFWFiYSEpKEl999ZXzs9GjR4vp06dLyr/11lvi3/7t30RYWJj4zW9+I95//33J5+78vXuLJ3Xv0aOHy2ubk5MjhBDi8uXLYuzYsaJ79+6iY8eOokePHiIjI8PjHz5v13v+/PnOstHR0WLcuHHi0KFDkuP56pp7+rdy/PhxAUB89NFHTY7lq+st9+/KUdfp06eL0aNHN9knMTFRhIWFid69e0t+zx2auxbu4lJZREREMjhOkoiISAaDJBERkQwGSSIiIhkMkkRERDIYJImIiGQwSBIREclgkCQiIpLBIElERCSDQZKIiEgGgyQREZEMBkkiIiIZ/x+9ZX67TPkfjAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "differential phase fit error=-69.106 mRad\n",
      "0.9022853875626692\n",
      "0.8257137521107901\n"
     ]
    }
   ],
   "source": [
    "plt.scatter(x, y)\n",
    "plt.ylim([-1,1])\n",
    "plt.xlim([-1,1])\n",
    "plt.gca().set_aspect('equal')\n",
    "plt.show()\n",
    "\n",
    "from alphautil.analysis import fit_ellipse_coeff\n",
    "from alphautil.ellipse import get_sci_params\n",
    "coeff = fit_ellipse_coeff(x, y)\n",
    "bx, by, Ax, Ay, phi_d_fit = get_sci_params(coeff)\n",
    "\n",
    "print('differential phase fit error=%0.3f mRad' % ((phi_d - phi_d_fit)*1e3))\n",
    "print(Ax)\n",
    "print(Ay)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}