Hydrostatic Equilibrium (rotating cylindrical polytrope).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Hydrostatic Equilibrium of Rotating Cylindrical Polytropes\n",
    "\n",
    "\n",
    "## Lane-Emden Equation\n",
    "\n",
    "The Lane-Emden equation is $$\\theta''+\\frac{\\theta'}{\\xi}=-\\theta^n+\\beta$$\n",
    "\n",
    "where β is a constant being correlate with the rotation of the star,\n",
    "\n",
    "n is the polytropic index from the polytropic EOS ($P=K\\rho^\\Gamma$ with $\\Gamma=1+\\frac{1}{n}$),\n",
    "\n",
    "and θ from the density's function $\\rho(\\theta)={\\rho_c}\\theta^n$\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### *Numerical Solution for the Lane-Emden Equation*\n",
    "\n",
    "#### Runge-Kutta Method:\n",
    "\n",
    "- The DE is $$\\theta''+\\frac{\\theta'}{\\xi}=-\\theta^n+\\beta$$\n",
    "\n",
    "\n",
    "- Assuming $\\theta'=\\frac{d\\theta}{d\\xi}=z$, $\\theta''=\\frac{d^2\\theta}{d\\xi^2}=z'$ and also replacing ξ with x, θ with y and β with b\n",
    "\n",
    "\n",
    "- The DEs are $$z'=-\\frac{z}{x}-y^n+b=f(y,z)$$ $$and$$ $$y'=z$$\n",
    "\n",
    "\n",
    "- Initial Conditions: The density at the center of the star is known to be equal to $\\rho_n$, so $\\theta(0)=y(0)=1$ and $z(0)=0$. Therefore, $z'(0)=-\\frac{z(=0)}{x(=0)}-1+b$ $\\implies(De L' Hospital)\\implies f(0)=z'(0)=\\frac{-1+b}{2} $\n",
    "\n",
    "\n",
    "- $\\underline{Numerical}$ $\\underline{Solution}$: $${y_{i+1}}={y_i}+\\frac{k_1 +2k_2 +2k_3+ k_4}{6}$$ \n",
    "$${z_{i+1}}={z_i}+\\frac{l_1 +2l_2 +2l_3+ l_4}{6}$$\n",
    "\n",
    "\n",
    "where \n",
    "\n",
    "$k_1=hy'_i,$\n",
    "\n",
    "$k_2=hy'(z_i+\\frac{l_1}{2}),$\n",
    "\n",
    "$k_3=hy'(z_i+\\frac{l_2}{2}),$\n",
    "\n",
    "$k_4=hy'(z_i+l_3).$\n",
    "\n",
    "$l_1=hf_i,$\n",
    "\n",
    "$l_2=hf(x_i+\\frac{h}{2},y_i+\\frac{k_1}{2},z_i+\\frac{l_1}{2}),$\n",
    "\n",
    "$l_3=hf(x_i+\\frac{h}{2},y_i+\\frac{k_2}{2},z_i+\\frac{l_2}{2}),$\n",
    "\n",
    "$l_4=hf(x_i+h,y_i+k_3,z_i+l_3).$\n",
    "\n",
    "\n",
    "## Rotation\n",
    "\n",
    "The dimensionless variable b is $$\\beta=\\frac{\\Omega^2}{2\\pi G\\rho_c}$$\n",
    "\n",
    "Every star has a maximum Rotation (Ω) depending from the polytropic index n. If b exceed the maximum limit $\\beta_c$, the star will be destroyed (mass shedding). In this case, the equilibrium not apply and therefore it is $\\frac{dP}{dr}=0$ at the surface of the star. So, at this boundary condition the Lane-Emden Equation (for $\\xi=\\xi_1$) takes the form: $$\\theta''(\\xi_1)=\\beta$$.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***************************************"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "import numpy as np\n",
    "from scipy import fftpack \n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import scipy as sp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set step h and DE:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "h=1e-5\n",
    "de= lambda x,y,z : -(z/x)-pow(y,n)+b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set constant b (*rotation*) and polytropic index:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "b=0.2738\n",
    "n=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set initial values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "x=[0]\n",
    "y=[1]\n",
    "z=[0]\n",
    "f=[(b-1)/2]\n",
    "\n",
    "k1=[]\n",
    "k2=[]\n",
    "k3=[]\n",
    "k4=[]\n",
    "l1=[]\n",
    "l2=[]\n",
    "l3=[]\n",
    "l4=[]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RK(i) :\n",
    "    x.append(x[i]+h)\n",
    "    \n",
    "    k1.append(h*z[i])\n",
    "    l1.append(h*f[i])\n",
    "    \n",
    "    k2.append(h*(z[i]+0.5*l1[i]))\n",
    "    l2.append(h*de(x[i]+h*0.5,y[i]+k1[i]*0.5,z[i]+l1[i]*0.5))\n",
    "    \n",
    "    k3.append(h*(z[i]+0.5*l2[i]))\n",
    "    l3.append(h*de(x[i]+h*0.5,y[i]+k2[i]*0.5,z[i]+l2[i]*0.5))\n",
    "    \n",
    "    k4.append(h*(z[i]+l3[i]))\n",
    "    l4.append(h*de(x[i]+h,y[i]+k3[i],z[i]+l3[i]))\n",
    "    \n",
    "   \n",
    "    y.append(y[i]+(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6)\n",
    "    z.append(z[i]+(l1[i]+2*l2[i]+2*l3[i]+l4[i])/6)\n",
    "    f.append(de(x[i+1],y[i+1],z[i+1]))\n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "i=0\n",
    "\n",
    "while np.real(y[i])>=0:\n",
    "      RK(i)\n",
    "      i +=1\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The root is ξ1=  3.478180000014319\n"
     ]
    }
   ],
   "source": [
    "print('The root is ξ1= ',x[i-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot: x-y (Radious-Density)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAF8CAYAAAAn/HmMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcjWXjx/HPNQtjK7JFZBdClskuk50RHlH8Ii3Sk7Q8laJ6iJSkIi2iIqkIZUkUQqQsUyEpssVEPfYtu+v3xzUmZjPDzLnP8n2/Xudlzjm3M9/bGec793ZdxlqLiIiIBKcwrwOIiIhI1lHRi4iIBDEVvYiISBBT0YuIiAQxFb2IiEgQU9GLiIgEMRW9iIhIEFPRi4iIBDEVvYiISBBT0YuIiASxCK8DZIYCBQrYkiVLeh1DRETEJ77//vvd1tqC6Vk2KIq+ZMmSxMXFeR1DRETEJ4wxv6d3We26FxERCWIqehERkSCmohcREQliQXGMXkRSd/LkSeLj4zl27JjXUUQkg6KioihWrBiRkZEX/RoqepEgFx8fT548eShZsiTGGK/jiEg6WWvZs2cP8fHxlCpV6qJfR7vuRYLcsWPHyJ8/v0peJMAYY8ifP/8l741T0YuEAJW8SGDKjP+7KnoRyXKTJ0+mYsWKVKxYkSlTpiQ+fvToURo1asTp06dZuHAhRYsWxRhD8eLFGTp0KACvv/4648aNu6jvO3PmTF544QUApk+fzrp16xKfi4mJueD4G6llWrx4MTVq1CAiIoKpU6dmKNPdd9/NddddR9WqVenYsSOHDx9Oddlt27aRO3duXnrpJQDWr19PtWrVEm+XXXYZI0aMAGDVqlXUqVOHatWqER0dzYoVKwD48MMPqVq1KlWrVqVevXqsXr06xe9lraVx48YcPHgwQ+uT1J49e7jxxhvJnTs3vXv3TnW5vXv30qxZM8qVK0ezZs3Yt2/fJX3flCxatIg2bdoke/y9995LM1ta0pN71apV1K1bl2uvvZaqVavy8ccfJz7XsGHDxPevaNGitG/fHoBZs2YxYMCAi8p0QdbagL/VrFnTikjK1q1b5+n3P3jwoK1Ro4Y9dOiQPXjwoK1evbo9dOiQtdba119/3Y4YMcJaa+2aNWtsvXr17Pbt28/7+0eOHLHVqlW75Bzdu3e3U6ZMSbzfqFEju3LlyjT/TmqZtmzZYlevXm27det23mumx4EDBxK//s9//mOHDBmS6rIdOnSwHTt2tMOGDUv23KlTp2zhwoXt1q1brbXWNmvWzM6ePdtaa+3nn39uGzVqZK21dunSpXbv3r3WWmtnz55ta9WqleL3mjVrln344YcztC4pOXz4sF2yZIkdNWqUvf/++1Ndrk+fPonrPmTIEPv4449f8vdOauHChTY2NjbZ4+PGjUszW1rSk3v9+vV2w4YN1lpr//jjD3vllVfaffv2JVuuQ4cOdvz48dZaa8+cOWOrVatmjxw5kmy5lP4PA3E2nR3p0y16Y8xYY8z/jDFrU3neGGNGGmM2GmPWGGNq+DKfiGS+L774glatWpE7d27y5MlD06ZNmTNnDuC2Ntu1awfA9u3byZs3L7lz5z7v7+fMmZOSJUsmbqGedfr0aUqXLo21lv379xMWFsbixYsBt9W0cePGxC23b7/9lpkzZ9KnTx+qVavGpk2bAJgyZQq1atWifPnyLFmyJFn21DKVLFmSqlWrEhaW8Y/Qyy67DHAbWUePHk111+z06dMpXbo01157bYrPf/XVV5QpU4YSJUoAbhfv2a3xAwcOULRoUQDq1atHvnz5AKhTpw7x8fEpvt6578XWrVupWLEi99xzD9deey3Nmzfn6NGj6Vq/XLly0aBBA6KiotJcbsaMGXTv3h2A7t27M3369GTLbN26lYYNG1KjRg1q1KjBt99+C7gt9ZiYGDp27EiFChW47bbbcN3nft4qVKhAgwYN+PTTT1P9/tu3b6dly5Zcc801DBw4MF3rlt7c5cuXp1y5cgAULVqUQoUKsWvXrvOWOXToEAsWLEjcojfGEBMTw6xZs9KdJb18fdb9e8DrwPupPN8KKJdwqw2MSvhTRDLDww/DqlWZ+5rVqkHC7uOUbNu2jXHjxiV+gO3atYvChQtz4sQJNm/ezNl5Klq1asXMmTO54oorKFq0KGvXriVv3rwAREdHs2TJEmrVqpX4uuHh4ZQvX55169axZcsWatasyZIlS6hduzbx8fGULVuWb775BnBl17ZtW9q0aUPHjh0TX+PUqVOsWLGC2bNnM3DgQObPn39e9rQypeTQoUM0bNgwxec++ugjKlWqBMCdd97J7NmzqVSpEi+//HKyZY8cOcLQoUOZN29e4m77pCZNmkSXLl0S748YMYIWLVrw2GOPcebMmcRSPNe7775Lq1atUny9pUuXMnr06MT7v/32GxMnTuTtt9/mlltu4ZNPPqFr164MGzaMDz/8MNnfv+GGGxg5cmSKr52Sv/76iyJFigBQpEgR/ve//yVbplChQsybN4+oqCh+++03unTpkni45ccff+Tnn3+maNGi1K9fn6VLlxIdHc0999zDggULKFu2LLfeemuq33/FihWsXbuWnDlzcv311xMbG0t0dDQNGzbk0KFDyZZ/6aWXaNq0abpyJ/0+J06coEyZMuc9Pm3aNJo0aZL4ix/883N+yy23pPmaGeXTorfWLjbGlExjkXbA+wm7JZYZY/IaY4pYa3f6JKCIZDprLffddx9PP/00AIMGDcIYw+7du88rzdmzZ7N8+XLWrl2bWIhnFSpUiF9//TXZazds2JDFixezZcsW+vXrx9tvv02jRo24/vrr05WtQ4cOANSsWZOtW7cmez6tTCnJkycPq9Lxi9S4ceM4ffo0DzzwAB9//DF33nnnec8PGDCA//znP8n2JJx14sQJZs6cyZAhQxIfGzVqFMOHD+fmm29m8uTJ3H333ef94rJw4ULefffdxF9+ktq7dy958uRJvF+qVCmqVasGnP/v06dPH/r06XPBdcwMJ0+epHfv3qxatYrw8HA2bNiQ+FytWrUoVqwYANWqVWPr1q3kzp2bUqVKJW5Nd+3alTFjxqT42s2aNSN//vyA+zn45ptvEos2s+zcuZNu3boxfvz4ZHt/Jk6cSI8ePc57rFChQuzYsSPTvv9Z/nYd/VXA9nPuxyc8lnbRb9sGjzwCkZHuli2bu11+ubvlzfvPrWBByJ8fLmKXm0jAS2PLO6uUKFGCBQsWJN7fvn07LVq0IEeOHOddNjR79mx69+6dYqEeO3aMHDlyJHu8YcOGvPXWW+zYsYNBgwYxbNgwFi1axA033JCubNmzZwfc3oFTp04lez6tTClJ7xb92e956623MmzYsGRFv3z5cqZOncrjjz+eeFgiKioq8QSyOXPmUKNGDQoXLpz4d8aPH8+rr74KQKdOnc4rkTVr1tCjRw/mzJmTWG5JRUREcObMmcRCOvtvczbr2V33mbVFX7hwYXbu3EmRIkXYuXMnhQoVSrbM8OHDKVy4MKtXr+bMmTPnHQ5Imu/s+5fes9STLnf2/oW26NOTG+DgwYPExsYyePBg6tSpc95ze/bsYcWKFUybNu28x1P7Ob9U/lb0Kb1DNsUFjekJ9ASoERYGY8bAyZPuZlP8K/+IjIQiReCqq6BoUShWDMqUgXLl3K1ECYjwt38akcDUqlUr/vvf/7Jnzx6stSxatIjhw4eTO3duTp8+zbFjx4iKiuKaa67h66+/pnv37kQk+f+3YcMG6tevn+y1a9euze23307p0qWJioqiWrVqjB49OsXjnHny5EnxAzwtaWVKyYW26K21bNq0ibJly2Kt5bPPPqNChQrJljt3q/KZZ55Jdgb7xIkTz9ttD+5Y8Ndff01MTAwLFixI3Krdtm0bHTp0YMKECZQvXz7Ndd28eTNly5ZNcx0za4u+bdu2jB8/nr59+zJ+/PjE8wPOdeDAAYoVK0ZYWBjjx4/n9OnTab5mhQoV2LJlC5s2baJMmTJMnDgx1WXnzZvH3r17yZEjB9OnT2fs2LEAF9yiT0/uEydO8K9//Yvbb7+dTp06JXt+ypQptGnTJtl5DBs2bKBy5cppfv+L4W+btfFA8XPuFwNS3I9hrR1jrY221kab6tXh8GE4fhzOnIFTp+DIEdi5E375BZYtgy++gIkTYeRIePRRaNQIcuaEdevg3XfhwQehVSsoWxZy5IDy5aF9exgwAD79FDZtcq8tIhmSO3duBg8eTN26dalbty5DhgxJ3CXdvHnzxF3J9913H/ny5aN69erUrl2byZMnJ77G0qVLadq0abLXzp49O8WLF0/cYjq7NValSpVky3bu3Jlhw4ZRvXr1xJPxLiS1TCtXrqRYsWJMmTKFe++9N9UT5pKy1tK9e3eqVKlClSpV2LlzJ/379wfcpYBnv07L33//zbx58xIPO5z19ttv8+ijj3Ldddfx5JNPJu6yHjRoEHv27KFXr16Jl96lJDY2lkWLFqVrPS6kZMmSPPLII7z33nsUK1Ys8bLGHj16JB5j79u3L/PmzaNcuXLMmzePvn37JnudXr16MX78eOrUqcOGDRvIlStXmt83KiqKMWPGEBsbS4MGDRJPVExJgwYN6NatG9WqVePmm29O9d8lqdRyx8XFJe5FmTx5MosXL+a9995LvJTu3F8Ak55fcdbChQuJjY1NV44MSe/p+Zl1A0oCa1N5LhaYg9uyrwOsSM9rXvLldWfOWLtzp7VLllg7dqy1/fpZ27GjtRUqWBsWZq3bR2BtnjzWNmhg7WOPWfvpp+7viPg5ry+vS8sPP/xgu3btmuzx+Ph4W7ly5TSX8bVzMwWjHTt22KZNm3odI2T9+eeftnHjxik+d6mX1/l0/7QxZiIQAxQwxsQDA4DIhF843gJmA62BjcDfwJ0pv1KmB4Mrr3S3Bg3Of+7vv2HtWli92t2+/97tFTh7JmypUlCvHtxwAzRtCqVL+ySySDCoXr06N954I6dPnyY8PDzx8ZdffpnbbrsNgN27d/Pss896FTHRuZmCUZEiRbjnnns4ePDgeWeCi29s27YtxSswMoOxFzqeHQCio6PthUa4ylTHj8MPP8C338J338HSpfDnn+650qVd4TdtCo0buxP/RDz0yy+/ULFiRa9jiMhFSun/sDHme2ttuo436Iyzi5E9O9St627gduyvXw/z57vbpEnu5EBj3DLt2kHbtpDCSTciIiJZyd9OxgtMxrgS790bpk+HPXvcln7//nDsGDzxBFSsCNdcA336uD0BOrFPfCgY9tyJhKLM+L+ros8KERFQpw4884w7pr9tG7zxBpQsCa++CvXru2P7TzzhRinTh7BkoaioqMRL20QkcFjr5qO/0HDCF6Jj9L528CDMnOku9Zs7110KWKECdOkCXbvqZD7JdCdPniQ+Pv6S57QWEd+LioqiWLFiREZGnvd4Ro7Rq+i9tHs3fPKJK/3Fi92WfePGcPfd0KEDXOJvcSIiEpwyUvTade+lAgXg3nth0SL4/Xd49lnYsgVuu82N2PfAA5k/AYmIiIQUFb2/KF4cnn4aNm50Z+63bAlvvw3Vq7vr9CdNcsP7ioiIZICK3t+EhUGTJvDRR7Bjh5uEZNcudwy/ZEkYPBguMC2iiIjIWSp6f3bFFfDQQ+4a/c8/hypV4L//dVv/d9wBP/3kdUIREfFzKvpAEBYGrVu7iXl++QXuuQemToWqVaFNGzcyn4iISApU9IGmQgV4/XV3bf6gQW5mvgYNoGFDt9UfBFdRiIhI5lHRB6orrnC78X//3U2ys22b27q/7jq3ta+R90REBBV94MuVy12Gt3EjvP++OzO/UyeoUQNmzNAWvohIiFPRB4vISOjWzU2pO2ECHDkC7dtDrVowZ44KX0QkRKnog014uBtK95df4N133eh7rVu78fW//trrdCIi4mMq+mAVEQF33eUuzXvrLXcMPybGTZf7yy9epxMRER9R0Qe7bNncMLu//QbPP++G261SBe67D/76y+t0IiKSxVT0oSJHDujXDzZtciX/zjtQtqwbX//IEa/TiYhIFlHRh5qCBeG11+Dnn6F5c+jf312bP3myTtgTEQlCKvpQVb68myJ3yRI3i96tt7ox9teu9TqZiIhkIhV9qGvQAOLiYNQoWL0aqlWDhx+G/fu9TiYiIplARS/ukrx//xs2bHDj6I8c6bb4x43T7nwRkQCnopd/5M/vtuzj4qBcOXd5XkwM/Pqr18lEROQiqegluRo13LH7d96BNWvc+PkDB8Lx414nExGRDFLRS8rCwuDuu93W/M03wzPPuOP3S5Z4nUxERDJARS9pK1wYPvrIjZd/7BjccIM7jr9vn9fJREQkHVT0kj4tW7pL7x57zJ2kV7kyzJ7tdSoREbkAFb2kX65cMGwYLF8OV1wBsbHuhD1diici4rdU9JJxNWu6M/OffBLGj3db91984XUqERFJgYpeLk727PDcc7BsGVx+ObRqBT16wIEDXicTEZFzqOjl0lx/PXz/PTzxhDt2X6WKmyFPRET8gopeLl1UFLzwAnz7rfu6cWPo2xdOnPA6mYhIyFPRS+apXRt+/NHtwh86FOrW1ah6IiIeU9FL5sqVC8aMgWnT4Pff3Sh7o0drzHwREY+o6CVrtG/vhs9t0MBNmNO+Peza5XUqEZGQo6KXrFO0qLvsbvhw92fVqrBggdepRERCiopeslZYmJvffuVKyJsXmjZ14+afPu11MhGRkKCiF9+oWtWVfdeubia8Zs3gzz+9TiUiEvRU9OI7uXO7kfTGjnUD7Vx3Hcyf73UqEZGgpqIX3zIG7rzTbd3nzw/Nm8OAAdqVLyKSRVT04o1rr3Vl3707DBrkjt1rV76ISKZT0Yt3cuVyw+a+956bEa9mTfjuO69TiYgEFRW9eK97d1fwUVHQqBG8+aYG2BERySQqevEP113npr5t1gzuvx/uuAOOHvU6lYhIwFPRi//Ilw8++8xdZz9hAtSrB1u2eJ1KRCSgqejFv4SFubPwP/sMtm51x+2//NLrVCIiAUtFL/4pNtbtyi9eHFq1giFDdNxeROQiqOjFf5Up407S69wZnnwSbrtNx+1FRDJIRS/+LWdO+PBDeP55mDQJGjaEP/7wOpWISMBQ0Yv/Mwb69YPp02H9eoiOdtfdi4jIBanoJXC0bet25efI4a63/+ADrxOJiPg9Fb0ElsqVYcUKqFsXunWDJ57QOPkiImlQ0UvgKVAA5s6Ff/8bXnwR2rWDgwe9TiUi4pd8XvTGmJbGmPXGmI3GmL4pPH+1MWahMeZHY8waY0xrX2eUABAZCaNGwRtvwBdfuJP0tm/3OpWIiN/xadEbY8KBN4BWQCWgizGmUpLFngYmW2urA52BN32ZUQJMr14we7YbQa92bfjhB68TiYj4FV9v0dcCNlprN1trTwCTgHZJlrHAZQlfXw7s8GE+CUTNm8O337qt/BtucKPqiYgI4Puivwo4d/9qfMJj53oG6GqMiQdmAw/4JpoEtMqVYdkyqFgR2reHkSO9TiQi4hd8XfQmhceSjmvaBXjPWlsMaA1MMMYky2mM6WmMiTPGxO3atSsLokrAKVIEFi1yl+E99BA8+KDOyBeRkOfroo8Hip9zvxjJd83fDUwGsNZ+B0QBBZK+kLV2jLU22lobXbBgwSyKKwEnVy6YOhUefRRee81t3R8+7HUqERHP+LroVwLljDGljDHZcCfbzUyyzDagCYAxpiKu6LXJLukXHg4vvQRvvulO1LvhBti50+tUIiKe8GnRW2tPAb2BL4FfcGfX/2yMGWSMaZuw2KPAPcaY1cBE4A5rNW2ZXIT77oNZs2DDBjfAzvr1XicSEfE5EwwdGh0dbePi4ryOIf4qLs5Ne3vqlCv+unW9TiQickmMMd9ba6PTs6xGxpPgFx3tLr/Llw+aNNHldyISUlT0EhrKlHFlf+217gS9d97xOpGIiE+o6CV0FCoECxdCixZwzz0wcCAEwaErEZG0qOgltOTODTNmwB13wDPPuIlxTp3yOpWISJaJ8DqAiM9FRsLYsXDVVfDcc/DnnzBxIuTM6XUyEZFMpy16CU3GwODBbva7zz5z4+Xv3+91KhGRTKeil9DWqxdMngwrVkBMDPz1l9eJREQylYpepGNH+Pxz2LgRGjSArVu9TiQikmlU9CIAzZrB/PmwZw/Urw/r1nmdSEQkU6joRc6qUwcWL3aX3DVs6Hbni4gEOBW9yLkqV4ZvvoG8eaFxY/jqK68TiYhcEhW9SFKlS7uyL1UKWreGadO8TiQictFU9CIpKVIEvv4aatRwJ+uNHet1IhGRi6KiF0nNFVe4E/SaNoW774YRI7xOJCKSYSp6kbTkygUzZ8LNN8N//gPPP+91IhGRDFHRi1xI9uwwaRLcdhs89RQ8/bQmwxGRgKGx7kXSIyICxo+HHDnc+PhHj8JLL7mhdEVE/JiKXiS9wsNh9GhX9q+84sr+9dchTDvGRMR/qehFMiIsDF591c10N3SoK/t33nG/BIiI+CEVvUhGGQNDhriyHzDAlf2ECW76WxERP6OiF7kYxkD//m43/uOPw7Fj8PHH7sQ9ERE/ooOLIpeiTx947TWYMQPatYO///Y6kYjIeVT0Ipeqd293nH7uXIiNhSNHvE4kIpJIRS+SGe6+Gz74wM1+17o1HD7sdSIREUBFL5J5/u//4KOPYOlSaNUKDh3yOpGIiIpeJFPdeqsr++++U9mLiF9Q0YtktltucUPmLlsGLVvCwYNeJxKREKaiF8kKHTu6y+1WrIAWLeDAAa8TiUiIUtGLZJWbb4bJkyEuTmUvIp5R0YtkpX/9C6ZOhR9+gObNYf9+rxOJSIhR0YtktXbtXNn/+CM0awb79nmdSERCiIpexBfatoVPP4U1a1T2IuJTKnoRX2nTxpX9Tz9B06awd6/XiUQkBKjoRXwpNhamT4e1a3XMXkR8QkUv4mutWv2zG1/X2YtIFlPRi3ghNhamTIHvv3fFr7HxRSSLqOhFvNKunRtBb/lyzXonIllGRS/ipZtvhg8/hG++cWfmaz57EclkKnoRr916K4wfDwsXQvv2cOyY14lEJIio6EX8QdeuMHYszJ8PHTrA8eNeJxKRIKGiF/EXd9wBY8bAnDnQqROcOOF1IhEJAip6EX/Sowe8+SZ89hl07gwnT3qdSEQCnIpexN/cdx+MHAnTpsFtt8GpU14nEpEAFuF1ABFJwQMPuK35Rx+FyEh4/30ID/c6lYgEIBW9iL965BFX9n37urIfOxbCtBNORDJGRS/iz554wp2U178/5Mjhjt8b43UqEQkgKnoRf/f0027UvKFDIVcuGDZMZS8i6aaiF/F3xsCQIW7UvJdfdmU/cKDXqUQkQKjoRQKBMTBihNuyHzTIlf3jj3udSkQCgIpeJFCEhbkBdY4edcfuc+aE3r29TiUifk5FLxJIwsPduPh//+0uwcuZE+66y+tUIuLHdK2OSKCJjISPP4YWLdxIehMnep1IRPyYil4kEGXPDp9+Cg0bQrduMH2614lExE/5vOiNMS2NMeuNMRuNMX1TWeYWY8w6Y8zPxpiPfJ1RJCDkzAmzZkF0tJvq9ssvvU4kIn7Ip0VvjAkH3gBaAZWALsaYSkmWKQf0A+pba68FHvZlRpGAkiePm+2uUiU3l/3XX3udSET8jK+36GsBG621m621J4BJQLsky9wDvGGt3Qdgrf2fjzOKBJZ8+WDuXChVCtq0gWXLvE4kIn7E10V/FbD9nPvxCY+dqzxQ3hiz1BizzBjT0mfpRAJVwYIwfz4ULgytWsGqVV4nEhE/4euiT2ncTpvkfgRQDogBugDvGGPyJnshY3oaY+KMMXG7du3K9KAiAadoUfjqK7c7v1kzWLfO60Qi4gd8XfTxQPFz7hcDdqSwzAxr7Ulr7RZgPa74z2OtHWOtjbbWRhcsWDDLAosElBIlXNlHREDTprB5s9eJRMRjvi76lUA5Y0wpY0w2oDMwM8ky04EbAYwxBXC78vVpJZJe5crBvHlw/Lgr+z/+8DqRiHjIp0VvrT0F9Aa+BH4BJltrfzbGDDLGtE1Y7EtgjzFmHbAQ6GOt3ePLnCIBr3Jldzb+rl3QvDns3u11IhHxiLE26SHywBMdHW3j4uK8jiHifxYtgpYtXfEvWACXXeZ1IhHJBMaY76210elZViPjiQSzmBiYOhVWr4abbnIT4ohISFHRiwS7Nm3g/fdhyRLo2BFOnPA6kYj4kIpeJBR06QKjRsHs2XD77XD6tNeJRMRHNE2tSKi49144cMDNZX/ZZTB6NJiUhrYQkWCiohcJJY8/7sr++efh8svhxRdV9iJBTkUvEmoGD3Zl/9JLkDcvPPWU14lEJAup6EVCjTEwcqQr+6efdrvxH3jA61QikkVU9CKhKCwMxo2DQ4fgwQfdbvzbb/c6lYhkAZ11LxKqIiJg0iRo0gTuvBOmTfM6kYhkARW9SCiLioLp06FWLejc2Y2RLyJBRUUvEupy53bX119zDbRvD99953UiEclEKnoRgXz5YO5cN6d969ZuyFwRCQoqehFxrrwS5s93W/jNm8OGDV4nEpFMoKIXkX+UKOHK3lo3l/22bV4nEpFLpKIXkfNdcw18+SUcPOjK/q+/vE4kIpcg3UVvjClqjGlmjOlsjLk14esiWRlORDxSvTp8/jnEx0OLFrB/v9eJROQiXbDojTHNjTErge3Al8CHwMSEr+ONMcuNMU2yNqaI+Fz9+u7a+nXrIDYWjhzxOpGIXIQ0i94Y0xn4HFgONAWuBgwQk/B1C2AtMMcY0ylLk4qI77VoAR99BMuWQYcOcPy414lEJIMutEX/DDDAWtvbWrsQ2Jnw+GFrbby1dr619m5gMPBsFuYUEa907Ahvv+0uv7vtNjh1yutEIpIBFyr6q4Gl6XidxQnLikgwuusuGD4cPvkEevaEM2e8TiQi6XShSW1+AG4Hvr7Acj2B7zMlkYj4p4cfdiflDRzoJsF55RXNZS8SAC5U9A8Ac40xPwJTgd8ACzQ0xpQHygIdgaK4Y/giEswGDIB9+2DECDeaXv/+XicSkQtIs+ittT8aY6oCjwF3AKVxJ+MNxxX+JmA68Iq19s+sjSoinjMmCXuYAAAXvklEQVTG7cI/cMCV/uWXw0MPeZ1KRNJwwfnorbU7gUeBR40xUUC+hKf2Wmt1Cq5IqAkLg3fecQPqPPww5M0L3bt7nUpEUpGhkfGstcestTsTbip5kVAVEQETJ0KzZu5EvU8/9TqRiKRCQ+CKyMXJnt0NqFO7NnTpornsRfyUil5ELl6uXG6o3AoVNJe9iJ9S0YvIpUk6l/2aNV4nEpFzqOhF5NIVLnz+XPa//eZ1IhFJoKIXkcxRooQ7Tn/6tJvedvt2rxOJCCp6EclMFSq4uez373dn5O/a5XUikZCnoheRzFWjBsyaBdu2udnvDhzwOpFISFPRi0jma9jQXVu/di20aQN//+11IpGQpaIXkazRsiV88AEsXeqmuj1xwutEIiFJRS8iWeeWW2DMGJgzB7p1cyfqiYhPXXCsexGRS9KjhztO/9hjcNllrvg1va2Iz6joRSTrPfqoOxN/8GA3492wYSp7ER9R0YuIbwwa5Mr+5ZfdaHpPPeV1IpGQoKIXEd8wBl591e3Gf/ppN73t/fd7nUok6KnoRcR3wsJg7Fg3l33v3m43fteuXqcSCWo6615EfCsiAiZNgsaN4Y47YMYMrxOJBDUVvYj4XlQUTJ8O0dHuEryvvvI6kUjQUtGLiDfy5IHZs6F8eWjXDpYv9zqRSFBS0YuId664ws1lf+WV0KoV/PST14lEgo6KXkS8VaSIm8s+Rw43l/2mTV4nEgkqKnoR8V7Jkm4u+5Mn3Vz2f/zhdSKRoKGiFxH/UKkSfPEF7Nnj5rLfvdvrRCJBQUUvIv4jOho++wy2bHGz3x086HUikYCnohcR/9KoEUydCqtXw003wdGjXicSCWgqehHxP7GxMGECLFkCnTq5Y/ciclFU9CLinzp3hrfegs8/h9tv11z2IhdJY92LiP/q2dPNePfEE25c/FGjNL2tSAap6EXEvz3+uCv7IUPcjHcvvOB1IpGAoqIXEf/33HOu7IcOdWXft6/XiUQChs+P0RtjWhpj1htjNhpjUv3faozpaIyxxphoX+YTET9kDLz+Ovzf/0G/fu7YvYiki0+36I0x4cAbQDMgHlhpjJlprV2XZLk8wIOAZrkQEScsDN57Dw4dgl694LLLXPGLSJp8vUVfC9hord1srT0BTALapbDcs8CLwDFfhhMRPxcZCZMnu2vtb7/dDa4jImnyddFfBWw/5358wmOJjDHVgeLW2lm+DCYiASIqCmbOhBo13DX2ixZ5nUjEr/m66FO6LsYmPmlMGDAcePSCL2RMT2NMnDEmbteuXZkYUUT8Xp48MGcOlCnjRs9btszrRCJ+y9dFHw8UP+d+MWDHOffzAJWBRcaYrUAdYGZKJ+RZa8dYa6OttdEFCxbMwsgi4pfy53fT2155pRsX/4cfvE4k4pd8XfQrgXLGmFLGmGxAZ2Dm2SettQestQWstSWttSWBZUBba22cj3OKSCAoUgS++spdcte8Oaxd63UiEb/j06K31p4CegNfAr8Ak621PxtjBhlj2voyi4gEiauvdmWfPbuby37DBq8TifgVY6298FJ+Ljo62sbFaaNfJKT9+qs7Gz8yEhYvhtKlvU4kkmWMMd9ba9M1zowmtRGR4FChgjtmf/QoNGkC27df+O+IhAAVvYgEjypVYO5c2LsXGjeGnTu9TiTiORW9iASXmjXdpXc7d7pj9rr8VkKcil5Egk+9ejBrFmze7M7G37fP60QinlHRi0hwiomBGTNg3Tp3nf3Bg14nEvGEil5Eglfz5jB1qhtMJzYWjhzxOpGIz6noRSS43XQTfPQRfPsttGvnzsoXCSEqehEJfp06uSluFyyAjh3hxAmvE4n4jIpeREJDt27w1lswezZ07gynTnmdSMQnVPQiEjp69oRXX4Vp06BrV5W9hIQIrwOIiPjUgw/C8ePw+OMQHg7vv+/+FAlSKnoRCT19+sDp09CvH4SFueP3KnsJUip6EQlNffvCmTPw1FOu7MeOVdlLUFLRi0joevJJt2Xfv78r+3ffdX+KBBEVvYiEtv/+15X9wIFui37MGJW9BBUVvYjIgAFuN/6zz7qSf+stlb0EDRW9iIgxbov+9Gl4/nm3Zf/mm+5xkQCnohcRAVfqgwe7sh861G3Rv/66yl4CnopeROQsY2DIEFf2L73ktuxffVVlLwFNRS8ici5j4MUX3TH7V15xW/bDh6vsJWCp6EVEkjLGbdGfPu226MPD3X2VvQQgFb2ISEqMcVvyZ7fsz/6pspcAo6IXEUmNMW6LPiwMRoxwk+CMHKmyl4CiohcRScvZLfvISLf7/uRJd+mdrrOXAKGiFxG5kLMn6EVGurPyT53SCHoSMFT0IiLpYQw895wr+0GD3Ja9JsKRAKCiFxFJr7Mj6IWHu2FzT52C8eMhQh+l4r/00ykiklH9+7st+yefdGX/wQfuvogfUtGLiFyMfv1cuffp48p+4kTIls3rVCLJ6EwSEZGL9dhj7oz8Tz+FTp3g+HGvE4kko6IXEbkUDz/sJr+ZORM6dIBjx7xOJHIeFb2IyKW6/34YPRpmz4abboIjR7xOJJJIRS8ikhl69oRx42DBAmjRAg4c8DqRCKCiFxHJPHfcAZMmwfLl0Lgx7N7tdSIRFb2ISKbq1AlmzIB166BRI9i50+tEEuJU9CIima11a5gzB7Ztg4YNYetWrxNJCFPRi4hkhZgYmD8f9uxxZb9+vdeJJESp6EVEskrt2rBokbu+/oYbYM0arxNJCFLRi4hkpeuugyVL3Kh5MTGwYoXXiSTEqOhFRLLaNde4ss+XD5o0ga+/9jqRhBAVvYiIL5Qs6cr+6quhZUuYNcvrRBIiVPQiIr5StKjbmq9cGdq3h/ff9zqRhAAVvYiILxUo4EbPi4mB7t3hlVe8TiRBTkUvIuJrefLA55+7wXUefRT69gVrvU4lQUrz0YuIeCF7djeHff78MHQo7NrlJsaJ0MeyZC79RImIeCU8HN58EwoXhoEDYe9eV/5RUV4nkyCiXfciIl4yBp55Bl57zY2R37KlZr6TTKWiFxHxB717w0cfwbffuhP1/vrL60QSJFT0IiL+onNnd339b79B/fqwaZPXiSQIqOhFRPxJ8+bw1Vewfz/Urashc+WSqehFRPxN7dpuF36ePG43/mefeZ1IApiKXkTEH5UvD999988oeqNGeZ1IApSKXkTEXxUqBAsXQuvW0KsX9OsHZ854nUoCjM+L3hjT0hiz3hiz0RjTN4XnHzHGrDPGrDHGfGWMKeHrjCIifiNXLpg2De69F154AW6/HU6c8DqVBBCfFr0xJhx4A2gFVAK6GGMqJVnsRyDaWlsVmAq86MuMIiJ+JyLC7bp//nn48EN3rf3+/V6nkgDh6y36WsBGa+1ma+0JYBLQ7twFrLULrbV/J9xdBhTzcUYREf9jjNt1P2ECfPMNNGwI27d7nUoCgK+L/irg3J/M+ITHUnM3MCdLE4mIBJKuXWHOHNi2DerUgR9+8DqR+DlfF71J4bEUp2wyxnQFooFhqTzf0xgTZ4yJ27VrVyZGFBHxc02auK36iAi3ZT99uteJxI/5uujjgeLn3C8G7Ei6kDGmKfAU0NZaezylF7LWjrHWRltrowsWLJglYUVE/FaVKrB8ubv8rkMHePFFTXUrKfJ10a8EyhljShljsgGdgZnnLmCMqQ6MxpX8/3ycT0QkcFx5JSxa5Oa1f+IJ6NFDZ+RLMj6dptZae8oY0xv4EggHxlprfzbGDALirLUzcbvqcwNTjDEA26y1bX2ZU0QkYOTI4aa2veYaePZZ2LwZPvkErrjC62TiJ4wNgl090dHRNi4uzusYIiLe+vBDuOsuKFHCTY5TvrzXiSSLGGO+t9ZGp2dZjYwnIhIsbrsNFixw19jXqeNG1ZOQp6IXEQkm9eu7k/SKFHEz4Y0apZP0QpyKXkQk2JQq5Wa/a9HCjZF/7706SS+EqehFRILR5ZfDjBnw5JPw9ttw443w559epxIPqOhFRIJVeDg89xx8/DGsWgXR0bBypdepxMdU9CIiwe6WW2Dp0n9G0nv/fa8TiQ+p6EVEQkG1am5rvm5d6N4dHnkETp3yOpX4gIpeRCRUFCwIc+fCAw/A8OFuulvNFRL0VPQiIqEkMhJGjoSxY93EODVrusvxJGip6EVEQtGdd7pL8MLD3XH7N97Q9fZBSkUvIhKqatRw89k3bw69e7u57o8c8TqVZDIVvYhIKMuXD2bOdJfhTZoEtWrBr796nUoykYpeRCTUhYW5gXXmznUn511/PUyZ4nUqySQqehERcZo0cbvyq1Rx194/9BAcP+51KrlEKnoREflHsWKwaBE8/LA7O79uXfjtN69TySVQ0YuIyPmyZXPX2c+YAb//7k7a++ADr1PJRVLRi4hIytq2hdWroXp16NbNjah3+LDXqSSDVPQiIpK6YsVgwQIYMAAmTHAD7Kxa5XUqyQAVvYiIpC0iAp55xhX+4cNQuza89poG2AkQKnoREUmfmBi3K79ZM3jwQWjVCnbs8DqVXICKXkRE0q9AAfjsM3j9dVi82F2K98knXqeSNKjoRUQkY4yB++9319yXLg0dO7oT9Q4c8DqZpEBFLyIiF6dCBTcxTv/+8OGHULUqfP2116kkCRW9iIhcvMhIGDjQTXmbLRvceCP06QPHjnmdTBKo6EVE5NLVqeMuu7v3XnjpJXft/XffeZ1KUNGLiEhmyZULRo2CL7+Ev/+G+vXhkUfc1+IZFb2IiGSu5s3hp5/g3/92Q+nq2L2nVPQiIpL5LrsM3nzTDbJjrbsGv3dvOHTI62QhR0UvIiJZ58YbYc0aN+Xtm2+66+5nz/Y6VUhR0YuISNbKlQtGjIAlSyBHDoiNhU6d4I8/vE4WElT0IiLiG/XruyF0Bw+GWbOgYkU35/3p014nC2oqehER8Z1s2eCpp2DtWqhXz+3Sr10b4uK8Tha0VPQiIuJ7ZcrAnDnw8cduF36tWvDAA7Bvn9fJgo6KXkREvGEM3HIL/Por9OoFb7wB5cvD6NHanZ+JVPQiIuKtyy93s+H98ANUquSuv69Z082OJ5dMRS8iIv6hWjVYtMjtzt+7Fxo1gltvhW3bvE4W0FT0IiLiP87dnf/MMzBzppslr39/DbZzkVT0IiLif3LmhAEDYP16aNsWnn0WypZ1x/FPnPA6XUBR0YuIiP+6+mqYNAmWL3fX3ffu7Y7jT57shtaVC1LRi4iI/6tVCxYuhM8/d6Pr3XrrP49JmlT0IiISGIyB1q3dvPfvvQd//QWNG0PTprB0qdfp/JaKXkREAkt4OHTv7o7fv/yymxK3QQNo1gy+/dbrdH5HRS8iIoEpRw545BHYvBmGDXPj6NevDy1awHffeZ3Ob6joRUQksOXKBY89Blu2wIsvuoF36tWD5s1hwYKQP2lPRS8iIsEhVy7o08cV/tChsGYNNGniTtqbMiVkh9VV0YuISHDJnRsefxy2bnXj5u/f7wbhqVDB3T92zOuEPqWiFxGR4BQVBT17ulH2pk6FfPncOPolSrjBeHbs8DqhT6joRUQkuIWHw803u0F3Fi6E6693I+2VKAGdO7tL84L4OL6KXkREQoMxEBMDs2bBb7/BAw/AF1+4S/Nq1oRx4+DoUa9TZjoVvYiIhJ4yZeCVV+CPP+Ctt9z4+XfdBUWKQK9eEBcXNFv5KnoREQlduXLBvfe6QXcWLYKbbnJb9tdf76bNHTkS9uzxOuUlUdGLiIgYA40awYQJsHMnvPkmZMsGDz0ERYtCx47wyScBeca+il5ERORcefPCfffBypVuXP377oNvvnFlX6iQG373iy/g5Emvk6aLz4veGNPSGLPeGLPRGNM3heezG2M+Tnh+uTGmpK8zioiIAHDddTBiBMTHw9y5ruxnzIBWrdyW/r//DXPmwPHjXidNlU+L3hgTDrwBtAIqAV2MMZWSLHY3sM9aWxYYDgz1ZUYREZFkIiLcpDljx7pZ86ZNczPnffCBm1GvQAHo1Mnd37fP67Tn8fUWfS1go7V2s7X2BDAJaJdkmXbA+ISvpwJNjDHGhxlFRERSlz07tG8PH38Mu3fD55/D//2f273frRsULAg33ujG3V+zxvOz931d9FcB28+5H5/wWIrLWGtPAQeA/D5JJyIikhFRUW6LfvRod6nesmVu+N29e+GJJ9yu/2LF4O+/PYsY4ePvl9KWedJfddKzDMaYnkDPhLvHjTFrLzFboCkA7PY6hI9pnUNHKK631jlY7djhLuNzMmudS6R3QV8XfTxQ/Jz7xYCkgw2fXSbeGBMBXA7sTfpC1toxwBgAY0yctTY6SxL7Ka1zaAjFdYbQXG+tc2jwYp19vet+JVDOGFPKGJMN6AzMTLLMTKB7wtcdgQXWBsnwRCIiIj7m0y16a+0pY0xv4EsgHBhrrf3ZGDMIiLPWzgTeBSYYYzbituQ7+zKjiIhIMPH1rnustbOB2Uke63/O18eAThl82TGZEC3QaJ1DQyiuM4TmemudQ4PP19lor7iIiEjw0hC4IiIiQSygij4Uh89NxzrfYYzZZYxZlXDr4UXOzGKMGWuM+V9ql0saZ2TCv8caY0wNX2fMCulY7xhjzIFz3uf+KS0XKIwxxY0xC40xvxhjfjbGPJTCMkH1XqdznYPqfQYwxkQZY1YYY1YnrPfAFJYJqs/udK6z7z67rbUBccOdvLcJKA1kA1YDlZIs0wt4K+HrzsDHXuf2wTrfAbzuddZMXOcbgBrA2lSebw3MwY23UAdY7nVmH613DDDL65yZuL5FgBoJX+cBNqTwsx1U73U61zmo3ueEdTJA7oSvI4HlQJ0kywTbZ3d61tlnn92BtEUfisPnpmedg4q1djEpjJtwjnbA+9ZZBuQ1xhTxTbqsk471DirW2p3W2h8Svj4E/ELyUTKD6r1O5zoHnYT373DC3ciEW9KTw4Lqszud6+wzgVT0oTh8bnrWGeDmhF2bU40xxVN4Ppik998kGNVN2BU4xxhzrddhMkvCbtrquK2ecwXte53GOkMQvs/GmHBjzCrgf8A8a22q73WQfHanZ53BR5/dgVT0mTZ8bgBJz/p8BpS01lYF5vPPb8XBKtje4/T6AShhrb0OeA2Y7nGeTGGMyQ18AjxsrT2Y9OkU/krAv9cXWOegfJ+ttaettdVwo6HWMsZUTrJI0L3X6Vhnn312B1LRZ2T4XNIaPjeAXHCdrbV7rLVnJ0J+G6jpo2xeSc/PQdCx1h48uyvQurEoIo0xBTyOdUmMMZG4wvvQWvtpCosE3Xt9oXUOxvf5XNba/cAioGWSp4LtsztRauvsy8/uQCr6UBw+94LrnOSYZVvccb9gNhO4PeGM7DrAAWvtTq9DZTVjzJVnj1kaY2rh/u/u8TbVxUtYl3eBX6y1r6SyWFC91+lZ52B7nwGMMQWNMXkTvs4BNAV+TbJYUH12p2edffnZ7fOR8S6WDcHhc9O5zg8aY9oCp3DrfIdngTOBMWYi7szjAsaYeGAA7kQWrLVv4UZVbA1sBP4G7vQmaeZKx3p3BO4zxpwCjgKdA/mDEKgPdAN+SjiOCfAkcDUE7XudnnUOtvcZ3NUG440x4bhfXCZba2cF82c36Vtnn312a2Q8ERGRIBZIu+5FREQkg1T0IiIiQUxFLyIiEsRU9CIiIkFMRS8iIhLEVPQikmHGmEXGGJv05nUuEUkuYK6jFxG/8l9gLO5a74+AdcA+TxOJSIp0Hb2IZJgxZiqQA2hrrT3tdR4RSZ123YvIxaiEm5FLJS/i57TrXkQuxkKgjzHmILAY2AUcDILhWkWCjopeRC7Gw7g5w8fg5mEA+BdBMq2qSDBR0YvIxegCPIqbeOZTYDOwxtNEIpIinYwnIhmSMP3mH8Cb1to+XucRkbTpZDwRyahrgJzAZK+DiMiFqehFJKN2JfxZ3dMUIpIu2nUvIhlmjJkBNANeBZYC8cA6a+0JT4OJSDIqehHJMGNMduARoBNQHsgFbATqWmt3e5lNRM6noheRS5Zwgt564Dlr7Uiv84jIP3SMXkQyw1W4zxONdy/iZ7RFLyKXzBizGHctfQ9r7Smv84jIP1T0IiIiQUy77kVERIKYil5ERCSIqehFRESCmIpeREQkiKnoRUREgpiKXkREJIj9P9Nzc/o6f1svAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "name='Radious-Density Graph for n=%.1f and b=%.2f' %(n,b)\n",
    "fig=plt.figure(figsize=(8,6))\n",
    "plt.plot(x,y,'r',label='θ(ξ) with ξ1=%.4f (n=%.1f and b=%.2f)' %(x[i-1],n,b))\n",
    "\n",
    "##set plot legend\n",
    "plt.legend(loc='best')\n",
    "\n",
    "##set x and y label\n",
    "plt.xlabel('ξ', fontsize=15)\n",
    "plt.ylabel('θ',fontsize=15)\n",
    "\n",
    "##set O(0,0)\n",
    "plt.xlim()\n",
    "a=plt.gca()\n",
    "a.set_xlim(0)\n",
    "a.set_ylim(0)\n",
    "\n",
    "##save the figure\n",
    "plt.savefig('%s.pdf' %name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}