"
]
},
{
"cell_type": "markdown",
"metadata": {
"lines_to_next_cell": 0
},
"source": [
"# Solving Ordinary Differential Equations"
]
},
{
"cell_type": "markdown",
"metadata": {
"lines_to_next_cell": 0
},
"source": [
"## Objectives\n",
"In this lecture, we will review the main computational techniques to solve a first order differential problem (initial value problem) of the type:\n",
"\n",
"\\begin{equation}\n",
"\\left\\lbrace\n",
"\\begin{array}{lll}\n",
"\\displaystyle{\\frac{dy}{dt}} & = & f(y,t) \\\\\n",
"y(t_0) & = & y_0\n",
"\\end{array}\n",
"\\right.\n",
"\\end{equation}\n",
"\n",
"Specifically you will be able to:\n",
"\n",
"- link the ODE to a physical problem\n",
"\n",
"- develop different computational methods to solve the initial value problem\n",
"\n",
"- understand the concept of stability and accuracy of a computation methods\n",
"\n",
"- describe Euler forward and backward approaches\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"lines_to_next_cell": 0
},
"source": [
"## Introduction\n",
"\n",
"Ordinary differential equations (ODEs) arise in many physical\n",
"situations. For example, there is the second-order Newton’s Second Law of Mechanics $F=ma$ .\n",
"\n",
"In general, higher-order equations, such as Newton’s force equation, can\n",
"be rewritten as a system of first-order equations. So the generic\n",
"problem in ODEs is a set of N coupled first-order differential equations\n",
"of the form defined above.\n",
"\n",
"$$\n",
" \\frac{d{\\bf y}}{dt} = f({\\bf y},t)\n",
"$$\n",
"\n",
"where ${\\bf y}$ is a vector of variables. As an example, the Newton's law of motion could be written as\n",
"\n",
"$$\n",
"M\\frac{d^2 x}{dt^2} = F(x,t)\n",
"$$\n",
"\n",
"Which is a second order differential equation. But nothing can prevent us from describing the same problem through two first-order coupled ODEs:\n",
"\n",
"\\begin{equation}\n",
"\\left\\lbrace\n",
"\\begin{array}{lll}\n",
"\\displaystyle{\\frac{dx}{dt}} & = & v_x(x,t) \\\\\n",
"\\displaystyle{m\\frac{dv_x}{dt}} & = & F(x,t) \\\\\n",
"\\end{array}\n",
"\\right.\n",
"\\end{equation}\n",
"\n",
"That means that we can limit our next discussions to the first order problems."
]
},
{
"cell_type": "markdown",
"metadata": {
"lines_to_next_cell": 0
},
"source": [
"## Euler's methods\n",
"\n",
"In engineering, conceptual models will describe how a certain quantity (mass, energy, position, ...) will evolve: it is about identifying which processes will affect that quantity, and translate them into a mathematical form.\n",
"\n",
"### From the mass balance to the ODE\n",
"#### Source and sink terms\n",
"\n",
"For the sulfates TMF problem, we had identified 4 contributions to the sulfate concentration.\n",
"\n",
"1. $M_{\\text{pit}}$ is the mass of sulfates added over $\\Delta t$ from the mill (mg).\n",
"2. $M_{\\text{mill}}$ is the mass of sulfates added over $\\Delta t$ from the pit (mg).\n",
"3. $M_{\\text{pore}}$ is the mass of sulfates added over $\\Delta t$ from that exchange term with the tailings porewater (mg)\n",
"4. $M_{\\text{discharge}}$ is the mass of sulfates removed over $\\Delta t$ from the discharge flux (mg)\n",
"\n",
"\n",
"Therefore, if the mass of sulfates at a certain time $t_0$ is $m(t_0) = m_0$, the mass of sulfates at $t_0 + \\Delta t$ is:\n",
"\\begin{equation}\n",
"m(t_0 + \\Delta t) = m(t_0) + M_{\\text{pit}} + M_{\\text{mill}} + M_{\\text{pore}} - M_{\\text{dis}}\n",
"\\end{equation}\n",
"\n",
"\n",
"And we have identified these terms:\n",
"\n",
"\\begin{equation}\n",
"\\left\\lbrace\n",
"\\begin{array}{lll}\n",
"M_{\\text{pit}} & = & Q_{\\text{pit}}c_{\\text{pit}} \\Delta t\\\\\n",
"M_{\\text{mill}} & = & Q_{\\text{mill}}c_{\\text{mill}} \\Delta t\\\\\n",
"M_{\\text{pore}} & = & kA_{\\text{bottom}}(c_{\\text{pore}} - c_{\\text{TMF}}) \\Delta t\\\\\n",
"M_{\\text{discharge}} & = & Q_{\\text{dis}}c_{\\text{TMF}} \\Delta t\\\\\n",
"\\end{array}\n",
"\\right.\n",
"\\end{equation}\n",
"\n",
"So we can rewrite the previous equation as:\n",
"\n",
"\\begin{equation}\n",
"\\begin{array}{llll}\n",
"& m(t_0 + \\Delta t) & = & m(t_0) + M_{\\text{pit}} + M_{\\text{mill}} + M_{\\text{pore}} - M_{\\text{dis}} \\\\\n",
"\\Longleftrightarrow & m(t_0 + \\Delta t) & = & m(t_0) + \\left(Q_{\\text{pit}} c_{\\text{pit}} + Q_{\\text{mill}} c_{\\text{mill}} + A_{\\text{bottom}} k \\left( c_{\\text{pore}} - c_{\\text{TMF}} \\right) -Q_{\\text{dis}} c_{\\text{TMF}}\\right) \\Delta t \\\\\n",
"\\Longleftrightarrow & \\frac{m(t_0 + \\Delta t)-m(t_0)}{\\Delta t} & = & Q_{\\text{pit}} c_{\\text{pit}} + Q_{\\text{mill}} c_{\\text{mill}} + A_{\\text{bottom}} k \\left( c_{\\text{pore}} - c_{\\text{TMF}} \\right) -Q_{\\text{dis}} c_{\\text{TMF}} \\\\\n",
"\\end{array}\n",
"\\end{equation}\n",
"\n",
"The link between concentration and mass of sulfates is the volume of water\n",
"\\begin{equation}\n",
"c_{\\text{TMF}}(t) = \\frac{m(t)}{V(t)}\n",
"\\end{equation}\n",
"\n",
"But, since the volume of water is constant ($V_0$), we have:\n",
"\\begin{equation}\n",
"c_{\\text{TMF}}(t) = \\frac{m(t)}{V_0}\n",
"\\end{equation}\n",
"\n",
"Therefore, if we divide the latter equation by $V_0$, we have:\n",
"\\begin{equation}\n",
"\\frac{c_{\\text{TMF}}(t + \\Delta t)-c_{\\text{TMF}}(t)}{\\Delta t} = \\frac{Q_{\\text{pit}} c_{\\text{pit}} + Q_{\\text{mill}} c_{\\text{mill}} + A_{\\text{bottom}} k \\left( c_{\\text{pore}} - c_{\\text{TMF}} \\right) -Q_{\\text{dis}} c_{\\text{TMF}}}{V_0}\n",
"\\end{equation}\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"lines_to_next_cell": 0
},
"source": [
"So, after each time step, we can compute the evolution of $c_{\\text{TMF}}$. How should we choose the timestep?\n",
"\n",
"\n",
"1. It is irrelevant\n",
"2. It should be as small as possible\n",
"3. The bigger the faster\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"lines_to_next_cell": 0
},
"source": [
"The choice of the timestep matters. Look at the previous equation:\n",
"\n",
"\\begin{equation}\n",
"\\frac{\\color{blue}{c_{\\text{TMF}}}(t + \\Delta t)-c_{\\text{TMF}}(t)}{\\Delta t} = \\frac{Q_{\\text{pit}} c_{\\text{pit}} + Q_{\\text{mill}} c_{\\text{mill}} + A_{\\text{bottom}} k \\left( c_{\\text{pore}} - {\\color{red}{c_{\\text{TMF}}}} \\right) -Q_{\\text{dis}} \\color{red}{c_{\\text{TMF}}}}{V_0}\n",
"\\end{equation}\n",
"\n",
"To compute the solution $c_{\\text{TMF}}(t+ \\Delta t)$, we use the value of $c_{\\text{TMF}}$. But, at what time? Intrinsically, at every time between $t$ and $ t + \\Delta t$, the value of $c_{\\text{TMF}}$ changes.\n",
"\n",
"So the timestep should be small!"
]
},
{
"cell_type": "markdown",
"metadata": {
"lines_to_next_cell": 0
},
"source": [
"#### ODE\n",
"\n",
"Going to the limit $\\Delta t \\rightarrow 0$ (i.e. $\\Delta t= dt$, $c(t+\\Delta t)-c(t)= dc$), the equation becomes:\n",
"\\begin{equation}\n",
"\\frac{dc_{\\text{TMF}}}{dt} = \\frac{Q_{\\text{pit}} c_{\\text{pit}} + Q_{\\text{mill}} c_{\\text{mill}} + A_{\\text{bottom}} k \\left( c_{\\text{pore}} - c_{\\text{TMF}} \\right) -Q_{\\text{dis}} c_{\\text{TMF}}}{V_0}\n",
"\\end{equation}\n",
"which is a 1$^{st}$ order linear ODE.\n",
"\n",
"#### Asymptotic behavior\n",
"\n",
"Let us build a bit of physical sense is often a good idea. One should realize that this problem is bound to have a steady solution (a solution which does not evolve in time). This is not always true. But, in this case, at some point, the concentration in the TMF will be so that the sink term will exactly compensate the source terms.\n",
"\n",
"This happens when the derivative is zero!\n",
"\n",
"\\begin{equation}\n",
"\\begin{array}{lll}\n",
"& \\frac{dc_{\\text{TMF}}}{dt} & = & 0 \\\\\n",
"\\Longleftrightarrow & 0 & = & \\frac{Q_{\\text{pit}} c_{\\text{pit}} + Q_{\\text{mill}} c_{\\text{mill}} + A_{\\text{bottom}} k \\left( c_{\\text{pore}} - c_{\\text{TMF}} \\right) -Q_{\\text{dis}} c_{\\text{TMF}}}{V_0} \\\\\n",
"\\Longleftrightarrow & c_{\\text{TMF}} & = & \\frac{Q_{\\text{pit}} c_{\\text{pit}} + Q_{\\text{mill}} c_{\\text{mill}} + kA c_{\\text{pore}} }{Q_{\\text{dis}}+kA} \\equiv c_\\infty\n",
"\\end{array}\n",
"\\end{equation}\n",
"\n",
"The asymptotic/final concentration (the steady-state concentration) $c_\\infty$ corresponds to the final solution. \n",
"\n",
"Can you try to guess which values $c_\\infty$ could be?\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us analyse this by investigating limit cases:\n",
"\n",
"1. Very fast production from the porewater: $k \\rightarrow \\infty$\n",
"\\begin{equation}\n",
"\\begin{array}{lll}\n",
"c_\\infty & = & \\frac{Q_{\\text{pit}} c_{\\text{pit}} + Q_{\\text{mill}} c_{\\text{mill}} + kA c_{\\text{pore}} }{Q_{\\text{dis}}+kA} \\\\\n",
" & = & \\frac{\\color{blue}{Q_{\\text{pit}} c_{\\text{pit}}} + \\color{blue}{Q_{\\text{mill}} c_{\\text{mill}}} + kA c_{\\text{pore}} }{\\color{blue}{Q_{\\text{dis}}}+kA} \\\\\n",
" & \\approx & \\frac{kA c_{\\text{pore}} }{kA} \\\\\n",
" & \\approx & c_{\\text{pore}} = 2000 \\text{ mg/L}\n",
"\\end{array}\n",
"\\end{equation}\n",
"\n",
"2. Very slow production from the porewater: $k \\rightarrow 0$\n",
"\\begin{equation}\n",
"\\begin{array}{lll}\n",
"c_\\infty & = & \\frac{Q_{\\text{pit}} c_{\\text{pit}} + Q_{\\text{mill}} c_{\\text{mill}} + kA c_{\\text{pore}} }{Q_{\\text{dis}}+kA} \\\\\n",
" & = & \\frac{Q_{\\text{pit}} c_{\\text{pit}} + Q_{\\text{mill}} c_{\\text{mill}} + \\color{blue}{kA c_{\\text{pore}}} }{Q_{\\text{dis}}+\\color{blue}{kA}} \\\\\n",
" & \\approx & \\frac{Q_{\\text{pit}} c_{\\text{pit}} + Q_{\\text{mill}} c_{\\text{mill}}}{Q_{\\text{dis}}} \\\\\n",
" & \\approx & 256.82 \\text{ mg/L}\n",
"\\end{array}\n",
"\\end{equation}\n",
"\n",
"Therefore, the asymptotic concentration should never be outside these two bounds!"
]
},
{
"cell_type": "markdown",
"metadata": {
"lines_to_next_cell": 0
},
"source": [
"### Analytical Solution\n",
"\n",
"The ODE can be rewritten in the following way:\n",
"\n",
"\\begin{equation}\n",
"\\frac{dc_{\\text{TMF}}}{dt} + \\lambda c_{\\text{TMF}} = Q\n",
"\\end{equation}\n",
"with\n",
"\\begin{equation}\n",
"\\left\\lbrace\n",
"\\begin{array}{lll}\n",
"Q & = & \\frac{Q_{\\text{pit}}c_{\\text{pit}} + Q_{\\text{mill}}c_{\\text{mill}}+ kAc_{\\text{pore}}}{V_0} \\\\\n",
"\\lambda & = & \\frac{Q_{\\text{dis}}+kA}{V_0}\n",
"\\end{array}\n",
"\\right.\n",
"\\end{equation}\n",
"\n",
"You can notice that\n",
"\\begin{equation}\n",
"c_\\infty = \\frac{Q}{\\lambda}\n",
"\\end{equation}\n",
"\n",
"The associated homogeneous problem is\n",
"\\begin{equation}\n",
"\\begin{array}{ll}\n",
"& \\frac{dc_{\\text{TMF}}}{dt} + \\lambda c_{\\text{TMF}} = 0 \\\\\n",
"\\Longleftrightarrow & \\frac{dc_{\\text{TMF}}}{dt} = -\\lambda c_{\\text{TMF}} \\\\\n",
"\\end{array}\n",
"\\end{equation}\n",
"\n",
"In other words, we are trying to find a function which is proportional to its opposite.\n",
"\n",
"Any idea?\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"lines_to_next_cell": 0
},
"source": [
"The solution of the homogeneous problem $c_H$ is\n",
"\\begin{equation}\n",
"c_H(t) = A \\mathrm{exp}(-\\lambda t)\n",
"\\end{equation}\n",
"\n",
"The solution to the particular problem should arise from our physical sense. We know that the solution should reach a steady-state, so a particular solution is a constant solution $c_p(t) = c$.\n",
"\n",
"By injecting that solution into the differential problem, we find:\n",
"\n",
"\\begin{equation}\n",
"c_p(t) = c = \\frac{Q}{\\lambda} = c_\\infty\n",
"\\end{equation}\n",
"\n",
"which corresponds to our physical sense!\n",
"\n",
"So that the general solution to the ODE is:\n",
"\\begin{equation}\n",
"c_{\\text{TMF}}(t) = c_\\infty + A \\mathrm{exp}(-\\lambda t)\n",
"\\end{equation}\n",
"\n",
"The last thing to do is to find the value of A so that the initial concentration meets our given problem. Usually, initial means at $t=0$, but you can take any other value of time $t_0$, as long as you are consistent with that value later.\n",
"\n",
"We know that:\n",
"\\begin{equation}\n",
"c(0) = c_0 (=93) \\Longleftrightarrow c_\\infty + A = c_0 \\Longleftrightarrow A = c_0 - c_\\infty\n",
"\\end{equation}\n",
"\n",
"So the exact solution of the problem:\n",
"\\begin{equation}\n",
"c_{\\text{TMF}}(t) = c_\\infty + \\left( c_0 - c_\\infty \\right) \\mathrm{exp}(-\\lambda t)\n",
"\\end{equation}\n",
"\n",
"So, what would happen if you arbitrarly chose $t_0$ to be something else than 0? For example, if you want $t_0$ to be the age of the universe, then, $c_{\\text{TMF}}$ has to match 93 mg/L after 4.343 $\\times$ 10$^{17}$seconds (13.772 billions of years!). That would give you a value A. As long as you stick to that value of A and to that initial $t_0$, the following solution will be the same. But computationnally, the difference between 13 billions of years and 13 billions + 1 year might be lost in rounding erros... So it's probably not the greatest idea!\n",
"\n",
"The fact that we know the exact solution to a problem allows us to assess the accuracy of a computational approach. Let us take a look at the solution!"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"6.358024691358025e-09\n",
"3.2469135802469134e-06\n"
]
}
],
"source": [
"# Parameters definition\n",
"import numpy as np\n",
"from matplotlib import pyplot as plt\n",
"\n",
"c_pit = 50\n",
"c_pore = 2000\n",
"c_mill = 700\n",
"Q_pit = 30\n",
"Q_mill = 14\n",
"Q_dis = 44\n",
"V0 = 8.1e9\n",
"c0 = 93 # initial concentration\n",
"k = 2.5e-5\n",
"Area = 3e5\n",
"#\n",
"# Eqn 13\n",
"#\n",
"c_inf = (Q_pit * c_pit + Q_mill * c_mill + Area * k * c_pore) / (Q_dis + Area * k)\n",
"A = c0 - c_inf\n",
"lam = (Area * k + Q_dis) / V0 # units are here in s^{-1}\n",
"#\n",
"# Eqn 15\n",
"#\n",
"Q = (Q_pit * c_pit + Q_mill * c_mill + Area * k * c_pore) / V0\n",
"print(lam)\n",
"print(Q)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Asymptotic concentration is: 510.6796116504854 mg/L,\n",
"exact concentration after 50 years is 510.6611233781182\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAEKCAYAAAD+ckdtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8FdX9//HXJwuEQIAQwiJbkH0JSSBsIioqatWqiFtrVVRcfrVq7bdWuqm1rfVbtFarX1sRxBV3lqptVYSCIiCrLIbVAIlAIBAIZCHL+f1xJyFglhuSm5vl/Xw84p05M3Pmk5HkkzNz5hxzziEiIiKNS0iwAxAREZHapwQvIiLSCCnBi4iINEJK8CIiIo2QEryIiEgjpAQvIiLSCCnBi4iINEJK8CIiIo2QEryIiEgjFBbsAGqiffv2Li4uLthhiIiI1ImVK1fud87F+rNvg07wcXFxrFixIthhiIiI1Akz2+HvvrpFLyIi0ggpwYuIiDRCSvAiIiKNkBK8iIhII6QELyIi0ggpwYuIiDRCSvAiIiKNUIN+D75WfbsGUj4IdhQiUgscDueg2DkAih045ytznLhc7O3rnPO2HV8udoBXT+lxvhP4Pr2y0vO6kvOXs63McSXlrszG0mO8leN1lJaetN/x/5y4b0nx8XOfcK6KrlkVO1S22VV1cFV1uSq2f2d3V3al7s51CiwsgpE3PVqjOk6VEnyJPetg0dRgRyFS79Ts11vg6vKnwpCTPkXqWrZFAkrwwTX0Bt+XSANWXOw4cqyQ7LxCsvMKTvg8nFfIkbxCco4VknOsiJxjReR6y7kFRSeUlaznHiuisLjmablZWAjNQ0MICzVCQ0IIDzXCQo2wkBDCQoyw0JJPIzwkhNCS5VDfcvgJ+x6vIzTECAsxQkKMEDNCzQgxjq97nyEGoSGGmRH6ne14+5i3D2WOK7NPmboM335WdhkwK7sMnLSt7LFw0vHfqcurj7Kf5Rx/0rYQA07ar8RJq9hJO3x3+8nHV1xhoM918vaaHFvT2Kqj9akfWmMBTfBmlgpkA0VAoXMu2czaAW8CcUAqcI1z7qD5rvBTwMVADjDJObcqkPGJ1FfFxY7svEIO5BzjwNFjHDx6jAM5J30eLeBQ7jEO5x5P4keOFVZ5q9UMIsNDadEsjMhmoUQ2CyUi3PcZHRnuKw8PpYW3rWR7RHioL1F7X77l42Xlr4fQLDTkO79ARSTw6qIFP845t7/M+hRgvnPuMTOb4q0/AHwP6ON9jQSe8z5FGo2ComL2ZeeTkZ3P3sN5ZGTnk3E4j4zD+ezN9n1mZOdx4OgxKmo4NwsNoV3LZrSNDCc6shk9YiKJigindYsw32dEGFERvuUTP8OIah5ORLgSrkhTEIxb9JcD53jLLwEL8SX4y4GXna+3xlIza2tmnZ1zu4MQo8gpySsoIj0rl/SDuaQdzCXtYA7pWceXM7Lzv9PCDjFo36o5HVo3p3ObCBK6tSGmZXOiWzajXUtfEm/XslnpZ2SzUCVoEalSoBO8Az4yMwf8wzn3PNCxTNLeA3T0lrsAu8ocm+aVnZDgzex24HaA7t27BzB0kfIVFhWTdjCXb/YfZdu+I3yz/yjb9x1l+/4j7D2cf8K+YSFG57YRdG0bydg+sXRp24KOrSPo2Lo5HaJ8n+1aNiMsVN3ARKR2BTrBn+mcSzezDsDHZpZSdqNzznnJ32/eHwnPAyQnJ9d6p1yREs450rNy+Xp3Nim7D5OyJ5tNe7PZkXmUgqLj//TatAjn9NiWnNk7lh4xkXSNbkHXaN9nx9YRhIaotS0idS+gCd45l+59ZpjZbGAEsLfk1ruZdQYyvN3TgW5lDu/qlYkEXHGxY9u+I6zZlcW69EOk7M7m6z2Hyc4rLN2nR0wkfTtGcf6Ajpwe25LT27fk9NhWtGvZLIiRi4iUL2AJ3sxaAiHOuWxv+QLgEWAecBPwmPc51ztkHvATM3sDX+e6Q3r+LoGyLzuf1TsPsmZXli+ppx0iO9+XzFs1D6N/pyiuSOxC/85R9O/Umn6domjVXG+VikjDEcjfWB2B2V5noDDgdefcv83sS+AtM7sV2AFc4+3/Ib5X5Lbie03u5gDGJk1MRnYey7YfYOn2TJZuz2TbvqOA7xl5/85RXJ50Gondokns1pbT27ckRLfVRaSBC1iCd85tBxLKKc8Eziun3AF3BSoeaVpyjhWyZGsmCzdn8MW24wm9VfMwhsdFc3VyN5J7RDO4SxsiwkODHK2ISO3TPUdpNL7Zf5QFKRks2JTBsu0HOFZUTGSzUEb2bMc1yd0YdXoMg05rrR7rItIkKMFLg+Wc4+vd2fxr/W4+XLe7tJXeK7YlN47uwbj+HUiOi6Z5mFroItL0KMFLg+KcY+Puw3zw1W7+tX4P3+w/SojBqNNjuHF0HOf270C3dpHBDlNEJOiU4KVB2H8knzmr03lnZRope7IJDTHO6BXDbWNP58JBHYlp1TzYIYqI1CtK8FJvFRYVMz8lg7dXpLFwUwaFxY6Erm34/eWDuGTIaXr/XESkEkrwUu8cOHqMN77cyWtLd5KelUtsVHNuPbMnE4d1pW/HqGCHJyLSICjBS72x4dtDzPgslX9+9S3HCos5o1cMv710IOcP6KCe7yIi1aQEL0G3/JsDPLtgK//dvI/IZqFcm9yNG0f3oI9a6yIip0wJXoLCOcfCTft4dsFWVuw4SEzLZvzion78aFQPWkeEBzs8EZEGTwle6tySbfv58783sWZXFl3atuB3lw3imuRutGim99VFRGqLErzUmXVph/jzf1JYvGU/ndtE8NiV8Uwc1pVwPV8XEal1SvAScLsP5fKnD1OYt/Zb2kaG8+uLB3DD6B4aA15EJICU4CVg8guLeGHxNzzz6VaKnOOucb244+xeesYuIlIHlOAlIBZsyuB38zaQmpnD+IEd+e0lA+keoyFkRUTqihK81KoDR4/xyD83MGfNt5zeviUv3TKCs/vGBjssEZEmRwleaoVzjg/W7eahuRs4lFvAPef14a5xvTSTm4hIkCjBS40dOHqMX723jn9v2EN8lza8OnkkAzq3DnZYIiJNmhK81MiSrfu57601HDh6jAcu6s9tY3tqWFkRkXpACV5OSUFRMU98tJl/LNpGz/YtmX7TcAZ3aRPssERExKMEL9W2+1Au/+/VVazZlcUPRnTjt5cOJLKZ/imJiNQn+q0s1bJ0eyZ3vbaKvIIinv3hUC4Z0jnYIYmISDmU4MUvzjlmfJ7Kox9+TY+YSN68YRS9O2i2NxGR+koJXqp0rLCYKe99xXur0hk/sCN/uSaBKI1GJyJSrynBS6UO5RRw56sr+WJ7Jj89vw/3nNuHkBALdlgiIlIFJXip0K4DOdw880t2ZB7lyWsTmJDUNdghiYiIn5TgpVzr0w8x6cUvyS8s4qVbRnBGr/bBDklERKpBCV6+Y+WOg0x6cTlRzcOYddsZ9OmoznQiIg2NErycYMnW/Ux+eQUdoprz2m2j6NK2RbBDEhGRU6AxRaXUpyl7mTTzS7pGt+CtO0YruYuINGBqwQsAC1IyuOOVlfTv1JqXbxlBdMtmwQ5JRERqQAle+Hzrfu54dSX9OkXx6uSRtGmhd9xFRBo63aJv4lakHmDySyvoGdOSV25RchcRaSyU4Juwr9KymPTil3RuE8Grk0fqtryISCOiBN9E7cg8ys0vfknbyHBeu20ksVHNgx2SiIjUIiX4JujA0WNMevFLipzjpVtG0LmNesuLiDQ26mTXxOQVFDH5pS9Jz8rl9ckj6RXbKtghiYhIAKgF34QUFzvufWM1q3dl8dS1iSTHtQt2SCIiEiBVtuDNrAMwBjgNyAXWAyucc8UBjk1q2V8+3sx/Nuzlt5cO5HvxnYMdjoiIBFCFLXgzG2dm/wE+AL4HdAYGAr8B1pnZ78ysdVUnMLNQM1ttZu976z3NbJmZbTWzN82smVfe3Fvf6m2Pq/m3JyU+XLebZxZs5ZrkrtwyJi7Y4YiISIBV1oK/GLjNObfz5A1mFgZcCowH3q3iHPcCXwMlfwz8L/Ckc+4NM/s7cCvwnPd50DnX28yu8/a7tjrfjJTv692H+Z+31pLUvS2/v2IwZprPXUSksauwBe+cu7+85O653Dk3xzlXaXI3s67AJcAL3roB5wLveLu8BFxRUqe3jrf9PFMmqrGDR49x+ysraN0ijH/8aBjNw0KDHZKIiNSBU+1k96Sf+/0V+AVQ8rw+BshyzhV662lAF2+5C7ALwNt+yNtfTlFxseO+t9aw91A+f//RMDq0jgh2SCIiUkdONcFX2bI2s0uBDOfcylM8R0X13m5mK8xsxb59+2qz6kZn2uLtLNy0j99eOoCk7tHBDkdEROrQqSZ458c+Y4DLzCwVeAPfrfmngLbeM3yArkC6t5wOdIPSZ/xtgMzvnNi5551zyc655NjY2FMMv/FbueMgU/+ziYvjO/GjUT2CHY6IiNSxCjvZmdk6yk/kBnSsqmLn3C+BX3p1nQP83Dl3vZm9DVyFL+nfBMz1DpnnrX/hbf/UOefPHxJykqycY9wzazWd20bwpyuHqFOdiEgTVFkv+ksDdM4HgDfM7A/AamC6Vz4deMXMtgIHgOsCdP5GzTnHA+9+RUZ2Hu/ceYZmhxMRaaIqS/DPA/8G/uWcS6nJSZxzC4GF3vJ2YEQ5++QBV9fkPALvrEzjPxv28quL+5PQrW2wwxERkSCp7Bn8TcBB4GEzW2Vmz5nZ5WbWso5ik2pKz8rlkX9uZETPdkw+8/RghyMiIkFUYQveObcHmAnMNLMQYCS+Ee1+YWa5wEfOuT/XSZRSpeJix/1vr6XYOZ64OoGQED13FxFpyvyaTc4bd/4L7+tBM2sPXBjIwKR6XvoilSXbMvnTlfF0axcZ7HBERCTI/Jls5m98tzf9IWBFQCKSavtm/1Ee+1cK4/rFct3wbsEOR0RE6gF/3oNvDiQCW7yvIfjeX7/VzP4awNjED845fj17Hc3CQnhsol6JExERH39u0Q8BxjjnigDM7DlgMXAmsC6AsYkf3l2VzpJtmfzhisF01FC0IiLi8acFHw20KrPeEmjnJfz8gEQlfsk8ks8fP9jIsB7R/HBE92CHIyIi9Yg/Lfg/A2vMbCG+UezOAh71Xpf7JICxSRX++MHXHMkv5E9XxqvXvIiInKDKBO+cm25mH3J8cJpfOee+9ZbvD1hkUqnPt+7nvdXp3H1ub/p2jAp2OCIiUs/4O9lMLL6e9KHAKDO7MnAhSVUKiop5aN4GesREcte43sEOR0RE6iF/XpObga+j3QaOz+vugPcCGJdU4uUvdrA14wgv3JhMRHhosMMREZF6yJ9n8KOccwMDHon4JfNIPn/9ZDNn9Y3lvAEdgh2OiIjUU/7cov/CzJTg64nHP9pE7rEiHrx0oN55FxGRCvnTgn8ZX5Lfg++1OAOcc25IQCOT71iffog3vtzFLWN60rtDq6oPEBGRJsufBD8duAHfoDbFVewrAeKc45H3N9Iushn3nNcn2OGIiEg950+C3+ecmxfwSKRSCzZlsPybA/z+8kG0aREe7HBERKSe8yfBrzaz14F/UmbkOuecetHXkaJix5//vYm4mEiu04h1IiLiB38SfAt8if2CMmV6Ta4OzV2TTsqebP72gyTCQ/0dukBERJoyf0ayu7kuApHy5RcW8cRHm4nv0oZL4jsHOxwREWkgKmwOmtlvzKxdJdvPNbNLAxOWlHh16U7Ss3J54KL+Gm9eRET8VlkLfh3wTzPLA1YB+4AIoA+++eE/AR4NeIRN2NH8Qp5dsJUze7fnzD7tgx2OiIg0IBUmeOfcXGCumfUBxgCdgcPAq8Dtzrncugmx6Xr5ix0cOHqM/7mgb7BDERGRBsafZ/BbgC11EIuUcTS/kGmLt3NW31iSukcHOxwREWlg1CW7nnp1qa/1fq8GtRERkVOgBF8P5Rwr5PlF2xnbpz3Deqj1LiIi1acEXw+9tnQnmWq9i4hIDfgzH3wscBsQV3Z/59wtgQur6co9VsQ/Fm1jTO8YkuMqfEtRRESkUv6MZDcXWIzvtbiiwIYj76zcxf4jx3j2XLXeRUTk1PmT4COdcw8EPBKhqNgxbfE3JHZry4iear2LiMip8+cZ/PtmdnHAIxH+vX4POw/kcOfZp2OmUetEROTU+ZPg78WX5PPMLNv7OhzowJoa5xzPL9pGXEwk4wd2CnY4IiLSwPkz0E1UXQTS1C375gBr0w7xxwmDCdWY8yIiUkP+PIPHzC4DzvJWFzrn3g9cSE3TP/67jZiWzZg4tGuwQxERkUagylv0ZvYYvtv0G72ve83sT4EOrCnZvDebBZv2cdMZcUSEhwY7HBERaQT8acFfDCQ654oBzOwlYDXwy0AG1pS8+HkqzcNCuGFUj2CHIiIijYS/I9m1LbPcJhCBNFWHcgqYszqdKxK7EN2yWbDDERGRRsKfFvyfgNVmtgAwfM/ipwQ0qibk7ZW7yC0o4sYz1HoXEZHa408v+llmthAY7hU94JzbE9ComoiiYsfLX+xgeFw0g07TjREREak9Fd6iN7P+3udQoDOQ5n2d5pVVyswizGy5ma01sw1m9juvvKeZLTOzrWb2ppk188qbe+tbve1xNf/26rf/bs5g54EcbhwdF+xQRESkkamsBf8z4HbgiXK2OeDcKurOB851zh0xs3DgMzP7l1fvk865N8zs78CtwHPe50HnXG8zuw74X+Da6n07DcvMJTvoENWciwZrYBsREaldFSZ459zt3uL3nHN5ZbeZWURVFTvnHHDEWw33vkr+MPihV/4S8DC+BH+5twzwDvCMmZlXT6Ozfd8RFm3ex33n9yU8VLP2iohI7fInsyzxs+w7zCzUzNYAGcDHwDYgyzlX6O2SBnTxlrsAuwC87YeAGH/O0xC9unQn4aHGD0Z2C3YoIiLSCFXYgjezTviSbgszS8LXgx6gNRDpT+XOuSIg0czaArOB/jULF8zsdnyPDujevXtNqwuK/MIi3ludxgUDO9EhqsqbISIiItVW2TP4C4FJQFfgL2XKs4FfVeckzrks7zW70UBbMwvzWuldgXRvt3SgG5BmZmH43rfPLKeu54HnAZKTkxvk7fuPNuwlK6eAa4er9S4iIoFR2TP4l4CXzGyic+7d6lZsZrFAgZfcWwDj8XWcWwBcBbwB3ATM9Q6Z561/4W3/tLE+f3/jy510aduCM3u3D3YoIiLSSPnzHvy7ZnYJMAiIKFP+SBWHdsb3B0Iovmf9bznn3jezjcAbZvYHfEPeTvf2nw68YmZbgQPAddX+bhqAnZk5fL41k5+N70uIZo0TEZEAqTLBe6+yRQLjgBfwta6XV3Wcc+4rIKmc8u3AiHLK84Crqw65YXtrxS5CDK5O1qxxIiISOP70oj/DOXcjvnfUf4fvOXrfwIbVOBUWFfP2yl2c068Dndu0CHY4IiLSiPmT4Evegc8xs9OAAny336WaFm7ax97D+VynznUiIhJg/kw280/vNbepwCp8g9VMC2hUjdRbK3YRG9Wccf07BDsUERFp5CpN8GYWAsx3zmUB75rZ+0CEc+5QnUTXiBw8eowFmzKYdEacRq4TEZGAqzTTOOeKgWfLrOcruZ+a99ftpqDIMSFJnetERCTw/GlKzjeziWamd7pqYPaqNPp1jGJA56hghyIiIk2APwn+DuBtIN/MDptZtpkdDnBcjUrq/qOs2pnFhKFd0N9JIiJSF/wZ6EZNzhqavTodM7g88bRghyIiIk1ElS14M5vvT5mUzznHnDXpnNErRu++i4hInalsNrkIfCPYtTezaE6cTa5LRcfJiVbtPMiOzBzuPrdPsEMREZEmpLJb9HcAPwVOA1ZyPMEfBp4JcFyNxnur0okID+GiwZ2CHYqIiDQhlc0m9xTwlJnd7Zz7Wx3G1GgUFhXz4brdjB/YiVbN/RlTSEREpHb408nub2Z2BhBXdn/n3MsBjKtR+GJ7JgdzCrh0iEb2FRGRuuXPbHKvAL2ANUCRV+wAJfgqvL92Ny2bhXJ239hghyIiIk2MP/eNk4GBzjkX6GAak4KiYv69YQ/jB3YkIjw02OGIiEgT489AN+sB9RCrps+37udQbgGXDNG77yIiUvf8acG3Bzaa2XIgv6TQOXdZwKJqBD74ajdRzcMY26d9sEMREZEmyJ8E/3Cgg2hsjhUW8x/dnhcRkSDypxf9f82sB9DHOfeJmUUCylqV+Hzrfg7nFXKJes+LiEiQ+DNU7W3AO8A/vKIuwJxABtXQvf/VbqIiwjhTt+dFRCRI/OlkdxcwBt8IdjjntgAdAhlUQ1ZQVMzHG32355uH6UaHiIgEhz8JPt85d6xkxczC8L0HL+VYtv0Ah/MKuXCQXjwQEZHg8SfB/9fMfgW0MLPx+OaG/2dgw2q4Pt64h4jwEM7qo8FtREQkePxJ8FOAfcA6fBPQfAj8JpBBNVTOOT7auJexfWJp0Uy350VEJHj8eU2uBTDDOTcNwMxCvbKcQAbWEK1PP8zuQ3n8bHzfYIciIiJNnD8t+Pn4EnqJFsAngQmnYft44x5CDM4b0DHYoYiISBPnT4KPcM4dKVnxliMDF1LD9dHGvSTHtaNdy2bBDkVERJo4fxL8UTMbWrJiZsOA3MCF1DDtzMwhZU82FwxU611ERILPn2fwPwXeNrNvAcM38cy1AY2qAfpo4x4ALhio1+NERCT4/Bmq9ksz6w/084o2OecKAhtWw/PRxr307xRF9xg9vRARkeDzpwUPMByI8/YfamY4514OWFQNTFbOMVakHuCucb2DHYqIiAjgR4I3s1eAXsAaoMgrdoASvGfRlv0UOxjXXyP4iohI/eBPCz4ZGOic0/C0FViYkkF0ZDgJXdsGOxQRERHAv1706/F1rJNyFBc7/rt5H2f3jSU0xIIdjoiICOBfC749sNHMlgP5JYXOucsCFlUDsi79EJlHj3FOP92eFxGR+sOfBP9woINoyBZsysAMzuqryWVERKT+8Oc1uf+aWUd8PekBljvnMgIbVsOxYNM+Eru11eh1IiJSr1T5DN7MrgGWA1cD1wDLzOyqQAfWEGQeyeertCzG6fa8iIjUM/7cov81MLyk1W5msfgmm3knkIE1BIu27MM5OKefbs+LiEj94k8v+pCTbsln+nOcmXUzswVmttHMNpjZvV55OzP72My2eJ/RXrmZ2dNmttXMvio7/n19tSBlH+1bNWPwaW2CHYqIiMgJ/Enw/zaz/5jZJDObBHwA/MuP4wqB/3HODQRGAXeZ2UBgCjDfOdcH31S0U7z9vwf08b5uB56r1ndSx4qKHYu27OPsvh0I0etxIiJSz/jTye5+M7sSONMret45N9uP43YDu73lbDP7GugCXA6c4+32ErAQeMArf9kbUGepmbU1s85ePfXOV2lZZOUUcLZuz4uISD1UYYI3s95AR+fc586594D3vPIzzayXc26bvycxszggCVjm1VmStPcAJfOrdgF2lTkszSs7IcGb2e34Wvh0797d3xBq3edb9wMwpldM0GIQERGpSGW36P8KHC6n/JC3zS9m1gp4F/ipc+6E+rzWerWGwHXOPe+cS3bOJcfGBq/1vHjLfgad1pqYVs2DFoOIiEhFKkvwHZ1z604u9Mri/KnczMLxJffXvLsAAHvNrLO3vTNQ0oEvHehW5vCuXlm9k3OskFU7D3Jm7/bBDkVERKRclSX4ymZOaVFVxWZmwHTga+fcX8psmgfc5C3fBMwtU36j15t+FHCovj5/X/7NAQqKHGOU4EVEpJ6qLMGvMLPbTi40s8nASj/qHgPcAJxrZmu8r4uBx4DxZrYFON9bB/gQ2A5sBaYBP/b/26hbn23ZT7OwEEb0bBfsUERERMpVWS/6nwKzzex6jif0ZKAZMKGqip1znwEVvT92Xjn7O+CuquqtDz7bup/kHtFEhIcGOxQREZFyVZjgnXN7gTPMbBww2Cv+wDn3aZ1EVk/ty84nZU82v7ioX7BDERERqZA/78EvABbUQSwNwpJtvtfj1MFORETqM39GspMyFm/ZT5sW4QzS8LQiIlKPKcFXg3OOz7fuZ0zvGEI1PK2IiNRjSvDVsH3/UXYfytPrcSIiUu8pwVfD0u2ZAIw+XcPTiohI/aYEXw3Lth+gQ1RzerZvGexQREREKqUE7yfnHEu3ZzLy9Bh8g/SJiIjUX0rwfkrNzCEjO5+RGr1OREQaACV4Py3znr+P0vN3ERFpAJTg/bTsmwO0b9WcXrF6/i4iIvWfErwfSp+/92yn5+8iItIgKMH7YdeBXHYfymPU6Xr+LiIiDYMSvB+WfuN7/j5Sz99FRKSBUIL3w9LtmbRr2Yw+HVoFOxQRERG/KMH7Ydn2A4yI0/N3ERFpOJTgq5B2MIf0rFw9fxcRkQalyvngm7ovUw8AMKKnnr+LSOAVFBSQlpZGXl5esEORIIqIiKBr166Eh4efch1K8FVYkXqQVs3D6NcpKtihiEgTkJaWRlRUFHFxcXos2EQ558jMzCQtLY2ePXuecj26RV+FlTsOktS9reZ/F5E6kZeXR0yM5rxoysyMmJiYGt/FUYKvxOG8AjbtzWZYj+hghyIiTYiSu9TGvwEl+Eqs3pmFc5DcQx3sRKTpCA0NJTExkcGDB/P973+frKysU64rLi6O/fv3f6f8yJEj3HHHHfTq1Ythw4ZxzjnnsGzZspqEXWMzZ87k22+/rfZxc+bMYePGjaXrDz74IJ988klthnZKlOArsTL1ACEGid3bBjsUEZE606JFC9asWcP69etp164dzz77bK2fY/LkybRr144tW7awcuVKXnzxxXL/EKhLlSX4oqKiCo87OcE/8sgjnH/++bUeX3UpwVdi5c6D9O/UmlbN1RdRRJqm0aNHk56eXro+depUhg8fzpAhQ3jooYdKy6+44gqGDRvGoEGDeP755yutc9u2bSxbtow//OEPhIT40lDPnj255JJLAPjLX/7C4MGDGTx4MH/9618BSE1NZcCAAdx2220MGjSICy64gNzcXAC2bt3K+eefT0JCAkOHDmXbtm0VxlpRPe+88w4rVqzg+uuvJzExkdzcXOLi4njggQcYOnQob7/9NtOmTWP48OEkJCQwceJEcnJyWLIIRHRMAAASuElEQVRkCfPmzeP+++8nMTGRbdu2MWnSJN555x0A5s+fT1JSEvHx8dxyyy3k5+cDvjsbDz30EEOHDiU+Pp6UlJQa/786mTJXBQqLilm9M4urhnUNdigi0kT97p8b2Pjt4Vqtc+BprXno+4P82reoqIj58+dz6623AvDRRx+xZcsWli9fjnOOyy67jEWLFnHWWWcxY8YM2rVrR25uLsOHD2fixInExJT/evGGDRtITEwkNDT0O9tKWvPLli3DOcfIkSM5++yziY6OZsuWLcyaNYtp06ZxzTXX8O677/KjH/2I66+/nilTpjBhwgTy8vIoLi6uMNbu3btXWM8zzzzD448/TnJycmk8MTExrFq1CoDMzExuu+02AH7zm98wffp07r77bi677DIuvfRSrrrqqhO+l7y8PCZNmsT8+fPp27cvN954I8899xw//elPAWjfvj2rVq3i//7v/3j88cd54YUX/Pr/4i+14CuQsiebnGNF6mAnIk1Obm4uiYmJdOrUib179zJ+/HjAl+A/+ugjkpKSGDp0KCkpKWzZsgWAp59+moSEBEaNGsWuXbtKy6vrs88+Y8KECbRs2ZJWrVpx5ZVXsnjxYsDXyk9MTARg2LBhpKamkp2dTXp6OhMmTAB8749HRkZWGmt59VTk2muvLV1ev349Y8eOJT4+ntdee40NGzZU+r1s2rSJnj170rdvXwBuuukmFi1aVLr9yiuv9CuGU6UWfAVW7jgIQHKcOtiJSHD429KubSXP4HNycrjwwgt59tlnueeee3DO8ctf/pI77rjjhP0XLlzIJ598whdffEFkZCTnnHNOpa94DRo0iLVr11JUVFRuK74izZs3L10ODQ0tvUVfnopiTU1NrVY9LVu2LF2eNGkSc+bMISEhgZkzZ7Jw4UK/Yy9PSRyhoaEUFhbWqK7yqAVfgRU7DtKpdQSntYkIdigiIkERGRnJ008/zRNPPEFhYSEXXnghM2bM4MiRIwCkp6eTkZHBoUOHiI6OJjIykpSUFJYuXVppvb169SI5OZmHHnoI5xzgS7wffPABY8eOZc6cOeTk5HD06FFmz57N2LFjK6wrKiqKrl27MmfOHADy8/NL/zApL9bKREVFkZ2dXeH27OxsOnfuTEFBAa+99lqVx/Xr14/U1FS2bt0KwCuvvMLZZ59daQy1SQm+AitTDzAsLlrvo4pIk5aUlMSQIUOYNWsWF1xwAT/84Q8ZPXo08fHxXHXVVWRnZ3PRRRdRWFjIgAEDmDJlCqNGjaqy3hdeeIG9e/fSu3dvBg8ezKRJk+jQoQNDhw5l0qRJjBgxgpEjRzJ58mSSkpIqreuVV17h6aefZsiQIZxxxhns2bOnwlgrM2nSJO68887STnYn+/3vf8/IkSMZM2YM/fv3Ly2/7rrrmDp1KklJSaUd/MD3uODFF1/k6quvJj4+npCQEO68884qr01tsZK/nhqi5ORkt2LFilqv99usXM547FMe+v5Abh5z6sMEiohU19dff82AAQOCHYbUA+X9WzCzlc655AoOOYFa8OVYtdP3/F0d7EREpKFSgi/Hmp1ZNAsLYUDn1sEORURE5JQowZdjbVoWg09rTXioLo+IiDRMymAnKSgqZl36IRK76fa8iIg0XErwJ9m8N5u8gmISurUJdigiIiKnTAn+JGt2+WZNSuymCWZERKThUoI/ydpdWURHhtO9XWSwQxERCZo5c+ZgZgGZBMUfqampvP7661Xut2bNGj788MPS9Xnz5vHYY48FLK6Kpr+tj5TgT7J21yESurXVADci0qTNmjWLM888k1mzZgXl/Kea4C+77DKmTJkSyNAajIAleDObYWYZZra+TFk7M/vYzLZ4n9FeuZnZ02a21cy+MrOhgYqrMkfyC9mcka3b8yLSpB05coTPPvuM6dOn88Ybb5SW7969m7POOovExEQGDx7M4sWLmTFjRunsaADTpk3jvvvuIzU1lf79+zNp0iT69u3L9ddfzyeffMKYMWPo06cPy5cvB+Dhhx/mhhtuYPTo0fTp04dp06YBMGXKFBYvXkxiYiJPPvkkeXl53HzzzcTHx5OUlMSCBQs4duwYDz74IG+++SaJiYm8+eabzJw5k5/85CcA7N27lwkTJpCQkEBCQgJLliw54fv8+9//zv3331+6XvbYqqa/TU1NZfDgwaXrjz/+OA8//DDgmw73oosuYtiwYYwdOzZod0ECOdnMTOAZ4OUyZVOA+c65x8xsirf+APA9oI/3NRJ4zvusU+vSDuEcJCjBi0h98K8psGdd7dbZKR6+V/kt7Llz53LRRRfRt29fYmJiWLlyJcOGDeP111/nwgsv5Ne//jVFRUXk5OSQlJTEH//4R6ZOnUp4eDgvvvgi//jHPwDfPO1vv/02M2bMYPjw4bz++ut89tlnzJs3j0cffbR0/PivvvqKpUuXcvToUZKSkrjkkkt47LHHePzxx3n//fcBeOKJJzAz1q1bR0pKChdccAGbN2/mkUceYcWKFTzzzDOAL0mXuOeeezj77LOZPXs2RUVFpePSl5g4cSKjR49m6tSpALz55pv8+te/BqjW9Lcnu/322/n73/9Onz59WLZsGT/+8Y/59NNP/Tq2NgUswTvnFplZ3EnFlwPneMsvAQvxJfjLgZedb9zcpWbW1sw6O+d2Byq+8pR2sOuqBC8iTdesWbO49957Ad8467NmzWLYsGEMHz6cW265hYKCAq644orSKVfPPfdc3n//fQYMGEBBQQHx8fGkpqbSs2dP4uPjAd8Mcueddx5mVrq9xOWXX06LFi1o0aIF48aNY/ny5bRte+Lv4c8++4y7774bgP79+9OjRw82b95c6ffx6aef8vLLvjZmaGgobdqc+HZUbGwsp59+OkuXLqVPnz6kpKQwZswYwDf97ezZswFKp7/1J8EfOXKEJUuWcPXVV5eW5efnV3lcINT1dLEdyyTtPUBHb7kLsKvMfmleWZ0m+LW7sugRE0l0y2Z1eVoRkfJV0dIOhAMHDvDpp5+ybt06zIyioiLMjKlTp3LWWWexaNEiPvjgAyZNmsTPfvYzbrzxRiZPnsyjjz5K//79ufnmm0vrKjsta0hISOl6SEjICdOjntznqS77QF133XW89dZb9O/fnwkTJmBmfk1/GxYWRnFxcel6yfbi4mLatm3LmjVr6ux7qEjQOtl5rfVqz3RjZreb2QozW7Fv375ajWnNriwS1HoXkSbsnXfe4YYbbmDHjh2kpqaya9cuevbsyeLFi9mxYwcdO3bktttuY/LkyaxatQqAkSNHsmvXLl5//XV+8IMfVPucc+fOJS8vj8zMTBYuXMjw4cO/MwXr2LFjS6do3bx5Mzt37qRfv36VTvF63nnn8dxzzwFQVFTEoUOHvrPPhAkTmDt3LrNmzeK6664D8Gv6244dO5KRkUFmZib5+fmljxJat25Nz549efvttwHfvPRr166t9jWpDXWd4PeaWWcA77Nkct50oFuZ/bp6Zd/hnHveOZfsnEuOjY2ttcD2HMpjz+E8dbATkSZt1qxZTJgw4YSyiRMnMmvWLBYuXEhCQgJJSUm8+eabpbfxAa655hrGjBlDdHT1RwEdMmQI48aNY9SoUfz2t7/ltNNOY8iQIYSGhpKQkMCTTz7Jj3/8Y4qLi4mPj+faa69l5syZNG/enHHjxrFx48bSTnZlPfXUUyxYsID4+HiGDRvGxo0bv3Pu6OhoBgwYwI4dOxgxYgSAX9PfhoeH8+CDDzJixAjGjx9/wvSxr732GtOnTychIYFBgwYxd+7cal+T2hDQ6WK9Z/DvO+cGe+tTgcwynezaOed+YWaXAD8BLsbXue5p59yIquqvzeli/71+D3e+upJ3/98ZmkVORIKmoU4Xe+mll3Lfffdx3nnnVeu4hx9+mFatWvHzn/88QJE1XPV2ulgzmwV8AfQzszQzuxV4DBhvZluA8711gA+B7cBWYBrw40DFVZGoiDDO7d+BQadpBjkREX9lZWXRt29fWrRoUe3kLoEV0BZ8oNVmC15EpD5oqC14qX31tgUvIiIiwaMELyJSzzTkO6tSO2rj34ASvIhIPRIREUFmZqaSfBPmnCMzM5OIiIga1VPXA92IiEglunbtSlpaGrU9zoc0LBEREXTt2rVGdSjBi4jUI+Hh4fTs2TPYYUgjoFv0IiIijZASvIiISCOkBC8iItIINeiBbsxsH7CjFqtsD+yvxfqaKl3HmtM1rDldw5rTNay52r6GPZxzfk3E0qATfG0zsxX+jhAkFdN1rDldw5rTNaw5XcOaC+Y11C16ERGRRkgJXkREpBFSgj/R88EOoJHQdaw5XcOa0zWsOV3DmgvaNdQzeBERkUZILXgREZFGSAneY2YXmdkmM9tqZlOCHU9DYGYzzCzDzNaXKWtnZh+b2RbvMzqYMdZ3ZtbNzBaY2UYz22Bm93rluo5+MrMIM1tuZmu9a/g7r7ynmS3zfqbfNLNmwY61vjOzUDNbbWbve+u6htVkZqlmts7M1pjZCq8sKD/PSvD4/lEDzwLfAwYCPzCzgcGNqkGYCVx0UtkUYL5zrg8w31uXihUC/+OcGwiMAu7y/u3pOvovHzjXOZcAJAIXmdko4H+BJ51zvYGDwK1BjLGhuBf4usy6ruGpGeecSyzzelxQfp6V4H1GAFudc9udc8eAN4DLgxxTveecWwQcOKn4cuAlb/kl4Io6DaqBcc7tds6t8paz8f1y7YKuo9+czxFvNdz7csC5wDteua5hFcysK3AJ8IK3buga1pag/Dwrwft0AXaVWU/zyqT6OjrndnvLe4COwQymITGzOCAJWIauY7V4t5bXABnAx8A2IMs5V+jtop/pqv0V+AVQ7K3HoGt4KhzwkZmtNLPbvbKg/DxrulgJGOecMzO9puEHM2sFvAv81Dl32Nd48tF1rJpzrghINLO2wGygf5BDalDM7FIgwzm30szOCXY8DdyZzrl0M+sAfGxmKWU31uXPs1rwPulAtzLrXb0yqb69ZtYZwPvMCHI89Z6ZheNL7q85597zinUdT4FzLgtYAIwG2ppZSSNGP9OVGwNcZmap+B5Rngs8ha5htTnn0r3PDHx/bI4gSD/PSvA+XwJ9vB6jzYDrgHlBjqmhmgfc5C3fBMwNYiz1nvecczrwtXPuL2U26Tr6ycxivZY7ZtYCGI+vL8MC4CpvN13DSjjnfumc6+qci8P3++9T59z16BpWi5m1NLOokmXgAmA9Qfp51kA3HjO7GN8zqFBghnPuj0EOqd4zs1nAOfhmS9oLPATMAd4CuuOb6e8a59zJHfHEY2ZnAouBdRx/9vkrfM/hdR39YGZD8HVcCsXXaHnLOfeImZ2OrzXaDlgN/Mg5lx+8SBsG7xb9z51zl+oaVo93vWZ7q2HA6865P5pZDEH4eVaCFxERaYR0i15ERKQRUoIXERFphJTgRUREGiEleBERkUZICV5ERKQRUoIXaaDMLMabsWqNme0xs/Qy60sCdM4kM5seiLr9OPcnmlVPxH96TU6kETCzh4EjzrnHA3yet4E/OOfWBqj+sDJjn5+87Sagq8aoEPGPWvAijZCZHfE+zzGz/5rZXDPbbmaPmdn13vzp68ysl7dfrJm9a2Zfel9jyqkzChjinFtrZiHe3Nax3rYQb87w2IrqMrMRZvaFN9/4EjPr55VPMrN5ZvYpMN/MOpvZIu9OxHozG+uFMA/4QeCvnkjjoAQv0vglAHcCA4AbgL7OuRH4pgW929vnKXzzfg8HJnrbTpaMb9hNnHPFwKvA9d6284G1zrl9ldSVAox1ziUBDwKPlql7KHCVc+5s4IfAf5xziV7sa7xzHgSae6OCiUgVNJucSOP3ZclUlWa2DfjIK18HjPOWzwcGlpnFrrWZtSozzzpAZ2BfmfUZ+MbU/itwC/BiZXUBbYCXzKwPvik1w8vU9XGZoTu/BGZ4k/DMcc6tKbNfBnAakFmN71+kSVKCF2n8yo4dXlxmvZjjvwNCgFHOubxK6skFIkpWnHO7zGyvmZ2Lb8asktZ8uXWZ2TPAAufcBDOLAxaW2Xy0TL2LzOws4BJgppn9xTn3src5wotDRKqgW/QiAr5WfcnteswssZx9vgZ6n1T2Ar5b9W97c7JXVlcbjk83OqmiQMysB7DXOTfNq3+oV25AJyDVn29IpKlTghcRgHuAZDP7ysw24ntmfwLnXArQpmQ6TM88oBXHb89XVtefgT+Z2Woqv3t4DrDW2+9afM/0AYYBSyvqZS8iJ9JrciLiNzO7D8h2zr3grSfj61A3tvIja+XcTwHznHPzA30ukcZALXgRqY7n8J7hm9kU4F3gl3V07vVK7iL+UwteRESkEVILXkREpBFSghcREWmElOBFREQaISV4ERGRRkgJXkREpBFSghcREWmE/j+crSKXKdoo8QAAAABJRU5ErkJggg==\n",
"text/plain": [
"