{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d7ecef97-cbe5-48f2-8a89-c60706d6d1b8",
   "metadata": {},
   "source": [
    "# Numeric interferometer precompute comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f62f9e4f-e643-4c11-981a-d2bd707db8f8",
   "metadata": {},
   "source": [
    "**What you'll learn:**\n",
    "- How to set up a numeric interferometer using `NumericBraggInterferometer`\n",
    "- How to plot interferometer arm trajectories from the numeric tree\n",
    "- How to define beamsplitter and free-evolution functions using `propagate`\n",
    "- How to simulate an atom cloud and compute interferometer output ellipses\n",
    "- How to sweep detuning and compute output port populations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "efe528bd-1cf4-4dac-a1ac-c074ebf4df76",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from mwave.numeric import NumericBraggInterferometer\n",
    "from mwave.integrate import make_kvec, make_phi, propagate, omega_fnc_gaussian, multi_omega_fnc, phase_fnc_constant, score_backends\n",
    "from mwave.utils.cloud import cloud_init\n",
    "import numpy as np\n",
    "from numba import jit, float64\n",
    "\n",
    "# Set interferometer parameters\n",
    "omega_r_conversion = 2066*2*np.pi\n",
    "nbragg = 4\n",
    "T = 4.5e-3*omega_r_conversion\n",
    "Tp = 8e-3*omega_r_conversion\n",
    "    \n",
    "# Initialize kvec and tree\n",
    "ifr = NumericBraggInterferometer(-2*nbragg, 4*nbragg, distance=0) # Should pass in wavefunction initialization functions\n",
    "\n",
    "ifr.split(nbragg) # Should pass in the relevant function here as well! Function will inspect input function arguments to ensure correct number and naming\n",
    "ifr.propagate(T)\n",
    "ifr.split(nbragg)\n",
    "ifr.propagate(Tp)\n",
    "ifr.split([3*nbragg, -nbragg])\n",
    "ifr.propagate(T)\n",
    "ifr.split([3*nbragg, -nbragg])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed0e46ae-3836-4343-9aea-b28c3c23872b",
   "metadata": {},
   "source": [
    "Plot the trajectories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c8bee6ae-a421-460a-bdf3-2bfa5d5b5bbc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGdCAYAAADnrPLBAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAmVRJREFUeJzs3QV0FNfbBvBnXeJuRAkaSEKQENzdocW1DqVAnbpD6dciLVDH3d2luIWEhCAB4u7ZJOu7852ZDbT8S1sISXY3eX8cDrOzk83lLrAPd+69L49hGAaEEEIIIVaKb+4GEEIIIYQ8DQozhBBCCLFqFGYIIYQQYtUozBBCCCHEqlGYIYQQQohVozBDCCGEEKtGYYYQQgghVo3CDCGEEEKsmhD1gNFoRFZWFuzs7MDj8czdHEIIIYQ8BnZf37KyMnh7e4PP59fvMMMGGV9fX3M3gxBCCCFVkJ6ejgYNGtTvMMOOyNzvDHt7e3M3hxBCCCGPQaFQcIMR9z/H63WYuX9riQ0yFGYIIYQQ6/JfU0RoAjAhhBBCrBqFGUIIIYRYNQozhBBCCLFqFGYIIYQQYtUozBBCCCHEqlGYIYQQQohVozBDCCGEEKtGYYYQQgghVo3CDCGEEEKsGoUZQgghhFg1CjOEEEIIsWoUZgghhBBi1SjMEEIIIaTKvt22D2+t3wmj0QhzoTBDCCGEkCqZt3kPvnH2wRqvABw6cwnmQmGGEEIIIU+EHYX5bNNuLHbz5R7PyE9D/y7tYS5Cs31nQgghhFhlkPlo0x786unPPZ5dkI53nx1i1jZRmCGEEELIYweZdzfuxmqvAO7x20WZeP2ZwTA3CjOEEEIIeawg88aGXdjgHcg9/qAkG6+OHAhLQGGGEEIIIf/KYDDgtY27sc07EDyjEZ9U5OOl4f1hKSjMEEIIIeQf6fV6TN+4B7t9TEHmS1Uhpg3pC0tCYYYQQgghj6TT6fDipr044BMIvsGABboSTBjUG5aGwgwhhBBC/kar0WLalv046hMIgUGPhcZyPNu/JywRhRlCCCGEPESj1mDy1gM46RMAoV6P7/lKDO/TDZaKwgwhhBBCHlCqVJi07TDO+ARApNNiuViHQd26wJJRmCGEEEIIp7y8AuN3HcNFH39ItBr8JDeiX+eOsHQUZgghhBCC8rJyjN5zAtHefpBq1PjNXoCeHSKtomcozBBCCCH1XGlpGUbvP4lYL1/I1SqscBaja2RrWAsKM4QQQkg9VlRUjGcPncN1T1/YqpRY7SFHh9bhsCYUZgghhJB6qqCgEKOOXcQtTx/YKSuwztsB7Vq1gLWhMEMIIYTUQ3n5BRh54gruuHvDoaIMG/xdENGyOawRhRlCCCGknsnOzsWIM7FIdvOEc5kCm4I90LJ5E1grCjOEEEJIPZKekY1RF+KR6uoBV0UJNjdtgOZNgmHNKMwQQggh9URqWgZGXr6FDBd3uJcWY2sLfzQODoK1ozBDCCGE1AP3ktMwKuYOsp1d4VVShK3hwWgY6Ie6gMIMIYQQUscl3k3CqOupyHNyQYOiAmxr2xT+fg1QV1CYIYQQQuqwhNt3MPpWJgocnOBXmI9t7VvAt4EX6hIKM4QQQkgdFX/jNkbfzUWRvSMCC3KxvVM4vLw8UNdQmCGEEELqoKvxNzA2tRCldvZolJ+Drd1aw8PdDXURv6pfOG/ePLRt2xZ2dnZwd3fHsGHDcPv27YeuUavVmDFjBlxcXGBra4uRI0ciNzf3oWvS0tIwcOBAyOVy7nXeeust6PX6h645efIkIiIiIJFIEBwcjJUrV1a12YQQQkiddynmOkanFaPUxg7NcrOwo0fbOhtknirM/PHHH1xQuXDhAo4cOQKdToc+ffqgoqLiwTVz5szBnj17sGXLFu76rKwsjBgx4sHzBoOBCzJarRbnzp3DqlWruKDy0UcfPbgmOTmZu6Z79+6IjY3F7Nmz8fzzz+PQoUNP8/smhBBC6qRz0bEYl61AmdwGLXIysb1vB7i6uqAu4zEMw1THC+Xn53MjK2xo6dKlC0pLS+Hm5ob169dj1KhR3DW3bt1Cs2bNcP78ebRv3x4HDhzAoEGDuJDj4WG6h/fjjz/inXfe4V5PLBZzx/v27cP169cffK8xY8agpKQEBw8efKy2KRQKODg4cG2yt7evjt8uIYQQYnH+uBiNqUVaKKUyhOVkYPOArnBwsIO1etzP7yqPzPwv9huxnJ2duV+jo6O50ZpevXo9uKZp06bw8/PjwgyL/bVly5YPggyrb9++XOMTEhIeXPPX17h/zf3XeBSNRsO9xl9/EkIIIXXZsXOXMblYxwWZ1tnp2Daom1UHmSdRLWHGaDRyt386duyIFi1M1TZzcnK4kRVHR8eHrmWDC/vc/Wv+GmTuP3//uX+7hg0oKpXqH+fzsEnu/k9fX9/q+G0SQgghFung6fOYVmaEWiJFZFYaNg3pAVs7W9QX1RJm2Lkz7G2gjRs3whLMnTuXGym6/zM9Pd3cTSKEEEJqxN6TZ/GCSgCNWIJOmanYMKI3bG1t6lVvP/XS7FdffRV79+7FqVOn0KDBn7sJenp6chN72bktfx2dYVczsc/dv+bSpUsPvd791U5/veZ/V0Cxj9l7ZzKZ7JFtYlc9sT8JIYSQumzHsVOYaZRDLxKiW2YKVo3qD4m0/n3+VXlkhp03zAaZHTt24Pjx4wgMDHzo+datW0MkEuHYsWMPzrFLt9ml2FFRUdxj9tf4+Hjk5eU9uIZdGcUGlebNmz+45q+vcf+a+69BCCGE1EebDp/Aq4wceqEQvTOTsfqZAfUyyDzVaqbp06dzK5V27dqFJk2aPDjPzlG5P2LyyiuvYP/+/dxyazagzJw5kzvPLsO+vzQ7PDwc3t7eWLBgATc/ZuLEidzS66+++urB0mx2Hg57K2vatGlccHrttde4FU7sRODHQauZCCGE1CVrDxzD22InGPl89M9Mxs+jB3EDCHXNY39+M1XEfumjfq5YseLBNSqVipk+fTrj5OTEyOVyZvjw4Ux2dvZDr5OSksL079+fkclkjKurK/PGG28wOp3uoWtOnDjBhIeHM2KxmAkKCnroezyO0tJSrm3sr4QQQog1+33PYcbzaDTjcTyGeX7N9r99ZtYlj/v5XW37zFgyGpkhhBBSF/y0+xA+sXEDw+djZFYylowZAoFAgLqq1veZIYQQQkjN+WHHAXxs58EFmbFZyfh+7NA6HWSeBIUZQgghxMJ9t20fvnD04o4nZafg27FDwefTR/h91BOEEEKIBZu/ZQ8WOPtwx8/npGL+mCEUZKp7nxlCCCGEVD92d/0vtuzFMnc/7vH0vDR8MHowBZlHoDBDCCGEWGCQ+XjzHvzi4c89npWfjrmjh5i7WRaLwgwhhBBiYUFm7sbdWOUVwD1+qygTbzw72NzNsmgUZgghhBALCjJvbNiFDd6mXfXfL8nGzJEDzd0si0dhhhBCCLEA7K74r23cjW2VQeZTRS5eGt7f3M2yCrSaiRBCCDEzvV6PVzaYggzPaMSX5Xl4aejjlewxt9vnTyPtepxZ20BhhhBCCDEjnU6HFzfuwW6fQPCNRizQFuO5wX2s4j0pzcvFoR+XYMvn7yE1PtZs7aDbTIQQQoiZaDVaPLdlH46wQcZgwEKjAqP797SK94MxGnFw+ULo1Cr4NG0O35CWZmsLhRlCCCHEDDRqDaZsPYATPoEQ6vVYwlNiRJ/uVvNeXD2wBxk3rkMkkaLfK3PA55uvtAKFGUIIIaSWKVUqTNp2GGd8AiDSabFUpMWQ7l2s5n0ozEzHmQ2ruOOuE6fB0dNUasFcKMwQQgghtaiiQolxO4/ioo8/xFotfpIZ0b9LJ6t5Dwx6PQ788B30Oi0CwiIQ2sv8K64ozBBCCCG1pLysHKP3nEC0tx+kGg1+teOhV8f2VtX/l3ZuQW7SHUhsbNDn5dfA4/HM3SQKM4QQQkhtKC0tw+j9JxHr5QuZWo0VTkJ0a9/Gqjo/N+kuLmzfyB33nPYK7JxdYQloZIYQQgipYcXFJXj24FnEe/rCVqXEanc5OrQJt6p+12u1OLD0OxgNBjSO7IimHbvCUlCYIYQQQmpQQUEhnjl6ETc9fWCnrMBaL3tERphvGXNVnd28FoUZaZA7OKLn89Mt4vbSfRRmCCGEkBqSl1+AkScu446HNxwqyrDezxmtQ0Osrr8zbl7Hlb07uOM+L82E3N4BloTCDCGEEFIDsrNzMfJMLJLcvOBcpsCGhu4IC2lqdX2tVatwcPkigGEQ0q0XGraOhKWhMEMIIYRUs/SMbIy6EI9UVw+4KEqxuak3Qpo0ssp+/mPNbyjNzYG9mzu6T34RlohqMxFCCCHVKDUtA8MvJiDVxR3upcXY3ryB1QaZ5JgriDt6kDvu+/JsSOTyh543Go04sno5fl34FcyJRmYIIYSQanIvOQ2jYu4g29kVXiVF2BLeEMGB/lbZv6ryMhz6aQl3HNF/CPxahD4UYnZe3In441fB6IQAA8Tu2oDwoWPN0lYKM4QQQkg1SLyXjGfikpHr5AKf4gJsa90EAf6+Vtu3x3//ERXFRXDyboBO4yZz5wxGA7ae2YqrF65CppRxMULA8OEpFMOhTVuztZXCDCGEEPKUbty+i2dvZaDA0Rl+hfnY2j4Efg28rbZfb58/g1tn/wCPz0f/GXMAPh9rj69F/KV4yNQysD+EDB/NDA0QZm9Aw9engic0X6SgMEMIIYQ8hfibiRh9JwdF9o4IKMjF9o5h8Pb2tNo+rSgpxtHflnHHbYaMwMGkc7i9ZQWkGikXYhgYEa4PQAu9P1w9MuA4czJ4QvNVzGZRmCGEEEKqKOb6TYxNLkCJnT0a5edga7fW8HB3s9r+ZBgGh39aAnWZAnxnB+y/dwPS23JIIYWOr0OAjRw98qMggQi23klweHUiN3pjbhRmCCGEkCq4HHsd4zNKoLC1Q9O8LGztGQlXVxer7surR/Yh6eplMDweyhwaQKqXQSvQwrWZK8aUu8F4K4C7ztY3GQ6vWEaQYVGYIYQQQp7QuehYTMpVotzGFiG5mdjSpwOcnZ2sth9LKkqwavtP0O4/y+3ZonHzgcqWD68QZ0zqMwn8zfuhuG2azGwXmAL7FyZYTJBhUZghhBBCnsCpS1cxpVADpUyOsJwMbBrQBY4O9lbZh8XlxVhxcAUKbuTDITkZQqMROrkcjr1aYna/KZCLZCj7fQMUd/246+2D02A3bbxFBRkWhRlCCCHkMR07dxnPKQxQS2VonZ2OjYO6wc7ezur6L0+Rh1UHVqHkdglERhFsikogVJYDQgGmfrYAHr4BYIxGKH5Zh7Jk060l+ybpsJ86HpaIwgwhhBDyGA6fuYAXKnjQSKSIzErDuqE9YWtrY1V9l1WUhdUHVqPsXhkXYtgfOqYUdgWZ3PO9prz0IMiU/rgW5WmB3HmHkEzYTRwHS0VhhhBCCPkP+/44h5c1QujEYnTMSsWaEX0gl7GbxlmH1PxUrD2wFqpkFYSMkAsxKqkKLdqEgHfsKnINRgSERSC0V39TkFm6FuWZpiDjGJoD23FjYMkozBBCCCH/Yufx05hpkEInEqFrZgpWjuoHmVRqFX12N+cuNhzYAG2aFgJGACGEUMlUiIiKwKhOo3Bpx2acS7oLiY0N+rz8GrvFL0q+X4uK3CC2aAEcIwpg++wzsHQUZgghhJB/sPnwSczh28IgFKJXZgp+f2YAxBKxxffXzYyb2HRwE4wZRvDBhwACqGxUiOwYiWHth4HP5yM36S4ubNvIXd9z2iuwtXdE8eK1UOabgoxTu2LYjBgJa0BhhhBCCHmEtQeO4W2RI4wCAfplJuOX0YMgEoksuq/iUuOw9eBWIBtciGF/qO3U6NS5Ewa0GcCFGJZeq8WBpd/BaDCgcWRHNGnbAcUL10FZyAYZA5yjFJAPHQZrQWGGEEII+R8r9h7B+zIXGPl8DM5MxvIxgyE0Y+2h/3Ll7hXsPLITglwBF2BYagc1unftjr4Rff92/dnNa1GYkQa5gyN6Tn4exYs2QFXMBhk9nLsoIR8wBNbEct8ZQgghxAx+3n0IH9u4geHzMTIzGUvGDoFAYN7aQ//k/O3z2HtkL0QFIm4+DEvrrEXv7r3RvWX3R35Nxs3ruLJ3B3fce+rLUP68H+pSNsjo4NJDB1mfgbA2FGYIIYSQSj/sOIAvHL244zFZyfhu3NAHt2YsyR/X/8Ch44cgLhJzK5NYOlcdBvYaiA5NO/zj12nVKhxcvogtwoSQTt3geDwJ6jI2yGjh0puBrOffR3GsAYUZQgghBMB32/ZhgbMP1xcTs1Pw9VjLCzJHYo/g+MnjkJRIIIYYDBgYPAwY2mso2jZq+59f/8ea31CamwM7F1e0LHSGWs0uv9bAtR8f0m49Ya0ozBBCCKn3vt6yBwtdTbWHnstJxedjhlhMkDEajTgQfQCnT52GtEwKCSQwwgh4ASP7jkRYQNhjvU5ybDTijh7kjqPsWoJRNwIPargMkkDaqQusGYUZQggh9drnm3djqZup9tAreWn4cPRgiwgybIjZdXEXLpy9AFm5DFJIuRDDb8DH6L6jEeIb8tivpSovw+EfF3PHTZ0bw4nXDjyo4DrMFpL2/3xbylpQmCGEEFIvsWHh48178IuHP/f4tfx0vDfa/Kt4DEYDtp7diqvnr0KmlIH9YeAZIPITYUK/CQj2Cn7i1zz++48oLy6CvdgeIXaDwEMFXEc6QdI2EnUBhRlCCCH1Msi8t2k3Vnqaiii+VZiJN54dbNY26fQ6bD69GXEX4yBTm0KMnqeHLECGKf2nIMDd1NYndfv8Gdw6+wd44KGd61CIBBq4jXaDOLwN6goKM4QQQupdkHlz4y6s9zLVHnqvOAuvjTLfcmStXov1J9bj1pVbkGqkphDD18O2oS0m9J+ABs4NqvzaFSXFOPrLD9xxM8f2cJPawnWcF8Qtw1GXUJghhBBSbxgMBszauBtbvU1B5hNFDl4eMcAsbVFr1Vh7bC3uxtyFVMvOiJFCx9fBsbEjJvWbBA9Hj6d6fYZhcGjJN1BXlMNR7IEWTi3gNtEPouYtUddQmCGEEFIv6PV6zNi4B7t8AsEzGvGFsgDPDe1X6+2o0FRg9eHVSItLg0Qn4UKMVqCFazNXTO47GS52LtXyfeJ3bEZyQhxXmSnKvQs8pzSCqGlz1EUUZgghhNR5Op0OL2/ai30+geAbjZivKcKkwX1qtQ0KlQIrD65EdkI2JHp2gbUEGqEGXiFemNRnEhxtHKvtexVdi8aJzaYikqHOrdH4lfYQBjdBXUVhhhBCSJ2m1Wjx/JZ9OMwGGYMB3xkUGDOgV619/+LyYi7E5N/Mh9ggNoUYkQa+ob6Y2Gsi7GR21fr9dEl3sG/hj9AzOrhJPdD5zTEQNmyMuozCDCGEkDpLo9ZgytYDOOETCKFej0W8Cozq++iaRdUtT5GHVQdWoeR2CURGEbdjr1qkRlCrIEzqNQlSsbTav6f+XiLOfLsZeapsCHlCDJj1EsR1PMiwKMwQQgipk1RqNSZuPYQzPgEQ6XX4QajB0O5da/z7ZhVlYfWB1Si7V8aFGPaHWqJGk9ZNMLb7WEhF1R9iWLrEm7j382VcK7rCPe4y8hm4RrRDfVDlLQ5PnTqFwYMHw9vbGzweDzt37nzo+SlTpnDn//qzX7+HJ1oVFRVh/PjxsLe3h6OjI5577jmUl5c/dE1cXBw6d+4MqVQKX19fLFiwoKpNJoQQUk9UVCgxdhsbZPwh1mrxk1iPod071ej3TM1PxZerv8Sy75dBfUfNBRmVVIXgLsH4/K3PMbXP1JoLMjfikfv7HZzPPQ0Do4d/0xCEjxqH+qLKIzMVFRUICwvDtGnTMGLEiEdew4aXFStWPHgskUgeep4NMtnZ2Thy5Ag3OWvq1Kl48cUXsX79eu55hUKBPn36oFevXvjxxx8RHx/PfT82+LDXEUIIIf+rvKwcY/acwBVvf0g1Gvxix0PvjlE11lH3cu5h/cH10KZqIWAEEEIIlUyFiKgIjOo0CgK+oEbfJG18LArWZyKh+AaKtNmQyOToO+stbhChvqhymOnfvz/389+w4cXT0/ORz928eRMHDx7E5cuX0aaNaRfC77//HgMGDMD//d//cSM+69atg1arxe+//w6xWIyQkBDExsbiu+++ozBDCCHkb0pLyzB6/0nEevlCplZjhZMQ3drXzE63tzJvYeOBjTBkGCCo/KGyUSGyYySGtR9WK/WdtLHRKNiUiwK1EgklZ7lzPZ97BXbOrqhPanTOzMmTJ+Hu7g4nJyf06NEDX3zxBVxcTOvnz58/z42w3A8yLHYEhn3zL168iOHDh3PXdOnShQsy9/Xt2xdff/01iouLudclhBBCWMXFJXj24FnEe/rCRqXEKjcZOrVtVe2dE5cah62HtgJZ7FwNdhcXAdR2anTq3AkD2gyotSKVmiuXULC1CHqjDJcK1oEBg0aRHdC0UzfUNzUWZthbTOztp8DAQNy7dw/vvfceN5LDBhSBQICcnBwu6DzUGKEQzs7O3HMs9lf26//Kw8PjwXP/FGY0Gg338z72dhUhhJC6q7CwCKOOXMBNTx/YKiuwzssekRHVu9Nt9L1o7Di8A4JcARdiWGoHNbp37Y6+EX1RmzQXz6FgRxkY2OBG+S6Uakshd3BEr+dn1KvbSzUeZsaMGfPguGXLlggNDUXDhg250ZqePXuiJs2bNw+ffvppjX4PQgghliEvvwAjT1zGHQ9vOFSUY72fE1qHhlTb65+/fR57j+6FMJ+dDWP62NQ4a9Cnex90b1k7y7z/Sn32NAr3qMFAjhLmMm4U3ebO93lpJuT2DqiPam1pdlBQEFxdXXH37l0uzLBzafLy8v621TS7wun+PBv219zc3Ieuuf/4n+bisObOnYvXX3/9oZEZdiUUIYSQuiUnNw8jTsUgyc0LTuUKbAxyR1hI02p57VMJp3Dw2EGIi8Tc8mqWzlWHgb0GokPTDjAH9amTKNyvBwMp+NI7uFBwnS3ChJBuvdCwdSTqq1oLMxkZGSgsLISXlxf3OCoqCiUlJYiOjkbr1q25c8ePH+eqmUZGRj645v333+dWOolEpj9I7MqnJk2a/Ot8GXbi8f+unCKEEFK3ZGRmY9T5eKS4esBFUYrNTb0R0qTRU7/ukdgjOH7yOCQlEm6jO3YuisHdgKG9h6Jto7YwF9Wxoyg8wt5CkkBql4wYZwVKb+bBztUN3SfX7xW+VQ4z7H4w7CjLfcnJydxKI3bOC/uTvc0zcuRIbgSFnTPz9ttvIzg4mJvAy2rWrBk3r+aFF17gll2zgeXVV1/lbk+xK5lY48aN416H3X/mnXfewfXr17F48WIsXLiwOn7vhBBCrFRaeiZGXrqBdBd3uJcWY0uIH5o0aljl12P/I33w6kFuDzWpQsqVHDDCCHgBI/qMQHhgOMxJdfgQCo+z/6kXQeqQDEXPRoj79kvuuX6vzIFELkd9xmPYGuFVwM596d797/cKJ0+ejOXLl2PYsGGIiYnhRl/YcMLuF/P5558/mMDLYm8psQFmz5493OxvNvwsWbIEtra2D22aN2PGDG4JN3ubaubMmVyweRLsbSYHBweUlpZyG/QRQgixXveS0zAq5g6ynVzgWVKEreENERzoX+UQs+viLlw4ewGycpnpHIzg+/DxTL9nEOJbfXNvqkq5fz+KTrFtE0LmmAT5S4Oweu4clBcXoVX/wegx5SXUVY/7+V3lMGNNKMwQQkjdkHgvGc/EJSPX0RnexYXY3roxAvyffE6kwWjAtrPbEH0+GjKlKcQYeAaI/EQY128cgr2CYQmUu/eg6BxbiFIAuXMSnGaPw/4fF+HW2T/g5OWDiV8vhkhSM7sKW9PnN9VmIoQQYhVu3L6L0bcykO/oDN/CfGxrHwK/BqZpCY9Lp9dh8+nNiLsUB5lKBvaHnqeHLECGKf2nIMA9AJaiYsdOFF90NAUZtyQ4zRqPxMsXuCDD4/PR/9XX63SQeRIUZgghhFi8+JuJGH0nB0X2jggoyMX2jmHw9v7nVa3/S6vXYv2J9bh15RakGqkpxPD1sG1oiwn9J6CBcwNYkoot21EczW4yy4eNRxIcZ06AslyBo78t456PHPYMvIKbmLuZFoPCDCGEEIsWc/0mxiYXoMTOHsH5OdjWrTU83N0e62vVWjXWHluLuzF3IdVKwf7Q8XVwbOyISf0mwcPxz3mclqJ8wxaUXDMFNRvvJDi+OhHg8XD4pyVQlyngHtAQ7Uf+uZcboTBDCCHEgl2OvY7xGSVQ2NqhSV42tvZoCze3/647VKGpwJrDa5AalwqJTsKFGK1AC9dmrpjcdzJc7EyldSxN2dpNKL1uunVm65sMh1cmcreU4k8cRtLVyxAIheg/Yw4EQtN2JcSERmYIIYRYpPPR1zAxpxzlNrZonpuFrX2i4Oz87zX5FCoFVh1ahazrWZDo2QXWEmiEGniFeGFSn0lwtGHnoFgmxcoNUNwy3e6yC0iB/YsTuCBTmpeLk6t+4c53eHYCXP0sZ14PS3n1KoQuLhD7V21FWXWgMEMIIcTinLp0FVMKNVDKbRCak4HNA7rA0eGfV7MUlxdzISbvRh7EBrEpxIg08A31xcReE2EnY1cEWS7Fb+uguOPHHdsHp8Fu2nguyDBGIw4tXwStSgXvJs3RZvBwWBJDSQkyZ82Gobwcfr/+AnnlJri1jcIMIYQQi3L8/GVMKzVALZUhIjsdmwZ1g539o8NIniIPqw+sRtHtIoiN7H69YqhFagS1CsKkXpMgFVv2ah82rCh+WY+yZNOohn2TdNhPHf/g+ZiDe5B+Ix5CiQT9ps8Gny+AJcn5/Avo8/MhDgyENMR8e/JQmCGEEGIxDp+5gBcqeNBIpGiXnYb1Q3rC1tbmb9dlF2dj1f5VKLtXBpFRZAoxEjUaRzTGuB7jIBVZdoi5H2RKf1yL8rRA7rFD80zYTRr34PnCzHScXr+KO+428Tk4eT7ZMvSapjhwAIp9+wCBAN5fzwdfar4+pzBDCCHEIuz74xxe1gihE4vRISsVa0f0gVxm2tDuvrSCNKzdvxbKZCWEjJArAKmSqtCyXUuM7jIaIiuZGMsFmaVrUZ5pCjKOodmwHffnCiWjwYCDS7+DXqdFQFgEQnv1hyXR5eUh55NPuWPXl16ELDTUrO2hMEMIIcTsdh4/jZkGKXQiEbpmpmDlqH6Q/eV/+vdy7mH9wfXQpmohYAQQQgiVTIWIqAiM6jQKAgu7/fJvGL0BJT+sRUVOEPfYMSIPts8++9A1F3duRs69O5DY2KDPy6+Bx2MLTFoGhmGQ8+FHMJSWQtK8GVxfftncTaIwQwghxLy2HP4Ds/k2MAiF6JmZghXPDIBYIuaeu515GxsPboQ+XQ9B5Q+VjQqRHSMxrP0wrq6fNWH0ehQvXgdlPhtkjHBqWwSbkSMfuiY36S4ubNvIHfec+jLsnP97KXptKt22DeV//AGeSATv+fPBE5veK3OikRlCCCFms+7AcbwlcoBRIEDfzGT8OnoQRCIRrqdex+ZDm4Esdg9cPhdi1HZqdOrcCQPaDLC6EMNitDoUL14PZSEbZAxwjlJAPvTh1Ul6rRYHln7H3WZqFNkBTTt1gyXRZmQi96t53LHb7FmQNm4MS0BhhhBCiFms3HcE70ldYOTzMTgzGcvHDEZcWhx2HN4Bfg4bYUyBRe2gRveu3dE3oq/VvlOMVoOihRuhKq4MMp0rIB845G/Xnd28FoUZaZA7OKLX8zMs6/aS0YjsuXNhVCoha90azlOmwFJQmCGEEFLrft1zCB/K3cDw+RielYwJrT3w8U8fQ5gv5EZhWBpnDfp074PuLbtb9TvEqNUoXLgJ6lI2yOjg0l0LWd+Bf7su41YCruzdwR33fnEm5PYOsCTFa9ZAefkyeHI5vOd9BZ7AcuYpUZghhBBSq5buPIDPHby4477pt+FYEosjmyTcyiSWzlWHgb0GokPTDlb/zjBqFQq/3QJ1GRtktHDpbYSsZ7+/XadVq3Bw2UJ2di1CuvVCcJtIWBLNvXvI+24hd+zx9tsQ+5k2+LMUFGYIIYTUmkXb92G+kw93HJFyHQGpd8GDBAwYGNwNGNp7KNo2alsn3hFGWYGC77ZDU84uv9bAtR8f0m69H3ntH2t+Q2luDuxc3dB98ouwJIxOh6x33gWj0cCmc2c4jn545ZUloDBDCCGkVny9ZQ8Wuvpyx63SEtE29S4YGMF4ASP6jEB4YHideSeM5QoUfrcbGmUAeFDDZaAY0s5dH3ltcmw04o4e5I77vTIHErkclqTg55+hvn4dfHt7eH3xuUXN47mPwgwhhJAapdPp8NqWldjhZRpxaZ1yCxGpN8Dz4ePZfqMR4mu+bfBrglFRgoKF+6BV+YMHFVyH2kAS1fGR16rLy3H4x8Xccav+g+HXwrybz/0v1fUEFCz/kTv2/PBDiDw8/nbNrTspiL4ci/HjhsFcKMwQQgipEVqtBl/vWYS99o2RWhlk2qfGIgJpGPfSJAR7Bde5njeWFKFg0UFo1X7gQQnXkY6QtP3n+S/Hfl+O8uIiOHn5oPPYybAkRo0GWe++A+j1sOvbF/aDHp60HH/jLnauXA1pagyMPD6SO7VDoJ95Si5QmCGEEFKt1BolvtyzGAccQ5DhbFpOLWHUeBbrMLhRAiLbbIJE4lbnet1YVID8xUeg0/iCxyuH2zOuEEf88/yfxAtncOvsH+Dx+Og/43WIJJZVTyp/8RJo796DwNUVnp98/OD2UnTsDexbsxayjHjIwXDnyp39UaooN1tbKcwQQgipFuXKUnyxbxkOOoUhx8VUS0jKKNGl9DxeDQ6EVnEVGk0OrsaMRatWayGVeNaZnjfk56Lg+5PQaRuAz1PAdaw3xKGt/vH6ipJiHPl1GXfcbtgz8GrUBJZEeeUKilas4I69PvsMQicnnL98DYfXrYNt9g3cL/1Z5tYY3caMRbdO5p20TWGGEELIUykuKcQXR37GEecI5LmaQoycKUfXkgt4M7QHQoLf4M4pla0REzMBSmUyrl4di4hW6yCVWlYl6Kow5GYj/4fT0Ou8weeVwm2CH0Qhof9a2+jwz99DXaaAW0AQokb9WWDSEhjKK5D17lxumbjDyBG4KrHHiZmvwy4vEbaV15R7NkPvcRPQITIMloDCDCGEkCrJL8rGF8dW4qhzWxRWhhg7RoFuxRfwTptBCPZ/86Hr5XJ/RERswNWYCVCp0hDNBZq1kMlMK5yskSErA/nLzkOv9wKfXwK3SQ0hatr8X7/m+skjSIq+BIFQiAEzXofAwip95y1YAF1GBvTOLvgtvxjyJR/DjqskxYOqQQsMnDgRrcP//fdY2yjMEEIIeSKZeSn46o+NOObSDiWVIcaeKUHPogt4t9Oz8Pd8+x+/ViZrgNYR63E1ZiJUqhREXx3DjdDI5QFW9y7o01NR8OMV6A2eEPCL4DqtKUTB/367qDQvFydX/cIdd3h2Alz9LOv3rTj5B0o2b+aOrzhJIC9N40KM2r8VhkyagLAWllGL6X9RmCGEEPJYkjNuY/6FXTjuHIkyV9Mutk7GQvQquoT3uk2Al9vjFUVkby21rhyhUSrvPRihsbFpaDXvhD7lHvJ/iYPB4A6BoABuz7eEMDD4P2sbHVq+CFqVCt5NmqPN4IeLTJqT0WjE3i170ODLjyBj32tXB+Tb2kAb1AYjp05Cs8aWFbr+F4UZQggh/+pW8jV8c/UwTjq1R4VLH+6cqzEPfYqu4IPeL8DZsecT96BE4o6IiPWIiZmIiorEykCzBra2ljUR9lH09xKR/9sNGIyuEAjy4fZiKwj92V1+/13MwT1IvxEPoUSCftNng883f20jg96A7bsOI2H/NnROuAaZVosyiQjxXQZg1IvT0DCgAawBhRlCCCGPFJd4Gd8mnMIphyionE3b8HsYs9G3KBbv9X0RjvamYFNVErErd4spJnYyystv4GrMeLQKXw07O8uaj/FXusSbKFh5BwajC4TCXLi90g4Cn/+uU1SYmY7T61dxx10nPAcnT/NOfNbp9diydT8SD+2AnTIfjUrK4V1SDiOPB9dvF+GdXj1gTSjMEEIIecil66ew+M5VnHFoD03lqIu3IQMDiuPx3qCZkMtM82Sqg1jszI3IxMROQVlZPHfrqVX4StjbW9ZOuCzdzevIX50CI+MEoTAHbjM6QOBlqjP1b4wGA1dEUq/Twj+0FcJ6V1//PSmtVocNG3ch5dhu2KqLuIm9PAMPLXJKuOfdXn4J7lYWZFgUZgghhHBOxxzB0rSbOGvXHjpH0/wXP30qBpbexLtDZkMiHlQjPSUSOXKBJjZ2KkoVMdzkYDbQODj88z4ttU0bH4uC9ZkwMg4QibLgOrMrBO6Pt0/OpZ1bkHM3ERIbG/R9eZZZahup1BpsWL8dGSf3wUZTwi2x1vAlkIZ1xZDEOOg0dyFp3gxur7wCa0RhhhBC6rkjl/bgp5w0XLBtB719F+5coP4eBpfdw1uD50AkGlrjbRAK7RAevhLXrj2PktLL3K2nsLDf4ORo/gra2thoFGzKhZGxh0icCbdZPcB3cX+sr81Nuovz2zZwxz2nvgw7F1fUpvIKFdat2YK8Mwcg15Vxm92pBTLYtumJKVPHQnDsCLJXLgdPJIL3/PngicWwRhRmCCGkntp3dgt+KSrEZZs2MNiZ9noJ1iVimDIdswa+BpGodvc/EQptER7+O67FvYji4vPcSE1Y2C9wdoqCuWijLyF/axEYxg5iSQZcZ/cB38nlsb5Wr9XiwNLvuNtMjSI7oGmnx1vtVR0UZUqsWbkBxRcOQ6avAFuHWyWUwymyD56fMgYO9rbQZmQi+at53PVus2dB2tgyl10/DgozhBBSz2w9uRarypW4Im8NxrYRd66p9gZGavIxc8gss7ZNIJAjLPRXxMW/jKKi07h27TmEtvwJLi6da70tmovnULCjDAxsIJalwXXWAPAdnR77689uXovCjDTIHRzR6/kZtXJ7qaikDGtXrEfZ5aOQGlTcMmulyA7unfrhpYnPwtZG9mCZePbcuTAqlZC1bg3nKVNgzSjMEEJIPbHuyO9YqwNipeFgbPjcuRbqODxrLMOLA2fAUggEUi7AXL/+KgoKj3MjNaEtl8HVtXuttUFz7jQKdqvBQA6JPBUusweDb+/w2F+fcSsBV/bu4I57vzgT8if42qrIKyjBut/XQhVzAhKjBmzJygqxA3y6DsL0CSMgk0oeur54zRooL18GTy6H97yvwBOYf5n406AwQwghddxvB37CJkgRJ4148K9+mDoW4wRaTO7/IiyRQCBBy5ZLcT1hFvLzDyMu/hW0aLEY7m6mKtw1SX3qJAr368FAColNClzmDAff9n5Vov+mVau41UtsbaOQrr0Q3CayxtqalVOA9b+vhS7+FMRGLdjIUi51RkCPIZg5dijE4r/fKtQkJSHvu4Xcscfbb0Hs999Lyy0dhRlCCKmjlu39HluFzrghNX2Y8hgjWquuYoJMiDH9Lf+2Ap8vRouQJUi48Qby8vbh+vWZCGm+EB4eA2vse6qOH0PhYfZIAqldMlzmjOJGL57EqbW/ozQ3B3aubug+5YUaaWdqRjY2/r4WzI1zEDE6sNN2y2SuaNx3OJ55ZiBEwkd/vDN6PbLeeReMRgObTp3gOHo06gIKM4QQUofodDos3fcDtsu8kGhjmmfCZwxoq7yCKfYOGD7weVgTPl+EkObfcb/m5OzE9YTZYBg9PD2rf4WV6vAhFB5nRzJEkNonw+X1Z8CTmuaYPK7k2GhcO3KAO+73ymxI5Oz6oeqTlJKJTb+vBj/xAoSMgTtXZuOB5gNGYuSwvhAI//12UcHPP0MdHw++vT28vvzCLMvEawKFGUIIqSMh5ru9i7DbJhD3HExzSwSMHpHll/CCmxf693gJ1orPF6J5swXg8UTIzt7CjdQYjTp4e4+qtu+hPLAfRX+wwUUImWMSnF8fA5744Xkm/0VdXo7DPy7mjlv1Hwy/FmHV1r5bd1KwbcVqiO5dgRhG7pzCzhuthjyLIYN6gM83zYH6N6rrCShYtpw79vzwQ4g8PFBXUJghhBArDzHzd3+HvfaNkepoKjkgYrSIKruIV3waonuP6agLeDwBmjX9ihuhycxcj5u33gHD6ODjM/apX1u5Zy+KzrJzYgSQOSfBefY48B4x1+S/HPt9OcqLi+Dk5YPOYyejOsTfuIudq1ZDkhIL6f0Q4+iHdsNGo3/fzo8VYlhGjQZZ774D6PWw69sX9oNq7ladOVCYIYQQK6TWKPHVniXY79gcGc6mSbFiRo2OiguYGRiKDj1moq7h8fho0vgzboQmI2MVbt3+AEZGB98Gk6r8mhU7dqL4oiMXZORuSXCaNR68f5hv8m8SL5zBrbN/cG3sP+N1iCTseqKquxp7E3vXrIEsIx5yMNw5hUsQOo4cjT49Oz7x6+UvXgLt3XsQuLrC85OP68ztpfsozBBCiBUpV5bii33LcNApDDku/bhzUkaJzqUXMKdpe0T0mI26jP0QbtzoQ26EJi3tVyQmfgqjUQt/vyefC1SxZTuKo9kN8Piw8UiC48wJ4P3HnJNHvk5JMY78uow7bjfsGXg1qnrl7/OXr+HwunWwyb7B7dbLKnNrjG5jxqJbp6rthqy8cgVFK1Zwx16ffQah0+PvlWMtKMwQQogVKC4txJeHf8Jh59bIczUVKpQz5ehacgFvhvZASPDrqC/YQBPc8F3weWKkpC7D3bvzwBi1CAh4/Ftq5Ru2oOSaqbaSjVcSHGdULcgwDIPDP38PdZkCbgFBiBo1BlXxx9lonNi4DnZ5iVzdJK6Nns3Qe9wEdIis+twbY0UFsua+xy0TdxgxAnY9am+vntpEYYYQQixYflE2vji2Eked26LQdQB3zo5RoFvxBbzTZhCC/d9EfcQGmoYN3+CWbyclL8K9pG+5ScGBga/95y2UsrWbUHrdmzu29U2GwysTwXvMuSf/K+HkUSRFX4JAKORuLwmETzbX5uiJ8zi9ZQPsC5O4CtZG8KBq0AIDJ05E6/DmeFq5C76BLj0dQm8veLw3F3UVhRlCCLFAmXkp+OqPjTjm0g4llSMx9kwJehZdwLudnoW/59vmbqJFCAycCR5fjHv3FiA5ZQk3h6Zh0Bv/GGjKVm5A6a0G3LFdQArsX5xQ5SBTmpeLE6t+5o47PDsBbn4Bj/V1RqMRh46cwYXtm2Bfkgr7yhCj9m+FIZMmIKxF9dRIKj91CiWbNnHH3l/Ng+AJNv6zNhRmCCHEgiRn3Mb8Cztxwrk9FK6mOTFOxkL0KrqE97pNgJdb7RUrtBYB/i+BzxPhzt0vkZq6nLvlFBw892+BRvHbOijumHa7tWuYBvvnxlc5yLC1jQ4tXwStSgXvJs3RZvDwxwoxu/cdR8zuLbBXZHIhxgA+tEFtMHLqJDRr/Hhh6HEYSkqQ/f4H3LHTpImwaV9zuxBbAgozhBBiARKT4/H11YM46dQeFS6m1Umuxnz0LrqC93tNg6tTT3M30aL5+U0Djy9CYuInSEv/DUZGi8aNTKt22OCh+HU9ypL8uWvtG6fDftr4p/p+MQf3IP1GPIQSCfpNnw0+/5/n2xj0BuzYfRjx+7bBvjyHCzF6ngDGRu0wauokNAoyVSyvTjmffwF9fj7EgYFwf73uz6eiMEMIIWYUl3gZ314/hVOOUVA5m/aJ8TBmo29RLN7r+yIc7U3nyH/zbTCRG6Fhl2xnZKzh5tA0afQpFD+tR3laIHeNQ/NM2E0a91TdWZiZjtPrV3HHXSc8BydP0/yb/6XT67Fl634kHtoBO2U+F2J0PCF4zTpg9NSJCPDzqpG3VXHwIBT79rEFruD99XzwpU+3TNwaUJghhBAzuHT9FBbficYZhyhoKkddvA0ZGFAcj/cGzYRcZponQ56Mj88Ybtn2jZvvICtrIzQJt+CS9hrYG06OodmwHVe11Ub3GQ0GroikXqeFf2grhPX++/uk1eqwYeMupBzbDVt1ETexV8sXQdiiMyZOmwAfL/cae1v1+fnI+eRT7tjlxRcgCw1FfUBhhhBCatHpmCNYmnYT5+wioXU0LZP106diQOlNzB0yGxLxIHo/npKX10iAEeDGzTdR6BQLbcuf0Uw0BLbPPn1RxUs7tyDnbiJXc6nvy7MempejUmuwYf12ZJzcBxtNCbfEWsOXQBreDc9NmwB3t5rd34VhGGR/+BE3X0bSvBncXnkF9QWFGUIIqQVHL+3FjzmpuGDbDnr7Lty5QP09DC67h7cGz4FIVP2FE+srtjK0ZGMpvHnTkdXyR5R5XUCKuzNCjCO4UZuqyk2+h/PbNnDHPaa9DDsXV+64vEKFdWu2IO/MAch1Zdxmd2qBDLZtemLK1LFwcXJAbSjdvh3lJ0+CJxLBe/588MRsLe36gcIMIYTUoH1nt+DXogJcsmkLg51pSXCwLhHDlGmYNXAWRKKqf7iSRweZ4oXroCwMgh380bjRZNyxWYO8vP1gjDq0aLEYfP6TFZBk6bVaHPjhW+42U6PIDmjWqRsUZUqsWbURRecPQa6vgJwdnRHK4RTZB89PGQMH+9pbCq3NyETuV/O4Y7dZr0HauHqWd1sLCjOEEFIDtp1ch5XlFbgibw3GthF3rqn2BkZq8jFzyCzq8xrAaDUoWrQRqqIgbtGzc+cKyAe+D1lBB8Rfn478giOIi5+Oli2WQSB4skBzdvNaFGakQe7giDbPTMb3i35G2eWjkBpUXIhRiuzg3qkfXpr4LGxt2OrbtYcxGpE9dy63268sIgLOU6eivuEx7E22Ok6hUMDBwQGlpaWwt2fnkxNCSM1Yd+R3rNUBsdJwMDzTHiYt1HF4xqjASwNfpW6vIYxajaKFm6AqZYOMDi7dtJD1M+3TwyosOoO4uJdgNKrh7NQJoaE/QiB4vNCRcSsBmz55lysJwAuOgDr5BiQGNfdchdgB3l0HYvyEkZBJn3zEpzoUrVqF3HnzwZPLEbRzB8R+pr106tPnd9V2CwJw6tQpDB48GN7e3twEqJ07dz70PJuRPvroI3h5eUEmk6FXr164c+fOQ9cUFRVh/PjxXAMdHR3x3HPPoby8/KFr4uLi0LlzZ0ilUvj6+mLBggVVbTIhhNSY3w78hD4HVuENYQRiZBFckAlTx+Jr3SUc7T+JgkwNYtQqFH67+c8g08vwUJBhuTh3QnjYbxAI5CgqPoNr156HwaD8z9fWqlXYt+T/uCDDbnDH3L3KBZlyqRNcB0zBu7+txPPPjzNbkNEkJSHvu4Xcscfbb9WpIPMkqhxmKioqEBYWhqVLlz7yeTZ0LFmyBD/++CMuXrwIGxsb9O3bF2q1Kc2y2CCTkJCAI0eOYO/evVxAevHFFx9KZH369IG/vz+io6PxzTff4JNPPsHPP5u2jyaEEHNbtu8H9Di4Du9LIxEnDQOPMaKN8goWMddwqP8UTO7z579ppPoxygoU/N9WqMvYfWQ0cO3HQNbr0XvzODm1R3jYCggEtiguuYCY2KnQ68v+8bXTMnKxePZrKC/M5x4LYESZzBWeQ5/He7+twOTJoyAWi8w6PyjrnXfBaDSw6dQJjqOffrVWvb7NxI7M7NixA8OGDeMesy/Jjti88cYbePNNUxE0dojIw8MDK1euxJgxY3Dz5k00b94cly9fRps2bbhrDh48iAEDBiAjI4P7+uXLl+P9999HTk4OxJWzst99911uFOjWrVuP3T66zUQIqU46nQ5L9/+AHVIv3BY35c7xGQPaKq9gir0Dhnd5ur1MyOMxlitQ+N1uaJT+bGUjuAwUQ9q5639+XWlpLGKvTeGCjL19K4SH/Q6R6M9bGEkpmdj0+2oIbp/nAgxLKXVCk8FjMXJYXwiqUF27JuQvW4aCJd+Db2+PoD27IfLwQF3zuJ/fNTIBODk5mQsg7K2l+9jGREZG4vz581yYYX9lby3dDzIs9no+n8+N5AwfPpy7pkuXLg+CDIsd3fn6669RXFwMJ6dHr9nXaDTcz792BiGEVEeI+W7vIuy2CcA9e9MeMQJGj8jyS3jBzQv9e7xEnVxLjIoSFCzcB62KDTIquA61gSSq42N9rYNDOFqFr0FM7BQoFDGIiZ2IVuGrcC+lBNtWrIbo3hWIK0MMy75JGOZ88jn3+WQpVAkJKFi2nDv2/PCDOhlknkSNhBk2yLDYkZi/Yh/ff4791d394V0QhUIhnJ2dH7omMDDwb69x/7l/CjPz5s3Dp5+adkAkhJDqCDHzd3+HvfaNkepouoUhYrSIKruIl32C0KPHdOrkWmQsKULBooPQqv3AgxKuIx0gadv+iV7D3r4lIlqtRUzsJJSVXceBfb1wd4c7pGpTYNGJ5BDplHDy8sHE9z+0qCBj1GiQ/e67gF4Puz59YD+INlq0nHenGs2dO5cbkrr/Mz093dxNIoRYIY1WjY+2zkfUHwew1LkvUoWBEDNqdC89iU0u5dg8dCZ6tBlo7mbWK8aiAuQvPASt2hc8XgXcnnV+4iBz3517wMWj7aBTCmBjX4zGg1Og8vaGX6/hXJDh8fjoP+N1iCSWVdsof/ESaO7chcDVFZ6fmIpp1nc1MjLj6enJ/Zqbm8utZrqPfRweHv7gmry8vIe+Tq/Xcyuc7n89+yv7NX91//H9ax5FIpFwPwkhpCrKlaX4ct8yHHQKRbaLaVWMlFGiS+l5zG4ahYges6ljzcCQn4uC709Ap20APq8MrmO9IA5t9cSvc+FyPA6tWwOb7Bvc/+jvZvsjaGgWZM5qdB6UjOvrTVMT2g17Bl6NmsCSKK9cQdGKFdyx12efQejsbO4m1d0ww94aYsPGsWPHHoQXdt4KOxfmlcpaEVFRUSgpKeFWKbVu3Zo7d/z4cRiNRm5uzf1r2AnA7BDv/V0y2ZVPTZo0+cdbTIQQUlXFpYX48vBPOOwcgTxXUwFBOVOOriUX8GZoD4QEv0GdayaG3GwU/HAaOp0P+LxSuE3wgyjkyYoonjobjeMb18Mu7zZXN4lV7tkMUeMmILylA6KvTsCt3Vqoy8vh6ueLqFGWNZGb3RQva+573DJxhxEjYNfDNG+LPEWYYfeDuXv37kOTfmNjY7k5L35+fpg9eza++OILNGrUiAs3H374IbdC6f6Kp2bNmqFfv3544YUXuOXbbGB59dVXucnB7HWscePGcXNf2P1n3nnnHVy/fh2LFy/GwoWmNfWEEFIdCopz8fnR33HUuS0KXQdw5+wYBboVX8DbbQaikb9pVSYxD0NWBvKXnYde7wU+vwRukxpC1LT5Y3/9sZMXcGrzBtgX3uMqWLNLeCt8WmLAhIloG/Hn69goX4MidTV4fCO8O9+CVpcDmdAXliJ3wTfQpadD6O0Fj/fmmrs5dWNp9smTJ9G9+99T4eTJk7nl1+zLfvzxx9yeMOwITKdOnbBs2TI0/ku9CPaWEhtg9uzZw02uGjlyJLc3ja2t7UOb5s2YMYNbwu3q6oqZM2dyweZJ0NJsQsijZOal4Ks/NuKYSzuU8EzD9fZMCXoWXcC7nZ6Fvye7CRsxJ316Kgp+vAK9wR0CfhFcpzWFKPi/b/2wo/yHjpzBhe2bYF+SajrHLuD2D8fgSRMR3uLh2kWK/DysemsGtCoV/Dpq4dziHiQSL26SsFweAHMrP30a6S+Y9izyW7kSNu1NdzDqOsVjLs2mcgaEkHonJSMR8y7swAnn9lDwTBWNnYyF6FV0Ce91mwAvN8v533h9pk9NQv7PsTAY3CDgF8LthRYQBgb/Z4jZve84YnZvhb0igzvH7tyrC2qN4VMmoXmTh1fI3q9ttOWLD5CeEAfvxs0w7L3XEXttMpTKe5CIPdCq1RrY2DSEuRhKS5E0eAj0eXlwmjgRnu+/h/pCYc59ZgghxBIlpsRjQfRBnHBqjwqXvtw5V2M+ehddwfu9psHVqae5m0gq6ZPuIP/XBBiMbhAI8uH2YisI/f8eRO4z6A3Ysfsw4vdtg315DtiPPT1PAGOjdhg1dRIaBf1zQI05uIcLMkKJBP1mzIFM5oWIiPWIiZmAioo7uBozDq3CV8PW1jyTgXM+/4ILMuKAALi/PscsbbB0FGYIIXVeXOJlfHv9FE47tofS2bRPjIcxG32LYvFe3xfhaP/o7e+JeegSb6FgZSIMRhcIhXlwe6UtBD6Prjmk0+uxZet+JB7aATtlPhdidDwheM06YPTUiQjw+3NF7aMUZqbj9PpV3HHXCc/BydM0Z1MidkVEq3WIiZ2M8vKbuBozgQs0dnbNUJsUBw9CsXcvwOfD++v54MtqtyK3taAwQwipsy4nnMHixMs449Ae6spRF29DBvoXx+P9QTMhl5lWLBHLobt5Hfmrk2FknCAU5sBtRgcIvHz+dp1Wq8PGTbuRfHQXbNVF3MReLV8EYYvOmDhtAny8Ht6U9VGMBgMOLlsIvU4L/9BWCOv98J8HsdilcmO9KSgri8fVmPHcTsHshnu1QZ+fj5xPTBvAurz0ImRhYbXyfa0RhRlCSJ1zOuYIlqbdxDm7SGgdTQsV/PSpGFB6E3OHzIZETDumWiJtfCwK1mfCyDhCJMqG68wuELg/vKeYSq3BhvXbkXFyH2w0JdwSaw1fAml4Nzw3bQLc3R5/245LO7cg524iJHIb9H151iM3nxOJHBHRii19MPVB6QO2WKWDw5Pvb/Mk2EU02R9+BENJCSTNm8GtclsT8mgUZgghdcbRS3vxY04qLti2g96+C3cuUH8Pg8vu4a3BcyASDTV3E8k/0MZGo2BTLoyMPUTiTLi+1gMC1z9HV8orVFi3divyTu+HXFcGGwBqgQy2bXpiytSxcHEyTeR+XLnJ93B+2wbuuMe0l2Hn4vqP1wqFdmgVvhKx155DaekV7tYTW5zS0fHP2oLVrXT7dpSfPAmeSATv+fPB+0uNQvJ3FGYIIVZv39kt+LWoAJds2sJg14A7F6xLxDBlGmYNnPVg001imbTRl5C/tQgMYweRJANus/uA7+TCPacoU2LNqo0oOn8Icn0F5OzojFAOp8g+eH7KGDjY/7mVx+PS63Q4uPQ77jZTo3Yd0KxTt//8GqHQFq3CV+DatRdQXHKBu/UUHvYrnJyqVkrh32gzMpH71Tzu2G3Wa5D+ZUsT8mgUZgghVmvbyXVYWV6BK/LWYGwbceeaam9gpCYfM4fMMnfzyGPQXDyPgh0KMLCBWJoO19n9wXd0QlFJGdau3ICyy0cg1au4EKMU2cKtY3+8NOkZ2NqwZ6rm3Oa1KEhPhdzBEb1emPHYtY0EAjnCwn5FXPwrKCo6zY3UhIb+BBfnTtX2XrPLxLPfe4/b7VcWEQHnqVOr7bXrMgozhBCrs/7o71irYRAjawXGxlQvt4U6HqOMpXh54Kvmbh55TJpzZ1CwWwUGckjkqXCZPRgFWgbrvlkKVcwJSAxqsCUeK8QO8O46ENMnjIRM+nR19zJuJeDynu3cce8XXoXc/sluTwkEMoS2/Anx12egsPAE4uJeQMsWy+DqWj2lBYrXrIHy0iXwZDJ4z/sKPIGgWl63rqMwQwixGr8f/AkbGQnipBHg/qsOIEwdi3ECLSb3N+2OSqyD+vRJFO7Tg4EUEpsU6CZ0wXc/rIAu/hTERi3YyFIudYJ/jyGYMWYopJKnnzOiVau41UtsbaOQrr0QXMVq2wKBBKEtl+H69deQX3CEG6lp2eJ7uLk93RJ/TVIS8r4zlevxePstiP39n+r16hMKM4QQi7d87w/YInTEDYlpC3ceY0SE6iomSIUY23+KuZtHnpD6xDEUHGIr6UggkCVjtSEburmzIWJ0YCNLmcwVjfoMw7PPDoJIWH0fU6fW/o7S3BzYubqh+5QXnuq1+HwxWrT4Hgk3Xkde3n7EX38VISEL4eFuqu31pBi9HlnvvAtGo4FNp05wHGNZRS4tHYUZQohFYovPLt3/A3ZIvXDbxjQngc8Y0FZ5BVPsHTB84PPmbiKpAtWRwyg8xn70iFFkSMShm/vAZ7Rgp2iX2Xig+YCRGDmsLwTC6r29khIbjWtHDnDH/V6ZzS3Hflp8vgghzReCzxMjJ3cnrl+fBaa5Hp6eQ574tQp+/hnq+Hjw7e3h9eUXjz2Ph5hQmCGEWFyIWbh3MXbZ+OOevWkegoDRI7L8El5w80L/Hi+Zu4mkipQH9qPoD3YHWyHSK27hfN4e8GGEws4b4UNGYeigXlzR4eqmLi/HoR8Xc8et+g2GX4vq23yOzxeiefMF4PGFyM7eyo3UGBktvL1GPfZrqBISULBsOXfs+eEHEHl4VFv76gsKM4QQiwkx83d/h332jZDi2Is7J2K0iCq7iJd9gtCjx3RzN5E8hXu/r4f4tjd4PAFSy2/gYv5elDr6ot2w0ejft3ONhJj7jq/4EeXFRXDy8kHncZOr/fXZ31OzpvPA54mQmbUBN2++A8aog4/P2P/8WqNGg+x33wX0etj16QP7QbShY1VQmCGEmJVGq8aXuxdhv2NzZDibij+KGTU6Ki5gZmAoOvSYSe+QFbuy+SD0Z3PhIw3kbp0kl8XjuPEWol58F316dqzx75944QxunjkJHo+PftPnQCRh10dVP/b1mzT5HDy+CBkZq3Hr9gcwMnr4Npj4r1+Xv2QJNHfuQuDiAs9PPqbbS1VEYYYQYhblylJ8uW8ZDjqFItulH3dOyqjQueQ85jRrj4ges+mdsWIXVu0GL6YU3hJ/8GRB3Lm75XeR5OSJ9z56CSJxzS85rigpxpFfl3HH7YY9A+/GTWv0+7FhrXGjj7gRmrT035CY+AkYoxZ+fs898nrllSso+n0Fd+z1+WcQOjvXaPvqMgozhJBaVaIoxBeHfsJh5wjkuZoK+8mZcnQtuYA3Q3sgpMfr9I5YsbM/b4X4hgYNpH6A1FQnKUuVivIgOe4kBUFfYMC+H65hwPRQiKXCGq1tdPjn76EuU8AtIAhRo2pndRAbaIKD53KrnVJSl+PO3a9gNOoQEPDyQ9exm+JlzX2PWybuMHw47Hr0qJX21VUUZgghtaKgOBdfHP0dR5zboNDVtHzVlilDt+LzeLt1fzQOeJPeCSt2+vuNkN9j4C9tAHanOzZMZGpSIGjnjHbjJnDXNLtbgj0/XENmYgn2fn8Ng14Ng1hWMx9DCSePIin6EgRCIfrPeB0CYe2VtGADTVDQG+DxJUhOXoR7Sd/AyOgQGPDqg9tIuQu+gS49HUJvL3i8N7fW2lZX8Rj2T1wdp1Ao4ODggNLSUtjb25u7OYTUK5l5Kfjqj4045tIOJTzTMLo9U4KeRRfwboeR8Pc2lSEg1oetcXT2+42wTxfDTeLNnTMyRmSqkyHu7I3Wo/r87WtykxXYvSQWWpUeHoH2GDwzDBJ59QYNRX4eVr01A1qVCp3HTUG7oY+/sqi6paT8yIUZVoD/K1zIqThzBukvmDZ59Fu5EjbtTfsnkap/ftPIDCGkRqRm3cG8c9tw3DkKClfTnBgnpgg9Cy/i/W4T4OX238X9iOWGmNPfrodzrg0aSgLYve9gYAzI0CTDtlcQogb980aGbIAZNqcVdi2O4YLNrkWxGDIrHFIbUbXVNjq4fBEXZLwbN0ObwcNhTuztJXY/GvZ2E3vbyaAoA+/9k9xzThMnUpCpJjQyQwipVokp8VgQfRAnnNqjgmfHnXM15qN30RW832saXJ1oDw1rDjGn5q+Ba7EDnMXu3DmDUY90TRKchjRHSM+ox36tgoxy7F4cA1WZDi4+thg6Oxwyu6cvWXD1wG6cWPkzhBIJJi34Hk6ephEjc0vPWI3ExE/huEIA+WUBxAEBCNyxHXwZu+8O+Sc0MkMIqVXxd6LxbfwJnHKMgtLZVKPGw5iDvkUxeK/vi3C0f7q6NcR8dCoNTn29Bh4KFzQWN2I374XeqEOaNgmeI8PRqfOjV+v8G9cGthg2JwI7F8WgMLMcO76L4QKNjUPVC0kWZqbj9LqV3HHX8dMsJsiwfBtMguHUHagubwXDY6CdHgjeUxbNJH+ikRlCyFO5cuMsFt2+hDMO7aHmmf6X6W3IRP/iOLw/aCbkMlvqYSulKqvAua/XwUvlAXuRab6TzqhBmi4ZvuMiEdQ29Km/R3FOBXYtjEFFqRaOHnIMnd0Ktk5P/iFvNBiw4aO3kHM3Ef6hrTDyvc8sas8WfX4+kgYPgaGkBGX9DCgbYoCX1yg0a/oVt+keebqRGQozhJAqOR17FEtTEnDOvj20PNOHj68+FQNKb+K9IbMhEdfM5mSk5pUVleDSN5vgrfWGnciRO6c1qJGqT0bw1C7wDa3e/VpK85XYuTAG5UUa2LvJuDk1ds5P9ufnwraNOLt5LVdzafL/LYWdiyssBbvOJuOV6Sg/eRKSZs0gWzIeN+68y0YweHoMRbNmC7iyCOTv6DYTIaRGHL20Fz/lpOC8bST0Dl25cwH6JAwuu4u3B8+BSDSUet5KleQVIPrbbWhgaIAmwuZgqz+qDUqkG1PQ5KW+6N64Zm4VOrjJMfyNCG6ERpGvwo5vr3KBxt718eaT5Cbfw/ltG7jjHtNetqggwyrdvp0LMjyRCN7z50Pq2xh8iQwJCXOQk7uLW7Yd0vw7bqIwqRqKgoSQx7Lv7Bb8VlSAizZtYbBrwJ0L1iVimDINswbOgkhE/xBbq6KMbMQu3o0GjJ8pxAgBlaEC6UwqWs4chGB/U5mJmmTvIuMCDTtCU5pnCjTsLSf21tN/TUo+uPQ77jZTo3Yd0KyTZa2S02ZkIveredyx26zXIG3SmDv2cB/A7RQcf30m8vL2g2H0aBGymNtsjzw5us1ECPlX206uw8ryClyRtwZTeW+/qfYmhmtyMWsIlRywZrl3U5Gw/CD8+AGQCkyhQakvQzo/HRGzhsPJx7RiqTZVlGq4EZriHCXkDmJuhMbJ0+Yfrz+1bgUu794GuYMjd3tJbu8AS8EuE0+bMhXKS5cgi4iA/5rV4Akenh9TUHAC8denw2jUwsWlO1q2WAqBgCYG30dzZv6CNs0j5MmtP/o71moYxMhageGZKhq3UMdjlLEULw98lbrUimXcuIc7vx6FnyAQEoHpVk65rhSZoky0fWMU7N1czNo+pUKLXYtiUJRVAZmdiBuhYZdv/6/MWzew8ZN3uJIAQ9/8AMFt28OSFK1ezY3K8GQyBO3cAbG//yOvKyw8jbj4l2A0auDs3BmhLX+EQEBzzlgUZv6Cwgwhj2/FwZ+wgZEgThr+4FyYOhZj+VpM6WvatZRYp7SYG0hafRZ+okCI+ab//ZfpipElzUH7t8bAxtFydkhXlWuxe3EsCtLLuQ31hswOh5uvad8illatwpq3X0NJbjZCuvbkKmJbEk1SEpKHjwCj0cDz44/gNHbsv15fVHwe1669AKNRBSfH9ggL+wWCytGy+kxBq5mevDMIqc+W7/0BW4SOuCFpwT3mMUZEqK5iglSIsT3/eUdXYvnuXohB1sYr8BUHQVQ5J6NUW4Qc2zx0fHcCpHLL3LhNXaHDniWxyEstg0Qu5HYKdvc3/Rt+9NeluHbkAOxc3DD5/37gVjFZCkavR8rYcVDHx8OmY0f4/vrLYy0TLym5gthr02AwVMDRoS3Cwn6FUFi/tzZQUJh58s4gpL7R6XRYeuAH7JB44ra4GXeOzxjQVnkFk+1tMaLLeHM3kTyF2ycvIW9nPHzFgRBWrpQp1uYjz7EIXd6aCJHM8udmaFR67P0+FjlJCoilAgx+LRxqxV1sm/cx9/wzH34JvxZhsCT5y5ahYMn34NvZIWjPbog8PR/7a0tLYxB7bSr0+jLY27dCq/AVEAr/HJGqbxQUZp68MwipTyFm4d7F2GXjj3siU6FHAaNHZMUlPOfigYEdRpq7ieQpJBw6jeL9d+ArDYSAZ1q0WqjJRZFbGTq/NQFCK1t5plXrsW9pHLLulEAo0sGgWgtVWTFa9RuMHlNfgiVRJSQgZfQYQK+H9zcL4DB48BO/hkIRj5jYydDrS2Fn1xKtwldCVLnfT32joDDz5J1BSH0IMV/vXoi99sFIEQZx50SMFlFlF/GyTxB6tBlo7iaSpxC78xiUJ9LQQBoIfuXKswJNNkq81Og0Z6zVhZi/0mkM2LcsDslX18GovQU7F09MXfgDRBLLmShr1GiQMmoUNHfuwq5PH/gsXlTlXYjLym4iJnYSdLoi2No2R6vwVRCLTbsw1ycKqppNCLlPo1Xjy92LsN+xOTKc+3DnxIwGHRUX8GpgC3TsMZM6y4pd2XwQ+rO58JYGwFUWzJ3LU2eiPMCILrPHoS4QSQRo0q4c9y7cYmd0QWfsgex7FfBrbjlhJn/JEi7ICFxc4PnJx09VTsHOrhkiWq1DTOxElJffwNWY8WjVag0kYsvaENBS0D4zhNRh5cpSfLlvGQ46hSJb4MOdkzIqdC45j1lNI9GmeUdzN5E8hQurdoN3tRTeUv8HH5w56nSomgjR+ZVn61TfVpQUY+WbM6AuU8DFtysqyluDL+Sh/0stEdDS/B/wyuhopE6YyC0Tb7BsKex69KiW162ouIerMROg1eZBLm+ICDbQSOpP5XkF3WZ68s4gpK4oURTiy0M/4ZBzBPL4psmHcqYCXUrO463QHggJjjB3E8lTOPfLdogSVPCS+j04l6VOhSFUjqhpI+pc37K1jXZ+8zmSoi/BzT8QYz77PxxbmYik2HzwBTz0fb4Fglq5ma19xooKJA0bDl16OhyGD4f3vK+q9fWVyhQu0Gg02ZDJ/BHRai2kUsupCF6TKMxUoTMIsXYFxbn44ujvOOLcBoV80z/utkwZuhWfx9ut+6NxQEtzN5E8hdPfb4T8HgMPaYMHH/KZmhTw2zij3YRBdbZvr584gkM/LoZAKMT4rxZygcZgMOLoihu4eyUPPD4Pvac1R6M25hmxyP74E5Rs2gShtxeCdu2CwK76Vx+pVOlcoFGrMyCV+nKBRiYz/Tmoy2jODCH1SGZeGuadWo9jzu1Q7NqfO2fPlKBH0UXM7TAC/t5vm7uJpIrY2kNnftgEhzQRAiU+gBQwMkZkqpMh7uyN9qMm1em+VeTn4cSqn7njDs9O4IIMSyDgo/fU5tyvty/m4MhvCTAaGDSJfPxl0NWh/PRpLsiwvL/6qkaCDEsm80XriA3c3BmVKg3RV8dwc2rk8kfvKlzfUKFJQqxYatYdzDu3Dcedo6Bw6cedc2KK0LPwIt7vNgFebpZVdI88WYg5/e16OOXKESzxBySAgTEgQ5MM215BiBpU9zcyZGsbHVy+CFqVCt6Nm6HN4OEPPc8X8NFjcjPuVtPNc9k4uvIGjAYjmnWonVswhtJSZL//AXfsNGECbNrXbDkF9tZS64iN3AiNUpmEq1fHolWrtbCxMa1MrM8ozBBihRJT4rEg+gBOOkWhvDLEuBrz0bvoCt7vNQ2uTtUz+ZCYJ8Sc+noNXAsd0EgSZAoxRj3StUlwGtgUHXtPrTdvS8yhvUhPiINQIkG/GXPA5z9cpJHF5/PQfUJT8IV8JJzKxPHVt2DQM2jRxTThvSblfP4F9Hl5EAcEwP2N11Eb2Mm/EREbEBMzARUVd3A1Zixaha+Bra2pGnd9RWGGECuScPcqvok7jlOOUVBWLrH2MOagT9FVvN/3JTja9zZ3E0kV6VQaLsR4KFzQWNyICzF6ow5p2iR4jgxHp87P1au+LcrKwOl1K7njruOnwcnzn0db2DkzXcc2hkDAQ9yJDPyx/jY3QhPa3bfG2qc4eAiKvXvZNAXvr+eDL6u9khDs8mzTsu3JKC+/aVq2Hb6aW85dX1GYIcQKXLlxFotuX8QZhyionXpx57wNmehXfA0fDHoNcplpdIZYH1VZBc5+vQ7eKg80ETUBxIDOqEGaNhm+4yPRpW39G2UzGgw4sPQ76HVa+Ie2QlifAf/5NezS9E7PNoJAyEfMkTSc3nSHG6Fp1fvPFV/VRZ+fj5xPPuGOXV58AbKw2i+nIBa7cJOA2UBTVnadu/XUqtUq2NuZaqvVN7TPDCEW7Oy1Y/g++TrO2beHlmeqo+OrT8WA0pt4b8hsSMSWs2EYeTLlxaW4uGAjvLXesKvcql5rUCPNkIyGU7rAN7Rpve3SC9s34eymNVzxyMn/txR2Lo+/jwy7wuvi7iREH0jlHkcODUKb/gHV1jb29TOmz0D5iROQNGuGwE0bwRObineag06n4IpTKhQxXA2n8LCVcHD4s+K9taPVTIRYsaOX9uKnnBSct42E3qErdy5An4TBZXfx9uA5EImGmruJpIoU+YW4/H9b0MDgiybC5oAIUBuUSDemoPELvdGtaf2+VZibfA/nt67njtm6S08SZO6P0LQf2pAbobm0JxkXdyXBqDei7aDAp9qR977S7du5IMMTieA9f75ZgwxLJLLnilHGXnsepaVXuJGa8LDf4OjYBvUJ3WYixILsO7cNvxXm4aJNWxjsTHtIBOsSMbQiDbMHzYLIimvr1HdFGdmIXbwbDRg/NBGGcP/6qgwVSGdS0XLmIAT790V9x05+Prj0O+42U6N2HdCsc/cqv1bbgYFcoDm/4x4u70uBwcCg/dCgpwo02oxM5H41jzt2fW0mpE0sY9KtkBuR+R1xcS+iuOQCV3U7LPQXODnV7OoqS0JhhhALsP3UOqxQVOCKvDUY24bcuabamxiuycWsIbPN3TzyFHLvpiJh+UH48QPQWNCcO6fUlyGdn46I2cPRyIfmO913bvNaFKSnQu7giF4vzHjqkZSIvv7csu2zW+/i6sFUGPRGdBwZXKXXZZeJZ7/3Hrfbr6xVK7hMmwZLIhTaICzsV8TFvYyi4jOIvfYcwkJ/hrNz/ShZQmGGEDPacGwl1qgNiJG1AmPD5861UMdjlLEULw98ld4bK5Z5KwmJPx+BnyAQjUWmEFOhVyBDmIG274xCY7f/ntRan2TeuoHLe7Zzx71feBVye4dqed3wXn7cCM2pjYm4djQdRj2Dzs824lZAPYnitWuhvHQJPJkM3vPngSf4+zJxcxMIZAgN/Rnx16ejsPAkrsU9j5YtlsHVteojXNaCwgwhZrDy0M9YbxQjThoOyE3nwtSxGMvXYkr/F+k9sWJpMTeQtPos/EQBaCQ2hZgyXQmypNlo/+4YNHGkkir/S6tW4eCyhVyRxpCuPRHctnpvj7Ts1oAboTm5/jbiT2ZwpRC6jW3y2IFGk5SEvG+/44493n4LYn/L3XVXIJAgtOUyxF9/DQUFRxEX/wpatvgebm51ey4WhRlCatGP+37AFoEjEiTtuMc8xohWqhhMkPAwrr9lDVuTJ3P3QiyyNl2GrygIwRLTSqRSbRFybPLQ8aMJaCavvX1IrM2ptStQkpsNOxc3dJ9SM2E+pLMPN0JzfPVN3DidxU0K7j6xGbfp3r9h9HpkvTsXjEYDm44d4ThmDCwdny9ByxY/IOHG68jL24/4668iJGQRPNxNpU7qIgozhNQwnU6HZQeWYrvEA7flnbhzfMaANspoTLG3wYiB9WsztLrm9slLyNsZD19xIIIqQ0yJtgC5joXo8vFEhMhMS+rJo6XERuPakf3ccd9XZnHLsWtK0ygvboTm6MqbuHU+h9uHptcUthyC6RbvoxT+8gvUcXHg29nB68svqmVFVG3g80UIab4QfJ4IObm7kJAwC4xRB0/PIaiLKMwQUoMhZuHexdht44e7dqYaSQJGj8iKS3jOxQMDe9DtJGuWcOg0ivffga80EIFS06qWIk0eCtwU6PLWBLSglWf/SV1ezlXDZoX3HQT/ljW/P0rjdp5ceGELU965nMsVp+z9nKlg5f9SJSQgf+ky7tjzww8g8qzdIpZPi88Xonnzb8DjCZGdsw0JN94Aw+jg5TUSdQ2FGUJqIMR8vXsh9toHI8XRtFuviNGifdklvOITiB49plOfW7HYncegPJGGBtJAOMgacecKNNko8VKj05yxEFKIeWzHV/yI8uIiOHn5oMv42iucGdzanRuhOfTLddy7mseVPuj7fAsIRH8GGqNGg+x33wX0etj17g37wYNhjXg8AZo1mw8eX4SsrI24cfMdGBk9fLxHoy7557G1avDJJ59wQ3J//dm06Z+7WqrVasyYMQMuLi6wtbXFyJEjkZub+9BrpKWlYeDAgZDL5XB3d8dbb70FvV5fk80mpEo0WjU+2jofUX8cwA/OfZAiDIKY0aB76R/Y6FKGLUNfRY82A6l3rdSVzQdxftYqOJ8Xwk8WDD5PgDx1JpI80xC+8Fl0e3sSBZknkHjxLG6eOQkej49+0+dAJKnd3ayDwt3Q/+WW3Dya5GsFOPBTPPQ6w4Pn85csgebOXQhcXOD5qemzzFrxeHw0bfIFGjSYyM4Cwq1b7yE9Yw3qkhofmQkJCcHRo0f//IbCP7/lnDlzsG/fPmzZsgUODg549dVXMWLECJw9e5Z73mAwcEHG09MT586dQ3Z2NiZNmsRtHPbVV1/VdNMJeSxKVTm+2LsEB53CkFVZwVrKqNCp9DxmN4lEmx6zqCet2IU1e8C7UgJvqT94siDuXI46HerGAnSabvmTQS1RRUkxjvyylDtuN2wUvBubp3RDQEtXDJweiv3L45B6vRD7l8Wh/yuh0MXHouj3Fdw1Xp99CqGzM6wdj8dD40Yfg88TIy39NyQmfsLNofHzqxsLD2q0NhM7MrNz507Exsb+7bnS0lK4ublh/fr1GDVqFHfu1q1baNasGc6fP4/27dvjwIEDGDRoELKysuDh4cFd8+OPP+Kdd95Bfn4+xI+5jfTj1nYg5EmUKArx5aGfcNg5Arl80710OVOBLiXn8VZoD4QER1CHWrFzv2yHKEEFL+mfhQqz1KkwtJQj6rkRZm2bNWM/cnb93xe4d+Ui3PwDMf6r7yAQmndn68zbxdi7LA56jQENAqVofuRj6DMy4DB8OLznfVXn+v9e0rdITV3OPW7Y8G0E+L8ES/W4n981epuJdefOHXh7eyMoKAjjx4/nbhuxoqOjubkFvXqZ5hSw2FtQfn5+XJhhsb+2bNnyQZBh9e3bl/vNJSQk/OP31Gg03DV//UlIdSkozsXsLfPQ6fJVrHEdwAUZW6YMg4oO42CQPVaOeJOCjBU7vXQjomdvgN89Ny7IsP/4Z6pSkNWiFO0WTaAg85QS/jjGBRm+QIj+M143e5Bh+TRxwpCZYRBJBZAfWsEFGaGnFzzem4u6hsfjoWHQGwgMNI0Y37u3AMnJ38Pa1ehtpsjISKxcuRJNmjThbhF9+umn6Ny5M65fv46cnBxuZMXR0VQt9j42uLDPsdhf/xpk7j9//7l/Mm/ePO57EVKdsvPT8eUf63DMuR2KXU37NdgzJehRdBFzO4yAv/fb1OFWXBPozA+b4JAmQqDEB5ACRsaITHUyxJ29ETmKnWtAnpYiPw8nVv7MHXd4djw3MmMpvIId0b+TGuUHz3CP74ZNhp9ACsvb57d6Ak1Q4Gvcsu17Sf+HpORFMBq1CAp63WrnBtVomOnf/88NekJDQ7lw4+/vj82bN0Mmq7kNpObOnYvXX3/9wWN2ZMbX17fGvh+p21Kz7mDeue047hwJReWcGCemCD2LLuH9ruPh5WZadk2sM8Sc/m49nHLkCJb4AxLAwBi4EGPT0x9RQ2pvhU1dx9Y2Orh8EbQqJbwaN0XbIZZ1q85QWgr1EtMtpayA7kjS+KJsUSyGzAqH1Mb8o0c1ISDgFW6V092785CSugxGRofghu9YZaCp1aXZ7ChM48aNcffuXfTu3RtarRYlJSUPjc6wq5nYCb8s9tdLly499Br3Vzvdv+ZRJBIJ95OQp5GYEo8F0Qdw0ikK5S6misauxnz0KrqCD3pNg6tTD+pgKw4xp75eA9dCBzSSBJlCjFGPdG0SnAY2RYfeU83dxDon5tBepCfEQSiRoP/0OeDzLWvMI+fzL6DPy4M4IAChP3yKlB9vIT+tDDsXxmDorHDI7B5vjqa18fd7nhuhSbzzGdLSfuEmBTdq9IHVBZoanQD8v8rLy7k5MezE4MmTJ3MTgDds2MAtyWbdvn2bmzfzvxOA2VtU7LJs1s8//8wtz87Ly3vswEITgMmT2HtmB9Zl3cNN5+bQ82grJkKqBQMYDTwY9XLweBZW2oH9GDQYTcfs5nnsBznDjiZVfjzy8MSFKa0Ow4CBqQ94pt/wE7/EkfZN0cC5sthcNXncz+8a/Zf6zTffxODBg7lbS+yKpI8//hgCgQBjx47lGvfcc89xt4OcnZ25Rs6cORNRUVFckGH16dMHzZs3x8SJE7FgwQJunswHH3zA7U1DIy+kum08tg47inNx0TkKajfLuZdPSJ3BDsZY7ACHZY0UWWMf6O+HPzOo0TCTkZHBBZfCwkJuFKZTp064cOECd8xauHAh+Hw+NzLDrkBiVyotW2baOprFBp+9e/filVde4UKOjY0NN6Lz2Wef1WSzST3z+/5fsF9VgUtOUdC6hHDnGuhT0Tr7Klq6+YFf84v+SA0x6PUovpMFGSOBgGea92CEASqjCja+zpDa21Lf1xI9yqDknYezuMj0PjBAsd4HDZzGwcs+1CzvA/vRW7x6NTS3bkHo5Q3X6S+DJ3j4Y1Gp0CLmcCo0Sj13q6lVH786O4eGVVp6FZlZm7m/KQ72EfD2fobbRfhxeNhL6sdtJnOh20zkUZbt/AFHGAZXHCKh45n+uxigS0LrrBh88+xcbtdpYp0U+YW4/O1WNNA3gI3QNDStMaiQZkxG4xd6w6epafM7Uvuupm7Djbv/By9e3oNz2fBCaKN3Eeo7qFbbUrJtO7Lffx88kQgBW7dC2sRUY+t/KQpU2LUoBooCNexcpBg2pxXsXS3sVlk1ys3dh4Qbc8AwBri7D0RI82+5wpWW/PlNYYbUO99t/RZ/CGW4Yt8Whsr/rTfUJqJ1Vjzmj6YQY82KMrIRu3g3GjB+kAtNoy4qQwXSmVS0mD4A7kG0qtFSxGXsQ1ziPHgh+8G5bMYNzYPfQoR/zRdC1GVmImnIUBgrKuD2xutwfeGFf72+rEiNXQtjUJqvgq2TBEPntIKje939D09+/mHEX3+NK0zp5tYHLUIWg8+v/XuEFGaq0Bmkbpu3cT7Oyp1w1bYNjJXDpk00NxGRcwNf0kiMVctLSsf1ZfvhxwuAVGD6gFHqy5HBT0OrWcPh5GNaQEAsz42so7h6+wt4GNNxf45tjtEZwUGzEBk0ocaWiadNmQrlpUuQtWoF/7VrwBP8962UihINN0JTnKOE3EHMjdA4edqgriooOIG4+OlgGC1cXXqgZcsfwOfX7q0kCjNV6AxSN32y7gtcdPBErDwCTOUM/eaqeLTOu4NvptD8K2uWeSsJiT8fgZ8gEBKBadi/Qq9AhjADbd8YBXs3F3M3kTymxNzTuHTzM7gbkiCoDDW5Rgf4B0xHVNA0bn5ldSlavRq5X80DTyZD0M4dEPv7P/bXsnNo2EBTlFUBmb2YW7bt4lN3514VFp5GXPxLMBo1cHbujNCWP0IgqL2ioBRmqtAZpG55f/WnuOzihzh5qwfnWipj0aYwFfMmfWzWtpGnkxZzA0mrz8BPFAgx3/QPa5muBFnSbLR/awxsHOnvubVKLryCc9c/hKs+EcLKUJNntIVXg+fRpfGMpw41mqQkJA8fAUajgcdHH8J53Lgnfg1VuRa7F8eiIL0cUlsRhs4Oh2sDO9RVRUXncC3uRRiNKjg5RSEs9GcIKkdAaxqFmSp0Bqkb3l75Ma64B+OGrCX3mMcYEa68isiSbHwy4UNzN488hXsXryFz4yX4ioIgqrx/X6otQo5NHjrOnQCpvO5Oyqxv0ovicDrhAzhrEyCqDDX5BjncfCajW9PZEPCffDEuo9cjZdx4qOPiYNOhA3x/+7XKm8OpK3TYsyQWeallkMiF3E7B7v519/OluOQyrl17DgZDBRwd2iIs7FcIK+el1SQKM1XoDGK9lEol3t88D1c9muG2tDl3jscYEFEejQ4VhXh/bN0rGFef3D55CXk74+ErDoSwclVFibYAeQ6F6Pz2RIhktON3XZVdehsn49+DgzoWkspBmQKDFI6e49Cz+VsQCh5/UmrB8uXIX7wEfDs7BO3ZDdG/7CT/ODQqPRdocpMVEMuEGDwzDJ5BDqjLy7ZjYqfCYCiHg30rhIevgFBYsyNSFGaq0BnEOkPM3M3zEO3VEnfFpmWVAkaPNopL6KxX4c1Rb5i7ieQpJBw6jeIDd+ArCYSgcjfmIk0eCtwU6PLWBAhFdXe/D/Kw/LJkHIubCzvVFUj5ph1Figxi2LiPQu+Q9yES/vs8DvWNG0h+djSg18N7wddwGDKkWrpYq9Zj7w/XkH23FCKJAINmhsE7+OECynWJQhGHmNgp0OtLYW8XivDwlRCJai7AUZipQmcQ6woxb7EhxjscKaKG3Dkho0Pb0ovoyQNeHfaquZtInkLszmNQnkhDA2kg+JUrzwo0OSj2UqHznLEUYuqxImUmjlx7F/KK85BXhpoSgwhi1yHo0+JjSER/X11k1GiQMuoZaO7cgV3v3vBZsrhaaw/pNAbsW3YNmbdLIBTzMXBGGBo0cUJdVVZ2AzGxk6HTFcHWtjlaha+CWOxcI9+LwkwVOoNYPkVZCd7evhDR3hFIF5pWIIgZDdqVXMAAqRzTBvz7XhHEsl3ZfBC6s7nwkQaAX7nyLE+dhfIAPbrMHm/u5hELUqLMwZH4DyAp+wM2fFNNoVKDAHzn/ugb+jlkoj//rc/95hsU/fY7BC4u3O0loXP1f/DqtAYc+DEe6TeKIBDxMeCVlvBrXndX05WX38bVmInQ6QphY9MYEa3WQCx2rfbvQ2GmCp1BLFd+cR7e370UV7zbIkvYgDsnZVRoV3QBwxzdMK5XzexHQWrHhTV7wLtSAm+p/4P/Meeo06FuLECn6aPpbSD/qExdgMNxH0JQegx2AoPpnEEAg0NP9An9HIKEVKROmMgVUmyw9AfY9exZY72p1xlw8OfrSI0vhEDIR7+XWiCgZfV/wFuKiop7uBozAVptHuTyYC7QSCTVu6cThZkqdAaxPGk5yfjs8Gpc9m6HXIEXd07OVCCy4DxGewRgWNdR5m4ieQrnftkOUYIKXlK/B+ey1WnQt5Qh6rkR1Lfksal0ChyK+wjGov1wqAw1FQY+XE/I4LBPC6f+w+E9f16N96hBb8ThXxOQFJsPvoCHvi+0QFC4qR5hXaRUJnOBxmBQISJiPexsm1br61OYqUJnEMtxOzUBX5/ejktekSjgm5K+DVOG9vnnMCWwJXq3G2DuJpKncHrpRsjvMPCQmkbZ2BJxWepUMG0c0X7iYOpbUmUaXQUOX/8U2oLdcBToTCeVQJF9JHq2+QbOcp8a712DwYijv9/A3eg88Pk89H4uBMGt6+4u1CpVOnTchOAW1f7aFGaq0BnE/K4lRmPh5QO46BmFYr7pfrM9U4LI3PN4JaQjOoR2MXcTSRXpdTqc+WET7FNFcJd6c+eMjBGZ6hSIO3mi9TN9qW9JtSn54zju/fYyyvsZoK9cga028lAma4OeofPgZhdYo71tNBhxbPVNJF7MBXvntNfU5mjc7umWgtdHCio0+eSdQcznXNwpLE84i4vuUVDwTcsanYyFaJdzAbNa90VE0zb09lhxiDn93Xo45cjhKjH9Y25gDMhUJ8Ompz/ChvQwdxNJHWMoLUXS4CHQ5+XBYfwYxA8XoiR7PVwFau55jREolYajW8uv4OXQpMbaYTQyOLH2Fm6dywZ4QI+JzdCsg+l2OXk8FGaq0Bmk9h25fACrkuJwwS0K5TzTe+NizEdk9gW803kEmviH0NtixSHm1Ndr4FroAOfKSYEGox7p2mQ4DWyCkN4dzd1EUkdlvvU2FHv2QBwQgMAd28GXybg/eydvLUJ+5iq4CZTcdToGKBQ1R5cWX8LXObRG2sIYGfyx4TYSTmdxj7uNb4KQzjV/q6uuoDBThc4gtWfvmR1Yl3UPF13bQ8kzbYntbshBu+yL+Kj3JPh51uwQMKk5OpWGCzEeChc4Vi7V1Bt1SNcmwWNEGBp3aUvdT2qM4uAhZM6eDfD5CFi/DrLw8IeeNxqNOJW4FNkZv8KdX27688kA+cJG6BDyOYJcq//PJzsn7PTmO4g/kcE97jKmMVp2M80XI/+OwkwVOoPUvI3H1mFHcQ4uOneAmmeqo+NlyESbrMv4bMDz8HKlv+DWSlVWgXML1sNL6Q57kWkfD51RizRtEhqMbYeGkWHmbiKp4/T5+dztJUNJCVxeegnuc2b/47VsqDmf9DtSU5bBg1/KnTMwQJ4gEJHNP0Uj947VHmjObb+H2CNp3OOOo4IR3uvPVXzk0SjMVKEzSM1ZdeB37FEqcMkpClqeqY5OA30q2mRdxYIRc2BvV3e3/67rKkoUOL9gI3w0XrATmd5HrVGNNH0yGk7pAt/Q6l2qScg/hYWM6TNQfuIEJE2bInDzJvDEj1e36WLSWtxNXgxPXhH32MgAuXxfRDT5AM29e1VrGy/uSkL0wVTucfthQWjdL4De0H9BYaYKnUGq37KdP+AIw+CKQyR0PNM/LP66JLTJjsU3z7wLubx2ysiT6qfIL8Tlb7eigb4BbISmv1cagwppxmQ0fqE3fJoGUbeTWlOybTuy338fEIkQuHULpE2efGJvdOpW3Lz7Lbx4eQ/OZcMLoY3eRajvoGoLNFf2p+DSnmTucdtBgWg7MKBayyvUJRRmqtAZpPos3PodTgqluGLfFgaeqRhgQ20iWmfFY/7ouRRirFhRRjZiF+9GA8YPcqFpvpPKUIF0JhUtpg+Ae5CvuZtI6hldZiaShgyFsaICbm+8DtcXnq6sSVzGPsQlzoMXsh+cy2bc0Dz4LUT4j6yGFgPRB1NwYWcSd9y6nz8ihwZRoHkECjNV6Azy9OZtnI+zcidctW0DY2WBwCaam4jIuYkvn6WRGGuWl5SO68v2w48XAKnANKKm1Jcjg5eGVrOHw8mn7m4KRiwXYzQibeo0KC9ehKxVK/ivXQOewPRvz9O6kXUUV29/AQ9jOviVAyc5jDMaBc5Gu6CnrxUWezQNZ7fe5Y7De/miw8hgCjT/g8JMFTqDVN0n677ARQdPxMojwFQWCGyujkfr3Dv4Zspn1LVWLPNWEhJ/PgI/QSAkAtOk7Qq9AhnCDLR9YxTs3epuMT1i+YpWr0HuV1+BJ5MhaOcOiP1NBWirU2LuaVy6+SncDckQVIaaXKMD/AOmIypoGvh80795VRF/MgOnNiZyxy27N0DnZxtRoPkLCjNV6Azy5N5f/Skuu/ghTt7qwbmWyli0KUzFvEkfU5dasfS4W0hacQq+okCI+VLuXJmuBFnSbLR/awxsHOnvEjEvTVISkoePAKPRwOOjD+E8blyNfr+kgss4n/ARXPWJEFaGmjyjLbwaPI8ujWdUOdQknM7EyfW3AQYI6eyNrmObgHd/KKieU9AOwE/eGeTxvb3yY0S7N0SCzLTRFI8xIkwZg/YlWfhkwofUlVbs3sVryNx4Cb6iIIj4pknbpdoi5NjkocPb4yCzszF3EwkBo9cjZdx4qOPiYNOhA3x/+7XWRjTSi+Jw6vr7cNHdgKjyW+Yb5HDzmYxuTWdDwBc+8WvePJeN42tucoGmaQcvdJ/QlKvrVN8pKMw8eWeQ//b6yg8R7dEMt6XNucc8xoCI8mh0qCjE+2PnUhdascRTl5G7PQ6+4kAI+aZJ2yXaAuQ5FKLz2xMhkpmW1BNiCQqWL0f+4iXg29khaM9uiDxrv+5RdultnIx/D47qWIgrB2UKDFI4eo5Dz+ZvQSh4vKXh9yVeysHRFTfAMEDjdh7oObkZ+IKq38KqCyjMVKEzyKMplUrM3TwP0Z4tcFdiWu4oYPRoXXYZXbRKvPnMG9R1Vizh0GkUH7gDX0kgBDzT/yiLNHkocC1Fl7cnQigyBRtCLIX6xg0kPzsa0Ovh/fV8OAwdatb25Jcl41jcXNiprkDKZ7hzRQYxbNyfQe+Q9yASmm7TPg620vaR3xK4uk5spe1e05pDUI8DjYJGZp68M8jfQ8xbm+fhqlc4ksUNuXNCRoe2pRfRkwe8OuxV6jIrdm3XcZQfT4WvNBD8ypVnBZocFHup0HnOWAoxxCIZtVqkjBwFzZ07sOvdGz5LFlvMhNnCinQcjXsP8orzkFeGmhKDCGLXIejT4mNIRI93izYpNh+HfrkOo4FBULgb+jwfAoGwfgYaBYWZJ+8MUtlfZSV4e/tCRHtHIF1oWhkgZjRoV3IB/aVyPDfg6fZwIOYVveUQtGdy4CMNAL9y5VmeOgtl/np0fPVZCjHEouX93/+h8NffIHBx4W4vCZ1NpTMsSYkyB0fi34ek7BRs+EbuXKlBCL5zf/QN/Qwy0X9/DqXEF+DgT9dh0Bvh39IF/V5sAaGoepacWxMKM1XojPouvzgP7+9ZiitebZAlNG18JmHUiCw6j6EOzhjfe7K5m0iewoU1e8C7UgJvqf+D/8nmqDOgasRD5xljqG+JxVNGRyN1wkR2G100WPoD7Hr2hCUrUxfgcNyHEJQeg53AYDpnEMDg0BN9w76EreTfg1j6jSLsWx4Hg84I3+bOGPBySwjF9SvQKGhk5sk7o75Ky0nGZ4dX47J3O+QKvLhzcqYCkQXnMdojAMO6jjJ3E8lTOP/bdgjjVfCS/lnULludBl2IDB1eGEF9S6wCu7tv0rDh0KWnw2HYMHjPnwdrodIpcCjuIxiL9sOhMtRUGPnQ2HVF39CvYC/75w0nM24XY9/Sa9BrjfBp4oiB08MgktSfQKOgMPPknVHf3E5NwNent+OSVyQK+Ka/TDZMGdrnn8ekgOboG1k9tUiIeZxeuhGyOww8pQ0e1ITJUqeCaeOI9hMH09tCrEr2J5+gZOMmCL28ELR7FwR2drA2Gl0FDl//FNqC3XAU6LhzSiMPSpso9A6bD2e5zyO/LutuCfb+cA06tQFewQ4Y9GoYxNInX/5tjSjMVKEz6otridFYePkALnpGoZhv2r3VnilBZO55vNQ8Cp3Cupm7iaSK9Dodzv6wGXapQrhLvblzRsaITHUKxJ080fqZvtS3xOqUnz6D9Mp6S34rfodNVBSsmU6vxpGEL1GRtxXOAi13Tm3koUzWBj1D58HNLvBvX5OTVIo931+DVqWHR6A9Bs8Mg0Re91caKmhk5sk7o647F3cKyxPO4qJ7FBR8R+6co7GICzGvRfRC66aR5m4ieYoQc/q79XDKkcNVYtpvw8gYkK5Ohm1Pf4QN6UF9S6ySobQUSYOHQJ+XB6fx4+H54QeoK/QGLY7d+AalOevhIlBz5zRGoFQajm4tv4KXw8OVv/NSFdi9OBYapR7u/nYY/Fo4pDZ1O9BQmKlCZ9RVJ6KP4re70bjgFoVynun372LMR2T2BbzTeQSa+IeYu4nkKULMqQVr4FrgAGeJ6VahgdEjXZMMp/6NENK3M/UtsWqZb70NxZ49EAcEIHDHdvBlpvpgdYnBqMfJW4uQn7kKbgIld07HAIWi5ujS4kv4Opt2WmcVZJRh16JYqMt1cPW1xZBZ4ZDZPtnmfNaEwkwVOqOu2XtmB9Zl3cNF1/ZQ8my5c+6GHLTLvoiPek+Cn+ffhzKJddCpNDj19Rp4KFzgKHblzumNOqRrk+ExIhSNu7Q1dxMJeWqKQ4eROWsWwOcjYP06yMLD63SvGo1GnEpciuyMX+HOL+fO6RmgQNgYUSGfIcjV9Pe6MKucCzQqhRbO3jYYOrsV5PZ1M9BQmKlCZ9QVG4+vx46ibFx0joKaJ+fOeRky0SbrMj4b8Dy8XE0TQon1UZVV4NyC9fCscIeD2LSsU2fUIk2bhAZj26FhZJi5m0hItdDn53O3lwwlJXB56SW4z5ldb3qWDTXnk35HasoyePBLuXMGBsgTBCKy+ado5N4RxTkV2LkwBspSLZw85VygsXGseyVHKMxUoTOs3aoDv2OvUoGLTu2h5Zm2z26gT0WbrKtYMGIO7O1M82SI9akoUeD8go3w0XjBTmR6H7VGNdJ1yQic1Al+rUy1sgipC9iVdxnTZ6D8xAlImjZF4OZN4Inr5sjDf7mYtBZ3kxfDk1fEPTYyQC7fFxFNPoC3sAN2LYxBebEGDm4yDJ3TCnbOj186wRpQmKlCZ1irn3YvxyGDHpcdIqHjmf7C++uS0Do7Fv/3zLuQy02jM8T6KPILcfnbrWigbwAboenPrsagQpohGY1f7A2fpkHmbiIh1a5k23Zkv/8+IBIhcOsWSJs8PBG2PopO3Yqbd/8PXrz8B+ey4YVGXm/i9lp3lBWqYe8q5UZo7F3rzrwiCjNV6Axrs3DrdzgplOKKfVsYeKYZ7Q21iWidFY/5o+dSiLFixZl5iFm0Aw0YP8iFpvlOKkMF0pkUtJg+EO5Bph2aCalrdJmZSBoylNskz+311+H6IpVP+au49L2IuzMfXsh+cC7L6IbS6yPAu9UGtk4SDHu9FRzc6sZ/YinMVKEzrMX8jV/jjNwRV23bwFhZILCx5iZaZ9/El6NpJMaa5SWl4/qy/fDl+UMmMBWlU+rLkcFLQ6vZw+Hk8887hRJi7RijEWlTp0F58SI32dd/3VrwBPVnt9sncSPrKK7e/gIexnTwK+tsplbYozB+MNxK+3C3nJw8H6+wpSWjMFOFzrB0n6z7AhcdPBErjwBTWSCwuToerfPu4JvJn5m7eeQpZN5KQuLPR+AnCIREYBoirtArkCHMQNs3RsHezbS5ISF1WdHqNcj96ivwZDIE7djOLccm/y4x9zQu3fwU7oZkCCpDTVqFDUqThmD8yLfh1sB6P/NYFGaq0BmW6oM1n+KSsx/i5K0enGupjEWbwmTMm/SpWdtGnk563C0krTgFX1EgxHzTxL0yXQkyJdmIensMbByt788rIVWhSUpC8vARYDQaeHz0IZzHjaOOfAJJBZdxLuFDuOnvQFgZajKVMri5T0PfNrPB55v+A2xtKMxUoTMszTurPsYVt4ZIkJk2TOIxRoQpY9CuJAufTfjQ3M0jT+HexWvI2HAJfuIgiPimSdul2iLk2OShw9vjILOz/uFhQh4Xo9cjZdx4qOPiYNOhA3x/+/VBZXfyZNKL4vBH/Ptw1d2AqDK/5OllcG8wBd2azoaAb101nSjMVKEzLMXrKz/EVY+muCU17czLYwyIKI9G+4pCfDh2rrmbR55C4qnLyN0eB19xIIR806TtEm0Bcu0L0eWdiRDJ6t4+EYT8l4Lly5G/eAn4dnYI2rMbIk9TSQ5Sdam517H31JsIcroDcWWoKTBI4eg1Dj2bvQWhwDqWulOYqUJnmJNSqcS7m+fhqmcL3JWYliEKGD1al11GF60Sbz7zhrmbSJ5CwqHTKD5whwsx9/9nVKTJQ4FrKbq8PRFCUd2ur0LIP1HfuIHkZ0cDej28v54Ph6FDqbOqiVatx5YfDyDPeSkae9+BtDLUFBnEsHF/Br1D3oNIaNn70lCYqUJnmCvEvMWGGK9wJIsbcueEjA5tSy+iO2PAayNmmbuJ5Clc230c5cdS4SsNBL9y5VmBJgfFnkp0fn0chRhSrxm1WqSMHAXNnTuw690LPkuW0O2laqbTGLBv2TWkpt6FMWItgr1vQM5nuOdKDCKIXYegT4uPIRFZ5q1tCjNV6IzapCgrwTvbv0O0VwTSRKYZ+yJGg3YlF9BfIsPzA180dxPJU4jecgjaMznwkQaAX7nyLE+dhTJ/PTq++iyFGELYvxP/938o/PU3CJydEbR3D4TOphIdpHrptAYcWB6H9JvFgE0J7AfshANzHjZ8I/d8qUEIvnN/9A39DDKRZXxG3kdhpgqdURvyi/Pw/p6luOLVBllC08ZnEkaNyKLzGOrgjPG9J5u1feTpXFq7F8zlYnhL/R/8DzNHnQFVIx46zxhD3UtIJeXVq0gdP4GtXYAGP3wPu169qG9qkF5nwMGfryM1vhACIR9dnvPGTd13ECiOwa4y1JQZBDA49ETfsC9hK7GMYElhpgqdUZPScpLx2eHVuOzdDrkCL+6cjFEisvAcxrgHYFjXUWZpF6ke53/bDmG8Cl5SvwfnstVp0IXI0OGFEdTNhPwFu7tv0vAR0KWlwWHYMHjPn0f9UwsMeiMO/XIdydcKwBfw0PeFFvAKkeBQ3IcwFh2Ag8DAXVdh5ENj1xV9Q7+Cvcy8G3VSmKlCZ9SE26kJ+Pr0dlzyikQB3/SHwoYpQ/v885gU0Bx9IwfVantI9TqzbBOkiUZ4Shs8KJCXpU4F08YR7ScOpu4m5BGyP/kEJRs3QejlhaDduyCws6N+qiUGgxFHfruBe1fzwOfz0Pu5EAS3dodGV4HD1z+FtmA3HAU67lqlkQelTRR6h82Hs9zHLO8RhZkqdEZ1upYYjYWXD+CiZxSK+abdW+2ZEkTmnscLTdujS6vutdIOUv30Oh3O/rAZdqlCuEu9uXNGxohMdQpEHT3Q5tl+1O2E/IPy02eQ/oKp3pLfit9hExVFfVXLjAYjjq26icRLueDxeeg1pRkatzMth9fp1TiS8CUq8rbCWaDlzqmNPJTJ2qBn6Dy42QXWalvrXJhZunQpvvnmG+Tk5CAsLAzff/892rVrZ3Fh5mLCOSyNO4mL7lEo5Ttx5xyNRVyIeS2iF1o3jazR709qNsScXrgBTtkyuEpMf/GNjAHp6mTYdPdD+LCe1P2E/AtDaSmSBg+BPi8PTuPHw/PDD6i/zMRoZHBizU3cOp8D8ICek5qhaZRpCgRLb9Di2I1vUJqzHi4CNXdOYwRKpeHo1vIreDnUTiXzOhVmNm3ahEmTJuHHH39EZGQkFi1ahC1btuD27dtwd3e3iDBzIvoofrtzBRfdO6CMZ/oeLsZ8RGZfwBtRQxHS0LSLL7HOEHNqwRq4FjjAWWL682Zg9EjXJMOpfyOE9O1s7iYSYhUy33obij17IPb3R+DOHeDLTHXIiHkwRgYnN9zGjdNZXKDpNq4JQjo/fDvJYNTj5M1FyM9aBTeBkjunY4AicQg6h3wBX+ea/WyrU2GGDTBt27bFDz/8wD02Go3w9fXFzJkz8e6775olzLDVXXXKUhy6fAgb8tNxwTUKSp4t95y7IQftsi/ivY7D0cDdtGKJWB+dWovzP2yBm8IZDmLTzH69UYcsbRrcBjVHw/bh5m4iIVaj7PhJZM39EODz4ffzMshamnY4J+bFMAwu7ErGjbNZ3OOo4Q3RvMOfIzT3sZ+751JWIDd3DdwF5dw5PQPk84MRGTwX/k5hEMkdwKvmGlB1JsxotVrI5XJs3boVw4YNe3B+8uTJKCkpwa5du/72NRqNhvv5185gw091hpkTF/bhp5x7uGDfBmqenDvnacxCZPY1eNxVg8/GXEIIIaROYeDomAM/vzg4OOZxZ4xGHvLyghDV5SMENehkljBj8RWnCgoKYDAY4OHh8dB59vGtW7ce+TXz5s3Dp5/WbDVpRxs7xNiHckGmgSEd7bLi4JqkAR9sKqUgQwghpC7ioaTEi/tp75ALP794ODllw8MjCVqjymytsvgwUxVz587F66+//reRmeoUHtIJw39dALEOaHdNABHYe790/7dO4QEyVweIZJZdu4QQa2DQCyDVuMKFMX3slMOITD8pWg8OgoMTFVi1BAzDIOZIBmKOpHKPI/r4o1Xv//7sjM3cg/TCM+jZwHyLICw+zLi6ukIgECA3N/eh8+xjz3+orCqRSLifNYm9Lzj/xf+er0MIIeTPbfXP7UmE7Go+vA0COKUZULQ0EVf95GjzTDO4uFlmfaD6pP0IZ4jk9riwMwlXDhWA4dsgckjQv9bMatdkEtphEsypemfq1ACxWIzWrVvj2LFjD01EYh9H0f4EhBBiNURiAbqObIaITzshuYMHMoQMbMBDSJoKRd9G49D3l5GbXWbuZtZ7rfsFoOOoYK4fog+k4tz2e9yojSWz+AnA95dmsxN+f/rpJ25vGXZp9ubNm7k5M/87l8ZSyxkQQgj5+/b6F48kgTmXDX/TprPQgsEdTylajGgCHz8H6jIzijuRgdObErnj0O4N0OnZRrVe1bzOrGa6j12WfX/TvPDwcCxZsoRbsv04KMwQQojlYkfbLx9PgfZUJgJNm85CBwaJbhI0H94YvkGmDUhJ7Us4nYmT625zxyFdfNB1TGNu1+DaUufCzNOgMEMIIdYRamJOp6PsRBqCTZvOQg8Gt13EaDQkGEFNXM3dxHrp5rksHF9zi12VjWYdvNBtQlOurlNtoDBThc4ghBBiGWLPpaPoaCoaK03/3zawocZRhMDBDdEoxLyVnOuj2xdzcGzlDbDDH40jPbjyB3xBzU+7pTBThc4ghBBiWeIvZyLvYAqaVBgfnLtpJ0CDgUFoFv7oFa2kZty5kosjv9/gyiAEt3FHr6nNIajhQENhpgqdQQghxDLdjM1Bxr4kNCszPDh324YP934BaNn24XpCpOYkxebj0C/XYTQwCAp3Q5/nQyAQ1lygoTBThc4ghBBi2e4k5CF5zz00KdFBULnb+h05D069/BHegWrh1YaU+AIc/Ok6txotoKUL+r7YAkKRoEa+F4WZKnQGIYQQ65B0uwB3dt9Fk0IthJWh5p4UsOnmh4guvuBXc8FD8rC0G4XYvzweBp0Rfs2d0f/llhCKqz/QUJipQmcQQgixLulJxbixIxGN8zUQVYaaZDEg7uKDtj0CKNTUoIxbRdi3LA56rRE+TZwwcHooRJLqDTQUZqrQGYQQQqxTVnop4rfdRnCOGpLKUJMqAhDlhfZ9gmp0Xkd9lnWnBHt/uAaDwYihs1vBO9ixWl+fwkwVOoMQQoh1Y8shxG67jaAMJWSVoYYtm6Br54moAcEQUqipdjlJpVCV6xAYWv37AFGYqUJnEEIIqRsK8ytwZestBKRWcPWfWFkCBqoIN3QY3JirE0UsH4WZKnQGIYSQuqWkUIWL227BL6kMdpWhJpfPQBHmgg5Dm0AiFZq7ieRfUJipQmcQQgipmxQlKlzYfhs+iQo4VIaa/2/vzmOiuvo+gH+HfZd9k92yuOBWlUJj0z7yuPX1VYuJtUtqY2rc2rq0aU0XTf/RamramkbTf2qbNNYaij76WN/Hur4KrogWK6OIiMqioggiy8CcN+cwzMvUDa4MzDDfT3Id594LmZwf997vnHvuvbd0AtVDApAxPQWeXnKADdkahhkNjUFERH1b/b0m5GUXIbToLgJFW6i5rRO4kdIP6Vkp8PZx7+2PSB0wzGhoDCIicgwN9w3IzSlCUOEdBJtCTQ0EypP88NwryfDz9+ztj0hgmLHAMENERA/T1NiC3O16+J2pRpixLdTUQaAswQdpWQPhH8RQ05vYM6OhMYiIyDEZmluRu+MCPPNvIrK1LdTUQ6A01hujZqQgKMS7tz+iQ6rt5PFbJ4R8oHffxjBDRESd0dJiRN6uYrger0RUS1uoaYBASX9PDJ+RgrAIXzZkD2KY0dAYREREknyI4tH/lAB5FYg1tLVJEwSKwz2QmpWMyOh+bKgewDCjoTGIiIg6MhqNOLGvFM2HriO+uW2eAQIXQtwxaHoSohMC2GBWxDCjoTGIiIgeFWryD11F/YEyDGhsm9cCAX2QGxL/+xkkJHf/rfwJDDMdMcwQEVF3OX2kDDV7y5B4v23IaasMNf6uiJ8yAImDQ9nQ3Yg9Mxoag4iIqLP+PHEdN3aXIrneaJ533tcZUS8nYODwcDZkN2CY0dAYREREXXW+oBLX/l2CgXWt5nl6byeETYrHkFGRbNCnwDCjoTGIiIi0unjuBi7vuITkGgOcTc9/uuilQ0BmLIZnRLNhNWCY0dAYRERET6tEfwsX/1WM5OpmuJhCTbEH4PtSDEaMjYaTkxMbuZMYZjQ0BhERUXe5WnIHf+VcQNLNJriaQs1lN8Dthf4Y/Y84hppOYJjR0BhERETd7XrZXRT+pkdiZSPcTKHmiiugy4hA2j8T4OzCnppHYZjR0BhERETWUlVei4JsPRKuN8DTFGquuQgYxoQjffIzcGGoeQDDjIbGICIisrZbVfdwKrsIcWX34W0KNeXOAg0jQ5AxJQmubs4sggnDTAcMM0REZGtqqhtwLPs8YkruwdcUaqqcBGqHBSFjajLcPVzg6Gr51OyuNwYREVFPq61pwNHf9Oh/oRb9TKHmlk6gekgAMqanwNPL1WGLUssw0/XGICIi6i3195qQl12E0KK7CBRtoea2TqAqxR8ZWcnw9nF3uOLUMsx0vTGIiIh6W8N9A3JzihBUeAfBplBTA4HyJD+kzxgIXz/HCTW1DDNdbwwiIiJb0dTYgtztevidqUaYsS3U1EGgLMEXaVkp8A/yRF9XyzDT9cYgIiKyNYbmVuTuuACv/JuIaG0LNfUQKI31xqgZKQgK8UZfxTCjoTGIiIhsVUuLEXm7LsL1WBWiTKGmAQIlUV4YnpWMsAhf9DUMMxoag4iIyNa1thhx9D8l0OVVIMbQNq8JAsXhHkjNSkZkdD/0FQwzGhqDiIjIXhiNRpzYexnN/1uO+Oa2eQYIXAhxx6DpSYhOCIC9Y5jR0BhERET2GGryD11F/YEyDGhsm9cCAX2QG5KmJSI+Mai3P6JmDDMaGoOIiMienT5Shpq9ZUi8L9T7Vhlq/F0RP2UAEgeHwt4wzGhoDCIior7g7PHruPk/pUiuN5rnnfdzRtTkBAwcHg57wTCjoTGIiIj6kvMFlbj27xIMrGs1z9N7OyFsUjyGjIqErWOY0dAYREREfdGFwhso3XkJKTUGOJme/3TBS4fAzFgMz4iGrWKY0dAYREREfVmJ/hYubi9G8u1muJhCTbEH4PtSDEaMjYaTkxNsCcOMhsYgIiJyBGWX7uD8Nj2SbjbD1RRqLrsBbi/0x+h/xNlMqGGY0dAYREREjuR62V0U/qZHYmUj3Eyh5ooroMuIQNo/E+Ds0ruhhmFGQ2MQERE5oqryWhRk65FwvQGeplBzzUXAMCYc6ZOfgUsvhRqGGQ2NQURE5MhuVd3DqewixJXdh7cp1JQ7CzQ8G4KM/0qCq5tzj34ehhkNjUFERERATXUDjmWfR0zJPfiaQk2Vk0DtsCBkTE2Gu4dLjzQTw4yGxiAiIqL/V1vTgKPZevS/WIt+plBzUydwOzUAGdNS4OnlCmtimNHQGERERPSg+ntNyMsuQmjRXQSKtlBzWydQleKP57NS4OXjBmtgmNHQGERERPRoDfebkZujR1DhHQSbQk0NBMqT/JA+YyB8/dzRG8dvqw1PjouLg06ns5hWr15tsc7Zs2cxduxYeHh4IDo6GmvWrHng92zduhUpKSlqndTUVOzatctaH5mIiIgew9PLDeNeT8XAFRm4OCJQjaPxhw6DLtThavFt9BarXmv1xRdfoKKiwjy9++67Fmlr/PjxiI2NxalTp7B27VqsXLkS33//vXmd3NxczJo1C3PmzMHp06cxbdo0NRUWFlrzYxMREdFjyAHAL80cjNSVGbg0OhjnwtwxaGQEeotOCNH2nHAr9MwsXrxYTQ+zYcMGfPLJJ6isrISbW9u5to8//hjbtm1DUVGRej9z5kzU19dj586d5p977rnnMHz4cGzcuLHTn4WnmYiIiOxPr59mkuRppaCgIIwYMUL1vLS0tJiX5eXl4YUXXjAHGWnChAnQ6/W4c+eOeZ3MzEyL3ynXkfOJiIiIJKtdKP7ee+9h5MiRCAwMVKeLli9frk41rVu3Ti2XPTLx8fEWPxMWFmZeFhAQoF7b53VcR85/nKamJjV1THZERETUN3WpZ0aeBvr7oN6/T+2niJYuXYoXX3wRQ4cOxbx58/DVV19h/fr1FiHDWlatWqW6pdonObiYiIiI+qYu9cwsW7YMs2fPfuw6CQkJD52flpamTjOVlpYiOTkZ4eHhqKqqslin/b1c1v76sHXalz+K7AWSYapjzwwDDRERUd/UpTATEhKiJi0KCgrUI8VDQ0PV+/T0dDUA2GAwwNW17Q6Ce/bsUUFHnmJqX2fv3r0Wg4jlOnL+47i7u6uJiIiI+j6rDACWA3S//vprnDlzBiUlJfj555+xZMkSvPHGG+ag8tprr6nBv/Ky63PnzmHLli345ptvLHpU3n//fezevVudopKnr+Sl2ydPnsSiRYus8bGJiIjIDlnl0uz8/HwsWLBABRA5RkYO9H3zzTdVUOnYYyJvmrdw4UKcOHECwcHB6j40H3300QM3zfv000/V6anExER1Y73Jkyd36fPw0mwiIiL7w8cZaGgMIiIish02cZ8ZIiIiImtjmCEiIiK7xjBDREREdo1hhoiIiOwawwwRERHZNas9m8mWtF99zmc0ERER2Y/24/aT7iLjEGGmrq5OvfKRBkRERPZ5HJeXaPfoTfNsjdFoRHl5OXx9fdXDMLtL+zOfrl69yvvX2BDWxTaxLraLtbFNrAtUj4wMMpGRkeqRSA7dMyMbICoqymq/X97Ihzfjsz2si21iXWwXa2ObHL0u/R7TI9OOA4CJiIjIrjHMEBERkV1jmHkK8qGZK1assHh4JvU+1sU2sS62i7WxTaxL5znEAGAiIiLqu9gzQ0RERHaNYYaIiIjsGsMMERER2TWGGSIiIrJrDDNP4bvvvkNcXBw8PDyQlpaG48ePd19l6IlWrlyp7ujccUpJSTEvb2xsxMKFCxEUFAQfHx9kZWWhqqqKLdvNDh06hClTpqg7dMoabNu2zWK5vMbg888/R0REBDw9PZGZmYmLFy9arHP79m28/vrr6sZg/v7+mDNnDu7du8daWbEus2fPfmD7mThxIutiZatWrcLo0aPVHelDQ0Mxbdo06PV6i3U6s+8qKyvDyy+/DC8vL/V7PvzwQ7S0tMBRMcxotGXLFixdulRdmp2fn49hw4ZhwoQJuHHjRvdWiB5r8ODBqKioME+HDx82L1uyZAl27NiBrVu34uDBg+qRFq+88gpbtJvV19erv38Z7h9mzZo1+Pbbb7Fx40YcO3YM3t7ealuRO+x2MsicO3cOe/bswc6dO9WBeO7cuayVFesiyfDScfvZvHmzxXLWpfvJfZEMKkePHlV/7waDAePHj1f16uy+q7W1VQWZ5uZm5Obm4scff8SmTZvUlwaHJS/Npq4bM2aMWLhwofl9a2uriIyMFKtWrWJz9pAVK1aIYcOGPXRZTU2NcHV1FVu3bjXPO3/+vLwNgcjLy2ONrES2b05Ojvm90WgU4eHhYu3atRa1cXd3F5s3b1bv//rrL/VzJ06cMK/z+++/C51OJ65fv85aWaEu0ltvvSWmTp36yJ9hXXrGjRs3VH0OHjzY6X3Xrl27hJOTk6isrDSvs2HDBuHn5yeampqEI2LPjAYyDZ86dUp1l3d8/pN8n5eX151Zk55Anq6Q3egJCQnqW6TsepVkfeQ3no41kqegYmJiWKMedPnyZVRWVlrUQT5nRZ6Wbd9W5Ks8tTRq1CjzOnJ9uU3JnhyyngMHDqhTFMnJyZg/fz6qq6vNy1iXnnH37l31GhgY2Ol9l3xNTU1FWFiYeZ0JEyaoB1PKHk5HxDCjwa1bt1Q3X8c/JEm+lztu6hnygCi7Vnfv3o0NGzaoA+fYsWPVE1ZlHdzc3NRBkjXqPe3bw+O2FfkqD6gdubi4qJ07tyfrkaeYfvrpJ+zduxdffvmlOp0xadIktW9jXXqG0WjE4sWL8fzzz2PIkCHmdn/Svku+Pmybkhx1m3GIp2ZT3yR3vO2GDh2qwk1sbCx+/fVXNdCUiB7t1VdfNf9ffsuX29CAAQNUb824cePYdD1Ajp0pLCy0GOtH2rBnRoPg4GA4Ozs/MLpcvg8PD9dYCnpa8ptMUlISiouLVR3k6cCamhrWqBe1bw+P21bk698HzsurMuQVTtyeeo48VSv3bXL7YV2sb9GiRWqw+/79+xEVFWWe35l9l3x92DYlOeo2wzCjgewCfPbZZ1X3bMfuQvk+PT29O+tDXSAv5b106ZK6BFjWx9XV1aJG8vJHOaaGNeo58fHxaufasQ7yvL4cC9NeB/kqd9xyrEC7ffv2qW1K9rZRz7h27ZoaMyO3H9bFeuR4bBlkcnJy1N+53EY66sy+S77++eefFl8C9uzZo25tMGjQIDik3h6BbK9++eUXdUXGpk2b1Kj/uXPnCn9/f4vR5WRdy5YtEwcOHBCXL18WR44cEZmZmSI4OFhdHSDNmzdPxMTEiH379omTJ0+K9PR0NVH3qqurE6dPn1aT3KWsW7dO/f/KlStq+erVq9W2sX37dnH27Fl1BU18fLxoaGgw/46JEyeKESNGiGPHjonDhw+LxMREMWvWLJbKSnWRyz744AN1dYzcfv744w8xcuRI1e6NjY2sixXNnz9f9OvXT+27KioqzNP9+/fN6zxp39XS0iKGDBkixo8fLwoKCsTu3btFSEiIWL58uXBUDDNPYf369eoPzs3NTV2qffTo0e6rDD3RzJkzRUREhGr//v37q/fFxcXm5fJguWDBAhEQECC8vLzE9OnT1U6Dutf+/fvVwfLvk7z0t/3y7M8++0yEhYWpLwDjxo0Ter3e4ndUV1er8OLj46MuL3377bfVAZesUxd54JQHQnkAlJcBx8bGinfeeeeBL2OsS/d7WE3k9MMPP3Rp31VaWiomTZokPD091Zc4+eXOYDAIR6WT//R27xARERGRVhwzQ0RERHaNYYaIiIjsGsMMERER2TWGGSIiIrJrDDNERERk1xhmiIiIyK4xzBAREZFdY5ghIiIiu8YwQ0RERHaNYYaIiIjsGsMMERER2TWGGSIiIoI9+z8g5Lpesn6sPAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "\n",
    "nodes = ifr.get_nodes()\n",
    "nA, nB, nC, nD = nodes[4*nbragg], nodes[2*nbragg], nodes[0.0*nbragg], nodes[-2*nbragg]\n",
    "Ts = [0, 0, T, T, T + Tp, T + Tp, T + Tp + T, T + Tp + T]\n",
    "\n",
    "plt.figure()\n",
    "for nn in [nA, nB, nC, nD]:\n",
    "    for p in nn:\n",
    "        for node in nn[p]:\n",
    "            t, x = node.get_trajectory()\n",
    "            plt.plot(t, x)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5e27341-ce00-4943-8823-5eb55571628a",
   "metadata": {},
   "source": [
    "Create function for simulating the full interferometer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7e015dbc-a33c-45e5-adc1-3fc54e7dde85",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading single frequency Bragg precompute table, this could take a while...\n",
      "Precompute table loaded! Performing checks...\n",
      "Checks passed!\n",
      "Loading multifrequency Bragg precompute table, this could take a while...\n",
      "Precompute table loaded! Performing checks...\n",
      "Checks passed!\n"
     ]
    }
   ],
   "source": [
    "from mwave.integrate import propagate, omega_fnc_gaussian, multi_omega_fnc, phase_fnc_constant\n",
    "from mwave.precompute import load_precomputed_gbragg\n",
    "from mwave.utils.cloud import cloud_init\n",
    "\n",
    "Omega0 = 19.5\n",
    "w0 = 6.2e-3\n",
    "sigma = 0.259658916\n",
    "tau_factor = 5\n",
    "\n",
    "def deltalookup(vz):\n",
    "    return 4*nbragg + 4*(vz/0.0035) # 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",
    "def make_gbragg(sigma, use_precompute=True, tau_factor=5, tol=1e-8):\n",
    "    if use_precompute:\n",
    "        return load_precomputed_gbragg('../precomputation/gbragg_single_sig0.260.h5',\n",
    "                                       '../precomputation/gbragg_multi_sig0.260.h5',\n",
    "                                       table_sigma=sigma,\n",
    "                                       table_modulation_frequency=8*4)\n",
    "\n",
    "    tfinal = 2*tau_factor*sigma\n",
    "    t_center = tau_factor*sigma\n",
    "\n",
    "    def gbragg_fn(kvec, k_init, sigma, omegas, deltas, delta_phase=0.0, mod_freq=None, mod_phase=0.0):\n",
    "        natoms = len(omegas)\n",
    "        phi0_batch = np.zeros((natoms, len(kvec)), dtype=np.complex128)\n",
    "        phi0_batch[:, np.argmin(np.abs(kvec))] = 1.0  # start at k=0 frame\n",
    "\n",
    "        delta_eff = deltas - 4*k_init  # transform detuning to k=0 frame\n",
    "\n",
    "        if mod_freq is None:\n",
    "            result = propagate(\n",
    "                kvec, phi0_batch, tfinal, delta_eff,\n",
    "                omega_fnc_gaussian, np.array([1.0, sigma, t_center]),\n",
    "                phase_fnc_constant, np.array([0.0]),\n",
    "                omegas=omegas, tol=tol,\n",
    "            )\n",
    "        else:\n",
    "            result = propagate(\n",
    "                kvec, phi0_batch, tfinal, delta_eff,\n",
    "                multi_omega_fnc, np.array([1.0, sigma, t_center, mod_freq]),\n",
    "                phase_fnc_constant, np.array([0.0]),\n",
    "                omegas=omegas, tol=tol,\n",
    "            )\n",
    "\n",
    "        # Apply post-hoc laser phase in k=0 frame (same convention as load_precomputed_gbragg)\n",
    "        phi = result.phi_final * np.exp(-1j * np.outer(delta_phase, kvec/2))\n",
    "\n",
    "        # Shift output by Deltan = k_init//2 positions (same as load_precomputed_gbragg)\n",
    "        Deltan = k_init // 2\n",
    "        if Deltan == 0:\n",
    "            return phi\n",
    "        phi_out = np.zeros_like(phi)\n",
    "        if Deltan < 0:\n",
    "            phi_out[:, :Deltan] = phi[:, -Deltan:]\n",
    "        else:\n",
    "            phi_out[:, Deltan:] = phi[:, :-Deltan]\n",
    "        return phi_out\n",
    "\n",
    "    return gbragg_fn\n",
    "\n",
    "gbragg = make_gbragg(sigma, use_precompute=True)\n",
    "\n",
    "def cpops(x0, y0, z0, vz, vx, vy, cphase=0.0, injected_dphase=0.0):\n",
    "    bs_lookup_dict = {}\n",
    "\n",
    "    # Define beamsplitter function\n",
    "    def calc_bs(x0, y0, z0, vz, vx, vy, k_init, k_final, klattice, t, z, idx):\n",
    "\n",
    "        # Check if this is a multifrequency beamsplitter\n",
    "        multifrequency = isinstance(klattice, list)\n",
    "\n",
    "        # Load cached result\n",
    "        if idx in bs_lookup_dict:\n",
    "            if k_init in bs_lookup_dict[idx]:\n",
    "                if int(k_final) in bs_lookup_dict[idx][k_init]:\n",
    "                    return bs_lookup_dict[idx][k_init][int(k_final)]\n",
    "\n",
    "        # Compute omegas and deltas\n",
    "        x = x0 + vx*t\n",
    "        y = y0 + vy*t\n",
    "        omegas = omegalookup(x, y, z0 + z)\n",
    "        deltas = deltalookup(vz)\n",
    "\n",
    "        # Set sigma\n",
    "        sigma = 0.259658916\n",
    "\n",
    "        # Compute phases\n",
    "        phases = deltas*t\n",
    "\n",
    "        # Apply common mode phase if provided and at pulse 3\n",
    "        if idx == 3:\n",
    "            phases += cphase\n",
    "\n",
    "            # Apply differential phase if provided at pulse 3. The way I am injecting this is unphysical/different from how we do in the experiment. But I'm not sure if there is a better way.\n",
    "            if k_init == 0.0 and k_final == -2.0*nbragg:\n",
    "                phases += injected_dphase\n",
    "\n",
    "        # Compute effect of Bragg beamsplitter\n",
    "        if not multifrequency:\n",
    "            phi = gbragg(ifr.kvec, int(k_init), sigma, omegas, deltas, delta_phase=phases)\n",
    "        else:\n",
    "            phi = gbragg(ifr.kvec, int(k_init), sigma, omegas, deltas, delta_phase=phases, mod_freq=8*4, mod_phase=0.0)\n",
    "\n",
    "        # Determine index of k_final state\n",
    "        kf_idx = np.argmin(np.abs(ifr.kvec - k_final))\n",
    "\n",
    "        # Save wavefunction to cache\n",
    "        if idx not in bs_lookup_dict:\n",
    "            bs_lookup_dict[idx] = {}\n",
    "        if k_init not in bs_lookup_dict[idx]:\n",
    "            bs_lookup_dict[idx][k_init] = {}\n",
    "        if int(k_final) not in bs_lookup_dict[idx][k_init]:\n",
    "            bs_lookup_dict[idx][k_init][int(k_final)] = phi[:,kf_idx]\n",
    "        else:\n",
    "            raise RuntimeError('Array should not have been created but it was!')\n",
    "\n",
    "        # Return\n",
    "        return phi[:,kf_idx]\n",
    "\n",
    "    # Define free evolution function\n",
    "    propfnc = lambda x0, y0, z0, vz, vx, vy, t, k: np.exp(-1j*t*k**2)\n",
    "\n",
    "    # Apply functions to interferometer\n",
    "    fnclst = [lambda x0, y0, z0, vz, vx, vy, k_init, k_final, klattice, t, x: calc_bs(x0, y0, z0, vz, vx, vy, k_init, k_final, klattice, t, x, 1),\n",
    "              propfnc,\n",
    "              lambda x0, y0, z0, vz, vx, vy, k_init, k_final, klattice, t, x: calc_bs(x0, y0, z0, vz, vx, vy, k_init, k_final, klattice, t, x, 2),\n",
    "              propfnc,\n",
    "              lambda x0, y0, z0, vz, vx, vy, k_init, k_final, klattice, t, x: calc_bs(x0, y0, z0, vz, vx, vy, k_init, k_final, klattice, t, x, 3),\n",
    "              propfnc,\n",
    "              lambda x0, y0, z0, vz, vx, vy, k_init, k_final, klattice, t, x: calc_bs(x0, y0, z0, vz, vx, vy, k_init, k_final, klattice, t, x, 4)]\n",
    "    ifr.set_operation_funcs(fnclst)\n",
    "\n",
    "    # Load population functions\n",
    "    popfnc = ifr.get_population_func([4*nbragg, 2*nbragg, 0*nbragg, -2*nbragg], lambda x1, x2, x3, x4, x5, x6: np.zeros_like(x0), lambda x1, x2, x3, x4, x5, x6: np.ones_like(x0,dtype=np.complex128), lambda x1, x2, x3, x4, x5, x6: np.zeros_like(x0,dtype=np.complex128))\n",
    "\n",
    "    # Evaluate populations and return\n",
    "    return popfnc(4*nbragg, [x0, y0, z0, vz, vx, vy]), popfnc(2*nbragg, [x0, y0, z0, vz, vx, vy]), popfnc(0*nbragg, [x0, y0, z0, vz, vx, vy]), popfnc(-2*nbragg, [x0, y0, z0, vz, vx, vy])\n",
    "\n",
    "x0, y0, z0, vz, vx, vy = cloud_init(natoms=1000, sigma_cloud=1e-3, sigma_transverse_v=3.5e-3, sigma_vertical_v=0.1*3.5e-3)\n",
    "pA, pB, pC, pD = cpops(x0, y0, z0, vz, vx, vy, cphase=np.pi/4)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "397b42b8-5ff5-4b93-a17a-a2377c959ea2",
   "metadata": {},
   "source": [
    "Make ellipse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f2910c10-5690-4d4d-9db8-c827a6a048cf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAckAAAGiCAYAAACMDD3oAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAK2JJREFUeJzt3Q9UVWW+//EvoII6gRoqUij+mTG9/g0TqabpJkv8MzO68k7SnzG5pjdTG/M/3dTUJtIcb9eJGcubqWsyzX7qWNOQptM0U6SGY2YXXWkaaoKSAYGGCvu3vs9d58w5wIMgHODA+7XWXoe9z3P2ec5mcz48ez/P3gGO4zgCAADKCSy/CAAAEJIAAFSCliQAABaEJAAAFoQkAAAWhCQAABaEJAAAFoQkAAAWhCQAABaEJAAA9RGSH3zwgfzsZz+TyMhICQgIkO3bt1/zNe+//77ceuutEhwcLD169JB169aVK5OamirR0dESEhIisbGxsm/fPh99AgBAU+bTkCwqKpL+/fubUKuKEydOyKhRo+Rf//Vf5eDBgzJjxgx55JFH5N1333WX2bx5s8ycOVMWLVokBw4cMOtPSEiQc+fO+fCTAACaooC6usC5tiS3bdsmY8aMsZaZN2+e/OlPf5LDhw+7lyUmJkpeXp6kpaWZeW053nbbbfLiiy+a+dLSUomKipLp06fL/Pnz6+CTAACaimbSgKSnp0t8fLzXMm0laotSXb58WTIyMiQ5Odn9fGBgoHmNvtamuLjYTC4arBcuXJAbb7zRhDcAwL84jiPfffedOZ2nOdAkQjI7O1s6duzotUznCwoK5NKlS/Ltt99KSUlJhWWOHDliXW9KSoosXrzYZ/UGANSPU6dOyc0339w0QtJXtOWp5zFd8vPzpXPnzmbjhoaG1mvdAADVp40nPdV2ww03iC81qJCMiIiQnJwcr2U6r0HWsmVLCQoKMlNFZfS1NtpTVqeydL2EJAD4rwAfnzJrUOMk4+LiZPfu3V7Ldu3aZZarFi1aSExMjFcZPb+o864yAAD4RUgWFhaaoRw6uYZ46M9ZWVnuw6Djx493l3/00Uflyy+/lLlz55pzjL/73e/kjTfekCeeeMJdRg+brlmzRtavXy+ZmZkyZcoUM9QkKSnJlx8FANAE+fRw6yeffGLGPLq4zgs+/PDD5iIBZ8+edQem6tq1qxkCoqH43//93+Zk7P/8z/+YHq4u48aNk/Pnz8vChQtNR58BAwaY4SFlO/MAAOA34yQb2gnfsLAw04GHc5IA4H8K6uh7vEGdkwQAoCEhJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAqM+QTE1NlejoaAkJCZHY2FjZt2+ftezdd98tAQEB5aZRo0a5y0yYMKHc88OHD6+LjwIAaEKa+foNNm/eLDNnzpTVq1ebgHzhhRckISFBjh49Kh06dChXfuvWrXL58mX3/DfffCP9+/eXX/ziF17lNBRfffVV93xwcLCPPwkAoKnxeUty5cqVMmnSJElKSpLevXubsGzVqpWsXbu2wvLt2rWTiIgI97Rr1y5TvmxIaih6lmvbtq2vPwoAoInxaUhqizAjI0Pi4+P/+YaBgWY+PT29Sut45ZVXJDExUVq3bu21/P333zct0Z49e8qUKVNMi9OmuLhYCgoKvCYAAOo1JHNzc6WkpEQ6duzotVzns7Ozr/l6PXd5+PBheeSRR8odat2wYYPs3r1bli1bJn/9619lxIgR5r0qkpKSImFhYe4pKiqqhp8MANAU+PycZE1oK7Jv374yePBgr+XasnTR5/v16yfdu3c3rcuhQ4eWW09ycrI5L+qiLUmCEgBQry3J8PBwCQoKkpycHK/lOq/nEStTVFQkmzZtkokTJ17zfbp162be69ixYxU+r+cvQ0NDvSYAAOo1JFu0aCExMTHmsKhLaWmpmY+Li6v0tVu2bDHnEh966KFrvs/p06fNOclOnTrVSr0BAKiT3q16mHPNmjWyfv16yczMNJ1stJWovV3V+PHjzeHQig61jhkzRm688Uav5YWFhTJnzhz5+OOP5eTJkyZwR48eLT169DBDSwAA8JtzkuPGjZPz58/LwoULTWedAQMGSFpamrszT1ZWlunx6knHUP7973+XnTt3llufHr49dOiQCd28vDyJjIyUYcOGydKlSxkrCQCoVQGO4zjSxGjHHe3lmp+fz/lJAPBDBXX0Pc61WwEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAR86m39JPjqeax4B+J86CcnU1FSJjo6WkJAQiY2NlX379lnLrlu3TgICArwmfZ0nx3Fk4cKF0qlTJ2nZsqXEx8fLF198UQefBKi6zfuz5I7n9sgDa/aaR50H4F98HpKbN2+WmTNnyqJFi+TAgQPSv39/SUhIkHPnzllfExoaKmfPnnVPX331ldfzy5cvl1WrVsnq1atl79690rp1a7PO77//3tcfB6gSbTkmb/1MSp3/m9fHJ7cepkUJ+Bmfh+TKlStl0qRJkpSUJL179zbB1qpVK1m7dq31Ndp6jIiIcE8dO3b0akW+8MIL8tRTT8no0aOlX79+smHDBvn6669l+/btvv44QJWcyC1yB6RLiePIydyLbEHAj/g0JC9fviwZGRnmcKj7DQMDzXx6err1dYWFhdKlSxeJiooyQfj555+7nztx4oRkZ2d7rTMsLMwcxrWts7i4WAoKCrwmwJe6hreWwADvZUEBARId3sovNjznUoE6CMnc3FwpKSnxagkqndegq0jPnj1NK/OPf/yj/OEPf5DS0lK5/fbb5fTp0+Z51+uqs86UlBQTpK5JwxfwpU5hLSXl3r4mGJU+PntvH7O8oeNcKvBPzaSBiYuLM5OLBmSvXr3kpZdekqVLl17XOpOTk815URdtSRKU8LVxt3WWu37U3hxi1RakPwSk7Vyqfg5/qD/gVyEZHh4uQUFBkpOT47Vc5/VcY1U0b95cBg4cKMeOHTPzrtfpOrR3q+c6BwwYUOE6goODzQTUNQ0WfwqXys6l+tPnAPzicGuLFi0kJiZGdu/e7V6mh0913rO1WBk9XPvZZ5+5A7Fr164mKD3XqS1D7eVa1XUCaJznUgG/692qhznXrFkj69evl8zMTJkyZYoUFRWZ3q5q/Pjx5nCoy5IlS2Tnzp3y5ZdfmiEjDz30kBkC8sgjj7h7vs6YMUOeeeYZ2bFjhwlQXUdkZKSMGTPG1x8HaNT8+Vwq4JfnJMeNGyfnz583g/+1Y40eEk1LS3N3vMnKyjI9Xl2+/fZbM2REy7Zt29a0RD/66CMzfMRl7ty5JmgnT54seXl5cuedd5p1lr3oAICmcS4V8JUARwceNjF6eFZ7uebn55sLFwAA/EtBHX2Pc+1WAAAsCEkAACwISQAALAhJ+B0umQagyV5xB7jWJdNcV4TR8Xw6XEF7YwKAL9CShN/g9lMA6hohCb/B7acA1DVCEn6DS6YBqGuEJPwGl0wDUNfouAO/wiXTANQlQhJ+x99uPwXAf3G4FQAAC0ISAAALQhIAAAtCEgAAC0ISQIPGtXpRn+jdCqDB4lq9qG+0JAE0SFyrFw0BIQmgQeJavWgICEkADRLX6kVDQEgCaJC4Vi8aAjruAGiwuFYv6hshCaBB41q9qE8cbgUAwIKQBADAgpAEAMCCkAQAwIKQBADAgpAEAMCCkAQAwIKQBADAgpAEAMCCkESt4ga5ABoTLkuHWsMNcgE0NrQkUSu4QS6AxoiQRK3gBrkAGqM6CcnU1FSJjo6WkJAQiY2NlX379lnLrlmzRn784x9L27ZtzRQfH1+u/IQJEyQgIMBrGj58eB18Ethwg1wAjZHPQ3Lz5s0yc+ZMWbRokRw4cED69+8vCQkJcu7cuQrLv//++3L//ffLX/7yF0lPT5eoqCgZNmyYnDlzxquchuLZs2fd0+uvv+7rj4JKcINcAI1RgOM4ji/fQFuOt912m7z44otmvrS01ATf9OnTZf78+dd8fUlJiWlR6uvHjx/vbknm5eXJ9u3bq1SH4uJiM7kUFBSYOuTn50toaOh1fzZUfG7yZO5FiQ5vZYITAHxBv8fDwsJ8/j3u05bk5cuXJSMjwxwydb9hYKCZ11ZiVVy8eFGuXLki7dq1K9fi7NChg/Ts2VOmTJki33zzjXUdKSkpZmO6Jg1I+IYGY1z3GwlIAI2CT0MyNzfXtAQ7duzotVzns7Ozq7SOefPmSWRkpFfQ6qHWDRs2yO7du2XZsmXy17/+VUaMGGHeqyLJycnmvw3XdOrUqRp+MgBAU9Cgx0k+99xzsmnTJtNq1E4/LomJie6f+/btK/369ZPu3bubckOHDi23nuDgYDMBANBgWpLh4eESFBQkOTk5Xst1PiIiotLXrlixwoTkzp07TQhWplu3bua9jh07Viv1BgDA5yHZokULiYmJMYdFXbTjjs7HxcVZX7d8+XJZunSppKWlyaBBg675PqdPnzbnJDt16lRrdQcAwOdDQHT4h459XL9+vWRmZppONkVFRZKUlGSe1x6res7QRc8xLliwQNauXWvGVuq5S50KCwvN8/o4Z84c+fjjj+XkyZMmcEePHi09evQwQ0sAAPCbc5Ljxo2T8+fPy8KFC03YDRgwwLQQXZ15srKyTI9Xl9///vemV+y//du/ea1Hx1k+/fTT5vDtoUOHTOjqMBDt1KPjKLXlyXlHAIBfjZNsyuNrAAC+0SjGSQIA4M8ISQCNyqenvpU1fztuHoFGPU4SAKpj1hsH5f8d+Od1nsfeepP85r4BbERcN1qSABoFbTl6BqTSeVqUqAlCEkCjsO/khQqXf3KSw664foQkgEZhcLT3TRBcBkW3rfO6oPEgJAE0Cv2j2ppzkJ50XpcD14uOOwAaDe2kMz6uiznEqi1IAhI1RUgCaFQ0GAlH1BYOtwIAYEFIAgBgQUgCAGBBSAIAYEFIAgBgQUgCAGBBSAIAYEFIAgBgQUgCAGBBSAIAYEFIAgBgQUii1pzNvyQfHc81jwDQGHCBc9SKzfuzJHnrZ1LqiAQGiKTc21fG3daZrQvAr9GSRI1py9EVkEofn9x6mBYlAL9HSKLGTuQWuQPSpcRx5GTuRbYuAL9GSKLGuoa3NodYPQUFBEh0eCu2LgC/RkiixjqFtTTnIDUYlT4+e28fsxwA/Bkdd1ArtJPOXT9qbw6xaguSgATQGBCSqDUajIQjgMaEw60AAFgQkgAAWBCSAABYEJIAAFgQkgAAWBCSAABYEJIAGiTuKoMmE5KpqakSHR0tISEhEhsbK/v27au0/JYtW+SWW24x5fv27SvvvPOO1/OO48jChQulU6dO0rJlS4mPj5cvvvjCx58CQF3eVeaO5/bIA2v2mkedBxplSG7evFlmzpwpixYtkgMHDkj//v0lISFBzp07V2H5jz76SO6//36ZOHGi/OMf/5AxY8aY6fDhw+4yy5cvl1WrVsnq1atl79690rp1a7PO77//3tcfB4CPcVcZNCQBjjbLfEhbjrfddpu8+OKLZr60tFSioqJk+vTpMn/+/HLlx40bJ0VFRfL222+7lw0ZMkQGDBhgQlGrGxkZKbNmzZLZs2eb5/Pz86Vjx46ybt06SUxMLLfO4uJiM7kUFBSYOujrQkNDffTJAVwPvXG3tiDLen3SEInrfiMbFe7v8bCwMJ9/j/u0JXn58mXJyMgwh0PdbxgYaObT09MrfI0u9yyvtJXoKn/ixAnJzs72KqMbSsPYts6UlBRTxjVpQAJomLirDBoSn4Zkbm6ulJSUmFaeJ53XoKuILq+svOuxOutMTk42/224plOnTtXocwHwHe4qg4akSVzgPDg42EwA/AN3lUGTaEmGh4dLUFCQ5OTkeC3X+YiIiApfo8srK+96rM46Afhni1LPQXJnGTTakGzRooXExMTI7t273cu0447Ox8XFVfgaXe5ZXu3atctdvmvXriYMPcvoCVzt5WpbJwAADfJwqw7/ePjhh2XQoEEyePBgeeGFF0zv1aSkJPP8+PHj5aabbjKda9SvfvUr+clPfiK/+c1vZNSoUbJp0yb55JNP5OWXXzbPBwQEyIwZM+SZZ56RH/7whyY0FyxYYHq86lARAAD8JiR1SMf58+fN4H/tWKNDOdLS0twdb7KyskyPV5fbb79dNm7cKE899ZQ8+eSTJgi3b98uffr0cZeZO3euCdrJkydLXl6e3HnnnWadevEBAAD8ZpxkUx5fAwDwjUYxThIAAH9GSAIAYEFIwm9wVwgAda1JXEwA/k/vApG89TMpdUQCA0RS7u1rBpwDgC/RkkSDx10hANQXQhIN3oncItOC9FTiOHIy92J9VQlAE0FIosHjrhAA6gshiQaPu0IAqC903IFf4K4QAOoDIQm/alFyRwgAdYnDrQAAWBCSANy4YAPgjcOtAAwu2ACUR0sSABdsACwISQBcsAGwICQBcMEGwIKQBMAFGwALOu4APughqteb1cvp+dO4Ti7YAJRHSAK1yN97iHLBBsAbh1uBWsItvYDGh5AEagm39AIaH0ISqCXc0gtofAhJoJZwSy+g8aHjDlCL6CEKNC6EJFDL6CEKNB4cbgUAwIKQBADAgpAEAMCCkAQAwIKQBADAgpAEAMCCkAQAwIKQBADAgpAEAMCCkAQAoD5C8sKFC/Lggw9KaGiotGnTRiZOnCiFhYWVlp8+fbr07NlTWrZsKZ07d5bHH39c8vPzvcoFBASUmzZt2uTLjwIAaIJ8eu1WDcizZ8/Krl275MqVK5KUlCSTJ0+WjRs3Vlj+66+/NtOKFSukd+/e8tVXX8mjjz5qlr355pteZV999VUZPny4e15DGACA2hTgOI4jPpCZmWmCbv/+/TJo0CCzLC0tTUaOHCmnT5+WyMjIKq1ny5Yt8tBDD0lRUZE0a/Z/ma4tx23btsmYMWOuq24FBQUSFhZmWqjaygUA+JeCOvoe99nh1vT0dNO6cwWkio+Pl8DAQNm7d2+V1+PaAK6AdJk6daqEh4fL4MGDZe3atVJZ1hcXF5sN6jkBAFBvh1uzs7OlQ4cO3m/WrJm0a9fOPFcVubm5snTpUnOI1tOSJUvknnvukVatWsnOnTvlscceM+c69fxlRVJSUmTx4sU1+DQAgKao2i3J+fPnV9hxxnM6cuRIjSumrb1Ro0aZQ7ZPP/2013MLFiyQO+64QwYOHCjz5s2TuXPnyvPPP29dV3JysmmRuqZTp07VuH4AgMav2i3JWbNmyYQJEyot061bN4mIiJBz5855Lb969arpwarPVea7774znXJuuOEGc+6xefPmlZaPjY01LU49rBocHFzueV1W0XIAAGo1JNu3b2+ma4mLi5O8vDzJyMiQmJgYs2zPnj1SWlpqQq2yFmRCQoIJtR07dkhISMg13+vgwYPStm1bghAA4B/nJHv16mVag5MmTZLVq1ebISDTpk2TxMREd8/WM2fOyNChQ2XDhg2mA44G5LBhw+TixYvyhz/8wauTjQZzUFCQvPXWW5KTkyNDhgwxAarDS5599lmZPXu2rz4KAKCJ8uk4yddee80Eowah9modO3asrFq1yv28BufRo0dNKKoDBw64e7726NHDa10nTpyQ6Ohoc+g1NTVVnnjiCdOjVcutXLnShDEAAH4xTrIhY5wkAPi3An8fJwkAgL8jJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAqI+QvHDhgjz44IMSGhoqbdq0kYkTJ0phYWGlr7n77rslICDAa3r00Ue9ymRlZcmoUaOkVatW0qFDB5kzZ45cvXrVlx8FANAENfPlyjUgz549K7t27ZIrV65IUlKSTJ48WTZu3Fjp6yZNmiRLlixxz2sYupSUlJiAjIiIkI8++sisf/z48dK8eXN59tlnfflxAABNTIDjOI4vVpyZmSm9e/eW/fv3y6BBg8yytLQ0GTlypJw+fVoiIyOtLckBAwbICy+8UOHzf/7zn+WnP/2pfP3119KxY0ezbPXq1TJv3jw5f/68tGjR4pp1KygokLCwMMnPzzetXACAfymoo+9xnx1uTU9PN4dYXQGp4uPjJTAwUPbu3Vvpa1977TUJDw+XPn36SHJysly8eNFrvX379nUHpEpISDAb7PPPP69wfcXFxeZ5zwkAgHo73JqdnW3OF3q9WbNm0q5dO/OczQMPPCBdunQxLc1Dhw6ZFuLRo0dl69at7vV6BqRyzdvWm5KSIosXL66FTwUAaEqqHZLz58+XZcuWXfNQ6/XSc5Yu2mLs1KmTDB06VI4fPy7du3e/rnVqa3TmzJnueW1JRkVFXXcdAQBNQ7VDctasWTJhwoRKy3Tr1s10rDl37pzXcu2Bqj1e9bmqio2NNY/Hjh0zIamv3bdvn1eZnJwc82hbb3BwsJkAAPBpSLZv395M1xIXFyd5eXmSkZEhMTExZtmePXuktLTUHXxVcfDgQfOoLUrXen/961+bAHYdztXes3riVjsKAQDQ4Dvu9OrVS4YPH26Gc2jL78MPP5Rp06ZJYmKiu2frmTNn5JZbbnG3DPWQ6tKlS02wnjx5Unbs2GGGd9x1113Sr18/U2bYsGEmDH/5y1/Kp59+Ku+++6489dRTMnXqVFqLAAD/uZiA9lLVENRzijr0484775SXX37Z/byOndROOa7eqzp847333jNBqK/TQ7tjx46Vt956y/2aoKAgefvtt82jtiofeughE6Se4yoBAGjQ4yQbMsZJAoB/K/D3cZIAAPg7QhIAAAtCEgAAC0ISAAALQhIAAAtCEgAAC0ISAAALQhIAAAtCEgAAC0ISAAALQhIAAAtCEgAAC0ISAAALQhIAAAtCEgAAC0ISAAALQhIAAAtCEgAAC0ISAAALQhIAAAtCEgAAC0ISAAALQhIAAAtCEgAAC0ISAAALQhIAAAtCEgAAC0ISAAALQhIAAAtCEgAAC0ISAAALQhIAAAtCEgAAC0ISAAALQhIAgPoIyQsXLsiDDz4ooaGh0qZNG5k4caIUFhZay588eVICAgIqnLZs2eIuV9HzmzZt8uVHAQA0Qc18uXINyLNnz8quXbvkypUrkpSUJJMnT5aNGzdWWD4qKsqU9/Tyyy/L888/LyNGjPBa/uqrr8rw4cPd8xrCAAD4RUhmZmZKWlqa7N+/XwYNGmSW/fa3v5WRI0fKihUrJDIystxrgoKCJCIiwmvZtm3b5L777pMf/OAHXss1FMuWBQDALw63pqenmyBzBaSKj4+XwMBA2bt3b5XWkZGRIQcPHjSHacuaOnWqhIeHy+DBg2Xt2rXiOI51PcXFxVJQUOA1AQBQby3J7Oxs6dChg/ebNWsm7dq1M89VxSuvvCK9evWS22+/3Wv5kiVL5J577pFWrVrJzp075bHHHjPnOh9//PEK15OSkiKLFy+uwacBADRF1W5Jzp8/39q5xjUdOXKkxhW7dOmSOXdZUStywYIFcscdd8jAgQNl3rx5MnfuXHPe0iY5OVny8/Pd06lTp2pcPwBA41ftluSsWbNkwoQJlZbp1q2bOV947tw5r+VXr141PV6rci7xzTfflIsXL8r48eOvWTY2NlaWLl1qDqsGBweXe16XVbQcAIBaDcn27dub6Vri4uIkLy/PnFeMiYkxy/bs2SOlpaUm1KpyqPXnP/95ld5Lz1u2bduWIAQA+Mc5ST2XqEM0Jk2aJKtXrzZDQKZNmyaJiYnunq1nzpyRoUOHyoYNG0wHHJdjx47JBx98IO+880659b711luSk5MjQ4YMkZCQEDO85Nlnn5XZs2f76qMAAJoon46TfO2110wwahBqr9axY8fKqlWr3M9rcB49etQcVvWkvVVvvvlmGTZsWLl1Nm/eXFJTU+WJJ54wPVp79OghK1euNGEMAEBtCnAqGzvRSOkQkLCwMNOJR68GBADwLwV19D3OtVsBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCwICQBALAgJAEAsCAkAQCo65D89a9/Lbfffru0atVK2rRpU6XXOI4jCxculE6dOknLli0lPj5evvjiC68yFy5ckAcffFBCQ0PNeidOnCiFhYU++hQAgKbMZyF5+fJl+cUvfiFTpkyp8muWL18uq1atktWrV8vevXuldevWkpCQIN9//727jAbk559/Lrt27ZK3335bPvjgA5k8ebKPPgUAoCkLcLT55kPr1q2TGTNmSF5eXqXltBqRkZEya9YsmT17tlmWn58vHTt2NOtITEyUzMxM6d27t+zfv18GDRpkyqSlpcnIkSPl9OnT5vUVKS4uNpOLrrdz585y6tQp0yIFAPiXgoICiYqKMtkSFhbmuzdyfOzVV191wsLCrlnu+PHjGtbOP/7xD6/ld911l/P444+bn1955RWnTZs2Xs9fuXLFCQoKcrZu3Wpd96JFi8y6mdgG7APsA+wDjWsfOH78uONLzaSByM7ONo/acvSk867n9LFDhw5ezzdr1kzatWvnLlOR5ORkmTlzpnte//Po0qWLZGVl+fY/EB/95+RvLWB/rbc/1516s70b837ieURQv/99qVohOX/+fFm2bFmlZfSQ6C233CINSXBwsJnK0oD0tx1DaZ2pN9ucfaXh4W+z7gUGBjackNTzhRMmTKi0TLdu3a6rIhEREeYxJyfH9G510fkBAwa4y5w7d87rdVevXjU9Xl2vBwCgXkKyffv2ZvKFrl27mqDbvXu3OxT1UID2cnX1kI2LizOHSjMyMiQmJsYs27Nnj5SWlkpsbKxP6gUAaLp81k7V830HDx40jyUlJeZnnTzHNOph2W3btpmfAwICTC/YZ555Rnbs2CGfffaZjB8/3vRYHTNmjCnTq1cvGT58uEyaNEn27dsnH374oUybNs30fLX1bK2IHnpdtGhRhYdgGzLqzTZnX2mY+NtsxNvcVz2CHn744Qp7Iv3lL39xl9F57f3qUlpa6ixYsMDp2LGjExwc7AwdOtQ5evSo13q/+eYb5/7773d+8IMfOKGhoU5SUpLz3Xff+epjAACaMJ+PkwQAwF9x7VYAACwISQAALAhJAAAsCEkAAJpSSPrzbbqq+x4nT540w2cqmrZs2eIuV9HzmzZtqrd6q7vvvrtcnR599FGvMjqEaNSoUeZ3qZcknDNnjrmARH3VW8tPnz5devbsafYTvSzW448/bi6R5am2t3dqaqpER0dLSEiIGROsQ6Aqo797HWKl5fv27SvvvPNOtff32lKduq9Zs0Z+/OMfS9u2bc2k9SpbXi9oUnbb6tCw+qy33oShbJ30dfWxzatT74r+BnXSv7m63N4ffPCB/OxnPzND+XT927dvv+Zr3n//fbn11lvNEJAePXqY30FN/24q5DRCCxcudFauXOnMnDmzShdXV88995wpu337dufTTz91fv7znztdu3Z1Ll265C4zfPhwp3///s7HH3/s/O1vf3N69OhhhqPUpuq+x9WrV52zZ896TYsXLzZDZDyHxriG23iW8/xsdV1v9ZOf/MSZNGmSV53y8/O9PlufPn2c+Ph4c+H7d955xwkPD3eSk5Prrd6fffaZc++99zo7duxwjh075uzevdv54Q9/6IwdO9arXG1u702bNjktWrRw1q5d63z++edmm+mF/nNycios/+GHH5qL/i9fvtz53//9X+epp55ymjdvbupenf29NlS37g888ICTmppqft+ZmZnOhAkTTD1Pnz7tNbxMf2+e2/bChQv1Wm/9XeuQNM86ZWdne5Wpi21e3XrrkDrPOh8+fNjsO55D8x6ug+2tf9v/+Z//aW5UoX8727Ztq7T8l19+6bRq1cp8x+s+/tvf/tbUOy0t7bq3hU2jDMnq3oFEx2dGREQ4zz//vHtZXl6eGav5+uuvm3n9Regvb//+/e4yf/7zn52AgADnzJkztVLf2nqPAQMGOP/+7//utawqO15d11tD8le/+lWlfziBgYFeXza///3vzZdRcXFxvdW7rDfeeMP8MeodaXyxvQcPHuxMnTrVPV9SUuJERkY6KSkpFZa/7777nFGjRnkti42Ndf7jP/6jyvt7balu3cvSf5RuuOEGZ/369V5f2qNHj3Z8qbr1vtZ3TV1t85pu7//6r/8y27uwsLBOt7enqvztzJ071/mXf/kXr2Xjxo1zEhISam1buDTKw63VdeLECXMXET384Xnxc22ep6enm3l91MNxrvtYKi2vF9fVS+fVhtp4D71kn17ZSA8bljV16lQJDw+XwYMHy9q1a83hn/qu92uvvWbq1KdPH3O3losXL3qtVw8Vet4ZRm/CrZcr1Btv12e9PemhVj1cq3ekqe3trTcv19+p576p9dN5175Z0efyLO/abq7yVdnfa8P11L0s3R+uXLlS7k4PeqhND7/rYW+9bOU333xT7/XWw/R6dyG9q8bo0aO99tG62Oa1sb1feeUVcwUzveF9XW3v63Gtfbw2toVLg7lVVmO9TVd161HT99CdXC/fp+dkPS1ZskTuuecec25v586d8thjj5k/aj2fVl/1fuCBB8yXip6HOHTokMybN0+OHj0qW7duda+3ot+J67n6qren3NxcWbp0qUyePNkn21vXr5d1rGg7HDlyxPq5rrUvu5bZytSG66l7WbpP6P7h+WWn58Puvfdec73n48ePy5NPPikjRowwX35BQUH1Um8ND/1HqF+/fuafphUrVpi/QQ3Km2++uU62eU23t56vO3z4sPkO8TTcx9v7etj2cf0H+tKlS/Ltt9/WeN/zu5D019t0VafuNaU7x8aNG2XBggXlnvNcNnDgQCkqKpLnn3++0i9tX9fbM1i0xagdGoYOHWr+ELt3797gt7f+QWoHh969e8vTTz9d4+0Nb88995zp7KStGM9OMNrS8dxvNJh0f9Fyuv/UB735gk4uGpD6z+pLL71k/onyBxqOuj31yIenxAa4veuS34SkP9+mq6p1r+mtwN58801zeEovDH8tephH/3iLi4utFwiuq3p71kkdO3bM/BHqa8v2RtPfiapsvXVR7++++878h33DDTeYi/Q3b968xtu7Inq4Vv9bd31uF5231VGXV1a+Kvt7bbieurtoS0xD8r333jNfytf6Xep76X5TG1/aNam3i+4P+s+R1qmutnlN6q3/xOk/JHoE5Fq61fL2vh62fVxPe2jPYd0ONf0dujmNWHU77qxYscK9THtZVtRx55NPPnGXeffdd33Sced630M7wpTtZWnzzDPPOG3btnVqQ21tm7///e9mPdrzz7PjjmdvtJdeesl03Pn+++/rrd66bwwZMsRs76KiIp9vb+2AMG3aNK8OCDfddFOlHXd++tOfei2Li4sr13Gnsv29tlS37mrZsmXmd5yenl6l9zh16pT5nf3xj3906rPeZTsc9ezZ03niiSfqdJtfb731u1LrkpubWy/b+3o67mjPd0/aK71sx52a/A7d9XEaoa+++sp0IXcNhdCfdfIcEqE7sHY39uyerd2D9Rd/6NAh05uroiEgAwcOdPbu3Wu+0LXrvy+GgFT2HtoVXuuuz3v64osvzI6rvTPL0uEKa9asMUMAtNzvfvc7031ah8rUV711+MSSJUtMQJ04ccJs927dujl33XVXuSEgw4YNcw4ePGi6d7dv377Wh4BUp976xaY9Rfv27Ws+g2e3eK2vL7a3dmXXL7B169aZYJ88ebLZV129fn/5y1868+fP9xoC0qxZM/OFrMMoFi1aVOEQkGvt77WhunXXemlP4TfffNNr27r+dvVx9uzZJkB1v3nvvfecW2+91fzeauMfp+utt37X6D9Yx48fdzIyMpzExEQnJCTEDD2oy21e3Xq73HnnnaZ3aFnf1dH21vdxfU9rSOoQPv1Zv8uV1lnrXnYIyJw5c8w+rsOGKhoCUtm2aNIh6c+36brWe+iOWvazKA2OqKgo899SWRqcOixE19m6dWszLnD16tUVlq2remdlZZlAbNeundneOj5Rd3jPcZLq5MmTzogRI5yWLVuaMZKzZs3yGmpR1/XWx4r2LZ20rK+2t44D69y5swkQ/Q9Zx3W6aItW9/myw1J+9KMfmfLaVf5Pf/qT1/NV2d9rS3Xq3qVLlwq3rQa9unjxovmnSf9Z0uDX8jr+rbpffLVd7xkzZrjL6jYdOXKkc+DAgXrZ5tXdV44cOWK28c6dO8ut62IdbW/b35WrrvqodS/7Gv0708+p/2B7fp9XZVtUFbfKAgDAgnGSAABYEJIAAFgQkgAAWBCSAABYEJIAABCSAABUDy1JAAAsCEkAACwISQAALAhJAAAsCEkAAKRi/x/E0UecJOYARwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def calc_xy(a, b, c, d):\n",
    "    return (a-b)/(a+b), (c-d)/(c+d)\n",
    "\n",
    "x0, y0, z0, vz, vx, vy = cloud_init(natoms=1000, sigma_cloud=1e-3, sigma_transverse_v=0.1*3.5e-3/omega_r_conversion, sigma_vertical_v=0.1*3.5e-3)\n",
    "\n",
    "cphases = np.linspace(0, 2*np.pi, 12)\n",
    "x, y = np.full_like(cphases, np.nan), np.full_like(cphases, np.nan)\n",
    "for i in range(len(cphases)):\n",
    "    pA, pB, pC, pD = cpops(x0, y0, z0, vz, vx, vy, cphase=cphases[i], injected_dphase=np.pi/8)\n",
    "    x[i], y[i] = calc_xy(np.sum(pA), np.sum(pB), np.sum(pC), np.sum(pD))\n",
    "\n",
    "plt.plot(x, y, '.')\n",
    "plt.gca().set_aspect('equal')\n",
    "plt.xlim(-1,1)\n",
    "plt.ylim(-1,1)\n",
    "plt.show()"
   ]
  }
 ],
 "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.12.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}