Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Go to file
 
 
Cannot retrieve contributors at this time
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"from numpy import sin,cos,pi\n",
"from scipy.linalg import *\n",
"from scipy.optimize import fsolve,root\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import pretty_plots"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/ryan/bin/pretty_plots.py:6: MatplotlibDeprecationWarning: \n",
"The text.latex.unicode rcparam was deprecated in Matplotlib 3.0 and will be removed in 3.2.\n",
" plt.rcParams['text.latex.unicode'] = True\n"
]
}
],
"source": [
"pretty_plots.setdefaults()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from __future__ import print_function\n",
"from ipywidgets import interact, interactive, fixed, interact_manual\n",
"import ipywidgets as widgets"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Intro to functions and m-files"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Building our own functions."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Functions are saved in memory (or better yet) in a folder in your path or current directory\n",
"\n",
"Example of storing a function in memory\n",
"\n",
"$f(x,y) = (xy^{3}-x^{3}y)$"
]
},
{
"cell_type": "code",
"execution_count": 75,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 75,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f= lambda x,y: (x*y**3-x**3*y)\n",
"f(0.1,0.1)"
]
},
{
"cell_type": "code",
"execution_count": 84,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.contour.QuadContourSet at 0x7f44495caf60>"
]
},
"execution_count": 84,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEhCAYAAAAwMdReAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAciElEQVR4nO3dz44c13mG8feY8pAwIJGWQkNGKFEmF6GArMZaBdmZvAPK3mZjapVFNnK8y85Q7oD0HZjcZkUayAVIXAWQNhyTEoMIkC1TImCIAwkni+4Sa5rV3fXnfOd8p+r5AYQ542Z3zWi6nvmqTneFGKMAACjtB6U3AAAAiSABAJwgSAAAFwgSAMAFggQAcIEgAQBceKn0BowRQjgn6ZKkJzHGo9LbAwCYrtYJ6bakj9b/CwCYgeqCtJ6OrpbeDgBAWtUFSdJvS28AACC9qoIUQrgq6f3S2wEASK+KIIUQDkMINyXdLb0tAAAbblfZhRBua7WS7rD0tgAA7LkNkqTrpTcAAJCP5yAdSXq14/Pncm8IAMCe2yDFGC93fT6EcFcs+waA2aliUQMAYP4IEgDABYIEAHDB7TmkXEIINyTdkKRw+uDnP/zpeYXjsPffnTru/xinvoljN0+SFL4Z8GALE88cmN33d2f2/xz0up8EmxgPpv0MSdLBwbfTN0TSKy99k+R+Nv341N9M7neu/vrdj0b/26+/PdP7tsfH+zNx/PB//xxjPD96g9YWH6QY4y1JtyTp9M8uxJ/+x79Kkk5/2m8v8vKjYTuKsw+eDdvADgcffzb5Pmp3/PYbpvf/1eXTye7r6cXpYXv25vRfSt668MXk+2hce/2TZPfVdv2V+yb3Oyd3vp720sy7n18ZdPuHj/d35tG//PujsdvTtvggbdPsAPaF6enFMChKzY5uSpg2d8ZLCZR1hKS0IZLSxCiVh4/PJ4vS3c+vmESpvbMlTitTA9SwCFHfX9z7Ikh7PHvzuFeUpGHTUoowNeYcqBwRaqSOkUc1RKmx1DilCpA0PEJSvxBJ6WMkEaRehkxL0rgwSWniJHXvxGuJVM4AtVnEyNN0ZKXZ4VmGSZp3nFIGqFFbiBoEaYA+05I0LkxS2qlp064dfYlYlQrPJqupyHOMUk5JDetpqW1zB15ToCzi0xgTIal/iCTbGEkE6YSDg2/11oUvdv4H6jstSdPDJNnEadPQOGwLmJfI9LWEQ3Q55ZqWNm3byZcKlWV0uswhRA2C1GFflKT+05I0PkxS/jj1UVt4NlmHyPN01LCYkho5p6VdhoZhW8ByB6aPsRGShoVIyhcjiSBt1TdKUv//YFPCJPmMU01yTEQ1xCiHUtPSFB7D0zYlQpJNiMbuy7YhSDv0iZI0bFqSpodJIk5DcGium+WU1KgxTJ5MjZBkNxGljpFEkPZqnrCppyUpTZgk4rRN7hDVOB3liJJEmIZIESFpeIiksjGSCFJvVtOSlC5M0os74SUGqsREVGOMGrmiJBGmbVJFSLINkWQXI4kgDTIkStLwk4Epw9RYSqBKHparOUalLD1MKQPUGBMiqfxU1EaQWvq8aWTfKEm+wtSYU6A4N5ROzimprb1jnnucLCIk2YdIyhMjiSC94Nrrn+z9wel7Xqkx5jCeZBumRtdO3WukPAZoTtNRqSg15hQnq/i0jQ2R5DNGEkHq1CdKUp5pScoTprZtO/5cofIYni5zilGjdJQam88/z4HKEZ+2XCGS+u1zUu4XqgtSjPFajsexiJI0flqSTu4Ac/7W0qglFDnMMUYNL1Fq2/ZczBmq3OHZNCVEks1UlPqX1OqClJNllKRpr4DOPTXhuTnHqOExSl1KRyKHnCGSysVI4hLme/X9DWzMkzfFRdeeXgyL2EF6sLTv9dQdIcZ7+Pj893+mqClGEkHqZUiUhobp2ZvHScO0pB1mTkv9vhKlvFJESFqFqLYYSQSptyHHqktNSw3ClNbSv5dEyV7JEEk+YiQRpEFyRMkiTEvfoU7B924l1Q4Tz6U6LNcYe07aS4wkgjSYdZSktNNSgzANw/erG1GaLnXcp0xFnmIkscruhB+f+luv2/VdfScNX4HXSLESr0vppePeEaH9mp/nGlbheWEVcuupKDcmpA19rzKZY1KSbKalBof0TuL7MAzT0m6pD8m1jZ2KpOExyvnOLUxImYydlCS7aaltyZMTIRqPaemkHJGesh/wHCOJCamTxZQkTX/Spl70sM0SJqclfI05LXnRg+Uk1DZlKpLsYnTw8WdjNqcTE9IW11+53+uSxkPOJ0nD35i1y5S3Hxpqc4dd+/REgGy1f67nPDXlju/U53sNMZIIUhJDoyRNO4Qn5TmM16XGQBGhMuZ0OK/U9Jfi+e39MF0bQdqh75QklYmSlHda6tK1sy8dKQLkS41Tk4fDjyViNETq6UgiSMWlipKUf1raZlcQUj5BCE99Nn/WPQTKQ3zaUj2PxzzXSl8LjSDtYT0lSWmiJJWflvogImjb9nOfOlTeorNNLTGymI4kgpSchyhJfqYlYIxaApJKyudrrTGSWPbdS99l4FOl/K0wx/JwANNMXcq9qfT526kIkgEvl1vO9bolAMOlPooxNkZepiOJIPU2dEoaGyWLk7xECfAj9VQk5YlRDgTJIasoESagLItzu7UfpmsjSAPkmpIku+WwhAnIz2IqkqbFaOh0ZH24TiJIrlm+RoMoAXlYrXid02TUIEjGpi5wsI4SYQJsWE1F0vQYeZyOJII0WK4l4G3Wr2YnSkBavA5wHIKUgZdl4LswLQHTWU5FjdzTUU4EqRK53vOLKAHj5JiK5njeqI0gZZJiSsoZJcIE9JNjKpLSxGjMdJTr/JFEkKqT892RiRKwG+eK0iJII5RY2FAK0xLwolxTUaPUdJQbQcoo1eKGEteQIUpA/hBJ8z9v1EaQ0BvTEpaMw3P2CFKlSl5pkyhhSUpMRY1U01ENh+skglS10lEiTJi7klPRkg7VNQhSZjW8SHYIwoQ5KjkVLRlBqlzJKamNKGEuPIRoidORJL1UegMwH02UPDyhgaH4uS2PCQnJMS2hNp5ilHo6qmVBg0SQZsHLYbs2zi2hBpwr8oUgwRRRgleEyB+CNMKdrw9Lb8ILPE5JDaYleOJ5KlrqYoYGQUI2RAmleQ2RZ8dvv5HtsQgSsmJaQgmepyJrX10+XXoTeiNIKIIoIYeaQrT0w3USQRrM4/mjWjEtwVItIapBrsN2BGlGPC9s2IUoIaWapiKcRJAGSDEd3f38SoItmR+mJaRAiOzkmJIIElwhTBiDqWi3WhY2EKSeOHeUF1FCX4QoH+spiTdXzYjDdcPwZq3YhZ+L+WFC6oHpqCymJWwiRsOlOmxnOSURpD2IkQ+cW4LEuSIvrKJEkDLJcbju4ePz5o9RGlFaLkI0nffFDQRph1TTEeeO0mJaWhamorQ8H7ojSFtwqM4/ojRvSwvR04uh9CYMljpKBKlDyhgxHdliWpqnJYWohJSH7lJGiSBtqDVGSzh/tAtRmoelTUWbck5JHs8nEaSWv373o9KbgAmYluq25BCV4i1KBMkI01E5RKkuS5+KNuU+l+QpSgTJAOeNymNaqgMh8sFLlAhSYrljxHS0G2HyialotxIr7jxEiSAlxGTkF1HygxD1s8QoEaRESsSI6WgYpqWymIrqUDJKBGmiu59fIUaVIUr5EaJxSr1YtlSUCNIEpQ7REaPpmJbyYCqarmSUcoeJII1EjOaBKNkhROmUfFuhnFHiAn0DlVy4QIxscCHAtPg+2nh6MejlR7HIYzdROvvgmenjEKSeSq+gI0b2nr15zM50Ar539ppJaa5h4pBdD8RoOTi3NA4xyqv0O4NbnV9iQtqBEC0X01I/fI/KKXkIr/F9lP47zf0RpA6lQyQRIw84t7Qb35fySh/CS40gtXz97ZniMSJE/hCmk/g++DOXMBEkJwiRfxzGI0be1R4mguQAMarHUqelpX29tas1TASpIEJUryVNS0v5OueotjARpMyI0HzMfVqa69e1RLWEiSBlQojma47T0ty+Hqy0X7/kMU4EyRARWo65TEu1bz/68xgngpQYEVq2mqelWrcb03mJE0FKgAihrbZpqZbtRB6bb0uUM1DVBSmEcCjpPUlXJV2S9ETSl5LuSboZY7xvvQ0ECH14n5Y8bxv8yBmoaoIUQjgn6bZWIWo7t/5zQ9KNEMIdSb+OMT5J9dgECGN5nZa8bQ/qYRmoKoK0jtFHWk1E+1xf3+7nYx+PACE1L9OSh23AvKR85/FaLj/xR52M0ZGkdyVdXv95b/25xmEI4fbQBzk+fokYwUzpS1sQI3jnfkIKIVyXdNj61L0Y47WNm90KIfxBq3A1t70eQjjMcU4JGCL3YTxChFrUMCF90Pr7E60moxeszxlt/n8fdN0W8CDHtESMUBPXQQohXNLJQ3W/27VYIcZ4JOlW61NX1+efAJesDuOd/vSAGKE6roOk1bmhtludtzrp5sbHv0y0LYCZlFEiRKiV93NI7XNHR32WcscY74dwYtXH6NV2QE5Tzy0RItTO+4T0TuvvQxYntG/7ztZbAQ51TUv7YkOMMAfeg9Q+/3O09VYvat+2z2uXAFe6zi11RYdzRZgTt4fs1gsa2h4M+OftILGoAbNBfDBnboOkF0Py5YB/+5f2ByGEcynfSgiwRniwRJ6DtMkkKCGEG1q9D55OvcYwhfKGxKjkOz8AqXkO0qsT/u1mvF7t+JwkKcZ4S+vl5Gf+/o3Ib6aoASHCHHkO0hSMOqhaOzj8koSl8LzKbvOc0ZTIDDn/BBS1Of2MWQYO1MhzkLoOu/X1WvsDFjSgBrveRogoYQncHrKLMR5tvOPC5QH/nNceoRpTzgc1UeKcEubA84QknZyShkSmfVsuPwG3UoWEaQlz4D1IH7b+frj1Vi9q3/bDrbcCCrF4l2/etQG18x6ku62/X+pzKYmOd3i423lDoBDrw2tECbXyHqR7Gx/3uZTE5iUrNu8DKCLnJcyZllAj10FaX368/b50v9l1+/UEdaP1qXussIMHpRYdECXUxO0qu5bfSLq9/vulEMIHMcZtYfq9Tr5eaWfAAGseVr+xEg+pvXXhixMfP0p0v+6DFGO8E0K4r+cLFd5fT0IfrC9ZrhDCoaQPJF1t/dM76wmrt1PH0suP4vcfP70Ydtwa2M1bAAgTxtoMkBX3QVr7haQ/6fn0c0PSjY3XKbXdjzG+O/VB23GSCBT68b7DP/3pgfttRFm5ArSpiiDFGJ+EEH6m1aG7q3tufidFjLoQKOxTy46eaQmbSkWorYogSd+//c+1EMJVrVbSHer5C2CPtFpNd3PoYbopOLyHRq07dqalZfMQobZqgtSIMd6Tw6XcxGmZ5rAzZ1paFm8RaqsuSDUgTsswtx0409J8eY5QG0EyRpzmZ847baaleaklRA2ClFETJ8JUr6XsqJmW6lZbiBoEqQDCVJ8l7pyZlupTa4gaBKkgwlSHpe+QmZb8qz1EDYLkAGHyiZ3wc0xLPs0lRA3Xb666NJsvvEU57Hi78WatfniI0bXXP9G11z9Jdn9MSC2nvok6++CZvrp8utg2MC2VRYj2Y1oqy0uILBCkDmcfPJOk4mEiSnmxgx2GMOVXOkZWIWoQpB1Kh4lpKQ92qNOw6CGPkjGyDlGDc0g9nH3w7Ps4lcC5JTvsSNPgCrW2lhAjiQlpkJLnlziElxYhssG0lF6pGOUMUYMJaaCS0xKT0nTP3jxmh2mMaSmdJcVIIkijEaX6EKK8iNI0JWKUehn3UARpAqJUB6aicpiWxikVo9II0kQlFztgP0LkA1Hqb6kxkghSEiXOKzEl7cZU5A/Tkk9eYiQRpKSIkg+EyDeitF3u6chTjCSClByH8MphKqoH01J53mIkEaTqMSWtEKI6EaXnck5HHmMkESQTTEn5MBXVj2mpbtdfuZ/svghSS/gm3Y4tZ5SWOiURonkhTHmknI5SxkgiSC84+PizZPfFpGSDqWjelhilXIfrPMdIIkidUkYJaRGiZWBa8s0iRhJBMseUlAZT0TIRpXRSTUdWMZII0la1TUlzPY9EiMC0tBwEaYdUUWJKGocQoY0ojVfDdCQRpL1qm5TmgKkI2zAtlWMdI4kgZcOU1A8hQh9EaZ4IUg9MSfaYijAU01I/KQ7X5ZiOJIIEBwgRpiBK80GQekoxJVketnt6MZjdtxWmIqTCtGQn13QkESQUQohggTDVjSAhK6Yi5ECU6kSQBmBxwzSECDnVNC09fHze5H69XmZiG4I0A97PHzEVoaRaouRRzvNHEkGCMUIED2qalpaMIA3EYbt+mIrgEVHyjSBVztvhOkIE77xOS1bnkWpCkJAMIUJNPEYptbufXym9CYMQpIp5mY6YilArb9PS0qckgoRJCBHmwFOUlowgZfTV5dPJ7qv0dMRUhLnxMi0teUoiSBiMEGHOPERpqQhShUpNR0xFWIrS09JSpySChF4IEZaodJhSqGmlHUHKJNX5o9zTEVMRUOYw3hKnJII00PHbbxR77BIxArBSYlpaWpQIUgYpV9flwFQEbFfjIbxaDtsRpErkmo4IEbBfzmlpSVMSQapAjhgxFQHD1RSlMVPSna8PJz/uEARpgDHnj6YerrOOESECpsk1LS1hUiJIjuWIEYA0ajy35A1BMuR1MQNTEWDDelqaOiV5X9xAkHrKvdzbajoiRIC9OUUp53kkgmRkynRkESOmIiAvy2lprueTCJIBjzECUIbHKHmdkghSD7kO16WOEVMR4IPVtDS3SYkgJTZ2OrKIEQBfLMI0NkoeFzgQpD2GTEceYsRUBPhXY5RyHLZ7yfwRsFPqGNXmrQtf9L7t3A5PYNmaKKV63j58fH7Q88mjEGMsvQ1unP3hT+I/vXb9+4+tp6NUMfIaotxPjrkHa8kvvHz5Uf/9VKkLWE6R8jk85nl37fVPet/2+iv3X/jclTf/76MY4zuDH3gDE1ICpWLkJURefivr2o65R2pOhkRn6v14i1bKaWnMpHT38yuDomSFIG3RdzpaYoy8BKiP9rYSJz9SxSfl43uI1OlPD4pFqa87Xx92TkkpcMiupTlkZ3mobuoPfYkQ1RSgvmqMU+2H7EpHaIySkUrxXB/63B176I5Ddg7MOUZzjFBb8/XVGKaa1Bihts3tzxmoFNPS0Emp9KE7JqSWsz/8SXznn/+t121zxihHiOYeoD68x6mmCan2EPWRM05T9wFWk1IzJaWakHgdUks80+8JP6cYvXXhC2K0xvdhupcfxUXESHr+teb4eqf+MjL0l62+r09K/dokDtkNlCtGliFix7sdh/LGWUqEtml//VaT09SVeDUcviNIA9QcIyI0DGHqZ+kh6mIdpynnlixW362mpP9Kcl8csjPw9GIY9YOY+m1/msNxxGg8vnfbEaP9rA7rTXlPvCG/ZOV+vzuC1FPf6cjDVESE0uL7edKSzhOl5ClMXqNEkHqwjFGqqYhpyB7fW6aiFCymprlEiSDtYR2jqYhQXkv9XjMV2Uj5fR0zLT18fL53mHJEiSDtYBWjFFMRISpnad93QmQv5dRkOS1ZR4lVdltYxKjEW4FYSb0c1OPFwnZ568IXi1iBR4zya77nU1bojVki3ncFnuVycILUoU+MxkxFY5WKUM7XIOx6LK+xmnuUaozR2QfPkt3X2AtuppIqTDVFiSBtSB2jGkLk4W3nd+naPi+RmmuUPMcoZXTGPk7OWE0N09BpqWSUCFLLd2f2/wfPESPrEHkPUB/tr6F0nOYWJU8xyhWfobZtl2Wopr7gdsi01Pw879sXpY4SQRqg7w+BtxDNIUC7eIjTXKLkIUZeI9TH5rZbBWrs1GQxLa2ec2neqYEg9WA9FaUO0dwDtEvJOM0lSiXUHKFdrAM1JUypD+GlUG2QQgjnJF2S9CTGeGT1OJZTUcr/yEuO0DbN96T0Ib1a5J6O5hqhXawCNSZMQ6alXFGqNkiSbku6Kum+pJ9bPIBVjFL8hyVA/eUMU61TUs4YLTFE27S/FyniNDZMfaMk2Z7jrjJI6+noquVj9PkPmjtERGiaXGGqLUq5YkSIdksZp6Fh8jItVRkkSb+1umOLqYgQ+cKhvLwI0XCp4mQVJqsoVRekEMJVSe9b3Pd3Pd5xI1eISkWouSSxldRXmJzi2uufmEWplinJcjoiRGmkiNOYMJWIUjVBCiEcSnpP0o1S29A3RmP/I+WIkHVwxj5+qVBZTkveo7SkGB18/Jnp/R+//Ybp/TemxmlImPpMS6l/vl0HKYRwW6uVdEV/ra4xRKXDM1TX9uaMlOW0tCSlQ2QdnjGPaxWrKXEaGiaLq1h3cR0kSddLPrhliFJHqLYA9dH+mnLEyWJa8jolWUxHuWNUKj5DdW1n6kg13/sxYUo1LaXgPUhHkl7t+Pw56we2ilGqEM0xQLvkjBPT0nC5YlRLhPaxitSYqcnTtOQ6SDHGy12fDyHcldGyb4sQpYjQ0gK0S444pYyStykp5XSUI0RzidA+m1/n1EANnZr6hslyWnIdpNziwf4nas4QEaH9mu+RRZiYlHazjNFSIrRLqkBZhil1lEKM5d9McaiNCel+jHH0OzWEEG7o+cq9f5T0PxM3D7Dyd5L+XHojgA7/EGN8eeqdLH5CijHeknRLkkIIH8YY3ym8SUAnfj7hVQjhwxT384MUdwIAwFQECQDgAkE66VbpDQB24OcTXiX52Vz8ogYAgA9MSAAAF3auslu/oemvZPPOCE8k3bS82isAID2rK3bvW/ZtdqmHtQfiuDjwgta721/V+okv6UtJ97T6RY5XTaMkkyt2cw4JcGT9m2fzZN/ljqRfxxif2G8V8Nz6Z/Sv6w+T7n8X/8JYwIv1E/0jrSaifa6vb8cvY8jN7IrdLGqYKIRwLoRwGELosxMBdvmjTsboSNK7ki6v/7y3/lzjcH3NMCALyyt2SwQphdta/VbLjgGjhRCu6+SFKO/FGC/HGO/EGI/Wf25pNRG1zx9dX59vAsysf+m+Kemu5eNwyG6C9SEWk8tgYHE+aP39iVaT0QtijE9CCO9qtSCo/W+vGW4bFqjEFbuZkKYxO5aK5Vgf7m0fqvvdrsUK62W27dWpV9e/HAEpbU7t5gjSSNbHUrEo72183OelEDc3Pv5lom0BGkdaTeubf8wQpIFyHUvForR/Cz3qs5S743VIrLZDUutzmD/e/KPVa+FMcA6phxLHUrEo7WscDXnB6309/5nkOkmoXpVBijHmPoF7PfPjYVna53+GvA3LkZ4HiZcdoHpVBqmAI0mvdnyeE8mYpOP1aw86b9itHS9+FlE9gtRDjPFy1+c33sIIGGMzJF8O+Ld/aX8QQjjHWwmhZixqAHwhKFgsggSU1XUouK/NeE25L6A4ggTUi/NGmBWCBJS1ec5oSmSGnH8C3CFIQFlTDru91v6ABQ2oHUECCuq4/HPnis4teO0RZoUgAeW1J5shkWnflkuao3oECSjvw9bfh7w9Vfu2H269FVCJKl4Yu74A2a9ks6roiaSbHYdOgFzaL7C+1OcFrh3v8MCb/aJ6VQRJqyer5aUeHqjfW/4DFjbfPfmX2v/zuHnJCrN3YAZyqeKQXYzxP2OMwfAPMUIx60tJtCf03+y6/fpifDdan7rHCjvMQRVBAhagHaFLIYQPtt5S+r1OHr7eGTCgFgQJcCDGeEcnV8q9H0K42T5XtL445F2dvBzKnY6L9QFVquUcErAEv5D0Jz2ffm5IuhFC2Hb7+zHGd3NsGJADExLgxPo80M/Ub4HCnRgjly3HrDAhAY6so3QthHBVq5V0h3r+AtgjrWJ1k8N0KMXyit0ECXAoxnhPLOXGwnDIDgDgAkECALgQYoyltwEAACYkAIAPBAkA4AJBAgC4QJAAAC4QJACACwQJAOACQQIAuECQAAAuECQAgAsECQDgwv8DVVi+CilPQvAAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"x=np.linspace(-1,1,101); # 10 data points from -1 to 1\n",
"y=np.linspace(-1,1,101); # 10 data points from -1 to 1\n",
"X,Y=np.meshgrid(x,y); # 100 data points from -1 to 1 on the x- and y-axes\n",
"Z=f(X,Y);# Z=f(X,Y) evaluated at the 100 data points\n",
"\n",
"plt.contourf(X,Y,Z) # also pcolor, countour, ..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**For more plotting examples and options check out the [Matplotlib gallery](https://matplotlib.org/gallery/index.html)**"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Here we will save a function called `my_ode` as `my_ode.m` For a single-DOF forced spring-mass-damper:\n",
"\n",
"$\\ddot{x}+2\\dot{x}+9x=\\cos(2t)$"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"def my_ode(t,r):\n",
" \"\"\" Help documentation for \"my_ode\"\n",
" input is time, t (s) and r=[position (m); velocity (m/s)] and time\n",
" output is dr=[velocity (m/s); acceleration (m/s/s)] at time, t\n",
" the ODE is defined by:\n",
" \n",
" a = -2*v-9*x+cos(t/2)\"\"\"\n",
" \n",
" dr=np.zeros(np.size(r))\n",
" dr[0]=r[1]\n",
" dr[1]=-2*r[1]-9*r[0]+cos(2*t);\n",
" return dr"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given any `[position, velocity]` and `time`, the defined `my_ode` returns the velocity and acceleration:"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([ 2., -12.])"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"my_ode(0,[1, 2])"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## Euler approximation\n",
"\n",
"The simplest integration routine is the Euler approximation. A first-order Taylor series expansion about the current timestep:\n",
"\n",
"$r_{i+1}=r_{i}+\\frac{dr_{i}}{dt}\\Delta t$\n",
"\n",
"We can integrate the equation as such"
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"dt=0.01\n",
"t=np.arange(0,6*pi/2,dt) # (start, stop, by some incr)\n",
"x0=0.1\n",
"v0=0\n",
"r=np.zeros((len(t),2))\n",
"r[0,:]=np.array([x0,v0])\n",
"for i in range(1,len(t)):\n",
" dr=my_ode(t[i-1],r[i-1,:])\n",
" r[i,:]=r[i-1,:]+dr*dt\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 68,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'position (m)')"
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAAE9CAYAAADanxEeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO29eXQc53nm+7yNfV8IruAKcN9EgdDqRbJFOrbjcRKblLJP5toi5dw5NzkniSjlnJkzd86MZdKZLXOvbVKZJDfJJJZIOY4dO5YJWt4kayFAiqS4AAS4ggtA7PvS/d4/qqr7q2Kj0d1V3bX0+zunD7qB6uoP3UA937sTM0MQBEEQBHcJub0AQRAEQRBEkAVBEATBE4ggC4IgCIIHEEEWBEEQBA8ggiwIgiAIHkAEWRAEQRA8QL7bCxD8SV1dHa9evdrtZQiCIPiK1tbWe8y8MN7PRJCFtFi9ejVOnjzp9jIEQRB8BRFdm+tn4rIWBEEQBA8ggiwIgiAIHkAEWRAEQRA8gAiyIAiCIHgAEWRBEARB8AAiyIIgCILgAaTsSRCEpBgcn8b/fuc63r3Sj4K8ED66vg5PN69AcUGe20sThEAggiwIwry8e6Ufz/1dK/rHpqPfa7lwF3/7i2s48rvNWFNX5uLqcpuekUn85FIvxqZmsaW+Cs2rakBEbi9LSIOMCzIRVQKoBdDPzMOZfj1BEJzl9I1B/M7/egdTs5H7ftbRM4rffPltHPvS46ivLnFhdblLJML47yc68I2fdGJa+WyaVlbjz/Y+gIaF5S6uTkgHx2LIRFRJRJ8joq8T0XtE1EdEYQADADoBDBBRWP9+BxG9TkR/TEQ7nFqDIAjOMjg+jX1/czIqxnXlhTj4+W34009vRHGBdvm4PTSJP/iHU5gN3y/YQmaIRBh/+Mpp/PmJDpMYA0Db9UHs+cYvcOG22D9+w5aFrIvpfgC7ADRYfzzH02oAVOvH79LPAwAtAF5l5v9lZ02CIDjHl79/AT0jUwCAqpICHH3u8ah7evPSKvzrv3oX4Qjj5LUB/NWbV/HsR62XASET/PeWdnzn/VvRx9vqq7BuUTm+e+YWZsKM/rFpfPH/O4nv/V8fRnVpoYsrFVIhLQtZt4TfA9AKYB80cSXLLeEp4tx2ATiiW9FfI6LV6axNEARnONc9hFdP3ow+/uqe7aZY8YfX1eEPn1oXffznP+owxZiFzHDq+gD+5xuXo49/85GV+Pb/+SH812d24JX9j6GiSLOzugcn8O/+6QO3limkQUqCrLuY+wAcBdAEs6AOQrNyDwHYC2A3gJ0AGqFZxTX6/d36z/cDOAKgyzi9ctsPoJOIfkBED9j4/QRBSJM/P9ERvb9782J8YsuS+4557slGNOgiPTI5i8M/7cza+nKRSITxp/94Dsza48cbF+A/fnYL8kKaDdS0sgZ/9nTskvnd92/hzcv33FiqkAZJCTIRfZyIOgAchCashnC2QBPPRmauZeZPMPMLzPwaM59g5lPMfIWZh/TbFf37rzHzy8z8HDOvZeYQNKE+BE2gjfPvBtBGRF9z/lcXBGEuzt8axg/P340+/qNPrI97XEFeCM9/cmP08d+/fR0jkzMZX1+u8r2zt6Ox4ZKCPBz8/Hbk55kv47+0ZQl+Zcey6ONDr18CGwoueJqEgqwnar0C4Dg065YAtEET4RpdgF9m5it2F6IL9QvMvBaaZf0ygCH9NfcT0T0i+pjd1xEEYX7++q3Yv/Snti7BxiWVcx77ic2L0bBQt5KnZvHNd29kfH25SDjC+B+K1+ILH16DFbWlcY/9009vQmG+dnl//8YgftYhVrIfmM9CvgJgDzRRPAbNEm7WRXgoU4vSLev9zFwL4DlowlwLoIWIvpip1xUEARienMF3378dffzFjyRO1AqFCM8qx/zN21cRiYhF5jTHz9/F5Z5RAEBFUT6++JE1cx67uLIYv/HQiujj/+dHl+c8VvAO8wlyDYATAHYy89NOWMKpwsxHdGF+AZowV2d7DYKQS3z7VDcmZsIAgI1LKtC0cv5/uV97sB5VJQUAgBv9E3j3an9G15iL/N3bsbn2v/PYqnmzp/c/0YiCPC22/O7Vfly8I2VQXmc+Qd6pu6VPZWU1CWDmQ7owv+b2WgQhyBxVMqt/65GVSXV9Ki7Iw2cfiMUt1XMI9rlybww/15OzQgT81qOr5n3OsuoS/JKSiPcP71zP2PoEZ0goyF4QYituWOmCkCtcvTeGs91aNKowL4TP7qhP+rl7di6P3v/+2dsYm5p1fH25yjffjYnpxzcuSror2m8+vDJ6/1unujExHXZ8bYJzyLQnQRCifO9sLHb80fV1UTd0MmxfrjWnAICJmTB+0t7r+PpykUiE8V2lCchvKCI7H482LMDqBVri18jkLFou3J3nGYKbiCALghBFvfB/ZvuyBEfeDxHhU9uWRh//y7k7jq0rl2m7PoBbQ5MAgOrSAnxk3cKknxsKET7XFPNc/POZWwmOFtzGseES+hCJZthMumLmbzmzIkEQUuFG/zgu3hkBABTmh7Br8+KUz/GprUuiDUV+dOEuJmfCMp7RJv98Jua1+OSWJdFypmT5zPal+K/H2wEAb1zqxcjkDCqKk/d8CNnDtiDrZUgH4Uz2M0NGQgqCK/zoYk/0/ocaF6C8KPV/xY1LKrBqQSmu9Y1jbDqMNy/fw1ObUhd2QSMcYVMY4V89kJrXAgAaFpZj89JKnL89jOnZCI6fv2uymgXvYMtlTURfB3AY5u5daZ/O5vMFQbDBG5digvyxjYvSOgcRmTJ7f3xJ4sh2aL02gF59uEddeREebViQ1nk+80AslHD8vMSRvUragkxED0Lr2AVolq3RCSDe4IhkboIguMTEdBi/6OyLPv7YhvQEGQCe3BCLcf64vUfaNtpA9Vrs3rw42rM6VXYrXoqfddy7b2Sj4A3sWMgvKvcJwClovadrmDmU5k2CTYLgAr/ouhedebxuUfmcLRmToXlVLUoLtX/lG/0TuNo37sgac5EfK16Lj6fptQCAtYvKsaJWK5UanZrFSWnc4knsCHITYlbxcb2l5olMttQUBCEzvHEx5lq2c+EHtISwxxtjrtWfKKIiJE/34EQsyS7P/J6mChHh44rX48RF+Uy8iB1BNmYgAzHXtSAIPuStztjwgSfWJ19WMxfqOX4qgw3SQrWOH2moRVkaSXYqal7AGyLInsSOIA8ad5j5qv2lCILgBj3Dk+jsHQOgWWJNq2psn/OJ9bGL/1udErNMB9VrYSemb/BowwKU6CVoXffGcF1CCSkzPRvJaE6EHUGWIIQgBIBfdMWSuR5cWe1I3fDKBaVYqcehJ2ciOHNzcJ5nCCqz4QjeVj6XdLPeVYoL8vBoQ230seoVEZLjL9+8gke+fAL/9u/b8NZl598/O4IcHfIgc4r9ARFVE1ETESWepyfkFOqF/zEbcUorj6yJXfzfuSL791Q4d2sYo3ov8GVVxdH2l3b50Nq66P23lKx6ITne6epDz8gU/vnMbdwcmHD8/HYE+RvK/cN2FxI0dOE7TESdRMRENKDfP0xETS4t6yiAVv2rIAAA3u6KieVjada5xuMR5Vyq6Avzo5agPdqwIKmJW8mgbrje6uyTkrQUCEcYJ68ORB8/rGw4nSJtQdanLj0HLbGrkYh+QETzzwQLOLoVehya8O2DlvwGaJ3MGvTvtRLRUSLK2mxn/bV2Zev1BH9wZ2gSV+5p8eOi/BB2JDH7OFlUC7n12gBmwhJHThZ1A/Oog16LTUsqUV2qtc28NzqFyz2jjp076Jy/NYwR3WuxuLIIqxzyWqjYSttj5iP6hf4r0GqQu4joCDQxOpnmOU/bWZOb6O9FK2IinIg9+nE7M7qoGC/Of4iQa7yr1KM2raxBUb5zrQBW1JaivroE3YMTGJ8O42z3EJpW2k8YCzoz4Qjeu5oZr0UoRHisYUF08MdbnX1Yt7jCsfMHmXeuxDZJj6xxzmuh4kTf6E7lPkGzANPF772sT8Asxl0ADgBo0x/v0h8bxzQR0VFm3pvJRRHRLgDPZ/I1BH/Sdi3mgntotfNi+UhDLb7V1g0AeKerXwQ5Cc52D2Fcn1tcX11iq0lLPB5rVAX5Hv7146sdPX9QUfMgHmlw3l0N2O9l/RUAryLWIEQNSKTaNtPXLTSJaA+0ZikGLczcyMzHmLlLvx2BZhG3KcftyVRM2YhjAzieifML/qftekyQH3Sg3MnKo2ti1p10h0oONX7sZJJd9JwN6mcyIHHkJIhE2OS1eCQD8WPAfi/r52EW0XR7U/tWiBUOKvcHAcS1epk53s8Oxjs2HfTYdCsRMWJxbEG4j8mZMM7fGo4+blrhvCCrNc2nbgzKxT8JTPFjB93VBo0Ly1FVosWR+8ampbVpErT3jGBwfAYAUFdeiMaF5Rl5HSd6WTM0Qb0CrWPXTgCN0CZApXrLzLYjw+hlRKqr+iVdeOPCzF0Ajijf2uVggpfVUheEuJy5OYTZiCaQaxeVo6rU+Rm5DXVlqCzWolD9Y9O43i8X/0REIozT12OXjkxYYqEQoUlJ3mtVwhZCfN5RKhEeXlObkfgx4Fwv61ZmXsvMLzPzKWa+wsxD6dwc+J3cwNo69Ejco8xYS8WedmgtXdAsdOtNEEyoF+ImB7OrVUIhwoNK3Fh1kQv3c7l3NJrJW1dehOU1JRl5nebVaga8hBLmQ01+fHh15uxGp3pZP+vAWvyMapF2JbKODZi5zfItR7Kt9bh1jfUGoMWJ8wvBQRXHTCZbPaiI/anrsjdMRJtlk5QpS2ynEkoQC3l+VK9Fs0cFWe1l7dtSJYdoVu5bhTYR6rHNcx4lCA7DzKaL/84MJHQZqGIvgpwY9f1xoqf4XDywvBr5+mzl9rujGNLjo8L99IxMontQ68pVXBDCxiWZKxOzI8hdjq3C/6j+vlTeF/VYaWdpgZkRibAkAmWA6/3j6BubBgBUFudnLEkFAB5YEfv3uHB7GBN6SY9wP6as9xWZ6xtUUpiHLcsq476uYEa1jrfVVyE/z1ZxUkLs1Py+Ct1VS0QfY+Y3nFmSv4jTF7oz7oHxUQU5a127vMjY1Cwu3R3BxdsjuHhnGBfvjODi7WEMT86CCFhaWYxty6vwS1uW4NPbljoyACGXOXMzlq7xwIpqhEKZK3SoKinAukXl6OgZxWyEcbZ7KCNtB/3O0MQMOvTOWXkhwvblmb0kNK2qwfv630HrtQFHBlgEkVM3YoL8YIbr6O0I8mEALwCogla287AjK/If1v+aVDIkTA1+iag6mfizXxmdmsU/v38LY9NhjE/Non98GlfujaGzdxQ3+udu1M4M3BqaxK2hSbz+wV189fVLeOFTG/ErO+qzuPpgcbY7Jsjb6qsy/noPrqyOik3b9QER5Di8r1z4Ny2tQElhZjedzatq8VdvXgUgFnIiVAt5Rwa9FoANQWbmISJ6GsAPAewkoh8A2MvMI46tzp8EVlDtMjo5ixe+ddb2eW4PTeIPvnkaP+u4hy//2jYU5mfOhRRUzt7MtiDX4NWTN+97bSGGKX6chY5mD6yIfe5nu4cQiXBGPSV+JBxh0+hQzwoyADBzCxHthtYJajeAQSI6CK2PdRvSmJnMzMPzH+Up7Gz1reJdG+d7noGI9kFvNLJy5cqUn19aNPeOPy9EaFxYho1LKrFxaQU2LqnAxiWVWFpVjJkw48q9MbRcuIu/evMK7o1qsc9jrTcxMDaNr/12k6M9mIMOM+PcrZgobs2CIKuir1rnQgxT/DhDZWgq9dUlWFBWiL6xaYxMzuJq3xgaMphL4Ecu94xiTM95WFRRhKVVxRl9PVuCTEQd+t1BxFy3B2yc0u+9rFPFV3FjvfXnEQBobm5OOdOqtCAPTzcvR2lhPsqK8lBRXIBVtaVYs7AMa+rK5hTVwnzChiUV2LCkAr/3+Gr8u386F+2PfOJiDw4cO4P/9syOjJWIBI1rfeMYmdRqXWtKCzJW66qyfnEFCvNCmA5HcL1/HEPjMxlpROJXIhHG6RvZtZCJCNuWV+HHl3oBaBslEWQzpyybpExfY+yKXyPm7mOdK1i9AHZENtAV+vl5IRza84Ctc5QV5eO/7H0AS6uK8f++oeXPffv0LWxbXo0vfHiNE8sMPKqFurW+KisbmcL8EDYsqYi+9rlbQ/jQ2rqMv65fuN4/jqEJrfSotqwQKx0eKDEX25dXRwX5/RtDkpdhQd0k7chAa1krTgTf5hoWkSvEczsni6lRbZATupyEiPDHn9iAX39oRfR7B//lIjru5nr6QnKcy3JCl8HW+liZjbitzbixSQKA7aZQglx+rJgFOfMOTbsWcrZm+XoWZu6y/PM0pvB0qT1OEyLC//0rW3D+9jDO3BzCdDiC5187g2PPPY48SUxJSLYzrA20WPWN+9YgwBzTV+qDM8325bHP/1z3MGbDkYzW2fqJiekw2vVNPpH5vcoUdpO6Tjm1EJ+jxtBTEVn12FQ6fAkAivLz8Gd7H8Av//nPMBNmnLo+iH881Y09O5e7vTTPwswmCzkbCV0Gqvh/IIJs4oPuWC5rNjdJiyqLsaSyGHeGJzExE8bl3lFsXJK9DYGXuXhnGPrsFTTUlaGsKPPpTbIVcoaTyv1UJi2px56c8yhhTtYvrsCXnog5Jf7LDy9hckY6Qc3F9f5xDGc5octgw5KKaLvGq33jGJ6Udo2AtkmyuqyzyTbF8jsjJWlRzimjSbP1mYggO8Nx5X5DMqMU43T4Oh73QGFe9j3RiLryIgBajfLf/uKayyvyLm7FKgHNo7F+cawP8DmxkgEANwcmogldVSXZ3SQBwAMmQZY4ssF5JYywJUthhISCTEQfz8oqkoSIKoloh9vriIN1klIyoxStIxtlGlOalBfl4w92rYs+/oufd2FqVqzkeHzgwq5fRXXHiiBrmEMIlVkv39umtOiUpi0xzilhhK3LvGEhtxDRu0T0saysZg50If4KgAEAu9xcSzz0UYpqX+qEtdi6Bb1P+VbLXBnWRFRtvdlfcfB4unk5FlVoVvLd4Sl8+1S3yyvyJhdvxy4ym5dmP1a4dbma1eu3HkCZwZzQlf1NkpppfeH2CGbCkayvwWvMhCO4dCdWtbHZCxYygCFoYwFbiOgHRGSviDRFdCH+Y2hC/CfQSqq8OmVKFeEGvWPZXLwMc71yIgF/Edrvb9wSnTdnKcrPM9Uhv/yzKzIlKg4XlYvMpqWZGyM3F2oGsSR2aagbky0ueC1qygpRX625yafDEXT2jmZ9DV6j4+4opvWNyfKaElSXFmbldRMKsj7Y/mVoQrgbQBsRdRDRHxHR6kwtiog+R0SvICZABOAKgJ3M/K1Mva4dmPkYzJnSzxPRYTVWTERNRHQcwB7luGO6hS3Y5DcfWYkyvSH/5Z5RvHdVGuarDI5P4/bQJACtUcfqBWVZX8PGJZUwqtKu9o3l/ChGZjZtTLKZYa2ySfGWnL8lnotzLsSPgSSSuph5PzQr+So0YWwEcAhAJxG9R0R/bDeuS0SrdRH+OhH1ATgKTbSMRiMHmXmtD8qsnoK5Ucg+aO8TExEDaIXZ5d7GzHuzucAgU1FcgM8qnYb+4d3rLq7Ge1y4HbOO1y8ud6XetKQwD6vrtI1AhBGt88xVbg9NRudSlxflY1WWOnRZUV2yIsjm9yCbYYSkCqt0C65RHy7wPGL1s02IzUQGNAvxJLSZwF3QxKkfMZFqgOaqrdW/7oYm9qr7Vs1oOAJNjK+k+Hu5AjMPEtEaaBuK+WLdx0SMnee3HlkZFeLvnb2Nf/+Zzagpy467yetcvBO7yGxysdZ009JKdPWOAQAu3B7GA1nogORV1ISuLcsqXZu2pOYTnL8tguxWrX5Klc7GcAEi2gUt7vmU5ZCoQKeA9S9wEMBLAI4ws++CTHpy1m79PdoP7f0wNjBd0LKpDyfrpmbmA7A3sMM4z2675/ADW+ursK2+Cme7hzA9G8G3T3fj33xIelwDwEXFQt7oQkKXwaYlFfjemdsAzDHtXMStC78V1S174fYwmDlnh7VEImzalHjKZR0PZm7RL/CN0MTCEJd4fa2tPa7jfW8QmjW8l5lrmfmrfhRjFf092svMjcxM+q2RmfdLzDizPKP0uP7O+7dcXIm3MFvI2U/oMlA7QV3IcWvsvLJJyuaF38rymhJU6J2oBsZncGd40rW1uM2VvjGM67kNCyuKsKgysyMXVWwFkZj5ii6ezcwcguaCPgTNCmzD/YMXjCzpNgDHoFmQjboIP8fMr9lZjyAAwKe3LY12hDp1fRA3+sddXpH7hCOMS0q8doOLgrwpjjWWq3T0eOMzISJJ7NKxhhGyiaNZHcx8gplfYOZP6CJdqwt1DTOH9Nta/WdPM/PLfokPC/6htqwQH14XG+333TNiJV/tG8PkjFbGsaiiCAv0zmZusKyqGJXFmjU2PDkbzfzONSamw7iubxZDBDS6PItYErs0LrjotchKmqXf3c+C//jsA8ui979zWgTZK/FjQLPG1DXkqtv6cs8oDOfA6gVlKC7Ic3U9ktiloWb+Z3vQhvSyFgLJ7s2LUZSv/XlfvDOCK/fGXF6Ru3glfhxvDbma2KWGENQe326xeZlskgCzIGf7cxFBFgJJRXEBPqK4rU9cuOviatxHdcNtctlCtq4hVy/+HXfNdeFus3ZReXSW+NW+cYxOzbq8ouwzNjWLmwMTAID8EGFNXXab54ggC4Fl9+bF0fvHz+e2IKsW8kYXWmZaEZe1xUL2gNeiuCAPa5U49sUc/Fw6emJtQ9fUlaEwP7sSKYIsBJaPb1wMo5Ty5LUBDOgdkXKN4cmZ6K6/II/QUOe+NbZhcUX0s7lybywnZ1h33I1d/L3gsgYsiV05KMhuuqsBEWQhwCysKMIOvQtUOMJ441KPyytyB3VqTePC8qzv+uNRUpiHNQtyt4XmyOQMugdjmyQ3+orHY3OOey7alf+VdS6EEdz/zxSEDLJrU8xtfeJCbgqymjS10QOuUQM1jqxmgecC7XfddY3OhVoLfSkHk+3aFZf1BrGQBcFZntq0KHr/zc57CEdyrwlFp3KR8UKs0sB08c8xC7nDYxnWBupn0n53NOeatqifyzoRZEFwlg2LK7CwQmuCMTg+Y+rCkyuo3aDWLfLOxV/NLM41l7Wpa5qHBHlRRRGqSgoAAKNTs7iVQ01bhiZmYuNJ80JYvSD7k7dEkIVAQ0T4yNpY+dPPL99zcTXuoCYPrVvkfkKXgWqBqGvMBUyfiYcEmYjMVnIOua0vKxvXhoVlrownFUEWAo/aRvOn7b0uriT7DI3PoGdkCgBQmB/CCpfm7cZjVW0pCvWL3p3hSQxNzLi8ouzhlb7i8VAt9lxq2tLugax3EWQh8HxYsZDbrg9gLIcaHlzuVXb9dWXRxg9eID8vhEbFYu/IEbf1wNg0evVNUlF+CCs9tEkCzHkGuRRKUJPY3NokiSALgWdRZXE0u3gmzHjnSp/LK8oeXnWNGpjjyLnhtlZFTu2O5RU25mimtTnXwp3QjgiykBOobTR/1pE7ceTLPd6MHxuorsFcscbU0hovZVgbrFcS/y73jmI2HHFxNdnj0h2l5EksZEHIHI8rbut3r/S7uJLs0uFxQVbXlDOCfMebJU8GVaUFWFJZDACYno3gal/w54kPjE3j3qgWRiguCGFFjTthhHynTkRElQB22T0PM3/LgeUIgonmVTUIkdYV6vztYQxPzqCyuMDtZWUck4XsgQEGVswWcu65rL0wVCIe65dU4M6wVgLUfncEaz24mXMSaxgh5FIYwbYgE9HHARwFUG1/OWA4uEkQBIOK4gJsWVaFs91DYAZarw7gYxsXzf9EHzM6NRttz5gfIqzySHtGlRW1pSguCGFyJoJ7o1PoH5tGbVmh28vKGMzser/kZNi4pCJakXDxzgg+vW2pyyvKLF75TGy5rInoWQDHAdQAIP0G5X6yN1juC4LjPLymNnr/nRxwW6sdulbXlaHAhbrK+cgLkcn6Crrb+t7oNAbGtfKu0sI81FeXuLyi+Jg8FzmQ2OWFkifAhiATURWAw9BElPUb9MeDKd6GlPuCkBFUQX43BzKtvR4/NlCTiIJe+tRuac3olmt0PjbkWLLdJY+EEey4h1/UvzI0EW4D8Cwzn7K9KkHIAA+tjgnymZtDmJgOo6Qwz8UVZRavZ1gbrMuhOLLJNerpz6QcRAAzcLVPG49ZXBDM/xVm9kxvcTs+rCblfiszN4sYC16mtqwwuvudjTBOXR9weUWZRW0FuNajsUoA2LAkd1zW7R7u0KVSXJAXHQkZYfPmLmioYYQyl8MIdgS5Wbn/rN2FCEI2eGTNguj9oMeR/eKyVgdetN8dCfSEoXaPN2pRUV23QW4QYg0jELkXRrAjyNGsamY+7cBaBCHjmBO7ghtHnpwJ43q/Vj8aIm3mrlepry5BqR46GBifwb3RaZdXlBmY2ZQg5aUpT/HYsCQ2rzrIngsvlaHZEeQux1YhCFlCFeTTNwYD24Wos3cUhqG5akGZp+N/oRCZLPigJnbdGZ7EiN5HvaI4H4sri1xeUWJyZciEV0qeAIcEWW8KIgieZ3FlcTRGNDkTCeyFRo35NS70rrvaQL0QXgqoIF+yWMduukaTQY3tB3WTBHin5AmwJ8iHlPsSQxZ8w46VsR42QU3s8nqHLiu50LHL64M+rKxaUIaCPG3TcGtoEsOTwRuP6bVGLWkLMjO3ADgBreTpEBE94NiqBCGDNK2sid4/dT2Ype+mi7+HE7oM1E1DUK0x0wxkH2ySCvJCJu9KRwA3SneGJzEy6Z0wgt3WPXsBXIEmyj8ioi/YX5IgZJYHVQv5RkAF2TRKzvvWmHXqUxAzrU21rh4ueVJZF/AGIao3xgthBFuCzMyDzNwI4Bi09plHiKiPiF4ios8R0Q4iWk1ElcneHPmtBCEBW5ZVolBvI3nl3hgGxoKV1Wud0NO4yLsZ1gZLq4pRUaT1KRqenEXPyJTLK3KWSIQ9FatMlg2Lg10j3mEpeXIbW4MciKjD+i1owvx8mqeU4RJCxinKz8PmZZU4rVvHp28MBmrQxNW+MYQjmoW5vKYEpYXe/5ciIqxdXB4NIbTfHcFifQRgELg5MIGJmTAAYEFZIerKvZ1hbeH1U5EAACAASURBVBB0C9mcaOd+GMGuy7oRQIN+Y9zf0zqdWyAgoiYiOkxEnUTERDSg3z9MRE3znyEYa/Aqqtu6LWCJXX6LHxuYe1oHK15pbj7ho88k4Ml27T3e8lo4Mf4l0MKaKkRUTUTHAbQC2AdtswJojVQa9O+1EtFRInJiZKUn1+B1HgxwYpeaYe2nObamxK6eYFlj5oQu9y/8ybKythRF+ZpM9I5MBSq8E4lwsFzWAHY6soqAoItbK2ICmIg9+nGOvodeWIMfeHBFbB9y+sYgwhFGnkcn76SK3xK6DNQLYtAsZK9d+JPFGI/5wa1hAJql/0jDgnme5Q+6BycwPq2FEWrLClFX7v4cbluCLMMk7uMEzELYBeAAtElYALBLf2wc00RER5l5b8DW4HmW15SgrrwI90anMDo1i87eUU+4rJzAZCH7yD26zjIXmZldz3p1iktqNq9PMqwNNiyuiAlyz2hgBNm8cS33xN+a9yaW+xQi2gPzBKwWZm5k5mPM3KXfjkCzRtuU4/Y4Fc/1whr8AhGhSY0jXwtGHHk2HEFX71j0sZ9c1kurilGuZFr3BiTTejYcQWevEqv0kdcCsCR2Baiz3aU73oofAyLITnJQuT8IrUb7Ppg53s8OxjvWp2vwDWrHrjPdQy6uxDmu949jWu/PvaSyGJXFBS6vKHmIyLSBCEoS0bX+cUzPap/J4soiVJX65zMBgjse04t14RkVZL0GeQcRfVytSc7ka7oBERmZ5gYv6aIXF2buAnBE+dYuu8lVXliD39heH/t1z94MhiB3+KxlppX1AUzs6vBQa8Z0COp4TDXRbr1HPEmOCrIuul8nog4iCgPohJZgZGT8dgIY0H/+dSLa4eTru8h+y+MjcY8yc9jy+OkArMFXbKuvit6/eGcYU7NhF1fjDH7NsDYIYpmNF12jqRDE8ZjhCJv+V7zyuTgiyLrl+zpiZTaNSFxr3IhY6c27RLTKiXW4iBp/7UpkmRowc5vlW3Yznb2wBl9RVVqAlbWlAICZMKP9jv8FwO+CvDaAYxjbe1QL2X+fSShElgx4/38uN/rHMaWHERZWFKGmzP0Ma8ABQSaip6BZvruQev0xAWgG0EVEH7O7FhdpVu5bRS4R6rHNcx7lnzX4jm3LY1bymW7/1yP7teTJQLVUOnpGA+EeVROhvGKJpYrq0g3CeEyTu9pDmyRbgkxED0JzR6tCTNAu8oeguVF3Q7O8duuPD+k/N57D+v0WH0+MUmOvXXMedT/qscnUDXt9Db5ju+K29nscOWJxw/mpS5eBmmk9NDHj+0zr6dkIrtyLZb37qQZZRS3VCkIowatxfbuNQY7qX6OiCuBAgvrkE8YdvczmIICnlOcfBbDe5pqyip5MpdKZwtNVMUw7ocoLa/ArJgvZ54LcPTiByRnNDVdXXugZN1wqGJnWRp/xjp5RLPJxT+sr98Ywq/cVr68uiW42/EbQelpf8uigj7QtZCL6PGI9rAHgIDN/ItlmIczcxsy7AXwVMWu5kYh+Ld01uYRVxPpTeG6f+sBGlrMX1uBLtioWcvvdEUzO+DexS3VX+zF+bGBtEOJn2j3qGk2V9Yvvb9riZ7xqIdtxWT+j3G9j5hfTOQkzq12kAODXbazJC3ghEOmFNfiCyuICrKnTxhPORhgXfdz4wOyu9s5FJlWscWQ/0+7BWtd0WFJZjIpizbofmZzFneFJl1eUPjOW5jleKg+0I8hqVu9LNtfxlTnO6wdqbTzXKpzpnisrayCifUR0kohO9vb22nhJb7HNFEf2717GNOXJQxeZVDENmQiShezjTRIRBaYk7VrfWLR5ztIqbzXPsSPI6oU7lazeeLTqXwm5lVjkBfdw0mtg5iPM3MzMzQsXLszkmrLK9oDEkVVrcu1CPwuy+cLvZ/dou497WFtZH5DSp3aPxo8Be4KsXshTiVnGw+7z3cS6djsim+774IU1+BaThezTFprM7NuhElaWVRWjTG9EMTQxg95Rf2ZaT86Eca1Pc40S+TuuD5jjyJd8HNrxclzfjiCr2bl261dVq9hvPkM7bmfT2JRkmnl4eA2+ZUt9FYxBLx09o5iY9l9i153hSYxOzQIAqkoKsLC8yOUVpQ8RYW0ARjFe7hmFnmCNVbWlKC7Ic3dBNlHnOLf7OLbf7uFRmE4J8j6b6zDaPjKAkzbPlVX0ntAqjSk83RH3vBfW4GfKi/LRoCd2hSOM87eHXV5R6pjixx4ZJWeH9QHo2NXR481M3nSxduuKRPwZSjCFETz2udgR5GP6VwKwN91OW3qnr32IlU8dTXC4V1GtylQETj3WbhzeC2vwLduXq4Mm/OckuOzzoRJW1N/Br9aY33tYW6krL0SNPqlqfDqM7sEJl1eUOlOzYVOjFq+FEewI8ivQREDttPWFVE5ARM8C+KFyjkEAr9pYk1uoVn0qWeLqsXY9A15Yg29R48h+HMVoSujycTavgWqNXfapy9qL4/3scH+mtf88F1fujSGsW/YraktQ5rFGLWkLMjMPAXgBmpAagnqEiPqI6Mv65CfTqEUiqtS//xIR9QH4BmIDJxhaly//+Qu19qEGDck014jTXet43AP9tQbfonbsOudDQb5s6mHtrV1/Opgu/D3+bETh1X7JdvB76ZMpw9qDG1dbvayZ+QhinbYMUa4BcABaKdMAEYV1kQ4DGNC//7x+nPE8ADjEzH9hZz0u0mJ5nMwYQ+u4ROs5/LgG37J5aWU0seuyzxK7mNl0ofGaGy4d1EzrQR+O/BubmsXNAc2lmx8iNNT5/zMBzJa+H2P76qAPryV0AQ5Me9I7bT0HwDArDIFVxy3WWB5bj9ufbqcvL6CPMVQTqw4kOl63XtVEuJa5spuJqNp6y/YacoGyonw06rW7EYavErvujU5jaGIGAFBWmIelVf7t/Wxwf6a1vy7+aghhdV0ZCvMdHT3vGn6f+qS62Tcs8d4myZG/Et1SXgPNhX0F849hJGgCfgBADTO/7MQ6XEYVwAYiOpjg2JdhrhVOJJ4vQvMsGLdE583UGnKCrctiERY/ua1NPawXV/g+w9pAdb37rYWmaol5LZPXDqrL+nLPaDQe6xdMJU8edFk7FtHWY8qHABwioipo85EboNW5VkNL2OqDZsW16McHBmY+RkRtiCVJPa9boQeNsiRlwtUu5anHdOs2EGvwM1vrq/Dt07cA+EuQ1Qzr9QFwVxtYBxr4CXOta3A+k5qyQiysKELvyBSmZiO43j8e7QXvdSZnwrjWPw4ACHm0UUtGUsx0sX0tE+f2OE9B8xAYluc+APsSWCxtzLw3gGvwJX7t2BXUi/86HzcHUd25QbKQAW2jZMypbr874htBvtwzCiM3cNWCMk82aglGYMMj6DHYNUguOeoYM+8M4hr8ymbFZd3RM+qbUYzmpiDBufibxjD6LNPaPOgjOJ8JYMm09lELTbO72psbVxFkh2HmQX3O825ozVPURKsuAEcA7EzWKmXmA8xMys2aGZ3xNeQKFcUFpo5dfhnFaK5B9uaFJh3qq0t8mWk9NDETHU9YmBfC6gWlLq/IWdb7tIWmHwZ9eKsqOkAwcwtcLiPywhr8xtb6KnTpnXzOdg9hxwovDOSam77RKfSPaUJVWpiH+uoSl1fkHESEtYvK8b4+gaujZwQLK7zfo1vNCG9YWIb8vGDZPX6d+uTlHtYGCf9S9Bpi4zZLRKstP79sOcbubTaTv6wgzMfWeiXT2gejGK3WcSgUjAxrAz/GkU3xY49aYnZQ8xQ6e0cxo88W9jpenvJkMN/WjZSvc/2nk8M3QXCNrUpi17lbPhBk5SITJHe1gbn0yR/WWIeH5+06QWVxAZbpte4zYY6OmPQyfmnUkowvRURSyBm2LIsJcvvdEUzNejuxS7WQg3jx92OrRnVWcBA/E8DsuVCHaHgVvzRqmS+GbE36sQ6v3w9zcwlB8DVVJQVYtaAU1/rGMRNmXLozYpoE5TWsYxeDhuoeveyTBCLz2MXgfSaA9nv9pL0XgLZx/WUsdXlFifGDuxqYR5CZOWEtMTOfcHY5guA+W+urcK1PayBwrnvY24Lc4+3OQ3ZZVlWC0sI8jE+H0T82jXujU6gr925iV9/oVDQbvLgghBU1wcqwNvDb1CdTD2sP/594024XBBfxS4MQTaBiF//lNcHJsDYIhchcj+zxi3+7pSY8aEl2Bn4TZDXRbqOHE+1sCbI+TrHSOmbR7XMJgh22LvPHKMbLAc+wNlDnO3vdbX3pTmwoSRAzrA3UUMLVvnHP51qocX0vfy52LeRBaAMPXnVgLca5TjpwLkFIG7X06dKdEUzPerOsw+uN8p3CTz2tL6nNJwKa0AUApYX5WFGreWTCEUZXr3czrQfGptGjt/oszA9h1QLvtvp0ymXtxNZ8UD9PowPnEoS0qS4tjF5spsMRz4qAai0GqYe1FfV383otcq5YyACwfpE/3NbWlpl5HvYkeSKGrE+H8m7mjJBz+MFtHfSELgP1d/PyGEZm9kV7RqdYv8QfguynRi1Jtc4koj8G8FCCQ5qJ6JU011ANoFl5PJjmeQTBMbbWV+Ffzt0B4N0GIe0BL3kyqK82Z1r3jU5hgQczrbsHJzA6pTUbrC4twCIftPm0gzmU4N2N0iUfzaZOtpf1J6CN9YsHQRPVPTbWQQCMUS7Se1lwHXOm9XCCI91hcHw6OgKvKD+EFbXBLK8BtEzrtYvKcUZvZdp+dxSPeVCQrQ1BEow8DQR+ybQ21SB73EJOxWWdqM2lE+0yja8H0vg9BMFR1BaaF24Pe65frxo/blzo7biYE6wzZVp78+Lvl9Iap2hcWA7jz+56/zgmpr2Xac1sntrm9c8lWQv5OMwj/Az2QbNsBwEctbmWTmjzea/aPI8g2Ka2rBD11SXoHpzA9GwEl3tGsWmpdyryTO7qACd0GazzgXs0F1pmqhQX5GH1gjJ03RsDszZoQt3IeoE7w5MYmdTCCBXF+VhSWezyihKTlCAz81fjfZ+I9ul3TzLzc46tShA8wNb6SnQPag3pz3YPeUqQze0Zg3/x90Pp0yUfWWJOsW5xeXRc6aU7I54TZOtn4vUwghNZ1t7+DQUhTbycaa2KUuPCHLCQPd4cZCYcQWev6rXIDUE2xZE9GErwm9ciWZf1XBjx3njubEHwNVuXe1eQ1QvNpqXev9DYpb66BCUFeZiYCaPPg5nWV++NYSas5aUuqypGVUmByyvKDiZBvuNBQfZRyRNgU5DncmULQhBQLeTzt4cxG44gP8/90v3ekdgAg9LCvMAOMFAJhQjrFscyrTt6Rj0lyGrikNczeZ3E6+MxVU+S10ueAI80BhEEL7KwoiiaBDI5E0GnR9oDqtbxusXBHWBgZe0itWOXt6yxdp9ZYk6xpq4M+frfn1qH7QXCETZ1dvODy1oEWRASoCapeMVtfVFpz7gphy7+6gXVax27Lvqo+YSTFOaHsKYu1hvaSxula31jmNL70C+qKEJNWaHLK5qfhIJMRGHlNktEqy0/v2w5xu7NO9srQYA3RzFe9MnkGqfx8hjGXLWQAbOL3ku9xv0y4UllPgtZbdoxl1/MblOQuZqNCILrqJOfvGIh+/FC4wQmC9lDF/7x6Vlc7x8HAOSFKCey3lXUIROXPLRRuuSz+DGQnMtaRFLIWVQL+fztYYQjnODozBOOsMka27jEO7XRmcboaQ0AfWOx1qFu03F3FKz/WaxeUIrigjx3F5RlvFojfsmHiXbzZVnvtTzutzzeD5nSJASYRZXFWFRRhJ6RKYxPh3Hl3ijWujhZ6aolLlbrg7iYU4RChPWLK3D6hjZ/5sLtYSysWOjyqqzNJ3Jng2Tg1alPF27Hci380qgloSAz82vz/PyEs8sRBO+xtb4KP7rYA0CLI7spyKaLv4c6h2WLzcsqTYL80fXuC/IFJcnOD5m8TrOqthSFeSFMhyO4OzyFoYkZ1+uwR6dmcU0JI/jlc5Esa0GYB3OmtbuTn/zUKD8TqO1LVQvITc7fiq1j87Lc2yTl54XQ6LGStEt3hqNhhLULy30TRhBBFoR58FKm9UVFhPySqOIkm5WuZBduu3/hZ2acv53bggyY48heSOxSN0l+6mTniiAT0Woi2mEtoxIEL6JmWp+/NYyIi4ldphF/PrrQOMUGJUZ7uXcUkzPujvy7OTARnSZUXVqAZVXeniaUKbyWAe/XTZKjgqyL7I4EP3+JiMLQRi22Augkog4i+oKT6xAEJ1lSWYy6ci15anRqFlf73OnYNTZlLq9RO1flCuVF+Vi9QGsVGo6w64MmTBf+pZWenyaUKVRBVhvXuIUpjLDUWxOoEmFbkImokohe14W2FcAzcxz3HoDncX/dcSOAI0T0TbtrEYRMQESmOLJbbuv2uyPRuNiaujIU5fsjLuY0ahz5vMtx5A9umQU5V1HzGc7fGgaze16k2XDElGuRMy5rInoQwBUAu5CgXpmIXgKwUz/G+kmx/v29YikLXsULoxjP+7CMIxN4KbEr1xO6DJbXlEQzq4cnZ3FzYMK1tailgUsqiz01hGQ+7FrILwOoUR7fJ8pEVAVtTCMrxxwAsBvACwCGEBPlQzbXIwgZwQuZ1urrem0QfDYxWci33BXkCz6NVTqN5kWK/f4f3HIv+fEDH2+S0hZkIvo8gCbExLQNmsh+xXLo08ZT9GP3MvNXmfkEMx8C0Kz8rFqsZMGLbFNnI98acsUld165yKkWe66xyZRp7Z57dHB8Gt2DmiVYmBfKuZaZVrYsc3/TCpg9SX5yVwP2LGQ1VtzKzM26yFq3Rmq3ry5rsxFm7oJmaRvW9SdsrEkQMsKyqmLUlGouuZHJWVzrG8/q68+EI7igxMW2+Gzn7yT11SWoLNZ6Gg1PzuLW0KQr61Av/OuXlKPAA7Oy3UT9mzznooXs14QuwJ4gNyj3X0pw3C5o1i8DODbHMa/OcV7BAYiomoiaiEje2zSxJnZl+4LT2TuKaT0uVl9d4otRcpmCiMxxZJfc1uclocuE+v/xgUufCTP7Oq7vlCC3xTtAT/oCYtbv8TnO1aUc52vR0IXvMBF1EhET0YB+/zARNbm0rKPQMuCPuvT6gcDNBiGqC9BvF5lM4IVMa2vJU66zZkEZyvThH70jU+gZzr7nondkCn1j0wCA0sI8rKotzfoa7GBHkNWhEtahEwa71AfM/KMUz+sbdCv0ODTh24fYxqJav78PQCsRHSWirP2O+mvtmvdAYV7MiV3ZFeQPJH5sQt2UuJb1brLE5DMJhcyeCzfc1h+Y4seVCIX8VRduR5BVq3guq9aIMzPmsKLjPH/QxppcQRe9ViQnfHsAZHMox4tZfK1As82SaZ3NZCLVBZjL8WODbS5ujgBgajZsakrit+ShTGFyW7uQ2HX2ZuxvwY//J3YEuUu5/7T1h3q5k5GFDQCvJDjXbuX+XNa2lzkB86aiC1oyW6N+2w/z+9VERBl3HxPRLmjNWAQHUGsthyZmslZrGYmY42K5XPJksG5ROYrytcvXraFJ3BvN7mzkC7dHMKu3UF21oBQVxe5ON/IKbid2nbkZs+e2L/efs9WOIJ/UvxKA/US0yvLzl5WfA3MkdOnC/TxiiV+JLGnPQUR7oG08DFqYuZGZjzFzl347Aq0xivq77clUTNmIY2PumL2QBtZay2zFka/3j2N0SuuXvKCsEIsr/dPoIFPk54VMF/9sx/T9fuHPFG6XPp1RLOTty/23cU1bkPUa4kFoIloDoE3vVf1FInodmms2KrLMfNV6DiL6OMzCDgCH012TSxxU7g/CXOYVhZnj/exgvGPTQY9NtxIRIxbHFhxGvfieuj6QlddULY0t9VU52y/ZivpZqK7KbPD+jdjrPeDDC3+mWLe4HIV6+Vf34AQGx6ez9tp3hibRM6J5SkoL83xZF263cO4riDX1qIFm6R5GrNTJuHIcUJ9ERA/qva+PQ3PpGsLdkmTilyfQy4hM5V+68MZFr7k+onxrl4MJXlZLXcgAO1fGGtOdvJYlQe6W+HE8VNf9mSwLsljI8SnIC5mmkGXTSlY/k63LqpDns4QuwKYg61bya4jfo9p4Nw7FEdla5eeGcA/BItw+YL/l8ZG4R5mxegDui7+nSRc0C916ExykaVVMkM91D2Vl/N/pGzHh3y7x4yiqS/Jsd/b+1EenZnG5V0voCpF5PKdg3ii9fzN7n4vf3dWAA9OemHkvgOcAnIJ5ilMXgN3MnCjL1zj2GIA1zHza7nqyjGqRdiWyjg2Y2Roj3+nEQvS4dY31BqDFifMLGrVlhWioKwMAzIQ54xm+4QibLjQ7Voo1ZtC4sBwlBVrd693hKdzNUt3rue6h6NStdYsqUFqYn5XX9Qs7FI/B6RtZFGTlf3FbrgoyADDzEb11ZghADTOHmHktM89V3tMFbZDEXgCNzPx0nJabfqBZuZ9KMpp6bPOcRwmeRLWSWzPstm6/O4Lxac0KX1xZhKVVJRl9PT+RF7Ik2WXJbW12V/vzwp9JHlyp5lkMZqU8kJlNn8sDPg0jON58NRlhZeYrzPwCM7/GzFecXkMWUT/1rjmPuh/1WF93JstFmrMoyKqFsWOFPy8ymWRbvZLYlaVM6/dV16h8JvfRuLAcFUWa1+De6FR0AEcmuTkwgcHxGQBAZXE+Vi3wV4cug9zuhm6DOH2hO1N4uirIvvmPJqJ9RHSSiE729va6vRzX2KkIctv1gYxaAKevq4Jck+DI3MQcR86+hSwZ1vcTChG2r4i9L9lwW79vSbLzayVCRgWZiFYT0Q4i+rj+dTURBSUDwiqkqTQ06VMfZLOVph2U0ETzwoUL3V6OazQuLI9OG7o3Op3RyU9iISdGjRWevpF592j/2DRu9MdGLm5cEpTLmbM8qGweT13PvCCreRZ+jR8DDguyLrpfJ6IOvaypE1pNrNHjuRPAgP7zrxPRDidf32UkozlHCIUoK3Hk0alZtPdoIxdDJPHKeKxZUBbtntY/ltnNEWC2jjctq0RhvjgZ46FuHrNhIav/g37euDry16Rbvq8j1pCiEeaMa+utEbFhC+/G6fLlB2ptPNf6F2rnXIILZKMe+czNwWg27/rFFSgrkmxeK6EQoUlJIsp0TL9NvfDLBmlO1GqAs91D0dGhmWByJmxK6FNDSn7DtiAT0VPQLN9diNUWJ/10aFnGXUT0Mbtr8RH+3cIJAIDm1bE91DtdfQmOTJ9T18VdnQzqBbg1w93T1M3XztWyj56LuvIirKjVKgKmZyO4eCdzDULOdQ9hOqwJ/pq6MtSV+7e1rK0ttz7v2OiXbARvCFpZTws0oe6CFl+thZZR3AhNvJuU5xGAFiJqYub37awpi1hjxnaumH4cqJHTPLiyGoX5IUzPRtB1bwx3hyexuLLY0dd490rsz6LJx7v+TKO+N20ZtJBnwhGT+/Wh1fKZJGLHippovP30jcGMdTRTvSJNK/39mdj1gRkTi6KiCuAAM5+a4/hoXbI+WOEggKeU5x8FsD7VRejnegaZsTwHARzW215av6+SynZ5gfogmYYigrcoLshD08pqvN2lieYvOvvwqw/WO3b+cIRNF5pH1yxIcHRu88DyauSFCOEI49LdEQxPzqAyA9OXLtwejtaE11eXSE34PDy4ohrfff8WAE00f/ex1Rl5HdVr0ezzTVLagkxEn4dm8RqW8cF5unKZ0DtW7SaigwD+RP92IxH9GjP/Y4rLyfSYwU5Y2mIyc5cltb4xhfNJ7XEAeKyhLmOCfOH2cHTC05LK4qj7T7ifsqJ8bFxSgQ9uDYNZKxX76HrnqwBOXlXc1eKxmBdVHN/p6gczO16OxMwmr0izzz8XOzHkZ5T7bamIsQozH4C5c9Wvp3GOQ8xMGbzN1aNatWxTEVn1WF+NmxRiPNYYs1p/4XAc+R3FXf3Qmlrf1lVmC2tteCZoDZAllg02L61EuZ6IeGd4Etf7nc+Av9o3jr4xbaJUVUmBLyc8qdgRZLWP80s21/GVOc7rdU4q91NZt3rsyTmPEjzNAyuqUFyg/Qtd7x93tCPRu1diAv/wGkkemo+dGS5DY2a8dzW2SWpeJZ/JfOTnhe6zkp3mpPKZNK2sRsiHE55U7Aiy+hdp18pr1b8S/OXOPa7cb0imwUecDl/H4x4oeJ6i/DzThfkXnc5YydrFPyYqj4ggz4sqyCevDjheZnNzYCI6a7eiKB8bllTM8wwBAB5Rch/evuJ8NYLqSWoOQNa7HUFWxcfu1sevWcbWSUrJjFK0jmyUaUw+RnVbv3n5niPnvHR3BP26G666tABrfe6GywbLa0qxslbrXzwxE3Z87N/bSkhix8pqX87adYNHGtTyQGcv88xs2gSr/4t+xY4gq1nHdicWqVajbzKO9cQ09X1IOM9Zt6D3Kd9qmSvDmoiqrTf7Kxac5iPr6qL3f9rei0jEfuvGn7bH+oR/aG2d791w2eJx5YL81mVnrbG3lAv/4411CY4UVLbVV6G0UBuR2T04gZsDzsWR1TBReVF+IGaFOyXI++Y8KjkMq5Hhv5iqKsINetb4XLwMs2chkYC/CGBAuSU6r+ASW5dVoa68EADQNzbtyICDn3XELO2PrpOLf7KoFtJbnc54KwDNElO9Hx9a639LLFsU5IVM4QS1tt4u6ibp4TW1yM/zfxtTO7/BMf0rAdibbqctvdPXPsTKp44mONxzMPMxmGPozxPRYTVWTERNRHQcwB7luGO6hS34mFCI8MT6RdHHb1zqsXW+yZmwKS72kXW5O8QjVVRBPnV9EJMzYUfO29k7Go0fVxbnY8sy/1ti2UTNgXjTQc+F2WsRjE2SHUF+BZp7We209YVUTkBEzwL4oXKOQQCv2liTWzwFs6t9H4BOImIiYmhJa7uUn7cx895sLlDIHE9uiInmG5fsjaV850p/NCFp7aJyLKuW+uNkWVRRjHWLtHj7dDjiWLa1KiKPNS6Q+HGKfGitEtbp6HVkIpcWP455LYIQPwZsb8yekgAAFydJREFUCDIzDwF4AZqQGoJ6hIj6iOjL+uQn02wyIqrUv/8SEfUB+AZiAycYWpevzDU9zRB6HHgNkkvQOsbMOzO8JCGLfHTdQhjX6DM3B9E3OpX2uX6mxI8/Iu7qlMmE29rsrpbPJFW2L69GdanWOa13ZAoXbo/YPuf528O4NxpLfNwUkDGYtpzuesOMr8IsyjXQYqOt0EYthnWRDkOLhbZC66pVozwPAA4x81/YWY+bMPMgM+8GsBuaO1+NsXdB6/S1M1nLmJkPWJqTWLOzk13XbuUcshHIAFWlBdE4GTPww/N30zoPM6PlQuy5HxV3dcqorsufttsX5OnZiCmTVxK6UicvRKbQy0/a7XmRAOCNi7HQ0EfWLQxM4qPtKLjeaes5AEY2izpkwrjVWB5bj9ufbqcvr8HMLcy8l5kbFSFsZOb9EjMOLp/cujR6//tnb6d1jva7o7iqz/MtK8wLjBsumzzWWId8/eJ8tnsId4cnbZ3vvav9GNFbmC6vKUHjwjLba8xFnlivCrK9PAsA+JEiyE9tXJTgSH/hSFqabimvgebCvoL5xzASNAE/AKCGmV92Yh2C4Baf3rYkev+tzr5oHXEq/ODcnej9JzcuQnFBniNryyWqSgrwkNIgQr1wp4PqsXhq4yJpYZomarXAyasDGJmcSftcfaNTOKVP3QqRWez9jmN54sw8pPeUXgvNIt4LTXAPQXPXHtIf74UmwrXM/FU9Fi0IvmZpVUnUbR2OMF7/4M48z7gf9Tm/tGVJgiOFRDy1KWYxnbiQXvgA0EIIJy4oltimxbbWlcssqizG5qVanHc2wvixjeTHn7T3wsgLa1pZg5qyQieW6AkyUrili/NruuC+wMzP6V+/qn9fRFgIHJ/eFnNb//OZWyk9t+PuCM7f1vIZC/NC+NiG4Oz6s80uRTh/fvle2uVPl3tGowMRygrzTF2nhNT55NbYJvNfzqUX1gFg2iR9fFNw3NVAhgRZEHKRT29bAsOj+eblPtxIYbrNsbab0ftPbVqEigzM880VVteVRWO9kzMR/LwjveQuNTnvo+sXoihfQgh2UMM6b1zsxfj0bMrnGJ+eNYUhdgXMa+G4IOulTZ/TS5teIaLXieg9/fY6EX1d/3kw8tQFQWdpVYkpnvXKezeSel44wvj2qe7o4883LXd8bbnG7s2xi/933k/NW2HwXeV5uzcH68LvBmsXVUTrxCdmwvhJGm7rExd6MKF7PNYvLsf6xcEa8uGYIOv1xa9DK206Cq20aQ+0hhhN+m0XtKYZR6GVRP0g3Q5fguBFfv2hldH7r568gZnw/FOHfnypB3eHtdrlBWWFeELc1bb57APLovd/eP4ORqdSs8Y67o7g4h2tXrYoPySC7BCfUsM6aVQjqJukz2xfluBIf+KIIBPRK4h1o1JLm4D7y53U7+2G1uHrm06sQxDc5qlNi7CwoggA0DMyhX86Pb919vLPYiXrn9+5HAUB6MnrNpuWVmD9Ys0am5yJ4IcpJtmpVrWEEJzjlxVBPn7+LgbHk69GGBibNiWDfWb70gRH+xPb//lE9B40SzheffEQtD7PLfrXoTmO20tE79pdiyC4TUFeCL/3+Oro46/9+DLCCSZAnb05hLf1sXT5ITI9V0gfIsKv7KiPPv52Ehsjg3CE8a22WAhBtbYFe2xYUoFt+lSm6dmIKVQzH6+13cS07nF6YHkVGgI4ltSWIBPR1wEY3Z+MTl1XAOxl5pBe2tTMzJ/Qv9Yycwja3OBTMHf42klEX7OzHkHwAr/96CpUFOUDALp6x/BPp+e+6Py3lvbo/V/evlR6VzvIr+yICenPO3qTTrJ742JPdKxfbVkhntwQrExet3nmoRXR+99870ZSva2ZGX//7vXo4994eGWCo/1L2oJMRA9CG5toCCpBE+K1zPxaoucy8zFmbgbwCZhFeT8R7Uh3TYLgBapKCvC7j6+KPv7y9y9iaOL+Rgg/ae+NZowSAc890Zi1NeYCy2tKo/3AIwz8zS+uJvW8v3vnWvT+3ubl0qDFYT67YxmKCzTpuXhnxNSadC7evNyHrt4xANrs438VUK+FHQv5Gf2rIag75xNiK8zcAqBZOQcQm40sCL7lS0+uxeJKLZZ8b3QK/+E7H5gsgYGxaRw4dib6eE/TcmxaKoUHTvNvPrQ6ev+b793A2DzJXZfujER7LRMBv/XwqoTHC6lTWVxgqiT42o87533O//xRR/T+55rqUaZ7oIKGHUE2ZvsygCPMfCqdk+j9nV9GzMrelfgZguB9yovy8e8/syX6+B9PdeM/f+8CZsIR9IxM4vf++j3c0fss15YV4k8+ucGtpQaaJ9cvwpo6rSZ5ZHIWf/nzKwmP/x8n2qNdoJ7auBgrF5Rmeok5yb6PNkQnpP388r2EozLf6eqLzgjPDxGe/UjDnMf6HTuCrLatOWpzHeoM5OC+20JO8eltS/BMcyxe9hc/v4JHv3wCT371x3j/Rmx89sHPb8eiimI3lhh4QiHCc0/ELimHf9o153jMc91D+P7ZWDb2H+5al/H15SqrFpSZ3M7/4TsfxE1+DEcY/+l7F6KPP9dUjxW1wd0k2RHkauV+15xHJYfd5wuC5yAi/Mdf3WLqJtQ3No3x6Vgrx//0q1ulxjXDfL5pOdbqDSlGp2bxn79/4b5jZsMR/Ok/no0+3r15Mbbq2cBCZvij3RtQmK9J0NnuIfz1W1fvO+Yvf34FZ7u1TsuF+SH8248Fe5NkR5BVEa2e86jUYIg4CwGiKD8Ph39nJ1741EbUKk3w1y0qx//+4iP47UclRplp8vNCeOGTG6OPv9XWjWOtN03HfOVfLuLMzdiF/8VPbYSQWVYuKMXvPxlLZHzp+xdMs5LfuNiDr/zgYvTx7z/ZGPgQgp3IeAu0rlsA8BSA0zbOtUe5LzODhUCRFyI890QjvvDhNbjWN4bigjzUV5fIKL8ssmvzYvzqjmXReuQDr53BvdEp7N68GH/15hX83duxkpo/eGpdIGtcvchzTzTixIUenO0ewmyE8YW/fg+/8fBKMBjffPdG1I29Y0U1vvRk8KsQKJkasLhP1MqeWqFZtYMAVjPzSBrnqYJWu1ytn2s3M/8orUUJWaO5uZlPnjzp9jIEIWmGJ2fw9Dd+EW2JGY/dmxfj8G/vRCgkm6Vs0T04gT1ffwu3hybj/ry+ugTHvvQYllYFo0afiFr1st/7SNtlrWdVG9nR1QBOEFFKnb71ARMtiIlxi4ixIAiZoLK4AH/zfzwc7RRl5ZNbluDPf/1BEeMsU19dgm/9/uNoWnl/5PPh1bV47UuPB0aM5yNtCzl6AqJWAA/qDxnAYQCHmPlqguesBnAAMZc3AeiEVss8bGtBQlYQC1nwK9OzEfzt29fwndPduDM8idULyvDbj67CZ7YvlTCCizAzftLeGy1xeqxhAT6yri5wn0kiC9mOy7oSWokSQbOUm/QfGSfsUm6D0KzgBuUGmPtaH4EmyvMxyMx/kdaiBccQQRYEQUidRIJsJ6nrGQDfUB6rwyIAoBHxa4rV7Y66G9hnPXAO2gCIIAuCIAiBwm7/sXjiajW5rcfMZZIHyy8hCIIgCClgR5D74U7NsNQpC4IgCIEjbUHWB0mkNExCEARBEIT42JqHLAiCIAiCM4ggC4IgCIIHEEEWBEEQBA9guzGIkJsQUS+Aa2k+vQ7APQeXIwjpIH+HghusYuaF8X4ggixkHSI6OVdhvCBkC/k7FLyGuKwFQRAEwQOIIAuCIAiCBxBBFtzgiNsLEATI36HgMSSGLAiCIAgeQCxkQRAEQfAAIsiCIAiC4AFEkIWMQkRNRHSYiDqJiIloQL9/mIia5j+DIKQHER3X/+ZSvQ24vXYhN5EYspARiKgawFEAu+Y59BiAZ5l5MPOrEnIJIupE/Jns8zHIzDVOr0cQ5kMEWXAcXYxbkfzFsI2Zd2ZwSUIOQkTpXtxEkAVXsDMPWRDm4gTMYtwF4ACANv3xLv2xcUwTER1l5r3ZW6IQZPRNocohAH1urEUQkkUsZMFRiGgPNFe1QQsz745zXDU04VbjyDuZuc16rCCkip6f0Go8ZmZycTmCkBSS1CU4zUHl/iCAuFavHjO2/uxgvGMFIQ1UD43kJwi+QARZcAwiaoD5QvhSomQtZu6CuVvSrjiuRkFIh4eU+12urUIQUkAEWXCS/ZbHybQmPGx5/LRDaxFyG2sOgyB4HhFkwUnUeHBXMqVMcWLGkm0tOIEIsuA7RJAFJ1Fny6aSnKUeK/NpBSdQBfk96w8lNCJ4ERFkwUnUi1wqVol6bDqNHATBiulvkYga1I5xAAb0rlydRHRQz38QBFcRQRYcIc4FrTOFp6uCLJaLYIs4f4svQvt73If7N3wNAJ4H0ElE1nwGQcgqIsiCU1iFtD+F55oaNog7UbCJtUf6niSft4+IWuc/TBAygwiykCmk9lNwi3ju5zZode81epOQRv3xMctxTUR0PMPrE4S4iCALTlFr47lW8bZzLkFotDw+xMw7mfmYkfnPzF364724v0HNLiKabyiKIDiO9LIWvIC4qAUnaYXWuxoAOpk5YT08Mx8jov0w18QfhJTgCVlGBFlwCmvM2I7IphJ/FgQT8wnwXM/RRdmIPzcRUbWMBRWyibisBaew43ZeoD6Qi6DgEi2Wx1ITL2QVEWTBEfS+1CrWOF4ipAZU8ALWBiISShGyigiy4CSqZZuKyKrHyvhFwS2sm0pJLhSyigiy4CQnlfvWWtBEqMeenPMoQcgs1k2k9MAWsookdQlOchyAUS7SkExSTJyuSlIDKqQNEe0B8IzxWC9rShYRZMFVxEIWnMSaFJPMKEXryEbrOQQhVfYYN12gk+UZ5X5XnLwIQcgoIsiCY+ijFNWL2IFEx+stMvcp32qRDGvBDsxs7bx1MJnnEdE+mEMn1vMIQsYRQRacRhXhBiJKdEF8GeZM1oQCLghJotYhN8zXClO3otW/00EAL2ViYYKQCGJmt9cgBAy9Qb9qbRwBcNBwARJRE7QLoNqe8FiK8T5BiIvuebkC82bPENkW/X41tDrjvTD/HQLAbmaW0ImQdUSQBceZ44KYiDZmljaFgmPoyYKtSL2WeH86nb4EwQnEZS04jh4HXoPkErSOiRgLTqN7Y9Yg+VhwF4CdIsaCm4iFLGQUfWqO0SPYKCvpgibWh/VEMEHIGLq1bJRDNUCzmgeh9UxvAXBUXNSCFxBBFgRBEAQPIC5rQRAEQfAAIsiCIAiC4AFEkAVBEATBA4ggC4IgCIIHEEEWBEEQBA8ggiwIgiAIHkAEWRAEQRA8gAiyIAiCIHgAEWRBEARB8AAiyIIgCILgAUSQBSGHIaI9RMTGze31+BUiGtDfQ+soRydfo0l/jQF9opoQMESQBUEQbEBEh6ENrGjL5JAKfRBLi/5aL2fqdQT3EEEWhABCRM8T0XH99rzb6wkqRNQEYJ/+8EAWXtJ4jT2ZtMYFdxBBFoRg0ghgl357yOW1BJmj+teWbIxw1K1kY8bz4Uy/npBdRJAFIbdpgSbexk1IEiLah9iM72xYxwbGazWI9yNYiCALQg7DzIPM3GXc3F6PzzCEsUu3XLOC/jkZr/ditl5XyDwiyIIgCCmix28N69gN17HxmtVEtMeF1xcygAiyIAhpQUTVevJYq1L200lER/VkJ+vxDUR0UD9m3uOTeP1dRHTYcr5W/XyZFinVRX1szqMs6O/XUX2dRglTOmt+VbkvVnJQYGa5yU1uAbgBeB4AJ3GrVp6zC0CncZvjvE3KczuV5w3M8zoHlXPsS2Jdh5P8PRsAHE/ifJ0A9mTgfW5QXqM1hc9mvveLAbQCaErynOp7kNRz5Obtm1jIgpDbVEMTGOM2L7q79rj+3EQ8rzceeR7JuXX36YlSiV67AZpoJVPy0wDg6HznTAPVkn1lvoOJ6CiAg5j//QK0zc+JJBt/HFXu70/ieMHjiCALQnDogpY13QJgUPn+oPJ9u6U5tdDEGNASi3YzMzEzAdiJWLKRgSFGxvrU4xvjrOcg5kAXqVaYha0Nmhjt1M+3O845Dzvswn5GuZ/w/dQ3A+prDwI4BG29NdDWvBfm960aZrGdC/W1pSY5AJDu+hAEIUDo3aMMy/AYM++d47g9UC7+ulBaj2mCJoQqLcy8e45zduJ+azvR8a3QLEODGmYejHPcUZjF7RAzxy030q1yVdy7mNmRsi61xWi898tyrPq7DQLYyXNks1vfh/nOrT9nALENStz3TfAPYiELgpAOcQVeJ56Vm8ilahXVZusBuqtaFeNjc4kxADDzIWiWqIEjNbuW7ljJlDqpG40jc4mxjun3STLR7aRyX6xknyOCLAhCqhyaxxI7aXk8nxBZj48XP7WK77MJzgcA0AVbXWdcCz1F1HNY120iThy4c55zn4S2ibBuJhKhbgqc+P0EF8l3ewGCIPiO9+b5uVWsre5uE8w8SDSvd1a1/o6l4Jo9Ai3D2XqOdFGt1oQCG+f32quvZ87jkXrHL/WzuM+zIPgLsZAFQUiV+Vy1/ZbHCS3J+dAtTTUmfXyuY+Ng2jzorm871Cr3k+lsph6zS683dnJ0oroxsfu7CS4jgiwIQqpYBTfTWC2/VATeunmwK1rq85Ox0q2x8z0AjGYgz6fTEMWC+lnIjGSfI4IsCILXsQqN0eVq3hvudyvbFS31+fNuTFibABUvHtwELfnN+F2OE9G+NKxn06bAYetbyDIiyIIgeJ3a+Q/J/LniiF1ScWw9uWw3Erv6d0FrnjJARHPWYsfBuilw8r0SsowkdQmC4DfsNDdJe6JVkslncz23BcBOpXzrGZgTxFSeJ6JdzLwzjZfKdjhBcBARZEEQvI5VZPZ6pAFGA1IUeL38K1rWpNc174Ym0mp8uomIDiaqtdYxWcQeeV+ENBGXtSAIXscqem66ZdW12I7XMnMLMx/Qu4hZm60k08gk1SQzwcOIIAuC4GmY2Rp7TbqeWB8RuUu/2c1oBsyil3BjoI+bTPq1mfkY7u/WNV9WeEpJZoK3EUEWBMEPqKKcSvOMo9Dqlo8jweCKFEjFQt6vvHbC5igK1s3HfK+Ral204GFEkAVB8AMvKfcbkhmpqFulqjWdzAjI+VCbkszXqtLalCQZy95qSc/nhlYTv1JpmCJ4EBFkQQg+vq9N1d25qgWYcKSi7uo9EeccdlEzvOdrVWm1dg8nqhPW1/yi8q2ueXqAA+YNh93RmoLLiCALQvBpDkjDCGvS01G9ocYePV5rxIsPQ2sIUp3guWmhC6RhtVYnivHqx6qbgAbodcb6OhuUOHO8NSfjYo++fpxYu+AzpOxJEIKJ2qGqGpoQDOr3fTk3l5nbiGg/zK7nXZg/yeuAQ9axwauIzZpuQuLY7bPQ1qcK7fOYP4P6GDPPOYgCuM8FLtZxABALWRCCSTwB8r2VrIvUbiQ/2GGvPhvZSZKOI+sbn6eQvGAOAtjPzMlY9OprH03y/IKHIWZ2ew2CIGQA3Z16EGYLrQvATj9ayFb0GLLR8aoBmph16bfj81mYNl97ANp7OsjMNUk+Zxe0zOsmaNnR1bCxZiLqhO6yZub0WogJnkIEWRAEIUWI6HnEYrx7HXaJJ/P6TYiVUh1hZutUKcGHiMtaEAQhdVRL9hkXXl8VYCfqqwUPIIIsCIKQIrrL3xDlPS5ksT+tf21JojRK8AkiyIIgCOmhWqbzNipxCr0pirEBSKVrmeBxJIYsCIKQJvrsYqOEKSvlZEpCmcSOA4ZYyIIgCGmij0c0XMYvJjrWCfRkMiO7W8Q4YIiFLAiCYANLxnNGrWTFOt7NzNIMJGCIIAuCINhEr4lugJZklZEWlnpd+R5o1nHGaqwF9xBBFgRBEAQPIDFkQRAEQfAAIsiCIAiC4AFEkAVBEATBA4ggC4IgCIIHEEEWBEEQBA8ggiwIgiAIHuD/BwJ9RGT0MSD7AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(t,r[:,0])\n",
"plt.xlabel('time (s)')\n",
"plt.ylabel('position (m)')"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Runge-Kutta\n",
"\n",
"\n",
"Using `solve_ivp` in Python you can use more advanced integration routines, such as [Runge-Kutta 5(4)](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#r179348322575-1), [Runge-Kutta 3(2)](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#r179348322575-3), [Adams/BDF method](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html#r179348322575-7), etc. \n",
"\n",
"Numerical integration algorithms for differential equations:"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [],
"source": [
"from scipy.integrate import solve_ivp # import the ordinary differential equation integrator in Python"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'position (m)')"
]
},
"execution_count": 70,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAusAAAE9CAYAAACydiySAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOy9e3hT15X3/926Wja25Ttgm4uNIZgQYiApEEiahNQkk6YOE9KBTPu+0GnizO/XduY3OAkD00kb+OVi+s4k7RsI7RtPmhZaSBM3T6aB5tImOEASDEkAcTE2GNtgY9mWLXyRZGm/fxzp6Ej4oqt1jrQ+z3Me7320zzlLPrb0PWuvvRbjnIMgCIIgCIIgCPmhirUBBEEQBEEQBEGMDIl1giAIgiAIgpApJNYJgiAIgiAIQqaQWCcIgiAIgiAImUJinSAIgiAIgiBkCol1giAIgiAIgpApmlgbQCiT7OxsPmPGjFibQRAEoSjq6+vNnPOcWNtBEIRyILFOhMSMGTNw9OjRWJtBEAShKBhjzbG2gSAIZUFhMARBEARBEAQhU0isEwRBEARBEIRMIbFOEARBEARBEDKFxDpBEARBEARByBQS6wRBEARBEAQhU0isEwRBEARBEIRModSNBEEERK+tF/vO7cPRjqPQqrRYNnUZVpeshl6tj7VpBEEQBBG3kFgnCGJc6jvq8c9/+Wf02HrEfX9t+St+f/BpvNhhxvThYSAlF6hqiKGVCUh1CcxDZtQZDOhXMcy1OVBms4HRvSAIgogboi7WGWNpADIBdHPO+6J9PYIgIssvPvkQrzRsBJjjutcadTp8b0ouXr/cgSn9V2NgXeLi4i7s0NrwanY+7Com7l8wZMPWzi7MiJ1pBEEQRASJWMw6YyyNMbaaMbaDMfY5Y6yLMeYE0AOgEUAPY8zp3t/AGDvAGNvIGLs5UjYQBBFZdn9+BjvP/Jso1F3Dk/CTzi78S1cPklwuAECHRoMnc7MwHEtDEwwXd+Gpg09hZ0a6j1AHgC+T9Pju1Dyc7T4bI+sIgiCISBKWZ90ttB8DsBJAkf/LoxyWAcDoHr/SfR4AeB/AXs75/wnHJoIgIsf2o9vBkoUJMe40YOBiJVarfwgAmGO34/HJuXAyhuNJSfhtWirSj7ehoiw/liYnBC9/8TLevfCu2C+12VBsH8a7k5IxzBh61Gqs/sM/IL37CTzxjYV0TwiCIBRMSGKdMbYawCYAC6W7gznFCPtWAljJGNsF4BUAL3DOL4ZiH0EQ4XO66zQcyZ+K/aHLD4E7sgG10F86ZMPjll78IsMIAHjFmI6e2iMAlpA4jCI//+QD7GrYJX6KrumzYnNXD9QAvm21onJyLq6pVFDpLOg2/B6b3tQCAN2TBKC+vp4BmKNSqVap1epvAcjjnCfH2i6CILwwxgYAdDidzj+6XK79AM4uWrSIj3kM52O+7n+BjRBEutGzS/KyBcBRAMcAfO7ud0t+AkLsepH7+EwAizCyV95j1HsAnuScfxmwkcSEsHjxYn706NFYm0FEkR99+CN82PIhAMBhLcVQ63cBABeT1oljHABW50/BRZ0gCG3mO5DjWI1Pnrprwu1NBN481oIfH/0emP4KAOBrg0PY2X7Vx+vyQbIB/5SXI/YHmv8Bk3Xz6Z7IBMZYPed8caTPW19fn6vRaF7TarXFGRkZ6rS0tEG9Xm9Xq9Uu9+w1QRAxhnMOp9Opstlsur6+PkNPT4/L4XCcHx4e/h+LFi0adeFXQJ51xthdELzdRfAV6O8D2Afgfc75hQBO1QtAOu6XkmvcDeAeAA/BK97vAXAPY+wVzvk/BmIrQRDhc7b7rCjUAcDe+Q2x3cnTkcN6AQBaAP/UYxHFoS7jCC6fv3NCbU0knvv492AZglDnLi1+2Nl33Yf43QODuNtqxwepOgCAPvcALl8snmBLiYnELdT3TZkyZWpOTk4viXOCkCeMMWg0GpdGoxlKSUkZmjx5Mjo7O2dduXJlX319/ZrRBPuYYt2dyeWXEAS057//GAThvpdz3hupN8A5/wDABwCeYoyVAagEsAaCF/4xxtjDANZwzv8SqWsSBDEyvz39W7Fdmr4crYaZuGwbxFSjAZ+UH0ZFWT5ue+5DtFkGgSEXko3/AbW+E0xtQ9aU4wD+NnbGxylOlxP9ye+KWQHs3cvxQH+5+DoDoGIMTs7B2vuQkvICmGoYakMLcnMvxcRmIvrU19czjUbz2pQpU6bm5uZG7DuZIIjowxhDbm5uL+d8ant7+2v19fX3jRQSM142mAvwCvU3ABRzzhdzzn8ZSaHuD+f8OOf8Mc55JgTR3gshbOZ9xtg/ROu6BEEAVrsV+y/uF/ublj2GT566Cxee+xt88tRdYuxzVfkcGLRqACo4uleI45OzP4WLuyba7Ljnry1/hUovOF24Uw97l/d3nm804MJzf4OfPbwABq0afDgNDsst4us50w5OuL3EhDFHq9UW5+TkkFAnCIWSm5vbq9VqiwHMHun18cR6BgRv9yLO+cMBhrpEFM75LrdofwqCaDeOcwhBEGHwTtM7GBweBADMzpiNBTkLRhxXUZaPZ1fPR77RgOHeMsAlrGPrsbejvqN+wuxNFH539ndi296zVPx9G7RqVJXPAeB7TxxddwBcWA3c3H8S53rOTbzRRNRRqVSrMjIy1BT6QhDKhTEGo9GoVqlUq0Z6fTyxvohz/g3O+fEo2BYUnPMX3KL9D7G2hSDiluoS1H7872J3TcMRsJ8YgeqSEYdXlOULXvdnK/Dtud8U99eer426qYlEc18zjlw5AgBgUCHLeQcYBI/6s6vn+2R68dyTpq2PYNXMe8T9b5x7Y6LNJiYAtVr9rbS0tMFY20EQRHikp6cPqtXqipFeG1Osy0Gk+xML7z5BJAqXbN0w6fUAAC3nuK+/X3ghgOqkFbO8nzHvNb+HAcdAVGxMRP5wzuujuKPgdhx+4qHrwpJG4qHZD4ntdxq9MyZEXJGn1+vtsTaCIIjwcP8f5430WsQqmBIEoWxqj7fhQIo3JfNtA4NIcwWe2nVe1jwUpwtZRwaHB1HXVhdxGxMRF3fh3YveAkhSAT4et0y+BZm6qQAAq8OKO37xv1F7vC3iNhKxg3OerFaraZEIQSgctVrtGq0uAol1giAAANUHzmK/RKyX9wfnGWeMoVC/ROw/uf+3JAwjwJedX6K9vx0AkK5Px7KpywI+9u0vrqDzyo1iv099FJvePEH3Jc6geHWCUD5j/R+HVMF0lIukAViMMBeAcs7fjIxFBEEEw5X+NvRNEXJz61wcdw4EFzJRe7wN7x/NhWaa0HfoTmLTW8cAUPXMcNh/wZuZZ+W0ldCqtQEfW33gLAYH5mNS1p8BAJpJZ3HtSj+qD5yle0IQBKEQwhbr7lSKzyMyWVo4IvgAQRBE4GTlNGLI3f7a0BBSgqhuDLiFYX8uUuxZUOm6wNR22HVnUX0ghYRhiDhdTvy5+c9if9XMERMFjMplyyA4cuAcmgJ10hUw1TA0k0y4bFkYaVMJgiCIKBFWGAxjbAeEAkkZEHKxhzMXF+7xBEGEQWGBt3DOCn+vekruuMdftgwCYBi2zhP3aVLOuvcTofBF5xcwD5oBAFlJWbgl75ZxjvBlqtEAABju86bf1KSaxP0EQRCE/AlZrLurjD7m7nL3BnhFd7AbQRAxYnB4EJcGvhL7my3P4bakt1D7LRPwdC9Q1TDuOURheG2OuE8z6RymGJMib3CC8HHrx2L7zml3Qq1SB3W8p3DVsHWuuE8zqQH/fE9RxGwkCIIgoks4ISebJG0G4BiAJwEcjWZ1U4IgIs/n7Z/D5rQBAIrTi1H7zHeCPkdV+RxsevMEBgemg7t0YCo7VLpufO/rqZE2N2H47/MfettHjFiQ1BZUSJFn7AsHtOi1Z0Kl6wZT2VA4tQPAjAhbSxAEQUSDcMJgFsLrTX+Pc76Yc/4BCXWCUB5SD+7tBbeHdA5v9cxUOPuLxf1J6eN75Qk/qktw5ZlMdAxdBCDkvP9r37/jttqlQWdyqSjLx6Gn7sYjw95c+R/texh4Ol3YRil4RRAEQciDcMR6EbzhK4+NNZAgCHnzWftnYvu2/NtCPo+neua/r1wt7jt0+VBYtiUk/VdxMNkbV37L4BCSOUcO60X1gbMhnfL2vm6xfdAgiVkPoOAVQRAEETvCEesWT4NzfjF8UwiCiAWdA5240CsUBtaqtFiQs2CcI8ZHKvg/u/IZHE5H2OdMNA4avLH+KwaHxHaoC3ZvGRqCwSXUzrmo06JFE1z8OwE4nA7wILMkEQRBhEs4MevdiEy6RoIgYsjn7Z+L7QU5C5CkCX9BaGFqIQomFaD1WiuGnEM42XUSZbllYZ83URgG8LlUrEuy84SayUXPgcVDNtFj/1lSEgqv9YdlZ0JRXYLXNYP4TVoqFg3Z8JD1Gr42ZBMyJQWwAJsgokFtbW3qgw8+ODuUY1NTU52FhYW2BQsWDPzjP/5j5/LlywOqhFdYWHhja2urHgAKCgpsLS0tJ4O9dl1dXfKKFSvmSvdVVla279ixI+SKbbW1takAMHv2bFtpaak91PPIkXA863/wNBhjd0bAFiLKMMaMjLGFjDFKBUGISENgbp18a8TOu3jyYrF9tP1oxM6bCJzW6dCvEj6eJw8PY9rwsPhaVfmc0Q4bl69JPPSfGihLT1D0X8XRpCR0ajTYPykFlzUacT9BKBGr1ao2mUzJe/bsyV6xYsXcefPmzTWZTLpoXzcaQv3ee+8tevDBB2c/+OCDs3/+85/nBHLMsmXLShhji4Ld0tLSbg7VzlAJR6zvlLRfCdeQeMMtil9hjDUyxjhjrMfdfoUxFquKJPsA1Lt/EgQA4GiHV0jfMjm4PN5jsThPItY7SKwHw2cGvdi+ZdDmk9s2nAJTXxvyivXPkpJAAR2B4wRwPMl7XxYN2WJnDEFEAZPJlLxkyZLSaAr2kYT62rVrzeEI9e3bt2fv378/I9jjWlpa9OOPkgchh8Fwzi8wxiohiPZixth+AI9xzpsjZp0CYYwZIYjhlX4vGd3bowAeZYy9AeD7nHMLJgC3Xf42EQlOR38HmvuEf1m9Wo+bcm6K2LmlnvXjV4/D4XJAq9JG7PzxzOcpaWL7FonADqQ41aik5GJ2/1WkO53oVavRpVGjSatBsS4zDEsThzM6La65Zztyh4dRKJntIAi5UFlZ2V5eXt433rienh5NY2Oj7rXXXsvxhLQAgre9vLx8diihLeMxmlDfvXt3yLrRZDLpqqqqpodyrPR9y51wYtbBOd/lFoHPAbgHQBNjbBcE721IrjTO+Rfh2BRL3L+LegiZcsbjIfe4RVE1ysum8YcQicaxq8fE9oKcBdCpI+dQyd/xdUzN0OCyVoPB4UGYnp+KBTY7xfiOg8PlwLGUScCwEKd+a+XnwKTQvekiVQ1QAbj1r/8f3mt+DwDw6eqfo3juuvDPnQDUJ3nDhhYN2aiSHyFLiouLbRUVFdZAx2/durVjy5Ytedu2bSvw7GttbdXX1NRkrF+/vidSdtXV1SXfd999PrH14Qp1AFizZk3x+KOux2w2+6ywr6ysbM/KypLtE3hYYt1No6TNIHiOQ4UjMjbFig/gK9SbIBSK8iiile6+Z8xCxtg+zvmaaBrFGFsJ4IloXoNQJl92fim2I74AtP8qFhsy8bZ2EgDgaJJeEOsU4zsmpi4TBt1CfWrKVORHQqhLuGXyLaJY/6z9M6wjsR4QRyUhMIspBIaII7Zu3dpRX1+fIg0l2bt3b8TEukeoW61WUSBHQqg//vjj+SaTKTmUY8+cOePjVQ8nDGciCCdmHYyx5wDshbc4kjQEkgWxwa+tOBhjD0EoFOXhfc55Mef8Dc55k3vbBcGTfkwy7qFoxbB74uYBvBeN8xPK54ur3omsSKRs9Ecqao4n0YLGQJBm54nkGgIP0kXEx68ep1SEAeDiLhyT5KZfHKnQJIKQCZs3b26X9i9duhSREBGTyaSLhlCvra1N3blz5+RQj29oaBDfX2pqqjMcWyaCkMU6Y6wMgrdWKrD9BXjApwvVDhnxvKRtATCit9wdo+7/2vMjjQ0Fxtg+xlg9Y4xDCMkJZ6aDiGOGhodwtttbYCeS8eoeFti8Yv0rvY4WNAZAtMX6zPSZSNMJMfHdQ924ZL0U8WvEG+ct59GrEr6mMpMyMfNfzcDTvcJGIV1EHHDDDTf4TBdFYvGlyWTSLVmypDTSQt1sNqu/+93viuEv1dXVQZ/vs88+Ez3yhYWFsp8qC8ez7omB5hDE9gUIlUwXASgGkBHCpsiVTu5UiNLwl2fHWjjKOW8CsEuya6U73j0S+Hv4CWJETnWdwjAXQvSK0ouQrk+P+DVmOIaR6hQK8fSo1WjVKDnKLfq4uAtfdX4l9hflRX5Ji4qpfGZRpLMrxMhIU48uylsExuLBv0QQXrKzs328y+EK2GgJdQB44IEHijznXbp0ad/GjRvNwZ7j4sWL4sPItGnT4lqsL4Q37KWecz6Lc/5LzvlxzvkFznlvKFsE3lMseMyvv2vEUb74p7t8OEK2NEHw7PtvBOHDb7/4WGy3XMlB7fHIh+ypANwk8a5/oY96Cl9F02RpwjXHNQBAVlJWxOPVPUjXJzz5ztu47bkPo3L/4wXpQuxoPEARRKzxT9cYjoCNplDfvn179uHDh9MAIXzl7bffbgrlPNIwnxkzZsS1WC+CN3zl+xGwRclIPdlNgaRj5Jwf89sVkW8Ad5x8hv8G4P1InJ+ID2qPt+HA+U/FvrU3H5vePBFZweaO5fUJhUnSU4zvGEgX/C7IWRA1D+61XjHxA1SGZrRZBiN//+MI6WzHwlyauCTij927d/vkKX/44YdDWlwaTaHun6bx17/+daP/jECgSMN8br311usqt/pni4k14Yh1UZAqOd1ihFgsafuL8LGQjl086iiCiDAvHDgDJF0U+86B6Rh0OFF94OzoBwVLVQPwdC8WrPm9uOu32htwm/MVEoWj8JXZKwoX5EZ+wa+HfZ8wcC58/KuTOgDVQOTvfzxQXQLzMxm40n8FAJDkcqHk50uA6pIYG0YQkcNkMumkqRsLCgpsoWSCiaZQB4Dy8nIx9WNlZWV7MCkq/ZHaWFJSYjOZTLp169ZNLywsvJExtignJ+dmxtiiwsLCG91ZZ2I6LRxOAGkTKDbagzTePJgpGenvMJDc7HFJ7fE2VB84izbLINSMwck5Ptc/jmzWKy6I8PgXh/RZSNoU0qwXIaG9vw0pU/oBANyZBJddqM582TIY8Ws1X8kS26qkK2jr7cWmN08ACK8aZzzy5VWvZ/2m7Mgv+PVwxeKEwTgVakMrAEBtuARn/w1Ruf+Kpv8qvkr2ZoEptduFL01KPxp1fnOkOfOlDxryO602XU6q3v7Du0va/n7J9O5Y2xUvmEwmXXd3t+bll1/O2bNnT7Znf2pqqvPAgQPngj2f2WxW+wt1ALh48WJERO66deume4oYlZaWDoSTatFfeG/btm3yaBVQW1tb9Tt37py8c+fOyZF88AiWcMT6XriFJmPsTs75XyJjkrJwLy6V0jjiwJGRqs5ILTCVJ9Ulo37BVQD4BmNoSNfimDYNz6nvRpVBj3O6AljVKjDOked0Yp7NjrsHBvCDTbXIT09DVfkcEnshkp11FR5Z5hwshGeSbarRMOoxofKL99vgzMyFWn8VjLmgTmrD4OBMVB84S/dPQp+9D429wseHmqkxL3te1K411WiAeXC6V6wnN8PZf0NU7r/S+UqyzuKmIXsMLUkcfnOkOfOZd0zTbcMuFQBctdp0z7xjmg4AJNjHpqqqanqoFT2XLl3a96tf/aq5tLQ0qD90s9msLisrm+sv1AHg8OHDadu3b88OZRGoh5qamgzpA8W+ffuC0VnX8emnn6ZI+6MJdX/27NmT/eWXXyafOnXqdDjXD4VwwmBegTcUJmKpBxWIv8gO5oOkS9qJYEYY+dF/Ff2M4Q+TUvB6WipeMabh+UwjHs/LwaqCqfjajEL8/dTJ+F85ydBlHka9IQlWtfDnyRlDu0aDD1KS8a852Ugp3o4O12GKsQ2DhSXXxLZzSBDMBq0aVeVzIn6ty5ZBOAeniX214ZK4n/BystNb3Xt2xmwYNNETzlXlc6C2zRT7asOlqN1/pXNC781gN98m+3VoccFLHzTke4S6B9uwS/XSBw30dB8lli5d2nfo0KGGYIV6b2+vpqysbK7H6w1cn0qxqqpqeqhhJGazWf2jH/1IfPiorq4O+mHCn8bGxutsKS0tHXj11VebOjs7v+Cc1586derEq6++2rRq1SqfcCCTyZS8bNmyCY+DC1msuzO3PAwhQmERY2w/Yyw1YpYpF8q8MgrXVCo8nZOFF7Iy8IsMI36Tnoa6ZAPatMFN8Ki0vTDk/w48+3d44cCpKFkb3zg0LWLbNZSPfKMBz66eHxVP91SjAS6JWFe5vbnkxfXlS7Pv4tJoUlGWjyfuvEfsawyXse3BeTTT4YcTwEmpZ91GnvWJoNNqG1HYjbafCJ/Dhw+nFRYW3lhXVxdURVCr1aqWCvWDBw+e3rhxo3nz5s2t0nFr1qwpvv7o8bnjjjvEgkqrVq3qCcdD76GpqcmnQl9lZWX7qVOnTq9fv77Hs2C1tLTUvn79+p5333236dVXX/WJvT18+HBabW3thOrdsJIec87fZ4zdA6FC5j0ALIyx5wEchbB4MujpKs55Xzg2xYBwcsP7C/vMEfbJBsbYo3AXWZo2bdo4o68n2eUa9TU155jpcKDE7kDtwN/iV65fY7bdgTynE8MALmq1+GuyAb9JT0W3Wphp0xrrsVBXB/vTnRA/wVNyqUjJOHDOcbrLO4v30Y++E7UUgYDgxd30397POnVSG3lxR0CaCSYaBar8+R+3lqHmYia6h7oB1SAWFsu+iN+E06TVYkAl+LRyhoeR56Tf0USQk6q3Xx1BmOek6ulpaRwqKyvby8vLx9VRPT09msbGRt1rr72W4xHbra2t+vvuu2/2kSNHTKF4rw8ePHh6+fLlAwCwdevWjrfeeivTZDIlA4JHesuWLXlbt27tCPR8W7ZsyfMcn5qa6nz99dcjEi++cOHC/vT09GEAKC4uto33ALB+/fqerq6uZml40ebNmwsqKiomLBwmLLHOGPOoIgu84SBPhnFKHq5NCkNRYS+c811w55BfvHhx0MUoDZzjQes1JLs4krkLk1wuFDqGMd0xjBkOhyi49w59HbcneVPVawGUOBwo6XXgkT4rFhq/Ba1RSKTzUbIBP87JwrOdXcIiVFr4NS4t1hZYHcIieqPeiKkpU6N6vYqyfAy7VuInJ34BMCdUum5seWAGeXEluLgLJzpPiP2bc26O+jUZYyjNKkVdWx0AoUjWjPQZUb+ukvgqTQyTxU02u7fUNqUfjSo/vLukTRqzDgB6jcr1w7tLKO5xHIqLi23BZEnZunVrx7p166Z7YsKtVqt669atU4JdSCkV6h727dvXOG/evPme/rZt2wpWrVpl9R83EnV1dcnSDDXhpGn0JxTv/MaNG801NTU50ocPs9msjpRN4xFOzDogVCotApAOQWh7BBwLY1Ma/rMH4QjwuF44owHwU3M3nuruwQ97erGh14p7BgYxWyLUPXTykatp9rvSMHRlDb5v8dbP+u9JKfhNGkVgBYqpyyS2S7NKJ6Qa40OLZqA02+tJn1kg2wmkmNBqbUWfXXCGZegzUJBaMM4RkeHG7BvF9knzyTFGJiYnFv2d2J6/4l+Bp3uFjWbvosrfL5ne/W/3lzbnpurtDEBuqt7+b/eXNtPi0uiwe/fu5tTUVFF0vvPOOwEtuPQwklAHhFAS//j1tWvXBpT5Tjou3DSNkeL222/3mbEINmQoHMIV60B8CO5wGCmUJVCypJ1AiikplkA9USm5uPjc3yDnJ5e8X4zurfZbJlQY/gsAww96evG3Vu8iyf/INKIxyNj3RMVfrE8UczPnjmgDAbxWXye2+3rz8McvLk/IdedleTPO0D25nrqWerG9671hWtA+gfz9kundn21eeeLCc39T/9nmlSdIqEeXRx55pNPTtlqt6kCLAhUUFNjG8pRv3LjRvHTpUlHktra26tetWzdutpre3l7xC33nzp2TGWOLxtqkx/qPv/feeyOSGtu/eFJPT8+EiY5wL5TwdZc5501+nslgFlEkTm71CHiiKsryUVGWj9ue+xBsCPhXczfO6LQ4pdfDwRh+nJ2FX1/pgKzKjsmQWIn10qxS/KHhD9fZkOjUHm/D7748DJXblzVwbfKE5aGXivXT3acx7BqGRkUPvQCwt74R7YPNYAzgnKGjM4fqAxBxS3FxsU+qozNnzugDCVcJhLfffrupqKhovmeh6J49e7Iffvjhbjl4y4OhpKTE53fU1dU1YXIjLM865/x4pLdIvbEJRuoRD0aAS8cGU/k0ofEsTNQB2NrZDQ0Xoq++StLjnUkpuO25D8kDNgqcc5i6YyPWfYRh14SnqZUt1QfOwqX1Jk5wDeVPWDXRnOQc5CYLs16Dw4No6qWCYx5+9tePwJjw2eKyZwNcT1Veibhl1qxZUctLmp2d7XzxxRd9wmG++93vFgfqvZcLDQ0Nemk/mr8zf8iFEhmOAljpbgdT1VU69mjkzIlvKsryMbQ/C0m2LsxyOPA9Sx9eyRBi3F/MyEDb+T7ygI1Cq7UVVvvELS6VUpJRAg3TYJgP45L1Eqx2K1J1tNbgsmUAKbneh0tP3vuJykM/L2serg4IC7NPmU9hdsbscY5IDLodF+DJ7+Ya8n6OUH0AIhGItJBev359z969e3s8BYisVqv6O9/5zvR33313RA/BkSNHgpp+lS5kXbt2rXnLli1XPP3c3FwxHr+mpiZj7969Ykz+aNcfCf/87LNnzyaxrjDeg1esFzHGjOPFn49Q+fS9qFgWpyRt8v5/7X7+T3ANPwOV5ho6NSpoMw5jsPt2qpA5Aqe6vXnpJ2pxqQedWodZGbNwpvsMAMG7fuuUWyfs+nJlctYgrqkFAcidBnCH8D0yUXnob8y+EX9pEQpQn+o6hQdLHpyQ68qd1PQOONxt55D3oZbqAxCJQDTisV9//fXmoqKiNE84zP79+zNqamoy1q9f3+M/NpzCR4E63t4AACAASURBVOnp6cNjHS+tWDra9UfirbfeEtckFhQU2MItzhQMY4bBMMbumihDAoExlsYYi35Os+B536//cADHPDbOOYgAudLDYe9cKfZ1mXUAGyYP2Aic6TojtqULPicKadgNxa0L3H+Lt/6A4FVnE5qHXhqedMpMRcY85GR5s7u53GKd6gMQiUI04rGzs7Odv/71rxul+zZs2FA0keEw/sL8xz/+cUAeve3bt2d70jYCwP333x+QwI8U48Wsv88Y+4wxdueEWDMKbpH+HIAeeD3YsoFzfgyAdCplzFzzjDEj3MWF3Lw/mieeMWb038K3OL6YajTA0bsYLocQUqHS9kGbdpw8YCNwruec2L4h84YJv35pJol1f1LTvDVCXENTo1pNdiSkYv1sz1k4XI4xRicGDpcDnbaLYj8W94UgJhL/fOGNjY360caGQ0VFhXXt2rU+ec4feOCBCU22Ib1+a2urftmyZSVjja+pqcn46U9/KubTTU1NdT7zzDPt0bTRn/HEei+AxRBE+37GWHRrYPvhFukbIYj0KghpIeW6Akoq0IvclVxH45fwzcc+lrjfBOH9e7axzpuQVJXPgUGjh717ubhPn30QG79Bsbf+SMV6LGKT52Z5vfmnu2mRKeD70PKfD96PT566a0IFoTHJiCkpUwAIIvVC74UJu7ZcabI0iQ8t+ZPyceH/XzPh94UgJpIbbrjBJ/76448/TovWtXbv3t1cUFAgXu/w4cNpW7ZsyYvW9fx56aWXWqV55Q8fPpyWlpZ285YtW/Lq6uqSTSaTrq6uLnn79u3Zy5YtK9mwYUORJ3QHiGyBpkAZU6xzzjMgCEsG4B4AxxhjDYyxf2GMzYiWUYyx1Yyx38MrThmACwAWcc7fjNZ1w4Fz/gZ8M7o8wRh7RRqbzhhbyBh7D8BDknFvuD3zRIhUlOXj2dXzkcO/Du4U1n8w3VVMyw+4qnFC0GvrRceA8DvRqXSYljZtwm0oySiBigkfO5eslzA4nNihSpxzn4eWiczOI2VOpje042w3ZTuRPkDFYgaKICaa7Oxsp1TAmkym5O3bt2ePdUw47Nmzx8fxum3btgKTyeRfHzEqZGdnO48cOWKSvl+r1aretm1bwYoVK+bOmzdv/ooVK+ZWVVVNP3z4sM9DS3V1dXMsUk6Om7qRc/4YBO/6RQiiuRjACwAaGWOfM8Y2hhtHzhib4RboOxhjXQD2QRC0niJLz3POZykgtePd8E3j+CiE3xNnQg6weviG8RzjnK+ZSAPjlYqyfBx68j6smfuAuO+Nc2/E0CL5IfWqFxuLY5JP26AxYFqq8JDg4i40WhrHOSK+6RjoQPeQUOslRZuCwtTCmNghFaSeBcCJjPR3EIu1HQQRC/zjsKuqqsYtXhQqy5cvH6isrPQJJSkvL5+w6d7S0lJ7U1PTiVWrVgUUe15QUGA7ePDg6Y0bN5rHHx15Avq2dnt+ixljjwJ4At784AvdmyerxDEIKQgbIYSrWAB0wytgiyCEf2S6f94D4UFAGhIiTU+xC4JQV8S8LOfcwhibCeFhY7zY+jdIqEeeNbPXiCL9zxf/jCdveRLGJArzB3zFutSTOtHMyZyDi30XAQheXGnJ+0RD6sGdmzlXnHWYaG7I8Ip18qz7hmhJQ7cIQm5UVFRYOef1448cn927dzfv3r27ebxxLS0tJyNxvR07drTt2LEjIkVRQvkdZGdnO999990mk8mk2717d8Zbb72V2dLSordarerU1FRnenr68IoVK6xyKOAUlGuNc74LwC7G2EoIcdZ3+w0RxXsQ+OeOswB4FsAuznlvkOeKOe6Fove4f0ePQfh9eB5umiBkfXkl0NAXzvmTGGfBaoDnuSfccyiB0qxSlGaVwtRlgt1lx60vvoAcvhJV5XMSPt401vHq0msfuHgAgK9NiYiPWI+hKPQJg+k5C875hKb1lBMu7iLPOkEkEKWlpfatW7d2bN26VbaxsyHNg3PO34ew6HQmhHCVb0MQpeN9uvNRxlgA7AXwHuf8D6HYJDc8v6NY25GIFCfdDRMEEaRJ+xJtzbdRkSQA57rlI9Y9JLpYP9vj9WLHUhTmT8rHJO0kXHNcg8VmQcdAByanTI6ZPbGkua9ZXEuRbchGTnJOjC0iCCLRCWvOlXN+gXNezTlfzDlXQQhreQGCSD0G3/htwJvN5RiANyB4nos555mc88p4EepEbPlL/WQIf46AOvkSmLY74cuEO11OnLecF/slGWNmqooqczKu9+ImKtKY/VjeE8aYz0NUIofCnO7yhsDQ4lKCIORARAMkOecfcM6f4px/wy3gM90iPoNzrnJvs9yvPcw5/6VS4tEJ5dDeo4az3yt8tGlfAUjsMuGvfvo5hpxDAADmTMPHp2P3u5icMhmpOiEnvtVuFTPUJBqDw4NotbYCAFRMhZnpM2NqDy0yFZDLbAdBEISHCVnNpMTYc0K5CEWSvCUBNGlfiPsTkdrjbfiPjz4S+47Bydj05gnUHo/Iup6gIS+uQFNvEziEWYVpqdOgV0elBknASMW6VLAmGtIZqFiGixEEQXiITeoBgogiVeVzoB2aD+4SlmSok9phSO5O2DLh1QfOYlh7Wey7bJNjHhZEcevA+R6vKJxlnBVDSwQo17qANDSp2FgcQ0sIgiAESKwTcUdFWT6effAWaO1eT+H9S3oSdnHpZcsgVPorYt85NEXcHyukceuJKtblJgqLjcVQM6FI3yXrJfQ7+mNs0cQz4BhA2zVhxknDNJiRNiO2BhEEQYDEOhGnVJTl49/u8qaxN7vkXk8rSlSX4ELSOsxIOiHu+hP/T3yufzymYUE+YTAJGnLRYGkQ27MyYu9Z16v1yNJ7izLd9dKemIVKxQrpA9T0tOnQqrUxtIYgCEJg4ksYEsQEcfsfnwTL1oMzhuPtn8Py0wwYXS4gJReoahj/BPFA/1VYGUObVvhX13COmQ4HtKw3pmFBszJmgYGBg6O5rxlDw0NI0iTFzJ5YIBWGs9JjL9Zrj7eh/WomVGkXAQDdwxcTLuWpNF5dDrMdBEEQAHnWiTgm+9pVzLfZAQBOxnAw2S0G+6/G0KqJp0GnE9tFdgc8vsJYCjCDxoBM3VQAQhGar7/4u4Ty4l6zX8OVfiE0SaPSYHpa1Kp6B0z1gbOwD3pzq6v0V2K+tmGikYp1OawjIAiCAEisE3HOnQPeuOy/GhIzG0yDzjuVX+JwxNASL7XH29DZlSX2zfaLMc1QM9FIReGMtBmyCLe4bBmEy+YV62p9u7g/UfCZ7ZBBaBJBEARAYp2Ic26XiPVPDUlwxtCWWNGk9QrBWXZ5iPXqA2fhGMoT+2p9R0J5cX1EoUw8uFONBh+xrtJ3AOAJlfKUwmAIgpAjJNaJuKbE4UD2sCDRe9VqnJaEhCQKjRLPerFMPOuCF9cr1gVhmDhe3AMNX4jtv5xQyWJGoap8DpKYEdwpiHOmtsFgsCZMytM9R896C3RxNY430pIugiDkAYl1Iq5hAJYODon9w4bEWsSIlFw0ab2io9jjWU/JjZFBAlONBjhHEOuJ4MWtPd6GTy6dEvuW3kxZhABVlOXj2dU3QTM8Vdy3/utJCbG4tPZ4G3564C9i32nLwZa3TDG/JwRBEACJdSKecQvSpUNeb+0hQ1LMhepE0vvDo+jUCGJdp9Ih/1+vAk/3xjwbTlX5HCTxHHCXkNdbpe2DQW9PCC9u9YGzgK5d7DuHYl+kykNFWT4eummx2M/K7I6hNRNH9YGzGFZ7axG4bHmyuScEQRAk1on4paoBeLoXS/7xK3HX50kpWOZ6KWE8Zhd6L4jtGekzoFapY2iNF8GLezPUTm+MdOU9KQnhxb1sNUOlsQIAuEsD7sgU9sskBKjEWCK2pTHc8YxQOMz7AOWJ3ZfLPSEIIrEhsU7EPZ+ctYO7v3wZc6LDYZJF2MFE4FMlM11eC+YqyvKxavYCsT81tzeG1kwcuVkWse2y58LzMSyXEKCSjMQT61ONBqj03pSuLluuuJ8gCCLWkFgn4p7qA2fhuOYVIJqUhoSZ4m7qbRLbRcaiGFoyMtJMKNIHi3jm7pu42PYssjVo1bIJAZJmQWmyNGHYNRxDayaGqvI5YqpKAHDaJsvqnhAEkdiQWCfinsuWQQz3e0WhOvmCuD/eaeyVeNZlmIpOatP5nsTw4iZP6hTbLlse8o0GPLt6vmxCgNJfWozcYUGg2112XHo2D3g6HaguGedI5fL1uSlgmmsAAO7SYkrKVFndE4IgEpuI5aZijKUBWBnueTjnb0bAHIIQmWo0oK1vOjhnYIxDpb8CqIYwNS0j1qZFnSaL17MutzAYwNeznighF9IZhF/+3f24o/COGFozAv1XUTIpB1fdC5PPa7UocgzHdeVf6d9eafYs7F0f9lcZQRBExAhbrDPG7gKwD4AxfHPAEcEHCIIAhCnuTW+egGtoKtSGNjDGkZx6CVXlS2JtWlTpd/R7S9ozDQrTCmNs0fXkT8pHkjoJQ84hdA11oWeoBxlJ8fsQxTn3LWkv0yqZJXYHPkkW4rUbdDp8YyC+Z6Gk90Qas08QBCEHwgqDYYx9H8B7ADIgpLRmnpeC3ODXJoiIIWQemQ+D0/sl/PUF1+J+iluaCWZa2jRoVbEvae+PWqXGzPSZYj/evetdQ12w2IQFpgaNAVNSpsTYopGZJSmedV4nv7+bSEOVSwmCkDMhi3XGWDqAVyAIbO7e4O5bgtx6JW2CiDgVZfl44ZsVYr8P8b+41CcTjIwFSCItMvXxqhtnQcXkuWxolt0uts9rE0usS/8eCYIg5EA4ISeb3D85BIF+DMD3OefHw7aKIKLAwtyFYvuU+RQGhwdh0MRvajafTDDp8ssE48FnkWmce9aV8gBV7BgG4xycMVzSajDEGJI4H/9ABcI597kvJNYJgpAb4bh1Fkra9ZzzxSTUCTmTkZQhfhEP82F81fnVOEcoG5/FpTIWhtIY4Xj3rDf0eCvHylYUpuQiiXNMc2eEcTGGJq0mbiv/SkOTkjXJsg1NIggicQlHrC+WtL8friEEMREsylsktus76mNoSXSpPd6Gjy6eFPvNVybF0Jqx8fes8zj14AJQhgfXXfl3VvEqcdf5b78q7I9D/ENgGKOlUwRByItwxLqY/YVz/kUEbCGIqLM4z/uMebTjaAwtiR61x9uw6a1jcKm7AACcM7y43yLbiq1TUqaI4UgWmwVdQ10xtig6XJcJRq5i3Y1PJdM4zoGvlNAkgiASl3DEetP4QwhCXizM80Zvneg8EZfVGasPnIWNtYMxwUPN7VkYtKtkW7FVxVQ+OeDjNRSmY6AD1xxC4Z1UbSpyk+UdViJ9mDhnORdDS6KLNDSJxDpBEHIkImLdXRCJIGRPbnKuGJM65Bzy+aKOFy5bBqHSewvYOO054n65Is03HpeLTKtLcP4XC8TuLGsn2E+Msq4KWmJMjLUE0vcmfc8EQRByIRyx/oKkTTHrhGK4Kecmsf1l55cxtCQ6TDUafMS6y5Yn7pcrcV/JtP8qGiUpEIvtDnG/XClMK4RGJSQMa+9vh9VujbFFkcc/Ewx51gmCkCMhi3XO+fsAPoCQtvEFxtiCcQ4hCFmwIMf7pxqPGWGqyudAmyQV67kwaNWoKp8TQ6vGRiqS4tWL2yApLiQtOiRXtCqtT8GqeLwvHQMdsDqEhxAlhCYRhD+1tbWpjLFFoWxpaWk3z5s3b+66deum19XVJQd6zcLCwhs95ygsLLwxFLvr6uqS/e15/PHH47tSYRiEW5FjDYALEAT7h4yx74VvEkFEF6ln/Stz/In1irJ85Gb1iv1s3TQ8u3q+rCu2+nvW4zEjTKNErJfY5S/WAWBWenzPePhk58mgTDBEYmG1WtUmkyl5z5492StWrJg7b968uSaTSRft69bV1SWvWLFirnRfZWVl+44dO8bNglBTU5Nx7733FhUWFt6YlpZ2s/Sh4/HHH88Pxf66urrkxx9/PH/evHlzPef0PIjce++9RTU1NRnBnjPShCXWOecWznkxgDcAZADYxRjrYow9yxhbzRi7mTE2gzGWFugWkXdFEGMwN3MutCpBODX3NcMyFF+Fcx1OB3rsl8X+hz9aI2uhDgB5yXmYpBXSS1rtVnQOdsbYosjiAkYOg5E50rUEm//7fdz23IeyzSoUCrWnvOlbT14wxNV7I4hgMZlMyUuWLCmNpmAfSaivXbvWPJ5Qr6mpyUhLS7t5w4YNRfv3789obW3VW61WNeB96Ni5c+fkefPmzb/33nsDqgJoNpvVy5YtK1mxYsXcnTt3TjaZTMmecwJAa2urfv/+/RkbNmwomjdv3txgZh8iTTgVTMEY81+dxyCI9idCPCUP1yaCGA+dWoe5mXNFr/pX5q9we8HtMbYqcjT3NcPJnQCA/En5SNbG7PMlYBhjKDIWiWFJ5y3n4yYkofZ4GxZp1BhUCb6RTKcTWS5XjK0KDHO316Gk0negrWMQm948AQCyfwAcj9rjbfjT2S+gcruIrNasuHlvROJSWVnZXl5e3jfeuJ6eHk1jY6Putddey2ltbdV79lutVnV5efnslpaWk2MdHwqjCfXdu3c3j3Xcli1b8rZt21YQ6HX279+fUVhYeOPx48dPZ2dnO0caYzab1WVlZXOl730sTCZT8ooVK+ZWV1c3b9y40RyoLZEiXGFcDEFgQ/LTQ0LPJzLGFgJ4DMBKAEUALAC6AbwP4BXO+bFEsEGupKAIgCAM/7n2j/i324rj5gu6sdc7tV+UHpCDIfZUl2CWYRhfpQre9abfPYxlfVahaqbCi/FUHziLJ7RiWQpfr7rMq4K+/bkTcJuo0ncAAAYdTlQfOKv4/5fqA2fBje1i32XLi5v3RiQuxcXFtoqKioBXg2/durXDXwy3trbqa2pqMtavX98TKbvq6uqS77vvvtnSfYEI9Zqamgx/ob506dK+73//++avfe1r/QDw6aefprz33ntpe/bsyZa+hzvuuGP2qVOnTo903gceeKBIKtRTU1OdjzzySOcjjzzSc8MNN9jOnDmjb2ho0G/fvn2yyWQSPV5VVVXT77vvvr7S0lJ7cL+B8Ag3Zh0QRPlIW0LCGDMyxt4DUA/gUQgiGRCKSBW599UzxvYxxoyjnEbxNsiZ2uNtqDuZIvYHWBM2vXkibqbAm3q9JRAUI9b7r6JIImLPe0JGZJwtJVAuWwZRpVot9j8ZvBMzhnZj5tBu2T+IdHSlgLsEn45KYwXU/QDknQY0UC5b+sUHEMCbNSke3htBBMPWrVs7Vq1a5SPM9+7dG7E4bY9Ql4aYBCLUAeDHP/6xz5NzdXV186FDhxrWr1/fU1paai8tLbWvX7++Z/fu3c2nTp06UVBQYPOMNZlMydu3b8/2P2ddXV3y4cOHxbDr0tLSgaamphM7duxoW758+UB2drZz+fLlA+vXr+85derU6bVr1/p40v/lX/4lYC9/pAhXrC+K8LYYCsYtfOsheLLH4yEI2XTizga5U33gLAavef//1YYWDDocsi0aFCxNFq9YV1IqOmmGlCZd/ETDCak0rxeFck6l6WGqMQUuu9f7r9Z1uPfL3/bxmJw1BKYSnGOu4RRwpzCrEw/vjSCCZfPmze3S/qVLlwIKDxkPk8mkC1Wo19bWpkq935WVle1jhaCUlpba9+zZ41Ow88UXX5zsP+7ll1/OkfY/+uijc6OFywDA7t27m6UPAfv375/wBafhLjA9HuktUm8sRnwArxcbEApHrYEQLlQMISRF+oe0kDG2Lw5tkDWXLYPgjgy4hoUvZ6a2QaXrjBuPmk8YjFEhnnX4hoec1+qui6tTKlXlc6BJkor1ybJPpemhqnwOmN37XafSX1WM7eNRcauoHeCy5QJgcfPeCCJYbrjhBpu039LSErZYN5lMuiVLlpSGItQB4JVXXvER1c8880z7aGM9LF++fEA6S9Da2qo3m81q6ZiDBw+metpLly7tG0uoe7j//vt9Zh4mImuOlEiEwRAAGGMPAVgo2fU+57yYc/4G57zJve2CMIMgjRV/yB1bHhc2KAHBc8bgHJwm7lMnN8eFR23YNYyLvRfFvmLCYADkOZ1IcS+8tKpVMKvj4+Pp/gV50CZ5nUF5SdNln0rTQ0VZPspn3yz2U9PMirF9PLIzvd+9Llse8o2GuHlvBBEs/oK1sLDQNtrYQAhXqAPAyZMnxVjx0tLSgUBENQDcfffdPgtsz5w54/PgIfXW33XXXeMuxgWAW2+9dUDaP3fuXERmHgIlPr4N5cHzkrYFgjf7OjjnI732/EhjFWqD7KkqnwODVg2XRKzrki/HhUet1doKh0vwUOcm5yJVlzrOEfKBwd+7rh19sIJosbZgmLvviSEXh568X1GCsGLeIrG9oGhIUbaPhTTH+k/vXYlPnrorbt4bQQSLv6d42rRpIYv1SAh1AOjt7RXjIYOxZ9asWT5jpZ51fy+70WgM6AHAn0AfHCJFVANDGWMzICxqzISQhcQCoJtzHtCTjFJgjBXBN/TkWbcgHhHOeRNjbBeEhZ4AsJIxZhzrGCXYoBQ8X8jP/uUCPI/KBZO74uKLWhoCU5yunHh1pOQC/VdR7HDgqyTBYdGk02KpOj3GhoWP0svZS232FKyKh+JB0iJPSrwvBBFJdu/e7ROH/fDDD4eUCSZSQh0Abrzxxn5P299bPhbnz5/38XrPnj1bFO/Z2dnOzs7OL6T9QM752Wef+eRA9g8bijYRFeuMsZvhmypwtHFN8KYP/GK0cQriMb/+rgCOeQVeoQwADwd4nJxtUAwVZfm4s3Qdlv/uZQDAVdsF2J126NQTGoYWcS70XhDbihIg7qwoxadeA45uBwCcX/oosPTfY2lVRGiweDO+SIsMKYUpKVNg0BgwODwIi82CrqEuZBuuS7CgKJwup0/WJGkFXYJINEwmk06aHrGgoMAWStrGSAp1ADh06FBI6bJqamp8Yt390ywG6xU3mUy63/72t+I5V61a1aNIz7rbg/4KvBlIxnO7FMOdQpAxVg9gDec8pJspE6Tx3k2BeKc558f8vFOLRhurIBsURbo+HQWTCtB6rRXDrmE0WBowL2terM0KC6kXd2b6zBhaEhrSBwzpe1EyPiXtFSgKVUyFWcZZOGEWCgY1WhoVL9bbrrXB5hQcY9mGbBiTEi6Drbx5fuYCDHZfr08MmcN48sKXMbAo7jCZTLru7m7Nyy+/nCPNT56amuo8cODAuWDPZzab1f5CHQAuXrw4oR6wmpqaDGle9MrKynEXpY6G2WxW/9d//VfGT3/60wLp+/rZz37WGq6dwRK2WGeM3Q3gz54uri+ONObhENI1NjHGVnLO/xKuPTFCmnIymEJDx+AV2eGmrZSDDYpjXvY8tF4T/u9OmU/FlVhXlGfdjVTMNloa4yLk4nyPN9xCiWIdEP6WPGL9vOU8vjblazG2KDyksx1K/D+Je0YS6mPtJ0SqqqqmV1VVTQ/l2KVLl/b96le/ag624I+nGqi/UAeAw4cPp23fvj17Iqp+1tbWpm7YsMEnqiOQDDIeTCaTbs2aNcWAkA3H//2kpqY6//SnP52b6IJIQJgLTBljZQDeg68nnUEQgC9ACM24B4LH9h53/wX3655juLv9PmNsQTj2xBCpW6Zp1FHXIx0bbtoOOdigOKTi3NRliqEl4ePiLt8wGCXFrLvJS85DilYoWNVn74N5cMKrOkcUh9OB5j7vpKFShaH0IUMa661UpA+1JcaSGFpCEPJg6dKlfYcOHWoIVoj29vZqysrK5kozrFRXV/tESlRVVU2PdqrDxx9/PP/BBx/0qZD61ltvjZk/3Z/u7m6NyWRKNplMyf5CvbS0dKCvr++L5cuXD4x2fDQJNxuMJz+3KLgBLOKcL+acP8U5/yXn/AN3DvUP3P2nOOeLIXhxP4DXG88k51MM7oWdUoKZu5cK5ZDnYeVgg1KRivVTXadiaEn4XL52GUPOIQBAZlKmIqf2GWM+DxnSBbNK5GLfRQzzYQDA1JSp4oOI0vAR6z3KF+vS96DUByiCiCSHDx9OKywsvLGuri55/NFerFarWirUDx48eHrjxo3mzZs3+4SKeDzWkaampiYjLS3t5p07d/oUP3r11VebKioqrJG6jslkSk5LS7u5pqZmwgsiAWGEwTDG/haCJ9YT9vI853xToMdzzo8BuIcx9jyAKvfuYsbYg5zzt0K1Kwb4K6LuII7tknbCyMYiBxsUydysuWL7fM952Jw26NUTmj41YkgXzClZgBQZi/CV+SsAggd0yZQlMbYodJQeluTBfy2B0sOTzvcqPzSJIEaisrKyvby8fNzMKT09PZrGxkbda6+9luMR262trfr77rtv9pEjR0yhhHocPHjwtMfzvHXr1o633nor0xM/bjKZkrds2ZK3devWjrHPEhi1tbWpmzdvLpDGpwPC4tg9e/Y0heIBX758+cDBgwdPA0Joz/nz5/UffPBBmqdiqdVqVW/YsKGosbGxNVLvI1DCif/6tqR9LBihLoVz/qQ77t0TN/13AJQk1v2Rg9CVgw2KIFWXiulp09Hc14xhPoxz3ecwP2d+rM0KiSaLV6wrqRiSP/5x60pG6ZlgPOQl5yFVmwqrwwqrw4qOgQ5MTrmuircicLgcPoXDlPwQRRD+FBcX24LxKG/durVj3bp10z2LTK1Wq3rr1q1Tgs3gIhXqHvbt29c4b9488Qt127ZtBatWrbKGE0piNpvV3/nOd6Z7BLSUcDLPePCzzbpx40azyWTSlZeXz/Y81ETifQRLOGEw0uwjz4Zpx3OjnFcJZIZxrL+oDvVcE2IDY+xRxthRxtjRzs7OMC4pL0qzSsW2kkNh/tJ0Qmz/8TMnao+3xdCa0ImnjDBKzwTjgTEWN/elpa9FLByWl5ynqMJhCYMhczio/URY7N69uzk1NVWM7X7nnXeCCvUYSagDQspE//j1tWvXhuxJqqmpycjJybnZX6ivWrWqvLJiFgAAIABJREFU59SpUyfCFeqjUVpaat+zZ4/PWsAnnnhiQguzhCPWpaIumOwjI1Hv/smQWIsc5RBUHLANnPNd7vUIi3NycsY/QCFI49afPrAftz33oeKEbu3xNtRfOSP2uy1GbHrzhOLeB+C7MNZThEepSBdjKlmsA74zA0peZOpzTxQ82xHXPHnhSzzdW3/dRmkbo8YjjzwieuCsVqvav9LnaBQUFNjG8jBv3LjRvHTpUjEsp7W1Vb9u3bqgs9WsW7duun+ml6VLl/YdPHjw9LvvvtsU7Qwty5cvHygtLRXf5+HDh9OieT1/whHrUpEXTIz0SIR7fCzxtz0cAR7q70EONiiWnu48sa1KakObZVBxQveFA2cArTeEzmXLw6DDieoDZ2NoVWhMTpmMZI0Qhthn70PXUNc4R8iToeEhtFhbAAAMTNGhSUD8ZITxme1IJ7FOEIAQPiPtnzlzJmKLt95+++0mqed+z5492bW1tQFPad17771F/rngX3311aZDhw41hBKKYjabA34YkbJgwYKYZIIBwhPr0imBcPNzS7/FlBZvHU4oS5a0E8bCTjnYoFj2HnKBc2GxnEp/FWB2xQnd9msdYGrhs5Y7DeDOSQCAy5bBWJoVEv4hF0oVhhd6L8DFXQCAwtRCJGmSYmxReEjF+h9O1CtvBqq6BHg6HQ2Htou7ij/6X8J+gkhwZs2aZRt/VGhkZ2c7X3zxRZ8Qle9+97vFgQjmLVu25EnDXkpLSweamppOhFJhFRCEf05Ozs05OTk3FxYW3hjMsenp6T5hWKEI/lCJlFh/dNRRgfGY+ycHcDTMc00onHP/nObBrFaKiKtNDjYomSs9HC678NDOmAuqpCsAlCV0czK9z1hOWy48ZQymGg0xsig8tM4pYvsH+/YrSxS6iacQGABoaPGmnVTpO9Bm6VfWDFT/VQBAo1Yr7iqxO8T9BEF4ibQQXb9+fc+qVatEgW21WtXf+c53xgyHMZvN6m3bthV4+kuXLu07derU6WByp/szY8YM8aFEmnIyEC5evOgzPhw7giUcsf6G+ycDsIYxdmcoJ3FngnkU3hSQisu1Dl/PdjDiVzo23Lh/OdigSKYaDXANiZ8HUCe1ivuVwtfne+O6XTYhrMegVaOqfE6sTAqZ2uNt+LzBK6j6XK3KEoVu3j37hdiuM2kUZ78/L3/YAdewEJ7EVHYwrUVxM1B2AM1abxK0IocjdsYQhIzp6emJeLXY119/3Wch6/79+zPGylv+n//5n9nS/ttvvx1MwccRufXWW31CWYIJxzl58qSYJrKgoCBqMxEjEY5Y/z0EgSitQPq9YE7AGPs+gD9LzmEBsDcMm2KFdDYgmGw20rHhzijIwQZFUlU+Byp7odhXJ7UpTugmT/JW+uS2XOQbDXh29XxUlE3ogvWIUH3gLGwDuWJfpe9QnCisPd6Gjy56s/P09mUp8oFDyhXLkPggCAj3BVDWDNRFrRZOd374fMcwkhW8eJkgoklXV1fEQzyys7Odv/71r31SSW3YsKFoNC/+Sy+9JE6xrl271hwJT/Y3v/lNnxz0mzdvLhhtrJTt27dnSz3x999/f0hhOKES8pMT57yXMfYUgJ3wiu1djLEXALwCQXQ3cc7FXwxjLA2CJ/fbELzpRnjm64VzPCkdryDeA7DS3S4KpLDQCFVH34sDGxRJRVk+mvtvx68a3wYA6FMu4ye3K0voSnOs/9fffxPL8pfF0JrwuGwZBDReUajWXQXAFSUKqw+cBbLaxb50wa+S/q6kTDUaYLblASkXAABqfQec1+YqagaqUeedsZlFXnWCEPEXwo2NjVGpDlhRUWFdu3atWbpg9IEHHig6dOhQg3RcXV1dstVqFUV8T0+PesuWLXkIgX/6p38ShX52drZz1apVPZ44eJPJlPz444/n79ixY1RPSm1tbWpVVZVPyM4PfvCDCc1fHdY0B+d8F2OsGEIFUo9gzwDwpHvzVLmz4PoMJVKRzgC8wDn/VTj2xJD3/foPA9g1zjGP+fX9z6FEGxTL94/8v/g/2XpwxgDNFZS/PQ/4IwdScoGqhvFPEEM45z7x0UVGZS9DmGo0oM3CwZ06MLUdTDMApr6Gqam54x8sEy73WjBpiuB44VwFl034XlLSA4c/VeVz8K/ve9cSqHQdipuBapCI9WJ7VDO9EYSiuOGGG3zCOj7++OM0AFGZCty9e3fzwYMHUz2e6sOHD6f5Vzc9cuSIT2XS/fv3Z4xUCCkQ1q1b1yN9GPnZz37W+sknn6R5HgZ27tw5+Z133sn40Y9+1D5r1izb7Nmzbd3d3ZqGhgb93r17r7vu5s2bW6OdKtKfcMJgAAgVSAFUAuj17HL/ZJItw6/vP+6xUCugygHO+TH4Lrh9cqzxjDEjfBflvj+aF5wxZvTfJtqGRCD52lXMdAgLvV2M4aznS10Bi8+6hrrQZxcmpJI1ychLDsn5IBuqyufAoNXAZfe+D0OKWVGiMC+7V2wLi5cFv4iSvND+VJTl49El3hkbfUqnokKtumD0WVw6y+4Q9xNEopOdne2UxpObTKbk7du3Z491TDj4Fxnatm1bgclk0nn60fLsA0KRI/9wnNbWVn1VVdX0Bx98cPa8efPmr1ixYu6GDRuK/IX62rVrzdKHiokibLEOCB52ADMBPAXgAryCfDQYBHH/JIAMzvkvI2FHjJGK4yLG2PNjjP0lfGcaxhLWmwD0SLaxzhstGxKCuRJPm0mnG2OkvJCGwBQbiz2zWYqloiwfz66ejyTu9eI+cItaMaIQAO65OX4W/Er5n7csEdvapE58c8HkGFoTHAe/dQgfaKeJ/R9d+wnmOn+Pg986FEOrCEI++Mdh+4d+RJLly5cPVFZWtkv3lZeXz/a0e3t7I77AVUpFRYX14MGDpwNdKJqamuqsrq5ujlaV1PGI2C+Dc94L4AUALzDG0iHETxdByONthBAK0wXB+/u+e3zcwDl/gzF2DN4Fm0+4vdfPe1IrMsYWQhDbKyWHvuH2iseFDUqm1GbHf08S0tOd1usAa4wNCpDGXq+DQJqfXMlUlOXDol2Bn9V/BgBIkSygVQJpad5CTi5bHvKNBlSVz1HUA8dIGJOMyDZkwzxohs1pQ+u1VkxPi9r3eURZNT8LP/6qCxwA5wyTDdPwRLlyZgYIYiQqKiqsnPP68UeOz+7duwMSoy0tLScjcb0dO3a0jRYrHqgt4bB8+fKBlpaWk7W1tal79+7N/PLLL5NbWlr0VqtVnZqa6kxPTx++8cYbBx5++OGeUPO6R4qoPLm4hfgfonFumXM3hJkFj8f6UQCPjuHpPMY5XxOHNiiSUoV61qUVGYvT40OsA1B0YSSpvf97zX24Z/pdMbQmshQbi2EeFB6ezlvOK0asX+i9AO6OvpyRPg3v/M9VMbaIIAg5UFFRYa2oqJC1ey4iYTCEgDvmeyYCW6j5Bud8UTzaoFRusHnFeqNOiyGFhJNIxbrSF5dKkYr1xt5GcAWl2ZOK9XiZ7fBQYvRW/Dzfo5yHKJ97EkcPtQRBxD8k1iMM59zCOb8HwD0QCkdJF1E0QcjQsihQbzbn/EnOOZNs/hlcom5DQpCSi0mcY4Z70ZmTMTRotUI2GJnj41mPI2E4JWUKkjVCQoBeWy+6hrrGOUIe9Nn7cHVAWJisVWkxLXXaOEcoC5+HKEvjGCPlhU9F2QzlV5QlCCJxiGoAfyLDOX8fMU6FKAcbFIM7PePcj5/AxQvvAgBMq3+B+Td8O5ZWjUv3UDd6bEIonUFjwJSUKeMcoRwYYyhKL8LJLiE8ssnShGxD1JITRAypgJ2ZPhMaVXx9zM4yeoXu+V5letal74EgCELujOlZZ4w5JdswY2yG3+vn/caEuw1H880SxHiUZpaK7X8/sB+3PfehrKtO+oTApBdBxeJrskyJcesNPd68/PEoCqX35ELvBThcyiguFK8zUARBxD/jfbMzyc/RAnhZhDeCiBndPZIy90mX0WYZlHWZ+HgXINL31NTbNMZI+SC9J/Eo1lN1qZicIqRsHHYNo6WvJcYWjc+AYwBt14T/YQ3TYGbazBhbRBAEETiBuOFIQBMJwxuHXGJbpe8A2LBYJl6OJJJYV4pnPRHCLaT3pcEi7wq/gO//ybS0adCqtWOMJgiCkBfjBVP6L0Ds9us/BlD5NyJ+uNLDkJyZBZWuC4w5odK3wzVUINsy8VJvczxmuJCKXWnxJzmTCGJ9VvosfNL2CQBlLDKN5+w8BEHEP2OKdc75mLnSOecfRNYcgogtU40GdA3lQ6UTMo+ok9rgGiqQbZl4qQiJp7SNHianTIZBY8Dg8CB6bD3oGuxCliEr1maNSvdQN7qHBJ9GkjoJ+anxWXBHmk1FCTMeifAARRBE/BJfq9EIIkyqyudAZS8Q+6qkNtmWie8Z6vEVhpPiTxiqmMpnxkDuXlz/nPfxtuDXg09GGIWJ9ZKMkjFGEgRByI+wvkkYY2meLVxDInkuggiVirJ8fO+W5WI/KeUKnl0tz5Lkrx39TGzbh7Lx9hdXYmhN9NC6Jovt/2ffu7Jd7AsA53rOiW1p8aB4oyjdO4vTZGnGsuf+LOv7Is3QE8/3hSCI+CRct48FQA+AvRGwxXOuoxE4F0GEzD987XaxrU5qx9/cJL/CSLXH27Dr8CGxbxvIkXXWmlCpPd6Gow16sd/napX1+0wUD+6fT/aAOzIBAIy50D54Sbb3xTJkQedgJwBAp9KhMLUwxhYRBEEER6TmaCORMcbiPg+t/iFiSro+XQwpcbgcspzmrz5wFk5Nu9h32fJknbUmVKoPnIVtMEfsq3Qdsn6fieLBrT5wFs4hSZpTvXzvi//iUrVKHUNrCIIggkcWAZWMsXRQVhlCRpRmeYsjmbpMMbRkZC5bBqHSXxX7TluuuD+euGwZhMsmFYVXxf1yg3OeMCXtL1sG4bR5w5NU+g5xv9yQppaM59kOgiDil4DqYDPGNgK4ZYwhixljvw/RBiOAxZK+JcTzEMT/be/ew6M47zzRf3/dugtJrQs3iau4C2MDIrbjy8SJwRCfPYlwwJ5NJpNjT4Bkn53N7FkLzHEm4+xjFoN8JuNzdtYBz7EnyWadGGLLmRnbGOwkNja+ILCNLQxY8k1CRghdEEKXvrznj6rurm5a6lt1d3X39/M8/VDVeqvqlVrAr976vb/XNHWVdTj46UEAwMnekynuzZWqHYXo1wMkAL6A1qpVa2JV7ShEZ78DypMHsY3BljMEsV9Cdcnk8AcnWddQF4acQwC0pzOTC63XR7NUOwpxbnSqb98brFvx9++jPlaCIaL0FlGwDuA2ALeO8zWBFnBviKMfAkDp24fiOA+RKeoqLDyy3rQAz4314KYcrWpNvseDd+z/Cb22Mry29kiKO2euxrWLsP3pE/CMToG9sAMAUFjcg8a1N4U5MvmMKTDzHfMhkrnryTWuXYTt//qxb9+ef86yVZNYtpGI0l00aTAS4jXR16J5wfDnthi+DyJTLalc4ts+1XsKTo8zhb0JMtSN9lz/ffZcpwt2AJNlwJJVa+LRsKIGO+9YhgJV7Xvvm1+yW/L7DEi3yOB8dUD7XP7r7V8DlPbPti23Fz9tWGC5z0UpFTiPgGkwRJSGIh1ZPwgg1PKBm6GNiPcD2BdnX9oA7FdKfRLneYjiVl5QjunF09E11IUxzxja+9uxqMI6o4Yf5fmXS691WuhGIgEaVtSgN/dG/KxFK1VZXNKT4h5dqfl4J/7xzcOAngUydClzU2C8NtbPxS8/m41PLn4CiMJVc0ZT3aUr/PKtdzHoHNR2PIU4csqF9StT2yciomhFFKwrpZpCvS8im/XNo0qpH5jWKyILqKusQ9eQVru89UKrpYL19lx/sD5/LLODdSAwfcFqCyM1H+/E9qdPQGZ0wltn5Jk33bi2stNyI81mm+eYpwXr0NJNjBOzU635eCd2vfxH2PWHMq6Rqfi/nnkfIpLxnwsRZRYzqsFkbmImZbUlFf5UGKvlrZ/Jy/NtZ/rIOqAFhV7tA6Ee8qVO04FTGHaOwZZ/3vfe8NBkS5YxNJvxc7FaidOmA6fgtJ/17WdqeVMiynyRpsGMx5tfbq3/PYlMYBwltFpFmDOGNJiFWTCyPr14OgpzCjHsGkbvSC96R3pRUVCR6m4B0Mto5vVAxA0A8DjLAE+hJcsYms2Ym2+sumIFZ/uHkT/dWDFpmu99IqJ0EtfIulKqSX/9zqwOEVlF8CRTl8eVwt749Uyagl67lnBR6PGgxqX3q9h6K62axSa2gCXurZQKU+0ohC3fuEDVNN/7mc44sm6lzwTwfi6GYH1kqu99IqJ0YolFkYisqKqwCqW5VQCAEfcIbvr731hiOfUz3/6lb3v+lGtge2AAeGAAaDwzwVHpz6qBYePaRcgvMixQNTLNsmUMzTandA5yRHtAe3borK/OvBX8l9vmBwTr7rGpWfO5EFFmYbBONI7m453o7/OPVvc427H96RMpD9iNpegWli9MYU+Sy6rBesOKGtTNuezbL8uZiZ13LMuKSYy59lzMLp3t27fS57Jyngdi0546eZwlqCmpyprPhYgyy4TBuoi4DS+XiMwJ+vpHQW3ifVkjz4AI+gS1YX99b3tBpyUmqJ3uO+3bzqa60fPKDMH6gHWCQgAYUh2+7V9995tZFRDOL7dmpR5j3fubZi/Da/d9Las+FyLKHOEmmHpXFp2o4gurwVBGOts/DNsk/3/utoJO3/uplE2L7xhZdWT9svMyOga1YN0udswtm5viHiWX8XMx/m6mmnHCK1cuJaJ0FkkaDINxykrVjkJ4RvzBur2gC4AnpRPU3B53QKCaTSPr1ZOqUZij/ex7R3rRM2yNxZHaB9qhoAAAs0pnId+en+IeJZdVa+AbbxwYrBNROgs3sr4xaL83aH8LAId53SGyjsa1i7D96TF4nCWw5Q5CbGMoLOpF49rULYH42eBnGHVrK0VOLpyM8oLylPUl2WxiwwLHArzX8x4A4HTvaVTVVKW4V4FzCLLpSYeXMRC2UvnGU73+dLVsmttBRJlnwmA9XElGpdRL5naHyDq8+a0/fWsmXLnaokh33ZTa1Q+zdXKp16KKRb5g/VTfKdxQc0OKexQ4h8CYv50tZpbMRK4tF06PE93D3bg4dhGleaUp7dOQcwifD34OQEtNysbPhYgyB6vBEE2gYUUNvn/dn/n2i0q+mKB14gXkq2dRCozXonJ/2b1TfdZYifLD3g9924vLF6ewJ6mRY8uxXA38M31nfKlJc8vmZl1qEhFllnhXMCXKeHUV/pVMWy+0prAnWuqHV1YG6xWGYL039cG6UiqgH4srsi9YB7RJpt6bpzN9Z7BiyoqU9sd4A2X8nSHKNA8//HBVY2Pj7PAtI7Nu3bq+559/Pmmr0jc3N5esX7/e95hYKdWSrGunk5SMrIvIHBFZHlwKksiKjCuZftj7ITzKk7K+GEfWszENxniD8vHAx778/VQ5O3QWg85BAEBZfhmmFU9LaX9Sxfi5WGFkPdufdhBRZjE1WNcD8OUTfH2niLgBtAFoAdAmImdE5K/M7AeRmaYWTUVFQQUALRf2s4ufpaQfwSUCjakH2aI4txgVeVrte7dy4ys/S+2qssFBoUh2Fs8y1sD/ZcsbuPGhl1P6uRifdnBknYjSXdxpMCJSCmAfgNX6W7sBvBOi3dsAVuLKUpDzAOwVkTVKqT+Ptz9EZhMR1FXW4XDnYQBaKsycsjlJ78dH/R/58nBnl85Gnj0v6X1ItebjnejprYRt0lkAQI/zE2x/+gQApGTiL9MtNLN/831gahEAoLzgIxzuXg95Fhh5oRIF25P2RB0A4PK4Ap5AZfPnQtnn/vvv71i1atXl8C1DW7hwYWofV1JIcQXrIrICwCFo5Ru9CyiFarcTQL2+G7zIknd/o4gcVEr9f/H0iSgRllQsCQjWb6+9Pel9MAaG2ZgCA+iryuZMQ/4kLUC3F3RheEBbVTbVwXq25qsDwNzLPSh11+Ci3Y5Buw2dOXbMcLlRMHoh6X357KK/vOmUoim+p2JE2WDVqlWXGxoaBlPdDzJXvGkwjwEwFnq+4hmwiJQB2AZ/IC/6/hoA9wEYgD9g3x1nf4gSYmnlUt/2yd6TSb9+8/FOPPSSv1KqcbGmbHK2fxjukem+fVv+Wd/7qcB0C40AWDLm9O1/mJe6pz68gSKiTBNzsC4i34KW1uINtI9BC8AfCmp6p/cQve1GpVSTUuolpdRuAKsMX3Mwf52sqK7SXxHm5IWTUCrkQ6SEaD7eie1Pn8CIzZ8r/3yLLaU5wamirSrrD9a1VWVVSlaVHRgdQNdQFwAg15aLuWVzk94HK1kyOubbbs1PYbDeZ0hNKs/eGygiyhzxjKzfZdhuUUqt0gPwgaB2xlVQ24MXWlJKtUMbofeOyt8WR5+IEmJa8TQ48rXFegedg74FV5Kh6cApDDvHYMv313i/PDQNTQdSX7ow2RrXLkKBVEK5CwAAYh9BYeEgGtcmPygzjqrPd8xHri036X2wkiVj/mD9ZApH1vm0gyjzNTc3l3z729+ePXPmzKtKS0uXi0j90qVLl3z961+vfeKJJ+Je2vvHP/7x1JkzZ14lIvU//OEPU/4oO55g3ViKYucE7VZDGzVXAPaP0+apcc5LJhARh4isFBH+bGPknWTq1dqbvHrrZ/uHYcs/D7G5AAAepwNwF6cs9SOVGlbUYOcdVyPHNcP33vduyUt+vnrTAnz45B2+3cWfvg08UAY0ZV/tey9jsJ6qNBilFNNgiBLg4YcfrhKRehGpLy0tHbfqn1FPT4/de4yI1Le2tsb9D0Nra2veDTfcsGD9+vULn3zyyaqOjo78wcFBu/61ohdeeKH8nnvuqZ05c+ZV4YL25ubmEm/fli5dusTb56VLly7ZsWPHjI6OjnwAeOWVV1K7JDPMC9aPhWqgT0AF/KPmB8c5V7uhXVoHlHpQvEdE2kREiUifvr1HRFamqFv7oJXK3Jei62eEgGA9iYsjVTsKYSvwp7y4R6p972ejhhU1+PPl1/n2HY7zye/EUDdO5flH0hd5A9Wh7uT3xQqKp2C204Uij7YGQU+OHeftNqB4SlK70TPcg96RXgBAYU4hZpbMTOr1iShxWltb866//vq6I0eOhA2eOzo68u+5557ahx9+uCqaa3zlK19Z2NraWhR7LxMjnmDdYdjuHafNauOOUurlKM+bNvTR64PQguLN8N90OPTtzQBaRGSfiCTte9SvtTpsQwprSYV/caRkBuuNaxchv6jLt+8ZqUZhrj0lqR9WYcxFPnkh+RN+AeCkIS978ahzgpZZoPEMbA8MYNG0et9bJ+9+Fmg8M8FB5gsopVm+CDZJybp/RGSynp4e+/XXX1/nHUUHgLq6ustNTU2fvvrqqyc/+OCDE88888zpL3/5yxeNxzU2Ns6ONC3mhz/8YU2oQL2srMwV/3cQn3j+JTOOpo83Gu7Na1cYZ/Q9xPH9cfQpJfSAuAWRBcUbALwUtpV5tifxWhktVZNMG1bUoLbGPxXEYZ+LnXcsS0mpQqtIVUqS1xiAj3NDjKxnOeMNbSpuoj648IFvmykwRJnju9/97mxjoP6DH/zgiw8++ODkvffe23PTTTddrqurG2toaBh8/fXXz9x///0dxmN/8pOfhP3P8vPPP8//+c9/Pg0ASkpK3Pfff3/H+fPn31FKtbz++uvJHXUIIZ466+3QqsEAWsWXgIWQ9JKN3moxAPDbCc61xrA93ii9lb2EwBuOdmjlKb03KKv1fW+blSKyTyllnHxrOhFZDWBrIq+RTWom1aA0rxQXxy7i4thFdF7qxIySGeEPjJNHeXB+zL+wzL9svgtTi6cm/LpWNs8xD/n2fIy6R/HF0Be4MHwBlYWVSbv+qbw8uPTVSmc6nZiUxOpAVrak0hCsp6DE6Qc9/mD9qqqrkn59olQ7evRoUVVVlTvW42+66aaYF1RKlNbW1rwXXnjBNzq+bt26vkcffXTccmgPPvjguQsXLuR4g++Ojo78H//4x1MffPDBc+Md470RmDFjxujx48dPxvMzTIR4gvWj0EaJBcAWEdmjlPrU8PXH9D+9ZRlDTi7Vg/qt8Af1E43AW46IbID/pgUADiml1gQ12ysiT0EL6r1tN4jISqWU6d+vnhu/BVrqDZlERLCkcgne7HoTAHDL//M/MdV+LRrXLkroKHfHYAeGnEMAgIqCCkwpSm4esBXl2HKwuGIx3j3/LgAtLenmGTcn7frvG1JgrhrlqLqXlUbWjWsjEGWLHTt2zNixY0fMx58/f/4dqwWqDz744HTj/q9+9atPx2vr9eijj3b++te/nuwNwl9++eVSAOMG614HDhw4bbXvH4gjWFdK7RaR+wCUQVsY6ZiI7AXQBq1c4xoYAnCl1CfB5xCRrwHY493V2+8Jbmdxuwzb/QgsVemjlOoXkY3Qfj7GY4MD+5iIyD5oI/epmsSaFfJcswFowbq98DN0di9L+HL3xjSPJRVLIHLF2mNZaWnlUl+w/sGFD5IbrBeX+fthDNaTPKHSamodtci15cLpceLs0FkMjA6gLL8s/IEmODd0DueHtcnGhTmFWV/3PhLLfrGsPnwrisSJ751oSXUfMtWrr75a4t1et25dX6TB9He+853z3tH1SCalrlu3rq+urs6Soy/xzr55CP4guxzaCPke+Ms1eqOKbcaDRGSFiLihVYeZB39px0MRTkK1BL0UYkAJS6XUuDn3ek35vYa3Vps42TR4hJ8S4Ogpf+BhL9Ru7oed7oTWPDeOUDIP18+Yt24cUU2GD6Yt9G1f9Z1ngQcGtFeSJ1RaTa4tFwvL/T+bZE7ENv4OLKlYArvNPkFrIkoHPT09dm8JRQC49dZbL07U3ujaa68NSOkJVzryzjvv7Iu+h8kRV7Cur0D6O/gDdiNvoL47RABeYfi6N6gfQFBQnwa2BO3vDdkqUPCTgztDtopeO7QmQF2AAAAgAElEQVSR/eAXmeh8zzTftq2wExCtCkgia56f6Dnh215axUf7Xkt/3+jbbv34kFbnPAm1zn979Aza+vU5BErQ1pGWBawSJhU3Uc3HO9H4+3/z7Rd45iTlukRW88wzz5xWSrXE+rJaCsjhw4cDqrNcf/31EefUX3fddUPG/dOnT+eP1xYAFixYMBpd75Innpx1AIBSaqOIbIaWH20c2W0HsEUpNVHlE2/Avh/AphCrn1pdwPc70ai6l1LqWFAagymPIZVS80K9r5eTZOlGk1SXVKF/tAq2/B6IuGEv6IR7eE7Cap67PW683/O+b39Z1bKEXCcdzb14DoXlMzBss6E7JwfddjumuN0JrXXefLwTDxx4ETkztLEJ9+gU/F3zGeTZCrK6Oo/Rsqpl2HdaW9LhvfPvJfx6zcc7tVS0aZ/4/kN79f0CNFd38jMJg6kbZHV9fX0BcerNN9+8ZLy20Z4r2OLFiy0brJtShFYptVcptUopZQNQrpSyKaXmTxCotwPYDS2/e55S6s40DNQBYJVhO5qJosa2q8ZtRZbTuHYRMDrHt28r/DShNc8/6v8Iwy5t1H5K4RRMK54W5ojsYUfgqpmtSVg1s+nAKbhyP/Pte0ZmJDwNKt1cPflq3/Z7599LeInTpgOnMOx0wV7gr9Y2fKmanwlRBrhw4YJp+WxmnivZ4h5ZDxZJ0K2U+hjAfWZfOwWMz7/bx211JWPZy7ResTXbNKyowdELN+LZzqMAgEllHfi7NYmreW5MgVk2maPqwepGx3CsoAAA0Jqfh1uGE5eOBGjpTvk1/qDQPTzT9z5p5pbNxaTcSbjkvIQLIxfQNdSF6knVCbve2f5hSG4fJEd7Oq7cBVDOSn4mRCnW3d1tenAcvOhRNObPn2/ZkfNwTA/Ws4U+udSoLWTD0IyBfdokvBrSnTBr1qwU9yZ17l71VTzb+QgAoNTRiW8uT1wgEhCsMwXmCsZKLB/kJ35kvdpRiH7DCK57ZIbvfdLYxIalVUt9JU7f63kvocF6taMQ59zGG6gZAISfCVGK9fb2xh1jVlZWBuTQ//73v2+3Wl59MiR0LWYRmSMiy0Xka/qfc0QkbPmcNBEcZEezmNMF446JFWESypDutGry5Mmp7k7KzC2bi5I8rZJU70gvPh/8PGHXMub8GtMLSLPUkAZzIj/vilnuZvsPt06DLU/7q648dnhGpiU0DSpdXV0VmAqTSI1rFyG/2L8+intkBj8TIgs4c+bMhBM6IxE86TMRo/XpwNRgXQ/IHxWRM3ppxjYALdBKNLbo+3361x8VkeVmXj/FWHklS9jEhmsmX+Pbf+f8OxO0jt2Qcwht/W2+a3KRlyDFUzDb6UKpWxtk6bPb8XlOTkJrnc+u9t+Te0arUeMowc47EpcGla6MN5Ynzp+YoGX8GlbUYE7Ned++wzaPnwmRBbS1tcX9uDN4RdXnnnsu4gHfnp4ee3Nzc0lzc3NJcFWZdGNKsK6PmB+AFpBvhlY7XSZ4zdPbtYjIWyIy24x+JFlFHMcGB/bxnItSYPlk/33m8e7jCbnGBz0fQOljxfMc81CUm9b/1piv8QxsDwzgmlm3+N56598/kdBa5+90+2/M/mLFTXjtvq8xKAzBmLLVeqEVTrczMRdqWoDRn5aha8hfIvL5wZ+i4dAtibkeEcG7Kmg4v/jFL0x5BF9XV+cL2B955JGIqyx84xvfqF2/fv3C9evXL9y6dWta/0Mdd7AuIrdCGzFfDX8pxogPh1YNpV1EvhpvX9JIWqS90PhW/uHvfdtH338yITW+3+sxpMBUMQVmPMun+G+cjMF0IhjPv2LKioReK51VFlaiZpL2f+OYZwyn+04n5kJD3WjNy4NTL4c72+lEpceT0PKdRNkoOHc83AJDTzzxRLlxMaN43HvvvV94tzs6OvIffvjhqnDHHD58uMi4aummTZt6zOhLqsSV/C8iK6CluAD+RZEEWmnCQ9CC+HZo+dwV0CqfzIMW2K80HCcADonISqXUu/H0KYmCc9TjCcCjyXcnC7h64Avklc3EmE3wSV5uQmp8Hz131LdtTLuhQManHIlKSQIAp8cZcAPFYH1iV1ddjc5LWi75ez3vJWxBr3cK/PHANSNpW+yByBRHjx4tincC5uLFi0eDzxG8wND3v//92a+//nrIx5iHDx8u+tGPfmRaxsTdd9/d95Of/GTUG/w3NjbOrqysdN99990hVxxtbW3Nu/322xcGn8Os/qRCvDN19+l/+gJuANuUUuPlBfjqrovISgC7ANxqOH4fgIWhDx2ffq67kJgR634Ae5RSwaUZ40llqTTuRLKYEllLvgKuGR3F24Va2cC3CvLx74YiXlgtLLfHHTCKu2oay/GP56qqq2CDDR54cLr3I9yw69+w9bblpqennO497at5P714Omveh3H15Kvx/CfPAwB++uJz+O/NU9G4dpHpn8vxfH+wvmKUwTpltx07dszYsWNHXOdoamr69N577w0Yia6rqxsrKSlxe1Ngjhw5UnrDDTcs2L17d6c3r7ynp8f+D//wD1U7duyYAQAzZswYBbTR8Lg6BODJJ59sNy6IdM8999Q+9thjFzdt2tRz3XXXDU2ZMsV9+PDhoqeeeqriySefDBh5f/zxx6MprW1JMQfrIvItaCPl3hH1XUqp7ZEer5Q6BmCNiOwC4F03fJ6IrFdKPRNld1YD2BrlMdFoA7DX+IZSqj1oJdKQK4iOg7XVM8CXRkZ8wfrbhQWmBuun+k5hyKkNZEwpmoIZk2aYdu5M8+L7fXCPTIcUdEJE4dzYaWx/WkupNDMwNM5NMKbeUGgDff6fvb2oHZ1nL2srjcK8z0UBeNcwsr6CI+tECfPII498es899/jilyNHjpTefPPNpQBgDOS9+wcOHDi9cePGaGKjcd10002Xm5qaPm1sbPSN2B85cqTUmOoSyv3339+R7qPqQHw563cZto9FE6gbKaW2IXBFzz+P4Ry7lVKSwNfecS5tHBGPJgA3to1m5VOykGuH/YHBWwWmpOYB0JZP/8tf/8a3PzV3CYJuDMmg6cApOC/76/7bCz9LyKqixmCdKTDh/fpVF5Rb+3thy70Iye01/XP5LCcHvXYtPih1uzHX6TLt3EQU6O677+67//77O0J9zRioz5gxY/SNN95oraurGwvVNlb33ntvzzPPPHPaO2I/kRkzZow+/vjj7Q8++OA5M/uQKvGkwaw0bO+Msx8PAXgqxHmt7ii0UX0gun4b2x4dtxVZ2rLRURR4PBix2dCRm4suux3T3fGt1dB8vBPbnz4BNeU0cvX33jlTgebjnaw6Mo6z/cOwl84GcAQAYC/81Pe+WZRSAcH6yinp9M9UanT1j6GgZA5yJmnBub2oHa4Bc1cWPV7mLzZxzeiYf/QpgeU7iazk3nvv7QlOWUmkBx988Nzf/M3f9Pzt3/7ttFdeeaX0888/zx8cHLSXlJS4r7rqqqFNmzb1GEey//SnP/lml4fKpW9oaBhUSrVEev2GhobBhoaG95944onyp556qvz9998v6ujoyC8pKXHPnDlzdNasWaO33nrrxUh+JtFeO5XiCdaNOdrxjg57f1iC9EoROQh/sF4rIo5w+echVj49GLIhWVvxFOQNdWP56CjeKNRWSnyrsADfVMVxnbbpwCkMO10oLvrE997I4Bw0HTjFYH0c1Y5CnL3kn8tkL/oEgAvVjhLTrtF5qRPnh7Va3pNyJ2G+Y75p585U1Y5CnL881xes5xR9DNfAl0xdWfToig1A27MAgBU3bgOu3mTauYkotKqqKvejjz7aCaAzkraJ6MPdd9/dlwnpLZGKJw3GOJkz3mom6VoN5VDQ/p0RHLMlzDkoHTSeAR4YwLU3bPO99cZ134u7xvfZ/mHY8s/BlqPlqytXETxjk00djcw0jWsXoUCq4BnTxg/E5kRhyVlTV7A0Vua5evLVsNuychG9qDSuXQT7qP+mxl7UburKokopvPXFW779a6dfa8p5iYisJp5g3Ti7Nt5SFcbR5rSpjKJPkjX+HLaN1xYARMQBbTEor0PjjcSLiCP4FX+PyWw3VN/g23797OvwKE9c56t2FMJe7K9J7bo8H4DN1NHITNOwogY771iGfJc/CLytftDUJxFvdr3p2752GoPCSDSsqMGDt68DPFo5ZlteP7b9uymmfS4dgx3oGuoCABTnFnOFXyLKWGYF65vHbRUZ72izQvrlcBsD9Fq9us14HkPgE4mJgvvtAPoMr4nOSymypHIJKgq0Ed3ekV4s+Ole3PjQy2g+Hvbp4JWaFuC1kfX4aunvfG/9t9E/4O38H5o6SpyJGlbU4MG1Db79izhp2rmVUgHB+vXTrzft3JnuWytn48s19b79yqoY/l6M480v/J9J/dR65NjirURMRGRN8QTr+/U/BcDGWFcg1VdA3Qx/Cch9EzS3HKXUfgTm7G8VkT3G3HQRWSkiBwFsMLTbr4/MUxqziQ2zC/2TDe2TTqGzfxjbnz4RfcA+1I0REbTkF/jeumF4BJNlgPnqETCmQbx7/l2MuEZMOe/HAx/78tVL8kqwuGKxKefNFsY1At7oesO0877VZUiB4dMOIspg8QTrv4WWsmJcgfSvojmBiGwC8KLhHP3wV4VJJ7ciMH1nM4A2EVEioqBNoF1t+PoxpdTGZHaQEufDdn8g7Z1MF2uJupaCfIzZtDKNtWNOTIuzukw2qSqswrwyraSv0+M0bTVTY4B57bRrma8eJeOTiNfPvg6l1AStI6OUChhZv276dXGfk4jIqmIO1pVSAwDugxZke4PtvSJyQUT+m4gsF5GAYvUiUqq/v1NELgD4uX6c9xzblFIXY+1Tquh553MR2WTR/Uqp+vDNKF2c754NpbQA21bQAbFfAhBb6cDXC42j6pxUGq0vTfuSb9s48hoPYwoMg8LoLa1cirL8MgBAz3APTvedDnNEeKf6TqF3RKtLUJZfhoXlUS98TUSUNuIZWYe+WFATAgP2cmi52C0A+kTErQfwbmi51y3QVhstNxwHALuVUv8UT39SSSnVr5RaA2ANtBQhY05/O7QVUOsjHVFXSm0LWpgpuIpMpP1aYzgHbxISoLq0Eu5hrXSgiEJOSav2fpSTQhWAPxb5j/nysDlpHNnEGEy/dva1uM+3v+UTvPzp6779S31z4z5ntrHb7Lhhun8i9uHOw3Gdr/l4J/7if/2zb39WwXLYJK7/yoiILC3uf+H0FUh/AGDA+5b+pxhe5UH7we22xLoCqtUopQ4ppTYqpeYZguR5SqktzFHPTI1rF0GGrvbt55SciKlE3Ue5ufgsV1sKqcjjwXUjDNajde30a5Ej2kTD1gut6L7cHfO5mo934icH/gWwaZ+DZ6wcf/9cX2yTh7PcjTU3+rbjuYnyLhp22X7C9947p6r5mRBRRjNlOEIfYZ8LLS3mY/gD8vEItOB+G4BypdRjZvSDKBUaVtRg283+ByY5xW348TdmRT0p9JCjyrd98+Vh5HtvZ7kaY8RK80qxcqp/wu8rHa/EfK6mA6fgLvzAt++6tBjDTk9McxGynbHE6fFzx3Fp7FJM52k6cAojngHYCj8HACgluHxxPj8TIspoptW60nPYdwPYLSJl0CZU1gKohFausB/ABWgpIYf09kQZ4XvXLcehC8u1SY3iQX5ZK4C6qM7x8qxlQO+HAIBbb/9HYO7XE9DTzPeVGV/xLZbzp8//hA0LN4Q5IrSz/ZdRNM9fAtJ1aYn+PucSRGty0WQsrliMD3s/hEu5cLjzMNbNXRf1ec72D8NeehravH3APTwLcBfzMyGijJaQRD+l1IBS6ndKqSal1H1KqR/ofzbp7zNQp4xz25zbfNsvfPJCxMc1H+/E9Q//L3yoB+p2ycXNNTeb3r9sccvMW3zbb3S9EXMJx6lVF2HL0yYxKnce3Je1aqxcoCo2t8661bf94qcvxnSOakchckr8N1Bu/QaKnwkRZTLOyiEyyZrZayB6BtibXW+iY7Aj7DHeHNxe8U9idA4uxqEPeD8bq1mlszC3TJsIOuIeibm295eXdfm2XUMLAZUT01wE0tw2238ze7jzMC47L0d9jh+tmY2cSR/69l2DS/iZEFHGMz1Y18sz3qGXZ/ytiBwQkbf11wEReVT/emn4sxGlj2nF0wIm0j195umwxzQdOIVhpxO5Zcd97430rWAObpy+OtO/RttzHz8X0zk+G/XfQLkH61DjKMTOO5ZxgaoY1TpqfXXwh13DMU00La04A7E5AQDukamYXjSHnwkRZTzTctZFZDmAXQhc/Gc8m/VjDgLYpZT6g1n9IEqlDQs2+ErT/dM7v8XPnpqDasckNK5dFDKgONs/DPukU7DlassLeFzFcF9ahLNgDm48bp97Ox5//3EAwB8++wOGnEMozi2O+Pi2/jZfPfB8ez6O3PufMSlvUkL6mk3WzFmDtnfbAABbn/8V+j8ZQ7WjcNy/H8Ge//h53/Z/un4DfnDN1xLWVyIiqzBlZF1Efgv/Kp3G8ozAlSUbje+tgbby6W/M6AdRqv3ZzD/DpJxyAICyX4S97B109g9j+9MnQpaXq3YUIq/iVd++a6AegJ05uHFaWL4Q8x3zAWipMMub/h43PvRy+BJ/TQuAB8rw3D/7R+b/7GIfJj2yIpHdzRrGVBhb/jG8U/QXeG1kPRqercPIztoJj+0f6cernf6/K+vmRD9BNVOZsSosEaXWRH+P4w7WReRtABsQun76AIBj0Fb2PKbvh2q3UUTMWW6QKIVybblw9vrL1OVV/hGAB8NOd8jUlm/fbENOsbZ+llI2jPXewBxcE4gI5hTc5NvPCXPT5DPUDTeAfynxj8LffmkIGIq9Xjv5LShfgLrRUQDAmE3wr5OKfF8rGL0w4bHPtj0Lp0dLgbmq8irMKZuTsH6mExG57Ha7Of+MKM253W6biISczBPXX3AReRSAd1VM7wqmHwPYqJSyKaUqlFKrlFK36X9WKKVsAO4EcByBK5/Wi8j/iKc/RFZwoWsVlDsfAGDPP4+c0ncABJb8az7eiRsfehn//fg/+t5zXbwa1ZOmMwfXJG+8N9u3bS8+A8ntHfemyejVwkJ05WgZguVuN24eZkqSmb41OOTbfnrSJEQyJqyUwv7T+337sZbjzFDnRkdH81LdCSKKj/73+Fyor8UcrIvICgBb4A+2BVqQPl8p9buJjlVK7VdKrQJwGwID9i167jtR2qourcBYn390fcbUJ/Fe0XfwccG3gQfKgAfKcGPzl/GF813klGiVLZQS/M2qzXjtvq8xUDfJFxcK4bq0AAAgopBXfgTA+HXSvSPuvy3156Y3DA75F6ciU3z90hAKPB4AwOn8PLxVkD9h++bjnbj+Zz/HJxc/AQDk24rwda5B4ON2u5+9ePEi8+aI0tzAwECh2+1uDvW1eEbW79L/9Abb9eGC9GBKqUMAVhnOAWg3AERpq3HtItgHvgaPUyt4dCHHjp2VFQEjiLn2QRRM948UugZW4levOJPc08xW7SjEWJ+/Ok+u421ARkPOB/CW0DyTm4vXCgsAAKIUNg4OJq2/2aJEKXzjkn90/Z8c4xcG834uFwv8E0uHe5fjxff7EtrHdOLxeF7o6+tzM2+dKH0ppdDf3+/2eDwhF2mJJ1j3PodUAPYqpY5P1Hg8SqljAB6Df3Q+kmoyRJbVsKIGO9d/CcWX7vC996+TivFwhQNOAD12G/7D1MkBFWBGu9dxFUaTNa5dhLzRJfCMVgEAxD6C7VN/hNdG1vuecKBJG3nXSmi68Wh5GZRo02q+cnkYM13ulPU/k/0fA4Ow6cHlG4WFeCc/RBZH0wI0PFuHXxX/JXKKPwYA5CiF31x8gaVNA51yOp1t58+fL0t1R4goNt3d3WVOp7MNwOlQX48nWK8wbO+L4zwA8JRhe+KSAERpoGFFDd780f+JOwYv+d77ZVkpVs+qwe0zqnHC8Oh/pOtbUO4SVoAxWcOKGuy84xoUDfvv/58oK0WvzfDP3lA3VvzXF9HZPwxbQScOFvsnPP6w37AwVfGUZHQ5OxRPwUyXC+uG/POodlaWYyi/Ejc+9DLm3vdvWP7TF32TfZsqHL52//ulISxz9/PG1qC+vl65XK7vdXV1ne3u7i7jCDtR+lBK4dy5c2VffPHFWZfL9b36+vqQf4HjCdYdhu32OM5jxvFEliMiuL+nF7cYgpJeux3DhmBxpKsB7kt1rACTIA0ranDkr7ehdkxLMRqy2fB/VzgC2vRddgJwo2DaM773nIN1uP3iL3BjwTPAAwNA45lkdjuzNZ4BHhjAf/zLPyHPpo2ot+bn40u5/x6d/cNQAPqHtc/rf5aWoDVfu7HN8yhs0m+geGMbqL6+vtvlcm3s6ur66OTJk6Vnz54tHxoaKnC5XDYG70TWoZSCy+WyDQ0NFZw9e7b85MmTpV988cVHLpdrY319/bhlx+JZFKkd/lFwx0QNo6DAwJ0ySB6Af+juwS/KSvDPZaXos9sBAPPGxvBx31/D1T8XNVEsCkPRy7Hl4D/39uOvp00GAPy+ZBKuHRnFNw150/lTnoe9sAMAoDw5GO3+Om+gEmxmyUx8f9n38T/e1YqA2ar+FfbhKriHFgIAXikswM8MN1bfHxjwpSXxc7lSfX19d0tLy+0ul2vh6Ojoup6engYAU5VSRWEPJqKk0csznnO73c16jvrp8UbUveIJ1g9BX4kUwK0A3onjXMY6XMfiOA+RtRRPgX2oG/cMDOK7A4P4PDcHBUqh3O5A4fb/mOreZY1bhofxv10awr9N0uqn/11VBS7Ybfjq5WHkO5qRV/6Gr+1Yz62oLprFG6gkuGfZPfhjxx/ReqEVIh4UzvxnOPuuBQT4kWMy3Pr8gatHRvFX/Rd9x/FzCU3/D/+U/nokxd0hIpNIrI/I9NKNLdBGw/sBzFFKRV06QUTKoNVmd+jnWqOUejmmTlHSrFq1Sh09ejTV3SCKzANlGBTB96qn4kze+CWpnYN1qLi0Ca/fx3nuydJ1qQtrn7oTyt4f8uvVThd+0XUO09yGyb4PDIRsmw5EpEUvXUxEFJGYc9b16i/eKi4OAC+JSEk05xCRUmgj9N5A/RADdSIyXfEUlCiFPV90+1bQDOa8uBTS/R1sXbskyZ3LbtMnTcd/Wfr/Qg3PvuJrK0dG8MvgQJ2TfYkoy8Q8su47gUgLgBX6rgKwB8BupdQnExwzB8A2+NNoBEAbtFrtF8c5jCyEI+uUrva3fIKHXnscl/Pehj33IuCajJGeazHZfh22rl3MFIsUeeZYBx7607Po93yI/BwbZHQ+LvbORbWjKKNSkjiyTkTRiicNphTaBFOBNsK+Uv+S94Tthlc/tNHzWsML+rHeY/ZCC9jD6VdK/VNMnSbTMFgnIooeg3UiilY8E0zvAvBzw743SPcG4PMQuma6GLaNdwqbgxuO4xgAButERERElPHiCdaB0IF38FB9cJvxhvJlnPeJiIiIiLJSPMF6L1JTE5112ImIiIgoK8QcrCulfgfgdyb2hYiIiIiIDGIu3UhERERERInFYJ2IiIiIyKIYrBMRERERWVTciyJRdhKR8wA+jfHwKgA9JnaHKBb8PaRUmK2UmpzqThBR+mCwTkknIke5KAilGn8PiYgoHTANhoiIiIjIohisExERERFZFIN1SoW9qe4AEfh7SEREaYA560REREREFsWRdSIiIiIii2KwTkRERERkUQzWKaFEZKWI7BGRNhFRItKnb+8RkZWp7h9lLhE5qP/ORfvqS3XfiYiIvJizTgkhIg4A+wCsDtN0P4BNSqn+xPeKsomItAGojeHQfqVUudn9ISIiigWDdTKdHqi3IPJA6ZhSqj6BXaIsJCKx/uPGYJ2IiCwjJ9UdoIz0EgID9XYA2wAc0/dX6/veNitFZJ9SamPyukiZTL9hNNoN4EIq+kJERBQPjqyTqURkA7T0F69DSqk1Ido5oAX1xrz1eqXUseC2RNHS50O0ePeVUpLC7hAREcWME0zJbLsM2/0AQo6W6znqwV/bFaotUQyMT3Y4H4KIiNIWg3UyjYjUIjBI2jnRxFGlVDsCV5FcHSJ9gSgWXzJst6esF0RERHFisE5m2hK0H8ly7nuC9u80qS+U3YLnTBAREaUlButkJmP+eXsk5RhD5KizKgyZgcE6ERFlBAbrZKZVhu1oJooa264atxVR5IzB+tvBX2S6FRERpQsG62QmYwAUzWimsW0si9gQBQv4XRSRWuNKugD69NVK20Rklz7fgoiIyHIYrJMpQgQ7bVEcbgzWOeJJcQnxu7gd2u/jZlx5M1gLYCuANhEJnj9BRESUcgzWySzBQXZvFMcGLFbDFAWK08qg/Q0RHrdZRFrCNyMiIkoeBuuUKKxtTakSKqXlGLS6/uX6Aknz9P39Qe1WisjBBPePiIgoYgzWySwVcRwbHNjHcy6ieUH7u5VS9Uqp/d4KRUqpdn1/I65cnGu1iKxOSk+JiIjCyEl1B4jAPHUyVwuA3fp2m1Jqwnr/Sqn9IrIFgTX/d4FlRImIyAIYrJNZgnPU4wnAo8l3JwoQLjgf7xg9YPfmu68UEUckawUQERElEtNgyCzxpLJUGncYIFGKHAraZ81/IiJKOQbrZAqlVHBd9eC84YmwxjVZQfDiSUzPIiKilGOwTmYyjohHE4Ab20az8imRmYJvODnRmYiIUo7BOpnpqGE7uNb1RIxtj47biiixgm8wo1mFl4iIKCE4wZTMdBCAt+RdbSQT9EKsNska1xQzEdkA4C7vvl6aMVIM1omIyHI4sk5mCp6gd2cEx2wJcw6iaG3wvvTgPVJ3GbbbQ8zDICIiSjoG62QapdQxBI5GbpuovYg4AGw2vHWIlWAoHkqp4BVJd0VynIhsRmA6VvB5iIiIUoLBOpnNGKDXishEwdJjCKy4MWFwTxQhY531WhGZMLVKH303/p72A9iZiI4RERFFS5RSqe4DZRgRaaQhSdIAAAkzSURBVEHgKOVeALu8aQUishJacGRc0n1/lPnFRCHpT2w+RuCNoDcAP6RvO6DVUd+IwN9DAFijlGI6FhERWQKDdTLdOMHSRI4ppbi0O5lGn7jcguhrpW+JZQVUIiKiRGEaDJlOzzufi8gmi+5noE5m05/izEXkueftAOoZqBMRkdVwZJ0SSkRWQ6v4shL+0njt0AL5PfqkVKKE0UfZvSUda6GNtvcD6IX2e7iPaS9ERGRVDNaJiIiIiCyKaTBERERERBbFYJ2IiIiIyKIYrBMRERERWRSDdSIiIiIii2KwTkRERERkUQzWiYiIiIgsisE6EREREZFFMVgnIiIiIrIoButERERERBbFYJ2IiIiIyKIYrBNlMRHZICLK+0p1f9KViPTpP8PVCbzGSv0afSLiSNR1iIjIWhisExHFQUT2AHAAOKaUOpSo6yiljgE4pF/rsURdh4iIrIXBOlEGEpGtInJQf21NdX8ylYisBLBZ392WhEt6r7EhkaP4RERkHQzWiTLTPACr9deXUtyXTLZP//NQIkfVvfTR9f367p5EX4+IiFKPwTpRdjsELbD3vihCIrIZQK2+m4xRdS/vtWr51ISIKPMxWCfKYkqpfqVUu/eV6v6kGW/Q3K6PeCeF/jl5r7c9WdclIqLUYLBORBQlPV/cO6qeinQU7zUdIrIhBdcnIqIkYbBORDEREYc+kbXFULqwTUT26RMvg9vXisguvU3Y9hFcf7WI7Ak6X4t+vkQHsMa0l/3jtgqi/7z26f30lmGMpc9PGbY5uk5ElMFEKZZWJsoEev7yrgialiul+vVjVsMwMqyUuiJvXQ+kW/TddqXUPP24fdDKCI5nt1Jqm36OzQg/Ar1XKbUlXOdFpFY/V7hqKO0AtimlIg6mI6Ffv03fPaaUqo/gmK3Qgupw9dGPAdgUSVqNiByE/2dQn8xUHCIiSh6OrBNlNwe0dA7vKyw9UD+I8IHnVn3Rpa2ILFVksx7UT3TtWmg3DpGULawFsC/cOWNgHAH/bbjGIrIP2k1UJAsZrQTwUoSLHu0zbIe9ySEiovTEYJ0oc7RDq+5yCEC/4f1+w/vxlhesgBaoA9oo8BqllCilBEA9/BMfvbyBqrd/xvbzQvRn3CcDegDbgsCg9xi0QLVeP9+aEOfcY3JazF2G7Ql/nvqNgvHa/QB2Q+tvObQ+b0Tgz82BwEB8PMZrs+Y6EVGGYhoMUQbSV9X0jijvV0ptHKfdBhgCQz2IDm5jTIPxOqSUWjPOOdtw5Sj9RO1boI0oe/nSdILa7UNg4OtLswnRNjglqD1Uik8sRMT3j2aon1dQW+P31g8tXSVk1Z3gn0O4c+vH9MF/8xLy50ZEROmNI+tEFIuQwb8u1Oj4RGkawQH3quAGevqLMVDfP16gDgBKqd3QRrC9TKlJHrRqaCQ54sabkL1hymMGfD8RTro9atjm6DoRUQZisE5E0dodZgT3aNB+uCA1uH2ofO3gwHzTBOcDAOjBvLGfIUf2o2Q8R3C/A4TIO28L2TDwfLsNr0gYbxjM+P6IiMhiclLdASJKO2+H+XpwIB+cQhNAKdUvEjbjwzhqvD+KdI+9ALwj6maMPBtHuycMvkN8Xxv1/ozbHtGvhGr8LK54IkFEROmPI+tEFK1w6R+9QfsTjkCHo49QG3PgD47XNoSAGws9nSYeFYbtSFZ8NbZZrddTj6TSS6SMNy3xfm9ERGRBDNaJKFrBwXiiBY8YRxP8B99YxBvQGo+PZHQ/OFd/AwDvQkhbY1kMKojxszDzJoCIiCyCwToRWV1wEOpd/TPsC1emqsQb0BqPD3vTopQ6hND55yuhTcT1fi8HRWRzDKPuATcMJo/aExGRBTBYJyKrqwjfJPHnChEIR5Q3r090XYOJ04e8K8n2iUgkq9B6Bd8wmPmzIiIiC+AEUyJKN/Es7BRJnnlIEU6EHe/YQwDqDSUo70LgZFWjrSKyWilVH8Olkp2iRERECcZgnYisLjgA3WiRxX9qEWXwr5ew9JVm1Ou2r4EWwBvz4VeKyK6JasnrAkbSLfJzISIiEzENhoisLjggTmWqh7EvceeHK6UOKaW26aurBi80FckiTtFOeCUiojTDYJ2ILE0pFZzrHXG9dBFxiMhq/RVv5RUgMCCe8KZBRGqjubZSaj+uXMU0XPWaqCa8EhFR+mGwTkTpwBiwR7Nw0D5oddkPQqu+Eq9oRta3GK494cJQBsE3JuGuEW3ddyIiSjMM1okoHew0bNeKyOZwB+ij2cZR+D0m9MO4INOaMG2DF2SK5IlA8Ah8uNQW4yTUaBaLIiKiNMFgnSjzpX3tbT1FxDhyvEdENozXXk8feSnEOeJlrEQTvFhTsOBR8j0T1UHX+7zd8Fa7PiF1IsYbgHiq5BARkUUxWCfKfKsyZLGc4AmY+/TFhDbo+eHe/PQ90BZDckxwbEz04Nk72u2YKKdcb2u8QaiFXkdd72etIa89VJ8jSdvxXT9Ebj8REWUAlm4kykzGlTsd0ILEfn27PB1L/CmljonIFgSms6xG+Amn20waVfd6CoA3DWclJs4V3wStf8YgfCvCV3rZr5TaO1GDoLQajqoTEWUojqwTZaZQwWnaj67rAewaRDaZsh1aTfbdJncj4rx1/aboVkQeTPcD2KKUiuRJgPHa+yI8PxERpRlRSqW6D0SUAHqKxi4Ejuy2A6hPx5H1YHrOuncl0FpogW67/joYbmQ6zmv3QfuZ9iulyiM8ZjW0CjEroVVxcSCOPotIG/Q0GKVUbEurEhGR5TFYJyKKkohshT+nfKPJaTaRXH8l/OUg9yqltiTz+kRElDxMgyEiip5xBPyuFFzfGJybUT+eiIgsisE6EVGU9DQib8C+IQXVdu7U/zwUQXlHIiJKYwzWiYhiYxzRDrtIk1n0BaG8NwfRrOZKRERpiDnrREQxEpFd8JdhTEpJTMPkVuaqExFlAY6sExHFSCm1Df4yktsnamsGfWKrtwoNA3UioizAkXUiojgEVWZJ6Oi6YVR9jVKKCyEREWUBButERHHSa77XQpvweSxB16gFsAHaqHrCasgTEZG1MFgnIiIiIrIo5qwTEREREVkUg3UiIiIiIotisE5EREREZFEM1omIiIiILIrBOhERERGRRTFYJyIiIiKyqP8fh4dHvt4jRNIAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"\n",
"r23=solve_ivp(my_ode,[t[0],t[-1]],[x0, v0],method='RK23');\n",
"r45=solve_ivp(my_ode,[t[0],t[-1]],[x0, v0],method='RK45'); # default = 'RK45'\n",
"plt.plot(r23.t,r23.y[0],'o',label='RK45')\n",
"plt.plot(r45.t,r45.y[0],'s',label='RK23')\n",
"plt.plot(t,r[:,0],label='Euler')\n",
"plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
"#axis([0 10 -0.2 0.3])\n",
"plt.xlabel('time (s)')\n",
"\n",
"plt.ylabel('position (m)')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}