Class_Work/NUCE_2100/HW5.ipynb
2024-10-01 20:40:44 -04:00

742 lines
81 KiB
Plaintext

{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"id": "adafa2ac-0f1b-4766-b221-37d33b5246e9",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import sympy as sm\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"id": "7eb9455d-669d-40aa-9c71-2146373f1ff3",
"metadata": {},
"source": [
"Homework 5 - NUCE 2100\n",
"\n",
"**Dane Sabo**\n",
"\n",
"*October 1st, 2024*"
]
},
{
"cell_type": "markdown",
"id": "a6e08ac9-30bf-4c98-842d-f7ac8ec40641",
"metadata": {},
"source": [
"# Problem 1:\n",
"\n",
"The neutron flux in a bare spherical reactor of radius 50 cm is given by \n",
"\n",
"$$ \\phi(r) = 5x10^{13} \\frac{\\sin(0.0628r)}{r} \\frac{\\text{neutrons}}{\\text{cm}^2 \\times \\text{sec}} $$\n",
"\n",
"where r is measured from the center of the reactor. The diffusion coefficient for the system is 0.80 cm.\n",
"\n",
"## What is the maximum value of flux in the reactor?"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "d918e366-30f5-492e-90e8-b257de968b1b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7ef4d8152d20>]"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAHJCAYAAABpOFaGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbUUlEQVR4nO3deVhUZfsH8O+wr+MCCCg0oJi7uOSCaZKKpGlhZWqZqOSuuZSmLS5ZYlkulamYivlGlhZWvpkSob4mmoq4oYYGKgq4MyyBI3N+f/ibiZEZmH1hvp/r4qpz5pznPHPPoLfPuc/ziARBEEBERERkhxws3QEiIiIiS2EiRERERHaLiRARERHZLSZCREREZLeYCBEREZHdYiJEREREdouJEBEREdktJkJERERkt5gIERERkd1iIkRUhyQmJkIkEiE3N9fSXSEjsJXP88iRI+jRowc8PT0hEomQmZlp6S5h4cKFEIlEKvtsJZ5kXkyEyCoo/oByc3PD1atXq70eGRmJtm3bmrQPBw8exMKFC3H37l2TXkcfivio+5k7d66lu2fVrD121vy904ZMJsPQoUNx+/ZtrFixAlu2bIFEIlF77MOfhZOTE5o0aYLRo0er/b0nMgcnS3eAqKqKigosXboUn332mdmvffDgQSxatAijR49G/fr1zX59bbz33nsIDQ1V2WfqBLGusNbY1fS9e+WVVzB8+HC4urpapnNauHjxIi5duoT169fj1Vdf1eocxWdRXl6OQ4cOITExEQcOHMDp06fh5uZmsr7aQjzJ/JgIkVXp0KED1q9fj3nz5qFx48aW7o5GpaWl8PT0NPt1BwwYgMcee8zs160LbDF2jo6OcHR0tHQ3anT9+nUA0OkfD1U/i1dffRW+vr748MMP8dNPP+HFF180RTcBWFc8y8vLTZr0kfZ4a4ysyltvvYXKykosXbq01mOvXr2KsWPHwt/fH66urmjTpg02btyocszo0aMREhJS7dyH6wcWLlyI2bNnAwBCQ0OVQ/e5ubnKY7OysvDSSy+hQYMG6Nmzp/Lc48ePY8CAARCLxfDy8kLfvn1x6NAhtde7cOGC8l/+9erVw5gxY1BWVqZLiHSi7fv/559/0LJlS7Rs2RL//POPcv/t27cRGBiIHj16oLKyUuN1Ll26hMmTJ6NFixZwd3eHj48Phg4dWq0Wo7i4GDNmzEBISAhcXV3RqFEjREVFISMjo8b3oW37htA2VlX3afN5Xr16FXFxcWjcuDFcXV0RGhqKSZMm4d69ezV+7wD1NS3afN907aM62lxn9OjR6N27NwBg6NChEIlEiIyMrLXth/Xq1QvAg9ElBV0+8wMHDqBLly5wc3NDs2bNsG7dOrXXURdPbT93fb+7ABAVFYUePXrgf//7H3r37g13d3dMnz691vPIPDgiRFYlNDQUo0aNwvr16zF37lyNo0KFhYXo3r07RCIRpk6dCj8/P+zatQtxcXGQSqWYMWOGTtd97rnn8Ndff+Gbb77BihUr4OvrCwDw8/NTHjN06FA0b94cS5YsgSAIAIAzZ86gV69eEIvFmDNnDpydnbFu3TpERkZi37596Natm8p1XnzxRYSGhiI+Ph4ZGRn48ssv0ahRI3z44Yda9bOoqAg3b95U2afoqyHc3d2xefNmPP7443j77bexfPlyAMCUKVNQVFSExMTEGv8lfeTIERw8eBDDhw9HUFAQcnNzsWbNGkRGRiIrKwseHh4AgIkTJ2L79u2YOnUqWrdujVu3buHAgQM4e/YsOnXqZHD7NTFF7Gr7PK9du4auXbvi7t27GD9+PFq2bImrV69i+/btKCsr0+p7V5Wu3zdt+mjIdSZMmIAmTZpgyZIleO2119ClSxf4+/vrHEdFYtKgQQPlPm0/81OnTqF///7w8/PDwoULcf/+fSxYsECvftRE3+8uAJw8eRL169dHTEwMxo0bh5deeglhYWFG7R8ZQCCyAps2bRIACEeOHBEuXrwoODk5Ca+99pry9d69ewtt2rRRbsfFxQmBgYHCzZs3VdoZPny4UK9ePaGsrEwQBEGIjY0VJBJJtestWLBAePjrv2zZMgGAkJOTo/bYESNGVGsnJiZGcHFxES5evKjcd+3aNcHb21t44oknqrUxduxYlfOHDBki+Pj4aIjKvxTxUfej7jjFe9Dl/QuCIMybN09wcHAQ9u/fL2zbtk0AIKxcubLW/iniXVV6eroAQPjqq6+U++rVqydMmTKl1vb0bV8dbWOnS6y0/TxHjRolODg4CEeOHKnWrlwuFwRB8/euat8Vr2n7fdOlj+rocp20tDQBgLBt27Ya26z6fn777Tfhxo0bwpUrV4Tt27cLfn5+gqurq3DlyhXlsdp+5jExMYKbm5tw6dIl5b6srCzB0dGx1t8PQdD+c9f3u1tYWCgAELy8vISzZ8/qfD6ZHm+NkdVp2rQpXnnlFSQkJCA/P7/a64Ig4Pvvv8fgwYMhCAJu3ryp/ImOjkZRUZFWw9W6mjhxosp2ZWUl9uzZg5iYGDRt2lS5PzAwEC+99BIOHDgAqVRaYxu9evXCrVu3qh2nyerVq5GSkqLyY0wLFy5EmzZtEBsbi8mTJ6N379547bXXaj3P3d1d+f8ymQy3bt1CWFgY6tevr/JZ1K9fH4cPH8a1a9d06pe27dfEFLGr6fOUy+XYsWMHBg8erLY26eHbbbXR5/tWWx+NeR1d9OvXD35+fggODsYLL7wAT09P/PTTTwgKClIeo81nXllZid27dyMmJgaPPPKI8vhWrVohOjraoD4+TN/v7smTJwE8uO3fsmVLo/aJjIOJkJb279+PwYMHo3HjxhCJRNixY4dO55eXl2P06NFo164dnJycEBMTU+2YH374AVFRUfDz84NYLEZERAR2795tnDdgY9555x3cv39fba3QjRs3cPfuXSQkJMDPz0/lZ8yYMQD+LeA0poefOLpx4wbKysrQokWLase2atUKcrkcV65cUdlf9Q9r4N9bAXfu3NGqD127dkW/fv1UfozJxcUFGzduRE5ODoqLi7Fp0yat/sL+559/MH/+fAQHB8PV1RW+vr7w8/PD3bt3UVRUpDzuo48+wunTpxEcHIyuXbti4cKF+Pvvv43Wfk1MEbuaPs8bN25AKpUa7ck0fb5vtfXRmNfRhSIp3b59OwYOHIibN29We5JLm8/8xo0b+Oeff9C8efNq11DXf0Po+909deoUAGDYsGFG7Q8ZDxMhLZWWliI8PByrV6/W6/zKykq4u7vjtdde0/gH8P79+xEVFYVffvkFx44dw5NPPonBgwfj+PHjhnTdJjVt2hQjR45UOyokl8sBACNHjqz2L3zFz+OPPw5A87+6ayr81aTqv1D1panORvj/miNj0+f9K5Lv8vJyZGdna3WdadOm4YMPPsCLL76I7777Dnv27EFKSgp8fHyUnxfwoF7l77//xmeffYbGjRtj2bJlaNOmDXbt2mWU9g2hT6zM/Xnqwxr7qEhKn3/+efz0009o27YtXnrpJZSUlCiPMcdnDmj/uev73T158iQCAwNVRtfIurBYWksDBgzAgAEDNL5eUVGBt99+G9988w3u3r2Ltm3b4sMPP1Q+QeHp6Yk1a9YAAP744w+1k6etXLlSZXvJkiX48ccf8fPPP6Njx47Geis245133sF//vOfakWdfn5+8Pb2RmVlZa3/qm/QoIHaWF+6dKnaPl1vVfj5+cHDwwPnz5+v9tq5c+fg4OCA4OBgndo0Nl3eP/DgD+333nsPY8aMQWZmJl599VWcOnUK9erVq/E627dvR2xsLD755BPlvvLycrXXDgwMxOTJkzF58mRcv34dnTp1wgcffFDj75cu7etL11jVRjGye/r06RqP0/Z7Z67vm7m/146OjoiPj8eTTz6Jzz//XDnJpTafuZ+fH9zd3dUm7Or6r44un7s+392TJ08iPDxcq76QZXBEyEimTp2K9PR0bN26FSdPnsTQoUPx1FNPaf0vanXkcjmKi4vRsGFDI/bUdjRr1gwjR47EunXrUFBQoNzv6OiI559/Ht9//73av2Ru3Lih0kZRUZHyPj0A5OfnIzk5udp5inmBtP3L1dHREf3798ePP/6o8jhuYWEhkpKS0LNnT4jFYq3aMhVd3r9MJsPo0aPRuHFjrFq1ComJiSgsLMTMmTNrvY6jo2O1EYbPPvtM5V/VlZWV1W5jNWrUCI0bN0ZFRYXB7RtKl1hpw8HBATExMfj5559x9OjRaq8r3o+23ztzfd8s8b2OjIxE165dsXLlSpSXlyv7Udtn7ujoiOjoaOzYsQOXL19W7j979qzWZQXafO76fncrKyuRlZXFRMjKcUTICC5fvoxNmzbh8uXLyse933jjDfz666/YtGkTlixZole7H3/8MUpKSkw6wZi1e/vtt7FlyxacP38ebdq0Ue5funQp0tLS0K1bN4wbNw6tW7fG7du3kZGRgd9++w23b98GAAwfPhxvvvkmhgwZgtdeew1lZWVYs2YNHn300WpFtp07d1Zec/jw4XB2dsbgwYNr7N/777+PlJQU9OzZE5MnT4aTkxPWrVuHiooKfPTRR0aOhu50ef/vv/8+MjMzkZqaCm9vb7Rv3x7z58/HO++8gxdeeAEDBw7UeJ1BgwZhy5YtqFevHlq3bo309HT89ttv8PHxUR5TXFyMoKAgvPDCCwgPD4eXlxd+++03HDlyROVf/fq2byhdYqWtJUuWYM+ePejduzfGjx+PVq1aIT8/H9u2bcOBAwdQv359jd87dRN2muv7Zonv9ezZszF06FAkJiZi4sSJWn/mixYtwq+//opevXph8uTJuH//Pj777DO0adNGJbnRRJvPXd/vbnZ2NsrLy5kIWTuLPa9mwwAIycnJyu2dO3cKAARPT0+VHycnJ+HFF1+sdn5sbKzw7LPP1niNr7/+WvDw8BBSUlKM3HvrVPXx+YfFxsYKAFQenxeEB4+lTpkyRQgODhacnZ2FgIAAoW/fvkJCQoLKcXv27BHatm0ruLi4CC1atBD+85//aHx8fPHixUKTJk0EBwcH5WO2imNv3Lihtu8ZGRlCdHS04OXlJXh4eAhPPvmkcPDgQZVjNLWh7nFeXeNTW3vavP9jx44JTk5OwrRp01Tau3//vtClSxehcePGwp07dzRe986dO8KYMWMEX19fwcvLS4iOjhbOnTsnSCQSITY2VhAEQaioqBBmz54thIeHC97e3oKnp6cQHh4ufPHFFzW+J23bry0mtcVOELT/rujyeV66dEkYNWqU8hHxpk2bClOmTBEqKiqUx6j73mlqT5vvm659VEfb6+jz+Ly6z6KyslJo1qyZ0KxZM+H+/fs6feb79u0TOnfuLLi4uAhNmzYV1q5dq/Zz0/Tea/vc9f3ufvfddwIA4fTp07XGhixHJAhWVNVnI0QiEZKTk5VPfn377bd4+eWXcebMmWqFiV5eXggICFDZN3r0aNy9e1fjk2dbt27F2LFjsW3bNjz99NOmeAtEREQE3hozio4dO6KyshLXr19XThWvr2+++QZjx47F1q1bmQQRERGZGBMhLZWUlODChQvK7ZycHGRmZqJhw4Z49NFH8fLLL2PUqFH45JNP0LFjR9y4cQOpqalo3769MqHJysrCvXv3cPv2bRQXFyMzMxPAg4VGASApKQmxsbFYtWoVunXrpiwQdnd3r/WpHSIiItIdb41pae/evXjyySer7Y+NjUViYiJkMhnef/99fPXVV7h69Sp8fX3RvXt3LFq0CO3atQMAhISEqH0kU/ERKNbx0XQNIiIiMi4mQkRERGS3OI8QERER2S0mQkRERGS3WCxdA7lcjmvXrsHb21vn5ReIiIjIMgRBQHFxMRo3bgwHh5rHfJgI1eDatWsWXyuKiIiI9HPlyhUEBQXVeAwToRp4e3sDeBBIY6+tI5PJsGfPHvTv3x/Ozs5GbZv+xTibB+NsPoy1eTDO5mGqOEulUgQHByv/Hq8JE6EaKG6HicVikyRCHh4eEIvF/CUzIcbZPBhn82GszYNxNg9Tx1mbshYWSxMREZHdYiJEREREdouJEBEREdktJkJERERkt5gIERERkd1iIkRERER2i4kQERER2S0mQkRERGS3mAgRERGR3WIiRERERHaLiRARERHZLSZCFpInzcOp4lPIk+ap7EvLSVPue3hb0z4iIiLSDxddtYANGRswfud4yAU5FqxegIRBCQCg3OcgcsAr7V/BlpNblNvqjkkYlIDosGhk38pGc5/mCBIHAXiQLD28j4iIiKpjImRmedI8ZTIDAHJBjvE/jwdEUNm3+cRm5Tmajhn38ziIRCImS0RERHpiImRm2beylcmMghxyQKj5PHXHCBAgCA92miJZIiIiquuYCJlZc5/mcBA5qCRDDnBQSWDU0eYYYyZLcZ3iOGpERER1HoulzSxIHISEQQlwFDkCABxFjkgYnFBtX2x4bI3HOMABIohU2naAAxxENX+kcsirJVMCBJXEaMLOCfj44MeQrJSgz1d9IFkpwYaMDQBYrE1ERHULR4QsIK5THPpI+uDrXV/j5QEvI9QnFAAQHRaNC7cvIKxhGILEQXi/z/sq2w8fs/vCbkzYOQGVQiUcRY5YN2gdACj3OcDhwYhQlWEibUaWKoVKvPnbm9WSozvld5T7eUuNiIjqAiZCFhIkDkI773YqyUOQOKjG7Yf3xXWKq5Y8AcZJlh5OlNQlR5puqREREdkKJkI2zhTJUnzfeMxNnVutjkndLbWq9UcTdk5AdFg0AHCUiIiIbAITITugT7LU0L1hrcnRwyqFSqw6vArL05fz9hkREdkEJkIEoHqyVFtypOmWmiIJAnj7jIiIrB+fGiONgsRBiAyJVBk5yp2Ri7TYNFyaeQnrB69XebJtVsQsrZ5Iy5Pm8ekzIiKyChwRIp3UdEsNAJYfWs7bZ0REZDOYCJFBHr6lljAogbfPiIjIZtjMrbE1a9agffv2EIvFEIvFiIiIwK5du2o8Z9u2bWjZsiXc3NzQrl07/PLLL2bqrf3i7TMiIrIlNjMiFBQUhKVLl6J58+YQBAGbN2/Gs88+i+PHj6NNmzbVjj948CBGjBiB+Ph4DBo0CElJSYiJiUFGRgbatm1rgXdgP0x1+4yjREREZGw2MyI0ePBgDBw4EM2bN8ejjz6KDz74AF5eXjh06JDa41etWoWnnnoKs2fPRqtWrbB48WJ06tQJn3/+uZl7TlWLrh9eYkTTUiEP3z6bsHMCjlw9whEiIiIyKpsZEaqqsrIS27ZtQ2lpKSIiItQek56ejlmzZqnsi46Oxo4dOzS2W1FRgYqKCuW2VCoFAMhkMshkMsM7XoWiPWO3awtGtRuFPpI+uHjnIpo1aIaUv1Mweddk5ZxF07tOx/LDy1XOqRQq0f3L7pDjwQjRmgFrMKbDmFqvZc9xNifG2XwYa/NgnM3DVHHWpT2bSoROnTqFiIgIlJeXw8vLC8nJyWjdurXaYwsKCuDv76+yz9/fHwUFBRrbj4+Px6JFi6rt37NnDzw8PAzrvAYpKSkmaddWnMRJ+MMf61qtQ35FPgJdA4FiQASRSpE18GDBWODBCNGkXybBMffBqJLiPF8XX43Xsfc4mwvjbD6MtXkwzuZh7DiXlZVpfaxNJUItWrRAZmYmioqKsH37dsTGxmLfvn0akyFdzZs3T2UUSSqVIjg4GP3794dYLDbKNRRkMhlSUlIQFRUFZ2dno7ZdF1SGVCpHiRzgoEyCFOSQI8s7Cyv/XKmsI1I3SsQ4mwfjbD6MtXkwzuZhqjgr7uhow6YSIRcXF4SFPSi47dy5M44cOYJVq1Zh3bp11Y4NCAhAYWGhyr7CwkIEBARobN/V1RWurq7V9js7O5vsF8GUbduy8V3GY2CLgbhw+wI8nT3RfUP3amufKZIg4MEo0eRdkzGwxUAA/6515u/+YFSQcTYPxtl8GGvzYJzNw9hx1qUtmymWVkcul6vU9FQVERGB1NRUlX0pKSkaa4rI+iiKrLs06aJSYK3pMXzF02aSlRL0+aoPJCsl2JS5yRJdJyIiG2EzI0Lz5s3DgAED8Mgjj6C4uBhJSUnYu3cvdu/eDQAYNWoUmjRpgvj4eADA9OnT0bt3b3zyySd4+umnsXXrVhw9ehQJCQmWfBukJ20ew1f3tNnkXZMRHxYPj1wPtPJvxVmriYhIhc2MCF2/fh2jRo1CixYt0LdvXxw5cgS7d+9GVFQUAODy5cvIz89XHt+jRw8kJSUhISEB4eHh2L59O3bs2ME5hGxYTY/h1zRK9Gb2m+if1B+SlRJsyNhgia4TEZGVspkRoQ0bav4LbO/evdX2DR06FEOHDjVRj8jStJ2sUfH0mWI+ouiwaADg2mZERGQ7I0JE6tQ2WePD1NURcZSIiMh+2cyIEJE2qo4SaXraTN2s1e3926PkXglHiIiI7AxHhKjOqfq02ZoBa5QjQzXVEXXf0J0jREREdogjQlSnjekwBo65jpB0lKBlo5YA1NcRPTxCFB0WzZEhIiI7wBEhqvN8XXzRW9JbpzqiC7cvIE+ax0VeiYjqOI4Ikd2prY7IUeSIo9eOou9XfZXLdyQMSkBcpzgL9pqIiEyBI0Jkl2qatTq+bzze/O3NarfLjlw9whEiIqI6hiNCZPceno8o+1a2xoJqjhAREdUtTISIAGX9kIKDyIEF1UREdoC3xogewoJqIiL7wREhIjVYUE1EZB84IkSkAQuqiYjqPo4IEWmBBdVERHUTEyEiLbGgmoio7uGtMSI9sKCaiKhu4IgQkZ5YUE1EZPs4IkRkAH0KqjkyRERkPTgiRGQk2hZUX7h9AQCQfSsbzX2as4aIiMiCmAgRGVFtBdW8XUZEZF14a4zIRB4uqObtMiIi68MRISIT0vZ2WfqVdPh6+PJWGRGRmTERIjKx2m6XiSDC8O+H81YZEZEF8NYYkRlpmn+It8qIiCyDI0JEZlb1dtn10usYtn2Yyut8soyIyHyYCBFZgOJ2WZ40j0+WERFZEG+NEVkQnywjIrIsjggRWRifLCMishwmQkRWgE+WERFZBm+NEVkZPllGRGQ+HBEiskJ8soyIyDyYCBFZKT5ZRkRkerw1RmTl+GQZEZHpcESIyAZo+2TZhdsXeIuMiEgHTISIbERtT5Y5ihzh6eyJtJw01gwREWmJt8aIbJC622Uj249E9w3d0eerPpCslGBDxgYL95KIyPpxRIjIRlW9Xebp7InuG7pXqxmKDovmyBARUQ2YCBHZMMXtsrScNI01QwAfsSci0oSJEFEd0NynOR+xJyLSA2uEiOoAPmJPRKQfjggR1RF8xJ6ISHdMhIjqED5iT0SkG94aI6qj+Ig9EVHtbCYRio+PR5cuXeDt7Y1GjRohJiYG58+fr/GcxMREiEQilR83Nzcz9ZjI8uI6xSF3Ri7SYtOQHpeOLSe3sGaIiKgKm0mE9u3bhylTpuDQoUNISUmBTCZD//79UVpaWuN5YrEY+fn5yp9Lly6ZqcdE1iFIHITIkEiU3CvRWDOUJ81DWk4akyIisjs2UyP066+/qmwnJiaiUaNGOHbsGJ544gmN54lEIgQEBJi6e0RWj4/YExFVZzOJ0MOKiooAAA0bNqzxuJKSEkgkEsjlcnTq1AlLlixBmzZt1B5bUVGBiooK5bZUKgUAyGQyyGQyI/Ucyjar/pdMg3H+l7+7P9YMWIPJuyajUqiEo8gRH0R+oPYR+z6SPjoVUjPO5sNYmwfjbB6mirMu7YkEQRCMenUzkMvleOaZZ3D37l0cOHBA43Hp6enIzs5G+/btUVRUhI8//hj79+/HmTNnEBRU/Q/5hQsXYtGiRdX2JyUlwcPDw6jvgchSbt67ifyKfAS6BiK/Ih/vXny32jGLmy1GO+92FugdEZHhysrK8NJLL6GoqAhisbjGY20yEZo0aRJ27dqFAwcOqE1oNJHJZGjVqhVGjBiBxYsXV3td3YhQcHAwbt68WWsgdSWTyZCSkoKoqCg4OzsbtW36F+NcszxpHsJWh1W7XZY9JRsAlHMS1TY6xDibD2NtHoyzeZgqzlKpFL6+vlolQjZ3a2zq1KnYuXMn9u/fr1MSBADOzs7o2LEjLly4oPZ1V1dXuLq6qj3PVL8Ipmyb/sU4qxfqE4qEQQmYsHOC8nbZukHr8Pul3zF+53id64YYZ/NhrM2DcTYPY8dZl7ZsJhESBAHTpk1DcnIy9u7di9DQUJ3bqKysxKlTpzBw4EAT9JDINj08IzUASFZKuJI9EdkFm0mEpkyZgqSkJPz444/w9vZGQUEBAKBevXpwd3cHAIwaNQpNmjRBfHw8AOC9995D9+7dERYWhrt372LZsmW4dOkSXn31VYu9DyJrVHVGak0r2adfSYevhy9npCaiOsVmEqE1a9YAACIjI1X2b9q0CaNHjwYAXL58GQ4O/06NdOfOHYwbNw4FBQVo0KABOnfujIMHD6J169bm6jaRzVH3mL0IIgz/fjgfsSeiOsdmEiFtarr37t2rsr1ixQqsWLHCRD0iqpsUS3Mo6oYc4AABAm+VEVGdZDOJEBGZT9W6oeul1zFs+zCV17mKPRHVFUyEiEgtRd1QnjRP7YzUYQ3DkCfNQ/atbISIQyzXUSIiA9jMWmNEZBnqVrFfN2gddl/YDclKCfp81Qdhq8OQcivFwj0lItIdR4SIqFbaPGK/5soavC59HaE+uk9tQURkKUyEiEgrtT1iL4ccF+9cZCJERDaFt8aISGeKR+yrcoADPJw9kJaThjxpnoV6RkSkGyZCRKQzdXVDvRv0Rq/NvdDnqz6QrJRgQ8YGC/eSiKh2Rr815ujoiMrKSmM3S0RWpmrdkIvIBT0Te0LAg/m+ONcQEdkKo48I2eBi9kSkpyBxECJDIlF6r1SZBCko5hoiIrJmRk+ERCKRyvacOXNw9+5d5fadO3cwd+5cY1+WiCworGEYRFD93a861xDrhojIWpm8RiglJQX169dXbjdo0AB79uwx9WWJyIyCxEGYHDy5xrmGWDdERNZIp0QoJSUFJSUlAIAvvvgC48ePx7lz52o8Ry6Xo7i4WLktlUohk8n06CoRWbMonyhkT8lGWmwacmfkIjosGuN3jq+2RhlHhojImuiUCL3xxhvw8vLCoUOH8PXXX6Nfv36Ii6t5Berp06ejZ8+eWLJkCZYsWYJevXph5syZBnWaiKyTomYoSByE7FvZ1eYaYt0QEVkbvW6N7dixAxMnTsSLL76IsrKyGo8dO3Ysvv76a3h7e8Pb2xtJSUkYO3asXp0lItuhbq4hR5EjPJ09WTNERFZDp0SocePGeOWVV/DNN99g0KBBqKio0OpR+StXrgAApk2bBh8fH5w9e1a/3hKRzVA319DI9iPRfUN31gwRkdXQKRHavn07hgwZgt9++w0NGjTA7du38fHHH9d4zhtvvIGtW7di9erVAB7MMzR69Gi9O0xEtiOuUxxyZ+QiLTYN6XHp2HJyC2uGiMiq6DShoqenJ5577jnldmBgIAIDA2s8JzU1FcePH0fHjh0BAH5+figvL9ejq0RkixRrlKlbn0xRM8RJF4nIUnSeWfqPP/5Abm6uyi2xUaNGaTze2dkZcrlcOb/Q7du34eDAlT2I7I2iZqhqMlR1rqHsW9lo7tOcSRERmZVOidCIESNQUFCAjh07wtHxwX3/hydQfNhrr72GYcOG4ebNm1i8eDG+/fZbvP322/r3mIhskqJmaMLOCagUKlXmGlI8Zu8gckDCoATEdar5aVQiImPRKRE6ceIEsrKydLrAyJEj8dhjjyE1NRVyuRzbtm1Dq1atdGqDiOqGquuThTUMAwBIVkqq1Q1xjTIiMhed7lF17doV58+f1+kCKSkpCAoKwpQpU+Dg4IAVK1bUOgkjEdVdnGuIiKyJTolQZmYmwsPDER4ejq5du6JLly7o2rVrjedUnYQxKSlJq0kYicg+cK4hIrI0nW6N/fjjj3pfqOokjPHx8Xq3Q0R1h7q6IcVcQ6wZIiJz0CkRkkgkKCgowJEjRwA8uFXm7+9f4zmKSRj379+PzMxMrSdhJCL7ULVuyNPZU5kEAawZIiLT0+nWWFJSEnr27In//ve/2LlzJ3r16oWtW7fWeI4+kzASkX1R1A2V3CthzRARmZVOI0Iffvghjhw5ggYNGgAA7ty5g8jISAwfPlzjOfpMwkhE9olzDRGRuek0IiSXy+Hl5aXc9vLyglwur+EMIiLtqVufTDHXkGSlhGuUEZHR6TQiNHLkSPTo0QPPP/88AOCHH36ocVZpIiJdca4hIjInnUaEZsyYgTVr1sDd3R3u7u5Ys2YNZs+erfH4nJwczJkzR2VfUlIStm3bpl9vicgucK4hIjIXrUeEBEFAx44dkZWVhccee0yrc0JDQwEAL7zwApKSkrBmzRokJycb9Bg+EdkXTXVDirmGWDNERIbQekRIJBIhPDwcZ86c0ekCH330Ebp3747WrVvj999/x6+//op69erp3FEisk/q6oYUcw2xZoiIDKVTjdCZM2fQsWNHPProo/Dw8IAgCBCJRPjzzz81niMIAv766y80a9YMhYWFKC4uhpubm8EdJyL7wbmGiMhUdEqEhg4dqlIcLQgCtmzZovH48vJyDB8+HO3bt0dCQgJ27dqFfv364YcffkCzZs307zUR2Z0gcRCCxEFIy0nTWDPERIiIdKVTsXRycjIkEonyJyQkBNu3b9d4/N27d/H888/jvffeAwAMGDAA69evR14e1w8iIv1oWp9MMdcQ1ygjIl1oNSK0fv16JCQk4Pz58yqLrBYXF6Njx44azwsICMArr7yisq+2RVqJiGqibn0yxVxD43eO5xplRKQTrRKhF198EVFRUXjnnXfwwQcfKPd7e3ujYcOGJuscEZE6nGuIiIxFq0SoXr16qFevHqKiorBv375qr3NSRSIyN0XNEADWDRGR3nQqlj59+rTy/ysqKpCSkoL27dszESIii+IaZUSkL50SoWXLlqlsl5SUICYmRmWfIAjIysrCL7/8gpdeegmNGzdGcnIyhgwZYnBniYjUYd0QEelLp6fGHiYSiXDp0qVq+xctWoSoqCgsWrQIJ06cQEpKiiGXISKqVVynOOTOyEVabBpyZ+QiOixamQQB/9YN8YkyIqpKpxGhLl26QCQSAQAqKyuRn59fbS0xAGjQoAHCw8Oxdu1azJw5EydPnjROb4mIasC6ISLSlU6JUNU5g5ycnNCoUSM4OztXOy4qKgrAgxGjFStW4OOPPzawm0REuqmpboiISEGnW2MSiQRZWVn46aef0KRJE9y+fRtnz55VOUYkEuH5559Xbv/999+4evUqhgwZgmeeeUb5Q0RkSurWKFs3aB0AcNJFIlLSKRF64403sHXrVqxevfrByQ4OGD16dI3nDBkyBG3atMGcOXPw9ttvK390FR8fjy5dusDb2xuNGjVCTEwMzp8/X+t527ZtQ8uWLeHm5oZ27drhl19+0fnaRGSbHq4bAh7MN8TFWolIQadEKDU1FZs3b4a7uzsAwM/PD+Xl5TWe4+HhgXHjxiEiIgLdunVT/uhq3759mDJlCg4dOoSUlBTIZDL0798fpaWlGs85ePAgRowYgbi4OBw/fhwxMTGIiYlRmQaAiOq2IHEQIkMiAYDF00RUjU41Qs7OzpDL5cqC6du3b8PBoeZcat68eXjzzTfRr18/uLq6Kvc/8cQTOnX0119/VdlOTExEo0aNcOzYMY1trVq1Ck899RRmz54NAFi8eDFSUlLw+eefY+3atTpdn4hsW/atbBZPE1E1OiVCr732GoYNG4abN29i8eLF+Pbbb2u9zbV7927s3bsXFy5cUCZNIpFI50ToYUVFRQBQ4xIf6enpmDVrlsq+6Oho7NixQ+3xFRUVqKioUG5LpVIAgEwmg0wmM6i/D1O0Z+x2SRXjbB62EOcQcYja4mmJtwQ5t3KUy3VYe1JkC7GuCxhn8zBVnHVpTyQIgqBL4+fOnUNqairkcjn69u2L1q1bq7zu6OiIyspK5fajjz6K8+fPK0eRjEEul+OZZ57B3bt3ceDAAY3Hubi4YPPmzRgxYoRy3xdffIFFixahsLCw2vELFy7EokWLqu1PSkqCh4eHcTpPRBaTcisFa66sgRxyOMABk4InAQC+uPIFBAgQQYTJwZMR5RNl4Z4SkSHKysrw0ksvoaioCGKxuMZjdRoRAoCWLVuiZcuWWh/ftWtXXLx4EWFhxntkdcqUKTh9+nSNSZA+5s2bpzKCJJVKERwcjP79+9caSF3JZDKkpKQgKipK7RQEZByMs3nYSpwHYiBel76Oi3cuolmDZgCAsNVhEPDg34MCBKzNW4vXn33dakeGbCXWto5xNg9TxVlxR0cbOidCf/zxB3Jzc1VGfWpaa+z06dNo27YtWrRoAVdXVwiCAJFIhD///FPXSwMApk6dip07d2L//v0ICqr5D6qAgIBqIz+FhYUICAhQe7yrq6tKHZOCs7OzyX4RTNk2/YtxNg9biHOoTyhCfUIBaJ508VLxJeUx1soWYl0XMM7mYew469KWTonQiBEjUFBQgI4dO8LR8cHcHLXd8vrxxx91uYRGgiBg2rRpSE5Oxt69exEaWvsfUhEREUhNTcWMGTOU+1JSUhAREWGUPhGRbdM06aKnsyfSctK4UCuRHdApETpx4gSysrJ0uoBEItHpeE2mTJmCpKQk/Pjjj/D29kZBQQEAoF69esrH+UeNGoUmTZogPj4eADB9+nT07t0bn3zyCZ5++mls3boVR48eRUJCglH6RES2Td1irSPbj0T3Dd25UCuRndBpHqGuXbtqNYlhVXPmzMHdu3eV23fu3MHcuXN1agMA1qxZg6KiIkRGRiIwMFD58+233yqPuXz5MvLz85XbPXr0QFJSEhISEhAeHo7t27djx44daNu2rc7XJ6K6qeqki+lx6dhycgvnGiKyIzqNCGVmZiI8PFynep+UlBR89NFHyu0GDRpgz549WLp0qU4d1ebhtr1791bbN3ToUAwdOlSnaxGRfVEs1sqFWonsj06JkD71PnK5HMXFxfD29gbwoJKb8zIQkTXiQq1E9kenREifep/p06ejZ8+eGDZsGADg22+/xcyZM3Vuh4jI1NTVDK0btA5B4iDkSfOQfSubBdREdYzOj8/X5uFbWGPHjkXXrl2RlpYG4MHkhG3atDH2ZYmIjCKuUxyiw6JVZprekLFBuU4ZC6iJ6hajJ0JyubzavrZt27JAmYhshqJmCADypHlqF2uNDovmyBBRHaDTU2NERPampsVaicj2MREiIqqBooC6qqqTLvLReiLbxkSIiKgGigJqR9GD2fSrTrrY56s+kKyUYEPGBgv3koj0ZZQaIZlMhoKCApSVlcHPzw8NGzY0RrNERFahagG1p7OncuZpgDVDRLZO7xGh4uJirFmzBr1794ZYLEZISAhatWoFPz8/SCQSjBs3DkeOHDFmX4mILCZIHITIkEiU3CthzRBRHaJXIrR8+XKEhIRg06ZN6NevH3bs2IHMzEz89ddfSE9Px4IFC3D//n30798fTz31FLKzs43dbyIii9BUM8RJF4lsk163xo4cOYL9+/drnA+oa9euGDt2LNauXYtNmzbhf//7H5o3b25QR4mIrAEnXSSqW/RKhL755hsAQGVlJX7++Wf07dtXuYRGVa6urpg4caJhPSQisjKcdJGo7jDoqTFHR0eMGDECN27cMFZ/iIhsgqJmSDESpG7SRT5aT2T9DH58vkuXLsjJyTFGX4iIbBInXSSyXQYnQtOmTcNbb72FK1euGKM/REQ2h5MuEtkug+cRUqwq36ZNGzzzzDOIjIxEx44d0a5dO7i4uBjcQSIia6eugFox6SJrhoism8GJUE5ODk6cOIHMzEycOHEC8fHxyM3NhZOTE1q0aIGTJ08ao59ERFaNky4S2SaDEyGJRAKJRIJnnnlGua+4uBiZmZlMgojIrihWrU/LSdNYM8REiMi6GGWJjYd5e3ujV69e6NWrlymaJyKyaoqaoarJECddJLJOehVLX758Wafjr169qs9liIhskrqFWqtOusgCaiLroVci1KVLF0yYMKHGtcSKioqwfv16tG3bFt9//73eHSQiskVxneKQOyMXabFpyJ2Ri7hOcdiQsQGSlRKuWk9kRfS6NZaVlYUPPvgAUVFRcHNzQ+fOndG4cWO4ubnhzp07yMrKwpkzZ9CpUyd89NFHGDhwoLH7TURk9RQ1QwA0TrrIAmoiy9JrRMjHxwfLly9Hfn4+Pv/8czRv3hw3b95ULq768ssv49ixY0hPT2cSREQETrpIZK0MKpZ2d3fHCy+8gBdeeMFY/SEiqpNqKqDmYq1ElmPwzNIAcPLkSaxevRrr169HVlaWMZokIqpTNBVQ776wm3VDRBZk8OPzq1atwsyZMyEWi+Ho6Ig7d+6gXbt22Lx5Mzp06GCELhIR1Q0Pr1oPAJKVEtYNEVmQXiNCGzduREZGBioqKvDBBx9g6dKluHPnDm7duoW///4bAwYMQK9evXDw4EFj95eIyKZVXbWedUNElqfXiNDHH3+sLIyWy+U4cuQIVq1ahY4dO6JDhw5YunQpgoOD8cYbbzAZIiLSgBMvElmeXiNCWVlZKC4uxsGDB+Hs7AwHBwds3boVAwcORMOGDdG0aVMkJyfj2LFj+O9//4vc3Fwjd5uIyPZpqhsCwEkXicxE7xohNzc3dOnSBY8//jjCw8Px7bffQi6X49y5c8jMzMT+/fvx+++/Y9SoUbhz5w68vLwglUqN2XciIpv3cN2Qoniaq9YTmYfBxdKffPIJIiMj8ffff2PixIkIDw9HcHAwMjIy0LhxY+Tl5SEvLw+nT582Rn+JiOocxcSLnHSRyPwMToQ6dOiAY8eOYeLEiejevTsEQXjQsJMTNm7cCAAICgpCUBB/iYmIalJT8TQTISLTMMrq882aNUNKSgoKCwtx6NAh3Lt3DxEREUx+iIh0oM2kiyHiEMt1kKgOMkoipODv749nn33WmE0SEdkNRfH0hJ0TUClUqky6qLhl5iBywKSgSRgILl9EZAxGTYSIiMgw2ky6uObKGrwufR2hPqGW7CpRncBEiIjIylRdtT4tJ61a3ZAccly8c5GJEJERGGWtMSIiMg1F3VBVDnBAswbNLNQjorrF4EQoLS1N42vr1q0ztHkiIrumbtLFScGTAHDSRSJjMDgReuqppzB79mzIZDLlvps3b2Lw4MGYO3euoc0TEdm9uE5xyJ2Ri7TYNGRPebC8UdjqMK5YT2QERhkRSk5ORpcuXZCVlYX//ve/aNu2LaRSKTIzM43QRSIiUizWCgBfXPmi2qSLHBki0o/BiVCPHj2QmZmJtm3bolOnThgyZAhmzpyJvXv3QiKRGKOPRET0/y7cvgABgso+rlhPpD+jFEv/9ddfOHr0KIKCguDk5ITz58+jrKzMGE0TEVEVYQ3DIIJIZV/VSRdZN0SkG4MToaVLlyIiIgJRUVE4ffo0/vzzTxw/fhzt27dHenq6MfoIANi/fz8GDx6Mxo0bQyQSYceOHTUev3fvXohEomo/BQUFRusTEZG5BYmDMDl4crUV6xWLtbJuiEg3Bs8jtGrVKuzYsQMDBgwAALRt2xZ//vkn3nrrLURGRqKiosLgTgJAaWkpwsPDMXbsWDz33HNan3f+/HmIxWLldqNGjYzSHyIiS4nyicLrz76OS8WXNE66yMVaibRjcCJ06tQp+Pr6quxzdnbGsmXLMGjQIEObVxowYIAy2dJFo0aNUL9+faP1g4jIGgSJg5QTKqqbdJGLtRJpx+BE6OEkqKrevXsb2rzBOnTogIqKCrRt2xYLFy7E448/rvHYiooKlREsqVQKAJDJZCrTAxiDoj1jt0uqGGfzYJzNR12sQ8QhahdrlXhL+Jnoid9p8zBVnHVpTyQIglD7YZq99957Nb4+f/58Q5pXSyQSITk5GTExMRqPOX/+PPbu3YvHHnsMFRUV+PLLL7FlyxYcPnwYnTp1UnvOwoULsWjRomr7k5KS4OHhYazuExEZXcqtFKy5sgZyyOEAB0wKnoSO3h2RX5GPQNdA+Lpo/kcrUV1TVlaGl156CUVFRSrlMeoYnAh17NhRZVsmkyEnJwdOTk5o1qwZMjIyDGleLW0SIXV69+6NRx55BFu2bFH7uroRoeDgYNy8ebPWQOpKJpMhJSUFUVFRcHZ2Nmrb9C/G2TwYZ/OpKdZ50jxcvHMRzRo0Q8rfKZi0a5Jyxfo1A9ZgTIcxFuq17eF32jxMFWepVApfX1+tEiGDb40dP35cbQdGjx6NIUOGGNq8UXXt2hUHDhzQ+LqrqytcXV2r7Xd2djbZL4Ip26Z/Mc7mwTibj7pYh/qEItQnFHnSPGUSBDwonp68azIGthjImiEd8TttHsaOsy5tmWTRVbFYjEWLFuHdd981RfN6y8zMRGBgoKW7QURkUtm3sjUWTxORKoNHhDQpKipCUVGR0dorKSnBhQv//hLn5OQgMzMTDRs2xCOPPIJ58+bh6tWr+OqrrwAAK1euRGhoKNq0aYPy8nJ8+eWX+P3337Fnzx6j9YmIyBopVqx/uHha8ag9Ef3L4ETo008/VdkWBAH5+fnYsmWLXo+7a3L06FE8+eSTyu1Zs2YBAGJjY5GYmIj8/HxcvnxZ+fq9e/fw+uuv4+rVq/Dw8ED79u3x22+/qbRBRFQXKVasn7BzAiqFSuWki0HiIORJ85B9KxvNfZrzNhkRjJAIrVixQmXbwcEBfn5+iI2Nxbx58wxtXikyMhI11XUnJiaqbM+ZMwdz5swx2vWJiGxJXKc4RIdF48LtCwhrGIYgcRA2ZGzA+J3jlQXUCYMSENcpztJdJbIogxOhnJwcY/SDiIiMLEgcpBz1yZPmKZMggLNPEymYpFiaiIisCwuoidTTa0RIUZ+jjeXLl+tzCSIiMiJNBdSezp5Iy0ljzRDZLb0SIXVzB6kjEon0aZ6IiIxMXQH1yPYj0X1Dd9YMkV3TKxFKS0szdj+IiMjEqhZQezp7KpMggDVDZL/0rhH6+++/a3yKi4iIrE+QOAiRIZEouVfCmiEiGJAINW/eHDdu3FBuDxs2DIWFhUbpFBERmZaiZqgqTrpI9kjvROjh0aBffvkFpaWlBneIiIhMT1Ez5ChyBIBqky6m5aQhT5pn4V4SmZ7JltggIiLrxkkXiQwYERKJRNWeCuNTYkREtkVRM6QYCVI36SJHhqgu03tESBAEjB49Gq6urgCA8vJyTJw4EZ6enirH/fDDD4b1kIiIzKKmSRf5JBnVVXonQrGxsSrbI0eONLgzRERkOVy1nuyR3onQpk2bjNkPIiKyMK5aT/aIxdJERKTEAmqyN1x0lYiIVLCAmuwJEyEiItKIq9ZTXcdEiIiINNI0A7Vi1XqODJGtYyJEREQaqZuBWrFqfZ+v+kCyUoINGRss3Esi/RmlWDo1NRWpqam4fv065HLVIdSNGzca4xJERGQhXLWe6jKDR4QWLVqE/v37IzU1FTdv3sSdO3dUfoiIyPZx1XqqqwweEVq7di0SExPxyiuvGKM/RERkxTjpItU1Bo8I3bt3Dz169DBGX4iIyMpx1XqqawweEXr11VeRlJSEd9991xj9ISIiK8dJF6kuMTgRKi8vR0JCAn777Te0b98ezs7OKq8vX77c0EsQEZGVCRIHKYujNU26yAJqsgUGJ0InT55Ehw4dAACnT59WeU0kEhnaPBERWTmuWk+2zOBEKC0tzRj9ICIiG8UCarJlnFCRiIgMoqmAGgCLp8nqGWVCxbt372LDhg04e/YsAKB169aIi4tDvXr1jNE8ERFZuYcLqHdf2A3JSgmLp8nqGTwidPToUTRr1gwrVqzA7du3cfv2baxYsQLNmjVDRkaGMfpIREQ2QDHpIgCuWE82w+ARoZkzZ+KZZ57B+vXr4eT0oLn79+/j1VdfxYwZM7B//36DO0lERLaDxdNkSwxOhI4ePaqSBAGAk5MT5syZg8cee8zQ5omIyMbUVDydJ81D9q1sNPdpzqSIrILBt8bEYjEuX75cbf+VK1fg7e1taPNERGRjNBVPK+qGuGo9WRODR4SGDRuGuLg4fPzxx8qlNv744w/Mnj0bI0aMMLiDRERkex4ungagLJ4GOOkiWQ+DE6GPP/4YIpEIo0aNwv379wEAzs7OmDRpEpYuXWpwB4mIyDZVnX06LSeNdUNklQxKhGQyGQYMGIC1a9ciPj4eFy9eBAA0a9YMHh4eRukgERHZPk66SNbKoBohZ2dnnDx5EgDg4eGBdu3aoV27dkyCiIhIBSddJGtlcLH0yJEjsWEDC96IiKhmcZ3ikDsjF2mxacidkQsALJ4mizO4Ruj+/fvYuHEjfvvtN3Tu3Bmenp4qr3P1eSIiUlDUDXHFerIWBidCp0+fRqdOnQAAf/31l8prXH2eiIjU4aSLZC24+jwREZkdJ10ka2FwjdDly5chCILG14iIiB7GSRfJWhg8IhQaGor8/Hw0atRIZf+tW7cQGhqKyspKQy9BRER1ECddJGtgcCIkCILaWqCSkhK4ubkZ2jwREdVhnHSRLE3vRGjWrFkAHhREv/vuuypzB1VWVuLw4cPo0KGDwR1U2L9/P5YtW4Zjx44hPz8fycnJiImJqfGcvXv3YtasWThz5gyCg4PxzjvvYPTo0UbrExERGQ8nXSRL0LtG6Pjx4zh+/DgEQcCpU6eU28ePH8e5c+cQHh6OxMREo3W0tLQU4eHhWL16tVbH5+Tk4Omnn8aTTz6JzMxMzJgxA6+++ip2795ttD4REZHxcNJFsgS9R4QUT4uNGTMGq1atglgsNlqn1BkwYAAGDBig9fFr165FaGgoPvnkEwBAq1atcODAAaxYsQLR0dGm6iYRERng4bohRfG0XJDDQeSAhEEJiOsUZ+luUh1icI3Qpk2bjNEPo0tPT0e/fv1U9kVHR2PGjBkaz6moqEBFRYVyWyqVAniwpppMJjNq/xTtGbtdUsU4mwfjbD72EGt/d3/4N/HXOOliH0kfk9cM2UOcrYGp4qxLewYnQu+9916Nr8+fP9/QS+iloKAA/v7+Kvv8/f0hlUrxzz//wN3dvdo58fHxWLRoUbX9e/bsMdn6aSkpKSZpl1QxzubBOJuPPcT6VPEptcXTX+/6Gu2825mlD/YQZ2tg7DiXlZVpfazBiVBycrLKtkwmQ05ODpycnNCsWTOLJUL6mDdvnrIIHHgwIhQcHIz+/fsb/dafTCZDSkoKoqKi4OzsbNS26V+Ms3kwzuZjT7FuL22PBasXVCuefnnAy2YZEbKXOFuSqeKsuKOjDYMToePHj6vtwOjRozFkyBBDm9dbQEAACgsLVfYVFhZCLBarHQ0CAFdXV7i6ulbb7+zsbLJfBFO2Tf9inM2DcTYfe4h1qE8oEgYlYMLOCagUKpXF06E+oWabfdoe4mwNjB1nXdoyOBFSRywWY9GiRRg8eDBeeeUVU1yiVhEREfjll19U9qWkpCAiIsIi/SEiIt09XDwdJA7ChowNytohFlCToQxeYkOToqIiFBUVGa29kpISZGZmIjMzE8CDx+MzMzOVy3jMmzcPo0aNUh4/ceJE/P3335gzZw7OnTuHL774At999x1mzpxptD4REZHpBYmDEBkSWeOq9Xy0nvRl8IjQp59+qrItCALy8/OxZcsWnR53r83Ro0fx5JNPKrcVtTyxsbFITExEfn6+ytpmoaGh+O9//4uZM2di1apVCAoKwpdffslH54mIbBhXrSdjMzgRWrFihcq2g4MD/Pz8EBsbi3nz5hnavFJkZKTGxV0BqJ28MTIyUm0NExER2SauWk/GZnAilJOTY4x+EBER1Uox+/TDBdS7L+xm3RDpxSjF0v/73/+wbt06/P3339i2bRuaNGmCLVu2IDQ0FD179jTGJYiIiABw1XoyLoOLpb///ntER0fD3d0dGRkZypmZi4qKsGTJEoM7SERE9LCqBdQ11Q0R1cbgROj999/H2rVrsX79epXn9h9//HFkZGQY2jwREVGNFHVDVXHVetKWwYnQ+fPn8cQTT1TbX69ePdy9e9fQ5omIiGrEVevJEAbXCAUEBODChQsICQlR2X/gwAE0bdrU0OaJiIhqxVXrSV8GjwiNGzcO06dPx+HDhyESiXDt2jV8/fXXeOONNzBp0iRj9JGIiKhWirohAJx0kbRm8IjQ3LlzIZfL0bdvX5SVleGJJ56Aq6sr3njjDUybNs0YfSQiItIaJ10kXRicCIlEIrz99tuYPXs2Lly4gJKSErRu3RpeXl7G6B8REZFOOOki6cJoa425uLigdevW6Nq1K5MgIiKyGE3F04q6oT5f9YFkpQQbMjZYuKdkDfQeEXJwcIBIJKrxGJFIhPv37+t7CSIiIr1w0kXSlt6JUHJyssbX0tPT8emnn0Iul2s8hoiIyJSCxEHKJCctJ411Q6SW3onQs88+W23f+fPnMXfuXPz88894+eWX8d577xnUOSIiImOoqW6I7JtRaoSuXbuGcePGoV27drh//z4yMzOxefNmSCQSYzRPRERkEE66SJoY9NSYYj2xzz77DB06dEBqaip69eplrL4REREZDSddJHX0HhH66KOP0LRpU+zcuRPffPMNDh48yCSIiIisGiddpIfpPSI0d+5cuLu7IywsDJs3b8bmzZvVHvfDDz/o3TkiIiJT4KSLpKB3IjRq1KhaH58nIiKyRrVNuni28Cxu3rtpwR6SueidCCUmJhqxG0REROajKJ6esHMCKoVKlUkXFbfMRBChMqQS47uMt3R3yYQMXmKDiIjIFtU26aIAAZN3TcbAFgN5u6wOYyJERER2i5MuktHWGiMiIrJlirqhqjjpYt3HRIiIiAjVJ110gAO+GPAFAE66WJfx1hgREdH/U9QNnbt+DpeOX0IlKjnpYh3HESEiIqIqgsRB6C3pDQCYtGsSJ12s45gIERERqZFfka+xeJrqDiZCREREagS6Bmosns6T5rFuqI5gIkRERKSGr4sv1gxYU23FesVirX2+6gPJSgk2ZGywcE/JECyWJiIi0mBMhzEY2GKgxkkXFXVD0WHRnGvIRjERIiIiqgEnXazbeGuMiIhIS5omXfR09mTNkI1iIkRERKSlhydddBQ5YmT7kei+oTtrhmwUb40RERHpoOpirZ7Onui+oTtrhmwYR4SIiIh0FCQOQmRIJErulXCuIRvHRIiIiEhPXKjV9jERIiIi0pO6mqF1g9YhSBzESRdtBGuEiIiIDFC1ZiisYRiCxEHYkLEB43eO52KtNoAjQkRERAZS1AwpRoIUSRDAxVqtHRMhIiIiI8q+lc0CahvCRIiIiMiIaiqgZt2Q9WEiREREZESaCqi5WKt1YrE0ERGRkT1cQA1wsVZrZXMjQqtXr0ZISAjc3NzQrVs3/PnnnxqPTUxMhEgkUvlxc3MzY2+JiMheVS2gZt2Q9bKpROjbb7/FrFmzsGDBAmRkZCA8PBzR0dG4fv26xnPEYjHy8/OVP5cuXTJjj4mIiLhYqzWzqURo+fLlGDduHMaMGYPWrVtj7dq18PDwwMaNGzWeIxKJEBAQoPzx9/c3Y4+JiIi4WKs1s5kaoXv37uHYsWOYN2+ecp+DgwP69euH9PR0jeeVlJRAIpFALpejU6dOWLJkCdq0aaP22IqKClRUVCi3pVIpAEAmk0EmkxnpnUDZZtX/kmkwzubBOJsPY20epojzqHaj0EfSBxfvXISHswd6be5VrWaoj6SPXdUMmer7rEt7NpMI3bx5E5WVldVGdPz9/XHu3Dm157Ro0QIbN25E+/btUVRUhI8//hg9evTAmTNnEBRU/YsWHx+PRYsWVdu/Z88eeHh4GOeNPCQlJcUk7ZIqxtk8GGfzYazNw1RxPlR8SG3N0Ne7vkY773YmuaY1M3acy8rKtD5WJAiCYNSrm8i1a9fQpEkTHDx4EBEREcr9c+bMwb59+3D48OFa25DJZGjVqhVGjBiBxYsXV3td3YhQcHAwbt68CbFYbJw3UqUvKSkpiIqKgrOzs1Hbpn8xzubBOJsPY20epo5znjQPYavDVJIhR5EjsqdkA4DKch11maniLJVK4evri6Kiolr//raZESFfX184OjqisLBQZX9hYSECAgK0asPZ2RkdO3bEhQvqq/RdXV3h6uqq9jxT/YFjyrbpX4yzeTDO5sNYm4ep4hzqE4qEQQmYsHMCKoVK5VxDv1/63S7XKDN2nHVpy2aKpV1cXNC5c2ekpqYq98nlcqSmpqqMENWksrISp06dQmBgoKm6SUREpJW4TnHInZGLtNg05M7IRXRYNNcoswCbGRECgFmzZiE2NhaPPfYYunbtipUrV6K0tBRjxowBAIwaNQpNmjRBfHw8AOC9995D9+7dERYWhrt372LZsmW4dOkSXn31VUu+DSIiIgAPniZT3P5Ky0nTONdQXb9FZkk2lQgNGzYMN27cwPz581FQUIAOHTrg119/VRZQX758GQ4O/w5y3blzB+PGjUNBQQEaNGiAzp074+DBg2jdurWl3gIREZFairmGHq4bUsw11NynORMiE7CpRAgApk6diqlTp6p9be/evSrbK1aswIoVK8zQKyIiIsMo5hqqWjekmGvI3mqGzMnmEiEiIqK6quoaZZ7OnsokCOD6ZKZiM8XSRERE9kCxRlnJvRKuT2YGTISIiIiskKb1ycIahiFPmsc1yoyEiRAREZEVUrc+2bpB67D7wm5IVkq4RpmRsEaIiIjISlWtGQprGAYAkKyUsG7IiJgIERERWTHONWRavDVGRERkIzTVDSnmGmLNkO6YCBEREdkIdXVDirmGWDOkH94aIyIisiGca8i4OCJERERkYzjXkPEwESIiIrJRnGvIcEyEiIiIbBTnGjIca4SIiIhsGOcaMgwTISIiIhvHuYb0x1tjREREdQjnGtINEyEiIqI6hHMN6Ya3xoiIiOoYzjWkPY4IERER1UGca0g7TISIiIjqMM41VDMmQkRERHUY5xqqGWuEiIiI6jjONaQZEyEiIiI7wLmG1OOtMSIiIjvDuqF/MREiIiKyM6wb+hdvjREREdkh1g09wESIiIjITrFuiLfGiIiICPa7RhkTISIiIrLbNcp4a4yIiIgA2OcaZRwRIiIiIiV7W6OMiRARERFVYy9zDTERIiIiomrsZa4h1ggRERGRWvYw1xATISIiItJIm7mG0q+kw9fDF819mttcQsREiIiIiLSiqBuqmgyJIMLw74dDLsjhIHJAwqAExHWKs2AvdcMaISIiItLKw3VDDv+fRjx8q8yWiqg5IkRERERaq1o3dL30OoZtH6byuq0ty8FEiIiIiHSiqBvKk+ZVu1VW9RH77FvZVl83xFtjREREpJe68Ig9R4SIiIhIb7b+iD0TISIiIjKILT9iz0SIiIiIjMbWHrFnjRAREREZja09Ym9zidDq1asREhICNzc3dOvWDX/++WeNx2/btg0tW7aEm5sb2rVrh19++cVMPSUiIrJPcZ3ikDsjF2mxafjmhW8gQFB53ZpWsbepROjbb7/FrFmzsGDBAmRkZCA8PBzR0dG4fv262uMPHjyIESNGIC4uDsePH0dMTAxiYmJw+vRpM/eciIjIvgSJgxAZEokewT1qXMX+VPEpi44O2VQitHz5cowbNw5jxoxB69atsXbtWnh4eGDjxo1qj1+1ahWeeuopzJ49G61atcLixYvRqVMnfP7552buORERkX2q6RH7sNVhePfiuwhbHWaxR+xtplj63r17OHbsGObNm6fc5+DggH79+iE9PV3tOenp6Zg1a5bKvujoaOzYsUPt8RUVFaioqFBuS6VSAIBMJoNMJjPwHahStGfsdkkV42wejLP5MNbmwTgb16h2o9BH0gcX71xEswbNAABhq8Oq1Q31kfQxyhNlunxuNpMI3bx5E5WVlfD391fZ7+/vj3Pnzqk9p6CgQO3xBQUFao+Pj4/HokWLqu3fs2cPPDw89Ox5zVJSUkzSLqlinM2DcTYfxto8GGfjO4mTOFV8Su0j9l/v+hrtvNsZfI2ysjKtj7WZRMgc5s2bpzKCJJVKERwcjP79+0MsFhv1WjKZDCkpKYiKioKzs7NR26Z/Mc7mwTibD2NtHoyzabWXtseC1QuqLc3x8oCXjTIipLijow2bSYR8fX3h6OiIwsJClf2FhYUICAhQe05AQIBOx7u6usLV1bXafmdnZ5P9IpiybfoX42wejLP5MNbmwTibRqhPKBIGJWDCzgmoFCqVdUOhPqFGaV+Xz8xmiqVdXFzQuXNnpKamKvfJ5XKkpqYiIiJC7TkREREqxwMPhjk1HU9ERETmEdcpDtlTsrG42WJkT8m22ASLNjMiBACzZs1CbGwsHnvsMXTt2hUrV65EaWkpxowZAwAYNWoUmjRpgvj4eADA9OnT0bt3b3zyySd4+umnsXXrVhw9ehQJCQmWfBtERESEB0+UtfNuZ9ElN2wqERo2bBhu3LiB+fPno6CgAB06dMCvv/6qLIi+fPkyHBz+HeTq0aMHkpKS8M477+Ctt95C8+bNsWPHDrRt29ZSb4GIiIisiE0lQgAwdepUTJ06Ve1re/furbZv6NChGDp0qIl7RURERLbIZmqEiIiIiIyNiRARERHZLSZCREREZLeYCBEREZHdYiJEREREdouJEBEREdktJkJERERkt5gIERERkd2yuQkVzUkQBAC6rWKrLZlMhrKyMkilUi7oZ0KMs3kwzubDWJsH42wepoqz4u9txd/jNWEiVIPi4mIAQHBwsIV7QkRERLoqLi5GvXr1ajxGJGiTLtkpuVyOa9euwdvbGyKRyKhtS6VSBAcH48qVKxCLxUZtm/7FOJsH42w+jLV5MM7mYao4C4KA4uJiNG7cWGUNUnU4IlQDBwcHBAWZdkVcsVjMXzIzYJzNg3E2H8baPBhn8zBFnGsbCVJgsTQRERHZLSZCREREZLeYCFmIq6srFixYAFdXV0t3pU5jnM2DcTYfxto8GGfzsIY4s1iaiIiI7BZHhIiIiMhuMREiIiIiu8VEiIiIiOwWEyEiIiKyW0yEiIiIyG4xEbKA1atXIyQkBG5ubujWrRv+/PNPS3fJ5u3fvx+DBw9G48aNIRKJsGPHDpXXBUHA/PnzERgYCHd3d/Tr1w/Z2dmW6awNi4+PR5cuXeDt7Y1GjRohJiYG58+fVzmmvLwcU6ZMgY+PD7y8vPD888+jsLDQQj22TWvWrEH79u2Vs+1GRERg165dytcZY9NYunQpRCIRZsyYodzHWBtu4cKFEIlEKj8tW7ZUvm7pGDMRMrNvv/0Ws2bNwoIFC5CRkYHw8HBER0fj+vXrlu6aTSstLUV4eDhWr16t9vWPPvoIn376KdauXYvDhw/D09MT0dHRKC8vN3NPbdu+ffswZcoUHDp0CCkpKZDJZOjfvz9KS0uVx8ycORM///wztm3bhn379uHatWt47rnnLNhr2xMUFISlS5fi2LFjOHr0KPr06YNnn30WZ86cAcAYm8KRI0ewbt06tG/fXmU/Y20cbdq0QX5+vvLnwIEDytcsHmOBzKpr167ClClTlNuVlZVC48aNhfj4eAv2qm4BICQnJyu35XK5EBAQICxbtky57+7du4Krq6vwzTffWKCHdcf169cFAMK+ffsEQXgQV2dnZ2Hbtm3KY86ePSsAENLT0y3VzTqhQYMGwpdffskYm0BxcbHQvHlzISUlRejdu7cwffp0QRD4fTaWBQsWCOHh4Wpfs4YYc0TIjO7du4djx46hX79+yn0ODg7o168f0tPTLdizui0nJwcFBQUqca9Xrx66devGuBuoqKgIANCwYUMAwLFjxyCTyVRi3bJlSzzyyCOMtZ4qKyuxdetWlJaWIiIigjE2gSlTpuDpp59WiSnA77MxZWdno3HjxmjatClefvllXL58GYB1xJirz5vRzZs3UVlZCX9/f5X9/v7+OHfunIV6VfcVFBQAgNq4K14j3cnlcsyYMQOPP/442rZtC+BBrF1cXFC/fn2VYxlr3Z06dQoREREoLy+Hl5cXkpOT0bp1a2RmZjLGRrR161ZkZGTgyJEj1V7j99k4unXrhsTERLRo0QL5+flYtGgRevXqhdOnT1tFjJkIEZFepkyZgtOnT6vc6yfjadGiBTIzM1FUVITt27cjNjYW+/bts3S36pQrV65g+vTpSElJgZubm6W7U2cNGDBA+f/t27dHt27dIJFI8N1338Hd3d2CPXuAt8bMyNfXF46OjtWq4QsLCxEQEGChXtV9itgy7sYzdepU7Ny5E2lpaQgKClLuDwgIwL1793D37l2V4xlr3bm4uCAsLAydO3dGfHw8wsPDsWrVKsbYiI4dO4br16+jU6dOcHJygpOTE/bt24dPP/0UTk5O8Pf3Z6xNoH79+nj00Udx4cIFq/g+MxEyIxcXF3Tu3BmpqanKfXK5HKmpqYiIiLBgz+q20NBQBAQEqMRdKpXi8OHDjLuOBEHA1KlTkZycjN9//x2hoaEqr3fu3BnOzs4qsT5//jwuX77MWBtILpejoqKCMTaivn374tSpU8jMzFT+PPbYY3j55ZeV/89YG19JSQkuXryIwMBA6/g+m6Ukm5S2bt0quLq6ComJiUJWVpYwfvx4oX79+kJBQYGlu2bTiouLhePHjwvHjx8XAAjLly8Xjh8/Lly6dEkQBEFYunSpUL9+feHHH38UTp48KTz77LNCaGio8M8//1i457Zl0qRJQr169YS9e/cK+fn5yp+ysjLlMRMnThQeeeQR4ffffxeOHj0qRERECBERERbste2ZO3eusG/fPiEnJ0c4efKkMHfuXEEkEgl79uwRBIExNqWqT40JAmNtDK+//rqwd+9eIScnR/jjjz+Efv36Cb6+vsL169cFQbB8jJkIWcBnn30mPPLII4KLi4vQtWtX4dChQ5buks1LS0sTAFT7iY2NFQThwSP07777ruDv7y+4uroKffv2Fc6fP2/ZTtsgdTEGIGzatEl5zD///CNMnjxZaNCggeDh4SEMGTJEyM/Pt1ynbdDYsWMFiUQiuLi4CH5+fkLfvn2VSZAgMMam9HAixFgbbtiwYUJgYKDg4uIiNGnSRBg2bJhw4cIF5euWjrFIEATBPGNPRERERNaFNUJERERkt5gIERERkd1iIkRERER2i4kQERER2S0mQkRERGS3mAgRERGR3WIiRERERHaLiRARERHZLSZCRGT1IiMjMWPGjGr/b2mRkZEQiUQQiUTIzMw02XVGjx6tvM6OHTtMdh0ie8REiIiMoupf1s7OzggNDcWcOXNQXl5u1Ov88MMPWLx4sVHbNMS4ceOQn5+Ptm3bmuwaq1atQn5+vsnaJ7JnTpbuABHVHU899RQ2bdoEmUyGY8eOITY2FiKRCB9++KHRrtGwYUOjtaWL+/fvw8mp+h+ZHh4eCAgIMOm169Wrh3r16pn0GkT2iiNCRGQ0rq6uCAgIQHBwMGJiYtCvXz+kpKQoX//111/Rs2dP1K9fHz4+Phg0aBAuXryo0kZpaSlGjRoFLy8vBAYG4pNPPlF5/eFbYyEhIVi5cqXKMR06dMDChQuV29u3b0e7du3g7u4OHx8f9OvXD6WlpRrfR25uLkQiEb777jv06tULrq6u+Omnn7SKgVwux0cffYSwsDC4urrikUcewQcffKDS/2nTpmHGjBlo0KAB/P39sX79epSWlmLMmDHw9vZGWFgYdu3apdX1iMgwTISIyCROnz6NgwcPwsXFRbmvtLQUs2bNwtGjR5GamgoHBwcMGTIEcrlceczs2bOxb98+/Pjjj9izZw/27t2LjIwMvfuRn5+PESNGYOzYsTh79iz27t2L5557DjWtN33ixAkAwLJlyzB//nycOXMGffv21ep68+bNw9KlS/Huu+8iKysLSUlJ8Pf3Vzlm8+bN8PX1xZ9//olp06Zh0qRJGDp0KHr06IGMjAz0798fr7zyCsrKyvR+30SkHd4aIyKj2blzJ7y8vHD//n1UVFTAwcEBn3/+ufL1559/XuX4jRs3ws/PD1lZWWjbti1KSkqwYcMG/Oc//1EmHps3b0ZQUJDefcrPz8f9+/fx3HPPQSKRAADatWtX4zmZmZnw9PTEtm3bEBISovW1iouLsWrVKnz++eeIjY0FADRr1gw9e/ZUOS48PBzvvPMOgH8TJ19fX4wbNw4AMH/+fKxZswYnT55E9+7dtb4+EemOI0JEZDRPPvkkMjMzcfjwYcTGxmLMmDEqyU92djZGjBiBpk2bQiwWK5OMy5cvAwAuXryIe/fuoVu3bspzGjZsiBYtWujdp/DwcPTt2xft2rXD0KFDsX79ety5c6fGc06cOIFnnnlGpyQIAM6ePYuKiopaR4/at2+v/H9HR0f4+PioJGeKEaTr16/rdH0i0h0TISIyGk9PT4SFhSE8PBwbN27E4cOHsWHDBuXrgwcPxu3bt7F+/XocPnwYhw8fBgDcu3dP72s6ODhUu80lk8mU/+/o6IiUlBTs2rULrVu3xmeffYYWLVogJydHY5uZmZmIjIzUuS/u7u5aHefs7KyyrXjSruo2AJVbhkRkGkyEiMgkHBwc8NZbb+Gdd97BP//8g1u3buH8+fN455130LdvX7Rq1arayEyzZs3g7OysTJAA4M6dO/jrr780XsfPz0/l0XKpVFotyRGJRHj88cexaNEiHD9+HC4uLkhOTlbbnlQqRW5uLjp27Kjze27evDnc3d2Rmpqq87lEZBlMhIjIZIYOHQpHR0esXr0aDRo0gI+PDxISEnDhwgX8/vvvmDVrlsrxXl5eiIuLw+zZs/H777/j9OnTGD16NBwcNP9R1adPH2zZsgX/+9//cOrUKcTGxsLR0VH5+uHDh7FkyRIcPXoUly9fxg8//IAbN26gVatWats7ceIEHB0da60jUsfNzQ1vvvkm5syZg6+++goXL17EoUOHVEbFiMi6sFiaiEzGyckJU6dOxUcffYRJkyZh69ateO2119C2bVu0aNECn376abVbUMuWLUNJSQkGDx4Mb29vvP766ygqKtJ4jXnz5iEnJweDBg1CvXr1sHjxYpURIbFYjP3792PlypWQSqWQSCT45JNPMGDAALXtnThxAi1atICbm5te7/ndd9+Fk5MT5s+fj2vXriEwMBATJ07Uqy0iMj2RUNMzpEREpFFkZCQ6dOhQbR4jUxGJREhOTkZMTIxZrkdkD3hrjIjIAF988QW8vLxw6tQpk11j4sSJ8PLyMln7RPaMI0JERHq6evUq/vnnHwDAI488ojJ5pDFdv34dUqkUABAYGAhPT0+TXIfIHjERIiIiIrvFW2NERERkt5gIERERkd1iIkRERER2i4kQERER2S0mQkRERGS3mAgRERGR3WIiRERERHaLiRARERHZLSZCREREZLeYCBEREZHd+j/5jUkXKr+WtgAAAABJRU5ErkJggg==",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"neutron_flux = lambda r: 5e13 * np.sin(0.0628*r)/r #neutrons/(cm^2-s)\n",
"\n",
"r_range = np.linspace(1e-5, 50, 100)\n",
"flux_range = neutron_flux(r_range)\n",
"\n",
"plt.figure(); plt.title(\"Neutron Flux as a Function of Radius $r$\")\n",
"plt.xlabel('Radius $r$ [cm]'); plt.ylabel('Neutron Flux $\\\\phi(r)$ $\\\\frac{\\\\text{neutrons}}{\\\\text{cm}^2 \\\\times \\\\text{sec}}$')\n",
"plt.grid('both')\n",
"plt.plot(r_range, flux_range, '.g')"
]
},
{
"cell_type": "markdown",
"id": "17b06d8d-fc16-47f0-b5c0-8bbf49198a17",
"metadata": {},
"source": [
"Flux technically isn't defined at the center of the reactor, but neutron flux is greatest as r approaches 0."
]
},
{
"cell_type": "markdown",
"id": "080f3f81-dc35-46fd-8cf8-dda1749e5fa9",
"metadata": {},
"source": [
"## Calculate the neutron current density as a function of position in the reactor."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "263c4561-9a00-44e1-8732-a442979aed5b",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle J = - \\frac{0.05024 a \\cos{\\left(0.0628 r \\right)}}{r} + \\frac{0.8 a \\sin{\\left(0.0628 r \\right)}}{r^{2}}$"
],
"text/plain": [
"Eq(J, -0.05024*a*cos(0.0628*r)/r + 0.8*a*sin(0.0628*r)/r**2)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"'neutrons/cm^2/s'"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"J, r, a = sm.symbols('J, r, a')\n",
"D = 0.80 #cm\n",
"\n",
"psi = a*sm.sin(0.0628*r)/r\n",
"\n",
"#Phi only depends on r.\n",
"#Therefore, we can take the gradient with only the partial w.r.t. r.\n",
"diffusion_approx = sm.Eq(J, -D*psi.diff(r))\n",
"display(diffusion_approx, 'neutrons/cm^2/s')"
]
},
{
"cell_type": "markdown",
"id": "b5a31a40-dd14-40b0-ad0d-af8a4de4baeb",
"metadata": {},
"source": [
"Where $a = 5\\times 10^{13}$ and J is in the $\\hat e_r$ direction."
]
},
{
"cell_type": "markdown",
"id": "1e2203b1-19bf-4348-a2fb-c59410abde46",
"metadata": {},
"source": [
"## How many neutrons escape from the reactor per second?\n",
"The neutrons escaping from the reactor per second is equal to the neutron current density multiplied by the area:\n",
"\n",
"$$ \\int_A \\left( J(\\vec r, t) \\cdot n \\right) dA $$\n",
"\n",
"We can see that our neutron current density is always normal to the surface, and our surface is at a constant value of r. This lets us simplify our expression:\n",
"\n",
"$$ \\int_A \\left( J(\\vec r, t) \\cdot n \\right) dA \\rightarrow \n",
"A J(r,t)|_{r = 50 \\text{[cm]}}$$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "8f13df8c-04c4-4833-b12e-05d24d0dd2d0",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"There are 1.579e+15 neutrons leaking per second from the reactor.\n"
]
}
],
"source": [
"a = 5e13\n",
"J = lambda r: -0.05024*a*np.cos(0.0628*r)/50 + 0.8*a*np.sin(0.0628*r)/r**2\n",
"A = lambda r: 4*np.pi*r**2\n",
"\n",
"print(f'There are {J(50)*A(50):.3e} neutrons leaking per second from the reactor.')"
]
},
{
"cell_type": "markdown",
"id": "45b572f6-2cfd-43d3-a2ee-6effafd2cbc2",
"metadata": {},
"source": [
"# Problem 2\n",
"A bare slab of thickness 2a contains uniformly distributed sources emitting Q neutrons/cm^3/sec. Given the expression for the scalar flux below verify the equation of continuity by computing per unit are of the slab the total number of neutrons\n",
"\n",
"$$ \\phi(x) = \\frac{Q}{\\Sigma_a} \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L) } \\right) $$\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "efef4bcb-e47b-4d6a-8a84-54d6c934b894",
"metadata": {},
"outputs": [],
"source": [
"Q, sigma, x, a, L = sm.symbols('Q, Sigma_a, x, a, L')\n",
"\n",
"phi = Q/sigma*(1-sm.cosh(x/L)/sm.cosh(a/L))"
]
},
{
"cell_type": "markdown",
"id": "906a04a7-def0-4eba-a8ea-d9859fa5b813",
"metadata": {},
"source": [
"## Produced within the whole slab"
]
},
{
"cell_type": "markdown",
"id": "a2ebdb98-99cf-4413-97fd-e1bcc8855322",
"metadata": {},
"source": [
"$$ \n",
"\\int_V s(\\vec r, t) dV = \\int_V Q dV \n",
"$$\n",
"$$\n",
"= \\int_A 2aQ dA\n",
"$$\n",
"but assuming per unit area...\n",
"$$\n",
"= \\frac{d}{dA} \\int_A 2aQ dA\n",
"$$\n",
"$$\\int_V s(\\vec r, t) dV = 2aQ \\text{ per second per unit area}$$"
]
},
{
"cell_type": "markdown",
"id": "4687ccaf-dbce-4474-91a4-7fb1412298b0",
"metadata": {},
"source": [
"## Absorbed per second within the slab"
]
},
{
"cell_type": "markdown",
"id": "777059b4-e89d-4062-bd1e-9781067682e3",
"metadata": {},
"source": [
"$$\n",
"\\int_v \\Sigma_a(\\vec r, t) \\phi(\\vec r, t) dV = \n",
"\\int_v \\Sigma_a(\\vec r, t) \\frac{Q}{\\Sigma_a} \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L)} \\right) dV\n",
"$$\n",
"$$\n",
"= \\int_{A} \\int_x Q \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L)} \\right) dxdA\n",
"$$\n",
"and doing things per unit area...\n",
"$$\n",
"\\frac{d}{dA} \\left(\\int_v \\Sigma_a(\\vec r, t) \\phi(\\vec r, t) dV \\right) = \\int_x Q \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L)} \\right) dx\n",
"$$\n",
"Now we can do this integration for the thickness x using SymPy:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "71824f97-740c-4ca8-a857-cffd085d3aa6",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle Q \\left(- \\frac{L \\sinh{\\left(\\frac{2 a}{L} \\right)}}{\\cosh{\\left(\\frac{a}{L} \\right)}} + 2 a\\right)$"
],
"text/plain": [
"Q*(-L*sinh(2*a/L)/cosh(a/L) + 2*a)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"absorption_per_second = sm.integrate(Q*(1-sm.cosh(x/L)/sm.cosh(a/L)), (x, 0, 2*a))\n",
"absorption_per_second"
]
},
{
"cell_type": "markdown",
"id": "c8a5504a-08b5-4743-bde5-650a325f5c67",
"metadata": {},
"source": [
"Which is our final expression for the absorption per second per unit area!"
]
},
{
"cell_type": "markdown",
"id": "577476c8-0bb8-46b9-896e-4681bb3bb05a",
"metadata": {},
"source": [
"## Escaping per second within the slab\n",
"$$\n",
"J = -D \\nabla \\phi(\\vec r, t) \n",
"$$\n",
"but since we're dealing with a large slab, all current is in the thickness (x) direction...\n",
"$$\n",
"= -D \\frac{d}{dx} \\left( \\phi(x) \\right)\n",
"$$\n",
"$$\n",
"= -\\Sigma_a L^2 \\frac{d}{dx} \\left(\\frac{Q}{\\Sigma_a} \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L)}\\right)\\right)\n",
"$$\n",
"\n",
"Now we know that the leakage rate is the neutron current out of each face. We can see by inspection that we won't have leakage out of the sides that aren't in the thickness direction. This is because we're assuming an infinitely large slab.\n",
"\n",
"$$\n",
"\\int_a \\left(J(x,t) \\cdot n \\right) dA = \\int_a \\left(-\\Sigma_a L^2 \\frac{d}{dx} \\left(\\frac{Q}{\\Sigma_a} \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L)}\\right)\\right) \\right)dA\n",
"$$\n",
"but once again, unit area...\n",
"$$\n",
"\\frac{d}{dA}\\left(\\int_a \\left(J(x,t) \\cdot n \\right) dA\\right) = \\frac{d}{dA}\\left(\\int_a \\left(-\\Sigma_a L^2 \\frac{d}{dx} \\left(\\frac{Q}{\\Sigma_a} \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L)}\\right)\\right) \\right)dA \\right)\n",
"$$\n",
"$$\n",
"J(x,t) =-\\Sigma_a L^2 \\frac{d}{dx} \\left(\\frac{Q}{\\Sigma_a} \\left( 1 - \\frac{\\cosh(x/L)}{\\cosh(a/L)}\\right)\\right) \n",
"$$\n",
"\n",
"And once again just doing the manipulation in SymPy"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "7ab93d22-cde3-4408-9ade-b9aeeab96ceb",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{L Q \\sinh{\\left(\\frac{x}{L} \\right)}}{\\cosh{\\left(\\frac{a}{L} \\right)}}$"
],
"text/plain": [
"L*Q*sinh(x/L)/cosh(a/L)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"J = -sigma*L**2*sm.diff(Q/sigma*(1-sm.cosh(x/L)/sm.cosh(a/L)),x)\n",
"J"
]
},
{
"cell_type": "markdown",
"id": "9dd7ea6e-dceb-44ec-bfd5-d183a2ee3f52",
"metadata": {},
"source": [
"But! There is leakage out of two faces: leakage at x = 0, and x = 2a. The total leakage is the addition of these two."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "f1a07559-c152-4185-923e-d3885ff63dcc",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 2 L Q \\sinh{\\left(\\frac{a}{L} \\right)}$"
],
"text/plain": [
"2*L*Q*sinh(a/L)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"total_leakage = J.subs(x,0) + J.subs(x,2*a)\n",
"sm.simplify(total_leakage)"
]
},
{
"cell_type": "markdown",
"id": "c9b8927d-ce50-4848-9772-d7e674f7c97a",
"metadata": {},
"source": [
"which is the total leakage of neutrons per unit area, per second."
]
},
{
"cell_type": "markdown",
"id": "ea81b9e6-fe62-4bad-ab6b-126d5d5b960e",
"metadata": {},
"source": [
"# Problem 3\n",
"Consider a critical bare slab reactor with the following properties: "
]
},
{
"cell_type": "markdown",
"id": "9c0f0535-d5e6-4f8c-b12f-e8499c073d9f",
"metadata": {},
"source": [
"## Calculate $k_{eff}$ for this reactor"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "cb4ff251-c19a-4f89-bcb5-b11f4969322e",
"metadata": {},
"outputs": [],
"source": [
"a = 40 #cm\n",
"Sigma_a = 0.066 #1/cm\n",
"D = 0.9\n",
"L = (D/Sigma_a)**(1/2)\n",
"nu = 2.6\n",
"nuSigma_f = 0.070 #1/cm\n",
"e = 200 #MeV / fission\n",
"L = (D/Sigma_a)**(1/2)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "d6fd9637-b0bc-43ec-a7ea-0013321db089",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.978318710166636"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"B_g = (3.1415/a)\n",
"k_eff = nuSigma_f/(D*B_g**2 + Sigma_a)\n",
"k_eff"
]
},
{
"cell_type": "markdown",
"id": "93c7fd3e-918d-44d7-a792-ab51ea9495ab",
"metadata": {},
"source": [
"## Calculate the width, a, that yields k_eff = 1"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "97ec39bf-c1e9-440c-a32e-cd7573bf524a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[-47.1210000000000, 47.1210000000000]"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"a = sm.symbols('a')\n",
"B_g = (3.1414/a)\n",
"k_eq = sm.Eq(1, nuSigma_f/(D*B_g**2 + Sigma_a))\n",
"sm.solve(k_eq)"
]
},
{
"cell_type": "markdown",
"id": "0a4c2ab9-174b-4af6-8f24-c417b5326a72",
"metadata": {},
"source": [
"The thickness that yields $k_{eff} = 1$ is 47.121 [cm]"
]
},
{
"cell_type": "markdown",
"id": "58847762-ba38-4764-9a0e-f89ea2ad3a4d",
"metadata": {},
"source": [
"## If the reactor is operating at 100 kW, calculate the flux distribution in the reactor, $\\Phi(x)$ using the width found in b."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "214e6142-4dee-4d85-b71f-1b0f20d1fa3c",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\Phi{\\left(x \\right)} = 1.22995582046794 \\cdot 10^{15} \\pi \\cos{\\left(0.0212219604847096 \\pi x \\right)}$"
],
"text/plain": [
"Eq(Phi(x), 1.22995582046794e+15*pi*cos(0.0212219604847096*pi*x))"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Phi = sm.Function('Phi')(x)\n",
"Sigma_f = nuSigma_f/nu\n",
"E = e*1.60218e-13 #J\n",
"P = 100e3 #W\n",
"a = 47.121 #cm\n",
"phi_eq = sm.Eq(Phi, P * sm.pi/2/a/E/Sigma_f * sm.cos(sm.pi*x/a))\n",
"phi_eq"
]
},
{
"cell_type": "markdown",
"id": "0fb10451-32b4-4cc7-b74d-3d1aefa4edb2",
"metadata": {},
"source": [
"Where x is in [cm]."
]
},
{
"cell_type": "markdown",
"id": "92b90596-1991-41f4-84f6-b46ad711a5f3",
"metadata": {},
"source": [
"## Calculate the power density as a function of x (Px)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "47a3f8c6-5fde-405c-8747-8f5591387a07",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\operatorname{Power Density}{\\left(x \\right)} = 1.22995582046794 \\cdot 10^{20} \\pi \\cos{\\left(0.0212219604847096 \\pi x \\right)}$"
],
"text/plain": [
"Eq(Power Density(x), 1.22995582046794e+20*pi*cos(0.0212219604847096*pi*x))"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Px = sm.Function('Power Density')(x)\n",
"Power_Density_eq = sm.Eq(Px,phi_eq.rhs*P)\n",
"Power_Density_eq"
]
},
{
"cell_type": "markdown",
"id": "0df4cd19-0dc2-4f0a-8490-3c9c33f273d8",
"metadata": {},
"source": [
"# Problem 4\n",
"Jezebel is a bare, fast, spherically shaped critical reactor constructed of pure $^{239}Pu$ metal ($\\rho = 15.4$ g/cm^3)\n",
"\n",
"## Calculate the critical radius and critical mass of the reactor using the following:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "61536d51-6bdd-40d2-a371-f58e3e9090d2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The critical radius is 1.912e+0 cm.\n",
"The critical mass is 2.926e+1 g.\n"
]
}
],
"source": [
"rho = 15.4 #g/cm^3\n",
"nu = 2.98\n",
"sigma_f = 1.85\n",
"sigma_a = 2.11\n",
"D = 1.26\n",
"R = sm.symbols('R')\n",
"B_g = (3.1415/R)\n",
"\n",
"critical_radius_eq = sm.Eq(1, nu*sigma_f/(D*B_g**2 + sigma_a))\n",
"critical_radius = sm.solve(critical_radius_eq)\n",
"print(f'The critical radius is {critical_radius[1]:.3e} cm.')\n",
"\n",
"critical_mass = 4/3*3.1415*critical_radius[1]**3\n",
"print(f'The critical mass is {critical_mass:.3e} g.')"
]
},
{
"cell_type": "markdown",
"id": "dc787a55-8bc4-4f20-8bd1-b12d12552b00",
"metadata": {},
"source": [
"## There has been considerable interest in the possibliilty of super heavy nuclei with mass numbers >300. Such nuclei would be characterized by large values of v (6-10). using the results from part a) plot the critical radius and mass as a function of v in the range [2-10]."
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "8e1535af-f4a3-471a-86bf-a260747463fc",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 35.263250041651 \\sqrt{\\frac{1}{185.0 \\nu - 211.0}}$"
],
"text/plain": [
"35.263250041651*sqrt(1/(185.0*nu - 211.0))"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nu = sm.symbols('nu')\n",
"critical_radius_eq = sm.Eq(1, nu*sigma_f/(D*B_g**2 + sigma_a))\n",
"R_eqn = sm.solve(critical_radius_eq, R)[1]\n",
"R_eqn"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "b49aa3e1-e284-4b0e-8491-38e0b923094d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7ef4b0bbbef0>"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAHHCAYAAAC2rPKaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABaF0lEQVR4nO3deVxUZdsH8N9h2HcVZRFEEHcUE5dHTEVTedQsW9wfQ7RExdza9C0V31JL39wKNX1IeiwrtbSyJxXXMveU0sQ1NxQXTERBEGbO+wfNxMg2B85wZs78vp+PH50zh3Oue2aEi/u+7vsWRFEUQURERKRCdkoHQERERGQuTHSIiIhItZjoEBERkWox0SEiIiLVYqJDREREqsVEh4iIiFSLiQ4RERGpFhMdIiIiUi0mOkRERKRaTHSoRqWkpEAQBFy8eFHWc80dC1UuMTERgiAoHYaRsmLi+162w4cPIyoqCm5ubhAEAWlpaUqHVIq1vHfW8FraEiY6JMn58+cRHx+P0NBQODs7w9PTE507d8aSJUvw4MGDKl1z3759SExMRHZ2trzBykT/zVUQBOzdu7fU86IoIigoCIIg4Mknn1QgQutS8vUUBAH29vaoX78+Ro4ciatXryodnk0qLCzEwIED8eeff2LRokVYs2YNgoODFYvH0r8nVMTSXksC7JUOgKzH999/j4EDB8LJyQkvvPACwsPD8fDhQ+zduxevvfYafv/9d6xcubLCa4wYMQJDhgyBk5OT4di+ffswe/ZsjBw5Et7e3hWeqyRnZ2esXbsWjz/+uNHxPXv2ICMjw2LitBb/+7//i5CQEOTn5+PAgQNISUnB3r17ceLECTg7O5vtvpb2ubIE58+fx6VLl7Bq1Sq8+OKLSodjNd8TymJpryUx0SETXbhwAUOGDEFwcDB27twJf39/w3MJCQk4d+4cvv/++3K/Pjc3F25ubtBoNNBoNCbdU8q5NaFv375Yv349li5dCnv7v//rrF27FpGRkcjKylIwOuvTp08ftGvXDgDw4osvwsfHB++99x6+/fZbDBo0yGz3tbTPlSW4efMmABglFZbIGt47a3ktbQmHrsgk8+fPx/3795GcnGyU5OiFhYVh0qRJAP6uizh58iSGDRuGWrVqGXpBHh1jT0xMxGuvvQYACAkJMQxnXLx4sdzx+KtXr2L06NEICAiAk5MTQkJCMG7cODx8+BAAcOnSJYwfPx5NmzaFi4sL6tSpg4EDB1Z7XH/o0KG4ffs2UlNTDccePnyIDRs2YNiwYaXONzWOe/fuYfLkyWjYsCGcnJxQr1499OrVC0ePHpV0TllMjUH/np07d87wW7SXlxfi4uKQl5dX6rp79+5F+/bt4ezsjEaNGuGjjz4y4RWsWJcuXQAU/0YsNX4pMZX1uRo5ciQaNmxY6txHa3zM/T6Y+/plGTlyJLp16wYAGDhwIARBQHR0tMmvScljpn5+Kvo/XJXvCceOHUOfPn3g6ekJd3d3PPHEEzhw4EC1YixLZfcp77Usz3vvvQdBEHD27FnDsRMnTkAQBGzYsMGkmKhy7NEhk3z33XcIDQ1FVFSUyV8zcOBANG7cGHPnzoUoimWe8+yzz+LMmTP4/PPPsWjRIvj4+AAA6tatW+b5165dQ4cOHZCdnY0xY8agWbNmuHr1KjZs2IC8vDw4Ojri8OHD2LdvH4YMGYLAwEBcvHgRy5cvR3R0NE6ePAlXV1fpLwCAhg0bolOnTvj888/Rp08fAMAPP/yAu3fvYsiQIVi6dKnR+abGMXbsWGzYsAETJkxAixYtcPv2bezduxfp6elo27atyeeUReprMWjQIISEhGDevHk4evQo/v3vf6NevXp47733DOccP34cvXv3Rt26dZGYmIiioiLMmjULvr6+VXpd9fQ/vGrVqiU5fnPF9Chzvw819T6XFB8fj/r162Pu3LmYOHEi2rdvD19fX3z22WeSXx9TPj+V/R+W+j3h999/R5cuXeDp6YnXX38dDg4O+OijjxAdHY09e/agY8eOkmOs6n3Key3LExAQYHhNGjduDAAIDw9Hy5YtsXnzZjz//POVv+hUOZGoEnfv3hUBiE8//bRJ58+aNUsEIA4dOrTUc6tXrxYBiBcuXDAcW7BgQalj5Z37wgsviHZ2duLhw4dLXVun04miKIp5eXmlntu/f78IQPzPf/5T4fXLoj/v8OHD4ocffih6eHgY7jFw4ECxe/fuoiiKYnBwsNivXz/D15kah5eXl5iQkFBhDKacUxZTY9C/Z6NGjTI695lnnhHr1KljdGzAgAGis7OzeOnSJcOxkydPihqNRjTlW4r+9dy+fbt469Yt8cqVK+KGDRvEunXrik5OTuKVK1ckxy8lprLe99jYWDE4OLjUvfSvi5653wdzX788u3btEgGI69evNxwz9TUpecyUz48p/4elfE8YMGCA6OjoKJ4/f95w7Nq1a6KHh4fYtWvXKsVYFlPvU9ZrWZ7t27eLAMS1a9caHe/Xr58YFRVV6deTaTh0RZXKyckBAHh4eEj6urFjx8oah06nw6ZNm9C/f39DbUdJ+u50FxcXw7HCwkLcvn0bYWFh8Pb2rnQIoDKDBg3CgwcPsHnzZty7dw+bN28uc9hKShze3t44ePAgrl27Vu59TTmnOjHoPfqedenSBbdv3zZ8BrRaLbZu3YoBAwagQYMGhvOaN2+OmJgYSbH17NkTdevWRVBQEJ5//nm4ubnh22+/RWBgoKT45YypMuZ+H2rqfTaXyj4/pv4fNpVWq8W2bdswYMAAhIaGGo77+/tj2LBh2Lt3r+HepsYo131MUbJHp6Q7d+5YdMG1tWGiQ5Xy9PQEUFw/IEVISIiscdy6dQs5OTkIDw+v8LwHDx5g5syZCAoKgpOTE3x8fFC3bl1kZ2fj7t271Yqhbt266NmzJ9auXYuvv/4aWq223O5lU+OYP38+Tpw4gaCgIHTo0AGJiYn4448/jK5lyjlyvBYlEwXg72GkO3fuACh+Dx48eGDoZi+padOmlcZTUlJSElJTU7Fhwwb07dsXWVlZpb65mxK/nDFVxtzvQ029z+ZiyufHlP/Dprp16xby8vLKfJ+bN28OnU6HK1euSIpRrvuYoqxE5+rVqzh48KCh1oeqj4kOVcrT0xMBAQE4ceKEpK8r+VtmTXr55ZcxZ84cDBo0COvWrcO2bduQmpqKOnXqQKfTVfv6w4YNww8//IAVK1agT58+5c6uMDWOQYMG4Y8//sAHH3yAgIAALFiwAC1btsQPP/wg6Rw5XovyZrSI5dRYVUeHDh3Qs2dPPPfcc/j2228RHh6OYcOG4f79+1WOv6rK60nQarVGj839PtTU+2wKU1+Tkmry81NVlhSjl5cX3NzcjNaPmjZtGlxcXDBmzJgaj0etWIxMJnnyySexcuVK7N+/H506dZL12qZ2V9etWxeenp6VJlwbNmxAbGws3n//fcOx/Px82RYfe+aZZxAfH48DBw7gyy+/lCUOf39/jB8/HuPHj8fNmzfRtm1bzJkzx1D0bOo51YnBFHXr1oWLi4vRLBG906dPV+maQPEPn3nz5qF79+748MMPMW3aNACmxS9HTLVq1SrzNbl06VKpY+Z+HyzhfQakvSamMvX/sJTvCa6urmW+z6dOnYKdnR2CgoKqFGtN3cff39/Qo7NmzRp8+umnWLVqVZmzW6lq2KNDJnn99dfh5uaGF198ETdu3Cj1/Pnz57FkyZIqXdvNzQ0AKv2mbGdnhwEDBuC7777DkSNHSj2v/41Mo9GU+u3sgw8+qPA3USnc3d2xfPlyJCYmon///uWeZ0ocWq221NBCvXr1EBAQgIKCApPPqU4MUmg0GsTExGDTpk24fPmy4Xh6ejq2bt1apWvqRUdHo0OHDli8eDHy8/MN96ssfjliatSoEe7evYvffvvNcCwzMxMbN240PDb3+2BJ7zNg2msilan/h039nqDRaNC7d2988803RlPOb9y4YVjcUz/0Xh3mvE9AQACuXbuGn376CfHx8Zg2bRoXGpQZe3TIJI0aNcLatWsxePBgNG/e3Ghl5H379mH9+vUYOXJkla4dGRkJAHjzzTcxZMgQODg4lJtAzJ07F9u2bUO3bt0wZswYNG/eHJmZmVi/fj327t0Lb29vPPnkk1izZg28vLzQokUL7N+/H9u3b0edOnWq2vxSYmNjKz3HlDju3buHwMBAPP/884iIiIC7uzu2b9+Ow4cPG347N+Wc6sQg1ezZs7FlyxZ06dIF48ePR1FRET744AO0bNnS6IdiVbz22msYOHAgUlJSMHbsWJPjr25MQ4YMwRtvvIFnnnkGEydORF5eHpYvX44mTZoYinnN/T5Y2vtsymtSFab8H5byPeGdd95BamoqHn/8cYwfPx729vb46KOPUFBQgPnz51c5zpq6T0BAAA4cOIAnn3wSI0aMwLx582SLmf6i0GwvslJnzpwRX3rpJbFhw4aio6Oj6OHhIXbu3Fn84IMPxPz8fFEU/57GeevWrVJfX96U7rffflusX7++aGdnZ3i+vHMvXbokvvDCC4bpyKGhoWJCQoJYUFAgiqIo3rlzR4yLixN9fHxEd3d3MSYmRjx16pQYHBwsxsbGVhpLeTGXNR22pEenl5sSR0FBgfjaa6+JERERooeHh+jm5iZGRESIy5YtM1zHlHPKY+prUd57Vt5rtGfPHjEyMlJ0dHQUQ0NDxRUrVpQ57bgsFb2eWq1WbNSokdioUSOxqKjI5PilxFRem7Zt2yaGh4eLjo6OYtOmTcVPP/3U6OvN/T7UxPtcnvKmRFf2muhJ/fxU9n9YFKV9Tzh69KgYExMjuru7i66urmL37t3Fffv2VSvGsphyHynTy0VRFKdOnSoCEMeOHWuYXk/yEkTRgqrEiIiIiGTEGh0iIiJSLSY6REREpFpMdIiIiEi1mOgQERGRajHRISIiItViokNERESqZfMLBup0Oly7dg0eHh6Sd84lIiIiZYiiiHv37iEgIAB2duX329h8onPt2jVZ9kIhIiKimnflyhUEBgaW+7zNJzoeHh4Ail8oOfZE0SssLMS2bdvQu3dvODg4yHZdS6H29gHqb6Pa2weov41sn/VTexvN2b6cnBwEBQUZfo6Xx+YTHf1wlaenp+yJjqurKzw9PVX74VVz+wD1t1Ht7QPU30a2z/qpvY010b7Kyk5YjExERESqxUSHiIiIVMvqE50rV64gOjoaLVq0QOvWrbF+/XqlQyIiIiILYfU1Ovb29li8eDHatGmD69evIzIyEn379oWbm5vSoRER2SStVovCwsJqX6ewsBD29vbIz8+HVquVITLLo/Y2Vqd9Dg4O0Gg01Y7B6hMdf39/+Pv7AwD8/Pzg4+ODP//8k4kOEVENE0UR169fR3Z2tmzX8/Pzw5UrV1S7zpna21jd9nl7e8PPz69ar43iic6PP/6IBQsW4JdffkFmZiY2btyIAQMGGJ2TlJSEBQsW4Pr164iIiMAHH3yADh06lLrWL7/8Aq1Wy3VxiIgUoE9y6tWrB1dX12r/4NbpdLh//z7c3d0rXBDOmqm9jVVtnyiKyMvLw82bNwHA0KFRFYonOrm5uYiIiMCoUaPw7LPPlnr+yy+/xNSpU7FixQp07NgRixcvRkxMDE6fPo169eoZzvvzzz/xwgsvYNWqVTUZPhERoXi4Sp/k1KlTR5Zr6nQ6PHz4EM7OzqpMAgD1t7E67XNxcQEA3Lx5E/Xq1avyMJbiiU6fPn3Qp0+fcp9fuHAhXnrpJcTFxQEAVqxYge+//x4ff/wxpk2bBgAoKCjAgAEDMG3aNERFRVV4v4KCAhQUFBge5+TkACgeR5RjTFlPfy05r2lJ1N4+QP1tVHv7APW30ZLaV1BQAFEU4ezsDJ1OJ8s1RVE0/C3XNS2N2ttY3fY5OztDFEU8ePAATk5ORs+Z+rkXRH0UFkAQBKOhq4cPH8LV1RUbNmwwGs6KjY1FdnY2vvnmG4iiiGHDhqFp06ZITEys9B6JiYmYPXt2qeNr166Fq6urTC0hIrIt9vb28PPzQ1BQEBwdHZUOh1Ti4cOHuHLlCq5fv46ioiKj5/Ly8jBs2DDcvXu3wgV/Fe/RqUhWVha0Wi18fX2Njvv6+uLUqVMAgJ9//hlffvklWrdujU2bNgEA1qxZg1atWpV5zenTp2Pq1KmGx/olpHv37i37ysipqano1auXale7VHP7APW3Ue3tA9TfRktqX35+Pq5cuQJ3d3c4OzvLck39po1q3nRZ7W2sbvvy8/Ph4uKCrl27lvpc6UdkKmPRiY4pHn/8cUndYU5OTqW6v4DiaWzm+EZhrutaCrW3D1B/G9XePkD9bbSE9mm1WgiCADs7O9lqTfTf2/XXtTS7d+9G9+7dcefOHXh7e5d7XsOGDTF58mRMnjy51HNVaWN0dDTatGmDxYsXVyHq6klMTMSmTZuQlpYGABg5ciSys7MNHQ2Pqu57aGdnB0EQyvyMm/qZt7xPTgk+Pj7QaDS4ceOG0fEbN27Az89PoahMlJEBn+PHgYwMpSMhIqIKXL9+HS+//DJCQ0Ph5OSEoKAg9O/fHzt27Kjw66KiopCZmQkvLy8AQEpKSpkJz+HDhzFmzBhzhF6mlJQUCIJgSC78/f0xePBgXL58WfZ7LVmyBCkpKbJfV04Wneg4OjoiMjLS6MOm0+mwY8cOdOrUScHIKpGcDPuwMHSeMQP2YWFAcrLSERERURkuXryIyMhI7Ny5EwsWLMDx48exZcsWdO/eHQkJCeV+XWFhIRwdHU1a46Vu3bo1XgPq6emJzMxMXL16FV999RVOnz6NgQMHyn4fLy+vCnuzLIHiic79+/eRlpZm6Aa7cOEC0tLSDJnn1KlTsWrVKnzyySdIT0/HuHHjkJuba5iFZXEyMoAxYyDou+t0OiA+nj07REQmysgAdu2qmW+b48ePhyAIOHToEJ577jk0adIELVu2xNSpU3HgwAHDeYIgYPny5Xjqqafg5uaGOXPmYPfu3RAEAdnZ2di9ezfi4uJw9+5dQ2+KfoJMw4YNjYaZsrOzER8fD19fX7i6uqJTp07YvHkzAOD27dsYOnQo6tevD1dXV7Rq1Qqff/655HYJggA/Pz/4+/sjKioKo0ePxqFDh4zqWt544w00adIErq6uCA0NxYwZM0rNZHr33Xfh6+sLDw8PjB49Gvn5+UbPjxw50miy0KNtBYAuXboYJgGJoojExEQ0aNAATk5OCAgIwMSJEyW3TwrFa3SOHDmC7t27Gx7rC4VjY2ORkpKCwYMH49atW5g5cyauX7+ONm3aYMuWLaUKlC3G2bPAozVDWi1w7hwQGKhMTEREViI5GRgzpvjbqJ2dgMWLHVFBx0q1/Pnnn9iyZQvmzJlT5mr6j/ZUJCYm4t1338XixYthb2+PP/74w/BcVFQUFi9ejJkzZ+L06dMAAHd391LX1Ol06NOnD+7du4dPP/0UISEhOHLkiGGNmPz8fERGRuKNN96Ap6cnvv/+e4wYMQKNGjUqc6FcU9y8eRMbN26ERqMxWovGw8MDKSkpCAgIwPHjx/HSSy/Bw8MDr7/+OgBg3bp1SExMRFJSEh5//HGsWbMGS5cuRWhoaJXiAICvvvoKixYtwhdffIGWLVvi+vXr+PXXX6t8PVMonuhER0ejshnuEyZMwIQJE2S9b1JSEpKSkuTfW6RxY8DOzjjZ0WiAsDB570NEpDJ/dYgbvn3qdAKmTHHB00+LaNBA/vudO3cOoiiiWbNmJp0/bNgwo9GEkomOo6MjvLy8DD0p5dm+fTsOHTqE9PR0NGnSBDqdDj4+PoZZv/Xr18err75qOP/ll1/G1q1bsW7dOkmJzt27d+Hu7m5YYRgAJk6caJTQvfXWW4Z/N2zYEK+++iq++OILQ6KzePFijB49GqNHjwYAvPPOO9i+fXupXh0pLl++DD8/P/Ts2RMODg5o0KBBlRM4Uyk+dKWUhIQEnDx5EocPH5b3woGBwMqVEP/KmkWNBvjoI/bmEBFVouwOcQHnzpnnflKXkWvXrl2175mWlobAwEA0adKkzOe1Wi3efvtttGrVCrVr14a7uzu2bt0quZDYw8MDaWlpOHLkCN5//320bdsWc+bMMTrnyy+/ROfOneHn5wd3d3e89dZbRvdJT09Hx44djb6muvWxAwcOxIMHDxAaGoqXXnoJGzduLLU+jtxsNtExq9GjUXT2LPa+/TaKzp4F/sqGiYiofPoO8ZI0GtFsHeKNGzeGIAiGddkqI8dm0fptDcqzYMECLFmyBG+88QZ27dqFtLQ0xMTE4OHDh5LuY2dnh7CwMDRv3hxTp07FP/7xD4wbN87w/P79+zF8+HD07dsXmzdvxrFjx/Dmm29Kvk9Z9300gSyZyAQFBeH06dNYtmwZXFxcMH78eHTt2tWsq3sz0TGXwEDcbtWKPTlERCb6q0Mc+jISjUbEokUPzPZttHbt2oiJiUFSUhJyc3NLPS91F3ZHR8dKyyFat26NjIwMnDlzpsznf/75Zzz99NP417/+hYiICISGhpZ7rhTTpk3Dl19+iaNHjwIA9u3bh+DgYLz55pto164dGjdujEuXLhl9TfPmzXHw4EGjYyULtMtSt25dZGZmGh7n5OSUuq6Liwv69++PpUuXYvfu3di/fz+OHz9eneZViIkOERFZjNGjgYsXi2dd/fGHiBEjqtfDUBl9rWaHDh3w1Vdf4ezZs0hPT8fSpUslD9M0bNgQ9+/fx44dO5CVlWWojSmpW7du6Nq1K5577jmkpqbiwoULSE1NxZYtWwAU9zKlpqZi3759SE9PR3x8fKm15KoiKCgIzzzzDGbOnGm4z+XLl/HFF1/g/PnzWLp0KTZu3Gj0NZMmTcLHH3+M1atX48yZM5g1axZ+//33Cu/To0cPrFmzBj/99BOOHz+OkSNHGhVAp6SkIDk5GSdOnMAff/yBTz/9FC4uLggODq52G8vDRIeIiCxKYCAQHV0zHeKhoaE4evQounfvjldeeQXh4eHo1asXduzYgeXLl0u6VlRUFMaOHYvBgwejbt26mD9/fpnnffXVV2jfvj2GDh2K8PBwzJo1y9AT9NZbb6Ft27aIiYlBdHQ0/Pz8jKZvV8eUKVPw/fff49ChQ3jqqacwZcoUTJgwAW3atMG+ffswY8YMo/MHDx6MGTNm4PXXX0dkZCQuXbpkNPxVlunTp6Nbt2548skn0a9fPzz99NNo2LCh4Xlvb2+sWrUKnTt3RuvWrbF9+3Z89913su14XxaL2tRTCTk5OfDy8qp0UzCpCgsL8d///hd9+/ZVfGl2c1B7+wD1t1Ht7QPU30ZLal9+fj4uXLiAkJAQ2fa60ul0yMnJgaenp0VuASEHtbexuu2r6HNl6s9v9b2qRERERH9hokNERESqZbOJTlJSElq0aIH27dsrHQoRERGZic0mOmZbMJCIiIgshs0mOkRERKR+THSIiIhItZjoEBERkWox0SEiIiLVYqKjpIyM4nXOMzKUjoSIiEiVmOgoJTkZCA4GevQo/js5WemIiIhIgt27d0MQhEo3/2zYsCEWL14s232jo6MxefJk2a6ndkx0lJCRAYwZA+h0xY91OiA+nj07REQKuH79Ol5++WWEhobCyckJQUFB6N+/P3bs2FHh10VFRSEzMxNeXl4Aijes9Pb2LnXe4cOHMWbMGHOEXqaUlBQIgoDmzZuXem79+vUQBMFo/ym1s9lER9EFA8+e/TvJ0dNqgXPnaj4WIiIbdvHiRURGRmLnzp1YsGABjh8/ji1btqB79+5ISEgo9+sKCwvh6OgIPz8/CIJQ4T3q1q0LV1dXuUOvkJubG27evIn9+/cbHU9OTkaDBg1qNBal2Wyio+iCgY0bA49ubqbRAGFhNR8LEZGlqcH6xfHjx0MQBBw6dAjPPfccmjRpgpYtW2Lq1Kk4cOCA4TxBELB8+XI89dRTcHNzw5w5c4yGrnbv3o24uDjcvXsXgiBAEAQkJiYCKD10lZ2djfj4ePj6+sLV1RWdOnXC5s2bAQC3b9/G0KFDUb9+fbi6uqJVq1b4/PPPJbfL3t4ew4YNw8cff2w4lpGRgd27d2PYsGFG554/fx5PP/00fH194e7ujvbt22P79u1G5yxbtgyNGzeGs7MzfH198fzzzxue27BhA1q1agUXFxfUqVMHPXv2RG5uruSYzcVmEx1FBQYCK1cWJzdA8d8ffVR8nIjIlpWoXxRCQuC4Zo3ZbvXnn39iy5YtSEhIgJubW6nnHx2GSkxMxDPPPIPjx49j1KhRRs9FRUVh8eLF8PT0RGZmJjIzM/Hqq6+WuqZOp0OfPn3w888/49NPP8WJEycwa9YsaP76eZCfn4/IyEh8//33OHHiBMaMGYMRI0bg0KFDkts3atQorFu3Dnl5eQCKh7T++c9/wtfX1+i8+/fvo2/fvtixYweOHTuGf/7zn+jfvz8uX74MADhy5AgmTpyI//3f/8Xp06exZcsWdO3aFQCQmZmJoUOHYtSoUUhPT8fu3bvx7LPPQhRFyfGai73SAdis0aOBmJji4aqwMCY5RESP1C8KOh1cpkyB+PTTgBmGW86dOwdRFNGsWTOTzh82bBji4uIMj//44w/Dvx0dHeHl5QVBEODn51fuNbZv345Dhw4hPT0dTZo0gU6ng4+PDzw9PQEA9evXN0qQXn75ZWzduhXr1q1Dhw4dJLXvscceQ2hoKDZs2IARI0YgJSUFCxcuNIobACIiIhAREWF4/Pbbb2Pjxo349ttvMWHCBFy+fBlubm548skn4eHhgeDgYDz22GMAihOdoqIiPPvsswgODgYAtGrVSlKc5sYeHSUFBgLR0UxyiIiAMusXBTPWL0rtdWjXrl2175mWlobAwEA0adKkzOe1Wi3efvtttGrVCrVr14a7uzu2bt1q6F2RatSoUVi9ejX27NmD3Nxc9O3bt9Q59+/fx6uvvormzZvD29sb7u7uSE9PN9yzV69eCA4ORmhoKEaMGIHPPvvM0EsUERGBJ554Aq1atcLAgQOxatUq3Llzp0qxmgsTHSIisgxl1C+KZqxfbNy4MQRBwKlTp0w6v6zhLalcXFwqfH7BggVYsmQJ3njjDezatQtpaWmIiYnBw4cPq3S/4cOH48CBA0hMTMSIESNgb196IOfVV1/Fxo0bMXfuXPz0009IS0tDq1atDPf08PDA0aNH8fnnn8Pf3x8zZ85EREQEsrOzodFokJqaih9++AEtWrTABx98gKZNm+LChQtVitccmOgQEZFleKR+UdRo8GDRIrP1eteuXRsxMTFISkoqs3i2svVxHuXo6AitVlvhOa1bt0ZGRgbOnDlT5vM///wznn76afzrX/9CREQEQkNDyz3XFLVr18ZTTz2FPXv2lKorKnnPkSNH4plnnkGrVq3g5+eHixcvGp1jb2+Pnj17Yv78+fjtt99w8eJF7Ny5E0BxoXbnzp0xe/ZsHDt2DI6Ojti4cWOVY5YbEx0iIrIco0cDFy8Cu3ZB/OMPPBwxwqy3S0pKglarRYcOHfDVV1/h7NmzSE9Px9KlS9GpUydJ12rYsCHu37+PHTt2ICsryzC8U1K3bt3QtWtXPPfcc0hNTcWFCxeQmpqKLVu2ACjuZUpNTcW+ffuQnp6O+Ph43Lhxo1ptTElJQVZWVrm1SI0bN8bXX3+NtLQ0/Prrrxg2bBh0JYYQN2/ejKVLlyItLQ2XLl3Cf/7zH+h0OjRt2hQHDx7E3LlzceTIEVy+fBlff/01bt26VeYaPkphokNERJalBusXQ0NDcfToUXTv3h2vvPIKwsPD0atXL+zYsQPLly+XdK2oqCiMHTsWgwcPRt26dTF//vwyz/vqq6/Qvn17DB06FOHh4Zg1a5ahJ+itt95C27ZtERMTg+joaPj5+WHAgAHVaqN+2nd5Fi5ciFq1aiEqKgr9+/dHTEwM2rZta3je29sbX3/9NXr06IHmzZtjxYoV+Pzzz9GyZUt4enrixx9/RN++fdGkSRO89dZbeP/999GnT59qxSwnQbSkOWAKyMnJgZeXF+7evWuoepdDYWEh/vvf/6Jv375wcHCQ7bqWQu3tA9TfRrW3D1B/Gy2pffn5+bhw4QJCQkLg7OwsyzV1Oh1ycnLg6ekJu0fXHlMJtbexuu2r6HNl6s9v9b2qRERERH+x2URH0S0giIiIqEbYbKKj6BYQREREVCNsNtGxKjW47wsREZGaMNGxdCX2fUFwcPFjIiILZePzW0hmcnyemOhYskf2fYFOB8THs2eHiCyOftZXWWvHEFWV/vNUnVmF3NTTkpWx7wv0+75wfywisiAajQbe3t64efMmAMDV1RWCIFTrmjqdDg8fPkR+fr4qp14D6m9jVdsniiLy8vJw8+ZNeHt7G3Z3rwomOpZMv+9LyWTHjPu+EBFVh37Xbn2yU12iKOLBgwdwcXGpdtJkqdTexuq2z9vbu8Ld4E3BRMeS6fd9iY8v7snRaICPPmJvDhFZJEEQ4O/vj3r16qGwsLDa1yssLMSPP/6Irl27Kr4gormovY3VaZ+Dg0O1enL0mOhYutGjgZiY4uGqsDAmOURk8TQajSw/oDQaDYqKiuDs7KzKJABQfxstoX1MdKxBYCATHCIioipQX+UTERER0V+Y6BAREZFqMdEhIiIi1WKiQ0RERKpls4kOdy8nIiJSP5tNdLh7ORERkfrZbKKjOtzhnIiIqBQmOmrAHc6JiIjKxETH2nGHcyIionIx0bF2Fe1wTkREZOOY6Fg7/Q7nJXGHcyIiIgBMdKyffodz/QZ63OGciIjIgJt6qgF3OCciIioTEx214A7nREREpXDoioiIiFSLiQ4RERGpFhMdIiIiUi0mOkRERKRaTHRsCffDIiIiG8NEx1ZwPywiIrJBNpvoJCUloUWLFmjfvr3SoZgf98MiIiIbZbOJTkJCAk6ePInDhw8rHYr5cT8sIiKyUTab6NgU7odFREQ2iomOLeB+WEREZKO4BYSt4H5YRERkg5jo2BLuh0VERDaGQ1dERESkWkx0iIiISLWY6BAREZFqMdEhY9wmgoiIVISJDv2N20QQEZHKMNGhYtwmgoiIVIiJDhXjNhFERKRCTHSoGLeJICIiFWKiQ8W4TQQREakQV0amv3GbCCIiUhkmOmSM20QQEZGKcOiKiIiIVIuJDhEREakWh65IuowMCOnpcM7KUjoSIiKiCtlsopOUlISkpCRotVqlQ7EuycnAmDGw1+nQWxCKX78xY5SOioiIqEwmJTo5OTmSL+zp6Sn5a2pSQkICEhISkJOTAy8vL6XDsQ6PrJ4siCI048cDffuygJmIiCySSYmOt7c3BEEw+aKCIODMmTMIDQ2tcmBkgcpYPVnQr57MRIeIiCyQyUNXGzZsQO3atSs9TxRF9O3bt1pBkYXSr55cItkRNRoIXD2ZiIgslEmJTnBwMLp27Yo6deqYdNHQ0FA4ODhUKzCyQPrVk+PjAa0WOjs76JYtgz17c4iIyEKZlOhcuHBB0kVPnDhRpWDICvy1enLRqVPYcekSerzwgtIRERERlYvr6JB0gYEQu3VDvo+P0pEQERFVSPL08qVLl5Z5XBAEODs7IywsDF27doVGvzkkERERkUIkJzqLFi3CrVu3kJeXh1q1agEA7ty5A1dXV7i7u+PmzZsIDQ3Frl27EBQUJHvAZEUyMopnajVuzFlZRESkCMlDV3PnzkX79u1x9uxZ3L59G7dv38aZM2fQsWNHLFmyBJcvX4afnx+mTJlijnjJWiQnA8HBQI8exX8nJysdERER2SDJic5bb72FRYsWoVGjRoZjYWFh+L//+z9Mnz4dgYGBmD9/Pn7++WdZAyUr8sjCgtDpimdqZWQoGxcREdkcyYlOZmYmioqKSh0vKirC9evXAQABAQG4d+9e9aMj61TGwoLQLyxIRERUgyQnOt27d0d8fDyOHTtmOHbs2DGMGzcOPXr0AAAcP34cISEh8kVJ1kW/sGBJGg3AhQWJiKiGSU50kpOTUbt2bURGRsLJyQlOTk5o164dateujeS/6jDc3d3x/vvvyx4sWQn9woL6mXcaDfDRRyxIJiKiGid51pWfnx9SU1Nx+vRpnD59GgDQtGlTNG3a1HBO9+7d5YuQrNNfCwvi3LninhwmOUREpADJiY7eo8kNUSmBgUxwiIhIUSYNXU2dOhW5ubkmX3T69On4888/qxwUERERkRxMSnSWLFmCvLw8ky+alJSE7OzsqsZEtiQjA9i1i1PPiYjILEwauhJFEU2aNIEgCCZdVErvD9mw5OS/19uxsysuYB49WumoiIhIRUxKdFavXi35wr6+vpK/hmxIeYsKxsSwroeIiGRjUqITGxtr7jjI1lS0qCATHSIikonkdXSIZMFFBYmIqAYw0SFlcFFBIiKqAVVeR4eo2rioIBERmRkTHVIWFxUkIiIzqvbQVU5ODjZt2oT09HQ54iEqjWvtEBFRFUlOdAYNGoQPP/wQAPDgwQO0a9cOgwYNQuvWrfHVV1/JHiDZuORkIDgY6NGj+O+/No4lIiIyheRE58cff0SXLl0AABs3boQoisjOzsbSpUvxzjvvyB4g2bDy1tphzw4REZlIcqJz9+5d1K5dGwCwZcsWPPfcc3B1dUW/fv1w9uxZ2QMkG1bRWjtEREQmkJzoBAUFYf/+/cjNzcWWLVvQu3dvAMCdO3fg7Owse4Bkw7jWDhERVZPkRGfy5MkYPnw4AgMDERAQgOjoaADFQ1qtWrWSOz6zSUpKQosWLdC+fXulQ6HycK0dIiKqJsnTy8ePH48OHTrgypUr6NWrF+z++o07NDTUqmp0EhISkJCQgJycHHh5eSkdDpWHa+0QEVE1VGkdnXbt2qFdu3YAAK1Wi+PHjyMqKgq1atWSNTgiAKattZORUVzT07gxkyEiIjKo0tBV8l9TfLVaLbp164a2bdsiKCgIu3fvljs+ospxCjoREZVDcqKzYcMGREREAAC+++47XLhwAadOncKUKVPw5ptvyh4gUYU4BZ2IiCogOdHJysqCn58fAOC///0vBg4ciCZNmmDUqFE4fvy47AESVYhT0ImIqAKSEx1fX1+cPHkSWq0WW7ZsQa9evQAAeXl50OhnxxDVFE5BJyKiCkhOdOLi4jBo0CCEh4dDEAT07NkTAHDw4EE0a9ZM9gCJKsQp6EREVAHJs64SExMRHh6OK1euYODAgXBycgIAaDQaTJs2TfYAiSrFKehERFSOKk0vf/7550sdi42NrXYwRFXGKehERFSGKiU6ubm52LNnDy5fvoyHDx8aPTdx4kRZAiOSVXLy37Oz7OyKh7tGj1Y6KiIiMjPJic6xY8fQt29f5OXlITc3F7Vr10ZWVhZcXV1Rr149Jjpkecqbgh4Tw54dIiKVk1yMPGXKFPTv3x937tyBi4sLDhw4gEuXLiEyMhL/93//Z44YiaqHU9CJiGyW5EQnLS0Nr7zyCuzs7KDRaFBQUICgoCDMnz8f//M//2OOGImqh1PQiYhsluREx8HBwbCRZ7169XD58mUAgJeXF65cuSJvdERy4BR0IiKbJblG57HHHsPhw4fRuHFjdOvWDTNnzkRWVhbWrFmD8PBwc8RIVH1SpqDrZ2c1bFhj4RERkXlI7tGZO3cu/P39AQBz5sxBrVq1MG7cONy6dQsrV66UPUAi2QQGAtHRFSc5JTYItQ8LQ4PU1BoLj4iI5Ce5R6ddu3aGf9erVw9btmyRNSAixTwyO0vQ6RCxfDm0r7wChIQoHBwREVWF5B4dItUqY3aWnU4H4fx5hQIiIqLqMrlHp0ePHiadt3PnzioHQ6Qo/eysEsmOzs4OYqNGCgZFRETVYXKis3v3bgQHB6Nfv35wcHAwZ0xEytDPzoqPB7RaiBoNfh07FuGcnUVEZLVMTnTee+89rF69GuvXr8fw4cMxatQozrIi9SkxO6soOBiXf/sN/JQTEVkvk2t0XnvtNZw8eRKbNm3CvXv30LlzZ3To0AErVqxATk6OOWMkqlmmzM7KyAB27Sr+m4iILJbkYuROnTph1apVyMzMREJCAj7++GMEBAQw2SHbUWIKOoKDix8TEZFFqvKsq6NHj2LPnj1IT09HeHg463bINpS3QSh7doiILJKkROfatWuYO3cumjRpgueffx61a9fGwYMHceDAAbi4uJgrRiLLwQ1CiYisisnFyH379sWuXbvQu3dvLFiwAP369YO9veT1BomsWxlT0LlBKBGR5TI5U9myZQv8/f1x+fJlzJ49G7Nnzy7zvKNHj8oWHJHFeWQKeoUbhOr3zGrcmBuIEhEpxOREZ9asWeaMg8h6mLJBaHLy37U8dnbFydHo0TUfKxGRjWOiQ1QVgYHl99KUV7AcE8OeHSKiGsa9rojkxoJlIiKLwUSHSG76guWSWLBMRKQIJjpEctMXLGs0xY8rKlgmIiKz4vxwInMwpWAZ4MwsIiIzY48OkblUtmcWt5IgIjI7k3p0li5davIFJ06cWOVgiGwGZ2YREdUIkxKdRYsWGT2+desW8vLy4O3tDQDIzs6Gq6sr6tWrx0SHyBQVzcxiokNEJBuThq4uXLhg+DNnzhy0adMG6enp+PPPP/Hnn38iPT0dbdu2xdtvv23ueInUgTOziIhqhOQanRkzZuCDDz5A06ZNDceaNm2KRYsW4a233pI1OCLVkjIzKyMD2LWLO6QTEVWB5EQnMzMTRUVFpY5rtVrcuHFDlqCIbMLo0cDFi8VJzMWLZW8RwYJlIqJqkZzoPPHEE4iPjzfavPOXX37BuHHj0LNnT1mDI1K9imZmlVewzJ4dIiKTSU50Pv74Y/j5+aFdu3ZwcnKCk5MTOnToAF9fX/z73/82R4xEtolbSRARVZvkBQPr1q2L//73vzhz5gxOnToFAGjWrBmaNGkie3BENk1fsFwy2WHBMhGRJFVeGblJkyZMbojMSV+wHB9f3JNT2VYSXGWZiKgUyYmOVqtFSkoKduzYgZs3b0L3SNf6zp07ZQvOVM888wx2796NJ554Ahs2bKjx+xOZjalbSSQn/13PY2dXnCCVVdxMRGRjJCc6kyZNQkpKCvr164fw8HAIgmCOuCTHNGrUKHzyySdKh0Ikv8DAintouMoyEVG5JCc6X3zxBdatW4e+ffuaI54qiY6Oxu7du5UOg0gZXGWZiKhckmddOTo6IkzGYsgff/wR/fv3R0BAAARBwKZNm0qdk5SUhIYNG8LZ2RkdO3bEoUOHZLs/kdXjKstEROWSnOi88sorWLJkCURRlCWA3NxcREREICkpqcznv/zyS0ydOhWzZs3C0aNHERERgZiYGNy8eVOW+xNZPa6yTERULslDV3v37sWuXbvwww8/oGXLlnBwcDB6/uuvv5Z0vT59+qBPnz7lPr9w4UK89NJLiIuLAwCsWLEC33//PT7++GNMmzZNavgoKChAQUGB4XFOTg4AoLCwEIWFhZKvVx79teS8piVRe/sAK2vjCy8APXpAOH8eYqNGxUnOI3ELq1dDM24cBJ0Oop0ddB9+CAQEWEf7qsiq3sMqYPusn9rbaM72mXpNyYmOt7c3nnnmGckBVcXDhw/xyy+/YPr06YZjdnZ26NmzJ/bv31+la86bNw+zZ88udXzbtm1wdXWtcqzlSU1Nlf2alkTt7QOssI2//Vb8pwTnrCz0HjsWwl89sYJOB4cJE+C8cqX1ta8K1N5Gts/6qb2N5mhfXl6eSedJTnRWr14tOZiqysrKglarha+vr9FxX19fw2KFANCzZ0/8+uuvyM3NRWBgINavX49OnTqVec3p06dj6tSphsc5OTkICgpC79694enpKVvshYWFSE1NRa9evUr1eqmB2tsHqKuNwu7dhiRHz06ng1tmJroMHWr17SuPmt7DsrB91k/tbTRn+/QjMpWp8oKBlmT79u0mn6vftuJRDg4OZvmQmeu6lkLt7QNU0sbmzUutsixqNMj19y/dPhUuPKiK97ACbJ/1U3sbzdE+U69XpURnw4YNWLduHS5fvoyHDx8aPVdys8/q8vHxgUajKbUr+o0bN+Dn5yfbfYhUr4xVlrXLliHfx8f4PC48SEQqI3nW1dKlSxEXFwdfX18cO3YMHTp0QJ06dfDHH39UWFRcFY6OjoiMjMSOHTsMx3Q6HXbs2FHu0BQRlWP0aODixeJZVxcvQvyrwN+Au6UTkQpJ7tFZtmwZVq5ciaFDhyIlJQWvv/46QkNDMXPmTPz555+SA7h//z7OldiN+cKFC0hLS0Pt2rXRoEEDTJ06FbGxsWjXrh06dOiAxYsXIzc31zALi4gkKLnK8qMzFrjwIBGpkORE5/Lly4iKigIAuLi44N69ewCAESNG4B//+Ac+/PBDSdc7cuQIunfvbnisLxSOjY1FSkoKBg8ejFu3bmHmzJm4fv062rRpgy1btpQqUCaiauJu6USkQpKHrvz8/Aw9Nw0aNMCBAwcAFPfEVGURwejoaIiiWOpPSkqK4ZwJEybg0qVLKCgowMGDB9GxY0fJ93lUUlISWrRogfbt21f7WkSqwIUHiUiFJCc6PXr0wLfffgsAiIuLw5QpU9CrVy8MHjy4xtbXkUNCQgJOnjyJw4cPKx0KkeV4pI6nzELk5GQgOBjo0aP47+Tkmo6SiMhkkoeuVq5cCd1fXdsJCQmoU6cO9u3bh6eeegrx8fGyB0hENayi3dK5UzoRWRnJiY6dnR3sSmwgOGTIEAwZMkTWoIjIQrFgmYisjOShKyKyYVJ2SmcdDxFZACY6RGQ6UwuWWcdDRBaCiQ4RSVNZwTIXHiQiC6KKva6IqIZVVLDMOh4isiDs0SEieUmp4wFYy0NEZmVSj07t2rUlXVQQBBw9ehTBwcFVCqomJCUlISkpCVqtVulQiNSljA1Ey114kJuIEpGZmZToZGdnY/HixfDy8qr0XFEUMX78eItPIBISEpCQkICcnByT2kVEEoweXby2zrlzxT055a2uzDV5iMjMTK7RGTJkCOrVq2fSuS+//HKVAyIilaiojgdgLQ8R1QiTEh3do9+MKqHf6JOIqFxSNhHNyChOjBo3ZhJERJKwGJmIlME1eYioBlRpevm1a9ewd+9e3Lx5s1Rvz8SJE2UJjIhsQGW1PKzjIaJqkpzopKSkID4+Ho6OjqhTpw4EQTA8JwgCEx0ikoZr8hCRGUlOdGbMmIGZM2di+vTpRpt7EhHJjnU8RFRNkjOVvLw8DBkyhEkOEZkf63iIqJokZyujR4/G+vXrzRELEVFp3FuLiKpB8tDVvHnz8OSTT2LLli1o1aoVHBwcjJ5fuHChbMGZE1dGJrIiVa3j8fU1Ps7hLSKbU6VEZ+vWrWjatCkAlCpGthZcGZlIJUyt4+F2E0Q2SXKi8/777+Pjjz/GyJEjzRAOEZFEFe2tVVhYfA6nqRPZLMmJjpOTEzp37myOWIiIqqaS9XiEc+c4TZ3IRkkuRp40aRI++OADc8RCRFR1gYFAdHSZiYsYFlY8XFVSRdPUd+1iMTORSkju0Tl06BB27tyJzZs3o2XLlqWKkb/++mvZgiMikkVFw1slsY6HSHUkJzre3t549tlnzRELEZH5cLsJIpskOdFZvXq1OeIgIjI/Obeb4FR1IqsguUbnwoULOHv2bKnjZ8+excWLF+WIiYio5umnqZdUXh0PV2ImshqSE52RI0di3759pY4fPHiQU86JyHqZut0EV2ImsiqSE51jx46VOb38H//4B9LS0uSIiYhIGZVtNwFUPMT1KM7gIlKc5ERHEATcu3ev1PG7d+9a1XYKSUlJaNGiBdq3b690KERkSSqYpg7A9CEuDm8RWQTJiU7Xrl0xb948o6RGq9Vi3rx5ePzxx2UNzpwSEhJw8uRJHD58WOlQiMiamDLExeEtIoshedbVe++9h65du6Jp06bo0qULAOCnn35CTk4Odu7cKXuAREQWp7Kp6lJncBGR2Uju0WnRogV+++03DBo0CDdv3sS9e/fwwgsv4NSpUwgPDzdHjERElqeiIS4pM7hYx0NkVpJ7dAAgICAAc+fOlTsWIiJ14ErMRBbDpB6d3377DbpHu2Er8Pvvv6OoqKjKQRERWb3KZnCxjoeoRpiU6Dz22GO4ffu2yRft1KkTLl++XOWgiIhUoaLhLU5TJ6oRJg1diaKIGTNmwNXV1aSLPnz4sFpBERGpnr6Op2SyU940dQ5vEVWZSYlO165dcfr0aZMv2qlTJ7i4uFQ5KCIi1TOljocbjRJVm0mJzu7du80cBhGRDZJzmjo3GSUqk+Tp5UREJCMZpqkLq1dzFWaicjDRISKyVCaswuyclQXNuHGcvUVUjiqto0NERDWkkuEt98xMCFJWYeYQF9kYm+3R4aaeRGQ1Khjeuu/vD9HUVZi50SjZIJtNdLipJxGpQb6PD7TLl1e8ySjABQrJZklOdD755BN8//33hsevv/46vL29ERUVhUuXLskaHBERVU6Mi6t4FWaACxSSzZKc6MydO9ewRs7+/fuRlJSE+fPnw8fHB1OmTJE9QCIiMkFFs7cA0zca5fAWqYzkROfKlSsI++s/xqZNm/Dcc89hzJgxmDdvHn766SfZAyQiIhmYMINL0vAWe33ISkhOdNzd3Q37Xm3btg29evUCADg7O+PBgwfyRkdERPKpbKNRU4e32OtDVkTy9PJevXrhxRdfxGOPPYYzZ86gb9++AIp3LG/YsKHc8RERkZwCAysf3qpo/y1uS0FWRnKPTlJSEjp16oRbt27hq6++Qp06dQAAv/zyC4YOHSp7gEREVENMGd5iUTNZGck9Ot7e3vjwww9LHZ89e7YsARERkYIq23+Lu66TlZGc6Pz4448VPt+1a9cqB0NERBagouEt7rpOVkZyohMdHV3qmCAIhn9rtdpqBURERBaOu66TFZFco3Pnzh2jPzdv3sSWLVvQvn17bNu2zRwxEhGRpZFh13XO3qKaILlHx8vLq9SxXr16wdHREVOnTsUvv/wiS2BERGSlzDG8xZ4fqiLZ9rry9fXF6dOn5bocERFZM7nW7AEgrF7Nnh+qMsk9Or/99pvRY1EUkZmZiXfffRdt2rSRKy4iIrJ21V2zB4BzVhY048axsJmqTHKi06ZNGwiCAFEUjY7/4x//wMcffyxbYEREpGKmDG8BcM/MhMDCZqoGyYnOhQsXjB7b2dmhbt26cHZ2li2ompCUlISkpCTOEiMiUkpls7cA3Pf3h2hnZ5zscN0ekkByohMcHGyOOGpcQkICEhISkJOTU2aBNRER1YCKhrcA5Pv4QLt8OezHj5ensJm9PjbHpERn6dKlGDNmDJydnbF06dIKz504caIsgREREQGAGBcH9O1b/XV72Otjk0xKdBYtWoThw4fD2dkZixYtKvc8QRCY6BARkfy4GSlVkUmJTsm6nEdrdIiIiBRlSmEzV2u2WbKto0NERKSYytbt4WrNNsukHp2pU6eafMGFCxdWORgiIqIqq8nNSNnrYzVMSnSOHTtm9Pjo0aMoKipC06ZNAQBnzpyBRqNBZGSk/BESERHJQa7NSFnUbFVMSnR27dpl+PfChQvh4eGBTz75BLVq1QJQvNFnXFwcunTpYp4oiYiI5MCiZpsjuUbn/fffx7x58wxJDgDUqlUL77zzDt5//31ZgyMiIqox+uEtjab4sdSi5rJkZBTXDWVkmCdmqpTkRCcnJwe3bt0qdfzWrVu4d++eLEEREREpQq6iZoCFzRZCcqLzzDPPIC4uDl9//TUyMjKQkZGBr776CqNHj8azzz5rjhiJiIhqTmAgEB1d9lCUKb0+QPlDXGX07DhnZUHYvZu9PmYieQuIFStW4NVXX8WwYcNQWFhYfBF7e4wePRoLFiyQPUAiIiKLYsIeXaYWNgurV6P32LEQRJGFzWYiOdFxdXXFsmXLsGDBApw/fx4A0KhRI7i5uckeHBERkUWqZI8uUwubNePGFSc5AKezm0mVFwx0c3ND69at0bp1ayY5REREJZlY2CyYUtjMWp9qkdyjAwBHjhzBunXrcPnyZTx8+NDoua+//lqWwIiIiKxaZUNcjRtDtLMzTnaqM52dvT5lktyj88UXXyAqKgrp6enYuHEjCgsL8fvvv2Pnzp3w8vIyR4xERETWqZLCZu3y5dDpZ3FVZzo7e33KJblHZ+7cuVi0aBESEhLg4eGBJUuWICQkBPHx8fD39zdHjERERKokxsUhVaPBE8HBsG/WrMxeH1kXMbTBXh/JPTrnz59Hv379AACOjo7Izc2FIAiYMmUKVq5cKXuAREREapbv4wOxW7eqT2dnr0+FJCc6tWrVMiwMWL9+fZw4cQIAkJ2djby8PHmjIyIisnVyLGIoYV0fw/kqWdFZcqLTtWtXpKamAgAGDhyISZMm4aWXXsLQoUPxxBNPyB4gERGRzavuIoZStq5QWc+P5BqdDz/8EPn5+QCAN998Ew4ODti3bx+ee+45vPXWW7IHaC5JSUlISkqCVqtVOhQiIqLqMWGGV6W1PoAq630kJzq1a9c2/NvOzg7Tpk0zPH7w4IE8UdWAhIQEJCQkICcnh7PFiIjI+lW0iKG+1yc+vrgnp7ytK0xc0RnJyX8nRBa+onOVFwwsqaCgAAsXLkRISIgcl1OFjAzg+HEfNQxvEhGRGlRW6wPIX++TkQGf48cVrfUxOdEpKCjA9OnT0a5dO0RFRWHTpk0AgNWrVyMkJASLFi3ClClTzBWnVUlOBsLC7DFjRmeEhdlb+/AmERGpRUW1PvrnZZzlZR8Whs4zZsA+LEyxWh+Th65mzpyJjz76CD179sS+ffswcOBAxMXF4cCBA1i4cCEGDhwIjf6FsWF/J7oCgOK/yxveJCIisjhy1Pv89cNQv+qzUFGtj5mZnOisX78e//nPf/DUU0/hxIkTaN26NYqKivDrr79CEARzxmhVTB3eJCIisljVrfexoB+GJic6GRkZiIyMBACEh4fDyckJU6ZMYZLzCFML24mIiKyWXLO8aoDJNTparRaOjo6Gx/b29nB3dzdLUNbs7+FNEUDx32UVthMREVk1E9b2Ef8qaRHLm+VVA0zu0RFFESNHjoSTkxMAID8/H2PHjoWbm5vRedy9vDjR7dGjCJ99dhDDh3dESIhDmedZyRIERERE0o0ejaIePXDws8/QcfhwOCg0M9vkRCc2Ntbo8b/+9S/Zg1GTwECgVavb5SYwVrQEARERUdUEBuJ2q1aK/jZvcqKzevVqc8ZhU6QsPElERERVJ8uCgSSNlC1HiIiIqOqY6CjAlIUniYiIqPqY6CjAlIUnS8rIKF6xm9tJEBERScNERyGmbDkCFBctBwcDPXoU/83tJIiIiEzHREdBlW05ImXfNCIiIiqNiY4FY9EyERFR9TDRsWAsWiYiIqoeJjoWTErRMguWiYiISmOiY+FMKVpmwTIREVHZmOhYgYqKllmwTEREVD4mOlaOBctERETlY6Jj5aQULLOOh4iIbA0THStnasEy63iIiMgWMdFRgcoKllnHQ0REtspe6QBIHoGB5a+wXFEdT3lfQ0REpAbs0bEBrOMhIiJbxUTHBrCOh4iIbBUTHRvBOh4iIrJFNlujk5SUhKSkJGi1WqVDqTFy1fFkZADp6QKyspzNEygREZFMbDbRSUhIQEJCAnJycuDl5aV0OIrT1/GUTHbKquNJTtb3/NhDEHpDq9VizJiajZWIiMhUHLoiAKbV8Tw6vCWKAsaP13B4i4iILJbN9uhQaaNHAzExxcNVYWGlh6zKHt4Syp2mnpFR/DWNG3MaOxERKYM9OmSkog1Ey56mLpY5TZ0zuIiIyBIw0SGTPTq8ZWenw7Jl2jKLlTmDi4iILAETHZJEP009NbUIK1emIi5OLHWOlB3VuUAhERGZExMdkiwwEOjWTYSPT36Zz5u6EjOHt4iIyNyY6JDsqjKDi8NbRERkDpx1RWZRtRlc5S9QyNlbRERUFezRIbORPoOLw1tERCQvJjqkCLmHt1jUTEREZWGiQ4qpbKNRU2dvsdeHiIjKw0SHFFXd4S32+hARUUWY6JDFMmV4i70+RERUESY6ZNEqG96Su9eHiIjUhYkOWbyKhrfk7PXR4xAXEZF6MNEhqydHr49eySGusDB7pKY2MFfYRERUA5jokCpUt9cHKGuIS8Dy5REsbCYismJMdMgmVNbrA5Q9xKXT2eH8ecHoGAubiYisBxMdshkV9foAZQ9x2dnp0KjR3zu0czo7EZF1YaJD9JfSQ1wixo371azT2ZkMERGZFxMdohJKDnGdPVuEXr0uGz0v53R2DoEREZkfEx2iR9TEdHYOgRER1QwmOkQSyTGdnSs6ExHVDCY6RFVQ3V4fc6zozJ4fIqLSmOgQmUFlvT5yr+jM4mciorLZKx0AkVoFBpY/lR0oTn5iYooTl7Cw0ufqe31KJjtlrehcXs9PTIzxNZOT/z7Pzg5YvlyAr2/12khEZOnYo0OkIDlWdK5q8fP48RpkZTmXui97fYhITZjoEFkwU1Z0rnrxs4DMTDejYxwCIyK1YaJDZOEqW9G56sXPIvz9cw2Puf4PEakREx0iFahK8fOyZVr4+OQbzuH6P0SkRkx0iFSisp6fR5OhuDjR6Hml1v9hQkRE5sREh8iGWNr6P6wJIiJzY6JDRAY1uf4Pa4KIqCZwHR0iMlJT6/9UlBDpr2nqGkH6c8+eLb5/RfETkW1hjw4RSSbH+j9K1ARxCIzI9jDRISLZmbL+T03XBDEZIrJNTHSIyCwqmwUG1FxNkDnqgZgQEVkHJjpEpCip0+IfTYZM6fU5d06QdY0g9g4RWQ8mOkRk8apbExQWJspWD8ShMiLrwkSHiKyeHENgpvQMAcoMlTEZIqo6JjpEpArVHQKr6dliTIaIagYTHSKyGdVNhvTXkKN3SM5kaPVqgUXUROVQRaKzefNmNG3aFI0bN8a///1vpcMhIitWU7PF5EqGsrKcMW6chkXUROWw+kSnqKgIU6dOxc6dO3Hs2DEsWLAAt2/fVjosIlK5mhgqMyUZysx0h04nGJ3DImqiv1l9onPo0CG0bNkS9evXh7u7O/r06YNt27YpHRYRUY0kQ/7+92FnZ7wTPYuoif6meKLz448/on///ggICIAgCNi0aVOpc5KSktCwYUM4OzujY8eOOHTokOG5a9euoX79+obH9evXx9WrV2sidCKiaqtuMuTjk4/ly7UsoiYqh+KJTm5uLiIiIpCUlFTm819++SWmTp2KWbNm4ejRo4iIiEBMTAxu3rxZw5ESESmjsmQoLk60yiJqrkRNNUHx3cv79OmDPn36lPv8woUL8dJLLyEuLg4AsGLFCnz//ff4+OOPMW3aNAQEBBj14Fy9ehUdOnQo93oFBQUoKCgwPM7JyQEAFBYWorCwsLrNMdBfS85rWhK1tw9QfxvV3j5A/W0s2T5fX8DXV3+87PNfeKE4qTh/XkCjRiICA43P9fUFli8XMH68BlqtAI1GxLJlWvj6iobzGjYE7OzsjeqCNBoRwcFFhnPS0wXodMY/XrRa4NSpIvj6Fg+zFSdDf1+nOBkS0aNHkSH50rfv3//WYcIEETqdADs7EcuXaxEXZzxcp7/muXMCwsJEq9nB3pY+o+a6dmUEURRLf1oUIggCNm7ciAEDBgAAHj58CFdXV2zYsMFwDABiY2ORnZ2Nb775BkVFRWjevDl2794NLy8vREZGYt++fahTp06Z90hMTMTs2bNLHV+7di1cXV3N0SwiIquSleWMzEw3+Pvnwscnv9TzqakNsHx5BHQ6O9jZ6TBu3K/o1euy0de/9FJviOLfyZCdnQ4rV6Yarnf8uA9mzOhc6tpvv70XrVrdlnQtfUzLlrWBKAoQBBHjx6cZxWTcNnf4+98vs22mnkPKy8vLw7Bhw3D37l14enqWe57iPToVycrKglarha/+15S/+Pr64tSpUwAAe3t7vP/+++jevTt0Oh1ef/31cpMcAJg+fTqmTp1qeJyTk4OgoCD07t27whdKqsLCQqSmpqJXr15wcHCQ7bqWQu3tA9TfRrW3D1B/G5VqX9++wCuvaHH+vO6vnqFwAOFG52i12kd6hnR44YUehudbtwZmzRJL9QwNH97RqEdn4cJjRkkOAOh0dggOfgLduv3dO/Tss/aG80RRwIoVbfDKK+FGPTurVwt/TcUvv2fIlHP095Sj94if0arTj8hUxqITHVM99dRTeOqpp0w618nJCU5OTqWOOzg4mOVDZq7rWgq1tw9QfxvV3j5A/W1Uon0hIcV/yjNmTHFCdO4cEBYmIDDQ+MdNSEhxzVB8fPGwVnHNkICQEON26GeVGSdEQLNm9tA3+eLFsuqGBFy65GCIMSMDGDeuZN2QgPHj7dG37991SqacAxTXEulrkOzsittRVl2U/ppnzxbXP1WUEPEzWrVrmkLxYuSK+Pj4QKPR4MaNG0bHb9y4AT8/P4WiIiIiU8ixErUps8rkKqKWs9AaML3YOivLGbt3CxUWWrMYu+osOtFxdHREZGQkduzYYTim0+mwY8cOdOrUScHIiIhIDqasRF3ZrDK5ZpTJlTAB0rbveOml3ujd277cZIhT9atH8UTn/v37SEtLQ1paGgDgwoULSEtLw+XLxUVkU6dOxapVq/DJJ58gPT0d48aNQ25urmEWFhERqV9NLL6oxC7348ZpDLVFZSVDXLeo+hSv0Tly5Ai6d+9ueKwvFI6NjUVKSgoGDx6MW7duYebMmbh+/TratGmDLVu2lCpQJiIi2xYYWPkeZTEx+rqhss+t7Bx9MmRcW1T+Ao0lk52ye4fK3r5Df72KEqaStUVlJUMxMeatLTK1/khpivfoREdHQxTFUn9SUlIM50yYMAGXLl1CQUEBDh48iI4dOyoXMBERWS1Thspqdpf7irfvqOlFHAHTeoek9CAdP+6jaA+S4omOUpKSktCiRQu0b99e6VCIiMjKyLXL/fLlWtjZFWcf5hxOk7O2SMpwWliYPWbM6IywMPsKV742J5tNdBISEnDy5EkcPnxY6VCIiEilTNm+Y+XKVKSmFpXbOyRH/ZGctUXSepD09UdCuT1I5qZ4jQ4REZEt8/HJR7duIipaFqa69Udy1haZXn9kfO1Ha4tqis326BAREalJTdUWydmDVBPYo0NERGQjKusZAuSenSYatgH56CNBkdlZTHSIiIjIiCkJkSnDaT16FOGzzw5i+PCOpbb3qCkcuiIiIiKzCAwEWrW6reg6O0x0iIiISLVsNtHhOjpERETqZ7OJDtfRISIiUj+bTXSIiIhI/ZjoEBERkWox0SEiIiLVYqJDREREqsVEh4iIiFSLiQ4RERGpFhMdIiIiUi0mOkRERKRaNrupZ1JSEpKSklBUVAQAyMnJkfX6hYWFyMvLQ05ODhwclNnIzJzU3j5A/W1Ue/sA9beR7bN+am+jOdun/7ktimKF5wliZWeoXEZGBoKCgpQOg4iIiKrgypUrCKxg11CbT3R0Oh2uXbsGDw8PCIIg23VzcnIQFBSEK1euwNPTU7brWgq1tw9QfxvV3j5A/W1k+6yf2ttozvaJooh79+4hICAAdnblV+LY7NCVnp2dXYWZYHV5enqq8sOrp/b2Aepvo9rbB6i/jWyf9VN7G83VPi8vr0rPYTEyERERqRYTHSIiIlItJjpm4uTkhFmzZsHJyUnpUMxC7e0D1N9GtbcPUH8b2T7rp/Y2WkL7bL4YmYiIiNSLPTpERESkWkx0iIiISLWY6BAREZFqMdEhIiIi1WKiI7N58+ahffv28PDwQL169TBgwACcPn1a6bBks3z5crRu3dqw+FOnTp3www8/KB2W2bz77rsQBAGTJ09WOhTZJCYmQhAEoz/NmjVTOixZXb16Ff/6179Qp04duLi4oFWrVjhy5IjSYcmmYcOGpd5DQRCQkJCgdGiy0Gq1mDFjBkJCQuDi4oJGjRrh7bffrnRPI2ty7949TJ48GcHBwXBxcUFUVBQOHz6sdFhV9uOPP6J///4ICAiAIAjYtGmT0fOiKGLmzJnw9/eHi4sLevbsibNnz9ZIbEx0ZLZnzx4kJCTgwIEDSE1NRWFhIXr37o3c3FylQ5NFYGAg3n33Xfzyyy84cuQIevTogaeffhq///670qHJ7vDhw/joo4/QunVrpUORXcuWLZGZmWn4s3fvXqVDks2dO3fQuXNnODg44IcffsDJkyfx/vvvo1atWkqHJpvDhw8bvX+pqakAgIEDByocmTzee+89LF++HB9++CHS09Px3nvvYf78+fjggw+UDk02L774IlJTU7FmzRocP34cvXv3Rs+ePXH16lWlQ6uS3NxcREREICkpqczn58+fj6VLl2LFihU4ePAg3NzcEBMTg/z8fPMHJ5JZ3bx5UwQg7tmzR+lQzKZWrVriv//9b6XDkNW9e/fExo0bi6mpqWK3bt3ESZMmKR2SbGbNmiVGREQoHYbZvPHGG+Ljjz+udBg1atKkSWKjRo1EnU6ndCiy6Nevnzhq1CijY88++6w4fPhwhSKSV15enqjRaMTNmzcbHW/btq345ptvKhSVfACIGzduNDzW6XSin5+fuGDBAsOx7Oxs0cnJSfz888/NHg97dMzs7t27AIDatWsrHIn8tFotvvjiC+Tm5qJTp05KhyOrhIQE9OvXDz179lQ6FLM4e/YsAgICEBoaiuHDh+Py5ctKhySbb7/9Fu3atcPAgQNRr149PPbYY1i1apXSYZnNw4cP8emnn2LUqFGybkyspKioKOzYsQNnzpwBAPz666/Yu3cv+vTpo3Bk8igqKoJWq4Wzs7PRcRcXF1X1rupduHAB169fN/p+6uXlhY4dO2L//v1mv7/Nb+ppTjqdDpMnT0bnzp0RHh6udDiyOX78ODp16oT8/Hy4u7tj48aNaNGihdJhyeaLL77A0aNHrXq8vCIdO3ZESkoKmjZtiszMTMyePRtdunTBiRMn4OHhoXR41fbHH39g+fLlmDp1Kv7nf/4Hhw8fxsSJE+Ho6IjY2Filw5Pdpk2bkJ2djZEjRyodimymTZuGnJwcNGvWDBqNBlqtFnPmzMHw4cOVDk0WHh4e6NSpE95++200b94cvr6++Pzzz7F//36EhYUpHZ7srl+/DgDw9fU1Ou7r62t4zpyY6JhRQkICTpw4oboMvWnTpkhLS8Pdu3exYcMGxMbGYs+ePapIdq5cuYJJkyYhNTW11G9balHyt+LWrVujY8eOCA4Oxrp16zB69GgFI5OHTqdDu3btMHfuXADAY489hhMnTmDFihWqTHSSk5PRp08fBAQEKB2KbNatW4fPPvsMa9euRcuWLZGWlobJkycjICBANe/hmjVrMGrUKNSvXx8ajQZt27bF0KFD8csvvygdmupw6MpMJkyYgM2bN2PXrl0IDAxUOhxZOTo6IiwsDJGRkZg3bx4iIiKwZMkSpcOSxS+//IKbN2+ibdu2sLe3h729Pfbs2YOlS5fC3t4eWq1W6RBl5+3tjSZNmuDcuXNKhyILf3//Ukl38+bNVTU8p3fp0iVs374dL774otKhyOq1117DtGnTMGTIELRq1QojRozAlClTMG/ePKVDk02jRo2wZ88e3L9/H1euXMGhQ4dQWFiI0NBQpUOTnZ+fHwDgxo0bRsdv3LhheM6cmOjITBRFTJgwARs3bsTOnTsREhKidEhmp9PpUFBQoHQYsnjiiSdw/PhxpKWlGf60a9cOw4cPR1paGjQajdIhyu7+/fs4f/48/P39lQ5FFp07dy61pMOZM2cQHBysUETms3r1atSrVw/9+vVTOhRZ5eXlwc7O+MeTRqOBTqdTKCLzcXNzg7+/P+7cuYOtW7fi6aefVjok2YWEhMDPzw87duwwHMvJycHBgwdrpL6TQ1cyS0hIwNq1a/HNN9/Aw8PDMP7o5eUFFxcXhaOrvunTp6NPnz5o0KAB7t27h7Vr12L37t3YunWr0qHJwsPDo1Q9lZubG+rUqaOaOqtXX30V/fv3R3BwMK5du4ZZs2ZBo9Fg6NChSocmiylTpiAqKgpz587FoEGDcOjQIaxcuRIrV65UOjRZ6XQ6rF69GrGxsbC3V9e38v79+2POnDlo0KABWrZsiWPHjmHhwoUYNWqU0qHJZuvWrRBFEU2bNsW5c+fw2muvoVmzZoiLi1M6tCq5f/++Ua/whQsXkJaWhtq1a6NBgwaYPHky3nnnHTRu3BghISGYMWMGAgICMGDAAPMHZ/Z5XTYGQJl/Vq9erXRoshg1apQYHBwsOjo6inXr1hWfeOIJcdu2bUqHZVZqm14+ePBg0d/fX3R0dBTr168vDh48WDx37pzSYcnqu+++E8PDw0UnJyexWbNm4sqVK5UOSXZbt24VAYinT59WOhTZ5eTkiJMmTRIbNGggOjs7i6GhoeKbb74pFhQUKB2abL788ksxNDRUdHR0FP38/MSEhAQxOztb6bCqbNeuXWX+7IuNjRVFsXiK+YwZM0RfX1/RyclJfOKJJ2rssyuIooqWmiQiIiIqgTU6REREpFpMdIiIiEi1mOgQERGRajHRISIiItViokNERESqxUSHiIiIVIuJDhEREakWEx0iIiJSLSY6REREpFpMdIiIiEi1mOgQker07dsXsbGxhse7du2Cj48PtFqtglERkRKY6BCR6tSvXx9Xr141PO7WrRsePHiAAwcOKBgVESmBiQ4RqU79+vWRkZFheGxnZwcXFxfcvHlTwaiISAlMdIhIdR7t0UlLS0N2djY6deqkYFREpAQmOkSkOvXr18f9+/eRk5MDnU6HKVOmYPjw4fDz81M6NCKqYfZKB0BEJLf69esDADIyMrB69Wpcv34d33zzjcJREZESBFEURaWDICKSU1ZWFurWrYt//vOfOHPmDH788UdD8kNEtoWJDhGpkrOzM0JDQ7Ft2zYEBgYqHQ4RKYSJDhEREakWi5GJiIhItZjoEBERkWox0SEiIiLVYqJDREREqsVEh4iIiFSLiQ4RERGpFhMdIiIiUi0mOkRERKRaTHSIiIhItZjoEBERkWox0SEiIiLVYqJDREREqvX/lykY6jyUqfAAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"critical_radius = lambda nu: 35.26325*(1/(185*nu-211))**(1/2)\n",
"critical_mass = lambda nu: 4/3*3.1415*critical_radius(nu)**3\n",
"\n",
"nu_range = np.linspace(2,10,50)\n",
"critical_radius_range = critical_radius(nu_range)\n",
"critical_mass_range = critical_mass(nu_range)\n",
"\n",
"plt.figure\n",
"plt.title('Critical Mass and Radius as a function of $\\\\nu$')\n",
"plt.xlabel('$\\\\nu$'); plt.ylabel('Radius [cm] and Mass [g]')\n",
"plt.plot(nu_range,critical_radius_range,'.b', label = \"Critical Radius\")\n",
"plt.plot(nu_range,critical_mass_range,'.r', label = \"Critical Mass\")\n",
"plt.yscale('log')\n",
"plt.grid('both')\n",
"plt.legend()"
]
},
{
"cell_type": "markdown",
"id": "f1e97331-57df-4abb-a9e0-5e9cf65df4c4",
"metadata": {},
"source": [
"## What is the minimum value of $\\nu$ necessary to support criticality with this material? (Consider an infinite system)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "87a18296-a54c-4607-9f5a-ad4e1468b3cd",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The minimum value of nu for k_inf = 1 is 1.141e+00.\n"
]
}
],
"source": [
"nu_minimum = 1*sigma_a/sigma_f;\n",
"print(f'The minimum value of nu for k_inf = 1 is {nu_minimum:.3e}.')"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}