{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a0e822ab-6601-47f5-8ec2-be0ede781531",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import astropy.units as u\n",
    "import astropy.constants as c\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddd134bf-1a4f-496b-beea-305f49b73aef",
   "metadata": {},
   "source": [
    "# Stellar Formation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db8d121c-e00f-434f-b0da-220342fe2a0a",
   "metadata": {},
   "source": [
    "In the last lecture we covered 4 important concepts for stellar formation\n",
    "\n",
    "- The equation of motion for collapsing cloud of gas:\n",
    "    $-\\frac{d^2r}{dt^2}=g(r)+\\frac{1}{\\rho (r)}\\frac{dP}{dr}$.\n",
    "- The freefall time:\n",
    "    $t_{\\rm FF} = \\left[\\frac{3 \\pi}{32 G \\rho }\\right] ^{1/2}$.\n",
    "- Hydrostatic equilibrium:\n",
    "    $<P> = -\\frac{1}{3}\\frac{E_{\\rm GP}}{V}$\n",
    "- The Jeans Mass/Density:\n",
    "    $M_{\\rm J} = \\frac{3 k T}{2 g \\bar{m} } R$; \n",
    "    $\\rho_{\\rm J} = \\frac{3}{4 \\pi M^2}\\left[\\frac{3kT}{2G\\bar{m}}\\right]^3$\n",
    "    \n",
    "Using the Jeans Density, we found that for a 300 solar mass cloud, a density of $1\\times 10^6$ particles/m$^3$ is enough to drive a collapse. Clouds with this density typically follow a hierarchical collapse. In this model, one large cloud begins to collapse, driving up the density at various parts of the cloud until the cloud seperates into smaller clouds, with densities high enough to cause graviational collapse. This process continues until a single cloud produces several stars.\n",
    "\n",
    "## An adiabatic collapse\n",
    "To relate our current problem to some of the concepts that were covered in thermodynamics last year, let's consider what happens if the cloud collapses adiabatically (that is, there is no exchange of heat with the surrounding environment). First, if we consider $M$ and $\\bar{m}$ to be fixed, the temperature dependance on temperature is given by\n",
    "$$\n",
    "    T \\propto \\rho^{1/3}\n",
    "$$\n",
    "Now, imagine we have an ideal gas that is collapsing adiabatically. In such a scenario, the equation of state gives\n",
    "$$\n",
    "    T \\propto \\rho^{\\gamma -1}\n",
    "$$\n",
    "where $\\gamma$ is known as the adiabatic index. For a monatomic gas, $\\gamma = 5/3$. The dependence of such a gas is shown in the below figure, and mainly highlights two behaviours. First, if a monatomic gas is initially located in the \"contraction\" phase of the diagram, it will slowly contract until it reachs the solid line - at which case, contraction ceases as the gas is in equilibrium. Alternatively, if a monatomic gas is in the \"expansion\" phase, it will expand until it reaches this equilibrium. This presents us with a problem, as it suggests that a monatomic gas will naturally cease expanding or contracting with relative ease, which is not conducive to star formation.\n",
    "\n",
    "So what if instead $\\gamma \\leq 4/3$? In this case, such a gas will never cross that equilibrium line, which means a gas cloud that starts in the contraction region will continue to contract indefinitely, while a cloud in the expansion region will expand indefinitely. So then the question becomes - what sort of gas do we need such that $\\gamma \\leq 4/3$? Fortunately, this is trivial to do. If the cloud is not monatomic, or if all the hydrogen is either dissociated or ionized, then the adiabatic index is $\\gamma \\leq 4/3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "1bea07ee-9a82-4841-9010-4a59681fd8f6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAGOCAYAAAAZykA4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAABcSAAAXEgFnn9JSAACG90lEQVR4nO3dd3hTR9Y/8O9Vtyxbci+4UUIN1fTeMWAbY5sNBBKS3bRN2GTTgCS/d5O82X130/tGptcAgdBMh9BrMD30ZrngIrk39fn9YSxbXGGMm1zO53n8ZHdm7tWRsKWjuXPncIwxEEIIIYQ8jMDZARBCCCGkaaNkgRBCCCHVomSBEEIIIdWiZIEQQggh1aJkgRBCCCHVomSBEEIIIdWiZIEQQggh1aJkgRBCCCHVomSBEEIIIdWiZIEQQggh1aJkgRBCCCHVomSBEEIIIdWiZIEQQggh1RI5O4Ca4jguE4AcQKqzYyGEEEKamWAApYwx/9oczDWXEtUcxxVKpVK39u3bOzsUQgghpF7oTRbcK9DDYqn8LHaRCBGglEEo4OrtcW7fvg2DwVDEGHOvzfHNZmYBQGr79u27Xr582dlxEEIIIXX2S1Iq/t+mP+BnsdraZg4IwYdR3SAR1e8qgW7duuHKlSu1nplvTskCIYQQ0uyZLFb8a/tVLDuebGsTCTh8PKUbZg4IdV5g1aBkgRBCCGkkuSVGvLb6LE7cybG1eblK8NOscPRv6+nEyKpHyQIhhBDSCK5mFOLFFUlIyyuztT3Zxh0Jz/RFG5WLEyN7NEoWCCGEkAa241IG3v7lAspMFltbdM9AfBrXAy4SoRMjqxlKFgghhJAGYrUyfL3vBr7ff8vWxnHA/IjOeGl4O3Bc/d3x0JAoWSCEEEIaQJHehDfXXcC+q1m2NjeZCN/N6I1RnXydGNnjo2SBEFIvGGNoLvu2ENLQknUlePXns7iTXQzJ/bsg2/q44r9P90FbHwWsVmv1J6gGx3GNPiNByQIhpNasVisKCgqQl5cHg8Hg7HAIaRKK9Gak5Jbg+W5SoJsUAOAuEyHYUw5jbjqu59b9MaRSKTw8PKBUKiEQNHzlBkoWCCG1whhDZmYmCgoKnB0KIU2G2WKFSMChnbfC1iYSchAL6/cD3WAwIDMzE3q9Hv7+/g0+00DJAiGkVoqKimyJgq+vL9zd3SEUNv1V3YQ0BKuVIT2/DPllRlT8FQg4Dm1ULlDJJfX6WBaLBYWFhcjOzkZ+fj5cXV3h7l6rXZxrjJIFQkitFBYWAgA8PT3h5eXl5GgIcR6j2QpNTinKTBZwXPkMgkQoQKiXHC6S+v+YFQgE8PLygtlsRm5uLoqKiho8WaAS1YSQWiktLQUAuLm5OTkSQpynxGDGrexiu/0TXKUidPBVNEiiUFXF315JSUmDPg5AMwuEkFpgjMFiKX9zlEqlTo6GEOfIKTHgXr7e7i4gL4UUAUoZBI1wt0LF357FYgFjrEHXLVCyQAh5bFXfHJvLpjKE1BcrY8go0COnuPIOII7j0EYlg6dr4yXPVf/2KFkghBBCmgizxQpNbilKDGZbm0hQvj7BVdpyP1Jb7jMjhBBC6lGZ0QxNTimMlsoNlVwkQoR6ukIiatlLAClZIIQQQh4hv9SItLwyWKtcgvOQS9BG5QKBoOVfimvZqRAhhDQyZ2zF62wjR460PW+O4yAUCuHh4YEOHTogLi4OP/74Y7PbvCssLAwcx5VvPlZQhpTcUluiwAEIULogyKN1JAoAJQuEEELqyYQJEzB79mw888wzGD16NDw9PZGYmIg5c+YgODgYy5Ytc3aIj02TU4rsovKFjOmpKegZ7IFXn54CHzdpq0oK6TIEIYSQejF//nyMHDnSrq2goABff/01/vnPf+L555+HyWTCiy++6JwAH8OOXXug0RahUG+ytVWsSxC2ktmEqmhmgRBCSINRKpX46KOPbLMKr7/+OrKysqo/yMkK9SYwdz+0adve1uYuEyPMS+7EqJyLkgVCCGkEqampmDNnDtq3bw+ZTAZPT09ERkbi+PHjvLGMMaxZswbTp09Hx44d4erqCjc3N/Tv3x///e9/HZY3/uijj8BxHJYtW4ZLly4hOjoaHh4ecHV1xYgRIxw+DgDs2LED48aNQ5s2bSCVShEYGIihQ4fi448/rtfnP2vWLAwdOhR6vR4LFy7k9ZeWluLf//43evfuDYVCAYVCgYEDB2L58uUOz8dxHMLCwmCxWPDpp5+iY8eOkEqlCA4Oxrx58xxWQdVqtZg/fz66du0KhUIBpVKJjh074tlnn8Xvv/8Oxhiyi/RI1pVg/IDu6BnsAQDwdZNh6fefoUP78uTh0KFDdms0nnvuOWRkZEAsFiM4ONi2YdmDfv75Z3Ach9mzZ9f2ZXSeihr0Tf0HwOWuXbsyQojzWSwWduXKFXblyhVmsVicHU6TAoCVv7VWOn78OPPw8GAAWKdOnVhsbCwbNmwYE4lETCgUsrVr19qNLysrYwCYl5cXGzZsGHvqqafY2LFjmVwuZwDY7NmzeY/74YcfMgDstddeY3K5nHXv3p099dRTrGfPngwAk8lk7NKlS3bH/PDDDwwAEwqFbPjw4WzGjBls3LhxLCgoiPccqjNixAgGgB04cKDacRWPN3r0aLv2rKws1qNHDwaA+fv7s0mTJrGJEycypVLJALA5c+bwzgWAhYaGsj/96U9MoVCwyMhIFhkZaTtm5syZduMLCwtZ27ZtGQAWHBzMYmJiWHx8POvfvz8Ti8XsH//4B9PklLALqXnsQmoeCwwKZgBYfomBMcbYpk2bWFxcHAPA/Pz82OzZs20/CxcuZIwxFhsbywCwbdu2OXz+I0eOZADY0aNHa/rSVutx/g67du3KAFxmtf0Mru2Bjf1DyQIhTcej3qSsVivLy8trdj9Wq7XOr82DyUJBQQELCAhgQqGQrVq1ym7s6dOnmYeHB1MoFCw7O9vWbjKZ2KZNm5jRaLQbn52dzfr27csAsEOHDtn1VSQLANi3335r1/f3v/+dAWDPPPOMXXtISAjjOI6dPn3art1qtT7yg7+qmiYLR48eZQBYQECAXfukSZMYAPbGG28wvV5va8/MzLQ93507d9odU/Fcu3TpwjIyMmztd+7cYSqVigFgt27dsrUvWbKEAWDR0dG839m0e5ls24GTtkShPFkI4SVMd+/eZQDYiBEjHD6/PXv2MABsypQpvL6bN2/a4q0vlCxQskBIk/aoN6m8vDzbm3lz+snLy6vza/NgsvD1118zAOztt992OP6rr75iANhXX31Vo/Pv3buXAWBvvfWWXXtFsjBkyBDeMTqdzvZNvCoXFxfm4eFRo8etTk2ThWvXrtlmOSqcO3eOAWD9+vVz+Lt09uxZ24d8VRWv8969e3nHzJkzhwFgS5cutbV9+umnDAD75ptv7MYW603scnqBXaJwO7uIhYaGPnayYLVaWYcOHZhIJGL37t2z65s3b95j/TvXRGMmC3Q3BCGENKA9e/YAAGJjYx32Dxs2DADw+++/8/rOnz+PPXv2QKPRoLS0FIwxFBUVAQBu3rzp8Hzjx4/ntXl5ecHT0xMZGRl27eHh4Th69Cj+8pe/4K233kK3bt1q/sRqgbH7+xRUueWw4vWJiYmBQMBfRlexhsHR6yMWizFq1Chee8eOHQHA7vmGh4cDAD7//HP4+flh8uTJMHIS3CuwLwTlrZDCXymrzdMDx3F46aWXMHfuXCxduhTvv/8+AMBkMmHZsmWQSqV49tlna3VuZ6MFjoQQ0oCSk5MBAEOGDLFbFFfx069fPwCATqezHWM0GvH000+jd+/emDdvHv773/9i2bJlWL58OTZu3AgAtqThQUFBQQ7b3dzcYDQa7dp+/PFHtG3bFkuWLMGTTz4Jf39/PPXUU1i3bt1DF+nVRcVz9PT0tLVVvD4ffPCBw9eH4zgUFxfbvT4V/P39IRQKee0VpZurLnIcM2YM3nzzTdy7dw8zZsyAp6cnRg4bgu8//QRpmmRwHIcgDzkCVS51qhj5/PPPQyqVYvHixbYkJDExEVlZWYiNjYWXl1etz+1MNLNACKl3SqUSeXl5zg7jsSmVyno/Z8WdC/Hx8XB1dX3ouM6dO9v+91dffYU1a9age/fu+Oyzz9CnTx94eHhALBbjxo0b6NSpk9234aocfTt/mB49euDKlSvYtWsXduzYgYMHD+KXX37BL7/8gkGDBuHgwYOQSCQ1Pt+jnDt3DgDQtWtXW1vF6zN06FC0b9/e4XEP8zjPFSh/Xf/8wotYuXYDDh88gPNJp/DH+TNYqv4Oy1asQvfpf3qs8zni7e2NuLg4/Pzzz/jtt98wduxYLFq0CACaxf4SD0PJAiGk3nEcB5VK5ewwmoSgoCBcv34d8+fPt02FP8qmTZsAAGvWrOFdGrhz5069xieTyRATE4OYmBgAwOXLl/H000/jxIkTWLRoEV599dV6e6x169YBgN2lg4qZkJiYGLz99tv19liOlBnNEHq0wcyX/oaZL/0NBr0ev65ajE8//n948/U5mFkPyQIAvPLKK/j555+xcOFCdOzYEbt378YTTzzh8JJJc0GXIQghpAGNGzcOQGUCUBMVszKOLin88ssv9RPYQ3Tr1g2vvfYaAOCPP/6ot/OuXLkSx44dg1wuxwsvvGBrr83rUxv5pUbc1pbAVKVipL+nO/7vH+8jICAAWq0W2dnZ1Z6jYpbFbDZXO27YsGHo1q0bNm/ejM8++wxWq9XuOTdHlCwQQkgDevnll+Hr64vPPvsMCxYs4G2oZDabsXv3brsP5ooFemq12m7shg0bsGLFinqJq7S0FN999x3y8/Pt2q1WK3bt2gUACA4OrvPjFBQU4OOPP8bzzz8PAPjhhx/g4+Nj6x8wYADGjRuHY8eO4bXXXkNhYSHvHBcuXLDF9LgYY8goKMOyn9fj/JnyRZJVC0GdO3cWWVlZUCgUj5wN8/b2hlgsxu3btx+5puPll1+G0WjEjz/+CLFYjOeee65W8TcVdBmCEELqWdXV/iqVClu2bEFUVBRefvll/POf/8STTz4JDw8PZGZm4uzZs8jPz8emTZvw5JNPAgDmzp2LXbt2Yf78+Vi/fj06duyImzdvIikpCe+88w6++OKLOsdoNBrxxhtv4J133kF4eDjCwsJgNBpx+vRppKamIiwsDC+99NJjnfM///mPbVvn4uJipKWl4dy5czAajXB3d8cPP/yAZ555hnfcqlWrEBERgf/+97/4+eef0atXLwQGBqKgoAAXL15Eamoq3njjDURERDxWPFbGkJxTiiK9CUknj2L1YjX8/APRu08veKpUuHfvHo4cOQKr1YqPP/74keszJBIJIiIikJiYiJ49e6JPnz6QSCQYMmSILRmq8Oyzz2L+/PkoLS3FlClT4Ovr+1ixNzWULBBCSD3R6/UAwFvIOHDgQFy6dAlff/01tm/fjkOHDgEAAgICMGLECEydOhVjx461jR8+fDiOHj2KDz74AOfOncONGzfQvXt3/Prrr+jTp0+9JAsKhQI//vgjfvvtN1y4cAEXL16ERCJBSEgIXnjhBcyZM8furoWa2L17N4DyhYdubm7w8vJCZGQkxowZg1mzZsHd3d3hcb6+vjh+/DgWLlyItWvX4ty5czh+/Dj8/PzQrl07vP7665g+ffpjP8fcEiOK7heCip72NKRiMS4mncDZpCQUFBTA398fkyZNwhtvvIExY8bU6JyLFi3CO++8g7179+Lnn3+GxWKB2WzmJQtKpRJ9+vTB0aNHm/XCxgrcw1bUNjUcx13u2rVr18uXLzs7FEJaPavViuvXrwMAOnXq9Nir0luqK1euoFu3bujatSvovcp5CstMSM0thaXK55u7TIxgT3mjVYxMTU1F27ZtERwcjDt37jRIOevH+Tvs1q0brly5coUxVqvNNGhmgRBC6gFjDD/88AMA8Mo0k8bBGIO22IDMAr1du5+7DL5u0gb5wH6Y//znP7BYLHjttdca9XEbCiULhBBSBzqdDm+88QYuXbqES5cuwc3NrcFvASR8VitDWl4p8stMtjYBxyHYUw6li7hRYrh+/To+//xz3L17F/v370dQUBBeeeWVRnnshkZzh4QQUgfFxcVYs2YN0tLSEB0djSNHjqBdu3bODqtVMZotuK0ttksUJCIB2vsqGi1RAMq3l168eDFOnDiB4cOHY/v27VAoFI32+A2JZhYIIaQOwsLCeLdDksZTbDAjJacU5ir/BgqpCCGecoiEjft9eOTIkQ/dWbO5o2SBEEJIs8MYQ26JEffy9WCwLwQVoJS1iHUCTQklC4QQQpoVK2O4l1+G3JLKwlgcx6GNygWervVXy4JUomSBEEJIs2GyWJGSU4oSY+WWy2KhAKFecsgl9JHWUOiVJYQQ0iyUGs3Q5JTa1XeQS0QI9ZJD3MjrE1obShYIIYQ0eXmlRqTnlcFaZQGhp1yCQA8XCGh9QoOjZIEQQkiTxRhDZqEe2iKDrY0DhwCVDF6uElrI2EgoWSCEENIkmS1WpOaV2eo7AIBIwCHEUw6FrPH2TyCULBBCCGmC9CYLNDklMJgr1yfIxEKEeckhEQmdGFnrRMkCIYSQJsVRISilixhBHo1XCIrYo+WjhBBSTziOs/2cOHHioeN++eUX27iwsLDGC7CJY4whu1CP5JwSu0TBz12GkCoVIz/66CNwHIdly5Y5KdLWh5IFQghpAKtXr35o36pVqxoxkvrz3HPPgeM4HDx4sN7PbbEypOSWIrOwsmKkkOMQ5uUKP3fakdHZKFkghJB6JBQK0b17d6xbtw5ms5nXn5OTg127dqFPnz5OiK5pqigEVeCgEJS7g0JQc+bMwdWrVzF16tTGDLNVo2SBEELq2cyZM6HT6bB7925e37p162AymTBr1iwnRNb0FOtNuJVdDL3JYmtTSEXo4KOATOx4IaO3tzc6d+4MpVLZWGG2epQsEEJIPXv66afBcZzDyw2rVq2CQqHAlClTqj3Hjh07MG7cOHh4eEAmk6FTp06YP38+8vPzeWOrXsO/dOkSoqOj4eHhAVdXV4wYMQLHjx/nHaPX67F48WJMmTIF7dq1g4uLC1QqFYYPH461a9fyxnMch+XLlwMARo0aZbc+Izk52TbObDbj+++/R3h4OBQKBRQKBfr374+ffvoJFktlQsAYg67YgFGjRqFbGxXSU1Owf9d2PD91PLqH+cHXxxszZsxAWlpatc/3QSaTCWq1GkOHDoVKpYKLiws6dOiA559/HmfOnKn2Na9Kp9Phr3/9KwIDA+Hi4oInn3wSP/74IxhjDteaPO7rCQBGoxH//e9/0a9fP3h5eUEulyMsLAyRkZEPPcZZ6G4IQgipZ8HBwRg+fDi2bt2K4uJiKBQKAMCdO3dw4sQJPPPMM5DL5Q89/t///jfef/99iEQijBgxAt7e3jh27Bg+/fRTbNq0CYcPH4afnx/vuKSkJLz22mto3749JkyYgGvXruHw4cMYM2YMTp8+jSeffNI2Njk5GS+88AICAwPRqVMn9O/fH5mZmTh+/DiOHDmCa9eu4aOPPrKNnz17No4ePYrbt29jwoQJ8Pf3t/VVPD+LxYIpU6Zgx44dcHd3x7hx48AYw/79+/Hqq69i79692LBhA8BxuJdXhtzSykJQv6xYjJULf8SwYcMwadIknDp1CmvXrsWZM2dw4cIFuLi4PPJ1LykpwaRJk3D48GG4urraEobk5GSsXr0aSqUS4eHhjzyPTqfD4MGDcfPmTQQGBiI6Ohp5eXl48803cfPmTYfHPO7rCZTPQG3YsAFubm4YNmwY3N3dkZ6ejqNHj6K4uBjTp09/ZKyNhjHWLH4AXO7atSsjhDifxWJhV65cYVeuXGEWi8XZ4TQZAJhQKGSMMbZw4UIGgC1fvtzW/7//+78MANu9ezfLyMhgAFhoaKjdOX7//XcmEAiYQqFgJ0+etLXr9Xo2bdo0BoDFxcXZHfPhhx8yAAwA+/bbb+36/v73vzMA7JlnnrFr1+l0bO/evcxqtdq137lzh4WFhTGBQMDu3r1r1zd79mwGgB04cMDh8//iiy8YANatWzeWmZlpa7937x7r1KkTA8C++fY7djOriF1IzWMXUvNY34FDGAAml8vZ8ePHbceUlJSwwYMHMwBs8eLFDp/v0qVL7dr/8pe/MABs+PDhLDs7264vMzPT7vWsTsV5oqOjWVlZma39zJkzTKlUOvx3e9zX886dO7bz6HQ6u2PKysrsXouHeZy/w65duzIAl1ltP4Nre2Bj/1CyQEjT8ag3KavVyvJLjc3u58E3+sdVNVnIy8tjUqmUjR8/3tbfqVMnFhAQwMxm80OThWeffZYBYO+99x7v/FlZWczFxYUJBAKWkpJia6/48BwyZAjvGJ1O5/BxqlOR6Hz33Xd27Y9KFkJCQmzJ0IO2bt3KALCQtu1sicKF1DzWf/BQBoB98MEHvGM2bNjAALDZs2fbtTtKFtLT05lQKGRSqZQlJyfX+Lk+qKioiMlkMiYUCnnJEmOMffDBB/Xyep46dYoBYDExMbWOtTGTBboMQQipd4V6M3p+vMfZYTy2Cx+Oh9LB6vvaUKlUmDx5MrZs2YLMzEykpqbi+vXrePPNNyEUPnwHwiNHjgAon6J+kK+vL8aPH48tW7bg2LFjvGnq8ePH847x8vKCp6cnMjIyHD7e0aNHcfDgQaSnp0Ov14MxZhv7sCl3R1JSUpCSkgIfHx+HcQwZNR5uSiVS7t6BLjsL3r5+8JRL4HJ/EaOjYzp27AgAD429qoMHD8JisSAyMhKhoaE1jvtBZ86cgV6vx8CBAx3ugfHUU0/hX//610OPr+nr2blzZ7i6umL79u34/PPPMXPmTAQGBtY67oZGyQIhhDSQWbNmYePGjVi7di3u3r1ra6vOvXv3AOChmzVVtKenp/P6goKCHB7j5uaG3Nxcu7aCggLExsZi//79D42lqKio2lirqoj7wQ9qxhgyC/TQFhsQ2CYY1wsKoM3MRI+OYfBylVQbu5ubGwDAYDDw+h6UmpoKAGjfvn2NY3ak4oM9ODjYYX9ISIjD9sd9Pd3d3bFw4UK89NJLmDt3LubOnYuOHTti1KhReOaZZzBkyJA6PIv6R3dDEEJIA5k0aRJUKhVWrFiBdevWoUuXLnXeX6G6zYkEgpq/pc+bNw/79+/HiBEjcPDgQeh0OpjNZjDGbLd8MsYecZbqmS1W3NWVQFtsqAgeABCoksFbIbV7Lo8Te1NUm9dzxowZuHPnDhYuXIhp06YhPz8fCQkJGDp0KN5++21nPI2HopkFQki9c5eJcOFD/rRyU+cuq9+3RKlUimnTpmHhwoUAgNdff/2RxwQGBuLu3bvQaDTo2rUrr7/iNsU2bdrUKbZNmzZBKBRi69atcHd3t+u7c+fOY5+vYgpdo9EAKC8ElZxTAmOVQlCZ6eXf/p9oV/vLBA9TMRNw+/btOp0nICAAQOVMxYMe1l7b19PHxwcvvPACXnjhBVti8dRTT+Grr77Cn//8Z3Tr1q2Wz6R+UbJACKl3HMfV27X/5u6ZZ57Bxo0bwXGcw3UIDxo2bBju3r2LNWvW4JNPPrHr02q12L17NziOq/M0dV5eHtzd3XkfbEB57QpHJJLyywaOdqYMCQlBSEgIUlJSsHXHbrTrORDWKt+kzxz5DQX5+ejQoYPdbZf1ZeTIkRAKhdi9ezdSU1MfehnhUcLDwyGTyZCUlISUlBTeZYeHvTa1eT0fxHEcIiIiMHnyZKxZswaXL19uMslC8573IYSQJm7YsGHQ6XTQarU1Wnj32muvQSAQ4LvvvkNSUpKt3Wg04m9/+xvKysoQGxtb6w/DCh07dkReXh7WrVtn1/7111/jwIEDDo+pmD24fv26w/45c+YAAN59523odFpbu6AsH59/8j8AgDfeeKNOcT9MYGAgnn32Wej1esyePRs5OTl2/dnZ2Th16tQjz6NQKDBz5kyYzWa88cYbduslLly4gO+//97hcY/7ep47dw4bN26E0Wi0a8/NzbXFWdd/43pV29soGvsHdOskIU0G7bPgGKrcOvkoD7t1kjHG/vWvfzEATCQSsbFjx7Lp06ez4OBgBoA98cQTdnsYMPbwfQcqhIaGsvK3+0qrVq2y7c0wbNgwNmPGDNa1a1cmEAjYm2++6fCWxaSkJMZxHJPJZGzKlCnsL3/5C/vLX/7CdDodM1us7HZWARs6aiwDwNzc3dmYiEgWGRXN3NzcbLcJPvj7MmLECAbA4W2Kd+/eZQDYiBEjavR8CwsLbXszuLq6sokTJ7KnnnqKDRw4kEkkEvbGG284fH0epNVqWYcOHRgA1qZNG/bUU0+xCRMmMLFYzObMmWP7d6jL67lp0yYGgCmVSjZmzBg2c+ZMNnnyZNtrFRUV9cg4G/PWSZpZIISQJub999/Htm3bMGLECJw+fRobN26EVCrF3LlzcerUKYe7Nz6umTNnYvv27Rg4cCDOnz+PnTt3IjAwEPv370d0dLTDY8LDw7Fq1Sp07doVe/bsweLFi7F48WLo8vJxW1uMYqMV3y5Zg7kf/wfBIW1x4vB+7P9tHzp16oQff/wRGzZsaNCFjG5ubjhw4AC+/fZbdOvWDUeOHMHWrVuh1Woxc+ZMPPvsszU6j7e3N44fP46XX34ZFosFmzdvRkpKCj7//HPMnTsXQPktqVU97us5cOBA/POf/0R4eDiuX7+O9evXIykpCT169MCSJUvw66+/1v0FqUccq+Nq18bCcdzlrl27dr18+bKzQyGk1bNarbap6E6dOjX7leyk9or1JmhyS2GxVn6WuMnECPZ0gagF/l6sXbsWM2bMwCuvvIKffvrJqbE8zt9ht27dcOXKlSuMsVotgmh5/5KEEEIaHGMMuiID7ursEwUfNynCvOTNPlFwVHTq/PnzePfddwE8er+MlobuhiCEEPJYrFaG9Pwy5FUpBCXgOAR5uEAll1RzZPMxZMgQ+Pv7o0uXLnB3d8fdu3dx5swZWK1WzJkzp8ltmtTQKFkghBBSYyaLFZqcUpQaK2+fFAsFCPWSQy5pOR8p7733Hnbs2IGkpCTk5+dDoVBg+PDheOGFF2p0C2xL07zniRpI1TrtD/t57rnnnB1mk1BdXXlCSMtSYjDjVnaxXaLgKhGhg6+iRSUKAPDhhx/i1KlT0Gq1MJlMyMvLw4EDB1plogDQzEK1Zs+e/dC+oUOHNmIkhBDiXLklRqTnl9ltWezpKkGgygWCaragJi0DJQvVoG/LjzZnzhxMnz7dtkUqIaRlYYwho0APXXHl5kQcOASqZPBSSJ0YGWlMlCyQOvH29oa3t7ezwyCENACzxYqU3FIUGyovO4gEAoR4yaGQ0sdHa0JrFurBvHnzwHEc/vSnP/H6dDodAgMDIRQKcfToUVv7c889B47jcPDgQezcuRNDhw6FQqGAh4cHYmNjce3aNd659Ho9Fi9ejClTpqBdu3ZwcXGBSqXC8OHDsXbtWoexVX2cw4cPY/To0XBzc4O7uzsmT56MK1eu8I5hjGH16tUYOnQo/Pz8IJPJEBwcjLFjx+LHH3+0G1vdmoWcnBy8++67eOKJJyCTyeDp6YmIiAjs2bPHYawcxyEsLAwWiwWffvopOnbsCKlUiuDgYMybN69GZWoJIfWjzGTBLW2xXaLgIhaig6+CEoVmqK57KlGyUA8++eQT9O7dG+vXr8fy5cvt+l588UVkZGTgvffec7jOYf369Zg8eTKMRiOioqIQGBiITZs2YeDAgbhw4YLd2OTkZLzwwgtISkpCWFgYpkyZgl69euHkyZOYMWMGPvroo4fGmJiYiNGjR6O0tBSTJk1CQEAAduzYgeHDhyMzM9Nu7Ny5czFr1iwkJSWhZ8+eiI2NxRNPPIGLFy/i888/r9Frkp6ejv79++OLL76A0WhETEwMevfujX379mHChAn4+uuvH3rs008/jX/+85/o1KkTxo8fj6KiInz22Wf4y1/+UqPHJoTUTUGZEbezi+0qRqpcJGjvo4BERB8bzcndu3fx3nvv4caNG3U7UW33iW7sHzRibQjc39/7cVy5coW5uLgwNzc3dufOHcYYYwsWLGAAWL9+/ZjJZLIbP3v2bNvjLFiwwNZutVrZvHnzGADWq1cvu2N0Oh3bu3cvs1qtdu137txhYWFhTCAQ8PZXr3gcgUDANm3aZGs3m80sLi6OAWD/8z//Y2svKytjUqnU7nlUMJlM7PDhw3ZtD9ujPTIykgFgTz/9NDMYDLb2I0eOMLlczoRCITt37pzdMRWvR5cuXVhGRobd81OpVAwAu3XrFiPOR7UhWiar1coyCsrYhdQ8u5+swjLe+w5xvof9HZpMJrZ582YWERHBOI6zvbcCuM2oNkT9q+7Wyc2bN9uN7dKlC7744gsUFRVh5syZuHr1Kt58803I5XKsWrUKIpHjabvBgwfjxRdftHvMTz75BEFBQTh//rzdpQsvLy+MHTsW3AMrj9u2bYsPPvgAVqsViYmJDh9nxowZiImJsf1/oVCI9957DwBw+PBhW3thYSEMBgPat2+Ptm3b2p1DJBJh2LBhD3/B7rtz5w62bdsGhUKB77//3lbWFii/i+SVV16BxWLhXdKo8N1339mVsG3btq1tt7QjR4488vEJIY/PYmXQ5JQiu1BvaxMKOIR5ucLXTcZ73yFNT1paGj766COEhYUhJiYGu3btevDyQ63rxtOFp2pUd+vkgzXOAeDVV1/Fjh07bMVESkpKkJCQgI4dOz70PNOnT+e1icVixMfH45tvvsGRI0d4ly+OHj2KgwcPIj09HXq9vny1ckYGAODmzZsOH2f8+PG8toq4Ko4FAF9fX1uiMn/+fLz00kto167dQ+N3pCLBiYiIgKenJ6//mWeewVdffeXwg18sFmPUqFE1ipUQUj8MJgs0uaXQmyy2NqlIiFAvOWRioRMjIzVRUFCAP/3pT9i0aROsViuvPyAgAGVlZcjPz691xkfJQjVqc+vk4sWLERYWhsLCQkycOBEvvfRSteMfVt8+LCwMAHDv3j1bW0FBAWJjY7F///6Hnq+oqMhhe1BQEK/Nzc0NAHgLB5cvX47p06fj008/xaefforQ0FCMGDEC06dPx8SJE6t9PlVjrngOD6poT09P5/X5+/tDKOS/OT0sVkJI3RTpTUh5oBCU+/1CUMJmXt+hJTOZTNBqtSgrK8PNmzexbds2u0RBKBRiwIABiI2NxbBhw/DUU08hPz+/1o9Hvwn1bOvWrdDry6fxrl+/juLi4no797x587B//36MGDECBw8ehE6ng9lsBmMMu3fvBvDwFa+PUxVw9OjRuHXrFlavXo1nnnkGVqsVK1aswKRJkxAfH1/n51HddCZVLySkcTDGoC0yIFlXYpco+LpJEeolp0ShCWKMobCwELdv38bFixdx79493kyCj48Pnn32WWzatAnfffcdRo4c6fAL2OOimYV6dPPmTds6hYiICGzcuBGvv/46lixZ8tBjNBpNte2BgYG2tk2bNkEoFGLr1q1wd3e3G3/nzp16eAaV3N3d8fTTT+Ppp58GAJw8eRLTpk3Dr7/+ih07dmDSpEkPPbYi5oc9t+TkZABAmzZt6jVmQpztwUSY4zi4ubmhS5cumD59Ol577TWIxbW+bFwnYWFh0Gg0YIy1ikJQLYnZbEZOTg60Wq3ty+iD+vTpg6ioKIwaNapBfscodawnZrMZs2bNQklJCb766iusWrUKnTt3xtKlS/Hrr78+9LhffvnF4bkqjqm6XiEvLw/u7u68ROFh56lPAwcOxDPPPAMA+OOPP6odWxHzrl27HE57rVq1CgBqtFiSkOZo9uzZmD17NmbOnIkePXrg9OnTePPNNzFx4kSYzeZHn6ABmcxW3NEV2yUKEqEA7X1cKVFoQhhjKC4uxt27d3HhwgWkpqbyEgWBQACRSIT27dvj+++/x/jx4xssGaVkoZ787//+L37//XdERUXh5ZdfhouLC1atWgWxWIyXXnrJbu1BVUePHuXNPHz44YdISUlBjx497D5QO3bsiLy8PKxbt85u/Ndff40DBw7Uy/NISUnBsmXLUFpaateu1+ttjxEcHFztOdq1a4fJkyejqKgIb7zxBkwmk63vxIkT+OmnnyAUCvHaa6/VS8yENDXLli3DsmXLsHLlShw5cgTHjx+HTCbDb7/99tAN1BrLTW0xSo2VCxldpeWFoFxaWCGo5spisSA7OxtXrlzBtWvXkJOTw7u87OLigsDAQHTo0AFisRhSacNvu02/HdWorrJkSEgI/vd//xdA+Qfg//3f/8HPzw+LFy+2jQkPD8dHH32EDz74AM899xx2797Nm6b861//ihdeeAEJCQlo3749Ll68iMuXL8Pd3Z23wPK9997DrFmzMH36dPz4448ICgrChQsXcO3aNbz55pvVbnRUU7m5uXj++efx2muvoW/fvggKCkJJSQmOHz8OrVaLvn37IjY29pHnSUhIwLBhw7BixQocOnQIgwYNglarxcGDB2GxWPDll1+iV69edY6XkOZgwIABeO6556BWq7F7927brcCNqWJZgtlSeY3bSyFFgFJGhaCagJKSEmi1WuTm5jq8o0EoFMLd3R0eHh625MDRuIZCMwvVWL58+UN/tm7dCgAoLi7GrFmzYLFYsGTJEvj4+NidY/78+Rg6dCj27t2Lb775hvcYf/rTn7B161YIhUJs2bIFaWlpmDJlCk6cOIHevXvbjZ05c6bttszz589j586dCAwMxP79+xEdHV0vz7l9+/b48ssvMXLkSKSkpGDjxo04evQoQkND8fXXX+PQoUM1ymLbtGmD06dP4+2334ZIJMLGjRtx5swZjBkzBrt378Zbb71VL/ES0lx069YNAJCdnc3rY4xhzZo1GD16NDw8PCCTydClSxd89NFHvFk+ALh16xY++ugjDBo0CP7+/pBIJAgKCsKzzz7L26nPyhju5ZfBUuWDheM4tPFwQZtHVIwsKSnB/PnzERYWBplMhg4dOuCTTz6ByWRCWFgY78tPxfOYPn06OnbsCFdXV7i5uaF///7473//6/DDjT3G9vItjcVigU6nw9WrV3H16lXodDreaySTyeDv748OHTrA39+/UWYRHOEetnq+qeE47nLXrl27Xr582dmh1IvnnnsOy5cvx4EDBzBy5Ehnh0PIY7Farbh+/ToAoFOnTnQXy30VH56O3lf//e9/4/3338esWbOwcuVKW7vVasWsWbOwZs0aKBQK9O3bFx4eHkhKSkJqair69++PgwcPwsXFxXbM/Pnz8dlnn+HJJ59EaGgopFIprly5gqtXr8Ld3R1HjhxBjx497ApBTRzUA/fSUnE5vQChXnK4PqK+g8FgwMiRI3Hy5El4enpi1KhRMBgM2L9/P8aPH49z587ZFkxW0Ov1cHFxgZeXF7p27YrAwEDk5OTg+PHjKC0txezZs3kzpu+++y6++OILSKVSDB8+HN7e3sjMzMQff/wBuVxuWxDdkpSVlUGr1SInJwcWi4XXLxAI4ObmBg8PD7t/9wdZrVakpKQ8dNFjVfHx8UhOTk5jjFV/Hfkh6DIEIaT+MQboC5wdxeOTKYEGmpLftWsXgPLNyqr68ssvsWbNGowcORJr1qyx7V5qNBrx6quvYvHixfj444/xn//8x3ZMTEwMXn75Zd4uq0uXLsWf//xn/P3vf8f2XXuhySmB0WL/TbWDb83qO3z++ec4efIk+vfvj927d0OlUgEov5tp+PDhSE1N5R0jEomwadMmTJ482W6hnVarxaRJk7B8+XL8+c9/xvDhwwGUJxfff/893NzccOHCBbvnYzabceLEiUfG2VxYrVbk5eVBq9U+9JZ6qVQKlUoFlUrV5BJwmllwEppZIM3ZI2cWyvKBTx1vONakzdMALqpaH/7gzILVasXdu3fxxRdfQK1WY8qUKdiwYYNt+3ez2WzbXe/27dvw8/OzO19ZWRnatWsHg8EAnU5Xow+QoUOH4vjx4zh2ORmubpV3Tk0e3BNpqSk1rj4YFBSE9PR0h7vILlq0yLZNfU3Pt2/fPowbNw5vvfUWvvzySwDll2T8/PzQq1cvnDt3rkbnaW70ej10Op1tX5wHVdxe6+HhAblc/ljnppkFQghpxhxtPPbiiy8iISHBru/s2bPQ6XQYN24cL1EAyle9h4eHY/v27bh58yY6depk6ysuLkZiYiLOnz+P3Nxc211Haen3wBhDSvJddOneExwAf6UMQkHNZ0w0Gg3S09Ph7+/vsFruU089ZVfT5kHnz5/Hnj17oNFoUFpaCsaYbXfZqlvS18f28k2R1WpFQUEBtFotCgsLHY6RSCS2WYT62DSpoVGy4CQVt1YRQlqeiroyer3edsfSwoULMXjwYLu7rCqux+/du/eRhZp0Op0tWdi/fz+mT58OrVb70PElxcUQCjiEeMrhJnu8e+8rarA87DZpNzc3qFQq3j4qRqMRzz33HNasWfPQcz+4JX1dt5dvSipmgHQ6nd0t4xU4joNCoYCHhwdcXV2dEGHtUbJACKl/MmX5lH5zI1PWy2ke/CLw+eefY+7cuXjttdcwatQoW02YipXvHTp0wJAhQ6o9p5eXF4DyGYU//elPyM3NxT/+8Q9Mnz4d/oFByCplMFqsmD/nBezc8itEQg4dfBSQNmIhqK+++gpr1qxB9+7d8dlnn6FPnz7w8PCAWCzGjRs30KlTJ95li4rt5bdt24Zdu3bh4MGDWLFiBVasWIG4uDhs2LCh0eKvjYotmLOzs1FQ4Hidjlgsts0iPKwCcVPXPKMmhDRtHFena/8tzbvvvot9+/Zhz549+Pjjj20bsVUUeOvcuXONZxqPHDmCnJwcxMfH4+OPP+YVgkpLKU/S2qhktU4UAgICAMDhIkagfHbA0e6smzZtAgCsWbPGdqtoheq2pK/L9vLOYjKZoNPpoNVqYTQaef0cx8HV1RUeHh5QKBROiLB+Na3lloQQ0kJV3M2wcuVKW92Ufv36QalU4tChQ8jNza3RefLy8gCU72WiLdLbFYJKuXsH1/64AAB1KgQVGhqKNm3aIDMzE8ePH+f1r1+/vtrYHFW5fZwt6R9ne/nG9GAhp/T0dF6iIBKJ4OXlhfbt2yM4OLhFJAoAJQuEENIoevfujZiYGJjNZnz22WcAym+Vmzt3LoqKihAbG+vw23d6errdvgwdO3YEAGz4dSMu30lDxaR+cWEB/v3+3x1eK6+NV155BQDw9ttv202vazQa2+61D6qITa1W27Vv2LABK1as4I2vj+3lG4PZbLbt/XDjxg3k5eXxLqe4urqiTZs26NChA3x9fZ1WMKyh0K2ThJDHRpsyOVbdpkwAcOHCBfTu3RtSqRR3796Fv78/rFYrnnvuOaxcuRISiQS9e/dG27ZtYTQacf36dVy5cgU9evTA+fPnAQBGsxWjx47FsUMH4KZUou/AoRBwHM6cPApvb29069YNW7Zs4d2WXbXqZE1U3ZTJy8vLblOmMWPG4Pz588jIyLD7Zn348GGMHj0aFosF4eHh6NixI27evImkpCS88847+OKLLzBixAgcPHgQQPldE71794ZcLn/o9vJHjx51yq6FjDGUlJQgOzvbYXIAlM8iVGzBLJE0fhGuxrx1kv7CCSGkkfTs2RNTp06FXq/HV199BaB8t74VK1Zgy5YtGDduHO7evYtff/0VR48ehUwmw7vvvmtb41BiMONWdjG+WrgaL/7tbXh4euPYwX24cfkCpk+fjpMnT9o2T6orqVSKvXv3Yu7cuXB1dcXWrVtx+fJlvP3221i3bh2ysrJsiy4rDB8+HEePHsXo0aNx584dbNu2DRKJBL/++qvDwnH1tb18fTKbzXaFnHJzc3mJglwutxVy8vPzc0qi0NhoZoEQ8thoZqHx5ZQYcC9fb/fB5a2Qwt8JhaBOnjyJQYMGISIiAjt37mzUx24otSnk5Gy0KRMhhBAA5YWgMvL1yCkx2No4jkMblQyerg37oXXu3Dn07NnTLhm8c+cOXn75ZQBwSvXM+mSxWJCbmwutVuuwYBdQvjGWSqWCUql85F4YLRklC4QQ0kSZLVZocktRYqjcJlgkENSoEFR9mDFjBgoKCtC9e3d4eXkhJSUFZ86cgcFgQHR0tO1Wx+amtLQUOp2u2kJOFbMIMpnMCRE2PZQsEEJIE1RmNEOTU2pXCEouESLE07VGhaDqw9/+9jf88ssvuHjxInJzcyGTydCrVy/MmjULr7zySrP6pl2TQk4ymcw2i0CX1uxRskAIIU1MfqkRaXllsFZZn+Ahl6CNygWCx6jxUFevvfaaw4WJzYler7eVg3ZUyKmm5aBbO0oWCCGkiWCMIatQj+yiKusTAPgrXeCtkDSrb/LOZLVakZ+fD61Wy6tFUUEqlUKpVDabQk7ORskCIYQ0ARarFam5ZSjUV26qVNtCUK1VTQo5VRTBam6FnJyNkgVCCHEyg8mC5JxSGMyVi+1kIiFCveSNWgiqOWKM2cpBP6yQk0Qisc0iNNdCTs5Grxoh5LFVnQ63Wq20GKwOCvUmpFYpBAUA7jIxgj3lEDbi+oTmxmg02mYRWkMhJ0cYYzXekbOuKFkghDw2juMgkUhgNBpRUlICpbJ+Sju3JowxaIsNyCyw31DH110GPzcprU9woKKQk1ardVj1EigvB61UKuHh4dHiZxHKysooWSCENG1ubm7IyclBVlYWgPJCOjTDUDNWK0N6Xhny9ZXfiAUchyAPFyhdJI36jbE5MJlMts2THM0iAOVbMHt4eMDV1dWWaDnaibElYIyhrKwMubm5Du/waAiULBBCasXLywslJSXQ6/W4d++es8NpNhhjMJitqJoLcBwgEgmQWcQh03mhNTkWiwVmsxkWi8Vh8sRxHIRCIUQiERhjyM3NrXGp7+asIpk0m82NlhBRskAIqRWhUIiQkBDk5OSgqKjood/4SKUSgxma3FKYLZUffApZ+UZLjV3foakym83IycmBVqt9aM0DV1dXqFQquLm5AYDDXRhJ/aJkgRBSa0KhEL6+vvD19aWp80dY83sK/rntit1CxmcHhWHuyE4QCVv35RvGGJKSkrBw4UKsX78eBoOBN0apVGLcuHGYOnUqQkNDnRBl60bJAiGkXnAcR4vyHDCarfgo8TJ+PpVia5MIBfjX1CcxrW+tCgC2GEVFRfj555+hVqtx/vx5h2N69OiBKVOmYMKECVSnwYkoWSCEkAaiLTLg1dVncDo5z9bm6yaF+plw9AnxcGJkznX+/HkkJCRg1apVDus0uLm5YcyYMZg2bRo6derkhAjJgyhZIISQBvBHegFeWpGEe1VujewVrELCM+Hwc29935DLysrwyy+/4KeffsKpU6ccjunSpQuio6MxefJkyOXyRo6QVIeSBUIIqWdbzqdj7oaLMJgrV6rH9QnCv6Y+CVkr25Hx6tWrSEhIwPLlyx3ujSCXyzFq1ChMmzYN3bp1o0tZTRQlC4QQUk8sVobPdl9DwqE7tjahgMP7k7rgz0PCWs0HocFgwKZNm6BWq3Ho0CGHY9q3b4+oqChMmTLFdlcDabooWSCEkHpQUGbC62vO4dANra1NJRfjx6f7YEgHbydG1nju3LmDBQsWYMmSJdBqtbx+qVSK4cOHIy4uDuHh4a0meWoJKFkghJA6upVdhBdXnMFdXYmtrZOfGxY+2xchXi372rvZbMa2bdugVquxe/duh2NCQkIQGRmJmJgYeHp6NnKEpD5QskAIIXXw29UsvLH2PIoNldvuRnTzx5d/6glXact9i01LS8PChQuxaNEihzt4isViDBkyBHFxcejfvz+Ewta1VqOlabm/yYQQ0oAYY/jvwdv4Ys91u62b3xrXEXNGdYCgBVaMtFqt2LNnD9RqNRITEx1uNRwQEIDIyEhMnToVvr6+ToiSNARKFggh5DGVGs14d/1FbL+UYWtzlQjx1VO9MKGbvxMjaxhZWVlYunQpEhISkJyczOsXCoUYMGAAYmNjMWzYMJpFaIEoWSCEkMeQmluKl1aewdWMQltbqJccC5/ti45+LWdVP2MMBw8eREJCAjZu3AiTycQb4+Pjg4kTJyI+Ph6BgYFOiJI0FkoWCCGkhk7czsGrq88gr7Tyg3PYE974fkZvqOQSJ0ZWf3Jzc7F8+XIkJCTg+vXrvH6O49C3b19MnToVo0aNglgsdkKUpLFRskAIIY/AGMPKkxp8nGhfCOqFoW0xf2LnZl8IijGGkydPQq1WY926dQ4LOXl4eCAiIgJxcXEICwtr/CCJU1GyQAgh1TCYLfjH5stYl5Rqa5OIBPj31O6ICw9yYmR1V1hYiNWrV0OtVuPixYu8fo7j0KNHD8TExGDcuHFUyKkVo2SBEEIeIrtIj7+uOoszmspCUH7uUix4pi96BqucF1gdnT17FgkJCVi9ejVKSkp4/e7u7hg3bhzi4+PxxBNPOCFC0tRQskAIIQ5cTMvHSyvOILOwshBUnxAV1LPC4dsMC0GVlJRg3bp1UKvVOH36tMMx3bp1w5QpUzBx4kS4uLg0coSkKaNkgRBCHrDpXBrm/3rJrhDUU32D8b8x3SAVNa/bAi9fvoyEhASsWLECBQUFvH5XV1eMHj0a8fHx6Nq1K23BTByiZIEQQu6zWBk+3XUNCw7bF4L6R2RXPDsotNl8kBoMBvz6669Qq9U4cuSIwzFPPPEEoqOjERkZSYWcyCNRskAIIQAKSk2Ys+YsjtzU2do85GL8OLMPBrdvHoWgbt68iQULFmDp0qXIycnh9ctkMowYMQLx8fHo1atXs0l+iPNRskAIafVuZhXhxRVJSM4ptbV19i8vBBXs2bQLQZlMJmzduhVqtRr79u1zOCYsLAzR0dGIjo6GSqVq3ABJi0DJAiGkVdt7JQt/X3sOJUaLrW1Sd398Ma0n5JKm+xaZkpJiK+SUmZnJ65dIJBgyZAji4+PRr18/CATNey8I4lxN9y+BEEIakNXK8MOBW/hq7w279nfGd8Rrozo0ySl6i8WCXbt2Qa1WY8eOHQ4LObVp08ZWyMnbu3lcPiFNHyULhJBWp8RgxjvrL2DnH5XfyBVSEb5+qhfGdfVzYmSOZWRkYMmSJViwYAFSUlJ4/SKRCIMGDUJsbCwGDx5MhZxIvaNkgRDSqqTmluLFFUm4lllkawvzkmPR7L7o4Nt07gqwWq3Yv38/1Go1tmzZArPZzBvj6+uLyZMnIy4uDv7+La/aJWk6KFkghLQax2/p8OrPZ5FfpRDU8I4++H56byjlTaMgkk6nw7Jly5CQkIBbt27x+gUCAfr27Yu4uDgMHz6cCjmRRkHJAiGkxWOMYfnxZHyy/apdIaiXh7fD3IjOEAqcuz6BMYZjx45BrVZj/fr1MBqNvDFeXl6IiIhAfHw8goODnRAlac0oWSCEtGgGswX/b9MfWH8mzdYmFQnwaVwPxPRu48TIgPz8fKxatQpqtRqXL1/m9XMch969eyMmJgZjxoyBVCp1QpSEULJACGnBsgv1eHnVGZxLybe1BShlSHgmHD2CVE6JiTGGpKQkqNVqrFmzBmVlZbwxSqUS48ePR3x8PNq3b++EKAmxR8kCIaRFOp+aj5dXJiGr0GBr6xvqgf/O6gNft8YvBFVcXIw1a9ZArVbj7NmzDsd0794d0dHRmDhxIpWDJk0KJQuEkBZnw5k0vL/pEoxVCkHN6B+Mj6OfhETUuJsTXbx4EQkJCVi5ciWKiop4/QqFAmPHjkV8fDw6d+7cqLERUlOULBBCWgyzxYr/23ENS47dtbWJBBw+jO6GWQNCGm2jpbKyMmzYsAFqtRrHjx93OKZTp06YMmUKJk2aBIVC0ShxEVJblCwQQlqEvBIj/rbmHI7eqiwE5ekqwU8z+2BAO69GieH69etISEjAsmXLkJeXx+t3cXHBqFGjEB8fj+7duzfJXSIJcYSSBUJIs3c9s7wQVEpuZSGorgHuWPBsOII8GrYQlNFoxObNm6FWq3HgwAGHY9q1a4eoqChER0dDqVQ2aDyENARKFgghzdquPzLx1i/nUVqlEFRkjwB8Ht8TLpKG2/b47t27WLhwIRYvXozs7Gxev1QqxbBhwxAfH4/w8HCaRSDNGiULhJBmyWpl+G7/TXyz76atjeOAdyd0wl9HtG+QD2ez2YwdO3ZArVZj165dYIzxxgQFBSEqKgoxMTHw8mqcyx+ENDRKFgghzU6xwYy3fzmP3ZezbG1uUhG+ndELozvXfyGo9PR0LF68GAsXLkRaWhqvXywWY9CgQYiLi8PAgQOpkBNpcShZIIQ0Kyk55YWgrmdV3obYztsVC2f3RXuf+rurwGq1Yu/evVCr1UhMTITFYuGN8ff3x+TJkxEbGws/v6ZXrZKQ+kLJAiGk2Th2S4fXHigENaqTD76Z3htKl/opqJSdnY2lS5diwYIFuHPnDq9fKBSif//+iI2NxbBhwyAS0dsoafnot5wQ0uQxxrDkWDL+b4d9Iai/jmyPd8Z3qnMhKMYYDh8+DLVajV9//RUmk4k3xtvbGxMnTkRcXByCgoLq9HiENDeULBBCmjS9yYIPNv2BX89WrhWQicsLQU3pVbdCULm5uVixYgUSEhJw7do1Xj/HcQgPD0dMTAxGjx4NiURSp8cjpLmiZIEQ0mRlFerx0sozuJCab2sLVMqw4Nm+eLJN7fYrYIzh1KlTUKvVWLduHfR6PW+MSqXChAkTEB8fj7Zt29Y2fEJaDEoWCCFN0tmUPLyy8gyyiyoLQfUL88BPs8LhrXj8Us1FRUVYvXo11Go1Lly44HBMjx49EBMTg/Hjx1MhJ0KqoGSBENLk/JKUiv+36Q8YLZWFoGYOCMGHUd0euxDU+fPnoVarsXr1ahQXF/P63dzcbIWcOnXqVOfYCWmJKFkghDQZZosV/9pxFUuPJdvaRAIOH0V3w6yBoTU+T2lpKX755Reo1WqcOnXK4ZguXbrYCjnJ5Q27JTQhzR0lC4SQJiGvxIjXfj6L47dzbG1erhL8NCsc/dt61ugcV69eRUJCApYvX478/Hxev1wux+jRozFt2jR07dqVtmAmpIYoWSCEON21zEK8uCIJqblltrZuge5Y8GxftFG5VHuswWDApk2boFarcejQIYdjOnTogOjoaERFRcHNza1eYyekNaBkgRDiVDsvZeDt9RfsCkFF9wzEp3E9qi0EdefOHSxYsABLliyBVqvl9ctkMgwfPhxxcXHo06cPzSIQUgeULBBCnMJqZfhm3w18t/+WrY3jgPkRnfHS8HYOP9zNZjMSExOhVquxZ88eh+cNCQlBVFQUpkyZAk/Pml2+IIRUj5IFQkijK9Kb8Oa6C9h3tUohKJkI383ojVGdfHnjU1NTsWjRIixatAj37t3j9YvFYgwdOhSxsbHo378/FXIipJ5RskAIaVTJuhK8uCIJN7Mrb2Ns7+OKhc/2RbsqhaAsFgv27NkDtVqNbdu2wWq18s4VGBiIyZMnY+rUqfD15ScZhLR6zAJl1ilI9fxLdY+DkgVCSKM5fEOLOT+fRaHebGsb3dkX30zvBXdZeSGorKwsLFmyBAsWLEBycjLvHEKhEAMHDkRsbCyGDh1KswiEOCDW6+CVshM+mu2Q6LMhtJQ9+qBqULJACGlwjDEsOnIX/955FVXqQOG1Ue3x1rhOEHDA/v37oVarsWnTJpjNZt45fH19MXHiRMTHxyMgIKARoyekmWBWuOnOwic5EaqsY+AYfzautihZIIQ0KL3Jgvc2XsKmc+m2NplYgM/je2JwkBTffvM1EhIScOPGDd6xAoEAffv2xdSpUzFy5EiIxfVThpqQlkRkyIdX6i54a7ZBVspf08M4ISxCOQD+DqY1fow6xEcIIdXKKCjDyyvP4GJaga0tUCXDnJ5irP9yPqb98gsMBgPvOE9PT0RERCAuLg6hoTXfuZGQVoMxKHIvwVuTCI+MwxBY+WXVDS6+0IVEQhcyEYZlL8MpyQLHcXcArGeMzav1oxNCWqyk5Fy8suosdMWVyUCwVI+8Df/EzPdO8sZzHIeePXsiJiYGY8eOpUJOhDggNBXDM20PfJIT4VKs4fUzCFDgNwDa0CgU+vYDuPpZ01OXmYUwAD71EgUhpEVZ+3sK/mfLHzBZKhcolF3chaO7fwKsFrux7u7uGDduHKZNm4YOHTo0dqiENH2MQZ5/DT6aRHimH4DAyp+NM0q9oAuZBF3IJJjkfvUeAl2GIITUG5PFik+2XcGKE5XfeJjFhNw9P6H4ov0mSt26dcOUKVMwceJEuLhUv6UzIa2RwFwGz/Tf4JO8FfLCWw7HFHqHQxsWhXy/wYCg4T7SKVkghNSLnGIDZi84gj+yK7/1WIrzoN38fzCkXwUAuLq6YsyYMYiPj0eXLl1oC2ZCHHApuA1vTSK80vdBaC7l9ZvF7tCFTIQuNBIG1zaNEhMlC4SQOtHr9fh+1Wao/7DCIlPa2g0ZN6Dd9C9YinLQsWNHREdHY/LkyVTIiRAHOIsBHvcOwUezFYq8Kw7HFHn2gDY0CvkBw8CEkkaNr67JQi+O4/5RmwMZY/9bx8cmhDjRzZs3sWDBAqw8eBmS4X+BoMqCxOI/9qPk4CKMGzYY8fHx6NmzJ80iEOKAtCgFPppEeKXtgchUxOs3i1yREzwButBI6N3CGj/A++qaLPS8//M4OAAMACULhDQzJpMJW7ZsgVqtxm+/7Ydq2Cwox7xm62dWC7iLW/B0JwWmvLkRKpXKecES0kRxVhNUGUfho0mEW855h2NKVJ2hDY1CbuAoMJHz7wyqa7JwG8Cx+giEENJ0aTQaLFy4EIsXL0ZmZiY4iRw+cf8DeYf+tjECsx7xgYX405/+SrMIhDggKbkH75Tt8E7ZBbExj9dvEcqQ22YstKGRKFN1dEKED1fXZOEoY+zP9RIJIaRJsVgs2LlzJ9RqNXbs2AHGym+DFHkEwjfufyD2CraNDZADH4wIQoCClkERYsdqgTLrBHw0iXDXJoED4w0pdW9XPovQZiysYlcnBPlo9JdNCLGTkZGBxYsXY8GCBUhNTbXrk7XtA5/oeRDIKt/Q+gVK8Xp/JeRiQWOHSkiTJS7T3p9F2AGJXsfrtwrEyAscCW1oNEo8ugJNfDaOkgVCCKxWq62Q0+bNm2GxWHhjAsfMhjg83u5NLb6LK57qpoCgib/REdIomBXu2iT4JCdCmXUCHPiFnPSuwdCGRSEnaDwsEncnBFk7lCwQ0opptVosW7YMCxYswK1b/E1fhEIh+g4YBLcxr+CmUWVrlwo5/K2/EoOCnL/wihBnExly4ZWyCz4p2yEtzeD1WzkR8gOGQRsahWKvnk1+FsERShYIaWUYYzh69CjUajU2bNgAo9HIG+Pl5YWJEydidGQcVt51wc38ypLRvnIh5g1RIUxFFSBJK8YYFDnn4aNJhCrjKASMX1bd4OIPbWgkckIiYJZ6OiHI+lPrZIExRhcoCWlG8vPzsXLlSqjValy5wt/0heM49O7dG1OnTsXo0aNxu5DDFyfyUWCofBN80keCtwep4C6lP3/SOgmNhfBK3Q0fzTbISlJ5/eWFnAZBGxaNQp9wgGsZfys0s0BIC8YYQ1JSEtRqNdasWYOysjLeGKVSiQkTJiAuLg7t27cHAOy5XYpF5wpRpQ4UJnWQY3ZPN4gEzW8KlZA6YQyueVfgo0mEx72DEFj5s3FGmTd0IZPLCzm5tLwai5QsENICFRcXY82aNVCr1Th79qzDMd27d8eUKVMQERFhKwdtsjIsOVeIPXcqkwqRAHi5jztGt5U3SuyENBUCUwk80/fBR5MIeeEdXj8Dh0KfftCGRaHAdyAgqJ9y0E0RJQuEtCAXL15EQkICVq5ciaIi/taxCoUCY8eORXx8PDp37mzXV6C34PMT+biqM9naVDIB5g1WoaNX4+5DT4gzuRTchE9yIjzT90Fo0fP6TRJVeSGnkMkwugY6IcLGR8kCIc1cWVkZ1q9fD7VajRMnTjgc07lzZ0yZMgWTJk2Cqyt/05c7eSZ8eiwPurLKW72e8BTj3cEqeLm03G9LhFTgzHp43jsAH00iXPOvORxT5NWzvJCT/9BGL+TkbJQsENJMXb9+HQkJCVi2bBny8vhbx8rlcowaNQpxcXHo3r37Q7dgPppShh+TCmCssrXCyFAZXg5XQiKk9QmkZZMVJcNbsw1eqbshMpfw+s1iBXKCKgo5hTohwqaBkgVCmhGj0YjNmzdDrVbjwIEDDse0a9cO0dHRiIqKglKpdDgGACyMYc0fxdh0rfINUgBgdk83TH5CTvUdSIvFWYxQZRyBj2Yr3HIvORxT7NEV2tAo5AWOBBNKGznCpoeSBUKagbt379oKOWVnZ/P6pVIphg8fjri4OISHhz/yg77EZMU3JwtwNtNga1OIObw1SIWefvTGSFomaUn6/VmEXRAbC3j9FqELcoLGQRcahTJleydE2HRRskBIE2U2m7Fjxw6o1Wrs2rXLVsipquDgYERFRSEmJgaenjXb9CWt0IxPj+XhXnHldYdgdxHmD1HBnwpBkZbGaoYq6zh8khPhrjvjcEipe4fyQk5BY2AV0V0/jtA7AyFNTHp6OhYtWoSFCxciPT2d1y8WizF48GDExcVhwIABEAprvgDxTIYe35wsQKm5MvHof78QlAsVgiItiLg0C94pO8oLORlyeP1WgRS5bUZBGxqFUlXnZrkFc2OiZIGQJsBqtWLv3r1Qq9VITEx0WMgpICAAkydPxtSpU+Hn5/dY52eMYdP1Evx8qdiuQO5TXRWI7+pKhaBIy8AscM8+DR9NIpRZpxwWcipThN6fRRgHi8TNCUE2T5QsEOJE2dnZWLp0KRISEnD37l1ev1AoRP/+/REXF4dhw4Y91ixCBYOZ4cekAhxLrbxfXCbk8PoAJQa0oUJQpPkT6XPhnboT3pptkJZl8fqtAnFlISfPHjSLUAuULBDSyBhjOHToENRqNTZu3AiTycQb4+Pjg4kTJyIuLg5t2rSp9WNll1jw2fE83K1SCMrPtbwQVKiSCkGRZoxZ4aY7Dx/NVqgyj4Fj/Nk4vTwQutBI5ARHwCxVNX6MLQglC4Q0ktzcXKxYsQIJCQm4do2/6QvHcQgPD0dMTAxGjx4NiaRum75c1hrxxfE8FBorLzx095Xg7YEquFEhKNJMCY0FVQo5pfH6GSdAvt8QaMOiUOTdp8UUcnI2ShYIaUCMMZw8eRIJCQlYt24d9Hr+1rEeHh4YP3484uPj0bZt23p5zN13yrDkgUJQkU/I8WwPNwipEBRpbhiDa95l+CRvhUfGIQis/Nk4o8znfjnoiTDJvJ0QZMtGyQIhDaCwsBCrV6+GWq3GxYsXHY7p2bMnYmJiMG7cOFshp7oyWRkWnyvE3iqFoMQC4OVwJUaFudTLYxDSWASmYnillRdycinir+lh4FDo2w/a0GgU+A5o0YWcnI2SBULq0blz56BWq7F69WqUlPC3jnV3d8e4ceMQFxeHjh071utj5+st+Px4Pq7lVH7r8pAJMJcKQZFmRp5/Hd6aRHim73dcyEnqAV3wJOhCJ8Mo93dChK0PJQuE1FFpaSnWrVsHtVqN33//3eGYrl27YsqUKZg4cSLk8vrf9OVWrgmfHc9DzgOFoOYOVsGTCkGRZkBgLoNH+v7yQk4FNxyOKfTuA21oFAr8B4MJaIFuY6JkgZBaunLlChISErB8+XIUFPC3jpXL5Rg9ejSmTZuGrl27NlithcOaMvyUVABjlVvKR4e54MU+7lQIijR5ssI78NEkwittH4QOCzm5QxdcXsjJoAh2QoQEoGSBkMdiMBiwceNGqNVqHD582OGYJ554AlFRUYiKioKbW8Nt+mJhDKsuFmHrjVJbm4ADnu/phokdqBAUabo4ixEeGYfgk5wIRd4fDscUezwJbVgU8gJGtLpy0E0RJQuE1MDt27exYMECLFmyBDqdjtcvk8kwYsQIxMfHo1evXg3+QV1stOLrk/k4n2W0tSkkHN4ZpEJ3XyoERZomaXEqvDXb4J26GyJTIa/fInJFTtA4aEOjoHev+51BpP5QskDIQ5hMJiQmJkKtVmPv3r0Ox4SGhiIqKgpTpkyBh4dHo8SVWmjGf47lIbNKIagQpQjzBlMhKNL0cFYTlJnH4KNJhLvunMMxJcpO5bMIgaNgFdFdO00RvbMQ8oCUlBQsWrQIixYtQkZGBq9fIpFgyJAhiI+PR9++fWu1BXNtnb6nx7enClBWpRDUwDZSzOmvhIuINp8hTYekNBPemu3wTt0BsSGP128RypDbZjR0oVEoVXVyQoTkcVCyQAgAi8WC3bt3Q61WY/v27bBa+QVoAgMDERkZiZiYGPj6+jZqfIwx/HqtBGv/sC8ENb2bAnFdqBAUaSKYBcqsU+WzCNm/gwO/rHqZW9vyzZOCxsEqVjghSFIblCyQVi0zMxNLlizBggULoNFoeP0ikQgDBw5EbGwshgwZ0qizCBXKzFb8eLoAJ9IMtjaZiMMb/ZXoT4WgSBMg1uvglbITPprtkOizef1WgRh5ASOgDYtGiUc3KuTUDFGyQFodq9WKAwcOQK1WY/PmzTCbzbwxvr6+mDx5MmJjYxEQEOCEKMtllZjx6bF8aAoqY/RXCDFvsAohVAiKOBOzwk17Bj6aRKiyjoNj/Nk4vWsQtKFRyAkeD4tE6YQgSX2hZIG0Gjk5OVi2bBkSEhJw8+ZNXr9AIEDfvn0RGxuLESNGQCx27ofxH9kGfHEiH0VVCkH19JPgrYEqKCS0PoE4h8iQD6/UXfDRJEJayl/Twzgh8vyHQhcWhSKv3jSL0EJQskBaNMYYjh8/DrVajfXr18NgMPDGeHp6IiIiAvHx8QgJCXFClPYYY9h5qxRLLxTBWuWSb3RHOWZ1p0JQxAkYgyL3YvksQsYRh4WcDC5+0IVMhi5kEswyTycESRoSJQukRSooKMCqVaugVqvxxx/8TV84jkOvXr0QExODMWPG1Fshp7oyWRgWnC3E/mT7QlB/7avEiFC6pYw0LqGxCJ5pe8sLORXz1/QwCFDgNwDa0CgU+vYDONpavKWiZIG0KElJSVCr1VizZg1KS0t5/UqlEuPGjUN8fDw6dOjghAgfLq/Mgs9P5ON6lUJQni4CzBvsgQ6etD6BNBLGIM+/Bh/NVnimH4DAauQNMUq9oAstn0UwuTTunUHEOShZIM1eSUkJ1qxZA7VajTNnzjgc061bN8TExCAiIgIuLk3vG/rNXCM+O5aPXH3lIrFOXmK8O1gFDxl9WyMNT2AuhWfab/DRJEJeeMvhmAKfvtCFRiHfbxAgoI+P1oT+tUmzdenSJSQkJGDlypUoLORvHatQKDBmzBjEx8ejc+fOTbZWwsHkMqjPFMBUZTH5mLYueLG3O8RUCIo0MJeC2+WzCGn7ILSU8fpNEiVygidCFzoZBtc2ToiQNAWULJBmRa/XY8OGDVCr1Th27JjDMZ06dUJ0dDQmT54MhaLpbvpisTKsvFSExCqFoIQc8HwvN0S0p0JQpOFwFgM87h2EjyYRirwrDscUefaANiwK+f7DqJAToWSBNA83btxAQkICli1bhtzcXF6/TCbDqFGjEB8fjx49ejT5D9qi+4WgLlQpBOUu4fD2IBWepEJQpIHIijTw1iTCK20PRKZiXr9Z5Iqc++Wg9W5hjR8gabIoWSBNltFoxJYtW6BWq7F//36HY8LCwhAdHY3o6GioVKrGDbCWUgpM+M+xfGSVVBaCClOKMG+ICr6u9CdJ6hdnMUKVeRQ+mkS45VxwOKZE1Rna0CjkBo4CEzWNO4NI00LvTKTJSU5OxsKFC7F48WJkZWXx+iUSCYYNG4a4uDj069evyc8iVHUqXY/vfi+AvkohqEFBMszp5w4ZFYIi9UhScg8+KdvglbILYmM+r7+8kNNYaEMjUabq2PgBkmaFkgXSJFgsFuzYsQNqtRo7d+4EY/wCNEFBQZg8eTKmTp0Kb29vJ0RZe1bGsOFqCdZdrpz65QBMf1KBuM6uzSrhIU2Y1QJl1gn4aBKh1J52OKTUvV35LEKbsbCKXRs5QNJcUbJAnOrevXtYvHgxFi5ciNTUVF6/SCTC4MGDMXXqVAwePNgphZzqqsxsxfe/F+BUeuXukS4iDm8MUKJfIE35kroTl2nhnbId3inbIdHn8PqtAglyA0dBFxaFElUX2oKZPDZKFkijs1qt+O2336BWq7FlyxZYLBbeGD8/P1shJ39/fydEWT8yi8sLQaUUVhaCClAIMX+IB4Lc6c+P1AGzwl2bBJ/kRCizToCDo0JOwdCGRSEnaDwsEncnBElaCnq3Io1Gq9XaCjndvn2b1y8UCtGvXz9MnToVw4cPd3ohp7q6mGXAlyfzUVylEFRvfwneHKCCKxWCIrUkMuTCK2UXfDTbIC3L5PVbORHyA4ZBGxqFYq+eNItA6gUlC6RBMcZw9OhRqNVqbNiwAUYjf+tYLy8vTJw4EXFxcQgODnZClPWLMYbtt0qx/IFCUDGdXPF0dwWE9OZNHhdjUORcqCzkxPhl1Q0u/tCGRiInJAJmKRVyIvWLkgXSIPLz87Fy5Uqo1WpcucLf9IXjOPTp08dWyEkiaRmbvhjvF4I6UKUQlEQAvNpPiWEhTW+badK0CY2F8ErdDR/NNshK+Gt6GATI9x8EXWg0Cn3CAY5mrEjDoGSB1BvGGE6fPg21Wo21a9eirIy/daxKpcL48eMRHx+Pdu3aOSHKhpNbZsFnx/NxM7eyEJSXiwDzhnigvUfzvqRCGhFjcM27Ah9NIjzuHXRcyEnmbSsHbXLxcUKQpLWhZIHUWXFxMX7++Weo1WqcO3fO4Zju3btjypQpiIiIaDLloOvTjRwjPjuej7wqhaC6eIvxziAVVFQIitSAwFQCz/R99ws53eH1M3Ao9OkHbVgUCnwHAgL6vSKNh5IFUmsXLlxAQkICVq1ahaKiIl6/QqHA2LFjbYWcWqr9yaVIOFMIc5XF6OPbueDPvd0hFtD6BFI9l4Kb8ElOhGf6Pggtel6/SaKCLmQidCGTYXQNdEKEhFCyQB5TWVkZ1q9fD7VajRMnTjgc06VLF1shJ7lc3sgRNh6LlWH5xSJsv2lfCOovvd0xoX3Lfd6k7jizHp73DsBHkwjX/GsOxxR69YIuLBr5/kPABHQZizgXJQukRq5du4aEhAQsX74ceXl5vH65XI5Ro0YhLi4O3bt3b/E7EhYZrPjyZD4uZVcpBCUV4J1BKnTzaRmLNUn9kxXdhY9mGzxT90BkLuH1m8VuyAkaD21oFAxuIU6IkBDHKFkgD2U0GrFp0yao1WocPHjQ4Zj27dsjKioKU6ZMgZubW+MG6CSa+4WgsqsUgmqrEmHeEA/4yOk6MrHHWYxQZRyBj2Yr3HIvORxT7NEV2tBo5AWOABNS1VHS9FCyQHju3LmDhQsXYsmSJcjOzub1S6VSDB8+HHFxcQgPD2/xswhVnUzT4/vfC6C3VG6gMCRYhtf6KiEVtZ7XgTyatDgN3inb4JW6G2JjAa/fIpIjp81Y6EKjUKZs74QICak5ShYIAMBsNmPbtm1Qq9XYs2ePw0JOwcHBiIqKQkxMDDw9W9emL1bG8MuVYqy/Ujl1zAGY2V2BmE5UCIrcZzVDlXkcPppEuOvOOBxS6t4B2rBo5LYZA6uI9t4gzQMlC61cWloaFi1ahEWLFiE9PZ3XLxaLMXjwYMTGxmLgwIHNspBTXZWZrPj29wKcvldZCEou5vD3ASqEB9CUMQHEpVnwSdkO75QdEBtyef1WgRS5bUZBGxaNUmUn2oKZNDuULLRCVqsVe/bsgVqtRmJiIqxWfgGagIAAREZGYurUqfD19XVClE1Dxv1CUKlVCkEFupUXgmrjRn8+rRqzwD37NHw0W6HM+t1hIacyRSi0YVHIDRoPi1jhhCAJqR/0bteKZGVlYenSpViwYAHu3r3L6xcKhRgwYABiY2MxbNiwVjmLUNX5TAO+OpmPElPlJZnwACneGKCEq5i21W2tRPoceKfshHfKNkjL+Gt6rAIx8gKGQxcahWLP7jSLQFoEShZaOMYYDh06BLVajY0bN8JkMvHG+Pj4YOLEiYiPj0dgIG36whjDtpulWHGhyO67YmxnV0x/kgpBtUrMCjfdufJCTpnHwDF+WXW9PBC6sCjogiJgkSqdECQhDYeShRYqNzcXK1asgFqtxvXr13n9HMehb9++mDp1KkaNGtXsy0HXF4OFIeFMAQ5pKnfSkwiBOf2UGBJMi9FaG6GhAN5pu+Ct2QZZCX9ND+MEyPcfCm1oJIq8+1AhJ9JiUbLQgjDGcPLkSajVavzyyy/Q6/lbx3p4eCAiIgJxcXEICwtr/CCbsJwyCz47lo9beZWzL95yAeYN9kA7KgTVejAG19w/ygs5ZRyCwMqfjTPKfKENnYyckIkwybydECQhjYuShRagsLAQq1evhlqtxsWLF3n9HMehR48eiImJwbhx41pkIae6uqYz4vMT+civUgiqq7cY7wz2gFJK3xZbA4GpGF5pe+GjSYRLUTKvn4FDge8A6MKiUODbH+Ba95oe0rpQstCMnTt3Dmq1GqtXr0ZJCX/rWHd3d4wbNw7x8fF44oknnBBh87DvbikWnrUvBBXRXo7ne7lBRIWgWjx5/vXyWYT0/Y4LOUk9oQuZBF3IJBjl/k6IkBDno2ShmSktLcW6deugVqvx+++/OxzTrVs3TJkyBRMnToSLC11nfxizlWHZhSLsvFVZCErEAS/0cce4dlQIqiUTmMvgmf4bvDWJcC246XBMoXc4tKGRyPcfAgjorZK0bvQX0ExcvnwZCQkJWLFiBQoK+FvHurq6YvTo0YiPj0fXrl1pR8FHKDRY8cWJfFzWVhaCUkoFeHewCl28qRBUS+VSeBvemm3wStsLobmU128Wu0MXEgFdSCQMiiAnREhI00TJQhNmMBjw66+/Qq1W48iRIw7HPPHEE4iOjkZkZGSrKeRUV8n5JvznWB60pZXXHdp7iDB3sAe8qRBUi8NZDPC4dwg+mkQo8i47HFPk2R260CjkBQwHE1KySMiDKFlogm7evIkFCxZg2bJl0Ol0vH6ZTIYRI0YgPj4evXr1olmEx3A8VY8fThfAUKUQ1LAQGf7aVwmpkF7HlkRanAIfTXkhJ5GpiNdvEbkiJ3g8tKGR0Lu1dUKEhDQflCw0ESaTCVu3boVarca+ffscjgkLC0NkZCRiYmKgUqkaN8BmzsoY1l4uxq9XKxeCCgDM6uGG6I5ySrhaCM5qgirzGLyTt8I957zDMSWqTvfLQY+kQk6E1BAlC06WkpKChQsXYtGiRcjMzOT1SyQSDBkyBPHx8ejXrx8EArqN73GVmqz49lQBkjLsC0G9NVCF3v5UCKolkJRmwFuzHd4pOyE25vH6LUIZctuMgS40CqWqjk6IkJDmjZIFJ7BYLNi1axfUajV27NjhsJBTYGCgbRahNRdyqqt7RWZ8eiwPaUWV2/O2uV8IKpAKQTVvVguU2SfLy0FnnwYHfln1Mre20IZFI6fNGFipkBMhtUbvlo0oIyMDS5YswYIFC5CSksLrFwqFGDhwIOLi4jBkyJBWX8iprs5lGvD1A4Wg+t4vBCWnQlDNlrhMC++UHfBO2Q6Jnr+mp7yQ0whow6JR4tGNCjkRUg8oWWhgVqsVBw4cgFqtxubNm2E2m3ljfH19MWnSJMTFxSEgIMAJUbYsjDFsuVGK1RftC0HFd3HFU90UENCHR/PDrHDXnoG3ZitUWSfAMf5snN41CNrQKOQEj4dFQoWcCKlPlCw0EJ1Oh2XLliEhIQG3bt3i9QsEAvTt2xexsbEYMWIEFXKqJwYLw09JBTiSUrkTn1TIYU4/JQYH0zbXzY3IkAev1F3w0WyDtDSD1884IfL8h0IXFoUir940i0BIA6FkoR4xxnD8+HGo1WqsX78eBoOBN8bT0xMRERGYNm0agoODnRBly6UrteCz43m4nVc5e+MjF2D+EA+EqSgZazYYgyLnInw0W6HKOAIB48/GGVz8oAuNhC54IswyTycESUjrQslCPSgoKMDKlSuhVqtx+TJ/0xeO49CrVy/ExMRgzJgxVMipAVzTGfHZ8XwUGCqnp7v5SPDOIBXcqRBUsyA0FsErbQ+8NYlwKeav6WEQoMBvALSh0Sj07UuFnAhpRJQs1EFSUhLUajXWrFmD0lL+1rFKpdJWyKlDhw5OiLB12HunFIvOFsJcZTH8xA5yPNeTCkE1eYxBnn8NPpqt8Ew/AIHVyBtilHpBFzoZupBJMLnQnUGEOAMlC4+ppKQEa9asgVqtxpkzZxyOefLJJzFlyhRERERQIacGZLYyLD1fhF237QtBvRjujrFtqRBUUyYwl8Iz7Tf4aBIhL+Sv6QGAAp++0IVGId9vEBVyIsTJ6C+whi5duoSEhASsXLkShYWFvH6FQoExY8YgPj4eXbp0cUKErUuBwYovT+ThstZka1PdLwTVmQpBNVkuBbfLZxHS9kFoKeP1myRK5ARPhDZ0MoyubZwQISHEEUoWqqHX67Fhwwao1WocO3bM4ZhOnTohKioKkZGRUCho05fGcPd+ISjdA4Wg5g32gBcVgmpyOLMenhkH4a3ZBkXeFYdjijy7QxsWjXz/YVTIiZAmiJIFB27cuIEFCxZg6dKlyM3N5fW7uLhg5MiRmDZtGrp37051BRrRsdQy/HC6AMbKDRkxPESGV6gQVJMjLUqBjyYRXmm7ITIV8/rNYgVygsZDFxoJvVtY4wdICKkxShbuMxqN2LJlC9RqNfbv3+9wTLt27RAVFYXo6GgolbTpS2OyMIa1fxRj4zX7QlDP9nRD5BNUCKqp4CxGqDKPwkeTCLecCw7HFKu6QBcWhdyAkWAiujOIkOag1ScLycnJWLhwIRYvXoysrCxev1QqxbBhwxAXF4e+ffvSh5ITlNwvBHWmSiEohZjDW4NU6OlHhaCaAknJPfikbINXyi6Ijfm8fovQBblBY6ENjUSZ8onGD5AQUietMlmwWCzYsWMH1Go1du7cCcb4BWiCgoIQGRmJqVOnwsvLywlREgBILzLjP8fycK9KIaggdxHmD1EhQNEqf32bDqsFyqwT8NEkQqk97XBIqXt7aEMjkRs0DlYR3aFCSHPVqt5t7927h8WLF2PhwoVITU3l9YvFYgwaNAixsbEYNGgQFXJysjMZBnxzMh+lVTZQ6BcoxRv9lXChQlBOIy7LrlLIKYfXbxVIkNtmFLSh0ShVdaYtmAlpAVp8smC1WrFv3z4kJCRgy5YtsFgsvDH+/v6YPHkyYmNj4efn54QoSVWMMWy+XoLVl4rtig5P6+qKP3WlQlBOwSxwz04qn0XIOgkO/EJOZYoQ6EKjkBM0HhaJmxOCJIQ0lBabLGi1WixduhQJCQm4c+cOr18oFKJfv36IjY3F8OHDIRK12JeiWTGYGX5MKsCx1MpCUDIhh7/1V2JgEC2Ga2wiQy68U3bBW7MN0rJMXr+VEyE/YBi0odEo9upBswiEtFAt6hOSMYYjR45ArVbj119/hdHI3zrW29sbERERiI+PR1BQkBOiJA+jLbXg02N5uJtfWTjI11WI+UNUCFVSIahGwxgUOefho0mER8YRcIw/G2eQB0AbMhk5IRNhlno4IUhCSGNqEclCXl6erZDT1atXef0cxyE8PBwxMTEYPXo0JBLa9KWpuaw14osT+SisUgiqu68Ebw9UwY0KQTUKobEAXql74KNJhKwkjdfPOAHy/QZDFxqFQp9wgKN/F0Jai2abLDDG8Pvvv0OtVmPt2rXQ6/W8MSqVChMmTEB8fDzatm3rhChJTey+XYrF5wphqbJAIfIJOZ7t4QYhFYJqWIzBNe9y+SzCvYMQWE28IUaZT3khp+CJMLn4OCFIQoizNbtkoaioCD///DPUajXOnz/vcEz37t0RExODCRMmUDnoJsxkZVhyrhB77lTWCBAJgJfD3TE6jG6za0gCUwm80vfBOzkR8iL+mh4GDoW+/aANjUKB70BAQHcGEdKaNatkISMjA23atEFRURGvz83NDWPHjkV8fDw6derkhOjI48jXW/DFiXxc1VV+k/WQCTB3sAodvegyUUOR59+AtyYRnum/QWjhz8aZJB7QhUyELnQyjPIAJ0RICGmKmlWykJeXx2vr0qULoqOjMXnyZMjl9G20ObiTZ8Knx/KgK6tcn/CEpxhzB6vg6ULfYOubwFwGj3sH4JOcCNeC6w7HFHr3Li8H7T8ETECLSQkh9ppVslBBLpdj1KhRmDZtGrp160ZbMDcjR1LK8N/TBTBWuU1/VJgLXurjDgkVgqpXssK79ws57YXQXMLrN4vdkBMcAW3oZBgUIU6IkBDSXDSrZEEikeC5557DjBkz4OZGm740JxbG8POlYmy+XqUQFAc819MNkzpQIaj6wlmM8Mg4DG9NItxyLzkcU+zRDdrQKOQFjgATUm0NQsijNatkITAwEKNHj6ZEoZkpNlrxzal8nMus3PdCIeHw9kAVelAhqHohLU6Dd8o2eKfsgshUyOu3iOTICRoHXWgkytzbOyFCQkhz1qySBdL8pBWWF4LKKK7c2CfEXYR5Q1Twp0JQdWM1Q5V5HD6aRLjrzjgcUureAdqwaOS2GQOryKWRAySEtBT0bk0azOl7enx7qgBlVQpBDWgjxd/6K+Eiog19aktcmgWflO3wTtkBsSGX128VSJHbZjS0YVEoVXaiLZgJIXVGyQKpd4wxbLxWgjV/2BeCmt5NgbgurlQIqjaYBcrs3+GtSYQy6/eHFHIKhTYsCrlB42ERK5wQJCGkpaJkgdQrvdmKH04X4kRalUJQIg6v91diQBvaIOtxifQ5tnLQ0rJsXr9VIEZewAhoQ6NQ4vkkzSIQQhoEJQuk3mSXmPHpsXwkF1QWgvK7XwgqhApB1Ryzwk13Dj6aRKgyjzks5KR3bQNdaCR0QRGwSJVOCJIQ0ppQskDqxR/ZBnxxIh9FxsoLDz18JXhrkApuElqfUBMiQz68UnfBO2U7ZCXpvH7GCZDvPxTa0CgUefemQk6EkEZDyQKpE8YYdt0uxZLzRbBSIajHxxgUuZfgrUmER8bhhxRy8oU2tLwctEnm7YQgCSGtHSULpNZMFoaF5wrx293KQlBiAfBKuBIjw+g2veoITcXwTNsDn+REuBRreP0MHAr8BkAbGo1C334AR9tgE0Kch5IFUit5egs+P56P6zmV34Q9ZQLMHaLCE55UCMohxiAvuA6f5ER4pu+HwGrgDTFJPaELmQRdyCQY5f5OCJIQQvgoWSCP7VauCZ8ez0NulUJQHe8XgvKgQlA8AnMZPNN/g7cmEa4FNx2OKfTuA21oNPL9BwMC+rMkhDQt9K5EHsshTRl+SiqAqcpt/qPvF4ISUyEoOy6Ft+GdnAiv9H0Qmkt5/WaxO3QhEdCFRMKgCHJChIQQUjOULJAasTCGVReLsPVG5YeegAOe7+WGie2pEFQFzmKAx71D8NFshSLvisMxRZ7dy2cRAoaBCemSDSGk6aNkgTxSsdGKr07m40JWZSEoNwmHtwep0N2XCkEBgLQ4BT6abfBK3Q2RqYjXbxa5Ijd4PLShkdC7tXVChIQQUnuULJBqpd4vBJVZpRBUqLK8EJSfa+v+9eGsJqgyj8E7eSvcc847HFOi6gRtaDTyAkdSISdCSLPVut/tSbVO39Pjm1MF0FcpBDUoSIo5/ZSQteJCUJLSDHhrtsM7ZSfExjxev0UoQ26bMdCGRqFM1dEJERJCSP2iZIHwMMbw69USrLlcbNc+434hqFa5PsFqgTL7ZHk56OzT4OxKZJUrdWtXXsipzVhYxa5OCJIQQhoGJQvETpnZih9+L8DJ9Mo9AFxEHN4YoES/wNZXCEpcprUVcpLodbx+q0CMvMBR5YWcPLpSISdCSItEyQKxySox4z/H8pFSpRCUv0KI94Z4IMi9Ff2qMCvctWfgrdkKVdYJcIxfDlrvGgRtaBRygsfDIqFCToSQlq0VfQKQ6lzKNuDLBwpB9fKT4M2BKihaSSEokSEPXqm74KPZBmlpBq+fcULkBQyDLjQSRV69aRaBENJqULLQyjHGsPNWKZZesC8ENaWTK2Z2V0DY0j8QGYMi50J5OeiMIxAwM2+IwcUf2tBI5IREwCz1dEKQhBDiXJQstGImC8OCs4XYn1xZCEoiAP7aV4nhoS37Nj+hsRBeaXvgrdkGl+IUXj+DAAV+g6ANi0KhT18qB00IadUoWWil8sos+Ox4Pm7kVhaC8nIRYO5gD3TwFDsxsgbEGFzzr8I7ORGe9w5AYDXyhhhlXtCFTIYuZDJMLj5OCJIQQpoeShZaoZu5Rnx2LB+5+sqFe529xHhnsAoespZXCEpgLoVn2j74aBIhL7ztcEyBTz9oQ6NQ4DcIELS814AQQuqCkoVW5mByGdRn7AtBjW3rghf6uEMsaFnrE1wKbsFHkwjPtH0QWsp4/SaJCjkhEdCGRMLoGuiECAkhpHmgZKGVsFgZVlwswrablYWghBzw597umNDOpcVstMSZ9fDMOAjv5EQo8q86HFPk1Qva0Ejk+w+lQk6EEFIDlCy0AkUGK748mY9L2ZXX6N0lHN4Z7IFuPi3jw1JWpIG3JhFeaXsgMhXz+s1iBXKCJkAXGgm9W6gTIiSEkOaLkoUWLqXAhP8cy0dWSWUhqLYqEeYO9oCva/O+Ns9ZjFBlHoGPZhvcci44HFPs0RW60EjkBowEE7W+HSgJIaQ+ULLQgp1K1+O7UwXQWyo3UBgcJMNr/dybdSEoSUk6fDTb4ZW6E2JjAa/fInRBbtDY8kJOyg5OiJAQQloWShZaICtj2HClBOuuVE7HcwCeflKBqZ2baSEoqxmqrBPw1iRCqU1yOKTUvT20odHIDRoDq0jeyAES0jSUlZXh119/xZEjR3D37l0UFhbCxcUFYWFhGDBgAGJiYuDv7+/sMJudhIQELFy4EB9++CGioqKcHU6jo2ShhSkzWfHd7wX4/V5lISi5iMPfByoRHtD8puHFZdnl5aBTd0Ciz+H1WwUS5LYZBV1oFEpUXWgLZtKqXbhwAXPnzkVOTg5kMhm6d+8OT09PFBcX48qVK7h06RJWrFiBr7/+GgMGDGi0uF566SWcPXsWW7duRWBg07zzqDnE6EyULLQgmcXlhaBSCyu3LA5UCDGvuRWCYha4Z5+GjyYRyqxT4MAv5FSmCIEuNAo5QeNhkbg5IUhCmpbr16/j1VdfhcFgwOzZs/HCCy/AxaVyJ1ar1YqDBw/iu+++Q3Z2thMjbZ6eeuopTJgwAd7e3s4OxSma0ScIqc6FLAO+OpGPYlPl+oTe/hK8OUAF12ZSCEqkz4V36k54a7ZBWpbF67dyIuQHDoc2NArFnj1oFoGQ+xhj+Mc//gGDwYCXXnoJL730Em+MQCDA6NGj0b9/f2RmZjohyuZNpVJBpVI5OwynoWShmWOMYdvNUqy4UGT3/XtqZ1fMeLIZFIJiDG455+CdnAiPzKPgmIU3xCAPuF8OegLMUg8nBElI03b8+HHcvn0bfn5++POf/1ztWIVCgQ4dKhf+6vV6rFq1Cnv27EF6ejpEIhE6duyI+Ph4TJgwgXd8VFQUMjIykJSUhM2bN2Pt2rVITU2FQqHA8OHD8frrr8PNrXy27969e4iOjrYdW/V/A0BSUvn6o48++gjbtm2DWq2GyWTCsmXLcP36dRQXF+PAgQNwc3PDuXPnsGfPHpw7dw5ZWVkwGo3w9/fHyJEj8dxzz9ke80F3797FypUrcfr0aeh0OigUCgQHB2PUqFGYMWMGsrOzaxRjdWsW8vPzsXz5chw8eBBZWVmQyWTo1q0bZs6ciYEDB/Ji6tu3LwICArB582asXLkSW7duRUZGBjw9PTFhwgS88sorkEia1m3tlCw0Y0YLQ8KZAhzU6G1tEiHwal8lhoU07UJQQmMBvFJ3w0ezDbKSNF4/4wTI9xsEbWg0inzCqZATIdU4evQoAGDMmDEQiWr+tl5SUoJXXnkFV69ehYeHB4YOHYqysjIkJSXh3LlzuHTpEt555x2Hx3777bdYu3YtwsPDERwcjAsXLmDTpk1ITk7GggULwHEc5HI5IiMjceLECeTk5GD06NGQyx+++Hj37t3YvHkzunTpgsGDByMtLc22IPvbb7/FzZs30aFDB/Tr1w9GoxHXrl3D8uXLceTIESxbtox37n379uEf//gHjEYj2rZti5EjR6KkpAS3b9/Gt99+i5iYmMeO8UHZ2dl48cUXkZ6eDn9/f4wYMQL5+fn4/fffceLECbz55puYOXOmw2P/3//7fzh27BjCw8MRGhqKc+fOYcWKFdBqtfjkk09qHENjoGShmcops+Dz4/m4WaUQlLeLAPOGeKCdRxMtBMUYXPMuwyd5KzwyDkFgNfGGGGU+98tBT4RJ1jqvDRLyuG7cuAEA6Ny582Md9+OPP+Lq1avo27cvvvzyS7i6ugIAkpOT8dJLL2Ht2rUYMGAAhg0bxjt2x44dWLNmDcLCwgCUf7t+/vnnce7cOSQlJaFfv35QqVT46KOP8NJLLyEnJwd///vfq108uGnTJvzf//0fxo8fz+t78cUX0bNnTygUClub0WjE559/jk2bNmH16tV48cUXbX0pKSn48MMPYbFY8M9//hMRERG2PsYYTp06BalUColE8lgxPujf//430tPTERERgQ8//BBicfn77/nz5zFnzhx899136Nu3Lzp16mR3XEZGBmQyGTZu3GhbB5Geno5Zs2Zh586dePnllxEUFFTjOBoafV1rhm7kGDFvX45dotDFW4zPxno1yURBYCqGz93N6HroBXQ+9jq80vfZJQoMHAp8B+BWv3/i0pifkdnxGUoUCHkMBQXl+414eNT8Ml1ZWRm2bNkCgUCA+fPn2xIFAAgLC7NdzlizZo3D4//617/aEgWg/Jp+XFwcAODs2bOP+xQAAEOHDnWYKADAkCFD7BIFAJBIJHj77bchFApx6NAhu76ff/4ZBoMBMTExdokCAHAch4EDB9Z5qj8tLQ1HjhyBXC7Hu+++a0sUAKBXr16Ii4uDxWLB+vXrHR7/zjvv2C2YbNOmDSZOnAgAOHfuXJ1iq280s9DM7L9bioSzhTBXWaAwvp0L/ty76RWCkuffgLdmKzzT90No0fP6TVIP6IInQRc6GUY53fdNSGO6evUqDAYDunbtavehX2Hy5Mn44osvcOHCBVitVggE9t8tHV2LDwkJAQDodLpaxTR8+PBq+7Ozs3H48GFoNBoUFxeDsfIF3WKxGKmpqXZjf//9dwBAbGxsrWKpifPnzwMABg0aBKVSyeufNGkSVq9e7fCDXyQSoW/fvrz2ur6GDYWShWbCbGVYfqEIO27ZF4L6S293TGjfdDYgEpjL4JG+Hz6abXAtuO5wTKF3b2hDo1HgPxhM0PRmQghpbio+qPLy8mp8jFarBQAEBAQ47Hdzc4NCoUBxcTEKCwt5dwL4+vryjqm41m8y8S8x1kR1m0WtWrUKP/zwA8xm80PHVJWVVX5HVUNO5Vd8oD/sskVFe8VrXZWXlxeEQv6W+xUzPEajkdfnTJQsNAMOC0FJBXh3kApdm0ghKFnhHfhoEuGVtg9Ccwmv3yx2hy64vJCTQRHshAgJabk6duyICxcu4Nq1a5g0aVK9nbe63V4fnGmoDw+7LHDp0iV88803UCgUeOeddxAeHg4vLy/b+IiIiCb3TRxo/NevIVGy0MRp7heCyq5SCKqdSoS5QzzgI3duISjOYoRHxiH4JCdCkfeHwzHFHk9CGxaFvIARVA6akAYydOhQrF+/Hr/99htef/31Gt0R4ePjA6B8oZ0jxcXFKCoqglQqhbu7e73G+7gOHDgAAHj11VcRGRlp16fX65GTw9/d1c/PDykpKUhLS+MtLqwvFesNHvYa3rt3D0Dla92cNa/UppU5kabH+7/l2iUKQ4Nl+OcoL6cmCtLiVLS5/BN67P0T2p77Ny9RsIjkyA6bgssjFuH60O+QGzSOEgVCGtDgwYPRrl07ZGVlYcmSJdWOLS4uxu3bt9GlSxdIpVJcu3YNKSkpvHE7duwAAPTs2bPO34IrFv5ZLPx9VGqiqKgIgONLH/v27bOtXaiqf//+AMrvsGioGHv16gUAOHHihC3Gqnbu3AkA6N27d43P2VRRstAEWRnDmj+K8MWJfFvFSA7AM90V+PsAJaSixl/IyFlNUN07iCdOvIMnD8yG/531EJkK7caUKJ+ApsdbuDhuPVK7vwG9e7tGj5OQ1ojjOHzyySeQSqVYsGABfvjhB5SVldmNYYzh0KFDePbZZ3HlyhW4uLggOjoaVqsVn376qd14jUaDxYsXAwCmT59e5/gqvllrNJpaHV+x6G/Lli12axbu3LmD77//3uExM2bMgFQqxaZNm7Bnzx67PsYYTp48abcuoDYxBgUFYejQoSgpKcEXX3xhF9vFixexYcMGCIVCTJs2rcbnbKroMkQTU3q/ENTpqoWgxBzeHKBCnwBpo8cjKc20FXISG/iLp6wCKXKDxkAbGoVSVcNM9RFCHq1Tp0748ccfMW/ePCxbtgxr165Fjx49bIWkrl69ipycHEilUvj5+QEA5syZgz/++AOnTp3ClClT0KdPH9umTAaDAdOnT3/kHQo1MXz4cGzbtg0ffPABBg4caLsF8n/+539qdHxUVBRWrVqFI0eOIC4uDl27dkVBQQHOnj2LkSNH4vLly7xLAaGhofjHP/6BDz/8EO+//z4WLVqEDh062GZWsrKycODAAdu6h9rG+P777+PFF1/E9u3bcfbsWXTv3h35+fk4c+YMLBYL/v73vzfYZZDGRMlCE5JxvxBUWtVCUG5CvDfEA4FujfhPxSxQZp2CjyYR7tm/gwN/iq/MLQza0CjkBo2DRaxwcBJCSGPr1asXNm3ahI0bN+Lw4cO4efMmCgsLIZfLERoaitjYWMTExNiSBVdXVyxYsACrVq3C3r17cfjwYYjFYnTp0gXx8fG8/Qlqa/To0XjrrbewefNmHDlyxPaNvqbJgkqlwooVK/D999/j7NmzOHz4MAIDA/HKK6/gmWeeQUxMjMPjJkyYgHbt2mHFihU4c+YM9u/fD3d3dwQHB2PGjBl2hbZqG6Ovry+WL1+OZcuW4eDBgzhw4ABkMhn69euHWbNmObzFtDniHF3raYo4jrscFhbW9T//+Y/dvuYtxflMA74+aV8IKjxAijcGKOEqbpyrRWK9Dl4pO+Gj2Q6Jnl+VzioQIy9gBLRh0Sjx6EaFnAghpJmIj49HcnJyGmOsVrej0cyCkzHGkHizFCsfKAQV29kV0xujEBSzwk13Fj7JiVBlHQPH+OWg9a5tygs5BU2ARcrfeIQQQkjLRsmCExksDOqkAhxOsS8ENaefEkOCG7YQlMiQD6/UXfDWbIOs9B6vn3FC5PsPKS/k5N2LCjkRQkgrRsmCk+SUWvDp8Tzczqtcn+AjLy8E1VbVQLsaMgZF7iV4axLhkXHYYSEng4svdCGR0IVMhFnm1TBxEEIIaVYoWXCCazojPj+ej3xD5ZR/Nx8x3h7kAaW0/r/BC03F8EzdAx9NIlyK+bcFMXAo8BsAXWgUCnz7A5xzN3sihBDStFCy0Mj23SnFwrOFMFdZVxrRXo7ne7lBVJ+FoBiDPP8afDSJ8Ew/AIHVwBtilHohJ2QidCGTqJATIYSQh6JkoZGYrQxLzxdh1+3KQlAiDnixjzvGtqu/QlACcxk803+DT/JWyAtvORxT6B0ObVgU8v0GAwL6FSCEEFI9+qRoBAUGK748kY/L2srdwlRSAd4drEJn7/rZBtml8Da8kxPhlb4PQnMpr98sdocuZGJ5ISfXNvXymIQQQloHShYa2N18Ez49lgdtaeX6hPYeIswb7AGvOtZ34CwGeNw7BB/NVijyrjgcU+TZHdrQaOQHDKP6DIQQQmqFkoUGdCy1DD+eLoTBUrlAYXiIDK/0VUIqrP36BGlxCnw02+CVuhsiE794iVnkitzg8dCGRkLv1rbWj0MIIYQAlCw0iPJCUMXYeK3E1iYA8EwPN0R1lFdb4/xhOKsJqoyj8NYkwj3nvMMxJarO5VswB44CE8lqGT0hhBBij5KFelZisuLbUwU4k1F594GrmMNbA1Xo5f/4haAkpRnw1myDd8ouiI38Qk4WoQy5bcZCGxqJMlXHOsVOCCGEOELJQj26V2TGf47lIb2osh56kLsI84eoEKB4jJfaaoEy+yR8khPhrj3tsJBTqXu78lmENmNhFbvWR/iEEEKIQ5Qs1JNzmQZ8dTIfpVUKQfULlOKN/kq41LAQlLhMC++UHfBO2Q6JXsfrtwokyA0cBV1YFEpUXaiQEyGEkEZByUIdMcaw5XoJVl0qtvv+H9/FFU91U0DwqA90ZoW7Ngk+yYlQZp0AB0eFnIKhDYtCTtB4WCTu9fsECCGEkEegZKEODGaGn84U4EiVQlAyIYc5/ZUYFFT9AkORIRdeKbvgk7Id0tIMXr+VEyE/YBi0oVEo9upJswiEEEKchpKFWtKWWvDZsTzcya8sBOUrF2LeEBXCHlYIijEoci7AR5MIVcYRCJiZN8Tg4g9taCRyQiJglno2VPiEEEJIjVGyUAtX7xeCKqhSCOpJHwneHqSCu4NCUEJjIbzS9sAnORGyklReP4MABX6DoA2LQqFPXyoHTQghpEmhZOEx7bldisXn7AtBTeogx+yeDxSCYgyueVfgo0mEx72DEFiNvHMZZd7QhUyGLmQSTC4+jRA9IYQQ8vgoWaghk5VhyblC7LlTZmsTCYCX+7hjdNvKQlACUwk80/fBR5MIeeEd3nkYOBT69IM2LAoFvgMBAZWDJoQQ0rRRslADBXoLPj+Rj6s6k61NJRNg3mAVOnqV11twKbgJn+REeKbvg9Ci553DJFEhJyQC2pBIGF0DGy12QgghpK4oWXiEO3nlhaB0ZZXrE57wFOPdwSp4i03wTNkJH00iXPOvOTy+yKsXtKFRyA8YCiZ4yMJHQgghpAmjZKEaR1LK8N+kAhgrN2TEyFAZ/t4xD4G3f4Zn6h6IzCW848xiN+QET4A2JBIGt5BGjJgQQgipf5QsOGBhDD9fKsbm65WJgAwm/KvtRYw37IHbkUsOjyv26AptaDTyAkeACR+/DgQhhBDSFFGy8IASoxVfn8rHuczyuxdCuUzMFu/HDMlhuGQU8sZbhC7IDRoLbWg0ypTtGztcQgghpMFRslBFWqEZnx7LQ3axARMEZzFT+BuGC+/PIjywf1Kpewdow6KR22Y0rCI5/2SEEEJIC0HJwn1nMvRYd/I2puM3TJcehC+XzxtjFUiR22YUtKFRKFV1pi2YCSGEtAqtPllgVjMunzuKdqnbsEd4DkKOXw66TBFaXg46aBwsEjcnREkIIYQ4D8cY/8OxKeI4rlAsFrv5+flBIpHU/YRWC4SmIgiMRRDC4mAAB4tIDovYDVahlGYRCCGENFupqakwmUx6xphLbY5vTslCJgA3ALkA+HsnE0IIIeRh/AEYGWMetTm42SQLhBBCCHEOKm9ICCGEkGpRskAIIYSQalGyQAghhJBqUbJACCGEkGpRskAIIYSQalGyQAghhJBqUbJACCGEkGpRskAIIYSQalGyQAipM47j/DmOW8pxXAbHcWaO4wo4jtvEcVyos2MjhNQd7eBICKkTjuPGAfgFgBzAIQB5AEYA8ANwB0APxliJ8yIkhNQVzSwQQmqN47hwAIkAdACeZIyNZ4w9BaAjgKsA2gGY4cQQCSH1gJIFQkitcBwnAPATykvdT2eM3azoY4wVAlh5//8OdkJ4hJB6RMkCIaS2hgLoB2A9Y+yMg/6K5EHeeCERQhoCJQuEkNqKuv/fzQ/p97z/38yGD4UQ0pAoWSCE1Nao+/899pD+Pvf/e64RYiGENCC6G4IQ8tg4jhMDKAYgAeDCGNM/0C8AkAbAB0AwY4xmFwhpxmhmgRBSG0+iPFEAgGAH/W8CCACwgRIFQpo/ShYIIbXRp8r/ns9xHFfxfziOiwbwCcr3W5jX2IERQuqfyNkBEEKapYpk4UcArwHox3HcHwA6oPwOiXwA0YyxFOeERwipTzSzQAipjYpk4V8AngMgBBCL8ksSiwH0ZIwddU5ohJD6RgscCSGPheM4IYBCAIWMsQBnx0MIaXg0s0AIeVydUb7R0llnB0IIaRyULBBCHlfFJQhKFghpJShZIIQ8LkoWCGllaM0CIYQQQqpFMwuEEEIIqRYlC4QQQgipFiULhBBCCKkWJQuEEEIIqRYlC4QQQgipFiULhBBCCKkWJQuEEEIIqRYlC4QQQgipFiULhBBCCKkWJQuEEEIIqRYlC4QQQgipFiULhBBCCKkWJQuEEEIIqRYlC4QQQgipFiULhBBCCKkWJQuEEEIIqdb/B7LTt3rVxIbOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 600x450 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rho = np.logspace(-1,1,50)\n",
    "T = rho**(1/3)\n",
    "T_mona = rho**(5/3-1)\n",
    "T_poly = 0.5*rho**(4/3-1)\n",
    "\n",
    "fig,ax = plt.subplots(1,1,figsize=[4,3],dpi=150)\n",
    "ax.set_ylim(0.4,1.6)\n",
    "ax.axes.get_xaxis().set_ticks([])\n",
    "ax.set_xlim(-1,1)\n",
    "ax.axes.get_yaxis().set_ticks([])\n",
    "ax.set_xlabel(r\"$\\rho$\")\n",
    "ax.set_ylabel(\"T\")\n",
    "    \n",
    "ax.plot(np.log10(rho),1+np.log10(T),'k',label='Jeans Density')\n",
    "ax.fill_between(np.log10(rho), 0, 1+np.log10(T), color='k', alpha=.2)\n",
    "ax.plot(np.log10(rho),1+np.log10(T_mona),label='Monatonic gas')\n",
    "ax.plot(np.log10(rho),1+np.log10(T_poly),label='Real gas')\n",
    "ax.text(0.4,0.5,'Contraction')\n",
    "ax.text(-0.9,1.3,'Expansion')\n",
    "plt.legend(loc=1)\n",
    "plt.savefig(\"Figures/Temp_density.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67e7ca11-94c1-426d-95ac-57c0495797a0",
   "metadata": {},
   "source": [
    "## A collapsing protostar\n",
    "So let's assume we are at the final stage of this collapse - that is, we have a 1 solar mass dust cloud of molecular hydrogen at a temperature of 20 K, and an average density of:\n",
    "$$\n",
    "    \\rho = \\rho_{\\rm J} = \\frac{3}{4 \\pi (1 {\\rm M_{\\odot}})^2}\\left[ \\frac{3k (20 {\\rm K})}{2 G (2\\times m_{\\rm H})}\\right] ^3\n",
    "$$\n",
    "The radius of such a cloud is $\\sim 10^{15}$m. This cloud will collapse unopposed as long as the gravitational energy released is not converted into thermal motion of the particles. This will happen as long as the energy goes into dissassociation of hydrogen molecules and the ionization of H atoms.\n",
    "\n",
    "## The radius at which pressure becomes important\n",
    "So how small is the gas cloud when this happens? In order to estimate it, we need to know how much energy can go into dissociating and ionising H.\n",
    "\n",
    "The energy to dissassociated a H molecules is $\\epsilon _D=4.5$eV, while the energy needed to ionise a H atom is $\\epsilon _I=13.6$eV. As such, the energy needed to dissociate all of the H molecules, and ionise of the H atoms, is approxmiately\n",
    "$$\n",
    "    \\frac{M}{2m_H}\\epsilon _D+\\frac{M}{m_H}\\epsilon _I\n",
    "$$\n",
    "Assuming that this energy comes from the graviational collapse of a cloud from $R_1$ to $R_2$, we can then equate these\n",
    "$$\n",
    "    \\frac{GM^2}{r_2}-\\frac{GM^2}{r_1}=\\frac{M}{2m_H}\\epsilon _D+\\frac{M}{m_H}\\epsilon _I\n",
    "$$\n",
    "If we plug in $10^{15}$m as the initial radius, and assume a 1 solar mass cloud, we find that $R_2\\sim 10^{11}$m (or 0.6 AU), and that the free fall time is $\\sim$ 100,000 years.\n",
    "\n",
    "## The core temperature once freefall has been halted\n",
    "So at this stage, the gravitational energy released by the collapse can no longer go into dissociating or ionising hydrogen. From this point onwards, any further contraction equates to an increase in temperature of the cloud, which (assuming an ideal gas) means the pressure increases, slowing the collapse. Our next step is to calculate the temperature at this stage. \n",
    "\n",
    "We are in hydrostatic equilibrium, so we can use the virial theorem:\n",
    "$$\n",
    "    <P>V = -\\frac{1}{3} E_{\\rm GP} = NkT\n",
    "$$\n",
    "to estimate the temperature. Now, $E_{\\rm GP} \\sim -\\frac{GM^2}{r_2}$. As such, T is $\\sim 60,000$ K.\n",
    "\n",
    "This is below the temeperature required to drive fusion $(\\sim 10^{7})$ K, so the cloud must collapse further to become a star.\n",
    "\n",
    "Below is a code snippet showing where all of these numbers come from."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "19e6a43f-1bdc-4160-ac53-c453b5be4195",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The critical density is 3.854896e-16 kg / m3\n",
      "The initial radius is 1.540682e+06 solRad/1.071852e+15 m\n",
      "The final radius is 125.635490 solRad/8.740461e+10 m\n",
      "The freefall time is 1.072130e+05 yr\n",
      "The internal temperature assuming hydrostatic equilibrium at R2 is 61315.537409 K\n"
     ]
    }
   ],
   "source": [
    "M = 1*u.solMass\n",
    "m = 2*c.m_p\n",
    "T = 20*u.K\n",
    "ed = 4.5*u.eV\n",
    "ei = 13.6*u.eV\n",
    "\n",
    "rho = (3/(4*np.pi*M**2)) * ((3*c.k_B*T)/(2*c.G*m))**3\n",
    "r1 = (3*M/(4*np.pi*rho))**(1/3)\n",
    "r2 = c.G*M**2/((M/(2*c.m_p)*ed+M/(c.m_p)*ei+c.G*M**2/r1))\n",
    "teff = ((3*np.pi)/(32*c.G*rho))**(1/2)\n",
    "\n",
    "E_gp = -c.G*M**2/r2\n",
    "T = -E_gp/(3*M/c.m_p*c.k_B)\n",
    "\n",
    "print(\"The critical density is {:e}\".format(rho.to(u.kg/u.m**3)))\n",
    "print(\"The initial radius is {:e}/{:e}\".format(r1.to(u.solRad),r1.to(u.m)))\n",
    "print(\"The final radius is {:f}/{:e}\".format(r2.to(u.solRad),r2.to(u.m)))\n",
    "print(\"The freefall time is {:e}\".format(teff.to(u.year)))\n",
    "print(\"The internal temperature assuming hydrostatic equilibrium at R2 is {:f}\".format(T.to(u.K)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43733f01-3984-4fc1-80f8-5a075eb21d7c",
   "metadata": {},
   "source": [
    "## Degeneracy Pressure\n",
    "\n",
    "So, will all dust clouds collapse enough to create a high enough temperature to ignite fusion. The short answer is no, but to understand why we must start considering quantum mechanical effects.\n",
    "\n",
    "We will start by assuming that there is a minimum total kinetic energy for a collection of electrons, and that this degeneracy occurs when the average distance between electrons becomes comparable to the De Broglie wavelength, $\\lambda _{\\rm DB} = h/p$ where $h$ is Planck's contant and $p$ is the momentum. The kinetic energy of a classical electron is $E_{\\rm KE} = kT$, so the momentum is given by $p\\sim \\sqrt{2 m_{\\rm e} kT}$ (please excuse the mixing of classical and quantum here, but we just want to get a rough handle on the picture). As such, the De Broglie wavelength is\n",
    "$$\n",
    "    \\lambda_{\\rm DB} = \\frac{h}{\\sqrt{2 m_{\\rm e} kT}}\n",
    "$$\n",
    "So, a classical treatment of  the particles is appropriate when the distance between the electrons is much larger than this wavelength. Assume we have a cloud of electrons with a density of $\\rho = \\frac{\\bar{m}}{\\Delta x^3}$ (that is, each particle occupies a volume of $\\Delta x^3$). For a classical treatment, we need $\\Delta x >> \\lambda$. This is satisfied when:\n",
    "$$\n",
    "    \\rho << \\frac{\\bar{m}}{\\lambda ^3} = \\bar{m} \\frac{(2 m_{\\rm e} kT)^{3/2}}{h^3}\n",
    "$$\n",
    "\n",
    "Now, using our expression for the gravitational potential energy ($E_{GP}$) and kinetic energy $E_{KE}$, and assuming hydrostatic equilibrium is achieved and maintained, then we can show that the temperature of the star increases as\n",
    "$$\n",
    "    kT=\\frac{GM\\bar{m}}{3R} \\simeq G\\bar{m}M^{2/3}\\rho^{1/3}.\n",
    "$$\n",
    "As such, the temperature at the core will rise with $\\rho^{1/3}$ as long as the inequality on $\\rho$ holds.\n",
    "\n",
    "So, what happens when we reach the critical density? By comining the above two expressions, we arrive at an expression which tells us what the temperature of the core will be when this happens:\n",
    "$$\n",
    "    kT \\sim \\left(\\frac{G}{h}\\right)^2 \\bar{m}^{8/3} m_e M^{4/3}\n",
    "$$\n",
    "Around this temperature, the electrons become generate, and exert degeneracy pressure, halting collapse. So, if we want fusion to occur, we need the temperature at which this occurs to be greater than $1\\times 10^7$ K. The above approximateion predicts this should happen somewhere between 0.04 and 0.1 ${\\rm M_{\\odot}}$. A more detailed treatment shows this occurs when $M>0.06 {\\rm M_{\\odot}}$. So, for a dust cloud with an initial mass higher than this, the core temperature will get high enough to start fusion. \n",
    "\n",
    "Smaller than this, and electron degeneracy presssure will halt collapse before the temperature has gotten high enough to ignite. Such stars are called brown dwarfs, and produce very low luminosities, making them very hard to observer."
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}