{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to programming for Geoscientists through Python\n",
    "\n",
    "# Lecture 6 solutions\n",
    "\n",
    "## [Gerard J. Gorman (g.gorman@imperial.ac.uk)](http://www.imperial.ac.uk/people/g.gorman), [Nicolas Barral (n.barral@imperial.ac.uk)](http://www.imperial.ac.uk/people/n.barral)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Read a two-column data file**</br>\n",
    "The file *data/xy.dat* contains two columns of numbers, corresponding to *x* and *y* coordinates on a curve. The start of the file looks like this:\n",
    "\n",
    "-1.0000   -0.0000</br>\n",
    "-0.9933   -0.0087</br>\n",
    "-0.9867   -0.0179</br>\n",
    "-0.9800   -0.0274</br>\n",
    "-0.9733   -0.0374</br>\n",
    "\n",
    "Make a program that reads the first column into a list *x* and the second column into a list *y*. Then convert the lists to arrays, and plot the curve. Print out the maximum and minimum y coordinates. (Hint: Read the file line by line, split each line into words, convert to float, and append to *x* and *y*.)</br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## STAGE 1: input\n",
    "\n",
    "# Open data file\n",
    "infile = open(\"data/xy.dat\", \"r\") # \"r\" is for read\n",
    "\n",
    "# Initialise empty lists\n",
    "xlist = []\n",
    "ylist = []\n",
    "\n",
    "# Loop through infile and write to x and y lists\n",
    "for line in infile:\n",
    "    line = line.split() # convert to list by dropping spaces\n",
    "    xlist.append(float(line[0])) # take 0th element and covert to float\n",
    "    ylist.append(float(line[1])) # take 1st element and covert to float\n",
    "\n",
    "\n",
    "# Close the filehandle\n",
    "infile.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First ten elements in x and y lists:\n",
      "\n",
      "X = -1.00, Y = -0.00\n",
      "X = -0.99, Y = -0.01\n",
      "X = -0.99, Y = -0.02\n",
      "X = -0.98, Y = -0.03\n",
      "X = -0.97, Y = -0.04\n",
      "X = -0.97, Y = -0.05\n",
      "X = -0.96, Y = -0.06\n",
      "X = -0.95, Y = -0.07\n",
      "X = -0.95, Y = -0.08\n",
      "X = -0.94, Y = -0.09\n",
      "\n",
      "The min and max y are: -0.95 and 0.95\n"
     ]
    }
   ],
   "source": [
    "## STAGE 2: view file\n",
    "\n",
    "# Sanity check -- Let's see what our lists look like for the first 10\n",
    "print(\"First ten elements in x and y lists:\\n\") # N.B. \"\\n\" adds a newline\n",
    "for i in range(10):\n",
    "    print(\"X = %.2f, Y = %.2f\" % (xlist[i], ylist[i]))\n",
    "\n",
    "# Max and min y coords\n",
    "print(\"\\nThe min and max y are: %.2f and %.2f\" % (min(ylist), max(ylist)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl83HWd+PHXO3ea+2qbs016Jz3SNi0F5JCzgLYIKGVR\nAUHUXXRdfrri+vuJy8pP1N8uuivqdhGooCLishRaKJRDzpamB23TM2naNHeapDmbcz6/P+Y7dQhJ\nc83Md4738/GYR2a+5zvfTOY9388pxhiUUkqpyQqzOwCllFLBQROKUkopj9CEopRSyiM0oSillPII\nTShKKaU8QhOKUkopj9CEopRSyiM0oSillPIITShKKaU8IsLuAHwpPT3dzJw50+4wlFIqoOzcufOU\nMSZjtO1sTSgi8hjwKaDRGLNwmPUC/By4FugGbjfG7LLW3Qb8b2vTHxpjNox2vpkzZ1JaWuqp8JVS\nKiSIyImxbGd3kdcTwOpzrL8GmGM97gZ+BSAiqcD9wHnASuB+EUnxaqRKKaXOydaEYox5C2g5xyZr\ngd8ap21AsohkAlcDrxpjWowxrcCrnDsxKaWU8jK771BGkw2cdHtdbS0bablSSimb+HtCmTQRuVtE\nSkWktKmpye5wlFIqaPl7QqkBct1e51jLRlr+McaY9caYEmNMSUbGqI0UlFJKTZC/J5SNwBfFaRXQ\nZoypA7YAV4lIilUZf5W1TCmllE3sbjb8B+BSIF1EqnG23IoEMMb8GtiMs8lwOc5mw3dY61pE5F+A\nHdahHjDGnKtyXymllJfZmlCMMbeMst4AfzfCuseAx7wRl1LBYmDQwWuHGqlo6mRWRjyXzZ9KZLi/\nF0yoQBVSPeWVCiUnmrv4ypM7OVTfcXbZ/OkJ/Orzy8lPj7MxMhWs9KuKUkGoru0M69Zvo769h1/e\nuoy9P7iKX926jMaOXm5Zv42TLd12h6iCkCYUpYLMoMPwzaf30Hamn9/ftYprF2WSGBPJNYsy+f2X\nz6Orb4C/f3o3A4MOu0NVQUYTilJB5pnSk2yvbOEHa4oozEr8yLr50xN58DOL2FV1msffPW5PgCpo\naUJRKoh09Q7wr68coWRGCp9dnjPsNmuWZHHpvAz+4/WjnO7u83GEKphpQlEqiPx+exWnOnv57rUL\ncA7WPbzvXrOAjt4B/uvtYz6MTgU7TShKBYm+AQePvVvJqoJUls849+Db86YnsLpoOk++f4KOnn4f\nRaiCnSYUpYLElrJ66tp6+MrFs8a0/VcvmUV7zwDPlFZ7OTIVKjShKBUknt5RRU5KLJfMHduYdUty\nk1mal8zvtp/A2YdYqcnRhKJUEKhq7ubd8mZuLsklLGzkupOhbj1vBseauth2TEcuUpOnCUWpIPDC\n3loAbhihZddIPrU4k/joCP5n97CDdSs1LppQlAoCL+6tY1leMtnJsePaLyYynMsXTGXLgXr6taOj\nmiRNKEoFuMpTXRysa+e6xVkT2v+6RZmc7u7n/YpmD0emQo0mFKUC3NYDDQBcXTRtQvtfPDeD+OgI\nNu2t82RYKgRpQlEqwL15pJG50+LJSZkyof1jIsO5Qou9lAdoQlEqgHX2DvBBZQufnDd1Use5Vou9\nlAdoQlEqgL1Xfor+QcMl88bW92QkF8/NIDYynK0HGzwUmQpFtiYUEVktIodFpFxE7htm/cMissd6\nHBGR027rBt3WbfRt5Er5hzcONxEfHUHJjNRJHScmMpzzZ6Xx1pEmD0WmQpFtMzaKSDjwCHAlUA3s\nEJGNxpgDrm2MMf/gtv3XgaVuhzhjjCn2VbxK+RtjDH853MgnZqcTFTH574YXz0nn9UONVDV3k5c2\nsfoYFdrsvENZCZQbY44ZY/qAp4G159j+FuAPPolMqQBwpKGT2rYeLp1kcZfLxdaQLX85qncpamLs\nTCjZwEm319XWso8RkRlAPvC62+IYESkVkW0icr33wlTKP71bfgr4ayKYrPz0OHJSYrXYS02YbUVe\n47QOeNYYM+i2bIYxpkZECoDXRWSfMaZi6I4icjdwN0BeXp5volXKB3YcbyEnJZascfaOH4mIcPHc\nDDbuqaV/0EFkuLbZUeNj5zumBsh1e51jLRvOOoYUdxljaqyfx4A3+Wj9ivt2640xJcaYkowMz3yT\nU8puxhh2HG9h5czJVcYPdfGcDDp7B9h1otWjx1Whwc6EsgOYIyL5IhKFM2l8rLWWiMwHUoD33Zal\niEi09TwduBA4MHRfpYJV5akuTnX2sSLfswnlwtlphIcJ71jFaUqNh21FXsaYARG5B9gChAOPGWPK\nROQBoNQY40ou64CnzUcnbFgA/KeIOHAmxYfcW4cpFex2HHcON7/Cw3coCTGRLMxKZHulDmevxs/W\nOhRjzGZg85Bl3x/y+gfD7PcesMirwSnlxz6obCUtLopZGXEeP/Z5BWk88d5xevoHiYkM9/jxVfDS\nWjelAtCO4y2UzExBZOyTaY3Vypmp9A042HPy9OgbK+VGE4pSAaa+rYeqlm6PF3e5rJiZigh8oMVe\napw0oSgVYLxVf+KSNCWS+dMT2V6pA0Wq8dGEolSA+fDkaaIiwijMSvTaOc7LT2XniVb6BnQ4ezV2\nmlCUCjB7q9soykr0asfD8/JT6el3sK+mzWvnUMFHE4pSAWRg0PkhvyQn2avnWWn1b9F6FDUemlCU\nCiDlTZ2c6R9kSW6SV8+TFh9Nfnocu6q0x7waO00oymuMMVSe6qL0eAuN7T12hxMU9p50FkEt9vId\nCkBxbjJ7Tp7mo32KlRpZoAwOqQLMlrJ6fvLyISqaus4uu3ReBvd/uoj8dM93xgsVH1afJiEmgvw0\n71/D4txknttdQ21bD9keGoBSBTe9Q1Ee5XAYHnjhAF95cicRYWE8+JmFPHHHCr55xRx2nmjlun9/\nm3eO6jhRE/Vh9WkW5yQRFub5Do1DLc1z3gXtqdIOjmpsNKEojzHG8P2N+3ns3Upuv2AmL3z9E9x6\n3gwunTeVb14xl1f/4RLyUqdw54Yd7K3WD6nx6ukf5FBdh0+KuwDmT08kKiKMPSe1HkWNjSYU5TFP\nbjvBU9uq+MolBdz/6cKPTUs7PSmG3911Hunx0XzlyZ20nem3KdLAdLCunQGH8XoLL5eoiDAWZiXq\nECxqzDShKI84UNvOD188yGXzp/Kdq+ePOMZUWnw0v7x1GY0dvfzLizpA9Hi4+oQszvFuCy93xbkp\n7Ktpo39QOziq0WlCUZM2MOjgW3/6kKQpkfz0psWjlu8vyU3mKxcX8OzOanae0H4OY3Wgtp3UuCgy\nk2J8ds7ivGR6+h0cru/w2TlV4NKEoiZtw/snOFDXzgNrikiLjx7TPvdcNpupCdH8cNNBbZY6RmW1\n7RRmJnplhOGRLM11Fq/t1mIvNQaaUNSknO7u4+dbj3Dx3AxWL5w+5v2mREXwD1fOZXfVad7SVl+j\n6h90cLihw6vjdw0nJyWWtLgobemlxkQTipqUX/2lgo7eAb57zcj1JiO5cVkOmUkxPPJGuZeiCx7H\nmrroG3BQ5OOEIiIszUvWll5qTGxNKCKyWkQOi0i5iNw3zPrbRaRJRPZYj7vc1t0mIketx22+jVwB\ntHT18dv3TrB2SRYLMsf/QRcVEcadn8jng8oW9lXrIITnUlbrvD6FE7jOk1Wcm0xFU5e2ylOjsi2h\niEg48AhwDVAI3CIihcNs+kdjTLH1eNTaNxW4HzgPWAncLyIpPgpdWR5/t5Iz/YP83SdnT/gYn1uR\ny5SocDa8f9xjcQWjA7XtREeE2TLKQHGu819L+w6p0dh5h7ISKDfGHDPG9AFPA2vHuO/VwKvGmBZj\nTCvwKrDaS3GqYXT09PPEe8dZXTSdOdMSJnycxJhIbliWzcYPaznd3efBCIPLgbp25k9PIMKLQ9aP\nZHFuEiLaY16Nzs6Ekg2cdHtdbS0b6kYR2Ssiz4pI7jj3VV7y5LYTdPQMTOruxOXmkjz6Bhxs3lfv\ngciCjzHG2cIry3f9T9wlxkSSnx6nc6OoUfl7pfwLwExjzGKcdyEbxnsAEblbREpFpLSpqcnjAYai\n/kEHG947zkVz0lnkgU52C7MTmZURx//sqfFAdMGntq2HtjP9Pm/h5W5hVhJlte22nV8FBjsTSg2Q\n6/Y6x1p2ljGm2RjTa718FFg+1n3djrHeGFNijCnJyMjwSOChbuuBBhrae7n9gpkeOZ6IcH1xNh9U\ntlDd2u2RYwaTshr7KuRdFmYnUnP6DC1dWiypRmZnQtkBzBGRfBGJAtYBG903EJFMt5drgIPW8y3A\nVSKSYlXGX2UtUz7w2/dPkJ0cy6XzpnrsmGuLnSWWGz+s9dgxg8Uhq5f6/OkTr6uarCKruM3V2kyp\n4diWUIwxA8A9OBPBQeAZY0yZiDwgImuszb4hImUi8iHwDeB2a98W4F9wJqUdwAPWMuVl5Y0dvH+s\nmVtX5RHuwSHU89KmsHxGCv+zu0Z7zg9xpKGDvNQpxEXbN32Rq/+LFnupc7F1gi1jzGZg85Bl33d7\n/l3guyPs+xjwmFcDVB/z1LYqosLD+FxJ7ugbj9Pa4iy+/3wZ5Y2dk2o5FmyONHQwd1q8rTEkT4ki\nJyWW/Voxr87B3yvllR/pHRjkud01XL1wOuljHLNrPK4snAbA1oONHj92oOobcHCsqYu5fpBgi7IS\n9Q5FnZMmFDVmrx9spO1MPzctz/HK8TOTYinKSuS1gw1eOX4gqjzVxYDDMM/G+hOXhVlJVJ7qoqNH\ne8yr4WlCUWP2513VTE2I5hOz0712jisWTGNnVSvNnb2jbxwCDjc4K+T94Q5lYbazYv5gnQ5lr4an\nCUWNyanOXt483MRnlmZ7tDJ+qCsWTMMYeOOw9hkCONrQQXiYUJDh+yFXhirKdlbMaz2KGokmFDUm\nG/fUMuAw3Oil4i6XhdmJTEuM1mIvy+H6DmamTSE6ItzuUJiaEENGQjT7temwGoEmFDUmf95VzaLs\nJK8XvYgIly+YxltHmugdGPTquQLBkYYOv6g/cVmYlcgBrZhXI9CEokZ1qL6dstp2bljmm+HSLps3\nla6+QXaeCO05OM70DXKipdsv6k9cFmYncbSxk55+Tfbq4zShqFH9964aIsKENUuyfHK+8wpSCQ8T\n3itv9sn5/FV5YyfGwDw/SihFWUkMOszZ3vtKudOEos7J4TBs3FPLpfOmjnm++MlKiIlkcU4S71aE\n9tTAZ1t4+VGRl6vHvFbMq+FoQlHntKuqlfr2Hj69JHP0jT3owlnp7K1uC+k+D0cbOoiKCGNG6hS7\nQzkrJyWWpNhI7eCohqUJRZ3Ti3vriIoI4/IF03x63gtmpzHoMGw/FrpDtB1u6GBWRrwtk2qNRERY\nmJ2og0SqYfnPO1X5HYfD8NL+Oi6dm0G8jwcmXJaXQnREWEgXex2p72CezWN4DacoK4lDdR30Dzrs\nDkX5GU0oakQ7q1ppaO/lusW+Le4CiIkMZ8XM1JCtmG/v6ae2rcev6k9cirIS6Rt0UN7YaXcoys9o\nQlEj2mRTcZfLBbPTONzQQVNH6A3DctSqkPenFl4uriFYtGJeDaUJRQ3LVdz1yXm+L+5yWVWQBkDp\n8dCrRznS4Pz27099UFzy0+KIiwrXinn1MZpQ1LBcxV3XLvJ9cZfLwqwkoiPC2HE89Do4Hm3oJDYy\nnOzkWLtD+ZiwMGFBplbMq4/ThKKGtWlvHdE2FncBREWEUZybTOmJ0LtDKW/qpCAjjjAvDsQ5GYXW\nECwOh86uqf7K1oQiIqtF5LCIlIvIfcOsv1dEDojIXhF5TURmuK0bFJE91mPj0H3VxDkchs376rjU\nxuIulxUzUymrbaerd8DWOHytorGT2VP9r4WXS1FWIl3W0DBKudiWUEQkHHgEuAYoBG4RkcIhm+0G\nSowxi4FngZ+4rTtjjCm2HmtQHlN6opXGjl6uW+yboVbOpWRmCoMOw56Tp+0OxWe6+waoOX2GWRn+\nnFCcFfNa7DU6h8NQ1dzNrqpW9tcEd2ddO79+rgTKjTHHAETkaWAtcMC1gTHmDbfttwGf92mEIWrz\nPqu4a/5Uu0Nh2YwURGDH8RYu9OLEXv7kWFMXgF/focyZFk9EmFBW286n/OCLhz862tDBY+8e55Wy\nepq7+s4uDw8TFucksbpoOrecl0diTKSNUXqWnQklGzjp9roaOO8c298JvOT2OkZESoEB4CFjzP94\nPsTQ4yru+uS8qcTZXNwFkBgTyfzpiZSGUMV8RZOzhZc/36FER4QzZ1qCDmU/jN6BQR566RAb3jtO\nVEQYVxVO5/xZaUxPiqG330FZbRtvHz3Fj146xC9eL+fOi/L52qWz/GLOm8my/xNjDETk80AJcInb\n4hnGmBoRKQBeF5F9xpiKYfa9G7gbIC8vzyfxBjJXcde1NnRmHMmKmSk8u7OagUGHXw1D4i0VjZ2E\nCcxM958xvIZTlJXImzqz5kfUtZ3ha0/tYs/J03x+VR73XjmP1Lioj2yzeuF0/tdV89hf08Yv3yzn\nZ1uPsnlfHf/vs0tYnJNsU+SeYed/Zw2Q6/Y6x1r2ESJyBfA9YI0x5mwPN2NMjfXzGPAmsHS4kxhj\n1htjSowxJRkZGZ6LPkht2lvrN8VdLstnpNDdNxgyc5lXNHWRl+ofszSeS2FmIqc6e2ls77E7FL9Q\neaqL6x95l6MNHfzy1mX88PpFH0sm7hZmJ/HLW5fz+O0raD8zwE2/ep9nSk+OuH0gsDOh7ADmiEi+\niEQB64CPtNYSkaXAf+JMJo1uy1NEJNp6ng5ciFvdi5qYQYfhpf31flPc5bIsLwWAPSdDo9ir3M9b\neLm4hrLXDo5wqrOXzz+6nf5Bw7Nfu2Bc/bc+OX8qL3/zIlbmp/KPz+7lgRcOBGxzbNsSijFmALgH\n2AIcBJ4xxpSJyAMi4mq19VMgHvjTkObBC4BSEfkQeANnHYomlEkqPd5ite7yn+IucA6Znh4fze4Q\naOk16DBUnury6/oTl8KzCSW0W3r1DTj42lM7ae7qZcMdK1mQmTjuYyRPieKJO1Zw+wUzeezdSr79\n7F4GAnDwTVu/hhpjNgObhyz7vtvzK0bY7z1gkXejCz2u1l2X+VFxFziHTC/OTWZPVfAnlJMt3fQN\nOpgVAHcoCTGRzEibEvJ3KP/8Qhk7jrfy77csZVFO0oSPExEexv2fLiQ1Lop/e/UIPf2D/HxdcUDV\nG/pPuYay1aDDsHl/PZfN96/iLpeleclsPdjA6e4+kqeMXC4d6AKhhZe7oqxE9teEbkJ541Ajv9te\nxVcuLvDIFNkiwjcun0NsZDgPbj5IbFQ4P7lxsd+OmDBU4KQ+5VWlx1to6rB37K5zWZrrbP0S7B0c\nXUPCzw6YhJJEVUs37UHcWW8k7T39/NNz+5g7LZ57r5rr0WN/+eICvnnFHJ7dWc0DLx7AmMCoU9GE\nogDYtK+OmEj/K+5yWZybjAjsDvJir4qmTtLjo0maEhid3Qqt+oKDIVjs9dBLh2ho7+EnNy3xSou8\nv798Dl+6MJ8n3jvOz7Ye9fjxvUETivLb1l3u4qMjmDs1ISTuUGZPjbM7jDEL1ZZe+2va+P32Kr50\nYT7Fud7pOyIi/O/rFnDT8hx+/tpRntnh/02KNaEotlc20+SHrbuGWpqXzJ6TpwPm9n+8jDFUNAVG\nCy+XqYkxpMdHh1RCMcbwfzcfJDUuim9cMcer5woLE350wyIumpPOPz23j3eO+veU2JpQFC/urSM2\nMtxvi7tcluYl03amn8pTXXaH4hWnOvtoO9MfEH1Q3BVlhdbcKG8ebuK9ima+cdlsn4zDFRkexiO3\nLmNWRjxfe2onh+v9t4OvJpQQNzDo4OX99Vy+YCpTovyzuMvFNSzFviCdejbQWni5FGUlUt7YSe/A\noN2heN2gw/Cjlw6Snx7H35w3Y/QdPCQxJpLH7lhBTFQ4X3pih9+OTqAJJcS9f6yZlq6+gBgxds7U\neKIjwthbHZwJ5WwLr4C7Q0liwGE4Ut9pdyhe9/L+eo40dHLvlXOJivDtx2d2ciyP376C1u4+7txQ\nSnef/80RpAklxL34YR3x0RFcOs//xzmLCA9jQWZiUN+hTIkKJzMpxu5QxsVVMX+gLjj/Li7GGH7x\nRjkFGXG2Na9fmJ3Ef9yylLLaNr7xhz0M+tkQLZpQQljfgIOXy+q5snAaMZH+PRChy6LspKCdera8\nsZNZGfGIBEYnNpe81CnER0cEfcX864caOVjXzt9dOptwGzsaXr5gGvd/uoitBxv44Sb/GnFKE0oI\ne7f8FG1n+vmUn7fucrcoJ4nO3gEqm4OvYv5YUxezMgKnybBLWJiwIDMhqBOKMYb/eL2c3NRY1hTb\nXzx82wUz+dKF+Tz+7nGeeLfS7nDO0oQSwl7YW0tCTAQXzfH/4i6XRdnOsZL2B1mxV1evc9rfQKs/\ncSnKSuJgXbvfFcF4yq6qVvacPM3dF88i0k/G1vredQu4snAaD7x4gK0HGuwOB9CEErJ6+gd5tayB\nq4um+7xycTJcFfP7gqxi3tUUOtBaeLkUZiXS3TfI8SC8cwR4/N3jJMZEcOOybLtDOSs8TPj5umIW\nZifx9T/s9ov/icD5JFEe9daRJjp6BwKquAv+WjG/N8juUAK1hZdLMPeYr2/r4aX99dy8ItfvmtZP\niYrg0dtKSI2L4s4NO6g9fcbWeDShhKgX99aRMiWSC2en2x3KuAVjxXxFUyfhYUJemn9P+zuSOVMT\niAyXoOzg+LvtJ3AYwxdWzbQ7lGFNTYjh8TtWcKZvkDse30Fbt30DdWpCCUFdvQNsPdjA6oXT/aY8\neDwWZQdfxXxFU2dATPs7kqiIMOZMTeBAkN2h9A4M8vvtVVw+f5pfJ/u50xL49ReWU3mqi9uf+ICu\nXnv6qATep4matC1l9XT3DfKZpTl2hzIhC4OwYr6iMbDG8BpOUVYiB2rbg2qstdcONtLc1cfnV+XZ\nHcqoLpydzr/fspS91W18+bel9PT7fuQCWxOKiKwWkcMiUi4i9w2zPlpE/mit3y4iM93WfddaflhE\nrvZl3IHuz7uqyU2NZcXMFLtDmZA504KrYv7stL8BNMrwcIqyEmnu6qOhvdfuUDzmjztOkpkUEzAt\nIVcvnM5Pb1rMexXN3PP7XT4fDmfEhCIim90/wD1NRMKBR4BrgELgFhEpHLLZnUCrMWY28DDwY2vf\nQmAdUASsBn5pHU+Novb0Gd6raOaGpTkB14HOJTLIesxXt1rT/gb6HYp15xgs9Si1p8/w1tEmblqe\nY2tHxvG6YVkO/3L9QrYebOTLv93JmT7fJZVz3aE8DrwiIt8TEW8MqbkSKDfGHDPG9AFPA2uHbLMW\n2GA9fxa4XJyfgmuBp40xvcaYSqDcOp4axXO7azAGblwWmMVdLouykygLkor5QB0UcqgFmYmIBE9L\nr//eVY0x8NnluXaHMm5fWDWDH9+4iLePNvHFx7b7bEbNEROKMeZPwDIgESgVkW+JyL2uhwfOnQ24\nzxhTbS0bdhtjzADQBqSNcV+PefTtY/x0y6GALxs2xvDnndWsnJnq1xWMYxFMFfOuJsOB2EveXXx0\nBDPT4oLiDsXhMDxTWs35BWkB+79y84o8/uOWpew5eZpb1m/jVKf3iyJHq0PpA7qAaCBhyCMgiMjd\nIlIqIqVNTU3j3t8Yw7FTXTzyRgX/9Nz+gO4JvPvkaY6d6uLG5f7TOWuiFp4tXgn8b8MVjV2kx0eR\nPCXK7lAmrTAzkQN1gf832V7ZQlVLN59bEdh38p9anMX6L5bQ0TNAZ4/3W36N2EtHRFYD/wZsBJYZ\nY7o9fO4awP1eMsdaNtw21SISASQBzWPcFwBjzHpgPUBJScm4s4GI8OD1C0mOjeSXb1bQ3tPPz28u\nJiIAm9v+eWc1MZFhto2U6kmzp8af7fewZon9YytNRkVTJwUBXtzlUpiVyKZ9dbSd6Scp1vuTT3nL\nn0pPkhATwTULA/9/5ZPzpnLhvek+GRHjXGf4HvBZY8x9XkgmADuAOSKSLyJROCvZNw7ZZiNwm/X8\nJuB14yx32giss1qB5QNzgA+8ECPgTCr/uHo+371mPpv21nHvMx8G3J1KT/8gL3xYy9VF00nwwSxz\n3hYVEcbcacHR76GiqTPg609czg5lH8B/l46efjbvr2PNkqyAGYV7NL4aXmnEOxRjzEXePLExZkBE\n7gG2AOHAY8aYMhF5ACg1xmwEfgM8KSLlQAvOpIO13TPAAWAA+DtjjNebMnzlklk4DPz45UNEhAk/\n/eySgGn98eLeOtp7Brh5ReBVMI6kKCuRrQcbMcYEbIu1lq4+WrsDb9rfkRRl/bWl1/mz0myOZmJe\nO9hIT7+DzywN/KJhX7N1YBpjzGZg85Bl33d73gN8doR9HwQe9GqAw/japbMYGHTwr68eITI8jIdu\nXBQQH2ZPbTvBrIw4zi8IzH/y4RRlJfFMaTX17T1kJsXaHc6E/LWFV2BXyLtkJEQzNSE6oO9QXtxb\nR2ZSDMvyArOflp0CryLAD3z98jl8/bLZ/LH0JP/6yhG7wxnV/po29pw8za3nzQiI5DdWZwckrAnc\nD6+KxuBoMuyuKCsxYBtLtJ3p560jTVy3KJOwACl98CeaUCbo3ivnsm5FLr94o5yntp2wO5xzemrb\nCWIiw7hxeWC3WBlqvtXvIZBbFVU0dRIdEUZ2cmDeYQ2nKCuJ8qZOW4b+mKytBxroG3RwXYCNwu0v\nNKFMkIjww+sXctn8qXz/+f28UlZvd0jDOt3dx/N7alm7JDugW90MJxj6PZQ3Olt4BdO34cKsRAYd\nhiMNHXaHMm6b9tWRnRxLcW6y3aEEJE0okxARHsYv/mYpi6wJbnaeaLU7pI/53fYqzvQPcvuFM+0O\nxSsKA7h4BaAiQKf9PZdAnRulrbuft482cd3izKAqGvYlTSiTNCUqgt/cvoLMpBju3LDj7Mx7/qBv\nwMGG945z0Zx0FmQm2h2OVxRlJVLdesbWOSAmqqd/kJOt3UFVfwKQmzKFhOiIgLtz3HKgnv5BE3CT\nzvkTTSgekB4fzRN3rCRMhNsf/4BmHwxxMBYbP6ylsaOXL19UYHcoXnO2mWpdYH14ARxv7sIYmBUk\nTYZdwsLrkfafAAAdgElEQVSEwqxE9gVYY4lNe+vITY1lkTUKgxo/TSgeMjM9jv/6Ygl1bT22zUXg\nzhjDo28fY/70BC6aE3izMo5VYWbgdqSraHTNIx9cRV4Ai3OSOFjXTt+Aw+5QxqS1q493y09x3aIs\nLe6aBE0oHrR8Rgo/u7mY3SdPc+8ze2wdCfeVAw0cqu/gyxcVBPU/SCD3e6ho6kQECtKD6w4FYElu\nMn0DDg7XB0bF/JayegYcWtw1WZpQPOzaRZn80zUL2LyvnodePmRLDA6H4d9eOUJBehxriwN7nKux\nCNR+DxVNnWQnxxIbFRzDe7hbkuNsJfVh9WmbIxmbTfvqmJk25WyDAjUxmlC84K6L8vnCqhmsf+sY\nT75/3Ofn37SvjsMNHfz9FXMCchDL8QrUfg/BNIbXUDkpsaRMiWRvACSU5s5e3qto1tZdHhD8nzY2\nEBHu/3Qhl8+fyv0by3jtYIPPzj0w6OBnW48wd1o8n14c/Hcn4LxDGXSYgCleAeddZDDMIz8SEWFx\nTjJ7A2Ca5pfL6hl0GK5bFBr/L96kCcVLIsLD+PdbllKYlcg9v9/ts/nPn9x2goqmLv7XVfOCqrPc\nufx1QMLAKfaqa+/hTP9gwM8jfy5LcpM50tBBd5/35+GYjE176yhIj2NBZsBM8+S3NKF4UVx0BI/d\ntoLUuCi+tGEH1a3emAXgrxrbe/i3V45w8dwMriqc5tVz+ZPc1FgSoiM4EEBNh4NxDK+hluQk4TCw\n34+bDzd19LLtWDOf0uIuj9CE4mVTE2N4/I4V9PQP8qUndtB2xnsd8H700iF6Bxz885qikPrnEBEW\nBFjFfLDMI38ui62KeX+uR3m5rB6HgetCpHjY2zSh+MDcaQn85+eXc6ypi689tdMrlcevHWzgud01\n3H1xAfnpwVuMMpKirEQO1XUEzMRnFU2dJMZEkB4f+NP+jiQjIZrs5Fj2nPTfhPLih7XMnhrP3GnB\nm9h9SROKj1wwO52f3LSY9481c9eGUs70eS6pNLb38O1n97IgM5F7LpvtseMGkqKsJM70D1J5qtPu\nUMakorGL2VPjg/5Ocklukt8mlMb2Hj443qLFXR6kCcWHbliWw09vWsJ7Fae4/fEP6OqdfGXloMNw\n7zMf0t03wH/cUhw0U5aOV6ANSBjMTYbdLctLobr1DI3tPXaH8jEv7a/HGLhukXZm9BRbEoqIpIrI\nqyJy1Pr5sanRRKRYRN4XkTIR2SsiN7ute0JEKkVkj/Uo9u1vMHE3Lc/h4ZuLKT3Ryhd+s52Wrr4J\nH8sYw/95fj/vlJ/iB58uYvbU0G2lMntqPFHhYQGRUNp7+mns6A26MbyGs2yG8197V5X/jcT94t5a\n5k1LYM600P2/8TS77lDuA14zxswBXrNeD9UNfNEYUwSsBn4mIu6TFHzbGFNsPfZ4P2TPWVuczS9u\nWcr+2nbW/OIdDk5wgqh/feUIv99exdcuncW6lXkejjKwRIaHMXd6fEAMwRIKLbxcirISiQoPY1eV\nfxV71bf1sON4q06k5WF2JZS1wAbr+Qbg+qEbGGOOGGOOWs9rgUYgw2cRetk1izJ55ivn0z/o4IZf\nvscT71aOeeyvgUEHD7xwgF+8Uc66Fbn849XzvBxtYCjKTKKstg1j/LtivqIpeAeFHCo6IpyF2Yns\n8rO5gjbvqwPQhOJhdiWUacaYOut5PXDOThMishKIAircFj9oFYU9LCLRXorTq4pzk3nhnk+wMj+V\nH7xwgJvXv8/2Y83n/EA8XN/Bzeu38di7ldx+wUwe/MwirVC0FGUn0trdT12b/5XXu6to6iQyXMhN\nnWJ3KD6xfEYKe2va/Grk4U376liQmRgSd4m+FOGtA4vIVmD6MKu+5/7CGGNEZMRPUBHJBJ4EbjPG\nuN6R38WZiKKA9cB3gAdG2P9u4G6AvDz/KxaamhjDE3es4Nmd1fz45UPcvH4bS3KTuapwGqsKUpma\nEEN33yAH6tp4aV89Ww82kBATyc/XFbO2ONvu8P2Ke8V8lh/P0X60oZOZaXFEhsA4a+CsmP+vtysp\nq21jad7Hqkt9rvb0GXaeaOXbemfvcV5LKMaYK0ZaJyINIpJpjKmzEkbjCNslApuA7xljtrkd23V3\n0ysijwPfOkcc63EmHUpKSvyyLERE+GxJLp9eksWfdlbz++1V/HTL4Y9tlx4fxVcvmcWXLyogJS54\n+y9M1PzpiYhAWW0bV/rxSAHljR1nh4sJBa6K+Z0nWv0ioZwt7tLWXR7ntYQyio3AbcBD1s/nh24g\nIlHAc8BvjTHPDlnnSkaCs/5lv/dD9r6YyHC+sGoGX1g1g1OdveyuOk1rdx/REWHMyoinMDMxZMbn\nmoi46Ajy0+P8uqVXT/8gJ1q6Q+ruclpiDNnJsez2k4r5F/fWsTA7kZkh2AHY2+xKKA8Bz4jIncAJ\n4HMAIlICfNUYc5e17GIgTURut/a73WrR9TsRyQAE2AN81cfxe116fLRff8v2V4WZiX7zwTWciqZO\njHGOnhBKSmam8F6Fs37Qzjq/ky3d7Dl5mu+snm9bDMHMloRijGkGLh9meSlwl/X8KeCpEfa/zKsB\nqoBVlJXEi3vrON3dR/IU/ysWPNrgbDI8J8SG+jgvP43n99RSeaqLAhsrwrW4y7tCo1ZQhQxXxby/\n9kc50tBBRJgwMy20iltWFaQCsO1Yi61xbNpXx+KcJPLSQqOFna9pQlFBxd+HYDna2El+ehxREaH1\nr5efHsfUhGi2HWu2LYaq5m72VrfpvPFeFFrvahX00uKjmZ4YQ1mtf86NcrShI+SKu8DZkvG8gjS2\njdLPypte3FcLwLVa3OU1mlBU0Cn007lRevoHqWrpZk6Ijrm2qiCVxo5eKk912XL+TXvrKM5NJidF\ni7u8RROKCjpFWYlUNHV6Zd6Zyaho6sRhQq9C3mVVQRoA2yt9X49SeaqLstp2Le7yMk0oKugUZSXi\nMHCovsPuUD7C1cIr1JoMuxSkx5Eeb089yvN7ahDRsbu8TROKCjquXuj+Vo9ytDE0W3i5iAjnz0o7\n2x/FV4wxPL+nllX5aWQm+e+QPMFAE4oKOjkpsSTGRPhdPcqRhk5mhmALL3cXzUmnqaOXg3W+u3v8\nsLqNylNdXL9U5433ttB9Z6ugJSJ+WTFf3tgZ8nOXXzLXOQPFX440+eyc/7O7hqiIMFYv1OIub9OE\nooJSYWYSh+raGRj0jyHTe/oHOdHcFdKzaoJzXK/50xN4y0cJZWDQwYt7a7l8/lSSYiN9cs5QpglF\nBaWirER6Bxy2NVEdytXCK9TvUMB5l1J6ooWu3gGvn+ud8lOc6uzj+qWhMxinnTShqKBUlO1fPebL\nG0O7hZe7i+dm0D9oeL/C+629nt9TS2JMBJfOC5rJXv2aJhQVlGZlxBMVEeY3Lb1CdQyv4ZTMTCE2\nMpy3jnq32Ku7b4AtZfVctziL6Ihwr55LOWlCUUEpMjyM+dMT/OYORVt4/VV0RDjnz0rjjcONXm0+\n/NK+err7Brm+WFt3+Yq+u1XQKsx0tvSya+wod4fq25k/XYu7XK4qnMbJljNeTfhP76giPz2Olfmp\nXjuH+ihNKCpoFWUl0namn5rTZ2yNo6Onn5MtZ1iQmWhrHP7kysJphAm8vL/eK8cvb+xgx/FW1q3I\ntXVCr1CjCUUFrUKrx7zdc6McaXB24tM7lL9Ki49mVUEam/fXeeUO8g8fnCQyXLhxeY7Hj61GpglF\nBa0FmQmI2N/Sy9UrfL7eoXzENYsyOdbUxVGrBZyn9PQP8t+7qrmqcDrp8dEePbY6N1sSioikisir\nInLU+pkywnaDIrLHemx0W54vIttFpFxE/igi/jfXq7LdlKgICtLjbE8oh+rbSYiJICspxtY4/M3V\nRdMQcVaee9JL++to7e5n3cpcjx5Xjc6uO5T7gNeMMXOA16zXwzljjCm2Hmvclv8YeNgYMxtoBe70\nbrgqUBVlJXHA5qbDh+o6WDA9Ucvyh5iaEMOKGals2lfrsWIvYwy/eaeSWRlxXDgr3SPHVGNnV0JZ\nC2ywnm8Arh/rjuL8r7wMeHYi+6vQUpiVSG1bD61dfbac3xjDofoO5mdq/clwrl+azZGGTvacPO2R\n431Q2cL+mna+9Il8wsI0gfuaXQllmjGmznpeD0wbYbsYESkVkW0i4koaacBpY4xr3IZqYMRxFUTk\nbusYpU1NvhuQTvkH1xzzB+rsKfaqbj1DZ+8A86dr/clw1hRnMSUqnD98UOWR4/3mnUpSpkRyw1Kt\njLeD1xKKiGwVkf3DPNa6b2ec97oj3e/OMMaUAH8D/ExEZo03DmPMemNMiTGmJCNDh18INQutll57\nq+0p9nJN8qV3KMOLj45gbXEWL3xYR3tP/6SOdaK5i1cPNnDreTOIjdKe8XbwWkIxxlxhjFk4zON5\noEFEMgGsn40jHKPG+nkMeBNYCjQDySISYW2WA9R46/dQgS0lLooZaVPYc7LVlvMfsu6M5ukYXiO6\nZWUeZ/oHeX735P6Nf/F6OZHhYXzx/BkeikyNl11FXhuB26zntwHPD91ARFJEJNp6ng5cCByw7mje\nAG461/5KuRTnJnusjH68Dta3k5c6hbjoiNE3DlGLspMoykrkd9urJlw5X97YyZ93VfOFVTOYmqit\n6exiV0J5CLhSRI4CV1ivEZESEXnU2mYBUCoiH+JMIA8ZYw5Y674D3Csi5TjrVH7j0+hVQCnOTaah\nvZe6Nt/3mN9f087CbK0/ORcR4bbzZ3KovmPCE289vPUIsZHh/O2l4y4VVx5ky9cmY0wzcPkwy0uB\nu6zn7wGLRtj/GLDSmzGq4FGcmwzA7qrTZC7y3Zzibd39VLV0c/MK7Q8xmuuXZvPw1iM88kY5l8zN\nGFcT6/01bWzaW8fXL5tNmnZktJX2lFdBrzArkajwMJ8Xe7mGzl+UneTT8waiqIgw/vbSWew43srr\nh4atUh2Ww2H4P8/vJzUuirsuKvBihGosNKGooBcdEU5hViJ7qnybUPbVOBPKQk0oY7JuZR4F6XE8\nuPkgPf2DY9rnt+8fZ3fVab537QKd4tcPaEJRIaE4N5l9NW0+nWN+f2072cmxpMbpyEBjERkexv1r\nijjW1MXPXzs66vaH6tv5vy8d4pPzMrhhmU7x6w80oaiQsDQvmTP9gxy2Rv71hbKatrMdK9XYXDI3\ng5tLcvn1Xyp44/DIRV/Nnb18+belJMVG8pObluiwNn5CE4oKCa6KeV/Vo3T09HPsVJfWn0zA/WsK\nmT89kXt+t2vYeedrTp9h3fptNLb3sv4Ly8lI0Ip4f6EJRYWEvNQppMZFsdtH9SiuOVi0/mT8pkRF\n8PjtK8hKjuULv9nOg5sOcKi+nWNNnTz69jFW/+wt6tp6eOKOlSzNG3agcmUT7W2lQoKIsCQnyWd3\nKPuthFKkfVAmZHpSDM9+9QJ+uOkAj75TyX+9XXl23QWz0vjRDYuYkRZnY4RqOJpQVMgozk3hzSNN\ntPf0kxjj3RZB+2vamJoQzdQE7bU9UUlTIvnpZ5fwD1fOZcfxFvoHDYuyk5inM1/6LU0oKmQszUvG\nGNh7so1PzPHuXBm7q1rP1tuoyclKjmVtsbbiCgRah6JCxpLcZERg5wnvDhTZ3NnL8eZuls3Q8n0V\nWjShqJCRFBvJvGkJfHD84y2HPGmXVfG/TCuMVYjRhKJCyqqCNHaeaKVvwHsdHHdVtRIRJizO0RZe\nKrRoQlEhZWV+Kj39DvZ7cZ75XSdaKcxKJCZSJ3lSoUUTigopK/NTAdh+rMUrx+8fdPBh9Wkt7lIh\nSROKCinp8dHMyojjg0rv1KMcquugp9+hFfIqJGlCUSHnvII0So+3MuiY2OyA57KrytmCbFmeNhlW\noUcTigo55+Wn0tE7wEFrvndP2lXVytSEaLKTfTeRl1L+wpaEIiKpIvKqiBy1fn6sfEBEPikie9we\nPSJyvbXuCRGpdFtX7PvfQgWq8/LTANh2zLPFXsYYth1rZkV+qo5+q0KSXXco9wGvGWPmAK9Zrz/C\nGPOGMabYGFMMXAZ0A6+4bfJt13pjzB6fRK2CwvSkGGZlxPH20VMePW5FUxcN7b1cOMu7vfCV8ld2\nJZS1wAbr+Qbg+lG2vwl4yRjT7dWoVMi4eG4G2yubxzwz4Fi8X+FMUBfMSvPYMZUKJHYllGnGmDrr\neT0wbZTt1wF/GLLsQRHZKyIPi8iIEyKIyN0iUioipU1NTZMIWQWTi+dm0NPv4INKzzUffq+imezk\nWGakTfHYMZUKJF5LKCKyVUT2D/NY676dMcYAIza3EZFMYBGwxW3xd4H5wAogFfjOSPsbY9YbY0qM\nMSUZGRmT+ZVUEFmVn0ZURBhvHfHMlwyHw/D+sWbOn5Wm9ScqZHlttGFjzBUjrRORBhHJNMbUWQlj\n5Lk+4XPAc8aYfrdju+5uekXkceBbHglahYzYqHBWzkzlraOeSSgH6to53d3PhbO1uEuFLruKvDYC\nt1nPbwOeP8e2tzCkuMtKQojzq+D1wH4vxKiC3MVz0znS0Eld25lJH8s1Ve0FWiGvQphdCeUh4EoR\nOQpcYb1GREpE5FHXRiIyE8gF/jJk/9+JyD5gH5AO/NAHMasgc/FcZxGoJ4q93ik/RUFGHNMSdUIt\nFbpsmWDLGNMMXD7M8lLgLrfXx4GPzaxjjLnMm/Gp0DBvWgLZybG8UtbAzSvyJnycjp5+3q9o5ovn\nz/BgdEoFHu0pr0KWiHDtoum8dbSJtjP9o+8wgjcON9E36ODqhdM9GJ1SgUcTigpp1y7KpH/QsPVA\nw4SPsaWsnvT4aB1hWIU8TSgqpBXnJpOdHMumfXWjbzyMnv5B3jzUyJWF0wgP0+bCKrRpQlEhzVXs\n9fYEi73eLT9FV98gVxeN1jdXqeCnCUWFvOsWZ9E/aHilrH7c+24pqychOkKbCyuFJhSlWJKTxIy0\nKfxpZ/W49jvTN8hL++q5snAaURH6r6SU/heokCci3LIyjw8qWyhv7Bjzfpv31dHRO8DNK3K9GJ1S\ngUMTilLAZ5fnEBku/Pb9E2Pe53fbT1CQHnd2nnqlQp0mFKWAtPhoPrM0m2dKT9Lc2Tvq9qXHW9hV\ndZovnj9DB4NUyqIJRSnL3RfPonfAwaPvVI667a//UkHylEg+p8VdSp2lCUUpy+yp8axZksVj71Se\nc8DIHcdb2HqwkTsvzGdKlC2jFynllzShKOXmW1fNwwA/2FiGc6qej+obcPDPL5QxLTGauy4q8H2A\nSvkxTShKuclNncK3rprLlrKGYSvof/LyIfbXtPPPa4qIjQq3IUKl/Jferys1xF2fKGD7sRZ+8EIZ\nPf2D3HVRAf2DDv71lcM8+k4lXzx/BqsXZtodplJ+R4a7rQ9WJSUlprS01O4wVADo6R/k75/ezZay\nBhJjIhhwGLr7Bvmb8/L4l7ULddwuFVJEZKcxpmS07fQORalhxESG8+vPL+fVAw28eaSJqPAwri6a\nzvmzdIpfpUaiCUWpEYgIVxVN56oinedEqbGwpVJeRD4rImUi4hCREW+jRGS1iBwWkXIRuc9teb6I\nbLeW/1FEonwTuVJKqZHY1cprP3AD8NZIG4hIOPAIcA1QCNwiIoXW6h8DDxtjZgOtwJ3eDVcppdRo\nbEkoxpiDxpjDo2y2Eig3xhwzxvQBTwNrxTnOxWXAs9Z2G4DrvRetUkqpsfDnfijZwEm319XWsjTg\ntDFmYMjyYYnI3SJSKiKlTU1NXgtWKaVCndcq5UVkKzBcbeb3jDHPe+u8Qxlj1gPrwdls2FfnVUqp\nUOO1hGKMuWKSh6gB3Efey7GWNQPJIhJh3aW4liullLKRPxd57QDmWC26ooB1wEbj7In5BnCTtd1t\ngM/ueJRSSg3PrmbDnxGRauB8YJOIbLGWZ4nIZgDr7uMeYAtwEHjGGFNmHeI7wL0iUo6zTuU3vv4d\nlFJKfVRIDb0iIk3A2Kfk+6h04JQHw/EUjWt8NK7x0bjGJ1jjmmGMyRhto5BKKJMhIqVjGcvG1zSu\n8dG4xkfjGp9Qj8uf61CUUkoFEE0oSimlPEITytittzuAEWhc46NxjY/GNT4hHZfWoSillPIIvUNR\nSinlEZpQ3PjrsPoikioir4rIUetnyjDbfFJE9rg9ekTkemvdEyJS6bau2FdxWdsNup17o9tyO69X\nsYi8b/2994rIzW7rPHq9Rnq/uK2Ptn7/cut6zHRb911r+WERuXoycUwgrntF5IB1fV4TkRlu64b9\nm/oorttFpMnt/He5rbvN+rsfFZHbfBzXw24xHRGR027rvHK9ROQxEWkUkf0jrBcR+Xcr5r0issxt\nneevlTFGH9YDWADMA94ESkbYJhyoAAqAKOBDoNBa9wywznr+a+BrHorrJ8B91vP7gB+Psn0q0AJM\nsV4/Adzkhes1priAzhGW23a9gLnAHOt5FlAHJHv6ep3r/eK2zd8Cv7aerwP+aD0vtLaPBvKt44T7\nMK5Pur2HvuaK61x/Ux/FdTvwi2H2TQWOWT9TrOcpvopryPZfBx7zwfW6GFgG7B9h/bXAS4AAq4Dt\n3rxWeofixvjvsPprreON9bg3AS8ZY7o9dP6RjDeus+y+XsaYI8aYo9bzWqARGLXj1gQM+345R7zP\nApdb12ct8LQxptcYUwmUW8fzSVzGmDfc3kPbcI6b521juV4juRp41RjTYoxpBV4FVtsU1y3AHzx0\n7hEZY97C+eVxJGuB3xqnbTjHQczES9dKE8r4eWRY/XGaZoyps57XA9NG2X4dH38zP2jd8j4sItE+\njitGnFMIbHMVw+FH10tEVuL81lnhtthT12uk98uw21jXow3n9RnLvt6My92dOL/pugz3N/VlXDda\nf59nRcQ1iKxfXC+raDAfeN1tsbeu12hGitsr1yrk5pQXPxlWf6hzxeX+whhjRGTEpnnWt49FOMdA\nc/kuzg/WKJzNB78DPODDuGYYY2pEpAB4XUT24fzQnDAPX68ngduMMQ5r8YSvVzASkc8DJcAlbos/\n9jc1xlQMfwSPewH4gzGmV0S+gvPu7jIfnXss1gHPGmMG3ZbZeb18JuQSivHTYfXPFZeINIhIpjGm\nzvoAbDzHoT4HPGeM6Xc7tuvbeq+IPA58y5dxGWNqrJ/HRORNYCnwZ2y+XiKSCGzC+WVim9uxJ3y9\nhjHS+2W4bapFJAJIwvl+Gsu+3owLEbkCZ5K+xBjT61o+wt/UEx+Qo8ZljGl2e/kozjoz176XDtn3\nTQ/ENKa43KwD/s59gRev12hGitsr10qLvMbPjmH1N1rHG8txP1Z2a32ouuotrgeGbRHijbhEJMVV\nZCQi6cCFwAG7r5f1t3sOZ/nys0PWefJ6Dft+OUe8NwGvW9dnI7BOnK3A8oE5wAeTiGVccYnIUuA/\ngTXGmEa35cP+TX0YV6bbyzU4RyMH5135VVZ8KcBVfPRO3atxWbHNx1nJ/b7bMm9er9FsBL5otfZa\nBbRZX5i8c6082eIg0B/AZ3CWJfYCDcAWa3kWsNltu2uBIzi/YXzPbXkBzn/4cuBPQLSH4koDXgOO\nAluBVGt5CfCo23YzcX7zCBuy/+vAPpwfjE8B8b6KC7jAOveH1s87/eF6AZ8H+oE9bo9ib1yv4d4v\nOIvQ1ljPY6zfv9y6HgVu+37P2u8wcI2H3++jxbXV+j9wXZ+No/1NfRTXj4Ay6/xvAPPd9v2SdR3L\ngTt8GZf1+gfAQ0P289r1wvnlsc56L1fjrOv6KvBVa70Aj1gx78Ot9ao3rpX2lFdKKeURWuSllFLK\nIzShKKWU8ghNKEoppTxCE4pSSimP0ISilFLKIzShKGUTEckV56jGqdbrFOv1THsjU2piNKEoZRNj\nzEngV8BD1qKHgPXGmOO2BaXUJGg/FKVsJCKRwE7gMeDLODtX9p97L6X8U8iN5aWUPzHG9IvIt4GX\ngas0mahApkVeStnvGpzDZyy0OxClJkMTilI2Euf0wlfinE3vH4YMfKhUQNGEopRNrNGMfwV80xhT\nBfwU+H/2RqXUxGlCUco+XwaqjDGvWq9/CSwQkUvOsY9SfktbeSmllPIIvUNRSinlEZpQlFJKeYQm\nFKWUUh6hCUUppZRHaEJRSinlEZpQlFJKeYQmFKWUUh6hCUUppZRH/H9DPxgX5SWeqAAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1116a8da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## STAGE 3: plotting\n",
    "\n",
    "# Convert to array\n",
    "from numpy import *\n",
    "xarray = array(xlist)\n",
    "yarray = array(ylist)\n",
    "\n",
    "# plot\n",
    "from pylab import *\n",
    "plot(xlist, ylist)\n",
    "xlabel(\"X\")\n",
    "ylabel(\"Y\")\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Read a data file**</br>\n",
    "The files data/density_water.dat and data/density_air.dat contain data about the density of water and air (respectively) for different temperatures. The data files have some comment lines starting with # and some lines are blank. The rest of the lines contain density data: the temperature in the first column and the corresponding density in the second column. The goal of this exercise is to read the data in such a file and plot the density versus the temperature as distinct (small) circles for each data point. Apply the program to both files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-10.0, -5.0, 0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0]\n"
     ]
    }
   ],
   "source": [
    "## STAGE 1: work with one file\n",
    "# precede exactly as before but with modified for loop:\n",
    "#  line[0] -- refers to the first character in a string\n",
    "#  != -- boolean operator for 'not equal to'\n",
    "#  and -- boolean operator, both T/F statements must be met\n",
    "#  commented lines will start with a #, blank lines will consist of just a newline character '\\n'\n",
    "infile = open(\"data/density_air.dat\", \"r\")\n",
    "temp = []\n",
    "dens = []\n",
    "for line in infile:\n",
    "    try:\n",
    "        t, d = line.split()\n",
    "        t = float(t)\n",
    "        d = float(d)\n",
    "    except:\n",
    "        continue\n",
    "    temp.append(t) # N.B. we're now filling out temp and dens lists\n",
    "    dens.append(d)\n",
    "infile.close()\n",
    "\n",
    "\n",
    "# Sanity check:\n",
    "print(temp)\n",
    "\n",
    "# ... double check the original file to make sure this is right."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-10.0, -5.0, 0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0]\n"
     ]
    }
   ],
   "source": [
    "## STAGE 2: \n",
    "# Let's write it up as a function\n",
    "# We'll use exactly the same code as before\n",
    "\n",
    "\n",
    "def readTempDenFile(filename = \"data/density_air.dat\"):\n",
    "    infile = open(filename, \"r\")\n",
    "    temp = []\n",
    "    dens = []\n",
    "    for line in infile:\n",
    "        try:\n",
    "            t, d = line.split()\n",
    "            t = float(t)\n",
    "            d = float(d)\n",
    "        except:\n",
    "            continue\n",
    "        temp.append(t) # N.B. we're now filling out temp and dens lists\n",
    "        dens.append(d)\n",
    "    infile.close()\n",
    "    return temp,dens\n",
    "            \n",
    "# function defined, let's run it\n",
    "temp,dens = readTempDenFile() # unpack into two objects temp and dens\n",
    "print(temp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAELCAYAAADz6wBxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VVX69vHvnUYXpAlSpArSlQgi3YKIKCgoVmzIKBYU\nnRnH+dl1nHGUYkVAFBwrKqIIFhghNJGAdES6gkiVKp3n/ePszGR4SQgkJ/skeT7XdS5P1i65cww8\n7L3WXktmhnPOOXei4sIO4JxzLm/zQuKccy5bvJA455zLFi8kzjnnssULiXPOuWzxQuKccy5bolpI\nJA2XtFHSwgy2d5E0X9JcSamSWqXbdihonyvp03Tt1SXNlLRc0vuSkqL5MzjnnMucovkciaQ2wC5g\npJk1OMr24sBuMzNJjYAPzKxusG2XmRU/yjEfAB+b2XuSBgPzzOzVqP0QzjnnMhXVKxIzSwG2ZrJ9\nl/23khUDMq1qkgScB3wYNI0AuuZAVOeccyco9D4SSZdL+gH4HLgl3abCwe2ubyWlFYsywDYzOxh8\nvRaolItxnXPOHSEh7ABmNhoYHdwGexK4INh0mpmtk1QD+LekBcD2rJ5XUm+gN0CxYsWa1q1bN4eT\nO+dc/jZ79uzNZlbuWPuFXkjSmFmKpBqSyprZZjNbF7SvlDQJOBP4CCglKSG4KqkMrMvgfEOAIQDJ\nycmWmpqaKz+Hc87lF5LWZGW/UG9tSaoV9Hsg6SygELBF0smSCgXtZYGWwOKgP+UboHtwihuBMbmf\n3DnnXJqoXpFIehdoB5SVtBZ4FEgEMLPBQDegp6QDwB6gRzCC6wzgNUmHiRS7v5vZ4uC0fwbek/QU\n8D3wejR/Buecc5mL6vDfWOG3tpxz7vhJmm1mycfaL/RRW8455/I2LyTOOeeyxQuJc865bPFC4pxz\nLlu8kGRi1uqtDJuyksOH8/+ABOecO1FeSDLx6dxfeOrzJdw6YhZbdu0LO45zzsUkLySZeKJLfZ7o\nUp9py7dw8aApTF+xOexIzjkXc7yQZEISPVtUY/Sd51K8cALXDZtJ/6+WcvDQ4bCjOedczPBCkgX1\nTy3JZ3e1ottZlXnh38u5Zui3/LJtT9ixnHMuJnghyaJihRJ47srGDOzRhMW/7ODiQVP4atGvYcdy\nzrnQeSE5Tl3PrMTYe1pTpXQRer81m0fHLGTvgUNhx3LOudB4ITkB1csW46M7zuXWVtUZMWMNl78y\nnRWbdoUdyznnQuGF5AQVSojn4c71GH5TMr9u38OlL07lw9lrKQiTYDrnXHpeSLLpvLqnML5vGxpW\nKskDo+bR74N57Np38NgHOudcPuGFJAdUKFmYd247h/suOJ0xc9fR+YUpLFyX5VWBnXMuT/NCkkPi\n40TfC2rz7m3nsO/gYS5/ZRrDp67yW13OuXzPC0kOa16jDOPuaU3b08vzxNjF9BqRytbd+8OO5Zxz\nUeOFJApOLpbE0J5NeezSekxZtpmLB6Xw7cotYcdyzrmo8EISJZK4qWV1Pu5zLsWSErh26LcM+PpH\nn17FOZfvRK2QSBouaaOkhRls7yJpvqS5klIltQram0iaIWlRsL1HumPelLQqOGaupCbRyp9TGlQq\nyWd3t6LrmZUYNHEZ1w6byfrtPr2Kcy7/iOYVyZtAx0y2TwQam1kT4BZgWND+O9DTzOoHxw+UVCrd\ncX80sybBa24Ucue4YoUS6H9VE56/sjEL123n4kFTmLB4Q9ixnHMuR0StkJhZCrA1k+277L9DmooB\nFrT/aGbLgve/ABuBctHKmZu6Na3M2LtbUalUEXqNTOXxzxax76BPr+Kcy9tC7SORdLmkH4DPiVyV\nHLm9GZAErEjX/HRwy2uApEK5FDXH1ChXnI/7nMvNLavxxrTVXPHKdFb69CrOuTws1EJiZqPNrC7Q\nFXgy/TZJFYG3gJvNLK2H+i9AXeBsoDTw54zOLal30PeSumnTpqjkP1GFEuJ59NL6DO2ZzLpte+j8\n4lQ+nrM27FjOOXdCYmLUVnAbrIaksgCSTiJylfJXM/s23X7rLWIf8AbQLJNzDjGzZDNLLlcuNu+M\nXVjvFMb3bU2DSiXp98E8+n0wl90+vYpzLo8JrZBIqiVJwfuzgELAFklJwGhgpJl9eMQxFYP/ishV\nzFFHhOUlFUsW4d3bzqHv+bX55Pt1dH5xqk+v4pzLU6I5/PddYAZQR9JaSbdKul3S7cEu3YCFkuYC\nLwM9gs73q4A2wE1HGeb7tqQFwAKgLPBUtPLnpvg4cd+Fp/N2r3P4ff9BrnhlOm9O8+lVnHN5gwrC\nX1bJycmWmpoadows2bp7P38cNY+JP2zkgjNO4Z/dG3FysaSwYznnCiBJs80s+Vj7xUQfifuv0sWS\nGHZjMo90rsfkHzfS6YUpzPTpVZxzMcwLSQySxC2tqvPxHS0plBDHNUO/ZdCEZRw6nP+vHp1zeY8X\nkhjWsHJJxt7Tmi5NKjFgwo9cO/Rbft2+N+xYzjn3P7yQxLjihRIY0KMJz13ZmPlrt3PxoBT+/YNP\nr+Kcix1eSPKI7k0rM/aeVlQoWYRb3oxMr7L3gE+v4pwLnxeSPKRmueKM7nMuN50bmV7lkhem8P1P\nv4UdyzlXwHkhyWMKJ8bz2GX1eevWZuzZf4hur07n2S9+8MkfnXOh8UKSR7WuXY4v7mtD96aVeWXS\nCrq8NM2fiHfOhcILSR52UuFEnu3emOE3JbN19366vjyNgRN+5ICvwuicy0VeSPKB8+qewlf3teHS\nxqcycMIyur48jR9+3RF2LOdcAeGFJJ8oVTSJAT2aMPj6pmzYsZdLX5zKy98s9zXinXNR54Ukn+nY\noAJf3tuGDvUq8M8vl9Jt8AyWb/SFs5xz0eOFJB8qU7wQL193Fi9ecyZrtuym0wtTGJqy0qdYcc5F\nhReSfOzSxqfy1X1taFO7HE+PW0KP12awevPusGM55/IZLyT5XPkShRnasynPX9mYpRt2cvGgKYyY\nvprDfnXinMshXkgKAEl0a1qZr+5rw9nVS/Pop4u4bthMft76e9jRnHP5gBeSAqRiySKMuPls/n5F\nQxas207HgSm8M/MnX4nROZctXkgKGElc3awqX9zbmsZVSvHQ6AXc+MYs1m/fE3Y051we5YWkgKp8\nclH+dWtznuhSn1mrttJhQAqjUn/2qxPn3HGLaiGRNFzSRkkLM9jeRdJ8SXMlpUpqlW7bjZKWBa8b\n07U3lbRA0nJJL0hSNH+G/CwuTvRsUY3xfVtTt0IJ/vjhfG4bmcrGHb54lnMu66J9RfIm0DGT7ROB\nxmbWBLgFGAYgqTTwKNAcaAY8Kunk4JhXgduA2sErs/O7LKhWthjv9W7B/11yBlOWbabDwBTGzF3n\nVyfOuSyJaiExsxRgaybbd9l//7YqBqS9vwj42sy2mtlvwNdAR0kVgZPM7NvguJFA1+j9BAVHfJzo\n1boGn9/TmmplitH3vbn0eXsOW3btCzuacy7Ghd5HIulyST8AnxO5KgGoBPycbre1QVul4P2R7S6H\n1CpfnA9vb8GfOtZh4pKNdBiQwhcL14cdyzkXw0IvJGY22szqErmyeDKnziupd9Dvkrpp06acOm2B\nkBAfR592tfjs7lZULFWY2/81h3ve/Z5tv+8PO5pzLgaFXkjSBLfBakgqC6wDqqTbXDloWxe8P7L9\naOcbYmbJZpZcrly5KKXO3+pUKMHoPi2574LTGbdgPRcOSGHikg1hx3LOxZhQC4mkWmmjriSdBRQC\ntgBfAh0knRx0sncAvjSz9cAOSecEx/UExoQUv0BIjI+j7wW1+eTOlpQplsStI1J5YNQ8duw9EHY0\n51yMSIjmySW9C7QDykpaS2QkViKAmQ0GugE9JR0A9gA9gk70rZKeBGYFp3rCzNI67fsQGQ1WBBgf\nvFyUNahUkjF3teSFict4ddIKpi3fzD+6NaLN6X6151xBp4IwxDM5OdlSU1PDjpFvzP15G/d/MJcV\nm3ZzbfOqPNTpDIoXiuq/SZxzIZA028ySj7VfzPSRuLyjSZVSfH5Pa25rXZ13v/uJjgNTmL5ic9ix\nnHMh8ULiTkjhxHj+ekk9Rv2hBQlx4tqhM3ns00X8vv9g2NGcc7nMC4nLluRqpRnXtzU3nVuNN6ev\nptOgKXz/029hx3LO5SIvJC7biiYl8Nhl9XnntuYcOGR0HzyDQROWcfDQ4bCjOedygRcSl2POrVmW\ncX1b07lRRQZM+JGrXpvBmi2+tK9z+V2mhURSC0kvBzP0bpL0k6Rxku6UVDK3Qrq8o2SRRAZdfSaD\nrm7Cso276DRoCh/49PTO5WsZFhJJ44FeRB4O7AhUBOoB/wcUBsZIuiw3Qrq8p0uTSozv25oGlUry\npw/n0+ftOfy226dYcS4/yvA5EkllzSzTMZ1Z2ScW+HMk4Tl02Bg6ZSXPf7WU0sWSeO7KxrSu7Q8x\nOpcXZPs5kqMViGCdkEz3cS69+Dhxe9uajO7TkhKFE7nh9e944rPF7D1wKOxozrkcktmtrZaSlkha\nJKm5pK+BWZJ+ltQiFzO6fKBBpZJ8dlcrbmxxGsOnraLLS9NYsn5H2LGcczkgs872AcBVRPpJPgce\nN7OaQBfguVzI5vKZIknxPN6lAW/cfDZbdu+ny0vTGDZlJYcPe0e8c3lZZoUk0cwWmNkMYJOZTQUw\nszlEJkx07oS0r1OeL+9tTds65Xjq8yXcMHwm67fvCTuWc+4EZVZI0m/7yxHbkqKQxRUgZYoXYsgN\nTfn7FQ2Zs2YbHQdO4fP5vhKjc3lRZoXkYUlFAczsk7RGSTWJrJXuXLZI4upmVRnXtzXVyhbjznfm\n0O+Duez0tU6cy1N8GnkXEw4cOsyLE5fx0jfLObVUEQb2aEJytdLHPtA5FzU5No28pIY5E8m5jCXG\nx9GvQx1G3d4CCa56bQbPf7WUAz5fl3Mx71hTpJwPvJJLWZyj6WmlGd+3Dd3OqsyL/15Ot1ens3LT\nrrBjOecykdlzJNcB/wSuyL04zkHxQgn888rGvHLdWfy09XcueWEq78z8yefrci5GZXZF8gbQ3cw2\n5VYY59Lr1LAiX/RtQ9PTTuah0Qu4bWQqm3ftCzuWc+4ImRWSR4HXJZ3QMyOShkvaKGlhBtuvC2YV\nXiBpuqTGQXsdSXPTvXZIujfY9pikdem2dTqRbC7vqFCyMCNvacYjneuRsmwzHQem8M0PG8OO5ZxL\nJ7O5tp4BhgOfZLTPMbxJZNbgjKwC2ppZQ+BJYEjwfZeaWRMzawI0BX4HRqc7bkDadjMbd4LZXB4S\nFyduaVWdT+9qSdnihbj5zVk8/MlC9uz3+bqciwWZdrab2VvAsydyYjNLAbZmsn26maWtyfotUPko\nu50PrDCzNSeSweUvdSucxCd3tuS21tV569s1dH5xCgvXbQ87lnMF3jGH/5rZxFzIcSsw/ijtVwPv\nHtF2V3BLbLikk6MfzcWSwonx/PWSerzdqzm79x2i68vTeGXScg75fF3OheaYDyRKigcuAaoBCWnt\nZtb/mCeXqgFjzaxBJvu0JzLEuJWZbUnXngT8AtQ3sw1B2ynAZsCI3A6raGa3ZHDe3kBvgKpVqzZd\ns8YvavKbbb/v56+jF/L5gvU0q16a/lc1pvLJRcOO5Vy+kWMPJAKfATcBZYAS6V7ZJqkRMAzokr6I\nBC4G5qQVEQAz22Bmh8zsMDAUaJbRuc1siJklm1lyuXK+kFJ+VKpoEi9deybPX9mYxb/s4OKBU/jk\n+3Vhx3KuwEk49i5UNrNGOf2NJVUFPgZuMLMfj7LLNRxxW0tSRTNLm9nvcuCoI8JcwSGJbk0r06x6\nae57fy73vj+Xf/+wkSe7NqBkkcSw4zlXIGTlimS8pA7He2JJ7wIzgDqS1kq6VdLtkm4PdnmEyFXO\nK8FQ3tR0xxYDLiRSaNJ7NhguPB9oD9x3vLlc/lSldFHe630OD3Q4nXEL1nPxwBRmrDjyItc5Fw1Z\n6SO5HPgXkaJzABBgZnZS9OPlDJ+0sWCZ9/M27nt/Lqu27KZ3mxr0u/B0CiXEhx3LuTwnJ/tI+gMt\ngKJmdpKZlchLRcQVPI2rlGLsPa24tllVXpu8kstfns6yDTvDjuVcvpWVQvIzsNB8oiOXhxRNSuDp\nyxsyrGcyG3bs5ZIXp/LKpOU+m7BzUZCVzvaVwCRJ44H/THSUleG/zoXtgnqn8EWVNjz66UKe/WIp\nY+et59nujWhQqWTY0ZzLN7JyRbIKmEhked0cHf7rXG4oV6IQr1zXlMHXN2XTrn10eXkafx//A3sP\n+BQrzuWEDK9IJP0F+MLMHs/FPM5FTccGFWhRowx/G7eEwZNX8OWiX3nmioacU6NM2NGcy9MyuyJZ\nCfSV9L2kNyX18ClJXF5Xsmgi/+jeiHd6NefQYePqId/y0OgF7PB14p07YVlas13SmURm8u0AxAMT\niFytfBfdeDnDh/+6o9mz/xADJvzIsCkrKV+iME92bcCF9U4JO5ZzMSMnh/9iZt+b2TNm1h7oDCwC\nemUzo3OhKpIUz0OdzmB0n5aUKprIbSNTueudOb54lnPHKSsPJB5tqd3twAIzyxMrDPkViTuW/QcP\n89rkFbz47+UULRTPI53rcfmZlZAUdjTnQpPVK5KsFJLPiTyQ+E3Q1A6YDVQHngjWLIlpXkhcVi3f\nuJM/f7SA2Wt+o83p5fjb5Q18RmFXYOXkra0E4Awz62Zm3YB6RKZxbw78OXsxnYsttcqXYNQfWvD4\nZfWZvXorHQak8Oa0Vb7eiXOZyEohqZJ+KndgY9C2lcjcW87lK3Fx4sZzq/FVv7Y0q16axz5bzJWD\nfZoV5zKSlUIySdJYSTdKuhEYE7QVA7ZFN55z4alUqghv3HQ2A3o0ZtXm3VzywlQGTVjG/oM+zYpz\n6WWlj0TAFUCroGmamX0Y7WA5yftIXHZt3rWPJz5bzKfzfqHOKSX4R/dGNKlSKuxYzkVVTna2X2xm\n449ou93MBmczY67xQuJyysQlG/jr6IVs3LmXm1tW5/4Op1M0KStT1jmX9+RkZ/vDks5Ld+I/AV2y\nE865vOr8M07h635tuLZ5VV6fuoqLBqYwbfnmsGM5F6qsFJLLgL9Jai3paSKjtbyQuAKrROFEnura\nkPd7n0NiXBzXDZvJH0fNY/vvPvbEFUzHLCRmtplIMXkZOBXobmb7ox3MuVjXvEYZxvVtTZ92Nfn4\n+3VcMGAy4xesDzuWc7kuw0IiaaekHZJ2AMuB04ErgbQ25wq8wonx/KljXcbc2ZLyJQpxx9tz+MNb\nqWzcsTfsaM7lmgwLSdqSuulehc2seFaX2pU0XNJGSQsz2H6dpPmSFkiaLqlxum2rg/a5klLTtZeW\n9LWkZcF/fTZiFxMaVCrJmDtb8ueOdZm0dBPn95/M+7N+whcWdQVBZlck1TI7UBGVM9nlTSIzBmdk\nFdDWzBoCTwJDjtje3syaHDFi4EFgopnVJrLY1oOZZXQuNyXEx3FHu5p8cW8b6lU8iT9/tIDrhs1k\nzZbdYUdzLqoy6yP5p6SPJPWUVF9SeUlVJZ0n6UlgGnBGRgebWQqwNZPt083st+DLb4HMilKaLsCI\n4P0IoGsWjnEuV1UvW4x3bzuHpy9vwIK127loYApDU1Zy0NeLd/lUZre2rgQeBuoQ6WifQuSp9l7A\nUuA8M/s6h3LcCqR/VsWAryTNltQ7XfspZpbWm/kr4ItHuJgUFyeua34aX/drS6ta5Xh63BKueHU6\nS9Z796LLf7K0sNUJnzxye2ysmTXIZJ/2wCtAKzPbErRVMrN1ksoDXwN3m1mKpG1mVirdsb+Z2VH7\nSYIC1BugatWqTdesWZNTP5Zzx8XM+HzBeh4ds4jtew5wR7ua3HVeLQolxIcdzblM5ejCVtEiqREw\nDOiSVkQAzGxd8N+NwGigWbBpg6SKwbEViUwgeVRmNsTMks0suVy5ctH6EZw7Jkl0bnQqE/q15bIm\np/Liv5fTadAUUldneOfXuTwltEIiqSrwMXCDmf2Yrr2YpBJp74ks75s28utT4MbgfdoEks7lCScX\nS6L/VU0YcUsz9h44zJWvzeDRMQvZte9g2NGcy5ao3dqS9C6RRbDKAhuAR4FEADMbLGkY0A1Iu+d0\n0MySJdUgchUCkbVQ3jGzp4NzlgE+AKoGx10VTGefKZ9ry8Wa3fsO8s8vlzJixmoqnlSYp69oSPs6\n5cOO5dz/yMlJGz8GXgfGm1meHHbihcTFqtlrfuPPH81n+cZdXH5mJR7uXI/SxZLCjuUckLN9JK8A\n1wLLJP1dUp1sp3POAdD0tJP5/J5W3HN+bT6b9wsX9p/MmLnr/EFGl6dkZa6tCWZ2HXAWsBqYEDyJ\nfrOkxGgHdC6/K5QQT78LT2fsPa2oXLoofd+bS68RqazfvifsaM5lSZY624O+iZuIPEPyPTCISGHJ\nqedInCvw6lY4iY/vOJf/u+QMpq/YwoX9U3jr2zUc9vXiXYw7ZiGRNJrIw4hFgUvN7DIze9/M7gaK\nRzugcwVJfJzo1boGX97bhiZVSvHwJwu5esi3rNy0K+xozmUoK53tncxs3BFthcxsX1ST5SDvbHd5\nkZkxavZanhq7mL0HD9P3/Nr0blODxPhQH/9yBUhOdrY/dZS2GccfyTl3PCRxVXIVJtzflvPrluef\nXy6ly0vTWLhue9jRnPsfmc3+W0FSU6CIpDMlnRW82hG5zeWcywXlSxTm1eubMvj6pmzatY8uL0/j\nmfFL2HvgUNjRnAMiD/xl5CIiHeyVgf7p2ncCD0Uxk3PuKDo2qECLGmX427glvDZ5JV8u/JVnrmhE\ni5plwo7mCris9JF0M7OPcilPVHgfictvpi/fzIMfL+Cnrb9zTbOq/KVTXU4q7KPxXc7K9pPtkq43\ns39Jup/ItO7/w8z6H+WwmOSFxOVHe/YfYsCEHxk2ZSXlShTiqa4NubCer6zgck5OdLYXC/5bHChx\nlJdzLkRFkuJ5qNMZjO7TkpOLJnHbyFTufGcOm3bmmQGVLp+I6nokscKvSFx+d+DQYV6bvIIXJi6n\nSFI8j3SuxxVnVUJS2NFcHpZjw38lPSvpJEmJkiZK2iTp+pyJ6ZzLCYnxcdx1Xm3G9W1FrfLFuX/U\nPHoO/46ft/4edjRXAGTlOZIOZrYD6Exkrq1awB+jGco5d2JqlS/BqD+04PHL6jNnzW9cNDCF4VNX\nccinWXFRlJVCkjZE+BJglJn501DOxbC4OHHjudX4ql9bmlUvzRNjF9N98HSWbdgZdjSXT2WlkIyV\n9APQFJgoqRywN7qxnHPZValUEd646WwG9GjM6s276fTCFAZO+JH9B/PkskIuhmWps11SaWC7mR2S\nVBQ4ycx+jXq6HOKd7a6g27JrH49/tphP5/1CnVNK8PduDTmz6slhx3IxLifn2gKoC/SQ1BPoTmQd\ndedcHlGmeCFeuOZMXr8xmR17D3DFq9N54rPF/L7f14t32ZfZFCkASHoLqAnMBdIm9zFgZBRzOeei\n4PwzTqFZ9dL844sfGD5tFV8t/pW/X9GIVrXLhh3N5WFZuSJJBlqaWR8zuzt43XOsgyQNl7RR0sIM\ntl8nab6kBcGKi42D9iqSvpG0WNIiSX3THfOYpHWS5gavTln9QZ1zESUKJ/JU14a83/sckuLjuP71\nmTwwah7bft8fdjSXR2WlkCwEKpzAud8EOmayfRXQ1swaAk8CQ4L2g8D9ZlYPOAe4U1K9dMcNMLMm\nwWsczrkT0rxGGcb1bU2fdjUZ/f06LuifwrgF6329eHfcslJIygKLJX0p6dO017EOMrMUYGsm26eb\n2W/Bl98SmWUYM1tvZnOC9zuBJUClLOR0zh2nwonx/KljXT69qyUVShaiz9tz+MNbs9mwwwdmuqw7\nZh8J8Fi0QwC3AuOPbJRUDTgTmJmu+a6g0z+VyJXLb0ce55w7PvVPLcknfVoybOoqBnz9Ixf0n8yf\nOtbl2mZViY/zaVZc5o55RWJmk4k80Z4YvJ8FzMmpAJLaEykkfz6ivTjwEXBv8GQ9wKtEOv6bAOuB\n5zM5b29JqZJSN23alFNxncu3EuLjuL1tTb64tw2NKpfk4U8WcsWr031FRndMWZlr6zbgQ+C1oKkS\n8ElOfHNJjYBhQBcz25KuPZFIEXnbzD5OazezDWZ2yMwOA0OBZhmd28yGmFmymSWXK1cuJ+I6VyBU\nL1uMf93anEFXN2Hdb79z2UtTefyzRezceyDsaC5GZaWP5E6gJbADwMyWAeWz+40lVQU+Bm4wsx/T\ntQt4HVhy5Jonkiqm+/JyIgMBnHM5TBJdmlRiYr92XNu8Km9OX80F/Sd7Z7w7qqwUkn1m9p9xgZIS\nOMpCV0eS9C4wA6gjaa2kWyXdLun2YJdHgDLAK8FQ3rRHz1sCNwDnHWWY77PBcOH5QHvgviz9lM65\nE1KyaGSo8Og+LSlTLNIZf9Mbs/hpi88q7P4rK0vtPgtsA3oCdwN9gMVm9tfox8sZPkWKc9l38NBh\nRs5Yw/NfLeXgYePu82pxW5saFEqIDzuai5JsL7Wb7kRxRDrDOwACvgSGWR66vvVC4lzO+XX7Xp4Y\nu4hxC36lZrliPNW1IS1qlgk7louCHCskwcnKAZhZnhz+5IXEuZz3zQ8beeTThfy8dQ9XnFWJhzqd\nQdnihcKO5XJQtidtVMRjkjYDS4GlweqIj+RkUOdc3tS+bnm+urctd7WvxWfzfuH85yfzzsyfOOyL\naBU4mXW230ek4/tsMyttZqWB5kBLSd7J7ZyjSFI8D1xUh/F9W1O3QgkeGr2A7oOns/iXHcc+2OUb\nmRWSG4BrzGxVWoOZrQSuJ9Lx7pxzQGSJ3/d6n8PzVzZmzZbfufSlqTz9+WJ27/Np6guCzApJoplt\nPrIx6CdJjF4k51xeJIluTSsz8f62XJVchaFTVnFB/8l8sfBXf/Ykn8uskGQ2p7TPN+2cO6pSRZN4\n5oqGfHRHC0oWSeT2f82m14hUft7qz57kVxmO2pJ0CNh9tE1AYTPLM1clPmrLuXAcOHSYN6etZsCE\nHzlsRt/zT6dX6+okxmd1cVYXpmyP2jKzeDM76SivEnmpiDjnwpMYH8dtbWrwdb+2tKldjn988QOX\nvDCF71bWhyYsAAAR30lEQVRluMKEy4P8nwXOuairVKoIQ3omM6xnMrv3HeKq12bwx1Hz2Lrb75Ln\nB15InHO55oJ6p/B1vzbc3jayKuN5z0/ig1k/+7MneZwXEudcriqalMCDF9fl83taU7t8cf700Xx6\nDJnB0l93hh3NnSAvJM65UNSpUIL3e7fg2e6NWL5xF5e8MIVnxi/h9/3+7Ele44XEOReauDhxVXIV\nJt7fjivOqsRrk1dyYf8UJizeEHY0dxy8kDjnQle6WBLPdm/MB39oQbFC8fQamUrvkams27Yn7Ggu\nC7yQOOdiRrPqpRl7d2v+3LEuKcs2cWH/yQxNWcmBQ4fDjuYy4YXEORdTkhLiuKNdTb6+ry0tapTh\n6XFLuPTFqcxe48+exCovJM65mFSldFGG3ZjMazc0ZfueA3R7dQYPfjTfnz2JQV5InHMxSxIX1a/A\nhH5t6d2mBqNmr6X9c5MYOWM1B/12V8yIaiGRNFzSRkkLM9h+naT5khZImi6pcbptHSUtlbRc0oPp\n2qtLmhm0vy8pKZo/g3MufMUKJfBQpzMY37c19U89iUfGLKLzi1P5duWWsKM5on9F8ibQMZPtq4C2\nZtYQeBIYAiApHngZuBioB1wjqV5wzD+AAWZWC/iNyHryzrkC4PRTSvB2r+YMvv4sdu49yNVDvuWu\nd+bwi4/uClVUC4mZpQAZ9pCZ2XQz+y348lugcvC+GbDczFaa2X7gPaCLJAHnAR8G+40AukYlvHMu\nJkmiY4OKTOjXlnsvqM3Xizdw/vOTeenfy9h74FDY8QqkWOojuRUYH7yvBPycbtvaoK0MsM3MDh7R\n7pwrYIokxXPvBacz8f62tKtTjue++pELB0zmq0W+kFZui4lCIqk9kULy5xw8Z29JqZJSN23alFOn\ndc7FmMonF+XV65vydq/mFE6Ip/dbs+k5/DuWb9wVdrQCI/RCIqkRMAzoYmZpPWfrgCrpdqsctG0B\nSklKOKL9/2NmQ8ws2cySy5UrF53wzrmY0bJWWcb1bc0jnesx9+dtdByYwtOfL2bn3gNhR8v3Qi0k\nkqoCHwM3mNmP6TbNAmoHI7SSgKuBTy1yvfoN0D3Y70ZgTG5mds7FrsT4OG5pVZ1vHmhH96aVGTZ1\nFe2fm8yoVJ+qPpoyXGo3R04uvQu0A8oCG4BHgUQAMxssaRjQDVgTHHIwbVlHSZ2AgUA8MNzMng7a\naxDpfC8NfA9cb2b7MsvhS+06VzDN+3kbj322iO9/2kaTKqV4/LL6NK5SKuxYeUZWl9qNaiGJFV5I\nnCu4Dh82Rn+/jmfG/8DmXfu4Krkyf7yoLuVKFAo7WszL9prtzjmXH8TFiW5NK/PNA5Gn4z+es47z\nnpvE61NX+WSQOcQLiXOuQChROJGHOp3BF/e24czTTubJsYu5eNAUpi7bHHa0PM8LiXOuQKlVvjgj\nbj6boT2T2X/wMNe/PpPb35rNz1t/DztanpVw7F2ccy5/kcSF9U6hde2yvD51FS/9eznfLN3IH9rW\n5I62NSmSFB92xDzFr0iccwVW4cR47mxfi4n3t6VD/Qq8MHEZF/SfzLgF6/3p+OPghcQ5V+CdWqoI\nL15zJu/3PocShRPo8/Ycrhs2k6W/7gw7Wp7ghcQ55wLNa5Rh7N2teLJLfRb9soNOL0zhsU8Xsf13\nfzo+M15InHMunYT4OG5oUY1JD7TjmmZVGDljNe2fn8R73/3EIX86/qi8kDjn3FGcXCyJp7o25LO7\nW1GrXHEe/HgBXV+exuw1vx374ALGC4lzzmWi/qklef8P5zDo6iZs2rmPbq9Op9/7c9m4Y2/Y0WKG\nFxLnnDsGSXRpUomJ97elT7uajJ2/nvbPTeK1ySvYf9CfjvdC4pxzWVSsUAJ/6liXr+5rQ4uaZXhm\n/A90HJjCN0s3hh0tVF5InHPuOFUrW4xhN57NGzefDcDNb8yi14hZrN68O+Rk4fBC4pxzJ6h9nfJ8\ncW8bHry4LjNWbKHDgBSe/eIHdu87eOyD8xEvJM45lw1JCXHc3rYm3zzQjs6NKvLKpBWc//xkxsxd\nV2CejvdC4pxzOaD8SYXp36MJH93RgrIlkuj73lyuem0Gi37ZHna0qPNC4pxzOajpaaUZc2crnrmi\nISs27ebSF6fyf58s4Lfd+8OOFjVeSJxzLofFx4lrmlXlm/vb0bNFNd797mfaPTeJt2aszpdPx3sh\ncc65KClZNJHHLqvP5/e0ol7Fk3h4zCI6vziVmSu3hB0tR0WtkEgaLmmjpIUZbK8raYakfZIeSNde\nR9LcdK8dku4Ntj0maV26bZ2ild8553JK3Qon8c5tzXnlurPYsecAPYZ8y93vfs/67XvCjpYjFK1R\nBZLaALuAkWbW4CjbywOnAV2B38zsuaPsEw+sA5qb2RpJjwG7jrZvZpKTky01NfUEfgrnnMtZe/Yf\nYvDkFQyevII4ibvOq8WtrapTODH2FtOSNNvMko+1X9SuSMwsBdiayfaNZjYLyGx+5vOBFWa2Jqfz\nOedcGIokxXPfhaczoV9b2p5ejn9+uZQOA1L4evGGPDtcONb7SK4G3j2i7S5J84NbZyeHEco557Kr\nSumiDL6hKf+6tTlJCXHcNjKVm96YxYpNu8KOdtxitpBISgIuA0ala34VqAk0AdYDz2dyfG9JqZJS\nN23aFNWszjl3olrVLsv4vq15uHM95qz5jYsGpPC3cUvYuTfvLKYVs4UEuBiYY2Yb0hrMbIOZHTKz\nw8BQoFlGB5vZEDNLNrPkcuXK5UJc55w7MYnxcdzaqjrf/LEd3c6qzNApKznv+cl8NHsth/PAcOFY\nLiTXcMRtLUkV0315OXDUEWHOOZcXlS1eiH90b8QnfVpSqVQR7h81j+6DpzN/7bawo2UqmqO23gXa\nAWWBDcCjQCKAmQ2WVAFIBU4CDhMZ4VXPzHZIKgb8BNQws+3pzvkWkdtaBqwG/mBm64+VxUdtOefy\nmsOHjY+/X8ffx//Alt376JFchQcuqkPZ4oVyLUNWR21FrZDEEi8kzrm8asfeA7w4cRlvTFtNkaR4\n+l14OtefcxqJ8dG/oRT68F/nnHPZd1LhRP56ST2+uLcNTaqU4vHPFnPJC1OYvnxz2NH+wwuJc87l\nAbXKF2fkLc0YckNT9hw4xLXDZnLHv2bz89bfw47mhcQ55/IKSXSoX4Gv72vLAx1OZ9LSTVzQfzID\nJ/zI3gOHQsvlhcQ55/KYwonx3HVebSbe35YL653CwAnLOP/5yYxfsD6Up+O9kDjnXB51aqkivHTt\nWbzX+xxKFE7gjrfncP3rM/lxw85czeGFxDnn8rhzapRh7N2teKJLfRau28HFg6bw+GeL2L4nd56O\n90LinHP5QEJ8HD1bVOObB9px9dlVeHP6as57bhLTV0R/dJcXEuecy0dKF0vi6csb8tldrahfqSQ1\nyhaP+vdMiPp3cM45l+saVCrJyFsynI4wR/kViXPOuWzxQuKccy5bvJA455zLFi8kzjnnssULiXPO\nuWzxQuKccy5bvJA455zLFi8kzjnnsqVArJAoaROw5gQPLwvEzgoy/+W5jo/nOj6e6/jEai7IXrbT\nzKzcsXYqEIUkOySlZmWpydzmuY6P5zo+nuv4xGouyJ1sfmvLOedctnghcc45ly1eSI5tSNgBMuC5\njo/nOj6e6/jEai7IhWzeR+Kccy5b/IrEOedctnghyYCkKyUtknRYUvIR2/4iabmkpZIuCjHjY5LW\nSZobvDqFlSXI0zH4TJZLejDMLOlJWi1pQfAZpYaYY7ikjZIWpmsrLelrScuC/54cI7lC/92SVEXS\nN5IWB38W+wbtoX5mmeQK9TOTVFjSd5LmBbkeD9qrS5oZ/Ll8X1JSjn9zM/PXUV7AGUAdYBKQnK69\nHjAPKARUB1YA8SFlfAx4IOzPKsgSH3wWNYCk4DOqF3auINtqoGwM5GgDnAUsTNf2LPBg8P5B4B8x\nkiv03y2gInBW8L4E8GPw5y/UzyyTXKF+ZoCA4sH7RGAmcA7wAXB10D4YuCOnv7dfkWTAzJaY2dKj\nbOoCvGdm+8xsFbAcyJ1lyGJbM2C5ma00s/3Ae0Q+KxcwsxRg6xHNXYARwfsRQNdcDUWGuUJnZuvN\nbE7wfiewBKhEyJ9ZJrlCZRG7gi8Tg5cB5wEfBu1R+by8kBy/SsDP6b5eS7i/RHdJmh/cnsj12yLp\nxNrnkp4BX0maLal32GGOcIqZrQ/e/wqcEmaYI8TK7xaSqgFnEvlXdsx8ZkfkgpA/M0nxkuYCG4Gv\nidwl2GZmB4NdovLnskAXEkkTJC08yitm/iV9jIyvAjWBJsB64PlQw8auVmZ2FnAxcKekNmEHOhqL\n3HuIlWGUMfO7Jak48BFwr5ntSL8tzM/sKLlC/8zM7JCZNQEqE7lLUDc3vm9CbnyTWGVmF5zAYeuA\nKum+rhy0RUVWM0oaCoyNVo4syNXP5XiY2brgvxsljSbyBywl3FT/sUFSRTNbL6kikX9Jhs7MNqS9\nD/N3S1Iikb+s3zazj4Pm0D+zo+WKlc8syLJN0jdAC6CUpITgqiQqfy4L9BXJCfoUuFpSIUnVgdrA\nd2EECf4QpbkcWJjRvrlgFlA7GCGSBFxN5LMKlaRikkqkvQc6EO7ndKRPgRuD9zcCY0LM8h+x8Lsl\nScDrwBIz659uU6ifWUa5wv7MJJWTVCp4XwS4kEj/zTdA92C36HxeYY0wiPUXkV+EtcA+YAPwZbpt\nfyVy73EpcHGIGd8CFgDzifzhqhjyZ9aJyAiWFcBfw/5/GGSqQWQE2TxgUZi5gHeJ3PI4EPxu3QqU\nASYCy4AJQOkYyRX67xbQishtq/nA3ODVKezPLJNcoX5mQCPg++D7LwQeCdprEPnH7nJgFFAop7+3\nP9nunHMuW/zWlnPOuWzxQuKccy5bvJA455zLFi8kzjnnssULiXPOuWzxQuLyLUll0s3E+usRM7Pm\n/AyoOUDSLZIqRPH8xSRNkhQXfF1X0vhgJt05kt6TVF5SE0mvRyuHy18K9JPtLn8zsy1EpqtA0mPA\nLjN7LtRQkSzxZnYog823AHOIzCGV1fOlPbWcFb2AUWZ2WFJR4HPgbjMbF5zrfKCMmc2VVENSJQtm\nBnAuI35F4gokSTcGazfMlfSKpDhJCZK2SeofrOfwpaTmkiZLWpm2voSkXpJGB+3LJP1fFs87UNJ8\noJmkxyXNCuZNG6yIHkQK3/tpV02S1qZ7WvkcSROC909JGilpGvBm8D36B997vqReGfzo1/HfJ5uv\nByanFREAM5toZkuCL8cCPXLuU3f5lRcSV+BIakBk5oJzLTLBXQKRKV0ASgLjzaw+sJ/IGhPnA1cC\nT6Q7TTMi03E3Aa4NbgUd67wpZtbIzGYAg8zsbKBhsK2jmb1P5CnpHmbWxCLT8WemLnC+mV0P9AY2\nmlkz4Gwik1NWPeLnLgxUNrO1QVMDYHYm508FWh8jg3N+a8sVSBcQ+cs2NTJtEkX47xT4e8zs6+D9\nAmC7mR2UtAColu4cX5rZbwCSPiEybUZCJufdD4xOd/z5kv4IFAbKEvkLffxx/hxjzGxv8L4DcIak\n9IWrNvBTuv3Lc3zrjmwETj3OTK4A8kLiCiIBw83s4f9plBKI/IWf5jCRudbS3qf/83Lk3EJ2jPPu\nsWA+oqBv4iUiq+ytk/QUkYJyNAf5752DI/fZfcTP1MfMJmZwHoA9R5xjEdA8k/0LB8c4lym/teUK\nognAVZLKwn9Gd1U9xjFH6iCpVFAUugDTjuO8RYgUps3BzMTd0m3bSWT51jSrgabB+/T7HelLoE9Q\ntJBUJ5gB9j/MbBNQJN2ItbeAdpI6pu0jqb2kM4IvTye2Zkp2McqvSFyBY2YLJD0OTAiGwR4Abgd+\nOY7TzCLSaX0qMMLM5gJk5bxmtkXSCGAxkVl3Z6bb/AYwTNIeIv0wjwFDJW0j8zVUXgOqAnOD22ob\nOfpSxxOAc4FJZva7pM7AAEkvBnnnAn2DfdsTI9Pau9jms/86d5yCEVENzOzesLMcL0lnE7kFdvMx\n9itCZB2LlpkMVXYO8FtbzhUoZjYLmJr2QGImqgJ/8iLissKvSJxzzmWLX5E455zLFi8kzjnnssUL\niXPOuWzxQuKccy5bvJA455zLFi8kzjnnsuX/AQ2+ym0V1ftSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1115d7470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VPXZ//H3nT1hC0tQCJtsCrITNhVwadWqFRdAVBRR\nUJT2cXl+7VPbp+tja7Vq1VpZBNzqAu5LFXdBZQ2Cyk4AlZ2whAAJEOD+/TEHTRGyQJKTTD6v65or\nM985M3PPdSCfnO855z7m7oiIiJRUTNgFiIhI1aLgEBGRUlFwiIhIqSg4RESkVBQcIiJSKgoOEREp\nFQWHiIiUioJDRERKRcEhIiKlEhd2AeWhQYMG3qJFi7DLEBGpUubNm7fF3dOKWy4qg6NFixZkZmaG\nXYaISJViZt+UZDlNVYmISKkoOEREpFQUHCIiUioKDhERKZVyCw4zm2Rmm81sYaGxemb2npmtCH7W\nDcbNzB42sywz+9LMuhV6zbBg+RVmNqy86hURkZIpzy2OJ4DzDxv7FfCBu7cBPggeA/wEaBPcbgTG\nQCRogN8DvYCewO8PhY2IiISj3ILD3acD2w4bHgA8Gdx/Erik0PhTHjELSDWzRsB5wHvuvs3dtwPv\n8cMwEhGRClTR53Gc4O4bgvsbgROC++nAmkLLrQ3GjjZe7tydvfsPsmvvfnbv3R/8PMCuvQXs2nuA\n3cH47r0HqFsjniZ1k2lSN4X01GRqJEbl6TEiIkCIJwC6u5tZmV3w3MxuJDLNRbNmzY7pPRau28FN\nT8/7Liz2Hzy28uqmxH8XIpFASSa9bsp392slxR/T+4qIVAYVHRybzKyRu28IpqI2B+PrgKaFlmsS\njK0Dzjxs/OMjvbG7jwfGA2RkZBzTb/w6yfH0blmfmomx1EiMo2ZSHDUT46iREEeNxDhqJUV+Hnq+\nRmIcKfGxbMvbx9rt+azdns+67fms3Z7H2u35ZGXv4uPlm9lTcPAHn/N9qKSQXvf7gGmSmkLt5DjM\n7Fi+gohIuavo4HgdGAb8Nfj5WqHxn5nZ80R2hO8IwuUd4C+FdoifC9xZXsU1rZfC/YM7l/p1DWsl\n0bBWEt2a/XC/vbuzdfe+H4TKupx8vt66m0+ztpC378B/vKZWYlyhMPl+SyU9NXI/NSVewSIioSm3\n4DCz54hsLTQws7VEjo76KzDFzG4AvgEGB4u/BVwAZAF5wHAAd99mZv8HzA2W+5O7H77DvVIzMxrU\nTKRBzUS6NE39wfPuTk5eQbDFkse6nPzv7q/dns+sVdvYtXf/f7ymRkJsECwpQaAk0yqtJr1b1aem\n9q+ISDkz9zLbzVBpZGRkeLQ0OXR3cvP3szYn77vpsLXb84Ktl8j93D2RYImLMbo3r0v/k9Po1yaN\n9o1qExOjLRMRKRkzm+fuGcUup+Co+nL3FLBoXS7TV2QzbVk2izfkAtCgZiL92jSg/8lpnNG6AfVr\nJoZcqYhUZgqOahQch9u8cw+fLN/CtOXZfLIim+15BZhBx/Q69G+bRv+2aXRpmkpcrDrOiMj3FBzV\nODgKO3DQWbhuB9OXZzNteTbz1+Rw4KBTKymOM1o3oF8QJI1Tk8MuVURCpuBQcBzRjvwCZmRFtkam\nL89m/Y49APRoUZdLuqZzUcfG1EnReSYi1ZGCQ8FRLHcna/Mu3l28iVfnr2PF5l0kxMZw1ilpXNq1\nCWedkkZiXGzYZYpIBVFwKDhKxd1ZtD6XV+av47UF69myay91kuO5sFMjLuuaTvfmdXXuiEiUU3Ao\nOI7Z/gMH+WzlVl75fC3vLNpEfsEBmtZL5tIu6VzSNZ2WaTXDLlFEyoGCQ8FRJnbv3c87izbyyvx1\nfJa1hYMOnZumclnXdC7q1EiH+IpEEQWHgqPMbcrdw+sL1vPy/HUs2ZBLXIzRv20a153egjNaN9BU\nlkgVp+BQcJSrpRsj+0Ne+Xwdm3fupXPTVH5+VmvOaddQASJSRSk4FBwVYu/+A7w0bx1jpmWxZls+\n7RrVZvRZrfhJh0bEqt2JSJWi4FBwVKj9Bw7y+hfreeSjLFZl76ZlWg1Gn9mai7s0Jl5nqItUCQoO\nBUcoDhx0pi7cyD8+XMHSjTtpWi+ZUf1bMbB7E50TIlLJKTgUHKFydz5Yspl/fJTFF2tyOLF2Ejf2\na8mVPZuRnKAAEamMFBwKjkrB3fksayv/+HAFs1dvo36NBG7oexLX9G6uS+iKVDIKDgVHpTP36208\n8mEW05ZnUzspjuGnn8Tw01uQmpIQdmkigoJDwVGJfbk2h0c+zOLdxZuokRDL0D7NGXFGS9Jq6WRC\nkTApOBQcld7Sjbk8+tFK3vxyPfGxMVzZsxk39W9Jozpq8S4SBgWHgqPKWJW9izEfr+SV+eswg4Hd\nm3Jz/1Y0q58Sdmki1YqCQ8FR5azdnse4aauYnLmGAwedAZ0bc8tZrWjdsFbYpYlUCwoOBUeVtSl3\nD49NX8Uzs79lz/4DXNChEbec1YpTG9cJuzSRqKbgUHBUeVt37WXSZ6t5asY37Ny7n3NOacjos1vT\nrVndsEsTiUoKDgVH1NiRX8BTM75m4meryckr4PTW9fnZWW3o3bKeGiqKlCEFh4Ij6uzeu59nZn/D\n+Omr2bJrLxnN6/Kzs1vTv22aAkSkDCg4FBxRa0/BAaZkrmHsxytZv2MPHdPrMPqs1pzb/gRi1JFX\n5JgpOBQcUW/f/oO8Mn8tj368km+25tH2hJqMPqs1F3VqrJbuIsegpMGhftdSZSXExXBFj2Z8cEd/\nHhrSBXe49fkFnPfgdOas3hZ2eSJRS8EhVV5cbAwDuqTzzm39GHN1N/YUHGDwuJnc+fJX7MgvCLs8\nkaij4JCoERNj/KRjI969vR8j+57E5Lnf8uMHpvH2VxuIxilZkbAoOCTqpCTE8ZsL2/P6z84grVYi\nNz/zOSOfmseGHflhlyYSFRQcErU6pNfhtdGn85sL2vFpVjY/fmA6T874mgMHtfUhcjwUHBLV4mJj\nGNmvJe/d3p+uzVL5/euLGDh2Bss27gy7NJEqS8Eh1ULTeik8dX1P/n5FZ77ZmseFD3/Cfe8sY0/B\ngbBLE6lyFBxSbZgZl3Ztwvt39OfiLo155KMsLnjoE2at2hp2aSJVioJDqp16NRJ4YHAXnr6hJ/sP\nOkPGz+J/XvySHXk6dFekJBQcUm31bZPGO7f146b+LXnx87Wc88A03vhivQ7dFSmGgkOqteSEWO78\nSTteG306jeok8fPn5nPDk5msy9GhuyJHo+AQIXLo7iu3nMb/XtiOmSu38uMHpjHp09U6dFfkCBQc\nIoG42BhG9G3Ju7f3o0eLevzpzcVcNmYGSzbkhl2aSKUSSnCY2a1mttDMFpnZbcFYZzObaWZfmdkb\nZla70PJ3mlmWmS0zs/PCqFmqj6b1UnhieA8eGtKFtdvy+Ok/PuXeqUt16K5IoMKDw8w6ACOBnkBn\n4CIzaw1MAH7l7h2BV4BfBMu3B4YApwLnA4+aWWxF1y3Vi5kxoEs6H/x3fy7tms6jH6/k/AenMyNr\nS9iliYQujC2OdsBsd89z9/3ANOAyoC0wPVjmPeDy4P4A4Hl33+vuq4EsIqEjUu5SUxL426DOPDui\nFwBXTZjNL174gu2794VcmUh4wgiOhUBfM6tvZinABUBTYBGRkAAYFIwBpANrCr1+bTAmUmFOa92A\nqbf145YzW/HK/HX86IFpvLZgnQ7dlWqpwoPD3ZcA9wDvAlOBBcAB4HrgFjObB9QCSvUnnZndaGaZ\nZpaZnZ1dxlWLQFJ8LL88/xTe+PkZNKmXwq3PL2D4E3NZsy0v7NJEKlQoO8fdfaK7d3f3fsB2YLm7\nL3X3c929O/AcsDJYfB3fb30ANAnGDn/P8e6e4e4ZaWlp5f0VpBpr16g2L998Gr//aXvmrN7GuX+f\nzoRPVrH/wMGwSxOpEGEdVdUw+NmMyP6NZwuNxQD/C4wNFn8dGGJmiWZ2EtAGmFPxVYt8LzbGGH76\nSbx3R39Oa1Wfu/69hEsfncGi9TvCLk2k3IV1HsdLZrYYeAMY7e45wJVmthxYCqwHHgdw90XAFGAx\nkamt0e6u4yKlUkhPTWbCsAweuaorG3bs4eJHPuPut5eQv0//RCV6WTTu3MvIyPDMzMywy5BqZkde\nAXe/vYTn566hWb0U/nxpB/q20bSpVB1mNs/dM4pbTmeOi5SROinx/PXyTjx/Y2/iYoxrJs7hjskL\n2KZDdyXKKDhEyljvlvV569a+/Pzs1rz+xXp+9MA0Xpm/VofuStRQcIiUg6T4WP773JP593/1pXn9\nFG6f/AXXTprDt1t16K5UfQoOkXJ08om1eHHUafxpwKnM/zaHcx+cxsufrw27LJHjUmRwmFkfM/un\nmX1pZtlm9q2ZvWVmo82sTkUVKVKVxcYY1/ZpwXt39KNr07rcMeULJn26OuyyRI7ZUYPDzN4GRgDv\nEGku2AhoT+QciyTgNTO7uCKKFIkGjeok88T1PTj/1BP505uLeeDdZdrvIVXSUQ/HNbMG7l5kK9CS\nLBMGHY4rldmBg86vX/6KyZlruKZ3c/548anExFjYZYmU+HDcuKM9caRAMLN67r6tqGVEpGixMcZf\nL+9Iao14xk1bxY78Au4f3Jn4WO1ylKqhqKmq081sSXCxpV5m9h4w18zWmFmfCqxRJOqYGXf+pB3/\nc/4pvP7FekY+lamzzaXKKOpPnL8Dg4ns5/g38Ed3b0Wk9fl9FVCbSNS7+cxW3H1ZR6Yvz+aaibPZ\nkV8QdkkixSoqOOLd/St3nwlku/unAO7+OZBcIdWJVANX9mzGI1d144u1OQwZP4vNO/eEXZJIkYoK\njsLP3XnYcwnlUItItXVBx0ZMuq4H32zdzaCxM3WND6nUigqO3wZX6MPdXz00aGatgKfKuzCR6qZv\nmzT+NaIXOXkFXD5mBss27gy7JJEjOmpwuPvr7v6DP3vcfaW731u+ZYlUT92a1eWFUX0wg8HjZvL5\nt9vDLknkB4o9/s/MOlZEISIS0faESJuS1JR4hk6YzScrdClkqVyKazlyDvBoBdUiIoGm9VJ4YVQf\nmtevwfVPzOWtrzaEXZLId4o6j+Nq4G9ELu0qIhWsYa0knr+xN52bpPKzZz/nuTnfhl2SCFD0Fsfj\nwEB313aySEjqJMfz9A296N82jTtf/ooxH68MuySRIoPj98BEM9M5GyIhSk6IZfy1GQzo0ph7pi7l\n7reWqDmihKqoXlV3m9la4FXgvIorSUQOFx8bw98Hd6FOcjzjpq8iJ6+AP1/agTj1t5IQHDU4ANz9\naTNbX1HFiMjRxcQYf7z4VFJTEnj4gxXsyC/goSu7kBgXG3ZpUs0U++eKu39QEYWISPHMjDt+3Jbf\nXdSeqYs2cv0Tc9m1d3/YZUk1U+QWB4CZxQIXAi0KL+/uD5RfWSJSlOvPOIk6yfH88qUvuXrCbJ64\nrgd1a6gTkFSMkkyQvgFcB9QHahW6iUiILu/ehLFDu7NkQy6Dxs1kw478sEuSauKoVwD8bgGzL929\nUwXVUyZ0BUCpTmat2sqIJzOpkxzPv0b04qQGNcIuSaqokl4BsCRbHG+b2bllUJOIlIPeLevz/I29\n2VNwgEFjZ7Bw3Y6wS5IoV5LgmAW8Ymb5ZpZrZjvNLLe8CxORkuuQXocpo/qQGBfLleNnMWf1tuJf\nJHKMShIcDwB9gBR3r+3utdy9djnXJSKl1CqtJi+M6kPD2olcM3E2Hy7dFHZJEqVKEhxrgIWuU1VF\nKr3Gqcm8MOo0Tj6xFiOfmser89eFXZJEoWIPxwVWAR+b2dvA3kODOhxXpHKqVyOBZ0f2ZuSTmdw2\neQE5efu47vSTwi5LokhJtjhWAx8QuVysDscVqQJqJsbx+PAenNv+BP7wxmIefH+5+ltJmTnqFoeZ\n3QlMdfc/VmA9IlJGkuJjefTqbvzq5a948P0V5OQV8LuL2hMTY2GXJlVcUVNVq4Bbzawz8AXwNvCu\nu+taliJVRFxsDPde3onU5HgmfLqaHfkF3DuwE/FqjijHoajuuJOByQBm1hU4H3g5aEHyPpGtkTkV\nUqWIHLOYGOM3F7ajbo0E/vbOMnLzC/jn1d1IildzRDk2Jfqzw93nu/vd7n4WcBGwCBhRrpWJSJkx\nM0af1Zo/X9qBD5dt5tqJc8jdUxB2WVJFlaTJ4ZEuHZsD/G/ZlyMi5enqXs2pnRTPHVMWMGTcLJ68\nvidptRLDLkuqmJJscdwATACuDm6PAf8DfGZm15RjbSJSDn7auTEThvVg9ZbdDB43kzXb8sIuSaqY\nkgRHHNDO3S9398uB9oADvYgEiIhUMf3bpvGvET3Zumsvg8bOZMWmnWGXJFVISYKjqbsX7l2wORjb\nBhzTJKmZ3WpmC81skZndFox1MbNZZrbAzDLNrGcwbmb2sJllmdmXZtbtWD5TRP5T9+b1mDKqDwfc\nGTRuJgvW5IRdklQRJQmOj83sTTMbZmbDgNeCsRpE9nWUipl1AEYCPYHOwEVm1hq4F/iju3cBfhc8\nBvgJ0Ca43QiMKe1nisiRnXJibV4adRq1k+K56rFZfJa1JeySpAooSXCMBh4HugS3p9z9FnffHRxl\nVVrtgNnunufu+4FpwGVEpr8ONU+sAxy61vmA4DPd3WcBqWbW6Bg+V0SOoFn9FF4c1Ydm9VIY/vhc\npi7cEHZJUsmVJDjOd/eX3P324PaimY06js9cCPQ1s/pmlgJcADQFbgP+ZmZrgPuAO4Pl04k0Wjxk\nbTAmImWkYe0kJt/Yhw7ptbnlmc+ZMndN8S+SaqskwfFbMzv70AMz+yWRrYBj4u5LgHuAd4GpwALg\nAHAzcLu7NwVuByaW5n3N7MZg30hmdnb2sZYnUm3VSYlcQfCMNmn88qUvGT99ZdglSSVVkuC4GPiL\nmfU1sz8TOZrqmIMDwN0nunt3d+8HbAeWA8OAl4NFXiCyDwRgHZEtkkOaBGOHv+d4d89w94y0tLTj\nKU+k2kpJiGPCtRlc1KkRf3lrKfdMXarmiPIDxQaHu28hEh7/BBoDA9193/F8qJk1DH42I7J/41ki\n+zT6B4ucDawI7r8OXBscXdUb2OHumoQVKScJcTE8NKQrV/dqxpiPV/LrVxZy4KDCQ75XVHfcnUR2\nWB+SALQEBpqZH+dVAF8ys/pEDucd7e45ZjYSeMjM4oA9RI6gAniLyH6QLCAPGH4cnysiJRAbY9x1\nSQfqpiTwyEdZ5OYX8MAVnUmMU38rKbrJYbldc8Pd+x5h7FOg+xHGnciRXSJSgcyM/3feyaSmxHPX\nv5eQu6eAsUO7UyOxJNd/k2h21KkqM2tR1AuDqaMmZV2QiFQuI/q25G8DO/FZ1haGTpxNTt5xzVRL\nFChqH8ffzOwlM7vWzE41s4Zm1szMzjaz/wM+I3JOhohEuUEZTRkztDuL1udyxbhZbMrdE3ZJEqKj\nBoe7DwJ+C5xMZMf4J0TOGh8BLAPOdvf3KqJIEQnfeaeeyBPDe7B2ex6Xj5nB11t2h12ShMSi8VC7\njIwMz8zMDLsMkaj05dochk2aQ2xMDE9d35P2jY/nOBmpTMxsnrtnFLecrh8pIqXSqUkqL4w6jfhY\n44rxM8n8elvYJUkFU3CISKm1bliTF28+jbSaiQydOJuPlm0OuySpQAoOETkm6anJTBnVh9YNazLy\nyUxeW/CDhg4SpYoNDjN72cwuNDOFjIj8hwY1E3luZG+6N6/LbZMX8PTMr8MuSSpAScLgUeAqYIWZ\n/dXMTi7nmkSkCqmVFM+T1/fknFNO4LevLeIfH6xQf6soV5JeVe+7+9VAN+Br4H0zm2Fmw80svrwL\nFJHKLyk+lrFDu3FZt3Tuf285//fmEg6qv1XUKlHvgKCv1FDgGmA+8AxwBpGOtmeWV3EiUnXExcZw\n38DOpCYnMOmz1eTk7+PeyzsRF6tZ7mhTbHCY2StETgJ8Gvhpoc60k81MJ0uIyHdiYozfXtSOuinx\n3P/ecnLz9/PIVV1JildzxGhSkj8FHnP39u5+96HQMLNEgJKcKCIi1YuZ8fNz2vB/A07lg6WbGDZp\nDjv3FIRdlpShkgTHXUcYm1nWhYhIdLmmTwsevKIL877ZzpWPzWLrrr1hlyRlpKjuuCeaWXcg2cy6\nmlm34HYmkFJhFYpIlTWgSzqPXZtB1uZdDBo3k3U5+WGXJGWgqC2O84D7iFyq9QHg/uB2B/Dr8i9N\nRKLBWac05OkbepG9cy8Dx8wga/OusEuS41Rsk0Mzu9zdX6qgesqEmhyKVD6L1+dy7aQ5HHTnieE9\n6NQkNeyS5DDH3eTQzIYGd1uY2R2H38qsUhGpFto3rs2Lo/qQkhDLleNnMWPllrBLkmNU1FRVjeBn\nTaDWEW4iIqXSokENXrr5NNLrJnPd43N5Z9HGsEuSY6DrcYhIhcvJ28fwJ+byxZoc7rm8E4MymoZd\nklCG1+Mws3vNrLaZxZvZB2aWXWgaS0Sk1FJTEvjXDb04vXUDfvHil0z4ZFXYJUkplOQ8jnPdPRe4\niEivqtbAL8qzKBGJfjUS45gwLIMLOzbirn8v4b53lqk5YhVRkl5Vh5a5EHjB3XeYWTmWJCLVRWJc\nLA9f2ZXayXE88lEW2/P28acBHYiN0e+YyqwkwfGmmS0F8oGbzSwN2FO+ZYlIdREbY/zl0o6kpiQw\n5uOV7Mgv4IHBXUiIU3PEyqrY4HD3X5nZvcAOdz9gZruBAeVfmohUF2bG/5x/CqnJ8dz99lJ27tnP\nmKHdSEkoUQNvqWAlXSunEDmfo/DyT5VDPSJSjd3UvxWpKfHc+fJXXDNxDpOG9aBOii77U9mU5Kiq\np4m0HjkD6BHc1BVXRMrFFT2a8ejV3fhq7Q6uGD+TzbmaGa9sStJyZAnQ3qvQ4Q46j0Ok6vt0xRZu\nfDqTBjUT+dcNvWhWX71Vy1uZnccBLAROPP6SRERK7ow2DXh2ZG9y9xQwcOwMlm7MDbskCZQkOBoA\ni83sHTN7/dCtvAsTEenSNJUXbupDjBmDx85k3jfbwy5JKNlUVf8jjbv7tHKpqAxoqkokuqzZlse1\nk+awcccexl7Tnf5t08IuKSqV2VRVEBBfA/HB/bnA58ddoYhICTWtl8KUm/pwUoMajHhyLm98sT7s\nkqq1khxVNRJ4ERgXDKUDr5ZnUSIih0urlcjzN/Wma9O6/Nfz83lm9jdhl1RtlWQfx2jgdCAXwN1X\nAA3LsygRkSOpnRTPUzf05OyTG/KbVxbyz4+y1N8qBCUJjr3uvu/Qg+AkQK0pEQlFUnwsY6/pziVd\nGvO3d5bxl7eWKDwqWEnOHJ9mZr8Gks3sx8AtwBvlW5aIyNHFx8bwwOAupKYk8Ngnq8nJK+DuyzoS\nF6v+VhWhJMHxK+AG4CvgJuAtYEJ5FiUiUpyYGOP3P21Pako8D76/gh35BTx8ZVeS4mPDLi3qlaTJ\n4UEzexV41d2zK6AmEZESMTNu+1FbUpPj+cMbi7n+ibmMvzaDmolqjliejrpdZxF/MLMtwDJgWXD1\nv99VXHkiIsW77vSTePCKLsxevY2rHpvFtt37in+RHLOiJgRvJ3I0VQ93r+fu9YBewOlmdvvxfKiZ\n3WpmC81skZndFoxNNrMFwe1rM1tQaPk7zSzLzJaZ2XnH89kiEp0u6ZrO+Gu6s2zjTgaNncH6nPyw\nS4paRQXHNcCV7r760IC7rwKGAtce6weaWQdgJNAT6AxcZGat3f0Kd+/i7l2Al4CXg+XbA0OAU4Hz\ngUfNTJOYIvID57Q7gaeu78nm3L0MGjuTldm7wi4pKhUVHPHuvuXwwWA/x/E0yG8HzHb3PHffD0wD\nLjv0pEWuSzsYeC4YGgA87+57gxDLIhI6IiI/0KtlfZ67sTd79x9g8NiZLFy3I+ySok5RwVHUJOHx\nTCAuBPqaWX0zSwEuAJoWer4vsCk40RAiZ6qvKfT82mBMROSIOqTX4YVRp5EUH8uQ8bOYtWpr2CVF\nlaKCo7OZ5R7hthPoeKwf6O5LgHuAd4GpwALgQKFFruT7rY0SM7MbzSzTzDKzs3Xwl0h1d1KDGrx4\ncx9OrJPEsElzeH/xprBLihpHDQ53j3X32ke41XL347qWo7tPdPfu7t4P2A4sh+/OSr8MmFxo8XX8\n5xZJk2Ds8Pcc7+4Z7p6RlqbOmSICjeokM+WmPpxyYi1u+tc8Xv58bdglRYVQTrM0s4bBz2ZEguLZ\n4KkfAUvdvfDafR0YYmaJZnYS0AaYU5H1ikjVVa9GAs+M7E3vlvW4Y8oXTPp0dfEvkiKFdZbMS2ZW\nHygARrt7TjA+hMOmqdx9kZlNARYD+4PlDyAiUkI1E+OYdF0Pbn1uAX96czE5efu4/cdtiRyLI6VV\n7IWcqiJdyElEjmT/gYP85pWFTM5cw7V9mvOHn55KTIzC45CSXshJ5+WLSLURFxvDXy/vSGpKPOOm\nryInr4D7B3cmXs0RS0XBISLViplx5wXtSE1J4J6pS9m5p4BHr+5OcoLOKy4pxayIVEs3n9mKuy/r\nyMfLs7l20mx25BeEXVKVoeAQkWrryp7NeOTKbixYk8OQ8bPYvHNP2CVVCQoOEanWLuzUiInDevD1\nlt0MGjuTNdvywi6p0lNwiEi1169tGs+M7EVOXgEDx85g+aadYZdUqSk4RESAbs3qMuWmPrjD4HEz\nmf/t9rBLqrQUHCIigZNPrMVLN59GneR4rp4wm09WqO/dkSg4REQKaVovhRdG9aFZvRSuf2Iub321\nIeySKh0Fh4jIYRrWSmLyTX3o3CSVnz37Oc/N+TbskioVBYeIyBHUSY7n6Rt60a9tGne+/BVjPl4Z\ndkmVhoJDROQokhNiGX9NBhd3bsw9U5dy91tLiMb+fqWlliMiIkVIiIvhwSu6UCf5+/5Wf7msI7HV\nuDmigkNEpBgxMcafBpxK3ZR4Hv4wi9w9BTw4pAuJcdWzv5WmqkRESsDMuOPck/ntRe15e+FGbngi\nk91794ddVigUHCIipXDDGSdx/6DOzFy1lasmzGb77n1hl1ThFBwiIqV0efcmjB3anSUbchk8biYb\nd1Sv5ogKDhGRY/Dj9ifw5PCebNixh8vHzGD1lt1hl1RhFBwiIseoT6v6PDeyN/kFBxg0dgaL1u8I\nu6QKoeBVDA0tAAAMk0lEQVQQETkOHZvU4YVRfUiIjWHIuFnMWb0t7JLKnYJDROQ4tUqryYs3n0bD\n2olcM3E2Hy7dFHZJ5UrBISJSBhqnJjPlpj6cfGItbnxqHq/OXxd2SeVGwSEiUkbq10zk2ZG96dGi\nHrdNXsCTM74Ou6RyoeAQESlDNRPjeHx4D37c/gR+//oiHnp/RdT1t1JwiIiUsaT4WMZc3Y2B3Zvw\n9/eX88c3FnPwYPSEh3pViYiUg7jYGO69vBOpyfFM+HQ1O/ILuHdgJ+Jjq/7f6woOEZFyEhNj/ObC\ndtStkcDf3lnGzj0FPHJVN5Liq3ZzxKoffSIilZiZMfqs1tx1SQc+WLqZayfNIXdPQdhlHRcFh4hI\nBRjauzkPD+nK599s58rxs9iya2/YJR0zBYeISAX5aefGTBiWwars3QwaO5O12/PCLumYKDhERCrQ\nmSc35F8jerJ1114GjpnJik07wy6p1BQcIiIVrHvzeky+qQ8H3Bk8biYL1uSEXVKpKDhERELQrlFt\nXhzVh1pJ8Vz12Cw+y9oSdkklpuAQEQlJ8/o1eHFUH5rWTWH443OZunBj2CWViIJDRCREDWsnMeWm\nPnRIr80tz8xjytw1YZdULAWHiEjI6qTE868RvTijTRq/fOlLxk9fGXZJRVJwiIhUAikJcUy4NoOL\nOjXiL28t5d6pSyttc0S1HBERqSQS4mJ4aEhXaifH8+jHK9meV8Bdl3QgNsbCLu0/hLLFYWa3mtlC\nM1tkZrcVGv+5mS0Nxu8tNH6nmWWZ2TIzOy+MmkVEKkJsjPHnSzow+qxWPDfnW/7rufns238w7LL+\nQ4VvcZhZB2Ak0BPYB0w1szeBpsAAoLO77zWzhsHy7YEhwKlAY+B9M2vr7gcqunYRkYpgZvzivFOo\nm5LAXf9eQu6eAsZd052UhMoxSRTGFkc7YLa757n7fmAacBlwM/BXd98L4O6bg+UHAM+7+153Xw1k\nEQkdEZGoNqJvS+4d2InPsrZw9YTZ5OTtC7skIJzgWAj0NbP6ZpYCXEBka6NtMD7bzKaZWY9g+XSg\n8PFpa4MxEZGoNzijKY9e3Z1F63K5YtwsNuXuCbukig8Od18C3AO8C0wFFgAHiEyb1QN6A78ApphZ\nifcImdmNZpZpZpnZ2dllX7iISEjO73AiTwzvwdrteQwcO4Ovt+wOtZ5Qdo67+0R37+7u/YDtwHIi\nWxIve8Qc4CDQAFhHZIvkkCbB2OHvOd7dM9w9Iy0trfy/hIhIBTqtdQOeHdmbXXv2M3DsTJZsyA2t\nlrCOqjq047sZkf0bzwKvAmcF422BBGAL8DowxMwSzewkoA0wJ4y6RUTC1LlpKi+M6kN8rDF43Ewy\nv94WSh1hnQD4kpktBt4ARrt7DjAJaGlmC4HngWHB1sciYAqwmMjU1mgdUSUi1VXrhrV48ebTSKuZ\nyNCJs/lo2ebiX1TGrLKemXg8MjIyPDMzM+wyRETKzZZde7nu8Tks3bCT+wd3ZkCX4z9myMzmuXtG\nccup5YiISBXUoGYiz43sTffmdblt8gKenvVNhX22gkNEpIqqlRTPk9f35JxTGvLbVxfyjw9WVEh/\nKwWHiEgVlhQfy5ih3bmsazr3v7ecv7y1pNw/s3Kcvy4iIscsPjaG+wZ1JjUlgVZpNcv98xQcIiJR\nICbG+N1P21fMZ1XIp4iISNRQcIiISKkoOEREpFQUHCIiUioKDhERKRUFh4iIlIqCQ0RESkXBISIi\npRKV3XHNLBs4no5fDYhcC6S6qG7fF/Sdqwt959Jp7u7FXgkvKoPjeJlZZklaC0eL6vZ9Qd+5utB3\nLh+aqhIRkVJRcIiISKkoOI5sfNgFVLDq9n1B37m60HcuB9rHISIipaItDhERKRUFRyFmdr6ZLTOz\nLDP7Vdj1lAcza2pmH5nZYjNbZGa3BuP1zOw9M1sR/Kwbdq1lycxizWy+mb0ZPD7JzGYH63qymSWE\nXWNZM7NUM3vRzJaa2RIz6xPN69nMbg/+TS80s+fMLCka17OZTTKzzWa2sNDYEderRTwcfP8vzaxb\nWdSg4AiYWSzwT+AnQHvgSjOrmKuiVKz9wH+7e3ugNzA6+J6/Aj5w9zbAB8HjaHIrUPiamvcAf3f3\n1sB24IZQqipfDwFT3f0UoDOR7x+V69nM0oH/AjLcvQMQCwwhOtfzE8D5h40dbb3+BGgT3G4ExpRF\nAQqO7/UEstx9lbvvA54HBoRcU5lz9w3u/nlwfyeRXybpRL7rk8FiTwKXhFNh2TOzJsCFwITgsQFn\nAy8Gi0TV9wUwszpAP2AigLvvc/ccong9E7miabKZxQEpwAaicD27+3Rg22HDR1uvA4CnPGIWkGpm\njY63BgXH99KBNYUerw3GopaZtQC6ArOBE9x9Q/DURuCEkMoqDw8CvwQOBo/rAznuvj94HI3r+iQg\nG3g8mKKbYGY1iNL17O7rgPuAb4kExg5gHtG/ng852notl99rCo5qysxqAi8Bt7l7buHnPHKoXVQc\nbmdmFwGb3X1e2LVUsDigGzDG3bsCuzlsWirK1nNdIn9dnwQ0Bmrww+mcaqEi1quC43vrgKaFHjcJ\nxqKOmcUTCY1n3P3lYHjToU3Y4OfmsOorY6cDF5vZ10SmH88mMvefGkxpQHSu67XAWnefHTx+kUiQ\nROt6/hGw2t2z3b0AeJnIuo/29XzI0dZrufxeU3B8by7QJjgKI4HIjrXXQ66pzAXz+xOBJe7+QKGn\nXgeGBfeHAa9VdG3lwd3vdPcm7t6CyDr90N2vBj4CBgaLRc33PcTdNwJrzOzkYOgcYDFRup6JTFH1\nNrOU4N/4oe8b1eu5kKOt19eBa4Ojq3oDOwpNaR0znQBYiJldQGQ+PBaY5O5/DrmkMmdmZwCfAF/x\n/Zz/r4ns55gCNCPSWXiwux++A65KM7Mzgf/n7heZWUsiWyD1gPnAUHffG2Z9Zc3MuhA5ICABWAUM\nJ/LHYlSuZzP7I3AFkSMH5wMjiMznR9V6NrPngDOJdMHdBPweeJUjrNcgRB8hMm2XBwx398zjrkHB\nISIipaGpKhERKRUFh4iIlIqCQ0RESkXBISIipaLgEBGRUlFwSNQws/pmtiC4bTSzdYUeV8quqGZ2\nvZmdWI7vX8PMPjazmODxKWb2dtBF9XMze97MGppZFzObWF51SHSJK34RkarB3bcCXQDM7A/ALne/\nL9SiIrXEuvuBozx9PfA5kf5CJX2/uEL9l4ozAnjB3Q+aWQrwb+Dn7v5W8F7nAPXdfYGZtTSz9KDv\nk8hRaYtDqgUzG2Zmc4Ktj0fNLMbM4swsx8weCK7j8I6Z9TKzaWa2KjghFDMbYWavBOMrzOx/S/i+\nD5rZl0BPM/ujmc21yLUixgZn8l5BJOgmH9oqMrO1ZpYavHdvM3s/uH+XmT1lZp8BTwSf8UDw2V+a\n2YijfPWr+f4s4qHAtEOhAeDuH7j7oXbzbxI5gU6kSAoOiXpm1gG4FDjN3bsQ2dIeEjxdB3jb3U8F\n9gF/INKuYhDwp0Jv05NIq+ouwFXB1E5x7zvd3Tu5+0zgIXfvAXQMnjvf3ScDC4Ar3L1L0M6/KKcA\n57j7UCLXVtjs7j2BHkSuq9LssO+dBDRx97XBUAciHWOPJhPoW0wNIpqqkmrhR0R+uWZGOjCQzPet\npvPd/b3g/ldEevnsN7OvgBaF3uMdd98OYGavAmcQ+f9ztPfdB7xS6PXnmNkvgCQirSLmAW+X8nu8\n5u57gvvnAu3MrHBQtSHSs+mQhvzwug1F2Uyks6xIkRQcUh0Ykd5jv/2PwUjX1MJ/5R8E9ha6X/j/\nx+G9ebyY980P2lsT7Ft4BOjm7uvM7C4iAXIk+/l+JuDwZXYf9p1ucfcPjvI+APmHvccioFcRyycF\nrxEpkqaqpDp4HxhsZg3gu6OvmhXzmsOda5FreKcQue7DZ6V432QiQbTFzGoBlxd6bidQq9Djr4Hu\nwf3Cyx3uHeCWIKQws5PNLLnwAu6eTeSKeIeOKHsaONPMvrtOhZmdZWbtgodtgYWIFENbHBL13P2r\noHPq+8FhqQXAKGB9Kd5mLpGdzI2BJ919AXzXkbXI93X3rWb2JJE23xuIdCI+5HFggpnlE9mP8gfg\nMTPLAaYXUc84Ip1QFwTTZJs58qWO3wdOAz529zyLXNjq72b2j6DeBUSuxw5wFtHbdlzKkLrjihQj\nOGKpg7vfFnYtpWVmPYhMaQ0vZrlkIteuOL2IQ4dFAE1ViUQ1d58LfHroBMAiNAN+qdCQktAWh4iI\nlIq2OEREpFQUHCIiUioKDhERKRUFh4iIlIqCQ0RESkXBISIipfL/AV8gyDbwYb4CAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1115dd828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "([0.0, 4.0, 15.0, 20.0, 25.0, 37.0, 50.0, 100.0],\n",
       " [999.8425, 999.975, 999.1026, 998.2071, 997.0479, 993.3316, 988.04, 958.3665])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## STAGE 3: plotting\n",
    "# same as before but with the plotting lines of code in our function\n",
    "\n",
    "# remember to have %pylab inline for ipython\n",
    "# ... and to import pylab and numpy if not already imported\n",
    "from pylab import *\n",
    "from numpy import *\n",
    "\n",
    "def readTempDenFile(filename = \"data/density_air.dat\"):\n",
    "    infile = open(filename, \"r\")\n",
    "    temp = []\n",
    "    dens = []\n",
    "    for line in infile:\n",
    "        try:\n",
    "            t, d = line.split()\n",
    "            t = float(t)\n",
    "            d = float(d)\n",
    "        except:\n",
    "            continue\n",
    "        temp.append(t) # N.B. we're now filling out temp and dens lists\n",
    "        dens.append(d)\n",
    "    infile.close()\n",
    "    plot(array(temp), array(dens))\n",
    "    xlabel(\"Temperature (C)\")\n",
    "    ylabel(\"Density (kg/m^3)\")\n",
    "    show()\n",
    "    return temp,dens\n",
    "\n",
    "# run function\n",
    "readTempDenFile()\n",
    "readTempDenFile(\"data/density_water.dat\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Read acceleration data and find velocities**</br>\n",
    "A file data/acc.dat contains measurements $a_0, a_1, \\ldots, a_{n-1}$ of the acceleration of an object moving along a straight line. The measurement $a_k$ is taken at time point $t_k = k\\Delta t$, where $\\Delta t$ is the time spacing between the measurements. The purpose of the exercise is to load the acceleration data into a program and compute the velocity $v(t)$ of the object at some time $t$.\n",
    "\n",
    "In general, the acceleration $a(t)$ is related to the velocity $v(t)$ through $v^\\prime(t) = a(t)$. This means that\n",
    "\n",
    "$$\n",
    "v(t) = v(0) + \\int_0^t{a(\\tau)d\\tau}\n",
    "$$\n",
    "\n",
    "If $a(t)$ is only known at some discrete, equally spaced points in time, $a_0, \\ldots, a_{n-1}$ (which is the case in this exercise), we must compute the integral above numerically, for example by the Trapezoidal rule:\n",
    "\n",
    "$$\n",
    "v(t_k) \\approx v(0) + \\Delta t \\left(\\frac{1}{2}a_0 + \\frac{1}{2}a_k + \\sum_{i=1}^{k-1}a_i \\right), \\ \\ 1 \\leq k \\leq n-1. \n",
    "$$\n",
    "\n",
    "We assume $v(0) = 0$ so that also $v_0 = 0$.\n",
    "Read the values $a_0, \\ldots, a_{n-1}$ from file into an array, plot the acceleration versus time, and use the Trapezoidal rule to compute one $v(t_k)$ value. Experiment with different values of $\\Delta t$ and $k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8XOWV8PHf0ajZsrrkJlnFtty75YLpJWASgiEQAiF1\nISxJSEjbTfZ9d5MsWfYNm93sZrPJBqdAEnogAe+GugQwa+Mi9y43WcVFXVZvc94/ZqTIxmVszZ07\n5Xw/H32suXM1cwZGOnOf8zznEVXFGGOMAYhzOwBjjDHhw5KCMcaYQZYUjDHGDLKkYIwxZpAlBWOM\nMYMsKRhjjBlkScEYY8wgSwrGGGMGWVIwxhgzKN7tAC5UTk6OFhUVuR2GMcZElE2bNtWrau75zou4\npFBUVERZWZnbYRhjTEQRkSOBnGfDR8YYYwY5lhRE5FciUisiO89z3iIR6ROR252KxRhjTGCcvFJ4\nHFh+rhNExAM8ArzuYBzGGGMC5FhSUNXVQON5TvsS8AJQ61QcxhhjAudaTUFE8oBbgf90KwZjjDGn\ncrPQ/G/AN1XVe74TReQ+ESkTkbK6uroQhGaMMbHJzSmppcAzIgKQA3xQRPpU9cXTT1TVlcBKgNLS\nUtsqzhhjHOLalYKqFqtqkaoWAc8DXzhTQjDGRJ/efi9Pra+kor7d7VDMaRy7UhCRp4GrgBwRqQa+\nAyQAqOrPnHpeY0x4q23t4oEnt7ChopGURA//+JHZrJiX53ZYxs+xpKCqd13AuZ9xKg5jTPjYdKSJ\nLzy5iZbOXh5aMZNVW4/y4DNbee9gA9+9eSbJCR63Q4x5EdfmwhgTmQ7Xt3PXz9cxNi2ZP3xhMdPH\npXHX4gL+9Y1yfvr2QQ7VtfOrzy5iVJL9WXKTtbkwxoTEP7+2j/g44fn7L2H6uDQAEjxx/PXyafz4\nrvlsqmzik79cT0tnr8uRxjZLCsYYx22rauaPO45x7+UTGZ2W/L77Pzx3PD/5+AJ21rRw9y/W0dTe\n40KUBiwpGGMcpqo88upeslIS+dzlxWc9b/mssaz8ZCnlx9v4/it7QxihGcqSgjHGUe/ur2ftwQYe\nuHoyqckJ5zz36mmjub00nxe31tDcYVcLbrCkYIxxjNerfP+VveRnjuDupQUB/cynLimku8/Lc2VV\nDkdnzsSSgjHGMWsPNrD72Em+ct0UkuIDm246bWwaS4qz+O26I/R7rYFBqFlSMMY45nebqkhLjuem\nOeMu6Oc+vayIqsZO3t5nDZRDzZKCMcYRJ7t6eXXncW6eN/6CF6V9YMYYxqYl8/jaCmeCM2dlScEY\n44g/bj9Gd5+X2xdOuOCfTfDEcfeSAt7dX8/BujYHojNnY0nBGOOI5zdVUzJ6FHPz0y/q5+9cXECC\nR3h6fWWQIzPnYknBGBN0h+ra2HSkidsX5uNvj3/BclOTuLwkl1d3HUfVCs6hYknBGBN0L2yuJk7g\n1vnD6356w8wxVDd1svvYySBFZs7HkoIxJqj6vcrvN9dw5ZTcM7a0uBDXTh+DCLy+60SQojPnY0nB\nGBNUmyubONbSxUcW5A/7sXJGJVFamMnruy0phIolBWNMUP3v/nriBK4oyQ3K490wcyx7jp2kqrEj\nKI9nzs2SgjEmqNYerGdWXjrpI8/d5yhQ188YC8Bru44H5fHMuVlSMMYETXt3H1sqm1k2KSdoj1mQ\nPZJpY1OtrhAilhSMMUGzsaKRPq9y6eTsoD7u9TPHUnakkfq27qA+rnk/SwrGmKBZe7CBRE8cpYVZ\nQX3cG2aOwavwpz3WC8lpjiUFEfmViNSKyM6z3H+3iGwXkR0islZE5joVizEmNNYcqGd+QQYjEi+s\n19H5zBiXxoSsEfzX9qNBfVzzfk5eKTwOLD/H/YeBK1V1NvA9YKWDsRhjHNbU3sPuYye5dHLw6gkD\nRIRb5+Wx5kA9x1u6gv745s8cSwqquhpoPMf9a1W1yX9zHTD8Sc3GGNe8d6gBVYJeTxhw28J8vAq/\n31LtyOMbn3CpKdwDvHK2O0XkPhEpE5Gyurq6EIZljAnUmgP1pCR6mJOf4cjjF2ansLgoi+c3VVsv\nJAe5nhRE5Gp8SeGbZztHVVeqaqmqlubmBmdBjDEmuNYebGDJxGwSPM79WbltYR6H6trZUtXs2HPE\nOleTgojMAX4BrFDVBjdjMcZcvKPNnRyub2fZJGeGjgZ8cPY4khPieGGTDSE5xbWkICIFwO+BT6pq\nuVtxGGOGb8NhX/lw6URnk0JqcgLLZ45l1bajdPX2O/pcscrJKalPA+8BU0WkWkTuEZH7ReR+/ynf\nBrKBn4rIVhEpcyoWY4yzNlY0Miopnunj0hx/rtsXTqC1q483rEmeI+KdemBVves8998L3OvU8xtj\nQqesookFhZl44i5uQ50LccmkbMakJfHqzuN8eO54x58v1rheaDbGRLaWjl72nWhlUWFmSJ7PEycs\nLMxk59GWkDxfrLGkYIwZlk2VvnpCaVFwW1ucy4xxaRxp6KC1qzdkzxkrLCkYY4ZlY0UTCR5h3gRn\n1iecyczx6QDsPd4asueMFZYUjDHDUlbRyKy89KD3OzqXGeN9Be1dNTaEFGyWFMwZ1bZ28fXntvHr\ntRUca+l0OxwTprp6+9lW1cKiEA4dAYxOTSI7JZHdx06G9HljgWOzj0zkUlX+6nfbWb2/jhc2V/Od\nVbtYUJDBD++YR1FOitvhmTCys6aFnn4vpSEqMg8QEWaMT7Ok4AC7UjDv88T6St4pr+Ohm2fyP1+7\nkr+6YSoH69r50tNb6Onzuh2eCSMbK3w9LReGOCmAbwip/Hgbvf32ngwmSwrmFAfr2nj4j7u5Ykou\nn1hayOTRo/ji1ZN55LY57Khp4V/e2Od2iCaMlFU0Mik3hexRSSF/7hnj0ujp93Kgti3kzx3NLCmY\nQb39Xr727FaSEzz84PY5iPx5IdLyWWP5+JICVq4+xJoD9S5GacKF16uUHWkKeT1hwEx/sXn3URtC\nCiZLCmbQw3/cw7bqFh6+ZTZj0pLfd//ffWgGE3NS+OqzW2lq73EhQhNODtS10dLZG9L1CUMV54wi\nOSHO6gpBZknBAPDMhkoeX1vBPZcV86E54854zohEDz+6cz61rd08W1YV4ghNuNl8xL16AvhWNk8b\nm8YuW9kcVJYUDBsON/J3L+3kiim5/M2N08557qy8dGbnpfParuMhis6Eq61VzWSMTKAoe6RrMcwY\nn8buoydt050gsqQQ4463dPH5JzYxIXMkP75rPvEBbJCyfNZYtlQ22/qFGLe1qpm5+Rmn1J5Cbca4\nNE529VHTbO/FYLGkEOOeWn+Exo4eVn5qIekjEgL6meWzxgLw+i5rXRyr2rv7KD/RGtLWFmcyUGze\nZcXmoLGkEMNUlZe2HWXZpGwmj04N+Ocm5Y6iZPQoXtl5zMHoTDjbXt2CV2FegbtJYdrYNOLEZiAF\nkyWFGLatuoUjDR2smJd3wT9746yxbDjcSENbtwORmXC31b9H8rx8d5PCiEQPxTkpVmwOIksKMezF\nLTUkxscNDgddiBtmjcWr2O5XMWprVRNF2SPJTEl0OxTmTchkc2WzFZuDxJJCjOrr9/Lf249xzdTR\npCUHVksYasa4NCZkjeBVm4UUk7ZWNbteTxhQWpRJY3sPh+rb3Q4lKlhSiFHvHWqgvq2bFfMubjtD\nEeHGWeNYc6Celk7b6CSWHGvp5MTJ7vBJCv51Epv86ybM8DiWFETkVyJSKyI7z3K/iMi/i8gBEdku\nIgucisW834tbjpKaFM/V00Zf9GPcMHMsvf1qQ0gxZmulv55Q4M6itdNNyh1F+ogENlVYUggGJ68U\nHgeWn+P+G4ES/9d9wH86GIsZoqu3n9d2HWf5rLEkJ1z8xijzJ2RQnJPC0xsqgxidCXdbq5pJ9MQx\nfVzgM9acFOffs7nsSKPboUQFx5KCqq4GzvV/aQXwG/VZB2SIyJn7K5igentfLW3dfRc162iouDjh\n7iUFbDrSxB7rPxMztlQ1M318Gknxodtp7XwWFmZysK7denIFgZs1hTxgaAOdav8x47C399WRmhTP\n0onDb2R2+8J8kuLjeGLdkSBEZsJdX7+XHdUtzA+TesIAqysET0QUmkXkPhEpE5Gyuro6t8OJaKrK\nu/vrWTY5O6CWFueTMTKRm+aM58UtNbR19wUhQhPOyk+00dnbHzZF5gFzJ2SQ4BHKLCkMm5tJoQaY\nMOR2vv/Y+6jqSlUtVdXS3NzckAQXrQ7Vt1PT3MnlJcH77/iJpQW09/Tz4pYz/u8zUWRHja/IPDfM\nkkJygoeZ49PZZHWFYXMzKawCPuWfhbQUaFFV65vgsHfLfVdaVwQxKcybkMHM8Wk8se6ILSCKcjtq\nWkhNiqcwy73OqGdTWpjJtuoWuvv63Q4lojk5JfVp4D1gqohUi8g9InK/iNzvP+Vl4BBwAPg58AWn\nYjF/9u7+egqzR1IQxHbHIsInlhay93grmyvt8j2a7aw5ycy8NOLi3OuMejalRZn09HnZWWOTHobD\nydlHd6nqOFVNUNV8Vf2lqv5MVX/mv19V9YuqOklVZ6tqmVOxGJ+ePi/vHWrg8pKcoD/2zXPHkxgf\nxys7bIVztOrr97Ln2ElmjU93O5QzWljomzhhQ0jDExGFZhMcmyub6OjpD2o9YUBKUjylhZmsOdgQ\n9Mc24eFAXRvdfV5m5YVnUshNTaIweyRltohtWCwpxJB399fhiRMumZTtyONfOjmHPcdOUm+dU6PS\nwLDMrLw0lyM5u4WFmWw60mS1rWGwpBBD3t1fz4KCjItqgBeISyf7hqXes6uFqLSzpoWRiR6Kc0a5\nHcpZLSrKoqG9h8PWHO+iWVKIEY3tPeyoaXFk6GjA7Lx0UpPjWXOg3rHnMO7ZdbSFGePS8IRhkXnA\noiLfIjYbQrp4lhRixJoD9ajiSJF5gCdOWDoxmzUHLSlEm36vsuvoybCtJwyYlDuKzJEJ1gdpGCwp\nxIj1hxtISfQw2+Ff6ssm51DV2EllQ4ejz2NC63B9Ox09/YN7IocrEWFhYZZdKQyDJYUYsfFwEwsK\nM4PS2uJcLp3sK2Lb1UJ0GdjucnZ+eF8pgG+9wqH6dpvwcJEsKcSA5o4e9p1oZXHR8Bvgnc+k3FGM\nSUuyukKU2VHdQlJ8HJNzw7fIPMDqCsNjSSEGDHSOXFTsfFIQES6dlMN7Bxvwem1aYLTYebSFaePS\nHL/SDIZZeekkxsdRVmF1hYsR/v+HzbBtqGgkwSMh62y5bHIODe097D3eGpLnM87yepVdNSeZFeb1\nhAFJ8R7m5WdYx9SLZEkhBmw83MjsvPRh7bJ2IQbqCmutrhAVqpo6aO3uc3ySQjCVFmWys6aFzh5r\njnehLClEua7efnbUtIRk6GjAuPQRFGaPZMNhu3yPBjtqfEXmmWHa8+hMSosy6fMqW6ua3Q4l4sSf\n7wQRSQJuA4qGnq+qDzkXlgmWLZXN9PZrSIrMQy0qyuLNPSfwejUsO2qawG0+0kxSfBxTx4bHnsyB\nWFjge7+XVTQ61tYlWgVypfASvv2U+4D2IV8mAmysaEQESgtDmxQWF2XR1NHLwbq2kD6vCb7NlU3M\nzc8gMT5yBhbSRyYwdUwqG6zYfMHOe6UA5KvqcscjMY7YWNHI1DGppI90pt/R2Sz2D1dtqGikZEzk\nfMI0p+rq7WfX0Rb+4rJit0O5YJdMyuaZjZV09faHrJ4WDQJJ/WtFZLbjkZig6+v3svlIE4tCPHQE\nUJg9ktzUJDZaXSGi7TraQm+/sqAg0+1QLtiVU3Pp6vVabesCBZIULgM2icg+EdkuIjtEZLvTgZnh\n23Oslfae/pAWmQeICIuLsthoC4gi2sAal0hMCkuLs0mMj+Md/xa0JjCBJIUbgRLgeuDDwE3+f02Y\n2+gfTx1Y4Rlqi4oyqWnupLrJ+iBFqs1HminI8l31RZoRiR6WFGfx9r5at0OJKOdNCqp6BMjAlwg+\nDGT4j5kwt7WqmXHpyYxLH+HK8w9coWy0Yl9EUlU2VzaxoCA0ix6dcNXU0Rysa6eq0T6YBOq8SUFE\nHgSeBEb7v54QkS85HZgZvq1VzSFbxXwm08amkZocz4bDNoQUiWqaO6lt7WZBYeQNHQ24copv/5DV\n+20IKVCBDB/dAyxR1W+r6reBpcDnAnlwEVnur0UcEJFvneH+AhF5S0S2+OsVH7yw8M3ZNLb3UNnY\nwVwXk4InTigtzLQrhQgVyfWEAZNyU8jLGMHb+ywpBCqQpCDA0LXi/f5j5/4hEQ/wE3w1iRnAXSIy\n47TT/hZ4TlXnA3cCPw0kaHN+26p9Kznn5rt76b+oOIsDtW00WBvjiLOlspkRCR6mRdCitdOJCFdN\nzWXtgXp6+rxuhxMRAkkKjwHrReS7IvJdYB3wywB+bjFwQFUPqWoP8Ay+RXBDKTDQZSsdOBpQ1Oa8\ntlY2Eycwx+X+9wMrqW0WUuTZXNnE3AnpEdEZ9VyunJJLe0//4JWPObdACs0/BD4LNPq/Pquq/xbA\nY+cBVUNuV/uPDfVd4BMiUg28DFitIki2VTdTMjqVlKRA1ic6Z3Z+OknxcTZXPMJ09vSz++jJiB46\nGrBscg4JHuHtcpuFFIizJgURSfP/mwVUAE/4v474jwXDXcDjqpoPfBD4rYi8LyYRuU9EykSkrK7O\nxgbPR1XZVtXM3AnuNzBLivewoCCT9Ycb3A7FXIDt1c30eSNz0drpRiXFU1qYxepy69obiHNdKTzl\n/3cTUDbka+D2+dQAE4bczvcfG+oe4DkAVX0PSAbet7O8qq5U1VJVLc3NzQ3gqWNbZWMHTR29zJsQ\nHr/QSyZmsfvYSVo6e90OxQRoc6WvJjU/gqejDrVsUjZ7j5+kuaPH7VDC3lmTgqre5P+3WFUnDvkq\nVtWJATz2RqBERIpFJBFfIXnVaedUAtcCiMh0fEnBLgWGaaBdcDhcKQAsnZiNKtbyIoJsrmyiOCeF\n7FGRt2jtTBYXZ6FqW3QGIpB1Cm8Gcux0qtoHPAC8BuzBN8tol4g8JCI3+0/7OvA5EdkGPA18RlVt\nD8dh2lbVQnJCHFPCpBHdvAm+Dps2hBQZVJUtlU1Rc5UAMNfegwE7axVSRJKBkUCOiGTy52moaby/\nYHxGqvoyvgLy0GPfHvL9buDSC4zZnMfWqiZmjU8nIUxmjSQneJg/IYN1h+xKIRJUNXZS39bD/Cio\nJwxITvAwb0KGTXgIwLn+avwlvvrBNP+/A18vAf/hfGjmYvT2e9l59KSrK5nPZMnEbHYdbeFkl9UV\nwt3myoFFa+H1HhquJcVZ7Dx6krbuPrdDCWvnqin8SFWLgW8MqSUUq+pcVbWkEKb2HW+lp8/r6krm\nM1k6MQuv+nbCMuFtc2UTIxM9TA2T4cdgWVKcTb9Xbb3CeQSyTuHHIjJLRO4QkU8NfIUiOHPhBorM\n4XalsKAgk0RPHOttCCnsDey0FumL1k63oDCD+Dhhg9UVzimQQvN3gB/7v64G/gm4+Zw/ZFyzraqZ\nrJRE8jPd6Yx6NgNjuusO2S9kOOvo6WPPsVYWFIbXh4pgGJkYz6y8dKsrnEcgHwVuxzdt9LiqfhaY\ni68lhQlD26tbmJOfjsh521OF3JKJvjHdVqsrhK3t1S30R8mitTNZMjGLbVUtdPX2n//kGBVIUuhU\nVS/Q51/lXMupi9JMmGjv7mN/bStzXG6CdzZLJ/rGdMtsTDdsDRSZo2nm0VBLirPo6feyxb84z7xf\nIEmhTEQygJ/jm320GXjP0ajMRdlZ04JXYV6YLFo73YKCTBI8YnWFMLb5SDPFOSlkpSS6HYojFhZm\nIYINIZ3DObuliW8M4v+pajPwMxF5FUhTVdujOQxtr24BCNsrhRGJHmaMTx/8NGrCy8CitSunRm8r\nmfQRCcwYl+ZfxFbidjhh6ZxXCv7VxS8PuV1hCSF8batuJi9jBDlh3Jpg/oQMdlS30Ndvve3DTWVj\nBw3tPVFbTxiwbFI2ZRVNVts6i0CGjzaLyCLHIzHDNlBkDmfzCzLo7O1n34lWt0Mxp/nzorXoTgo3\nzBxLT7+Xt2w3tjMKJCksAd4TkYP+LTN3iIhdLYSZJv/2m+E6dDRg4A+OFfrCz8aKJlKT4pkawTut\nBWJBQSa5qUm8tvO426GEpUB2YLnB8SjMsP15+83wvlLIzxxBdkoiW6ua+cTSQrfDMUNsONxIaVEm\nnrjwm84cTHFxwg0zx/D7zTV09faTnOBxO6SwEsiK5iP4pqBe4/++I5CfM6G1vboFEZgV5klBRJhf\nkMEWKzaHlfq2bg7UtrG4ONvtUEJi+cxxdPT0s7rchpBOF+iK5m8Cf+M/lIBvBzYTRrZXNzMxJ4W0\n5AS3Qzmv+QWZHKxrp6XDCn3hYmCK5pKJwdpUMbwtmZhF+ogEXrUhpPcJ5BP/rfjaWrQDqOpRILoH\nHSOMqrKtuoW5YV5PGDDQl2lrtdUVwsWGw42MSPAwOy+8rzSDJcETx3XTx/A/e07Q02cz4YYKJCn0\n+KemKoCIpDgbkrlQx092UdfaHfYzjwb42nDAVis2h411hxpYWJgZNntwhMKNs8ZysqvP+nGdJpB3\nwHMi8iiQISKfA/4H3+pmEya2VfkXrYVZZ9SzSU1OYMroVLZUWV0hHDR39LDvRCuLi2Nj6GjAZSU5\njEz08IoNIZ0ikELzPwPPAy8AU4Fvq+qPnQ7MBG5bdTPxccKMcWluhxIwX7G5Gdt91X0bK5pQ9fUF\niiXJCR6unjaa13cdp9cWUw4K6FpRVd9Q1b9S1W+o6htOB2UuzPbqZqaNS42oqXXzJmTQ0tnL4fp2\nt0OJeesPNZAYHxd2GzOFwm0L8mho7+HNPbVuhxI2zpoURKRVRE6e4atVRE6GMkhzdl6vsr0qcorM\nAwa6cA5sCmTcs6GikXkTMiLqQ0WwXFGSy9i0ZJ7ZWOl2KGHjXNtxpqpq2hm+UlU1oHEKEVkuIvtE\n5ICIfOss59whIrtFZJeIPHWxLyRWHapvp7W7L+I+5U0ePYpRSfG2NaLLWrt62VnTEnNDRwPiPXF8\ntDSfd8rrONrc6XY4YSGg4SMRuUxEPuv/PkdEigP4GQ/wE+BGYAZwl4jMOO2cEnzrHy5V1ZnAVy4w\n/pi3LUy33zwfT5ywsDDTWhi7bNORJrzq2784Vt1R6tse5ndl1S5HEh4uZvFaIoEtXlsMHFDVQ6ra\nAzwDrDjtnM8BP1HVJgBVtYG9C7StupmURA+Tcke5HcoFWzIxi/21bdS3dbsdSsx6aetREuPjonL7\nzUBNyBrJZZNzeK6sin6vTXxwcvFaHlA15Ha1/9hQU4ApIrJGRNaJyPIzPZCI3CciZSJSVldny9KH\n2lbVzOz89IjsV7N0ou/TqV0tuGPTkUb+sKWGey8rZmRiIG3Qotediwqoae7kfw/Uux2K69xevBaP\nb6eLq4C7gJ/7d3k7haquVNVSVS3NzY3eDUAuVHdfP3uOtUZcPWHA7Lx0RiZ6WG+Lh0Ku36t8d9Vu\nxqQl8cWrJ7sdjuuumzGarJREnrWCs6OL12o4dS/nfP+xoaqBVaraq6qHgXJsO6SA7T3WSk+/N+Jm\nHg1I8MSxsDCTdbY9Z8j9rqyKHTUt/J8PTiclKbavEgCS4j3ctiCP13ed4HhLl9vhuMrJxWsbgRIR\nKRaRROBOYNVp57yI7yoBEcnBN5x0KODoY9xgu+wIvVIA3xDSvhOtNLb3uB1KzGjp7OWfXttHaWEm\nN88d73Y4YeNTlxThVeU371W4HYqrAik0FwPvDixeA/5XRIrO93Oq2gc8ALwG7AGeU9VdIvKQiNzs\nP+01oEFEdgNvAX+lqjaWEKCtVc3kjEpifHqy26FctKX+rpwbDtv/9lB5fE0FTR09fPfmmfi2YTfg\nKzhfP2MsT22opKOnz+1wXBPI8NHvgKFrwPv9x85LVV9W1SmqOklVH/Yf+7aqrvJ/r6r6NVWdoaqz\nVfWZC30BsWxbVTPzJqRH9C/27LwMRiR4bAgphMqONDJjXBqzYqQj6oW49/Jimjt6eWHz6SPdsSOQ\npBDvn1IKgP/7ROdCMoE42dXLwbr2iK0nDEiMH6gr2JVCKKhqROzl7ZaFhZnMzU/nsf89jDdGp6cG\nkhTqhgz3ICIrAJu35bId1b7OqJFcTxiwdGIWe4+30mR1BcdVNXbS0tnL7LzIf984QUT4i8uKOVTf\nztvlsblsKpCkcD/wf0SkUkQq8S1k+0tnwzLnM9AzKBo+8S3xr1dYb+sVHDcwOSEa3jdO+eDscYxL\nT+YX7x6OyS6+gcw+OqiqS/G1qpihqstU9YDzoZlz2X30JBOyRpAxMvJH8ubkp5OcEGdDSCGwo6aF\nxPg4poyxzRPPJsETx72XT2TtwQZ+tabC7XBCLpDZR/8oIhmq2qaqbSKSKSL/EIrgzNmVn2hl6pjI\n2T/hXJLiPSwqymKNrSZ13PbqZqaPSyMxPnZ2WLsYn11WxPKZY/mHP+7mf3afcDuckArknXGjqg72\nN/b3KfqgcyGZ8+np83K4vp0pYyKv39HZXFGSy/7aNutU6SCvV9lZc5I5NuvovOLihH/92Dxm56Xz\n5We2sLOmxe2QQiaQpOARkaSBGyIyAkg6x/nGYYfr2+nzalQNAVw+JQeAd/dbbyunHKpvp627j9lW\nTwjIiEQPv/hUKekjErjn1xupiJENoQJJCk8Cb4rIPSJyL/AG8GtnwzLnUn6iFSCqksLUMamMSUti\n9X4bQnLKjhr/CvgIn8YcSqPTknnss4vo7Vc+tvI9Dta1Be2xG9t7wrKQHUih+RHgH4Dp+NpcvAYU\nOhyXOYfyE63ECUzMDWZvQneJCJeX5PK/++utfbFDtlW1MCLBw6Qoet+EwrSxaTz9uaX0e5U7V67j\nQG3rsB/zj9uPUfoPb/AXj2+ktvX9vZb6+r3srGnh2Y2VVDZ0DPv5LkSgnbBO4OuS+lHgML4+SMYl\n5SdaKcpJibrtEy8vyeH5TdVsr24e3K7TBM+OmhZmjk8j3mNF5gs1dWwqz9y3lLt+vp47Hl3H3900\nnVvm5V1UN4ENhxv56nNbmZg7irUHG1j+b+/y8C2zSB+ZQFlFExsrGtlS2Uxbt6/VRkqih4dvnc0t\n80/fecD5d169AAAXTUlEQVQZZ00KIjIFXzvru/AtVnsWEFW9OiSRmbMqP9HG1CgaOhpweUkuIvDu\n/npLCkHW1+9l19EWPr7YLvIv1uTRqTx731K++tw2vvrsNp5cV8l3b555Qe1CDtS28rnflJGfOYLf\n/eUl1Ld18+AzW/n8k5sBEPENpd4yfzyLirIozknhH/57D195diur99fx0IpZjHK4q+25Hn0v8C5w\n08C6BBH5qqPRmPPq6u3nSEM7H47C7pZZKYnMzktndXkdX77WOqgH04G6Nrp6vbZobZgm5o7iD59f\nxvObqnnk1b2s+MkaHr5lFncuLjjlvJ4+L82dPTR39FLf1s3h+nYO1Lbxyo7jJHji+PVnF5OZkkhm\nSiIvfvFS/mvbUbJGJbKgIJP0EQmnPNZTn1vCf7x1gH9/cz8jEnxXDU46V1L4CL5212+JyKv4ttOM\n3M5rUeJgXRteJaqmow51RUku//nOQU529ZKWnHD+HzAB2V7lm1JpM4+GLy5OuGPRBG6YNZYvPb2F\nb/1+BydOdvPlaydz/GQX//7mfn5XVk3fabWxkYkepo5N5XsrZjEha+Tg8cT4OG5bmH/W54v3xPGV\n66Zw6eQcJuY4Xw86a1JQ1ReBF/07ra0AvgKMFpH/BP6gqq87Hp15n4GZR9E4fAS+usJ/vHWAtQca\nWD5rrNvhRI11hxtITYqnONuKzMGSPiKBX366lG++sJ1//Z9y3jtUz5bKZryq3LFoAtPHpZExIoGs\nlESKclIYl5ZM3DC2zV1UlBXE6M/uvINTqtoOPAU8JSKZ+IrN3wQsKbig/EQbCR6hKASfGNywoDCT\nUUnxrN5fZ0khSFo6enl5xzFunZ83rD9K5v0SPHH8y0fnMjYtmZ+9c5Bb5ufx1eumnHIlEGkuqGLh\nX8280v9lXFB+vJWJOaNIiNIZJAmeOC6ZlM3q8jpUNaL3iggXv99STVevl7uXWJHZCSLCXy+fxpev\nLYmKGYHR+ZclipXXtlISpfWEAVeU5FDd1ElFiOdnRyNV5cn1lcydkGGb6jgsGhICWFKIKB09fVQ1\ndkZtPWHAFVNyAWt5EQwbDjdyoLaNu5cUnP9kY7CkEFH2n/AtsS+J8qRQmJ1CQdZIVpdbUhiuJ9ZX\nkpYcz4fnRN8UZuMMR5OCiCwXkX0ickBEvnWO824TERWRUifjiXSDM4/GRndSALhiSg7vHWygp897\n/pPNGdW3dfPqzmPctjCfEYnRMbRhnOdYUhARD/AT4EZ8G/TcJSIzznBeKvAgsN6pWKJF+YlWkuLj\nKIjgmQ2Burwkl/aefjZXNrkdSsR6rqyK3n61oSNzQZy8UlgMHFDVQ6rag2/x24oznPc94BHg/V2h\nzCBV5a19dczNz8ATA9MKl03KxhMnVle4SO3dffzy3cNcXpLD5NHRf2VpgsfJpJAHVA25Xe0/NkhE\nFgATVPWPDsYRFbZXt3Cgto1bF4SmKZbbUpMTWFCQwepya6V9MR5bc5iG9h6+9oEpbodiIoxrhWYR\niQN+CHw9gHPvE5EyESmrq4vNT44vbK4mKT6OD80Z53YoIXN5SS47j7bQ0NbtdigRpaWjl0dXH+K6\n6WOssaC5YE4mhRpgwpDb+f5jA1KBWcDbIlIBLAVWnanYrKorVbVUVUtzc3MdDDk8dff1s2rbUa6f\nOTam+gFdMSUXVVhzsMHtUCLKo6sP0tbdx9evt6sEc+GcTAobgRIRKRaRRHzN9VYN3KmqLaqao6pF\nqloErANuVtUyB2OKSH/aU0tzRy+3n6NpVjSanZdOxsgE3t5X63YoEaOutZvH1lRw05zxTB+X5nY4\nJgI5lhRUtQ94AN9ObXuA51R1l4g8JCI3O/W80eiFzdWMSUvissk5bocSUp44YfnMsfxx+zFqT9o8\nhED8/N1D9PR7+ep11nrcXBxHawqq+rKqTlHVSar6sP/Yt1V11RnOvcquEt6vvq2bt/fVccv8vJiY\ndXS6z181iT6vsnL1IbdDCXuqyn9vO8o100YzMTe6W6EY59iK5jD30taj9HmV2xfE1tDRgMLsFFbM\nG88T649QbwXnc9p3opWjLV1cO22026GYCGZJIcy9suMYM8enRX1ri3N54OrJ9PR5+fm7drVwLn/a\n66u9XG1JwQyDJYUw1tLZy5aqZq6J8V/yibmjuHnueH773hEa23vcDidsvbW3lpnj0xiTlux2KCaC\nWVIIY2sP1NPvVa6cEnvTcE/3wDWT6eztt6uFs2ju6GHTkaaY/wBhhs+SQhh7p7yO1OR45k3IcDsU\n100encqHZo/jN2sr7GrhDN4pr8OrNnRkhs+SQphSVVaX13HZ5Bzio3SXtQv14LUldNjVAuAbWjzZ\n1Tt4+629tWSlJDI33z5AmOGxvzZh6kBtG0dbugY3nDG+fSRumjOeX6+tiOnWF6vL67j6n9/m+h+u\nZt/xVvq9yjvldVw1JTcmpy2b4LKkEKbe8W8wY0nhVA9e66strIzBq4V+r/LDN8r59GMbyBmViKLc\n/rO1rFx9iKaOXhs6MkFhSSFMvVNex+TRo8jLGOF2KGFl8uhUbp47nt+sja11C1sqm/jYo+/x72/u\n5yPz83npi5fx+y9cyujUJB55dS+eOLEPECYoLCmEoa7efjYcbuSKEvslP5MvX1tCd18/j75z0O1Q\nHFfZ0MEXntzErT9dy+H6dn5w+xz++aNzGJHoIS9jBM/fv4xLJ2dz46yxpI+InWaJxjnxbgdg3m/d\noQa6+7xcOdWSwplMyh3FRxbk89iaClbMy2NWXrrbITnC61Xu/uU6Gtp6ePDaEj53xURGJZ36K5uZ\nksiT9y51KUITjexKIQytLq8nKT6OJcVZbocStv72Q9PJSknka89tpau33+1whq2334uqnnJs3aEG\nqho7+f5tc/jqB6a8LyEY4wRLCmFGVXljz3GWTswmOcE2Wz+bjJGJ/NPtcyg/0cYP3yh3O5wL1tPn\nZdORRn785n4+9uh7zPj2q3z/1b2nnPPC5hpSk+K5fsYYl6I0scg+eoSZXUdPUtXYyQNXT3Y7lLB3\n1dTR3L2kgJ+/e4hrpo1m6cRst0M6p86efn67roJ3yuvYfKSZzt5+RGDm+DRmjEvj8TUV3HNpMaPT\nkmnv7uOVnce4ee54+3BgQsqSQpj5445jeOKE62eMdTuUiPB/PzSdNQfq+eqzW3nxi5eGTd+f9u4+\nkhM8eOLEd/W3+wR//1+7qWnuZPq4ND62aAJLirNYMjGbrJREKhs6uOZf3uY/3znIdz48k9d2Haej\np5/bYmxjJeM+SwphRFV5eccxlk3KJjMl0e1wIsLIxHh+cvcC7vjZe3z2sY08d/8lIR97V1WqGjvZ\nWNFI2ZFGNlY0caC2jQSPkJcxgpGJ8ew+dpKpY1J55r6lZ7yiKcgeyUcW5PHU+ko+f+UkXthczYSs\nEZQW2h7LJrQsKYSRXUdPcqShg/uvnOR2KBFl5vh0fvqJhfzF4xv5/BOb+NVnFpEQgtYg7d19fGfV\nLlaX11Hb6lszkZYcT2lRFjfPHU9HTz9VTR3Unuzibz80nU8vKzpnXA9cXcILm2v4zqpdrD3YwJev\nKUHEViib0LKkEEZe2ekbOrphpg0dXagrp+Ty/26dzV+/sJ0Hn9nCt5ZPpyB7ZFAeu727j2MtXUwe\nfepuZk+tr+T5TdXcNGccSyZmU1qYydQxqcRdZKuJguyR3LYgj+fKqgG4LUY3VjLusqQQJnxDR8e5\nxD/GbC7cHYsmUN/ezQ9fL+eVnce5bvoYPrm0kKUTs0mMf/8n9BMnu/jHl/fw9r46rpqay63z81g2\nKYfKxna2VbWwrbqZzZVN7Dnm6y/06CcXDibsvn4vj605zJLiLP7j4wuC9hoGrhYWFGQELakZcyEs\nKYSJPcdaOVzfzr2XF7sdSkT7wlWT+cj8fJ5Yd4SnNlTyxu4TpCR6uKwkhyXF2eSkJpE5MoFdR0/y\n4zf30+tVrp02mrf31fHS1qPECXj9ywVGJnqYNyGDL1w1iZd3HOMfX97DVVNzSYr38PLO4xxt6eLv\nV8wKavwF2SP56d0LKMiyhGDc4WhSEJHlwI8AD/ALVf3+afd/DbgX6APqgL9Q1SNOxhSuXt5xjDjB\nho6CYGx6Mt+4YSoPXDOZd/fX89a+Wt7aW8tru06cct5100fzdzfNoDA7he6+ft7aW8eWyiYmjx7F\n3AkZTModNdh1dGFhJp95bCO/WXuEey8v5hfvHqI4J8WR/ZDtPWDc5FhSEBEP8BPgA0A1sFFEVqnq\n7iGnbQFKVbVDRD4P/BPwMadiClf/te0oK989xOUlueSMSnI7nKiRnODhAzPG8IEZY1BV6tt6aO7o\nobmzl6T4OOYM2XsgKd7D8lljWT7rzH+Qr5o6miun5PLvf9rPhKyRbK9u4Xu3zLro+oEx4crJKRqL\ngQOqekhVe4BngBVDT1DVt1S1w39zHRBTlTVV5cdv7udLT29hbn46//qxeW6HFLVEhNzUJErGpLKo\nKOuUhBCov/3QdDp6+vny01vIGJnA7VYINlHIyaSQB1QNuV3tP3Y29wCvOBhP2Hn4j3v4lzfKuXV+\nHk/cu8QKzGGuZEwqH19cQE+/l08sKWREoq00NtEnLArNIvIJoBS48iz33wfcB1BQUBDCyJzT1dvP\nE+uPsGLeeH54x1ybjx4hvn79FBLj47jnMpsQYKKTk1cKNcCEIbfz/cdOISLXAf8XuFlVz7hriqqu\nVNVSVS3NzY2OdtIbDjfS1evllnl5lhAiSMbIRP7uphm24txELSeTwkagRESKRSQRuBNYNfQEEZkP\nPIovIdQ6GEvYeae8jsT4uLBv4maMiS2OJQVV7QMeAF4D9gDPqeouEXlIRG72n/YDYBTwOxHZKiKr\nzvJwUeftfbUsKc6ycWljTFhxtKagqi8DL5927NtDvr/OyecPV1WNHRysa+fjSwrdDsUYY05hm+y4\n4J3yOsDXr8cYY8KJJQUXvL2vjvzMEUzKTXE7FGOMOYUlhRDr6fOy9mA9V07JtVlHxpiwY0khxMoq\nGuno6eeqqcHvmWOMMcNlSSHE3i6vI8EjLJtkU1GNMeHHkkIIdfX288rOYywqyiIlxFtGGmNMICwp\nhNAjr+6lqrGTL1w12e1QjDHmjCwphMiaA/U8tqaCzywr4rKSHLfDMcaYM7KkEAItnb1843fbmJib\nwjeXT3M7HGOMOSsb2HaYqvLtl3ZS29rNC59fZm0tjDFhza4UHKSqfP/Vvby09ShfubaEeRMufGMX\nY4wJJUsKDvrp2wd59J1DfGJpAQ9cY8VlY0z4s+GjIDpxsouDdW2c7Oxje3UzP337ILfMG89DN8+y\n1cvGmIhgSSFIdlS38NFH19LV6x08dsPMMfzgo3Ntc3djTMSwpBAEda3d3PfbMrJTknjktjlkpiSQ\nlpxAfuYIu0IwxkQUSwrD1Nvv5YtPbqapo4fn71/GrLx0t0MyxpiLZklhGNq6+/jef+1mQ0UjP7pz\nniUEY0zEs6RwFvuOt/LNF7Zz1dRc7r18IqP8vYpUle3VLTyzsYpVW2to7+nnL6+cyIp5eS5HbIwx\nw2dJ4QwO1bVx9y/W09nTx9aqZn7z3hHuuayYutZu3th9gprmTpIT4vjwnPHctaSA+bb+wBgTJSwp\nnKayoYOP/3w9qspLD1xKW3c/P3htLz94bR9J8XFcXpLLg9eWsHz2WNKSE9wO1xhjgsrRpCAiy4Ef\nAR7gF6r6/dPuTwJ+AywEGoCPqWqFkzGdSWdPP1uqmth4uIlnN1bS1dfP059byuTRqQA8ee9SDtW1\nMTY9mZGJlkeNMdHLsb9wIuIBfgJ8AKgGNorIKlXdPeS0e4AmVZ0sIncCjwAfcyqmAT19XlaX17Gh\nopGNFY3srGmht18RgRnj0vjZR2YzfVzaKT8zMXeU02EZY4zrnPzYuxg4oKqHAETkGWAFMDQprAC+\n6//+eeA/RERUVZ0IqLOnn2c2VvLz1Yc42tJFoieOOfnp3HPZRBYXZ7KwIIv0kTYkZIyJXU4mhTyg\nasjtamDJ2c5R1T4RaQGygfqhJ4nIfcB9AAUFBRcVzJ/2nuAbv9tOY3sPi4oy+d4ts7h0cg7JCda1\n1BhjBkTEALmqrgRWApSWll7UVURRdgrzJmTw+asmsagoK6jxGWNMtHAyKdQAE4bczvcfO9M51SIS\nD6TjKzgH3cTcUfzqM4uceGhjjIkaTrbO3giUiEixiCQCdwKrTjtnFfBp//e3A39yqp5gjDHm/By7\nUvDXCB4AXsM3JfVXqrpLRB4CylR1FfBL4LcicgBoxJc4jDHGuMTRmoKqvgy8fNqxbw/5vgv4qJMx\nGGOMCZztvGaMMWaQJQVjjDGDLCkYY4wZZEnBGGPMIEsKxhhjBkmkLQsQkTrgyEX+eA6ntdCIAfaa\nY4O95tgwnNdcqKq55zsp4pLCcIhImaqWuh1HKNlrjg32mmNDKF6zDR8ZY4wZZEnBGGPMoFhLCivd\nDsAF9ppjg73m2OD4a46pmoIxxphzi7UrBWOMMecQM0lBRJaLyD4ROSAi33I7HieIyK9EpFZEdg45\nliUib4jIfv+/mW7GGGwiMkFE3hKR3SKyS0Qe9B+P2tctIskiskFEtvlf89/7jxeLyHr/e/xZf8v6\nqCEiHhHZIiL/7b8d7a+3QkR2iMhWESnzH3P8fR0TSUFEPMBPgBuBGcBdIjLD3agc8Tiw/LRj3wLe\nVNUS4E3/7WjSB3xdVWcAS4Ev+v/fRvPr7gauUdW5wDxguYgsBR4B/lVVJwNNwD0uxuiEB4E9Q25H\n++sFuFpV5w2Zhur4+zomkgKwGDigqodUtQd4BljhckxBp6qr8e1LMdQK4Nf+738N3BLSoBymqsdU\ndbP/+1Z8fzTyiOLXrT5t/psJ/i8FrgGe9x+PqtcsIvnAh4Bf+G8LUfx6z8Hx93WsJIU8oGrI7Wr/\nsVgwRlWP+b8/DoxxMxgniUgRMB9YT5S/bv9QylagFngDOAg0q2qf/5Roe4//G/DXgNd/O5vofr3g\nS/Svi8gmEbnPf8zx97Wjm+yY8KKqKiJROd1MREYBLwBfUdWTvg+SPtH4ulW1H5gnIhnAH4BpLofk\nGBG5CahV1U0icpXb8YTQZapaIyKjgTdEZO/QO516X8fKlUINMGHI7Xz/sVhwQkTGAfj/rXU5nqAT\nkQR8CeFJVf29/3DUv24AVW0G3gIuATJEZOCDXjS9xy8FbhaRCnxDv9cAPyJ6Xy8Aqlrj/7cWX+Jf\nTAje17GSFDYCJf7ZCon49oJe5XJMobIK+LT/+08DL7kYS9D5x5Z/CexR1R8OuStqX7eI5PqvEBCR\nEcAH8NVS3gJu958WNa9ZVf9GVfNVtQjf7+6fVPVuovT1AohIioikDnwPXA/sJATv65hZvCYiH8Q3\nLukBfqWqD7scUtCJyNPAVfg6KZ4AvgO8CDwHFODrLnuHqp5ejI5YInIZ8C6wgz+PN/8ffHWFqHzd\nIjIHX5HRg++D3XOq+pCITMT3SToL2AJ8QlW73Ys0+PzDR99Q1Zui+fX6X9sf/DfjgadU9WERycbh\n93XMJAVjjDHnFyvDR8YYYwJgScEYY8wgSwrGGGMGWVIwxhgzyJKCMcaYQbai2Ziz8E//e9N/cyzQ\nD9T5b3eo6jJXAjPGQTYl1ZgAiMh3gTZV/We3YzHGSTZ8ZMxFEJE2/79Xicg7IvKSiBwSke+LyN3+\n/Q52iMgk/3m5IvKCiGz0f13q7isw5swsKRgzfHOB+4HpwCeBKaq6GF+b5y/5z/kRvt7/i4Db/PcZ\nE3aspmDM8G0caGcsIgeB1/3HdwBX+7+/DpgxpHtrmoiMGrIvgjFhwZKCMcM3tN+Od8htL3/+HYsD\nlqpqVygDM+ZC2fCRMaHxOn8eSkJE5rkYizFnZUnBmND4MlAqIttFZDe+GoQxYcempBpjjBlkVwrG\nGGMGWVIwxhgzyJKCMcaYQZYUjDHGDLKkYIwxZpAlBWOMMYMsKRhjjBlkScEYY8yg/w+1cB5OFLk2\nKAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x111511eb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## PART 1: read and plot\n",
    "# specify libraries\n",
    "%pylab inline\n",
    "from pylab import *\n",
    "from numpy import *\n",
    "\n",
    "dt = 0.5\n",
    "\n",
    "# read in acceleration\n",
    "infile = open(\"data/acc.dat\", \"r\")\n",
    "alist = []\n",
    "for line in infile:\n",
    "    alist.append(float(line))\n",
    "infile.close()\n",
    "aarray = array(alist)\n",
    "time = array([e*dt for e in range(len(alist))]) # time is specified by dt and the number of elements in acc.dat\n",
    "\n",
    "# plot\n",
    "plot(time, aarray)\n",
    "xlabel(\"Time\")\n",
    "ylabel(\"Acceleration\")\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0820085591101\n"
     ]
    }
   ],
   "source": [
    "## PART 2: trapezoidal function\n",
    "\n",
    "def trap(dt, k, alist):\n",
    "    if not (1 <= k <= (len(alist) - 1)):\n",
    "        raise ValueError\n",
    "    return dt*(.5*alist[0] + .5*alist[k] + sum(alist[:k]))\n",
    "\n",
    "dt = 2\n",
    "k = 4\n",
    "print(trap(dt, k, alist))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Make a dictionary from a table**</br>\n",
    "The file *data/constants.txt* contains a table of the values and the dimensions of some fundamental constants from physics. We want to load this table into a dictionary *constants*, where the keys are the names of the constants. For example, *constants['gravitational constant']* holds the value of the gravitational constant (6.67259 $\\times$ 10$^{-11}$) in Newton's law of gravitation. Make a function that reads and interprets the text in the file, and thereafter returns the dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'speed of light': 299792458.0, 'gravitational constant': 6.67259e-11, 'Planck constant': 6.6260755e-34, 'elementary charge': 1.60217733e-19, 'Avogadro number': 6.0221367e+23, 'Boltzmann constant': 1.380658e-23, 'electron mass': 9.1093897e-31, 'proton mass': 1.6726231e-27}\n"
     ]
    }
   ],
   "source": [
    "def read_constants(file_path):\n",
    "    infile = open(file_path, \"r\")\n",
    "    constants = {} # An empty dictionary to store the constants that are read in from the file\n",
    "    infile.readline(); infile.readline() # Skip the first two lines of the file, since these just contain the column names and the separator.\n",
    "    for line in infile:\n",
    "        words = line.split() # Split each line up into individual words\n",
    "        dimension = words.pop() # pop is a list operation that removes the last element from a list and returns it\n",
    "        value = float(words.pop()) # Again, use pop to obtain the constant itself.\n",
    "        name = \" \".join(words) # After the two 'pop' operations above, the words remaining in the 'words' list must be the name of the constant. Join the individual words together, with spaces inbetween, using .join.\n",
    "        constants[name] = value # Create a new key-value pair in the dictionary\n",
    "    return constants\n",
    "\n",
    "print(read_constants('data/constants.txt'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Explore syntax differences: lists vs. dictionaries**</br>\n",
    "Consider this code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{0: -5, 1: 10.5}\n"
     ]
    }
   ],
   "source": [
    "t1 = {} # this specifies a dictionary\n",
    "t1[0] = -5 # create a new key called 0, give it the value -5\n",
    "t1[1] = 10.5 # create a new key called 1, give it the value 10.5\n",
    "\n",
    "print(t1) # shows the contents of 0 and 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explain why the lines above work fine while the ones below do not:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "list assignment index out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-22-f5b601b7e131>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mt2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;31m# this specifies a list\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mt2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0mt2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m10.5\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;31m# cannot work, raises and IndexError\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;31m# [] notation is for indexing in lists\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: list assignment index out of range"
     ]
    }
   ],
   "source": [
    "t2 = [] # this specifies a list\n",
    "t2[0] = -5\n",
    "t2[1] = 10.5\n",
    "# cannot work, raises and IndexError\n",
    "# [] notation is for indexing in lists\n",
    "# this is trying to specify the 0th element is -5\n",
    "# IndexError is raised because the intial list has no dimensions, and therefore [0] is out of index\n",
    "# Dictionaries behave differently, the [] notation will create a new key and add the specified element to the dictionary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What must be done in the last code snippet to make it work properly?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-5, 10.5]\n",
      "[-5, 10.5]\n"
     ]
    }
   ],
   "source": [
    "## To make it work properly we have two options\n",
    "\n",
    "# 1. specify the list's dimension before asignment\n",
    "t2 = [0]*2 # a list containing two 0s\n",
    "t2[0] = -5\n",
    "t2[1] = 10.5\n",
    "print(t2)\n",
    "del t2\n",
    "\n",
    "# 2. start with a dimensionless list, but append new elements\n",
    "t2 = []\n",
    "t2.append(-5)\n",
    "t2.append(10.5)\n",
    "print(t2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Compute the area of a triangle**</br>\n",
    "An arbitrary triangle can be described by the coordinates of its three vertices: $(x_1, y_1), (x_2, y_2), (x_3, y_3)$, numbered in a counterclockwise direction. The area of the triangle is given by the formula:\n",
    "\n",
    "$A = \\frac{1}{2}|x_2y_3 - x_3y_2 - x_1y_3 + x_3y_1 + x_1y_2 - x_2y_1|.$\n",
    "\n",
    "Write a function *area(vertices)* that returns the area of a triangle whose vertices are specified by the argument vertices, which is a nested list of the vertex coordinates. For example, vertices can be [[0,0], [1,0], [0,2]] if the three corners of the triangle have coordinates (0, 0), (1, 0), and (0, 2).\n",
    "\n",
    "Then, assume that the vertices of the triangle are stored in a dictionary and not a list. The keys in the dictionary correspond to the vertex number (1, 2, or 3) while the values are 2-tuples with the x and y coordinates of the vertex. For example, in a triangle with vertices (0, 0), (1, 0), and (0, 2) the vertices argument becomes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n"
     ]
    }
   ],
   "source": [
    "## STAGE 1: area function on list\n",
    "\n",
    "def area(vertices):\n",
    "    # nb. vertices = [[x1,y1],[x2,y2],[x3,y3]]\n",
    "    x2y3 = vertices[1][0] * vertices[2][1]\n",
    "    x3y2 = vertices[2][0] * vertices[1][1]\n",
    "    x1y3 = vertices[0][0] * vertices[2][1]\n",
    "    x3y1 = vertices[2][0] * vertices[0][1]\n",
    "    x1y2 = vertices[0][0] * vertices[1][1]\n",
    "    x2y1 = vertices[1][0] * vertices[0][1]\n",
    "    return .5*(x2y3 - x3y2 - x1y3 + x3y1 + x1y2 - x2y1)\n",
    "\n",
    "print(area([[0.0,0.0],[1.0,0.0],[0.0,2.0]]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n"
     ]
    }
   ],
   "source": [
    "## STAGE 2: area function on dict\n",
    "\n",
    "def area(vertices):\n",
    "    # nb. vertices = {v1: (x,y)}\n",
    "    x2y3 = vertices[2][0] * vertices[3][1]\n",
    "    x3y2 = vertices[3][0] * vertices[2][1]\n",
    "    x1y3 = vertices[1][0] * vertices[3][1]\n",
    "    x3y1 = vertices[3][0] * vertices[1][1]\n",
    "    x1y2 = vertices[1][0] * vertices[2][1]\n",
    "    x2y1 = vertices[2][0] * vertices[1][1]\n",
    "    return .5*(x2y3 - x3y2 - x1y3 + x3y1 + x1y2 - x2y1)\n",
    "\n",
    "print(area({1: (0,0), 2: (1,0), 3: (0,2)}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Improve a program**</br>\n",
    "The file *data/densities.dat* contains a table of densities of various substances measured in g/cm$^3$. The following program reads the data in this file and produces a dictionary whose keys are the names of substances, and the values are the corresponding densities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def read_densities(filename):\n",
    "    infile = open(filename, 'r')\n",
    "    densities = {}\n",
    "    for line in infile:\n",
    "        words = line.split()\n",
    "        density = float(words[-1])\n",
    "    \n",
    "        if len(words[:-1]) == 2:\n",
    "            substance = words[0] + ' ' + words[1]\n",
    "        else:\n",
    "            substance = words[0]\n",
    "        \n",
    "        densities[substance] = density\n",
    "    \n",
    "    infile.close()\n",
    "    return densities\n",
    "\n",
    "densities = read_densities('data/densities.dat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One problem we face when implementing the program above is that the name of the substance can contain one or two words, and maybe more words in a more comprehensive table. The purpose of this exercise is to use string operations to shorten the code and make it more general. Implement the following two methods in separate functions in the same program, and control that they give the same result.\n",
    "\n",
    "1. Let *substance* consist of all the words but the last, using the join method in string objects to combine the words.\n",
    "2. Observe that all the densities start in the same column file and use substrings to divide line into two parts. (Hint: Remember to strip the first part such that, e.g., the density of ice is obtained as *densities['ice']* and not *densities['ice     ']*.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'air': 0.0012, 'gasoline': 0.67, 'ice': 0.9, 'pure_water': 1.0, 'seawater': 1.025, 'human_body': 1.03, 'limestone': 2.6, 'granite': 2.7, 'iron': 7.8, 'silver': 10.5, 'mercury': 13.6, 'gold': 18.9, 'platinium': 21.4, 'Earth_mean': 5.52, 'Earth_core': 13.0, 'Moon': 3.3, 'Sun_mean': 1.4, 'Sun_core': 160.0, 'proton': 280000000000000.0}\n"
     ]
    }
   ],
   "source": [
    "## STAGE 1: using string function join\n",
    "\n",
    "def read_densities(filename):\n",
    "    infile = open(filename, 'r')\n",
    "    densities = {}\n",
    "    for line in infile:\n",
    "        words = line.split()\n",
    "        density = float(words.pop()) # pop is a list operation that removes the last element from a list and returns it\n",
    "        substance = \"_\".join(words) # join the remaining words with _\n",
    "        densities[substance] = density\n",
    "    infile.close()\n",
    "    return densities\n",
    "\n",
    "densities = read_densities('data/densities.dat')\n",
    "print(densities)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## STAGE 2: Using string slicing and operations\n",
    "\n",
    "def read_densities(filename):\n",
    "    infile = open(filename, 'r')\n",
    "    densities = {}\n",
    "    for line in infile:\n",
    "        density = float(line[12:]) # column 13 onwards\n",
    "        substance = line[:12] # upto coumn 12\n",
    "        substance = substance.strip() # remove trailing spaces\n",
    "        substance = substance.replace(\" \", \"_\") # replace spaces with _\n",
    "        densities[substance] = density\n",
    "    infile.close()\n",
    "    return densities\n",
    "\n",
    "densities = read_densities('data/densities.dat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Write function data to a file**</br>\n",
    "We want to dump $x$ and $f(x)$ values to a file named function_data.dat, where the $x$ values appear in the first column and the $f(x)$ values appear in the second. Choose $n$ equally spaced $x$ values in the interval [-4, 4]. Here, the function $f(x)$ is given by:\n",
    "\n",
    "$f(x) = \\frac{1}{\\sqrt{2\\pi}}\\exp(-0.5x^2)$</br></br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# libs\n",
    "from numpy import *\n",
    "from math import pi\n",
    "\n",
    "# define our function\n",
    "def f(x):\n",
    "    return (1.0/sqrt(2.0*pi))*exp(-.5*x**2.0)\n",
    "\n",
    "# let's make our x\n",
    "xarray = linspace(-4.0, 4.0, 100)\n",
    "fxs = f(xarray)\n",
    "\n",
    "# let's zip them up for a simple for loop when writing out\n",
    "data = zip(xarray, fxs) # this combines each element into a tuple e.g. [(xarray1, fxs1), (xarray2, fxs2) ...]\n",
    "\n",
    "# write out\n",
    "outfile = open(\"ex5_out.dat\", \"w\") # w is for writing!\n",
    "for x,y in data:\n",
    "    outfile.write(\"X = %.2f Y = %.2f\" % (x, y))\n",
    "    outfile.write(\"\\n\")   # ensure newline\n",
    "outfile.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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}