{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Foster and Chevalier setup " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Foster and Chevalier (1993) sets up an isothermal sphere that has constant external pressure on its outer edge. The density on the outside is decreased so that the sound speed in the outer region is 1000 times larger than the sound speed inside the inner region. \n", "\n", " $$P =\\rho C_s^2 \\rightarrow \\rho_{out} = 10^{6} \\rho_{edge}$$ \n", " \n", " where $\\rho_{edge}$ is the density at the edge of the cloud. Since the Lane-Emden solution is monotonically decreasing, this is simply the minimum density of the sphere, which is $7.31\\times 10^{-22}$ for our setup. By the ideal gas law $T_{out}=10^6T_{in} = 10^7K$. " ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from yt.units import g,cm,K,s\n", "from yt.utilities.physical_constants import kboltz , mp,G" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "rho_out = 7.31207875416e-16 g/cm**3\n" ] } ], "source": [ "rho_edge = 7.312078754161586E-022*g/cm**3\n", "print \"rho_out = \",rho_edge*1e6" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Therefore, to find P_out\n", " $$P_{out}=P_{edge} = \\Bigg(\\frac{\\rho_{edge} k T_{in}}{m_p}\\Bigg) $$" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "P_out = 6.03166654669e-13 erg/cm**3\n" ] } ], "source": [ "T_in = 10*K\n", "print \"P_out =\", rho_edge*(kboltz*T_in/mp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The inside density $\\rho(r)$ is simply given by the Lane-Emden equation times $\\rho_c =1.1\\times10^{-19} g/cm^3$ corresponding to the initial cloud density in Larson (1969). The inside pressure is obtained by ideal gas law using the observational value for $T_{in}=10K$, \n", "$$P_{inside}=\\Bigg(\\frac{\\rho kT_{in}}{m_p}\\Bigg)$$" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "P_in = rho*8.25e+08\n" ] } ], "source": [ "print \"P_in = rho*%.2e\"%(kboltz*T_in/mp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the linear interpolation part of the code, we compute the values using the non-dimensional radius scheme: \n", "\n", "$$\\xi =\\Bigg[ \\Bigg(\\frac{\\sqrt{4\\pi G\\rho_c}}{a}\\Bigg)\\Bigg]r $$ where a = 28731cm/s. is sound speed (computed from $T_{in}$=10K)\n", "\n" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "rc = 1.057e-17 rr\n" ] } ], "source": [ "rho_c = 1.1e-19\n", "c_s = 28730*cm/s\n", "print \"rc = %.3e rr\"%(sqrt(4*pi*G*rho_c)/c_s)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Sphere Fattening\n", "\n", "Using the Foster and Chevalier setup, we find that the code takes a long time to track the full evolution of the sphere, because the extremely high temperature outside ($10^7K$) cause the sound speed to be very high, so the code has to be more careful by taking smaller timesteps to resolve fast-moving waves. So, we decided to make outside temperature cooler to speed up the calculation. Since we still have to maintain the condition that $P_{edge} = P_{out}$, so that implies that if the outside temperature decreasess then the density must increase accordingly. This is why we call the process \"fattening\", as we multiply the density everywhere in the simulation by a ``fattening_factor``." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "%pylab inline\n", "import yt\n", "yt.mylog.setLevel(50)" ] }, { "cell_type": "code", "execution_count": 77, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "After fattening, edge of sphere is: 7.30195258565e-20\n" ] } ], "source": [ "# Loading in the non-dimensionalized density values from solving the Lane-Emden equation\n", "D = np.loadtxt(\"../data/density.txt\")\n", "\n", "xi_max= 16.9 # radius of sphere\n", "rho_c = 1.1e-19 #central density of sphere\n", "fattening_factor =100 #factor to enhance the density everywhere by\n", "\n", "rho_edge = rho_c*dens[int(xi_max*100)] #physical density at the edge of the sphere (1 xi =100 index array values)\n", "print \"After fattening, edge of sphere is: \",fattening_factor*rho_edge" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Outside the sphere, the density should drop by a factor of 150 (rather than $10^6$ as in the FC case): " ] }, { "cell_type": "code", "execution_count": 80, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Density outside of the sphere is: 4.86796839043e-22\n" ] } ], "source": [ "rho_out = fattening_factor*rho_edge/150.\n", "print \"Density outside of the sphere is: \",rho_out" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulation Results:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is a radial density plot that illustrates how the fattened sphere setup looks like: " ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Loading in data from the simulation and creating a spherical radial profile \n", "pf = yt.load(\"../../data_astroSim/data/sphere_hdf5_chk_0000\")\n", "sp = pf.sphere(pf.domain_center, (1.5,\"pc\"))\n", "rp = yt.create_profile(sp,'radius','density')" ] }, { "cell_type": "code", "execution_count": 76, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 76, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEZCAYAAACjPJNSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd0VNX2wPHvDk0SkkASOiGAiigI2ECahmehSH0qSAmC\nXdEnKogokkQQsYCKKBZ40osIgoCCNRQR4flTUKQLoSM9NEPJ/v1xkzEJCWmTmUmyP2vNkjlzy567\nYnbOueeeLaqKMcYYk9/8vB2AMcaYosESjjHGGI+whGOMMcYjLOEYY4zxCEs4xhhjPMISjjHGGI+w\nhGOMMcYjLOEYY4zxiOLeDsBdRKQm8AIQpKpdktsEGAoEAatVdbIXQzTGmCKt0PRwVHWbqj6Qrrkj\nUA04A+zyfFTGGGNS+FzCEZHxIrJfRNama28tIhtEZJOIDMzm4a4AflDV/sBjbg/WGGNMtvlcwgE+\nBlqlbhARP2BMcntdoJuI1En+LEpERolI5ZTNU+26EziS/O/z+Rq1McaYi/K5hKOqy/knSaRoBGxW\n1XhVPQvMwBkuQ1Unq+rTQKKIjAUapuoBfQa0FpG3gSWe+QbGGGMyUlAmDVTF6a2k2IWThFxU9TDw\naLq200D6+zppiIgtl22MMbmgqpL1Vv/wuR6ON0RHR/P999+jqkX6FR0d7fUYfOVl18KuhV2LjF/f\nf/890dHRufpdW1B6OLuB6qneV0tuc4uYmBh3HcoYYwq1yMhIIiMjiY2NzfG+vtrDEdLe/F8NXCYi\nESJSErgH+NxdJ4uJiSEuLs5dhzPGmEIrLi4u13+ki6pv3cIQkWlAJBAK7AeiVfVjEWkDvIWTJMer\n6gg3nU997Rp4S1xcHJGRkd4OwyfYtfiHXYt/2LX4h4igObyH43MJx9NERKOjo13dRGOMMZmLi4sj\nLi6O2NhYSzg5ZT0c4yk1atQgPj7e22EYkyMRERFs3779gnbr4eSCJRzjKcn/g3o7DGNyJLOf29wk\nHF+dNOBRNmnAGGOyp1BNGvA06+EYT7EejimIrIdjjDGmwLGEgw2pGZMXr7zyCg899JDbt82Kn58f\nf/75p1uO5Ul9+vRhyJAhud6/Xr16LF261I0R5YwNqeWBDakZTykIQ2oTJkxg1KhRbN26leDgYDp1\n6sQrr7xCcHCwt0O7QLFixdi8eTO1atW64LOWLVsSFRXFfffdd9FjTJw4kfvvvx9/f38AVBURYdOm\nTVSqVClf4u7Tpw/h4eG89NJLbt02v9iQmjHG7UaOHMmgQYMYOXIkCQkJrFy5kvj4eG677TbOnTuX\n4T7nz3uv6oe7knfTpk1JSEggISGB48ePk5CQkG/JpqizhIMNqRlz/PhxYmJiGDNmDLfddhvFihWj\nevXqfPLJJ2zfvp0pU6YAEBsby913301UVBRly5Zl4sSJxMbGEhUV5TrWpEmTqFGjBuXLl2fYsGHU\nrFmT7777zrV/yrbx8fH4+fkxadIkIiIiqFChAsOHD3cdZ/Xq1TRt2pRy5cpRtWpVnnjiiUwTnzcc\nOnSI9u3bU65cOUJDQ7n55ptdn23YsIGWLVtSrlw5rr76aubPn5/hMSZOnEiLFi3StKUMFX700UdM\nnTqV1157jaCgIDp27AiQ5nqeOXOGfv36UbVqVapVq8ZTTz3F2bNnAViyZAnh4eGMGjWKihUrUrVq\nVSZMmJDn752XITVLODgJx1YZMEXZihUrSExMpHPnzmnaAwICaNu2LV9//bWr7fPPP6dLly4cPXqU\n7t27A87wCsAff/xB3759mT59Onv37uXYsWPs2bMnzTFTtk3xww8/sHnzZr755hteeuklNm7cCDhD\nZm+99RaHDx/mxx9/5LvvvuO9995z+3fPrZEjRxIeHs6hQ4f466+/XMny3LlztG/fntatW3PgwAFG\njx5Njx492Lx5c4bHSX89Ut4/+OCD9OjRg2effZaEhATmzZt3wb7Dhg1j1apVrF27ljVr1rBq1SqG\nDRvm+nzfvn0cP36cPXv2MG7cOPr27cuxY8fy9L0jIyMt4RhTIInk/ysbDh48SFhYGH5+F/5KqFy5\nMgcPHnS9b9KkCe3btwfgkksuSbPt7Nmz6dChA02aNKF48eJZ3nsQEWJiYihZsiT169enQYMGrFmz\nBoBrr72WRo0aISJUr16dhx56iCVL3F9H8ccffyQkJISQkBDKlSvH5Zdfnq39SpQowd69e9m2bRvF\nihWjWbNmAKxcuZKTJ08ycOBAihcvTsuWLWnXrh3Tp0/P1nFzMlQ4bdo0oqOjCQ0NJTQ0lOjoaCZP\nnuz6vGTJkrz44osUK1aMNm3aUKZMGVdC9wZLOMZ4k2r+v7IhLCyMgwcPkpSUdMFne/fuJSwszPU+\nPDw80+Ps2bMnzeelS5cmNDT0oueuWLGi69/+/v6cOHECgM2bN9O+fXsqV65M2bJleeGFF9IkPndp\n0qQJhw8f5vDhwxw5ciTTnkh6zz77LJdeeim33347l112Ga+++ipw4TUAZ3mY3bvdVlHFZc+ePVSv\n/k/lloiIiDQ9ytDQ0DR/RKS+vt5gCccYQ5MmTShVqhRz5sxJ037ixAm+/PJLbr31Vldb+iGg1CpX\nrsyuXbtc70+fPs2hQ4dyFdOjjz7KlVdeydatWzl69Cgvv/yyT83yCwgI4I033mDr1q18/vnnjBo1\niu+//54qVaqwY8eONNvu2LGDqlWrZniMU6dOud7v27cvzecXu9YAVapUSbM+X3x8PFWqVMnN1/EI\nSzjYpAFjgoKCGDJkCE888QSLFy/m3LlzbN++na5du1K9enV69uyZrePcddddzJ8/n5UrV3L27Nks\nx/ovlkCOHz9OUFAQ/v7+bNiwgbFjx+bkK3H27FkSExNdr8wmHOQ2iS1cuJCtW7cCEBgYSPHixfHz\n86Nx48YEBATw2muvce7cOeLi4liwYAHdunW74BgNGjRg3bp1rF27lsTERGJjY9MkmYoVK170WaNu\n3boxbNgwDh48yMGDBxk6dGiaCRz5wSYN5JFNGjAGBgwYwPDhw+nfvz/BwcE0adKEiIgIvvnmG0qU\nKJGtY1x11VW88847dO3alSpVqhAUFESFChUoVapUhttndsMc4I033mDq1KkEBQXx8MMPc88991x0\n3/Qee+wx/P39Xa+UZ3ICAwP54YcfXNutXLmSoKAggoKCCAwMJCgoiJ9//hmAtm3bMmJExqW3Nm/e\nzK233kpgYCDNmjWjb9++3HzzzZQoUYL58+fzxRdfEBYWxuOPP87kyZNd94ZSx3355ZczZMgQbrnl\nFmrXrn3BjLX777+fdevWERISwr///e8L9h88eDDXX3+96/7X9ddfzwsvvJDpNcnqmmVHXiYN2IOf\n9uCn8ZCC8OCnu508eZKyZcuyZcsWIiIivB2OyQV78NMY47MWLFjA6dOnOXnyJM888wz169e3ZGMA\nSzjGGDebN28eVapUoVq1amzdupUZM2Z4OyTjI2xIzYbUjIcUxSE1U/DZkJoxxpgCxxIONi3aGGOy\ny8oT5IENqRlPsSE1UxDZkJoxxpgCxxKOMcYYj7CEY4wpVF588cUsK336oq1bt2a4Wnd2TZo0iTvu\nuMONEbmfJRxjDAA1atTA398/zRIv6ReTzEhUVNQFZQjCw8NZunRpfoWaa99++y01a9bM1rbVqlW7\n4Ho8/fTT+RpfdpeeySg59erVi4ULF+ZHWG5T3NsBuIuI1AReAIJUtUtyW3OgB873vFJVm2e077p1\nUL06BAZ6LFxjfI6IsHDhQlq2bOntUPJVdn+piwiLFy++YH0zX6CqblkXzdMKTQ9HVbep6gPp2par\n6qPAAmBiZvvWqwdly8J118HLL0M2/qgzxi1iYjKumZbZrNOcbp9TGc1GUlXuvvtuKleuTEhICC1b\ntnQV8Ro7diwzZ85k+PDhBAUFceedd9K9e3f27NlDmzZtCAoK4q233gKcyp5NmjShXLlyXHvttSxb\ntsx1jhYtWhATE0OzZs0ICgqibdu2HD161PX5xfbdtm0bN910E8HBwbRp0ybX5RCyez2yY/z48dSs\nWZOgoCAuu+wyPvnkE9fxXnrpJWrUqEGlSpW47777Mq1Pk76XmHqoMKWcderFRsePH5/mj4Xly5dz\nww03UK5cOW688UZWrVrl+iyr651vVNWnXsB4YD+wNl17a2ADsAkYeJH9P8mgbSYQkMn2F1SsKlFC\n9amnVI8eVWPcxvnfLa3o6IyrpkVHZ3yMnG6fEzVq1NBvv/32gvakpCSdOHGinjx5UhMTE/WJJ57Q\n66+/3vV5z549NTY2Ns0+1apV06VLl7re79y5U0NDQ/Xrr79WVdXFixdrWFiYHj58WFVVmzdvrrVr\n19atW7fq6dOntUWLFvriiy+qquqOHTsuuu8NN9ygAwcO1DNnzmhcXJyWKVNG+/Tpk+F3/Oabb7Rm\nzZrZuh7VqlXTJUuWZGvb1BISEjQ4OFi3bt2qqqr79u3T9evXq6rqBx98oFdccYXGx8friRMntGPH\njq5Yt2zZon5+fpmef/DgwZluq6o6btw4bdmypaqqHjhwQIODg3XmzJl6/vx5nTx5soaGhurR5F9q\nF7ve6WX0c5uqPUe/332xh/Mx0Cp1g4j4AWOS2+sC3USkTvJnUSIySkQqp2yebt9w4KiqnsxuAGfP\nwrhx8PffefgWxhRAnTp1cpVbTr0cfq9evfD396dkyZIMGTKEn3/+mdOnT1/0WJqqdzBp0iQ6duzo\nKuR2++2306BBAxYtWuTa5v7776dWrVpccskl3H333fz6668ATJ48OdN9t23bxtq1a4mJiaFEiRLc\nfPPNtG3b1m3Xo127dq7S0yEhIUycmOlASRp+fn789ttvJCYmUrFiRerUqQM4JaH79+9P9erVCQgI\nYPjw4UybNs1t8aZYsGAB9erVo0uXLvj5+dGzZ09q1aqV5h5PZtc7P/lcwlHV5cCRdM2NgM2qGq+q\nZ4EZQMfk7Ser6tNAooiMBRqKyMBU+96Pk8Ry5PnnIVXlW2OKhHnz5rnKLadU/0xKSnKVUy5btiyX\nX345IpKjcs/x8fFMmzbNlczKlSvHTz/9xN69e13bVKpUyfXv1KWQM9t3z5497Nmzh9DQUC655BLX\nvu5cmXrhwoWu0tOHDx/m3nvvzXKfwMBApk+fzpgxY6hUqRIdOnRgy5YtgFMSOnV8ERERnDlzhgMH\nDrgt5ozOk3Ku1GWuM7ve+amgTBqoCuxM9X4XThJyUdXDwKPpd1TVmKwOHhwcw7FjKe8iqVEjkn79\ngEWLoFkzm01giozUvZIUkyZNYtGiRcTFxREeHs6hQ4coX768a9uMbl6nbwsPD+e+++7j3XffzXFM\nF9v3zz//5NChQyQmJrqKvO3YsQN/f/8cnycjGV2P7GjVqhWtWrUiMTGR5557jocffphvv/02w5LQ\npUqVonz58iQkJKQ5xsXKT2en9HT6GWs7duygc+fOufo+4Cxpk9clwHyuh+MNR4/GsG1bDM89F4O/\nfySvvQaXXAJ8+ik88YS3wzOFWExMRndkLj5pICfbu8Px48cpVaoU5cqV4+TJkzz//PNZlkGuVKlS\nmraoqCg+++wzvvnmG5KSkvj777+Ji4vL9rTrzPatVasW9evXJyYmhrNnz7J06dIspwYnJSWlKT2d\nmJiYwytycfv27XPVBCpevDgBAQGuKczdunVj1KhRxMfHc/z4cQYPHkz37t1d+6ZOcA0bNmTGjBmc\nP3+eVatWuXqcABUqVEBE2LZtW4YxtGvXjj/++INZs2Zx/vx5pk2bxtatW/P0nE5Kpc+UV24UlISz\nG6ie6n215Da3iImJYfv2OF55BbZuhbvuSv7g7behaVNo3x7WruXYMSf/HD/urjMb4zsy+6u5T58+\nVK5cmSpVqnD11VfTvHnapwseeOABfv31V0JDQ+nSpQsAgwYNYsiQIYSEhDB69GgiIiL47LPPGDp0\nKOXLl6dGjRqMGjWKpKSki54byHLfGTNmsHz5ckJDQ3nllVfo1avXRb/nzp07XWWnS5cujb+/Pzt2\n7GDYsGF07NgxzbYpM+1SXl27dgVgyZIlhISEZHj88+fP8/rrr1OlShXKly/Pjz/+6OqdPfjgg3Tt\n2pUWLVpw2WWXERwc7JrFl/46vPzyy6xfv55y5crx8ssv06NHD9dnZcqUYdCgQTRu3JiQkBD+7//+\nL00MYWFhfP7554wYMYKwsDDefvttFi5cSHBwcJbXOyuFbvFOEakBzFfVq5PfFwM2ArcAe4FVQDdV\nXe+Gc2l0dDSRkZFERkZmvNEnn5D0nye5t/yXTPm9Ic2bO6NtAQF5PbspSmzxTlMQpf+5TRlai42N\nzfHinT6XcERkGhAJhOJMj45W1Y9FpA3wFk6vbLyqjnDT+TSra3D0KLx2/Sc8srU/DVjDUcpxxx0w\ndy4ULyh3wYzXWcIxBZE7V4v2uYTjaVklHFW46SZYvhxG8jTX8AvtWMApAnjwQfjgA+fBO2OyYgnH\nFERWnsDNLlaATQQGDYKSJWEAr7ONmnxOB0rxNx99BOPHezZWY4zxpkJ3D8eTsluAbfZsuPtuED3P\nFHpShhOMbzuHqZ+UsHs5Jlush2MKIuvhuFl2SkzfeSe89RYkUYxeTCIiXJlTvZ8lG2NMkWI9nDzI\naYnpZ591ngMd/MQx5JqGThZKN5XSmIxYD8cURDZpwI1ymnBUU00S+PFHJ9ksXw61a+dPgKbQsIRj\nCiIbUnOz7AyppUgzI61JE6eeQceOcOQIZ8/ClClOUjLGmMLIhtTyIKc9nAw99RRnV//KHcUX8/WS\nkrz5Js5abMakYj2cjC1fvpwHH3yQ9evz/Bx3hmrWrMn48eP517/+lS/HT61ly5ZERUUVyBLXmbEe\njo/5vfcbxP0SxL+X/AeAZ56Br7/2clDG5FBKieng4GBCQkJo3rw5H3zwQb4nyebNm6dJNjVr1uS7\n777L13MWBA8//DB16tShWLFiTJo0Kc1nEydOpHjx4mnKX6cu1nbkyBE6d+5MmTJlqFmzJtOnT/d0\n+BmyhJNHO3dCk+bFuPPUZG5mCfczjqQk6NoVklckN6ZASCkxfezYMeLj43nuued49dVXuf/++70d\nWr45f/68t0PIVMOGDRk7dizXXXddhp83bdqUhIQEjh8/TkJCAjfddJPrs8cee4xLLrmEAwcOMGXK\nFB599NF860HmhCUccnYPJ73wcOjdG44TRCfmMpznacxKjhyBDh0g3YrjxqQRExeDxMoFr5i4GLds\nn1MpvZnAwEDatWvHzJkzmThxIn/88QcAZ86coX///kRERFC5cmUee+wx12rLS5YsITw8nFGjRlGx\nYkWqVq3KhAkTXMf+4osvqFu3LkFBQa7tUu8H0KtXL3bs2EG7du0ICgrijTfeoF27dheUJmjQoAHz\n5s3L8DtMnjyZGjVqUL58eYYPH57ms9jYWO6++26ioqIoW7YsEydO5MyZM/Tr14+qVatSrVo1nnrq\nKc6ePZsmtldeeYXy5ctTq1atLAumbdmyhcaNGxMcHEznzp1dpZtz+j0effRRWrZs6Sq7kF2nTp1i\nzpw5DBs2jNKlS9OsWTM6duzI5MmTc3SczOTlHo4lHJyEk+nCndkwahS0bAmbuIIHGMen3EVF9hER\nAckL2hpTIN1www1Uq1aNZcuWATBw4EC2bNnC2rVr2bJlC7t37+all15ybb9v3z6OHz/Onj17GDdu\nHH379uVYcrGpBx54gI8++oiEhAR+//33NPdUUlYvnjRpEtWrV2fhwoUkJCTQv39/7r333jS/LNes\nWcOePXsyXGr/jz/+4LHHHmPq1Kns2bOHQ4cOpSk6BvD555/TpUsXjh49Svfu3Rk2bBirVq1i7dq1\nrFmzhlWrVjFs2LA03+nw4cPs2bOHCRMm8NBDD7F58+ZMr9nkyZOZMGEC+/bto1ixYvznP85Qe06+\nR3b88ssvVKhQgTp16jBs2DDX6tmbNm2iRIkSXHrppa5tGzRowLp163J1nvRSyhTkhiUcNyhRAmbN\ngpo1YT4dGMcDrKjWlQWfnaVsWW9HZ0zeVKlShcOHDwPw0Ucf8eabbxIcHExAQADPPfdcmvsDJUuW\n5MUXX6RYsWK0adOGMmXKsHHjRtdn69at4/jx4wQHB9OwYcNMz5n6vlGHDh3YvHkzW7duBWDKlCl0\n7dqV4hmsnDt79mzat29Ps2bNKFGiBEOHDr1gKf4mTZrQvn17AC655BKmTZtGdHQ0oaGhhIaGEh0d\nnSYxiAhDhw6lRIkS3HTTTdxxxx188sknmcYeFRXFlVdeSenSpRk6dCgzZ85EVXP0PbJy88038/vv\nv/PXX38xe/Zspk+fzuuvvw7AiRMnCAoKSrN9UFAQx32groolHDcJDYV58yAsDC6bOIRa9fwpNniQ\nt8MyJs92795NSEgIBw4c4NSpU1x33XWucs9t2rTh0KFDrm1DQ0NdxcYgbeni2bNns3DhQiIiImjZ\nsiUrV67M1vlLlSpF165dmTJlCqrK9OnTiYqKynDbPXv2uIbnUs4fGhqaZpvUn6fsU736P+W2IiIi\n2LNnj+t9uXLlLihhnfrz9FIfPyIigrNnz3Lw4EFKlSpFly5dsvU9slKjRg1XCem6desyZMgQPv30\nU8CplZO+euixY8cI9IHKxZZwyNs9nNSuvhq2bYOevfycB3LmzIGpU/MeoDFesnr1avbs2UOLFi0I\nCwvD39+fdevWcfjwYQ4fPszRo0ddQ2ZZue6665g7dy4HDhygY8eOrmJt6WVUHKxXr15MmTKFb7/9\nloCAABo3bpzhvpUrV2bnzn+q0Z86dSpNQszo+FWrVr2g7HOVKlVc748cOcLp06dd73fs2JHm8/RS\nnz8+Pp6SJUsSFhYGOMNq2fkeuZHSK6xduzbnzp1z9aTAGb6rW7euW86Tl3s4qGqRfjmXIJ+sXasa\nFqa6apWqqv79t+qhQ/l3OuPb8vVnzQ1q1Kih3377raqqJiQk6Pz58/XSSy/V3r17u7bp16+fdunS\nRf/66y9VVd21a5cuXrxYVVXj4uI0PDw8w2OeOXNGp06dqseOHVNV1XHjxmmNGjUy3K9Jkyb60Ucf\nXRBf7dq1tX79+jp06NBMv8O6des0MDBQf/jhBz1z5ow+88wzWqJECdf3iomJ0aioqDT7DB48WJs1\na6YHDhzQAwcOaPPmzXXIkCGu2IoXL64DBgzQM2fO6NKlS7VMmTK6cePGDM8fGRmp4eHhun79ej15\n8qTefffd2rNnzxx/D1XVM2fO6OnTp7VZs2b60Ucf6d9//61JSUmqqvrll1/q/v37VVV1/fr1Wq9e\nvTTH69atm3bv3l1Pnjypy5Yt07Jly+off/xx0fNlJrOf2+T2nP2+zekOhe2V778E5s5VrVpV/1q1\nTZs0UW3eXDUxMX9PaXxTQUg4/v7+GhQUpGXLltWmTZvq2LFjXb/kVFUTExP1+eef11q1amlwcLBe\nddVV+s4776hqxgmnZs2aroTTunVrDQkJ0eDgYG3UqJGuWLEiw/3mzZun1atX13LlyunIkSNd7cOG\nDVM/Pz/dtm3bRb/HpEmTtHr16hoWFqbDhw93xaCaccL5+++/9cknn9TKlStrlSpVtF+/fpqY/D9p\nSmzDhw/XsLAwjYiI0KlTp2Z67pYtW+rzzz+vjRo10uDgYO3YsaMeSvdXZna/R2RkpIqI+vn5uV5L\nlixRVdX+/ftrxYoVtUyZMnrppZdqTEyMnjt3zrXv4cOHtVOnThoQEKARERE6Y8aMi57rYtyZcGyl\nAXesNJCFHf1Hk/j2WBqf+4EjhPDQQ/D++1a4raixlQbyZvLkyXz00UdpHnDMb0uWLCEqKoodO3a4\n7Zje+B55YSsNFCCzZsEV7/6Hz861Zx4dKUkiH34IY8Z4OzJjCo5Tp07x3nvv8fDDD3s7lDwpLN8j\ntyzh5LPy5eHcOXiOEfxFBd7nEUB58kmYO9fb0Rnj+7766isqVKhA5cqV6datm7fDybXC8j3ywobU\nPDCk9t570Lcv+HOS5TRnKj0YSX/69YM338zXUxsfYkNqpiBy55Bazp84KoRSVhrIy2oDF/Poo7B2\nLXzwQQAd+JyV3EiTqMv59ygr3GaMKVji4uJy/RiJ9XA80MMBZ1itQwf45hv4/MXVtB7dFhYvhmuv\nzfdzG99gPRxTEFnFTzfyVMIBOHHC6ek0bQp8+ik89RSsWOGsAGoKPUs4piCyhONGnkw4Fxg1Cj74\nAJYuhYoV2bEDKlSAVKtomELEEo4piGxadGHx9NPQowfceisbfjjEjTdCt27O8JsxxhQ2lnC87cUX\n2XttW07d3Ibje48zdy706WNlDYxvmTZtGq1bt86XY/fp04chQ4bkev/AwEC2b9/uvoBMvik0CUdE\naorIOBH5JFVbuIh8ltw+0JvxZWbpMuGKOSNYff4aPqMzJUlkyhR4/HGw0RfjacuXL6dZs2aULVuW\nsLAwWrRowc8//0z37t1ZtGiRt8OjZcuW/Pe//03Tdvz4cWrUqOGdgEyOFJqEo6rbVPWBdM1XA7OS\n2zMvvuFFu3bBiZPCY7zHEcoxlR74cZ6xY+GZZ7wdnSlKjh8/Tvv27XnyySc5cuQIu3fvJjo6OscV\nJ43JjM8lHBEZLyL7RWRtuvbWIrJBRDbloLeyEnhARL4BvP/nWQa6d3ceDE2iGD2ZQjDH+JCH8COJ\na67xdnSmKNm0aRMiQpcuXRARSpUqxa233kq9evWYOHEiLVq0cG3r5+fH2LFjqV27NsHBwQwZMoQ/\n//zT1Tu65557OJd8MzL9vin7//nnnxfEcPToUdq3b0+FChUIDQ2lffv2rtozgwcPZtmyZTz++OME\nBQW5KmmmPlZCQgK9evWiQoUK1KxZk5dfftl17JQ4BgwYQEhICJdeeqlP9NqKEp9LOMDHQKvUDSLi\nB4xJbq8LdBOROsmfRYnIKBGpnLJ5ql37AENU9VagXb5HnkuPPAKvvgpnKEVnPuMKNvL7v54gqqeN\nqRnPqV27NsWKFaN3794sWrSIo0ePpvk8fR2Zr776il9++YWVK1fy2muv8fDDDzNt2jR27tzJb7/9\nlqYSaPpP1QQPAAAgAElEQVR9M6p5A5CUlMR9993Hzp072bFjB/7+/vTt2xeAYcOG0aJFC8aMGUNC\nQgKjR4++4FiPP/44x48fZ/v27cTFxTFp0iQ+/vhj1+erVq3iyiuv5NChQwwYMID7778/F1fK5JbP\nJRxVXQ4cSdfcCNisqvGqehaYAXRM3n6yqj4NJIrIWKBhqh7QIuDJ5PZtnvkGufPss/Dii3CSMmx9\neyFXJqyCAQPsRk5hJ5L/r2wKDAxk+fLl+Pn58dBDD1G+fHk6derEX3/9leH2AwcOJCAggCuvvJJ6\n9epx++23ExERQWBgIG3atOGXX37J9FyZTQ8PCQmhc+fOlCpVioCAAAYNGpTlqsopx0pKSmLmzJmM\nGDECf39/IiIieOaZZ9KUi46IiOC+++5DRLj33nvZt29fpt/PuF9BWdqmKrAz1ftdOEnIRVUPA4+m\na1sH3J3VwVNXr8vPJW6yEhvrrEZw/fXB0HMx/Otf0L8/vPGG6xeHqpU1KFR87A+KK664wnVTftOm\nTfTo0YN+/frRqlWrC7atUKGC69+lS5emYsWKad7v378/x+c/ffo0/fr1Y/HixRw9ehRV5cSJE04t\nlSx+8A8ePMi5c+cuKBe9e/du1/tKlSqliTHl+Km/i8lYXpa0SVFQEk6+ynW5VDcTgeuvT34TEgLf\nfQdt28IDD8AHH/Dl18V57z2YPh3KlPFqqKYIqF27Nr179+bDDz/MMOFkV0BAAKdOnXK937dvX6bb\nvvHGG2zevJnVq1dTvnx51qxZw7XXXutKOBdLOmFhYZQoUYL4+Hjq1KkDOCWeq1atmuvYzT/S/zEe\nGxub42P43JBaJnYD1VO9r5bc5hYxMTF5ztz5IiTEWXxtxw4O3daVbv9OZMECuOUWOHjQ28GZwmbj\nxo2MGjXK1SPYuXMn06dP58Ybb8zTcRs0aMC6detYu3YtiYmJxMbGZpo4Tpw4QenSpQkKCuLw4cMX\n/DFYsWLFDCcbgDN5oEuXLrzwwgucOHGC+Ph43nzzTaKiovIUv0krLi4u13+k+2rCEdLe/F8NXCYi\nESJSErgH+NxdJ0tZLdonlSnDH68t4IcfhNl/tyWQBFatgubNwY1FCI0hMDCQn376icaNGxMYGEjT\npk2pX78+I0eOvGDb7E4CALj88ssZMmQIt9xyC7Vr175gxlpq/fr149SpU4SFhdG0aVPatm2b5vMn\nn3ySWbNmERoaSr9+/S449+jRo/H396dWrVrcdNNN9OzZkz59+mR6vqyG6cyFIiMjc51wslxLTUTa\nASG5OjocVtUFOQpIZBoQCYQC+4FoVf1YRNoAb+EkyfGqOiKXMaU/n0ZHR3v13s3F/P031KkDO+PP\nM4bHacxPtOUL9lOJSpVgwQK47jpvR2myw9ZSMwVR+p/blHs5sbGx7l+8U0RmAe+QtseRXX1VtUsu\n9vMYry7emU1z58I990BiovIiQ7mXibRiMX8FXsaKFVCvnrcjNNlhCccURJ4uwLZcVS8+LzETIuKT\nT/enl98F2PKqUyf46ito314YmjCEfVRiGS3YHjOHevWaeDs8Y0wRYgXY8qAg9HBSrFkDbdrA3r0w\n7+Ev6DD7Xnj3Xeji051Ik8x6OKYg8lh5AhEpLyKPi0j95Pf3iMhCERkpIoWmaovPzlJLp0EDWLUK\nXnsN2o9t68xg698fhg93Pc+RlAS73TZ/zxhj0srLLLWL9nBEZBxwDLgUWIEzNXkNUAGoqqqP5eqs\nPqQg9XAytHs3tG/v3Mh5/32GjvRn1CiYMgXuuMPbwZnUrIdjCiJPFmBbpKrPqGon4JiqPq6qH6nq\ny8CPOTmRySdVq8KyZXD+PAlXN2XCkD85ehTatYMhQ+D8eW8HaIwxjqwmDQSKyGfAvcDMlEYReQ1Y\nm+leBYyvTxrIUkAA24dN4YN677CCJvRmAotow9ChsHw5TJ7s5CXjXREREfbchylwIiIi0rzPt0kD\nIlIcuE5Vf0rXPhCYoarxuTqrDynwQ2rJ7rwT5syBZixnJl0ZxwMM5UXOU5zevSHVgrnGGJNnuRlS\ns1lqhSThHD4MPXrAokVQib1MoheX8DfPVJzKoj+qE5LbR3eNMSYD+XEPJ7snbioiD4pICXccz+Rc\nSAgsXAjR0bBfKtOKxSykHcsTryckbo63wzPGmNwlHBG5RkRGi8hLIhKJM4FgFtDXncF5SkGZFp0V\nPz+IiXEeEq1U2Y9izw+k5KL5Tl2dhx+G48dd227d6iybY4wxOZFv06Iz3UlkKrAMCAda46zevAgo\npar35CoSLyksQ2rpHTwIwcFQogRw7Bj06wdLlsB//8vpxpE0bAjFijn3dho39na0xpiCxmP3cETk\nEVV9P9X72jiJ5wdV/TnHB/SiwppwMrRgATz8MMsr3sntv4zgNP74+cHTTzs9o4AAbwdojCkoPHkP\nJ0lEyqe8UdVNqjq6oCWbIqddO1b/9zfifznMGhrQgqUkJTkFRa+6ylm4wBhj8ktuE84SYIGIPCYi\nddwZkDcUlns4WVGFB54NoSdTGMDrTKUH47mPEA6xY4dzD8gYYy7GG/dwZgN/AhFAM5wHSJcB81R1\ncq4i8ZIiNaQG/PYb9OkDP/8MgSQwlBfpykzm3vgqj6zo5dS5NsaYLHjyHs7DqvpBqveXATcBtVX1\nuRwf0IuKWsIBOHcOXn8dYmMhMRFuDvgfX9d6mBKhQc7q01dd5dr21CkoXdrykDEmLU/ewymderVo\nVd2iqv8taMmmqCpeHAYNcno7t98Od424nhL/9xN07gw33wxPPAGHDgHQt69TznrVKi8HbYwp8HLb\nw2kAvAm8hDMz7ay7A/OUotjDSU3Vebnu3xw86Dw9OmsWO+8dzKVvPMpZnOd5u3eHl1+GGjW8Fq4x\nxkd4ckhtNpAANMa5j7MKZyLBQlVdneMDelFRTziZ0d9+5383P03gkR08w0i+oC0glCgBTz0Fr77q\n7QiNMd7kySG1n4FoVb0KqAGMAUKAIbk8nvEx87bWo9GRxTzDSF5nAEu5iWYs5+zZ5IdJjTEmh7Iq\nT5CZEUB7EamkqquA2cmvAqnAlyfIB+XKwVVXCV/8cQeLaE1PpjCFnmwqXpcmrV4GGno7RGOMF+Rb\neYJMdxKpCNygqguS3xcDrk9fxqAgsCG1zJ07B+PGOYXcDhyAkiQyv92H3P6/4RAZ6Uxzq10bcAq9\nvfaaM+W6UiXvxm2MyX+eHFJ7EpguIi0AVPU8cJmI2J+9hUjx4vDII7BlizOPoPplpWg6/QnYvNkp\nad2sGURFwYYNzJkDzz8PNWvC449DfIGvlGSMcbfc9nAeAyap6ol07f1U9S13BecJ1sPJvvPnnQU/\nXRIS4J130LffZtHZW3nm6GDW4zzDU7y4U59n0CC44grvxGuMyT+e7OFcDmQ0FfpULo9nCoA0yQYg\nKAheeIGv3tvKkqP1+Z6WzKArdfmdc+dg4kTYtMkroRpjfFBuE8504AcR6ZSu6FpEZjuYwuu1sYG8\nynNcylb+x/V8zW0spC1RVb7ljrbWezTGOHKVcJJnpj0DvAEcFZFfRWQD8Js7g8sJEakpIuNE5JNU\nbVeKyEwReVdE7vRWbIXdtGkwcCD4BZbhDQZQk218yl28mfQf/K67BiZNgjNnXNsfPgzPPuvcCjLG\nFB25uofj2llEgOZAFeD/VNXrv0JE5BNV7ZL876eBn1T1BxGZp6odM9je7uG4yZEjMGYMvPWWc79n\n106lzPJFMHIkrF8P//kPPPQQr35YjueSF0G65RZnYkLHjvZ8jzEFSb6sNCAiLwLfq+ryvASX7YBE\nxgPtgP2qWj9Ve2vgLZxe2XhVzfBZ93QJpzzOw6ingSaq2iKD7S3huNnJk7BmDTRtmqrx119h1Ch0\n/nxmnunMyFOP8j9ucH1csaIzBbtdO8/Ha4zJufxKOMFAZ5wyBIeBz1R1Za6jzCogkebACZxZcPWT\n2/yATcAtwB5gNXCPqm4QkSjgGuB1Vd0rIrNU9e50x/QDZqtq5wzOp9qsGVx9tbOSZbt29qd2Pprz\n/l/89OjHPMwHHCaE93mE6XTjFAH8/jvUrevtCI0x2ZEvs9RU9ZiqTlDVB4FXgXrJ90peEZHrcxvs\nRc63HDiSrrkRsFlV45MXCp0BdEzefrKqPg0kishYoKGIDAQQkQgR+QCYCLye6Ulfftl5gPHNNyEi\nAt55B84W2PVIfdpb0yrwGgO5jC28yFDaM58dVOfTyk9QV3+/YHtV+PFH57/GmIIt1/dwkoer7gRu\nAPYCn6rqr24JSiQCmJ+qh3Mn0EpVH0p+3xNopKr/ccO5NDo62vU+slo1ImfOhH37YPp05wFH4zaH\nDjnTpT/44J8p0+HsYGGncVz90zgID4feveGee6BcOVascJ4vrVULevVyXjVrevUrGFMkpV/SJjY2\n1jOrRV9wEJFKwF3AtcBOYJZqBn+uZv94Hk04F1wDVee34oABTqWy3r3zehqTjirExcH778OyZfDn\nn3BJ8XPw9dfOtV+0CFq1YnRCb55edBvnUy37d+ONzqy4Tp28F78xRZ0nH/xMQ1X3qeoYVb0P+BC4\nVUQ+FhF3FWTbDVRP9b5acptbxMTEpF2MTsRJMsuWwdChzmJiNqbjViLQsiXMnJmcbC7BWZ6gTRuY\nMQO2beNMs5Y0+TqWHVTnVZ6lNhsBWLnSKdtjjPG8uLg4YmJicrVvbpe2uRQ4pKpHs9guVFUP5eL4\nNXB6OFcnvy8GbMSZNLAXp/5ON1Vdn9NjZ3Cui89S++svaN8errkG3nsvVaUyk99mznRG1uqwnnuZ\nyH38lxE8xzvFnmL/X0JIiLcjNKbo8mQP5zXgoIj8IiJvi8i/RSQs/Ua5TDbTgBVAbRHZISJ9khcH\nfQL4ClgHzHBHsklxQQ8ntQoVnGGetWudVSmtp+MxFSrAbbfBRrmSQYygEavoykyWh3YkhMPeDs+Y\nIskbPZzBwFwgCGiR/LoJ2IFT+fMDVV2bq4g8LNvP4SQkQKtW0LixM5tNcpTYTR7s3AlTpjgrGmz8\n/QxrWj/HlX/MdobemjTxdnjGFEme7OGcVNXfVXWFqr6qqu1wyk0vBLYBE0WkZS6P7XEX7eGkCAqC\nL7+Eb7+FN97wSFzGER7urDr922/wf7+VJHzWKBg92pk18MYbkJTk7RCNKTK80cMZDSxT1Vnp2u9X\n1fEicgnwiqo+lauoPCjHKw3s2uU8Qj9iBHTvnn+BmazFx0PXrlC+PEyYAKGh3o7ImCLDkz2cF4BH\nReQPEXlDRLqJSDfgZgBV/RvYkMtje1y2ejgpqlWDL76Afv3g++/zNS6ThYgIWLrUKbhz7bWwYoW3\nIzKm0PN4D8e1s0hX4H6cacprgGdVdaeIDALCVfWxXB/cQ3K9ltp33zk9nBUrnKcSjXfNnw8PPujc\n6PnXv7wdjTGFXr6spZbLQB4HglX1Zbcf3M3ytHjnu+/C2LFO0gkKcm9gJudefx127HCWJjLG5Cuv\nPfiZXvJDoD6fbFLkaEgttccec9Zd6dnTblz7ghYtYLlHFjU3psjK1yE1EXlJVYfk6uB52NdT8lye\n4MwZ52GR5s2dRUCN95w540wc2LULgoO9HY0xhVpuejjFs96EBiLSKzfxAPWz3KqgK1kSPv0UGjVy\nFvrs1s3bERVdJUvC9dc7Q5xt2ng7GmNMOtnp4VwHlMnl8U+o6s+53Ncj3FaAbe1ap3zlokVw3XV5\nP57JnSFDnHKj1ts0Jl/lSw/H1xOGO8TExBAZGUlkZGTuD1K/vrP0cefOsGoVVKrktvhMDrRo4Sy4\naozJF+nLFOREvsxSK0jcXmI6Ohq++caZNl2qlPuOa7Ln+HGoXNkpvGPX35h84zOz1Iq06GioWBH6\n9rWFPr0hMBDq1IH//c/bkRhj0rGE425+fjBpEvz0E4wZ4+1oiqYWLZxaRsYYn2IJJz+UKQPz5jk3\nrr/91tvRFD3Nm1vCMcYHWcIhDw9+XkytWs4yK927OyUtjec0b+5MjbaHcY1xO6+tpVYYuH3SQHpj\nxjiz13780bm/YDzjiitg1ixn9qAxxu1s0oAv6tvXKRIWFWV/cXtS8+a2zI0xPsYSTn4TcRb5PHgQ\nctkNNblgEweM8TmWcDyhZEmYPRsmTnSGeUz+S0k4RXzI2BhfYgnHUypWhM8+c1aYXr3a29EUfrVq\nOUvcxMd7OxJjTDJLOOTTLLWMXHstjB8PHTvazLX8JmLDasbkA5ullgf5PkstI++9B2+/7UzdDQ31\n7LmLktGj4fff4cMPvR2JMYWOx2apieMuEYkWkd4iYuUuc+Kxx6BTJ+jQAU6f9nY0hZcVZDPGp+Sq\nhyMi7wKVcBJWQyAA6K+qk9wbXv7zSg8HnCnSPXrA2bPwySfOkjjGvc6fh5AQ2LoVwsK8HY0xhYon\nn8PZrKp3qmpnVa0JtAG6iEj3XB6v6PHzgwkTnOnSTz5ps6nyQ7FizjNQ1ssxxifkNuGEiogrsyXX\nzGkPNHBLVEVFqVLOmmsrVsCLL3o7msLJhtWM8Rm5TThfAnEi8q+UxJM8LrXLbZHlkIh0FJEPRWS6\niNyW3OYvIhNE5AOf7X0FBztVQmfPhtde83Y0hY/NVDPGZ+R6lpqIXA+8C9QAlgEJwO+qOspt0eUu\nrrLA66r6oIj0BI6o6kIRmaGq92SwvXfu4aS3e7fzy/HZZ+GRR7wdTeFx+rRz/+avvyAgwNvRGFNo\neHQtNVX9n6o2Bm4GvgVKA8+IyAER+VREonJzXBEZLyL7RWRtuvbWIrJBRDaJyMCLHGIwkFKIphqw\nM/nf53MTj8dUrQpffw3DhsHUqd6OpvAoXRoaNnTqExljvCrPU6NUdYOqjlXVbqpaFWiCM+R2dS4P\n+THQKnWDiPjhJJFWQF2gm4jUSf4sSkRGiUgVERkBfKGqa5J33YmTdABylIm94tJLYfFiGDDAKW1g\n3MOG1YzxCcXdfUBV3QJsycP+y0UkIl1zI5yZcfEAIjID6AhsUNXJwGQReQK4BQgSkctU9UPgM2CM\niNwBzM9tTB5Vty588w3cdhucOwe9enk7ooKveXPnQVtjjFe5PeHkk6r8MzQGzuSERqk3UNV3gHfS\ntZ0C7svq4KmXaYiMjCQyMjL3kbrDVVc5lUJvvdVJOvdl+RXMxTRr5hTCO3cOiheUH3ljfEtcXFye\nlwDzyaVtkns481W1fvL7O4FWqvpQ8vueQCNV/Y8bzuUbkwYysmkT3HKLM2X6oYe8HU3BVr++s47d\nDTd4OxJjCoXCXIBtN1A91ftqyW1u4bHFO3Oqdm34/nsYPhxef93b0RRsVpDNGLcodIt3ikgNnB7O\n1cnviwEbce7R7AVWAd1Udb0bzuW7PZwUu3ZBq1bQti28+qotg5Mb06c7tYjmzPF2JMYUCoWihyMi\n04AVQG0R2SEifVT1PPAE8BWwDpjhjmSTwmd7OCmqVXNmWf3wg3M/5+xZb0dU8KSsOODrf1wY4+MK\nXQ/HkwpEDyfFqVNw991OrZdPPgF/f29HVLDUqOFMO7/iCm9HYkyBVyh6ON7g8z2cFP7+MHeu8+T8\nTTfBnj3ejqhgsedxjMkz6+HkQYHq4aRQhREjnEJu8+Y5lURN1j780BmWnDjR25EYU+BZD6eoEIFB\ng+DNN53JBHPnejuigsFmqhnjVfYUHM6Qmk888JlTd93l3Jfo1AnWr4fnnnOSkcnYlVfCsWPOUGSV\nKt6OxpgCKS8PgNqQWkEcUktv1y4n+VSu7BR1Cw72dkS+q2NHZ9WBrl29HYkxBVpuhtQs4RSGhAOQ\nmAhPP+2sOD1nDtSr5+2IfNNbbzkrN5Qp4+1IjPGccePgjjvcekhLOLkgIhodHV0wh9QyMnmyk3hG\nj4Zu3bwdje9JSoL9+70dhTGeVbasU6rDDVKG1GJjYy3h5FSh6eGktmaNM8TWsqUzscAKjxlj3Mxm\nqRlHgwbw88/OMNt118Evv3g7ImOMsYQDBejBz5wICnKeNxkyxJk6/eabznCSMcbkgT34mQeFckgt\nvW3boEcPCAx0bh6Gh3s7ImNMAWdDaiZjNWvC0qXOcjjXXuskncKeZI0xPsd6OEWhh5Pa779D794Q\nEuIknurVs9zFGGPSsx6OyVq9erBypTOD7brr4P337d6OMcYjLOFQSCcNXEzx4s5abHFxMGkSNG1q\nM9mMMdlikwbyoMgNqaWXlAQffwzPP+88KPrSS84MN2OMuQgbUjM55+cH998P69bBiRPOApczZtik\nAmOM21kPp6j3cNL74Qfo29dZa2zUKGjUyNsRGWN8kPVwTN41a+asUnDffU7Zg549YedOb0dljCkE\nLOGYCxUr5iScjRudZ3gaNnRWWD5+3NuRGWMKMEs4FMFZatkVGAhDh8Kvv8L27XD55c4w2+nT3o7M\nGOMlNkstD+weTg789puzNtvq1TB4sNMLKlnS21EZY7zA7uGY/HX11fDZZ/+86tRxFgg9d87bkRlj\nCgDr4VgPJ/eWLnXu7ezaBQMHwr33QqlS3o7KGOMB1sMxnnXTTbBkCUyY4PR4Lr3UKYNw8qS3IzPG\n+CBLOCbvWrSAL7+Ezz+HFSucmW3DhsHhw96OzBjjQyzhGPe59lqYNcsZatu61enxPPoobNjg7ciM\nMT6g0CQcEekoIh+KyHQRuS25raaIjBORT7wdX5FSp46zPtv69VC+PNx8M7RtC19/bUvmGFOEFbpJ\nAyJSFnhdVR9M1faJqnbJZHubNJDfTp+GadOc+zsATz7pLBRapox34zLG5FqhmDQgIuNFZL+IrE3X\n3lpENojIJhEZeJFDDAbezd8oTY6ULu0sEPrbb/DWW7BggVP4rW9fWLs26/2NMYWCzyUc4GOgVeoG\nEfEDxiS31wW6iUid5M+iRGSUiFQRkRHAF6r6a7pj5igLm3wiArfeCvPmOYmmfHlnqK1pU6cuj61g\nYEyh5nMJR1WXA0fSNTcCNqtqvKqeBWYAHZO3n6yqTwN3ArcAd4nIQwAiEiIiY4GGWfSKjKdVqwYx\nMc6SOQMHOiURwsOd4bZffrF7PcYUQsW9HUA2VQVSL1m8CycJuajqO8A76doOA49mdfDU6wJFRkYS\nGRmZ+0hNzhQvDh07Oq/t253JBv/+t7OO2733Qo8eUKmSt6M0psiLi4vL85qTBSXh5DtLND6gRg2I\njYXoaGdq9cSJTkG4pk2d5NOhA1xyibejNKZISvkdmZfE45Oz1EQkApivqvWT398IxKhq6+T3zwGq\nqq+64Vw2S82XnTwJc+Y4yeeXX6BzZ+jaFVq2dHpHxhivKBSz1JIJaW/0rwYuE5EIESkJ3AN87q6T\nWXkCHxYQAFFR8M03TpmEq66CF16AqlXhscecnlBSkrejNKbIKFTlCURkGhAJhAL7gWhV/VhE2gBv\n4STJ8ao6wk3nsx5OQfTnnzBzpvM6cAC6dHF6Po0agZ+v/h1lTOGRmx6OzyUcTxMRjY6Otns4BdmG\nDf8kn4QEZwJCp07OCgdWr8cYt0q5hxMbG2sJJ6esh1PIbNwIc+c6r40boU0bJ/m0bu3MfDPGuIX1\ncHLBEk4htnevs4L13Lnwww9OOYUOHZwkFB7u7eiMKdAs4eSCDakVEceOOSUUFiyARYugShVnlYM2\nbZxp1yVKeDtCYwoEG1LLA+vhFA4xcTHELom9oD365mhiImPSNp4/T8z0h4ndOv7C7a99mpj2I/N2\nfNu+0G9vCte0aGPyT7FiztI6GXnvPbjmGnjpJc/GZEwRYE/O4TyHY0NqBoABA6Dk7c60a2PMBfKy\n0oAlHMj1Q0ymEPLzg+bNnZcx5gIpf5zHxl44BJkVG1IzxhjjETZpwGapGWNMttkstTywWWrGGJNz\nNkvNGGOMz7KEY4wxxiMs4RhjjPEISzhYPRxjjMmuQlUPx9Ns0oAxxuScTRowxhjjsyzhGGOM8QhL\nOMYYYzzCEo4xxhiPsIRjjDHGIyzhYNOijTEmu2xadB7YtGhjjMk5mxZtjDHGZ1nCMcYY4xGWcIwx\nxniEJRxjjDEeUdzbAbiTiHQE7gACgfGq+k26tv+q6tfejNEYY4qqQtXDUdV5qvoQ8CjQNYO2Lt6M\nz9fZ1PB/2LX4h12Lf9i1yBufTDgiMl5E9ovI2nTtrUVkg4hsEpGBFznEYODdbLSZVOx/pn/YtfiH\nXYt/2LXIG59MOMDHQKvUDSLiB4xJbq8LdBOROsmfRYnIKBGpIiIjgC9U9ddU+17QZowxxrN88h6O\nqi4XkYh0zY2AzaoaDyAiM4COwAZVnQxMFpEngFuAIBG5TFU/zKjNg1/FGGNMMp9daSA54cxX1frJ\n7+8EWiXfj0FEegKNVPU/eTyPb14AY4zxcTldacAnezielNMLZowxJnd89R5ORnYD1VO9r5bcZowx\npgDw5YQjya8Uq4HLRCRCREoC9wCfeyUyY4wxOeaTCUdEpgErgNoiskNE+qjqeeAJ4CtgHTBDVdd7\nM05jjDHZ57OTBjxBRFoDb+Ek3vGq+qqXQ/IYERkPtAP2p5qYUQ6YCUQA24EuqnrMa0F6iIhUAyYB\nFYEk4CNVHV0Ur4eIlAKWAiVx7vF+qqqxRfFagOtxjP8Bu1S1Q1G9DgAish04hvP/yFlVbZTT6+GT\nPRxPuNhzPUXEBc86Ac8B36jqFcB3wCCPR+Ud54CnVbUu0ATom/yzUOSuh6omAi1V9RqgIdBGRBpR\nBK9FsieBP1K9L6rXAZxEE6mq16hqo+S2HF2PIptwSPVcj6qeBVKe6ykSVHU5cCRdc0dgYvK/JwKd\nPBqUl6jqvpSHglX1BLAeZ1JKUb0ep5L/WQqnl6MUwWuR3PNtC4xL1VzkrkMqwoU5I0fXoygnnKrA\nzlTvdyW3FWUVVHU/OL+EgQpejsfjRKQGzl/2K4GKRfF6iIifiPwC7AO+VtXVFM1r8SYwACfhpiiK\n1yGFAl+LyGoReSC5LUfXo8g/h2Muqkjd4BORMsCnwJOqeiKDh4KLxPVQ1STgGhEJAj4Tkbpc+N0L\n9Yo/PJ0AAAKvSURBVLUQkTtw7m/+KiKRF9m0UF+HdJqp6l4RKQ98JSIbyeHPRVHu4dhzPRfaLyIV\nAUSkEvCXl+PxGBEpjpNsJqvqvOTmIns9AFQ1AYgDWlP0rkUzoIOI/AlMB/4lIpOBfUXsOrio6t7k\n/x4A5uLclsjRz0VRTjj2XM+Fzzp9DvRO/ve9wLz0OxRi/wX+UNW3U7UVueshImEiEpz879LAbTj3\ntIrUtVDV51W1uqrWwvnd8J2qRgHzKULXIYWI+CePACAiAcDtwG/k8OfCpkXD2/wzLXqEl0PymORn\nnSKBUGA/EI3zV8ssIByIx5nieNRbMXqKiDTDmQr8G86QgALPA6uATyhC10NErsa5+euX/Jqpqi+L\nSAhF7FqkEJGbgWeSp0UXyesgIjWBz3D+3ygOTFXVETm9HkU64RhjjPGcojykZowxxoMs4RhjjPEI\nSzjGGGM8whKOMcYYj7CEY4wxxiMs4RhjjPEISzjGGGM8whKOMcYYj7CEY4yPE5E+3o7BGHewlQaM\n8WEi0g3wx6lNs1JV/8/LIRmTa1aewBgfJSLvAL+p6ofejsUYd7AejjE+KLkGzc9AgKqe93Y8xriD\n3cMxxjddDiRYsjGFifVwjPFBIlIB+B2YA6wF9qnqHO9GZUzeWA/HGN9UDPgCOIZzr/WMd8MxJu+s\nh2OMjxGRUsBU4L7kMs/GFArWwzHG9zwCfGbJxhQ2Ni3aGN9zGc6QmjGFig2pGeNjROQa4EtgKfCB\nqn7r5ZCMcQsbUjPGx6jqL0BdYD3wXxHZIiKR3o3KmLyzHo4xPk5EHgRaq+qd3o7FmLywezjG+CgR\naQA0B6oDb3s5HGPyzHo4xhhjPMLu4RhjjPEISzjGGGM8whKOMcYYj7CEY4wxxiMs4RhjjPEISzjG\nGGM8whKOMcYYj/h/S8MyJq1/2/cAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "xi = np.arange(1e-6, 20, 0.01) #start at small poisitive number to avoid div-by-0\n", "plt.semilogy(xi,rho_c*D,label ='Original L.E. solution',color=\"red\")#$\\rho$')\n", "plt.semilogy(xi[:int(xi_max*100)],fattening_factor*rho_c*D[:int(xi_max*100)],linewidth=5,linestyle=\"--\",label ='Fattened L.E. solution')\n", "xi_out = np.arange(xi_max,40)\n", "plt.semilogy(xi_out,rho_out*ones_like(xi_out),linewidth=5,linestyle=\"--\",label ='Density drop by 150')\n", "sim_xi = rp.x.value*1.05e-17 #xi values from the simulation\n", "sim_dens = rp[\"density\"].in_units(\"g/cm**3\").value #density values from the simulation\n", "plt.semilogy(sim_xi,sim_dens,label=\"Simulation\")\n", "\n", "plt.xlabel(r\"$\\xi$\",fontsize=15)\n", "plt.ylabel(r\"$\\rho \\quad [g/cm^3]$\",fontsize=15)\n", "plt.legend(loc='upper right',prop={'size':12},numpoints=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Since the Lane-Emden solution is self simmilar, multiplying the inside density by a constant factor would not matter. The original outside density is already very low ($\\approx 10^{-21} g/cm^3$ ), so multiplying it by a ``fattening_factor`` should not affect the sphere's behaviour by much." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/global/u2/d/dorislee/astroSim-tutorial\n" ] } ], "source": [ "cd .." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/anaconda2/lib/python2.7/site-packages/matplotlib/__init__.py:1350: UserWarning: This call to matplotlib.use() has no effect\n", "because the backend has already been chosen;\n", "matplotlib.use() must be called *before* pylab, matplotlib.pyplot,\n", "or matplotlib.backends is imported for the first time.\n", "\n", " warnings.warn(_use_error_msg)\n" ] } ], "source": [ "from scripts.plotSim import *" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/global/project/projectdirs/astro250/dlee/FLASH4.3/object\n" ] } ], "source": [ "cd ~/proj/dlee/FLASH4.3/object/" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_dens(0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The temperature outside is only a factor of ``fattening_factor`` more than the inside of the sphere, which is less extreme than the $10^7$K proposed by Foster and Chevalier: " ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_var(0,\"temperature\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For a visual check, in the pressure plot the outer region of the dense core should be the same color as the ambient, since $P_{out}=P_{edge}$." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_var(0,\"pressure\")" ] }, { "cell_type": "code", "execution_count": 90, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plot_dens(65,\"sphere\",grid=True,velocity=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Free fall time \n", "\n", "To estimate how long your simulation is going to take, use the free fall time as an order-of-magnitude estimate: \n", "\n", "$$t_{ff} = \\sqrt{\\frac{3\\pi}{32G\\rho}}$$\n", "\n", "We could check that this agrees with the ``ds.current_time`` of the last timestep in the simulation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reference: \n", "\n", "[1] P. N. Foster, R. A. Chevalier. Gravitational Collapse of an Isothermal Sphere. ApJ 416, 303-311 (1993).\n", "\n", "[2] R. B. Larson. Numerical Calculations of the Dynamics of a Collapsing Proto-Star. MNRAS 145, 271–295 (1969). [Link](http://mnras.oxfordjournals.org/content/145/3/271)" ] } ], "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 }