2d_Cole-Cole_IP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# OUTDATED, the examples moved to the manual\n",
    "## See https://empymod.readthedocs.io/en/stable/examples\n",
    "\n",
    "----\n",
    "\n",
    "# Cole-Cole\n",
    "\n",
    "There are various different definitions of a Cole-Cole model, see for instance Tarasov and Titov (2013). We try a few different ones here, but you can supply your preferred version.\n",
    "\n",
    "The original Cole-Cole (1940) model was formulated for the complex dielectric permittivity. It is reformulated to conductivity to use it for IP,\n",
    "\n",
    "$$\n",
    "\\sigma(\\omega) = \\sigma_\\infty + \\frac{\\sigma_0 - \\sigma_\\infty}{1 + (i\\omega\\tau)^C}\\ . \\qquad\\qquad\\qquad (1)\n",
    "$$\n",
    "\n",
    "Another, similar model is given by Pelton et al. (1978),\n",
    "\n",
    "$$\n",
    "\\rho(\\omega) = \\rho_\\infty + \\frac{\\rho_0 - \\rho_\\infty}{1 + (i\\omega\\tau)^C}\\ . \\qquad\\qquad\\qquad (2)\n",
    "$$\n",
    "\n",
    "Equation (2) is just like equation (1), but replaces $\\sigma$ by $\\rho$. However, mathematically they are not the same. Substituting $\\rho = 1/\\sigma$ in the latter and resolving it for $\\sigma$ will not yield the former. Equation (2) is usually written in the following form, using the chargeability $m = (\\rho_0-\\rho_\\infty)/\\rho_0$,\n",
    "\n",
    "$$\n",
    "\\rho(\\omega) = \\rho_0 \\left[1 - m \\left(1- \\frac{1}{1 + (i\\omega\\tau)^C} \\right)\\right]\\ . \\quad (3)\n",
    "$$\n",
    "\n",
    "In all cases we add the part coming from the dielectric permittivity (displacement currents), even tough it usually doesn't matter in the frequency range of IP.\n",
    "\n",
    "#### Note: Requires `empymod v1.7.4` or newer.\n",
    "\n",
    "### References\n",
    "\n",
    "**Cole, K.S., and R.H. Cole, 1941**, Dispersion and adsorption in dielectrics. I. Alternating current characteristics; *Journal of Chemical Physics*, Volume 9, Pages 341-351, doi: [10.1063/1.1750906](https://doi.org/10.1063/1.1750906).\n",
    "\n",
    "**Pelton, W.H., S.H. Ward, P.G. Hallof, W.R. Sill, and P.H. Nelson, 1978**, Mineral discrimination and removal of inductive coupling with multifrequency IP, *Geophysics*, Volume 43, Pages 588-609, doi: [10.1190/1.1440839](https://doi.org/10.1190/1.1440839).\n",
    "\n",
    "**Tarasov, A., and K. Titov, 2013**, On the use of the Cole–Cole equations in spectral induced polarization; *Geophysical Journal International*, Volume 195, Issue 1, Pages 352-356, doi: [10.1093/gji/ggt251](https://doi.org/10.1093/gji/ggt251)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from empymod import bipole, versions\n",
    "from empymod.utils import _check_var"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot-style adjustments\n",
    "plt.style.use('ggplot')\n",
    "plt.rcParams['figure.dpi'] = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use `empymod` with a user-defined function to adjust $\\eta$ and $\\zeta$\n",
    "\n",
    "In principal it is always best to write your own modelling routine if you want to adjust something. Just copy `empymod.dipole` or `empymod.bipole` as a template, and modify it to your needs. Since `empymod v1.7.4`, however, there\n",
    "is a hook which allows you to modify $\\eta_h, \\eta_v, \\zeta_h$, and $\\zeta_v$ quite easily.\n",
    "\n",
    "The trick is to provide a dictionary (we name it `inp` here) instead of the resistivity vector in `res`. This dictionary, `inp`, has two mandatory plus optional entries:\n",
    "- `res`: the resistivity vector you would have provided normally (mandatory).\n",
    "- A function name, which has to be either or both of (mandatory)\n",
    "    - `func_eta`: To adjust `etaH` and `etaV`, or\n",
    "    - `func_zeta`: to adjust `zetaH` and `zetaV`.\n",
    "- In addition, you have to provide all parameters you use in `func_eta`/`func_zeta` and are not already provided to `empymod`. All additional parameters must have #layers elements.\n",
    "\n",
    "The functions `func_eta` and `func_zeta` must have the following characteristics:\n",
    "- The signature is `func(inp, p_dict)`, where\n",
    "    - `inp` is the dictionary you provide, and\n",
    "    - `p_dict` is a dictionary that contains all parameters so far calculated in empymod [`locals()`].\n",
    "- It must return `etaH, etaV` if `func_eta`, or `zetaH, zetaV` if `func_zeta`.\n",
    "\n",
    "Dummy example\n",
    "```\n",
    "def my_new_eta(inp, p_dict):\n",
    "    # Your calculations, using the parameters you provided\n",
    "    # in `inp` and the parameters from empymod in `p_dict`.\n",
    "    # In the example below, we provide, e.g., inp['tau']\n",
    "    return etaH, etaV\n",
    "```\n",
    "And then you call `empymod` with `res={'res': res-array, 'tau': tau, 'func_eta': my_new_eta}`.\n",
    "\n",
    "## Define the Cole-Cole model\n",
    "\n",
    "In this notebook we exploit this hook in empymod to calculate $\\eta_h$ and $\\eta_v$ with the Cole-Cole model. By default, $\\eta_h$ and $\\eta_v$ are calculated like this:\n",
    "\n",
    "$$\n",
    "    \\eta_h = \\frac{1}{\\rho} + j\\omega \\varepsilon_{r;h}\\varepsilon_0 \\ , \\qquad (4)\\\\\n",
    "    \\eta_v = \\frac{1}{\\rho \\lambda^2} + j\\omega\\varepsilon_{r;v}\\varepsilon_0 \\ . \\qquad (5)\n",
    "$$\n",
    "\n",
    "With this function we recalculate it. We replace the real part, the resistivity $\\rho$, in equations (4) and (5) by the complex, frequency-dependent Cole-Cole resistivity [$\\rho(\\omega)$], as given, for instance, in equations (1)-(3). Then we add back the imaginary part coming from thet dielectric permittivity (basically zero for low frequencies).\n",
    "\n",
    "Note that in this notebook we use this hook to model relaxation in the low frequency spectrum for IP measurements, replacing $\\rho$ by a frequency-dependent model $\\rho(\\omega)$. However, this could also be used to model dielectric phenomena in the high frequency spectrum, replacing $\\varepsilon_r$ by a frequency-dependent formula $\\varepsilon_r(\\omega)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cole_cole(inp, p_dict):\n",
    "    \"\"\"Cole and Cole (1941).\"\"\"   \n",
    "    \n",
    "    # Calculate complex conductivity from Cole-Cole\n",
    "    iotc = np.outer(2j*np.pi*p_dict['freq'], inp['tau'])**inp['c']\n",
    "    condH = inp['cond_8'] + (inp['cond_0']-inp['cond_8'])/(1+iotc)\n",
    "    condV = condH/p_dict['aniso']**2\n",
    "\n",
    "    # Add electric permittivity contribution\n",
    "    etaH = condH + 1j*p_dict['etaH'].imag\n",
    "    etaV = condV + 1j*p_dict['etaV'].imag\n",
    "\n",
    "    return etaH, etaV\n",
    "\n",
    "\n",
    "def pelton_et_al(inp, p_dict):\n",
    "    \"\"\" Pelton et al. (1978).\"\"\"\n",
    "\n",
    "    # Calculate complex resistivity from Pelton et al.\n",
    "    iotc = np.outer(2j*np.pi*p_dict['freq'], inp['tau'])**inp['c']\n",
    "    rhoH = inp['rho_0']*(1 - inp['m']*(1 - 1/(1 + iotc)))\n",
    "    rhoV = rhoH*p_dict['aniso']**2\n",
    "\n",
    "    # Add electric permittivity contribution\n",
    "    etaH = 1/rhoH + 1j*p_dict['etaH'].imag\n",
    "    etaV = 1/rhoV + 1j*p_dict['etaV'].imag\n",
    "\n",
    "    return etaH, etaV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "\n",
    "Two half-space model, air above earth:\n",
    "\n",
    "- x-directed sourcer at the surface\n",
    "- x-directed receiver, also at the surface, inline at an offset of 500 m. \n",
    "- Switch-on time-domain response\n",
    "- Isotropic\n",
    "- Model [air, subsurface]\n",
    "    - $\\rho_\\infty = 1/\\sigma_\\infty =$ [2e14, 10]\n",
    "    - $\\rho_0 = 1/\\sigma_0 =$ [2e14, 5]\n",
    "    - $\\tau =$ [0, 1]\n",
    "    - $c =$ [0, 0.5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAGKCAYAAADNBgH+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzs3XecXPV97//XOdO3V23Rane16h0J0QVIiGaDwYDDNdiOMbHjGOcm5tqx4wcmgO048b0mdmLin0uMwY5xbIILLqEJEL0jCTXU+66296nnfH9/nNVIKwlQ2dXMzr6fj8c8dubM2ZnvaFcz7/2Wz9cyxhhEREREspSd6QaIiIiIvBuFFREREclqCisiIiKS1RRWREREJKsprIiIiEhWU1gRERGRrKawIiIiIllNYUVERESymsKKiIiIZDWFFRHJOMuyWLp06ag9/k033YRlWezYsWPUnuO99PX1ceuttzJ58mQCgQCWZbFq1SoAkskkX/3qV5k+fTqhUAjLsvjtb3+bsbaKZBuFFRE5guM4/OhHP+LCCy+krKyMQCDAhAkTmD9/Pp/85Cd5+OGHR70Nd955J5Zl8fTTT4/6c50KX/rSl/jOd77D3Llz+fu//3vuuOMOqqurAfj2t7/NHXfcQU1NDV/4whe44447mDlzZoZbLJI9/JlugIhkF8dxuPLKK3nkkUcoKSnhiiuuoK6ujs7OTrZs2cLPfvYzNm7cyFVXXTViz7lhwwby8vJG7PGy0cMPP8z06dP5/e9/f9T7CgoKePzxxwkGgxlonUh2U1gRkWF+8Ytf8Mgjj7BgwQJWrlxJcXHxsPu7urp4/fXXR/Q5x0Mvwr59+7jgggve8b7y8nIFFZF3oGEgERnm+eefB7x5HocHFYDS0lIuvvhiAB599FEsy+IrX/nKsHMee+wxLMvCsix279497L7rr78ey7LYtm1b+tjhc1YaGxu56667AFi2bFn6sSzLGvZYg4ODfPOb32Tx4sUUFhZSUFDArFmz+Ju/+Rv2799/1Nf3gx/8gHnz5hEOh6mqquJTn/oU3d3dx/ivc9C+ffu45ZZbaGxsJBgMUllZyTXXXMOrr7467LylS5diWRbGGFauXJl+HUuXLk3Ppdm+fTs7d+5M39fY2Hjc7RHJZepZEZFhKisrAdi0adN7nnv++ecTDAZZsWIFX//619PHn3zyyfT1FStWcNNNNwFgjOHpp5+msbGRpqamd3zcz33uc/z2t79l5cqVfPzjHz/qh3dXVxfLli1j9erVzJw5k5tvvplgMMiWLVu49957ufbaa6mqqhr2PV/84hd59NFH+cAHPsCll17KU089xX/8x3+wadMmVq5c+Z6v94Bt27axZMkSmpubWb58OTfccAO7d+/mwQcf5I9//CMPPvggV199NeCFvqVLl3LXXXfR0NCQ/rdobGykpKSExsZGvvOd76RfN0BJSckxt0VkXDAiIodYtWqVCQQCxrIs85GPfMT86le/Mtu2bXvH888//3zj8/lMT09P+tjixYvN4sWLTWlpqfnoRz867LEBc/PNNw97DMBceOGFw47dcccdBjBPPfXUUZ/3hhtuMID5q7/6K+M4zrD7ent7TVdXV/r2xz/+cQOY+vp6s3PnzvTxZDJpzj//fAOYl1566R1f4+EuueQSA5h//ud/Hnb82WefNbZtm9LSUtPb2/uer/GAhoYG09DQcMzPLzLeaBhIRIZZsGABDzzwANXV1fz85z/n+uuvp6mpiYqKCq677jr+9Kc/DTt/+fLlOI6T7pno7u7mjTfe4NJLL2Xp0qWsWLEife6B68uXLz+pNra2tvLLX/6Smpoa7r77bmx7+FtZYWHhUXsn/uEf/oH6+vr0bb/fzyc+8QmAI4Zv3smePXt4/PHHaWho4POf//yw+5YsWcKHP/xhurq6+M1vfnO8L0tE3oHCiogc4UMf+hA7d+7k0Ucf5fbbb+fKK6/EcRx+/etfc8UVV/AXf/EXGGMAuOiii4CDQz9PP/00ruuyfPlyli9fTnNzMxs2bBh2zoHvOVGvvvoqrutywQUXHNcqosWLFx9xbNKkSYA3rHTAfffdx5133jnscmAJ9Ztvvgl4Q2B+/5Ej6Qfm87zxxhvH3C4ReXeasyIiRxUIBLj00ku59NJLAW9J80MPPcTNN9/Mvffey1VXXcXVV1/N2WefTX5+frrXZMWKFYTDYc4991xqa2vTx6ZNm8YzzzzD7Nmz0/VFTtSBCbETJ048ru872oThA4HDcZz0sfvuu++oc1iWLl1KT08PwDu+hpqaGoD0eSJy8tSzIiLHxOfzcf3113PrrbcCB4d0AoEAS5YsYe3atbS2trJixQrOO+88wuEwM2fOpLa2lieeeIJXXnmFvr6+k+5VgYMTUPfu3XvSj3U0Tz/9NMaYYZc777wTOBh4Wlpajvq9zc3Nw84TkZOnsCIix6WwsBAgPQwE3hwUYwwPPPAAGzZsGDYnZfny5Tz99NM8/vjj6dvHwufzAcN7PA4488wzsW2bZ599lsHBwRN+LSdi4cKFADz33HOkUqkj7n/qqacAWLRo0Sltl0guU1gRkWF+8Ytf8Pjjj+O67hH3tbS08KMf/QhgWIGzA70l//RP/wRwRFjp6enhe9/7HrZtH/MeQOXl5QBH1GkBb3n1hz/8Yfbt28eXvvSlYcEJoL+/f9SGYerq6rjkkkvYsWNHesnxAS+//DIPPPAApaWlXHPNNaPy/CLjkeasiMgwL7/8Mv/6r/9KdXU1S5YsYfLkyQBs376dP/7xj0SjUa6++mo+9KEPpb9n4cKFlJWV0draSnFxMaeffnr6vgPBpbW1lcWLFx9zDZFly5Zh2zZf/vKXeeuttygtLQVIF6C75557WLt2Lffccw8rVqzg0ksvJRgMsn37dh599FEefvjhUdsc8fvf/z7nnXcef/d3f8djjz3G4sWL03VWbNvmJz/5SboHSkROnsKKiAzz+c9/nmnTpvHEE0+wZs0aHn30UWKxGOXl5SxdupQbb7yRG2+8cVg12QM9Jr/+9a+58MIL00M44PVETJ8+nU2bNh3XfJVZs2Zx//33861vfYvvfe97xGIx4GBYKS0t5YUXXuA73/kOv/zlL/nhD3+Iz+dj0qRJ3HzzzcyePXuE/kWO1NTUxGuvvcbXv/51/vSnP/H0009TVFTE5Zdfzm233cYZZ5wxas8tMh5Z5vD+UxEREZEsojkrIiIiktUUVkRERCSrKayIiIhIVlNYERERkaymsCIiIiJZTWFFREREsprCioiIiGQ1hRURERHJaqpge4K6urqOuonZiaqsrKStrW3EHk9Gh35O2U8/o7FBP6exYaR/Tn6/P711xnF934i1YJxJpVIkk8kReawDZctTqdQRG7JJ9tDPKfvpZzQ26Oc0NmTTz0nDQCIiIpLVFFZEREQkqymsiIiISFZTWBEREZGsprAiIiIiWU1hRURERLKawoqIiIhkNYUVERERyWoKKyIiIpLVVMFWREaFMQaMC4Yjv3LofWbo9oHjh10/cN+B60N3ka6oecj5lkXKNpj2VozrHt6gg+cf8mX4fbzDOYfdf/j5Ry3u+R4VP4+pIuhxVg0diSKjp6JSqWWRGOzFtLePcmVUVcc9KZZF0k2CHch0SxRWRDLBuA4kEpCIQzLhXU8eckklIZXEJJOQSg3dPuSrkwTHGTqWAifl3T7w1XXAcTAHrg/dxnWHbruH3Ha94HDg+oH7jTn41Rx6+8D1Q46nvx5yyZDmjD2zHI/9mW6AHJOueYvgb+/KdDMUVkSOlUkmcDrbMft2Ywb6IDqAGRyA6CDEBoe+RiEWxcS9r8RjQ5c4JGJeOEnEvYAhx86yAAusA9c5ePvA9aHzLMvCYIYdO3j90Mc74kmGn8Nh5xz+PUc8xFEe86jP8y6O5fzjfcyjtWs0HOfT+Hw+HMcZnbbIcXiXH5wFdlEJ2fBTUliRccskE9DTlb6Y3m7o6/Eu/b2Yvh7o74OBPhjohUSCfaPRkEAQAgEIhIa+Dt32H3IJBMDnx/IHwOdL3/YuPvAfcj391Qe2D2z74HWfD8u2wTpwbOj6gXMs65BjPrCtg/cPu28oABy4bXHw9oH7rAOXQ+8/5L50+Djkfssadt06zg9my7Koqamhubk54xuvyTvTz2lssCyLiqGfU6YprEhOMqkUdLVDZxumvRU69kNXB6a7A7qGLoP9x//AlgWRfMgvgEiedz2SjxWJeNfDEe8SCkM4ghWKQCgEB74GQxAMQzDoXfcHvPAgIiLvSGFFxizjOtDeCs17MPv3QlszprUZWpuho82bR/Fe/AEoLvUuRSVYhcVQWAyFRVBQjFVQBPmFUFCIVVBIzeQptOzfr78GRUROIYUVyXrGGOjuhD07MHt2eF/37YT9+7zJqO/EH4CySqiYgFU+AcoqoLQCq6QcSsuhpMzrFTnGoQbLstQLIiKSAQorknVMdyfs2IzZsRmzcwvs2AL9vUc/2R+A6olYVROhqhYm1GBV1sCEaigqVbgQEckBCiuSUcZ1oXk3ZvN62LrB+9rReuSJtg1VE7HqGmHSZKzaBqip83pNbN8pb7eISC5xHUM8bojHXBJxQzxmCEdsamoy3TKPwoqccqatBbNhNWxc4309vNfEsqF2ElbjVGiYhtU4DeoasALBTDRXRGRMi8dcerod/H6LsgrvY991DS+tHCAec4nHDMnEkfPwqmoDLFh0qlt7dAorMupMKgmb12PWvIZZ8yq0HrYAOBiCphlYU2djTZ3lXY/kZaaxIiJjgOMYYlGXWNQQG3SJRl1ig0O3oy6zT4ukg0nL3iRrXosyocbPWRcUAGDbFt2dKZxDSj5ZFoTCFsGQTShsUVyaPb3WCisyKkw8Dmtfx7z2HGbdG17BtAN8Ppg8A2vWfKyZC6Bpulc/REREhjGuYdumONFBl+igGfrqDdW8m8F+l7IK73pegU1hsU1e/vA5fIvOzsfvh1DYCyeB4PDaRsdb52g0KazIiDHJJLz1mhdQ1rzqVW49oLAYa/5irPlnwOzTsMLqORGR8ckYQyJuGBxwycu3CYW9ELF/X5J1q6IUFftYfF6+d7IFm9bFjlr02vZBJGITzrMJRyzCEXvoYlFafvDjvbIqwNLLj/yDsHri2PkjUWFFTprZtRXz/ArMKyu9iq8HlE/AWnwe1qJzoXGaVuaIyLhhXK8XZKD/4GWw32Wg32FwwE0Pv5x2Vh6TGr35eJYFA30uh75VWpZFw9QQAJE8e+hiEcmzj+gJyWUKK3JCTDyGeelpzMr/gd3bD95RUoZ1xvlYZ5zvBZRx8h9JRGTjW1F6ux0G+lwGB1wO3/j7cOGIhXEPDueUlPs4Z1k+efnD54rMXhAZjeaOKQorclxMazPmqT9hnn8CogPeQb8fa8FZWOddDHNO01JiEckprmtwXfD7vT++erpSrHktis9vce6ygvR5LXuT9PUcTCi27c0XyS+wySvwDX21yc+3ieTb+HzD/5gLBm0qJqgH+mgUVuSYmO2bcf/0IKx+GQ6Ump9Qg7Xs/VjnXISVX5jZBoqInCTjGgYGXPp6HPp6hr72ej0l0+eEmTY7DIDPZ9Hd6eDzefNPDvQgN00P4bqQX2CTX+gjErGw7LHTuxxNuuzuibOrJ87O7jgTi0LcnCWFVhRW5B0ZY2DTWi+krF918I65i7AvuhLmLNI8FBEZk+Jxl95uh95uh75ul96hYOI6Rz9/oO9gj0legc3i8/IoKBzei1zfFBrNJo8Y1xia+5Ls6I4xvyqfwpD3On69voNfre1Inze/Oo+bM9XIwyisyFGZLetxH/opbFnvHbBtrLOWYr3vOqyaSZltnIjIMTLGEB00w5btvv7iAPt2JY96vu2DwiIfhcU2hcW+9PVI3sHvt22LmrqxUaQy4bjs7I6zrTPOtq4Y27vi7OyOEUt5PeT/sLSO0yd6Q1n1xSFKwz4mlYRoKA4xvSJ75soorMgwZu9O3N/8DFa/4h3wB7CWXIJ12TVYFVWZbZyIyLswxpBMGoJBL1g4KcPjv+8lmTBc9sEigiHv+IHgkVdgU1Tso6jEpqjER2Gxj/x8e0wN3RwqnnLZ0R2nIs9PeZ63LPmZHb1896WWI84N+izqi4f3BC1pKOT8xqL07WxaIKGwIoC3eaD5zc8wLz7pzUmxbS+kfODD3i7FIiJZJhZ16epI0d3p0N3p0NPpUFTqS0969fm9QmdOyjDQ56bDytSZIabPDuMPZM+H8fFyXMPunjibOmJs7oiyuSPGzu44roG/XFzFFTNKAZhSFqYo5KOpLExTaYjJpWEml4aoLQziOyyUZVM4OZzCyjhnHAfz1B8wv3sAYlHv4OnnYn/wo1jVdZltnIjIEMcx9HQ5dLWn6Opw6OpIEYseWcW1r8cZNun1nKUFhMMW9iErbw6ElrHEcU06XGztjPHlx3YSd458/cVhH6lDlkM3loT46XVTszqIHAuFlXHMbFqL+8APYO9O70DjNOwb/hKraUZmGyYiAvT1OuzamvB6T7oczOF1SywoLLIpLfNTUu6jpMwbyjn0g/nwEvNjwYFekw1tUTa2RdnQHuXsugJuPt0biq8tDJJ0DRG/zbTy8NAlwrSKMOURf9aWzD8ZCivjkIkOYh68F/PsY96B/EKsa//cG/bR6h4RyYDBfof21hRFJT5KyryPpnjUZdumePqcYMiitMJHWbmfknI/JaW+MT2Uc6ikY/jdxk7Wtw6ysS3KQHJ4Mnu7/eD2JZGAzfc+0ERVQQA7R8LIe1FYGWfM+lW4938XOtsAsC64DOuaj2EVFL3Hd4qIjAxjDN2dcVzXcOCzdsvGODu3JmiaHkqHlZIyPw1TgpRV+Cmt8JGXb+dET0HCcdnUHqMnluK8Bu+912/D7zd20h3z1k6H/TYzKsLMqowwszKPGRXhYY9RUzg2ViONFIWVccLEBjH/fR9m5SPegYoq7Jv+BmvGvMw2TETGhQM9J+37U7S3pojHullycSGl5V6Nj4oqP709DvmFB3t3/QGL+YvH/qanKdewuSPKWy2DrNnv9ZwkXUNJ2Me59YVYlrfHzzWzy7AtizkT8mgsCR0xAXY8U1gZB8yubbg/+L/Qug/Aqzp77cexwtmzhl5Ecksi7tLemqKtxQsogwPDhzV8PouBficdVmonBamdlHu9BT9+fT+Pbekhlhr++kvDPuZW5RFNueQFvH+DD87Syst3orCSw4wxmJX/g/nljyGVhLIK7E98Dmvm/Ew3TURyjOsaujocWpuTtLWk6OkaXgrWsryN+iom+KmsCjBrbh1tbfu9Stk5oDeWYlXLIG82D/DpM6oI+70eIr9tEUu5FIZ8zKvKY35VHvOq8phYFMyJIa1TRWElR5nBAcxP78G8/rx3YMGZ2J/4W+3hIyKj4qWVA3S0poYdKyyyqagOUFnlp7zSn54Ma1kWfv/YnszvuIYtnTFe29vPm80DbOmIcSB2nd9QyKJar9bL5dNKOL+hiMbS0LiZDDsaFFZykGneg3vP171hH58P67qbsC6+SileRE6a6xjeXhejrSXFucsK0gGkrMJHX49DZZWfyuoAldV+wpGxHUjeyWt7+/nXF5vpjQ/vPWosCbGwJp8JBYH0saqC3BvaygSFlRxj1r6O+8NvQXQAyiqxP/1F1U0RkROWTBr6ehzKKryPC8uGfbuTDPa7tO1PpvfImTorzIy54Zz7o2hvb4JX9/bRWBLmtJp8AKoKAvTGHfIDNqfV5LOoNp+FNfnpEvcy8hRWcoQxBrPiYcyvfgLGhamzsT/z91hFJZlumoiMMbGoS8veJC17k7S3pvD54LKri7F93qqV6XPCWED5hIMfIX5/boQU1xg2d8R4eXcfL+/pZ09vAoDz6gvTYaWuKMg/X1LPtIoIfq3YOSUUVnKASaUwD3w/XeTNOm851kduwQoo5YvIsRnoc2jek6R5T5LuzuHDG6F8m2jUJb/AW7UyqTH3hjYc1/Cj1/bz0p5+uqIH5974bZgzIS8dVMCbczNrwthfUj2WKKyMcSYex/3BN+Gt18Cysf7sE5qfIiLHpK9nKKDsTtDbc9jS2nIf1RMDVE0MUFjky1ALR0/ScdnWFWdGhVfCwWdbbOmM0RVNEfHbnD4xn7PqCjm9Np/8YO69/rFGYWUMMwP9uPd8DbZsgEDQm5+y4MxMN0tEspgxhs3r4+zblaCv92BAsSxvWKemLkD1xEBOTo5NOobVLQM8t7OXV/b0E3dc7r9uGgVDYeTG+RUAzKvKI+DLvdc/limsjFGmuwP3O3d6mxBG8rH/9+1Y02ZnulkikoXiMZdQ2PvwtSyLtpYkfb0ulg2VVX5qJwWoqg2Myd2I34vjGt7aP8gzO3p5aU8fA4mDAa004qe5L8G0cq935cByY8k+CitjkGlrwb37K9DRCsWl2J+7E6tucqabJSJZJpU0vPh0Pz1dDpdcXURoKIxMnRUmETdUT/QTCOZeQDnUI5u7+eFr+9O3SyN+zq0vZEl9ITMrI6p9MkYorIwxpn0/7rdu8zYirKzGvvWrWJXVmW6WiGSBVNLQ3eVQMbRKxx+wcF0wQFe7Q/VEL5hU1ebm5Ps9PXFW7uilqSzMOZO8AphnTyrgl2vbObuukAsai5g9QQFlLFJYGUNMR+vBoFI1EfsL/4hVUpbpZolIBhnX0LY/xZ4dCZr3JgG49OpiAkPF2k47M0I4YqeHgXJNbyzFszv7eGp7D5s7YgDMr8pLh5XyvAA/uWaqNgUc4xRWxgjT0eYFlY5WmFCD/fmvK6iIjGN9PQ67dyTYsyNBPHZwf538QpvBfpfiUm/SaHFpbr7Nv7y7jxXbenhtbz/O0Mu3LVhUk8/SycXDzlVQGfty87c4x5jOdty7b4P2/d7Qz+f/EatUu3OKjDfJhGHvrgS7tyeG1UIJBC0m1geoawxSUuYbF6ULHn67i7X7BwGYUhZi2eRizm8soiSsj7VcpJ9qljMDfbjfuQPaWryg8oV/xCqryHSzROQUMcbQ2eawa1ucfXuSuEMZxbJgQo2fSZODTKgJ4PPlZkAZTDo8s6OXJ7f18OUL6iiNeB9bV84oZWpZmIuaimkoCWW4lTLaFFaymEnEvQ0Jm3dDSbk39FNWmelmicgptOqVQfbsSKZvFxbZTJocpK4xmLPzUIwxbGyP8tiWHp7f2Ut8aJzn6e09XDPb61U+Z1Jhel6K5D6FlSxlHAf3R9/yCr7l5XvLk8snZLpZIjKKjDG0t6YoKvalg8iE6gDNe5JMrA9S35TbwzzRpMsTW7t5bEs3u3oS6eOTioNcMqWEpZOLMtg6ySSFlSxkjME88H1Y9TL4A9if/QrWxIZMN0tERtkbLw6yb3eSmfPDTJsVBqCmLkBVbTH+QG4GlEMlXcP9b7aRdA0hn8WShiIunVrCjIrc281Zjo/CShYyf/wl5plHwbKxP/UFrOlzMt0kERlhxhi6Ox0KCu10YbbKaj+tzUnMIdv02D4LOwe3poklHR7b0s3m9ii3nOXViioK+bh2ThklYT8XNhZpTx5JU1jJMua15zC/ewAA68ZPYy06J8MtEpGR5DiGfbsSbN+coKfLYc7CCE3TvQmiExuC1E4K5nQvSnNfgv/Z1M2T2zfTF/d2N37/9BIaS72epBvna16eHElhJYuYXVtxf/IdAKxLP4i99H2ZbZCIjJjBAZedW+Ps3JogmfAmjNo2JOIHu1F8PgtysDPBNYZVzQP84e0u3tg3wIGqMFUFAS6fVkJFXm5W1JWRo7CSJZyuDpx7vg6JBMxdhHXdxzPdJBE5ScYYOtsdtm+Ke9Vlhz6lI3kWjVNDTGoKpvfryWUv7+7nn5/dm759em0+HztnCo3hBLnbhyQjSWElC5hkkva7b4POdq+M/qe+gJWLg9Qi44TrGPbtTrJtU5yeroPF28on+Jk8LUhVbQA7h6uq7u9P0DaQYm5VHgCLJxZQVxRkYW0+759WysTiEDU1FTQ3N2OMeY9HE1FYyThjDOY/v4ezYQ1E8rH/+jasPG1TLjIWpVKG7ZvibN8cT5fAt31Q1xBk8rQQRSW5/UfIxrYov9vYyUu7+6jI8/P9q6bgsy0CPovvXjlZGwjKCVNYybTN6zDPPwG2jf3pL2JV12W6RSJygiwLtm2Kk4gbwhFvqKd+Sm4P9Tiu4aU9ffxuQydvt8fSx2sLg/TFHUqGKs4qqMjJUFjJMGv6XOxP/C1FwQD9cxepS1RkDOnq8HY7nrsogmVZ+HwWs+aHsWyLiZMC2DlaAv+AN/b184NX99PS71XY9dsWSycXcdXMMpXAlxGlsJIF7PMuprCmhv7m5kw3RUSOUSpleHnlAMmkobI6QPVEb0VLfdP4+ZDOC/ho6U9SGPLxvmklvH96aXrvHpGRpN8qEZFj4LqG/fuSVE8MYFkWfr/F5OlBogOGgqLcHeY5YH9/gt9u6CTst/n4Qm/rj5mVEf5uSS1nTCwg5M/9fwPJHIUVEZF3kUoZdm1LsO3tGNFBw9kX5lNZ7fWizJgbyXDrRt+OrhgPre/kuZ29uAaCPotrZ5dTGPImCy9p0H49MvoUVkREjiKZcNm+JcH2oQmzAKGwlS7olus2tA3y32s7eG3fQPrYwpp8rp1dRkFQvShyaimsiIgcIh5z2bYpzo4tcVLevFHy8m2mzAwxaXLQqzKb4/74dhc/fG0/ALYF50wq5ENzymkqC2e4ZTJeKayIiACxqMuWjXF2bo3jDtVxKyyymTo7TO2k3C7iZoxhIOFSMDS0c/akAn62qo0lDYVcO7uc2qJghlso453CioiMa4MDLls3xti1LYE7tE1PcamPabND6cm0uco1hhd39fHgug5Kw37uuGgSAOV5Ae69dgp5gdwuYidjh8KKiIxbibjLU//Tm+5JKa3wMX12mMpqf06HFMc1PL+rj1++1c6e3gQALf4kXdFUeumxgopkk3EdVv7whz/w5JNPYoxh3rx5fOITn8jpNygR8QJKcKiibDBkU1MXIBY1TJ8Torwy90PKszt7+dXaDvYOhZT8gM0HZpZy5Yyy9AofkWwzbsNKb28vjz76KHfffTd+v5877riDzZs3M3369Ew3TURGgXENa16PsmdHggsvK6SgyPtgXnBG3rgajU5cAAAgAElEQVSYNAuwckcv//qiV3yyIGhz9cwyrphRSn5QIUWy27gNKwCO45BMetP9U6kURUWqFyCSqyzbIh5zcV1o2Ztk6lBYyeWg4riG9sEkVQXeBNnzGwp5eGMnS+qLeP+MEg31yJhx3GHlN7/5Da+88gp79+4lGAwyffp0PvrRj1JbWztijVq/fj0PP/ww27dvp6uriy984QuceeaZR5z36KOP8vDDD9Pd3U1dXR033XQTs2bNOqbnKCoq4gMf+AC33HILtm1zySWXUF1dPWKvQUQyKxZ12bIhxpSZYSJ53rDPrPkRps40lFXm9t9prjE8v7OP/3qrHdcY7rmyaWj3Y5tvv68xp4e6JDcd9//Y9evXc9lllzFlyhQcx+G//uu/+PrXv86//Mu/EA4fuQZ/48aNTJ06Fb9/+FPt3buX/Px8SkpKjvieeDxOY2Mjy5Yt4+677z5qO1544QXuu+8+PvnJTzJjxgyeeOIJvvGNb/Dtb3+biooKAL70pS+RSqWO+N7bbruNYDDIG2+8wb//+78TDAb5xje+wfr165k9e/bx/pOISBaJDqZY++YgO7Z4S5BdF+YvzgOgsDi3exKMMby0u59frGlnZ08c8IZ79vQm0hsLKqjIWHTcYeW2224bdvuWW27hk5/8JNu2bTvig951XX784x9TU1PD5z73OWzb++tm37593HXXXVxxxRVcffXVRzzHwoULWbhw4bu24w9/+AMXXXQRy5cvB+Cmm25i9erVPPbYY9x4440AfPOb33zH73/xxRepqqqioKAAgEWLFrF582aFFZExKpFw2bYxwfYtm0klvSqzpRU+ausDGW7Z6DPG8Pq+AX6+uo1tXV5IyQvYXD2rjA9oTorkgJOumTw4OAiQ/tAf9uC2zZe//GW2b9/OPffcg+u6tLS08NWvfpXFixcfNagci1QqxbZt21iwYMGw4/Pnz+ftt98+pscoLy9n06ZNJBIJXNdl3bp17zqU9cgjj3Drrbe+Y0+PiGRGKmnYtD7Gij/0snlDjFTSUFLm46wL8jnvogIqJuR+WNnUEeNrT+9hW1ecsN/mz+aU86Orp/DheRUKKpITTmrg1hjD/fffz8yZM6mvrz/qOWVlZdxxxx3ccccd/Nu//RubNm1i7ty5fOpTnzrh5+3t7cV1XYqLi4cdLy4upru7+5geY/r06SxcuJAvfelLWJbF3LlzWbx48Tuef/nll3P55ZefcJtFZGQ5jmHnljibNxzcu6ew2ObcCycSjPRluHWjrzuaomSoJsr08jCn1+ZTXxzi2tllFIVze06OjD8n9Rv94x//mF27dvHVr371Xc+rqKjgs5/9LHfeeSdVVVV85jOfGZFx06M9xvE87g033MANN9xw0u0QkVPHdQ17diTYtM7bBRkgv8BmxtwwExuC1NYW0tzcjzG5ueHgjq4YD6xp5639g/zg6ikUhXxYlsXtS+s0H0Vy1gmHlXvvvZfXX3+du+66i/Ly8nc9t7u7mx/+8IecfvrpbN26lfvvv5+bb775RJ+aoqIibNs+ohelp6fniN4WEcktr78wSMter+RAOGIxfU6YSZOD2LaV0x/WzX0JfrGmnWd29GLwNhhc0zLAkgav5EIuv3aR4w4rxhjuvfdeXnnlFe68804mTJjwruf39vbyta99jYkTJ/J//s//obm5mbvuugu/38+f//mfn1ij/X6amppYs2bNsCXNa9as4YwzzjihxxSR7GWMSX8Y1zUG6GhLMW1WiMapIXz+3P6Q7oym+NVb7Ty2pRtnqLPovPpCbpxfQV1xKLONEzlFjjus/PjHP+a5557ji1/8IpFIJN27kZeXRzA4fGdO13X5xje+QUVFBbfeeis+n4+6ujpuv/127rrrLsrKyrjyyiuPeI5YLEZLS0v6dmtrKzt27KCgoCC9LPnKK6/ku9/9Lk1NTUyfPp0nnniC9vZ2LrnkkuN9SSKSpXq6HDasiVJZ7WfKDK80QvXEABVVAQKB3A4pAINJh8/+fhuDSW+HxYU1+Xx0QSVTy48sEyGSy447rDz22GMA3HnnncOO33LLLSxdunTYMdu2ufHGG5k5c+awOiv19fXcfvvtR11BBLB161buuuuu9O2f/vSnAFx44YV89rOfBeDcc8+lr6+Phx56iK6uLiZNmsSXv/xlKisrj/cliUiW6u5M0daSorfbYfK0UHqoJ5DDC3wc1+CzvSCWF/BxQWMR27tifOy0SuZV5We4dSKZYZlcnYU2ytra2tKl+k+WZVnU1NTQ3Nycs5MCc4F+TqMvHncZ7HcpLff+uHFdw4bVMRqnBckveO8luGP5Z+S4hqe29/Bfa9r5ytI6Gku93pN4yiXoy635OGP55zSejMbPKRAInFCngta3iUjGpVKG7ZvibNkYIxCwWPb+Inw+C9u2mLMwkunmjSpjDK/s6eenq9rYM7QT8u82dvK353h1n0L+ky6HJTLmKayISMYY17B7R4K318aIRb2/3PLybeJRl7xj6EkZ6za0DnL/qjY2tEUBKAzafGhuOe+fXprhlolkF4UVETnljDG0tqTYsDpKX483eTSSZzFzXoSJDYGcGvJ4J999qZkntvYAEPRZXDWzjGtml1GgirMiR1BYEZFTqqcrxfrVMdr3e5uMBoIW02YPLUP25X5IOaChJIRtwfKmYm6YX0F5Xg7PGhY5SQorInJKRAddNr4VZc8Ob2K6bUPjtBDTZocIBnN7XsZg0uE36zuZVh7mzLpCAN43rYSFNflMUq0UkfeksCIioyqVNGzeEGPbpjiu4x2bWB9g5rxwzs9LSbmGx7Z0819r2umJO9QWBllUW4Dftgj4bAUVkWOksCIio8pxDDu2eEGlrNLHnAURSspz+63HGMNLe/r56Ztt7OvzVvjUFgb42GkVjKORLpERk9vvGCJyyhlj6O5wKK3w3l5CYZu5CyP4AxbVE3N/8uzWzhg/em1/eoVPccjH/5pXwWXTSvDbuf3aRUaLwoqIjBhjDK8+N8D+fSnOuiCfCTXepNFJk8fPcEdXNMWGtihBn8XVM8u4dk4ZeYHcHu4SGW0KKyIyYizLIr/Qh2Wn6O9zmVCT6RaNvt64w46uGPOrvVL4p9fm8/GFlVzQWESFVviIjAiFFRE5YY5j2LYpTmWVn5Iy7+1k+uwQjVOPrTz+WJZwXP74dhcPru0A4PtXNVEU9mNZFtfOLs9w60Ryi8KKiBw3Ywz7difZsDpKdNDQWunj3GUF3iaDQZtA8L0fY6xyjeG5nX38bFUrrQNerZjGkhDdMYeisN5SRUaD/meJyHHp7kix9s0oXR3eOuRwxKKhaXzMSVnXOshP3mhlc0cMgLKIn48uqGDp5OL0TskiMvIUVkTkmEQHXTauibJnp1fUzeeDqbPCNM0I4ffn/gd1+2CSrzyxC9dA2G9z3ewyrp5Vpo0GRU4BhRUReVeplGHrxjhbN8Zwhoq61TUGmDkvQiQvtz+o4yk3HUYq8gK8b1oJjoEb5lVQEtHbp8ipov9tInJUxhj27vLmpRzYEbmswsechZH0ZNpclXBc/rCxi1+v7+AfL2mgocQb5vrU4qqcrxMjko1y+x1HRE5Ib7fDmtcG0/NSInkWsxdEqJmU20XdjjZ59tEt3fzl4iqAnH7tItlMYUVEjmAMdHU4+PwwbVaYpukhfDk+L+XwybPlET8fPa2SpZOLMtwyEVFYERGclKGzPUVltVfErLjUx2ln5lFZ7Sccye15KQDffamZJ7b2AJo8K5KNFFZExrlE3OWZx/uJDbpceHkhhUVeMbdJk3O4WMph6otD2BZcMqWEG+ZXUKrJsyJZRf8jRca5YMimqNjGuIZ41E2HlVx1oPJsfXGI0ycWAPD+6aUsrM2nvnh81IsRGWsUVkTGmXjcZfO6GFNnhdNDPPMX5+EPWDldL8UYw7M7+/jZqjZaB5JMLAqyoCYfv20R8FkKKiJZTGFFZJxwXcOOLQk2rY2RTBpSKTjtzDyAnJ+XcrTKs9fNLiN3o5lIblFYERkHWluSrHszSn+vC0BRiY9Jjbk/J6W5L8F9b7by0u5+QJNnRcYqhRWRHDbQ77DuzSj793k1Q4Ihi5nzwtRPDmKNg71s9vYmeGl3P7YFl04tUeVZkTFK/2tFclAqadi8Ica2t+O4LlgWNE4LMWNOiEAwd3sU4imXtft6KB/KYafX5nP93HLObyzSnBSRMUxhRSSHGGPYuzPJhjUHS+RXVvuZszCS06t8XGN4ensvP1/dRsyB7181mcKgD8uy+MiCykw3T0ROksKKSI7o6Urx1uvRdIn8vHybOQsjVNX6c7pM/KrmAe57s5XtXXEAqgpDtPQlKSzP3XAmMt4orIjkiIF+d3iJ/BkhfL7cDSk7umLc/2YbbzQPAJAfsPnQ3Ao+tXQWnW2tGGMy3EIRGSkKKyJjlOsa+ntdikq8HoSaugAz54WpawwSycvdeSkA3dEUn39kJynX4LfhfdNKuX5uOcWRACG/elREco3CisgYNDjg8sqz/cSihouuKCQYtLEsi2mzw5lu2qhJOi4BnxfCSiJ+ljcV059w+NhpldQU5v4ybJHxTGFFZAwKRywYGuXo63Epr8zdnpSkY3hkcxcPruvga8vraSjxVvV8+owqfONg+bWIKKyIjAmOY9i1NUH9lCA+n4VtW5x+bj6hsEUwlJtBxTWG53b28fPVbbT0JwH406YuPnNmNYCCisg4orAikuX27/Oqzw70uziOYeosb6insDh352asaRngvjfb2NrplccvDfu4YX4lF08pznDLRCQTFFZEstTh1WdDYYu8/NzsRTnUvzy/j5U7egGvPP41s8v44KwywiqPLzJuKayIZJlUyrBlQ4ytG4eqz9rQND3E9Nlh/IHcH/qYUhbmuZ29XD6thOvnVVAS1tuUyHindwGRLGGMoWWvN+QTHfRmz1ZU+Zm7KHerz/bEUjy4toP51XmcWVcIwPunl3BmXYFW+IhImsKKSBbo73NY+0aUthZvyCecZzHntAg1dYGcrD4bTbo8vLGT36zvJJpyebN5gNNrC/DZFgGfraAiIsMorIhk0OFDPrYNTTNCTJsdxu/PvZCSdAyPbenmV2vb6Y552wI0lYb4+MIJWt0jIu9IYUUkg3ZujbN5vbenTWW1N+RTUJibQz5v7OvnB6/uTy9Dri4I8JEFlSxpKMTOwd4jERk5Cisip5gxJj200zg1RGtzisapQaon5uaQzwGugZb+JCVhHx+eV8HFU0oI5PDeRSIychRWRE4RxzFs3Rhn/74k5y0vwLYtfD6Lc5YWZLppo2J96yDtgykuaCwC4PTafP732dUsaSjSMmQROS4KKyKnSCpl2LYpTjJhaNmTpLY+NyeRbuuM8Z+r23h93wD5QZtFNfkUhHxYlsXFU0oy3TwRGYMUVkRGUTzuEgxaWJZFKGQz7/QIADWTAhlu2cjb0xPngTXtPL+rDwCfBUvqi3CNyXDLRGSsU1gRGQWua9ixOc7GtTEWnpVHTZ3XizIxB3tTOgaT/Ofqdp7e3oNrwAIuaCzihvkVWoIsIiNCYUVkhHV3pFj9WpTebm9p7t5dyXRYyUWxlEkHlTPrCvjI/AoaS8OZbpaI5BCFFZERkkwYNr4VZceWBACBgMWsBWHqm3IrqPTEUqxqHuDCyd6mghOLgnxi0QRmVESYURHJcOtEJBcprIicJGMMzXuSrH0jSjzmzc+Y2BBgzmkRQuHcWfXSH3f47YZOfv92J/GUYXJZmPriEABXzSzLcOtEJJcprIichMEBl7VvDKZ3Rs4vsJm3OEJlVe5MoB1MOvx+Yxe/29DJQNIFYEpZiHjKzXDLRGS8UFgROQHGNWwfmkDrpLydkafO9Mrk+3Kk0Fk85fL7jV38dkMHfQkvmDQUh7hxQQVn1RXkdAE7EckuCisix6m/z+GNFwfo7vQm0JZV+Ji/OI/C4twqk+8Ykw4qE4uCfHhehUrji0hGKKyIHCe/36K/z8EfgNkLItQ3BXOilyGecnluZy8XNRVjWRZ5AR83LZqA37Y4v6FIGw2KSMYorIgcg/4+h8Ii779LOGKz+Nx8Cot9hCNjfwJtPOXy6JZuHlrXQXfMoSDk46y6QgBVnBWRrKCwIvIujDGsezPK9s0JzrqggJoa73hl9difQBtPuTw2FFK6Yt6Q1oT8sf+6RCT3KKyIvAvLstJDPJ3tqQy3ZmQ4ruFPm7oOCyl+PjSngouairUTsohkHYUVkcMkE4ZUyhDJ84Z4ZswLU1Xrp7I6N4q72Ras2NZDV8xhQr6fP5tbwbLJCikikr0UVkQO0dqSZPUrg0TybM67qADLtvD7LSrGcN2U2NBwzyVTSogEbCzL4s9Pq6R9MKWQIiJjgsKKCJBMGtavirJrm1cq37YN0aghL3/sfpAPJh3+tKmb323opDfukHIN184uB2BRbUGGWycicuwUVmTca29NserlAaKDXqn8ydOCzJwfwe8fm0FlIOHwx7e7eHhjZ7qYW3VBQJNnRWTMUliRcctxDBvXxNi2KQ5AJN/mtDMjVEwYmx/qxhh++VYHD7/dycBQSKktDPBncyu4sFF1UkRk7FJYkXGppyvFGy8N0t/rfajXNwWZc1oEf2DsfqBblsWunjgDCZe6oiDXzy1niYq5iUgOUFiRccW4hq1ve3v6GBdCYYsFZ+RRVTv2elM6oyl+u76D900vpabQW6n04fkVnNdQyDmTVBZfRHKHwoqMG9FBlzdfHqSj1auXUj0xwPzFEULhsVWFtm0gyUPrOnhiaw9J1zCYdPnrs71qdfXFIeqLQxluoYjIyFJYkXFj68YYHa0pfH6YuzDCpMlja0+f5r4E/72ug6e29eB4c4GZVRnhvIaizDZMRGSUKazIuDFzfoR43DBjbpiCwrG1Q/IPXm3hkc3duEMhZX51HtfPLWfuhLwxFbhERE6EworkrL4eh51b48xZGMGyvOJup5+Tn+lmnZCCoA/XwOLafP5sbgUzKyOZbpKIyCmjsCI5KZU0PP9kP8mEIb/Ax+TpY2cex4a2QR5c28GVM0rTxduumlnGOZMKaSoLZ7h1IiKnnsKK5CR/wGLmvDAte5PU1mf/Sh9jDGv2eyHlrf2DgLcr8oGwUhjyURgaW0NXIiIjRWFFckZfr4PrGIpLvV/rhilBGqZk9yRaYwyv7xvgV2vbebs9BoDfhmWTi7luTnmGWycikh0UViQn7N2VYPWrgwSDFhdcWkgwZGd1SDng7uf38ezOPgCCPotLppZwzawyKlUaX0QkTWFFxjTX9TYg3L7Z24Awr8yHMRlu1LtwXINrSO90vHhiAa/u7ed900q5elYZpRH9lxQROZzeGWXMisdcXn9hgI42B4Cps0LMmBvGzsLy8inX8PT2Hv57XQdXTC/lAzPLADi/oYhFtQUUaT6KiMg7UliRMam7M8Wrzw8QGzT4/bDw7HyqJ2bf0EnScXliaw+/Xt9B64BXOXfFth6unFGKZVn4bEtBRUTkPSisyJize0eCNa8O4rqQX2hzxpJ8Couy6wM/nnJ5fGs3v17XSUfUCyklYR/XzC7jsqmlY2I+jYhItlBYkTHDGMPGt2Js2RAHoKrWz8Kz8gkEs++D//uv7ufJbT0AlEf8XDunjEumlBDyj619iEREsoHCiowJqZThzZcGadmbBLz5KTPnhbOmhyKeckk6hoKhIZ0rppfyVssA180p5+IpxQR8CikiIidKYUWyXjzu8tLTA/R2O9g2LDgjj7rGYKabBXgh5ZHN3Ty0voPz6gv59BnVAEwtD/ODq6fgy8LJviIiY43CimS9QMAiFLYIhizOOC+fssrM/9oeGlJ6Yt5qpDUtg6Rcg38ooCioiIiMjMy/64u8B9v2NiBMJl3y8jM7kTaecvjDxk4eXNtO11BIqSoIcP3ccpZOLk4HFRERGTkKK5KVdm9P0NOVSu+YHAhaBIKZX/Fz/8s7+dFr+wGYkO/n+rkVLGtSSBERGU0KK5J1+vscVr06CAbKJ/ipqcvc/BTHNfQnHIrD3n+VP1tYx/+s3ceVM0q5eEpJuhKtiIiMHoUVyToFhT7mnBYhFnUzVujNGMPLe/r5+eo2SiJ+vra8HoDSvCD//oGmjLRJRGS8UliRrOC6hmTSEAp5S3ybpocy1pa1+we5/81WNnV4uyB3RFO0DyapzPd6eCzLwmTzBkQiIjlGYUUyznEMb7w4yECfw7kXFRAMZaYmyY6uGD9d1cbr+wYACPksrppZxgdnl1GQBfNlRETGK4UVyahU0vDq8wO0709h29Db7VBRderDyqrmAe58cjcG8Flw6dQS/te8Cu2CLCKSBfROLBmTShpeeqafrnYHnx/OXJJPRdWpm6NijElXwJ0zIY+awiCTS0N8dEEltUXZUXROREQUViRDkknDyyv76epwCAQszrown9LyU/Pr6LiGFdt6eGpbD19dXk/AZxHwWfzL+xqJBFQWX0Qk2yisyCmXTBhefmYoqAQtzr4wn5KyU/OruG7/ID96fT/bu7zNEJ/a3sOlU0sAFFRERLKUwoqcUsmE4aWV/XR3ntqg0tqf5L43W3l+Vx8A+UGbG+ZVsGxy8ag/t4iInByFFTllUsnhQeWcpfkUl47ur6DjGh5a18GD6zpIOAbbgsumlnDj/AqKwvr1FxEZC/RuLaeEkzK88tzAKQ0qALYFa1sHSTiGuVV5fPL0CUwuDY/684qIyMhRWJFR5zqG114YoKM1hd8PZ18wukGlfTBJ2G9TEPRhWRZ/eUYV2zrjnN9QmF79IyIiY4dmFMqoSyQMA30utg/OPL+AklFa9eMawx/f7uKvf7+d/1zVlj5eVxTigsYiBRURkTFKPSsy6sIRm3MvKqC/z6F8wuj8yu3qjnPPyy283R4FYEd3nJRrtBuyiEgOUFiRUdPf51BQ6JWpD0dswpGR78hLOob/XtfOf6/rIOVCxG/zsdMqed/0Emz1pIiI5ASFFRkVO7bEWftGlPmLI9Q3jc6mhHt7E/y/5/ama6acWVfAp8+ooiIvMzs1i4jI6FBYkRFnjKG708EYiA6O3u7EeQGb9oEkhSEfn15cxRJNoBURyUkKKzLiLMtiwRkRqmr9VE8c2V6O7liKkqH6KKURP39/QR11RUFKtOGgiEjO0mogGTHxmIsxXk+KZVnU1AVHrKfDGMNjW7r5y99u5aXdfenjc6vyFFRERHKc3uVlRCQTLs8/2U9hkY+FZ+XhD4zccExvLMU9L7fw8p5+AJ7Z0cvZkwpH7PFFRCS7KazISXNdw+svDjLQ5+I4hlTKjFhYWbt/kG89v4+uaAq/DR87rZKrZpaNyGOLiMjYoLAiJ239qihtLSl8PjhzSf6ILFF2jeHX6zr5+Zo2XAN1RUE+f14tTWUqlS8iMt4orMhJ2bk1zvbNCQAWnp03YmX017UO8rPVXhXapZOL+MyZ1YT9mmIlIjIeKazICWtvTfHW617F2Blzw9TUBUfssedV5fPBWWXUFQW5eEqxliSLiIxjCityQqKDLq+/MIAxUFsfYNrsky/89syOXuZX56WXJn9i0YSTfkwRERn71K8ux811DK+/MEAibigqtllwRt5J9Xw4ruE/Xt/P3c/v45vP7CXpjF4hORERGXvUsyLHbf3qKF0dDv4ALF6Sj99/4kGlN+7w/57by5qWQcCrm+JThBYRkUMorMhx2bszcXBC7Vn55Bf4TvixmvsSfPWp3ezrSxL2W3zunFrOqVf9FBERGU5hRY5ZX4/D6le9HpCps0InVUp/Y1uUf1y5h964Q2Wen68sraOxVMuSRUTkSAorcsySCUMgaFFa5GPm3BMPFq4xfO/lFnrjDlPKwty+tI5SlcwXEZF3oE8IOWZllX4uuMwbprHsE5+nYlsWf3/BRH61tp2/Uv0UERF5D/qUkPfkHrI6JxSyCYWO/9fGGMP2rlj6dm1RkM+dW6ugIiIi70mfFPKuBgdcVvyxl907Eif8GI5r+P6r+/n8/+zgtb39I9g6EREZDxRW5F1t3xQnFjXs3BLHuMdf/yTlGr79wj4e2dyNa6AzmhqFVoqISC7TnBV5V7MXhAmGLCbWB457nkrCcfnmM3t5bd8Afhs+d04t5zcWjVJLRUQkVymsyLuybItps49/5U/CcfmnlXt5o3mAoM/iyxdMZFFtwSi0UEREcp2GgeQIqZRh8/oYzgmWvU86Lv/8jBdUQj6Lf1hWp6AiIiInTGFFjrBxTZSNb8V45dmBE/p+n21RGvET9FncvqyOeVX5I9xCEREZTzQMJMO070+my+lPmXFiOynblsVnz6rm6pll1Jec/G7MIiIyvqlnRdJSScOqV7xy+g1TgkyoOfZy+o5reGRzF87QiiHbshRURERkRCisSNq6VVGig4a8fJvZCyLH/H3GGL7/agv/3yv7+fYL+0axhSIiMh4prAgAbS1Jdm3zhn8WnJmHP3Dsy5R/tqqNx7b0YFtwrnZNFhGREaawIqRShjWvRQGYPC1IxYRjn8r0m/UdPLS+E4DPnFnNufWqoyIiIiNLYUXYtC7G4IBLOM9i5rxjH/55clsP973ZBsCfn1bJpVNLRquJIiIyjimsjHM9XSm2vR0HYP7pxz78s6ZlgH9/uRmAD84q47o55aPWRhERGd+0dHkcc13D6lejGAO1kwJU1R776p94yuCzLM5uKODjCytHsZUiIjLeKayMY9s3x+npcggELOYsPPbhH4Az6gr4v5c1UFsUxLaOb88gERGR46FhoHHKGEPzniQAsxaECUfe+1ch4bi0DSTTtxtLwwR9+hUSEZHRpU+accqyLM5dVsBpZ+VR3xR8z/ONMdzzUguff2QHb7dHT0ELRUREPBoGGsds22JS43sHFYCHN3axckcvPgviKXeUWyYiIvL/t3fncVFX++PHX8MAAoLiAiKbwggoIGVIiukFM69aplm4pKVGblmZispVC7fSy08tRTPLbpnmNZduZj09AakAACAASURBVO6RqV31m5amuaO4oKIossk62+8PrpMjA4ICM8j7+XjwqPl8zud83jPHmXnPOedzPn+RnpVaRqvVcz6psEJ3VD56LZflh9MAiA51JcRNbkwohBCi+kiyUsskny7k2KF8ft19u1zlb+Sqmfvfq+j00NmnHs/5N6jiCIUQQghjkqzUMg6OVtSxU+Dte/+bDBZqdMzZc4XsQi2qhnV440k3FHLljxBCiGomc1ZqGQ9vW5o0tUFZjpZffzydc7cKcKqj5B+dPKljLbmtEEKI6ifJSi1U3lVqXwxsxLUcNc+0qI+rY/kXjBNCCCEqkyQrtYBOp+fAL7l4NLPFs5lNuYdy7G2siOnoXsXRCSGEEGWTfv1a4FJyETeuaTh+OB+1uuyrgAo1On46l4leX/6rhYQQQoiqJMnKI06t1nP6WAEA/kF22NqW3eRfHEoj4f+usfjXa9URnhBCCHFfkqw84s6eLKCoUE9dRyuaq8peAG7vxWy2JWUC0LFZveoITwghhLgvSVYeYfl5OpLPFALF9/+xUpY+VyXttpqP/9eb8mJgQ9o0lYXfhBBCWAZJVh5hJ4/mo9NCQxclbh6lX82j0+tZ+H+p5Kp1BDS2Y9BjLtUYpRBCCFE2uRroEZWZruHKxeI7JAc9bl/mFUCbTmdw7HoedZQKxnVwx9pKFn4ToqI0Gg15eXnmDqPGyM/Pp6ioyNxhiPt4kHZycHDA2rpy0wtJVh5Ber2e40eK74zs2cwG54alN/PtQi2rjtwE4LUnXGnqVL4bGwoh/qLRaMjNzcXJyQkrK+mwLg8bGxvUarW5wxD3UdF20ul05OTkULdu3UpNWORd9QhKS9Vw64YWKyW0DLEvs6xjHSWzunjRw8+Z7n7O1RShEI+WvLw8SVSEAKysrHBycqr0XkbpWXnE6PV6Tv1Z3Kvi06IO9g73//D0b2yPf+OykxohRNkkURGiWFW8F+Td9Yi5mqImO1OHtQ20aFX6zQovZBRwKbOwGiMTQgghHowkK4+Yc6eKExBVgB22dUw3r0anZ8H+VMZvvcDBy7erMzwhhBCiwiRZecS0j6iLX2AdfP1L71XZcPIW5zMKsbNW4NfYrhqjE0II09q1a8eyZcvMdv6xY8cSHR1dqXWuWbOGVq1aVWqdtZXMWXnE2NaxomXr0uefXM0u4pujxVf/RIc2wdlO/gkIURuNHTuWdevWAaBUKmnSpAldunThH//4B87Oj85k+5SUFNq3b294bGNjg4eHB3379uWdd94xLOswc+ZMuSeaBZNvqkdEfp4OO3tFmeup6PR6Pv41FbVOz+NuDnT2kSX1hajNOnfuzIcffohGoyEpKYnx48eTnZ3NkiVLzB1ahWm1WhQKRamTO7/55hsCAgIoKiriwIEDTJw4kSZNmvDyyy8DUK+efB5aMhkGegSo1Xp2b89h/8+3KSzQlVou8VwWx9LyqaNUMLqdW5mJjRDi0Wdra4urqyvu7u5ERETQq1cvdu/ebVQmOzubSZMmERISQkBAAH379uX48eNGZRYsWEBISAj+/v5MmDCB2bNn07VrV8P+qKgo4uLijI4ZPHgwY8eOLTW2Tz/9lC5dutCiRQvatm3L5MmTyc3NNey/M8Ty448/EhkZiY+PD5cvXy61vgYNGuDq6oqnpycvvvgioaGh/Pnnn4b99w4DRUVFMXXqVKZOnUqrVq0ICgoiPj7eqPclMzOTMWPGEBgYiEql4pVXXiE5ObnUGAB27NhB9+7d8fX1JTw83JAsirJJsvIIyLipQavRU1igx9bWdAKSka9h+aE0AAY95kITR1n8TYiqotfr0RcWVP/fQwxjXLx4kV27dmFj89etOfR6PYMHDyYtLY2VK1eydetWWrduTf/+/cnIyADgP//5D4sWLWLq1Kls3boVDw8PVqxY8dCvoZWVFTNnzmTnzp0sWLCAvXv38v777xuVyc/PZ/HixcydO5edO3fSuHHjctV95MgRjh07Rps2bcost27dOpRKJT/88AOzZs1i2bJl/Pvf/zbsHzduHEePHuXLL79k48aN6PV6Xn311VIXUdu1axdjxowhOjqan3/+mfj4eNauXUtCQkK54q7NZBjoEeDa1Iann6tHYb4ORSlL5TvVUfJiUCN+v3KbngENqjlCIWqZokJ0b/Wr9tNaLV4Ldco/aT4xMRE/Pz90Oh0FBQUATJs2zbB/7969nDp1iiNHjlCnTvGk/bi4OLZv387mzZt55ZVX+OKLLxgwYAD9+/cHir/Ad+/ebdQL8iCGDx9u+H9vb28mTpzI5MmTmTNnjmG7Wq1m9uzZBAUF3be+3r17Y2VlhVqtRq1WM2jQIPr27VvmMe7u7syYMQOFQkGLFi04deoUy5YtY9CgQSQnJ7Njxw42bNhAWFgYAIsWLSIsLIxt27bx/PPPl6gvISGBN998k379iv9tNGvWjIkTJ/LBBx8wfvz4cr0utZUkK48IewerMheAs7ZSEBXUiBcDG2Ilwz9CCKBDhw7MmTOH/Px8Vq9eTXJystFQyJ9//klubi7BwcFGxxUUFHDx4kUAkpOTGTJkiNH+xx9/nL179z5UbHv37mXRokUkJSWRk5ODVquloKCAvLw8HBwcgOJhrMDAwHLV98knn+Dn54dGo+HkyZPExcXh7OzMlClTSj3miSeeMBouDw0N5dNPP0Wr1XL27Fmsra154oknDPsbNmyISqXi7NmzJus7evQoR44cMepJuZMo5ufnY28vi3OWRpKVGkxdpCf3trbMe/9odMXdwnduTiiJihDVwLZOcS+HGc5bEQ4ODvj4+AAwa9YsoqKi+PDDD5k0aRJQ/EXq6urK+vXrSxxbv359w//fO//t3uEoU/PjypqncfnyZQYPHswrr7zCxIkTcXZ25uDBg8TExBgNsdjZ2ZV77p27u7vhufr5+XHp0iXmzp3L+PHjsbOr+BIOpQ25lTUUp9friYmJoUePHiX23em5EqZJslKDnU8q5PSxAnz86xDcxnRGvvHkLXadz2bUk00IdHWo5giFqJ0UCkWFhmMsxfjx43n11VcZPHgwbm5utG7dmhs3bmBtbY2Xl5fJY3x9ffnjjz+IiooybDt69KhRmUaNGnH9+nXDY61Wy6lTpwgPDzdZ55EjR9BoNEybNs1wdc8PP/zwsE/PiFKpRKPRoFarS01WDh06VOKxj48PSqXS0Etz6NAhwzDQrVu3SE5Oxs/Pz2R9wcHBnDt3zpA0ifKTCbY1lLpIT/Lp4tVqGzRSmixzM0/NmmM3uZhVyLXbcndTIUTZOnTogL+/P4sWLQKgU6dOhIaGEh0dza5du0hJSeHgwYPEx8dz5MgRAKKjo1m9ejVr164lOTmZBQsWcPLkSaMej6eeeoqffvqJxMREzp49y+TJk8nKyio1jmbNmqHRaPjiiy+4ePEi69evZ+XKlQ/13DIyMkhLS+Pq1avs3LmTzz//nA4dOuDk5FTqMVevXmX69OmcPXuWDRs28MUXX/D6668DxUlat27dmDRpEgcOHOD48eOMGTMGNzc3unXrZrK+cePGsX79eubPn8/p06dJSkri+++/Jz4+/qGeW20gPSs11PmkQtRqPY71rHD3tDFZ5ovf0yjQ6GnZ2J5IWVNFCFEOI0aMYPz48YwePRoPDw9WrlxJfHw8MTExpKen4+LiQvv27Q1X3rz44otcvHiRWbNmUVhYyPPPP0+/fv04fPiwoc4BAwZw4sQJ3nnnHaytrRk+fDhPPfVUqTEEBwczbdo0lixZwpw5c2jfvj2TJ0/mnXfeeeDnNWDAAKC4R8XV1ZUuXboQGxtb5jFRUVEUFBTQs2dPlEol0dHRvPLKK4b9H374IXFxcQwZMoSioiLat2/PypUrja6oultkZCRfffUVH330EUuWLMHGxoYWLVoY1noRpVPoZcm+B3Ljxo1SL0+rKIVCQdOmTUlNTS3XpYfqIj0/bcpGrdbzRLgDHt4lL0M+ci2XuJ9SsFLA/O7N8W1Y87qkLU1F20lUP3O1UXZ2tiwqdpcBAwbg4uJi6KExxcbGptI+Q6tCVFQUgYGBzJw509yhmNWDtlNp7wkbGxtcXFwqXJ/0rNRA58+W3aui0en57GDx+HAPP2dJVIQQVSY/P58VK1YQGRmJUqlkw4YN/PLLL6xevdrcoYlHiCQrNYxGo+f8meK5Kn6BdibXVdlxNpPL2UXUq6Nk4GMVz2CFEKIidu7cSUJCAoWFhahUKpYtW8bf/vY3c4clHiGSrNQwl5KLKCrU4+BohbuX6XHRI9eKF2MaGNIYR1vTk2+FEKIy2Nvbs2bNGnOHUelMXa4tzEeSlRpEq9Vz7lTxKpMtWtbBqpTVav/RyYPfruTyhHvd6gxPCCGEqBKSrNQgly8UUZCvx85egWfz0u/to1AoCPN0rMbIhBBCiKoj66zUEDqdnrOniueq+AbUQaks2avyc3IWuUXa6g5NCCGEqFKSrNQQqZfV5N3WYWOroJlvyWWZT6TlsWB/KqN/SCZPLQmLEEKIR4ckKzVExs3i+2j4+tfB2sa4V0Wn1/PFoTQA2nk64WAjk2qFEEI8OmTOSg0R/IQDns1scXAqmV/uvZhDUnoBdtZWDAxpbIbohBBCiKojPSs1iHMja2xtjZtMo9Oz6ugNAPoENsTZXvJPIUT1mT9/Pl27djV3GCatWbOGVq1aVcu52rVrx7Jlyx66ntWrV1vs8vsnT54kNDSUvLy8aj+3JCsWLve2lsICXan7fzybSWqOmvp1lPRu2bAaIxNC1HRpaWm8++67hIeH4+PjQ9u2bRkyZAi//PKLuUOrVnv37uXVV18lKCgIlUpFZGQkM2bMIDU1tVrjKCwsZN68eYwdO9aw7fTp0wwfPpx27drh4eFhMiG6ffs2cXFxPPnkk6hUKnr16sUff/xR6nkmTZpksq6FCxfSq1cvVCqVySSvVatWtGnThs8+++whnuWDkWTFwp34o4DEH7JJuVBUYl+hRseaP28C0K91I+xtpDmFEOWTkpJCjx492Lt3L1OnTiUxMZFVq1bRoUMHpk6dau7wqs3KlSsN9zJatmwZu3bt4p///Cc5OTl8+umn1RrLli1bcHBwoF27doZt+fn5eHt7M2XKFFxdXU0eN2HCBH755RcSEhJITEwkIiKCAQMGmEy2tm3bxuHDh3FzcyuxT61W07NnTwYPHlxqjP369WPlypVotdV7IYd8u1kwnVZPYYEOnQ6cG5acNFuo1dPG3RE3Rxu6tWhghgiFEDXVlClTANi8eTM9e/ZEpVIREBDAyJEj+eGHHwzlrly5wmuvvYafn59h/40bN8qse82aNURERODr68vf/vY3li9fXmb5n3/+mRdeeIFWrVoRFBTE4MGDuXDhgmF/SkoKHh4ebNmyhaioKFQqFc888wy//fZbifOGhYWhUql4/fXXycjIKPO8V69eJS4ujujoaD788EM6dOiAl5cX7du3Z968eYwbN85QdvPmzXTu3BkfHx/atWvH0qVLy6w7OzubSZMmERISQkBAAH379uX48eNlHvP999/z97//3Wjb448/znvvvUfv3r2xtS25vlZ+fj5btmxh6tSptG/fHh8fH2JiYvDy8mLFihVGZVNTU5k6dSqLFy/G2rrklIEJEyYwYsQIWrZsWWqMkZGRZGRksH///jKfS2WTZMWCWSkVPNXFkcjuTjjVK5ms1Kuj5J3wpiQ854ONiXVXhBDmVaDRlfpXpNWVu2yh5v5lKyIjI4Off/6ZoUOH4uDgUGJ//fr1AdDr9URHR5OZmcm3337L6tWruXjxIm+88Uapda9atYr4+HhiY2PZtWsX//jHP5g7dy5r164t9Zi8vDxGjBjB5s2bWbNmDVZWVgwbNgydzvh5xcfHM2rUKHbs2IGvry9vvvkmGk3xlZKHDh0iJiaGIUOGsGPHDjp06EBCQkKZr8OmTZsoKipi9OjRJvffeR2OHj3KqFGj6NWrF4mJiYwfP565c+eWepsBvV7P4MGDSUtLY+XKlWzdupXWrVvTv3//MhOoAwcOEBISUmbM99JqtWi1WurUMV7Sws7OjoMHDxoe63Q6xowZwxtvvEFAQECFznE3W1tbAgMDOXDgwAPX8SBkNqaFUygUONUv+1LkOtaScwphifqvOVPqvlD3usR19jI8Hrw+iUKt3mTZYFd7PujazPB4+IZzZBcad8N/P6j0X8P3unDhAnq9nhYtWpRZ7pdffuHkyZPs378fDw8PABISEujcuTN//PEHjz/+eIljFixYQFxcHM8++ywA3t7enDlzhq+//pp+/fqZPM9zzz1n9Hj+/PmEhIRw5swZo1/5o0aN4plnngGKewE6d+7MhQsXaNGiBf/617+IjIzkrbfeAkClUvHbb7+xa9euUp/f+fPncXJyokmTJmW+Dp999hkdO3Y09LSoVCqSkpJYunQp/fv3L1F+7969nDp1iiNHjhiSiLi4OLZv387mzZt55ZVXShyTlZVFVlaWyeGZsjg6OhIaGsrChQvx8/PDxcWFDRs2cPjwYXx8fAzlPv74Y6ytrXn99dcrVL8pbm5upKSkPHQ9FVErk5VNmzaxc+dO9Ho9rVu35rXXXkOhsKyeiZtpGuo5W5W4+gcgPU/N8kM36B/SCM96JReIE0KIsuj1xUnR/T73kpKScHd3NyQqAP7+/tSvX5+kpKQSyUp6ejpXr14lJiaGiRMnGrZrtVqcnJxKPc+FCxeYO3cuhw4d4tatW4YelStXrhglK3dP+rwzf+PmzZu0aNGCpKQkevToYVRvaGhomcmKXq8v12d/UlIS3bp1M9oWFhbG559/jlarRak0/kH5559/kpubS3BwsNH2goICLl68aPIcBQXF9327t4ekPBISEoiJiSE0NBSlUknr1q3p06cPf/75J1DcM/Svf/2Lbdu2Vcp3nZ2dHfn5+Q9dT0XUumQlOzub7du3M3/+fKytrZk2bRpJSUn4+/ubOzQDjVrPb//NRafX07GLE/Wcjd8I646ls+diNrfy1Ua/toQQlmVN/9I/V+69D+mKKL9Sy9779bLsBdVDRAU+Pj4oFAqSkpLo3r17qeVK+zIvbfudJGPu3Lm0adPGaN+9X+h3Gzp0KO7u7vy///f/cHNzQ6fT8fTTT6NWq43K3T3P4s7575zzTgJWEb6+vmRnZ3P9+vUye1dMPd+yzqfT6XB1dTV55+Y7Q0v3atCgAQqFgqysrHJG/5fmzZvz7bffkpeXR05ODk2aNGHUqFF4eRX33P3666/cvHmTJ5980nCMVqtl5syZfP755/z6668VOl9mZibNmlXvd0+tHD/QarWo1Wo0Gg0ajYZ69eqZOyQjl5ILUav12NlZ4VTPuIlu5Kr58VwmAC+HuJgjPCFEOdlZW5X6Z6u0KnfZe4d6TZWpiAYNGhAZGcny5ctNrplx5wvT39+fK1eucOXKFcO+M2fOkJ2djZ9fyeTKxcUFNzc3Ll68iI+Pj9Gft7e3yVhu3bpFUlIS77zzDp06dcLPz++BvrD9/f05dOiQ0bZ7H9/rueeew9bWliVLlpjcf/frcO8cjd9++w1fX1+TSVjr1q25ceMG1tbWJV6Hhg1NLzFha2uLv78/Z86UPnR4Pw4ODjRp0oTMzEx2795t6A166aWXSExMZMeOHYY/Nzc33njjDVatWlXh85w6dapEr1FVs7ielRMnTrBx40bOnz9PRkYGEyZMMMoGAbZv387GjRvJzMzE09OToUOHlnvhn3r16vH8888zevRorKys6Nq1a4XHCKuSTqcn+cxfNyxU3PPza92xdDQ6CGniQHCTkhPjhBCiPGbPnk3v3r157rnnmDBhAq1atUKr1bJnzx5WrFjB7t276dSpE61ateLtt99mxowZaDQapkyZQnh4OI899pjJemNiYnjvvfdwcnKic+fOFBUVcfToUTIzMxk5cmSJ8s7OzjRo0ICvv/4aV1dXrly5wpw5cyr8fKKjo+nduzdLliyhW7du7Nmzp8whIAAPDw+mTZvGu+++y+3bt4mKisLLy4vU1FTWrVtH3bp1mTZtGiNHjuTZZ5/lo48+olevXvz+++98+eWXzJ4922S9nTp1IjQ0lOjoaKZOnYpKpeLatWvs3LmT7t27l/raRUREcODAAYYPH27YVlRUZEhg1Go1165d49ixY9StW9cwJ2XXrl3o9XpUKhUXLlxg1qxZqFQqw3yahg0blkiSrK2tcXFxMZq3dOXKFTIyMrh69SparZZjx45hbW2Nl5cXdevWBYqvzLp27RqdOnUq87WtbBaXrBQWFtK8eXM6d+7M/PnzS+zft28fy5cvZ9iwYQQEBJCYmMjs2bP56KOPaNy4eKn52NhYwwzxu02dOhVbW1sOHTrExx9/jK2tLbNnz+bEiRMEBgZW+XMrj2tX1OTn6bGto8CrufFlatdvF5Fo6FWRZfWFEA/O29ubbdu2kZCQwMyZM0lLS6Nhw4aEhIQYkgWFQsEXX3zBu+++y4svvoiVlRWRkZG8//77pdY7cOBA7O3t+eSTT/jggw9wcHCgZcuWDBs2zGR5KysrlixZQlxcHF26dMHX15dZs2YRFRVVoecTGhrKvHnzmDdvHvPnz6dTp06MGTOGhQsXlnnc0KFD8fX15dNPP2XYsGEUFBTg6enJM888w4gRI4DinpKlS5cyb948Fi5ciKurKxMnTjQ5uRaKX7eVK1cSHx9PTEwM6enpuLi40L59e8P3lCmDBg2iW7duZGdnG3r8r1+/bjRfZunSpSxdupTw8HDDMFN2djb//Oc/SU1NxdnZmWeffZbY2FhsbGwq9BrOnTuXdevWGR7fOe+6devo0KEDABs2bCAiIgJPT88K1f2wFPoHGeirJv369SvRszJlyhR8fHyMMs9x48YRFhbGwIED71vn/v37OX78uOGNs3HjRvR6Pb179zZZXq1WG42bKhQK7O3tuXHjhsmE6EEoFArc3Ny4du0ae3ZkkZGuxT/Ijpat7Y3KLdp/lR/PZfF407rM7GK6S1VUnbvbyYLfNrWaudooKyvL4oaTLZ2NjU2JOSkCRo4cSXBwMG+//ba5QwGM26mwsJCOHTuyZMkSwsLCyjwuOzvb5PycOz06FWVxPStl0Wg0JCcn88ILLxhtDwkJ4fTp0+Wqo1GjRpw5c4aioiKsra05fvy44VI4U7777jujSVI+Pj7Ex8c/0It9X7p6ZKTfwkqpoN1T3jjU/at5LmfksTP5FABvd25J06amJ2mJqmdJw4bCtOpuo/z8/Ar/ihXIa2bCjBkz2L59u0W9NndiuXTpEuPGjTP0spTF1taWpk2bVloMNSpZyc7ORqfTlcjW6tevT2ZmZrnq8Pf3p02bNsTGxqJQKAgODqZt27allu/Tpw89e/Y0PL4zI7wqelYO7rsKgIe3DVnZN8jK/qtMkVrHoMcaczGzEBerPFJTq/9GUrWd9KxYPnO1UVFRkfQSVJD0rJjm5ubGkCFDLOa1ubudvL29GThwYLliKyoqMrncf63oWbnD1CVzFbl2/OWXXy73XS1tbGxKzXAr88MwJ6uIq5eL7//j61+nRN121gpeCmpU6ecVFafX66UNLJy0kRDmV5nvwRp16XK9evWwsrIq0YuSlZVV6rXrNcWff9wCPTRuYl1iXRUhhBCiNqtRyYq1tTW+vr4cPXrUaPvRo0cf6l4H5qZW6zl1rDgB8w0wXr3wWk4RE7dd4P9ScswRmhBCCGF2FjcMVFBQwLVr1wyP09LSuHDhAo6OjjRu3JiePXuyaNEifH198ff3JzExkZs3b9K1a1czRv1wLiUXoi7S4VjPClc34yb59kQ6Z9IL2JaUSXuv0perFkIIIR5VFpesnDt3jhkzZhge37nFdUREBG+++SYdOnQgJyeHb7/9loyMDLy8vJg8eXLVXJ1TTVJTiueqqALsjObe3MhVszO5eAXF/sGNzBKbEEIIYW4Wl6wEBQWVeStxKF6o5t6bStVkHZ52Ij+nLnZ1c422f3fyFhodBDdxoJWrrFYrhBCidqpRc1YeVVZWClq0rI/S+q9elYx8DT+eLZ7H0k96VYQQQtRikqxYqA0nb1Gk1RPQ2J4QuQeQEMICRUVFERcXZ+4wqtyDPs/Vq1eXe5mM6nby5ElCQ0NN3sjSEkmyYoGyC7VsS8oAintVKrKGjBBClMfYsWPx8PDAw8ODZs2aER4ezsyZMx/qy6tdu3YsW7asEqN8cPv27cPDw+OB7uBcGQoLC5k3bx5jx441bDt9+jTDhw+nXbt2eHh4mHytbt++TVxcHE8++SQqlYpevXrxxx9/GJW50273/n3yySeGMufOneO1114jODiYgIAAevfuzd69ew37W7VqRZs2bfjss8+q4NlXPklWLJCjrRXjOrjztG89Qt3rmjscIcQjqnPnzhw+fJh9+/YxadIkvvrqK2bOnGnusB4JW7ZswcHBgXbt2hm25efn4+3tzZQpU3B1dTV53IQJE/jll19ISEggMTGRiIgIBgwYYLQa7OHDh43+PvzwQxQKBc8++6yhzODBg9FoNKxdu5atW7cSFBTEkCFDSEtLM5Tp168fK1euRKvVVsErULkkWbFAVgoF7b2ceCfcXXpVhBBVxtbWFldXVzw8POjTpw8vvvgi27dvN+w/c+YMr776Kn5+fjz22GO8/fbb3Lp1y2RdUVFRXL58menTpxt+6d+xefNmOnfujI+PD+3atWPp0qVGx7Zr146EhATGjx+Pv78/YWFhfP3112XGrtfrWbJkCeHh4ahUKp555hk2bdoEQEpKCn379gUgMDAQDw8Pox6Ou926dYvRo0cTGhqKSqWiS5cubNiw4f4v3n18//33/P3vfzfa9vjjj/Pee+/Ru3dvbG1tSxyTn5/Pli1bmDp1Ku3bt8fHx4eYmBi8vLwMV8YCuLq6Gv1t376dDh060KxZM8NzunDhAm+99RaBgYH4+voyZcoU8vPzje6jFxkZHrMsYwAAE/ZJREFUSUZGBvv373/o51vVJFmxMLJEuBCPDo1GX+E/ne6vzwCdrnibVqO/b72Vwc7OznDfl+vXr/PSSy8RGBjI1q1bWbVqFTdv3mTkyJEmj122bBlNmzZlwoQJhl/8ULxo56hRo+jVqxeJiYmMHz+euXPn8s033xgd/+mnnxISEsL27dsZMmQIkydP5uzZs6XGGh8fz5o1a5gzZw47d+5k+PDhjBkzhv379+Pu7m4YYtmzZw+HDx8utceosLCQkJAQvvrqK3bu3MmgQYMYM2YMhw4dqvDrd7cDBw4QEhJSoWO0Wi1arZY6dYwXB7Wzs+PgwYMmj7lx4wY//fST0dyYBg0a4Ofnx/r168nLy0Oj0fD111/j4uJiFJOtrS2BgYEcOHCgQnGag8VdulybFWl1TNx2gY7N6vF8QAPqWEsuKURNtvXbis+XCO3ggLtX8a/ua1fU/L4vj0YuSjo8/deikD9tyqao0DhBeb6/80PFevjwYb777js6duwIFK9x1bp1ayZPnmwoM3/+fMLCwjh37hwqlcro+AYNGqBUKnF0dDQa4vjss8/o2LEj48aNA0ClUpGUlMTHH3/MSy+9ZCj39NNPM3ToUADefPNNli1bxr59+2jRokWJWPPy8li2bBlr1qwx3Ii2WbNmHDx4kK+//prw8HCcnYtfj8aNG5d5O5amTZsyatQow+Po6Gh+/vlnNm3axBNPPFGu1+5eWVlZZGVlVfju346OjoSGhrJw4UL8/PxwcXFhw4YNHD58GB8fH5PHrFu3DkdHR3r06GHYplAoWL16NdHR0fj7+2NlZYWLiwtff/11idfCzc2NlJSUij/JaibJigXZeS6L8xmF5BRm0LtlQ3OHI4R4xCUmJuLn54dWq0WtVtOtWzfef/99oLhHZN++ffj5+ZU47uLFiyWSldIkJSWVWBcrLCyMzz//HK1Wi1JZfC+0wMBAw36FQoGLiwvp6ekm6zxz5gwFBQUlrrRRq9UEBweXK647tFotixcv5ocffiA1NZWioiKKiopwcHjwqzALCgoASvSQlEdCQgIxMTGEhoaiVCpp3bo1ffr04c8//zRZ/ptvvqFPnz7Y2dkZtun1eqZMmULjxo357rvvsLOz49///jdDhgxhy5YtNGnSxFDWzs6O/Pz8CsdZ3SRZsRAanY7/nCh+Y77QqiE2SpmrIkRN1+Olit9g1equDlU3Dxt6vFSfez8NuvSs93CB/U+HDh2YM2cONjY2NGnSxOgO83q9nq5duzJlypQSx939ZXc/er2+xNw7U8Pd1tbGX0cKhQKdTmeyzjvbV6xYUaL3wtRckLJ8+umnLFu2jBkzZtCyZUscHByYNm2aYTjsQTRo0ACFQvFAVyI1b96cb7/9lry8PHJycmjSpAmjRo3Cy8urRNlff/2Vc+fOGV0FBPDf//6XxMRETpw4gZNTcY/cnDlz2LNnD+vWreOtt94ylM3MzDTMdbFkkqxYiJ1nbnDtthqnOkq6tni47lwhhGWwtn64Hx1WVgqj5KWy6r3DwcGh1OGF4OBgtmzZgpeXV4lEojQ2NjYlrizx9/cvMSfit99+Q6VSGXpVKsrf3586depw5coVwsPDS40FuO+VLr/++ivdunUzDEnpdDrOnz9vskepvGxtbfH39+fMmTNEREQ8UB0ODg44ODiQmZnJ7t27mTp1aokyq1evJiQkhKCgIKPtd3pKrO75x2NlZVUiATx16hTPPffcA8VYnWRShAXQ6/Us/7+LAPQMaICdzFURQpjZ0KFDyczMZPTo0Rw+fJiLFy+ye/duxo8fX2oC4OXlxa+//kpqaqrhqqGRI0fy3//+l48++ohz586xdu1avvzyS0aPHv3AsTk6OjJy5EimT5/O2rVruXDhAseOHWP58uWG27V4enqiUChITEwkPT2d3Nxck3U1b96cPXv2cPDgQZKSkoiNjeXGjRtlnn/OnDmMGTOmzDIRERElkrSioiKOHTvGsWPHUKvVXLt2jWPHjnH+/HlDmV27dvHzzz9z6dIl9uzZQ9++fVGpVPTv39+orpycHDZt2mRy0bm2bdtSv359xo4dy/Hjxzl37hyzZs0iJSWFLl26GMqlpKRw7do1OnXqVOZzsQTSs2IBDqXmknTjNnbWCp7zb2DucIQQAjc3NzZs2MDs2bMZNGgQhYWFeHp6EhkZWeIX+x0TJkwgNjaWp556isLCQq5cuULr1q1ZunQp8+bNY+HChbi6ujJx4kQGDBjwUEMtkyZNonHjxixevJhLly5Rr149Wrduzdtvvw0UT5yNiYlhzpw5jB8/nqioKBYsWFCinrFjx5KSksKgQYOwt7dn0KBBdOvWjZycnFLPff36da5evVpmfHfqyc7Opl69eobj7p6/s3TpUpYuXUp4eDjr168HIDs7m3/+85+kpqbi7OzMs88+S2xsrNEQHRRfGq3X63nhhRdKnLthw4asWrWK+Ph4+vXrh0ajwd/fny+++MKoF2bDhg1ERETg6elZ5nOxBAq9XCv7QG7cuPFQb7S7Tf3xEsfS8ujdqiHRT5heKEiYn0KhoGnTpqSmpsol5hbKXG109xeSKB8bG5tK+wy1VCNHjiQ4ONiQQFmSwsJCOnbsyJIlSwgLCyu13IO2U2nvCRsbG1xcXCpcn4w3WIARYU14NtBNrgASQohHyHvvvUfdupa5Cvnly5cZM2ZMmYmKJZGelQdUmT0r8ou9ZpB2snzSs1Jz1IaelUeB9KwIIYQQQpSDJCtCCCGEsGiSrAghhBDCokmyIoQQlaC01VaFqG2q4r0gyYoQQjwkBwcHcnJyJGERtZ5OpyMnJ+eh7q1kiiwKJ4QQD8na2pq6dety+/Ztc4dSY9ja2lJUVGTuMMR9PEg71a1bt9y3aCgvSVbKadu2bWzfvh1PT09iYmLMHY4QwsJYW1vL5cvlJMsA1AyW1E6SrJRT9+7d6d69u7nDEEIIIWodmbMihBBCCIsmyYoQQgghLJokK0IIIYSwaDJn5QFV9kznqqpTVD5pJ8snbVQzSDvVDJXZTg9al9zIUAghhBAWTYaBKsn8+fPLvf3ebfn5+cTGxpKfn18lsZVHafFXVz3lPe5+5craX5E2MrXd3O1UWW30oHVV5BhztZO52+jeeMxRV2W9l+5XRj7zHq4e+cyrGElWKsnly5fLvf3ebXq9nvPnz5v1OvbS4q+uesp73P3KlbW/Im1karu526my2uhB66rIMeZqJ3O30b3xmKOuynov3a+MfOY9XD3ymVcxyunTp083dxCPihYtWpR7+93b1Go1GzZs4IUXXsDGxqbK4ruf0uKvrnrKe9z9ypW1vyJtdO92S2inymqjB62rIseYo50soY3ujsdcdVXWe+l+ZeQz7+Hqkc+88pM5KxYgLy+PoUOHsnz58kq/n4KoPNJOlk/aqGaQdqoZLKmdpGfFQlhZWREUFIRSqTR3KKIM0k6WT9qoZpB2qhkspZ2kZ0UIIYQQFk0m2AohhBDCokmyIoQQQgiLJsmKEEIIISyaJCtCCCGEsGiSrAghhBDCosldpGqQmzdvsnjxYrKyslAqlbz00kuEh4ebOyxhwty5czlx4gTBwcHExMSYOxzxP7///jsrVqxAr9fTu3dvunTpYu6QhAny/rFs5vgukmSlBlEqlQwdOpTmzZuTlZVFbGwsbdq0wc7OztyhiXv06NGDzp07s3v3bnOHIv5Hq9WyYsUKpk2bhr29PbGxsbRr1w5HR0dzhybuIe8fy2aO7yIZBqpBGjRoQPPmzQGoX78+jo6O3L5927xBCZOCg4Oxt7c3dxjiLmfPnsXT05OGDRtib29PmzZt+OOPP8wdljBB3j+WzRzfRdKzUolOnDjBxo0bOX/+PBkZGUyYMIEnn3zSqMz27dvZuHEjmZmZeHp6MnToUFq1alXhc507dw69Xk/jxo0rK/xaozrbSVSeh223jIwMGjZsaCjbqFEjbt26Va3PoTaQ95flq8w2qq7vIklWKlFhYSHNmzenc+fOJm/BvW/fPpYvX86wYcMICAggMTGR2bNn89FHHxkaOjY2Fo1GU+LYqVOnGj5oc3JyWLx4MaNGjaraJ/SIqq52EpXrYdvN1GLdCoWiOkKvVSrj/SWqVmW1UXV+F0myUonatGlDmzZtSt2/adMmnn76acOkvqFDh3LkyBF27NjBwIEDAYiPjy/zHGq1mnnz5tGnTx8CAgIqL/hapDraSVS+h223hg0bGvWkpKen4+fnV+Vx1zaV8f4SVasy2qi6v4tkzko10Wg0JCcn89hjjxltDwkJ4fTp0+WqQ6/X8/HHHxMUFMTf/va3qgiz1quMdhLVrzzt1qJFC1JSUrh16xb5+fkcPny4RHlRteT9ZfnK00bm+C6SnpVqkp2djU6no379+kbb69evT2ZmZrnqOH36NPv378fb25uDBw8C8Pbbb+Pt7V3p8dZWldFOAB988AHJyckUFhYyatQoJkyYQIsWLSo7XPE/5Wk3pVLJ4MGDmTFjBjqdjt69e+Pk5GSOcGut8r6/5P1jPuVpI3N8F0myUs1MjZGXd9y8ZcuWrFmzprJDEiY8TDtB8dwVUf3u125t27albdu21RmSMOF+7STvH/Mrq43M8V0kw0DVpF69elhZWZX4dZ6VlVUigxXmI+1UM0m71QzSTpbPUttIkpVqYm1tja+vL0ePHjXafvToUZkoa0GknWomabeaQdrJ8llqG8kwUCUqKCjg2rVrhsdpaWlcuHABR0dHGjduTM+ePVm0aBG+vr74+/uTmJjIzZs36dq1qxmjrn2knWomabeaQdrJ8tXENlLoTS0+IB7I8ePHmTFjRontERERvPnmm8BfC+1kZGTg5eXFkCFDCAwMrO5QazVpp5pJ2q1mkHayfDWxjSRZEUIIIYRFkzkrQgghhLBokqwIIYQQwqJJsiKEEEIIiybJihBCCCEsmiQrQgghhLBokqwIIYQQwqJJsiKEEEIIiybJihBCCCEsmiQrQgghhLBokqwIISrd6dOnWbt2Lbm5uSX2TZ8+nenTp1d/UPexe/duXn/9dfLz88t9zNWrV3n55ZdJTk6uwsiEEJKsCCEq3enTp1m/fr3JZGXYsGEMGzbMDFGVrrCwkNWrV9O7d2/s7e3LfZy7uzudOnXiq6++qsLohBCSrAghqpWnpyeenp7mDsPI7t27ycnJoUuXLhU+tnv37pw8eZLTp09XQWRCCABrcwcghHi0rF27lvXr1wPw1ltvGbZPmzaNoKAgwxDQnf+mpaXx1ltv8corr6DT6fjxxx/JzMzEz8+PESNG4Orqypo1a9izZw8FBQUEBwczYsQI6tevb3Teffv2sXnzZi5dugRAy5YtGThwID4+PveNeceOHYSGhlK3bl2j7fv37+eHH37gypUraLVanJ2dCQwMZPTo0YYyvr6+eHh48OOPPxIQEFDh10sIcX+SrAghKlWXLl24ffs227ZtY8KECTg7OwPctzdl+/bteHt78/rrr5Obm8uKFSuIj4/Hz88PpVLJG2+8wY0bN1i5ciVLly4lNjbWcOx//vMf1qxZQ2RkJC+99BIajYaNGzcSFxfHnDlzyjx3eno6ly5d4u9//7vR9jNnzrBgwQLCw8Pp27cvNjY23Lx5k2PHjpWoIygoiP3796PX61EoFBV5uYQQ5SDJihCiUjVq1IjGjRsD0Lx5c1xdXct1nIODAxMnTsTKqnh0Oicnh+XLl+Ph4cGkSZMM5a5cucKWLVvIy8vDwcGBmzdvsm7dOrp160Z0dLShXEhICGPGjGHdunWMGzeu1PPeGb65twfm9OnT6PV6RowYgYODg2F7ZGRkiTp8fHzYsWMHV69excPDo1zPVwhRfjJnRQhhEdq0aWNIVADDl/4TTzxhVO5OL8nNmzcBOHLkCFqtloiICLRareHPxsaGwMBATpw4UeZ5MzIyAEoMK6lUKgA++ugj9u3bx61bt0qt486xZZURQjw46VkRQlgER0dHo8fW1tZlbler1QBkZWUBMHnyZJP13m9YpqioCAAbGxuj7YGBgUycOJGtW7fy8ccfo1ar8fLyok+fPnTs2NGo7J1j79QlhKhckqwIIWo0JycnAMaPH4+Li8sDH3/79m0aNGhgtC8sLIywsDDUajVJSUl89913JCQk4Orqir+/v6Hc7du3jeoSQlQuSVaEEJWuOnsaHnvsMZRKJdevX6d9+/YVPv7OcNP169fx8vIyWebOkJKDgwNHjhzh/PnzRsnK9evXUSgUuLu7P9iTEEKUSZIVIUSl8/b2BmDLli1ERkaiVCpxd3ev0IJr5eXq6kq/fv345ptvuH79Oo8//jiOjo5kZmZy9uxZ7Ozs6NevX6nH+/n5YWtry5kzZ2jbtq1h+5o1a0hPT6d169Y0atSI3NxctmzZglKpJDAw0KiOpKQkmjdvXmLISghROSRZEUJUuqCgIF544QV2797NTz/9hF6vN6yzUhX69OmDp6cnW7ZsYe/evWg0GpydnVGpVHTt2rXMY62trWnfvj2//fYbAwcONGz38/MjOTmZVatWkZ2dTd26dfH19WXatGlGPTAFBQUcO3aM/v37V8lzE0KAQq/X680dhBBCmNO5c+eYPHkyH3zwAX5+fhU6dufOnXz55Zd88skn0rMiRBWRS5eFELWeSqUiPDycb7/9tkLHabVavv/+e/r06SOJihBVSJIVIYQABg8ejEqlqtBdl9PT0+nUqRM9e/aswsiEEDIMJIQQQgiLJj0rQgghhLBokqwIIYQQwqJJsiKEEEIIiybJihBCCCEsmiQrQgghhLBokqwIIYQQwqJJsiKEEEIIiybJihBCCCEs2v8HjuQv/mfwa9sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Times\n",
    "times = np.logspace(-2, 2, 101)\n",
    "\n",
    "# Model parameter which apply for all\n",
    "model = {\n",
    "    'src': [0, 0, 1e-5, 0, 0],\n",
    "    'rec': [500, 0, 1e-5, 0, 0],\n",
    "    'depth': 0,\n",
    "    'freqtime': times,\n",
    "    'signal': 1,\n",
    "    'verb': 1\n",
    "}\n",
    "\n",
    "# Collect Cole-Cole models\n",
    "res_0 = np.array([2e14, 10])\n",
    "res_8 = np.array([2e14, 5])\n",
    "tau = [0, 1]\n",
    "c = [0, 0.5]\n",
    "m = (res_0-res_8)/res_0\n",
    "\n",
    "cole_model = {'res': res_0, 'cond_0': 1/res_0, 'cond_8': 1/res_8,\n",
    "              'tau': tau, 'c': c, 'func_eta': cole_cole}\n",
    "pelton_model = {'res': res_0, 'rho_0': res_0, 'm': m, \n",
    "                'tau': tau, 'c': c, 'func_eta': pelton_et_al}\n",
    "\n",
    "# Calculate\n",
    "out_bipole = bipole(res=res_0, **model)\n",
    "out_cole = bipole(res=cole_model, **model)\n",
    "out_pelton = bipole(res=pelton_model, **model)\n",
    "\n",
    "# Plot\n",
    "plt.figure()\n",
    "plt.title('Switch-off')\n",
    "plt.plot(times, out_bipole, label='Regular Bipole')\n",
    "plt.plot(times, out_cole, '--', label='Cole and Cole (1941)')\n",
    "plt.plot(times, out_pelton, '-.', label='Pelton et al. (1978)')\n",
    "plt.legend()\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.xlabel('time (s)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table style='border: 3px solid #ddd;'>\n",
       "  <tr>\n",
       "     <td style='text-align: center; font-weight: bold; font-size: 1.2em; border: 2px solid #fff;' colspan='8'>Sat Oct 27 08:41:39 2018 CDT</td>\n",
       "  </tr>\n",
       "  <tr>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>Linux</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>OS</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>4</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>CPU(s)</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>1.15.3</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>numpy</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>1.1.0</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>scipy</td>\n",
       "  </tr>\n",
       "  <tr>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>1.8.0</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>empymod</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>7.0.1</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>IPython</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>2.6.8</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>numexpr</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>3.0.1</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>matplotlib</td>\n",
       "  </tr>\n",
       "  <tr>\n",
       "     <td style='text-align: center; border: 2px solid #fff;' colspan='8'>3.6.6 |Anaconda custom (64-bit)| (default, Oct  9 2018, 12:34:16) \n",
       "[GCC 7.3.0]</td>\n",
       "  </tr>\n",
       "  <tr>\n",
       "     <td style='text-align: center; background-color: #ddd;border: 2px solid #fff;' colspan='8'>Intel(R) Math Kernel Library Version 2019.0.0 Product Build 20180829 for Intel(R) 64 architecture applications</td>\n",
       "  </tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "versions('HTML')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}