{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forward Euler, continued\n",
    "\n",
    "First, revisit the simulations with large time steps. On Tuesday/Wednesday, we determined that, in order for the simulation to be constantly shrinking, that the time step should be < 1. \n",
    "\n",
    "Let's try that with our single-state degradation model. Here is copy-pasted code from Tuesday/Wednesday."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6F0lEQVR4nO3deXxU5b348c83O1sgQxaWsAXZQQQTtKUqClqtVixqxV57qdXa3lu7b7b9tfZae69t/V3b/tpbsS7F1qXWFW9dKgSXVkGiopBEloR9SQKBsCck8/z+eObAECaTmcyZOTPM9/16zWtmzpxz5hmW+c55nu/zfcQYg1JKqfSV4XUDlFJKeUsDgVJKpTkNBEopleY0ECilVJrTQKCUUmkuy+sG9ERhYaEZOXKk181QSqmU8s477+w2xhR13p6SgWDkyJFUVVV53QyllEopIrI51HbtGlJKqTSngUAppdKcBgKllEpzGgiUUirNaSBQSqk050ogEJEHRaRRRNZ08bqIyG9EZIOIfCAi04NeWyAi6wO3BW60RymlVOTcuiL4I3BpmNcvA8YEbrcAvwcQER9wO3AOMAO4XUQKXGqTUkqpCLgSCIwxrwPNYXaZCzxsrOXAABEZDHwceMUY02yM2Qu8QviAEpN3/nchK574ZbxOr5RSKSlRYwRDga1Bz7cFtnW1/RQicouIVIlIVVNTU48akVm7mKKaP+L36xoMSinlSJnBYmPMfcaYcmNMeVHRKTOkI9J70BiGml2s2Rbu4kUppdJLogLBdmBY0PPSwLautsfF0LJJ5Eo7K94POaatlFJpKVGBYDHwr4HsoXOBFmPMTuBl4BIRKQgMEl8S2BYXfQaPAaB+7ep4vYVSSqUcV4rOichjwCygUES2YTOBsgGMMfcCLwCfADYAh4EbA681i8hPgZWBU91hjIlfv42vDACzdyM79h1hyIBecXsrpZRKFa4EAmPM9d28boAvd/Hag8CDbrSjW/lDMRk5jJQGln7YyGfPHZGQt1VKqWSWMoPFrsjIhIIRTMhtYklNg9etUUqppJBegQAQXxnjc3bzVt0eDrW2e90cpZTyXNoFAnxlFLZtp62jgzfW7/a6NUop5bm0DASZHUcoyzvI0lrtHlJKqbQMBABXDmul8sNGOnSWsVIqzaVhIBgFwAVFB9lzqI1VW/d52x6llPJY+gWCAcNBMpmQu5usDNHuIaVU2ku/QJCZDQOGkbd/ExUjfSzRQKCUSnPpFwjAjhM01zNnYgnrGg6ytfmw1y1SSinPpHEg2Mic8baKqV4VKKXSWfoGgtYWRvRu44ziviytbfS6RUop5Zn0DQQAzfXMnlDM8vo97D96zNs2KaWUR9IzEBTYFFKa67l4QgntfsPr63q26plSSqW6NA0EIwGB5nqmDS/A1ydHu4eUUmkrPQNBdh7kD4XmejIzhFnjiqj8sJH2Dr/XLVNKqYRLz0AAdoZxcz0AF08ooeXIMd7ZvNfjRimlVOKlcSCwKaQA540tIiczg6UfaveQUir9pHcgOLwbjrbQNzeLc8p0lrFSKj2lcSBwMofsVcGcCSXUNx2ivumgh41SSqnEcyUQiMilIrJWRDaIyG0hXr9HRFYFbutEZF/Qax1Bry12oz0RCZpLADB7QjGAZg8ppdJOzIvXi0gm8DvgYmAbsFJEFhtjapx9jDHfCNr/K8C0oFMcMcacFWs7ohY0lwCgtKA34wf1Y0ltA184vyzhzVFKKa+4cUUwA9hgjKk3xrQBjwNzw+x/PfCYC+8bm9y+0LfkeNcQ2O6hqs172Xe4zcOGKaVUYrkRCIYCW4OebwtsO4WIjABGAZVBm/NEpEpElovIVV29iYjcEtivqqnJpVnAvjLYeyIQzJ5QTIff8OpanWWslEofiR4sng88aYzpCNo2whhTDnwG+JWIjA51oDHmPmNMuTGmvKioyJ3WBMpRO6aWDqCwb65mDyml0oobgWA7MCzoeWlgWyjz6dQtZIzZHrivB17l5PGD+CoYBQd2QtshADIyhNnji3ltXRNt7TrLWCmVHtwIBCuBMSIySkRysF/2p2T/iMh4oAB4K2hbgYjkBh4XAjOBms7Hxo2TQrp30/FNsycUc+BoOys3NSesGUop5aWYA4Exph24FXgZqAWeMMZUi8gdInJl0K7zgceNMSZo2wSgSkTeB5YBdwVnG8VdpxRSgI+NKSQ3K0O7h5RSaSPm9FEAY8wLwAudtv240/OfhDjuTWCKG23oEd/JKaQAvXOymHlGIUtqG/jxFRMREY8ap5RSiZG+M4sBehVAL99JKaRgu4e2Nh9hfaPOMlZKnf7SOxDAKZlDALPHlwC6lrFSKj1oIAiqQuoY1D+PKUP7a7kJpVRa0EDgGwUtW6G99aTNsycU8+6Wvew52NrFgUopdXrQQOArAwzs3XzS5jkTSjAGKnWNAqXUaU4DQYgUUoBJQ/IZlJ+XOt1DW1bAbyvgsM5/UEpFRwOBEwj2njxOICLMnlDM6+ubOHqsI8SBSWbDEti9DrYs97olSqkUo4Gg90DIzT/ligBs99Dhtg6W1+/xoGFRagzMw9u20tt2KKVSjgYCkZMWsg/2kdED6ZWdmRrdQw3V9l4DgVIqShoIwBafCxEI8rIzOW9MIUtrGzi5MkaSaTtk6yVlZMH2d8GfAl1ZSqmkoYEA7DjBvi3QceyUl+ZMKGFHy1Fqdx7woGERavoQMDDuE3DsEDTWet0ipVQK0UAANhD42+18gk4uHF+MSJLPMm4IjA9MX2DvtXtIKRUFDQQQlEK68ZSXivrlMrV0AEuTORA01kBWLxh9oR383lbldYuUUilEAwF0OZfAcfHEEt7f1kLj/qMJbFQUGqqheDxkZEJphV4RKKWiooEAoN8g+4s6xBUB2HITAEuTdZZxYy0UT7KPS8th91o4ss/TJimlUocGAgibQgowrqQfQwf0Ss7uoUO74VAjlEy0z0sr7P32d7xrk1IqpWggcIQoR+0QES6eWMI/NuxOvlnGzvyB4gn2fsh0QHScQCkVMQ0EDt8om4vvD71o/ewJxRw95uefG3Yntl3dcWYUO11Defk2KOg4gVIqQq4EAhG5VETWisgGEbktxOufE5EmEVkVuN0c9NoCEVkfuC1woz094iuDjlY4sCPky+eMGkjf3KzkSyNtqLaZQn2LT2wrLbeBIJknwSmlkkbMgUBEMoHfAZcBE4HrRWRiiF3/Yow5K3C7P3CsD7gdOAeYAdwuIgWxtqlHuskcysnK4PyxhSytbcTvT6Iv2MZaKJ5oxzkcpRVwdB/sqfOsWUqp1OHGFcEMYIMxpt4Y0wY8DsyN8NiPA68YY5qNMXuBV4BLXWhT9LoJBGBnGTceaGXNjpYENaobfv+JQBDMGTDW7iGlVATcCARDgeApudsC2zq7WkQ+EJEnRWRYlMciIreISJWIVDU1NbnQ7E7yh0JGdthAcOG4YjIEltQkSffQvs22pERJp0BQOM5WVNVAoJSKQKIGi58HRhpjzsT+6l8U7QmMMfcZY8qNMeVFRUWuN5CMTCgYGTYQFPTJ4ewRBSxJlmqknQeKHRkZMHS6BgKlVETcCATbgWFBz0sD244zxuwxxjiL/94PnB3psQkVYiH7zmZPKKFm53527DuSoEaFcTwQjD/1tdIKO5DcdiixbVJKpRw3AsFKYIyIjBKRHGA+sDh4BxEZHPT0SsApj/kycImIFAQGiS8JbPOGEwjCZNvMmVACJMks44YaGDACcvud+lppBZgO2LEq4c1SSqWWmAOBMaYduBX7BV4LPGGMqRaRO0TkysBuXxWRahF5H/gq8LnAsc3AT7HBZCVwR2CbN3xlts/9YNdf8qOL+jByYO/kGCdorDl1oNgxtNzea/eQUqobWW6cxBjzAvBCp20/Dnr8feD7XRz7IPCgG+2IWXDmUL+SkLvYtYxL+NNbmznU2k6fXFf+CKPX3gq718P4y0O/3meg/TwaCJRS3dCZxcF8o+x9mAFjsN1DbR1+3ljv4Szj3ets109XVwRwohKpTixTSoWhgSBY/2Egmd0GgvKRBeTnZXlbhM5ZhaxkUtf7lFbAwQZo2ZaYNimlUpIGgmBZOTBgWLeBIDszg1njiqn8sJEOr2YZN1TbeQ8Dz+h6n1IdJ1BKdU8DQWe+MtgbPoUUbBG6PYfaWLV1X/zbFEpjDRSOhczsrvcpmQxZealfidQYqLzzRKVVpZSrNBB05iuDPfXd9qvPGltMVoZ41z3UUHPqjOLOMrNhyLTUvyJorIHXfwmv3+11S5Q6LWkg6MxXBq0tcGRv2N36986mYqTPm2qkR/bB/m3hB4odpeWw832bZZSq6irt/bqXdIKcUnGggaCzgsgyh8B2D61rOMjW5sNxblQnTR/a+3ADxY7SCltee9ea+LYpnuoq7VKixw7bYKCUcpUGgs4iqELqcGYZJ/yq4PiqZJFcEaR4JdJjR2DzmzD9s9B3EKx52usWKXXa0UDQWcFIQCIKBCML+3BGcV+WJroIXWONrS7av7T7ffOH2MqqqRoItrwF7UfhjIth0qdg/StwdL/XrVLqtKKBoLPsPPvF2U3xOcfsCcUsr9/D/qPH4tywIA01djnK4MVownFWLEtFdZU2TXbkTJg8z3ZzrX2h++OUUhHTQBCKb1REVwRgu4fa/YbX18VhjYRQjIHG6si6hRylFXbtgjA1lJJW3TIYfi7k9LGfo/9wWPOU161S6rSigSAUX1nEgWD68AIKemcnrnvowE442hLZQLHj+DhBis0nOLALGtbA6IvscxGYdJW9SjjsXW1CpU43GghC8Y2Cw7vtF243MjOEC8fbWcbtHf74t63BWYNgQuTHDJ4KGVmp1z1U/6q9dwIBwOSrwd8Otc970iSlTkcaCEI5njkU2TjBnAkltBw5xjubw889cEVjFBlDjuxeMGhK6gWCukroPRAGnXli2+Cp9u+nWrOHlHKLBoJQokghBThvTCHZmZKYxWoaaqDfYOjti+640grY/i74O+LTLrf5/XZ8oOxCu/SmQwQmzYONr6fmmIdSSUgDQSjOpLIIag4B9MvL5tyygYmZTxDtQLGjtMIuuuNULU12jdVwqPHkbiHH5KvB+KHmucS3S6nTkAaCUHL7Qt+SiK8IwHYP1Tcdor7pYPza1dEOTeu6rzEUSqpVInXKSoy+8NTXSiZC0XidXKaUSzQQdCWCheyDzZ5QDBDf7KHmeptH35MrgoJRtr89VTKH6pZB0QQ7IS6UyVfbyWYt2xPbLqVOQ64EAhG5VETWisgGEbktxOvfFJEaEflARJaKyIig1zpEZFXgtrjzsZ4piHwuAUBpQW/GD+oX3+6hngwUO0ROrFiW7JyyEqG6hRyT5gEGap5NVKuUOm3FHAhEJBP4HXAZMBG4XkQ6f1O9B5QbY84EngR+EfTaEWPMWYHblSQLX5nN2Y+i2uWcCSVUbd7LvsNt8WlTQw1IBhSN69nxpeWwe62tXprMNr9pr3zCBYLCM2w20enQPbTzA9j0T69bodKYG1cEM4ANxph6Y0wb8DgwN3gHY8wyY4xTonM5EEGRHI856xfv3RTxIbMnFNPhN7y6Nk6zjBtrwDfapoP2hDOxbPs77rUpHuoqITMHRnw0/H6T58H2qqj+jpKO3w9//Rw8+mk4mKDZ6Up14kYgGApsDXq+LbCtKzcBLwY9zxORKhFZLiJXudAed0Q5lwBgaukACvvmxq97qDGCxWjCGTIdkOQfJ6hbBsM/Ajm9w+83aZ69r34m/m2Kl7pKaK6DtoPwhi68o7yR0MFiEbkBKAd+GbR5hDGmHPgM8CsRGd3FsbcEAkZVU1MCfjn5Il+XwJGRIcweX8xr65poa3d5lnHbIRuUejI+4MjLtzOSk3mc4MAuOxYSrlvIUTAChpandu2htxfaDLWpn4GVD8DezV63SKUhNwLBdmBY0PPSwLaTiMgc4IfAlcaY48tlGWO2B+7rgVeBaaHexBhznzGm3BhTXlRU5EKzu9GrAHr5ogoEYLuHDhxtZ+Uml2vhNH0ImNgCAZyoRNrNUpyeqVtm7yMJBGCzh3atht3r49emeNlTB+v/DmffCLN/BBmZsOw/vW6VSkNuBIKVwBgRGSUiOcB84KTsHxGZBizEBoHGoO0FIpIbeFwIzARqXGiTO6KoQur42JhCcrIy3O8ecmoMRVNsLpTSCji6z34JJaO6SuhTBCWTI9t/0lWApOag8dt/sCW2y2+0abLnfBE++EtqryanUlLMgcAY0w7cCrwM1AJPGGOqReQOEXGygH4J9AX+2ilNdAJQJSLvA8uAu4wxSRQIoptLANA7J4uZo+0sY+Pmr+7GWrtcY8HI2M6TzCuW+f1QH6KsRDj5Q+yg8pqnkvcqJ5TWA7DqERvI+g2y22Z+3XbfVf7Uy5apNOTKGIEx5gVjzFhjzGhjzM8C235sjFkceDzHGFPSOU3UGPOmMWaKMWZq4P4BN9rjGl8ZtGyNeuH3ORNL2Np8hPWNLs4ybqyG4vG2+yAWhePs6mbJGAga1sChpsi7hRyTPmXTYhuT5zdEt95/HFr3wzlfOrGtt88Gg3Uvwea3PGuaSj86szgcXxlgoh7Amz0+DmsZN9TEPj4A9pf20OnJGQjClZUIZ+JVdn5FqnQP+f3w9n02i8sp/eE450u2qOCS21PrCkelNA0E4TgppBEWn3MM6p/H5KH57pWbOLTbFmBzIxCA7R5qqI5qslxC1FVC8aQTXSWR6lsEo85Pne6h+mWwe50dE+gspzdc8D3YusJeGSiVABoIwomyHHWwORNKeHfLXvYcjK5bKaSGQGmJWOYQBCutANMBO1a5cz43tB22tYOivRpwTL7aBuydq1xtVly8fZ8dEJ/0qdCvT7vBThxc8h+pUzY8nCi7VlXiaSAIp/dAyOnX40BgDNz993UciHVhe6d0dHGMGUOOoUlYiXTzm9DRFv34gGP8FXYVtmTvHmquh3Uv25TRrNzQ+2Rm23TSplqbRZTK3vkj3DUCNizxuiUqDA0E4Yj0KIUUYNKQfOZXDOOxt7dwwS9fZdGbm3o+yayx2galvsU9O76zPgPt1U4yBYK6SsjM7b6sRFd6+2D0bDvLOJm7h96+3w74l38+/H4Tr4LBZ9l5BceOJqJl7tu7CV76AbQfhb9+3pZQV0lJA0F3oljIPpiIcNfVZ7L41pmMLenL7YurueSe13hh9c7o00qdgWKRqNvRJacSabJ8adZV2iDQ0zpKYGsPtWxNrgAXrPUgvPdnmDgX8geH31cE5vzEfp6qBxPSPFf5/fDcrTboff4lyMqx9ZQOuzzRUrlCA0F3fGWwb4tdFKYHziwdwGNfOJeHPldBTlYG//7Iu8z7/ZuRzzz2+23XkFsDxY7SCjjYAC3b3D1vT+zfYbtBetot5Bj3CXtVkawlJz74C7S2wIwQg8ShjL4QymbZGkRH98e1aa6regA2vQEf/xkMPxfmPwr7t8MT/wrtcarOq3pMA0F3fGXgb7e/zHpIRLhwfDEvfu18fn71FHbsO8K1977FLQ9XsaG7uQb7NtslJt0aKHYk04pl0ZaV6EpePoy5GKqfTb5BVmPsIPHgs2DYjMiPm307HN4Db/6/uDXNdc0b4ZUf2666aZ+124bNgCt/a4PDi99JnitRBWgg6F4MmUOdZWYI11UMZ9m3Z/HtS8byZt0ePv6r1/nhM6tpPNBFP7DbA8WOksmQlZcclUjrKqFPcezlM8BmDx3cZQefk8nG12y9qHO+GF0X39Dpdrzgrd/BwTiufueW411CWXDlb07+rFOvg4990w4gr1joWRPVqTQQdKcHVUi70zsni1svGsOr35nFDecM5y8rtzLrl6/y6yXrOdTaqQvq+Kpk4117f8BmpgyZ5v0VgVNWYvRF7oyBjP04ZPeG6iTLHlpxnx3wd0pnR+OiH9kB19d/2f2+Xlt5P2z+h+0S6h9i2ZGLfmQzvF7+PqxP8UwiY2D1k8lbtysKGgi603eQrfETZc2hSBT2zeU/5k7mlW9ewAVji7hnyTpm3f0qj67YQntHIMOooQYGDIfcfq6/P6XlsPN9b/O8d31guz5i7RZy5PSBcZdBzXM9Htdx3d5NsPYFmzKanRf98YVnwPR/haqH4vLv0DXNG+2M6DPmnOgS6iwjAz610F7hPnkjNK1NbBvd4vfDC9+Gp26ChRfY7sgUpoGgOxkZPU4hjdSowj78/oazeerfPsoIX29+8MxqLv31GyypacA01rjfLeQorbBLQnpZ7dIpK1E2y71zTppng8vG19w7ZyxW3m9LYHSXMhrOBd+z3S3JWqY6uEvok78Jf3WX2xeuf8x2TaZiJpG/Axbfav9eZ9xil4796wKbKtsR45whj2ggiEQPU0ijdfaIAv76pY+w8LNn4/cb/v3ht/A3rWNXr7L4vGEyVCKtq4SSKdCvxL1znjHHFtZLhsllbYfg3Ydhwiehf7iF+7qRPxjO/RKs/qtdfyHZHO8S+s/IPueAYYFMop3wl8+mTiZRxzF46mZbOXbWD+CyX8CNL9oaUct/B3+8HFpOWY4l6WkgiIRvlL2897u86lgIIsLHJw3i5W+cz/+d3YdM/PxspfDlR99l8x6XawPlD4H8od4FgrZDsGV5z8tKdCU7D8ZfDh8+7315gw+egKMtJ1cZ7amZX4e8/rb0RDJprg90CV1sy2NEalgFzP2tDSAvfDv5M4mOHbXpr9VPw8U/hVnfs1c+WTlw2c/hmodsOZiF55240k0RGggi4SuzXSgHdiTsLbMzM/hkyV4AymfMpLK2kTn//Rr/8Xw1zYdc/PXkrFjmhU3/BP8x98YHgk2aZ7+AvfwP6aSMDppic+lj1WsAfOwbsOEV2PSP2M/nhuNdQtnwyV9HP+B/5qfhvG/Du4tg+e/j00Y3tB2Gx+bbsZ5P3A0zv3rqPpPnwReW2Qy4P82D136RkB+PbtBAEIkC9zOHItJQDRnZLPjkxbz2nVlcc/YwFr25iQt+sYz/eXUDR4+5kCtfWmHnKniRmlhXafuJh3/E/XOXzbLLjXrZPbTpDbtGwjlfcm9W+DlfhH5DYMlPkuMX9Mo/wOZ/wqURdgmFcuEPbSbR338I619xt31uaD0Aj1xjx5zm/g5mfKHrfYvGwheWwpRrYdnP4NFrU2IMRANBJFycSxCVxhooHAuZ2RTn5/Ff86bw8tfP55wyH794aS0X3v0qf63aSoc/hi+E4+MEHswnqKuEETN7lknTnawc2y+/9gU4dsT980dixUK77vXkq907Z3YvmHWbvYr78G/unbcnmuttQDrjYjjrX3p+nowMmHefnUfy1xtPzJ1JBkf2wsNzbRfmvD9E1vWV08d+nivugY2vw73nJcd8nTA0EESif6m99E14IKg9ZUbxmJJ+3L+ggsdvOZfifrl858kPuPw3b/DauqaeLY05eKrN9Eh091DLNruqWDy6hRyTr4a2g3aB+ETbtyWQMrogtvpJoZz1LzBwDCy9w7sU2Vi7hDrL6QPXP27/rB69Dg7tcaedsTi0GxZ90g7OX/cnmHJN5MeK2Cyxm/5uA92Dl9q5JMlwFReCBoJIZGTatYITGQiOttiyFl3UGDq3bCDPfnkmv/3MNA63dbDgwbe59t63+MEzq/nvV9bxp7c28eLqnVRtambT7kMcbG0PHSiye9k+7EQHArfKSoQz4mO27r8XtYdW3g8IlN/k/rkzs2yZ6t1r4YPH3T9/JNzoEuqsf6lNKz2wC57wOJNo/0546BOwe71t0/jLe3aeIdPgltfsv/MXv2PnHbS6uIStS7LcOImIXAr8GsgE7jfG3NXp9VzgYeBsYA9wnTFmU+C17wM3AR3AV40xL7vRJtf5yqB5U+Le73hpia5rDIkIV5w5hEsmDuKRFZv5y8qtvLxmF82H20L+8MjLzqCwby6FfXMp6he475vDJ7LGccbWxbxX30Rhfm8K++bQNzcLcbPaaWd1lXayXvGE+L1HZpYtz/Den20/bzwm5YVy7IhNGR1/uU2TjIcJV8LQs2HZf8Hka+LTvdYVp0tozCWxdQmFUlpu++Gfvhn+9k248v+5W3U3Evu2wKIr7frZNzwFIz8W2/l6++zVzj/vgco77RXGp//kfrWAGMQcCEQkE/gdcDGwDVgpIouNMcErid8E7DXGnCEi84GfA9eJyERgPjAJGAIsEZGxxpgkqxiGDQSb/mEv7RLxDzOKVclysjK4ceYobpxpB7XbO/w0H2qj6WAruw+2sftAq318oJXdgW1bmw/z3pa97DnUxmbx8eucw/zoD0/yoRkOQG5WIGj0y6Woby5F/XKOBxF7yyG/VzZZGUJm0C0rIyNwL2Rm2vsMkeP7iYidkFO/DMZeFv8/y8nz7K/XtS/BmdfG970cq/9q+5ZDLUXpFqdM9aJP2s/30a/E772Cud0lFMqZ19qrndd/aX8ofOTL7r9HV/bU2SDQdgA++6xNcXVDRgac9y07Jvfk5+EPF9qJd4n6N9kNN64IZgAbjDH1ACLyODAXCA4Ec4GfBB4/CfxW7M/NucDjxphWYKOIbAic7y0X2uUu3yhbBfRgo7uTn7rSWGMnRfWP/hdlVmYGxfl5FOd3/yuxvcNPy47R8MD/8Jvz2lkzaCpNQQFj98FWtu09zKqt+2g+1Eos49JgC+9Nzajn6ay9fG9VIS+v/vvxIOEEkRNB5cR9RiCgAAgnvn+E4w+cRye9lgH8JqOQTS/czz1vDz/pdfv4xDmd16TTa+Gcsocx/GTHPZBdxk+W5YCsjOK7Mtov1d58vVc5I5f8nO+vncKRjL4RH9nT7++LWp7m+j3/5KGi7/DmszuA+KRUi7mEL/ZZzrSX/w+/fd+wurcL6bfdGNy2iW/u/C6Zpp17Bv+crcsEcHuQtzf9B/6WWxrvZOzTN7NsyWKeGPhvtEtOxGe4Y+4kBvd3d9zJjUAwFAiu0bwNOKerfYwx7SLSAgwMbF/e6diQHY4icgtwC8Dw4cNdaHaUgjOHEhIIau2voTj/Ys7KzGBg6TjoPZCxbR8ydnqIQmEBHX5D86G2QJBo5eDRdjqMocNvaO8I3PuN3dbht4+dbUG3iq2vwzYonHoJV2YWHN8evL/fb2j3+0/abgwYzPFur+P3wduc7QYMNof7zbzzuPTw80hrC4ekz8n7BR47JzScfO5wTux9wqRjaxh2rJ7f9PkKDQdbIx4b7OkY4r1Zn+W3R77GzIZHebhXF/V9Or9Xz96KwR07+NT++1mRXcETx86HvfHNxvpp9te4O3MHN++6k2/k382WzBFxe6/R7XV86+CP6CCT7/X7L7YcLoXD8fl82+jDN/Lu5HNmEdfuf5qhh2r5WZ/baMyM7HulvcP9AWdXxggSwRhzH3AfQHl5eeKH3oMDwYg45L0HM8Z2DXW1uLnbRE6sWBZGZoZQ1M+OL8TkodUw6Ey+M++82M4TqW25cP8z/OmjTTDtkvi+1xN/gI0FfPVrP+CrOb3j+14AnAdPvs78tc8z//N3QL9B8Xkbv9+WT2jN45wvP8yL+UPi8z6dtUyEP1zIwoxfwM2V0KfQ/ffYVgV//jH0y4cFi1k4cLT77xHSLKi9mrHP/huL2r5lU07HfjxB730yN7KGtgPB/RelgW0h9xGRLKA/dtA4kmOTQ/9hIJmJyRw6sBOO7nN/VbJwSsttv+yRffF9n9YDsHVFfLOFOhs6HQaMiH/2UMs2qP1fWyk0IUEg4MIfQkdbfMtUv70QtrwJl/6XLU2SKP2HwvzH7Gp68ahJtOkfdp5ALx98/kVIWBAImHAFfPE1m1Tw6Kdh6U89WVTJjUCwEhgjIqNEJAc7+Lu40z6LgQWBx9cAlcbmMi4G5otIroiMAsYAb7vQJvdl5di/rL0JKAPcEBhecXtVsnCciWU73o3v+8SzrERXROygcf2r8c1PX/kAYKDi5vi9RygDR8P0BXbBl3jUxt9TZ+sbjfk4nPUZ98/fndKzbSbRljfhf7/hXi7+hiXw52tsva0bX7Tl3r3gK4ObXrE/IN64G/50VcJn+sccCIwx7cCtwMtALfCEMaZaRO4QkSsDuz0ADAwMBn8TuC1wbDXwBHZg+SXgy0mZMeRIUBXSE4vRJDAQDJkOSPxnQNZV2vUd3Ki9E41J88B0QO1z8Tn/sSP2i3jcJ7z5Qrngu5CZY8sauMnvh+e+bM/9yV8lPpXTMeUaW4p71Z/hrd/Gfr4P/waPXQ8Dz4AbX7DVXb2U3cumys79H9j6Niw8HzYnLmfGlQllxpgXjDFjjTGjjTE/C2z7sTFmceDxUWPMtcaYM4wxM5wMo8BrPwscN84Y86Ib7YmbglGwpz7+swMba6HfYJt/nCh5+XZwOt4Ty+oqbV52VozjDNEaNMXOxo1X7aE1T8GR5vimjIbTbxCc+2+2HTtWuXfetxfClrfgsrsS2yUUygW3wcS58Pcf2XTgnlrzlO1mGjQFPvd8fMYdemrav8DNS2xg+OPldq3qBMxG1pnF0fCVQWuLzRGPp4bqxF4NOJxKpPH6h7dvC+xZn9huIYfTPbTpH3bmqpuMsXWFiifCyAQNgIcy82u20N7SO9w5X3CX0NTr3TlnLDIy4Kp7bVmUp2460YUajfcesesJDDvHzhPoVeB6M2M2aArc8iqM/wT8/f/AX26wlQbiSANBNBJRfK6j3S7fF88Zt10prbBBLl6fLxFlJcKZNA8wdhlLN21ZbpfcnHGLd10nYNcqOO9bULfUFjuLhdMllJUTv4ljPZHT25Z8yOkLj11n6wFF6u0/wHP/DqMusDOG8/Lj185Y5fW3s48//p+w7iW4b1ZcFyTSQBCNRASC5nq79kFJnJanDCfeK5bVVdoSykXj4nP+7hSPt8t+up099PZC+x/3zE+7e96eqPiCHfyMtUz1inttl9ClP/e+/7yz/CFw/aN2QPUvN0S2+NA/f2MXvxl7mS33kMisrp4SsbOqP/c3OwZ1/xxbLiUONBBEo2AkIPFdQNyLgWJH4Tg7mzkegcDfYbN2Rl/k7a/LyfNs+uq+rd3vG4mW7VCzOJAy2sedc8YiOw9mfR+2vwO1z/fsHHvqbPfS2Eth6nx32+eWoWfDVb+3wSpcJpEx8OrP4ZUf2Xk51/0psXWZ3DD8XPjiG7Y767kvw1b3/39qIIhGdp79tRXPK4LGWrvQuRe/mjMybM59PALBjlV2boTby1JGa/I8e1/9jDvnq3oQjD/xKaPhTL3eBvWelKn2d8Cz/267hK74VfJ0CYUyeZ4dQF71CLz5m1NfN8Yuofnqf8LUz8DVD0BmduLb6Ya+RfDZZ+zVjFv1j4JoIIiWb1R8A0FDte2CcruGfaRKK2DXGrs0n5vqKgGBMo8Dga/MlgaudiF76NjRQMroZYGrxSThlKnes95+SUZjxb2wdXlydgmFcsH37C/9V26HtUFJh34/vPhd+OevbSnwub+z5eRTWUam/bcWj1PH5ayns3gHgsYab7qFHKUVNt9+5yp3z1tXabM9+gx097w9MWke7Hgv9slX1c/A4d3epYyGM/4K+3f56l2Rr9CWCl1CnWVk2Nz7IWfZbKCGantV8/xX7HrRH7kVLv+/dj/VJf3TiZavzP7nj0c6V9shO/7gxUCxY2i5vXeze+joftj2tnfZQp05NZxi6R4yxv56Lhpvs1CSjVOm+sAO+4XYneNdQrnJ3yXUWU5vmP+oXW/i0fm2zPN7f7ZXC5fcmVqfxSMaCKJ1PHMoDgPGTWsB4+0VQZ+B9jO6GQg2/QP87ckTCAYMswNvsUwu27bSXjXN+ELyftGM/BicMQfe+O/ua0g5XUKX/SI1uoQ6yx9ig8GhJqh51gbBC3+QvH83SUYDQbScQBCPmkONTo0hD68IwHYpbHVxYlldJWT3gWEz3DmfGyZfbTO0Gj/s2fErFkJufzgzybtQZt9uB+n/+auu99m9IdAldBmceV2iWua+odPt/IBrF8HHvuF1a1KKBoJoFdhVwOIyTtBQY+vweD3wWFoBB3fBfpcKwXpVViKciXMB6dmg8f6d9lfntBsgN/LFYDwx+EyYci0sv9e2uzN/h51klZXrbS0ht4ycCZOu8roVKUcDQbRy+0Kf4vgEgsZqmzbqdXZDqYvjBHs3QXNd8nQLOfoNssFpzdPRX/m885D9Ap2RRCmj4Vz4Q1vx9bWfn/ra8t/beRWX/SJ+axmopKeBoCd8ZfEZI2io8b5bCKBkMmTluVOJ1OuyEuFMnmdTLBvWRH5MeytUPWQXEHG6CZOdbxScfSO8+7DtBnLsXg+VP039LiEVMw0EPRGPctSHdsOhRm8Hih2Z2TbX3o0rgrpKyC+FwjGxn8ttE+baxYaiKTlR/az9e5pxS9yaFRcXfNcG92V32uf+jkAtobzTo0tIxUQDQU/4yuwqYm5Oujo+UJwEgQBs99COVbGtCNXRDhtfs7OJk/GLps9AKJsVXffQ2wuhcGxyXuGE07fY1q2pfsbOodAuIRVEA0FP+AIDxns3uXdOp6RuMlwRgB0w7miFhhgqHu54z863SOYvzcnzYN9m2B7BymzbqmwNH6+rjPbUR79il2R8/uu2S2jcJ5KjUJ7ynAaCnohHFdLGavuftG+Je+eMxfFKpDGMExwvKzHLjRbFx/grICM7suyhFQshp1/qzLrtLC8fzv+2nf+QlQdX3JOaAU25TgNBT/jikELqDBQny3/M/CG2wF4s4wR1lXasIZErrUWr1wA76WrN07Y+TVcONNhulWk32Bmsqar8Jph4FXxqoXYJqeNiCgQi4hORV0RkfeD+lOV+ROQsEXlLRKpF5AMRuS7otT+KyEYRWRW4nRVLexKmV4G9uRUI/H5o+jB5uoUczoplPXG0xR6bzN1CjslX21IMW1d0vc87D9kUzBlfSFy74iE7Dz69CMZd6nVLVBKJ9YrgNmCpMWYMsDTwvLPDwL8aYyYBlwK/EpEBQa9/xxhzVuC2Ksb2JI6bmUMtW6DtYPIMFDtKK+w4yMGm6I/d+IYtXpcKgWDcpbarpKvsofY2W256zCUwcHRi26ZUAsQaCOYCiwKPFwFXdd7BGLPOGLM+8HgH0AgUxfi+3nNzLkGyDRQ7nHGC7T0YJ6irtMsJlrpfO911uf3svICaZ0PX769dDAcbYEYSVhlVygWxBoISY4wzb30XEHakU0RmADlAcP3fnwW6jO4RkSSqQdANXxns3xbZMnndOb4qmQfrFIczeCpkZPWse6iu0i7knpXjfrviYdI8W7Bs8z9OfW3FveAbnRpXN0r1QLeBQESWiMiaELe5wfsZYwzQZTK2iAwG/gTcaIxxRuW+D4wHKgAf8L0wx98iIlUiUtXU1IOuCrf5yuzKVPu2xH6uhhoYMDz5BiGze8GgKdEHguZ6W5Qvlb44x1xir2A6VyTd/q79/Od8UWvaq9NWt/+yjTFzjDGTQ9yeAxoCX/DOF31jqHOISD7wN+CHxpjlQefeaaxW4CGgy/KUxpj7jDHlxpjyoqIk6Flys/hcY61dVD0ZlVbYL0N/R+THJHNZia7k9LarP9Uuho5jJ7a/fZ8NEFOv965tSsVZrD9xFgMLAo8XAM913kFEcoBngIeNMU92es0JIoIdX4ii6IvH3JpL0N5m690k20Cxo7TCDmQ3RVGuuX4Z9B+eegOrk6+GI3uh/lX7/GCTHUA+6zM2B1+p01SsgeAu4GIRWQ/MCTxHRMpF5P7APp8Gzgc+FyJN9BERWQ2sBgqBO2NsT+L0KbSTi2INBLvX2UVbkm2g2BFtJdKOdqh/PXnLSoQz+iK7xoCTPfTOH6GjLfXqCikVpaxYDjbG7AFmh9heBdwcePxn4M9dHJ9CfQediLizfnFjkmYMOQpGQe+BNhCc/bnu99/xLrQmeVmJrmTlwoRP2u6h1gNQ9QCMnp2cBfOUcpGOfsXCjRTShmpb4iBZv2xEbPdQpKUm6ipBMmDU+fFtV7xM/hS07ofnv2YLC57zJa9bpFTcaSCIha/MFiwLlXseqcZaW80yM9u9drmttNyOEXS37i0EykpMT+6yEuGMusBeAa15yv79njHH6xYpFXcaCGLhG2X791u29vwcjTXJN3+gM2dS2I5uKnQe2WevHFKxW8iRmQ0TrrSPK76gKaMqLei/8ljEmjl0tMUGkWTNGHIMmQ5I991Dm5yyEhcmpFlxc84XbVXSaTd43RKlEkIDQSxiDQSNtfY+WecQOPLyoWh895lDqVRWIpziCTD/EU0ZVWlDA0Es+g6CrF49HzBOtlXJwnEqkYZbyauu0g4SJ/N4h1LqFBoIYpGRYccJ9vYwEDTU2LkI/Ye52654KK2wk626uvpprreVSlN5fECpNKWBIFaxlKN2BopTYeLV8RXLuugeqqu09xoIlEo5GghiVTDSdg2FW90qFGPsHIJU6BYCKBpnr166DATLbOE8Z9xEKZUyNBDEyldmF3k/sCO64w7shKP7kn+g2JGRCUOnhw4EHcdg4+v2aiAVrm6UUifRQBCrnmYOpdJAsaO0AnatgbbDJ2/f/o6djavdQkqlJA0EseppIEjWVcnCKa2w8wR2rjp5e6qXlVAqzWkgiFX/UlsrKNoU0sYam36aSqUYuqpEWlcJQ8+GXgWJb5NSKmYaCGKVkRkYMI72iiCFBoodfQptNdLgQHBkr+0a0m4hpVKWBgI3+EZFd0XQ0Q5Na1OrW8hRWgFbgyaWbXzdLtmpgUCplKWBwA3OXIJws26D7d1oM41KUiRjKFhpBRzcBfu32+d1lZCbb7uGlFIpSQOBG3xlcOwQHAy5ZPOpGqrtfUpeEQSNExgDG7SshFKpTgOBG6LNHGqssVk2RePi16Z4KZkMWXm2EmlzPbRsSf1qo0qlOQ0EbnACQaQ1hxqq7THZveLXpnjJyoHBZ9krAi0rodRpIaZAICI+EXlFRNYH7kPmD4pIR9DC9YuDto8SkRUiskFE/iIiObG0xzP9h4FkRndFkIrdQo7SctixCta9ZDOmtKyEUikt1iuC24ClxpgxwNLA81COGGPOCtyuDNr+c+AeY8wZwF7gphjb442sHDufIJJA0HbYZhil4kCxo7TCDnZvWKJXA0qdBmINBHOBRYHHi4CrIj1QRAS4CHiyJ8cnnUirkDZ9CJjkX54ynOCFZzQQKJXyYg0EJcaYnYHHu4CSLvbLE5EqEVkuIlcFtg0E9hljnJXftwFDu3ojEbklcI6qpqamGJsdB74y2BNBCqlTYyhVis2F0n8o9Btiu8NGnud1a5RSMcrqbgcRWQIMCvHSD4OfGGOMiHT1LTjCGLNdRMqAShFZDbRE01BjzH3AfQDl5eURJuwnkK8MWlvsTNtwZSMaauyqZr5RiWtbPIy/3FZQ7TXA65YopWLUbSAwxszp6jURaRCRwcaYnSIyGAiZSG+M2R64rxeRV4FpwFPAABHJClwVlALbe/AZksPxFNKN4QNBY7VNG83ITEy74uXyu71ugVLKJbF2DS0GFgQeLwCe67yDiBSISG7gcSEwE6gxxhhgGXBNuONTRqRzCRprU3ugWCl12ok1ENwFXCwi64E5geeISLmI3B/YZwJQJSLvY7/47zLGBDrK+R7wTRHZgB0zeCDG9ninYIS9DxcIDu2Bgw2pPVCslDrtdNs1FI4xZg8wO8T2KuDmwOM3gSldHF8PzIilDUkjuxfkDw0fCBpTuLSEUuq0pTOL3dRdCqmzGI12DSmlkogGAjf5RnVzRVADvXzQt6ssW6WUSjwNBG7ylcHh3XB0f+jXG2vs1YAu8K6USiIaCNwUrvic328zhnSgWCmVZDQQuKkgMEksVPdQyxZoO6gDxUqppKOBwE2+MIFAB4qVUklKA4GbcvtBn+LQgeB4jSHtGlJKJRcNBG7zlYVeyL6xBgYMt8FCKaWSiAYCt3UVCBpSfDEapdRpSwOB23xlcGCHXYDG0d4Ge9ZrIFBKJSUNBG5zBoz3bjqxbfc68LfrQLFSKilpIHBbqMyhxlp7r1cESqkkpIHAbaHKUTdWQ0Y2DDzDmzYppVQYGgjc1qvA3oIDQUMNFI6xi9wrpVSS0UAQD76yk8tMNGrGkFIqeWkgiIfgctRHW6BlK5RoIFBKJScNBPFQMApatkF7KzR+aLcVa8aQUio5aSCIB18ZGD/s2xK0KpmWllBKJaeYAoGI+ETkFRFZH7gvCLHPhSKyKuh2VESuCrz2RxHZGPTaWbG0J2kEZw411EBOP1teQimlklCsVwS3AUuNMWOApYHnJzHGLDPGnGWMOQu4CDgM/D1ol+84rxtjVsXYnuQQHAgaa+zVgC5Go5RKUrEGgrnAosDjRcBV3ex/DfCiMeZwN/ultj6F9iqguR4aqnWgWCmV1GINBCXGmJ2Bx7uA7hbjnQ881mnbz0TkAxG5R0RyuzpQRG4RkSoRqWpqaoqhyQkgYmcYb34Lju7TgWKlVFLrNhCIyBIRWRPiNjd4P2OMAUyY8wwGpgAvB23+PjAeqAB8wPe6Ot4Yc58xptwYU15UVNRds73nGwUNq+1jHShWSiWxrO52MMbM6eo1EWkQkcHGmJ2BL/rGMKf6NPCMMeZY0Lmdq4lWEXkI+HaE7U5+zjgBaLE5pVRSi7VraDGwIPB4AfBcmH2vp1O3UCB4ICKCHV9YE2N7kocTCPoOgt4+b9uilFJhxBoI7gIuFpH1wJzAc0SkXETud3YSkZHAMOC1Tsc/IiKrgdVAIXBnjO1JHk4g0IFipVSS67ZrKBxjzB5gdojtVcDNQc83AUND7HdRLO+f1JxAoDWGlFJJLqZAoMLoNxgu/CFMvMrrliilVFgaCOJFBC74rtetUEqpbmmtIaWUSnMaCJRSKs1pIFBKqTSngUAppdKcBgKllEpzGgiUUirNaSBQSqk0p4FAKaXSnNjq0alFRJqAzT08vBDY7WJzUkU6fu50/MyQnp9bP3NkRhhjTqnjn5KBIBYiUmWMKfe6HYmWjp87HT8zpOfn1s8cG+0aUkqpNKeBQCml0lw6BoL7vG6AR9Lxc6fjZ4b0/Nz6mWOQdmMESimlTpaOVwRKKaWCaCBQSqk0l1aBQEQuFZG1IrJBRG7zuj3xJiLDRGSZiNSISLWIfM3rNiWKiGSKyHsi8r9etyVRRGSAiDwpIh+KSK2IfMTrNsWbiHwj8G97jYg8JiJ5XrcpHkTkQRFpFJE1Qdt8IvKKiKwP3Bf09PxpEwhEJBP4HXAZMBG4XkRO9wWF24FvGWMmAucCX06Dz+z4GlDrdSMS7NfAS8aY8cBUTvPPLyJDga8C5caYyUAmMN/bVsXNH4FLO227DVhqjBkDLA0875G0CQTADGCDMabeGNMGPA7M9bhNcWWM2WmMeTfw+AD2i2Got62KPxEpBS4H7ve6LYkiIv2B84EHAIwxbcaYfZ42KjGygF4ikgX0BnZ43J64MMa8DjR32jwXWBR4vAi4qqfnT6dAMBTYGvR8G2nwpegQkZHANGCFx01JhF8B3wX8HrcjkUYBTcBDgS6x+0Wkj9eNiidjzHbgbmALsBNoMcb83dtWJVSJMWZn4PEuoKSnJ0qnQJC2RKQv8BTwdWPMfq/bE08icgXQaIx5x+u2JFgWMB34vTFmGnCIGLoKUkGgT3wuNggOAfqIyA3etsobxs4D6PFcgHQKBNuBYUHPSwPbTmsiko0NAo8YY572uj0JMBO4UkQ2Ybv/LhKRP3vbpITYBmwzxjhXfE9iA8PpbA6w0RjTZIw5BjwNfNTjNiVSg4gMBgjcN/b0ROkUCFYCY0RklIjkYAeVFnvcprgSEcH2GdcaY/7b6/YkgjHm+8aYUmPMSOzfcaUx5rT/lWiM2QVsFZFxgU2zgRoPm5QIW4BzRaR34N/6bE7zAfJOFgMLAo8XAM/19ERZrjQnBRhj2kXkVuBlbHbBg8aYao+bFW8zgc8Cq0VkVWDbD4wxL3jXJBVHXwEeCfzQqQdu9Lg9cWWMWSEiTwLvYjPk3uM0LTUhIo8Bs4BCEdkG3A7cBTwhIjdhy/J/usfn1xITSimV3tKpa0gppVQIGgiUUirNaSBQSqk0p4FAKaXSnAYCpZRKcxoIlFIqzWkgUEqpNPf/AZwTdW5VpvK8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Run the 1D model to see where the steps break down..\n",
    "import numpy as np\n",
    "import scipy.integrate\n",
    "\n",
    "def degradation_model1D( t, y, params ):\n",
    "    dydt = -params[0] * y;\n",
    "    return dydt\n",
    "\n",
    "def forwardEuler1D( fn, t, y0, params ):\n",
    "    \"\"\" [t,y] = forwardEuler1D( fn, t, y0, params )\n",
    "    Given a function handle fn for dydt, find the value of\n",
    "    y at timepoints t, given y(0) = y0.\n",
    "    fn takes t, y, and params as input.\n",
    "    t is 1-dimensional (NT,)\n",
    "    y is 2-dimensional (1,NT) (NOTE WE ASSUME WE HAVE JUST ONE VARIABLE)\n",
    "    \"\"\"\n",
    "    y = np.zeros( (1,t.size) )\n",
    "    y[0,0] = y0;\n",
    "    for i in range(1,t.size):\n",
    "        h = t[i]-t[i-1];\n",
    "        yprime = fn(t[i-1],y[0,i-1],params);\n",
    "        y[0,i] = y[0,i-1] + h*yprime;\n",
    "    \n",
    "    return (t,y)\n",
    "\n",
    "dt = 1\n",
    "y0 = 1\n",
    "params = (2,)\n",
    "t = np.arange(0,10,dt)\n",
    "tru_y = np.exp( -params[0]*t)\n",
    "(t,y) = forwardEuler1D( degradation_model1D, t, y0=y0, params=params)\n",
    "\n",
    "plt.plot(t, tru_y)\n",
    "plt.plot(t, y.T);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rewrite the forward Euler method to work for multi-state models\n",
    "\n",
    "Rewrite the method so that y0 can have more than one value and fn can return a vector of more than one value.\n",
    "\n",
    "The return value should be NY x NT (number of states by number of time steps) big, so it gets filled in one column at a time. Adjust the 1D code to account for multiple states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "def forwardEuler( fn, t, y0, params ):\n",
    "    \"\"\" [t,y] = forwardEuler( fn, t, y0, params )\n",
    "    Given a function handle fn for dydt, find the value of\n",
    "    y at timepoints t, given y(0) = y0.\n",
    "    fn takes t, y, and params as input.\n",
    "    t is 1-dimensional (NT,)\n",
    "    y is 2-dimensional (NY,NT)\n",
    "    \"\"\"\n",
    "    # Write code here\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, make a 2-state degradation model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def degradation_model( t, y, params ):\n",
    "    # Write code here\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the 2-state degradation model with forward Euler.\n",
    "\n",
    "Use 2 and 2.5 as the degradation rates. Simulate from t=0 to t=4 with a time step size of 0.1.\n",
    "\n",
    "Then plot the true solution and the approximate solution (indicate simulated with dashed lines). And in another plot, show the difference (between the true and approximate solutions) over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determine how the error depends on the step size\n",
    "\n",
    "When we compute the difference between the true and approximate solutions, we have one value at each state for each time step. I would like to summarize this as one number, so that we can plot \"error\" versus step size.\n",
    "\n",
    "Let's explore different measures for the error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.linalg\n",
    "\n",
    "# Write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine how the error changes with the time step. Plot two figures - one with the simulations and one with the errors. Compute the error using two methods:\n",
    "- Compute the Euclidean distance between the true solution and the approximate solution (2-norm of the difference) at each time point, then take the mean across time.\n",
    "- Compute the infinity norm (max value) of the difference between the true solution and the approximate solution  at each time point, then take the mean across time.\n",
    "\n",
    "Make an array of time steps. For each time step, plot each solution in a separate subfigure, and keep track of the error in an array. Then plot the errors in a second figure, with each error computation in a separate subfigure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'np' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-3-5380e110293d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# range from smallest timestep that doesn't take way too long to run to largest time step with qualitatively correct answer.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mdts\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1e-4\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1e-1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0merr1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mdts\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0merr2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mdts\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mparams\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2.5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'np' is not defined"
     ]
    }
   ],
   "source": [
    "# range from smallest timestep that doesn't take way too long to run to largest time step with qualitatively correct answer.\n",
    "dts = np.linspace(1e-4,1e-1,10)\n",
    "err1 = np.zeros( (dts.size,) )\n",
    "err2 = np.zeros( (dts.size,) )\n",
    "params = (2,2.5)\n",
    "y0 = np.array((1,1))\n",
    "\n",
    "for i in range(dts.size):\n",
    "    t = np.arange(0,4,dts[i])\n",
    "    sol = scipy.integrate.solve_ivp(fun=degradation_model, t_span=(0,150), y0=y0, method='LSODA', rtol=1e-8, atol=1e-8, t_eval=t, args=[params] )\n",
    "    (t,y) = forwardEuler( degradation_model, t, y0, params)\n",
    "    plt.subplot( dts.size, 1, i+1 )\n",
    "    plt.plot( t, sol.y.T )\n",
    "    plt.plot( t, y.T, '--')\n",
    "    err1[i] = np.mean(scipy.linalg.norm( sol.y-y, 2, axis=0) )\n",
    "    err2[i] = np.mean(scipy.linalg.norm( sol.y-y, np.inf, axis=0) )\n",
    "    \n",
    "plt.figure()\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot( dts, err1, '.-' );\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot( dts, err2, '.-' );\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus Content\n",
    "\n",
    "Here, we use Forward Euler to solve the Goodwin oscillator (coded up in Lecture 8). This code is hear as an example to follow when working on Project 3. Aside from being more interesting than the degradation model, the most important difference is that we don't have an analytical solution. So, we use solve_ivp's LSODA method with really tight error controls (i.e. we allow only tiny little errors) to generate the reference \"true\" solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here is the model, copied from a notebook from earlier in the course\n",
    "def goodwin_oscillator( t, statevars, params ):\n",
    "    X = statevars[0];\n",
    "    Y = statevars[1];\n",
    "    Z = statevars[2];\n",
    "    alpha_1 = params[0]\n",
    "    K = params[1]\n",
    "    n = params[2]\n",
    "    gamma_1 = params[3]\n",
    "    alpha_2 = params[4]\n",
    "    gamma_2 = params[5]\n",
    "    alpha_3 = params[6]\n",
    "    gamma_3 = params[7]\n",
    "        \n",
    "    dydt = np.zeros( (3,) )\n",
    "    dydt[0] = alpha_1*K**n/(K**n+Z**n) - gamma_1*X\n",
    "    dydt[1] = alpha_2*X - gamma_2*Y;\n",
    "    dydt[2] = alpha_3*Y - gamma_3*Z;\n",
    "\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAACW0klEQVR4nOydd3hUVdrAf3fSe++FFHovSehNkKICYseKva51XXvdddVVd+29IAqiCCjSkd4hQIAU0nvvvUw53x+XxITMTGaSmQT45vc8PMrcc899Sea+55y3SkIILFiwYMHCxY+irwWwYMGCBQumwaLQLViwYOESwaLQLViwYOESwaLQLViwYOESwaLQLViwYOESwbqvHuzt7S3CwsL66vEWLFiwcFFy/PjxMiGEj7ZrfabQw8LCiI2N7avHW7BgwcJFiSRJ2bquWUwuFixYsHCJYFHoFixYsHCJYFHoFixYsHCJ0KVClyTpW0mSSiRJitczZoYkSXGSJCVIkrTHtCJasGDBggVDMGSHvgyYp+uiJEnuwKfAQiHEMOB6k0hmwYIFCxaMokuFLoTYC1ToGXIzsFYIkXNufImJZLNgwYIFC0ZgChv6QMBDkqTdkiQdlyTpdhPMacGCBQsWjMQUcejWwDhgFuAAHJIk6bAQIuX8gZIk3QfcBxAaGmqCR1uwYMHChU1TXR1ndm3DOziU8DFRZn2WKXboecBWIUS9EKIM2AuM0jZQCPGlECJKCBHl46M10cmCBQsW2hAaTdv/nz24l/K8nD6UpnvY2NsRt3UjhWnJZn+WKRT678AUSZKsJUlyBMYDSSaY14IFCz2krvIv99eJzX9QXVLUh9IYR2NtDT888yh5SfEIITj060/8+fWnXCxNeZTNTQghsLK24c7/fkb0ousAaG5oMNszDQlb/Ak4BAySJClPkqS7JUl6QJKkBwCEEEnAFuA0cBT4WgihM8TRgoWLjYqC/L4WoVsUpaXw3RP3U5KVgdBoOPbHGg7+sqKvxTKYxtparO3ssLGzR5Ikrn3uNRY/+wqSJPW1aAaxb+X3/Paf1xFCYG1ri42tHWW52Xzz2L0UZ6ab5ZlSX612UVFRwlLLxcKFzqntm/jzm8948IsfcHRz72txjKKiII/YDeuYevNSHJxdqCjIx83XFytrm74WzWCEEJ0UuEatpqasFHc//z6SyjDitm6ktryUqTcvbftMrVJRkplOwIBB3Z5XkqTjQgitxniLQrdgQQ8NNdXEbd1I1ILFtDQ0sPmT95h68534Rw7oa9G6jUqppLm+Did3j74WRSd5ZxPw7ReOrYNjp2vr//tvqgoLuO0/H100u3VTok+hW1L/LVjQQutGx9HVjUnX34ytvQNWtrY01ddTVVTQx9J1zek/t1CpRU4hBCuff4LtX33cB1IZhkatZv27b/Dn159qvT585uVMuPYmuIBt6QUpZ9Fo1L3+XItCt2BWmurr+O2df5J9Oq6vRTGKU9s3s/bNV2hp/MuB5eDswq1vvs/gydP7ULKuqa+q5M9vPiXl0P5O1yRJImrhtYyee1UfSGYYkkLB1f94magF12i9HjEmmoETpiApLkz1VVmYz08v/Z24rZt6/dl9Vg/dwv8PmuvraaipprGuBoCTW/5A2dxMzDmP/4WK0KixsbPHxt6hw+eSJCGEIDfhNMFDh6NQWPWRhLpxcvfgvk++w8pGu6186NSZvSyRcUiSRODAwXrHNDc0kHxwL/1jJuLo6tZLkhmGk4cnC598nsBBQ3r92RfmEmfhksHN148lr7/D4EnTEEKQn5xEQUpSh/jiC5Ex8xaw4MnntNpos+KOs/qfL5B+/GgfSGYYzp5eOLi46rzeUFPNic3rUauUvShV1wghiN2wjvL8XL3jaspK2P7Vx2ScONZLkhmOrb0DA8ZP6hMfhUWhWzALapWSI+t+obGutk0pSpLEvIeeYNHfX7xgj8sgmyz0BQv0GzWGKx99mvDR5s366w4N1VVs/vg9ynKy9I4rSkth17IvyT+b2DuCGUhtWSl7fviGvET9kc8+oWHc/s7HDJs+q5ckMwyVUsmZXdtoqK7qk+dfuG+VhYua/LNJ7F+1nILkjjlm1jY2SJJEQ031BRvf/dPLT7Pti490XlcorBg8eTrWOkwafUlFQR4ZccdRq1R6x4WNGsvS9z4jdLjWpO4+w9XHlwe/WsHgydO6HOsTGnbBRbkUJCex7fMPeyUrVBuWsMWLjLRjhwkcNOSCsxtqo6qoEFdf3052ZiEE3z35IK7ePlz3wj/7SDrtaDRq4ndux9Xbh7DR43SOE0IQv2s7tg4ODJo4tRcl7BqNRo0kKS44ZWdq1Col+1f9gF9EfwZP6noB6A2ERkNFQR6uPr7Y2Nmb5RmWsMWLnPbhT4fX/syG99/uQ2kMx90/QKvTUJIkZt5+D9NvvasPpNKPQmHFyNnz9CpzkP8NZ3ZuJWn/hdfPRaGwMkiZN9bWsOPbz8hNON0LUnWNEIKd331BTvwpg8YrrKzJPBnbpXmpN5EUCryCQ82mzLvCotAvcJQtzfz43BOkHj0IwNwHHmXW3Q/2sVT6STlygK2ff0BLU6POMeFjovDpF96LUhlGYVoyzQ31Bo1d/MwrLPr7C2aWyHDqqypZ+eJTBitoG3sHUg4foDxPvwOyt2iqq+XsgT0GF+CSJInb//MRU266MCp2CyE4uHolxRlpfSaDRaFf4DTX1eHo6tYWPufTLxyvoBCERkPW6ZN9LJ12astKKUxN7nKXUlGQx4FfVvRJAoY2hBCsfes1dn3/lUHjHVxcLyizRlNdLZLCCmtbO4PGW9vYcP9n3zN67pVmlswwHFxcefCrFYycPd/gexRWF07YaF1lOUfW/UxRemqfyWCxoV+knN6xhe1ffsytb76PX0T/vhanE9pqcJzP2QN72PzJf7n1rQ/wCQ3rHcH0IDQachPjsXd2xjcswqB7Tmz6neKMNOY/8pSZpbOgi40fvoOThyczbru7r0VB2dIMGoGNvflMLvps6JbEoguYstxsnNw9tMYTD5o4DVt7B3zCLjyzBWDQzrV/9EQe+PJHHJxdekGirpEUCkKHjzTqnpbGRprq69Bo1BdkklFXqFVK1r39OhFjohh7xaI+lWXXsi9x8/Nn7PyFRt1n7+yMnWPnmi99gY2BpyNzYVHoFzDbv/wYZXMTt/+ncwidnaPjBZmCfnzj76QcOcD1L73RZVifta0t1ra2vSRZ1+QmnMbWwdGoE8+Ea28yo0SGIzQavn38fsZeuYgxRqT1W1nbYOvggJVN3/8eKgrzkbphQpl114XhUzq8ZhUu3j59GhtvUegXMJfd9QDKRt2ORZCLMKlamvt8d9WKnZMTLl7eBsdo15SWsOv7r4heeG2X6d7mZs+P32Lr4MgNL//b6HuFRtOnyVLK5iaChw7Hxcv4TmALn3zeDBIZz7XPvdaj+1UtLX26QUiLPYJveIRFoVvQjl94ZJdjss/E0dxQf8Eo9OEzZjN8xmyDx9s5OVGak0l9VUXXg83MoqdfpLnesAiX9hz5bTXxu7Zx1/tf9pmT1NbBkbkPPNbt+4UQqJXKC+rEZAy/vPYcDi6uLHjyuT6T4dY3/9fnDn5DOhZ9K0lSiSRJenNxJUmKliRJJUnSBV11qaGmmu1ffczB1Rd255bTO7ZSkpXR5bh5Dz1+wSTnqFUqo9uD2Tk6cc+HXzMgZpKZpDIcF09vvEP6GX2fZ1Aw4WOiULU0m0Eqw9Cou69IhEbDd088wP5V35tQIuOI27qRtW+92m2F2D9mEv1GjTGxVMbT134UQ86Iy4B5+gZIkmQFvA1sM4FMZsXO0ZH6qipcffzaPuvJy2AO1Colu5Z9SeLenV2ObQ0NvBD6LCbs2cFn997SoY+lMfTlv6EwLZmEPTu6VaxqQPRELlt6f58lkwD8/u6/WPvmK926V1IoGDx5OoGDhppYKsMRyL/77irEsfMXMHKWXjVlVpL27WLr5x/0uS7p0uQihNgrSVJYF8P+BqwBok0hlDloDaM7U5HAhpGZlFQdZerREwyPs0M0K7nib3/vaxHbsLK24f7Pv0etlJVLQV0BH5/8mJTKFEb5jOKRMY/gYf9XJbfkQ/vY99P33PGfj80aLtUVHgGBDJwwBSc3d+qV9Xx+6nNii2KJcI/ggZEPEOIaovW++qpKVv/zBaIWXGOUucaUJB/cx6k/N3PMM5ddebvwcvDi3hH3MsTLsBKoQggaqqv6pMLepoxNHLBPxkayxi//AJODJhs9x6TrbzaDZLo5VnSM5QnLaVI3sSByAQvmLDDKmasNZXMTzQ0NOHt4mkjKjqRWpvLV6a8obSxlRsgMbhlyC9YKWYXWlJdRkpnR53HxPfbiSJIUBCwGPuu5OOYj5fB+vn7+YR5cfzfZtTkEOwfz49kVbCvfg62bS5+vrOdj7+SMk7sHebV53LbpNv7M+RNPe0/Wpa1jycYlFNX/1b3d2cMLv/D+NNXX9aHEEDJ0BLPveYhGdRN3b72b5YnLsbO2Y0fODpZsWkJiufbKfo5u7ngFh/ZpfZopt9xBwQ0hvHfiPZrVzcQWxXLzppvZk2tYav/+Vcv55tF7uyyKZWo+i/uMZ/Y9Q3JoLccDi3jgzwdYkdQ9c2JzQwMNNdUmlrAzmzM3c8+2e0isSKS4oZgX9r/Am0ff7PEJbdlTD7Pnh29MJGVH4kriuGXTLRwoOEC9sp53Y9/lbzv/hlItb7rGX309t739gVmebRRCiC7/AGFAvI5rq4EJ5/5/GXCdnnnuA2KB2NDQUNGbnNizRfzjkbli0ZqFoqqpSgghxM7snWLU96PE33f/XWg0ml6VRx8Hf10pUo4eFCq1Sty26TYxceVEkVyRLIQQ4nTJaRHzY4y46Y+bhFKt7GNJ/0KjVouG2hohhBCvHnxVjFg2QuzM3imEECKnOkdcvvpyMXv17Laf/YXGj4k/iuHLhouvT38thBCisrFS3PjHjSL6x2iRVZ3V5f35yUnixJY/REtTo7lFbeNA/gExfNlw8fyuZ0VTU6NoVjWLR3c8KkYsGyEO5B8wai5lc7N4/5bFYu/KZWaSViazKlNE/xgtbt90u6hvqRdqjVr8Z+Mr4qU754hlWz7o0dyJ+3aJnPhTJpL0L2qba8Xs1bPF/DXzRUl9iRBCiFVJq8TwZcPFvw//2+TP6wogVujSsbouCMMVeiaQde5PHVACXN3VnOPGjeudf/053jrylhj5/UiRVJ7U4fMvTn0hhi8bLjYeWCXykxN7VSZtqNUq8dUjd4ldy78W69PWi+HLhovfUn/rMGZz5mYxfNlwsTxheYfPG2trhVql6k1x2yjPzxPv3nCl2LFlpRi+bLh45+g7Ha7Hl8aLUd+PEq8ffF3nHCpli2hp7D2F2EpBSba468WZ4pHVd3VY2Ivri8XEFRPFnVvuvKAWfCGEUGvUYtG6ReLKtVeKE39uFO/duEBUFReJ+pZ6sXDdQjH317miUWnczzJu2yZRmJZiJollHt3xqBi/Yrwori9u+6w4K1288sx14rLPJ4iiuiKzPr87/C/2f2LEshHiVEnHxeLfh/8thi8bLg6n7BGrXn3GLIuJNvQp9B6bXIQQ4UKIMCFEGPAr8JAQ4reezmtKsjISWZOyhqsirmKwZ8dY5zuH30mYSz+Of/MD+3/+sY8k/AuFwoq7P/yaSTfcwjdnvqG/e38WRC7oMGZuv7lMCZrCJ3GfUN0sH5FzE07z6T03k5/cNw0L7BwdmXrzUjY3H8DJxol7R97b4fow72EsGbyE1SmrSa7oXCu6sa6Wj++8idM7tvSWyG2sPfIjg1MduCX4ug5hh76OvjwR9QTHio6xM7drB3VzQ0OvFWbalr2N9Op0Hh79MIERgxh/zY24eHvjaOPIC+NfIL8un+WJy42ac9Tl8/GPHGAmieFsxVl25u7kzmF34uvo2/a5b78I7nnhf1Q5tvDBie6bLYQQlOflmLTOfk1LDauSVzE3bC4jfTpmET829jH8HP347OjHaNTqbiVFmRpDwhZ/Ag4BgyRJypMk6W5Jkh6QJOkB84vXc1qaGlnz8nMMi7fj9qGdq7LZKGx4POoJto8uwOHaC8OnK0kSR0qPkV6dzj0j7kEhKTpdf2LcE9Qr6/kxSV6E/CL6M37x9bh6G59YYgqc3D0InT2FrZV7WTJ4CW52ne3hD4x6AEcbR76J72zndHB2IXrhtQQMGNQb4rbRrG7ml/qtZC0NJTrq8k7XF/dfTD/Xfnwa9ykaob9t3u7lX/Prv17slfZ6KxJXEOYaxpx+c/ALj2TyDbe0RYjEBMQwPXg6PyT+QIOyoYuZ/kKjUVOcmU5NWYlZZP45+Wfsrey5aXDH7FohBMEuwdwy5BY2Zm4kt6ab1R+F4KeXnyb2jzUmkFZmTcoa6pX13D2ic50YRxtHHhr9EMdVSQTft5DgwcNM9tzu0qVCF0IsEUIECCFshBDBQohvhBCfCyE+1zJ2qRDiV/OI2k0kiYQoNbaj+jHIU7uymBkyE7fgIH5IX9nn4X9H1v3CgV9W8Ef6H7jbuTOn3xyt4wZ6DGRW6CxWJK6gQdmArYMjk2+8DTdf/16WWKaiIJ/1Z39DILhmgPZu7W52btww8Aa2Zm0lt7bzSzvp+psJHNi7jXX35+2noqmCm0fcipV156Ava4U194+8n5TKFA4WHNQ71+i5V3LV48+2heCZi8zqTOJK47hmwDVYKayoLinq5NS/Z8Q9VDVXsS5tncHztjQ28uNzjxO/609Ti0xdSx0bMzYyL3xep8X+p5efZvuXH3Pb0Nuwkqz4PrF78fCSQsFVjz/L+MU3mEJkhBCsT1/PaJ/RnU72rSyIWECgUyBfn/naJM/sKZd8+dyM+iyOeWQze5zuTEqFpOCOYXdQlp7B92/8Xa6Y1kdUFhZQkpfJztydzA2bi42V7hT624feTq2yli1ZsplCaDQUpaf2SqTC+az+5/Ok/rKBcX7jCHHRHp4IcNvQ21BIClYmrex0TQhBVVFhr/Zj3Ji5kQmZ/rim6t7Jzgubh5e9F6vOrtI7l194JP1GjjZ7csn69PVYSVZcFXEVqpYWvnn0Pg6t6SjbaN/RjPUdyw+JP3R5smjF3smZxf94mZGz5ppc5j15e2hUNWpd7MNGjsEvcgA+jj4siFzAb2m/tZkSjSVs5BiTbWqSK5NJq0rrZPJsj42VDQtS+mOzM4uUyhSTPLcnXNIKXQjBpo3LcGyxYXao/vjmqyKuwllyoCQrg+qiwl6SsDPzHnocu6vH0qxu5qoI/XG5Y3zH0N+9P78k/wLIu+QVzz9B2rHDvSFqG0IIBl+3gBP+RVwRfoXesT6OPlwWchkbMjbQom7pcK2+qpJvHrvXoIQqU1CvrGdv3l4GVnhSeF7v0/bYWNlw7cBr2Zu3l7zaPL1zluVkkXnSvGWhd+TsINo/Gh9HHwSCOQ88yoCYiZ3G3TjoRvLr8jladNTguSPGRuPs6WVKcQFZZl8HX0b5dO5hOun6W9oWkSWDl9CsbmZT5qZuPUfV0kLq0YMmadqxOXMz1pI1c8P0L3ADA4ehsqPtPexLLmmFXpqdidgQz9TmIR0ScbRhb23PuOhZ/DojH4eAvrFDt7Ivfx+e9p6dnDDnI0kS1w+8noTyBJLKk/AMCuaqx5+hf/SEXpL0LznSfKsp8WxmenDXFSCvHXAtVc1VnRyNzh6ezHvoCQaM750yAIcLD9OsbmbG83/n8nsf0Tv2+oHXA/BH+h/651z3C9u++NBsprvcmlwyqzOZETIDkMu1Dp8xW2v99ln9ZuFq68ra1LUGz9/cUM/ZA3uorSgzlcg0qhrZn7+fy0Iv6+QPOr9cxGDPwQz2HMxvab9161kajZo//vsWZw/u7YnIgPwejvUbq9Uf1J45Sx8icO5kNmRsoEnV1OPn9oRLWqGrvexZP7mAYRNmGDR+Uf+raVQ3sjVza684ts4nftd2fv33SxzKOcCkwEmdvvzauCL8CqwlazZlbkKSJAZNnNrrCTo1pSUcObOTwR6D8HPy63L8hMAJBDgF8Hva752uDZs+q9f8AAfzD+Jo7chon9FdVkr0d/In2j+ajZkb9SrryTfcwpJ/vWe2Il1782VFNS1IbopcU1pCdUmR1rF2VnZcFXEVf2b/SW1LrUHz11dVsvHDd8iKO2EagYHYolgaVY3MDJ3Z6Vrivp18eMd11Jb/tYBc3f9qEssTu2XCsLV34Na33u+xHb2ovojUylSmBhnWAHxh/4VtJ76+5JJW6IeLjlDhpmRK5AyDxo/yGcUgq36c/c/3JB/aZ17htCCEoLaphnJVpcFfJHd7dyYGTmRz5mY0QkNzQwNJ+3aZLVJBG0e3rKP/+iqmBhoms0JSMC9sHocLDneylSqbm8g8GWvSHaI2hBAcKDjAdNUIdnzxCU11XWfZXhlxJdk12TozXgE8AoLMGmm0J3cP4W7hbWUUDq9dxcoXdZetuCLiCpQaJbtzdxs0v4d/ILe++T5Dp11mAmllDhcexlZhy1jfsZ2ueQWFMnLWvA4lE64IvwJrhTUb0jd063m+YREGl2/Wxf78/QBMCZqid1zascN8/eg9DBDBeDt4szlzc4+e21MuWYUuhODU+vVEtPjS392whgWSJDF9yFwKnWtR2vV+GdQRl82h6drBSEhMCjTc7DA/fD7FDcWcLDlJU10tmz5+j/TjhttNe0rLUG92jS1laug0g++ZGzYXlVCxM6ej2aW+qoq1b71Khpnlz67JJr8un8HWYWTHx2Fj33Wnmdn9ZmOjsGFDhn5Fk3LkAClHDphK1DZa1C2cKDnB5MC/arWMmnOlXnPRSO+R+Dv5sy3LsLp5kkKBX0R/rRE/3eVI4RFG+47G3rpznaHAgYOZece9HWqgeNh7MD5gPNuzt3fLdFVXWcGR31ZTXVLcbZkP5B8gwCmASHf9JaztnZ0J6D8IZ3cP5obNZW/eXoNPQ+bgklXolcWFOB8tYawYYNTxd074HHaPKSXZtfd2uO05WXySQZ6DcLd3N/iey0Ivw97Kns2Zm3Hz9eO2tz9k1OWGN9rtKfGadEoDNQzzMjwOd6jXUIKcg9iavbXD526+ftzwypsMnWq6HaI2jhQeAeDyBXdw/6fLsLLuekfnauvKlKApXSqak1v+4MSm9SaTtZXE8kSa1c2M8xvX9plfeKRen4kkSVze73IOFBwwWNFUFhVw6NefaG4wvjb8+VQ0VZBcmcz4gPFarzfW1Wr9WV4eejl5dXkkV3ZOQuuK5oZ69v/0PYVpxt8L8mbwRMkJov2ju9QdwUOGc+WjT2Pn6MT88Pm0aFo6bVJ6k0tWoVc5NLFiTi7DJxmnGAZ6DCTYOZgd6dt6NfyvprSErx+9h+KEREb7jDbqXicbJyYGTmRP3h6EEPiGRfRaXWa1SkVG7FFGOQ3F1srw5giSJDEnbA5HCo50MLtIkkTI0BFmrxoZVxqHt4M3wS7BRt13WehlFDcUk1ih2+xy1WPPdKvrUVccLz4OyNFNIDsAcxNO01hbo/e+uWFzjTK71JSUcPDXlZTmZPVAWpnWCBttCl1oNHz5wB3s+6lz3PnM0JkoJAXbs7cb/UzPgCAe+nolgycZfmJsT3ZNNhVNFVpNROfTfjEa6T0SXwdfg3/O5uCSVeinSk6hshaMCdHaHFsnkiQxO2QWISvy2PNL7xX8V6tVOAT5UmXVyFi/rr9I5zMzZCZF9UUkVyZTX1XJgV9+NEnoVlcU5WUwYHcTI+oCjb53bj/Z7HL+C1BXUc7xjb+bdUGNK4ljtPco1vz7ZaP8JdODp6OQFOzK2aVzjJO7h1nKqJ4oOUGYaxheDnJYYV15Ob+8/jypR/UnPLWaXf7MNixhKHjoMB759meTZD6eKjmFvZU9Q70611pXq9VMvfkOIsfGdLrmae/JOL9xBsvcHkmh0NpY3VBOlpwE/lo49fHjs4/z59efyM+VJKaHTOdAwQGa1X2Ty3LJKvSkdRsZXOZFmGuY0fdeHj6H2EGV1EY4mF4wHXj4B8LCYVS4tRj0RTqfqcFTkZDYlbsLodFwZO0vFGekmkHSjmRJxWyYVMjoGOP7KA7xGoK3g3enyICaslJ2L/+KwtSzphKzA2WNZeTV5THSbRjKpiajSt562Hsw2mc0u3J1K3SNRs3B1StIPaJf0RqDRmg4WXKyg7nFwdWV6178F+Gj9W9aJElievB0DhUe6hT7rw0raxvsHB17LDPA6bLTDPUaio2is0nL2saGsVcsImiw9sYas0Nnk1GdQWZ1ptHPLUpPZctn79PSaHjpg1ZOlpzE3c6dcLfwLsdGRsUQMOCvLNIZITNoVDVyrOiY0c81BZekQtdo1IizRQxsCTQo9O98hnsPp3ygLSds0s0gnXaEEJwsOUmgUyD+TsaH7Xk7eDPCZwR7cvfg7OnFI8t+Nmmkgi7iKk9T4aEiKrxzYktXKCQFU4OmcrDgIErNX52C/CL6c99ny4gcp93u2lPiSuIAGBsaw5J/vsPQqZ3D6fRxWehlpFSm6EwyUiisSNq326SF0lIrU6ltqe2g0G3s7Ok3YjQuXt5d3j8teBqNqkZiiw1LesqJP8XWzz/sUfhui7qFpPIknfkUjXW1NNbW6PRHTA+RcxpaI06MobG2hozjR7vlGD1ZcpLRvqMN8r1Nuv6WDk2hxweMx8Haoc/MLpekQq9T1fPz9ByC5xkWRnc+CknBpICJnE46THFW7yj1n176O8rN8Yz2Hd3tOWaGzCShPIGShhJs7XvndJF66hhjmyJwtnXu1v3TgqdRp6xrU7IAVtbWuHh2raS6S1xJHLYKW4Z4dq9uzMwQeQHYk6e7+cXS/37GjNvv6db82jhTdgagQ6ZlSVYGeWcTDLo/2j8aOys79uUZZl6qLikmKy6W+h6UYUiuSEapUepU6HFbNvDpPTejUmo/NQQ5BxHuFt4thR42cgwPfrUCn35d77LbU9lUSVZNlkGnZI1a3WnBs7OyY1LgJHbn7u6TulCXpEI/U3oGgWC03+huzzEleCrRh+zY+kOnGmRmwW/4UHKcqxjhPaLbc7Rmae7J20N5Xg5/vP82lUUFphKxE0II7I8WMTih+53iJwZOxFph3UnRFGek8efXn5qlrk6rGSB+2xZ+feMlo1+8UNdQQl1C9RbrMmXYH0BSeRIuNi4d6uTE/rGWTR+9a9D9DtYORPtHG5z4MmzGLO777PsetXM7XXYaQOd3OmJsNLPveQgbW90ho1OCprQlJhmDpFB0K7mrNcfAkPcw5cgBPrjtGioKOp7UpgdPp7ihuFsROj3lklToJ//4nXFn3RnuPbzbc0wKnMSBkRXUTuudrEXbKQNJC6nXWdXNEPq798fP0Y9DBYdAkihISaK2zHwJOkX1RewYVUjgDcbbz1txsnFinN+4ToqmtqKcpP27TF5XRwhBamUqgzwHobCyxtrWrlsv/qTASRwrOqbTJt1QU82GD/5DVtzxnooMyIpmiNeQDrJOWXI7C594zuA5pgVPI6c2h6zqrC7HKhRWPc52PV16Gl9HX50mRL+I/oy6XH/tnylBU2jRtHTLJp10YA/r3n7NqAU7qUKu6WPIe+gREMTYKxZ1qn3TmozUVXVOc3BJKvSqkkJ8WlxwsXXp9hwe9h74DRzIwTrTpUDrQq1Stu0MeqLQJUlOSDpSeAQ3/wDu/3QZocP114PpCYkViTTZaRg5sGe27mlB00ivTqeg7q/TRMSYKB7+dhXeoWE9lLIjRfVF1CnrGOgxkNFzruDqp1/s1jyTAifRqGpsi4g4HztHR4ozUntksmhFqVGSUpnSKVLE1dsX//4DDZ5nWrAcxrcv3zCzy8mtG9j88XuGC3oeieWJene6pTlZXca6R/lF4WDtYLCpqD2q5maa6upoaTR8d59YnkiIS4hBusMvPJJpt9zZybzp4+jDAI8BFoVuKk6ObqB+fliP55kSNIWqhDROH9Ed0WAKTm3fTO3bG4i069dtW3QrEwMnUtNS07bTMCcJOScZku1KoOhZdb7WrNjDhX9ViVRYWZkllj61So78MTR7WBcxATFYS9YcKNCeEWplbcPdH3zVwWHWXTKqMmjRtHSw+ataWkjcu5OaUsMT4IKcgwhzDevwc9ZHS0MD9dVV3XKMNqmayKnNYaCH9gVHo1bzwzOPcmy9/mYUtla2xPjHsD9/v9GmsRGXzWHJP98xKmInqTzJYN9KU32dTpkmBkzkZPFJo01FPeWSU+hKjZLM6kz6e/TshQWYHDSZkWmuHFz/swkk041veCRZA9QM8O/+7ryV1gSOgwUHyYk/xfJnHqW+qrLH82ojJzWB8QketFToT2zpikj3SLwdvDspmtQjB/n93TdM6lxKrZQVej/HELlU777uLdZONk6M9h3NwXzz78JaT2/td+jVJcVs/uS/BjtFWxkfMJ5jRcc6RBXpHLv4Bq574Z9dFi7TRmZ1Jhqh0blwCiFY8PizDDIg+WdK0BTy6vLIrsk2Wg5jqG6uJq8ujyFehin0lS88xcYP39F6bVLgJFo0LZwoNv8Jvz2GtKD7VpKkEkmS4nVcv0WSpNOSJJ2RJOmgJEmdCx73Iod3/8HMQx6EWxuXAaiNYV7DODq+kaJ5vl0P7gEu4cHsichliLf2eFxj8LT3ZIjnEA4VHsLO0QkHZxeDCk8ZixCCWIdMKu8cil8P+1BKksSEgAkcKTzSoRlDU30dNWUlNNWbTv60qjT8HP1wELb4hffvUQLK5KDJJFcmU9ao3U9RnJnOD88+1uM+ownlCTjZOBHqGtr2mbt/AHf+73MixhjXNnFCwAQaVY2cKT3TI5m6Iq1K/jfr2lhZWVszYPwkfAwwqU0MlENijanr3sqeH7/l93f/ZdDY1l63Qz0New+jFixmyJQZWq+N9RuLjcJG9mf1IoYsvcuAeXquZwLThRAjgH8CX5pArm5TUFeAlUZiUGDPs9ysFFYMDx/HkRLzJgmcyo4FQbfD6M5nYuBETpWewiUkkOtfegOvYN0dhLpLaWMp5c3lDOk3qseV7UDeOVY0VbTtoEE+Mt/21gc4OHffF3I+aVVp9Pfoj6OrG1c9/gzho8d1fZMOWk1Ful5aBxdX7J2cjUpc0kZSRRKDPQd3yKmwsrbGMzAYe2fjTHTR/tFISG21bLpiy2fvs/XzD416BsimLRuFDaEuoVqv15SVUJKVgUaj1nq9PaEuofg5+hksc3sc3dxx8jDMJNjmEPUy7KQ8ctY8Isd1znIFOaporN9YDhb2rh3dkJ6ie4EKPdcPCiFaz/SHgZ5vjXtAaQhsnVRKpGfPTS4g20rdTtewb8tqk8x3PiqlkmMv/49RaW4mVegqjaoticQctd3PVpwlMs8JrxzTzD0hQC4w1Z2X1lBUGhUZVRkMcDdNZ/vBnoNxt3PXKbOrtw/Xv/QGgQO7b0rTCA2plamdvhvZp+O6VeLZzc6NoV5DDbajO3t44uyhvzmMNtIq04hwi8BaoT1888zO7fzw7GMGfTclSWozFRnaTq+V6AXXMPvuBw0am1ieiJ+jH572XYdqtjQ1Ul9VqdccODFgIqmVqTpPcObA1Db0uwGdBYElSbpPkqRYSZJiS0tLTfxombSqNEJdQrGz6rocqiGM9x9PvyIHEk8Yn9xgCEKjpm5qAA0h9kZVWNTHGN8x2FvZc7DgIPG7tvPJPUtoaTKtcyatKo2hWa7UxfXMnNCKv5O/Vofd4TWrWPf2ayZ5Rk5tDi2aFgZ4DGDfT9+z/JlHezSfQlIQ5RfVZUhdTxbU/Lp8GlWNnWzRp7Zv4uAvK7o154SACZwuPU2Dsuu0+Ck33c7kG28z+hlpVWl6S88Omz6Lq59+yaAqlwAx/jFUNld2OMEZgyEngZTKFIOjzNJjj/D5/bdRka+7XlJXJzhzYDKFLknSTGSF/oyuMUKIL4UQUUKIKB8f0zcBaGlswHt5GiPLuu6aYyiR7pEcna4kc6p5Mi9t7OxJCK3GN8I0JwqQs9XG+o3laNFRPAKDGTx5Bqpm0ybopFelc2S2mgV/0/nrNprxAeOJLY5Fqf7LYWfr4ICdk7NJHKNplefsuu798QoKIWRI9/MUWon2j6agvkBnGYCkA3v45O4lXVZE1EV6lZypfL5ynP+3p7juRcNsw+czIXACKqEyuAwAGKYQW6lrqaOwvpABHrpPQu5+/jrNFdpodfZ35wT300tPs+2Lj/SOUWqUZNVkdVn/vBX/yAFcduf9uPrq1jWDPAfhZufWLdt/dzGJQpckaSTwNbBICFFuijm7Q21dNYVu9fj79jPZnJIkER00nqOFR82SyltTXkpeSabBXyRDifaPJrUyFcd+/sy++0Ec3dxNOn9aVRqRHpEmK+IE8hG1UdXYlmEIMPaKRVzxyFMmaemWWpWKQlIQ4RbB0GmXMXPpfT2eM8ZfVkq6dumeAUEMnjwdtbLrqBJttDoXz/9+2NjaGVTDRRujfUZjq7A1yOyi0aj59vH7OLDqB4Pnb3OI6gkNzYw7blTIpb+TP/1c+3VLOYaPidJZAKyVnJocVBqVweGsHgFBjJm3QG+Wq6EnOFPSY4UuSVIosBa4TQhhfBNAE1KsqGTvmDIGjTJtUado/2j6H1WzdfXXJp0XYMeqb7hyl3eP46LPJ8pPrsDXWkO7uztEbWiEhrK8HAacUXToBdlTovyjUEgKrbswU/gB0ir/MseZyq8Q6R6Jp72nTkXjF9Gf2fc81Cmb0FDSq9Lxc/TrkOjS0tjAkXW/dLs8sr21PWP8xhik0BUKKwbETMLPiBNkVwpdrVKy7q3XOLPLuFrnMf4xxBbHotIY52SecM2NjJg5R+8YQxah9lQU5Bn0TkX7R5Nfl98hac6cGBK2+BNwCBgkSVKeJEl3S5L0gCRJD5wb8jLgBXwqSVKcJEmGn+NMTKt9Td9RrztMCJiAfYsV2cWmsRe3x2pkMEeGVZh8hz7MexgO1g4cKzrG3pXL+Ppvdxt1bNZHfl0+DlUaFLH5qExYa8XNzo3BnoM77Wh+f/cNfn+v5w0j0qrS6O/en6a6Wj64/VrijVQo2pAkiWj/aI4W6T/BdTf0slXm9lQVF7F/1XLK83O6NSfIvqHUylQqm7rOUZh681IGTtDfW7M9aVVpOFg7EOisvUa+pFBw87/eZfiM2QbPCXKAQr2yXm9PV12olEq9pXTTq9JRSAqDSuYCrHv7NXZ881mX46L95bBSY8xbPcGQKJclQogAIYSNECJYCPGNEOJzIcTn567fI4TwEEKMPvfHuI4SJiTtm7VcdsK3QwEjUxDsEkzyNBuSR3Tv2KyPIrcGMgMbTK7QbRQ2jPEdw7HiY0SMjWbyjbf1OHyulfSqdLIDGrjsf6/h7hdgkjlbifKL4nTp6Q4NAoKHDCNkaM/s3a2Zi/09+iOEYMy8BXiFaA+pM5YY/xhKGkrIqdWuYPf8+C3f/O0eo08Fao2azOrO5jjfsAge/f5Xo2PQ2xPlL7+mhia+KJuaUDY1GTQ2rVJehHSVrlYorPDvPxA3PfZnbbSat4y1o7c0NfLxnTdwcovuXrBpVWkEOwdr7XuqjZlL72PMvAVdjuvv3h93O/deM7tcUpmipb4qNMGuOkOluoskScQExHC06KjRYVP6UCmVZCWfIcQuECcbJ5PN20qrHd0pLJCx8xfqtfcZQ9vx1HNAt7II9RHtH02LpoXTpX/Z0cddeTXjrry6R/NmVGegERoGuA/A0dWN6bfeRUD/QT2UVqZ1F6bL7BI+OooJ1y5BrTbuhJRXl0ezulmrGcDG3h5r2+5XuRzuNRx7K3uOFXetaGrKSvjozhtIOrDboLm1nSraU5abTfrxI0ZvMDztPRnoMZAjRcYpdFt7ByZdfwvBQ3XXlUmvSjdqUxUxJrpLuzz0vh39klLoscHFOE/oefq8NqJ9o5m4x54/ln1gsjkrC/Lw+jWTYdXmqejYakePLY6lpamRslzTpE6nV6UzJTOI7AOGxTIbw1i/sUhInY6oGrW6R02L22cuKluaTergDnMNw8fBh2OF2l/a0OEjGXflIqMTsHQ5RJMO7OHE5j+6J+w5bKxsGOU7itiirk0BLl4+TLx2Cf6RXRcCq2iqoLypXK9CP3tgL7+/+0a3HN0x/jHElcQZ1Hmpw32LriNokPY8jxZ1Czk1OQbbzxtqqilKT0VloKM7yj+K/Lp88uvyDZa3u1wyCr2yvpySumKTOxdbiQmMocpZSZFVlcnmdPTyZE9UBf6DzLMItbejb/3sA5PFc6dXpRNc5kBRuul94K62rgzyHNRB0Qgh+OrhO9m3svs9XlMrU7FV2BLqEsruZV/x1SN3mUJcoOMJTtdC0dLUaHRtel0hi+nHDpOw2/hem+cT5RdFSmVKhybd2pAkiYnXLcE3LKLLOVtl1ldLKeqqxdzy7/91q+9qtH80zermDic4QxAaDRUFeTQ3dLajZ9VkoRKGR7hkxR1nxfNPGByl02ZHN2Dx7CmXjEI/sG0tt24LpZ/CPLtdfyd/CiY6Ex+sM2nWaErUFWT61tI/uOc1XLTRZkcvOsa4Kxcx+56He7wzVWvUZFRnYHvHJGbd/ZCJJO1IlF8Up0pPte3CJEkiZvEN9I/qfvRSalUqEe5y5mL42GjGXbHIVOIC8s6xvKlcZ//L3999gw3/e9uoOdOq0ghwCuhkjrvq8WdY8i/DGlvoI9o/GoFoi4TSh0ajpjQnq8uGI22BCXqyce2dnfEL757PaJzfOCQko00YJVkZfPfEA2Sd6uwz0LVw6iJ0xGgW/v0FXH0Mq/HUm3b0S0ahV7oqSQyrYUiI+WqDRftHc7zwOE3daDyrjTNnDuBea2O2UwXIMqdVpWEf6kf46HE9jufOr8tvs+uaIjZcG1H+UTSrm9vargGMmXsVYT2ou9LqqAPoHzW+xzb58+nKjh6z6DqmLrndqDn12XVNUT9nhPcI7KzsDIrAyD4dx/KnH6EwRX/j7rSqNNzs3PB20B0jf2bnNkqyMoyWF/6KhDI2Ht07tB9z7n+UQC2n4bSqNKwkK4MjXJw9PBkQPdHg30GrHb03Il0uGYWe61zF2RFKAl2CzPaMKJ+xzN/qyuafug5XMoT0dVuZkOBp8BepO7SPRy/LyTK63Or5pFal4lduR8Pvx81Wlnecr6y4zze71JSVdOuZ1c3VFDfI5jih0dBQU23yJLFg52D8nfx17sL6jRht1IKk0qjkMtDnLfYNNdX8+fUnFGf2vNetrZUtI31GGmQKCBw4mHkPPYF3iP6kvbSqNCLdInUu9qqWFrZ98SHpx7tfsyfaP7pTJFRXWFnbMOKyOVp71aZXpRPiEoKtlWFO5rzE+E5t57qit+zol4xCz8pPpr+eL5IpGB80kbP96ij1NE34YvEUD7LH2uBoY7psy/Npb0ffuexLdi3rWTHM9Kp0HJqtqMvMx8oEu0RtuNu7M9BjYIcdTXNDPV89fBdndm4zer7WI/UAjwHUVpTx2b23dGsefUiSRLRfNLHFsToXi9LsTEqztZtkzie3NhelRtlJoddVlJN8aD8NJlpMo/2iOVtxlpoW/Ukydo5ODJs+S2/GsRCCtMo0vXkgVjY2PPDFD122ntNHjH8MLZoWTpWcMuq+pro60o8f7RRdk16VbtQpefOn/+PwmlVGPburjGJTcUko9ObGBgb9XM7QtJ51++kKH0cfasZ6EOfc/YSO9iQrco3KwOsO7e3oM26/h4VPGt6DUhtpVWk0D3Dj/k++w97JfD/vKL8o4kri2uq62Ds5M//hJ41KcGmlNVpkgPsAbOzsmbn0PoNCzowl2j+aiqaKtgXkfH5/918GK4I25+J5isY3LIKHv/mpR+an9kT5RyEQnCzW3kqvPQ3VVaQdO6xzwSpuKKZWWatXOUqShJO7B46ubt2WeazfWBSSwmizS9ap4/z2n9c7RHs1q5vJqc0xKmRx0d9fYPziG4x6dqR7JB52HhaFbggVTRUcHF5O4MiuO3X3lBj/GBJyT1Jd0bOU98qyYtTJRUTYmya5RR+tdnRrf3fcfHvmNDY2Xre7RPlH0aRuIqH8LxPR0GmX4RlovEktpTIFZxtn/J38cXBxZez8hXgFmb5GfKsdXVds97yHn2TarXcaNFfrIqTLHGeqk+gI7xHYKGwMsu8mH97P7+/+i9py7ZVSDUmfL0g5S9zWjQaH/GnDxdaFoZ5DjVaO/UaO4cZX3+rwu8+qztLbWUkbvmEReAUb994qJAXj/MYZ5IDuCZeEQs9qzCUltI4hg7ufOWco4zxHc9VWT3asW9ajeU6f3Mv0416EYPqqk+fTPmwq5cgBrZ5+Q2iz6+5tJm7bJlOK2Ilxfufs6O0UjbK5iazTJ2kwsvFya6KLJEnUVVaYtK5Ne4KcgwhwCtCpaIIHDzN4QU2vSifIOaiTOe7Ept858LPhhbK6wt7anpE+Iw1SjgPHT+bmf72Hk7v2euHtq1nqIv34EXZ9/xUKq56pnuiAaE6XnTaqZ6eDiyvBQ4Z3SMgytoZLTVkJqUcP6i0joItWO3phXaHR9xrKJaHQkzPjsG9WmDVapJWY0IkcGFFOWVjPGhg39HNg/eRChvQ3zdFZH0O9hrbZ0Q+tXsnJLd1LSsmtzUWpVuKgskbVbFgaeHfxtPekv3v/DoqmuqSYNW+8RGac4bscIURblyKAvSu+48fnHje1uMBfdV1ii2K1ZhSrlEqSD+03yI6uK9uyLC/HJA7R9kT5RZFUkURdi/56M07uHgQMGISVtfZM7NSqVHwcfPTW9Z9y423c/9myHjcAj/GPQaVRcbKka1NRe8pysjr4T9Kr0rGWrOnnaliF1pwzp1j/3r9pqNYfu6+N9ol+5uKSUOhFG/ZzxdFAvBx61n3eELwdvGGYPyc0yT2aJ7Mhh0o3JZFepi0kpg0bhQ1j/cZyrOgYi599hYVPvdCtedKr0kGCaU8+StSCa0wsZWfG+Y3jZMnJtobGXkEhXPvCP+kfPdHgOUobS6lurm5TjiNnzWXaLaZLKjqfaP9oKpsrddrRN3/8Lkn7d+udQ19t7jn3/Y1rnn3VBJL+RZR/FBqhMUg5FmekkbBnh9ZrXaX8g1yYyxSlnMf4jsFKsjI6WSf12CG2fflRW4JRWlUa/Vz7YWNlmIN/4ITJ3Prm+wbHoLdngMcAXG1dLQq9KzIGqqiaYH5l3kq0bxQ5SWcoK+xe+VKA/KMnGNYcYnAxoJ4S7RdNenU6SicrnTusrmg9nka4dZ0xaAqi/aNpVDW2VdeTFArCRo4xqgZ7qxlgoIecth48ZDiDJhrvWDUUffHo1jY23Pb2h0y+8Va9cxhbm7unjPIZhbXC2qC6Lkn7d/Hn15+iOa8ujUZoyKjK0JshCnBw9UpyE4zL8tSGk40Tw7yHGe0YHXX5FTzw+fK275CxPiFbB0f8Ivp3K8tVISkY6zfWrBmjF71C1wgNp2yzCBjZ8+4zhjLOazQzD7qzZ8vP3Z7DZV8Rgwq633HeWNoUTf4RDq35iYwTxnvb06vSGVcezB9vvN6juiqG0mZHb/cC1JaXcXLrBoNtmKlVcuZipHskapWK4ow0k7fja0+QcxBBzkE6X1qv4NAu267pquFSU1bCurdfo6CL5B5jcbB2YIT3CI4XdW3Kil54Hfd9tqyTQsuvzadJ3aQ3Q1SlVHL099XkJyf1WGY4F6BQlmBQK71WHF3dcHKXe6Q2qZrIrc01auE8e3Bvj+SP8osipzaHkgbDm3sYw0Wv0LOKU3EqVRPhFNZrzxzfbzLboksoH9q9tnRKtZJ10wvwmDXWxJLpZqjXUBytHTleepzTf24hN/FM1zedR1pVGr5OfkiShK29eVrytcfbwZtwt/AOR9Ty3Gx2fvs5RemG1aZPrUzFy94LT3tPqooL+fG5x0k/ZvqiYu1pzQrUZkevr6rk8Nqf9dZ1Sa9KR0LqdBJqrq+ntrzMLE2/o/yiSCjvWjk6uXvg4OzS6fPWhVOfcrS2seGx5WuIumpxz4Q9R7R/NCqh4kSJcU7+1GOHiNu6kczqTASCCHfDT5x7ln/NmZ1bjRW1jdayxebapV/0Cv1U7B6uPORPQFPv7XY97D1wGRDCseq4bt2fXZNNg42SyCDt1d/MgbXCWrajF8dy9wdfMf1W4+zIrXbdgKhR3PDKmyYvm6uLaL9oThSfaOtSEzR0OPd+/C0hwwwLUU2r+ivRxdnDk4VPPU+wgfd2l2j/aKqaq7Q2NFY2NXHg5x8oStNd2CytKo0Ql87mOJ9+4dz+n4/MEkMf5R+FWqgNsqMn7dvF6R0dlZquU8X5SApFj8r+tmeM7xisFdZGm13Sjh7i5JY/jI5wAbjtPx8x5SbjSji0Z7DHYJxtnDu0WTQlhnQs+laSpBJJkuJ1XJckSfpQkqQ0SZJOS5LUe9tOoNRHyZ9RJYwYNKE3H0uM5zhqYs9SmG18F6Mz8QcYkulCmIPpY6H1Ee0fTUZ1BlVq48P2etuu20qUfxQNqgaSyuVjro2tHa4+vgbFYWuEpkMWoJ2jEwNiJmlN/zYl+rrUuPn58/A3qxgyZYbO+w1xLpqa0T6jsZasDY5Hj9/VMdM2rTJNa5hle7JOn2TfymWoWowrfasLB2sHRnqP1Fm2WBez7nqApe99Kke4KKwJdTU8ptzR1Q1nD+1hm4ZgpbBi9YLV/CP6H92eQx+GbLOWAfP0XJ8PDDj35z7ANIVODCS9OQd1uDvuzt3/IXeHsT5jiDrlyqG9xocA5pyOIybJg3AjjnqmoDX9+HD6PjZ9/J5R4X+pVakgoOCj37od9tgdtDkZC9OS2fPjt12aHvJq82hSN7U5RMtysgxOve8Jgc6BBDkHaY3tliQJe2fdGbZttbm1OBd3//ANO741z+vlaOPIcO/hHC3serd7xSNPseSfHas9plaldrkIFaenErdtY7ed8tqI8o8isSKxy5DL9tg6OCIpFKRXpRPmGoaNwrAIl/K8XE5sXk9TXffaCbYS7BKss5tTTzGkBd1eQF/N2EXAciFzGHCXJMm0fcn0UB6fzBAMiyE1JePDJrNuRgFFI4zvApQ9FA4slnBy6GyLNCeDPQfjZOPEyarTFKQkUVdZbvC96VXp2AgrAiMG4ejmYUYpO+LtIDfQbt92rCw3m7htG6ntIlv3fLvugV9WsOGD/5hP2HZE+0frtKMXZ6Sx+ZP/au0zmlmdiVqodSpHE9cU60BMQAzx5fHUttTqHWfr4NjhhKRUK8mqzupSoY9ffAOPfPuzSc11Mf4xaITGaDv6wdUraDxuXIRLXtIZdi370qR9dE2NKZbKIKB9/F7euc86pUNJknQf8i6e0NCep7y3qJoZuE+J1TjTtpwzBDc7N4KC+3crlTetOp1wX/Onz5+PtcKasb5jOVoayx8fGrfLTqtMI8g9hAV3mueoqI8Y/xjWpq6lRd2CrZUtQybPYNi0WV2GjrWGLLa+tFNuuo3GOv3KylTE+MfwW9pvpFamMsizY6u7xrpask6doKa0pFM9HH123Rm33W0+gZGboX95+ktOFJ9gesh0neOEEOxd8R3eIf0YNn3WXw0iughZBEzuexnlMwobhQ1HC48yLXiawfdlJ5zCqs64Xr4jZ89nQMwkHFx6z19nLL3qFBVCfCmEiBJCRPn49DzlPacmh/VTCgifMdUE0hlPtMso2J1ObqrhXcibVc14HakmvKb3drntifGPIasmi9IG7fU4dNEXdt1WxgeMp0nd1NalxtrW1qA44NSq1A52Xa/gUIIHDzOrrK20ZgVqc9j1Gz6KB774QWsHoNbMxTDXMHOL2ImRPiOxs7LrsmenJEnkJZ5pK3KVUik7eAd56O/Ruu3Lj0gzcYSRvbU9o3xGGe0YHfHI7Rwcob9V3vlIkoSjm3uvBQR0B1NIlg+09+4Fn/vM7KRUpVLtrGJYhPnT57UxNjCKgVmOxCXsN/ietKKzDMp2wqvGNJ5+Y2m1Se8/vY0Vzz9hUPhia0W6sET45tF7OyWVmJso/ygUkqKDokk6sIftX32s977kiuQ2JdPS2EDqsUM01Bifst0dApwDCHYO1m5HVyh0OnVTq1K1Zi4WZ6bLDSbSepahrA87KztG+47uYN7Sxc1v/LctUiqlMgVrhTVhbmE6xytbmsmJP0VVsenrmMT4x3C24myXrfTak1EtN9gwZoceu2Ed6ceNWzh6G1Mo9PXA7eeiXSYA1UII81WfaUdy3GHCi517LXPxfMaHTWbVnHxywgyvHJerLGTFnFyi5pm2BZqhDPY8FzbVeBZrWzuDYpozqzPRCA0BoZH0GzW2W1lyPcHV1pWhnkM7OOyqi4soTE3W2Tm+UdVITm1Om7mjNDuL9e++YZY+qLqICYjhePFxrXb0pAN7WPPvlzuVok2vStdpunDx9jFryWKA8f7jSalMoaJJf6vF9gtSamUqEW4Rep2LNrZ23PPh1yaLQW+PMa30WkktOcus4740nDbcSX58wzoyT5q/61BPMCRs8SfgEDBIkqQ8SZLuliTpAUmSHjg3ZBOQAaQBXwHmaTSphZrDSUSnehtch8HUuNi6MMhniEGRAa0kVyZjrbDulRou2rBSWDHObxyxVXHc+OpbhA7vumVfazz1uKnzmX33g+YWUSsxATGcLj3dlvgyfvEN3P6fj3RGTKRVpqERmrYdum9EJLf8+38EDTJP/1ZtRPlFUdNSQ3JF5121RqVC2dzcIeO2QdlAXm2eVjOAX3gki595BY8A83XkAtm8BV03YmioqWbd26+ReuwQKZUpbZFEfUGrqciYcrqpdem4qxzRNBu+Gbvvs++Zfrt5/Rg9xZAolyVCiAAhhI0QIlgI8Y0Q4nMhxOfnrgshxMNCiEghxAghRK8tYfvHVtB8tX67nbmJth+O5+Z8MpMM655ScCCWaTmhBre7MgfR/tFk1WRR0lCC0Gi6bMeWVpWGtWRNkGNgL0nYmfEB4ztkBXYVh95q1x3oKSsaG1s7/CMHYOfopO82k9JWH12Lohk2fRY3vfZ2hx13RnUGAtFnvgqQM4qdbJy6NLvYOztTX1VFTV0lxQ3FXSr0M7u2sfnj98yS5WprZcton9FGKfSzFckobx7JyNn6IrI7IkkSNrbGR7X1Jheudb8LqpqqKGouoX+/3qvhoo0xodG41FmRlG9Y5ldTXgkB5b1TkEsXrYpm1741fHLPkg4dXLSRVJ7EMNtIPlt6E0kH9vSGiJ0Y4zumLZqhlf2rlrP18w+1jk+uTMbR2pEgZ3lHm378CHmJWnPjzIa/kz+hLqEGFb0C2oqQDfHqnEG89s1X2PaF9n+rKbFWWBPlF9Wlk1GhsOLWN/+HZogc3NCVQm+oqqKyuNBsDsVo/2iSK5OpaqrqcmxZYxnlTeUM8TQ8Uzsn/jR7fvy2W3XQe5OLVqHHnT3IiDRXwm3MewTtivHhU1g/o5g0z657PJY1lrFzaAGB91zZC5LpZrDnYDzsPDitSmPghCl6Ez2EECRWJDLAayBRVy3uskmwuXCwdmCUzygOF54fJaH9dJFckcxAj4FtCRz7f1pO7MZ1ZpayM9H+0RwvOt5WuqA9e1cuY82br7T9PbE8EVdbVwKdOp+E/CIH9NrPPsY/huyabIrqi7oc23YS6kKhj198Azefl4xkSmIC5KQ5QzJdW7OOg1s8WfbUQ+TEd326Ls3O5NT2zVjZ9N3J2hAuWoWekhjLuBQPwpzM38JNH042TgzzGmbQce9shVwlT9sOrDdRSAomB01mX+0xZt/7MJ6BwTrHFtYXUt1czeCQUUy9eSk+oWG9J+h5xATI0QyVTfLiOeWm25n7wGOdxgkhOsV/3/DKm1x25/29Jmsrk4MmU6usJa4krtM1Zw8v3P3820xeSRVJDPUaqtWcNPmGWxl7Re840icGyvXm9+frj94qSk8l/39riGzwlvsE9CHDvYbjYO1gUPhi63s4vN9YXH18URiQuTruykX87bufTZrlag4uWoWeEdLAtoWNhAb1nTOmlShpMOG/lpCeoD9bLSHhMDOOe+Pf0veJCVOCplDZXCmXH9UTytdqBhjgEG4W+6cxTA2aikB0UjTn+wDyavOoVdZ2UOgOLq64ehvflKCnTAyYiLXCmr35eztdGzt/AbPuehBJklCqlaRWpmpd7DUadZd+DlPS370/AU4B7M3rLHN7nD08qXRREuEeqden0VhXy08v/4MsI0pNGIuNldwM3ZAAhaSKJIKdg/H28OeaZ181ODfhQo4/b+XCl1AHZ0rPMCRwhMma5faE6IFTqHFUklCuP8EoqzAVn3pH3J17rxmHLiYFTkJCYtemlXx27y3UlGqvz5xYnoiVZEXKsjX8/NpzvSxlR4Z6DcXL3quDovn1jZfY8c2nHca1VrIb6T0SgNKcLE5u+aNXarifj7OtM+N8x7Evb5/OMcrmJtKq0lBqlAz17ByFk7B7B5/evYS6CsNLNfQESZKYFjyNw4WHaVHrLqQludizcWQ2A4bozwNpaWhAoVCAmRXipMBJpFenk1ebp3dcUnlSh4VTpVSi0ejOrRAaDRvef/uCD1mEi1ShF1fkE7q/niFK3aaC3mR8xBSOTm7iiEJ34XshBIcdUmhYOsLs1f4MwcPegxE+Izhhl8m0W+7E2k679z6xPJFI90hGX34lo+dc0ctSdkQhKZgWPI0D+Qfa2tL5hkfiHRLWYdzp0tM4WDu0RYvknIlj53dfIDS9t8ttz9TgqaRVpVFQ17kG+qaP3+PnV58loTwB0G6O8wgIZNCkaTi6u5tb1DamBk2lUdWo1yadWJ6IQDDEZaDe05ubrx83vvoWYSPHmEPUNmaGzARgd+5unWPKGsvIq8tjlI8crpubcJpP7rqJ4gzdVVMb62opycqgvqprP1lfc1Eq9JOpBwkqdSDc9sJQ6DYKG6YHT+dA1l6aW7R3wymoL6CssYzRPqN7Vzg9TAmawonmRCIun4Gjq1un62qNmlOlpxjlM4rBk6YxeLLu+h69xfTg6dQqazlZLNftnnbzUkbP7ehkPl16muHew7E614h47BWLuP/z5XqrHJqT1hoj2nbpEWOiGDp1JieLT+Bp70moS2efUPCQ4cy+56EeN1Y2hpiAGOys7PSeLOLL4gkstefkyx9RkpXRa7LpItQ1lEi3SHbl7tI5ptWX0arQvUL6MXLWXOwcdX83HF3duOv9Lxg+83KTymsOLkqFnmZbzOrZBUyeeFVfi9LGBKvhzN/owq49a7ReP55/jIX7AvDP690sS33MCJ6BQLAz809y4k91yrpMrkymTlnHKMch1FWU96odVxcTAidgo7Bhd97uts80GnWbH6BZ3czZyrNt5haQTQg9qWHdU8JcwwhxCWFPXueQz8GTpzP2ikUcLznBOL9xWk2IPS3X2h0crB2I8Y9hV+4unb/3E8UncAj0Zez8hXrj+9e/9292LvvCXKJ2YGboTI4XH9dZBiCuJA5bhS1DvWTTlqOrGzOX3odnYN9Gy5mKi1KhHy06KidA2PXNjksbM8dcRWpYA6dVnbvUAJzOPUGTgyDEM6x3BdPDYM/BhLqEcmjfH6z+5wud6oS0plI7xFfwxUNLUTX3fdlQJxsnJgRMYHv29raU+pUvPNUWo32y5CQqjYqxfnKfFZVSyc5lX1CUrv330htIksTs0NkcKjikNU46tyKbmuLith6q7WmoqeaTu28ibtumXpC0I5f3u5z8unziyzrH76s1ao4XH2d0WDTTb70Ld3/dFbNdfXxw9ugdv9HMkJmohVqnQ/dk6UmGew/vkNgnhKCyqAC1SnvW6L6Vy9j9wzdmkdfUXHQKvbKmjMDfColp6rtsOm24OrhhP3s4W+v2oz7PwSKE4ED1MRoXDSBybHQfSdgZSZKYGzaXvdbxzHrscfzCOhYqii2KJdg5mLFT5zLvwcexse/bhKhWroy4kqL6orYFZ8y8BQybNguAgwUH25JjAOrKy4jf9SdVenp49gZXRFyBSqjYlr2t07WNH7zN5cd8GevbudmXpFAw7ZY7CTZD27mumNVvFjYKGzZldl5MzlacpVZZS4x/DEIIynKzUSm1K8QZt99LzKLrzC0uAMO9h+Pn6KdV5urmahLKEtr6eraSdeoE3z52H/lntfvAWpoaL/iEolYuOoV+NH0/Gkkw1LdvM0S1saj/IppLKtgVv7nD59k12eTX5hlVr7m3mBc+j2ZrNQkexR0UdrO6mUOFh5gQOAHv0DCGTZ/Vh1J2ZGbITBysHdiYsRGQ0+gHjJ8EwKGCQ4zyGdVWMtfdP4BHvl1F/5hJfSYvyKVlI9witCqanEGQOFrJAPfO9X0cnF2IXngt3n0Q/+9q68rUoKlszdraaZOyN28vEhIxATFknTrB939/mNyEztnS+qJHzIFCUrCo/yIOFhzslBh1sOAgaqHu9B4GDRrCzKX34x2iPadl1l0PMue+v5lNZlNy0Sn0/uEj8L/3CmZMvrqvRenEZJ+JLDgYyME1Kzt8vit7F9fvDMI7XrvDtC8Z6DGQ4V7D+TV+FSe3bKCiQK58fKjgEI2qRqa5TyA34bTO3Vdf4GjjyKzQWWzN2trWeqyhpprjx3dwtuIsM4JndBivsLLC2qZvCri1IkkSV0VcxfHi422NN0BeOHeI4wyLnoa1VeeklbLcbJP14OwOC/svpLSxtJOjcXvOdsb4jsHbwZuQoSO4/N5H8IvofGret/J7vn3igV5V7FdHXo1GaPgjvWMTl125u3C3c2e4V8fNoK2DI2PnL8DRzb3XZDQXF51Cj3SP5MlxT2JndeEVyXF0dMZq0Ug2ByWTWys3cRJCsCl1PbWRjvTvP7KLGfqGm4fcTF5lDjuXfUHmSTnjdVPGJlxsXXDNbuGX15/vtRhoQ7l1yK3UKetYnbIagB1ff8ruTz5BISTmh89vG7fxw3dIPmR4vXpzct3A67C3smd54vK2z3bn7qZB1cA0lxjOnlcnR2g0/PzKM+z87vNelvQvZgTPIMg5qIPMqZWppFamcnk/OerD2taWkbPnaY2U8o8cwICYib0aoRPiGsJ4//GsSl5Fs1r2+9S21LIzZydz+s1pi35qj0qpJP34UapLijt8fuS31fz00tM6yzRfaFx0Cv1C57YrH6fFAT6Nk5NdDhUc4mxdGtFLlhA2um8acXTF3LC5uHv6cnKRHaPnLyCvNo/t2dtZ3H8xI2deztX/eAk3X7++FrMDw7yHMT5gPMsSllHdXM3oqxeze0IVEwIn4Ocky9rS1EhFft4FEz/sYe/Bov6L+CPjj7byvssTlhPsHIxDSjUbP3q3w8IphGDew08wcvZ8PbOaFyuFFbcNvY2TJSfb4ruXJSzD3sqeqyL+ijJTKZWkHDnQqdDboIlTmbrkjl6UWObekfdS0lDCqrOrAPg5+Wea1c1cM+AareOb6+v47Z1/cvZgR2eqo5sbHgGBF3zKfytSX4WiRUVFidjYCz/zqjt8sOUNkrduZ9BNi9iWuRWHWsHKu/7o05K5XbElawtP73mahZELyavNI6kiifVXr8ffyb+vRdPJ2YqzLNm4hCi/KFxtXdmWvY3l85czxrdjAosQ4oLIKAaoaKrg6t+uxsvBiwkBE/gx6Udem/Qa8/xm0dxQj4d/35Uo1oVSreSGDTdQ2VTJksFL+DjuY5YOW8pTUU+1jWlpbODz+29nxKy5zLzjXgDqqyqxsbfH1t6hT+R+eMfDHCk8wv0j7+erM18R4x/Dx7N0d7nKS4oncOCQXm/gYiySJB0XQkRpvWZR6KYnKyGOX/7zMhvH5RFZ5s6wZEeWvvcZXsEhXd/cRwgheP/E+xz+YzUhJY7MeOIxwgscUatUjLhszgWjEM9nXeo6Xj/0Oiqh4uH+9zIg2YboBdfi4u2DWqnE2vbCW0QPFx7m8V2PU6+sZ374fN6a+lZbVchWhBAk7NlBv5GjL4jM4oyqDO7dLu96R/uM5ss5X+Jg3VFRl2Zn4hUS2mZe2fLp/8g4GcsDXyzvVZNLK+WN5Tz454MkVSQR6hLKsnnL8HE0vJdxY20Ndk5OfSK7PvQpdIQQXf4B5gHJyF2JntVyPRTYBZwETgNXdDXnuHHjxKVMU32dSKtME2XVxSL58P6+FsdgDu9aL355+2XRWFcr1v3nn+KXf77Q1yJ1SWlDqcirzRPVpSXig9uvE2cP7hMZJ46JT+5eIgrTUvpaPK1UN1eL9Kp0odFo2j5Tq1Ri2xcficPrfhHFmeni3RuuFHHbNvahlB1pVDaKlIoUodao9Y5Tq1RCCCEK01JE/O4/e0M0nSjVSpFSkSKaVc0GjU85ckD89s4/hUatFj+9/LT48fknzCyh8QCxQpeu1nVB/KWsrYB0IAKwBU4BQ88b8yXw4Ln/HwpkdTXvpa7QLwXUapWor67qazGMoqG2RgghRG15mVj79mtC2dLSxxIZx/r3/i0O/LJCCCFEaXbmRSd/UXqq+OKhpaI0O7OvRekWCXt3iuXPPCoa62pFWuxhkXL4QF+L1Al9Cr1Lk4skSROBV4UQc8/9/blzO/s32435AsgQQrx9bvx7Qgi9gb+XssnFgoXuolGrL3gbrj6ULc3s+OYz+o0cw5ALoPaPsQghLlhTXSv6TC6GuG6DgNx2f88Dxp835lVgmyRJfwOcgNk6BLkPuA8gNLRvG1NYsHAhcjErc5B7t8578PG+FqPbSJJ0QSvzrjBV2OISYJkQIhi4AvhBkqROcwshvhRCRAkhonx8DHdOWLBgwYKFrjFEoecD7cMzgs991p67gV8AhBCHAHug713zFixYsPD/CEMU+jFggCRJ4ZIk2QI3AevPG5MDzAKQJGkIskIvNaWgFixYsGBBP10qdCGECngE2AokAb8IIRIkSXpdkqSF54Y9BdwrSdIp4CdgqejK22rBggULFkyKQfmsQohNwKbzPnu53f8nApNNK5oFCxYsWDAGSy0XCxYsWLhEsCh0CxYsWLhEsCh0CxYsWLhEsCh0CxYsWLhEsCh0CxYsWLhEsCh0CxYsWLhEsCh0CxYsWLhEsCh0CxYsWLhEsCh0CxYsWLhEsCh0CxYsWLhEsCh0CxYsWLhEsCh0CxYsWLhEsCh0CxYsWLhEsCh0CxYsWLhEsCh0CxYsWLhEsCh0CxYsWLhEMEihS5I0T5KkZEmS0iRJelbHmBskSUqUJClBkqSVphXTggULFix0RZcdiyRJsgI+AS4H8oBjkiStP9elqHXMAOA5YLIQolKSJF9zCWzBggULFrRjyA49BkgTQmQIIVqAVcCi88bcC3wihKgEEEKUmFZMCxYsWLDQFYYo9CAgt93f88591p6BwEBJkg5IknRYkqR52iaSJOk+SZJiJUmKLS0t7Z7EFixYsGBBK6ZyiloDA4AZwBLgK0mS3M8fJIT4UggRJYSI8vHxMdGjLViwYMECGKbQ84GQdn8PPvdZe/KA9UIIpRAiE0hBVvAWLFiwYKGXMEShHwMGSJIULkmSLXATsP68Mb8h786RJMkb2QSTYToxLViwYMFCV3Sp0IUQKuARYCuQBPwihEiQJOl1SZIWnhu2FSiXJCkR2AU8LYQoN5fQFixYsGChM5IQok8eHBUVJWJjY/vk2RYsWLBwsSJJ0nEhRJS2a5ZMUQsWLFi4RLAodAsWLFi4RLAodAsWLFi4RLAodAsWLFi4RLAodAsWLFi4RLAodAsWLFi4RLAodAsWLFi4RLAodAsWLFi4RLAodAsWLFi4RLAodAsWLFi4RLAodAsWLFi4RLAodAsWLFi4RLAodAsWLPQp1aUNpJ8oQdms7mtRekTe2QrK8mr7VAaLQrfw/xpls5rkI0W0NKn6WhSjUas0xO/Jo7ViavzefPasTL7oFGP6iVK2fBnf9u9orGtBaPqmCqwxqJRqTu/Ka/v7kfUZ7PrhbB9KZFHoFkxARWF92/831SnRqDV9KI1xpBwt4s/vEqkpawQg81Qp+39NRdVy4SvF+D357PkpheoSWfb6qmaqyxqxtr3wX+vmRhVVJQ0ADJsWxA0vRGNrb40Qgs2fn+HPZYl9LGHXnNqRy75fUqgqlv8d8x8Yyaw7hgKgVmv6ZLdu3etPtHBJcXJ7DgfXpnH3u1Oxd7IhdksWOfHlXPdMFLYOF/7Xa+jkQLyCnfEKdAagNKeW/ORK+qhNgFGMnBmMb5gr7n6OAIxfGIEQAkmSaGlSkZ9cSfioC7N378FfU8lOqODWf07AzsEanxCXtmsRo31w9rDvQ+kMY+zcfoQO9Wr7+Tu62uLoagvAsQ2ZxP2Zy5KXx+Pm49BrMhm0lEuSNE+SpGRJktIkSXpWz7hrJUkSkiRpLb5uQTtCI0g+UkRRZnXbZwWpVfRV8xFjGDzBnwmLIrC2kb9KoUM8GTol8IJW5qoWNVu+OENpbi2SQsI/3A1JIQEQsyCC6/4RhY2dFRqNoKle2cfSdqY0p5bGuhYkhURApFuHa5Ik/ztiN2Wx5ct4aiua+kLELhl9eSgxC8KxtrHq8LkkSYyeHUr/cb4AlOfXob7ATnxN9UpaGlVIkoRPqIvWMaNmhTDj5kG9qszBAIUuSZIV8AkwHxgKLJEkaaiWcS7AY8ARUwt5qdPcqGL/L6mkxZYA8gu77r0TJB0s7GPJdFNRWI/QCBxcbBk3LwxrW/nFDB3mxejZoQCUF9RxfEtWH0qpncqiBoqzamisbdF63erc4rT3p2TWvXcC5QVkfhFCsHvFWTZ9ekbvuJirwln42GhcPC/Mna6HvxNDJwfqHVNT3sjqt2KJ3ZjVO0IZyNH1Gax49bDe74WDsy2DJwYAUFvRREON9u+aqTFkhx4DpAkhMoQQLcAqYJGWcf8E3gYumC2B5iJwrADYO9lw3bPjmHxdfwA8A52YedtgBkb79bFk2mmsbWHN27EcWJumd1zy4SJO78zTqTj7Cp9QF27910RCh3rpHRc51pfBEwLaTh8XApIkcdntQxi/KELvOGtbK4IGegBQlFFNTXljb4jXJfkplez9KZnmxq6d0K5eDky9YQAjZwb3gmSGMyDaj5irwrGxtepyrEqpZu07x9m9onecpYaci4OA3HZ/zwPGtx8gSdJYIEQIsVGSpKd1TSRJ0n3AfQChoaHGS9sFQgiyTuQTrtgNRfEcPjuUvKogFv1jCnZOdiZ/Xk+J35uPsknN6MsCcCvbCSePgMIKq8AxDI2ZAzZWNDco2fZNIuMXhuPbz7WvRQbkBWjqjQPxs02DnT+DugV8h8LAeeDg3jZuwtWRjJoVgoOLbd8Jex61KadxLtyGVXMVeEbAoCvBRfvCGTLEk5AhngColZq2nXuvU5ULib9DbSG4BuE1cC4ERRp0q7JFzabPThMQ6c78B0aYWdB2NFRA4m9QkQmOnhB5GQSMojy/joy4Uiad27x0xbCpQYD8bjfXq7B3tjGfzC31kLQBSpPAxhH6TYJ+k+GcGauVgP7uBPR3N2hKaxsrpt44EM8AJzMIrOV5PZ1AkiQF8F9gaVdjhRBfAl+C3CS6p88+n4Tf9rFnq4obvN7Dx7GYgPoRiKZB2C17Hua/jbrfNKysLpzdVlFGNQ1lVYzOuBWKToGVLQgNaFTgEgizX0UZspCaskbqKpvx7dfXEstITZUMznwcUraAZAUKK1mp27nBtKdg4iOgsEKhkHBys0MIQfyefIIHe+Dh3ztf7E6oWhCbn+W3LVH429RyudeXoG6Gzc9AzH0w83mw1S5bRUE9f3wUx2W3DSFkqGfvySwE7P8f7HoDNCpqCOZ0zeWMcXoHpzGXw5w3wFm/09PG1or594/Ao5cUCgDxa+GPx6C5Bqzs5J/zn69CxAxGXvlfhk2dhJW1ce/h7h/PUpRZw/XPRXWyu5uE7EOw5m6oyQeFjfwOIsB/BFzxHoTKe9j4vfmEDvXE1dtw23jE6L9+Ry2NKrP6lwz5qeYDIe3+Hnzus1ZcgOHAbkmSsoAJwPped4xmH2Ro/I3MDV6O970fw/P5hL/2G5PvnAZqJTXf3seqZzdTmlPTq2LpY/ZCB67Q3IdUnQvXfgPPF8ILxXDrGnANgHX34bznCW56bkSHL0Vfkrw/k8T3X0Ok7YQ5/4LncmWZ7/4TQifA9pfhh6uhvqztnqY6JUf/yCRhX0HfCK3RwK93ImK/I2Z0GYPvuBNeLIaHDsPIG+HQx/DlDKjI0Hq7i7c93iEu5t0dauPPV2HHazD4KnjsNIXzdxDfsghN1L2Q8Bt8PkVWRF0Q0N8deycbhEaY3/QS9xP8eif4DoH798k/56fTYc4biLyT8NlkrM7+ZvS0EWN9GTo50Dwbsqz98nfW2h5uXw8vlsjf60WfyieN7+bBgQ+pKWtgz8pkss6UdTmlNmI3ZfHzG0dpMcDc1G2EEHr/IO/iM4BwwBY4BQzTM343ENXVvOPGjROmojIjVzS/OViID8cJUVfWeUBznaj88WnxyyOfi8of/yGESmmyZ3eHlmaVaKypF+LzqUK8GSpEUULnQWqVEDv/LcQrbkJ8d6UQTTUi60yZ2Pz5aaFWqXtd5lY2v/StWP3IZ0Ikb+18UaMR4sQPQvzTV4iPooSoymu7VFVSLzRqTS9K2o5dbwnxiqsQBz/Rfj19lxBv9RPirTAh8k/0pmS6SfhNlnn9o/LP9RwtzSr5fwrPCPHhWCFe9xHi7CaDpty94qz47h/7RFODmb7/+SeEeM1LiGULhFA2d7q84f0jYueL78r/rgMfmUcGY6kuEOI/kfL3ta608/WmWiF+vl2W+Y8nRE1pvWisa+nWowrSqsS+1SlC2aLqkchArNChV7tc7oQQKuARYCuQBPwihEiQJOl1SZIWmmORMQYhBFs/PcSmgofgxh/BSYujy9YJ95vf5rprKnFP/Rx+e5Dm+r5z1J09WMj3zx6gJq8AFn0Mfp2ChmQzxsznYPEXkH0QfriGhopaqssaaW7oo6zGpA3MlR5nwcIGGDin83VJgjG3wq1roaYQvpsv/xdw83FEUkg01SspTKvqPZlLk2HvO2iG3UCyYjHNDVrCECNmwD07wNYZll8NRdojSDRqDUc3ZJKwL1/rdZPRWAl/PA6BY2D+Ox1suG2OOP/hcPd2+bvz861wdlOX0w6ZHED0VeHY2pvBZKFWyWYWRy+44Xuw7uw38QrzwmPKYhi6CLa9AIc/N/oxBWlVrH33uOl2uVufh+ZauOEHcPLufN3OGa5fBpMfg9hvcDn0HPaO3TOZBES6MeW6AeYxGZ3DoPOLEGKTEGKgECJSCPHGuc9eFkKs1zJ2hhAi1tSC6kIqOMFUm/cYO9EKfAfrGSghzXwWZr3MqcP1rHrxT+qr+iYgxz8IxjmvxWXoeBiyQP/gUTfCdd9CfixDcv/BdU+P6hsno1oFf76K5DsEu8v+pn9s2GS443doKIcV10FjVdulncuT2PJlPCplL4UCbnsJbJ0oGvgCf36XRN7ZSu3jvCJh6R/nlPoiKE/vNERSSBSlV1GaY+YMwH3vyUp94UdtivHM7jw2fnq648/N0VM2EQSMks0cOYf1Tuvbz5VhU4OQJMn0OQ6nfoLCUzD/LXDw0DpkwqJIRs+NkM2Lg6+CLc/AyRVGPUahkGiqU1Jf3dxzmbMPQcJamPx4l7qjPup5djl+SuWR7bDzXz1/tpm4cDyE3WX3WwS6FRG6+FbDxk95kqBxQwmXdmJ/5mvzyqYDn8wviHL4CWnWS4bdMOxquOJdSN2K1eYnUbWoOLQujfoqE3ypDSX+V/5IvZF4r5fBygBbctA4uPEHeYe86hZQySeiSdf056pHRpl1l9JG4SlI3QqTHiFgRD+ufWYcocP0hCp6hMEd62XH9MobZaXaDkmSuOKhkcy4Rc/L31PqSuHIlzDqJtkhdw6NRiBJdP652bvCzavBLRhW3iD/vLsgL7mSX9+KNV3SlFoF+/8LAaNh6NVahzTUtPy1iFjZwHXfQcRM+ONRyNxn8KP8I9y46eXxpnGu73sXnHzk3XcXVJY0klYQjGrQNfJ9x7/v+fPNwEWt0ItPxHHspDvN4x6Vv9iGIEl43/A80yaUYrXzJdTJf/ZqRmZefAGVhzbD8Gv17wrOJ/pumPY0nPyRul3fc3p3fredM0YjBKq9H6Owd+qgZLok8jK4+jPI3i/vxgB3P8e27LrGOjObvfb9F+xcIfretoxQG7suFhKvSLhpJVRmwc+3tS1ErbQq1NqKJjJPlZpe5uPfyVEhU57o8PGoy0K44sGR2u9x8pLNXFa28NOSDicibdjaWyGE7Kw2CYm/yQ7laX/vFOIHsln0l38fY/eP7WKxrW1lU4ZnpGwyKtOf09AehUJCrdZwakdu94uqFSdA2p8w/n6wdexyePAgD+5+byreN70M/WfDhicgY3f3nm1GLmqFnr/vAHH1C5HG3WbcjQoFXP0pTZ7jWPNJNmc2xplFPm3s/fE0e8pvhfEPGH/zjOdh0JW4H36aW+/+K0bX7OQcxrr8DFfe5MLwGSFdj2/PyOvP2R+/hePL2j5OOVrEDy8corKoXve9PaGhAs5ugLG3U9/iwJE/MgxPg+83SfZtZO2D7dpPUQfXpLHrx7OmzSJVtcCxr2WF4TOo7WODEuQ8+sENy6EqG9beCxrdcvn2c+X656LaapD0mJM/gns/OaZfC0IjiL4yjP5R58X7O7jDLb+AwhpWLZFt2QZSnlfH/l9TSTte0j2Zj3wux5pH3W3wLQorBZLNuYXIeyD8epecI3ABcfEq9JZ6xja8xe0ztmDr0Y2QPlsn7JZ8iYdNIc6nP4TmOtPLeD5CsDDwfaaG74HgbkR1KhSw+HPwisRpy51QlUNlUb35bbqx36Kx9YDh13Tv/lmvyLv1jX+HHLkyRNAgDwZP8MfJzUwJX4m/y7HEI2+gNLuW45uyjDMxjLoJJjwkv/hnfu10ecr1A7j2H1EGZQsaTOo2qCuW4+LbcWxDJitfPdx1TZN+k2D+2/I8u/6td6gkSaiVGmKN/bmcT10JZO6BEdfJ308tKKwUDJsa1Jak1QGPMLj+OyhPg98fwdCqaL79XLnpxZguywdoRdUsfz+GLpL9EF1QW9HEuvdO/FVryc5FDsBQK+GX20B5wSTHX8QKPXE9NNdgN8FA27kWJO8ILn9oMhHK3+D3hwz+MnWb4nicKw7hNeVKrUdTg7B3hZt+ArUSseo2tnxxmt0rk81nNmqqQZOwnu+LP+XE7m6aeBRWsiPMLVh+AWqLcHKzY9qSQeZLsohfA14DwH8kYSO95eNykLNxc1z+OoROhPV/g5KkDpec3O3aCi+ZrABWwlpwOJdV2Q6vIGdCh3sZFoMddTeMuU228yZ2ilnoQGVxPUc3ZJIR1wPTUcJvss9h+HU6h5Rk1+hfNMKnyYt+4m9w+FODH+117vdZV9lsnOklfSc0VctmTwNorG1B2azu6L/w7i9vrgpOwmadyfG9zkWr0Lf+WsuhlofkZJaeEDkTZr9K6vFydryz1qz29JNrjpDbMhqGLu7ZRN794ZqvkIrimB3xB/PvH95WZc/kJG9GpYKBo12MV4jtcfSEm1bIx+rVS+XdDXIN79/fP9mh0mSPqSmQk0VGXNe2cNo52rRVVDSYVuedrbNs523qnJQWvzefFa8c7rnpqKUBkrfA0IWdnM79x/ky5boBhs0jSXDlexA4Fn57CMpSdQ71Dnbh5lfHd2+X20r8r3LZB22ht8j2842fnGbvqhT980x+TI582faSHKZrIE31Slb96wiHf9OeFKZd5jVyJE7EDIOG+/Zz5Ybno/EOPu/7P/hKmPoUnFh+wThJL0qFLhoqsGvIxMY/ovs73fZMepQ6n2lU5xahTN7b8/m0oFapOZnoQ67DVdpj5Y1l0DyY9jQ+6Z/gnLEKIURbkwaTkrAWW3dvJt8xWX+EiCH4DZND8XIOyS8uYG1nRUNNi2kjdhJ+AwQMv5aq4ga2fh3ffYXrGiCbBCoytZ7iIkb7MHZOKK5ePSyTmroVlPUwrKNZq6leaXyzDWs7OcLI2lZeiPSYE919ZTt6TVmj8U7SqhzIPaJ/pytg7r3DGD27C9+LJMHVn8ommNVLobbYIBHsnWyYsCiSkZcZWMCrpUGO2R+6yLBoLdC/yZv5gnyi2vS0vFvvYy5KhS4lb2KG66dEXTvORBNKjH7kQa4euBzb9XdBtekTR6xKTrPUeylRl/ubbtIZz0HkLNj0NCd+Ociqfx41bWp3YxWk7aAm7CaEqU4AI66D8Q/Ckc/gzK/YOVhz44sxRI7xNc38IO8a/UeC9wBqy5soTK0yunZIB8KmwOxXIekPOPhRh0uOrrbELIjAykbRs9Nd/Fpw8pWf1Y4TW7P59h/7UauMrAnuFiyfLspS4PeH9ZoTW5pU/PLmMQ6s0b2b1y7zGvm/ehS6pJAIHOBhWGE5ezd5IWqqkePq1YYtMMOnBbUtTF22rkvdpnXh1EVNWSPLnjlATkK59gEKK7jmazn88efbZWd8H3JRKvTGk5vklTxwjMnmlOxdUdz0A8oWwYH/fk9zjYmjL+LXorBSYDtSeyRAt1BYwbVfg7M/A7OeIWq2r2k7vSRvolllw4/bJ3N8S7bp5p3zTwidJNumixNQnDOFZMSVkpvUwxeiIgPyj8sLBxAy1JM73ppsVDElrUz6GwxZKNdXydrf6XJ1aQOr34ylJLsbtYKaa2VFM+xq+XfajvBRPkxYFNm9BSli+l+26UOf6Bxma2/N9JsGEX1VuHHzx6+BoCjw1H1fRlypca3Y/IbBwg8h+4D8szYQIQR7fkruukxt/BqtC6cuNGpB0CAPXLz0vFdOXnKEUV1RlxFG5uaiU+iNJSUsO7qUM3YPmcbc0h6fgVRO/oAzxWPJXaX7BTAWodawcasPGe53GeRVNwpHT7hxOS4tKYytfBYFGtN1eIlfi+QWxPSbBxE+UktadHexspFDv+xc5aSjxio0ag1H1mdwakcPw8Bad43tdmAm8S9IEiz6RC65u/rOtpIGrdg5ysd3Q+p8dyJ5M6iatO4aAyLdelYPfPJjcjby9pf1JvAMiPZrMxsZdBooTZFLJIzQ7QwV55pxnNyeY5zMI2+A6HvkomkJvxl0iyRJ2DpYY+too/uk1FRzbuFc3Gnh1IW7nyNz7h7WdSJT8DiY95Yc277nPwbNbQ4uOoUupW0jxvkngiZFm2V+35kLuPWKWPoXvGF0WrIuGlNjqW+2R+k/ySTzdSJwjOwIy9hN1e/v8dOrR8g728OdbkMFZOzCdsR8hk0NaosoMBkufnLNj+pcWHc/CgmuemRUz2t2x6+VI1PcQyjJruHnN46aLqzT3lU2CbTUdXDsgmzLvf65KEIGd2PBjl8rl0sO6dBmgJqyRiqL6ntmypEkuWqgZ4RsxqjRX/Fy94qzbP0qvutnxv8KkkJWjjofLXHjizHEXKW/GYdW5v5b3v3//rBex257JiyKYPK1/XUv4K0Lp4HRLUII7bV/dBF1F4xaAnvehtTtht9nQi46hW4//ibGPXIPniNMZ245H+crn4XwaZSu/YCCQ8d6PJ9j9u/c4PcCAxfMNIF0Ohh7G4xbitPJd3FzrMOqp6n1Zzci1CqybK/q3q7TEEInyC9uyhbY9y4unvZYWStQKdXdM10UJ0JJYtsLq1EL7ByscXQzYe0b3yGyYzf3sLzrbUdrjZTE/QXkJ+uoGXM+jVXyrm7Y4k5x3Kd25vLzG8fQqHoYeWXvKsdNtzTAL7d3yn5tj4e/k7x463ukEHJsftgUcNHvE3Jys+teX01ru3NFvuy6dOy20qrIywvqOLQuvfOiFL8G3EIg2LDNYHVJI18/tY/UWMMctHKE0X9ls9Gae6DShGZKA7noFDpW1nLxJ3OF6Z17hrj2W3bWPMz+nxMQ9TocIoag0SDif4PIWUjtuvmYhfn/wSZ4GAu4iwC3HqalJ6yl2imajavqSTlSZBr5tBFzH4y4QU6EOber2fNTCus/iDNudwTndo1WbfVE/CPcuPrJsaZPXhpxHcTcL8dMx6/tcEml1HBye47h/WDPbgSNUmvS1qhZIcy7b7hpOiX5DoarP4G8Y3KFQR2MmhXC+IUR+kM8C+OgIl1v7DnA6V15ZJ7uQXkKt2C5MF1ZilzzxcCTSnZ8OUkHCzpGTjVUQPoOrQunLmzsrIi5Khy/MCM6hdk6yqc4IeTFs5eTji4+hd5LSM4+zL17CAvcX0fqgaOjIeUo36W/TobjjSaWUAvWdrJzRmGNWHUbJzandq/Ua305ZOzBZfRUFj81hsixJoxAOR9JggUf/LWrqcgk+oowZt85tM0ubRBCyDuwiOng7INGI9CYs1v8nH9BcIyc3diuIJaNrRVXPzmGWXcMMWye+DXgHioXMzsPVy8HwkaY0HcxbLHs3D32FZxapXdoaW4tmz47rb20wZlf5a4+Q3VXzxZCcHpnLpk9SVoCOVb8shfln9ORLwy6ZczsUJa8PL5jgEDSH3LmsBHZzk7udkRfGW68Q90zQk46Kozr9aQji0LXg/vIGByuegmRtoPcVYZnsLVHmbCNYLtEXEdMNLF0OnAPOberSSZ//2EK07qRsHP2DxBqrEYsJnCAB46uZi7X27qrQcAvt+HqRpsiMzg+Pe+YXFDr3K6xKL2Kr5/aR2G6CROW2mNtK5sEbB3hp5vkRfAcTm52SAqJ5gYlGSf1KLS6EsjYJct83omzNKeWjLhS0zm4W5n1KoRNlWuXF57WOaypVklpTm3n3AaNWlau/WfrLJMLsvnjltcmMPl6AxOi9DH5CblOzNbn5SzPLpAUEg4utgghSD5SJC9Kp38Br/5yRUgDEBpBcWZN9zcFg6+AKU/KSUex33Zvjm5gUehdMW4pyb7Ps37vMHK3bjHuXrUKt6yVzIlJxDvShPHnXRE5E2nWi8y1+juzBnb9AnQifg1K98Ecj3M1XVp7V3hGwDVfyZETG54EISjNqeXHlw+RbIjJJ26FXGzp3K7RzsmGAVF+eAaasZema6BcmbE6H1bd3Ol4fXRDJtu/TaChRofNOn6tnDY/svPpLWF/ATuWJWJyw6KVtbzgO3jKZRh0xE2HDPXkltcn4BV4njM8Y5fcrHrUTV0+SlJI2JmitENrDSOfwfDLHZ3KMOiiLK+OP79LJGl7glzxc+RNBptqy/Lq+PXtWFJju1n8C+STRf/L5RpGveQkNUihS5I0T5KkZEmS0iRJelbL9SclSUqUJOm0JEk7JEm6QNoZmwBJYuA9jzInbDXBx+6SY5wNpCXxTxqqGmD0zWYUUAdTnsB22Byk7S/SeGwt+35OMaypREUGZO6lOHAph3/PoKqkwfyytjJwLkx/Fk6thL3v4hXkxIjpwQQP1r0TBEDZCPHr5DhxO7k0r1egMzNvHWwahaKPkBhZ2eQehvUdi0uNXxDBwsfH6D7hnPpJLkespYzy1BsHcM3T41CYo4ems69smqspkE8XLdp/x9Y2VgghiPszh+LMc07quJ/A3h0Gzdf7iF0/JJF00IQ9ZO1d5cqMNo6w4gaDMkl9QlxY9PhoRjhvlT8YeYPBj3PzdeDyu4dqLyhmKAorOcvYb6gcFaXnRGQquvy2SJJkBXwCzAeGAkskSTq/cMNJ5D6iI4Ffgb4LxDQDCjsHBjz8CpKzF/XL7qb8TIJB96VsO8Ky0m+pcp9uZgm1IElwzZcQOpGidV+QuC+P0hwDKkoe/x4kK4LnLeL2f08iaIC72UXtwPRn5B3rrn+hOPEdk67tj5ObHUII3QWeEn+H5moYvQSQiymZvRlye4ZfIyfwnFkN215sU+q2DtYERLoBUJhW1bFtWu4x2cY65natU1pZKUwfKtqekGg5KS33aKcQzPYom9Sc2Z3H2cOFcuPvsxvkKCJr3Y5mtUpDZVGD6RuwuAXDzaugoQxWXt9l3XeA4AGuSCe+pyn4coqq3A1+lK29NQOj/XtubrRzkRuQ2LvDiuvNHvliyPIfA6QJITKEEC3AKmBR+wFCiF1CiNZl/jDQg0yICxQXf8St69ha+hCbvzyDpixL//iyVEIrf2TCiGzc/F16RcRO2DjAkp8ID67lVr9HCLDtoptNcy2c+F7efbkG4uJpb54doj4UCjmBZ8Bc2PhUW+navatSWPfeic5OOo0G9r8PPkMgbBogR1f8+OIh3aYOczDlCTli59DHsOP1Djv1hpoW1n8Qx9ENmX+NP/KZnFh1bhFqz+lduT1PsDKEoYvgqv/KdWR+e1DuPHQetg7WXPP0OKbdOFB2Sqqa5KYQerCyVnDN0+MYNz/M9DIHjpFPF8WJ8OM1ctVEfST+BtW57Ky4ny2fn0Gt7NomXlfZRNLBwu43zzgf1wC4ZTWoGuH7q8yq1A15W4OA9t+uvHOf6eJuYLO2C5Ik3SdJUqwkSbGlpWbo9mJmJO/+XHZPFJd7f47i+3lQoifNePdbuNrXMvaOhearhGgIDu5w6xo5dG/51WRt28mJrTq+UEc+h8ZKTts/zI5licbXDzEVrZmk/SbJkS/HlxEx2ocB0X5Ynx/Cd/YPKE2SFeq5cLShUwKZtXSo+Z257ZEkmPc2jFsqt2P781V5sUGu93L53cOIvjJMHlt0Rrafj1vaZiJqT3Z8hfmcuecTdZdcp+bMap1hdk5udkiNFTQd+pH9Vq+jcu+vd8pW057ZvvcDLpeVeuEp+PFa3fVT1ErY+w549WfK0knMucewEND0k6Xs/CHJdB2dQDa73P67vAB9f5XswDcDJt1+SZJ0KxAFvKPtuhDiSyFElBAiysenG00pLgDch4/F74FPQGg4++Hr5O/8s9MYTfKf7D/oQtWIp8H5Avh3ugbAnVvAI4z0LbtJ35+IWnXeTrc0Bfa+C4OvosU2gMY6Zc8KWvUUW0e45Vc5muKPxwjJ/5CoOcFIkkRtRRPKZrX8Im95DnyHdQhHc/G0Z9D4XnRCt6JQwJX/g3F3woH3Yc1d0CLXBIoY7YOdow2a5mbOLvsWYecOU5/UOs2Cv41i9p0Ghj2agilPwPx3IHkj/LC4U1kDhICtz5FfF0588QjKcnWb7oqzavjuHwcMT6zqLoOvgOu/l+3SX12mvZfqgQ+g9Cxc/jquvk4EnjMf5iSW6zUHjZwZzE0vxvS8/s/5BI45p9Rr4PDnpp37HIa8sflA+9qXwec+64AkSbOBF4CFQohe7F7cB/gNQ33HZuLqriLu91j443H5GCUEpG6nYtXLJDbNoSygF2LPDcXFD5ZuYMaoRBZId2P1yy20ZJ+ipVEJBXGw4lrZ4XTFu0RdEc6VD+voX9mb2DrKUSRj74D9/4Nv56BK2c1v751g2+fHZJtkfamcuWllg7JZza4VZynP74XuU7pQKOCq/8nNMRJ+g08nyP9VK6GhgvQv/82OrAXkjPywU9ifqkXdlkzVK0202zP+PrkJSWEcfD5ZboXXXCc7TLe9CKd/JnLOVG7712T8I2S/gLbKhnaO1oSN8MI7tBfMjEOugqUb5FIMX0yTa6g0VMg/6yNfws5/yfVxBl3RdktTvZJtXydwYI32HqYatQZJksznvwgcA/ftkgvUmQGpq5oNkiRZAynALGRFfgy4WQiR0G7MGGRn6DwhhEGFF6KiokRsbGx35b4gaK6qQdr7NrZxn1PZ4svxhpuY5vwZtr6h1F/9C07BYX0tYmc0atnOu+cddpXeRqFyCDd5PYZw8uVk6HcERQ9rc+RdUJz5Fba+AHVFZCknYC/K8XcqQL3oCxh0JVY2CsoL6lj7zgnm3z+c4O7UVDE12YfkipLlqXIijlAjNILC0f8lcPFdgNyo2d5ZTqDa/0sqyUeLWPLy+N41F7WnNFmOUc85JGfdSpKckBN1N1zxbptZKyehnP2/pnHFAyPaepNq1Jre97kA1BbBpr/LyUMgN8tWt8ghgzcs79QEuiyvFhdPe+wcbWhuUGJrb42kkKgorOePj+KYdpOJi9GZGEmSjgshtPaw7FKhn5vgCuB9wAr4VgjxhiRJrwOxQoj1kiT9CYwAWs9qOUII3WlkXBoKvY3KbBLW7yc2zoOrF1XjNnkx2JiwjK05qC+nYOdWCjLqiJogECOuZ8VbyQQNcGfmbb143DeGlgZI3iTvIh08YcR1HD+s4cTWHO79n+wQbahp6TtlqA21CtK2Q85hsLaXS+T6yj/fusomfnrtCBOujmTEjGAqCurJSSxn9OzQvpVZCLlrUMZuWZkPmAP9OibGFaRWErs5m/kPjMDG1ood3yfSVK9i/gMj2soh9zpFZ+SuT801coG2gfO6TPPf/PkZVEo1C/42muYGJZs/P8PsO4fh7GGmXrcmoMcK3RxcUgr9EqGpXom9kxHp9hcAuWcryEmoYPK1+h11FyItjSqSDhbi28+FgP7ufS1OjzixNRtJITHqsuC+2aV3k5RjRdRXtTDmcnkRFUL0bRCDAVgUugULFixcIuhT6BfPUmrBggULFvRiUegWLFiwcIlgUegWLFiwcIlgUegWLnmEUolQmjDrrxcRmj7K1jUxDbGxaJovvvQUodFQf+hQX4thMBaF3ov0lQO6pzScPEne40/QcNzwSpMXCk2JiZwdO47Kn/Q3dLgQacnNJXPxNahr5EqHjfEJNKWk9LFUxiM0GnLuvY/Cl17qa1GMpvR/75Nz731tf2+MT0DT1LtdiIzholXoF5NyFC0t5D3xBEUv/9WDUl1rosbFvYBtWBjNyclIVnL2oqal5aL5+dsNHozXPXfjMudyABrPxFOzbVsfS2UYqtJSNE2NqKuqACh56y3yHvkbQmWmHq8mRAhB4cuvoDpXsyn4ww/wfVIudaCuq0NTX9+X4hmM59I78H/5JYQQ8nv88MMUPq+7hV9fc9GFLR5IKyPlH88TmRXPv29/i4WjA7ltfCguDhdQMsl5nM6r4vRTL3DU1pdTI6dxbZAVc957isA3/oXrvLl9LZ5WcisaWPPut/xgE4GtvR0zBvlyz9RwIn2cKXjhBdTlFQR/+gmSgf0Ze4Oqhhbe3ZbM1oRiwkuzCBsSzk3zxzI29K8U+/wnn6ThZByRmzehsO/75K8WlYZPdqXx6/E86ppVjAl159bx/Zg1xFduOq1SIVnLNd1VlZWoiouxHzwYodGgrqjA2rv3MxqFEKw4ksOyg1kU1zQxxN+V66OCWTwmCOtzMejNGRlkXXc9fi+/hPvVV3e4P//vT9N0NomItWuRbHvvvd2eWMxHO1PJLK0n1MuRq0cHccuEUBxtO9fMV9fVo3By7BSTXrt7N7ZBQdgNGNBnMeuXVBz6yZxK1n20gqGaGtYOmMnRzHLeOfIVbtOmctkrT7Z9oS4UDqUUc8f3J/BwsmHuMH/yKhs5eTqdexI3Y3v/Q9x2VRT2vV23owtSimt58p3feWfDGxy5cinHxl3OzqQSWtQalkQH80DlSeya6vF56KG+FrWNsrpmbvnqCGmldVw53J8lnz5NCbY8NuURLh/qx/NXDCHc2wnR0oKysBDbfv0QGg3K/HxsQ0K6foAZUGsES787yr7UMi4b7Iu/mz27zpZQW1rBTeocFv/9ToaFaC9hUPH995R9+hlha37FNrj3qlULIXjht3hWHslhbKg7w4PcOJxRTkpxHQN8nXn+yiHMHCT3oFVVVGDl5tZ2smul/shRWjLS8VjSuXSwufh6Xwb/2phEf19nJkd6cSa/mhM5Vfi42PHU5QO5PioEq3YZrrkPPQxqNcGff6ZTaZe89180DQ34vfB8r25s9Cl0+SjRB3/GjRsnTMGJlEKxatFScd+NL4nZ7+0Wu88WC41abZK5e0r+jj3iz+jp4oaXfhYVdc1tn2eX1Yv7lh8T/Z7ZIKa+vVMce/AJUfbdsj6U9C8amlVi1nu7xbh/bhcpm3YIdV2dEEKI0tom8eK6MyLiuY1ixCtbxDf7MkSLSi0ak5JE4WuvCXV9fZ/JrFJrxI1fHBSDXtwk9qeWCiGEaM7OFhUn4sSHf6aIoS9tFv2f3yje2Jgoqhtb2u4rX/6DSBo1WjSlpfWJ3P/dliz6PbNB/Hg4q+0zpUottr36PxE/aIiY/tDX4oV1pzt8d1ppSs8Qxe++JzQaTW+KLH46ki36PbNB/HtTYtuzNRqN2HymUMx8Z5fo98wG8cAnO0R6Sa1B8zXExYncvz0qVJWVZpP5QGqpCHt2g3jwx1jRovpLNxzLLBfXfnpA9Htmg5j//l5xKL1MCCH/eypW/SwqVq7UOadGoxFFb/9HFLz8itnk1gVyyRWtevWiV+hCyD/cLfGFYtp/dopbbn1D7Jg+T6QkpJts/u6gVmvEC698L1ZMWyQSMoq1jtmbUiLm/OdP8fWM68Rndz8vUosNewnMyXNrT4uwf6wXe1NKtF5PKaoRt359WPR7ZoOY/d5ucfTdT0XKlKlCVVPTy5L+xWe700S/ZzaIn4/laFVwxdWN4u+/xImwZzeIcf/cJlYeyRYqtUYoS0pE6Weft93Tunj1BgfTykT4sxvEEz+f7HRNo1aLkkPHxCu/x4uI5zaKka9uFcsPZgqlSvtGRVleLrJuvU00nk02q8ypxbVi8Iubxc1fHRJqdeefc7NSLb7enij2jp4gXl/woPjXhoQOC6g2KtesFWlz5wlVdbVZZC6vaxYxb2wXl727S9Q3Kztd12g0Yn1cvpj05g7R75kN4v7lsSK7zPDNSet3pzkzUxS/+55QNzSYTHZdXPIKvZUmpUqs/mCF+HHGNaL/s3+Il387I8qr+mbn+Pk5JbP8UJbecS0qtfh6X4YY8dImEfncRvHBe6tEygMPC2VpaS9J+hebzxSIgU+tFYdmzBVVv/+uc5xGoxHbEorEtP/sFP2e2SAe+my3SC6qkXctb74lGk6d6jWZT+dWif7PbxQP/hgrNBqNKPnoY5H/3PNCo1JpHdt+R3YwraztmrKiQiRPmSLKl/9gdplblcyMd3aJuqaOSkaj7Pj3s4U1YsmXh0S/ZzaIOf/dIzacKuikTBvi40XqzMvaThrm2LU3KVVi/vt7xejXtoqi6kad49RNTSLro8/Ee//5qW0BXXE4WzQpO/8+WtE0yycQjVotCl55RTScPmMSmTUajbjn+2NiwPObxJm8Kr1jG1tU4sM/U8T8h78Sd9z+pnhzQ4Ko6WIxak/599+Ls+OihLKioqdid4k+hX5hGZx7iJ21Fdc9ejPz/ljJjeP78dOBdGIvv5KfXvgfuRW91+z45MFTnP3ka64Y6sut4/VXzrOxUnD3lHB2/uMybogOIeHEWTJjz/DS9kzOFtW0hayZm7zKBp5Zc4ZxvnYEDAzDJiBA51hJkrh8qB/bnpjGM/MGs7ugkTn/28tT72+kbN1vNCXLoXVCrTZrHHVds4q//XQCb2c7/r14hGzr1KhBrepktwUYEezG6gcm8tGSMVQ3Klny1WGu/ewgm88UopEUuF1xJQ5jxwByhIk5QgQ1GsHTq09RWa/k45vH4GT3l0OuKSWFtMvn0BgX1/bZIH8XVtwzns9uGYtSo+HhlSeY8/5eVsfm0nSuM5DDsGFEbtmMXWQkAMVvvknhSy+bNBLpzU1nSSys4Z3rRuHnqtuZrLCzo98jD/Dk0zex/uEphHk58fy6M0x9exef70mnpqlzPkCrY1SZm0vt9j9pyZJb9QmVqkffnx8OZ7M9sZh/zBvE8CD9JaHtbaz426wBfOyex2Nn1vD9ziQmv7WTNzcnUVDVdX9az9tvJ2LTRqw9ZAd84auvUr1+fbdl7y4XnVPUGM6m5JHw0r/4xWkAx/wGc2WoI9eKfGJuWYSjm3kK8JfWNvPtXc8yK3kPEZs24RlsXOecgqpGvtqdysrYfJqVar458CH2w4cz4v23cbU3TyXEFpWGG744RFpJHRv+NoUwbyej7q+ob+GHQ9l8fyiL+upawvzduSa6H1cUxdH87Zf0+/57bPxN20FICMGTv5zi97h8frp3AuMjvDpc6yr6oEmpZuWRHL47mEluRSNB7g4sHhPE4rFBRPo4U/K/9yn/6isG7N1j0kiSD3ek8t/tKby2cBh3TArrKFNyMiXvvkfQO//Byt29071qjWBzfCEf70zjbFEtrvbWXDM2mOvGBTMs0FWOihGCknfeBbUav+eeBaB25y4co6Owcuned37zmUIeXHGCpZPCeHXhMJ3jyr76CocRI3GaML7tMyEE+9PK+GJPBvvTyrC3UTB/eADXRwUzIdyrU6ldTWMjko0NkrU1latXU/7V14St+BFrIzucxedXc82nB5kywJtv7ogyOBpFqNU0p6WR4uTH53vS2RJfhCRJzBrsyzVjg5gxyLfLIAZ1XT05d9+Fy6zZeN93L0KtpikxEfvhw00SFXNJRbl0h4KqRn44nE3JqtXcc3glT855moGTxjDbR8H4SG98w/S1SDWc6kYld3x7lLOF1axZ1I9h0bq//F1RUd/C2mPZ5C1fQWaLNQf7jWVquDv37P+BwKW3ETZjkklkVqk1PLYqjt3HM/haOkXMy3/Hyq17DS4aW9T8FpfP6thcTuRUMbY0hRvLTlH5+AvMHuaP196tiKZmPG+9pUcyCyH47/YUPtqZxhOzB/LY7AFUr1+PbXgEDiOGGzWXWiPYllDEyqM5HEgrQyNgRJAbc0PsmVqfw/Ali7FSSBS88AKStQ0Br73abbnXnczjiZ9PsXhMEP+9YVS3X24hBIczKlh5NIet8UW0qDWEeDowb5g/s4b4MSbUHTtrWekoCwtJm3kZPo8/hvcDDyCUSlSVldj4+hr0rKOZFdz+7REG+buy+v6J2OpoS6iuqyNz8TW4zJ6N3zP/0DomPr+an47msP5UAbVNKnxc7Jg9xJfZQ/yYFOmNg21HRVm3bz81mzYR8O83kCSJylWrEGo1nrfo//5kl9dz/eeHUEgSGx+dgpdz17XNlcXFKBwdOy16eZUNLD+UzbqT+ZTWNuNib83lQ/yYPsiHqQN88HTSHXYp1GokKyvqDx8mZ+mdBH/yMS6zZqFpaUHRg3DN//cKvRWVSk3slv2sa3RnW1Ix1x5dy4KMA7x8zweMjvQlhkoGBXsSMWao0UX6E1MKOPjsa3wRfhlv3DWdOcNMsyMVQnAip4qtCUWc3H+SRzZ/xKcjF1MwLIrZHhpmndyC+9KlDBw7VOfLpgvZzHKaA2nlvBtSz/Cv3yLk889xGh/TY7nTS+v47WQ+fyaVkFQom43+eeIH/BQq8l96h7GhHvhtXYtTRDgul800eN7qBiX/3JjIr8fzuDEqhLeuHQFCkLloEbbhEQR/+EG3ZS6paWL9qQI2nikkLrcKIcDD0YZx/Ty5LnYdvh5ODHzhHzjbWZP3xBM4T52G+zWLAf2ngmaVmk92pfPRzlQmRnjx7dLoDru8ipUraUlLw+/FF40Of6uob2FbQhFbEoo4kFaGUi2ws1Ywrp8HEyK8GBHkysCqXHwjQrHx9aX+8BFyli4l9NtvcJo0CVVlJZraWmxCQjrIr9YIfj6Wy2t/JBDk4cAv90/EW4dibP23q+vqkGxtu1RWTUo1WxOK2JZYzJ7kUuqaVdhYSYwIciM6zJOoME+GBroS6GbfQabchx5GqFWEfvEFAGVffIld/0hcZs1qG7PrbAlPrT6FRgh+uX8iA/26PpUIIci64UaQJMJ+XqX196jWCA6ml/HbyQJ2JZdQUd+CJMHwQDfGhrozJtSD0SHu9PPqHLuurq2ldtt2XK+Yj8LBgYYTJ3E8Z9rrDqboWDQP+AC5Y9HXQoi3zrtuBywHxgHlwI1CiCx9c/Z1PXSNRhB/8CTJ+4+zxnM4p/OqeXrvVwTWlfHkFc8R7uPEwvQDeLo5wMJr8XK2xbu2HBcPVxTuHmiEoCwlnayyejZXWHPmxFm+2vEOzU+9SMxdN5hFZiEEZ4tqOJRaxuGsSpT79/LwwR94fPqjFLv7sag+natOb+HUPc/iHRlKcEs1bnWVOI0ahZOjHc0qDdWNStJK6tifWkry3mME1Jcx72+3cUNUCKrycqy9vLoWxEgKqhrZebaEo5kVnEkrIrNeA0Lw7Z9vkRo8hNhr7yfCx4kZ37yBatosnK++Gmc7a2zSklEEBVNj40BWeT2H0sv541QBdc0qHp0Syp02RTiNGY21jw8t2dlYeXlj5WycuUgXlfUt7E0tZV9qGSdyKsko/SuzMcTZiuf3fknViCgaFy/B10Ei/G+3YnPPA9gtWoykVlO7dy9Fgf05UiX47WQ+JbXNXDM2iH8vHoG9jRXq2loUdnZItraUfPABzYlJBH/0YY8SbWqalBxOL/+/9s48uKr6iuOfQ5YHCVkIhCSIkLAIAnVhE6y0bmxWpUyxRaxLwTIddaqOy4gUB+lUoLRaO1UqtdgZZWwVrVC0Iog6qIggewghhDWBkDyyLy9vO/3jXuQ1JBAI5N7E32fmTu5vycs3573fufed3/39Dhv2n2BD/gn2FJ1ajZya4GFAWgIDOtQydO9GmDyFlG5dSFz7Hzq8sJC4t1fg757B8S+/xvv1Fl7pchU5pX7G9E7k+buGk5rYeMLkonnz0HCYjLlzz0uzPxjmq/0n+CLfy+aDZewoKCcQsnxSQsdoBqQl0K97Z3okd7KOTkJScmc6RgmBqZOR637AifseJLuggozfPs6K1Cs4OOJGXrl7GBnH9hObmUlUYuJpf1dDIfyHDhGTnk6HuDiq139Oh44e4kaMOKvmUFjZWVjBp7nFfLX/BDsKKqj1W/MZcbFR9E3tTJ/UePqmdiY9qSOpCR5SO3vo7IkmqoMQLisjPq1bkxfIs9Eihy4iUVg5RccCBVg5Re9U1d0RfR4ArlDVX4nIVGCyqp4xQ7LTDr0hobCSt2kn+blH2JScyX5vDbctW0B1SHh2lJX/8S+fPI+3YxJzR88AYOlH88ntcimv33w/U4b15J7+cXTLar1FKuGwcqCkipyiKrKPVVH/+Xq+9+X7zBnxC6olmjv2rmP67g+YfOvv8EV7mJq7lml71vDj2+eTGO9hzuE1DNq0lsvWf9boh/5icbS8ju1HytlbVMWBghKyy4IUFlcy+4tX+bLHEN7PuhZP0M97q57mtUETeeuym+gU8LF07QK2jruTHz46g/6hSvLHjiPtmTmkTJt20TWX1fjZVlBOdmEF+SU17CuuJr+kmlp/iAR/DdOz32d9jyvZkjaAjGovS9cu4I9Df8ZnWSO5JSXEzLeeo9eC50i44QZqt2zl0LRp394layAAUVEXfHFKRV2AnGOVZB+tJPuopftIaS2lNf5v+3SvKeWqkjw+6j0SRPh5zmp+mreOeQ+8zL1j+jFyzZuUvfEGA7Z8g3TogHfxYio/XE2fFe8BcHzRIut1Hn/8gsSHfYEQuworyCmqIreoktyiKg54a/BW+0/vrIonFKA+OhZPsJ7525YRO3Y84x67n6jaavaOvIbuTzxB1xnTCZaVcfAnU+j+5BMkTpiAL3cvByZNosei35N0220t0hwMhckrrmbr4XL2Hq9iv7eG/OJqCs8wmfrA9X15csLA8/p7LXXoo4G5qjreLs8CUNX5EX1W23022Emli4BUPcOLu82hN0VdfYCS6gDemnpqPv2MaqKoHWx9XcrI2056r3SyRg91VdqqUFg5XumjaN9h6vLy8F5+NXX+EF1yt5Gcs52ujzxKZtc4KCtFYmNb1Zk3hapSWuPnWIWP45U+amvqiNm6mfruGXTq15eMqCApy/5G8o8mEj9qFBoKUbtpM3HDhiIxzqTNU1XKawOUVNdTUlWPt7qeYEgJ19fT9dhBkvv0YuDgLGK8xXiXLCF5yh10GjKYYGkp5cvfIXHcWGIzM1tdd3V9kCOltZTV+qmsC1LpCyCAJyaKjKSO9O0UJiXNmgiu2bCBup276DbzlwCUv/Muvt27SZ/zm1bV7AuEKKrwcbSijsq6IPXBEL5AiE6x0XSJi2FgeiKpCafueMM+HzUbNuDp04fY3r0JlpZSvHAhSZMmEX/ttYT9fipXriR+zBhi0tIumuaSqnqKq+opqfJR6w8Rtj3i5RkJDO5xfnNVLXXoU4AJqnq/Xb4buEZVH4ros8vuU2CX8+0+3gavNROYCdCrV69hhw4dOq9/yGAwGL6ruCYFnaouUdXhqjo89RwfQzIYDAbDmWmOQy8EIgPDPe26RvvYIZckrMlRg8FgMLQSzXHom4D+IpIlIrHAVKDhEqiVwL32+RRg3Zni5waDwWC48Jy+EXADVDUoIg8Bq7EeW1yqqtkiMg9rT4GVwN+B10VkH1CK5fQNBoPB0Iqc1aEDqOoHwAcN6p6JOPcBd1xYaQaDwWA4F9rV5lwGg8HwXcY4dIPBYGgnGIduMBgM7QTHNucSkRLgfFcWdQO8Z+3lLEZjy3G7PnC/RrfrA6PxXOmtqo0u5HHMobcEEdnc1Eopt2A0thy36wP3a3S7PjAaLyQm5GIwGAztBOPQDQaDoZ3QVh36EqcFNAOjseW4XR+4X6Pb9YHReMFokzF0g8FgMJxOW71DNxgMBkMDjEM3GAyGdkKbc+giMkFEckVkn4g85QI9l4rIJyKyW0SyReRhuz5FRNaISJ79s4sLtEaJyFYRWWWXs0Rko23Lf9m7aTqpL1lElovIHhHJEZHRbrKjiDxqv8e7RORNEenotA1FZKmIFNtJZk7WNWozsfizrXWHiAx1UOMi+33eISL/FpHkiLZZtsZcERnvhL6ItsdEREWkm112xIbNpU05dDu/6UvARGAQcKeIDHJWFUHgMVUdBIwCHrQ1PQV8rKr9gY/tstM8DORElBcCL6hqP6AMmOGIqlO8CHyoqgOBK7G0usKOInIJ8GtguKoOwdp5dCrO2/AfwIQGdU3ZbCLQ3z5mAosd1LgGGKKqV2DlLJ4FYI+dqcBg+3detsd9a+tDRC4FxgGHI6qdsmHzUNU2cwCjgdUR5VnALKd1NdC4Aiuhdi6QYddlALkO6+qJNbhvBFYBgrXyLbox2zqgLwk4gD1RH1HvCjsClwBHgBSsXUpXAePdYEMgE9h1NpsBr2AleD+tX2trbNA2GVhmn//fmMbatnu0E/qA5Vg3FgeBbk7bsDlHm7pD59SgOkmBXecKRCQTuBrYCKSp6jG7qQi4OJlom8+fgCeBsF3uCpSratAuO23LLKAEeM0OC70qIvG4xI6qWgj8Aetu7RhQAXyDu2x4kqZs5tbxMx34r33uCo0iMgkoVNXtDZpcoa8p2ppDdy0i0hl4B3hEVSsj29S6lDv2fKiI3AoUq+o3TmloBtHAUGCxql4N1NAgvOKkHe049CSsC08PIJ5Gvqa7Dac/e2dDRGZjhS2XOa3lJCISBzwNPHO2vm6jrTn05uQ3bXVEJAbLmS9T1Xft6uMikmG3ZwDFTukDvg/cLiIHgX9ihV1eBJLtHLDgvC0LgAJV3WiXl2M5eLfY8WbggKqWqGoAeBfLrm6y4Umaspmrxo+I3AfcCtxlX3jAHRr7Yl24t9tjpiewRUTSXaKvSdqaQ29OftNWRUQEKwVfjqo+H9EUmWf1XqzYuiOo6ixV7amqmVg2W6eqdwGfYOWABec1FgFHRGSAXXUTsBv32PEwMEpE4uz3/KQ+19gwgqZsthK4x35SYxRQERGaaVVEZAJWCPB2Va2NaFoJTBURj4hkYU0+ft2a2lR1p6p2V9VMe8wUAEPtz6hrbNgoTgfxz2Py4hasWfF8YLYL9FyH9ZV2B7DNPm7BilF/DOQBa4EUp7Xaeq8HVtnnfbAGyz7gbcDjsLargM22Ld8DurjJjsCzwB5gF/A64HHahsCbWDH9AJbjmdGUzbAmwl+yx85OrCd2nNK4DysWfXLM/DWi/2xbYy4w0Ql9DdoPcmpS1BEbNvcwS/8NBoOhndDWQi4Gg8FgaALj0A0Gg6GdYBy6wWAwtBOMQzcYDIZ2gnHoBoPB0E4wDt1gMBjaCcahGwwGQzvhfwY5SEG3e4l1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Simulate the model for a given time step size.\n",
    "# SRT reran this many times, adjusting the time step.\n",
    "import scipy.integrate\n",
    "\n",
    "alpha_1 = 1\n",
    "K = 0.61\n",
    "n = 11\n",
    "gamma_1 = 0.15\n",
    "alpha_2 = 1\n",
    "gamma_2 = 0.15\n",
    "alpha_3 = 1\n",
    "gamma_3 = 0.25\n",
    "\n",
    "params = (alpha_1,K,n,gamma_1,alpha_2,gamma_2,alpha_3,gamma_3)\n",
    "\n",
    "# Find a point on the limit cycle and start all subsequent simulations from that point.\n",
    "sol = scipy.integrate.solve_ivp(fun=goodwin_oscillator, t_span=(0,300), y0=(.1,.1,.1), method='LSODA', rtol=1e-8, atol=1e-8, args=[params] )\n",
    "y0 = sol.y[:,-1]\n",
    "\n",
    "# Run this with a time step of 0.1 (which is the biggest that looks like it is still qualitatively correct)\n",
    "# Run this with a smaller time step - the smallest you are willing to wait for, as well. \n",
    "t = np.arange(0,150,.1)\n",
    "sol = scipy.integrate.solve_ivp(fun=goodwin_oscillator, t_span=(0,150), y0=y0, method='LSODA', rtol=1e-8, atol=1e-8, t_eval=t, args=[params] )\n",
    "(t,y) = forwardEuler( goodwin_oscillator, t, y0, params)\n",
    "\n",
    "plt.plot( sol.t, sol.y.T );\n",
    "plt.plot( t, y.T, ':' );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine how the error changes with the time step. Since we don't know exactly what the region of absolute stability is, we want to make sure we plot every solution (and not just the error we compute from it). So plot each solution in one figure, but keep track of the error in an array. Then plot the error in a second figure. Here, I compute the error by first computing the Euclidean distance between the true solution and the approximate solution at each time point, then take the mean across time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'np' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-4-5a16665a452a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# Determine how the error changes with the time step.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;31m# range from smallest timestep that doesn't take way too long to run to largest time step with qualitatively correct answer.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mdts\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1e-4\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1e-1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m \u001b[0merr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mdts\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'np' is not defined"
     ]
    }
   ],
   "source": [
    "# Determine how the error changes with the time step.\n",
    "# range from smallest timestep that doesn't take way too long to run to largest time step with qualitatively correct answer.\n",
    "dts = np.linspace(1e-4,1e-1,10)\n",
    "err = np.zeros( (dts.size,) )\n",
    "\n",
    "for i in range(dts.size):\n",
    "    t = np.arange(0,150,dts[i])\n",
    "    sol = scipy.integrate.solve_ivp(fun=goodwin_oscillator, t_span=(0,150), y0=y0, method='LSODA', rtol=1e-8, atol=1e-8, t_eval=t, args=[params] )\n",
    "    (t,y) = forwardEuler( goodwin_oscillator, t, y0, params)\n",
    "    plt.subplot( dts.size, 1, i+1 )\n",
    "    plt.plot( t, sol.y.T )\n",
    "    plt.plot( t, y.T, '--')\n",
    "    err[i] = np.mean(scipy.linalg.norm( sol.y-y, 2, axis=0) )\n",
    "    \n",
    "plt.figure()\n",
    "plt.plot( dts, err, '.-' );"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}