utf-8''AnalysisDropcounter.ipynb

{
  "cells": [
    {
      "metadata": {
        "trusted": false
      },
      "cell_type": "code",
      "source": "# Import libraries\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport scipy.stats as sp\nimport scipy.special as sps\nimport urllib as ul\nimport requests\nfrom datetime import datetime, date, time",
      "execution_count": 1,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "# Rainfall intervalometer: Theory and data processing"
    },
    {
      "metadata": {
        "slideshow": {
          "slide_type": "slide"
        }
      },
      "cell_type": "markdown",
      "source": "The rainfall intervalometer is a very simple instrument. Basically, it registers the times when raindrops fall on a small drum. If you want to build this instrument, please have a look at: \nhttps://github.com/nvandegiesen/Intervalometer/wiki/Intervalometer.\nHere, we concentrate on the underlying theory and on how to convert arrival times to rainfall rates."
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "## Theory"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "You can skip the theory if you are not familiar with integral calculus or the Gamma function. We will end up with a simple polynomial that describes the relation between the number of drops measured per time interval and the rainfall rate. This will be sufficient to process measured date. \n\nBasic idea is to first calculate what the arrival rates are as a function of rainfall rate according to:\n\\begin{equation} \n\\rho_A^* (\\Lambda)   = \\frac{\\alpha N_0}{\\Lambda^ {(\\beta +1)}} \\Gamma({\\beta+1, D_{min} \\Lambda})\\ \\  \n\\label((Eq. 1)\n\\end{equation}\n\nIn which:\n\n\\begin{align}\n\\Gamma(a,x) &: Upper\\ incomplete\\ Gamma\\ function\\\\\n\\rho_A^* (\\Lambda) &: Arrival\\ rate\\ (m^{-2} s^{-1}) (measured)\\\\\n\\Lambda &: Constant\\ related\\ to\\ rainfall\\ rate\\ R\\ (mm/h)\\ according\\ to\\ \\Lambda =4.10 R^{-0.210}\\\\\nN_0 &: Constant,\\ here\\ taken\\ as\\ N_0=8000\\ (-)\\\\\n\\alpha &: Constant\\ associated\\ with\\ terminal\\ velocity,\\ here\\ taken\\ as\\ \\alpha=3.25\\ m\\ s^{-1}\\ mm^{-\\beta}\\\\\n\\beta &: Exponent\\ associated\\ with\\ terminal\\ velocity,\\ here\\ taken\\ as\\ \\beta=0.762\\\\\nD_{min} &: Minimum\\ detectable\\ dropsize\\ diameter,\\ here\\ taken\\ as\\ D_{min}=0.8\\ mm\n\\end{align}"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Eq. 1 can be derived from the emprical relations found in 1948 by [Marshall & Palmer](https://doi.org/10.1175/1520-0469%281948%29005%3C0165:TDORWS%3E2.0.CO;2), shown in Figure 1 below, taken from [Uijlenhoet & Stricker, 1999](https://doi.org/10.1016/S0022-1694%2899%2900032-3)."
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "![Clipboard01.png](attachment:Clipboard01.png) \n <p style=\"text-align: center;\">*Figure 1: Probability density of drops per cubic meter of air of certain size for four rainfall intensities (1.0mm/h, 2.8mm/h, 6.3mm/h, 23.0mm/h).*</p>",
      "attachments": {
        "Clipboard01.png": {
          "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAE0CAIAAACwwpMFAAAgAElEQVR4nO2d77lrztfGZ5/r14CUQAmUQAmjBCmBEiiBEiiBEqKEKIESPC/uJ/OdjdjCGJKsz4tz7ZM/SMIys+Ze9/rp+57tSVEUpmnatr3rXgiC+FT+7br1uq5932+aZte9EATxwewbpK7Xq+u6u+6CIIjP5n/7bTqKoiAImqbpum6/vRAE8dn87JST6rrucrm4rlvXtWEYZVmaprnHjgiC+GzWT/eiKKqqSvy3aRrf9z3Py7KMMWYYRtu2eZ5zzsMwpAhFEMQ61gSpruuyLEuSRMzjuq6zLKvrOtd1r9drFEWMMcMwDMMIgoDSUgRBrGZNTgrxiDEmglQURaZplmXJGDNN0/f9MAwNw2CMkfiAIIgtrAlSbdsyxizLEo/UdS2GS5xzxljTNDuFp67rEP705OOR+D9P7l8cyeBLkL8WwzDGBzz54DPGLx5sHP8VB4D/Dh4cv3H8lvNg27b86QzDaJrGNM2u6/Y7mQXytyp/VwTbsronf5VN0wRBMHh28l2O40ye/ezxO4lTRJzK8qkjX4GDC2n8G0/+6jj5mqZZMgkVu+i6Dufrgi/m6UbwucTjRVHYtr0iWzcZQeTgJYRpfwamJV/FfMgb/CiTz85vsOu6oigG54/8rPwrjGPlzI5mHh+/oK5rRKKZc29mO6+eipORSD5Pxjud/6RLEPcSeVM4t+u6ZtIdZcnWcOqqCqnj2wDnHGfmpiA1898Z7vf76p0qATl+TE4PpK7rPM+PvWdmWdY0TRzHBx5D13VVVaVpeuAxMMY8zzv8lIiiiHN+SIZEvsk5jlOW5eDON/n3eCNymNsy/6iqKkmSrUFKHlnYto1IzBjDbfzMy3knmb4N5k2HcJJiAJrdCI46M+WfQAz5Bw9O/j25EZzbM0HgzzMfty78vV6C0DSN2A3nvCgKfL9Ioj87vqZpLpdLURSr97sdwzDOkM63bfvwi9N13cO/CsMwzrD+e/hv8WHMf5/zz2ZZBoUAWC/mvFwucRyLVILneVVVYd+32+1ZkDrDoJogiBkOH+BD5IR8CNsy3RuklsqyrOu667r5AcIZ5lkEQcxw+KDSMAzTNMV0b32QGn+SJROHrusOzA4SBHF+xCo//ruvC8Lk7pMkOUm+liCIE4KCFjGO0R2ksHDg+/71etW8a4Ig3oIwDOM4Fqkh3UFKVM9kWWZZlhAuEARByBwWpJqmgc8U/nYcR15rJAiCqOu6rmuhENAdpBhjYRimaSq0C0mSWJZFWSqCIEBd17IN1AHTPSj30jQV2Xvbth3HSZJE88EQBHFCgiAIw1D894DVPfG3UFqhvrSuaxpSEQQx4IDpnlyjj/GUYRhJkpimmaap4zjHFs0QBHEqDhhJeZ4nFvXg22mapmEYWZZlWXa/35Mk8X2ftOkE8Z0Id19wgE7qfr/LcvM8z7uui+MYRki+76P0z7IsGlIRxBciL6yxY3NSgrIsoyhK0xT2Do7jhGF4v9+v1+v1eqUhFUF8G7JDwQE5qTGmaYZhmGVZHMdVVbmuCwNP+BQ7jiOvRxIE8fHIC2inCFKMsTAMYSaNZb44jh3HaZomTdM8zwdzVIIgvodfQQo2LlEUobYuyzLlUy2ozCerYfI8h2GoaZpZluV5jjGUbdv3+71pGiqjIYhvoCiKX5Onvu/7vm/bFl1e4Frpuq7rusajJ9XtdusVYZpm27bPnr3f74ZhtG0bBIFt2+h7nOc5noXvchiGqg6GIIgTEsex67pxHOO/DA+JmGXbthxE2rZN09Q0TfGGjZimOf+CNE1d1+37PggCzvn9fod+ShwPvKju97uS4yEI4oTEcfxfkMJ6P+JUnuciRgxQFRT+DFJ939u2naZp27au64Zh2LbtIEqiyYqquEkQxNmQg9S/uq4552EYuq7LOQ+CYDIPpbP7S1mWMMbL87woiqIobrdbURTCggojrKqqPM8jgQJBfDb/XNeF1Lssy6ZpPM873NjXMIw8zz3PMwzjdrtFUVRV1e12w+GJ15RlyTknzSdBfDb/TNMUEvCu64Ig2LVVJBYQ/xz+2LYdhiHiFFSdRVEgj+44jnhZEAS32y1JEhpSEcSn8o9Jfeg8z+Oc790rYmE0gVdDkiSIUxhPIWVmWZbYiGmat9uNhlQE8TH4vv/LuElOVi3Jam/kpV0gZQ4BBNQJ+BveCYNcPtYBOedqD5ggCP2kafpf4vy4cPk3IjmFmr6yLCEEDYIAUk9Z24lwZprm5XKhMhqCeGvk+dYBLggvvR7JKd/32UNW6nle0zSu6+JveYoHXcLtdkNlsuJDJwjiCH4FKRSmeBLKL/UVCS+RnGKMQUKF8mMMnaIoGqSi8DhjjHw+CeJ9EbHiVwfj6/Va1/XhEoQxeZ5blsU5R9ap6zrHcTC5u91uCEayKTI8P4uigOuL/BRBEG/BdJCq6/p2u+0apNYJBaCKchwH6XMoThGnkEHHHDBNU/ldnHPXda/Xq+M4eZ7r1KMSBLGR6b57tm3vnXJerW+wbTsIAjH9DMMwCAJoOxHC6rpG6mqwuzzPwzCkbjQE8a7Iy37oLWyaJowQ8K/alUW2zVZBro3uH0XI8n+fHTAqAV3XnfFgIAjiDMRxbNu2sD/5NZK6Xq9I+sDUCdV8amMiItTqGWWe51mWCeUB5ndiAAXDBqTVB2/EaAtCUNJ8EsSZQfz5b9VLDmCMMYXWUZNsH5ohDzXYZhAE4r8IVc9GTG3b7jFCJAhCIb9cEOQAFgTB3omb7RV2SE6JSmPGGBJSwl8Y5YfP9AeQrXPOSfNJEG/BryBV13VRFD8/P5cHcjWvEpQUBqIEWnY9h5eLiLCc87EkXQaVydSNhiDOz/8YY3VdG4bRdR2mgqDrOvQWVrs/VREBMQi5cDwCAToECowx13XLsvQ8L89z8RoZ2D9EUQSBwgnVYQRBMMYYVvQYY4hNmEztN9VUmAwShujiERQkl2Upv+ZP72OIQsk6nSDOQxzHwjT8n+u6aZqiGef9fr/dbgNJpFoUzq0QWeTkFPJNvu+LZBMk6bKr5xgsOKIbDZXREMTZ+McYC4JAW+GI2vljGIamacrJfgjQfd8X4QaRq+s6OZyNjwr+7o7jUIM/gjgDQ8U5YofyNPnMjlWRpikM0cUjooZGfhC9G+YNPF3XRYO/wXsJgtCPGND8Wt3TcGUqz8SLkCRHH8hVBw8iOz4fpzCkQv9k0nwSxFFgKQ9/6/aTQi91tdu0bTuO40HhHnJtstEwYwxNB/80RMeQKssysk4niEOAuuD//yNn1AdK7j0wDAMhQPmWB2V9AOYHA/U53IeXFPHFcSw8iwmC0AP0QKJ271eQQkvO+2+U717tBgWyIboMWgoOQhI6ti+JPhAxBEFAlckEoY04jqcLjA3DuFwulsQ5FefPtiwM0eXH4zg2TXOwtCeknn/OPYXP54x+nSCI/Ri6IIRhWErg+lTIrikeJKfGUoM0TTnng4Br2zbi1J8VfPD5RAQkgQJB6EYeYjHG9sgWyew33RNwzicl5kEQjMX0mM0JbeufBEEwOakkCEIh8nTvl30wBNxyB2PkuVUFRD3iI2gyYdo3fjxJElm5iogDscISRWuaplVV+b6vUwFLEN/MryDVNE3TNPJaPqpwVe1sj4rlyb2kaer7Ppyn5KcQj0QRsnj9M5f0SdBNC9bpSMCr/wwE8d087btXFMVg4qMwQo33vR+wwRtbnjPGyrJMkmRsm1WWZdd1k28ZI6zTyeeTIPZAvvcPV/c+xrFE7tYnA+VBlmXj4ILSmeULmpzz+/2OToWk+SQIhcgX1K/pHtTYyA03TYPZmTwz2ojmJXxUxsAsWH4cccqyrHHGDZWAlmWNp4qTYFNJkjiOE8exckt4gvhOft315cmdbdvmb9R6gaMFg8INLtnjM3E57Kgm1+kwpHppoVNYp5PmkyC2EwTBtOJcLejOgjmReERzkOr7HvV6k0+J3qLPnpL985aAghvx5RIEsY4wDKcV53VdOyNWj9Z834faW2SjD1kIQ3JqUoQJPeekMQuc8K7X60t5cVinJ0ni+z5lqQhiNU9X96DVtiVWJ6QwoGC/q5mPum7zPC+KYjIjNhOn4EQcRVGWZcv3BeGVbduWZVEZDUEoQB5iMdWK8zRN5SZ3+nNSAhiZP0sYIQk1+SySTSsc0KFlJ+t0gljB0+mebdtRFMnjndVjn6ZpoigKggBN8baEUSUg0DxzEEYlzbPWx7fbra7rV2dw0MF2XUfW6QTxKk+ne67rFkVxuVx+HszkpAZTpCzLkMOCNMkwjCzLrter7/tCe3VsmgbSimcVwihzefZ5oSyfcUl/BiqTxddCEMTLyEMs9kj9wv/gdrtNzv5Elzqxfo9qkjiOkaWGbL1t2ziO5QWyA6d7YNzzakAcxzNHiHWAFSKDtm055xhbvfpegvhC5Oneyzmp+/0OORVjTFztsIUTWzdNc/K9hwepfqpb34AwDPeIU33fl2X5ZxNAgiD63zqpX9M9zrnv+1gIq6qqqqpxOgmrV1i5EzRNI8TWUKtPjtq6roPHeV3XR6VpECZmJm4D2cSAMAzjOF6XZkJlMnQelKUiiKXI0Wv87LMxBYIUpntt2zJp6oeWyM8Uktimci37q3DO59fdOOec82fP4jO+KvWU3w6rhnVvJ4hPxXVd1LrIUuqVinM5EiFIiS0ifj171+HTPcF8cqrve875TMd5TBtXi8vbtnVdl8poCGISCAPw9z/G2Lxeseu68QvkdTpoNcUjdV0jYzUGvbTOoEhgjOV5Pq8qyPO8rutnq3LIgkdRtM5QGDcKzjmZvRDEGLnvHusf+WwYHuR5Llb3UPWGxwc3fAyXxJwuCAIUwaHk7dlMCtOc8wwfoDWdP5g/lZy2bc9MDP+EKpMJYsDTlla4scvBbKaVE7JL4ilcaXjXzPV2htW9AZPG5zJQLcxP6zjnGz8XsvWrk1wE8WE8lSAIVt/V5994qpyU4M8YtCROQXixZTSEMpotgzKC+BieShDkCeHqmeT8s/BQP5VDANwOZjQByB9dr9eZ5lcQsm7RFiDJhcrkP7tsEcQHk2XZr8y1zuiIkdQJky/Q0M8fWNu2fzpMoav7Rlk5VSYTX06e567r/jGS2gmMs07YXgVJpfl1OlQa+74/szrJOUfjrC0rmJDLojL5JCuhBKET1JCJ/2oNUqea5Q1I07Su63nrKOS255uziwbuW+ILpJ5ozEU9k4lvZzDKwpgCQgE5d6UErAluTDDvByZ0f07WIKSYf9lGqad8SPhFzvmNEcQewLBkqDiH65uITWEYcs4hklJYwHHO1T2ZhUeI3NN84EBqSUk5MX4dKqMhvoT7/c45/yVBCIIALrqTb4C4SckVckKd1BgE6D9ftiROrXb1nNwUldEQ38OwW4x83u96DbxFkOr73rbtJTO1mT40MvO1yi+Rpumfii2C+ACGOikst0VRdLlchC2n7/t72Il0XfdSX4NDQP+FPz8+nKf+bF+MyZqSLseiGw31TCY+GHg6if/+/+qe53lJkmAFPc/zIAhge6T8SoCf1MkvMCSAlgSCefMpQZqmrusq+T6hhEBlMmk+iY9kKJzuH4tug16+yKfM17W9yvkT5zILZ3N93wdBsOSVSqSeAvxAnHPKUhGfx7BbTFVVhmGICmFgGAbnXO29+oQyzhng175kcorGxdfrdf5lSqSeAgypTNN0HIeGVMSHMTHdm3SA+mXpooL3ClKMsTzPlySn8Mqu6/6MU3AQ9n1fVVhBq4vr9UqaT+JT+UNxrjasIOSdPCElg2K9hbmkeZM8AVRqCuOUXEZD1unEZ/Ar8mD6xxgzR+BBhfNMJL+OdTdfQRzHCzUEyBMt0ZQplHoKyDqd+AwGpnf/Y4xxzuH523Ud/kUkAwqjY9d1UI0q3KYG4MFSFMXAFHCMGHmh08TMKzH8wQJiHMdKjtN13fv9jn4/0D0o2SxBaOZ2u/1KX+gMkPf7/e2GUWBhWZ/84sFq6bNXog5p8wH+Aol80nwS78tQzIlm6M/SGVVVOY6jRIE50Gi9ERgiLZRQIkItMULAZruuW9HAfYYgCO73e5IkZ/MXJIgV/GOMobe4ZVmWZV2v1+SB53mXy8X3/TAMgyDYvrO3vmAgSlq4iCay4wtXBqFc33yM/4FAads2CRSId2SYOAdt2+JqQWc+0zSRA1aoFYQ9rqqtHcJLpXPoq75wkohbhXJlJr7zZw01COKcDAuMn4G+VZzzJemVJYjVPVWqa/28lJzqH3FqYYBA/fAe0QR9It73aye+ClSbLQpSfd+3bYs7PLrvbTzLMQF591s6dN7LP8USUxcBZAR7fEVoiUgCBeL8YEq3NEgJoBvYOFn7gOkeWOg5JUCcWvhixClVo1eZsc9nGIbvfs8gPpJh7d4MVVX5vn+5XKIoStP07SROOxHHcdM0y1c8OedBEDiOs+TFKJ3xPE95whu6BCH7Yoyd35GC+E7k0/J/44cEURRlWYaOAKQMHADRpm3bg8LsZ4Rh2HWd7/vI9M2D/JHjOLBy3nywv4AxtGVZRVHAOee0XXwIgkGCkCTJZYosy9BpXdXp+0mVZRiVvGQ+t9B8Smz/drtVVfVn0fIK6rqGer6ua8/zrtfrHnshiNXIMecfe/g9yoiiij8LQVbv+AOwbTuO44VBB6ACZmFEeEnq2TTNcgcY13XhbogUFWMsTVOa9xHnYWjVAjMpKKRM00QeCjMOVWVl4x1/BhhmvjQMWWjqIr8emsz5b6+ua6SZXqLrOnSLpUakxGn5lTjvus5xnCiK4jjGQrvancGg6sME0HmeV1X10odK07RpmuUOUHEcB0Hw59RyxUAVQRZD6SiKFo4KXxo8EsSrDOvnsMjXti0KX3ZVWmJxPQzDD1MVohXoS2v5y01dBJNST9HqCnUC+Hu1fGGh5vNNq8SJdwHCzF86KVxjjDGoZu73e9u24g+F+8aVqXCD5wGTspfeAvH6S14Fk82T7/f7/X5P0xS1Lxt/tSWmVJ/6IxLnYVgWk6bps3GXYRgKd/wxYs5JOOev9gFdEaewpjEeK+V5rspC71lvdwRi+GThD+UmMwQBZDHn/5vePRP77GEf/KnkeQ7LgXm7OxmEG5jkLdRbcc5N0/Q8DwlveVNrDvrJUeV5XhSF4ziyAQa6nDLGHMdBGSZB7MRETkoP79LBeDUrklPiXS8lkuA+vHch3rPe7pSTIvbmhbKYHaPjJ4LAsbD8RX4XWsgsF7viLVmW7dokBkKtIAggTxePq1XPEcSYiZZWevgwMeckKDp5NXYgur2kX8fgq2katW55Yzjn8PkUh6fEAZEgFqI1SH0JaZoWRfGqHMx13TiOLct6aby5UOq5Ebm3+wrJKEG8yrAs5hlZll2v1yzLVF0AHz/dE8A7+NXPK1ocv/RGdNx6aRS2DmGdrmFfxJezdLpXFEVd113XKenCwL5jugfE9O3VN0LH8OobwzCEAcveJdwYUrmuS9bpxK7IseJ/M69T7h71VbdfznlVVagxeumNQRCgPumlZX7OuWEYeJfyeqYB8PxDgecS5xmCeJW5kVRd11VVIaXyaknanxiGAW+QL4lWcRyvSE4xxsIwtG17eRFyURSe50VRhK4/GsY4WF40TZMqkwnlZFn266SStQli1Ub0LlYra3rVHfwDgKZ83UeGePLPlwVBYNs2PHYgDtDZGRS/6atSe4KYAULlaY/zvY36y7L8bDHnJCvK+gSu685f/4gR8iPIvmvuYAwd/IfVjRMHMuxgLOi6bnlJxwq+J3Eugwt4neoSPjBJkjx7QVVVA2ml67pN08B6ZVeppwwqB1V1uiaIpxKEMAx936+qqn6gdrXoS1JRoGka3/cdx/E8LwiCl5JTYkKO1bQsy55d/KZpDn6jruvQ23U/9+FJoPlExd9X/dDEHjyt3RuPdJY3YlrCx9fuCVCOl+d527boDxrH8fLk1OBrb9v2Wedk5LxE3d+4d6nrui9139oOJptlWercKfFhyNO9XxKE2+02uAd+5wRtO9frVbgUuK5blqXneVBOrfAPQBiyLGtsloBwADUAY6yua3jjiRfgWc/ztLUjC4LAdV3P89DfUc9OiQ/j6XQPJ3eSJEmSZFkm1vhU8T2zgKZp5Owevli4yM8nmLIsS5IEAtosy8QMEXHK87zxer9t2/f7HZ4qsIIavECUzkwe5x76TyTRDcMggQKxjqc6qSzLoFqGoMmyLLVn8FeNy8YRGRNAOfQMgDANDTshWJO/f8z4PM+b/FEwrXt2MCidkQsD0ccMzayU/9AAHRt930fYpWhFLOdXrJDngYwxOfERBIHanNT3SBDSNJVtmDDMwd9LPKdmvvZJB+GFwNXzfr/jhxDHgKzZig0uAYM7wzA0p8aIt2ZoHyxgjMlnP85phTv+nsR5/2hqgBYJnHM5KqVpOv89zIeMLZrY2+2Gtb9BYhty0BUbXAgsqlUZHBMfj2x69ytIyZcTVpQ2mlijR4D471cFKfCsLQL6pz97159LY5OdY5YfEmNssAtk91ds7U+CIEAqXdQw7BoNic/g6UhKtI0BYyv+l4BXv23b4mr8nunen+AesOVyjeN49ZeJmmdRXYCD2alc6X6/3263OI7RawtSDBIoEPOEYSiujgmPc7RZ31jigPO+f7hxiy1TkBKsM0SXwZW/4o34dTD/QkJq7+K7sixFTgpNgwZTYIKQQcNa/P2PMVbXNXTJ1+sVazFFUaCf7cx6+Z/JeaiurterWHWCGHrdBj8PKDy3mP+GYSiavA+bvs6CXyeO4yzLIOB61U/mVeTfHUNIGMvQkh8xiWEY/60491IeF0IemS05qTzPcR2KR2gkNQauweK/cDJAS7uFAw345K3uu2fbtp72eeOPo2cQR7wj8kjq13RvfBrNXCdlWQ6eLctS5BpEuQaaIYsXUDekMSI5hYsWc20MbZbEKfSFtm179dqZ/tIZmclGpMSXM5GTQnaTMZam6e0BFtEnN4EXi7wV0iuCtm0xaALiAsAwYedP935g7iOcp8TlCrHVn2/P8xwL/JzzPM9RMPjqMWD49vKhKwJiF53eMsTJmRhJPWufO44pSITjWXExYBUPf89MEuXsKSGDGfdgLiwWH2ZA/0743mEmvzohDaHAUSMajAcpm04AjKRwMgxHUnmei5HU5AJf27Z5nodhKI+kmCRVnx9/0XTvGZB9yo9gFWzh2/M8R4XAFlkDbIsPDBNxHG/8CMRn8DQn1T+SUOLfZ+crqtvbB4wxsUU8NfkuzGvgiqsnWftGYLonT3lc111ulIoVOnzDW+QjWO440GOT/Ii/GeQuMLH4JUEQ/Pz8XC4X+d9nHcOxojwpKRBPjd+F16P/EnXBHSBMVxzHQdEvJs4L3+66LibdM0XIS+CcQxgxuYWqqhzHwYmxU7sHmDp0XbdT2TNxZlBXj056/z0qh7G2be8PUMj6LJcpRlJ4F5PKLOB5NvkuSpz/CeY745XTl8Aq4ZYtYEQ2mHZhkIUHMd7ZVTguzAL32wVxWuamezIoEJt8Cg5qYlIgj8+xovzsXYZhLFcAfScvzfKeMfPbycxM64S5qHgEYxz5Bft5JwAsCxyYzicOAcnrv4MUHDZmwo0cpDCwwjyOSfmp8b5t26a6rXnGFsDrWFI047run3FKRMxxSNIjzUXApdPme8Di23SQGk8Rn0337vf7YECEXJfcLWsMiTkXosrgKQiC+fk1gtTMOAXDJcSpwRRyiUJCFTgMWmz5HuTp3k8vxSaRr0Ly2zRNtaV2dV1HUaTNbPutiaKo6zqoNLfg+z5GQ/KDaAUEk2J03GKMPUvSd13neR7ckJumQc6x6zrf99EWcOMRLieKoizL4KWhbafEIVyvV6H++7W6Z9u2YRhFUaDGWPmO4SH7PW3WtxDHcdM029vY5XkuCsgFuPfId6CZuxEy6AhqqAp2HMeyLJjwbDy8l4Bnw5a6d+It8DzvV/yRh1g451DJhburWlkdiTlfYrvnlMC27clk/HxOasCxpTMykK1SNv2DmV7dG2gy+9l1unWQBOFVVCWnEO/Gy/lYOFt+tR8uSRdgpXj7MihxToZ+UkDkocQjymf+Xdc1TbPHRPJTwajW9/2N28GUbTyLx7Ls8swjus6cYcKOMWBRFGc4GEI5E2JO3E6xPIflHqj11N6psE2S570KdOTbtwNJwfb5IzQBB5bOyJBA4cNAIGLSrI71fT9Te6E2B0EuCOvYbjQswH1ie5w6vMRPBpXYJFD4JGQ/qf8xxqD6kwdaYvwvz/62g6VrhRv8EkRb0BUt2gcgg+553qDpxquglZ7jOGcQBCDsRlFkWRayeMceD7GdYQdjRCWcdmjKNl6iVrVj8jhfB4qHB0qC1ZtK09RxnI03DBQueJ63U6UxYyyKouUm6HEc53nuOE4URTsdD6GNYQfjIAhQEjx+qdrpHspWT7KM/Y6oSk71jwqB7VPIcYmfQuA1+tJb2raFCPAkU1FiBQOrFtY/suZt28KQSHT0fGlxegmYGtDZsxphMaxka1s698ngqPZQA6wIUgCTPhIovClYDHmqk9p1lYR0UttBxlrV1sIwVDKwxYqMqnVbfEaklpB/WCHdJAeFt+apVQuKj3MJtYpzKjBWwnyL9ldRpSMXbuvbNyWQe22vA/GOWjy8HdNiTsZYURRo5ylQWySFxLzCDX4naZpWVaUqXY1OLdtT8nAWReGxkgNjT/xdX4JzDhWr7/u0svxGyIHiV5CCjFMGimRVoMB4e9Eskabp9XpVddXled40jZIbEkLeQAVe17Xv+yvWAdGcYuMhQb9q2zY1TH4Xmqb55RwtD7EGyz1lWaodvUNJuK4xHDFA9J1WgtqkkpySl3MISBIp2cUKYEpFLR7OD+RQw5yUsG5BZQz+xX8V7ptcENTCOVe7gKWwvkSUzgzS/Idbs8JBgZaYT46ck/ofxlOIUHVdw1IKDyIPqnAUR2JOteR5Dl8n09MIIB8AACAASURBVDRROYxb0OoN3m43y7Iw5Nl4bDhzLMsaeE4FQVBV1fbtrwYZPcdxhNs1cULk0/j/gxROKXizoHilrmvlzpwUoZQDjbVpmpxzxphlWXEcr761YIzteV6e59uLXYIggPknBi948Aw3KlTRX69XfNLDj4cY87SlVRzHkJ6LG51yCQLppNQCwbeYRCtReyqsZ4b4jj16dsDW6jxTLUQoJElJ+XkqnkoQoihK07Rpmqqq7vc751xtGRQVGCsnyzLcWrA2h96rGx27kOe2LKvruo2/F3QJhmFYluX7PmZ/excALz9mzvn9fk+ShBb+zoa8uvdv8KjrukVRQObruq7aFrJd1yEdQKFKFfgm0zTNskws8G+fv6Drz7P+1a9uCvZkVVXp6VxtWdbyF2OGGwRBURT7VUoTL+H7/q97hjzEMgwDi8doVwVbSIVDOCqLUQ4UQP1D3qF2PoXM1Hjet27NFyoHDbOql8qGsOwgBnf7FUsTy0HwmS6LEWaymOsxykm9A1juQMkRY0yt8elkh9HVxYOIUzsplcqyzPMczbvwx/KzFyd8/1B1nSdr9rU8zUnhd5JN+9X6mdEyyh7keR7HcVEUpmnKChIloPMC5n3dFC9tDSmqqqqUGGMNQC4VCYq6rpumWZ5mEmuOcRzDbIu6Zh2LfBr/9H1fFEXTNEEQjNPkmPqp2nFVVVmWqS21IWS6rnMcR4mAQEbEFFz20NMxxjjn66RGnudhYqXwIAWXy0WsKi4E/UHkz+J5HntU+Sg+PmIBw+aguAUZhlGPUJs4p9W9vcGVr7yYNk1TjDVut5twSb/dbqvFkPB7Ok+jF5TLyI+gJsyyrJM0N/rqxcdnc8I9yusoJ6UHNJ5SvllR3KfK0GqnLn4Ky27u9zuWjw6vNlVboHZ+hjkp3/d/Rlwul5+fHyWL0AIaOeshDMM9+huWZVkUxWBatAWI45VLUhTW3KCY0XVdx3GOFSh89bXT9/3tdkvTFFX1UAOmaYoTUW0tKHR90Dco3CwxBtJz5atUWFRRu0hfliV7SNJPCya5+h0UoCxDSWYQBEEQfIMyHklVcUoM7YPl8KHKW1YAF4STn44fA/I+yu8HiFNqtSkoxDl/g89Bi4cwDPcWVaHnAH7KsizLsvyGywetYWBv188HKfwkCvdNOSnN7JScQkxRe7W8S5xCxyPk5rSNayBu1LCj8yB3i/lPJ2UYhm3blmVlWVYUhe/7WZap9bKg1T3NBEGwR3IKd3W1uSSMzq7X60lW056BkVRVVZ7nNU2zTi+2gm+7cH5ZjcvRCz3L8LTruspva1QWox+1XbBk9phOYi558rQLUukiPQ+P42+YhekE5vT4+6kEYQ9ouncIartgycB+c484pba4Zyegc6aFoD2Qa/f+PRtu7TSE+7ZR6xlA06o9KlHQXEthexj2sCUoimKPA1ZL13VpmhqGQU4vypEDhdYgdQZXxu8kTdO6rvdQ+mB1BUUkqkCcaprmPJL0SVDaHcdxnuee56k1XyMEukdSOndHyJRl6fu+2jonANc9teOpoijquu667nK5nDZUIUixRza9aRrLsvb4hr8T4Z8zDFJoo3i9XmUTNVXQ73cghmGkaSpCCbrAep7n+/72HxqlwqomaFVVRVF0v99vt1scx1hKU7Ll/UDV5BIHhcvlou2o3pr/7kxysgoGYPg3jmO2g+L82+QeZwOqZYiS0jRt2xaLdNtFiQq9ojjn8omHU/Fdls/atkUS8Fk2HW0EiHnkos4JMafoOonvWuGOaXXvDMCIUo5K4+54q7Fte3u8G1/hKNh6lzjV9z0S6oMvGZJx3PvLslSr2v8w5CD1bzy4Emtwtm0rzwUURfHz86N2m8RLlGXZNI1chYveZUom42VZRlG0cf5o23aWZfIjXdfBc/1dFtFQxJ8kiZgCVw8YY9BLn1y2ehRIRDxts84YQ3WxcMtWW1GZ5zmNpM6AbduDgieFQirIRwfDBIwdlm8Bks62beVZJKp8z186IwM/YvnboOnen2AcOjHd6x8uBSJ+zcyr10FB6iTACFwIJieNzLcwLu4Lw/AlfSYyOwimsgB9DxuGvcG3Ie73FKSWMJ2TAug+lOf5HhNmuCAo3yyxAmSjUdKh/G7UP2544ixCT3MlW0YGbTLkoZ8CRl5nU4GjjfP9fj950c9JeJo4H9y14EmmcMeUOD8VKOjfL32LOIUgiD8wMtqe/4YrwODkFMuUEC6cMNGOI3yLip/DQZCatmphjInfPggC5X5SFKROhVjJ3QlRhBzHsXLHuEFfyEGVb57nJxy2i2ns2QZ6Z0O+mQ3FnOiFC42v8lK7c+qGvxnMQfZr3+S6bhiGcDVRvnEs54t2WxisiWc55yc836BLiOPYsizqmjWP+PmGQYpzfr/f67rGiaW2kIUKjE9Inud7VBcIIB/Nskxtly2ApRgRp+Sn9Ng8rQP+tFDSn/Ygz8NE7R6GzV3XFUWhNkhRgfE5Kcvyer3ud7XALGEgfVJFGIacczTyk3cRRZHClpHKgQ0p59yyrGNbPJyWadO7wZoINFMK55mUkzotylUIYzjnatdhZESfUbSzNU3zXc40LFbu9828KUP7YDQHxQ9cVVXxAGsxCkMjjWxPSxiGhmHsmiXJ87zrup38TGDkCPdxnN/v0igbExfDMCzLehc9vW76vodYZrIz9ca7q1hEBKSTOjN7tIEZIzqM7oEQVbzj2pnc4oEYjqRQ6wgd+oAttyPU4Mj2HZSTOjM7tWgfIDqM7rFxnNlKzGf0g0X3uq4dxyFTI5n/T5yL2IHUI9YdkiTZElOyLBsMXylCnRxotXf17UXCOIqineIUivuiKNopT78ruE+EYeg4zjse/078Wt1zHAcrer7vm6aZZdmWDAKar8qBiXJS5wcrYrs64WJGdr1ed0rBYNKaZdmbCpEgA6qqSnkD+jdCjhu/rFqapsFonDGGiK522Ew6qbcgTdOiKHadcWA85XneTvMyxMGqqtSaGmtDDKnQB/PowzmAPxoxVFUFs6H5FNLg9IKuqqqqQRgajKRoxnd+cIXsrTO0bRu26/staUGSfv6uM8/AkAqder/w7j6tk4JnI2Msz3OkzJ+tNcDrQxTX3O939mg6KjdiK8tS3gLppN4IDcqpfp8OowOUG8zqB/XSX+XkKdfuTXQwRliBt8b47IHxGAKceNa2bRF9ZpRp5ILwXgxcMXZi4HC2B2EY7mFHoxNcdye0oNmJaReEhcBwKgxDMZKCwEr4JYZh+KywnnRS7wU8NjUYnqDH8q6XXxzHe+9CA8KU6ugD2Z0507uFYLaPLwtzPfHFYSb47F2GYUDITkOqtwA3cA3XNpTiu+4CofDdr3BIGpX73pwEWPoMBtcrm4MijSdnwZ/9PciEoZsIOpSt2zWhE9xONOSeMSPbtb8e5/y9ujlMIhqRfqTmE0ltZLen++6BQS3LJEirC2E6k6Z7SPJNvov67r0pruvqMb2FpcGuu0Ch3Ht1c5gE85JPNSOWc1K/RlIIz5cHM1o4jJUgKQDi7lQUhdwuaUBd128qXflm8jyPokjDfRvO67ueIZB6Xq9XIXlPksSyrMvl8l4r/TClKori8zSfSZI0TfOfYkmOXoZhuK57u93gEs0YexanMZIa5KFwG5xJtd5uNwxW3z1/+YVoS071fe+6LpaY99sdHP3jOOacc86xozzP915n3INxI9K3Bv2J5TDyX5BC/nvgJ/Usl3m/3+M4ll8Mj0Q4ZjzbPUkQ3hoNczEB1A9QM+60C8SpQWpiDy92DeCziGj7AUzrpJBaks8J5Y0Y8jyHzljhNgmdwFIuCIK9UyHwjdl7XIMJgRx5Mdjfb4+78kmaT3kd9r+cFOZ6lmVdr9coinzfz7LMMIwoipIkUTLpNU2zrus3LfskmqZpmgapnKqqLMvaKRWCEqsgCLquS5IETcn32JFt2wiFYlWxaRq5m8N7EYYhnGqwGouv8dWN1HV9bIar67qfn59fDRZE6MK96xlKwjOVxbw1GAWjkKXfoS2jAGkB5DcZY6hh2GlIhXklBBDjxstvCjSf65p6cc7PMNf5W8y50wlBEoS3BrGp73vM+HBX23WPuGB2VWCiz6i4aX/Mij40nyumrkEQnCEHLwepf4OBlu/7Pz8/l8tlp75g64agxBmQr+S6rvW4X2Kks58lAzRTWHFO01RVZuNYHMeBI1hd15ZlLbHEgTbI9/2iKJIkwd9HfRXQH0yLOTFFhyIc1XlqYyq1mX5r5EUPrAXvPXkvyxK3U3Ra1rB0BRXC8oHbOVfTIINEqSyW85eMEOENzzlP07Qsy6PGU2JsOzHdw+qe/KUHQaB2pYMkCG8NZhBY2uOcm6apcyFs747wAtHN4c/X2A/OGapEPT8ECgt9IM4+3QNilAUpuZ4BHnF+IF1B3zMI4mzb1mYpB0GMhpa/6Obged6zvETXdZ7n3R7AkvyEk0TRdx4xF2v3f877znjVy9ELC3l5npdlucd0D4JPhRskDgcnjLbdBUGg5xRCz87JjxaG4SBl8S7qP5QNvMU1OGd6J8ru9qhdhMHLB9j6EAJ4Tun8QTnnekTh+Gjjq2A8IYrj+I0yrWEYamiwuBoRhaZzUvupUQAkCBShPgxt2SLBrh1GZZDNGcSpsUDsXUZSAjgonDawDkylh4nzXb9rSpx/Kq7rai5502Nt3D9u7IOoJBapcWt/R/UfFvJOm/WfTpyjLAaa+uv1mmWZbGehhM/z6CIAOqHp7BtclmWSJBo0dxBSIVkuHkS/LMdx8OCWRt9HAeOEIAjeoGuWCF2TZTFqb49Yula4QeI8oKZE520ZOSNt8yyYgpxz3LGFlwQK2lDgcb4OKOV2KvgiDidN070LZQZottmEEfupLmZV4Lc7g0KqHzUQ+scYi6LIsizHccb/qm233XUdikUVbpM4D1D/6mzGiVWqXTuMyoRhiPnRCVVRGwmCADPow+1J0UpdfuQfYwzSYei+YFUB6rq2bVvh7qEuU7tN4lSkaVpVlc5OB6ZpohOynownBPeWZa34jHCe8TzvnBWC+CZN01z36VSBZUeEo///luRRFvttejde1Ng+iqOc1McDGaTmOZFm51/s7iWpEWamYRiiMPbMnjDQfB6blpnOSd1uNzayD1buzPmOi7XEq+hp0T4AWRVtcerVdJjrunLGB00A9zk0BUCgcGCbwqeJc8Mw0JusLMv5RgzroJHU96CtC5aM5sw9lhcXJpvHB3bmIAWg+Twkm/40SGGgLs8M1e4YHuen1eMTCsEFrP+3Vj78nwfCnSXheDwbPX+Q6o8TKPwhQdiv5RRWFt/ityG2Az8T/Qv2Orva9I849WcGB3am4r9pmr5R6gMtHvQMqcRQSZw5P33fi7w6uiTI6w4od1CVty+KoiiKd5TnEutIkqSqKhSW68TzPAxw9OwOlra4jOdf03Wd67p1XTdNgyCu5wi30zSN53m2beu5fh3HwWSTQYIgP1HXtfEbhTs+4bIrsSu4w+mvuijLsmkabX2J5NKZZyc5XoPJRBiGyLvrOTwlIIl+jEBBHmWx36t7yqHVvS8EsyGdySmRT9BmliD41NIZGQwA964nf1pgvHdA/K+5O/E1YBJ0iI4ZZc86G39Ae6HBPvRAYEfXNI3jOHoEtP8Tf8Gz2bIsuWzFdV1hg7cdilDfCTqAX6/XXdMZor8ImpiyxwzLcRzGGOd8v13LxHGcJInjOO+VcnoJ3HiKonAcB+5ae+zlvzGNGF9hzRiPin9JcU6oYm8HKCzq4Z4KjydMSZD90SyGQNL3tJpyVewnUDjSBYFyUl8L7oIartuxxuWQ1sSIU9+gCsRqgFqBwjBICV2J/RvlPfLwYShOfS0DW9idmBTi7d0JeRLM+N7LWXgdaCPEOd/+47ZtC1nWr8Q5QhIUHCJCIRslBOhKQNpLv2qGOAmu6yI5teteUNE13nWappp7T+Fsh8+ttp0eAhZwTdOEjGn71n79fyak3W43tUM4ykkRyGIcZa6m2SwBYJ572q4HalHloPC0OShsm60HjuOoda02DOODl2aJSbquw0l1uVywNo/BxSGG95zzOI41j6cQFoui0GkHeBQQKDDGFGo+/wtS0MuKRVOMuBTWxDDSSX0fiFBBEKAglHMuOhcc1fUX7UV939e8X4jg9e/3ENI0TdMU9n7bt/YrSDHG0HXKdV38oXYuTRHq2yiKwrZtoaOBv3BVVfBTVGtOvRxYlUM/pQ1k0BljcteZD8Z13fv9jsnZxrvRv8H/u66zbRvjNDgIb9n6eOMKt0acn6ZpBipK13Ux0YvjuK5rnXJwGawf6Z9/wevOsizN+z0ExOUwDC3L2vJD/xek4D5+uVwQpHzfT5JEodycUU7q+zBNc3B2VlUl1m7yPI+i6KhujGVZ1nWtfzSXpmkQBEfNdvXDOb/dbjB3X/mRB0l1aDqgulReQ0ire1+I3Chp7JwJ+dIRx9X3j6XGQ9bdBl2bvgH8+gtVY9OKcw3NHfM8pyD1bSAQoP8Q8hSDFwzc4DQDk4ZDJBGIUx9fOiMDs5clP/ecC8KuOYIvGd8SMigNud/v9/sd46bBC9AFS2eL9sHhQRKhPzsWBAFWOY/67PqB5pMxZlnWC9N8OXrhBIIuHMXGam9xVLtHTAK544Fzn0OKkAFKZ07SOlgbKGycmWjLI6n/yQFrbGehNnHedV1d157nUWUMIWMYBmpWoAPUD3IlOAC1pWB/gtYkSCp/T3NvTPx934efuPydQ1v3a5ylK3T2PeWkiFmOTU71uuqfJ0FqbG+7yxOCQqWBhw++jemymL2hnBQxQxzHByan2KMI2bIs/ScqJpt1XX9D6YwM5/x+v2dZJgsUBt0VtAYp0kkRYPI0wIV6iNGwAMobnWJ0+cqEdOtLSmcE+N055880n7pHUlQZQ1wul8vl8ixOxXG8a+HI5XL5+fkR/44VhmEYQmy53zEwxuARjmO4XC6ixg2LXx9WOtM0TfQgSZLJSBQEwe12i6Jo4i6lef5JOakvR6yZzKxncc73y86whwYaTkSMscl9wYB4p2NAJEK2HkWyjDE5L4Mix4+ReqIBIlzqkCOfMSAMw9A0TXm1V2uQStP0C1ODhAx872Cv+Ow1u3bBwk1b/BcS08lXcs53SuRD3SrHIM75QHmPa/Uz4hSClPgvFCczJwCC+GEe57ta8RMnp21bjKEwhJm5Avdr0S5GUmjVOTOmg1Z+j9sqY2xJLc7HSNLxc8uP4Juf+X2H9sHaEE2HiO8kyzLDMERPlxkjIESHnVLI6MXked71ep1p042EblVVe3RCXqLGCoIA8jHdHYNVM77qcQLMfC75+6HVPUIfRVF0XXe5XJCWnncrg+Gi8gAxqKOwbXt+PVF5h1HRHHDw4KT2gnNelqXneW9dOrNxuYxW9whNNE1T13UQBFg+C4JAdPF8Rp7nWZapvT4H8SgIgq7rZg4DyoAoilTFKYzdBh8KTiaTr4ckHeJsJQegH3zn8jePL3ypuH+nWWjf9/f7HZbS4pE8z4+VFBMHEoahYRjyI+x3DnsS5ckpeafIOg2OahLkelXl8tHPRpiWIE88n/xCvc6blviNc1JYOph5i7xosFeQwo+KimIRp8hP6psZ16sHQbAkQMRxrFANMBjLL19GVFuEDNkBnCYZY0s+IELqO97mkSYXg2h85PlvUlZg/PR9r2A8NwLpRqGgReFoURRN06ht7kC8NQszAJ7nQRCgZI/yvOOlcmJRHo/LbCNwT+66Tqwk/EnXdb7vm6aJy/5dqOtazi2iSnH+d3ccB04JjO023RMzOwxTRUB9x/sAcThqZ1tbOKQT8gDO+cfPSHQUGBuGgdxY0zQiZBqGodkHg/gMkL1eb5KtDlHcd+CRIMVzhm9jPxQUGKNJn/wIlic8z4P2wbZtUbAj5ncn+U4x6zz2GLqu20N98ypnOIaiKJbogKCc2skkoGma5WuIQRCgw6jyw1h+DHmeoyuX8mvqJBepzJogVVWVCEYAdYOYqOOLE13GwjAUXnozwjmdJElyuDqu6zq1PQ3XcYZjWG7PEoZh0zR7HPOrzbXQYVR5nHqprjiO4z26zkRRdAZNlvyh/jfzukmqqoL4TZ64JUmC4mHGmGVZ1+sVbntYapV3fMI4TRzO8ltXWZaWZYk61QOJ4ziKIrU2s6/ewlHfZ1mWqlw+O+VI6uUg5bpu27ZZlomGZRiViOUJzvmzSPyS4nz8SrEStOUPBMqmafAHHhEvkB8RCjT5EVX/4gDqup55zeSXhmfHD7LHKY4k4JLDwB9YYDJNU7wR17+8TfmrEA8ahlHXNeb1pmmu+BKYlLv88+4t3hIEged5YnlL3tSff0yCzw6tKT7Is10PPkIQBL7v+76fpumW80E+zcZn+zycc8MwFK45nqQsRP7sKyUISZIkSYJ60aIofN/H+gtjLIqiLMvw1ACs4K6Y8a3+7cV1C83en7t49cD+ZHwwgx0tCa9Ljnb58S/Z7Pz2F36EZ/8dJAQRI8R/T3KRIGSLu5oIpoP7B/sdwTeGKoH81T07wj9DHhvdV9gj5TJzBxJa8Mm4/GxHk5/iVeR7JGIINvjySEpsbsWzEGit2yNBEDsxc1NZMvr7c7ws72jhG2Uh1cogJe8bEz350A/PFxAEsZzxaOvPfxe+bDBeWTfUWilBkHeGASTWXLquy7Js3BqLIAhiHet1UvJ/y7JMkgRuza7rfk/7MIIg9mZ97d44e1pVlWEYavuJEgTx5awvixnPNufrJNFT7Hq96rTF6brO933LsiZrCJqmsSQ0KxufNfDZFcdxJkXbRVGIziWXy0WPnE/umKKt31zXdZ7nXS4Xy7LGgvuu6xzHEaeEtqNCHT5+HW2Lm0jOXK/XJEnGO+26LkkS6LSjKNKvf4Y4hjFd3WLggoDCSPa7McaumKYJF55J/xrUQEE9DAtEPUeF49H5PWCnuIVM1nhzzm3bTtM0TdM4jvX4/8OoRHRtWWL7rWSnaNMCpdVgp+JExVmh55CwU5ga4fA07LTvezRlgdLKGJl2weXKMAy8THMLFfmU0BSkXNcVjjlj87OdkHtOQPowCENhGO7XtmgG9ErRHKRw1bEn5mqy7ZcecGUKOwHYDO29U5wGYqcIzfILYA+y92EMQMMF/I2TVoPfA4Ty8tUxcNRD0cjehzGJECrhkDTZBzdNI0vSIZDbe6fIkYmlUNu2B9Mr6Iw9z3McR2epbVmWOBd1kqYpRnCT33xd11mWYeYlagl2xTAMXAbiADZKAReC0Qr+hlZefhZSRpwS2mZeruuK2hqcohpEPBg4y3rOwfePywe1JZoL8j3Pw1BG63SPSaOGwd1sPwb3SXR8k18gJoPwadA8oGV6R1LgmbUjhv33+x2TIJ0HdrvdUA+g06SpbVuMKwdjFpwJcRzHcYxVIG2HJApddXquwW6UMTYexmJUgTwJm23mqpYgCDCCw22s1zbdk897jPM1ZD0GDRfH16d8DNomoWLXpwpSMvOdOxXSti3UKq7r6oxQiETPvOvEWSFPiDTQtu3tdsOxaUuPtm0LK8rx2ZjnuTiMhUbP25G/cxGkNE335LpNbUNHzCvFf8dj+/+WD7TbyOjcl8xkgVXTNINlLD2Hhw4oaDWurUrB87wsy9I0hWfs4NnBKcG0uAKgFQ0yEhjBaVhKq6oKO0V3v3EyBHlk/K0tRRNFkWEYvu9jOR5ri5pGUmII1/c9WsJr2KmcDpRb5oZhiEE+ewyt0ddbsyUrO3ok1bYtsiEY2wZB0LYtbmUaDgxZOdyu0U9Yw2AKn7Qsy7Zt27a93+84JbCs2T8WHPFiLHvtfUj9Y8iPI8ERaphbIRqKcaJpmsLvG1+OfBpom2ekaRo+QKDM81xTkEIUAExLQgrIhlbiGzcMAytZiFxA9lTWA9M4qhfIq3iio0//WGDCV6EnJzLwGmPLOqYo3ylCkrhxIkYI9LiqIzGEkZSe76F//PqITRAZ4KpEdrLve8wBhWxIfzctcXLu1S1mTPewT3JdV+dkR1htiJ0Opjyv9Sn8aHbyqzkPkxOWyflv13WqbOQWUlUVivO17RdiTqRBkHVij0U9HINo3RwEgf5Kkuv1CmGKviB1IB9/7RHEB/MVQYogiPdF0+oeQRDEOihIfRRVVS1UeMhr7aqATfhOG99IXddVVeHfP1+5cV96Vuu/CM0Ze2JX2OJVmOWvXI6Q+BuGoVbEkOf5loW2cUMXzvnkYi4UQxuOtO8fBeQbN0IIaCT1Udzv9+WrMHvc7bFAcbvd1LqzQu64cSNC+JLneVEUkyWKURSJPjSrgfTxDN3rPgMKUu+E7/tZlsHnCKXRA1+k6/Va13Vd147jRFGEZ4VPVpZlsIsS1yfsk+CqitmQZVk/Pz+O42DWgzaLvu9jU+MLD/ZMPz8/nueJmabv+6JRleM4l8vF8zy8F72FrtcrTL7gozRwbsKz4nMlSYKeoIhT6E3rOA76P2Kb8Dx61lwTByOCMuc8TdOxfViWZWL13fd9fHYcZFEUnueJz1UURZZlSZL4vh9FUdd11+sVbxG70FOn/RUcO5AjXgLjFCjmoT8U1V4YJjDGIJ5mjNm2XZYlRjT9Q6aI94qSUTxommYYhsJK6Xa7YTiGwi5s6na7Qd03mCWhALUsSyER7h86VaiWgyBASzhIliGndF1XyEdh28Qe5kGQ0cVxjA2maQr7FBwDJm54i1CH48BQUzL5veFTyBLiSWE3Ornjb1wd0OkwxlA+IvaIbxXyfbxSOD5iWqqtQPUboCD1TuBSwd+oIxG+BbjekAyCkh4vE7ZN6HaLB1GSIhxjcF2hTkLeF8wA5IttkGxCGBLZIqHbZg8xfZ7nqD4RsRIbxOvRiFh+r7D9wh4hd+77HmEUjwhBtggEiBQz39uzIDVInMlhSwRNUUnTP0oU+t9lIqZpikMSW5j0LyPWsbKlFXEIYjLCHjMOVDaIFwy09TKyMPtdqQAAAjVJREFUolX0KJbFzWhHLP6L1sQYAclvHK8eigOQy8ixwpVlme/7YgsD3yJ5y/ILMGszpvpN1nXdNM3lchFPYVq6RKw7qDbHZ5RfMPjq5OzeuCBBPDLonDTuu0lshHJSbwbOe+RoUAs6WLp6Jq8fRJDBH4wx13XllJNwoZNfMwhkQLxLfiX80qqqQrUqBlCDA5MPVX4KQx4MweT6Svaw0hd+knEc/9mbmk31gMMxDz7L6pgiV1xNPk5sgYLUG5NlWZZlGHeIJuCD14jLJggCWLIgdS0/JXJDcMlAkpgxhqwQY8zzPPEueXwBWzgsvUVRJKKVIbWuRdYZefHxKGxg4IOkDzLQRVE4jiMW9fCCIAhQUIbPgsy6HBqQwx7sCC+oqgqbxRvHq3i2bcsiKWPUgPcZg29S3sKStxN/cPR8k3gB5JvxtzDth2mZqOZHpknodOA/gaEHbJuQ38ErUQqPwU7f97BYElnq/pHNQTJbPCgDi0s8GwQBUmbYeP9IMIteGJxzWUMUx7GwxwnDUNjmiFy1eFZumiBCp2maIpsmskKTboJIvYtzHs65468XXwv+Zr9T4BjcYdA6eCU+F/4WhfvilcR2qHaPmAPjDmGMfyDjPo/Pxime56VpOs4iPRvvyC+wLAurn/L2J/+ePADxIEaO2yVXBKMCY2IeaIXOEKQWgmy9MB55FYieNgaXrusulwtGqVu2QwDKSRFzmKY5Noo7M4ZhQB+w7u1b3iuo6zrPc4pQqvg/idzlwFLRZUYAAAAASUVORK5CYII="
        }
      }
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "We see that:\n\n\\begin{equation}N_V(D) = N_0 exp(-\\Lambda * D)\\ \\ \\label((Eq. 2) \\end{equation} \n"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "When we, as Marshall and Palmer did, assume that a drop falls with a terminal velocity:\n\n\\begin{equation} v(D) = \\alpha D^{\\beta}=3.25D^{0.762}\\ \\ \\label((Eq. 3) \\end{equation} "
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "then the probability density of drops falling per second on one square meter of ground surface is:\n\n\\begin{equation} N_A = \\alpha N_{0} D^{\\beta} exp(-\\Lambda D)\\ \\ \n\\label((Eq. 4) \\end{equation} "
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "If we want to know the total number of drops falling per second on one square meter, then we have to integrate over all possible raindrop diameters:\n\n\\begin{equation} \\rho_A = \\alpha N_0 \\int_{0}^{\\infty} D^{\\beta} exp(-\\Lambda) dD\\ \\ \n\\notag\\\\\n\\rho_A = \\frac{\\alpha N_0}{\\Lambda^ {(\\beta +1)}} \\Gamma({\\beta+1})\\ \\ \n\\label((Eq. 5) \\end{equation} "
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "Unfortunately, due to the limited sensitivity of the intervalometer, we can only measure raindrops from a certain diameter upwards ($D_{min}$) so the measured arrival rate becomes:\n\n\\begin{equation} \\rho_A^* = \\alpha N_0 \\int_{D_{min}}^{\\infty} D^{\\beta} exp(-\\Lambda) dD\\ \\ \n\\notag\\\\\n\\rho_A^* = \\frac{\\alpha N_0}{\\Lambda^ {(\\beta +1)}} \\Gamma({\\beta+1, D_{min} \\Lambda})\\ \\ \n\\label((Eq. 1) \\end{equation} \n\nwhich can readily be calculated in Python.  "
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "One can find many parameterizations in the literature, both for $\\Lambda$ and for the terminal velocity (Eq. 3). A very good overview is given by [Uijlenhoet & Stricker, 1999](https://doi.org/10.1016/S0022-1694%2899%2900032-3), which also consolidates many of the different possible combinations.\n\nWe use the following parameters, which can be changed if needed. $D_{min}$ may vary depending on the sensor used for the intervalometer. $D_{min}$ can also be seen as an effective parameter as small drops falling in the middle of the drum will generate a similar signal strength as a larger drop near the edge."
    },
    {
      "metadata": {
        "trusted": false
      },
      "cell_type": "code",
      "source": "# Constants\nalpha = 3.25   # (m/s)\nbeta = 0.762   # (-)\nN0 = 8000     # (-)\nDmin = 0.8   # (mm)",
      "execution_count": 2,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": false
      },
      "cell_type": "code",
      "source": "# Forward model (see Eq. 1), calculates arrival rates for rainfall rates between 0 and 200 mm/h\nd = 0.1 #Rainfall intensity increment (mm/h)\nR = d  #Rainfall intensity (mm/h)\n\nnr = 2000     #Number of points in forward model\nfw = np.zeros((nr,2))   #Values forward model\n\nfor i in range(1,nr):\n    R = i * d    #Calculate rainfall intensity\n    Lambda = 4.1 * R ** -0.21    #Calculate Lambda\n    pre1 = alpha*N0/(Lambda**(beta+1)) \n    pre2 = sps.gamma(beta+1)\n    rho = pre1 * pre2 * sps.gammaincc((beta+1),Dmin*Lambda) #Note: gammaincc is normalized wrt gamma(x)\n    fw[i,0]=rho   #Arrival rate (1/s*m^2)\n    fw[i,1]=R     #Rainfall rate (mm/h)\n\n",
      "execution_count": 3,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": false
      },
      "cell_type": "code",
      "source": "# Plot inverted results (rainfall rate as function of arrival rate) \nplt.plot(fw[0:nr,0],fw[0:nr,1])\nplt.axis([0,9000,0,200])\nplt.title(r'Rainfall rate vs Arrival rate')\nplt.xlabel(r'Arrival rate $\\rho_A^*\\ (m^{-2}s^{-1}$)')\nplt.ylabel('Rainfall rate (mm/h)')\n",
      "execution_count": 4,
      "outputs": [
        {
          "data": {
            "text/plain": "Text(0,0.5,'Rainfall rate (mm/h)')"
          },
          "execution_count": 4,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEdCAYAAADNU1r0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xd8VuX9//HXh703yCYgCAIygwJucWsdrQsX4sBZtbXW1X7VWn+t2m2r1g2ogOKsWhWtW0D2HrIJI+ywAmR8fn+cE72JIblJcufcSd7Px+N+5JzrrM89cn/uc13nXJe5OyIiIsVRJeoARESk/FISERGRYlMSERGRYlMSERGRYlMSERGRYlMSERGRYlMSkVJnZu3NbKeZVY1z/fPNbHW4Td8i1r3KzL6KmXcz61zSmCsiM3vKzH5bCvt50cx+XxoxScWjJCIFMrMVZpYZfrGvD79I6sWzrbuvcvd67p4T5+H+BNwSbjOj+FGXjJmdYGZpUR0/PzP7zMy2mlnN4mzv7je4+0OlHVdJhZ+tk6OOQ0qHkogU5ifuXg/oA/QF7knQcToA8xK0bwAsUG4+72aWAhwLOHBOEev+6Iwv3rPA0mZm1aI4rkSn3PxTSXTcfT3wIUEyAcDMzjKzGWa2PayKeiBmWUpYzVQtnP/MzB4ys6/NbIeZfWRmzcysppntBKoCs8xsabj+3Wa2NFx3vpmdX5y4w+M+bGZfA7uBTmY23MwWhPteZmbXh+vWBf4LtA7PvnaaWWszqxITz2Yze9XMmhzgeAvM7OyY+WpmtsnM+plZLTN7KdzHNjObYmaHFBL+lcAk4EVgWL7jvGhmT5rZ+2a2CzixkLLfFxVbOP9aeMaZYWZfmFmPOF/jq8L39a9mtgV4wMwONbP/hc91k5m9bGaNwvVHA+2B/4Sv8a/D8oFm9k342swysxPiOb5ET0lEimRmbYEzgCUxxbsIvugaAWcBN5rZeYXs5lJgONACqAH8yt33hmc6AL3d/dBweinBr/CGwIPAS2bWqpjhXwGMAOoDK4ENwNlAgzCev5pZP3ffFT7HtWG1Wj13XwvcCpwHHA+0BrYC/zrAscYAQ2PmTwM2uft0gkTQEGgHNAVuADILiftK4OXwcVoBCedS4OHweX1VSFk8sUGQQLsQvD/Tw+PG6yhgWbjtw4ABfyB4vQ4neM4PALj7FcAqwrNcd3/UzNoA7wG/B5oAvwJeN7PmBxGDRERJRArzlpntAFYTfPnen7fA3T9z9znunuvuswm+pI4vZF8vuPtid88EXiXmrCY/d3/N3deG+x4HfAccWczn8KK7z3P3bHfPcvf33H2pBz4HPiJIWAdyPXCfu6e5+16CL8MLDlBt8wpwjpnVCecvDcsAsgiSR2d3z3H3ae6+vaADmtkxBFV8r7r7NIKkemm+1d5296/D12hPIWXxxIa7P+/uO2KeY28za3jgl2U/a9398fA1znT3Je4+IfyRsBH4C4V/Ni4H3nf398PYJwBTgTPjPL5ESElECnOeu9cHTgC6Ac3yFpjZUWb2qZltNLMMgl/WzQreDQDrY6Z3AwdspDezK81sZli1sQ3oWcS+C7M6377PMLNJZrYl3PeZRey7A/BmTCwLgBzgR1VR7r4kXP6T8Mv6HH74oh5NUCU41szWmtmjZlb9AMccBnzk7pvC+VfIV6WV/3kVUlZkbGZW1cz+GFbZbQdWhJvF+5rnf41bmNlYM1sT7u+lIvbVAbgw7zUOX+djgOKefUoZUiOYFMndPzezFwmuosqrsnoF+CdwhrvvMbO/Ufwv+u+ZWQfgGWAIMNHdc8xsJkEVSXF83021BVc5vU5QVfS2u2eZ2Vsx+y6oS+vVwNXu/nWcx8urNqoCzA+/vHH3LIKquQctaDR/H1gEPBe7sZnVBi4CqppZXuKtCTQys97uPquQWIvqkrvA2AjOSs4FTiZIIA0Jqu3ifc3zH/cPYVkvd98cVnP+s5D1VwOj3f26OI8nSURnIhKvvwGnmFleNVR9YEuYQI7kx9UtxVWX4EtmI4CZDSc4EykNNQi+kDcC2WZ2BnBqzPJ0oGm+apyngIfD5IaZNTezcws5xthwnzcSU11kZiea2REWXDW1naB6q6BLoM8Ly7sTVPn1IWhX+JIg+ZVEgbERvJd7gc1AHeD/lfA49YGdwLawvePOfMvTgU4x8y8RnCGdFp4V1bLgcuu2JYxDyoCSiMQlrNseBeTdvHYT8LuwzeT/CNo5SuM484E/AxMJvmyOAOI9Cyhq3zsIGspfJfilfSnwTszyhQS/1peF1Sqtgb+H63wUPtdJBA3JBzrGujD2wcC4mEUtgfEECWQB8DnBl2d+wwjaj1a5+/q8B8Ev+csO0BYTl0JiG0Vw0cEaYD7BcyyJB4F+QAZBg/kb+Zb/AfhN+Br/yt1XE5wJ3UuQ4FcTJB59P5UDpkGpRESkuJTpRUSk2BKWRMysXXj1zgIzm2dmt4XlTcxsgpl9F/5tHJabmf3DzJaY2ey8m6BERCR5JfJMJBu4w90PBwYCN5tZd+Bu4BN37wJ8Es5DcKNXl/AxAngygbGJiEgpSFgScfd1eXfDhg2aC4A2BA1oI8PVRvLDJaPnAqPCm8AmEVzSqOvERUSSWJm0iYTXxfcFJgOHhFeJ5F0t0iJcrQ3737SUFpaJiEiSSvjNhhZ0H/46cLu7bzc74P1LBS340aVjZjaCoLqLunXr9u/WrVtphSoiUmGty9jDpp17adu4NssWzNnk7qXSN1lCk0jYrcPrwMvunneteLqZtXL3dWF11YawPI2go7Y8bYG1+ffp7k8DTwOkpqb61KlTExa/iEhF8M6stdw6ZgZ3D+rAg+f2xMxWlta+E3l1lhF06bDA3f8Ss+gdfugHaBjwdkz5leFVWgOBjLxqLxERKZ4F67Zz1/jZDEhpzH1ndS/1/SfyTORogm6454R9H0FwR+ofgVfN7BqCLqEvDJe9T9AZ3hKCDvqGJzA2EZEKL2N3FtePnkb9WtX412X9qFGt9M8bEpZE3P0rDtyB25AC1nfg5kTFIyJSmeTkOreNm8G6jEzGjhhEi/q1EnIc3bEuIlIB/e3jxXy2aCP3/6QH/Ts0TthxlERERCqYD+et5/H/LeHC/m257Kj2CT2WkoiISAWyZMNO7nh1Fr3aNuSh83pSyG0VpUJJRESkgti+J4vrR0+lRrUqPHl5f2pVr5rwYyqJiIhUADm5zu1jZ7Jy827+dWk/2jSqXSbHVRIREakA/vzRIv63cAP3/6Q7gw5tWmbHVRIRESnn/jNrLU98tpShR7bj8oEdyvTYSiIiIuXY3DUZ3Dl+FqkdGvPgOYlvSM9PSUREpJzatHMvI0ZNpXGdGjx5ef+E3JFelIT34isiIqVvX3YuN740jc279jH+hsE0r18zkjiUREREyqEH/jOPKSu28vdL+nBE24aRxaHqLBGRcualSSt5ZfIqbjj+UM7tE+3YfUoiIiLlyKRlm3ngnXmc2LU5d57WNepwlERERMqLtK27uenl6bRvWoe/D+1L1SpleyVWQZRERETKgV17s7lu1DSycnJ55spUGtSqHnVIgJKIiEjSy811bh83k0Xrt/P40L4c2rxe1CF9T0lERCTJPfLhQibMT+e3Z3fnhK4tog5nP0oiIiJJ7LWpq/n358u4fGB7rhqcEnU4P5KwJGJmz5vZBjObG1M2zsxmho8VeWOvm1mKmWXGLHsqUXGJiJQXk5dt5t4353BM52bc/5MeZd6lSTwSebPhi8A/gVF5Be5+cd60mf0ZyIhZf6m790lgPCIi5cbKzbu44aVptGtSh39d1o/qVZOz4ihhScTdvzCzlIKWWZBOLwJOStTxRUTKq4zMLK4ZORUHnh82gIa1k+NKrIJEldqOBdLd/buYso5mNsPMPjezYyOKS0QkUtk5udzyynRWbNrFk5f1J6VZ3ahDKlRUfWcNBcbEzK8D2rv7ZjPrD7xlZj3cfXv+Dc1sBDACoH37xA5ALyJS1n737ny+/G4Tj/zsiDIdXKq4yvxMxMyqAT8FxuWVufted98cTk8DlgKHFbS9uz/t7qnuntq8efOyCFlEpEyMmriCURNXct2xHbl4QPn4kRxFddbJwEJ3T8srMLPmZlY1nO4EdAGWRRCbiEgkvli8kQf/M5+TD2/B3WccHnU4cUvkJb5jgIlAVzNLM7NrwkWXsH9VFsBxwGwzmwWMB25w9y2Jik1EJJksWr+Dm1+eTpcW9fj7JcnRJ1a8Enl11tADlF9VQNnrwOuJikVEJFmlb9/D8Be+pXaNqjx/1QDq1ixfwzyVr2hFRCqQXXuzuWbkFLZlZvHq9YNo3ah21CEdtOS8e0VEpILLzsnl52NmMH/tdv51aT96toludMKS0JmIiEgZc3ce/M98/rdwA78/rycndkuuThUPhs5ERETK2HNfLWf0pJVcf1wnLh/YIepwSkRJRESkDP13zjoefn8BZx7RkrtO7xZ1OCWmJCIiUkamrdzK7eNm0rddI/5yUR+qlKNLeQ9ESUREpAys3LyL60ZNpWXDWjxzZSq1qleNOqRSoSQiIpJgW3ftY/gLU8h154WrBtC0Xs2oQyo1ujpLRCSB9mTlMGL0VNK2ZfLKtUfRKYnGRy8NOhMREUmQnFzn9rEzmbJiK3++sDepKU2iDqnUKYmIiCRAcC/IPD6Yt57fnt2dn/RuHXVICaEkIiKSAE98tpRRE4N7Qa45pmPU4SSMkoiISCl7bepqHvtwEef1aV0h7gUpjJKIiEgp+nThBu5+Yw7HdmnGoxf0rhD3ghRGSUREpJTMXL2Nm16ezuGt6vPk5f2pUa3if8VW/GcoIlIGlm3cydUvTqF5/Zq8cNWR1Ctn44IUl5KIiEgJbdixh2EvfAvAyKuPpHn9inMzYVGURERESmDHniyGvzCFTTv28cJVA+jYrG7UIZWpynG+JSKSAHuzc7jxpeksXL+DZ4el0rtdo6hDKnMJOxMxs+fNbIOZzY0pe8DM1pjZzPBxZsyye8xsiZktMrPTEhWXiEhpyMl1fjFuJl8t2cQjP+vFiV3L78BSJZHI6qwXgdMLKP+ru/cJH+8DmFl34BKgR7jNE2ZWMbq4FJEKx9257805vD9nPb8563Au6N826pAik7Ak4u5fAFviXP1cYKy773X35cAS4MhExSYiUhKPfLCIsVNWc8uJnbn22E5RhxOpKBrWbzGz2WF1V+OwrA2wOmadtLDsR8xshJlNNbOpGzduTHSsIiL7eerzpTz1+VIuO6o9d5x6WNThRK6sk8iTwKFAH2Ad8OewvKBbOr2gHbj70+6e6u6pzZs3T0yUIiIFGPvtKv7434Wc3asVvzu3J2YV+270eJRpEnH3dHfPcfdc4Bl+qLJKA9rFrNoWWFuWsYmIFOa/c9Zx75tzOP6w5vzloj5UreDdmcSrTJOImbWKmT0fyLty6x3gEjOraWYdgS7At2UZm4jIgXz53UZuGzuTvu0b8+Tl/SpFdybxSth9ImY2BjgBaGZmacD9wAlm1oegqmoFcD2Au88zs1eB+UA2cLO75yQqNhGReM1YtZXrR0+jU/O6PD9sAHVq6Pa6WOZeYNPDj1cMGsFbA5nAirBKKlKpqak+derUqMMQkQpqcfoOLnxqIg1rV2f8DYNo0aBW1CGVCjOb5u6ppbGvQlOqmTUEbgaGAjWAjUAt4BAzmwQ84e6flkYgIiLJZOXmXVz+7GRqVqvCS9ccVWESSGkr6rxsPDAKONbdt8UuMLP+wBVm1sndn0tUgCIiZW3NtkwufWYy+3JyGTdiEO2b1ok6pKRVaBJx91MKWTYNmFbqEYmIRGjD9j1c9swktmdm8cp1A+nasn7UISW1uFuIzKwN0CF2m/CudBGRCmHLrn1c9uxkNuzYy+hrjuSItg2jDinpxZVEzOwR4GKCq6fyrppyQElERCqEjMwsrnhuMqu27OaF4QPo36FJ1CGVC/GeiZwHdHX3vYkMRkQkCjv3ZnPVC9+yOH0HT1+ZyuBDm0UdUrkR7x0zy4DqiQxERCQKe7JyuHbkFGanZfD40L6Vtkv34irqEt/HCaqtdgMzzewT4PuzEXe/NbHhiYgkzt7sHK4fPY3Jy7fwt4v7cHrPVkVvJPspqjor706+aQRdk4iIVAhZObncOmYGny/eyCM/O4Jz+xTYcbgUoagkcjTwX+Bjd99RBvGIiCRcTq5zx6uz+HBeOg/8pDsXD2gfdUjlVlFtIs8DvYH3zewTM7vLzHqXQVwiIgkRJJCZvDNrLXed3o2rju4YdUjlWlE3G04CJgEPmFlT4FTgDjPrBUwHPnD3VxMfpohIyeXkOne+Nou3Zq7lztO6cuMJh0YdUrkX982G7r4ZGBM+8ro9KWgMdRGRpJOT69w5fhZvzFjDr049jJtP7Bx1SBVCvDcbNgKuBFLY/451XZ0lIkkvN9e56/XZvDF9Db885TBuOalL1CFVGPGeibxPUK01B4i8C3gRkXjl5jp3vzGb8dPSuP3kLtw6RAmkNMWbRGq5+y8TGomISCnLzXXufXMOr05N49YhXbj95MOiDqnCifeO9dFmdp2ZtTKzJnmPhEYmIlICubnOfW/NYeyU1dx6Umd+cbLOQBIh3jORfcBjwH0Ed7AT/u2UiKBEREoiN9f5zdtzGfPtam45sTO/OOUwzCzqsCqkeM9Efgl0dvcUd+8YPgpNIGb2vJltMLO5MWWPmdlCM5ttZm+GDfaYWYqZZZrZzPDxVPGfkohUZrm5zm/fnssrk1dx0wmHcsepSiCJFG8SmUfQf9bBeJEfXwI8Aejp7r2AxcA9McuWunuf8HHDQR5LRIScsBH95cmruPGEQ7nztK5KIAkWb3VWDkEHjJ8SZweM7v6FmaXkK/soZnYScEHckYqIFCI7J5c7x8/mzRlruG1IF24/uYsSSBmIN4m8FT5K09XAuJj5jmY2A9gO/MbdvyxoIzMbAYwAaN9e/d2ISNCZ4u3jZvLe7HXceVpX3UhYhuJKIu4+sjQPamb3AdnAy2HROqC9u28O74R/y8x6uPv2AmJ5GngaIDU11fMvF5HKZW92Dj9/ZQYfzU/nvjMP57rjdL1PWYqrTcTMzjazGWa2xcy2m9kOM/vRF3yc+xoGnA1c5u4O4O57w25VcPdpwFJAF3SLSKH2ZOVww+hpfDQ/6I1XCaTsxVud9Tfgp8CcvC/+4jCz04G7gOPdfXdMeXNgi7vnmFknoAvBaIoiIgXK3JfDiNFT+fK7Tfy/84/g0qNUvR2FeJPIamDuwSQQMxsDnAA0M7M04H6Cq7FqAhPCBq9J4ZVYxwG/M7Nsgkb8G9x9S9zPQkQqlV17s7lm5BQmL9/Coxf04qLUdlGHVGnFm0R+TTCmyOfsf3XWXw60gbsPLaD4uQOs+zrwepyxiEgltmNPFsNfmML0VVv560V9OK+vRiSMUrxJ5GFgJ1ALqJG4cEREDmzLrn1c9cK3zF+7nceH9uOsXhoTPWrxJpEm7n5qQiMRESlE+vY9XP7sZFZu2c1Tl/fn5O6HRB2SEP8d6x+bmZKIiERi1ebdXPDUN6zdlsmLwwcogSSReM9EbgZ+bWZ7gSzAAHf3BgmLTEQEWJy+g8ufncy+nFxevm4gfdo1ijokiRHvzYb1Ex2IiEh+s1ZvY9gL31KjahXGjRhE15b6Kko2hVZn5e/7qoDlZmZtSzMgERGAiUs3c+kzk6hfqxqv3aAEkqyKOhN5zMyqAG8D04CNBFdodQZOBIYQ3P+RlsggRaRy+d/CdG58aTrtm9Rh9DVH0bJhrahDkgMoNIm4+4Vm1h24jKDDxFYEXcIvIBh3/WF335PwKEWk0nh75hrueHUW3Vs34MXhR9Kkru4qSGZFtom4+3yCEQ1FRBJq9KSV/N/bczkypQnPDkulfq3qUYckRYj36iwRkYRxd/46YTH/+N8STj68Bf+8tB+1qleNOiyJg5KIiEQqOyeX34bjoV+c2o6Hz+9Jtarx3sImUVMSEZHI7MnK4edjZjBhfjq3nNhZ46GXQ3ElEQve1cuATu7+OzNrD7R0928TGp2IVFgZu7O4dtQUpq7cyoPn9GDY4JSoQ5JiiPec8QlgEJDXM+8O4F8JiUhEKrz1GXu46N8TmbU6g8eH9lUCKcfirc46yt37hWOg4+5bzUzX3YnIQVuyYSfDnv+WjMwsXhw+gMGdm0UdkpRAvEkky8yqAg7fj0SYm7CoRKRCmr5qK1e/OIVqVaowdsRAerZpGHVIUkLxVmf9A3gTaGFmDwNfAX9IWFQiUuFMmJ/Opc9MomHt6rxx42AlkAoi3g4YXzazaQTdnBhwnrsvSGhkIlJhvPj1ch58dz692jTk2WEDaF6/ZtQhSSmJ60zEzEa7+0J3/5e7/9PdF5jZ6Di2e97MNpjZ3JiyJmY2wcy+C/82DsvNzP5hZkvMbLaZ9Sv+0xKRZJCb6/z+3fk88J/5nHz4IYwZMVAJpIKJtzqrR+xM2D7SP47tXgROz1d2N/CJu3cBPgnnAc4AuoSPEcCTccYmIkloT1YON78ynWe/Ws5Vg1N46vL+1KmhW9MqmqK6gr/HzHYAvcxsu5ntCOc3EPTsWyh3/wLYkq/4XGBkOD0SOC+mfJQHJgGNzEwDKIuUQ5t37uXSZybxwbz1/Pbs7jxwTg+qVtFNhBVRoUnE3f8QDkj1mLs3cPf64aOpu99TzGMe4u7rwv2vA1qE5W2A1THrpYVl+zGzEWY21cymbty4sZghiEiiLN+0i58++Q3z1m7nycv6cc0xHaMOSRIo3ob1e8K2iy4E44nklX9RirEU9DPFC4jlaeBpgNTU1B8tF5HoTF2xhetGTcXMGDNiIP3aN446JEmweLs9uRa4DWgLzAQGAhOBk4pxzHQza+Xu68Lqqg1heRrQLma9tsDaYuxfRCLw3ux1/OLVmbRpVJsXhw+gQ9O6UYckZSDehvXbgAHASnc/EehLMMphcbwDDAunh/FD28o7wJXhVVoDgYy8ai8RSV7uzuOffMfNr0ynV5uGvHHjYCWQSiTeSyX2uPseM8PMarr7QjPrWtRGZjYGOAFoZmZpBEPp/hF41cyuAVYBF4arvw+cCSwhGD1x+ME9FREpa3uycrjr9dm8PXMt5/dtwx9+eoTGAalk4k0iaWbWCHgLmGBmW4mjqsndhx5g0ZAC1nXg5jjjEZGIbdixhxGjpjFz9TbuPK0rN51wqLpxr4TibVg/P5x8wMw+BRoCHyQsKhFJavPXbufakVPYujuLpy7vx+k9dTV+ZVVkEjGzKsBsd+8J4O6fJzwqEUlaE+anc9vYGTSoVZ3XbhikPrAquSKTiLvnmtksM2vv7qvKIigRST7uzr+/WMYjHyzkiDYNeebKVA5pUKvoDaVCi7dNpBUwz8y+BXblFbr7OQmJSkSSyt7sHO57cy7jp6VxVq9W/OmC3tSuoQZ0iT+JPJjQKEQkaW3YsYebXprO1JVbuW1IF24b0oUq6sJEQvE2rKsdRKQSmrV6G9ePnsa2zH38Y2hfzundOuqQJMmoS00RKdD4aWnc++Ycmteryes3DqZHazWgy48piYjIfrJycnn4vQW8+M0KBh/alH9e2o8mdWtEHZYkKSUREfne5p17ueWVGUxctplrjunIPWd0o1rVeHtHksqo0CRiZnMooCddgh533d17JSQqESlzc9dkcP3oaWzcuZe/XNSbn/ZrG3VIUg4UdSZydplEISKRenvmGu56fTaN69Rg/A2D6NW2UdQhSTlRaBJx95VlFYiIlL2snFwe+e9Cnv1qOUemNOFfl/XTGOhyUIqqztpB4dVZDRISlYgkXPr2PdzyynSmrNjKsEEduO+s7tSopvYPOThFnYnUL6tARKTsTFy6mZ+PmcGuvdn8/ZI+nNvnRyNRi8TloK7OMrMW7D88rvrSEilH8vq/evSDhaQ0q8sr1x3FYYfot6IUX7zD454D/BloTTCcbQdgAdAjcaGJSGnavieLX706i4/mp3NWr1Y88rNe1Kupq/ylZOL9BD1EMK76x+7e18xOBA404JSIJJn5a7dz08vTSNuayf+d3Z3hR6doACkpFfG2omW5+2agiplVcfdPgT4JjEtESsn4aWmc/8TXZGblMHbEQK4+pqMSiJSaeM9EtplZPeAL4GUz2wBkF+eA4djs42KKOgH/BzQCrgM2huX3uvv7xTmGiMDufdnc//Y8XpuWxqBOTfnH0L66fFdKXVGX+NZ0973AucAe4BfAZQTD4/6uOAd090WEZzFmVhVYA7wJDAf+6u5/Ks5+ReQHC9dv55ZXZrB0405+flJnbhvSRd2XSEIUdSYyEegHPOXuV4RlI0vx+EOApe6+UqfXIiXn7oydspoH3plH/VrVeemaozi6c7Oow5IKrKgkUsPMhgGDzeyn+Re6+xslPP4lwJiY+VvM7EpgKnCHu2/Nv4GZjQBGALRv376EhxepOHbsyeLeN+fyn1lrObZLM/5yUR9VX0nCmXtBN6SHC82OIai+ugh4J99id/eri31gsxrAWqCHu6eb2SHAJoI75B8CWhW1/9TUVJ86dWpxQxCpMOakZXDLmOmkbc3kl6ccxo3HH6rRB+WAzGyau6eWxr6KumP9K+ArM5vq7s+VxgFjnAFMd/f08FjpeQvM7Bng3VI+nkiF4+68+M0K/t/7C2hWrybjRgwkNaVJ1GFJJRLv8LjPmdlgICV2G3cfVYJjDyWmKsvMWrn7unD2fGBuCfYtUuFt2bWPX4+fzccL0jn58BY8dkFvGmvwKClj8d6xPho4FJgJ5ITFDhQriZhZHeAU4PqY4kfNrE+43xX5lolIjC8Wb+SO12aRsTuL357dnat186BEJN77RFKB7l5YA8pBcPfdQNN8ZVccYHURCe3JyuHRDxbx/NfLOeyQeowcfiTdW6szbYlOvElkLtASWFfUiiKSGAvXb+f2sTNZuH4HVw1O4e4zulGretWow5JKLt4k0gyYb2bfAnvzCt39nIREJSLfy80NGs//+MFCGtSqxgtXDeDEbi2iDksEiD+JPJDIIESkYBu27+FX42fzxeKNDOnWgkcu6EWzerr3Q5JHvFdnfZ7oQERkfx/OW889b8xh975sHjqvJ5cf1V6N55J0iuo76yt3P6bUse9vAAAXSElEQVSAYXI1PK5IgmTszuLB/8zjjRlr6NG6AX+/pA+dW2jgKElORd1seEz4V59gkTLw+eKN3DV+Nht37uXWIV245cTOGvdckpqGxxVJAjv3ZvPwewsY8+0qurSox9NX9qdX20ZRhyVSJA2PKxKxiUs3c+f4WazZlsn1x3XiF6ccpkt3pdzQ8LgiEcncl8OjHy7kha9XkNK0Dq9dP0j9Xkm5E28SyXL3zWb2/fC4ZvZIQiMTqcCmrtjCr8fPZtmmXQwb1IG7zuhGnRoHVbsskhTKfHhckcps595sHvtgIaMmraR1w9q8cu1RDNagUVKOxZtEzgUyKYXhcUUqq08XbeC+N+awbvserhqcwq9O7Urdmjr7kPIt3psNd4WTucDIcGz0S4CXExWYSEWxddc+Hnp3Pm/MWEPnFvUYf8Ng+ndoHHVYIqWiqJsNGwA3A20IRjacEM7fSdAtvJKIyAG4O+/NWcf9b88jIzOLW0/qzM0ndaZmNV15JRVHUWcio4GtwETgWoLkUQM4191nJjg2kXJrfcYefvPWXD5ekE6vtg156dqjOLyVOniQiqeoJNLJ3Y8AMLNnCcZAb+/uOxIemUg5lJPrvDRpJX/6cBFZubncd+bhDD86hWpVdde5VExFJZGsvAl3zzGz5UogIgWbk5bBvW/OYc6aDI7t0oyHzu1JSrO6UYclklBFJZHeZrY9nDagdjivDhhFQtv3ZPGXjxYzauIKmtaryeND+3J2r1bqcVcqhaI6YExYC6CZrQB2EIzZnu3uqWbWBBgHpBCMs36Ru29NVAwiJZHXcP67/8xn4869XDGwA3ec2pWGtatHHZpImYn6IvUT3X1TzPzdwCfu/kczuzucvyua0EQObOXmXfz27Xl8sXgjPVo34JkrU+ndTh0mSuUTdRLJ71zghHB6JPAZSiKSRPZk5fDvz5fxxGdLqF61Cvf/pDtXDOyghnOptKJMIg58ZGYO/NvdnwYOcfd1AO6+Lux6XiRy7s6H89L5/XvzSduayVlHtOK3Z3enZcNaRW8sUoFFmUSOdve1YaKYYGYL49nIzEYAIwDat2+fyPhEAFiyYQcPvDOfr5Zsoush9XnluqMYfKj6uxKBCJOIu68N/24wszeBI4F0M2sVnoW0Ihi7JP92TwNPA6Smpnr+5SKlZfueLP7+8XeM/GYFdWpU5YGfdOdyVV2J7CeSJGJmdYEq7r4jnD6VoEPHd4BhwB/Dv29HEZ9Ubrm5zvjpaTz6wUI279rHxantuPO0rjStVzPq0ESSTlRnIocAb4bX0VcDXnH3D8xsCvCqmV0DrAIujCg+qaSmrdzK796dz6zV2+jbvhHPXzVAw9SKFCKSJOLuy4DeBZRvBoaUfURS2a3avJtHPljIe3PW0bx+Tf58YW/O79uGKlV0w6BIYZLtEl+RMpWxO4t/fvodI79ZSZUqcOuQLlx/XCeN8yESJ/2nSKWUlZPLS5NW8vdPviMjM4sL+rXljlO76pJdkYOkJCKVirvz0fx0/vjfhSzftIujOzfl3jMPp0frhlGHJlIuKYlIpTFt5RYe+WAR3y7fQucW9Xj+qlRO7NpCHSWKlICSiFR4C9dv508fLuLjBRtoVq8mD53Xk6ED2ul+D5FSoCQiFdaqzbv568eLeWvmGurVrMadp3Vl+NEp1Kmhj71IadF/k1Q4G3bs4Z//W8KYb1dRxYzrjzuUG47vRKM6NaIOTaTCURKRCiNjdxZPf7mU579aQVZOLhcPaMetQ7pwSANdcSWSKEoiUu5lZGbx/FfLef7r5ezYk805vVvzy1MO09C0ImVASUTKrfzJ4/QeLbl1SBe6t9aozSJlRUlEyh0lD5HkoSQi5YaSh0jyURKRpLdp515e+Ho5oyauZMeebE7rcQi3Dumiu8xFkoCSiCSt1Vt288yXyxg3ZTX7cnI5vUdLbjmps5KHSBJREpGkszh9B099tpS3Z62lisH5fdtw/fGHcmjzelGHJiL5KIlI0pi+aitPfLqUjxekU7t6Va4anMK1x3akVcPaUYcmIgegJCKRys11Plm4gWe/XMbk5VtoVKc6tw3pwlWDU2hcV3eYiyQ7JRGJxO592YyflsYLX69g+aZdtG5Yi9+cdThDj2yvAaFEypEy/281s3bAKKAlkAs87e5/N7MHgOuAjeGq97r7+2UdnyTWuoxMRn6zkjHfriIjM4ve7Rrx+NC+nNGzpXrVFSmHovjJlw3c4e7Tzaw+MM3MJoTL/uruf4ogJkmw2WnbeO6r5bw3ex257pzesyXXHNORfu0bazwPkXKszJOIu68D1oXTO8xsAdCmrOOQxNubncMHc9czauJKpq3cSr2a1bhyUArDj06hXZM6UYcnIqUg0spnM0sB+gKTgaOBW8zsSmAqwdnK1uiik+Jasy2TVyavZNyU1WzauY+UpnX4zVmHc9GAdjSoVT3q8ESkFEWWRMysHvA6cLu7bzezJ4GHAA///hm4uoDtRgAjANq3b192AUuhcnOdr5duYtTElXyyIB2Ak7odwpWDOnBM52ZUqaIqK5GKyNy97A9qVh14F/jQ3f9SwPIU4F1371nYflJTU33q1KkJiVHis3XXPl6fnsbLk1exfNMumtatwcUD2nHpUe1p21hVViLJyMymuXtqaewriquzDHgOWBCbQMysVdheAnA+MLesY5P45J11jJ2ymgnz0tmXk0v/Do25bUgXzjiiJTWrVY06RBEpI1FUZx0NXAHMMbOZYdm9wFAz60NQnbUCuD6C2KQQa7dl8trUNF6dupo12zJpVKc6lw1sz8UD2tGtpXrSFamMorg66yugoApy3ROShPZl5/LxgnTGTVnNF99txB2O6dyMu8/oxindD6FWdZ11iFRmujVYfsTdmb5qG2/NWMO7s9eydXcWrRrW4ucndubC1Ha6PFdEvqckIt9bsWkXb85Yw1sz17By825qVqvCqT1a8tN+bTiuS3Oq6gorEclHSaSS27JrH+/NXssbM9YwY9U2zGBQp6bccmJnTu/Zkvq6r0NECqEkUgllZGYxYX46789ZxxeLN5Kd63RrWZ97zujGOX1aq+t1EYmbkkglEZs4vvxuI1k5TptGtbnmmI6c17cNh7fS1VUicvCURCqwjMwsPp6fznv5Esfwozty5hGt6N22oTo/FJESURKpYNZlZPLx/HQmLNjAxKWblDhEJKGURMo5d2f+uu18PH8DExasZ+6a7QB0bFaX4Ud35IyeLenTrpESh4gkhJJIObQ3O4dvl2/h4/npfLxgA2u2ZWIG/do35q7Tg5sAO7eoF3WYIlIJKImUEys27eLzxRv5fPFGJi7dTGZWDjWrVeHYLs25dUhnTup2CM3r14w6TBGpZJREktSuvdlMXLr5+8SxastuADo0rcOFqW05rktzju7cjNo11O2IiERHSSRJ7M3OYeaqbUxctpmJSzczfdVWsnKcOjWqMqhTU649tiPHdWlOSrO6UYcqIvI9JZGIZOXkMjstg0nLNvPN0k1MW7mVPVm5mEGP1g24+uiOHH9Yc/qnNFbX6iKStJREysjufdnMWp3BtJVbmLJiK1NWbGH3vhwAurWsz9Aj2zOoU1OO6tiUhnXU1YiIlA9KIgmyLiOTqSu2Mm1l8Ji/bjs5ucEokl1a1ONn/doy6NCmHNWxCU3rqUFcRMonJZFSsH1PFnPXZDAnLYPZazKYsXIrazP2AFC7elV6t2vIjccfSv8OjenXvrHONESkwlASOUg792Yzb00Gc9ZkMDst+Lt8067vl7dtXJu+HRpzXYfG9O/QmMNbNaB61SoRRiwikjhKIgeQk+us2LyLRet3sHDddhau38Gi9B2s2rIbD2qlaN2wFke0bcgF/dtyRJuG9GzTkCZ1a0QbuIhIGar0SSQrJ5dVW3azfOMulm3ayaL1O1mUvp3v0neyNzsXgCoGKc3q0qN1A37Wry1HtG3IEW0a0kxtGSJSySVdEjGz04G/A1WBZ939jyXdZ3ZOLuk79rJy0y6WbdrF8pjHqi27v2/wBmhevybdWtbnioEd6NaqAd1a1qdzi3oaS1xEpABJlUTMrCrwL+AUIA2YYmbvuPv8wrbbk5XD2m2ZrNmWyZqtP/xNC/+u375nv0RRq3oVOjarR/dWDTjriFZ0bFaXjs3r0qlZXRrVUXWUiEi8kiqJAEcCS9x9GYCZjQXOBQpMIks37mTAwx+zccfe/cqrGLRqWJs2jWpzZMcmtGlUmzaNa9O+SR06Na/LIfVrUUXjhYuIlFiyJZE2wOqY+TTgqNgVzGwEMCKc3Tv1N6fMLWhHyxMSXtyaAZuiDeFHFFN8FFP8kjEuxRSfrqW1o2RLIgWdHvh+M+5PA08DmNlUd08ti8AORjLGpZjio5jil4xxKab4mNnU0tpXst3AkAa0i5lvC6yNKBYRESlCsiWRKUAXM+toZjWAS4B3Io5JREQOIKmqs9w928xuAT4kuMT3eXefV8gmT5dNZActGeNSTPFRTPFLxrgUU3xKLSZz96LXEhERKUCyVWeJiEg5oiQiIiLFVm6TiJmdbmaLzGyJmd2d4GM9b2YbzGxuTFkTM5tgZt+FfxuH5WZm/wjjmm1m/WK2GRau/52ZDSthTO3M7FMzW2Bm88zstqjjMrNaZvatmc0KY3owLO9oZpPD/Y8LL5rAzGqG80vC5Skx+7onLF9kZqcVN6aY/VU1sxlm9m4SxbTCzOaY2cy8Sy6T4HPVyMzGm9nC8LM1KOLPVNfw9cl7bDez25PgdfpF+Bmfa2Zjws9+pJ8pM7stjGeemd0eliX+dXL3cvcgaHRfCnQCagCzgO4JPN5xQD9gbkzZo8Dd4fTdwCPh9JnAfwnueRkITA7LmwDLwr+Nw+nGJYipFdAvnK4PLAa6RxlXuO964XR1YHJ4rFeBS8Lyp4Abw+mbgKfC6UuAceF09/A9rQl0DN/rqiV8D38JvAK8G84nQ0wrgGb5yqL+XI0Erg2nawCNoo4pJraqwHqgQ8Sf8zYE9zPXjvksXRXlZwroCcwF6hBcMPUx0KUsXqcSvalRPYBBwIcx8/cA9yT4mCnsn0QWAa3C6VbAonD638DQ/OsBQ4F/x5Tvt14pxPc2QZ9jSRFX+GGeTtDjwCagWv73juAqvEHhdLVwPcv/fsauV8xY2gKfACcB74bHiDSmcB8r+HESiez9AxoQfDlassSUL45Tga+jjokfetZoEn5G3gVOi/IzBVxI0GFt3vxvgV+XxetUXquzCuoepU0Zx3CIu68DCP+2KCK2hMUcnh73JfjlH2lcYbXRTGADMIHg19U2d88uYP/fHztcngE0Le2YgL8R/EPlhvNNkyAmCHpj+MjMplnQnQ9E+/51AjYCL1hQ9fesmdWNOKZYlwBjwunIYnL3NcCfgFXAOoLPyDSi/UzNBY4zs6ZmVofgTKMdZfA6ldckUmT3KBE6UGwJidnM6gGvA7e7+/ao43L3HHfvQ/Dr/0jg8EL2n/CYzOxsYIO7T4stjjKmGEe7ez/gDOBmMzuukHXLIq5qBNW2T7p7X2AXQRVIlDEFBwraF84BXitq1UTHFLYrnEtQBdUaqEvwHh5o/wmPyd0XAI8Q/HD7gKCaLLuQTUotpvKaRJKhe5R0M2sFEP7dUERspR6zmVUnSCAvu/sbyRIXgLtvAz4jqG9tZGZ5N7bG7v/7Y4fLGwJbSjmmo4FzzGwFMJagSutvEccEgLuvDf9uAN4kSLpRvn9pQJq7Tw7nxxMklWT4TJ0BTHf39HA+yphOBpa7+0Z3zwLeAAYT8WfK3Z9z937ufly4/+8oi9eppPWUUTwIfjEtI/glkNew3iPBx0xh/zaRx9i/werRcPos9m+w+jYsb0JQ39w4fCwHmpQgHgNGAX/LVx5ZXEBzoFE4XRv4Ejib4NdjbIPjTeH0zezf4PhqON2D/Rscl1HCRuxwvyfwQ8N6pDER/HqtHzP9DXB6EnyuvgS6htMPhPFEGlO4z7HA8CT5nB8FzCNo9zOCixF+ngSfqRbh3/bAwvB5Jvx1KtE/ZZQPgjq/xQR17vcl+FhjCOo+swgy9TUEdZqfEGT7T/Je6PBN+VcY1xwgNWY/VwNLwsfwEsZ0DMFp5mxgZvg4M8q4gF7AjDCmucD/heWdgG/D/b8G1AzLa4XzS8LlnWL2dV8Y6yLgjFJ6H0/ghyQSaUzh8WeFj3l5n+Ek+Fz1AaaG7+Fb4RdJ1DHVATYDDWPKoo7pQYIv6rnAaIJEEPVn6kuCsZdmAUPK6nVStyciIlJs5bVNREREkoCSiIiIFJuSiIiIFJuSiIiIFJuSiIiIFJuSiIiIFJuSiEgEzKy7mV1lQZf+9cvgeOeZ2TNm9raZnVrMfXQys+fMbHxpxyfll5KIlDkzO9/M3My6HeR23xTzeDuLs13M9o3M7KaS7KMA1Qnucj4fKDA+M6ttZp+bWdWSHszd33L36wi6LL+4mPtY5u7X5Iuxhpl9EdPdh1QySiIShaHAVwRdQOwnHCynSkHz7j44EcHkP2YBGhGMCVGa2gEvENwVfKAzkauBN9w952B2bGZHmNm7+R55vbf+huBO5eJuvx9330dwJ3SxEpNUACW5zV4PPQ72AdQD1gCHAQvDshRgAfAEQbcpx+eb7xCut5Ogp9KbYvb3AHBHOP0WQZfc84ARMevsLCCO/MfsUMj2Y4FMgq5lHgvLLifowmImwZgLBfZ5RJB85gIrgZ/nW9ayiNfqGyAlZv414J8ECXglQdc3owi6/3muiH1Z+NqdnK98WPicZwNfxvkejs833xt4P+rPlh7RPCIPQI/K9Qi/fJ8Lp78h6CU2hWCsj4Fh+X7zMdvuJBg35fOYsvlA+3A6r1+g2uEXd9O87QqI40fHKGT7FPbvfPNw4D9A9XD+CeDKAo7xM4LeeasTDPjz/aBFcbxONYD1+coWAr8Mp3/HDwMJ1QC2EvbVdID93Romi6eAG8Ky+uHrVyOcb1RETE3D7Zey/2BKVYGNUX+29IjmoXpMKWtDCbpih+AX/lCC6pWV7j4pZr388wC4+wwza2FmrQl6Dd7q7qvCxbea2fnhdDuC4UE3FxJL/mPEu/0QoD8wxcwgSDobCljvVuA6D7oLX2dmWcRfhdwM2JY3Y2a1CKrV8l67TIJkvC5cvhvYd6Cdufs/gH/kK84JY/+zmY1096mFBeTum4EbCijPMbN9Zlbf3XcU+cykQlESkTJjZk0JxvPoaWZO8AvWCX7J78q3ev75WOOBC4CWBIkIMzuBYJyHQe6+28w+I+g9tTDfH+MgtzdgpLvfc6Adh2O99HL3xeF8K2CTB20I8cjMd/weBONp5I3O2Bt4Mtx3W2Ctux9Ub6rh8+wJ/AR42syedfcnDmYfMWoCe4q5rZRjaliXsnQBMMrdO7h7iru3IxivoO1B7mcsQaP8BQQJBYKBfraGX4zdCMZIOBiFbb+D/Ru/PwEuyGtsNrMmZtYh3/66Aw3Dy2KrAH8AHs9baGbjzOyOAwXj7luBquEZCMARBF185+lF0I4BQUKZzUEysy7uvsvdxxKME15U0j3QfpoSVGdlFWd7Kd+URKQsDSVoI4j1OnDvwezE3ecRfKmvyavOIRgStJqZzQYeAn5UFVaEA24fVuN8bWZzzewxd59PcJXTR+H6EwjaJmL1BV4mGItmNrDK3Z8GMLNzCb60exUR00cEjecQJJGZ4fa1gNphooH9E8rBuM/MFpnZdIJBkYp7FnIi8H4xt5VyTuOJiCSAmf0NmOju4/KV1wKedffLzewtdz+vkH30JWhIvyLB4ZaImb1B0NC+KOpYpOzpTEQkMfoQnjnkcydQz8yeAnqYWe0D7cDdZwCflsbNholiZjWAt5RAKi+diYiUETNrD9zv4V3fZnY/8IG7T442MpHiUxIREZFiU3WWiIgUm5KIiIgUm5KIiIgUm5KIiIgUm5KIiIgUm5KIiIgUm5KIiIgUm5KIiIgUm5KIiIgU2/8HhkYa0xQnN5IAAAAASUVORK5CYII=\n",
            "text/plain": "<matplotlib.figure.Figure at 0x7f1ca716abe0>"
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "We have now calculated the arrival rate as function of rainfall rate,  whereas we need to know the rainfall rate as function of the measured arrival rate. Although the inverse of the upper incomplete Gamma function is available in Python, we prefer to use a polynomial approximation for its ease-of-use. We fit a 3rd order polynomial through this curve. The polynomial is forced through (0,0), which, strictly speaking, need not be true because there may be very light rainfall with all drops below the minimal detectable dropsize."
    },
    {
      "metadata": {
        "trusted": false
      },
      "cell_type": "code",
      "source": "weight=np.zeros(nr)+1  #Equal weights to all points\nweight[0]=10000        #Huge weight to (0,0)\n\nfit=np.polyfit(fw[0:nr,0],fw[0:nr,1],3,w=weight)   #Fit 3rd order polynomial\n\ny = np.polyval(fit,fw[0:nr,0])     #Calculate values of polynomial for given arrival rates\n\n#Plot results\nplt.plot(fw[0:nr,0],fw[0:nr,1], label='Forward model', linewidth=3, linestyle='--')\nplt.plot(fw[0:nr,0],y, label='Polynomial', color='red')\nplt.axis([0,9000,0,200])\nplt.title(r'Rainfall rate vs Arrival rate')\nplt.xlabel(r'Arrival rate $\\rho_A^*\\ (m^{-2}s^{-1}$)')\nplt.ylabel('Rainfall rate R (mm/h)')\nplt.legend()\n\n",
      "execution_count": 5,
      "outputs": [
        {
          "data": {
            "text/plain": "<matplotlib.legend.Legend at 0x7f1ca67d8240>"
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEdCAYAAADNU1r0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzs3XWYlOXXwPHv2aW7WxpBOlZKkJBUaVCUNhAE6RYVUBSVkhAE6VKUkFeRFAQUUEpAGiQWkG5Y2DjvH/PAb8Fld1h2djbO57rmYuaeJ84Ec/Z+7hJVxRhjjIkMH28HYIwxJvayJGKMMSbSLIkYY4yJNEsixhhjIs2SiDHGmEizJGKMMSbSLImYKCciOUXkuoj4url9IxE54exTKoJt24rIhlCPVUTyP27McZGITBSR96LgONNF5KOoiMnEPZZETJhE5KiI3HJ+2P91fkhSuLOvqh5X1RSqGuzm6YYDnZ19tkc+6scjIlVFxN9b53+QiKwVkUsikjgy+6tqB1X9MKrjelzOd6uGt+MwUcOSiAlPPVVNAZQESgH9PXSeXMDfHjo2AOISa77vIpIbqAwoUD+Cbf9T43O3FhjVRCSBN85rvCfW/Kcy3qOq/wLLcSUTAETkBRHZLiJXnUtRg0I9l9u5zJTAebxWRD4Ukd9E5JqIrBCRDCKSWESuA77AXyJy2Nm+n4gcdrbdIyKNIhO3c96hIvIbcBPIKyLtRGSvc+wjIvKWs21y4Gcgm1P7ui4i2UTEJ1Q8F0Rkvoike8j59orIi6EeJxCR8yJSWkSSiMhs5xiXReRPEckcTvitgU3AdKDNA+eZLiITRGSpiNwAqoVT9lFEsTmPv3NqnFdEZJ2IFHHzPW7rfK6jROQiMEhE8onIL85rPS8ic0QkjbP9LCAn8H/Oe9zHKS8vIr87781fIlLVnfMb77MkYiIkIjmAusChUMU3cP3QpQFeADqKSMNwDvMq0A7IBCQCeqnqbaemA1BCVfM59w/j+is8NTAYmC0iWSMZfiugPZASOAacBV4EUjnxjBKR0qp6w3mNp5zLailU9RTQBWgIVAGyAZeA8Q851zzglVCPawPnVXUbrkSQGngCSA90AG6FE3drYI5zqx1GwnkVGOq8rg3hlLkTG7gSaAFcn88257zuKgcccfYdCgjwCa736ylcr3kQgKq2Ao7j1HJV9TMRyQ78BHwEpAN6AQtEJOMjxGC8xJKICc9iEbkGnMD14/vB3SdUda2q7lLVEFXdietHqko4x5qmqgdU9RYwn1C1mgep6neqeso59rfAQaBsJF/DdFX9W1WDVDVQVX9S1cPq8iuwAlfCepi3gHdV1V9Vb+P6MWz6kMs2c4H6IpLMefyqUwYQiCt55FfVYFXdqqpXwzqhiFTCdYlvvqpuxZVUX31gsx9U9TfnPQoIp8yd2FDVqap6LdRrLCEiqR/+ttznlKqOdd7jW6p6SFVXOn8knANGEv53oyWwVFWXOrGvBLYAz7t5fuNFlkRMeBqqakqgKlAIyHD3CREpJyJrROSciFzB9Zd1hrAPA8C/oe7fBB7aSC8irUVkh3Np4zJQNIJjh+fEA8euKyKbROSic+znIzh2LmBRqFj2AsHAfy5Fqeoh5/l6zo91ff73Qz0L1yXBb0TklIh8JiIJH3LONsAKVT3vPJ7LA5e0Hnxd4ZRFGJuI+IrIMOeS3VXgqLObu+/5g+9xJhH5RkROOsebHcGxcgHN7r7HzvtcCYhs7dNEI2sEMxFS1V9FZDquXlR3L1nNBcYBdVU1QERGE/kf+ntEJBcwGXgO2KiqwSKyA9clksi4N021uHo5LcB1qegHVQ0UkcWhjh3WlNYngNdU9Tc3z3f3spEPsMf58UZVA3FdmhssrkbzpcB+YEronUUkKfAS4CsidxNvYiCNiJRQ1b/CiTWiKbnDjA1XraQBUANXAkmN67Kdu+/5g+f9xCkrrqoXnMuc48LZ/gQwS1XfdPN8Jgaxmohx12igpojcvQyVErjoJJCy/PdyS2Qlx/Ujcw5ARNrhqolEhUS4fpDPAUEiUheoFer5M0D6By7jTASGOskNEckoIg3COcc3zjE7EupykYhUE5Fi4uo1dRXX5a2wukA3dMoL47rkVxJXu8J6XMnvcYQZG67P8jZwAUgGfPyY50kJXAcuO+0dvR94/gyQN9Tj2bhqSLWdWlEScXW3zvGYcZhoYEnEuMW5tj0TuDt47W1giNNm8j6udo6oOM8eYASwEdePTTHA3VpARMe+hquhfD6uv7RfBZaEen4frr/WjziXVbIBXzjbrHBe6yZcDckPO8dpJ/aKwLehnsoCfI8rgewFfsX14/mgNrjaj46r6r93b7j+km/xkLYYt4QT20xcnQ5OAntwvcbHMRgoDVzB1WC+8IHnPwEGOu9xL1U9gasmNABXgj+BK/HY71MsILYolTHGmMiyTG+MMSbSPJZEROQJp/fOXhH5W0S6OuXpRGSliBx0/k3rlIuIjBGRQyKy8+4gKGOMMTGXJ2siQUBPVX0KKA90EpHCQD9gtaoWAFY7j8E10KuAc2sPTPBgbMYYY6KAx5KIqp6+OxrWadDcC2TH1YA2w9lsBv/rMtoAmOkMAtuEq0uj9RM3xpgYLFraRJx+8aWAzUBmp5fI3d4imZzNsnP/oCV/p8wYY0wM5fHBhuKaPnwB0E1Vr4o8dPxSWE/8p+uYiLTHdbmL5MmTlylUqFBUhWqMMXHSmasB+Jw6RcYblwh5IifbTxw/r6pRMjeZR5OIM63DAmCOqt7tK35GRLKq6mnnctVZp9wf10Rtd+UATj14TFWdBEwC8PPz0y1btngsfmOMie1W7jnDwgGjmXBwGHNK1uHSiDG889yTx6Lq+J7snSW4pnTYq6ojQz21hP/NA9QG+CFUeWunl1Z54Mrdy17GGGMe3eFz1xk/djHDl45ma7ZCrGw/gA5V8kW84yPwZE3kGVzTcO9y5j4C14jUYcB8EXkd15TQzZznluKaDO8Qrgn62nkwNmOMidOu3w6ix6S1jP52MDcSJWVIm8FMb12OBL5RW3fwWBJR1Q08fAK358LYXoFOnorHGGPiC1Wlz7fb6TJtCDmunKV1y2EM7VSbtMkTRfm54twsvoGBgfj7+xMQ8OByCiYmS5IkCTly5CBhwofNjm6McdeEXw9T8KtRPHf4TwbW7Eizrs0pmt3d5WEeTZxLIv7+/qRMmZLcuXMTTk8wE4OoKhcuXMDf3588efJ4OxxjYrV1B86xfdxMJv8+j/nFapCg09s0Lu25CZHj3NxZAQEBpE+f3hJILCIipE+f3mqPxjymExdvMmLcEkb83wj+ylKAH958l3dfLOzRc8a5mghgCSQWss/MmMfne+M6Y+d/yJ0ECRnYahBT2pYnYRQ3pD8oztVEYgJfX19Klix573b06FFvh8TatWt58cUXPXqOo0ePUrRo+OtHubONMSYSQkLI1q0jT5z355uewxnUqS6ZUibx+GnjZE3E25ImTcqOHTsi3vABQUFBJEgQNR9JcHAwvr6+UXIsY0ws8NFHsHAhMnIknbu/EW2ntZpINAkICKBdu3YUK1aMUqVKsWbNGgCmT59Os2bNqFevHrVq1eLtt99myRLXYnuNGjXitddeA2DKlCkMHDgQgIYNG1KmTBmKFCnCpEmT7p0jRYoUvP/++5QrV46NGzeybNkyChUqRKVKlVi48MHF5bh3/oYNG1KvXj3y5MnDuHHjGDlyJKVKlaJ8+fJcvHgRgB07dlC+fHmKFy9Oo0aNuHTpEgBbt26lRIkSVKhQgfHjx987bnBwML179+bpp5+mePHifPXVV1H8jhpjAK4GBMKiRfDBB9C6NXTrFq3nj/M1kVErD/DF6oNubftK2Sf4pHHx+8r6L9zJvD/+Ny9k1+cK0L3mk+Ee59atW5Qs6VqKPE+ePCxatOjeD+yuXbvYt28ftWrV4sCBAwBs3LiRnTt3ki5dOr755hvWr19P/fr1OXnyJKdPuwbtb9iwgebNmwMwdepU0qVLx61bt3j66adp0qQJ6dOn58aNGxQtWpQhQ4YQEBBAgQIF+OWXX8ifPz8vv/zyQ+PdvXs327dvJyAggPz58/Ppp5+yfft2unfvzsyZM+nWrRutW7dm7NixVKlShffff5/BgwczevRo2rVrd6+8d+//LaU9ZcoUUqdOzZ9//snt27d55plnqFWrlrV9GBOF9v97jXcHz2HejJ4kLFcOvvoKovn/mNVEPODu5awdO3awaNEiwJUEWrVqBUChQoXIlSvXvSRSs2ZN0qVLB0DlypVZv349e/bsoXDhwmTOnJnTp0+zceNGKlasCMCYMWMoUaIE5cuX58SJExw86EqSvr6+NGnSBIB9+/aRJ08eChQogIjQsmXLh8ZbrVo1UqZMScaMGUmdOjX16tUDoFixYhw9epQrV65w+fJlqlSpAkCbNm1Yt27df8rvvj6AFStWMHPmTEqWLEm5cuW4cOHCvTiNMY/v8s079JqwmlHzBnHBNwmTug+HJJ5vA3lQnK+JxBThrWWfPHnye/ezZ8/OpUuXWLZsGc8++ywXL15k/vz5pEiRgpQpU7J27VpWrVrFxo0bSZYsGVWrVr3XNTZJkiT3tYO4+1d/4sSJ79338fG599jHx4egoKBwX9PDzqGqjB07ltq1a99XHhM6GRgT2wUFh9B15p+8O+09Ml2/SKu2nzOkSvGId/SAOJ9Eutd8MsLLT+H5pHHx/1ziioxnn32WOXPmUL16dQ4cOMDx48cpWLAg27Zt+8+2FSpUYPTo0fzyyy9cuHCBpk2b0rRpUwCuXLlC2rRpSZYsGfv27WPTpk1hnq9QoUL8888/HD58mHz58jFv3rxIx546dWrSpk3L+vXrqVy5MrNmzaJKlSqkSZOG1KlTs2HDBipVqsScOXPu7VO7dm0mTJhA9erVSZgwIQcOHCB7dlsexpioMHTpXmpM+pjyJ3bT9cWevNbtZQplSeWVWOJ8Eokp3n77bTp06ECxYsVIkCAB06dPv68GEFrlypVZsWIF+fPnJ1euXFy8eJHKlSsDUKdOHSZOnEjx4sUpWLAg5cuXD/MYSZIkYdKkSbzwwgtkyJCBSpUqsXv37kjHP2PGDDp06MDNmzfJmzcv06ZNA2DatGm89tprJEuW7L5axxtvvMHRo0cpXbo0qkrGjBlZvHhxpM9vjHGZv+UEgeO+pNX2pUwo15Q8XdtTp2gWr8Uj4V1mienCWk9k7969PPXUU16KyDwO++yMCd/WY5cY+e5kps8bwLo8pfnuvbF82aosPj6P1pguIltV1S8qYrKaiDHGxAL/XglgyJgfmbZwKMfSZGXca4OY3bzMIyeQqGZJxBhjYrhbd4J5Z/I6hs16H9+QYHq0HMK4t6qSPLH3f8Kti68xxsRw249d4M3Jg3jy/HG6NOxHv3fqkTN9Mm+HBVhNxBhjYryKM8bAgY0Mr9OBmt1aUTF/Bm+HdI8lEWOMiclmzYKhQ+H113l99GjSpgi7V6e3eOxylohMFZGzIrI7VNm3IrLDuR29u/a6iOQWkVuhnpvoqbiMMSbW2LAB3ngDqlWDL7+McQkEPNsmMh2oE7pAVV9W1ZKqWhJYAISeFfDw3edUtYMH4/K4u1PBFy1alGbNmnHz5s1wt0+RIkU0RXa/LVu20KVLl3C3iY4p5I0x9ztx8SZdhy4guGEjyJ0bFiyARFG/PnpU8FgSUdV1wMWwnhPXXBkvAZEfRh2D3Z07a/fu3SRKlIiJE2NmxcrPz48xY8Z4OwxjTCjXAgLpOnENnUd15/rN2+ycOBvSpvV2WA/lrd5ZlYEzqhp6Rr48IrJdRH4VkcpeiivKVa5cmUOHDgEwcuRIihYtStGiRRk9evR/tm3VqhU//PDDvcctWrRgyZIlTJ8+ncaNG1OnTh0KFChAnz597m0zb948ihUrRtGiRenbt++98hQpUtC3b1/KlClDjRo1+OOPP6hatSp58+a9N9V86FrGH3/8QcWKFSlVqhQVK1Zk//79Hnk/jDEPFxyi9Ji9hW6T3yPXpdN0avwut/Pk83ZY4fJWw/or3F8LOQ3kVNULIlIGWCwiRVT16oM7ikh7oD1Azpw5wz9Lt24QicWhwlWyJISRAMISFBTEzz//TJ06ddi6dSvTpk1j8+bNqCrlypWjSpUqlCpV6t72b7zxBqNGjaJBgwZcuXKF33//nRkzZjB79mx27NjB9u3bSZw4MQULFuSdd97B19eXvn37snXrVtKmTUutWrVYvHgxDRs25MaNG1StWpVPP/2URo0aMXDgQFauXMmePXto06YN9evXvy/WQoUKsW7dOhIkSMCqVasYMGAACxYsiNK3zhgTvo9/2kOV8R/y7NHt9K7blUbdW/B07nTeDitc0Z5ERCQB0Bgoc7dMVW8Dt537W0XkMPAksOXB/VV1EjAJXNOeREfMjyr0eiKVK1fm9ddfZ8KECTRq1OjejL2NGzdm/fr19yWRKlWq0KlTJ86ePcvChQtp0qTJvZUOn3vuOVKnTg1A4cKFOXbsGBcuXKBq1apkzJgRcNVc1q1bR8OGDUmUKBF16riapIoVK0bixIlJmDDhvendH3TlyhXatGnDwYMHERECAwM99v4YY/5r7ubjBH8xhpY7fmZiuSZk6NKBJmVyeDusCHmjJlID2Keq/ncLRCQjcFFVg0UkL1AAOPLYZ3KzxhDVwloe1905ylq1asWcOXP45ptvmDp16r3y0JM1+vr6EhQUFO4xEyZMeG+adnemd3/vvfeoVq0aixYt4ujRo1StWtWteI0xj++3Q+f5ZcRUvvrla5YXKM+Ojn34slZBb4flFk928Z0HbAQKioi/iLzuPNWc/zaoPwvsFJG/gO+BDqoaZqN8bPXss8+yePFibt68yY0bN1i0aNG9mXlDa9u27b32kiJFioR7zHLlyvHrr79y/vx5goODmTdv3r0Foh7VlStX7k3VPn369Egdwxjz6A6cucaIEd8z+ofP2JMpD5PaD2bkK6W9PieWuzxWE1HVVx5S3jaMsgW4uvzGWaVLl6Zt27aULVsWcLV/hL6UdVfmzJl56qmnaNiwYYTHzJo1K5988gnVqlVDVXn++edp0KBBpOLr06cPbdq0YeTIkVSvXj1SxzDGPJqz1wLoMWY5E+e+z/VESenf5iMmv1mZZIlizzhwmwo+hrl58ybFihVj27Zt99pA4ovY/tkZ86iGL9pKzbdfocCF47RuM5zBg1tRJJvn/99H5VTwNgFjDLJq1SoKFSrEO++8E+8SiDHxTnAwPacOptjZw3Rt0JdOvV6KlgQS1WJPnSkeqFGjBsePH/d2GMYYT1OFrl2RH/8Pxo6lZ9PWXlve9nFZEjHGmOg2ciSMHw+9eiGdO1PI2/E8hjiZRFT1XvdWEzvE5rY5Y9y1bPe/3Jozj0bDekGzZvDpp94O6bHFuTaRJEmScOHCBftRikVUlQsXLpAkSRJvh2KMx2w/fokZn8/m+c/7cujJEgRMmQY+sf8nOM7VRHLkyIG/vz/nzp3zdijmESRJkoQcOWL+6FxjIuP4hZsMGbGYKd99yMlUGenxyiCm+yQkLvzZFOeSSMKECcmTJ4+3wzDGGAAu37xD13ErGD1zACEidG09lDGda5Iuecyc2v1RxbkkYowxMUVAYDCdvt7A+1/1I/P1i7RqOYwPutUnT4bk3g4tysT+C3LGGBMDBYcoXedsodX4gZQ4fYCu9XrRpvvL+MXwWXkflSURY4yJYqrK+4t3UX7cUOoc2MiHz73B091f58Xi2bwdWpSzJGKMMVFs3C+HSP3FCNpt/T8mP92QBN278UblvN4OyyMsiRhjTBS6dOMOF8Z9RZ91M1lUuCp/dxtI/7pxd044a1g3xpgolHbtSj74v9FsLuDH4nc+ZPJLpWLNtO6RYUnEGGOiyqZN0KwZUrIkeX9Yyrj0aUiUIG5f8LEkYowxj0lVkf374YUXIFs2WLqUjJkyejusaBG3U6QxxnjY2asBtBu6mDs1akKCBLB8OWTK5O2woo0lEWOMiaRrAYF0Gv8Lfcf0IPDcBf74cjbky+ftsKKVJRFjjImEgMBgOk35nV4T+pLvgj8dGr/L9SLFvR1WtPNYEhGRqSJyVkR2hyobJCInRWSHc3s+1HP9ReSQiOwXkdqeissYYx5XUHAIXeZs4ZUxAyh3Yjc9X+hOvZ6tqV4os7dDi3aerIlMB+qEUT5KVUs6t6UAIlIYaA4Ucfb5UkR8PRibMcZESkiI0vf7nVQdPYi6B35nSPU3KdzjLV7ye8LboXmFx5KIqq4DLrq5eQPgG1W9rar/AIeAsp6KzRhjIkNVGbp0L/lGD+XVv5YxvnwzEvbsTseq8asdJDS3k4iIpBWRIiKSV0QeJ/l0FpGdzuWutE5ZduBEqG38nbKw4mgvIltEZIutGWKMiU7j1xwiwYjPeXvT98wq9TzHe7xLv7qxeXHbxxduMhCR1CIyQER2AZuAr4D5wDER+U5Eqj3i+SYA+YCSwGlgxN1ThbFtmEsTquokVfVTVb+MGeNHP2xjjPfN2nSMk5+Nof/a6fzwVBV+7zqIoY2LxfuluCMabPg9MBOorKqXQz8hImWAViKSV1WnuHMyVT0Tav/JwI/OQ38g9AXFHMApd45pjDHRIeXiBQxdPp5f8vqxoMtHTG5RmgS+1sE13CSiqjXDeW4rsPVRTiYiWVX1tPOwEXC359YSYK6IjASyAQWAPx7l2MYY4zHLltFgZD9OFy/DhLbDmNauPIkTWN8feIRpT0QkO5Ar9D5O4/nDtp8HVAUyiIg/8AFQVURK4rpUdRR4yznO3yIyH9gDBAGdVDX4UV+MMcZEuQ0boHFjpGhRsq1ZxZzkKeP8fFiPQlTDbHq4fyORT4GXcf3I3/1xV1Wt78HYIuTn56dbtmzxZgjGmDjqxMWbZP1nHwmeqw6ZM8P69XFmOhMR2aqqflFxLHdrIg2Bgqp6OypOaowxMdmhs9fpPfQ7ZkzvRcrUqZCVK+NMAolq7tbJjgAJPRmIMcbEBMcu3KD7iP9j7PS+BAYG8UHnkZAzp7fDirHCrYmIyFhc7Rc3gR0ishq4VxtR1S6eDc8YY6LPqcu36DzyZ0Z/3ZtUATdo1/pTBrz8nLfDitEiupx1t8FhK64eVMYYEyedvRZAxy9WMHxST7Jcu8Drr35Er36vUCZX2oh3jsciSiLPAD8Dq1T1WjTEY4wx0e7ijTt0GLOKj7/sQa7L//LGS4Po0L8VFfKl93ZoMV5EbSJTgRLAUhFZLSJ9RaRENMRljDHR4sqtQN4a/wvvj+9J/gvH6dj4XVr2b0fVgtaQ7o6IBhtuwjXdySARSQ/UAnqKSHFgG7BMVed7PkxjjIl6N24H0WHir/Qe25MiZw7TqWF/GvZ/nTpFs3g7tFjD7cGGqnoBmOfc7k57EtZU78YYEysMnr+FTl/0oszJfXSp15sa/drToGSYc7+ah3AriYhIGqA1kJv7R6xb7yxjTOx05w6DZw0i8fGd9Hq+G2X7duClp+PnmiCPw92ayFJcl7V2ASGeC8cYY6JBYCA0b07Slcu5OnocZSvUo3lZGwsSGe4mkSSq2sOjkRhjTHQIDoY2bWDRIvjiC1J16URzb8cUi7k7Yn2WiLwpIllFJN3dm0cjM8aYKBQYHELXOVs40/RVmDcPhg2DLnZF/nG5WxO5A3wOvMv/FotSIK8ngjLGmKgUGBxCl9lbqDG8P5l3/8I/XfqSp29fb4cVJ7ibRHoA+VX1vCeDMcaYqHa3BlJjeH+a7P6FEZVaEFy3LX28HVgc4e7lrL9xzZ9ljDGxRlBwCN3nbqX65wPuJZBbfQfQu3ZBb4cWZ7hbEwnGNQHjGmwCRmNMLBAUHEL3eVup+ll/mu5ezYhKLbjRZwDvvfBUvF8XPSq5m0QWOzdjjInxgoJD6PHNNp79dABNd69mZKUWXOvdnw9etAQS1dxKIqo6w9OBGGNMVAgOUXp9u43Kw/rTbPcqRlZqwZVe/RhUr7AlEA9wq01ERF4Uke0iclFErorINRG5GsE+U0XkrIjsDlX2uYjsE5GdIrLIGQmPiOQWkVsissO5TXy8l2WMiY9UlR7ztlHRSSCjnnmVSz37Mah+EUsgHuJuw/pooA2QXlVTqWpKVU0VwT7T+e/cWiuBoqpaHDgA9A/13GFVLencOrgZlzHG3COqvD13GC/tciWQC736MaSBJRBPcjeJnAB2q6pGuKVDVdcBFx8oW6GqQc7DTUAOd49njDHhCgqCtm0puPR7/mjThcu9+/Nhg6KWQDzM3Yb1PrjWFPmV+3tnjXyMc78GfBvqcR4R2Q5cBQaq6vqwdhKR9kB7gJy27rExBlxzYbVoAd99B0OHUnbAAJ5WtQQSDdytiQzFNU4kCZAy1C1SRORdIAiY4xSdBnKqailcAxvnikiYl8tUdZKq+qmqX8aMGSMbgjEmDrh5J4jPl+wkuElTVwIZMQIGDACwBBJN3K2JpFPVWlFxQhFpA7wIPHf38piq3sap4ajqVhE5DDzJ/9Z4N8aY+1y/HcRbX23gzdG98P1nK4FjxpLwnc7eDivecbcmskpEHjuJiEgdoC9QX1VvhirPKCK+zv28QAHgyOOezxgTN125Fcib49fQaURXnv1nG73rdmFl1SbeDitecrcm0gnoIyK3gUBAAA2vh5aIzAOqAhlExB/4AFdvrMTASqequcnpifUsMEREgnCNju+gqhfDPLAxJl67fPMOb437hV5jelL61D66v9iDYn3e5vliWb0dWrzk7mDDR27/UNVXwiie8pBtFwALHvUcxpj45fz123Qcu5p3x3ajyJkjdK7fhwp936J1hdzeDi3eCvdylojkjuB5ERHrpmuM8biTl2/x5vCfGTT6HZ46+w8dGw2gyoCOlkC8LKKayOci4gP8AGwFzuHqoZUfqAY8h+sylb8ngzTGxG9Hzl2n68ilDP+6N7ku/0v7Ju9Rv89rNCljf8N6W7hJRFWbiUhhoAWucR1ZcXX13Ytr3fWhqhrg8SiNMfHWgTPX6P3pIsZN60v6m1d486XBtOjfhjpFrQ0kJoiwTURV9+Ba0dAYY6Jd5mOHmDKlB7537vBai094p/+rVC5gY8RiCnd7ZxmDo/62AAAgAElEQVRjTPTbvJnUdesSnDwJnd8eSd8OL1ImVzpvR2VCsSRijImZVq+GBg0gc2Z8V61ifK7c+PjYKPSYxt3BhsYYEy2+3+rPllFT4PnnIW9e2LAB8uSxBBJDRSqJiEhBEZkc1cEYY+K36b/9w6b3R1CqZ3uuPFUM1q6FrNaAHpNFNE6kuIisEJHdIvKRiGQWkQXAamBP9IRojInrVJXhy/dz7INPGL50NL/lKkGHNsPQtGm9HZqJQEQ1kcnAXKAJrjEi23DNaZVfVUd5ODZjTDwQFBxCv+93knDIYD5YPZmlT1ZkfPeRTOxQxWbijQUialhPrKrTnfv7RaQX0E9Vgz0bljEmPrh1J5guc7ZQdcwgWuxYxrfFarKy62CmtypL0kS+3g7PuCGiJJJERErhmnAR4DpQXJw/D1R1myeDM8bEXZdv3qHj5N9oO34AtQ9uYlyFlzjarT8TmhQnoa/1+YktIkoip4HQqxf+G+qxAtU9EZQxJm47dfkWncetov/EfpQ5uZf3a7xFip7d+Lx2QbuEFctENO1JtegKxBgTPwQEBtP50x8YNrkvuS6folODvjzdqz2vVcrj7dBMJNhgQ2NMtEqyfy+zp3Yn6NpVXmv+IS/3bk39Etm8HZaJJEsixpjos24dNGhAsmTJWDDxOzpWKkulAhm8HZV5DJZEjDEed/NOEMl+XAKvvgp58sCyZTTJlcvbYZko4FYXCGfxqZYi8r7zOKeIlPVsaMaY2C4kRPl46V6mvtobbdoUSpd2TWNiCSTOcLcf3ZdABeDukrfXgPER7SQiU0XkrIjsDlWWTkRWishB59+0TrmIyBgROSQiO0Wk9CO+FmNMDBIQGMw7c7aQfvBAOi8YzfbizxC4fAWkT+/t0EwUcjeJlFPVTkAAgKpeAhK5sd90oM4DZf2A1apaANf0Kf2c8rpAAefWHpjgZmzGmBjm4o07tPtyHS981I23/ljIjNIv8FXXzwlOktTboZko5m6bSKCI+OIaG4KIZARCItpJVdeFsU57A6Cqc38GsBbo65TPVFUFNolIGhHJqqqn3YzRGBMDHD1/g25jl/P+1wMoeeoAH1Z/g+AuXfmyXhF8bSbeOMfdJDIGWARkEpGhQFPgvUieM/PdxKCqp0Ukk1OeHTgRajt/p+y+JCIi7XHVVMiZM2ckQzDGeMLWYxf5aPhCxswaSMYbl+nYqD/lerxhY0DiMLeSiKrOEZGtwHO4pkBpqKp7oziWsP5E0TBimQRMAvDz8/vP88YY71i66zRzPpvFtO8/JNAnAa1aDuONns2pUzSLt0MzHuRWEhGRWaraCtgXRtmjOnP3MpWIZAXOOuX+wBOhtssBnIrE8Y0x0UhV+WrdEfYNn8C0pV9wPE0WurcZypBu9SiV06Zyj+vcbVgvEvqB0z5SJpLnXAK0ce63AX4IVd7a6aVVHrhi7SHGxHwhIUrG0Z8x+scRbM3xFD26fsm4d5tYAoknwq2JiEh/YACQVESu8r9LTndwLilFsP88XI3oGUTEH/gAGAbMF5HXgeNAM2fzpcDzwCHgJtDuUV+MMSaa3b6Nb/v2NFk8k5VlajG17bvMaFeetMnd6bxp4gJxdYaKYCORT1S1fzTE80j8/Px0y5Yt3g7DmPjp7Flo1Ah+/x0GD+ZM196kTZ6YRAlsGveYTkS2qqpfVBzL3Yb1/s6gwAJAklDl66IiCGNM7LFm/1lOrdtMi486w7lzMH8+NGtGZm8HZrzC3Yb1N4CuuBq7dwDlgY3YeiLGxBuqytTfjrJ51FRG/d9wbqZNQ7L166FMZJtHTVzgbr2zK/A0cMxZY6QUrjXXjTHxQGBwCAMW7uT8wMFMXDiUgxmeoHX7sdwpUcrboRkvc3ewYYCqBogIIpJYVfeJSEGPRmaMiREu3bhD1+kbaTBhME3+XsOSp55lbvsPmPD6M9b+YdxOIv4ikgZYDKwUkUvYGA5j4rw9p67S78sVDJo2kNKn9jO8ckv8O3ZnetMSJEno6+3wTAzgbsN6I+fuIBFZA6QGlnksKmOM1/248xRTxy5kwrdDSBtwlQ4N+1P0nXaMqpbf1kE390SYRETEB9ipqkUBVPVXj0dljPGa4BBl+Ir9nBr/NXOXjeVi0lS0bjuCt7o1pWZh64Nl7hdhElHVEBH5S0Ryqurx6AjKGOM9Zy9dJ9vgd+m7aSGbnyjKsNeG8PnbNcifKaW3QzMxkLttIlmBv0XkD+DG3UJVre+RqIwx3nHuHFlffplWm9YwrUw9fnurH9NbPk3qpAm9HZmJodxNIoM9GoUxxvu2bXONQD9zBmbMoEjlF2mdK62tAWLC5W7DurWDGBMHhYQoY345SIPda8jTvxtkzAi//QZlylDW28GZWMHdmogxJo65eOMOPeZuofKkz8iz5QfuVKpMogXfQ6ZMEe9sjMOSiDHx0F8nLjPgq9W8O2sIFY/vZGqZ+pzpNYT+lkDMI7IkYkw8oqrM3nycpV9+x5RFn5Am4Do9XuhO1nfa06emTUJhHl1E64nsIowlanGtK6KqWtwjURljotzNO0G8u3AX6SeNY9ba6ZxIk5lOLT6iY5fG1LDxHyaSIqqJvBgtURhjPOrIuev0nLyOt2Z8RJ0DG/n5yYpMaTeQUW9WIWf6ZN4Oz8Ri4SYRVT0WXYEYYzxj2e5/mTx2IaPmf0T2q2f5sPobXOvQidkNi9n8V+axRXQ56xrhX85K5ZGojDFRJvU3s5gzdSCXk6SgVcthNOr8Ei8/ndPbYZk4IqKaSJTPc+BMIf9tqKK8wPtAGuBN/rdOyQBVXRrV5zcm3rh1Czp1osK0aRwoWpaejfrwyVs1KJo9tbcjM3HII/XOEpFM3L887iPPpaWq+4GSzvF8gZPAIqAdMEpVhz/qMY0x/3PjdhDJjxyE5s1h504YOJAnBrzHbBVSJ7PpS0zUcnd53PrACCAbcBbIBewFijzm+Z8DDqvqMZta2pjHExAYzEc//k2a+XPpuWQskiwZLF0KdeuSFEjq7QBNnOTusmQf4lpX/YCq5sH14/9bFJy/OTAv1OPOIrJTRKaKSNqwdhCR9iKyRUS2nDtnK/QaA3Do7HVeGbESv4Fd6PXNpxzJWwT++gvq1vV2aCaOczeJBKrqBcBHRHxUdQ3OJanIEpFEQH3gO6doApDPOe5pXDWf/1DVSarqp6p+GTNmfJwQjIkTFmz1p9+70xj52WvU27ueEZVaMLLPOAIzZ/F2aCYecLdN5LKIpADWAXNE5CwQ9JjnrgtsU9UzAHf/BRCRycCPj3l8Y+K0G7eDeH/xblJP+pK5a6dxPnkaWrUcxvNvv0SPcjlt9UETLSLq4ptYVW8DDYAAoDvQAtfyuEMe89yvEOpSlohkVdXTzsNGwO7HPL4xcdaOE5f5YMpa3pn9MTUO/8nK/OUY36o/H79RjcLZrOe9iT4R1UQ2AqWBiarayimb8bgnFZFkQE3grVDFn4lISVzjUo4+8JwxBtfStRPWHuK3qQuZ+MPnpLt1hUHPtefqmx2Y07AYyRPbdHgmekX0jUskIm2AiiLS+MEnVXVhZE6qqjeB9A+UtXrI5sYYx+I//iHBgP7M2byQY2mz0PKVUbzSoSGNS+fwdmgmnoooiXTAdfkqDVDvgecUiFQSMcZEwp49NO7YAtmxg7klavNTm16MbFORJ9LZ3FfGeyIasb4B2CAiW1R1SjTFZIwJTRXGj4fevZEUKbgw51suZS3NjGfzksDX3Q6WxniGu8vjThGRikDu0Puo6kwPxWVMvPfn0Yv8sGwrQxaPwGf5cteYj6lTSZ8lC528HZwxDndHrM/CNYZjBxDsFCtgScSYKBYQGMyolQf4Z8pcPvl5DCHBt/EZPx46dgTrtmtiGHe7cvgBhVU1rBl9jTFRZKf/ZQbO3Ejz+V/Q/6/l7M6cj/Yv9WfG621IYQnExEDuJpHdQBZcI8mNMVHsTlAI49YcYuv0BYz/6QuyXznLhHJN2dyuK+OalyGFdd01MZS738wMwB4R+QO4fbdQVet7JCpj4pF9/15lwKxN1P9mLHO2/cg/abPSqu3n1H2rKdNs5LmJ4dxNIoM8GYQx8VFwiDJp3RHWT1nIyB9HkfvyaaaWqc8vrbvxcYuy5Eqf3NshGhMhd3tn/erpQIyJb+as3UfCvv2ZvWUJ/qkz0aLlp1Rr35SZz+TBx8dqHyZ2iGjurA2qWimMZXJteVxjHsfGjbTs0BafgweYWeoFfmrRhaGtKpA/U5QvJmqMR0U02LCS8699s415TCEhis/tAPjgAxgxAp8cOTgybzFXMxVmTpV8NnDQxErRvjyuMfFNQGAwY1YfhDVr6L1gBHL4MLz5JgwfTt5Uqejs7QCNeQzeXh7XmDjtj38u8tGs33h1wTia71zB9Ry5SLF6NVSv7u3QjIkS7tZE7i6Pu0pVS4lINVzrgRhjwnAtIJDPlu3n3Iy5fL1yIuluXmFiuSb8/WZ3xlZ/xtvhGRNl3E0igap6QUTuLY8rIp96NDJjYqlVe84wZtZa3l74BXUObGR35nx0enUIjV6vzxdPP+Ht8IyJUt5cHteYOOXk5VsMWbyLtN/MYvbaaSQKDuSTqm050rI9Y5qWJGvqpN4O0Zgo524SaQDcImqXxzUmzvh6/RGWzFnBuz+No5z/3/yeszifNurO6+1q0a94Vht1buIsdwcb3nDuhgAzRMQXaA7M8VRgxsQa169TYMSHLPi/WVxPlIw+dbrg8/prTK/7FGmTJ/J2dMZ4VESDDVMBnYDswBJgpfO4N65p4SOdRETkKHAN19TyQarqJyLpgG9xrVtyFHhJVS9F9hzGeJQqLF4MXbtS5cQJFpWuw9yGHenbqhJ+udN5OzpjokVENZFZwCVgI/AGruSRCGigqjui4PzVVPV8qMf9gNWqOkxE+jmP+0bBeYyJEqrKDztOUS3RdVL37Qk//QTFisG8eTyVrxhzM6YgoQ0aNPFIREkkr6oWAxCRr4HzQE5VveaheBoAVZ37M4C1WBIxMcTuk1f4cMF2/L6fygub5kOiBDBiBLzzDiRMSCFvB2iMF0SURALv3lHVYBH5JwoTiAIrRESBr1R1EpBZVU875zvtjJA3xqsu3rjD58v3c3L+Yj5eOYl8F/1Z+mRFsk+fSIkKxbwdnjFeFVESKSEiV537AiR1HkfFBIzPqOopJ1GsFJF97uwkIu2B9gA5c+Z8jNMbE76g4BDmbD7Ot9+sofvPX1Hz0GaOpsnK6y8PJn/rZlQpXcDbIRrjdRFNwOjrqROr6inn37MisggoC5wRkaxOLSQrrilWHtxvEjAJwM/Pz5brNR6x8fAFPvvuT2r+MIVFfy4m0Dchw6q05dCrr/Nuo5LkzZjC2yEaEyN4Zc1NEUkO+KjqNed+LVzjTpYAbYBhzr8/eCM+E3+dunyLj3/8m4TfzGPir9PJfP0iC4pWZ06DDnRq8Sx9C2WyMR/GhOKthZszA4uc/4wJgLmqukxE/gTmi8jrwHGgmZfiM/GU//Jfea1vN0qf2s+OrAXo2nQglVvXY17lPCRO4LGKuTGxlleSiKoeAUqEUX4BeC76IzLx3smTMHAgZadP53KqdPSu25XAlq0Y9UJhm67EmHB4qyZijFepKmsPnCPxzRtUXDjV1VU3OBh69eJs+640S5CMsnlswKAxEbEkYuKdPaeuMuz/dpFz4Vx6/D4Prl+C5s3h448hTx6e9HaAxsQilkRMvHH6yi1GrdjPpW8X8v6aaeS/6M/mHEU4M2Iq9ds39HZ4xsRKlkRMnHfpxh2+XHuIrQtW0mfV15Q/sZvD6bLTvvFAMrRoRveaBb0dojGxliURE2fduB3ElA3/sGLROt5aPYN3963nQtJUDKzZkdPNWtK3fjGezJzS22EaE6tZEjFx0uYjFxgyfhmtVs1k8a5V3E6QiHEVXmJDo3Z0buhHpQIZvB2iMXGCJRET95w9S/ERH7Loq4koMLP0iyx9sQ1vNK1ApyJZbLCgMVHIkoiJ9YJDlJt3gkgZcMPVVXfUKJLeusWuWo15v3gjXmlamW9KZyeBTdFuTJSzJGJireAQZemu00xa+hcddy/j+Z9nwaVL8PLLMHgw+fPkZ55AkoQ20twYT7EkYmKdu8lj8k9/8cyKb5n+xyLS37rKzZp1SPbZJ1CyJAA2ztwYz7MkYmKNu8nj6x+3U2n5t8z4czFpA66xJm8ZJldtSevuzalTNIu3wzQmXrEkYmK8oOAQftp1mmk/befZZd8w88/FpL59g1X5nubrqi3xa1ab8ZXykDZ5Im+Haky8Y0nExGir9pxh5PzN1F75DTO3LiHV7RusKFCeSVVaUL5pTSZY8jDGqyyJmJjL35+Cnw9l/twZpLhzi5+frMjXVV6lQpMaTLbkYUyMYEnExBjnr98mddKEJDywHz7/HGbPJkdICGtLVuPLsk14pnF1vq6Q25KHMTGIJRHjdccu3GDKhn/Yt3glo48tJ9ua5ZA0Kbz1FtKzJzmTZ2B6qiQkT2xfV2NiGvtfabxCVfnz6CWmrD/M7R+X0nHT9ww5sZurSVOi772HvPMOZMwIQD4vx2qMeThLIiZaBQaHsHTXaWb9spcCK36g55YlPHnhOKdSZmBI9TfZWbcJX75VhUwpk3g7VGOMG6I9iYjIE8BMIAsQAkxS1S9EZBDwJnDO2XSAqi6N7viMZ1y5Gci8P4+z9Oct1P51AZN3LCNtwDV2Z85Hjxe6c6l+E9pVfZL3CmSwua2MiUW8URMJAnqq6jYRSQlsFZGVznOjVHW4F2IyHhQUHELXXpNo+Ot3LNi3AR9VVhQoz6xyDclZvxYdKue1KdmNiaWiPYmo6mngtHP/mojsBbJHdxwmGgQGwqJFJPjiC6b//jvXEiVlRukXWfJsY56rW4Ex5XOSIUVib0dpjHkMXm0TEZHcQClgM/AM0FlEWgNbcNVWLnkvOvOo/r0SwLw/jhN07Di9j/0KX38N//4L+fJx6ePPecOnKC8/V4T5JbLZpIjGxBGiqt45sUgK4FdgqKouFJHMwHlAgQ+BrKr6Whj7tQfaA+TMmbPMsWPHojFq8yBVZdORi8z+/Qg3f/yZV7ctpfrhP/FBkeefh44doU4d8LWkYUxMISJbVdUvSo7ljSQiIgmBH4HlqjoyjOdzAz+qatHwjuPn56dbtmzxSIwmfGevBbBw20mWr9lF2V9/4NUdy8h1+V/OJUvD/OI1SdnlbVq/UtXbYRpjwhCVScQbvbMEmALsDZ1ARCSr014C0AjYHd2xmfAFBYfw64FzfLv5GLeWr6TpXyv55sBvJA4OYtMTRRleuRXna79Ii8r5qVXYZtM1Jj7wRpvIM0ArYJeI7HDKBgCviEhJXJezjgJveSE2E44Rk1eQePZM3t+1mhxXz3IlcXLmlajD4rIvUrR2RTqXz03BLNbLypj4xBu9szYAYQ0EsDEhMcidoBASJfCBmzdh0SKYOpW+v/xCCMKG3CUZVrUtl2rUpXHF/MwrlpWkiazNw5j4yEasm3sCg0NYd+Aci7b5c33tBibzNwnnfwtXr0KePIQMHsxbPkXIW6Yw3f2eIF/GFN4O2RjjZZZE4jlVZceJyyzefpLdqzdTdctK+uxZS84rZwhKkgSaNYPXXoNnn8XHx4fJ3g7YGBOjWBKJpw6eucaPO0/z+7q/KPHbMprt+ZXBZw4TLD78lqsEXzzzKmlaNOO9Vyt4O1RjTAxmSSSe+WnnaaYu+ZN8v62m4Z61dD2+Cx+UHVkLMPi5N/nDrwaVqhTjzVLZKZQllbfDNcbEcJZE4ouzZ2HxYp6eNofamzeQQEM4kjYbXzzzCqtKVqdwFT8alcrOwLzp8fWxCRCNMe6xJBLHqCq7T15l5d4zHNixn/FJj+O7cAH8+iuEhJAhX34mVWjKyqeeIXPVijxfPBsdCmW23lXGmEixJBIH3A4KZtORi6zc8y/71m+nxI711DqwkW7+e/BB4amn4N13oWlTfIoVo8yxS7TOlopkiezjN8Y8HvsViaUu37zD2v3nWL3rJFd+WU/Fvb/T9tAf5L/oD8DejLkZVelVErzUlK7vNLxv36dzp/NGyMaYOMiSSCw0YsGf/DN7IdUObmbwka2ku3WVOz4J2PxEUWaVfoE/i1SkYIXi1CycmSpPZvR2uMaYOMySSAx2/vpt/C/domT2VLB9O6xYAcuX0+233/ENCuRi0lSsyVuG1fnLcaT0MzxTOi91nsrMe7nTksDXx9vhG2PiAUsiMcitO8FsOXaRjYcv8PeWfWTYuI46J/9CT+5Ezp93bVSyJNc6dubNC5kJKleeGsWy061wZgpkSmHLyhpjop0lES8KCAxm2/FLbDp8gb92HyXJpt/xO7qTekd30OfcUQDOJUvD9efrkLLBC1CjBmTJQmpVJt64Q3pbFdAY42WWRLwgKDiEDuNW4/vbb/gd/Yuax3fR7cwRfFBu+yZgS47CfFK1Lb/lLU3yMqUZUK8IJZ5Ic29/EbEEYoyJESyJeNDpK7fYeuwST2ZOyZN6AzZuhA0bSLB2LZO2b8dHldu+CdmWvRBfPPMKm3MW5VqJMpQplI2K+dLTKX8GUiVJ6O2XYYwxD2VJJIrcuhPMntNX2el/mR1HznHtj23k2LeD0qf2kfnCIThz0rVh4sRQsSLrX3mbLzUHl4uVpkyhrFTIm56WedOTMaXVMIwxsYclkUg6duEG6w6cY+eJy5zdtZ/kf++k2KkDlD61j+anD5I06DYA/6ZIx74nS1CuTw+oUAFKl4bEiSl4JYAvgCypk3j3hRhjzGOwJBKBc9duc/rKLYrnSAOqcOQIbN3K9R/XknvdRl48c5i0AdcACPTx5e/MeZlXojbbshdiT64iZC5SgEpPZqRctfz3HdeShzEmLrAk4rhyK5Aj565z4Mw19v17jeP//Evwrt1kOXGIYpf9KZbsCrJ9O1y+DMBTCRPyd7qc/FywIn9nzsfurPm5XagIT+XJRPEcqXkrVzqeyprSxmsYY+K0eJtErt8OYsj//c1p/3OEHDhIxhOHKXjuGE+eP8br546S4+q5e9veSJiEwOJFSdS8uetyVJkyBBUqzJQleymaPTX1s6emf/bUpEgcb99OY0w8Jarq7RjuIyJ1gC8AX+BrVR32sG39/Px0y5Yt9x4HBAZz6vItTl6+xclLzr+Xb3Hm7BU+LJaUvJdOwcGDcOAAevAgZ7fuIvP1i/f2v+OTgMPpc7A/Yy4OZMjF/oy5OJYtHykL5OXTl0ryZOaUHnzlxhgTPURkq6r6RcWxYtSfziLiC4wHagL+wJ8iskRV94S1/bnrt+k0dxtnzl4h6Nhxkp32J/vVs+S4co7sV89S8epZnrh8hmxXz7lms70rY0bkySfZ8VRZ/kqaiX/SZuNoltyQPz95s6WlUOaUlMySkuZZUpEjbVJ8bH0NY4wJU4yqiYhIBWCQqtZ2HvcHUNVPwtq+WOKkuipRsvtqEwAhCP+mTM+pVBnxT52Jo2mzUaLa01SrXwkKFIA0roF7a/adJYGvkCdDcrKltmRhjIkforImEtOSSFOgjqq+4TxuBZRT1c6htmkPtHceFgV2R3ugEcsAnPd2EA+wmNxjMbkvJsZlMbmnoKpGyfX5GHU5CwirKnBfllPVScAkABHZElXZNCrFxLgsJvdYTO6LiXFZTO4RkS0Rb+WemNb/1B94ItTjHMApL8VijDEmAjEtifwJFBCRPCKSCGgOLPFyTMYYYx4iRl3OUtUgEekMLMfVxXeqqv4dzi6ToieyRxYT47KY3GMxuS8mxmUxuSfKYopRDevGGGNil5h2OcsYY0wsYknEGGNMpMXaJCIidURkv4gcEpF+Hj7XVBE5KyK7Q5WlE5GVInLQ+TetUy4iMsaJa6eIlA61Txtn+4Mi0uYxY3pCRNaIyF4R+VtEuno7LhFJIiJ/iMhfTkyDnfI8IrLZOf63TqcJRCSx8/iQ83zuUMfq75TvF5HakY0p1PF8RWS7iPwYg2I6KiK7RGTH3S6XMeB7lUZEvheRfc53q4KXv1MFnffn7u2qiHSLAe9Td+c7vltE5jnffa9+p0SkqxPP3yLSzSnz/PukqrHuhqvR/TCQF0gE/AUU9uD5ngVKA7tDlX0G9HPu9wM+de4/D/yMa8xLeWCzU54OOOL8m9a5n/YxYsoKlHbupwQOAIW9GZdz7BTO/YTAZudc84HmTvlEoKNz/21gonO/OfCtc7+w85kmBvI4n7XvY36GPYC5wI/O45gQ01EgwwNl3v5ezQDecO4nAtJ4O6ZQsfkC/wK5vPw9zw78AyQN9V1q683vFP8beJ0MV4epVUCB6HifHutD9dYNqAAsD/W4P9Dfw+fMzf1JZD+Q1bmfFdjv3P8KeOXB7YBXgK9Cld+3XRTE9wOuOcdiRFzOl3kbUA7XaN0ED352uHrhVXDuJ3C2kwc/z9DbRTKWHMBqoDrwo3MOr8bkHOMo/00iXvv8gFS4fhwlpsT0QBy1gN+8HROuJHIC1w9tAuc7Vdub3yng/9s7/xi5qiqOf76ULl1o07WrFQhtN00kQKG0aPwBptZAjFUEGmu0Cf4AgiElNEaiiVYF9Q+j1aRBLQTbEGoIiy1ljYaQkkawWgUpreu2tIjFNi2VrbVVWmpa6vGPc6d9M8zs7LzdmbdLzye5mXfvvHvud967+8579+4791N4wNpS/pvAV1txnEbrcFbpJJbYk8payTvNbB9A+pxcR1vTNKfH49n4nX+hutKw0RagH3gSv7s6ZGZvVLF/su30/b+BzuHWBCzD/6D+l/KdI0ATeDSGdZI2ycP5QLHnbzqwH3hAPvS3QtI5BWvK8hng4bRdmCYz2wv8ENgN7MP7yCaK7VN9wBxJnZLOxp80ptCC4zRanUjd8CgFUktbUzRLGg88CnzJzP5TtC4zO2Fms/C7//cCFw9gv+maJF0L9JvZpmxxkZoyXGVmVwDzgNslzRlg31boOhMftr3XzAi0CssAAAbmSURBVGYDR/AhkCI1eUM+v3AdsLrers3WlOYVrseHoM4HzsHPYS37TddkZi8A38dv3J7Ah8neGKDKsGkarU5kJIRHeVXSeQDps7+OtmHXLGks7kAeMrO1I0UXgJkdAp7Cx1s7JJVebM3aP9l2+n4i8K9h1nQVcJ2kvwPd+JDWsoI1AWBmr6TPfuAx3OkWef72AHvM7JmUX4M7lZHQp+YBz5vZqylfpKZrgJfNbL+ZHQfWAldScJ8ys5VmdoWZzUn2/0orjtNQxymLSPgd0078TqA0sT6jyW12UT4nspTyCasfpO2PUz5h9Wwqn4SPN78tpZeBSUPQI2AVsKyivDBdwDuAjrTdDmwArsXvHrMTjovS9u2UTzj+Im3PoHzCcSdDnMROdudyamK9UE343euEzPZG4KMjoF9twCO8Atyd9BSqKdnsBm4aIf38fcBWfN5P+D8j3DEC+tTk9DkV2J5+Z9OP05D+KItM+Jjfi/iY+5Imt/UwPvZ5HPfUt+Bjmutxb7++dKDTSflp0vUX4D0ZOzcDL6V00xA1fRB/zOwFtqT0sSJ1ATOBzUlTH/CtVD4deDbZXw2clcrHpfxL6fvpGVtLktYdwLxhOo9zOeVECtWU2v9zSltLfXgE9KtZwHPpHPakC0nRms4GDgATM2VFa/o2fqHuA36OO4Ki+9QGYFvqU1e36jhF2JMgCIIgN6N1TiQIgiAYAYQTCYIgCHITTiQIgiDITTiRIAiCIDfhRIIgCILchBMJgiAIchNOJAgKQNIlkr4gD+k/oQXt3SDpZ5J+KekjOW1Ml7RS0prh1heMXsKJBC1H0nxJJumiButtzNne4Tz1MvU7JC0aio0qjMXfcp4PVNUnqV3S05LGDLUxM+sxs1vxkOWfzmljp5ndUqGxTdJvM+E+gtOMcCJBESwEfoeHgCgjLZZzRrW8mV3ZDDGVbVahA18TYjiZAjyAvxVc60nkZmCtmZ1oxLCkyyT9uiKVord+A39TOW/9MszsGP4mdC7HFLwFGMpr9pEiNZqA8cBe4EJgeyrrAl4AluNhUz5UkZ+W9juMRypdlLF3N3Bn2u7BQ3JvBb6Y2edwFR2VbU4boH43cBQPLbM0ld2Ih7DYgq+5UDXmEe58+oBdwB0V351b51htBLoy+dXAT3AHvAsPfbMKD/+zso4tpWN3TUX559Nv7gU2DPIcrqnIXw48XnTfilRMKlxApNMrpYvvyrS9EY8S24Wv9fH+VF6Wz9Q9jK+b8nSmbBswNW2X4gK1pwt3Z6leFR1vamOA+l2UB9+8GPgVMDbllwOfq9LGJ/HovGPxBX9OLlo0iOPUBvyjomw78OW0/R1OLSTUBhwkxWqqYW9xchb3Abelsgnp+LWlfEcdTZ2p/t8oX0xpDLC/6L4VqZgU45hBq1mIh2IHv8NfiA+v7DKzP2b2q8wDYGabJU2WdD4eNfigme1OXy+WND9tT8GXBz0wgJbKNgZb/2rg3cCfJIE7nf4q+y0GbjUPF75P0nEGP4T8duBQKSNpHD6sVjp2R3FnvC99/zpwrJYxM7sHuKei+ETS/iNJD5rZcwMJMrMDwG1Vyk9IOiZpgpm9VveXBW8pwokELUNSJ76ex6WSDL+DNfxO/kjF7pX5LGuABcC5uCNC0lx8nYcPmNnrkp7Co6cOxMk2Gqwv4EEz+1otw2mtl5lm9mLKnwf803wOYTAcrWh/Br6eRml1xsuBe5PtC4BXzKyhaKrpd14KfAK4X9IKM1veiI0MZwH/zVk3GMXExHrQShYAq8xsmpl1mdkUfL2CCxq0041Pyi/AHQr4Qj8H04XxInyNhEYYqP5rlE9+rwcWlCabJU2SNK3C3iXAxPRvsWcA3wN+XPpS0iOS7qwlxswOAmPSEwjAZXiI7xIz8XkMcIfSS4NIepeZHTGzbnyd8HpOt5adTnw463ie+sHoJpxI0EoW4nMEWR4Fvt6IETPbil/U95aGc/AlQc+U1At8F3jTUFgdatZPwzi/l9QnaamZbcP/y2ld2v9JfG4iy2zgIXwtml5gt5ndDyDpevyiPbOOpnX45Dm4E9mS6o8D2pOjgXKH0ghLJO2Q9Dy+KFLep5APA4/nrBuMcmI9kSBoApKWAX8ws0cqyscBK8zsRkk9ZnbDADZm4xPpn22y3CEhaS0+0b6jaC1B64knkSBoDrNITw4VfAUYL+k+YIak9loGzGwz8JvheNmwWUhqA3rCgZy+xJNIELQISVOBuyy99S3pLuAJM3umWGVBkJ9wIkEQBEFuYjgrCIIgyE04kSAIgiA34USCIAiC3IQTCYIgCHITTiQIgiDITTiRIAiCIDfhRIIgCILchBMJgiAIchNOJAiCIMjN/wEAXlx7Nb0gsgAAAABJRU5ErkJggg==\n",
            "text/plain": "<matplotlib.figure.Figure at 0x7f1ca6d8bf28>"
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "The correlation coefficient between the polynomial and real curve is almost equal to one:"
    },
    {
      "metadata": {
        "trusted": false
      },
      "cell_type": "code",
      "source": "corr=np.corrcoef(fw[0:nr,1],y)\n(corr[0,1])**2",
      "execution_count": 6,
      "outputs": [
        {
          "data": {
            "text/plain": "0.9999998803855714"
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "We now store the coefficients of the polynomial in a simple text file for later use, together with the constants used.\n"
    },
    {
      "metadata": {
        "trusted": false
      },
      "cell_type": "code",
      "source": "poly = open('polyconstants.txt','w') \n \npoly.write('alpha: '+str(alpha)+'\\n')\npoly.write('beta: '+str(beta)+'\\n')\npoly.write('coeffs: '+str(fit))\n\npoly.close() \n\nfile = open('polyconstants.txt','r') \nprint(file.read())",
      "execution_count": 7,
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "alpha: 3.25\nbeta: 0.762\ncoeffs: [ 3.48030475e-11  1.75845637e-06  3.70090706e-03 -4.20322042e-08]\n"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "## Data processing"
    },
    {
      "metadata": {
        "trusted": false
      },
      "cell_type": "code",
      "source": "# Note: If your data file is already stored locally, comment out lines below\n# Copy remote file to local file\nfile_name='http://iahs-icwrs.org/ival0998.txt'\nresponse = requests.get(file_name)     # Get remote file\nfile_content = response.text           # Copy file content into variable\n\nlocal_file='temp.txt'      # Temporary file to stire data locally\n\nfh = open(\"temp.txt\", \"w\")\nfh.write(file_content)\nfh.close\n# Note: If your data file is already stored locally, comment out lines above\n\n# Open local file and split into time stamps and interval times\nread_file = open('temp.txt','r')\ncontent= read_file.readlines()\n\nfirst_t = True            # Set flag for first timestamp\nold_millis = 0            # Zero milliseconds at start\n\nfor line in content:\n  if line[0] == 't':       # timestamp\n    first_t = False       # Reset flag first timestamp\n    \n    ftime = line.split(',')     # file time\n    millis = ftime[7].strip('\\n')    # Extract milliseconds without line return \n    \n    rdate=datetime(int(ftime[1]), int(ftime[2]), int(ftime[3]), \n                   int(ftime[4]), int(ftime[5]), int(ftime[6]))     # Date & time from timestamp\n    \n    print(rdate)\n    \n  elif not first_t:        # Consider times only after first time stamp\n    print('ms')          # time\n\n    ",
      "execution_count": 8,
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "2018-05-18 08:38:55\n2018-05-18 09:08:55\n2018-05-18 09:38:55\n2018-05-18 10:08:55\n2018-05-18 10:38:55\n2018-05-18 11:08:55\n2018-05-18 11:38:55\n2018-05-18 12:08:55\n2018-05-18 12:38:55\n2018-05-18 13:08:55\n2018-05-18 13:38:55\n2018-05-18 14:08:55\n2018-05-18 14:38:55\n2018-05-18 15:08:55\n2018-05-18 15:38:55\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\n2018-05-18 16:08:55\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\n"
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "ms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\n2018-05-18 16:38:55\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\nms\n2018-05-18 17:08:55\n2018-05-18 17:38:55\n2018-05-18 18:08:55\n2018-05-18 18:38:55\n2018-05-18 19:08:55\n2018-05-18 19:38:55\n2018-05-18 20:08:55\n2018-05-18 20:38:55\n2018-05-18 21:08:55\n2018-05-18 21:38:55\n2018-05-18 22:08:55\n2018-05-18 22:38:55\n2018-05-18 23:08:55\n2018-05-18 23:38:55\n2018-05-19 00:08:55\n2018-05-19 00:38:55\n2018-05-19 01:08:55\n2018-05-19 01:38:55\n2018-05-19 02:08:55\n2018-05-19 02:38:55\n2018-05-19 03:08:55\n2018-05-19 03:38:55\n2018-05-19 04:08:55\n2018-05-19 04:38:55\n2018-05-19 05:08:55\n2018-05-19 05:38:55\n2018-05-19 06:08:55\n2018-05-19 06:38:55\n2018-05-19 07:08:55\n"
        }
      ]
    },
    {
      "metadata": {
        "trusted": false
      },
      "cell_type": "code",
      "source": "# 1 Check span during which logger was active\n\n# 2 Extract millisecond intervals\n",
      "execution_count": 9,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": false
      },
      "cell_type": "code",
      "source": "",
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "python36",
      "display_name": "Python 3.6",
      "language": "python"
    },
    "language_info": {
      "mimetype": "text/x-python",
      "nbconvert_exporter": "python",
      "name": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.6",
      "file_extension": ".py",
      "codemirror_mode": {
        "version": 3,
        "name": "ipython"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}