NSM_test_suite.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<module 'chem_evol' from 'chem_evol.pyc'>"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sygma as s\n",
    "reload(s)\n",
    "import omega as o\n",
    "reload(o)\n",
    "import chem_evol as c\n",
    "reload(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SYGMA run in progress..\n",
      "   SYGMA run completed - Run time: 0.56s\n",
      "SYGMA run in progress..\n",
      "   SYGMA run completed - Run time: 0.68s\n",
      "SYGMA run in progress..\n",
      "   SYGMA run completed - Run time: 0.61s\n",
      "SYGMA run in progress..\n",
      "   SYGMA run completed - Run time: 0.62s\n",
      "SYGMA run in progress..\n",
      "   SYGMA run completed - Run time: 0.63s\n"
     ]
    }
   ],
   "source": [
    "# Do a SYGMA run for each NuGrid metallicity\n",
    "\n",
    "s_02 = s.sygma(iniZ=0.02, imf_type='salpeter')\n",
    "s_01 = s.sygma(iniZ=0.01, imf_type='salpeter')\n",
    "s_006 = s.sygma(iniZ=0.006, imf_type='salpeter')\n",
    "s_001 = s.sygma(iniZ=0.001, imf_type='salpeter')\n",
    "s_0001 = s.sygma(iniZ=0.0001, imf_type='salpeter')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.109535502807\n",
      "0.108478909105\n",
      "0.108056271737\n",
      "0.107633634269\n",
      "0.107633634271\n"
     ]
    }
   ],
   "source": [
    "# Show the number of neutron star mergers at each timestep; should not be negative, and should be roughly equal!\n",
    "\n",
    "print np.sum(s_02.history.nsm_numbers)\n",
    "print np.sum(s_01.history.nsm_numbers)\n",
    "print np.sum(s_006.history.nsm_numbers)\n",
    "print np.sum(s_001.history.nsm_numbers)\n",
    "print np.sum(s_0001.history.nsm_numbers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "## IMF check\n",
    "\n",
    "The Salpeter allows one to calculate the number of stars $N_{12}$ in the mass interval [m1, m2] with\n",
    "\n",
    "(I) $N_{12} = k_N \\int _{m1}^{m2} m^{-2.35} dm$ \n",
    "\n",
    "Where $k_{N}$ is the normalization constant, which can be derived from the total amount of mass in the system $M_{tot}$.\n",
    "Since the total mass $M_{12}$ in the mass interval above can be estimated with\n",
    "\n",
    "(II) $M_{12} = k_N \\int _{m1}^{m2} m^{-1.35} dm$\n",
    "\n",
    "the mass interval of [8, 100] for neutron star progenitors, [0.1, 100] for the total mass, and $M_{tot}=1e4$ will yield for $k_N$:\n",
    "\n",
    "$1e4 = \\frac{k_N}{0.35}(0.1^{-0.35} - 100^{-0.35})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Name the relevant quantities\n",
    "mtot = s_001.mgal\n",
    "nsm_l = s_001.transitionmass\n",
    "imf0 = s_001.imf_bdys[0]\n",
    "imf1 = s_001.imf_bdys[1]\n",
    "\n",
    "# Compute the normalization constant as defined above\n",
    "k_N = (mtot*0.35) / (imf0**-0.35 - imf1**-0.35) #(I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The total number of NS merger progenitors $N_{12}$ is then:\n",
    "\n",
    "$N_{12} = \\frac{k_N}{1.35}(8^{-1.35} - 100^{-1.35})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Compute the total number of neutron star merger progenitors as defined above\n",
    "N_nsm = (k_N/1.35) * (nsm_l**-1.35 - imf1**-1.35)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compared to a SYGMA run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Theoretical number of neutron star progenitors:  74.2177993298\n",
      "Number of neutron star progenitors in SYGMA run:  74.2177993298\n",
      "Ratio (should be ~1):  1.0\n"
     ]
    }
   ],
   "source": [
    "# Compute the number of NS merger progenitors in one of the SYGMA runs (normalize to mgal)\n",
    "A_imf = mtot / s_001._imf(imf0, imf1, 2)\n",
    "N_sim = A_imf * s_001._imf(nsm_l, imf1, 1)\n",
    "\n",
    "print 'Theoretical number of neutron star progenitors: ', N_nsm\n",
    "print 'Number of neutron star progenitors in SYGMA run: ', N_sim\n",
    "print 'Ratio (should be ~1): ', N_sim / N_nsm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ensure r-process yields are being read in properly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Obtain the fractional isotope yields from a SYGMA run\n",
    "l = len(s_01.history.ism_iso_yield_nsm)\n",
    "y = s_01.history.ism_iso_yield_nsm[l-1]\n",
    "n = np.sum(y)\n",
    "yields = y / n\n",
    "\n",
    "# Exclude isotopes for which there are no r-process yields in the yield table\n",
    "nonzero = np.nonzero(yields)\n",
    "yields = yields[nonzero]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Obtain the mass numbers for the isotopes (x-axis ticks)\n",
    "massnums = []\n",
    "\n",
    "for i in s_01.history.isotopes:\n",
    "    massnum = i.split('-')[1]\n",
    "    massnums.append(float(massnum))\n",
    "\n",
    "# Again exclude zero values\n",
    "massnums = np.asarray(massnums)\n",
    "massnums = massnums[nonzero]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Hacky text parser to get all the fractional isotope values from the r-process yield table\n",
    "r_yields_text = open('yield_tables/r_process_rosswog_2014.txt')\n",
    "r_yields = r_yields_text.readlines()\n",
    "lines = []\n",
    "\n",
    "for line in r_yields:\n",
    "    lines.append(line)\n",
    "\n",
    "newlines = []\n",
    "\n",
    "for line in lines:\n",
    "    if '&' in line:\n",
    "        new = line.strip()\n",
    "        new = new.split('&')\n",
    "        newlines.append(new)\n",
    "        \n",
    "massfracs = []\n",
    "rmassnums = []\n",
    "        \n",
    "for ind, el in enumerate(newlines):\n",
    "    if ind is not 0:\n",
    "        massfracs.append(float(el[2]))\n",
    "        rmassnums.append(float(el[1].split('-')[1]))\n",
    "\n",
    "# Array of r-process yields to compare with simulation yields\n",
    "massfracs = np.asarray(massfracs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x574e210>"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtkAAAHwCAYAAABzBnP9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XGcm3WZ7/3vVZDusC8KtbUo0plRCitCWfXsAXcWJF18\npGDZ1oUqQRCBHj3lKbyKPH1A1zpoGU+1HjrYCivbisJKcMvuFjmy6qrcnoX49GHP6qEP0AVcZwZQ\nWwaQUkgLttfzx510Mpkkk2TuJHeSz/v1yqu9f7lz378kk8w1v/v6XT9zdwEAAACIzrRmdwAAAABo\nNwTZAAAAQMQIsgEAAICIEWQDAAAAESPIBgAAACJGkA0AAABEjCAbAAAAiBhBNgAAABCxQ5vdgWqY\n2eGSbpG0T9JP3f2uJncJAAAAmKDVRrL/UtIWd/+kpL9odmcAAACAYpoaZJvZZjPbaWaPFLQvNLMd\nZvaEmV2Xd9exkp7O/n9/wzoKAAAAVKHZI9m3Szo7v8HMpknamG0/SVLSzN6RvftphYG2JFmjOgkA\nAABUo6k52e7+oJn1FDSfKulJdx+WJDO7W9JiSTsk/aOkjWb2QUn3FTummXkduwwAAAAc5O5FB36b\nPZJdzFs1lhIiSc9k2+Tur7r75e7+f7p7qtQB3J1bDbf+/v6m96FVb7x2vHa8dq1z43XjteO1a61b\nnF+7cuIYZAMAAAAtLY5B9rOSuvO2j822AQAAAC0hDkG2afwkxoclzTOzHjM7TNKFkr7blJ51mEQi\n0ewutCxeu9rx2tWO1642vG6147WrHa9d7Vr1tbPJ8knqenKzuyQlJM2StFNSv7vfbmbnSBpU+EfA\nZndfW8UxvZnPCQAAAJ3BzOQlJj42NciuB4JsAAAANEK5IDsO6SIAAABAWyHIBgAAACJGkA0AAABE\njCAbAAAAiBhBNgAAABAxgmwAAAAgYgTZAAAAQMQIsgEAAICIEWQDAAAAESPIBgAAACJGkA0AAABE\njCAbAAAAiBhBNgAAABAxgmwAAAAgYgTZAAAAQMQIsgEAAICIEWQDAAAAESPIBgAAACJGkA0AAABE\njCAbAAAAiBhBNgAAABAxgmwAAAAgYgTZAAAAQMQIsgEAAICIEWQDAAAAESPIBgAAACLWUkG2mb3N\nzDaZ2d81uy8AAMRSJjNxu1gbgLpqqSDb3X/l7sua3Q8AAGIpk5Hmz5eCYGx73rzwlgusgyDch0Ab\nqKtDm3FSM9ssaZGkne5+Sl77QkmDCoP/ze7+pWb0DwCAltTVJW3aJC1dKm3ZErbt2SOZSdu2hdu5\n+7q6mtdPoAM0JciWdLukDZLuyDWY2TRJGyWdJenXkh42s3vdfYeZXSLp3ZLWuftvJFkT+gwAQPwl\nEkr39yu5YJEkKbVhrSSN2+5LJJrVO6BjmLs358RmPZLuy41km9l7JfW7+znZ7eslef5otpm9UdKA\npPdL2lRspNvMvFnPCQCAOOg5+jiN7FojSeqes1qSxm0P7/zl+AdkMoxsAzUwM7l70cHfZo1kF/NW\nSU/nbT8j6dT8Hdz9BUnLJztQIpFQb2+vent7lUgklOAvdgBApwgCaXR0bPu558bfv2uXNDgorVw5\ntv+yZdL27QTawCSCIFAQBBoaGtLQ0FDZfeMUZEcmyE34AACgk2Qy0rJlSt08oOS6Acldqcx0SVKy\n60bJTKkll0vXXiudeKI0fTo52kAVCgdvzUpnMMcpyH5WUnfe9rHZNgBRK7w0nKsyUNjGL12gtXR1\nSdu3q6+rS8MrVoRt2c/3cN7nOd3bq+TC8yWRow3USzNL+JnGT2B8WNI8M+sxs8MkXSjpu03pGdDq\nipXmyrVR4gtob4V/HHd1TWhLrt2oEd2mEd2m5Jr1Dewc0DmaEmSb2V2S0pJOMLMRM7vM3fdLukrS\nDyU9Kulud3+8Gf0DWlphEC2ND5rzS3wFQVjWa88e6ZVXwv8HQXjfpk2MZAOtoNwf1cUU5myPjo7/\nvgAQiaZVF6kXqougIxWmdgSBdMEF0j33hNu5nMu8S8LpjRuVvOp6SXklvvK2+3KXmgHEV+6P6k2b\nxj7f5SYyZvdPr1yp5LpbJUmpVcvVNzjIxEegBuWqixBkA62uxC/Z9NKlSo6Go1nFguaqS3wBiKfc\n1afc4jNF/qgep9h8C+ZgADUpF2S31LLqAIooTP/I/sJN7u8qnXNZrMRXfpkvLh8DrSO7+EzPgkXq\nWbBI6f7+0gG2VDyYJsAGIsdINtAmxqV/XH2Fkhu/qZED4eXgbn1Sw+vXhHVxsxMd06+9pmTXm8L9\nM2GAnex6U1jii8vHQEspvDLFlSigMUgXATrAuF+yh1yp1OCNYzmXS96vvo0bpfvvD+vi5vK1c6Nd\nlPADWlcQqOesxWN/VE9bruEf31t+NBtAJFplxUcAtSpM/3BX38kna3j40YNNZevicvkYaE2Fi89I\nSq0aYAVHIAYYyQZaXYXVAricDLQpJjICTUO6CNDuJvsly+VkAAAiR3URoN2VS/fIu5zc3R3eUjdn\nLyeXWrCisD2TKd4GoPVVu5gNgIowkg10gmouJxfW3c4tuy5JTz0VPqbcYhcAWke1i9kAGIeJj0Cn\nq2ZiY37d7dziFnv2SGbhsuvS2H38AgZaW7HPO59vIBIE2QAmyi5ukVywSFLesut5233kcwPtocjn\nnc83MHWkiwCtqs4VBVh2HegcVB8CasPER6Dd5PIo85c+D4KwLYoJSyy7DnSOws87n28gEoxkA60q\nCIrnUU71Mm9h3W13ll0H4iaqK1kV1tmPFeqCI0aokw20qfTGjUpedb2kbB7lihXRHLjwFxbLrgPx\nEXVFkFYKWqmGgpghyAbaFHmUQIeq15WsVtDJz73B0um0ksn/IklKpf5GfX19Te5R/FDCD2hHpfIo\n+UUDtL9OrgjSyc+9wZLJ/6KRkb86+P/h4Ueb3KPWwkg20IpaMY8SQKQ6+UpWJz/3RurpOelgkN3d\nPUCQXQQj2UC76eqStm9XX1eXhvPzsK+4ggAb6ASdfCWrk597g6VSfzMuXQTVYSQbAIBW0slXsjr5\nuSOWmPgIAEAro+LPmFaqhtJhOnGiJIvRAADQqootPrVt28TFpzolyCz2PDvlucdcbqLkyMhfHQy2\nOxlBNoDqFVtVMoqVJgFM1NUV1oVeujQMtHMl7DZtIrhES0un0+rpOUk9PScpnU43uzuRI8gGUJ16\nL+kOYKJs2bqeBYvUs2CR0v39TPRD7KRSf6Pu7gF1dw9UNFGy3Ue+yckGUD0WgwAarm3L1pFj3bHa\noUQgOdkAosWoGtBYpcrWtTqujHW0ake+W01LjWSb2WJJH5R0hKRvuPs/F9mHkWygAdp2VA2Im3Yv\nW8eVMbSwtivhZ2ZHSVrn7hMSeAiygQYIAvWctVgjB8Jf+N3Tlmv4x/fySxGolzZPqUhv3KjkVddL\nyi6Tnr/IFhBjsQuyzWyzpEWSdrr7KXntCyUNKkxj2ezuXyrx+K9I+lt3/0WR+wiygXpq91E1AA3H\nlTG0qjgG2adL2iPpjlyQbWbTJD0h6SxJv5b0sKQL3X2HmV0i6d2SviLpakk/dPeflDg2QTZQb20+\nqgaggbgyhhYWu4mP7v6gpBcLmk+V9KS7D7v765LulrQ4u/+d7v4pSecrDMIvMLNPNLLPAPKwGASA\nKGQy0rJlSt08MDYB7uYBadkyJj6i5TUtJ9vMeiTdlzeSfb6ks939E9ntiyWd6u5XV3lcP/PMM9Xb\n26ve3l4lEgkl+GsYAIB44soYWkgQBAqCQENDQxoaGtJPf/rTkiPZhza6c40QtENZIwAAmqHRQS9X\nxtBCCgdvzYrG15LiVSf7WUndedvHZtsAAEAjULcaiEwzR7Ite8t5WNK8bBrJbyRdKCnZjI4BANCR\nurqkTZuK161mdBmoSlNGss3sLklpSSeY2YiZXebu+yVdJemHkh6VdLe7P96M/gEA0LFY0RWIREsu\nRlMOJfyAJmDiEtBWqFsNVCZ2JfwA1KBYPmQcciTJ4QTaSxBIo6Nj26Oj4z/f9Vb4vZHJFG8DYo4g\nG2gFcQ5k83M4gyC8LV0atjGSDbSWZtetLvyuy2SkefPCW+78cfnuAyZBugjQKnLBa+FkpJjkSqY3\nblTyquslSakNa9W3YkWTe4QJSOtBJZr9c1L4Xbd4sWQmbd0absfsuw+dLXbLqtcTQTbaWZwDWXI4\nYy43Qrhp01hwEgThCOX27QTaiJXC7zpJsf3uQ2cjyAbaRGwD2SBQz1mLNXLgVklS97TlGv7xvYw0\nxU3Mr4Z0tMKR4lwqRGFbh/wxVPhdJyme333oeEx8BNpBsycjldLsHE6MmWzCWCKh9Gc+Q2m2uCEP\nebzC77rnngtvOXH57gMmwUg20Aqyv4TTK1cquS4cLU6tWq6+wcF4XOpvdg4nJqaD5AI1SXrqqfC9\nCAL1vH+JRvbfIokRwVgJAqWXLFFy+ixJUmrPb9X3hjd0Xh5y4Xedu1KZMMBOdr1JMovXdx86Huki\nQDsgkMVkikwYS+/fr+QfHi0pDNySmUM04n8tqURaDz9TTVOYIpFafU1n5iGTOoMWQpANAB2icMJY\ncs36scDtkCuVGrxx3Ahh32GHjRvpHjcRkmCnochDBloPOdkA0CGSa9ZrRLdpRLcpuWb9+DtnzVLf\nihUaHn5UwyOPqW/LFum116Rt2ybWNydPuLGCQKl9z6t7zupwFHvPb8lDBlrcoc3uAACgCuVGl4tM\nGEsdfkDJ7Khoat/z4T659JBEQun+fiUXLArv37BWfbn78hcZyqWf7NkTpp/0nnzweH1btzKSPVXZ\nycN9W7dqOC+fPnX4a0p23ZjNQx6g3CLQYkgXAYBWUW5y4/bt0qmnKv3BDyp5yx3Sm9+s1N7R8ukg\nmrwsZNn0E1IYokNqDtCSyqWLMJINRIWJiai3EqPLMpMeeUT62tfUd/HFGv7nfxwLwnOPk8K2/JHQ\nUmUh8yZC5tJPwv+vruOT63CF3xPFvjf4LgFaCjnZQK3y81BzI4w/+MFY2w9+MDFfldxVTFU2xeNg\nreuBAaVvvDHcXnj++NrXXV2lg7di9c2/csP4+ubF0k9e2TmWN5xLP4mDYp8tPm8Amoh0EaAWxZao\nHhyUrr1Wuv/+cPvcc6X//t+llSvDbZawRkQiq0KRf6Ul9zP9ta9JZ599MBUl/dprYX1iafJqJFNV\n69UglowH0CSkiwBRK3bZfmBA6RUrlFx4viQpdfUK9Q0MSO96V3h/bt9O+IVP6kz9FFsNL1+RlI+S\n8t+PYj/T+/ap7557wsl40uTpJ1N536cSKBfreyd93gDEk7u31S18SkBjPLRhg3frD71bfxj+f87b\nXfq2S9/27jlvn3B/R3j1VffjjnN/4IGxtgceCNtefbVZvWoP2df2oQ0bvLv7nd4990R/aPZsf2j2\nbO+ee6J3d78z/Dmbwmtd889sFO/7Aw+4z54d/pv//wp15OcNQFNl486iMSnpIsAUjLtsP/N66aWX\nNHIgXPa8e9py6cgjNfLi2nC7kyoxFK482ClLQjdCnatQTFZtpKzC9/2CC6R77hl73wv7WqTv6fXr\nlfxUtuRglascTqnvAFADFqMB6qHwsv2LLyq14uNjk8hWfFx68cWx+ztpMYnCyXn5k/EwNcUmMlZS\nmaISpaqNVKrwfbe83zuFi9mUWNwmuerzpRfTKZzImDtOJX1nYiSARis1xN2qN5EugkYovGzf/U5/\n6Kabxi6N5+6/6aax+6d4Gb/VFKbOoEbFfl7q8TNU7Ge6hp/Zce/7zLnj0z9mzHA/8siy2912xNjj\np80YSxcpTEd59VX3Y44Jb88/H/b96qu9+9CjwjSa/L6TwgSgTkS6CBCBai/Td/LkvyBQz1mLx6XO\nDP/4Xkazq9XoqhnV/swW3lfkfU+t+LiSX90sKUz/kDRucZtx27O7pP5+JdeFj0+tWq6+wcGx51qY\njrJ4cVgjfOtWad8+6eKLx9KSCj+fk6UwdfLnFUDNyqWLNH3kOeqbGMlGPTASVrmIRkQ7VuFrNMXJ\ngHVTYmS5cBJm9yFHjruiUXiFY9IrHgWvR+HkxqKTHUt8Xh+aPbv4xEg+3wBqJEaygRoUGaVjMl+F\nGBWsTYmR6/TSpUqOhiOz1U4GrKvJPhOFI9v2XyVJI/7XxbcruOJRcY3wIn3r2d9VeiIyn28ANWAk\nGx3joYceGhs9feih2g9U7UgYEJUiI9fdM+fGNr+9ZNm8ycoNVlJ+sMiofve0GWOvhR1ROoe7sG9X\nXz3+sQX7uns4h4LPN4AqiJFsdIqenpM0MvJXkqTu7gENDz9a+8GqHQkDIpLeuHEsT/nqK5Tc+M3Y\n5reXLZs32TyGcvMaCkf1S6xAKSncNpuQwz2ub4dcqdTgjWG+t3vR1St73r9EI/tvKf5cosbVHqAt\nUMIPqEVhObKLLpJeemns/k4qyYeGSq5ZP1bG7mvfUurmgbHSkDcPhBMf41B+brKyeZOVGyxXfjB/\nFccgkLZtk157TX1btmh45DENjzymvpER9Y2MhNvDj4ZpNPmTJPP75q6+k0/W8PCj4WO3bJFeey08\nbhCEkygPHCj9XKYq//3K/QHxgx+MtQVB2BaH9xVANEoNccfxJukdkm6V9B1JV5TYJ6ILAGhFkaWL\nZI2blHXIkUzmQ/0VpkQUSWuIxc9cgya41rSKY4V9G3fs2bPr91yKpZ+tX+8+bZr7978frwmtAKqi\ndksXMTOTdLe7f6TIfd6KzwkxVEkZOi7vIkrZEc70ypWly9jFSQNSHmpexbGCvk167ChLchZJP0tf\ndNG48oaxmdAKoGKxSxcxs81mttPMHiloX2hmO8zsCTO7rsRjz5P0PUl3N6Kv6CCFl3OXLVPqKzeU\nv0xfTTDBinOYTFeXtH27+lasCNMaClMg4qZYn6Ls51RWoJysb5UcOz/Anj9//P3VpncUWQU1eff/\nKL26JYDWV2qIu543SadLepekR/Lapkl6SlKPpDdI+oWkd2Tvu0TSTZLekrf/vSWOHfGFAHSEYpdz\nv//9iZeLa710TB1eoDr1TEep5dgR1CsvXA1z0rQgALGnOKaLmFmPpPvc/ZTs9nsl9bv7Odnt67Md\n/1LeY86U9JeS/kDS4+4+WOS43qznhBbT6DrY1OEdQ2UFVKKePyc1HHtc1Zdq0zsK08/0SaWuvlzJ\nrT8KjxfntCAAJZVLF4lTkH2+pLPd/RPZ7YslneruV1d5XD/zzDPV29ur3t5eJRIJJToxiEF5TVr0\no6pf0u0aiDZ6qfC4maysXbF9EAtVlSvMbyuWa79ymfq+9rXxP/O870DsBUGgIAg0NDSkoaEh/fSn\nPy0ZZDezUkiPxqeLnC/ptrztiyV9tYbjTm3cH52jCYt+TLqEdE67p5fEdanweiuxFLkfc8zY+9pO\n73Mrm2whnPz0jko+r8XeT95joOWpTLpInOpkPyupO2/72GwbUB+NroNdzSSuwhrBuVSTTZvaY6Sr\nyCSwtkmbKTfBtVjt5z17pFdeGavX3E7vc6sqnOiYyUgf/ahSbzxM3XNvnDgRupLPa70niQKIn1LR\nd71vknolbc/bPkRjEx8PUzjx8cQajhvx3yhoZw2rg13jJK6aagS3iIpH9VtJhVcgCt/Xdn6fW1bh\n1ZZZsyatV877CHQexW3io5ndJSkhaZaknQonPN5uZudIGlRYaWSzu6+t4djejOeEFtToOtg15FjX\nXCM47ip57eNqsvexcILrBRdI99wz7rkVvq+S2vN9bnHVTnRs288rgJJiVyfb3S9y92Pcfbq7d7v7\n7dn2f3L3P3L342sJsIGK5epgT7ZcdZSXc6u9XDyVGsFxVulrHxfVLoddmApjBd+9g4PSrl1j2889\nF95y2uV9bgPjlrefrI51u35eAdSsJVd8LIeRbFQsztU7Wm3lv2rF+bXPV6wSyuCgdO210v33S9On\nTz5SPfN6DR+SCUe29+2Tzj1X6RUrwtJt7kplwgA72fUmyay93udWVs3Vlnb/vAIoKZYl/OqFIBtt\no1UC0Uq00nOpoH76uOWwZ3epL7/e+eCgeq5ZrRF9XVIYnKVWfHxs/5vWqO+aa8afT6KEX5zUEjS3\n0s84gMgQZANS/H4Jxq0/9dJKdbErrJ+eXLO+spFqhcFZcuVnNbL/lnB/cnVbQ6d8PgFMSexysoGG\nKyzJJU3Mpe3k/tRTK5UjLNHX5P6usdzcz315fO7tSy8pfdFFYQ72wvOV/spX1HfzzRoeflTDw4+q\n7+STpfw//MnVbQ2U3AMwRYxko3PEbVnzuPWnzqa0JHWDjevr1VcoufGbZZfDLjlSTa4uALQ10kWA\nrLgFeuP6U5irK7XV5enYlDerIA1gXF8PuVKpwRtLL4c92QQ50g7Q6vgZBkoiXQTIqqokV6P7s+rz\n7Zs+EpfyZpWk6RT21V19J588lv5xzTVjAXYl5QhJO+hM5Vb+bCWdlNoGRK3UKjWtehMrPqKUBx7w\n7mkzxlYZnDZj4gpuzeyPHeF+5JFjK8zlVptrdZWsdlls1csoVtospnAlv/zXuZaVORvZd8RHufe9\nwpU/W0a5zwzQ4VRmxcemB8VR3wiyUVSNy5o3uj8PzZ7dnssyNzMgKXKMh266qfTrTNCMyVTyM9tm\ngSlLxgPFlQuyyclG54hbXmGRc8cmb7nRgkDpJUuUnD5LkpTa97z6tm6d+hL3JUry9bx/SfmSeuSb\nYjIVTFyO2xyQqejY7yZgEuRkA1L8cmMLzx2XvOVmSCSUnD5LI7vWaGTXGiX3FNxfaw5osZJ8S5ZI\nBw6M7bNrV7iK41TPhc6SSCjd3x+WblywSOn+/gmVgeI2B6RmnfzdBEzBoc3uAFA3cRu5LidvAl1y\n3YAkKbVqIJ4LtjTCjBnFRwlreR2ywVBywSJJ4QqNqa/2j73OSy4Pl0k/8cRwmfSpnAsdJRdEh/9f\nreH8kepSgWmrlejkuwmoGekiaE+ttMpgTiv9URC1IukiuvTSsaXIp3ipfbJL3en165X81OpIzoUO\nUa50Y7vVR+/k7yZgEtTJRmfqsMVeWlatedOVmqyOtcg3RZUqCaIJTIGOUC7IJl0E7aswTWDDWvUR\nYMdPV1fx0T0vshR5te9fJZe62+WyPhon+zPb19U1PkXkiivGfo7jNgcEQMMxko22xghlC4r6Unu5\nEcV2u6wPAGgo0kXQmSpIE0BMTeVSe7WP5bI+AKBGlPBD56lkuWvEVyVBb6llq6tdAprL+gCAOmAk\nG+2LEcr2MFmlGGn8e8qEVwBAgzCSjc7ECGV7KLagzNKlYZs0ceRaUloqu0gIAAD1xkg2gJZQconq\nIjW2k9NmaOTFtZKY8AoAqB8mPgJoeeUqxYy7b+b10ksvMeEVAFB3pIsAaG2lalkXs3s3E14BAE3H\nSDaAeCusZe2uVOY59R12mPTUU9K2bUovWqTk3kOlWbOU2juqvnsLRq6Z8AoAqAPSRdCeqB7SOQrf\n19zkx7/9W+nyy6U9e6RcYJ1feYSfBQBAHZEugvZTSz1ktK7CYDmRULq/Xz0Lz1fPr19SemBgbOQ6\nkSDABgA0XcuNZJvZ4ZJ+Kqnf3e8vcj8j2Z2CesgdrdxESAAAGqHdRrKvk/SdZncCTZI/Sp0bzaQe\ncuepZiIkAABN0JQg28w2m9lOM3ukoH2hme0wsyfM7Loij3u/pMckPSep6F8NaGNFUkSSn/uyRnSb\nRnSbkmvWN69vaJxMRlq2jAoiQKso/FxmMsXbgDbTlHQRMztd0h5Jd7j7Kdm2aZKekHSWpF9LeljS\nhe6+w8wukfQeSTMkvSTpJEmvuvuHihybdJF2lp8i8otfqOea1RrR1yVRD7mjMOkVaA25wZFNm8Lv\n5kxGmjcvvO+pp8LPLJOV0cLKpYsc2ujOSJK7P2hmPQXNp0p60t2HJcnM7pa0WNIOd79T0p25Hc3s\nY5JGhc6TTRFJLlgkSfrM0nP1xW0DkqTUqgG+qDtFsfeX9xxojnJ/9HZ1hQF2/vyZPXskM2nbtnA7\ndx+fYbSZpgTZJbxV0tN5288oDLwncPc7yh0okUiot7dXvb29SiQSSjCy2VaSa9ZrRLdJkr7404IJ\nb1dcwRc1ADRK4Ui1NHFkumBwJLVhrSSN2+7j9zRaRBAECoJAQ0NDGhoaKrtv06qLZEey78tLFzlf\n0tnu/ons9sWSTnX3q6s8Luki7aRIfeSesxazZDYAxEUFlZ4KqwFJojoQ2kKrVBd5VlJ33vax2TZ0\nqsKJjpmM9NGPKvXGw9Q990YmvAFAHExW6amwGtBzz4W3HKoDoU01cyS7V+FI9vzs9iGS/l3hxMff\nSPp/JSXd/fEqj8tIdjspHCG54ALpnntYMhsAYqRk3frsYEl65Uol190quSuVCQPsZNebJDOlVi1X\n3+Ag82nQkmK3rLqZ3SUpIWmWpJ0KF5a53czOkTSocIR9s7uvreHYBNltJr1xo5JXXS8pm7u3YkWT\newQAOGiyNL7CgZDclcfCNgJstKDYBdn1RJDdfljZDwBiqnCkWmJkGh2FIButi4mOABBv1K1HByPI\nRmtihAQAAMQYQTZaB7l7ANDaGNlGB2mVEn7odIUl+6RwRbD588eX6OOLGgCap1jJ1Fxbse/xIJj4\nPQ50AIJsxEf+8rtBMFa+b9MmAmsAiIPJgujs93h6yRL1HH2ceo4+TuklS/geR0ciXQSxQ8k+AIix\nIFB6yRIlp8+SJKX2Pa++rVvLrvBIVSi0K3Ky0VL4cgaAeJvse5rvcXQKcrLROgqX32W5XQBoLUGg\n1L7n1T1ntbrnrFZq3/N8j6MjMZKN+KBkHwDEX7l0kVzO9qZNY+kjQSAtW8b3ONoS6SJoHZR+AoD4\nqiSI5nscHYQgGwAARIMgGjiIIBsAAACIGBMfAQAAgAYiyAYAAAAiRpANAAAARIwgGwAAAIgYQTYA\nAAAQMYJsAAAAIGIE2QAAAEDECLIBAACAiBFkAwAAABEjyAYAAAAiRpANAAAARIwgGwAAAIgYQTYA\nAAAQMYJsAAAAIGIE2QAAAEDEWirINrMzzex/mtmtZva+ZvcHAAAAKKalgmxJLullSdMlPdPkvgAA\nAABFNSVfm8nvAAAgAElEQVTINrPNZrbTzB4paF9oZjvM7Akzu67wce7+P939g5Kul/SFRvUXAAAA\nqEazRrJvl3R2foOZTZO0Mdt+kqSkmb0je98lZnaTmb0lu/vvJB3WwP4CAAAAFTu0GSd19wfNrKeg\n+VRJT7r7sCSZ2d2SFkva4e53SrrTzD5kZmdLOlJhQF5UIpFQb2+vent7lUgklEgk6vNEAAAA0DGC\nIFAQBBoaGtLQ0FDZfc3dG9OrwhOHQfZ97n5Kdvt8SWe7+yey2xdLOtXdr67yuN6s5wQAAIDOYWZy\ndyt2X6tNfAQAAABiL05B9rOSuvO2j822AQAAAC2lmUG2ZW85D0uaZ2Y9ZnaYpAslfbcpPQMAAACm\noFkl/O6SlJZ0gpmNmNll7r5f0lWSfijpUUl3u/vjzegfAAAAMBVNm/hYL0x8BAAAQCNMaeKjmf2l\nmT1pZi+Z2W4ze9nMdkffTQAAAKA9TDqSbWZPSTqvVVI3GMkGAABAI0y1hN/OVgmwAQAAgDioZMXH\nfzWz70jaKmlfrtHd/6FuvQIAAABaWCVB9gxJr0r6QF6bSyLIBgAAAIqguggAAABQg6lWFznWzP7R\nzHZlb39vZsdG300AAACgPVQy8fF2hSsvHpO93ZdtAwAAAFBEJSX8fuHu75qsLS5IFwEAAEAjTLWE\n3/NmdrGZHZK9XSzp+Wi7CAAAALSPSoLsyyV9WNJvJf1G0gWSLqtnpwAAAIBWRnURAAAAoAbl0kVK\n1sk2s//b3b9sZhsU1sUex92vjrCPAAAAQNsotxhNbin1f21ERwAAAIB2UTLIdvf7sv991d235N9n\nZkvr2isAAACghVVSwu/f3P09k7XFBTnZAAAAaIRac7LPkXSupLea2Vfz7poh6ffRdhEAAABoH+Vy\nsn+tMB/7LyT9r7z2lyVdU89OAQAAAK2sknSRGZJecff92e1DJE1391cb0L+qkS4CAACARpjqio8/\nlNSVt90l6UdRdAwAAABoR5UE2X/g7ntyG9n/H16/LkUgkwlvhW0AAABAA1QSZL9iZgcriZjZf5IU\n74h13rzwlgusg0CaP59AGwAAAA1RbuJjzkpJW8zs15JM0pslfaSuvZqqPXskM2nbtnB76VJpyxap\nq6v84wAAAIAITBpku/vDZvYOSX+Ubfp3d3+9vt2amvTAgCQpuWCRJCm1Ya36Eokm9ggAAACdZNLq\nIpJkZidLeqekP8i1ufsddexXzczMu+e8XZI0smuNJKl7zmoN7/xlM7sFAACANjOl6iJm1i9pQ/a2\nQNKXFdbOjq/nngtvOaOjYV52KcVytcnfBgAAQI0qmfh4gaSzJP3W3S+T9MeSjqxrr0qw0I1m9lUz\nu6TUfqlZ05WaNV3dc29Ud/eAUjcPSMuWlQ6m588fH4QzURIAAABTUMnEx4y7HzCz32cXptklaW6d\n+1XKYknHShqV9EypnfpGRiRJw/kTHa+4ovjEx64uadMmpZcsUXL6LElSat/z6tu6lYmSAAAAqEkl\nI9n/amZHSfobhcur/5ukn03lpGa22cx2mtkjBe0LzWyHmT1hZtcVeegfSXrI3f8vSVeWPEFX18QA\nuVzAnEgoOX2WRnat0ciuNWGwzURJAAAA1KjsSLaZmaT/5u6/k/TXZvZ9STPc/ZFyj6vA7QpzvA9O\nnjSzaZI2KkxN+bWkh83sXnffkU0Nebekn2usRvf+KfYBAAAAqIuyQba7u5ndL2l+dnsoipO6+4Nm\n1lPQfKqkJ919WJLM7G6F6SE73P1OSXeaWZekDWZ2hqSfRtEXSVIQKLXveSXnrJYUposoCBjNBgAA\nQE0qycn+NzP7z+7+cJ378lZJT+dtP6Mw8D7I3TOSlk12oEQiod7eXvX29iqRSChRLljOZKRly9S3\ndauGc/sFQThRcvt28rIBAAAgSQqCQEEQaGhoSENDQ2X3nbROtpntkDRP0rCkVxSu+ujufspUOpkd\nyb4vdxwzO1/S2e7+iez2xZJOdferqzyuT1r7O5MZHzznqogUtuW2C/cv1QYAAICOUVOdbDN7W/a/\nZ0s6TtKfSzpP0qLsv1F7VlJ33vax2bZoFSvZt23bxJJ9+QE2Jf4AAABQhXLpIvdI+k+SvuHuZ9Xh\n3Ja95TwsaV52hPs3ki6UlIz8rNmSfVq6VNqyJWzL/b9Mib+K9wcAAEDHKxdkTzOzz0g6wcw+VXin\nu99U60nN7C5JCUmzzGxEUr+7325mV0n6ocIR9s3u/nit5ygrkVC6v1/JBYskSakNa9WXn7ddmAqS\nSCj9mc+U3h8AAADIUzIn28z+SNISSSsl/XXh/e7++fp2rTZFc7KL5E/3HH2cRnatkSR1z1mt4Z2/\nHNt3/vxw9DpvImTP+5doZP8tE/cvcw4AAAC0r5pyst393939S5Iud/fPF97q1tuoFcupHhyUdu0a\n2x4dHbs/Pz0kCMLbkiXSgQNj++/aFR4jhxxtAAAA5Jm0ukirKTqSHQRjedT79knnnqv0ihVKbv2R\nJCm1arn6BgfHlexLb9yo5FXXh/fP7pL6+5Vcd2u4veT96tu4Ubr/fmn69LFjk0JSPSq3AACAFlVu\nJLszgmwVBM03rVHfNdeM36EgsCuZTpI73vr1Sn4qu3jNhrXqW7EiyqfRGUqk5lCjHAAAtAKCbE0e\nNI8TBOo5a7FGDoQj193Tlmv4x/eOG6mu6ngoLf8qg8RVAQAA0DJqysnOe/BSMzsi+//Pmtk/mNl7\nou5kXQVBmHedk5+DXSi7AmTq5gF1d4e31M0D4ehqLue6muOhvGyll54Fi9SzYJHS/f0E2AAAoOVV\nsuLjI+5+ipmdLulGSeskfc7dT2tEB6s1YSQ7m5KQXrlyLKe6SA72OOXyhGs5HsriqgAAAGhFU0oX\nMbOfu/u7zey/Sdru7nfl2urR2amqtITflCbXVbsse9TaabJgBak5AAAAcTSldBFJz5rZ1yV9RNL9\nZja9wsfFR6mVHKM4Xm7y3rZtY231LOnXTsu8V5KaAwAA0IIqGck+XNJChaPYT5rZWyTNd/cfNqKD\n1So18bGuGj15r50mC7bTqDwAAOgoU00XOU7SM+6+z8wSkk6RdIe7/y7ynkagKUG2CkoENqCkX6PP\nBwAAgPGmGmT/QtKfSOqVdL+keyWd5O7nRtzPSDQryG705L2mTRZsdj46AABATEw1J/uAu/9e0l9K\n2uDuqyS9JcoOtrxGl/Qrd75iucxR5TcX5oNnMtK8eeEtv7xhfn54PfsDAAAQU5UE2a+bWVLSxyT9\nj2zbG+rXpRYTxeS9agLRcud74YXJJ0VWG/Tm39fVFa7OuHRpeNxt26RXXpH27An/n8sV37RpXLnD\ntpikCQAAUIVK0kXeKem/SvqZu6fM7G2SPuzuX2pEB6vVlHSRqUzeq2Vp8XLnKzcpstpzldg/vXSp\nkqNhkJzasFaSSueHF/bngguke+4ZP0mT9BIAANCCWFa92SYLwiOuFlJ2UmS15yqyf8/+Lo28GAbX\n3XNWS1LZ/PD0+vVKfircLzW7S33555vsDwoAAICYmuqy6seb2T1m9piZ/UfuFn0321QlKRMRLy2e\nXLNeI7pNI7pNyTXrx99ZeK7PfGbiufJTOQr3v+gi6aWXxu5/7rnwlrNrlzQ4OO65Jld9fqw/+7vG\n0k0K00sAAADaRCXpIg9K6pe0XtJ5ki6TNM3dP1f/7lUvliPZFYwe11wtpHCUvIIVFMed65ArNfyj\nrWVHlgv3Tw3eGC4p765UJgywk11vksyUWvJ+9W3cKN1/vzR9urR4sXpedo34Xx/sT2rFx5X86mZJ\nlB8EAACtq9xI9qEVPL7L3X9sYfQ6LOkGM/tfkmIZZMdSdjQ4uWCRpGxgmT96XKpayGSj2YU505mM\n9NGPKvXGw5TsujEMelcNjA+aC8914IC0eLF0773hdu4PgPxUlvz93dV38skaHn50rA+ShvMC/XRv\nr5ILzw+f6+wupb7ar+S6gXB71YCSKz+rEd0mSUquWa1hgmwAANBu3L3sTVJaYVrJP0haIelDkv59\nssc16xY+pSZ79dUJTd1z3u7St136tnfPefv4fY87zh/asMG7u9/p3d3v9Ic2bHA/7riix5nggQf8\noSOP9O45b/fuOW/3h2bMcH/ggeL9KXGuh2bP9m79oXfrD8NzT7FvJZ9rtr/d02aM3T+tSH8BAABa\nQDbuLBqTVpIu8p8lPS7pKElrJB0p6cvu/v/UK/CfiqanixSryDE4qJ5rVmtEX5dUJIVjikuLl001\nqWDxmEkfX7B/2b6VS1fJvjbplSvDdBNJqVXL1Tc4yMRHAADQcqY08dHdH3b3Pe7+jLtf5u5/GdcA\nOxYKa0n/4AfStdcqdfXlpetoFwsuowg4i0263LZt/KTLwnSQwomLhfuX69tkNcO7uqTt29W3YoWG\nhx/V8PCjYT42ATYAAGgzJUeyzey75R7o7n9Rlx5NUdNHsrPGldG7aY36rrlm/A7lRoOrGdkOAqWX\nLFFy+qzwXPueV9/WgomMk9TNHjeyXDhxsdpyglMclQcAAGgVNdXJNrPnJD0tKSVpm6RxB3D3n0bc\nz0jEJcieUrWQSheMqXDfsnWziwTA4+paU/0DAACgqFqD7EMk/R+SkpJOkfQ9SSl3f7ReHY1CLILs\nycroRbk4TQUjx9UG/DX/gQAAANBBasrJdvf97v59d79U0nslPSUpMDOGNcuZLC856sVpJsvnLlUe\nsJRq9wcAAMAEZauLmNl0SR9UOJrdK+m7kr7h7s82pHc1iMVIdgQj1ZGMJldbzYPqHwAAABWrNV3k\nDkknS7pf0t3u/v/Vr4uVMbPTJX1U4SI6J7r76UX2aX6QXYGyedIVrNpY8WTCaiciMnERAACgIrUG\n2QckvZLdzN/JFBbenhFpL6tgZoslzXH3vylyX0sE2SVHqgtHk7NLl/cddpj01FNjqzYWmwgJAACA\nhqlpWXV3n7SG9hQ6tFnSIkk73f2UvPaFkgYV5opvdvcvlTjERZIur1f/6q7cMuq5WtJdXWPLjQeB\ndMEFYc1qaeLS5wAAAIiVugXSk7hd0tn5DWY2TdLGbPtJkpJm9o7sfZeY2U1m9hYzmyvpd+7+SuFB\nW8JkEyOlicFzIqH0DTdUNhESAAAATTfpsup1O7FZj6T7ciPZZvZeSf3ufk52+3qFaSlfKnjcDZK+\nX2rVyZZIF6kh75myegAAAPFSU7pIE7xV4eI3Oc9IOrVwJ3e/YbIDJRIJ9fb2qre3V4lEQom4jfpW\nu4x6ufQSAAAANEQQBAqCQENDQxoaGiq7b5yC7MgE7VTXOS+9JLluQJKUWjXAxEcAAIAGKxy8NSs6\niB3eF7N0kRvcfWF2u2i6SAXHjX+6SLUoqwcAABA7Na342ACWveU8LGmemfWY2WGSLlS4+A2qTS8B\nAABAUzUlyDazuySlJZ1gZiNmdpm775d0laQfSnpU4QI4jzejfwAAAMBUNC1dpF7aMl0EAAAAsRPX\ndBEAAACgLRFkAwAAABEjyAYAAAAiRpANAAAARIwgGwAAAIgYQTYAAAAQMYJsAAAAIGIE2QAAAEDE\nCLIBAACAiBFkAwAAABEjyAYAAAAiRpANAAAARIwgGwAAAIgYQTYAAAAQMYJsAAAAIGIE2QAAAEDE\nCLIBAACAiBFkAwAAABEjyAYAAAAiRpANAEAny2QqawNQFYJsAAA6VSYjzZ8vBcFYWxCEbQTawJQc\n2uwOAACAJunqkjZtkpYulbZsCdty/+/qam7fgBZHkA0AQCdLJJTu71dywSJJUmrDWvUlEs3tE9AG\nzN2b3YdImZm323MCAKCeeo4+TiO71kiSuues1vDOXza5R0BrMDO5uxW7j5xsAAA6WRBIo6Nj26Oj\n43O0AdSEkWwAADpVduJjeuVKJdfdKklKrVquvsFBaft28rKBSZQbyW6pINvM3ippg6QXJD3p7l8q\nsg9BNgAAlcpkJgbTxdoATNBO6SKnSLrH3ZdJelezOwMAQMsrFkwTYANT1pQg28w2m9lOM3ukoH2h\nme0wsyfM7LoiD31I0ifN7EeSvt+QzgIAAABVakq6iJmdLmmPpDvc/ZRs2zRJT0g6S9KvJT0s6UJ3\n32Fml0h6j6Rdkv7F3R80sy3uvrTIsUkXAbj8WzteOwBAhWKXLuLuD0p6saD5VIV51sPu/rqkuyUt\nzu5/p7tfI+l7klaa2a2SftXIPgMtgxXcasdrBwCISJwWo3mrpKfztp9RGHgf5O6PSLpgsgMlEgn1\n9vaqt7dXiURCCYrqo5OwglvteO0AAGUEQaAgCDQ0NKShoaGy+8YpyI5MQH1PdDpWcKsdrx0AoITC\nwVuzopki4X3Nyl82sx5J9+XlZL9X0g3uvjC7fb0kL1amb5LjkpMNiBXcpoLXDgBQidjlZGdZ9pbz\nsKR5ZtZjZodJulDSd5vSM6DVsYJb7XjtAAARaFYJv7skpSWdYGYjZnaZu++XdJWkH0p6VNLd7v54\nM/oHtLRMRlq2TKmbB9TdHd5SNw9Iy5YxeW8yvHYAgIi01IqPlSBdBBBl6KaC1w4AUKG2WVa9EgTZ\nAAAAaIS45mQDAAAAbYkgGwAAAIgYQTYAAAAQMYJsAAAAIGIE2QAAAEDECLIBAACAiBFkAwAAABEj\nyAYAAAAiRpANAAAARIwgG6iXTKayNgAA0HYIsoF6yGSk+fOlIBhrC4KwjUAbAIC2d2izOwC0pa4u\nadMmaelSacuWsC33/66u5vYNAADUHUE2UC+JhNL9/UouWCRJSm1Yq75Eorl9AgAADWHu3uw+RMrM\nvN2eE1pXz9HHaWTXGklS95zVGt75yyb3CAAARMXM5O5W7D5ysoF6CQJpdHRse3R0fI42AABoW4xk\nA/WQnfiYXrlSyXW3SpJSq5arb3BQ2r6dvGwAANpAuZFsgmygXjKZicF0sTYAANCSCLIBAACAiJGT\nDQAAADQQQTYAAAAQMYJsAAAAIGIE2QAAAEDECLIBAACAiBFkAwAAABEjyAaAVpLJVNYGAGiqlgqy\nzexEM/uOmX3NzM5vdn8AoKGyK4kqCMbagiBsI9AGgFg5tNkdqNI5kr7q7g+Z2b2S/r7ZHQKAhunq\nkjZtkpYulbZsCdty/2clUQCIlaaMZJvZZjPbaWaPFLQvNLMdZvaEmV1X5KF3SrrQzL4s6Y0N6SwA\nxEkioXR/v3oWLFLPgkVK9/dLiUSzewUAKNCUZdXN7HRJeyTd4e6nZNumSXpC0lmSfi3pYUkXuvsO\nM7tE0rslrXP332T3/Xt3/1CRY7OsOoC21nP0cRrZtUaS1D1ntYZ3/rLJPQKAzhS7ZdXd/UFJLxY0\nnyrpSXcfdvfXJd0taXF2/zvd/VOSDjOzr0v6lqR1jewzAMRCEEijo2Pbo6Pjc7QBALEQp5zst0p6\nOm/7GYWB90HuPizpk5MdKJFIqLe3V729vUokEkpwKRVAO8hkpGXLlLp5QMl1A5Kk1KoBadkyaft2\n8rIBoM6CIFAQBBoaGtLQ0FDZfZuSLiJJZtYj6b68dJHzJZ3t7p/Ibl8s6VR3v7rK45IuAqB9ZTIT\ng+libQCAuotdukgJz0rqzts+NtsGAMgpFkwTYANA7DQzyLbsLedhSfPMrMfMDpN0oaTvNqVnAAAA\nwBQ0q4TfXZLSkk4wsxEzu8zd90u6StIPJT0q6W53f7wZ/QMAAACmomk52fVCTjaApiJnGgA6Rqvk\nZANAa2PZcwBAVpxK+AFAa2PZcwBAFkE2AEQpu+x5csEiSVJqw1r1UasfADoOOdkAqkfecVksew4A\nnYGcbADRIe+4PJY9BwCIkWwAtQiC4nnHnZ4Wkf0DJL1ypZLrbpUkpVYtV9/gIMueA0AbKjeSTZAN\noCbpjRuVvOp6Sdm84xUrmtyjmCCVBgA6BkE2gMiRdwwA6HTkZAOIFnnHAACUxUg2gOqQdwwAgCTS\nRQBEjbxjAAAIsgEAAICokZMNAAAANBBBNgAAABAxgmwAAAAgYgTZAAAAQMQIsgEAAICIEWQDAAAA\nESPIBgAAACJGkA0AAABEjCAbAAAAiBhBNgAAABAxgmwAAAAgYgTZAAAAQMQIsgEAAICIEWQDAAAA\nESPIBgAg7jKZytoAxEZsg2wze5uZbTKzv8trO9zMvmlmXzezi5rZPwAAGiKTkebPl4JgrC0IwjYC\nbSC2zN2b3YeyzOzv3P3D2f9fLOlFd/+emd3t7hcW2d/j/pwAAKhKEEhLl0pbtoTbuf8nEs3sFdDx\nzEzubsXuq/tItpltNrOdZvZIQftCM9thZk+Y2XUVHu5YSU9n/7+/4k5wmQ0A6o/v2vpJJJTu71fP\ngkXqWbBI6f5+Amwg5hqRLnK7pLPzG8xsmqSN2faTJCXN7B3Z+y4xs5vM7C253fMe+rTCQLuwvTQu\nswFA/fFdW3fJNes1ots0otuUXLO+2d0BMIlD630Cd3/QzHoKmk+V9KS7D0uSmd0tabGkHe5+p6Q7\nzeyNZnarpHeZ2XXu/iVJ/yhpo5l9UNJ9pc6ZSCTU29ur3t5eJRIJJTZtKn6Zrasr8ucLAB2pq0vi\nu7Z+gkAaHR3bHh0N2xjNBhoqCAIFQaChoSENDQ2V3bchOdnZIPs+dz8lu32+pLPd/RPZ7Yslneru\nV0dwrqI52emNG5W86npJUmrDWvWtWDHVUwEACtT1uzaTmRiwF2trN9mrBOmVK5Vcd6skKbVqufoG\nB6Xt29v/+QMxVi4nu2OC7J6jj9PIrjWSpO45qzW885dTPRUAoEDdvmtz6SibNo2N3gaBtGxZZwSa\nnfoHBhBzTZ34WMKzkrrzto/NttVHqctsAIDo1PO7Nj8dJQjGqm1s2tQZgWax59gJzxtoYY0Ksk3j\nJyo+LGmemfWY2WGSLpT03bqcOZORli1T6uYBdXeHt9TNA+HoB5NxACAajfiupcIGgBZS93QRM7tL\nUkLSLEk7JfW7++1mdo6kQYWB/mZ3XxvR+Sami3CZDQDqrwHftaT+AYiTcukijaguUnRlRnf/J0n/\nVO/zS+IyGwA0Qr2/a6mwAaCFxHZZdUASi1sACJH6B6DFxH5Z9WqxrHobqXc1AdKIgNbCZxZAzMSx\nuggwuXpWE2B1OqD1kPoHoIXUPScbmJJsNYHkgkWSsotbRJF/yep0AACgjkgXQezVs5oAK4ECAIBa\nkS6C1lXnhYSSa9ZrRLdpRLcpuWZ9ZMcFAACdjSAb8VXvagKsBAoAAOqEdBHEW72qCWQnPqZXrlRy\n3a2SpNSq5eobHIymcgkAAGh75dJFCLLRuSgHBgAApoAgGwAAAIgYEx8BAACABiLIBgAAACJGkA0A\nAABEjCAbAAAAiBhBNgAAABAxgmwAAAAgYgTZAAAAQMQIsgEAAICIEWQDAAAAESPIBpolk6msDQAA\ntJxDm90BoCNlMtL8+dKmTVIiEbYFgbRsmbR9u9TV1czetbdMZuLrW6wNqIPe3l4NDw83uxsAqtTT\n06OhoaGqHmPuXp/eNImZebs9J7SpIJCWLpW2bAm3c//PBd2IHn/coMnMTPyOAlpPqc9utt2KPYaR\nbKBZEgml+/uVXLBIkpTasFZ9BNj11dUVBtjF/rghwAYARIiRbKCJeo4+TiO71kiSuues1vDOXza5\nR50hvXGjklddLyn7x82KFU3uEToFI9lAa6plJJuJj0CzBIE0Ojq2PToatqHukmvWa0S3aUS3Kblm\nfbO7AwBoQwTZQDNkMtKyZUrdPKDu7vCWunkgzA2mwkh98ccNAKABYpsuYmZvk/RXkma4+4dLtRV5\nHOkiaA1UuWi87MTH9MqVSq67VZKUWrVcfYODTHxEQ7RiusjTTz+tk046SS+99JLMil4Vr9hll12m\nuXPn6gtf+EJEvUO7u+uuu3THHXfo+9//flP70VbpIu7+K3dfNlkb0LKKBXQEefXV1SVt366+FSs0\nPPyohocfDfOxCbARZ6+9Jv3qV+PbMhlpZCTS0/T29urwww/XjBkzdMQRR2jGjBn67W9/q7lz52r3\n7t1TDrCbacGCBfrGN77R7G5Ikt72trfpJz/5SbO70TIuuuiiKQfYw8PDmjZtmg4cOBBRrypT9yDb\nzDab2U4ze6SgfaGZ7TCzJ8zsunr3AwAk8ccN4ufZZ6XCX/7PPCPlRs0efFA6/XTp8cfD7UxGWrJE\nWh/tfAIz0/e+9z3t3r1bL7/8snbv3q03v/nNkZ4jrvbv39/sLlSs0YFi1OdrdP8lyd2bchWpESPZ\nt0s6O7/BzKZJ2phtP0lS0szekb3vEjO7yczektu9yDFb989pAADyfepT0ic+MRZo/+//Lf3Jn0i/\n+EW4/ed/Lq1dK73//dLPfx4G2LNmSevWhff//vfhMfLnGrz6qnTNNdLLL1fVlWJBSOEo4IIFC/S5\nz31Op59+umbMmKGFCxfqhRdeOLj/hz/8Yb3lLW/RzJkzlUgk9Nhjj1V07m9961s644wztGrVKr3x\njW/UcccdN24Ec/fu3Vq2bJmOOeYYzZ07V6tXrz7Y389//vO65JJLivb5s5/9rP7lX/5FK1as0IwZ\nM3T11VdLkqZNm6ZbbrlFJ5xwgk444QRJUjqd1qmnnqqZM2fqtNNO089+9rODx5zseed7/vnndd55\n52nmzJmaNWuWzjzzTEnSxz72MY2MjOi8887TjBkz9JWvfGXS1+yyyy7TlVdeqQ9+8IM64ogjFBSZ\nQ7JgwQKtXr1af/Znf6YjjjhCixcv1vPPP6+LL75YRx55pE477TSN5F352LFjhz7wgQ9o1qxZOvHE\nE7UlV9K0xPleeOEFnXfeeQePtXr1ap1xxhk1H69Qufc293NRybn27t2ra6+9Vr29vTrqqKP0vve9\nT3v37j34+h911FGaMWOGtm3bpv/4j//QWWedpdmzZ2vOnDm6+OKLtXv37qLvZ83cve43ST2SHsnb\nfq+kf8rbvl7SdQWPeaOkWyU9mbuvWFuRczkAAHFU9HfUyy+7n3GG+xVXuP/85+5HH+3+d383cb+v\nf2QTOUQAABvOSURBVN1dcn/Xu9xff32s/cAB9+uvd//jP3Z/7jn3V15xX7DA/WMfc//97yvuW29v\nr//4xz+e0D40NOTTpk3z/fv3u7t7IpHwefPm+VNPPeV79+71RCLhn/70pw/uf/vtt/srr7zir732\nml9zzTX+rne96+B9H//4x3316tVFz//Nb37TDzvsMN+8ebMfOHDAb731Vj/mmGMO3r9kyRJfvny5\nZzIZf+655/y0007z2267zd3db7jhBr/kkkvK9nnz5s3jzmdm/oEPfMB/97vf+d69e/2FF17wmTNn\n+re//W3fv3+/p1Ipnzlzpr/wwgsVPe98n/70p3358uW+f/9+//3vf+8PPvjguNf5Jz/5ybj9J3vN\njjrqKP/Zz37m7u779u2bcL5EIuHHH3+8/+pXv/Ldu3f7O9/5Tj/++OP9Jz/5ie/fv98/9rGP+eWX\nX+7u7q+88orPnTvXv/Wtb/mBAwf8F7/4hc+ePdsff/zxoufbu3evf+QjH/FkMul79+71xx57zOfO\nnetnnHFGTccr1v9y7+03v/nNis915ZVX+oIFC/w3v/mNHzhwwH/2s5/5a6+9dvDn4cCBAwfP+dRT\nT/mPfvQjf/311310dNTPPPNMv+aaa4q+n+4lPrtj7cXj31J3RHkrEmSfL+m2vO2LJX01onP5mWee\n6Zdeeqn39/f7Aw88UPIFAwCgkUoOBL38svuRR4a/lu++e+L9r77q/oEPhPcffbT7Y4+Nvz8XaB9/\nvPspp1QdYLuHwd8RRxzhM2fO9JkzZ/qHPvQhdy8esA4MDBx83C233OLnnHNO0WO++OKLbma+e/du\nd588yD7++OPznvKrbma+c+dO37lzp0+fPt337t178P5UKuULFixw99qD7CAIDm7feeedftppp43b\n50//9E/9W9/6VtXP+3Of+5wvWbLEn3rqqQn3lfpjJqfYa3bppZeW3D/Xty9+8YsHt6+99lo/99xz\nD27fd999/u53v9vd3b/zne/4+973vnGP/+QnP+lf+MIXip5v//79/oY3vMGffPLJg22f/exnDwa+\n1R6v0GTvbX6QXe5cBw4c8K6uLt++ffuEcxT+PBSzdetWf8973lPy/txn94EHHvD+/n6/9NJL/cwz\nzywbZLflio/FLkUAqAKVT4DG+uUvpVxe8D//c7gS6bRsRmcuB3vWLOn116VUKkwd+dGPpBNPDPcx\nk1avDtNKcsc45JCqu3HvvfdqwYIFk+6Xn6t9+OGHa8+ePZLCfNvPfOYzuueeezQ6Oiozk5lpdHRU\nRxxxRFXH7cp+3+zZs0f/f3v3Hh5Vee59/HvnIIckSMIxQAgI2uoriiJCpa9iRQ4VVKSIoIbSS3Zb\nLaJQBHRXRFsPmMLeIC99bQWFcKrKFhC6zVZBqpjigQpSxAByCkFCQiSKgZB59h+zMk6SmRBiyInf\n57pyZc2zTvfcM5N1Z82znpWbm0tRURGJif6epCVFTMeOHc/o+ZXVoUOHwPTBgwdJTk4uNT85OZms\nrKyQ8QU/77Ieeughpk2bRv/+/TEzxo4dy+TJoS8/q0zOkpKSTvtc2rRpE5hu0qRJucclse7du5eM\njAwSEhIAfy6Li4tJSUkJLB+8v5ycHIqLi0vlKnj+mW6vrL1791b6ta1oX0eOHKGwsJALLrigojQF\nHD58mPHjx/P3v/+dr7/+muLi4sB2K9K3b1/6Bt2duaILgmtrdJEsIDh7Hbw2Ealt3jB3pcaOXr/e\n36YxvEWq3yefwIABMH++vw/155+X7qO9aRO0bQsLF0JUFNx9t7+YXrjwu20cPw6DB0NKCkyaBP37\nl+6jXUnue14YtnjxYlavXs3bb79Nfn4+e/bsCf6mucqSkpJo3Lgxubm55OXlcfToUfLz89myxT+m\nQkxMDMePHw8sn52dXWr9cIVQcHu7du3Ys2dPqfn79u2jffv2ZxxvTEwMqamp7Nq1i1WrVjFz5kzW\nrVsXMpYlS5acNmfVObJLUlISffv2JS8vL5DLY8eO8dxzz4XcX6tWrYiKiuLAgQOBtv3791d5e6Hi\nqei1rWzsLVu2pEmTJuzaVf7OyaH2//DDDxMREcG2bdvIz88nLS2t2i+MrKki2yh9seIHQFczSzaz\n84A7gFU1FIuIVKRJE/jLX/xn0tav9/8MH+5v05lskeo3cybMmeP/nMXGwtq1sHcvlBQZ110HL73k\nL7BL3H03PPWUf7qoyF9gJyX5C/VnnoFBg/xnu6vpQq7KFh9ff/01jRo1Ij4+nm+++YapU6dWS4HY\ntm1b+vfvz4MPPkhBQQHOOXbv3s2GDRsA6N69Oxs2bGD//v189dVXPF1yRt/Tpk0bdu/eXeE+fvrT\nn5KZmcmyZcsoLi5m+fLlbN++nSFDhpxxvGvWrAkUe3FxcURFRRHpfbNQNpaCgoKzkrNwBg8ezOef\nf05aWhqnTp2iqKiIDz/8kB07doRcPiIigttuu43HHnuMb7/9ls8++4yFQf/gnen2yjrda1vZ2M2M\nMWPGMGHCBLKzs/H5fGRkZFBUVESrVq2IiIgoVYAXFBQQGxtLXFwcWVlZPFtyIXE1qokh/JYAG4GL\nzGyfmY1xzhUD44B0YBuwzDm3/WzHIiKV1LcvG6dNI/n6wSRfP5iN06ZB0NdjIlKNXnzRX2CXiI2F\n9HTo3r1y60dHwwMP+AvsyEh/15Enn4Tp0/3bqqSKCrvgeRUtl5KSQseOHWnfvj2XXnop11xzTaX3\nf7r9Lly4kJMnT3LJJZeQkJDA8OHDOXToEAD9+vVjxIgRXHbZZfTs2bNcYTx+/HhefvllWrRowQMP\nPBDyeSQkJPD666+TmppKy5YtSU1NZc2aNcTHx5/2eZeVmZlJv379iIuLo0+fPtx3331ce+21AEyd\nOpUnnniChIQEZs6cyejRo793zs4kttjYWNLT01m2bBnt2rWjXbt2TJkyhRMnToRdZ86cOeTn55OY\nmMjo0aMZNWoUjRo1qvL2yqrotT2T2FNTU+nWrRs9e/akRYsWTJkyBZ/PR5MmTXjkkUfo06cPCQkJ\nbNq0iWnTpvHRRx/RvHlzhgwZwrBhwyodb2XV2Ts+VpXu+ChSPZLbdGHf4ScA6Nj6d+z9svxXcCJy\nZurjHR9FypoyZQpffvklCxYsOOv7WrBgAYsXL+bNN9886/uqSIO646OI1KL160v35zxypHQfbRER\nOWfs2LGDrVu3ArBp0yZeeOEFbrvtthrZ97Zt2+jcuXON7Ku6NcjRRUTke/j2W7jnHpb+5x8Y+ewf\nAFg66Q9wzz26/biIyDmooKCAkSNHkp2dTZs2bZg0aVKV+qqfqaFDh7Jz585SN5ypT9RdRETK0xB+\nImeFuouI1E9V6S6iIltERKSGqMgWqZ/UJ1tEREREpA5QkS0iIiIiUs1UZIuIiIiIVDMV2SIiIiIi\n1UxFtoiIiEiQuLg49uzZc9rl3nnnHZKSksLOHzNmDI8++mg1Rib1icbJFhEREQlSUFBQ6WXP5Jbm\ncm7RmWwREZF6wDlHTk4ORUVFZ31fxcXFdWIbIvWZimwREZE6wOfzUVhYGHJecXExAwcOo0OHC0lK\nuqhSXRnOVOfOnZkxYwaXX345sbGx+Hy+cstMnz6d4cOHc8cdd9CsWTOuuuoqtmzZUuE2tm/fzvXX\nX098fDzdunVj9erVgeULCwuZOHEinTp1Ij4+nmuvvZYTJ04AkJGRQZ8+fYiPj+eKK67gnXfeCaz3\n4osv0qVLF5o1a0aXLl1YunQpALt27aJv3740b96c1q1bM3LkyJDPdfDgwcydO7dU2+WXX87KlSsB\niIiIYPfu3QCcPHmS3/72tyQnJ5OYmMi9994biLGszZs306NHD84//3zuuOOOUq9nbm4uQ4YMIT4+\nnhYtWnDdddeFfzGkQVCRLSIiUst27NhB27adiY09n/vum1hu/ocffsjGjf/i5MkccnJuY+7c/19u\nmcLCQmbN+g+eeuppvvrqqyrFsWzZMv72t7+Rn59PREToEmHVqlWMGDGCo0ePMnLkSG699dZSZ62D\nt+Hz+bj55psZOHAgOTk5zJ49mzvvvJPMzEwAJk6cyObNm8nIyCAvL48ZM2YQERHBwYMHGTx4MI8+\n+ihHjx4lNTWVYcOGkZuby/Hjxxk/fjxvvPEGx44dY+PGjXTv3h2A3/3udwwYMID8/HwOHDjAuHHj\nQj6H0aNHs2jRosDjTz75JLBPKN0FZPLkyezcuZMtW7awc+dOsrKyePzxx8tts6ioiKFDhzJ69Gjy\n8vIYPnw4r776amD+H//4R5KSksjNzeXw4cM8+eSTlX1ZpJ5SkS0iIlLLHn88lSNHxlJcfJj58xey\nb9++UvPbtGlDcXEusJzGjTPo1Kn8xXYpKb/i4YfTeeyxTxg48GdVimP8+PG0a9eORo0ahV2mR48e\nDB06lMjISCZMmEBhYSEZGRkht5GRkcE333zD5MmTiYqK4vrrr2fw4MEsXboU5xwLFixg9uzZtG3b\nFjOjd+/eREdHk5aWxk033cSAAQMAuOGGG7jqqqtYu3YtAJGRkWzdupXCwkLatGnDxRdfDEB0dDR7\n9+4lKyuL8847j2uuuSbkc7j55pvJzMxk165dAKSlpTFixAgiIyMBSt3Z789//jOzZs3i/PPPJyYm\nhilTpgTOnAd7//33OXXqFPfffz+RkZEMGzaMnj17BuZHR0eTnZ3NF198QWRkJH369KnUayL1l4ps\nERGRWpaY2IrzzvsYeAOzYmJjY0vN79SpE8uXz+cnP3mVBx8cwK9+9cty23j//Q8oLJzOyZOz+OST\nD6oUR4cOHQLTS5YsIS4ujmbNmnHTTTcF2oNH0zAzOnTowMGDB0Nu4+DBg+VG30hOTiYrK4sjR45Q\nWFjIBRdcUC6OvXv38te//pWEhAQSEhKIj4/nvffeIzs7m6ZNm7J8+XLmzZtHYmIiQ4YMYceOHQA8\n++yz+Hw+rr76arp168aCBQtCPs9GjRpx++23k5aWhnOOpUuXcvfdd5dbLicnh+PHj9OjR49ALIMG\nDSI3N7fcstnZ2bRv377ccy0xadIkunTpQv/+/enatSvPPPNMyNik4VCRLSIiUsumT3+EUaMSufrq\nv/Dqq4tJSEgot8yQIUN4663/4ve/fzRwxjXYvfeOoWnT24iJ+b+kpPy8SnEEd5MYNWoUBQUFHDt2\njDVr1gTa9+/fH5h2znHgwIFSxWXwNtq1a1dqeYB9+/bRvn17WrZsSePGjQNnk4MlJSWRkpJCXl4e\neXl5HD16lIKCAh566CEAbrzxRtLT0zl06BA/+MEPGDt2LACtW7fm+eefJysriz/96U/ce++9gb7V\nZaWkpJCWlsZbb71FTEwMvXr1KrdMy5Ytadq0Kdu2bQvEkp+fH7I7TmJiIllZWeWea4nY2FhSU1PZ\ntWsXq1atYubMmaxbty5kbNIwqMgWERGpZTExMcyfP5d//COdQYMGVWkbU6f+lnffXcWbby5i3rxZ\n1Rzhdz766CNee+01iouLmTVrFo0bNw5ZoAL06tWLpk2bMmPGDE6dOsX69et5/fXXGTlyJGbGL37x\nCyZMmEB2djY+n4+MjAyKioq46667WL16Nenp6YELQt955x0OHjzI4cOHWbVqFcePHyc6OprY2NjA\nPx2vvPJKoNBt3rw5ERERYfuW/+hHP8LMmDhxYsiz2OD/h2Hs2LE88MAD5OTkAJCVlUV6enrI7UVF\nRTFnzhxOnTrFihUr2LRpU2D+mjVrAv9QxMXFERUVFTY2aRj06oqIiDQQV1xxBb17967S2M2VXeeW\nW25h+fLlxMfHs3jxYlasWBEocstuIzo6mtWrV7N27VpatmzJb37zGxYtWsSFF14I+C8G7NatGz17\n9qRFixZMmTIFn89Hhw4dWLlyJU8++SStWrUiOTmZ1NRUfD4fPp+PmTNnBs6Gb9iwgXnz5gHwwQcf\n0KtXL5o1a8att97K7Nmz6dSpU9jnkpKSwqeffspdd90VNhfPPPMMXbt2pXfv3jRv3pz+/fvz+eef\nl9tWdHQ0K1asYMGCBbRo0YKXX36ZYcOGBeZnZmbSr18/4uLi6NOnD/fdd59GGGngLLhzf0NgZq6h\nPScREWkYzIz6fIyaPn06u3btYuHChbUdSrVIS0vj+eefZ8OGDbUditRx4T67XnvI/1B1JltERETO\nOcePH2fu3Ln88pflLyIVqQ4qskVEROSckp6eTuvWrUlMTAx7wxqR70vdRURERGpIfe8uInKuUncR\nEREREZE6QEW2iIiIiEg1U5EtIiIiIlLNomo7ABERkXNFcnJylcawFpHalZycfMbr1OkLH82sM/AI\n0Mw5d7vXdgtwExAHzHfO/U+ZdXTho4iIiIicdfX2wkfn3BfOuXvKtK10zv0b8Gvg9tqJrGFav359\nbYdQbyl3VafcVZ1yVzXKW9Upd1Wn3FVdfc1djRTZZvaCmX1pZlvKtA80s8/M7HMzm3yGm/13YG71\nRSn19U1cFyh3VafcVZ1yVzXKW9Upd1Wn3FVdfc1dTZ3JXgAMCG4wswjgOa/9/wAjzeyH3ry7zWym\nmSWWLF5m3aeBtc65f571yEVEREREzlCNFNnOuXeBo2WarwYynXN7nXNFwDLgFm/5Rc65CcAJM5sH\ndC85021m44AbgJ+Z2b/VRPwiIiIiImeixi58NLNkYLVz7jLv8TBggNe/GjO7C7jaOXf/99yPrnoU\nERERkRoR7sLHBjeEX7gnKiIiIiJSU2pzdJEsoGPQ4w5em4iIiIhIvVaTRbZR+gLGD4CuZpZsZucB\ndwCrajAeEREREZGzoqaG8FsCbAQuMrN9ZjbGOVcMjAPSgW3AMufc9pqIR0RERETkbKqp0UVGOefa\nOecaOec6OucWeO1/c879wDl3oXPu6aps28ymmtk2M9tiZovN7DwzizezdDPbYWZvmNn51fuM6qdQ\n45VXlCsvt5lmtt3M+tdO1HVDmNzN8HLzTzN71cyaBc1T7gg/Rr43b6KZ+cwsIahNefNUcH+BcV5+\ntnrDmZa0K3eeMJ/Xnma2ycw2e7+vCpqn3AFm1sHM3vaOqVvN7H6vXceJ0wiRu3Feu44TpxHufRc0\nv/4eK5xz9fYHSAZ2A+d5j5cDo4FngIe8tsnA07Uda134AX4MdAe2BLWFzBVwCbAZ/8WxnYCdeKPR\nnIs/YXLXD4jwpp8GnlLuTp83r70D8N/AF0CC13ax8nba91xf/N/+RXmPWyp3lc7dOqC/Nz0IWOdN\n6/P6XY7aAt296VhgB/BDHSe+V+50nKhi7rzH9fpYUadvq14Jx4CTQIyZRQFN8F88eQvwkrfMS8Ct\ntRNe3eJCj1ceLlc34+/Cc8o5twfIxD+2+TkpVO6cc28653zewwz8fwxAuQsI854DmAVMKtN2C8pb\nQJjc/Rp/gXPKW+aI167cBQmTu2yg5Axsc7670F6fV49z7pDzbvLmnPsa2I7/75qOE6cRJnftdZw4\nvXC582bX62NFvS6ynXNHgT8C+/D/wfzKOfcm0MY596W3zCGgde1FWee1DpOr9sD+oOWy+O5NL+X9\nAljrTSt3FTCzm4H9zrmtZWYpb6d3EXCtmWWY2Toz6+G1K3enNwWYaWb7gBnAVK9duQvBzDrh/zYg\ng/DHVOUuhKDc/aPMLB0nTiM4dw3hWFGvi2wzuwB4EH+3kXb4z2jfCZS9IY1uUFN5ytUZMrNHgCLn\n3NLajqWuM7MmwMPAtNqOpZ6KAuKdc72Bh4CXazme+uQFYJxzriP+48b8Wo6nzjKzWOAVYLx3ZlHH\n1EoKkbuSdh0nTiM4d0AxDeBYUa+LbOAq4D3nXJ7zj1byX8A1wJdm1gbAzNoCh2sxxrouXK6ygKSg\n5TSOeQhm9nPgp8CooGblLrwu+PvQfWJmX+DPzcdm1hqNnV8Z+4EVAM65D4BiM2uBclcZvZxzrwE4\n514Benrt+rwG8bpevgIscs6t9Jp1nKiEMLnTcaISQuSuQRwr6nuRvQPobWaNzcyAG4B/4R9v++fe\nMqOBlaFXPyeVHa88XK5WAXeYf7SWzkBXYFNNBVlHlcqdmQ3E31fsZufciaDllLvSAnlzzn3qnGvr\nnLvAOdcZOABc4Zw7jD9vI5S3Usp+Xl8DfgJgZhfhv+g7F+UulLK5yzSz6wDM7Ab8/ThBn9ey5gP/\ncs79Z1CbjhOVUy53Ok5UWqncNZhjRW1fefl9f/C/ebcBW/BfkBENJABv4i/C04HmtR1nXfgBlgAH\ngRP4+7GPAeLD5Qp/n8Wd+C9C6F/b8dfB3GUCe4GPvZ//p9ydPm9l5u/Gu2JceavUey4KWARsBT4E\nrlPuKp27Hvj7yG4G3sd/wFbuSuetD/6v6f/p5eljYGBFx1TlrsLcDdJxourvuzLL1MtjhXnBioiI\niIhINanv3UVEREREROocFdkiIiIiItVMRbaIiIiISDVTkS0iIiIiUs1UZIuIiIiIVDMV2SIiIiIi\n1UxFtohIHWBmPjNbGPQ40sxyzGxVbcZ1psxstJnNqe04RERqm4psEZG64RvgUjNr5D2+Ef9t1Ouj\nKt+Awcx0XBKRBkF/zERE6o61wE3e9EhgackMM+tpZhvN7CMze9fMLvTaLzGzf5jZx2b2TzPrYmZN\nzex1M9tsZlvMbHjZHZnZOjN72lv3MzPr47WXOhNtZqvN7FpvusDMZpjZp2aWbma9zGy9me00s8FB\nm+/obX+HmT0atK07g2KdZ2YWtN1UM9sM9K62bIqI1CIV2SIidYMDlgEjvbPZl+G/DXiJ7cCPnXM9\ngGnAU177r4D/cM5dCVwFHMB/K+ws59wVzrnLgP8Os89I51wv4EHgsTKxhBIDvOmcuxT4Gngc+Alw\nG/BE0HI9gaHA5cBwM7vSzH4IjACu8WL1AXcGbfd9L96NYfYtIlKvRNV2ACIi4uec+9TMOuE/i70G\nsKDZzYGF3hlsx3d/v98HHjGzJGCFc26nmW0FUs3sKWCNc+7dMLtc4f3+CEiuRIgnnHPp3vRWoNA5\n5/P2F7z+/zjn8gHM7FXgx0Ax0AP4wDuD3Rg45C1fHBSLiEiDoDPZIiJ1yyrgWYK6inieAN52znUD\nhuAvUnHOLfUefwusNbO+zrlM4Er8hfDvzezfw+zrhPe7mO+K9lOUPjY0DpouCpr2lazvnAsu+qH0\nmXALevyic+5K74z1xc65krPf33rbEBFpMFRki4jUDSVnrecD051z28rMPx/I8qbHBFYy6+yc+8I5\nNwdYCVxmZon4C9cl+Av2K89g/3uA7uaXBFwdYpmK1ge40cyam1kT4FbgPeBt4Gdm1sqLO97b/um2\nKyJSL6m7iIhI3eAAnHNZwHMh5s8AXvLOSq8Jar/dzO7Gf5Y5G/gD/sL4WTPzASeBX4fbX4j9v2dm\ne4Bt+PuBf1TBOuG2twl/94/2wCLn3McAXuzp3ggiJ4H78I+gorPYItLgmL6hExERERGpXuouIiIi\nIiJSzVRki4iIiIhUMxXZIiIiIiLVTEW2iIiIiEg1U5EtIiIiIlLNVGSLiIiIiFQzFdkiIiIiItXs\nfwGbv4IIch6UiAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x62e85d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot r-process yields against neutron star merger simulation yields (should be nearly identical)\n",
    "\n",
    "plt.figure(figsize=(12,8))\n",
    "plt.scatter(massnums, yields, marker='x', s=32, color='red', label='Final neutron star merger ejecta')\n",
    "plt.scatter(rmassnums, massfracs, s=8, label='r-process yields')\n",
    "plt.xlim(80, 250)\n",
    "plt.ylim(0.000000000001, 1)\n",
    "plt.yscale('log')\n",
    "plt.xlabel('Mass number')\n",
    "plt.ylabel('Mass fraction')\n",
    "plt.legend(loc=4)\n",
    "#plt.savefig('yield_comparison.png', dpi=200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check DTD fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Define the three functions which are fit to the DTD (Power law, 5th and 6th degree polynomials)\n",
    "\n",
    "def quintic(t, a, b, c, d, e, f): \n",
    "    y = (a*(t**5))+(b*(t**4))+(c*(t**3))+(d*(t**2))+(e*t)+f\n",
    "    return y\n",
    "\n",
    "def sextic(t, a, b, c, d, e, f, g):\n",
    "    y = a*(t**6) + b*(t**5) + c*(t**4) + d*(t**3) + e*(t**2) + f*t + g\n",
    "    return y\n",
    "\n",
    "def powlaw(a, t):\n",
    "    y = a / t\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x600e450>]"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeYlNX1wPHvoYiIURFBUBTELmhAFDRYRoRILGDDggVb\nokawxMSaBDSJPyUJiSVYwSDRIKEKqCDCmIAiKCBVQDoIiIoKCki5vz/OrK6wCzOz78x9y/k8zzzu\nzu7MHJd9z945995zxTmHMcaY+KjkOwBjjDHBssRujDExY4ndGGNixhK7McbEjCV2Y4yJGUvsxhgT\nM1kndhGpJiLvishUEZklIg9l7u8mIstFZErm1q5w4RpjjNkVyWUdu4js4Zz7RkQqAxOAO4E2wDrn\nXM8CxWiMMSYHOZVinHPfZD6slnns2sznEmRQxhhj8pdTYheRSiIyFVgFpJ1zszNf6iIi00TkORHZ\nO/AojTHGZC2nUsx3DxLZCxgN3A3MBj51zjkR+SNQzzl3fbBhGmOMyVaVfB7knPtKREYCJzjn3ir1\npWeB4WU9RkSsKY0xxuTBOZdTuTuXVTH7lZRZRKQ60BaYJiJ1S33bhcDMnQQX2Vu3bt28x5DU+KMc\nu8Xv/xb1+PORy4i9HtBXRAT9g9DPOfemiLwgIk2BbcBi4Ma8IjHGGBOIrBO7c24GcHwZ918daETG\nGGMqxHaeZimVSvkOoUKiHH+UYweL37eox5+PvFbF5PVCIq5Yr2WMMXEhIrgcJ0/zWhVjgucczJ4N\ns2bp5yJQqxacfDJUr+43NmNMtFhi9+ztt+HJJ2HMGNh9dzj+eKiUKZCtWAEzZkDLlnDRRXDddVCt\nmt94jTHhZ6UYT2bNgvvug6lT4a674OyzoVGjHb/vyy8hnYZnn9Uk3707XHUVVLE/ycYkQj6lGEvs\nRbZtG/zpT/D443DPPfDLX+pIPRsTJugfg/Xr4T//KfsPgTEmXiyxh9y6ddC5M6xaBYMGQb16uT+H\nc/DYY/DQQ9CnD5xzTvBxGmPCI5/Ebssdi2TZMp0IrV0bxo3LL6mDTqredhsMHgw33qgJ3hhjSrNK\nbRGsWQNt2+rk5113BfOcrVrBe+9B69Za3vntb4N5XmNM9FliL7CvvoJ27eDii4NL6iXq1oWxYyGV\n0snUe+4J9vmNMdFkib2ANm6E9u3hpJPgD38ozGuUTu7Vq2uZxhiTbDZ5WkC//CWsXq0rWCoVeDZj\n6VJd7/7CC1r2McbEg02ehsiAATBqlK5cKXRSBzj4YOjfH668EhYuLPzrGWPCy0bsBfDRR7oC5vXX\noXnz4r7244/rZqZ33oEaNYr72saY4Nk69hDYtAl+8hO45hro2rX4r++crr7ZulXLMsaYaLNSTAj0\n6KFr1Lt08fP6IvDEE9qDZsgQPzEYY/yyEXuA5s3T0fqUKVrz9mnCBF1i+cEHUKeO31iMMfmzUoxH\nzkGbNnDuuXDHHb6jUXffrfX+gQN1JG+MiR4rxXjUrx988YWfunp5HnxQ30W89JLvSIwxxZT1iF1E\nqgH/BXbL3IY55+4TkZrAy0AD9DDrS5xzX5bx+NiO2NeuhaOPhpEji78KZlcmT9ZNUnPmwD77+I7G\nGJOrgpdiRGQP59w3IlIZmADcCbQHPnPO9RCRu4GazrkdNrfHObHfdZeO1p95xnckZbvxRm0N/Oij\nviMxxuSqaDV2EdkDSAPXAIOB051zq0WkLpB2zh1VxmNimdiXLYOmTfUQjAMO8B1N2T77DI45Bt54\nA447znc0xphcFLzGLiKVRGQqsApN4LOB/Z1zqwGcc6uARK3B6NYNbropvEkd9OzUBx+EW27RSV5j\nTLzllNidc9ucc82A+sCpIpICtk8ViUkdM2dqXT3oro2FcMMNsGGDTaQakwR5dXd0zn0lIq8CJwCr\nRWT/UqWYT8p7XPfu3b/7OJVKkUql8nn50Lj3Xr3tvbfvSHatcmWtsV9xha5vt0OxjQmndDpNOp2u\n0HPksipmP2Czc+5LEakOjAIeAH4KfO6ceyRJk6fvvguXXKLLCaOUJNu318M5br/ddyTGmGwUdPJU\nRI4F+gKClnD6Oef+IiL7AgOAg4Al6HLHL8p4fKwSe4cO8NOfat06SmbOhDPP1D9IUXinYUzS2c7T\nIpkxQ5P6woV6uEXUXHst1K9fuMM/jDHBscReJJ066RLHKEyalmXpUmjWDGbN0hOYjDHhZYm9CEp6\nrS9YAHvt5Tua/P3qV9ra1zYtGRNultiL4Be/0La8DzzgO5KKWbUKGje2UbsxYWeJvcBWrtRkOH++\nbvqJuttug6pV4S9/8R2JMaY8ltgLrFs3WLMGevXyHUkwVqzQFgMffgi1a/uOxhhTFkvsBbRpEzRo\nAOPGaSfHuLj5ZqhZEx56yHckxpiyWGIvoBdegBdfhFGjfEcSrCVL4Pjjtby0776+ozHGbM8O2igQ\n5+Cxx+DWW31HErwGDeD88+Hxx31HYowJio3YszBhAnTurLs1K8XwT+GcOZBKweLF0dxwZUyc2Yi9\nQB57TI+8i2NSB50zaNFCj/czxkSfjdh3YeVKPaRiyZJob0jalbfe0jX6c+bE9w+YMVFkI/YCeP55\n6Ngx3kkd4LTT9P9xxAjfkZgkW7gQ/vlPeP99+PZb39FElyX2ndi2DZ57Tg+piDsR+PWvbbOSKT7n\n4M03taV0ixbw6qtwzTV6+Po558Dnn/uOMHosse/EuHHwox/BiSf6jqQ4LrpIG4S9+67vSExSOAf3\n368Hrp93nv7+DRigHVTXrIGjjtJ3kytW+I40WvI6QSkpnn0Wfv5zHc0mQZUq2mbg0UftCD1TeM7B\n3XfD6NEwcSLst98Pv16jhr6DrFMHTjlFv+/ww/3EGjU2eVqOTz+Fww6DRYt0Z2ZSfPEFHHIIzJ6t\nzc6MKQTn4M47ddL+jTd2vTnumWfgkUe0ad3uuxcnxrCwydMA9eunNb8kJXXQuuall+qFZEyh9OkD\nY8boLZsdz7/4Bfz4x9CzZ+FjiwMbsZfBOe3i+PTTcOqpvqMpvpkz9YSoxYtht918R2PiZvFinbca\nNw6aNMn+cQsX6uTqBx/AgQcWLLzQsRF7QCZPhs2bta6XRE2a6KTVoEG+IzFxs22brnj5zW9yS+oA\njRrpJOs99xQktFjJOrGLSH0RGSsis0Rkhoh0zdzfTUSWi8iUzK1d4cItjn794OqrkzNpWpauXeGJ\nJ3xHYeLmscd00HTnnfk9/t57daT/zjvBxhU3WZdiRKQuUNc5N01E9gTeBzoAlwLrnHM7rX5FpRTz\n7bf6Nm/SJJ1ETKotW+DQQ2HIEO3+aExFLV4MJ5ygK2AOOyz/53nhBejdWydek6CgpRjn3Crn3LTM\nx+uBOUBJpSs2Y9vXXtPeKUlO6qBLH2+6CZ580nckJi66dYNf/rJiSR3g8sv1cJh584KJK47yqrGL\nSEOgKVCylaWLiEwTkedEZO+AYvOipAxj4NprYeBA+Oor35GYqJs5E15/XXc3V1TVqnDVVdruw5Qt\n51UxmTJMGviDc26YiNQGPnXOORH5I1DPOXd9GY9z3bp1++7zVCpFKpWqSOyBW7sWGjbUhl/77OM7\nmnDo2BFat9aTlozJV/v2cMYZcMcdwTzf7NnQpo3uVK0Ss22W6XSadDr93ecPPPBAYU9QEpEqwAjg\nNefco2V8vQEw3Dl3XBlfC32N/emnYexYePll35GEx5gxOtE1bVqyJ5NN/iZMgE6dYO7cYDcXnXwy\n/Pa32k8mzoqx3LEPMLt0Us9Mqpa4EJiZ43OGxgsvWBlme61bw9dfW/8Ykx/ndCXLAw8Ev2P0uut0\no5PZUS6rYloB/wVmAC5zuw/ohNbbtwGLgRudc6vLeHyoR+yLFkHLltpsqGpV39GES48e2qfdapom\nV+m0rj2fPRsqVw72ub/6Cg4+WM/rrV072OcOEzvMugIeeUSXY9kqkB198gkccYT+fGzuweTivPPg\n3HM1uRdC587QrBncfnthnj8MbOdpBfTvrz1SzI7q1IF27eBf//IdiYmSDz/U/SCFLG9efbV1Ii2L\nJXb0F3D16mT2hcnWDTdYPdPk5m9/070QhTwg/dRT9fr99NPCvUYUWWJHV8FccknwNcA4ad1aT7KZ\nOtV3JCYK1qzRAzNuuaWwr7PbbnD66bp6y3wv8YndOS3DXHaZ70jCrVIl3bBko3aTjV69dA9EnTqF\nf62zzoJRowr/OlGS+MnT6dN188SiRbZOe1eWLIHmzWH58uQddmCyt3EjNGigK2KOPrrwrzd/PqRS\n+nsZx2vYJk/zUDJpGsdfiKA1aKANwYYM8R2JCbPBg/VQjGIkddDeM9WqadsCoxKd2J3T+rqthsne\n9ddbOcbs3LPP6olHxSKiq7asHPO9RCf2qVP1l6JZM9+RREeHDvpzW7zYdyQmjObN081s7dsX93Wt\nzv5DiU7sgwbBRRdZGSYXu++ubVP79vUdiQmjZ5/VE5KKfaTiGWdon/evvy7u64ZVYhO7c9qS9uKL\nfUcSPZ07a1+dEM6FG482bdI/+DfcUPzX3msvnf9JyuEbu5LYxD5rls7en3CC70iip3lz3XQyfrzv\nSEyYDB0Kxx1X8YM08mXlmO8lNrEPHGhlmHyJ6KjdyjGmtGeeKe6k6fbOOMMGGyUSu4792GPhqaeg\nVSvfkUTTypXQuLGuHd5jD9/RGN8WLtTuqMuX69JDHzZsgFq14LPPCtvGoNhsHXuW5s7Vf/yTT/Yd\nSXTVq6cXsq1pN6AN4i67zF9SB03mxxwDU6b4iyEsEpnYBw2CCy/UbfImf1aOMaCT6P366TmkvrVs\nqR0lky6Rqa1kmaOpmA4d4P339e23Sa6JE7WB3okn+o4EWrSw074ggYl9yRK9WYveiqteXf9Avvii\n70iMTyWj9TAsRGjZ0hI7JDCxDxump7rE7WRzX666Si/sEM2LmyLatEnbclxxhe9I1BFHwNq1eupX\nkiUusQ8dChdc4DuK+GjVCtavhw8+8B2J8WHkSF1h1rCh70hUpUpaEkp6nT3rxC4i9UVkrIjMEpEZ\nInJr5v6aIjJaROaKyCgR2btw4VbMZ59pTbhtW9+RxEelSnDllTpqN8kTlknT0qwck9uIfQvwK+dc\nY+Bk4BYROQq4BxjjnDsSGAvcG3yYwRg5Es48M15rXMPgyivh3/+GrVt9R2KK6fPPYezY8LXlaNHC\nRuxZJ3bn3Crn3LTMx+uBOUB9oANQsuitL3B+0EEGZehQOD+00UXXUUfBgQfCm2/6jsQU05Ah0KYN\n7B2y9+gliX3bNt+R+JNXjV1EGgJNgYnA/s651aDJHyjCYVi5++YbTTznnOM7kngqmUQ1yRHWIyXr\n1tWmYB995DsSf3JeGyIiewIDgducc+tFZPv1EOWuj+jevft3H6dSKVKpVK4vn7cxY7R5Va1aRXvJ\nRLnsMvj973Uidc89fUdjCm31apg8WVeZhVFJnf2II3xHkrt0Ok06na7Qc+TUK0ZEqgAjgNecc49m\n7psDpJxzq0WkLjDOObfDoVi+e8Vcdx00bQq33uothNg75xzo1Ck8S99M4fTqBRMmhHcPw1//qucY\nP/GE70gqrhi9YvoAs0uSesYrwDWZjzsDofsbvnUrjBihOyVN4VxxRXgvdBOssJZhSjRrluwluFmP\n2EWkFfBfYAZabnHAfcAkYABwELAEuMQ590UZj/c2Yp8wAW65BaZN8/LyibF+PdSvr8ej1QnlTIsJ\nwvLlelj1xx/7bfq1M6tX62Han30Wjh2xFZHPiD3rGrtzbgJQuZwvt8nlRYvtlVd0t6kprD331HLM\ngAHQpYvvaEyhDBigq8vCmtRBBxaVK2t76QMO8B1N8SVi5+nw4cU/XDeprrgCXnrJdxSmkMJehgEd\npTdpoielJVHsE/v8+do7onlz35EkQ9u2usxs4ULfkZhCWLRIb2ec4TuSXWvcGGbO9B2FH7FP7MOH\naxnGeq8XR9Wq0LGjjdrjatAg7bUUhSZ6NmKPMSvDFF/J6hjr+Bg/AweGr4VAeZI8Yo/1madr10KD\nBjpDbv1hisc5OPRQTQLHH+87GhOUpUv133PlSn1nFnaff65dJ7/8MtorY+zM0+289prWAi2pF5cI\nXH65NgYz8TFokO4FiUJSB9h3X12ptXSp70iKL9aJ3ZY5+tOpk66eSHIjpriJUhmmRFLr7LFN7Js3\nw6hR1vTLl8aNoWZNGD/edyQmCMuXw5w52vY6Spo0SWadPbaJffx4OPxwqFfPdyTJ1amTlWPiYvBg\nXYSw226+I8lN48Y2Yo+VESPg3HN9R5Fsl12mb9+//dZ3JKaioliGARuxx86IEVaG8a1hQ22b+sYb\nviMxFbFqFUyfHs0jJY85Bj78MHmne8Uysc+bB+vWaYc345eVY6Jv6FA4++xw94Ypz49+BLVr627Z\nJIllYh85UkfrttvUv44d9d3TN9/4jsTka/BguPBC31HkL4kblWKZ+qy+Hh516sBJJ+kOYBM9n38O\nEyfCz37mO5L8JXHJY+wS+5df6kG2UVuWFWe2WSm6hg/Xa6lGDd+R5O+YY2D2bN9RFFfsEvvo0XDK\nKXbuZphccAGMG6ctHky0RL0MA7rsOWkHW8cusVsZJnz22gvatNEkYaJj/Xr9gxz13duHHgoLFviO\norhildi3btX+MLbMMXysHBM9r70GP/kJ7LOP70gqpk4d2LQpWe8YY5XYJ0/Wf8SGDX1HYrZ3zjnw\n/vu6JtpEw6BB0S/DgDalO+ywZI3as07sItJbRFaLyPRS93UTkeUiMiVza1eYMLNTsszRhE/16vqW\nfsAA35GYbGzcCK+/rt0c4yBp5ZhcRuzPA2eVcX9P59zxmdvrAcWVF0vs4WblmOh480047jjYf3/f\nkQTjsMOSNYGadWJ3zo0HyqpShaKF/YoVsHix1gRNOLVpo6OmpO0CjKIhQ3Q1U1xYKSZ3XURkmog8\nJyJ7B/B8eXn1VTjrrGicxZhUVavCRRdpn3YTXlu2wLBh8Urshx6arBF7RdNgL+BB55wTkT8CPYHr\ny/vm7t27f/dxKpUilUpV8OW/N3JkNLvPJc3ll0OXLnDvvb4jMeWZMAHq14/XIoQolWLS6TTpdLpC\nz5HTmaci0gAY7pw7LpevZb5esDNPN23S1TALFsB++xXkJUxAtm3Tc2hff117eJjwuf12qFULfvc7\n35EEZ9s23T376afR20VbjDNPhVI1dRGpW+prFwJeWu289Zb2g7CkHn6VKmmfdptEDSfn4ldfB/29\nO+QQWLjQdyTFkctyx5eAt4EjRGSpiFwL9BCR6SIyDTgduKNAce6UrYaJlpLVMQV6A2cqYMoUPSUp\nju+molSOqaisa+zOuU5l3P18gLHkxTlN7IMG+Y7EZKtZM53knjQJWrb0HY0prWS0LqFY6xasJK2M\nifzO07lztcZ+XJmVfRNGIramPaziWIYpkaSVMZFP7CNH6ukucRxhxNnll+su1KQdWRZmc+dqP5W4\nvotKUikmFondujlGz5FHQr16UMFVXSZAJaP1uJ48ZqWYiPjyS3jvPWjd2nckJh92Hmq4DB4c3zIM\nwMEHw8cfa+k27iKd2EePhlatorcu1ahLL9VRYhIutLBbtkxHs6ef7juSwqlaFQ46SFuPxF2kE7st\nc4y2+vXh2GN1s5Lxa+hQLWlWreo7ksJKSjkmsol92zbtD2OJPdo6dYKXXvIdhRkyJB6913clKStj\nIpvYJ0+G2rV1N5mJrosu0hH7unW+I0muTz/VQ1B++lPfkRSeJfaQszJMPNSqBaedpqUA48fw4dC2\nrR6GEncNG8LSpb6jKDxL7MY7K8f4FffVMKUdfHAyEntO3R0r9EIBdndcsUIn3T75xPqvx8HXX8OB\nB8K8edql0xTPunX6s1+6NPqHVmfjk0/gmGO0/BQVxejuGAojR0K7dpbU46JGDX33ZeehFt/IkXDK\nKclI6qDzcl9/Dd984zuSwopkYh8xwnabxs0VV8CLL/qOInkGD07GapgSIrrMdtky35EUVuQS+4YN\nug29XTvfkZggtW2r64uTsMY4LDZsgFGjoEMH35EUVxLq7JFL7OPGadvXfff1HYkJUtWqcMklNola\nTKNHQ/PmWp5IEkvsITR8uJVh4qqkHGMHcBTHoEHJKsOUsMQeMs5ZfT3OTjoJNm/WzTKmsL79Vq+l\npCxzLO2gg6zGHirTp0O1anDUUb4jMYUgYpOoxTJunF5HBx7oO5LisxF7yJSM1u1Qjfi64gro398O\n4Ci0pJZhwBL7D4hIbxFZLSLTS91XU0RGi8hcERklInsXJkxl9fX4O/JIXY725pu+I4mvLVtg2LDk\nJvaSUkyc53JyGbE/D5y13X33AGOcc0cCY4F7gwpseytX6tFdce4XbdRVV0G/fr6jiK///leTW6NG\nviPxo0YNvUVp92musk7szrnxwNrt7u4A9M183Bc4P6C4djBihK5dj3u/aAOXXabvztav9x1JPA0c\nCBdf7DsKvw46KN7lmIrW2Os451YDOOdWAQXr9PHKK9C+faGe3YRJnTq6zX3IEN+RxM/WrbrbtGNH\n35H4Ffc6e9DdVnZaterevft3H6dSKVKpVFZP+vXX8NZb9vY8Sa66Cnr31v+a4Pzvf3DAAdqXPMnC\nnNjT6TTpCp7ynlN3RxFpAAx3zh2X+XwOkHLOrRaRusA459zR5Tw27+6OQ4fCE0/AmDF5PdxE0IYN\nuhRv5kxNRCYYXbroz/O++3xH4tef/wyrV8Nf/uI7kl0rRndHydxKvAJck/m4MzAsx+fLipVhkqd6\ndd08Yy0GgrN1qy5zTHoZBqzG/h0ReQl4GzhCRJaKyLXAw0BbEZkLnJn5PFBbt+rEqSX25LHVMcGa\nMAH23x8OP9x3JP6FuRQThKxr7M65TuV8qU1AsZRp4kSoV0+PtDLJctpp8OWXMG0aNG3qO5ros9Uw\n34t7Yg/9zlMrwyRXpUpw9dXQt++uv9fs3NatmtitDKPq1dN17N9+6zuSwgh1YndOl7wlrV+0+V7n\nzlpnj+sFWCz/+5+WYY480nck4VC5sk4ir1jhO5LCCHVinz0bNm7UntEmmQ49VJPRa6/5jiTa+vfX\njV/me3GeQA11Yi85tsuafiVb587wz3/6jiK6Nm/Wa+nSS31HEi5xrrOHOrEPGZLcRkXmex07apvZ\nNWt8RxJNY8fqOx9bgPBDltg9WLRI61+tWvmOxPi2115w3nnw73/7jiSa+ve30XpZmjeHmjV9R1EY\nOe08rdAL5bjztGdPmDMHnn22gEGZyHjzTfj1r2HqVN+RRMumTboCZMaMZB6qEQfF2HlaNCX1dWMA\nzjgDPv/cEnuuRo2C446zpJ40oUzsq1bBrFnQurXvSExYVKoE112njcFM9qwMk0yhLMU8/bQeBmBn\nX5rSli6FZs1g+XLtJWN2bt06XdI3fz7Uru07GpOv2JRiBg5M5unpZucOPhhOOEHLdGbXhg6FU0+1\npJ5EoUvsn3wCkyfD2Wf7jsSE0Q03WDkmW//6F1x5pe8ojA+hS+yDBmlS32MP35GYMGrfXnu0L1jg\nO5JwW7kSJk3SZaImeUKX2F9+2SZ7TPmqVdNRaJ8+viMJt/794fzzbYCUVKGaPP34Y2jSREcb1aoV\nJSwTQbNmQdu2OplaJejDHWOieXPo0QPOPNN3JKaiIj95OnCgvnW0pG52pnFjOOQQPYDF7Gj2bF0y\nnOWRwiaGQpXYrQxjsnXzzfDUU76jCKcXX4ROnbQ1rUmm0JRili3TNcoffwy77VaUkEyEbdyoa7Tf\nfRcaNfIdTXhs3arvZoYPhx//2Hc0JgiRLsUMGKCTPZbUTTZ2313b+T79tO9IwmXsWF23bkk92QJJ\n7CKyWEQ+EJGpIjIpn+d48UU7CMDk5he/gOef10ZXRvXpA9de6zsK41tQI/ZtQMo518w51yLXB8+Y\nob22zzgjoGhMIhxxhDa4sp2oau1aPWmqU3nHzpvECCqxS0Weq29fuOoqm+wxubv5ZnjySd9RhMO/\n/w1nnQX77us7EuNbUIndAW+IyGQR+XkuD9yyRcswnTsHFIlJlPbtdRfq9Om+I/Hv+eetDGNUUNs7\nWjnnVopIbTTBz3HOjd/+m7p37/7dx6lUilQqxejR0KCBnZ5u8lO1qo7aH3882YeyzJihG/vatvUd\niamodDpNOp2u0HMEvtxRRLoB65xzPbe7v8zljpdeqhspbr450DBMgnzyiQ4MPvoIatXyHY0fv/qV\ntjL+0598R2KCls9yxwondhHZA6jknFsvIjWA0cADzrnR233fDon9iy/0gN2FC60uaCrmmmvg6KPh\n7rt9R1J8GzZoS+OJE/XQahMvvtax7w+MF5GpwERg+PZJvTwDBuhbR0vqpqK6doV//EPnbJJm4EDt\nDWNJ3ZSocGJ3zi1yzjXNLHU81jn3cHaPg2eesckeE4zmzXXUOmyY70iK78kn4aabfEdhwsTbztPJ\nk/Vw4rPO8hWBiZtbb4VHH/UdRXF98IF2uTz3XN+RmDDxltj/8Q+dMLW16yYoF1wAS5booCEpnnoK\nfv5za19sfshLE7BPP4XDD0/2KgZTGH/7G7zzjs7fxN26dVp+mjkTDjzQdzSmUCLTBKx3b234ZUnd\nBO2GG7QRVhKOznvxRWjd2pK62VHRR+xbt+rs/cCBeuK8MUG7/35dSvuPf/iOpHCc09PGHn0U2rTx\nHY0ppEiM2F99Ffbf35K6KZyuXbVvypo1viMpnNGjdX7Kjr4zZSl6Yv/rX/XCM6ZQ6taFiy+O94i9\nZ0/dbSo5jeNMUhS1FJNOO66/Hj780GbxTWHNnQunnqq19h/9yHc0wZo5Uzf2LV5s5wMnQehLMX/4\ng9Y/LambQjvySC1T9OrlO5Lg/f3vcMstltRN+Yo6Ym/Y0DFvnnbkM6bQZs/Ww1sWLIA99/QdTTBW\nr4ajjoL582G//XxHY4oh9CP2++6zpG6K55hjtHNonEbtvXppR1RL6mZnijpi37TJ2WHVpqhmztTl\ngAsWQI0avqOpmK++0qXCb7+tG/xMMoR+xG5J3RRbkyZw2mnxOD7viSe0t5IldbMrXloKGFNMJaP2\nefNgr70rRiW7AAAIRklEQVR8R5Of9et1tJ5Oa995kxyhH7Eb40OTJjrS/fOffUeSvyef1IlgS+om\nGzZiN4mwdCk0a6Zngx5wgO9ocvPNN9CoEYwZo3+kTLLYiN2Ychx8MFx/PZQ6Tz0ynnoKWrWypG6y\nZyN2kxhr1+rGpXRal0JGwRdfaMxjxsCxx/qOxvjgbcQuIu1E5EMRmSciCTxO2ERBzZpwzz3ROvD6\n//4PzjvPkrrJTYVH7CJSCZgHnAl8DEwGLnPOfbjd99mI3Xi3aZMmyZ49w3+c3JIlcPzx0ZwXMMHx\nNWJvAcx3zi1xzm0G+gMdAnheYwJXrZp2fezaVSclw+z++6FLF0vqJndBJPYDgWWlPl+euc+YUGrb\nFlq2hIce8h1J+d57T0+C+s1vfEdioshWxZhE6tlTV5t8+OGuv7fYtm7VdxQPPhif5mWmuIJooLsC\nOLjU5/Uz9+2ge6m1ZqlUilQqFcDLG5O7Aw6A3/4Wbr4Z3nwTKoVoiNOrlzbLu+4635EYH9LpNOl0\nukLPEcTkaWVgLjp5uhKYBFzunJuz3ffZ5KkJlS1btI9Mx45wxx2+o1HLlulGqvHjtT2vMflMngay\njl1E2gGPoqWd3s65h8v4HkvsJnQWLICTToJx4/xvAHIO2reHFi3gd7/zG4sJD2+JPasXssRuQqp3\nb3jsMZg0ye+pRC+/rKeMTZlinVDN9yyxG5MH5+DCC+Gww/w1CluwAE4+GUaOhBNP9BODCSfrFWNM\nHkTgmWfgP//RUXOxbdyodf7f/96SugmGjdiNyfjgA+3b/tprcMIJxXvdm2+Gzz7TPyqS07jMJIGN\n2I2pgB//WEfuF1wAK1cW5zVfeEEbfD33nCV1E5wg1rEbExsXXACzZmkfmTFjtHFYoQwfDnfdpevo\no3qykwknG7Ebs53774fTT4czz9QSSSGMG6f94V95BRo3LsxrmOSyxG7MdkTgr3/VnjKtW8OaNcE+\n/9tvw6WXwoABumbdmKBZYjemDCLw8MPaC71VK5g2LZjn7dMHzj9fa+vWUcMUiq2KMWYXXnoJbrtN\nm3LddFN+k5ybN2vbgjfegGHDrF2AyZ6tijGmADp10vLJM8/opOrUqdk/1jmtozdtqgdqT5pkSd0U\nno3YjcnSxo3a6rdHD91IdPvtult09913/N4vvtDVLn//u561+vDDcM45tqTR5M5aChhTBBs26Lrz\nvn21n/vxx8PRR2sf9S1b4KOPdLPTKafAZZfBlVdC5cq+ozZRZYndmCJbtw7eeUd7vVSpon3UDzgA\nTj0Vqlf3HZ2JA0vsxhgTMzZ5aowxxhK7McbEjSV2Y4yJGUvsxhgTM5bYjTEmZiqU2EWkm4gsF5Ep\nmVu7oAIzxhiTnyBG7D2dc8dnbq8H8HyhlE6nfYdQIVGOP8qxg8XvW9Tjz0cQiT0Rm6Sj/ssR5fij\nHDtY/L5FPf58BJHYu4jINBF5TkT2DuD5jDHGVMAuE7uIvCEi00vdZmT+ex7QC2jknGsKrAJ6Fjpg\nY4wxOxdYSwERaQAMd84dV87XrZ+AMcbkIdeWAhU6zFpE6jrnVmU+vRCYGVRgxhhj8lOhxA70EJGm\nwDZgMXBjhSMyxhhTIUXr7miMMaY4CrLzVER6i8hqEZle6r6aIjJaROaKyKgwr6ApJ/4eIjInswJo\nkIjs5TPG8pQVe6mv3Ski20RkXx+xZaO8+EWka+bnP0NEHvYV366U87tzoohMEpGpmf+e4DPG8ohI\nfREZKyKzMj/nWzP3R+LaLSP+rpn7o3LtlvnzL/X17K9f51zgN+AUoCkwvdR9jwB3ZT6+G3i4EK9d\nwPjbAJUyHz8M/J/vOLONPXN/feB1YBGwr+84c/zZp4DRQJXM5/v5jjPH+McBP818/DNgnO84y4m9\nLtA08/GewFzgqKhcuzuJPyrXbpnxZz7P6fotyIjdOTceWLvd3R2AvpmP+wLnF+K1g1BW/M65Mc65\nbZlPJ6I/6NAp52cP8DfgN0UOJ2flxH8zmky2ZL7n06IHlqVy4l8JlIxy9wFWFDWoLDnnVjnnpmU+\nXg/MQX/PI3HtlhP/gRG6dsuMP/PlnK7fYjYBq+OcWw36PwDUKeJrB+064DXfQWRLRNoDy5xzM3zH\nkqcjgNNEZKKIjAtrKWMn7gF6ishSoAdwr+d4dklEGqLvPCYC+0ft2i0V/7vbfSkS127p+PO5fiu6\nKqYiIjlrKyL3A5udcy/5jiUbIlIduA9oW/puT+HkqwpQ0zl3koicCAwAGnmOKRe9ga7OuaEicjHQ\nhx/+e4SKiOwJDARuc86tL2MPSqiv3e3jL3V/JK7d0vEDW8nj+i3miH21iOwPuv4d+KSIrx0IEbkG\nOBvo5DmUXBwKNAQ+EJFF6NvQ90Uk9KOuUpYBgwGcc5OBbSJSy29IOWnpnBsK4JwbCLTwHE+5RKQK\nmlT6OeeGZe6OzLVbTvyRuXbLiD+v67eQiV344V+WV4BrMh93BoZt/4CQ+UH8mZbEvwHaO+c2eYsq\nO9/F7pyb6Zyr65xr5Jw7BFgONHPOhfbiZMffnaFAawAROQKo6pz7zEdgWdo+/vkicjqAiJwJzPMS\nVXb6ALOdc4+Wui9K1+4O8Ufs2v1B/HlfvwWa3X0J+BjYBCwFrgVqAmPQmd7RwD6+Z6FzjH8+sASY\nkrn18h1ntrFv9/WFhHtVTFk/+ypAP2AG8B5wuu84c4y/OVrrnQq8k7kwvcdaRuyt0Lf+0zKxTgHa\nAftG4dotJ/6fRejaLfPnv933ZHX92gYlY4yJGTsazxhjYsYSuzHGxIwldmOMiRlL7MYYEzOW2I0x\nJmYssRtjTMxYYjfGmJixxG6MMTHz/x3SaDUbuB44AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x60d6610>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solar metallicity fit, parameters from chem_evol (see fitting notebook to derive new parameters)\n",
    "a = -0.0138858377011\n",
    "b = 1.0712569392\n",
    "c = -32.1555682584\n",
    "d = 468.236521089\n",
    "e = -3300.97955814\n",
    "f = 9019.62468302\n",
    "t = np.linspace(10, 22.2987, 100) # Polynomial portion of solar metallicity DTD x-axis\n",
    "\n",
    "# Define the DTD fit and plot it\n",
    "y = quintic(t, a, b, c, d, e, f)\n",
    "plt.plot(t, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x44e6110>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEACAYAAACj0I2EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHMRJREFUeJzt3Xu81XO+x/HXJ2WYSBcRZZCZyK3ILZdpRSLGXZhChDGM\nxrgfOWPvOXOGOMalMHQmuRUnoYtbiRbdpYsuE8KMTLQHQ9Gk9m5/zx/flcm2L2uv23f9fuv9fDz2\nw9qrtdfvXfp9+q7P7/v9fc05h4iIRF+T0AFERCQ3VNBFRGJCBV1EJCZU0EVEYkIFXUQkJlTQRURi\nosGCbmYjzKzCzBbVeH6QmS0zs8VmNiR/EUVEJB1N03jNSGAY8OimJ8wsAZwE7OecqzKz7fMTT0RE\n0tXgCN05Nx34osbTlwFDnHNVqdd8lodsIiLSCJn20DsBPzWz2WY21cwOymUoERFpvHRaLnX9XCvn\n3GFmdjAwBuiYu1giItJYmRb0j4BnAJxzc82s2szaOOc+r/lCM9PNYkREMuCcs8a8Pt2Wi6W+NhkH\nHA1gZp2AZrUV881CRfarrKwseIZSzR/l7Mof/ivq+TPR4AjdzEYDCaCNma0AyoCHgJFmthhYD5yf\n0dFFRCRnGizozrl+dfzSeTnOIiIiWdBK0QYkEonQEbIS5fxRzg7KH1rU82fCMu3VpH0AM5fvY4iI\nxI2Z4fJ0UVRERIqcCrqISEyooIuIxIQKuohITKigi4jEhAq6iEhMqKCLiMSECrqISEyooIuIxIQK\nuohITKigi4jEhAq6iEhMqKCLiMSECrqISEyooIuIxESmm0SLSAn79FNYuhSWLYPly+Hzz2H1avj6\na9hqK2jeHFq0gI4doVMn6NwZ9toLmmgImVfa4EJEGvTNNzBpEkyZAq+8Ah9/DPvu6wt1p06w/fbQ\nsiVss41/7dq18OWX8P77vuAvXgxr1kCPHtCrF5xxBrRtG/p3Vdwy2eBCBV1EauUcJJPw2GMwbhx0\n7Qp9+sDRR/vHW2zRuPf7+9/9+734Ijz/PBx5JAwYAKedBk3VK/geFXQRydr69fDEE3DnnbBxIwwc\nCOecA+3b5+4YX3/t/5F48EFYuRKuuw4uuAC23jp3x4g6FXQRyVhVFTz0EPzXf8Hee8M110Dv3mCN\nKimNN2MG3HYbzJsHt9wC552nXjuooItIBpzzo+Ubb4Sdd4YhQ+CQQwqfY84cuPJKqK6GoUPhsMMK\nn6GY5GWTaDMbYWYVZraoll+7xsyqzax1Yw4qIsXhb3+Dn/0M/vM/4Z57/AXPEMUc4NBDYeZMGDTI\n99Wvv95fYJX0pfPBZiRwXM0nzawDcCzwYa5DiUh+bdwId9wBBx3kL04uWADHHZf/9kpDmjTxLZdF\ni+C99+Dgg2HhwrCZoqTBa8vOuelmtmstv3QXcB0wIeepRCRvVqyA88/3rZY5c2CPPUIn+r62beHp\np+Hxx30f//bb/UVTqV9Glx7M7GTgI+fc4hznEZE8GjPGj8r79IFXXy3OYr6JmR+tJ5P+YumgQVBZ\nGTpVcWv07E8z2xoYjG+3fPt0fT9TXl7+7eNEIkEikWjsYUUkC5WVcO21fv73iy9Ct26hE6Vv773h\njTd8cT/2WH8Bt2XL0KlyL5lMkkwms3qPtGa5pFouE51z+5vZvsAU4F/4Qt4BWAkc4pz7Ry0/q1ku\nIgGtWgVnnQXbbutbGK1ahU6UmepquOoqmDoVXnrJz8iJs7zMctn03qkvnHNLnHPtnHMdnXO7A38H\nDqitmItIWG+95WetHH00TJwY3WIO/oLp3XfDz38ORxwB77wTOlHxabDlYmajgQTQxsxWAGXOuZGb\nvcTRQMtFRArvhRf8hcT77oO+fUOnyQ0zP19+xx2hZ08/zbJz59CpiocWFonE0IMPQnk5PPMMdO8e\nOk1+PPIIDB7sL+7uuWfoNLmXSctFt8QRiRHn/IyQkSNh+vTinsWSrQED/Hz6Xr18X/3HPw6dKDwV\ndJGYcM7f5GrSJJg2DXbaKXSi/Bs40N+D5thj/SrTUvg910cFXSQGqqvhl7/09x1/7TVoXUI34/jF\nL/yGGyec4H/vLVqEThSOeugiEVdd7Yvau+/6C6HbbBM6UeE5B5ddBh98AM89B1tuGTpR9nS3RZES\nU10Nl17qp/CVajHfpKoKTj8d2rTxtwEOfV+abOVzHrqIFJlNo1IVc69pU78xx4IF/s6RpUg9dJEI\n2nQB9K23/D6fpV7MN2ne3N8a4LDD/J6nvXqFTlRYGqGLRNAf/gCTJ2tkXpvddoMnn4T+/f0m1aVE\nBV0kYu691y+qmTy5tGazNEYiAb/9re+pr1sXOk3h6KKoSISMHQu/+Y1fNLTbbqHTFDfn/ObWrVrB\nAw+ETtN4uigqEmOvvw6XX+6n5amYN8wM/vd//TWG//u/0GkKQyN0kQhYutTfMXH0aDjmmNBpomX+\nfL+93syZ8JOfhE6TPo3QRWJo1So48UT44x9VzDNx4IFQVgb9+sV/xyMVdJEitm4dnHoqXHghnHtu\n6DTR9atfwfbbw3//d+gk+aWWi0iRqq72F/WaNfM7DUV95WNoH38MBxwAEybAoYeGTtMwtVxEYqSs\nDFauhBEjVMxzYeed/ZTP886DtWtDp8kPFXSRIjRmDDz2GDz7LGy1Veg08dG3r9+S75ZbQifJD7Vc\nRIrMggXQuze8/DJ07Ro6TfysXu3/u912YXM0RDsWiUTcP/4Bp50G99+vYp4vxV7Is6ERukiRqKz0\nN5M66qj4z8aQhul+6CIRdvXV8PbbfiVoE13dKnlquYhE1BNPwPjxMHeuirlkTiN0kcCWLIGePf09\nR7p0CZ1GikVe5qGb2QgzqzCzRZs9d7uZLTOzhWb2tJmV8LasIplbswbOOAPuvFPFXLKXzoe7kcBx\nNZ6bDOzjnOsKLAduzHUwkbhzDi6+2I/OzzsvdBqJgwYLunNuOvBFjeemOOeqU9/OBjrkIZtIrA0b\n5nfUufvu0EkkLnJxUXQg8GQO3kekZMyZ47eRmzVLK0Eld7Iq6GZ2E1DpnBtd3+vKy8u/fZxIJEgk\nEtkcViTS/vlPOPtsGD4cOnYMnUaKRTKZJJlMZvUeac1yMbNdgYnOuf03e+4C4BLgaOfc+np+VrNc\nRFKc8ytBd98d7rordBopZvmch26pr00HOh64DvhpfcVcRL5r6FB/B8UxY0InkThqcIRuZqOBBNAG\nqADKgMHAlsDnqZfNds5dXsfPa4QuArz5JpxwAsyerVaLNExL/0WK1Jo1fiu0W2/1t3AVaYgKukgR\ncs7PM//hD/2FUJF06F4uIkXo0Uf9zvNvvhk6icSdRugiefTuu3DEEfDqq7DffqHTSJRoT1GRIrJh\nA/TrB+XlKuZSGBqhi+TJDTfAsmX+trja5FkaSz10kSLxyivw+OOwcKGKuRSOWi4iOfb55zBgADz8\nMLRtGzqNlBK1XERyyDl/f/Pdd4c//jF0GokytVxEAhsxAj74wG8pJ1JoGqGL5Mjy5XD44ZBMwj77\nhE4jUadpiyKBVFZC//5w880q5hKOCrpIDvz+99CmDVxxRegkUsrUQxfJ0qxZ/h4tCxZoiqKEpRG6\nSBa++grOPRf+9CfYaafQaaTU6aKoSBYuvhiqq+Ghh0InkbjRtEWRAho/HqZO9atBRYqBRugiGVi1\nCrp2haef9ndTFMk1bXAhUgDOwUknQZcu8Ic/hE4jcaWWi0gBDB8On3wCzzwTOonId2mELtIIy5dD\n9+4wbRp07hw6jcSZVoqK5FFVld8b9OabVcylOKmgi6Tp1luhRQutBpXipR66SBrmzoV77/WbPTfR\nMEiKVIN/Nc1shJlVmNmizZ5rZWaTzewdM5tkZtvlN6ZIOP/6l2+1DBsG7duHTiNSt3TGGiOB42o8\n9x/AFOfcnsCrwI25DiZSLK6/Hrp1g7POCp1EpH4Ntlycc9PNbNcaT58C9Eg9fgRI4ou8SKy89BJM\nnAhvvRU6iUjDMu2h7+CcqwBwzq0ysx1ymEmkKHz2GVx0kd/suWXL0GlEGpari6L1TjQvLy//9nEi\nkSCRSOTosCL54Rxceimccw707Bk6jZSCZDJJMpnM6j3SWliUarlMdM7tn/p+GZBwzlWYWTtgqnOu\n1pm5WlgkUfTII3DHHX52y1ZbhU4jpSifC4ss9bXJBOCC1OMBwPjGHFSkmP31r3DttTBqlIq5REuD\nI3QzGw0kgDZABVAGjAOeAnYBPgTOcs59WcfPa4QukbFxI/ToAaedBtdcEzqNlDLdbVEkS7fcAlOm\n+C8tIJKQVNBFsjBvHvTp4/+7yy6h00ip0825RDK0di306+dXg6qYS1RphC6Cn6K4bh08+mjoJCKe\nNrgQycD48fDyy9obVKJPI3QpaZ98AgccoL1Bpfiohy7SCNXVcMEFvt2iYi5xoIIuJevuu+Grr+C3\nvw2dRCQ31HKRkrRgARx3HMyZA7vvHjqNyPep5SKShrVr4ec/h3vuUTGXeNEIXUrORRdBZaWmKEpx\n07RFkQY8+SRMm+ZXg4rEjUboUjI++AAOPRQmTYIDDwydRqR+6qGL1KGy0vfNb7pJxVziSyN0KQnX\nXgtvv+33B7VGjXlEwlAPXaQWzz0HTz0F8+ermEu8qaBLrH30kZ/V8uyz0KZN6DQi+aUeusRWZaXf\n5Pnqq+Hww0OnEck/9dAltq67DpYu9S0X7T4kUaMeukjKs8/6vvm8eSrmUjo0QpfYef996N7dj8wP\nOSR0GpHMaB66lLx16+DMM6GsTMVcSo9G6BIbzvn7m1dWwqhRmqIo0aYeupS0++/3t8WdNUvFXEpT\nViN0M7sROBfYCCwGLnTObajxGo3QJe+mT4czzoCZM2GPPUKnEcleQXvoZrYrcAlwgHNuf/xo/5xM\n308kUytXwtlnw8MPq5hLacum5bIG2AA0N7Nq4IfAxzlJJZKmdevgtNPg8suhT5/QaUTCyrblcglw\nJ/AvYLJz7rxaXqOWi+SFc3D++f4i6BNPqG8u8VLQi6Jm1hG4CtgVWA2MNbN+zrnRNV9bXl7+7eNE\nIkEikcj0sCLfuuMOvxJ0+nQVc4m+ZDJJMpnM6j0yHqGb2VnAsc65S1Lfnwcc6py7osbrNEKXnJs4\nES69FGbPhh/9KHQakdwr9MKid4DDzGwrMzPgGGBZFu8nkpaFC2HgQL+8X8Vc5N8yLujOubeAR4F5\nwFuAAcNzlEukVh9/DCefDPfd57eTE5F/00pRiYyvv4ZEAk4/HQYPDp1GJL8yabmooEskVFbCKafA\nTjvBn/+si6ASf1r6X4/qar97zTff+Mdm0L49bLtt6GTSEOf8BVAzeOABFXORusS2oFdVwSuvwNix\nfi/Jt9+GVq2geXN/f+yNG/0Kw+bN/erCww+Hnj3hqKNgu+1Cp5fNlZXBkiUwdSo0axY6jUjxil3L\n5dNPYcgQf7e9XXf1W5B17w577w0tWnz3tc5BRQW88w5Mm+YLxhtv+KLev7//iL/NNgWLLrW45x64\n916YMQN22CF0GpHCKeke+jffwNChcPvt0K8f/PrX8OMfN/59vv4axo/3Kw9nzPj3e+25Z+4zS/0e\nfhhuvtn/Y7vrrqHTiBRWyW5wsWQJdOni77Q3c6Yv7JkUc/Aj8v79/W43S5dC69bw05/CCSf495bC\neOYZuPFGmDxZxVwkXZEfoT/1lL8x0113wbnn5ucY33wDjzwCt9wCe+0F5eW+jSP58dxzcNFF8NJL\ncMABodOIhFFSLRfn/Mfxxx+Hp5+GAw/M+SG+Z8MGX9h//3vo1s336tWKya3nn4cLL9R+oCIl1XIp\nK/Mn/dy5hSnmAFtuCZdcAu++60foRx4JV1wBn39emOPH3Qsv+GI+caKKuUgmIlnQ/+d/fKtl8mTY\nfvvCH3+rreD662HZMj8nunNnv/1ZVVXhs8TFs8/6/UAnTNCSfpFMRa7lMny4b3W8/jp06JCzt83K\n4sV+Jsw//+kL+xFHhE4ULY895v+BfP75wn3aEil2se+hz57tb8w0c2bms1jyxTkYMwauuQaOPRZu\nu03zptNx//1w663+01bnzqHTiBSPWPfQv/jCLxIaPrz4ijn41svZZ/s2TOvWsO++Pmt1dehkxam6\n2t9g6667/KctFXOR7EVihO6c39F9l138ysEoWLQIfvlLn/1Pf4KuXUMnKh7r1/uLn3/7m++Zh7gO\nIlLsYjtCv+8+WLHCrwKNiv3391ujXXQRHHccXHUVfPVV6FThffop9O7ti/orr6iYi+RS0Rf0jz7y\nC3meeAJ+8IPQaRqnSRO4+GK/knX1an8/maee8qP2UrRgARx8sL8R2pgxsPXWoROJxEvRt1z69oV9\n9vFFPeqmTfOrWtu1g2HD/KrTUjF6NFx5pb8I2rdv6DQixS92LZeXX4Z58+CGG0InyY2jjvK38j3x\nRL8o6YYb4t+GWbvWf0opK/MtFhVzkfwp2oK+fr1fhTl0aLw+mjdrBr/5jZ+7XlHhbx0wcmQ8Z8Ms\nWgQHHeRvmTB/vr+uICL5U7Qtl9tu87evnTAhD6GKyBtv+AK/fr2/6HvMMaETZa+y0v//u/tuuPNO\nOP/80IlEoic2C4u+/NLPNZ81C37ykzwFKyLO+Yulgwf73/eQIdGd5rhgAQwc6K8TPPgg/OhHoROJ\nRFNseuhDh8LPflYaxRz8oqSzzoK//MX31/v0gTPP9G2ZqPjsM7jsMjj+eP+J44UXVMxFCq3oCvrq\n1X4GyE03hU5SeFtuCYMGwXvvwWGH+VsInHkmzJkTOlnd1q3zbZXOnX3+ZctgwABt5CwSQlYF3cy2\nM7OnzGyZmS01s6zvkzdsmB+hlsrovDbNm8O118L77/uZMWef7XdNGjeueO7ouG6dX7W7xx5+6f5r\nr/nvW7cOnUykdGXVQzezh4HXnHMjzawp8EPn3Joar0m7h75mje8hT58OnTplHCt2qqpg7Fj/j92H\nH/rVpwMHhtmabcUKeOABGDHC3xO+rEy7ConkQ0F76GbWAjjKOTcSwDlXVbOYN9Z99/ll4Srm39W0\nqb8x2YwZ8OKLfkONbt38bXqHDYOVK/N7/DVrYNQoOOkkX7zXrvWj8nHjVMxFiknGI3Qz6wIMB/4C\ndAHeBK50zq2r8bq0RuiVlX7EOXmyv1Oh1G/DBpgyBZ580t9HfOed/QXJnj398vq2bTN/7+pqf0H2\n1Vf9YqBp03zLp29fOP10v5G2iORXQactmlk3YDbQ3Tn3ppndDax2zpXVeJ0rK/v3U4lEgkQi8b33\nGzPGLwtPJjOKU9I2boQ334RJk3zxnTvX97L339/3uPfYA9q3h5Yt/VezZr6NU1Xlp4iuWgWffOK3\n1lu82N97pl07OPpo/w9E797QqlXo36VIvCWTSZKbFcDf/e53BS3oOwKznHMdU98fCdzgnDupxuvS\nGqH36OFXhmppePaqq31xXrrUX1h97z2/KvXLL/195auqfBtniy18gW/XDnbc0V+/2G8//wmpTZvQ\nvwuR0pbJCL1ppgdzzlWY2Udm1sk59y5wDL790mhLlsDy5XDqqZmmkc01aeJv/FVKN/8SkSwKesqv\ngVFm1gz4ALgwkze5/374xS98K0BERDITfOn/mjWw225+lL7zznmNIiISGZFc+v/oo9Crl4q5iEi2\nghf0hx7y7RYREclO0IK+bJmffdGzZ8gUIiLxELSgjxrlV0BusUXIFCIi8ZDtLJeMOef3mRw7NlQC\nEZF4CTZCnz0bfvAD3QtERCRXghX0UaOgf3/dN1tEJFeCzEOvrPT3Fpk9Gzp2zOvhRUQiKTLz0F9+\n2d8wSsVcRCR3ghT0J5+Efv1CHFlEJL4K3nKpqvJ391u4EDp0yOuhRUQiKxItlxkz/EYWKuYiIrlV\n8II+YQKcfHKhjyoiEn8FLejOwfjxcMophTyqiEhpKGhBX7bM74XZpUshjyoiUhoKWtA3tVu0mEhE\nJPeCFHQREcm9gk1brKjwe1xWVMCWW+b1kCIikVfU0xafew5691YxFxHJl4IV9BdegBNPLNTRRERK\nT0FaLlVVjrZtYelS2GmnvB5ORCQWirblMn++3wRaxVxEJH8KUtCnTIFevQpxJBGR0pV1QTezJmY2\n38wm1PUaFXQRkfzLuoduZlcB3YAWzrnvzTI3M9e8ueOTT2DbbbM6lIhIySh4D93MOgAnAH+u73Vd\nu6qYi4jkW7Ytl7uA64B6h/lqt4iI5F/TTH/QzE4EKpxzC80sAdT50WDFinLKy/3jRCJBIpHI9LAi\nIrGUTCZJJpNZvUfGPXQzuwU4F6gCtga2BZ5xzp1f43VuwwZHs2ZZ5RQRKSmZ9NBzsrDIzHoA19R1\nUTTfi5dEROKmaBcWiYhI/hV8k2gREWmYRugiIiVMBV1EJCZU0EVEYkIFXUQkJlTQRURiQgVdRCQm\nVNBFRGJCBV1EJCZU0EVEYkIFXUQkJlTQRURiQgVdRCQmVNBFRGJCBV1EJCZU0EVEYkIFXUQkJlTQ\nRURiQgVdRCQmVNBFRGJCBV1EJCZU0EVEYkIFXUQkJjIu6GbWwcxeNbOlZrbYzH6dy2AiItI42YzQ\nq4CrnXP7AN2BX5nZXrmJVTySyWToCFmJcv4oZwflDy3q+TORcUF3zq1yzi1MPf4aWAa0z1WwYhH1\nvxRRzh/l7KD8oUU9fyZy0kM3s92ArsCcXLyfiIg0XtYF3cy2AcYCV6ZG6iIiEoA55zL/YbOmwHPA\ni865e+p4TeYHEBEpYc45a8zrsy3ojwKfOeeuzvhNREQkJzIu6GZ2BPA6sBhwqa/BzrmXchdPRETS\nldUIXUREikdOV4qa2QgzqzCzRZs918rMJpvZO2Y2ycy2y+Uxc6mO/Leb2TIzW2hmT5tZi5AZ61Jb\n9s1+7Rozqzaz1iGypaOu/GY2KPXnv9jMhoTK15A6/u4cbGZvmNmC1H8PCpmxLnUtEozKuVtL/kGp\n56Ny7ta7SLNR569zLmdfwJH46YuLNnvuNuD61OMbgCG5PGYB8vcCmqQeDwFuDZ0z3eyp5zsALwF/\nBVqHztnIP/sEMBlomvp++9A5G5l/KtA79bgPMDV0zjqytwO6ph5vA7wD7BWVc7ee/FE5d2vNn/q+\nUedvTkfozrnpwBc1nj4FeCT1+BHg1FweM5dqy++cm+Kcq059Oxv/B1x06vizB7gLuK7AcRqtjvyX\n4YtIVeo1nxU8WJrqyP8JsGlU2xJYWdBQaXK1LxLsQETO3Tryt4/QuVvfIs1Gnb+FuDnXDs65CvDB\ngR0KcMx8GQi8GDpEuszsZOAj59zi0Fky1An4qZnNNrOpxdqyqMd/AHea2QrgduDGwHkatNkiwdnA\njlE7d+tZ5BiJc3fz/Jmcv03zlKs+kbwKa2Y3AZXOudGhs6TDzLYGBgPHbv50oDiZagq0cs4dZmYH\nA2OAjoEzNcYIYJBzbpyZnQk8xHf/fxSVmosEa1lDUtTnbl2LHKNy7m6eH9hIBudvIUboFWa2I4CZ\ntQP+UYBj5pSZXQCcAPQLHKUx9gB2A94ys7/iP27OM7OiH2Vt5iPgGQDn3Fyg2szahI3UKIc658YB\nOOfGAocEzlOn1CLBscBjzrnxqacjc+7WkT8y524t+TM6f/NR0I3v/ksyAbgg9XgAML7mDxSZ7+Q3\ns+PxPayTnXPrg6VKz7fZnXNLnHPtnHMdnXO7A38HDnDOFe1Jyff/7owDjgYws05AM+fc5yGCpalm\n/uVm1gPAzI4B3g2SKj0PAX9x313xHaVz93v5I3bufid/xudvjq/WjgY+BtYDK4ALgVbAFPyV28lA\ny9BXlRuZfznwITA/9XV/6JzpZq/x6x9Q3LNcavuzbwo8hl+89ibQI3TORubvhu/lLgBmpU7I4Flr\nyX4E/iP+wlTW+cDxQOsonLt15O8ToXO31j//Gq9J6/zVwiIRkZjQFnQiIjGhgi4iEhMq6CIiMaGC\nLiISEyroIiIxoYIuIhITKugiIjGhgi4iEhP/D6H32MBxWYP1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x45c6950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = -2.88192413434e-5\n",
    "b = 0.00387383125623\n",
    "c = -0.20721471544\n",
    "d = 5.64382310405\n",
    "e = -82.6061154979\n",
    "f = 617.464778362\n",
    "g = -1840.49386605\n",
    "\n",
    "y = sextic(t, a, b, c, d, e, f, g)\n",
    "plt.plot(t, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}