{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "dd8ff36d", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pylab as plt" ] }, { "cell_type": "markdown", "id": "f661644f", "metadata": {}, "source": [ "Let's consider an luminosity profile that declines exponentially with radius with a scale length of h" ] }, { "cell_type": "code", "execution_count": 2, "id": "8509ca97", "metadata": {}, "outputs": [], "source": [ "h=4000 # scale length [pc]\n", "I0=100 # central luminosity density [Lsun/pc^2]\n", "R=np.arange(0,20000,10) # range of R: 0 to 20 kpc in steps of 10 pc\n", "I=I0*np.exp(-R/h) # luminosity profile" ] }, { "cell_type": "markdown", "id": "e288f360", "metadata": {}, "source": [ "Now let's say we want to know the radius at which I=30 Lsun/pc^2. In this case, there *is* an analytic solution so let's work that out." ] }, { "cell_type": "code", "execution_count": 3, "id": "84f5918c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The luminosity density hits a value of 30.00 Lsun/pc^2 at R = 4815.891 pc\n" ] } ], "source": [ "Iwant = 30 # Lsun/pc^2\n", "Rwant=-h*np.log(Iwant/I0) # note that we are using a natural log (np.log) since we are just solving an exponential.\n", "print('The luminosity density hits a value of {:.2f} Lsun/pc^2 at R = {:.3f} pc'.format(Iwant,Rwant))" ] }, { "cell_type": "markdown", "id": "52767493", "metadata": {}, "source": [ "But what if there *isn't* an analytic solution to the function? Below are two methods to work out the solution computationally. They work best when the function is monotonic.\n", "\n", "First lets plot the function to make sure it looks right and is monotonic. " ] }, { "cell_type": "code", "execution_count": 4, "id": "dc7a761e", "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEYCAYAAABRMYxdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+TklEQVR4nO3dd3yV9d3/8dcnJ4uEMEIChA2yZAqE6ag46qxoQQUX1oGo3VPb2vbX3va2te19a6ssxb3B1br3aFgBmTIFFGTvTUj4/P44B+9IEwghOdcZ7+fjcR7nXNe5zrneF4GLT67rO8zdEREREUkkKUEHEBEREalpKnBEREQk4ajAERERkYSjAkdEREQSjgocERERSTgqcERERCTh1FqBY2YTzWyDmc0vt+5uM1tkZnPN7AUza1DuvdvNbJmZLTazc2orl4iIiCS+2ryC8zBw7mHr3gK6uXsPYAlwO4CZdQGGA10jn7nfzEK1mE1EREQSWK0VOO7+IbDlsHVvuntpZHEq0CLyegjwtLvvd/cVwDKgX21lExERkcSWGuC+rwOeibxuTrjgOWR1ZN1/MLNRwCiA7OzsPp07d67NjCJJZebMmZvcPT/oHLEsLy/P27RpE3QMkYRRW+edQAocM/sVUAo8cWhVBZtVOIeEu48HxgMUFhZ6cXFxrWQUSUZm9nnQGWJdmzZt0HlHpObU1nkn6gWOmY0ELgTO9P+bCGs10LLcZi2ANdHOJiIiIokhqt3Ezexc4BfARe6+p9xbLwPDzSzDzNoCHYDp0cwmIiIiiaPWruCY2VPA6UCema0Gfku411QG8JaZAUx199HuvsDMngU+JXzr6lZ3L6utbCIiIpLYaq3AcfcRFax+8Ajb3wncWVt5REREJHloJGMRSThm1tLM3jOzhWa2wMx+UME2Zmb3RgYYnWtmvYPIKiK1I8hu4iIitaUU+Im7zzKzHGCmmb3l7p+W2+Y8wu39OgD9gTGRZxFJALqCIyIJx93XuvusyOudwEL+c2ytIcCjHjYVaGBmBUf77oNe4QgWIhJjVOCISEIzszZAL2DaYW81B1aVWz7iAKNmVmxmxUvWbmfjzv21klVEao4KHBFJWGZWF5gM/NDddxz+dgUfqXSAUXcvdPfCMlIYNraILzbvqWhTEYkRKnBEJCGZWRrh4uYJd3++gk2qNcBo27xstu89wLfHFLFgzfaaCSsiNU4FjogkHAsPtPUgsNDd/1bJZi8D10R6Uw0Atrv72qN9d1Z6iEmjB5IeMi4fN5WizzbVYHIRqSkqcEQkEZ0MXA2cYWazI4/zzWy0mY2ObPMqsBxYBkwAbqnql7dvnMPkWwbRrEEm106cwavzjloXiUiUqZu4iCQcd/+YitvYlN/GgVuru4+C+nV47qZBXP/IDG59cha/H9KNqwe0ru7XiUgN0xUcEZFqqp+VxmPX9+fMzo2548X5/O2tJbi6kYvEBBU4IiLHoU56iLFX9eGywhbc+85SfvnCfMoOqsgRCZpuUYmIHKfUUAp/GtqD/JwM7nvvM7bs3s89w3uRmRYKOppI0tIVHBGRGmBm/Oyczvz2W114Y8F6rpk4ne17DwQdSyRpqcAREalB3zm5LfeO6MUnX2zl8nFTWL9jX9CRRJKSChwRkRp2Uc9mPHRtP1Zt2cPQMUUs37gr6EgiSUcFjohILTilQx5PjRrA3pIyho2dwpxV24KOJJJUVOCIiNSSHi0aMOnmQWSlhxgxYSofLtkYdCSRpKECR0SkFrXNy+b5mwfRulE21z08g5dmfxl0JJGkoAJHRKSWNa6XyTM3DaBP64b84OnZTPx4RdCRRBKeChwRkSiol5nGI9f149yuTfn9vz7lT68v0qjHIrVIBY6ISJRkpoW478reXNm/FWPe/4yfTZpLadnBoGOJJCSNZCwiEkWhFOO/Lu5Gfk4G//v2UrbuLuEfV/SmTrpGPRapSbqCIyISZWbGD8/qyH9d3I13F2/gqgensW1PSdCxRBKKChwRkYBcNaA191/Rm3mrt3Pp2Cms2bY36EgiCUMFjohIgM7rXsAj1/Vj3fZ9DB1TxLINO4OOJJIQVOCIiARs4AmNeOamgZQedIaNncLMz7cGHUkk7qnAERGJAV2a1WPy6EE0qJPGlQ9M5d1F64OOJBLXVOCIiMSIVo2ymHTzIDo0zuHGR2cyaebqoCOJxC0VOCIiMSSvbgZPjRrAwHaN+Olzcxj7wWcaEFCkGlTgiIjEmLoZqUy8ti/f6tmMu15bxJ2vLOTgQRU5Isei1gocM5toZhvMbH65dblm9paZLY08Nyz33u1mtszMFpvZObWVS0QkHqSnpnDP5Sdx7aA2PPDxCn787GxKSjXqsUhV1eYVnIeBcw9bdxvwjrt3AN6JLGNmXYDhQNfIZ+43Mw3rKSJJLSXF+O23uvDzczvx4uw13PBoMbv3lwYdSyQu1FqB4+4fAlsOWz0EeCTy+hHg4nLrn3b3/e6+AlgG9DvaPpau36XfaEQkoZkZt5zenj8P7cHHSzdyxYSpbN61P+hYIjEv2m1wmrj7WoDIc+PI+ubAqnLbrY6s+w9mNsrMis2seF9pGQvWbK/VwCIiseCyvi0Zd3Uhi9bt5NKxU1i1ZU/QkURiWqw0MrYK1lXYos7dx7t7obsXAhoQS0SSxtldmvDEDf3ZtGs/Q8cUsXDtjqAjicSsaBc4682sACDyvCGyfjXQstx2LYA1R/uy9FAKxStV4IhI8ihsk8ukmweRYsZl46YwfcXhLQFEBKJf4LwMjIy8Hgm8VG79cDPLMLO2QAdg+tG+LCsjxMwvtmqMCBFJKh2b5DD5lkHk52Rw1YPTeGPBuqAjicSc2uwm/hQwBehkZqvN7HrgLuBsM1sKnB1Zxt0XAM8CnwKvA7e6e9nR9pGdnsrGnftZtUUz8IpIcmneoA6TRg+iS0E9bn58Jk9N/yLoSCIxJbW2vtjdR1Ty1pmVbH8ncOex7CMrPcQuoPjzLbRqlHWMCUVE4ltudjpP3tifW56Yxe3Pz2PTzv1894z2mFXUrFEkucRKI+NqyUwLkZORqobGIpK0stJTmXBNId/u1Zy/vrWE3768gDKNeixSe1dwoqVX64YqcEQkqaWFUvjLpT3Jz8lg3IfL2byrhL9d3pOMVI2XKskrrq/gABS2bsji9TvZvvdA0FFERAKTkmLcfv6J/Or8E3ll3lq+89AMdu7TeVGSV9wXOH1aN8QdPvlCV3FERG48rR3/c3lPpq/YwvDxU9m4U6MeS3KK+wLnpJYNCKUYs3SbSkQEgEt6tWDCyEKWb9zNsLFFfL55d9CRRKIu7guc7IxUTizIoVgFjojIVwZ3asyTN/Znx94DDB1TxPwvNa2NJJe4L3AAClvnMnvVNkrLNPGmiISZ2UQz22Bm8yt5/3Qz225msyOP30Q7Y23r1aohz40eREZqiOHjp1K0bFPQkUSiJiEKnN6tG7KnpIyFa3cGHUVEYsfDwLlH2eYjdz8p8vh9FDJFXfvGdZl88yCaN6jDtQ/N4JW5a4OOJBIVCVHg9GuTC8C0FZsDTiIiscLdPwQ0URPQtH4mz940kJ4t6/Pdp2bx6JSVQUcSqXUJUeA0rZ9J60ZZTNOkcyJybAaa2Rwze83Mula2kZmNMrNiMyveuHFjNPPVmPpZaTx2fX/O7NyE37y0gL++uVjz+ElCS4gCB2BA20ZMX7GFgxrBU0SqZhbQ2t17An8HXqxsQ3cf7+6F7l6Yn58frXw1LjMtxNirenN5YUv+/u4yfvnCPLVdlISVMAVO/3a5bN97gEXr1A5HRI7O3Xe4+67I61eBNDPLCzhWrUsNpXDX0O58d3B7npq+iluemMW+A0ed21gk7iRQgdMIUDscEakaM2tqkVkpzawf4fNhUpxAzIyfntOJ332rC28tXM81D07XaPCScBKmwGneoA4tc+swdXlSnJ9E5CjM7ClgCtDJzFab2fVmNtrMRkc2GQbMN7M5wL3AcE+yRinXntyWe4f34pNVW7l83BTW79gXdCSRGhP3k22W179tI95ZuJ6DB52UFAs6jogEyN1HHOX9fwD/iFKcmPWtns1omJXOTY8V8+37i3j0+n6ckF836Fgixy1hruAA9G+by9Y9B1i6YVfQUURE4sYpHfJ4etRA9peWMWxMEbNXbQs6kshxS6gCZ0CkHY5uU4mIHJvuLeozafQgcjLTuGLCVD5YEp/d4UUOSagCp2VuFs0b1FFDYxGRamiTl82kmwfSulE21z88gxc/+TLoSCLVllAFDoRvU01bvkUDWImIVEPjnEyeuWkAhW0a8sNnZvPAR8uDjiRSLYlX4LTLZfPuEpapHY6ISLXUy0zj4e/04/zuTfmvVxby368t1C+NEncSrsBROxwRkeOXmRbi7yN6c9WAVoz7YDk/fW4uBzTqscSRhCtwWkXa4RR9pgJHROR4hFKMPwzpxo/O6sjkWau56bGZ7C3RqMcSHxKuwDEzTm7fiKLPNlOmealERI6LmfGDszpw5yXdeH/xBq54YCpbd5cEHUvkqBKuwAE4uX0e2/ceYP6X24OOIiKSEK7s35r7r+zNgjU7uHTcFNZs2xt0JJEjStgCB+DjZZsCTiIikjjO7VbAo9f1Y/32fQwdU8TS9ZrcWGJXQhY4eXUzOLGgHv9WgSMiUqMGtGvEMzcNpPSgM2zsFGZ+viXoSCIVSsgCB+CU9o0oXrlVDeJERGpYl2b1eP7mQeRmp3PlA9N4Z+H6oCOJ/IeELXBObp9HSdlBivXbhYhIjWuZm8Wk0QPp2CSHUY/N5LniVUFHEvmahJpNvLx+bXNJCxkfL93EqR3yg44jIuWY2Y+rsNludx9X62Gk2hrVzeDJGwdw8+Mz+dmkuWzaVcLob7TDzIKOJhLMFRwz+5GZLTCz+Wb2lJllmlmumb1lZksjzw2PZx9Z6an0btVQDY1FYtPPgLpAzhEePwksnVRZ3YxUHhzZl4t6NuNPry/iD/9ayEEN0SExIOpXcMysOfB9oIu77zWzZ4HhQBfgHXe/y8xuA24DfnE8+zq1Qx5/eXMJW3aXkJudftzZRaTGPObuvz/SBmaWHa0wcnzSU1P438tPolHddCb+ewWbd+/n7mE9SU9N2FYQEgeC+tuXCtQxs1QgC1gDDAEeibz/CHDx8e7kUHfxos90FUcklrj7z2tiG4kdKSnGby7sws/P7cRLs9dw/SMz2LW/NOhYksSiXuC4+5fAX4AvgLXAdnd/E2ji7msj26wFGh/vvro3r09OZiofL1WBIxJrzKyzmZ1pZnUPW39uUJnk+JgZt5zenj8P60HRZ5u5YsJUNu/aH3QsSVJRL3AibWuGAG2BZkC2mV11DJ8fZWbFZla8cePGI26bGkrhlPZ5fLBko2bCFYkhZvZ94CXge8B8MxtS7u0/BpNKasplhS0Zf3UflqzfybCxU1i1ZU/QkSQJBXGL6ixghbtvdPcDwPPAIGC9mRUARJ43VPRhdx/v7oXuXpiff/TeUad3ymft9n0s1oibIrHkRqCPu18MnA7cYWY/iLynLjgJ4MwTm/DEDf3ZsruEoWOKWLh2R9CRJMkEUeB8AQwwsywL9yU8E1gIvAyMjGwzkvBvd8ft9E7hO13vLTry1R4RiaqQu+8CcPeVhIuc88zsb6jASRh9Wufy3OiBhFKMy8ZNYeryzUFHkiQSRBucacAkYBYwL5JhPHAXcLaZLQXOjiwftyb1MjmxoB7vL67wgpCIBGOdmZ10aCFS7FwI5AHdgwolNa9jkxwm3zyIJvUyuWbidF6fvy7oSJIkAulF5e6/dffO7t7N3a929/3uvtndz3T3DpHnGhuCeHCnfIo/38qOfQdq6itF5PhcA3ztfzp3L3X3a4DTgokktaVZgzo8d9NAujarxy1PzOTJaV8EHUmSQFIMUjC4c2PKDrp6U4nECHdf7e4V/irv7v+Odh6pfQ2z03nihv6c3qkxv3xhHve8vVSdP6RWJUWB06tlA+plpuo2lUiMMbPfBZ1BoicrPZVxV/dhaO8W/M/bS/jNSwso06jHUksqHck4keaKSQ2lcGrHfN5fHO4urnlSRIJlZinABCrpLSmJKy2Uwl8u7UFeTjrjPljO5t37+Z/LTyIjNRR0NEkwR7qCk1BzxQzu1JgNO/fzqboqisSCfwJb3P32oINI9JkZt593Ir++4ERenbeOayfOYKfaSEoNO9JcVAk1V8w3OobHzHl/8Ua6NqsfcBqRpFcI3Bl0CAnWDae2I69uBj99bg6Xj5vKw9f1pXFOZtCxJEFUegUn0eaKyc/JoHvz+ry7SFfERWLAYGCcmfUPOogE6+JezXnw2r6s3LybYWOmsHLT7qAjSYI4YiPjRJsrZnDnxnzyxVbNjSISMHf/FDgHuDvoLBK8b3TM58kbB7Bz3wGGjS1i/pfbg44kCaDSAicR54r5ZpcmHHR4R1dxRALn7muAC4LOIbHhpJYNmHTzIDJSQ1w+bgr/XqZhPeT4HOkKTsLNFdO1WT2a1c/krU/XBx1FRAB31yRx8pUT8uvy/C2DaJmbxbUPTedfc9cEHUni2JEKnISbK8bMOLtLEz5aupG9JWVBxxFJemb2lpk1KLfc0MzeCDCSBKxJvUyeuWkgvVo25HtPfcIjRSuDjiRx6kgFTkLOFXN2l6bsO3CQj5Zq8k2RGJDn7tsOLbj7VqBxcHEkFtSvk8aj1/fjrBOb8NuXF/CXNxZr1GM5ZkcqcBJyrpj+7XLJyUzVbSqR2HDQzFodWjCz1oD+JxMy00KMubI3I/q15B/vLeO2yfMoLTsYdCyJI5WOg+Puq4/wXtzOFZMWSmFwp8a8s2gDZQedUEpc3m0TSRS/Aj42sw8iy6cBo2rii81sIuGrzhvcvVsF7xtwD3A+sAe41t1n1cS+pWakhlL44yXdya+bwb3vLmPLnhL+PqIXmWka9ViO7qhzUSXiXDHf7NqELbtLmPn51qCjiCQ1d38d6A08E3n0cfeaaoPzMHCkIS3OAzpEHqOAMTW0X6lBZsaPv9mJ/3dRV95euJ6rH5zG9j0a9ViO7kjdxFPM7EEgI4p5ouIbHfNJCxlvfVrhZMYiEl2nAWcQHvzv1Jr6Unf/ENhyhE2GAI962FSggZkV1NT+pWaNHNSGv4/oxZxV27ls3BTWbd8XdCSJcUe6gpOwc8XkZKYx8IQ83vp0vRquiQTIzO4HRgPzgPnATWZ2X5R23xxYVW55dWTdfzCzUWZWbGbFGzeqg0JQLuzRjIe/05cvt+1l6Jgilm3YFXQkiWFHKnAKgReiFSTavtmlCSs372HROg3DIRKgbwDnuPtD7v4Q4fYwp0dp3xU1wKvwNx53H+/uhe5emJ+fX8ux5EgGtc/j6VED2F9axqVji5i9alvQkSRGHanASei5Ys7t1pQUg1fnrQ06ikgyWwy0KrfcEpgbpX2vjuzvkBaARpaLA92a12fyzYPIyUxjxPipvL9Yo9PLfzrSZJsJPVdMXt0MBrRrxCvz1uo2lUhwGgELzex9M3sf+BTIN7OXzezlWt73y8A1FjYA2O7u+o0nTrRulM3kmwfRLj+bGx4p5oVPKu34K0mq0m7iEJ4rxswSdq6Y87sX8OsX57No3U5OLKgXdByRZPSb2vpiM3uK8O2uPDNbDfwWSANw97HAq4RviS0j3E38O7WVRWpHfk4GT48awE2PzeRHz8xh864Sbji1XdCxJEYcscCBxJ4r5txuTfnNS/N5dd5aFTgiURSZjuF14DV3X1Qb+3D3EUd534Fba2PfEj05mWk89J2+/OiZ2fzXKwvZuHM/t53XmfAwR5LMqjIOTsLOFfPVbaq5uk0lEmUjga3A78xslpmNMbMhZlY36GASfzJSQ/x9RG+uGdiacR8u5yfPzeGARj1OekctcEjwuWLO717A8k271ZtKJIrcfZ27P+zuwwn32HwU6AO8YWZvm9nPg00o8SaUYvy/i7ryk7M78vysLxn1aDF7SkqDjiUBqkqBk9Bzxag3lUiw3P2gu09x99+4+8nAcODLoHNJ/DEzvndmB/772935YMlGrnxgGlt3lwQdSwJy1DY41OJcMbGg/G2qH5/dUfdtRaLIzPKBG4E2lDsfuft1QWWS+DeiXytys9P53lOfMGxsEY9e35/mDeoEHUui7KhXcGp5rpiYcEGP8G2qhWt1m0okyl4C6gNvA6+Ue4gcl3O6NuWx6/qxYed+ht5fxJL1Or8nm6rcooJamismVpzXrYDUFOOl2boqLhJlWe7+C3d/1t0nH3oEHUoSQ/92jXj2poEcdGfYmCKKVx5pajJJNFXpRRXkXDFRkZudzjc65vPS7DWUHUyY5kUi8eBfZnZ+0CEkcZ1YUI/JNw8ir24GVz4wjbc/XR90JImSqlzBCXKumKi5uFdz1u3Yx7Tlm4OOIpJMfkC4yNlrZjvMbKeZ7Qg6lCSWlrlZPDd6IJ2b5nDT4zN5dsaqo39I4l5VCpwg54qJmrO7NKFuRiovfKLbVCLR4u457p7i7nXcvV5kWaNuSo1rVDeDJ28cwKATGvHzyXO5771lGv8swVWlwKnxuWLMrIGZTTKzRWa20MwGmlluZFDBpZHnhtX57urKTAtxXremvDZ/HfsOlEVz1yJJy8xOq+gRdC5JTNkZqTw4si9DTmrG3W8s5vf/+pSDapaQsKrSTbw25oq5B3jd3YeZWTqQBfwSeMfd7zKz24DbgF/Uwr4rdUmv5jw3czVvL1zPhT2aRXPXIsnqZ+VeZwL9gJmEOzWI1Lj01BT+57KTyKubwYMfr2DzrhL+cmlP0lOr2udG4kWlBU5tzRVjZvUI98q6FsDdS4ASMxvC/7XteQR4nygXOP3bNaJpvUxemPWlChyRKHD3b5VfNrOWwJ8DiiNJIiXF+PUFJ5Kfk8Fdry1iy+4Sxl7dh7oZVfmdX+LFkUrW2porph2wEXjIzD4xswfMLBto4u5rASLPFU4HYWajzKzYzIo3btx4nFG+LpRiDDmpGR8s2cjmXftr9LtFpEpWA92CDiGJz8wY/Y0TuHtYD6Ys38wVE6aySef9hFJpgVOLc8WkEh44cIy79wJ2E74dVSXuPt7dC929MD8/v5oRKndxr+aUHnT+OWdNjX+3iHydmf3dzO6NPP4BfATMCTqXJI9LC1sy4Zo+LFm/k2Fjili1ZU/QkaSGVOmmYw3PFbMaWO3u0yLLkwgXPOvNrAAg8ryhmt9/XE4sqEfXZvV4tnh1ELsXSTbFhNvczASmAL9w96uCjSTJ5ozOTXjihgFs3XOAb48p4tM1GqkgEVRloL98M/ulmY03s4lmNhH4s7s/UZ0duvs6YJWZdYqsOpNwz6yXCd8WI/L8UnW+vyZc3rcln67dwfwvtwcVQSQpuPsjhx7Aq4DG05dA9GndkEmjB5KaYlw+bgpTNSZa3KvKFZzamCvme8ATZjYXOAn4I3AXcLaZLQXOjiwHYkjP5qSnpvCMBoMSqVWR4SfqmVku4VtTD5nZ34LOJcmpQ5McJt88iCb1M7lm4nRen7826EhyHKpS4NT4XDHuPjvSjqaHu1/s7lvdfbO7n+nuHSLPgU0aUj8rjfO6NeXF2V9qTByR2lXf3XcA3wYecvc+wFkBZ5Ik1qxBHSaNHki3ZvW45YlZPD7186AjSTVVpcBJyrliLu/bkp37Snl9/rqgo4gkstRIm7vLgH8FHUYEoEFWOk/cMIDTOzXm1y/O53/fXqJRj+NQVQqcpJwrZkDbRrTKzeLpGV8EHUUkkf0eeANY5u4zzKwdsDTgTCLUSQ8x7uo+DOvTgv99eyl3vDRfkzHHmaOOauTuOdEIEmtSUozLClvwlzeX8Pnm3bRulB10JJGE4+7PAc+VW15uZh8FGEnkK2mhFO4e1oP8nAzGvP8Zm3eV8D+Xn0RmWijoaFIFVelFlbRzxQzr05IUQ42NRaLrx0EHEDnEzPjFuZ2548IuvDZ/Hdc+NJ0d+w4EHUuqoCq3qH5W7nEH8E/gd7WYKWY0rZ/JGZ0b82zxKkpKDwYdRyRZWNABRA53/SltuWf4SRSv3Mrl46ayYce+oCPJURy1wHH3b5V7nE14GPX1tR8tNlw9sA2bdpXwmroLikSLGjpITBpyUnMmXtuXzzfvZujYIlZs2h10JDmC6kyfmlRzxZzaPo82jbJ4bIq6CorUlEOdFSp47AQ0063ErNM65vPUjQPYvb+MYWOKmLdaA8LGqqq0wUnquWJSUoyrBrSm+POtGr5bpIa4e46716vgkePumtJZYlrPlg2YNHogmWkhho+fwsdLNwUdSSpQlSs4ST9XzKV9WpKZlsJjU1cGHUVERGJAu/y6PH/LIFrmZvGdh6drguYYVJU2OEk/V0z9rDSG9GzOi5+sYftetZ4XOV5mNqsmthEJUpN6mTxz00B6tWrI95/+hIf/vSLoSFJOVW5Raa4Y4OqBrdl7oIxJMzXLuEgNONHM5h7hMQ/ICzqkyNHUr5PGo9f145tdmvC7f37K3W8s0qjHMaIq97rru/sOM7uB8Fwxv41MkplUujWvT+9WDXhsykquHdSGUIp6sooch85V2EYTwUlcyEwLcf+Vffj1i/O5773P2LSzhDsv6UZqqDr9eKSmVKXAKT9XzK9qOU9Mu/6Udtz65CzeXriec7o2DTqOSNxyd3VLlIQSSjH+eEk38nMyuPedpWzeXcLfR/SiTrpGPQ5KVcpLzRUTcU7XJrRoWIcJHy4POoqIiMQYM+PHZ3fkD0O68s6i9Vz94DS271G7zaBUpZHxc+7ew91viSwvJ9xVPOmkhlK4/pS2FH++lVlfbA06joiIxKCrB7bhvit6M3f1di4dV8S67Rr1OAjVvUGYtHPFXFbYknqZqTzwka7iiNQUM8s2M13Ll4RxfvcCHv5OX9Zs28fQMUUs27Ar6EhJp7oFTtK2sM3OSOXKAa15ff46vti8J+g4InHJzFLM7Aoze8XMNgCLgLVmtsDM7jazDkFnFDleg9rn8fSoAewvPcilY4v4RFf+o6q6BU5S94E71ItqosY8EKmu94ATgNuBpu7e0t0bA6cCU4G7zCypBhSVxNSteX0m3zyQenXSuGLCNN5bvCHoSEmj0gJHc8VUrkm9TC7q2ZxnZqxi6+6SoOOIxKOz3P0P7j7X3Q8eWunuW9x9srsPBZ45nh2Y2blmttjMlpnZbRW8f7qZbTez2ZHHb45nfyKVad0om0mjB9EuP5sbHynm+VkaTy0aKi1wNFfMkY06rR17D5TxkK7iiFRHjpnlVvYAcPdqdz+JtOe5DzgP6AKMMLMuFWz6kbufFHn8vrr7Ezma/JwMnh41gP7tcvnxs3MY/+FnQUdKeBqFqJo6Nc3hnK5NeKhoJTv2qRugyDGaydfnuSv/KK6B7+9HeGiL5e5eAjwNDKmB7xWptpzMNCZe25cLehTwx1cXcecrn3LwYFK3+KhVR7pFpblijuJ7Z3Rg575SHi1aGXQUkbji7m3dvV3k+fBHuxrYRXNgVbnl1ZF1hxtoZnPM7DUz61rZl5nZKDMrNrPijRs31kA8SVYZqSH+PrwXIwe2ZsJHK/jpc3M4UHbw6B+UY3akW00nHmVKBgPq13CeuNKteX3O6NyYBz9ewXdObkt2RtLfuROpEjNr4+4rj/C+Ac3dvbqNFSrq6Xn4r8qzgNbuvsvMzgdeBCrsveXu44HxAIWFhfqVW45LSorxu4u6kp+TwV/eXMKWPSXcf2VvstL1f0hNOtKfpuaKqYLvntGeb99fxONTP+emb5wQdByReHG3maUALxG+LbURyATaA4OBM4HfEr7yUh2rgZblllsAa8pv4O47yr1+1czuN7M8d99UzX2KVJmZ8d0zOpBXN4NfvjCPEROm8dC1fcnNTg86WsKotMDRXDFV07tVQ05pn8eEj5ZzzcA2mndEpArc/dJIo98rgeuAAmAvsBB4BbjT3Y9n+NcZQAczawt8CQwHrii/gZk1Bda7u5tZP8K37Dcfxz5Fjtnwfq3IzU7ne099wrCxRTx6XT9aNMwKOlZCUCPjGvC9M9qzaVcJT0xTTShSVe7+qbv/yt1Pd/dOkZ5MI9z98eMsbnD3UuC7hOfRWwg86+4LzGy0mY2ObDYMmG9mc4B7geHurttPEnXf7NqUx67vz8ad+xk2ZgqL1+0MOlJCsHj+91xYWOjFxTXR4eL4XTFhKovW7eTDnw+mrtriSJwys5nuXhjF/X27gtXbgXnuHpMjosXSeUcSy6J1Oxg5cTp7S8p48Nq+9G2TG3SkqKit846u4NSQn53TiS27S3jwI42LI3IMrgceIHyr6kpgAuG57v5tZlcHGUwk2jo3rcfkmweRl5PBVQ9M461P1wcdKa5VayRjM9tR2eeSVa9WDflmlyZM+Gi5RjcWqbqDwInuPjQyenEXYD/QH/hFoMlEAtCiYRaTRg+ic0E9bnqsmGdmfBF0pLhV3ZGM6x3vjs0sZGafmNm/Isu5ZvaWmS2NPDc83n1E20/P6cTuklLGfKARKkWqqI27l/81dQPQ0d23ABpBU5JSbnY6T97Qn1M65POLyfO4771lxHNzkqAEeYvqB4Qb/x1yG/COu3cA3oksx5WOTXK4pFdzHilaybrtx9VGUiRZfGRm/zKzkWY2EngZ+NDMsoFtwUYTCU52RioPjizkkl7NufuNxfy/f2rU42MVSIFjZi2ACwjfez9kCPBI5PUjwMVRjlUjfnRWRw66c887S4KOIhIPbgUeAk4CehH+t3+ru+9298FBBhMJWloohb9e2pMbTmnLw0Ur+f7Tn7C/NOmHn6uyoK7g/C/wc8L33w9p4u5rASLPjQPIddxa5mZxZf/WPDNjFYvWqamSyJFEumV/DLwLvA18qK7aIv8nJcX49YVd+OX5nfnX3LVc/3Axu/aXBh0rLkS9wDGzC4EN7j6zmp+P+TlhfnhWB3Iy0/ivfy3UfVORIzCzy4DphMekuQyYZmbDgk0lEntGnXYCf720J1OWb2bE+Kls2rU/6EgxL4grOCcDF5nZSsIz/J5hZo8D682sACDyXOEYGO4+3t0L3b0wPz8/WpmPSYOsdH54Vgc+XraJdxfF5FAeIrHiV0Bfdx/p7tcQngX8joAzicSkoX1a8MA1hSzdsJNhY4r4YvOeoCPFtKgXOO5+u7u3cPc2hIdPf9fdryLcuHBkZLORhOeoiVtXDWhNu/xs7nxlISWlmilWpBIphw3otxmNzyVSqcGdG/PkjQPYtvcAQ8cWsWDN9qAjxaxYOpHcBZxtZkuBsyPLcSstlMIdF3Rh+abdPDZVUziIVOJ1M3vDzK41s2sJz0P1asCZRGJa71YNmTR6IGkpxvBxU5nymaZQq0igBY67v+/uF0Zeb3b3M929Q+R5S5DZasLpnfI5tUMe97y9RPdLRSrg7j8DxgM9gJ7AeHfXAH8iR9G+cQ6TbxlE0/qZjJw4nVfnrQ06UsyJpSs4CcfM+O23urD3QBn//eqioOOIxCR3n+zuP3b3H7n7C0HnEYkXBfXr8NzogXRvUZ9bn5zF47pb8DUqcGpZ+8Y5jDqtHZNnrWbqcl1GFAFNBSNSUxpkpfP49f05o1Njfv3ifP7nrSXqvRuhAicKvju4Ay0a1uHXL85Xg2MRan8qGJFkUic9xLir+3BZYQvueWcpv3pxPmUa9VgFTjTUSQ/xhyHdWLZhFxM+Wh50HBERSTCpoRT+NLQHt5x+Ak9O+4Jbn5jFvgPJPeqxCpwoGdy5Med2bcrf313Kqi0au0BERGqWmfHzczvzmwu78PqCdYycOJ0d+5J3zloVOFH024u6EDLj9ufn6R6piIjUiutOacs9w09i1hdbuXzcVDbsSM7Jn1XgRFFB/Tr88oIT+XjZJp6c/kXQcUREJEENOak5E6/ty+ebd/PtMUWs2LQ76EhRpwInyq7o14qT2zfij68sZPVW3aoSEZHacWqHfJ4eNYA9JWUMG1PE3NXbgo4UVSpwoszM+NPQHgDcNlm3qkREpPb0aNGASaMHUic9xPDxU/loaWxOUl0bVOAEoEXDLN2qEhGRqGiXX5fnbx5Eq9wsrnt4Bi/N/jLoSFGhAicgV/RrxSnt87jzlYUs37gr6DgiIpLAGtfL5NnRA+ndqiE/eHo2Ez9eEXSkWqcCJyBmxl8u7Ul6ago/eHq2BgAUEZFaVS8zjUeu68e5XZvy+399yp9fX5TQzSRU4ASoaf1M/jy0B/O+3M5f3lwcdBwREUlwmWkh7ruyN1f0b8X973/GLybPpbQsMX/BVoETsG92bcpVA1ox/sPlfLgkeRp/iYhIMEIpxp0Xd+MHZ3bg2eLV3PTYTPaWJN6oxypwYsCvL+hCxyZ1+fGzc9i4c3/QcUREJMGZGT86uyN/uLgb7y7ewNUPTmPbnpKgY9UoFTgxIDMtxL0jerFz3wG+99SshL1cKCIiseXqAa2574rezF29nUvHTmHt9r1BR6oxKnBiROem9fjjJd2ZunwLf3p9UdBxREQkSZzfvYBHruvHuu37GHp/Ecs27Aw6Uo1QgRNDhvZpwdUDWjPhoxX8a+6aoOOIiEiSGHhCI56+aQAlZc6wsVOY9cXWoCMdNxU4MeaOC7vQu1UDfj5pLkvWJ0YVLSIisa9rs/o8f/MgGtRJ44oJU3lv0YagIx0XFTgxJj01hfuv7ENWeohRjxazdXdiNfoSEZHY1apRFpNuHkSHxjnc8Ggxk2euDjpStanAiUFN62cy9qo+rNm2j9GPz9QggCIiEjV5dTN4atQABrTL5SfPzWHcB58FHalaVODEqMI2ufx5WA+mrdjC7c9rUk4REYmeuhmpTLy2Lxf2KOC/X1vEna98ysGD8fX/UGrQAaRyF/dqzopNu7nnnaW0y8/m1sHtg44kIiJJIiM1xL3De5FXN4MJH61g064S/jysB2mh+Lg2ogInxv3wrA6s2LSbu99YTIuGdRhyUvOgI4mISJJISTF++60u5OdkcPcbi9m8u4QxV/YmOyP2y4f4KMOSmJnx52E96N82l588OyfuW7WLRIuZnWtmi81smZndVsH7Zmb3Rt6fa2a9g8gpEuvMjFsHt+dPQ7vz8dKNXPHANLbEQQcYFThxIDMtxAMjC+lckMPox2cyY+WWoCOJxDQzCwH3AecBXYARZtblsM3OAzpEHqOAMVENKRJnLu/binFXF7Jo7Q6GjS1i9dY9QUc6IhU4cSInM42Hv9OP5g3qcN3DM/h0zY6gI4nEsn7AMndf7u4lwNPAkMO2GQI86mFTgQZmVhDtoCLx5OwuTXj8hv5s2rmfoWOKWLQudv8vUoETR/LqZvDo9f2om5HKNROnaSBAkco1B1aVW14dWXes2wBgZqPMrNjMijdu3FijQUXiTd82uTw3ehAAl46dwvQVsXlXQQVOnGnRMIvHb+iPmTFi/FQWr1ORI1IBq2Dd4X1cq7JNeKX7eHcvdPfC/Pz84w4nEu86Nc1h8s2DyM/J4OoHp/HmgnVBR/oPKnDi0An5dXl61ABCKcaICVNj+hKhSEBWAy3LLbcADp/grSrbiEglWjTMYtLoQZxYUI/Rj8/k6elfBB3pa6Je4JhZSzN7z8wWmtkCM/tBZH2umb1lZksjzw2jnS2enJBfl2duGkh6KIUR46eqTY7I180AOphZWzNLB4YDLx+2zcvANZHeVAOA7e6+NtpBReJZbnY6T97Yn9M65nPb8/P4x7tLY2Zg2iCu4JQCP3H3E4EBwK2R3g23Ae+4ewfgnciyHEHbvGyeHjWAzLQQl4+P3fugItHm7qXAd4E3gIXAs+6+wMxGm9noyGavAsuBZcAE4JZAworEuaz0VCZcU8glvZrzlzeX8LuXF8TEqMdRL3Dcfa27z4q83kn45NOccI+GRyKbPQJcHO1s8ahNXjbPjR4Y0/dBRYLg7q+6e0d3P8Hd74ysG+vuYyOv3d1vjbzf3d2Lg00sEr/SQin89dKejDqtHY9M+ZzvPf0J+0vLAs0UaBscM2sD9AKmAU0OXR6OPDeu5DPqzXCYQ/dBO8fofVAREUl8KSnGL88/kV+e35lX5q7lOw/NYOe+A8HlCWrHZlYXmAz80N2r3IBEvRkqlpudzlM39ufUDuH7oH97c3FMXCIUEZHkMuq0E/jbZT2ZvmILIyZMZePO/YHkCKTAMbM0wsXNE+7+fGT1+kODbEWeNSfBMcpKT+WBkYVc2qcF9767jFufnMWektKgY4mISJL5du8WTBhZyGcbdjNsbBGfb94d9QxB9KIy4EFgobv/rdxbLwMjI69HAi9FO1siSAul8OdhPfjV+Sfy+oJ1XDp2Cmu37w06loiIJJnBnRrz5I392bH3AEPHTGH+l9ujuv8gruCcDFwNnGFmsyOP84G7gLPNbClwdmRZqsHMuPG0dkwc2ZfPN+/hon/8Wz2sREQk6nq1ashzoweRHjKGj59K0WeborbvIHpRfezu5u493P2kyONVd9/s7me6e4fIs/5HPk6DOzfm+VsGkZ0eYsSEqYz94DO1yxERkahq37guk28ZRLMGmVw7cQavzovOcFMayTjBdWySwz+/dwrndm3KXa8t4sZHi9m2J/anuRcRkcRRUL8Oz900iB4t6nPrk7N4bMrKWt+nCpwkkJOZxj+u6MX/u6grHy7dyAX3fsy05ZuDjiUiIkmkflYaj9/QnzM7N+GOlxbwtzcX1+qoxypwkoSZMXJQGyaNHkRqyBg+YSp/fHVh4AMxiYhI8shMCzH2qt5cVhju7fvLF+bX2r5U4CSZni0b8Or3T2VEv1aM/3A5F/3935rHSkREoiY1lMKfhvbg1sEn8FQtDkyrAicJZWek8sdLuvPQtX3ZsqeEIfd9zN/eXMy+A7qaIyIitc/M+Nk5nXn0un61tg8VOElscOfGvPHD07igewH3vruM8+75iKJl0evCJyIiye20jrU3I4EKnCSXm53O/w7vxWPX9+OgO1c8MI0fPzubTbuCGVpbRESkJqjAEQBO7ZDPGz88je8Obs/Ls9cw+O73GffBZ2qELCIicUkFjnwlMy3ET8/pxOs/PI2+bXP579cWcdbfPuDVeWtrtSufiIhITVOBI/+hfeO6TLy2L49d34+stFRueWIWl46dEtUhtkVERI6HChyp1Kkd8nnl+6dw5yXd+GLLHq6YMI0R46cyY6Vm0RARkdimAkeOKDWUwpX9W/Phzwfzmwu7sHTDLi4dO4WrH5zG9BVbdOtKRERikgocqZLMtBDXndKWj34+mF+e35kFa3Zw2bgpXHJ/Ea/NW0uZJvEUEZEYogJHjkmd9BCjTjuBf//iDP4wpCtb95Rw8xOzOPOv7/PY1M/Zvb806IgiIiIqcKR66qSHuHpgG979yemMubI3DbLSuePF+fT/4zvc8eJ8Fq3T9A8iIhKc1KADSHwLpRjndS/g3G5NmfXFNp6Y9jnPFK/isamfU9i6IVf0b8V53Qqokx4KOqqIiCQRFThSI8yMPq0b0qd1Q+64oAuTZ63miWlf8ONn53DHi/M5t1sBl/RqzsATGhFKsaDjiohIglOBIzWuYXY6N5zajutPacu0FVt48ZMveWXeWibPWk2TehkMOak5F/YooHvz+pip2BERkZqnAkdqjZkxoF0jBrRrxO8u6sq7izbw/KwvmfjxCsZ/uJxm9TP5ZtemnNO1KX3bNCQ1pCZhIiJSM1TgSFRkpoU4v3sB53cvYNueEt5ZuIHXF6zjqelf8HDRShpmpXHmiU34Rsd8Tu2QR4Os9KAji4hIHFOBI1HXICudoX1aMLRPC/aUlPLB4o28vmAdby5Yx6SZq0kx6NGiAad1zOcbHfPp2aK+ru6IiMgxUYEjgcpKT+W87gWc172A0rKDzFm9nQ+XbOTDpRv5x7tLufedpeRkpNKnTUP6tc2lf9tcujdvQHqqCh4REamcChyJGamhlK96Yv3o7I5s21PCv5dtpuizTUxfsYU/L14MQEZqCr1aNaBfm1x6tmxAjxYNyM/JCDi9iIjEEhU4ErMaZKVzQY8CLuhRAMDmXfuZsXIr01dsYfrKzfzjvWUcmiGieYM69GhRnx4tGtCzRX26Nq9P/TppAaYXEZEgqcCRuNGobgbndmvKud2aArB7fynzv9zO3NXbmbN6G3NXb+e1+eu+2r5Z/Uw6Nc2hU9N6dG6aQ6emOZyQX1e3t0REkoAKHIlb2Rmp9G/XiP7tGn21buvuEuZ+uZ0Fa7azeN1OFq/bycfLNnGgLHypJzXFaJefTbu8urTJy6ZdXjZt8rJpm5dNXt10jcsjIpIgVOBIQmmYnc43Ir2vDjlQdpDlG3ezaN0OFq/byZL1u1i6YSfvLFr/VeEDkJORSptIwdOiYR2aN6hD84Z1aBF5zkrXPxcRkXihM7YkvLRQSuRWVc7X1peWHWTNtn0s37SLFZt2s3LTbpZv2s3sVVt5bd5aSg/617ZvmJVG80jhU1C/Do3rZdA4J5PGORlfvW6YlaarQCIiMUAFjiSt1FAKrRpl0apRFqd3+vp7ZQedDTv38eXWvXy5LfKIvF6+cTf/XraZXftL/+M700JGft0M8uuFC5+8uuk0yEonNyudBllp5GZHlrPTaZiVRr3MNFI0N5eISI1TgSNSgVCKUVA/fKWmsJJt9pSUsnHnfjbs3M+GHfvZsHPf116v2rKH2au2sW1PydduhZWXYtAwUvzUr5NGTmYadTNTqZeZSk5mGjkZqeREXtfNDL+ul5lGTmYqdTNSyUpPJTMtRVeNREQOE3MFjpmdC9wDhIAH3P2ugCOJVCgrPZXWjVJp3Sj7iNu5O7v2l7JtzwG27C5hy54Stu0pYcvuA5HnErbuKWHH3lK27Slh1ZY97NhXys59B9hfevCoOcygTlqIrPQQddJDZKWlhp/TQ9RJC331Ois9sj6yLiM1hYzUEBlp//ecCMwsF3gGaAOsBC5z960VbLcS2AmUAaXuXlktKyJxKKYKHDMLAfcBZwOrgRlm9rK7fxpsMpHqM7Pw1ZjMNFrmZh3TZ0tKD7Jz3wF27S9l575Sduw7wM594de79h1gz4Ey9paUsaekjL1fvS4NL5eUsX3vga/e31NSyt4DZZVeTUogtwHvuPtdZnZbZPkXlWw72N03RS+aiERLTBU4QD9gmbsvBzCzp4EhgAocSUrpqSk0qptBo7o1N1LzgbKD7DtQxv7Srz/vO3CQwj/V2G6CNAQ4PfL6EeB9Ki9wRCRBxVqB0xxYVW55NdC//AZmNgoYFVncb2bzo5QtmvKARPutUscUHzodfZOY18Td1wK4+1oza1zJdg68aWYOjHP38ZV9YRKcdxLx77KOKX7Uynkn1gqcilpKfu16euQkNB7AzIoT8b55Ih6Xjik+mFlx0BmqwszeBppW8NavjuFrTnb3NZEC6C0zW+TuH1a0YaKfd3RM8SERjwlq77wTawXOaqBlueUWwJqAsohIjHL3syp7z8zWm1lB5OpNAbChku9YE3neYGYvEL5FXmGBIyLxJ9a6TcwAOphZWzNLB4YDLwecSUTiy8vAyMjrkcBLh29gZtlmlnPoNfBNINFuO4kktZgqcNy9FPgu8AawEHjW3Rcc4SOV3jOPc4l4XDqm+JAIx3QXcLaZLSXcI/MuADNrZmavRrZpAnxsZnOA6cAr7v56Fb8/Ef6MDqdjig+JeExQS8dl7gnfZVRERESSTExdwRERERGpCSpwREREJOHEbYFjZuea2WIzWxYZrTSmmdlKM5tnZrMPdYkzs1wze8vMlkaeG5bb/vbIsS02s3PKre8T+Z5lZnavRXESIjObaGYbyo8BUpPHYGYZZvZMZP00M2sT0DH9zsy+jPysZpvZ+XF2TC3N7D0zW2hmC8zsB5H1cf2zCprOOTrn1PJxxe15J2bPOe4edw/C81R9BrQD0oE5QJegcx0l80og77B1fwZui7y+DfhT5HWXyDFlAG0jxxqKvDcdGEh4zKDXgPOieAynAb2B+bVxDMAtwNjI6+HAMwEd0++An1awbbwcUwHQO/I6B1gSyR7XP6sgHzrn6JwTheOK2/NOrJ5zAv9HWM0/zIHAG+WWbwduDzrXUTJXdLJZDBSU+wuyuKLjIdyrbGBkm0Xl1o8gPAJrNI+jzWH/KGvsGA5tE3mdSnjETgvgmCo70cTNMR2W+yXCvYni/mcV1EPnnP/8OxDF40i4c04lx5Uw551YOefE6y2qiqZ0aB5Qlqo6NCz8TAsP+w6HDSkPHBpSvrLjax55ffj6INXkMXz1GQ8PGbAdaFRryY/su2Y2N3Ip+dBl1bg7pshl3F7ANBL3ZxUNOud8fX2QEvnvcdyfd2LpnBOvBc5Rp3SIQSe7e2/gPOBWMzvtCNtWdnzxdNzVOYZYOb4xwAnAScBa4K+R9XF1TGZWF5gM/NDddxxp0wrWxexxBSQej1fnnPj6exz3551YO+fEa4ETd1M6eLlh4YEXCA8Lv97CQ8ljXx9SvrLjWx15ffj6INXkMXz1GTNLBeoDW2oteSXcfb27l7n7QWAC4Z/V1/JFxOwxmVka4RPNE+7+fGR1wv2sokjnnK+vD1JC/j2O9/NOLJ5z4rXAiaspHazyYeErG1L+ZWB4pNV4W6ADMD1yiW+nmQ2ItCy/hgqGoY+ymjyG8t81DHjXIzdco+nQP8iIS/i/Ifzj4pgiGR4EFrr738q9lXA/qyjSOUfnnFoVz+edmD3nRKvRUS00YjqfcEvtz4BfBZ3nKFnbEW4xPgdYcCgv4fuH7wBLI8+55T7zq8ixLaZcrwWgkPBf/M+AfxDdhmNPEb50eoBwNX19TR4DkAk8Bywj3JK+XUDH9BgwD5gb+UdVEGfHdArhS7dzgdmRx/nx/rMK+qFzjs45tXxccXveidVzjqZqEBERkYQTr7eoRERERCqlAkdEREQSjgocERERSTgqcERERCThqMARERGRhKMCR46JmZVZeKbb+Wb2TzNrUME2bcxsr5nNruY+TjWzT63cTLsikpx0zpHqUoEjx2qvu5/k7t0IjyJ5ayXbfebuJ1VnB+7+EeExFEREdM6RalGBI8djClWYeC/y29UiM3skMpHcJDPLirzX18yKzGyOmU0/NPqqiEgFdM6RKlOBI9ViZiHgTKo+XH0nYLy79wB2ALdEhrx/BviBu/cEzgL21kZeEYlvOufIsVKBI8eqTuQ+92YgF3irip9b5e7/jrx+nPDQ3p2Ate4+A8Ddd7h7aQ3nFZH4pnOOVIsKHDlWeyP3uVsD6VR+P/xwh88J4oBVsF5EpDydc6RaVOBItbj7duD7wE/NLK0KH2llZgMjr0cAHwOLgGZm1hfAzHLMLLVWAotIXNM5R46VChypNnf/hPBsxcOrsPlCYKSZzSV8mXmMu5cAlwN/N7M5hC89Z9ZWXhGJbzrnyLFQ5SrHxN3rHrb8rSp+9KC7j67g+2YAA2oim4gkHp1zpLp0BUdqQxlQ/3gG3QL+CWyqyVAikrB0zpH/YO5qbyUiIiKJRVdwREREJOGowBEREZGEowJHREREEo4KHBEREUk4KnBEREQk4fx/QJkYQgxv91MAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "Fig,(LinPlot,LogPlot)=plt.subplots(1,2,figsize=(8,4))\n", "LinPlot.plot(R,I)\n", "LinPlot.set_ylim(0,120)\n", "LinPlot.set_ylabel('I [Lsun/pc^2]')\n", "\n", "LogPlot.plot(R,np.log10(I))\n", "LogPlot.set_ylim(-0.5,2.1)\n", "LogPlot.set_ylabel('log(I) [Lsun/pc^2]')\n", "\n", "for panel in [LinPlot,LogPlot]:\n", " panel.set_xlim(0,20000)\n", " panel.set_xlabel('R [pc]')\n", "\n", "Fig.tight_layout() " ] }, { "cell_type": "markdown", "id": "3c63839a", "metadata": {}, "source": [ "OK, so how do we work out computationally where the function hits I = 30 Lsun/pc^2?\n", "\n", "First way is to use a root finder. If we have a function f(x) and want to know where does f(x)=c, that's the same as wanting to know when f(x)-c=0. So we can make a new function g(x)=f(x)-c and find the root using scipy optimize." ] }, { "cell_type": "code", "execution_count": 5, "id": "a05d7c2f", "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The luminosity density hits a value of 30.00 Lsun/pc^2 at R = 4815.891 pc\n" ] } ], "source": [ "from scipy import optimize\n", "\n", "# define the function. Note that inside the function I have not defined the variables I0, h, or Iwant,\n", "# so it uses the values that are currently in memory. That may or may not be a good idea.....\n", "def gfunc(R):\n", " I = I0*np.exp(-R/h) \n", " g = I - Iwant\n", " return g\n", "\n", "sol = optimize.root(gfunc, 5000) # give it the function name, and a guess as to the root value\n", "Rwant = sol.x[0]\n", "\n", "print('The luminosity density hits a value of {:.2f} Lsun/pc^2 at R = {:.3f} pc'.format(Iwant,Rwant))" ] }, { "cell_type": "markdown", "id": "11370351", "metadata": {}, "source": [ "Or, we can use the simple, lazy Mihos method: if you can plot it, you can solve it. Since to make the plot, you had an array of (R,I) values, just find the index where I is closest to the value you want, and ask what R that corresponds to." ] }, { "cell_type": "code", "execution_count": 6, "id": "e77c006c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The luminosity density hits a value of 30.00 Lsun/pc^2 at R = 4820.000 pc\n" ] } ], "source": [ "idx=np.argmin(np.abs(I-Iwant)) # this gives you the index in the I array where I is closest to Iwant\n", "\n", "Rwant=R[idx] # this gives you the R value that corresponds to that index.\n", "\n", "print('The luminosity density hits a value of {:.2f} Lsun/pc^2 at R = {:.3f} pc'.format(Iwant,Rwant))" ] }, { "cell_type": "markdown", "id": "c932ab19", "metadata": {}, "source": [ "Since the R,I values were only calculated in steps of 10 parsecs, your answer is only accurate to that level. \n", "\n", "If you want a more accurate value, you need to have a more finely stepped array of (R,I). Make the array as fine as you want, but the finer it is the more work the computer does....." ] }, { "cell_type": "code", "execution_count": 7, "id": "7a6c5cdc", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The luminosity density hits a value of 30.00 Lsun/pc^2 at R = 4815.891 pc\n" ] } ], "source": [ "# first make a more finely stepped array of R,I values\n", "R=np.arange(0,20000,0.001)\n", "I=I0*np.exp(-R/h)\n", "\n", "idx=np.argmin(np.abs(I-Iwant)) # this gives you the index where I is closest to Iwant\n", "\n", "Rwant=R[idx] # this gives you the R that corresponds to that index.\n", "\n", "print('The luminosity density hits a value of {:.2f} Lsun/pc^2 at R = {:.3f} pc'.format(Iwant,Rwant))" ] }, { "cell_type": "markdown", "id": "29c97406", "metadata": {}, "source": [ "Let's plot it to make sure it worked!" ] }, { "cell_type": "code", "execution_count": 8, "id": "89364e1d", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEYCAYAAABRMYxdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABJBUlEQVR4nO3dd3xV9fnA8c9zb/YgEBIgQNhL9ggQnKhYwYpoHYCACAhUsWodrdpWW1dt1VatMyACgoA/tXUvrOCAQAAB2TJk7xE2ZHx/f5xLCRhCSO6959xznvfr9X2d+70595znkOThyRnfrxhjUEoppZRyE5/dASillFJKBZsWOEoppZRyHS1wlFJKKeU6WuAopZRSynW0wFFKKaWU60TZHUBlpKWlmQYNGtgdRvCsWGEtmze3Nw6lSpg3b95OY0y63XHYQXOMUqEXqhwT0QVOgwYNmDt3rt1hBM/EidZy4EB741CqBBFZZ3cMdtEco1TohSrHRHSB4zqadJRSoaQ5RnmI3oPjJIcOWU0ppUJBc4zyED2D4yRXXGEtp0+3NQyllEtpjlEeogWOk9x6q90RKKXcTHOM8hAtcJykb1+7I1BKuZnmGOUheg+Ok+TnW00ppUJBc4zykJAVOCIyVkS2i8jiEu+lisgXIvJjYFmtxNceEJFVIrJCRC4PVVyO1qeP1ZRSZ0VEMkXkKxFZJiJLROTOUtYREXk+kGcWiUhHO2K1leYY5SGhvEQ1DngBmFDivfuBL40xT4rI/YH+70WkJdAPaAXUBqaJSDNjTFEI43OeO+6wOwKlIlUhcI8xZr6IJAPzROQLY8zSEuv0ApoGWlfg5cDSOzTHKA8JWYFjjPlaRBqc8nYfoHvg9XhgOvD7wPtTjDFHgbUisgroAswqax/rdh2isKiYKL9LrrT96ld2R6BURDLGbAG2BF7vF5FlQB2gZIHTB5hgjDFArohUFZGMwGdLtW7XQfYfKSA5LjqU4YeP5hjlIeGuDGoeTyaBZY3A+3WADSXW2xh472dEZISIzBWRufuOFLBqx4GQBhxWO3daTSlVYYE/rDoAs0/5UrnyzMk5ppBBr81h35GCkMUbVppjlIc45dSHlPKeKW1FY0yOMSbLGJMFsGiDi26Yu+46qymlKkREkoB3gLuMMftO/XIpH/lZnimZY+qnJrBkcz6Dxswm/5ALihzNMcpDwl3gbBORDIDAcnvg/Y1AZon16gKbz7QxnwgLNu4Ndoz2ueceqymlzpqIRGMVN5OMMe+WsspZ55kq8dG8PKATy7bsZ8Bruew9dCx4AdtBc4zykHAXOO8DgwOvBwPvlXi/n4jEikhDrJsA55xpYwkxfha5qcDp3dtqSqmzIiICvAYsM8b84zSrvQ/cFHiaKhvIL+v+m+N6tKzJq4M6sXLbAW4cPZs9ByO4yNEcozwklI+JT8a6Sbi5iGwUkWHAk8BlIvIjcFmgjzFmCfAW1g2BnwKjyvMEVXyMn+Vb9nOkwCUPW23dajWl1Nk6DxgEXCIiCwLtChH5tYj8OrDOx8AaYBUwGritvBu/uEUNRt+UxaodB+g/OpddB44G/QDCQnOM8hCxHiiITM1atTPHej/Bu7edS8d61c78Aafr3t1a6jwxykFEZN7xe968Jisry8ydO/d//W9/3Mmw8XnUr57ApFuySU+OtTG6CtAcoxwoVDnGKTcZV0h8jPWU+8INe+0NJFjuv99qSilHOr9pGq/f3Jn1uw/Rf3Qu2/cfsTuks6M5RnlIRBc40X6hRnIsiza65Emqnj2tppRyrHObpDFuSBc27z1Mv5xctu2LoCJHc4zykIgucADaZVZ1zxmcDRusppRytOxG1Rk3pAvb8o/QLyeXrfkRUuRojlEeEvkFTt0U1uw8SP5hF4xRMWiQ1ZRSjtelYSoThnVhx/6j9M2Zxea9h+0O6cw0xygPifwCJ7MqAD+44TLVH/9oNaVUROhU3ypydh84Rt+cWWzcc8jukMqmOUZ5SMQXOG3rVAVgoRvGw+nRw2pKqYjRsV41Jt7SlfxDBfR9NZcNux1c5GiOUR4S8QVOSkI0DdMS3XEfzpo1VlNKRZR2mVWZdEs2B44W0vfVWazbddDukEqnOUZ5SMQXOABt66a440mqoUOtppSKOG3qpvDm8K4cLiii76u5rN3pwCJHc4zyEFcUOO3qVmXrviOR9bhmaf7yF6sppSJSq9opvDk8m2NFxfR9dRardxywO6STaY5RHuKOAidwo/H36/fYG0hlXXSR1ZRSEeucjCpMHp5NsTH0y8ll1fb9dod0guYY5SGuKHBa16lCjN/H/PV77Q6lclassJpSKqI1r5XM5OHZGAP9cnJZuc0hRY7mGOUhrihwYqP8tK5ThfnrIvwMzsiRVlNKRbymNZOZMiIbnwj9cnJZtmWf3SFpjlGe4ooCB6BT/Wos2pTPscJiu0OpuCeesJpSyhWa1Ehi6shuxPh93Dg6lyWbbX4YQnOM8hDXFDgd61XjWGGx/QmkMs4912pKKddomJbI1JHZxEf7uXH0bBZvsjFHaY5RHuKeAqd+NQDmRfJlqsWLraaUcpX61ROZOrIbSbFR3Dg6l0V2DUyqOUZ5iGsKnJpV4qhTNZ7vI/lG49tvt5pSynUyUxOYMiKblIRoBoyZbc9Tn5pjlIdE2R1AMHWqX405a3fbHUbFPfWU3REopULIKnK60T8nl0GvzWH80M50qp8avgA0xygPcc0ZHLAKnK37jkTGrL6l6dzZakop16pTNZ6pI7NJS4rhptfmkPdTGP8o0xyjPMRVBU7HehF+H86CBVZTSrlaRko8U0d2o2aVOAaPnUPuml3h2bHmGOUhripwWmQkEx/tj9wC5667rKaUcr2aVeKYMiKb2lXjGfJ6HjNX7wz9TjXHKA9x1T040X4fbeumRO6UDc8+a3cESqkwqlEljsnDsxkwJpeh4/IYc1Nnzm+aFrodao5RHuKqMzhg3YezZPM+jhQU2R3K2Wvf3mpKKc9IT45l8vBsGlRPZNj4PGas3BG6nWmOUR7iugKnY71qFBYbFm2MwAH/8vKsppTylOpJsbw5PJtG6UkMnzCXr5ZvD82ONMcoD3FfgRMY8C+sTyYEy333WU0p5TmpiTFMHt6VZjWTGPnGPKYt3Rb8nWiOUR7iugInNTGGZjWTInM8nBdesJpSypOqJsQwaVg2LTKSuXXSPD5fsjW4O9AcozzEdQUOQJeGqcxbt4eiYmN3KGendWurKaU8KyUhmjeGdaVV7RRumzSfTxdvCd7GNccoD3FlgdO5QSoHjhaybMs+u0M5OzNnWk0p5Wkp8dG8MawL7TKrMurN7/loUZCKHM0xykNcWeB0aWgNfT470i5TPfig1ZRSnpccF834oV3oWK8qd0z5nvcWbKr8RjXHKA9x1Tg4x2WkxFMvNYE5a3cx7PyGdodTfq++ancESikHSYqNYtyQLgwZl8dvpy6g2Biu6VC34hvUHKM8xJVncMA6izNn7W6MiaD7cJo3t5pSSgUkxkYxbkhnujaszt1vLeTteRsrvjHNMcpDXF3g7DlUwKrtB+wOpfxmzLCaUkqVkBATxdibO3Ne4zTue3shb+VtqNiGNMcoD3HlJSqAriXuw2laM9nmaMrp4Yet5fTptoahlHKe+Bg/YwZnMeKNefzunUUUFhtu7Frv7DaiOUZ5iGsLnHqpCdSsEsuctbsZmF3f7nDKZ+xYuyNQSjlYXLSfnEGduHXiPB789w8UGcOgs8lvmmOUh7i2wBERujSs/r/7cETE7pDOrFEjuyNQSjlcXLSfVwZ14raJ8/nTfxZTVFTMzeeV82EKzTHKQ2y5B0dEfisiS0RksYhMFpE4EUkVkS9E5MfAslpl99OlYSpb9x1hw+7DwQg79KZNs5pSSpUhNsrPywM7cVnLmvz5g6WM+WZN+T6oOUZ5SNgLHBGpA9wBZBljWgN+oB9wP/ClMaYp8GWgXykn7sPZVdlNhcdjj1lNKaXOICbKx0sDOtKrdS0e+2gZOV+vPvOHNMcoD7HrKaooIF5EooAEYDPQBxgf+Pp44OrK7qRJehLVEqLJXRMhA/698YbVlFKqHKL9Pp7v34Ffts3giY+X8/L0MxQ5mmOUh4T9HhxjzCYReRpYDxwGPjfGfC4iNY0xWwLrbBGRGqV9XkRGACMA6tUr+wkCn0/o1rg6s1bvjIz7cDIz7Y5AKRVhov0+nuvbHr8If/t0OUXFxdx+SdPSV9YcozzEjktU1bDO1jQEagOJIjKwvJ83xuQYY7KMMVnp6elnXP/cxmlszj/CT7sOVTjmsPn0U6sppdRZiPL7+Gff9lzToQ5Pf76SZ6etLH1FzTHKQ+x4iqoHsNYYswNARN4FzgW2iUhG4OxNBrA9GDs7t3F1AL5btZOGaYnB2GToPPmktezZ0944lFIRx+8Tnr6+HT4Rnp32I8XFht9e1uzkM9eaY5SH2FHgrAeyRSQB6xLVpcBc4CAwGHgysHwvGDtrmJZIRkocs1bvcv54OFOm2B2BUhFJRMYCVwLbAw8vnPr17lg5ZW3grXeNMY+ELcAw8fuEp65rS5RPeP6/qygsNtx3efMTRY7mGOUhdtyDM1tE3gbmA4XA90AOkAS8JSLDsIqg64OxPxHrPpyvlm+nuNjg8zn4PpxateyOQKlINQ54AZhQxjrfGGOuDE849vH5hL/+qg0+n/DS9NUUGcP9PVtYRY7mGOUhtgz0Z4x5GHj4lLePYp3NCbrzGqfx7vxNLN+6n5a1q4RiF8HxwQfWsndve+NQKsIYY74WkQZ2x+EUPp/w+NWtifIJr85YQ1GR4Q+/PAf58ENrBc0xygNcO5JxSec2se7Dmbl6p7MLnGeesZaafJQKhW4ishBrWIp7jTFLSlvpbJ7UdDKfT3ikTyv8PmHMt2spLDY8/I9nENAcozzBEwVORko8jdISmbl6F7dc4OChyt9+2+4IlHKr+UB9Y8wBEbkC+A9Q6rPUxpgcrMvmZGVlmbBFGAIiwsO9W+ITYex3a4kd+ii/79nCtgHQlAonz/ycd2tcndlrdlFQVGx3KKeXlmY1pVRQGWP2GWMOBF5/DESLyJl/2VasgHHjrNcFBdC9O0ycaPUPHbL6U6da/fx8q//uu1Z/506rf/zS89atVv/4Y9obNlj941MnrFlj9WfMOLHv7t1h5kyrv3ix1c/Ls/oLFlj9BQusfl6e1V+82OrPnAnduyMrV/KnK8/hsWq7uPiPt/HMO3MpLjbWfrt3t+IAK67u3a04wYq7e3frOMA6ru7dreME67i7d7f+HcD6d+ne3fp3AuvfrXv3E/+Wo0dDjx4n+i+9BL16neg/9xxcddWJ/tNPw7XXnug/+ST063ei/+ijMLDECCMPPQRDhpzoP/AAjBhxon/vvTBq1In+XXdZ7bhRo6x1jhsxwtrGcUOGWPs4buBAK4bj+vU78ZQaWLE//fSJ/lVXWcd4XK9e1r/BcT16WP9Gx3XvHvE/e6xYYfVnzLD6awJTipz6sxcinilwzmuSxsFjRSzamG93KKf37rsnfkCVUkEjIrUk8CiRiHTByn0RModL5YkIA7rWo7EcZu+/P+DBf/9gFTlKuZgYE7k/5FlZWWbu3LnlWnf3wWN0fPQL7rmsGb+59DSjfNrt+F8706fbGYVSJxGRecaYLLvjKIuITAa6A2nANqyHGKIBjDGviMjtwK1YT24eBu42xsw803bPJsdEAtO9Oxv3HOKCXg9zXae6/O3atvid/GSp8oRQ5RhP3IMDkJoYQ6vaVfhm1U7nFjjvBWXoH6U8xxjT/wxffwHrMXJPk/feIxO4K2/b/wYDfOr6dlrkKFfyTIEDcEHTdMZ8s4b9RwpIjou2O5yfS0mxOwKllJsFcsxdPVLwi/DMFyspMoZnrm9HlN8zdywoj/DUT/RFzdIpLDbMWu3QS+9Tp564aUwppYKtRI75zaVNue/y5ry3YDN3Tl3g7AcwlKoAT53B6VS/Gokxfmas3MEvWjlwRM+XX7aWffvaG4dSyp1OyTGjLm5ClE/46yfLKS42PN+/A9F6Jke5hKcKnJgoH90apzFj5Q6MMSdPQucEH39sdwRKKTcrJceMvKgxfp/w2EfLuP3N+fyrf0diorTIUZHPcz/FFzVPZ+Oew6zdedDuUH4uIcFqSikVCqfJMbdc0Ig/927JZ0u2cdukeRwtLLIhOKWCy3sFTtN0AGas3GFzJKWYOPHEQE5KKRVsZeSYm89ryKN9WjFt2XZ+/cY8jhRokaMim+cKnHrVE2iYlsjXTixwxoyxmlJKhcIZcsygbg144po2fLViByO1yFERzlP34Bx3UbN0puSt50hBEXHRfrvDOeGLL+yOQCnlZuXIMTd2rYffB/e/+wPDJ8wlZ1AW8TEOypNKlZPnzuAAXNgsjSMFxeT9tNvuUE4WHW01pZQKhXLmmL6d6/H3a9vy7aqdDBufx+FjeiZHRR5PnsHJblSdGL+Pr1fu4ILAPTmOcHxitZtvtjMKpcJORO4ux2oHjTGvhjwYNzuLHHN9ViZRfuGetxYyZNwcXhvcmcRYT/6XoSKUJ8/gJMRE0blhNaavcNh9OOPGnUhASnnLfUASkFxGu8e26NziLHPMNR3q8s++7ZmzdjdDXs/jwNHCkIWmVLB5thy/uHkNHvtoGet3HaJedYc8mq2TbCrvesMY80hZK4hIYriCca0K5Jg+7evg9wl3TlnA4LFzGDekszOnulHqFJ48gwPQ45yaAHy5fJvNkSiljDG/C8Y6KjSubFubF/p3YOGGvQx6bQ77jhTYHZJSZ+TZAqdBWiKN0xP5ctl2u0M5YfRoqynlQSLSQkQuFZGkU97vaVdMrlOJHNOrTQYvDujIks35DBozm/zDWuQoZ/NsgQPWWZzZa3ex3yl/jehkm8qjROQO4D3gN8BiEelT4stP2BOVC1Uyx1zeqhYvD+jEsi37GThmNnsPHQticEoFl6cLnEvPqUlBkeHrlTvtDsUybZrVlPKe4UAnY8zVQHfgTyJyZ+BrDps0LoIFIcf0aFmTVwd1YsW2/dw4ejZ7DmqRo5zJ0wVOx3pVqZoQzZfL9D4cpWzmN8YcADDG/IRV5PQSkX+gBY7jXNyiBqNvymLVjgP0H53LrgNH7Q5JqZ/xdIET5fdxcfMafLViO0XFxu5w4KWXrKaU92wVkfbHO4Fi50ogDWhjV1CuE8Qcc1GzdMYO7szanQfpPzqXHfu1yFHO4ukCB+DSc2qw51AB89fvsTsU+OADqynlPTcBW0u+YYwpNMbcBFxoT0guFOQcc37TNF6/uTPrdx+i/+hctu8/ErRtK1VZni9wLmyWTpRPmOaEy1SffGI1pTzGGLPRGLP1NF/7LtzxuFYIcsy5TdIYN6QLm/cepl9OLtv2aZGjnMHzBU6VuGi6Nkp11uPiSnmUiPzZ7hjU2ctuVJ1xQ7qwLf8I/XJy2ZqvRY6y32lHMvbS3DCXtqjJIx8uZe3OgzRMs3Gw1Oees5Z33ln2ekq5jIj4gNGA/qURSiHMMV0apjJhWBcGj82jb84sJg/PpnbV+KDvR6nyKusMjmfmhrm8dS0APlm8xd5AvvzSakp5zwfAbmPMA3YH4mohzjGd6ltFzu4Dx+ibM4uNew6FbF9KnUlZc1F5Zm6YOlXjaVc3hU8Xb+W27k3sC+T99+3bt1L2ygIetzsI1wtDjulYrxoTb+nKoNdm0/fVXKaMyCYz1SHz/SlPOe0ZHK/NDdOzdQaLNubrXxxK2eNi4FUR6Wp3IKry2mVWZdIt2Rw4WkjfV2exbtdBu0NSHlTmTcZemhumV+Ay1aeLS32QIzyeftpqSnmMMWYpcDnwlN2xuFoYc0ybuim8ObwrhwuK6JeTy9qdWuSo8DptgeO1uWEapCXSolayvQXOrFlWU8qDjDGbgV/aHYerhTnHtKqdwpvDszlaWEy/nFms3nEgbPtWqqwzOJ6bG6ZX6wzmrd/DdrvGcXjnHasp5VHGmP12x+BqNuSYczKqMHl4NkXFhn45uazart9iFR5lFTghmxtGRKqKyNsislxElolINxFJFZEvROTHwLJaZfZREb3a1MIY+GyJjWdxlPKwwO9+1RL9aiLymY0hqSBoXiuZycOzMQb65eSycpsWOSr0yipwQjk3zHPAp8aYFkA7YBlwP/ClMaYp8GWgH1ZNayTRKD2RT+y6TPXkk1ZTyrvSjDF7j3eMMXuAGvaF4zI25pimNZOZMiIbnwj9cnJZvnWfLXEo7yirwAnJ3DAiUiXw+dcC2zwWSGh9gPGB1cYDV1d0H5WIjV6tazF77W52HzwW7t3DggVWU8q7ikWk3vGOiNQHHDATrkvYnGOa1Ehi6shuxPh99M/JZelmLXJU6JT1mHio5oZpBOwAXheR70VkTGA8nZrGmC2B7W/hNH+1icgIEZkrInN37NhRiTBK16t1BkXFxp7LVFOmWE0p7/oD8K2IvCEibwBfAzr4X7A4IMc0TEtk6shs4qP93Dgml8Wb8m2NR7nXGeeiCsHcMFFAR+BlY0wH4CBncTnKGJNjjMkyxmSlp6cHOTRoVbsKDdMSeX/B5qBvWylVNmPMp1j5YWqgdTLG6D04LlO/eiJTR3YjMSaKG0fnsmjjXrtDUi5U1mPiPhF5DYgN8j43AhuNMbMD/bexEto2EckI7DsDm+akERGualeb3LW7wj8r7qOPWk0pb7sQuARr8L8LbI7FXRyUYzJTE5gyIpsq8dEMGDOb79fvsTsk5TJlncEJydwwgcteG0SkeeCtS4GlwPvA4MB7g7HG4LHFVe1rYwx8sDDMZ3FWrLCaUh4lIi8BvwZ+ABYDI0XkRXujchGH5ZjM1ASmjuxGtYQYBr02h3nrtMhRwSPGlH7/nohsA64xxswM+k6tp7PGADHAGmAIVrH1FlAPWA9cb4zZXdZ2srKyzNy5c4MdHgBX/usb/CK8d/v5Idm+UpFCROYZY7LCtK8lQGsTSEyBWcZ/MMa0Csf+TxXKHKNO2JJ/mP45uezYf5RxQ7vQuUGq3SGpMApVjinrDE7I5oYxxiwI3EfT1hhztTFmjzFmlzHmUmNM08CyzOIm1K5qV5uFG/N1eHGlwmsF1h85x2UCi2yKRYVJRko8U0d2o2aVOAaPncPsNbvsDkm5QFlPUXl6bpje7WojEubLVA89ZDWlvKs6sExEpovIdKzL1+ki8r6IhH4qbLdzcI6pWSWOKSOyqV01nptfz2Pm6p12h6QiXFRZXzTGbBYRT84Nk5EST+cGqby3YBO/uaQJImGYnWLDhtDvQylnc+b/vm7h8BxTo0ock4dnM2BMLkPH5THmps6c3zTN7rBUhCqzwAFvzw3Tp31t/vDvxSzdso9WtVNCv8PXXw/9PpRyoMB0DJ8Cnxhjltsdj2tFQI5JT44NFDmzGTY+j5ybsrioWfCHBFHuV55xcDw7N8wVrTOI8gnv6Zg4SoXaYGAP8GcRmS8iL4tIHxFJKs+HRWSsiGwXkcWn+bqIyPMiskpEFolIx2AGr4KrelIsbw7PplF6EsMnzOWrFbaMGqIi3BkLHDw8N0y1xBi6N6/Bv7/fRGFRceh3+MADVlPKY4wxW40x44wx/YAsYALQCfhMRKaJyO/OsIlxQM8yvt4LaBpoI4CXKx91BIqgHJOaGMPk4V1pVjOJkRPm8eWybXaHpCJMeQocT88Nc12nuuzYf5RvfgzDDW+7dllNKQ8zxhQbY2YZYx4yxpwH9AM2neEzXwNlPXnZB5hgLLlA1eMDi3pKhOWYqgkxTBqWTYuMZH49cR6f2zGFjopYZ7wHhxNzw8wI9C/E+gvIEy5pUYPUxBj+b94GLm4R4hNXOTmh3b5SDici6cBwoAEl8pMxZmglN10HKHmH7cbAe1squd3IEoE5JiUhmjeGdWXw2DncNmk+L9zYgZ6tvVebqrN3xjM4Xp8bJibKR5/2tZm2dDt77JhhXClveQ9IAaYBH5VolVXaY5ClnokO9YS+6uylxEfzxrAutMusyqg3v+ejRd6qS1XFlOcSFXh8bpjrO2VyrKiY90M9Js6991pNKe9KMMb83hjzljHmneMtCNvdiDVo4HF1gVJ/oUM9oa+tIjjHJMdFM35oFzrWq8odU74PfT5WEa88T1F5fm6YlrWr0DKjCm/P2xjaHR0+bDWlvOtDEbkiBNt9H7gp8DRVNpBvjPHeaYAIzzFJsVGMG9KFTvWrcdeU7/nP92XemqU8rjz34FzEyXPDjMcqdjzluk51eeTDpSzfuo8WtaqEZicveqpuVKo0dwIPishRoADr0pIxxpT5Sycik4HuQJqIbAQeBqKxPvwK8DFwBbAKOIQ1/533uCDHJMZGMW5IZ4aNm8tv31pAYbHhuk517Q5LOVB5LlHp3DDA1R3qEO0X3p4b4rM4SnmYMSbZGOMzxsQbY6oE+mf8i8IY098Yk2GMiTbG1DXGvGaMeSVQ3BB4emqUMaaxMaaNMUZn0IxgCTFRjL25M+c1TuO+txfyVp6zR2hW9ihPgaNzw2CNyXBJC2tMnKOFRaHZyV13WU0pjxKRC0trdsflGi7KMfExfsYMzuKCpun87p1FvDl7vd0hKYcpzyUqnRsm4Mau9flsyTY+X7KN3u1q2x2OUm50X4nXcUAXYB7WQw5KnSQu2k/OoE7cOnEeD/77B4qMYVB2fbvDUg5x2gJH54b5uQuapJGZGs+k2etCU+A8+2zwt6lUBDHG9C7ZF5FM4O82heM+LswxcdF+XhnUidsmzudP/1lMcbFh8LkN7A5LOUBZl6gqNTeMG/l8wo1d6pO7Zjerth+wOxylvGAj0NruIJSzxUb5eXlgJy5rWZOH31/Ca9+utTsk5QCnLXCCMDeMK12fVZdov4Tmeu+oUVZTyqNE5F+BSTGfF5EXgG+AhXbH5RouzjExUT5eGtCRXq1r8eiHSxn99Rq7Q1I2K9dAfxWZG8at0pJi6dk6g7fnbeBIQZBvNo6Pt5pS3jUX656becAs4PfGmIH2huQiLs8x0X4fz/fvwC/bZvD4x8t4efpqu0NSNjrjTcYhnBsmYg3oWo8PFm7mw0Vbgjv+wtNPB29bSkUgY8z4469FpBonjz6sKssDOSba7+O5vu3xi/C3T5dTVFzM7Zc0tTssZYPyPEX1HtZp4mlAiJ6PjixdG6bSOD2RSbPX6QBTSgVRYCiKq7By0wJgh4jMMMbcbWdcKrJE+X38s297/D7h6c9XUlhsuKtHM7vDUmFWngInwRjz+5BHEkFEhAFd6/PIh0tZvCmf1nVSgrPhEYFJ2iNwxl+lgiTFGLNPRG4BXjfGPCwinhtYNGQ8lGP8PuHp69vhE+HZaT9SXGz47WXNEClt3lXlRuW5BydUc8NEtOuy6pIY42fsd0G8W796dasp5V1RIpIB3AB8aHcwruOxHOP3CU9d15a+WZk8/99VPP35CgKzDikPKM8ZnArNDeN2VeKiuT4rk0mz13F/rxbUSI6r/Eb/+tfKb0OpyPYI8BnwrTEmT0QaAT/aHJN7eDDH+HzCX3/VBp9PePGr1RQWG+7v2ULP5HjAGc/gVHRuGC+4+dwGFBYbJs5aZ3coSrmCMeb/jDFtjTG3BfprsO4BVKrCfD7h8atbMyi7Pq/OWMPjHy3TMzkeUJ6nqEqdB8YY83Xww4ksDdISubRFTSbOXs9tFzchLtpfuQ0OCUxw/PrrlQ9OKfe4G3jW7iBcwcM5xucTHunTCr9PGPPtWoqM4aErW+qZHBcrzyUqnRumDEPPb8C00dt4f8FmbuhcySdaM/WJWKVKof8DBYvHc4yI8HDvlvhEGPvdWoqKDX+5qpUWOS51xgJH54YpW7dG1WlRK5mx363l+qy6lftFeeSR4AWmlHvotYRg0RyDiPCnK88hyi/kfL2GomLDo31a4/NpkeM25RrJ+BQ6N0wJIsLQ8xuyfOt+Zq7eZXc4SkUkEdkvIvtKafuBEMxsq7xMRHigVwtu7d6YSbPX8+C/f6C4WOtotynPPTj/4sRfUD6gPTo3zEmualebv3+6gldmrOa8JmkV39DAwIj0EycGJzClIoQxJtnuGDxBc8z/iAi/u7w5UT7hX/9dRWGx4W/XtsWvZ3Jcozz34Mwt8boQmGyM+S5E8USkuGg/w85vyN8+Xc4PG/NpU7eCA/81bx7cwJRSqiTNMScREe75RXP8vhODAT51fTstclyiPPfg6Nww5TAwux4vTV/FS9NX8fLAThXbyJ/+FNyglIoQIjLfGNOxsuuoM9AcU6q7ejTDL8IzX6ykyBieub4dUf6K3MGhnKQ8l6imo3PDnFFyXDQ3davPS9NXs3rHARqnJ9kdklKR5JwzTMkgQJDmRFHq535zaVN8PuGpz1ZQVGx4tm97LXIiXHkuUencMOU05LyGjPlmLa/OWM3fr2t39hvo189aTpkS3MCUcr4W5VhHJ/utLM0xZRp1cROifMJfP1lOsTE8168D0VrkRKzyFDgl54b5Q4jjiWhpSbH07ZzJ5Dnr+e1lzchIiT+7DbRvH5K4lHI6Y4wOBx4OmmPOaORFjfH7hMc+WkZR8Xz+1b8jMVFa5ESi8nzXjs8NsyqYc8OIiF9EvheRDwP9VBH5QkR+DCyrVXYfdhh+QSOKDeR8vebsP3z//VZTSqlQ0BxTLrdc0Ig/927JZ0u2cdukeRwt1JOHkag8c1GFam6YO4FlJfr3A18aY5oCXwb6ESczNYFrOtThzdnr2bbviN3hKKWUqoCbz2vIo31aMW3Zdm6dOJ8jBVrkRJqKnner1A3GIlIX+CUwpsTbfYDjT2yNB66uzD7sdMclTSkqNrz01aqz++C111pNKY8TkUQRqeTkbupnNMeclUHdGvD4Na357/LtjHxjnhY5EaaiBU5lBwl4FvgdUFzivZrGmC0AgWWNUncsMkJE5orI3B07dlQyjNCoVz2B67PqMnnOBjbvPVz+D3brZjWlPEZEfCJyo4h8JCLbgeXAFhFZIiJPiUhTu2N0Bc0xZ21A1/r87do2fP3jDoZPmKtFTgSpaIFT4TGtReRKYLsxZl6FdmxMjjEmyxiTlZ6eXtEwQm7UxU0wGF48m7M4995rNaW85yugMfAAUMsYk2mMqQFcAOQCT4rIQDsDdAXNMRXSt3M9/n5tW75dtZOh4/I4fEyLnEhw2qeoAnPAlFbICHCWjwed5DzgKhG5Amt28ioiMhHYJiIZxpgtgae2tldiH7arWy2Bfp3rMSVvPb++qDGZqQl2h6SUk/UwxhSc+qYxZjfwDvCOiESHPyylLNdnZRLlF+55ayFDxs3htcGdSYwtz4PIyi6nPYNjjEk2xlQppSUbYyr8XTXGPGCMqWuMaQD0A/5rjBkIvA8MDqw2GHivovtwilEXN0FEeOG/5TyLc9VVVlPKe5IDT1KW2gBKK4DUWdIcUynXdKjLP/u2Z87a3Qx5PY8DRwvtDkmVwUnl55PAWyIyDFgPXG9zPJVWKyWOAV3rMWHWOkZc1OjMoxtfeml4AlPKeeZhnTEu7f4+AzQKbzgupTmm0vq0r4PfJ9w5ZQE3j53D60M6kxynJxedSIwp/XaaSJgbJisry8ydO/fMK9po14GjXPTUdM5tXJ2cm7LsDkepsyYi84wxnvzhjYQco+zxyQ9b+M3k72lTN4XxQ7tQRYucCgtVjinrJuNzRGRRGe0HIC3YAblN9aRYbu3emM+XbmPO2t12h6OUI4lIgzN8XQLDSyjlCL3aZPDigI4s3pTPoDGzyT+sV1CdpqwCpwXQu4x2JXBuqAN0g6HnNaRWlTie+HgZpztjBkCvXlZTynueEpF3ROQmEWklIjVEpJ6IXCIijwLfAefYHWTE0xwTVJe3qsXLAzqxbMt+Bo6Zzd5Dx+wOSZVw2ntwdG6Y4ImP8XP3L5rxu7cX8dEPW7iybe3SV+zdO7yBKeUQxpjrRaQlMAAYCmQAh7FGO/8IeNwYo0ODV5bmmKDr0bImrw7qxMiJ87hx9Gwm3dKVaokxdoelKOMenEgQSdfHi4oNv3z+Gw4dK+KLuy8kNkoHaVWRQe/BiYwco+w1Y6U1EGCjtEQm3dKV6kmxdocUMey4B0cFkd8nPHDFOazffYgJM/XkmFKlEZFfldIuFZFSRzZXyikuapbO2MGdWbvzIDeOns3OA0ftDsnztMAJo4uapXNJixo89+WPbC9tIs4ePaymlHcNw5qjbkCgjcaa++47ERlkZ2CuoDkmpM5vmsbrN3dm3e6D9M/JZft+vapqp9MWOCKyX0T2ldL2i8i+cAbpJg9d2ZJjhcX89ZPlP/9i375WU8q7ioFzjDHXGmOuBVoCR4GuwO9tjcwNNMeE3LlN0hg3pAub9h6mX04u20r7Y1aFRUVHMq4SziDdpEFaIiMubMS/v9/088fGhw+3mlLe1cAYs61EfzvQLDBlgz6HW1maY8Iiu1F1xg3pwrb8I/TLyWVrvhY5dtBLVDYYdXET6lSN56H3FlNYVHzmDyjlHd+IyIciMlhEBmNN4fK1iCQCe+0NTany69IwlQnDurBj/1H65sxi897DdofkOVrg2CA+xs8ff3kOy7fuZ9Ls9Se+0L271ZTyrlHA60B7oAMwHhhljDlojLnYzsBcQXNMWHWqbxU5uw8co2/OLDbuOWR3SJ6iBY5NerauxQVN03j6sxUnTl/efLPVlPIoY41b8S3wX2Aa8LUp51gWItJTRFaIyCoRub+Ur3cXkXwRWRBoDwU3+gigOSbsOtarxsRbupJ/qIC+r+ayYbcWOeGiBY5NRITHrm5NQXExf/zPYmuEY00+yuNE5AZgDnAdcAMwW0SuK8fn/MCLQC+sG5P7BwYOPNU3xpj2gfZIEEOPDJpjbNEusyqTbsnmwNFC+uXksm7XQbtD8gQtcGxUv3oid1/WjGnLtvHxD1uhoMBqSnnXH4DOxpjBxpibgC7An8rxuS7AKmPMGmPMMWAK0CeEcUYmzTG2aVM3hTeHd+XQMavIWbtTi5xQ0wLHZkPPa0ibOik8/P5iCi/tAZddZndIStnJZ4zZXqK/i/LlqTrAhhL9jYH3TtVNRBaKyCci0qq0DYnICBGZKyJzd+zYUe7AI8Jll2mOsVGr2im8OTybo4XF9MuZxeodB+wOydW0wLFZlN/Hk9e2Yc+hAv6v3eVwyy12h6SUnT4Vkc9E5GYRuRlrHqqPy/E5KeW9U+/dmQ/UN8a0A/4F/Ke0DRljcowxWcaYrPT09PJHHgluuUVzjM3OyajC5OHZFBUb+uXksmr7frtDcq3ILnBWrIBx46zXBQXW0wETJ1r9Q4es/tSpVj8/3+q/+67V37nT6n/wgdXfutXqf/qp1d+wwepPm2b116yx+jNmnNh39+4wc6bVX7zY6uflWf0FC6z+ggVWPy/P6i9ebPVnzrT6K1bQqnYKj1XdScP3JjOrZjPr69OmWV/fEPij9NNPrf7WrVb/gw+s/s6dVv/dd61+fr7VnzrV6h8K3NA2caLVP356ety4k5+mGD365BFOX3rp5FmHn3sOrrrqRP/pp+Haa0/0n3wS+vU70X/0URg48ET/oYdgyJAT/QcegBEjTvTvvRdGjTrRv+suqx03apS1znEjRljbOG7IEGsfxw0caMVwXL9+VozHXXutdQzHXXWVdYzH9epl/Rsc16OH9W90XPfurvjZA6ztdu9u7Qd+/rMXRsaY+4AcoC3QDsgxxpRngL+NQGaJfl1g8ynb3meMORB4/TEQLSJpQQk8UgwcePLvpbJF81rJTB6ejTHQL2c2K7dpkRMKkV3guMi1WXVJjPLxzw8WsfvgMbvDUco2xph3jDF3G2N+a4z5dzk/lgc0FZGGIhID9MMaQ+d/RKSWiEjgdRes/LcrmLE73qFDJ/7oUbZqWjOZKSOy8Qn0z8ll+VadICDYdDZxBzl47gUs3pTP64+9zssDOxLIxUrZKhyziYvIfn5+SQmsS0+mPKOni8gVwLOAHxhrjHlcRH6NtYFXROR24FagEDgM3G2MmVnWNt2WY/531nb6dDujUCWs3WnNW3W0sIhJt2TTsrb3JgoIVY6JCvYGVcUl3nk7+5Zu49MlW3l3/iau7VTX7pCUCgtjTHIQtvExp9yvY4x5pcTrF4AXKrufiHbrrXZHoE7RMC2RqSOz6Z+Ty41jcpk4rCut66TYHZYr6CUqJ+nbl0se/g1dGqTy5/eX6KiXSqng0sk2Hal+9USmjuxGYkwUN47OZdHGvXaH5Apa4DhJfj7+/ft45oZ2GODOKQso0LmqlFLBkp9/4kEE5SiZqQlMGZFNlfhoBoyZzYINe+0OKeJpgeMkffpAnz5kpibw+DWtmbduD09/vsLuqJRSbhHIMcqZMlMTmDqyG9USYhg0Zjbz1u2xO6SIpgWOk9xxh9WAPu3rcGPXerw6Yw1fLttmc2BKKVcokWOUM9WpGs/UkdlUT4rhptdmk/fTbrtDilha4DjJr35ltYCHrmxJy4wq3PN/C9m097CNgSmlXOGUHKOcKSMlnqkju1GzShyDx85h9hpvjWYQLFrgOMnOnScG7gPiov28OKAjhUWG29+cz7FCvR9HKVUJp+QY5Vw1q8QxZUQ2tavGc/PrecxarUXO2dICx0muu85qJTRMS+Rv17bl+/V7efj9wKzjSilVEaXkGOVcNarEMXl4Npmp8QwZN4fvVmlxejZ0HBwnueeeUt/+ZdsMlmxuzEvTV3NORhVu6tYgvHEppdzhNDlGOVd6ciyTh2czYMxsho7LY/RNWVzYzGVzpIWInsFxkt69rVaKe3/RnEtb1OAvHyxlplbxSqmKKCPHKOeqnhTLm8OzaZSexC0T5vLViu12hxQRtMBxkq1bT0ymeQqfT3i2X3sapiVy25vzWb9LBwFUSp2lMnKMcrbUxBgmD+9Ks5pJjJwwT5+uLQctcJykX7+TZ+Q+RXJcNKNvyqK42DB0fB57D+mknEqps3CGHKOcrWpCDJOGZdMiI5lfT5zH50u0WC2LFjhOcv/9VitDw7REXh2UxfpdhxgxYR5HCorCFJxSKuKVI8coZ0tJiOaNYV1pVTuF2ybN59PFW+wOybG0wHGSnj2tdgbdGlfnmRvaMeen3dz91gKKivXJKqVUOZQzxyhnS4mP5o1hXWhbN4VRb37PR4u0yCmNFjhOsmGD1cqhd7va/PGX5/DxD1t59MOl+vi4UurMziLHKGdLjotmwrCudKxXlTumfM/7CzfbHZLj6GPiTjJokLWcPr1cq99yQSM27z3C2O/WkpoYwx2XNg1dbEqpyHeWOUY5W1JsFOOGdGHIuDzumvI9xcWGqzvUsTssxwh7gSMimcAEoBZQDOQYY54TkVRgKtAA+Am4wRjjrZnG/vjHs//IL89h7+Fj/OOLlcRF+xhxYeMQBKaUcoUK5BjlbImxUYwb0plh4+by27cWUFhsuK5TXbvDcgQ7zuAUAvcYY+aLSDIwT0S+AG4GvjTGPCki9wP3A7+3IT779Ohx1h/x+YS/X9uWo4XFPPHxcuKi/ToQoFKqdBXIMcr5EmKiGHtzZ4ZPmMt9by+kuNhwQ+dMu8OyXdjvwTHGbDHGzA+83g8sA+oAfYDxgdXGA1eHOzbbrVljtbMU5ffxbN/2XNayJg+9t4SpeetDEJxSKuJVMMco54uP8TNmcBYXNE3nd+8sYvIc/X/A1puMRaQB0AGYDdQ0xmwBqwgCapzmMyNEZK6IzN2xY0fYYg2LoUOtVgHRfh8v3NiBi5ql8/t3fuCN3HVBDk4pFfEqkWOU88VF+8kZ1ImLm6fzwLs/MNHj/w/YdpOxiCQB7wB3GWP2iUi5PmeMyQFyALKystz16NBf/lKpj8dG+Xl1UCduf3M+f/rPYg4dLWTkRXpPjlIqoJI5RjlfXLSfVwZ14raJ8/njfxZTVGwYfG4Du8OyhS1ncEQkGqu4mWSMeTfw9jYRyQh8PQPw3mQbF11ktUqIi/bz8sBOXNk2g79+spx/fL5CHyFXSlmCkGOU88VGWf8PXNayJg+/v4TXvl1rd0i2CHuBI9apmteAZcaYf5T40vvA4MDrwcB74Y7NditWWK2Sov0+nuvXgRuy6vL8f1fxlw+W6mCASqmg5RjlfDFRPl4a0JFerWvx6IdLGf219+69suMS1XnAIOAHEVkQeO9B4EngLREZBqwHrrchNnuNHGktgzBGhd8nPPmrtiTFRjP2u7VszT/Cs/3aExftr/S2lVIRKog5RjlftN/H8/07cNfUBTz+8TIKiw23dvfObQthL3CMMd8Cp7vh5tJwxuI4TzwR1M35fMJDvVtSp1o8j320lP6jcxlzUxbVk2KDuh+lVIQIco5Rzhft9/Fc3/b4Rfjbp8spNoZRFzexO6yw0JGMneTcc0Oy2WHnN6RO1TjunLKAX708k7E3d6ZxelJI9qWUcrAQ5RjlbFF+H//s2x6/T3jqsxUUFhnu7OH+ke91LionWbzYaiHQs3UGk0dkc+BIIVe/8B1fLtsWkv0opRwshDlGOZvfJzx9fTuu7ViXf05b6YkHULTAcZLbb7daiHSsV433bj+P+mkJDBs/l+em/Uix3nyslHeEOMcoZ/P7hKeua0vfrEye/+8qnnZ5kaOXqJzkqadCvou61RJ4+9fn8uC7P/DPaStZvDmfZ25oR5W46JDvWyllszDkGOVsPp/w11+1wecTXvxqNYXFhvt7tqC8Y9FFEi1wnKRz57DsJi7azzM3tKNN3RQe+2gZVz7/Lc/370D7zKph2b9SyiZhyjHK2Xw+4fGrWxPlE16dsYbiYsODV5zjuiJHL1E5yYIFVgsDEWHIeQ15a2Q2RcWG616eySszVuslK6XcLIw5Rjmbzyc80qcVN5/bgNHfrOWRD5e67nKVnsFxkrvuspZhHKOiU/1UPr7zAh54dxFPfrKcb3/cyTM3tKNmlbiwxaCUChMbcoxyLhHh4d4t8Ykw9ru1FBUb/nJVK9ecydECx0mefdaW3abER/PijR2ZmreBP3+whMv+MYM/X9WKazrUcc0PulIK23KMci4R4U9XnkOUX8j5eg1FxYZH+7TG54v83K8FjpO0b2/brkWEfl3q0aVhKr97exF3v7WQDxdt4Ylr2lArRc/mKOUKNuYY5VwiwgO9WuD3CS9PX01RseGJa9pEfJGj9+A4SV6e1WzUKD2JqSO78dCVLZm5eieX/XMGk+es13tzlHIDB+QY5Uwiwu8ub85vLmnClLwN/P6dRRE/h6GewXGS++6zljZfH/f7hKHnN+SSFjX4/TuLeODdH5iSt4HH+rSmTd0UW2NTSlWCQ3KMciYR4Z5fNMfvE56d9iNFxvDUde3wR+iZHC1wnOSFF+yO4CQN0hKZMiKb/yzYxOMfLeeqF79lQNd63PuL5lRNiLE7PKXU2XJYjlHOdFePZvhFeOaLlRQVG565vh1R/si74KMFjpO0bm13BD8jIlzToS6XnlOTf36xkvEzf+LDRVv4zSVNGZhdj9gonZ1cqYjhwByjnOk3lzbFF5i7qqjY8Gzf9hFX5ERWtG43c6bVHKhKXDQP927FR3dcQJs6KTz64VIufWYG7y3YpPfnKBUpHJxjlPOMurgJD/RqwYeLtnDHlO8pKCq2O6SzomdwnOTBB62lg6+Pn5NRhTeGdeXrlTt48pPl3DllATlfr+GeXzTj4uY19LFypZwsAnKMcpaRFzXG7xMe+2gZRcXz+Vf/jsRERca5ES1wnOTVV+2OoNwubJbO+U3SeG/hJp75fCVDx82ldZ0q3H5xU37RsmbEP16olCtFUI5RznHLBY2I8gl//mApt02az4sDOkTE7QmRUYZ5RfPmVosQPp91f85X93bn79e15cCRQn49cR5XPP8NHyzcTGGEnc5UkU1EeorIChFZJSL3l/J1EZHnA19fJCId7YjTVhGWY5Rz3HxeQx7t04ppy7Zx68T5HCkosjukM9ICx0lmzLBahIn2+7ghK5Npd1/Es33bU1BUzG8mf89FT00n5+vV5B8usDtE5XIi4gdeBHoBLYH+ItLylNV6AU0DbQTwcliDdIIIzTHKGQZ1a8Dj17Tmv8u3M/KNeY4vcvQSlZM8/LC1jNDr41F+H1d3qEPvdrWZtmwbr3+3lic+Xs6z037kuk51ufncBjRKT7I7TOVOXYBVxpg1ACIyBegDLC2xTh9ggrFmFMwVkaoikmGM2RL+cG0S4TlG2W9A1/pE+YT73/2B4RPmMvqmLOKinXm5SgscJxk71u4IgsLvEy5vVYvLW9ViyeZ8xn77E1PmbGDCrHVkN0qlX+d69Gxdy7G/FCoi1QE2lOhvBLqWY506wEkFjoiMwDrDQ7169YIeqK1ckmOUvfp2rodPhN+9s4ih4/J4bXBn4mOcl8/1EpWTNGpkNRdpVTuFZ25ox7f3X8x9lzdn894j3DV1AV0en8ZD7y1m8aZ8rD+olaqU0u5qP/UHqzzrYIzJMcZkGWOy0tPTgxKcY7gwxyh7XJ+VyT9uaEfuml0MGTeHQ8cK7Q7pZ/QMjpNMm2Yte/SwN44QqJEcx6iLm3DrRY3JXbuLqXkbmJJnndVpUiOJ3m1rc1X72jRMS7Q7VBWZNgKZJfp1gc0VWMfdXJxjVPhd06EuPhF+O3UBN4/NY+yQziTFOqeskEj+6zkrK8vMnTvX7jCCp3t3a+mR6+N7Dx3jg0Vb+GDhZuas3Q1Amzop9G6XQa/WGWSmJtgcoQIQkXnGmCy74yiLiEQBK4FLgU1AHnCjMWZJiXV+CdwOXIF1+ep5Y0yXsrarOUapM/tw0WbunLKADplVeX1IZ5Ljos/q86HKMVrgOMmGwO0BmZllr+dCW/IP89GiLby/cDOLNuYD0KJWMj3OqUmPljVpWydFx9axSSQUOAAicgXwLOAHxhpjHheRXwMYY14RaxTKF4CewCFgiDGmzASiOUap8vnkhy38ZvL3tKmbwvihXahyFkWOFjilcF3yUQCs23WQL5ZuY9qybeT9tIeiYkN6ciyXtqjBhc3S6daoOtUSdbLPcImUAicUNMcoVX6fLdnK7W/Op2XtFCYM7UJKfPmKHC1wSuG65PPpp9ayZ09743CQvYeOMX3FDr5Yto2vV+xg/9FCRKB17RTOa5LG+U3SyGpQTZ/ICiEtcDTHKFVe05Zu47ZJ82leK5k3hnWhasKZ/xjVAqcUrks+en28TIVFxSzcmM93q3by7aqdfL9+DwVFhpgoH+3qptCxfjWy6qfSqX41UvUMT9BogaM5Rqmz8dXy7YycOI+mNZKYOKzrGc+4a4FTCtcln61brWWtWvbGESEOHi1kzk+7mblqJ3PX7WHxpnwKiqyf50bpiWTVr0aHetVoUyeFZjWTI2aCOKfRAkdzjFJna8bKHQyfMJdGaYlMuqUr1ZNiT7uuFjilcF3yUZVypKCIRRvzmbtuN/PX7WHeuj3sOWRNExHtF5rXSqZ17RRa17Fai1rJemmrHLTA0RyjVEV8++NOho3Po0H1RCYN70raaYocLXBK4brk88EH1rJ3b3vjcAljDOt3H+KHTfks3rSPxZvyWbw5n72BokcE6qUm0LRGMs1qJtG0ZhJNayTTpEaSFj4laIGjOUapipq5aidDx+eRWS2BScO7UiM57mfraIFTCtclH70+HnLGGDbtPcziTfks27KfH7fv58dtB1i78yCFxdbvwvHCp2FaIvVTE6hfPZH61a1lZmo8sVHeKn60wNEco1Rl5K7ZxdBxeWSkxDF5eDY1qpxc5GiBUwrXJZ+dO61lWpq9cXjQscJi1u06yMptB1i5bT+rth/gp10HWbfrEAeOnhiCXARqp8RTv3oCdavFk5EST+2qcSctEx00kmcwaIGjOUapypqzdjdDXp9DzSpxvDk8m1opJ4ocLXBK4brkoxzHGMPug8f4adch1u8+yE87D7F+9yF+2nWQTXsOs+PAUU79FaoSF0XtqvHUrhpPzSpxpCfFkJYcS3pSLGnJsaQlxZKWFENSbBTW2HPOpgWO5hilgmHeut0MHptH9aQYJg/PpnbVeCB0OcZdf2pGunfftZa/+pW9caj/ERGqJ8VSPSmWTvWr/ezrxwqL2bbvCJv3HmZL/hE25x9my94jbMk/zOa9R1i0cS+7Dh77WREEEBftCxQ7VsFTNSGGqvHRVE2IJiUhhmoJ0VSNj7H68dFUS4whMcYfEUWRcijNMcpGneqnMmFYFwa/Noe+ObOYPDybutVCNyWP487giEhP4Dms4dbHGGOePN26rvvrSq+Pu1JhUTG7Dx1j5/5j7DxwlB37j7LzwNESr6339x0uYO/hAg4dKzrttqJ8QtWEaJLjokmM9ZMUG0VSbDTJcVGBvvXaej+KxNio//UTYvzERfuJj/ETH2299pdj+gs9g6M5RqlgWrhhL4Nem01yXDRTRmRTr3qi+8/giIgfeBG4DGvm3zwRed8Ys9TeyMLkvffsjkCFQJTfR43kuFKfHijNkYKi/xU7ew8VsPfQMWt5+PiygANHCjlwtJADRwrZvPew9TrQP1ZUXO7YYqJ8xEdbBU/88QIo2ndSEaRcRHOMcoB2mVWZdEs2A1+bTb+c3JDtx1EFDtAFWGWMWQMgIlOAPoA3CpyUFLsjUA4QFygsTn3SoLyOFhZx8GgRB44Usv9oAQePFrH/SAGHC4o4fKyIIwVFgdfFgWWhtSwoPunrew5an1EuojlGOUSbuim8Obwrw8eH7gyp0wqcOsCGEv2NQFebYgm/qVOtZd++9sahIlpslJ/YKH/QpquQ+4KyGeUEmmOUg7SqncJX93Un7sHQbN9pBU5pNwScdJOQiIwARgDUq1cvHDGFz8svW0tNPkqpUNAcoxwmlOOKOa3A2QhklujXBTaXXMEYkwPkgHUDYPhCC4OPP7Y7AqWUm2mOUR7itAInD2gqIg2BTUA/4EZ7QwqjhNA9LqeUUppjlJc4qsAxxhSKyO3AZ1iPiY81xiyxOazwmTjRWg4caG8cSil30hyjPMRRBQ6AMeZjwJvnUceMsZaafJRSoaA5RnmI4wocT/viC7sjUEq5meYY5SFa4DhJdLTdESil3ExzjPIQn90BqBLGjbOaUkqFguYY5SFa4DiJJh+lVChpjlEe4rjJNs+GiOwHVtgdR4ikATvtDiIE9LgiT3NjTLLdQdhBc0xE0uOKPCHJMZF+D84Kt85yLCJz3XhselyRR0RcNJ32WdMcE2H0uCJPqHKMXqJSSimllOtogaOUUkop14n0AifH7gBCyK3HpscVedx8bGfi5mN367HpcUWekBxbRN9krJRSSilVmkg/g6OUUkop9TNa4CillFLKdSK2wBGRniKyQkRWicj9dsdTHiLyk4j8ICILjj8WJyKpIvKFiPwYWFYrsf4DgeNbISKXl3i/U2A7q0TkeRGRMB/HWBHZLiKLS7wXtOMQkVgRmRp4f7aINLDxuP4sIpsC37MFInJFBB5Xpoh8JSLLRGSJiNwZeD/iv2ehpDlGc0wYjy2i84xjc4wxJuIa4AdWA42AGGAh0NLuuMoR909A2inv/R24P/D6fuBvgdctA8cVCzQMHK8/8LU5QDdAgE+AXmE+jguBjsDiUBwHcBvwSuB1P2Cqjcf1Z+DeUtaNpOPKADoGXicDKwPxR/z3LIT/ZppjNMeE89giOs84NcfY/gtZwX/MbsBnJfoPAA/YHVc54i4t+awAMkr8kKwo7ZiAzwLHnQEsL/F+f+BVG46lwSm/oEE7juPrBF5HYY3eKTYd1+kST0Qd1ymxvwdc5pbvWYj+jTTHlPJ9DvOxuDLHnObYXJVnnJJjIvUSVR1gQ4n+xsB7TmeAz0VknoiMCLxX0xizBSCwrBF4/3THWCfw+tT37RbM4/jfZ4wxhUA+UD1kkZ/Z7SKyKHBq+fgp1og8rsBp3Q7AbNz9PasszTE/f99ubv95dUWecVKOidQCp7TrwZHwvPt5xpiOQC9glIhcWMa6pzvGSDv2ihyHk47xZaAx0B7YAjwTeD/ijktEkoB3gLuMMfvKWrWU9xx9bCEQqcejOcYSaT+vrsgzTssxkVrgbAQyS/TrApttiqXcjDGbA8vtwL+BLsA2EckACCy3B1Y/3TFuDLw+9X27BfM4/vcZEYkCUoDdIYu8DMaYbcaYImNMMTAa63t2UowBjj4uEYnGSjyTjDHvBt525fcsSDTH/Px9u7n259UNecaJOSZSC5w8oKmINBSRGKwbjt63OaYyiUiiiCQffw38AliMFffgwGqDsa5dEni/X+DO8YZAU2BO4DTffhHJDtxdflOJz9gpmMdRclvXAf81gQuv4Xb8lzPgGqzvGUTQcQXieA1YZoz5R4kvufJ7FiSaYzTHhE2k5xnH5phw3ngU5JuYrsC6U3s18Ae74ylHvI2w7hpfCCw5HjPWNcQvgR8Dy9QSn/lD4PhWUOIpBiAL6xdgNfAC4b+BbDLWadQCrKp6WDCPA4gD/g9YhXVHfSMbj+sN4AdgUeAXLCMCj+t8rFO5i4AFgXaFG75nIf530xyjOSZcxxbRecapOUanalBKKaWU60TqJSqllFJKqdPSAkcppZRSrqMFjlJKKaVcRwscpZRSSrmOFjhKKaWUch0tcFSFiEiRWLPeLhaRD0SkainrNBCRwyKyoIL7uEBElkqJWXeVUt6gOUZVlhY4qqIOG2PaG2NaY40mOeo06602xrSvyA6MMd9gjaWglPIezTGqUrTAUcEwi3JMxhf4a2u5iIwPTCr3togkBL7WWURmishCEZlzfERWpZRCc4yqAC1wVKWIiB+4lPIPY98cyDHGtAX2AbcFhsKfCtxpjGkH9AAOhyJepVRk0RyjKkoLHFVR8YHr3ruAVOCLcn5ugzHmu8DriVhDfDcHthhj8gCMMfuMMYVBjlcpFVk0x6hK0QJHVdThwHXv+kAMp78+fqpT5wYxgJTyvlLK2zTHqErRAkdVijEmH7gDuFdEosvxkXoi0i3wuj/wLbAcqC0inQFEJFlEokISsFIqomiOURWlBY6qNGPM91gzGPcrx+rLgMEisgjrtPPLxphjQF/gXyKyEOtUdFyo4lVKRRbNMaoitIJVFWKMSTql37ucHy02xvy6lO3lAdnBiE0pFfk0x6jK0jM4KpSKgJTKDMIFfADsDGZQSinX0ByjTkuM0fuulFJKKeUuegZHKaWUUq6jBY5SSimlXEcLHKWUUkq5jhY4SimllHIdLXCUUkop5Tr/D+ARe9kOnNafAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "R=np.arange(0,20000,10)\n", "I=I0*np.exp(-R/h)\n", "\n", "Fig,(LinPlot,LogPlot)=plt.subplots(1,2,figsize=(8,4))\n", "LinPlot.plot(R,I)\n", "LinPlot.set_ylabel('I [Lsun/pc^2]')\n", "LinPlot.axhline(y=Iwant,ls=':',color='red')\n", "\n", "LogPlot.plot(R,np.log10(I))\n", "LogPlot.set_ylabel('log(I) [Lsun/pc^2]')\n", "LogPlot.axhline(y=np.log10(Iwant),ls=':',color='red')\n", "\n", "for panel in [LinPlot,LogPlot]:\n", " panel.set_xlim(0,20000)\n", " panel.set_xlabel('R [pc]')\n", " panel.axvline(x=Rwant,ls=':',color='red')\n", "\n", "Fig.tight_layout() " ] }, { "cell_type": "markdown", "id": "903015f5", "metadata": {}, "source": [ "Yay!" ] }, { "cell_type": "code", "execution_count": null, "id": "06c4687d", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.12" } }, "nbformat": 4, "nbformat_minor": 5 }