{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "7520724c-2be0-44a3-b657-4247cd9c497d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import configparser"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7b936622-1122-438e-a7d6-c302922aaadc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b'#Units: Masses in Msun / h\\n'\n",
      "b'#Units: Positions in Mpc / h (comoving)\\n'\n",
      "b'#Units: Velocities in km / s (physical, peculiar)\\n'\n",
      "b'#Units: Halo Distances, Lengths, and Radii in kpc / h (comoving)\\n'\n",
      "b'#Units: Angular Momenta in (Msun/h) * (Mpc/h) * km/s (physical)\\n'\n",
      "b'#Units: Spins are dimensionless\\n'\n",
      "b'#Units: Total energy in (Msun/h)*(km/s)^2 (physical)\\n'\n"
     ]
    }
   ],
   "source": [
    "dirname = '/data54/lavaux/VELMASS/halo_central/halos_new/'\n",
    "lam = 5\n",
    "R_max = 100\n",
    "Mmin = 3e12\n",
    "    \n",
    "with open(dirname + \"/halos_0.0.ascii\", \"rb\") as f:\n",
    "    for _ in range(30):\n",
    "        r = str(f.readline())\n",
    "        if \"Units\" in r:\n",
    "            print(r)\n",
    "\n",
    "# Get box size\n",
    "with open(dirname + '/auto-rockstar.cfg') as f:\n",
    "    data = [r for r in f]\n",
    "Lbox = [r for r in data if r.startswith('BOX_SIZE')][0].strip()\n",
    "Lbox = float(Lbox.split('=')[1])\n",
    "origin = np.array([Lbox/2, Lbox/2, Lbox/2])    \n",
    "omega_m = [r for r in data if r.startswith('Om')][0].strip()\n",
    "omega_m = float(omega_m.split('=')[1])\n",
    "\n",
    "halo_file = dirname + '/out.npy'\n",
    "halos = np.load(halo_file)\n",
    "\n",
    "# Cut by mass\n",
    "# Mmin = ???\n",
    "halos = halos[halos['Mvir'] > Mmin]\n",
    "\n",
    "xtrue = halos['X']  - origin[0] # Mpc / h\n",
    "ytrue = halos['Y']  - origin[1] # Mpc / h\n",
    "ztrue = halos['Z']  - origin[2] # Mpc / h\n",
    "rtrue = np.sqrt(xtrue**2 + ytrue**2 + ztrue**2) # Mpc / h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "2c5f507e-348f-4e12-904f-3e8dd8455a58",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "12449803\n",
      "2000.0\n",
      "\tMade 342 of 400\n",
      "\tMade 400 of 400\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x149298041100>]"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAGdCAYAAAASUnlxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbrElEQVR4nO3de1yO9/8H8NfdXd2FCkUHkkJObQ7l1CSHKTFzFiZnk3M1m9MOxrawRqL0c545tQkzMuVQQkiK0Bwjh1oKFdHd4fr90de93St0p7qq+/V8PK7H1nW97+t63R+pt+u+rs8lEQRBABEREZEa0BA7ABEREVFFYeNDREREaoONDxEREakNNj5ERESkNtj4EBERkdpg40NERERqg40PERERqQ02PkRERKQ2NMUOUJkUFBTg4cOH0NPTg0QiETsOERERlYAgCMjKyoKZmRk0NN58ToeNz788fPgQ5ubmYscgIiKiUrh37x4aNmz4xho2Pv+ip6cHoHDg9PX1RU5DREREJZGZmQlzc3PF7/E3YePzL68+3tLX12fjQ0REVMWU5DIVXtxMREREaoONDxEREakNNj5ERESkNtj4EBERkdpg40NERERqg40PERERqQ02PkRERKQ22PgQERGR2mDjQ0RERGqDjQ8RERGpDTY+REREpDbY+BAREZHaYONDREREaoONDxEREakNTbEDEBGVxMqw66V6nWdv6zJOQkRVGc/4EBERkdpg40NERERqo1SNT0BAACwtLaGjowNbW1tERka+sT4iIgK2trbQ0dGBlZUVAgMDlbZfuXIFQ4YMQePGjSGRSODr61tkH6+2/XeZPn26ombcuHFFtnfu3Lk0b5GIiIiqIZUbn6CgIHh4eGDhwoWIjY2Fg4MDXFxckJSUVGx9YmIi+vbtCwcHB8TGxmLBggWYNWsWgoODFTXZ2dmwsrLC0qVLYWJiUux+oqOjkZycrFjCwsIAAMOGDVOq69Onj1JdSEiIqm+RiIiIqimVL25esWIFJk6ciEmTJgEAfH19cfjwYaxduxbe3t5F6gMDA9GoUSPFWZyWLVvi/Pnz8PHxwZAhQwAAHTp0QIcOHQAA8+bNK/a49erVU/p66dKlaNKkCRwdHZXWy2Sy1zZPREREpN5UOuMjl8sRExMDJycnpfVOTk44ffp0sa+JiooqUu/s7Izz588jNzdXxbj/5Ni2bRsmTJgAiUSitC08PBz169eHtbU1Jk+ejNTU1NfuJycnB5mZmUoLERERVV8qNT5paWnIz8+HsbGx0npjY2OkpKQU+5qUlJRi6/Py8pCWlqZi3EL79u3D06dPMW7cOKX1Li4u2L59O44dO4affvoJ0dHR6NmzJ3Jycordj7e3NwwMDBSLubl5qfIQERFR1VCqeXz+e5ZFEIQi695WX9z6ktq4cSNcXFxgZmamtN7V1VXx/zY2NrCzs4OFhQUOHjyIwYMHF9nP/Pnz4eXlpfg6MzOTzQ8REVE1plLjY2RkBKlUWuTsTmpqapGzOq+YmJgUW6+pqQlDQ0MV4wJ3797FkSNHsGfPnrfWmpqawsLCAjdu3Ch2u0wmg0wmUzkDERERVU0qfdSlra0NW1tbxR1Vr4SFhcHe3r7Y13Tp0qVIfWhoKOzs7KClpaViXGDz5s2oX78++vXr99ba9PR03Lt3D6ampiofh4iIiKoflW9n9/LywoYNG7Bp0yYkJCTA09MTSUlJcHd3B1D48dGYMWMU9e7u7rh79y68vLyQkJCATZs2YePGjZgzZ46iRi6XIy4uDnFxcZDL5Xjw4AHi4uJw8+ZNpWMXFBRg8+bNGDt2LDQ1lU9WPXv2DHPmzEFUVBTu3LmD8PBw9O/fH0ZGRhg0aJCqb5OIiIiqIZWv8XF1dUV6ejoWL16M5ORk2NjYICQkBBYWFgCA5ORkpTl9LC0tERISAk9PT/j7+8PMzAx+fn6KW9kB4OHDh2jXrp3iax8fH/j4+MDR0RHh4eGK9UeOHEFSUhImTJhQJJdUKkV8fDy2bt2Kp0+fwtTUFD169EBQUBD09PRUfZtERERUDUmEV1caEzIzM2FgYICMjAzo6+uLHYeI/oUPKSWi11Hl9zef1UVERERqg40PERERqQ02PkRERKQ22PgQERGR2mDjQ0RERGqDjQ8RERGpDTY+REREpDbY+BAREZHaYONDREREaoONDxEREakNNj5ERESkNtj4EBERkdpg40NERERqg40PERERqQ02PkRERKQ22PgQERGR2mDjQ0RERGpDU+wApN5Whl0v1es8e1uXcZLqoTTjybEkInXCMz5ERESkNtj4EBERkdpg40NERERqg40PERERqQ02PkRERKQ22PgQERGR2mDjQ0RERGqDjQ8RERGpDTY+REREpDbY+BAREZHaYONDREREaoONDxEREakNNj5ERESkNtj4EBERkdrQFDsAEVF1sjLseqle59nbuoyTEFFxeMaHiIiI1AYbHyIiIlIbbHyIiIhIbbDxISIiIrXBxoeIiIjURqkan4CAAFhaWkJHRwe2traIjIx8Y31ERARsbW2ho6MDKysrBAYGKm2/cuUKhgwZgsaNG0MikcDX17fIPhYtWgSJRKK0mJiYKNUIgoBFixbBzMwMurq66N69O65cuVKat0hERETVkMqNT1BQEDw8PLBw4ULExsbCwcEBLi4uSEpKKrY+MTERffv2hYODA2JjY7FgwQLMmjULwcHBiprs7GxYWVlh6dKlRZqZf2vdujWSk5MVS3x8vNL25cuXY8WKFVizZg2io6NhYmKC3r17IysrS9W3SURERNWQyo3PihUrMHHiREyaNAktW7aEr68vzM3NsXbt2mLrAwMD0ahRI/j6+qJly5aYNGkSJkyYAB8fH0VNhw4d8OOPP2LEiBGQyWSvPbampiZMTEwUS7169RTbBEGAr68vFi5ciMGDB8PGxgY///wzsrOzsWPHDlXfJhEREVVDKjU+crkcMTExcHJyUlrv5OSE06dPF/uaqKioIvXOzs44f/48cnNzVQp748YNmJmZwdLSEiNGjMDt27cV2xITE5GSkqJ0LJlMBkdHx9dmy8nJQWZmptJCRERE1ZdKjU9aWhry8/NhbGystN7Y2BgpKSnFviYlJaXY+ry8PKSlpZX42J06dcLWrVtx+PBhrF+/HikpKbC3t0d6erriOK/2XdJs3t7eMDAwUCzm5uYlzkNERERVT6kubpZIJEpfC4JQZN3b6otb/yYuLi4YMmQI3nvvPXz44Yc4ePAgAODnn38udbb58+cjIyNDsdy7d6/EeYiIiKjqUelZXUZGRpBKpUXOoKSmphY50/KKiYlJsfWampowNDRUMe4/atasiffeew83btxQHAcoPPNjampaomwymeyN1xQRERFR9aLSGR9tbW3Y2toiLCxMaX1YWBjs7e2LfU2XLl2K1IeGhsLOzg5aWloqxv1HTk4OEhISFE2OpaUlTExMlI4ll8sRERHx2mxERESkXlR+OruXlxfc3NxgZ2eHLl26YN26dUhKSoK7uzuAwo+PHjx4gK1btwIA3N3dsWbNGnh5eWHy5MmIiorCxo0bsXPnTsU+5XI5rl69qvj/Bw8eIC4uDrVq1ULTpk0BAHPmzEH//v3RqFEjpKam4rvvvkNmZibGjh0LoPAjLg8PD/zwww9o1qwZmjVrhh9++AE1atTAqFGj3m2UiIiIqFpQufFxdXVFeno6Fi9ejOTkZNjY2CAkJAQWFhYAgOTkZKU5fSwtLRESEgJPT0/4+/vDzMwMfn5+GDJkiKLm4cOHaNeuneJrHx8f+Pj4wNHREeHh4QCA+/fvY+TIkUhLS0O9evXQuXNnnDlzRnFcAPjiiy/w4sULTJs2DU+ePEGnTp0QGhoKPT09lQeG6N9Whl0v1es8e1uXcRKif/D7kkh1Kjc+ADBt2jRMmzat2G1btmwpss7R0REXLlx47f4aN26suOD5dXbt2vXWXBKJBIsWLcKiRYveWktERETqh8/qIiIiIrXBxoeIiIjUBhsfIiIiUhtsfIiIiEhtsPEhIiIitcHGh4iIiNQGGx8iIiJSG2x8iIiISG2w8SEiIiK1wcaHiIiI1AYbHyIiIlIbbHyIiIhIbbDxISIiIrXBxoeIiIjUhqbYAah6WBl2XewIREREb8UzPkRERKQ22PgQERGR2mDjQ0RERGqDjQ8RERGpDTY+REREpDbY+BAREZHa4O3sRGom58VzPLiVgLSHSch4lIJrwS9x//59PHz4EC9fvoRcLkdOTg7kcjk0NDSgr68PAwMDxX/Nzc1hZWWltOjo6Ij9toiISoSND1E1JggC/r57E9djo3D/xmXcu34ZqfduQxCEEu8jJSXljdulUilat26NDh06KJb3338fmpr88UJElQ9/MhFVM/Kcl7h18Syung3H1XMRePL3gyI1+ob1YWLRDLXrmaBPp9Zo2LAhzMzMULNmTWhra0NbWxsymQz5+fnIzMxEZmYmMjIy8OTJEyQlJeH27du4ffs2bt26hczMTFy6dAmXLl3Cxo0bC/evr49evXrByckJTk5OsLKyquhhICIqFhsfompAEATExMTgt1U+uHD8AHKynyu2aWppo0mbjmjcsh3MrW3QsFlr6Netp9ju2dv6nY774MEDREdH4/z584iOjkZ0dDSePn2KvXv3Yu/evQCApk2bYujQoXB1dUWbNm0gkUhK/2aJiN4BGx+iKiwzMxO//PILNmzYgLi4OMX62kYmaNnJEa06dkfTtp0h061RLseXSCRo2LAhGjZsiEGDBgEA8vPzceHCBYSGhiI0NBSnT5/GzZs3sXTpUixduhTW1tYYPnw4PvnkE7Ro0aJcchERvQ4bH6IqKCMjA35+fli5ciWePHkCAJDJZGht3xudXIahyfsdoaEhzk2bUqlUca3PwoULkZWVhUOHDuHXX3/FwYMHcf36dXz33Xf47rvv4ODggE8//RRDhgyBrq6uKHmJSL3wdnaiKuTp06dYtGgRLCws8PXXX+PJkydo3rw5fH198fDhQ4ye/xOate0sWtNTHD09PQwfPhy7d+9Gamoqtm/fjo8++ghSqRSRkZFwc3NDgwYN4OHhgZs3b4odl4iqucrz05GIXisvLw+rV69G48aN8e233yIjIwOtWrXCzp07ceXKFcyePRt169YVO+Zb6enpYdSoUfjjjz9w9+5dLFmyBI0aNcKTJ0+watUqWFtbY+jQoThz5ozYUYmomuJHXaRWVoZdFztCifw7Z+KVGASvXoyHt/8CAJg0tobT6Gl4v6szkjU04HfsllgxS+Xf761mp+GYZTcE1y6cQuS+X/BX9AkEBwcjODgYlja26DlsElp17iHKxdBV5XuFiFTDxoeoksp6ko4DG30QHboHAKCrZ4B+473Q2WUYNKRSkdOVHQ2pFC07dEPLDt2QnHgdEXs2I+boH0i8HIONl2Ngbm0DZ7eZaNnRkXeDEdE7Y+NDVAldOXMcQT8twLOMxwCATn2Got/Ez1DLoPJ/nPUuTC2tMeIzb7iM88CJvVtxav8O3Lt+GRu+moJGzd+H85iZaGHnwAaIiEqN1/gQVSLZ2dmYPn06Nn7tjmcZj2FqaY1Zq4Lg6vV9tW96/s3A0Bj9J32OL7ceRY9hE6Et00XStUtYv3AyAj4fg3vX48WOSERVFBsfokoiLi4OdnZ2CAgIAAA4Dh4Hj9W70bhlW3GDiahW7broP/kLLNx6BI5DxkNTW4Zbl85h5Yyh2LZ0Dh4XMys1EdGbsPEhqgQ2bNiATp06ISEhASYmJpjyw0YMcJ8PLW2Z2NEqBb06RhgwZR7mb/oTdh8OAABcOPYHlk7ogwMbfZDz4vlb9kBEVIiND5GI8vLy4OHhgcmTJ0Mul+Pjjz9GfHw8mtt1FTtapVSnvhlGfbEcXv570LRNJ+TlynEsaD2WTnRBbHiISg9fJSL1xMaHSCRPnjxB3759sWrVKgDA4sWLsW/fPhgZGYmcrPJr2Kw1pi7/GRO/XQtDU3NkpP2NX37wRODccUi5y0kQiej12PgQieDatWvo3LkzwsLCUKNGDezevRtfffUV71ZSgUQiQesuPfHF+oNwHjMTmtoy3Ig7Ax/3ATiw8SfIc16KHZGIKqFSNT4BAQGwtLSEjo4ObG1tERkZ+cb6iIgI2NraQkdHB1ZWVggMDFTafuXKFQwZMgSNGzeGRCKBr69vkX14e3ujQ4cO0NPTQ/369TFw4EBcu3ZNqWbcuHGQSCRKS+fOnUvzFonKzblz52Bvb4/r16+jUaNGOHXqFIYMGSJ2rCpLS1sG59EzMHf9QbTu0hMF+Xk4FrQOPlM+xo04zgBNRMpUbnyCgoLg4eGBhQsXIjY2Fg4ODnBxcUFSUlKx9YmJiejbty8cHBwQGxuLBQsWYNasWQgODlbUZGdnw8rKCkuXLoWJiUmx+4mIiMD06dNx5swZhIWFIS8vD05OTnj+XPmixj59+iA5OVmxhISEqPoWicpNeHg4evXqhcePH6Njx444d+4c2rZtK3asasHQ1BwTv12L8d/4Q9+wPtIe3sXaL8Zi8uTJige5EhGp3PisWLECEydOxKRJk9CyZUv4+vrC3Nwca9euLbY+MDAQjRo1gq+vL1q2bIlJkyZhwoQJ8PHxUdR06NABP/74I0aMGAGZrPi7WP7880+MGzcOrVu3Rps2bbB582YkJSUhJiZGqU4mk8HExESxVIXnF5F6OHjwIFxcXPDs2TP07NkTR48ehbGxsdixqp33PvgQczeEwP6jEQAK75hr1aoV/vjjD5GTEVFloFLjI5fLERMTAycnJ6X1Tk5OOH36dLGviYqKKlLv7OyM8+fPIzc3V8W4/8jIyACAIo1NeHg46tevD2tra0yePBmpqamlPgZRWQkKCsLAgQPx8uVL9O/fHwcPHkStWrXEjlVt6dbUw9BZ32LGT9vRvHlzpKSk4OOPP8b48eMVPzuISD2p1PikpaUhPz+/yL9SjY2NkZKSUuxrUlJSiq3Py8tDWlqainELCYIALy8vdO3aFTY2Nor1Li4u2L59O44dO4affvoJ0dHR6NmzJ3JycordT05ODjIzM5UWorK2detWjBw5Enl5eRg1ahSCg4Oho6Mjdiy1YPWeHeLi4jBnzhxIJBJs2bIFNjY2CA0NFTsaEYmkVBc3//fOE0EQ3ng3SnH1xa0vqRkzZuDSpUvYuXOn0npXV1f069cPNjY26N+/Pw4dOoTr16/j4MGDxe7H29sbBgYGisXc3LxUeYheZ8+ePRg/fjwEQcCUKVPwyy+/QEtLS+xYakVHRwc//vgjIiMj0bRpU9y/fx/Ozs6YNm0asrOzxY5HRBVMpYeUGhkZQSqVFjm7k5qa+tprFUxMTIqt19TUhKGhoYpxgZkzZ2L//v04ceIEGjZs+MZaU1NTWFhY4MaNG8Vunz9/Pry8vBRfZ2ZmsvmhMnPt/EnM/cYdBQUFGD9+PAICAqChwRkkxPLBBx8gLi4O8+bNw5o1a7B27VocP34cO3bsQLt27cSOh5Vh18WOQKQWVPoprK2tDVtbW4SFhSmtDwsLg729fbGv6dKlS5H60NBQ2NnZqfQvX0EQMGPGDOzZswfHjh2DpaXlW1+Tnp6Oe/fuwdTUtNjtMpkM+vr6SgtRWbh9+Tw2fTsdubm5GDp0KNavX8+mpxKoWbMmVq9ejdDQUJiamuKvv/5Cp06d4OPjg4KCArHjEVEFUPknsZeXFzZs2IBNmzYhISEBnp6eSEpKgru7O4DCsyhjxoxR1Lu7u+Pu3bvw8vJCQkICNm3ahI0bN2LOnDmKGrlcjri4OMTFxUEul+PBgweIi4vDzZv/zMA6ffp0bNu2DTt27ICenh5SUlKQkpKCFy9eAACePXuGOXPmICoqCnfu3EF4eDj69+8PIyMjDBo0qNQDRKSq+zeuYMOXU5Cb8xJ9+vTB9u3bIZVKxY5F/9K7d29cunQJAwYMQG5uLj7//HM4OTnh4cOHYkcjonKmcuPj6uoKX19fLF68GG3btsWJEycQEhICCwsLAEBycrLSnD6WlpYICQlBeHg42rZtiyVLlsDPz09pwraHDx+iXbt2aNeuHZKTk+Hj44N27dph0qRJipq1a9ciIyMD3bt3h6mpqWIJCgoCAEilUsTHx2PAgAGwtrbG2LFjYW1tjaioKOjp6ZV6gIhUkZ58D+sWTMLL7GewsrFDcHAwtLW1xY5FxTAyMsLevXvxf//3f9DV1cXRo0fRtm1bXvhMVM2pdI3PK9OmTcO0adOK3bZly5Yi6xwdHXHhwoXX7q9x48Zvfbjg27br6uri8OHDb6whKk8vnmdhw1dT8CzjMRo0bYVJS/4PNWrUEDsWvYFEIsGnn36Kbt26wdXVFZcuXYKzszPmz5+Put1GQyot1Y9IIqrEeNEBURnIz8/D1u888HfSLegb1sfExYHQqcl5eqqKFi1a4MyZM5g6dSqAwjs+A+a44UlqssjJiKissfEhekeCIGBfwPe4FnMS2jJdTFociNpGnJG5qtHV1UVAQAB+/fVX6OvrI/HKBayYNhB/nX/zswiJqGph40P0jk7+vg2n/tgBiUSCT+b5oGGz1mJHoncwbNgwxMbGomGz1nie+RTrF07Gn1tXoyA/X+xoRFQG2PgQvYOE6BPYF/gDAKDfxDl474MPRU5EZcHKygozV+5El36uEAQBodvWYP2Xn+JZxmOxoxHRO2LjQ1RKj/9+gG1L50AoKEBH5yHoMWyi2JGoDGlpyzBs9mKM/HwZtGQ6uBZzEiumDUbStUtiRyOid8DGh6gU8uRy/PzdbLzIyoB58/cwdOaiUj+ChSq3Dr0HYrbfr6jXoDGePkrGGq9PcPbP3WLHIqJSYuNDVAq/r1uKe9fioatngLFfroIm5+qp1swsm8NjzW607tITeblyBK1YiN1+i5CXKxc7GhGpiI0PkYouHD+AU/u3AwA+mfsj6ho3EDkRVQTdmnoY/40/+oyZBYlEgtMHdiLg8zHITE8VOxoRqYCND5EK/k66hV9XfgUA6D1qKlp1dBQ5EVUkDQ0NOI2e/r95mvRw52osVkwfjLsJF8WORkQlxMaHqITkL19gy5JZkL/MRrO2neHsNlPsSCSSVp26w3PNbphYNEPm40fwnzMa0WH7xI5FRCXAxoeohP5Yvxx/370J/br1MHr+CmjwwaNqrV6Dxpi1ahdsuvRCXq4cO3+ci/3rliE/P0/saET0Bmx8iEog4VwETv2xAwAw6ovl0KtjKHIiqgx0atTCuG/WoPeowkddhO/ehA1fTkF2VobIyYjoddj4EL3Fs4zH2LViIQDAYdAYWLe3FzkRVSYaGhpwGeeBMV/6Kub78ZvtikcP7ogdjYiKwcaH6A0EQcBvvl8j6/EjGFs0Rb8Jn4kdiSqptt1cMHPlTtSuZ4rU+4nwnTUcN+LOiB2LiP6DjQ/RG0SH7kH8qTBINbUwep4PtGU6YkeiSqxh01bwWP0bGrVogxdZGfi/+RMRdTBI7FhE9C9sfIheIz35HvYGfAcA6DN2Fho0aSlyIqoK9OvWw3SfX9C+x0coyM/Db6u+xt613/Mhp0SVBBsfomIUFBRg108LkPMiG1Y2dugxlM/hopLT0pbhk3k+cBnnAQCI3LsVmxZNw8vsZ+IGIyI2PkTFOXvoN9y6dA7aMl2M/Hwpb10nlUkkEvQeNfV/jzSR4erZcKzx+gRPH6WIHY1IrbHxIfqPp2l/44/1ywEALuM9YGhqLnIiqsradOuD6T7boFfHCA9v/wXfmUNx7/plsWMRqS02PkT/IggC9qz+Fi+zn6FRizZwGOAmdiSqBixavI/Zfr/+M9PzZ6Nx+fQRsWMRqSU2PkT/cinyMC5HHYWGVBOunt/xIy4qM3WNG2Cm7040t+0Kec4LbP52Bk7s3Sp2LCK1w8aH6H+eZz7FHv8lAIBeIz6FqaW1yImoutGtqYdJ3/0fOvd1hSAI2Lf2e97xRVTB2PgQ/c8f65ch60kajBs1Qe+RU8WOQ9WUVKqJYbO/xUeT5gAovOOr8OG3L0RORqQe2PgQAbh58SzOHd4DiUSC4Z7fQVNbW+xIVI1JJBL0HD4ZbgtWQlNLG5dPH4H/527IepImdjSiao+ND6m9/Lxc7FlT+BFXl36usGzdXuREpC7ade8L92WbUUOvNu5di8eq2a5IvZ8odiyiao2ND6m9k/u3I+XuDdTUrw2XcZ5ixyE1Y2Vjh1mrdsHQ1ByPU+5jtccIJF65IHYsompLU+wARGLKfPwIh7euBgD0mzgHNfVrixuIytzKsOtiR3ir+g0tMcs3CBu+noJ71+Kxdu44jJ7ng/e7OokdrUoq7Z+5Z2/e0KAOeMaH1NqBDT54mf0M5s3fQ0fnIWLHITWmV8cQ05ZvRevOPZAnz8HPS2bxdneicsDGh9RW4pUYnD+yDxKJBENmfA0NDf51IHHJdGtg3DdrYP/RSMXt7n+sX46CggKxoxFVG/xJT2qpID9fcUFzxz5D0aj5+yInIioklWpiyMxv0He8FwDg+G8bsWPZ58jLlYucjKh6YONDainq4C48uJUA3Vr66Pe/XzBElYVEIsGHI6dg5OfLoCHVxIXjB7D+y0/x4nmW2NGIqjw2PqR2srMycOhnPwCAyzgP1KpdV+RERMXr0HsgJi/5P8h0a+BGbBTWfPYJMtL/FjsWUZXGxofUzpGdgcjOegoTi2bo0s9V7DhEb9Tcrium/7QdenXrIfn2Nfh5jMTfSbfEjkVUZbHxIbWSnnwPkb//AgDoP/lzSKWc0YEqv4ZNW2HWyp2o17Axnvz9AKs9R+HO1VixYxFVSWx8SK2EbF6J/NxcWLezR4sO3cSOQ1RihqbmmLliJxo1fx/ZWU+xdu44XDlzXOxYRFUOGx9SG2fPnkVs+EFIJBL0//QLSCQSsSMRqaRW7bqYuvxntOzoiNycl9i0aBrOHPpN7FhEVQobH1ILgiBgzpzCp2Hb9R6EBk1aipyIqHRkujUwYZE/OjoPhlBQgF9Xfomw7QEQBEHsaERVAhsfUgv79u3DyZMnoSXTgcvY2WLHIXonUk0tuHr9gA9HugMADv28Cnv8l6AgP1/kZESVHxsfqvbkcjm++OILAED3IeNRu56JyImI3p1EIkHf8Z4YNO1LSCQSnNq/Hb/84IlceY7Y0YgqtVI1PgEBAbC0tISOjg5sbW0RGRn5xvqIiAjY2tpCR0cHVlZWCAwMVNp+5coVDBkyBI0bN4ZEIoGvr2+pjisIAhYtWgQzMzPo6uqie/fuuHLlSmneIlUjGzZswM2bN1G/fn30GD5J7DhEZcphoBvcFqyAVEsLFyMPY93CSZzokOgNVG58goKC4OHhgYULFyI2NhYODg5wcXFBUlJSsfWJiYno27cvHBwcEBsbiwULFmDWrFkIDg5W1GRnZ8PKygpLly6FiUnx/xovyXGXL1+OFStWYM2aNYiOjoaJiQl69+6NrCz+EFBX2dnZWLKk8NEUX3/9NXRq1BI5EVHZa+vYF59+tx6yGjVx6+I5+H82GpnpqWLHIqqUVG58VqxYgYkTJ2LSpElo2bIlfH19YW5ujrVr1xZbHxgYiEaNGsHX1xctW7bEpEmTMGHCBPj4+ChqOnTogB9//BEjRoyATCYr1XEFQYCvry8WLlyIwYMHw8bGBj///DOys7OxY8cOVd8mVRP+/v5ISUlB48aNMXnyZLHjEJWbZu26YMZP26FXxwgPb/8FP8+RePTgjtixiCodlRofuVyOmJgYODk5Ka13cnLC6dOni31NVFRUkXpnZ2ecP38eubm5ZXbcxMREpKSkKNXIZDI4Ojq+NltOTg4yMzOVFqo+MjMzsXTpUgDAN998A21tbZETEZWvBk1aYubKnTA0a4THKfex2mMk7l2/LHYsokpFpcYnLS0N+fn5MDY2VlpvbGyMlJSUYl+TkpJSbH1eXh7S0tLK7Liv/qtKNm9vbxgYGCgWc3PzEuWhqmHlypV4/PgxmjdvjtGjR4sdh6hCGJk1wqyVO9GgaSs8y3iMgM/dcP1C8f/4I1JHpbq4+b8TvwmC8MbJ4IqrL259WRxXlWzz589HRkaGYrl3755KeajySk9Px08//QQAWLx4MTQ1+WgKUh96dYww/cdf0KxtZ+S8yMb6Lz9FbHiI2LGIKgWVGh8jIyNIpdIiZ1BSU1OLnGl5xcTEpNh6TU1NGBoaltlxX10UrUo2mUwGfX19pYWqh+XLlyMrKwtt2rTB0KFDxY5DVOF0atbC5O/Wo023PsjPy8U2by9E/r5N7FhEolOp8dHW1oatrS3CwsKU1oeFhcHe3r7Y13Tp0qVIfWhoKOzs7KClpVVmx7W0tISJiYlSjVwuR0RExGuzUfWUnJyM1atXAwC+++47aGhwuipST5ra2nCbvwIf9B8FQRCw138JDv28irM8k1pT+fy/l5cX3NzcYGdnhy5dumDdunVISkqCu3vhDKLz58/HgwcPsHXrVgCAu7s71qxZAy8vL0yePBlRUVHYuHEjdu7cqdinXC7H1atXFf//4MEDxMXFoVatWmjatGmJjiuRSODh4YEffvgBzZo1Q7NmzfDDDz+gRo0aGDVq1LuNElUpP/zwA168eIHOnTujX79+YschEpWGVIrBM75GrTqGOLx1NcK2B8C9Zh4CAgIglUrFjkdU4VRufFxdXZGeno7FixcjOTkZNjY2CAkJgYWFBYDCf23/e24dS0tLhISEwNPTE/7+/jAzM4Ofnx+GDBmiqHn48CHatWun+NrHxwc+Pj5wdHREeHh4iY4LAF988QVevHiBadOm4cmTJ+jUqRNCQ0Ohp6en8sBQ1fTgwQOsW7cOAPD999/zQaREKPyHofPoGdCrbYjg1d9i3bp1SEtLw/bt26GjoyN2PKIKJRF4zlMhMzMTBgYGyMjI4PU+KloZdr1Cj+fZ27rY9bNnz4afnx+6deuGiIiIItsrS863qcicpc1YWhX9Z0DKLkYexs5lcyCXy9G9e3fs27cPBgYGYscqU6X9HqvovwtUdlT5/c2LH6jaSElJUZzt+frrr0VOQ1Q5tXFwxp9//gk9PT2Eh4eje/fur53yg6g6YuND1YaPjw9evnyJLl26oGfPnmLHIaq0evTogYiICBgbGyMuLg4ffPABbt26JXYsogrBxoeqhUePHikeX/LVV1/x2h6it2jXrh1OnToFKysr3L59G/b29oiNjRU7FlG5Y+ND1cKKFSuQnZ0NOzs79OnTR+w4RFVCkyZNcOrUKbRt2xapqalwdHTE8ePHxY5FVK7Y+FCV9/jxY6xZswYAz/YQqcrExERxrU9WVhb69OmD4OBgsWMRlRs2PlTl+fr64tmzZ2jTpg369+8vdhyiKsfAwACHDh3C4MGDIZfLMWzYMAQGBoodi6hcsPGhKi0jIwN+fn4AeLaH6F3o6Ojg119/xZQpUyAIAqZOnYpvv/2WszxTtcMnN1KVtmbNGmRkZKB169YYNGiQ2HGqJM6ro37e9GduPdgTTs81EbrNH4sWLULo+WsYPP0raLzDLM+cH4cqE57xoSrrxYsXWLVqFYDCR6XwmVxE704ikaDPmFkYMuNrSCQSnD6wE1u/90SuPEfsaERlgr8pqMrasmULHj16BAsLC7i6uoodh6ha+eDjTzBmoS+kWlq4dPIw1i2chBfPs8SORfTO2PhQlZSXlwcfHx8AwGeffQZNTX5qS1TW2nTrg0+/Ww9ZjZq4dfEcAua4IfPxI7FjEb0TNj5UJQUHB+P27dswNDTEhAkTxI5DVG01a9cF0322Qa+OER7cSoCfxwg8enBX7FhEpcbGh6ocQRCwbNkyAMDMmTNRs2ZNkRMRVW8Nm7bCzJU7YWjWCI9T7mO1xwjcu35Z7FhEpcLGh6qc6xdOIzY2FjVq1MCMGTPEjkOkFozMGmHWyp1o0LQVnmU8RsDnbrh+4bTYsYhUxsaHqpxjv64HAEyaNAmGhoYipyFSH3p1jDD9x1/QrF0X5LzIxvovP0Xs8YNixyJSCRsfqlLuXb+MG7FRkEql8PLyEjsOkdrRqVkLk5esQ1tHF+Tn5eIXby9E7NkidiyiEmPjQ1XKq7M9I0eOhIWFhchpiNSTprY2Rs9fga4D3AAAvwd6448NP6KgoEDkZERvx8aHqoz05Hu4dDIUAPD555+LnIZIvWloaGDQtIXoN+EzAMDxXzdgl8885OflipyM6M3Y+FCVcWLfVggFBWhu1xXvv/++2HGI1J5EIkGvEZ9ixBxvaGhIcf7I79jwtTtyXjwXOxrRa7HxoSrhxfMsnPtzNwDAcfA4ccMQkZKOToMxcfFaaMt0ce38SfjPcUPWkzSxYxEVi40PVQlnD/2GnBfZMLZoiua2XcWOQ0T/0bKjI6b+uBU1Derg/o0r8PMcyYkOqVJi40OVXn5+HiL3/QKg8GyPRCIRORERFceixfuYuXIn6po0RPrDpP9NdBgvdiwiJXzAERWxMuy62BGUxJ8Kw5PUh6hpUAfte/YHUPkyvk5VyUlUVuo3tMQs313Y8OUU3L95Bf5zxqC3VQ24uLiovC/+/aHywDM+VOmd2PMzAMD+o5HQlumInIaI3ka/bj1M89kK6/YfQP4yG/3798eWLVvEjkUEgI0PVXJ3EuJw52ospFpa+KD/KLHjEFEJ6dSohUlLAmHb62Pk5+dj/PjxWLJkCQRBEDsaqTk2PlSpnfjfjLDte/SHft164oYhIpVoamlj1BfLMW/ePADA119/jSlTpiAvL0/kZKTO2PhQpfUk9SEuRRZOWNht0FiR0xBRaUgkEnh7e8Pf3x8SiQTr16/HoEGD8Pw55/ohcbDxoUrr5O/bUFCQj2btuqBBkxZixyGidzBt2jQEBwdDR0cHBw4cQI8ePZCamip2LFJDbHyoUpK/fIEz/5uwkGd7iKqHQYMG4ejRozA0NER0dDS6dOmC69d55xZVLDY+VCldOH4AL7IyYGhqjpYduokdh4jKiL29PU6fPg0rKyvcvn1b8TVRRWHjQ5WOIAiI/L1wwsIP+o+ChlQqciIiKkvW1taIiopChw4dkJ6ejl69emHPnj1ixyI1wcaHKp3b8eeRfPsatGW66Og8ROw4RFQO6tevj+PHj6N///54+fIlhg4dCl9fX7FjkRpg40OVzsn/ne1p3+tj1NAzEDkNEZWXmjVrYs+ePZg6dSoEQYCnpydmz56N/Px8saNRNcbGhyqVp49SEH/qCACg64BPRE5DROVNU1MT/v7+WL58OQDAz88PQ4cORXZ2tsjJqLpi40OVyukDO1FQkI8m73eEmWVzseMQUQWQSCT4/PPPERQUBJlMhn379qFHjx7IepIudjSqhtj4UKWRK8/BmZBfAQBdB4wWOQ0RVbThw4fjyJEjqFu3Ls6dO4dVs4bh76RbYseiaoaND1UacRGH8CzjMWobmcDGvpfYcYhIBF27dkVUVBSaNGmCx38/gJ/HCNyIOyN2LKpGNMUOQPTKyf3bAAD2/UdCKuW3JpG6sra2xpkzZ9DB0Ql3rsbi/+ZPhKvnd+jgNEjsaMVaGab6JIyeva3LIQmVRKnO+AQEBMDS0hI6OjqwtbVFZGTkG+sjIiJga2sLHR0dWFlZITAwsEhNcHAwWrVqBZlMhlatWmHv3r1K2xs3bgyJRFJkmT59uqJm3LhxRbZ37ty5NG+RKtjdvy7h3rV4SLW00NlluNhxiEhkRkZGmLr8Z7R17IuC/Dzs9JmHQz+v4tPd6Z2p3PgEBQXBw8MDCxcuRGxsLBwcHODi4oKkpKRi6xMTE9G3b184ODggNjYWCxYswKxZsxAcHKyoiYqKgqurK9zc3HDx4kW4ublh+PDhOHv2rKImOjoaycnJiiUsLAwAMGzYMKXj9enTR6kuJCRE1bdIIog6uBMA0M6xL2rVrityGiKqDLS0ZRg9/yd8ONIdABC2PQDbls5BrjxH5GRUlanc+KxYsQITJ07EpEmT0LJlS/j6+sLc3Bxr164ttj4wMBCNGjWCr68vWrZsiUmTJmHChAnw8fFR1Pj6+qJ3796YP38+WrRogfnz56NXr15Kk1nVq1cPJiYmiuXAgQNo0qQJHB0dlY4nk8mU6urW5S/Ryi47KwOx4YUNqv1HI0VOQ0SViYaGBvqO94Sr1/fQkGoi9vgBrP1iLJ49fSx2NKqiVGp85HI5YmJi4OTkpLTeycnptc9aiYqKKlLv7OyM8+fPIzc39401r9unXC7Htm3bMGHCBEgkEqVt4eHhqF+/PqytrTF58uQ3Pv03JycHmZmZSgtVvOiwvcjNeQlTq+awaNlW7DhEVAl16jMUU37YAN1a+rhzNRa+s4Yh5e5NsWNRFaRS45OWlob8/HwYGxsrrTc2NkZKSkqxr0lJSSm2Pi8vD2lpaW+sed0+9+3bh6dPn2LcuHFK611cXLB9+3YcO3YMP/30E6Kjo9GzZ0/k5BR/WtTb2xsGBgaKxdzc/LXvncqHIAiIOrALAPDBRyOLNLJERK80a9cFs3x3wdDUHI9T7sPPYwSuxZwSOxZVMaW6uPm/v5wEQXjjL6zi6v+7XpV9bty4ES4uLjAzM1Na7+rqin79+sHGxgb9+/fHoUOHcP36dRw8eLDY/cyfPx8ZGRmK5d69e699D1Q+bsadQer9RMh0a6B9z/5ixyGiSs64URPM9vsVlq3b4+XzLKxfOBmnD+wUOxZVISo1PkZGRpBKpUXOxKSmphY5Y/OKiYlJsfWampowNDR8Y01x+7x79y6OHDmCSZMmvTWvqakpLCwscOPGjWK3y2Qy6OvrKy1UsU4fLDzbY9vrY+jUqCVyGiKqCmoZ1IX7si2w7fUxCgrysdtvEfYGfIf8/Dyxo1EVoFLjo62tDVtbW8UdVa+EhYXB3t6+2Nd06dKlSH1oaCjs7OygpaX1xpri9rl582bUr18f/fr1e2ve9PR03Lt3D6ampm+tpYqXmZ6qeC4XL2omIlVoacsw6ovl6DveEwAQue8XbPzKHS+eZ4mcjCo7lT/q8vLywoYNG7Bp0yYkJCTA09MTSUlJcHcvvN1w/vz5GDNmjKLe3d0dd+/ehZeXFxISErBp0yZs3LgRc+bMUdTMnj0boaGhWLZsGf766y8sW7YMR44cgYeHh9KxCwoKsHnzZowdOxaamsoT3D179gxz5sxBVFQU7ty5g/DwcPTv3x9GRkYYNKhyTnql7s7+uRsF+Xlo3KodzKxaiB2HiKoYiUSCD0e6Y+xXftCS6eCv85Hw8xiB9GRetkCvp3Lj4+rqCl9fXyxevBht27bFiRMnEBISAgsLCwBAcnKy0pw+lpaWCAkJQXh4ONq2bYslS5bAz88PQ4YMUdTY29tj165d2Lx5M95//31s2bIFQUFB6NSpk9Kxjxw5gqSkJEyYMKFILqlUivj4eAwYMADW1tYYO3YsrK2tERUVBT09PVXfJpWzgvx8RP3vuVw820NE76KNgzNm/LQd+ob18ffdm/CdORQ3L50TOxZVUhKB02AqZGZmwsDAABkZGWp9vU9ppl9X1eWoo9j0zTTU1K+Nr3ecgJa2rNyPSUTiKO3jGVT9WfQ07W9sXjQN965fhoZUE0NnLUJnl2Fvf+H/VFTOdzkWFU+V3998SCmJIupgEACgg9NgNj1EVCZqGxljus82xWMufl35Jfau/Z4XPZMSNj5U4Z6kPsRf0ScAAF36uYqchoiqE20dXbgtWIE+Y2YBACL3bsWGL6cgOytD5GRUWbDxoQp39s/dEAQBzdp2Rr0GjcWOQ0TVjEQigdPo6Rj7lR+0Zbq4FnMSq2YNx99Jt8SORpUAGx+qUPn5eTj7524A4FPYiahctXFwxkzfnahT3wyPHtzBqlnDcfVchNixSGSaby8hKjt/RUciI+1v1NSvjfc+6C12HCKq5ho0aQmP1buxZclMJF6OwcavpqDfxM/QY9gktXpETmlvWqmOF2HzjA9VqDMh/1zUrKmtLXIaIlIHenUMMXXZFnR2GQ5BEHBggw+2LZ0D+csXYkcjEbDxoQrz9FGK4jRzpz5DRU5DROpEU0sbwzwWY/CMr6Eh1UTs8QNY7TUKj/9+IHY0qmBsfKjCnD28G0JBAaze6wDjRk3EjkNEakYikaDrx5/Afdlm1DSogwc3r2LljCG4efGs2NGoArHxoQpRkJ+Pc38GAwC69OVFzUQknqbvd4SX/x40aNoKzzOeIHDueJzYuxWcz1c9sPGhCnEt5iSepD6Erp4B3ndwFjsOEam5OvXNMHPFDrTv2R8FBfnYt/Z7jB07FtnZ2WJHo3LGxocqxKvnctl9OIAzNRNRpaCto4tP5v6IAVPmQ0NDil9++QVdu3bFnTt3xI5G5YiND5W7zPRUXD1zHADQhXP3EFElIpFI4DhkHKYs3YR69eohNjYWdnZ2OHLkiNjRqJyw8aFydy50LwoK8tG4VTuYNG4mdhwioiKate2MmJgY2NnZIT09Hc7Ozli6dCmv+6mG2PhQuRIEAecOF17UrMpTkomIKpq5uTkiIyMxYcIEFBQUYP78+Rg8eDAyMvicr+qEjQ+Vq1vx0Uh7eBcy3Rpo062P2HGIiN5IR0cHGzduxLp166CtrY19+/ahQ4cOiI+PFzsalRE2PlSuzh76DQDQrvtHkOnWFDkNEVHJTJ48GSdPnoS5uTlu3LiBzp07Y/v27WLHojLAxofKzYtnmbgYeRgAZ2omoqqnQ4cOuHDhAnr37o3s7GyMHj0a06dPR05OjtjR6B2w8aFyc+H4AeTJc2Bi0QyNWrwvdhwiIpUZGRnh0KFD+PLLLwEAAQEBvOW9imPjQ+Xm7KHdAIBOLkPV6inIRFS9SKVSLFmyBAcPHkTdunVx/vx5tG/fHgcPHhQ7GpWCptgBqHp6cCsB929egVRTC7a9PhY7DhHRO+vbty8uXLiAYcOGITo6Gh999BHmz5+Put1GQypV7dfpyrDrpcrg2du6VK+jf/CMD5WLs38Wnu1574MPUcugrshpiIjKhoWFBSIjIzFjxgwAgLe3N9Z+PhZP0/4WORmVFBsfKnPynJeIObofANDJmRc1E1H1IpPJsHr1agQFBUFPTw+3L5/HiqkDce38SbGjUQmw8aEyF38qDC+eZaJOfTM0a28vdhwionIxfPhwxMTEoEGTlniW8RjrFk5CyOaVyM/PEzsavQEbHypz5/73MVdH5yHQ0OC3GBFVX82aNcOsVUGw/2gEBEHAkZ2BhR99PUoROxq9Bn8rUZlKT76HG3FnIJFI0MFpkNhxiIjKnZa2DENnfQu3+Ssgq1ETty+fh4/7AFyJOiZ2NCoGGx8qU9FhewEAzdrZo65xA5HTEBFVnHY9+sHLfy8aNmuN7Kyn2PjNVOxb+wPy5HKxo9G/sPGhMlNQUIDo0MLGp6PzYJHTEBFVvHoNLDBr5S50GzQWAHBi78/w8xiB1PuJIiejV9j4UJm5GXcGT1IfQreWPmzsPxQ7DhGRKDS1tTFw6gJM/HYtaujVxv2bV7Bi2mCcC90DQRDEjqf22PhQmTl7OBgA0L7HR9CW6YichohIXK279MSc/9uPJm06Qv4yG7t85mOb92d48TxL7GhqjY0PlYnsrAzEnwwFwI+5iIheqW1kjKlLt6DveE9oaEgRG34QP7kPQOKVGLGjqS02PlQmYo8fQF6uHKaW1mjYzEbsOERElYaGVIoPR7pjxortqGvcAI//foA1n43Gn1v9OOePCNj4UJk4F7oHQOHcPXwgKRFRUY1btcNngb/DttfHEAoKELrNH2u8PkF68j2xo6kVNj70zh7e/gv3rl+GhlQTtj35QFIiotfRramHT+b+iNHzf4JOTT3cTYiDj/vHvPC5ArHxoXd27nDh2Z7WXXqiVm0+kJSI6G3a9/gIcwJ/h5WNHXJeFF74vGXxTDx7+ljsaNWeptgBqGrLy5Uj5tirB5IOETkNEVHprQy7XqHHq2vcANN+3Ipjv27A4V9WI/5UGO5cjYWr1/do1al7pchYHfGMD72Tq2fD8TzjCfTr1kNzu65ixyEiqlIKL3yegtl+v8LYoimynqRhw1dT8Nuqr5Hz4rnY8aolNj70Tl59zGX74QBIpTyBSERUGg2btoLnmmDFjM9RB4Pg4z4At+PPi5ys+mHjQ6WW+fgR/oo+AQDo6MS5e4iI3oW2TAcDpy6A+7ItqF3PFOnJ9+A/ZzT2r1uGXHmO2PGqjVI1PgEBAbC0tISOjg5sbW0RGRn5xvqIiAjY2tpCR0cHVlZWCAwMLFITHByMVq1aQSaToVWrVti7d6/S9kWLFkEikSgtJiYmSjWCIGDRokUwMzODrq4uunfvjitXrpTmLVIJnD/yOwoK8mHRsi2MGzUROw4RUbVg3a4LPl/3Bzo6D4YgCAjfvQkrpg/GveuXxY5WLajc+AQFBcHDwwMLFy5EbGwsHBwc4OLigqSkpGLrExMT0bdvXzg4OCA2NhYLFizArFmzEBwcrKiJioqCq6sr3NzccPHiRbi5uWH48OE4e/as0r5at26N5ORkxRIfH6+0ffny5VixYgXWrFmD6OhomJiYoHfv3sjK4vTgZU0QBES/mruHZ3uIiMqUbk09jPjMGxO+DYBeHSP8ffcmVs0ajpDNK/m093ckEVScOKBTp05o37491q5dq1jXsmVLDBw4EN7e3kXq586di/379yMhIUGxzt3dHRcvXkRUVBQAwNXVFZmZmTh06JCipk+fPqhTpw527twJoPCMz759+xAXF1dsLkEQYGZmBg8PD8ydOxcAkJOTA2NjYyxbtgxTpkx563vLzMyEgYEBMjIyoK+v//bBqKZKctfA3YSLWDV7OLRkOli06yR0a+pVQDIiqoo8e1uX6nW8g6nQs4zH2LNmCeIiQgAAJo2tMXLODzC3fq/cj13aP7uKpsrvb5XO+MjlcsTExMDJyUlpvZOTE06fPl3sa6KioorUOzs74/z588jNzX1jzX/3eePGDZiZmcHS0hIjRozA7du3FdsSExORkpKitB+ZTAZHR8fXZqPSOxdaeMbuvQ96s+khIipHtQzqYszClRj7lR9qGdRFyp3rWDXLlWd/SkmlxictLQ35+fkwNjZWWm9sbIyUlJRiX5OSklJsfV5eHtLS0t5Y8+99durUCVu3bsXhw4exfv16pKSkwN7eHunp6Yp9vHpdSbPl5OQgMzNTaaG3k+e8RGx44b88OnLuHiKiCtHGwRlfbDiIto59UVCQjyM7A/HT9EG4kxAndrQqpVQXN//3WUyCILzx+UzF1f93/dv26eLigiFDhuC9997Dhx9+iIMHDwIAfv7551Jn8/b2hoGBgWIxNzd/7Xugf8SfCsPL51moY9wATdt0EjsOEZHaUDr7U9sQf9+9idUeI/B7oDdyXmSLHa9KUKnxMTIyglQqLXIGJTU1tciZlldMTEyKrdfU1IShoeEba163TwCoWbMm3nvvPdy4cUOxDwAq7Wf+/PnIyMhQLPfu8UFxJRH9v7l7OvQeCA0NzohARFTR2jg4Y+6Gg7D7cCAEQUDEni3wcf8YN2KjxI5W6an0W0tbWxu2trYICwtTWh8WFgZ7e/tiX9OlS5ci9aGhobCzs4OWltYba163T6DwY6qEhASYmpoCACwtLWFiYqK0H7lcjoiIiNfuRyaTQV9fX2mhN3uS+hA34gr/YnXoPUjkNERE6qumfh2M+mIZJn+3TjHvz9q547DrpwV4nvlU7HiVlsr/XPfy8sKGDRuwadMmJCQkwNPTE0lJSXB3dwdQeBZlzJgxinp3d3fcvXsXXl5eSEhIwKZNm7Bx40bMmTNHUTN79myEhoZi2bJl+Ouvv7Bs2TIcOXIEHh4eipo5c+YgIiICiYmJOHv2LIYOHYrMzEyMHVs4y6VEIoGHhwd++OEH7N27F5cvX8a4ceNQo0YNjBo1qrTjQ/8RHboXgiCgSZuOMDTlR4NERGJr2dERX6w7gA/6F/6uO3c4GMsm9UVcRAif+F4MlZ8x4OrqivT0dCxevBjJycmwsbFBSEgILCwsAADJyclKc/pYWloiJCQEnp6e8Pf3h5mZGfz8/DBkyD8Xxdrb22PXrl348ssv8dVXX6FJkyYICgpCp07/XD9y//59jBw5EmlpaahXrx46d+6MM2fOKI4LAF988QVevHiBadOm4cmTJ+jUqRNCQ0Ohp8e7jspCQUEBosMKJ5bs6MSLmomIKgudmrUwZOY3aN/zI/y68iv8nXQLW7/3RKsj+zFk5teoU99M7IiVhsrz+FRnnMen0Ovmzrh56RwC5rhBplsDi3adgky3RgUnI6KqiPP4VKw8uRxHg9bhyM5A5OflQlunBvqMmQmHQWNUfqZidZzHh0+VrAJK+5e/rL9hX13U3NaxL5seIioxNjAVS1NbG85uM9CmWx/8tuprJF6Owf51yxBzdD+Gzl4Mixbvix1RVLwlh0rkZfYzXDzxJwCgozMfUUFEVNmZWDTFdJ9tGO75HWro1caDWwnwmz0cwau/xYtn6jtvHRsfKpGLJ/6EPOcF6jVsjMat2osdh4iISkBDQwOdXYZh3sZDilvfT/2xA94T+uD8kX1qefEzGx8qkXOhr+buGfzGySqJiKjyqVW7LkZ9sQxTl/+M+uZWePY0HTuWz4X/HDek3LkhdrwKxcaH3urRgztIvBwDiYYG7HoPFDsOERGVUrO2nTEn8Hf0m/AZtGQ6uB0fDZ+pA7F/3TK8zH4mdrwKwcaH3io6tPAW9ubtP0Bto9fPpk1ERJWfppY2eo34FHM3hMDG/kMU5OchfPcmLJ3QBzFH91f7j7/Y+NAbFeTn/zN3Dy9qJiKqNuoaN8CERf6Y9N06GJlZIPPxI2xf9jn8PxuNB7f+EjteuWHjQ290PfY0MtL+hq6eAVp36SV2HCIiKmOt/jfzc9/xnoUff10+jxXTB2G33yKkpaWJHa/MsfGhNzr3v7l72vf4CFraMpHTEBFRedDU1saHI90xb+MhtOnWB0JBAU4f2IlmzZrBz88Pubm5YkcsM2x86LWyszJw+fQRAEBHJ37MRURU3dWpb4axX67CtB+3wtSqOZ4+fYrZs2ejbdu2CA0NFTtemWDjQ68Ve/wA8nLlMLW0RsNmrcWOQ0REFaRpm074zH8vAgMDYWhoiKtXr8LZ2Rn9+vVDQkKC2PHeCRsfeq2zh4MBFD6QlHP3EBGpFw2pFFOmTMGNGzfg4eEBTU1NhISE4L333sOMGTOq7PU/bHyoWA9u/YX7N65AqqkF214fix2HiIhEUqdOHaxcuRJXrlzBgAEDkJ+fD39/fzRt2hTLly/Hy5cvxY6oEjY+VKxzoYVne1p37olateuKnIaIiMRmbW2Nffv24ejRo2jTpg0yMjIwd+5cNG/eHNu2bUNBQYHYEUuEjQ8VkSeX48LR/QCAjn2GiJyGiIgqk549eyImJgZbtmxBgwYNkJSUBDc3N3To0AHHjh0TO95bsfGhIq6cOYbnmU+hb1gfzW0/EDsOERFVMlKpFGPHjsX169fx/fffQ09PDxcuXECvXr3g7OyM2NhYsSO+FhsfKuLVRc0deg+CVKopchoiIqqsatSogQULFuDWrVuYMWMGtLS0EBoaivbt22PUqFG4ffu22BGLYONDSu7fv49rMScBAB2d+TEXERG9Xb169bB69WokJCRg5MiRAICdO3eiRYsWmD59OpKTk0VO+A82PqRk69atEAoKYPVeB9RrYCF2HCIiqkKaNGmCHTt24MKFC3ByckJubi4CAgLQpEkTzJ07F+np6WJHZOND/xAEAZs2bQIAdOLZHiIiKqV27drh8OHDOH78OLp06YIXL15g+fLlsLKywpIlS/Ds2TPRsrHxIYXIyEjcunULMt0aeN/BWew4RERUxXXv3h2nTp3CgQMH0KZNG2RmZmLp0qV4/vy5aJl45SopvDrb07Z7P8h0a4ichoiIqgOJRIJ+/frBxcUFv/32Gx49egRjY2PR8rDxIQBAZmYmfvvtNwD8mIuIiMqehoYGXF1dxY7Bj7qo0M6dO5GdnY2WLVvComVbseMQERGVC57xIQDA+vXrAQCTJ08G+EBSIiKqpnjGhxAbG4uYmBhoa2vDzc1N7DhERETlho0PYcOGDQCAQYMGwcjISOQ0RERE5YeNj5rLzs7G9u3bAQCTJk0SOQ0REVH5YuOj5nbv3o2MjAxYWlqiZ8+eYschIiIqV2x81Nyrj7kmTpwIDQ1+OxARUfXG33Rq7K+//kJkZCQ0NDQwbtw4seMQERGVOzY+amzjxo0AgH79+qFBgwYipyEiIip/bHzUlFwux88//wyAFzUTEZH6YOOjpvbv349Hjx7B1NQUffv2FTsOERFRhWDjo6bWrVsHABg/fjw0NTmBNxERqQc2Pmroxo0bCAsLg0QiKXxEBRERkZpg46OGXp3tcXFxQePGjcUNQ0REVIHY+KiZly9fYvPmzQCAqVOnipyGiIioYrHxUTO7d+9Geno6GjVqBBcXF7HjEBERVahSNT4BAQGwtLSEjo4ObG1tERkZ+cb6iIgI2NraQkdHB1ZWVggMDCxSExwcjFatWkEmk6FVq1bYu3ev0nZvb2906NABenp6qF+/PgYOHIhr164p1YwbNw4SiURp6dy5c2neYrW1du1aAMCnn34KqVQqchoiIqKKpXLjExQUBA8PDyxcuBCxsbFwcHCAi4sLkpKSiq1PTExE37594eDggNjYWCxYsACzZs1CcHCwoiYqKgqurq5wc3PDxYsX4ebmhuHDh+Ps2bOKmoiICEyfPh1nzpxBWFgY8vLy4OTkhOfPnysdr0+fPkhOTlYsISEhqr7FauvSpUs4ffo0NDU1MWHCBLHjEBERVTiJIAiCKi/o1KkT2rdvrzhzAAAtW7bEwIED4e3tXaR+7ty52L9/PxISEhTr3N3dcfHiRURFRQEAXF1dkZmZiUOHDilq+vTpgzp16mDnzp3F5nj06BHq16+PiIgIdOvWDUDhGZ+nT59i3759qrwlhczMTBgYGCAjIwP6+vql2kd5WBl2vVSv8+xtrfT1tGnTsHbtWgwdOhS//fZbmR+PiIiql//+HqmsVPn9rdIZH7lcjpiYGDg5OSmtd3JywunTp4t9TVRUVJF6Z2dnnD9/Hrm5uW+sed0+ASAjIwMAULduXaX14eHhqF+/PqytrTF58mSkpqa+dh85OTnIzMxUWqqrrKws/PLLLwB4UTMREakvlRqftLQ05Ofnw9jYWGm9sbExUlJSin1NSkpKsfV5eXlIS0t7Y83r9ikIAry8vNC1a1fY2Ngo1ru4uGD79u04duwYfvrpJ0RHR6Nnz57Iyckpdj/e3t4wMDBQLObm5m8egCps586dePbsGaytrdGjRw+x4xAREYmiVFP2SiQSpa8FQSiy7m31/12vyj5nzJiBS5cu4eTJk0rrXV1dFf9vY2MDOzs7WFhY4ODBgxg8eHCR/cyfPx9eXl6KrzMzM6tl8yMIguKjSXd39zf+WREREVVnKjU+RkZGkEqlRc7EpKamFjlj84qJiUmx9ZqamjA0NHxjTXH7nDlzJvbv348TJ06gYcOGb8xramoKCwsL3Lhxo9jtMpkMMpnsjfuoDqKiohAXFwcdHR2MHTtW7DhERESiUemjLm1tbdja2iIsLExpfVhYGOzt7Yt9TZcuXYrUh4aGws7ODlpaWm+s+fc+BUHAjBkzsGfPHhw7dgyWlpZvzZueno579+7B1NS0RO+vuvLz8wMAjBo1qsg1UUREROpE5dvZvby8sGHDBmzatAkJCQnw9PREUlIS3N3dARR+fDRmzBhFvbu7O+7evQsvLy8kJCRg06ZN2LhxI+bMmaOomT17NkJDQ7Fs2TL89ddfWLZsGY4cOQIPDw9FzfTp07Ft2zbs2LEDenp6SElJQUpKCl68eAEAePbsGebMmYOoqCjcuXMH4eHh6N+/P4yMjDBo0KDSjk+Vd//+fezevRsAMGvWLJHTEBERiUvla3xcXV2Rnp6OxYsXIzk5GTY2NggJCYGFhQUAIDk5WWlOH0tLS4SEhMDT0xP+/v4wMzODn58fhgwZoqixt7fHrl278OWXX+Krr75CkyZNEBQUhE6dOilqXl2j0r17d6U8mzdvxrhx4yCVShEfH4+tW7fi6dOnMDU1RY8ePRAUFAQ9PT1V32a1ERAQgPz8fDg6OqJNmzZixyEiIhKVyvP4VGfVbR4f967mMDc3R3p6Ovbs2VPiM1+cx4eIiADO40NVzI4dO5Ceng4LCwt8/PHHYschIiISHRufakoQBKxatQpA4e3/fC4XERERG59q6+bFs4iPj0eNGjUwceJEseMQERFVCmx8qqnIfYWPpxgzZgzq1KkjchoiIqLKgY1PNZSefA9Xoo4CKJzwkYiIiAqx8amGTu7fDkEQ0Lt3b7Rq1UrsOERERJUGG59qJjsrA2dCggBAaQJIIiIiYuNT7Zzavx05L7JhamkNFxcXseMQERFVKmx8qhH5yxc4sW8rAKCn66d8CjsREdF/sPGpRs4eDsbzjCeoa9wAbR15toeIiOi/2PhUE/l5uQjfvQkA0GPYREilKj+GjYiIqNpj41NNxIaH4MnfD1CrtiE6OA95+wuIiIjUEE8LVAMFBQU4FrQeANBt0Fhoy3QA8GGjRERE/8UzPtVAwtlwpNy9AVmNmvig/0ix4xAREVVabHyqOEEQcDRoHQDgg/6joFtLX+RERERElRcbnyru5sWzuHM1Fppa2ug2cIzYcYiIiCo1Nj5VmCAIOLhpBQCgk8sw6BvWFzkRERFR5cbGpwq7EnUUSX9dhLZMF71HTRU7DhERUaXHxqeKKsjPR8hmXwCAw6Ax0K9bT9xAREREVQAbnyrqwvE/kHL3BnRr6aPHsIlixyEiIqoS2PhUQXm5cvz5sx8AoKfrZNTQMxA5ERERUdXAxqcKigr5FY//fgD9uvXgMMBN7DhERERVBhufKibnRTaO7FgLAOj9yTRo6+iKnIiIiKjqYONTxUTu24qsJ2kwNDVHpz5DxY5DRERUpbDxqUKepD7EkZ3/BwDoM2YWNLW0RU5ERERUtbDxqSIEQcCeNYshf5kNSxtbtOvxkdiRiIiIqhw2PlVE/KkjuHLmOKSaWhg2ezE0NPhHR0REpCr+9qwCXj5/hr0BSwAAPYZNhIlFU5ETERERVU1sfKqAQz/7IiPtbxiaNcKHfDQFERFRqbHxqeSio6Nx8vdtAIChMxdBW6YjciIiIqKqi41PJZaXl4cpU6ZAEAS079kfzW0/EDsSERFRlcbGpxL7+uuvERsbC109AwyYMl/sOERERFUeG59K6tdff4W3tzcAYOiMb6BXx1DkRERERFWfptgBqKiLFy9i/PjxAIDPP/8cpj36iZyIiIioeuAZn0omLS0NAwYMQHZ2NpycnBRnfYiIiOjdsfGpRHJzczF8+HDcvXsXTZo0wa5duyCVSsWORUREVG2w8akk8vPzMWPGDBw/fhy1atXC77//jjp16ogdi4iIqFrhNT6VwJMnTzBy5EgcPnwYALB161a0bt1a5FRERETVT6nO+AQEBMDS0hI6OjqwtbVFZGTkG+sjIiJga2sLHR0dWFlZITAwsEhNcHAwWrVqBZlMhlatWmHv3r0qH1cQBCxatAhmZmbQ1dVF9+7dceXKldK8xQpz+fJldOjQAYcPH4auri527tyJQYMGiR2LiIioWlK58QkKCoKHhwcWLlyI2NhYODg4wMXFBUlJScXWJyYmom/fvnBwcEBsbCwWLFiAWbNmITg4WFETFRUFV1dXuLm54eLFi3Bzc8Pw4cNx9uxZlY67fPlyrFixAmvWrEF0dDRMTEzQu3dvZGVlqfo2K8SePXvQuXNn3Lp1CxYWFjh9+jRGjBghdiwiIqJqSyIIgqDKCzp16oT27dtj7dq1inUtW7bEwIEDi70Dae7cudi/fz8SEhIU69zd3XHx4kVERUUBAFxdXZGZmYlDhw4pavr06YM6depg586dJTquIAgwMzODh4cH5s6dCwDIycmBsbExli1bhilTprz1vWVmZsLAwAAZGRnQ19dXZVhKJC8vD2fOnEFISAhCQkJw8eJFAEDPnj0RFBQEIyOjYl+3Mux6mWchIiJ6G8/e1mJHKBFVfn+rdI2PXC5HTEwM5s2bp7TeyckJp0+fLvY1UVFRcHJyUlrn7OyMjRs3Ijc3F1paWoiKioKnp2eRGl9f3xIfNzExESkpKUrHkslkcHR0xOnTp4ttfHJycpCTk6P4OiMjA0DhAJal5ORkLFiwAEePHlUc45Xp06dj8eLF0NTUfO1xXz5/VqZ5iIiISqKsfx+Wl1c5S3IuR6XGJy0tDfn5+TA2NlZab2xsjJSUlGJfk5KSUmx9Xl4e0tLSYGpq+tqaV/ssyXFf/be4mrt37xabzdvbG99++22R9ebm5sXWlwd/f3/4+/tX2PGIiIhKaoHYAVSUlZUFAwODN9aU6q4uiUSi9LUgCEXWva3+v+tLss+yqnll/vz58PLyUnxdUFCAx48fw9DQ8I3vpzQyMzNhbm6Oe/fulcvHaOqK41o+OK7lg+NaPjiu5aMqjasgCMjKyoKZmdlba1VqfIyMjCCVSouc3UlNTS1ypuUVExOTYus1NTVhaGj4xppX+yzJcU1MTAAUnvkxNTUtUTaZTAaZTKa0rnbt2sXWlhV9ff1K/w1UFXFcywfHtXxwXMsHx7V8VJVxfduZnldUuqtLW1sbtra2CAsLU1ofFhYGe3v7Yl/TpUuXIvWhoaGws7ODlpbWG2te7bMkx7W0tISJiYlSjVwuR0RExGuzERERkZoRVLRr1y5BS0tL2Lhxo3D16lXBw8NDqFmzpnDnzh1BEARh3rx5gpubm6L+9u3bQo0aNQRPT0/h6tWrwsaNGwUtLS1h9+7dippTp04JUqlUWLp0qZCQkCAsXbpU0NTUFM6cOVPi4wqCICxdulQwMDAQ9uzZI8THxwsjR44UTE1NhczMTFXfZpnLyMgQAAgZGRliR6lWOK7lg+NaPjiu5YPjWj6q67iq3PgIgiD4+/sLFhYWgra2ttC+fXshIiJCsW3s2LGCo6OjUn14eLjQrl07QVtbW2jcuLGwdu3aIvv87bffhObNmwtaWlpCixYthODgYJWOKwiCUFBQIHzzzTeCiYmJIJPJhG7dugnx8fGleYtl7uXLl8I333wjvHz5Uuwo1QrHtXxwXMsHx7V8cFzLR3UdV5Xn8SEiIiKqqviQUiIiIlIbbHyIiIhIbbDxISIiIrXBxoeIiIjUBhufChAQEABLS0vo6OjA1tYWkZGRYkeqUry9vdGhQwfo6emhfv36GDhwIK5du6ZUIwgCFi1aBDMzM+jq6qJ79+64cuWKSImrJm9vb0gkEnh4eCjWcVxL58GDBxg9ejQMDQ1Ro0YNtG3bFjExMYrtHFfV5eXl4csvv4SlpSV0dXVhZWWFxYsXo6CgQFHDcS2ZEydOoH///jAzM4NEIsG+ffuUtpdkHHNycjBz5kwYGRmhZs2a+Pjjj3H//v0KfBfvQMQ7ytTCq/mH1q9fL1y9elWYPXu2ULNmTeHu3btiR6synJ2dhc2bNwuXL18W4uLihH79+gmNGjUSnj17pqhZunSpoKenJwQHBwvx8fGCq6trpZnDqSo4d+6c0LhxY+H9998XZs+erVjPcVXd48ePBQsLC2HcuHHC2bNnhcTEROHIkSPCzZs3FTUcV9V99913gqGhoXDgwAEhMTFR+O2334RatWoJvr6+ihqOa8mEhIQICxcuFIKDgwUAwt69e5W2l2Qc3d3dhQYNGghhYWHChQsXhB49eght2rQR8vLyKvjdqI6NTznr2LGj4O7urrSuRYsWwrx580RKVPWlpqYKABTzOBUUFAgmJibC0qVLFTUvX74UDAwMhMDAQLFiVhlZWVlCs2bNhLCwMMHR0VHR+HBcS2fu3LlC165dX7ud41o6/fr1EyZMmKC0bvDgwcLo0aMFQeC4ltZ/G5+SjOPTp08FLS0tYdeuXYqaBw8eCBoaGsKff/5ZYdlLix91lSO5XI6YmBg4OTkprXdycsLp06dFSlX1ZWRkAADq1q0LAEhMTERKSorSOMtkMjg6OnKcS2D69Ono168fPvzwQ6X1HNfS2b9/P+zs7DBs2DDUr18f7dq1w/r16xXbOa6l07VrVxw9ehTXr18HAFy8eBEnT55E3759AXBcy0pJxjEmJga5ublKNWZmZrCxsakSY12qp7NTyaSlpSE/P7/IQ1KNjY2LPHCVSkYQBHh5eaFr166wsbEBAMVYFjfOd+/erfCMVcmuXbtw4cIFREdHF9nGcS2d27dvY+3atfDy8sKCBQtw7tw5zJo1CzKZDGPGjOG4ltLcuXORkZGBFi1aQCqVIj8/H99//z1GjhwJgN+vZaUk45iSkgJtbW3UqVOnSE1V+N3GxqcCSCQSpa8FQSiyjkpmxowZuHTpEk6ePFlkG8dZNffu3cPs2bMRGhoKHR2d19ZxXFVTUFAAOzs7/PDDDwCAdu3a4cqVK1i7di3GjBmjqOO4qiYoKAjbtm3Djh070Lp1a8TFxcHDwwNmZmYYO3asoo7jWjZKM45VZaz5UVc5MjIyglQqLdIBp6amFumm6e1mzpyJ/fv34/jx42jYsKFivYmJCQBwnFUUExOD1NRU2NraQlNTE5qamoiIiICfnx80NTUVY8dxVY2pqSlatWqltK5ly5ZISkoCwO/X0vr8888xb948jBgxAu+99x7c3Nzg6ekJb29vABzXslKScTQxMYFcLseTJ09eW1OZsfEpR9ra2rC1tUVYWJjS+rCwMNjb24uUquoRBAEzZszAnj17cOzYMVhaWiptt7S0hImJidI4y+VyREREcJzfoFevXoiPj0dcXJxisbOzwyeffIK4uDhYWVlxXEvhgw8+KDLdwvXr12FhYQGA36+llZ2dDQ0N5V9ZUqlUcTs7x7VslGQcbW1toaWlpVSTnJyMy5cvV42xFu2yajXx6nb2jRs3ClevXhU8PDyEmjVrCnfu3BE7WpUxdepUwcDAQAgPDxeSk5MVS3Z2tqJm6dKlgoGBgbBnzx4hPj5eGDlyJG9jLYV/39UlCBzX0jh37pygqakpfP/998KNGzeE7du3CzVq1BC2bdumqOG4qm7s2LFCgwYNFLez79mzRzAyMhK++OILRQ3HtWSysrKE2NhYITY2VgAgrFixQoiNjVVMs1KScXR3dxcaNmwoHDlyRLhw4YLQs2dP3s5O//D39xcsLCwEbW1toX379orbsKlkABS7bN68WVFTUFAgfPPNN4KJiYkgk8mEbt26CfHx8eKFrqL+2/hwXEvnjz/+EGxsbASZTCa0aNFCWLdundJ2jqvqMjMzhdmzZwuNGjUSdHR0BCsrK2HhwoVCTk6OoobjWjLHjx8v9mfq2LFjBUEo2Ti+ePFCmDFjhlC3bl1BV1dX+Oijj4SkpCQR3o3qJIIgCOKcayIiIiKqWLzGh4iIiNQGGx8iIiJSG2x8iIiISG2w8SEiIiK1wcaHiIiI1AYbHyIiIlIbbHyIiIhIbbDxISIiIrXBxoeIiIjUBhsfIiIiUhtsfIiIiEhtsPEhIiIitfH/L4Tg2yOXlFoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Sample \n",
    "Nt = 400\n",
    "R_lim = 100\n",
    "print(len(rtrue))\n",
    "print(Lbox)\n",
    "\n",
    "# Get relative probability of each object according to selection function\n",
    "# lam = float(config[f'sample_{i}']['lam'])\n",
    "# R_max = float(config['mock']['R_max'])\n",
    "prob = np.exp(- lam * rtrue / R_max)\n",
    "prob /= np.sum(prob)\n",
    "\n",
    "accepted_count = 0\n",
    "rsamp = np.zeros(Nt)\n",
    "\n",
    "# Loop until we have Nt valid positions\n",
    "while accepted_count < Nt:\n",
    "\n",
    "    ids = np.random.choice(len(rtrue), size=Nt, p=prob, replace=False)\n",
    "    \n",
    "    # Here I would actually scatter\n",
    "    \n",
    "    valid_indices = rtrue[ids] < R_lim\n",
    "    ids = ids[valid_indices]\n",
    "    \n",
    "    # Calculate how many valid positions we need to reach Nt\n",
    "    remaining_needed = Nt - accepted_count\n",
    "    selected_count = min(len(ids), remaining_needed)\n",
    "    ids = ids[:selected_count]\n",
    "\n",
    "    # Append only the needed number of valid positions\n",
    "    rsamp[accepted_count:accepted_count+selected_count] = rtrue[ids]\n",
    "\n",
    "    # Update the accepted count\n",
    "    accepted_count += selected_count\n",
    "\n",
    "    print(f'\\tMade {accepted_count} of {Nt}')\n",
    "    \n",
    "    # Set up for next iteration\n",
    "    prob[ids] = 0\n",
    "    prob /= np.sum(prob)\n",
    "\n",
    "plt.hist(rsamp, bins=30, density=True, alpha=0.5)\n",
    "x = plt.gca().get_xlim()\n",
    "x = [max(0, x[0]), x[1]]\n",
    "x = np.linspace(x[0], x[1], 100)\n",
    "y = x ** 2 * np.exp( - lam * x / R_max) * (lam / R_max) ** 3 / 2\n",
    "plt.plot(x, y, color='k')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "73491c07-4ad8-4aa7-b277-192bc28f26d8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3, 12449803)\n"
     ]
    }
   ],
   "source": [
    "comb_x = np.array([xtrue, ytrue, ztrue])\n",
    "print(comb_x.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "7302fe50-8271-4d0f-bcf2-dc02bdae4f94",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['FILE_FORMAT', 'PARTICLE_MASS', 'MASS_DEFINITION', 'MASS_DEFINITION2', 'MASS_DEFINITION3', 'MASS_DEFINITION4', 'MASS_DEFINITION5', 'STRICT_SO_MASSES', 'MIN_HALO_OUTPUT_SIZE', 'FORCE_RES', 'FORCE_RES_PHYS_MAX', 'SCALE_NOW', 'h0', 'Ol', 'Om', 'W0', 'WA', 'GADGET_ID_BYTES', 'GADGET_MASS_CONVERSION', 'GADGET_LENGTH_CONVERSION', 'GADGET_SKIP_NON_HALO_PARTICLES', 'GADGET_HALO_PARTICLE_TYPE', 'RESCALE_PARTICLE_MASS', 'TIPSY_LENGTH_CONVERSION', 'TIPSY_VELOCITY_CONVERSION', 'PARALLEL_IO', 'PARALLEL_IO_SERVER_ADDRESS', 'PARALLEL_IO_SERVER_PORT', 'PARALLEL_IO_WRITER_PORT', 'PARALLEL_IO_SERVER_INTERFACE', 'PARALLEL_IO_CATALOGS', 'RUN_ON_SUCCESS', 'RUN_PARALLEL_ON_SUCCESS', 'LOAD_BALANCE_SCRIPT', 'INBASE', 'FILENAME', 'STARTING_SNAP', 'RESTART_SNAP', 'NUM_SNAPS', 'NUM_BLOCKS', 'NUM_READERS', 'PRELOAD_PARTICLES', 'SNAPSHOT_NAMES', 'LIGHTCONE_ALT_SNAPS', 'BLOCK_NAMES', 'OUTBASE', 'OVERLAP_LENGTH', 'NUM_WRITERS', 'FORK_READERS_FROM_WRITERS', 'FORK_PROCESSORS_PER_MACHINE', 'OUTPUT_FORMAT', 'DELETE_BINARY_OUTPUT_AFTER_FINISHED', 'FULL_PARTICLE_CHUNKS', 'BGC2_SNAPNAMES', 'SHAPE_ITERATIONS', 'WEIGHTED_SHAPES', 'BOUND_PROPS', 'BOUND_OUT_TO_HALO_EDGE', 'DO_MERGER_TREE_ONLY', 'IGNORE_PARTICLE_IDS', 'EXACT_LL_CALC', 'TRIM_OVERLAP', 'ROUND_AFTER_TRIM', 'LIGHTCONE', 'PERIODIC', 'LIGHTCONE_ORIGIN', 'LIGHTCONE_ALT_ORIGIN', 'LIMIT_CENTER', 'LIMIT_RADIUS', 'SWAP_ENDIANNESS', 'GADGET_VARIANT', 'ART_VARIANT', 'FOF_FRACTION', 'FOF_LINKING_LENGTH', 'INITIAL_METRIC_SCALING', 'INCLUDE_HOST_POTENTIAL_RATIO', 'TEMPORAL_HALO_FINDING', 'MIN_HALO_PARTICLES', 'UNBOUND_THRESHOLD', 'ALT_NFW_METRIC', 'EXTRA_PROFILING', 'TOTAL_PARTICLES', 'BOX_SIZE', 'OUTPUT_LEVELS', 'DUMP_PARTICLES', 'ROCKSTAR_CONFIG_FILENAME', 'AVG_PARTICLE_SPACING', 'SINGLE_SNAP', 'RAMSES_DP'])\n",
      "0.315\n",
      "<class 'float'>\n"
     ]
    }
   ],
   "source": [
    "dirname = '/data54/lavaux/VELMASS/halo_central/halos_new/'\n",
    "\n",
    "def parse_file_to_dict(file_path):\n",
    "    config_dict = {}\n",
    "    with open(file_path, 'r') as file:\n",
    "        for line in file:\n",
    "            line = line.strip()\n",
    "            if not line or line.startswith(\"#\"):  # Skip empty lines and comments\n",
    "                continue\n",
    "            key, value = line.split(\"=\", 1)\n",
    "            key = key.strip()\n",
    "            value = value.strip()\n",
    "\n",
    "            # Convert the value to the appropriate type (int, float, or leave as string)\n",
    "            if value.isdigit():\n",
    "                value = int(value)\n",
    "            else:\n",
    "                try:\n",
    "                    value = float(value)\n",
    "                except ValueError:\n",
    "                    pass  # Leave it as a string if it cannot be converted\n",
    "\n",
    "            config_dict[key] = value\n",
    "\n",
    "    return config_dict\n",
    "\n",
    "\n",
    "config_dict = parse_file_to_dict(f'{dirname}/rockstar.cfg')\n",
    "print(config_dict.keys())\n",
    "print(config_dict['Om'])\n",
    "print(type(config_dict['Om']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3480b9de-7dc7-4c1c-a5c3-9d3d9c8e45ca",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "borg_new",
   "language": "python",
   "name": "borg_new"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}