DAE simulation python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3d23b88f",
   "metadata": {},
   "source": [
    "# Simulating dynamic models in Python\n",
    "\n",
    "In general, many dynamic models encountered in engineering can be discretised into a system of differential algebraic equations (DAE) of the form:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "\\frac{d \\textbf{x}}{d t} & =\\textbf{f}(\\textbf{x}, \\textbf{y}, t) \\\\\n",
    "0 & =\\textbf{g}(\\textbf{y}, \\textbf{x},  t)\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "\n",
    "with the initial conditions $\\textbf{x}(t_0) = x_0$ and $\\textbf{y}(t_0) = y_0$. In general this system can be simulated with a `for` or `while` loop with the following three steps each loop:\n",
    "\n",
    "1. Solve the system of algebraic equations $0 =g(\\textbf{y}, \\textbf{x}, t)$ for $\\textbf{y}$ at the current timestep $t$ using $\\textbf{x}=\\textbf{x}(t)$.\n",
    "1. Compute $\\textbf{f}(\\textbf{x}, \\textbf{y}, t)$.\n",
    "1. Add the solution $\\textbf{f}$ to an _integrator_ (e.g. Euler, Runge-Kutta family of methods etc.) to advance the system for a timestep $\\Delta t$.\n",
    "\n",
    "When the timestep $\\Delta t$ is small enough many dynamic simulations without singularities will converge using simple numerical methods. \n",
    "\n",
    "The code block below is a template of the simplest python loop using a forward Euler method: \n",
    "```python\n",
    "t_f = 10  # Simulate for 10 seconds of time\n",
    "dt = 0.1  # The timestep in seconds\n",
    "t = 0  # Set initial time equal to 0 s\n",
    "x = [1.5, 2, ..., 1]  # Set the initial x vector equal to x_0\n",
    "y = [1, 0.5, ..., 0]  # Set the initial y vector equal to y_0\n",
    "x_sol = []  # List of solutions at each time step\n",
    "y_sol = []  # Usually this is not needed as it can be reconstructed from x, t\n",
    "t_sol = []  # List of time steps at each solution\n",
    "\n",
    "# Start main while/for loop:\n",
    "while t < t_f:\n",
    "    # 1. (solver) Solve algebraic equations 0 = g(y, x, t) \n",
    "    # using analytical solutions, iterative methods or a solver \n",
    "    # such as scipy.optimize.fsolve\n",
    "    y = ...  # New solution vector\n",
    "    # 2. Compute f(x, y, t):\n",
    "    f = ...  # solve the system with current x, y, t\n",
    "    dxdt = f  # \n",
    "    # 3. (integrator) Advance the time of the (Eulerian step):\n",
    "    x = x + dxdt * dt  # Advance the system using forward Euler step of the simulation\n",
    "    t += dt  # Advance the time\n",
    "    # Save the solutions:\n",
    "    x_sol.append(x), y_sol.append(y), t_sol.append(t)\n",
    "```\n",
    "\n",
    "In steady-state systems the independent variable $t$ is often replaced with another variable such as reactor volume $V$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "618b5492",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Basic imports \n",
    "import numpy as np\n",
    "import scipy.optimize\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "342063ee",
   "metadata": {},
   "source": [
    "## Example tank flow system\n",
    "\n",
    "We will use the following example from [AP Monitor](http://apmonitor.com/che263/index.php/Main/PythonDynamicSim) in order to demonstrate the three steps outlined above:\n",
    "\n",
    "Example 1: A system of two tanks have the following differential equations that relates the inlet flow to the height of the tanks:\n",
    "$$\n",
    "A_c \\frac{d h_1}{d t}=q_{\\text {in }}-q_{\\text {out } 1} \\quad A_c \\frac{d h_2}{d t}=q_{\\text {out } 1}-q_{\\text {out } 2}\n",
    "$$\n",
    "The system has two algebraic equations for the height dependent outlet flow rate according to Bernoulli's equation for each tank:\n",
    "\n",
    "$$q_{\\text {out } 1}=c_1 \\sqrt{h_1} \\quad q_{\\text {out2 }}=c_2 \\sqrt{h_2}$$\n",
    "\n",
    "Both tanks are initially empty when the inlet to the first starts to flow at a rate of $0.5 \\mathrm{~m}^3 / \\mathrm{h}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efcf1355",
   "metadata": {},
   "source": [
    "<img src=\"ant_ii_ex_.svg\"  width=\"240\">"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bac1173",
   "metadata": {},
   "source": [
    "###  Restate the problem in standard format\n",
    "\n",
    "\n",
    "First we start with the algebraic equations:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "q_{\\text {out } 1}=c_1 \\sqrt{h_1} \\\\\n",
    "q_{\\text {out2 }}=c_2 \\sqrt{h_2}\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "\n",
    "We need to put this into the format\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "0 & =\\textbf{g}(\\textbf{y}, \\textbf{x}, t)\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "\n",
    "First we replace the notation for the differential variables which are $x_1 = h_1$ and $x_2 = h_2$ and the algebriac variables  $y_1 = q_{\\text {out } 1}$ and $y_2 = q_{\\text {out } 2}$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "y_1=c_1 \\sqrt{x_1} \\\\\n",
    "y_2=c_2 \\sqrt{x_2}\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "\n",
    "Note that standard solvers usually need the algebraic equations to zero:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "0 = c_1 \\sqrt{x_1} - y_1 \\\\\n",
    "0 = c_2 \\sqrt{x_2} - y_2\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8c74db7",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "\n",
    "Next, we transform the differential equations into our notation:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "A_c \\frac{d h_1}{d t}=q_{\\text {in }}-q_{\\text {out } 1} \\\\\n",
    " A_c \\frac{d h_2}{d t}=q_{\\text {out } 1}-q_{\\text {out } 2}\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "\n",
    "becomes:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    " \\frac{d x_1}{d t}=\\frac{q_{\\text {in }}-y_1}{A_c} \\\\\n",
    " \\frac{d x_2}{d t}=\\frac{y_1 - y_2}{A_c}\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "\n",
    "Where $q_{\\text {in} }$ is called an input to the system and $A_c$, $c_1$ and $c_2$ are constant parameters. The overall system is then \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6805a424",
   "metadata": {},
   "source": [
    "\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    " \\frac{d x_1}{d t} & = f_1(\\textbf{x}, \\textbf{y}, t) = \\frac{q_{\\text {in }}-y_1}{A_c} \\\\\n",
    " \\frac{d x_2}{d t} & = f_2(\\textbf{x}, \\textbf{y}, t) = \\frac{y_1 - y_2}{A_c} \\\\\n",
    "0 & = g_1(\\textbf{x}, \\textbf{y}, t) = c_1 \\sqrt{x_1} - y_1 \\\\\n",
    "0 & = g_2(\\textbf{x}, \\textbf{y}, t) = c_2 \\sqrt{x_2} - y_2\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "\n",
    "with initial conditions:\n",
    "$$ x_1 (t=0) = x_2(t=0) = y_1(t=0) = y_2(t=0) = 0$$\n",
    "\n",
    "In our example, we will make the tank inflow time depedent so that:\n",
    "\n",
    "$$\n",
    "q_{in} = \n",
    "     \\begin{cases}\n",
    "       &   0.5  &, t \\leq 4 \\\\\n",
    "                 &  0  &, t > 4\n",
    "     \\end{cases}\n",
    "$$\n",
    "\n",
    "In code this is simply:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "dbdc994e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def q_in(t):  # Flow rate into first tank m^3 /h\n",
    "    if t <= 4:\n",
    "        return 0.5  # m^3/h\n",
    "    elif t > 4:\n",
    "        return 0.0  # m^3/h\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f784e8e",
   "metadata": {},
   "source": [
    "## 1. Solving algrebraic equations\n",
    "\n",
    "While the equations in this problem statement have exact analytical solutions, we will show how to solve any general system of equations when it is written in standard form $0  =\\textbf{g}(\\textbf{y}, \\textbf{x},  t)$ using a numerical equation solver like [scipy.optimize.fsolve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html):\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "51609a8c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0.])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Define the parameters \n",
    "Ac = 2  # m^2, cross sectional area of the tanks\n",
    "c1 = 0.13  # -, valve constant 1\n",
    "c2 = 0.20  # -, valve constant 2\n",
    "\n",
    "# Define a function that returns the system of algebraic equations:\n",
    "def g_sys(y, x, t):\n",
    "    # Return the vector [g1, ..., gn]\n",
    "    g1 = c1 * np.sqrt(x[0]) - y[0]  # = c_1 sqrt(x_1) - y_1\n",
    "    g2 = c2 * np.sqrt(x[1]) - y[1]  # = c_2 sqrt(x_2) - y_2\n",
    "    return [g1, g2]\n",
    "\n",
    "# Test the solver for the initial conditions:\n",
    "x0 = [0, 0]  # [x_1 (t_0), x_2 (t_0)]\n",
    "y0 = [0, 0]  # [y_1 (t_0), y_2 (t_0)]\n",
    "t0 = 0  # h, starting time\n",
    "\n",
    "scipy.optimize.fsolve(g_sys, y0, args=(x0, t0))  # Solves g_sys, outputs solution y\n",
    "# It is also worth testing g_sys(x0, y0, t0) while writing g_sys"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bf0f544",
   "metadata": {},
   "source": [
    "Next we wrap the entire solution method into a function that computes the solution to simplify the later loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "da3b4d5f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0.])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def g_sol(y, x, t):\n",
    "    # Note: Instead of returning the numerical solution,\n",
    "    #       the analytical can also be set here.\n",
    "    sol = scipy.optimize.fsolve(g_sys, y, args=(x, t))  # = [y1, y2]\n",
    "    return sol\n",
    "\n",
    "y = g_sol(y0, x0, t0)\n",
    "y"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d1bc2ea",
   "metadata": {},
   "source": [
    "## 2. Compute $\\frac{d \\textbf{x}}{d t}  =\\textbf{f}(\\textbf{x}, \\textbf{y}, t) $\n",
    "\n",
    "The next step is to write a function that returns the right-hand side of the system $\\frac{d \\textbf{x}}{d t}  =\\textbf{f}(\\textbf{x}, \\textbf{y}, t) $. This can then be used to advance the system forward in time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8c3e4fe5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.25, 0.  ])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def f_sys(x, y, t):\n",
    "    f1 = (q_in(t) - y[0]) / Ac\n",
    "    f2 = (y[0] - y[1]) / Ac\n",
    "    f = np.array([f1, f2])  # Wrap into array data type (vector)\n",
    "    return f  # Return vector of solutions f\n",
    "\n",
    "# Test for initial conditions (including y computed in the previous cell):\n",
    "dxdt = f_sys(y, x0, t0)\n",
    "dxdt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dc1b2b9",
   "metadata": {},
   "source": [
    "This is a good time to perform some basic sanity checks when possible to see if the system makes physical sense. For this example we expect the height of the first tank to increase while the flow out of the second tank should be zero (since the initial height of tank 1 is zero and there can be no flow out) as the output shows.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e2a8dea",
   "metadata": {},
   "source": [
    "## 3. Integrator\n",
    "\n",
    "Finally we add the functions to our integrator. First we will use a simple forward step Euler method `while` loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "272e8a65",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the parameters \n",
    "Ac = 2  # m^2, cross sectional area of the tanks\n",
    "c1 = 0.13  # -, valve constant 1\n",
    "c2 = 0.20  # -, valve constant 2\n",
    "\n",
    "# Define the numerical simulation parameters\n",
    "t = 0  # Set initial time equal to 0 h\n",
    "t_f = 10  # Simulate for 10 hours of time\n",
    "dt = 0.1  # The timestep in hours\n",
    "\n",
    "# Set the initial conditions\n",
    "x = [0, 0]  # Set the initial x vector equal to [x_1(t=0), x_2(t=0)]\n",
    "y = [0, 0]  # Set the initial y vector equal to [y_1(t=0), y_2(t=0)]\n",
    "\n",
    "# Define data containers to save each time step\n",
    "x_sol = []  # List of solutions at each time step\n",
    "y_sol = []  # Usually this is not needed as it can be reconstructed from x, t\n",
    "t_sol = []  # List of time steps at each solution\n",
    "\n",
    "# Start main for loop:\n",
    "while t < t_f:\n",
    "    # 1. (solver) Solve algebraic equations 0 = g(y, x, t) \n",
    "    # using analytical solutions, iterative methods or a solver \n",
    "    # such as scipy.optimize.fsolve\n",
    "    y = g_sol(y, x, t)  # New solution vector\n",
    "    # 2. Compute f(x, y, t):\n",
    "    dxdt = f_sys(x, y, t)  # solve the system dxdt = f(x, y, t) with current x, y, t\n",
    "    # 3. (integrator) Advance the time of the (Eulerian step):\n",
    "    x = x + dxdt * dt  # Advance the system using forward Euler step of the simulation\n",
    "    t += dt  # Advance the time\n",
    "    # Save the solutions:\n",
    "    x_sol.append(x), y_sol.append(y), t_sol.append(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5702f7b0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.025     , 0.04897226, 0.07253383, 0.09578324, 0.11877157,\n",
       "       0.14153146, 0.16408611, 0.18645312, 0.20864641, 0.23067735,\n",
       "       0.25255547, 0.2742889 , 0.29588468, 0.31734899, 0.3386873 ,\n",
       "       0.3599045 , 0.38100502, 0.40199286, 0.42287167, 0.44364481,\n",
       "       0.46431538, 0.48488623, 0.50536004, 0.52573927, 0.54602626,\n",
       "       0.56622318, 0.58633207, 0.60635487, 0.6262934 , 0.64614938,\n",
       "       0.66592446, 0.68562019, 0.70523804, 0.72477945, 0.74424574,\n",
       "       0.76363821, 0.78295809, 0.80220657, 0.82138478, 0.84049381,\n",
       "       0.83453472, 0.82859678, 0.82268   , 0.81678439, 0.81090995,\n",
       "       0.80505666, 0.79922454, 0.79341358, 0.78762379, 0.78185515,\n",
       "       0.77610769, 0.77038138, 0.76467624, 0.75899227, 0.75332946,\n",
       "       0.74768781, 0.74206733, 0.73646802, 0.73088986, 0.72533288,\n",
       "       0.71979706, 0.7142824 , 0.70878891, 0.70331659, 0.69786543,\n",
       "       0.69243544, 0.68702661, 0.68163895, 0.67627246, 0.67092713,\n",
       "       0.66560297, 0.66029998, 0.65501816, 0.6497575 , 0.64451801,\n",
       "       0.63929969, 0.63410253, 0.62892655, 0.62377173, 0.61863808,\n",
       "       0.6135256 , 0.60843429, 0.60336414, 0.59831517, 0.59328737,\n",
       "       0.58828073, 0.58329527, 0.57833097, 0.57338785, 0.56846589,\n",
       "       0.56356511, 0.55868549, 0.55382705, 0.54898978, 0.54417368,\n",
       "       0.53937875, 0.534605  , 0.52985241, 0.525121  , 0.52041076,\n",
       "       0.51572169])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Process data\n",
    "x_sol = np.array(x_sol)  # Put list into an array for easier processing\n",
    "y_sol = np.array(y_sol)  # \"\n",
    "x_sol[:, 0]#, y_sol "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "618157ca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAsbklEQVR4nO3deXhU9fXH8fcBBAQRlEWFsAoqSGUxgmt/rojaFttqQa27IioKdcGtdcHdiksVREQU1EqtpS4VN7St1roQ9s0FESEgsosiSEK+vz/OpAkxCUmYmzuT+bye5z6TzEyYE4V77v0u51gIARERyVy14g5ARETipUQgIpLhlAhERDKcEoGISIZTIhARyXB14g6gspo1axbatWsXdxgiImll2rRpq0MIzUt7Le0SQbt27cjJyYk7DBGRtGJmX5b1moaGREQynBKBiEiGUyIQEclwaTdHICIStby8PHJzc9m8eXPcoVRa/fr1ycrKYqeddqrwzygRiIiUkJubS6NGjWjXrh1mFnc4FRZCYM2aNeTm5tK+ffsK/5yGhkRESti8eTNNmzZNqyQAYGY0bdq00ncySgQiIqVItyRQqCpxKxFI2ps5E558EjZtijsSkfSkRCBpbe1aOPFEOPdcaNsWbr0V1qyJOyqR9KJEIGnt8sth1SoYOxZ69YIbb4Q2bfz5xYvjjk4kPSgRSNqaNAmeeQZ+/3s4/3z4xz9gzhw49VQYPRo6doTTToPp0+OOVKRqpkyZwplnnhn55ygRSFpauRIGDYKePeH664ue79rV5wsWLYLf/Q5eeQUOPBCOOQZeew3UmVXSyaxZs+jRo0fkn6NEIGknBLj4YvjmGxg/HkrbN5OVBX/8IyxdCvfcAx9/DCecAAcc4D+zZUv1xy1SWbNmzWLFihUcccQR7LnnnkyZMiWSz9GGMkk7Eyf6sNCdd/odQHkaN4arr4YhQ+DZZ+Hee+Gcc/wuYsgQuOgif49IWYYO9ZVpydS9OzzwwPbfN2vWLLp27cq7777LpEmTeOaZZzj22GOTGwy6I5A089VXcOml0Ls3XHVVxX+ubl04+2yYPRsmT4bOneGaa6B1a7jySliyJLqYRaoiLy+PtWvXclXiL3p+fj5NmjRh0aJFnH/++ZxyyilJ+6xI7wjMrC/wIFAbGBtCuKvE642Bp4E2iVjuDSE8EWVMkr5CgIEDfb/A+PFQpwp/e818iOiEE2DGDBgxAh580I/+/T0p9OyZ/NglfVXkyj0K8+fPp1u3btSq5dfrs2fPpmvXrnTo0IHHH388qYkgsjsCM6sNjAROALoAp5lZlxJvuxSYH0LoBhwJjDCzulHFJOltwgRfGXTHHbDvvjv+5/XoAU8/7RPLQ4fCyy/7xPLRR/skc0HBjn+GSFXNmjWLbt26/e/72bNnc8ABB0TyWVEODfUCFoYQFoUQtgATgX4l3hOARuZ7oncB1gL5EcYkaSo318f0Dz/c9wgkU5s2PnewdKlPMH/2GfzsZz7/8NhjkIYFKKUGmDVr1jYn/rlz59J1e5NiVRRlImgFLC32fW7iueIeBjoDy4E5wJAQwo+uw8xsoJnlmFnOqlWroopXUlQIcOGFkJcHTzwBtWtH8zmNG/u8w6JFfqdQv74PRbVpA7fc4hvXRKrLiBEj6N+///++X7RoETvvvDNr1qxh0KBBzJgxgzvvvDMpnxVlIiit8lHJVdzHAzOBlkB34GEz2/VHPxTCmBBCdgghu3nzUnsvSw02bpzvAbj7bt8kFrWddoIzzoBp0+Dtt33H8s03e0IYOBAWLIg+BpGyNG3alNGjR/P5559z3XXXJeXPjDIR5AKti32fhV/5F3cuMCm4hcAXwH4RxiRpZskS3xh25JFwySXV+9lmcNRRPi+xYAGcdRY89RR06QInnQRvvaUNalIzRJkIpgKdzKx9YgJ4APBSifcsAY4BMLM9gH2BRRHGJGkkBC8dEYLfFdSKcbHzfvvBo496Yrr5ZsjJgWOP9QnnCRO0QU3SW2T/tEII+cBg4HVgAfBcCGGemQ0ys0GJt90KHGpmc4C3gGtCCKujiknSy5gxMGWKT+BWotlSpJo3h5tugi+/hMcfh/x835/Qti3cfjus1t9eSUMW0uzeNjs7O+Tk5MQdhkRs8WL4yU/g4IPhjTd8mCYVhQBvvgn33+/zGPXre2IYOtTvIiQ9LViwgM6dO8cdRpWVFr+ZTQshZJf2fu0slpRTUOBDQmZeXjpVkwB4bH36wKuvwrx5cOaZXvSuc2fvk/DGG5pHkNSnRCApZ/RoX60zYoQPuaSLLl18OGvpUhg+3MtfH3+839mMHasOapK6lAgkpXzxBQwbBscdBxdcEHc0VdO8OfzhDz6P8OSTXgrjwgt9+ekf/uD1kkRSiRKBpIyCAjjvPF8dlOpDQhVRr57PF8yYAf/8Jxx2mE8ot23rS1GnTYs7QhGnRCApY/Ro+Ne/4L77/Oq5pjDzfRAvvACffuq9FP7+d8jOhiOOgL/9zVcficRFiUBSQuGQUJ8+PlFcU3Xs6JVOc3N9DiQ3F045Bfbe2+sdrVsXd4SSStSqUjJG8SGhxx5L/yGhimjcGK64AhYu9LuD9u29gU5Wlu+g/vjjuCOUVKBWlZIxCoeERoyoWUNCFVG7Npx8sv/+M2bAb37jG9U6d4a+fX1ZqsphZy61qpSMUHxIKF1XCSVL9+5eXfXuu72cxahRvhdhn33gsst84rlRo7ijzEAx9qpUq0qp8Qo3jmXSkFBFtGhRtPz0mWdgt908EWRl+Tlp4cK4I5TqUFaryhdeeIELL7yQfv368cYbbyTls3RHILEZM8aXVY4Zk3lDQhVRty6cfrofH34If/oTjBzpjyee6A16jjtOCTRyMfWqLKtV5cknn8zJJ5/MunXruOqqq+jTp88Of5buCCQWixf75Oixx2pIqCJ69/a7gy+/9LuFnBzftdyliyeHb7+NO0JJtu21qrztttu49NJLk/JZSgRS7UIoOvnXhI1j1allS++W9uWX3huhUSMYPLho2Oizz+KOUJKlrFaVIQSuueYaTjjhBHr27JmUz1IikGr32GPe1OWPf0yvWkKppF49+O1v4aOP4P33vcfyqFE+sXziiTB5slYbpbuyWlU+9NBDTJkyheeff57Ro0cn5bNUhlqq1ZIl3hT+oIO8fHOczWZqmhUrfLXR6NH+9d57w6WXwrnnQpMmcUeXXlSGWiQiIXjP34ICHxJSEkiuPfcsaprz7LOwxx6+aa1VKxg0CObOjTtCSVX6pyjV5okn4PXX4a67UqfjWE1Uty4MGADvveeF7QYMgPHjvRz2UUeptpH8mBKBVItly/zq9Kc/rf4m9JmsZ0/fqZybC/fc46u1TjkF2rWD226Dr7+OO0JJBUoEErkQ4KKLvMH7449rSCgOTZv6ct2FC+HFF2H//X0ZauvWvk/hvffUSa2kdJs/LVSVuPVPUiL39NPwyitei79jx7ijyWy1a8MvfuFDdJ984ndnkyfD4Yf73cNjj8HGjXFHGb/69euzZs2atEsGIQTWrFlD/fr1K/VzWjUkkVqxwjc97bcfvPuun4gktWzc6JvVRo6E2bO9Muo553jfhH33jTu6eOTl5ZGbm8vmzZvjDqXS6tevT1ZWFjvttNM2z5e3akiJQCITAvz6137FOXOmJwNJXSH4ENGoUfD885CX5zu/L7kEfv5zb7kp6UvLRyUWzz/vtfZvuUVJIB2Y+RDRn//s+z1uvdWHj371K1/ldeut6rdcU+mOQCKxerUPCbVt6ztfdTWZnvLzfX5n5EjfAFinjieGiy+G//s/lQdJJ7ojkGo3ZAisXw/jxikJpLM6daBfP3jjDe+3fPnlnhCOOspXHj30EHzzTdxRyo5SIpCke+klH1644QbfxCQ1Q6dO3kVu2TLfHNiokSeGli3hwgth+vS4I5Sq0tCQJNX69T4k1Lw5TJ3qu1yl5po2DR55xBP/pk1eQ+rii6F/f2jQIO7opDgNDUm1ueoqWLnSh4SUBGq+Aw/0ulHLl3vDnI0b4bzzvL7R0KGwYEHcEUpFKBFI0rz5pu8cvvpqP0FI5mjSxNtpzp0L//439O3ry1C7dPFJ5WefhR9+iDtKKYuGhiQpvvvOy0vXr+97Biq5sVFqoJUr4cknvTT2okXQrJnfLQwc6CWypXppaEgid911vvb88ceVBMS1aAHDhnnXtNdfhyOO8Mnmjh2913LhpjWJnxKB7LD//AcefthbJh52WNzRSKqpVQv69IFJk/xiYfhw36h26qle9O766+GLL+KOMrNpaEh2yKZN0L27VxadMwd22SXuiCQdbN0Kr73mw0avvOLlLY47zqvU/vznUKJMjiSBhoYkMsOH+0ajMWOUBKTiateGk07yPSeLF8ONN8L8+V6bqk0b34Oiu4Tqo0QgVTZ9ujegP+88v5oTqYrWreHmm/3E//LLkJ3tXew6dIDjj/e5hC1b4o6yZtPQkFRJXp5vHvr6a7+S2223uCOSmiQ31/eijB0LS5f6xPM55/gOZvW0qBoNDUnS/fGPMGuWrxVXEpBky8ry4aIvvvA5hEMP9RVHnTrB0Uf7voQ0bBWQspQIpNI+/thLS596Kvzyl3FHIzVZ7dpw4oleznzJEu9yt3ixt9ds1Qp+9zuYNy/uKNNfpInAzPqa2SdmttDMri3jPUea2Uwzm2dm/44yHtlxBQVwwQXQsKFXnhSpLi1b+lLThQu9Guqxx3p57K5d/Y5h3Di12ayqyBKBmdUGRgInAF2A08ysS4n3NAFGAb8IIewPnBpVPJIco0Z5F6v774c99og7GslEtWr54oS//MUrod57L6xbB+efD3vt5UtQp071JalSMVHeEfQCFoYQFoUQtgATgX4l3nM6MCmEsAQghLAywnhkB335pe8gPv54OOusuKMR8Sq3V17pCxb+8x9vmvPUU9CrF3Tr5oXw1q6NO8rUF2UiaAUsLfZ9buK54vYBdjOzf5nZNDMr9fRiZgPNLMfMclatWhVRuFKeEGDQIH989FF1ppLUYua72p980ttpjh7t1W+HDPEhpdNOgylTfGhTfizKRFDaqaLkzVod4EDgJOB44A9mts+PfiiEMSGE7BBCdvPmzZMfqWzXM8/4TtA77/T2kyKpqnFjHx7KyfECiAMHeq2j447zYnfDh/uSVCkSZSLIBVoX+z4LWF7Ke14LIWwMIawG3gG6RRiTVMHKlX5ldcghcMklcUcjUnGFw0PLl3vznL33hptu8ouZvn3hr39VeWyINhFMBTqZWXszqwsMAF4q8Z4XgSPMrI6ZNQB6A2plkWKGDPEy02PH+nI+kXRTv37R8NDnn8Pvf+/LTn/zGx86GjIEZs+OO8r4RJYIQgj5wGDgdfzk/lwIYZ6ZDTKzQYn3LABeA2YDHwFjQwhzo4pJKu8f/4CJE732S5cu23+/SKrr0MGHhxYvhldfhWOO8Xab3bp5eYtRo3wVUiZRiQkp04YNsP/+3n1q2jS1npSaa/VqnwcbN87vDOrV8xVI557riaJWDdh6qxITUiXXXefrtMeOVRKQmq1ZMx8emjnTL3ouuMAXR/TpA+3be7mLRYvijjI6SgRSqvfe81vkIUOgd++4oxGpHmbQs6c3Wlq+3IdFO3eG227zieYjj4Tx433OrCbR0JD8yA8/eLOZTZu8Gbn6DEimy82FCRPgiSe8xMUuu3itrXPPhcMPT499NRoakkq5/XYvLPfoo0oCIuDVUK+/3pswvfsu9O/vS09/+lOviHrrrb7zPl0pEcg25szxTWNnnumlJESkiJnfAYwdCytW+F1C27Y+h9CunZfInjAh/YrfKRHI/2zd6o0/mjSB++6LOxqR1NawoV8wvfWW90245Ra/Kzj7bC/IeO658K9/pUdZCyUC+Z+RI+HDD+GBB3wVhYhUTLt2flewcCG88w4MGAB/+xscdZRPMhe+lqo0WSyAN/3o0gWOOAImT06PyS+RVPb9995QZ/x439EcghfGO/tsn2hu0qR649FksZQrBLj4Yn985BElAZFkaNAAzjjDm+gsWeJzb2vWeBG8Pff0u4bJkyE/P+5IlQgEb/AxebKvlW7XLu5oRGqerCy49lrvm/DRRz4XN2UKnHSSv3bFFb6ZLa4BGg0NZbi1a33DTNu28P77KionUl22bPFaR+PHe02vvDxvu3nWWUU9mZNJQ0NSpquu8tvVxx5TEhCpTnXrQr9+MGmSN9MZNcr37QwbBq1be/+ECROqZxdzuYnAzLLM7Coze9HMpprZO2Y2ysxOMjMlkTT39tu+U/Lqq73yoojEo2lTn6d7/3345BMvk/3550VLUX/7W2+uE9V8QplDQ2b2BN5a8h9ADrASqI+3lzwK7yx2bQjhnWhCK52GhpJj0yY44AAfk5wzB3beOe6IRKS4ELzm11NPwXPPwfr1cOmlXgepKsobGqpTzs+NKKM3wFxgUqLZTJuqhSRxu+02X9c8ZYqSgEgqKtzFfPjh3mXtlVe8l0IUykwE22sQE0LYAqTwFgkpy5w5cM89ftt5zDFxRyMi21PYHyEq2x3nN7OfmdkMM1trZhvM7Fsz2xBdSBKl4mUkRoyIOxoRSQXlDQ0VegD4FTAnpNtaU/mR0aO9jMTTT/sElYhIRVb+LAXmKgmkv9xc7zrWp4+vUxYRgYrdEQwDJpvZv4EfCp8MIag+ZZq57DJffqYyEiJSXEUSwe3Ad/jSUXWuTVN//zu88IJPEke18kBE0lNFEsHuIYQ+kUcikfnmGxg82DeNDR0adzQikmoqMkcwxcyUCNLYDTf4FvYxY2CnneKORkRSTUUSwaXAa2a2SctH088HH3gNk8GDoVevuKMRkVS03aGhEEKj6ghEki8vz2uft2zpO4lFREpT5h2BmbUr7wfNZSU9Ikma++/3XcQPPwy77hp3NCKSqsq7I/hjosLoi8A0YBW+cqgjXnTuGOAmIDfqIKXyvvgCbr4ZTj7ZDxGRspRXa+hUM+sCnAGcB+wFfA8sACYDt4cQNldLlFIpha0na9eGhx6KOxoRSXXlzhGEEOYDN1RTLJIkzz3ntcsffNDb4ImIlEfNZWqYdetgyBDIzvba5SIi21ORDWWSRq67Dlat8l6oaj0pIhWhO4Ia5P334dFH/Y6gR4+4oxGRdFGRfgRvVeQ5iVfhnoHWrWH48LijEZF0UubQkJnVBxoAzcxsN6CwXuWuQMtqiE0q4f77Ye5cLyy3yy5xRyMi6aS8OYKLgKH4SX8aRYlgAzAy2rCkMorvGejXL+5oRCTdlLeP4EHgQTO7LISg1egpKgTvM1Crlje4FhGprIrUGnrIzA4F2hV/fwhhQoRxSQVNmgSvvAL33efzAyIilbXdRGBmTwF7AzOBrYmnA6BEELMNG+Dyy6F7d78rEBGpiorsI8gGulSlZ7GZ9QUeBGoDY0MId5XxvoOAD4D+IYTnK/s5merGG73PwN//DnW0I0REqqgi+wjmAntW9g82s9r4pPIJQBfgtETtotLedzfwemU/I5NNn+51hC6+WH0GRGTHlLd89GV8CKgRMN/MPmLb5vW/2M6f3QtYGEJYlPjzJgL9gPkl3ncZ8DfgoEpHn6G2boWLLoIWLeCOO+KORkTSXXkDCvfu4J/dClha7PtcoHfxN5hZK+CXwNGUkwjMbCAwEKBNmzY7GFb6Gz0acnLg2WehceO4oxGRdFfe8tF/7+CfbaU8V3Ke4QHgmhDCVrPS3v6/WMYAYwCys7MrPVdRk3z1FVx/PRx3HPTvH3c0IlITVGTV0Lf8+AT+DZADXFk49FOKXKD4gsYsYHmJ92QDExNJoBlwopnlhxBe2H7omemKK+CHH2DkSCgnd4qIVFhF1prch5/A/4xf5Q/AJ48/AcYBR5bxc1OBTmbWHliW+LnTi78hhNC+8GszexL4h5JA2d58EyZOhFtugU6d4o5GRGqKiqwa6htCeDSE8G0IYUNimObEEMJfgN3K+qEQQj4wGF8NtAB4LoQwz8wGmdmgpESfQTZvhksu8QRwzTVxRyMiNUlF7ggKzOw3QOH6/lOKvVbueH0IYTLe1rL4c6PLeO85FYglY919Nyxc6HcF9erFHY2I1CQVuSM4AzgTWAl8nfj6t2a2M37FLxFbuBDuvBMGDIBjj407GhGpaSpSa2gR8PMyXv5PcsORkkLwlpP16nk9IRGRZCtvQ9mwEMI9ZvYQpQwBhRAujzQyAeCvf4U33vDKonvtFXc0IlITlXdHsCDxmFMdgciPbdgAQ4dCz54+USwiEoXyNpS9nHgcD2BmDUMIG6srMIGbboIVK7zrmBrRi0hUKtKz+BAzm0/iDsHMupnZqMgjy3AzZ/pw0EUXqaiciESrIquGHgCOB9YAhBBmAT+NMKaMV1DgQ0FNm6qonIhEr0JV7EMIS0vUAtpa1ntlx40bB++/D08+CbuVuWVPRCQ5KpIIliZaVQYzqwtcTtFEsiTZ6tW+c/iII+Css+KORkQyQUWGhgYBl+JlpXOB7onvJQLXXQfffAOjRqmonIhUj4psKFuN7y6WiH3wAYwdC1deCV27xh2NiGSK8jaUlbqRrJA2lCXX1q0+QdyypS8bFRGpLuXdERTfSHYLoNNThB55BGbMgL/8BRo1ijsaEckkFsL2G36Z2YwQQo9qiGe7srOzQ05Ozdrs/PXXsO++cNBBXk5CcwMikmxmNi2EkF3aaxWZLIbtlJuWHTNsGHz/PTz8sJKAiFS/iiYCici778KECXDVVX5XICJS3cqbLC7eq7iBmW0ofAkIIYRdow6upsvP9xLTbdrADTfEHY2IZKryis5pyjJiDz8Mc+bApEnQsGHc0YhIptLQUEy++sqXifbtCyefHHc0IpLJlAhiMmyYN6T/0580QSwi8VIiiME778DTT8PVV0OnTnFHIyKZTomgmuXnw+DBPkF8/fVxRyMiUsEy1JI8I0cWTRA3aBB3NCIiuiOoVitWwI03wvHHa4JYRFKHEkE1uuYa2LRJE8QiklqUCKrJe+8V7SDeZ5+4oxERKaJEUA22bvUJ4qws7SAWkdSjyeJq8OijMHOml5jWDmIRSTW6I4jY6tXw+9/D0UfDqafGHY2IyI8pEUTshhvg22/hoYc0QSwiqUmJIEI5OfDYY3DZZdClS9zRiIiUTokgIgUFngBatFAPYhFJbZosjsiECfDBB/Dkk9C4cdzRiIiUTXcEEfjmG988dvDBcOaZcUcjIlI+3RFEYPhwWLUKJk+GWkq1IpLidJpKsgULvITEBRfAgQfGHY2IyPYpESRRCHD55bDLLnD77XFHIyJSMRoaSqIXXoApU/yOoHnzuKMREamYSO8IzKyvmX1iZgvN7NpSXj/DzGYnjv+aWbco44nSpk1wxRXQtStcfHHc0YiIVFxkdwRmVhsYCRwH5AJTzeylEML8Ym/7Avi/EMI6MzsBGAP0jiqmKI0YAYsXw9tvQx3dZ4lIGonyjqAXsDCEsCiEsAWYCPQr/oYQwn9DCOsS334AZEUYT2SWLoU77oBTToGjjoo7GhGRyokyEbQClhb7PjfxXFnOB14t7QUzG2hmOWaWs2rVqiSGmBzDhvlE8b33xh2JiEjlRZkISiuxFkp9o9lReCK4prTXQwhjQgjZIYTs5ik2C/vuuzBxoieDtm3jjkZEpPKiHM3OBVoX+z4LWF7yTWZ2ADAWOCGEsCbCeJJu61ZfLtq6te8kFhFJR1EmgqlAJzNrDywDBgCnF3+DmbUBJgFnhhA+jTCWSDz+uDecmTgRGjSIOxoRSRsh+FLD778v+9i0yY/Nm4seDz4Yjjkm6eFElghCCPlmNhh4HagNjAshzDOzQYnXRwM3Ak2BUebF+vNDCNlRxZRM69d7r4EjjoDf/CbuaEQkUnl5sGHDtse33xY9Fj++++7Hx8aNRY8bN/qJviquvjq9EgFACGEyMLnEc6OLfX0BcEGUMURl+HBYswYefFANZ0RSXl4erFvnx9q1/rh+vR/Fv16/3qtGFj4WHps3V+xzGjTw0gING0KjRv71rrtCy5ZFzxceDRoUPRY/dt5526NBA6hf37/eaadI/vNoxXsVfPyxdxy78ELo0SPuaEQySAh+Zb169Y+PNWv8WLt226/XrfMr9fLsvDM0aeI145s0gaZNoUMH/37XXf0o/nXh0ahR0bHLLlC7dnX8V0g6JYIquOIKT+S33RZ3JCI1wA8/wNdf+7FyZdFj8WPVKj9Wr/b3l6ZWLdh9dz+aNoWsLDjggKLndtut6LHwaNLEj3r1qvM3TjlKBJU0eTK8+qrvJE6xlawiqSMEvyL/6qsfHytWFB1ff+1DL6Vp2ND/kbVo4Sf1Hj38+2bN/Cj8umlTP5o0Ud33KlIiqIS8PL8b2GcfGDw47mhEYrJ5MyxbBrm52z4uX170+NVXsGXLj392l11gr71gzz39an3PPf3YY4+io0ULP7QUr9ooEVTCyJHwySfw8stQt27c0YhEID/fT+RLlhQdS5cWHbm5PjxTUqNGPiHaqpUvpWvZ0o+99tr2aNiw+n8n2S4lggpavRpuuQX69IGTToo7GpEq2rLFT+6LF/vx5Zfbfr1sGRQUbPszu+3muyZbt4bevX2YpnVrP+m3auXfN2pU/b+LJI0SQQXdeKMvPLjvPi0XlRQWgo+7f/45LFrkxxdfFD0uW+bvKVS7tp/M27WDI4+ENm28VkqbNn60bu3DOVKjKRFUwNy58Oij3mdg//3jjkYyXkGBD9MsXLjt8fnnfpTcrNSqlS+FPPpoaN/eT/rt2vnXrVqpbrooEWxPCD5BvOuuPjQkUi1C8OWSn37qxyef+ONnn/lJv/gSynr1/ETfsaPvOt17bz86dPCr+/r14/s9JC0oEWzH5Mnw5ptw//2+Qk0kqfLz/Sr+44/9WLDAHz/5xHe3Fqpb10/unTrBCSf4Y8eO/tiqlZZNyg5RIihHXh5ceaUvF73kkrijkbS2ZYtf0c+bB/Pn+wl//nx/Li+v6H0tW8K++8Jpp/njvvv6X8C2bdN216qkPiWCcjzyiJaLSiVt3eoTs3Pm+ORS4fHZZ371D3713qEDdO4MP/uZP3buDPvt52OQItVMiaAMa9fCzTfDscdquaiUYe1amD1722PuXC8ZDL68rEMHX2Hwy1/64/77+1W+xu0lhSgRlOHWW33nu5aLCiH4OvsZM7wBReGxtFgn1mbNfKfsRRfBT34CXbv6SV8bqCQNKBGU4tNP4eGH4YIL/N+0ZJCCAv8LMH160TFjRtHEba1afkV/+OHQvTt06+bHHnvoikHSlhJBKYYN8zv34cPjjkQiFYKv2Jk6FXJy/Jg+3cscgy/L7NYN+vf3gmc9eviVwc47xxu3SJIpEZTwz3/Ciy/CHXf4RZ7UIF9/DR99BB9+6I85OV6rHjzzd+8OZ58NBx7oR+fOkTUCEUklSgTFbN3qm8fatoXf/S7uaGSHbNniQzrvvw8ffOAn/8WL/bXatX0M/5RT4KCD/Nh/f530JWMpERTz1FM+B/jnP2tRR9pZuRL++9+iIyenaPdt69be9HvwYC+a1rOnShyLFKNEkPD9996MvndvGDAg7mikXFu3+mas4if+hQv9tbp1fVhn8GA45BA/WraMN16RFKdEkDBihJdhf+45Lf5IOd9840M7//1v0VDPhg3+WosWcOihMHAgHHaYJ4EMbzsoUllKBHjHvLvvhl//2s8lEqOCAt/O/f77Rcf8+b7Cx8zH9k8/3a/0Dz3U6+8oc4vsECUCvNfAli1w111xR5KB1q3zFTwffOAn/Q8/LFqzv9tuPrbfv7+f+Hv1UgkGkQhkfCKYNw8efxwuv9yLOUqE8vP9P/gHHxSt5FmwwF8z85U7p57qJ/2DD/aNW6qqKRK5jE8Ew4Z5l73f/z7uSGqg3Fw/2RceOTlFTVOaNfOZ+cJhnoMO0tW+SEwyOhG8/bb3G7jnHvUa2GEbNviJvviGreXL/bWddvJdueef7yf/Qw7x7lga2xdJCRmbCAoK4OqrvS3rZZfFHU2a+eEHmDXLSzN89JE/fvxxUS/cjh29/23v3n50766VPCIpLGMTwcSJXlbmqae0eaxc+fk+jj91atExe3ZRM5UWLXwS97TTfHinVy/Yffd4YxaRSsnIRLB5M1x/vY9WnH563NGkkK1bfelmTg5Mm+aPM2YU1dffdVfIzvY6HIWlGVq31hCPSJrLyEQwahR8+SWMHZvBi1Ly8304Z9o0vzWaNs3ra2zc6K83bOiZ8qKL/IR/4IHeHzdj/4OJ1FwZlwjWr4fbb4fjj/fuYxnhhx+8c9aMGUU19mfN8lsj8Lo7hZO5Bx7oV/377qseuSIZIuMSwd13e4fBGrt5bP36bbtozZjhO3ML++Xuuquf9AcNKiq3vM8+OumLZLCMSgTLlsEDD8AZZ/hClrRWUOBN0mfN8snbWbP8xP/ll0Xv2XNP/0VPOqmosUqHDhreEZFtZFQiuPlmP3/eemvckVTSunUwZ46f8Is/Fo7n16rlV/WHHAIXX+xdtbp390QgIrIdGZMI5s+HceO8lET79nFHU4aNGz3QefN8TL/wWLas6D277+5N0s8/39smduvmpRlUX19EqihjEsGyZX7RfMMNcUeC78JdsMCP+fOLji++KHpPvXrQpQscfbRX3DzgAD/xt2yp5ZoiklQZkwiOO87PtdV2Di0ogCVLfF1+4fHxx/5Y/Aq/bl1fodO7N5x7rl/dd+3qY/l1MuZ/j4jEKKPONElPAlu3+kn988+9Q9bChfDZZ/Dpp/5c4fJM8NU6nTv7mtX99vOvO3fWCV9EYhfpGcjM+gIPArWBsSGEu0q8bonXTwS+B84JIUyPMqZKKSjwXrhLlvhqnMWLffhm8WI/0S9e7I0MCtWt641SOnWCvn39Sn+fffxxjz00pCMiKSmyRGBmtYGRwHFALjDVzF4KIcwv9rYTgE6JozfwSOIxWlu3wpo1fpJfuRK+/hq++sqP5cu9fHLhUfxED94spX17H7M/+WS/ou/Y0Y+sLK3HF5G0E+UdQS9gYQhhEYCZTQT6AcUTQT9gQgghAB+YWRMz2yuE8FXSo3n1VRg6FFav9uWYhZUyi6tXzydjs7J8zP7Xv/bypG3b+mO7dtC4cdJDExGJU5SJoBWwtNj3ufz4ar+097QCtkkEZjYQGAjQpk2bqkWz++6+tr5ZM2je3B9btPCjeXNPAE2aaPhGRDJOlImgtDNqycvwiryHEMIYYAxAdnZ2KZfyFdC7N/zlL1X6URGRmizKWgO5QOti32cBy6vwHhERiVCUiWAq0MnM2ptZXWAA8FKJ97wEnGXuYOCbSOYHRESkTJENDYUQ8s1sMPA6vnx0XAhhnpkNSrw+GpiMLx1diC8fPTeqeEREpHSR7iMIIUzGT/bFnxtd7OsAXBplDCIiUj7VIxYRyXBKBCIiGU6JQEQkwykRiIhkOAullVpIYWa2Cvhyu28s0gxYHVE4qUq/c2bQ75wZkvU7tw0hNC/thbRLBJVlZjkhhOy446hO+p0zg37nzFAdv7OGhkREMpwSgYhIhsuERDAm7gBioN85M+h3zgyR/841fo5ARETKlwl3BCIiUg4lAhGRDFejE4GZ9TWzT8xsoZldG3c8UTOz1mb2TzNbYGbzzGxI3DFVBzOrbWYzzOwfccdSXRJtXZ83s48T/78PiTumKJnZ7xJ/p+ea2bNmVj/umJLNzMaZ2Uozm1vsud3N7E0z+yzxuFsUn11jE4GZ1QZGAicAXYDTzKxLvFFFLh+4MoTQGTgYuDQDfmeAIcCCuIOoZg8Cr4UQ9gO6UYN/fzNrBVwOZIcQuuJl7QfEG1UkngT6lnjuWuCtEEIn4K3E90lXYxMB0AtYGEJYFELYAkwE+sUcU6RCCF+FEKYnvv4WPzm0ijeqaJlZFnASMDbuWKqLme0K/BR4HCCEsCWEsD7WoKJXB9jZzOoADaiBnQxDCO8Aa0s83Q8Yn/h6PHByFJ9dkxNBK2Bpse9zqeEnxeLMrB3QA/gw5lCi9gAwDCiIOY7q1AFYBTyRGBIba2YN4w4qKiGEZcC9wBLgK7yT4RvxRlVt9ijs2ph4bBHFh9TkRGClPJcRa2XNbBfgb8DQEMKGuOOJipn9DFgZQpgWdyzVrA7QE3gkhNAD2EhEQwapIDEu3g9oD7QEGprZb+ONqmapyYkgF2hd7PssauDtZElmthOeBJ4JIUyKO56IHQb8wswW40N/R5vZ0/GGVC1ygdwQQuHd3vN4YqipjgW+CCGsCiHkAZOAQ2OOqbp8bWZ7ASQeV0bxITU5EUwFOplZezOri08uvRRzTJEyM8PHjReEEO6LO56ohRCuCyFkhRDa4f9/3w4h1PgrxRDCCmCpme2beOoYYH6MIUVtCXCwmTVI/B0/hho8OV7CS8DZia/PBl6M4kMi7VkcpxBCvpkNBl7HVxmMCyHMizmsqB0GnAnMMbOZieeuT/SOlprlMuCZxEXOIuDcmOOJTAjhQzN7HpiOr4ybQQ0sNWFmzwJHAs3MLBe4CbgLeM7MzscT4qmRfLZKTIiIZLaaPDQkIiIVoEQgIpLhlAhERDKcEoGISIZTIhARyXBKBJIRzKypmc1MHCvMbFni6+/MbFREnznUzM5KfP0vM6tyA3IzmxJV5UmRGruPQKS4EMIaoDuAmd0MfBdCuDeqz0sURzuPSuz4NbM6IYT8Ml5+CrgEuD0J4YlsQ3cEktHM7MjCPgZmdrOZjTezN8xssZn9yszuMbM5ZvZaonwHZnagmf3bzKaZ2euFJQBKOBqYXuLEfqqZfWRmn5rZEYk/6xwz+6uZvQy8YWZ7mdk7ibuVuYXvw3eYnhbdfwnJZEoEItvaGy9r3Q94GvhnCOEnwCbgpEQyeAg4JYRwIDCO0q/SDwNKFsOrE0LoBQzFd40WOgQ4O4RwNHA68HoIoTveZ2AmQAhhHVDPzJom4XcU2YaGhkS29WoIIc/M5uClSV5LPD8HaAfsC3QF3vSyN9TGSyOXtBc/rodTWARwWuLPKvRmCKGwDv1UYFwi4bwQQphZ7H0r8eqbayr9W4mUQ3cEItv6ASCEUADkhaIaLAX4hZMB80II3RPHT0IIfUr5czYBJdsp/pB43Mq2F2EbC79INCf5KbAMeKpwsjmhfuLPFUkqJQKRyvkEaF7YI9jMdjKz/Ut53wKgY2X/cDNri/dYeAyvJNsz8bwBewKLqxi3SJk0NCRSCSGELWZ2CvAnM2uM/xt6AChZ2fZVfKVPZR0JXG1mecB3QOEdwYHAB+WsKhKpMlUfFYmImf0dGBZC+CwJf9aDwEshhLd2PDKRbWloSCQ61+KTxskwV0lAoqI7AhGRDKc7AhGRDKdEICKS4ZQIREQynBKBiEiGUyIQEclw/w+RvRqbb02A+wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(t_sol, x_sol[:, 0],'b-', label='$h_1$')\n",
    "plt.plot(t_sol, x_sol[:, 1],'r-', label='$h_2$')\n",
    "plt.xlabel('Time (hrs)')\n",
    "plt.ylabel('Height (m)')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c7353d2",
   "metadata": {},
   "source": [
    "#### Using advanced integrators\n",
    "\n",
    "Most DAE systems you will encounter in Chemical Engineering have numerical challenges such as stiffness (especially common in reactive systems). It is recommended that you use more advanced scientific and numerical libraries to solve these problems. The development so far in this notebook has intentionally been working towards a common standard used by the `scipy` library among others. The `LSODA` method of `scipy.integrate.solve_ivp` can automatically detect stiffness while the solver is progressing which is particularly suited to the highly stiff coagulation systems that we deal with in this course:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "42e5b241",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "  message: 'The solver successfully reached the end of the integration interval.'\n",
       "     nfev: 117\n",
       "     njev: 0\n",
       "      nlu: 0\n",
       "      sol: None\n",
       "   status: 0\n",
       "  success: True\n",
       "        t: array([0.00000000e+00, 1.26491096e-04, 2.52982193e-04, 1.31490008e-03,\n",
       "       2.37681798e-03, 5.16432362e-03, 7.95182926e-03, 1.07393349e-02,\n",
       "       1.93900111e-02, 2.80406873e-02, 3.66913635e-02, 5.36606826e-02,\n",
       "       7.06300017e-02, 8.75993207e-02, 1.16958979e-01, 1.46318638e-01,\n",
       "       1.75678296e-01, 2.27160717e-01, 2.78643137e-01, 3.30125558e-01,\n",
       "       3.81607978e-01, 5.06043388e-01, 6.30478797e-01, 7.54914207e-01,\n",
       "       8.79349616e-01, 1.14792890e+00, 1.41650819e+00, 1.68508747e+00,\n",
       "       1.95366675e+00, 2.54202871e+00, 3.13039066e+00, 3.71875262e+00,\n",
       "       3.83642501e+00, 3.95409740e+00, 3.98706964e+00, 3.99944732e+00,\n",
       "       4.00096607e+00, 4.00248483e+00, 4.00552234e+00, 4.00855985e+00,\n",
       "       4.03893495e+00, 4.06931005e+00, 4.09968516e+00, 4.40343618e+00,\n",
       "       4.70718720e+00, 5.01093822e+00, 5.31468924e+00, 6.55454073e+00,\n",
       "       7.79439222e+00, 9.03424371e+00, 1.00000000e+01])\n",
       " t_events: None\n",
       "        y: array([[0.00000000e+00, 3.15765726e-05, 6.31340178e-05, 3.27364615e-04,\n",
       "        5.91165830e-04, 1.28259508e-03, 1.97220378e-03, 2.66038423e-03,\n",
       "        4.78909572e-03, 6.90893952e-03, 9.02153578e-03, 1.31482443e-02,\n",
       "        1.72548913e-02, 2.13442041e-02, 2.83839577e-02, 3.53836269e-02,\n",
       "        4.23476951e-02, 5.44816654e-02, 6.65292945e-02, 7.84989121e-02,\n",
       "        9.03971812e-02, 1.18890024e-01, 1.47049931e-01, 1.74913960e-01,\n",
       "        2.02509524e-01, 2.61248465e-01, 3.18990958e-01, 3.75846627e-01,\n",
       "        4.31897605e-01, 5.52160999e-01, 6.69363724e-01, 7.83855202e-01,\n",
       "        8.06404760e-01, 8.28859347e-01, 8.35143823e-01, 8.37501958e-01,\n",
       "        8.37411620e-01, 8.37321287e-01, 8.37140640e-01, 8.36960012e-01,\n",
       "        8.35154714e-01, 8.33351364e-01, 8.31549964e-01, 8.13643158e-01,\n",
       "        7.95931262e-01, 7.78414275e-01, 7.61092197e-01, 6.92408397e-01,\n",
       "        6.26971999e-01, 5.64783004e-01, 5.18592087e-01],\n",
       "       [0.00000000e+00, 4.34815864e-08, 1.04709894e-07, 1.23384805e-06,\n",
       "        2.73602362e-06, 7.56890757e-06, 1.39332805e-05, 2.14634134e-05,\n",
       "        5.03492364e-05, 8.60937474e-05, 1.27274283e-04, 2.20720062e-04,\n",
       "        3.28423592e-04, 4.48104757e-04, 6.78943240e-04, 9.36029400e-04,\n",
       "        1.21577991e-03, 1.75446280e-03, 2.34476371e-03, 2.98041489e-03,\n",
       "        3.65644112e-03, 5.43599750e-03, 7.39027055e-03, 9.49320729e-03,\n",
       "        1.17259153e-02, 1.69251930e-02, 2.25595337e-02, 2.85592669e-02,\n",
       "        3.48735304e-02, 4.96322900e-02, 6.54261578e-02, 8.20615043e-02,\n",
       "        8.54894501e-02, 8.89435510e-02, 8.99134446e-02, 9.02778257e-02,\n",
       "        9.03225195e-02, 9.03671971e-02, 9.04564878e-02, 9.05457140e-02,\n",
       "        9.14347630e-02, 9.23174135e-02, 9.31937186e-02, 1.01618254e-01,\n",
       "        1.09460338e-01, 1.16759867e-01, 1.23551263e-01, 1.46554874e-01,\n",
       "        1.63117451e-01, 1.74312269e-01, 1.79838326e-01]])\n",
       " y_events: None"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import scipy.integrate\n",
    "\n",
    "tspan = (t0, t_f) # Timespan of start and end time, h\n",
    "\n",
    "# Wrap the DAE system for scipy.integrate.solve_ivp solver\n",
    "def dxdt(t, x, y):\n",
    "    # 1. Solve algebraic equations\n",
    "    y = g_sol(y, x, t)  # New solution vector\n",
    "    # 2. Compute f(x, y, t):\n",
    "    dxdt = f_sys(x, y, t)  # solve the system dxdt = f(x, y, t) with current x, y, t\n",
    "    return dxdt\n",
    "\n",
    "# Simulate the system:\n",
    "res = scipy.integrate.solve_ivp(dxdt, tspan, x0, method='LSODA', args=(y0,))\n",
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "942ebbbb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA/XElEQVR4nO3deZzVY/vA8c/d1DStqmkXbbRoX2U3IlkeFVpsRXkK9cgu/CISlUJEhLI8tEgqQqihPEjFbJWSRKPSSkXTNDPX74/rnDln1mZqzpw5c67363Vec/Zzn2a6r+/3Xq7LiQjGGGPCV5lgN8AYY0xwWSAwxpgwZ4HAGGPCnAUCY4wJcxYIjDEmzJUNdgMKq2bNmtKoUaNgN8MYY0LKmjVrdotIrdweC7lA0KhRI1avXh3sZhhjTEhxzv2a12M2NGSMMWHOAoExxoQ5CwTGGBPmLBAYY0yYs0BgjDFhzgKBMcaEOQsEJuRVqgSNGkFsrO++s86CqCiYODFozTImZITcPgJjsqtZE379Fbp3h44d4Z9/YP16KFMGypaFYcOgaVO4775gt9SYksnOCExImzgRRo7UTl8E1qzRIABw883w6KMwZw506RLcdhpTktkZgQlpXbpAv35QvTrs2ZP1senTdXhozhyIiQlO+4wJBXZGYEJaTAz07JkzCHilpED//nDXXbB1a/G2zZhQYYHAhLS1a+G///Xdbtw46+POaTB45hmdUL7hBn2NMcbHAoEJWWlpcMUVvtudOsEvv+hwUJkyEBEB5ctDxYoaBDIyNGi0bq1DSsuX67wC6IojW2FkwpXNEZiQMnGiduIxMTB6NGzerEf9l1wC550HdevCihUwbpwGii5dYPZsaNJEr48bB8uWwerV+vx27aBPH5g6FebODfa3MyY4nHgPiUJE586dxdJQh6/YWJ0cnjBBVwWJQOXKsGiRb0I4NhZWrcp7uWh8PNx5Z9Z9B7Vq6QqjG2+EChUC/jWMKXbOuTUi0jnXxwIZCJxzPYEpQATwqoiMz/b4CcB/gZPRs5NJIjIzv/e0QGA+/xwuvliHeipWhA8/PLZVQSNHwnPP6RBSerreV7WqBokRI3R/gjGlRX6BIGBzBM65COAF4BLgNOAa59xp2Z42HFgnIu2A84HJzrnIQLXJlA5Ll2oQAO20jyUIxMbCO+/o8FK1ajBkCERHw/79emZQpw7861+waVPW19g8gimNAjlZ3BXYJCKbRSQVmA30yvYcAao45xxQGdgLpAWwTSbEJSbCeM955fDh8PLLWYd4CsI7vDR3Ljz2GLz7LixcCG+/DTNm6C7kjAw90zj1VJ2QfuEFfY1tTDOlUSADwYmA/8rtZM99/qYCLYFtQCIwUkQysr+Rc26oc261c271rl27AtVeU8KlpcFFF+n1wYN9E7z9+hUuGKxapa/znknExOjt+Hi46Sb46Sf46CNo1Uof/+ADHSqqUQP++ss3jGRMaRHIQOByuS/7hMTFQBxQH2gPTHXOVc3xIpHpItJZRDrXqpVr7WVTSk2c6Ovk/+//4I8/oHZtaNZM7/N24qtWFfw977sv53BSTIxvctm7CikpSSekvTZu1BVGjRvDSy9pTiNjSoNABoJk4CS/2w3QI39/NwHzRW0CfgFaBLBNJsR4U0jMmKFBoUwZOHIEunb1Pce/Ey9KsbGwYIHOI9SoAVddpSuKtm6FW2/VtBY33KDByf81No9gQk0gA8Eq4FTnXGPPBPAAYFG25/wGdAdwztUBmgObA9gmE2JiYmDWLPj3v3WpaFQUvPde4HMHZZ9HmDcPvvxS8xY99ZSelaSm6ga1+vV12enMmTaPYEJTwAKBiKQBI4AlwHpgroisdc7d4py7xfO0scCZzrlEYClwv4jsDlSbTGiKjfWtErrrruJJIJfXPML69XDPPZCcrKuOvDuW33hD5y1OOklvh9j2HBPmbEOZKdHWrtWUEAC33aadsX8HHWwiMGgQvPVW1vubN4cHH4QBAyDSFkSbEiAo+wiMOV7+q4QGDtQlnMeySiiQvvgCPv5Y5xGqV4fLL9eOf8MGDRD16ml21IULs77O5hJMSWKBwJRYjz0G27frGPyMGXrfsawSCpTs8wjvvQfffqvzCI8/roFh715YsgR699YNahs3+l5ncwmmpLChIVMibdgAp3n2oSck+Nb0lyT+CfC8/PMcHT6sCe+eeEIDgFe5crop7s47damqMcUhaLmGAsECQemXnq5BYONGeOAB7UhDmYhmPP33vzVNtlezZr55hPLlg9c+Ex5sjsCElCef1CDQqJHm/Ql1zun+hwMHNC1GVJSeFWzcqMtO69SBMWPANs2bYLFAYEqUTZu0U3RON3OVKxfsFh0//7mEqVM1fUXVqroKql49TVvx6KN6/brrdEezMcXJAoEpMTIydFI1PR3uuEOLxpQGue1JePddaNhQ9yN88AGccYZ+73fegTZtoFs3DRgZOTJvGVP0bI7AlBhPPaWTrA0awM8/h9/6+3XrdAL6nXc0jQboWcKDD+oQUuXKQW2eCXE2R2BKvC1btMNzDt5/P/yCAOgE+euvw44dOk9Ss6Yun/3Pf7SC2ogR+u9kTFGzQGCCTkSzeqal6bh551yPWcJHjRowapQGgXnzdIgsJUU31DVuDD16wPLllsbCFB0LBCbonn8e4uK08PykScFuTclRtqxmPI2Lg++/h759tazmZ5/BeefBKafoRruUlGC31IQ6CwQmqLZu1SRuoDtzo6KC256SqkMHnXDetg0efhhOOAE2b9YSm7Vq6dzKtuxJ3o0pIAsEJmhE4OqrdWJ0yBA488xgt6jkq11bl5ru3KmlNVu0gIMHdaK9QQMtq7lyZbBbaUKNBQITNK+9Bt99p0e0zz8f7NaElshIuPZaTYv97beax8g5XYrarZsGiP/+V2smGHM0FghMUOzYoatgQIc8KlQIbntC2emnw6JFOsx27726zHTDBq2eVqsWPPSQTjwbkxcLBCYoBgzQpGzXXw/nnx/s1pQO9evrPoRdu/Rsq2lT2L9fczU1aKCb9b791lYbmZwsEJhi9/bbWvaxenWYNi3YrSl9oqK0WtpPP+ky04sv1s5/4ULdwdyyJfTvr+mx/VmNhPBlgcAUqz17YOhQvT5nju2WDSTn4Jxz4JNPNOvpHXdAxYo6bDR3LlxyiZ6ZJSdbjYRwZ4HAFKvrr4d//tHVQt7qYybwGjaEZ57R1UbTpultEQ3GJ52kZw0PPWTDdOHKAoEpNgsX6tFp1ao6hm2KX6VKcMsteobw2Wea6ht0Ce+dd+rO5Rdf1CWpJnxYIDDFYv9+XcUCuqyxatXgtifcOae7lA8e1FxGFSro5ddftWZCzZowbFjWymqm9LJAYAJm4kRfkfnBg7Uwy+mn69p3E1z+NRKeew4WL9azhQcf1NxGhw/D9OnQvLn+zhYu1DTZpnSyQGACpksX7WwmTfKlj9i0ySYkS4LcaiTMnaupK+LitE70jTfqxrXvvtOlp7Vrw2OPWSW10sjqEZiA+vhjuOwynZisUkWPLP2LvZuS7cABePNNPbv77Te9r0wZ/Z0++KCeLTgX3DaagrF6BCZo5szxbWC64w4LAqGmShWdM9iyBb76SgMA+KqqNWkCL7+sK8FM6LJAYALm22/hjTf0+t1367JF75yBCS3OwVlnwYcfarqKMWMgOloDxC23aA2Fm2/WPQom9FggMAGRmuo7erz7bp0nmDtX5wwsGIS22rXhkUfgjz80MJxxhk4uv/aaJrvr2FEL6qSlBbulpqAsEJiAeOAB2LsXmjXTFMngm5BctSq4bTNFIyJCg/3XX2tthP/8R1ce/fCDFtGJjob774fffw92S83R2GSxKXJr10KbNtpR/PwznHxysFtkisvhw/Duu1pzed06vc85raj24IPQvbtONpviZ5PFptikp+tSQxFdamhBILyUL69pRNau1SWo112nJTe/+EJrLderp9lQ9+wJdkuNPwsEpkhNmKB7BU45RYcFTPhq00Z3ke/ZA1OmaE6jnTs1p1GdOtCrlw4rhdigRKlkQ0OmyPz6q+bAz8iAxERo1SrYLTIliQj87386bPTJJ/p3AnrWeM89sG+fZkv1X2IcG6tzSvfdF5w2lyY2NGQCTgSuvFKHhu64w4KAyck5OPtsTWfhvwT1t9/g9tt1KPGSS+Cll/T5lhq7+NgZgSkS06drkrJ69XRteWRksFtkQkF6uu4+Hz9ezxa8qlfXjKhz5sCllwavfaWJnRGYgNq5U4/oQFeMWBAwBRURAZdfrruWN2/WVNjlyukw0cGD0KePZq2Niwt2S0s3CwTmuF13na/+8FlnBbs1JlQ1bgz/+pcmvuvTR1cbpabqhHOHDron5eWXrVZCIFggMMfl/ffh88+hWjWrP2yOj39q7Pnz4dNPNXVF374aHH76SdNZREfDtdfqxjVTNCwQmGN24AAMGqTX33rL6g+b45Nbaux586BzZ9i9WzPXnn22zh3MmqWpLJo00YpqBw4Et+2hLqCTxc65nsAUIAJ4VUTG5/Kc84FngXLAbhE5L7/3tMnikuOaa2D2bOjZUyf8jCkO27bpyqIXX/RtTCtXToeVRo3SwGGpsXPKb7I4YIHAORcBbAQuApKBVcA1IrLO7znVgK+BniLym3OutojszO99LRCUDMuXa9qAChVg61Y9XTemOGVk6HDShAmwdKlvX0KDBrqEecgQHbI0KlirhroCm0Rks4ikArOBXtmecy0wX0R+AzhaEDAlw+HDOpYL8MILFgRMcJQpo7mLPv1UV66NHw/160Nysm5Qq1lTl56uWGG7l48mkIHgRGCr3+1kz33+mgHVnXNfOOfWOOcG5vZGzrmhzrnVzrnVu6xOXtDdf7+mIO7SRcsZGhNs3kynyclaB6N3bx0e+vhjOPdcTWnx6KNWZjMvgQwEuY3SZY/LZYFOwGXAxcBo51yzHC8SmS4inUWkc61atYq+pabA1q7VYudly2odYhuLNSWJc1o+8/33NQ361Kk6obxrl+5krlNHA8PHH+tmNqMCGQiSgZP8bjcAtuXynE9E5G8R2Q0sB9oFsE3mOGRk+DKLjh2rScSMKam8ZTZ//lkzoV5/vW52XLFCh4yio7VokrcWczgL5GRxWXSyuDvwOzpZfK2IrPV7TktgKno2EAl8BwwQkaS83tcmi4NnwgRdlXHqqfDjj4XPK3/kyBGSk5NJSUkJTANNwEVFRdGgQQPKlSsX7KYck5QUXaL61FOQ5NfLdOgA996r+bLKlw9e+wIpKKuGPB98Kbo0NAKYISLjnHO3AIjIS57n3AvcBGSgS0yfze89LRAER3Ky7vxMT9f/QKedVvj3+OWXX6hSpQrR0dE4G1MKOSLCnj17OHDgAI0bNw52c47bzz/D00/Dm2/6ditXqAD9+8Ndd2ka7dIkaIEgECwQFD8R6NYNvvtOl+U988yxvc/69etp0aKFBYEQJiL8+OOPtGzZMthNKTJpaVp7efx4WLnSd3+TJvr3XlrSY1vSOXNc3npLg0Ddujo8dDwsCIS20vj7K1tW576+/VbPfB94QOcPNm/WZIqPPqqbJp9+2rd3obSlx7ZAYPL155+a3wV0lZBlFjWl2YknainNnTvhs8+0vCZo8ru779Y0KpddBs8/n/UMIdRZIDA5TJyoRz2guYQOHdKNO199Fdx2GVNcypSBCy+EJUs0KEyYoKuQDh3SyzXXaK6jt9/WCehQZ4HA5NCli576PvUULFoEFStqPvjSdCpsTEFFR+vffvnycPPN+jMyUrOfXn+9prEYMEDnDEJsyjWTBQKTQ0yM5oD3Fp8vW1YLzpSWU+HPP/+cG264oUDPPXToEOeddx7p6em8/PLL1KtXj/bt22deEhMTi+yzAFJTUzn33HNJS0sr8GtMYPmnx37lFd2MVrUqPPIInHGGDhvNmQNdu2qKi8cegx07gt3qwrFAYHK1YIHv6GbkyNITBADi4+Pp0KFDgZ47Y8YMrrzySiIiIkhISODxxx8nLi4u89LmKGsM4+Pjadeu4HskIyIi6N69O3PmzCnwa0xg5ZYee+5cPVP++mvt9J94QpPd7dihAaJePQ0Mc+ZooCjpLBCYHNau9RUQv/tuLTjjnTMoDeLj49mxYwfnnHMOdevW5fPPP8/zuW+//Ta9emmuxMTERNq3b5/r87p168aWLVsA+P333+ncuXPmZ3lf88svv9CrVy86d+5M165d2bBhAwB9+/blrrvuIiYmhieffJLevXvz9ttvF82XNcftvvtyHgjFxPiWjtaurSuNtm7V4aKBAzVIrFqlQ0YnnKCFdL7/vuQOHZUNdgNMyZKR4VspMWIETJqkqyS8p8ZFdWZwxx1FX4e2fXt49tmjPy8+Pp7WrVuzYsUK5s+fz9tvv82FF16Y43mpqals3ryZRo0aAbB27Vpuuukmyni2VN92220MHToUEeG3336jYcOGACQkJGSeKXjPCI4cOcLNN9/M9OnTadq0KR999BHjx49n5syZJCYm0rJlS2I90TY9PZ1Vq1Yd97+HKX7t28Mbb8Crr+rehEmTdFnqrFl6qVcPhg3TS926wW6tjwUCk8XkyVr446STNLkc+E6FV60K/SGiI0eOsHfvXu655x4A0tLSqFatGps3b2bcuHH89ddfzJs3D4Ddu3dTzZPQfuvWrdSuXZuEhIQc77lp0yYaN26cucbeGwiOHDnC/v37qVWrFu+++y5r167lqquuyvzcc845h5SUFPbu3cvDDz+c+X4RERFERkZy4MABqlSpEsh/DhMg5cpp3eU+fbS62muvaSGd337T5HdjxkCnTnrG3acPREUFt70WCEym7dvhoYc0g+Onn2bNLBoTU7RBoCBH7oGwbt062rVrl3lUn5CQQOvWrWnSpAmvvfYaV199deZzK1SokJkXKSEhgRYtWuT6nomJiVnmClavXs2wYcNYt25d5g7c+Ph4xo0bx5AhQ7K8ds2aNZx++umULZv1v+Lhw4eJCnbvYIpEzZq68OL++zU9yzPP6NzBmjU6ZBQZCVdcoWktunULTkZfmyMwmfr313qw//kP5NHnhbzsk7cJCQm0bds21+dWr16d9PR0UlJSSExMzDMQ7N27lwoVKgCaRmPx4sW0adMmy/xAvXr1WLJkCRmeMlqJiYmICImJiTk+f8+ePdSqVStkE7uZvLVurWcHf/6pS7PPPVdTXMybB2eeqWmyH3hAL9nn5WJjdY9PIFggMADMn6/peWvWDNwfW0kQHx+fpeNNSkqidevWeT6/R48efPXVVyQmJvLWW29lLhvt0KEDBz2Zyi6++GKWLl1Kv379ePfdd4mOjqZOnTpZgs7gwYPJyMigZcuWtG/fngkTJuCcyzUQxMbGcumllwbg25uSomxZrbH85ZdaK2HyZE3quGuX5jwaPx4uukiHjg4cKIa0FiKS5wWtIXAPsBBNI70ceBEtJFMmv9cG6tKpUycxRevgQZGqVUVA5PPPA/c569atC9ybH6fdu3fLsGHDpEmTJvLEE09k3v/999/L9ddfX6xt6dOnj/z444/F+pmFUZJ/j6Fu3TqRW24RqVJF/z+CSJkyIuXLi0yYIJKWduzvDayWPPrVPOcInHMz0dKSHwITgJ1AFFpesifwkHNulIgsD1CMMsXktttg/35dHdS9e7BbExzR0dG85F0z66dDhw7ExMSQnp5OREREwNuRmppK7969ad68ecA/y5Q8LVvqcu2pUzXX0W23wS+/6LzB/ffrZPPUqUX/uflNFk+W3AvEJAHznXORwMlF3yRTnFav1nzsUVG67M3kNHjw4GL7rMjISAYOzLV0twkjERGayuLAARg9WoPDqFHg2dJS5PIMBHkEAf/HU4FNRd4iU2zS08GzmpFnn9WcKsaY4PNPa+Fdsdevn04uB8JRJ4udc5c7535wzu11zu13zh1wzu0PTHNMcRo3Tk81W7eGoUOD3RpjjFdeaS0Ctc+wIPsIngWuBBI9Ew6mFNi6VQvQlymjK4ZKYb0RY0JWbpXPinovj7+CLB/dCiRZEChd+vXT9ct33qnF6I0x4asgZwT3AR85574EDnvvFJGnA9YqE1CzZ2v+k9q1NWuiMSa8FSQQjAMOoktHrVBhiDtwAP79b70+Z46VnjTGFCwQ1BCRHgFviSkWt9wCBw/qMrTzzw92a4wxJUFB5gg+d85ZICgFVq2Cd96BChVg5sxgtyZ4rEKZMVkVJBAMBz5xzh2y5aOhy3/PwPPPQ/XqwW1PMFmFMmOyOmogEJEqIlJGRCqISFXP7arF0ThTdMaN0yWjbdtCMW6ULZGsQpkxWeWXa6iRiGzJ53EHnCgiyYFomCk6v//u2zPw/vslZM9AEEuUWYUyY7LKb7L4KedcGTTz6BpgF7py6BQgBugOPAJYICjhvHsG7r0XmjQJdmuCK68KZQsWLGDx4sXs3LmT4cOH06NHD6tQZsJGfrmG+jrnTgOuAwYD9YB/gPXAR8A4EUkpllaaYzZvHnz9tRa8GDcu2K3xE6QSZXlVKOvduze9e/dm37593HPPPfTo0cMqlJmwke8cgYisE5GHROR8EWkuIh1E5FoR+a8FgZLvn3/A2++8+67WUQ13R6tQ9vjjjzN8+HDAKpSZ8GEVykqx4cO1zkDv3nDOOcFuTcmQV4UyEeH+++/nkksuoWPHjpmPW4UyExbyqlhTUi9WoSxvEyaILFum1+PitLpRZKTIo48Gt11eJbmy1ZQpU6Rjx44ybNgwmTZtWub9VqEsp5L8ezR541gqlJnQ06WLTgzPnu1LIxEZaWcDBXH77bdz++2357jfKpSZcFCQegRLC3KfCT5vzvJevbS8XUQELFwYuNS14WLw4MHFEgTAKpSZ4MhvH0EUUBGo6ZyrDnhXn1cF6hdD28wxaNMGDh3S67fdBhdcENz2GGNKvvzOCIah+wdaeH56LwuBFwLfNHMsLr4YMjLg9NNh1iwteWeMMfnJMxCIyBQRaQzcIyJNRKSx59JORKYWYxtNAT31FHz/PVStCitW6DBRv34WDIwx+TvqZLGIPO+cOxNo5P98EXkzgO0yhXTkCIwZo9fnzdM9A/51Tm2ewBiTl6MGAufcW0BTIA5I99wtgAWCEuT++3UDWffucNFFvvsDWefUGFM6FGT5aGfgNM861EJxzvUEpgARwKsiMj6P53UBvgX6i8i8wn5OuNuyBaZM0bOA2bOD3RpjTKgpyM7iJKBuYd/YOReBTipfApwGXOPJXZTb8yYASwr7GUZdeaVOEI8dCzVrBrs1xphQk9/y0Q/QIaAqwDrn3HdkLV5/xVHeuyuwSUQ2e95vNtALWJftef8B3gO6FLr1hnfegR9+gJNP1uyixhhTWPmdEUwCJgNjgN7AE57b3svRnAhs9bud7Lkvk3PuRKAP8FJ+b+ScG+qcW+2cW71r164CfHR4+OcfGDZMr8+fr/UGzNEda6lKrxEjRmTWHsjLkCFDWLx48XG182iGDRvG//73vywlNE855RR69+5NamoqQJbH2rVrR9++ffnll18y3yN7+c1rrrnGSmWGofyWj36Z36UA751b+ZPs8wzPAveLSHouz/Vvy3QR6SwinWvVqlWAjw4Pt96qhej79YNOnYLdmqI1cWLOZa+xsXr/8TrWUpWgVca++OILUlNTOXDgQJ6vi4uLK1SJymOxcuVKunXrRkJCAk888QRxcXFs3LiRpKSkzLoJCQkJPPbYY8TFxREfH0/37t258sor8U75ZS+/OWvWLCuVGYYKkmLigKdWsf9lq3PufedcfmVOkoGT/G43ALZle05nYLZzbgtwNfCic6534b5CeEpKgjffhIoV4dVXg92aoufNm+QNBrGxertLEQwgHmupSoBHHnmE//u//+O0005j7dq1mfdv3LiRs88+mzZt2vDMM8+wY8cOGjRoAMC8efPo1q0b7dq14+yzz8Z7Vtu3b19GjBjB2WefTcOGDfnqq68YOHAgzZo1y1K34I033qBTp060bduWczyJo9avX0+zZs2IiIggMTExM7Bt2rQJEaFZs2aAprtu3bp15nvdcsst7Nixg61bt2Y+nr38ppXKDEN5ZaPzXoBH0V3GVdD0EkOBh4H+wBf5vK4ssBloDEQC8UCrfJ7/OnD10doTztlHvdlFMzJETjlFs4vee6/eHwr8s1aOHCly3nn5X9q2FSlXTuTkk/Vn27b5P3/kyIK1o23btjLB84/23nvvyY033pjr8w4fPix16tTJvJ2UlCQdOnSQjIwMGT58uLz66qsiInLkyBHp0qWLrFy5UkREbr31VrngggsyX7d79+7M62PGjJGpU6eKiEjz5s1l8uTJIiIyevRoadasmWzbtk0OHz4s1apVk5SUFNm/f7+0bNlSDh8+LCIi+/btExGRyZMny2uvvSYiIjVq1JA2bdpIixYtpHLlyvLll19mfl716tXlr7/+yvK9GjduLHFxcVle265dO+nevbuIiKSlpUnNmjXz/Pez7KOhiXyyjxZkVLmniLwsIgdEZL+ITAcuFZE5QPV8AkwaMAJdDbQemCsia51ztzjnbilswDK+o+SRI2HTJmjUCGbOLJqj5JKoenWoVw9++01/Vs/zr63g8ipVuXnzZoYMGcLVV1+d+Vz/UpUADz30EGPHjsU5R8uWLUlKSgJg/vz5tGzZkq5duwLQqlWrLEfZr7/+Ol27dqVdu3a8+OKLREVFkZKSwp9//skdd9wBQIUKFRgyZAj16tUjMjKSihUrEhkZSUREBIcOHeLuu+9m9erVme1ZsmQJPXv2zFJCc/369bzwwguMHTsW0PKaVapUoWrVqlm+//bt22nSpAlbt26lbt26JCQkEBcXl3lm5F8q04SHguwjyHDO9QO86/uv9nss370FIvIRWtbS/75cJ4ZF5MYCtCWsxcTAjBlwhWe91l9/wXvvheaGsYJUqvQOB40eDdOmwSOPHP93zatUZZMmTXjttdeyBAL/UpUrV65kyZIlxMXFMXz4cFJSUjILyiQkJNDJb5JmzZo1nH/++QC8+eabfPfddyxbtozKlStz7rnn0qpVK9auXUvHjh0z2xEfH8+tt94KQHJyMvXr18c5R8WKFUlKSuKDDz5g6NCh3Hzzzdx44438+eef1K9fn8WLF3Paab5V2e3atWPy5MmZ7WrVqlWW7z9z5kwuuOACqlSpwvLly3M87mWlMsNLQc4IrgNuAHYCf3iuX++cq4Ae8ZtiNGOG7/qIEaEZBArCGwTmzoXHHiu6vElHK1Xpz79U5YMPPsiHH37Ili1b2LJlC/Hx8ZlnBNHR0ZnX16xZw6xZszLPCBITEznzzDOpXLky7733Hl9//TVt2rQhMTExz3b4V1H76aefqFSpEgMGDODyyy8nJSWF2NhYYjy/+MTExMy6yCLCG2+8wYUXXpj5mP/8wKeffsqTTz7JpEmTMh/PLRBYqczwU5BcQ5uBf+Xx8FdF2xyTn+++gwULwDkYNUqPkktrColVq7Tz9363osqbFB8fnzmEA75SlXnp0aMHK1as4PDhw3Tv3j3z/jp16vD333+zd+9ebrjhBi699FLat29P8+bNqVatWmbnPGjQIHr16sW8efO49NJLadKkCZUqVSIxMTGzHSkpKRw6dIjqnrEv/6Awbtw4vvnmGypVqkSrVq145ZVXuOeeezLPXBITE/nyyy9ZvHgxzjm6deuWpaP/4osvWLp0KSJCy5Yt+eSTTzKL3iQmJtK7d+8c39lKZYahvCYPgPs8P58Hnst+yet1gb6E62RxerpIrVo6QTxmjN63bJlIzZq+8pQlXUmeZNy9e7cMGzZMmjRpIk888UTm/cEoVXk0HTp0kNTU1IC9/9FKZZbk36PJG8dYqnK95+fqwIUhU1BPPQW7dkGzZjpWDpZdtChFR0fz0ks5p6+Ku1RlQXz//fcBe28rlRmenBQwl5xzrpKI/B3g9hxV586dZfXq8IpNe/fqqpm0NC1BefLJwW7RsVm/fn3mkIkJXfZ7DE3OuTUi0jm3xwqyoewM59w6PGcIzrl2zrkXi7iNJh/XXgupqbpsNFSDgDGm5CrIqqFngYuBPQAiEg+cG8A2GT8rVsCSJVCjBkyYEOzWGGNKowKlKRORrdnuyjc3kCka6enQv79ef+MNrTdgjDFFrSAbyrZ6SlWKcy4SuB3fRLIJoHHjYPt2OOssuPzyYLfGGFNaFeSM4BZgOJpCOhlo77ltAmjXLi00ExEB774b7NYYY0qzowYCEdktIteJSB0RqS0i14vInuJoXLjxT718zTW6SqhfP3jrreC2yxhTuuVXoex58sklJCK3B6RFYcybVO7BB2HpUjjhBPjsM/j3v4PdMmNMaZbfHIH/Yv1HgUcC3JawFxMDs2ZBjx6++/zTLBhjTCDkV6HsDe8F2Od/23OfCYCvvwbvHr/bbw/jIBDAEmVWqtJKVZqsClrltmDbj81x2b0bHn1Ur995pyaVO95smyErgCXKrFSllao0WVm58xKkZ0/IyNCdxE8/XXSpl0ukO+6A88/P+/Loo1C/Plx8MTRsqD/r19f783qNp8jL0VipSitVabLKMxD41yoG2vrVKz7guc8UoRUrYM0aqFpVN49B1qRyYSkQJcrQQFCzZk1WrFjBiy++mGenl5qayubNm2nUqBEAa9euJSkpif79+9OyZcvMQJCWlsb111/P008/TWJiIj/99BMtWrTIfJ+YmBi+/fZb4uPjueiii5g7dy6gnXCTJk346quvGDRoEEOGDGHChAkkJSUxf/58Dh8+zIEDB5gwYQLffPMNCQkJfPDBBwB8/PHH9OzZM7NdAwcOpGXLlnTq1ImZM2dmViVLSkrKUXOgQoUK7Nu3L/O1N910E+3bt8+sY9C6dWtWhe0fXXjKc7JYRKoUZ0PCWUYGDBig1+fMgbJ+v5XSWm8gWCXK8ipVuWDBAhYvXszOnTsZPnw4PXr0OK5SlRUqVMh83euvv86cOXM4fPgwO3bs4Iknnsi3VCWQWaoyPT09s1TloEGD6NxZc4YtWbKEmTNnZilVCVoRbezYsXz22WeFKlXpz79UZZUq1g2EAxsaKgEmTIBt2+DMM3V4yBCwEmV5lars3bs3r7zySmanDbmXqhw+fDiNGjVi7NixmWcEuZWq9M4P+JeqjI+Pp3nz5nmWqjz99NOB3EtVnnXWWQwdOpQXX3yRf/75J7NUZUJCQo5SlTt37sxsV36lKnN73MtKVYYXCwRBtncvjBljO4hzyK9E2XE4WqnKxx9/nOHDdeO8laq05FbhoiC5hkwADRyoKabvvlvnQo3HffflvK8IxsnyKlUpIowaNYpLLrmEjh07Zj5upSpNWMirdFlJvZSmUpUrV2rpyRo1RAJYebDEKMklDqdMmSIdO3aUYcOGybRp0zLvt1KVOZXk36PJG8dYqtIEkIgOeQPMmGEppoPt9ttv5/bbc2ZNsVKVJhzYHEEx8t8s+9xz8OuvWoN4w4bgtsvkb/DgwSUmCARaZGQkAwcODHYzTDGzQFCMvJtlP/oIRo0C53Q3cRFsljXGmGNmQ0PFyLvwpWdPnSCOioJ580rpPgFjTMiwM4JiVru2BgGAu+4KvyAgYmmrQpn9/konCwTF7JJL9GefPjB9einNI5SHqKgo9uzZY51JiBIR9uzZYxvNSiEbGipGDzwAW7dCo0Ywf37WzbPhcGbQoEEDkpOTM5OumdATFRWVmUzPlB4WCIrJ4cOaURR0shiybpYNh0BQrlw5GjduHOxmGGOysUBQTO68U+cGBgwAz6ZToBQnlTPGhAybIygGW7fCyy/rKqFXXgl2a4wxJisLBMWgb19NNf3kk1C5crBbY4wxWVkgCLAlS2DlSjjpJBg5MtitMcaYnCwQBFBGhmYXBS0441xw22OMMbmxQBBAjz8OO3dC9+5wxhnBbo0xxuTOAkGA/PmnBoKyZWHWrGC3xhhj8maBoAj5ZxcdOBCOHIGrroKZM4PbLmOMyU9AA4FzrqdzboNzbpNzblQuj1/nnEvwXL52zrXL7X1ChTe76Kuvwgcf6AqhpUstu6gxpmQL2IYy51wE8AJwEZAMrHLOLRKRdX5P+wU4T0T2OecuAaYDpweqTYHm3Sl80UV6u0yZ8EkfYYwJXYE8I+gKbBKRzSKSCswGevk/QUS+FpF9npvfAiGfxOTXXyE9Xa+PHGlBwBhT8gUyEJwIbPW7ney5Ly9DgI9ze8A5N9Q5t9o5t7okJyw7cgRuu02vDx8O06aFV3ZRY0xoCmQgyG3VfK75h51zMWgguD+3x0Vkuoh0FpHOtWrVKsImFq1+/eDQIbjgApg6VYeF+vWzYGCMKdkCGQiSgZP8bjcAtmV/knOuLfAq0EtE9gSwPQG1YwcsXKhF6Bcu1Pv8s4saY0yB+S9B9IqN1fsDIJCBYBVwqnOusXMuEhgALPJ/gnPuZGA+cIOIbAxgWwLC/3c1YACIwODB8OKLvufExMB99wWnfcaYIhCITllEhw/27NGslBs2wA8/wP/+B59+Cn//DVdcAcOG6efceitcfjlERBzfd8lDwFYNiUiac24EsASIAGaIyFrn3C2ex18CHgaigRed5l9IE5HOgWpTUfMuF334YfjyS6hRA957T88CjDGlhPc/+owZ0L49fP651pkdNQreegt274a9e/Xy1196OXAADh7UDv3QIUhJ0Utqqk4mpqUV7LOnT896+48/ivzrAbhQKxvYuXNnWb16dbCbkWnZMl0umpEBJ5wA779vK4WMKXFEtGP2dtj79unPP/6Abdv05+7deoT+11+wf78+/9Ah7bwL2nHnJyJCx44jI6F8eb1UrAgVKkClSnqpXBmqVNHLCSdoxsqlS3VY4fHH9fXHyDm3Jq8DbStMc5w2bNAgAHD77RYEjDlmEyfq0bf/f6LYWJ1k8x9fTU/XTnz3br3s3AnJyXrZvl1v79mjnf2BA3pUnpKiwaAwIiK0sz7hBO2g//kHdu2Cxo2hY0eoVk0fq14doqP1Uq2ar0P379wrViz8sE5srJ4RjB6tSxB79gxYB2OB4DgcPqxniKB7BqZNs4pjxhSa92i9fn3o0wduukmPfL//HpYvh5NP1om3/fu1Mz58uHDvX66cHmFXrgxVq/o67tq1oU4d/dz69fW+GjV8lwoVfO/hLTDu7ZSHDw/sf/TsBc1jYgJa4NwCwXG4+mo90Lj4Ynj2WejVK7yK0RuTr/R0PYLevh1+/x1++gk2bYItW/S+Xbt0GObvv32n1aD/mfz9/LP+LFdOy/zVqOHr0GvV8nXmDRpA3bpZO/Pq1XUo5ngUc6cM6FmQ//sHuMC5BYJC8D9z3b4dFi/W7KLnnKOPh1sxehPCCjoMk5d//tHOPTlZO+off9Sfyck63r53rz7naMMxUVG+YZW6dbUz//VXWLECBg2CBx/UDr1aNf3PFgzF3CkDuf8OAjjcYIGgELyLB+bOhTFj9G+8fHk480zfc2xoyIQE/z/mmJisR71paXqks3kzJCXB+vV6fetWHX//66+jD8+ULesbgqlXT4d3mjSBFi2gUSM9gq9XL+fRevYhmEGDoFmzgP0zFEgxd8rBYIGgELwHAldeqfUGypSBRYtK1d+DCQdpadCwoa5Nv+wyraO6ebNOfF52ma6UyYtzOvHpHZJp0ACaNoXmzbWjP/FE7eSrVi18Sb5gDMEYwAJBoZ1/vu9s99prNZ2EMSXO3r06Hv/993r58Uf47Tcdl//nn6zP3ejZy/nnnzqhWreuXho2hFNPhVattLM/8USdYA3QpqagDMEYwAJBod13n54Z168Pn3yiBzH2N2qKnYh26hs3QlwcrF6tQzhbt+qSyiNHcr7GOe3oTzlFl0BWrqybo667DubNC/6RdxgMwZRUFggK4bPPYNIkvf7113o2bWeuJqD+/ls7+zVr4LvvYO1aXXWza1fenX2VKnoE37QptGkDp5+uPxs29E24eodhFi70DcHYH3PYskBQCI89pj+vu07/TzVsaGeupgCOtkJHRCdhV6/W1TJxcbrUcseOnMM4XlWq6FH9KadAu3bQtSucdppOxBZkuaQNwxg/lmLiKLz/h9u10+HRsmU1jURioiWTC3sFXYLpPfqeM0dXz7zwgq6IadTIl+ogtxQGFSvqhGzTptC2rR7Zt26tk7JRUQH/eqZ0sRQTx8G7yq5xY90fc+ONWpjeEsuZfJdg7tihG02WL4f4eN152L171tdv2KCbpGrW1NPLVq20s+/WTY/0K1YMzvcyYcfOCApgxgwYMsSXG8qGUcNAQY/2P/sM+vbV9fFr1uiQzYEDOY/wK1XSTSd79+oOxIcf1uGcqlWL5/uYsJffGUEg6xGUGpMn689Dh3TptQWBMOA92vfmoY+N1Q5/3z49JWzTRjvxHj10GdnKldr5//23Lr284AJNU7xsmWav/OAD3XgyerSu7omIsCBgSg4RCalLp06dpDgtWCACIhERIqNHi9SsKbJsWbE2wRS1CRNy/hKXLdP7vf74Q2TkSJHISJFq1fSPIPulalWRJk1EypcX6dVLJDo69z+OZcuy/uFkv21MMQBWSx79qp0R5ENEd7gDvPKKrhqyOsSlQPaj/YULtRrU4sU6Nh8VpZO0U6bo0fyff+p4/Vlnwd13a374tDRYsEAzYn78sV5/993c/zjyW6FjTAlgk8X5eOopPevv3Fkz44KtsiuxJk70ZakcMMA3eTt7tt7XtKmO7R85op13p05aUcg533j+8uU6fFOzpo75n3wyfPghjBgBL70EY8dm/aUXdAmmbZQyJZxNFmfjnSM84wzNmXXkiG663LjRlouWSBMnwkMPac6bvXt1aZfX33/rDH9Ghh7pb9umY/zZNWigm0OuvFJ/+c7lzHuT/bYxIcYmiwvBO2pw5ZW64u/yy+Hf/9b7TQmQvZB4ly7a+W/ZoqmKvRO2f/+tjx86pJky167VX+hpp8HNN8OTT+qR/+jRvqISXbv6EqXZcI4JJ3lNHpTUS6Ami/3nD+fN800QV6xoc3pBlX1id9kykRNOEBk6VGT7dpFbb9VfVG6TuSDStKnIlCkiO3ZkfQ+bvDVhBpssPjr/+cNp0/S+9HS46iobCQgq7y9m6VL45ht4/nldpz99uuaznzZNf1FlcvlTLl9eE7C1aaOTv152tG9MFjZH4Cc2Vjt+7zDy9ddrhlEbFg6CvXt1cubDDzUA7N6d9fHy5XXI59xzdahn0KCs1bDq1NFhofR0X14Q+yWaMGZzBEfhHXaOifHNNbZvrweStly0GKSnazrXO+6Ajh01PXJ0NAwbphux9u3zbb7q3VuXeVapouP7cXG+IFC9uhZWcU7LJdaooZPA/fvb0b4x+bDlo/hGH3r10pWFlSppepgbbrDlosclrzQNsbG6+/ajj+CHHzQvj3/x8ipVfK/r31/X8ffvrx3/lCn6eu8R/h9/6DBRo0bwyy++zxgwQGvqvvxycX5jY0KSBQK0Pxk1Cu65R29HRurmsSefhA4dbMn3MfNG2FGjtFrWkiW+ztorIkKXb3bpAv/6F/TpkzX1QmysBgHv+Nwff/j2BoB29AMGZD3i9z7PGFMgYR0I/A9YV6zw3d+xI9x1lwYBOxMoBBFdpjlvnnbga9fCnj2+CAt6utWmjebi6ddP0yvnV9s2+8Suf8fvP9lrvyRjjllYTxZ79wjNmqUHoykpuv8oMjLM5xYLmnnz99/1H2rJEh1L27Yt64au8uV1d2758pCUBPffD+PHF9/3MMZkssniPHjH/y+7TINAuXI6D/n++2E+QZxX5s2//tLJ12bNNB9Pgwbwn//oyp7ff9eVOpdcAs88o4XSU1L0CH7HDh3ff+21MP5HNabkCtuhIe9Bb9euvtKvffv6DnrDeoK4SxcYORIuvVTH63fu1PufeML3nOhofd555+ma23btcq7lz56WwVsb19bjGlOihG0g8B70NmigQ9utW+sQUadO+njYDDv/+acWV/nkE418v/wCBw/6Hk9J0bGy1q01+2avXlpYxeriGlNqhOUcgfdsYPlyGDNG9xtVqqT927ffhvAB69HG9nfv1pTJn36q1bR+/TVncfQqVbQmboMGOoN+223w6qsh/I9ijAGrWZzFxIna8ffr5zuoTUvT/vCuu/R2iT1gPVpH719D99RTYdIkTcFQsyY88oge3furWlWXRp1+ug4DXXCBRkTvkM6CBfpZPXrYkI4xpVjYBQJvXzl4sParXt6gUKKHhPIqlj5mjNbA/eYbfd4FF2R93fbtmpmzbVvNr33ppXr6U6FC7p9jQzrGhJWwGhryHlADXHihbzNrx4568Bz0g978jvjvvVeXZz73nO6urVZNJ3Gz//4iIrSDP3hQx/QnTdLKOmXDLuYbY/zY8lGPn3/WjasJCb4gEBEB69bp9YAnoMyeSx/0tvfUxHvE/8EHsGiRrsbp0UPX3kdG6rj9xImabO2PP7Rzb94crr4apk6FzZt1/D8qSpdrbtigidcsCBhj8pNXfuqSejnWegQTJohMnqz1xp3TVPVlyujPW28tonT0RyuKnj3v/fvva4OuukrkvPNETjxRpGzZnDn1y5UTqV9fn9O/v+bjHzUqZ6Mtz74xJg/kU48g6B17YS/HEggmTBC5/HLtP6+9Nmsf262b9pWTJ/v663zfyL9T9UYX/47eWzTFe7tmTZEZM0ReeEFk0CCR5s01EnmjkP8lIkKkdm3t9EGkT5/CFVQ5WiAyxoSt/AJBWMwRnHyybm6NjNTVQf5fefJk6LVhIgfifqb9eE8OG+9EwuzZsHWrTiikpfmGbh54QG+XLat5dCZO1PTIr74KTz+tH1C5sq7Rz02ZMjo21aABXHGFTtyed55m5PziC/2MW2/VFT/+kxYFTf1gjDHZ5DdHENCjd6AnsAHYBIzK5XEHPOd5PAHoeLT3LOwZwQEqyhHKyB6qyTpayCIuk23UkXSQIzj5lZPkwEktdPglMlJPEaKitEblCSfouJFzIqNHiyxYIHLxxXq0Hh2tR/W5DeV4L+XLi7RsKXLFFSKPPCLy2WciH3ygR/GjR9vQjjGm2BCMoSEgAvgZaAJEAvHAadmecynwsScgdANWHu19CxsI4iqfKRmQ7+WIi8i9Q4+MzLuTBw0QNWpoZw0iZ5whMmaM3vd//1f4jt6GdowxAZJfIAjkqqGuwCYR2SwiqcBsoFe25/QC3vS081ugmnOuXlE2ot2B/5FQ+cxcH3OeS1lJ16Ge7FJTdYNV69a62xZ081W1appJMzoaHnpI7x89WtMuP/OMpmEeOzZnebOj1cq9776ca1djYmzYxxgTUIEMBCcCW/1uJ3vuK+xzcM4Ndc6tds6t3rVrV6EbsnfR/0ilXGbH773kULasLr2MjNQc+TfcoGvyb7pJS5fdcAN89512+uPH61zBPffoz8ce0zz5/hMQ1tEbY0JAIANBbn1t9pnpgjwHEZkuIp1FpHOtWrUK3ZDy3c8ikiOI582zXzKlpWmnX6GCXhYt0pSk3s6+dWvdoPXkk3qUn5amt71nEy+/rGkZslfLso7eGFOCBXKnUTJwkt/tBsC2Y3jOcYmvchZnyNd5Pu5AV/GUKaNnBDNnwrhxmoNn9mzYssXX2Xs7dG/pstw6+BKdo8IYY3IKZCBYBZzqnGsM/A4MAK7N9pxFwAjn3GzgdOAvEdlelI1oejCODMqQUbkq5RrUhaZN2b1kNTXS/gAc7qQGUK8etG/ve1FaWv4dunX2xphSJGCBQETSnHMjgCXoCqIZIrLWOXeL5/GXgI/QlUObgH+Am4q6HZXlb/A0wKtmUX+IMcaEsIAmoRGRj9DO3v++l/yuCzA8kG0wxhiTv7BKOmeMMSYnCwTGGBPmLBAYY0yYs0BgjDFhLuSyjzrndgG/HuPLawK7i7A5oSDcvrN939It3L4vFN13bigiue7IDblAcDycc6slrzSspVS4fWf7vqVbuH1fKJ7vbENDxhgT5iwQGGNMmAu3QDA92A0IgnD7zvZ9S7dw+75QDN85rOYIjDHG5BRuZwTGGGOysUBgjDFhLmwCgXOup3Nug3Nuk3NuVLDbE0jOuZOcc7HOufXOubXOuZHBblNxcM5FOOd+cM59GOy2FAfnXDXn3Dzn3I+e3/UZwW5TIDnn7vT8PSc552Y556KC3aai5pyb4Zzb6ZxL8ruvhnPuM+fcT56f1Yv6c8MiEDjnIoAXgEuA04BrnHOnBbdVAZUG3C0iLYFuwPBS/n29RgLrg92IYjQF+EREWgDtKMXf3Tl3InA70FlEWqOZ5QcEt1UB8TrQM9t9o4ClInIqsNRzu0iFRSAAugKbRGSziKQCs4FeQW5TwIjIdhH53nP9ANpB5KgFXZo45xoAlwGvBrstxcE5VxU4F3gNQERSReTPoDYq8MoCFZxzZYGKFHE1w5JARJYDe7Pd3Qt4w3P9DaB3UX9uuASCE4GtfreTKeUdo5dzrhHQAVgZ5KYE2rPAfUBGkNtRXJoAu4CZnuGwV51zlYLdqEARkd+BScBvwHa0muGnwW1Vsanjrdzo+Vm7qD8gXAKBy+W+Ur9u1jlXGXgPuENE9ge7PYHinLsc2Ckia4LdlmJUFugITBORDsDfBGDIoKTwjIv3AhoD9YFKzrnrg9uq0iNcAkEycJLf7QaUwtNKf865cmgQeFtE5ge7PQF2FnCFc24LOux3gXPuv8FtUsAlA8ki4j3Tm4cGhtLqQuAXEdklIkeA+cCZQW5TcfnDOVcPwPNzZ1F/QLgEglXAqc65xs65SHSSaVGQ2xQwzjmHjh2vF5Gng92eQBORB0SkgYg0Qn+3y0SkVB8tisgOYKtzrrnnru7AuiA2KdB+A7o55yp6/r67U4onx7NZBAzyXB8ELCzqDwhozeKSQkTSnHMjgCXoaoMZIrI2yM0KpLOAG4BE51yc574HPTWkTenxH+Btz8HNZuCmILcnYERkpXNuHvA9uiruB0phugnn3CzgfKCmcy4ZeAQYD8x1zg1BA2LfIv9cSzFhjDHhLVyGhowxxuTBAoExxoQ5CwTGGBPmLBAYY0yYs0BgjDFhzgKBCQvOuWjnXJznssM597vn+kHn3IsB+sw7nHMDPde/cM4dcwFy59zngcg6aQyEyT4CY0RkD9AewDk3BjgoIpMC9XmexGiDKcRuX+dcWRFJy+Pht4DbgHFF0DxjsrAzAhPWnHPne+sXOOfGOOfecM596pzb4py70jk30TmX6Jz7xJO2A+dcJ+fcl865Nc65Jd7t/9lcAHyfrWPv65z7zjm30Tl3jue9bnTOveuc+wD41DlXzzm33HO2kuR9Hrq79JrA/UuYcGaBwJismqLprHsB/wViRaQNcAi4zBMMngeuFpFOwAxyP0o/C8ieBK+siHQF7kB3jHqdAQwSkQuAa4ElItIerTEQByAi+4DyzrnoIviOxmRhQ0PGZPWxiBxxziWi6Ug+8dyfCDQCmgOtgc805Q0RaFrk7OqRMxeON/nfGs97eX0mIt4c9KuAGZ6As0BE4vyetxPNvLmn0N/KmHzYGYExWR0GEJEM4Ij4crBkoAdODlgrIu09lzYi0iOX9zkEZC+leNjzM52sB2F/e694CpOcC/wOvOWdbPaI8ryvMUXKAoExhbMBqOWtD+ycK+eca5XL89YDpxT2zZ1zDdHaCq+gGWQ7eu53QF1gyzG225g82dCQMYUgIqnOuauB55xzJ6D/h54Fsmez/Rhd6VNY5wP3OueOAAcB7xlBJ+DbfFYVGXPMLPuoMQHinHsfuE9EfiqC95oCLBKRpcffMmOysqEhYwJnFDppXBSSLAiYQLEzAmOMCXN2RmCMMWHOAoExxoQ5CwTGGBPmLBAYY0yYs0BgjDFh7v8BKNWm1xLOgMcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(t_sol, x_sol[:, 0],'b-', label='$h_1 (Euler) $')\n",
    "plt.plot(t_sol, x_sol[:, 1],'r-', label='$h_2 (Euler) $')\n",
    "plt.plot(res.t, res.y[0, :],'b-x', label='$h_1 (Adams/BDF) $')\n",
    "plt.plot(res.t, res.y[1, :],'r-x', label='$h_2 (Adams/BDF) $')\n",
    "plt.xlabel('Time (hrs)')\n",
    "plt.ylabel('Height (m)')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a807085",
   "metadata": {},
   "source": [
    "# Discretising partial differential equations\n",
    "\n",
    "The next step is to understand how to turn the partial differential equation (PDE) from our coagulation model into a DAE system that we can simulate. This will in turn also allow us to understand the discretisation of the coagulation model term $\\frac{\\partial n\\left(v_p, t\\right)}{\\partial t}$ into a system of differential equations $\\frac{\\mathrm{d} N_i}{\\mathrm{~d} t}$ where $N_i$ corresponds to a bin interval of particle volumes $v_p$. The notebook `Coagulation 1` demonstrates this process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8dc5530",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}