Generate_configs.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9ddbb7d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['figure.figsize'] = [10, 6]\n",
    "plt.rcParams['font.size'] = 18\n",
    "plt.rcParams['pdf.fonttype'] = 42\n",
    "plt.rcParams['ps.fonttype'] = 42\n",
    "# plt.rcParams['mathtext.fontset'] = 'cm'\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4268d987",
   "metadata": {},
   "outputs": [],
   "source": [
    "def loc_1D_gradient(image_properties, b_array, b_fixed=0):\n",
    "    '''\n",
    "    output dimensions: (len(z) x len(x), 5)\n",
    "    '''\n",
    "    b_3d = np.zeros((len(b_array), 3))\n",
    "    b_3d[:, image_properties['fixed_axis']] = b_fixed * np.ones(len(b_array))\n",
    "    b_3d[:, image_properties['main_axis']] = b_array\n",
    "    \n",
    "    if image_properties['gradient_along'] == 'z' and len(b_array) == len(image_properties['z']):\n",
    "        loc_x = np.zeros(len(image_properties['z']))\n",
    "        loc = np.array([loc_x, image_properties['z']]).transpose()\n",
    "        b = b_3d\n",
    "        \n",
    "    elif image_properties['gradient_along'] == 'x' and len(b_array) == len(image_properties['x']):\n",
    "        loc_z = np.zeros(len(image_properties['x']))\n",
    "        loc = np.array([image_properties['x'], loc_z]).transpose()\n",
    "    else:\n",
    "        print(\"ERROR: dimension mismatched\")\n",
    "    return np.hstack((loc, b_3d))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e436675",
   "metadata": {},
   "source": [
    "# Generate config files: velocity samples and b-fields\n",
    "## 1. Basic $b_z$ gradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f7d15b1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "bz = np.linspace(0, 75, 151)\n",
    "bx = by = np.zeros(len(bz))\n",
    "b = np.array([bx, by, bz]).transpose()\n",
    "np.save(\"./configs/bz_linear_75.npy\", b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8ed879b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "bz = np.linspace(-25, 25, 51)\n",
    "np.save(\"./configs/bz_2Dmap_25.npy\", bz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "11662ce3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-25., -24., -23., -22., -21., -20., -19., -18., -17., -16., -15.,\n",
       "       -14., -13., -12., -11., -10.,  -9.,  -8.,  -7.,  -6.,  -5.,  -4.,\n",
       "        -3.,  -2.,  -1.,   0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,\n",
       "         8.,   9.,  10.,  11.,  12.,  13.,  14.,  15.,  16.,  17.,  18.,\n",
       "        19.,  20.,  21.,  22.,  23.,  24.,  25.])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22d79270",
   "metadata": {},
   "source": [
    "## 2. $b_z$ gradient with offset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "a658e558",
   "metadata": {},
   "outputs": [],
   "source": [
    "bz = np.linspace(0, 75, 151)\n",
    "by = np.zeros(len(bz))\n",
    "bx = 10 * np.ones(len(bz))\n",
    "b = np.array([bx, by, bz]).transpose()\n",
    "np.save(\"./configs/bz_linear_75_bx_10.npy\", b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec85aca0",
   "metadata": {},
   "source": [
    "## 3. Load Hongquan's data and divide it by 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "bb599864",
   "metadata": {},
   "outputs": [],
   "source": [
    "b_old = np.load(\"./configs/bz_frames_2.npy\")\n",
    "b_new = b_old / 2\n",
    "np.save(\"./configs/bz_frames_2_new.npy\", b_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94328d6d",
   "metadata": {},
   "source": [
    "## 4. $v_x$ sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "8b27b115",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoEAAAGGCAYAAAAEvludAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABaJElEQVR4nO3dd3gVVf7H8fc3vQEpdFBBRUFcAQEVXYoFAbFgRVy7gBXFXte2rmBBURQXYRV0URRFUUCKIohYQUFBQQVRBKQkoaQRkpzfH3PjL8a0C0kmyf28nuc+SWbOmfncJyT5MmfmHHPOISIiIiKhJczvACIiIiJS/VQEioiIiIQgFYEiIiIiIUhFoIiIiEgIUhEoIiIiEoJUBIqIiIiEoAi/A9Q2DRs2dK1atfI7hoiIiEi5li5dus0516ikfSoCg9SqVSuWLFnidwwRERGRcpnZL6Xt03CwiIiISAhSESgiIiISglQEioiIiIQgFYEiIiIiIUhFoIiIiEgIUhEoIiIiEoJUBIqIiIiEIBWBIiIiIiFIRaCIiIhICFIRKCIiIhKCVASKiIiIhCAVgSIiIiIhSEWgiIiISAiK8DuAiNRNre6YGXSfdSP7V0ESEREpia4EioiIiIQgFYEiIiIiIUhFoIiIiEgIUhEoIiIiEoJUBIqIiIiEIBWBIiIiIiFIRaCIiIhICFIRKCIiIhKCVASKiIiIhCAVgSIiIiIhSEWgiIiISAhSESgiIiISglQEioiIiIQgFYEiIiIiIUhFoIiIiEgIUhEoIiIiEoJUBIqIiIiEIBWBIiIiIiFIRaCIiIhICIrwO4CISF3U6o6ZQbVfN7J/FSURESmZr1cCzSzMzG40s1VmlmNm681slJnFV7D/eWb2opktN7M9ZubMrFUZ7RuY2Rgz2xA430ozu9rMrNLelIiIiEgt4Pdw8JPAE8B3wDBgKnA98K6ZVSTbNcD5QDawpqyGZhYFzAOuAl4LnG81MBa4by/zi4iIiNRKvg0Hm1l7vEJsmnPu7CLbfwaexivuXinnMBcDG51zeWb2DHBoGW0HA12B651zYwLbxpvZm8BdZvaic+6XvXw7IiIiIrWKn1cCBwEGjC62fTyQBVxY3gGcc7865/IqeL4LAscdX2z7aCASGFjB44iIiIjUen4WgV2BAuCLohudcznAssD+ShEYWj4S+Dpw/KK+COSotPOJiIiI1HR+FoHNgW3Oud0l7NsANAzcx1cZkoDYwHH/JHD+VKBFaZ3NbKiZLTGzJVu3bq2kSCIiIiL+8bMIjANKKgABcoq0qaxzUc75Sj2Xc+5551wX51yXRo0aVVIkEREREf/4WQRmAdGl7Isp0qayzkU556usc4mIiIjUeH4WgRvxhnxLKsxa4A0V51bSudLxppH5y5Bv4PwplDBULCIiIlJX+VkEfhk4/1FFN5pZDNARWFJZJ3LOFQBfAZ1KKDqPCuSotPOJiIiI1HR+FoGvAQ4YXmz7ELz78yYXbjCzZmbW1sz25R7BVwPHHVps+3AgD3h9H44tIiIiUqv4Nlm0c+5bM3sWuM7MpgGzgHZ4K4Ys5M8TRY8ALgGOBxYUbjSzHkCPwJddAh+vM7PtgXM8VOQY44HLgCcCS8t9D5wCnAk85Jz7uRLfnoiIiEiN5lsRGDAcWId3da4/sA0YA9wbGMItzwn8dcm3m4t8/kcR6JzLNbOTAtsG4d0HuAZv1ZJn9y6+iIiISO3kaxHonMsHRgVeZbW7FLi0hO33A/cHcb7twHWBl4iIiEjI8vOeQBERERHxiYpAERERkRCkIlBEREQkBKkIFBEREQlBfj8dLCIiQKs7ZgbdZ93I/lWQRERCha4EioiIiIQgFYEiIiIiIUhFoIiIiEgIUhEoIiIiEoJUBIqIiIiEIBWBIiIiIiFIRaCIiIhICFIRKCIiIhKCVASKiIiIhCAVgSIiIiIhSEWgiIiISAhSESgiIiISglQEioiIiIQgFYEiIiIiIUhFoIiIiEgIivA7gIhITZObm8snn3zC7NmzmTNnDhs3biQxMZHExESSkpJK/Vj4eXJyMq4gHwsL9/utiIiUSkWgiAiwbt06Zs+ezezZs/nggw/IyMj40/4tW7YEdbzweo1I6HAyCUecTES9lMqMKiJSKVQEikhIys7OZuHChX8UfqtXr/7T/sMPP5y+ffvSt29fDjvsMHbs2MH27dtJT08nPT29zM83btzI5s2b2fHxZHYsfpW4NseQ0OkUYg44AjPdhSMiNYOKQBEJGWlpabz00kvMnj2bhQsXkpOT88e+Bg0a0Lt3b/r27UufPn1o2bLlH/ta3TGz2JEaBF4HeL9FGwZeAdGugMa/fEPG17PI+vEzsn74hKwfPiEiqRkJHfqRcMRJhMfWr8J3KiJSPhWBIlLnFRQUMGnSJG699VZSU1P/2N6lS5c/rvYdffTRRERUzq9EszBiW3UktlVH8nalkvHNXDKWzyEvfRPbF7zA9kUvE9/27yR07Ed0i3aYWaWcV0QkGCoCRaROW7FiBVdffTUff/wxAD179mTw4MGcfPLJNG7cuMrPH1EvhcTjBtGg23lkr1nCrmWzyFn7FZkrPyRz5YdENmpFvY79iG9/PGHRcVWeR0SkkIpAEamTMjIyeOCBB3jyySfJz8+ncePGPPHEE1xwwQW+XHmzsHDi2hxNXJuj2bP9dzKWzybjm3ns2bqOtHnPkb5wIkk9LyGhU39dGRSRaqEiUETqFOccb7/9NjfccAPr16/HzLjmmmv497//TWJiot/xAIhMbEpSz0tJ/Ps/yPrhU3Z9PYvd61eQNu8/ZK9bRkq/GwiPred3TBGp41QEikid8fPPPzNs2DBmzvQe5OjcuTPPPfccXbt29TlZySw8kvh2PYhv14PM7xeROucZsn/8jE2//0TD024hZr/D/Y4oInWY5ioQkVovNzeXhx9+mPbt2zNz5kzq16/PmDFj+Pzzz2tsAVhcfLvuNL/saaKaH0r+rm1sfvUutn/8Cq4g3+9oIlJHqQgUkVrtww8/pEOHDtx9991kZ2czaNAgVq1axXXXXUd4eO1asSOiQROaXvAI9budB86xY/ErbJ5yN3k7t/odTUTqIBWBIlIr7dq1i4svvpgTTjiBVatWccghhzBv3jxeeeUVmjVr5ne8vWbhEST1uJjGA/9FeEIyu9evYNOL15P142d+RxOROkZFoIjUOlu3buWEE07g5ZdfJjo6mgcffJBvvvmGk046ye9olSa2VUeaXTaG2AO7UJCzi63THiJt3nO4vFy/o4lIHaEHQ0SkVvnll184+eST+eGHH2jdujW7T7yd/2a25L/3vV+h/utG9q/ihJUnPK4Bjc65j11LppO+YCK7vppJzvqVNDr9diIb7ud3PBGp5XQlUERqjZUrV3Lsscfyww8/cMQRR7B48WIiU1qW37EWMzPqdx1A04seJyKpGXu2rmPTpOHsWj4H55zf8USkFlMRKCK1wieffEL37t3ZuHEj3bt3Z+HChbX63r9gRTc9mGaXPEX84Sfg8naTNnsMF1xwAbt37/Y7mojUUioCRaTGmzVrFieddBLp6emcfvrpzJkzp8ZM/FydwqLjaNj/JlJOvRmLimXKlCmcc8455ObqPkERCZ6KQBGp0V5++WVOP/10srOzufzyy3nzzTeJjY31O5avEtofT9N/PEpycjIzZszgvPPOUyEoIkFTESgiNdaTTz7JxRdfTH5+PnfccQcTJkwgIkLPswFENW7N+++/T1JSEtOnT2fQoEHs2bPH71giUouoCBSRGsc5x5133slNN90EwKhRoxgxYgRm5nOymqVTp07MmzePBg0aMG3aNP7xj3+Ql5fndywRqSVUBIpIjZKXl8fgwYMZOXIk4eHhvPTSS38Ug/JXnTt3Zt68edSvX5+pU6dy0UUXqRAUkQpRESgiNUZ2djbnnHMOL7zwArGxsUyfPp2LLrrI71g1XteuXZkzZw716tVjypQpXHrppeTna81hESmbikARqREKdmfSt29fpk+fTlJSEu+//z79+9eeiZ39dswxxzB79mwSEhKYPHkyl19+uQpBESmTr0WgmYWZ2Y1mtsrMcsxsvZmNMrP4II5xipl9YmaZZpZmZlPNrHUpbfc3s3Fm9pOZZZvZBjN718x6VN67EpFgFezOYvOrd/HRRx/RokULFi1axLHHHut3rFrn2GOPZdasWcTFxfHSSy8xdOhQCgoK/I4lIjWU31cCnwSeAL4DhgFTgeuBd82s3GxmdhYwA4gFbgUeA3oAi82sebG2zYGlwHnAG8B1wH+ADsCHZqZLDiI+cAX5bHv3MXI3r+Gggw5i8eLFtG/f3u9YtVb37t2ZOXMmsbGxvPDCC1x11VUqBEWkRL7NtWBm7fEKv2nOubOLbP8ZeBo4H3iljP6RwBhgPdDdOZcR2P4eXrF3PzC0SJdLgIbAAOfc9CLHeRX4ERgCzKyM9yYiFZc+fwLZa74kLKYe7733HgcccIDfkWq9Xr16MWPGDPr378/48eOJiIjg2Wef1dPVIvInfl4JHAQYMLrY9vFAFnBhOf17As2BCYUFIIBzbhmwABgYKBQL1Q983FjsOL8DBUBmxaOLSGXYufRddi19F8IjaHTW3bRp08bvSHXGCSecwDvvvEN0dDTPPfcc119/vdYaFpE/8bMI7IpXfH1RdKNzLgdYFthfXn+AT0vY9xle0XdIkW1zAh/HmllPM2thZl2BV4EMYFRQ6UVkn2St+ZL0D8YDkNLvBmL2O9znRHVP7969efvtt4mKiuKZZ57hxhtvVCEoIn/wswhsDmxzzpW0+vkGoKGZRZXTv7BtSf0BWhRucM4tAK4FWuNdKfwNrwA9FDjGOfdVaScys6FmtsTMlmzdurWMSCJSEbmb17Jt+iPgCmhw3CAS2h/vd6Q6q2/fvkybNo3IyEieeuop7r77br8jiUgN4WcRGAeUVAAC5BRpU1Z/SjlGaf23AkvwHiI5I/CxATDTzPYr7UTOueedc12cc10aNWpURiQRKU/erm1seeMB3J4c4g/rRYPjLvA7Up3Xv39/3nzzTSIiIhgxYgSvvvqq35FEpAaolCLQzKL3olsWUFq/mCJtyupPKcf4S38zG4L3oMktzrnHnXPvOOceB04C9gNGVDS4iOydgtxstr75L/IzUolueRgp/W7QwwrV5LTTTuOpp54C4IorrmDZsmX+BhIR31W4CDSzfmZ2f7Ft15jZTiDTzF4p9iBGeTbiDfmWVMS1wBsqzi2nf2HbkvrDn4eK7wRWOedWFG3onPsWWIX3oImIVBFvKpjHyd28hojEZjQ6824sIphfGbKvrr76ai677DKys7M588wzSU1N9TuSiPgomCuBtwJtC78ws3bAU3jF2DxgIN49dxX1ZeD8RxXdaGYxQEe8Ydvy+gN0K2HfMcBO4Ici21oA4aUcKwIfp8sRCQXpH75A9k+fExaTQONz7iM8roHfkUKOmTF27FiOOuoo1q1bx8CBA7XOsEgIC6bwaQfMKvL1QCAbOMo5t9PMXsGbi290BY/3GnAXMBxYVGT7ELx7+SYXbjCzZnj37v3qnCsc4l0IbAIGm9mTReYJ7AD0Al50zu0pctzvgCPM7Bjn3GdFjt0N7ylizREoUkV2fT2LXUumQ1gEjc68i8iUliW2a3WHfgyrWkxMDG+++SadO3fmgw8+4M477+Sxxx7zO5aI+CCYK4FJwLYiX58EzHfO7Qx8vQDvydsKCQzDPgucZWbTzGywmY3CW0FkIX+eKHoE8D1FrhoGCrwb8O7nWxQYmr4DmIv3AMh9xU5Z+PU8M3ss8MTvY3hXMXOBByqaXUQqLnvtUtLm/QeAlL7DiNn/CJ8TScuWLXnjjTeIiIjg8ccf14MiIiEqmCJwG3AAgJnVw5un7+Mi+yMpfbi1NMOBW4D2eAXh+XirgJzqnCt3nSPn3FTgdLwnhB8Hbse7qnicc25DsbbvAL0DmS8HxgKX4s0f2M05tzTI7CJSjtyt69g6fSS4Aup3G0jC3070O5IEdO/endGjRwPegyLLly/3N5CIVLtghoM/Ba4ys5VAv0DfosPDB+MNz1aYcy4fb5LmMidqds5dilewlbRvBt76wRU533xgfjAZRWTv5GWksWXqA7jcbOLadiex+z/8jiTFXHPNNSxZsoSJEycyYMAAlixZQkpKit+xRKSaBHMl8L5A+9eBy4CXnHPfAZg3x8OZwOJKTygitU5mZqY3FcyurUQ3b0vD/jdi5ue0pFISM+O5556ja9eurFu3jvPPP18PioiEkAr/Vg4UfO3wJlnu5Zy7rMjuROBJKv5QiIjUUc45hgwZQu7vPxLRoAmNzv4nFlHW4j/ip8IHRRo3bsz777/PXXfd5XckEakmQf3X3DmX5px71zn3UbHt6c65p5xzuqlEJMRNnDiRV199FYuMoZGmgqkV9ttvP6ZOnUpERASPPfYYU6ZM8TuSiFSDvRqfMbM4M9vPzPYv/qrsgCJSe6xatYrrrrsOgOTeVxPVUL8SaosePXrw5JNPAnD55ZfrQRGREBDMiiHhZnaXmW0AdgHrgJ9LeIlICMrJyeH8888nKyuLf/zjH8QffoLfkSRI1157LZdccolWFBEJEcE8HfwEMAz4CpgKpFdJIhGplW677TaWL1/OQQcdxHPPPcff/v1R+Z2kRjEz/vOf/7By5UqWLFnCoEGDWN3xOiwsuNm/1o3sX0UJRaQyBVME/gOY5pw7p6rCiEjt9M477zBmzBgiIiKYMmUK9erV8zuS7KWYmBimTZtG586dmTdvHvV31CPp+Mv9jiUiVSCYewIj8VbjEBH5w2+//cZll3mTBYwYMYIuXbr4nEj2VdEHRXZ+MY3M1Zr9S6QuCqYI/AQ4rKqCiEjtk5+fz4UXXkhaWhp9+vThpptu8juSVJKePXvy+OOPA5A2ewx5u7aV00NEaptgisDbgAvM7IyqCiMitcvDDz/MwoULadKkCZMmTSIsTBNC1yXXX389MQd2piAng9RZT1GB1TxFpBYJZrLob4EhwJtm9quZLTSz+cVeH1RdVBGpST7++GPuv/9+AF5++WWaNGnibyCpdGZGw37DCYutT866r9m1tEIrdIpILRHMFDGn4C0ZFwbUB/YHWhd7HVgFGUWkhklLS+OCCy6goKCA22+/nd69e/sdSapIeEISKX2HAZC+4EVyt/7icyIRqSzBjN2MBNYDHZxzic651iW9qiiniNQQzjkGDx7M+vXrOfroo/nXv/7ldySpYnGHdCP+b70hfw/bZozC5e/xO5KIVIJgisA2wNOBYWERCVH/+c9/eOutt6hfvz6vvvoqkZGRfkeSapB84hAiEpuyZ8tati+a7HccEakEwRSBvwAxVRVERGq+b7/9lhtvvBGA559/ntatdfE/VIRFx5HS/2awMHZ+/iY561f4HUlE9lEwReDTwGAzS6iqMCJSc2VlZTFw4EB2797NFVdcwcCBA/2OJNUspmU76h9zLuDYNuMJCnZn+h1JRPZBMCuGZADbge/N7EW8dYLzizdyzr1UOdFEpCYZPnw433//PW3btuWpp57yO474JPG4QeT8vJTc338i7f1xNOyvuSFFaqtgisCJRT6/p5Q2DlARKFLHTJ06lfHjxxMdHc2UKVOIj4/3O5L4xMIjaHjqzWyaOJzMFfOJPego4tv+3e9YIrIXgikCj6+yFCJSY/3yyy8MGTIEgMcff5wOHTr4nEj8FpmyH0nHX07avOdIm/MM0S3aElGvod+xRCRIFS4CnXMLqzKIiNQ8zjmGDBnCjh07OP3007n22mv9jiQ1REKnU8ha8wU5a5eSOuspGp/3AGZaMUakNtnrn1gza2hm+q+fSB323//+l3nz5pGSksLzzz+PmfkdSWoIMyOl3w1aTUSkFguqCDSz5mY2ycy2A5uBzWaWbmYTzaxFlSQUEV/8+uuv3HSTd9P/mDFjtCyc/EVEQjIpfa8DYPvCieRu+9XnRCISjGCWjdsfWAJcBKwFXgm81gIXA1+Y2X5VEVJEqpdzjqFDh7Jr1y4GDBjA+eef73ckqaHiDjmW+L+dhMvLZdu7j2s1EZFaJJgrgf8CkoBTnXNHOucuCrw6A/2B5EAbEanlJk6cyJw5c0hKSuK5557TMLCUKfnEoUQ0aOKtJvKxVhMRqS2CKQJPBsY652YV3+Gcew94DuhbWcFExB+//fbbH6uCPP300zRt2tTnRFLThUXHkXJqYDWRz95k0aJFfkcSkQoIZoqYJODHMvb/CCTuUxoRqRat7phZ4nbnHFvfeIDsHTuIPfgo7v42kXsCbdeN7F+dEaUCSvs+lqWqvo8xLQ+j/jHnsvPT1zj+1HNofvkzhEXHVcm5CunfpMi+CeZK4G9ArzL29wi0EZFaKnPlfLLXLiEsOp7kk6/VMLAEJfG4QUQ1PZj8nVtIXzjJ7zgiUo5gisCpwLlmNsLMGhRuNLP6ZvYwcB7wWmUHFJHqkbcrlfT3nwcg6cShRNRL8TmR1DYWHkFKvxsgLJyMr2eSs36F35FEpAzBPhjyKXA7sM3MfjGzX4BU4A7gE+Chyo8oIlXNOUfa3Gcp2J1J7IFdiD/8BL8jSS0V1bg1DY45F4DU2WNwebk+JxKR0lS4CHTOZQE9gSuBuUAmkAXMAYYCxzvnsqsipIhUrczvFpD90xdYVBzJfa7TMLDskwbdBhKZsh95aRvYvvhVv+OISCmCmizaOZfvnBvvnOvvnDvMOdfOOXeqc26Ccy6vqkKKSNXJz0gn/f1xACSfOJiI+loISPaNRUSS3Pd6wNj5+Zvkbl7rdyQRKUEwk0XPN7MTy9h/vJnNr5xYIlIdnHOkzn2WgpwMYlofSfzfevsdSeqImJbtqNf5VHAFpL73FK4g3+9IIlJMMFcCewFlrRvVGG+4WERqiazvPyL7x8+wqFhS+g7TMLBUqsQeFxNevzG5m9ew84u3/I4jIsUENRxcjkRgdyUeT0SqUH5mOmmBYeCk468gon4jnxNJXRMWFUtKn2sB2P7xZPakbfA5kYgUVeZk0WZ2BNCxyKbuZlZSn2TgGuC7yosmIlUpbd5/KMjeScwBHUno0MfvOFJN9maC6X0Re2Bn4g8/kcwVH5A6ewxNBj2MWWVefxCRvVXeiiFnAvcFPnd4TwZfWUrbXcD1lZRLRKpQ5qqPyVq92BsG7ne9hoGlSiWdMJjstUvZvX4FGctmU6/TKX5HEhHKLwInAgsAA+YDDwPzirVxQAbwnXMup5LziUgl27p1K2lzxwKQ1OsyIho09jmR1HXhsfVI7n0V26aPJH3Bi8QedJSeQhepAcosAp1zvwC/AJjZZcBC59y6asglIlVk2LBhgWHgI0jo2NfvOBIi4g49jtg2x5D942ekzX2WRmffqyvQIj4LZrLoSSoARWq3d999l9deew2LjCG57/W6N0uqjZmR3PtqLDqe7DVfkvX9R35HEgl55Q0H/4WZdQGOBpL4axHpnHP/qoxgIlK5MjIyuPZa70nNxO4XEZnY1OdEEmoi6qWQdPzlpM0eQ9r744hp1ZHwuAbldxSRKlHhItDMYoFpwMl49wi6wEeKfO7w1hgWkRrmn//8J+vXr6dz585s7Xyq33EkRCUccTKZ3y1k96/fkP7BeBqedovfkURCVjBjQffiFYD/Bo7HK/ouAfoBi4AvgcMqO6CI7LslS5bw9NNPEx4ezvjx47GwcL8jSYgyM29i8ohoMr9bQNaaL/2OJBKygikCzwGmOufuBVYEtm1wzs0BTgKigEsrN56I7Ku8vDyGDBlCQUEBw4cPp1OnTn5HkhAXmdSMxO4XApA251kKdmf5nEgkNAVTBO4HLAx8XrgIZBSAcy4PeBU4P5iTm1mYmd1oZqvMLMfM1pvZKDOLD+IYp5jZJ2aWaWZpZjbVzFqX0f4wM3vFzDaZ2W4z+83M3jKzspbEE6m1Ro8ezbJlyzjggAN44IEH/I4jAkC9LqcT1awN+bu2kb5wot9xREJSMA+G7CrSfhdQADQvsn8HEOyd5k/iTTD9FjAKaBf4upOZneScKyirs5mdBbwBLAduBRoAw4HFZtbFObexWPs+wNvAGuBpYDPemsfdgPqBr0V8FeyKDutG9i9937p13HefN9/7c889R3x8hf9/JVKlLCyclH43sGniDWR8PYv4dj2I2e9wv2OJhJRgisA1wCEAzrl8M1uJN0T8gnmTPZ0FrK/owcysPTAMmOacO7vI9p/xCrTzgVfK6B8JjAmcs7tzLiOw/T1gKXA/MLRI+8aB4y0ATnfO7aloVpHayDnH1VdfTVZWFueffz79+vXzO5LIn0Q1akWDY85jxyevkjp7DM0vewaLiPQ7lkjICGY4+H3gbDMrvKN8HNDXzNYAP+LdF/jfII43CO/hktHFto8HsoALy+nfE+9K5ITCAhDAObcMr9AbGCgUC12Ft8bxbc65PWYWV2y/SJ3y2muvMXv2bBITExk9erTfcURK1KDbeUQktyQvbQM7PpvqdxyRkBJMETiS/38qGOfcWOAWvGHgdOAu4NEgjtcVb0j5i6IbA0vPLQvsL68/wKcl7PsMb3j3kCLbTgF2AolmtgzIBHLMbJGZlXcukVolPT2dG264AYBHH32UJk10y6vUTBYRSUrf6wDY8dnr7Emt8ICSiOyjYFYMyXDOrQ48BFK47Qnn3JHOua7OuUeccy6IczcHtjnndpewbwPQ0Myiyulf2Lak/gAtimw7FG/4ezZekXkOcBtwOLAgMDxdIjMbamZLzGzJ1q1by4gkUjPcdtttbNmyhe7du3PFFVf4HUekTDH7HU7CESdDfh6pc54luD8lIrK3/FwzKg4oqQAEyCnSpqz+lHKMkvrXC3z9lnPuUufcm865UcCZge33lnYi59zzzrkuzrkujRo1KiOSiP8++ugjJkyYQGRkJOPGjSMsTEvDSc2X2OsywuIasHv9CjJXfOB3HJGQ4OdfhywgupR9MUXalNWfUo5RUv/swMeJRRs65xYAvwK9yjiXSK2we/duhg71noe66667aNeunc+JRComPLYeSScMBiD9wxfIz9rhcyKRus/PInAj3pBvSUVcC7yh4txy+he2Lak//Hmo+LfAx99LaL8Jby1kkVpt5MiRrF69mkMPPZQ777zT7zgiQYk/rBcxB3SgIHsn6R++4HcckTrPzyLwy8D5jyq60cxigI7Akgr0B2+Ov+KOwXsI5Ici2wofQGlZQvuWwJZyzidSo61atYqHH34YgHHjxhEdXdqFdpGaycxIPvkaCI8kc8UH5Pzyjd+RROo0P4vA1wCHN7lzUUPw7tGbXLjBzJqZWVszK3qP30K8K3iDzSyhSNsOeEO7U4vNBfhy4ONVRU9mZqfhXTmctS9vRsRPBQUFDB06lNzcXK644gp69uzpdySRvRKZ3ILEY73Fp1LnPovLK2tASET2RamTRZvZC8A459znga97AN875yrl8Vjn3Ldm9ixwnZlNwyvCClcMWcifJ4oeAVyCN0XNgkD/PWZ2A14xucjMxuNNC3MjsBW4r9j53jezV4FBZjYLmAEcgDdh9Sa8yaVFaqUXXniBRYsW0bhxYx59NJiZmkJPsCuySPWrf/RZZH63gD2p69nx6VQSu//D70gidVJZVwIvBQ4q8vWHQO9KPv9wvLkG2wPP4q0SMgY4tbwl4wCcc1OB0/GeEH4cuB1YBBznnCtp6piLgTvw3tdo4DK8ZeeOLr7EnEhtsXnzZm699VbAWyc4OTnZ50Qi+8bCI0nucy0AOz6bqrkDRapIWUXgNqDoDLNW2Sd3zuU750Y55w51zkU751o4524qugJIoN2lzjkLPMlb/BgznHPHOOfinHNJzrlznHNrSjlfXmA+w0Odc1HOucbOuYudc/oNI7XW8OHD2b59O3379uX888/3O45Ipfhj7sACzR0oUlXKWjv4E+AeM9sfb0UQgLPM7OAy+jjn3L8qLZ2IlCl77VKmTJ1CbGwsY8eOxVvGW6RuSOx1GVk/fe7NHfjt+yQcUdmDUSKhrawicDgwCe8ePcN7iOOswKs0DlARKFINCvbkkDp3LAAPPPAArVu39jmRSOUqnDswdcYo0j98gdiDjyI8roHfsUTqjFKHg51z65xzPYFYoBVeITgcaF3G68CqjSsihXYsnkL+js106NCB4cOH+x1HpErEH9aLmFadKMjZRfqH//U7jkidUu4UMc65XOfcr3hXBT93zv1S1qvqI4tI7tZ17PzyLcAYN24ckZGRfkcSqRKFcwdaRBSZK+aT/ctyvyOJ1BkVnifQOXdZ4XQxIuIf5wpInT0GCvKpd+QpHH300X5HEqlSkUnNaBCYOzBtjuYOFKksQU0WbWbxZvaAmX1jZhmB1zdmdr+ZxVdVSBH5fxnLZpO7cTXhCckk9rjY7zgi1aL+UWcSmbI/eekb2fHp637HEakTKlwEmlky3tJr/wSaAl8HXk2Ae4EvAm1EpIrkZ6STvnASAEknDiUsWv/3ktBg4ZEk9y2cO/AN9mzTzF4i+yqYK4EPAm2B64BmzrnuzrnuQHPgWuBQtOqGSJVKmz8etzuT2AO7EHfocX7HEalWMS3bk9ChT2DuwGcoKCh3TQERKUMwReDpwATn3FjnXH7hxsCEz88BLwADKjmfiARkr11K1vcfYRHRJJ98teYElJCU2OsywuIS2f3bSl588UW/44jUasEUgU3whn9L8xV/XmFERCpJwZ4c0gJzAjb4+wVENNCPmoSm8JgEkk8cAsCtt97Kli1bfE4kUnsFUwRuBjqVsb9ToI2IVLIdn0whb8dmIhu1on6XM/yOI+KruHY9iGnVifT0dG6++Wa/44jUWsEUge8CV5jZlWb2Rz8zCzOzocDlwDuVHVAk1OVuXcfOL7w5AVP6XIeFl7XQj0jdZ2Yk97mWmJgY/ve//zFv3jy/I4nUSsEUgfcCa4GxwEYzW2hmC4GNwHOBffdVfkSR0OVcAWlznoWCfBI6nUJ0i7Z+RxKpESITm3Lffd6fnKuvvprs7GyfE4nUPsFMFp0KdAFGAqlA18BrGzAC6BpoIyKVJGP5XHZv+J7w+CSSempOQJGibr75Zg4//HDWrFnDQw895HcckVonqMminXM7nXN3O+faO+fiAq/DnXP3OOd2VlVIkVCUn5HO9gXe04+aE1DkryIjI3n++ecxMx599FFWrFjhdySRWiWoIlBEqk/a/AkU7M4k5sDOxLX9u99xRGqkbt26cdVVV5GXl8eVV16puQNFgqA7zEVqoOyfvyLr+4XenIC9y54TsNUdM6slU3WdR6SiCv9NFkT1JDx+Cp988gmNTrmeeh37ldh+3cj+e32OYOzNeUT8oCuBIjVMwZ7dReYEHERkYlOfE4nUbGExCSSddCUA6QsmkpeR5nMikdpBRaBIDbPj09fI2/57YE7AAX7HEakV4g49jtiDuuJ2Z5L+wXi/44jUCioCRWqQlStXsvPzNwFI6XOt5gQUqSAzI7n3VVhkNFmrFpG95ku/I4nUeCoCRWqIgoICrrzySm9OwI79iG7Rzu9IIrVKRIMmJP79HwCkzn2OgtwcnxOJ1GwqAkVqiAkTJrB48WLC4hNJ6nmJ33FEaqV6Xc4gsvGB5O/cwo7Fr/gdR6RGq5Qi0MxizayHmfWojOOJhJqNGzdy2223AZB84lDCYhJ8TiRSO1lYOCl9rwOMnV++Te7mtX5HEqmxKutK4P7AAmB+JR1PJKQMGzaMHTt2cOqppxLXtrvfcURqtehmh1Cv86ngCkidMwZXkO93JJEaqbKKwB3AS8DLlXQ8kZDx1ltvMW3aNBISEhg7dmyZcwKKSMUkdr+I8IQUcjf9yK6vZ/kdR6RGqpQi0Dn3u3PuUufcZZVxPJFQsX37dq699loARowYwX777edzIpG6ISw6juTeVwGw/aOXyNu1zedEIjWPHgwR8dEdd9zBpk2b6NatG1dffbXfcUTqlLhDuhHb5hhcbjZp74/zO45IjVPhItDMOpnZtWXsv9bMOlZKKpEQsGjRIsaNG0dkZCTjx48nPDzc70gidU7ySVdiUbFk//Ap06dP9zuOSI0SzJXA+4CyFkTsB9y7b3FEQkNOTg5DhgwB4M4776R9+/Y+JxKpmyLqNyKx+4UAXHfddezatcvnRCI1RzBFYFdgYRn7FwJH7VsckdDw73//m9WrV9O2bVvuuusuv+OI1Gn1jjyVqKZt+O2337jzzjv9jiNSYwRTBDYEylqVe3ugjYiU4dtvv2XkyJEAjB8/nujoaJ8TidRtFhZOSr9hRERE8Oyzz7Jo0SK/I4nUCMEUgVuAssasDqfsIlEk5OXn5zNkyBDy8vK4+uqr+fvf/+53JJGQENX4QO644w4ABg8eTHZ2ts+JRPwXTBH4PjDYzP5SCJrZYcAVgTYiUoqxY8fy+eef07x5c0aMGOF3HJGQcs8999CuXTt++OEHHnzwQb/jiPgumCLwISAf+NLMnjGzwWZ2hZk9AywB8oB/VUVIkbrg119//eN+pLFjx9KgQQOfE4mElujoaP773/9iZjz22GN89dVXfkcS8VWFi0Dn3BrgRGAtcA3wPDA+8Pka4CTn3I9VEVKktnPOcfXVV5OZmck555zDGWec4XckkZDUrVs3rr/+evLz87niiivYs2eP35FEfBPUZNHOuSXOucOBI4GBwPlAJ+fc35xzS6oioEhd8NprrzFr1iwaNGjA008/7XcckZD20EMP0apVK5YtW8Zjjz3mdxwR30TsTSfn3DJgWaUmEanhWt0xc6/65WfvZOMEbzWQ8G4X0+0pDUGJ+CkhIYHx48fTu3dvHnjgAc466yzatm3rdyyRaqdl40SqWPr8FyjI2kH0foeTcMTJfscREeCkk07i8ssvJzc3lyuuuIL8/Hy/I4lUu6CKQDM7zsxmmNlWM8szs/xir7yqCipSG2WvW0bmivchPJKUvsMwM78jiUjAqFGjaNasGZ988gljx471O45ItQtm7eAewIfA0cDngb4fAl8CBqwAXq6CjCK1UsGeHNLmPANA4nGDiExu4XMiESkqMTHxj+LvzjvvZN26df4GEqlmwVwJvBvYBBwGXBrY9rBz7higL9AamFCp6URqsR0fv0Le9t+JbNSK+ked5XccESnBgAEDOPfcc8nMzGTo0KE45/yOJFJtgikCjwImOOe2AgVF+zvn5uJdBdQ8gSLA7o2r2fnl22Bh3jBw+F49gyUi1WDMmDEkJyczb948Jk2a5HcckWoTTBEYDWwIfL478LFekf3LgM7BnNzMwszsRjNbZWY5ZrbezEaZWXwQxzjFzD4xs0wzSzOzqWbWugL9jjCzPWbmzOycYHKLlKVgz262zXwSXAH1u5xBdPND/Y4kImVo0qQJo0ePBuDGG2/k999/9zeQSDUJpgjcBLQEcM5lAtvx1gsu1BJv1ZBgPAk8AXwHDAOmAtcD75pZudnM7CxgBhAL3Ao8BvQAFptZ8zL6heFNdJ0TZF6Rcu34eDJ5ab8RkdySBt0v9DuOiFTAhRdeSN++fdm+fTvXXnut33FEqkUwReCXwHFFvp4L3GhmF5vZpcB1eA+MVEhgDeJhwDTn3FnOufHOuZuAm4Dj8SaiLqt/JDAGWA90d86Ndc6NAPoATYD7y+g+DGiPVzSKVJqc375j5xdvgYXRsP+NhEVG+x1JRCrAzBg3bhwJCQlMmzaNN9980+9IIlUumCLwv8A2M4sNfH0XkA1MBF7AGyK+LYjjDcJ7qnh0se3jgSygvEsoPYHmePcpZhRuDExkvQAYGCgU/8TM9sNbB/l+4Ncg8oqUqSA3h9RZTwKO+kefrWFgkVpm//3355FHHgHg2muvJS0tzedEIlUrmLWD5znn/uGcyw58vRY4BBgAnAa0c86tCOLcXfEeMPmi2Hly8O4v7FqB/gCflrDvM6B+IF9xY/HWPx5d8agi5dv+0STy0jcR2agVicdd4HccEdkLV111Fd27d2fz5s3cdNNNfscRqVIVKgLNLDYw7Ht00e3OuUzn3DvOuZnOuR1Bnrs5sM05t7uEfRuAhmYWVU7/wrYl9Qf408RsZjYQ6A9c5ZzTxNZSaXJ+/YZdS9+FsHBSThmORfzlIrSI1AJhYWFMmDCB6OhoJk2aRPbapX5HEqkyFb0SuBtvDsBOlXjuOP7/KePicoq0Kas/pRzjL/3NLBHv6t9451xJVw9LZWZDzWyJmS3ZunVrMF0lBBTkZrNt1lMANOh2HtFND/Y5kYjsi0MOOYQHHngAgNTZz1CwO9PnRCJVo0JFoHOuAO/+ufqVeO4svGlnShJTpE1Z/SnlGCX1fxzv/d5R0YCFnHPPO+e6OOe6NGrUKNjuUself/gC+Ts2E9XkIBp0G+h3HBGpBDfffDNdu3Ylf9dW0t4f53cckSoRzIMhk4CLzKyyHnfciDfkW9LxWuANFeeW07+wbUn9ITAsbGZHApcDzwApZnawmR0MNA60axrYpkc5JSjZP39NxrL3ICyClP43alJokToiIiKCl19+GYuIJnPFfDJXL/Y7kkilC6YI/ARvHsBlZjbMzPqaWY/iryCO92Xg/EcV3WhmMUBHYEkF+gN0K2HfMcBO4IfA1/vjPYn8IPBjkdcjgf1jAl//LYj8EuIKdmeS+t7TACT+/QKiGrXyN5CIVKpDDz2UxF6XApA251nyMvS0sNQtwVy2mFfk86eA4gssWmBbeAWP9xreNDPDgUVFtg/Bu5dv8h8HNmsGNAB+dc4VDvEuxJvAerCZPVk4TYyZdQB6AS865/YE2n4BnFtChl7AtcAovCeK11QwuwhpH0wgf9dWopq1of7RZ/sdR0SqQL0j+5P90xfkrPuatPeeptE592FmfscSqRTBFIGXVeaJnXPfmtmzwHVmNg2YBbTDWzFkIfBKkeYjgEvwJpFeEOi/x8xuwCsmF5nZeLx7Fm8EtgL3FTnXRuCN4hnMLCHw6WfOub/sFylN9povyfx2HoRH0vCUm7Cwiv7fR0RqE7MwUk65gU0vXEf22iVkLJ9NvY79/I4lUinKLALN7CjgJ+dcmnOuKlbVHg6sA4biTd2yDW9o9t7Awyhlcs5NNbNs4B68Bz92Ax8AtzvnSpo6RmSf5edkkDp7DACJ3S8isuF+PicSkaoUUa8hySdfw7Z3HiV9/gRiDuhAZFKpK5OK1Brl3RP4KdC38AszSzCzV8zssMo4uXMu3zk3yjl3qHMu2jnXwjl3U9EVQALtLnXOmXNuQQnHmOGcO8Y5F+ecS3LOneOcq9CwrnNuYuC4ugooFZb+/jjyM9KIbtGO+l3P8DuOiFSD+HY9iGvXE7dnN9tmjMIV5PsdSWSflVcEFr/xIRpvTd+mVRNHpGbL+uFTMld+iEVEe5NCaxhYJGQkn3w14Qkp5G5czc7PdO1Aar9gng4WCWn5WTtInfMsAIk9LyEyuaTZiUSkrgqPSSDllOEAbF/8Crt//8nfQCL7SEWgSAWlzfsPBVnbid7vcOp1PtXvOCLig9jWnajX+TQoyCd1xigK9pS28JVIzaciUKQCXn/9dbJWLcIiY7xhYNOPjkioSux5CRHJLdmTup7tH73kdxyRvVaRKWJOMbPCewDj8OYCPNfMOpbQ1jnnnqyscCI1wa+//sqVV14JQNLxlxOZqFtiRUJZWGQMDU+9md//dwu7lkwn9qCuxLbq6HcskaBVpAi8IPAq6spS2jpARaDUGXl5eVxwwQVs376d2IO6kqD5wUQEiG7WhgbHns+OjyeTOms0zS9/hrCYhPI7itQg5RWBx1dLCpEa6sEHH2Tx4sU0b96csFOGa6UAEflDg27nkb1mCbmbVpP2/jgannqz35FEglJmEeicW1hdQURqmg8//JCHHnoIM2Py5MlcOjvT70giUoNYWDgNT72JTS9eT+bKD4k9+Gji2/7d71giFaa720VKsG3bNi688EKcc9xzzz306tXL70giUgNFJrcg6YQrAEib8yx5u1J9TiRSccGsHSxSY7W6Y2alHcs5x9Y3HyR740aiWx7GpJwuvFSJxxeRuiWhYz+yfvycnJ+XkvreUzh3kW4dkVpBVwJFitm15B2y13xJWEwCDU+7RauCiEiZzIyUftcTFlOPnJ+/4tlnn/U7kkiFqAgUKWL37z+RvuBFAFL6XU9E/cY+JxKR2iCiXgrJfa4F4Oabb2bJkiU+JxIpn4pAkYCC3Vlse+dRKMgjoVN/4g451u9IIlKLxLf9Owmd+pObm8u5555Lenq635FEyqQiUCQgbd5z5KVvJLJRK5KOv9zvOCJSCyWfMJjOnTuzbt06LrvsMpxzfkcSKZWKQBEgY8V8Mld+iEVG0+j02wmLjPY7kojUQhYRydSpU2nQoAHTp0/niSee8DuSSKlUBErI25O2gbS5YwFIOvFKIhvu53MiEanNWrduzaRJkwC4/fbbWbx4sc+JREqmIlBCmsvbw7Z3HsXtySGubXcSjujtdyQRqQPOOOMMbrnlFvLz8xk4cCBbt271O5LIX6gIlJCWvnAiuZvXENGgCSl9r9PcXiJSaR5++GGOPfZYNmzYwEUXXURBQYHfkUT+RJNFS8jKWvMlu5ZMh7BwGp5+G2HR8X5HEpEqVJmTylf0PHkdhhL29QrmzJlDcq9LSDz2/L+0Xzeyf7XkEilOVwIlJOXt2kbqzCcBSOxxEdHND/U5kYjURRH1G9LwtFsAY8fHr5D9y3K/I4n8QUWghBxXkM+2GaMoyN5JTKtO1D/qLL8jiUgdFtv6SBocOxBcAdveeYy8jDS/I4kAKgIlBG3/+BV2//otYXGJNOx/E2b6MRCRqtXguEHEHHAEBVnbvYfRCvL9jiSiIlBCS+b3i9j56WtgYTQ87RbCE5L8jiQiIcDCwml42q2Exyexe/0Kti/6n9+RRFQESujI3byW1PdGA5B0/BXEturoax4RCS3h8Uk0PP02sDB2fjaV7DVf+h1JQpyKQAkJ+Znb2TLtX7g9u4k//ETqdTnd70giEoJi9v8bid0vBGDbjCfI27nF50QSylQESp3n8vew9e0R5O/cSlSzQ0npc63mAxQR39Q/5hxiD+xCQc4utr79CLm5uX5HkhClIlDqvLT3n2f3bysJT0im0Zl3YRFRfkcSkRBmFkbKqTcRXq8RuZtWc/vtt/sdSUKUikCp03Z9PYuMZe9BeCSNzrybiHopfkcSESE8tj6NzrgdwiIYPXo0//ufHhSR6qciUOqsnPUrSHt/HAApfa/ThNAiUqNEt2hL8omDAbjiiitYtGiRz4kk1KgIlDopb8cWtr49Agryqdd1AAmHn+h3JBGRv6h35KkMGzaM3NxczjzzTH766Se/I0kIUREodU5Bbg5bpj1EQdYOYlp1IqnXZX5HEhEp1RNPPMEpp5xCamoq/fv3Jy1NK4pI9VARKHWKc47U955iz5a1RCQ1o+EZt2Nh4X7HEhEpVUREBFOmTOGII47ghx9+4Oyzz9YTw1ItVARKnbLzs6lkrVqERcXS6Kx/Eh6T4HckEZFy1atXjxkzZtCsWTMWLFjAlVdeiXPO71hSx6kIlDoj66fP2f7Ry4DR8LRbiGq4v9+RREQqbL/99uPdd98lLi6OiRMnMnLkSL8jSR2nIlDqhNxtv7Lt3ccBR2L3C4k7+Gi/I4mIBK1z585MnjwZM+Ouu+5i6tSpfkeSOkxFoNR66enpbJ32L1xuNnFtu1O/23l+RxIR2WsDBgzg0UcfBeDiiy/ms88+8zmR1FUqAqVW27NnD+effz556ZuIbHwgKf1u0JJwIlLr3XzzzQwdOpScnBzOOOMM1q1b53ckqYNUBEqtVVBQwGWXXcbcuXMJi61P47PuISwqxu9YIiL7zMx45pln6N27N1u2bKF///7s2LHD71hSx0T4HUDqvlZ3zAyq/bqR/ctt45zjhhtuYPLkySQkJJBw5n1ENGi8txFFRGqcyMhIXn/9dY477ji+++47zj33XGbOnElkZGTQv1f3RkV+F0vtpiuBUivdf//9PPPMM0RFRTF9+nQtCScidVJiYiIzZsygUaNGzJs3j2HDhmnqGKk0KgKl1hk9ejQPPvggYWFhvPbaa5xwwgl+RxIRqTKtW7fmnXfeITo6mnHjxvHkk0/6HUnqCBWBUqtMmjSJG2+8EYAJEyYwYMAAfwOJiFSDY445hpdeegmAW265hawfPvU5kdQFKgKl1nj77be54oorAG+tzcsu05rAIhI6zjvvPP7973/jnGPbu4+RvW6Z35GklvO1CDSzMDO70cxWmVmOma03s1FmFh/EMU4xs0/MLNPM0sxsqpm1LqFdTzN71sy+NbNdZrbVzBab2SDTnCI13vz58xk4cCD5+fncc889f1wNFBEJJXfeeae3pFxeLlvf/Bc5v37jdySpxfy+Evgk8ATwHTAMmApcD7xrZuVmM7OzgBlALHAr8BjQA1hsZs2LNX8EOANYANwMPAyEA68Az1fCe5Eq8uWXX3LGGWeQm5vLNddcw4MPPuh3JBERX5gZY8eOJf5vvXF5u9nyxoPkrF/hdyyppXybIsbM2uMVftOcc2cX2f4z8DRwPl6BVlr/SGAMsB7o7pzLCGx/D1gK3A8MLdLlduBj51x+kWM8BXwIDDazp5xz+kmqYb777jv69etHRkYGF1xwAWPGjNFk0CIS0sLCwkjpNwxcPpkr5rPljQdofO6DxLRs53c0qWX8vBI4CDBgdLHt44Es4MJy+vcEmgMTCgtAAOfcMryrfQMDhWLh9oVFC8DAtgLgjcCXhwf9DqRKrVu3jpNPPpnU1FT69+/PxIkTCQvz++K1iIj/zMJI6XcDcYf1xOVms2XqvezeuNrvWFLL+PkXtStQAHxRdKNzLgdYFthfXn+Akh6R+gyoDxxSgRwtAx83V6CtVJPNmzfTu3dvNmzYQPfu3Xn99deJjIwsv6OISIiwsHAa9r+JuLbdcbnZbH79XnZv+tHvWFKL+FkENge2Oed2l7BvA9DQzKLK6V/YtqT+AC3KChC4b/BKYC3wcdlxpbps376dPn368NNPP9GxY0feffdd4uLi/I4lIlLjWFg4DU+9mbhDjsXtzmTL6/8kd/Mav2NJLeFnERgHlFQAAuQUaVNWf0o5Rrn9zSwOeAuIBy51zu0po+1QM1tiZku2bt1aRiTZVwV7cjj11FNZvnw5hxxyCHPmzKFBgwZ+xxIRqbEsPIKGp99GbJtjKMjJYPOUe8jd8rPfsaQW8LMIzAKiS9kXU6RNWf0p5Rhl9jezGOBtoAtwmXNuUVlBnXPPO+e6OOe6NGrUqKymsg8KcrPZ+uZDLF68mJYtWzJ37lwaN9Z6wCIi5bHwCBqdfjuxB3WlIGcXm6fcTe7WdX7HkhrOzyJwI96Qb0lFXAu8oeLccvoXti2pP5QwVFykADwJGOKc+1+FE0uVyc/eyeYpd5PzyzIaNWrE3LlzOeCAA/yOJSJSa1hEJI0G3ElM684UZO9k85R72LNtvd+xpAbzswj8MnD+o4puDBRpHYElFegP0K2EfccAO4Efih07Gm8I+GRgqHPuhaBTS6XL27mNzZNvJ3fTD4TXb8zixYtp105THYiIBMsiomh05l3EHNCRgqztbJ5yF3tSf/M7ltRQfhaBrwEOGF5s+xC8e/kmF24ws2Zm1jZwH1+hhcAmvDn+Eoq07QD0AqYWvc8vUAC+DfQBrnLOTajMNyN7Z0/qb/w++Vb2pK4nsuH+NL3wUdq0aeN3LBGRWissMppGZ99D9P5HkJ+Z7hWC6RvL7yghx7ci0Dn3LfAscJaZTTOzwWY2Cm8FkYX8eaLoEcD3FLlqGCjwbgD2AxaZ2TVmdgcwF9gK3FfslJOBvsAHQJaZXVjsdUTVvFMpze7ff+L3ybeRv3MrUc0PpckFjxBRr6HfsUREar2wyBgan30v0fsdTn5GGptfvYs923/3O5bUMH7PvDscuAVoj1cQno+3CsipgYmcy+ScmwqcjveE8ON4q4IsAo5zzhW/H7BL4ONJwMslvM7ax/ciQcj55Rs2v3onBdk7iWndmSYD/014bD2/Y4mI1BlhUYFCsMVh5O/axu8v36IJpeVPfFs2DiCwgseowKusdpcCl5aybwbe+sHlnatV0AGlSmT98Alb33kU8vOIa9eDhv1vxMI1EbSISGULi46j8bn3s/Wth8n5ZRmbX72Lhqd58wqK+FoESula3TEzqPbrRvavoiSVK+ObuaTOfgZcAfWO7E/SSVdi5vcFaRER/wT7+z5YhYVg6pxnyfx2HlvfGkHSCVdQr8sZWos9xOmvr1SbHZ+/Qep7T4MroMFxg0g66SoVgCIi1cDCI0jpdz2JPS4GHOnzJ5D+/n9wBfl+RxMf6S+wVDnnHOkfvsD2BRMBSDrpShL//g/9D1REpBqZGQ26nUfD026F8Ah2fTWTrdMeoiA32+9o4hMVgVKl8vLySH3vaXZ+MQ3Cwml42i3U73ya37FEREJW/GE9aTLwIcJi6pG95ks2v3IHebtS/Y4lPlARKFUmJyeHc889l8xv52ER0TQ+65/EH9bL71giIiEvZr/DaXrR40QkNSN38xp+f/lmLTMXglQESpVYu3Ytxx57LG+//TZh0fE0HvgQsQd1Kb+jiIhUi8jkFjS98HGiW7TzppD5361k//yV37GkGqkIlEo3Y8YMOnfuzNdff82BBx5Ik388QkxLLQMnIlLThMc1oMn5/yaubXdcbjZbpt7PrmWz/Y4l1URFoFSa/Px87r77bk477TS2b9/O6aefztKlS4lq1MrvaCIiUgqLiKLh6bdS/5hzwRWQNucZ0hdMpKCg3DUbpJZTESiVYsuWLZx88sk8/PDDhIWFMXLkSN566y0SExP9jiYiIuUwCyOp5yUk9x0GFsbOz99g0KBBZGZm+h1NqpCKQNlnixcvplOnTsyfP58mTZrwwQcfcPvttxMWpn9eIiK1Sb0OfWh87v1YVCyvv/46nTt35quvdJ9gXaUVQyQoRWe2d86xa8k7pC94AQryiW55GBGn386lszNhdtXOgC8iIlUjtvWRNL3wcba98wirV6+mc9ejSep5MfW6Dih3gv/asnqVeHSpRvZKwe4stk1/hPT546Egn/pdz6TJ+Q8TUS/F72giIrKPohodQNOLn6Tekf2hII/0D19gy+v3kZeR5nc0qUQqAiVouVvXsemlG8la/TEWFUvDAXeSdMIVWLguLIuI1BVhkdEk976aRmf/k7DY+uSs+5pNL1xH1k+f+x1NKomKQAlKxsoP+f3lm8lL20BkwwNodslo4g89zu9YIiJSReIOPppmlz9DTKtOFGTvZOub/yJt3nMU7NntdzTZRyoCpUIyMjK4+uqrSZ0xCrdnN/Htj6fpRaOITG7hdzQREaliEQnJND7vAZKOvxzCvHWHf3/pRnK3/Ox3NNkHKgKlXDNnzqR9+/b85z//gfAIkvtcS0r/mwiLivE7moiIVBOzMOofdZa33FxyS/Zs+5VNL93EziXv4JzzO57sBRWBUqrff/+dgQMHcuqpp/Lrr7/SqVMnml00inod+2FmfscTEREfRDc9mGaXjCahQx/I30P6B8+z9Y0HyM/c7nc0CZKKQPmLgoICnn/+edq1a8frr79OXFwcjz/+OF988QVRTQ7yO56IiPgsLCqGlL7DaDjgTsJiEsheu4SNL17HrFmz/I4mQVARKH/y/fff07NnT6688kq2b99Ov379WLlyJTfffDMREXr6V0RE/l/8ocfR7LIxRO//Nwoyt9O/f38GDBjAmjVr/I4mFaAiUADIycnhvvvuo0OHDnz88cc0btyYKVOmMHPmTFq1auV3PBERqaEi6jeiycCHSOx1OfHx8UyfPp3DDjuMu+66i4yMDL/jSRl0aUdYuHAhV155JatXrwZg8ODBPPLIIyQnJ/uSp+iqJCIiUvNZWDgNjj6LvMN6wsKJZK78kBEjRvDomHEk9rqM+MN6lXovuVYZ8Y+uBIaw9PR0hgwZQq9evVi9ejWHHnooCxcuZPz48b4VgCIiUntF1Euh4ak30/TCx4hq2ob8jDRSZ4xi8+Tb2P37T37Hk2JUBIagvLw8Jk6cSNu2bZkwYQJRUVHcd999LF++nB49evgdT0REarnoFu1oevEoUvrdQFhcIrs3fM/vk24k9b2n9RRxDaLh4BCSn5/PlClTeOCBB/jxxx8B6N69O+PGjaNdu3Y+pxMRkbrELIyEI3oTd+ix7Fg8hZ1L3yHjm7lkrvqYxL9fQL0jT9Vyoz7TlcAQUFBQwJQpUzj88MO58MIL+fHHHznwwAOZNGkSCxYsUAEoIiJVJiw6nqQTrqD55c8S07ozLjeL9PkT2PjCdWSvXep3vJCmIrAOKygo4I033uCII45g0KBBrFq1ilatWvHf//6XVatWcfHFFxMWpn8CIiJS9SJTWtL43PtpdPa9RCQ1Iy/tN7ZMvY8TTzyRuXPnatURH6gCqIOcc7z11lt06tSJc889l5UrV7L//vvz/PPPs3r1ai6//HIiIyP9jikiIiHGzIg7+CiaXz6WxF6XYlGxzJ8/nz59+nDkkUfy6quvkpeX53fMkKEisA5xzvHOO+/QuXNnzjrrLL755htatmzJc889xw8//MCQIUOIioryO6aIiIQ4i4ikwdHn0PLqFxkxYgRNmjRh2bJlXHDBBbRp04ZnnnmGrKwsv2PWeSoC6wDnHDNnzuSoo47ijDPO4Ouvv6Z58+aMGTOGH3/8kauuuoro6Gi/Y4qIiPxJWEwCd9xxB+vWrWPcuHG0adOGdevWMWzYMA444AAefPBBUlNT/Y5ZZ5nG4IPTpUsXt2TJkio/T0UmTM7P2kHmig/YtXwOeWkbAAiPT6L+MeeQ0KEvYZEq/EREpPZwBflk/fgZOz9/g9xN3iwWFhlNwhEnU7/rmUQ0aLxXx62OCan3ZqGD6shlZkudc11K2qdns2sZ5wrI+eUbMpbPIeuHT6HAu3ciPCGZ+l3PJKFTP8IiY3xOKSIiEjwLCyf+0OOIO+RYdq//lh2fvUnOz0vZtfRddn01k/h2Pah/9NlENW7td9Q6QUVgLZGfmU7Gtx+QsXwOeds3eRstjNiDupLQsS+xB3bBwsL9DSkiIlIJzIyY/Y8gZv8jyN2ylp2fTyPz+4/I/G4Bmd8tIKpZG+IPO574dj0Ij0/0O26tpSKwBnOugJx1y8hYNpusnz6HgnwAwus1IqHDyST8rTcR9Rv6nFJERKTqRDU+kIan3UJij4vY+eXbZHz7PrmbfiR304+kz59AbOsjiW9/PLFtjtZIWJBUBNZAGzduZMcnr7Hrm7nk79jsbbQwYtscQ70OfYhpfaSu+omISEiJaNCE5JOuJLHnpWT/9DmZKz8k++evyF67hOy1S7CoWOIOOY749r2I2f9v+jtZASoCa5hZs2Zx+umnk58fuOpXvzH1OvQh/m8nEVEvxed0IiIi/gqLjCa+XQ/i2/XwHpD8/iMyVy4gd9NqMle8T+aK9wlPSCH+sJ7EH34CUY1a+R25xlIRWMMcd9xxJCQksKdJexI69CGmdSfMNJOPiIhIceFxDajf+TTqdz6NPWkbyFz5IZkrPyRvx2Z2fjGNnV9MI7JRK0Y0+Ia+ffvSoUMHrZRVhKaICVJ1TBGTmZlJ+38tqNJziIiI1EXOOXZv+J7MlR+StWoRBTkZf+xr3LgxvXv3pk+fPvTu3ZumTZtW2nk1RYxUivj4eL8jiIiI1EpmRkzLw4hpeRjJJw4l++elnJ60iTlz5vDbb78xefJkJk+eDECHDh04+eST6dOnD8cddxwxMaH1YImKQBEREamTLCKSuDbHMGFkf5xzrFq1ijlz5jB37lwWLFjA8uXLWb58OY899hixsbH06tXrj6Kwbdu2mJnfb6FKqQgUERGROs/MaNeuHe3atWP48OHk5OSwePHiP4rC5cuX89577/Hee+8BEBZbn6imbYhu1oaopm2IataGiIRkn99F5VIRKCIiIiEnJiaGE088kRNPPJFHH32UTZs28f777zNnzhzmzZvHli1byPl5KTk/L/2jT3hCClHN2hDd7BCvMGx6MOGx9Xx8F/tGRaCIiIiEvGbNmnHRRRdx0UUX4Zxjv2smsnvTD+T+/iO7N/1I7u8/kZ+RSvaPqWT/+Nkf/SISm3mFYdODiUzZj4jkFkQ0aFIr5ilUESgiIiJShJkR0aAxEQ0aE9/274C3ilde2gZ2//4TuZt+8FYt2bKWvO2byNu+iazvP/r/A4RFEJnUjIjkFkQGXhHJLYlMbkFYbP0ac6+hr0WgeRPg3QBcCbQCtgKvA/c65zIreIxTgHuADsBu4APgNufczyW0bQA8BJwFpABrgGeA/zjNlSMiIiKlMAsjMmU/IlP2g/bHA+Dy89iz7Vd2//4juZvXkpf2G3vSNpC/axt7UtezJ3U92cWOExaT8Edx+HD95QwdOpSGDf1ZAtbvK4FPAtcDbwGjgHaBrzuZ2UnOuYKyOpvZWcAbwHLgVqABMBxYbGZdnHMbi7SNAuYBnYAxwPdAP2As0AS4vzLfmIiIiNRtFh5BVJMDiWpy4J+2F+TmkJe+kT2BojAvbQN70jawJ+03CnIyyN24mtyNq7n77vlcfvnlPqX3sQg0s/bAMGCac+7sItt/Bp4GzgdeKaN/JF4xtx7o7pzLCGx/D1iKV9QNLdJlMNAVuN45NyawbbyZvQncZWYvOud+qaS3JyIiIiEqLCqmxOLQOUdB5vY/isMhHeNp0qSJTynBz7VTBgEGjC62fTyQBVxYTv+eQHNgQmEBCOCcWwYsAAYGCsVCFwSOO77YcUYDkcDAYMKLiIiIBMPMCE9IImb/v1GvY19Gjhzp6/2BfhaBXYEC4IuiG51zOcCywP7y+gN8WsK+z4D6wCHwx72HRwJfB45f1BeBHOWdT0RERKTO8LMIbA5sc87tLmHfBqBh4D6+svoXti2pP0CLwMckILaktoHzpxZpKyIiIlLn+flgSBze07wlySnSJreM/pRyjJxibcpqW9g+rpR9mNlQ/v/+wgwzW11a2xqkIbDN7xBSYfp+1S76ftUu+n7VLpX+/bJHKvNolaeach1Q2g4/i8AsoHEp+2KKtCmrP0B0BfqX1bawfanncs49DzxfRpYax8yWOOe6+J1DKkbfr9pF36/aRd+v2kXfr+rj53DwRrwh35IKsxZ4Q8WlXQUs7F/YtqT+8P/Dv+lAdkltA+dPoeRhZREREZE6yc8i8MvA+Y8qutHMYoCOwJIK9AfoVsK+Y4CdwA8AgfkGv8Kbf7B40XlUIEd55xMRERGpM/wsAl8DHN7kzkUNwbs/b3LhBjNrZmZtzazofXsLgU3AYDNLKNK2A9ALmOqc21Ok/auB4xadO5DA+fPwViqpS2rV8LXo+1XL6PtVu+j7Vbvo+1VNzM/V0sxsDHAd3oohs/j/FUMWAycUrhhiZhOBS4DjnXMLivQ/F6+YXI43/1994Ea84rKzc25DkbZRwCd4y8s9jbdiyCnAmcBDzrl/VuFbFREREalR/F42bjiwDu/qXH+8p4HG4K0dXOaScQDOualmlo23dvDj/P/awbcXLQADbXPN7CS8tYMH8f9rBw8Dnq2k9yMiIiJSK/h6JVBERERE/OHnPYFSBczsSjObbGarzCzfzMqs8s2suZm9ZGZbzSzbzJYEhtnFR2a2zsxcKa+GfucLNWYWZmY3Bn6ucsxsvZmNMrN4v7NJycr4+ckov7dUFTO708ymmtnawPdjXTntDzWzt80s3cwyzWyRmZ1QTXHrPF0JrGMCP1ApwNdAa6Clc67EhQnNLBnvqejGwBPAb3hrLPcELnfOvVgdmeWvAt/HbODfJeyeWspKO1JFzOwpvPuV3wLew7t/eRiwCDipIrevSPUK/Ad4EX99yGCPc+41HyIJf3xf0vBm7OgM7HTOtSql7UF4S7vmAaOBHXgPjx4O9HPOvV8Nkes0FYF1jJm1An51zhWY2QygfxlF4KPArcDpzrl3A9vC8dZjPgg4wDmn/zX7IFAErnPO9fI5Ssgzs/bAt8Bbzrmzi2wfhveQ2T+cc6/4lU9KFig2JjnnLvU7i/w/MzvQObc28PkKIKGMIvB14Gy8Bz2XBbYlACvxVvpq61TE7BMNB9cxzrl1QVyVuABYU1gABvrn4z2ck4z39LT4yMwizKy+3zlC3CDA8K5EFDUeb6WhC6s7kFScmUUVnUZM/FVYAJYncKvF6cCCwgIw0D8DmAAcAnStioyhREVgiDKzZngrqHxWwu7CbfoB89fReEXGDjPbbmaTzKy536FCUFegAG9Y6g/OuRxgGfo5qcnOwfsZ2mVmW8xsjJk18DuUVMgReEu9flrCPv2NqiR+TxEj/iksJkpaLq9wW0lL8kn1WIn3v91VeD+nvYDBwIlmdpRzbmMZfaVyNcdbxrKk+zA3AMeaWVQ5y1xK9fsCmAr8hDeH7Cl489L2NLNjdatLjae/UdVARWANZGaJ/HUllbI87ZxLC/I0hauvlPSHLadYG9kL+/J9dM71L7Zvipl9hLeSzgN4N0dL9Yij5J8T+PPPiorAGsQ5d3SxTS+Z2Td4D1vdQMkPXUnNob9R1UBFYM2UCNwXRPv/4T1tFYyswMfiaykDxBRrI3snkUr8PjrnXjGzf+NNrC7VJwvvCfqS6GeldnkM72eyPyoCazr9jaoGKgJrIOfcOrwb0atS4XBiSZfTC7eVdBleKqiKvo/rgOMq+ZhSto3AYWYWXcKQcAu8oWJdBawFnHN7zGwjoLk2az79jaoGejAkRDnnNuH9AB1Twu7CbUuqL5FU0MHAZr9DhJgv8X5XHlV0o5nFAB3Rz0mtEfietUQ/Q7XBt3hDwd1K2Ke/UZVERWBoexU4yMxOK9wQmCdwGLAdmOVTrpAWmMS7pO3X4v0Be7ek/VJlXgMcf72/cwjePUmTqzuQlM3MUkrZ9S+8ETD9DNVwgQd33gV6mVmHwu2B6X4GAz9S7Il9CZ4mi65jAgVd4Q/MhcChwD8DX293zj1TpG0KsBRvhZEn8K4MDiLwJKpz7r/VFFuKMLPhwBXAbLzh38KngwcAa4Buzrmt/qQLTWY2Bu/J0rfw/nPUDm8FkcXACVoxpGYxsyfxrhZ9CPwKJOA9HXw88DlwvHMu27+EocvMLgIOCHw5DIgCRgW+/sU593KRtgfjFXp7gCeBnXj/+fob3kIIc6ord12lIrCOMbOJwCWl7P6l+MzsZtYCGAn0w/tF+R3wiJZV8o+ZHQfcjjfU2AjvvsKfgenASOfcdt/ChajAFfLhwFCgFbAN7wrhvZpqpOYxszOAa/CWF0sB8vGuHL0OPBGY41F8YGYL8JYmLcnC4qskmVk7vL9RPfEKxq+A+7VkXOVQESgiIiISgnRPoIiIiEgIUhEoIiIiEoJUBIqIiIiEIBWBIiIiIiFIRaCIiIhICFIRKCIiIhKCVASKiIiIhCAVgSIiIiIhSEWgiIiISAhSESgiUk3M7HAzyzOz3n5nKcrMBphZrpm18TuLiFQfLRsnIlJNzGwuEO2cK23t1Mo4x7PAAKClC+IXvJktxVtf/KyqyiYiNYuuBIqIVAMz6wb0Bp6ownMYcAYwPZgCMOAp4Ewza1/5yUSkJlIRKCJSPa4BUoFZVXiOrkAL4O296DsNyAKuqsxAIlJzqQgUkZBnZi+YmTOzRiXsOzhwv9xz+3D8CLwh2nnOuT3F9l0aOPeJZnavmf1iZtlm9rmZHRNo09PMPjazTDPbZGb/LOVUZwI7gA8D/WLM7H4zW21mWWa23cy+NbPHind0zmUAi4Bz9/Z9ikjtoiJQRAS+CXw8vIR9jwA5wH37cPzOQALwRRltRuIVik8BDwAHAnPMbADeVbpFwC3AKuBBM7uwhGOcCcwsUmg+G8j9GXATcDfwAXBCKRk+BZqYWduKvjERqb0i/A4gIlIDFC0CPyzcaGbHAmcBdznntuzD8Q8LfFxTRptw4BjnXG7g3N8B04E3gG7OuS8D2/8L/AJcC/yvSNZ2wKHAPUWOeSbwnnPukgrmLMzXHq/YFJE6TFcCRURgeeBj8SuBjwPrgSf38fiFw8xpZbR5rrAADFgU+PhZYQEIEGjzBVB8OpcBwG5gdpFtO4D2ZlbSFc6SpAY+Nq5gexGpxVQEikjIc86lApvwroABYGbnAd2AO51zOft6isLDltFmbbFM6YFPfy6hbTqQUmzbmXj3HGYU2TYcSAK+NbM1ZjbBzM4ws9J+9xfm09xhIiFARaCIiOcbAkWgmUUBI4AvgVcq4dhbAx+Ty2iTH+T2P5hZS6ALxZ4Kds5NB1oBFwHzgRMDbRYE3mNxhfm2lrBPROoYFYEiIp5vgMRAQXUt3oMZNxedb8/MbjezZYH5+DCzfma2uQJz660IfKyqFTkG4F29e6f4DudcmnPuf865IXjv6VGgO958gsUdHPi4ooR9IlLHqAgUEfEUPhzSHe/himnOuUXF2jyFNwx7bmDy50nAGc65leUc+2tgJ3BMJeYt6kzgY+fcH1fwzCzczBKLNgoUtF8HvizpquQxwGbn3OoqyikiNYieDhYR8RQWgaOBesDtxRs453LM7F686VxigYucc5+Vd2DnXL6ZTQPOMLNo59zuygptZklAD+C2YrvqAZvM7B28wm8L0Bq4Gu+ewneLHScBrwB+obKyiUjNpiuBIiKe74FcvCdjn3XO/VRKu6/w7rN7zTk3J4jjP4f3kMap+xKyBKfh/Yf+7WLbs/AK2tbArYHzX4Q3ZHy0c25jsfZnA3HAuErOJyI1lAW/vKSISGgys9bAR8BbwPnAQc65XUH0nw3EO+e6V2Kmt4DWzrmO+3icpcAvzrmzKiWYiNR4uhIoIlIBZtYEmAs87Jy7HvgO7wpbMG4GupnZyZUY7VPgrn05QGBVkr9RwhC4iNRduhIoIlIOM6sPLASmO+fuD2w7DpgDtHHObfIxnojIXlERKCIiIhKCNBwsIiIiEoJUBIqIiIiEIBWBIiIiIiFIRaCIiIhICFIRKCIiIhKCVASKiIiIhCAVgSIiIiIhSEWgiIiISAj6Pxn6G+E5KrnxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "k_B = 1.38e-23\n",
    "T = 450 + 273.15\n",
    "m = 174 * 1.67e-27\n",
    "\n",
    "sigma = 2 * 0.017 * np.sqrt(k_B * T / (2 * m))\n",
    "# assuming width of the transverse velocity is 2 x 17mrad x width of longitudinal v\n",
    "\n",
    "# sigma = 4.93 # estimate from the distribution of v_r and v_z\n",
    "\n",
    "vx = np.random.normal(0, sigma, 500)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "count, bins, ignored = plt.hist(vx, 40, density=True)\n",
    "plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) *\n",
    "               np.exp( - (bins - 0)**2 / (2 * sigma**2) ),\n",
    "         linewidth=2, color='k')\n",
    "ax.set_xlabel(\"$v_x$ (m/s)\")\n",
    "ax.set_ylabel(\"Frac. of atoms\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ef123bb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "vx = np.array([0, 0.1]) # 2 atoms\n",
    "np.save(\"./configs/vx_input_2.npy\", vx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8588f366",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0. , 0.1])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84604aca",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save(\"./configs/bz_frames_2_new.npy\", b_new)"
   ]
  }
 ],
 "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}