tutorial/3-Lane Emden Numerical Solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our goal is to simulate a isothermal sphere undergoing gravitational collapse. We can write down the force balance equation for the sphere, which describes the competition between gravitational infall versus thermal pressure support.\n",
    "$$ -\\frac{1}{\\rho}\\nabla P -\\nabla \\Phi_g =0$$\n",
    "Using Poisson's equation ($\\nabla \\Phi_g^2 =4\\pi G \\rho$) and the isothermal equation of state ($P = \\rho a_T^2$), we can rearrange the equation to the form of the isothermal Lane-Emden equation: \n",
    "\n",
    "\n",
    "$$\\frac{1}{\\psi^2} \\frac{d}{d\\psi}\\Bigg(\\xi^2 \\frac{d\\psi}{d\\xi}\\Bigg)=e^{-\\psi}$$\n",
    "Using product rule, we can rewrite the Lane-Emden equation as: \n",
    "\\begin{align}\n",
    "\\frac{d^2\\psi}{d\\xi^2}+\\frac{2}{\\xi}\\frac{d\\psi}{d\\xi}=e^{-\\psi}\n",
    "\\end{align}\n",
    "with boundary conditions as: $\\psi(0)=0$ and $\\frac{d\\psi}{d\\xi}(0)=0$.\n",
    "\n",
    "The Lane-Emden equation describes a family of solution based on the single parameter $\\xi_{max}$, which specifies the cloud radius and determine the density contrast $\\rho_c/\\rho_0$ (ratio between central density and edge of the sphere). Analytically, the sphere is marginally stable (i.e. just at the verge of collpasing) when $\\rho_c/\\rho_0$=14.1 [1]. This marginally-stable state is called the Bonner-Ebert sphere. Numerically, we can see that a cloud radius of  $\\xi$=6.451 corresponds to a density contrast of $\\rho_c/\\rho_0$=14.038. The singular isothermal sphere case corresponds to $\\xi_{max}=\\infty$. [1]\n",
    "\n",
    "We want to solve the Lane-Emden equation to figure out what is the density distribution that yields a marginally-stable sphere, but the Lane-Emden solution can not be solved in closed form, so instead we need to do numerical intergration to get the solution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start integrating at a small poisitive number (1e-6) to avoid div-by-0 error and integrate up to xi_max, with interval of $\\Delta \\xi$ = 0.01."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xi_max = 6.451\n",
    "from scipy import integrate\n",
    "def solvr(Y, t):\n",
    "    return [Y[1], exp(-Y[0])-2/t*Y[1]]\n",
    "xi = np.arange(1e-6, xi_max, 0.01)\n",
    "IC = [0, 0] # initial conditions for psi and psi prime\n",
    "psi = integrate.odeint(solvr, IC, xi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a spherically symmetric cloud, $\\rho(r)=\\rho_c e^{-\\psi}$, so we can use the numerically integrated results for $\\psi$ to compute the density profile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f3b6692ba50>"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEXCAYAAABbKnTjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8lNXZ//HPxb6vsu8iCiKCLG6goFRB24q1WgFbFTew\ntlZrXVr9VVqtbZ+nj60WNcUFoeKGVAVFWZQUBWQPexCRJYBElpAAAbKd3x9nBMRAJslk7lm+79dr\nXrPdM3NNxG9Ozn3f1zHnHCIiEp8qBV2AiIiUnUJcRCSOKcRFROKYQlxEJI4pxEVE4phCXEQkjpUY\n4mZW3cwWmNkyM1tpZo+eYLunzWy9maWZWY/IlyoiIserUtIGzrnDZnaJcy7XzCoDc83sA+fcwm+2\nMbMrgI7OuU5mdh6QApxfcWWLiAiEOZ3inMsN3ayOD/7jzxAaAkwIbbsAqG9mzSJVpIiIFC+sEDez\nSma2DNgBzHTOLTpuk1ZAxjH3t4UeExGRChTuSLzIOXcO0Bo4z8zOrNiyREQkHCXOiR/LOZdjZrOB\nwcCaY57aBrQ55n7r0GPfYmZq1CIiUgbOOSvu8XCOTjnFzOqHbtcELgPSj9tsCnBjaJvzgb3OucwT\nFJIwl0cffTTwGvR99H3i9aLvE/7lZMIZibcAxptZJXzov+Gcm2ZmI30mu7Gh+1ea2RfAAWBEGO8r\nIiLlFM4hhiuBnsU8/q/j7v8ignWJiEgYdMZmOQwYMCDoEiJK3ye26fvEtqC+j5U03xLRDzNz0fw8\nEZFEYGa4E+zYLNXRKRWpffv2bN68OegyYla7du3YtGlT0GWISIyJmZF46DdN1GqJN/r5iCSvk43E\nNScuIhLHFOIiInFMIS4iEscU4iIicUwhLiISxxTiIiJxTCEeAYsWLWLDhg1BlyEiSUghHgHp6el0\n7Ngx6DJEJAkpxEVE4phCvJT27t3LmDFj+MlPfkJBQQEFBQVUrVqVwsJCOnTowK5du4IuUUSSSEyH\nuFn0LuGaOnUqo0aNYvny5eTl5bFw4UL69OlD5cqV6devH7m5uSW/iYhIhMR0iDsXvUu4rr76ahYs\nWEC3bt2oVasW69evPzIf3rdvX9q2bVtBPw0Rke+K6RCPRXXr1uXDDz9kyJAh33o8Ly+PqlWrBlSV\niCQrhXgZHD58mLp165KXl0e1atUoKioiJSWFYcOGBV2aiCQZtaItg927d/PYY49Rr149GjduTO3a\ntbn22mtp0KBBhX1mPP18RCSyTtaKViFeDhMmTODGG2+MymfF489HRCJD/cQrSFFRUdAliEiSU4iX\nQ+XKlYMuQUQSXEbGyZ/XdEqc0M9HJLlkZ8Nf/gJjx8KePZpOERGJC3l58PTTcPrp8PXXsHz5ybeP\nmdXuRUSSmXPw1lvw29/6AJ81C7p1K/l1CnERkYB9+in85jd+FP6vf8HAgeG/ViEuIhKQdevgwQdh\n2TL4059g+HCoVMpJbs2Ji4hEWWYm/Pzn0K8f9O3rw/ynPy19gEMYIW5mrc3sYzNbbWYrzezuYrbp\nb2Z7zWxp6PJI6UsREUlsubnw2GPQtSvUqAHp6XD//f52WYUznVIA/No5l2ZmdYAlZjbDOZd+3HZz\nnHNXlb0UEZHEVFQEEyfC734HF14ICxfCqadG5r1LDHHn3A5gR+j2fjNbC7QCjg/xUnTlFhFJDnPn\nwr33+tuvv+6nTyKpVDMwZtYe6AEsKObpC8wszczeN7MzI1CbiEjc2rQJrr8ehg6Fu++Gzz6LfIBD\nKUI8NJXyFvAr59z+455eArR1zvUAxgDvRK7E+DVv3rygSxCRKMvJ8cd69+rl577Ls9MyHGEdYmhm\nVfAB/m/n3LvHP39sqDvnPjCzZ82skXNuz/Hbjh49+sjtAQMGMGDAgDKUHfv27t3L7NmzufDCC4Mu\nRUSioLAQXnoJfv97GDwYVqyAVq3K9l6pqamkpqaGtW1YvVPMbAKwyzn36xM838w5lxm6fS7wpnOu\nfTHbJU3vlBdffJG+ffvSuXPniLxfov18RBLJRx/Br38NDRrA3/8OPXtG9v1P1oq2xJG4mfUFbgBW\nmtkywAG/A9oBzjk3FrjWzO4E8oGDwPWRKj5erVmzhltvvTXoMkSkAn3+uT/TcvVq+N//hR/9qHQL\nr0dCOEenzAVO2nPVOfcM8Eykiop3mZmZtGzZMugyRKSCZGXBH/7gDxt88EGYNAmqVw+mFp2xWUaF\nhYU8+eSTjB8/nlGjRpGdnX3kuUmTJnHttdeW6jUiEvsKC31r2C5d4NAhWLPGj8SDCnCI8d4p9ofo\n/V3iHi3dfPODDz7IFVdcwcCBA9m3bx9jxozh4YcfBmDjxo20a9euVK8Rkdg2dy788pdQpw58+CH0\n6BF0RV5Mh3hpgzVavvrqKz799FP+9re/AZCdnc2OHTsA2LRpE6cWcyrWyV4jIrFr2zZ44AGYM8fP\ne19/ffTnvU9G0yllMG3aNAYPHnzk/uzZs7nssssAP5Vy3XXXhfWa733vexVfrIiUyaFD8Oc/Q/fu\n0KGD73MydGhsBTgoxMvkk08+oUuXLgCkpaVRWFjIVVf5tjGZmZk0bdo0rNcMGTKE8ePHM2vWLB54\n4IHofQEROSHnYMoUOOss3+Nk4UJ4/HGoXTvoyooX09MpsWr+/Pn07duXl19+mTVr1jB16lQA1q5d\nS9euXcN+zeTJkzlw4ADf+973NCoXiQHp6XDPPbB5Mzz7LFx+edAVlUwhXkrbt2+nfv363H777d95\nbvLkyfzyl78M+zWpqancc889AOzfv586depUTNEiclLZ2fDHP8KECb7T4C9+AVWrBl1VeDSdUkqf\nfPIJfU/QxSY7O5v69euH/Zphw4bx6aefMmPGDBYuXBjxWkXk5IqKYNw46NzZB/mqVb7jYLwEOGgk\nXmrp6elcXszfWKtXr+aCCy4o1WsuvPBC9VYRCciyZXDXXT7Ip0yBPn2CrqhswuqdErEPS6LeKZGm\nn49IZOzdC//v/8Gbb8ITT8CIERXXYTBSTtY7JcZLFxGJDOf8nHeXLn5V+TVr4NZbYz/AS6LpFBFJ\neCtX+oWJDx6Ed9+Fc88NuqLIifPfQSIiJ5aT41vEDhwIN9wACxYkVoCDQlxEEpBz8NprfuokO9u3\nih01CiqftB9rfNJ0iogklDVr/HHee/b4FrGJfgCYRuIikhD27/e9vfv394szLF6c+AEOMTQSb9eu\nHRZrnWViSHGtbUXEe+cdv6L8gAF+J2bz5kFXFD0xc5y4iEhpbdnie3yvWwcpKT7EE5GOExeRhFJQ\nAE8+6Rck7t0bli9P3AAvScxMp4iIhGPRIrjjDmjcGObPh06dgq4oWBqJi0hcyM72UydXXeXXtZw5\nUwEOCnERiXHOwVtvQdeucPiwP+b7hhtib4WdoGg6RURi1qZN/pjvjRvh9dehX7+gK4o9GomLSMzJ\nz/eLEvfuDX37+raxCvDiaSQuIjHls89g5Eho0cL3OunYMeiKYptCXERiwr59fmm0yZP94YPXX695\n73BoOkVEAjdtml9dPjfXL5E2dKgCPFwaiYtIYHbu9KvLL1gAL73kW8ZK6WgkLiJR5xy88gp06wYt\nW8KKFQrwsipxJG5mrYEJQDOgCHjeOfd0Mds9DVwBHABuds6lRbhWEUkAmzf7HZc7dsD770OvXkFX\nFN/CGYkXAL92znUFLgDuMrPOx25gZlcAHZ1znYCRQErEKxWRuFZYCE895UN7wAB/+rwCvPxKHIk7\n53YAO0K395vZWqAVkH7MZkPwo3WccwvMrL6ZNXPOZVZAzSISZ1atgttug+rVYd48OP30oCtKHKWa\nEzez9kAPYMFxT7UCMo65vy30mIgkscOH4fe/h0sugVtugdmzFeCRFvbRKWZWB3gL+JVzbn9ZP3D0\n6NFHbg8YMIABydo/UiTBzZ3rR9+dO/tWsS1bBl1R/EhNTSU1NTWsbcNaFMLMqgDvAR84554q5vkU\nYLZz7o3Q/XSg//HTKVoUQiTxHTjgT9qZNAmefhp+/GMd811ekVgU4iVgTXEBHjIFuDH0YecDezUf\nLpJ8UlPh7LMhK8vPg197rQK8opU4EjezvsAcYCXgQpffAe0A55wbG9puDDAYf4jhCOfc0mLeSyNx\nkQT0zSLF777rl0n7wQ+CriixnGwkrjU2RaRcPvrIz31fconvedKgQdAVJZ6ThbhOuxeRMsnJgfvv\nhw8+gLFjYfDgoCtKTjrtXkRKbfp0f8q8c7BypQI8SBqJi0jY9u6F++7zUygvvACXXRZ0RaKRuIiE\n5f33/ei7enU/+laAxwaNxEXkpLKyfLvYTz6BCRP8DkyJHRqJi8gJfTP6rl/ft4tVgMcejcRF5Dty\ncuDee32vk4kToX//oCuSE9FIXES+5eOP/VmXVar4nicK8NimkbiIAH59y4cegrffhuef12GD8UIj\ncRFh/nzo0cPvxFyxQgEeTzQSF0lihw/D6NHw8svwzDNwzTVBVySlpRAXSVLLlsGNN0KnTn7uu2nT\noCuSstB0ikiSyc+Hxx6DQYN858HJkxXg8UwjcZEksmYN3HQTNG4MS5dC69ZBVyTlpZG4SBIoLPRt\nYvv3h9tv950HFeCJQSNxkQS3aZMffQMsWACnnhpoORJhGomLJCjnfK+Tc8+FH/7Qn32pAE88GomL\nJKDdu2HkSFi3DmbN8mdgSmLSSFwkwUyfDt27Q7t2sGiRAjzRaSQukiByc48uVjxhAlx6adAVSTRo\nJC6SAJYuhV69/DTK8uUK8GSiEBeJY4WF8MQTvtfJ738Pr74KDRsGXZVEk6ZTROLUl1/60+arVYMl\nS6BNm6ArkiBoJC4SZ5yDcePgvPN8w6pZsxTgyUwjcZE4snOnP3Rwwwa/eEO3bkFXJEHTSFwkTkyb\n5g8dPO00WLhQAS6eRuIiMe7gQXjgAZgyxe+4HDAg6IoklmgkLhLDVq6EPn3g66/9oYMKcDleiSFu\nZi+aWaaZrTjB8/3NbK+ZLQ1dHol8mSLJxTn45z/98d6/+Q28/jo0aBB0VRKLwplOGQf8E5hwkm3m\nOOeuikxJIsnt669hxAi/E3P+fD8HLnIiJY7EnXOfAlklbGaRKUckuX34oV+wuHt3mDtXAS4li9SO\nzQvMLA3YBtzvnFsTofcVSQqHDsFvf+uXStPOSymNSIT4EqCtcy7XzK4A3gFOP9HGo0ePPnJ7wIAB\nDNC/Vklya9bAsGF+weK0NGjUKOiKJGipqamkpqaGta0550reyKwdMNU5V2JTSzPbCPRyzu0p5jkX\nzueJJAPnICXF9zz5y1/gllvANDEpxTAznHPF/usIdyRunGDe28yaOecyQ7fPxf9i+E6Ai8hRO3fC\nbbfB1q3w6adwxhlBVyTxKpxDDF8F5gGnm9kWMxthZiPN7I7QJtea2SozWwb8A7i+AusViXszZ/qd\nl507+6NPFOBSHmFNp0TswzSdIkns8GF4+GF/zPf48TBwYNAVSbyIxHSKiJRDejoMH+6XTFu+HBo3\nDroiSRQ67V6kAjkHL74IF10Eo0bBf/6jAJfI0khcpILk5Pi2satXw3//C2eeGXRFkog0EhepAIsX\nQ8+evt/JggUKcKk4CnGRCHIO/v53uPJK+POf4bnnoGbNoKuSRKbpFJEI2bULbr7ZXy9YAB06BF2R\nJAONxEUi4L//hXPOga5d4ZNPFOASPRqJi5RDYSE8/rg/fX7cOBg8OOiKJNkoxEXKaNs2uOEGqFwZ\nliyBli2DrkiSkaZTRMrg/fehVy/43vdgxgwFuARHI3GRUsjL832/J03yl4suCroiSXYKcZEwbdgA\nQ4f6UfeyZTrzUmKDplNEwvDGG3DBBfCzn8E77yjAJXZoJC5yErm5cM89MHu2X/+yZ8+gKxL5No3E\nRU5g9Wo491wf5EuXKsAlNinERY7jHDz/vF+s+L774N//hrp1g65KpHiaThE5Rna27zy4Zg3MmQNd\nugRdkcjJaSQuErJokZ8yadTI9z5RgEs80Ehckl5Rke88+Ne/wrPPwrXXBl2RSPgU4pLUdu70nQd3\n74aFC6F9+6ArEikdTadI0kpN9Z0Hu3XznQcV4BKPNBKXpFNYCH/8I4wdCy+/DIMGBV2RSNkpxCWp\nbN3qOw9WqeKP/W7RIuiKRMpH0ymSNN57D3r3hssv950HFeCSCDQSl4SXlwcPPQRvveUv/foFXZFI\n5CjEJaF98YXvPNiqFaSl+WPARRKJplMkYb32mu88eNNNvvOgAlwSkUbiknAOHIBf/cqfNj99uhpX\nSWIrcSRuZi+aWaaZrTjJNk+b2XozSzOzHpEtUSR8q1ZBnz5w6JBf91IBLokunOmUccAJj6Q1syuA\njs65TsBIICVCtYmEzTl/3Pcll8ADD6jzoCSPEqdTnHOfmlm7k2wyBJgQ2naBmdU3s2bOucxIFSly\nMtnZcPvtkJ6uzoOSfCKxY7MVkHHM/W2hx4q1O3c3OYdzOFRwiMKiwgh8vCSzhQv9qfOnnKLOg5Kc\nor5js/WQ1hQWFVLkiihsW0iNTjVoWKMhjWo2omHNhjSs0ZCmtZvSul5r2tRrQ5v6bWhTrw3tGrSj\nVtVa0S5XYlRRETz5JPzP/8Bzz8GPfxx0RSKRk5qaSmpqaljbmnOu5I38dMpU59zZxTyXAsx2zr0R\nup8O9C9uOsXM3LGf55wjNz+XrENZZB3MOnKdeSCTjOwMMnIy2JqzlYycDDKyM2hSuwmdT+lMl1O6\nHLnu0bwH9WvUD+vLSmLYudMfNpiV5Q8jVOMqSXRmhnPOinsu3JG4hS7FmQLcBbxhZucDe8OdDzcz\nalerTe1qtWldr/VJty0sKmRz9mbSd6Wzdudalmxfwvjl41mZuZKWdVvSq2UverXoRe+WvenTsg+1\nq9UO86tJPJk92684/9OfwmOPQdWqQVckEqwSR+Jm9iowAGgMZAKPAtUA55wbG9pmDDAYOACMcM4t\nPcF7uXBG/qVRUFRA+q50lmxfwpKvlrB4+2KWZy6nW9NuXNzuYi5qexH92vajYc2GEf1cia6CAt95\n8IUXYNw4dR6U5HKykXhY0ykRLCTiIV6cg/kHWbBtAZ9s/oQ5W+bw2dbP6NiwI5d3vJzBpw2mb5u+\nVK9SvcLrkMjYssV3HqxRwx862Lx50BWJRFfShfjx8gvzWbR9EdO/mM70DdNZu2stF7e7mMEdB3Nl\npyvp0LBD1GuS8Lzzjl+4+N57/fHfldQoQpJQ0of48Xbn7mbWl7P4cMOHTFs/jRZ1WnB156v5Uecf\ncXazszE70fS/RMuhQ3D//TB1Krz6Klx4YdAViQRHIX4ShUWFzMuYxzvp7/B2+tsAXN35aq7pcg0X\ntrmQSqahX7StWwfXXw+nnQbPPw8NtTtDkpxCPEzOOVZ+vZK3177NpDWTyDmcw9CzhjK823C6N+uu\nEXoUjB8P990Hjz/up1H0IxdRiJfZysyVvLbqNV5b9Ro1qtRg+FnDGdZtGKc1Oi3o0hLOvn3w85/7\nplVvvOEXLxYRTyFeTs45FmxbwKsrX+XN1W/Spn4bbup+E8O7DadRTTWpLq+lS/30Sf/+8NRTUFuH\n+It8i0I8ggqKCvh448eMSxvHB+s/4PKOlzOixwgu73g5lStVDrq8uOIcPP00/OlP/nro0KArEolN\nCvEKknUwi9dXvc64tHFs27eNG8++kRHnjOD0xqcHXVrM27ULbrkFvvoKXn8dOnYMuiKR2KUQj4JV\nX6/i5bSXeWXFK3Rs1JHbzrmN68+6Xk27ijFnjj955/rr4YknoFq1oCsSiW0K8SjKL8xn2vppPL/0\neeZvnc8N3W5gZK+RdG3aNejSAldY6I86SUmBl16CK64IuiKR+KAQD8jmvZt5YekLvLjsRU5teCoj\ne43k2jOvpWbVmkGXFnUZGb5pVeXK8Mor0LJl0BWJxA+FeMDyC/N57/P3SFmSwpLtS7ix+43c0esO\nOp/SOejSomLSJLjrLvj1r/1ZmJW1/1ekVBTiMeTLrC95fsnzjEsbR+dTOjOy10iu6XJNQjbk2r//\n6Krzr77qFzAWkdJTiMegvMI83k1/l5QlKaz6ehUjeoxgZK+RCdOMa8kSGDYM+vb1hw9q0WKRslOI\nx7h1u9YxdslYJqyYQJ+WfRjVexTf7/T9uDzuvKgI/vY3f/nnP/0RKCJSPgrxOHEw/yBvrn6TlCUp\nbMvZxh297uDWc26lRd0WQZcWlm3b/LJphw/7nZft2gVdkUhiOFmIq0VfDKlZtSY39biJ+bfO592h\n75KRncGZz57JdZOu46MvPyKWfwG++y706uVPnZ89WwEuEi0aice4nMM5vLLiFZ5b/Bx5hXmM6jWK\nm3rcFDM9W3JzfdfB6dNh4kS44IKgKxJJPJpOSQDOOeZlzOO5xc/x3ufvcXXnqxnVexTntTovsBa5\naWkwfDj07AnPPAP16wdShkjCU4gnmJ0HdvJy2sukLEmhXvV63Nn7ToZ3G06danWi8vlFRb7b4BNP\nwN//7k/iEZGKoxBPUEWuiJkbZpKyJIX/bvovw84axqjeo+jWrOKacW/ZAjff7Hde/vvfcOqpFfZR\nIhKiHZsJqpJVYtBpg3j7+rdZcecKmtRuwuCJg+n3Uj8mrpjIoYJDEfss5/wRJ717w2WX+RN4FOAi\nwdNIPMHkF+Yz9fOppCxOIW1HGjf3uJmRvUbSsVHZe73u2QOjRsHq1T7IzzknggWLSIk0Ek8iVStX\n5Zou1zDjZzOYe8tcilwR5794PoNeGcQ76e9QUFRQqvebPh3OPhtat/ZnYSrARWKLRuJJ4FDBISat\nnkTKkhQ2793M7T1v57aet9GqXqsTviY3Fx580B//PW4cDBwYxYJF5Fu0Y1OOWL5jOSmLU3h99etc\n2uFSRvUaxcBTB1LJjv5RtmgR/Oxnfv57zBho0CDAgkVEIS7fte/wPiaunMhzi58jNz+Xkb1G8tOz\nbmbsP07hmWd80yr1PRGJDQpxOSHnHJ9t/Yy/fPQc738xhcY5A/nr0Ju54bzBVK1cNejyRIQI7Ng0\ns8Fmlm5mn5vZg8U839/M9prZ0tDlkfIWLdFRWGjMefUC5t0/gf9rvZnHbhrM2LV/ps3f2/CbGb9h\n1dergi5RRE6ixJG4mVUCPgcGAtuBRcBQ51z6Mdv0B+5zzl1VwntpJB5D1q71J+7UrQsvvADt2x99\nbt2udYxfPp4JyyfQom4Lbu5+M8O6DYuZni0iyaS8I/FzgfXOuc3OuXzgdWBIcZ9TjholigoL4X/+\nBy6+GEaMgJkzvx3gAGeccgZPDHyCzfds5vFLHueTLZ/Q4akO/GTST5i2fhr5hfmB1C4i3xbOSPzH\nwCDn3B2h+z8FznXO3X3MNv2BycBWYBtwv3NuTTHvpZF4wNLT/ei7Vi148UXoUIqFhLIOZvH6qtcZ\nv3w8X2Z9yU+6/oTh3YZzQesLAmvCJZIMonGyzxKgrXOuBzAGeCdC7ysRUlgI//u/0K8f3HgjzJpV\nugAHaFizIXf2uZPPbvuM+bfOp3md5tw65VZOffpUHv7oYVZ/vbpiiheRE6oSxjbbgLbH3G8deuwI\n59z+Y25/YGbPmlkj59ye499s9OjRR24PGDCAAQMGlLJkKa3Vq+G226BGDVi4MDI9Tzo26sgjFz/C\nwxc9zPLM5UxcMZFBrwzilFqnMLzbcIaeNZS29duW/EYi8h2pqamkpqaGtW040ymVgXX4HZtfAQuB\nYc65tcds08w5lxm6fS7wpnOufTHvpemUKDp82LeLffZZ+OMfYeRIqFSBjRaKXBGfbP6EV1e+yuS1\nk+nSpAvXnXkd13S5htb1WlfcB4skuHIfJ25mg4Gn8NMvLzrn/mJmIwHnnBtrZncBdwL5wEHgXufc\ngmLeRyEeJXPn+tH3GWf4BRtanfgM+wqRV5jHjA0zmLx2MlPWTeH0xqdzbZdr+fGZP6Z9g/bRLUYk\nzulknySSkwMPPeR7njz9NFxzDQS9zzGvMI/ZG2czee1k3k5/m/YN2h8J9NManRZscSJxQCGeJKZM\ngbvugiuugL/+FRo2DLqi7yooKmDO5jm8teYt/rP2PzSv05whZwzhh2f8kJ4ten6rh4uIeArxBLd9\nO9xzj1/zcuxYiJd9xYVFhczLmMeUdVOY+vlU9uXt4wedfsBVZ1zFpR0upWbVmkGXKBITFOIJqqDA\ndxl8/HG/0/KRR6BmHOfe57s/Z+q6qUz5fArLvlrGpR0u5Yen/5AfnP4DmtVpFnR5IoFRiCegefPg\nzjuhSRMf5J07B11RZO3O3c0HX3zAlHVTmLFhBp0ad2JQx0EM6jiI81ufr+ZcklQU4glk1y6/WMOH\nH8L//Z9vFxv0jsuKlleYx7yMeUz/YjrTN0zny6wvubTDpT7UTxuko10k4SnEE0BRkW9S9cgjcMMN\n8Ic/QL16QVcVjMz9mcz8cibTN0xnxoYZNKjRgEEdB3F5x8u5qO1F1K9RP+gSRSJKIR7nFiyAu++G\nypX9iTs9egRdUewockUs37Gc6RumM/PLmSzctpAup3ThkvaXcEmHS+jXth91qtUJukyRclGIx6lt\n2/wx3x9/DH/6k+95UpFnXCaCwwWH+WzrZ6RuSmX2ptks3r6Ybs26+VBvfwl92/alVtVaQZcpUioK\n8Thz8CD87W/wj3/AqFHw299CHQ0my+Rg/kHmb53P7I2zmb1pNmk70jir6Vn0bdOXvm370rdNXx35\nIjFPIR4nnIM334QHHoBzz/U9v0vbaVBOLjc/l0XbFjE3Yy5zM+YyP2M+DWs29KEeCvYzm5ypk44k\npijE48CCBfCb38D+/X4E3r9/0BUlhyJXxNqda4+E+twtc9l9cDfntz6fc1ueS59WfejTso9G6xIo\nhXgMS0+Hhx/2LWJHj/YLNlSuHHRVyS1zfybzt85n0bZFLNq+iMXbF1OnWp0jgd6nZR96t+yto2Ak\nahTiMWjbNn+Y4Ntvw/33wy9/Gd9nWyYy5xwbsjYcCfVF2xex7KtltKrXit4te3NO83Po3qw73Zt3\np2ntpkGXKwlIIR5D9u71zanGjvWtYh96KDYbVcnJFRQVsHbnWhZvX8zyzOWk7UhjeeZyalapSY/m\nPejerLuIcrh5AAAKN0lEQVS/bt6dTo06UbmS/rySslOIx4CcHPjnP/1891VX+amTNm2CrkoiyTnH\nluwtRwI9bUcaaTvS+PrA13Rt2pWzmpxF16ZdObPJmXRt0pXW9VprbVIJi0I8QMeG9+DB/ozLM84I\nuiqJpuxD2azIXMHqnatZs3PNkesDeQfo0qQLXZscDfYzm5xJ2/ptFe7yLQrxAOTk+MZU//gHXH65\nD+9Ea1Il5bPn4B7W7Fxz5PJNuOcczqFTo050atyJTo06cXrj04/cb1yzsQI+CSnEo2jnTr8c2rPP\nKrylbPYe2sv63etZv2c9n+/+nPV71rN+t79dySodCfcjAd+4Ex0adKBRzUYK+ASlEI+CL7+EJ5+E\niRPhuuvgvvs0bSKR5ZxjV+6uo+F+TNBv3LsR5xwdGnagfYP2dGhw9Pqbx+pVT9KOaQlAIV6Bli71\nZ1bOmgV33OEPFWzRIuiqJBllHcxi095NbNy70V9nbTx6e+9GalSpcSTY29VvR+t6rWlTvw2t67Wm\ndb3WtKjTQkfRxCiFeITl5/uFiMeMgS++gHvv9QFet27QlYkU75tR/Ma9G9mYtZEt2VvYmrOVrfu2\nsjVnKxnZGezK3UWzOs18uNc7Gu7H3m9WpxnVKlcL+uskHYV4hGRmwvPPQ0oKnHoq/OIX8KMfQVUt\nMiMJIL8wn6/2f0VGdoYP+NAlI+fo/a8PfE296vVoXqc5zes0p0XdFjSv3fzI/SOP1WlOwxoNNUcf\nIQrxcigqgjlz/IIM77/v57vvugu6dw+6MpHoK3JF7M7dzY79O9ixfwdf7f/qyO3j7+fm59KsdjNa\n1G1B09pNaVKrCafUOuXode1v369XvZ5C/wQU4mWwZQuMHw/jxvk2sLfcAjfdpLMrRcJ1qOCQD/Z9\nX7Ezdyc7D+xkZ+5OduXuOnp94Oj9wwWHiw33JrWa0LhWYxrWaEijmo1oWLMhDWs0PHKdDOutKsTD\ntG8fTJ0KL78MS5bA0KEwYgT06pX461iKBO1QwaHvBPs393cf3E3WoSz2HNxD1sEssg5lkXUwi72H\n9lKzas0jod6oZiN/+5iQPzb469eoT/3q9alXvR71qtejdrXacdF2WCF+EgcOwHvv+T7es2ZB375+\nBZ2rr4YaNYKuTkROxjnHvrx93wn3b673HNzjb4fuZx/OJudwzpFLbn4udarVoV71et8K9+MvJ3qu\nTrU61K5WmzrV6lC9cvUKmw5SiB8nKwumT4f//MdfX3CBXzV+yBBo1Cjo6kQkWgqLCtmXt+9bwZ5z\nOIfsQ9nfeSzncM63fglkH87mQN4BDuQfYH/efgqKCnyoV639rXD/1v2qRx8/0ba1q9WmdtXa1Kpa\n68ilWpVqyR3izsGaNX7H5HvvQVoaXHyxD+1rroHGjaNekogkmPzCfA7kH+BAng/1b8K9xPuh62Of\nO1hwkNz8XHLzczmQd4DCRwvLF+JmNhj4B1AJeNE599ditnkauAI4ANzsnEsrZpuohLhzsGkTpKb6\ny+zZfk77Bz+A738fLrlEvbtFJD4456hUqdIJQ7zEGX0zqwSMAQYBXYFhZtb5uG2uADo65zoBI4GU\ncldeCocO+eXNxozx89nt28OFF8KMGf565kwf6s88A1deGbkAT01NjcwbxQh9n9im7xPbKur7lDTP\nHs5u2XOB9c65zc65fOB1YMhx2wwBJgA45xYA9c0s4osSOgdbt/pQfvppGDkSevb089ijRsGKFXDR\nRT68t2+H117z25xxRsUcXaJ/hLFN3ye26ftERpUwtmkFZBxzfys+2E+2zbbQY5nhFuKcP1IkK8uf\nGZmR4S9btvjrjRv9epS1akGXLv7SrZtfk7JHD02PiEhyCifEI+qyyyAvz/cfycvzx2ZnZflly6pW\n9SfTNG3qV7355tKzp58i6dxZR4+IiByrxB2bZnY+MNo5Nzh0/yHAHbtz08xSgNnOuTdC99OB/s65\nzOPeK/jjC0VE4tCJdmyGMxJfBJxmZu2Ar4ChwLDjtpkC3AW8EQr9vccH+MmKEBGRsikxxJ1zhWb2\nC2AGRw8xXGtmI/3TbqxzbpqZXWlmX+APMRxRsWWLiAhE+WQfERGJrKh1fjGzwWaWbmafm9mD0frc\nimBmL5pZppmtCLqWSDCz1mb2sZmtNrOVZnZ30DWVlZlVN7MFZrYs9F0eDbqmSDCzSma21MymBF1L\neZnZJjNbHvpvtDDoesrLzOqb2SQzWxv6f+i8qH5+NEbioROGPgcGAtvx8+xDnXPpFf7hFcDM+gH7\ngQnOubODrqe8zKw50Nw5l2ZmdYAlwJA4/u9TyzmXa2aVgbnA3c65uA4LM7sX6AXUc85dFXQ95WFm\nXwK9nHNZQdcSCWb2MvBf59w4M6sC1HLO5UTr86M1Eg/nhKG44Zz7FEiIf4AAzrkd37RJcM7tB9bi\nj/OPS8653NDN6vj9PnE9Z2hmrYErgReCriVCjCjOAlQkM6sHXOScGwfgnCuIZoBD9H6QxZ0wFLch\nkcjMrD3QA1gQbCVlF5p6WAbsAGY65xYFXVM5/R24nzj/ZXQMB8w0s0VmdnvQxZRTB2CXmY0LTXeN\nNbOonnqYEL8NJTJCUylvAb8KjcjjknOuyDl3DtAaOM/Mzgy6prIys+8DmaG/lCx0iXd9nXM98X9d\n3BWanoxXVYCewDOh75QLPBTNAqIV4tuAtsfcbx16TGJEaC7vLeDfzrl3g64nEkJ/1s4GBgddSzn0\nBa4KzSO/BlxiZhMCrqlcnHNfha53Am/z3TYe8WQrkOGcWxy6/xY+1KMmWiF+5IQhM6uGP2Eo3vey\nJ8qo6BsvAWucc08FXUh5mNkpZlY/dLsmcBkQlztoAZxzv3POtXXOnYr//+Zj59yNQddVVmZWK/QX\nH2ZWG7gcWBVsVWUXOqkxw8xODz00EFgTzRqi0jvlRCcMReOzK4KZvQoMABqb2Rbg0W92bMQjM+sL\n3ACsDM0lO+B3zrkPg62sTFoA40NHRFUC3nDOTQu4JjmqGfB2qAVHFWCic25GwDWV193ARDOrCnxJ\nlE921Mk+IiJxTDs2RUTimEJcRCSOKcRFROKYQlxEJI4pxEVE4phCXEQkjinERUTimEJcRCSOKcRF\nQsxMywpK3NEZmyKAmQ0DauF7kH/mnFsacEkiYYlK7xSRWGZm/wRWOufGBl2LSGlpJC5Jzcy64pej\nq+2cKwy6HpHS0py4JLtOQI4CXOKVRuKS1MysKb6f9X+AFcAO59x/gq1KJHwaiUuyqwxMA7Lx+4jy\ngi1HpHQ0EpekZWbVgYnALdFeoVwkUjQSl2Q2CnhbAS7xTIcYSjI7DT+dIhK3NJ0iScvMzgE+AOYA\n/3LOfRRwSSKlpukUSVrOuWVAV2At8JKZfWFmA4KtSqR0NBIXCTGz24HBzrkfB12LSLg0Jy5Jz8y6\nA/2AtsBTAZcjUioaiYuIxDHNiYuIxDGFuIhIHFOIi4jEMYW4iEgcU4iLiMQxhbiISBxTiIuIxDGF\nuIhIHPv/bXGLA8Ug8v0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3b668c80d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rho = np.exp(-psi[:,0]) \n",
    "plt.plot(xi,psi[:,0],label ='$\\psi$')\n",
    "plt.plot(xi,rho,label =r'$\\rho/\\rho_c$') \n",
    "plt.xlabel(r\"$\\xi$\",fontsize=15)\n",
    "plt.xlim(0,xi_max)\n",
    "from matplotlib.legend_handler import HandlerLine2D\n",
    "plt.legend(loc='upper left',prop={'size':13},numpoints=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that since ``rmax`` is the only parameter that needs to be adjusted for different models of the Lane-Emden solution. To make sure that the boundary of the sphere is properly filled with the Bonnert-Ebert values, we extend the density profile in ``density.txt`` to $\\xi$= 20. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "xi = np.arange(1e-6, 20, 0.01) \n",
    "psi = integrate.odeint(solvr, IC, xi)\n",
    "rho = np.exp(-psi[:,0]) \n",
    "np.savetxt(\"FLASH_4.3/density.txt\",rho)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inside ``Simulation_initBlock.F90``, we can read in the values of the array \n",
    "~~~fortran\n",
    "  open(12,file=\"../density.txt\")\n",
    "  read(12,*,IOSTAT=io) dens_arr\n",
    "  close(12)\n",
    "~~~"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have a 1D array of non-dimensional density values from numerical integration, we need to assign each cell in the simulation box with the density values that corresponds to the right $\\rho(r)$. Since we are working with a radial distribution on a cartesian grid, we don't have the exact radius value at each cell, so need to do linear interpolation to map the cells to somewhere between their closest density values. As explained by [Wikipedia](https://en.wikipedia.org/wiki/Linear_interpolation), linear interpolation basically weighs what the y values are  based on how close the x values lie to the previous and subseqeunt datapoints. In our case, \n",
    "\n",
    "$$ \\rho(r) = \\rho_0 + (\\rho_1-\\rho_0)\\frac{r-r_0}{r_1-r_0}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Translating this into our Simulation_init.F90 and picking dr ($\\Delta \\xi$) as 0.01 to correspond with np.arange used in the numerical integration: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    if (rc <= rmax) then\n",
    "       rc0 = int(rc/dr)+1 !to prevent from hitting index 0 which is yields zero density\n",
    "       rho0 =  rho_c*dens_arr(rc0,1)\n",
    "       rho1 = rho_c*dens_arr(rc0+1,1)\n",
    "       rc0 = rc0*dr\n",
    "       rc1 = rc0+dr\n",
    "       rhoZone = fattening_factor*(rho0+(rho1-rho0)*(rc-rc0)/(rc1-rc0) )!linear interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulation_init.F90 reads in values up to ``rmax`` $\\xi$= 16.90 for initialize the sphere, note that this has a much higher density contrast than the Bonner-Ebert sphere. This ensures that the sphere in our toy problem will actually collapse in a short time. We also multiply the sphere's density by a \"fattening_factor\" of 100 for the same reason."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "At xi =  16.900001\n",
      "Density is:  7.2918475616e-20 g/cm^3\n",
      "density contrast is  150.853400418\n"
     ]
    }
   ],
   "source": [
    "idx = np.where(np.isclose(xi,16.9))[0][0]\n",
    "print \"At xi = \", xi[idx]\n",
    "print \"Density is: \", rho[idx]*rho_c, \"g/cm^3\"\n",
    "print \"density contrast is \", 1./rho[idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now run our simulation to check that the sphere is initialized with the right density distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/global/u2/d/dorislee/dlee/FLASH4.3/object\n"
     ]
    }
   ],
   "source": [
    "cd ~/dlee/FLASH4.3/object/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import yt\n",
    "from matplotlib.legend_handler import HandlerLine2D\n",
    "from mpl_toolkits.axes_grid1 import AxesGrid\n",
    "from yt.mods import *\n",
    "yt.mylog.setLevel(50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def compareplot1Dprofile(timestep):\n",
    "    fattening_factor=100\n",
    "    rho_c = 1.1e-19*fattening_factor\n",
    "    pf= yt.load(\"sphere_hdf5_chk_{}\".format(str(timestep).zfill(4)))\n",
    "    sp = pf.sphere(pf.domain_center, (0.32,\"pc\"))\n",
    "    rp = yt.create_profile(sp,'radius','density')\n",
    "    plt.plot(rp.x.value*1.05e-17,rp[\"density\"].in_units(\"g/cm**3\").value/rho_c,label=\"Simulation\")\n",
    "    xi_max = 6.451\n",
    "    def solvr(Y, t):\n",
    "        return [Y[1], exp(-Y[0])-2/t*Y[1]]\n",
    "    xi = np.arange(1e-6, xi_max, 0.01) #start at small poisitive number to avoid div-by-0\n",
    "    psi = integrate.odeint(solvr, [0, 0], xi)\n",
    "    rho = np.exp(-psi[:,0])  \n",
    "    plt.plot(xi,rho,\"--\",label ='Integration',color=\"red\",linewidth=5)\n",
    "    plt.xlabel(r\"$\\xi$\",fontsize=15)\n",
    "    plt.ylabel(r\"$\\rho/\\rho_c$\",fontsize=15)\n",
    "    plt.legend(loc='upper right',prop={'size':12},numpoints=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEXCAYAAACDChKsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucVfP++PHXey6lyXRXVJoumiR0XEPRuByVI7l8UR2l\ncpJD5HBSQsqJ5Lim4xKKkCgcjh+pOIMOXShK6UK60/0y3Zvm8/tjzZ5Zs+977evs9X4+HvuxZ9Ze\ns9Znd9nv+Xzen8/7I8YYlFJKqUhlJLsBSimlKicNIEoppRzRAKKUUsoRDSBKKaUc0QCilFLKEQ0g\nSimlHEmJACIir4jIJhFZFOScsSKyUkS+F5E/JLJ9SimlfKVEAAEmAp0CvSgiXYAWxpiWwADghUQ1\nTCmllH8pEUCMMbOBHUFO6QZMKj13LlBTRBokom1KKaX8S4kAEoZGwDrb9xtKjymllEqSyhJAlFJK\npZisZDcgTBuA423fNy495kNEtLiXUko5YIyRSM5PpR6IlD78+RDoDSAi5wA7jTGbAl3IeD8uuABT\nUoIxptI/HnzwwaS3Qd+fvj+3vTc3vD8nUqIHIiKTgQKgroisBR4EqgDGGDPeGPOxiFwmIj8De4G+\nQS/41VfwzTcwZQosWACjR4NEFFiVUkqFkBIBxBjTM4xzBoZ9wQ4drMfgwVYAOf30Ci/v2gWTJsHt\nnVfCX/8Kr78Oxx0XecOVUsrFUmkIKz68ggfAzp0wdfTPcOGF8Nln1vO2bUloXOQKCgqS3YS40vdX\neaXze4P0f39OiNOxr1QlIibUe9rw1Srkwo40PLK+/GCHDjBzJhx1VJxbqJRSqUdEMBEm0d0XQEpK\nOHTK6VRZ+oPva927w+TJmi9RKgxNmzZlzZo1yW6GilBeXh6rV6/2Oa4BhPB6IL9/9C3VrriEmmaX\n74svvww33RSn1imVPko/cJLdDBWhQH9vTgJI+udA/Djc9kz6130PsrMrvnDDDVYvRCmlVEiuDCAi\n8L8qF1m9DYCqVeHFF62pWdWrJ7dxSilVSaTENN5EK0tx9O4NGzfCpZf6na2llFIqMFf2QADKhgCH\nDtXgoZTLTJ48mc6dO8fl2n379mX48OGOfz43N9dvkjsVuTKAiNgCSDAlJVBUFPf2KKXiY/bs2bRv\n355atWpRr149zj//fL777jt69uzJ9OnTk908LrzwQiZMmFDhWFFREU2bNk1OgyLk2gAS0qpV1gLD\nfv3i3h6lVOwVFRXRtWtXBg0axI4dO9iwYQMPPvggVatWTXbT0oYrAwgE6YEYAy+8AKeeCl9+CdOm\nQQr8pqKUisyKFSsQEa677jpEhKpVq3LJJZdw8skn89prr3H++eeXnZuRkcHzzz9Pfn4+NWvWZPjw\n4axataqs99K9e3eKi4sBfH7W8/OrVq3yacPOnTvp2rUr9evXp27dunTt2pWNGzcCcP/99/PVV18x\ncOBAatSowR133OFzrd27d9O7d2/q169Ps2bNePjhh8uu7WnH4MGDqVOnDi1atEh4r8qVASToEFav\nXlZ9rL17y4/ddhvs35+QtimlYiM/P5/MzEz69OnD9OnT2blzZ4XXxWsoYsaMGSxcuJA5c+bw2GOP\nMWDAACZPnsy6detYvHgxb731VsCf9f7eo6SkhH79+rFu3TrWrl1LTk4Ot912GwCjRo3i/PPPZ9y4\ncezevZuxY8f6XGvgwIEUFRWxevVqCgsLmTRpEhMnTix7fd68ebRu3Zpt27YxePBgbkrwGjbXBpCA\n/CXWVq2yKvoqpYISif8jXLm5ucyePZuMjAxuvvlmjjnmGK688ko2b97s9/whQ4ZQvXp1Wrduzckn\nn8yll15KXl4eubm5dOnShYULFwa8V6AFlXXq1OGqq66iatWqVK9enXvvvZcvv/wyaLs91yopKeHt\nt9/m0UcfJScnh7y8PO6++25ef/31snPz8vLo168fIsKNN97I77//HvD9xYMrAwgE6YH8+c/gr2ja\nmDHwyy/xbJJSlZ4x8X9EolWrVkyYMIG1a9eyZMkSNmzYwJ133un33Pr165d9Xa1aNRo0aFDh+z17\n9kT857F//34GDBhA06ZNqVWrFh07dmTnzp1hreDfunUrxcXFNGnSpOxYXl4eGzaU76V37LHHVmij\nMcZRO51yZQAJ+luMCDz3HMa+Sj0jA/r0gdzceDdNKRUn+fn59OnThyVLlkR1nerVq7Nv376y73//\n/feA5z7++OOsXLmS+fPns3PnzrLehyeABBr6AqhXrx7Z2dkV6o2tWbOGRo0aRdX+WHJlAIEQv8m0\nbs22Pn+3vu7cGX74wVqpbvsNRSmV2pYvX86TTz5Z9hv7unXreOuttzjnnHOium7btm1ZsmQJixYt\n4uDBg4wcOTJgINizZw/VqlWjRo0abN++nREjRlR4vUGDBn6T72Al06+77jruu+8+9uzZw5o1a3jq\nqafo1atXVO2PJVcGkHDWgWzpP4w/HzMDPvkETj45MQ1TSsVMbm4uc+fOpV27duTm5nLeeedx6qmn\n8sQTT/icG25SHKBly5YMHz6ciy++mPz8fJ8ZWXZ33nkn+/bto169epx33nlcdtllFV4fNGgQU6dO\npW7dumVDa/Z7jx07lpycHJo3b84FF1zADTfcQN++gTdkDdbueHBlNd6tW+HEE63nQH78Ebp0gXXr\nYtxApdKEVuOtnLQabwyE+ndvjLUQXSmllH+uDCDhDGGVlAQ5Z+lSmDEj5u1SSqnKxLUBJBS/UwZ/\n/tlaaHjyydC3ry4uVEq5misDCEQ4hHXkCNx8s5U4eeMN68WNG62SJ0op5VKuDCDhDGFV6IFkZsL2\n7VYgsRs9GhK4aEcppVKJawNIKD5DWCNH+v7gli0wblxM26aUUpWFKwMIOEiit2kDPXv6nvjkk3Do\nUEzbppRSlYErA0i4Q1g+03gffJCSjEzr6ypV4JZbYO5c62ullHIZ1waQUPzOwmrZkqUFt/Kf/Lvg\n11/h+eehWbO4tFEp5V6XXXZZhaq7qcqVAQQiTKLbfH39WF475Qlo2DA+DVNKxUSzZs34/PPPQ57n\nb1vZRBo5ciS9e/eucOzjjz9OqZpXgbgygITbA/G3Er2kRFeoK6XCc8R75maacWUAAecr0YOuUFfK\n7UaM8L8LlFcVWsfnO+Bv69dPP/0UCLyt7LJly7j00kupW7curVu3ZurUqWXX2759O127dqVmzZq0\na9eOBx54wGd73Oeee478/Hzy8/MBq6hikyZNqFmzJmeddRazZ88G4NNPP+WRRx7h7bffJjc3l9NO\nOw2o2CsyxjBq1CiaNm3KscceS58+fdi9ezdglXfPyMhg0qRJ5OXlUb9+fR555JGY/dmF4soAEvE6\nEK/jfnsgCxZYW9+W7puslEodc+fOrbD1a79+/QD/28ru27ePSy+9lBtuuIGtW7cyZcoUbr31VpYt\nWwbArbfeSm5uLps3b+bVV1/ltdde86mC+8EHHzB//nyWLl0KwNlnn82iRYvYsWMHPXv25Nprr+XQ\noUN06tSJYcOGcf3111NUVOR318OJEycyadIkvvjiC1atWkVRUREDBw6scM7//vc/Vq5cyaxZs3jo\noYdYvnx5PP4Yfbg2gIQSKIB4D2GZr2az5awucMYZ8NxzYPtNRSmVGpo2bVph69fffvst4NavH330\nEc2aNaN3796ICG3btuWaa65h6tSplJSU8N577/HQQw9RtWpVWrduzY033uhzjWHDhlGzZk2qVq0K\nQM+ePalVqxYZGRn87W9/4+DBg2F/yE+ePJm77rqLvLw8cnJyGD16NFOmTKGk9INIRBgxYgRVqlTh\n1FNPpW3btvzwww8O/6Qi48oAAs6r8VYIIEOGIBeczzHfTi8/YcwYHeNSKsV4b/0KBNz6dc2aNcyZ\nM4c6depQp04dateuzeTJk9m0aRNbtmyhuLiYxo0bl51//PHH+1zD/jpYOxOedNJJ1K5dm9q1a7N7\n9262BttPwmbjxo3k5eWVfZ+Xl0dxcTGbNm0qO2bffjcnJydh29q6MoBEU423Qs/Ea3MYwNq98LPP\nom6jUioxvIefjj/+eAoKCti+fTvbt29nx44d7N69m3HjxnHMMceQnZ3N+vXry85f52fTIPs1Z8+e\nzT//+U+mTZvGjh072LFjBzVq1AhrW1uAhg0b+mxrm52dXSFoJItrA0goYQ1hXXABB8841/ekxx+P\nqn1KVVojRpT/57E/giXRIzk/Dry3lb388stZsWIFb7zxBsXFxRw+fJhvv/2W5cuXk5GRwdVXX82I\nESPYv38/y5YtY9KkSUGvX1RURHZ2NnXr1uXQoUM89NBDFBUVVbj/6tWrA27O1aNHD5566ilWr17N\nnj17uO++++jevTsZGdbHdzI39XJlAIEYDWGJsKv/331PKiyE336LtolKqSgE+83e/pr3trJHH300\nM2bMYMqUKTRs2JCGDRsydOhQDh48CMCzzz7Lzp07Oe6447jxxhvp2bNnWa7D3307depEp06dyM/P\np1mzZuTk5FQY9rr22msxxlC3bl3OPPNMn2v069ePXr16ccEFF9CiRQtycnIYO3ZswPslclvblNnS\nVkQ6A09jBbVXjDFjvF6vAbwBNAEygSeMMa/6uU7ILW0PHICaNaH034Nfn30GnTr5Tqr65z/h009h\n1izr+19/PoLJb0Vz8wvUqQO33mrNxrKNuSqVjnRLW8vQoUPZtGkTEydOTHZTwhLLLW2zYtaqKIhI\nBjAOuBjYCMwXkQ+MMctsp90GLDHGXCEi9YDlIvKGMSbiebOxnIVVIpn8o8ooXn1iG/TpA9WrR9oc\npVQlsnz5cg4dOsQpp5zCvHnzeOWVV5K6kj2ZUiKAAGcDK40xawBEZArQDbAHEAPkln6dC2xzEjzK\nLhajhYQlJfC2dOfV25y2RClVmRQVFdGjRw9+++03GjRowODBg+natWuym5UUqRJAGgH2qQzrsYKK\n3TjgQxHZCBwNXO/0ZrFcSKilTZRylzPPPJOVK1cmuxkpIVUCSDg6AQuNMReJSAtgpoicaozxmfA8\nwjaDo6CggIKCggqvhzuE5Xm2n+8zhOUngBgD334LZ50V+j5KKZUMhYWFFBYWRnWNVAkgG7CS4x6N\nS4/Z9QVGAxhjfhGRX4ETgW+9LzYijCmA4fRAPM+RBpDNm+FPlxk2T/kcXnsNXnkFsrNDtkkppRLF\n+5frkSNHRnyNVJnGOx84QUTyRKQK0B340OucNcAlACLSAMgHVuFAJD0Qf70L7xxIhXMOH6bKtMnM\n3HEGXHIJvP46vP++k2YqpVRKS4keiDHmiIgMBGZQPo33JxEZYL1sxgOjgFdFZFHpj91jjNnu/J7B\nX/cEBe/z/PVAPOeJALfcQu0JE6ht/6Gnn4brrnPaVKVSUl5eXkLXHKjYsJdFiVZKBBAAY8x0oJXX\nsRdtX/+GlQeJWrhJdPuzR6AAcuQIZGUBvXqB95S+b76B+fM1KaLSyurVq5PdBJVkqTKElVCRJtG9\nj/sLIGXHOnbkYOu2vhd85pmI26mUUqnMlQEEwu+BeOdA/K0DqXCeCDt7D/K94PTpsG+fo7YqpVQq\ncmUACacHEmkOxH5s12U92Mwx1jctWlg5kFWrICfHeaOVUirFpEwOJJHiOoQFHMk+isE8zmsf1II/\n/QkyM6NrsFJKpSBX9kA8gg1jRZpErxBAjsAkemO6XgGZmcybBz16xKbNSimVKjSAhHgt3BzIkSO+\nxzzPmzfDv/9tVQFWSql04doAEmoYK9whLE/g8HfM/nzgAHz9tfP2KqVUqnFtAIHgPZBokujeQcXz\n/eczimHqVHjgAeeNVkqpFOHKJDqEXkwYbAgrVADxDhwZu7YzsvrL3PzUODi0zrp5377QvHl0b0Ip\npZLItT0Qp0NYIdeB4DWEZQwX3t+e4XuHcOyhdeUXffZZx21XSqlU4NoAAs5mYYUzjbfCMRFWnd/H\n9wYTJkBRUaRNVkqplOHaABJNDyTYcJX9a8/z8gv6czCzWsUL7d4Nr74aUZuVUiqVuDaAQHhJ9Fjk\nQPZXq8PXLXr73uTddyNrsFJKpRDXBpBwk+j+hrDCzYHYX/uszR3lJ7RrB5Mnw8yZzhqvlFIpwNWz\nsIKJZiW6dw/kyBHYXO8kNtz2CMNmXshrc85x3nCllEoRrg0gEL9pvP4WEmZmwtb+97Lwy+jarJRS\nqUKHsAIItJAw4llYlAeQzMyKyXallKrMXB1Agol0HUiwWVhHjkBGhrVjYXGx8zYrpVQq0SGsEK/F\nKgcSsAeyZw+89hpUqwb9+kX8HpRSKllcG0CcljIJaz8QP7WwMjO9eiC//grjxrH32VeofngXNGpk\n7aeene34PSmlVCK5OoAEE9ceyNat0KoVHD5Mdc8PbdgA770H11/v4N0opVTiuTYHAvGvxusdQLKy\nSo/Vqwddu/re9OmnI2q/Ukolk2sDSDRDWKEWEgabhVU2hDVokO9N58yBuXPDfg9KKZVMrg4gwcSy\nFpbfJPr558Npp/neWMubKKUqCdcGEIjNLKxwdyT0SaKLlPVCjAj/L+sKdr3/OYwZ4+zNKKVUgrk2\ngMRqS9uoFhJ2784LtYayZuZKJl39AW/9fmHohimlVIpwbQAB59V4I9pQigA9EICqVRlTazSmeQt6\n97a2CFFKqcrCtQHEaTXeWC8kLC62AkvnzvD777BgQeTvRSmlksHVASSYaIawgi0kDBRAMjPh5pvh\nxRcjex9KKZUsrg0gENtqvPbAEKgHkpHhG4A8AQTgppvgnXegaMVvMHw4fPVV5G9KKaUSxNUr0Z0u\nJPQcF4ksBwLlvZCMjPLXPAHkuI3f8e8aT5Nz0ttw5DB8/7013VcppVKQqwNIMMGGsDzPgQJIoFlY\nUL6Y0FPyqri49LXPPoNLLqGj/WYffQS//AItWkTwzpRSKjF0CCvEa4F6IIGeIXgPxDuRXjaE1bEj\nNGni24hnnw3rvSilVKK5NoA4LWUSTgDxzoGUlFQcwrJP5S0LIFlZcPvtvg2ZMAF27w75fpRSKtFc\nHUCCCWcIC4KXMgk0hOW3BwJWFj0np+INi4rg00+DN1YppZIgZQKIiHQWkWUiskJEhgQ4p0BEForI\njyLy32jvGU0SPZIeiL8kuv08T0Kd2rWhTx8ADmcdxZet+sPixXDttWG/J6WUSpSUSKKLSAYwDrgY\n2AjMF5EPjDHLbOfUBP4FXGqM2SAi9aK7Z/yGsELlQDxDWPbjZQYNgsaN2X11f648tx5L68GxYb8r\npZRKnFTpgZwNrDTGrDHGHAamAN28zukJvGuM2QBgjNkazQ2jHcIKpwfibwjL3gOpMHzlkZ8P995L\n3Vb16NkTxo4N6+0opVTCpUoAaQSss32/vvSYXT5QR0T+KyLzRaRXtDeNZhZWoEAC/nsgnmEqew/E\nbwCxuesuGD8eduwI/j6UUioZUiWAhCMLOB3oAnQGHhCRE5xeLJr9QII9278ONY03VABp3hyuvBKe\neCJ4W5VSKhlSIgcCbADsiyAalx6zWw9sNcYcAA6IyJdAW+Bn74uNGDGi7OuCggIKCgr83tRJNd5A\nPY9wZ2HZp/HaV6EHMny4te/U3y5eRN2Jj8MDD0DLlsF/SCmlQigsLKSwsDCqa6RKAJkPnCAiecBv\nQHegh9c5HwDPikgmUBVoBzzp72L2ABJINNV47c+RbCgFvj0QnyS6VyOarPiM2dX/Sd2LZljHqleH\n558P8kNKKRWa9y/XI0eOjPgaKTGEZYw5AgwEZgBLgCnGmJ9EZICI3Fx6zjLgU2ARMAcYb4xZ6vSe\n0Q5hBcuBhDuNN9QQFhMnwh//SOt1M8qPvfoqbNkSvPFKKZUAKRFAAIwx040xrYwxLY0xj5Yee9EY\nM952zuPGmDbGmFONMVHX+HDSAwlnFlaohYThJtG55hrIza147MAB+Ne/gvyQUkolRsoEkEQLtxpv\nLEqZOO6B1KwJAwb4Hh83DvbtC/KDSikVf64OIMHEs5hi0IWE3gYN8o0ye/bAvHkhflAppeLLtQEE\nYjuE5b2hVHa2w4WE3ho3hp49ATiUW4eXGz4Aa9dCgJllSimVKKkyCyvhoi1lEmohYXZ26B5IWAEE\nYMgQOOssMnv35R+nVOeMDXBa/TB+Timl4si1PZBYDmFlZPjmQIIFkIh6IAAnnQQDB5JZozp/+Qu8\n9FIYP6OUUnHm2gACzqrx+ut5ZGX574GEWkgYdgCx6dcPpkyBvXsj+zmllIq1mAQQEbnY9nW+iGTH\n4rrxFKuFhN75Ds+xKlVC90DCSqJ7adTI2ib97be9GqqUUgkWqx7I1bagsRa4LkbXjZtwh7DCyYEE\n6oFEPY03gFtugWfHGsxH/w/at4fZsyO/iFJKRSlWAeTH0jLslNaqqhT1Y53OwhLxHcIKdxaWoyS6\nXUkJnYqm8eayM5Cul8M338DDD0d4EaWUil6sAshPIvKJiNwgIu2Bs2J03biJJomemRleDiTQnuhR\n9UDef5+M66/lpIMLy49Nnw7ffRfhhZRSKjoxCSDGmELgNuBE4ApgUiyuG29OV6LbA0agHEhMp/Ha\ndetm1Xn39sgjEV5IKaWi4yiAiOX/RORBEekjIjWMMauMMfcbY4YYY36NdUNjzWkS3RgrCISTAwk0\nhBVNEp2sLBg61Pf4e+/BUse1JZVSKmJOeyDjsMqt/wF4EPhZRHrHrFUJEMshrEh6INFO4wWgd29r\nOpbdccdZK9SVUipBnAaQlcaYa4wxVxljmmHtEnidiPSMYdvizuk0Xu8hLH89kHCm8ToOIFWrwuDB\nAGzLzePDLs/DqlXQubODiymllDNOA0hdkfLf4Y0x3wFdsXYIrBScVuP1DGE5nYUVi2m8APTvD5Mm\nsX3OSm6afwt7io9yeCGllHLGaQD5BCgUkYs8gcQYY7C2na0UohnCysoqP+7ZljbUOpCM0j/pqJPo\nHjk50KsXLU/K5sILtbyJUirxHAUQY8zXwN3AaOB3EZkmIhOAlF+Bbud0CCvSabyBeiCOkuh+DBkC\nTz4Jhw5Ffy2llAqX42m8xphvjTHtgI7AZ0A14G4R2VIaUHrFqpHx4LQab6AhrEhKmcSkB2JzxhnQ\nujW8+WbpgYMHYevW6C+slFJBRL0OxBizzBjzvDGmhzGmEXAu1hDXKVG3Lo6iHcIKNgsrVDHFmORA\nvAwdCk88epiSF1+C/Hy4447YXFgppQKI+X4gxpifgZ9jfd14cFKN1zsHEqgHEreFhP4UF3Phmjf4\ndM1DZNxSugRn3Tq47z5o0yYGN1BKKV8heyAi8oCIdEhEYxIp2oWEwWZhhdpQKuY9kEOHkHuH0uig\nbf2mMTByZAwurpRS/oUzhDUWOEFEXhKRMSJyTrwblQjRVOMNJwcSzhBWrJLo5OT4X50+dSr88EMM\nbqCUUr5CBhBjzC5jzKvGmP7AGOBkEXlZREaLyJnxb2L8xHMhoacHYkzFYopxGcICGDAAGjb0Pf7g\ngzG6gVJKVRRREt0Ys90Y87Ix5i/Ak8CZIvKKiIwSkT/Ep4nx4TSJHk4tLPssrJIS616e+8UriU61\najBsWMVjZ50FgwbF6AZKKVVRNNN4txhjXjDG3IRVG6uDiEwQkZEicnLsmhg/TqvxhqqFZZ+F5T1M\nFbceCMBf/gKNG7O/0QncWu8dDnwxFy68MIY3UEqpcrEq5/67MWacMaYfMB64REQmioifgfnUEM2W\ntqGGsOyzsLwDSNx6IGDVyJoxg2q/LmX9udfy/AshullKKRUFRx9fpeVLrgHaAGuA94wxuwGMMRuA\np0vPqxujdsZcOENY9qEq7+OhZmF5hrCC9UBilkS3a90asLYHuegi6NcPataM8T2UUoo4l3M3xmyL\nom1xF6oHkpERXQ8k0BBW3HogNiefDJdfDqNHx+f6Sinl2nLu4VTjzcgInAOxJ9ED5UASPoTlZdQo\nePllWLYMq7TJAw/A/v3xu6FSylVcXc49mEA9EH9DWME2lLJP4YU4J9G9NGwII4fs46s/jca0aGFF\nlGefjd8NlVKu4tpy7hB6CMs7B2KM9Qh3HYi/IaxE9kCYNYtbn25J/1XDkN27rWOPPKKFFpVSMeHa\ncu7hzMKy9zQ8xyD8lehJSaLbNW6MbNpU8diuXbq4UCkVE64u5x6MvyEsY8oXBXovJAx3FlZCeyAn\nngg33eR7/IUXYPHiON5YKeUGri3nDqGT6N5DWJ7Euj25HmxHwoQvJPRn5EjIza14TARmz47zjZVS\n6S4mCwntjDE/G2NeMcbcE+trx1I4Q1jePRB/ASSlFhL6c+yxcP/9Zd9+W/sSzMLv4a9/jfONlVLp\nLuYBpLIIdyGhdw5ExDeARDKNN+E9ELDqYXXqxJH3PqB3gxn8vzWVotKMUirFpUwAEZHOIrJMRFaI\nyJAg550lIodF5Opo7xlNDyScDaU8Q1gZtj9l+0LCuCfRPapWhenTybzqCv75uDB4cHkQU0opp1Ii\ngIhIBtbq9k5Y5VF6iMiJAc57FPg0+nsGf93fNF5PABEJv5x70oewvFx2mbU+5KWXEntfpVT6SYkA\nApyNtbp9jTHmMDAF6ObnvNuBacDmWNw0nJXo/mZh+cuB2Gdh2cu5p8QQlo0IPP64lVv3LA1h3rzE\nNkIplRZSJYA0AtbZvl9feqyMiDQErjTGPA9EXWbWyTqQcJPoKbOQMIDTToNOneDFYWugWzdo1w4+\n/jjxDVFKVWpJ+Phy7GnAnhuJKohEM4QVbg4knB5IQnIg3g4f5plGT5E9eiSwzzo2cCD8+KO1Pa5S\nSoUhVQLIBqCJ7fvGpcfszgSmlJZOqQd0EZHDxpgPvS82YsSIsq8LCgooKCjwe9NIk+j2hYShZmGF\ns5DQs4Yk4f7zH2qN9pqn8OuvVpmTUaOS0CClVKIVFhZSWFgY1TVSJYDMB04QkTzgN6A7Vrn4MsaY\n5p6vRWQi8B9/wQMqBpBAnFTjjXQdSEosJPTnqqusnQr/+9+Kxx97DG64wVrBrpRKa96/XI8cOTLi\na6REDsQYcwQYCMwAlgBTjDE/icgAEbnZ349Ee89oh7Aq6ywswHrzzz1nNdLu8GF45pkkNEgpVRml\nSg8EY8xv114tAAAXmUlEQVR0oJXXsRcDnNsvNvcM/pq/aryBciCeoOCp2BssB5L0AAJWL2PIkLIh\nqyNZVcgc+SAMHpykBimlKpuU6IEkg9NSJv5yIN7l3UXKZ3ClZBLdY9gwaN6cotZn06neAvbdOcy3\nV6KUUgG4OoAE4wkgkeZAPEUYMzJSOInuUa0azJpF7uKvqXtBG554IoltUUpVOq4NIBB5NV77EFag\nWVie0iWeoaqUTKLbNWsGmZmMGWOlP9auTXJ7lFKVhmsDSDRDWN45EHsA8QSeQENYKZFE96NpUysl\ncv31cOhQ6cHt28F7QyqllCrl6gASTLCV6MFyIJ4eiGcIK5l7okfq7ruhfn245x7g/ffhpJOgb9/g\nkVYp5VquDSDgbCFhoById4XdQENY3j2QpCbRvWRkwGuPb+GSl6+Hq6+2eh+ffGJN+VVKKS+uDSBO\n14F4F1P03pHQ00sJNgsrZZLo3oyh1jUXc/nedyoe//vfYenS5LRJKZWyXBtAIPJqvOHUwvIEjECz\nsFJ5CAsRq5yJtwMH4M9/hoMHE98mpVTKcm0AcVKN1zOEFawWlj2JnrIr0YO5/HL/293++CP873+J\nb49SKmW5OoAEE2oWVqgkergLCVMugIC1YYitHlbRcfnwzTdw0UVJbJRSKtW4NoBA5KVMgi0k9PQq\nIllImGpJ9DI5OfDmm5CdzeYrb+bU4gWsP/bMZLdKKZViXBtAnFTj9Z6F5fl5+1CXfRqvMb5Bwt4D\nSbkkut3pp8OyZdR//0X+Mqg6N95Y8c9CKaVcHUCCCTYLyxN8AvVIMjPLh7oOH65kORC75lYF/aFD\nrcWFTz6Z5PYopVKKawMIRD+E5T1lF8p7IFAeQDJsf8qeHogxvsNbqSozE15/3douZOHC0oMTJsD0\n6Ultl1IquVL599+4clLKxHsIK1gPBKznQ4d8h7BKSsqHtkL1hFJF06bw9NPQr/s+5p19G9lvvAp1\n6lgRpUmTUD+ulEpDru2BOK3Ga5+F5S+A2Hsg/gKI5+e9j1cGPc9YzodbzrGCB1i1sq68EvbtS2q7\nlFLJ4doAApFX47XXwvLOgXjPwgL/ORCwvj94MMXzH/7cfz/H71hc8djChVovSymXcm0ACXchYaRD\nWKF6IGAFjkoZQJ5/Hho39j3+zjuaYVfKhVwdQIJxOoQVKgfiOX7gQCUMIPXqWVV6jzqq4vEzzrDq\nwCulXMW1AQSc7QfiZBZW2vRAAM4805qBVWpGne4cnPml/56JUiqtVcaPsJhwsg7EXgsr1DoQCN0D\nqWxJ9DI9esCSJZhqOYxfcC9T7xHGj688M8qUUrHh2gACzqrxhjMLK1QAycqyjlfKHojHqFEIMLEI\nOnaEESNg5MhkN0oplUiuHcJyUo3X3xCWve4VhDeEVWlnYfmRm2utJ5wyxVonAljR8dFHrW6WUipt\npcFHmDPhJtE9gcF+zBNA7HWvXJFED6B+fZg5E84/H+pVLeKG96+xDnz3nRVZKu1YnVIqGNf2QMBZ\nKZNQtbC8p/GmXRI9gCZN4LM3fqPtHR2t4AEwbRoMGKBVGJVKU2n0ERaZcKvxRjoLK9yFhJU6ie7P\noUOc8JcCKF5R8fgrr1jl4Z95RrPsSqUZ1/ZAwp2FFaycu+sWEgZTpQoMH+7/tWefhUmTEtsepVTc\nuTaAgPMhLO8A4glGnmGtcHIgaRdAwNo3/amnfI9fe6019VcplVZcG0CcVOMNVAsLKu5AmNYLCUO5\n8074xz/Kvn3nqF68e81kq4eilEor6fgRFhYnpUwCDWGB79RecGEPxOO++2DvXtixg/ybn6PLnzLY\nexB69052w5RSsZSuH2FhcVKN198QFvhO7YXgOZC0S6LbicAjjwDwBxE+/xz++EfYv9+alAXAzp1Q\no0bF3baUUpWKa//3OqnGG2gWFpQn3NO6nHskPPOdgdat4YsvrLWFY8aA2bIVOnSA/v3LN4hXSlU6\nrg4gwYTakdBfDiTcHogrAoiXFi1g9mz44NUdbGzzR1iyxCrKeNVV1nCXUqrScW0AgfCS6P7KuYtE\nlwNJ6yR6EI1qFPFl9S402vJ9+cGPPoKLLoLNm5PXMKWUI64NIE6q8QYbwrLPwvIewvIe5k/LhYTh\n+Pe/yfpuru/xefOgfXvtiShVybg2gEDkSfRwZ2G5ciFhOHr1grFj/Ufv3r2hevXEt0kp5VjKBBAR\n6Swiy0RkhYgM8fN6TxH5ofQxW0ROie5+zoewQs3Ccv003mBuvx2mToWqVcsOvVOzP7/0uD+JjVJK\nOZESAUREMoBxQCegDdBDRE70Om0VcIExpi0wCngpunsGfz3ShYT2WViuXkgYjmuugVmzoHZtuPxy\ntv3jOdp3EL74ItkNU0pFIiUCCHA2sNIYs8YYcxiYAnSzn2CMmWOM2VX67RygUbQ3jXQabzhDWNoD\nCVOHDjBnDkyZwl9vz+KNN+C666zai2W++AI2bkxaE5VSwaVKAGkErLN9v57gAeIvwCfR3DCcarxO\nFhJ690rsAcXDtUl0b/n5ZXmPSy6BL7+01oncdRccWbQELr8c/vAHmDEjyQ1VSvlT6X4HFpELgb5A\nh0DnjBgxouzrgoICCgoK/Fwn+H0C5UD87UgIgWdhgQ5hhatVK6tTctNV29n0Yjca7tsDe/ZA584w\nbJi1b67+oSkVE4WFhRQWFkZ1jVT537gBaGL7vnHpsQpE5FRgPNDZGLMj0MXsASSYSIspeo55ei/e\nhRP9LSS0P3voEFZgdWoU826V7mTs+6X8oDHw8MNWF+Wtt6BR1KOXSrme9y/XI0eOjPgaqTKENR84\nQUTyRKQK0B340H6CiDQB3gV6GWN+8XONiCSilEmgAKI9kCBWrSJj8SL/r61YUWH2llIquVIigBhj\njgADgRnAEmCKMeYnERkgIjeXnvYAUAd4TkQWisi8aO7pdAgrkmKKgYaw0m1P9JjKz4fvv4eLL/Z5\naf8z46FevSQ0SinlT0oEEABjzHRjTCtjTEtjzKOlx140xowv/bq/MaauMeZ0Y8xpxpizo79n4Nei\nWUgYbg/E9Un0QI49Fj79FB56qOwPeHbLPpz24BUsXpzktimlyqRMAEm0aIawAm0opTmQGMrMhAce\ngM8/h/bt6TD/ae67zyqbNX687e9l3z646SZYtSqpzVXKjVwdQIIJtJAw2BDWkSO+M7NAcyBR6dgR\nvvoKatakVy/ry5despaRLFiA1UuZMAHatLHqxR8+nOwWK+Uarg0gEHkpEycLCe3PHtoDiZAt2p94\nIsydC/36wd2XLubIY09YLxw4APfeC6edBjNnJqmhSrmLawNIrKvx+itlEiyAaBLduYwMuKlvCTOb\nDyDTeG1ItWQJdOoEK1cmp3FKuYhrAwg4X4keKgcSzhCWvxpZKgJr15L12zr/r914I7Rsmdj2KOVC\nrg0gTqrxOinnbn/28HyvPZAoNG0KS5fC3/5WYcOVvVKdHqse5vvvA/+oUio2XB1AgglWjTfcHEiw\nHoj9WTmUmwtPPgnffgvnnAPAUcOH0OG6hnTubG0/smaN7fwbbrAS7bpxlVIx4doAAs6m8YYzC0t7\nIAl22mnw9dfwzjtk3nM3t91mpUCaNYPTT4e//x12f/I/ePNNK9HevDk8/jgUFSW75UpVaq4NIE7W\ngXjXwnI6C0t7IHEgAtdeCzk5gNU5eegh+PFHK0780O2B8nM3b4bBg6FJE2uHRKWUI64OIMF4gkOo\naryhNpTyvGYXKLCo2DvuOHjx+s85//B/fV/cudN3w3qlVNhc/b8n0mq8Tra09bxmp0NYCfbdd34D\nRVHVeiw7r18SGqRUenBtAIlmCCvaWVg6hJVggwfDTz9Bz54Vup4L2t9Oxy45dOtmrXAv+7veu9fK\nq4wZA1u3JqfNSlUCrg4gwXgCiL8hrGjXgWgPJAny860k+k8/Qf/+UKcOHd+5jdWrrXWH/ftb1VCe\nfBKKXnjTqgg8dCg0bAj/93/w0UdQXBzyNkq5iWsDCATugXiOe/dSwpmFFU4xRe2BJFGrVlY1xnXr\noG5dqlWDW2+14sr48fDD94Y19/yr/PzDh+Hdd6FrV6sno5Qq49qPsGBDWJ7jgXYkDNTbCLecuybR\nU0DpbC0PEatAYwczG14PsKHVddcloGFKVR6u7YEEG8IypmLJEo9wZmHpQsJK7ptv/P7j+Dkjn4Kh\n5/DSS7DDezPlHj2sfImWlFcu49oAAsF7IJ5ch78dCT3HtZRJGrrnHvj5Z2u4qkGDssNNR/Thb3cJ\nM2daVVSuugqmTYODC5fClClWvqRFCzj1VLjvPqtksFJpzrUBJJweSLAhLF1ImMaaN4fHHoP16+E/\n/4FrriGrX2+6dYN33oG1a+GKK+CFF+Cp86ZW/NnFi+GRR+Duu5PTdqUSyNUfYYF6IN5Vd+3Hs7LC\nn8ars7AquawsuPxy62FTsyb07Ws9Drd6B1b4/mhRQVdy/V3zhx9g1y6rdleVKnFptlKJ4uoeSLAh\nrEABJJxZWOEuJNQkeiW3ciXZK5b6fenSsZfTrh384x/WjOCyf0dPP23tslinDlx2GTzxBMyfrzsp\nqkrJ1QEkkEA5EO/j4eRAPOfb6RBWmmjZElasgIcfthYeejRrxhdbTmL0aCvhfu21cPzx0LePYe+H\ns6xz9u6FTz6xKj2efba1klGpSsa1AQRC90CClXP3zoEEmoXlr5ehQ1hppGVLGDbM2qB9/XorMXLv\nvVSpKlx0kbUwccUK+Pxz+GPeCqpvX+9zCZOZCe3a+b/+++/H+Q0o5ZxrP8LiMYTlrwfiL4BoDyRN\nNWoEAwb4HBaxFsLn15/l98cWmNO4/Y/VufhiaN/eSo/UqgVs2wa//BLnRivlnGt7IMGGsLyn63qE\nUwvLOweiPRBVpmVLazyrTp0Kh9ve2oFRo6x/P489Zg13nXQSPN1zHh9tacfSpRX/HSqVKlz9EeZ0\nCCvcHEigISxP4NAkustceqn1KCmxMuuzZsHs2WR1voSLLoKLLrJOKy62ZgMfGv4t49fezRdXWJ2R\ndu3g3HPhrLPglFOgcePQNd2UiifXBpBQQ1j2ooke9iEsf+tA/M3C0h6I8pGRYW2VePrp1sJFL1lZ\npTn5k/fRbrRVcmXzZmuR/DffWBO5Fi2CgwetQOL9qFEjwe9HuZZrP8KclDIJNoTlb0MpDSAqKn36\nlH1Zvz5062Y9PLZssXoqixZZM4EnTIAlS+CYY6xA0qZNae6l9HHMMdpjUbHl6o+wcBYShrsjYaAC\ni5pEV461ahX05WOOocLQF1j//latsgLLTz/BF1/ASy/B8uXWv097QPE8Wra0tgBWKlKu/QgLZxaW\nkx0JvavxBuuBaA5ExVpmphUQWrb0fW3bNmtKsefx7rvW88qVUL26VePL3yMvD44+OpHvQlUWrg4g\ngQTKgXgf95Q2gYo9EJ3Gq1JR3bpWEv7ccyseN8bKsaxeXf748UdrD63Vq2HNmooBpnFja8Zyw4bW\ns+dRrVqi35FKNld/hIWzDiQeQ1iaA1GpRMQqPNyggf/1jN4BZv162LABvv3Wet6wATZutAKIPaA0\nbGg9PNdu0MDK5dSoobmYdOHajzAn1XjDqYXlr5SJN+2BqMokVIAB6//Jtm1WIPEElQ0brNqRmzfD\npk3lz4cPW4HEE1C8n+vVs3pLdepYzxpwUperP8Iircbrb0fCYKVMdBaWcgsR64O/Xj1rS5Rg9u2z\ngol3YFmzBubNswKR57F9O+zfD7VrW8HEHljsX9eubVVJrlXLevZ8fdRRGnziybUfYU5LmXjnQPxt\naasLCZUKLCenPJ8SjkOHrEDiCSjeAebXX2HnTuuxa1fFZ2PKg4r3c82a1uwz78fRR/seq17d/2iC\n27k6gAQSakdCLWWiVOJUqQLHHms9InXggG9Q2bWr/OuiImvYbc8e62vvh+f4vn1W4PMElLvv9lv2\nzHVS5iNMRDoDT2PV53rFGDPGzzljgS7AXqCPMeb7aO4Z6TTecGph6UJCpVLHUUdZD9vuxI4cOWJV\n4PcEltq1Y9O+yi4lOmUikgGMAzoBbYAeInKi1zldgBbGmJbAAOCF6O7pfAgr3A2lgg1heQJUpAoL\nCyP/oUpE31/llc7vLTMTFiwopFEjOPHE6ANSukiJAAKcDaw0xqwxxhwGpgDdvM7pBkwCMMbMBWqK\niOO/RifVeL2T69H0QJz2PtL5Pyno+6vM0vm9Qfq/PydSJYA0AtbZvl9feizYORv8nBORSIsphlML\nK5wcSFaWJtCVUpWfa0fhMzOtjeRmzPB9bd++8l7C1q3Qtat1fMEC6NLFChrFxdZK3VNOsV7LyoJJ\nk6y5757eRXa2lQD0lp0NVavG530ppVSiiAn0a3giGyFyDjDCGNO59PuhgLEn0kXkBeC/xpi3S79f\nBnQ0xmzyulby35BSSlVCxpiIVs2kSg9kPnCCiOQBvwHdgR5e53wI3Aa8XRpwdnoHD4j8D0AppZQz\nKRFAjDFHRGQgMIPyabw/icgA62Uz3hjzsYhcJiI/Y03j7ZvMNiullNulxBCWUkqpyidVZmHFhIh0\nFpFlIrJCRIYkuz2xJCKNReRzEVkiIotF5I5ktynWRCRDRBaIyIfJbkusiUhNEZkqIj+V/h0GKEtY\nOYnI30TkRxFZJCJvioif6SOVh4i8IiKbRGSR7VhtEZkhIstF5FMRqZnMNkYjwPt7rPTf5/ci8q6I\nhNwcOW0CSDiLESu5YuAuY0wb4FzgtjR7fwCDgKXJbkScPAN8bIxpDbQFfkpye2JGRBoCtwOnG2NO\nxRoa757cVkVtItZnid1QYJYxphXwOXBvwlsVO/7e3wygjTHmD8BKwnh/aRNACG8xYqVljPndU7rF\nGLMH6wMoqnUwqUREGgOXAS8nuy2xVvqb3PnGmIkAxphiY8zuJDcr1jKB6iKSBeQAG5PcnqgYY2YD\nO7wOdwNeK/36NeDKhDYqhvy9P2PMLGOMZ+n0HKBxqOukUwAJZzFiWhCRpsAfgLnJbUlMPQUMBtIx\nKdcM2CoiE0uH6MaLSNrs32eM2Qg8AazFWuC70xgzK7mtiov6npmfxpjfgfpJbk889QM+CXVSOgUQ\nVxCRo4FpwKDSnkilJyJ/AjaV9rCk9JFOsoDTgX8ZY04H9mENh6QFEamF9dt5HtAQOFpEeia3VQmR\njr/sICL3AYeNMZNDnZtOAWQD0MT2fePSY2mjdHhgGvC6MeaDZLcnhtoDV4jIKuAt4EIRmZTkNsXS\nemCdMebb0u+nYQWUdHEJsMoYs90YcwR4DzgvyW2Kh02e+nsiciywOcntiTkR6YM1lBzWLwDpFEDK\nFiOWzgDpjrX4MJ1MAJYaY55JdkNiyRgzzBjTxBjTHOvv7XNjTO9ktytWSoc91olIfumhi0mvyQJr\ngXNE5CgREaz3lw6TBLx7wx8CfUq/vhGo7L/EVXh/pVtqDAauMMYcDOcCKbGQMBYCLUZMcrNiRkTa\nA38GFovIQqzu8zBjzPTktkyF6Q7gTRHJBlaRRgthjTHzRGQasBA4XPo8Prmtio6ITAYKgLoishZ4\nEHgUmCoi/YA1wHXJa2F0Ary/YUAVYKb1ewBzjDG3Br2OLiRUSinlRDoNYSmllEogDSBKKaUc0QCi\nlFLKEQ0gSimlHNEAopRSyhENIEoppRzRAKKUUsoRDSBKKaUc0QCiVAKJSNqsQFdKV6IrlSAi0gNr\nr4yqWGUiFiS5SUpFJW1qYSmVykTkWWCxMaZS14hSyk57IErFmYi0Ab4DqpeWO1cqLWgORKn4awns\n1uCh0o32QJSKMxGpD/yItdHSIuB3Y8x7yW2VUtHTHohS8ZcJfAzswso7Hkpuc5SKDe2BKBVHIlIV\neBPoZ4zZnez2KBVL2gNRKr5uAd7X4KHSkU7jVSq+TsAawlIq7egQllJxJCKnAZ8AXwIvGmM+S3KT\nlIoZHcJSKo6MMQuBNsBPwAQR+VlECpLbKqViQ3sgSiWQiPQHOhtjrkl2W5SKluZAlEoAEWkLdACa\nAM8kuTlKxYT2QJRSSjmiORCllFKOaABRSinliAYQpZRSjmgAUUop5YgGEKWUUo5oAFFKKeWIBhCl\nlFKOaABRSinlyP8HYMajjOD9NU0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3b6809b850>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "compareplot1Dprofile(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reference: \n",
    "\n",
    "[1] Stahler, Steven William., and F. Palla. *The Formation of Stars*. pg. 243-248. Weinheim: Wiley-VCH, 2004.  \n",
    "[2] Typanski,Nathan, *Solving a second-order ODE with NumPy and SciPy*. Online. [Link](http://www.nathantypanski.com/blog/2014-08-23-ode-solver-py.html)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}