{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###### Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 L.A. Barba, C.D. Cooper, G.F. Forsyth."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Riding the wave"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Convection problems"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Welcome to *Riding the wave: Convection problems*, the third module of [\"Practical Numerical Methods with Python\"](https://openedx.seas.gwu.edu/courses/course-v1:MAE+MAE6286+2017/about). \n",
"\n",
"In the [first module](https://github.com/numerical-mooc/numerical-mooc/tree/master/lessons/01_phugoid), we learned about numerical integration methods for the solution of ordinary differential equations (ODEs). The [second module](https://github.com/numerical-mooc/numerical-mooc/tree/master/lessons/02_spacetime) introduced the finite difference method for numerical solution of partial differential equations (PDEs), where we need to discretize both *space* and *time*.\n",
"\n",
"This module explores the convection equation in more depth, applied to a traffic-flow problem. We already introduced convection in [Lesson 1 of Module 2](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/02_spacetime/02_01_1DConvection.ipynb). This hyperbolic equation is very interesting because the solution can develop *shocks*, or regions with very high gradient, which are difficult to resolve well with numerical methods. \n",
"\n",
"We will start by introducing the concept of a conservation law, closely related to the convection equation. Then we'll explore different numerical schemes and how they perform when shocks are present."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Conservation laws"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You know from (non relativistic) physics that mass is _conserved_. This is one example of a conserved quantity, but there are others (like momentum and energy) and they all obey a _conservation law_. Let's start with the more intuitive case of conservation of mass."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Conservation of mass"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In any closed system, we know that the mass $M$ in the system does not change, which we can write: $\\frac{D\\,M}{Dt} =0$. When we change the point of view from a closed system to what engineers call a _control volume_, mass can move in and out of the volume and conservation of mass is now expressed by:\n",
"\n",
"![massconservation-CV](./figures/massconservation-CV.png)\n",
"\n",
"Let's imagine the control volume as a tiny cylinder of cross-section dA and length dx, like in the sketch below.\n",
"\n",
"![1Dcontrolvolume](./figures/1Dcontrolvolume.png)\n",
"#### Figure 1. Tiny control volume in the shape of a cylinder.\n",
"\n",
"If we represent the mass density by $\\rho$, then mass is equal to $\\rho\\times$ volume. For simplicity, let's assume that mass flows in or out of the control volume only in one direction, say, the $x$-direction. Express the 1D velocity component by $u$, and conservation of mass for the control volume is translated to a mathematical expression as follows:\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial}{\\partial t}\\int_{\\text{cv}}\\rho \\, dV + \\int_{\\text{cs}}\\rho \\, u\\, dA =0\n",
"\\end{equation}\n",
"$$\n",
"\n",
"where \"cv\" stands for control volume and \"cs\" stands for control surface. The first term represents the rate of change of mass in the control volume, and the second term is the rate of flow of mass, with velocity $u$, across the control surface."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since the control volume is very small, we can take, to leading order, $\\rho$ as a uniform quantity inside it, and the first term in equation (1) can be simplified to the time derivative of density multiplied by the volume of the tiny cylinder, $dAdx$:\n",
"\n",
"$$\n",
"\\frac{\\partial}{\\partial t}\\int_{\\text{cv}}\\rho \\, dV \\rightarrow \\frac{\\partial \\rho}{\\partial t} dA dx\n",
"$$\n",
"\n",
"Now, for the second term in equation (1), we have to do a little more work. The quantity inside the integral is now $\\rho u$ and, to leading order, we have to take into consideration that this quantity can change in the distance $dx$. Take $\\rho u$ to be the value in the center of the cylinder. Then the flow of mass on each side is illustrated in the figure below, where we use a Taylor expansion of the quantity $\\rho u$ around the center of the control volume (to first order)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![1Dfluxbalance](./figures/1Dfluxbalance.png)\n",
"#### Figure 2. Flux terms on the control surfaces."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Subtracting the negative flux on the left to the positive flux on the right, we arrive at the total flux of mass across the control surfaces, the second term in equation (1):\n",
"\n",
"$$\n",
"\\int_{\\text{cs}}\\rho \\, u\\, dA \\rightarrow \\frac{\\partial}{\\partial x}(\\rho u) dA dx\n",
"$$\n",
"\n",
"We can now put together the equation of conservation of mass for the tiny cylindrical control volume, which after diving by $dA dx$ is:\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial \\rho}{\\partial t} + \\frac{\\partial}{\\partial x}(\\rho u)=0\n",
"\\end{equation}\n",
"$$\n",
"\n",
"This is the 1D mass conservation equation in differential form. If we take $u$ to be a constant and take it out of the spatial derivative this equation looks the same as the first PDE we studied: the linear convection equation in [Lesson 1 of Module 2](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/02_spacetime/02_01_1DConvection.ipynb).\n",
"But in the form shown above, it is a typical _conservation law_. The term under the spatial derivative is called the _flux_, for reasons that should be clear from our discussion above: it represents amounts of the conserved quantity flowing across the boundary of the control volume."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### Dig deeper"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can follow the derivation of the full three-dimensional equation of conservation of mass for a flow on this screencast by Prof. Barba (duration 12:47)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"image/jpeg": "/9j/4AAQSkZJRgABAQAAAQABAAD/2wCEABALDA4MChAODQ4SERATGCgaGBYWGDEjJR0oOjM9PDkzODdASFxOQERXRTc4UG1RV19iZ2hnPk1xeXBkeFxlZ2MBERISGBUYLxoaL2NCOEJjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY2NjY//AABEIAWgB4AMBIgACEQEDEQH/xAAaAAEAAwEBAQAAAAAAAAAAAAAAAwQFAgEG/8QASRAAAgICAAMEBgcEBwcDBAMAAQIAAwQRBRIhEzFBURQVImFxgTJTkaGiwdEjUpKxQlRicpPh8AYkM0NVc4I1wvFFY4OyJTRk/8QAFQEBAQAAAAAAAAAAAAAAAAAAAAH/xAAXEQEBAQEAAAAAAAAAAAAAAAAAEQEh/9oADAMBAAIRAxEAPwD7H1Th/VfiMeqcP6r8Rl2IFH1Th/VfiMeqcP6r8Rl6IFH1Th/VfiMeqcP6r8Rl6IFH1Th/VfiMeqcP6r8Rl6IFH1Th/VfeY9U4f1f3mXogUfVOH9V95nvqnD+r+8y7ECl6pw/q/vMeqcP6r7zLuvfECj6pw/qvvMeqcP6r8Rl6IFL1Vhj/AJX4jHqrE+q+8y7ECl6qxPqvvM99V4n1X3mXIgUvVWH9V95nnqnD+q/EZeiBR9U4f1X4jHqnD+q/EZeiBR9U4f1X4jHqnD+q/EZeiBR9U4f1X4jHqnD+q/EZeiBR9U4f1X4jHqnD+q/EZeiBS9U4f1X3meeqcP6r7zL0QKPqnD+q+8x6pw/qvvMvRAo+qcP6r7zPfVOH9V95l2IFL1Th/VfeY9U4f1X3mXYgUvVOH9V95j1Th/VfeZdiBR9U4f1X4jPfVOH9V+Iy7ECl6pw/qvxGeeqsP6r8Rl6IFH1Vh/VfiM99VYn1X4jLsQKPqnD+q/EY9U4f1X4jL0QKXqnE+r/EY9U4f1X4jLsQKXqnD+q/EZ56pw/qvxGXogUvVOH9V+Ix6pw/qvxGXYgUvVOH9V+Ix6pw/qvxGXYgUvVOH9V95j1Th/VfeZdiBS9U4f1X4jPG4Vh8p/ZfiMvTxvomB7ERAREQEREBMoW8UyGtOK+KES1kHaqwPQ68JqTMx7vR+DZOWP8A713X4sYHdQ4wN9r6E390sPykh9Z66LiD/wA2/SZeLxjNbCwltVPSucDJ6dAOYL0952Pvk2Zxu/G4jfj+j1MqBSvNcFY7HgO8wLDLxvfsNgAe/nM85ePfv8P+x5pqSVBI0SJkZnFL6V6viUE2BAHfmOie/wB3TrAnrHGB/wAQ4Tf3Sw/KTG7MQDeIrnzS39RIxxWs8R9DWuzYJUuV0vNrfLv4Rh5OVZm2V3qgHIG5F69md9AT4k98Dr0rL/qDf4onvpOX/UD/AIolyU7uI10Wcj0ZHfra17H3QHpOX/UD/iiPScv+oH/FE74hc1GDbbX9ML7JPgT03IGybOG8Ovuzn7QUk8rAdXHh08/CB02ZlqN+rnPuFiyH1nm/9HyP41/WSpxjFaytH56jbrs+0XXP8Ptk2Vn4+LjrdZZtHPKnIOYsfIagVPWeb/0bI/jX9Y9Z5v8A0fI/jX9ZZr4niWVCztgoOthxog71oj4zPwc/Ju/2gzKewbsVCAlnHsdD119kCf1nm/8AR8j+Nf1nQ4jmH/6TcPjYsqZnEXy0zvQcxRXjIVdez6liDrR989wOJUJw/HxMh351x92uT9EgdQT5wLR4jmD/AOlXf4izuvOyrAebh7VAeNloAEpUcdpqrxlNVppuUdhYWDM48SfLXvknFeJYb4uXj9shNXKLevTRbqv2AwLjNxBhutMbR7t2MfylZvXezqzh4A8+fpLHCQVwVABCczdmG7+TmPL90zs/Fsu/2iCJ2fJZigkWglSVbyB79EQJ0s4y5IXI4Yx8gWMmrXjGv2luGP7qMfzlHIXD4fmYISumtw5e96xoKORgCfLZPSW6uOY7pWbFatnFjcrf0eUnv951AmCcU37V2L8q2/WdFOIFdrkUH41H9ZTrzxi8JrsQ2XvYpZEtYcygDqGPu6/yk3+z7OeEVLavLbWWRx5EEwOtcUHVr8QD+4f1nfJxP67G/wAI/rK3GK67c7h6Xf8AB5nNgJ6EBd9ZFVxS2mrGxkpfIyHp7RVZtMRvps/AHZ/WBoJXnn/iZFIH9mv9TI7E4qG/ZXYpX+1WwP8AOR3cewKK6WttINqghQNkfHy75JxLjGJwsJ6Qx2/UBRvp5wPCnF9dLsPf/bb9ZGV474WcP+YeS53GcLAsau+wiwJz8oG9+74zhOM8PyWantV0RYG2egC9Dv3dYHPLx79/h/2POq140D+0OCw/slh+Ukp4rhWDSXAa6AH+kPMeY98uI62IHRgysNgjuMDPyLeKY+NbcUxG7NC2g7eA35TQqbnqV9a5gDqQcSIHDconuFL/AP6mSY3TGqH9gfygSxEQEREBERAREQEREBERAREGAiIgIPUREBERAREQEREDi5uSl3/dUmZ64lmRwPGxhpeZK+02PDoW/OWeKAnhuSqkBmrYDZ11IkZ4rg0Fa3yE5tDu6wGRw2s87Y6KllttbufPlI/ISJuH3PnX3D0dSxBSxq+dhoAa93jLS8QxGUMMmrR/tToZmMe7Iq/jECYbCgE7PnM7JwceziOMvoVRQc1jv2Q1vWgCfPqfsl30rH+vq/jEelY/19X8YgQtgVektkIWW070d7CnWt689Svh4OXhWKq5C3o7l7ndNMen6/ZL3pWP/WKv4xHpWP8A1ir+MQJYkQyKGOhdWfgwnrZFCnRurB97AQPbqkvpeqwcyOCpHulVuGU241lF7WWrYQW5m69O7X2Sf0rH/rFX8Ynj5mKg22RVr+8IEGXwnEzDX2yEitORNHXKNg9D/wCIlDO4fZXXh4uPlGoHJJqYICahysdDzlp/9oOGI2jefkjH8pFZxzhNz1M1zlq25l/Zt36I8vfA6o/2dwKkUOrWvykMztsuT4n3/wApbp4dVRYGqJC8hRgepbZ3sn7ftkA47gMdCywn/tN+kl9a4v8A97/Bb9IHNPBcChqmrp0au7r3+W/PXhudpwrDRSvYKwNhsPMN9T3zxuK4wGyLgPPsW/SRjjWETrns3/2m/SBPZw7EsZCaEHIjIoA0AG7+nyni8LwVqFfolJGtbZASfifOeJxGiz6CXN/+Jp6vEK2JCVXsR3jsiP5wLFNYppSsEsFGgWOzPDRWchbyv7RVKA+4/wDxK751gHsYOS/yA/ORpxDJP0uGZC/NT+cCzdh496utlKEWa5/Z6tru3OTgYjM7NjVMXIZtqDsiR+nW/wBQyvsX9Z0uaSPbxMhPioP8jA9s4di2WVu1Q3WxZR4bPU7Hj1k9dSVlyg0Xbmb4yE5qd/Y36/7ZnA4njluVVuLeXYt+kDvLwcfMao3pz9kdr16fPzkfEiK61NSqt9hFK2aG1B7+vw2YPE8cNykXBvLsW/Sc5TU5dXZ20ZBAII0hBBgVeG8KxL+F4fbVBiqlgT3tsEbPn0MsJwTBF7XPW1tjLyE2sW6fAydMlK0VEx71VQAAK+4T30zyx8g/+H+cDivheHXUK+xDgb6v7R6jR6n3dJ36uwuXXolAH/bHWeemf/58j+D/ADnNnEaatdqlyb86if5QOU4Vj0VLXi7oAPtMv0iv7u+8CWqKa8ahKal5a0AVR5CV04jTZ1rS5h/2miziKJWXFGQwHfqowHF9+rblHe4Ca89kDX3y2BoaHQTB4jx3h9tVa12uSt1bH9m3cHBPhNyuxba1sQ7VhsdNQO4iICIiAiJ5A9iIgIiICIkGcrvg3pV9M1sF+OoEOPndqoudQtVtgSnzYeZ+M64hxGnAFZtV2LnoEG+nTZ+8SvlEvw7EtwqWuVHR1RSAeUfGVOKUvxFsBrMGznDP+zZtBegHMxHhA2jkUrclLWKLXBKpvqQJxVkCzLvp0AKuXR337E+exeDcQquttd+a5qXpWwtsg7BVvhoARl8EzRlMcfZrC1libCDY+zzH74G1nFqGXNrYlK+lqg9Cm+p+I75dBBAIOwZncNwxwzhL49zBkr5iWPdo9ZY4aHHDcYWb5+yXe+/ugWoiICIiAiIgc21pahSxQynvBEjrtrS/0VF1yoHGh0A3qTblHFdDxHNdmAbmSoAnyG//AHQLnKu98o38JxZVTyMXrQrrZ2omfdxmuviT0jkOPRWXvt39A+AlHiuZnNgW2v8AsUsrIpqrb2n2N7Y+Gh16QNReGcMyUW70Slw4DAlO8Ge+puGjuwaP4JSpHDuGWOqhv92oFrszluUdwHX4GXL+M4dKXE2gtUASviSRsD4wOjwjhx78Kn+ATz1Lwz+o0fwCd8OvutwarcsV12uNlVOwN90tMdKT5CBSHCOGBtDCoB9yz08H4c3fhUn4rMjHpvs5LqctxmXJ29trfQrQ9y6+75EzexLjk4lNxUqbEDEHw2IECcI4dWdphUqfMIJZrx6azuupF+CgSPMbLWsHDrqsffUWMQAJT7Xjn9Vw/wDFP6QNSeO6oNuwUb1szmo2GpTaqrZr2gp2AZk/7VPy8JAAJL3Vrpe/6XhA0MLOqzlsNXN7DaOx39Ng/YZZ1Pn0z/Q+HY9mOFa3Iygjqe9N79n4gACX+J8VHDha1igqtPPX16u29cv3r9sDRiZtN/ZnKyb+ZGqqBtTe16LvY8vEfKMTiyXU4jWoKrMhWbl5vohd9fh3fbA0Z70mPVmXcVes42QcWpq1tU8gZnB6Ede7X5y9ThV47vctlrWMNEu5I+zwgWtSJsilL1oa1Baw2qE9T/rRkfD7/ScCm3tUtLL1dPok9x185l8SspXJzmsZBeqVin94sNsAPHvgbmwASfCU6OKYl9LWi3s0XWzYOXoe49fAyNMj/d8y2ytjarFCqDZ0Po9PgZRu9GTB4ZfeorUtWtgsHcoRtbHxgbqMtiB0YMrDYIPQzrUycTiuBTQrKrY+PZay1sw0rdN8w/syX11jcmTaA5ox05mt10Y+Q8z+sDQJCjZOgPEzwMpdkBBZe8b7pgcezMfN4OxqZ3HahAEbW9a5j08ANy9VxWg+lWmvkFYDB/G1eoBHzBAgaNjpUhexlRB3knQE6GiOkwM7PGfwS/8AZjtVZQtY2RYdggDzB/WbGE6vh0lHDjkHXz6QJLbK6a2stdURRssx0BOlIZQykEEbBHjMzjeTj11U0XuurblBTvLAHfd8pxg8WxG4dWuGWsZSKUrPRt9w35DpvfkIGvGpjDj61bryaGFqWGuwV+0FP9EfE9Jewcq29ra8ikU2oQQoO9qe4/zHygW4iZDcdxXXIXm7LkDBGc65iOh++Bfws2jPo7bHYsnMV6jUnZlRSzEBR1JPhM/gl2JdhIcSxXPInacvgeUD8pxxrItXCvqrx3Ysmg51y9fvgacEhRskAeZmNg56YiGrKt7XItyXCqo6n2tb14AflNPMFrYzLRVVax17Nh0pECZWDb5SDo6Op7PmeEcLa9e2yqQiM7NuvIbp7R8Jp4GVXjYFPbOdWdo4YnY5dk7+wiBpxKd9vpLviY2V2OQqrYWChtAny9+pQ4Rm22caz8e11dVVVR17mZfpdPPqu4G2TobPcJyrK6BlIZSNgjxlPjhZeDZZQEnsj9HvmbmW5fDqRkileSrWlfJPL16fRA98D6CJ8/Vx29MfKN1PNkKOeqsAhWXR6g+IGu+aFXFsdxcWPKKWVCf3mIHQeffr4wLlNFdAcVjQdi5HvPfJJVwc5c0WfsbaWrOmS1dHulqAiIgRX0V5ChbV5lB3rfQyWIgIiICIiAiIgJnUcJqF2XdkhbmyHJ6j6K+AmjECivB+HqUIxU2gKr393595nVfDMSoMErPtJ2e2YtpfIblyIGavBMTX7YNe5+m7nq/lvXfqT0cNxcexHrqG0QoN9ehOz98txApVcJwabEsro01Z2vtMQD8N6lzWxo909iBl18FrTmTt7TQx9qroAw1oAnv0BNMDQAHcJ7EBERATN4jxLh2Nelea+rEIdNoTo+fSaUhycXHylC5FKWgdwYb1Axcnif8As/kNvIZXYj6tt/dPbeKf7P311La/OtJBTaOSpmiODcNU7XBoB8wkep+HA7GFT/DAq+tuE5PaqH5+2XkcchHMOvTr8TKxr4UEC9lbYy9Btxsr3cu992ppng3DT1ODRv8AuTz1Lwzexg0b/uCBBw8hs3tV4dbSHB/aFwVG+/oD46E1p4ihFCqAFHcBPYEIRcXHYUU9BtgidNk9Zmhz6Wco8Eu7dgAXJUnp85sRAynzcopYq8JvBcHrzL16a85njHua7HtPB7RdWVL2NYCGAGta3PpYgfPvjHTK/Cb7V/oKbVIrGwdD5gT2xXsud7OC5Lqx2UNq8u9a3rflN+IGDjj0ewtVwC5CV5T7anp8N+6Wr3ryKjVZwu91IAK8oHQHeu/zmpEDHaxlyBkJwS43BeQNzKND7YxbXputuTg19T29XIZTv75sRA+euW1mHZcEuXd3auS6kt5+M6xwcY0mrgN4akEK3aLvr3+M34gY65eQru44HcGsYMx516kDW+/yE69NyRcbRwa/tCvKW517vLvmtEDM9Z5v/R8j+Nf1lMGyqq9aeA3I1oOyXU9T8TN+IGVg5V9OLXW/C70ZFC+zy9dD4zjJy8rJr7JuE5AXmBPtr10d/lNiIGI1fMliNwZyLX53/ajZO/OW6b8lK1rTh7IqjQ3aJoRAy0vyaazTXwyzlJJ32i6Gzv8AOU7cS+zBxse/BtsFKcjCu0AONaIPun0EQMXEpupPOuBd2pBDWPcNtvXf9gnuLjvgUVivh7O1WyG7UczE95+c2YgY2XmZuRiXUeqMhe0QrsuvTY+MqV15SvWx4dd2+uV8m5g5C+5Qe+fSRA+exezr4hfy05WW7UKtjOmierdOutDXlLeJRw5ciuurFejISsqgZDtV+Pd4y/Xj8mbfkc2+1VF15cu/1k8CrXgpXiPQHfdgPPZv2mJ8d+csgaAHl0nsQEREBERAREQEREBERAStflCmxdFWRSBbo9U33H4Sw45kIBKkjvHhMumujGo4gCCyD2WLdS55evzJMDVle7Mrqy6sdgeazxHcvlv46Mkxwy49av8ASCgH4ynaalfNyL9ctXLrfhyjmH3mBoRMdM3KHDmzLrE3U3tJWu1YbHj85O3ELlotd6OVlI0O/lUjoTCVoxMz1jYbcOmtqrHt9qxh0AUeQk/plj5VuPTTtqtczMdDr3fn9kFXImeOJdwcCs1BmyAevKB06fHvHunOPxBMjHyLLSionQoCeYdO4/5Qq3Xl02JWwbQsYqm/Ej/4k8wuE5yY2MmNZWR2bOHbwQ9Tr7AZOOMWitrXwbBWH5OjDfXWunzhK1okGJeb6iXTkdWKsu96M6yrHqxLrKwC6IWAPiQIV7kWijHstI2EUtqeY72vX+2rFbg6IB2D7xM++/KzOGNZQKQhRg4ck92wdal3Aa9sRHyTWWYAjkBA1qBYkZvqVLH5wVr3z6PdqSTFSi26+x6/aoybiLNdwCkDfz0RCNiqxbaksX6LgMJ1MuvK7DgmM6kKzKqgnuHmfsBg8Qd8BXJFLg6uJGzWOvXXnA1IlHIsBx8W2u0tX2iEnf0wen8yDL0KREit2Lqf2nKCSCv73SBLERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBI3oqdlZkB5TzD4+ckiAlZ8GqzKF7M++m037LEdxI85ZiByaqygQovKCCBroNQK1FhsA05GiZ1ECnTwzGqD8ym5nbnZ7epJloIiuzhQGYAE+eu7+c6iBF6PV2tlhQFrFCv7wIbGoexXapC6DSnXdJYgQeh0ctq9mNWsWYeZPfOqseqmoVqvQHfXqd+fxksQPAqrvlAGzs+8z2IgcitAnZhQE1rlA6anqqFUKo0ANAT2ICc1olSBK1CqO4CdRAjOPSa1rNalFIIXXQanNmHjWljZQjFiC2x3kd0miBG+PS7q71qWT6JI7pJOLLUr12jqvMdDZ753ATxlVipYAlTse4zizIqqsrrdwGs2FB8dDZjGvryaFuqO0bej89QJIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiIDwmWubl3tYqejUMlnKRa+zrQPcPjNJ1V1KsNgzJs4TXdi5YWkLZZYeU60ddOm/fo/bCNQWKKO0axSoGy47p1W4trWxDtWGwZlY+BmCmpHdFqNvO9IA6De9bmvA5LorqhZQzdw31M6mMMPIt4n6aa/2lLkLs/SU9ND3AdfiTNmBy1iIwVnUE9wJ751MriyX4r+mYlfa2vpCpXm5feJzaMo41ovteuugMe13oueuvkBqCteQ15Ae69COUUkDZPfsA/nMqvOyr763rS42kAirl1XynxLefj908fBzLbMntiLVNqMa1HILAFHj/AK7oG5G5jXYPEcm1nN7UKwUhFb6BB7vs++X8KhscW1/8sNtNnfTQ39+4VagnXfE8ZQylWGwRoiAVlcbVgw8wdz2VOHY5xqGq7FKlVzyhD0I30MtwERECjldk2cEyCnZdg2wx13kfpJeHOz4VbPs9+i3eV30J+WpXzmVb1svx8fs1ICW2P3Hv8pNc+SiCztqVQkD6BPedDx98IrZ+NZm5n7IhHx0BRmHTmY9fuH3znhl3oeLVidmd1B2cnpyrzNo++alXaBB2pUt5qNCcW41Vt1drrt6wQPn3wIl4hjuyrWWdmOgFHuB/kRLUy8Xh7KlRrIoCWs2gvUjfd9mpqQKFvEmTJtqWkOKyBvtAvhvxlvHt7ehLeXl5hvW9z22iq9QttauB4MNzpVVFCqAFHcBCqeYLL8ivGquenpz2MnfrwH+vKS4Vr2VuthDNW5QsO5ovxBdYXFtlbFeQlDrYktVSU1LXWNKIHcSOrtee3teXl5v2ZH7uh3/PckgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgJS4ullmCVrrNvtDmQDZYbl2IEeNX2OLTV+4gX7BJIMQEREBERAREQERECnmhLMjGrsClAWsbm7ug1/wC6ULksbDd8a9qsZLlZDyc3s9Nke4HZmy9aWa51Da7tie6XWgBryhEOG4NZr9J9IdPpN08fhLE5VFT6KgfATqB5E9nkK8ns8ZlTRZgvh1Op7AREQEREBPAytvlIOjo68J7PnMcXcO/2gyrbrz6Nk6YJrYU9Bs+W4H0cREBESlZxTFrsatjZzKdHVbH8oF2JQGTacgWUiy2lyAUZCpXw2N+EvwE8gEHuIOunSewEREBERAREQEREBERAREQEREBERAREQERECtxKw1YVjrcKWUbDHWt+Uq46VX2Le2TkswQMRzFU/lLmdSb8O2tUV2ZSFDd24yaC+C9FRAJTlHhA6tyK6qRaSWVtcvL1Lb7tSStudA2iu/BhozNxeG3UHFU389VI3yH+idEdD5dfumnCK6ZtNmUcdNswBJOunTv6yxM/D4WMZ67Huax6xyodcul8Rrx85oQpM/FQ22pnNe4FmwtZPslT9Hp5+M0D1GpSownRqe1uFiULpFC66+Z9+oFuwOa2FbBXIPKSN6MylF/Wn0iw1G9q3s2OcHQ7vdvc061sFthdtoSOQeXTrK1WHYMktY6mpbDYgA6lj4n4bMIs0UrRUK1ZmA8WbZPzkkjrV1ewu/MrNtRr6I13SSFJDj5NeSrtXvSOUOxrqO+TSHGrtrFnaurlnLLpdco8B74E0hzMg41PMqdpYx5UTeuY+UmlbOK0p6Tyc71ghQe4b8f84EuPcMihLVBAcb0fCSSlwpXTE/auSzuzAHpoE9OkuwKudt2oxwxXtX9rXfygbP5SxXWtSBEGlHcPKV79esMXf7r6+6WoCV8nNpxWVbebZ6nS7CjzPkJYlDJoSzPrFoJSxdd/QsvUA/In7IFjNya8XFeyxwvsnXmekoLZlDExD2v7J0pTY6ly3RjLGTgHIyS7WDs2AUgrsgeIHx8Z36CRiJQLddm4atuXuAOwPyhEOBi1PUVvBuspsKlrDzdR3EfdNGRYtPo9XJzczEks2vpE95ksKTwsBrmIGzoT2ZeTbm81Iux6gpuABWzr4ny90DUicUtYyk2qFOzoA76TuAlXh45qmyG+leef5eA+ySZrcmFe3lW38p3SnZ01oO5VAgdyLJu7CrmCl3JCqoOtmSyDNpa6j9n/AMRCHT4jw/KB5hZLZNTF6+zdGKsobmG/j857nGxcZrKieav29D+kB3j7J7iGo4yNSgRCN8utakxG+hgUe37Xi1NaN7C0NYdeOyAPzkzZNVnpFNdwFlQ9vXXl6SDh/D3xcrItd+YMFSsa+ig7h95lo41W7WVFV7RpmA6mEVuHJRU9ldFoZSqsV37QPiT8ekuWWJVW1jnSr1JlTEwBijHKvt61Ku37+9dT8xLGVSb6goYKQwYEjY6HcCtw/L9IyMpSGAVlZQ/QgED9JBi5xbP5ms5qcgla1A+iVPT7Rs/KWF4fyq5FxF1isHsC/SJ8de6TNiVtjLQvsqgHKR3rrugTyCzNxqbOztyKkf8AdZgDJ5Qu4RhX2NZbWWZu8ljCrtdtdoJrdX15Hcq5mRki9KcRFZ+jOW7lXf8A8/ZJcTEow6ymPWEUnZ985ysV7mFlF5otA5eYLzAjyIhDHvsN749/L2iKGDL3MD/LqDLJkGNjdiWd7GttYAM7eXl7pPCkREBERAREQEREBERAREQEREBERAREQEREBERAREQEREBKXEq2asOE7StOticxBK+7XjLsq52PfkCvsL+y5W2djfNAylPDrMtDi4r5GlLDRbewR4k++blLtZUrPWa2I6qTvUiOIB6PyOymkcoPmvkfsEndgiM7b0o2dQK+SP8AfMRvHmYfap/SWZTe1L8rEFTBwOazYPhrX5yPmut7bLrZmNZZKqd6U66HfzgaEqrhH0kWvkW2KrFlrbWlP+tyGjKyMmygKa61atbW6b5v3gPLX5zQgIiICIiAnFlSWlCw2a25l+OiPzncQERECvxHrw7I/wC238pYXuEzsrOR6MqhhyXHdaVk9X2OhElfIvGR2OPStgq5e0JbXf5fLrCLkq5OLZdZzJl3UjWuVNairKe3Pto7NQlf9Lm9o/LylqFcU1LTUtab0PMzuIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAnLtyIzaJ5RvQ7zOpXy6HtUNUwW1N8u+4+YPugZbW8SuyKGXINFd5AC9mCF6Ejr4np1mviWm/GrsYaYj2h5HxmN2mTSi05Diims6RhWXcdPAjoJcr4tw+lK6UZ1H0VHZt+kI0olZM1LAezS1teHIR/ORtm5AOl4defeWUfnCrsSh6fkc3L6C3NrfL2q7nQ4gyIz5WLbQijZY6YfdAs10VVEmqpELd/KutzxKK67XtQaZ/paPQ/KV1tzr1FtSUIjdVDsd6+QnVl+VUB/ugtPiUf9YEtOLTRY9ldYVn79SaURnZO+vDbh8GU/nOjm3gdMC/7V/WBciUTnZPhw28/wDkv6zwZ2V48OsA99iwL8TMqzrM8IuOhpJIfb76qCN66TTgIiIFfOzK8KkWWAkE8oC98z8LiWTavb3Iq1llBX90N3Hf2b+MvcRx7cilfR3VLq25kLDY3rXX5GY6YeZj2qjY9t1Q0TqwHZ3v3dCYTX0BrRnV2RSy9zEdRIbcKm63tG5w2gDysRsDznD2Z1lW6aaqn8rW3/Ke1pmcn7e6oHzRP1gc4+DZXldvdkG4gFU9gAgHzPjLsqCi2wBkzX5T1BVV0funjrnoy9lZTYg+kHBBPwIhVuJSxrsjL7RyDjqrcoQqCdjv3PGXios9izEZP7SsDAvRK6NmAe2lJP8AZc/pOubI+rq/jP6QJolW3JsoCmxK/aIUe33n7J2Lcg/8lR8X/wAoE8So5z+9RjKv9pj+k43xQgFDhkHuPM36QL8jqtW3n5T1RirDyMpBeME9XwlHuDGWceh63eyxwz2a5uUaHSBPETiu6u13VGBNZ0wHgYHcREBERAREQEREBERAREQEREBERAREQEREBERAREQEgzrmx8DIuTXNXWzLvzA6SeRZVAycd6S3KG8YGVhZFOJi25zkhH0oB+lYw7yR57/lPMniiW14D2VuC9qv7K8w8QeomhdwzDvt7SyhWO9+7fnrznnq2hQez5kYkEHe+XTc2gPCE654fmWZllzdny0KSqEjRJB0ZNdm00XrS5bmbyG/P9DOlxaUyDeqadu8g989bGqbIW9kBsUaBgUeG2pkZeVZ2L83a6DsmtLyjpL+RemPXzPs7PKFA2SfKK6hW9rD/mNzH3dAPynGVjekchWxq2Q7DL98KoPxEYuTlF9lUrTVQOyGO+n2Dcs38Q7PJRFrDJzhHctrTHuAHjPU4ZjJYrhCSqlep3zb7yfMw3C8RiSUO9DrzHoRrr8eg6wiRs6hctcYse1Y6A17tyPH4jVkZjY6K3RSwY9x0dGdXYlHbduzcljMg5vPXcPnsytZhUVX42PWhWthZzaPU9POBbws2nOV2oLFVbl5iOh+EjzclGx7qhzc781QCjrvlkmLi04ukqYgAcqpzdAPhOvRKfSfSOX9p8enx15wK2DkZV1dX+7LXVyj2i+yfkJfkePSKEZQdqWLAa7t9dfzkkKREQEREDD4hxC7Huya3YrWLa+zb+Esv3/zmnk5VKtZjs+n7IufcvdOGwK8hblykWwWWc4HloAD7hDcOrbLa9mYq2iyHuJHd/8AEIr8Bt//AIzGq7OwEJvZXp3900LMiqp+WxuU8pfZ7tDW+vzntVSUghBoFi3zPUyPKxK8p6Ws/wCU/MPI9O4/68IFDN4m12NlJhIWetSC5YKVOu/XfLGNm/7hj3ZPR7WVCFHcxOv5ztuG4rLo16PNzFgfaPxPjOLOGrZbs2utXN2nZr00/fvcC9I8i+rGqNtzhEBA2Z7VX2aaLM5J2SxnGZjJl0Gp9gHqCPAwqK0q3EsYb3+zcgfwy3KWXw9MlsfY6VnTHeiV1rW/jqWK3qXVKWAlBrl5tmBW4xbVXhctzey7qCviw5hsa8ekrYGX6Fg6vDlO27KlQvtAf0QZezcZb6mYVhrQByk+HXchyOF15Nxeyx+TmD8g6AMPGESNmFqb+VezurGgr+etiWaX7SlH/eUGVsjhuNk2my1WJPkxHXWt/ZLSqEUKo0ANAQIKs1LMy3GIKuncT3N0Hd9sp42dXTZfUabmdLWBK176E7H85dGHSMo5IB7Q+/p5d3yE7SlUussHfZrY94GoCi4X186q6jycaMqY/F8fIyuwTm2V5gxHQ67/ALJflVOH4tZrK1AGvm5T/e74FdeKdtQzVKFdLFDKeulY9D9nX3S2LLGzmrBArSsMenUkk/pIvVuMMYUIrIn9lu/4+c7ox7EyGuttDkoE0F13Enf3wLMREKREQEREBERAREQEREBERAREQEREBERAREQKd/E8fHssS3mBrIB6b6aB38Os89Z0Hs2Uk1MoYuemtnQHxM99XUnKyL3HMbuXYPgBrp906ThuKgA7Pm0djm8PL7IRDZm3LdawFRprYIVJPNs66nyHWaAPSVXwKmNzKzK1w053vp5DfdLFaCqta1JIUaG++FVeKXLTjDnZ1V3C81Z9oHw1JcTKqyaga32QPaB7wffOrsdLzWbNns35wPfOVxakvW1F5CAwIXoDs76/P+cCZmCqWPcBsynhrku6ZNli8ltezXr6PiNfI9ZdI2NHqDK2JRdj/s2tD0qNICvtAe8+MCTKqF+M6eJHsnfcfA/bMevGvRbLrcq2yyqxVOzr2SBzDp7z9wmzRSuPUK0JKjff8ZC2Er5BsZ25Swcp4FgNb/lCFfDsSrIGQlX7UdzFif5mWoiFIiICIiAiIgV8q1wyU0kC2zrs/wBFR3n/AF5yJsxnw63xyrWWN2as/wBHfXqdfAyXJwkyLA7PYp5eUhTrY2Do/ZC4WOqsgrHIxB5PAEeQ8IFZuItTfViOouymUs3ZnSj7ZyOKm/h6ZNCKhazkPanovvOvl9suJiY9bq6VKrDZ6e/v/lPPQcbshWKVCgkgDprffCPbbit1CrorYSCfkT0k8rjDpD1sFIFYARdnlGu7pLEKhzLmx8S21QCyqSNznGqNQLvkvdzddsRofCWD1Gj3SsmBQnMqhgjLyGsMeUD3CBZ2CNg9POZqVVC8UWKrI7m2iwHrveyN/M/KX6q0pqSpBpFHKB5CcV4tFThkqVSO4gd0CVnVF5nYKo8SdCe7BGx1E4uprvqaq5FdG71YdDOlAVQqjQA0BA9iIgIiICIiAiIgIiICIiAiIgIiICIkGbe2LiPcqc5TR5fPrAniQJlpZUjqrEM3KenVT75Fbk3jiK0V1q1YQM/X2upI6fZAuRKleTe+a6CsGhW5ObfUHW9/DrqW4CIkdtorNYI3zty/CBJE9kTXBclKdHmdSwPh01+sCSJy7BEZjvQGzoblM59jANRh3Wof6Wwv84F6I7hsxAREQEREBERAREQEREBERAREQEREBERAbiIgIjwkK5WO9YsW1ChblB34+UCaIiAiJHTdVkJz0uHXZGx5wJIiICIkZuAyFpIPMylgfgR+sCSIiAiVcjiGNjWMlrMCoBOkJ1vu7vhLKsHUMvUEbED2JXvy66bezYPzcvN0XfSSUXJkVC2o7Q9xI1AkiVs662mpDSELM4XbkgDf+eh85LR23J+3FYb+wSR98CSIiAiIgIiICIiAlTiv/p13y/mJbkWS4rodjU1wA+go2T8oFJXL33WYx3W1lY2vcT/SP2aHykqoa+KvZZZvnq6eAAB/znFXEQFAHD8useXYn8p62enNzHCyiQNb7EwinbxBlterGYLzZA57nHsKpG++Sesxa+Tegbkxa20NdHJ7j933yT0urT74dlMH6sDSSDOkz00R6vylBGiOw7xA4XPenhovtsDMhDOCOvIT36HdIr8+1sQM6ft670PIvgDogE+HQ63LC5NC8yrw68BhptY/f8Z02XV7X+45B5xpv2B6/GBXyMm+pnuZuxtQqhTfMhB8fj3/AGSTDvtvuxjkVmu4KyuNaG+h6TtuIqT/AOn5ZP8A2TOfWRd1PqzMJXuJrA1AvXra1RFLqj+BYbEw0wsxcOm2i6/t+0bac2k+k3h5E6mwuRdYu1xnU+TsBIns4j/QxqPnaf0geWU238QNd3XFVQ4AH0j3aPnrW/nKpuz8nNbGVhXjh3RrFXroaPQ+B66+2Ww/EietWMB/eJhjxBf+HXinfU9SIEuBa1uIjOdsCyk+eiRv7pYnNYIQbAB8QvdudQpKmTxCrHt7Nq7nbv8AYrJH2y3ECLHvGQpIrsQD99dSWIgJSzXdcqjVrVoAW3rak7HQ/Lepdg9YGNTm5uXnZC0OiV1qGrR06t578tzluPbw3t5Urt59Ilh1tfE/cRNoKoYsFAY9511MjbGx3GnoqYa1ooD074SKtWdZkOorr7Ou1CarH7ydfu/67pSGZa/Eq677EarGs09qjlHOQdD/AF5ibXZp7PsL7H0end8Jz2NPKyGpOVjthyjRPmYIz/WLXY+VyMisLOypIO9++aNVa01LWm9KNde8zwY1CsGWmsMvcQo2JJCnhMmrPyrMwWNX2eHzchDLo7Pcf5fbNaCARojfxgJBmC04r9iTz68O/Xjr36k85dBZWyN9FgQYGTTxJa8M3Y1VluJX9O2x/a9/TvMnw+KrkNcl1ZosqPVWPU9N/pLVOHTRU1artWOzzddyRqKmbmapCd72VHfCMuziiX8PuBBVnrblKg8vd0G/OUMyynKxeap1XnuQhfFVA5SzDw1+U+k5EChQqhR3ADoJyMegFyKawX+kQo9r4+cCSBOa0WutUXfKo0Nnc6hVNswvxFsOtNhay1lm/oHwEqcLyBVwe3LKjlG25Qf3VA/KXbcGu213LOFs0XUHQYjunlnD6m5lUtXW/R0TorfpCOb+IpVj4+Qf+HaN+/6OwPjI04qmRhZF+Op3UOnN49O/4b39knXh2Oq1KQzLUQa1Zthdd0kXEoRlZK1UqCo10GjAy6uNHkRCDc7KTz1jvO+g15666nO82ziFdoyQKutSM1XeT1PT5TZSmqsAJUi8vUcqganRVW1sA6Oxvwgj0A6GzsxEQrKybL34lU2KVRVbsndhvn8da92vvkz5dw4hVi+xyk7Zx4dCQvxOpOuHQuScgJ+0Pv6fHXnFuFRarKyn2m5yQdHfnuEVuJ3EWoKlBsp1cWJ+iu9fPfUTv1ki3BGUnmt7NdDuHQbPzOpNZhUWtWXUk19B17/j5yrmrV+3qqoYX2lRvl+l1HXfkIHPF8hij41bKrFOYk958gvv3Os3KyUzKcehwrPrYNRI9533d00GqrZ1dq1Lr3MR1EgycGrKtWy1rNqpUcrle/4QLMTmtBVWta70oAGzszqFIiICIiAiIgIiICIiAiIgexPJUyeI4+NatTFnsJA5UG9E92/KEXJ5OKLlvqWxN6PmNEe6dwpE4tsWmprHOlUbPTcgw8i2yyxL1CNoOoHgp8/f0gWp4zBRtiAPMz2U8mtMnProuXmrWtnKnuJ2APzgXO+JxTUtFQrTfKvds7ncBE5tsWmp7HOlRSxPuE4xshMqkW175T00w0QfeIEsSkufvjD4Whpag+/He+v5SBeKX2mrsMQMlpcKzWaBAP5jrA1IkdRtI3aqr7lO51bYKqmsbfKo2dDZgdRKeJxTHzXC0C0j94oQB85cgIiZ7ZV9F722kei9p2YBGivhze/rA0IkORl045UWtonyG9DzPkJNAREgzb2xsZrUr5yCBregPefdAniQ0i8MxusRlP0Qq61GVkDGqFhUsNgHXh7z7oE09kFNzWdDUy7G97BB+Yk0AYmRRjvZmZNduXeHoYNXptDlPUbHj12PlNKnJovYrVdXYwHUKwOoEsREBBIA2ToCJX4h/wD0bvLlO/hA8x8izIIcVctB+ixPVh4HXlLM8GgoA7tT2AiQ5OVViqrWk+0dAKCSfkJ0cikUi42qKz/SJ6QJInFdqWjaHY3retQbaxz7dRyfS693xgdxPFdWTnUgrrex1lN817DX6LXzluYkOCpIGu7fxgXYgd0QETmyxKwpdgvMeUb851AREQEREBERAREQEREBERAREQB7pgV14mNepsY+kpcxYNss466IHjvYm/PCq83NyjfnrrCIMFHFb2OpQ2uX5D/R/wA+n3yxEQrw611185mYYyH4k197sqhCoQgAd/TXn3ffJ+J4gyq1bkrdk66sJ0R490r4b4jZdKY9dQVqTZ06kHY/zhGrKtp5eJY/T6Vbgn+Ey1KuVsZWJy/WHfw5T/lCrUSKvKpstsrVwWrIDe7c5GVU+SaEJZ1+lodF+JgdX9m6mixgDapABPePH+cpYva1cGY1lnu222C7JbZBOvlLuRjU5Sclycyg7HUjXzE7qrSmta61Coo0AIGVhYoq4kL1xXWt07IM30uhJLN8ZJi02h6cc0lFxnLdoe5u/WvtmnEJCcW013cvaLzBTsCdxCuURa1CooVR4AaE6iICZXEAbbb67K7Hbk5KFCEjbDRbf+tfOarAlSAdE+PlM3kYXPU/Eb1ZRze0qgEeY6QIs/FtpLZByWNJRUtrAG+Qd+j/AK75r613T52urHzM0rbnXvkV2FVIHTQ7u4am/UnZ1qhZnIH0m7zCY7lDiTMStRW3sSrFuyXZc+C/fv5S/EKx7jk1Y1OPWzl1q3agGyF+Pn5S3lXD0BWoPsWaUPy75QfHUu6G966xoa1rpCMmntsPEqNdN3Z9rvsgOZgmj0Pl1+zclpyMlMPJ7cBsiok6Xr0PUfz18pdyGsSktSgdx15SdbkWGlha2+5ORrSNJvfKAPH398DPw0Zcxr0rsKPUVNzj6bDrvXl16SLhRtGZ2xpvTtv6K0BEUeHWb3wiCEREKSPJQWY1qMNhlI+6STi+uu6h67RuthpuuukDjEc24dFjd7VqT9kmlPhaCvBQKNVklq1/dUnoPs1LkCnxOmx6FtocrdSedTrfTWj9xM49BryMOivtS1OuZug/ab8fd37lrIrssUdlb2ZHf03sTqmoU011LshFCjfugc4+OmNUK6+blH7xJmDxGrLqPEUx8e17Mllbn1tAgA2Pj39J9HKVmRZ6VapsWuqkAkAbZvHfw8ITWImPmpmm7Dd/RCzmpF3rmKE93lvpNa/0o5eIwKIpUpzEbOyu+75Tv0l8l6q6D2KuH5iR7QKkDQ8PGWceg1Jp7XuO9hn1sQMuv1qeYO7kBnI2mt8p9kfA/lLWHl519w7fE7KokjZPXu2P0+M0IgV8+prsUooJ2y71362NyDDszmy3TIrUVKCQ3i3tHX3CX4hSIiAiIgIiICIiB5PYiAiIgIiICIiAiIgUeI1X2mtUrayrTF1DAcx6aB93fIcfH9WWJyY5YPWEPZD+kD/nNSIHjDmUqe4jXSV6cJarBYbbbHVORS53yiWYgZHYPfhJjBLBdiuHK9y2EHff7++X8Kt0qbtE5Czs3LvetmWIgIiICIiAiIgJ5PYgJlXrlZGagavSpZ7I5Rrl31JPvHhNWIFPkyMfKbsahZVawLe0ByHuP8pciICIiAiIgIiICIiAiIgJX4gdYN/Uj2D1HhLEEBgQQCD4GB4mgihe7Q1PZWqwhSy9nbYK1O+z5un+vdLMBERATMx8XIPE2yLqwCCw7QtvanuUDwHj9s04gUacJ683tNr2Sl2Xz22t/wAvvl6IgIiICIiAiIgIiICIiAiIgIiVMux2uXGqJUspZ3/dX3e+BbiYGPklK3sW8M9LAKFbfaDpsHzPXv8AObObecXEtvC83Zrza84E0Sjwuy61LbbnJDEaHgDrqB7vD5SHifEajw256SzgA6dBsAg+MJWpEjrvrtd0Rtsmtj4jYnVti01Pa50iKWY+4Qri/JqoKq5PMQSAqknQnVNqX0rbUdow2DKFuXSvEMe4FnWypkHIpPeQR/IzTgeRKl/EqMa1q7SQV1vQ30Pj8J0ufQ+M96liity/R6k+6EWYlZs/GU07s6XAlW8PtneNlVZSM1ROgddRqFTEgDZOgPGcpYlgJR1YDptTuVuLELwzILa0EO9ytw8VLxDKehOypZEOtaDHZ9qEakTxWVxtSCPMSOq7tLbqyvKa2A+II3uFSxIr8hMcIbA2nYKNLvRMlgIiICJHY2nrUOFYnej/AEgO+dCxDYawwLqNld9QIHURObGCIxJA6EwOomfwnNGVUE3zmutOZ972xHWX1YMNqQR3dIHsRK+ZmJiVsWBLcu1GuhPgNwLER3Dr0nLutaM7sFVRsk+EDqJQo4kt+E19fKT2jIgb2eYg6A6+MmGbWWrUh1Z1LEEfRHv8oFmJXrzse2hbls9hjyqSNbPulPB4hk3s/NQnYi5k7U2AePTpA1IkJyUW81uyA9NdeuzJVZXG1II7ukD2IiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgJXysXt9MlhrsUEBtb6HvBHiJYiBQo4WiWLZeUtdD7GqwgX5CWcyn0nEtpH9NeWTRAzn4YaqWqwrmoVxyuDs6HmPIzyjhFdCdiLObF3zGorvZ14ny8ZpRCRQxuHPjvzjKcnoD7IG1HcD+su3VrdS9TfRdSp15GdRCqleD2NONXU51S/MSe9tg7/nLcRAq24FV1mS77JvrFZ9w6932yGzhK2WMTfYE5gyovQKdaPxmhEJGW/BlbEahb2H0QpZd8qg7C6l3Hoap7HdwzProBoAASeIVVzsOrJrdnTnfkKqD1Hd5TPzMPLylrRAFWyhO05unVTvR+O5tRCRj18KyAqILuyUKrHk/fUnXy7vjqX8Wu4PZbkcgsfQ0m9ACWYhUdqWM9RSzlVW24I3zDXdJIiAiIgVM2nIeyuzGKB1DIebwB11Hw1IuH8NOC4Pam1ihFjsPaY76fdNCICVuIUi7CuHZh3Fbcg1s71LMQMxeFt2y212mgLUqoqDXK3XZI8Za4fiehY3Yc5fTE8x7zvr1lmICU87FbNYUt0pClifNu4fZ3/ZLkQIFo7UVWZCL2qDwOwDI+I0XZFda0lOj7YP3f6B6/KW4gZR4P+0QC39kjLYF1/TGuvzAP2y1dgrb6T+0ZTkAKSB3AeX3y3EJGf6r/Y0V+kuBTsDQA2p6a+PvkN3CrGvPZWKlKlXrTXQMNAk/Z981ohWUvB3XGuobJZxdathcjTdCCevyl7FxlxVdUPss3MAB3dJPEBERAREQEREBERAREQEREBERAREQETyewEREBERAREQEREBERAREQEREBERAREQERG4CI2JRtsszLWpocpSnSywd5P7o/WB02eGvejHqe6xBs6IA+0z0ZGYQd4LD/wDIv6yarHppINahdLy9PKS7HmIFVMjJJ9vCdR7nU/nPLcnJ5QaMJ3/vOqy1seYjY8xAoek8T1/6cn+OJ1XdxGzo2JXT/aazm18hLxYeYmdkZxdHNNgppT6WQw2P/EeMCZsxcdqqLrBZkOQNIPPx9wluZ1dFdedSqDelaxmJ2xJ6Df2maOx5iAiNjzEbHmICI2PMTzY8xA9iebHmI2PMQPYjY84gIiICIiA8IiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgIiICIiAiIgZ/HMyzCwOero7tyg+X+tT5J8m92LNc5J7/aMRAibIuZuRbrB+8eYzpbbEUKtjgDwDGIge9vb9a/8Rjt7frX/AIjEQHb2/Wv/ABGO3tH/ADX/AIjEQInybXXrZZyt0A5j7U8DuSqu5PL1VN9FiIEnb3dtvtX6rr6R8512931r/wARnsQPO3t+tf8AiMdvb9a/8RiIDt7frX/iMdvb9a/8RiIDt7frX/iMdvb9a/8AEYiB3Xl5FTh0ucMPHmn2PC8pszAruce0eh+IiIFuIiAiIgIiICIiAiIgIiICIiAiIgIiICIiB//Z\n",
"text/html": [
"\n",
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.display import YouTubeVideo\n",
"YouTubeVideo('35unQgSaT88')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### General conservation laws"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"All conservation laws express the same idea: the variation of a conserved quantity inside a control volume is due to the total flux of that quantity crossing the boundary surface (plus possibly the effect of any sources inside the volume, but let's ignore those for now).\n",
"\n",
"The _flux_ is a fundamental concept in conservation laws: it represents the amount of the quantity that crosses a surface per unit time. Our discussion above was limited to flow in one dimension, but in general the flux has any direction and is a vector quantity. Think about this: if the direction of flow is parallel to the surface, then no quantity comes in or out. We really only care about the component of flux perpendicular to the surface. Mathematically, for a vector flux $\\vec{F}$, the amount of the conserved quantity crossing a small surface element is:\n",
"\n",
"$$\n",
"\\vec{F}\\cdot d\\vec{A}\n",
"$$\n",
"\n",
"where $d\\vec{A}$ points in the direction of the outward normal to the surface. A general conservation law for a quantity $e$ is thus (still ignoring possible sources):\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial}{\\partial t}\\int_{\\text{cv}}e \\, dV + \\oint_{\\text{cs}}\\vec{F}\\cdot d\\vec{A} =0\n",
"\\end{equation}\n",
"$$\n",
"\n",
"To obtain a differential form of this conservation equation, we can apply the theorem of Gauss to the second integral, which brings the gradient of $\\vec{F}$ into play. One way to recognize a conservation law in differential form is that the _fluxes appear only under the gradient operator_.\n",
"\n",
"Recall the non-linear convection equation from [Lesson 1 of Module 2](https://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/02_spacetime/02_01_1DConvection.ipynb). It was:\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial u}{\\partial t} + u \\frac{\\partial u}{\\partial x} = 0\n",
"\\end{equation}\n",
"$$\n",
"\n",
"If we look closely at the spatial derivative, we can rewrite this equation as\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial u}{\\partial t} + \\frac{\\partial}{\\partial x} \\left(\\frac{u^2}{2} \\right) = 0\n",
"\\end{equation}\n",
"$$\n",
"\n",
"which is the *conservation form* of the non-linear convection equation, with flux $F=\\frac{u^2}{2}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Traffic flow model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We've all experienced it: as rush hour approaches certain roads in or out of city centers start getting full of cars, and the speed of travel can reduce to a crawl. Sometimes, the cars stop altogether. If you're a driver, you know that the more cars on the road, the slower your trip will flow.\n",
"\n",
"Traffic flow models seek to describe these everyday experiences with mathematics, to help engineers design better road systems.\n",
"\n",
"Let's review the [Lighthill-Whitham-Richards](http://en.wikipedia.org/wiki/Macroscopic_traffic_flow_model) traffic model that was offered as an exercise at the end of Module 2. This model considers cars with a continuous *traffic density* (average number of cars per unit length of road) rather than keeping track of them individually. If $\\rho(x)=0$, there are no cars at that point $x$ of the road. If $\\rho(x) = \\rho_{\\rm max}$, traffic is literally bumper to bumper.\n",
"\n",
"If the number of cars on a bounded stretch of road changes, it means that cars are entering or leaving the road somehow. _Traffic density obeys a conservation law_ (!) where the flux is the number of cars leaving the road per unit time. It is given by $F=\\rho u$—as with mass conservation, flux equals density times velocity. But don't forget your experience on the road: the speed of travel depends on the car density. Here, $u$ refers not to the speed of each individual car, but to the _traffic speed_ at a given point of the road. \n",
"\n",
"You know from experience that with more cars on the road, the speed of travel decreases. It is also true that if you are traveling at fast speed, you are advised to leave a larger gap with cars ahead. These two considerations lead us to propose a monotonically decreasing $u=u(\\rho)$ function. As a first approximation, we may consider the linear function:\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u(\\rho) = u_{\\rm max} \\left(1-\\frac{\\rho}{\\rho_{\\rm max}}\\right)\n",
"\\end{equation}\n",
"$$\n",
"\n",
"![velocityvsdensity](./figures/velocityvsdensity.png)\n",
"#### Figure 3. Traffic speed vs. traffic density.\n",
"\n",
"The linear model of the behavior of drivers satisfies these experimental observations: \n",
"1. All drivers will approach a maximum velocity $u_{max}$ when the road is empty.\n",
"2. If the road is completely jampacked ($\\rho \\rightarrow \\rho_{max}$), velocity goes to zero. \n",
"\n",
"That seems like a reasonable approximation of reality! \n",
"\n",
"Applying a conservation law to the vehicle traffic, the traffic density will obey the following transport equation:\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial \\rho}{\\partial t} + \\frac{\\partial F}{\\partial x} = 0\n",
"\\end{equation}\n",
"$$\n",
"\n",
"where $F$ is the *traffic flux*, which in the linear traffic-speed model is given by: \n",
"\n",
"$$\n",
"\\begin{equation}\n",
"F = \\rho u_{\\rm max} \\left(1-\\frac{\\rho}{\\rho_{\\rm max}}\\right)\n",
"\\end{equation}\n",
"$$\n",
"\n",
"We can now use our numerical kung-fu to solve some interesting traffic situations, and check if our simple model gives realistic results!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Green light!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's say that we are examining a road of length $4$ where the speed limit is $u_{\\rm max}=1$, fitting $10$ cars per unit length $(\\rho_{\\rm max}=10)$. Now, imagine we have an intersection with a red light at $x=2$. At the stoplight, traffic is bumper-to-bumper, and the traffic density decreases linearly to zero as we approach the beginning of our road. Ahead of the stoplight, the road is clear.\n",
"\n",
"Mathematically, we can represent this situation with the following initial condition:\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\rho(x,0) = \\left\\{\n",
"\\begin{array}{cc}\n",
"\\rho_{\\rm max}\\frac{x}{2} & 0 \\leq x < 2 \\\\\n",
"0 & 2 \\leq x \\leq 4 \\\\\n",
"\\end{array}\n",
"\\right.\n",
"\\end{equation}\n",
"$$\n",
"\n",
"Let's see what a plot of that looks like."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy\n",
"from matplotlib import pyplot\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Set the font family and size to use for Matplotlib figures.\n",
"pyplot.rcParams['font.family'] = 'serif'\n",
"pyplot.rcParams['font.size'] = 16"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def rho_green_light(x, rho_light):\n",
" \"\"\"\n",
" Computes the \"green light\" initial condition.\n",
" It consists of a shock with a linear distribution behind it.\n",
" \n",
" Parameters\n",
" ----------\n",
" x : numpy.ndaray\n",
" Locations on the road as a 1D array of floats.\n",
" rho_light : float\n",
" Car density at the stoplight.\n",
" \n",
" Returns\n",
" -------\n",
" rho : numpy.ndarray\n",
" The initial car density along the road as a 1D array of floats.\n",
" \"\"\"\n",
" rho = numpy.zeros_like(x)\n",
" mask = numpy.where(x < 2.0)\n",
" rho[mask] = rho_light * x[mask] / 2.0\n",
" return rho"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Set parameters.\n",
"nx = 81 # number of locations on the road\n",
"L = 4.0 # length of the road\n",
"dx = L / (nx - 1) # distance between two consecutive locations\n",
"nt = 30 # number of time step to compute\n",
"u_max = 1.0 # maximum speed allowed on the road\n",
"rho_max = 10.0 # maximum car density allowed on the road\n",
"rho_light = 10.0 # car density at the stoplight\n",
"\n",
"# Discretize the road.\n",
"x = numpy.linspace(0.0, L, num=nx)\n",
"\n",
"# Compute the initial traffic density.\n",
"rho0 = rho_green_light(x, rho_light)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEYCAYAAAB2qXBEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4VdW9xvHvLzNhFAgggoIiCMggEGxVbLRqrVWr4lCnMqj0Wqu9tlqtdazaOvV6a1utVCE413m4zrXGuZIAMskgMs+TQBLIcHLW/eOcBESQhHOSvc/K+3mePHncyUneZxvyZq+91trmnENERGRvpQUdQEREUpuKREREEqIiERGRhKhIREQkISoSERFJiIpEREQSoiIREZGEqEhERCQhKhIREUlIRtABGlu7du1cr169go7hhfLyclq2bBl0DG/ofCaXzmfyTJkyZb1zLq++n+99kXTu3JmSkpKgY3ihqKiIgoKCoGN4Q+czuXQ+k8fMljTk8zW0JSIiCVGRiIhIQlQkIiKSEBWJiIgkREUiIiIJUZGIiEhCVCQiIpIQFYmIiCRERSIiIglRkYiISEJUJCIikpBAi8TM9jWzN8zMBZlDRET2XmBFYmanA58AB+3h8zLN7FYzm2tms8zsYzM7qmlSiojIngR5RXItcDzw0R4+7y/AOcAI59yhwATgbTMb3Mj5RESkHoIskiOdc1982yeYWR9gHHCHc24dgHPuIWAhcHvjRxQRkT0JrEicc5F6fNrpgAHv7nT838AJZtYq6cFERKRBwj5rayAQBZbudHwRsYdy9WvyRCIi8jVhf0JiR2Crc65mp+Nb4u877OpFZjaO2JAYeXl5FBUVNVrA5qSsrEznMol0PpNL5zM4YS+S3bFv+6BzbjwwHqBPnz5Oj99MDj3KNLl0PpNL5zM4YR/aWg/kmln6Tsdbx99vaOI8IiKyk7AXyQxiGbvvdLwnEAHmNHkikb3gnNbcir/CXiQvAA4o2On4McBbzrnSJk8k0kDzVpdScE8Rf3prXtBRRBpFqIvEOTeP2L2O35pZRwAzG0tsNfzvgswmUh+rNm9j1ITJLNmwlXfnrQ06jkijCOxmu5ndTWxl+/7x//4s/qHhzrmqHT71cuAm4CMzqwZKgROcc58hEmKbt1UzekIxq7dUABCp0fCW+CmwInHOXV3Pz6sGro+/iaSEykgNP3u0hHlrSmnfMouN5VVU1USDjiXSKEI9tCWSiqJRx9XPzOA/CzeS1zqbv553GADVKhLxVKquIxEJrTvfmMvL01fSMiudwjH5tG+ZBUB1RENb4iddkYgkUeFHi3jw/YVkpBkPXDCU/l3bkpEW+2emKxLxlYpEJEnemLWKW/7vcwDuGDmQo3vnAZCVHvtnpnsk4isViUgSlCzeyC+f+gzn4KoTenPm0G51H8vMiO3oo1lb4isViUiCFqwt4+JHSqiMRDnv8P257JheX/t4ZrqGtsRvKhKRBKwtrWD0xMls2lrNcX078ftT+2P29T1FM9LiVyRRRzSqqxLxj4pEZC+VVUYYM7GY5V9tY1D3dtx37mFkpH/zn5SZ1d0nqY7qqkT8oyIR2QvVNVF+/vhUZq/cwgEdcnl41DBys3Y/mz4j3eKv0xWJ+EdFItJAzjl++/xM3p+/jg4ts5g0ZjgdW2V/62tq75NEdJ9EPKQiEWmge9+ez7NTlpOTmcbDo/Pp0bHlHl+TqSnA4jEViUgDPPHpUu779wLSDP567hAGd29Xr9dlaWhLPKYiEamnd+as4foXZwJw22kDOK5f53q/NjMjfrM9oisS8Y+KRKQepi/bxC+emEbUweXH9uK8w/dv0OtrpwBrLYn4SEUisgeL15cztrCYbdU1jBzSjV8d37vBX2P7okQNbYl/VCQi32JDWSWjJ05mQ3kVIw7uyB0jB3xjwWF9ZGVodbv4S0Uishvbqmq4aFIJizdspX/XNjxwwdC6K4uG0jYp4jMVicguRGqiXP7kVD5bton92rVg4uh8WmXv/eN7MuOztjT9V3ykIhHZiXOOm16ezb/mrKVti0wmjc2nU5uchL6m7pGIz1QkIju5v+hLHv90KVkZaTw8ahi9OrVO+GvWFYmm/4qHVCQiO3h+6nLufnMeZvDncwYzrEf7pHzd2qGtiDZtFA+pSETiPvhiHb95dgYAN57cjx8O2DdpX3v7Fika2hL/qEhEgNkrN3PpY1OJRB3jjj6QMUf2TOrXz9LQlnhMRSLN3vKvtjJmYjFllRFOGdSVa088JOnfQ9N/xWcqEmnWNm+tZvTEYtaWVnJ4z/bcc9ZA0tIavuBwT7Y/j0RFIv4JfZGY2TAze93M5pjZTDObbGZnBZ1LUl9FdQ2XPFLCgrVl9O7civE/HUZ2RnqjfC9N/xWfhbpIzKwH8A6wHhjgnBsATACeNrNTAowmKS4adfz66elMXryRLm1yKBwznLYtMhvt+2mLFPFZqIsEOAloA/yPcy4C4Jz7O7AFOC/IYJLabn9tDq/OXEXr7AwKx+bTtV2LRv1+mRraEo+FvUgi8fd1e1NYbMe8NKBxxiDEew99sJCHP1xEZrrx4IVDOaRLm0b/npr+Kz4Le5E8BcwFrjezVmaWBlwHZAN/DzSZpKRXZ6zitlfnAHDPWYM4olfHJvm+mrUlPtv7XeiagHNui5l9H5hI7D5JGbAZON45997uXmdm44BxAHl5eRQVFTVBWv+VlZWl9Lmct7GGu4srADi7dyZtN31BUdEXTfK9ly6uBmDR4iUUFa0GUv98ho3OZ3BCXSRm1ofYzfbXgPZABXA28LyZXeCce31Xr3POjQfGA/Tp08cVFBQ0TWDPFRUVkarncv6aUq544GMiDkZ99wBuPrX/Xj1XZG8tylwE8z6nS9duFBT0B1L7fIaRzmdwwj60dSvQDvilc26rcy7qnHsKeB+YZGahLkIJhzVbKhg9YTJbKiL8oH9nbjylaUsEdrxHoqEt8U/Yi2QAsNw5t22n4/OBPCC5+1iId0orqhk1YTIrN1cw9IB9+PNPDiO9ERYc7om2SBGfhb1I1gL77uLK4wDAAV81fSRJFVWRKJc+NpW5q0s5sGNLHvrpMHIyg5nsl5mh6b/ir7AXyV+IrSP5fXzaL2Z2DHAG8E/n3Pogw0l4Oee45rkZfLhgPR1bZTNp7HD2aZkVWJ6MNK1sF3+F+h6Dc+5ZMzsRuBb43MxqgCjwO+C+QMNJqN395jxemLaC3Kx0Jo7Op3v73EDzaPqv+CzURQLgnHsTeDPoHJI6Hv3PEu4v+pL0NOP+84cwoFvboCORpaEt8VjYh7ZEGuSt2au56aVZAPzxjAEU9OkUcKIYbdooPlORiDemLPmKy5+cRtTBlcf15uxh3YOOVEfTf8VnKhLxwsJ1ZVw8qZjKSJRzhnXniu/3CjrS12jTRvGZikRS3rrSSkZNnMxXW6s5pk8et59+aJMvONyT2iuSiIa2xEMqEklp5ZURxhYWs2zjNgZ2a8tfzxtCRnr4fqw1a0t8Fr5/cSL1FKmJ8osnpjJzxWb2b5/LhNH5tMwO50RE3SMRn6lIJCU557j+xVm8O28d++RmMmnscDq2yg461m5l6YpEPKYikZR03zsLeKp4GTmZaTw8Op+eHVsGHelbZdTebI/oHon4R0UiKefp4mXc+6/5pBn85dwhDNl/n6Aj7ZHukYjPVCSSUormreW3L8wE4Pc/PpTj+3UOOFH9aGhLfKYikZQxc/lmfv74VGqijp8XHMQF3zkg6Ej1tn33Xw1tiX9UJJISlm3cypjCYrZW1XD6Yftx9Q/6BB2pQTS0JT5TkUjofVVexagJk1lfVslRvTpy58iBoVtwuCcZ8YdpRaKOaFRXJeIXFYmEWkV1DRdNKmbh+nIO6dKaBy4YQlZG6v3Ymtn2bVKiuioRv6Tev0hpNmqijiuenMbUpZvo2jaHwjHDaZ2TGXSsvaZtUsRXKhIJJecct7wym7c+X0ObnAwKxw6nS9ucoGMlRPdJxFcqEgmlB99fyCOfLCErPY3xPx1G786tg46UMG2TIr5SkUjovPTZCu54fS4Afzp7EN85sEPAiZIjK11TgMVPKhIJlY8XrOeqZ6YDcP2P+nLKoK4BJ0qe2l2JqyO6IhG/qEgkNOas2sLPHp1CdY1j7JE9uXjEgUFHSqraWVsRzdoSz6hIJBRWbtrGmInFlFZG+NGAfbn+R32DjpR0dfdItHGjeEZFIoHbvK2a0RMns3pLBcN7tOdPZw8iLS21FhzWR+36F83aEt+oSCRQlZEaxj1Swvw1ZfTq1IrxPx1KTmZ60LEahab/iq9UJBKYaNRx1TMz+HTRRjq1zqZwTD7tcrOCjtVoardJ0fRf8Y2KRAJzxxtzeWX6SlplZzBxTD7d9skNOlKj2j60pXsk4peUKBIzG2lm75vZFDNbaGYlZnZh0Llk7038aBHj319IRprxwAVD6N+1bdCRGt32LVJ0RSJ+CX2RmNmVwO+A85xzQ4E+wHzg+4EGk732+sxV/P7/PgfgrjMHMuLgvIATNY26TRtVJOKZjKADfBsz6wHcARzlnFsO4JyrNrOrAH9WqjUjkxdt5Jf//Azn4Oof9OGMId2CjtRktm+RoqEt8UuoiwS4ENjknCve8aBzbiWwMphIsrcWrC3lkkdKqIpEOf/w/fl5wUFBR2pSWVrZLp4K+9DWEcDi+D2SD8xsrpl9bGZjgw4mDbN2SwWjJhSzeVs1x/XtxC2n9k+5h1MlKkNDW+KpsF+RdAd6AFcBpwNrgZHAk2a2r3Pu9l29yMzGAeMA8vLyKCoqapKwvisrK9urc7kt4vjjpxWsKI1yYNs0zupWxocfvJ/8gCG3bk0lAJ/PnUfR1oV7fT5l13Q+gxP2IskBWgJXO+dWx489Y2Y/Aa4zs3udc1t3fpFzbjwwHqBPnz6uoKCgqfJ6raioiIaey+qaKGMLi1laupUeHXJ55tIj6NAqu3EChlzRltmwbDE9DuxFwVE99+p8yu7pfAYn7ENbpfH3n+10fBqQC/Rr2jjSEM45rn1uJh98sZ4OLbOYNHZ4sy0R0BYp4q+wF8nc+Pudc9bs5riEyL1vz+e5qctpkZnOhNH5HNChZdCRAqXpv+KrsP8ifiX+fuBOxw8FtgGzmzaO1NcTny7lvn8vID3N+Nv5hzGoe7ugIwUuI03Tf8VPYS+SfwLFwG1m1grAzEYAZwK3O+fKgwwnu/bOnDVc/+JMAG477VCOPaRzwInCoXZoSyvbxTehvtnunKsxsxOBO4HZZlYBVAK/cM79I9h0siufLdvEL56YRtTBFd8/mHOH7x90pNDQ0Jb4KtRFAuCc2whcEnQO2bPF68sZW1jMtuoazhrajSuPOzjoSKGyfRt5DW2JX8I+tCUpYkNZJaMmTmZjeRVH987jD2cMaHYLDvdk+xYpuiIRv6hIJGFbqyKMnVTCkg1bOXS/Ntx//pC6X5qynbZIEV/Ve2jLzM4Cfgy0AhYBLzjnmt/yZPmaSE2Uy5+YxvRlm+i2TwsmjM6nVXboR0wDoS1SxFf1+rPRzG4iNoPqFOBAYpspFpnZVDPr04j5JMScc9zw0mzembuWdrmZTBo7nE6tc4KOFVp190iiukcifqnv+MNlwDNAB+fcQOdcR2AEUAZMNrO+jRVQwutv7y7gyclLyc5I4+FRwzgor1XQkUItU0Nb4qn6Fklb4GHnXKT2gHPuI+B7wFTgrkbIJiH23JTl3PPWfMzgzz85jKEHtA86UuhlZWhoS/xU3yJZTmwn3q9xzjngL0BBEjNJyL0/fx3XPDcDgJtP6c+Jh3YJOFFq0PRf8VV9i+QB4CYz2283H69IUh4JuVkrNnPpY1OIRB0/+96BjDqiR9CRUsb2LVJ0RSJ+qW+R3EtsX6tZZnarmX3XzLqbWQFwG6BV5s3A+m1RxhQWU15Vw48Hd+WaHxwSdKSUUju0pS1SxDf1KhLnXA2xGVt3Entg1IfAYuAdYs8LWWxmh5mZ5n16atPWKv5UUsG60kq+e2AH7jpzIGlpWnDYEBraEl/Ve9WYcy7inLsD6ELsEbi/Ah4nNnPrfqAEKDWzyY0RVIJTUV3DJY+UsKrccUiX1jz406FkZ6QHHSvlbC8SXZGIXxp8BRG/wf5p/A0AM8sFDgOGAUOSlk4CF406rvznZxQv/or2OcbEMfm0yckMOlZK0hYp4qukDEXFH3f7UfxNPOGc49ZXP+f1WatpnZPBr4Zmsm/bFkHHSlna/Vd8pQ2RZLce+mAREz9aTFZ6Gg9eOJRurfXjkojaK5KI7pGIZ/SbQXbp5ekruf21OQDcc/YgjjioY8CJUp/ukYivVCTyDZ98uYGrnp4OwHUnHcKpg7oGnMgPtbv/VmmLFPGMikS+Zt7qUsY9WkJVTZTRR/TgkhEHBh3JG5l1W6RoaEv8oiKROqs3VzB64mRKKyKc2L8LN5zcTw+nSiINbYmvVCQCwJaKakZPnMyqzRUMO2Af/vcng0nXgsOkyoifz0jUEZtFL+IHFYlQFYnyX49OYe7qUg7Ka8lDo4aRk6kFh8lmZjtMAVaRiD9UJM1cNOr4zbPT+fjLDeS1zqZwzHDa5WYFHctbGt4SH6lImrm73pzHi5+tpGVWOhNH59O9fW7QkbymIhEfqUiasUc+Wczf3/uSjDTjgQuGcuh+bYOO5D1tkyI+UpE0U2/OXs1NL88G4I9nDODo3nkBJ2oedI9EfJRyRWJmH5iZM7MeQWdJVVOWfMUVT07DOfjV8b05a9g3Hn4pjWT7Nim6IhF/pFSRmNlI4Kigc6SyL9eVcdGkYiojUc4d3p3Lj+0VdKRmRRs3io9SpkjMLAv4I/Ba0FlS1drS2ILDTVurOfaQTtz640O14LCJ1d0jiWhoS/yRMkUCXEbs4VnFQQdJReWVES4qLGHZxm0M6taWv553GBnpqfS/3w9ZGZq1Jf5Jid8kZtYeuBq4Lugsqai6JsplT0xl5orNHNAhl4dH55ObpaciB6F2dbuKRHySEkUC3Ag85pxbHHSQVOOc4/oXZlE0bx3tW2YxacxwOrbKDjpWs6XntouPQv9nqZn1As4G+jbgNeOAcQB5eXkUFRU1TrgU8OKCKl5cUE1WGlw2II3Fs4pZvJdfq6ysrFmfy2Qo27INgCnTptEjp0LnM4n08xmc0BcJcBdwh3Nuc31f4JwbD4wH6NOnjysoKGikaOH2dPEyXlwwgzSD+y8YxnH9Oif09YqKimiu5zJZHllczOwNa+nbfwDpa+bofCaRfj6DE+oiMbMRwKHAOUFnSTXvzl3Lb1+YCcCtpx2acIlIcuw4/VfbYoovQl0kwPFAOlC8wzTVLvH3r5lZFXCdc05TgncwY/kmfv74VGqijsuOOYjzDz8g6EgSt32LFEdOwFlEkiXUReKcu5HYjfY6ZnYzcBNwkm6+f9PSDVsZW1jMtuoazhiyH1ed0CfoSLKDupvtetyueCRVZm1JPWwsr2LUxMmsL6viqF4dueOMgVpwGDK1Q1uRqIpE/JEyRWJmJ5nZZ8B/xQ+9Fv9vAbZV1XDRpGIWrS+n775teOCCIXWL3yQ8dhzaEvFFqIe2dhS/D6J7IbtQE3Vc8dQ0pi3dxH7tWlA4Jp/WOZlBx5Jd0NCW+Eh/sqY45xw3vzybtz9fQ5ucDCaNzadzG93GDSttkSI+UpGkuAfe+5JH/7OErIw0HhqVT69OrYOOJN9CW6SIj1QkKeyFacu56415mMH/njOY4T3bBx1J9kBbpIiPVCQp6qMF6/nNszMAuOFH/ThpwL4BJ5L60NCW+EhFkoI+X7mFnz06heoaxyUjejL2qJ5BR5J60oOtxEcqkhSzYtM2xhROpqwywskD9+W3P6z3XpYSAhraEh+pSFLI5q3VjJ4wmTVbKjm8Z3v+dPYg0tK04DCVbF9HoisS8YeKJEVUVNdwyaMlfLG2jN6dWzH+wmFkZ2jbv1RTt7JdRSIeUZGkgGjU8etnpjN50UY6t8mmcMxw2uZqwWEq0tCW+EhFkgL+8NocXp2xilbZGRSOGU7Xdi2CjiR7SUNb4iMVScg9/OEiHvpwEZnpxoMXDqXvvm2CjiQJ0BYp4iMVSYi9OmMVt736OQB3nTmQI3t1DDiRJCorQ9N/xT8qkpD6dOEGrvznZzgHvzmxD6cf1i3oSJIEGWm6RyL+UZGE0BdrSrnkkRKqaqJc+J0DuPR7BwUdSZJk+812XZGIP1QkIbNmSwWjJxazpSLC8f06c/Op/fVwKo9oaEt8pCIJkdKKakZPLGbFpm0ctn877vvJYaRrwaFXNP1XfKQiCYmqSJRLH5vKnFVb6NmxJQ+PyqdFlhYc+kZDW+IjFUkIOOe49rkZfLhgPR1bZTFpzHDat8wKOpY0gtqV7VpHIj5RkYTAPW/N4/lpK8jNSmfC6Hz275AbdCRpJLVXJBENbYlHVCQBe+w/S/jbu1+Snmb87fwhDOzWLuhI0og0tCU+UpEE6O3P13DjS7MA+MPph3JMn04BJ5LGpiIRH6lIAjJ16Vdc/uRUog5++f2DOSd//6AjSRPIqt1rS1ukiEdUJAFYtL6ciyeVUFEd5exh3fjv4w4OOpI0kcy6dSS6RyL+UJE0sfVllYyaMJmN5VV8r3cet58+QAsOm5HtW6ToikT8kRF0gD0xs8HAZcAQYnkzgX8Btzrn1gWZraG2VkW4qLCYpRu3MmC/ttx//pC6MXNpHuoebBV1OKerEvFDKvwWewpoDxztnBsEHA+cAHxkZinzYI5ITZRfPDGN6cs30719CyaMzqdlduh7XJLMzOrKRKNb4otUKBKAa5xz5QDOuRXA3cDBwEmBpqon5xw3vDSLf89dyz65mUwaM5y81tlBx5KA1K0l0eiWeCIV/iQe6Jyr2unYyvj7fZo6zN74y78X8OTkZWRnpPHQqHwOzGsVdCQJUKxIalQk4o3QX5HsokQAegMOeL+J4zTYMyXL+J+355NmcN+5hzH0gJToPmlEdfdJdI9EPBH6ItmZmaUDY4GHnXPzg87zbd6bv47fPj8TgJtP7c8P+ncJOJGEQe3QliZuiS9SYWhrZzcAEeDK3X2CmY0DxgHk5eVRVFTUNMl2sHhzDXdMriAShZN6ZrJ/5WKKihY3eY5kKisrC+Rc+iZSVQnAlrKtOp9JpJ/P4KRUkZjZGOBsoMA5V7a7z3POjQfGA/Tp08cVFBQ0TcC4ZRu3cvUDH1NRA6cN7sr/nD2YNA+eK1JUVERTn0sftZlSxLpt5WS3yNX5TCL9fAYnZYa2zOxC4NfAsc65tUHn2Z2vyqsYNXEy60orOeKgDtx15iAvSkSSZ/usLd0jET+kRJGY2QXANcBxzrnV8WMnx4ewQqOiuoaLHylh4bpyDunSmr9fOJSsjJQ4xdKEan8mIuoR8UToh7bM7HzgH8TujRy3w3YiI4BVQeXaWU3U8cunpjFlyVfs2zaHwjHDaZOTGXQsCaGM+BWqbraLL0JfJMBfgBxiixB3dksTZ9kl5xy/f2U2b85eQ+ucDCaNHU6XtjlBx5KQ0oJE8U3oi8Q51z7oDHsy/v2FTPpkCVnpaYy/cBi9O7cOOpKEWO3QVo3WkYgnNICfoJc+W8EfX58LwD1nD+K7B3UIOJGEna5IxDcqkgR88uUGrnpmOgC/O6kvpw7qGnAiSQXbdwAOOIhIkqhI9tLc1VsY92gJ1TWO0Uf04OIRPYOOJCkiI12ztsQvKpK9sGrzNkZPKKa0IsIPD+3CDSf308OppN6y6rZIUZOIH1QkDbSloprRE4pZvaWC/B77cO85g0nXgkNpgO2bNgYcRCRJVCQNUBmp4WePTGHemlIOymvJP346jJzM9KBjSYrRpo3iGxVJPUWjjqufmcEnCzeQ1zqbwjHDaZebFXQsSUGatSW+UZHU051vzuXl6StpmZVO4Zh8urfPDTqSpCg9j0R8oyKph0kfL+bB9xaSkWb8/cKh9O/aNuhIksI0tCW+UZHswRuzVnPzK7MBuHPkQEYcnBdwIkl1GtoS36hIvsWUJRv55VPTcA6uOqE3I4d2CzqSeGD7FikBBxFJEhXJbny5royLJpVQGYly3uH7c9kxvYKOJJ7YvrJdTSJ+UJHswtrSCkZNmMymrdUc17cTvz+1vxYcStJoaEt8oyLZSVllhLGFxSz/ahuDurfjvnMPq9vSQiQZan+eNLQlvtBvyB1U10S57PGpzFqxhR4dcpkwahi5WaHfaV9STJY2bRTPqEjinHP87oWZvDd/HR1aZlE4ZjgdWmUHHUs8VDe0pXUk4gkVSdy9//qCp0uWk5OZxsOj8+nRsWXQkcRTWkcivlGRAE9OXsp973xBmsFfzx3C4O7tgo4kHtPNdvFNsy+Sf89dw/UvzgLgttMGcFy/zgEnEt9p91/xTbMukunLNnHZ49OoiTouP7YX5x2+f9CRpBnI1PNIxDPNtkiWbChnbGEx26prGDmkG786vnfQkaSZ0NCW+KZZFsmGskpGTZjMhvIqRhzckTtGDtCCQ2kyWRmxnzWtIxFfNLsi2VZVw0WTSli8YSv9u7bhgQuG1v2FKNIUdEUivmlWv0EjNVEuf3Iany3bxH7tWjBxdD6tsrXgUJpWRlptkeiSRPzQbIrEOcfNr8zmX3PW0LZFJpPG5tOpTU7QsaQZ0tCW+Cb0RWJmnczscTObF3971swavJ/7/UVf8th/lpKVkcbDo4bRq1Prxogrskca2hLfhLpIzCwLeBvIAvoD/YBy4F0za1Xfr/P81OXc/eY8zODP5wxmWI/2jRNYpB4ytWmjeCbURQKMAgYC1zjnIs65GuAa4EDg0vp8gW0R+M2zMwC48eR+/HDAvo2VVaRetl+RqEnED2EvkpHAUufcwtoDzrnVwOfxj+3R2q1RIlHHuKMPZMyRPRsppkj9ZWloSzwT9ilLA4H5uzi+CPh+fb6AA04Z1JVrTzwkmblE9lpGfIuUihq4/MlpAafxx9o1FTy3SuczCGEvko7AlF0c3wLkmlkL59y2nT9oZuOAcQC5nXtwaudNvP/+e42btBkoKyujqKgo6BgpL+ocLTOhvBpemb4y6Dh+Wa3zGYSwF8nufOsydOfceGA8wMG9+7jjjz2mSUL5rqioiIKCgqBjeOGl/mU8+6//0Ldf36DmN30QAAAF7klEQVSjeGPO53N0PpPktDsb9vlhL5L1wK7m6bYGtu7qamRnadr5REKoV6dWfKdrBgWD9ws6ijfabvpC5zMgYb/ZPgPosYvjPYGZTRtFRER2JexF8jxwgJn1qD1gZp2BvsBzAWUSEZEdhL1IColdedxpZhlmlgbcQWzW1gNBBhMRkZhQF4lzrgo4HqghtnZkDtAGONY5VxZkNhERiQn7zXacc2uA84LOISIiuxbqKxIREQk/FYmIiCTEnPN74zgzKwXmBZ3DEx2Jre2R5ND5TC6dz+Tp45yr97M2Qn+PJAnmOeeGBR3CB2ZWonOZPDqfyaXzmTxmVtKQz9fQloiIJERFIiIiCWkORTI+6AAe0blMLp3P5NL5TJ4GnUvvb7aLiEjjag5XJCIi0ohUJLJHZravmb1hZrp8FWkmzOwDM3M7bpq7O14WiZl1MrPHzWxe/O1ZM+sWdK5UZGanA58ABwWdJdWZ2WAz+4eZTTGz6Wb2uZndZ2Z5QWdLRWZ2kJndEz+fU8xsfvyX34+CzpbqzGwkcFR9P9+7IjGzLOBtIAvoD/QDyoF3zaxVkNlS1LXENs78KOggHngKaA8c7ZwbROy8ngB8ZGYtAk2Wmn4I/AQ4xzk3FDiE2B89L5vZ9wJNlsLiv0P/CLxW39d4VyTAKGAgcI1zLuKcqwGuAQ4ELg00WWo60jn3RdAhPHKNc64cwDm3ArgbOBg4KdBUqWkFcLNzbgGAcy4K/IHY77UfBxksxV0GlADF9X2Bj0UyEljqnFtYe8A5t5rYNvQjA0uVopxzkaAzeGRg7S+9HayMv9+nqcOkOufcC865h3Y63Cb+fl1T5/GBmbUHrgaua8jrfCySgcQefLWzRcCAJs4iUif+fJ2d9QYc8H4Tx/GOme0H/A2YGn8vDXcj8JhzbnFDXuRjkXQESndxfAuQq7FoCQszSwfGAg875+YHnSdVxW+6LwCWA+nAac65LQHHSjlm1gs4G7i9oa/1sUh2x4IOILKTG4AIcGXQQVKZc+5L51wvoC0wH5huZvWecSR17gLucM5tbugLfSyS9cCutj9uDWx1zm1r4jwi32BmY4j99fdDPTY6OeJXIVcCa4D7A46TUsxsBHAo8MDevN7HbeRnEJsGuLOewMwmziLyDWZ2IfBr4Fjn3Nqg86Sq+DB1hdthnyfnnDOzmcCZZpbtnKsMLmFKOZ7YsGCxWd3gTZf4+9fMrAq4zjm3yynBPl6RPA8csONqTDPrDPQFngsokwgAZnYBsenox8VnE2JmJ5vZuGCTpaTXge/s4ngPYvdEdzW5QXbBOXejc+4g59zg2jfg7/EPnxQ/ttt1JT4WSSGxK487zSzDzNKAO4jN2tqryzaRZDCz84F/EPsZPc7MLogXyylA1yCzpbBbzKwDgMVcDuQD9+14pSKNy8vdf+NXIPcCw4hNrZwF/LdzblmgwVKQmd1N7LJ3f2JrHabHPzR8N9NZZTfMbCO7Xy9yi3Pu5iaMk/LM7EjgYmLFEQFygA3E7o88oSLZO2Z2ErGFnV2AzsAcoCp+lbLr1+hci4hIInwc2hIRkSakIhERkYSoSEREJCEqEhERSYiKREREEqIiERGRhKhIREQkISoSERFJiIpEREQSoiIREZGEqEhERCQhKhKRJmJmvcys2sxu2en4A2ZWambDgsomkggViUgTcc4tAB4CrjSzjgBmdiOx57af7pwrCTKfyN7S7r8iTcjMugBfEtvqfC4wHjjXOfd0oMFEEuDjo3ZFQss5t9rM/pfYo3YzgCtUIpLqNLQl0vS+ALKBT5xzfws6jEiiVCQiTcjMjgUeBD4BjjSzQQFHEkmYikSkiZjZEOBFYjfcC4ClxB5pKpLSVCQiTcDMegGvA28Bl8efd38LcJKZHR1oOJEEadaWSCOLz9T6mNgVyA+cc5Xx4+nALOAr59wRAUYUSYiKREREEqKhLRERSYiKREREEqIiERGRhKhIREQkISoSERFJiIpEREQSoiIREZGEqEhERCQhKhIREUmIikRERBLy/7mXPz6jAa7sAAAAAElFTkSuQmCC\n",
"text/plain": [
"