Coagulation_lm1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[@lmaedler](https://www.uni-bremen.de/mvt/team/prof-dr-ing-habil-lutz-maedler) Jan. 3rd 2021"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Particle Coagulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's begin by importing some python libraries we need for this notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to declare a few parameters here that become relevant later. The most important one is \n",
    "\n",
    "        t_total\n",
    "It is the real time in seconds our coagulation is taking place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# variable declarations\n",
    "t_total = 2e-0                           # [s], change from 2e-4 to 2e-0 or any \n",
    "                                         # other number   \n",
    "nt = 1000                                # number of iteration steps \n",
    "dt = t_total / nt                        # [s], small time step during the interation\n",
    "dp_min = 0.5e-9                          # [m]\n",
    "dp_max = 1000e-9                         # [m]\n",
    "vp_min = np.pi/6 * (dp_min)**3           # [m^3]\n",
    "vp_max = np.pi/6 * (dp_max)**3           # [m^3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Population balance equation for aggregation "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The population balance equation for pure coagulation can be writte as (Smoluchowski type aggregation model):\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "& \\frac{\\partial n(v_p,t)}{\\partial t} =\n",
    "\\frac{1}{2}\\int_0^{v_p} \\beta(\\tilde{v_p},v_p-\\tilde{v_p}) \\, n(\\tilde{v_p},t) \\, n(v_p-\\tilde{v_p},t) \\, d\\tilde{v_p} \n",
    "- \\int_0^{\\infty} \\beta(v_p,\\tilde{v_p}) \\, n(v_p,t) \\, n(\\tilde{v_p},t) \\, d\\tilde{v_p} \n",
    "\\end{split}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"Hounslow et al. (1988) used a geometric grid with factor 2 based on volume $(v_{i+1} = 2 \\, v_i$, with $v$ the particle volume) and were the **first** to develop a set of equations conserving both numbers and mass for purely aggregating systems\" (Nopens et al. 2005):\n",
    "\n",
    "$$\n",
    "\\dfrac{\\mathrm{d} N_i}{\\mathrm{d} t} = \\,  N_{i-1} \\sum_{j=1}^{i-2} 2^{j-i+1} \\, \\beta_{i-1,j} \\, N_j\n",
    "+ \\frac{1}{2} \\, \\beta_{i-1,i-1} \\, N_{i-1}^2 \\\\\n",
    " - N_{i} \\sum_{j=1}^{i-1} 2^{j-i} \\, \\beta_{i,j} \\, N_j\n",
    "- N_{i} \\sum_{j=i}^{\\infty} \\beta_{i,j} \\, N_j\n",
    "\\label{eq1}\\tag{1}\n",
    "$$\n",
    "\n",
    "Please consult the paper (Nopens et al. 2005) for more information of the strenghts and weaknesses of this model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialization\n",
    "vp_num = int(np.log(vp_max / vp_min) // np.log(2) + 2)\n",
    "\n",
    "vp = np.zeros((vp_num))\n",
    "for i in range(vp_num): \n",
    "    vp[i] = vp_min * 2**i \n",
    "    \n",
    "n = np.zeros((vp_num))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coagulation in the free molecular regime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we want to consider coagulation in the free molecular regime with \n",
    "\n",
    "$$\n",
    "\\beta_{i,j} = \\left( \\frac{3}{4 \\pi} \\right)^{1/6} \\sqrt{\\frac{6 k_B T}{\\rho_p}}\n",
    "\\left( \\frac{1}{v_i} + \\frac{1}{v_j} \\right)^{1/2}\n",
    "\\left( v_i^{1/3} + v_j^{1/3} \\right)^{2}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 293                                  # [K]\n",
    "rho_p = 1.0e3                            # [kg / m^3]\n",
    "k_B = 1.380649e-23                       # [J / K = N m / K = kg m^2 / s^2 / K]\n",
    "beta_pre = (3/(np.pi*4))**(1/6) * (6 * k_B * T / rho_p)**(1/2)\n",
    "# [kg m^2 / s^2 / K * K * m^3 / kg = m^5 /s^2]**(1/2) = [m^5/2 / s]\n",
    "\n",
    "beta_n = np.zeros((vp_num, vp_num))\n",
    "for i in range(vp_num): \n",
    "    for j in range(vp_num): \n",
    "        beta_n[i,j] = beta_pre * (1/vp[i] + 1/vp[j])**(1/2) * \\\n",
    "                      (vp[i]**(1/3) + vp[j]**(1/3))**2\n",
    "        # [m^5/2 / s * 1/m^3/2 * m2 = m^3/s]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have all necessary equation, but still need to define our starting conditions. This is done with *n_0*. We select a bin (let's say number 4) and fill this bin with all our particles (let's say 2e18 #/m^3). You can use any other bin or several bins to start the coagulation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial conditions\n",
    "n_0 = np.zeros((vp_num))\n",
    "start_bin = 4                 # this is an easy way to select a starting size\n",
    "n_0[start_bin] = 1e17         # [#/m^3]\n",
    "#n_0[6] = 1e17                # try to use these additional particles as starting condition\n",
    "n = n_0.copy()                # copy the initial condition to our interation array "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, in order to solve equation ([1](#mjx-eqn-eq1)) we need to implement a numreric scheme for the time derivertive. Here, we use the forward in time decritization. From our declarations above we use \n",
    "\n",
    "            nt = 1000                                # number of iteration steps \n",
    "            dt = t_total / nt                        # [s] \n",
    "      \n",
    "With this in place we can compute the following equation. Note, that the summations are done individually. This take a lot of time, but I have not come up with a vecor based way yet..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{split}\n",
    "N_i^{n+1} = N_i^{n} + \\Delta t \\cdot \\left(\n",
    "N_{i-1}^n \\sum_{j=1}^{i-2} 2^{j-i+1} \\, \\beta_{i-1,j} \\, N_j^n\n",
    "+ \\frac{1}{2} \\, \\beta_{i-1,i-1} \\, {(N_{i-1}^n)}^2 \\\\\n",
    "- N_{i}^n \\sum_{j=1}^{i-1} 2^{j-i} \\, \\beta_{i,j} \\, N_j^n\n",
    "- N_{i}^n \\sum_{j=i}^{\\infty} \\beta_{i,j} \\, N_j^n\n",
    "\\right)\n",
    "\\end{split}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn = np.empty_like(n)\n",
    "for t in range(nt - 1):\n",
    "    nn = n.copy()\n",
    "    for i in range(1, vp_num, 1):\n",
    "        sum1 = np.zeros((vp_num))\n",
    "        sum2 = np.zeros((vp_num))\n",
    "        sum3 = np.zeros((vp_num))\n",
    "        for j in range(0, i - 1, 1):\n",
    "            sum1[i] = sum1[i] + (2**(j-i+1) * beta_n[i-1,j] * nn[j])\n",
    "        for j in range(0, i - 0, 1):    \n",
    "            sum2[i] = sum2[i] + (2**(j-i) * beta_n[i,j] * nn[j])\n",
    "        for j in range(i, vp_num - 1, 1): \n",
    "            sum3[i]= sum3[i] + (beta_n[i,j] * nn[j])\n",
    "        n[i] = nn[i] + dt * ( \\\n",
    "                             nn[i-1] * sum1[i] + \\\n",
    "                             1/2 * beta_n[i-1,i-1] * nn[i-1]**2 - \\\n",
    "                             nn[i] * sum2[i] - \\\n",
    "                             nn[i] * sum3[i]     \n",
    "                             )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, all the calculations are done. Let's look at the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conservation of mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "During coagulation the total volume of all particles must be consered. This can be expressed by the thrid moment of particle ditribution.\n",
    "$$\n",
    "V_{total} = \\int_0^\\infty v_p \\cdot n(v_p) \\, d v_p\n",
    "$$\n",
    "\n",
    "Let's calculate the ratio of $V_{total}$ at the beginning and the at the end. This must be 100%!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ratio of total particle volumes: 100.0 %\n"
     ]
    }
   ],
   "source": [
    "print(\"Ratio of total particle volumes:\", round(np.sum(vp * n) / np.sum(n_0 * vp) ,3) * 100,\"%\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "During coagulation the total number of all particles should decrease. \n",
    "$$\n",
    "N_{total} = \\int_0^\\infty n(v_p) \\, d v_p\n",
    "$$\n",
    "\n",
    "Let's calculate the ratio of $N_{total}$ at the beginning and the at the end. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ratio of total particle numbers: 0.619 %\n"
     ]
    }
   ],
   "source": [
    "print(\"Ratio of total particle numbers:\", round(np.sum(n)/np.sum(n_0),5) * 100,\"%\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can we check our simulation results further? We know that for a simple **monodisperse coagulation** with a constant $\\beta$ (let's say $10^{-15} \\, m^3/s$):\n",
    "\n",
    "$$\n",
    "\\dfrac{\\mathrm{d} N(t)}{\\mathrm{d} t} = - \\frac{1}{2} \\cdot \\beta \\cdot N(t)^2; \\quad \\quad t=0, \\,N=N_0\n",
    "$$\n",
    "integration leads to\n",
    "$$\n",
    "N(t) = \\frac{N_0}{1 + N_0 \\cdot \\frac{1}{2} \\cdot \\beta \\cdot t} \\\\\n",
    "$$\n",
    "using the fact that the total volume is conserved during coagulation\n",
    "$$\n",
    "v_p \\cdot N_p = v_{p,0} \\cdot N_{p,0} \\\\\n",
    "$$\n",
    "it follows \n",
    "$$\n",
    "\\frac{d_p}{d_{p,0}} = \\left( \\frac{v_p}{v_{p,0}} \\right)^{1/3}  = \\left( {1 + N_0 \\cdot \\frac{1}{2} \\cdot \\beta \\cdot t}  \\right)^{1/3} \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated factor of diameter increase: 4.66\n"
     ]
    }
   ],
   "source": [
    "print(\"Estimated factor of diameter increase:\", \\\n",
    "      round((1 + np.sum(n_0) * 1/2 * 1e-15 * t_total)**(1/3), 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's compare to our numeric solution. Here, for simplicity we just take the particle volume bin with the largest number concentration. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computed factor of diameter increase: 5.04\n"
     ]
    }
   ],
   "source": [
    "maximum_index = np.argmax(n)\n",
    "maximum_index_start = np.argmax(n_0)\n",
    "print(\"Computed factor of diameter increase:\", \\\n",
    "      round((vp[maximum_index]/vp[maximum_index_start])**(1/3), 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We don't expect these two values to be the same, but they should be in the same order of magnitude.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we would like to plot the orginal particle size distribution with the one after the coagulation. For this we need to initialize the particle diameters from our particle volumes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "dp_start = np.zeros((vp_num))\n",
    "dp_end = np.zeros((vp_num))\n",
    "dp_start = (6/np.pi * vp[start_bin])**(1/3)\n",
    "dp_end = (6/np.pi * vp)**(1/3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plotting looks a bit more complicated than usal, mainly because we want a plot with two axes. Please also note that particle size distributions are normilzed to the total particle number, respectively.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABBSElEQVR4nO3deZgU1dX48e8RxIERFJdEBaV5E7UKZABlcxeNiDbuCm5BYpSgcXt/0djGqBhj0om8Ro1RRMU1UYkKKq2JcSG4gIIGFahCFltBMSIoiIgycH9/3Jqxp6eX6p6te+Z8nqef6aq6det0dU+frqpb94oxBqWUUqrUbNXSASillFKZaIJSSilVkjRBKaWUKkmaoJRSSpUkTVBKKaVKkiYopZRSJUkTlCqYiBgR+WHwfKKIXN3I9Y8RkVcas84Ct3+wiCwKUa7R4xSRM0XkucasM8M2DhORFSnTC0TksEaqu078qZ+VRqp/vYj8T2PVp0pb+5YOQJU3Y8y4lo6hoUTEAHsaY5YAGGNeBvZuiViMMX8F/trM2+ydr4yIRID3ga2NMdU56mq0+EVkBvCQMebulPq3bYy6VXnQI6hWTkT0R0gWum8al+5P1dg0QZUoEUmKyGUi8o6IrBWRR0WkImX5eSKyRETWiMhTIrJbyjIjIj8XkcXA4ppTOiLySxH5VERWisgJInKMiLwX1PGrlPUHicgsEfkiKHubiHTIEud9IvLb4PnTwSmYmscWERkTLHNE5F/BthaJyMiUOnYMXsM6EXkD+EGO/RIJXt9YEfk4iO8XYWPPsG9mBoveDmIeleEU2O4i8oSIrBKR1SJyW5bYsr7GDGXHiMgyEflSRN4XkTNT5r8SPP9l2v7cJCL3Bcu2E5F7gtf4kYj8VkTaZdlWx+B9+lxEFgID05YnReRHKftvbvBe/FdEbgqK1eynL4JY9g9ifVVE/iQia4Dxkvm05zHBa/1MRG4Uka2CbY0XkYdS4qh5b9uLyA3AwcBtwfZuS3n/ak4vbyciDwTvywci8uuUuseIyCsiMiF43e+LyNHZ3g9Voowx+ijBB5AE3gB2A3YAPGBcsOxw4DNgX2Ab4M/AzJR1DfCvYL2OwGFANXANsDVwHrAK+BvQGegNbAT+J1h/P2AI9hRwJNj2pWn1/zB4fh/w2wzxDwc+BnYHKoHlwE+COvcN4u8dlH0EmBKU2wf4CHgly36JBNt/OCjfJ3gtPyog9tp9k/56gunDgBXB83bA28Cfgu1VAAcFy8bUxJnvNaa9hkpgHbB3ML1ryr6orTNtnd2D/XlMMD0NuDOo63vYz8rPsuyzOPBy8Jp3B+bXvL6Uz1rN/psF/Dh4vi0wJG2/t09Zbwz2c3VR8Jo7pscfrPNSsO09gPeAc4Nl47Gn8NLf2/bB9Iyaslk+ew8AT2I/w5Gg7p+mxLYJ+1lvB5wf7D9p6f9tfYR/tHgA+sjyxtgvjbNSpv8ITAye3wP8MWXZtsE/YySYNsDhKcsPA74G2gXTnYMyg1PKvAmckCWWS4GpKdM5ExSwF/ApcHAwPQp4Oa3MncC1wZfHJsBJWfY78ieo1PJ/BO4pIPbD08rkSlD7YxNg+wx1134Z53qNGdarBL4ATiZIkpnqTJnXMXh/rgimvw98k7oucDrwUpZ9sAwYnjI9luwJaiZwHbBTlv2enqA+zBV/sE7qti8AXgiej6fIBBV8br4BeqUs+xkwIyWOJSnLOgXr7tIU/6/6aJqHnuIrbZ+kPN+ATURgj6o+qFlgjFkPrAa6pZRfnlbXamPM5uD518Hf/6Ys/7qmfhHZS0Smi8gnIrIOmzB2ChOwiGyH/VV7tbGNDQB6AIOD025fiMgXwJnALsDO2F/fqfF+QH7p5XcrIPb0fZPL7sAHJkfDgECu11iHMeYrbEIbB6wUkYSIODnqvgdYZIz5Q8q2tg7WrdnWndgjqUx2I/z+/Sn2B4YvInNEZESOshBuX2Z8rxpoJ6ADdV/LB9T9H6j9/zHGbAieaiOLMqIJqjx9jP2SAkBEKoEdsafGajSkm/o7AB/bsq0L8CtA8q0UnP//G/aX/J0pi5YD/zbGbJ/y2NYYcz726KQamwhq7BEixvTyHxcQeyH7Zjmwh+RvAJDrNdZjjPmnMeZI7Ok9H7grUzkRiWFbFP40bVvfYI9yarbVxWRvjbeSkPvXGLPYGHM6Ntn9AXgs+Hxl22dh9mW29+or7JFNjfRknqvuz7BH3j1S5u1B3f8BVeY0QZWnvwE/EZF+IrIN9ijhdWNMspHq74y9RrI++GWf8Us2gxuwp68uSZs/HdhLRH4sIlsHj4Ei4gZHdU9gL7B3EpFewNkhtnV1UL439rrPow2I/b9Atntr3sB+wcdFpFJEKkTkwAzlsr7G9IIi8n0ROS744v8GWA9szlDuaOBi7KnXmqNejDErgeeA/xORLiKylYj8QEQOzfIapgBXikhXEemOvWaUkYicJSI7G2O2YE9DEsS2CthC9v2Uy+XBtnfHfjZq3qt5wCEiskdw5H1l2npZ35fgczMFuEFEOotID+D/AQ9lKq/KkyaoMmSMeQG4Gngc++X5A+C0RtzEZcAZwJfYX/aP5i5e63RsA4XPU1qenWmM+RIYFsT4MfbUyx+wDTwALsSeevkEe03r3hDb+jewBHgBmGCMqbk5tJjYxwP3B6fL6rS8C74Ij8Ve9/gQWIE9PUdauXyvMdVWwC+CcmuAQ7HXZtKNwp4C9VL258Rg2WjsKa6FwOfAY9ijsUyuw57+eh+b2B7MUg5s45YFIrIeuAU4zRizMThFdgPwarCfhuSoI92T2Gto84AE9pQlxph/Yd+fd4Ll09PWuwU4JWiFd2uGei/CHoUtA17B/nCbXEBcqsSJMTpgoSofEvKGUaVU+dMjKKWUUiVJE5RSSqmSpKf4lFJKlSQ9glJKKVWSNEEFJKVPuZDlG32YCVV+svQ915zbb7GhQZqCiMwQkXNbOg5VGsomQUlKh5ZNUT5PXfX+uY0x44wx1zdG/U1BRKJBZ5lfBL0q3CUinZt4myIifxDboepqEfmjiEjK8oiIvCQiG0TEb6z3J09Mhf7wKKh8c5O08ZWMMS8bY1pkaJC2pNDProjsKyIzg1sD/isil6Qsu15E3hWRahEZn7aeiMhVIvKh2A57HxGRLjm2c3PQDH+WiHRLmX+miNzSgJdcEsomQamCbQf8FtutjAt0B25saKVie/qekWXxWOAEoC9QBYzA9o9W42HgP9heL67C9lKwc0NjagtEh7IoSPBF35jfb6E/uyKyE/APbPdTO2LvoUsdhHIJ8EvsPWHpRgM/Bg7E/u92xHYGnWk7g7CdI++CvQ/symD+dtj7Aa8p5AWWpJbuDDDMA3tj4RZsf3HrgV8G848DFmDveJ8BuHnK/x17A+VabKeYvVO2cR+Ze+V2sT19bw7q+iK9PEHnotgP3afYm2dPAI7B9rC8BvhVSp1bATFgKbYPvSnADk28D08C3s2x3MH28r0GWASMzFLuMIIOOTMsew0YmzL9U2B28HwvbK8JnVOWv0zQQ3uGugYBc7G9QvwXuKmI1zwW2x3Ot8F793TKezoj+NwsAI7LU77mvfoSe2PsiSnbGEP+jm3HYm/KXQn8Iu01zgriWAncBnRIWW6AnwOLsfd+zQzmfRXEN4qUjm2DdXbH9syxKvhs3ZYpzrDvd1B2BnA98GqwD54j6Ew2ffvBvCTfdT47Hvt/91Cw7rvBZ+FK7P/KcmBY2rZ+j+3BYy32Jt8dUpYPwX7OvsD2Mn9Y2ro3BHF+TUoHwA383yn0s/s74MEQ9T4EjE+b9xhwecr0Adjvn04Z1h8F/D54Phx4Jnh+G3BGY7z2ln6UxRGUMebH2Lv4jzW2f7M/ishe2F81l2Lvtn8GeFpEOmQqH1T1LLAntp+xtwgx8qcxxsN26jkrqGv7LEV3wQ7F0A37y+Uu4CzsL5yDgWvku6GqL8YmsEOxv5I+B/6SqdKgG5gvcjzOyPcaAodgv4wzbaMS+2X1N+y+OR24XWw3QoXojf3SqPF2MK9m2TJje1zItDzdLcAtxvan9wNsEi+IMWYS9j3+Y/DeHSsiWwNPY79kv4ftjeCvIrJ3pvJBVUux7+F22F4ZHhKRbL02ZDIU+7kbBsRSTg9tBv4X2/Hp/sAR1O9R4gRgMLbX7kOCeX2D+Or0kiF2PKjp2F4jItjP4iPpwRT5fp+B7VLqe9geLC7L96JTHIv90dgVexTyT+yPtG7Ab7BHGqlGA+dg/zeqgVuDuLthjzp+ix2+4zLg8bQjmR9jfxB0JkOnuGI7Es72v5Tek0WNQj+7Q4A1IvKa2PHXnhaRMP1Lgu03UtKmt8F+ftItAA4WkY7Yz84CERmAHcblbyG3V9LKIkFlMQpIGGP+ZYzZBEzAHg4fkG0FY8xkY8yXxphvsL/s+gaHw41hE3BDEMsj2C+dW4LtLcB+mKqCsj8DrjLGrEiJ5ZRMp3GMMR+auh2Qpj/yfhBF5Ehs/3bZDvlHAEljzL3GmGpjzFvYbpROKWgP2O6K1qZMrwW2Da5DpS+rWZ7tutgm4IcispMxZr0xZnaBsWQzJIglboz51hjzIvZL/fRsKxhj/m6M+dgYsyVICouxRz9hXWeM+coY8y62G6fTg3rfNMbMDvZ5EvtFnd6f3u+NMWtMSl98OQzCfqlfHmxvozEmU8OIYt7ve40x7wVxTAH6hYinxsvGdo5bjT2a2hm7/2v+VyIisn1K+QeNMfON7fX9amBkkHzPwh4lPBO8F//CHmUfk7LufcaYBcHr2pQeiDFmRI7/pWw9txf62e2O/X+7BNuB7fvYH9NhPAucG1zz2g64IpjfKb2gMWY+9n2bHWznD9gfdheLyMXBNbC/pu3bslLOCSp9yIkt2NMF3TIVFpF2IhIXkaVih2FIBotCDSMRQujhLLA9ME+V74ZK8LC/pr/fSLHUEttn2t+AU4wx72UplnOoCBGJpcyfDhyUVrbGeiD1gm4XYL0xxmRYVrP8SzIrdNiHsHYDlgeflxrpwzTUISKjRWReyuvdh8I+N2U7NEiKbEO/hJH+f/BZhv+V1PrS99fW2P3SAzg1Le6DqNsHYSH7K6xCP7tfY8cgm2OM2Yg96j4g5I/hydhkNgP7o/alYP6KTIWNMX8yxvQ1xowiGJMM+70+FntU5WFPUZelckpQ6XcUpw85Idh/0I+ylD8DOB74EfZUTaRm1SK23VDLgaPTfr1VGGPqDRUQnOJbn+NxZraNiEh/4CngHGM7mM0VT9ahIowx8Zr52F/fr6SWTalnAbaBRI2+fHdacQHwP1K3JWHq8jpM9mEfCpXpc7N72gX01GEa6pQX20v2XdgObXcMXu98wn1uapT10CB51BkyIzjSaWjDl/T9tQk7vMZy7NFVatyVxph4Svmc+0tEns3xv/RsltUK+uxiO79NjaPmed7PTHBkeK0xJmKM6R5s4yPyDCMiIt/Hnpn5DfYH1DvBEeQcvjtzU3bKKUGld70/BYiKyBHBdYVfYC9kvpalfOdg+WrsP9TvCtx2dxHpUGTs6SZihwnoASAiO4vI8ZkKBqf4ts3xyHgdTUT2wbYkusgY83SeeEIPFZHHA8D/E5FuIrIb9j25L3gd72F7s75W7JAVJ2L/cR7PEn+2YR9qbiEYEzKm9M/B69gv1V8Gr/Mw7DWSR7KUrxkLaVWw7Z9gvwAKUbZDg4TwHlAh9raGrYFfk7kH90KcJSK9RKQT9gv3seCI6yHgWBE5KjgjUiG2VWn3sBUbY47O8b90dJZ1CvrsYk/jnih2OJytsacpXzHGfAEQ7O8K7Pdv+6DOdsGyHcQOnSJih565CfhN2hF/JjdhR2/egD2lOFBEtsU2YlkWcveUnHJKUL8Hfh0c2l9mjFmEPSf9Z+yvq2OxjSK+zVQe++X5AfaXyELseduwXsT+kvlERD5rhNdyC/bI5jkR+TKIZXAj1JvqF9hfsvek/EKs/cUn9kbjiVDwUBG53IltgPAu9igjQd0L4KcBA7CNQuLY046rstSVcdiH4EfCjoR//+4BegWfg2nB5+M44Gjs5+Z2YLQxxs9SfiHwf9jWdv8F+mBbiRWinIcGyckYsxbbsONu7P/WV2Q5HVWAB7E/bD7BNjy6ONjWcuxZkF9hfzAsBy6neb7Hsn52xd4svb6moLHXNX+F/fx/in0/Uhsz3YU9DXg6tsn619jGHWBPZT6D3Y/PApONbbyTlYgMBbY3xkwNtv9GsO3l2AY68RyrlzTti0+VFRE5CPh5cPqvpIkODaJUg2iCUqqJaIJSqmHK6RSfUkqpNkSPoJRSSpUkPYJSSilVksquA8qtttrKdOzYsaXDUEqpkrVhwwZjjCn7A5CyS1AdO3bkq6++aukwlFKqZIlImK6xSl7ZZ1illFKtkyYopZRSJUkTlFJKqZJUdteglFLf2bRpEytWrGDjxo0tHYpqARUVFXTv3p2tt966pUNpEk2WoCKxxGRsz9efJuPRep1rRmIJwfaxdgy2+/4xyXj0raaKR6nWaMWKFXTu3JlIJILt0F+1FcYYVq9ezYoVK+jZs2dLh9MkmvIU333YDj+zORo7SuSe2LFL7mjCWJRqlTZu3MiOO+6oyakNEhF23HHHVn303GQJKhmPzgTW5ChyPPBAMh41yXh0NrB9JJYoZBhtpRRocmrDWvt735LXoLpRd/TLFcG8lS0TTtP62+sf8uS8nGOOAXB8v26cMXiPZoioeX3+6BTWTZ+et1yXESPoOmpk3nJKqdavJVvxZUr9GTsGFJGxIjJXROZWV5dnp9BPzvuIhSvX5SyzcOW6UEmsHK2bPp2Nvp+zzEbfD5XEVOlIJpPss09h4zf6vk+/fv3o378/S5cubaLImtf48eOZMGFCzjLTpk1j4cKFtdPXXHMNzz//fFOHVtZa8ghqBXWHdu7Od0Nh1xEM2DUJoLKysmx7t+21axce/dn+WZePunNWM0bT/Cochx4PPpB1+Qc/Ht2M0aiWMm3aNI4//niuu+66UOWNMRhj2Gqr8r4rZtq0aYwYMYJevXoB8Jvf/KaFIyp9LfmOPwWMjsQSEoklhgBrk/Foqzy9p1RrVl1dzdlnn01VVRWnnHIKGzZsAODNN9/k0EMPZb/99uOoo45i5cqVPPPMM9x8883cfffdDB06FICbbrqJffbZh3322Yebb74ZsEdmrutywQUXsO+++7J8+XJuvPFGBg4cSFVVFddee23GWP7xj3+w77770rdvX4444ggA1qxZwwknnEBVVRVDhgzhnXfeAeCNN97ggAMOoH///hxwwAEsWrQIgA0bNjBy5EiqqqoYNWoUgwcPZu7cuQBsu+22tdt67LHHGDNmTL0Y7rrrLgYOHEjfvn05+eST2bBhA6+99hpPPfUUl19+Of369WPp0qWMGTOGxx57DIAXXniB/v3706dPH8455xy++eYbACKRCNdeey377rsvffr0wc9zFqK1acpm5g8DhwE7RWKJFcC1wNYAyXh0InZY42OwQ2FvAH7SVLEo1RZc9/QCFn6c+zRyoXrt1oVrj+2ds8yiRYu45557OPDAAznnnHO4/fbbueSSS7jooot48skn2XnnnXn00Ue56qqrmDx5MuPGjWPbbbflsssu48033+Tee+/l9ddfxxjD4MGDOfTQQ+natSuLFi3i3nvv5fbbb+e5555j8eLFvPHGGxhjOO6445g5cyaHHHJIbRyrVq3ivPPOY+bMmfTs2ZM1a2wbrWuvvZb+/fszbdo0XnzxRUaPHs28efNwHIeZM2fSvn17nn/+eX71q1/x+OOPc/vtt9O1a1feeecd5s+fT79+/QraZyeddBLnnXceAL/+9a+55557uOiiizjuuOMYMWIEp5xySp3yGzduZMyYMbzwwgvstddejB49mjvuuINLL70UgJ122om33nqL22+/nQkTJnD33XcXFE85a7IElYxHcw7JnYxHDfDzptq+Uqp57L777hx44IEAnHXWWdx6660MHz6c+fPnc+SRRwKwefNmdt21fiPdV155hRNPPJHKykrAfrm//PLLHHfccfTo0YMhQ4YA8Nxzz/Hcc8/Rv39/ANavX8/ixYvrJKjZs2dzyCGH1N4TtMMOO9Ru4/HHHwfg8MMPZ/Xq1axdu5Z169Zx9tlns3jxYkSETZs21Za/5JJLANhnn32oqqoqaH/Mnz+fX//613zxxResX7+eo446Kmf5RYsW0bNnT/baay8Azj77bP7yl7/UJqiTTjoJgP32248nnniioFjKnfYkoVQrke9Ip6mkN3UWEYwx9O7dm1mzcl9XzTVgak3Sqil35ZVX8rOf/SxnXZmaXWfahohw9dVXM3ToUKZOnUoymeSwww7LG1Nq/dnuPxozZgzTpk2jb9++3HfffcyYMSNrffm2B7DNNtsA0K5dO8q1kVixyvuqo1KqxX344Ye1iejhhx/moIMOYu+992bVqlW18zdt2sSCBQvqrXvIIYcwbdo0NmzYwFdffcXUqVM5+OCD65U76qijmDx5MuvXrwfgo48+4tNPP61TZv/99+ff//4377//PkDtKb5DDjmEv/71rwDMmDGDnXbaiS5durB27Vq6desGwH333Vdbz0EHHcSUKVMAWLhwIe+++27tsu9///t4nseWLVuYOnVqxv3x5Zdfsuuuu7Jp06ba7QJ07tyZL7/8sl55x3FIJpMsWbIEgAcffJBDDz00Y91tjSYopVSDuK7L/fffT1VVFWvWrOH888+nQ4cOPPbYY1xxxRX07duXfv368dprr9Vbd99992XMmDEMGjSIwYMHc+6559aexks1bNgwzjjjDPbff3/69OnDKaecUu/Lfuedd2bSpEmcdNJJ9O3bl1GjRgG2CfjcuXOpqqoiFotx//33A/DLX/6SK6+8kgMPPJDNmzfX1nPBBRewatUqqqqq+MMf/kBVVRXbbbcdAPF4nBEjRnD44YdnPGUJcP311zN48GCOPPJIHMepnX/aaadx44031mteX1FRwb333supp55Knz592GqrrRg3blzY3d+qSb7Dy1JTWVlpynHAwpom5GGamecqU65qmpCHaWaeq4yqy/M8XNdt6TBalc2bN7Np0yYqKipYunQpRxxxBO+99x4dOnRo6dAyyvQZEJENxpjKLKuUDb0GpRosTC8RG32fipRfk0qVqg0bNjB06FA2bdqEMYY77rijZJNTa6cJSjVYTS8RuRJQhePQZcSIZoxKqeJ07ty59r4n1bI0QalGka+XCKWUKpQ2klBKKVWS9AhKKaXaIM9xh2MHjW0H3O36Xjxt+ZnAFcHkeuB81/feDrNuY9EjKKWUamM8x20H/AU7cGwv4HTPcXulFXsfONT1vSrgeoIOu0Ou2yg0QSmlmsTf//53XNdl6NChzJgxI+N9UOXisMMOy9tw4uabb67tKBfgmGOO4YsvvmjiyIo2CFji+t4y1/e+BR7BDiJby/W911zf+zyYnI0dcSLUuo1FE5RSqkncc8893H777bz00ktFJahy69YnPUE988wzbL/99i0VTvuaMfSCx9i05dkGjM3mp8CzRa5bNE1QSqkGOeGEE9hvv/3o3bs3kyZNAuxYR6+88grjxo3j1FNPZeLEifzpT3+iX79+vPzyy6xatYqTTz6ZgQMHMnDgQF599VXA9vowduxYhg0bxujR9ccH++Mf/0ifPn3o27cvsVgMgHnz5jFkyBCqqqo48cQT+fxz+6M/07AXAEuXLmXIkCEMHDiQa665pnYIjRkzZjAi5VaICy+8sE4XSDXOP/98BgwYQO/evWuH/bj11lv5+OOPGTp0aO0wIpFIhM8++wzIPaTIeeedR+/evRk2bBhff/11g96LFNXGmAEpj0lpy0MPGOs57lBsgqq5HhV63YbSRhJKtRKf/O53fOM17nhB27gOu/zqVznLTJ48mR122IGvv/6agQMHcvLJJ3PNNdfw4osvMmHCBAYMGMD48eNrh9gAOOOMM/jf//1fDjroID788EOOOuooPM8D7DhSr7zyCh07dqyznWeffZZp06bx+uuv06lTp9q+9kaPHs2f//xnDj30UK655hquu+46br755qzDXlxyySVccsklnH766UycOLHgfXLDDTewww47sHnzZo444gjeeecdLr74Ym666SZeeukldtpppzrlcw0psnjxYh5++GHuuusuRo4cyeOPP85ZZ51VcExFCDVgrOe4VcDdwNGu760uZN3GoAlKKdUgt956a23HqcuXL2fx4sXsuOOOOdd5/vnn6wx/vm7dutq+9Y477rh6yalmnZ/85Cd06tQJsMNprF27li+++KK2c9Wzzz6bU089Fcg+7MWsWbOYNm0aYBNlTdIMa8qUKUyaNInq6mpWrlzJwoULcw7JkWtIkZ49e9aON7XffvuRTCYLiqUB5gB7eo7bE/gIOA04I7WA57h7AE8AP3Z9771C1m0smqCUaiXyHek0hRkzZvD8888za9YsOnXqxGGHHZZ1GIpUW7ZsYdasWRkTUeowG6myDaeRTaHDXrRv354tW7bUTmd6He+//z4TJkxgzpw5dO3alTFjxuR9vbn6O60ZSgPscBqNeIovJ9f3qj3HvRD4J7ap+GTX9xZ4jjsuWD4RuAbYEbjdc1yAatf3BmRbtyni1GtQSqmirV27lq5du9KpUyd832f27NkZy6UPNTFs2DBuu+222ul58+bl3dawYcOYPHly7bWkNWvWsN1229G1a1defvlloO5QFdmGvRgyZEjtAIaPPPJI7fwePXqwcOFCvvnmG9auXcsLL7xQL4Z169ZRWVnJdtttx3//+1+effbZ2mXZhtMIO6RIc3N97xnX9/Zyfe8Hru/dEMybGCQnXN871/W9rq7v9QseA3Kt2xQ0QSmlijZ8+HCqq6upqqri6quvrh0BN92xxx7L1KlTaxtJ3HrrrbVDYPTq1SvUtaDhw4dz3HHHMWDAAPr168eECRMAuP/++7n88supqqpi3rx5XHPNNUD2YS9uvvlmbrrpJgYNGsTKlStrh9LYfffdGTlyJFVVVZx55pkZh/3o27cv/fv3p3fv3pxzzjm1IwkDjB07lqOPPrq2kUSNsEOKqPp0uI1m0pqH22isYTJ0uI3C6XAbhduwYQMdO3ZERHjkkUd4+OGHefLJJ1s6rKLpcBtKKdVKvPnmm1x44YUYY9h+++2ZPHlyS4ekstAEpZRqUw4++GDefvvtlg5DhaDXoJQqc+V2ml41ntb+3uc8gorEEoMAk4xH50RiiV7AcMBPxqPPNEt0SqmcKioqWL16NTvuuGNBTbBV+TPGsHr1aioqKlo6lKw8xx0EGNf35gQdyg4HfNf3QuWQrAkqEktci+2ttn0klvgXMBiYAcQisUT/ZDzaZE0LlVLhdO/enRUrVrBq1aqWDkW1gIqKCrp3756/YAvwHLc2h3iOWyeHeI7bP0zz9FxHUKcA/YBtgE+A7sl4dF0klrgReB3QBKVUC9t6663p2bNnS4ehVCb1cojre+s8xw2dQ3Jdg6pOxqObk/HoBmBpMh5dB5CMR78GtuRYTymllKp2fW+z63sbgKWu760DcH0vdA7JlaC+jcQSnYLn+9XMjMQS24WtXCmlVJv1ree49XKI57ihc0iuU3yHJOPRbwCS8WhqZVsDZxcYqFJKqbblENf3vgFwfa+oHJI1QdUkpwzzPwM+C1N5JJaoM259Mh6Npy3fDngI2COIZUIyHr03TN1KKaVKV01yyjA/dA7Jex9UJJYYX1hYtevVG7c+aKqe6ufAwmQ82hc4DPi/SCzRoZjtKaWUKj2e444vdt1czcy3Au4CPi2y7kHAkmQ8uiyor2bc+oUpZQzQORJLCLAtsAYor3GelVJK1eM5bkNzSM4jqKeBNcl49Moi6w4zbv1tgIsdjfFd4JK0610AiMhYEZkrInOrqzV/KaVUGXgaWOP6XrE5JGeCGgBMLbZiwo1bfxQwD9gN217+tkgs0aXeSsZMMsYMMMYMaN9euw9USqky0NAckjNBDQXujMQSg4usO8y49T8BnkjGoyYZjy4B3gcclFJKlbuhwJ2e4xabQ7InqGQ8uhB7hHNjkXXPAfaMxBI9g4YPpwFPpZX5EDgCIBJLfB/YG1hW5PaUUkqVCNf3GppDcrfiS8ajHwPRYipOxqPVQM249R4wJRmPLojEEuMiscS4oNj1wAGRWOJd4AXgiqAZu1JKqTLn+l7ROQRCjAeVjEe/LLbyoNfzZ9LmTUx5/jEwrNj6lVJKlTbX94rOIXkTVCSWGABcBfQIygt2CI6qYjeqVDYbfb926PdsuowYQddRI5spIqVUQ3iOmzGHuL6XN4eEaRL3V+BybDNw7YNPNZkuI0bkLbPR9wE0QSlVPorOIWES1KpkPJreuEGpRtd11Mi8iSff0ZVSquSscn2vqBwSJkFdG4kl7sY2YqjtWykZjz5RzAaVUkq1Kdd6jlsvh7i+lzeHhElQP8Hem7Q13x2eGUATlFJKqXyKziFhElTfZDzap/jYlFJKtWF9Xd8rKofk7c0cmJ2hF3KllFIqjNme4xaVQ8IcQR0EnB2JJd7Hnj/UZuZKKaXCOgg423PcOjmksZqZD29gcEoppdquonNIrvGg5gKvAs8CM5Lx6MZiN6KUUqpt8Ry3Tg5xfa/gHJLrCGoI9tBsOHBdJJZYje1X79lkPPpeEfEqpZRqO+rkEM9xa3OI63uhckjWBBV09jojeBCJJXbFDt/+20gs8UNgdjIevaAh0SullGqdXN+rk0M8x63NIZ7j/hCY7fpezhwSevS/ZDy6EpgMTA6Gg9+/uLBVOfn80Smsmz49Z5mNvk+Fo8N4KVVOPMcdDtwCtAPudn0vnrbcAe4F9gWucn1vQsqyJPAlsBmodn1vQL7tub5Xm0OC4eDz5pC8zcwjscS0lOcHAiTj0S3JePTVfOuq8rdu+vTa/u+yqXCcUP3oKaVKg+e47YC/YI9oegGnZ2gKvga4GJhAZkNd3+uXLzl5jjst5fmBAK7vbXF9L28OydVIIgG8AvSKxBKdkvHoBuDP2Gyq2pAKx6HHgw+0dBhKqcYzCFji+t4yAM9xHwGOBxbWFHB971PgU89xixrPyXPc2hziOW4n1/cKziG5jqDOBBYA2wFPRWKJWUCPYMBBt5iAlVJKNYv2IjI35TE2bXk3YHnK9IpgXlgGeM5z3Dc9x02vu0adHOI57iygh+e44zzHDZVDciWoi4H1wCfJePRHyXh0f+AT4GvgF2FfhVJKqWZXbYwZkPKYlLZcMqxjCqj/QNf39sWeIvy557iHZChTm0Nc3/uR63sF55BcjSRmA0OBSCSWeAVYCnQG3gH0fI9SSpWvFcDuKdPdgY/DrhwM5Y7re596jjsVe8pwZlqx2hziOW5ROSTrEVQyHn0uGY9eDSwDDgZ+h22x8VPseUWllFLlaQ6wp+e4PT3H7QCcBoQas8lz3ErPcTvXPAeGAfPTy7m+95zrew3KIWGamf8+GY8aYFEklvgsGY9eGKZipZRSpcn1vWrPcS/E3jjbDpjs+t4Cz3HHBcsneo67CzAX6AJs8Rz3UmyLv52AqZ7jgs0hf3N97x85Nvd71/cMsMhz3M9c3wudQ/ImqGQ8OiXl+cCwFSullCpdru89AzyTNm9iyvNPsKf+0q0D+hawnSkpzwvKIWGG21BKKaWanSYopZRSJUkTlFJKqZIUqi++SCzRDeiRWj4Zj6Y3KVRKKaXq8Ry3Xg5xfS9vDsmboCKxxB+AUdguMDYHsw3127wrpZRSdXiOW3QOCXMEdQKwdzIe/abYAJVSSrVZJwB7u75XcA4Jk6CWAVtjx5IvSCSWqNOdezIejWcocxhwc7CNz5Lx6KGFbkcppVTJKjqHhElQG4B5kVjihdQNJOPRi3OtFIklarpzPxLbrcacSCzxVDIeXZhSZnvgdmB4Mh79MBJLfK/QF6CUUqqkbQDmeY5bJ4e4vpczh0C4BPUUIbvASDMIWJKMR5cBRGKJet25A2cATyTj0Q8BkvHop0VsRymlVOkqNoeE6kni/mIqJnN37oPTyuwFbB2JJWZgOxG8JRmP1utEMOgqfixAhw4digxHKaVUc3N9r9gcknPAwinJeHRkJJZ4lwzdsCfj0ao8dYfpzr09sB9wBNARmBWJJWYn49H36qxku4qfBFBZWVlIl/BKKaVagOe4U1zfG+k5bsYc4vpevhyS8wjqkuBvsWN5h+nOfQW2YcRXwFeRWGImto+n91BKKVXOGppDsieoZDy6Mvj7QZF1zwH2jMQSPYGPsN25n5FW5kngtkgs0R7ogD0F+Kcit6eUUqpEuL63MvhbbA5puq6OkvFoNVDTnbsHTEnGowuCIePHBWU84B/YAazewDZFrzeuiFJKqbYnVFdHxUrGo/W6c0/GoxPTpm8EbmzKOJRSSpWfnEdQkViiXSSWeKi5glFKKdV6eI7bznPconNIzgSVjEc3AztHYglt262UUqogru9tBnYOhpUvWJhTfEng1Ugs8RTwVe3MePSmYjaolFKqTUkCr3qOWyeHuL6XN4eESVAfB4+tsDfTKqWUUmEVnUPC9CRxHUAklqgM7ldSSimlQnF97zoAz3ErXd8rKIeEGQ9qf+AeYFtgj0gs0Rf4WTIevaCYYJVSSrUdnuPWySGe4/YFfub6Xt4cEuY+qJuBo4DVAMl49G3gkKKjVUop1ZbcTEoOcX0vdA4JdaNuMh5dnjZrc8aCSimlVBrX94rKIWEaSSyPxBIHACZobn4xtmcIpZRSKp/lnuMeAJiguXnoHBLmCGoc8HPs8BkfAf2CaaWUUiqfonNImFZ8nwFnNiA4pZRSbZTre0XnkDCt+P4HuAUYgh3TYxbwvzUj5SqllFLZeI6bMYe4vpc3h4Q5xfc3YAqwK7Ab8Hfg4aKjVUop1ZYUnUPCNJKQZDz6YMr0Q5FY4sKCQ1RKKdUWiet7dXKI57ihckiuId93CJ6+FIklYsAj2MOzUUCi2EiVUkq1fp7j1uYQz3GLyiG5jqDeDCqTYPpnKcsMcH1B0SqllGpLGpxDcg353rNBoSmllGqzXN9rcA4J04qvHRAFIqnldbgNpZRS+XiOmzGHNNZwG08DG4F3gS3FhaiUUqqUeI47HNv8ux1wt+t78bTlDnAvsC9wlet7E8Kum6boHBImQXVPxqNVhVSqlFKqdAVHNX8BjgRWAHM8x33K9b2FKcXWYLslOqGIdVN1d32vqBwS5j6oZyOxxLBiKldKKVWSBgFLXN9b5vret9gWdsenFnB971PX9+YAmwpdN82znuMWlUPCHEHNBqZGYomtgkAFMMl4tEsxG1RKKdXk2ovI3JTpScaYSSnT3YDUHsZXAIND1l3ourOBqZ7j1skhru/lzSFhEtT/AfsD7ybjUROivFJKqZZVbYwZkGO5ZJgX9vu90HVrc4jrewXlkDCn+BYD8zU5KaVUq7EC2D1lujvwcROtuxiYX2hygnBHUCuBGZFY4lngm5qZ2sxcKaXK1hxgT89xe2KHwDgNOKOJ1l0JzPAct04Oaaxm5u8Hjw7BQymlVBlzfa866A/vn9im4pNd31vgOe64YPlEz3F3AeYCXYAtnuNeCvRyfW9dpnVzbK7oHBJmPKjrCqlQKaVU6XN97xngmbR5E1Oef4I9fRdq3RzbKTqHhOlJ4iUyXABLxqOHF7tRpZRSbYPnuBlziOt7eXNImFN8l6U8rwBOBqrDBBaJJercbZyMRzPebRyJJQZimyKOSsajj4WpWzXc549OYd306TnLbPR9KhynmSJSSrVCReeQMKf43kyb9Woklvh3vvWCPvzq3G0ciSWeSsajCzOU+wP2fKZqRuumT8+bgCochy4jRjRjVEqp1sT1vXo5xHPcvDkEwp3i2yFlcitgP2CXEHUPApbUDA0fiSVq7jZO7w7jIuBxYGCYgFXjqnAcejz4QEuHoZRqpVLGhYLCckioU3ypY3pUY1tj/DTEennvNo7EEt2AE4HDyZGgRGQsMBagQwdtSKiUUmWk2BwS6hRfsWN6hLnb+GbgimQ8ujkSyz7AYtBFxySAyspKvWFYKaXKREPGhQpzBEUkljiA+uNB5TsvFOZu4wHAI0Fy2gk4JhJLVCfj0Wlh4lJKKVX6PMetl0Nc38t7bSHMNagHgR8A84DNwWwD5Kt8DrBnJJbIerdx6tFZJJa4D5iuyUkppVoPz3GLzSGhjqAGAL0K7YsvGY9WR2KJOncbJ+PRBZFYYlywfGLOCpRSSrUGA7A9UDRJX3zzsS0uVhZaeTIerXe3cbbElIxHxxRav2qbNvo+H/x4dM4yXUaMoOuokc0UkVIqh6JzSJgEtROwMBJLvEHdzmKPK3RjSjVUmHuyNvo+gCYopUrDTsBCz3Hr5BDX9/LmkDAJanzxcSnVuLqOGpk38eQ7ulJKNavxxa4Yppl5qDt+lVJKqXSu7xWdQ8IMWKiUUko1O01QSimlSlKoBBWJJTpGYom9mzoYpZRSrY/nuB09xy04h+RNUJFY4ljsDVb/CKb7RWKJpwqOUCmlVJvjOW6dHOI5bj/PcUPlkDBHUOOxPZN/AZCMR+dhu6xQSiml8hlPSg5xfW8eIXNImARVnYxH1xYXl1JKqTau2vW9onJIqJ4kIrHEGUC7SCyxJ3Ax8FoxG1NKKdXmzPcc9wygnee4BeWQMEdQFwG9sXcAPwysAy4tLk6llFJtTNE5JMyNuhuAq4KHUkopFZrre0XnkKwJKhJLPE39AQZraV98SimlsvEcN2cOaWhffBOKCUoppZSiEXJI1gRV0wdfJJaoBL5OxqNbgul2wDYN3bBSSqnWq6YPPs9xK4GvXd/bEkyHziFhGkm8AHRKme4IPF9YqEoppdqoonNImARVkYxH19dMBM875SivlFJK1ahwfa82hwTPQ+WQMAnqq0gssW/NRCSW2A/4uuAQlVJKtUVfeY5bm0M8xw2dQ8LcqHsp8PdILPFxML0rMKrQCJVSSrVJlwJ/9xy34BwS5j6oOZFYwgH2BgTwk/HopiIDVUop1Ya4vjfHc9w6OcT1vVA5JOspvkgscXjw9yTgWGAvYE/g2GCeUkoplZHnuIcHf+vlkGBeXrmOoA4FXgwqTmeAJwqKVimlVFvS4ByS6z6oa4Onv0nGo++nLovEEj0LCFIppVSJ8Rx3OHAL0A642/W9eNpyCZYfA2wAxri+91awLAl8CWzG9lY+IL1+1/dqc4jre3VyiOe4oXJImFZ8j2eY91iYypVSSpWe4GbZvwBHA72A0z3H7ZVW7GjsKbk9gbHAHWnLh7q+1y9TckpTdA7J1Refg+2Bdru0a05dgIowlSullCpJg4Alru8tA/Ac9xHgeGBhSpnjgQdc3zPAbM9xt/ccd1fX91aG2UDQMKI3sF3aNafQOSTXNai9gRHA9tQ9h/glcF6YypVSSrWI9iIyN2V6kjFmUsp0N2B5yvQKYHBaHZnKdANWYq8hPec5rgHudH1vEvU1OIfkugb1ZCSWmA5ckYxHfxemMqWUUiWh2hiT69SbZJiX3vN4rjIHur73see43wP+5Tmu7/rezNSCru896TnudOAK1/eKyiE574NKxqObI7HEkUBRlUdiiToX4ZLxaDxt+ZnAFcHkeuD8ZDz6djHbUkopFdoKYPeU6e7Ax2HLuL5X8/dTz3GnYk8ZzkxbH9f3NnuOW3QOCdOTxGuRWOI24FHgq5qZyXj0rVwrBb2e/wU4EvtC50RiiaeS8WjqOc73gUOT8ejnkVjiaGAS9Q8zlVJKNa45wJ5Ba7qPgNOAM9LKPAVcGFyfGgysdX1vZdA7+Vau730ZPB8G/CbHtl7zHLdeDqlpEZhLmAR1QPA3NQADHJ5nvUHAkmQ8ugwgEkvUuwiXjEdTx6Wfjc3QSimlmpDre9We414I/BN7hmuy63sLPMcdFyyfCDyDbWK+BNvM/CfB6t8HpnqOCzaH/M31vX/k2FyxOSRUV0dD85XJIsxFuFQ/BZ7NtEBExmKbOdKhQ4ciw1FKKVXD9b1nsEkodd7ElOcG+HmG9ZYBfQvYTrE5JNQRFJFYIoptLljbNDAZj+Y6pINwF+Fq6h+KTVAHZVoetD6ZBFBZWZl1CGH1nc8fncK66dNzltno+1Q4TjNFpJRqqzzHrZdDXN/Ll0Py36gbiSUmYnuevQibdE4FeoSIKcxFOCKxRBVwN3B8Mh5dHaJeFcK66dPZ6Ps5y1Q4Dl1GjGimiJRSbZHnuMXmkHDXoJLxaFUklngnGY9eF4kl/o9w/fDNAfYMukXKeBEuEkvsEdT142Q8+l6YgFV4FY5DjwcfaOkwlFJt2wGu71V5jvuO63vXeY4bNoeE6uqoZmCpDZFYYjdgE5C3H6VkPFoN1FyE84ApyXh0QSSWGBeJJcYFxa4BdgRuj8QS8yKxxNws1SmllCpPtTnEc9zQOQTCHUFNj8QS2wM3Am9hryPdFabyZDxa7yJcMh6dmPL8XODcMHUppZQqS9M9x92eInJImFZ81wdPHw96lqhIxqNriwxUKaVUG+L6Xm0OCXqWqHB9L1QOyZugIrFEBXABtoWdAV6JxBJ3JOPRjcUGrJRSqm3wHLdeDvEc9w7X9/LmkDCn+B7Adu7352D6dOBBbEsMpZRSKpeic0iYBLV3Mh5NvSnrpUgsof3lKaWUCmNv1/fq5BDPcUPlkDCt+P4TiSWG1ExEYonBwKsFBqiUUqpt+o/nuLU5xHPc0DkkzBHUYGB0JJb4MJjeA/AiscS7gEnGo1WFRquUUqrNGAyM9hy3Tg7xHPddwLi+lzWHhElQwxshQKWUUm1T0TkkTDPzD4qtXCmlVNvm+l7ROSTMNSillFKq2WmCUkopVZI0QSmllCpJmqCUUkqVpFADFipVbjb6Ph/8eHTOMl1GjKDrqJHNFJFSqlCaoFSrE2YQxprBHDVBKVW6NEGpVqfrqJF5E0++oyulVMvTa1BKKaVKkiYopZRSJUkTlFJKqZKk16DKzOePTmHd9Ol5y230fSocpxkiUkqppqFHUGVm3fTptS3QcqlwnFCt2ZRSqlTpEVQZqnAcejz4QEuHoZRSTUqPoJRSSpUkTVBKKaVKkiYopZRSJUkTlFJKqZKkCUoppVRJ0gSllFKqJGkz8xKy3zsz6LPodT54pUvWMnoDbuPRITmUKm1NmqAiscRw4BagHXB3Mh6Npy2XYPkxwAZgTDIefaspYyplfRa9zi6rPoRd98laRm/AbRw6JIdq6zzHrfP97PpePG15ve9n1/feCrNuY2myBBWJJdoBfwGOBFYAcyKxxFPJeHRhSrGjgT2Dx2DgjuBvm/XJznvQX2/CbXI6JIdqyzzHrff97DnuU67v5f1+Drluo2jKI6hBwJJkPLoMIBJLPAIcD6S+iOOBB5LxqAFmR2KJ7SOxxK7JeHRlYwfz5LmX0SG5pLGrDe3Ub6qp3KZ9ztN3u6z6kKVddmPUnbOaMTKVzZiV69hl1Ycs/NEJLR2KaoO+jfyQ4++e0FTVDwKWuL63DMBz3Kzfz67vGWC257jbe467KxAJsW6jaMoE1Q1YnjK9gvpHR5nKdAPqJCgRGQuMBejQoUOjB9ocKrdpz47bbpOzzKaee7KsW7/mCUjl9e7ebfpgXpW39iIyN2V6kjFmUsp0Q76fw6zbKJoyQUmGeaaIMgQ7dhJAZWVlveVhNOEvkUbTA/uzRpWK/Vs6AKWKVW2MGZBjeUO+n0N9bzeGpkxQK4DdU6a7Ax8XUUYppVTjasj3c4cQ6zaKpkxQc4A9I7FET+Aj4DTgjLQyTwEXBtenBgNrm+L6k1JKqTrmAHt6jpv3+zm4xjQYWOv63krPcVeFWLdRNNmNusl4tBq4EPgn4AFTkvHogkgsMS4SS4wLij0DLAOWAHcBFzRVPEoppSzX9+p9P7u+t8Bz3HGe4+b8fs62blPEKcY0yanDJlNZWWm++uqrlg5DKaVKlohsMMZUtnQcDaVdHSmllCpJmqCUUkqVJE1QSimlSlLZXYMSkS3A1xkWtQM251k9V5lsy8LOT59uD1TniacxhXn9jbV+c+7rMPOac183534OU77Y5bqvCy9bTt8fHY0x5X8AYoxpFQ/sndJFl8m2LOz8DNNzS+31N9b6zbmvw8xrzn3dnPs5TPlil+u+Lrxsa/7+KNVH+WfY7zzdwDLZloWdH2b7Tamh2y9k/ebc12HnNZfm3M9hyhe7XPd14WVb8/dHSSq7U3zlQkTmmtxdjahGovu6+ei+bh66n63WdARVaiblL6Iaie7r5qP7unnofkaPoJRSSpUoPYJSSilVkjRBKaWUKkmaoJRSSpUkTVBKKaVKUlOOB6WyEJFewHhgNfCCMeaxlo2o9RKRPYDbgM+A94wx8RYOqVUSkYOBM7HfKb2MMQe0cEitlohsBVwPdMHe0Ht/C4fUZPQIqkAiMllEPhWR+Wnzh4vIIhFZIiKxPNUcDfzZGHM+MLrJgi1zjbSv9wISxphzgF5NFmwZa4z9bIx52RgzDpgOtNovzIZqpM/08UA3YBN21NtWS5uZF0hEDgHWAw8YY/YJ5rUD3gOOxH5g5gCnY/vX+n1aFecEf68FNgAHGGMObIbQy04j7evNwGOAAR40xtzbPNGXj8bYz8aYT4P1pgDnGmPWNVP4ZaWRPtPnAJ8bY+4UkceMMac0V/zNTU/xFcgYM1NEImmzBwFLjDHLAETkEeB4Y8zvgRFZqvp58MF8osmCLXONsa9F5DLg2qCuxwBNUGka6zMdnE5dq8kpu0b6TK8Avg0mG9KZbsnTU3yNoxuwPGV6RTAvIxGJiMgk4AHgxiaOrbUpaF8D/wAuFpGJQLIJ42ptCt3PAD9FfwAUo9B9/QRwlIj8GZjZlIG1ND2CahySYV7Wc6fGmCQwtsmiad0K3dfzgVZ7CqQJFbSfAYwx1zZRLK1doZ/pDdgfA62eHkE1jhXA7inT3YGPWyiW1k73dfPQ/dx8dF9noQmqccwB9hSRniLSATgNeKqFY2qtdF83D93PzUf3dRaaoAokIg8Ds4C9RWSFiPzUGFMNXAj8E/CAKcaYBS0ZZ2ug+7p56H5uPrqvC6PNzJVSSpUkPYJSSilVkjRBKaWUKkmaoJRSSpUkTVBKKaVKkiYopZRSJUkTlFJKqZKkCUqVDRE5IRhLq2b6NyLyoxzlDxOR6Q3Y3n0ickrw/O7UbTc2EdleRC5ohHqaNE6lmpMmKFUWRKQ9cAIpYzoZY64xxjzfHNs3xpxrjFnYhJvYHigoQYlV53+4GeJUqtloglLNIujB3ReR+0XkHRF5TEQ6BcuuEZE5IjJfRCaJiATzZ4jI70Tk38AVwHHAjSIyT0R+kHaEM1BEXhORt0XkDRHpnLb9ymCwuDki8h8ROT5DjCIit4nIQhFJAN9LWTZDRAYEz+8QkbkiskBErkspkwzinRUs31dE/ikiS0VkXEq5y4M43klZPw78IHhtN2YrF+xHT0RuB96ibh9u6XGuF5Ebgn0yW0S+n+E1jw/ek+eC+E8SkT+KyLsi8g8R2TrcO6xU49MEpZrT3sAkY0wVsI7vjhhuM8YMDAZw60jdMXC2N8Ycaoy5Ads/2eXGmH7GmKU1BYL+yx4FLjHG9AV+BHydtu2rgBeNMQOBodhEV5lW5sQgxj7AeUC2YcuvMsYMAKqAQ0WkKmXZcmPM/sDLwH3YntSHAL8JYh0G7IkdA6gfsJ/YQexiwNLgtV2eo1zNfnzAGNPfGPNBlhgBKoHZwT6ZGbymTH4ARLEjtT4EvGSM6YPdh9Ec9SvVpDRBqea03BjzavD8IeCg4PlQEXldRN4FDgd6p6zzaIh69wZWGmPmABhj1gX9m6UaBsREZB4wA6gA9kgrcwjwsDFmszHmY+DFLNsbKSJvAf8JYk295lPTyee7wOvGmC+NMauAjSKyfRDHsGDdtwAHm4jS5Sr3gTFmdpbYUn2LHYId4E0gkqXcs8aYTUHM7bBjaNW8hmzrKNXkdDwo1ZzSO340IlIB3A4MMMYsF5Hx2ORR46sQ9UqGujOVOdkYs6jAGOtWItITuAwYaIz5XETuo2683wR/t6Q8r5luH8Txe2PMnWn1RjLEm61cmH0CsMl819nmZrL/v38DYIzZIiKp62zJsY5STU6PoFRz2kNE9g+enw68wndf7p+JyLbkHlzwS6Bzhvk+sJuIDAQQkc5Bo4pU/wQuSrm+1T9DPTOB00SknYjsij0VmK4LNkGsDa7pHJ0j3kz+CZwTvFZEpJuIfC/Da8tWTqk2Q38dqebkAWeLyJ3AYuAOY8wGEbkLezopiR0bJ5tHgLtE5GJSEpkx5lsRGQX8WUQ6Yq+dpDc/vx64GXgnSFJJ6l7rApiKPcX4LvAe8O/0AIwxb4vIf4AFwDLg1fQyuRhjnhMRF5gV5Mr1wFnGmKUi8qqIzMeecrs8UznskZBSbYIOt6GaRXBqanrQEEIppfLSU3xKKaVKkh5BKaWUKkl6BKWUUqokaYJSSilVkjRBKaWUKkmaoJRSSpUkTVBKKaVK0v8HyqD9yPrvyt4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax1 = plt.subplots()\n",
    "color1 = \"tab:blue\"\n",
    "ax1.semilogx(dp_end, n_0 / np.sum(n_0), drawstyle=\"steps-mid\", \\\n",
    "             label=\"before coagulation\", color = color1) \n",
    "ax1.set_ylabel(\"particle number in #/m^3\", color = color1)\n",
    "ax1.tick_params(axis=\"y\", labelcolor = color1)\n",
    "title_string = \"normalized particle size distribution\" + \"\\n\" + \\\n",
    "    \"total time = \" + \\\n",
    "    str(np.format_float_scientific(t_total, precision=0)) + \" s,  \" + \\\n",
    "    \"total particle number = \" + \\\n",
    "    str(round(np.sum(n)/np.sum(n_0),5) * 100) + \" %\"\n",
    "ax1.set_title(title_string)\n",
    "ax1.set_xlabel(\"particle diameter in m\")\n",
    "\n",
    "color2 = \"tab:red\"\n",
    "ax2 = ax1.twinx() \n",
    "ax2.semilogx(dp_end, n / np.sum(n), drawstyle=\"steps-mid\", \\\n",
    "             label=\"after coagulation\", color = color2) \n",
    "ax2.set_ylabel(\"particle number in #/m^3\", color = color2)\n",
    "ax2.tick_params(axis=\"y\", labelcolor = color2)\n",
    "\n",
    "# ask matplotlib for the plotted objects and their labels\n",
    "lines, labels = ax1.get_legend_handles_labels()\n",
    "lines2, labels2 = ax2.get_legend_handles_labels()\n",
    "ax2.legend(lines + lines2, labels + labels2, loc=0)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Things to try"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Change t_total from 2e-0 to 2e-2 to 2e2 or any other value.\n",
    "* What happens if t_total becomes too large and why ?\n",
    "* Change the start_bin from 4 to 8 or any other value.              \n",
    "* Change the value of the starting number concentration (n_0[start_bin] = 1e18). What happens if it is too large and why? \n",
    "* Start with two (or more) start bins such as  \n",
    "        start_bin = 4 \n",
    "        n_0[start_bin] = 1e17 \n",
    "        n_0[8] = 1e17 \n",
    "        \n",
    "* Write your own jupyter skript. Try to implement a $\\beta$ that is valid in the continuum regime."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Singh, M., Kaur, G., De Beer, T., Nopens, I. \"Solution of bivariate aggregation population balance equation: a comparative study\", Reaction Kinetics, Mechanisms and Catalysis, 123 (2) 385-401 (2018). \n",
    "\n",
    "Nopens, I., Beheydt, D., Vanrolleghem, P.A. \"Comparison and pitfalls of different discretised solution methods for population balance models: a simulation study\", Computers & Chemical Engineering, 29 (2) 367-377 (2005). \n",
    "\n",
    "Lister, J., Smit, D., Hounslow, M. \"Adjustable discretized population balance for growth and aggregation\", AIChE Journal, 41 (3) 591-603 (1995). \n",
    "\n",
    "Hounslow, M., Ryall, R., Marshall, V. \"A discretized population balance for nucleation, growth, and aggregation\", AIChE journal, 34 (11) 1821-1832 (1988). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": 4
}