Star_Formation_History.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 530,
   "metadata": {},
   "outputs": [],
   "source": [
    "#this cell has what you will need to import\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import optimize as optimise\n",
    "from scipy import signal, ndimage\n",
    "import re\n",
    "import fnmatch\n",
    "\n",
    "#I use these to loop over the directories\n",
    "import os\n",
    "from pathlib import Path\n",
    "\n",
    "from astropy.io import fits #to read a fits file\n",
    "from astropy.convolution import convolve, Gaussian1DKernel # to smooth histogram \n",
    "\n",
    "\n",
    "#I use this package to convert from redshift to year\n",
    "#Renyue implied doing this conversion is trivial\n",
    "#but it is not. I suggest reading up on this\n",
    "#package and then using it.\n",
    "from astropy.cosmology import FlatLambdaCDM\n",
    "cosmo = FlatLambdaCDM(H0=70, Om0=0.272)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_sun(director):\n",
    "    \n",
    "    \n",
    "    masses=[] \n",
    "    sfrs=[]\n",
    "    metals=[]\n",
    "    ages=[]\n",
    "    \n",
    "    for filename in os.listdir(director):\n",
    "        \n",
    "        if filename.startswith('sun'):\n",
    "            \n",
    "            sfile=director+'/'+filename\n",
    "                    \n",
    "            hdul = fits.open(sfile)\n",
    "            data=hdul[4].data                \n",
    "            hdul.close()\n",
    "                    \n",
    "            #here I calculate the star formation rate as the total\n",
    "            #mass of star particles less than a Myr old\n",
    "            good_age=data['age']<1e6\n",
    "            sfr = sum(data['mass'][good_age])/1e6\n",
    "            \n",
    "            masses.append(np.array(data['mass']))\n",
    "            sfrs.append(sfr)\n",
    "            metals.append(np.array(data['metallicity']))\n",
    "            ages.append(np.array(data['age']))\n",
    "                    \n",
    "\n",
    "    masses=np.array(masses)\n",
    "    sfrs=np.array(sfrs)\n",
    "    metals=np.array(metals)\n",
    "    ages=np.array(ages)\n",
    "    \n",
    "\n",
    "    return masses,sfrs,metals,ages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 534,
   "metadata": {},
   "outputs": [],
   "source": [
    "mass,sfrs,metals,ages=read_sun('out01461')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 535,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_mass_hist(galaxy):\n",
    "    \n",
    "    \"\"\"\n",
    "    Plots a histogram of star particle age, weighted by mass\n",
    "    Parameters:\n",
    "        galaxy: galaxy number (0-585)\n",
    "    Returns figure object\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    max_age = int(max(ages[galaxy]))\n",
    "    digits = len(str(max_age))\n",
    "    max_bin = (int(np.ceil(max_age / 10**(digits-1))) + 1) * 100\n",
    "    fig = plt.figure()\n",
    "    plt.hist(ages[galaxy], bins = np.arange(max_bin)*1e6, weights = mass[galaxy])\n",
    "    plt.title(\"Histogram of star particles ages in galaxy %i, weighted by mass\" % galaxy)\n",
    "    plt.ylabel(\"Total star particle mass ($M_\\odot$)\")\n",
    "    plt.xlabel(\"Star particle age (years)\")\n",
    "    return fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 536,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_mass_hist(galaxy):\n",
    "    \n",
    "    \"\"\"\n",
    "    Calculates a histogram of star particle age, weighted by mass\n",
    "    Parameters:\n",
    "        galaxy: galaxy number (0-585)\n",
    "    Returns mass value of each bin (n) and the left edges of the bins + the right edge of the last bin (bins)\n",
    "    \n",
    "    \"\"\"\n",
    "    max_age = int(max(ages[galaxy]))\n",
    "    digits = len(str(max_age))\n",
    "    max_bin = int(np.ceil(max_age / 10**(digits-1))) * 100\n",
    "    n, bins = np.histogram(ages[galaxy], bins = np.arange(max_bin)*1e6, weights = mass[galaxy])\n",
    "    return n,bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 537,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sfr(galaxy):\n",
    "    \n",
    "    \"\"\" \n",
    "    Calculates star formation rate by dividing the mass value of each histogram bin by bin width\n",
    "    Parameters:\n",
    "        galaxy: galaxy number (0-585)\n",
    "    Returns array of star formation values, and the (n,bins) output of the histogram calculation\n",
    "      \n",
    "    \"\"\"\n",
    "    \n",
    "    max_age = int(max(ages[galaxy]))\n",
    "    digits = len(str(max_age))\n",
    "    max_bin = (int(np.ceil(max_age / 10**(digits-1))) +1) * 100\n",
    "    n,bins = np.histogram(ages[galaxy], bins = np.arange(max_bin)*1e6, weights = mass[galaxy])\n",
    "    bin_width = bins[1] - bins[0]\n",
    "    sfr_array = n/bin_width\n",
    "    return np.array(sfr_array), n, bins\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n",
      "29\n",
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n",
      "35\n",
      "36\n",
      "37\n",
      "38\n",
      "39\n",
      "40\n",
      "41\n",
      "42\n",
      "43\n",
      "44\n",
      "45\n",
      "46\n",
      "47\n",
      "48\n",
      "49\n",
      "50\n",
      "51\n",
      "52\n",
      "53\n",
      "54\n",
      "55\n",
      "56\n",
      "57\n",
      "58\n",
      "59\n",
      "60\n",
      "61\n",
      "62\n",
      "63\n",
      "64\n",
      "65\n",
      "66\n",
      "67\n",
      "68\n",
      "69\n",
      "70\n",
      "71\n",
      "72\n",
      "73\n",
      "74\n",
      "75\n",
      "76\n",
      "77\n",
      "78\n",
      "79\n",
      "80\n",
      "81\n",
      "82\n",
      "83\n",
      "84\n",
      "85\n",
      "86\n",
      "87\n",
      "88\n",
      "89\n",
      "90\n",
      "91\n",
      "92\n",
      "93\n",
      "94\n",
      "95\n",
      "96\n",
      "97\n",
      "98\n",
      "99\n",
      "100\n",
      "101\n",
      "102\n",
      "103\n",
      "104\n",
      "105\n",
      "106\n",
      "107\n",
      "108\n",
      "109\n",
      "110\n",
      "111\n",
      "112\n",
      "113\n",
      "114\n",
      "115\n",
      "116\n",
      "117\n",
      "118\n",
      "119\n",
      "120\n",
      "121\n",
      "122\n",
      "123\n",
      "124\n",
      "125\n",
      "126\n",
      "127\n",
      "128\n",
      "129\n",
      "130\n",
      "131\n",
      "132\n",
      "133\n",
      "134\n",
      "135\n",
      "136\n",
      "137\n",
      "138\n",
      "139\n",
      "140\n",
      "141\n",
      "142\n",
      "143\n",
      "144\n",
      "145\n",
      "146\n",
      "147\n",
      "148\n",
      "149\n",
      "150\n",
      "151\n",
      "152\n",
      "153\n",
      "154\n",
      "155\n",
      "156\n",
      "157\n",
      "158\n",
      "159\n",
      "160\n",
      "161\n",
      "162\n",
      "163\n",
      "164\n",
      "165\n",
      "166\n",
      "167\n",
      "168\n",
      "169\n",
      "170\n",
      "171\n",
      "172\n",
      "173\n",
      "174\n",
      "175\n",
      "176\n",
      "177\n",
      "178\n",
      "179\n",
      "180\n",
      "181\n",
      "182\n",
      "183\n",
      "184\n",
      "185\n",
      "186\n",
      "187\n",
      "188\n",
      "189\n",
      "190\n",
      "191\n",
      "192\n",
      "193\n",
      "194\n",
      "195\n",
      "196\n",
      "197\n",
      "198\n",
      "199\n",
      "200\n",
      "201\n",
      "202\n",
      "203\n",
      "204\n",
      "205\n",
      "206\n",
      "207\n",
      "208\n",
      "209\n",
      "210\n",
      "211\n",
      "212\n",
      "213\n",
      "214\n",
      "215\n",
      "216\n",
      "217\n",
      "218\n",
      "219\n",
      "220\n",
      "221\n",
      "222\n",
      "223\n",
      "224\n",
      "225\n",
      "226\n",
      "227\n",
      "228\n",
      "229\n",
      "230\n",
      "231\n",
      "232\n",
      "233\n",
      "234\n",
      "235\n",
      "236\n",
      "237\n",
      "238\n",
      "239\n",
      "240\n",
      "241\n",
      "242\n",
      "243\n",
      "244\n",
      "245\n",
      "246\n",
      "247\n",
      "248\n",
      "249\n",
      "250\n",
      "251\n",
      "252\n",
      "253\n",
      "254\n",
      "255\n",
      "256\n",
      "257\n",
      "258\n",
      "259\n",
      "260\n",
      "261\n",
      "262\n",
      "263\n",
      "264\n",
      "265\n",
      "266\n",
      "267\n",
      "268\n",
      "269\n",
      "270\n",
      "271\n",
      "272\n",
      "273\n",
      "274\n",
      "275\n",
      "276\n",
      "277\n",
      "278\n",
      "279\n",
      "280\n",
      "281\n",
      "282\n",
      "283\n",
      "284\n",
      "285\n",
      "286\n",
      "287\n",
      "288\n",
      "289\n",
      "290\n",
      "291\n",
      "292\n",
      "293\n",
      "294\n",
      "295\n",
      "296\n",
      "297\n",
      "298\n",
      "299\n",
      "300\n",
      "301\n",
      "302\n",
      "303\n",
      "304\n",
      "305\n",
      "306\n",
      "307\n",
      "308\n",
      "309\n",
      "310\n",
      "311\n",
      "312\n",
      "313\n",
      "314\n",
      "315\n",
      "316\n",
      "317\n",
      "318\n",
      "319\n",
      "320\n",
      "321\n",
      "322\n",
      "323\n",
      "324\n",
      "325\n",
      "326\n",
      "327\n",
      "328\n",
      "329\n",
      "330\n",
      "331\n",
      "332\n",
      "333\n",
      "334\n",
      "335\n",
      "336\n",
      "337\n",
      "338\n",
      "339\n",
      "340\n",
      "341\n",
      "342\n",
      "343\n",
      "344\n",
      "345\n",
      "346\n",
      "347\n",
      "348\n",
      "349\n",
      "350\n",
      "351\n",
      "352\n",
      "353\n",
      "354\n",
      "355\n",
      "356\n",
      "357\n",
      "358\n",
      "359\n",
      "360\n",
      "361\n",
      "362\n",
      "363\n",
      "364\n",
      "365\n",
      "366\n",
      "367\n",
      "368\n",
      "369\n",
      "370\n",
      "371\n",
      "372\n",
      "373\n",
      "374\n",
      "375\n",
      "376\n",
      "377\n",
      "378\n",
      "379\n",
      "380\n",
      "381\n",
      "382\n",
      "383\n",
      "384\n",
      "385\n",
      "386\n",
      "387\n",
      "388\n",
      "389\n",
      "390\n",
      "391\n",
      "392\n",
      "393\n",
      "394\n",
      "395\n",
      "396\n",
      "397\n",
      "398\n",
      "399\n",
      "400\n",
      "401\n",
      "402\n",
      "403\n",
      "404\n",
      "405\n",
      "406\n",
      "407\n",
      "408\n",
      "409\n",
      "410\n",
      "411\n",
      "412\n",
      "413\n",
      "414\n",
      "415\n",
      "416\n",
      "417\n",
      "418\n",
      "419\n",
      "420\n",
      "421\n",
      "422\n",
      "423\n",
      "424\n",
      "425\n",
      "426\n",
      "427\n",
      "428\n",
      "429\n",
      "430\n",
      "431\n",
      "432\n",
      "433\n",
      "434\n",
      "435\n",
      "436\n",
      "437\n",
      "438\n",
      "439\n",
      "440\n",
      "441\n",
      "442\n",
      "443\n",
      "444\n",
      "445\n",
      "446\n",
      "447\n",
      "448\n",
      "449\n",
      "450\n",
      "451\n",
      "452\n",
      "453\n",
      "454\n",
      "455\n",
      "456\n",
      "457\n",
      "458\n",
      "459\n",
      "460\n",
      "461\n",
      "462\n",
      "463\n",
      "464\n",
      "465\n",
      "466\n",
      "467\n",
      "468\n",
      "469\n",
      "470\n",
      "471\n",
      "472\n",
      "473\n",
      "474\n",
      "475\n",
      "476\n",
      "477\n",
      "478\n",
      "479\n",
      "480\n",
      "481\n",
      "482\n",
      "483\n",
      "484\n",
      "485\n",
      "486\n",
      "487\n",
      "488\n",
      "489\n",
      "490\n",
      "491\n",
      "492\n",
      "493\n",
      "494\n",
      "495\n",
      "496\n",
      "497\n",
      "498\n",
      "499\n",
      "500\n",
      "501\n",
      "502\n",
      "503\n",
      "504\n",
      "505\n",
      "506\n",
      "507\n",
      "508\n",
      "509\n",
      "510\n",
      "511\n",
      "512\n",
      "513\n",
      "514\n",
      "515\n",
      "516\n",
      "517\n",
      "518\n",
      "519\n",
      "520\n",
      "521\n",
      "522\n",
      "523\n",
      "524\n",
      "525\n",
      "526\n",
      "527\n",
      "528\n",
      "529\n",
      "530\n",
      "531\n",
      "532\n",
      "533\n",
      "534\n",
      "535\n",
      "536\n",
      "537\n",
      "538\n",
      "539\n",
      "540\n",
      "541\n",
      "542\n",
      "543\n",
      "544\n",
      "545\n",
      "546\n",
      "547\n",
      "548\n",
      "549\n",
      "550\n",
      "551\n",
      "552\n",
      "553\n",
      "554\n",
      "555\n",
      "556\n",
      "557\n",
      "558\n",
      "559\n",
      "560\n",
      "561\n",
      "562\n",
      "563\n",
      "564\n",
      "565\n",
      "566\n",
      "567\n",
      "568\n",
      "569\n",
      "570\n",
      "571\n",
      "572\n",
      "573\n",
      "574\n",
      "575\n",
      "576\n",
      "577\n",
      "578\n",
      "579\n",
      "580\n",
      "581\n",
      "582\n",
      "583\n",
      "584\n",
      "585\n"
     ]
    }
   ],
   "source": [
    "# Saves Histograms in folder\n",
    "for i in range(len(mass)):\n",
    "    print(i)\n",
    "    fig = plot_mass_hist(i)\n",
    "    fig.savefig('Star Formation Histograms/galaxy_%i.png' %i)\n",
    "    plt.tight_layout()\n",
    "    fig.clf()\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 566,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n",
      "29\n",
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n",
      "35\n",
      "36\n",
      "37\n",
      "38\n",
      "39\n",
      "40\n",
      "41\n",
      "42\n",
      "43\n",
      "44\n",
      "45\n",
      "46\n",
      "47\n",
      "48\n",
      "49\n",
      "50\n",
      "51\n",
      "52\n",
      "53\n",
      "54\n",
      "55\n",
      "56\n",
      "57\n",
      "58\n",
      "59\n",
      "60\n",
      "61\n",
      "62\n",
      "63\n",
      "64\n",
      "65\n",
      "66\n",
      "67\n",
      "68\n",
      "69\n",
      "70\n",
      "71\n",
      "72\n",
      "73\n",
      "74\n",
      "75\n",
      "76\n",
      "77\n",
      "78\n",
      "79\n",
      "80\n",
      "81\n",
      "82\n",
      "83\n",
      "84\n",
      "85\n",
      "86\n",
      "87\n",
      "88\n",
      "89\n",
      "90\n",
      "91\n",
      "92\n",
      "93\n",
      "94\n",
      "95\n",
      "96\n",
      "97\n",
      "98\n",
      "99\n",
      "100\n",
      "101\n",
      "102\n",
      "103\n",
      "104\n",
      "105\n",
      "106\n",
      "107\n",
      "108\n",
      "109\n",
      "110\n",
      "111\n",
      "112\n",
      "113\n",
      "114\n",
      "115\n",
      "116\n",
      "117\n",
      "118\n",
      "119\n",
      "120\n",
      "121\n",
      "122\n",
      "123\n",
      "124\n",
      "125\n",
      "126\n",
      "127\n",
      "128\n",
      "129\n",
      "130\n",
      "131\n",
      "132\n",
      "133\n",
      "134\n",
      "135\n",
      "136\n",
      "137\n",
      "138\n",
      "139\n",
      "140\n",
      "141\n",
      "142\n",
      "143\n",
      "144\n",
      "145\n",
      "146\n",
      "147\n",
      "148\n",
      "149\n",
      "150\n",
      "151\n",
      "152\n",
      "153\n",
      "154\n",
      "155\n",
      "156\n",
      "157\n",
      "158\n",
      "159\n",
      "160\n",
      "161\n",
      "162\n",
      "163\n",
      "164\n",
      "165\n",
      "166\n",
      "167\n",
      "168\n",
      "169\n",
      "170\n",
      "171\n",
      "172\n",
      "173\n",
      "174\n",
      "175\n",
      "176\n",
      "177\n",
      "178\n",
      "179\n",
      "180\n",
      "181\n",
      "182\n",
      "183\n",
      "184\n",
      "185\n",
      "186\n",
      "187\n",
      "188\n",
      "189\n",
      "190\n",
      "191\n",
      "192\n",
      "193\n",
      "194\n",
      "195\n",
      "196\n",
      "197\n",
      "198\n",
      "199\n",
      "200\n",
      "201\n",
      "202\n",
      "203\n",
      "204\n",
      "205\n",
      "206\n",
      "207\n",
      "208\n",
      "209\n",
      "210\n",
      "211\n",
      "212\n",
      "213\n",
      "214\n",
      "215\n",
      "216\n",
      "217\n",
      "218\n",
      "219\n",
      "220\n",
      "221\n",
      "222\n",
      "223\n",
      "224\n",
      "225\n",
      "226\n",
      "227\n",
      "228\n",
      "229\n",
      "230\n",
      "231\n",
      "232\n",
      "233\n",
      "234\n",
      "235\n",
      "236\n",
      "237\n",
      "238\n",
      "239\n",
      "240\n",
      "241\n",
      "242\n",
      "243\n",
      "244\n",
      "245\n",
      "246\n",
      "247\n",
      "248\n",
      "249\n",
      "250\n",
      "251\n",
      "252\n",
      "253\n",
      "254\n",
      "255\n",
      "256\n",
      "257\n",
      "258\n",
      "259\n",
      "260\n",
      "261\n",
      "262\n",
      "263\n",
      "264\n",
      "265\n",
      "266\n",
      "267\n",
      "268\n",
      "269\n",
      "270\n",
      "271\n",
      "272\n",
      "273\n",
      "274\n",
      "275\n",
      "276\n",
      "277\n",
      "278\n",
      "279\n",
      "280\n",
      "281\n",
      "282\n",
      "283\n",
      "284\n",
      "285\n",
      "286\n",
      "287\n",
      "288\n",
      "289\n",
      "290\n",
      "291\n",
      "292\n",
      "293\n",
      "294\n",
      "295\n",
      "296\n",
      "297\n",
      "298\n",
      "299\n",
      "300\n",
      "301\n",
      "302\n",
      "303\n",
      "304\n",
      "305\n",
      "306\n",
      "307\n",
      "308\n",
      "309\n",
      "310\n",
      "311\n",
      "312\n",
      "313\n",
      "314\n",
      "315\n",
      "316\n",
      "317\n",
      "318\n",
      "319\n",
      "320\n",
      "321\n",
      "322\n",
      "323\n",
      "324\n",
      "325\n",
      "326\n",
      "327\n",
      "328\n",
      "329\n",
      "330\n",
      "331\n",
      "332\n",
      "333\n",
      "334\n",
      "335\n",
      "336\n",
      "337\n",
      "338\n",
      "339\n",
      "340\n",
      "341\n",
      "342\n",
      "343\n",
      "344\n",
      "345\n",
      "346\n",
      "347\n",
      "348\n",
      "349\n",
      "350\n",
      "351\n",
      "352\n",
      "353\n",
      "354\n",
      "355\n",
      "356\n",
      "357\n",
      "358\n",
      "359\n",
      "360\n",
      "361\n",
      "362\n",
      "363\n",
      "364\n",
      "365\n",
      "366\n",
      "367\n",
      "368\n",
      "369\n",
      "370\n",
      "371\n",
      "372\n",
      "373\n",
      "374\n",
      "375\n",
      "376\n",
      "377\n",
      "378\n",
      "379\n",
      "380\n",
      "381\n",
      "382\n",
      "383\n",
      "384\n",
      "385\n",
      "386\n",
      "387\n",
      "388\n",
      "389\n",
      "390\n",
      "391\n",
      "392\n",
      "393\n",
      "394\n",
      "395\n",
      "396\n",
      "397\n",
      "398\n",
      "399\n",
      "400\n",
      "401\n",
      "402\n",
      "403\n",
      "404\n",
      "405\n",
      "406\n",
      "407\n",
      "408\n",
      "409\n",
      "410\n",
      "411\n",
      "412\n",
      "413\n",
      "414\n",
      "415\n",
      "416\n",
      "417\n",
      "418\n",
      "419\n",
      "420\n",
      "421\n",
      "422\n",
      "423\n",
      "424\n",
      "425\n",
      "426\n",
      "427\n",
      "428\n",
      "429\n",
      "430\n",
      "431\n",
      "432\n",
      "433\n",
      "434\n",
      "435\n",
      "436\n",
      "437\n",
      "438\n",
      "439\n",
      "440\n",
      "441\n",
      "442\n",
      "443\n",
      "444\n",
      "445\n",
      "446\n",
      "447\n",
      "448\n",
      "449\n",
      "450\n",
      "451\n",
      "452\n",
      "453\n",
      "454\n",
      "455\n",
      "456\n",
      "457\n",
      "458\n",
      "459\n",
      "460\n",
      "461\n",
      "462\n",
      "463\n",
      "464\n",
      "465\n",
      "466\n",
      "467\n",
      "468\n",
      "469\n",
      "470\n",
      "471\n",
      "472\n",
      "473\n",
      "474\n",
      "475\n",
      "476\n",
      "477\n",
      "478\n",
      "479\n",
      "480\n",
      "481\n",
      "482\n",
      "483\n",
      "484\n",
      "485\n",
      "486\n",
      "487\n",
      "488\n",
      "489\n",
      "490\n",
      "491\n",
      "492\n",
      "493\n",
      "494\n",
      "495\n",
      "496\n",
      "497\n",
      "498\n",
      "499\n",
      "500\n",
      "501\n",
      "502\n",
      "503\n",
      "504\n",
      "505\n",
      "506\n",
      "507\n",
      "508\n",
      "509\n",
      "510\n",
      "511\n",
      "512\n",
      "513\n",
      "514\n",
      "515\n",
      "516\n",
      "517\n",
      "518\n",
      "519\n",
      "520\n",
      "521\n",
      "522\n",
      "523\n",
      "524\n",
      "525\n",
      "526\n",
      "527\n",
      "528\n",
      "529\n",
      "530\n",
      "531\n",
      "532\n",
      "533\n",
      "534\n",
      "535\n",
      "536\n",
      "537\n",
      "538\n",
      "539\n",
      "540\n",
      "541\n",
      "542\n",
      "543\n",
      "544\n",
      "545\n",
      "546\n",
      "547\n",
      "548\n",
      "549\n",
      "550\n",
      "551\n",
      "552\n",
      "553\n",
      "554\n",
      "555\n",
      "556\n",
      "557\n",
      "558\n",
      "559\n",
      "560\n",
      "561\n",
      "562\n",
      "563\n",
      "564\n",
      "565\n",
      "566\n",
      "567\n",
      "568\n",
      "569\n",
      "570\n",
      "571\n",
      "572\n",
      "573\n",
      "574\n",
      "575\n",
      "576\n",
      "577\n",
      "578\n",
      "579\n",
      "580\n",
      "581\n",
      "582\n",
      "583\n",
      "584\n",
      "585\n"
     ]
    }
   ],
   "source": [
    "# Saves plots of star formation rate in folder, on a log scale\n",
    "\n",
    "for i in range(len(mass)):\n",
    "    print(i)\n",
    "    sfr_array, n, bins = sfr(i)\n",
    "    fig = plt.figure()\n",
    "    plt.plot(bins[:-1], sfr_array)\n",
    "    plt.title(\"Star formation rate of star particles in galaxy %i\" % i)\n",
    "    plt.ylabel(\"Star Formation Rate ($M_\\odot yr^{-1}$)\")\n",
    "    plt.xlabel(\"Star particle age (years)\")\n",
    "    plt.yscale(\"log\")\n",
    "    fig.savefig('Logged SFRs/galaxy_%i.png' %i)\n",
    "    plt.tight_layout()\n",
    "    fig.clf()\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 538,
   "metadata": {},
   "outputs": [],
   "source": [
    "def smoothing(galaxy, resolution):\n",
    "    \n",
    "    \"\"\"\n",
    "    Smooths star formation rate\n",
    "    Parameters:\n",
    "        galaxy: galaxy number (0-585)\n",
    "        resolution: standard deviation of Gaussian used for smoothing\n",
    "    Returns the smoothed array and the (n,bins) output of the histogram   \n",
    "    \"\"\"\n",
    "    sfr_array, n, bins = sfr(galaxy)\n",
    "    smoothed_sfr = convolve(sfr_array, Gaussian1DKernel(resolution))\n",
    "    return np.array(smoothed_sfr), n, bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 539,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extrema(galaxy, resolution):\n",
    "    \n",
    "    \"\"\"\n",
    "    Finds peaks of smoothed SFR, requiring > 10 Myr separation between neighboring peaks\n",
    "    \n",
    "    Parameters:\n",
    "        galaxy: galaxy number (0-585)\n",
    "        resolution: standard deviation of Gaussian used for smoothing\n",
    "    \n",
    "    Returns a numpy array of indices of relevant maxima\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    smoothed_sfr, n, bins = smoothing(galaxy, resolution)\n",
    "    maxima = np.array(signal.find_peaks(smoothed_sfr)[0])    \n",
    "\n",
    "    # Checks 10 Myr separation between peaks\n",
    "    maxima_separation = np.diff(maxima) - 10\n",
    "    ignored_peaks = np.where(maxima_separation <= 0)\n",
    "    final_maxima = np.delete(maxima, ignored_peaks)\n",
    "    \n",
    "    return np.array(final_maxima)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 549,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sfr_end(galaxy, resolution):\n",
    "    \n",
    "    \"\"\"\n",
    "    Find the indices of the ends of starbursts, defined as the point at which SFR is a tenth of peak value\n",
    "    \n",
    "    Parameters:\n",
    "        galaxy: galaxy number (0-585)\n",
    "        resolution: standard deviation of Gaussian used for smoothing\n",
    "        \n",
    "    Returns maxima, indices of end times, and array of peak widths that were used for the search window\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    \n",
    "    smoothed_sfr, n, bins = smoothing(galaxy, resolution)\n",
    "    maxima = extrema(galaxy, resolution)\n",
    "    t = 0.1 * smoothed_sfr[maxima] # points that are searched for \n",
    "    idx = np.zeros(len(maxima), dtype = int) # array to hold indices for ends of star formation\n",
    "    \n",
    "    # Finds widths associated with each peak, to define a window for search\n",
    "    widths = np.array(signal.peak_widths(smoothed_sfr, maxima, rel_height= 1.0)[0], dtype = int)\n",
    "    \n",
    "    # Condition to ensure that first peak is identified\n",
    "    if ((maxima[0] < 20) & (widths[0] <= 65)):\n",
    "            widths[0] = 65\n",
    "            \n",
    "    # Defines a search window and populates array of indices   \n",
    "    for i in range(len(maxima)):  \n",
    "        if (widths[i] == 0):\n",
    "            widths[i] = 50\n",
    "        search_array = np.asarray(smoothed_sfr[maxima[i]:maxima[i] + widths[i]])\n",
    "        idx[i] = (np.abs(search_array - t[i])).argmin()\n",
    "\n",
    "    return maxima, maxima + idx, widths\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 562,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_tend(galaxy, resolution):\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    Plots end times for SFRs as vertical black lines\n",
    "    \n",
    "    Parameters:\n",
    "        galaxy: galaxy number (0-585)\n",
    "        resolution: standard deviation of Gaussian used for smoothing\n",
    "    \n",
    "    \"\"\"      \n",
    "    smoothed_sfr, n, bins = smoothing(galaxy, resolution)\n",
    "    maxima = extrema(galaxy, resolution)  \n",
    "    maxima, times, widths = sfr_end(galaxy, resolution)\n",
    "    \n",
    "    \n",
    "    # An array to hold height ratios to check for shallow peaks\n",
    "    height_differences = np.zeros(len(times))\n",
    "    \n",
    "    # Removes maxima where identified end times are too shallow\n",
    "    for i in range(len(times)):\n",
    "        valley = smoothed_sfr[times[i]]\n",
    "        if (valley == 0):\n",
    "            valley = 1e-10\n",
    "        height_differences[i] = smoothed_sfr[maxima[i]] / valley       \n",
    "    ignore_peaks = np.where(np.ceil(height_differences) < 9)\n",
    "    maxima = np.delete(maxima, ignore_peaks)\n",
    "    times = np.delete(times, ignore_peaks)\n",
    "    widths = np.delete(widths, ignore_peaks)\n",
    "    \n",
    "    \n",
    "    # Repopulates widths adjusted for the new maxima  \n",
    "    t = 0.1 * smoothed_sfr[maxima]\n",
    "    for i in range(len(widths)-1):\n",
    "        new_width = np.argmin(smoothed_sfr[maxima[i]:maxima[i+1]])\n",
    "        widths[i] = new_width\n",
    "        search_array = np.asarray(smoothed_sfr[maxima[i]:maxima[i] + new_width])\n",
    "        times[i] = (np.abs(search_array - t[i])).argmin() + maxima[i]\n",
    "    \n",
    " \n",
    "    # Removes peaks that are less than a factor of 10 from the subsequent valley\n",
    "    forward_differences = np.zeros(len(times))\n",
    "    for i in range(len(times)):\n",
    "        if (i == len(times)-1):\n",
    "            forward_valley = np.min(smoothed_sfr[maxima[i]:])\n",
    "        else:\n",
    "            forward_valley = np.min(smoothed_sfr[maxima[i]:maxima[i+1]])\n",
    "        \n",
    "        if (forward_valley == 0):\n",
    "            forward_valley = 1e-10   \n",
    "        forward_differences[i] = smoothed_sfr[maxima[i]] / forward_valley      \n",
    "    ignore_peaks = np.where(forward_differences < 10)[0]\n",
    "    maxima = np.delete(maxima, ignore_peaks)\n",
    "    times = np.delete(times, ignore_peaks)\n",
    "    widths = np.delete(widths, ignore_peaks)  \n",
    "     \n",
    "    # Removes peaks that are less than a factor of 4 from the previous valley\n",
    "    backward_differences = np.zeros(len(times))\n",
    "    for i in range(1, len(times)):\n",
    "        if (i == 0):\n",
    "            backward_valley = np.min(smoothed_sfr[:maxima[i]])\n",
    "        else:\n",
    "            backward_valley = np.min(smoothed_sfr[maxima[i-1]:maxima[i]])\n",
    "        if (backward_valley == 0):\n",
    "            backward_valley = 1e-10\n",
    "        backward_differences[i] = smoothed_sfr[maxima[i]] / backward_valley\n",
    "    ignore_peaks = np.where(backward_differences[1:] < 4)[0]\n",
    "    maxima = np.delete(maxima, ignore_peaks + 1)\n",
    "    times = np.delete(times, ignore_peaks + 1)\n",
    "    widths = np.delete(widths, ignore_peaks + 1)  \n",
    "    \n",
    "    \n",
    "    # Plots end times and SFR on a log scale\n",
    "    plt.plot(bins[:-1], smoothed_sfr)\n",
    "    plt.scatter(bins[maxima], smoothed_sfr[maxima], color = 'r')\n",
    "    plt.yscale(\"log\")\n",
    "    [plt.axvline(bins[times[i]], color = 'k') for i in range(len(times))]\n",
    "    plt.title(\"Smoothed Star formation rate of star particles in galaxy %i\" % galaxy)\n",
    "    plt.ylabel(\"Star Formation Rate ($M_\\odot yr^{-1}$)\")\n",
    "    plt.xlabel(\"Star particle age (years)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 559,
   "metadata": {},
   "outputs": [],
   "source": [
    "def t_end(galaxy, resolution = 5):\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    Finds the end times for star formation peaks\n",
    "    \n",
    "    Parameters:\n",
    "        galaxy: galaxy number (0-585)\n",
    "        resolution: standard deviation of Gaussian used for smoothing\n",
    "        \n",
    "    Returns array of star formation peak values and array of end times in years\n",
    "    \n",
    "    \"\"\"   \n",
    "    \n",
    "    \n",
    "    smoothed_sfr, n, bins = smoothing(galaxy, resolution)\n",
    "    maxima = extrema(galaxy, resolution)  \n",
    "    maxima, times, widths = sfr_end(galaxy, resolution)\n",
    "    \n",
    "    \n",
    "    # An array to hold height ratios to check for shallow peaks\n",
    "    height_differences = np.zeros(len(times))\n",
    "    \n",
    "    # Removes maxima where identified end times too shallow\n",
    "    for i in range(len(times)):\n",
    "        valley = smoothed_sfr[times[i]]\n",
    "        if (valley == 0):\n",
    "            valley = 1e-10\n",
    "        height_differences[i] = smoothed_sfr[maxima[i]] / valley\n",
    "\n",
    "        \n",
    "    ignore_peaks = np.where(np.ceil(height_differences) < 9)\n",
    "    maxima = np.delete(maxima, ignore_peaks)\n",
    "    times = np.delete(times, ignore_peaks)\n",
    "    widths = np.delete(widths, ignore_peaks)\n",
    "    \n",
    "    \n",
    "    # Repopulates widths adjusted for the new maxima  \n",
    "    t = 0.1 * smoothed_sfr[maxima]\n",
    "    for i in range(len(widths)-1):\n",
    "        new_width = np.argmin(smoothed_sfr[maxima[i]:maxima[i+1]])\n",
    "        widths[i] = new_width\n",
    "        search_array = np.asarray(smoothed_sfr[maxima[i]:maxima[i] + new_width])\n",
    "        times[i] = (np.abs(search_array - t[i])).argmin() + maxima[i]\n",
    "    \n",
    " \n",
    "    # Removes peaks that are less than a factor of 10 from the subsequent valley\n",
    "    forward_differences = np.zeros(len(times))\n",
    "    for i in range(len(times)):\n",
    "        if (i == len(times)-1):\n",
    "            forward_valley = np.min(smoothed_sfr[maxima[i]:])\n",
    "        else:\n",
    "            forward_valley = np.min(smoothed_sfr[maxima[i]:maxima[i+1]])\n",
    "        \n",
    "        if (forward_valley == 0):\n",
    "            forward_valley = 1e-10\n",
    "    \n",
    "        forward_differences[i] = smoothed_sfr[maxima[i]] / forward_valley\n",
    "        \n",
    "    ignore_peaks = np.where(forward_differences < 10)[0]\n",
    "    maxima = np.delete(maxima, ignore_peaks)\n",
    "    times = np.delete(times, ignore_peaks)\n",
    "    widths = np.delete(widths, ignore_peaks)  \n",
    "    \n",
    "    \n",
    "    \n",
    "    # Removes peaks that are less than a factor of 4 from the previous valley\n",
    "    backward_differences = np.zeros(len(times))\n",
    "    for i in range(1, len(times)):\n",
    "        if (i == 0):\n",
    "            backward_valley = np.min(smoothed_sfr[:maxima[i]])\n",
    "        else:\n",
    "            backward_valley = np.min(smoothed_sfr[maxima[i-1]:maxima[i]])\n",
    "        if (backward_valley == 0):\n",
    "            backward_valley = 1e-10\n",
    "        backward_differences[i] = smoothed_sfr[maxima[i]] / backward_valley\n",
    "    ignore_peaks = np.where(backward_differences[1:] < 4)[0]\n",
    "    maxima = np.delete(maxima, ignore_peaks + 1)\n",
    "    times = np.delete(times, ignore_peaks + 1)\n",
    "    widths = np.delete(widths, ignore_peaks + 1)  \n",
    "    \n",
    "    return bins[times], smoothed_sfr[maxima]\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 558,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.00000000e+00 2.08783456e+01 4.09486242e+03 9.44498987e+02\n",
      " 4.68820509e+04 1.73830797e+02 4.39265500e+02]\n",
      "[0.00000000e+00 2.08783456e+01 4.09486242e+03 9.44498987e+02\n",
      " 4.68820509e+04 1.73830797e+02 4.39265500e+02]\n",
      "[ 98  42  16  64  23  27  95 227  27 193  45  45]\n",
      "[ 54 114 134 161 205 244 318 448 387 510 484 513]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEWCAYAAABSaiGHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9d5hkZZX4/zmVOqeZnunJkQkMUUEyCAgKSjCsqxhZUFZ2Dbvrz11Xdl1Mq2zQNbusAVEEcb8q6IIoSgbJSBoGJk9P6Omcq6uq6/z+uPd2V3dXuLfqVldV9/t5nnq669at9566des994T3HFFVDAaDwWAoNoFSC2AwGAyG+YFROAaDwWCYFYzCMRgMBsOsYBSOwWAwGGYFo3AMBoPBMCsYhWMwGAyGWcEonCIhIpeLyIM+jbVGRFREQn6M5wcicrWIdIjIkIgsLLU80xGRO0Xk/aWWo9iIyCYReVpEBkXko6WWp9iIyLtF5Lcu9rtWRH7s87E/JSLf9XNMl8f1/bOUiopXOCJyhog8LCL9ItIjIg+JyGtmWYaSKoRs58BPxZdyvDDwZeD1qlqvqt1+jp+HPDN+kKp6oar+sFQyOczCtfH3wL2q2qCqX/Mgl+/Xhd+kO3eqepOqvr4U8qjqv6rqB0px7NlGRI4WkbtEpEtEdNprVSLyPRHZY9/oPC0iF7oZt6IVjog0Ar8Gvg4sAJYDnwHGSinXbFLsc5BhomwDqoEX8hhPRKSir7tURCRYYhFWk8f3UCjFvrkqJ2t+nhIHbgWuTPNaCNgHvBZoAv4ZuFVE1uQcVVUr9gGcCPRlef1y4CHgK0AfsBM4zd6+DzgMvD9l/ybgRqAT2AP8ExCwXwvYz/fY77sRaLJf2wsoMGQ/TrWP8SDwH0AvsAu4cNqxvgccBPYDnweC9mtB+31dtsx/bY8f8nIOgCOBKDBuy9Vnb38T8DQwYJ+Ha1Pes8Y+1pX257p/2pgbgeGUz/sHe/tpwONAv/33tJT33At8wf4uRoEj7G2fBx62x/kVsBC4yZbrcWBNyhhftWUdAJ4EzrS3XwDEsH4gQ8CfUo75ARffnfN5329/3i7gmizX1A3At4E77PNwXo7zOePasLdfAWy1r427gNVZjnkJllLpsz/Xkfb2P9jfbdQee2OG38BOYBDrGnx3Ma4Le5+zgXbgU/Z53A28O+V1T+OnO3f253kw5X1HAb8DeoAO4FP29muBH6fsdwrWtdYH/Ak4O9s5yvA9TIyZx3WzEOsad67tz0/7HGmv7wyf5WfAIazf2v3AUfb2CPAM8JGUeeQh4NPAEmAEWJgyzglYc104i9xHAOpiLn4WeFvO/XLtUM4PoBHoBn4IXAi0pPmxJYC/sE/+5+2L45tAFfB6+yKrt/e/EbgNaLAvqJeBK1MmiO3AOqAe+Dnwo2kXX2jasePAB+1jXw0cAMR+/ZfAfwN1wGLgMeAv7dc+BLwErMSyWu6ZPr7Hc/DgtG1nA8dgTcTHYv1Q3zzts9xoy1aT5phTPq8tYy/wXqy7n8vs5wvt1++1z/tR9uthe9t2YD2W8n3RPt/n2fvcCPwg5ZjvwfrRhoCPY/3gqtP9IFOO+QEP393/ADXAcVjW4ZEZrrkbsH7op9vnr9rl+Uy9Nt5sy3Ok/Xn+CXg4w/EcBX++fd7+3n5vZPrnTPPeOqwJbJP9fCmTk1MxrouzsX5vX8b6fb3Wln1TPuNnOHcTcmP9Tg/a10O1/fzk6dcEltXfDbzRPvb59vNF2c5Rms+XOqYjm9vr5hb7UQtswVIuqQrH9fWNdT032Of4v4BnUl47Guu3dyRwDfBHJm9k7wCuTtn3K8DXc8yxORUOlscjCmzOOWd7meDL8WGf2Buw7qwSwO1AW8rF+UrKvsfYF0lbyrZu4HgspTAGbEl57S+x/OMAvwf+KuW1TVgKJZTlh7E95Xmtvc8S+wsaI+VHizVJ32P//wfgQymvvX76+B7PwYOZzp+9z38BX5n2Q1qXZf8pnxdL0Tw2bZ9HgMvt/+8FPjvt9XtJuSME/hO4M+X5xak/pDQy9ALH2f9fS3aF4+a7W5Hy+mPAOzMc9wbgxjzOZ+q1cSf2jYz9PIB197k6zVj/DNw6bd/92Hfo5FY4fcDbmKYginRdnG1ff3Up224F/jmf8TOcuwm5sX4zT2cYe+KaAP4B+wYj5fW7sKyTjOcox5iurxusuSWOrdTsbVMsHK/Xd8p+zbYcTSnbPo51w9oLbEjZ/g7goRSZDgEn5fjMWRUO1k3Q3cB/ZxvHeVS8L11Vt6rq5aq6Aku7L8O6kB06Uv4ftd8zfVs90Iplku5JeW0P1t0R9rjTXwthKY9MHEqRc8T+tx7L7x4GDopIn4j0YVk7i1OOtW/asTLi4hxMQUROFpF7RKRTRPqxLKrWabvtS/PWTEw/N47My1Oepxtv+veQ7ntxZP64iGy1EyP6sKyi6TK7lS/dd3co5f+R1GOnYcpncXk+U1kNfDXlu+8BhKnnK63sqpq0j59u3ymo6jDWJPMhrGvt/0Rkc6b9fboueu3jOuyxP0MxrruVwA4X+60G3u6cb/ucnwEs9XqO0uDmulnEZNzDYfo15Or6FpGgiHxJRHaIyACW25Jp+/4QSyHeoaqvpGy/DdgiIuuwrLx+VX3MxWdMix2L/RGWS/vDbt5T8QonFVV9CesO9Og83t6FdReyOmXbKqy7SbDcYdNfS2BNkurxWPuwLJxWVW22H42qepT9+kGsH1PqsVyR5hykk+0nWFbQSlVtAr6DNeFNGcrtMZl5bmDqufM63hRE5Eysu9Q/x3IZNmO5tRyZc42d7bvLh+nHy3Y+08m2D8t92pzyqFHVh3PJLiKCdW3sT7PvTEFV71LV87FcRS9huYAyyeXHddEiInUpz1fZnyGf8XMdax+WSzYX+7AsnNTzXaeqX4Ks58gvOrGutxUp2yZ+3y6u71TeBVyK5XpuwlIsTNv3W1iJRG8QkTOcjaoaxbI4343llfhRvh/Ivg6/h3XT9jZVjbt5X0UrHBHZbN8ZrLCfr8Qys//odSxVHcf6Mr4gIg0ishr4O8BJt70Z+FsRWSsi9cC/Aj9V1QTWBZXEihG4OdZB4LfAf4pIo4gERGS9iLzW3uVW4KMiskJEWoBPZhrLxTnoAFaISCTlbQ1Aj6pGReQkrIu4EO4ANorIu0QkJCLvwPJT/7rAcR0asH6wnUBIRD6NFbty6ADWZMl+y/bd+SVfpvOZ7tr4DvCPInIUgIg0icjbM4x9K/AmEXmdnY7+cayblXTKaQoi0iYil9gKYAwr8D5uv1zM6+IzIhKxJ9KLsILc+Yyf63f1a2CJiPyNnarbICInp9nvx8DFIvIG20KoFpGz7d9XtnPkC/bc8nPgWhGptS2o96Xskuv6Ztq+Y1ihgFqsa3kCEXkvVjLA5cBHgR/a17zDjfZrlzA5t83AziatxvL6YJ+zqpRdvo3lyr9YVUczfvhpVLTCwQr4nww8KiLDWJPs81g/ynz4CFaQcydWhtlPgO/br30f647gfqxMlqi9v+Mu+wLwkG2yn+LiWO/D+jJfxPK1/i/WHRZYd1h3YWXTPIV1sWYi1zn4A1aG0yER6bK3/RXwWREZxMpgudWFvBlRax3ORfYxu7EC2xepalfWN7rnLqy4x8tYLpooU10SzoTWLSJPpXl/xu/OJzKez3TXhqr+ArgOuMV2izyPlfAxA1XdhhVQ/jqWFX4x1o885kKuANZ3cgDLbfdaW1Yo3nVxCOt6PoCVcfgh2+r2PH6u35WqDmK5hi62j/sKcE6acfZhWQWfwprU9wGfwDo/2c6Rn3wYyyI5hHUt3szk0oVc13cqN9r77MeaOyZurkVkFZYr/X2qOqSqPwGewEoOAEBVH8JS4k+p6u4s8q7Gcms7KfejwDb7OKux4tvHY10/Q/bj3blOgpMxZTAYDAUhImdjBbdX5Np3viMi1wFLVPX9JTj2H4CfqOqsV02odAvHYDAYyh7b9X2s7ao6CWu90S9KIMdrgFcDP53tY4OVOWEwGAyG4tKA5UZbhrX4+D+xssZmDRH5IdYasI/Z7shZx7jUDAaDwTArGJeawWAwGGaFeeVSa21t1TVr1pRaDLZt2wbApk2bymqsSjp2uVGJ56JYMvs9biWe27nGk08+2aWqiwodZ14pnDVr1vDEE0+UWgzOPvtsAO69996yGquSjl1uVOK5KJbMfo9bied2riEiWauduMW41AwGg8EwK8wLhSMiF4vI9f39/aUWxWAwGOYt80LhqOqvVPWqpqamUotiMBgM85Z5oXAMBoPBUHrmhcIxLjWDwWAoPfNC4RiXmiEvbroJ1qyBQMD6e9NNpZbIYKho5lVatMHgmptugquughG7b96ePdZzgHfnLIprMBjSMC8sHIPBM9dcAyMjHGho5cHVxxELhCzlc801pZbMYKhYjIVjMKRj714eX76F977js0TD1ZzQ/iI3/fSfqN67t9SSGQwVy7ywcEzSgMEryVWr+eSFH6FtqIdP3309T67YwrdP/jNY5brbt8FgmMa8UDiznjRQjGCzCWDPKg988kvsWLiSv33wJ1zx5O286aUHuP7kt9L1mX/N/WaDwZCWeaFwZhUn2LxnD6hOBpsLURDFGNOQlVvqj6A1lOSNo3tBhL/dcQ+j4Wp+uuqkUotmMFQsRuH4zTXXMJhQvnL6u/jaae9kJFxVeLD5mms4LBHuWXci0WDY2mYC2EUjPp7kwVe6OP/Vq4ns3AHJJEc89ygnr13AzY/tJZk0PaQMhnwwSQM+k9y7jyvf+QUeW3UMAE8u38wNP7sWKSDYvK9/jEuu+Ca9tU2ctO95fnLzpwhpEvbutdxrBl95Zl8fg2MJztowtRr7ZSet4m9++gyP7+7h5HULSySdwVC5VKyFIyJvFpH/EZHbROT1pZbH4benXcRjq47hS3d+jWt/9x3uW3civ91wSkHB5i9deDXxYJgPP3wLj608mp8f/TrrBRPALgoPvNJFQOC09a1Ttp+/pY2acJDb/nSgRJIZDJVNSRSOiHxfRA6LyPPTtl8gIttEZLuIfDLbGKr6S1X9IHA58I4iiuuJb51/BWt7D/D25+7mPU/fwfL+Dn58wkXwhS/kNd6BvlF+s+rVvOuFu/n4Az/myI6d/OhVb4Ta2rzHNGTnT/v62NjWQFNteMr2uqoQ529p447nDhJLJEskncFQuZTKwrkBuCB1g4gEgW8CFwJbgMtEZIuIHCMiv572WJzy1n+y31dyth8e4tnREO/b0kxw1UpCKG9tf4qHVh/HoYvelteYv3h6P+MI773sbGT1at764j08t3QDB77xP7Oz4n2eZcepKs/v7+eY5ekzGi85bhl9I3Ee3N45y5IZDJVPSRSOqt4P9EzbfBKwXVV3qmoMuAW4VFWfU9WLpj0Oi8V1wJ2q+lSmY4nIVSLyhIg80dlZ3EnizucOAnDh5RfD7t2QTPKW73+JJMJdLxzKa8yHtnexeUkDK//iMti9mzNvu8HafsyZPkmdhXmYHXegP0r3cIxjVqRXOGdtXERTTZjbnjFuNYPBK+UUw1kO7Et53m5vy8RHgPOAPxORD2XaSVWvBz4DPBWJRPyQMyO/fbGDV69qZklT9cS2dYvqWdFSwyM7uj2PF42P88SeXk4/YjKWsKmtgdb6CA9t7/JF5qxccw3JkVF+dvTr+Mapf05/Vd2cz457rt1aHHx0BgsnEgrwxmOW8rsXOxiJJWZTtPyZZ1aqoXwpJ4UjabZlzD9V1a+p6gmq+iFV/U62gWdj4edANM4LB/o5Y1pmE8Cp6xbyx13dntNpn9jdSyyR5PQjJjOiRITT1rfy6K7pBmIR2LuXr5z5bj7xpr/lP856H+99x+eIB4JWdlwF0d47wp99+2FO/9IfuO2Z/Vn33XZoEBE4ckljxn0uPX4ZI7Fx7t562G9R/WceWqmG8qWcFE47sDLl+QrAF7/FbJS2eWJ3D0mFU9YumPHaqesX0jcSZ+uhAU9jPrSji1BAOGnt1BTcY5Y3cbA/Ss9wrCCZc7Frywl85+S38dbn/8DXb7uOZ5du5OdHn1tR2XHR+DgfvPFJth0apLk2zMdv/RMvHsj8PezoHGJ5cw01kWDGfU5as4AljdXcnkN5lQXXXMP/HPUGLrz8a/zo+AutbXPcSjWUL+WkcB4HNojIWhGJAO8Ebi+xTK55dGcP4aDwqlUtM157tb3Ncde45eHtXRy/spn6qqnLpbYss+6+tx70psCmEx9PMjyW2S303ff9I8Fkkk/e+30ueukBNnXu5mfHvaGg7Lgn9/Tyuxc7GJ+lxZM/fXwfWw8O8OV3HM9PPnAKdVUhvvr7lzPuv6NziPWL6rOOGQgIlxy/jPte7qRvpLhKv1B+H27jC+deSWd9C//8hr/mkZXW+rBKs1INc4NSpUXfDDwCbBKRdhG5UlUTwIeBu4CtwK2q+oIfx5sNl9qju3o4bkVz2jvjVQtqqa8K8UKWO+vp9I/EeW5/P6cd0TrjtSOXWgon2516Ltp7Rzjr3+7hVZ/7HXe/2DHj9Wh8nF+N1HHBImFxaxMiwiUHn+WJZZvZ/6a35nXMXz69n7d9+2E+eOMT/NMvn8/9hgJJjCe5/v6dnLi6hfOOXExTbZh3vGYlf3jpMP0j8Rn7J5PKzs7hnAoHrGy1+Lhy5/P5JYPMFtef+U5W9R7knuuvYulAJ984zV5BUEFWqmHuUKostctUdamqhlV1hap+z95+h6puVNX1qurbIpNiu9RGYgme29/PyetmutPAuiPesrSRFw64P/4jO7tJKpyRRuEsqIuwpLGaFwuwcL54x0v0j8ZZ0ljNJ3/+HNH4+JTX793WyUA0wdvedsZExt2bbvwyAL/LI+MuPp7kC3ds5VWrmnnvKau5+bG9PL23Ny/ZB6NxV5bFY7t62N83yhVnrEXEChG+6ZilxMeV37448zMcHIgyGh9n/eK6nGMftayRdYvqcsaESsmurmEeXbyRd2y9h4bYKO/80108tOZ4OhYtN2u4DCWhnFxqRaPYFs7z+wcYT+qE6ywdW5Y1svXgoGtX0sM7uqgJBzl+ZXPa1zcuaeDljsG85O0cHOM3Lxzivaes5gtvOZquoTHu3jrVyrnv5cM0VIU4NaWEy5rWOpY1VfP4Hu+K4p6XDtM5OMZfnX0En7xwM/VVIX78R+9une8+sJMTPnc3r/rc7/jSnS+hmvl83vH8QWrCQc7ZNLls69gVTSxvruF3aay6HYeHAFjXmtvCEREuPnYZj+7qoXtozPPnmA3usNP033bFm2D1as7f/igA9/3zV0zXUkNJmBcKp9gWznP7rXEzrd0AS+GMxsfZ2zPiasyHtndx0toFRELpv6K1C2vZ2+1urOn8fqsVQ3nzq5Zz2vpWljZV8/+ebJ94XVW5/+UuTjtiIaHg1OOfsGYBT+7uzTrRp+PnT+1nUUMV52xaRF1ViEuOX8avnz3AQHSmaysTt//pAJ//v62ctXERbzl+Od+5bwc/fjS90komlbte6OCczYumuDlFhFPWLeTJPTM/g/PdrGmtdSXPeUe2oQr3vVyei0D/uLObTW0NLLn8XbB7N0ce2kFbYxX3LlhXatEM85R5oXCKbeE8197HksZqFjdUZ9zHiQvs7BzKOd7+vlF2dA6ndac5rF5Yx+BYgvi49xIrf9zZzaKGKjYvaSAYEC45bhkPvNI1Mfnv7Bpmf98oZ22cmeJ94uoWDg1E2d836vp4yaTyyM5uzt20eEKBvfVVyxlLJHngZXfriUZj43z2Vy9y/MpmvvXuV/Mfbz+OMze0ct2dL9GVxsJ4cm8vnYNjXHD00pmfYU0L3cMxdnUNT9m+r3eESDBAW5bvMZWjljXSWl/FH14qv/ToWCLJE7t7OSXFzSsivHbjIh7a3u35hsFg8IN5oXCKzXP7+zMuFHRYv8iKC0yf5NLhTGDnbF6ccR/nLjwa965wntrbxwmrWibiGq87so1EUnnwFWvyv9++Y59eLRkmM+6e9ZBxt61jkP7ROCelpIwfv7KZppow92xzN1nf9OgeuobG+NQbjyQSChAICNdechTR+DhfvfuVGfvf+dwhIqEA56Y5hyeutj7DE7unugbbe0ZZ3lJDIJBuSdhMAgFnAp+FRbgeeW5/H6PxcU6ZVtX6+JUt9I/G2dfj/obBYPCLeaFwiulSGxpLsLNrmGOzuNMAmmsjtNSG2dHpQuFs7WD1wtoJJZWONQut16YH+3PROTjG3p4RXr16Mjb06lXNNFaHJhTdfS93sra1jpULZrqWjlhcjwie4keP7rSqLKQmVYSCAc7auIh7t3XmXBA7GhvnO/ft5PQjFk5RWusX1fOO16zklsf30t476V5UVX7z/EHO2tA6I6XceV9DdYhn2vumbN/XO8KKlhrXnwvgpLUt9I7EGfX4PRQbJwV/epq+UyPOcQMbDLPJvFA4xXSpvbC/H1UyFntMZd2i+pwutcMDUR7c3sXrt7RNWCDpWNFSS0BgLOFtonMyw1ITHCYn/8MMjSV4eEd3WssAoCYSZPWCWk8K5+l9fbQ1VrGiZaoCO3vjIrqGxnIuiHWsm785b+OM1z587hGICF///fYpxzvQH03rTgPLMjlySSPbDk39DPt6RtIq2WycsNpSgIPR8ipzs/XgIAvqIrQ1Vk3ZvnFJPZFggGf392V4p8FQPOaFwikmztqao5ZnLoXisK61LqdL7cZH9pBIKu85ZXXW/SKhAMuaazy71JxJ1lk86nDu5sV0DcX4xh+2E0skeV0Wd97GtgZe7sgdi0o9prN2KJUzN1gxqgdeyeySsqybHZxxRCuvWTMz7XxpUw3vPnkV//tUO7vtc/u/T7ZTHQ7whqPaMo67aUkDLx8anIhlDI0l6B2Js7LFm8JZv6iOltowgx6SH2aDFw8OsGVp44yblqpQkE1LGni+nCwcU+tt3mAUToG83DHIwrpI1oQBhzWtdRweHMuYGj0aG+fHj+7h9VvaWL0w91qQZU01nvuybLdLt9RGprqazrObi33nvh201IZ5TZoSPQ6bljSwq2uYpIvAc3w8yY7OITYtaZjx2uLGajYvaeCBVzJneVnWTYyPnbch4z5Xn72ecFD40p0vMTSW4Fd/OsAbj1lKQ3U443s2LWlgcCwxkfywz85Q8+pSExGOW9nM8FhxXWrjSeUPL3XwkovySInxJNs6Bjly6cxzDtYNw/bD7m8Yiopd663/UBftDa2m1tscZ14onGLGcLZ1DLKhLfe6DbAqDgCMZVAS//tUO30jcT5wpru01aXN1Yx5zFLbfniI9YtnyttYHeZj520gGBD+/oLNhIOZL40jFtcznlTGXFhXOzuHiY8rm9MoHLCsnMd39zIamzlh94/E+eY92zNaNw6LG6r56Os28JsXDnHhV+9nMJrgfaeuySqXMxm/dNCy+ByF49WlZo1lpbwXM/HrUz9/jitueII3fe3BtGuIUtndPUwskWRzhgKk6xbV0TGQ+cZnVrnmGp5qXsnpV/+AM67+AV85/V2m1tscZl4onGLFcFSVVzqG2NSWfjKdzsoJhTNzck0mle8/uIvjVjZPZFHlYqlHCyeZVHZ0DnFEhtItH3rtel787Bu47KTsZU8c6yvqIn7k3JFvaks/+Z25YRGxRJLHdk+tfq2qfO7/XqR/NM6n3nhkzuNc/dr1XH7aGqLxJJ964+aMC2YdNtrf2TY7FrWv17J0Vnq0cMBSOKpatMSB59r7+ekT+3jXyavY1NbAp297Puv3vtNOTEl3YwGTKfpeE04A391f4/va+cQb/4bm0QHe+NKDfPWMd/HM0o2m1tscJS+FIyJ1dofOec2B/ihDYwk2Zrh7n86EhZPGMnh4Rze7uoa54vQ1WZMFUlnaVI2qul6Lc6B/lGg8yREZJiKwfPy5yPY5prPj8BABIWO5GGdx6wMpiye7h8b4+M/+xP8+2c5fn3PEjHhTOkSsNOnHrzmPq85an3P/huowS5uqeaVj0sKpjQRZUOe9Z9IW21oqVn+cmx/fS004yD9euJlPXLCJg/3RrCV1dndbCmdtBresk/3oWUEWodXBvSdfwI6FK/mH+37Iv935VZpGB/nGqX9uar3NUVwpHBEJiMi7ROT/ROQw8BJwUEReEJF/F5HMDvY5zMt2AH6jSwunpTZMfVUorWVwy+N7aa4N84ajlrg+/lK70VsmF910nISFdVnSrd2Q7XNMZ0/PCEubajIqsupwkFPXLeTO5w+RGE9y4yO7Ofvf7+X2Zw7wkXOP4O/On5mZ5hcb2hp4xY5l7OkeZvXCOtfKPpU1C+sIiBQljqOq3P1iB+ceuZiG6jBnb1zEhsX13PrEvozv2dU1zIK6CE216WNYqxbWEgyId4VzzTX8qXE5F73/vzjzL7/LLce+vmD3152XXkFjdIgLXn6Y+tgo73/q19y94RQOXvuveY9pKF/cWjj3AOuBfwSWqOpKVV0MnAn8EfiSiLynSDIWTLFiOE5q8MbF7hSOiLCipWaGZTA8luC3L3Zw6XHLqA67NxyXNVvuH7dutXbHbZRHnCIVEWHlglpXFs7enhFWL8x+vPecspr9faO89t/v5dO3vcDxq5r5zd+cxcdfvykvBeCWjYvr2X54iPGksqd7hDU55MxEKBigOhzMz0WVg11dwxweHOP09VZGn4jwpmOX8sSeXg4PRjO+J9tnqQoFWdniPcOxq6ufy99+Lb01jbQNdfPJCz/KT457Q97ur2RSuTfewFltEcIrV4AIl/ZZrSPu2nBqXmMayhu3Cuc8Vf2cqj6rqhNXqar2qOr/U9W3AT8tjoiFU6wYzraOQdoaqzLeSaZj1YLaGTGcB7d3EUskecPR7q0bmLRw3CucEUIBoa2hKvfOOVi1oMaVhbO3e2TCBZeJ845czNVnr6cqFOBfLt7CjVeclNXt5xcb2uoZSyTZ3T3M3p4R1rTmb/lVhwNFUTh/3GnFtlJL1Fx49FJUyZg8sKtrmLU5CpCuaKllzKO8PzjnPfRX13PDz/6Fm2/+FGfueorPn/sBDmw+1tM4Di8cGKBrKMbrLnjNREXy9c89xsa2+rJv+2DID1cKR1UnFhmIyD/k2me+sKNz2PPEuHJB7Yw7y3tesiozZxYW1HUAACAASURBVMvESseCuggBEWIuYzj7e0dZ2lw9oyBnPqxeWJfTwhkaS9A9HGNVDstBRPiHCzbzh//vbP7i9LVFtWpS2WC7Qu/d1kkiqRljHm6oDgeJJpIk8qhtl43Hdll179amKMONbfUsb65JW1JnNDZOx8BYTmttRUuNa1csWIkutxz7es7d/RQbuvcRTo7zr7/5BuOBIF+//NPuP1AKT+5xlOnU8juvO7KNJ/ekz1w0VDY5Zx4RuTXl8TPgA7MgV0Wwp3t4osSMW1YtqCU5LdD/wCtdnLGhNWsqcjpEhHBQXCcNtPeOsrzZexZWOpY1Vc/4HNPZYwevVy8oLGZULDbYNwu/s3vjFGLh1ISDqKqnoqZu2HpwkGOXN01RwiLCyWsX8NiunhlFOJ3jr1iQ/Xte3lxDfDzpai0VWNXLu8cDvPvcI2H1ahBhZUsNlyxSbhuqzWvh67Pt/SxqqGJJ49Q1bCetXUAiqXn3SzKUL25muAFV/XP78Xbg7mILVQn0jcToG4l7Vjgr7YnAsXIO9o+yv2/Us3XjEA4GiI+7mzTae0dnlJfJlyVNueNHztqWXC61UuFkqjluq3xjOMBE7M1NcVa3xBKZF82evG4BXUMxdkwrleQonOXNOSwc+zp0a+Xc/3IX1eEAp175ZxPuL3bv5j3vOpeR2Di/fvagq3FSeXZ/P8etaJph0Z6wuoWAWF10DXMLNwpnemvAsliRJSJHish3ROR/ReTq2T7+nm6nd4p3Cwcm1+I8tceqafVql2tvphMOBlzFcGKJJB2DUc8r6TMxET/KYuEc7LeC2kub3ZX7LwXHrbDW6yxuqGJRAbGt6rD1U9rto8LZ2TVEIqlpFc5Jay031OPTKl7vtxNDluf4nh2F5FrhvNLJyWsXzkhqcRra/X6rtxYNQ2MJdnQOcczymeulGqvDHLm0kcd3G4Uz18ipcFR1F4CItNrPC74KROT7InJYRJ6ftv0CEdkmIttF5JM55Nqqqh8C/hw4sVCZvOKsdfB6V+xYGE7846m9vVSFAmxJU2vMDeFQwJVLrWMgiupkZluhuElYODQQJRwUFtR6X9syW7zhaKve2lkbFxUUOwoHAwREJjIB/cCpe5dO4axZWEtDdWhGTbT9fSMEXSSGODceMReZah0DUXZ2Dk/UvktFRDh382Ie2t7lqZDstkODqFo9hdJx7IpmXjgwYPr2zDG8BA2+7+NxbwAuSN1gLyT9JnAhsAW4TES2iMgxIvLraY/F9nsuAR4Efu+jbK7Y3TWCiPcU4+pwkEgwMJHh9dTeXo5d0ZSxs2cuIkEhMZ7MWabESaFta/TH2lhYX4WIZFU4Hf1RFjdUu+4vUwouPW45v/zr0/nMJUcVPFZVKOBrDGfboUFCAUnb8lpEOHpZ00yF0zvKksbciSFtjdWIiCsl4fQ+ylS94ZzNixiNj8/oL5QNp2p6pmoIRy1rpH807ntMzFBavMxyvs0aqno/MN1SOgnYrqo7VTUG3AJcqqrPqepF0x6H7XFuV9XTgFlv0L6ne5iljdWe1s04VIWDjCWSROPjPL+/f0qrAK+EgwEU6BmOZd3v8IDVFXNRfeEp0QDBgBAJBrK61A4NRFnSVL7uNLBaFRy/spm6NH1zvBLxWeHs6bb682S6GTl6eSNbDw1OsXD3943mdKeB9f2Fg+4yHJ/f309AZlYYd3BaNHgJ8u/oHCYclIylhJxjvXggd7FSQ+XgReEU27ZdDqQun263t6VFRM4Wka+JyH8Dd2TZ7yoReUJEnujs9K/3/G57ZXo+VIUCjMWTPL+/n/i4zmiS5QUns61zcGab5VQO268vbvRH4QBEQjksnIGxGRlIc5mqcJADPiqcvT0jrMpyjR29vIlYIskrKa0i9veOssKl2zTiMv73/P5+1i+qn1Fh3KGpJsy6RXU8s899j52dnUOsWViX0RI7ckkjAZls/2GYG5TEwvEwfkYlp6r3qupHVfUvVfWbWfa7HvgM8FQk4l8sYXf3yESbZ69Uh4PEEuM8vMPqhPmaNYUrnK6h7Aqnc3CMYMDfeEokGMw4YakqHQNR31x4lUBVMEDXUMy3BaB7uodZlSW92Yn7beuwJuX4eJJDA1FXFg5YFpmbDMfnD+RuoX78ymae2dfnOuayo3Moa4mlmkiQNa11M5rkGSobLwrnH4smhUU7sDLl+QrgQJGPmRdDYwl6hmOsynN9SVXIcoP9/Kl21i+qY2EBbq5w0NLTuS2cKK31EV/jKeFQZpfM4FiCkdg4S5r8s6jKnSo7U80PK6d/JM5ANJF1DdOa1jrCQeEle1I+1B8lqbheaxXO4RIF6B+N0zEwljZxIZVXrWymayjGgf705XZSSYwn2dszwroMVcsd1rXWs7OrTPr2GHzBi8J5UUQ+VTRJ4HFgg4isFZEI8E7gdj8G9ru0jTOhuL2TnI4TL9jdPZL3+hsHxyXRO5IjhjM4VlDabzrCwQDjSU17R9/R72+SQiXgxFr8iOPs6bGyILNVaQgHA6xfVD9RRHa/x+syEgqQGE9mtcgmgvs5lMNm29py03r8YH+U+LjmzPBcv6iO3d0jRfflG2YP1wrHrqF2nh8HFZGbgUeATSLSLiJXqmoC+DBwF7AVuFVVX/DpeL4W75xcXJffZFobmUw08Fo/bTqhgCAiORVO5+CYq66kXnDced1pEhYODVgKZ17FcOyK2Af7ct/l52Kvy0Wzm5ZMtvt2UrLdLu6N2N+fk1CSDqe3Tq4K406Ppe0uWo9PVEPIIee6RXXEEklXqduGysBrLu7TIvIvIlJQMS5VvUxVl6pqWFVXqOr37O13qOpGVV2vqtMXnBZyvKJYOIWsaTl2RTOfvfQozt64qGB5QgGhdyR7aZHDg2Ms9t3Csdxz3WniR4dsC6fcs9T8xDkfmao4e2H/hPLIfo1tbGtgf98og9E4e7uHCYh7l5pjkXVkkXdn1xChgORUfC11EVrrq3jlcG4Lx1GMueR0XG7FamxnmH28Ko6VWK6uAyJym4h8TkTeXgS5fMVvC+dA3yjBgBRkMdRGgrzvVPfN1rIRCgp9WSyc8aTSPVQclxqkT1joGJh/LrWACE014YmMwEI42B+loSpEQ3X2SuROt9mXOwbZa/cecrumy/n+DmWJu+zsHGbVglpXdf42LK6f6C+UjfZeaw1brgoU61rzbBRnKFs8KRy7ntqRwGqszK/tWOtnyhr/LZwoSxqrCZbJgsZwIJB1HU738BhJpQgWjqNw0rvUmmvDea1TqmTaGqsmlG0hHOp3t4bJCeZvOzTEHhe9h1KJuMhw3NU17Lph34a2encutd5R2hqqc3aXXVAXoaEq5KmCgaG88aRwRORuETlOVcdU9SlV/aGqfqJYwvlFMWI4flVd9gPLwsnsUptY9Fkkl1q6CatrMObbItNKYnFDNR1ZYiJuOehy0ezy5hrqIkFe7hhkX0/u3kOphIKCkHnRsKqyr2fEdTWNta11DI4lcpZasorI5v79iAjL0zQsNFQuXl1qfw98RUR+ICJLiyFQMfDbwjnYP8qyMipIGcph4Tgp04t8ThoIiBAMCN1pLJyekRgL6sq3hlqxWNxYlTNF3Q2H+kcn6tVlIxAQNrQ18NTeXrqGcvcemk4oGEib9AFWSvRwbNz1zdXKFncFQdv7Rlxn0q1oqfXUt8dQ3nh1qT2lqucCvwZ+YycQlM+t/iwwnlQO9UdZWoYWTqZFd84E6LdLDSy3WjoLp3d4fiqctsZqDg9GCyo6GR9PcnjQfZWGTW0NE/XOvLaCCAclbdIHpGa9uVQ4C6YWpk2HqtLRP+Y6mcRqFGdcanMFz9lmYkW5twHfBj4CvCIi7/VbMD/x06XWNTRGfFx9q7rsB6GAtQBzOEOHxK5ha0JZWO+/ArAmrJl3yL0jcVrmocJZ3FBFfFxzZg1mo3NwDNXJnkO5eNWqyaKaXuvyhYKZrWOvadaOYsrWerx/NE5sPOk64WZFSw3jSSWRozitoTLwGsN5ENgPfAWrztnlwNnASSJyvd/C+YWfLrVC1+AUAyd435th4ugbiVMVCmSshVUIoTQWjqrSOxKjpTZ7htVcxMnKKyRxwFnD5MalBnDeljYiwQAb2+o93wiFg4G0NwxgZZOB+zTruqoQrfWRrBaOE99qc1nTb7Kdh7Fy5gJeZ6APAS/oTH/BR0Rkq08ylTXOor6lLu8+Z4OQHbzvG4mzMk3hgt7hGC1F6kkTtuuHpTIQTTCe1KIds5xptRMlctW2y4bXNUyt9VU8+qnXURPxnhEYDkjGGM7+vlHqIkGaPdw4rGipZV8WC8dZo+TFwgH3jeIM5Y3XGM7zaZSNw5t8kKco+OlS6yjDFfShgPU19mRYi9M7EvM0aXghHBR6hsdIprg8HEtrPsZwnM+cq11ENiY6pXpYNNtSF8krBT0UDNA/Gk+bWba/12p14GWt2MoF2YP8Tsak23iiY7G5qWptKH8KqhggIiEROU5EFqjqTr+E8hs/XWodg1EiwUDRJvB8mLRwMimceNEm/3AwQFKn1nJzFN98jOEstD9zJjeVGw71j1IVCtBUU/xrzEltT+eO7RzyXg5pWVN19pYVjoXj0qXWUhsmIO769hjKn7wVjp0WfQOwDrhGRN7il1DlTOeAtWLfjwoBfhF2LJwMd9VWPKV4CgemLv6csHDmoUutqSZMMCAFWzhLm6pn5RrLtni3a2iMVo+JJm2N1SRVSWRoe3B4YIyGqpDreKKIWFWtjYUzJyjEwvk74BqgXVU/DrzHH5HKm8ODY742MfODYFAQIWNmVN9IvHguNbvaQqqF48gxH11qgYDQUhvOGBdxw6H+2esj5FjH6RRk91DMc+sMJ+6UySLpHBxjkcffTySUu42CoTIoROHUAsPA8T7JUhF0DESLsp6lEATrzjqdWySZVPqKuAgzlCZDzvm/nNyOs8mCugg9wwUkDQxEPcVvCmGy4vdUeUdiVj+jVo8Kx1GUsQyJA/n8fiKhAHFj4cwJClE4XwX+DXhMRP4OuNsfkfzHz6SBw4NjZVmQsqU2krZFwUA0TlKhuUjurYk75GkxnHBQqK/yPw27ErAUTn4WjqpyeGCMttlSOLY7dnrMyXnude2WoyhjiQwutTzaZETsRnGFLKY1lAd5KxxVfRn4BFYM57eq+m3fpPIZv5IGovFx+kfjZWfhADTWhBmIJmZsd9xbxVoTEwrMDDo7adjlFOeaTRbWV+WdNNAzHCM2npy1LMhQ0C5PNM3C6bTTur3Ww3Pq9aVzgakqhwejrtfgOERCwnhSGRqbeX0bKouCbkFVtRv4hU+ylD2TJWLKz8JprA4xMDozhtNb5IyxgFiWTGr8qGeelrVxWFgXyTuGU4rGdS21My2yfC2ccDCQMcg/EE0QjbuvMpA6JljuuFztGgzljSsLR0TOFZHCO4VVOIc9pnTOJpaFk0bh2BNJMRdhNtdOjR8VMyuuElhQF8m4tiUXE32EZrFx3cK6yIwsNWfhqtcYDmQO8nfm+fuJTCicwouiGkqLWwvnbuCwiCSB54Fngefsvy+q6ry4EiYXrZWjhRNmsAQuNbBjFqkxnOEYm5c0Fu145Y6zFqd3JOb5WjnUb11js2rh1IVnrOFyCnrmY6k6MZfpdA5ax/CqxMKhzK3MDZWF2xjOR4EDwNeAz2MV7zwB+DKwpzii5UZE6kTkSRG5aDaO59x9lqWFk8Gl1jcLizCthIXJYxczDbsSWFBnXR/5JA4cGogi4n/vomxM//7AWpfTUB3Kq3pBOChprTvHvetViU0mNsyL+9o5jSuFo6rfAE4HFPgvIA58TFXPUdUlXg8qIt8XkcMi8vy07ReIyDYR2S4in3Qx1D8At3o9fr4cHhwjFJCyXNDYWBNmLJEkOq3IYc9wjFBAaChixlhLikstmbQKd87nGM5EeZs8Egc6+qO01le5aunsF821kRkWjrXoMz+lFwoGiI/rjKwyx0JZ6PHayNUozlA5uL6qVXVUVa/Dqg59BFY69Ml5HvcG4ILUDSISBL4JXAhsAS4TkS0icoyI/HraY7GInAe8CHTkKYNnDg9aVQYCZdJaOpXGakuhTHer9drWRjEzxlrqIhMKx0nDns8xHCfQno8L6NBAdNbr9LXUhmf0U8qnyoBDOBhAVWdkTU6uz/I+bihNkVhD5eHqtldEzgSOBDbbfxcDg8DCfA6qqveLyJppm08Ctjs12UTkFuBSVf0iMMNlJiLnAHVYymlURO5Q1Rl2vIhcBVwFsGrVqnzEnaAcF306NNp1twai8SnumL6RWNHW4DgsqI1MtBbumceFOx0KKeDZMRB13X/GL1pqIySSyuBYgkY7C6x7KMb6RfV5jefUZ+seGptSD65nOEZDVYhIyLv1ZrVRMC61Ssetn+U+4E/AzcDXVHV3EWRZDuxLed4OZLSgVPUaABG5HOhKp2zs/a4Hrgc48cQTC1o51jk45rq/+2zjTBTTLZy+kTjNRS4C2ZwSJC92GnYlYK1Byi/mcGggyolrvDVRKxQn3tY3HJ+4jrqGxjh5XZpeFy6YrF4QY11KbmvPcIwFeVtNhdWnM5QHbm81rgYewmpB8KiIvCgiPxWRfxKRN/skSzqfT04Foao3qOqvsw7sU6UBa5V0uVo41r3D9MSBwbH4hPVTLJyYVu9wnJ7h+JRt85FgQGisDtOfJokjGyOxBH0j8VnvteS4P52bhcR4kt6ROAvr8rvWJxTONIVbSLp8KBgwWWpzAFcWjqr+d+pzEVkBHAscA7wN+KUPsrQDK1Oer8DKjCsLYgnLXVSOKdHAxIK46WtxBqMJ1i8qbomZljrr2KkWznzOUgOrtl2fR4Wzu8vqsLlmYV0xRMpI6vcHk67A1jxvrhyX2vSYS/dQLO8aceGgFNTUzlAeuF34OcX6UNV2Vb1DVa9T1fem2ycPHgc2iMhaEYkA7wRuL3BMwJ/SNs7FPpvpql5wXCEDo1NdakPRRNFrmk3cIQ/H5nXztVSaarxbOLu6hgFY2zq7CseJ8fXZqdGOomjN8zsMZajP1jsSy9vVGg4EGIwmTJuCCsetS+0eEfmIiEyJuotIxK5C8EPg/W4PKiI3A48Am0SkXUSuVNUE8GHgLmArcKuqvuB2zBzHK9il5vx48s3cKTYTLrU0Fk6xy4FMBMlHYvSMxIiEAtTm0e54LtFcm4/CGQJgTevsJw3ApIUzUWUgz5srEccFNmmRqCo9wzHPKdEO2dooGCoHt7e+FwBXADeLyFqgD6jBUli/Bb6iqs+4PaiqXpZh+x3AHW7HmU26hgv7ERabmnCQUEAYTFE40fg4sfEkDdXFtXAmgs4jcXqHYyyYx4U7HRprwuzvG/X0np1dwyxprHbdnMwvmmrCU/opOYoi33U4YPVJSrVwRmLjjCWS+Vs4wQCjtmxLZrHsj8Ff3MZwosC3gG+JSBhoBUZVta+YwvmFqv4K+NWJJ574wXzH6LILd7bmGUgtNiJi1VNLcak5GWvFVjhVoSD1VSF6hmP0DMfndYaaQ1NNmP4MDfEysbtreNatG5hMcnAWf3YN5le4M5VwMDAl5jKRLp9n0sBkIoKxcCoZzwnxqhpX1YOVomz8YmKVdJm61MAub5Ni4Tjl3IutcGCygKdVZWB+JwzAZAzHbQ+X8aTySsdQ3mtfCqWlNjxh4XQNjxEJBQqqThEOypSsskLXZ02s7SmgsZ2h9Mxe/YwS4k8MZ4yacJC6Mm4q1lAdnpIW7bjXGqqKrwCcAp69s7DQtBJorgmTSCojsfSdL6fz4oEBBscSnLQ2v7UvhZJa3qZrMEZrXWFu0dC0hZo9Ba7PChkLZ04wLxSOH1lqVn/38p5IG2tCU8qJOC61+lmwcJwCkE4MZ77jrLB3mzjwx53dAJyyLq/iHQVjWTjWZN49PFZwrDIcDNA7EidhF/HszbOOmkMoIHZqtFE4lcy8UDh+WDidQ2MsLCCIOhs0zrBwZs+l1lIbpmtwjL5RE8OBSYXT5zKO88jObta11pWsfXlLbYTeYScteixvxeDguMAcN53jUivk2miqidA/ahROJeNJ4YjFe0Tk0/bzVSJyUnFE8w+/LJxF5W7hTOuJ47jUGmehS2JLXYT9faOowoJ5vugTvFk4ifEkj+/q4ZT1pbFuYKpLrXsoVlCGGkyuxelLWUwaCshEkdn8ZPSeam4oL7xaON8CTgWctOZBrArPc57u4bG8S33MFg3VoSk/yAmX2izEnVLdaMbCgaZa9wrnBTt+Uyp3GlgW6nBsnLHEuO0+LlDhTLNw/Kha3lQTdm0xGsoTrwrnZFX9ayAKoKq9wJyfXZJJrYgYTkN1mNH4OONJKzNqNmM4zSlKZj63JnBwLJx0TfGmMxG/KVHCAEx+f/t6RoiNJwte4BwKTF2o6UfV8uY8qjcYyguvCidu961RABFZBJR9rYlCYzgD0TiJpBbsZig2dVXW6n4nHXpoLE5NODgrzbxSLZxSxSHKiYkYjouYwyM7u1m/qI7FJTxvTgvy7YetagcFu9SCU11qVuHOwlytxsKpfLzORF8DfgEsFpEvAA8CX/RdKp8pNIbjZMaUv4VjWTKOwrHK2sxOGndLytqbZc1G4dRXhQgGJOcd+UT8poTuNJi0Sl/p8EfhTE8asNqOF/b7aaoNu7IYDeWLp9lIVW8SkSeB12G1E3izqm4timRlxERtqbK3cKyvczhF4cyGOw2mutGKXbutEhARVwU8X+4YYjg2XrL1Nw5OeaLtnZbCKfTmKiBCVSgwxcI5dkX+STtgWTiDYwkS48kJC8pQWXjNUrtOVV9S1W+q6jdUdauIXFcs4cqFycKd5a1wnOQAJ3YzOFb8wp0OTkn9TW0Ns3K8SsBSOIms+zzbbhXsOHZF82yIlBG/LRxnzN6RGKpK70i84Nhe80RX2+zn1FC+eL1NOD/Ntgv9EKScccppVIpLbdLCiReUhuqFmkiQZz59Pj+6suyz5GeNxprJ+mSZeHZ/Pw3VIVaXuJOsowy2dw4hQsHxFrDLHY3EGY2PE0skfXGpgfvFtIbyw9VsJCJXA38FrBORZ1NeagAeLoZgfiIiFwMXH3HEEXm9v2twjICUf/aV41JLjeEsmcVAtClpMxWrgGcOhdPex7ErmggESltduyYSpDocIBpPsrAu4ovLqsVe2+PEcfxIGgAnEWF2ewYZ/MHtVfUT4GKshmgXpzxOUNV3F0k23yg4aWA4xoK6CMESTwq5qJ+mcIZmMWnAMJOmmnBW909iPMm2Q4Mcvayw2IZfODdUflnyLXXhiXJHUPgNSVON9X5j4VQubtsT9AP9wGUi0gJsAKrBCo6q6v3FE7H0dA+V/6JPSFE40UmXmgngl476qtCUyg/T2dszQnxcOWJxaSpET6e5NsLB/qhvXW2d6gV9Pls4RuFULp5uf0XkA8DHgBXAM8ApWJ07z/VftPKhEhZ9wtQstfGkMhwbn5UqA4b0NFSHGBrLPDnu6LRaSpeLwnEUwsoWf+JJTsuDQitFOzSbGE7F49VR+zHgNcAeVT0HeBXQ6btUZUZAhKVNNaUWIyfhYICqUIChscSElWNcaqWjoSpENJ4kPp5+bfQOOwV5fZkonEjImg5W+pTA0FIbYTyp7OsZASYVRr54LYhqKD+8zkZRVY2KCCJSpaovicimokiWAxE5G/gc8AJwi6reW6xj3fqhU4s1tO9Yd9UJBsdmr3CnIT31KVmD6eIXOw4Psbihqmy+o5PXLuTebZ2+KRznM+/qsiy55ppCK1AHqIsEjYVTwXhVOO0i0gz8EvidiPQCB7weVES+D1wEHFbVo1O2XwB8FQgC31XVL2UZRoEhrFhSu1cZ5ip1VbbCmcU6aob0pK6LSqtwOodYt6h8sq0+9Np1HLuiiZN9WoTquOh2dQ1TXxWasKAKwZS3qWy8Vhp4i/3vtSJyD9AE/CaP494AfAO40dlg12j7JtZan3bgcRG5HUv5TC+fcwXwgKreJyJtwJeBss+Wmw3qq0IMpygc41IrHc65z5Q4sLdnhPO3tM2mSFkREU4/otW38VItnELdaQ5NtRFj4VQwec9G9mQfBN4J3OTxvfeLyJppm08CtqvqTgARuQW4VFW/iGUNZaIXyJhWIyJXAVcBrFq1youYFUmdnRnlBKtNllrpcM69k6aeymhsnK6hGCt8CtCXIwvsJIGe4RjLm/1J/W6qCZl6ahWMKxtXRBpF5B9F5Bsi8nq7EduHgZ3An/sky3JgX8rzdntbJpneKiL/DfwIy1pKi6peD3wGeCoSKf9Ms0JpmOZSMxZO6ZhcFzVzgtzfZwXSV7SUfzJKvqSmV/tl4TTXRFxV4DaUJ25nox9hWRKPAB8APoHVB+dSVX3GJ1nSrarUTDur6s+Bn7sZWFV/BfzqxBNP/GCeslUM9dUhhjsTEwsOG0xadMmoz+JS29c7CsxthVNfFaLevgFa5FMdwsaaEAM56tMZyhe3s9E6VT0GQES+C3QBq1R10EdZ2oGVKc9XkEdCQjoKLW1TSUwmDRiXWqlpqMqscNonFM7cdamB5VYbGkuwwqfMt/qqcFoXpaEycJs2MuETUNVxYJfPygbgcWCDiKwVkQhWbOh2n48x53FcakPRBKGAUB02ZdxLRbYYTnvPCJFgwLc7/3LF6T7rlyXnpP0nkxmdH4Yyxu1sdJyIDNiPQeBY538RGfB6UBG5Gcs9t0lE2kXkSlVNAB8G7gK2Areq6gtex05HobXUKok6e7Fh70ic+upQQT3kDYVRHQ4QDMjEItxU2ntHWdZcXfKincXm/aetBmDzEn/aVkxURI8ZK6cScVtLLejnQVX1sgzb7wDu8PNYML9cak6g+lD/qEkYKDEiYtdTm5k0cGggWhHVKwrlqrPWc/Fxy3z7rKlrm4y7uPKYF/6W+WThOD/Ig/1RGqrMD7LU1FeFGEzjUjs8GKWtcW670xz8VKz109qoGyqLeaFwRORiEbm+v7+/1KIUHecHn62o8wAAGdVJREFUeWggaqoMlAEN1aEZLjVVpWNgjMWz2KtoruBYNdmqcBvKl3mhcOaTheNUjO4bmb1un4bMOEHuVAZGE8QSSRb71AZgPjHpUjOLPysRr+0JqoC3AWtS36uqn/VXLEO+pLYjMD7u0lNfFaJ7eOpCxY7BKABtxsLxTINxqVU0Xi2c24BLgQQwnPIoa+aVSy1F4ZheOKWnvjo8w/1zeGAMwFg4eTChcIxLrSLxOiOtUNULiiJJEZlvlQYcTJZa6Wmontn1s2PAWDj5Up9lMa2h/PFq4TwsIscURRKDL9RHjEutnLAW4k6NNxwetC2ceZKl5id1kRAipM38M5Q/Xm+BzwAuF5FdwBhW/TNV1WN9l8xH5tM6nLqqySVTJkut9NSndP0MB637u46BKA1VIWoj5vvxSiAg1EdmZv4ZKgOvV/yFRZGiyMwnl1ooOGm0miy10lOfEnNoscv1dw6OschYN3lTX51+Ma2h/PHkUlPVPUAzcLH9aLa3GcoQE8MpPZMtCibvyDsGorQ1mPhNvjgVqA2VhyeFIyIfw2q2tth+/FhEPlIMwQyFY2I4pcdROKm1vw4Pjs2bKgPFIN3aJkNl4PUW+ErgZFUdBhCR67CKcH7db8EMhWPSoktPraNw7AnSqjIQNVUGCqC+Omy6flYoXmckAcZTno+TvnFaWTGfkgYA/ud9J/L03l42tvlTodeQP/V2EsfwmPWzGY6NM5ZI0lo/97vPFouGqhD7e0dKLYYhD7wqnB8Aj4rIL+znbwa+569I/jOfkgYAzt/Sxvlb2kothgEmMtFGbJdar111oLnWKJx8MS61ysWTwlHVL4vIfcDpWJbNX6jq00WRzGCYA9RFHJeaZeH0jlgKp8UonLyprzJp0ZWKZye/qj4JPFkEWQyGOUet7VKbsHBGrNjDgjqT0JEv9dUhhmPjjCeV4BxvYDfXcJWlJiIP2n8HUzp/5t3x02CYLzgWzpBj4RiXWsFka91tKG/cdvw8w/5bNlFoEQkAnwMagSdU9YclFslgmEF1OEBAUi0cS+EsMAonbxpS1jY11RhLsZLwug7nOjfbXIzzfRE5LCLPT9t+gYhsE5HtIvLJHMNcCiwH4kC7VxkMhtlARKiLhCZjOMMxRKDRTJR541RvMNUGKg+vxTvPT7Mtn3I3NwBTqk6LSBD4pj3eFuAyEdkiIseIyK+nPRYDm4BHVPXvgKvzkMFgmBVqq4JTYjjNNWETeyiA2ogTFxvPsaeh3HDlUhORq4G/AtaJyLMpLzUAD3s9qKreLyJrpm0+CdiuqjvtY94CXKqqXwQuSiNTO+B0tsp45YnIVcBVAKtWrfIqqsFQMHURK8gNlkvNZKgVhtPVdmTMKJxKw22W2k+AO4EvAqmurkFV7fFJluXAvpTn7cDJWfb/OfB1ETkTuD/TTqp6PXA9wIknnqg+yGkweKK2KjhRaaB3JDZRxNOQH46FY5IGKg+3SQP9QD+Wm6sF2ABUg+WjVtWME74H0vkYMioIVR3BKrWTe+B5VmnAUF5YMRxn4WecZc2mrE0hOCWbRmJG4VQantbhiMgHgI8BK4BngFOwaqmd64Ms7cDKlOcrgAM+jGswlJS6qhCddtO13pEYRy1rLLFElY1TvWHYxHAqDq9JAx8DXgPsUdVzgFcBnT7J8jiwQUTWikgEeCdwux8Dq+qvVPWqpqYmP4YzGDxRGwlOVIs2LrXCcZoMjhiXWsXhVeFEVTUKICJVqvoSVraYJ0TkZizLaJOItIvIlaqaAD4M3AVsBW5V1Re8jp3heBeLyPX9/f1+DGcweKIuEmJkbJzR2DjReNIkDRRIdSiIiLFwKhGvpW3aRaQZ+CXwOxHpJQ+3l6pelmH7HcAdXsdzcbx5VbzTUF44SQOTddTMGpxCCASE2nDQWDgViNfinW+x/71WRO4BmoDf+C6Vz5ikAUMpqa8KMRxL0GOXtTEutcKprQoZC6cC8epSm0BV71PV21U1lnvv0mJiOIZSUhsJkVQ41B8FTKVoP6iLBE2WWgXiNUvtROAaYHXqe1X1WJ/l8hVj4RhKiRPkbrebhplK0YVTm1IuyFA5eI3h3AR8AngOSPovTnEwMRxDKXHSeNt7RwFTKdoP6qsm1zYZKgevCqdTVX1JVZ5NjIVjKCV1EcfCsRWOKdxZMLVVwYlWD4bKwWsM519E5LsicpmIvNV5FEUyHzExHEMpcWp/7e8bpbE6RCiYd+jUYJNan85QOXi1cP4C2AyEmXSpKVZdM4PBkIbUGI7JUPOH2ohJi65EvCqc41T1mKJIYjDMUZwYTu9InNUL60oszdygzqRFVyRebfs/isiWokhSREylAUMpcdpMg1n06Re1Ji26IvGqcM4AnrG7cj4rIs9N649TlpgYjqGU1NouNTCLPv2iripEfFwZSxgrp5Jw7VITEQH+EthTPHEMhrlHQ3WqhWMUjh/UhC0lHo0lqQoFc+xtKBdcKxxVVRH5iqqeUEyBDIa5RlUoSCggJJLKAmPh+EKNnWoeTYzThHFTVgr5xHBeUxRJDIY5jNjtBZtNDMcXqsPW1DVqEgcqCq9ZaucAHxKR3cAwVpdONaVtDIbsJJJW89qFxsLxBcelNho3CqeS8KpwLiyKFEXGlLYxlBq1m6VvbGsorSBzhGqjcCoSTy41Vd0DNAMX249me5vBYMhCxK4usMasw/GFyaQBo3AqCa/Voj8GfJDJygI/FpHrVfXrvktmMMwh7vybMznYFyUQkFKLMicwFk5l4tWldiVwsqoOA4jIdVitoo3CMRiysH5RPesX1ZdajDmDk6VmFE5l4VXhCJD6DY/b22YdETkTeDfWZ9iiqqeVQg6DwTD7TLjU4hXTJcWAd4XzA+BREfmF/fzNwPe8HlREvg9cBBxW1aNTtl8AfBUIAt9V1S9lGkNVHwAeEJE3A497lcFgMFQuxqVWmXhSOKr6ZRG5Dzgdy7L5C1V9Oo/j3gB8A7jR2SAiQeCbwPlAO/C4iNyOpXy+OO39V6jqYfv/dwEfyEMGg8FQoUws/DRJAxWFK4UjIqtUdS+Aqj4JPFnIQVX1fhFZM23zScB2Vd1pH/MW4FJV/SKWNZRWLqBfVQeyyH4VcBXAqlWrChHbYDCUCU7WX9RYOBWF27ToXzr/iMj/K5Isy4F9Kc/b7W3ZuBLLzZcRVb0e+AzwVCRiFt0ZDHOBkJ3t5yyoNVQGbhVOamLAumIIQvrkg6xXk6r+i6o+nGtgUy3aYJhbBAJCQGDcKJyKwq3C0Qz/+0k7sDLl+QrggB8Dm344BsPcIxQMGAunwnCbNHCciAxgWSE19v8wWUut0QdZHgc2iMhaYD/wTqyEAIPBYJhBKCAkxk1adCXhysJR1aCqNqpqg6qG7P+d556VjYjcjLVgdJOItIvIlaqaAD4M3AVsBW5V1Re8jp1BfuNSMxjmGE7LB0Pl4HUdji+o6mUZtt8B3OH38Uy1aINh7mG51IyFU0l47YdTkRgLx2CYe4QCYpIGKgzXCkcsVubes/wwSQMGw9wjFBDi40bhVBKuFY6qKinrcSoJY+EYDHOPYNBYOJWGaTFtMBgqknDApEVXGvm0mP5LEdmDaTFtMBhKSNCkRVccpsW0wWCoSMzCz8rDa7XoPSLSAmwAqlNeMm2mDQbDrGIWflYeXltMfwD4GFbZmWeAU7AWcJ7rv2gGg8GQmaBZ+FlxeE0a+BjwGmCPqp4DvAro9F0qnzFp0QbD3CNsstQqDq8KJ6qqUQARqVLVl4BN/ovlLyYt2mCYe1hJA0bhVBJekwbaRaQZaz3O70SkF58qOhsMBoMXwsEAw4lEqcUweMBr0sBb7H+vFZF7gCbgTt+lMhgMhhyYGE7l4cmlJiLXOf+r6n2qejvwed+l8hkTwzEY5h6hQMC41CoMrzGc89NsK/u1OSaGYzDMPUzxzsrDlUtNRK4G/gpYJyLPOpuBeuChIslmMBgMGQkG5f9v7+6D5arrO46/P/fepBECZCBoIyENJQ+UsZRAhNJMI4oyYHmwU2KTCi0lA8UWMAVraX2gDzphRsbpVB40BoxUCQIpFQEDTkkEmUSBhIQEDEYI5ZbUAMpDFArBb//4/TbZ3tzde2+ye885ez+vmZ27e/ac3/ne3Xv3u7+H8/vxppcnqJTB9uHcROqrWQhcXrf91Yj4WcujMjMbwCjXcCpnsE1q00hDoudFxDPAe4B/JQ0eOLBt0ZmZNdDtPpzKGWzC+TLwBoCk2cCVwI3Ay8Ci9oTWnKRJku6QdIOkywc+wsw6yahuecXPihlswumuazr7Y2BRRCyLiE8DQ56COSeJbZI29Nl+iqRNkjYPIolMA+6KiPOAI4cag5lVmy/8rJ5BJxxJtf6ek4D76p4b6sWjAEuAU+o3SOoGriGNejsSmCfpSEm/LenOPre3A2uBuZLuA1bsQQxmVmE9vg6ncgabLJYC35P0AvAa8ACApCmkZrUhiYj7JU3us/k4YHNEPJXLvhk4MyIWAqf1LUPSx4Erclm3AV8dahxmVl093V0eNFAxg0o4EfE5Sf8JTADuzctNQ6ohXdyiWA4Bnq173Asc32T/5aRBC38CbGm0k6QLgAsAJk2atPdRmlkp9HSJN708QaUMujksIlb3s+3JFsai/k7bJJ4NwFkDFRoRiyRtBU4fPXr0sXsRn5mVSI9ni66coc400E69wKF1jyfiiUHNrIHurrTi564GFyu7MiWch4Cpkg6TNBqYC9zRioI9tY1Z5xnVlRpFXMupjkISjqSlpJVCp0vqlTQ/InYAFwH3AE8At0TExhadz5N3mnWY7u6UcDxSrTr2ZEjzXouIeQ223w3c3YbzfRv49syZM89vddlmVoyeLiecqilTk1rbuIZj1nl6utLH11u++LMyRkTCcR+OWefpyU1qnjG6OkZEwnENx6zz7KzhuEmtMkZEwnENx6zz1PpwfPFndYyIhGNmnafbw6IrZ0QkHDepmXWeHg+LrpwRkXDcpGbWeWp9OF6ioDpGRMIxs86zq4bjPpyqGBEJx01qZp1n54WfruFUxohIOG5SM+s8Pd25Sc19OJUxIhKOmXWeHo9SqxwnHDOrpO6dTWruw6kKJxwzq6RRHhZdOSMi4XjQgFnn6a4Ni/YotcoYEQnHgwbMOo9HqVXPiEg4ZtZ5atfheNBAdTjhmFkl7Zy80wmnMpxwzKySupQSToQTTlVUNuFIOlLSLZKuk3RW0fGYmVlzhSQcSTdI2iZpQ5/tp0jaJGmzpMsHKOZU4IsR8VHgT9sWrJmZtURPQeddAlwN3FjbIKkbuAb4ANALPCTpDqAbWNjn+POAfwOukHQGcNAwxGxmZnuhkIQTEfdLmtxn83HA5oh4CkDSzcCZEbEQOK1BUX+VE9W/NzqXpAuACwAmTZq0l5GbmdmeKlMfziHAs3WPe/O2fkmaLGkRqZb0+Ub7RcSiiJgZETMPPvjglgVrZmZDU1STWn/Uz7aGw08iYgu55jJgwdLpwOlTpkzZs8jMzGyvlamG0wscWvd4IvBcQbGYmVmLlSnhPARMlXSYpNHAXOCOVhTsqW3MzIpX1LDopcAqYLqkXknzI2IHcBFwD/AEcEtEbGzR+Tx5p5lZwTSSrtKV9DzwzB4ePh54oYXhtFuV4q1SrFCteKsUKzjedtqbWH8jIvZ61NWISjh7Q9LDETGz6DgGq0rxVilWqFa8VYoVHG87lSHWMvXhmJlZB3PCMTOzYeGEM3iLig5giKoUb5VihWrFW6VYwfG2U+Gxug/HzMyGhWs4ZmY2LJxwzMxsWDjhDGCIa/QUqtE6Q2Ul6VBJKyQ9IWmjpI8VHVMjksZI+qGkdTnWfyw6psGQ1C1praQ7i45lIJK2SHpM0qOSHi46nmYkjZN0m6Qf5b/fE4qOqRFJ0/NrWru9ImlBIbG4D6exvPTBk9St0QPMi4jHCw2sAUmzge3AjRHxrqLjGYikCcCEiFgjaT/gEeBDZXx9JQnYNyK2SxoFfB/4WESsLji0piRdCswE9o+IRst8lIKkLcDMiCj9hZSSvgY8EBGL81Rc+0TES0XHNZD8mfbfwPERsacXwe8x13Ca27lGT0S8AdwMnFlwTA1FxP3Az4qOY7AiYmtErMn3XyVNadRwSYoiRbI9PxyVb6X+tiZpIvAHwOKiY+kkkvYHZgPXA0TEG1VINtlJwE+KSDbghDOQIa3RY3suL8g3A/hBsZE0lpunHgW2Ad+NiNLGmv0L8AngV0UHMkgB3CvpkbxwYln9JvA88NXcXLlY0r5FBzVIc4GlRZ3cCae5Ia3RY3tG0lhgGbAgIl4pOp5GIuKtiDiatHTGcZJK22wp6TRgW0Q8UnQsQzArIo4BTiWt5ju76IAa6AGOAa6LiBnAL4BS9+8C5Ka/M4Bbi4rBCac5r9HTZrk/ZBnwjYhouFR4meTmk5XAKQWH0sws4IzcL3Iz8D5JXy82pOYi4rn8cxtwO6lJu4x6gd66Gu5tpARUdqcCayLip0UF4ITTXNvW6LGdHfHXA09ExBeKjqcZSQdLGpfvvw14P/CjYqNqLCL+LiImRsRk0t/tfRFxdsFhNSRp3zxwhNw8dTJQytGWEfE/wLOSpudNJwGlG+jSj3kU2JwG5VpiunQiYoek2ho93cANrVqjpx3yOkMnAuMl9QJXRMT1xUbV1CzgHOCx3DcC8PcRcXeBMTUyAfhaHuXTRVqvqfRDjSvkHcDt6TsIPcBNEbG82JCauhj4Rv4i+hTw5wXH05SkfUijbf+i0Dg8LNrMzIaDm9TMzGxYOOGYmdmwcMIxM7Nh4YRjZmbDwgnHzKxDDWVCX0mT8mS6ayWtl/TBVsfjhGOlJOmTeVbm9XmG2+Pz9gV5iGeRsY2T9Jd1j98p6bYBjlkpaWb7o2tOyX15PrDhPvdFkko9fLgDLWHwFyh/ijTcfwbp2q1rWx2ME46VTp7q/TTgmIg4inSRZW1OuwXAkBJOvnamVbF1A+OAnQknIp6LiLNadY42+yCwrp1TCDV5vW8ALmnXeW13/U3oK+lwScvznHUPSDqitjtQ+yJyAG2YVcUJx8poAvBCRPwvQES8EBHPSboEeCewQtIKAEnXSXq47xo1eW2Vz0j6PjCnvnBJSyR9Kf+zPZnnHUPS5LxtTb79Xt5+Ym5quAl4DLgSODzXvD6fj9uQ9+2WdJXSui7rJV3c95eTdLKkVfkct+a55Pruc76kh5TW31lWq9XlD4vV+bl/krS97pi/ydvXq/F6PR8BvpX3/2fVrUEk6XP5NW5YlqT/yB9UG1U3waak7TmeHwAnSLpS0uP5+Kvy+/hLYIuksk5ZM1IsAi6OiGOBj7OrJvMPwNn5ovG7SRe3tlZE+OZbqW7AWOBR0lpE1wLvqXtuCzC+7vGB+Wc3aX6zo+r2+0SD8pcAy0lfuKaS5sYaQ6o5jcn7TAUezvdPJE3QeFh+PBnYUFfezsfAR0lzw/X0iW8laV2a8cD9pLV1AP4W+Ew/MR5Ud/+zpA8IgDtJazIBXAhsz/dPJn2QKP9edwKz+yn3GWC/urjX5PtdwE+Ag5qVVff7vI009cxB+XEAH67tA2xi14Xl4+rO/0ngsqL/xkbSrc/f51jgtfz/Vbs9kZ+7tPbeACeQpuvpamUsntrGSifSImfHAr8PvBf4pqTLI2JJP7t/OH/T7iHVjI4E1ufnvtnkNLdExK+AH0t6CjgCeBq4WtLRwFvAtLr9fxgRTw8i/PcDX4qIHfl36bs+0e/mGB/M07iMBlb1U867JH2W1Hw3ljS9EqQPgg/l+zcBV+X7J+fb2vx4LClp3t+n3AMjrT1ERGyR9KKkGaSpZdZGxIuSmpV1iaQ/zNsPzdtfJL1ey/L2V4DXgcWS7iIlrJptpNfaitEFvBRp1vO+5pP7eyJilaQxpC9I21p1ciccK6WIeItUK1gp6THgz0g1k50kHUZqEnh3RPxc0hJSTaXmF81O0c/jvwZ+CvwO6R/z9UGW9f/C6qfsvs9/NyLmDVDOEtLqp+sknUuqZQ103oUR8eUB9tshqSsnW0iLs50L/Dqpj6VhWZJOJCXUEyLil5JWsuv1fj2/Z0Sag/A40qSWc4GLgPfl/caQvmFbASLiFUlPS5oTEbcqfes5KiLWAf9Fes+WSPot0nv1fCvP7z4cKx2lNdin1m06mtQUBPAqsF++vz8pEbws6R2k6dcHa46kLkmHkxbU2kTqKN2aP4zPITXT9ac+hr7uBS6U1JN/lwP7PL8amCVpSn5+H0nT2N1+wFal5Rs+0uf4P8r359Ztvwc4r9YfJOkQSW/vp9xNpN+35nbSt9p3s6sW1aisA4Cf52RzBKm2tpt83AGRJmFdQHr/aqZR0lmgO5HShL6rgOmSeiXNJ/09zZe0DtjIrlWMLwPOz9uXAudGbl9rFddwrIzGAl9UWg5gB7AZqHVQLwK+I2lrRLxX0lrSP81TwINDOMcm4HukpqQLI+J1SdcCyyTNAVbQoFaTm50ezAMFvgNcU/f0YtKH6npJbwJfAa6uO/b5XGNZKunX8uZPkfqr6n2atPrpM6SBCrUEtwD4uqTLgLuAl3O59+ZvpatyU9124Gx2bw65i1Rb2pyPe0NpAMZLdTWURmUtJyXT9fn1W93f65Nj/VZukhGp5lgzC2g0oMFarElNereh0hHxOOn9aRvPFm0jTm56uzMiml47U0Z5tNprERGS5pIGEJw50HF1x08AboyID+THXcAaYE5E/LgtQe869wzg0og4p53nsfJyDcesWo4lDWwQ8BJw3lAOjoitkr6idOHnRFKH/u3tTjbZeFLNzUYo13DMzGxYeNCAmZkNCyccMzMbFk44ZmY2LJxwzMxsWDjhmJnZsPg/1knBdOOtZfMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "maxima = extrema(12,5)\n",
    "maxima, times, widths = sfr_end(12,5)\n",
    "#maxima, times, widths = filter_peaks(maxima, times, widths)\n",
    "#print(maxima)\n",
    "#print(np.diff(maxima))\n",
    "plot_tend(12,5)\n",
    "smoothed_sfr, n, bins = smoothing(12,5)\n",
    "\n",
    "t, maxes = t_end(12)\n",
    "#[plt.axvline(t[i], color = 'r') for i in range(len(t))]\n",
    "\n",
    "#for i in range(len(times)):\n",
    "    #print(smoothed_sfr[maxima[i]:maxima[i]+ widths[i]][times[i]])\n",
    "#print(smoothed_sfr[maxima])\n",
    "print(widths)\n",
    "print(times)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 563,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n",
      "29\n",
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n",
      "35\n",
      "36\n",
      "37\n",
      "38\n",
      "39\n",
      "40\n",
      "41\n",
      "42\n",
      "43\n",
      "44\n",
      "45\n",
      "46\n",
      "47\n",
      "48\n",
      "49\n",
      "50\n",
      "51\n",
      "52\n",
      "53\n",
      "54\n",
      "55\n",
      "56\n",
      "57\n",
      "58\n",
      "59\n",
      "60\n",
      "61\n",
      "62\n",
      "63\n",
      "64\n",
      "65\n",
      "66\n",
      "67\n",
      "68\n",
      "69\n",
      "70\n",
      "71\n",
      "72\n",
      "73\n",
      "74\n",
      "75\n",
      "76\n",
      "77\n",
      "78\n",
      "79\n",
      "80\n",
      "81\n",
      "82\n",
      "83\n",
      "84\n",
      "85\n",
      "86\n",
      "87\n",
      "88\n",
      "89\n",
      "90\n",
      "91\n",
      "92\n",
      "93\n",
      "94\n",
      "95\n",
      "96\n",
      "97\n",
      "98\n",
      "99\n",
      "100\n",
      "101\n",
      "102\n",
      "103\n",
      "104\n",
      "105\n",
      "106\n",
      "107\n",
      "108\n",
      "109\n",
      "110\n",
      "111\n",
      "112\n",
      "113\n",
      "114\n",
      "115\n",
      "116\n",
      "117\n",
      "118\n",
      "119\n",
      "120\n",
      "121\n",
      "122\n",
      "123\n",
      "124\n",
      "125\n",
      "126\n",
      "127\n",
      "128\n",
      "129\n",
      "130\n",
      "131\n",
      "132\n",
      "133\n",
      "134\n",
      "135\n",
      "136\n",
      "137\n",
      "138\n",
      "139\n",
      "140\n",
      "141\n",
      "142\n",
      "143\n",
      "144\n",
      "145\n",
      "146\n",
      "147\n",
      "148\n",
      "149\n",
      "150\n",
      "151\n",
      "152\n",
      "153\n",
      "154\n",
      "155\n",
      "156\n",
      "157\n",
      "158\n",
      "159\n",
      "160\n",
      "161\n",
      "162\n",
      "163\n",
      "164\n",
      "165\n",
      "166\n",
      "167\n",
      "168\n",
      "169\n",
      "170\n",
      "171\n",
      "172\n",
      "173\n",
      "174\n",
      "175\n",
      "176\n",
      "177\n",
      "178\n",
      "179\n",
      "180\n",
      "181\n",
      "182\n",
      "183\n",
      "184\n",
      "185\n",
      "186\n",
      "187\n",
      "188\n",
      "189\n",
      "190\n",
      "191\n",
      "192\n",
      "193\n",
      "194\n",
      "195\n",
      "196\n",
      "197\n",
      "198\n",
      "199\n",
      "200\n",
      "201\n",
      "202\n",
      "203\n",
      "204\n",
      "205\n",
      "206\n",
      "207\n",
      "208\n",
      "209\n",
      "210\n",
      "211\n",
      "212\n",
      "213\n",
      "214\n",
      "215\n",
      "216\n",
      "217\n",
      "218\n",
      "219\n",
      "220\n",
      "221\n",
      "222\n",
      "223\n",
      "224\n",
      "225\n",
      "226\n",
      "227\n",
      "228\n",
      "229\n",
      "230\n",
      "231\n",
      "232\n",
      "233\n",
      "234\n",
      "235\n",
      "236\n",
      "237\n",
      "238\n",
      "239\n",
      "240\n",
      "241\n",
      "242\n",
      "243\n",
      "244\n",
      "245\n",
      "246\n",
      "247\n",
      "248\n",
      "249\n",
      "250\n",
      "251\n",
      "252\n",
      "253\n",
      "254\n",
      "255\n",
      "256\n",
      "257\n",
      "258\n",
      "259\n",
      "260\n",
      "261\n",
      "262\n",
      "263\n",
      "264\n",
      "265\n",
      "266\n",
      "267\n",
      "268\n",
      "269\n",
      "270\n",
      "271\n",
      "272\n",
      "273\n",
      "274\n",
      "275\n",
      "276\n",
      "277\n",
      "278\n",
      "279\n",
      "280\n",
      "281\n",
      "282\n",
      "283\n",
      "284\n",
      "285\n",
      "286\n",
      "287\n",
      "288\n",
      "289\n",
      "290\n",
      "291\n",
      "292\n",
      "293\n",
      "294\n",
      "295\n",
      "296\n",
      "297\n",
      "298\n",
      "299\n",
      "300\n",
      "301\n",
      "302\n",
      "303\n",
      "304\n",
      "305\n",
      "306\n",
      "307\n",
      "308\n",
      "309\n",
      "310\n",
      "311\n",
      "312\n",
      "313\n",
      "314\n",
      "315\n",
      "316\n",
      "317\n",
      "318\n",
      "319\n",
      "320\n",
      "321\n",
      "322\n",
      "323\n",
      "324\n",
      "325\n",
      "326\n",
      "327\n",
      "328\n",
      "329\n",
      "330\n",
      "331\n",
      "332\n",
      "333\n",
      "334\n",
      "335\n",
      "336\n",
      "337\n",
      "338\n",
      "339\n",
      "340\n",
      "341\n",
      "342\n",
      "343\n",
      "344\n",
      "345\n",
      "346\n",
      "347\n",
      "348\n",
      "349\n",
      "350\n",
      "351\n",
      "352\n",
      "353\n",
      "354\n",
      "355\n",
      "356\n",
      "357\n",
      "358\n",
      "359\n",
      "360\n",
      "361\n",
      "362\n",
      "363\n",
      "364\n",
      "365\n",
      "366\n",
      "367\n",
      "368\n",
      "369\n",
      "370\n",
      "371\n",
      "372\n",
      "373\n",
      "374\n",
      "375\n",
      "376\n",
      "377\n",
      "378\n",
      "379\n",
      "380\n",
      "381\n",
      "382\n",
      "383\n",
      "384\n",
      "385\n",
      "386\n",
      "387\n",
      "388\n",
      "389\n",
      "390\n",
      "391\n",
      "392\n",
      "393\n",
      "394\n",
      "395\n",
      "396\n",
      "397\n",
      "398\n",
      "399\n",
      "400\n",
      "401\n",
      "402\n",
      "403\n",
      "404\n",
      "405\n",
      "406\n",
      "407\n",
      "408\n",
      "409\n",
      "410\n",
      "411\n",
      "412\n",
      "413\n",
      "414\n",
      "415\n",
      "416\n",
      "417\n",
      "418\n",
      "419\n",
      "420\n",
      "421\n",
      "422\n",
      "423\n",
      "424\n",
      "425\n",
      "426\n",
      "427\n",
      "428\n",
      "429\n",
      "430\n",
      "431\n",
      "432\n",
      "433\n",
      "434\n",
      "435\n",
      "436\n",
      "437\n",
      "438\n",
      "439\n",
      "440\n",
      "441\n",
      "442\n",
      "443\n",
      "444\n",
      "445\n",
      "446\n",
      "447\n",
      "448\n",
      "449\n",
      "450\n",
      "451\n",
      "452\n",
      "453\n",
      "454\n",
      "455\n",
      "456\n",
      "457\n",
      "458\n",
      "459\n",
      "460\n",
      "461\n",
      "462\n",
      "463\n",
      "464\n",
      "465\n",
      "466\n",
      "467\n",
      "468\n",
      "469\n",
      "470\n",
      "471\n",
      "472\n",
      "473\n",
      "474\n",
      "475\n",
      "476\n",
      "477\n",
      "478\n",
      "479\n",
      "480\n",
      "481\n",
      "482\n",
      "483\n",
      "484\n",
      "485\n",
      "486\n",
      "487\n",
      "488\n",
      "489\n",
      "490\n",
      "491\n",
      "492\n",
      "493\n",
      "494\n",
      "495\n",
      "496\n",
      "497\n",
      "498\n",
      "499\n",
      "500\n",
      "501\n",
      "502\n",
      "503\n",
      "504\n",
      "505\n",
      "506\n",
      "507\n",
      "508\n",
      "509\n",
      "510\n",
      "511\n",
      "512\n",
      "513\n",
      "514\n",
      "515\n",
      "516\n",
      "517\n",
      "518\n",
      "519\n",
      "520\n",
      "521\n",
      "522\n",
      "523\n",
      "524\n",
      "525\n",
      "526\n",
      "527\n",
      "528\n",
      "529\n",
      "530\n",
      "531\n",
      "532\n",
      "533\n",
      "534\n",
      "535\n",
      "536\n",
      "537\n",
      "538\n",
      "539\n",
      "540\n",
      "541\n",
      "542\n",
      "543\n",
      "544\n",
      "545\n",
      "546\n",
      "547\n",
      "548\n",
      "549\n",
      "550\n",
      "551\n",
      "552\n",
      "553\n",
      "554\n",
      "555\n",
      "556\n",
      "557\n",
      "558\n",
      "559\n",
      "560\n",
      "561\n",
      "562\n",
      "563\n",
      "564\n",
      "565\n",
      "566\n",
      "567\n",
      "568\n",
      "569\n",
      "570\n",
      "571\n",
      "572\n",
      "573\n",
      "574\n",
      "575\n",
      "576\n",
      "577\n",
      "578\n",
      "579\n",
      "580\n",
      "581\n",
      "582\n",
      "583\n",
      "584\n",
      "585\n"
     ]
    }
   ],
   "source": [
    "# Saves Plots of End Times on log scale\n",
    "\n",
    "for i in range(len(mass)):\n",
    "    print(i)\n",
    "    fig = plt.figure()\n",
    "    plot_tend(i,5)\n",
    "    plt.title(\"Smoothed star formation rate of star particles in galaxy %i\" % i)\n",
    "    plt.ylabel(\"Star Formation Rate ($M_\\odot yr^{-1}$)\")\n",
    "    plt.xlabel(\"Star particle age (years)\")    \n",
    "    fig.savefig('End Times Logscale/galaxy_%i.png' %i)\n",
    "    plt.tight_layout()\n",
    "    fig.clf()\n",
    "    plt.close()"
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}