mgxs-part-i.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multigroup Cross Section Generation Part I: Introduction\n",
    "This IPython Notebook introduces the use of the `openmc.mgxs` module to calculate multi-group cross sections for an infinite homogeneous medium. In particular, this Notebook introduces the the following features:\n",
    "\n",
    "* **General equations** for scalar-flux averaged multi-group cross sections\n",
    "* Creation of multi-group cross sections for an **infinite homogeneous medium**\n",
    "* Use of **tally arithmetic** to manipulate multi-group cross sections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to Multi-Group Cross Sections (MGXS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many Monte Carlo particle transport codes, including OpenMC, use continuous-energy nuclear cross section data. However, most deterministic neutron transport codes use *multi-group cross sections* defined over discretized energy bins or *energy groups*. An example of U-235's continuous-energy fission cross section along with a 16-group cross section computed for a light water reactor spectrum is displayed below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAAI8CAYAAAAz5idmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlcFdX/P/DXXJTlsigqpeICon4k9z1MwyVJJVHcEFMTl8o1rMT8lAp+3D/lj9K0VNyKQDEXzCWzML9oH0tRK5PS3HPFBQVEtvP7g+6N693mLnNnDvf9fDx4lDNnzrxmHOR9hzNnBMYYAyGEEEIIIcRiKrkDEEIIIYQQwisqpgkhhBBCCLESFdOEEEIIIYRYiYppQgghhBBCrETFNCGEEEIIIVaiYpoQQgghhBArUTFNCCGEEEKIlaiYJoQQQgghxEpUTBNCCCGEEGIlKqYJIQTAhg0boFKpsHHjRrmjcI/OJSHEmVAxTQinDh48CJVKhR49ehhtc/HiRahUKgQGBoru96+//sLy5cvRt29fBAYGwt3dHbVq1UJYWBi2b99ucJtffvkF48ePR9u2beHn5wc3NzfUr18fPXv2xBdffIGysjKD250+fRoxMTFo2rQpPD09Ubt2bXTp0gVr1qxBcXGx6Mzx8fFQqVRQqVQYP3680XZbtmzRtuvWrZvOOkEQtF/ENo44l2VlZdi6dSsGDx6M+vXrw8PDA15eXnjmmWfw2muv4ciRI5Lt2xieP0QkJSVBpVKhX79+RtuEh4dDpVJh7dq1Osvv37+POXPmoE2bNvDy8oK7uzvq1auHkJAQvP322zh58qTU8QmRVRW5AxBCbCOmYLGkqFm+fDmWLl2KwMBA9OzZE7Vr18bFixexbds2HDhwALGxsVi2bJnONllZWdi5cydCQkLQtWtXVKtWDdevX8euXbswcuRIbNmyBTt27NDZ5ttvv0W/fv1QWlqKvn37YvDgwcjNzcWuXbvw2muvYevWrdi3b59F2atUqYItW7bgww8/hKenp976NWvWoEqVKigpKdHrNzIyEiEhIahdu7bo/RHDpD6XN27cwJAhQ3DkyBH4+Pigd+/eCAoKAmMMZ8+exebNm7FmzRosX74ckydPliSDKTx+IBs3bhx27dqF9PR0rFy5EpMmTdJZv2rVKuzduxf9+/fX+cB67do1PPfcc7h06RKCgoIwatQo1KpVC/fu3cPRo0exbNkyqNVqtGnTxtGHRIjjMEIIlzIyMpggCKxHjx5G21y4cIEJgsACAwNF97tt2zZ28OBBveVnzpxh1apVY4IgsGPHjumse/z4scG+Hjx4wIKDg5kgCOz777/XWRcaGsoEQWCbNm3SWZ6fn8+aN29ucBtj5s6dywRBYAMGDGCCILC1a9fqtblw4QJTqVRs4MCBTBAE1q1bN1F9E2XJz89nrVu3ZoIgsBEjRrD79+/rtcnLy2MJCQls4cKFDs22fv16JggC27Bhg0P3ay+3bt1iTz31FPP09GS///67dvnvv//O1Go1e/rpp9mtW7d0thk3bhwTBIGNGzfOYJ8XLlzQ+/eCkMqGhnkQQnRERkYiNDRUb3mzZs0QFRUFAPj+++911rm6uhrsy9vbGy+++CKA8ruJFeXk5EAQBEREROgsV6vV6NmzJwDgzp07FmXv168f6tati6SkJL11SUlJYIwZHQZi6lf0V69exbRp09CkSROo1WrUrFkTnTt3xvz583XaBQQEIDAwEA8ePMAbb7yBhg0bwtXVFQkJCdo233zzDV588UXUqFED7u7uaNq0Kd555x3k5ubq7ffcuXMYP348goKC4OHhgRo1auCZZ57B66+/jrt37+q137x5M3r16oUaNWrAw8MDgYGBGDFiBI4fP67TrrCwEIsWLULLli3h6emJatWq4fnnn8fmzZv1+tQMFYqJiUF2djYGDhyIGjVqwMvLC926dcM333wjybk0ZtmyZfj555/RtWtXJCcno1q1anptPD09MWfOHLz11lvaZdeuXcO8efPw3HPPoXbt2nBzc4O/vz9GjBiB3377zebj7t69O8aOHQsAiImJ0Q4nUqlUuHz5MgBgzJgxOn+uSDNsq+K1oulXpVKhqKgIc+bMQZMmTeDm5oaYmBidczplyhQ0atRIOyxrwIABOHbsmKhzquHn54c1a9agoKAAI0eORGlpKUpKSjBy5EgUFhZi9erV8PPz09nmyJEjEAQB06ZNM9hnQEAA2rdvb1EOQnhDwzwIIaJVrVpV57/mFBQU4LvvvoO7uztCQkJ01r3wwgv47bffsHPnTowePVq7PD8/H99++y28vLzQpUsXi/K5uLhgzJgxWLhwIX777Tc888wzAIDS0lKsX78ezz77LFq0aGGyjyd/RX/s2DG8+OKLuHfvHrp3744hQ4YgPz8fp0+fRkJCAt577z2dbR8/fowePXogNzcXffv2hZeXl3bM+sqVKzFlyhR4e3tj2LBh8PPzw3fffYelS5ciPT0dR44cQfXq1QGUF3+dOnVCXl4ewsPDMWzYMBQWFuL8+fNITk7GtGnTUKNGDQAAYwwxMTHYtGkT/Pz8MGTIEPj5+eHy5cs4ePAgmjVrpi1oioqKEBYWhszMTDRv3hxTpkxBfn4+0tLSEB0djRMnTmDx4sV65+XChQvo0qULWrVqhYkTJ+LatWvYvHkz+vbtiy+++ALDhg2z67k0Zs2aNQCA2bNnm21b8UPeoUOHsGTJEvTs2RPt2rWDp6cn/vjjD2zduhXp6ek4fPgwWrdubfVxx8TEwNfXFzt37sTAgQN1hjVULPjNDQExtn7QoEE4fvw4+vXrh1q1auHpp58GUD7EKiwsDPfu3UPfvn0xZMgQ3L59Gzt27EDXrl2xfft29O3b1+y50oiIiMDYsWOxbt06zJs3D4wxHDt2DGPHjtX74AuUF+DZ2dn4/fff0apVK9H7IaRSkfnOOCHESlIN8zAmNzeXPf3008zFxYVlZ2cbbHP27Fk2d+5c9t5777EJEyawunXrsjp16rAdO3botc3Ly2PR0dGsatWq7KWXXmJxcXFs4sSJzN/fn9WrV4/t27dPdDbNMI+kpCR2/vx5plKp2Jtvvqldv2vXLu16zTl5cpiH5lf0Gzdu1C57/PgxCwgIYCqViqWmpurt9+rVqzp/btiwIRMEgfXu3ZsVFBTorLtw4QKrWrUqq169Ojt79qzOutdff50JgsAmTJigXfbhhx8yQRDYhx9+qLffgoIC9ujRI+2fP/30UyYIAnv22WfZgwcPdNqWlpay69eva/+8YMECJggCi4iIYKWlpdrlN2/e1ObPzMzUyS0IAhMEgcXFxen0fezYMVa1alXm6+urs197nEtDLl26xARBYK6urkaHFhlz69YtlpeXp7c8KyuLeXp6sj59+ugst9dxV/TKK68wQRDYpUuX9NZpvp8TEhJ0lmuGQ7Vu3ZrduXNHZ11xcTELCgpiarWaHT58WGfdtWvXmL+/P6tduzYrLCw0mMeYhw8fskaNGrEqVaqwKlWqsKCgIIPnjjHGVq1axQRBYD4+PmzGjBls3759ekNBCKnsqJgmhFOOLKbLysrY0KFDmSAIbPLkyUbb7d27V1uACILAqlSpwqZNm8Zu3LhhtH2HDh10tnFzc2MzZ85k9+7dE52vYjHNGGO9evVifn5+rLi4mDHG2IABA5iPjw/Lz8+3qJjeunUrEwSBDRw4UFSOhg0bMpVKxU6dOqW37j//+Q8TBIHNnj1bb93du3eZt7c3U6vVrKioiDHG2EcffcQEQWCrV682u98WLVowlUrFTp48abZtUFAQc3Fx0SvoGWNs7dq1TBAENnbsWO0yzfny9fU1WFCNGTNG77zZ41wacvToUSYIAqtTp47VfRjy0ksvMXd3d1ZSUqJdZq/jrsiWYnrnzp162+zYsYMJgsBmzpxpcH+JiYlMEAS2e/duwwduguZYVCoV+/rrr022nT17NlOr1TrfxwEBAey1115jv/76q8X7JoQ3NGaaECdz8uRJxMfH63x9+OGHJreZPn06tm7diq5du+rN5FFRnz59UFZWhuLiYpw7dw5z5szBJ598gg4dOuDWrVs6bT/++GP069cPKpUKmZmZyMvLw5UrV5CQkIAPPvgAnTt3NjiOWIzx48cjJycHO3bswPXr17F7925ERUVBrVZb1M///vc/ALDo1+Tu7u4Gf9194sQJADA4laGvry/atm2LR48e4cyZMwCAAQMGwMvLC5MnT8awYcOwevVqg2N7NcMknn76aYPDFCp6+PAhzp8/D39/fzRu3Fhvfa9evXSyVqQZGvEkzfh6c9OfWXMu7W337t3o378/6tSpA1dXV+2Y5t27d6OoqAg5OTl629h63PYgCAI6d+6st/yHH34AUD4U5cnv6fj4ePz4448AgOzsbIv29+jRIyxZsgRA+RCitLQ0k+3nzZuH69evIzU1FdOnT0doaChu3ryJ1atXo23btli3bp1F+yeENzRmmhBOqVTln4WNzeFccZ2mLQCcOnUK8+bN02kXEBCAN954w2Afb775Jj766COEhoZi9+7dRh82rMjFxQWNGjXC7Nmz4ebmhnfeeQfvv/8+li5dCgDIy8vDzJkzoVarsWvXLjz11FMAyh8+nDlzJm7evInExEQkJiZi7ty5Zvf3pMjISNSoUQNr167F2bNnUVpaanL+aWPu378PAPD39xe9jeZYnqT5YGBsurg6derotGvQoAF+/PFHxMfHY9++fdi6dSsAoH79+oiLi9NO+WZJRnMZNMsNfYjRjNG1ZJuKrDmXT6pbty6A8gdTHz9+DDc3N9Hbfvjhh5g+fTpq1KiB3r17o0GDBlCr1RAEAdu3b8epU6fw+PFjve1sPW57MZRD84CuqWJXEATk5+dbtK+4uDj8/vvviI2NxcGDB5GUlITIyEiTc1D7+Phg2LBh2jHkBQUFWLx4MebPn4/JkyfjpZdeMvq9QQjv6M40IZzSPNRkasYLzZ02zUNtAPDKK6+grKxM5+v8+fMGt582bRoSExPRs2dP7N271+I7uwC0s3n88ssv2mXZ2dkoKChAcHCwwR+w3bt3BwC9WSjEcnNzw8svv4wDBw5gxYoVaNGiBTp16mRxP5rzdvXqVdHbGHuATPP3df36dYPrNcsrPqzWrFkzpKam4s6dOzh27BgWL16MsrIyTJ06FRs2bNDJ+Ndff5nNpun7yZlVTGXQuHnzpsFtNH0Z2qYia87lk+rVq4cGDRqguLgYhw4dEr1dSUkJ4uPjUadOHZw+fRopKSlYsmQJ5s6dizlz5pgs8mw97oo0H2pLSkr01mk+bFhCs+/09HS972nNV2lpqaiHNTX279+Pjz/+GK1atcKSJUvw2Wefwc3NDePHjzc4g4wxarVaO3vK48ePcfjwYYuPjxBeUDFNCKeaNWsGV1dX/PHHH0Z/yGl+DWzpU/aMMUycOBErVqxAWFgYdu/eDXd3d6tyaoo8Hx8f7TLNHUVDv1YHgNu3b+u0s8b48eNRVlaGGzduYNy4cVb1oZmB5Ouvv7Y6h0a7du0AlE+B9qT79+/j5MmT8PDwQHBwsN56FxcXtGvXDnFxcUhJSQEA7UtwPD090aJFC9y4cQOnTp0ymcHb2xtBQUG4evUqzp07p7c+IyNDJ2tFWVlZyMvL01uuOZ62bdua3Le9zuWrr74KAJg/fz4YYybbFhUVASi/znJzc9GlSxe9O7x5eXnIysoy+iHIkuN2cXEBUD57jCG+vr4AYHBqPEunsQP+OaeWfLAw5e7du4iJiYGbmxs+//xzVK1aFc2bN8d//vMf3LhxAxMnTrS4T29vb7tkI0TRZB6zrThHjhxhgiCw+fPnyx2FELNGjx6t98CYxpUrV5i/vz9TqVQGX8JiTFlZGRs/fjwTBIGFh4eLmjXhp59+Mrj81q1brGXLlkwQBJaSkqJdXlpayurUqWPwBSv37t1jzZo1Y4IgsFWrVonK/OQDiBp79+5lO3fu1JlxwZIHEIuKilhgYCATBIFt2bJFb79XrlzR+XPDhg2NPux58eJF5urqyqpXr87OnTuns27KlClMEAT26quvapcdP37c4AtJ0tLSmCAILCoqSrtszZo1TBAEFhISojebR0lJic5sHgsXLtQ+CFhxNo/bt29rZ9uoODNExVktZsyYodP3Tz/9xKpUqcJ8fX3Zw4cPtcvtcS6NKSgoYG3atGGCILCRI0caPEcPHz5kc+fOZQsWLGCMlV9vnp6erGHDhjoPExYVFbGxY8dqH7Sr+GCgNce9e/duJggCi4+PN5h9y5Yt2pfNVPTzzz8zLy8vow8gqlQqg/0VFxezxo0bM7Vazfbs2WOwzZEjR/RmljFG85Dx+++/r7O8rKyMPf/883rfx4wxtnTpUnb69GmD/f3f//0fc3d3Z66urjrXICGVDY2ZrqCsrAzTp09HSEgIl6+DJc5n2bJl+PHHH7F+/Xr88MMPeOGFF+Dj44NLly5h586dyM/Px9tvv23wJSzGzJs3D0lJSfDw8EDr1q2xcOFCvTZt27bFgAEDtH/W/Aq4U6dOqF+/PlxcXHDx4kXs2bMHhYWFeOWVVzB8+HBte5VKhRUrViAqKgoTJkxAamoq2rRpg3v37iE9PR05OTkICQmx+o6yRp8+fWzavmrVqkhLS0NYWBiioqLwySefoGPHjtoHBTMyMlBcXCyqr4YNGyIxMRGTJ09Gu3btMGzYMNSqVQvff/89/ve//yE4OFj70BcAbNq0CatXr0bXrl3RqFEj+Pr64s8//8SuXbvg7u6uM8Z9/Pjx+L//+z989tlnaNy4MSIiIuDn54e//voLBw8exLhx4zBnzhwAwNtvv429e/di586daN26Nfr27YuCggKkpaUhJycHcXFxBuf3fv7557F27VocPXoUXbp0wfXr17Uvefn000/h5eXlkHPp4eGBffv2YciQIUhOTsauXbvQu3dvNGrUCIwxnDt3Dt9++y3y8vKwYsUKAOXX27Rp07B48WK0bNkSERERKCoqQkZGBu7fv48ePXpo78rbctxdunSBWq1GYmIi7ty5ox0+Mm3aNPj4+GDAgAH417/+hZSUFFy9ehWdOnXC5cuXkZ6ejgEDBmDLli0GMzAjd+CrVKmCbdu24cUXX0R4eDi6dOmC1q1bQ61W48qVK/jpp59w4cIF3LhxAx4eHibP62effYatW7ciNDRU52U3QPnQpY0bN6JVq1aYPHkyQkNDtWP8v/jiC8ycORPNmjVD586dUadOHe1Dsd999x0EQcAHH3wg2avlCVEEuat5JVm5ciV788032ZgxY+jONOHGw4cP2YIFC1iHDh2Yj48Pq1q1Kqtduzbr378/++qrryzub8yYMUylUjGVSqUz1ZXmS6VSsZiYGJ1tPv/8czZkyBDWqFEj5uXlxVxdXVm9evXYwIEDWXp6utF9HT16lA0bNozVrVuXVa1alXl7e7MOHTqwJUuWWDSPcHx8PFOpVHp3pg0xdmd6w4YNTKVSGZzW7PLly2zSpEksMDCQubq6slq1arFnn32WLVq0SKddQECA2WkI9+/fz8LCwpivry9zc3NjTZo0YTNnzmS5ubk67Y4ePcomTpzIWrduzWrUqME8PDxYkyZN2NixY43eCUxOTmahoaGsWrVqzN3dnTVq1IiNHDmSnThxQqddYWEhW7hwIWvRogXz8PBgPj4+rFu3bgbnf9acr5iYGPb777+zAQMGMF9fX+bp6cm6du3K9u/fr7eNPc6lOWVlZSwtLY0NGjSI1atXj7m7uzO1Ws2Cg4PZhAkT2A8//KDTvqSkhC1btow988wzzMPDg9WpU4eNHj2aXb58WXvNG7ozbclxM8bYvn37WEhIiPZO85P9/vXXX2z48OHav9NOnTqx7du3s4MHDxq8M929e3ejd6Y1bt26xd555x3WokULplarmZeXF2vatCkbOnQoS05O1pnyz5BLly6x6tWrs+rVq7PLly8bbaeZOrFfv37aZSdOnGDz589nPXv2ZIGBgczDw4O5u7uzxo0bs5EjR+rNf01IZSQwZmbQmZO4c+cOunbtiqNHj+KNN95A48aN8e6778odixBCZHXx4kU0atQIY8aMcaopzpz1uAkhlqMHEP82a9YsvPXWW9qHpGiYByGEEEIIMYeKaZRPv3XixAnt+ExW/mZImVMRQgghhBCl47KYzsvLQ1xcHMLCwuDn5weVSoWEhASjbWNjY+Hv7w8PDw+0bdtW+/CIRmZmJn777Tc89dRT8PPzw+bNm7Fo0SKMGTPGAUdDCCGEEEJ4xeVsHjk5OVizZg3atGmDyMhIrF271uiwjEGDBuHYsWNYsmQJmjZtiuTkZERHR6OsrAzR0dEAyp+EHzp0KIDyu9JvvvkmAgMDMXPmTIcdEyGEKFFAQIDJt2xWVs563IQQy3FZTAcEBODevXsAyh8cXLt2rcF2e/bswYEDB5CSkoKoqCgAQGhoKC5duoQZM2YgKioKKpUKnp6e8PT01G6nVqvh4+OjnWCfEEIIIYQQQ7gspisyNbZ5+/bt8Pb21t511oiJicGIESNw9OhR7RukKlq/fr2ofV+/ft3oq4EJIYQQQoj86tSpo50bXQrcF9Om/PrrrwgODoZKpTs0vGXLlgCA06dPGyymxbh+/ToaN26MgoICm3MSQgghhBBpVK9eHb/99ptkBXWlLqbv3LmDxo0b6y2vUaOGdr21rl+/joKCAnz++ecIDg62uh9H6tq1KzIzM+WOIRpveQH+MlNeafGWF+AvM+WVHm+ZKa+0eMt75swZjBw5EtevX6diWqmCg4PRrl07uWOIolKpuMkK8JcX4C8z5ZUWb3kB/jJTXunxlpnySou3vI7A5dR4YtWsWdPg3ee7d+9q1zuTf/3rX3JHsAhveQH+MlNeafGWF+AvM+WVHm+ZKa+0eMvrCJW6mG7VqhXOnDmjN73RL7/8AgBo0aKFHLFk4+/vL3cEi/CWF+AvM+WVFm95Af4yU17p8ZaZ8kqLt7yOUKmL6cjISOTl5WHr1q06yzds2AB/f3907tzZ5n3ExsYiIiICKSkpNvdFCCGEEEJsl5KSgoiICMTGxkq+L27HTO/duxf5+fl4+PAhgPKZOTRFc3h4ODw8PNCnTx/07t0bEydOxIMHDxAUFISUlBTs378fycnJRl/0YonExERuxg699NJLckewCG95Af4yU15p8ZYX4C8z5ZUeb5kpr7R4yRsdHY3o6GhkZWWhffv2ku6L2zvTkyZNwrBhwzBu3DgIgoC0tDQMGzYMUVFRuH37trbdtm3bMGrUKMyZMwd9+/bFTz/9hNTUVO3bD53JV199JXcEi/CWF+AvM+WVFm95Af4yU17p8ZaZ8kqLt7yOIDBTbz0hRmk+6Rw/fpybO9NZWVncZAX4ywvwl5nySktpeXfvBjw8gJ49jbdRWmZzKK/0eMtMeaXFY16p6zUqpq3EYzFNCHFuvXsDvr7Ali327Tc9HQgKApo3t2+/xD7Onj2rHRJJSGXi7e2NJk2amGzjiHqN2zHThBBCLKNSAVLcPpk6FRg1Cpg/3/59E9ucPXsWTZs2lTsGIZL5448/zBbUUqNimhBCnIRKBTwxU6hd2OFZbiIRzR1pnt7WS4gYmjcbKuG3LlRMO5GkpCSMGzdO7hii8ZYX4C8z5ZWW0vKKKaaVltkcyisOT2/rJYQ33M7moRQ8zTOdlZUldwSL8JYX4C8z5ZWW0vK6uAClpabbWJtZrqdvlHaOzeEtLyG8cuQ80/QAopXoAURCCG8GDCgvetPT7dtvYCAwYgSwYIF9+yW2o59VpLISe2074nuA7kwTQoiTkPLWidi+P/wQePBAuhyEEOJoVEwTQogTkeJhQUv6jI0Fzp2zfwZCCJELFdOEEEIIIYRYiYppJxIRESF3BIvwlhfgLzPllRZveQHrM8v19A1v55i3vEQ+Bw8ehEqlQkJCgtxRiBlUTDuRKVOmyB3BIrzlBfjLTHmlpbS8YgpeazJbOnTEnoW30s6xObzlrWyys7MxdepUtGjRAtWqVYObmxv8/f3x0ksvYd26dSgqKnJYlosXL0KlUiEmJsZkO4Emclc8mmfaiYSFhckdwSK85QX4y0x5paXEvOZ+Ljsisz2LaSWeY1N4y1uZzJs3DwkJCWCMoUuXLnjhhRfg7e2NGzdu4NChQxg/fjxWrVqFn376ySF5NEWysWK5c+fOyM7ORq1atRySh1iPimlCCHEiUg3HsKRfmpCVONqCBQsQHx+PBg0aIC0tDR07dtRr8/XXX+O///2vwzJpZiY2NkOxh4cHvQqeEzTMgxBCnIRUvy2m30ITJbtw4QISEhLg6uqKPXv2GCykAeDFF1/Enj17dJZt3rwZ3bp1Q7Vq1aBWq9GyZUssWrQIjx8/1ts+ICAAgYGBKCgowIwZM9CgQQO4u7ujSZMmWLJkiU7b+Ph4NGrUCACwceNGqFQq7dfGjRsBGB8z3b17d6hUKpSWlmLhwoVo0qQJ3N3d0aBBA8TFxekNVTE3nETT35PKysqwcuVKdOzYEd7e3vDy8kLHjh2xatUqvQ8A1u5j3bp1CAkJgZ+fHzw8PODv74/evXtj8+bNBvtRKiqmbcTTGxB37NghdwSL8JYX4C8z5ZWW0vKKuSNsbWa57jYr7Rybw1veymDDhg0oKSnB4MGD8cwzz5hs6+rqqv3/mTNnIjo6GmfPnsWoUaMwdepUMMbw7rvvIiwsDMXFxTrbCoKA4uJihIWFYdu2bQgPD8eECRPw6NEjzJo1C/Hx8dq2PXr0wBtvvAEAaNOmDeLj47Vfbdu21evXkOjoaKxYsQKhoaGYNGkSPDw88P777+PVV1812N7U2GtD60aMGIEpU6YgJycHEyZMwGuvvYacnBxMnjwZI0aMsHkfM2fOxPjx43H79m0MHz4cb731Fl588UXcuHEDX375pdF+xHLkGxDBiFWOHz/OALDjx4/LHUW0YcOGyR3BIrzlZYy/zJRXWkrL+9JLjA0YYLqNNZkbN2YsLk5cW4Cxo0ct3oVRSjvH5jg6L48/q+ytR48eTBAElpSUJHqbzMxMJggCCwwMZLdv39YuLykpYeHh4UwQBLZgwQKdbRo2bMgEQWDh4eFFg54qAAAgAElEQVSssLBQu/zWrVusevXqrFq1aqy4uFi7/OLFi0wQBBYTE2MwQ0ZGBhMEgSUkJOgsDw0NZYIgsA4dOrB79+5pl+fn57PGjRszFxcXdv36de3yCxcumNxPaGgoU6lUOsuSk5OZIAisU6dOrKCgQGcf7du3Z4IgsOTkZJv24evry+rVq8cePXqk1z4nJ8dgPxWJvbYd8T1Ad6adCG+/NuEtL8BfZsorLd7yAvxlprzEnBs3bgAA6tWrJ3qb9evXAwDee+89nQcAXVxcsGzZMqhUKiQlJeltJwgCli9fDjc3N+0yPz8/RERE4MGDB/jjjz+0y5mNv85ZunQpqlevrv2zWq3Gyy+/jLKyMmRlZdnU97p16wAAixYtgoeHh84+NENWDB2/JVQqFVxdXQ0O/6hZs6ZNfTsaPYBICCHEZvQAYuUzcSLw11+O25+/P7BqleP2Z8qJEycgCAJ69Oiht65p06bw9/fHxYsX8eDBA/j4+GjXVa9eHYGBgXrb1K9fHwBw7949u+QTBAEdOnTQW675wGDrfk6cOAEXFxeEhobqrQsNDYVKpcKJEyds2sfLL7+M5cuXo3nz5hg2bBief/55PPvss6hWrZpN/cqBimlCCHESUhWxgkDFdGWklMLWVnXq1EF2djauXr0qepvc3FwAQO3atY32efXqVeTm5uoU08YKwSpVysut0tJS0RnM8fb2lmw/ubm5qFmzJlxcXAzuo1atWsjJybFpH//v//0/NGrUCOvXr8eiRYuwaNEiVKlSBeHh4Vi2bJnBDyVKRcM8CCHEiUgx8wbN5kGUrFu3bgCAb7/9VvQ2mqL4+vXrBtdrlvNwF1UzjKKkpMTg+vv37+stq1atGu7evWuwKC8pKUFOTo7Ohwhr9qFSqfDGG2/g5MmTuHnzJr788ktERkZi586d6NOnj94DnkpGxbQTMfeWJaXhLS/AX2bKKy3e8gLWZ7bkbrM9i2/ezjFveSuDmJgYVK1aFV9++SXOnDljsq1mWrl27dqBMYaDBw/qtTl37hyuXr2KwMBAnYLSUpq7vva8W22Ir68vAODKlSt6654cx63Rrl07lJaW4vvvv9dbd+jQIZSVlaFdu3Y27aMiPz8/REZGYvPmzejRowfOnj2L06dPmz4wBaFi2onw9uYt3vIC/GWmvNLiLS9gXWY5h3nwdo55y1sZNGzYEPHx8SgqKkJ4eDiOHz9usN3evXvRp08fAMDYsWMBAPPnz9cZzlBaWoq3334bjDGMGzfOplyaAvTy5cs29WOOt7c3goODkZmZqfNhorS0FG+++SYKCwv1ttEc/6xZs/Do0SPt8oKCArzzzjsAoHP8lu6jqKgIhw8f1ttvcXEx7t69C0EQ4O7ubuUROx6NmXYi0dHRckewCG95Af4yU15pKS2vmCLW2sxyDfVQ2jk2h7e8lcWsWbNQUlKChIQEdOzYEV26dEH79u3h5eWFmzdv4tChQzh37pz2hS4hISGIi4vD0qVL0aJFCwwZMgRqtRp79+7F6dOn0a1bN8yYMcOmTF5eXnj22Wdx6NAhjBo1Co0bN4aLiwsGDBiAli1bmtzW0plAZs6ciTFjxuC5557DkCFD4O7ujoyMDJSWlqJ169Y4deqUTvvo6Gjs3LkTW7ZsQfPmzTFgwAAIgoAdO3bg4sWLGD58uN61bMk+CgoK0K1bNzRu3Bjt2rVDw4YNUVhYiG+++QbZ2dno378/mjVrZtExyomKaUIIIQ5FDyASOcyePRtDhw7FypUrkZGRgQ0bNqCwsBC1atVCmzZtMGvWLIwcOVLbfvHixWjbti1WrFiBTZs2obi4GI0bN8aCBQvw1ltvaR/20zD3whJD6z/77DNMnz4de/fu1c7A0aBBA5PFtLG+TK0bPXo0ysrK8P7772PTpk2oUaMGBgwYgAULFmDw4MEGt0lJSUFoaCjWrVuH1atXQxAEBAcHY8aMGZg4caJN+/Dy8sKSJUuQkZGBH374ATt37oSPjw+CgoLwySefaO+M80Jgtk506KSysrLQvn17HD9+XGfcECGEKNVLLwFVqwLbt9u332bNyvt+/33zbQUBOHIECAmxbwZiGP2sIpWV2GvbEd8DNGbaiWRmZsodwSK85QX4y0x5pcVbXsD6zHI9gMjbOeYtLyHEPCqmncjSpUvljmAR3vIC/GWmvNJSWl4xBa8jMtvz96FKO8fm8JaXEGIeFdNOJDU1Ve4IFuEtL8BfZsorLSXmNXdX2NrMcj2AqMRzbApveQkh5lEx7UTUarXcESzCW16Av8yUV1q85QWszyzX0ze8nWPe8hJCzKNimhBCiE3oDYiEEGdGxTQhhDgJ3uZuKiwEHjyQOwUhhJhGxbQTsXWCeUfjLS/AX2bKKy0l5jV3F9kRmcUW9bNnA+Hhptso8RybwlteQoh5VEw7kQYNGsgdwSK85QX4y0x5pcVbXkBZmQsKgLw8022UlFcM3vISQsyjNyDaKDY2FtWrV0d0dLTiXxM7depUuSNYhLe8AH+ZKa+0lJjX3F1hazNLNYREqrxy4S0vIbxKSUlBSkoK7t+/L/m+qJi2UWJiIr1VihBCCCFEQTQ3OTVvQJQSDfMghBAnItXMG7z1Swgh9kLFtBPJzs6WO4JFeMsL8JeZ8kqLt7yA9ZnlmimEt3PMW15CiHlUTDuRuLg4uSNYhLe8AH+ZKa+0eMsLWJdZzrvHvJ1j3vISQsyjYtqJrFixQu4IFuEtL8BfZsorLaXlFXP3WGmZzaG8hBC5UTHtRHibkom3vAB/mSmvtJSY19xdZCVmNuXwYb7y8nZ+K4OtW7di6tSp6NatG3x8fKBSqTBq1CiT25SWlmLt2rV4/vnn4evrC7VajaCgIAwfPhxnz561KkdRURHWrVuH/v37w9/fHx4eHvDy8kLjxo0RFRWF5ORkFBUVWdU3kRfN5kEIIcSh7Dm+esQIQOGzkhKZzZ8/Hz///DO8vb1Rr149ZGdnQzDxqTIvLw8DBgxARkYG2rZti5iYGLi7u+Pq1avIzMzE2bNn0aRJE4sy/Pbbbxg0aBD++OMP1KpVC7169ULDhg0hCAIuXbqEgwcPIi0tDUuWLMHPP/9s6yETB6NiuoLhw4fj4MGDKCgoQJ06dfD2229jwoQJcscihBDFk2ueaULMSUxMRP369REUFITvv/8ePXr0MNn+1VdfRUZGBj799FODNUBJSYlF+7927RpeeOEF3LhxA3FxcUhISICbm5tOG8YYdu7ciWXLllnUN1EGGuZRwdy5c3H16lU8ePAAn3/+OaZNm4YLFy7IHctulixZIncEi/CWF+AvM+WVltLyiilMlZbZPL7y8nd++de9e3cEBQUBKC9aTcnKykJqaiqGDx9u9GZalSqW3Yf897//jRs3bmDMmDFYvHixXiENAIIgYODAgcjIyNBZfvDgQahUKiQkJOB///sf+vTpA19fX6hUKly+fBkAUFhYiEWLFqFly5bw9PREtWrV8Pzzz2Pz5s16+6nYnyEBAQEIDAzUWbZhwwaoVCps3LgRX331Fbp06QIvLy/UqFEDQ4cOxblz5yw6H5UR3ZmuIDg4WPv/Li4u8PHxgbe3t4yJ7KugoEDuCBbhLS/AX2bKKy3e8gLWZ5Zvnmm+zjGP14Qz+eKLLwCUv/AjNzcXu3btwpUrV1CzZk306tVLW5SLVVBQgJSUFAiCgNmzZ5tt7+LiYnD5kSNHsHDhQjz//POYMGECbt26BVdXVxQVFSEsLAyZmZlo3rw5pkyZgvz8fKSlpSE6OhonTpzA4sWL9fozNczF2Lpt27Zh7969GDRoEHr27IkTJ07gyy+/REZGBo4cOYKmTZuaPb7KiorpJ7z88svYtm0bACA1NRW1atWSOZH9GPskqlS85QX4y0x5paW0vGIKXmszWzIcw75tlXWOzVHaNUF0/fTTTwCAS5cuISgoCHfv3tWuEwQBEydOxEcffQSVStwv9o8dO4bi4mI0aNBA746vJb755huDw04WLlyIzMxM9O/fH9u3b9fmmjNnDjp16oSlS5eif//+eO6556zet8auXbvw1VdfoV+/ftplH330EWJjYzFp0iQcOHDA5n3wiorpJyQnJ6OsrAzp6emIiYnByZMn6elrQggxgd5SWAl16ADcuOH4/dauDRw75vj9/u3WrVsAgOnTpyMyMhLz589HvXr1cOTIEbz++utYuXIl/Pz8MHfuXFH93fj7HNatW9fg+k8++UTbBigv2EePHq1XeLdt29bgsJN169ZBpVLhgw8+0Cnwn3rqKcyePRsTJkzAunXr7FJM9+rVS6eQBoApU6bgo48+wnfffYfLly87bb1ExbQBKpUKAwcORFJSEtLT0zFlyhS5IxFCiM14fJiPCnWZ3LgB/PWX3CkcrqysDED5sM/Nmzdrhzy88MIL2Lp1Kzp06IBly5bh3//+N6pWrYqDBw/i4MGDOn0EBgbilVdeEbW/Tz/9FKdOndJZ1q1bN71iulOnTnrbPnz4EOfPn0f9+vXRuHFjvfW9evUCAJw4cUJUFnNCQ0P1lqlUKnTt2hXnz5936puP3BbTeXl5mDdvHk6ePIkTJ07gzp07mDt3rsFPi3l5eXjvvfeQlpaGu3fvolmzZnjnnXcQFRVlch8lJSXw8vKS6hAcLicnh6thK7zlBfjLTHmlpcS85opTR2S27zCPHADKOsemKPGaMKh2befa79+qV68OAOjfv7/e2OE2bdqgYcOGuHjxIrKzs9GyZUt8//33mDdvnk677t27a4vp2n8fz7Vr1wzur2KhGxMTg40bNxpsV9vAecnNzTW6ruJyTTtbPf300w7ZD4+4nc0jJycHa9asQXFxMSIjIwEYHzQ/aNAgbNq0CfHx8di3bx86duyI6OhopKSkaNvcvHkTW7duRX5+PkpKSrBlyxYcPXoUvXv3dsjxOMLYsWPljmAR3vIC/GWmvNJSYl5zxakjMtv3DrnyzrEpSrwmDDp2DLh61fFfMg7xAIBmzZoB+KeofpKvry8YY3j06BGA8lnAysrKdL6+++47bfuOHTuiatWquHLlCv7880+T+zY104ih+qZatWoAoDNMpKLr16/rtAOgHQpibHq/+/fvG81w8+ZNg8s1+6+4H2fDbTEdEBCAe/fuISMjA4sWLTLabs+ePThw4ABWrVqFCRMmIDQ0FKtXr0bv3r0xY8YM7a90gPKB9P7+/njqqaewYsUKpKenw9/f3xGH4xDx8fFyR7AIb3kB/jJTXmkpLa+YIRPWZpbqAUTz4u3ZmeSUdk0QXS+88AIA4JdfftFb9/jxY5w9exaCICAgIEBUfx4eHhgxYgQYY5g/f749o8Lb2xtBQUG4evWqwenpNNPstWvXTrvM19cXALTT6lV07tw5PHjwwOj+nhzOApS/KTIzMxOCIKBt27aWHkKlwW0xXZGpT3Pbt2+Ht7c3hg4dqrM8JiYG165dw9GjRwGU//ri0KFDuH//Pu7evYtDhw6ha9eukuZ2tIrfUDzgLS/AX2bKKy2l5RVTxFqTWb5p8QBAWefYHKVdE0TX4MGDUbduXWzevFk7s4dGfHw8Hj58iB49euCpp54S3eeCBQtQu3ZtbNy4ETNnzkRhYaFem7KyMpOFrDFjx44FY0zv5mBOTg7+85//QBAEnd+GBAcHw8fHBzt37sTt27e1yx89eoRp06aZ3Nd3332H3bt36yxbsWIFzp8/jx49eqB+/foW568sKkUxbcqvv/6K4OBgvWlsWrZsCQA4ffq0Tf3369cPEREROl8hISHYsWOHTrv9+/cjIiJCb/vJkycjKSlJZ1lWVhYiIiKQk5Ojs3zu3Ll6E/5fvnwZERERyM7O1lm+fPlyzJgxQ2dZQUEBIiIikJmZqbM8JSUFMTExetmioqLoOOg46Dgq0XH88kuMXoFqj+O4dCkCjx6JOw4gApcuiTuO3bsjkJdXef8+HHkczmzHjh0YM2aM9qUpQPm8zZplFf/O1Go1NmzYAEEQ0K1bN4wYMQJvv/02unbtiiVLluDpp5/Gp59+atH+69atiwMHDqBp06b473//i/r162P48OGYOXMm4uLiMHr0aAQEBGDHjh0ICAhAw4YNRfetybZz5060bt0acXFxmDJlCpo3b47Lly8jLi4OXbp00bavUqUK3nzzTeTm5qJt27aYMmUKXn/9dbRs2RL5+fmoW7eu0RuUERERiIyMRFRUFP7973+jX79+mD59OmrWrImVK1dadE7s6YcfftB+f6SkpGhrscDAQLRp0waxsbHSh2CVwO3bt5kgCCwhIUFvXZMmTVjfvn31ll+7do0JgsAWL15s1T6PHz/OALDjx49btT0hhDjaiy8yNniw/ft95hnG3nhDXFuAse++E9d20iTGWrc23x9jjO3fz9gXX4jr15nQzyrG4uPjmSAITKVS6XwJgsAEQWCBgYF625w6dYoNGTKE+fn5MVdXV9awYUM2adIkdv36datzPH78mCUlJbHw8HBWt25d5ubmxtRqNQsKCmJDhw5lycnJrKioSGebjIwMo/WNRmFhIVu4cCFr0aIF8/DwYD4+Pqxbt24sNTXV6DZLly5lQUFB2mObOXMmKygoYAEBAXrnY/369UwQBLZx40a2e/duFhISwjw9PZmvry8bMmQIO3v2rNXnxBZir21HfA9U+jvT5B9P3sFQOt7yAvxlprzSUlpeMcMmrMls6TAPsWOmxfVbnjc1FfjoI8tyyEFp14Qz0DwkWFpaqvOleWDw/Pnzetu0atUKaWlpuHXrFh4/foyLFy/i448/Njpzhhiurq4YO3YsvvrqK/z1118oLCxEfn4+zp07hy1btmDEiBGoWrWqzjbdu3dHWVkZ5syZY7RfNzc3zJo1C7/88gsKCgqQm5uLQ4cOmZyxbMaMGTh37pz22BYvXgwPDw9cuHDB4PnQ6NevH44cOYK8vDzcvXsXaWlpBqflczaVvpiuWbMm7ty5o7dc81ajmjVrOjqSbLKysuSOYBHe8gL8Zaa80lJaXjFFrCMyiy2mxbVT1jk2R2nXBCHEdpW+mG7VqhXOnDmjMzAf+OdJ3RYtWsgRSxYff/yx3BEswltegL/MlFdaSsxr7m6vEjObxlde/s4vIcScSl9MR0ZGIi8vD1u3btVZvmHDBvj7+6Nz58429R8bG4uIiAidOasJIYQYZ99hHv+05fENj4QonSAIRt/joWSahxEd8QAit29ABIC9e/ciPz8fDx8+BFA+M4emaA4PD4eHhwf69OmD3r17Y+LEiXjw4AGCgoKQkpKC/fv3Izk52eYLJDExkaY6IoRwQ6qCU755pqXrkxACvPLKK6Jfj64k0dHRiI6ORlZWFtq3by/pvrgupidNmoRLly4BKP/klJaWhrS0NAiCgAsXLmjfEb9t2za8++67mDNnDu7evYvg4GCkpqZi2LBhcsYnhJBKQaoHEKXOQQgh9sD1MI8LFy5on8at+GRuaWmptpAGAE9PTyQmJuLatWsoLCzEiRMnnLKQNjRPqZLxlhfgLzPllZYS85orOK3JLO9d4fK8vAzzUOI1QQixDdfFNLHMlClT5I5gEd7yAvxlprzS4i0v4JjM9i16+TrHPF4ThBDTqJh2ImFhYXJHsAhveQH+MlNeaSkxr7lC1prM8g7zCJOgT+ko8ZoghNiGimlCCCEOJefDin8/ZkMIIXZDxTQhhDiRyvqQntjjCgiQNAYhxAlxPZuHEsTGxqJ69eraKViUbMeOHRg4cKDcMUTjLS/AX2bKKy3e8gKOyWzJ3WbzRfIOAPycY7muiTNnzjh8n4RIydw1nZKSgpSUFNy/f1/yLFRM24ineaZTUlK4+sHOW16Av8yUV1q85QWszyzV0A3zbVPAUzHt6GvC29sbADBy5EiH7ZMQR9Jc40+ieaaJJDZv3ix3BIvwlhfgLzPllZbS8oqZPs6azPI+gCg+rxIeUnT0NdGkSRP88ccf2pebEVKZeHt7o0mTJnLHoGKaEEKchSAAZWX279fSItW+wzyIOUooNgipzOgBREIIcRK8vNikIiny/vgjMGCA/fslhDgnKqYJIcRJSFVMK+F14mKOTbP+zz+B9HT7ZyCEOCcqpp1ITEyM3BEswltegL/MlFdaSssrpuB0RGb7FtMxEvQpHaVdE2LwlpnySou3vI5AxbQT4e3NW7zlBfjLTHmlpbS8Yu4gOyKz2MJX3B3vf/LyML5aadeEGLxlprzS4i2vI1Ax7USUPg/2k3jLC/CXmfJKi7e8gLIyiyu6y/PyUEgDyjq/YvGWmfJKi7e8jkDFNCGEEJvJ+Ypwsf3yMhSEEMIXmhrPRjy9AZEQQqQg1QOIvNxtJoQojyPfgEh3pm2UmJiI9PR0LgrpzMxMuSNYhLe8AH+ZKa+0eMsLWJdZqnmmxbX7Jy8PxbezXBNyorzS4iVvdHQ00tPTkZiYKPm+qJh2IkuXLpU7gkV4ywvwl5nySou3vIBjMtt3uMU/eXkY5kHXhPQor7R4y+sIVEw7kdTUVLkjWIS3vAB/mSmvtHjLC1iXWao7wuL6te0cT5oElJTY1IVFnOWakBPllRZveR2Bimknolar5Y5gEd7yAvxlprzS4i0v4JjM9r1DbFveVauA4mI7RRGBrgnpUV5p8ZbXEaiYJoQQ4lBKGG5BCCH2QsU0IYQQh5KrmKYinhAiBSqmnciMGTPkjmAR3vIC/GWmvNLiLS9gfWb5ClW+zrEzXRNyobzS4i2vI1Ax7UQaNGggdwSL8JYX4C8z5ZUWb3kB6zJLNc+0OA2syiAXZ7km5ER5pcVbXkcQGKNffFkjKysL7du3x/Hjx9GuXTu54xBCiFkREeVF586d9u23dWugWzdgxQrzbQUB+OILQMzU/FOnAocOAadOme6PMWDCBODnn4GjR423LSoC3NzK9z9iRPl2ggAUFAAeHsDNm8DTT5vPRQjhhyPqNbozTQghxKHkvoVj7C527dqOzUEIqRyomCaEEGITJQyxsCSD3MU8IaRyoWLaiWRnZ8sdwSK85QX4y0x5pcVbXsAxmS0pZiu2zc0FCgufbJFtsK09LF4MzJ5t3z7pmpAe5ZUWb3kdgYppJxIXFyd3BIvwlhfgLzPllRZveQHHZLak6K14x7lXL2DBgidbiM9rabF9/Djw44+WbWMOXRPSo7zS4i2vI1Ax7URWiHk6SEF4ywvwl5nySou3vIBjMlt7Z/rhQ0N3pv/Jq4ThJubQNSE9yist3vI6AhXTToS36Wx4ywvwl5nySou3vID1ma0tkG3fj/RT+dmTM10TcqG80uItryNQMU0IIcQmjipOze3HXJGuWf9kOx7uaBNClIuKaUIIIQ4l9s60o4pcmt2DEGILKqadyJIlS+SOYBHe8gL8Zaa80uItL+CYzPYtXv/Ja4/iu1Wrf/5fiiKbrgnpUV5p8ZbXEarIHYB3sbGxqF69OqKjoxEt5pVeMiooKJA7gkV4ywvwl5nySou3vIBjMtu3SP0nr9hhHqb88ouNccyga0J6lFdavORNSUlBSkoK7t+/L/m+6HXiVqLXiRNCeNO/P6BS2f914m3bAl26AB9/bL6tIADr1gExMebbTpsGZGT8U+A2awaEhwMffKDbH2PAa68BJ06Ynsru0SNArQZSUspfZ/7k68Q1d7Y1PxWHDgUePAC+/tp8VkKIMtHrxAkhhHDBUbN5EEKI0lAxTQghTkIps1ZIVUyL7Zdm8yCE2BMV004kJydH7ggW4S0vwF9myist3vIC1me2pCC1djYPw/vIMbPe/H4deafcma4JuVBeafGW1xGomHYiY8eOlTuCRXjLC/CXmfJKi7e8gPWZ5Xtpi3TnWIoi25muCblQXmnxltcRqJh2IvHx8XJHsAhveQH+MlNeaSktr5ji0JrM8g6TiNf+nxTzV9v72JR2TYjBW2bKKy3e8joCFdNOhLdZR3jLC/CXmfJKS2l5xRSGjshsy11s/WMQn9eaO832vjuttGtCDN4yU15p8ZbXEaiY/ltRURFiYmLQoEEDVKtWDSEhIfjhhx/kjkUIIYpnScGpmcrOEfvSOHwYKC21vA96MJEQIgYV038rKSlBo0aNcOTIEeTm5mLixImIiIjAo0eP5I5mF2fOAAcOyJ2CEFJZiS08bSmmTe3D1LquXYFr16zblhBCzKFi+m9qtRqzZ89GvXr1AACjR49GWVkZzp07J3My+6hXD3j//SS8/DJw9arcacRJSkqSO4LFeMtMeaXFW17A+szyjVcWn9eWO+KZmUBQkO6yPn0ASyc2cKZrQi6UV1q85XUEKqaNyM7OxqNHjxD05L+enPL2BoKCsvDvfwMTJgDJyXInMi8rK0vuCBbjLTPllRZveQHrMlt6Z9e+45Btzysmz+3bwPnzusu+/hrw87Ns385yTciJ8kqLt7yOQMW0AQUFBRg1ahRmz54NtVotdxy7+fjjj9G8ObBrF/DHH+Wv883NlTuVcR+LeTexwvCWmfJKi7e8gPWZLbkzbd9iWrpzLDZnWRkg9pksZ7om5EJ5pcVbXkegYvoJxcXFGDp0KFq0aIFZs2bJHUcSVaoACQnA+PFAZCSwf7/ciQghPLOkQLa0mH6yraltzfVrTREv5q47Y8CJE5b3TQipHLgtpvPy8hAXF4ewsDD4+flBpVIhISHBaNvY2Fj4+/vDw8MDbdu2xebNm/XalZWVYdSoUXB1dXWKMUHPPVd+l/rrr8vvUtNLjQgh1rC0mJablC+Y+eUXy9oTQvjHbTGdk5ODNWvWoLi4GJGRkQAAwci/0oMGDcKmTZsQHx+Pffv2oWPHjoiOjkZKSopOu9deew03b95EamoqVCpuT41FPD2BDz4AJk0Chg8Htm+XOxEhhDfyjpm2LoOlfYrN3KqV/XMQQpSN24oxICAA9+7dQ0ZGBhYtWmS03Z49e3DgwAGsWrUKEyZMQGhoKFavXo3evXtjxowZKCsrAwBcunQJSUlJ+PHHH1GrVnldjv8AACAASURBVC14e3vD29sbhw8fdtQhSS4iIsLouo4dgd27gR9/BMaOBe7fd2AwI0zlVSreMlNeafGWF7A+sxTDPJ5sa7hgFp9XiiIeAEaMEN/Wma4JuVBeafGW1xG4LaYrYib+hdy+fTu8vb0xdOhQneUxMTG4du0ajh49CgBo2LAhysrKkJ+fj4cPH2q/nnvuOUmzO9KUKVNMrndzAxYtKh9LPWgQ8M03DgpmhLm8SsRbZsorLd7yAtZllnLM9JP0t7Uurz23SUsT34+zXBNyorzS4i2vI1SKYtqUX3/9FcHBwXrDNlq2bAkAOH36tE399+vXDxERETpfISEh2LFjh067/fv3G/w0N3nyZL3x2VlZWYiIiEDOE4OY586diyVLlugsu3z5MiIiIpCdna2zfPny5ZgxY4bOsq5duyIiIgKZmZk6y1NSUhATE6P9c5cu5WOpJ02KQr9+O5CfL89xhIWFGTyOgoICUcehERUV5bC/j2bNmon++1DCcYSFhdl8XTnyOMLCwiT7/pDiOMLCwgweByDd97mp4zh50vxxhIWFWXxdnT0bgUePxB1HUVEEbtwQdxzp6REoKNA9jt9/f/Lvo/wcf/NNFO7dE3ddrVs3GRXnp2ZMM91XBIAcneXnzhn/+wD0jwMw/fehuSaU8O+V2OsqLCxMEf9eiT0OzTmW+98rscehySv3v1dij0OT98nj0JDzOFJSUrS1WGBgINq0aYPY2Fi9fuyOVQK3b99mgiCwhIQEvXVNmjRhffv21Vt+7do1JggCW7x4sVX7PH78OAPAjh8/btX2vPjmG8Z69GDs8GG5kxBCbFFaylhEBGP9+9u/786dGRs3TlxbDw/GEhPFtZ0+nbHg4H/+/MwzjMXG6rbR/BR7/XXG2rUz3A/A2KVLjN2/X/7/X3zxz3YAY/n5//x/xZ+Kgwcz9uKL5f//5Ze66yq2FwTd/gghyuGIeq3S35kmtnnhBWDbNiApCZg5E3j8WO5EhBBrlJYCLi7SzaYh3zAP+Wky3b4tbw5CiDwqfTFds2ZN3LlzR2/53bt3teudxZO/GhGrevXyYrpLFyA8HPjpJzsHM8LavHLiLTPllZaS8mqKaXOsyeyoMdOGPwiIy2vthwhLsj7zjPk2SromxOItM+WVFm95HaHSF9OtWrXCmTNntLN2aPzy92SgLVq0kCOWLJ6cCtBSAwYAmzcDK1YAs2ZJf5fa1rxy4C0z5ZWWkvKKLaatyezIl7boS9H2K7YvsYW1pe0KCsr/27698bZKuibE4i0z5ZUWb3kdodIX05GRkcjLy8PWrVt1lm/YsAH+/v7o3LmzTf3HxsYiIiKCi4vL0ItqLFWzJrBxY/lUeuHhwMmTdghmhD3yOhpvmSmvtJSUV2wxbU1mS+762n+YyWaHDP3Q5DY1K5gmR1aW8TZKuibE4i0z5ZUWL3k1DyM64gHEKpLvQUJ79+7VTmUHlM/MoSmaw8PD4eHhgT59+qB3796YOHEiHjx4gKCgIKSkpGD//v1ITk42+qIXsRITE9GuXTubj4U3gwYBXbsC06YBzZsD77wDVK0qdypCiDFii2lrSflWQXtta8u+NP+/a5fj9k8IsV50dDSio6ORlZWF9qZ+XWQHXBfTkyZNwqVLlwCUv/0wLS0NaWlpEAQBFy5cQIMGDQAA27Ztw7vvvos5c+bg7t27CA4ORmpqKoYNGyZnfO499RSQklL+FR4OLFsGONGoGUK4oimmpXr7oFQvbTH1Zw0x/Wnm3jDU3tT29riTXlAAqNW290MIUSaui+kLFy6Iaufp6YnExEQkJiZKnMj5CEL527969ACmTgU6dQLeekvaO2CEEMtJeWfaUQ8gmttOrpk+xBTjSpyFhBBiH5V+zDT5h6EJ0O2lTp3yt4DVqgX07w/8+aftfUqZVyq8Zaa80lJS3opT45kq7KzJLO+Y6RhRhWrF/Uo1PaAYSromxOItM+WVFm95HYGKaSdS8a1FUhAEYOxYYOVKIDa2/L9PTKJiEanzSoG3zJRXWkrKqymmXVxMf19am1mqMdMV2xougsNMrLN+v7ZsY8iKFcDQocq6JsTiLTPllRZveR2BimknEh0d7ZD9BAQAO3eW/8COjAT+HtZuMUfltSfeMlNeaSkpb8ViuqTEeDtrMjtqmAdgaNtoyYZQVCzQS0vNt6+Y48kPLMuXA1u3KuuaEIu3zJRXWrzldQQqpokkVCpgyhTggw+A118H1q+nMYOEyKliMS2mMLSEI4tpaxmamcOSbSx9Xv3Jd4X98Ydl2xNC+EHFtI14mmdaDo0bA199Bdy6Vf4rzmvX5E5EiHPSFNNVqpi+M20NS8dMWzubh7FlxmbpMNZO7HJj+7Ol7c2b4vsjhFjPkfNMUzFto8TERKSnp3Pxa4/MzExZ9uviAsycCSQklI+p3rhR3A9TufLagrfMlFdaSsor9s60tZltKZBt24/leeV8ANHfXznXhFhKuo7FoLzS4iVvdHQ00tPTHTKTGxXTTmTp0qWy7r958/K71DduAFFR5f81Re681uAtM+WVlpLyir0zbU1m+78i3Ph+9C0VPZuH1MNLxPRfWrrU7L99SqOk61gMyist3vI6AhXTTiQ1NVXuCKhSpfwu9ezZwKhR5dPpGaOEvJbiLTPllZaS8oq9M21NZinHTJtvazxvfr74/dib8dypqFPHkUlsp6TrWAzKKy3e8joCFdNORK2gV3C1bAns3g38/DMwejRw965+GyXlFYu3zJRXWkrKW1pa/mHWXDFtTWapxkyLowZj+hkuXAC8vGzv3ZKs4s5D+fktLLRt6lBHUtJ1LAbllRZveR2BimkiG1dX4D//ASZPLn848euv5U5ESOUl5QOIgOPGTIvNUFSkv87SIt5YVnsM02jeHNi82fZ+CCHyo2KayK5zZ2DXLmDPHmDSJCAvT+5EhFQ+SpkaD7Ct8Da0rVRjoY31a2yYhiU5Ll8Gzp0z3UbOByUJIeJRMe1EZsyYIXcEo9Rq4MMPgcGDgYgI4LvvlJ3XGN4yU15pKSmv2Je2WJNZ3nmm/8lrbfFpr6nxxCnPW1ICzJkDPHpk7/7tT0nXsRiUV1q85XUEKqadSIMGDeSOYFavXuVvT9yxAzhypAEePJA7kWV4OMcVUV5pKSlvxWEepu5MW5NZ3jHT/+Q1VxQ78mUxxvele355GH6qpOtYDMorLd7yOgIV005k6tSpckcQxdsb+OgjYPHiqRg4EPj2W7kTicfLOdagvNJSUl6xd6atzeyIO9OGi/apinm7qrgPFcq5JsRS0nUsBuWVFm95HYGKaRvRGxCl060bkJ4ObNsGTJ0q7zRXhPBO7J1pa8j1inBLMjgin9zngBDyD0e+AbGK5Huo5BITE9GuXTu5Y1RaXl7Axx8D33wD9O8PLFgAhITInYoQ/ijlAUQpHlYU2x8Vu4Q4j+joaERHRyMrKwvt27eXdF90Z9qJZGdnyx3BIhXz9u5dfof644/LC2p7FwP2wvM55gHltV5xcfl0lOaGeViTWapiWtywiX/yKmn2C+NZDJ/fixeVW+wr6ToWg/JKi7e8jkDFtBOJi4uTO4JFnsxbvTrw2Wfl01JFRABnz8oUzATez7HSUV7rFRWVF9PmhnlYk9mRDyDqbxunXVZxnalM9ii6L1823a/xYzR8fgMDDfepBEq6jsWgvNLiLa8jUDHtRFasWCF3BIsYyisIwNixwMqVwIwZwOLF5XfclKIynGMlo7zW0xTT5u5MW5NZEMS/zc/SQtZ8gWw475PFrKki3priPihIf5m4ae6Mn99LlyzP4QhKuo7FoLzS4i2vI1Ax7UR4m87GVN6GDYHt24H69YF+/YAzZxwYzITKdI6ViPJaT+ydaWunxpPipS3iNLDruGqxrH8VuPHzu3QpUFBgbb/SUdJ1LAbllRZveR2BimnCLUEAXn4Z2LSp/C71Z5/JnYgQ5RJ7Z9oa8s4zrUyHDwMdO1q2ze7dgKenNHkIIdKhYppwr06d8pe8/PknMGwYcOWK3IkIUZ6KxbScD/Da/wHE8v7EtHVkEf/XX8CxY+LaXrsmbRZCiLSomHYiS5YskTuCRSzJW6UKEB8PzJsHvP56+UtfrP81rPUq8zlWAsprPbHDPKzJbEmRauvDf/r7WmJwnT0eMpSm+NY/v8nJun8eN06K/VpPSdexGJRXWrzldQQqpp1IgRIH45lgTd5mzYBdu8qLhsGDgdu3JQhmgjOcYzlRXus9fixumIc1mW15qNB2BQb7M/dqcfnon98vv9T987p1QEKCg+KIoKTrWAzKKy3e8joCFdNOJEFJ/zqLYG1elar87nRCAhAVBXz1lZ2DmeAs51gulNd6Fe9Mmyqmpc5s/4cVxeW1Zqy2NEW3ft6jR/VbxcdLsW/rKOk6FoPySou3vI5AxTSptFq1Ki+k//c/YPjw8jGMhDiroiLAze3/s3fm4VGVZ///TEgCCWFfBIIoIFRQqQRc+6pdEBBwFBRp3MFd1KZLqCuLSgtobVS0WkCtFQc3QFSwuLRWXvtaSfRX2UQtomxKAIEQAlnO748nh8xMzsycM5kzZ57M/bmuuWZy5syZ73nyzMk399zPfUNWVuIXINrFMNQ/u/FGpiOZW6fVPMKPkw4LIgVBcA8x00KzJjcX7r8f7rkHrr0WnnvOa0WC4A1mZDo7Wz32gro6lWbiRo51+H7RXpcM8ywGXRDSBzHTaUR5ebnXEhyRSL0nnKByqT/7DK6/Hg4cSNihQ0jnMU4Gojd+7JppNzXbrboRvH9syi07IMZ/PLdJnTlhl1Sax3YQve6im95kIGY6jZg0aZLXEhyRaL2ZmXDffap8nt8Pf/tbQg8PyBi7jeiNH7tm2k3N8aR5xDbf1noTsQCxKeY78nukzpywSyrNYzuIXnfRTW8yEDOdRkxPpRUtNnBL77BhsGwZvPUWXHVVYit+yBi7i+iNH7tm2qnmujr75rSuzpmZtlelY7qrEefEL0KcnugDuk4qzWM7iF530U1vMhAznUYUFBR4LcERbupt3RoefBB+8Qu49FJYvDgxx5UxdhfRGz92zbRTzbW19vOga2vVAsh4I9PWxtZar9W+5vsmPtXECc7nhHmN+uyzRGuxRyrNYzuIXnfRTW8yEDMtpDUFBariR2mpilLv3u21IkFwB7cWIJpmOtH7QmMjG8nYOjW8do/blKh0Ik34q6/C55+rOvqCIKQeYqaFtKdlS5g5E26+GcaPh+XLvVYkCInHbTNtx3jW1qq1C251TEyUgV29Gr79tmnHSKSZrqyE7dsTdzxBEBKLmOkmUlRUhN/vJxAIeC0lJgsWLPBagiOSrfe001SU+p13VMWPffucH0PG2F1Eb/zU1CjTG8tMO9XsNDLtxEyH72dtrBdYPh8t3zqWQR8xAp591pbEOHA2vpddpu4ffljde1EjPJXmsR1Er7voojcQCOD3+ykqKnL9vcRMN5GSkhKWLVtGYWGh11JiUlZW5rUER3ihNycH/vAHuPxyuPBCZaydIGPsLqK3afh8sc20U81ummmI3mBFPS5rcppHcnE2vs8/H/rzSSclUIpNUm0ex0L0uosuegsLC1m2bBklJSWuv5eY6TTiscce81qCI7zUe/bZquLHK6/AbbfZr0stY+wuorfpxDLTTjUnMzJt/by13kRU4XCnNF7T5sSGDU16eVyk4jyOhuh1F930JgMx04IQgbw8ePxxGDMGzj8fPvzQa0WCED+mMfR6AWJWliqRZ5donQ3tNmsJ39+J0Y7XlLsZ/d60yb1jC4LgHDHTghCD4cNVhPrRR+Hee73JWRSEROGlma6pUQt+a2vt7W+vzrS7pEbXxFD69PFagSAIwYiZFgQbdOgAf/2r+iM2Zow3X7UKQlMwI6xeR6azs+2baYheZzreyLTd7eb7bdli7/iCIKQnYqbTCL/f77UER6SaXp9PLUycPx/uvBN+9zuorg7dJ9U0x0L0uksq6o1lpp1qdtNM28uZbtCbqG6FhhH63kcfnZjjKhIzJxLfmTEyqTiPoyF63UU3vclAzHQaccstt3gtwRGpqrdnT5X20bs3jBoFn3zS8Fyqao6E6HWXVNSblRXdTDvVnMzIdDjK8N7SyPymNqk3J2KRivM4GqLXXXTTmwzETAfxpz/9iYKCArKzs5kxY4bXchLO8OHDvZbgiFTW6/NBYSEsXAizZsH996t80FTWbIXodZdU1JuREd14OtXs1Ey3bGl/3YE9g2xfr900D59P3dwx6ImbE8mKTqfiPI6G6HUX3fQmAzHTQfTo0YN7772XCy+8EF8yv0MTtKVrVwgE4NhjJZdaSE+SHZluXGfa3ai0Dn8Kqqpg82avVQhC+pLptYBU4oILLgDg1VdfxdDnO0PBY8xc6h//GG65BX7yE7j1VhUBFISUYOhQFqzdAT3Vj0/t4shjS7p1Uz21beC2mY5GtMu0F23I3T6mFTt3wtKlqmur/NkSBG+QP/dpxNKlS72W4Ajd9PbsCVddtZTsbLjgAvjqK68VxUa3MRa9cbJ1K52rtsLWrY0eh9+Wbt0KO3bYPrT3CxCX2i6hF4/ZTLxBTeyc6NpVGWo3SZl5bBPR6y666U0GYqbTiEAg4LUER+imF2DRogA33QR//CPcfDMsWJDa0SLdxlj0xofRvoN60KUL5OdT3iof8q1vgZwcFZm2iRel8b79Vi2iVJ8t98a4KSkekT/3idd7110JP2QIqTKP7SJ63UU3vclAzHQa8cILL3gtwRG66YUGzccdB6+9Brt2wcUXq6BfKqLbGIve+Kj883PqwZtvwpYtTDp3iyqebHF7obLSdooHeBOZPu44ePZZ86cXHP/DGi0P2/zZvX+CU2NOOCFV5rFdRK+76KY3GWhrpisqKpgyZQrDhw+nS5cuZGRkRKzAUVFRQVFREfn5+eTk5DB48OCYk0EWIApNpUULmDJFVfq45hp4+unUjlILzZdDh9w7drCZjjW/nRhvE6tLcWWlOpbTz1OimrwIgiAEo62ZLi8vZ968eVRXVzN27FggsgEeN24czz77LNOnT+fNN9/klFNOobCwsNFXFbW1tVRVVVFTU0N1dTVVVVXU1dW5fi5C82bAAHjjDSgvh/HjYft2rxUJ6UZVlXvHrq5WtaszMiDW5dKpmbaXM20dSW5KPMQsjdeU43gRj/n22+S/pyAIGlfzOPbYY9mzZw8Au3btYv78+Zb7LV++nLfffptAIMCECRMAOOecc9i8eTPFxcVMmDCBjPqyC/fddx/33nvvkdfOnDmTZ555hiuvvNLlsxGaOy1aQHExrFkDV12lItX101EQXMdNM11To8x0ixaxzXI8ZtpOaTyT4H2dVPpwIwL9y18m/pix6NZNoumC4AXaRqaDiVbGbsmSJbRp04bx48eHbJ84cSLbtm3jww8/PLJt+vTp1NXVhdyak5GeOHGi1xIcoZteiK35xBNVlHrdOrjiCti9O0nCIqDbGIve+KhykObhVHN1NWRmNpjpaNTWqn3tRm3DzbT5ODRdY6JmaRvuzokDBxJ/zFSZx3YRve6im95k0CzMdDTWrFnDgAEDjkSfTU466SQA1q5d26Tjjxo1Cr/fH3I744wzGpWOWblypWU/+8mTJ7NgwYKQbWVlZfj9fsrLy0O2T5s2jdmzZ4ds+/rrr/H7/WwI6xby6KOPUlxcHLLtnHPOwe/3s2rVqpDtgUDA8sMxYcIET89j+PDhludRWVmZsudRUFAQ8/eRlQUzZsD111fygx/4efBB785j+PDhTZ5Xyfx9DB8+3LXPhxvnYXYKS+bn3Oo8qg6qGhITp08HQk1lyHls387wfftYGQjYnlfr15cxb56fmprykDQPq/PYtu1rnnzSz7599s5j+XI/Bw6E/j6++CKAYQT/PtQYv/nmBPbuDS/ZZT2v5s2bDDSch2Go3wf4gdDfx+efTwNCzwO+pq7OD4R3aXoUKA7bVll/XPM8zO5xAayN9QQal89bWX+McELPA6CwMPHzavjw4Sl93Q0/D/Nz5/X1yu55mHq9vl7ZPY/gDoipdt0N1F+7/H4/vXv35uSTT6aoqKjRcRKO0QzYuXOn4fP5jBkzZjR6rl+/fsZ5553XaPu2bdsMn89nzJo1K673LC0tNQCjtLQ0rtcLgmEYRmWlYRQVGcaNNxrG/v1eqxGaK6v+vNbY1W2gYaxdaxiGYYwZE2HH0lKVfuzgurZ0qWHMn28Yl11mGN9/H33fV181jD//2TDOP9/esW+80TAGD274eehQw7j+esPIyDCMJ55Q7weG8c03hnHLLYYxaFDDvhs3qucMQ91/9ZVhfPGFevzyy6HPffttw2MwjI4dDWPOHMO46CLDGDGiYXvwzeez3p4Kt5077Y2vIKQDyfBrzT4yLQipTE6Oqkk9fjz4/bBypdeKhObIzi4DWTpzLQwcCNhbLGgXM80jK0s9jkaiqnmYOEnb8PmcVfGIVR4vlQs+/eQnXisQhPSi2ZvpTp06sWvXrkbbd9cnq3bq1CnZkgShET/9KSxbBitWwLXXwvffe61IaE5UVkJubsPPOTlw8GBijm0uQHTDTMfqYhit1J0To62raY7EmjXwzjteqxCE9KHZm+lBgwaxfv36RiXuPv30UwBOPPFEL2R5QnhOUqqjm15omua8PBWlnjgRxo1LTpRatzEWvfFx8GComc7NVQbbCqeK44lMOzGosfe1pzjRiw/jj+wnZ04MGwZPPpmYY6XKPLaL6HUX3fQmg2ZvpseOHUtFRQUvv/xyyPZnnnmG/Px8TjvttCYdv6ioCL/fr0V7zTlz5ngtwRG66YXEaP7Rj1T3xNdfVy3JKyoSICwCuo2x6I2PAwcam+lIkWmnis0605mZKkodDXfqTM+x3NfKhEc6ntVr3YtIJ29O3Hhj6M/hv/P//tfecVJlHttF9LqLLnrNxYjJWICobZ1pgBUrVnDgwAH2798PqMocpmkePXo0OTk5jBw5knPPPZebbrqJffv20bdvXwKBACtXrmThwoVN7nRYUlJCQUFBk88lGSxatMhrCY7QTS8kTnPr1vDII/DuuyqX+v774cwzE3LoEHQbY9EbH/v3Q9u2DT9Hi0w7Vew0zaNlS+e5zuGPzZJ56jiLbB0vNcrigfMRThy5uaHj0LevvXFJlXlsF9HrLrroLSwspLCwkLKyMoYMGeLqe2ltpm+++WY2b94MqO6HL730Ei+99BI+n49NmzbRq1cvABYvXsxdd93F1KlT2b17NwMGDGDRokVccsklXspPOrnBoSkN0E0vJF7zT38KQ4bArbeqHMg773S+gCsauo2x6I2PvXtDzXROTmQz7VSxmwsQ7UWmEzPGdqPWTSc15oST80qVeWwX0esuuulNBlqneWzatOlIc5Xa2tqQx6aRBmjdujUlJSVs27aNqqoqPv7447Qz0oK+tGsHf/kLHHMMjB6tGr4IghP27bMZmW7VSlX8aNXK9rHNNI9kVvOwbt4SvQNitGoehqEW7QX/bNV9UUeGDm28bedOd7tiCkK6oXVkWhDSBZ8PrrxSRap//Wvo319FqXNyvFYm6ICVmbbMmR44EBw2sqqpcW8BYrToqdOIcaz9Bw1q+nukIqWl8M03KnUM4Lvv4KijYMsWb3UJQnNC68i04IzwzkOpjm56wX3NPXvCCy+o1I/Ro+GTT5p2PN3GWPTGR8cd62h92glHvtaIlubhVLObCxCtaFwar9jS9EYz7ImtJuKU5M+JXr1g2jT12Gz5/sAD9l+fKvPYLqLXXXTTmwzETKcRwakvOqCbXkie5gsvhBdfhHvvhccei79Ml25jLHrjI6u2Ct+6dUe+24+2ANGp5njSPOxGfMNTM6wXINrXG6kudXIj0N7MifDf9yuvqPvPP4/92lSZx3YRve6im95kIGY6jbj11lu9luAI3fRCcjV37gwvv6yM9OjRoTmfdtFtjEVvYohWGs+p5mnTVIq1GznT4Wba2gTfauufScOw/0+nu6XxvJkTTz0V+rOZ5tG/f+zXpuo8joTodRfd9CYDyZluIkVFRbRv3/5ICRZBSCYZGarSx7hx8JvfwHHHwd13q/JjgmASbkJzcxObM9uihT0zbeZXO8mZtrOv1QJEJ/s1h9xoQRBCCQQCBAIBvk9CS2GJTDeRkpISli1bJkZa8JT8fAgE4OSTVZT644+9ViSkMtFypuPBrDXtRmk8qzrTwc+DvYhzcJQ7Ua3GdeTNNxtvmzkz+ToEwW0KCwtZtmwZJSUlrr+XmOk0YsOGDV5LcIRuesF7zRddpEz1nDkwY0Zsc+O1XqeIXudYmcFoaR5ONXfpAgMGuNcB0cpAh5bG22C7aUuk/cI7BUbbt+l4OycmTWq87e674cMPI78mFeaxE0Svu+imNxmImU4jpkyZ4rUER+imF1JDc5cu8PzzKhdy9Ojolc5SQa8TRK9zqqoal42OtACxvBwKC51pPv10+wsQq6shO9u+UY2U5hEaYZ5iOxfa3C/8mK+/Hvk1ic+d9n5OWDFuXOTnUmEeO0H0uotuepOBmOk0Yu7cuV5LcIRueiF1NPt8UFiomr1Mnw6zZzeUxAomVfTaRfQ6Z98+yMsL3RYpzeOzJet44JMNcXUGsmOmDx9W+9nFykw3bswy1/YCRKfR5miNXuLH+znhlLlz59Kundcq7JMKnzsniF79ETOdRuhWzkY3vZB6mrt3VyX0jjoKxoyBjRtDn081vbEQvc7Ztw9qu3ZXZTe6dwciR6azaqsYxudxtcdzEplO1AJEZaJ72Y5Mp0b+s/dzIhaDB6v7w4fVmPXq1Yt9+7zV5IRU+Nw5QfTqj5hpQWjm+Hxw9dXw5z/DlClQUhJ/XWpBP/buBV+P7uoriiAzbZUzHU9Kg2lQkxWZDsb8tsVNM50aBjy5mM2gevaEt97yVosg6IDt0nilpaX44rjSDhgwgBzpeSwInnP00bBkCTz5JFxwAcydC8cc6ZddggAAIABJREFU47UqwW3CW4lD5DSPpuQH21mAGByZrqtTpR2jYS8yre5jmd546ky7V2taD3buhP37vVYhCKmP7cj0KaecwtChQx3dTjnlFNavX++mfsEBs2fP9lqCI3TTC6mv2edTlQv++Ee46SYYP362VpG3VB/fcFJBr5WZbtHCOofeMMCpYtNwOolMR3r/cOrqIlf/MA05NMzhaGX0IFWizN7PiWh06KDug6t+BM/jHTuSLCgOUuFz5wTRqz+Omrbcfffd9OnTx9a+dXV1XHvttXGJEtyhMpGFZZOAbnpBH83HHQevvQYjRlQyfjw88gj06OG1qtjoMr4mqaDXykxD5KhrvIqd5ExnZiozHSvlIzx6Hd4NUZnpyoSnebhbGs/7ORENs7/FkiXq/qOP4O9/b9DcvXuq/FMSmVT43DlB9OqPIzM9ZswYTj31VFv71tTUpIWZ1qkD4owZM7yW4Ajd9IJemlu0gLffnsHatTBxIlxyiYpGpfJX2zqNL6SG3n37oFs3+/vHq9iumTYj07FSQkAZ7mipIMpEz3Ctmoc7eD8nrNi2LfRn01T/5S+wY0dqao5EKnzunCB63SGZHRBtm+nFixfTv39/+wfOzGTx4sX07ds3LmG6UFJSQkFBgdcyBCFuTjgBli+Hhx+GCy+ERx8FWazdfNi3D8uyZlbGcsDvLlcPRo5UIWQbPLUL6AlDq+HEKuCx+ie6dYPVq0P2PXxYHTY7O7bxhsZpHuH/6AXnTNvB6cJbd0rj6YUOaR2CYIUZ5CwrK2PIkCGuvpdtM33hhRc6Png8rxEEIfm0aAG/+pUqn3fDDSpKffXVqR2lFuyxd691mocVmRV71IOdO20fvzPAVshG3TBLqG3frspBBPH4Lmh5HPxxD7RdifWqnSATHmmRojkv6+qCc6dj49QYf/cdvPees9ekA/ffr7omCoKgcJTmIWjM0KGUb91KZye9fD2mvLZWK72gn+Zwvf2B5UDFB7B7MrRvDy3iKaBpEZVMBOXl5XTu3Dnhx3WLVNC7axd0yj0Ia/8LffqoUh5Y/6N0uHM+G3fDoHz7c7h8F3TuBNU1cOAAtK/crtxtXR1s3Rqyr2m8O4Ct1OFoFT/MnOmMjHLq6mKPcTxpHps3O9vfHuXUj0TK0bp1pGdCNUdrPZ4KpMLnzgmiV3/iMtPvvPMOu3fvZvz48QB8++23XH311Xz88cece+65zJs3j1bh/WsFb9mxg0k7drDMax0OmARa6QX9NFvp9QFtzB8sahHbwqVC1pMmTWLZMn1GOBX07toFHb9dD6cNgdJSqE9LM81lsKle/afV/Oxnfowt9jVP8sOyZbBlk6oS88gHQyPmBpjG+/u9qitjppVnD0rwjrYA0XzeMCZhGI31NqWah7tl8VL3KhF5XZnSbKbKpnrqRyp87pwgevUnLjM9bdo0hg0bdsRMT5kyhVWrVjFs2DBeeeUV+vXrx9SpUxMqVGgi3box3UxY1ATd9IJ+mmPpNVBpAoYB7ds5MBhOVrw5YPr06a4c1y1SQW9traqeEY7ZuCU3N/yZ6XG9T8uWcOgQUb+RMI33zGK49lr4wQ+iH9NM44j2fMuW0xNeZ9pdpnstIA6mA/DKK+onF750Siip8LlzgujVn7jM9MaNG/ntb38LQHV1NUuWLGHWrFlMnjyZBx98kKeeekrMdKqxejW6LZPUTS/opzmWXh/QHnjjDXjoIdWR+uyzkyAsArot9k1lvW3bqsWJwWZaGVL7moPN7hEzHQVz3+xstRjRDlZmOjhnumXLgqgmeffuhsdOS+O5E51O3TkRGaX50089lmGTVP7cWSF69SeuduL79u2jQ31l99LSUioqKrjgggsA1dxlszuJZoIgeMTo0Soq9eKLcN11KnVA0Ju2bZve3S74iw07ZtrEThm9mhpVC90Ks8pGXZ2Kukcz08GtEdK9MocgCO4Ql5nu2rUrn332GaDyp4855hh61q/a3r9/P1mxKvELgqAd7durFuTXXgsTJjQ0dRBSl2h1ms3IdDBOI7GHDikTDc7MtJ3SeHv3xj5OuJm20m8eJzjNQ6rUCIKQSOIy0yNHjuTOO+/k17/+NX/4wx9CSuB99tlnHHvssYnSJySQBQsWeC3BEbrpBf00x6P3tNNU2sfq1XDVVaFfo7tNOoxvIikvh0iL7tu0aWymVeTWvuaqKjDXmmdm2mvEAioyHSvNI5LhDd5eVwdVVQssI9PRFiA6WYiYePSawwq9NHv9uXOK6NWfuMz0zJkzGTx4MPPmzaOgoIC7gwpOPv/885x55pkJEygkjrKyMq8lOEI3vaCf5nj1tmwJM2fC5Mkwfnzkr+MTTbqMb6LYsSPyWlCryLTCvubgyLQd42maWMvI9Lp1qoPQunW237+uDmpqyhzlQnuPXnNY0Vjz44+ra4C5GHHHjlQZX+8/d04RvfoT1wLELl268Oabb1o+9+6775JTX8dUSC0ee+yx2DulELrpBf00N1XvqaeqKPXUqSrt46GHVDqIW6Tb+DaV+My0fc3BZtoJlpHpqiplpKuqYr4+OGf6qKMes4xMhxu74DrT3qZ56DWHFY01T56s7nfsgKFDoXt3+Pe/4ZRTkizNAq8/d04RvfrT5KYtO3fu5ODB0GK0e/fupZf0IxaEtKBVK5gzBz74AMaNg+JiOO88r1UJEGSmBwyANWtCVuO1bQtfftm04zs106aJtbMA0Q52FiCG7y+4h900H0FobsRdzeOaa64hNzeXo446imOPPTbk1rt370TrFAQhxTnzTHjjgXUMvvwEpo1fZ2sBmeAu335bb6ZzclQKRdC3homo5hGcM+0EOwsQ7USPg810rBSDAQOclcYT7LF3Lzz6qHos4yakK3FFpouKiggEAlxzzTWcdNJJtIznez5BEJodOb4qcnav44IRVYwdC3fcAeee67Wq9GXr1shpHpEXINrHiZkOTrOwswAx2nFMws10pMol4a/9/nt77yVVP2Lz17+qG0hkWkhf4opML1++nN///vfMnTuXG264gauvvrrRTUg9/H6/1xIcoZte0E+zW3oLCtSixNdeg5tuanoE1ETG1xlbt0J91dJGRM6Ztq+5okKZcjvU1CgTDfYi09Ewc6Zra+Grr/xHzHQs82ua6UmT7L9P4tFrDivsaT7nHJdl2MTrz51TRK/+xGWmq6qqGDRoUKK1CC5zyy23eC3BEbrpBf00u6m3dWt45BFVk/qCC+Ddd5t+TBlfZxw6FDlyHNlM29dcUQF5efb2DY5iJ6I0nrkAsVu3Wyw7FlpF2e3mTLubrqDXHFbopdnrz51TRK/+xJXmcd555/H+++/z05/+NNF6BBcZPny41xIcoZte0E9zMvT++Mdqtf9vfwuLF8OsWfYNWDgyvhEYOlStNgzj6V1AcGS6W7cjtczatGn8jYEykfY1V1QoU26HcDPdKDJ9+eXqfuRIyM6mbR18A7R6q+EcXvsOstfB1MPQ/m5o0QL+uR8yi7ux4qerbUem7eKOqdZrDiv00izXCXfRTW8yiMtM33PPPVx00UXk5eXh9/vp1KlTo306duzYZHGCIDQP8vLgscfg7bfh/PPhnntA/hdPIFu3WprpzgBbrV+SmanSJIJxWu1i/37o0aPh52jms6qqYf2jZZrHnj3qfudOQH1t2hOgiiPncBRANXQAqM97zgEO7q6jrs5+mofJtm3R9xcEQbBDXGb6xBNPBKC4uJji4uJGz/t8PmrDr9LNlKKiItq3b09hYSGFhYVeyxGElGbYMNVB8Y474KWXVEk9uzm3QhQ6dIAdO9jbqgu1Gdl07AA1tcrsdgiu+x22GtGqFrMTwnOmzVxmK1N78GBDZDo726KcdH6+CjXXU1cH27ar13Suj9d8+5167YEDqp55plFNq73fcbh1e0c50ybbt0ffXxYgCoK+BAIBAoEA39tdcdwE4jLTU6dOjfq8L42uQCUlJRQUFHgtwxZLly4Naf2e6uimF/TTnHC9YV/VW9EGmIvK5933NGS3sV+reOnBg1zYu3dD27UUJ2nz4bnnYMgQZv7oTT7JKGDlSvjX+/C//wu3327/MMpsLgXsaQ7PmW7ZMnKednCaR6tWar8Qwn6ne3bB0Z3h/HNh2TK17fxT4eST4dln4cH7YGBVGfuKh9Dypucwvohtfp3mTLvzp8z++KYO9jXPmwfXXeeumlik/XXYZXTRawY5y8rKGDJkiKvvFZeZnj59eoJlCMkgEAho8QEw0U0v6Kc54XrDvqqPRkugC0C4qYpCALhQow6ryZ4PGRkN5ck2b4Zjj3X2emU2A8RrpnNzQyPQwYSb6ViNDq2i5OHb6uqU2qtsLkB0Gnn/5htn+9vD/vimDvY1X3+9Sp+ZNs1dRdFI++uwy+imNxk4NtOVlZX069ePJ554gvPPP98NTYJLvPDCC15LcIRuekE/zQnXG/ZVvV2qqmB/BbRrGzGgDcALELlwcgqS7PnQogXU1v9zsmkT2FknFGxCldm0r9nKTFdWqqyTcILNdE5OfGY6/PlTHrmcYcDBOSM5vTZb/SNRv1jx2Bq1gDGYTteGbus6qvE+AL790HYmXBq7s3lc7GAop6DHtysKZ/N4+nRvzXTaX4ddRje9ycCxmc7NzeXgwYO0bt3aDT2CIOhMnOkXrYCDe+C6IlUXeepUZ22qBUVGRkMqw6ZNId3DLWndWplf83IezwLEYDOdk6OOZ0VwxNpOZDoaR9qSH1DfhOTs38mR7yvqFytmEVrIBIDdYdu+s9gHwAD2qZQkN8ig+fc1X71aFZkRhHQgrjSPn/zkJ7zzzjtSGk8QhITRoQP85S/wyiswerRanKjJcoSUIdhM79wJnTtH379LF7Wfaaab2gHRjExH2tfM0Ik3zSOcgx3yOXioxZEc7Joa6NpFPVddo9qpB9OxI+ze3fBz167w3XeNj+vzQds2sNeyDnf8ZFJDJ3axm+Zf7WrFCjHTQvoQl5m+++67GTduHNnZ2Vx00UV079690aJDKY0nCEI8XHQRnHWWqkvdpg3ce6+q3CDEpkULVe7OqomJFZ07KzNt5lbHU1c5+D3MnGkrEpEzDfDVVw3vueLe1axYocosvveeeu6f/6zf73Po3z/0tS/8STURMlm93Nrw5bWGe+5Sc1CIj/vuUyUwBSEdiKsD4pAhQ9i8eTMzZsxg0KBBdOnShc6dOx+5denSJdE6hQQwceJEryU4Qje9oJ/mVNXbtSs8/TRcfLEy10uXAuvWMbF9e1i3zmt5tkna+LZqBQMHUpfditpaFZG1k1repQuUlzf8rKLa8WuOluaRKDP91lsNz9fVwYcfTjwSja+rgxkzIr/eyT8L7nVBTM3PXHSca25Ku/imkqrXtUiIXv2R0nj17Ny5k6uvvpr33nuP/Px8HnvsMYYNG+a1rISiW9ci3fSCfppTXe/ZZ8Py5SqH+pOnqhi+d2/Tkm2TTNLGd+BAWLuW726G2i9g7Vo44YTYLzPTPEycdkAMJ1aah5l2kplp32xFq11dVwc9ew4/YqZra6MvfnOSE+7en7HU/sxZo5fmVL+uhSN69UdK49UzefJkevToQXl5OW+99RaXXHIJX3zxRbNKV9GtqYxuekE/zTrobdkSZs+G//c0/PA1+L//g9M1yaVO9viaaR5r1tgz0507Q1lZw8/KTNvXHG44o6V5HDignrd6nRVmZNjMA7cqElNXB/37F9qOIid6v/hI/c9cY/TSrMN1LRjRqz9xpXk0NyoqKnj11VeZMWMGrVq14vzzz+eHP/whr776qtfSBEGo54c/VPf//CdMnqzKsgmhZGQoI+gkMh2c5mGayKoqy+7kjQg3ndHSPPbvd9btMtxMB2Oa8bq6hrboZgTbid5E7SsIQnoTV2QaYOPGjTz55JNs2LCBg0GhCMMw8Pl8vPvuuwkRmAw+//xz8vLy6NGjx5FtJ510EmvXrvVQlSAIVkyZAm/vBr9fLU78n//xWlHqUF2tzOW2bRB0OYuIdZoHPPkk3H9/9N47dXXWkelIr4nHTN9wA+zaZZ2eYeZMZ2crM233mHbQMFNREAQPiSsyvWbNGgYPHszrr7/OihUr2LNnDxs3buQf//gHX375JYZm/9JXVFTQtm3bkG1t27alopmFvlatWuW1BEfophf006yd3vr7YcNgyRLVRds0XKlIsse3uhqystRjO4awY8fQsVOmdRW1tbFT0ysrG9I2TKLlTMdjpn0+68i0iVpsuYqamsaRaas/Q07bibuDXp85hV6atbuuiV7tictM33nnnYwYMYI1a9YAMH/+fLZs2cJrr73GoUOHmDlzZkJFuk1eXh779oUWFN27dy9tnFz5NWDOnDleS3CEbnpBP83a6Q163K4dPPEETJwIhYWwYEHqfTWf7PGtrlbpGfn59vbPympoPw7m+M2xlTKxezd06hS6LS8vcvpNPGY6IyO2mS4tnXMkMu00zcObCLRenzmFXpq1u66JXu2Jy0yXlZVx9dVXk5GhXm5GokePHs1vfvMb7rjjjsQptKCiooIpU6YwfPhwunTpQkZGBjPMekgW+xYVFZGfn09OTg6DBw9u1AqzX79+VFRUsG3btiPbPv30U06wk3SoEYsWLfJagiN00wv6adZK7+WXswhg5EjVJrH+dvrFPfnbup5c8uue7G7dk9oePUOet3VzqbtEsse3uhq+/hrOPDO+16tL+SJbZrq8vLGZbtMG9kVodHLgQENzGDuYaSQtWoQafmgwwTU18POfL2r0fCTCzynSOdo5/3gYwDo+4nMGoE95R4VG1wk0u64hepsDceVM79mzhw4dOtCiRQuysrLYs2fPkeeGDBkS0dgmivLycubNm8fJJ5/M2LFjmT9/fsRyfOPGjWP16tXMnj2b/v37s3DhQgoLC6mrqzuyIjUvL48LLriAadOm8eijj/LWW2/xn//8B7/f7+p5JJvc8O9kUxzd9IJ+mrXSu2cPuWCZlOsjqPVzhGoSXpDs8TXLzf3kJ/G9XhlIe5p37WrcYbFtWxWBjnTsDAfhGzPNw1xgaPV8bS20bZvLwYP2mtR4Xc2jFVUMZQOt0Ke8o0Kj6wSaXdcQvc2BuMx0fn4+39b3ae3bty/vvfce5557LqAiunl5eYlTaMGxxx57xMDv2rWL+fPnW+63fPly3n77bQKBABPq216dc845bN68meLiYiZMmHAkuv74449z1VVX0alTJ3r27MmLL77YrMriCYL25Odb10cLwzBgfwUcPqxSQbLsXOXsdDhJZdatg/Hj6XrUS1x33UB69bL/0qyshlxr00TajUyHm+k2bSKbaacEm+lINalralS0e/9+e2baSc60LEIUBMEucZnpH/3oR/zf//0fF198MZdffjlTp05l+/btZGdn88wzz3D55ZcnWmdEoi12XLJkCW3atGH8+PEh2ydOnMill17Khx9+yBlnnAFA586deeONN1zVKghCE1i92tZuPqAtsHkz3Ho79O6t2hrn5LiqzluqqmDdOjI7VTH3z85e2qmTijJ36xZajs6OmT7++NBtrVsnrmShaWiD87rDNd1+u8qbr6lpXF3ETgfEaIY51fLvBUFIXeLKmb7rrruOpEBMmTKFm2++mSVLlvDSSy8xYcIEHnzwwYSKjJc1a9YwYMCAI9Fnk5NOOgkgIaXvRo0ahd/vD7mdccYZLF26NGS/lStXWqaNTJ48mQULFoRsKysrw+/3Ux5cABaYNm0as2fPDtn29ddf4/f72bBhQ8j2Rx99lOLi4pBtRUVF+P3+RitxA4GAZXvQCRMmeHoexcXFludRWVmZsudxww032P59pMJ5FBcXN3leJfM8iouLbf8+jjkGZs/+mnfe8fPjH28guFpnss7DfI9kfs6dnseqVRN44QV1HipyW8yGDSs5dCj678PMmQ4+j2BzGus8MjMbTHKk83jtNT/ffbcqLCc6wKFDDeexZEkxNTXw/vsT2LMn9PcBK4GG82gwyJOBBXzwQfC+ZfX7hv4+YBowO2zb1/X7bgjb/ihQHLatsn5f9ftoeDaAdZvuCUD082hAnUcobpxHMeHn0UDk8/Dq74c5l7y+Xtk9D1Njql53w88jWEuq/f0IBAJHvFjv3r05+eSTKSoqanSchGNozs6dOw2fz2fMmDGj0XP9+vUzzjvvvEbbt23bZvh8PmPWrFlxv29paakBGKWlpXEfI9k88sgjXktwhG56DUM/zemid98+w7jtNsO45hrDKC8PemLtWsMYOFDdu0DSxre01DDA+MVZzq9Hf/iDYbz7rnr89NOGAY8Yjz1mGFlZ0V83ebJhbN7cePuYMdb7n39+6M+XXmoY+/dHPv7nnxtGUZG6ffGF2jZkiGGAYbRpYxgPPaQeX3HFI8avfmUYP/+5YZx5ptpmGIaxfr16HHx74onQn1u0aLwPGEZurmH8/vfWzzXlNphS4xEwBlOa8GO7e3skrtd5Rbpc17xCN73J8GvSATGNuPXWW72W4Ajd9IJ+mtNFb5s28PDDcO21MGECBAL1Ucr69IiYRZXjRIfxDW7coiK3t0YtR2dilTMdCasGL61aqWGPlE4RK2f6s8/U/Zgxt7JvH3z7bew8Z7vNXcz3d4PUnxFW6KVah89dMKJXf+LugFhTU8OLL77IP/7xD3bt2kWnTp348Y9/zCWXXEJmZtyHTSidOnVil0U3h927dx95XhCE9OH002HFCnjgAbjoIpg7CWw0CtSGeNZRdukCX36pHjtZgHjgQOOmLQCvv67MeZcuDdsqKlQN6mBMM52VBZ9+CgMGhD5vlTMdbJaffFLdZ2WBuQY9VjfMSG3JrTh0KPqx4uE51HqiFYykmuzEv4GL7KAbp2Bv3YIgpBtxud7y8nJGjBjBxx9/TGZmJh07djxSVePBBx9k5cqVdLYbsnCRQYMGEQgEqKurC8mb/vTTTwE48cQTvZImCIJHZGXBnXfC55/DnGugBFX5Qy9rY01BgfPXdO8O77+vHgcvQLRb+cKKPXtCzfT+/apsXjCmma6ttV60GByZjrQAERo6PkLsBYjhkelo/zD82eFCTju0R1WhOooofdpTlAyaMCEEoZkTl5n+5S9/ycaNG1m4cCHjx48nMzPzSKT6hhtuoKioiOeeey7RWh0zduxY5s2bx8svv8wll1xyZPszzzxDfn4+p512WpPfo6ioiPbt21NYWHikbnWqsmHDBo4PX36fwuimF/TTnM56+/WDP/4RGAq33gqj74Hzz09sSbRkj+8ppzh/Tc+esHWreqwM9AZ8vtiaI43TL38JB8NqfVt1PzTNdCSipXkEv/eOHRuA46NqMnHyD4LdRjBO2EY+G6njOLJi75wiZFHNHr6jjvZeS7FNOl/XkoEuegOBAIFAgO+//97194rLTL/22mvcd999IeYxMzOTSy+9lO+++45p06YlTGAkVqxYwYEDB9hfX9R07dq1vPzyy4DqxJiTk8PIkSM599xzuemmm9i3bx99+/YlEAiwcuVKFi5cGLHRixNKSkooiCcc5AFTpkxh2bJlXsuwjW56QT/N6a7XvASUlMCslfCXv8CsWcpoJ4Jkje/ult15e8A0LhnY3fFrO3ZUpfHAjNROISMjuuZoEd1u3Rp3QYxkpsNNd/h7hEemrXjqqSmA0hurKYyTnGk3UGkSfky9OjCYMnoyhC14HyCzS7pf19xGF71mkLOsrIwhQ4a4+l5xmWnDMCKmSJxwwglRaz8niptvvpnNmzcD4PP5eOmll3jppZfw+Xxs2rSJXvVdCxYvXsxdd93F1KlT2b17NwMGDGDRokUhkep0Ye7cuV5LcIRuekE/zaJXkZMDM2bAf/+rahcfdxzccUdjA+iUZI3v83/vTvf7poNzL90oNcLnmxvTlFZWWudLg3UXRCsz3bq1yruOhFXOtBW//vVcLr5YPU7UAkR3/4Tp9ZmrohW30Y/baOW1FNvIdc1ddNObDOIy0z/72c94++23GTZsWKPn3n77bX4Sby9bB2zatMnWfq1bt6akpISSkhKXFaU+vZy0RUsBdNML+mlOe71mg6mRIyE7mz7Ai0DVu7D/IfDlKNMX0aN16xa1mUyyxvf11+HVV5t+HMOAjIxeMc10tEoebdrYi0wHm+lI+c0tWsSOTB93XMMYJ7Kah3vo9Zlbz0DOZaPXMhyR9tc1l9FNbzKwbabNChgAU6dOZezYsdTU1HDZZZfRrVs3tm/fzsKFC1myZAmLFy92RawgCEJC2aMWhB2pDVdPq/obh4Bo6XZNWaWXID79FPr2hZYt4z+GaYDNaHAsU7prV2Qz3bYt7NgRum33bpVOEkxeXsPCQyszXV0N2dnR24mfcUbjLowmkQz6FVfAX/9q/ZpYrxcEQbDCtpm2qs7x0EMP8dBDDzXaPmTIEGpTIwQgCIIQmfx8Ff6MQp2hIqs11dC2HWRlotzdd99Be+8XZf3udzBzZtOO0bMnbNmi/jfIyIhtpnfuVN0PrWjbFjaGBTLLy1XqTDCtWyuTHYnDh1WKR7TIdOvW6nm71NaG7p/IxaaCIKQvti9DU6dOtX3QRCzsExLP7Nmz+e1vf+u1DNvophf005z2eqOkaJhkAO1Q+dST71T+e7q/jDY/HgIxqha5Pb5lZcrP9+nTtOP06aPOzzCgtnY2hhFd8zffwNFHWz9nleZhlRaSlwdffx35PczItFXOdHDU+IEHZgP2xjjcTHsTfbavN3XQS3PaX9dcRje9ycC2mZ4+fbqLMoRkUFlZ6bUER+imF/TTLHrt06cPLFoE774LxcXwBLHrU7utd+ZM1dmxqfTvD598YtZsroyZvbJ5Mwwdav2c1QJEKzMdawGincg0wMGDDWMcyxyHm+louGe09frMKfTSLNc1d9FNbzKQduJNpKioCL/fTyAQ8FpKTGbMmOG1BEfophf00yx6nfPTn4K5mP222+CRRyKXeHNT7+LFygT37Nn0Y/Xvr1IzDAOys2dmkQzcAAAgAElEQVTENJJffQXHHGP9nFVkevdu6NAhdFvr1tbNWkyscqbDdRmGszF2YqYtmucmCO/nsHP00pwK1wkniF53CAQC+P1+ioqKXH+v1Oj7rTE61ZkWBCExZF6tqoD86auRHJyeTcUUqG1VX/kjPMstRsWPeFi3DhYsgKVL6zccPKjyNPr0UXX+HNKjh2rcMmCA0h8rMv3995HTxa0i0zU1jU1sXp69yHRwmodppiOZ/URFpmXxoSDoTzLrTNuOTA8aNIg1a9bYPnBtbS2DBg1i/fr1cQkTBEFIWeqrgPh27iR3z1a6HNpK3t6t+LZtVa40+LZtW0Lfeu9e1bFxwYKgVtrr18OJJ6r7OMjIUAbSNL12zGSkpTFWzVis9m3duiFFJVY1j2jtxK247rrI1TyybDQfFDMtCIITbJvpNWvWOM6TWbNmDQejtbgSkkp5ebnXEhyhm17QT7PojZP8fMubkZ9PVed8duXks7dlV8oBI4EVP77/Hq66Cu6/XwW8E0lenvofITOzPGpkuro6enTXNM5vvNGwqDHS+5lY7XP4cGMzbeoKjlBbzYn5863f06xdHQt319CnyBy2yQDW8T4/YADrHL/2jTe8+cckZa4TNhG9+uMozWPs2LFkZ0dbbqPw+XxJ6YIoOGPSpElatAA10U0v6KdZ9MZJhLQNHw01qqv7DeSiL77jL1+Uk9m+J7m50KIJq1QOHoTqA/BCW2g5PuxJM6n48stVDkgc9OsHa9fCgQOTMIxlEQ3lli2x87QNQ0m54w646SbrTpKxslGqq1WXRauc6eA/L5MmTcJue+7USPOwrzcVaEUVc9hIK6ocv/bf/4YnnlDlG086yQVxEUiZ64RNRK/+2DbTV155peOD+3w+OkUqRiokHd0qsuimF/TTLHrdI6tiD9OBDtU7YS/q1gRy6m/sjLLT99E6zESnf3948UXo0GH6kXrTVvz3v9C7d/RjZWWpyDKoTJcePRrvE6vLolVpPKsFiNOnT+e116yfD6euzlldaneY7rUAx0wHro3jdTNmqCynKVPUNynTpln/Y5VodLpOgOhtDti+rDzzzDMuyhCSgW4LJXXTC/ppFr0ukp9PQVhOQZ0Bhw5BVVVDfnJmpjKWGb763GWUkaw+rMyfLwPatrGX69uU3A+zokf//gX1bcWt9/vsM/jBD6Ifq0cPqKxU5nbTptjm2zAad0m0Ko1nlX7iZE7YjUy7m+ah0RyupymKe/RQ5dhXroQLLoDJk2HcOHfHWKvrBKK3OeD5/+iCIAjNEotUkAwaIsym0fzmG2Uk9+xR9xkZcPLJ6ta5Y6NDuIbZlrtbt3oTH8HsfPYZDB8e/Vj5+ereMBqKjETDMFRHxeDIcvACxF/+Ek4/vXGaR7RItNVzNTWh/5REOkfJUkw8w4fD2WfDnDmwcCE88AD07eu1KkFIDGKmBUEQPMDnUyazqd0LE0W7dnDnnfDpp0SNTH/1FRx7bPRjmTnVdXXqH4Zhw5zrCY5MA3z+ufUCxGCclsYT05xcWrWCqVPhiy/gN7+BggK4/Xab37oIQgojTVvSiAULFngtwRG66QX9NIted9FN78yZ8PXXC6JGpq1qRodzwgkN+27aFLnBi0m00njBRstqPydj7KRpi3voNScg8YqPO041HTruOBgzBr78MrHH1+1zJ3r1R8x0E9GpA2JZWZnXEhyhm17QT7PodRfd9ALs3VtGTY11ZLq8XKVjxGLgQHVfW6vyw1u1ir5/pNJ4wZFpw7COSDsZY7sLEN2NWOs1J57jcsqAFYzkG3o6utEz8s13dE8Ki3vyxn960u6EnhzoGH1/J7eyX/wicr/7FES364QuepPZAdFnSA27uDA76pSWlkoyviAIzYarroJevaCkpHEnw+XLVQ70LbfEPo7Pp5qn1NTAU09F3gfgvffgnHNCTeztt8OVV6ptJ54If/mLqq/9+efQsqVayHn22eq15nHOPhv++U/1+NNPG5dju+QSlXJy/fXq5xYtlOEXIrOV7vRgh9cynJOfr+o4CmlPMvya5194CYIgCKlDhw6wc6d1Huu//w2jRtk7zqmnwiefwDXXxN7XKqRz8KCqRW2aXcNoMM12cqYjLUD0Ps1DL7aRTx02Ot1Y0DPf2f6HDsO+fapVfVZTf0+J7mokCFFI2GXl73//OwMGDKCbTGBBEARt6dBBmWAr0/nJJ6oRix0CAVWt4c03I+8zaVLkqLVppsMbtjSF8Hbi8r1sbE7BukGRHQyHgeGWQM12GH0VPPSQ+kZCEHQgYTnTixcv5tNPPwWUsRYEQRD04wc/gP/3/xpHpmtqVB5zy5b2jtOnjzLCHaOU97v8cnUfLTJt5ltbRaadVoEwFzU6QQxdcuneXdWlvu02VS5SEHQgYWa6W7duPPnkk8yaNYs333yT75vQiUtwB7/f77UER+imF/TTLHrdRTe9APPn+9m0qXFk+oMP4IwznB0rVkqFaYajmWnTvBuGynEO3j87O3SMY0Waq6vt1Zl2F/3mRLI1d+2q6lBPnNjQsMcJun3uRK/+JMxM33XXXcyZM4cePXrw4YcfMm7cOAoKCpgwYQIPP/wwmzdvTtRbCXFyi51VQymEbnpBP82i11100wtQXKw054flu776qupgl0iimW3T+AZHws0KI8Fm2ukYt4gv/TeB6DcnvNA8ZAj8/OeqHrUjtm/nlnbtYPt2V3S5gW7XCd30JoOELsXo06cPffr0oXfv3px11lkYhsHnn3/ORx99RElJCccddxyTJ09O5FsKDhgeq21ZiqGbXtBPs+h1F930AowYMZzNm1XbZxPDgP/8p3F1jKYSK3Lt8zWkZQRHpk2ysiKPsVWU2jBSoUGIfnPCK80//zmsXQt//nNDBZaYbN/O8OeeU20zu3d3VV+i0O06oZveZOBKnemzzjoLAJ/PR//+/bnsssswDIOKigo33k4QBEFIIL16hf782mvwk58kPi0iuIZ0OOa24Dxp00yb28KNcaw0D8NwnjMteMuMGfD3v6ubIKQqcZvp+fPn46RE9VVXXcWIESPifTtBEAQhifTuDR9/rJquPPoo3Hpr4t8jN1fd2/lTEtziPN4FiD5ffK8RvCMjA+bPh9//XrUhF4RUJG4zff3117PFQUH0wYMHc/LJJ8f7dkICWLp0qdcSHKGbXtBPs+h1F930QoPm3/4WfvUr8PtV3mqbNol/L7PSR7iZtiqXZ7UAMSvL2RjHE5n+6itn+8dGvznhtebWrdWcuP56sFPbQLcR1u06oZveZNCkNI+dO3fy6quvsnjxYr7++utEaRJcQoeW58Hophf00yx63UU3vdCgOT8fXn4Z/vpXcOtLxQ4d1H1dXej2a65pHBGOZKYDgQDHHKN+ttOQxWlkOrwLZNPRb06kguaePWHWLHsVPrxX6wzdrhO66U0GTVqAeMYZZzBo0CCqq6tZv349I0aM4LHHHuPoo49OlD4hgbzwwgteS3CEbnpBP82i11100wuhmjt1cve9THN81VWx9w1O8zDJzlZ6DxxQz51/fuj+VnifM63fnEgVzaeeqlrCFxfDH/8Yeb/UUGsf3a4TuulNBk2KTM+bN4+PPvqITz75hF27djFq1CjGjBkjUWpBEATBNjt2qPtoudORItOg0gBycmJHpsOreZitygV9KCyEvDxV4UMQUoW4zXTbtm0ZPHjwkZ/z8vK48cYbee2115g5c2ZCxAmCIAjpQ3i6RzCRSuMFYyfNw/vItNBUolb4aNUKBg5saJ0pCEkgbjN98cUX8+yzzzba3qtXL4466qgmiRIEQRDSgyeeaHgcbKbDTa9Vmke4mbaTDy1mWn/MCh+/+51FhY+BA1Vx6oEDPdEmpCdxm+kHHniAv/3tb1x77bVs2rTpyPba2lq2bduWEHFCYpk4caLXEhyhm17QT7PodRfd9ELyNQcb5OA0j5yc0P2sItOtWoXqdZrm4Q36zYlU1Ny6NTz9NNxwA+zdG/qcbp870as/cZvpDh06sGrVKrKysjj++OPp3bs3P/rRj+jXrx+jRo1KpEYhQejWtUg3vaCfZtHrLrrpheRrDjbTwZHp4DbiYG2mc3JC9Qab6dRdgKjfnEhVzT17qvrTl14KwT3hdPvciV798RlOOq9EYNeuXfzv//4vhw8f5qyzzkqLNI+ysjKGDBlCaWkpBQUFXssRBEHQkqefhkmT1OMDB1QjF58PbrlFNYsB9fNDD8G778Lrrze89vHH4aabGn6+4gp47jn1+OOPIWhZD6C6OC5b5k7NbKExTXcX9vjnP2HOHHjhBRWxFoRgkuHXmlQaz6RTp074/f5EHEo7ioqKaN++PYWFhRQWFnotRxAEQSuC0y6CzVdeXuh+hw41rj0dvsYsODJdWdn4vQxD1qU1R84+W91PmCCGWmggEAgQCAT43k6nnyaSEDOdzpSUlEhkWhAEIU46d254HJzmERw9njkTevWCDz4IfW14XnWwmZ4yxfr97FT8EPTDNNQ//zksWiSGWuBIkNOMTLtJk+pMC3qxatUqryU4Qje9oJ9m0esuuumF5Gvu2rXhcXBkOthM9+1rnTKQmxuqNzjKfehQ5PeM1UHPXfSbE7poPvts+M1vYPjwVRw44LUa++h2ndBNbzIQM51GzJkzx2sJjtBNL+inWfS6i256Ifmag5fYRIpMZ2RY16A+5ZRQvbEWIJrbwhcyJhf95oROms85B2AO48fD7t1eq7GHbtcJ3fQmAzHTacSiRYu8luAI3fSCfppFr7vopheSr7lLl4bHTs109+6hemOZaW9NtIl+c0IrzevW8dbuz3lw0jrGjwcdGjLrdp3QTW8yEDOdRuTm5notwRG66QX9NIted9FNLyRfc3Y2jBmjHkdK84hkpiFUbywz7X2NaQD95oRWmquqyN2wgYF9qpg/X1WKefddr0VFR7frhG56k4GYaUEQBMFTzIix08h0OJEqg5jI4sP0ondvePVVeP55mDYNamu9ViQ0V8RMC4IgCJ5imuBgA3zccQ2P7Zrp4DSO1I1MC8mkdWvVerxvX/D7LdqPC0ICEDNdz5/+9CcKCgrIzs5mxowZXstxheLiYq8lOEI3vaCfZtHrLrrpBW80d+qk7k3DfP75obnUGRkNUcVf/CL0tcF6Y6V5pEaNaf3mhG6ardReeSXMnQt33gm/+hXs2pV0WRHR7Tqhm95kIGa6nh49enDvvfdy4YUX4gvvDNBM6NWrl9cSHKGbXtBPs+h1F930gjea+/ZV96aZDjfCeXmwb59q2hJ+eQ7WGyvNI7xFuTfoNyd00xxJbe/e8OKLMHYsXHYZPPAAVFUlVZolul0ndNObDMRM13PBBRcwZswY2rVrRwI6rKckt956q9cSHKGbXtBPs+h1F930gjeazZSOSJfeDh1g505llsPNdLDeWJHp1DDT+s0J3TTHUnvWWbB8ORx9NIwaBc8+610+dV0dXH21ZuOr4XXNbWQ5hiAIguApppmuq7M2wR06wHffKbMc7YtDPdI8BFe5/HJ1P3KkKhUTgQzg58AE4MAvYPf1avfWrSEzwSUUDVSjoKyg+XnoMOzfXz9PDcg7sRu+0tWJfWMhaYiZFgRBEDylTx91bxiq0UaHDqHPd+wI33yj9muKmU6NyLTgKnv2qPudO23t7gPy6m8cAvYnXpIPCF/72rL+ZlLx33oNgpakpZleuHAhN954IwBnn302b7zxhseKksOGDRs4/vjjvZZhG930gn6aRa+76KYXvNHcurVadFhXB59/Dv37hz6fmws7dsAPf9hgptu1a6xXj5zpDYBec0Irzfn5bKir4/gmlG6pM6CyUuVT+1DzpmWr+m9GbB6jtg727lXztV07yPBBdY2KULdoAdmmvLo61u/dS8uaThzerf5xTHV0vK65jRY50xUVFUyZMoXhw4fTpUsXMjIyIlbcqKiooKioiPz8fHJychg8eDAvvPBCyD6XXXYZ+/fvZ//+/ZZGurkuQJwyZYrXEhyhm17QT7PodRfd9IJ3mvv0UWZ648bGZtrnUwsQ27ZVPx9/PHz1lXocrFePNA/95oRWmlevZsppp8GWLXHfMrZuIW/PFjof3ELLnVt4b+EWfjNhC8MHbOHOK7fw/97YgvFN49ftX7+Ffz6/hTm3bWH0oC189f4WOlSo47FlC1k7tpBTvoXsb4Ne9/rr/Laykp1znuaPf/R68Oyh43XNbbSITJeXlzNv3jxOPvlkxo4dy/z58yMa3nHjxrF69Wpmz55N//79WbhwIYWFhdTV1VFYWBjxPWpra6murqampobq6mqqqqrIzs4mI0OL/zdsMXfuXK8lOEI3vaCfZtHrLrrpBe80d+yo8qI3boSLL278/PffKzN98KAy3e3bq+3BevWITOs3J3TTnMg5nJcHo0erm2HARx+pBYtr1kCPHmoe7tqlbnl5UFAAJ5wAv/mNKuloSy9w9Kkw4x6V5pTq0Wkdr2tuo4WZPvbYY9lTnwe1a9cu5s+fb7nf8uXLefvttwkEAkyYMAGAc845h82bN1NcXMyECRMimuP77ruPe++998jPM2fO5JlnnuHKK69M8Nl4h27lbHTTC/ppFr3uopte8E7zj34E778PGzZAv36Nnz90qKGaR7BRDtYbHHlOXTOt35zQTbNbc9jng1NPVTfDUGnZ+/erf/KC66I7pReAT9VQf+wxuOeeRCl2Bx2va26jXdg1Wtm6JUuW0KZNG8aPHx+yfeLEiWzbto0PP/ww4munT59OXV1dyK05GWlBEIRU5n/+p6FEWevWjZ+/6y444wxlaCJ1Q8zJaXicumkeQnPA54OuXVWN9KYY6WCGD4f33lN51YJeaGemo7FmzRoGDBjQKPp80kknAbB27dqEv+eoUaPw+/0htzPOOIOlS5eG7Ldy5Ur8fn+j10+ePJkFCxaEbCsrK8Pv91NeXh6yfdq0acyePTtk29dff43f72fDhg0h2x999NFGXYoqKyvx+/2sWrUqZHsgEGDixImNtE2YMEHOQ85DzkPOIynncffdxdxxB5inE34e06fDKafA+vUBvv3W+jxWr244D2WmVwIN52FGpidPngyEngeU1e9bHrZ9GjA7bNvX9ftuCNv+KI3771XW77sqbHsAaHweqljb0rBtoefRQGqfRyrMK50+H9988zU7dvh5/HG9z8PL30cgEDjixXr37s3JJ59MUVFRo+MkHEMzdu7cafh8PmPGjBmNnuvXr59x3nnnNdq+bds2w+fzGbNmzUqYjtLSUgMwSktLE3ZMt0nk+ScD3fQahn6aRa+76KbXMFJf8x13GMaxxzb8HKz33XcNQ9low+jbt+GxeXvhhYbXtWjR+Pnk3GZ59L7J1+wVqT6HQygtNWaBYdR7iT17DGPkSMOoq/NYVxS0Gl8jOX6tWUWmhehUVlZ6LcERuukF/TSLXnfRTS+kvuZ27VTJMZNgveecA0OHqsexcqa9W1ue2uNrjV6aU30OhxOstn17NY8DAc/kxES38U0GzcpMd+rUiV27djXavnv37iPPpzORygmmKrrpBf00i1530U0vpL7mTp0a+nJAqN6MDFVhAVLZTKf2+Fqjl+ZUn8MhtGrFjIEDQxL6f/1reOopVdkjFdFqfJNEszLTgwYNYv369dSFrU759NNPATjxxBO9kCUIgiAkiJNOgjvuiL1frAWIzbSdgKAbAwfC2rXqvp6sLLj/frXoVtCDZmWmx44dS0VFBS+//HLI9meeeYb8/HxOO+20hL9nUVERfr+fQCp/JyMIgtBMOO00+N3vYu9nZaaDq4Q0oxYCQjPk9NOhogI++8xrJfpiLkZMxgJELepMA6xYsYIDBw6wf/9+QFXmME3z6NGjycnJYeTIkZx77rncdNNN7Nu3j759+xIIBFi5ciULFy50pbNhSUkJBQUFCT+uG5SXl9O5c2evZdhGN72gn2bR6y666QX9NEfSa1U+b9CghsdN6DbdRMoBfcZXoZfm5jKHp02DqVPh+ec9EBUFXca3sLCQwsJCysrKGDJkiKvvpc3/5jfffDOXXHIJ11xzDT6fj5deeolLLrmECRMmsHPnziP7LV68mCuuuIKpU6dy3nnn8dFHH7Fo0aKo3Q/ThUmTJnktwRG66QX9NIted9FNL+inOZLeWDnTV1zhkqCY6DW+Cr00N5c5fNxxasFtaWmSBcVAt/FNCq7VCWnm6FgaTyethqGfXsPQT7PodRfd9BqGfprD9fr9qixb9+7RS7WtW+dVmbnSFCh1lxzNXqH7HA7mm28MY9y4JIqxgY7j67Zf0yYyLTQdXdJRTHTTC/ppFr3uopte0E9zuF4zIh2pS6KJdwsQ9RpfhV6adZ/DwfTsCUcfDf/6VxIFxUC38U0GYqYFQRCEZodVmocg6Mjtt8OsWV6rEKIhZloQBEFodtTWRn9eSuMJutCtGxx/PLz7rtdKhEiImU4jFixY4LUER+imF/TTLHrdRTe9oJ/mcL2mSU7dNA+9xlehl2at5vC6dSzo0QPWrYu62+23w+9/H/ufxGSg1fgmCTHTTUSnOtNlZWVeS3CEbnpBP82i11100wv6aQ7Xa6Z3hJuOAwdCf/bOTOs1vgq9NGs1h6uqKNu+Haqqou7WoQNcdJHqjOg1uoxvMutM+wxDMsviwaxbWFpaKsn4giAIKcK4cbBkiWrQEmygw//SffEF9OuXXG3phrgLG5SVwZAhqv5dDC9RWwvDh8Mbb4R28xSikwy/JpFpQRAEodlgdjmsqfFWhyAkmhYt4Oqr4emnvVYihCNmWhAEQWg25OWpe7sLEAcPdlePICSSwkJ46SWorvZaiRCMmGlBEASh2WCa6ViR6RYt1P2UKe7qEYREkpkJF14Iy5Z5rUQIRsx0GuH3+72W4Ajd9IJ+mkWvu+imF/TTHK7XTPOIhblf8hci6jW+Cr00azeHHe5/xRXw17+6IsUWuo1vMsj0WoCQPG655RavJThCN72gn2bR6y666QX9NIfrNSPTsTD3S76Z1mt8FXpp1moOX365Gt2RIyE729ZLOgBP7YGa7pDZwk1x1txSVaVaM8ZDt26wenViBaUAUs0jTqSahyAIQurxxBNw002Nt4f/pTMMyMiAV15RJceExCPuwgbdu8OOHV6rSB75+bBlS1LfMhl+TSLTgiAIQrOhXTt7+5kR6eC0kIyM2M1eBCGh5Oc3JPA7wADKy6FL58RLcpVu3bxW4ApipgVBEIRmQ6dOzvYPTgvp3Ru+/DKxegQhKnGmPPiAGbfA5MkwYEBiJQnOkQWITUSnDohLly71WoIjdNML+mkWve6im17QT3O43o4dnb2+R4+Gx8lJS9BrfBV6adZ9Dttl/Hh4+eUEi7GBLuObzA6IYqabSElJCcuWLaOwsNBrKTHRwfAHo5te0E+z6HUX3fSCfprD9VpFpiOVvxsxQn3LbnLddXDNNfHpsFNFRBkfvcZXoZdm3eewXX70I/jXvxIsxga6jG9hYSHLli2jpKTE9feSBYhxIgsQBUEQUo99+xrnTa9erTo2R8LMnz58GJ59Fq691vn75uVBRUX0fb79Fo46yvmxdUXchftceik88gh01i13OolIO3FBEARBcECbNo232S1/l5Xlbqm85JfhE5o7I0fC3/7mtQpBzLQgCILQbLAyrBkO/tKJmRZ0YsQIWLHCaxWCmGlBEAShWZMMM20npUHMtJBojjpKlcirrfVaSXojZjqNmDhxotcSHKGbXtBPs+h1F930gn6a7ei122Ic4je8dupTq2PrNb4KvTQ3xzkcjdNOg3//O0FibKDb+CYDMdNpxPDhw72W4Ajd9IJ+mkWvu+imF/TTbEevEzOdGWf3BftmWq/xVeiluTnO4WiMGpXcVA/dxjcZSDWPOJFqHoIgCKnJWWfBqlUNP+/fH9qcJRwzGm0YsGgRxFPpNDtbVQOJxp490KGD82PririL5FBbC2PGSO50JKSahyAIgiA45IEH1L1Z2SMZkWnJmRa8okUL6NoVtm3zWkn6ImZaEARBaFZkZ6v7115T905MbIsW8b2n/TQPQUg8Y8fC4sVeq0hfxEynEauCv/fUAN30gn6aRa+76KYX9NNspffEE+F//gfOOcd5qoH7OdN6ja9CL83NYQ47ZcQIWLkyAWJsoNv4JgMx02nEnDlzvJbgCN30gn6aRa+76KYX9NNspTc7G95/P77juZ/modf4KpxrvvFGF2TYpDnMYafk5EBuriqT5za6jW8ykAWIcaLjAsTKykpyc3O9lmEb3fSCfppFr7vophf005wIvcELEFeuVFE+N6iqglatKgF9xlfhXPOtt6o2116QjnMY4Pnn4dAhcLtynW7jKwsQhYSi0+QH/fSCfppFr7vophf005xovdEi02YudtOOrdf4KvTSnK5zeNQoWL48IYeKim7jmwzETAuCIAhCPdHMdLTvcTt1in1sJ50YBcEp7dtDRYV0Q/QC+WgLgiAIQj3RqnlEW2S4eXPsY0s1D8FtCgrg44+9VpF+iJlOI4qLi72W4Ajd9IJ+mkWvu+imF/TTnGi98Uam7dey1mt8FXppTuc5/LOfwTvvJOxwlug2vslAzHQa0atXL68lOEI3vaCfZtHrLrrpBf00J0LvzJkNj4Mj03/6U+h+dsrfxUav8VU41+xlFD4d57DJmWfCv/6VsMNZotv4JgOp5hEn5urQs846i/bt21NYWEhhPD1oBUEQBE/55BMYPFhFns3HAE88Yb/Em2HENpB29jHp2xe+/NLevqnIbbfBww97rSI9GTUK3nhD0ooCgQCBQIDvv/+e999/39VqHnFW1BRMSkpKtCmNJwiCIDRm0CC4+mr1ODjNI93NiKAnvXvDV1+p+3TGDHKawU83kTQPQRAEIa3JyICnn1aPg9M8fD7n5fDkG3CF/CPiHaeeCv/+t9cq0gsx02nEhg0bvJbgCN30gn6aRa+76KYX9NOcaL3hkWmnZjo8z7oxSgzr8jMAACAASURBVG/nzs6O6y3Ox9jLBNJ0n8Num2ndxjcZiJlOI6ZMmeK1BEfophf00yx63UU3vaCf5kTrDS+N59RMx64lrfT26RN9r9SK7Kb3nHCbROv9wQ9g7dqEHjIE3cY3GYiZTiPmzp3rtQRH6KYX9NMset1FN72gn+ZE6w2PTGdlOXt9LDM9dGiD3rPPjrxfapUGSO854TaJ1puRAV27wo4dCT3sEXQb32QgZjqN0K2cjW56QT/NotdddNML+mlOtN5gM92ihXMzbRVR7toV7r9fPR40qNeR/VLLMEdDSuO5iRt6hw2Dt99O+GEB/cY3GYiZFgRBEIR6gs1zy5ZupHko9DLTgm64aaaFxoiZrufw4cNMnDiRXr160a5dO8444wz+5Xblc0EQBCGlCO5kGI+ZbtnSenu4cY5lplMrZ9o5uuvXnR49VJqH/MOWHMRM11NTU0OfPn344IMP2Lt3LzfddBN+v5+DBw96LS1hzJ4922sJjtBNL+inWfS6i256QT/Nidabk9PwOB4zHfx6k8xMqK1Vj//znwa9+hgd52N8zDEuyLBJus9hk759YdOmxB9Xt/FNBmKm68nNzeWee+6hZ8+eAFx55ZXU1dXxxRdfeKwscVRWVnotwRG66QX9NIted9FNL+inOdF6gyOq8Zjp3NzG2zIzoaZGPa6pqTzyPvqYaedjbDfdxQ3SfQ6bnH46fPhh4o+r2/gmA2knHoENGzZQUFBAeXk5uRZXR7OjjpvtKQVBEITkYxrq996D22+HaBl/ublQWamM8ZYtkJ/f2Ej26QPjx8OsWXDNNfDUU3DmmSpaHcns6N5O/OGHVUtxwTs2boTHH4eSEq+VeEsy/JpEpi2orKzkiiuu4J577rE00oIgCELzJytL1ewFuO46632CF3n17GmdK9yiBdTVqcfBRrs5h7IGDfJagdCvnzLUgvukrZleuHAhbdq0oU2bNowePfrI9urq6v/f3p3HN1Xn+x9/J6xdUlqWCpRV1ipl0UGmFEZUqEDHCswtWBm4lGW8IEi5XnDQEaiMV+HO405BvY4ggmAnMDBYyyCI3NGr/FB0KCiyKMoqKNKWTpsi0OX7+yOmNE3TJmm+Ofm07+fjkUfbk5PkdQ4hfHs4C1JSUtCvXz8sXrzYwEIiIjJSs2bAwoW1z+PJlQzN5pv7TP/hD/avJpP9jAvuSD+A7847jS4gx3nSb9wwuqThEzOYttlsWLRoERITE9GuXTuYzWZkZGS4nTc9PR0xMTEICQnBoEGDsGXLFqd5Jk+ejOLiYhQXF2Pnzp0AgIqKCkyZMgXNmzfHunXrtC9ToOXl5Rmd4BVpvYC8ZvbqJa0XkNeso/ezz+xfQ0NvblGuD7P55vOUltp7TSbgqaeAESPq//z68T2hk87e224D/H31b2nrNxDEDKbz8vKwdu1alJaWYvz48QAAk5tf3SdMmICNGzdi2bJl2L17NwYPHozU1FRYrdZaX+ORRx7BpUuXsHnzZpiNPHpCk+nTpxud4BVpvYC8ZvbqJa0XkNeso9dx4ZawsLoH055s9as6mHb0mkz2wfqgQTU/JriOfed7QiedvQMH3vzl0F+krd+AUALl5eUpk8mkMjIyXO7buXOnMplMavPmzU7TExMTVUxMjCovL6/xOc+cOaNMJpMKDQ1V4eHhlbd9+/bVOP/BgwcVAHXw4MH6L1CASGpVSl6vUvKa2auXtF6l5DXr6P3qK6UApa5dU+rIEfv3M2bYv1a/ffqp/WtV1eeJi1Nq7tybvYBSw4fbf05Pr/l5g+t20OvHFBX5/Y/FY3wP33T8uFL//u/+fU6J61f3eE3k5ldVy1Ebb775JiwWC1JSUpymp6Wl4eLFizjg5tDprl27oqKiAiUlJZW7fxQXFyMhIcGv7UaSdtYRab2AvGb26iWtF5DXrKPXsWW6RQvg9tvtB9M5Tm1XXXQ0MGZM7c9Xdcu0o9fxH6s1bfneuNGHaK34ntBJZ2+vXsDJk/59TmnrNxBEDqZr88UXXyA2NtZlN424uDgAwNGjR/36emPHjkVycrLTLT4+HtnZ2U7z7dmzB8nJyS6Pf/TRR132z87NzUVycrLLfklLly51OVn6uXPnkJycjBPVdop64YUXsLDakTNXr15FcnIy9u3b5zTdarUiLS3NpW3SpElcDi4Hl4PL0eiWo+olxZctWwqbbQVKS+3n7f1pSQAkAziBLl2At992Xo6HHnJaEpw+nYzvvnNejvJy+3I4Dkx0iIiYBIslG1eu3Jz27/++56fXuyk2FgAeBVD9+J7cn+atvl/rUrhefOXmcjh7AUD1Iy+v/jTvvmrTrQBc/zymTeP7KhiWo0kT+8Gv0pfDoa7lsFqtlWOx7t27Y+DAgUhPT3d5Hr/Tts1bo8uXL7vdzaNXr15qzJgxLtMvXryoTCaTev755/3SIHE3DyIiqltxsfOuGyNHKjVpkv376rsz1OSNN5zneeQRpV599eb9f/2rUrm59u/nznWet1+/m/M5pn3xhevr9u9v9K4fwbubBzn7zW+UOn/e6ArjcDcP8itpZyiR1gvIa2avXtJ6AXnNOnrDw4Fz527+7O5qhdHRNT+++rx/+pP9Yi2AvXfChJsHHs6c6Vtjkya+Pc433q9jI0/tx/ewswED/HsQorT1GwgNbjDdpk0b5Ofnu0wvKCiovL+xys3NNTrBK9J6AXnN7NVLWi8gr1lXb+fON793N5h+9dWaH1t1P2jH/tcO1XtDQm5+780VAwM7mOZ7QifdvQMHAocP++/5pK3fQGhwg+n+/fvj+PHjqKh2VMeRI0cAAP369TMiKyi89NJLRid4RVovIK+ZvXpJ6wXkNQei191g2t3W16r//LRq5Xxf9V5fr4IY2LO38j2hk+7euDjg88/993zS1m8gNLjB9Pjx42Gz2bBt2zan6Rs2bEBMTAyGDBni19dLT09HcnJyneewJiIimcrLnQevju/dDWirDpCff97z1+nTBxg8uO757rgDePZZz5/XCA35UunSWCyAzWZ0ReA5DkYMxAGITeueJXjs2rWr8tR1gP3MHI5Bc1JSEkJCQjB69GiMGjUKs2fPRlFREXr06AGr1Yo9e/YgKyvL7YVefJWZmcnTxBARNWDXr9tPk1edJ1um69onuuqgc86cuucBgG7dar4UuVLyL0NOeoSH2wfU4eFGlwROamoqUlNTkZubizs1X99e1GB6zpw5OHv2LAD71Q+3bt2KrVu3wmQy4fTp0+jSpQsAYPv27XjqqaewZMkSFBQUIDY2Fps3b8bEiRONzCciIoFatHDdXQNwv2Xam0uQcwsuBUJcHHDkCBAfb3RJwyRqN4/Tp0+joqICFRUVKC8vd/reMZAGgLCwMGRmZuLixYu4du0aDh06xIE0UOP5JIOZtF5AXjN79ZLWC8hrDkTv9u3AypU3f3YcVFj9HNEOffu6f6769L75pv1r4Afg3jcbuYWc72FX/jwIUdr6DQRRg2mqn7lz5xqd4BVpvYC8ZvbqJa0XkNcciF6LxfmsG8OH27+6GzA67l+71vW++vRW/y/6AQNqnq/aBYD9gO8JnQLRO2CA/wbT0tZvIHAw3YgkJiYaneAVab2AvGb26iWtF5DXHOje4cOBjh2BU6eAUaNqn7dTJ9dp1Xu92crcrp39q2MQ/8ADNc/Xo4fnz+kZvid0CkRvp07At9/657mkrd9A4GCaiIjIQ9262W/du7ueQ1q3/v2dfx49OrCvT3KZTPZ9/L3Zn588J+oARCIiIiO9/rp/n8+bLdMmEzBs2M2fExKAK1eAqCj/NlHD1LUrcPas/RdB8i9uma4nSeeZzs7ONjrBK9J6AXnN7NVLWi8grznQvSZT/Q6uq957yy11P6bqgDsiwvm+wFwJke8JnQLVe/vtwNGj9X8eKes3kOeZ5mC6njIzM5GTk4PU1FSjU+okYcBflbReQF4ze/WS1gvIa5be27Zt3Y8pLa17ngcf9DHII96vYyPP5iH9PaGLvwbTUtZvamoqcnJykJmZqf21TErxLJe+cJwE/ODBg7xoCxEROTGZgF27PNuv2THwrP6vcVQUUFgIfPwx8POf2+//4Qf7/I6DEYuL7Vur58wBXnoJmDDBfiVFb668qIPNBoSFGdtAzgoKgPnzgU2bjC4JrECM17hlmoiIKAj17Gn/WnUrb3T0zYF09fsA+zmxgwE30wWf1q3tA2ryPw6miYiIglBWlv3r4MH2LdKeGjFCSw41AE2berbbEHmHg2kiIqIg5NjqXHW3jupq2gJ8//36mki2Pn2Ar74yuqLh4WC6EUlLSzM6wSvSegF5zezVS1ovIK+5ofR26ODb8zkG09Wvjlhdy5a+Pb9dw1jHwSqQvXFxwJEj9XsOaes3EDiYbkSkXbVIWi8gr5m9eknrBeQ1N5Tems5+4c0FNpYvr/3+Dz7w/LlcNYx1HKwC2Xv77cCxY/V7DmnrNxB4Ng8f8WweRETkjskEvPsuMHKkZ/PGxLhe7vnECSA2tvaD+YqKgFatXOepPjj/5BPgrrs8a/cHns0jOP3zn/Yzvzj2x28MeDYPIiIioby53HhNW6Y92dTFzWHkjVat7ANq8i9eTrye0tPTERkZidTUVBEXbiEiosCoaz/mqsw1bNriYJrId1arFVarFYWFhdpfi1um60nSFRD37dtndIJXpPUC8prZq5e0XkBec7D2FhUBP/uZ6/Saejt3rnnLdLNmdb9OYAbTwbmO3QnW94Q7ge6NiKjf1mkp6zeQV0DkYLoRWblypdEJXpHWC8hrZq9e0noBec3B2mux1Dy9pt4pU2oeTPfqZb/CYW0CM5gOznXsTrC+J9wJdG/PnsA33/j+eGnrNxB4AKKPJB6AePXqVYSGhhqd4TFpvYC8ZvbqJa0XkNfcEHoXLwa2bAFOnfL++QoKgDZt6j4A8fp1oEUL75/f7iqAm83DhwMffuh+7ieeAH7/e+/2GfenhvCe0GnTJqB5c2DSJN8eL2398gBE8itJb35AXi8gr5m9eknrBeQ1N4RepWreMu0JT3YFAeyDJ4eNG93PV/Plpp2b6zoryPPPGzeQBhrGe0Knvn3tZ4rxlbT1GwgcTBMRERmorMzzQXF1Fov3u3r06+f8c2Tkze+9OWiSZOrTB/jyS6MrGhYOpomIiAx044bzluNAePXVm99PmOCf5+RlzGWIiLAfIEv+w8F0I7Jw4UKjE7wirReQ18xevaT1AvKaG0Lv9et6B9O+7EIyfrz9q33QbW8eNsw+7a67gMOHXR+zdKnzVm6jNIT3hG5NmwKlpb49Vtr6DQQOphuRLl26GJ3gFWm9gLxm9uolrReQ19wQem/cqM/BgXV7+umaOoB77nH/mI4d7V9nzAAAe7PjoMOmTYEBA1wfEx8P5OXVK9UvGsJ7QrfevYGvvvLtsdLWbyDwbB4+kng2DyIiCj7TpgGXLgG7dvnvOatujV6yBMjIuDnt0CFg4MCb802fDrz2mv1nxy4njz4KvPSS88GRju//+lf7riHVt3hzNCHHpk32/fQfesjoEv0CMV7jFRCJiIgMtGqV3oFo587OP3t7MoZ27YDLl2/+7OuZRyh4xMUB27YZXdFwcDBNRERkoFat9D13SAgwc6Zn89a06wYA/PAD8Oab/msi49X39HjkjPtMNyInhP3NkdYLyGtmr17SegF5zeyt3S9+Ufv9v/71ze/dbbE+ceJE5QGJEvA9UbeWLe0HvvpC2voNBA6mG5FFixYZneAVab2AvGb26iWtF5DXzN7a7d7tOq3qbhqbNgELF7qee7oqd82HDtUzThO+Jzzny+5F0tZvIHAw3Yi8+OKLRid4RVovIK+ZvXpJ6wXkNbO3Zg8+CMyf7zzt7bdrnrdvX2DwYPfP5a7ZcRBjsOF7wjPV94X3lLT1GwjcZ7oRkXY6G2m9gLxm9uolrReQ18zemmVnu04bM8b9/NW3UFb92ZPmvn09DAsAvic806MH8M03QHS0d4+Ttn4DgVumiYiICCaT/RzSABATU/t81R0/rqeJ9Ln1VvtgmuqPg2kiIqJGpLZT2znua9uWp8Br6Pr1Az77zOiKhoGD6UZkxYoVRid4RVovIK+ZvXpJ6wXkNbNXH8euHpKaAfZ66vbbgSNHvH+ctPUbCNxnup7S09MRGRmJ1NRUpKamGp1Tq6tXrxqd4BVpvYC8ZvbqJa0XkNfMXj2qbpWW0uzAXs+YzUCnTsC5c/bLy3tKyvq1Wq2wWq0oLCzU/lq8nLiPeDlxIiKSxmQCTp4EevZ0np6WBnz1FfD//p99njVrgN/8xvXARJPJfgGXceNu/gzwUuJSZWcDFy8Cc+YYXaJPIMZr3M2DiIioEfH3vtBVL/xCsowcCbz7rtEV8nE3DyIiokZu5UrPr4hXdTBuNgPDhulpIv3Cw+1/nsXFgMVidI1c3DLdiOTl5Rmd4BVpvYC8ZvbqJa0XkNfMXv9o186+/2xNgrXZHfZ6Z8yYmq+U6Y7RvcGIg+lGZPr06UYneEVaLyCvmb16SesF5DWz1zuJiUBUlHePMbrZW+z1zgMP1HyRH3eM7g1GHEw3IsuWLTM6wSvSegF5zezVS1ovIK+Zvd555x2gdWvvHmN0s7fY65327YEbN4AffvBsfqN7gxEH042ItLOOSOsF5DWzVy9pvYC8Zvb63+TJzj9LaK6Kvd6bORNYt86zeYOhN9hwMP2Thx56CO3bt0dERAT69OmDtWvXGp1EREQUcG+84f6+++4DOne++fP27cDYsfqbSK9Ro4C9ez0/CJWccTD9k6VLl+Lbb79FUVER3njjDTz22GM4ffq00VlERERBY+9eoOqGyQcfdB5ck0xmMzBjBvDCC0aXyMTB9E9iY2PRtKn9TIFNmjRBREQELA3sPDHrPP0/nCAhrReQ18xevaT1AvKa2auftGb2+iY11f7LUl0n6wiW3mDCwXQVkydPRkhICBISErBmzRq0bdvW6CS/ys3NNTrBK9J6AXnN7NVLWi8gr5m9+klrZq9vTCZgyRIgI6P2+T75JBf79gHl5YHpkoCXE6+moqICOTk5mD59Og4fPowubi5Yz8uJExFRQ7V2bc2XE6eGb+pUYNEioF8/1/uuXwfuvx8YOhT49FPg2WeBu+4KfKM3eDlxTbKysmCxWGCxWJCUlOR0n9lsxrhx45CQkICcnByDComIiIgC79lngSefrPkXqVdesZ/54z//E9i6FVi8GCgpCXxjsBExmLbZbFi0aBESExPRrl07mM1mZLj5fwibzYb09HTExMQgJCQEgwYNwpYtW5zmmTx5MoqLi1FcXIydO3fW+DxlZWUIDw/3+7IQERERBavOnYF77wVeftl5elER8Le/AQ8/bP85MhJ4/HHguecC3xhsRAym8/LysHbtWpSWlmL8+PEAAJPJVOO8EyZMwMaNG7Fs2TLs3r0bgwcPRmpqKqxWq9vnv3TpErZt24aSkhKUlZXhL3/5Cw4cOIBRo0ZpWR4iIiKiYPXYY8C77wIHDtyctnKlffBsrjJyHDsWOHYMOH8+8I3BRMRgulu3brhy5Qree+89PFfLr0Bvv/029u7di5dffhmzZs3C3XffjTVr1mDUqFFYuHAhKioq3D529erViImJQXR0NF588UXk5OQgJiZGx+IYJjk52egEr0jrBeQ1s1cvab2AvGb26vGznwG//a39eynNDuytP7MZeP11+8GIK1cCK1YAly/bL0dfvXfxYuC//sug0CAhYjBdVW3HS7755puwWCxISUlxmp6WloaLFy/iQNVfsaq45ZZb8MEHH6CwsBAFBQX44IMPMGzYML92B4O5c+caneAVab2AvGb26iWtF5DXzF49Bg26+d/3Upod2OsfERFATo79QMSYGOBPf7Kf8aN67+DBwNmzQEGBQaFBQNxgujZffPEFYmNjYTY7L1ZcXBwA4OjRo35/zbFjxyI5OdnpFh8fj+zsbKf59uzZU+Nvn48++qjLORtzc3ORnJyMvGone1y6dClWrFjhNO3cuXNITk7GiRMnnKa/8MILWLhwodO0YcOGITk5Gfv27XOabrVakZaW5tI2adIkQ5cjMTGxxuW4evVq0C5H3759Pf7zCIblSExMrPf7KpDLkZiYqO3vh47lSExMrHE5AH1/z+u7HImJiUHxeeXpcjjWsdGfV54uh6M3GD6vPF2OxMTEoPi88nQ5HOvY6M8rT5fD0Wv051VNy9G0qX1Xjttuy8WDD9qXw9FbdTlSU+1XwzR6OaxWa+VYrHv37hg4cCDS09NdnsffxJ0aLy8vD9HR0Vi2bBmWLFnidF/v3r3Rs2dPvP32207Tv/vuO8TExOC5557DE0884ZcOnhqPiIiIyH5Gj4cfBt56y+gSVzw1HhEREREFtbAwIDwc+P57o0uM0aAG023atEF+fr7L9IKfduRp06ZNoJOCSvX/Ggl20noBec3s1UtaLyCvmb36SWtmr17uepOTATdnG27wGtRgun///jh+/LjLWTuOHDkCAOhX0+V8GpHaTg8YjKT1AvKa2auXtF5AXjN79ZPWzF693PWOHAns3RvgmCDRoPaZ3r17N8aOHYvNmzdj4sSJldNHjx6No0eP4ty5c27PT+0txz44w4cPR2RkJFJTU5GamuqX5yYiIiKSZswY+4VdmjQxusQ+6LdarSgsLMSHH36odZ/pplqeVYNdu3ahpKQExcXFAOxn5ti2bRsAICkpCSEhIRg9ejRGjRqF2bNno6ioCD169IDVasWePXuQlZXlt4F0VZmZmTwAkYiIiBq9O+4ADh8G7rzT6BJUbuR0bPzUScxges6cOTh79iwA+9UPt27diq1bt8JkMuH06dPo0qULAGD79u146qmnsGTJEhQUFCA2NtZlSzURERER+VdCArB/f3AMpgNJzGD69OnTHs0XFhaGzMxMZGZmai4iIiIiIoef/xx44w1g3jyjSwKrQR2ASLWr6QTowUxaLyCvmb16SesF5DWzVz9pzezVq7be1q2BK1cCGBMkOJhuRKpetUgCab2AvGb26iWtF5DXzF79pDWzV6+6ejt1As6fD1BMkBB3No9gwSsgEhERETlbv95+EZdgOVSNV0AkIiIiIjGGDrUfhNiYiDkAMVilp6fzPNNEREREAHr3Br76yugK5/NM68Yt0/WUmZmJnJwcEQPpffv2GZ3gFWm9gLxm9uolrReQ18xe/aQ1s1evunpNJqBlS6CkJEBBbqSmpiInJycgZ3fjYLoRWblypdEJXpHWC8hrZq9e0noBec3s1U9aM3v18qT3rruATz8NQEyQ4AGIPpJ4AOLVq1cRGhpqdIbHpPUC8prZq5e0XkBeM3v1k9bMXr086X3/feDAAeCJJwLTVBsegEh+JekvKyCvF5DXzF69pPUC8prZq5+0Zvbq5UnvnXcCBw8GICZIcDBNRERERH5jsQA2m9EVgcPBNBERERH5Vfv2wHffGV0RGBxMNyILFy40OsEr0noBec3s1UtaLyCvmb36SWtmr16e9jamgxA5mG5EunTpYnSCV6T1AvKa2auXtF5AXjN79ZPWzF69PO296y7gk080xwQJns3DRxLP5kFEREQUCKWlwLhxwM6dxnYEYrzGKyDWE6+ASEREROSsWTMgLAy4cgWIigr86wfyCojcMu0jbpkmIiIicu/VV4GICGDiROMaeJ5p8qsTJ04YneAVab2AvGb26iWtF5DXzF79pDWzVy9veseOBXbs0BgTJDiYbkQWLVpkdIJXpPUC8prZq5e0XkBeM3v1k9bMXr286e3YESgoAK5d0xgUBLibh48k7uZx7tw5UUcNS+sF5DWzVy9pvYC8ZvbqJ62ZvXp527tqFdC9O5CcrDGqFtzNg/xK0l9WQF4vIK+ZvXpJ6wXkNbNXP2nN7NXL295f/QrYvl1TTJDgYJqIiIiItOjUCbh0CSgvN7pEHw6miYiIiEibiROBy5eNrtCHg+lGZMWKFUYneEVaLyCvmb16SesF5DWzVz9pzezVy5fetDSgfXsNMUGCg+lG5OrVq0YneEVaLyCvmb16SesF5DWzVz9pzezVS1pvIPBsHj6SeDYPIiIiosaEZ/MgIiIiIgpiHEwTEREREfmIg+lGJC8vz+gEr0jrBeQ1s1cvab2AvGb26ietmb16SesNBA6mG5Hp06cbneAVab2AvGb26iWtF5DXzF79pDWzVy9pvYHQZNmyZcuMjpDou+++w5o1a/DII4+gQ4cORud4pE+fPmJaAXm9gLxm9uolrReQ18xe/aQ1s1cvab2BGK/xbB4+4tk8iIiIiIIbz+ZBRERERBTEOJgmIiIiIvIRB9P1lJ6ejuTkZFitVqNT6rRu3TqjE7wirReQ18xevaT1AvKa2auftGb26iWl12q1Ijk5Genp6dpfi4PpesrMzEROTg5SU1ONTqlTbm6u0QlekdYLyGtmr17SegF5zezVT1oze/WS0puamoqcnBxkZmZqfy0egOgjHoBIREREFNx4ACIRERERURDjYJqIiIiIyEccTBMRERER+YiD6UYkOTnZ6ASvSOsF5DWzVy9pvYC8ZvbqJ62ZvXpJ6w0EXk7cRxIvJ96mTRv06NHD6AyPSesF5DWzVy9pvYC8ZvbqJ62ZvXpJ6+XlxA3w0UcfISEhAcuXL8dTTz3ldj6ezYOIiIgouPFsHgFWUVGBBQsWID4+HiaTyegcIiIiIgpyTY0OCCavvPIKEhISUFBQAG6wJyIiIqK6cMv0T/Lz87F69WosXbrU6BRtsrOzjU7wirReQF4ze/WS1gvIa2avftKa2auXtN5A4GD6J4sXL8bjjz+OiIgIAGiQu3msWLHC6ASvSOsF5DWzVy9pvYC8ZvbqJ62ZvXpJ6w2ERjmYzsrKgsVigcViQVJSEg4ePIhDhw5hxowZAAClVIPczaNdu3ZGJ3hFWi8gr5m9eknrBeQ1s1c/ac3s1UtabyCIGEzbbDYsWrQIiYmJnr6b4gAAFHxJREFUaNeuHcxmMzIyMtzOm56ejpiYGISEhGDQoEHYsmWL0zyTJ09GcXExiouLsXPnTuzbtw/Hjh1DdHQ02rVrhy1btuC5557DtGnTArB0RERERCSViMF0Xl4e1q5di9LSUowfPx6A+90wJkyYgI0bN2LZsmXYvXs3Bg8ejNTUVFitVrfPP3PmTJw8eRKfffYZDh8+jOTkZMydOxd//OMftSyPUS5cuGB0glek9QLymtmrl7ReQF4ze/WT1sxevaT1BoKIs3l069YNV65cAWA/UPDVV1+tcb63334be/fuhdVqxaRJkwAAd999N86ePYuFCxdi0qRJMJtdf38ICwtDWFhY5c+hoaGIiIhAVFSUhqUxjrS/ANJ6AXnN7NVLWi8gr5m9+klrZq9e0noDQcRguqra9mV+8803YbFYkJKS4jQ9LS0NDz/8MA4cOID4+Pg6X2P9+vUe9xw/ftzjeY125coV5ObmGp3hMWm9gLxm9uolrReQ18xe/aQ1s1cvab0BGacpYS5fvqxMJpPKyMhwue/nP/+5GjJkiMv0L774QplMJrV27Vq/dVy8eFFFRkYqALzxxhtvvPHGG2+8BektMjJSXbx40W9jwOrEbZmuTX5+Pnr27OkyvXXr1pX3+0uHDh1w7NgxfPfdd357TiIiIiLyrw4dOqBDhw7anr9BDaYDTfcfDhEREREFNxFn8/BUmzZtatz6XFBQUHk/EREREZG/NKjBdP/+/XH8+HFUVFQ4TT9y5AgAoF+/fkZkEREREVED1aAG0+PHj4fNZsO2bducpm/YsAExMTEYMmSIQWVERERE1BCJ2Wd6165dKCkpQXFxMQDg6NGjlYPmpKQkhISEYPTo0Rg1ahRmz56NoqIi9OjRA1arFXv27EFWVpbbC70QEREREflCzJbpOXPmYOLEiZgxYwZMJhO2bt2KiRMnYtKkSbh8+XLlfNu3b8eUKVOwZMkSjBkzBp9++ik2b96M1NRUA+uB1157Db169YLFYsFtt92GU6dOGdpTmxEjRiAkJAQWiwUWiwUjR440OskjH330EcxmM5599lmjU+r00EMPoX379oiIiECfPn2wdu1ao5PcunHjBtLS0tClSxe0atUK8fHx+Oijj4zOqtXLL7+MO+64A82bN0dGRobRObW6fPkykpKSEB4ejj59+mDv3r1GJ9VK0rqV+N6V9NlQnZTPYIn/xkkaQwBAeHh45fq1WCxo0qRJUF9V+ujRo/jFL36ByMhI9OjRA+vWrfPuCbSddI8q5eTkqAEDBqjjx48rpZT65ptv1JUrVwyucm/EiBEqKyvL6AyvlJeXqyFDhqihQ4eqZ5991uicOh07dkyVlpYqpZT65JNPVMuWLdWpU6cMrqpZSUmJeuaZZ9T58+eVUkq9/vrrqm3bturq1asGl7mXnZ2tduzYoVJSUmo8J30wSUlJUTNnzlQ//vijysnJUVFRUSo/P9/oLLckrVuJ711Jnw1VSfoMlvZvnLQxRHUXL15UTZs2VWfOnDE6xa0777xTLV++XCmlVG5urrJYLJXr2xNitkxLtnz5cvzxj39E3759AQC33norIiMjDa6qnarlSpPB6JVXXkFCQgJ69+4toj02NhZNm9r3smrSpAkiIiJgsVgMrqpZaGgonn76aXTq1AkAMHXqVFRUVODrr782uMy9Bx98EL/85S/RqlWroH4/2Gw2vPXWW8jIyEDLli3xwAMPYMCAAXjrrbeMTnNLyroFZL53JX02VCXtM1hCo4PEMURVWVlZGDp0KLp27Wp0ilvHjx+v3INh0KBBiI2NxZdffunx4zmY1qy8vByHDx/G/v370blzZ9x666145plnjM6q04IFCxAdHY2RI0fis88+MzqnVvn5+Vi9ejWWLl1qdIpXJk+ejJCQECQkJGDNmjVo27at0UkeOXHiBH788Uf06NHD6BTxTp48ifDwcHTs2LFyWlxcHI4ePWpgVcMl5b0r7bNB4mewlH/jpI4hqtq0aROmTp1qdEatEhMTsWnTJpSVleHAgQM4f/484uPjPX48B9OaXbp0CWVlZfjoo49w9OhRvPfee8jKysLGjRuNTnNr5cqVOHPmDM6fP4+kpCSMGTMGRUVFRme5tXjxYjz++OOIiIgAADEHmmZlZaGkpARWqxVpaWk4d+6c0Ul1unr1KqZMmYKnn34aoaGhRueIZ7PZKt+3DhEREbDZbAYVNVyS3rvSPhukfQZL+jdO4hiiqs8//xwnT55ESkqK0Sm1WrlyJdavX4+QkBAMGzYMzzzzDKKjoz1+PAfTfpaVlVW5w31SUlLlh/YTTzyBiIgIdO3aFY888gh2795tcKld9V4AGDx4MEJDQ9GiRQssWLAAbdu2xf79+w0utavee/DgQRw6dAgzZswAYP+vu2D777ua1rGD2WzGuHHjkJCQgJycHIMKnbnrLS0tRUpKCvr164fFixcbWOistvUb7MLDw13+Ef/nP/8p4r/1JQnW925tgvGzoSYSPoOrC+Z/46oLCQkBELxjiLps2rQJycnJLhsNgklJSQnuu+8+/OEPf8CNGzfw1VdfITMzE3/72988fo5GP5i22WxYtGgREhMT0a5dO5jNZrdHqNtsNqSnpyMmJgYhISEYNGgQtmzZ4jTP5MmTUVxcjOLiYuzcuRORkZFO/4Xr4Otv7rp7/U137759+3Ds2DFER0ejXbt22LJlC5577jlMmzYtaJtrUlZWhvDw8KDtraiowJQpU9C8eXPvj3I2oLcqf24l83d7r169YLPZcPHixcppR44cwe233x6UvdX5ewukjl5/vncD0VtdfT4bAtGs4zNYZ69u/u6Niory6xgiEM0OFRUVsFqtmDJlit9adfQeO3YMZWVlSElJgclkQvfu3fHAAw/gnXfe8TzKv8dDynP69GkVGRmpRowYoWbNmqVMJpPbI9RHjRqloqKi1Jo1a9T7779fOf+f//znWl/jqaeeUr/85S9VcXGxOn/+vOrbt6/PRxLr7i0sLFR79uxR165dU9evX1erVq1St9xyiyosLAzKXpvNpi5cuKAuXLigvv32WzVx4kT1xBNPqIKCAp96A9H8/fffq61btyqbzaZKS0vVli1bVFRUlPr222+DslcppWbOnKlGjBihrl275lNjoHvLysrUjz/+qKZNm6Z+97vfqR9//FGVl5cHZbvOs3no6NW1bnX1+vO9q7vX358NgWjW8Rmss9ff/8bp7lXKv2OIQDUrpdSePXtUdHS03z4fdPXm5+ersLAw9de//lVVVFSoM2fOqNjYWLVmzRqPmxr9YLqqvLw8t38oO3fuVCaTSW3evNlpemJiooqJian1zXLjxg01a9Ys1apVK9WpU6fK068EY+/ly5fVz372M2WxWFTr1q3Vvffeqw4ePBi0vdVNmzbNr6dl0tH8/fffq+HDh6tWrVqpqKgoNXz4cPXhhx8Gbe+ZM2eUyWRSoaGhKjw8vPK2b9++oOxVSqmlS5cqk8nkdHv99dfr3auj/fLly2rs2LEqNDRU9e7dW7377rt+7fR3byDWrb96db53dfTq/GzQ1Vydvz+D/d2r8984Hb1K6RtD6GxWSqmpU6eq+fPna2v1Z++OHTvUgAEDlMViUR07dlT/8R//oSoqKjzu4GC6isuXL7v9Q5k5c6aKiIhwebNYrVZlMpnU/v37A5VZib36SWtmb+BIa2evXtJ6lZLXzF79pDUHS2+j32faU1988QViY2NhNjuvsri4OAAIulNZsVc/ac3sDRxp7ezVS1ovIK+ZvfpJaw5kLwfTHsrPz0fr1q1dpjum5efnBzqpVuzVT1ozewNHWjt79ZLWC8hrZq9+0poD2cvBNBERERGRjziY9lCbNm1q/C2moKCg8v5gwl79pDWzN3CktbNXL2m9gLxm9uonrTmQvRxMe6h///44fvw4KioqnKYfOXIEANCvXz8jstxir37SmtkbONLa2auXtF5AXjN79ZPWHMheDqY9NH78eNhsNmzbts1p+oYNGxATE4MhQ4YYVFYz9uonrZm9gSOtnb16SesF5DWzVz9pzYHsbeq3ZxJs165dKCkpQXFxMQD7EZ6OlZ+UlISQkBCMHj0ao0aNwuzZs1FUVIQePXrAarViz549yMrK8vuVwNhrXK/EZvYGjrR29rJXejN72Rz0vX47yZ5g3bp1q7z4gNlsdvr+7NmzlfPZbDY1f/581aFDB9WiRQs1cOBAtWXLFvY2sF6JzewNHGnt7GWv9Gb2sjnYe01KKeW/oTkRERERUePBfaaJiIiIiHzEwTQRERERkY84mCYiIiIi8hEH00REREREPuJgmoiIiIjIRxxMExERERH5iINpIiIiIiIfcTBNREREROQjDqaJiIiIiHzEwTQRkR9t2LABZrPZ7e2DDz4wOlGbM2fOOC3r9u3bvXr86tWrYTab8c4777idZ+3atTCbzcjOzgYAjBs3rvL14uLi6tVPROQLXk6ciMiPNmzYgOnTp2PDhg3o27evy/2xsbGwWCwGlOl35swZ3HrrrXj66aeRlJSEXr16ISoqyuPHX7lyBR07dkRycjK2bNlS4zxDhw7FqVOncOHCBTRp0gQnT55EQUEB5syZg9LSUnz++ef+WhwiIo80NTqAiKgh6tevH+644w6jM1BaWgqz2YwmTZoE7DV79OiBu+66y+vHRUVFYdy4ccjOzsaVK1dcBuInTpzAxx9/jMcff7xyeXr16gUAsFgsKCgoqH88EZGXuJsHEZFBzGYz5s2bh02bNiE2NhZhYWEYOHAgdu7c6TLvyZMn8fDDD+OWW25By5Ytcdttt+F//ud/nOZ5//33YTab8cYbb+Dxxx9HTEwMWrZsiW+++QaAfReJ3r17o2XLlrj99tthtVoxbdo0dO/eHQCglEKvXr0wevRol9e32Wxo1aoV5s6d6/PyerIMM2bMwPXr15GVleXy+PXr11fOQ0QULLhlmohIg7KyMpSVlTlNM5lMLluId+7ciX/84x/4/e9/j7CwMKxcuRLjx4/Hl19+WTnIPXbsGIYOHYpu3brhv//7v9G+fXvs3r0bjz32GPLy8rBkyRKn51y8eDGGDh2KNWvWwGw2o127dlizZg3+7d/+Df/yL/+CVatWobCwEBkZGbh+/TpMJlNl37x587BgwQJ8/fXX6NmzZ+Vzbty4EcXFxT4Ppj1dhvvuuw9du3bFa6+95vRa5eXl2LRpE+Lj42vcfYaIyDCKiIj8Zv369cpkMtV4a9asmdO8JpNJdejQQdlstspply5dUk2aNFHPP/985bT7779fdenSRRUXFzs9ft68eSokJEQVFhYqpZR67733lMlkUiNGjHCar7y8XLVv317Fx8c7TT937pxq3ry56t69e+W0oqIiFRERodLT053mve2229R9991X67KfPn1amUwm9frrr7vcV9cyXLlypXJaRkaGMplM6tChQ5XTduzYoUwmk3r11VdrfO27775bxcXF1dpHRKQDd/MgItJg06ZN+Mc//uF0O3DggMt899xzD8LCwip/jo6ORnR0NM6dOwcAuHbtGv73f/8X48ePR8uWLSu3eJeVlWHMmDG4du0aPv74Y6fn/NWvfuX085dffolLly5h4sSJTtM7d+6MhIQEp2kWiwXTpk3Dhg0bcPXqVQDA3//+dxw/ftznrdLeLkNaWhrMZjNee+21ymnr169HeHg4HnroIZ8aiIh04WCaiEiD2NhY3HHHHU63QYMGuczXpk0bl2ktWrTAjz/+CADIz89HeXk5Vq9ejebNmzvdkpKSYDKZkJeX5/T4Dh06OP2cn58PALjllltcXis6Otpl2rx581BUVFS53/KLL76ILl264MEHH/Rw6Z15sgyORsA+yB85ciT+/Oc/o7S0FHl5edixYwdSUlKcfvEgIgoG3GeaiCiIRUVFoUmTJpg6dSoeffTRGufp1q2b08+OfaAdHAP277//3uWxNU3r2bMnxowZg5deegmjR49GTk4Oli9f7vK8OpdhxowZ2LNnD7Kzs3HhwgWUlZVh+vTpPr0+EZFOHEwTEQWx0NBQ3HPPPcjNzUVcXByaNWvm9XP07dsX7du3x1/+8hcsWLCgcvq5c+ewf/9+dOrUyeUx8+fPx/33349//dd/RfPmzTFr1qyALsO4cePQpk0bvPbaa7h48SL69OnjsksKEVEw4GCaiEiDI0eO4MaNGy7Te/bsibZt29b6WFXtWlqrVq3CsGHDMHz4cMyePRtdu3ZFcXExvv76a+zYsQN///vfa30+k8mEjIwMPPLII0hJSUFaWhoKCwuxfPlydOzYEWaz6x5/o0aNQmxsLN5//31MmTKlzua6eLsMzZo1w69//WusWrUKALBixYp6vT4RkS4cTBMR+ZFjV4i0tLQa71u7dm2duytU350iNjYWubm5WL58OX73u9/hhx9+QGRkJHr37o2xY8fW+liHWbNmwWQyYeXKlZgwYQK6d++O3/72t8jOzsb58+drfMzEiRORkZFRr3NL+7IMDjNmzMCqVavQtGlTTJ06td4NREQ68HLiRESNVGFhIXr37o0JEybgT3/6k8v9d955J5o1a+ZythB3HJcTX7duHaZMmYKmTfVvr1FKoby8HPfddx8KCgpw5MgR7a9JRFQVz+ZBRNQIXLp0CfPmzcP27dvxf//3f9i4cSPuuecelJSUYP78+ZXzFRcXY//+/XjyySdx6NAhPPnkk16/1owZM9C8eXNs377dn4tQo/Hjx6N58+b48MMPfT5AkoioPrhlmoioESgsLMTUqVPx6aefoqCgAKGhoYiPj0dGRgYGDx5cOd/777+Pe++9F23btsXcuXNdrq5Ym9LSUqctw7feeisiIyP9uhzVnTp1CoWFhQCAkJAQxMbGan09IqLqOJgmIiIiIvIRd/MgIiIiIvIRB9NERERERD7iYJqIiIiIyEccTBMRERER+YiDaSIiIiIiH3EwTURERETkIw6miYiIiIh8xME0EREREZGPOJgmIiIiIvLR/wcwWUbzha5yjwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 1,
     "metadata": {
      "image/png": {
       "width": 350
      }
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(filename='images/mgxs.png', width=350)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A variety of tools employing different methodologies have been developed over the years to compute multi-group cross sections for certain applications, including NJOY (LANL), MC$^2$-3 (ANL), and Serpent (VTT). The `openmc.mgxs` Python module is designed to leverage OpenMC's tally system to calculate multi-group cross sections with arbitrary energy discretizations for fine-mesh heterogeneous deterministic neutron transport applications.\n",
    "\n",
    "Before proceeding to illustrate how one may use the `openmc.mgxs` module, it is worthwhile to define the general equations used to calculate multi-group cross sections. This is only intended as a brief overview of the methodology used by `openmc.mgxs` - we refer the interested reader to the large body of literature on the subject for a more comprehensive understanding of this complex topic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introductory Notation\n",
    "The continuous real-valued microscopic cross section may be denoted $\\sigma_{n,x}(\\mathbf{r}, E)$ for position vector $\\mathbf{r}$, energy $E$, nuclide $n$ and interaction type $x$. Similarly, the scalar neutron flux may be denoted by $\\Phi(\\mathbf{r},E)$ for position $\\mathbf{r}$ and energy $E$. **Note**: Although nuclear cross sections are dependent on the temperature $T$ of the interacting medium, the temperature variable is neglected here for brevity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spatial and Energy Discretization\n",
    "The energy domain for critical systems such as thermal reactors spans more than 10 orders of magnitude of neutron energies from 10$^{-5}$ - 10$^7$ eV. The multi-group approximation discretization divides this energy range into one or more energy groups. In particular, for $G$ total groups, we denote an energy group index $g$ such that $g \\in \\{1, 2, ..., G\\}$. The energy group indices are defined such that the smaller group the higher the energy, and vice versa. The integration over neutron energies across a discrete energy group is commonly referred to as **energy condensation**.\n",
    "\n",
    "Multi-group cross sections are computed for discretized spatial zones in the geometry of interest. The spatial zones may be defined on a structured and regular fuel assembly or pin cell mesh, an arbitrary unstructured mesh or the constructive solid geometry used by OpenMC. For a geometry with $K$ distinct spatial zones, we designate each spatial zone an index $k$ such that $k \\in \\{1, 2, ..., K\\}$. The volume of each spatial zone is denoted by $V_{k}$. The integration over discrete spatial zones is commonly referred to as **spatial homogenization**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### General Scalar-Flux Weighted MGXS\n",
    "The multi-group cross sections computed by `openmc.mgxs` are defined as a *scalar flux-weighted average* of the microscopic cross sections across each discrete energy group. This formulation is employed in order to preserve the reaction rates within each energy group and spatial zone. In particular, spatial homogenization and energy condensation are used to compute the general multi-group cross section $\\sigma_{n,x,k,g}$ as follows:\n",
    "\n",
    "$$\\sigma_{n,x,k,g} = \\frac{\\int_{E_{g}}^{E_{g-1}}\\mathrm{d}E'\\int_{\\mathbf{r} \\in V_{k}}\\mathrm{d}\\mathbf{r}\\sigma_{n,x}(\\mathbf{r},E')\\Phi(\\mathbf{r},E')}{\\int_{E_{g}}^{E_{g-1}}\\mathrm{d}E'\\int_{\\mathbf{r} \\in V_{k}}\\mathrm{d}\\mathbf{r}\\Phi(\\mathbf{r},E')}$$\n",
    "\n",
    "This scalar flux-weighted average microscopic cross section is computed by `openmc.mgxs` for most multi-group cross sections, including total, absorption, and fission reaction types. These double integrals are stochastically computed with OpenMC's tally system - in particular, [filters](../usersguide/tallies.rst#filters) on the energy range and spatial zone (material, cell or universe) define the bounds of integration for both numerator and denominator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multi-Group Scattering Matrices\n",
    "The general multi-group cross section $\\sigma_{n,x,k,g}$ is a vector of $G$ values for each energy group $g$. The equation presented above only discretizes the energy of the incoming neutron and neglects the outgoing energy of the neutron (if any). Hence, this formulation must be extended to account for the outgoing energy of neutrons in the discretized scattering matrix cross section used by deterministic neutron transport codes. \n",
    "\n",
    "We denote the incoming and outgoing neutron energy groups as $g$ and $g'$ for the microscopic scattering matrix cross section $\\sigma_{n,s}(\\mathbf{r},E)$. As before, spatial homogenization and energy condensation are used to find the multi-group scattering matrix cross section $\\sigma_{n,s,k,g \\to g'}$ as follows:\n",
    "\n",
    "$$\\sigma_{n,s,k,g\\rightarrow g'} = \\frac{\\int_{E_{g'}}^{E_{g'-1}}\\mathrm{d}E''\\int_{E_{g}}^{E_{g-1}}\\mathrm{d}E'\\int_{\\mathbf{r} \\in V_{k}}\\mathrm{d}\\mathbf{r}\\sigma_{n,s}(\\mathbf{r},E'\\rightarrow E'')\\Phi(\\mathbf{r},E')}{\\int_{E_{g}}^{E_{g-1}}\\mathrm{d}E'\\int_{\\mathbf{r} \\in V_{k}}\\mathrm{d}\\mathbf{r}\\Phi(\\mathbf{r},E')}$$\n",
    "\n",
    "This scalar flux-weighted multi-group microscopic scattering matrix is computed using OpenMC tallies with both energy in and energy out filters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multi-Group Fission Spectrum\n",
    "The energy spectrum of neutrons emitted from fission is denoted by $\\chi_{n}(\\mathbf{r},E' \\rightarrow E'')$ for incoming and outgoing energies $E'$ and $E''$, respectively. Unlike the multi-group cross sections $\\sigma_{n,x,k,g}$ considered up to this point, the fission spectrum is a probability distribution and must sum to unity. The outgoing energy is typically much less dependent on the incoming energy for fission than for scattering interactions. As a result, it is common practice to integrate over the incoming neutron energy when computing the multi-group fission spectrum. The fission spectrum may be simplified as $\\chi_{n}(\\mathbf{r},E)$ with outgoing energy $E$.\n",
    "\n",
    "Unlike the multi-group cross sections defined up to this point, the multi-group fission spectrum is weighted by the fission production rate rather than the scalar flux. This formulation is intended to preserve the total fission production rate in the multi-group deterministic calculation. In order to mathematically define the multi-group fission spectrum, we denote the microscopic fission cross section as $\\sigma_{n,f}(\\mathbf{r},E)$ and the average number of neutrons emitted from fission interactions with nuclide $n$ as $\\nu_{n}(\\mathbf{r},E)$. The multi-group fission spectrum $\\chi_{n,k,g}$ is then the probability of fission neutrons emitted into energy group $g$. \n",
    "\n",
    "Similar to before, spatial homogenization and energy condensation are used to find the multi-group fission spectrum $\\chi_{n,k,g}$ as follows:\n",
    "\n",
    "$$\\chi_{n,k,g'} = \\frac{\\int_{E_{g'}}^{E_{g'-1}}\\mathrm{d}E''\\int_{0}^{\\infty}\\mathrm{d}E'\\int_{\\mathbf{r} \\in V_{k}}\\mathrm{d}\\mathbf{r}\\chi_{n}(\\mathbf{r},E'\\rightarrow E'')\\nu_{n}(\\mathbf{r},E')\\sigma_{n,f}(\\mathbf{r},E')\\Phi(\\mathbf{r},E')}{\\int_{0}^{\\infty}\\mathrm{d}E'\\int_{\\mathbf{r} \\in V_{k}}\\mathrm{d}\\mathbf{r}\\nu_{n}(\\mathbf{r},E')\\sigma_{n,f}(\\mathbf{r},E')\\Phi(\\mathbf{r},E')}$$\n",
    "\n",
    "The fission production-weighted multi-group fission spectrum is computed using OpenMC tallies with both energy in and energy out filters.\n",
    "\n",
    "This concludes our brief overview on the methodology to compute multi-group cross sections. The following sections detail more concretely how users may employ the `openmc.mgxs` module to power simulation workflows requiring multi-group cross sections for downstream deterministic calculations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Input Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import openmc\n",
    "import openmc.mgxs as mgxs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a model object to tie geometry, materials, settings, and tallies together\n",
    "model = openmc.Model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We being by creating a material for the homogeneous medium."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate a Material and register the Nuclides\n",
    "inf_medium = openmc.Material(name='moderator')\n",
    "inf_medium.set_density('g/cc', 5.)\n",
    "inf_medium.add_nuclide('H1',  0.028999667)\n",
    "inf_medium.add_nuclide('O16', 0.01450188)\n",
    "inf_medium.add_nuclide('U235', 0.000114142)\n",
    "inf_medium.add_nuclide('U238', 0.006886019)\n",
    "inf_medium.add_nuclide('Zr90', 0.002116053)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our material, we can now create a `Materials` object that can be exported to an actual XML file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate a Materials collection and export to XML\n",
    "model.materials = openmc.Materials([inf_medium])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's move on to the geometry. This problem will be a simple square cell with reflective boundary conditions to simulate an infinite homogeneous medium. The first step is to create the outer bounding surfaces of the problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate boundary Planes\n",
    "min_x = openmc.XPlane(boundary_type='reflective', x0=-0.63)\n",
    "max_x = openmc.XPlane(boundary_type='reflective', x0=0.63)\n",
    "min_y = openmc.YPlane(boundary_type='reflective', y0=-0.63)\n",
    "max_y = openmc.YPlane(boundary_type='reflective', y0=0.63)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the surfaces defined, we can now create a cell that is defined by intersections of half-spaces created by the surfaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate a Cell\n",
    "cell = openmc.Cell(cell_id=1, name='cell')\n",
    "\n",
    "# Register bounding Surfaces with the Cell\n",
    "cell.region = +min_x & -max_x & +min_y & -max_y\n",
    "\n",
    "# Fill the Cell with the Material\n",
    "cell.fill = inf_medium"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OpenMC requires that there is a \"root\" universe. Let us create a root universe and add our square cell to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create root universe\n",
    "root_universe = openmc.Universe(name='root universe', cells=[cell])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now must create a geometry that is assigned a root universe and export it to XML."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Geometry and set root Universe\n",
    "model.geometry = openmc.Geometry(root_universe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we must define simulation parameters. In this case, we will use 10 inactive batches and 40 active batches each with 2500 particles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# OpenMC simulation parameters\n",
    "batches = 50\n",
    "inactive = 10\n",
    "particles = 2500\n",
    "\n",
    "# Instantiate a Settings object\n",
    "settings = openmc.Settings()\n",
    "settings.batches = batches\n",
    "settings.inactive = inactive\n",
    "settings.particles = particles\n",
    "settings.output = {'tallies': True}\n",
    "\n",
    "# Create an initial uniform spatial source distribution over fissionable zones\n",
    "bounds = [-0.63, -0.63, -0.63, 0.63, 0.63, 0.63]\n",
    "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n",
    "settings.source = openmc.IndependentSource(space=uniform_dist)\n",
    "\n",
    "model.settings = settings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to generate multi-group cross sections! First, let's define a 2-group structure using the built-in `EnergyGroups` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate a 2-group EnergyGroups object\n",
    "groups = mgxs.EnergyGroups(group_edges=np.array([0., 0.625, 20.0e6]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use the `EnergyGroups` object, along with our previously created materials and geometry, to instantiate some `MGXS` objects from the `openmc.mgxs` module. In particular, the following are subclasses of the generic and abstract `MGXS` class:\n",
    "\n",
    "* `TotalXS`\n",
    "* `TransportXS`\n",
    "* `AbsorptionXS`\n",
    "* `CaptureXS`\n",
    "* `FissionXS`\n",
    "* `KappaFissionXS`\n",
    "* `ScatterXS`\n",
    "* `ScatterMatrixXS`\n",
    "* `Chi`\n",
    "* `InverseVelocity`\n",
    "* `NuFissionMatrixXS`\n",
    "\n",
    "Of course, we are aware that the fission cross section (`FissionXS`) can sometimes be paired with the fission neutron multiplication to become $\\nu\\sigma_f$. This can be accomodated in to the `FissionXS` class by setting the `nu` parameter to `True` as shown below.\n",
    "\n",
    "Additionally, scattering reactions (like (n,2n)) can also be defined to take in to account the neutron multiplication to become $\\nu\\sigma_s$. This can be accomodated in the the transport (`TransportXS`), scattering (`ScatterXS`), and scattering-matrix (`ScatterMatrixXS`) cross sections types by setting the `nu` parameter to `True` as shown below.\n",
    "\n",
    "These classes provide us with an interface to generate the tally inputs as well as perform post-processing of OpenMC's tally data to compute the respective multi-group cross sections. In this case, let's create the multi-group total, absorption and scattering cross sections with our 2-group structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate a few different sections\n",
    "total = mgxs.TotalXS(domain=cell, energy_groups=groups)\n",
    "absorption = mgxs.AbsorptionXS(domain=cell, energy_groups=groups)\n",
    "scattering = mgxs.ScatterXS(domain=cell, energy_groups=groups)\n",
    "\n",
    "# Note that if we wanted to incorporate neutron multiplication in the\n",
    "# scattering cross section we would write the previous line as:\n",
    "# scattering = mgxs.ScatterXS(domain=cell, energy_groups=groups, nu=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each multi-group cross section object stores its tallies in a Python dictionary called `tallies`. We can inspect the tallies in the dictionary for our `Absorption` object as follows. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OrderedDict([('flux',\n",
       "              Tally\n",
       "              \tID             =\t1\n",
       "              \tName           =\t\n",
       "              \tFilters        =\tCellFilter, EnergyFilter\n",
       "              \tNuclides       =\ttotal\n",
       "              \tScores         =\t['flux']\n",
       "              \tEstimator      =\ttracklength),\n",
       "             ('absorption',\n",
       "              Tally\n",
       "              \tID             =\t2\n",
       "              \tName           =\t\n",
       "              \tFilters        =\tCellFilter, EnergyFilter\n",
       "              \tNuclides       =\ttotal\n",
       "              \tScores         =\t['absorption']\n",
       "              \tEstimator      =\ttracklength)])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "absorption.tallies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Absorption` object includes tracklength tallies for the 'absorption' and 'flux' scores in the 2-group structure in cell 1. Now that each `MGXS` object contains the tallies that it needs, we must add these tallies to a `Tallies` object to generate the \"tallies.xml\" input file for OpenMC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate an empty Tallies object\n",
    "tallies = openmc.Tallies()\n",
    "\n",
    "# Add total tallies to the tallies file\n",
    "tallies += total.tallies.values()\n",
    "\n",
    "# Add absorption tallies to the tallies file\n",
    "tallies += absorption.tallies.values()\n",
    "\n",
    "# Add scattering tallies to the tallies file\n",
    "tallies += scattering.tallies.values()\n",
    "\n",
    "model.tallies = tallies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we a have a complete set of inputs, so we can go ahead and run our simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/pshriwise/.pyenv/versions/3.9.1/lib/python3.9/site-packages/openmc/mixin.py:70: IDWarning: Another Filter instance already exists with id=3.\n",
      "  warn(msg, IDWarning)\n",
      "/home/pshriwise/.pyenv/versions/3.9.1/lib/python3.9/site-packages/openmc/mixin.py:70: IDWarning: Another Filter instance already exists with id=4.\n",
      "  warn(msg, IDWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                %%%%%%%%%%%%%%%\n",
      "                           %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                    %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                     %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ##################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ###################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ####################     %%%%%%%%%%%%%%%%%%%%%%\n",
      "                #####################     %%%%%%%%%%%%%%%%%%%%%\n",
      "                ######################     %%%%%%%%%%%%%%%%%%%%\n",
      "                #######################     %%%%%%%%%%%%%%%%%%\n",
      "                 #######################     %%%%%%%%%%%%%%%%%\n",
      "                 ######################     %%%%%%%%%%%%%%%%%\n",
      "                  ####################     %%%%%%%%%%%%%%%%%\n",
      "                    #################     %%%%%%%%%%%%%%%%%\n",
      "                     ###############     %%%%%%%%%%%%%%%%\n",
      "                       ############     %%%%%%%%%%%%%%%\n",
      "                          ########     %%%%%%%%%%%%%%\n",
      "                                      %%%%%%%%%%%\n",
      "\n",
      "                 | The OpenMC Monte Carlo Code\n",
      "       Copyright | 2011-2022 MIT, UChicago Argonne LLC, and contributors\n",
      "         License | https://docs.openmc.org/en/latest/license.html\n",
      "         Version | 0.13.1\n",
      "        Git SHA1 | 33bc948f4b855c037975f16d16091fe4ecd12de3\n",
      "       Date/Time | 2022-10-04 12:30:35\n",
      "  OpenMP Threads | 2\n",
      "\n",
      " Reading settings XML file...\n",
      " Reading cross sections XML file...\n",
      " Reading materials XML file...\n",
      " Reading geometry XML file...\n",
      " Reading H1 from /home/pshriwise/data/xs/openmc/nndc_hdf5/H1.h5\n",
      " Reading O16 from /home/pshriwise/data/xs/openmc/nndc_hdf5/O16.h5\n",
      " Reading U235 from /home/pshriwise/data/xs/openmc/nndc_hdf5/U235.h5\n",
      " Reading U238 from /home/pshriwise/data/xs/openmc/nndc_hdf5/U238.h5\n",
      " Reading Zr90 from /home/pshriwise/data/xs/openmc/nndc_hdf5/Zr90.h5\n",
      " Minimum neutron data temperature: 294 K\n",
      " Maximum neutron data temperature: 294 K\n",
      " Reading tallies XML file...\n",
      " Preparing distributed cell instances...\n",
      " Reading plot XML file...\n",
      " Writing summary.h5 file...\n",
      " Maximum neutron transport energy: 20000000 eV for H1\n",
      " Initializing source particles...\n",
      "\n",
      " ====================>     K EIGENVALUE SIMULATION     <====================\n",
      "\n",
      "  Bat./Gen.      k            Average k\n",
      "  =========   ========   ====================\n",
      "        1/1    1.18505\n",
      "        2/1    1.17297\n",
      "        3/1    1.16184\n",
      "        4/1    1.14929\n",
      "        5/1    1.09928\n",
      "        6/1    1.18675\n",
      "        7/1    1.19772\n",
      "        8/1    1.17470\n",
      "        9/1    1.17208\n",
      "       10/1    1.09993\n",
      "       11/1    1.14342\n",
      "       12/1    1.10127    1.12234 +/- 0.02107\n",
      "       13/1    1.19914    1.14794 +/- 0.02834\n",
      "       14/1    1.18411    1.15698 +/- 0.02199\n",
      "       15/1    1.14556    1.15470 +/- 0.01718\n",
      "       16/1    1.20337    1.16281 +/- 0.01621\n",
      "       17/1    1.13853    1.15934 +/- 0.01413\n",
      "       18/1    1.18208    1.16218 +/- 0.01256\n",
      "       19/1    1.11842    1.15732 +/- 0.01210\n",
      "       20/1    1.15248    1.15684 +/- 0.01083\n",
      "       21/1    1.14903    1.15613 +/- 0.00982\n",
      "       22/1    1.23456    1.16266 +/- 0.01110\n",
      "       23/1    1.18876    1.16467 +/- 0.01040\n",
      "       24/1    1.13591    1.16262 +/- 0.00985\n",
      "       25/1    1.19559    1.16481 +/- 0.00943\n",
      "       26/1    1.16947    1.16511 +/- 0.00882\n",
      "       27/1    1.13198    1.16316 +/- 0.00851\n",
      "       28/1    1.15329    1.16261 +/- 0.00805\n",
      "       29/1    1.16538    1.16275 +/- 0.00761\n",
      "       30/1    1.18229    1.16373 +/- 0.00729\n",
      "       31/1    1.15060    1.16311 +/- 0.00696\n",
      "       32/1    1.15460    1.16272 +/- 0.00665\n",
      "       33/1    1.13875    1.16168 +/- 0.00644\n",
      "       34/1    1.13479    1.16056 +/- 0.00626\n",
      "       35/1    1.21125    1.16258 +/- 0.00634\n",
      "       36/1    1.15914    1.16245 +/- 0.00609\n",
      "       37/1    1.10457    1.16031 +/- 0.00624\n",
      "       38/1    1.17215    1.16073 +/- 0.00603\n",
      "       39/1    1.18462    1.16155 +/- 0.00588\n",
      "       40/1    1.15361    1.16129 +/- 0.00568\n",
      "       41/1    1.14983    1.16092 +/- 0.00551\n",
      "       42/1    1.14087    1.16029 +/- 0.00537\n",
      "       43/1    1.18725    1.16111 +/- 0.00527\n",
      "       44/1    1.19094    1.16199 +/- 0.00519\n",
      "       45/1    1.17371    1.16232 +/- 0.00505\n",
      "       46/1    1.18552    1.16297 +/- 0.00495\n",
      "       47/1    1.14194    1.16240 +/- 0.00485\n",
      "       48/1    1.12045    1.16130 +/- 0.00484\n",
      "       49/1    1.18476    1.16190 +/- 0.00476\n",
      "       50/1    1.17063    1.16212 +/- 0.00464\n",
      " Creating state point statepoint.50.h5...\n",
      "\n",
      " =======================>     TIMING STATISTICS     <=======================\n",
      "\n",
      " Total time for initialization     = 7.5024e-01 seconds\n",
      "   Reading cross sections          = 7.4182e-01 seconds\n",
      " Total time in simulation          = 1.9857e+01 seconds\n",
      "   Time in transport only          = 1.9835e+01 seconds\n",
      "   Time in inactive batches        = 2.4766e+00 seconds\n",
      "   Time in active batches          = 1.7380e+01 seconds\n",
      "   Time synchronizing fission bank = 1.1568e-02 seconds\n",
      "     Sampling source sites         = 1.0777e-02 seconds\n",
      "     SEND/RECV source sites        = 7.5224e-04 seconds\n",
      "   Time accumulating tallies       = 3.8921e-04 seconds\n",
      "   Time writing statepoints        = 3.0584e-03 seconds\n",
      " Total time for finalization       = 1.9732e-04 seconds\n",
      " Total time elapsed                = 2.0619e+01 seconds\n",
      " Calculation Rate (inactive)       = 10094.3 particles/second\n",
      " Calculation Rate (active)         = 5753.6 particles/second\n",
      "\n",
      " ============================>     RESULTS     <============================\n",
      "\n",
      " k-effective (Collision)     = 1.16239 +/- 0.00461\n",
      " k-effective (Track-length)  = 1.16212 +/- 0.00464\n",
      " k-effective (Absorption)    = 1.15435 +/- 0.00325\n",
      " Combined k-effective        = 1.15666 +/- 0.00304\n",
      " Leakage Fraction            = 0.00000 +/- 0.00000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Run OpenMC\n",
    "statepoint_filename = model.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tally Data Processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our simulation ran successfully and created statepoint and summary output files. We begin our analysis by instantiating a `StatePoint` object. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the last statepoint file\n",
    "sp = openmc.StatePoint(statepoint_filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to the statepoint file, our simulation also created a summary file which encapsulates information about the materials and geometry. By default, a `Summary` object is automatically linked when a `StatePoint` is loaded. This is necessary for the `openmc.mgxs` module to properly process the tally data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The statepoint is now ready to be analyzed by our multi-group cross sections. We simply have to load the tallies from the `StatePoint` into each object as follows and our `MGXS` objects will compute the cross sections for us under-the-hood."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the tallies from the statepoint into each MGXS object\n",
    "total.load_from_statepoint(sp)\n",
    "absorption.load_from_statepoint(sp)\n",
    "scattering.load_from_statepoint(sp)\n",
    "# Close the statepoint file now that we're done getting info\n",
    "sp.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voila! Our multi-group cross sections are now ready to rock 'n roll!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extracting and Storing MGXS Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first inspect our total cross section by printing it to the screen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Multi-Group XS\n",
      "\tReaction Type  =\ttotal\n",
      "\tDomain Type    =\tcell\n",
      "\tDomain ID      =\t1\n",
      "\tCross Sections [cm^-1]:\n",
      "            Group 1 [0.625      - 20000000.0eV]:\t6.81e-01 +/- 2.69e-01%\n",
      "            Group 2 [0.0        - 0.625     eV]:\t1.40e+00 +/- 6.00e-01%\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "total.print_xs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the `openmc.mgxs` module uses [tally arithmetic](../examples/tally-arithmetic.rst) under-the-hood, the cross section is stored as a \"derived\" `Tally` object. This means that it can be queried and manipulated using all of the same methods supported for the `Tally` class in the OpenMC Python API. For example, we can construct a [Pandas](https://pandas.pydata.org/) `DataFrame` of the multi-group cross section data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>cell</th>\n",
       "      <th>group in</th>\n",
       "      <th>nuclide</th>\n",
       "      <th>mean</th>\n",
       "      <th>std. dev.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>total</td>\n",
       "      <td>0.668083</td>\n",
       "      <td>0.001798</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>total</td>\n",
       "      <td>1.292060</td>\n",
       "      <td>0.007737</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   cell  group in nuclide      mean  std. dev.\n",
       "1     1         1   total  0.668083   0.001798\n",
       "0     1         2   total  1.292060   0.007737"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = scattering.get_pandas_dataframe()\n",
    "df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each multi-group cross section object can be easily exported to a variety of file formats, including CSV, Excel, and LaTeX for storage or data processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/pshriwise/.pyenv/versions/3.9.1/lib/python3.9/site-packages/openmc/mgxs/mgxs.py:2007: FutureWarning: As the xlwt package is no longer maintained, the xlwt engine will be removed in a future version of pandas. This is the only engine in pandas that supports writing in the xls format. Install openpyxl and write to an xlsx file instead. You can set the option io.excel.xls.writer to 'xlwt' to silence this warning. While this option is deprecated and will also raise a warning, it can be globally set and the warning suppressed.\n",
      "  df.to_excel(filename + '.xls', index=False)\n"
     ]
    }
   ],
   "source": [
    "absorption.export_xs_data(filename='absorption-xs', format='csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code snippet shows how to export all three `MGXS` to the same HDF5 binary data store."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "total.build_hdf5_store(filename='mgxs', append=True)\n",
    "absorption.build_hdf5_store(filename='mgxs', append=True)\n",
    "scattering.build_hdf5_store(filename='mgxs', append=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing MGXS with Tally Arithmetic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we illustrate how one can leverage OpenMC's [tally arithmetic](../examples/tally-arithmetic.rst) data processing feature with `MGXS` objects. The `openmc.mgxs` module uses tally arithmetic to compute multi-group cross sections with automated uncertainty propagation. Each `MGXS` object includes an `xs_tally` attribute which is a \"derived\" `Tally` based on the tallies needed to compute the cross section type of interest. These derived tallies can be used in subsequent tally arithmetic operations. For example, we can use tally artithmetic to confirm that the `TotalXS` is equal to the sum of the `AbsorptionXS` and `ScatterXS` objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>cell</th>\n",
       "      <th>energy low [eV]</th>\n",
       "      <th>energy high [eV]</th>\n",
       "      <th>nuclide</th>\n",
       "      <th>score</th>\n",
       "      <th>mean</th>\n",
       "      <th>std. dev.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0.000</td>\n",
       "      <td>6.250000e-01</td>\n",
       "      <td>total</td>\n",
       "      <td>(((total / flux) - (absorption / flux)) - (sca...</td>\n",
       "      <td>1.110223e-15</td>\n",
       "      <td>0.011435</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0.625</td>\n",
       "      <td>2.000000e+07</td>\n",
       "      <td>total</td>\n",
       "      <td>(((total / flux) - (absorption / flux)) - (sca...</td>\n",
       "      <td>-2.997602e-15</td>\n",
       "      <td>0.002567</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   cell  energy low [eV]  energy high [eV] nuclide  \\\n",
       "0     1         0.00e+00          6.25e-01   total   \n",
       "1     1         6.25e-01          2.00e+07   total   \n",
       "\n",
       "                                               score      mean  std. dev.  \n",
       "0  (((total / flux) - (absorption / flux)) - (sca...  1.11e-15   1.14e-02  \n",
       "1  (((total / flux) - (absorption / flux)) - (sca... -3.00e-15   2.57e-03  "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Use tally arithmetic to compute the difference between the total, absorption and scattering\n",
    "difference = total.xs_tally - absorption.xs_tally - scattering.xs_tally\n",
    "\n",
    "# The difference is a derived tally which can generate Pandas DataFrames for inspection\n",
    "difference.get_pandas_dataframe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, we can use tally arithmetic to compute the ratio of `AbsorptionXS` and `ScatterXS` to the `TotalXS`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>cell</th>\n",
       "      <th>energy low [eV]</th>\n",
       "      <th>energy high [eV]</th>\n",
       "      <th>nuclide</th>\n",
       "      <th>score</th>\n",
       "      <th>mean</th>\n",
       "      <th>std. dev.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0.000</td>\n",
       "      <td>6.250000e-01</td>\n",
       "      <td>total</td>\n",
       "      <td>((absorption / flux) / (total / flux))</td>\n",
       "      <td>0.076103</td>\n",
       "      <td>0.000658</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0.625</td>\n",
       "      <td>2.000000e+07</td>\n",
       "      <td>total</td>\n",
       "      <td>((absorption / flux) / (total / flux))</td>\n",
       "      <td>0.019441</td>\n",
       "      <td>0.000093</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   cell  energy low [eV]  energy high [eV] nuclide  \\\n",
       "0     1         0.00e+00          6.25e-01   total   \n",
       "1     1         6.25e-01          2.00e+07   total   \n",
       "\n",
       "                                    score     mean  std. dev.  \n",
       "0  ((absorption / flux) / (total / flux)) 7.61e-02   6.58e-04  \n",
       "1  ((absorption / flux) / (total / flux)) 1.94e-02   9.26e-05  "
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Use tally arithmetic to compute the absorption-to-total MGXS ratio\n",
    "absorption_to_total = absorption.xs_tally / total.xs_tally\n",
    "\n",
    "# The absorption-to-total ratio is a derived tally which can generate Pandas DataFrames for inspection\n",
    "absorption_to_total.get_pandas_dataframe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>cell</th>\n",
       "      <th>energy low [eV]</th>\n",
       "      <th>energy high [eV]</th>\n",
       "      <th>nuclide</th>\n",
       "      <th>score</th>\n",
       "      <th>mean</th>\n",
       "      <th>std. dev.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0.000</td>\n",
       "      <td>6.250000e-01</td>\n",
       "      <td>total</td>\n",
       "      <td>((scatter / flux) / (total / flux))</td>\n",
       "      <td>0.923897</td>\n",
       "      <td>0.007833</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0.625</td>\n",
       "      <td>2.000000e+07</td>\n",
       "      <td>total</td>\n",
       "      <td>((scatter / flux) / (total / flux))</td>\n",
       "      <td>0.980559</td>\n",
       "      <td>0.003729</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   cell  energy low [eV]  energy high [eV] nuclide  \\\n",
       "0     1         0.00e+00          6.25e-01   total   \n",
       "1     1         6.25e-01          2.00e+07   total   \n",
       "\n",
       "                                 score     mean  std. dev.  \n",
       "0  ((scatter / flux) / (total / flux)) 9.24e-01   7.83e-03  \n",
       "1  ((scatter / flux) / (total / flux)) 9.81e-01   3.73e-03  "
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Use tally arithmetic to compute the scattering-to-total MGXS ratio\n",
    "scattering_to_total = scattering.xs_tally / total.xs_tally\n",
    "\n",
    "# The scattering-to-total ratio is a derived tally which can generate Pandas DataFrames for inspection\n",
    "scattering_to_total.get_pandas_dataframe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, we sum the derived scatter-to-total and absorption-to-total ratios to confirm that they sum to unity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>cell</th>\n",
       "      <th>energy low [eV]</th>\n",
       "      <th>energy high [eV]</th>\n",
       "      <th>nuclide</th>\n",
       "      <th>score</th>\n",
       "      <th>mean</th>\n",
       "      <th>std. dev.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0.000</td>\n",
       "      <td>6.250000e-01</td>\n",
       "      <td>total</td>\n",
       "      <td>(((absorption / flux) / (total / flux)) + ((sc...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.007861</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0.625</td>\n",
       "      <td>2.000000e+07</td>\n",
       "      <td>total</td>\n",
       "      <td>(((absorption / flux) / (total / flux)) + ((sc...</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.003731</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   cell  energy low [eV]  energy high [eV] nuclide  \\\n",
       "0     1         0.00e+00          6.25e-01   total   \n",
       "1     1         6.25e-01          2.00e+07   total   \n",
       "\n",
       "                                               score     mean  std. dev.  \n",
       "0  (((absorption / flux) / (total / flux)) + ((sc... 1.00e+00   7.86e-03  \n",
       "1  (((absorption / flux) / (total / flux)) + ((sc... 1.00e+00   3.73e-03  "
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Use tally arithmetic to ensure that the absorption- and scattering-to-total MGXS ratios sum to unity\n",
    "sum_ratio = absorption_to_total + scattering_to_total\n",
    "\n",
    "# The sum ratio is a derived tally which can generate Pandas DataFrames for inspection\n",
    "sum_ratio.get_pandas_dataframe()"
   ]
  }
 ],
 "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}