{ "cells": [ { "cell_type": "markdown", "id": "03aa22e7", "metadata": {}, "source": [ "# Stochastic Response\n", "Richard M. Murray, 6 Feb 2022\n", "\n", "This notebook illustrates the implementation of random processes and stochastic response. We focus on a system of the form\n", "\n", "$$\n", " \\dot X = A X + F V \\qquad X \\in {\\mathbb R}^n\n", "$$\n", "\n", "where $V$ is a white noise process and the system is a first order linear system." ] }, { "cell_type": "code", "execution_count": 83, "id": "902af902", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy as sp\n", "import matplotlib.pyplot as plt\n", "import control as ct\n", "from math import sqrt, exp\n", "\n", "# Fix random number seed to avoid spurious figure regeneration\n", "np.random.seed(1)" ] }, { "cell_type": "markdown", "id": "d020a2ec", "metadata": {}, "source": [ "We begin by defining a simple first order system\n", "\n", "$$\n", "\\frac{dX}{dt} = - a X + V, \\qquad Y = c X\n", "$$\n", "\n", "and a (scalar) white noise signal $V$ with intensity $Q$." ] }, { "cell_type": "code", "execution_count": 84, "id": "60192a8c", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEGCAYAAACQO2mwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABNGklEQVR4nO2dd7wdRfn/P89t6b2Twk0gCWkkJDGAEDoYCV9KBAUBsSIqyle+3x8iIN8girGiKCoBERUBkSbSq9QUUklCIAkhCSmkkn6TW878/jhnz9kyszvb95z7vHnxyj3bZnZ3dp55yjxDQggwDMMwjIyqtCvAMAzDZBcWEgzDMIwSFhIMwzCMEhYSDMMwjBIWEgzDMIySmrQrECU9e/YU9fX1aVeDYRimrJg/f/42IUQv2b6KEhL19fWYN29e2tVgGIYpK4horWofm5sYhmEYJSwkGIZhGCUsJBiGYRglLCQYhmEYJSwkGIZhGCUsJBiGYRglLCQYhmEYJSwkGKYMefLtTdi5vzHtajCtABYSDFNmbNzZgG/dtwDf/PuCtKvCtAJYSDBMmXGwOQcA2LCzIeWaMK0BFhIMU2ZQ2hVgWhUsJBimTOGVh5kkYCHBMAzDKGEhwTAMwyhhIcEwZQaxU4JJEBYSDFOmCLBTgokfFhIMU2YQxzcxCcJCgmHKFI5uYpKAhQTDMAyjhIUEwzAMo4SFBMOUGRzdxCQJCwmGKVPYJ8EkQSaEBBHdTURbiGipadt0ItpARIsK/5+ZZh0ZhmFaI5kQEgDuATBFsv1WIcS4wv9PJVwnhmGYVk8mhIQQ4lUAO9KuB8OUA+yTYJIkE0LChSuJ6O2COaqb7AAiupyI5hHRvK1btyZdP4ZhQpDLCfx9zlocaGpJuyqMgiwLiT8AOAzAOACbAPxSdpAQYqYQYqIQYmKvXr0SrB7DMGF5eulHuP7Rpbj1hRVpV4VRkFkhIYTYLIRoEULkANwJYFLadWKYLCEqILxp78EmAMDH+3i97qySWSFBRP1MP88DsFR1LMO0JqiCnBIVIOcqnpq0KwAARHQ/gJMA9CSi9QD+D8BJRDQOgACwBsDX06ofw2SRSuhfjXvgpIXZJRNCQghxkWTznxKvCMMwqVBBylHFkVlzE8MwDJM+LCQYpsyoBIe1QQXdSsXCQoJJnV0NTVixeU/a1Sg7KqmDZXNTdmEhwaTO5+6YhTNufTXtapQNlSQcmOzDQoJJnXc/Yi0iCJWwxnUl3EOlw0KCYZgMwPamrMJCgmEYhlHCQoJhygz2STBJwkKCYcqUShAWlXAPlQ4LCYYpMyrR2cshsNmFhQTDlCmVJyqYLMJCgmHKjEoy0VTQrVQsLCSYzFBJ6SYYf7C1KbuwkGAyA8sIPSrqMfFLzzwsJJjMkOMOwxeV8LiK60mwKpFZWEgwmaEC+rxEqESzHC86lF1YSFQo+xubcaCpJe1q+KIC+76Y4QfGxA8LiQpl5I3P4vifvpx2NXxRifH/cZC1pySEwL2z12LX/qa0q8LEAAuJCmbb3oNpV8EXrEmUJ0s37MYNjy3F/z602Pe5/M6zDwsJRsoNjy3BObe/kWiZ3GHokbXn1JTLAQC27gk+KMma43r9x/uxZtu+tKuRCWrSrgCTTe6dvS7xMtnc5I+sCIuaqnwP35LLSIUkrNi8B0N7dwRpSiPDVLtmxtQ4q1UWsCbBZIasdHrZJ1sPqrogJJpacr7PTSJS6/WV23DGra/iH299GHtZlQgLCSYz8DwJf2TladVU5buRMJpEnNam1dv2AgCWbdwdYymVSyaEBBHdTURbiGipaVt3InqeiFYW/u2WZh2Z+MlKp5d1siZLDU2iJWsVs+Flznzz/W2Y+8GOhGpTPmRCSAC4B8AU27ZrAbwohBgK4MXCb6aCyXgfwyiozrhPQldL+fydc/DZO2bFWpdyJBNCQgjxKgC7CD8HwF8Kf/8FwLlJ1olJgWz2MZkja4/J6ISbW/zXLGv3wjjJhJBQ0EcIsQkACv/2lh1ERJcT0Twimrd169ZEK+jGR7sO4Kklm2IvZ8XmPXh22UeWbX/4z/uxlxsH7JPwR1bScxi1COWTSCAGNiOPq+zIspDQQggxUwgxUQgxsVevXmlXp8jnZs7CN/++AM0BIj78cMatr+Lrf5tf/L33YDN++sy7sZYZF1F8wwebW7C/sTmCK2WXrHV2hrBqzqi5KXOTMMqMLAuJzUTUDwAK/25JuT5SWnICOcnHsXb7/hRqk53RZRCiqPuZv3kNI298Vrpv484GLN9UOREuWXnTJU0iSAhstHVhoifLQuJxAJcV/r4MwL9SrIuSw657Cpf9ea5yf9LfQDl8c43NOWlMfRR1f3+repbsJ2e8hE//5rUISkmXrE46DKJJJJkqPJtPLftkQkgQ0f0AZgEYTkTriegrAGYAOJ2IVgI4vfA7k7y2clvaVSiS5MisobEF33vobezc3+jrvGE3PI1Tfvkfx/akR5Wf+cOb+NeiDQCAXE7gsrvn4vUMvctywXhvWY9uYq0lGJlIyyGEuEix69REKxIDldwwH5z3If4x70O0ra3CTeeM9nXuhzsaHNuSNpXNX/sx5q/9GOeM64+9jc14ZcVWLFj7MZbc9KlE6+GX7LWp4D6JJN45uyTCkQlNopIJYxrYsLMBH+7w6dtIoQOJagCZub4v42RFWEShSfCiQ9mFhUTMhPmQj5vxEib/zN+aEH6E0rPLPkL9tU/iqgcWov7aJ/2VIwQemr/ed5nu14zkMpplyQvLSL/rSlaEg52smptKJFe/hsbyWvDLDRYSFYZXBzJ79Xa8/F4+UOze2WsBAP9atNF3OYvX78KSDbu0ytQlSYesvc7lOI7NSiRbNmqhJmktZcXmPRhx4zNFf1e5w0KilXHhzNn40p/fCn0d80gpqk4iyYGoqqigHe/D89cntjKbWZg2t+RwsDndUWtGZJUUWXh63LxTSCT40ruZjNr3DQuJmEn6A0qquCrT4CwyTSLkhfys6R3l7O6Vm/fgf/65GFc/uCiya+rymT/OwvAbnkm8XDNRaIBxOJf/tWgDhlz3FNbuyIdFZ1mYZRkWEjGTdEx71CYIIQSeW/aRxd48/fFl+NzM2eajIior3Pl3vLI6dFnmzfsbm/H2+p2e12ooCKfNew5olx8Gc90Xf7gzkTLLkSffzqfFeXfTnlTK1zXjjpn+LKY/vizm2gSHhUTMJD16iVq7/vfbm3D53+bjz298UNx2z5trLMdkZYQmM7vkcgI/eGwp3t+617LdLryNX+Z7+fZ9C3H2797AngNqM9JfZ60J3VG/smIr3vvIf0cme+zNLTk0NsebCsZRjxDvP862Uy6hr3sONDu+qSyRiXkSWeGBuetw7SNLMPe6U9G7c9u0qxOIqDWXzbvyo+NNu+IfJYc1Ack6hVVb9+Jvs9di1urtlu32omRFL1j3MQC4dro3/iv8CPCyu/Mz9qNYKnPKb17Dqi17sfqWM9GcE6iriX8cGEVHH2d/npExTNnCmoSJhxfkQzo/0FgAXQihlW018bQcERdoCB23jzg6n0S482VRLNqdj6TsUsoIvauYyz/v92/g+48skR73r0Ub8MaqCGZ2S+q8akteY7r6wUUYdsPT4cvQqkb8rfzW51fgl8+95+sce3uIQ2sZdv3T+O2LK6O/cIZgIWFCN1TujVXbcN2jS7SyrSYdpuinOOdo2nmyscmtn4zKCRz2Kn7MC457L5Ru7vCK9657TdO5C9ftxP1z10mPu+qBRbj4rjnadXWUo/GgHgsQ1pxlfvPiSvz2pVW+zjHaQ1Tf4IrNe3DfHOs7bWzJ4ZfPr4jk+lmFhQSA97fuxRf/PBeLNZyUAHDxXXNw/1y9RdWjaJ4PzvMuy/gQ/Izq7MfmRP46jxQ0qvwxedxG00Hv8SdPL7f4EZIUqHbBZsz5MJO0gH9j1Ta8siI7a6LoEsonEaMWErVPYuptr+G6R5d4tguvci/90xz8+oXyESwsJADsPdCM/7y3FQcLtues2TCveehtz2OMdhvGcS2EwIvLt+DqBxc79sVhbrrjldV48K2SAAzrdHevo9xRbXDpn+YWjnOeq6sphZ20dfFdc4r+CTeSjJjbfaBJy/waFB1NNbKyQj63psLKe40h14h5beU2/PqF8jFRsZBAvA00qcFoKToneIEC+U7Bsq2oSridVypzw84GLNvoHJWraLQseRnWcS3xSSjqrfOcjCOymm0iiWqd/4c3cfIv/uNejygc17YXtfdgM3Y15Nvix/v8ZRkuXrPQaKPK0lxbnb/e/oOVk3JDBxYSiNnBlVAHY4x2w4Yj2jvVkuPaVUoUOW7GS5h62+uByw+DL5+Ej4Oykv7CThL1WrF5r/dBmgghHDmNVHcw/ubnMfam5wAAR//kxWAF2ttyyMfVtrYaALDfx6TNSoCFBJydS1h10kxSpgEhrP8GugaEUmC6dcBR3WFox7WLILNfW2cyXdY0iT0HmrB80+7A73jz7gOO+SK6uAkk3TZ+75x1GHHjM1j/sTOzsf3NmcOO45738caqbVoJLtvXFYTEweiWx83qAMQMCwk4O8DL7p6LRxeulx/sk6TaQFGTCNHVyjQJAz/2/qDYbf9zVm/HF+6eGzK7qBHiYt2qZW4qHJMTAg/PX4/6a590mONkfNLHyPe1lfqO6i/cPdeysp7OU7nHNAny6FtexKm/fEW7PDNuAyelwBUC0x9fhiXr8+bHpwozoJNa2tfeZlXP67llH2ldr31dflrZ/oI2FEW7123b/3hrHeqvfVK6omPcsJCAfAT67NLNKdQkOEEc17IJZfboGuNDCOK3mfX+ds8Ea+YPzV6fK+9fiFdXbMX2vQe1ypPVUe2TUFXIuSknBGa+mk/5sXFnabEkVSex0cfEQ8NhrsPCdTtVVVQy/d/v+DjainktE7fRvKo++xpbcM+ba3DhzFmF45xzbrIwkNatQruCuWlfY16TiKLuut/rjKfz4fZ7DkSnxejCQgL+Ohe/JPUNFOP8QzmuBR5ZYE1vbDRiP6YcAHhx+WZcdOds3G0aycowz+QO+9FVRe2TKJBT9I/2UWBSaSDOvf0NX8d/QSNiSsYiU7qRgwFMPsq2KHtOGU2hYb4Hw3Hd3GJo7fps2NmAz985u+iMN9CNnKsqNK4oE1PqwkICMQuJgC/Vb4pj43D7afPX7pAe7zZxTrbNfTKdc5sx4vYKnzRHQjnyKfl8dLozo/PX1o9uOuHnL+O9zc7cSlnxVXjxasC5FzUmqWvWJKbe9pplrQTveQNUOK7wOyGJYG8PQT5Fyzu2ddR+vu3fvbQSb76/vZh00EDX3GTcCguJlKiShU5G1JCDvtIWn42h1GCt533mD7Owv7HZkYQuP3HOdg3ZdbXScgRvuOZvRHmZCF6Fw3GtPE5t/rKTxAe7aVcDPnfHLKnJLe7iq01CwqxJLNu4G1c9sKhUD8X5xva9B5vx5vvbTBMz9evgJ/27He2Z8i7P0fyOjesFMe0aAsv+Xet+50WBZzu8JSewfNNu/YoEgIUEFI0pYMfU0NgSOK7bjF9nrXG4rM19+76FOOf2NyxOV9n1XTs9nzOutTsEDSGxcN1O1F/7JD4KkGRQVbxOB+8VBJDEoO72l1dhzgc78Pji5NNsVEs0CV0N1L7983fOkTYUr2f8g8eWelc0RsztpMo2mtcJEjGel3GuY1KnphWvVLZ1+y+eew+f/s1rWCnRdKOChQSitSVP/e1rOOrm54u/vTqSPQea8Lpkso/viB6X0c3CghZhNhnIru9qbtIo24yhVntpZOaPUPXR/W1WfpnVN993nxTlGqZrKufXL6ywCSd5uVnQJD7ckTfbDezW3rEvaCTblj0HMNuWFVeGVZPIj+h1m+UzSzfh6SWbpPvkLgn5y1v+kfsoWQhhcbBbrhmFBmq6Xyqam5z7vM4v+hRsD9CsSaxw6ehVPon5az8GAGyPYGCqLDu2K0cEEa0hoiVEtIiI5sVUisYWPVZvtdrgZR/yg/M+LH5AVz2wCJf8aQ627D5gO88fbqOb0mimdFctQjiPlY708rjPkxCYt2YHrnpgYXHbnA/kvhDV9fP1zP+7eutevLjcf3SZNAuspOK/fmEl1pjCMM3fXVOLQP21T2LeGu/6JyEkdhQ+/naFGP0oOPd3b+BCy6JRcmqqSt2DYW6S37Nz2xX3LsC1tiy4RdOlj97ba5Bx52urMflnL0vX43CGwPp/XzJzk5/3bhxZZRMwBubB2hm3voo1Ch+eqmzj/Go/URs+KZf1JE4WQkQzt16C3HFdmILf2IwqouJsS99I2pORi2nNjKnF1M4NNturXzu/0XhkkThGO9xnmgTU0iIRJlIpYfgkXMxNAvjSn9/CngCTjKyaRJ5TbLH8uh+3vyywpWvKPvofPbncM6LH/sH7/UybNWLe/cyk36f5/I0Q3Xve+AA79jfh6tOHSY8zyYhiZ+SmbXohC4IIK2fffD+vEW3YuR/D+3ZyPdZYe9oP5ndszyqrIyyWbdyFIwd0VTqe7d/5x/sbUY8OjuvYnf8GRhuKU0hkXpNIApnj+t+LN+LBtz7EyBufxeSfvRz42vZmZHfEyWyNzS05/H2OPM20VzlumsQtTy0vbpM5zNZIJjlpaRICyh7Sq+M2Yv8B9Uen25G4Otftx5oqNleiNSxSrDa3ZptZA7F98F4VtHFAI6y0OJNe4+q3veQvadz0f7+D21zWQjBrEjmXjlH3vottSbJPOYnTo/0Y343xDTc254rfmF1jeVex+p/bs7VqEv7NTWf/Lh+uXK3o5O3fodHZ/+n1D7B2e0mrKAko6/nNhcrUtHIhIQA8R0Tziehy+04iupyI5hHRvK1bg4X6qR7vNQ/nR/xb9+hN5tJhty1OmiS2xnveXFOcPKOL24jT2LbXNNKUjWJl8fdaPgkfuKZ3UOzSzsIqqaSqvO89XMqs+/k79dd2uOLe+bjjlfcL9dI+TcpBjcidovD3GMFfdvdcX2t8u9HckoMQQjri13kVsmVk8+ca5qawNSxh2PiriDB/7ccYdsPTOOIHz+TLieD6Fseyw3GtT1WV3Kdg9w1WEWFXQxNufuIdS7tU+SRaikIivq68HITEcUKI8QA+DeBbRHSCeacQYqYQYqIQYmKvXr0CFZBkFlj7yEE2QtgdZFaly0csa8xb9hwsmrrcyGl82PmcT3IcdmHXL0u+U7czlpnEVKeGSX/9k6ffxb6DzVhYWN60VL4/dIIThGaHFOU6FIdf/zSufnCxxcl68V1zsGFng0JTtf7e5yNL6r89ora8nmmLaSRtTjsfFUb7b2hswdyCny0n8tkEFpm0YC+M78f+/R//U6uVoqaais/dPKgrRkfZrmvcf5x9WOaFhBBiY+HfLQAeBTAp6jLinNxj/6hkIwfAOuINojkWVWDJRyzTMs767evYttc7IqJkbnL3SejiphWoQylLO7btPYilkgWCAJUmoV83P3z3H4vwlb/4j6Mw34tXjPzb63cWTSTPLHVGCsVxa4ap5tGFGxzXX7phl1Rg27U1lW/Efuq67fuL90eF8xbYBK/spT5UyKP18b7G4jOsqiJn2wr4WZ9pyo9lXNMcWSeEwEV3zsYlf9LXQKsU5iY7NVWE1dv2Fs4pbZdZHAB3X1FUZFpIEFEHIupk/A3gDACRB04nqUnYHculqIXStuoAFTIazyvvOUeUYaJwdE71lS8qwHXMdfjsH2fhrN++7sOxH89HpBJUXpjv0UuTMOzZALRXQgzL+o/zIbcd6qolTlY9zdJ7Alz+ZLOFhAj4zv0LMe33b2LXfvckin+dtQYAsG7Hfou5yeF3Cigl3jFNTjOEkPlVBfme/jkvnzDUK5NCFRE+84d8rivzwEy1FKvhk4gz0i7TQgJAHwCvE9FiAHMBPCmEeCbqQqISEvZc+TLsL7M4wjA18aoAqoRxtmy9XTctQxf3ZySUmoZ9u7smIaQzi83nrC6YiXY1NEEIgddXbit+eO5LrEb7EdVUB/t0zPeiygmVJkZuoa7t67QjmeybvIT9ys178MDcdY5rGcsHmzPOOs2VAm+vLwlo43lWVznblu537da/Gu9IhHxv2wrt2vxs/umxLLG5+vYQ2qaWXF6Tau2ahBBitRBibOH/UUKIH8dRjk7ctj0xl4wRNzrll8OGqGjI5oYXJJzNdYQSogHpLDrk1kDveXONJe2627ECwDOStM3mU3p0qAOQH/E+teQjXPKnObh3ztpCHdV127jT/2xtN2qqg40szB2Z39QrSeAeGixKzmfLOdajlFFqhX+vfWQJrn1kieO4UqSSun4HmkwChAAjkpuIArXzf7y1Ds+9o56T05wz5oeY6xn8vZnP/X+SZYmfN9XFoknYzv+fBxfjqJufN9Wv9WoSiaDzuZvtlDJU5o9/vPUhVm0phd7ZO3OZrVFmbtrfWLLzLvpwZ3FkokOoVeuKH6H7IW77v/uPxY66SK8jgCZJWKj5lI5t81N7GppaiovXGDNupT6Jwr/h1qRwUhswmsQ8GAhdp4huyWxvz5net/1V5c1Npd+G78HL72a5gOI4QsmnUCUxsxTLNH0HBCp+dwcaWxxzjXS+6+89vMQ1erE0Ui/V1es7cvffuZ/8E1NUo1lYGs/EqI+RpsXISMtCImZk8yTsbDCtIyBjp8KOetuLK3Har17Fzv15J7F99FjK6ZJvANc8tLhoUjFz3u1vFv8+9/Y3MOXXVqHl1kiiGvmoFtzxM/HPVeERQrq4jdWpb/1YgNJHeZNp7YRRNz6Dg80tsanhUk1Cox1ZzE0Z0SRufqI0f8ao0/qPGxSpW0o2+lH/9yweWeBcnEsZgGD7bb98yWyorqt9fWmjjp+/aw6eXuq9eNB//fZ1zHrfOyWJQVNL6X6L9Qz1Pekfa34OqnkSxvcSZ0ZiFhKIxifhNbIf98N8Pie7PdPo9JpzOazcsgcPzluP++c6J9K9t3kP1m7fVwzzs5fnalcV3sd4nXvnq6tx5PTnpEn2mnNCKSSd11NXokUI6eI2llmvhX+fXrLJdSC9r7El77eIJQZI4ZOQ3JtdgJoHCWE1iajuzZxF1Px+ZOGp9iq//N5WR++vMqM5NRPhuh9wagMWTYLUz/Dzd84uOnXNLNmwC9c/tkRyhhzj+lbhrn2683o+PkKziVc2OAJQXKkuzmVQyyUtR6xE4bfer+G0bmzOKZ1rOSGKqqOKE3/+H+W+KLKaul3XSCC2cZdTo3pNkqBQxQWFyA0Z+w42o1HyDCz3Vnhef5m1Ft+bcoR5k4O8OUK7ar6oVRjO7ebEnADMSoeI0twUIY3NOdTVVFme115bKKuAszOqqbJGFQkhtEfa5g7T7FKwpsKwPmd7nVRlvfn+dvTv2k66z0/0oMzmP0+xRouKPSYN3I8WUiXRJOznN0s0nahhTQKIREroCIkZT78rWc0sX/gPn1iOZzXX2pWh00iCdJgyu3QYZIv3GOxuaJZqEhIZkd9e6FbmrtkhTc+eEyI2ISEdERI5Rq/2TjWL5iYAuOahxdi1v8lSp1dX2peydZqMaqqsgnjw959SBlF4+S50FvPZYzN5ugnaj/fL5wH5CQyRRQ/ZV280s7+xGd+5f6Fl25jpzxX/9vPKzQJSNeM6iRBY1iSg55PwoqHJe5b06m17HS/TWAxo8Yc7HQsDeWFNfxBPI/FKSBYluw80yYWE6W/zuzKqsnDdTlx0pzOraZwj9YWK2baO52XbbxYucxWZcvccaEJdTbLjt8cWbcSyjbtxw1kji9vMkUQGCwqpqQ1kvhlNv7XD9CpLd+EwN5l8EjLnuhnVV+3ne2/y6Rhe5pFE0O+Kkwaq9SSK12UhES9JmZuqyDorNOziRMNvKIXcamkSAcqwN777JP6SqNjV0ITGFudzNH9Yqu9blrytJSdJhx4zdk3CGeZZ+v2jJ5dDxpjpz2HCod08y4q6X1i5xTmIsZQHgW/83TpKrqmqcjxj3USNTTYpoeM7azIFNvzw3+9IgzwMVFFGfgLT/M5D8JoI6McnYamnYj0Kg1Y7TyIp/OS3V6EnJABz8E6UcfJalwpQnD3YyE3VDsueA81oapY4f00V9zOLdsG6j2P9eGTYO8j7bdl8dSdizbeN2GXIHLNh8ZuAsaba6ffRHS2bNT0hSmW7CXazpum1ZonqXvz4JIyoRN2Rutd8Kj+ahNVxnf9X1WdwCGzMRKFJ2Bc4l5ZDFJsJJC7Htc6aB1HRkhOO0SVgd2SW/vYyfV31wCLfKdfDYl+OcropLBfIlh9Cht/ZxHbHNeBibrL9tggJiNJ5pgPn2YSlLERahepJ+8locPnf5uev5fHaGptz2L73oDJM3MCPYDe3dcNEtruh2bFAGRCv45rNTYjGJ7Fso3cun2qi2Gz6Op3PW2u8R6d2mhIUEs05iWcUamGg8yjdloSMA6/3kHkh4WpuciILBVaHwMojc+xluz0hrwhAHYw10/3g9d6+ce98vPjuFlx/5ohQ1zFjTcuR//eKe+eHvq5fWJMAtFWJM259RZnhUuZwtVNdRbGlYoirjTQlGKbZohjGWhyZZse1xjXjdLTbIaC4dvSAbqXwy4vvml0MSshiviYzsmVADWTPMh/dZI9a0rvJZtNxeXNT/u+L71JnV2328QCjfPVen8GL724B4G1C9lOnNdv3FzV5LzNrnO3cU0gQ0RGxlZ4RdBWJFZv3WpKLmdFRg+2Tf6J8r0IAGz1mhQdBliYjCB/va8SX/jzX9ZjmnNwgpgyB1Xh+SU5FEAC+8fcFAKxhlm+s2l5c5CgtTUIVSWVHliDSjbzj2orOgAmwm5tKuK1zcstT+otxRRW0cPvLq7Tfm9diYQ/4XPNinUvKGTNxDj50NImniOhuIhoUXzXSxY+xSfWydD4Me3RTlNJfQOCh+c4UCWGJytz0l1lr8rNzXVD5a4xntm3vQUsaZ51OIElN4h2TydHuHDU0oLSS+t38xDveBwXg7fU78aU/v2XZplob3BHdZDId6c7YT4OfP/uepd0libHinJdJPG1z0xEAFgJ4hYh+TUTBln/LMH6im1RH6oxY536wAx/uKI32DadYFLTk1KvDhSGqCBqdqI5mxT0YH8DL7zond3mWm2CfbO707M5R41eSQitqZFU3zCxmVAMmt8l0slQ0YYnyUd+XcACEgREGq7vWdyx18DpACNEohPgtgBEA1gOYQ0Q/NBYDqgT8ZOYOstaDwUe7D+D/Hl9W/L1IMXkuiB89jnBIQN904IVO9Vpa5LqBoUoHGS0lPU/CwL4wfXH5yoz7JPzSs2MbxzbdCCQ//oXWijGA9dIkUvVJmCpxQAjxCwBjABwAsICI/je2miWIn9j7EDJCmyDvu6klF8sKe35CDt3QMbN4CTqv2czSc1Lqh+wftSr3TlIsCbiSnhlZ3eskM66VmoTt9CgildwoX52thE5mXCDvC3u7sGhT1GgLCSKqJ6IpAL4KYBCAPQBuiaVWCeOnc41i4l0cNLeoV4cLQ2SahIYqoYxuKvQu9hBendFTWp2yPT+QMRDJegisGzKf0UZJVmCVH8ttnkQcBE2BkSVWbdmL42a8pLUe/W0vroylDjrRTW8T0Q4AjwH4IoCuAF4CcBmAjrHUKsNkU0TEp7pHZ27y/mD3HmzGowudM7p1cwEFPSYOHELC0CTK2MKiK+BUbcZ+/jUPO1dmi5KwQQJnjOwTUU2Cc9frq7FhZ4MlnbuaeHonncl05wFYLcrZ4+aBn8l0WdUkmloEfv7se5Ffd6VLOKIfdKxWquRo6uUwy0mTyFPOmoTuwFwV3WRfMChuwj7qoEvURokfk1xcXZOnkBBCvB9P0dnBl7kpvmqEIm77blhUkxDNqHwSYTSJrAgJFGbbry3EvZcjuuYhlSbhlbIia0SRiSEsQWdoRwnPuIa/h5vVrjjJ9BlB+Mc870lE6jxR8qe+ebf3Ot9pvS/HPAkAM19d7VhroJzQ7bCUmoRGEswo6VBXHep8e4RaGvjx28Ql1FhIwJ8Jye2lpdmosi4kdPCrMTwsWV/ZeXLw+oRB5pOYtVp/beUsoqut3vPmmngroklYv3WYcPeo8CUkYurNWUjAX1irm2smTRtmXPMkskC4heezYW4iJDuxLw6ytNyqDkGXCjXIgibhK2tsTAanzAsJIppCRO8R0SoiujamMrSPdXtntZKMmEmRZErvpAnTN6XVr9mFxLodDXh1hXtakqCcckTvWK5rR5bGPcv4GR/USIbhfpY5jZruHeoAeK90ZyYuF0qmhQQRVQO4HcCnAYwEcBERjXQ/K17cRidJConp/zUSE02rl72xKjpTRgb8dRbKUZOw24e37fX2nwSlNiENNuvBEXb8hMDKBEKaQuLUAII/rsjLTAsJAJMArBJCrBZCNAJ4AMA5aVboA5flEpNsVIf27IDDepWmqTyz7KPIrp0xGRHKr5DePInkypKt6RAH+xq9I9SyhJ8Bgsy05GcFu6gJ8k5ba3RTfwDmsJj1hW1FiOhyIppHRPO2bo1HnTfz/UeWKPfVJdgzVBNlwrGWBOWoScjMF3Ghandf/GR9pOXsbigvIeHL3CTRxqoTfId2gmiHrdLcBLlwtLx6IcRMIcREIcTEXr3STVCbpOO6uopiyyOVtQmDYbr51MxNCQpwVYcStRkqK/Mc4hiMyQRCii7GQIOM1qpJrAcw0PR7AICNKdXFkySjIaqIApm3fvXZsZ7HZEtEhDMZ6cylUDH9v4K7v5IMdFP5wqI2Q+ms454EbWqj77ak5qYy0yTiGg5lXUi8BWAoEQ0mojoAFwJ4POU6OTA66yQd1706tQk0eWba+AGex0StSISVnWloAz86dzTOHtff+0AFSXYwqnZXq/Hgw044S4O2tdHXWW5uirwYbYJYJeIKg8+0kBBCNAO4EsCzAJYDeFAIscz9rORpU5N/jEmamwZ0a1cUEmeO6RvptaOOtw47ok3DYjRtfH9Px+WAbu0w7Si5IEmygxnUvb10u85zryLCwO7tPI+Liz9eMsH3Oe3iEBIZ0ySClN0SU/RZpoUEAAghnhJCDBNCHCaE+HHa9ZFRZwiJhBrViH6d0ba2ujjiP7x3xOs/2b6XsDbgsGa4NDQJAnna3dw0uSQ7mAsmDkCfzs7Ff3Q0Wy/fSdzuqUmDu/s+p32E2o/RtqUhsCn65nS0QDs9O9XFUJMyEBLlgNHQuravdT1u2vjg5gszv/7cOAAlM07Ujdl+NUNTCkpYM1warmci7w7SbX+SmkS72mr89DNHOrbr2LWryF1TC/vuvQjSciM1NxUqIBvgpZlBQVf7NrfBG88aFUtdWEhEgKFJGLMkZVx58uH45kmHRVKe0TCMkWzU/nJ75xfWURg2yiYNTaKKyNPn47Y3yVFoPtLNWZ6OcPZ6sm1q4vVZBPGrRWpuKjwAmUBIM8jP65s5r2DmNGvprTUEtiwwhEQPk5Cwp0qoouiyNBrtwhyqKlOXn75qcqDr230ScXcUXqThk6jytja5pvxI0txECoGmMxIWwluT0ImIsyMzf8mgAI8pSnOTsSaJzCRq/g5ka3nHiZeJliQaEGeBzTBG4rPuHUoN6TLbRCYiitBnQYVroviv0Ty+c+rQ4lGd29VicM8O/q9ua2tJpX1ImrEDuij3EZHnyGzdjv3KkXjS8xxl5UXRmbatrca08QN8+6Xa18mXqrnrCxMtvwNpEgHv64ufrMehPaxOfuPblQ2yzFVLOtLJy9xkPDeLJhFXXWK6bqukrcksY39hRNGl8nWuZ1Pq0MxmDp3RsF554a6S1cXY3Jy2+WcX/L6TNlXIJkC2qw3/eQf1Saju3/4NBBGmQQMpbjxrJA40tWDt9tLCT7miucn9mkkn7/TSJIzdZm2RzU0ZxhiFVLvYB6si1CSMDsFchNGhmdtWFQWTEvZTdGcPf23yYP+FpYjbh6ijSWQJ2a3ojLi9ViU2/FF+n0VOYYuzDziCCOKgGQGI1CY4WVswb+rSzj0oxatcv3gJLaOvMZuCW2uCv7LAGL27zYCuilCTsJdi/tCrbIIqSLO5ZdoYy2/d0V7cqSjMprQoUL2vCybkJxymJSTqAozeZR1EFOYmoxPy+yzWbJcv0+oQEgk+YyJyaAQzpo1B385tpdqyWYB1bhtcSAQJAdfV4NyCZaKChUQEFDUJU0Ozj5Ci9EkYDdooTgigsbCehLmtVxH5Hl1MOLQbxg/qJi3PC9WoMApr08Du7XCVhpDQOcZAJSRuOCufjsPPKNc+IS2MiS1IpyI7Jczo18DorKKaYGl/5kmvI203VV04aRBmX3eq47gvfrLe8i2FeZZug8eLJg3E108c4tjuJSR2NeTzaPXoyEKiLDBUWC9zkyFEwg64iw5ryYd7wLSOcD6M0//17RqBtpCI8XvXuZeZl05A/276s4dl9/XEt48vdgg6z+4rxw9GpzY1OGV4dAv/+BESp4/sA0Dlk6jGUYO6hqqLMfKO6t3an3nSDn6Vb0GYhjJfmzwY08+2zjno2Da4f8ft+zljVF9p9GAbjzBfYxByeO+OrsdFAQuJCJCZmxx2fQKqJcIkCG4fbKNpan5Q56u9eqHrG+rsPDpaEZG/u5V1xmYNQEcLG92/C5bc9Cn07tzWR8klOrWpcUyyDJLGRPaK2tZWe2o0Au5+CSOyLaq+3Ln2dxCfRPDy/cxCN9ctTISf25wZlQDx0iROOaI3fnzeaHxvyhGB66ULC4kIkDmu7V+VWZMIG+1j/7BUl/Nyvt731aMd24QQgUd7cY4KdS/tp9ORCT/zEp32vZOH9nQpV7tYCxPqu+FXnx3nWS8vZJ1NW5+TIDu1cY6WqyUdZhjst5aUJjGkVz4UvLbGu8CiOde0Lcz9V7sIGAKkHYLXrPIqIlx89KGxJDt0lBV7Ca0Aw9fg9nETlfaHtdEbpZh9Egbm2cn26lz4iYGW3wMlieFyQmIq0/yS4/RJ6HyjBH8jXlnH2thsEhI+LqZrkvv80YPw08+McWx/+BvHFv8O5pNwnlNdRb6e/U8k9TLadlR9ub0tERG6eaSzsROkLkcP7gEAqC3cz0nDe+GBy4+RHmsM5qz+PVsdNCrRsSB03TQJIvmkTC9NIklXDguJApPq/ScaA/Kx16XZj2rHtXn9h8E9O+CPl0zAE98+HmtmTPVdpmykY9CSE5bjjBHQ/V87BjNs+X1kDS0n0ST0R/GaBwZApxMm8lcHWThkU4tZSOhfzN6JqDrnnh3qcK4kc+yEQ7tjZL/OAPxpEsaRshFyFZEvtVVmGy/WJaJ3K+swjztcraFFhRGSawip+h4dcMyQHtJjZbcs+569+ER9N0uZMghk8YcYxJ0zyw/ZqUnKBLW7n3Vkv+Lf5objdFzny7j7ixNx39eOxpTRfTG6v3rGrxtu7dMhJEzlO6/j3Njc4hQSUaWpnzqmn3T74hvPwM/PdyaoMyO752+dfJjjGD9CQpY6wywk/KDv3FfnhDLenR9Nwng1qsg5z1cnSsfIijXq4vfrUN2DNNtqAjanloKwNIpy88PIfBIOTcJH2V6ahAwvMxJrEing1lBHFEZ4UkynWWY/2g4zGt4pR/RB707BnJzFaztGraUGbxYS5k5TNpqR3XFekyj9nv39Uz0nXBmMUQg94/xrpgzHsD6laIyTh/fCuzdPQZf2teja3hrK988rjrX8tnesj37zk/jGSYdbtlHhP11kpuLGZvW9njhMvTyurtZhT8ZnPsvoyIJ0mrJZyH7DS2WHG3XxOwdGdby0DEU9Lz56kPY1vLCn32ixtWnzT3uIef5va6ETC1qCDq5maMiVPdYkMojbi3QLiOjYpsaSlXXqmH7ShGhRzoYsdoQSR3izw9xk/C25jsLcZK5r3y5tHR+UijNG9cWr/+9kZU4ke6fVvq6mOGKyp1m359ixc9SgblJxEEaTmDS4uyMxo5k+nduif1d5iK0f57752B37Got/54qahP/PUmZu0rE25aObjOOd16gJGN2kGljIcyRZt/Xr0hZrZkzFd08f5rNUNYaQMNqge3JG5z2b2+4T3z4ed132Cc8yjftyfZ0Kn4TXhMqoFwZzg4VEATcV301dNKctrqmqwu0Xj8e08QMcDT9KjdrNTJyzCAlzY5JpEhJzU0446urHAjPI1Lkf0bcT/n3l8aaPxXphs9Zjd17q2IAdm8hvdJP19y3njXb9ON36W3v9VKGWxx7Ww1JHc+jsJwp+sS4+HbmAXLD41iQk24Kagox3O/FQ64hb9i3ZizDegSoluJ8O8raLjrLUpygkXKRE8bkpzMej+3dBxzY1WHTj6a6+TCOIxNXcpLgXTyHB5qbkcU2p4ZXjx5aV1f43EGxmqWpZSbcQ2Bahr0nIyOUkPgmfTolte/Oj42nj+2PMgC7FUaW9Ds8s+6j4d5d2JXOTbNa0bDTmmNUuKcMN5zt3P7lNTRV27m+U7rO/3jpF2OOEQ62dyi8uKGmdPzx3FF64+gT0CpCWWmZuIkDqFFUhTTdeePB+26/RZC6cZDUZyb4l1SzsKMI7jcFfydwEy28ZOUl7lbWrru3r0FaS+sRYN7xGw1RHJH9HtSkunWonOzVJGbfc++aRwATbyEiFwycRQPKrRhlF3UCy2+GTcLu+1Nzk7BB0zU0GG3Y2AAB2NzRbtruNqAxN4rozj8B3Tx/mMFfInoX9cmSai6KDs3NSH3v9mSNw+og+2Gea0W4v24zuhDhzuoc2NdW+l6I1Sg1sbhLW9mKnNE/CV7WK2DV0uUZoE/ZkLTsMdh+EUZbbuKeh8I7Nbc6PkLzylKFYM2NqsSw3K0UVkfQdxZ0HzQ8sJAq4NQLzrhOG5p2Xk+q7489fUtslZR2YX1SnuNW12SIkqOS3kF1fsq0lJxzl+tUkvlhYS8MeKUS2D+KXplF0TXUV1syYistPyEcs2YvU+WYI/j4uu0Bxe0dfO2EIqqoI3zl1qDSu315s3y4hghNcbmH1LWfiiL5OQSIzb+m2OWMkq5pr4VYpr6SL9sGX7PU4J9hF10Ea79how8XfLtJzvyEkLJYB/ToZhxpFGPczvE8n3HzOKKz68ae1rwUAlxzjdODHlfFVBguJAl5hh93a1+LGQuI3ADh6SHecXMjXo/O+grxS5Tn2HaYGb+/Qq8hxiGsBLRJzk19NwnDuGgLLHGLZz+T4PUThBAacH7Hso3AKYn/Lhtod1zry5erTh2HhjWc4ttuf2fnjB+A3F47TrosZt2qohKCq/fp5dQTg+e+eYNnmpknU92iPq23O5a9NHownvn28qV72Z+xtbvJ6DUTQfrZ2TcKojr19mX82NDU76uGnTy4dao2o+tSoPrj02HrUVFcVfRlCCM/oQXsUX9KwkCjgNgIVAlh44xn48vGDXRuL9V3rj1L9UvSrwaklODrOYt2cDVFmwrGHwALu9lsZxQ/TIbAIt5k+bj/5cGSPT1Z/e4fTz2VEbx/lhokYcYyGqwjnjHNOmosCc1sSkm1m/Lw5IsLQPlYtpTj7WPMa3z51qGX+j/0dy0xIbunDn/3vE3DlyYc79us+W3tb1Ilu2ncwr0mMMUXp+bH+2O/HqIOlTGPwBmu/0bV9Lf7zvydZzpcNAMaHTNzoh8wKCSKaTkQbiGhR4f8z4yzPbyoEayI45/4oZILqwy8aACS77dkrjWtIFQnJ+bIZ137NTUbnaxcSAnlnnzFXws3m3K9LW3zzpNJkOZ3oJgL5ymAbZm2DH507GrecV0pjEe0gwP1afkryGqVaQ2Cd+43nqZ8u3ordNyN758ceZp35bC5reN9O6K25XrYMpZDI2TXV0t/DC+a8UYd0KbZB1QBixjRnKhOVucmsvRhXywmr27pz21rU25Yctj+zNTOmYkA39xDxKMmskChwqxBiXOH/p+IsyG3heskAQHGcsxEkTdgF22U+Cb/mJqNRG+Ym+7MozSxWP3MispgyZM/Tvk1mbnLr2+xmfD/9/CXHHIrPmyZ7RWlH1zG3xIHUFGRoEtKBkLdWUOvI1eS8zqdG9cVpI3qbjvF8Ah77Sxj5k3p3yn8Xhu/G3vEaTfz6M0fgihNLgxMjFFU1njmkazvlWhMXTMznSjMmmbbIggSEfSKf8zpBcnlFSdaFRGLYUzyY6STJJe8ntBBwb9b1HhPHnGXbfrtUxc0nIatTXkjYNAmX6z991WTHtlLYYd5xbT+9KCQ8zE3mDkcrKgb+nKBOn0R05qY4sThUPY796Wfc050A7mk5DEFq3vXYt47LnyczYdq1O7vpRfGMzTPu7fUwJ130y9iBXfGrz47FzeeOBgCcNqI3vjZ5MP7vv0ZKjx95SGeLAMkVtSz9F2y0oymj+xZG/Xnfm1l7MZuKLYNLDZ9N0mRdSFxJRG8T0d1E1E12ABFdTkTziGje1q1bAxckU99mXjoB1585whKFc0IhNcPJHovMqML6ZDz0jU8qM1K6llH4V9gcZNZjnKquqo7545zluJmb5GGT+WZlaBLGyKxDm2rLdq8RkjW6RLJfssFpblJf34/W4YX5XK/Z4mFRmT6G9OzgWOVs7MCu0rxYh5h8Nc2FKDRj1PzS/5yI0f3zqWh01lYwL5PqZcJTCWJZWgyD97fuk56jy7TxA9CpsPxoTXUVrp86Ej10Ne7inAl145AFUJiR+UHMJinzvfeRmNaiWtEyKMGXW4oAInoBQF/JrusB/AHAzcgL25sB/BLAl+0HCiFmApgJABMnTvQ3vPfgjFHOqo0d2NWRufXH547BT595F5MGlyZKyUwhKnp2bCM1E8lOmTy0JzrW1TiuufSmT8m7jgCahGObb3OT9VpXnHiYRYUvhiN6Cgky/S3bLyvbJiR8pFsJp0nkzx3Rr7NUu5KhMlMErcZLBYfnHa+s9jy26KsSQFNhoSrDFDOkV0ccO6QHlm7YLTU3lSaP5v/t0KamGDbqKMf2WykkTKNp+ys7f0J/3D93nec9RYG9dsbn4Gcwbz+0fWFwZF7foygkTPc9un9n/P7iCY7rpT2vLlUhIYQ4Tec4IroTwBMxVycw9T074A+XWF+uzKnqF1nX/LevyBYKKtle7di1Dcs+SZVkAmHn/iZlHeURRlZNwo6xXWuVMLJO8PvBWSMxtLBko9Pc5Myw6hYSazc3hVHq/Tqun/zO8cpEj0n4JMy33mTTJIBSKhZZCKy9/FOP6I0H3vowf93Czl9cMBaPLFjvfC4adbefM+HQ7vjrlyfhC3fPlZYfJ4YG7isE1nbwZycOxPa9jfja5JKWVzQ3iZLZ7jPjB6B7B+ea1X7CuuMgs+YmIjLnlT4PwNK06hIFYd5zyZ7ax3ZN74uaHWSOfZIv1og86t6hDtd++gjLNl0MM5LKTFXMaaMxPCvmfSr8+5XjBxdNfs5j1Ske3OppLysIfs3Gow7pgl6dggUZRNFlGM9FQBSFtjm1R87FzGLfZrRPoNTezp8wAPd97RhPU0wRD+et5Rruu/GpUX08jtDHK8RYhv3Q2uoqfOfUoWhnMsuZjyn6PRTXS9snkaom4cHPiGgc8u9pDYCvp1obn0SZpfG4w3pg9vdPdWRK1aG+RwfMXr1DeyH3v345r6ks+MHpxW0PXH4sxt/8vPR4N5OPlyahE7VRRUCLohxHXeDPz+BcJc27DOW1iuab8BZPzw7JtD9oaeYSDLOgVZOwCnJze5Z1ggZeJjudRxzG7Pf0VZMxpFcHDL/hGV/nqQJRotAk3Mv1Pi/J2dUyMiskhBCXpl2HMETxXs0djluaB7eOYvrZo3DKEb1x5ICuzp2SOsrUXdk2l0s4kqrZaWkxopu8FVkqpKnT05qca3r70SSyEt2URJdQEmqlbebOvjRD2Vkbt/o5fHGKcgHgmf8u+W6E4hhpGS67+3RuK11hLyj2uQ7S+nj8dr++KAqolGWBksyamyqNMKMB9Qgj/6/b4LVtbbXUAQ84OzZZqK8Mu9nLjpcm4WdxHeMetXI3ycxNthb+ly9PKmbpdGgS3kW4lJ2clCDF336Qdf5mTcKe6+j0kaV3btyqrGz7Y7CnvDbvP6JvaTEvr0SDukTdzxaXPfWlSegcY5j73NfzyAIsJGIigL8uccwjx5+dfyReuPpErfPa1FShUxtnhJVBaZarPL7dzzKdbut5y/DySZw4rFdx1Go3TUUR3ZQEQYuqraai/0B2CXNoqz219nVnjjDXwKVu1n1j+nfBcFOqDx0zbJhnGfo12M4vzSHRv7DWeuymAoR9W4GHv3GsdqRcnLCQiInObW2L6Gi0sYU/OB1zrz81pho5MTsqJxzaDX0662Uu/fLx9a69dlGTaFFoEpqT6QCzJqH34Tnj9J3nGSM3h5AK5ZMIfq4dr47UvPcixRKfZoy3cPbY/hhSSPkge551MnNT4Tiz8JW9ir9+eRLOHnuIs65E+fbicq65joA85FPX9xL1im0XTRqIjm1q8GnF+uzyOmgcY1gBIIoDLvsiSxMO7e6+dHJCZNYnUe4Ysyz90M3F9u+G39nfBmaTg+7Uf/scEdknMah7fjLZyYqlQJs10nIU6yiJ01dB5Oysf37+kfjVcyssCxwVr+1j4p1O2UA0kSi6g9aHrjgWE11WRnO7dqmTMu8zOcQLO3Tv54RhvZRRZ155zhzHeApJF00m4mHv4b07YelNn/J1jm6QBZC/7/8+bRi6d6jDuUfFkxAyLKxJxAQR4c1rT8Hxh/fM/47B4FS8ZgRTCIOq+LLTBnRrjwU/OB1fP2GIcydKaxDorhGRL0drfOa4j2F9OuGPl1rnsBhCNcoQWJKMuOPCbz3POrIfzhzTF9dMGV5sM17v254Qz4zfAC5z4IPOdxCFT+KCCQNwky3ZpRtfOLYeQL69+C4zgNnSHDjQrq4aXz/xsNRDXVWwkIiRQ7q2K05yi8NkHeU1dUw/lrI99nfvUKfszK4+fZhl5S7Xcnw6rnWeSU4xSg7zjRoffRS+CfsVjhmS1xYev/I4y37dvrp9XQ1+f/EE9OnctrTAkOnLnyIJbHCLbvKLzOltRycU1Osa5nN/fsFYXFZY/EqHM8f0w5oZU0MnyPziJ+tx9jin2c2OcQ9uix9lBTY3xUwxvC3WMsITdBRDyDvYdh9o9jw20PWLjms9n4T5uIsmDZQfqBASUawnEYe56YIJA3HXZZ9wDDgC9S+SkM7fXzzeEYlmj24KgyW9isbxMlOtbi3SHovbU/WrKGWXTbvG3rCQiBm3XP2e53rsj7J5+e0MzB/+hEP92cX9YPS5uvlrjGr179oOP5nmngE1zHoSqmvZn+OIfp2xfNPu4BdG/t5laVeCTNxzzCAWeW2hzibc7NFNUgI8L+XiSIV7OXl4L8sKkP6vH/jU0Pzly5O0j735nNEY0K09Thou9+NkCTY3ZZAzx1jV/yTafdBMk3HHdpcc1zqmqdKiQ24dqCprbhS2cPtjjCKE0S7MSik1wlzTfb9bWo4wqMo17mXqkYegba1zMpzFHOVy/Th8f14YJfbyYarq0bENrjtzhNaE0rRhTSKD/O6i8Wi5UOC0X73ielzJ7OCvu3ji28c7PkT7pCcvkhqxuU3cchyLUifktgaGKgQ2TIfoZ4KgF/aOzi4gf/nZsfjjK+9j4qHdfF9bZwYx4EzLERVewj5scWlqEmVgOQoEC4mYMc0j1T6nqopQBfK0OZszSfrBvAaxgTn5mIxfXjAWH+9vlNQhXuwJ/tyPlS8Vaac4QUpjxnWntjXFRevdcIsG8ovXRMwB3drjR+c6l83Uoei4JutvO8YUlygc11poCi8guzOTM1qt0LCQiJkwPgmDLDS+z0wYYPltVCnuuvkpRwhz1IjbcfIQWFkHtWS6Xox8TnFNM1N9TMgyY1/kJwrM60nIcHdcRx+RE1X+ojgdwV8/YQgO7dFBuT8NU1cSsJCInfijm8JQV12Fxhb/y0MmNZrzG1pqzHT/0nH1ymOKaTki9EkYs8vdzDO3Xzxe61r2ekS5Mllx0OJxnI65KY4WELaTj7NZft+SmiTZstOEhURChOlUvcxJYcZ1C2483fsgF+IePRl9lO7ja1tbLZkVbkU1mzjMO4rS0fvZiQNx/9wPMah7e6zbsR+1Pv1FbpSim9yPs6flcLtWFOj6SrxIo5+uVOFgkH3XepkT51yZUPHyBTq2qVGuaudadvAi/ZXjY56EX6I0TdhXcgvDUYO6Yc2MqcX08LUx+AUM7cRYWMqO7hKzURP2laTpr6hUWcFCImZUGR51uGbKcFSR+1oSaRO7T8KYJ6FRjltd+nd1TtDSzVelQ5QzlIvXLCZCjNLcVLL9r5kxFV+dLE+dUorWUl8rylfvNdAx+2VcZ1xHVJ8gVKpGweammDF/lH4568hDcNaR3lP87REqn6jvhh37nJFI5UhUZoh/f/t4bNrVYNkW5Sg5yhnKBs0+suXq4phMp2DaUf2xcN1ODOrudNT265IXuN846fAI6+Xuuzt2SA/07FiHbXvd23W6HXVlSgkWEjGjawMOgupD/+cVn4y+MEfZsRcBANiws0G7PLfRaPcOdY4V9qI0Nxlrc3RuF90n1VywYdVG6Lju3j7/DA7rpY7SAYBLjjkUFx99qFQz6tCmxtPvExTVKyEiXHHiYfjRk8vds8CmaW6qTBnBQiIpKjU8LrlJdfKCnvvuCfjyPW9h/ccN0v1uRDlCnzK6L/73jGH44nGDI7umn3U3dBk7sCv+9pVJOHpwD9fjZEvBxkkQv9pVpw7F8L6dMOqQzjjx5/+JvE76VOa3bcBCIma6tMuHZLaJMELFwGia6SSSzIZDc1ifTujeoS6QkIjS3FRdRbjylKGRXQ8wr7sR7bOePDR7+YKCTDr97unD4qhKYCpVVLCQiJkfnjMaY/p3wbGHuY/cgpAF9TYp9d41F1NAIRml/yAO4kqNwcRDVmeCh4WjmyTc+rmx+Mflx0RyrS7tavHVyUMqrgEldTvfm5IP0VQthWrGb52y3vk2FXwSUU6mywpjB1hTw+hkJnA75oi+/hcLippst6bgsCYh4byjBngfxACI/8MwQh9bNNQFvxqFubP52fnuacXTwBCMUfokssCqH39aGTQQ9E4fuPwYrN2+P3ilGCWpDlGI6AIiWkZEOSKaaNv3fSJaRUTvEZG/RWZbGUHWFQhLUrmbDHu8WySS3zrcMHUE6mqqitpd57Y1+OxExQJFKRJHCGwWqKmukkRNebdhYw1teyp9AOjavg5jB3aNoHbBqTBjQZG0NYmlAKYBuMO8kYhGArgQwCgAhwB4gYiGCSFakq9idik6rtMoO6EPwphIFqVl6KuTh+Crk4cU501k1ezUkqtcc5OdkilJ/S6G9+0UW+htFFRqBGOqrU8IsVwI8Z5k1zkAHhBCHBRCfABgFQD9ZZ9aCZXm55BRXPEtho7c+KizuoTkuMLIuG1t5QsJg2y+CT0y2oxCk7YmoaI/gNmm3+sL2xwQ0eUALgeAQYMGxV+zDJJGCGxSo6ZcccZ6fOXFrUnccekEtJOstubF7z4/Hh9s24f2dVn9TKMjlShuRovYWx8RvQDAaUQErhdC/Et1mmSbtB0JIWYCmAkAEydObFVtrVJHLmZK2VXVxwQVkjmRTIjpp0bJmr83HdrUSBeIqkTCpK9Jm7TrfPM5o2JdHCp2ISGEOC3AaesBmD2JAwBsjKZGal78nxOx/yC7PfwQtxaTi3DFNztRribHuPPC1Sdg9dZ9aVcjVuJsRq9dczIm/+xl6b5Lj62Pr2Bk19z0OID7iOhXyDuuhwKYG3ehh/XqGHcRkVJyXKcQ3ZRQv2qsMBdldFPp2jxZLSkO790Jh/dWz2WIM8dZUsRpEh3YvX1s1/Yi7RDY84hoPYBjATxJRM8CgBBiGYAHAbwD4BkA3+LIJgkeS1AmQdxFR7mYjx3Dr9K1fW3k12aCUc4RQuVbc3dS1SSEEI8CeFSx78cAfpxsjcqLEYVZpp+o75542Ul9EKVU4dFfe1CP9rhh6gitdOxMvKQ50AlLpQoHg6yamxgNJtZ3x9zrTkXvzuktShT3RL5cDIv5mFEtusMkS7EVlXGPW86mMjdaTwB2hZKWgEhqjsbIQzoDAI5ymU1bzqNQxko597PlbCpzgzUJJtNMHtoLb1x7inT5UUbOY986LhbzXJykkVomaipVk2AhwYQiiW/bS0BU6scZlHEp5zAKQzlnEagAOSeFzU0MwzAhKGO5pgULCSYQxodRqaMnJlmKCf7SrQYjgc1NTCAqffTEJEu7unxuq6iXaq0kfnb+kRiQgm+OhQTDMKkzY9oYjOzXGccMiX6Z30ohrTVPWEgwoUgjJQhTefTo2AbfPX1Y2tUIRaV+C+yTYAIxqT4/4uvQJv1xxleOHwwAqO/RIeWaVD7t6/ynPK90KnV+hEH6XzhTltwybTS+fuIQ9OzYJu2q4Jxx/XHOOOlyI0yEzL/hNNTW8LiytcFCgglEm5pqDOujzurJVB49MjAgYJKHhwUMwzARUKnh4CwkGIZhQlDp4eAsJBiGYRglLCQYhmEYJSwkGIZhQtCuNh8WXKlmJ45uYhiGCcE9X5qExxdvQN8UF/+KExYSDMMwIRjUoz2uPGVo2tWIDTY3MQzDMEpYSDAMwzBKWEgwDMMwSlIVEkR0AREtI6IcEU00ba8nogYiWlT4/49p1pNhGKa1krbjeimAaQDukOx7XwgxLtnqMAzDMGZSFRJCiOVAeS9+zjAMU8lk2ScxmIgWEtErRDRZdRARXU5E84ho3tatW5OsH8MwTMUTuyZBRC8A6CvZdb0Q4l+K0zYBGCSE2E5EEwA8RkSjhBC77QcKIWYCmAkAEydOrNA8jAzDMOkQu5AQQpwW4JyDAA4W/p5PRO8DGAZgntt58+fP30ZEawNVNE9PANtCnF9utLb7BfieWwt8z/44VLUjbce1FCLqBWCHEKKFiIYAGApgtdd5QoheIcudJ4SY6H1kZdDa7hfge24t8D1HR9ohsOcR0XoAxwJ4koieLew6AcDbRLQYwEMArhBC7EirngzDMK2VtKObHgXwqGT7wwAeTr5GDMMwjJksRzelwcy0K5Awre1+Ab7n1gLfc0SQqNSFWRmGYZjQsCbBMAzDKGEhwTAMwyhhIQGAiKYQ0XtEtIqIrk27PnFDRHcT0RYiWpp2XZKCiAYS0ctEtLyQVPKqtOsUN0TUlojmEtHiwj3flHadkoCIqgvZGp5Iuy5JQURriGhJISGq63wy39du7T4JIqoGsALA6QDWA3gLwEVCiHdSrViMENEJAPYC+KsQYnTa9UkCIuoHoJ8QYgERdQIwH8C5Ff6eCUAHIcReIqoF8DqAq4QQs1OuWqwQ0dUAJgLoLIQ4K+36JAERrQEwUQgR+QRC1iSASQBWCSFWCyEaATwA4JyU6xQrQohXAbSqeSdCiE1CiAWFv/cAWA6gf7q1iheRZ2/hZ23h/4oeFRLRAABTAdyVdl0qBRYS+Y7iQ9Pv9ajwzqO1Q0T1AI4CMCflqsROwfSyCMAWAM8LISr9nn8N4BoAuZTrkTQCwHNENJ+ILo/ywiwkAFme8ooebbVmiKgj8hM1/1uWMLLSEEK0FNZlGQBgEhFVrHmRiM4CsEUIMT/tuqTAcUKI8QA+DeBbBZNyJLCQyGsOA02/BwDYmFJdmBgp2OUfBvB3IcQjadcnSYQQOwH8B8CUdGsSK8cBOLtgn38AwClEdG+6VUoGIcTGwr9bkM9iMSmqa7OQyDuqhxLRYCKqA3AhgMdTrhMTMQUn7p8ALBdC/Crt+iQBEfUioq6Fv9sBOA3Au6lWKkaEEN8XQgwQQtQj/x2/JIS4JOVqxQ4RdSgEY4CIOgA4A/lVPyOh1QsJIUQzgCsBPIu8M/NBIcSydGsVL0R0P4BZAIYT0Xoi+kradUqA4wBcivzo0lg7/cy0KxUz/QC8TERvIz8Yel4I0WrCQlsRfQC8XkiIOhfAk0KIZ6K6eKsPgWUYhmHUtHpNgmEYhlHDQoJhGIZRwkKCYRiGUcJCgmEYhlHCQoJhGIZRwkKCYRiGUcJCgmEkEFEP03yKj4hoQ+HvvUT0+xjKu4eIPiCiK1yOmUxE77SmFO9M+vA8CYbxgIimA9grhPhFjGXcA+AJIcRDHsfVF46r2BxMTLZgTYJhfEBEJxmL2RDRdCL6CxE9V1j0ZRoR/ayw+MszhVxRIKIJRPRKIUPns4W1LbzKuYCIlhYWDHo17vtiGBUsJBgmHIchv37BOQDuBfCyEGIMgAYAUwuC4rcAzhdCTABwN4Afa1z3RgCfEkKMBXB2LDVnGA1q0q4Aw5Q5TwshmohoCYBqAEbOnCUA6gEMBzAawPP5HIOoBrBJ47pvALiHiB4E0Koy1jLZgoUEw4TjIAAIIXJE1CRKTr4c8t8XAVgmhDjWz0WFEFcQ0dHIaymLiGicEGJ7lBVnGB3Y3MQw8fIegF5EdCyQX9OCiEZ5nUREhwkh5gghbgSwDdY1TxgmMViTYJgYEUI0EtH5AG4joi7If3O/BuCVjv7nRDQUeU3kRQCLY60owyjgEFiGyQAcAstkFTY3MUw22AXgZq/JdAD+jbz5iWESgTUJhmEYRglrEgzDMIwSFhIMwzCMEhYSDMMwjBIWEgzDMIyS/w+fbXCnQ9RFvAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# First order system\n", "a = 1\n", "c = 1\n", "sys = ct.ss([[-a]], [[1]], [[c]], 0)\n", "\n", "# Create the time vector that we want to use\n", "Tf = 5\n", "T = np.linspace(0, Tf, 1000)\n", "dt = T[1] - T[0]\n", "\n", "# Create the basis for a white noise signal\n", "Q = np.array([[0.1]])\n", "V = ct.white_noise(T, Q)\n", "\n", "plt.plot(T, V[0])\n", "plt.xlabel('Time [s]')\n", "plt.ylabel('$V$');" ] }, { "cell_type": "markdown", "id": "b4629e2c", "metadata": {}, "source": [ "Note that the magnitude of the signal seems to be much larger than $Q$. This is because we have a Gaussian process $\\implies$ the signal consists of a sequence of \"impulse-like\" functions that have magnitude that increases with the time step $dt$ as $1/\\sqrt{dt}$ (this gives covariance $\\mathbb{E}(V(t_1) V^T(t_2)) = Q \\delta(t_2 - t_1)$." ] }, { "cell_type": "code", "execution_count": 85, "id": "23319dc6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "mean(V) [0.0] = 0.17348786109316244\n", "cov(V) * dt [0.1] = 0.09633133133133133\n" ] } ], "source": [ "# Calculate the sample properties and make sure they match\n", "print(\"mean(V) [0.0] = \", np.mean(V))\n", "print(\"cov(V) * dt [%0.3g] = \" % Q, np.round(np.cov(V), decimals=3) * dt)" ] }, { "cell_type": "markdown", "id": "3196c60d", "metadata": {}, "source": [ "The response of the system to white noise can be computed using the `input_output_response` function:" ] }, { "cell_type": "code", "execution_count": 86, "id": "2bdaaccf", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEGCAYAAAB2EqL0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABRYklEQVR4nO2deZwcZZ3/P9+uPqfnPnKH3AQIkAAhXAoegAFcgxcLKurqLmbXY1d/6uKxri66sKvruazIqqsurogggohyKshNAiQQchDINUkmk5nJnH1W9fP7o+qprqquvma6u3qmv+/XK6/MdNf0PNPH83m+NwkhwDAMwzDF8Hm9AIZhGGZ6wILBMAzDlAQLBsMwDFMSLBgMwzBMSbBgMAzDMCXh93oB1aS7u1ssXrzY62UwDMNMGzZv3jwghOhxu29GC8bixYuxadMmr5fBMAwzbSCiffnuY5cUwzAMUxIsGAzDMExJsGAwDMMwJcGCwTAMw5QECwbDMAxTEiwYDMMwTEmwYDAMwzAlwYLB1B1jiTTu2NwLbr3PMPXFjC7cY6Ynn/7VFty37QhOWdCG42e3eL0chmEM2MJg6o5n9x4DAIzE0x6vhGEYKywYTN0yOJ70egkMw1hgwWDqjqCivy0HxlMer4RhGCssGEzdIaAHu8eTqnmbluEAOMN4DQsGU3ekNV0cYoZgjCbSWP2V+/HNB3Z5uSyGaXhYMJi6YzyhC0UspQEAjowkMJ5U8d2HXsHTrw16uTSGaWhYMJi6IpMRSGkZAMCEIRjWbKkv//ZlT9bFMAwLBlNnpDMZ8+t4KuuSMu/XMjk/wzBMbWDBYOoKGb8AshbGaDwb/I6GuNaUYbyCBYOpK1SLBTGRzLUwIgF+yzKMV/Cnj6krUhbBGJrQ6zBGjRjGnNawzQJhGKa2sGAwdYUUBKKsYIzE0wgHfDh5fisSac3L5TFMQ8OCwdQV0iU1pzWMY7EUhBAYjatoiwQQCiiIs2AwjGfUhWAQ0Xoi2klEu4no2gLXnUlEGhG9q5brY2qHzIKa1RJCWhOYSGkYTaTRGg4gElCQSLFgMIxXeC4YRKQAuBHAJQBOAnAVEZ2U57p/A3BfbVfI1JKUqruk2pqCAPTA99BECu1NAYT8PiRVTqtlGK/wXDAArAOwWwjxmhAiBeBWABtcrvs4gDsA9NdycYydnX1jVa22lhZGeyQAQBeMwyMJzG2LIKD4uA6DYTykHgRjPoADlu97jdtMiGg+gLcDuKnYgxHRNUS0iYg2HT16tKILZYC3fPtR/OXNT1Xt8VWjcK+jSReM8aSKvpEE5raH4fcRVG5CyDCeUQ+CQS63OXeFbwP4RyFEUQe2EOJmIcRaIcTanp6eSqyPqSFOl9S+wRhSWgbz2yPwKz6onFbLMJ5RD2WzvQAWWr5fAOCQ45q1AG4lIgDoBnApEalCiN/UZIUNwnX3vIzF3VFcffYiz9bgdEm90j8OAJjXFsHRsaStdQjDMLWlHgTjWQAriGgJgIMArgTwHusFQogl8msi+gmAe1gsKsuR0QR+9NgeAChJMNJaBgGl8gaqFIyuZt3CeHz3AABgfkcELx0agRD6bAzF52aYMgxTTTx3SQkhVAAfg579tB3AbUKIbUS0kYg2eru6xuHoWHnjUCcsw40qSSKtC8aCjggAYPO+Y5jfHsHK2S2mQKlsZTCMJ9SDhQEhxL0A7nXc5hrgFkJ8sBZrajSOxQqPQ93RN4rlPc3m92MJFe1GnKGSyErunuYwAgohrQm88/T58PkIfsOqUDUB7kHIMLWHP3YMAOCVI3qsoNllJz44HMf6b/8ZHzx3sXlbtSquE6r+uOGAD01BP0biaVOY/NLC4MA3w3iC5y4pxnuSqoav/k4fTNRmBJutSPfTo7uyacrxKlVcS5dUKKBgYWfEtqaAolsYHPhmGG9gwWBw8Fgcsrwh6M99S6SM6uphy+S7ajUBlI8b8vtMF5i0Ovw+tjAYxktYMBj0HosDABZ2RkxxsCJna8vusUD1XFLJtAYiXTA+YLjAzljUAQDwSwuDq70ZxhNYMBhTCOa1RVx7NbmJQyUtjP6xRPZx1QxCfh+ICKcd14G9N1yGE+a0Asi6pLjam2G8gQWDwYQxO7szGkRKzRUCOVvbdluFBOP+bX1Y97WH8OSren+qRFpDOKC4Xpt1SbGFwTBewILBmEHtjmjQNvFOEnMJcMdTldm0f/9SHwDgwFAMADA4njKrvJ2YabVsYTCMJ7BgMJhI6oLQHgnkxDAyGYFP3bYl52cqZWHIgsHmsJ7O+9rABJZ0R12v5bTaqaFlBH616UDVii6ZmQ8LRoPzq00H8OD2I4gEFIQDCjLC7vKZcHFHAZWLYcgWH/LxDo/EMd+o8nbi57TaKfHzp/fhM7dvxS+e2e/1UphpChfuNTifuX0rAKC7OWim1Ka0jHmaz1dvUSnBkIHseFqDEAJjCRWtYXeXVIDTaqfEzr4xAHqVPsNMBrYwGAB6hXfQEAmrW8otfhEO+CpWuCcD2fGUhlhKg5YRaM0Xw5BZUhz0nhSHhvX06b6RRJErGcYdFowGRojsSb2rOZS1MIoIRjTor1gMQ4pAIq2ZJ9+WsLvhm630ZgujXIQQePHgCADg0Ejc49Uw0xUWjAbGWnPRFc26pJI2wci6Lz76xmW46X2nIxxQKl64F09rGEvoleQteVxSnFY7eQYnUhgY1+ttpKXBMOXCMYwGxpot094UQMgSw5BIC+NHH1iLN584GwDwjft3VSyGoRnWwpYDI5jTpge7W/NYGNlKb7YwykW+1m2RAI6MltfKnmEkLBgNzLhFMIJ+X8EYxpy2sHlbJKBULIYhN//Hdg/gMWNYUlELg7Okyka+1nNaw9jVP4ZMRsDHQ6iYMmGXVANjdSsFFF+eGIa+0TQFs2eLSAVdUm6bfzELQ+MYRtnIWpuelhCEyDZ0ZJhyYMFoYKzCELQKhotLKhrMtuuIBEu3MHqPxbB537G897ulyObLkpJptfXuknrlyBhu39xruy2parYkg1qgahl88TcvYv9gzHRJzWoJAXBPZmCYYrBgNDBWwTh7WZerS0oKQ8QiGE1BpeQN5y3fehTv/P4Trl1wAffOs/mypKZLWu2l3/0zPv2rLaZAHBqOY+UX/4AfP763puvYfngMtzy1H3//y+dNl1SPIRjVmmfCzGxYMBoYuYl//72n440rZ7m6pCZcXFJNQX/JgjFhXJdvg3JzL0XyNR+cJmm10gIaMzbp7YdHAQDX3fMy7t/WV7t1GO6+gfGk6VrsbmYLg5k8LBgNTNI4qc82AtpuabXxlIaQ32e28ACAaEjJ2zIkH/l85s7N/y9WzwORezA2MM3Sao8ZbeOtz+c1/7u5Zr9/1Bh4FU9pGLfEMID8LV8YphAsGA1M0hiHKl1RIVMwspt7LKWhKWg/8ZdjYUjyWRjWzf/TFx+Pb16xOu9jZF1S3lgYY4m07bkphpwzks8dV21GDMEIKj4zhsEuKWYqsGA0MDK4HQ7ob4OgoguD0yVldUcBegwjpWbKmnyXz8Kwbv7dzSEElPxvSXmfV80HT/ny/Xjvfz9d8vWycr0ckakEmYzAgy8fyQqGXxeMoN9n9ulilxQzGbgOo4GRwiCFwi1LKu5iYcig9Gg8jS7DJ+5GxuJuyneitW7+0VDht6M5D8PDLKlNBTK+nMjn1yrA0aB7fKaS3P9yHzbe8hzajGyzgOLDREpFc8hvJi/E2CU1bRBC4OEd/Vg5pwULOpo8XQtbGA2MKRiGULgHvXMFY3arHvMoVjFstSry1W1YN/98k/YkyjQboCQtMBnDuPrsRZhIaVWPwcjnR1oYmhCYSOqvYzQkBYMtjOnCtkOj+PBPN+Erv33Z66XUh2AQ0Xoi2klEu4noWpf7NxDRViJ6gYg2EdHrvFjnTEOOY5VCEXIRjHhKtaXUAtmq777Rwj2JrJuSjJc4sWZJFatTICL4feRJ0LvUYsHbnj1gfp1yCMbcdv15i1W4D5cTp6swmc5gPKlbGE0B3YpjwZg+PH9gGADQe8z7HmCeCwYRKQBuBHAJgJMAXEVEJzkuewjAaiHEGgAfAvDDmi5yhpIswcKIpTREHTGMOYaF0TdS2MKwuqHyWRjWzU0poVWFXyFPLIxS4xCfu/NFy8/YBaPFcLlVOwguq7ol8bSG0XgareGAKf7xlIp7XzyM8254GK8cGavqepipccRoR98ZdS9orSWeCwaAdQB2CyFeE0KkANwKYIP1AiHEuMgeP6MApodPos6RG1nYEAq/j0CUW+nttDB6WkLwEdBXpE229RT7wMtH8Mtncye9qRmBt546F59403K8YeWsomv2+3w2N5YQoib++HwWkpVXjozZLBFrDCPo9yHk15/HZNUFw/58xFMaRhMqWiN+BP0++H2EiZSGW589gIPDcbx6dLyq62GmxqjRxbkerMJ6EIz5AA5Yvu81brNBRG8noh0AfgfdymCmSDytIaj4zOl6RISg4svpJeW0MAKKD93NIfSNFh7EY93I73z+IP7xjhdzrklrGXQ3h/Cpi1eWYWFk1/ejx/bgpC/dh4Hx6nZgLaX30k+e2Gv7PmVaGBpCig8hIxstWYJLSgiB2zYdmFRX4AnHxpJQsxYGkG3tItdRbQFjpoaMRcWSLBgA4LZL5FgQQog7hRAnALgcwHV5H4zoGiPOseno0aOVW+UMJJ7SzJRaSdDvs8/DSOZaGAAwty2Mw0Umt5WS669qwhyMVAp+n90lJXs2vdpf3VNywsXC6B9L2E7zS7qjAIAfXH0GAHsMIxTwubaPz8fvX+rDZ2/fiv/64+6y1+q0MITQq71lj65IQEEirZmvs1d1IkxpyALMeii2rAfB6AWw0PL9AgCH8l0shHgUwDIi6s5z/81CiLVCiLU9PT2VXekMI+7ibgr5feaGlskIjKdU1+6xs1vDOFLUwsgVDGdgW8sI08IpBd0lld3gZG3G/qFYyY8xGawnfVXL4Pn9x7Duaw9h1T/fZ448lS1BzluuvzXlRpxIaQgHlGwlfQnurX2D+t/jtBZKwc1Fl1QzaDZiKJGg3m1YCka9N3NsdEaNeh52Sek8C2AFES0hoiCAKwHcbb2AiJaT0S+CiE4HEAQwWPOVzjDiaS2nb1NQ8eH/nt6PbYdGEEtrEAJodhGMuW3horOh3bKBnBtgOpMx6ytKQfG5B72HY+mSH2MyWAUjoWZsKcVSOOU1TQEFPsoKhnyey4lhjCf1v6c1z2yQwj+bXevs1mydjKxzkfNMZJZcilud1zXSwqjU0LKp4LlgCCFUAB8DcB+A7QBuE0JsI6KNRLTRuOydAF4iohegZ1T9pah1r+gZSDyt5dQ+yNP+W7/3WMGRqbPbwhhNqAUDzgmXE5F88wO6dSFEdjBSKQQUsgW9M8bbQAYGq4V1k9fdOdm/TboKEqqGoN8Hn48QtFhqibRuyQVdWq/kw2f00zoyVliU3YhZXFIyow2AWU8TDihIqBm2MOqQTEZgt8O9KmMY9eA69FwwAEAIca8Q4nghxDIhxNeM224SQtxkfP1vQohVQog1QohzhBCPebvimYHcyKzILB8hgHHDFG52qcCeK2sxClgZbmKytXfETKWV//vLiGEoPrJlIkl//Ui8dhaG3sxPtX0P6K4mGaewJg9IYXarc8mHFJv/e3o/DpTpbrP6ujuiQfPriCkYPiRSlhjGNGnm2Aj84NHXcOE3HzE7HAPZw5CaEbbuCV5QF4LBeEMsleuSsp6kZXtut/kUs1tLEAwXE3rjLZvxwf95BkII07VUTtA7oPhstRty4x6tumBkf2dS1WyBZelmS1gstqBfMZ/LeDpTtksqrWY3hpctm0cpWOsw2i3DqGS2m5yYmOAsqYpx86OvYvO+oSk/jnwMeUhIqhoS6YzZUsZrcWfBaGAmjOpfK9bNWDbPcxOMuW0RAHr16Zfv3oY9AxM518RTGqydyv9i9Tws64ni8d2D2HlkzAxel+OS8it2C0OucajKMQyrGymeytjiBPGUbDKYMbPOQn6fLegdCSjm8zhYQgqw9XX43dbDZa11Iqmacy+stS3SJSWD3lL0ymkiyeQyEkvjX+/dgXd+/0ns6CtP3J3I1v7yLT5ivK9ll2GvxZ0Fo4EZS6g5AW3b6d0UjNwYhvSN3/rsfvzkib34xn07bfcLIbCjbwxNAQWXnTIXAPC9q07DTz+0DgCwae8x08IozyXlM2dopCx++EPD1W2bYM1sSqiaGd8BstkribSGsD/byFGeBuOG629uWxid0SBeOlh8U5H1KQBw95ZD2OsiyPk4Fkth/cmzsfeGy8xUXyDrkmoNB7C7f9zclOrBNz6dsVqA67/957JdiFYUkv3S9NdEZv8t62kG4P1rxYLRwEykVLNdhcQqGAeH9TerWwwjEtRPzM/tHwaQPQFJ7nz+IB54+QgmUhq+e9Vp2HHdegDA/PYI2iIBbDs0agavywp6+wia8WGSJ2S/j3BgKFZyv6fJYGukmNJweDiBRV1659Cv/PZlaBlhd0kpPjP7SI9h+EBEmNUSwlAsVfT3pdQMIkEflvboG760pIqR1jI4FkubYmNNapAuKWdFPVsYU0POPZEcnMLhRX4U5OstBWP5LF0wat0q3wkLRoMihMC4q4WR3XT/9d4dAPLP2LZWgEeCCjIZYX5YXjAapgF6oFpuXESEk+a24uVDI5MOess1yvjFucu7kVQzePHgSMmPUy62tNq0ht7hmO303jeaQCKddUkFHS4p+fdHgkpJ6ZEpLYOA4sNXN5wMoPR25HLKnxQMa4xKWhjrT55j/11sYUwJ5wFgKvE06ZL63K9fxOfvfBGfum0LAGB+h+4C9vq1YsFoUJJqBmpGoDlUPM/f2RpEYq0ST6sZ3LXlIM674WE8+eqga/m+ZOWcFrzSP24KRrlBb2lJyOyR85Z1AUBOOmIlsQa942kNR0aTmNMaxhcuPRGAPo0voVqD3g6XlBQMowaiGGktg6DiQ1OovO6yu47oz4HcYMLB7GskW5s78XoTmu4cc1gYA+PFLci8WIzk/3s623utK6ofADjozXiCTEPNZz1I2iIB+PIU1lndHbG0hj0Duvn88I4j2H5Y74D6nSvX5Pzcku4oYinNzLBSynBJKZb25jLGIs31YpXnU8GZVjs0kUJXcxAnzWsFoAcnE+mMmQkVVHxIqwJpTRdmKRhNRsC5GGlNIOj3mYHqUgXj4R39CPp9OGtJJwD7ayRbmzvxehOa7hxzWBhTeR/mey3KScmuJiwYDYps1tfdHHS9X8YtZL2FG7IQDYCtmd1YQsWu/jG896zjsGFNTh9J0/e/y2irHfKXWbiXsbukuptDaIsEcNjRPXcknsbia3+Hu144WPLj58OaRntkNAktI9AZDZmV2KMJFcm0ZnNJJbWMKQ4RS9FcKRZGStVdUlIwSu0j9MzeQZy5uMMcq+vmknL7XczkcVoYu6fQ/TdfFpTb6AEvYMFoUGSgrjNqD1Zf/45T8OYTZplxhUKCYY13xFIqjo7pInRwOI6RAuNbF3fpvv+dhmAUm7RnRY9hGBaGsYk3h/3oigZz2oPIzKIf/nlPyY+fj/GkhvYmXRxkMkB3cxCtkey4Wnsdhh7DkNXuYatLquQYBpkbfykiI4TAvoGYmVEDwDYjPegizB1NAQ56TxFnSvfvth6edIFdvk7G5Vqa1YIFo0E5PCyHstgtjKvWHYcfffBMc/M9a2lX3sew9iAaiafNDXzv4ASEgDlT2smCjgj8PsKOvvItjJBfMU9ZZp1IyI+WsD8nk0h+ZKn0EEleJpIquoznSk4+64wGzb9xJJ5GQrVUevv1LCnTwgjYayCKkdbKtzBGEyrGkioWljD3+TtXrsHK2S2Y3xGxuUFueWofrr1ja9GfZ7IMu2S9TdbNl+/nZHLKeLI0S7NasGA0IH0jCfz6+V7MawtjcVfhzWWV4aN3Q765g34f+kYSpjl9YEjfUNvzCIZf8WFBRwS7+sq3MMIBHxLpDJKqZrMwWsKBnH5SMjheAb3QU5DDATSH/GaRYlc0ZLru+kYTGJpImX9L2K+7npwuqVKD3vFUtp2Ij0qzMMYLVOY72bBmPu775PloCQVsVeVf/M1LuNUyZpYpztBEynSzSiYrGG6djFfMajZrocZLTK+uFsXfWcyM4+zrHwIAvOO0+UVbi7sV7UnkSf+EOS3Y0TdmtguR5LMwAGBRVxR7jRbe5VoYfaMJrPrSfZjVEoLiI0QCCloj/pyBTtLVQhUwMeRM7O7moLnuruYg/IoPzSE/bn70NQAwWzjMaQuhfyxptumIWFxVpbiA4mkNTUEFRIRo0J8zdtWNhEOcJPd8/HU5r40k4Pch5pIG6tb6nnHn2EQKpy/qMFvSA8bGn9+bmxdrncX89gg2XrAUl5wy13QtVrvJZjHYwmhgls9uLnpNodPqFWv1MSYnzW1FSs3kuISkz98Nq2VTroUB6I3YDo0k0Bzyg4jQEgrYqq+B7BzxMrqn52UiqSIaUsy4jI+y7jyrMJ52XAcAYH57E9SMwL7BCWPd+t8YUHzICBQtMoyltGzgOqggni5+sow74iWSk+e35RRWSqxNEq1rKpTpM5pI45k9U++bNBNIqhomUhoWdjbl3D65x9Nfi0+8eQUev/ZNuPqcxehuzlqy7JJiaoo1PXR+e6To9YUE4/9dvBLPfOHNZpB1LGnfsItZGJJyLQwr8kTvFsOQAefn9g/jC3fmjocth4mkhmjIjzlGEkBGZAPK8jlavaANbzxBr6KWdRCyNkSe1uXPFLMy4inNjF9EQ6VZGHKTKkeAg/5sEsFrluyeQhvTX3zvMVzxgyc9z9ipB2Ss7ziHYEz2uUmpGVy1biE+ddHxttsVH6EzGnTt2VZLWDAaiJ8/vQ8n/NMfzO8XWzZtJ1952yoAQHvEPe1WMqslbGbfODfstkIWRvfULAyJ9BW3hAOIpTTbRmxt5/FzSxHUZJAuqb+9YBkAYJ4le0yOPj17WTZBQIqxKRimhUG2dbshhEAspWabBQaUkrJj4qmM7XeVgtXCsE4tLHRClq4Xp0XXaCTSGs76V929u7CjCVetOw7XnL8UwOSbBCbVTM6hSHLe8m5s2ntscoutECwYDcQX7nzJ/PoD5yzCqQva8l77gXMXY+8Nl7mmYjqRp+bxhGqzWqphYch4xOtX6GNQ5RQ5md5qDQrKDXSqCCEMl5Qfq+a14jtXrsFtG88x75cbtLUv14KOCBQf4cnXBm3XyOczXWBDSaoZZETWKomGFAxOFO9wK61Hp6gWYkffGPYPxbCjb9RMiwbyj5G1zh3x2j3iNVa3XWc0iOvfcQrOXqoXTE7GwhBCIJ7SEMrz+s1tC2NgPJkz5riWsGA0CM4smzedOLsiwWAguwmqGYFls7JxkXwnJcBuwpdzIpZdaS84vgdfe/vJuPE9pwPIBuetQcFS0ldLIdtGRY+XbFgzHwssqatyxGzUIhjhgILXLe82rS65iUtxdRszK3lw+xEA+ukfAE5d0I7n9w/nDXgOjCfxru8/gdcG7NZMKcg6m1ue2mcXjDwb3mn/cr/5dakNEWcCWkbggZeP2OorrE0HZTwrqJQ+88TJcCyNlJbB7Bb3aHlXNIikmpnUnPdKwYLRIAw4ZjC0lpB6WSrWXlDLeqJ4+vNvxp7rLy3yMz6cMr8NZy/tzNt6xA3ppjlpXivee9YinDxft5JkLONbD+wyr63UDGRZ5R3NkzUk1+/s6ruwM2tthY2fleJS6AS67ZDeLvsSoy38iXP11ObhCXfB+OWzB7Bp3zHc9IieqVWOi+/G956OWS0hPP3aEPptguH+3Fl1rpEE41sP7MLf/GwTHtl11LzN2hKkw3C/SutgMhbGYaNVTr5iWZlwUco8lWrBabUNwlHHm6yQu6hcrC6lpqCSN4XTyW8//rqyzetPXnQ8zljUiXOXddtuP2ORnp00anNJVUYwZDO5jqh7PGd2q/5Bdv4l1pNijkuqQAxj78AElvZETfeetTjQDVk4FjYeuxzBaAr6cfGq2bjlqf0YGE/C79Nbr7idkJ2vVSO5pKRr0Wq1DoxZBMN4b8jPwmSypPpGdet5Th7BmNeu3957LG5z6dYStjAahEFHB81KCoY1ztHRVDhI7qRct1hT0J/TnhsAZrWGsaQ7anMLVcolJdum5ytivPaSE/E3r1+Ci0+abbu909KnS7qisllS+YVyOJY2q8qBrDWYzyUl/055+i+3fuKz608AAByLpbHAyO5yi2E4RaSRgt7ShWhtXd57LAYfAa987RLzdQ2agjF5CyOfYEiRsNZ71BoWjAbB6p8GChfklYu1X1GgSCFgNWkN+22ncKtLSini9hJC4PbNva6n5v4x/YM8v929Kr455McXLjsJ7Q6xbHLZuEtJq9UzpLLCJ7PN8s1ZkLEiGRgPl5FEAOjBevn0yNiM2wlZisjfvUHPFGskCyPjaKkP6Fllc9sitve8fC2KuaRueWof/vfJvbbb+kYS8BHQk6cH29xWPSNR1vZ4AQtGA7C7fwyfd9QhlJL9VCpBywfGy0Z2rZGAbVONlyEYz+0fxqd/tQX/fNe2nPvGEyr8Pior+wgAIi7txEtJq51IabbZFTI2MpZng5aeorQmEFCoaPW+EyIyf4e0MOQY2bu3HMLnfv0iVC1jpinL4UyNFMOQAjoaz/7NQ7E0uh0FkcEiLqmkqkEIgS/+5iX8k+O9dngkgVkt4byvn89HWNgRYQuDqS6vHKneYCFAby8hOW95d4Erq0trJIBjsZTZ5jyR1rB8VjP+4cIVSKmZgtXVUmicyQFAdvZ5+e6zXAtDimuhtNpYUrUNrZIxiXydTOXIWuu15SID91IwfrnpAHb3j+E/7t+JXzyzH1sPjpgWRmskgIBCDSUYx4wCPauFMZZI54w4Lja3YsN/Po4Lvv4n1/v6RhJ53VGShZ1NUxoBO1VYMBqAck+c5WK1MGRGjxdEgwr2DcZwzvUPYzyp6v2QAoq5+RYac5oqMP1vPKmW1NDPiatglODj1i2MXMHIF5NJW4RwsoIhEwTmWepofre1DxnDfNnZN2ZaGOGADy3hAMaTjRHDEEKYMy9sdSiJ3PdFqMDre3A4bta9SA4MxfClu15CLKXi0Egcc4okjPQ0h1wPNbWCs6QaAGke/+4Tr0NaE2jOM6pzslQygD4VrH7/RFozR6NGLLME8sVuVCMI7XeZ/jeWUEsaZevELfgsf3+h07m1yhvIxiQSeYrpNEsAvZwaDCtyg+tpCWFBRwS9x+L41oO7zA3xTzv7zaB/2K+gKaggVkK7kpnAWFI1g953vXAIX3v7KWgO+Y33hX0LLXQgkBMmrfzosT342ZP7EA35sX8whresyk3osNLdEsLhkYRt9kotqQsLg4jWE9FOItpNRNe63P9eItpq/HuCiFZ7sc7pijSPm0N+rFnYjuWzWir6+M6ZGl5h9fvrgpFBOKiYtxdqr5HS9Pv8LhaGm+uhFJpcZqEXS5FNqRmkNWGzMPyKD34f5a0rSdtcUpP7SEsrsaclhAc/dYEpFFLYnnx10Hz+wgHFnCjYCDjrXx7fPQBAWp72g4R8Ht0Ewy2rTNZtPLy9H2pGYOXswp9NGRC/8/mpT5GcDJ4LBhEpAG4EcAmAkwBcRUQnOS7bA+ACIcSpAK4DcHNtVzm9kW/eSga6rbi5XrzAbmHo0+4iAZ8ZfL7h99vz/qxs7ueW5TVZl5TsEGvt2mtO6MuTkiqLBJ3PaSSg5LUwVIuF4czUKpUl3XrKZndzCOGAgi1fuhjnWnpjjSZUMzYUDvhsPahmOs7hVaom9HYxKdV2SAH0BIKg3+ca9HbLKpPp7rI1f75aH8l7zjoOQHXn1xfCc8EAsA7AbiHEa0KIFIBbAWywXiCEeEIIIbtuPQVgQY3XOK2RH+xCrTqmggwGL+vxpphIYq3ETqQ1JFTdJSUrZ+/bdgQHhtwzTOTp3S2uLYPe5dIc8uP//vos3PLhs8zbIgEFAYXyWhhyc4o6g6kBxdZM0Yo1mJ+vjXkxfvxXZ+K6DavMDCifj7DYEBGZYdY/aqTtGhZGo4x2lZ+ft5+mz6fvH0tAzQgIYY/fSUJ+dzF1G34k4xHy/VDs8BUOKOhuDuLIqDdxjHoQjPkArCO+eo3b8vFhAL/PdycRXUNEm4ho09GjR/NdNmN59eg4XjFmZUvkaadaFgYAPPuFC3H3x15XtccvhWaLeyCR1sypdasXtuOz61cCAJ7b797tUwa93QrPJ2thAMC5y7vN9iWALq4t4dzZHRLp9ok63FnhgM9s1+7EunGX08jRyvz2CK4+Z7Httk7DWplliJBsHdJoFoZ8b2xYMw+AvrmnCljtIb/i6pJyszCc9VGlWOtd0ZBn7UHqQTDcchVd8x+J6I3QBeMf8z2YEOJmIcRaIcTanp6eCi1xeiCEwJv/4xFc9K1Hbbebb+4qZkv1tIRyTsW15pKT5+AkI0srkc4gbgkMfui8JQCQ18KQbh231NvxSQa986GPaXXfbE2XVCjXJZUvS0rLCBDp/Yz+4tR5FVunjLfI11UKRsivGDPLG0MwZAq0HJkbT2kFBaM5pOQkNezsG8NXf5frEnWmvDsPCm6EAt5Zd/UgGL0AFlq+XwDgkPMiIjoVwA8BbBBCDNZobdMKaxdLa652Us2AyD1ldCYRDflx/TtOAWC4pNLZMaPSjZIvO0l+ANWM/YOYSGtIaZlJWxhuRIJK3gB2PgujMxrMae8iSWeE3tH2SxebA5wqgXTPyTYl/YbfXD6Xk51bPd1IWmbXNwX12STWefZOuppDGHK0o5cdiIGsEM9uDeU8h86DghsBxVewtUw1qQfBeBbACiJaQkRBAFcCuNt6AREdB+DXAK4WQuxyeQwG9sKuRy1dNVNqBkHFV7F25vWMFIiJlIq0Jmxppq3hgK05oRX5AUyp7g32KikY+qbjvo58Qe/ZrWEcGXMPdGqZjNkFt5JIcZWt6OV0uVDAh0AjuaQsFnpT0K8LRgGrvauAuAPAff9wPh781PmuEy9LsTC8dAd6LhhCCBXAxwDcB2A7gNuEENuIaCMRbTQu+xKALgD/RUQvENEmj5Zb11hPK0++qhthQggci6Um7dueboSNwL4M0FrbuLdG/Hmzk9Q8FobcNJ359lNaYwH3UtLi/rAyqyWUNzMmrYmqCMYHzl2My9fMw0cu0KfIyecu7G8sCyOdY2Go+MNLfeZtTrqaQ2aHY4n1MDenLYzls1pcM9pKqaMJePjc10XhnhDiXgD3Om67yfL1XwP461qva7ph7TAqg2lfvnsbbtvUa/brn+nIOgTZMLDVUlTYEg7kNPC74Ot/xDlLu0xBVR2m/ngVBKMpqJiVw06SZkabfSNqbwogkc4gqWo52W5aRky6/qIQndEgvn3laeZaxxIqfIZrM9SwFobukvravXo8wu0g1t0cxNBEEpmMMFuuyAPCZcaMEwBodyl4LWU2TFChho5hMBVgOJbCP96x1fxepun99Ml9ANyDuTORkHFCkwHaVkvmVGvYb4thqFoG+wZjuPXZA2Z7DecHccxof1HJ7r6FAtgyo805plP+frfUTFXLuFaoVwp5ih5PqggHFLPWoOEEw+9DJKiYbkN5m5OuaBAZAQw7GmF2RYO48b2nm7fJLsTlWv96DIMFg5kCX/zNS3h6zxCA3DbfQJ60sxmIPGnv7tezT6wWRmskYHNJHRrOunhkJoxzdKoUmEoHvfNVnUsr0WlFmB1r3QQjI6qa0GDdFKWrzMtNq9Zk+4z50B4J2npByZGsVuRkPGvPp3gqk+NmXGi0kp+cYDRu0JupADv6srUXs1rDeWcnzHRkEFKOObUWsrWG/Y721Fm3kBQK1bEJmuNZK+iSypenD+R3STlbdVhRNVG0fftU8PvIzJiSfa0a1cLoag7aRtm65ZHIQlFrCnc8reb0Flu9UK/POb5IOxAnXiYcsGDMEGQGC6D7UMeSqm1gfaOYGM5MsMVd2aFHepZU9nn6/YuHza/lB9B5csumuVauSj4c8OVtVW66pHIEQzYtzD0IJNTqNqIjIlOI5e8pJeg9MJ7MEeDpiPw7Q36faT1IZAaZlRWGAOyy1FjEUlpO5tsZizpx0/tOx9ffXV5rvKCfcHA4jm2HRsr6uUrAgjFDsFoUshlgrEIjSqcrqxe02QSkNRJASs0gkdbw6tFx/ODR18z7ZPqsM0tKtv0ud+xpIUJ+BYkCFkbQn5sCLX+/W3uQZDpT9Sw46ZaSMaKg4RbJ5ImNxVIq1n71QVx3z8tVXVctkIeJgOLDPMu8isevfRMWughGWySASECx1WIcm0i5ji9ef/Jcs49XqUjxvuy7j5X1c5WABWOGYD3tScGwBucaxMCw8YtrzrZ9H7W0OXemqJqC4bAwZG8nt86zk0X2GhIufUjybf4ybdbNd11tCwPIWjwyRiQFJJ+VMWRkVv3a0lV1cDyJP+7oL/h7nnh1ALc+s3/K660kKTUDxUdQfISLjLnts1pCNvFwEgr4bG7HgfEUuprzNxb86YfW4baPnFPSerwcg1wXabXM1HCa/bIH0EQDzVx2w7nJy+/Pvv4hcy61RGYfpV0sjJDfV9EYgcyASqq5gVC3tFkgu0k4Ba33WAzDsXTVBUOeaqVwmJMDtdy/AcgKhlyvEAKfuX0rHt7Rj6c//2bMzjMo6D3//TQA4Mp1x1X2D5gCaS1j/r1dzSH8+bNvRCSoFCyEDfuz1fxCCAxOJM3Gjm5ccHzpbYy8PPwVFQwiep8Q4pZaLIaZHIOOnP4O08LIui/cTrMzlQc+eb5rkDpsWBgpNWMWXklkbMO5Ibv5nqeKLC50E4x4WnO3MIwsKKfL7O9+/hwAfR50NQn4c2MYQP5RpFIwtIzATY+8iv9+9DWzzfvW3hFcdFLhyXJCiLrpTJAy3IQSNzeUk1DAZ7ajj6U0JNIZs8XKVKlGkWaplGLbXE1E3zHmVjB1iHPGr3RJDcezQlLNLJp6Y8XsFtuoUUmTZXO2ZpUB2RRIp8tnIqVW1B0FWCwMlxhTPKXlzFgAgIBPnujt65Pp09U+EMgKZCl2xVxSx4wMNCJ9qtzgRAqvHp0AANz1wsGi683XwsULUlqm7E7PYb9iJjDI95YzYD5ZvJi0JynlWVgPIA7gYSKqXGczpmJsPzxq+14G14YslscN7zy1pmuqR9wsBdndVm7EbkHvSlsYIYuF4WQipbkKlLQwnLUPcuZ2tdMsm8wmjvqWIV1k+X6v7KWk+AjOs8o9Ww/j7i05/UVt5KuE94KUKsru9By2WBiyTUihGEY5eDmwrOizIHSuBfAdAI8a8ybWEVFxu4ypCb95/qAtACctDNke5LrLT8allpYEjYoz0ykc8OHaS06w3VYLl5R0OblNZYslc6e4AdYYhkMwjG/zZV1VCunic7qk8hXv7R3UrQkh3OekO3stOXG6Wb1kMhZGyGJhyNkV3dHKWBgb1ujjglorWExaKiU9C0T0Vui9nFIATgfwDQAHiGh3FdfWkLxwYBi3PLXP9b59gxNY+9UHcmY6vNI/bmtrLX3F0q/d5tKzphFxfuiFyDXvnRtgPKVVNKUWyP5O58jV/rEENu07hoGx3M0yoLhnScmOpxcb2TvVosnSJh4oPLv63hcP45an9EyneFpzdT8152njLf/O+rIwtLItDGsMQ4pfd0tlLIw5bWFceOIsLOio/Zm96LNARK8B+FsA3xJCnCqE2CiEOF8I0QXgDdVeYKNx+Y2P44u/eck1v/2O5w5iYDyF2zf3mreNxNMYjqWxyFqgZgjEISO2wYKhM69N31xft7wbgL7ZhWxtL3LbXUyk1JJaTpdDPgvjni16IeFOx8REAPAbG9bB4Tj2Gad3QD8cLOpqwnvPqm5WkXSTSesnVCDo7TzwuMUj8vW+ki1QhupKMDII+MuLAYYD2SwpaWF0VijoDdgtmFpSimxeKoS4TAjxgPMOIUSv2w8wU8dt9oGsI7COGR02goudFnNXBncPsYVhoyMaxN4bLsMHzl1s3mZt8tfTEkJG2N0+1bAwTMFwWBiF+kHJzJgfPbYHb/qPR8zbJ1IquqLBqmcUyaQJ6ZqSz5vTSgKy4iaxjiaVdQz55pPL2enWti1ek9bKj2G0WBpdDoyn0BL2u6ZLT5aQRy3OS4lh7KjFQhg7buMcpVvgz68M4GWjV5LsjWT1Z/oVH8IBHw4eYwvDjc5o9vmwfoh7jCyWhJrB4ZE47ny+tzoxjIB70Puo4de/66Pn5fyMtVjL2nl4NK7aGixWCylH0tpqCeVvVZLP3QQA5y7rApCtoHcis8HqzcIoN4bR0xLC0bEkhBAYGC9cgzEZgn5fzoED0AenvXBguKK/ywpXetcBsrOqld9tPZwT4LSe3KQFIj+wzvbbzSG/JZ2vcqbwTMA6uCbssDAA4LFXjuJj//c8PvnLLegbTVQ+rTaPS0qfHe7H6oXtOT+TLy16JJ6uyYFAWhbSZSebIY67FIcWSvuU6c75RtTKU3M9CYZ1Nnyp9DTr41dHEyoGx1MVq8GQ5LMw3v/jZ3D5jY9X9HdZYcHwmM37hnDhNx/BTx7fk3Pffkdw25q3P2I0GxzN0367y3BRRQIKWirYaXUm0GkRDJuFYQjGxlues1XJV9L3DOQGvf/wUh8ODMWQVLWSByFJsRmJp10H8VQaZ3t1+X7rG01g094h27Uy/va/H16HhZ32epjZrWH4fZR3HogUpHoKek8k1bK7Fcv30tGxJAYnkhU/tOWzMKoNC4bHyPTCO184lJNN4pxpkbD1ptGtB2lhtEUC+NmH1uH6d5wCIPuGndUaqpuK2XpBnsg//LoltqB3T3M2NXlxV7YhXDXTaocmUth4y2Z85H83I5HOlOznHomnkckIjCZqY2FIV5JsyS0t2n//w06866YnzQMMAKQzAkt7onj9ih6cflyH7XG6okF9gFTKfbOTQfStB0fqptPteFJFc5lWpnRvfvnubdh1ZLziLqmQX8mxMGrxfPHR02PMCWpC5AQQndklVjNeugKsFsb5ln40sj3IqnmtFV/zdMfnI+z66iXw+8gsfAMcszMi2Y9GpdMXs4KRMVtUD8dSSKhazqS9fIzE0vARQYjsa11Nzl3ejc1fvNCsVnb69A+Pxs0JcqqWMWMRCx3P3bz2CMJBJW/QO60J+Eg/mQ/FUpjVUriFSLXYNziBC77+J9z10fOmZGE8tnsAQGVH/AL6869lhD5t0XBVH4tVfwYOWxgeY53PIDujvv+cRQBcLIy0HnwL+X1m0FBaGM435LvPWICFnRFcc769yR6jE/T74PORLS7UbXEbjMZVLOmO4jcfPQ9vWVXZGgcz6J3OmK6vjNC/D5doYQzH02Y9jttMhmpQqLVF/2i2lbeWEWZlujzErF7QhstOnQvFR3oVdJ6gd0rNmJ0KnEWUteRho6vubZsOYCKlmdlbpVKpNiD5CLm0ZnGLJ1UatjBqxOqv3I+/Om8x/uHC4223yyynLb0j+PMrRwEAc416AefUvERaQ9jonBozBUNFU1DJSWU8//ge/Pmzb6rK3zJTsZ4ih2IpRAIK1rgEoKeKnFoXT2umLz+e1sqKYQzH0mYcw4sCLgC46X2nY+MtevPDo5YpdGlNmO/HdUs68fw/XWSzgvLNNBdCIKVl0BRSMDiRv4q8FshENJnJVijzy41KWxROgpbUbBmSq0V3arYwakBS1TAST+PbD76Sc994MisKn/zlFgDAws4IIgEFOx0N8pKqhlBAQSSg4H+f2of/fvQ1jMbTaA1z2mwlsIZ6hiZSFa+/kPgVH5qCCkbjadOXP5ZII5F2b23uxnAsZbogOqLevP5nLOo0v44Z1rEQAmomY+uo6nSZ5RMMOSZXpu56NbcayDZzlEWy5X7GAgrZMtvefvr8yi0O7v3I8iUSVBIWjBpgNdeduL3Ic9vCWL2wLWcEYyKdQTjgMwPlX7t3OzbtO5aTIcWUxxcuPRFXn70Iaxd1YmmPHuze3T9e1Q6/7ZEAhuNp8/XPCN3aLNXCGImn8dsX9AZ+Xh0YelpC+NZf6uNFJwyL94Kv/wmP7x4s2II7HFBc6zBkwFsKtbcWhl2synUxEZFpSX78TctxwpzKxhLd2svH8rj5KgkLRg3oN2om3NJb3bJFZrWE0RkN5sQwkqqGsCM7Ys/ABJo4bXZK/M35S3Hd5Scj6PfhMxevNG/fPxgr8FNTozUSwHAsjXgq60bYeWSs5Hz/oYkUnjHSWb1sd71h9XwQ6U0TgWwqeKGpcJGg+4haKRDSwvAyhuH8XHZOwoqTbq1qZLHJGMYVP3jStO5ijeKSIqL1RLSTiHYT0bUu959ARE8SUZKIPu3FGqeCzH7yuZy6EqqGZT1R/Pf715q39bSE0BYJYiTuzJJyn26252hu4R8zOcIWN1Q1Lbf2pgCGY6kcC7MU3zdRthus1/h8hKaAgomUZsviK2Sdhf2Ka9BbHoRkGrMXrS8kua9L+Zu+bJVfDQtQCkbfaALP7RsGYLcwqjUfxXPBMAYz3QjgEgAnAbiKiE5yXDYE4BPQu+ROO+Qb362hYNKoIrUWOIUDCtoiAYzEU7YXPpHWg6LOvjb/xrMuKkbEIsjVLF+Z1RJG32gC8VTGVudRSpuPRZ1N2HVEPyR4ORtB0hTyI5bSzGaXQO5cESuRYDaGcdcLB7GjT29zI90rMvnAyzoMq6vn8jXzsHJOS9mPIWMwlS78BIAzF3dilpG6u7t/DHe9cBAvHsy6sN26CFcCzwUDwDoAu4UQrwkhUgBuBbDBeoEQol8I8SyA6icaV4G08eI5/aKAfpKJBBQs7W7G21bPwx1/qw+Cb4sEkNYEfrVJD7qNJtJ44tVBxNPZiWzRoIJrLzkBl/Csi4phFYxqDrFb1NWEQ8NxjCbSiIb8pk+60GlUtjLviAbNdjI//dC66i2yRKJBBbGUit5jFsEo4E4KKj6k1Ax294/h7299AV++exuA7AbbZMYwvHNJpbTsaf3tpy+Y0mPNaat8LUlHNIinP/9mBBTClt4R/P2tL+AnT+w175/JgjEfwAHL973GbZPCGPC0iYg2HT16dMqLqwQy+0Nz2YGkmyno9+G7V51mZp5Iv+dn79iKeErDjQ/ro0deOjiKT16kp+Zu+uJF2HgB11lUEmtm1OcuPaHAlVNj1bw2ZITeM6wpqJiWpLVg0Mk9H38d/vTpN2CuZQOa1VLdfP9SaAr6cdcLh/CF37xo3qa5WNOSUMCHpKrh4LAe25Mjhp0WRrqAlVJN0loGr/ZnXX6TTZFdbIwcqIZgAHpgPa0J3Pn8wZz73Mb/VoJ6iJa6Gf6TPloIIW4GcDMArF271rsjigUZzHP7EMVTGjqack+V7Zbbdh4Zs7X3eP85i/H+cxZXfqGMzcJ40wnVG0p04Ymz9AJMw8KUh4pCFkZHNIiOaBBzWrPuS68qoa1Ii+DAUNbCSBcSDL9uYcjW/PI5TztiGOkqTxHMx/X37sCTrw2a3/dMsgjvlx85B4/sPFrxtiBWWsN+13kjM9nC6AWw0PL9AgCFB/5OM1KmSyr3vkSeTpjWIjKrb/i049orvj4mS7VqL5z4FR9ONOaJhwOK6f4q5TQqJ7f5fVSz9RbCGuBeu0jvHVUorTbo9yGpZswGg1IwUs4sqQKiU03+tKvf9v2s1slt+LNbw7jizIXFL5wCv7jmbNfb83UDnir1IBjPAlhBREuIKAjgSgB3e7ymiiJ9sW4Wxmgi7RrotAa7B8aTpuj84m/c3yBMZYjUMEVV9vmy/s65JQiGdJGU29+oWljrJWQPpUIDh0J+3aKS9URmEZqRTdgU8rYOI+CYBuhl2nIxTprrXt8xYy0MIYQK4GMA7gOwHcBtQohtRLSRiDYCABHNIaJeAJ8C8EUi6iWiadNVL98bXwiB0bjqmqe9al6b+fXAeAppLYPOaLCu37wzAfn8WkfeVosVs5oB6FXS//7OU9HdHML8jkiRn8qOS612+4lSsWYDXXaqnoBRaOCQTAndZYyilbGKeFp3rcjPg1dBb79l8uG9n3i9J2soFSLC//zVmfjMW1bikpPn4BNvXgGgehZGXbzjhBD3ArjXcdtNlq/7oLuqpiVWwVj++Xux9csXoynox59fGUBKy7j6rXtaQth7w2U47V/ux9BEEuokxkQy5aP4CD/+4FqcbBHsanGcIUqHRxK44syFJbsv5KjeehEMObP65qvPMHtKhQoIhhQTObtc1inJOoL2iO5y88zCsHzOTpoG3Z7fuHIW3rhyFgDgmT1D+O5Dr+CnT+7D6cd1uNZ+TQXegWqA9aSkZgQOjySwu38M7//xMwAKV4K2hAMYT6hIaeUPomcmx5tOmI1ZrdUPJp80Vxel16/oKXKlHRm3iJbZEK9aHG/UKKya12ZaSGct7cp7vXRB9Rkz52VGTyyp/y8/D17VYVSzJUy1kRXp92w9VJXCx/o4osxwnCclLSPwzJ7slLJC07iaQ36MJ1WEA0rBdgvM9GNOWxh/+vQbzHbepSLDW/USw/jOlafh0HAckaCCUxe048FPXYBlPdG811vngQBZ94ls7y8FI+WRS6pa7pxasHxWC779l2tw/OyWqriv6+MdN8NxCsZYQsWrlnYey3qa8/5sc9iPsYQKxUfskpqBLO7Ov7HmQ7p0ls/K/76pJW2RgM1KLrYuZ3xDNi6ULimvLQxZhf65S6pXh1NNLj+tsp1xrbBg1ABn8O4Tv3jeLFZ6y6rZZoGPGy0hP46MJRAN+dnCYADo41K/ecVqXDpNK/yt8Q3FRxhN6ONmJ5IqAgp5niWVSGl41xkL8BEuis2Bd6AakHKkuEmxWLe4Ez+4em3O8CMrzWE/RuMq9gxMFMw8YRoHIsI7Tl8wbTPmQpZ1z2+PQAjg9ud6EUtpaAr6zRoOr1xSspiSyYV3oBqQzDO/+Kylna63W4mG/Ng/FMOegQkcMypjGWY6Y3WtSnH47O1b9dnZQQVEhIBCnrmkEulMXRRE1iMsGDXArYjmslPm4qNvXF70Z60zNAbHWTCY6U/IMiRqXnu27iSW0syNOqD4PHFJCSEQz9N9gWHBqAlJNWN2GpWsXthW0pvSmmvv1u2WYaYb1hjGNecvRUvYjwUdEcRSqpn55feRJ4V78nDHLil3WDBqQDKtoSXsx86vrjdvKzWVstkyxGf9qjkVXxvD1BqrYLRFAnjXGQswEktjIqWZjQeDfm8sDDk6ttRRuY0GPys1IKlmEAooNt9todoLK1YL42tvP6Xia2OYWiML9wD9QNQeCWIsqWI0njYbD/p9Pk9GtMqUWrYw3GHBqAFJVUPI77O1KO+MltYBU44JPXVBG2dJMTMCq4XRHPKbrfxl8R8ABPxkszCEENjaO1y10aMSUzA46O0K70A1IKlmcnrrdJU4tnFhp16jMWD06GGY6Y714NMVDZqCMZpQzXGmAZ/PNlPjlqf24W3/+ThueXp/VdeWdUmxYLjBglEDkumMzQwHSp/zK6vA1y0pnoLLMNMBuRm3hv3wKz60W+J5XYblHVB82DuQnXr30kF97vf2w6NVXVuCXVIF4UrvGpBQNVsqIZCdKlaMcEDBI595A2bXoBkew9SCcEDBL685GyfM0TvBtlvaikiXkOxke9uzB3DFmQvRP6Y3KpxI5k6XqyTskioMWxg1YDSu5rQwt8YzirGoK8omMjOjOGtpF9oMV5R1HPH6k+2ZgA9sPwIA6DdcslUXjBRbGIVgC6PKZDICx2IpM2bx/feeXjddRhmmHpBT+q67/GSsWdhuu280ngYAc5zreJUFQzZAZAvDHd65qsxYQoWWEegwBOOSadowjmGqRVPQj703XOZ6n+w7NZbQhUJu6NVCtt9pLzCjppFhl1SVOTqu+17lYBOGYYrzllWzAeh9pzIZgXFjVka1LYxjMd2iKTTUrJFhwagyMrvj+NktHq+EYaYP37nyNMxqCSGlZTCeUs2hURNJFcOxFP7ie4/hNctMmUoxHEuZ2VtMLvysVJnbN/eipyWElSwYDFMy4YCCpT1RJFKa6Y5qDvkxkdTw4PZ+vHhwBN97ePeUfsfPntyLc69/yFYMuHcwhq7m0opqGxEWjCoyEk/jsd0DeM+64/jEwjBlEgkoSKgaxhK6m2hOWxgTKdWslXDOmSmXL921DYdGEqab65k9Q3h011HsG5wo8pONC+9iVUS+0Z2dahmGKU44oCBusTDmtoUhBPDF37wEwH1swGQ4MqrHGZ/dOwQAWFpgZHKjw4Lhwkgsjb6RxJQfR2Z0yJGTDMOUTo6F4SheTRm9pt73w6fx1XtenvTvOTKq13gcHtEnYd6+8ZxJP9ZMhwXDhTd84484+/qHpvw4sshIduBkGKZ0IkEFsaSG0XjWwrCSMiZZPrZ7AD98bM+kf8+IUevRN5LAiXNbba1KGDt1IRhEtJ6IdhLRbiK61uV+IqLvGvdvJaLTq7kemVo3VUwLg4uAGKZsuppDGIqlzNqIBR1NtvtTamZS3WtTagbLPn+v+f14QsWR0QQe3N6P7hLHDjQqnh99iUgBcCOAiwD0AniWiO4WQlhtzEsArDD+nQXg+8b/VSWTEfD5Sm/h4cS0MLiym2HKZnZrCEIArxrpsytm22MLz+0fxq4j5afWHouloFk64f5xZz9+/ozeBfdSLqwtSD1YGOsA7BZCvCaESAG4FcAGxzUbAPxM6DwFoJ2Iqv7KymDbZJkwio24zQDDlI+MWdzylL6ZL+6K5lxzx3O9ZT/usMOD8PuX+rDlwDDmt0dw1brjJrHSxqEeBGM+gAOW73uN28q9BgBARNcQ0SYi2nT06NGyF2M1caUpPFkODeuB81ktnNfNMOXi7NBsbVIoyVgsBavVUIh8n2t2RxWnHgTDzefjfOVLuUa/UYibhRBrhRBre3p6yl8MEX5w9RkApt6GYEffGGa3htAS5jYDDFMus1rtBy0iwn++5zTbbUOWzT+WKu3zKi2MOa1hBJTs1qKWKDiNTD0IRi+AhZbvFwA4NIlrKoacoz0VwRBC4P5tfThveXellsUwDUWXZYzxi1++GADw1lPn2a45NBw3vy61MeGwITJ3/N25ttv9U4hXNgr1IBjPAlhBREuIKAjgSgB3O665G8D7jWypswGMCCEOV2tBMkg9ld77sZSGpJrhliAMM0kUH+Ga85fiirULbFb6bR85Bz98/1oAwEGLYJT6eR020mjbIwGkNd2qmNMaxjfevbpSS5+xeJ6+I4RQiehjAO4DoAD4sRBiGxFtNO6/CcC9AC4FsBtADMBfVXNNzUah3VQsDBkwZ3cUw0yez196Ys5t65Z0YsRwK8k4IVC6hXEslkJQ8dnS3W+6+gys4MNdUTwXDAAQQtwLXRSst91k+VoA+Git1pO1MCbfe19Wp7aE6+IpZpgZRTioO0esge5SBWMklkZbU8A29XJRZ1OBn2Ak9eCSqjukVSDnCE+GUdPCYMFgmEoTVHyQ+72MOco54MUYS6jm53JRly4UcsAZUxjezVxoDvlxyvw2PPHqIPYPvoCzl3bhijMXFv9BC1kLg11SDFNpiAiRgIJYSsOahe14+fAoXj40WtLPJlUNYb/ujrrz786bcvp8I8EWRh7mtoUxGk/j188fxGfv2Fr2z8sYRitbGAxTFSLG+NZoSEF7JGAe0oqRVDMIBfStrzMaxDLuTlsyLBh5aA77MTSRe/LYNziBj//iecSL+Es56M0w1SUsBSPoR0skYLqBi5FUMwjyfJpJwc9aHlrDAfSPJXNu/7ufP4ffbjmErb3Drj+3u38Mj+46ykFvhqkyYcNKaAopaA37y7QwuF3PZODdLA/NeRoGylbI+fyef/OzzdgzMIH3nX0cfMSdahmmWsgebdGgH63hgK2IrxDJtIYQt+uZFGxh5KE5j2UgZ1vkG7C0Z0Af73hgKI7mkN+WuscwTOUIGG6lpqAfLWF/Sc1Cn9kzhERaQ9DPW99k4GctD52ONDvZlNBv9J5xi29Y6RtJcPyCYaqIbDwYDSklCcYTrw7gih88ib2DMYRYMCYFP2t56HGYrHJ+sCwOku0FrFibnx0eiXP8gmGqiAxy97SE0BoOIJ7WkNbyz/neNxgzvw752VU8GVgw8vD65d3YeMEyfOSCpQCARFoXChlYc5vK1z+aDZKPJlS0soXBMFVDentntYTNw9nOvjFby3MryXQ2s5EtjMnBz1oe/IoP115yAo4zWgYk0vrJRZq9wy5Bb2dWFVsYDFM9rnm9fphb3N2E1oh+OHvr9x7Dv9+30/V66SUAWDAmCz9rRZAVoeNJFaqWMd90zqldADAwzoLBMLXiynXHYfu/rMfctgi6m7Mu5F/nmcJn7TXFUzAnBwtGEWRx0IXffAR9o9nMqOF4roXhbK/MQW+GqS5y4z95fpt522CehJS4xSXF6e6TgwWjCNFQ9o21u18fOB8O+HBgKG4b5woACdUecGMLg2FqgzWrMd+oVmtSSoQL9yYFC0YRFnREzK9l730Zz/jZk/ts1yYc7ULYwmCY2pNvcp7dJcWHucnAglGE+e3ZPvmykvRtq/UxkQ9uP2K71mryAmxhMIwX+BW7YAghIIQwMx0BtjAmCwtGESJBBddtWAUgOw7yyjMX4qwlnTkDW+JpDQGFIA84LBgMU3sUR3eFJZ+7F//vti22z6tTVJjSYMEogXev1WdhHDymC0ZTyI/5HREcdvSuSaQ1hAMKFEMxOpp4KAvD1Bq3tj6/fv6gTTDy1WowhWHBKIFwQEE44DMtjGhQwdy2MI6MJW1vvERaQySgYH67HvdYNa/Vk/UyTCNy89VnAAD8vuy2ploqv+MpDa9f0Y3vXLkGb1k1p+brmwmwYJRIR1PQFIymkB9d0RC0jDC71wL6GzIcUPA/f7UON7zjFHQ1c0dMhqkVF6+ag6vWHWcr0ItZ4hbxtIaWsB8b1syHL09gnCkMO9lLpC0SwGGjQ200qKCrWXc3DU6kzHnAccPCWNIdxZLuqGdrZZhGJeT3IalmRSKWtH6tmnVVzORgC6NEXjs6YX7dFPSbed/WuRiJdAZhLghiGM8I+X1IWSyMcUsx7eBEigv2pggLRolY3KII+n1mQPuYpao0ntYQ5h41DOMZuoWRMYtqrd0XkmoGTVx/MSX42SuRxV1R7Ogbw7y2MIDsRL6xhIo7Nvci6PchkdZy5mgwDFM72uRBLpZGZzQbd5Rw/cXU8FQwiKgTwC8BLAawF8AVQohjLtf9GMBbAfQLIU6u5RolP/7gmXh89wDedcYCANkai/Gkin++e5t53SUnc/YFw3jFQqMzQ++xGDqjQew6Mma7n5sOTg2v/SfXAnhICLECwEPG9278BMD6Wi3KjXntEbx77UJz5KrM9R51DFIaTxYfE8kwTHVY0KF3ZjgwpFsWu46MYVFXk2lZcAxjangtGBsA/NT4+qcALne7SAjxKIChGq2pJEJ+BUG/D/uHYrbbnScahmFqx8JO3cL45gM7MRxL4b5tR3DS3FZzrk1bhPu7TQWvBWO2EOIwABj/z5rqAxLRNUS0iYg2HT16dMoLLERr2J8zNGntos6q/k6GYfIjG36+enQCa/7lAWgZgXeevsCsu2jn7gtTouoxDCJ6EICbY/8L1fh9QoibAdwMAGvXrq1q/X9LOGATjLZIAN949+pq/kqGYYpw/Oxm7Doybn5/5uJOs78bWxhTo+oWhhDiQiHEyS7/7gJwhIjmAoDxf3+111NJmkN+HLUIxty2MAfVGMZjrjB6vwHAZ96yEm1NAbNdTzjgtVNleuP1s3c3gA8YX38AwF0erqVsWsL+nLGsDMN4S6tlDs0pxiS+f3/XqbjhHadg5ewWr5Y1I/BaMG4AcBERvQLgIuN7ENE8IrpXXkREvwDwJICVRNRLRB/2ZLUOnO3LebA8w3iP9XPZ3hQw/g/iynXHmVmOzOTwtA5DCDEI4M0utx8CcKnl+6tqua5SaQ7Z/aEhLgpiGM+xTrrkQtrKwkfiKTCnzd6NlqtIGcZ7rPMw5rVFClzJlAsLxhRY3GXvSHvWUk6pZRivsbqkuI15ZeFeUlNgfkf29PLLa87GmYtZMBjGa3g0cvVgC2MKdFsGJK05rp1PMwxTB1izpJjKwoIxBayCEVT4qWSYekBmKx4/u9njlcw82HabAu2WqlFO12OY+oCI8NuPvc7sK8VUDhaMKcAuKIapT05Z0Ob1EmYk7EdhGIZhSoIFg2EYhikJdklNkZ9+aB2GY6niFzIMw0xzWDCmyAXH93i9BIZhmJrALimGYRimJFgwGIZhmJJgwWAYhmFKggWDYRiGKQkWDIZhGKYkWDAYhmGYkmDBYBiGYUqCBYNhGIYpCRJCeL2GqkFERwHsm+SPdwMYqOBypgP8N898Gu3vBfhvLpdFQgjXiuQZLRhTgYg2CSHWer2OWsJ/88yn0f5egP/mSsIuKYZhGKYkWDAYhmGYkmDByM/NXi/AA/hvnvk02t8L8N9cMTiGwTAMw5QEWxgMwzBMSbBgMAzDMCXBguGAiNYT0U4i2k1E13q9nlpARD8mon4iesnrtdQCIlpIRH8kou1EtI2I/t7rNVUbIgoT0TNEtMX4m7/i9ZpqBREpRPQ8Ed3j9VpqARHtJaIXiegFItpU0cfmGEYWIlIA7AJwEYBeAM8CuEoI8bKnC6syRHQ+gHEAPxNCnOz1eqoNEc0FMFcI8RwRtQDYDODymfw6ExEBiAohxokoAOAxAH8vhHjK46VVHSL6FIC1AFqFEG/1ej3Vhoj2AlgrhKh4sSJbGHbWAdgthHhNCJECcCuADR6vqeoIIR4FMOT1OmqFEOKwEOI54+sxANsBzPd2VdVF6Iwb3waMfzP+tEhECwBcBuCHXq9lJsCCYWc+gAOW73sxwzeSRoeIFgM4DcDTHi+l6hiumRcA9AN4QAgx4/9mAN8G8FkAGY/XUUsEgPuJaDMRXVPJB2bBsEMut834U1ijQkTNAO4A8A9CiFGv11NthBCaEGINgAUA1hHRjHY/EtFbAfQLITZ7vZYac54Q4nQAlwD4qOFyrggsGHZ6ASy0fL8AwCGP1sJUEcOPfweAnwshfu31emqJEGIYwJ8ArPd2JVXnPABvM3z6twJ4ExHd4u2Sqo8Q4pDxfz+AO6G72isCC4adZwGsIKIlRBQEcCWAuz1eE1NhjADwjwBsF0J80+v11AIi6iGiduPrCIALAezwdFFVRgjxOSHEAiHEYuif5YeFEO/zeFlVhYiiRiIHiCgK4GIAFct+ZMGwIIRQAXwMwH3QA6G3CSG2ebuq6kNEvwDwJICVRNRLRB/2ek1V5jwAV0M/cb5g/LvU60VVmbkA/khEW6EfjB4QQjREmmmDMRvAY0S0BcAzAH4nhPhDpR6c02oZhmGYkmALg2EYhikJFgyGYRimJFgwGIZhmJJgwWAYhmFKggWDYRiGKQkWDIZhGKYkWDAYpgSIqMtSs9FHRAeNr8eJ6L+q8Pt+QkR7iGhjgWteT0QvN0pbesZ7uA6DYcqEiL4MYFwI8Y0q/o6fALhHCHF7kesWG9fN6L5QTH3AFgbDTAEieoMczENEXyainxLR/cYQm3cQ0b8bw2z+YPSvAhGdQUSPGN1E7zPmcxT7Pe8mopeMAUiPVvvvYhg3WDAYprIsgz5/YQOAWwD8UQhxCoA4gMsM0fgegHcJIc4A8GMAXyvhcb8E4C1CiNUA3laVlTNMEfxeL4BhZhi/F0KkiehFAAoA2cfnRQCLAawEcDKAB/QeiFAAHC7hcR8H8BMiug1AQ3XXZeoHFgyGqSxJABBCZIgoLbJBwgz0zxsB2CaEOKecBxVCbCSis6BbLy8Q0RohxGAlF84wxWCXFMPUlp0AeojoHECfy0FEq4r9EBEtE0I8LYT4EoAB2Oe2MExNYAuDYWqIECJFRO8C8F0iaoP+Gfw2gGJt9L9ORCugWygPAdhS1YUyjAucVsswdQin1TL1CLukGKY+GQFwXbHCPQC/he6iYpiqwxYGwzAMUxJsYTAMwzAlwYLBMAzDlAQLBsMwDFMSLBgMwzBMSfx/JaPB7GXhAn4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Response of the first order system\n", "T, Y = ct.input_output_response(sys, T, V)\n", "plt.plot(T, Y)\n", "plt.xlabel('Time [s]')\n", "plt.ylabel('$Y$');" ] }, { "cell_type": "markdown", "id": "ead0232e", "metadata": {}, "source": [ "This is a first order system, and so we can compute the analytical correlation function and compare this to the sampled data:" ] }, { "cell_type": "code", "execution_count": 87, "id": "d31ce324", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "* mean(Y) [0] = 0.165\n", "* cov(Y) [0.05] = 0.0151\n" ] } ], "source": [ "# Compare static properties to what we expect analytically\n", "def r(tau):\n", " return c**2 * Q / (2 * a) * exp(-a * abs(tau))\n", " \n", "print(\"* mean(Y) [%0.3g] = %0.3g\" % (0, np.mean(Y)))\n", "print(\"* cov(Y) [%0.3g] = %0.3g\" % (r(0), np.cov(Y)))" ] }, { "cell_type": "markdown", "id": "28321bee", "metadata": {}, "source": [ "Finally, we look at the correlation function for the input and the output:" ] }, { "cell_type": "code", "execution_count": 88, "id": "1cf5a4b1", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEHCAYAAACjh0HiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdf0lEQVR4nO3deZRcZZ3G8ecHAWRYRiCRffEwKAITg7aRRRZZhrAIAzI4MIyRQeIIRCKMyCIDOggiB3BlCRCIioAQEVAwKFvkgCEdskAgkJBACMSkYwIhIQl092/+ePuee6u6qrqqU1W3u9/v55z31K33bm9t97l7mbsLABCv9fJuAAAgXwQBAESOIACAyBEEABA5ggAAIkcQAEDkBjVrRma2o6RfSNpGUqekse7+YzPbUtLdknaR9Jqkk9x9eaVpDR482HfZZZeGthcABpKpU6cudfchpfpZs64jMLNtJW3r7s+Z2WaSpkr6V0lfkbTM3X9gZhdI2sLdv11pWi0tLd7a2troJgPAgGFmU929pVS/pu0acvdF7v5cV/e7kl6StL2k4ySN7xpsvEI4AACaJJdjBGa2i6S9JU2WtLW7L5JCWEj6SJlxRplZq5m1trW1Na2tADDQNT0IzGxTSRMkjXH3FdWO5+5j3b3F3VuGDCm5mwsA0AtNDQIz20AhBO5w9992VS/uOn6QHEdY0sw2AUDsmhYEZmaSbpX0krtfm+n1gKSRXd0jJd3frDYBAJp4+qik/SX9p6TnzWx6V91Fkn4g6TdmdrqkBZL+rYltAoDoNS0I3P0pSVam96HNagcAoBBXFgP18swz0syZebcCqFkzdw0BA9t++4VH/uwJ/QxbBAAQOYIAACJHEABA5AgCAIgcQQAAkSMIACByBAEARI4gAIDIEQQAEDmCAAAiRxAAQOQIAgCIHEEAAJEjCAAgcgQBAESOIACAyBEEABA5ggAAIkcQAEDkCAIAiBxBAACRIwgAIHIEAQBEjiAAgMgRBAAQOYIAACJHEABA5AgCAIgcQQAAkSMIACByBAEARI4gAIDIEQQAEDmCAAAiRxAAQOQIAgCIHEEAAJFrWhCY2TgzW2JmL2TqLjOzN81selc5qlntAQAEzdwiuF3SiBL117n7sK7yUBPbAwBQE4PA3SdJWtas+QEAqtMXjhGcbWYzu3YdbVFuIDMbZWatZtba1tbWzPYBwICWdxDcIGlXScMkLZJ0TbkB3X2su7e4e8uQIUOa1DwAGPhyDQJ3X+zuHe7eKelmScPzbA8AxCjXIDCzbTNPj5f0QrlhAQCNMahZMzKzOyUdLGmwmS2UdKmkg81smCSX9JqkrzWrPQCAoGlB4O4nl6i+tVnzBwCUlvfBYgBAzggCAIgcQQAAkSMIACByBAEARI4gAIDIEQQAEDmCAAAiRxAAQOQIAgCIHEEAAJEjCAAgcgQBAESOIACAyBEEABA5ggAAIkcQAEDkCAIAiBxBAACRIwgAIHIEAQBEjiAAgMgRBAAQOYIAACJHEABA5AgCAIgcQQAAkSMIACByBAEARI4gAIDIEQQAEDmCAAAiRxAAQOQIAgCIHEEAAJEjCAAgcgQBAESOIACAyDUtCMxsnJktMbMXMnVbmtmfzGxO1+MWzWoPACBo5hbB7ZJGFNVdIOlRd99N0qNdzwEATdS0IHD3SZKWFVUfJ2l8V/d4Sf/arPYAAIK8jxFs7e6LJKnr8SPlBjSzUWbWamatbW1tTWsgAAx0NQeBmW1iZus3ojGVuPtYd29x95YhQ4Y0e/YAMGD1GARmtp6ZnWJmfzCzJZJmS1pkZrPM7Goz220d5r/YzLbtms+2kpasw7QAAL1QzRbB45J2lXShpG3cfUd3/4ikAyT9VdIPzOzUXs7/AUkju7pHSrq/l9MBAPTSoCqGOczdPyiudPdlkiZImmBmG/Q0ETO7U9LBkgab2UJJl0r6gaTfmNnpkhZI+rca2g4AqIMegyAJATN72t33qzRMD9M5uUyvQ3saFwDQOLUcLP5QcYWZHVDHtgAAclDNrqHEx83sPkmzJL0gabGkWxSOHwAA+qlagmC+pCsk7SXp05K2k/TdRjQKANA8tQTB++4+RdKURjUGANB8tRwjOKhhrQAA5KaaC8pMktz93Z6GAQD0P1VdUGZmo81sp2ylmW1oZoeY2XilF4UBAPqZao4RjJD0X5LuNLOPSnpb0sYKIfKIpOvcfXqjGggAaKxqLihbI+l6Sdd3XUE8WNJqd3+7wW0DADRBTXcf7bqC+GZJBzamOQCAZuvN/xEcJ2knM7vDzD5W7wYBAJqr5iBw9w53/5mksyWdYWbfr3+zAADNUssFZZIkMztG4eri3RXuP7Sm3o0CADRPVUFgZutJusDdr5C0haSHJV1TzV1HAQB9W1W7hty9U9JhXd2/dPcZhAAADAy1HCOYZmaXchUxAAwstRwj2FHSP0v6uplNljRT0kx3v6chLQMANEXVQeDuJ0mSmW0kaU+FUBguiSAAgH6s5rOG3H2tpOe6CgCgn+vNBWUAgAGEIACAyBEEABA5ggAAIkcQAEDkCAIAiBxBAACRIwgAIHIEAQBEjiAAgMgRBAAQOYIAACJHEABA5AgCAIgcQQAAkSMIACByBAEARI4gAIDIEQQAELma/7O4EczsNUnvSuqQ1O7uLfm2CADi0SeCoMvn3X1p3o0AgNiwawgAItdXgsAlPWJmU81sVKkBzGyUmbWaWWtbW1uTmwcAA1dfCYL93f1Tko6UdJaZHVg8gLuPdfcWd28ZMmRI81sIAANUnwgCd3+r63GJpPskDc+3RQAQj9yDwMw2MbPNkm5J/yLphXxbBQDx6AtnDW0t6T4zk0J7fu3uf8y3SQAQj9yDwN3nSfpk3u0AgFjlvmsIAJAvggAAIkcQAEDkCAIAiBxBAACRIwgAIHIEAQBEjiAAgMgRBAAQOYIAACJHEABA5AgCAIgcQQAAkSMIACByBAEARI4gAIDIEQQAEDmCAAAiRxAAQOQIAgCIHEEAAJEjCAAgcgQB4uEuTZ4cHvubjg5pypS8W4EBiiBAPH71K2mffaR77sm7JbW7/HJp+HDp2WfzbgkGIIIAjff++3m3IHjllfA4e3a+7eiNqVPD46JF+bZDCltUfeUzRV0QBGisefOkjTaSbr8975ZI668fHjs7821HbyRtTl5Dnr7znfCZrl6dd0tQJwQBGmvWrPB47735tkOS1uv6uvcUBJ2d0plnSjNnNr5NTz4pXXxxz8MlbV6vD/xkb745PK5YkW87UDd94FuFAS05MGuWbzuk6oPgzTelG26Qjjqq8W06+GDpiit6Hq4vBUHyWfbHg+4oqQ98qzCgJQuLZAE2Y4b097/n05Zqg6AvhFaxvIPgrbekl18ubANBMGAQBCjte9+TPvax2sa59VZp003DqY6JZAGWLFyHDZM++9m6NLFmyQIs275K+tKCLmlzXkGw/fbS7ruH7uSzzAbq8uWhfuLE6qfZ0SFttpk0blz92oleIQgGsunTpQ8+qG7Yq6+WzjknfX7ppdKcObXN7+yzpVWrpDVr0roTTgiP2bXsV1+tbbr1si5bBGvXSgsWrHsb5s2rPoiy8t4iyEren//+77QuOZ5y+eXVT2fNGmnlyvC9SZxyijR+fHXjr1olvfhi9fNDWX3gW4WGmDNH2ntv6fzzqxv+/POln/yk5+H+/Gfp2GNLry0ndaUWdPVYgD36qNTeHvbhH3ig1NYWDlj+4hflx/ngg+4L0WrPGnrrrfSA6Je/LO28c/XBWsrcudKuu4atLam2U0GLX0NHR3gvyg17882hrXPnSgcdFF7HqlXSU0/1vv2JpA2//31at8EG4bFUm1avDsdCpk8vPb3s9+XOO6WvfKW6dpx4orTnnuv2mUASQTBwtbWFx8mTax+33AJGCiHw4IPSe+9177d2bXgs9cOcNq32dmRNmiQddpj03e9KP/qR9Je/SLfdJo0aJY0cKT39dOljDxtuKB1/fOgutW/7yisr78446aTw+MAD4bH4/PnHHpPuu6/7eOPHS62thXVvvJG+FnfpU58qP9877pDGjk2fFx9r+eQn04Vv1pIl6fvywx+GUz0nTZIeeigsYA84QFq8uPx8q7FwYfe6QYPC4/Ll3ftNmRLOjvrGNwrrk3Dr7TUJjz0WHnuzhYUCBEF/9tRT0jvvVB6m0n7u66+XjjiisG7JksrT3Gij8Jjd/SNJTzyRdq9c2f3HPX9+4fPsWvmzz6btnDkzLCDHjAm7DJJz1ZOF15VXFq59JgulQw+VBg+WnnlG+vznC+efLMT/9rfweNddYUG7dq100UXSiBHlX+/EiSH4kvP329ul666TXnstbJkcemjY/bVqVTpOe3tY6H7mM2FNfObMcPwkCchBg0KYJe2Run9Op54qfe1r4XX8/Och+CRp6dLwmJyWK0lvvx2Ou0ybJm29tfTVr4b6trbCYzTJKbxJW5culUaPlk4/PWw1zJ0b6j/4QHruuXT6PW1BrVmTLoxffjns/soqd3yjeLrZXYadnWGLdsKE8vNNdlGVC4LFi7uHMUpz935XPv3pT/uA0dHhPnKk+5QphfWdnaFk/fGP7ied5P744+6TJ7tL7occUjit995znzPH/dprQ3/J/Yorus/zwQfT/u5p99Ch7n//e2G/jg73tWvdH3vMff31Q/0bb4R27L+/+/vvu192WTpOUl55pfD5T35S+Nzd/f770+eXX959GpL73Lnu997bvf6yy8L8s3V77BEen3++8HVlu0uVjo70/VmwoLDfeee5/+M/hu6pU8Pjbru5n3ZaOsyXvpR2/9//pd033JB2P/BAeBwxwv3zny8//zVr0vqNN+7e1uLXddttofvwwwuHO+cc9xNPLP1+traWfh/OOy/tnj07fO7Z/j/+ceHz6dO7T+Ohh9zfftt92DD3F15w//Of0/fs6afD6+vsdF+2rPzns2pVeFx/ffeHH+7+2zj//HTY++5zf/NN95UrC38z22yTfmaPPBI+o6ee6v5bKP6dTZwY3ocBRlKrl1mm9nphnGfpt0Fw1VXud99dWPfii+Fj2GabwvrddnP/8IcL60r9cDffPO3/T/9UfkGX9bOflV+wSN2D4Ctf6T69e+9132mn0D1lSvn5Viru7tdc0/NwgweHH3tx/cUXd68bOjQ8Tprk/vrr5V9jcbnoonQBPmxYYb/zznPfcsvu42QXvHvumXZnA2LMmO7jfeELIcCzdbvvHuZx8cUh7Ht637IL9Z/+NHQfeWT3YU84oXvdSy9V9/n85S/uy5fX/rlut5377beH7hNPdL/yyu7DXH65e1tb+c/nnXcKnz/+ePr9fe+98vM+5pjKv5dkXtlhjjyy9O9s5crC+htucL/5Zu+vCIK+IvmCdXa6X3114Zdzq61KD+vuvnCh+4oVpb/Um2zS/UdT6cs/e7b7GWd075d0DxpU+QeaLRttVPtCIlseeqj6YU8/vXvd8cdXP/7w4b1v55gxIYwqDWOWdu+zT+Vhjz/efYcdet+eUgt3Kd1qyZZkCylbsmvTlcq++7r/9rfr9hlvu23p+l13dV+8OH1e/D3LfgelsBB+9dUwXHZLolwp3pLJlhUr3P/2t+6/s+LfXlLuvLP8sP0IQdAoL77oPmqU+wcfuP/wh+5Llxb2z25yPvZY+kX62tdK/4jd3b/zncL6ataYq1kgPv982CIp1e/Xvy58fuaZPU+PQlmXUhyEye6tpIwbV3q8CRMKf0vlSjVhV/xbu/768BssNey3vpV2P/ts+vsu3q302mshtFatcj/1VPe33mrEkqdXCIJ6aWsLa+dnnlm4v/urXw2PRx+d7ut9/vm0f7K/s1L5h38I4+X9A6VQ8iif+Uz+bdhuu+p/g2+/nXYvWhTGa29Pd5eOHp32v+qqcDxr3ryw9Z6TPh8EkkZIelnSXEkX9DR83YJgwQL3Aw8M++3vuCPse3z6afezzgoHQJcvD4lfah9nubLzzt3X4nfeOf8vOYXSl8snPpF/GyT37bfv3XiXXuq+6aY9D7fhhuHxV78KwfHOO6GMGuU+bVrYpTV2rPsvfuF+6KEhcOqkUhBY6J8fM1tf0iuSDpe0UNIUSSe7e9lLBltaWry1N6eFLV8uHX10OF1w6tT0tEIA6Ku++MVwIeK0aeEivg037NVkzGyqu7eU6jdonRpYH8MlzXX3eZJkZndJOk5Sfa8df+cdacstQ/czz9R10gDQMNlrKTbeOFzMmVzPUyd94YKy7SW9kXm+sKuugJmNMrNWM2ttS66arUWd3zgAaLrOzvQq7jrqC0FQ6p6/3fZXuftYd29x95YhQ4bUPpcPfShcpTp/frjU/8ILpf/4j8JhLrmk9ukCQL3ceGPh84suCmXp0nD/sPb2xvxLXbmDB80qkvaVNDHz/EJJF1Yap24Hizs7w4Hd119Pn996a/nzwO+6K5wSlq379rfdP/e5ng8S7bNPGLan4Spd7HP22aXri6/YpVBiL+WuUbn77vLjrFzZ83Qvu8x9vfV6Hm7MmO4X+I0e7X7PPaWHP/5499/8Jl02zZzpfsst9VnOdVFfPmtI4TjFPEkflbShpBmS9qw0TlNOH33vvXDa5+zZ7u++G662TSQf3o03hvD43e/C8yuuCJewJ/3POy9cXj90aLis3t39iCNCv+QWEdmyZEkY5pJLwvNtt3U/7rjCL2GpL1Gli2colBhL9neaLX/4Q9p90UVp9x//GMaZM6f7OMnK2de/HoZZujRcxT9zZuFK4NSpYRgpPQMouf5gp53S5cfixe6rV4dTzN9/PyxrmqBPB0Fon45SOHPoVUkX9zR87heUHXZY+mVLZM8PPuCA0H/ixO7jdnSED7/42oJvfCMdJlmb+da33I86Kh3me99Lu/ffP9y7ZsKEMM4tt1T+YTzySAitPH+ckyZVN9wmm1Tuf8EF7uPHl+43alT189h8896/lrvvLh/M2fK737l///u9n0+y4lBN2X330vU//Wnpew5lS09XTkvhtifPPpvvd6i9PZzqXWmYyZPDb2LcOPcnnyx8bQ8/nHZfeGG61r5gQfr7++xnQ11ytbh7WKhn7weVuO66MMyYMaWXBcn1Bmee2eNipdH6fBDUWnIPgtWrQ6qX8/77YWHek4UL3Z95JnwMDz6Y1nd2hqt916wpvJ9N9qZsy5cXTiu5Z1GpkpXUJTcCq7Zkb2swYkQIuWrHvekm9//93zD/5LYU++3n/s1vlh4+Ode6XHnjjcLXki3Zm+1VKsl9nFasCOds1/JeLFsWxp02redhlywJC5Ds7RSqKS++GLZEy73OUuXkk8t//j2FcLlQ/Pd/d99rr9B9wAFhWmedVdvuyHvuKfz+ZG9sV02ZMKHwtRR/l4tL8dW+ixal/bLf2/PPD7+j++4rHP5HPwr9589Pb0VRTkdH2MoonmfWm2+Guw/kjCDo65LdRqUku5puvDEEh5ReAZlV6dhCVrYuO8zSpeE4xiWXhH7z54fL47ffPiyU3MPaWHt7Oq0ddyw/zzfeCG0+6KDC+SdB0NYWnl96afe1zAMPLD3Nb34z3GSt+LVkS7VBkL1ZX/Ea5he/mG7VFZdNNknHmzmz5/mU2qVYXEaPDvfeyda98krP4xWXL3+58uc/e3Z65WtxyV7Z29IStiCvuy6M98QTXhAE7mENea+9whZPpVs+HHFEOk57e7pW/cQT4TW/805YG+/sDFfoH3NMCMDi9kvuH/946e9ype978bArVoQtASm8H6V0dlb+TfZTBEEMkpvSnXtu5R/GwoXuM2aE7p5+PD1ZtizsU91333Q648aFO5OWk6ztF2/RZLcC3nqr8BYdSZk1q3CcY45J71lzxhlhV83q1YXjDBpUemGx6abpdIrvtXTuuaVvclc8Xqk2Fpck8Irf72y56aZwFWm2bs6cnscrFTrf/366a2PffcPCNeumm0K/7N1Up04N+8iT5yNHFo7z5JOh/nOfK/+5jh2bHiuTwhbArFnd7+BZreLv5uTJ3e/lVfz6k5spoiSCIFbTphWuWRZb1yDImjUrvRlXJRtsEOaX7PZIZINg1apQd+GFhbfrKCXZAjjnnLQuGf7VV8PCI7trISnZLYLkFtennRZu0f3ee93v0JrcPiAbBNkTA6SwZXHnnYV12V0L5RbgN98cAmPcOPe99w518+dXHu/RR0NgJseTPvWpdPjkPxomTer+fs2bF/rtskt4f6+9NtQ//ng67dNOKxwn2a20//6lP4OsJ54obHtvVfPdnDGj/Fo9uqkUBH3hymI0yrBh1Q13+OHrPq899qhuOPfwWHxRzNZbp3/lmFz8d8UV4fEjHwl/v1jL9CRpl13Cv2KdcEIoDzyQ/g1n9g/qv/CF8D/Co0dLH/5wqCv+N61Fi6TNNpO22Sat22cfafPN0/81PuWU0J6TTw7Pjz46tL1Um7PzX2+98O9qp50W/vHs3ntD28s591zpkEPSdkmF//aV3IKg1D937bBD+Oevq64q/Nz33DPt3mqrwnGStibvdSUHHdTzMNXYaSdpwYLKwwwdWp95gS2CqElhn3GTTl9z93RXzdq1hfWPPhrqhw6tbXorV4bdH9ldTfffHw5ol5KsZSa3/S4n+fOeL30p7DpxD2vuyTUnifb27muu5dZkX3ghXIviHtbAk1MPH3usclsqrR2vWRMOvD/zTFr3+uvhLJVaD1Am81ixorA+2fLZd9/aprcusv+xgboQu4ZQ0qpV4QynZpoyJZw1UnyWRfKXh7UGQa2S87qL//2tWGdndbu63MNBzOIgOPbYnsfr6KhuHskC8bbbqmtPbyXzyZ4Q4B4C5YwzCg/UN8OaNeGYD+qiUhDkfvfR3uj13UfRd82YEXZlDR0auhtlyZKwG2qLLaRly+ozzdWrwx+4b7FFeL50adhl1Mu7RHZTy66Zesyns7Nw1xUGhL5+91GgeZJ9//Vc0G28cSiJwYPrN+08EALR6Qs3nQOkj30sHCC95prGzidZyG26aWPn0x/98pfVn2CAAYUtAvQNG28c7gzbaFttFc5GOvHExs+rXlpbpb/+tfHzOfXUUBAdjhEAQAQqHSNg1xAARI4gAIDIEQQAEDmCAAAiRxAAQOQIAgCIHEEAAJEjCAAgcv3ygjIza5P0et7tqNFgSUvzbkST8ZrjwGvuH3Z29yGlevTLIOiPzKy13FV9AxWvOQ685v6PXUMAEDmCAAAiRxA0z9i8G5ADXnMceM39HMcIACBybBEAQOQIAgCIHEGQAzP7HzNzM+vnf27bMzO72sxmm9lMM7vPzD6cd5sawcxGmNnLZjbXzC7Iuz2NZmY7mtnjZvaSmc0ys3PyblOzmNn6ZjbNzH6fd1vqhSBoMjPbUdLhkhbk3ZYm+ZOkvdx9qKRXJF2Yc3vqzszWl/RzSUdK2kPSyWa2R76tarh2See5+yck7SPprAhec+IcSS/l3Yh6Igia7zpJ50uK4ii9uz/i7u1dT/8qaYc829MgwyXNdfd57v6+pLskHZdzmxrK3Re5+3Nd3e8qLBi3z7dVjWdmO0g6WtItebelngiCJjKzYyW96e4z8m5LTv5L0sN5N6IBtpf0Rub5QkWwUEyY2S6S9pY0OeemNMOPFFbkOnNuR10NyrsBA42Z/VnSNiV6XSzpIkn/0twWNV6l1+zu93cNc7HC7oQ7mtm2JrESdVFs8ZnZppImSBrj7ivybk8jmdkxkpa4+1QzOzjn5tQVQVBn7n5YqXoz+2dJH5U0w8yksIvkOTMb7u5/a2IT667ca06Y2UhJx0g61AfmhSsLJe2Yeb6DpLdyakvTmNkGCiFwh7v/Nu/2NMH+ko41s6MkfUjS5mb2K3c/Ned2rTMuKMuJmb0mqcXd+9sdDGtiZiMkXSvpIHdvy7s9jWBmgxQOhB8q6U1JUySd4u6zcm1YA1lYmxkvaZm7j8m5OU3XtUXwP+5+TM5NqQuOEaDRfiZpM0l/MrPpZnZj3g2qt66D4WdLmqhw0PQ3AzkEuuwv6T8lHdL1uU7vWlNGP8QWAQBEji0CAIgcQQAAkSMIACByBAEARI4gAIDIEQQAEDmCAAAixy0mgHVkZptLelLShgq3EXlF0hpJ+7n7gLo5GQYmLigD6sTMhivcaG9A34IaAw+7hoD62UvSQL+1BAYgggConz0kvZB3I4BaEQRA/WwnqV/fUhxxIgiA+pko6VYzOyjvhgC14GAxAESOLQIAiBxBAACRIwgAIHIEAQBEjiAAgMgRBAAQOYIAACL3//sZqnL0VHmJAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Correlation function for the input\n", "tau, r_V = ct.correlation(T, V)\n", "\n", "plt.plot(tau, r_V, 'r-')\n", "plt.xlabel(r'$\\tau$')\n", "plt.ylabel(r'$r_V(\\tau)$');" ] }, { "cell_type": "code", "execution_count": 89, "id": "62af90a4", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEGCAYAAAB2EqL0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABA50lEQVR4nO3deVhV1f7H8fdidEIBxREEUVRwwHnIedZMbbBrdbXUTJvtVrdsvL/qVtdsMLPBqSwbrK6aWg5pas4oToAMMimgKIggKDOs3x8gqRf1gGefc8Dv63l88uy9z17f8yR8zt5r7bWU1hohhBDiRuysXYAQQoiqQQJDCCGESSQwhBBCmEQCQwghhEkkMIQQQpjEwdoFGKVBgwbax8fH2mUIIUSVcuDAgbNaa4/y9lXbwPDx8SE4ONjaZQghRJWilDpxrX1yS0oIIYRJJDCEEEKYRAJDCCGESSQwhBBCmEQCQwghhEksFhhKqZFKqSilVIxSalY5+5VSal7p/hClVJfL9h1XSoUqpQ4rpWTokxBCWIFFhtUqpeyBT4FhQBKwXym1RmsdftlhowC/0j89gc9L/3vJIK31WUvUK4QQ4n9Z6gqjBxCjtY7TWucDy4FxVx0zDvhGl9gLuCqlmlioPiEsori4mCVLlpCdnW3tUoSoMEsFRjMg8bLXSaXbTD1GA78rpQ4opaZfqxGl1HSlVLBSKjg1NdUMZQthXjExMUybNo3PPvvM2qUIUWGWCgxVzrarV2663jF9tNZdKLlt9YRSqn95jWitF2qtu2mtu3l4lPtkuxBW5evri7OzMykpKdYuRYgKs1RgJAFel732BE6ZeozW+tJ/U4BVlNziEqLKOXPmDE2bNiU8PPzGBwthYywVGPsBP6VUC6WUE3AfsOaqY9YAD5aOluoFnNdaJyulaiulXACUUrWB4UCYheoWwqymT59OfHw8ERER1i5FiAqzyCgprXWhUupJYCNgD3yptT6qlHq0dP8XwDrgdiAGyAamlL69EbBKKXWp3u+11hssUbcQ5hYZGQlAfHw82dnZ1KpVy8oVCWE6i81Wq7VeR0koXL7ti8v+roEnynlfHBBoeIFCGCwnJ4f4+HjGjh3LiBEjKC4utnZJQlRItZ3eXAhbc+zYMbTWPPDAA0yYMMHa5QhRYTI1iBAWcqmjOyAggKioKOn4FlWOXGEIYSHdu3dn3rx5tG7dmk6dOuHv78/KlSutXZYQJpPAEMJCWrVqxVNPPQWUXGUcPXrUyhUJUTFyS0oIC9m4cSMnT54EwN/fn5iYGPLy8qxclRCmk8AQwgKys7MZNWoUixYtAkquMIqKioiOjrZyZUKYTgJDCAuIjIxEa027du2AksAA5AE+UaVIH4YQFnCpv6J9+/YAtG3bltWrV9O7d29rliVEhUhgCGEBYWFhODo60qpVKwBq1KjB2LFjrVyVEBUjt6SEsICjR4/Spk0bHB0dy7YdPnyYb7/91opVCVExEhhCWMC8efP48ssvr9i2fPlypk6dSkFBgZWqEqJiJDCEsABfX1+6d+9+xbYOHTpQUFBAVFSUlaoSomIkMIQwWHx8PPPnz+fqVSA7duwIQGhoqDXKEqLCJDCEMNi2bdt46qmnSE9Pv2J7mzZtcHBwICQkxEqVCVExEhhCGOzo0aM4OzvTsmXLK7Y7OTnh7+8vVxiiypBhtUIY7OjRo/j7+2Nvb/8/+1avXk2jRo2sUJUQFSdXGEIYLDQ0tOyBvau1aNFCVt0TVYYEhhAGSk9P59SpUwQGlr9o5KlTp3jhhRfktpSoEuSWlBAGcnNz4/z589dcjrW4uJg5c+bg4+NDhw4dLFydEBUjgSGEwVxcXK65r1mzZri6uspIKVElyC0pIQz08ccf8/77719zv1KKjh07yi0pUSVIYAhhoKVLl/LHH39c95gOHToQGhqK1tpCVQlRORIYQhikoKCA8PDwa3Z4X9KhQwccHR1JSUmxUGVCVI4EhhAGiYyMJD8//4rAKC7WfLv3BBfzCsu2Pfzww5w9e1aexxA2TwJDCIMcOXIE4IrA2B2bxqu/hPH2ur9W2nNwcEApZfH6hKgoCQwhDHL+/HkaNWpE69aty7YdT7sIwInS/17y0ksv8fzzz1u0PiEqSgJDCIM88cQTJCcn4+Dw1+j1qNNZAGTlFl5xbFxcHCtXrrRofUJUlASGEAa6+lZT1JmSwIg+c4Hi4r9GRXXu3Jn4+Pj/mdFWCFsigSGEAU6fPk2XLl3YtGlT2TatNcfOZGGnIKegiKT0nLJ9Xbp0AUqWbRXCVklgCGGAw4cPc+jQIZycnMq2pWblkZFdwJjApgAcK73agJIrDICDBw9atlAhKsBigaGUGqmUilJKxSilZpWzXyml5pXuD1FKdblqv71S6pBS6ldL1SxEZV26Uri0qh78dTvqjo5Nr3gN4OHhwZAhQ2TmWmHTLDKXlFLKHvgUGAYkAfuVUmu01uGXHTYK8Cv90xP4vPS/l8wEIoC6lqhZiJsRHByMr68vbm5uZdsikjMB6OrtRjPXmkRfFhgAmzdvtmiNQlSUpa4wegAxWus4rXU+sBwYd9Ux44BvdIm9gKtSqgmAUsoTGA0stlC9QtyUAwcO0K1btyu2HUk6j6dbTdxrO+HXqA5RZy78z/u01jJFiLBZlgqMZkDiZa+TSreZesxc4AWg/DmiSymlpiulgpVSwampqTdVsBCVVVBQQK9evRgxYsQV20OTztPRsx4AbRq5EJtygfzCv/5J7969Gw8PD4KCgixarxCmslRglPcY69Vfo8o9Ril1B5CitT5wo0a01gu11t201t08PDwqU6cQN83R0ZEffviBqVOnlm1Lv5hPwrlsOjRzBaBzczfyi4oJScooO8bT05O0tDTp+BY2y1KBkQR4XfbaEzhl4jF9gLFKqeOU3MoarJT61rhShbg5ubm5/7NtT1waAN18Svo0erRwByAo/lzZMV5eXtSvX59Dhw5ZoEohKs5SgbEf8FNKtVBKOQH3AWuuOmYN8GDpaKlewHmtdbLW+iWttafW2qf0fVu01hMtVLcQFTZx4kT69OlzxbatkSnUreFAZy9XANxrO9G6UZ0rAkMpRZcuXeQKQ9gsiwSG1roQeBLYSMlIp5+01keVUo8qpR4tPWwdEAfEAIuAxy1RmxDmduDAATw9Pctea63ZdiyVfq09cLD/60eul299go+fI6+wqGxb586dCQsLIz8/36I1C2EKiy3RqrVeR0koXL7ti8v+roEnbnCObcA2A8oTwizS0tI4fvw4jz/+1/edo6cySc3KY1Cbhlcc29/Pg2/2nCD4eDp9WjUAYNSoUQDk5ORc8dCfELZA1vQWwowOHCgZm9G1a9eybX8eKxmxN6D1lQMxbmtVHyd7O7ZFpZQFxsCBAxk4cKBlihWigmRqECHM6FJgXJobCkr6Lzo0q4eHi/MVx9ZycqBHC3e2RV05BDwnJ4e4uDjjixWigiQwhDCjPn368MYbb+Dq6grA+ewCDiakM7BN+cO8B7bxIDrlAicz/pqIcMKECYwZM8YS5QpRIRIYQphR//79ef3118te74hJpVjDwKv6Ly65FCR/XnaV0a1bNyIiIsjMzDS2WCEqSAJDCDPJysri4MGDFBQUlG3bGpmKay1HOpUOp71aS486NHOtyZbIlLJtPXr0QGtddntLCFshgSGEmWzbto2uXbuWTe1RXKz581gq/f08sLcrf81upRRD/BuyMyaV3IKS4bXdu3cHYN++fZYpXAgTSWAIYSZ79+7F3t6+rMP76KlMzl7Iu2b/xSVD/BuRW1DM7tizANSvX5+WLVtKYAibI8NqhTCTvXv3EhgYWLamxbaoFJSC/q2vHxi9fN2p7WTPpvAUBrdtBMDcuXNp2LD8fg8hrEUCQwgzKCoqYv/+/UyaNKls29aoFDo2q0eDOs7XeSc4O9jTz8+DLZFn0Lo9SinuuOMOo0sWosLklpQQZhAREUFWVhY9e5as+ZV+MZ/DiRkMuMboqKsN9m/Imcw8wksXWcrLy2P16tWEh4ff4J1CWI4EhhBm4Ovry6ZNmxg5ciQA26NLhtMOukH/xSX9/UqO2xVT0o9RXFzM+PHj+fZbmZhZ2A4JDCHMoFatWgwdOrSs32Fn9FlcaznS0dPVpPc3rleDVg3rsDOmZBr0mjVr0rFjR1lMSdgUCQwhzGD+/PkEBweXvd53/Bw9fNyvOZy2PH1bNWBffFrZ8NqePXuyf/9+ioqKbvBOISxDAkOIm5SZmcnTTz/NunUlkzGfPp/LibTsskWSTNW3VQNyC4o5mJAOlEwzkpWVRWhoqNlrFqIyJDCEuEn79+9Ha02vXr2AkqsLgJ4t6lfoPL1a1sfBTrEzuqQf49IiTHv37jVjtUJUngyrFeImXfqF3qNHDwD2xadRx9kB/yYuFTpPHWcHOnrWK1uFz9vbm8jISPz8/MxbsBCVJFcYQtykXbt24e/vXzZD7b74c3T1drtidT1TdW/hTkhSBrkFRSilaNOmDXZ28mMqbIP8SxTiJmitCQ0NpV+/fkDJ8xfHzlyocP/FJd293Sko0oQknQcgPDycadOmcfLkSbPVLERlSWAIcROUUsTHxzN79mwA9pf2X1Q2MLp6u11xnpycHJYsWcKOHTvMUK0QN0cCQ4ib5ODgcMXtKCcHOzp61qvUudxqO+HXsA7BpYERGBhI7dq12bVrl7nKFaLSJDCEuAmvvPIKb731VtnrfcfP0dnLFWcH+0qfs5uPO8En0ikq1jg4ONCzZ08JDGETJDCEqCStNV999RUREREAZOYWcPRUZqVvR13Sy9edrNxCjp4q6cfo06cPR44cISsr66ZrFuJmSGAIUUlxcXEkJyeXdXjviU2jqFjTt1WDmzrvbS1L3r+zdF6pvn374uvrS0JCws0VLMRNksAQopIudURfCowd0anUcrKnc3O3mzqvh4szbRu7lE1EOGzYMKKjo2nXrt3NFSzETZLAEKKSduzYgZubGwEBARQXa/6ISOG2lg1wcrj5H6v+rT3YF3+OjOx8lCqZj0prfdPnFeJmSGAIUUlubm7cfffd2NnZERR/juTzuYzt1NQs5x7TsSkFRZr1YacB+Prrr/H29iY3N9cs5xeiMiQwhKik999/n8WLFwPw/b4E6jg7MMy/kVnO3b5ZXVo1rMM3e05QXKxxd3cnMTGRPXv2mOX8QlSGBIYQlZCXl1f292V7T7D2yCkm9vKmplPlh9NeTinFk4NaEZGcyUebj9GvXz/s7OzYtm2bWc4vRGVIYAhRCTNnziQwMJCUzFz+/Ws4g9p48Nzw1mZtY2xgU8Z39eSTLTEkXIDOnTuzdetWs7YhREVIYAhRCVu2bKF58+b8fCCJvMJiXr0jAMdKTDZ4PXZ2in+NCaBuDQcW/BnHoEGD2Lt3L9nZ2WZtRwhTWSwwlFIjlVJRSqkYpdSscvYrpdS80v0hSqkupdtrKKX2KaWOKKWOKqXesFTNQpQnMTGR6Oho+g8YxLd7T3Bby/q09KhjSFsuNRz5WzcvNh49zYDho3jyySclMITVWCQwlFL2wKfAKCAAuF8pFXDVYaMAv9I/04HPS7fnAYO11oFAJ2CkUqqXJeoWojxbtmwBwM6zA8nnc3mkv6+h7f29lzdFWhOlm/Hhhx/SoMHNPRgoRGVZ6gqjBxCjtY7TWucDy4FxVx0zDvhGl9gLuCqlmpS+vlB6jGPpHxmQLqxmy5YtNGjQgPVJDrRt7MLA1h6GtteiQW2GBzTi693HOZt5kZCQEEPbE+JaLBUYzYDEy14nlW4z6RillL1S6jCQAmzSWgeV14hSarpSKlgpFZyammqu2oW4wvjx47n7kWeJOZvNowNalj1YZ6SnBvuRmVvIA489T9euXbl48aLhbQpxNUsFRnk/UVdfJVzzGK11kda6E+AJ9FBKtS+vEa31Qq11N611Nw8PY7/1iVtXvyEjCK7Zjc7NXRkbaJ4H9W6kfbN6DAtoRBSeFBYWyvoYwiosFRhJgNdlrz2BUxU9RmudAWwDRpq9QiFMEBoayvML15Kenc/bd3bAzs74q4tLXhsdgFOzAOwcHNm0aZPF2hXiEksFxn7ATynVQinlBNwHrLnqmDXAg6WjpXoB57XWyUopD6WUK4BSqiYwFIi0UN1CXOHZl15j6euPMrWPDwFN61q07eb1a/GP2zvg1CyAFWvWWbRtIcBCgaG1LgSeBDYCEcBPWuujSqlHlVKPlh62DogDYoBFwOOl25sAW5VSIZQEzyat9a+WqFuIyxUVFbF961ZcW3XmH8PaWKWG6f188QnszYmYSOJOyHTnwrIcLNWQ1nodJaFw+bYvLvu7Bp4o530hQGfDCxTiBn5ev4387EzG3jGa2s4W+9G5goO9He88/xiP1fZmW0Ievt5WKUPcouRJbyFMoLVm9uLloOx49ZG/WbWWO/t2oF///izYcYLcgiKr1iJuLRIYQpjgt9BkIoJ34RsQiK9XE6vWopRiXPMiIn9dxFc7Yq1ai7i1SGAIcQNaaz7bGkufJz5g7c/fW7scAFR6Apl7fuI/363n3MV8a5cjbhESGELcQNjJTMKTM5k80J8A/7bWLgeA4cOHA3A+OphfDp20cjXiViGBIcQN/LA/gYv7V5ASdPVIcOvx8PCgS5cuqKQjbDh62trliFuEBIYQ15FfWMyaQ0lkH1hDcJBtrXY3YsQIMk+EExSZQLrclhIWUOHAUErVLp19Vohq78CJdM4lxpB9Po2RI21rgoE77riDOi4u5J9NZEtkirXLEbeAGwaGUspOKfWAUuo3pVQKJU9ZJ5euTTFHKeVnfJlCWMf26FTy4vYDJd/obUmvXr1ITUmhuX8nNkecsXY54hZgyhXGVqAl8BLQWGvtpbVuCPQD9gL/UUpNNLBGIazmz6hU9IlgevbsSePGja1dzhXs7OxwcnJkSNuGbItKkWcyhOFMeVx1qNa64OqNWutzwApghVLK0eyVCWFlqVl5HE1Ko5FHfcaPv3r5FtsQFhbGsmfvouC2aeyJ68qgNg2tXZKoxm4YGJfCQim1W2t92/WOEaI62RGdirJ35KdVv9LBs561yymXj48PKclJ1I4LZlvkGAkMYaiKdHrXuHqDUqqfGWsRwqZsP5aKq0Mh7Sw8K21F1KlTh8GDB1MYv5+tUdLxLYxVkcBoo5RapZT6t1LqPqXUIGCpQXUJYVXFxZptYYkcff8+PvlknrXLua6xY8eSlZJIbPQxjp+VlfiEcSoSGPHAO0As0BWYBrxhRFFCWFvYqfOcOhpEYV4ugYGB1i7nuu644w4AsmOC2CZXGcJAFZmjOV9rvZ+SNSmEqNZ+DUkmJyYIVzc3+vbta+1yrsvLy4tZs2ax/lwDtkalMrlPC2uXJKqpilxhDDCsCiFsyIm0i3y1M5bC4we4Y/RoHByss/ZFRbz77rvcf9dodkSnEnk609rliGrKlAf3FIDWOutGxwhRHXy06Rj5CSHkXsjgnnvusXY5JhvUuAjnjOO88N8QStYjE8K8THpwTyn1lFKq+eUblVJOSqnBSqmvgYeMKU8Iy7qQV8iGo6e5a+htzJ071+ae7r6ehyfdh97zNSFJ5wlJOm/tckQ1ZEpgjASKgB+UUqeUUuFKqXggGrgf+EhrvdTAGoWwmA1hp8ktKGbSoI7MnDmTmjVrWrskk40fP55jIcHY52Sw4mCStcsR1dANA0Nrnau1/kxr3QfwAf4DdNZae2utH9FaHza4RiEs5pdDJ3HLTiJ06xpyc3OtXU6F3HvvvWit8Twfytojp8gvLLZ2SaKaqdBstVrrfKApcK8x5Qhxc3Lyi3j7t3AeXXaAxTviyM4vNPm9ZzJz2RV7FueYP3jqqaeqXD9AQEAAAQEBZIRvJz27gO3HUk1+b1GxZvm+BJ74/iAvrQwlNSvPwEpFVVXh6c211v8BCpRSHyml+iml6hhQlxCV8tm2GBbtiCc8OZN//xbBfQv3mjwp35rDpyguLORY0B+MHTu2St2OumT8+PFEhRyknn0BqyqwEt97GyOZtTKUwwkZ/BScyCurQg2sUlRVlVkP4zHgTiAM6AF8buaahKiUrNwCvt59nJHtGrP9hUF8/vcuhCSd5511ETd8b3Gx5sfgRJrmxJF+7hz33ls1L6KffvppkpOTubNHKzZFnCEz98bTvK0+fJIFf8bxQM/m7HxxEE8P9uP38DNEJMvwXHElkwKjdE2Ml0tfntZa36m1XqK1/kBrPcnA+oQwidaaD34/RmZuIY8NbAnAqA5NeLhvC77Zc4LPt8Ve9xbT5ogzxKRcoHZSEHXq1KlSo6MuV79+fdzc3LizczPyC4tvuN73D/sS+Od/Q+ju48b/jWmHUorJt/lQ28me2RsiKSySfhDxF5MCQ2tdDAwt/fsqQysSooIKi4p5bXUYS3cfZ2Kv5gR6uZbte3FkW0Z3aMLsDZFMWLj3mvf1F2yPw9OtBgXppxk/fnyVvB11yaFDh3jqgTG0rnmRBX/Gldv5HXU6ixnLgnlpZSg9W7izYFI3nBxKfh3Uq+XI8yPasC0qlUe/PUBOvqyzIUpU5JbUIaXUv+QhPWFLtNa8uCKUb/cmMGOAL2+Na3/FficHO+Y/0Jl37urA8bMXefDLfUz/JpjoM389h3ooIZ0DJ9KZ1teXP//cxuefV+27rA0aNGD37t00O3uAkxk5/HL4r6uM0+dzmbn8ECM/3s7umDSeHdaapVN64F7b6YpzTOnTgrfubM8fkSk8+GWQLM4kgIoFhhdwHyXLs65WSr2llKqaN3pFtfHZtlhWHEzimaF+vDTKn/K+zyileKBnc3a8OIh/DG3N3rg07v5sN2EnSx5uW/BnHHVrODC2Q8laEjVq/M9M/lWKl5cXAwcOZOeGlbRt7MLiHXForUnNyuOez3ez8ehpHh3Qku0vDOLpIX7Y25X/HXBSL2/m3deZ4BPpvLIqzMKfQtgikwNDa/03rbU/4E3JLLUxlHR6C2EVSenZzPsjmts7NGbmkBsvLe/sYM/MoX6sf6Y/dWs68sg3wXy79wQbjp7mbv86tPBqyvLlyy1QufEmTZpETEwMA90yOHbmAvO3xDBz+SHSLubx04zevDiyLW5XXVWUZ0xgU54c1IoVB5NkJlxRqWG1eVrrg1rrr7XW/zSiKCFM8fZvESgFr44OKPfK4lqaudbki4ldSc3K49Vfwuju40adk/vIzMykQ4cOBlZsOffccw81atQgIWgDYwKb8sGmY+yOTeOtce3p6OlaoXM9NdgP3wa1eXNtuDwMeIuz/Wk4hSjHryGnWB92mn+OaENT14p3UHfwrMfSKT0IPnGOKbe1YEj/p+ncuTPt2rUzoFrLq1u3Li+99BItWrTggQmd6ObtRt2aDtzV2bPC53JysOP1MQFM/mo/X+2KZ8aAlgZULKoCZamnWZVSI4GPAXtgcekDgJfvV6X7bweygcla64NKKS/gG6AxUAws1Fp/fKP2unXrpoODg838KYQtyMkvot97W2jmWpMVj92Gg32FL5SvcPToUdq3b8+HH37IP/7xDzNVWf1M+3o/e2LT2PL8QBrVrdr9POLalFIHtNbdytt3cz9pphdgD3wKjAICgPuVUgFXHTYK8Cv9M52/HggsBJ4r7T/pBTxRznvFLeTH/QmcvZDPq3cE3HRYACxatAhHR0cmTpxohupsS0ZGBmvWrDHLuV67I4CCIs3s9ZFmOZ+oeiwSGJR0jsdoreNK56NaDoy76phxwDe6xF7AVSnVRGudrLU+CGVrckQAzSxUt7AxZy/k8em2WHr4uNPdx90s55wxYwaLFi3Cw8PDLOezJR9++CF33XUXCQkJN30u7/q1eaR/C1YeOknw8XNmqE5UNZYKjGZA4mWvk/jfX/o3PEYp5QN0BoLKa0QpNV0pFayUCk5NNX3iNWHbktKz+WpXPCPnbqfH25tJv5jPrNvbmu38/v7+PPRQ9VzSZerUqWit+fLLL81yvscHtqJx3RpMWLiXwR9s4+3fwsuGJ4vqz1KBUd4Qlqs7T657TOkkhyuAZ7TW5U5yo7VeqLXuprXuVh2/Ld5qLuQV8u66CPrO3soba8NxdrTn8YGtWD+zH12au5mljTfffJM9e/aY5Vy2yMfHh+HDh7NkyRKKim7+4bvazg58/0hPHh/YEm/3Wny16zh3fLKT6d8Ec+5ivhkqFrbMUqOkkih58O8ST+CUqccopRwpCYvvtNYrDaxT2Ijk8znct3AvJ9KymdDNi2n9WuDXyMWsbURFRfGvf/0LZ2dnevfubdZz25Lp06dzzz33sGHDBkaPHn3T5/P1qMNzw9sAcO5iPj/sS2Du5mMMeG8rn03sQj8/+bJWXVnqCmM/4KeUaqGUcqLkifGre+LWAA+qEr2A81rr5NLRU0uACK31hxaqV1hJQVExi3fEMXLuDtIu5PPDI72YPb6j2cMCSjq7HRwcmDx5stnPbUvGjBlDo0aN2LVrl9nP7V7biScGteK3p/vRzK0mU77az7vrIriYZ/o6JKLqsOSw2tuBuZQMq/1Sa/22UupRAK31F6XBMJ+SJWGzgSla62ClVF9gBxBKybBagJe11uuu154Mq616DpxI55VVoUSezqKfXwNeHNmW9s3qGdJWTk4OXl5eDBgwgBUrVhjShi1JT0/Hzc08t/GuJTO3gDfWhLPiYBJN69Xg3Xs6MqC1XG1UNdcbVmuxwLA0CYyqZdWhJJ796QiN69bgjbHtGN6usaHtLVmyhGnTprF161YGDhxoaFu2JCcnx/CZeA+cOMesFaHEpF7gldv9mdbP19D2hHlZ/TkMIa6nqFjz/sZjdPR0ZdOzAwwPC4CioiKGDh3KgAEDDG/LVnzyySe0aNGC7OxsQ9vp6u3Omif7MrhNQ97bECXLvVYjEhjC6rYfS+VkRg6P9veljrNlxmFMnz6dTZs2VWgOqqquU6dOnDlzhmXLlhneVk0ne166vS35RcWsPmz6UrHCtklgCKtbc+QU9Wo6MsS/kUXa27t3r1mGmFY1ffv2pUuXLsybN++6qw+aS6uGLnT0rHfFehyiapPAEFaVW1DE70dPM6p947IV34wUGxvLbbfdxvvvv294W7ZGKcXMmTMJDw9n8+bNFmlzXKdmhJ3MJCYl68YHC5sngSGsaktkChfzixgT2NQi7X3yySfY29szadKtuRT9hAkTaNSoER9/fMP5O81iTGAT7BT8cujqx65EVSSBIaxq7ZFTNKjjTC/f+oa3dfbsWRYtWsT9999P06aWCShb4+zszDfffMNnn31mkfYautSgT6sG/HL4pEVugwljSWAIq7mQV8iWyBRGd2h8zWVCzWnevHlkZ2cza9Ysw9uyZcOHD6d58+YWa+/OTs1ISs/hwIl0i7UpjCGBIaxmc/gZ8gqLLXI7SmvNxo0bueuuuwgIkNnxo6OjGTlyJDExMYa3NaJ9Y2o42knndzUggSGsZu2RUzStV8NsEwlej1KK3bt3s2jRIsPbqgpcXFzYtm0b7733nuFt1XF2YFhAY34NSZYlXqs4CQxhFeezC9gencrojk2wM/h2VF5eHtnZ2djb21O/vvF9JVVB48aNefjhh1m6dCknTxr/zf+uzk3JyC5g+zFZdqAqk8AQVrHx6GkKirRFbkctXrwYb29vEhMTb3zwLeSf//wnxcXFzJ492/C2+vl54F7biVVyW6pKk8AQVrH6yEm869eig0GTC16SnZ3Nv//9bwICAvD09DS0rarGx8eHyZMns2DBArOsyHc9jvZ2jO7QhM3hZ8jKLTC0LWEcCQxhcWEnz7MrJo3xXTwNn5pj/vz5nD59mrfffvuWmgbEVK+//jqvv/467u7mWe72eu7s3Iy8wmJ+C0k2vC1hDJmtVlhUUbFm2Ed/cjGvkA0z++NW28mwts6fP4+vry89evRg/fr1hrUjTKO1Zsz8nSSey2Hr8wNxN/D/vag8ma1W2IyguDTiUi/yyugAQ8MCYOXKlZw7d45///vfhrZTHaxatYpXXnnF0DaUUswZH8j5nAJ+Dpb+pKpIAkNY1NqQU9R2smeYBSYanDJlCiEhIXTt2tXwtqq6oKAg3n33XUJCQgxtx79JXbr7uLF8f6I8+V0FSWAIi8kvLGZd6GmGBTSippO9oW2lp5c8VdyhQwdD26kuXnjhBVxdXXn22WcN/0U+vqsn8WcvEnryvKHtCPOTwBAWszMmlfM5BYztZOxQ2rCwMDw9PVmz5upl48W1uLu788Ybb/DHH3+wdu1aQ9sa0a4xjvaKX6Xzu8qRwBAWs/ZIMvVqOtK3lXHrPGutefbZZ3FycqJPnz6GtVMdPfroo7Rt25bnnnuO/Px8w9pxreVEPz8Pfj1yiuJiuS1VlVhmeTNxy8vJL1n3YkxgU0PXvVi3bh2bNm1i7ty58lR3BTk6OvLpp5+SnJyMg4OxvxrGBDZhS2QKhxLT6ept/JBeYR4SGMIitkaVrHsx1sAnu/Py8nj22Wdp3bo1jz/+uGHtVGeDBw+2SDtD/Rvh5GDH2iPJEhhViNySEhax+vBJGtRxpqeB6178+eefxMbGMm/ePBwdHQ1r51Ywf/58Hn30UcPO71LDkcFtGvJbaDJFcluqypDAEIY7mZHD5ogU7unazNB1L4YPH05sbCwjRowwrI1bRUpKCgsWLGDTpk2GtXFHYBNSs/IIik8zrA1hXhIYwnDf7DkOwIO9fQw5v9aaS0/1e3t7G9LGrebll1/Gz8+Pxx57jJycHEPaGNK2ES7ODizbc8KQ8wvzk8AQhsrIzueHoARGtGtEM9eahrSxbNkyunfvzubNmw05/62oRo0afPHFF8TGxhr2pHxNJ3um9G3B+rDThMkzGVWCBIYw1Md/RHMhr5Cnh/gZcv6TJ08yc+ZMevfubbEO21vF4MGDeeihh3j//fdJTjbmmYmH+7agbg0H5myMkie/qwAJDGGY30KS+WrXce7v0Zy2jeua/fxaa6ZNm0ZeXh5ff/01dnbyz9ncPvzwQ37//XeaNGliyPnr1XTk6SF+/HkslR/2yfxStk5+woQh1oUm84+fDtPV243X7jBmDe2FCxeyYcMG3nvvPfz8jLmCudW5u7szYMAAAE6cMKavYUqfFvTza8Crv4Sy6lCSIW0I85DAEGZ34MQ5nll+mA7N6rH4wW7UcDRm3ihnZ2fuvPNOeebCAlavXk3Lli3Ztm2b2c9tb6dYOKkbvXzr8+xPRyQ0bJishyHMKjUrj1Ef76C2sz2rn+iDay1jpzDXWsvCSBZw8eJFOnXqRH5+PocOHTJkwaWc/CKmLt1P8Ilz/DijN12au5m9DXFjNrEehlJqpFIqSikVo5SaVc5+pZSaV7o/RCnV5bJ9XyqlUpRSYZaqV1RcUbHmxRUhZOYWsGBSV8PC4uWXX2bp0qUAEhYWUrt2bb7//nuSk5N56KGHKC4uNnsbNZ3s+XxiFxrXq8Fj3x4gNSvP7G2Im2ORwFBK2QOfAqOAAOB+pdTVN7ZHAX6lf6YDn1+2bykw0vhKRWXlFRYxY1kwWyJTeHlUW0M6uQF++eUX3n33XQ4ePGjI+cW1de/enQ8++IBff/2VDz74wJA2XGs5sWBiN87nFDB16X7OXpDQsCWWmkuqBxCjtY4DUEotB8YB4ZcdMw74RpfcI9urlHJVSjXRWidrrbcrpXwsVKswkdaaz7bF8s2e46RfLCC/qJg3xrbjwd7GPDwXHx/PlClT6NatG3PmzDGkDXF9Tz75JEFBQdSoUcOwNgKa1uWzv3fh8e8OMvTDP6nlaM/QgEa8OjrA0IkrxY1ZKjCaAZePmUsCeppwTDPA5AHgSqnplFyd0Lx580oVKky3NiSZORuj6OfXAC/3Wgxo7cGIdo0NaSsrK4uxY8cC8OOPP+Ls7GxIO+L6lFIsW7as7FagUX1Ig9s24rtpPflq13Fy8ov4Zs8J3Gs78czQ1mZvS5jOUoFR3r+oq3vbTTnmurTWC4GFUNLpXZH3iorJLyxmzsZI/JvUZemUHobOEQWwZs0aIiIiWL9+Pb6+voa2Ja7vUkD89ttvzJkzh99++43atWubvZ2u3u5lM9k+9cMhPtsWy/iunni61TJ7W8I0lrq+SwK8LnvtCZyqxDHCRqw6lETiuRxeGNHG8LAA+Pvf/054eDjDhg0zvC1hGqUUO3bs4MEHHzSkE/xyL41qi52C2RuiDG1HXJ+lAmM/4KeUaqGUcgLuA65eP3MN8GDpaKlewHmttazhaIMSz2Xzn/WRBHq5MrCNcavnAXz77bfs2bMHgNat5XaELbn99tt5//33WblyJa+++qqhbTV1rcn0fr6sPXKKffHnDG1LXJtFbklprQuVUk8CGwF74Eut9VGl1KOl+78A1gG3AzFANjDl0vuVUj8AA4EGSqkk4F9a6yWWqL06KC7WRJzOZGf0WVKy8uju484Q/4Y42lf8+8Km8DP848fDKOCDewMNHda6du1aJk+ezJgxY1i1apVh7YjKe+aZZ4iIiODdd9+lSZMmPPXUU4a1NWNAS1YcPMnUpft5fnhrJvdpUeFzaK3ZGXOWoLhzpF3MZ0jbhvT1a2DYw6XVjTy4V839tD+R9zZGcvZCyRrNTg525BcW08+vAYsq+BR21Oks7vx0F36N6vDpA13wcjfuXvLOnTsZNmwY7du3Z8uWLbi4uBjWlrg5hYWFjB8/nqZNm/Lpp58a+iXiZEYOs1aEsCP6LJ/9vQu3d6jYHFfvrotgwfY47O0UNR3tuZBXSE1He6b29eG5YW2ws8DtVVt3vQf3JDCqsWV7T/DaL2F093Hjvu7N6evXALdaTvwUnMirv4QxpG1DPp/Y1aShipm5BYybv4sLeYX89lRfGtY1bljlkSNHGDhwIA0bNmTnzp14eBh720vcvPz8fBwdHVFKUVhYaOia4IVFxdz12W6S0rNZ+XgfWjS4cYe71pp5f8Tw0eZj3N/Di3+NaYedUuyNS+PnA0msPXKKB3t788bYdrf8w6A28aS3sJyo01n88+cjvPZLGIPbNuTbaT25p6snjerWwMnBjom9vHnrzvb8EZnCMz8eoqDo+h2WxcWa5386QuK5bD77exdDwwLg888/p06dOmzcuFHCoopwcnJCKUV0dDTt2rVjw4YNhrXlYG/Hx/d1QinFg18GkZKVe8P3vLcxio82H+Puzs14a1x7ajja4+RgR//WHsy7rxPT+/vyzZ4TzFh2gKjTWYbVXtVJYFQjWmveWRfBiLnbWXnoJDMG+LJwUlecHf73ttOkXt68OtqfdaGnmbn8EIXXCI2k9Gymfr2f38PP8NLt/nT3Mf8cQpfXDyXrSe/duxcfHx/D2hLGcHd3p3bt2owbN461a9ca1o6vRx2+nNydtAv5/H1REOGnMq957Bd/xvL5tlgm9mrOB38LxOGqvjulFC+NastLo9qyI/osI+Zu54nvDpJbUGRY/VWVBEY1smhHHAu3x3F/j+bsf2UoL43y/58fjstN6+dbFhqvrQ4jM7fgiv2J57IZN38XQXHneGNsO6b28TGs9n379tG3b1+Sk5NxcHCgWbNmhrUljFO/fn3++OMPAgMDufvuu1m5cqVhbXXycmXxg91Iz85n3Kc72RaVcsX+gqJiPvw9iv+sj2RsYFPeHNv+mreblFLMGNCSXbMGM3OIH+vCknlm+WGKiqvnLfvKkj6MKqy4WBOTeoHDiRmcysjh4z+iub1DEz65r3OFOu9mb4jk822xONorWjV0obuPG31bNeDDTcc4mZHDysduw6+RcZ3O69evZ/z48TRq1IgtW7bIlUU1cP78eUaNGsW+fftYu3Yto0aNMqyt9Iv5TFwSxPGzF/m/se2ISb3ArpizxKRcILegmHu6ePKfezpUaFTglzvjefPXcIb6N2KIf0NaNaxDJy/XSo0srGqk07saKSrW/Lg/kV8OnSTkZAa5BX/dSurh487XU3tQ06niQwQPJ2awPiyZiOQs9sWnkVtQjKO94svJ3ennZ1w/wtdff83DDz9Mx44dWbduHY0bGzO1iLC8rKwsXnvtNd566y3DR7mdPp/L2Pk7ScnKw8FO0c3HjXZN69HPrwED2zSs1DkvfZG6xMnBjvZN6zI0oBEP9fahtrOlJsqwLAmMaiIu9QKv/hLG7tg02jZ2oXfL+gQ0qUtXbzecHe1pUreGWYYF5uQXse/4OXzq18K7vvmnfLhk2bJlPPjggwwZMoSVK1dSt64xM9wK67t48SIfffQRL7zwAk5Oxkx7f/ZCHqEnz9PN2w2XGo5mO2d+YTFHEjM4mJDO/uPpHE7MwLt+LeZO6ETnarhmhwRGFVdYVMyHm46xcHscTg52/N/Ydtzb1bPKD/9LTU1l9uzZvPPOO4b9EhG24ccff+S+++6jT58+/Pzzz4atEW4JQXFpPPvTEU5n5jJziB+PD2x53b7CqkYCoworLtY89t0BNh49w/iunrw4si0eLlV3ptbY2Fjee+895s+fj6Ojeb4Fiqrhxx9/ZOrUqdStW5f//ve/9OnTx9olVVpmbgGv/xLGL4dP0d3HjaVTelSbW1TyHEYV9uWueDYePcMrt/vz/r2BVTosli9fTpcuXfjvf/9LdHS0tcsRFjZhwgSCgoKoU6cOAwcO5Oeff7Z2SZVWt4Yjc+/rzId/CyT4RDpvrg2/8ZuqAQkMG3YqI4c5G6MY6t+Qaf0qPm+Orbhw4QJTp07l/vvvp127dhw4cICAgKsXXBS3gvbt27N//37uvfdeunUr90tslXJ3F09m9G/Jj8GJHExIt3Y5hpPAsGEfbTqGBt4Yd+3x41XBpEmTWLp0Ka+++irbt2+XYbO3OFdXV77//ntatGiB1popU6ZU6cklnxrcioYuzry5Npzqeov/EgkMGxV9JosVB5N4sJc3zVxrWrucCsvKyuLChQsAvPLKK2zZsoW33nrL0DmGRNWTkZFBSEgId999NxMmTOD06dPWLqnCajs78M8RbTicmMG60KpXf0VIYNigS1N81HZy4IlBraxdToVorfn1119p3749L7zwAgDdunVj4MCB1i1M2CQ3Nzf27NnDW2+9xerVq2nbti0LFiwwfEEmc7u7iydtG7vw3sZI8gurVu0VIYFhgzYePc3WqFRmDvXDrXbVGW4aHh7OyJEjGTNmDLVr12bixInWLklUAU5OTrz66quEhITQpUsX3njjjbKr06rC3k7x4si2nEjL5ps9x61djmEkMGxMdn4hb64Np21jFybf5mPtcky2dOlSOnbsSFBQEB999BFHjhzhtttus3ZZogpp3bo1f/zxB3v27KFu3brk5+czc+ZM4uLirF2aSQa28WBQGw8++P0YieeyrV2OISQwbMziHfGcOp/LW3e2t/mHgc6dO0dCQgIA/fv3Z8aMGURHR/PMM8/IMxaiUpRSeHt7A3D48GEWL16Mv78/jz/+eNm/NVullOLfd3VAo/lo8zFrl2MI2/6NdIvJzi/kq13xDPVvaOg04jcrNTWV119/HR8fn7IlOX19ffn0009l/QphNj169CA6OpopU6awePFiWrVqxfTp07l48aK1S7umZq41ebC3D78cOklMStW6rWYKCQwb8tP+RNKzC3h0QEtrl1KuyMhIpk+fjpeXF2+99RbDhw/n7bfftnZZohpr2rQpX3zxBTExMTzyyCMcOXKEWrVKlgZOTEy0cnXlm9HfFycHOxZuj73xwVWMBIaNyC8sZuH2OLp5u9HNhq4uCgsLKSwsBOCHH35g2bJlPPTQQ4SHh/Pf//6X9u3bW7lCcSto3rw5n376Kbt370YpRWZmJu3ataNnz55899135OXlWbvEMvXrOPO3bl6sOnSS0+dvvBpgVSKBYSNWHkzi1PlcnhxsG8NoY2JiePnll2nevDlr1qwB4JlnniEhIYEFCxbg7+9v5QrFrcjevmTqfkdHR959910yMjKYOHEiTZo04cknnyQ+Pt7KFZZ4pJ8vxRqW7KwaHfamksCwARnZ+bz/exSBXq4MaG29PoDCwkLmzZtH37598fPzY/bs2XTt2rVsjQo3NzfpoxA2oWbNmjzxxBNERETw+++/M3LkSJYsWUJWVsl63NHR0cTGWu+WkJd7LcYFNuXrPSeIS60+fRkSGDbgnXURpGcX8M5dlp8CJD4+ns2bNwMl394++eQTsrKyeOedd0hISGDt2rUyPFbYLDs7O4YNG8b333/PmTNn6NixIwBvv/02rVq1IjAwkDfffJOwsDCLT9sxa1RbnB3seGllKMXVZKlXmd7cyvbEpnH/or08OqAls0a1Nby93Nxctm/fzvr161m/fj1RUVG4ublx5swZHB0dSU9Px82t+i0KI24tJ06cYNWqVaxYsYJdu3ahtWbQoEFs2bIFgPz8fIuswbJ8XwKzVoby7t0duL9Hc8PbMwdZD8NG5RUWMerjHRQUFfP7MwMqtbTqjeTk5BAUFETv3r1xdnZm1qxZzJ49G2dnZwYMGMDtt9/OqFGjaN26tdnbFsIWJCcns3r1apRSzJgxg+LiYpo1a0bLli0ZNGgQffr0oXfv3tSrV8/sbWuteWBREGEnz7P5uQE0qlvD7G2YmwSGjZq/JZr3fz/GV1O6M6iS6w5fLS0tjU2bNhEcHMyePXvYv38/BQUFbN++nX79+hEaGkpCQgIDBw6kdm3jll8VwlZlZ2fzzjvv8Pvvv3Pw4EGKiopQSjFnzhyee+45cnNziYuLo02bNmWd7Dfj+NmLDJ+7neEBjZj/QBczfAJjSWDYoJCkDMZ/sYeh/g357O9dK/z+ixcvEhkZSXh4OEeOHGHs2LH079+fXbt20bdvX5ydnenSpQv9+/enX79+9O/fHxcXFwM+iRBV14ULFwgKCmLXrl0MGzaM3r17s3PnTvr160fNmjUJDAykc+fOdO7cmTFjxpQNAKmoeX9E8+GmYyyY1JUR7Sp3DkuRwLAxBxPSeXjpfmo5ObDmyT7Ur1P+KnqFhYUkJSURHx9P/fr16dixI2fPnqV79+4cP3687DhnZ2fmzJnDU089RW5uLpGRkbRr106m5xCiElJSUtiwYQMHDx7k0KFDHD58mMzMTHbu3EmfPn349ddfmTt3Lm3bti3707JlS5o3b37NK5K8wiLGf76HqDNZLHqwm1VHQ96IBIaNyC0o4ufgRN5eF0FDF2fm39MGh7zzJCcnU7NmTW677Ta01owePZqoqCgSEhLKHpqbMmUKX375JVprJk+ejJ+fH/7+/gQEBNCqVSsJByEMUlxcTHx8PJ6enjg7O7Nq1Spmz55NREQEmZmZZcclJCTg5eXF999/z8aNG/Hx8cHHxwcvLy+aNGlC4+Yt+fuSfRw7k8WTg1rxcD9f6tW0vZ9bCQwLi4yMJDExkbS0NNLS0kg4dYaws0UkNulPZm4hOb/8i6yEcHJycsreM3z4cDZu3AjA+PHjcXR0pEWLFvj6+tKiRQvatGmDp6enVT6PEOJ/aa05c+YMERERxMXFMXnyZOzt7ZkzZw6ffPIJSUlJZUN57ezsyM/PJ6dQM+i+xwgJ2oFTXXd8m3sysLMf7Vv58MgjjwBw/PhxtNa4u7tTt25diw+1t4nAUEqNBD4G7IHFWuv/XLVfle6/HcgGJmutD5ry3vLcTGCEhYURGxtLVlYWmZmZZGVlobVm1qxZALzzzjts3bqVrKwssrKyyMjIwN3dndDQUAAGDx7M1q1brzinc9O2TJvzHQ/0aM7v387nwoULNG3alKZNm9KkSRO8vb1l6VIhqpH8/HySkpJITEzk3Llz3HXXXQDMnTuXH1euJvZ4IufOplCUk4V7wyacPX0SpRSjR49m3bp1QMmzUe7u7nTt2pX169cD8Prrr5OQkICLiwsuLi7UrVsXX19f/va3vwElX1jbtq38EH2rB4ZSyh44BgwDkoD9wP1a6/DLjrkdeIqSwOgJfKy17mnKe8tzM4Exbdo0lixZcsW2evXqkZGRAcCLL77I9u3bqVu3Li4uLtSrVw9vb29ef/11APbt28fJs5msicpkc1w2nfy8+Pj+rvh61KlUPUKI6intQh4zvz/An0cTGNa5Jf8c0YazMUeIjY0lLS2Nc+fOkZaWRr169Zg9ezYAEyZMICgoqOzLbGFhIf369WP79u0AzJgxgy+++KLSVya2EBi9gf/TWo8off0SgNb63cuOWQBs01r/UPo6ChgI+NzoveWpbGD8FpLMs4s3UJM8PD3q09LTg9aeHvRr24wOnvVM+p/wy6GTzFoZQnExTO7jw3PDW+PsYP5nLIQQVV9xsebLXfHM3RxNdn4hjw9sxdND/HByuPFEHAlpF9kclsSJlPOczXcgKT2b2NgYHhvbt9LLO18vMBwqdcaKawZcPhdxEiVXETc6ppmJ7wVAKTUdmA4ls1tWhqdbTe4f2oOM7HwSzmWzPfEiq6OOM+eP4wR61uPxQa0Y5t8IO7v/DY6iYs17GyNZ8GccPVu488HfAvF0q1WpOoQQtwY7O8W0fr7c29WLt34LZ/7WGA4lprNgUjfqOJf/K/pIYgYLt8exPiyZYg1ODnZ4utakmVtNxvTrQkCTuobUaqnAKO9r+dWXNtc6xpT3lmzUeiGwEEquMCpS4CWBXq4Eerlese3shTzWh51m0fY4Ziw7gF/DOjx4mw+j2jemQemQ2JMZOby8MpQ/j6UysVdz/jWmHY42vmKeEMJ21KvlyPv3BtLLtz4vrghh/Oe7mXtfJ9o2Lvnln5VbwI7osyzdfZx98edwqeHA9P4tub+HF15utcr9EmtulgqMJMDrsteewCkTj3Ey4b2GalDHmUm9vLm/uxe/hSazcHscr/0Sxmu/hNGkXg20htOZuTjZ2/HOXR14oGfVmDNGCGF7xnf1pKGLMzOXH2Lk3B00dHHGwU5xOjOXYl2yqt+ro/2Z0N0LlxqWHZZrqT4MB0o6rocAJynpuH5Aa330smNGA0/yV6f3PK11D1PeWx4jh9VqrQk7mcnu2LNEnc7Czk7h61GbMR2b4uUut6CEEDcvNSuPtUdOEZGciaYkKHr6utOzRX3sDbyasHofhta6UCn1JLCRkqGxX2qtjyqlHi3d/wWwjpKwiKFkWO2U673XEnVfi1KKDp716OBp/snKhBACwMPFmal9W1i7jCvIg3tCCCHKXO8KQ3plhRBCmEQCQwghhEkkMIQQQphEAkMIIYRJJDCEEEKYRAJDCCGESSQwhBBCmKTaPoehlEoFTli7jkpoAJy1dhEWJp+5+rvVPi9U3c/srbUudw3ZahsYVZVSKvhaD81UV/KZq79b7fNC9fzMcktKCCGESSQwhBBCmEQCw/YstHYBViCfufq71T4vVMPPLH0YQgghTCJXGEIIIUwigSGEEMIkEhg2TCn1vFJKK6UaWLsWIyml5iilIpVSIUqpVUopV2vXZBSl1EilVJRSKkYpNcva9RhNKeWllNqqlIpQSh1VSs20dk2WopSyV0odUkr9au1azEUCw0YppbyAYUCCtWuxgE1Ae611R0qW433JyvUYQillD3wKjAICgPuVUgHWrcpwhcBzWmt/oBfwxC3wmS+ZCURYuwhzksCwXR8BLwDVflSC1vp3rXVh6cu9gKc16zFQDyBGax2ntc4HlgPjrFyTobTWyVrrg6V/z6LkF2gz61ZlPKWUJzAaWGztWsxJAsMGKaXGAie11kesXYsVTAXWW7sIgzQDEi97ncQt8MvzEqWUD9AZCLJyKZYwl5IvfMVWrsOsHKxdwK1KKbUZaFzOrleAl4Hhlq3IWNf7vFrr1aXHvELJLYzvLFmbBalytlX7K0gApVQdYAXwjNY609r1GEkpdQeQorU+oJQaaOVyzEoCw0q01kPL266U6gC0AI4opaDk9sxBpVQPrfVpC5ZoVtf6vJcopR4C7gCG6Or7cFAS4HXZa0/glJVqsRillCMlYfGd1nqlteuxgD7AWKXU7UANoK5S6lut9UQr13XT5ME9G6eUOg5001pXxVkvTaKUGgl8CAzQWqdaux6jKKUcKOnUHwKcBPYDD2itj1q1MAOpkm89XwPntNbPWLkciyu9wnhea32HlUsxC+nDELZgPuACbFJKHVZKfWHtgoxQ2rH/JLCRks7fn6pzWJTqA0wCBpf+vz1c+s1bVEFyhSGEEMIkcoUhhBDCJBIYQgghTCKBIYQQwiQSGEIIIUwigSGEEMIkEhhCCCFMIoEhhBDCJDI1iBAWopSqC/wJOFEy/csxIBe4TWtdrSapE9WTPLgnhIUppXpQMulitZ7aXFQ/cktKCMtrD1T3KUFENSSBIYTlBQBh1i5CiIqSwBDC8poCVXaqenHrksAQwvI2AkuUUgOsXYgQFSGd3kIIIUwiVxhCCCFMIoEhhBDCJBIYQgghTCKBIYQQwiQSGEIIIUwigSGEEMIkEhhCCCFM8v+pN2UEIugU3AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Correlation function for the output\n", "tau, r_Y = ct.correlation(T, Y)\n", "plt.plot(tau, r_Y)\n", "plt.xlabel(r'$\\tau$')\n", "plt.ylabel(r'$r_Y(\\tau)$')\n", "\n", "# Compare to the analytical answer\n", "plt.plot(tau, [r(t)[0, 0] for t in tau], 'k--');" ] }, { "cell_type": "markdown", "id": "2a2785e9", "metadata": {}, "source": [ "The analytical curve may or may not line up that well with the correlation function based on the sample. Try running the code again with a different random seed to see how things change based on the specific random sequence chosen at the start.\n", "\n", "Note: the _right_ way to compute the correlation function would be to run a lot of different samples of white noise filtered through the system dynamics and compute $R(t_1, t_2)$ across those samples. The `correlation` function computes the covariance between $Y(t + \\tau)$ and $Y(t)$ by varying $t$ over the time range." ] }, { "cell_type": "code", "execution_count": null, "id": "bd5dfc75", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.13.1" } }, "nbformat": 4, "nbformat_minor": 5 }