notebooks/kr08/2_1/fortran_original.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Program 2.1 (SIR model) - original Fortran code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting Program_2_1.f\n"
     ]
    }
   ],
   "source": [
    "%%writefile Program_2_1.f\n",
    "C\n",
    "C This is the FORTRAN version of program 2.1 from page 19 of\n",
    "C \"Modeling Infectious Disease in humans and animals\"\n",
    "C by Keeling & Rohani.\n",
    "C \n",
    "C It is the simple SIR epidemic without births or deaths.\n",
    "C\n",
    "C This code is written to be simple, transparent and readily compiled.\n",
    "C Far more elegant and efficient code can be written.\n",
    "C\n",
    "C This code can be compiled using the intel fortran compiler:\n",
    "C ifort -Vaxlib  -o Program_2_1  Program_2_1.f  \n",
    "C \n",
    "\n",
    "\n",
    "C Main program starts here.\n",
    "\n",
    "      program main\n",
    "\n",
    "      REAL beta\n",
    "      REAL gamma\n",
    "\n",
    "      REAL S,I,R\n",
    "\n",
    "      REAL S0\n",
    "      REAL I0\n",
    "      REAL MaxTime\n",
    "      REAL EVERY, step, t\n",
    "\n",
    "      INTEGER GivesName\n",
    "\n",
    "      CHARACTER*2000 str, FileName\n",
    "\n",
    "      COMMON /parameters/ beta, gamma\n",
    "      COMMON /variables/ S, I, R\n",
    "\n",
    "      GivesName=iargc()\n",
    "\n",
    "      if (GivesName.eq.0) then\n",
    "         beta=1.4247\n",
    "         gamma=0.14286\n",
    "         S0=1 - 1.0d-6\n",
    "         I0=1.0d-6\n",
    "         MaxTime=70\n",
    "      else\n",
    " \n",
    "c     \n",
    "c     READ IN ALL THE VARIABLES\n",
    "c\n",
    "\n",
    "         call getarg(1,FileName)\n",
    "\n",
    "         open(1,file=FileName,STATUS='OLD',ACCESS='SEQUENTIAL')\n",
    "  \n",
    "         \n",
    "         read(1,*) str\n",
    "         read(1,*) beta\n",
    "         \n",
    "         read(1,*) str\n",
    "         read(1,*) gamma\n",
    "         \n",
    "         read(1,*) str\n",
    "         read(1,*) S0\n",
    "         \n",
    "         read(1,*) str\n",
    "         read(1,*) I0\n",
    "         \n",
    "         read(1,*) str\n",
    "         read(1,*) MaxTime \n",
    "         \n",
    "         close(1)\n",
    "      endif\n",
    "C\n",
    "C     Check all variables are OK & set up intitial conditions */\n",
    "C     \n",
    "      \n",
    "      if ( S0.le.0) then\n",
    "         write(*,*) \"ERROR: Initial level of susceptibles (\",S0,\") is \n",
    "     .        less than or equal to zero.\"\n",
    "         STOP\n",
    "      endif\n",
    "\n",
    "      if ( I0.le.0) then\n",
    "         write(*,*) \"ERROR: Initial level of infecteds (\",I0,\") is \n",
    "     .        less than or equal to zero.\"\n",
    "         STOP\n",
    "      endif\n",
    "\n",
    "      if ( beta.le.0) then\n",
    "         write(*,*) \"ERROR: Transmission rate beta (\",beta,\") is \n",
    "     .        less than or equal to zero.\"\n",
    "         STOP\n",
    "      endif\n",
    "\n",
    "      if ( gamma.le.0) then\n",
    "         write(*,*) \"ERROR: Recovery rate gamma (\",gamma,\") is \n",
    "     .        less than or equal to zero.\"\n",
    "         STOP\n",
    "      endif\n",
    "\n",
    "      if ( MaxTime.le.0) then\n",
    "         write(*,*) \"ERROR: Maximum run time (\",MaxTime,\") is \n",
    "     .        less than or equal to zero.\"\n",
    "         STOP\n",
    "      endif\n",
    "\n",
    "      if (S0+I0.ge.1) then\n",
    "        write(*,*) \"WARNING: Initial level of susceptibles+infecteds\n",
    "     .        (\",S0,\"+\",I0,\"=\",S0+I0,\") is greater than one.\"\n",
    "      endif\n",
    "\n",
    "      if (beta.lt.gamma) then\n",
    "        write(*,*) \"WARNING: Basic reproductive ratio (R_0=\",\n",
    "     .        beta/gamma,\") is less than one.\"\n",
    "      endif\n",
    "\n",
    "\n",
    "      S=S0\n",
    "      I=I0\n",
    "      R=1-S0-I0\n",
    "     \n",
    "C     \n",
    "C     Find a suitable time-scale for outputs \n",
    "C     \n",
    "      step=0.01/((beta+gamma)*S0)\n",
    "\n",
    "      Every=1.0/((beta+gamma)*S0)\n",
    "      if (Every.gt.1) then\n",
    "         Every=10.0**INT(log10(Every))\n",
    "      else\n",
    "         Every=10.0**INT(log10(Every)-1)\n",
    "      endif\n",
    "      DO WHILE (MaxTime/Every.gt.10000)\n",
    "         Every=Every*10.0  \n",
    "      ENDDO\n",
    "\n",
    "      open(1,recl=3000,file='Program_2_1_f.out',ACCESS='SEQUENTIAL')\n",
    "C     for F77 use\n",
    "C      open(1,file='Output_Risk',ACCESS='SEQUENTIAL')\n",
    "C\n",
    "\n",
    "C\n",
    "C     The main iteration routine\n",
    "C     \n",
    "\n",
    "      t=0\n",
    "      write(1,*) t,S,I,R\n",
    "      DO WHILE (t.lt.MaxTime)\n",
    "         \n",
    "\n",
    "         CALL Runge_Kutta(step)\n",
    "         \n",
    "         t=t+step\n",
    "         \n",
    "C     If time has moved on sufficiently, output the current data\n",
    "         \n",
    "         if( INT(t/Every).gt.INT((t-step)/Every) ) then\n",
    "            write(1,*) t,S,I,R\n",
    "         endif\n",
    "      \n",
    "      ENDDO\n",
    "      \n",
    "      END\n",
    "  \n",
    "\n",
    "\n",
    "      SUBROUTINE Runge_Kutta(step)\n",
    "\n",
    "      REAL InitialPop(3), tmpPop(3)\n",
    "      REAL dPop1(3), dPop2(3), dPop3(3), dPop4(3)\n",
    "\n",
    "      REAL S,I,R, step\n",
    "      COMMON /variables/ S, I, R\n",
    "\n",
    "C\n",
    "C     Integrates the equations one step, using Runge-Kutta 4 \n",
    "C     Note: we work with arrays rather than variables to make the\n",
    "C     coding easier\n",
    "C\n",
    "      InitialPop(1)=S\n",
    "      InitialPop(2)=I\n",
    "      InitialPop(3)=R\n",
    "\n",
    "      CALL Diff(InitialPop,dPop1)\n",
    "      do k=1,3\n",
    "         tmpPop(k)=InitialPop(k)+step*dPop1(k)/2.0\n",
    "      ENDDO \n",
    "\n",
    "      CALL Diff(tmpPop,dPop2)\n",
    "      do k=1,3\n",
    "         tmpPop(k)=InitialPop(k)+step*dPop2(k)/2.0\n",
    "      ENDDO\n",
    "\n",
    "      CALL Diff(tmpPop,dPop3)\n",
    "      do k=1,3\n",
    "         tmpPop(k)=InitialPop(k)+step*dPop3(k)\n",
    "      ENDDO\n",
    "\n",
    "      CALL Diff(tmpPop,dPop4)\n",
    "\n",
    "      do k=1,3\n",
    "         tmpPop(k)=InitialPop(k)+step*(dPop1(k)/6 + dPop2(k)/3 +\n",
    "     .        dPop3(k)/3 + dPop4(k)/6)\n",
    "      ENDDO\n",
    "\n",
    "      S=tmpPop(1)\n",
    "      I=tmpPop(2)\n",
    "      R=tmpPop(3)\n",
    "\n",
    "      RETURN\n",
    "      END\n",
    "\n",
    "\n",
    " \n",
    "C The Main Differential Equations.\n",
    "\n",
    "      SUBROUTINE Diff(Pop, dPop)\n",
    "\n",
    "      REAL Pop(3), dPop(3)\n",
    "\n",
    "      REAL beta\n",
    "      REAL gamma\n",
    "\n",
    "      COMMON /parameters/ beta, gamma\n",
    "\n",
    "\n",
    "C     Set up temporary variables to make the equations look neater\n",
    "\n",
    "      REAL tmpS, tmpI, tmpR\n",
    "\n",
    "      tmpS=Pop(1)\n",
    "      tmpI=Pop(2)\n",
    "      tmpR=Pop(3)\n",
    "\n",
    "C\n",
    "C     The differential equations\n",
    "C\n",
    "\n",
    "C     dS/dt =\n",
    "      dPop(1) = - beta*tmpS*tmpI\n",
    "\n",
    "C     dI/dt =\n",
    "      dPop(2) = beta*tmpS*tmpI - gamma*tmpI\n",
    "\n",
    "C     dR/dt =\n",
    "      dPop(3) = gamma*tmpI\n",
    "      \n",
    "      RETURN\n",
    "      END\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gfortran -O3 -o Program_2_1_f Program_2_1.f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " WARNING: Initial level of susceptibles+infecteds            (  0.999998987     +   9.99999997E-07 =   1.00000000     ) is greater than one.\r\n"
     ]
    }
   ],
   "source": [
    "!./Program_2_1_f parameters.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use(\"ggplot\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sir_out = pd.read_csv(\"Program_2_1_f.out\",sep=\"\\s+\",engine=\"python\",header=None,names=[\"t\",\"S\",\"I\",\"R\"],index_col=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAENCAYAAABUyo/yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd8VFX++P/XnZJMyqQOEBJKKIEQmqALimAvCCjWu4LY1o9tZdWvfmy/ddVF113LrmJZldXPghWuii4Ciu6Ci6IgolKkSSchgfTeZub+/rgzIYS0SSZTyPv5eNzH3HLunfeEIe+cc889R9F1HSGEEEJ0nCnYAQghhBDhTpKpEEII0UmSTIUQQohOkmQqhBBCdJIkUyGEEKKTJJkKIYQQnSTJVAghhOgkSaZCCCFEJ0kyFUIIITrJEuwAOkmGbxJCiI5Rgh3AiSTckymHDh3q0HkOh4OCggI/R9N1winecIoVwivecIoVwivecIoVOhdvamqqn6MR0swrhBBCdJIkUyGEEKKTJJkKIYQQnSTJVAghhOgkSaZCCCFEJwWkN6+qqv8HTAOOaJo2opnjCjAXmAJUATdomvZDIGITQgghOitQNdP5wORWjl8EZHiWW4BXAhCTEEII4RcBqZlqmrZaVdX0VopMB97UNE0H1qqqmqCqam9N03K7Ip64J57AsnkzyU4nKMpxi24yNbu/zWNWK7rdjttuR7fbcTkcuAYMwDlwIG6HwygnRDei6+B0Qn29Qm2t8VpXZ7x6151OBacT3G7j1eU6ft3loqFMa+vGcnTdG4PbraDr3nUa1o1Fwe0Gm81MZWUcumcoGKOc0sI5R6/Z3DHvdb3Xavyq68f/jJoeO/qqtFjm3HNN3HCDf/+9RMeFyqANacDBRtvZnn3HJVNVVW/BqL2iaRoOh8PnN7Ps2YPp66+J7FisHaL374/7vPNwX3cd+vjxPidWi8XSoc8aDOEUK4RXvP6K1eWC8nIoK4PycoXSUmO9okKhuhqqqozFWD+6z7vtXfcudXUK9fVQV3f8ouvhNEBAbLADaDeHQw+b7213ECrJtN00TZsHzPNs6h0ZAcRy//0k3ncfpcXFjf/MNMbWau7PV7cbAKWlP229x2trMVVUoJSXYyovx3T4MJa9e7Hs2oVp/37Mb7yB+Y03qD39dEqffBLn4MHtjjmcRmcJp1ghvOJtGmttLRQVmY5ZiouP3S4pMVFebqK8XGl4rawMXN9Ds1nHatWJiACrVcdqhYiIo/ssFh2LBUwmY91kArPZ93WzWfe8GuuKYlzzaGOS3tCQ1LiByVjXiY2NoaqqspnjRlXw+HOOHmt6TuNjXt6/n73Hm+5r+djx11EUGDo0TkZACiGhkkxzgL6Ntvt49nUJZ2YmusNBXaB+gbpcWLdswbZsGTHvvEPkmjU4LryQ4ldeofaCCwITgwgrug6FhSYOHjRz+LCZvDwTeXlmSkrM7N+fRF6esb+0tONJ0W53ExurY7e7sduN15gYnaioo4vNdux2c/siIrwLDa9G0tTp3dtBcXG4/KESRUFBZbDDaDeHQydM/gbsFkIlmS4BZququhAYD5R21f3SoDCbqR89mvrRo6m47TbiH3mE6I8+Iummmyh+/XVqLrww2BGKIHC5YP9+M7t3WzhwwML+/WYOHjRz4ICFAwfMVFW1lCjNDWsWi05SkpukJDeJie6G9cZLQoIbu91NXJxObKzxGhNj1Oi6mtncdhkhTgSBejTmPeAswKGqajbwKGAF0DTtVWA5xmMxuzAejbkxEHEFg56URMmLL+JKS8P+0kskzJ5N/ooVuAYODHZooou43bBnj5kdO6zs3Gnhl18s7NxpZc8eC7W1Ld87T0hw06ePk5QUNykpLlJSXAweHE1MTCm9erno3dtIoNKvTYjgC1Rv3hltHNeBOwIRS0hQFMoffBDLgQNELVlCwn33UfjBB9Lb9wTgTZybN0ewcaOVzZutbNlipaKi+Wpg794uBg92kp7upH9/J337uujf30Xfvk4SEo6fYdDhsFFQUNvVH0MI4aNQaebtfhSFkj/9iYg1a4hcuxbbihXUTG7tUVwRiurrYfNmK999F8HatZGsXx9BScnxiTMlxUVWVj1DhjgZMqSejAwnGRlO7HaZkleIE4Ek0yDSk5KouOsu4h95BPuzzxr3TqV2GtJ0HXbtsrBqVSRffhnJd99FUF19bPJMSXExalQdo0bVNyw9eriDFLEQIhAkmQZZ5TXXEPvii1i3bSPiu++oGz8+2CGJJmpr4auvIvniCxtffhlJdvax/20GDapn/Pg6xo2r49RT6+jTxyV/EwnRzUgyDTabjaoZM7C/8ALRb78tyTRE1NTA6tWRfPJJFF98YaO8/GjtMznZxZln1nL22bVMmlQrtU4hhCTTUFA1cyb2F17A9tlnxnAyUVHBDqlb0nVYvz6ChQujWbbMdkynoayseqZMqebcc2sZMaI+II+VCCHChyTTEODq25e60aOJ2LgR2+rV8txpgOXlwWuvxbJwYTR79hz9LzFiRB3TptUwdWo1Awe6ghihECLUSTINETVTphjJdPlySaYBsnmzlX/8I4YlS6zU10cA0KuXi6uuquKqq6oYPFgSqBCifSSZhoia884j7s9/JvLrr432RunB0iV0Hb74IpJXX41l3TpjqgOTSefCC6uZMaOKs8+uxSL/K4QQPpJfGyHCOXQoLocDc14e5t27cfkwCL5omzeJ/vWvdrZsMWqhdrubGTOquOeeCOz24iBHKIQIZ5JMQ4WiUDdhAlFLlhC5Zg1Vkkz9Qtdh5cpInn76aBLt1cvFbbdVMHNmFbGxumcmliAHKoQIa9InMYTUnnYaABHr1wc5khPD9u0WrrkmieuuS2bLlgh69XIxZ04pa9Yc5pZbKomNldGHhBD+ITXTEFI3diwAET/9FORIwltRkYlnnrHz9tvRuN0KcXFu7rqrnOuvr5SnjoQQXUKSaQhxDh2KbrNh2bsXpbgYPTEx2CGFFV2HxYujePTROIqLzZjNOjfcUMm995aTlCQDKwghuo4084YSq5X6rCwAIjZtCnIw4eXgQTOzZiVx552JFBebmTixln//O58//alUEqkQostJMg0xdSedBIBVkmm76DosXBjFOef04MsvbSQkuPnb34pZuLCQIUOcwQ5PCNFNSDNviHFmZgJg2bkzyJGEvpIShQcfTOCTT4wbodOmVfPEE6UyVq4QIuAkmYaY+qFDAbDu2BHkSELb+vVW7rgjkZwcCzExbv70p1KuvLJaxroQQgSFJNMQ4/QkU8uuXeBygdkc5IhCz9tvR/Pww/HU1yuMGVPHSy8Vk54uQ/8JIYJH7pmGGN1ux5mailJbi3nfvmCHE1Lq6uCBB+J54IEE6usVbrqpgo8+KpBEKoQIOqmZhiDn0KFYDh3CunMnrkGDgh1OSCgqUrjppiS++y6SyEidp54q4aqrqoMdlhBCAFIzDUnOgQMBpGbqcfCgmUsvdfDdd5H07u3io48KJJEKIUKK1ExDkHPAAAAse/cGOZLg27LFwrXXJnPkiJlhw+p5++1CUlKkt64QIrRIzTQEudLTAbB085rp+vVWrrzSwZEjZiZMqGXx4gJJpEKIkCTJNAQ5PcnU3I1rpuvXR3DNNcmUl5uYOrWat98uJC5OBqYXQoQmSaYhyNWnD7rZjOXQIajufvcG162LYObMJCorTVx2WRV//3sxkZHBjkoIIVomyTQUWa24+vYFwHLwYJCDCawNG6zMmpVEVZWJyy+vYu7cEixyZ18IEeIkmYYoV58+AJhzcoIcSeDs3GnhuuuSGxLp88+XyJgVQoiwIMk0RLlSUwEwHzoU5EgCIyfHzIwZyZSUmDj//Bqee04SqRAifEgyDVHdKZkWFSnMmJFEXp6ZceNqeeWVImnaFUKEFUmmIcqVlgac+M289fVw661J7N5tZdiweubPLyIqKthRCSGEbwL297+qqpOBuYAZeF3TtL80Od4PWAAkeMo8qGna8kDFF2q6S8300Ufj+eabSHr2dPHmm4XEx8vjL0KI8BOQmqmqqmbgZeAiIAuYoapqVpNiDwOapmljgKuBvwcitlDVkExP4Jrp/PnRLFgQQ2SkzhtvFJGaKgMyCCHCU6CaeccBuzRN26NpWh2wEJjepIwOxHnW44ETu0rWhoZkmpsL+olXW1u3LoJHHokH4JlnShg7tj7IEQkhRMcFKpmmAY0fmMz27GvsMWCWqqrZwHLgd4EJLTTpsbG44+NRamsxFRYGOxy/Kiw08dvfJuJyKdx+ewVXXNH9BqYQQpxYQqnP5AxgvqZpf1VV9TTgLVVVR2iadkzbn6qqtwC3AGiahsPh6NCbWSyWDp8bMH37QmkpyVVVmMMhXo/WfrZuN9xwg4W8PBMTJrh59tkILJbgfq6w+C54hFOsEF7xhlOsEH7xnugClUxzgL6Ntvt49jV2EzAZQNO0b1VVtQEO4EjjQpqmzQPmeTb1goKCDgXkcDjo6LmBktSrF7YtWyjfupXYsWNDPl6v1n62L74YyxdfxJGY6GLu3HxKSoJ/nzQcvgte4RQrhFe84RQrdC7eVM9tJOE/gUqm64EMVVUHYCTRq4GZTcocAM4F5quqOgywAfkBii8knWidkNavt/L003YAXnihRDocCSFOGAG5Z6ppmhOYDawAthm7tJ9VVZ2jquolnmL3AjerqroReA+4QdO0E6/njQ9cvXsDYDp8OMiRdF5VlcLddyfidhv3Sc85pzbYIQkhhN8E7J6p55nR5U32PdJofStweqDiCQfuHj0AMOeHfwX9ySft7NtnYdiweu67ryzY4QghhF/JCEghzOVJpqYwT6Zffx3BP/8Zi8Wi8/zzMp2aEOLEI8k0hLl79gTAfORIGyVDV3m5wj33JABw993ljBjhDHJEQgjhf5JMQ9iJUDN9+mk7OTkWRo2qY/bsimCHI4QQXUKSaQhze54hMxUUgMsV5Gh8t2mTlfnzYzCbdZ59tgSrNdgRCSFE15BkGsoiInAlJqK43RBGz7+BkfsffDAet1vhppsqGT5cmneFECcuSaYhznvfVAmz+6ZvvRXNxo0R9O7t4t57y4MdjhBCdClJpiHO+3gMeXnBDcQHhw/DU08ZcxbMmVNKbGy3flxYCNENSDINcS5vzTSMBm549FEzZWUmzjmnhosuqgl2OEII0eUkmYa4cKuZbt1qYf58ExaLzmOPlaIowY5ICCG6niTTEOd9PCZcaqZPPBGHritcd10lgwaFXw9kIYToCEmmIa6hZhoGyXTVqkj++18b8fE6/+//yTOlQojuQ5JpiHOHyT1TpxMef9zodPTAAy6SkmRGGCFE9yHJNMS5vJP/hvgoSB98EMWOHVb69HFyxx2SSIUQ3Ysk0xDnTkoCQCksDHIkLaurg+eeM+YpfeCBcmy2IAckhBABJsk0xLkTE42VggLQQ/N5zUWLosnOtpCRUc/06dXBDkcIIQJOkmmoi4zEHRuL4nSilIfeSEK1tTB3rlErveeecszmIAckhBBBIMk0DHibek1FRUGO5HjvvRdNbq6ZzMx6pk2TARqEEN2TJNMw0JBMQ+y+aW0tvPiiUSu9995yTPJtEkJ0U/LrLwyEas30gw+iycszM2xYPZMnS61UCNF9STINA95OSO1JptuKtvH5/s+pqq/q0phcLnjllVgAZs+ukFqpEKJbswQ7ANE2d3IyAKbi4lbLvfzTyzy5/kkA+tn78d6U90iPS++SmD77zMbevRb69XMybZr04BVCdG9SnwgD7Wnm3VK4hSfXP4lJMZEWm8aB8gNc99l1VDv9n+h0HV5+2aiV3nprBRb5k0wI0c1JMg0D7Ummr2x8BYAbs25k5RUrGZIwhN2lu5n741y/x7NmTQQbN0aQnOzi17/u2uZkIYQIB1KnCANtJdOC6gI+2fMJZsXMraNuJTYilmfOeIbpS6bzj83/4DfDf0PP6J5+i+e114xa6W9+U0lUlN8uK4QIPxOAx4H4YAcSAKXAH4BvmjsoNdMw0FYyXXVwFS7dxaS0SaTFpgFwSq9TuLD/hdS4api3eZ7fYtm928zKlTZsNp3rrqv023WFEGGpuyRSMD7n4y0dlGQaBtpKpisPrgTg3L7nHrP/dyf9DoD3drznt3un8+fHAHD55VUkJYXm8IZCiIDpLonUq8XPK8k0DLSWTN26m9U5qwE4p985xxwb03MMox2jKaktYcnuJZ2Oo7xcYdGiaABuvFFqpUII4SXJNAy44+PRFQVTSYnxgGcje0v3UlJbQkp0SrOPwVyfdT0AC7Yu6HQcmhZNZaWJ006rJSvL2enrCSFOPGaz+eTMzMysjIyM4RdddNHA8vJyv+aZF154Ifm6667r11qZpUuX2r/44osY7/bTTz/d46WXXkr2ZxxNSTINBxYLJCai6LqRUBvZWLARgNE9Rjd76iWDLsFutbOxYCO/FP/S4RDcbvjnP43v5m9+I7VSIUTzIiMj3du3b9/6yy+//Gy1WvW//vWvPQIdw8qVK+1fffVVrHf7/vvvz589e3aXjscqyTRceCYJb9rU+1P+T0DLyTTKEsXUAVMBWLxrcYfffvXqSPbutZCa6uSCC2ToQCFE2yZOnFixa9euSIDHHnusV0ZGxvCMjIzhc+bM6QmwY8eOiAEDBgy/5JJLBgwcOHD45MmTG2qyaWlpI3Nzcy0Aq1evjh43btzQptd/991340eNGpU5bNiwrAkTJgw5ePCgZceOHRFvvvlmj1dffbVXZmZm1meffRZ7zz33pD7yyCO9AL755puo0aNHZw4ZMiTr/PPPH5Sfn28GGDdu3NDbb789beTIkcPS09NHfPbZZ7FN36817UqmqqpaVFX9j6qqLfZkEl1LbyGZbsrfBMBJPU5q8dzLMy4H4KNdH+HW3R16/3feMe6VzppVJYM0CCHaVF9fz4oVK+JGjhxZ/dVXX0W/++67yRs2bNj2/fffb3vzzTd7rFmzJgpg3759ttmzZx/Zs2fPz3a73f3MM8+0uyZ7/vnnV/z000/bt23btvXKK68smjNnTsrQoUPrrrvuuvzbbrvt8Pbt27dOnjy5ovE5N9xww4Ann3wye+fOnVuHDx9e/cADD6R6jzmdTmXz5s3bnnrqqYNz5sxJPf4dW9auZKppmhPIBFJ8uXhjqqpOVlV1h6qqu1RVfbCFMqqqqltVVf1ZVdV3O/peJyTvkIKNkqmu6+wo3gFAVlJWi6ee1vs0esf05mDFQb4//L3Pb52fb+Lzz22YzboM0iCEaFVtba0pMzMza+TIkVl9+vSpu+uuuwq+/PLL2ClTppTExcW54+Pj3VOnTi1etWqVHSAlJaXuggsuqAS49tprC7/55pt21wj37t0bMWnSpIwhQ4ZkvfDCCynbt29v9cn3wsJCc3l5uXnq1KkVADfffHPh2rVrG97vqquuKgaYMGFCZXZ2doQvn9uXZt45wKWqqp6pqqrVlzdRVdUMvAxcBGQBM1RVzWpSJgN4CDhd07ThwN2+vMeJrrmaaX51PmV1ZSREJuCIcrR4rkkxcflgo3bakabe99+PxulUOPfcGlJSOlazFUJ0D957ptu3b9+6YMGCgzabrdVn6BRFaXbbbDbrbrfx+6a6urrZXDV79ux+v/3tb4/s3Llz60svvbS/tra2U7cuvbFaLBZcLpfSVvnGfHnjV4AkYCVQo6qqy7O0p1vnOGCXpml7NE2rAxYC05uUuRl4WdO0YgBN0474ENuJr5ma6S8lRoeiQfGDjvtCNnXJoEsA+HTfp7jcrlbLNqbr8O67RhPvzJlSKxVC+O7ss8+uWL58eUJ5ebmprKzMtHz58sSzzz67HCA3Nzfi3//+dwzAO++8kzRhwoQKgD59+tStWbMmGkDTtMTmrlteXm7u169fPcD8+fMbeuva7XZXeXm5uWn55ORkV1xcnMt7P/SNN95IPu200yqalusIX+9+Nfcbuz3ZOw042Gg7GxjfpMwQAFVV1wBm4DFN0z5reiFVVW8BbgHQNA2Ho+UaWWssFkuHzw0GpYdxGyGmuhqbJ+7c/bkAjEwZ2eZnOTP5TAYmDGRPyR521OzgjH5ntOt9//tfxdPxSOeqq+xYLPY2zwm3n204xRtOsUJ4xRtOsUJ4xTtx4sSqmTNnFo4dO3YYwLXXXpt/+umnV+/YsSMiPT295sUXX+x5yy23RGdkZNT87//+bz7AI488cui2225LnzNnjmvChAnlzV3397///aEZM2YMio+Pd06cOLH8wIEDkQBXXHFFyZVXXjno008/TXj++ecPND7nn//8597bb7+9/5133mnq169f7XvvvbfPH5/Rl2Q6wB9v2AoLkAGcBfQBVquqOlLTtGOeBdE0bR7gHR9PLygo6NCbORwOOnpuMPRMSMAE1B46RIkn7p+yjZ68faP6tuuzXNT/Il4ueZl3f3yXrOiW77E29uqrCYCVq66qoKSk2e/zccLtZxtO8YZTrBBe8YZTrNC5eFNTfepb45Oqqqofm9v/2GOPHX7ssccON91vsVj417/+tbfp/smTJ1fs27dvS9P9d955ZyFQCDBr1qySWbNmlTQtM2rUqNqdO3dubXwt7/qECROqN27cuL3pOd99990O73rv3r2dOTk5m5v9gC1odzOvpmn7NU3bD0QCw73bnn1tyQH6Ntru49nXWDawRNO0ek3T9gI7MZKrAHTvKEiN5jTdU7oHgIHxA9t1De8jMsv3LW9Xr96yMoXly437+TNmSBOvEEK0pN01U1VVkwANOBvQPR2GfgGe1DTtkTZOXw9kqKo6ACOJXg3MbFLmY2AG8E9VVR0Yzb572hvfCc/bAalRMt1fbvwd09/ev12XGOUYRZ/YPmRXZLPh8AZ+lfKrVssvX26jpkbhtNNq6du3/fdZhRCiPYYOHVr3yy+//BzsOPzBlw5IzwDnAHWA4qk9rgUubutEz6M1s4EVwDZjl/azqqpzVFW9xFNsBVCoqupWYBVwn6ZpXTpiRThpWjN1uV3kVBiV+772vi2e15iiKA2106V7l7ZZ/oMPjI5HV1zh/wnGhRDiROLLPdPJwGfADuBOz76twK/bc7KmacuB5U32PdJoXQfu8SyiKU9vXsWTTPOq8qh31+OIchBtjW73ZaYOmMprm19j2d5lPHrqo5iU5v+eyskx8+23kdhsOlOnSjIVQojW+FIzjQKKm+xzYNRURVdLNHqGm0pLweXiYLnRObq9tVKvMT3H0DumN7mVuQ1DETZn8WLjXun559cQFydTrQkhRGt8SaabgGl4HmlRVfVZjCbejV0Ql2jKYsEdH4/idqOUlnKg3Ojt3d77pV4mxcSUAVMAWLqn+aZeXYcPPzSS6RVXSMcjIYRoiy/J9GGMnrynYjxbeg/gBh7zf1iiOW5v7bS4uMM1U4BpA6YBsHzvcnT9+Frnli1WfvnFSlKSi7POqu1ExEIIEXwPPPBAyuDBg4cPGTIkKzMzM2vlypUxbZ/lm3bfM9U07WtVVU8Bbgf6A/uAeZqmbfJ3UKJ57sRE2LcPU3ExOZVG56O02DSfr3NKr1PoGdWTgxUH2VywmVE9Rh1z/KOPjFrp9OnVWH0aOFIIIULLv//975gVK1YkbN68eWtUVJSem5trqa2t9WmowPbwaRxDTdO2APcB/w+4XxJpYDWumeZV5gHQO6a3z9cxKSYuGnARAMv2LTvmmK7DJ5/YACOZCiFEOMvJybEmJSU5o6KidDAGZEhPT6/39/v48pxpIsbIQ5c32rcYuFXTtKIWTxR+405IAIxkmqsYQwl2JJmC0at3wdYFLN2zlAdPebBhbN8ffrBy6JCFlBQXJ5/s9++bEKI7U5STu+S6ur6hpUOXXnpp2Z///OfU9PT0ERMnTiybMWNGkXfWGH/ypWb6BnAFxv1S73I58Lq/gxLNO6ZmWtXxminA+JTxJNuS2Ve2j21F2xr2L11qNPFOm1aNSaaOF0KEufj4ePeWLVu2vvTSS/t79OjhvP766we98MILyW2f6RtfnjM9D9gLXAZsB4YBH3n2iwDwJtOq4sOUJZcRaY4kMbLZyRTaZDFZmJw+mXe2v8OyvcvISs5C12HpUqOJd9o0aeIVQvhZKzXIrmSxWJg2bVr5tGnTykeNGlX91ltvJXvG+PUbX+oeu4FVmqZt0jStTtO0jRgjFe3yZ0CiZd5kmldxCICU6JQ2p15rTcNYvXuNsTS8Tby9e0sTrxDixLBx48bIzZs3R3q3f/zxx6g+ffr4fXyEVmumqqo2nqfrTeBRVVV/4mjN9Ark0ZiAaUim1cZUrx1t4vWakDqBhMgEdpbsZGfxTpYuPQWAqVOliVcIcWIoKysz33nnnf3KysrMZrNZT09Pr12wYEF7JmjxSVvNvF8CTR9EnNtoXQGeBZ73Y0yiBd5kmltvtE6kxKR06npWk5UL+1/Iop2LWLpnGUuXTgTg4ouliVcIcWKYNGlS1Y8//njclGv+1lYyPcDxyVQEidsz2P0hvRTofDIFo6l30c5FfPjf/Q1NvGPHShOvEEL4otVkqmlaeoDiEO2ge2umpkqg8828ABPTJmK32tn3nTFww4UX1kgTrxBC+MiX3rwAqKqaDBwzFJOmaQf8FpFokbeZ91BEDWB0QOqsSHMk5/c/n8U7pgNGMhVCCOEbXwZtOB/jWdOm49fpvlxHdJweFYUeGUl2jDFerj9qpgDjIq5m8ZGRmGzlnHqqjMUrhBC+8qVB71WgD8cO2qD4eA3RGYqCOzGRXLux6Y97pgBlm88GwD14KXk10sgghBC+8qVGmQysAK7SNM3vQzGJ9nEmJnA4xhj9KNnmn0E8Vv3bk52HLmH5vgxuG3WbX64rhBDdhS+1yheBfkCaqqp+H3FftE9hj1icZog3xWCz2Dp9vaIihXXrIjCZ3TD40xbnOBVCiHAVHR09pqvfw5ea6YfAbGArgKqq3v26pmlyzzRA8hzG2Lk9lFi/XO8//7HhditMmFjNT/Z6fsz/kf1l++kf59uk40II0Z35UjN9B4hH7pkGVV5SBAA93VF+ud7nnxu12ymT65icPhmAj3d/7JdrCyFEd+FLjbIfsB64HyjpmnBEW/LizAD0rI9so2Tb6uvhq6+M65y7Is4lAAAgAElEQVR3Xi399EtZvGsxH+36iDtPurNT4/4KIURTyh+7Zgo2/dHgDKDfmC/JdB4wDvhW0zS/DxIs2uewp3W3V42509f64YcIystNDB5cT9++LlLcZ5BkS+KXkl/4uehnRiSP6PR7CCFEd+BLMj0XGAHkq6q6D3B59uuapnXNhK/iOEdsTqiDXn7oT71ypVErPess49lSq8nKxQMvZsHWBXy862NJpkIIvwqFGmRX8eV+5yhPeTswEjip0SIC5IjVaBRIKXN3+lpffmkk03POOTpQw2WDLgOM+6ZuvfPvIYQQ3YEvNdMbuywK0W5HzMaMLr2KOtfSfuSIiS1bIrDZ3IwffzSZntzrZPrE9iG7Ipt1ees4rfdpnXofIYToDtqdTDVNW9CVgYj2OaIb7bspBZ2bJs1bK50woQ5bo8dVTYqJSwdfyks/vcRHuz6SZCqECHtVVVU/dvV7+DI27/+1cEjXNO0mP8Uj2pDvLAag9+HKTl1n1Sojg5599vFj8V426DJe+ukllu5ZypzT5vhlcAghhDiR+dLMewPGoPbe5yW86zogyTQAdF2noNZIpim5ZRS7XGD2vVevywWrV3s7Hx0/S0xmUiYjHSPZXLCZFftXMH3Q9M4FLoQQJzhfOiD9EZjjeX0CWIKRSFuqsQo/K6ktod5dT1wtRDnBVFraoev89JOVkhIT6elOBg50NVvm6iFXA7Box6IOxyuEEN2FL/dM/9h0n6qq8zh+SjbRRQqqCwDoVWMBnChFRZCU5PN1vE283kdimjN90HT+uPaPrM5ZTU5FDmmx8s8shBAt8eWeab8mu+KAocDodp4/GZgLmIHXNU37SwvlrgA+AH6ladr37Y2vO8ivzgegV30k4MRUXEzz9crWeTsfNdfE65VoS2Ry+mSW7FnC+zvf5+6xd3fgnYQQonvwpZl3b5NlIzAR2NfWiaqqmoGXgYuALGCGqqpZzZSzA3cB63yIq9vwJtMeejQApuJin69RVqawcaMVi0VnwoTWH6+5eqjR1Kvt1OSZUyGEaIUvybTpAPfVwDcYHZPaMg7YpWnaHs9QhAuB5nq1PA48BbRcZerGvMm0p8mYf7QjyXTt2gjcboUxY+qIidFbLTsxdSK9Y3qzv3w/6/Lk7xshRHgym80nZ2ZmZmVkZAw/55xzBhcUFHR+PNYmfLln2pnZYdKAg422s4HxjQuoqjoW6Ktp2jJVVe/rxHudsBpqppZEoGPJ9OuvjSbeiRPbHvTBbDKjDlGZ++NcFu5YKM+cCiHCUmRkpHv79u1bAS6//PL0Z555psdTTz2V58/3aDOZqqra1m25Ts9nqqqqCfgb7ajlqqp6C3ALgKZpOByODr2nxWLp8LnBYLFYqHAbAzakxvcGILa2ligfP8PatcY/1ZQpNhyOtmeeuXX8rcz9cS7L9i7j5Wkvk2BLaFes4fazDZd4wylWCK94wylWCL94Q8Wpp55auWnTJv/MYdlIe5KgP+bhygH6Ntru49nnZccYRP9Lz6TjKcASVVUvadoJSdO0eRgz2ADoBQUFHQrI4XDQ0XODweFwcLDYqNzHW4yEVpOTQ6kPnyE/38TPP6dgs7kZNCif9pwaTzxnpJ3B6pzVvPLtK9w88uZ2xRpuP9twiTecYoXwijecYoXOxZuamurnaNpHUeiaKdh02jWAvtPpZNWqVfabbrrJ7//Q7UmmY5psxwKzARUj0f7UjmusBzJUVR2AkUSvBmZ6D2qaVgo0/ImlquqXwP9Kb95jeR+NccQZj6n42sy7Zo1REx0/vo5IH6ZDvT7relbnrGbB1gXcNOImTIrMBy+ECB+1tbWmzMzMrMOHD1sHDRpUc+mll5b5+z3aTKaapm0EUFXVBtwB3Af0BDYDf9Q0bXE7ruFUVXU2sALj0Zj/0zTtZ1VV5wDfa5q2pBOfods4Un0EAEdiH6AjyTQCgNNP922Q/PP6nUdqTCp7y/bydc7XnNHnDJ/OF0IIaH8N0t+890zLy8tNZ511VsZf/vKXng8//PARf75He+6ZRgK3A/cDvYBtwO80TXvflzfSNG05sLzJvkdaKHuWL9fuDnRdp7C6EIBkRzrgezI92vmo5cEammMxWZg1bBZPf/8087fOl2QqhAhLdrvd/cILLxy46qqrBj/wwANHrFar367dnva6PcBfMWqjGvAwUKuq6iXexW/RiBaV1JRQ564j1hpLpCMF8C2ZHjhg5sABC/HxbkaMqPf5/WcOnYnVZOWLA1+QU5HT9glCCBGCTj/99OrMzMzqefPm+T58XCvac8+0N0cHtVc9S2N6O68jOuFw1WEAHFEO3ImNHo3RdVDa7iPmvV962mm1HRkbnx7RPZg6YCof7/6Yt7a9xYO/etD3iwghRBA0nYJt5cqVu/z9Hu2pmR5oYznY8qnCX45UGs37PaJ6QFQUbpsNpa4OpaqqXed//bVxv9TXJt7Gbsi6AYB3t79LjVPG1RBCCK/2dEBKD0Acog2HK4yaaY/oHgDoiYmQm2uMzxsT0+q5ug7ffmvUTH3tfNTYKb1OYUTyCLYUbuHDXR9yTeY1Hb6WEEKcSOQZhzBxuNKTTKOMZHpMU28bDhwwc/iwmcREFxkZzg7HoCgKt426DYDXNr0m4/UKIYSHJNMwcaSqUTMvviXTdeuMJt5x4+rac3u1VdMGTiMtNo3dpbv5Yv8XnbuYECLcdWxS5fDV4ueVZBomGpp5myRTpR3JdP36o8m0s6wmKzePMEZBemXTK52+nhAirP2B7pNQSzE+b7OkF26Y8Dbz9ozuCfhWM/3uO/8lU4CZmTN57ofnWH94Pd8f/p5Tep3il+sKIcLON8C5wQ4iFEjNNEwcd880ORkAcxtjcxYWmti1y4rN1rHnS5sTY43h2qxrAXh106t+uaYQQoQzSaZhomnN1NXDSKqm/PxWz/M28Y4dW09EhP/iuWn4TUSYIvhs32fsLN7pvwsLIUQYkmQaBty6uyGZJtuMGqnbM/WSqY2aqb+beL16Rvfk6qFXo6Pz/I/P+/XaQggRbiSZhoGS2hKcbifxEfHYLDYA3J6aqbmNmmlXJVOA2SfNJsIUwZLdS6R2KoTo1iSZhoH8KiNhegdsAHC1o2ZaVaWwebMVk0nn5JP9n0zTYtOYkTkDHZ3nfnjO79cXQohwIck0DHinXvN2PoKjNdPW7pn+8IMVp1Nh+PB6YmP1LontjtF3EGGK4JM9n7CjaEeXvIcQQoQ6SaZhIL/aUzNtlEx1ux09MhJTVVWL4/P68/nSljSuncq9UyFEdyXJNAw018yLohxt6m2hdtqV90sbmz16dkPtdEvBli59LyGECEWSTMOAt2baM6rnMftba+p1OmHDhsAk09TYVK7Puh4dncfXPY6ud02TshBChCpJpmGgYVzexjVTjj4e09zADVu3WqmsNJGe7qRnz64fkP6uMXcRHxHP14e+ZsWeFV3+fkIIEUokmYaBlmqmrQ3cEKgmXq9EWyJ3jrkTgIdWPoTT3fHZaYQQItxIMg0DzXVAgtYHbmg8U0yg3Dj8RvrG9mVrwVa0nVrA3lcIIYJNkmkYaEimTZt5Wxi4Qdcb9+StDUCEhkhzJA+NewiAZ75/hoq6ioC9txBCBJMk0xDndDsprC5EQWkYStCrpd68e/eayc83k5zsYuBAV8BiBbhk4CWMSx3Hkeoj/O2HvwX0vYUQIlgkmYa4wppCdHR6RPfAYjp2xryWevM2fr60s5OB+0pRFOZeMBcFhde3vM62om2BDUAIIYJAkmmIa+h8FNPzuGOu3r0BMOflHbM/0J2PmhrbeyzXZ12PS3fx0NcP4da7vjexEEIEkyTTEJdXaSTK3rG9jzvmSkkBPMnUfTRhrVsXCQQvmQLcf8r9OKIcrD+8nvd/eT9ocQghRCBIMg1xuZW5APSJ63P8wago3AkJKPX1mAoLAcjPN7F3r4XoaP9NBt4R8ZHxPDL+EQCeWPcERTVFQYtFCCG6miTTEOetmabGpjZ7vKGpN9dIuo0nA7dYmj0lYC4ffDkTek+gqKaIP3zzh+AGI4QQXUiSaYjz1kzT4tKaPe5NpibPfdNgPF/aEkVRePaMZ4m2RPPx7o9ZtndZsEMSQoguIck0xLW7ZnroEBCc50tb0z+uP78f/3sAHvr6IQqrC4MckRBC+J8k0xDX6j1Tjm3mraxU2LLFitmsM3Zs8O6XNnXdsOs4PfV0CmsKeXDNgzIQvhDihCPJNMTlVbWzZpqXx4YNVlwuhREj6omJCZ2EZVJM/PWMvxJjjWH53uVov8hQg0KIE0vAuqioqjoZmAuYgdc1TftLk+P3AP8DOIF84Deapu0PVHyhqLK+krK6MmxmG0lRSRRWHd9E6m5UM12/PviPxLSkr70vj5/2OPesvoffr/k9Y3uMJSMxI9hhCSGEXwSkZqqqqhl4GbgIyAJmqKqa1aTYj8ApmqaNAj4Ang5EbKHM28SbEpOC0sJQRo2beUOp81Fz1CEqlw++nGpnNbevvJ1qZ3WwQxJCCL8IVDPvOGCXpml7NE2rAxYC0xsX0DRtlaZpVZ7NtUDzNwm7EW8y7R1z/IANXt5k6jqUzw8/WIHQTaaKovDn0//MwPiBbCvaxmPfPhbskIQQwi8C1cybBhxstJ0NjG+l/E3Ap80dUFX1FuAWAE3TcHgGe/eVxWLp8LmBUplbCUD/xP4tx+twoNvtbCofSjUmMjJ0MjOTAhzpsVr72TpwsPCKhUxaMIm3t7/N2YPPZtbIWQGO8Fjh8F3wCqdYIbziDadYIfziPdEF+bH+46mqOgs4BTizueOaps0D5nk29YJm5vJsD4fDQUfPDZSdeTsBSLIk4XQ6W4y3R9++fLV1EgAnn1xFQUFpwGJsTls/2zRzGo9PeJz7v7qf25ffTg9TD07udXIAIzxWOHwXvMIpVgiveMMpVuhcvKmpzXdoFB0XqGbeHKBvo+0+nn3HUFX1POD3wCWapoXGg5JBlF2RDbTck9fL2b8/X2Ek01Bt4m3qmsxruD7reurcdfzPF//T0KQthBDhKFDJdD2QoarqAFVVI4CrgSWNC6iqOgZ4DSORHglQXCHtQPkBAPrZ+7VaztnvaDIdPz48kinAH0/7I6f1Po0j1Ue46fObpEOSECJsBSSZaprmBGYDK4Btxi7tZ1VV56iqeomn2DNALPC+qqo/qaq6pIXLdRsHy43bzP3j+rda7ueosRTioHdUEf37B3Yy8M6wmqzMO28e/ez92Fiwkd+u/C1OtzPYYQkhhM8Cds9U07TlwPIm+x5ptH5eoGIJBy63i+xyo5m3T2zrHZu/LjsJgInRG1CU4V0emz8l2ZJYcOECLvvkMj7f/zkPff0QT096usVHgYQQIhTJCEghKrcyF6fupFd0L2wWW6tlvzmYDsCZ9f8JQGT+NyRxCPMvnI/NbOPdHe/y7IZngx2SEEL4RJJpiNpfbgz+1Nb9Ul2HbzcnAnBm2VKoqeny2LrCr3r9ilfOfQWzYub5H5/n1U2vBjskIYRoN0mmIcp7v7SvvW+r5Q4cMJOXZyHZVEQWW7FkZwcivC5xQf8LeOaMZwB4fN3jklCFEGFDkmmI2l9m1Ezb6ny0dq0xhOCExJ8xoWPet6+rQ+tSvx7ya56ZJAlVCBFeJJmGqPbWTL/7zkimpw80aqSWX37p2sACYGbmzGMS6nM/PCfTtgkhQpok0xDlvWfa395WzdSYKWb8eGOMC+u2bV0bWIB4E6qCwrMbnuXhbx7G5Q6fx36EEN2LJNMQpOs6e0r3AJAel95iubw8E/v2WYiNdTP0bGM8Xuv27YEIMSBmZs7k1XNfJcIUwfyt8/ntyt9S6+r2A2MJIUKQJNMQVFBdQEltCXarnV7RvVos9803Rq30V7+qg+FDAbDs2gXOE2fgg2kDp/HORe9gt9pZuncpVy69ksNVh4MdlhBCHEOSaQj6pcS47zk4cXCrgxd89ZWRTCdOrEW323H27YtSW4tl796AxBkoE1In8OHFH5Iak8oPR35gysdT+Cn/p2CHJYQQDSSZhiBvMs1IyGixjK7D6tVGMj3jDKPp05mZCYDlBLlv2tjw5OF8etmnjOs1jrzKPC7/5HLe3/l+sMMSQghAkmlI2lWyC2g9me7aZSEvz4zD4SIz02jWrfck0xOlE1JTjigHi6Yu4prMa6h11XL3f+/mri/voqKuItihCSG6OUmmIWh7kdGJqLVk6q2VTppUi8nzr1g/bBgA1i1bujbAIIowR/D0pKd5dtKz2Mw2PvjlAyZ/NJlN+ZuCHZoQohuTZBpidF3n58KfARjhGNFiucbJ1Kt+7FgArD/+aLQDn8BmZM7gs8s+Y1jSMPaW7eWSJZfw3A/PUecKnynohBAnDkmmISa7IpvSulKSbcmkRKc0W6auDr791hiswXu/FMDVpw+unj0xFxdj3rMnIPEGU0ZiBkunL+XGrBupd9fz7IZnmfLxFDbmbwx2aEKIbkaSaYjZUmA00Y50jGyxJ+/atZFUVprIzKynd2/30QOKQt3JJwMQsWFDl8caCmwWG0+c/gSLpiyiv70/24q2Me1f03j020cprS0NdnhCiG5CkmmI2VJoJNMRyS038a5YYUzJdsEFx88QU3fKKQBErl3bBdGFrolpE/nPlf/h1pG3AvD6lteZqE3krW1vychJQoguJ8k0xHx/+HsARvcY3exxXYfPPzful1544fHJtHbiRAAiV68+4e+bNhVlieKRUx/h00s/ZXzKeIpqinjw6weZ/NFkVmevlvF9hRBdRpJpCKl31/PDkR8AY37P5vz8s4VDhyz06uVi1Kj64447s7JwJSdjzs3Fsnt3l8YbqkY4RvDhtA955ZxXSItNY2vRVmZ8OoPLP7mcNYfWBDs8IcQJSJJpCNlauJUqZxUD4gbQI7pHs2WWLIkCjCZeU3P/eiYTtZMmARC5cmVXhRryFEXhkkGX8N+r/stDv3qIhMgEvjv8HeoylSuXXsnqHKmpCiH8R5JpCPk291sAfpXSfK3U7YaPPjKS6WWXVbd4nZoLLgDAtny5nyMMP1GWKGafNJu1V6/lvpPvIz4inm9zv2XG8hmcv/h8Fu1cRK1TBs8XQnSOJNMQsurgKgDOTDuz2ePr10dw6JCFtDSnMbh9C2rPPRfdZiNy/XpMubldEmu4sUfYuXvs3aydsZYHTnmAnlE92Va0jXv+ew8Zf8/gqfVPNUzILoQQvpJkGiIq6ipYl7cOk2LijD5nNFvm/feNWumll1Y338TrocfGUnPOOQBEL1rk91jDWVxEHHeOuZO1M9by/JnPMyxpGIcrD/PCTy8wYdEE1GUqH+/6mBrn8Z27hBCiJZJMQ8Sq7FXUu+sZ02MMSbak444XFh5t4lXVqjavV3nttQDEvPXWCTUlm79EmiO5ashVfHH5F3xxzRdcMfgKbGYbaw6t4Y5Vd3DS2yfxu1W/4/P9n8scqkKINlmCHYAwfLTrIwAuHnhxs8f/+U8TNTUmzj67hsGD235usm7SJOoHDcK6eze2zz6jZto0v8Z7olAUhTP6nUFWdBaPT3icj3d/zKIdi9hYsJHFuxazeNdi7FY7F/S/gMnpk5mUNgl7hD3YYQshQozUTENAUU0RKw+uxKSYmD5o+nHHa2rglVfMANx4Y2X7LqooVN54IwD2v/1NaqftEB8Zz/VZ17P8suV8rX7Ng796kOHJwymvL+fDXR9y879vZsSbI7hy6ZW8/NPLbC3cKj2ChRCA1ExDwptb36TeXc+5fc+lZ3TP44+/GUN2tsKwYfWcfXb7mxyrZswg9rXXsO7YQfR771HlafoVbRsQP4DfnfQ7fnfS79hTuodle5ex8sBKNhzZwLe53/Jt7rc8uf5JkmxJjE8Zz/iU8ZzW+zSGJQ3DbDIHO3whRIApYf6XtX7o0KEOnehwOCgoKPBzOL6rrK9kwqIJFFQXsGjKIiamTTzmeGmpwumn96S42MyCBYWcd55v9+9sn3xC0m234Y6PJ/+zz3D16+fP8JsVKj/b9vIl3tLaUlbnrGbVwVX8N/u/5FXlHXPcbrUzpucYRjlGMarHKEY7RpMWm9biOMtdGWsoCKd4wylW6Fy8qampAP75UgpAaqZB9/LGlymoLmBMjzGcnnr6cccffjie4mIzkya5Ofdc3zvC1EybRvUFFxD1+eck3norhR9+iB4d7Y/Qu6X4yHguHngxFw+8GF3X2V++n3W56/g291vW5a3jQPkBVuesZnXO6oZzkmxJjHKMYmjiUIYmDiUjMYMhCUOIjYgN4icRQviTJNMg2pS/ib9v/DsAj5766HG1l/ffj2Lx4miiotz8/e9OOlS5URRKnnsO60UXEbFpE0mzZlG0YAG6XTrRdJaiKKTHpZMel86vh/4agEMVh9hUsImN+RsbXotqivgy+0u+zP7ymPPTYtMYkjCEgQkD6W/vT/+4/vS396ePvQ9RlqggfCIhREdJMg2SnIoc/uff/0O9u54bsm44btSjzz+P5N57EwB47LEyhgyJpqMtUHpCAoXvvIPjqquIXLeOHpMnUzx3LvWeGWaE/6TGppIam8rk9MmAMdl7dkU2mws2s7N4JztLdrKzeCe7S3aTU5FDTkUOq7JXHXedlOgU+tn70cfeh17RvUiJSSElOoXMmkyinFH0jO6J1WQN9McTQrQgYMlUVdXJwFzADLyuadpfmhyPBN4ETgYKgV9rmrYvUPEF0vq89dz6n1s5XHWYk3uezB/G/6HhWH09vPxyLH/9qx23W+GOO8qZNasK6FzTrGvgQAoWLybp5pux/vwzPaZPp3rqVCpvvJG6U0+lY9Ve0RZFUehr70tfe1+mDJjSsN/pdrKvbB87i3eyv2w/+8v3s79sPwfKD5Bdnk1eVR55VXl8d/i75q+LQnJUMj2iepAYmUiSLYnkqGSSbEkkRRrriTZjf0JEAvYIO7HWWOkcJUQXCUgyVVXVDLwMnA9kA+tVVV2iadrWRsVuAoo1TRusqurVwFPArwMRXyC43C6+P/w9C7Yt4F+7/wXAqSmn8sYFb2Cz2CgsNLF0qY1//COWvXuNf5Z77innnnvK/RdD//7kL1mC/bnniP3HP4hatoyoZctw9e5NzZlnUjduHM5hw6jPyIAoaWbsShaThcEJgxmcMPi4Y063k9zKXPaX7edQ5SHyKo3EerjyMAV1BWSXZnOk+ggF1QUUVPvWXBFjjcEeYSfOGoc9wt6wxEXEEWuNJcoS1bBEW6KP2fYuNovt6LrZhtVsxaJY/NbJSohwFKia6Thgl6ZpewBUVV0ITAcaJ9PpwGOe9Q+Al1RVVTRN83t345yCCo5UQGFREW5dx+32LOjoOke3dR0dcOtu3O6jx3Q8ZXQ3bjcN57lcbqqd1VTUVVFRX0lZbTm5FXlkl+ewrWg7ldVOqI3DUjeTSY5pjCy+iP/v0wh27rSybdvRJrsBA5w88UQpZ53VBSPv2GyUP/QQlddfT8w77xC9cCHm3FxiFi4kZuFCAHRFwe1w4OrVC3evXriTk3HHxqLHxKDHxuKOiTE6MVmt6BbLca+Kw4G1ogIsFqPG61l07y9bk+mY/ceVaW3xRXvLV1VhKi7ummt3oHwE0F+x0t88GOIGQ9zRY0lJSRQVFeF0O8mvLaKwrpiiuhKK6ksprC2myLNdWFdCUV0xhbUllNWXU+asoMJZSWW9seSR1+L7d4SCQoTJSoQpAqvJ0rBus9ow6ybPthWrZ3+EyYJVsWI1WTArFsyKybMY6xbFjEkxYzGZMWPCbPKWMWNWzFgU8zHlvfu964qiYMKEoigogEkxoaBgUhTg6DETSkNZu91OZUXl0XMbnWdcRznmOo23vcePrhs/E6Dhj4yG7UadaJs7dnQfxx5rVFZBwV2Xjiki0a//jqLjApVM04CDjbazgfEtldE0zamqaimQDPi9r/pFs4op3DwESPX3pdvFCazyLF42m8748bWoajVTplQTEdG1MbhTUym/7z7K770Xy9atRH71FRE//YRlxw4se/Zgzs/HnJ8PW7Z06PrNTyAXulKCHYAPvLH28fE8twLlEVBqg7JIKI08dr0sEqqtUGmFqmaWyogm21aotUCdGVwmnVp3HbXulidgEP6llvbluf9dG+wwhEfYdUBSVfUW4BYATdNwOBw+XyPavpeiGCNHK4oOig54X5vuA8XzCrpRwfAca27dpJgwmxTMJhNmk4kIs5UISwTR1ijs0Vbi4iAuDuLjdVJSYPBgncGDdcaM0bHZTECMZzmWxWLp0Gdtl5494ayzMD4h1DudcPgwSm4uHDqEUlgI5eVQWYlSXm6sV1VBfT1Kfb1xo9e7OJ0oTid6XR24XKDrzS9uN+g6SivHGhY4+tpevpTvymuHUHkFo5Ib5y1f51n8wKXo1JmgzqxTa4Y6k06dybNuNo7VmnXqPMdqPeWcJh2XAk7FSMhOBWPbZFzTZfIcU3TPvkbnmI7f7/ZcS1eM77Jb8bQkHbOv0aui48Y45vb5+LHX9pbTPa9grINxHe/20X1NyzRab7WMsdXDbO+63wnCZ4FKpjlA30bbfTz7miuTraqqBYjH6Ih0DE3T5gHzPJt6Rx5aXvv+AByOuHY88Kw0ee0Ip2dpfXD6igpjaUnAHyiPjIT0dGPxUXd6+D3QwiVWMxAF9A2TeCF8frZefhi0QfhRoJLpeiBDVdUBGEnzamBmkzJLgOuBb4ErgZVdcb9UCCGE8LeADHSvaZoTmA2sALYZu7SfVVWdo6rqJZ5ibwDJqqruAu4BHgxEbEIIIURnydi8YSKc4g2nWCG84g2nWCG84g2nWEHG5g01MgWbEEII0UmSTIUQQohOkmQqhBBCdJIkUyGEEKKTJJkKIYQQnRT2vXmDHYAQQoQp6c3rR+FeM1U6uqiquqEz5wd6Cad4wynWcIs3nGINt3jDKVY/xSv8KNyTqRBCCBF0kkyFEEKITurOyXRe20VCSjjFG06xQnjFG4ZU5RQAAAYiSURBVE6xQnjFG06xQvjFe0IL9w5IQgghRNB155qpEEII4RdhNzm4P6iqOhmYizHt4uuapv0lyCEdQ1XV/wOmAUc0TRvh2ZcELALSgX2AqmlacbBi9FJVtS/wJtAL41GleZqmzQ3FeFVVtQGrgUiM7/4HmqY96pkacCGQDGwArtU0zU/TZneOqqpm4HsgR9O0aSEe6z6gHHABTk3TTgnF74GXqqoJwOvACIzv7m+AHYRYvKqqDvXE5DUQeATj/11IxdqddbuaqeeX08vARUAWMENV1azgRnWc+cDkJvseBP6jaVoG8B9CZ4o6J3CvpmlZwKnAHZ6fZyjGWwuco2naaOAkYLKqqqcCTwHPaZo2GCgGbgpijE3dhTFtoVcoxwpwtqZpJ2madopnOxS/B15zgc80TcsERmP8nEMuXk3Tdnh+picBJwNVwEeEYKzdWbdLpsA4YJemaXs8f9EvBKYHOaZjaJq2Gihqsns6sMCzvgC4NKBBtUDTtFxN037wrJdj/EJKIwTj1TRN1zStwrNp9Sw6cA7wgWd/SMQKoKpqH2AqRu0JVVUVQjTWVoTc9wBAVdV44AyMeZTRNK1O07QSQjTeRs4Fdmuatp/Qj7Vb6Y7NvGnAwUbb2cD4IMXii16apuV61vMwmlVDiqqq6cAYYB0hGq+nZWIDMBijhWI3UOKZwB6M70NakMJr6nngfsDu2U4mdGMF4w+Tz1VV1YHXNE2bR4h+D4ABQD7wT1VVR2N8J+4idOP1uhp4z7Me6rF2K92xZhr2NE3TCbGhFFVVjQU+BO7WNK2s8bFQilfTNJenuawPRitFZpBDapaqqt575huCHYsPJmqaNhbjFsodqqqe0fhgKH0PMCoSY4FXtP+/vfsJ0aoK4zj+LUhoJIpQQrLJiGjRSsIgJJDIAgUXEg+llESLCIqgTeQmgxbRIqh9YpAmT2EoEbYZ3QSl/WMyq0X/bIxGCCbCXUGLc15nFMnivs49w/1+Nvee+76L38Cdee6cc+45mWuBs1zQTdpYXiJiGbAFeOfCz1rLOkRDLKangZsWtFfXa62bjYhVAPV4puc850TEVZRCujczD9TLzeYFqF16R4C7gesiYtRL08r9sB7YUif17Kd0775Gm1kByMzT9XiGMqZ3F+3eBzPATGZ+UtvvUoprq3mhPKR8npmztd1y1sEZYjE9DtwWEbfUJ72HgEM9Z/ovDgE76vkO4GCPWc6p43hvAN9k5qsLPmoub0SsrDM4iYirgY2UMd4jwIP1a01kzcznM3N1Zq6h3KNTmbmdBrMCRMTyiLhmdA7cD5ygwfsAIDN/A36pM2WhjEWepNG81cPMd/FC21kHZ3Bjppn5V0Q8BXxIeTVmd2Z+3XOs80TE28AGYEVEzAAvAC8DGRGPAz8D0V/C86wHHgG+iogv67WdtJl3FfBmHTe9EsjMfD8iTgL7I+Il4AvqpJRGPUebWW8A3osIKH9X9mXm4Yg4Tnv3wcjTwN76UP0D8Bj1vmgtb31A2Qg8seByi79jg+UKSJIkdTTEbl5JksbKYipJUkcWU0mSOrKYSpLUkcVUkqSOBvdqjDQOEbGH8m7fi5m5q980kvpmMZUuoq48dPO/fGU7MAd8vCiBJDXNYipd3G7g+nr+JLCMsmTiTL12LDP39RFMUntctEG6hIiYA66l7NV5tF7bw4Ju3ojYRVmp6jDwO7CVspzeNsrOL9soGzg/OtqyLiImgVeAe4DlwDHg2cw8sUg/mqQxcQKSNF4PULZKOwWsAz6lbEs3DdwBvA4QERPAFGUJuGnKOqsbgKmIWLHoqSV1YjGVxut7YBNl3VSACeA+5tdUXVuPm4FbgV+B7yibwZ8CVjK/kL2kJcJiKo3Xt3Vvybnans3MP4A/a3uiHtfU442UTamfoRRXKBuXS1pCnIAkjdffl2iP/FSPnwHragGmbhF3xeWJJulysZhK/fgA+BG4E/goIqaBScq46SbgaG/JJP1vdvNKPcjMs8C9lM2eJykzg28H3qKMoUpaQnw1RpKkjvzPVJKkjiymkiR1ZDGVJKkji6kkSR1ZTCVJ6shiKklSRxZTSZI6sphKktSRxVSSpI7+AXHlYL8exIYUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sline = plt.plot(\"t\",\"S\",\"\",data=sir_out,color=\"red\",linewidth=2)\n",
    "iline = plt.plot(\"t\",\"I\",\"\",data=sir_out,color=\"green\",linewidth=2)\n",
    "rline = plt.plot(\"t\",\"R\",\"\",data=sir_out,color=\"blue\",linewidth=2)\n",
    "plt.xlabel(\"Time\",fontweight=\"bold\")\n",
    "plt.ylabel(\"Number\",fontweight=\"bold\")\n",
    "legend = plt.legend(title=\"Population\",loc=5,bbox_to_anchor=(1.25,0.5))\n",
    "frame = legend.get_frame()\n",
    "frame.set_facecolor(\"white\")\n",
    "frame.set_linewidth(0)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}