Py_lectures_04_Classes_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python-in-the-lab: introduction to the classes I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import scipy.integrate as integrate\n",
    "import matplotlib.pylab as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will reuse files already analysed, so we can concentrate on the logic!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFg9JREFUeJzt3X+onfV9wPH3J7fX7rYwU9TBcjUkDpc2VLqUYLdKi5V2ibWrwZWi7QrrMlILrmy40Bv6R9mgmBEKrVRq02qdjGqDhDRdZBnMFp11kEgm6lxGpu28N2Om1XSsXshN/O6Pe098cjzPuc+55zznPM953i+4eM8359znk+ea8znf7+f7I1JKSJKaZ9WoA5AkjYYJQJIaygQgSQ1lApCkhjIBSFJDmQAkqaFMAJLUUCYASWooE4AkNZQJQJIa6i2jDqCbSy+9NK1bt27UYUhSbTz11FM/TyldVuS5lU4A69at4+jRo6MOQ5JqIyJ+VvS5DgFJUkOZACSpoUwAktRQJgBJaqhKF4Gr5sCxOfYcPs7J0/OsWT3Fzi0b2LZpetRhSdKKmAAKOnBsjl37n2F+4RwAc6fn2bX/GQCTgKRaMgEUtOfw8fNv/i3zC+fYc/j4UBKAvQ9Jg2YCKOjk6fme2gfJ3oekMlgEzjhwbI5rdz/K+plDXLv7UQ4cmzv/Z2tWT3V8TV77IHXrfUjSSpkAlrQ+Zc+dnifxxqfsVhLYuWUDU5MTF7xmanKCnVs2lB7bKHsfksaXQ0BLlhvjbw21DGscPjvmvyqCcym96TnD6H1IGl9DSwAR8QHg00vX3JhSev+wrl1EkU/Z2UQAbwwZtRLCh955GT/691MdE0R7EXe552bH/Du9+Q+r9yFpfPWVACLiPuBjwMsppXdn2rcCXwcmgO+klHanlB4HHo+IbcCRfq5bhjWrp5jrkATyPmV3Ksz+3b/81/k/zxZqgcLP3bZpumNvBGAigtdTWrb34YwhSUX02wO4H/gG8ECrISImgLuBjwCzwJGIOJhS+relp3wK+NM+rztwO7dsuOBNGrp/ys57k87KFmqLPnfbpunc3sjrKfHi7hu7/hxnDEkqqq8icErpMeCVtuZrgBMppRdSSmeAh4CbACJiLfDLlNL/9nPdMmzbNM2dN1/N9OopAphePcWdN1+d+6ZZtAB78vR8T8+FYjOO8mYsOWNIUlFl1ACmgZcyj2eB9y19vx34brcXR8QOYAfA2rVrSwgv33Jj/NmhlLwho3atN+0iz10VwfqZQ1w8NcnkRLBw7o2x/2xvpNun/CK1DIeIJEE5CSA6tCWAlNKXl3txSmkvsBdg8+bNb65+DslyQymdhozaTa4KXjtzlldfWyBYugldtIq9p+cXmFwVvONtk5x+beFNb9J5n/Lv2Pd07jWyyeVXZ86eTy4OEUnNVcY6gFngiszjy4GTJVynVMsNpXQaMvqj3117/vHqqUkIePW1BWDxzb+VGdufOxFvzpkLryfedtFbeHH3jTwxc/0Fb855n/I7zRbK/lliMblkexbtfy9JzVFGD+AIcFVErAfmgFtYLPzWykqmhWZdu/tRTs8vXNCWWHzzf2Lm+gva188c6imGosNPsJhcuiWGlrnT86yfOeSQkNQgffUAIuJB4ElgQ0TMRsT2lNJZ4HbgMPA8sC+l9Fz/oQ5Xv1s/9LJ6t9drdVqV3EmwOHOoqE4roCWNr35nAd2aUvrNlNJkSunylNK9S+2PpJR+O6X0Wymlrwwm1OHqd+uHXt7Ue71W+/BTpyGk1rVWslq4VU/otCeSpPHhXkA5ep0W2q6XN/WVXGvbpmmemLmeF3ffyFc/+Z7ca3WKo1Vg7pw2FrVqBvYIpPEVqYchgmHbvHlzOnr06KjDWLFhTrfsdq1uf3bt7kcL1RM61S4kVU9EPJVS2lzouSaAZmuf7tpNgEViqeJ6SQDuBtpw7buc5u08ChcWibOvlVRP9gB0gaI9gqIb00kaLnsAWrH2HkHex4NWL8EegVRfzgLSm2RnGE0XmEbqSmKpnkwA6qroojOPp5TqxwSgroouOkvgojGpZqwBaFnZPY+6FYmtB0j1Yg9APcn2CDqxHiDVhwlAPWsVifO2kmjtLOqQkFRtJgCtWLeN5txHSKo+E4BWrMgMIYeEpOqyCKwVK7pozCmiUjWN/VYQHoA+PHk7i7pthDQ8vWwFMdZDQK0pi3NLn04dky5X3pCQZwtI1TTWCWC5g901WEUWjXn/peoY6xpAL+fyajCyi8Z6Pexe0nCNdQ8gb5riqgjnqQ+B91+qtrFOAI5Jj5b3X6q2sU4AjkmPlvdfqraxnwaatX7mUMe56gG8uPvGgV1HneXdf/C8YWlQnAaaI29MutuWBhoct46QqqVRCaDTmPTU5AQ7t2wYUUTN4tYRUrWM9TTQdu1bFzjkMFxFt45o7Sbq70cqV6NqAKqWvK0jsqYmJ7jz5qtNAlJB1gBUCw4JSaPVqCEgVYu7iUqjZQ9AI9U6XezF3TfmHjPpgfNSOUwAqoxuQ0JOEZUGb2gJICKui4jHI+KeiLhuWNdVfXjgvDRcfSWAiLgvIl6OiGfb2rdGxPGIOBERM0vNCfg/4NeA2X6uq/HlgfPS8PTbA7gf2JptiIgJ4G7gBmAjcGtEbAQeTyndAHwR+Ks+r6sx56phqXx9JYCU0mPAK23N1wAnUkovpJTOAA8BN6WUXl/681eBt/ZzXY0/p4hK5StjGug08FLm8Szwvoi4GdgCrAa+kffiiNgB7ABYu3ZtCeGpDpwiKpWvjATQafg2pZT2A/uXe3FKaS+wFxZXAg84NtVI9nSxvFXDrcNl3DZC6l0Zs4BmgSsyjy8HTpZwHTWIh8tIg1dGAjgCXBUR6yPiIuAW4GAJ11GDeLiMNHh9bQYXEQ8C1wGXAv8DfDmldG9EfBT4GjAB3JdS+spKfr6bwSmPh8tInfWyGVxfNYCU0q057Y8Aj/Tzs6Vu1qyeyt1JNDskBJgEpBxuBaFacpqo1D8TgGqpvSaQx5XDUj63g1ZtFZkmCg4JSXnsAWgsOCQk9a7RPYADx+Y8H3hMeN6w1LvGJoADx+bYtf8Z5hfOAQ4RjAOHhKTeNHYIaM/h4+ff/FscIhgfDglJy2tsAsjbRMzNxcaDs4Sk5TV2CChvIVG3fehVLw4JSd01tgfQaYhganKCnVs2jCgilckhIenNGtsDaJ814qyQ8eb5AtKb9bUZXNncDE5l6TYkNO2HAdVYL5vBNXYISM3WbUjIswXUFCYANVJ2llAn1gPUBCYANda2TdM8MXN97jRRp4hq3JkA1Hjdpv563KTGmQlAjecUUTWVCUCN56phNVVj1wFIWa4aVhPZA5DaFB0SumPf0/YIVGv2AKQ2RVcNn1taRGmPQHVlD0DqoDVF9MXdN+auFciySKw6MgFIyygyJATuI6T6MQFIy2ifJTQRnecKrYqwJqBasQYgFZCdJdR+nGiLNQHVjT0AqUdFegTWBFQHbgct9Wn9zKHcmUIBnjWhoXI7aGmI3EtIdWUCkPrkwjHVlUVgqU8uHFNd2QOQBsCFY6qjoSWAiHhXRNwTEQ9HxOeHdV1p2IouHJs7Pe9wkEaqrwQQEfdFxMsR8Wxb+9aIOB4RJyJiBiCl9HxK6Tbgk0ChCrVUR0UXjoEFYo1Wvz2A+4Gt2YaImADuBm4ANgK3RsTGpT/7OPDPwD/1eV2p0rJDQl/95Hu69ggcDtKo9JUAUkqPAa+0NV8DnEgpvZBSOgM8BNy09PyDKaX3A5/u57pSnSx3AD142IxGo4wawDTwUubxLDAdEddFxF0R8S3gkbwXR8SOiDgaEUdPnTpVQnjS8LV6BN2SgGsGNGxlJIBOA54ppfTjlNIXUkqfSyndnffilNLelNLmlNLmyy67rITwpNFxzYCqpIx1ALPAFZnHlwMnS7iOVDuuGVCVlNEDOAJcFRHrI+Ii4BbgYAnXkWrJNQOqin6ngT4IPAlsiIjZiNieUjoL3A4cBp4H9qWUnus/VGn8eNiMRsndQKURO3Bs7vyQ0KqI88M/WRMRvJ6SO4tqWb3sBupeQNKIediMRsUegFQx9gjUD3sAUo1lewTrZw51fI49Ag2Cu4FKFdbtsJkWZwlppUwAUoX1srOoC8fUK4eApAprXziWVxOAC7eSyL5WymMRWKqRvFlC7SwSN5dFYGlMuZWEBskagFQzbiWhQTEBSDXm8ZPqh0NAUo31UiR2OEjtTABSzRXZSqKlNRxkAhA4BCSNFY+fVC9MANKY8fhJFWUCkMaUx09qOdYApDHlmgEtxx6ANMZcM6BuTABSQ3j8pNqZAKSGyM4QChb3C+pkVYQ1gYawBiA1iMdPKsvdQKUG8/jJ8eNuoJIK8fjJZrMGIAnw+MkmMgFIAjx+sokcApIEePxkE5kAJJ3Xy86i4O6idecQkKSO2tcN5HFIqL7sAUjKle0RXLv7UeZyVgk7JFRP9gAkFVJ0d1FnCdWHPQBJhRTdXdS9hOrDHoCkworsLprAekBNmAAkrUi3ISFPG6uHoSWAiLgyIu6NiIeHdU1J5Vnu/GHrAdXXVwKIiPsi4uWIeLatfWtEHI+IExExA5BSeiGltL2f60mqltaQUN40UesB1dZvD+B+YGu2ISImgLuBG4CNwK0RsbHP60iqsLx9hDxboNr6SgAppceAV9qarwFOLH3iPwM8BNzUz3UkVVtePeBcShesETAJVEsZNYBp4KXM41lgOiIuiYh7gE0RsSvvxRGxIyKORsTRU6dOlRCepEErctqYNYHqKWMdQKfhwJRS+gVw23IvTintBfbC4oEwA45NUkmKnC1gTaBayugBzAJXZB5fDpws4TqSKiqvJlDkzAENTxkJ4AhwVUSsj4iLgFuAgyVcR1JFdaoJTK4KXjtz1qJwhfQ7DfRB4ElgQ0TMRsT2lNJZ4HbgMPA8sC+l9Fz/oUqqi/aawOqpSQh49bUFi8IV4qHwkkqXt5Po9Oopnpi5fgQRja9eDoV3KwhJpcsr/loUHi0TgKTSWRSuJhOApNJZFK4mzwOQVLr2swQunprkV2fO8uprC4CniY2KPQBJQ5E9S+Dtb30LC+cunIDiSuHhMwFIGjqLwtVgApA0dBaFq8EEIGnoOhWFpyYn2Lllw4giaiaLwJKGrr0ovGb1FB9652XsOXycv/j+v7Jm9RQ7t2ywIFwyE4CkkcjuHnrg2By79j/D/MI5wFlBw+IQkKSR23P4+Pk3/xZnBZXPBCBp5JwVNBomAEkj56yg0TABSBo5ZwWNhkVgSSPXaVaQs4DKZwKQVAnZWUEaDhOApEo6cGzOHkHJTACSKsd1AcNhEVhS5bguYDhMAJIqx3UBw2ECkFQ5rgsYDhOApMpxXcBwWASWVDmuCxgOE4CkSnJdQPkcApKkhrIHIKnyXBRWDhOApEpzUVh5HAKSVGkuCiuPCUBSpbkorDwmAEmV5qKw8pgAJFWai8LKM7QicERcCXwJuDil9IlhXVdSvbkorDyFEkBE3Ad8DHg5pfTuTPtW4OvABPCdlNLuvJ+RUnoB2B4RD/cXsqSmcVFYOYr2AO4HvgE80GqIiAngbuAjwCxwJCIOspgM7mx7/Z+klF7uO1pJ0sAUSgAppcciYl1b8zXAiaVP9kTEQ8BNKaU7WewtSJIqrJ8i8DTwUubx7FJbRxFxSUTcA2yKiF1dnrcjIo5GxNFTp071EZ4kqZt+isDRoS3lPTml9AvgtuV+aEppL7AXYPPmzbk/T5LUn356ALPAFZnHlwMn+wtHkjQs/SSAI8BVEbE+Ii4CbgEODiYsSVLZCiWAiHgQeBLYEBGzEbE9pXQWuB04DDwP7EspPVdeqJKkQSo6C+jWnPZHgEcGGpEkaSjcCkKSGsoEIEkNZQKQpIYyAUhSQ5kAJKmhTACS1FAmAElqKBOAJDWUCUCSGsoEIEkNZQKQpIYyAUhSQ5kAJKmhTACS1FD9HAkpSSNx4Ngcew4f5+TpedasnmLnlg1s25R7JLlymAAk1cqBY3Ps2v8M8wvnAJg7Pc+u/c8AmAR65BCQpFrZc/j4+Tf/lvmFc9yx72nWzxzi2t2PcuDY3Iiiqxd7AJJq5eTp+Y7t51IC7BH0wh6ApFpZs3pq2efML5xjz+HjQ4im3kwAkmpl55YNTE1OLPu8vJ6C3uAQkKRaaQ3rtGYBrYo4P/yTVaSn0HQmAEm1s23T9PlE0D4rCGBqcoKdWzaMKrzaMAFIqrX2HoHrAoozAUiqvWyPQMVZBJakhjIBSFJDmQAkqaFMAJLUUCYASWqoSB0WUFRFRJwCftbWfDHwyw5P79Te3pZ93On77H8ngZ/3EG5eXEVi6yWuXwKXDii2Qd+zXmMb1D3rFo/3bPnYvGe9x1b2PVsujry4AK5KKV1c6OoppVp9AXuLtre3ZR93+j77X+DoIOIqElsvcS39dyCxDfqe9RrboO7ZMvF4z5aJzXtWvXu2XBxF7lmRrzoOAf2wh/b2th8u8337fwcRV96fLRfLoOLq9rpxuWfd4vGe5T/2nlX3ni0XR15cPcVS6SGgUYqIoymlzaOOoxNj611V44LqxlbVuKC6sVU1rjx17AEMy95RB9CFsfWuqnFBdWOralxQ3diqGldH9gAkqaHsAUhSQ5kAJKmhTACS1FAmgB5ExJURcW9EPNzW/vaIeCoiPlaVuCLiXRFxT0Q8HBGfH0VcXWLbFhHfjogfRMTvVyy2jr/jCsT19oj426X79ulRxZaJZ2NE7IuIb0bEJ0YdT0tErI2IgxFxX0TMjDqerIj4wNK/ye9ExE9GHQ9Qv4VgK/0C7gNeBp5ta98KHAdOADMFf9bDbY//Gvgi8LEqxbXUtgq4t2r3bKntHRWO7U1to4wL+AzwB0vff38lsQ0yRuAO4ANL3x/sJ54Bx/Vh4HNL3z8wiLhK+L1ua8U46q+RBzC0vyh8EHhv9hcITAD/CVwJXAQ8DWwErgb+vu3rNzKvy/7D/DBwC/DHK0wApcS19PjjwE+AT1XpnmXavgq8t6KxrTQBlPX/2S7gd5a+/96o/y0sfd0N7AGeqMq/UeAS4EfAo8BnBxFXCb/XfcCvDyq2vv5eow5gqH9ZWNf2C/w94HDm8S5gV4Gfk/2H+RXga8A/Aj8AVlUhrrb2QxW7ZwH8DfDhqv0+l7uXI7xnn2HpAwbwUD/3bcAxTgA/6DeeQcUF/CXwwX5/h2XdM2At8O1BxtXPV9NrANPAS5nHs0ttHUXEJRFxD7ApInYBpJS+lFL6c+B7LP5iX69CXBFxXUTcFRHfAh4ZQEwDiw34MxZ7Tp+IiNuqFFtOvCOPC9gP/GFEfJOVbzswyBjXRcRe4AEWewFl6Sku4B+ALyzdv5+WGBf0HhvAduC7pUXUo6afCRwd2nJXxqWUfgF0fMNKKd0/oJhgAHGllH4M/HiAMbUMIra7gLsGHBcMJrbc3/GI4/oV8NkBx5XVa4w/BXaUFs0beo3rWWBYRemeYgNIKX25pFhWpOk9gFngiszjy4GTI4olq6pxgbGtRFXjyqpqjFWNC6odWyFNTwBHgKsiYn1EXMRiMffgiGOC6sYFxrYSVY0rq6oxVjUuqHZsxYy6CDGsL+BB4L+BBRYz9/al9o8C/8FiNf9LxmVs4xhXHWKsalxVj62fLzeDk6SGavoQkCQ1lglAkhrKBCBJDWUCkKSGMgFIUkOZACSpoUwAktRQJgBJaigTgCQ11P8DlTovgi9C23IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mainDir = \"/home/gf/src/Python/Python-in-the-lab/Bk\"\n",
    "filename = \"F64ac_0.01_S.dat\"\n",
    "filename = os.path.join(mainDir, filename)\n",
    "s, ps = np.loadtxt(filename, comments=\"#\", unpack=True)\n",
    "plt.loglog(s, ps, 'o'); # Nothing new..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We would like to work a little with the data\n",
    "\n",
    "The idea is the following:\n",
    "\n",
    "Can we construction something general and sufficiently abstract from which we can extract the values more easily?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your first class\n",
    "class Sdist: # Note the capital letter\n",
    "    def __init__(self, filename): # Self???\n",
    "        self.size, self.ps = np.loadtxt(filename, comments=\"#\", unpack=True)\n",
    "\n",
    "# Done?\n",
    "THIS_IS_A_CONSTANT = 131231.23131   \n",
    "filename = \"F64ac_0.01_S.dat\"\n",
    "filename = os.path.join(mainDir, filename)\n",
    "s001 = Sdist(filename)\n",
    "# What is \"self\" ???"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*self* is not a reserved word, but it is universality used a (self)-reference to the class.\n",
    "It simply indentifies the class itself.\n",
    "\n",
    "s001 is a call to the class, named a *instance* of the class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([7.09426e+05, 4.26445e+05, 2.41955e+04, 1.50198e+07, 7.44998e+06,\n",
       "       2.29251e+07, 1.12218e+07, 1.05320e+07, 6.66718e+06, 1.16563e+07,\n",
       "       1.09002e+07, 6.96781e+06, 1.24547e+07, 1.95546e+07, 7.87439e+06,\n",
       "       6.28456e+06, 1.00967e+07, 6.11875e+06, 7.34575e+06, 5.14724e+06,\n",
       "       4.00538e+06, 2.85764e+06, 2.30620e+06, 1.50140e+06, 1.02132e+06,\n",
       "       6.97577e+05, 4.54409e+05, 3.07353e+05, 2.12505e+05, 1.47726e+05,\n",
       "       9.67887e+04, 6.87491e+04, 4.69031e+04, 3.30666e+04, 2.30293e+04,\n",
       "       1.62495e+04, 1.15204e+04, 7.98031e+03, 5.68260e+03, 4.23867e+03,\n",
       "       2.97384e+03, 2.13527e+03, 1.55214e+03, 1.17525e+03, 8.23620e+02,\n",
       "       6.61500e+02, 4.52724e+02, 3.50047e+02, 2.60123e+02, 1.88037e+02,\n",
       "       1.34791e+02, 9.54933e+01, 7.68208e+01, 5.41981e+01, 4.18298e+01,\n",
       "       2.87398e+01, 2.15285e+01, 1.39254e+01, 8.90375e+00, 5.06903e+00,\n",
       "       2.97192e+00, 1.50779e+00, 5.68595e-01, 2.21021e-01, 7.63318e-03,\n",
       "       6.06325e-03])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s001.ps[s001.ps != 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "s001.<tab> ? \n",
    "s001.__???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's redifine the class\n",
    "# as I want to get rid of the zeros of ps\n",
    "class Sdist:\n",
    "    def __init__(self, filename):\n",
    "        self.size, self.ps = np.loadtxt(filename, comments=\"#\", unpack=True)\n",
    "        s_len_ori = len(self.size)\n",
    "        self.size, self.ps = self.avoid_zeros()\n",
    "        print(\"%i lines deleted\" % (s_len_ori - len(self.size)))\n",
    "    \n",
    "    def avoid_zeros(self):\n",
    "        is_not_zero = self.ps != 0\n",
    "        s = self.size[is_not_zero]\n",
    "        ps = self.ps[is_not_zero]\n",
    "        return s, ps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14 lines deleted\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFqpJREFUeJzt3V2MHed52PH/w9XaWsuAaFi+qFYiSBcqHUa+ICrIQIkUshCXdKJEguI2UoJcNCpYB3WKAi6RJVxAaYFAC/CmNixEZSyGdYFIVWWCYUu1vGEMGYIvSJU1ZNlhoCoJtMuiYmyvi0gLcEk9vdhd+uzhzO6cnTPna/4/YME9s3POPPvBec77Pu9HZCaSpPbZMewAJEnDYQKQpJYyAUhSS5kAJKmlTACS1FImAElqKROAJLWUCUCSWsoEIEktZQKQpJa6bdgBbOauu+7K3bt3DzsMSRobr7/++t9k5ieqnDvSCWD37t1cvHhx2GFI0tiIiL+ueq5dQJLUUiYASWopE4AktZQJQJJaaqSLwKPm9KVFjp27zJWlZe7eOcORg3t5bP/ssMOSpG0xAVR0+tIiR0+9wfLKDQAWl5Y5euoNAJOApLFkAqjo2LnLN2/+65ZXbnDs3OWhJwBbJpK2wwRQ0ZWl5Z6OD4otE0nbZRF4zelLixyYP8+eubMcmD/P6UuLG75+986ZwueVHR+UzVomkrQZEwA/exe9uLRM8rN30Z1J4MjBvcxMT2143sz0FEcO7h1wtBuNastE0uizC4hq/fvr/w67r727v3/nR6b5yfsrt5w37JaJpNE3sAQQEb8A/ObaNfdl5j8Y1LW3UvVd9GP7Z2+54RcVYGHrRLGd5xX190/vCKangpUbefO8UWiZSBp9tRJARJwAHgHezcz7O44fAr4KTAHfyMz5zPwO8J2IeAy4UOe6/Xb3zhkWC5LAVu+ii27IR/7L9yC4eUMuKspu93lFLZWVD5KdM9Pc8eHbem6ZOHpIare6NYCTwKHOAxExBTwLfB7YBzwZEfs6TvkN4IWa1+2r7fbvl92QO9+Nw61F2e0+r6yl8tPlFV6be5i/nP9lXpt7uPLNf6u6h6TJVisBZOarwI+7Dj8IvJWZb2fmNeBF4FGAiNgF/DQz/1+d6/bbY/tneebxTzO7c4YAZnfO8Mzjn97yRtpLobXz3O0+r+pIpK1GNIGjhyQ1UwOYBd7peLwAfGbt86eAP97syRFxGDgMsGvXrgbCK1a1f7/znLKuoyKdN+lenrcjgj1zZ7l75wyf/dQn+Nbrixtu3N0tlarzAqrWPewmkiZXZObWZ232AhG7gf+2XgOIiH8MHMzMf7b2+LeABzPzd3t97QceeCCHtSFM940UVm+2nS2DonOmd8SGvvz1Yx+9/TaW3l8pvZEXPa/bzPQUv/b3Z/mzP79aekM+MH++MLl01wnev3a9cPRQ53l3zkzz3rXrtxSYq7SOJA1HRLyemQ9UObeJFsACcG/H43uAKw1cp1F1hoZ2Hlu/ia7fbBeXlvnW64uFN/LO5+2I4EbeWhP4sz+/ymtzD5fGXfbOfml5haXln8VQNHpoekfw3rXrN89b/3ezn4Gk8dVEArgA3BcRe4BF4AlWC79jpc7Q0PXjsPqOvPtGutmNfP15e+bO9hTXuqrdS0Wjh8paBd0Wl5Y5MH/ebiFpzNUqAkfEC8B3gb0RsRART2XmdeBLwDngh8BLmflm/VAHq19LP2x3pu52r180oqlM9+ihpQo3f4AARw9JE6DuKKAnM/PvZOZ0Zt6Tmc+vHX8lM/9eZv7dzPyD/oQ6WP1a+qGfN/Iq1y8a0fSxj0xXiqFKcgugu0qxvHKD3z/z5pYjjySNFtcCKrHdoaHd+nkjr3r9x/bPbnhn//Sv/HylGIpind4RfOwj0zdjKCtRLy2v2CqQxkztUUBNGuYooH4ahaGUVWPY6ryyUUZFZnfObFqwltR/vYwCMgGoJ0VDXzcTYKFYGqBhDwPVBCsa+rrZ6KHOLqHO50saPlsAqq1qq2C7i9ZJqs4WgAaqu1WwWaG4czKarQJpuBwFpL7oHHk0W3GuhIvPScNlAlDf9TIZza0rpeGxC0h910uh+M6ZaZeVkIbEBKBGdK+RVLZyaufic9YFpMGyC0gDUTSz+aO337blLmiSmuMwUA3NnrmzpSOGnEAmbU8vw0BtAWhoNlt8zjWFpOaZADQ0VUYL2SUkNccisIam6gQyh4pKzZi4GsAorLyp7am6n7G/U6lca2sA60MNXZd+PJXtR/Detev+TqUGTFQC2Gwjd40+h4pKgzVRNYDt7r+r0dE9gWzP3NnC8/ydSvVNVAK4e+dMYR+yyw2ML3+nUnMmqgvIPuTJ4+9Uas5EJQD7kCePv1OpORM3DLRb2XIDAfzl/C/Xem0Nh0tISOVaOwy0SNlyA5stQ6DR5hISUn9MfAIo6kOemZ7iyMG9Q4pIdbmEhNQfEzUKqEjR5iR2D4y3qktILC4tO1JI2sTE1wA0+cqWkAjYkBxmpqd45vFPmwQ00awBqFWKuoS6b/5gt5DUzQSgsVc0VNSVRaWt2QWkieTKomoru4DUes4glrY2sAQQEQ9FxHci4rmIeGhQ11U7OYNY2lqtYaARcQJ4BHg3M+/vOH4I+CowBXwjM+dZrcn9LXA7sFDnulIVVVcWXVxaZs/cWbuE1Dp1WwAngUOdByJiCngW+DywD3gyIvYB38nMzwO/B/zbmteVeuYMYmmjWgkgM18Fftx1+EHgrcx8OzOvAS8Cj2bmB2tf/wnw4TrXlbbDGcTSRk3MBJ4F3ul4vAB8JiIeBw4CO4Gvlz05Ig4DhwF27drVQHhqKzehlzZqIgFEwbHMzFPAqa2enJnHgeOwOgy0z7Gp5TrrAmVDRd1sRm3RxCigBeDejsf3AFcauI5Ui0NF1XZNJIALwH0RsSciPgQ8AZxp4DpSLQ4VVdvVmgkcES8ADwF3Af8XeDozn4+IXwL+PavDQE9k5h9s5/WdCaxB22yzmdmdM3YLaeT1MhO4Vg0gM58sOf4K8Eqd15aGoWwT+oCbx9e7hQCTgMaaS0FIHVxZVG1iApA69LKy6PoM4gPz5y0SayxN/I5gUq+6l5AoGy4KG2cQrz9XGhe2AKQtOINYk6qVLYDTlxbdI1iVuQexJlXrNoQ5fWmRo6feYHnlxs1j7hWrXrgHsUaZG8Js4ti5yxtu/mDzXb1xpJAmReu6gMoW+nIBMFXV3SVUNncA7BbSaGtdAij7z7rZWvFSt6ojhZxAplHWui6goub7zPQURw7uHVJEmgR2C2kcta4FUNR8t1muunrpFrK7UaOidaOApEEp6xaaiuCDTN98qBGOApJGQNkEshuZ7jWgkWACkBrSva7QVNy6WZ41AQ2TXUDSgLjXgAbBLiBpBJUNNV4fKmq3kAbNBCANiENFNWpaNwxUGhZnEGvUmACkAXIGsUaJRWBpiIpWpy3qFgLYOTPNHR++zVaBNmURWBoTvWxBubS8YrFYfWUXkDRkvWxB2Wm9WGwrQNtlC0AaMVW2oFznukKqwwQgjZiibqGPfWS68NwdEeyZO8uB+fN2B6lndgFJI6i7W6ioWAyr6wqBI4W0PbYApDHgukJqgsNApTHkukIq4zBQacK5rpD6wRaANIacQKYytgCkCecEMvWDo4CkMeUEMtU1sBZARPxcRDwXES9HxO8M6rpSW/QygWxxadn5A6qXACLiRES8GxHf7zp+KCIuR8RbETEHkJk/zMwvAv8EqNQ/Jam6XiaQAXYJqV4ROCL+IfC3wDcz8/61Y1PAXwCfAxaAC8CTmfmDiPhVYA74emb+yVavbxFYqqdsAlm32Z0zvDb38ICiUpN6KQLXqgFk5qsRsbvr8IPAW5n59lowLwKPAj/IzDPAmYg4C2yZACTV070JTdnbPTegaacmisCzwDsdjxeAz0TEQ8DjwIeBV8qeHBGHgcMAu3btaiA8qV06i8VuQKNOTSSAW+eoQ2bmt4Fvb/XkzDwOHIfVLqC+Ria13JGDeyvNH1heucHvn3lzw/aVtgomTxOjgBaAezse3wNcaeA6knrk/AF1aqIFcAG4LyL2AIvAE8BvNHAdSdvg/AGtqzsM9AXgu8DeiFiIiKcy8zrwJeAc8EPgpcx8s36okprgBjTt5VpAkjh9aXFDf//7167zk/dXbjlvKoIPMq0JjLCBDQOVNBncgKadbAFIKtTZKtgRcfPm38mVRkePLQBJtXW2CvbMnS08Z2l5haXl1a4iWwXjx+WgJW2pbAOabm5LOV5MAJK21OtKowfmz7va6BiwC0jSlrrXFNpspJDLSowPi8CStqWXbSkdPjo4bgkpqXG9LCtxI9MlJUaQXUCStm07y0q4pMTosAUgqW+qFostFI8GWwCS+qa7WFw2gcxC8WgwAUjqq85uoaqFYruFhsMEIKkxRcNHy2oEbks5eCYASY2qWii2W2jwLAJLGqiiQnFZt9CXX/qeheIGmQAkDZTzB0aHXUCSBs75A6PBFoCkoas6f8AtKfvLFoCkoas6f+DOmWlHCvWRCUDSSNhq/sD0juC9a9fdgKaPXA1U0kiqulG921Ju5JaQksZed6HYbSn7zyKwpLHgtpT9ZwKQNBZ63ZbSCWRbswtI0ljoZVtKYMMEss7n62dMAJLGRnddoGi0UDcnkJUzAUgaW92tgrIxja40WswEIGmsdbYKXGm0NxaBJU2MXlYadaSQCUDSBOllpVHXFXImsKQJV9YtNKkziHuZCWwLQNJEK+oWWl9XaHGtcNzW/QYGlgAi4pMR8XxEvDyoa0pSUbfQR2+/jZUbG3s/2lgXqDUKKCJOAI8A72bm/R3HDwFfBaaAb2TmfGa+DTxlApA0aFXXFWpbXaBuC+AkcKjzQERMAc8Cnwf2AU9GxL6a15GkvilbV2hHRKuWkKiVADLzVeDHXYcfBN7KzLcz8xrwIvBonetIUj+VrSvUtj2Im6gBzALvdDxeAGYj4uMR8RywPyKOlj05Ig5HxMWIuHj16tUGwpPUdt11gamIW85pQ02giZnAt/4kITPzR8AXt3pyZh4HjsPqMNA+xyZJwMa6QFtrAk20ABaAezse3wNcaeA6ktQXZTWBqnsQjKsmEsAF4L6I2BMRHwKeAM40cB1J6ouimsDM9BSf/dQnODB/fmILw3WHgb4APATcFRELwNOZ+XxEfAk4x+ow0BOZ+WbtSCWpIUV7DXz2U5/gW68v3lxqehIXkXMpCEkqULaExOzOGV6be3gIEVXjUhCSVFNZAXiSCsMmAEkq0IbCsAlAkgq0oTDsjmCSVKANhWETgCSV6F5E7sD8+Vs2oB/nTeftApKkiiatMGwCkKSKJq0wbAKQpIrKCsNHDu4dUkT1WAOQpIqKCsPrN/8D8+fHbn9hE4Ak9aC7MHz60iJHT70xliOD7AKSpBqOnbtcOjJo1JkAJKmGcR4ZZAKQpBrGeWSQCUCSahjnkUEWgSWphrKRQaNeAAYTgCTV1j0yaFyYACSpz05fWhyLFoEJQJL6aJzmBVgElqQ+Gqd5ASYASeqjcZoXYAKQpD4ap3kBJgBJ6qNxmhdgEViS+mic5gWYACSpz8ZlXoBdQJLUUrYAJGkARnFymAlAkho2qpPD7AKSpIaN6uQwE4AkNWxUJ4eZACSpYaM6OcwEIEkNG9XJYQMrAkfEJ4GvAHdm5hcGdV1JGrZRnRxWKQFExAngEeDdzLy/4/gh4KvAFPCNzJwve43MfBt4KiJerheyJI2fUZwcVrUFcBL4OvDN9QMRMQU8C3wOWAAuRMQZVpPBM13P/+3MfLd2tJKkvqmUADLz1YjY3XX4QeCttXf2RMSLwKOZ+QyrrQVJ0girUwSeBd7peLywdqxQRHw8Ip4D9kfE0U3OOxwRFyPi4tWrV2uEJ0naTJ0icBQcy7KTM/NHwBe3etHMPA4cB3jggQdKX0+SVE+dFsACcG/H43uAK/XCkSQNSp0EcAG4LyL2RMSHgCeAM/0JS5LUtEoJICJeAL4L7I2IhYh4KjOvA18CzgE/BF7KzDebC1WS1E9VRwE9WXL8FeCVvkYkSRoIl4KQpJYyAUhSS5kAJKmlTACS1FImAElqKROAJLWUCUCSWsoEIEktZQKQpJYyAUhSS5kAJKmlTACS1FImAElqKROAJLVUnS0hJUk1nL60yLFzl7mytMzdO2c4cnAvj+0v3Vq970wAkjQEpy8tcvTUGyyv3ABgcWmZo6feABhYErALSJKG4Ni5yzdv/uuWV27w5Ze+x565sxyYP8/pS4uNxmALQJKG4MrScuHxG5nAYFoEtgAkaQju3jmz5TnLKzc4du5yYzGYACRpCI4c3MvM9NSW55W1FPrBLiBJGoL1bp31UUA7Im52/3Sq0lLYLhOAJA3JY/tnbyaC7lFBADPTUxw5uLex65sAJGkEdLcIBjEvwAQgSSOis0UwCBaBJamlTACS1FImAElqKROAJLWUCUCSWiqyYOLBqIiIq8BfF3zpTuCnFY5v9rjo8/V/7wL+Zhshl8VV5Zxxjb3Ksc3i7TzWVOxVf+ZlsXZ/bdxi3+rzYcc+6L/1zWLb6uvDjr3K/9P7MvPOSq+WmWP3ARyvcnyzx0Wfd/x7sZ9xTXLsVY5tFu8gYq/6M9/s5zzOsVf42xlq7IP+Wx/n2OvcY4o+xrUL6L9WPL7Z46LPy163qirPn7TYqxzbKt6mY6/6M+8+VvZ9jFvsVT7fjn7FPui/9SqvMaqx17nH3GKku4CGJSIuZuYDw45jO4x9OIx98MY1bhid2Me1BdC048MOoAZjHw5jH7xxjRtGJHZbAJLUUrYAJKmlTACS1FImAElqKRNARRHxyYh4PiJe7jp+R0S8HhGPDCu2rRTFHhE/FxHPRcTLEfE7w4xvMyWxPxYRfxQRfxoR/2iY8ZUpibvwb2jUlMR+R0T8x7Wf+28OM74qImJfRLwUEX8YEV8Ydjy9iIhdEXEmIk5ExFyjF9vuRIpx+gBOAO8C3+86fgi4DLwFzFV8rZe7Hv874PeAR8Yt9rVjO4DnxzT2jzUR+wDivuXYqMcO/BbwK2uf/+em4u/X9wB8GfiFtc/PNBlvA7H/IvDP1z7/ZpPxtmVDmJPA14Fvrh+IiCngWeBzwAJwISLOAFPAM13P/+3MfLf7RSPiF4EfALc3EzbQUOxrr/OrwNza6zfhJA3FvubfrL1Wv52k2bibdJJmYr8HeGPt8xsFX++nk9T8HoD/BDy99jf+8Ybj7XSS+rFfAr4SEb/O6vfRmFYkgMx8NSJ2dx1+EHgrM98GiIgXgUcz8xmganfOZ4E7gH3AckS8kpkf9CfqVQ3GTmaeAc5ExFngT/oT8YbXbyT2iAhgHvjvmfk/+xfxqiZ/5k1rMPYFVpPA/6LhruM+fg//Yu3me6qpWLv1I/aI+NfA02uv9TLwx03F2+YawCzwTsfjhbVjhSLi4xHxHLA/Io4CZOZXMvNfsXrz/KN+3/w3UTv2iHgoIr4WEf8BeKXRaDeqHTvwu6w2k78QEV9sLNKN+vEzL/peBqEfP/NTwK9FxB/Sn6UYetXr97A7Io6z+k78WMOxbaWn2IH/AfzLtd/BXzUYVztaACWi4FjprLjM/BFQeLPJzJN9iqmq2rFn5reBb/c1qmr6EfvXgK/1Oa6t9CPu0r+hhvUj9veAf9rnuHrR6/fwV8DhxqLpTa+xfx8YSOG6zS2ABeDejsf3AFeGFEuvjH3wxjVuGO/Y143z9zCysbc5AVwA7ouIPRHxIeAJ4MyQY6rK2AdvXOOG8Y593Th/D6Mb+6CGRw3zA3gB+D/ACqvZ+Km1478E/AXwv4GvDDtOYx+Nj3GNe9xjn4TvYdxidzE4SWqpNncBSVKrmQAkqaVMAJLUUiYASWopE4AktZQJQJJaygQgSS1lApCkljIBSFJL/X/E1eSRXAW9VQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "s001 = Sdist(filename)\n",
    "plt.loglog(s001.size, s001.ps, 'o');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.15285e+01, 1.39254e+01, 8.90375e+00, 5.06903e+00, 2.97192e+00,\n",
       "       1.50779e+00, 5.68595e-01, 2.21021e-01, 7.63318e-03, 6.06325e-03])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s001.ps[-10:] # good!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Improvements\n",
    "* Why don't we give just the frequency, so it takes care of all the rest?\n",
    "* Can we leave the possibility to avoid the zeros as a choice?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's redifine the class again\n",
    "class Sdist:\n",
    "    def __init__(self, freq, mainDir, is_avoid_zeros=True):\n",
    "        filename = \"F64ac_%s_S.dat\" % freq\n",
    "        filename = os.path.join(mainDir, filename)\n",
    "        self.size, self.ps = np.loadtxt(filename, comments=\"#\", unpack=True)\n",
    "        s_len_ori = len(self.size)\n",
    "        self._filename = filename\n",
    "        if is_avoid_zeros:\n",
    "            self.size, self.ps = self.avoid_zeros()\n",
    "            print(\"%i lines deleted\" % (s_len_ori - len(self.size)))\n",
    "    \n",
    "    def avoid_zeros(self):\n",
    "        is_not_zero = self.ps != 0\n",
    "        s = self.size[is_not_zero]\n",
    "        ps = self.ps[is_not_zero]\n",
    "        return s, ps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f8ce69486d8>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt0nPV95/H3dzSSJZmDBLZAN1PbCcdbStSYOLSbGJoNiU0iLoZSl0vTbJfUpFlApLt27JRjFJaAg3dJTAMh1KWBbRqi9XGMQUntFErAkOzaxqmBpk6Jk9aypFqykaktCV3mt3+MRp4ZPSPNaG7PzHxe5+hY88wzM1+N5Oc7v9v3Z845RESk9ATyHYCIiOSHEoCISIlSAhARKVFKACIiJUoJQESkRCkBiIiUKCUAEZESpQQgIlKilABEREqUEoCISIkK5juA6cyfP98tXLgw32GIiBSM/fv39zvn6pI519cJYOHChezbty/fYYiIFAwz+5dkz1UXkIhIiVICEBEpUUoAIiIlSglARKRE+XoQ2K92HDjK5l2H6B4YorG2irUrl7BqaVO+wxIRSYkSQIp2HDjKhu2vMzQ6DsDRgSE2bH8dQElARAqKEkCKNu86NHnxjxgaHWfzrkM5SQBqfYhIpigBpKh7YCil45mk1oeIZFLRDwJ3Hu5kxbYVtDzZwoptK+g83JnU43YcOMqHN73AovWdfHjTC+w4cBSAxtoqz/MTHc+k6VofIiKpKuoE0Hm4k/ZX2+k53YPD0XO6h/ZX22dMApFP2kcHhnCc+aS948BR1q5cQlV5Wcz5VeVlrF25JIs/SVg+Wx8iUnyKrguo83AnW17bQu/pXsyMkAvF3D88PsyW17bQurg14XNM90n7lfUfnTxnNv3wqfbhR58fMGPcuSnn5KL1ISLFJ2cJwMwuA26ZeM2LnHMfyvRrRD7xD48PA+A8LpYAvad7p32emT5pr1raFHPRjnQXRS7q/+k/1PH3/9Q35SLv1Ye/9v/8A1969k0GBkenJIT4870u/rlqfYhI8UkrAZjZE8BVwDHn3MVRx68EtgBlwFbn3Cbn3MvAy2a2CtibzusmsuW1LZMX/+nUz62f9v7G2iqOeiQBr0/aXhf1v/7Jv07e/4F3fsgHd/wx7pnj/Dbz+fj477GT5ZP3j4Ycbw+OTj42elDXqyUCUGZGyLkZWxCaMSQi00m3BfAt4OvAU5EDZlYGPAJ8HOgC9prZTufcP06ccjPwmTRf19NMn+wBKssqabukbdpz1q5cEnNRh8SftBNdpAGuCexhU/lWqm0EgHr62FS+FUZhZ2i55/nrrIPGZ47Di83cdmoJH634KY3WT7ebz4Njq9kZWk7IOX65KXEXFmjGkIjMLK1BYOfcS8CJuMOXAm855w4750aAp4FrAczsAuCkc+6ddF43kUSf7AMWwDAa5jbQ/qH2afv/AVaVvcL+s+7icOUt7Km4k5sqf0JleYDPf/enMTOCYPoB2HXBjsmLf0S1jbAu2DHl3EiyaA70E8DBySN8Kvh34dsGzYF+NpVv5ZrAnpiWSKLZSsnMGEr0WBEpDdkYA2gCjkTd7gJ+a+L7W4G/mu7BZrYGWANwwQUXpPTCbZe0xYwBQPgTfzIX/UkHO+DZO6keDV/YmwP9bHSPcXp4jJ0sn/JJOlF3UfDsA/yX8yroDS6gfmyctrcHaD09CECjHZ9yvleysLhzqm2Eh8ofo2z4G/DVZva+5w427P21yQt9dHfTd0PzeDCwekpLI5Kw1EIQEUs0UJr0E5gtBJ6LjAGY2e8BK51zn5m4/SngUufcHak+97Jly1yqG8JEzwKqn1tP2yVtyV/8Ab56MZw8MuVwV2g+y0cenrzdVFvFK+s/OnkhHa3ax5y6XVj5AG68CguMYIEzn8ArQyHa+0/Qenpw8rlqq8o5PTLG6Ljj8JybCcRf8WcwxBy+MHIrO0PLp3Q3AQy6CtaPfiYmCUTGDxLNKIr8XCJSmMxsv3NuWTLnZqMF0AUsiLrdDHRn4XU8tS5uTe2CH+9k1+S3nXOr2XJOLb3BMurHxgkeO8DYO0uB2BlB//D2C2z7l+0QCA/mWnBqi2A4EGDLObV87NS7VNswhytvIXB2+FP8Xf94Id2D82m2/pRCreJdHip/jK/xKCECBC12ymuku2nnyJkEELnoe138o38uESl+2VgIthe40MwWmVkFcCOwMwuvkx01zUD44t8+/1x6yoM4M3rKg1Q2bCd49gEgakbQwQ5e+dVXJi/+0+kNlmEY59qpyX7+D75+D698sp/mGx5g0FXEnB+5RnfOrWZFcyMtCxewormRzrnVk+cELUTAmHLxj2iyfg7PuZlX5tzJNYE9k8evCexhT8WdHJ5zM3sqztwXMNOYgEiJSCsBmNl3gB8DS8ysy8xudc6NAbcDu4CfAR3OuTfTDzVHrtgI5VVsOaeW4UDs22OBUebU7TozI2hivKA3yXexbsxRYWOxB0eH4Pl7oWU1D5Z/jq7QfELO6ArN56nxj/FU9flTElH7/HNjksB0zCBg4UQQGUSOGXCOG2Aed27K6mcRKU5pjwFk02zGADLiYActr92LM49OeQf/4zf/NjxQOjFesKK5kZ7y6XvTXKicB/p6uHpw0ONeo/MPv80DP3mIgZFjuNFa3u1bydg7SznrvZuw8oEpjwhMXKjjB5hn4hyMe3QXAYy5AAFCMVNOk11zICL+kMoYQFHXApq1ltXUn9XoeVfDWQ2TF8HOsRPhi3+w7Ex/zYRgKETt+DjmHPVjjqqTN7L0lPen9s66Ztpfbefk6LHwJ/aKASobtlNX/yZWftLzMSGzyRbB3fPP5bILmia7iJ6tnhsTTnQX0soFjew6q9LzOSPdSWoRiJQGJYAE2i5po7Is9kJZGQrR9m/dcLAjXHaiLtw1g1n4yzlwjobRMe7rP8HL/3qUg0f7+eEHN7L3rnU03/AAlMetJg6Us6XSTVnBbIFRzmn+OxpmWLUMMBYIMFBWNpkQNtTV8zdnhR/nNZaRTBdSZMpp9BiBKo+KFBclgARaF7fS/qF2GsprsImL+rX/footc8Zpee1evvjyBobju4jMaBgfZ/fxYVpPD0HNArj6YWhZHb6/ZXX4ds0CwKDqXDCjt8x7/mfv6V7PRDQTC4zyZFNdwrGMyIykmXi1CI4ODGmAWKRIFF010HR4rSHY/W8DcPLI5CfpyMU03DEyVW8wCF/4ZeIXaVl9JiF89WIYOkH92LjnGEL93PrJKa3TVTj10jN6kpbmusQF8YJlYGXgQmABcN7lLCIiLYKv8Sjdg/P52vduBD6nMQGRAqYWwISEeweMhStdeH2S9lI/tyH5F51Yc9D29gCVodiLenTNotbFrey+YTcHP32Q+5ffn3SLwEG4a8orzvEQXPcYtA+E/43vmvIQ3SL4sn2Dy3dcSuieGnrb38vend9MKiYR8Q+1AAAOdrBl370Mx3XFDI8Ps2XeubSePh3+xDyDZArNxahphpNHJmfwTC46C0HbZd7lK+JbBGdXnM3g2CCjoZnXIUQEnWNozlm0HLiPs998GDPjZHMd9eOOy0+d4qW5VfSWlU07w2iOjTOHU0C4yF3N/rvZC3zwmtuS//lFJK80DXRiLn9L03zPaZ8GHOzqY8X553h200xOxwxB20e+ktoq5InXZjRq9W15Vey4QRKiu64SdU2FfxZLOWFEl7CYSS911Le/lXTcIpJ5mgaaiufvhdEh6se8+8Dr5zbA1Q/T9m7Z1G6aUIj7+45z8FdH2P2vXSld/DsPd7Li51tpaa5jxQUL6Jw7d+qgcZKiu4gaEnRBNcxt4OCnD1JdXp1Sa2E4EGDLuTMPGAOc7/ro2vge2u+7R4PEIgVACSCZfviW1bT+1zdo/50HaRh3k7OCYj4ZT5SQSEbseAP0lBntDU10XvuVyYv/bDez95y+GtU1lcyeCfF6gkHPMhTxbGJ8YN3oo+z53qNKAiI+pwQwceFuPT1Ie/8JGkbHwhf4cTeljHTr4lZ2L9vIwaP97O7qPnPxL68Kl5BIktfOZZG9imH2m9lHYmz/UDsNcxs890CYaTe0RKLXEOw4ez4jLvHwUbWNcBdPa82AiM9pDGA2/fAHO8JdRye7wgnkio0pddu0PNni2VdvGAc/fZAV21bQc7pnyv0NcxvYfcPupF/HS/y+ybMRsAAhF+K80RCfP3GCqzzKWzgXnoUUqF2Q8vsjIrOX73LQhSVyYUrlgh49lz9J0QO1ZuY5Pz/y6TxRN81sum/iec0iMjNOvnuS+rn1XN58OS91vTTtgHJkHcKx8gBfOm8e1seUQWKziQ1tTh5haPvtvPGrtzVDSMRnlABgVhf0VMR/6va6+Ef309fPrfdsAcy2+yZesnsmJGqJRBs2Y8u550w7S6iKd2nc/yA7FlylhWMiPqIxgBzw6vOHxHsVzzSQmyvJlqHoCZZNDhI/V+09SNzAcY0JiPhMybQAdhw4yuZdh+geGOLTZ/0/1pV/l+qh3ln14acqUdeNc46Dnz445Xh8N82strbMgFTKUEQGiRN1CYUwXh66jq6N89la8Qe8v3WNWgMieVYSg8DRG6B77Z07m8VXqcjmoG4uJT2A7BwNUauInYutSDHoKtjo1rD8OtUSEsk0LQSLs3nXIYZGwwu91gU7Yi/+cGZXrizxS5dOuuKnmCYUNWV0Z/XcKeWINE1UxB9KIgFEb3TemGjj9ajN4DNtprn5hSSZVccRw4EAXz+3xvO+Juvn5aHrwhVRD3ZkI1QRmUFJjAE01lZxdCIJdLv5NHskgcGqepLbZXd2kp15U0jaLmmbsUsoMkAcX1guepro2DN3hP8QtVZAJKdKogWwduUSqsrD1TwfHFvNoKuIuX/QVbDx9O+qdEGKols2CUVtXZloJ7Lg+DCDP0h+JbWIZEZJJIBVS5t44Pr30VRbxbOh5Xxx7I/pCs3n2eq5XNHczG8tqucHC17iyz/6dr5DLTiRLqFNl22accrodDuRVQ6lv8hNRFJTEl1AEE4CkRkni9bDc3OrqazbjgXClTEDFQMMBZ+m8/BvFF1XTS7ETxlNvGOa974KIWdwTw3HrI4jl6zVqmGRHCiJFkC8xtoq5tTtmrz4R1hgdLIgm6QuuQHiwJTKos6d2W2snj4u3n+3dhgTyYGSTABrVy7Bygc878tEvR1JvIrYmYsZE/CaJlplIyx4bXOOIhUpXSWZAFYtbaK24jzP+zJVb6fUxU99DdjUP7Xppome7/o0RVQky4o3ARzsCF9A2ms9LyQbfvtPi2Jxlp9FdwklWnEeXUcoeoaQGZOVRNUdJJIdxZkAIjX+Tx4BXPjfZ++MSQLFtDirECRsWc0wTXSykqim6IpkXHHWAvrqxRMX/zg1C+Dzb6QfmKQs2TpCDaNj7O7qjjkWcsZlVdt5Zf1HsxmiSFFQLaBEZR2yWO5BppdsHSGvLqEzlUS14bxIJhVnAki0QXsKG7dL5iU1TTSuS+i56urJKaLacF4ks3KWAMzsI2b2spk9ZmYfyeqLXbExXOI5Woobt0t2JbPZzHAgwMPnxq4cViVRkcxJKwGY2RNmdszM3og7fqWZHTKzt8xs/cRhB5wCKoHs9sW0rA7X969ZAFj43yzW+5fUJdsl5LVyOFJJVF1CIulJaxDYzC4nfFF/yjl38cSxMuDnwMcJX+j3AjcB/+ScC5nZ+cBDzrlbZnr+TG0II/6XaNOcgAsXlYivJhqhzWVEYuVsENg59xJwIu7wpcBbzrnDzrkR4GngWucm9xJ8G5iTzutK8UnUJRSaYZqouoREZi8bYwBNQPQczC6gycyuN7NvAv8b+HqiB5vZGjPbZ2b7+vr6shCe+FGyK4e9qomqS0hkdrKRALw6dJ1zbrtz7jbn3O87515M9GDn3OPOuWXOuWV1dXVZCE/8KtmVw16rhjVLSCR12UgAXcCCqNvNQHeCc0U8TbdyeLrNZdQlJJK8bCSAvcCFZrbIzCqAG4GdWXgdKWIzTROdbnOZRjsesw+0iHhLdxrod4AfA0vMrMvMbnXOjQG3A7uAnwEdzrk30w81DTMUhhP/SWa7yUSF5EIYv5hzM73t71UhOZFpFGctoGiRwnCjUZ8Iy6u0LqCAJJoiGq0yFKK9/wSfPDUYs7/AkKvgjQ/cpx3GpGSoFlC05++NvfhD+Pbz9+YnHklZsquGt5xTq81lRFJQ/AlAheEKXjqF5CC8uUzonhp1CYnEKf4EoMJwRWE2heQiSSAyTVT7DYvEKv4EoMJwRSeVLqF46hISOSOY7wCyLjLQ+/y94W6fmubwxT9uAHjHgaNs3nWI7oEhGmurWLtyiWrL+FRk17Ytr22h93Qv4WpBU3kVkgM4z/WzaH2nfs9S8op/FlASdhw4yobtrzM0Oj55rKq8jAeuf58uDgUg1UJyzsFRN58Hx1bzw7Lf0e9ZiopmAaVo865DMRd/gKHRca0mLRCpFpKzibIRm8q38vHxH+n3LCVLCQASrhrVatLCkGwhuS/WzYuZJVRtI6wLdnB0YIhF6zv58KYXVENISkrxjwEkobG2iqMeF/vG2iqPs8WPWhe3To4NtDzZ4nlOaGKRQKRFAPDJU/0cnnMz3W4+D76zmg3bRwDUJSQlQS0AYO3KJVSVxw4YVpWXsXblkjxFJOlIWEguSvTCsYC6hKREKQEQ/rT3wPXvo6m2CgOaaqs0MFjAkpkmClNnCalLSEqNZgFJUeo83Dk5TdTMCE1uSHeG1ywh58KbV3drlpAUqFRmAWkMQIpS9JhA5+FO2l9tZ3h8+MwJznmOCbSeHsSAZgt3CTEKm3dVKAFIUVILQEpCMi2ChtExdnfF7l005gIECNHt5rO14g94f+saJQPxNa0DEImTzHaTXiuHgxbSdpNStJQApOQkmiXkYEol0WjablKKjRKAlJyEs4Rm2G8YwttNHh0Y0gwhKQpKAFJyZtpu0mvVcEQI4/Ccm/nu4B+rO0gKnhKAlKTImECiDWa86gg5FzsmcK89zk87H89x5CKZowQgJS3VVcPRqm2Ez4z8tRaNScFSApCSluyq4UTbTTZZP79Ql5AUKCUAKWnJVBIFZtxuUl1CUoi0EEwkiueqYQ8zbTajRWOSL1oIJjJL8S2CRGbabEaLxqQQKAGIxIleNZxoqmg0rw3otWhMCoESgMg0ZltaGsIDxC8PXUfXxvfQft89ag2I7ygBiEwj2UFigymzhKIHiNUlJH6kBCAyg+guofuX35/SBvQR6hISP9IsIJEUJVNaGudo8JghFNlsRrOEJFtSmQWkBCCShpYnWwhPCPVWGQrR3n9iMglEDLoKNro1LL/uc0oCklGaBiqSIzOVkvCaIQRnuoTu+u5PVUZC8iZnCcDMft3MHjOzbWb2J7l6XZFsSmaW0HRlJFRZVPIprQRgZk+Y2TEzeyPu+JVmdsjM3jKz9QDOuZ855z4LrAaSap6I+N1MpaUBlZEQ30q3BfAt4MroA2ZWBjwCfAK4CLjJzC6auO8aYA/wfJqvK+IbkVlCmy7bNGNrYLouobtHHyZ0Tw297e9l785vZitckUlpJQDn3EvAibjDlwJvOecOO+dGgKeBayfO3+mc+xBwSzqvK+JHyZaRSNQlFNlroJ4+Lt5/t5KAZF0wC8/ZBByJut0F/JaZfQS4HpgDfD/Rg81sDbAG4IILLshCeCLZ07q4ldbFrQCs2LaCntM9U08yw8FklxAwZZZQlY3QuP9BFr3aTGNtFWtXLtFsIcm4bAwCe330cc65F51zdzrnbnPOPZLowc65x51zy5xzy+rq6rIQnkhuJDNAPN32k43083LFnXzgnR+yYfvrGiSWjMtGC6ALWBB1uxnozsLriPhapCUQWTSWaL1AaGKrsfgWgRk0Wz+byrfCKGzeVaFWgGRU2gvBzGwh8Jxz7uKJ20Hg58AVwFFgL3Czc+7NVJ9bC8GkmCTsEorTMDrG7q7Yz0xjLkCAEIHaBXDFRmhZna0wpcDlbCGYmX0H+DGwxMy6zOxW59wYcDuwC/gZ0DGbi79IsUln+8nIADEnjzC0/XYNEEtGpNUF5Jy7KcHx7zPNQK9IKYrvEkpYR2iGQeIq3qVx/4PsWHCVuoQkLaoFJJInyW4/6dUlpMJykohqAYkUgHTWDWivAckEtQBEfCKZQeJE1UW7QvP5ndE/J+Sc1g2UOLUARApQsusG1tfN8yws988VN2ndgKQkG+sARGQWkl03QFRhOTizZsDQugFJjVoAIj4Svf3ktBVGmb6w3LpgB0cHhrTXgExLCUDEp7TXgGSbEoCIT2mvAck2zQISKQDJrhkIuPDIQX3chvSRUhLHrI4jl6zlg9fcloOoJR9SmQWkQWCRApBuYbmghVcc19NHzf672QtKAqIWgEghSrawXKIWgQOsRoXlipHWAYgUuWQLy4USjRGACsuJEoBIIYovIxGwmf8re00bnSwsp1lCJUldQCJFINlBYpzDiO0SUmG54qIuIJESk3SLwKNLSIXlSpdaACJFKJ1po12h+fx+9V/wyvqP5iZYyShNAxUpcelMG/3kqX5eHrqOro3qEip26gISKVKp1BWCM4PE6hIqHUoAIiUg2WmjvcGymNvVNsIme5RrnvkN+OrFcLAjWyFKHigBiJSAZAeJHUwpLBe0EAGc1g0UISUAkRIR3SV0//L7vVsEHovGomndQHFRAhApQTNVGh0OBPhi3TzPUtMNHGfzrkO5ClWySAlApERFWgSJNqRPVEYihPHy0HX0tr9X3UEFTglApMTVz62f8ZzIDCHnJsYELFxZ9OL9dysJFDAlAJESl+wMoZ5gGb+5KLZLqMpGWLZ/HV0b30P7ffdobKDAKAGIlLh0ykhMHNaagQKlUhAiEkNlJAqbSkGIyKylU0biE6eO0z0wlJtAJW3qAhKRKWZbRqLbzSNgxqL1nXx40wvqDvI5JQARmVYqZSSqbZh/rriJlyvu5APv/JAN219XEvCxnHUBmdli4M+AGufcDbl6XRFJT3yXkJkRcqEp5xnwkUXnUD92Nm1vD7Dp1FYYhc27KlRN1KfSGgQ2syeAq4BjzrmLo45fCWwByoCtzrlNUfdtSzYBaBBYxH88B4mdC08HmlAZCtHef4KVp4YJECJQqw3ocyWXO4J9C7gy7sXLgEeATwAXATeZ2UVpvo6I+ITntFGLXU0cGROILBpTITl/SisBOOdeAk7EHb4UeMs5d9g5NwI8DVybzuuIiL9EDxIn6kXoCZbFLhpTITnfycYgcBNwJOp2F9BkZvPM7DFgqZltSPRgM1tjZvvMbF9fX18WwhORTEpYSsJj0ZgKyflLNhKAV2Up55w77pz7rHPuPc65BxI92Dn3uHNumXNuWV1dXRbCE5FMmmmWUKQ7CM4UktPmMv6QjVlAXcCCqNvNQHcWXkdEfCB6llDP6R7Pc3qCZbQsXHBm1fDJI/DsneE7NTCcN9loAewFLjSzRWZWAdwI7MzC64iIT0TGBBIuGvOqIzQ6xOAPNuY2UImRVgIws+8APwaWmFmXmd3qnBsDbgd2AT8DOpxzb6Yfqoj4XTKLxqK7hCoHezUonEdpdQE5525KcPz7wPfTeW4RKTzJ1hGKdAnVjTne+dG3WbV0XS7DlAkqBSEiGZVUHaGJLqFj5QGGap6m83BnboMUQAlARLIomS4hC4yy5cUvaFZQHqgctIhkTewMoV6cc/GLhgHoDcDYM3eEL0iaFZQzagGISFZFuoRe//RBaivO8zynfmyc4PiwZgXlmBKAiOTMht/+0yldQpWhEG1vD4S/H+rNR1glS11AIpIzkS6hh15YR1/QYraTBOgOzaM5nwGWGLUARCSnWhe3clnvDfzkl//G7q5wkYAVzY20LFzA9b82TzOCckgJQERy7v2ta9jo1vBU9fm0zz+XnvIgzoyh4BDtP1pH5yOqFZQLSgAiknOrljax/LrP8T/n1TMciL0MDQcCbJkzztgzdygJZJkSgIjkxaqlTRAc8LyvN1imWUE5oAQgInmTaC+B+rFxQLOCsk0JQETyxmulcPS00O7QvHyEVTKUAEQkbyL7C1eNVWHO0TA6Rnv/CVpPDzLoKtha8Qf5DrGoaR2AiORV6+JWRk++nz3fe5S7eJpGG6LLzedr3Mjy1jX5Dq+oKQGISN6tWtoEfI7f33UFx0KvUnX+blzwWQ7+Yi/lNW2TC8gks5QARMQXVi1torzmp7S/+gzD48MA9Jzuof3VdgAlgSzQGICI+MaW17ZMXvwjhseH+dLff0k7h2WBEoCI+Ebvae9pn8Nlg+z53qNKAhmmBCAivjHduoC7eJrNuw7lOKLipgQgIr7RdkkblaFQzLHIuoBGO073wFCeIitOSgAi4huti1tpHzQaRsemrAvodvNorK3Kd4hFRbOARMRXWi/byMpn7iAYNRg86Cr4GjeyduWSPEZWfJQARMRfWlYTBAZ/sJHKoV7+uvo8vj6vliGtC8g4dQGJiP+0rKb6C//ED/7w2/x5Yw1DwSHATa4L0KYxmaEEICK+lWhdwJbXtuQpouKiBCAivpVoXUDPqR6tCcgAJQAR8a1E6wLmjFWxYfvrSgJpUgIQEd9qu6SNSiuPOVYZCvHFEz18fPxHWhiWJiUAEfGt1sWttP/7yJR1AdcPvsO6YIcWhqVJ00BFxNda+7poxU053mjHtTAsTWoBiIi/1TR7Hu5hnhaGpSlnCcDMFpvZX5rZtly9pogUgSs2QnnsJ/0h5tD9gXUTG8nIbCWVAMzsCTM7ZmZvxB2/0swOmdlbZrZ+uudwzh12zt2aTrAiUoJaVsPVD0PNAsCgZgFV13+dD15zW74jK3jJtgC+BVwZfcDMyoBHgE8AFwE3mdlFZvY+M3su7uu8jEYtIqWlZTV8/g06//DbrFjQSMuB+1ixbYVWBKcpqUFg59xLZrYw7vClwFvOucMAZvY0cK1z7gHgqkwGKSLSebiT9lfbtV1kBqUzBtAEHIm63TVxzJOZzTOzx4ClZrZhmvPWmNk+M9vX19eXRngiUkxUFiLz0pkGah7Hps7Vitzh3HHgszM9qXPuceBxgGXLliV8PhEpLYnKQiQ6LjNLpwXQBSyIut0MdKcXjoiIt4TbRY6H4GBHjqMpDukkgL3ddAdmAAAGrklEQVTAhWa2yMwqgBuBnZkJS0QkVtslbVSWVcYcqwyFaDt+HJ69U0lgFpKdBvod4MfAEjPrMrNbnXNjwO3ALuBnQIdz7s3shSoipax1cSvtH2qnYdxN2S6S0SF4/t58h1hwzDn/drMvW7bM7du3L99hiIiftNfiPdxo0D6Q62h8x8z2O+eWJXOuSkGISGFJUBoi4XFJSAlARAqLR2kIyqvCxyUlSgAiUlg8SkNw9cPh45ISlYMWkcLTsloX/AxQC0BEpEQpAYiIlCglABGREqUEICJSopQARERKlBKAiEiJUgIQESlRSgAiUrA6D3eyYtsKWp5s0RaRs6CFYCJSkLRFZPrUAhCRgvTATx7SFpFpUgIQkYKz48BRBkaOed6nLSKTpwQgIgVn865DuNFaz/vGR2r48KYX2HHgaI6jKjxKACJScLoHhni3byUuVB5z3IXKebdvJUcHhtiw/XUlgRkoAYhIwWmsrWLsnaUM91xPaKQWHJw3GuKBvh5eHP4rrgnsYWh0nM27DuU7VF9TAhCRgrN25RKqyssYe2cpV/xyOf/3l70839XF1YODNAf62VS+lWsCe+geGMp3qL6mBCAiBWfV0iYeuP59NNVWsS7YQbWNxNxfbSM8VP4Yv6i8Bb56MRzsyFOk/qYEICIFadXSJl5Z/1GaA8c97w9aiAAOTh6BZ+9UEvCgBCAihS2ZzeBHh+D5e7MfS4FRAhCRwua1SbyXk13Zj6XAKAGISGGL3yTeyrzPS6alUGJUC0hECl/0JvEHO8J9/qNRM4DKq8ItBYmhFoCIFJf4FkHNgvDtSIKQSWoBiEjxiW4RSEJqAYiIlCglABGREqUEICJSopQARERKlBKAiEiJMudcvmNIyMz6gH9JcHcNcDLJ+6a77fV9/L/zgf4kw04lrmRi8UNcycYzm7imiy3duOKPlWcorpniSBSXV4z6G5s+Tv2NpRYXwIXOuZqkXt05V5BfwOPJ3jfdba/vPf7dl424konFD3GlEE/KcU0XW7pxxR/LVFz6G9PfWCH8jSXzVchdQM+mcN90t72+j/83W3ElE4sf4ko2ntnENd3j0o0r0bF045opjkQxeMWjv7HU4tHf2PRxpRSHr7uA/MLM9jnnluU7jniKKzV+jQv8G5viSo1f40qkkFsAufR4vgNIQHGlxq9xgX9jU1yp8WtcntQCEBEpUWoBiIiUKCUAEZESpQQgIlKilABmwcwWm9lfmtm2uONzzWy/mV3ll7jM7NfN7DEz22Zmf5KPuKaJbZWZ/YWZPWNmK3wUl+fv1wdxzTWzJyfes1vyFdtELBeZWYeZfcPMbshnLPHM7AIz22lmT5jZ+nzHE2Fml038X9xqZq/mOx6gcBeCzfYLeAI4BrwRd/xK4BDwFrA+yefaFnf7XuALwFV+imviWAD4S7+9ZxPHzplNbDmIa8qxfMYFfAq4euL7784mtkzFB/w34LKJ73fONpYsxfYx4LaJ75/yS1xRj1kViS/fX3kPIOc/MFwOXBL9iwTKgF8Ai4EK4B+Ai4D3Ac/FfZ0X9bjo/5wfA24E/vMsE0BW4pq4fQ3wKnCzn96zqGP/C7jEh3HNNgFk629sA/D+ie//Jp//Bya+HgE2A6/46f8nMA/4e+AF4I/8ElfU4zqAszP1nqX1c+U7gLz80LAw7hf5H4FdUbc3ABuSeJ7o/5xfBr4G7AaeAQJ+iCvueKfP3jMDvgJ8zE9xzfQ+5vH9+hQTHy6Ap2cbW4bjKwOeSSeWTMcG/Hfg8nR/h9l4z4ALgL/I5PuVzpfGAMKagCNRt7smjnkys3lm9hiw1Mw2ADjn/sw5dxfwN4R/wSE/xGVmHzGzh83sm8D3MxBTxmID7iDccrrBzD7rl7gSxJr3uIDtwO+a2TeYfUmETMW30MweB54i3ArIppRiA/4WuHPi/fuVj+ICuBX4q6xFlCLtCRxmHscSrpBzzh0HPC9YzrlvZSgmyEBczrkXgRczGFNEJmJ7GHjYh3El/P3mOa7TwB9lOK6IVOP7FbAmS7HESzW2N4BcDEynFBeAc+6eLMUyK2oBhHUBC6JuNwPdeYolml/jAv/Gprhmx8/x+TU2v8aVNCWAsL3AhWa2yMwqCA/m7sxzTODfuMC/sSmu2fFzfH6Nza9xJS/fgxC5/gK+A/QAo4Qz+K0Txz8J/JzwqP6fKS7/x6a4ii8+v8bm17jS/VIxOBGREqUuIBGREqUEICJSopQARERKlBKAiEiJUgIQESlRSgAiIiVKCUBEpEQpAYiIlCglABGREvX/ARohiyxKzuBeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Is it better to pass a string or a float?\n",
    "s001 = Sdist(0.01, mainDir, is_avoid_zeros=False)\n",
    "s002 = Sdist(0.02, mainDir, is_avoid_zeros=False)\n",
    "s003 = Sdist(0.03, mainDir, is_avoid_zeros=False)\n",
    "plt.loglog(s001.size, s001.ps, 'o')\n",
    "plt.loglog(s002.size, s002.ps, 'o')\n",
    "plt.loglog(s003.size, s003.ps, 'o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/home/gf/src/Python/Python-in-the-lab/Bk/F64ac_0.01_S.dat'"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s001._filename"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['size', 'ps', '_filename'])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# s001.ps, s001.size: Nothing else?\n",
    "s001.__dict__.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### You seem to be able to do it for all the frequencies..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14 lines deleted\n",
      "12 lines deleted\n",
      "14 lines deleted\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{0.01: <__main__.Sdist at 0x7f8ce66be630>,\n",
       " 0.02: <__main__.Sdist at 0x7f8ce6931ef0>,\n",
       " 0.03: <__main__.Sdist at 0x7f8ce66bea58>}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Can we now do it for the 3 frequencies?\n",
    "freqs = np.arange(1,4) / 100\n",
    "# Can I make a dictionary?\n",
    "s_distributions = dict()\n",
    "for freq in freqs:\n",
    "    s_distributions[freq] = Sdist(freq, mainDir)\n",
    "\n",
    "s_distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([11991900., 15108900.,   338738.,  2429310.,  5685190.,  2692080.,\n",
       "        3621790.,  5915990.,  6721880.,  1033120.])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s_distributions[0.03].ps[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEACAYAAAC6d6FnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt4lPWZ8PHvPYeQg5oAoeYwUUC6rIqpCOpe1dqqFWyjFK1Freu6W626Wwq1uwi0LqZoCxVf21h1rcWuuq8VKaUcmir4apWiXQsKi4pLS1EvhiSFRBKUJCQz83v/mMwwM3kmmfMhc3+ui8vMM888cycmzz2/0/0TYwxKKaUKjy3bASillMoOTQBKKVWgNAEopVSB0gSglFIFShOAUkoVKE0ASilVoDQBKKVUgdIEoJRSBUoTgFJKFShNAEopVaAc2Q5gKJWVlWb8+PHZDkMppfLGG2+80W6MGRfLuTmdAMaPH8/27duzHYZSSuUNEfkg1nO1C0gppQqUJgCllCpQmgCUUqpA5fQYgFKqcPX39+N2u+nt7c12KDmpuLgYl8uF0+lM+BqaAFJk3Y4DrNi0h5bOHmoqSlgwczKzp9ZmOyyl8pbb7ebEE09k/PjxiEi2w8kpxhg6Ojpwu91MmDAh4etoF1AKrNtxgMVr3+JAZw8GONDZw+K1b7Fux4Fsh6ZU3urt7WXs2LF687cgIowdOzbp1pG2AFJgxaY99PR7w4719HtZsWlP1loB2iJRI4He/KNLxc9GWwAp0NLZE9fxdNMWiVKp8fzzzzN58mQmTZrE8uXLLc85duwY1157LZMmTeL888/n/fffB6Cjo4OLL76YE044gblz50Z9j8997nNh653ef/99pkyZktLvI5oRnwCa9zUzY80M6p+sZ8aaGTTva477Gut2HOCC5S8xYVEzFyx/adCNtKaixPJ10Y6n21AtEqVGquH+TuPl9Xr5xje+wXPPPcfu3bt55pln2L1796DzHn/8cUaPHs3evXu54447WLhwIeAfpL3nnnu4//77k4ojnUZ0Amje10zja420Hm3FYGg92krja41xJYFYPk0vmDmZEqc97HUlTjsLZk5O1bcSl1xrkSiVbulo9f7xj39k0qRJTJw4kaKiIq677jrWr18/6Lz169dz0003AXDNNdfw4osvYoyhrKyMCy+8kOLi4oRjuOWWWzj77LM5++yzGTduHN/73vcSvpaVETcG0LyvmaY3m2g72oaI4DO+sOd7vb00vdlEw8SGmK4XS/9+4L+p6HNPpO8+8jUVpU4Od/cPOi9bLRKl0i0d43AHDhygrq4u+NjlcvH6668PeZ7D4aC8vJyOjg4qKytjfq8bbriBkhL/32dfXx82m/+z+cqVKwH44IMPmDlzJv/4j/+Y0PcSTcYSgIh8Brhh4D3PMMZ8OtXvEfjE3+v1j4wbYyzPazvaFvM1Y/00PXtqbdgvWqA5Gnojh6GTROBTTOAX+UBnDwt++T98b+M7dHb3x/wap01w2oV+7/HvP5stEqXSLR2tXqv7h9XAa6znDeXpp59m+vTpgH8M4Iorrgg+19vby1e+8hUeeughTj311LiuO5ykEoCI/By4AjhojJkScvxyoAmwAyuNMcuNMb8Hfi8is4FtybxvNE1vNgVv/kOpKquK+Zo1FSUcsPglGurTdLQbOULwpjztyAucu+7rmPUdSLkLLl3Cik2Vgz7F9PtM8NN8oFkL4a0Oq9dUlDgpG+VIqiWhM4dUvkjk73Q4LpeL/fv3Bx+73W5qamqinudyufB4PHR1dTFmzJiE3zfS7bffztVXX83nP//5lF0zINkxgCeAy0MPiIgdeBj4AnAGcL2InBFyyleBZ5J8X0uxfLIvthcz/5z5MV8zkf79aDflwM1/lm0ry50rqZV2BANd+2HjPKYfeWHYeC7zvsLfrf8sNFbAj6Yw/cgLzLJtZWvRPPaN+ipbi+Yxy7aVrp5+Xl10Ce8tb+DVRZfEdPPXmUMqX6VjHO7cc8/lz3/+M++99x59fX2sWrWKWbNmDTpv1qxZPPnkkwCsWbOGSy65JGXTVx9++GE++ugjFi1alJLrRUqqBWCM2SIi4yMOnwfsNcbsAxCRVcCXgN0icgrQZYw5ksz7RlNVVkXr0dZBx21iwxhDVVkV88+ZH3P/P8Bs+6vMOGEJxT1ttPjGsrLo7/FMuYYVm/Zwx7M7LT8pD9fsvNOxmlLpCz/Y38Piol+yvvfCqK8LJI5SBl7btZ8VRT8FIxSJBwCXtLPcuZIxziLg+Pc53Kf7WPtQtZWgclEqx+ECHA4HDz30EDNnzsTr9fK1r32NM888E4AlS5Ywffp0Zs2axc0338yNN97IpEmTGDNmDKtWrQpeY/z48Rw5coS+vj7WrVvH5s2bOeOMM6K95SD3338/TqeTs88+G/C3Bm6//faEv6dIEq2fPOYL+BPAbwJdQCJyDXC5MeaWgcc3AucbY+aKyPeATcaY14a43q3ArQCnnHLKtA8+iLm09aAxAPB/4m/8dGNcN/2gXath4zzoP35D99iLWdR/C2v6jg9hlDjtLLv6rOAv2wXLX7Jsjgb8ZMwt/GRMBW0OO1UeL/MPd9JwtBuDcIZ31aAbccDWonm4bO0xhe4TGzZjoNzFttO+yT9sOzXsutcUvcbSsl9R2tMG5S7mH7qS9b7ByUeA95b7f3aRXVtW37tSqfLuu+9y+umnZzuMnGb1MxKRN4wx02N5fToGga3aPgbAGHP3cC82xjwGPAYwffr0uLJT4CYfmAWUyCf+MC8uDbv5Azi8vXyLVazheAKI/KS8YObk4I3ScdIORo3bhDg7MZ4KPB/9Ld+rGMsxm//H1Op00Fjp7y88++NSvjytlt/97yFaOnsoL3FytM8T7Dqqkdhu/gC2wOynrv1MefPfucx7Mxvw3+Bn2bayVFZS2nO8JbHMuRLTDxsikkB5iTM4mG0TwRvxgSHbK56VUolLRwJwA3Uhj11ASxrex1LDxIbEb/iRutzBL5vLSmkaffxTu+PgDjxHpgafD+32CdwMv//K0/SUr0Vs/kFccXZSNPoPHIvoH+y12fjx6Ar+vusM5u68intph5P9A8PrvBcEm7UHZRxVHIr72yjhGHc7nuJOVlMj7fiw4ZDw6bGl0sedjtVs6DueAJw24Wifh84ef/yRN3+r710plT/SkQC2AZ8UkQnAAeA6/AO/+afcBV37aS4rpbFyDL0Dc3NbnQ6Kq9fSC8EkEDbbYNdqZr+8lEdO9NJri/gRRxkcanM4+Ip9S1j/PhvnMfvKB5m9aM7AdY/S/atvhI0fHDN2BP8YQGSSCnQtAYyRjxH5GAAb4Tf/gBpp59VR86imnYMyjiau55nevws7Z5ZtK3c6/ImkxVRyn2cOW0ZdPGjKq7YIlMp9Sc0CEpFngD8Ak0XELSI3G2M8wFxgE/AusNoY807yoWbBpUvAWULT6IrgzT9AbP2MGrcJiJhtEBg36NpPm8MeecWoPuHxWQ4M8+LS44/r53Cf819w+yrxGcHtq2RB/238W/+tPF1WRWPlGFqdDoxIsGupuazUH2+MkxJqpR2bQBWH+HfzKLNsW4PPBQahXTb/OS6bf8D54v6XdfaQUnko6UHgdJo+fbrJ+qbwu1ZT/+ZSjOUCEChvawr/xPujKf5P78AMVw2tTotGljFhd2Tjc7LsUCtXdndbBCDQ2Blc4dx6tA3TX07vwZnB1keJ087Y01fQ1X9w0KvLPV5KjbFsFYTyGbBZJAmPsWHDR4uppFR6GTPQioh2zn2eOWzwXYhdBJ8x2iJQCdNB4OElOwg8omsBpUT9HKpOGLz4A6DmhOpBc+ybPR8yw1VD/fg6ukVwRiTYYp+Pa498RHW/BzEG8YzmWOvVnNtdZv3+5a6wmkZgEGcnJdVrcZ60g9qKEpZdfRZH+q3HBrrstrBWwV2VY/jMKbXUj69jhquGjaVluH2VYSP3zWWlwe/hi3VVPHdCKS5bO6MZfPMHcIgvrEUwy7YVrzHaIlAqx2kCiMH8c+ZTbA8v6FTs8zH/ry3+Lp8BzfuaaRx3vBumy2HHGEOF14sYQ3W/h8b2D7nrw042//Uwu85Zwq6bt/Dn79xF1dU/AGfEqkVnCXxyBk0vLxy8wtnWz4S/2RJMQFFXN0e0XDw2G512ezAhfGdcFRcX/xN/lXH+72FgvMOqKymWbqRS6eNux1NhC9Mu876ilUhVXkqmHPQLL7zAtGnTOOuss5g2bRovvfSS5euzWQ56xBWDS4fg9NL/XkZbX2ewKwVgxrZG2nbcQ1VZNT2eHnotbrglHi+/r5wBf94MR3ugvM4/vlA/5/iJga9fXOqffVTugk/OgP/5BW211kWlQlc+zz9n/qA1EDEZSCRVf/MD2DjPcryj12ajaXSFZdeRldAB58DCtEVH4ILl6ECxSp9dq8P/fiL/xuIUKAf9wgsv4HK5OPfcc5k1a9aghVyh5aBXrVrFwoULefbZZ6msrGTjxo3U1NTw9ttvM3PmTA4cyK2WsCaAKEKrigbWE2z+a2ewf3/QzCCLFcgBbQ47XPHA8G9aPyf8F/ZHU6C/hyqP13IsIfRTv9UaiB5PD53HOod929ajrcz400raXOOiF9ALDGiXjIGiMv8fmdjADF60FtlSKJU+HnA+iq3nEVqKKrnvyBwWr/UPeGsSUCkRuWhzYBYdkHASCC0HDQTLQUcmgPXr19PY2Aj4y0HPnTsXYwxTpx6fJn7mmWfS29vLsWPHGDVqVMwx3HLLLcHWwYEDB5g7dy533z3scqqYaQKwELmiOLCPAJ4Pg8UVrD4pR1NVVp1YIAPrEOYf7gxLNmBd0yhyDYTVyuhoggksSj9Plcfr75L6wg+P/0FZrJQ2WK8EDKw7CLQIpnn/xLnrduJb559yuv+cBZw767Zh41TKksWizeAsugQTQCrLQf/qV79i6tSpUW/+I74cdF4YaEI2neilN+ITd6+3l6axY2g4ehQg+hTPiBk+8RafCzOwDiHQ9RKc4++D+Z8ZvrxFZKvgpKKT6PZ00+8bvFfAUIp9Pi7y2pkxYRJtO+6l6k8ruch1EVvcW2hzjaPKa5jf8SGUjqGp2NBmlyFnHJVKHzfa/19w1lEVhyh/4y62gSYBlZiQRZsxHY9BqspBv/POOyxcuJDNmzdHfa+8LAc9ooR8mm0bXWd5Sptd/J+Ch+iWKff6jk+7PKEmuVIUly4JxtRwtNt/M3WWwJUPQozXtGoVhHYTDdV1JQhVZVVc5LqI9XvX09vfBfhbC8/ueTZ4XqtduOvkkxGRYHIJLXFhlQQip5yWSB91b64ATQAqEQMfliyPJygV5aDdbjdXXXUVTz31FKeddlpCceRyOeiRI6QJWeWxLsZWVVbtv/mW1zH/cBfFFlM8F394mM3uFna972bzNZsTvvk372tmxp9WUu8ax4xT6mguK/MPHl/5YFIDWw0TG9h8zWZ23bSLzddspjpK91R1WXXwnC3uLcN2I3mMZ1DLotdmo2lMRcyxnWwO4WuswL3kNBrvvVunjqrYDSzaDOMs8R9PULLloDs7O2loaGDZsmVccMEFCcWQ7nLQmgACQpqK8w93UuwLL5cQ7MqpnwN3vE3DghYaL/oh1V4TNsUz+Gk3iU8e4XsZ+z9hN1bX0vylH4bd/FOx4b3lFNeIbqt4dlCL1OpwBNccBFYl+6KsPRQBGwaXrZ07+x9h668f0SSgYlM/J/jhDCQlH5ZCy0GffvrpzJkzJ6wc9IYNGwC4+eab6ejoYNKkSTzwwAPB6aIPPfQQe/fu5Z577gnu63vw4ODFmkO5//77eeutt4Kvf/TRRxP+fqzoSuCAkBW8EFH8baiuHIuB0GA3TYK/fDPWzLDsmqkuq2bzNf5+xFSWvraa8RR6jWjxxKvY56OxW/jb0guoeX8tJZGlLyK4fZVcW/ozXl10SdLvrfKPrgQeXrIrgTUBBCRzI0/x/OP6J+vxr6MNJwi7btoFxJYkUiWW2UQOcYSNAURTXlROqbOU1qOtfKLfx7c+PMwV3UctZw75DLSYSly2jpT8XFV+0QQwvFzcDyA/WS3EivWGEzl/P06Rn8DLR5Vbzt8PnfcfrVsmme6aaKzWGARnAYW0GkLPsUpgAF19XXT1+QeTDzptLK2pxdbVQ8MhiwE8OL4BTtd+etbO5e33D+tMIaVSRBNAqCRv5ImwWnPgEAdOmzPs03Rkv3y0GTzxbHgfj1j3WQicE2u3Ua+3l6bRFcz88BCOkBaGVXG6Eo5R88Z9rKu7QheQKZUCOgicZU1vNg3qWvEYD6WOUqrLqhGE6rLqQX37sQzeZpNVfNG09ndxzimf4LK6OjaWDRSni1J3qJoOrSukVIoUVAsgckPzH5/xZ879y09S1nefiGhdNkf6jrD1+q2Wz0Eatr9MsXhLUxigzSEsra6l8cuNuNYvtJzX7UP4fc9VuJdUsrLo7zm74VZtDSiVoIJJAJEbmk878gJT3lgJEr4DF5DRJJBMV05Kt79Mg0RKU/R6e1n2+jKaTq6gbTRhK4qNCS8pcWf/Iyz5tQf4F00CSiWgYLqAVmzaE7z5A9zpWD14GmLkDlwZkOtdOanUMLGBxk83Bru2ounq66K1vyusHPWG0jLLInPfYpV2Cam0SaYc9B//+Mfg/P1PfepT/PrXv7Z8/fjx42lvbw8+fvnll8NKQaRTwSSAyI3La6Td+sQkaockIvKmaNXfP5KErkSOtgo5Uq/NxkNjyi2fq5F2nu3+OjRW+NdyhOzPoApLKhZGhgqUg37uuefYvXs3zzzzDLt37x50Xmg56DvuuIOFCxcCMGXKFLZv387OnTt5/vnnue222/B4PEnFlGoFkwDCNm3HP7/cSndJembRDCWyPMNIvflHimug2GEftKI4wD9V1ASnim7b8NM0RKtyWfjqeROs4JtMEggtB11UVBQsBx1p/fr13HTTTYC/HPSLL76IMYbS0lIcDn8ve29vr2UhueF88YtfDLYiysvLgyUnUqVgEsCCmZMpcR6v4HmfZw7dpijsnG5TxJKjX9byAxli1fqpGBWldpDIoB3Khpwqqv8PC4rVbLpeby9NbzYlfE2rctBWG7pEKwcN8Prrr3PmmWdy1lln8eijjwYTQqSLL744eKO/5ZZbgsd/+9vfsnPnTh5//HFOPfVUZs+enfD3Y6VgEsDsqbUsu/osaitKEOCNky7jHrkdt6+SjaVlXOpycf6EKp6r28L3X3k62+EWjMjWz6LzFg3bKui12Xhg9BidKqqC0rEwMhXloM8//3zeeecdtm3bxrJly+jttZ4A8bvf/Y6dO3eyc+fO4B4AAe3t7dx444384he/oLzcuis0UQUzCwj8SSB0tsiERT388qRRFI9bi9j8i65sRZ30OFbRvO/MgumKySWR00ejrSg+5LQj5XU6VVQB6VkYmYpy0AGnn346ZWVlvP3228G6/7Hwer1cd911LFmyJC37BBdMC8BKTUUJo8ZtCt78A8TWn1TTUSUnloHik4pOYsbJFYPGBQJTRW2CVhUtIOmYTZdsOej33nsvOOj7wQcfsGfPHsaPHx9XDIsWLaK+vp7rrrsu4e9jKAWdABbMnIw4rRcmpaOmjoqf1R+2Qxx0e7p1qqgKSsdsumTLQW/dupVPfepTnH322Vx11VU88sgjYdtExuL+++9n8+bNwfGBwHumSsFXA73wF5fS1T+4Rnc6qmqqxEQWy4u2ori638Nmd8ug41pVND9pNdDhJVsNdOS2AHat9s8LH2Z++OK/+3bBLMTKV5EDxV3HuizP06miSsVnZCaAQG3/rv0E/ujZOM8yCRTaQqyRIOrAnk4VVSouIzMBhOzvGzREmYdCXYiVr2JZQKZTRZUa3sicBhqtnEOGyzyo9Ih1quhBp42Zp9TRZgsvKgc6VTRfGGMSWkFbCFIxfjsyWwDRNmRPYqN2lVtirSnUah/cLaRTRfNDcXExHR0dKbnRjTTGGDo6Oigujq2USjQZawGIyOeAe4B3gFXGmJfT9maXLrHe3/fSJWl7S5U988+ZP2yZafB3CzWNrgi2AgJKpY9vmVVcu+lSbQXkEJfLhdvt5tChQ9kOJScVFxfjciX3oTapBCAiPweuAA4aY6aEHL8caALswEpjzHL8e358DBQD6e2LSWZ/X5V3rDafibYdZZvDbnk8UFXU19hBi2+sdgvlAKfTyYQJE7IdxoiW1DoAEbkI/039qUACEBE78CfgMvw3+m3A9cD/GmN8InIy8IAx5obhrp+JdQBqZIq2J7HN+EcMBo0JRMwW6jZFLDG3cuFVutmMyi8ZWwdgjNkCfBhx+DxgrzFmnzGmD1gFfMkY4xt4/jAwKpn3VWo40WYK+WKcKqoriFUhSMcgcC0QWqHLDdSKyNUi8lPgv4CHor1YRG4Vke0isl37/lSiItd32GTwr/pwU0VrpX1gptBpNN57tw4SqxEn6VIQIjIe+E1IF9BXgJnGmFsGHt8InGeM+Wa819YuIJUq9U/WR50uWu01llNFQ2mXkMoX2S4F4QbqQh67gMEFWpTKoKHKAltNFY2kXUJqJEpHAtgGfFJEJohIEXAdkNoSdkrFKdbtJwNTRa3USMegvaWVymdJJQAReQb4AzBZRNwicrMxxgPMBTYB7wKrjTHvJB9qkmIsDqdGJquaT9FEKyp32JTxavE8fI0VOi6gRoTCKAcdKA4XuTDsygd1bUABizZVNFSxz0dj+4d8/uNjCEKReILP6biAykXZHgPIPXEWh1OFIZ6icj1SGnbzBx0XUPlvZBaDi6TF4ZSFeIrKzXGdSJujYtBMoeNTRbWonMo/hdEC0OJwKoqYi8o5HZYzhUTQonIqbxVGArh0ib/PP5QWh1MRkp0ppF1CKt8URhdQHMXh1u04wIpNe2jp7KGmooQFMydrk75ApKaoXAcHOnu4YPlL+jukcl5hzAKK0bodB1i89i16+r3BYyVOO8uuPkv/gAtUvEXlOnwn0EMxNdJOi6nkPs8cXrB/Vn+HVMboLKAErdi0J+zmD9DT79UmfQGLp6jcMWPnROnFZWsPjgssd67kMu8r+jukcpImgBDRVnnq6s/CFWtRue+PGUNDXRXTJ1SHLSArlT7udKzmQGcPExY1c8Hyl3SQWOWMwhgDiFFNRQkHLG72NRUlFmerQtEwsSE4PlD/ZL3lOR85bHw08HWgVQDQcLSbWmln36iv+ruEjsxh8do+AO0SUlmnLYAQC2ZOpsQZPrhX4rSzYObkLEWkcs1QReVChc4UCp0qql1CKpdoAggxe2oty64+i9qKEgSorSjRwTsVJtapomA9Uyi0S+iC5S9pt5DKKp0FpFScmvc1h00V7fH00Hmsc9B5Q20/2WIqdaaQSot4ZgFpAlAqSc37mml8rZFeb2/UcwJF5RqOdlvuP7yo/xbeOOkyXl10SQYiViOZJgClMiy0VSAi+IJbYB9X3e/heXeLZb+rx9iw4aPFaE0hlRxdB6BUhoXWFIr2oarNYY/6B+cQn9YUUhmnCUCpFIs2U+ikURXMOKXOcrOZUFpTSGWKJgClUsxqppBDHHR7umPafxjCawrpTCGVLroQTKkUsyoqZzVTqNdm4zvjxrJ43NhBM4UOmzK2Fs2jpqedliJdQKbSQweBlcqA+ifro244EzDc9pM6U0jFQgeBlcoxsawgHm77Sa0ppFJNE4BSGRDrCuLA9pNWA8WBmkK/L5rHtCMvsHjtW5oEVFJ0DECpDIgcF4i2VgD8xeQC/w0tKicCArjEX1OIflixqUjHBFTCdAxAqSyIZfVwwLCbz9g6aPGN1QVkCtAxAKVyXuQ+A0NtSG+1+QzAGPnYv/kMRheQqYRoC0CpHBFt+8lI1f0eNrtbLJ9z+yq5tvRnOlOogGkLQKk8FOtAcZvDHnVCaa208/ueq3AvOY3Ge+/W1oAakg4CK5UjYh0oPsnAzDoXbXYZNC4QOlB8Z/8jLPm1B/gXHRdQlrQFoFQOCS0q94MLf2BdUsJRRKvDNmxJCa0ppIajYwBK5bCUbT6jM4UKhu4HoNQIFU9JiWibzywxt3LhVdotNFLpILBSI1SsJSWaRlcMuvnD8W6hf139P1pSQmUuAYjI6SLyqIisEZF/ztT7KjWSxDpTqNVh5/K6mqglJf5cdL2WlFDJJQAR+bmIHBSRtyOOXy4ie0Rkr4gsAjDGvGuMuR2YA8TUPFFKhYtcQGaTKH/CAwPEVgPFIgR3H1vuXMll3ld0oLhAJTUGICIXAR8DTxljpgwcswN/Ai4D3MA24HpjzG4RmQUsAh4yxvxiuOvrGIBSQ4unpES0BWRaUmJkydgYgDFmC/BhxOHzgL3GmH3GmD5gFfClgfM3GGM+DdyQzPsqpfziKSnR6rBbdglpSYnClfQsIBEZD/wmpAVwDXC5MeaWgcc3AucDa4CrgVHALmPMw1GudytwK8App5wy7YMPPkgqPqUKTSwlJUJnCllx+yr5bP9P8BlDTUUJC2ZO1hZBnsj2LCCxOGaMMS8bY+YZY26LdvMfOPExY8x0Y8z0cePGpSE8pUa2WAaKe202lo0ZzQyXDhQXsnSUgnADdSGPXYB15SqlVMpFlpSItm6gy26jS/yf13TvgcKUji4gB/5B4EuBA/gHgb9qjHkn3mvrILBSyYu1yigMP1DssnVAuQsuXQL1c1IdqkqBjHUBicgzwB+AySLiFpGbjTEeYC6wCXgXWJ3IzV8plRqxrh0A/0CxVbdQYKAYDHTtp2ftXLZt+Gkao1aZoKUglCoAsdYUwhj/QoEBQw0WHzCVbJu9RbuFcozWAlJKDSkV6wd8BlqppJp2Dso49p+zgHNn3ZaOcFUcsj0LSCmV41KxfgD8s4VsAlUcYsobd2m3UJ7RFoBSCohv/cAXPu4eVGgOoIVKvlL8M1o6e3T9QJZoC0ApFbdY1w98f8yYqIXmqo1/S0pdP5AfdEtIpRQQ+/qBI3YbHw0UodP1A/lNu4CUUpZSsX7A7avkwr4HEdAuoQzRLiClVNLiXT9g1SVUI+1sLZrHX0Z9lWe7v65F5nKMJgCllCWrmUIVoyqsTxaJukm9y9Ye3H9gqTzGzubHMvQdqOFoF5BSKmaxrh8o93gpNYY2h33QRvUeY8Pfb3zCAAAPYklEQVQmRvceSBPtAlJKpUVkqyCaLrst6o5kDvHp3gM5QlsASqmExTNQbDP+eUWRLYIO3wn020qool0LzaWAtgCUUhkRz0CxL8o4wRj5mCoOESg051n/Tdi1Oo1RqwBNAEqphMU1UByi12ajabT/PInoSXJ4e+l+bkk6wlURtAtIKZVSMReaM4Zqj9dyoNgABtGB4gRoF5BSKmsiWwU2iX6biTZQLKADxRmgLQClVFrFU3o62kCx21fJtaU/49VFl6Q52vynLQClVM6Ip/R0tIHiGmnn2e6v47u7nLbGSVp2OkW0BaCUyrhYp48Gagz5DGHlp3tMEW9Pu1c3oLGgLQClVE6Ldfpom8M+6OYPUCJ9TH1zETRWwI+m6LTRBGkCUEplXKwDxSd4TdS9Bxz40E3qk6NdQEqprLMaKHaIAxGh39cfPKab1A9Pu4CUUnnFaqD4hKITwm7+4F9A9p1xY613I6ODFZv2ZDr0vKYtAKVUTqp/sj7qrmQBoS2CDt8J9FBMjbRzUMax/5wFBTlIrC0ApVTeqyqrGvacQEmJY8bOidIb3HugikNMeeMuHRcYhiYApVROinWmUKvDTkNdFdMnVId1CwVmCvkaK3AvOY3Ge+/W1cQRNAEopXJSrDOFDMJfnXbrvQfQvQeGomMASqm8oCUlYqNjAEqpESc1JSU6aOnsyVTIOU9bAEqpvBVvSQm3r5LP9v8EnzHUVJSwYObkEbduQFsASqmCEE9JiW5TxIu+s3nF+U3+MuqrPNv99YIfE8hYAhCRiSLyuIisydR7KqVGtnhKSlxW5+L/TNrD104p4rkTSnHZ2lkqj7Gz+bEMR507kuoCEpGfA1cAB40xU0KOXw40AXZgpTFmechza4wx18Ryfe0CUkrFI5GSEh5jwyFmxGxIn8kuoCeAyyPe3A48DHwBOAO4XkTOSPJ9lFJqWPGUlAjsSeyQwi0q50jmxcaYLSIyPuLwecBeY8w+ABFZBXwJ2J3MeymlVCwaJjbQMLEh+Lj+yXrL81oddurH14VNFS3hGDVv3Me6uitG3OCwlXSMAdQC+0Meu4FaERkrIo8CU0VkcbQXi8itIrJdRLYfOnQoDeEppQpJ1JISUaaKFlJRuXQkALE4ZowxHcaY240xpxljlkV7sTHmMWPMdGPM9HHjxqUhPKVUIYllplBol9BhU8az3V8viM1mkuoCisIN1IU8dgEtaXgfpZQaVqA7qOnNJtqOtkWtMNrqsHOZq4a/OuwD3UIlNHTth43z/Cfk+eCwlXS0ALYBnxSRCSJSBFwHbEjD+yilVEwaJjaw+ZrN7Lpp15AriNucjsHdQv098OLSDEabOUklABF5BvgDMFlE3CJyszHGA8wFNgHvAquNMe8kH6pSSiXPskvIABLeex3aLWS63BmKLrOSnQV0fZTjvwV+m8y1lVIqHSK7hKrKqqKWkwjMFBrnMXx7X3PY7KKRIB1jAEopldMip4pGrSkkggEOOoXG1xqDrx0ptBaQUqrgxTRTyNtL08sLR9SsIE0ASqmCF1hBXO78BMZAtAo5bTZG1GphTQBKKYU/CWz96ovc+6nnsXlHW55T5fEeXy08AqqIagJQSqkQs6fWsuzihYO6hIp9PuYf7gRGzmphHQRWSqkIwZlCLy+kzTZ4a8kWM3ZE7CymCUAppSw0TGyg4eOj9KydSwnHAGguK+XHoytoczg40ftDmvfl96wg7QJSSqlo6ufw9jn3cMBUsrG0jLsrx9LmdICAcRym8ZU7aX44f+sFaQJQSqkhnDvrNrbN3sJ3x07imM1itfAob97ODNIEoJRSw5g9tRYcnZbPtTnseTszSBOAUkrFINq+AlUeL5CfM4M0ASilVAysVguHTg3Nx5lBOgtIKaViEJwa+t/LaOvrDJsa2m2KuM8zh5qKkixHGR9NAEopFaNAEbltG35KzRv3UU0PblPJfZ45vGD/LMtmTs52iHHRBKCUUnE6d9ZtrKu7gu+/8jTdZRuxOX/D2KI/4iz/Nv5t0fODJgCllEqAs3wnVP4Sm7cXgK7+g3lXMloHgZVSKgFNbzbRO3DzD+j19vLAS3fSeO/deTElVBOAUkoloO1om+XxQw7hzv5H2PrrR3I+CWgCUEqpBAy1LqBU+vgWq3J+XYAmAKWUSsBw6wJqpCPn1wXoILBSSiUglpLRub4uQBOAUkolKFAy2rP+mzhCBoS7TRE/5joW5Pi6AE0ASimVjPo5OIDu55ZQ3NPG/y39BA+NraDHsZFdf9mGs3x+zk4L1TEApZRKVv0cShf+L8/9w9P8pKacHkcPYGg92krja40072vOdoSWNAEopVSKRFsb0PRmU5YiGpomAKWUSpFoawNaPm7lguUv5dy6AE0ASimVItHWBpzs8THtyAssXvtWTiUBTQBKKZUi0dYGfPvwhyx3ruQy7ys5tThME4BSSqVIw8QGGj/dSLXXIMZQ3e+hsf1DGo52Uyp93OlYnVOLw3QaqFJKpVDDxAYanroBMIOeq5GOnFocpi0ApZRKtXKX5eFWxubU4rCMJQARmSgij4vImky9p1JKZcWlS8AZ/km/h1G0TLuT2VNzZ8OYmBKAiPxcRA6KyNsRxy8XkT0isldEFg11DWPMPmPMzckEq5RSeaF+Dlz5IJTXAQLldZRc/RDnzrot25GFibUF8ARweegBEbEDDwNfAM4ArheRM0TkLBH5TcS/T6Q0aqWUynX1c+COt2n+h6eZUVdD/Y57mbFmRk6tCo5pENgYs0VExkccPg/Ya4zZByAiq4AvGWOWAVekMkillMpHzfuaaXytMbg6OFAaAnJj28hkxgBqgf0hj90MsRuyiIwVkUeBqSKyeIjzbhWR7SKy/dChQ0mEp5RS2ZXrpSGSmQYqFscGz3sKPGFMB3D7cBc1xjwGPAYwffr0qNdTSqlcF600RLTjmZZMC8AN1IU8dgEtyYWjlFIjR9RtI70+2LU6w9EMlkwC2AZ8UkQmiEgRcB2wITVhKaVU/ou6bWRHB2ycl/UkEOs00GeAPwCTRcQtIjcbYzzAXGAT8C6w2hjzTvpCVUqp/DJUaQj6e+DFpVmNT4zJ3W726dOnm+3bt2c7DKWUSk5jBdZDpAKNnSl9KxF5wxgzPZZztRSEUkqlW5TSEFGPZ4gmAKWUSjeL0hA4S/zHs0gTgFJKpZtFaQiufNB/PIu0HLRSSmVC/Zys3/AjaQtAKaUKlCYApZQqUJoAlFKqQGkCUEqpAqUJQCmlCpQmAKWUKlCaAJRSqkBpAlBKqQxr3tfMjDUzqH+yPqvbROpCMKWUyqBc2iZSWwBKKZUh63YcYPHvfpgz20RqAlBKqQxYt+MAi9e+hc9+2PL5bGwTqQlAKaUyYMWmPfT0ezH9FZbPe/vKuWD5S6zbcSBjMWkCUEqpDGjp7AHg2KGZGJ8z7Dnjc3Ls0EwOdPaweO1bGUsCmgCUUioDair8+wF4jkylt/VqfH0VYOAT/T6WHWrl5d7/ZJZtKz39XlZs2pORmDQBKKVUBiyYOZkSpx3wJ4FL37uQ199r40W3myu7u3HZ2lnuXMks29ZgayHdNAEopVQGzJ5ay7Krz6K2ogQBvlP0S0qlL+ycUunjAeej/KX4BvjRFNi1Oq0x6ToApZTKkNlTa5k9tdb/oPEGy3Mc4vN/0bUfNs7zf52mjWS0BaCUUtkQy4bw/T3w4tK0haAJQCmlssFqo3grXe60haAJQCmlsiFyo3ixW58XS0shQToGoJRS2RK6Ufyu1f4+//6QGUDOEn9LIU20BaCUUrkgskVQXud/nKYBYNAWgFJK5Y7QFkEGaAtAKaUKlCYApZQqUJoAlFKqQGkCUEqpAqUJQCmlCpQYY7IdQ1Qicgj4wOKpcqArysuiPRd5fKjHVl8H/lsJtMcQfjwxR3t+uGMac2pitjqe7zEPF2vosUTi1piHjjn060zHfKoxZlxMVzfG5N0/4LF4n4s8PtRjq69D/rs91TFHe364YxpzamJO5Pch12MeLtZk49aYY/9dyUbMsf7L1y6gjQk8F3l8qMdWXw/1nrEY7vVWzw93TGOO7f1jeT7e34fQr3Mx5sjHVrEmE7fGbP041b8ficYck5zuAspFIrLdGDM923HEQ2POjHyMGfIzbo05NfK1BZBNj2U7gARozJmRjzFDfsatMaeAtgCUUqpAaQtAKaUKlCYApZQqUJoAlFKqQGkCSAERmSgij4vImojjZSLyhohcka3YorGKWUROF5FHRWSNiPxzNuOzEiXm2SLyMxFZLyIzshmflSgxW/6+5IooMZeJyJMDP2vr3cxzgIicISKrReQ/ROSabMcTKxE5RUQ2iMjPRWRRxt44mUUEI+Ef8HPgIPB2xPHLgT3AXmBRjNdaE/F4KbAQuCJfYh44ZgMez7OYR+dhzIOO5WrMwI3AlQNfP5vquFMVO/CvwGcGvt6QjjjTFPfngdsGvn4qE3EbY3RDGOAJ4CHgqcABEbEDDwOXAW5gm4hsAOzAsojXf80YczDyoiLyeWA3UJwvMQ9cZxawaOD6eRHzgLsGrpVKT5DemNPhCdITswt4a+Brb4pjDniCJGMH/gu4e+D3eGya4oz0BMnHvQP4rohci/97yIiCTwDGmC0iMj7i8HnAXmPMPgARWQV8yRizDIi1O+dioAw4A+gRkd8aY3w5HjPGmA3ABhFpBn6RingHrpuWmEVEgOXAc8aYN1MVL6T355wuaYzZjT8J7CRNXccpjP0bAzfgtemIM1Iq4haRfwPuHrjWGuA/0xu1n44BWKsF9oc8dg8csyQiY0XkUWCqiCwGMMZ81xjzLfw30Z+l6uY/hKRjFpHPiciDIvJT4LdpjdYv6ZiBb+JvPl8jIrenLdLjUvFztvo+0ikVP+e1wJdF5D9IvvRFPOKNfbyIPIb/0/iKNMc2lLjiBp4H5g383N9PY1xhCr4FEIVYHIu6Ys4Y0wFY3nyMMU+kKKbhJB2zMeZl4OWURjW0VMT8IPBgiuMaSipijvr7kiapiPko8E8pjisW8cb+PnBr2qKJXbxxvw1kfNBaWwDW3EBdyGMX0JKlWGKlMWeGxpxZ+Rp7XsStCcDaNuCTIjJBRIqA64ANWY5pOBpzZmjMmZWvsedH3JmabpSr/4BngFagH3/Wvnng+BeBPwF/Ab6b7Tg1Zo15pMac77Hna9zGGC0Gp5RShUq7gJRSqkBpAlBKqQKlCUAppQqUJgCllCpQmgCUUqpAaQJQSqkCpQlAKaUKlCYApZQqUJoAlFKqQP1/NnqJUNpkmDgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's plot it\n",
    "for freq in freqs:\n",
    "    sd = s_distributions[freq]\n",
    "    label = \"%.2f Hz\" % freq\n",
    "    plt.loglog(sd.size, sd.ps, 'o', label=label)\n",
    "plt.legend(numpoints=1); "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ok, ok, but this not so general enough!\n",
    "#### I do not remember the frequencies, and I am interested in the files ending with S, T, u, v. Can we do something more general? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# It would be nice to have something like...\n",
    "# d_T = Dist(\"Duration\", mainDir)\n",
    "# d_S = Dist(\"Size\", mainDir)\n",
    "# d_E = Dist(\"Energy\", mainDir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### It is better to stop for a second. This starts to be pretty complex. I cannot put everything here\n",
    "\n",
    "Let's make a separated script for the class. We will learn how to use it...\n",
    "\n",
    "[Work on cplot_distributions.py]\n",
    "\n",
    "Done? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cplot_distributions as d1\n",
    "mainDir = \"/home/gf/src/Python/Python-in-the-lab/Bk\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_size = d1.Dist('Size', mainDir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAETCAYAAAAVhSD9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8lPW1+PHPmUlCFjQBEiELCBRLVURQ1Fq91hVooxStdWnrbaMVvK0V299VobUY8Vpo9Vr1qnXBovZagYsK2OB2bZXiCogFl6uluBC2EHZIQpKZ8/vjmQmTyTPJTDKTmUnO+/XKK5nneeaZM2hy5rudr6gqxhhjTLx4kh2AMcaYnsUSizHGmLiyxGKMMSauLLEYY4yJK0ssxhhj4soSizHGmLiyxGKMMSauLLEYY4yJq4xkB9BVIjIEuA+oBT5R1TlJDskYY3q1lGyxiMgfRKRGRN4POz5RRD4WkfUiMj1w+MtAlapeCRzT7cEaY4xpRVKxpIuInAHsB55Q1VGBY17gE+A8oBpYCVwObAMWAQr8UVXndXT/wsJCHTp0aGKCN8aYHmj16tW1qloUzbUp2RWmqstFZGjY4ZOB9aq6AUBE5gPfApqAWwLPWQS4JhYRmQJMARgyZAirVq1KUPTGGNPziMjn0V6bkl1hEZQCG0MeVweOvQBcJyIPAp9FerKqPqyq41R1XFFRVEnXGGNMJ6RkiyUCcTmmqvo+cHF3B2OMMcZdOrVYqoHBIY/LgM1JisUYY0wE6dRiWQkcJSLDgE3AZcB3kxuSMSbdNDU1UV1dTUNDQ7JDSUnZ2dmUlZWRmZnZ6XukZGIRkaeAM4FCEanGGZx/VESuBV4EvMAfVPWDJIZpjElD1dXVHHbYYQwdOhQRtx723ktV2bFjB9XV1QwbNqzT90nJxKKql0c4vgxY1s3hGGN6kIaGBksqEYgIAwYMYPv27V26TzqNsRhjTFzEmlQufehNLn3ozQRFk1rikXAtsaSB3vQ/tTEm/VliSXGL12xizRe7efvTnZw25y8sXrMp2SEZ06sk4nfwhRdeYOTIkYwYMYI5c9zLGx48eJBLL72UESNGcMopp/DZZ58BsGPHDs466yz69u3LtddeG/E1zjzzzFYLwT/77DNGjRrV5dijYYklhS1es4kZz6yj0ecHYNPuemY8sy7pycVaUKa3SMTvoM/n4yc/+QnPP/88H374IU899RQffvhhm+seffRR+vXrx/r16/nZz37GTTfdBDiztm677TbuvPPOTseQaJZYUtgdL35MfZOv1bH6Jh93vPhxkiKyFpTpXRLxO/jOO+8wYsQIhg8fTlZWFpdddhlLlixpc92SJUv4wQ9+AMDFF1/MK6+8gqqSl5fH6aefTnZ2dqdj+NGPfsSYMWMYM2YMRUVF3HrrrZ2+l5uUnBWWLipeqABg3sQO615GFPzkv2DqqW3Obd5d7/qcSMcTLdKnN4DJY0uTEpMxiZSI38FNmzYxePChtd5lZWW8/fbb7V6XkZFBfn4+O3bsoLCwMOrX+t73vkdOTg4AjY2NeDxOW2Lu3LkAfP7550yYMIEf/vCHnX07rqzF0klVG6pYu30tq7atYvyi8VRtqIr5Hh19+i8pyHF9XqTjiRDa7ZWKLShjEikRv4NuFeXdZmJFe117nnzySd577z3ee+89li1rvVKjoaGB73znO9x3330ceeSRMd23I5ZYYlDxQgUVL1RQtaGKyjcqafQ3ArDlwBYq36iMKblE03d7w4SR5GR6Wz0vJ9NLdoanU2McsY6NhCe+TSnWgjIm0SL9Dt4wYWSn71lWVsbGjYfq6VZXV1NSUtLudc3NzezZs4f+/ft3+nXDXXPNNVx00UWce+65cbtnkCWWKIW2UH6x4hc0+FqXg2jwNXDPu/dEfb9oPv1PHlvK7IuOI8vr/GcqLcjh2yeWsnFXfZtWTkdJIzxJ3Lx4HV/+5fMMnV7l2lpyS3yRPitF8+nNBvxNOnL7HZx90XFd6vo96aST+Mc//sGnn35KY2Mj8+fPZ9KkSW2umzRpEo8//jgAixYt4uyzz47bos7777+fffv2MX369I4v7gQbY4lCeAvFr37X67Ye2Br1PaPtu508tpSn3vkCgMtPHuLayln1+U7WfLGbRp+f0+b8hRuylzA5fz1UOC0otyTx32990fIabmMlbolPcUpMhzbQo/n0FkxqLfFNGGljMiZthP4Ouo2FxiojI4P77ruPCRMm4PP5uPLKKzn22GMBmDlzJuPGjWPSpElcddVVXHHFFYwYMYL+/fszf/78lnsMHTqUvXv30tjYyOLFi3nppZc45pjoN9C98847yczMZMyYMYDTernmmmu6/N6CUnIHyUQbN26cxrLR1/hF49lyYEuH1xXnFfPSxS9Fdc9IXUulBTm8Pv3sVseCn/Srd9W7Pif0j/0kzwpuylxAiexA8gfDOTM5bVlhxG4sQp73i8z5DJKdkD+YadsvYIn/dNdrs7weGn1+SgtyuGHCyHZ/6YJJLTRJ5WR6u/ypz5jO+uijjzj66KOTHUZKc/s3EpHVqjoumudbV1gUommJZHuzmXbCtKjvefcx/+D1Pj9lQ5/vsiLrOiZ5VnQ4fhKplROaVOZkzqVUdjjdVns2wnPXMW7vy+3GEnzeINnpHNizkTuyHuLdPlNaxQdQmrGXZ/Nu55Rh/VsSYHsTEGIZ8LfuMmN6BkssURiUN6jd88V5xVR+rZLy4eXR3XDtQk5adwulsgOPQJmnlt9kPcqtwz5wHT8Jam8sI+PwNXwwfClfHTaQ8WUlVOXlOiea6pmR9T/thnNjxkJypbHVsSx89Jf9LfHNyZzL7VnzeDnz/3FM4zru2/avrFz6UJsuthXPPkDdrSVoZT7bZx0VMalt3l3fKpHY+hhjeg4bY4nCtBOmUflGZasB+2xvNsV5xQzIGRD7OpZXZkFT69ZHDgc57fMHaPTdC7Qe9wh2Mbl1K4GScfh7ZBc/Q01gjvqWzAwqC53ZI+UH6hhILTnSRL26769QIrUdhpwrjXxXXkYCzaMifw193/0V5/muYilOl9kkzwpmyVxytbHlmtmZc9EmWBrWrZafk9ky7jLm1pc40NhMk0/bvHfrLjMm/VhiiUKwJTLz9Zk0+hspzitm2gnTWPTJos7dcE91y49Vebnc06+ArRleBjX7yKhZQ/PescChLqPgH9fJY0vhtd9y+7ZT2HXYBrKPeAHJ2IMHxR82W6TB4+GefgWUH6hjl+bxRp+fkq972ayFzM36Ps2jLmbhymoafX5qpIhBdFwmO3w+Sg4HuSvzQe7mATZrIbnS0KblkyuN3JixkKWNhxJLpkdaJZLd9U1tXiv43uM5aGqM6R6WWKJUPry8JZEEWyidTiz5ZbBnI1V5uVQW9qchpKWRXfwMDdCSXFqNq6xdyOR9T5HRfyGVhQM46HH+1PsjTATemuHloHrpSwNZNINAmdRSKQ/B8GP5xzZnVe+gU39N3dM/aZMUopEhTjdYmdQSaR5IqdSyoc93W5LaEv9p7KprnUwmeVZwY8ZCSqSWzVrIb5svYenu09m+76DNJjPJNy/QzV0R+0Lo3sjGWLpg3sR5nSvncs5MyMzhnn4FLUklSDxN9Cl6seVxy7jK2oXw3HXgO8i9/Qpakkp7ipqVA+SQJc2tTzTVO91xARWbn+dXXEG1vxC/Cjv8fWnUQ585qvJyGV9Wwuihg1uP34SJNMVehJaxmkp5iDMa/trqfHDyQJmnttWYzrc8K9pMrT7nzldtgN+YFGeJJRlGXwIX3MvWDK/racncDbReI1L75+taxmUiPS+U+jP5vOZyCtjvfkGgO263523Wbl/LS8NfZuKQYr7c5985sfFh/r1pCpt0AH/Oy+WWwgFsycxARVrGb/6jf4Frsulw9npTPbdkPs6KrOtaZpxVZj7h2oX2n5kPtpqVVt/k45+1B2yA33SvtQuheiV8vgJ+N8p53EVdKZv/8ssvc+KJJ3Lcccdx4okn8pe//MX1+cksm29dYTHoSrHJNkZfwqB3f80WfG1OaVNByxqRyWNLqdpQxT1H9GNrRiGDmn3k+/3s9rZNLh5VFKel8nnNRTTvHcvmrELK3AbnvVmcXTeN+/v50MDCT8ncQ07x09QDqz3nsTK7jgcKXuOgtI6xweNhweGHtTRRtmRmcHNhf+YM6Mcej4eiZuX6nbs4/8CBVq2Y8PGkabtyKT8QuQsttJttTuZcCJkEYAP8pluE9BQALVP4AecDYicEy+a//PLLlJWVcdJJJzFp0qQ2CxxDy+bPnz+fm266iQULFlBYWMhzzz1HSUkJ77//PhMmTGDTptT6kGUtliSapvlka+v+o2y/n5v31PCvpbezZNvNLav+Q1sM+0XIDPtrnO338+vtO1j5xWbm7BhGv33DAZyBem9YeW1PJqiPpw/bh4aPq3iaOaz4ZV6ffjaT89ez1SXxAW36vZo9HnZ7vagINZkeflE0kD/mDWw5HxxPCm/5VOXlRuxCC5UrjdyS8USrls55vte44+nlfPDr0617zCSGywzO8K7kWHW1bP7YsWNbaosde+yxNDQ0cPDgwZhisLL5PVBLuX36gsLtspv96gt8ine6wf7QuJmt22p4t+bdNiVkmj0e8pt95KBs9WZQ7PMxbecuTmlSHhs2lqlX/JnhD73JcKBy6q2w9lhYcq3zqSunPzTuB19TxC615uBCSWAQXtdWVYc8zfy+sB+XVtfSx+9zHU8KnbkWjf6yHxGna69Markj8yEOkEPBwf3856bvU/kf32dM+RSbSWbiJ2QGZ1THoxDPsvlPP/00Y8eOpU+fPq6vlayy+ZZYulmwmGWjv5Hx4mWa5rN00zYKG51PReEzxSLVJdvr9XCFdyBvFQ1h3pYaKBxGRfERAEwNv3j0JbDaKWbH7s+h3kkcg5p9bMls+79AhgYqqFZUMS3QYmpVdDNYNKwD+z11nHXkYParj0hDLy3JLac/jfX7yKQJES9o22QW3rLpIz76cCjR3Nj0AL98uomPmseyjxybSWa6LjCD0/V4J8WrbP4HH3zATTfdxEsvRS4j9eSTTzJunFOF5bPPPuP8889vOWdl89Pd2oVw2xFU3VFM5fKbDpXbFx+VmXW8lXnofxa3T/ZuBvUt4a2iIbHHEvJJa9qu3WT7Wycu0SyOaL6w5XH58HIqv1ZJlicLgGL1cqnmkR3evRbBPvyoSMQpY4OafRz0eHmoZDhZR57CsiPH8LWhg9tOCojitXKlkZs9j7Msazr/7PM9FtRdzYpnH7BBftN5gRmcrWTmOMc7KR5l86urq7nwwgt54okn+NKXvtSpOBJZNt9aLIkWMvh3T78SGsIXMvoauLt/P84PdAdFM+MrW8V1gWbo5IKI3UAhn8CCXVAtA+p9S5h2wrQ2pWnc1vCM3VDlLBj1NZKPUOfNoMkfsjYlilZNtt/PSfVw1pGD2eev4Y8i1KE0AYSMw7yRXcDbuVDj9bZ0F0bqPgvvLpulD/Pbqgwmj41vH7LpJYID9MGu5EBh184O3EPrsvmlpaXMnz+fP/3pT22uC5bNP/XUU1uVzd+9ezfl5eXMnj2b0047rVMxJLpsvrVYEi1k8C9S0qjJ8LZ8KhrU7D6e4VFFVClu9lGp/SgfXh7bOpqKKqio4qHCgRz0HIqj/EAdz23exk+8A3np4pcoH17esqFZe8qHlzO6aDTj6MMKHcxtp912qFWTV9xhUilWL9/iMF7q15d9OK2mPaI0hT2vweNh6WEZbMtoO+jvJrxhlCuN/Kjxv22A33Te6Eug7CQ48nT42ftdSirQumz+0UcfzSWXXNKqbP7SpUsBuOqqq9ixYwcjRozgrrvuapmWfN9997F+/Xpuu+22lgH4mpqamGK48847WbduXcvzH3zwwS69p3BWNj/RKgsIduSMLytxHdMozivmpS//CJZcS1W2l8qiAa1aNtl+P5W1Ow99Sj/y9E6tAK4KtjL8jRQ3OwP+5RkDeKgwMFYzcV7rawKla4LJBlq3isKPhT6OtNVAlieL0UWjmbelhvGyiS3SiYkBwMDmZp7fuJVMDnXlqbr3uKk6/wWCK//HlE+xcZdezMrmd8zK5qe6kEE+tzGNlnL7gU9F5YUnUHnGbxjY7HNaKE3NrZNK/uBOJ5VW2ylneKksKqLqW79pGatpb8vlWKsMTDthWptxmGxvNqV9D/1BjziVOQrbMjI4YVgZ5w0u5bncXKr9heyir+u1oSv/b2x6wMZdjEkwG2NJtHNmOmMsTfUtyeHu/gVszcho1SIIVT68nOrX7+KHn66hj//QH9+DHi99Yhw0DLYiNu3f1HY7ZVHuefeels3Jxi8aH3HL5ai3BAh5D9B+4c5OT2UOsTXDS2VRIbkDvo4v42/s137tjsPkSiPX63wuXHgsRz3/AMf+YkWXXt8Y01baJxYR8QC3AYcDq1T18SSH1FrY4F95xgCqsw51PUUSbEVM/Wxdy6BhnxgHDUOnNkcSuolZpA3NIh3vqAXTbuHOCFOZMyQDRfGpj+K8Ys4oO4Ml65e0SXihGj3QKK85D0LGYQDX5FIitTyb+StKDtaytXIEG0+4gZMmtZmkbYzppJRMLCLyB+B8oEZVR4UcnwjcA3iBuao6B/gWUArsBDq/aimRQteRVFTxVqSB8bAurreKhjC1Odv1XEfCu7UiCd3EbFDeINdxkY42OguKpqss9JqOWjUtM9COGNtyTUQug/6/698v4uyxMo9T5mYQ28lffTMrwZKLMXGSqmMsjwETQw+IiBe4H/gGcAxwuYgcA4wE3lTVnwP/1s1xtki1bXXvefeedj/lQ9vtlCONi8Sy5XKsWmaXDRzXMiutvWuK84qjvve2DG+b9TB+dcZbQuVII4NXz7GZYyaiaGZKmkNSssWiqstFZGjY4ZOB9aq6AUBE5uO0VjYCwY+yETvsRWQKMAVgyJBOLCxsR3Bb3eC+IXcf8w9O+vuvWs97DxHTJ/t5sY1tBEXqvgpyG9+J1IKIdXwlVFwLd+K+m2d7Qqco79DDuOLANtfrBrKTIw7ubFUaxmaOGdM5qdpicVOKk0SCqgPHngEmiMh/AcsjPVlVH1bVcao6rqioKG5BBbcLDu4bcuLelxm1+ua21VD3xzbPvEVg/UmsInVfZXmyom4dRLommdpUAsgr5tKRl3ZYCaDB4+HOov6cO7jMdV8ZmzlmIgmOVa7atorxi8ZTtaHrm311pWz+O++807L+5Pjjj+fZZ591ff7QoUOprT1U2fzVV19tVdIlkVKyxRKB25I7VdU64KruDibojhc/brUH/Y0ZC8kJrxjcVA/N9d26+5zbJ/tsb3ZMXUndJbxV05lJAdGMw6gcoCbz0G6dkQb4gzPHLn16FJPf+2/bNbAXizQFH+j0h66uls0fNWoUq1atIiMjgy1btnD88cdzwQUXkJGROn/O06nFUg0MDnlcBmxOUiwtWm0djDPjyFUXqqF2htsn+8qvVTIgZ0C3xhEvHa2jiWocxmWA/55+Ba6XlkgtCzw3o3Hc3MmkH7exyuAU/M7qatn83NzcliTS0NDgWsCyI9/85jdbWj35+fk8/nh8J9OmTorr2ErgKBEZBmwCLgO+m9yQnK2DN4Ukl83qvrFWXc4g3IuQJI7bJ/tU69pKhFjGYbYEBvjd1r4EZ46xZyPNS37q/LJ0sZyHSS+xTsGPRjzK5r/99ttceeWVfP755/zxj3+M2Fo566yz8AY2Bdy/fz9f+cpXAFi2bBkAq1evpqKigsmTJ3f6/bhJyRaLiDwFvAmMFJFqEblKVZuBa4EXgY+Ahar6QTLjBLhhwkhyMg/V3vpt8yXUaVara+o0i5kHvp02ffaxrrJPhvZidGutFfRxb5kg0qYGmdvMsQxfA3XPd76irUlPkcYqo52C7yYeZfNPOeUUPvjgA1auXMns2bNpaHD/EPXXv/6V9957j/fee69lD5ag2tparrjiCv70pz+Rn5/fmbcSUUomFlW9XFWLVTVTVctU9dHA8WWq+mVV/ZKq3p7sOMHZFnf2RceR5XX+KVcffh63yTVU+wt5LjePc8rKOGXYIJ4fvJzbX3syydH2HuGTEKafPD2qAf67+vWPWD8zp24LWpnP9llHWddYL5GIKfjxKJsfdPTRR5OXl8f7778fUww+n4/LLruMmTNnMmrUqI6fEKOUTCzpZvLYUsYOKeCUYf15ffrZzG/4KmdmVzCjaJAzWCzgydpNff78uMwoMe5iacVE2uClJsPT7swxAYr8NdQ/cy0rlz4U53dgUk2kscqudCmHls1vbGxk/vz5TJo0qc11wbL5QKuy+Z9++inNzc2AswPkxx9/zNChQ2OKYfr06YwePZrLLrus0++jPek0xpI2Sgpy2D3gRcTT1Oq4eJo6VXerK1K9S6s7hY45bdq/ybXKABDVzLEcDlKy+rcsHny+rXfp4dzGKrsitGy+z+fjyiuvbFU2f9y4cUyaNImrrrqKK664ghEjRtC/f3/mz58PwIoVK5gzZw6ZmZl4PB4eeOCBVtsVR+POO+/k2GOPZcyYMQDMmjXLNbl1lpXNj5Pgiu0FU09l8ZpN3Pz3ia4l3AVh7Q/WxvW1TfSCq6cv/vLFUW+5XNzUzEvVbScg+lU4pfEBHjvMilmmEyub3zErm58M88rbXRE/eWwpBf4+rue6Muhn4setiyOSSBu07dI8ns38FUcfXMfWyhHWNWZMgCWWWK1dCNUrIWx9w4Kpp7baDnjGmbciYbPDEl13y8QmfIC/uG+k9S/SZszloHo5TBoo89TiEaeY5ajVN1tyMQZLLLEJ2b8eOFSuxWWGUPnwcoqbrkDUGcaKx6Cf6br2BvjdZgAR2H0ydEryE7kDOUAOWdLc6lIrZpk+euMQQLTi8W9jiSUWIfvXt2iqd467KPCfQo4OT9m6W6a18O4xj8uYS4PHw2/6DaSA/a73GMhOjj64LlDM8pa0WbvUm2RnZ7Njxw5LLi5UlR07dpCd3f7U/I7YrLBYRCrL0s3lWkz8hLdeQmcArd62Grd5yZ7M3Zw3uIztGdJmxX5wSnKZOMUsZz7bDPzYZo6lkLKyMqqrq9m+fXuyQ0lJ2dnZlJWVdXxhOyyxxCK/zOn+cjvuYsHUU6l44fAEB2USJdLGZxDdlGQrZpmaMjMzGTZsWLLD6NGsKywW58yEzJzWxzJz2uy3YnqGSGMuVszSmPZZiyUWYfvXt2zi1U5hQlugmL7cNj7bst+9BWPFLI05xBZIdkZwDYt1bfQKwUWV8ybOY/yi8RG7x4Ky/X4qa3fyjf11bYpZAtTlFJN70/8lIlRjEsYWSBqTIK7dY2GiKWbpvyWf6plfspljpkeyxGJMDOJVzNK2QTY9mXWFGdMJwe6xiMUswwb5g91j4TPHAKr9hVya+wivTz87QdEa03XWFdbdOqgdZnqu+Mwc28Hm3XW2Yt/0GJZYuipC7TDTO8RSzDI4cyy8a2yX5vF6n5/ain3TY1hXWFcEa4eFlnnJzIEL7rXppL1MZ2eOnbv/IIK0qjtWp1nM1CmcfqGt2Depw7rCukuMtcNM7xDLzLEDZLcpZpkrjVzPfH6+4F2GTq/itDl/sRaMSSu2QLIrrHaYcRG+sDLSBmI1GR4uG5zP1oz+bRZWlkot6/t8n81ayG/3XsKMZxoBrAVj0oIllq6IsXZY6C6Tpmdpr5hle9sgb8nMaPkeWnMstJjlnMy50AR3vJhlicWkBesK64oYaoctXrOJNV/s5u1Pd1rXRi/T1ZljudLIjRkL2bS7jpNm/Mm6x0zKsxZLV0RZO2zxmk3MeGYdjT4/AJt21zPjmXWAdW30Bm1qjqmXLfhcr41Uc6xEalmRNY0SqbXuMZPybFZYPHRQO+y0OX9h0+76NsdLC3JsUVwv0tmZY+UH6vArreqO1WkW05t+xAtyBmOHFFj3qkk4mxWWYja7JJX2jpueL9qZY/f0K2iTVOBQ91ijz2fdqyblWFdYNygpyHFtsZQU5LhcbXqDaGeObcnwMnFwCVszvK4zxzb0+Z51jZmUY4klHjoon3/DhJHMeGYd9U2H+tVzMr3cMGFkoiMzKSy6mWMS88yxp975ArDZhyZ5rCusG0weW8rsi44jy+v8c5cW5DD7ouPsk2UvM2/ivIgbv8Vz5tiGzz617jGTVD2ixSIiecBy4BZV/XOy43EzeWypfZI0EdnMMdOTpGSLRUT+ICI1IvJ+2PGJIvKxiKwXkekhp24CrPqjSWvlw8sZXTSacQPH8dIP36O4b4SCliKoSEvXWGhByzJPbcteL3My53Ke7zVuX/iaVU023SolEwvwGDAx9ICIeIH7gW8AxwCXi8gxInIu8CGwrbuDNKarYu4eC9PRzLG7Mh/k7ax/s6rJplulZFeYqi4XkaFhh08G1qvqBgARmQ98C+gL5OEkm3oRWaaq/vB7isgUYArAkCFDEhe8MXESj5ljGeL8KpSJs1vlzGebAauabBIrJRNLBKVAaGGuauAUVb0WQER+CNS6JRUAVX0YeBicBZKJDdWdja2YWLWaObZ1TYRxl8gzx0LlSiPX63wufXoURz1/P7MG3GH/T5qESNWuMDcun9UO7Tiuqo+l6sC9MdFqt2vsjNvjsFtlLQs8N9umYiah0imxVAODQx6XAZuTFIsx3a7NbpXqjXhtpN0qofUA/41ND7Di2QcsuZi4StlaYYExlj+r6qjA4wzgE+AcYBOwEviuqn4Q673jXivMmG7UUnNsSw3jZRNbxH1aclCw5tg39te1GeAHqPYXckbj3fjxUFqQww0TRtoYjGkj7WuFichTwJvASBGpFpGrVLUZuBZ4EfgIWNiZpGJMj1FR5d49FqbB4+H2/v2ZOLjEtRUT3FRsRdZ1nLj3ZWY8s85aMKZLUnLwXlUvj3B8GbCsm8MxJmVFO3Nsn9fDPnE+R1ppGJNoKdsVlkjWFWZ6mmD3WHu7VYYrbmrmpeq2w5TV/kJOb7yHItnLds237jED9ICuMGNM50SsOeYi0gB/sDTM21n/Zt1jplMssRjTg7jNHDvM7zZTHysNYxLGEosxPUx4zbFfnjnbSsOYbpWSg/fGmNhEWlQJVhrGdD9rsRjTC4S2YoqJtLDS6RaL1D0WlCt9u3GZAAAcqklEQVSNXM987nh6OR/8+nTrHjNt2KwwY3qJ4Myxi798MZVvVNLgazh0MkIrxqOKQpsWjF9hsxa27P0yN+v7jCmfYi2YHsxmhRljIoqlNIw/ygH+YGmYc+581VowxhKLMb1RqwF+LW2ne+yQjgb4r2c+n9bus22RjQ3eG9NbRBzgr6hi2oaqtt1jLtob4A+WhrFtkY21WIwxbbrHPBGHXiMP8IvQZu2LDfD3TjEN3ovIV3G2DP4qUALkALXAx8BrwGJV3ZWAOOPKBu+NcRfPAf4d/r7Uk20D/D1E3AfvReQHIrIOeAO4HsgF/gG8DewCTgHmAptE5DERGdapyI0xKSEeA/z9Zb8N8PdSHSYWEfk7MAenqvCJQD9VPUNVv62q31fVb6rq0UB/4GrgCOADEbk0kYEbYxIrfIC/0BehNEyI0N0rxQb4e61oBu/nAQ+qarujeqq6B3gSeFJEjgcGxSE+Y0w3am+A/99jGOAfX2YD/L1Zh4lFVe+O9aaq+nfg752KyBiTksJLw3gUItW33JKZ0fK9vb1fTvR9wlcXv4su2Umt5wiKJt8Ooy/pjrdjEshmhRljohbaPfbrM+a4l+gP6wML7R4LlSuNXOH9XwbJTgQo8tdQ/8y1rFz6UOLegOkWHbZYRGQScJSq/mc3xGOMSSGxFLcsVi9b8LleG9z7JbxrLHyhZQ4HKVn9W875ZCSFh/Wx3SvTVDQtlp8A48MPishVIrJKRN4XkTtF5PD4h2eMSWVRD/C3UxomXAm1vLzvW1aiP41Fk1hGAYtCD4jI6cDDwGigAPg58FcR6RP3CI0x6aGiin8/K7q9X2b378f4shLXHSxDF1oGpyhbckkv0SSWAThrVkJdjdOb+jVVLQO+BgwHrotveMaYVDNv4ryIXWTh618ibYu8x+uxEv09WDSJ5QDQ0s0lIgKUA++o6ioAVX0LuA+wtSvG9HKt9n7pW+x+kcsA/y+KBri2YEqklgWemzn64DrrHksT0SSWlcA3Qh6firMY8qWw694ARsYpLmNMmmivBTPthGnuM8dcWIn+niOaxHIfcLWI3CwipwKVOP9r/Dnsuroo72eM6SXcSsMURPFnwkr0p7eoilCKyK3AzcGHwKuqenbYNT8GblDVlK8TZkUojeleweKW8ybOoyrKFfyoUtzsc13Br+p8ut2shfy2+RJe9n6d2RcdZyv4EyjuRShV9Rac2WHXARXAN10uuxJ4L9ogjTG9k5Xo7/mi7rpS1Y9U9X5VfTy8bpiIlAJ7gWfjHaAxpucJX8EvmtX6ApcS/ZEG+HOlkVsynrAB/hQSlx0kVXUTcHaHFxpjeqV4reD3B2aThdcg6y/7EdkPOHXIbmx6gJmL6uG1g0y+/p54vhUThWjK5i8RkbHR3lBEskXk5yJyTddCi/r1JovII4E421QIMMakvs6W6A8utDx+WNtWzJyMR5i06zG2zzoK1i5M9FswIaLpCvsCeEtE3haR60TkBBFp1dIRkZLAH/hHgS044y3vdjYoEfmDiNSIyPthxyeKyMcisl5EpgOo6mJVvRr4IbaOxpj0F+UKfmh/oWWG+PFIWHHLeeXOl0moaGeFfQln58jvAfk4PaB7gYNAPyATp0f0HeD3wB9V1d/poETOAPYDT6jqqMAxL/AJcB5QjbO+5nJV/TBw/j+BJ1W1w4Rms8KMST2hM8cAqjZURVWiP1ykrZJ3al+yPUqOHrAS/Z0Qy6ywWPe8z8LZ7/5UoBjIBnYA/wcsV9XPYw834msNBf4cklhOBSpVdULg8YzApXMCXy+r6v+2c78pwBSAIUOGnPj553EL1RiTIMFkc/GWDVTKLhok5O+VywB/uGy/n8ranZQfqEO19YL/evrw/gm3cdKkqfEPvAeKJbFENXgvIl/HmWKcA6wCfqeqjZ0PsVNKgY0hj6uBU4CfAucC+SIyQlUfdHuyqj6MUziTcePGRZ9NjTFJV17xGoS0YIrVy36/n33e9n+Vgwstg5uMhQqW6F88+Hxb/xJn0ezH8q842xOH/me5UUTOVtWahEXmEorLMVXVe4F7uzEOY0wSlA8vZ9EnTqH1eVtqeGjPdu7v50Ol/c+47W2VXEItkxYfw/bnrGssnqJpsdwAvAVchTOuciZwF3A38N2ERdZWNTA45HEZsLkbX98Y080iTlOuqGIqUBblOExHWyUHB/jf/2yXdY3FQTSzwkYAs1T1/1R1s6r+CfgxcFE377+yEjhKRIYFxnouA5Z24+sbY1JMq4WW2p+s8ClDMWyV3NI1dvc0W8HfRdEklj7ArrBjfwGygITUBRORp4A3gZEiUi0iV6lqM3At8CLwEbBQVT9IxOsbY1JTu3vB0Jcpu3MQdVonxeqNeJ9g91j4Kv5iajl511Jbwd9F0a68Dx8h2xf4nhPHWA69mOrlEY4vA5Yl4jWNMWku0D32VnDa8pYazvJvpDbDfYDfrXvsG/vrKJGdQMgK/mebgR/bAH8Moq0V9qqIrAwsXJwGnEPEXRWMMSY1XLY3O0IdMveNxsa4rOC/nvncvvA1Pvj16d0UdfrrcB2LiPwIGBP4Gg305dAM8o0404/XBL9UNeUH1G2BpDE9V3sLLVvqkMWw/sWvTnn+EqmlRorYeMINvXKAP67rWFR1btjNj8JJMmMD308FLgxeDkTu2DTGmG4WPk25ve6xoGAdsnv6FYRMU86l/MB28lffzMqPniPXt5tZA+5gwdRTu+NtpJWYVt5HvInIEcAJwPGq+psu3zDBrMViTO/S2Y3GQrvMQlsxzerBg5/NWsjcrO8zpnxKjx+DiftGXx1R1RpVfSEdkooxpndrtdGYtrPRWDvTlIMFLss8zgD/imcfsNljIeLSYkk31mIxxlS8UAFb1/HVXZ62K/gj1SFTRWhb4LLaX8iluY/w+vSeuy1Vt7dYjDEmLQ06jqk/e4fZZ8xq2Sq5WL0cFmkJv4hrif4S2cHm3XU2cyzAEosxptcL32jsB3tcpimHCe0a26V5/C3rOo4+uI6tlSOcvV96sbhsTWyMMT2GSx2ySF1jWzK8nFdWwja3mWPQK6clgyUWY0wvFbHAZUDoNOXPNq2OOEV5a4QCl4NXz4EdS6GiKo5RpwdLLMYYQ/uJ5rK92TxcUE9j6OBBOwUuyw/UMZCd+D9bweaZX+o1U5KDbIzFGGM6MPVn7zDr63NaDfBHsiXDy+ihg5kwuITn++b2yinJlliMMSYK4QP8hb7oZ44Fa47d8eLH3Rt0klhiMcaYCCKW6a+oIst/ZUwzx0qklgV1V0NlPvxuFKxdmIiQU4IlFmOM6YQXr76+1fqXSPXeg3u/jBk2mCuHZDlrX/ZshOeu67HJxRKLMcZ0Umj3WMSuMZwZY20WVjbVwyuzujHa7mOzwowxJkqxzxzTdmeO6Z7qjir4pyVrsRhjTBxMLShiFv3DZo65p43gzLFzy0qp2tDz1rlYi8UYY+KhoopyYFE0WyOLoEBNplC5fAbgdKv1FJZYjDGmCzrqHmtTOTlMgyh3/fVGyvcfgNGXJCLEbmddYcYYkwgVVfwl9x6Km67ASw6qzpCLm+1eof6Za3tM8Urbj8UYYxKs4oUKtm4bTHXTK5C5t8354qZmXqrezCYtZOXk5SlZ+sX2YzHGmBQzaOBGrt3XRJa/9fFsv59pu3YDUMyOHrE63xKLMcZ0k9Prs5iyO4eBzT5EleKmZiprd7bsRLlZB7B5dx3MS++BfBu8N8aYbjJrwB0AfCfrZv51/VpyAoP6VXm53N2vgK0ZGQxvvok/1+7k/N+NgnNmpuWAviUWY4xJsODMsUsfehOAt4qGsHvPAL6/fRXv5jVwa2F/DnqcNS81mR5uLeyP1O5gwpKfOn+k0yy5WFeYMcZ0kwVTT2XB1FMB+L8ReazsdwG/6D+sJakEBVfnZ/gaqHt+ZjJC7RJLLMYYkyST89cjGW1niQFszXD2fMmu39qdIcVF2neFiUge8ADQCLyqqk8mOSRjjGlXy6LKeeUM8Inr6vxBzT4ANvsHUNadwcVBSrZYROQPIlIjIu+HHZ8oIh+LyHoRmR44fBGwSFWvBiZ1e7DGGNMFl+3NbrOvS3AKcp1mMTfr+0mKrPNSMrEAjwETQw+IiBe4H/gGcAxwuYgcA5QBGwOX+boxRmOM6ZqKKqb+7B1mnzELLzmgcESTn5nbd3H8vlxm6hTGlE9JdpQxS8muMFVdLiJDww6fDKxX1Q0AIjIf+BZQjZNc3iN1E6UxxkRUPrycRZ8sYuu2wdT+8yv8/LB/0mfIi0jmUtb+cyWZ+dPSqkhlOv0hLuVQywSchFIKPAN8W0R+DzwX6ckiMkVEVonIqu3btyc2UmOM6YRBAzfyoyPu5vDihUjmHgC2HNjCzNemM+GRu5McXfRSssUSgdvGBqqqB4CKjp6sqg8DD4NTKyzOsRljTFzMP7yh9WZhQKMHmuURKv9jF2PKp6RkLbFQ6ZRYqoHBIY/LgM1JisUYY+IqOFNs9GPHuZ7fniHc2PQAM59tBn6c0sklnbrCVgJHicgwEckCLgOWJjkmY4yJK6/2dz0+qNlHrjRyPfNTvlBlSiYWEXkKeBMYKSLVInKVqjYD1wIvAh8BC1X1g2TGaYwx8fYfX7+RbG92q2OhFZBLZAebd9cnI7SopWRXmKpeHuH4MmBZN4djjDHdJjj763d/vYEar4dBzT6m7drdqgJySUFOMkPsUEomFmOM6c3Kh5dT/fpd/GDDe2Rrc8vxOs3iTv/F3DBhZBKj65glFmOMSUFvFQ0B4IoNH5HtP8B/5w3kv/rn05C5jA//uSal17ZYYjHGmBT1VtEQpjZn89Du7dzfz4dKA+CsbZmx3Kl6nIrJJSUH740xxhwy//AGNLApWJBKIzP++hsWr9mUpKgisxaLMcakoNAKyDu87mu6/d5dzHhmHUBKrWuxFosxxqSyiqqIa1sGNvs5z/dayq1rscRijDEpLtLalp/v2smczLmM2/tykiJzZ4nFGGNSXPnwciq/VskRTX5EleKmZiprd1J+oI5caWRG1v8kO8RWbIzFGGPSQPnwcr7xxHddWwMDqe32eNpjLRZjjEkTO7PcV9xLfmptXmyJxRhj0sTTg4/loMfb6thBjxfOmZmkiNxZYjHGmDTxVtEQHhs2lkYy+XNeLucOLmXckaWM/2QuVRuqkh1eCxtjMcaYNBFc2/LQ707m4YL6lg3BthzYQuUblUBqrMS3FosxxqQZt10mG3wN3PPuPckJKIwlFmOMSTORVuJvPbC1myNxZ4nFGGPSzKC+xe7HfX5Yu7Cbo2nLEosxxqSZaSdMc99lcscOeO66pCcXSyzGGJNmgivxBzb72qzEp6keXpmV1PhsVpgxxqSh8uHlfPOJ7yJuJ/dUd3c4rViLxRhj0tSOCCvxSfJKfEssxhiTptxW4pOZk/SV+JZYjDEmTU294hX6TH4QvH2cA/mD4YJ7YfQlSY3LxliMMSadjb4EVj/u/FyRGmVdrMVijDEmriyxGGOMiSvrCjPGmHSXIl1gQdZiMcYYE1eWWIwxxsRV2icWEZksIo+IyBIRGZ/seIwxJpkqXqig4oWKpMaQ1MQiIn8QkRoReT/s+EQR+VhE1ovI9PbuoaqLVfVq4IfApQkM1xhjUlrVhirWbl/Lqm2rGL9ofNJ2lUz24P1jwH3AE8EDIuIF7gfOA6qBlSKyFPACs8Oef6Wq1gR+vjnwPGOM6XUmPHI3WzL/iEojkNxdJZPaYlHV5cDOsMMnA+tVdYOqNgLzgW+p6jpVPT/sq0YcvwGeV9V3u/s9GGNMsi1es4lN8nRLUglK1q6SqTjGUgpsDHlcHTgWyU+Bc4GLReSaSBeJyBQRWSUiq7Zv3x6fSI0xJskWr9nEjGfWQcZu1/PJ2FUyFROLWxVo9304AVW9V1VPVNVrVPXBdq57WFXHqeq4oqKiuARqjDHJdseLH1Pf5EObCtwvaC5g8ZpN3RpTKiaWamBwyOMyYHOSYjHGmJS2eXc9AAe3T0D9ma3OqT+Tum3jmfHMum5NLqmYWFYCR4nIMBHJAi4DliY5JmOMSUklBc6eLM17x9Kw5SL8jQWgcESTn9nbt/BqwzzO873GHS9+3G0xJXu68VPAm8BIEakWkatUtRm4FngR+AhYqKofJDNOY4xJVTdMGElOprMnS/PesZzz6em8/elWXqmu5oK6Oso8tczJnMu4vS/DvO6ZHSaqEYcveqxx48bpqlWrkh2GMcbExeI1m7hx0VoafX5WZP2UMs+ONtc0q4cM8Tt7tpwzM+Y9W0RktaqOi+baZK9jMcYY00WTx5by1DtfAFCypW1SAZykArBnIzx3nfNzgjYES8UxFmOMMZ3kyR/c8UVN9fDKrMTFkLA7G2OM6X7nzHT2ve/InuqEhWBdYcYY0wMsmHpq4Cfne+MzPyaTJkS8oL62T8gvS1gs1mIxxpieZvQl/CPrK3yYdRxc+GDbFkxmjtOySRBrsRhjTA90bHG+80NwgH7JteA72OlZYbGwxGKMMT3d6Etg9ePOz92wjbElFmOM6Ym6IYFEYmMsxhhj4spaLMYY0xt0YwvGWizGGGPiyhKLMcaYuLLEYowxJq4ssRhjjIkrSyzGGGPiyhKLMcaYuLLEYowxJq4ssRhjjIkrSyzGGGPiqlfueS8i24HPww7nA3vaeVqk8+HH23vs9nPweyFQG0X40cbV3vmOjlnMvTdmt+OR4nSLNfR8Z+Luzt9Dt2MWc+TzR6pqUVR3V1X7cpLrw505H368vcduP4d8X9VdcXd0zGLuvTG7HY8UZ4RYQ6+NOe7u/D2McMxijuJ8R1/WFXbIc508H368vcduP3f0uh3pTNwdHbOYo3v9js6nY8xuxyPFGfo4lWMOf+wWa1fi7k0xR6VXdoWlIhFZparjkh1HLCzm7pGOMUN6xm0xx4e1WFLHw8kOoBMs5u6RjjFDesZtMceBtViMMcbElbVYjDHGxJUlFmOMMXFlicUYY0xcWWIxxhgTV5ZYUpSIDBeRR0VkUdjxPBFZLSLnJyu29rjFLSJHi8iDIrJIRP4tmfG5iRDzZBF5RESWiMj4ZMbnJkLMrv/PpIoIMeeJyOOBf+vvJTO+9ojIMSKyUER+LyIXJzueaIjIEBFZKiJ/EJHp3friXVldaV8RV63+AagB3g87PhH4GFgPTI/yXovCHs8CbgLOT6e4A8c8wKNpFnO/NIy5zbFUjRm4Argg8POCeMcdr9iB/wf8S+DnpYmIMwExnwtMDfz8RKJjDv3KaD/tmE56DLgPeCJ4QES8wP3AeUA1sFJElgJeYHbY869U1Zrwm4rIucCHQHZiwk5M3IH7TAKmB+6fFjEH3By4Vzw9RmJjToTHSEzMZcC6wM++OMcc9BhdjB34I3BL4P/jAQmKM94xrwF+KSKX4sTfbSyxJICqLheRoWGHTwbWq+oGABGZD3xLVWcD0XZrnQXkAccA9SKyTFX98Yk6oXGjqkuBpSJSBfwpPhEnLmYREWAO8LyqvhuveCGx/86JksCYq3GSy3skqGs+jrH/JPDH/ZlExBkqHjGLyL8DtwTutQiYl9ioD7Exlu5TCmwMeVwdOOZKRAaIyIPAWBGZAaCqv1TV63H+MD8Sz6TSji7HLSJnisi9IvIQsCyh0Tq6HDPwU5yuhItF5JqERXpIPP6d3d5HIsXj3/kZ4Nsi8nu6XmcsFrHGPlREHsZpQdyR4NgiiSlm4AXgusC/+WcJjKsNa7F0H3E5FrHsgaruAFz/oKnqY3GKKRpdjltVXwVejWtU7YtHzPcC98Y5rvbEI+aI/88kSDxiPgBUxDmuaMQa+2fAlIRFE51YY34fSMpEA2uxdJ9qYHDI4zJgc5JiiUU6xm0xd490jDkoHWNPm5gtsXSflcBRIjJMRLKAy4ClSY4pGukYt8XcPdIx5qB0jD19Yu7OKWi95Qt4CtgCNOF8yrgqcPybwCfAP4FfJjvOnhC3xWwx98TY0zHm0C+rbmyMMSaurCvMGGNMXFliMcYYE1eWWIwxxsSVJRZjjDFxZYnFGGNMXFliMcYYE1eWWEyPFthXZbmI1IhIvYh8LiKLRWRiyDU/FBF1KfqX1Lg6eP63RWSbiOQGHn8sIn9xue73gff2Y5dzNcG9UQK1sFREftTB64qIrBGRG6J7p6Y3ssRieiwRuQ54FvgHcBVQDvxH4PTZIZdWAafiLEhLpbgiPT8D+DVwh6rWBQ4vB74aWJEd6gygLvA99B5HA0WB50VNnYVvs4BfiEj/WJ5reg9bIGl6LBH5Alitqhe6nPNo91SHbqOrcYnIt3EqXBer6s7AsStwKu+epqpvBI4V4mwWdS9wiaqWhNxjKvAgMEZV/x5orX0KXK2qczt4fS/OavDfqepvo3vXpjexFovpyfoDW91OhP7xDu8KE5HHAo/dvs4Med7x4mz9uivQnfW6iPxLvOJqx4+AF4JJJeC1wPfQlsm/AAeA/wKKRWREyLkzgN0c2mQrKEtE7gp0k9WJyJ/DuwhV1Qf8TyAOY9qwxGJ6sneAH4jIDSLy5RiedxtO11jo1+s4XUpfAIjICcAbOEniauDbwA7gf0XkxATFhYj0Ac4E/hZ6XFW/AD6ndWI5A3hDVf+J08IITzorXBLZDOAonFL2PwFOBF4Skcyw65bjFEQcHkv8pnewxGJ6smtw9gb/LfCxiNSKyFMiMr69J6nqP1X1reAXcDpOcvmeBnbvw9ns6QvgbFVdpKrLgAuBDcCvEhFXwBicran/7nJuOXCaiAR/r8/gUAL6W+AxIjIMp/z6a23uAPtwdiWsUtV5OPt5HAX8a9h1awLfvxpFzKaXscRieixV/QQYC3wduB1n+9sLgRdF5OZo7iEiFwC/AW5S1cWBYzmBe/4P4BeRjMCAugD/S9hAeZzjCo6TbHc59xpwODBGRA4HjudQYlkRElfwu9vA/aLQVoyqvo7T2jk17Lrg65dgTBhLLKZHU1Wfqi5X1ZtV9VxgOM64wi0i0q+954rI8TiD5I+q6p0hp/oDXpyWSVPY17VAv5BWQ7zjyg58P+hyLpgozgBOA3zA24FjfwOGiUhZ4Px+4F2Xe2yLcCx8C9z6wPecdmI1vZQlFtOrqOpmYC7OttxHRbpORAbibKL0FhC+BmQ34McZFD/J7SvWGWfRxoUzjgPQJvmo6j9wpkyfEfhapaoNgdPvA7tCzr2hqs0u9x8Y4dimsGPBqca17cRqeinb8970WCIyWFU3upz6SuC768wsEckGluDMqPpO+B9gVT0gIn/D6Wp6N9Yk0tm4Av4v8H04zuSBcMuBc3CSQcsAv6qqiLwBXAKMAB6LcP+LRaQy+J5E5DScLXDfDLtuWOD7x+3EanopSyymJ3tfRP6KsxjxU5zxh2/iDJ4vDMykcnM3cALwQ+ArIhJ67kNV3Qv8HOeP+Isi8ihOS6Ew8Dyvqk5PQFyo6hci8jlwMvDfLpcsBy4FBuAsogz1N2B2yHVuDgMWi8hDOAsoZ+Ms5Hwi7LpTcLr+3ooUq+m9LLGYnuwmnD/Ys3A+wftwtnWdjpM8IvkKkAk86XLuLOBVVX1XRE4CbsFZgJiPM6D9Ls7Cw0TEFbQA+A5wncu54EwvxZkiHepvOBMMGnCmPLuZzaEWTR7wV+BaVW0Ku+58YGnIyn9jWtjKe2PSjIh8CacL6kxVXZGE1y/BmWo9QVVf6e7XN6nPEosxaUhEHsEp6XJ+El77d8DxqtphXTPTO9msMGPS06+AlcHqxt1sC86qfGNcWYvFGGNMXFmLxRhjTFxZYjHGGBNXlliMMcbElSUWY4wxcWWJxRhjTFz9fwsU4YlZ2qTCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "d_size.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Classes can redefine operations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Rectangle:\n",
    "    \"\"\"\n",
    "    This class defines operations on Rectangles\n",
    "    \"\"\"\n",
    "    def __init__(self, base, height):\n",
    "        self.base = base\n",
    "        self.height = height\n",
    "        self.area = base * height\n",
    "        self.p = self._get_perim()\n",
    "        self._diagonal = (base*base + height * height)**0.5 \n",
    "    \n",
    "    def _get_perim(self, half=False):\n",
    "        \"\"\"\n",
    "        Get the perimeter of the Rectange\n",
    "        Parameters:\n",
    "        ----------\n",
    "        half: Bool.\n",
    "             If True, get the semiperimeter\n",
    "        \"\"\"\n",
    "        p = self.base + self.height\n",
    "        if half:\n",
    "            return p\n",
    "        else:\n",
    "            return 2 * p\n",
    "        \n",
    "    def print_perimeter(self):\n",
    "        return \"The perimeter is %f\" % self.p\n",
    "    \n",
    "    def __add__(self, other):\n",
    "        base = self.base + other.base\n",
    "        height = self.height + other.height\n",
    "        return Rectangle(base,height)\n",
    "    \n",
    "    def __repr__(self):\n",
    "        return \"Rectangle of {0} by {1}\".format(self.base, self.height)\n",
    "\n",
    "p0 = Rectangle(3,2)\n",
    "p1 = Rectangle(1,1)\n",
    "p2 = Rectangle(2,0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print p0.base, p0.height, p0.p, p0.area\n",
    "p0.print_perimeter()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct a larger Rectange\n",
    "p3 = p0 + p1\n",
    "print p3.base, p3.height, p3.p, p3.area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p3"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}