{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# solving rate equations for the multi-photon multiple ionization dynamics due to intense x-ray radiation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CFEL 2022, Ludger Inhester"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:19.502304Z",
     "start_time": "2022-05-31T09:13:18.941247Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import solve_ivp\n",
    "from matplotlib import pyplot as plt\n",
    "AU2SEC = 2.418884326500e-17\n",
    "FS2AU = 1e-15 / AU2SEC\n",
    "AU2FS = 1 / FS2AU\n",
    "AU2M = 0.5291772083e-10\n",
    "M2AU = 1. / AU2M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:19.517174Z",
     "start_time": "2022-05-31T09:13:19.504550Z"
    }
   },
   "outputs": [],
   "source": [
    "def solveREQ(pIni,\n",
    "             processData,\n",
    "             configurations,\n",
    "             fluence,\n",
    "             pulseEnv,\n",
    "             tStart=-150 * FS2AU,\n",
    "             tEnd=150 * FS2AU,\n",
    "             dt=0.1 * FS2AU,\n",
    "             method='LSODA'):\n",
    "    \"\"\"\n",
    "    solves the rate equation and returns results for the population\n",
    "    Mandatory arguments: \n",
    "        pIni: initial population, array of shape(nconfigurations)\n",
    "        processData: process data, dictionary obtained from readProcessData\n",
    "        configurations: array of configurations, array of shape (nconfigurations, number of occupations)\n",
    "        fluence: fluence value, scalar (number of photons per area in atomic units)\n",
    "        pulseEnv: pulseEnvelop, function f(t), t in atomic units  \n",
    "    optional arguments:\n",
    "        tStart: time to start integration (i.e., where p=pIni) in atomic units, default -100*FSAU\n",
    "        tEnd: time to end integration in atomic units, default 100*FSAU\n",
    "        dt: timestep in atomic units for which to provide results, default 0.1*FSAU.\n",
    "            The timestep is also used as maximal timestep for the integrator. \n",
    "            It must be equal or smaller than the pulse duration.\n",
    "        method: one of 'RK23', 'RK45', 'BDF', 'LSODA', 'Radau', 'DOP853'\n",
    "    returns:\n",
    "        res: return results object of scipy.integrate.solve_ivp\n",
    "        in particular\n",
    "        res['t'] time values\n",
    "        res['y'] population values at different times\n",
    "    \"\"\"\n",
    "    confDict = {str(c): i for i, c in enumerate(configurations)}\n",
    "    gammas = [{\n",
    "        'i': confDict[str(proc['iconfig'])],\n",
    "        'f': confDict[str(proc['fconfig'])],\n",
    "        'rate': proc['rate']\n",
    "    } for proc in processData if proc['type'] == 'A' or proc['type'] == 'F']\n",
    "    sigmas = [{\n",
    "        'i': confDict[str(proc['iconfig'])],\n",
    "        'f': confDict[str(proc['fconfig'])],\n",
    "        'pcs': proc['pcs']\n",
    "    } for proc in processData if proc['type'] == 'P' or proc['type'] == 'V']\n",
    "\n",
    "    def REQ(t, p):\n",
    "        \"\"\"\n",
    "        returns dp/dt (the right-hand-side of the diff. eq.)\n",
    "        \"\"\"\n",
    "        pDot = np.zeros(p.shape)\n",
    "        Ipulse = pulseEnv(t) * fluence\n",
    "        for s in sigmas:\n",
    "            pDot[s['i']] += -s['pcs'] * Ipulse * p[s['i']]\n",
    "            pDot[s['f']] += +s['pcs'] * Ipulse * p[s['i']]\n",
    "        for g in gammas:\n",
    "            pDot[g['i']] += -g['rate'] * p[g['i']]\n",
    "            pDot[g['f']] += +g['rate'] * p[g['i']]\n",
    "        return (pDot)\n",
    "\n",
    "    def REQjac(t, p):\n",
    "        \"\"\"\n",
    "         returns d/dpj dp_i/dt (jacobian of the right-hand-side of the diff. eq.)\n",
    "        \"\"\"\n",
    "        pDot = REQ(t, p)\n",
    "        pDotjac = np.zeros((p.shape[0], p.shape[0]))\n",
    "        Ipulse = pulseEnv(t) * fluence\n",
    "        for s in sigmas:\n",
    "            pDotjac[s['i'], s['i']] += -s['pcs'] * Ipulse\n",
    "            pDotjac[s['f'], s['i']] += +s['pcs'] * Ipulse\n",
    "        for g in gammas:\n",
    "            pDotjac[g['i'], g['i']] += -g['rate']\n",
    "            pDotjac[g['f'], g['i']] += +g['rate']\n",
    "        return (pDotjac)\n",
    "\n",
    "    times = np.arange(tStart, tEnd, dt)\n",
    "    if method == 'RK45' or method == 'RK23' or method == 'Radau' or method == 'DOP853':\n",
    "        res = solve_ivp(REQ, (tStart - dt, tEnd + dt),\n",
    "                        pIni,\n",
    "                        t_eval=times,\n",
    "                        method=method,\n",
    "                        max_step=dt)\n",
    "    elif method == 'LSODA' or method == 'BDF':\n",
    "        res = solve_ivp(REQ, (tStart - dt, tEnd + dt),\n",
    "                        pIni,\n",
    "                        t_eval=times,\n",
    "                        method='LSODA',\n",
    "                        jac=REQjac,\n",
    "                        max_step=dt)\n",
    "    else:\n",
    "        raise ValueError(f\"method={method} unknown\")\n",
    "    if not res.success:\n",
    "        print(res.message)\n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:19.631531Z",
     "start_time": "2022-05-31T09:13:19.519135Z"
    }
   },
   "outputs": [],
   "source": [
    "def readProcessData(filename):\n",
    "    \"\"\"\n",
    "    reads in and returns process data generated by xatom from file filename\n",
    "    \"\"\"\n",
    "    processData = []\n",
    "    with open(filename, 'r') as f:\n",
    "        lines = f.readlines()\n",
    "        for line in lines:\n",
    "            line = line.rstrip()\n",
    "            els = line.split()\n",
    "            if len(els) == 0:\n",
    "                continue\n",
    "            configs = line.split('#')[1].split('->')\n",
    "            if (els[0] == 'P' or els[0] == 'V'):\n",
    "                processData.append({\n",
    "                    \"type\":\n",
    "                    els[0],\n",
    "                    \"i\":\n",
    "                    els[1],\n",
    "                    \"f\":\n",
    "                    els[2],\n",
    "                    \"pcs\":\n",
    "                    float(els[3]),\n",
    "                    \"energy\":\n",
    "                    els[4],\n",
    "                    \"iconfig\":\n",
    "                    np.array([int(n) for n in configs[0].split()]),\n",
    "                    \"fconfig\":\n",
    "                    np.array([int(n) for n in configs[1].split()]),\n",
    "                })\n",
    "            else:\n",
    "                processData.append({\n",
    "                    \"type\":\n",
    "                    els[0],\n",
    "                    \"i\":\n",
    "                    els[1],\n",
    "                    \"f\":\n",
    "                    els[2],\n",
    "                    \"rate\":\n",
    "                    float(els[3]),\n",
    "                    \"energy\":\n",
    "                    els[4],\n",
    "                    \"iconfig\":\n",
    "                    np.array([int(n) for n in configs[0].split()]),\n",
    "                    \"fconfig\":\n",
    "                    np.array([int(n) for n in configs[1].split()]),\n",
    "                })\n",
    "        return processData"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:19.756117Z",
     "start_time": "2022-05-31T09:13:19.634111Z"
    }
   },
   "outputs": [],
   "source": [
    "def gaussian(x, sigma=10 * FS2AU, x0=0.):\n",
    "    \"\"\"\n",
    "    Gaussian function:\n",
    "    x0 is the center\n",
    "    sigma is the stdev\n",
    "    \"\"\"\n",
    "    return 1. / np.sqrt(2 * np.pi) * 1 / sigma * np.exp(-0.5 *\n",
    "                                                        ((x - x0) / sigma)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial 3-configuration example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:19.890692Z",
     "start_time": "2022-05-31T09:13:19.759189Z"
    }
   },
   "outputs": [],
   "source": [
    "pIni = np.array([1., 0., 0., 0.])\n",
    "sigma = 1.\n",
    "Gamma = 3 / 2\n",
    "processData = [\n",
    "    {\n",
    "        \"iconfig\": np.array([1, 1]),\n",
    "        \"fconfig\": np.array([0, 1]),\n",
    "        \"pcs\": sigma,\n",
    "        \"type\": \"P\"\n",
    "    },\n",
    "    {\n",
    "        \"iconfig\": np.array([0, 1]),\n",
    "        \"fconfig\": np.array([0, 0]),\n",
    "        \"pcs\": sigma,\n",
    "        \"type\": \"P\"\n",
    "    },\n",
    "    {\n",
    "        \"iconfig\": np.array([0, 1]),\n",
    "        \"fconfig\": np.array([1, 0]),\n",
    "        \"rate\": Gamma,\n",
    "        \"type\": \"A\"\n",
    "    },\n",
    "]\n",
    "configurations = [\n",
    "    np.array([1, 1]),\n",
    "    np.array([0, 1]),\n",
    "    np.array([0, 0]),\n",
    "    np.array([1, 0])\n",
    "]\n",
    "fluence = 1.\n",
    "pulseEnv = lambda t: 1. if t > 0 else 0.\n",
    "\n",
    "res = solveREQ(pIni,\n",
    "               processData,\n",
    "               configurations,\n",
    "               fluence,\n",
    "               pulseEnv,\n",
    "               tStart=0,\n",
    "               tEnd=10,\n",
    "               dt=0.1,\n",
    "               method='LSODA')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:20.275680Z",
     "start_time": "2022-05-31T09:13:19.893571Z"
    }
   },
   "outputs": [],
   "source": [
    "times = res['t']\n",
    "p = res['y']\n",
    "# plot populations and compare to analytical solution\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(times, p[0, :], label=r'$p_1$')\n",
    "ax.plot(times, p[1, :], label=r'$p_2$')\n",
    "ax.plot(times, p[2, :], label=r'$p_3$')\n",
    "ax.plot(times, p[3, :], label=r'$p_4$')\n",
    "\n",
    "ax.set_xlabel('time')\n",
    "ax.set_ylabel('population')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Below is the solution for exercise 3a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:20.281772Z",
     "start_time": "2022-05-31T09:13:20.277748Z"
    }
   },
   "outputs": [],
   "source": [
    "def analyticalSolution(t, sigma, Gamma):\n",
    "    p1 = np.exp(-sigma * t)\n",
    "    p2 = sigma / Gamma * (np.exp(-sigma * t) - np.exp(-(sigma + Gamma) * t))\n",
    "    p3 = sigma / (sigma + Gamma) - sigma / Gamma * np.exp(\n",
    "        -sigma * t) + sigma**2 / (Gamma * sigma + Gamma**2) * np.exp(\n",
    "            -(sigma + Gamma) * t)\n",
    "    p4 = 1 - p1 - p2 - p3\n",
    "    return np.array([p1, p2, p3, p4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:20.533888Z",
     "start_time": "2022-05-31T09:13:20.283099Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "pA = analyticalSolution(times, sigma, Gamma)\n",
    "\n",
    "# plot populations and compare to analytical solution\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(times, p[0, :], label=r'$p_1$')\n",
    "ax.plot(times, pA[0, :], 'x', label=r'analytical')\n",
    "\n",
    "ax.plot(times, p[1, :], label=r'$p_2$')\n",
    "ax.plot(times, pA[1, :], 'x', label=r'analytical')\n",
    "\n",
    "ax.plot(times, p[2, :], label=r'$p_3$')\n",
    "ax.plot(times, pA[2, :], 'x', label=r'analytical')\n",
    "\n",
    "ax.plot(times, p[3, :], label=r'$p_4$')\n",
    "ax.plot(times, pA[3, :], 'x', label=r'analytical')\n",
    "\n",
    "ax.set_xlabel('time')\n",
    "ax.set_ylabel('population')\n",
    "ax.legend(loc='upper center', ncol=4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### End of solution for exercise 3a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rate Equation Solving with Calculated Atomic Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# read an preprocess the process data (rates and cross sections calculated by xatom)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:20.543498Z",
     "start_time": "2022-05-31T09:13:20.535029Z"
    }
   },
   "outputs": [],
   "source": [
    "# read in process data: neon, 2000eV\n",
    "processData = readProcessData('results_O_750.dat')\n",
    "\n",
    "# get a list of involved electronic configurations\n",
    "configurations = np.unique(np.array([p['iconfig'] for p in processData] +\n",
    "                                    [p['fconfig'] for p in processData]),\n",
    "                           axis=0)\n",
    "Nc = configurations.shape[0]\n",
    "\n",
    "# compute number of electrons for each configuration\n",
    "Nelectron = np.array([np.sum(c) for c in configurations])\n",
    "\n",
    "# initial population we start with the configuration that has the maximum number of electrons\n",
    "startCfgIdx = np.argmax(Nelectron)\n",
    "pIni = np.zeros(Nc)\n",
    "pIni[startCfgIdx] = 1.\n",
    "\n",
    "# Assume that nuclear charge is identical to maximum number of electron,\n",
    "# i.e., the configuration with the most electrons is neutral\n",
    "nuclearCharge = np.max(Nelectron)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# solve rate equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:21.927178Z",
     "start_time": "2022-05-31T09:13:20.544949Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 534 ms, sys: 801 µs, total: 534 ms\n",
      "Wall time: 536 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# set pulse parameters:\n",
    "# the fluence number of photons per micrometer square\n",
    "nPhotonsPerSqmuM = 1.5e13 \n",
    "# the pulse width (FWHM)\n",
    "fwhmPulse = 25 * FS2AU\n",
    "# convert fluence into atomic units and pulse width FWHM into stdev (atomic units)\n",
    "fluence = nPhotonsPerSqmuM / (1e-6 * M2AU)**2\n",
    "sigmaPulse = fwhmPulse / (2. * np.sqrt(2. * np.log(2.)))\n",
    "# pulse Envelope: simple Gaussian with stddev 25fs and t0 = 0fs\n",
    "pulseEnvelope = lambda t: gaussian(t, sigma=sigmaPulse, x0=0.0 * FS2AU)\n",
    "\n",
    "res = solveREQ(pIni,\n",
    "               processData,\n",
    "               configurations,\n",
    "               fluence,\n",
    "               pulseEnvelope,\n",
    "               method='LSODA')\n",
    "\n",
    "# p is population, times is timepoints\n",
    "p = res['y']\n",
    "times = res['t']\n",
    "\n",
    "# compute the average charge per time\n",
    "avgCharge = nuclearCharge - np.dot(Nelectron, p)\n",
    "\n",
    "# compute the population in a given charge state per time\n",
    "chargePops = np.zeros((nuclearCharge + 1, p.shape[1]))\n",
    "for q in range(chargePops.shape[0]):\n",
    "    idx = np.where(nuclearCharge - Nelectron == q)[0]\n",
    "    chargePops[q, :] = np.sum(p[idx, :], axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# plot/analyze the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:22.158323Z",
     "start_time": "2022-05-31T09:13:21.932125Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlcAAAHACAYAAACPsE3iAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABrk0lEQVR4nO3deXhTVeI+8DdJk3QvdKELW8tuKUsF2YZFBIq4AKIOIrKIOqMjClR03AVE6wII6lcQxCKOUGdUHH+OIousKlspUigiUGihtJRC9yXr/f0RbmhoC1lumu39PE8eyc3NPSdXuHlzzrnnyARBEEBEREREkpC7ugJERERE3oThioiIiEhCDFdEREREEmK4IiIiIpIQwxURERGRhBiuiIiIiCTEcEVEREQkIYYrIiIiIgkxXBERERFJiOGKiIiISEIMV0RERD5m586duPvuuxEXFweZTIZvv/3WpvfX1dVh+vTp6NGjB/z8/DB+/PgG+xQWFuLBBx9E165dIZfLMXv2bEnq7gkYroiIiHxMdXU1evXqhQ8//NCu9xsMBgQEBODpp5/GyJEjG91Ho9EgKioKL730Enr16uVIdT2On6srQERERM1rzJgxGDNmTJOva7VavPzyy/jiiy9QVlaGpKQkvP3227j11lsBAEFBQVi+fDkA4JdffkFZWVmDY8THx2PZsmUAgE8//VTyz+DOGK6IiIjIwsMPP4wzZ84gIyMDcXFx2LBhA26//XZkZ2ejc+fOrq6e22O3IBEREZmdOnUK69evx3/+8x8MGTIEHTt2xNy5czF48GCkp6e7unoegS1XREREZHbw4EEIgoAuXbpYbNdoNIiIiHBRrTwLwxURERGZGY1GKBQKZGZmQqFQWLwWHBzsolp5FoYrIiIiMktOTobBYEBxcTGGDBni6up4JIYrIiIiH1NVVYWTJ0+an58+fRqHDh1CeHg4unTpgsmTJ2Pq1KlYvHgxkpOTUVJSgp9//hk9evTAHXfcAQDIycmBVqvF5cuXUVlZiUOHDgEAevfubT6uuK2qqgoXL17EoUOHoFKpkJiY2Fwf1SVkgiAIrq4EERERNZ/t27dj+PDhDbZPmzYNa9asgU6nw8KFC7F27VoUFBQgIiICAwcOxPz589GjRw8ApqkW8vLyGhyjfqyQyWQNXm/fvj3OnDkj3YdxQwxXRERERBLiVAxEREREEmK4IiIiIpKQRw9o1+v1yMrKQnR0NORy5kQiIiJPYDQaceHCBSQnJ8PPz6OjSKM8+hNlZWWhX79+rq4GERER2WHfvn245ZZbXF0NyXl0uIqOjgZg+p8TGxvr4toQERGRNQoLC9GvXz/z97i38ehwJXYFxsbGok2bNi6uDREREdnCW4f0eOenIiIiInIRhisiIiIiCTFcEREREUnIo8dcWctgMECn07m6GuRCKpXKa/v2iUga/K6QjlKphEKhcHU1XMarw5UgCCgqKkJZWZmrq0IuJpfLkZCQAJVK5eqqEJGb4XeFc7Ro0QIxMTGNri/o7bw6XIn/WFq1aoXAwECf/B9Mpsnqzp8/j8LCQrRr145/D4jIAr8rpCUIAmpqalBcXAwAPjlVkteGK4PBYP7HEhER4erqkItFRUXh/Pnz0Ov1UCqVrq4OEbkJflc4R0BAAACguLgYrVq18rkuQq8dhCL2mwcGBrq4JuQOxO5Ag8Hg4poQkTvhd4XziOfUF8exeW24ErF5lwD+PSCi6+M1Qnq+fE69PlwRERERNSeGK7Lb9OnTMX78eFdXg4iI7HDmzBnIZDIcOnTI1VXxOgxXRERERBLy2rsFyTMJggCDwQA/P/7V9HXN8XdBEATojQIMxiv/NQgQIFx5DVf+ZNpPuLINgGkf85/FfRt/H7k3vVYDvdEIrd4Aud63bnjRXvm8Or3B/OfGyGUy+CnYFmMLfoO5oY0bN2LhwoU4cuQIFAoFBg4ciGXLlqFjx44AgIEDB2LYsGF46623zO+5ePEi4uLisGnTJgwfPhyFhYV49NFH8fPPPyMmJgZvvPEGXnzxRcyePRuzZ8+2ui5Hjx7Fc889h127dkEQBPTu3Rtr1qwx1wUAFi1ahMWLF0Or1eKBBx7A0qVLzdMd/Otf/8LSpUtx/PhxBAUF4bbbbsPSpUvRqlUrAMD27dsxfPhwbNy4ES+99BIOHz6Mn376CX379sXjjz+Ob7/9FqGhoXjuuefw3//+F71798bSpUsBAFqtFi+//DK++OILlJWVISkpCW+//TZuvfVWx/4HkMuVlpZi//790Gg0iI+PR1JSksXgWKNRwIXKOpwpqUFxZR1KqrS4XK3BpSotymt1qNEaUKPVo0ZrQK3WgBqtAVqDEXqDEQajAN2VQGUwMvz4utYhCswb3grGkmrI/DzrrrZH7r8LHbveBAD4YcO/IZcr8NcpM/Dksy9BJpOhV9uWeG/Vv3Db7Xea3zO4e3s8+1oaxv31QRRcrAYAnL5UA3VRJSrKypD2yrP4bec21FRXIzo2Do/MTMX06Q+jXUQgCgoKkJqaik2bNkEul2Pw4MFYtmwZ4uPjXfHx3ZpPhStBEFCra/5fJgFKhU13TVRXVyM1NRU9evRAdXU1Xn31Vdxzzz04dOgQ5HI5Jk+ejHfffRdpaWnm43755ZeIjo7GsGHDAABTp05FSUkJtm/fDqVSidTUVPOEbtYqKCjA0KFDceutt+Lnn39GaGgofvnlF+j1evM+27ZtQ2xsLLZt24aTJ09i4sSJ6N27Nx577DEApgD0+uuvo2vXriguLsacOXMwffp0/PDDDxZlPffcc1i0aBE6dOiAFi1aIDU1Fb/88gu+++47REdH49VXX8XBgwfRu3dv83sefvhhnDlzBhkZGYiLi8OGDRtw++23Izs7G507d7bps5L70Ov1OHDgADQaDQAgN/c0implKDIE49DZchwpKMeZS9XQ6I0uq6NMBsjMf5ZBdmUbAMhgetF375PyLCqF3PT/UyaDvN512hXTttgzF9T/+yoDEx54CF/8v604+nsWFvxzNuLatMN9k6cBMP29lF/z/SNuk1/ZLL/y/P8WvYncE8ex/POv0CI8HPlnTkNTVwuZDKipqcHw4cMxZMgQ7Ny5E35+fli4cCFuv/12HD58mKtfXMOnwlWtzoDEV39q9nJzFoxGoMr6U33vvfdaPF+9ejVatWqFnJwcJCUlYeLEiZgzZw52796NIUOGAADWrVuHBx98EHK5HH/88Qe2bNmC/fv3o2/fvgCATz75xObA8X//938ICwtDRkaGuSWqS5cuFvu0bNkSH374IRQKBbp164Y777wTW7duNYerGTNmmPft0KED3n//ffTr1w9VVVUIDg42v7ZgwQKMGjUKAFBZWYnPPvsM69atw4gRIwAA6enpiIuLM+9/6tQprF+/HufOnTNvnzt3LjZu3Ij09HS8+eabNn1Wch9nz55FTU0t8sr1yC7zw479v6Ncexj+8b0tfqT4yWVo0zIAsWEBiAhWISJIhYhgNVoEKhGgVCBQ5YdAtQKBV/6s8pPDTyGDUi6HQiGDn1wGhdz0Xz+FHH5y05drU8HJl28r92Z1dXU4ffo0EqJD4O/vD8AUrK79Adgc7rjjDpsCVpDaD+3btcXaVR9BJpPh7qF9UVZwCv9eswLznnsaANA+IghJrcPM75HLZGjTMhBJrcMQrAsFAHRqFYKk1mGovlyEgf36YuIdt5p2HtDT/L5PP/0Ucrkcn3zyifnfQnp6Olq0aIHt27cjJSXFwU/vXXwqXHmKU6dO4ZVXXsGePXtQUlICo9H0Cz0/Px9JSUmIiorCqFGj8MUXX2DIkCE4ffo0fvvtNyxfvhwAcPz4cfj5+eHmm282H7NTp05o2bKlTfU4dOgQhgwZct0Zzbt3725xMYiNjUV2drb5eVZWFubNm4dDhw7h8uXLFp8lMTHRvJ8YAgEgNzcXOp0O/fr1M28LCwtD165dzc8PHjwIQRAahD2NRsNZlj1Ync6Aj777FT8eOoMy/1goQqNQpxGglhtwS7Qcg5I6omebFugSHYK4Fv4cB0I+b8CAARbBf+DAgVi8eLFdLW9PPPEE7r33Xhw8eBApKSkYP348Bg0aBADIzMzEyZMnERISYvGeuro6nDp1yrEP4YV8KlwFKBXIWTDaJeXa4u6770bbtm2xatUqxMXFwWg0IikpCVqt1rzP5MmTMWvWLHzwwQdYt24dunfvjl69egFoehCtrYNrxeULrufa4CWTycwBqrq6GikpKUhJScG//vUvREVFIT8/H6NHj7b4LAAQFBTUoJ7XthTUr7/RaIRCoUBmZmaDX3r1W8TIMxiNAv594CwWfX8I+TnHAZkcUW1iMDa5LW4a1gLB2ktIaB+L5OSuNz4YkYMUCgXuuOMOl5QrJZlM1uC6f73Z0seMGYO8vDz873//w5YtWzBixAg8+eSTWLRoEYxGI/r06YMvvviiwfuioqIkrbc38KlwJZPJbOqec4VLly7h2LFj+Pjjj81dfrt3726w3/jx4/H3v/8dGzduxLp16zBlyhTza926dYNer0dWVhb69OkDADh58qTNK7737NkTn332GXQ6nV3r8f3xxx8oKSnBW2+9hbZt2wIADhw4cMP3dezYEUqlEvv27TO/r6KiAidOnDCPKUtOTobBYEBxcbH5PJFnyr1Yhee/zsa+M5ehu1yE8EAVJg3rgbkPpSBI7YdLly7h119/RXFxMQRBYPccNQtPWQtvz549DZ537twZCoUCUVFRKCwsNL924sQJ1NTUXPd4UVFRmD59OqZPn44hQ4bg2WefxaJFi3DzzTfjyy+/RKtWrRAaGuqUz+JN2KbuZlq2bImIiAisXLkSJ0+exM8//4zU1NQG+wUFBWHcuHF45ZVXcOzYMTz44IPm17p164aRI0fib3/7G/bt24esrCz87W9/Q0BAgMUX09SpU/HCCy80WZeZM2eioqICDzzwAA4cOIATJ07g888/x/Hjx636LO3atYNKpcIHH3yA3NxcfPfdd3j99ddv+L6QkBBMmzYNzz77LLZt24ajR49ixowZkMvl5vp36dIFkydPxtSpU/HNN9/g9OnT2L9/P95++22XjJUg+/zvcCHu+mA39p25jECVAlN6t8SbE5IwfWRvBKlNP4TCw8OhVCqh1Wpt/oFA5O3Onj2L1NRUHD9+HOvXr8cHH3yAWbNmAQBuu+02fPjhhzh48CAOHDiAxx9//Lo/lF999VX897//xcmTJ3H06FF8//33uOkm092IkydPRmRkJMaNG4ddu3bh9OnT2LFjB2bNmoVz5841y2f1JAxXbkYulyMjIwOZmZlISkrCnDlz8O677za67+TJk/H7779jyJAhaNeuncVra9euRXR0NIYOHYp77rkHjz32GEJCrg7YBEzjnur/qrlWREQEfv75Z1RVVWHYsGHo06cPVq1aZXUrVlRUFNasWYP//Oc/SExMxFtvvYVFixZZ9d4lS5Zg4MCBuOuuuzBy5Ej85S9/wU033WRR//T0dEydOhXPPPMMunbtirFjx2Lv3r3m1i5yX4IgYNFPx/HkuoOo0RowoEM4fnx6MIa0D4RSobAYNyeTyczPL1++7KoqE7mlqVOnora2Fv369cOTTz6Jp556Cn/7298AAIsXL0bbtm0xdOhQPPjgg5g7d+51F6hWqVR44YUX0LNnTwwdOhQKhQIZGRkATIsw79y5E+3atcOECRNw0003YcaMGaitrWVLViNkggfPcnfu3Dm0bdsWZ8+eRZs2bSxeM98BkpBg8YXsq8RzJfaje5rq6mq0bt0aixcvxiOPPGLz+/n3wX0YjQJe++4oPt+TBwD429AOeG50V1SUl+GXX36BSqXC6NGWYyNPnjyJY8eOITY21uLmByJHefK14dZbb7WY+8/dXO/cXu/72xu49wAkspvY4tSjRw8UFhbiueeeQ3x8PIYOHerqqlklKysLf/zxB/r164fy8nIsWLAAADBu3DgX14wcIQgCXvnvEXyxNx8yGfDmPT0wqZ+p1VVslQoPD2/wPnEbW66IyBMwXHkpnU6HF198Ebm5uQgJCcGgQYPwxRdf2DUw3VUWLVqE48ePQ6VSoU+fPti1axciIyNdXS1ywEfbT5mD1ZK/9sI9yVd/sZaXlwNAo1OGhIWFQSaTQaPRoKam5rpdG0RErsZw5aVGjx7doGvFkyQnJyMzM9PV1SAJfX/4PN79yXQzxLy7u1sEK+BquAoLC2vwXoVCgZCQEFRUVKCiooLhigim5cPIPXFAOxE53emSavzzq8MAgEcHJ2DaoHiL1w0GA6qrTeucNTU4VtxeUVHhvIoSEUnA68OVB4/XJwnx74HraPQGPLX+IKq1BvSLD8fzY7o12EcMTGq1Gmq1utHjMFyRM/EaIT1fPqcuDVfx8fGmtbuueTz55JMOH1scW3SjCdPIN4gzwnvKxIDe5P9+PokjBRVoGajEskm9G12yRgxM17ulW3ytsrLSORUln8TvCucRz6knjfWVikvHXO3fv99i/aMjR45g1KhRuP/++x0+tkKhQIsWLVBcXAzANEcHZ3b2TUajERcvXkRgYCD8/DjMsDmduFCJ5TtM644tHN8DsWGNL6lkTbgS1zSrrq6GwWBgUCZJ8LtCeoIgoKamBsXFxWjRooVP/lt16TfNtesRvfXWW+jYsaN5iRNHxcTEAID5Hw35Lrlcjnbt2vGi2YwEQcCLG7KhMwgYeVMr3NEjpsl9q6qqAKDBorD1+fv7Q6VSQavVoqqqqtGB70T24HeFc7Ro0cJ8bn2N2/yM12q1+Ne//oXU1NQmvwA1Gg00Go35+Y26B2QyGWJjY9GqVavrLlZJ3k+lUkEu9/ohhm5l45Ei7D9TigClAvPHJV032IqD2W+06HZwcDAuX77McEWS4neF9JRKpU+2WIncJlx9++23KCsrw/Tp05vcJy0tDfPnz7f52AqFwqf/JxM1N53BiHeuTLvw2JAEtG7ReHcgYLpTsLa2FoBpzczrCQoKwuXLl81hjEhK/K4gqbjNT/nVq1djzJgxiIuLa3KfF154AeXl5eZHTk5OM9aQiKyVsf8sTpdUIyJIhb8N63jdfcWgpFQqoVKprruvGL4YrojInblFy1VeXh62bNmCb7755rr7XXubNm/JJnI/OoMRy7edBAA8dVsnBKuvf5kRx1vdqEuw/j7ie4iI3JFbtFylp6ejVatWuPPOO11dFSJy0HeHzuN8eR0ig9V44Mq6gdcjtkLdqEuw/j5suSIid+bycGU0GpGeno5p06bxNnkiD2c0Cvh4p2nqhRmD4+GvvPH4FVtarsRwpdPpzHOXERG5G5eHqy1btiA/Px8zZsxwdVWIyEE7/ryIPy9UIUTth4cGtLfqPba0XCkUCvj7+1u8j4jI3bi8qSglJcWnp8gn8iZf7M0DAPz1lrYI9bduVmZxFmdrwpW4X11dHaqrq9GyZUv7KkpE5EQub7kiIu9QUFaLn/8wTcL4YP8bj7UCTNMwiHPXBQYGWvUeMYRxuRIiclcMV0QkiS/35cMoAIM6RqBj1I3HTwEwz2/l5+dn9fpjAQGmObMYrojIXTFcEZHDDEYBGfvPAgAm97durBVwNSBZ22pVf18xmBERuRuGKyJy2G+nLqG4UoMWgUqMSoy2+n1iQBJbo6wh7stwRUTuiuGKiBz230MFAIA7esRC5Wf9ZcXRliveDENE7ojhiogcUqczYOPRIgDAuF5NL1/VGHtarvz9/SGTyWA0Gi0WcicichcMV0TkkO3HL6KyTo/YMH/cEh9u03vtabmSyWTmua7YNUhE7ojhiogc8v3h8wCAu3vFQS6X2fRee1qu6u/POwaJyB0xXBGR3bR6I3YcvwgAuD0pxqb3Go1G1NXVAbCt5ar+/my5IiJ3xHBFRHbbd/oyKjV6RAar0LtNC5veKwYjhUIBlUpl03vZckXkW+bNmweZTGbxiImx7Qddc3L58jdE5Lm2HLsAABjRLbrZugQBtlwR+aLu3btjy5Yt5ucKxY0XhncVhisisosgCNicYwpXtsxtJRK7BMXB6bbgXFdEvsfPz8+tW6vqY7cgEdnlj6JKFJTVwl8px186Rdr8fkfClfge8RhE5P1OnDiBuLg4JCQk4IEHHkBubq6rq9QkhisissuuE6aB7AM7RCBAZXvzvBiM7OkWFMOVTqeDwWCw+f1E5B4qKytRUVFhfjQ1d13//v2xdu1a/PTTT1i1ahWKioowaNAgXLp0qZlrbB2GKyKyyy8nTRc1e1qtgKtdeva0XCmVSvN4C7ZeEXmuxMREhIWFmR9paWmN7jdmzBjce++96NGjB0aOHIn//e9/AIDPPvusOatrNY65IiKbafVG7Dt9GQAwuLN94cqRbkHxfdXV1airq0NQUJBdxyAi18rJyUHr1q3Nz9VqtVXvCwoKQo8ePXDixAlnVc0hbLkiIptl5ZeiVmdAZLAKXaND7DqGFOGq/nGIyPOEhIQgNDTU/LA2XGk0Ghw7dgyxsbFOrqF9GK6IyGa/nCwBAAzqGAmZzLYpGADTnYbi2AqGKyK6kblz52LHjh04ffo09u7di/vuuw8VFRWYNm2aq6vWKHYLEpHNfssVx1tF2PV+jUYDQRAgk8ms/qV6LYYrIt9x7tw5TJo0CSUlJYiKisKAAQOwZ88etG/f3tVVaxTDFRHZRKM34Pdz5QCAfgn2hSsxEKnVartavgCGKyJfkpGR4eoq2ITdgkRkk6PnK6DVGxEepEJ8hG1rAoocmYZBJL6X4YqI3A3DFRHZ5GBeKQDg5nYt7W51cnQwe/33cpZ2InI3DFdEZJOD+aZw1ad9S7uPIWW4EsdvERG5C4YrIrKaIAjIzHM8XDkygahIHAhvNBqh1WrtPg4RkdQYrojIagVltbhQoYGfXIaebcLsPo4ULVdyudwcsJpaMoOIyBUYrojIaofOlgEAbooNhb/S9vUERVKEq/rv57grInInDFdEZLWc8xUAgKTW9rdaAXB4AlERW66IyB0xXBGR1Y5eCVfd40LtPobRaIROpwNg/TpiTWG4IiJ3xHBFRFbLKTSFq0QHwpUYhORyOZRKpUP1qX/HIBGRu2C4IiKrFFfW4WKlBnIZcFOM4+HK0Var+sdguCIid8JwRURWEbsEO0QFI0Bl/2B2hisi8nYMV0RkFXEwe2Ks/a1WgHPCFZfAISJ3wnBFRFY5JsF4K4AtV0Tk/RiuiMgqpy5WAwA6twp26DhiK5OU4Uqv18NgMDh8PCIiKbg8XBUUFOChhx5CREQEAgMD0bt3b2RmZrq6WkRUj8EoIPdiFQCgk4PhSsqWK6VSCblcbnFcIiJX83Nl4aWlpfjLX/6C4cOH48cff0SrVq1w6tQptGjRwpXVIqJrnC+rhUZvhEohR5uWgQ4dS6oJREVqtRq1tbXQaDQIDHSsbkREUnBpuHr77bfRtm1bpKenm7fFx8e7rkJE1KiTxaZWq4TIICjkMoeOJWXLFWAKaWK4IiJyBy7tFvzuu+/Qt29f3H///WjVqhWSk5OxatUqV1aJiBpxSqIuQUD6cMVB7UTkblwarnJzc7F8+XJ07twZP/30Ex5//HE8/fTTWLt2baP7azQaVFRUmB+VlZXNXGMi3yS2XHWMCnLoOAaDAXq9HoD04YrTMRCRu3Bpt6DRaETfvn3x5ptvAgCSk5Nx9OhRLF++HFOnTm2wf1paGubPn9/c1STyeWLLVUeJBrMrFAr4+Ulz+WHLFRG5G5e2XMXGxiIxMdFi20033YT8/PxG93/hhRdQXl5ufuTk5DRHNYl83plLNQBMY64cIeU0DCKGKyJyNy5tufrLX/6C48ePW2z7888/0b59+0b3V6vVFhfliooKp9aPiIBarQEXK03BpV24NHcKShmuuHgzEbkbl4arOXPmYNCgQXjzzTfx17/+Ffv27cPKlSuxcuVKV1aLiOo5V2pqtQpR+yEsQOnQsZwRrthyRUS20ul0KCoqQk1NDaKiohAeHi7p8V3aLXjLLbdgw4YNWL9+PZKSkvD6669j6dKlmDx5siurRUT1nL0SrtqGB0Imc69pGOofi+GKiK6nqqoKH3/8MW699VaEhYUhPj4eiYmJiIqKQvv27fHYY49h//79kpTl0pYrALjrrrtw1113uboaRNSE/EtiuApw+FhSTyAKXA1XBoMBOp0OSqVjrWtE5H3ee+89vPHGG4iPj8fYsWPx/PPPo3Xr1ggICMDly5dx5MgR7Nq1C6NGjcKAAQPwwQcfoHPnznaX5/JwRUTu7WxpLQCgrYMzswPOabkS7zzU6/XQarUMV0TUwK+//opt27ahR48ejb7er18/zJgxAytWrMDq1auxY8cOhisicp6zl692CzrKGeEKAFQqFfR6PTQaDYKCHLujkYi8z3/+8x+r9lOr1fjHP/7hcHkuX7iZiNxb/pVw5eidgsDVcKVSqRw+Vn1iWNNqtZIel4i8T21tLWpqaszP8/LysHTpUvz000+SlcFwRURNEgQB58RuQQnGXInhxxktV/WPT0TUlHHjxplXgikrK0P//v2xePFijB8/HsuXL5ekDIYrImpSRa0eVRrTcjWtWzjWcmU0Gs1L3zir5Yp3DBLRjRw8eBBDhgwBAHz11VeIjo5GXl4e1q5di/fff1+SMhiuiKhJRRWmGdVbBCoRoFI4dCyxVUkmk0m29I2ILVdEZK2amhqEhIQAADZt2oQJEyZALpdjwIAByMvLk6QMhisialJhualLMCbU8akTxOCjUqkcni/rWmy5IiJrderUCd9++y3Onj2Ln376CSkpKQCA4uJihIaGSlIGwxURNenClZaraAnClbMGs9c/JluuiOhGXn31VcydOxfx8fHo378/Bg4cCMDUipWcnCxJGZyKgYiaVFRuCkSxYdK2XEmNLVdEdD2HDx9GUlIS5HI57rvvPgwePBiFhYXo1auXeZ8RI0bgnnvukaQ8tlwRUZOKJGy5ctadgsDVwMZwRUSNSU5ORklJCQCgQ4cOUCqVSE5Ohlx+NQb169cP3bp1k6Q8hisialKROObKQ1qutFotBEGQ/PhE5NlatGiB06dPAwDOnDkDo9Ho1PLYLUhETSqqMLUEST2gXWriMQVBgF6v5xI4RGTh3nvvxbBhwxAbGwuZTIa+fftCoWj8Dujc3FyHy2O4IqImiQPapWi5cuaAdrlcbl5fUKPRMFwRkYWVK1diwoQJOHnyJJ5++mk89thj5ukYnIHhiogapdEbcLna1Nrk7i1XgKlrUFy8mYjoWrfffjsAIDMzE7NmzWK4IqLmV3ylS1DlJ0eLQMdbgpw5oF08bnV1NQe1E9F1paenO70MDmgnokaJdwrGhPpLMumns1uuONcVEbkLhisialRJpakFKCrE8ZYmQRCaLVyx5YqIXI3hiogaVXJlvFVEkONhSK/Xm6dIcOaYK4AtV0TkegxXRNSoS1WmFqBICVquxMDj5+dnMWmflNhyRUTuguGKiBpVIoYrCVquxMDjrMHs9Y/Nlisiqq+2thYFBQUNth89etRpZTJcEVGjLlVd6RYMlq7lylldgvWPzZYrIhJ99dVX6NKlC+644w707NkTe/fuNb82ZcoUp5XLcEVEjRLDVaSHhCu2XBHRtRYuXIiDBw/i999/x6effooZM2Zg3bp1AODUpbI4zxURNUrsFowIlq5bsDlarsT1BaWYPoKIPJtOp0NUVBQAoG/fvti5c6d5pnZnXiPYckVEjTKPuZIgXDVnt6AgCNDpdE4rh4g8R6tWrXD48GHz84iICGzevBnHjh2z2C41hisiakCrN6KiTg8AiAiSrlvQmQPa5XK5eU1Bjrsi8m5paWmQyWSYPXv2dff7/PPP0apVK4ttKpUK69evx44dO5xWP3YLElEDl6pN4cRPLkNYgHRL3ziz5Uo8vk6n47grIi+2f/9+rFy5Ej179rzhvm3atGl0e11dHZRKJb7//nsYjUaL18aOHetwHRmuiKgBcTB7eJAKcrn7L30j4vqCRN6tqqoKkydPxqpVq7Bw4UK7jrFx40ZMmTIFly5davCaTCaDwWBwtJrsFiSihq4OZpemG685BrTXPz5brog8Q2VlJSoqKsyPG/0wevLJJ3HnnXdi5MiRdpc5c+ZM/PWvf0VhYSGMRqPFQ4pgBTBcEVEjSszTMEgThpqzW7B+eUTk3hITExEWFmZ+pKWlNblvRkYGDh48eN19rFFcXIzU1FRER0c7dJzrYbcgETVw+cqYKynWFTQYDOZfg84c0F7/+AxXRJ4hJycHrVu3Nj9v6hpx9uxZzJo1C5s2bYK/v79DZd53333Yvn07Onbs6NBxrofhiogaKK0xTWXQIlC6aRjkcjn8/Jx7yWHLFZFnCQkJQWho6A33y8zMRHFxMfr06WPeZjAYsHPnTnz44YfQaDRQKBRWlfnhhx/i/vvvx65du9CjRw/zXcaip59+2rYP0QiGKyJqoKzGFE5aShiunN0lWL8Mhisi7zJixAhkZ2dbbHv44YfRrVs3/POf/7Q6WAHAunXr8NNPPyEgIADbt2+3mExUJpMxXBGRc5RWm1quWgZ5zjQM9cvg3YJE3iUkJARJSUkW24KCghAREdFg+428/PLLWLBgAZ5//nnI5c4Zes4B7UTUQFmtKRBJ2S3IlisicgdarRYTJ050WrACGK6IqBFl4pgrCSYQFVuRnD2YHWC4IvIl27dvx9KlS21+37Rp0/Dll19KX6F6XNotOG/ePMyfP99iW3R0NIqKilxUIyICgFIPH3MlzldjyzgMIvINBoMB77zzDn766Sf07NmzwYD2JUuWOFyGy8dcde/eHVu2bDE/58WQyLUEQah3t6Bnjbny8/ODXC6H0WiEVqtFQECA08skIs+SnZ2N5ORkAMCRI0csXqs/uN0RLg9Xfn5+iImJcXU1iOiKOp0RWr1prS0pwlVzzc4uUqlUqKurg0ajYbgioga2bdvm9DJcPubqxIkTiIuLQ0JCAh544AHk5uY2ua9Go7GYJr+ysrIZa0rkG8QuQT+5DMFqx39/NWfLVf1yOO6KiGxx8eJFrF+/XpJjubTlqn///li7di26dOmCCxcuYOHChRg0aBCOHj2KiIiIBvunpaU1GKNFRNISw1WLQJUkTeRiyGmOAe0AwxURXd+CBQsa3X7q1Cl88803mDRpksNluDRcjRkzxvznHj16YODAgejYsSM+++wzpKamNtj/hRdesNheUFCAxMTEZqkrka8ol3C8FdD8LVdcAoeIrmfDhg0Wzw0GA86ePYuKigq8/vrrkpTh8jFX9QUFBaFHjx44ceJEo6+r1WqLX78VFRXNVTUinyEOZm8pQbgSBIHdgkTkVrKyshps0+v1mD17NnJyciQpw+VjrurTaDQ4duwYYmNjXV0VIp9Vv1vQUTqdzvzna293dhaGKyKylZ+fH2bPno1vvvlGkuO5NFzNnTsXO3bswOnTp7F3717cd999qKiowLRp01xZLSKfJq4rKMUEomLAUSqVTp0NuT6GKyKyR15eHhISEiQ5lku7Bc+dO4dJkyahpKQEUVFRGDBgAPbs2YP27du7slpEPk2cnb1lkOMtV809DUP9shiuiKgx77//foNtRUVFSE9Px913323xur2LOLs0XGVkZLiyeCJqhDMmEG2uOwUBhisiur733nuv0e3+/v7YvHkzNm/eDMA0oahHhisicj9lHrr0jUgsS2w1IyKq7/Tp004vw60GtBOR65XVSrdosyvDlU6ngyAIzVYuEZGI4YqILJRfCVdhHh6uBEGAXq9vtnKJiEQMV0RkobLOFK5C/KVbV7A5x1zJ5XL4+ZlGPHDcFRG5AsMVEVmorDO19oT4e966giKOuyIiV2K4IiIzvcGIGq0BgHeEK7ZcEZHoxRdfxL59+5qlLIYrIjKr0lwdoyRltyDDFRG5WmFhIe666y7Exsbib3/7G/73v/85rXWb4YqIzMQuQX+lHCo/xy8PbLkiIneRnp6OCxcu4N///jdatGiBZ555BpGRkZgwYQLWrFmDkpISycpiuCIiswoJB7Pr9XoYjUYAzTugvX55DFdEVJ9MJsOQIUPwzjvv4I8//sC+ffswYMAArFq1Cq1bt8bQoUOxaNEiFBQUOFQOwxURmVXUSj+YXS6XQ6FQOHw8W7DlioiscdNNN+G5557DL7/8gnPnzmHatGnYtWsX1q9f79BxOUM7EZlJOQ2DK5a+ETFcEZGtoqKi8Mgjj+CRRx5x+FhsuSIiM3HMVagH3ylYv0yGKyJyBYYrIjK72nLFcEVEZC+7wtWpU6fw8ssvY9KkSSguLgYAbNy4EUePHpW0ckTUvMwTiKo9dxqG+mUyXBGRK9gcrnbs2IEePXpg7969+Oabb1BVVQUAOHz4MF577TXJK0hEzadSI/2AdleOudLpdOY7FomImovN4er555/HwoULsXnzZotfpMOHD8dvv/0maeWIqHk5Y0C7K1qulEolZDKZRT2IiACgtrYWNTU15ud5eXlYunQpNm3aJFkZNoer7Oxs3HPPPQ22R0VF4dKlS5JUiohco8IL1hUETHPZKJVKi3oQEQHAuHHjsHbtWgBAWVkZ+vfvj8WLF2PcuHFYvny5JGXYHK5atGiBwsLCBtuzsrLQunVrSSpFRK5hvlswwLNbruqXy3BFRPUdPHgQQ4YMAQB89dVXiI6ORl5eHtauXYv3339fkjJsDlcPPvgg/vnPf6KoqAgymQxGoxG//PIL5s6di6lTp0pSKSJyDSnvFhQHtLtizBXAcEVEjaupqUFISAgAYNOmTZgwYQLkcjkGDBiAvLw8ScqwOVy98cYbaNeuHVq3bo2qqiokJiZi6NChGDRoEF5++WVJKkVErlFR6x1TMdQvl+GKiOrr1KkTvv32W5w9exY//fQTUlJSAADFxcUIDQ2VpAybr6BKpRJffPEFFixYgKysLBiNRiQnJ6Nz586SVIiIXOfqJKKOdQsajUbodKagxnBFRO7k1VdfxYMPPog5c+ZgxIgRGDhwIABTK1ZycrIkZdj987Rjx47o2LGjJJUgIvdQKdGAdjFYATAPLG9uXLyZiBpz3333YfDgwSgsLESvXr3M20eMGNHoDXv2sPkKmpqa2uh2mUwGf39/dOrUCePGjUN4eLjDlSOi5qMzGFGrMwBwfCqG+l2C4pQIzY0tV0TUlJiYGMTExFhs69evn2THtzlcZWVl4eDBgzAYDOjatSsEQcCJEyegUCjQrVs3fPTRR3jmmWewe/duJCYmSlZRInKuqiutVoDjLVeunJ1dxHBFRE3ZunUrtm7diuLi4gYTDX/66acOH9/mAe3jxo3DyJEjcf78eWRmZuLgwYMoKCjAqFGjMGnSJBQUFGDo0KGYM2eOw5UjouYjdgn6K+VQKhxbdtTVg9nrly0GPSIiAJg/fz5SUlKwdetWlJSUoLS01OIhBZt/nr777rvYvHmzxYj60NBQzJs3DykpKZg1axZeffVV8+h7IvIMFU6Ynd1V0zAAbLkiosatWLECa9aswZQpU5xWhs0/T8vLy82LNdd38eJFVFRUADBNNMoLGpFnuXqnoOdPw1C/bF6LiKg+rVaLQYMGObUMu7oFZ8yYgQ0bNuDcuXMoKCjAhg0b8Mgjj2D8+PEAgH379qFLly5S15WInMhb1hUUiWUbjUYYDAaX1YOI3Mujjz6KdevWObUMm3+ifvzxx5gzZw4eeOAB6PWmX7p+fn6YNm0a3nvvPQBAt27d8Mknn0hbUyJyKqmmYQDcY0C7n58f5HI5jEYjNBoNAgMDXVYXInIfdXV1WLlyJbZs2YKePXs2mC5myZIlDpdh01XUYDAgMzMTb7/9Nt577z3k5uZCEAR07NgRwcHB5v169+7tcMWIqHmJLVeOTiAKuMeYK8AU7urq6qDVahmuiAgAcPjwYXNOOXLkiMVrUk0dY1O4UigUGD16NI4dO4aEhAT07NlTkkoQketVSNhy5Q7dgmL5YrgiIgKAbdu2Ob0Mm6+iPXr0QG5uLhISEpxRHyJyESkXbXaXcMVZ2omoMWVlZVi9ejWOHTsGmUyGxMREzJgxA2FhYZIc366Fm+fOnYvvv/8ehYWFqKiosHgQkWe6OubKOwa01y+f4YqIRAcOHEDHjh3x3nvv4fLlyygpKcGSJUvQsWNHHDx4UJIybP6JevvttwMAxo4da9E3KQgCZDIZ78oh8lBSrisoznjMcEVE7mbOnDkYO3YsVq1aBT8/0/VOr9fj0UcfxezZs7Fz506Hy7D5Kuqsvsq0tDS8+OKLmDVrFpYuXeqUMoioaVJNIioGGYVCAYVC4XC9HMFwRUTXOnDggEWwAkx3Fz/33HPo27evJGXYHK6GDRsmScH17d+/HytXruQAeSIXkmoSUXe5UxDgEjhE1FBoaCjy8/PRrVs3i+1nz55FSEiIJGXYfRWtqalBfn5+g1+EtgakqqoqTJ48GatWrcLChQvtrQ4ROUiqSUTdZbxV/Tqw5YqIRBMnTsQjjzyCRYsWYdCgQZDJZNi9ezeeffZZTJo0SZIybA5XFy9exMMPP4wff/yx0ddtHXP15JNP4s4778TIkSMZrohcSKoxVwxXROTOFi1aBJlMhqlTp5onQ1cqlXjiiSfw1ltvSVKGzXcLzp49G6WlpdizZw8CAgKwceNGfPbZZ+jcuTO+++47m46VkZGBgwcPIi0tzar9NRqNxZ2JlZWVtlafiJpwtVvQsZYrd5idXcRwReQdli9fjp49eyI0NBShoaEYOHBgk408N6JSqbBs2TKUlpbi0KFDyMrKwuXLl/Hee+9JNpzB5p+oP//8M/773//illtugVwuR/v27TFq1CiEhoYiLS0Nd955p1XHOXv2LGbNmoVNmzbB39/fqvekpaVh/vz5tlaZiG5AZzCiVmdqdfamlqv681yJdzQTkedp06YN3nrrLXTq1AkA8Nlnn2HcuHHIyspC9+7dbTpWfn4+2rZti8DAQPTo0aPBa+3atXO4vja3XFVXV6NVq1YAgPDwcFy8eBGAaXJRW+aHyMzMRHFxMfr06QM/Pz/4+flhx44deP/99+Hn59do9+ILL7yA8vJy8yMnJ8fW6hNRI6qutFoBQLAXDWivv2YYW6+IPNfdd9+NO+64A126dEGXLl3wxhtvIDg4GHv27LH5WAkJCebsUt+lS5ckmyDd5qto165dcfz4ccTHx6N37974+OOPER8fjxUrViA2Ntbq44wYMQLZ2dkW2x5++GF069YN//znPxu9hVutVltcsDlpKZE0xGkYApQKKBU2/+ay4E4tV3K5HEqlEjqdDlqt1i0CHxFdVVlZafFdfu33fGMMBgP+85//oLq6GgMHDrS5zKZasauqqqzuSbsRm8PV7NmzUVhYCAB47bXXMHr0aHzxxRdQqVRYs2aN1ccJCQlBUlKSxbagoCBEREQ02E5EziXVYHbAvcZcAaZ6iOGKiNxLYmKixfPXXnsN8+bNa3Tf7OxsDBw4EHV1dQgODsaGDRsavP96UlNTAZgWZ37llVcsFnM3GAzYu3eveUFnR9l8JZ08ebL5z8nJyThz5gz++OMPtGvXDpGRkZJUioiaV4UXrisoUqvVqK6uZrgickM5OTlo3bq1+fn1Wq26du2KQ4cOoaysDF9//TWmTZuGHTt2WB2wsrKyAJharrKzsy2uUSqVCr169cLcuXPt/CSWHL6SBgYG4uabb5aiLti+fbskxyEi25jvFAyQbl1Bd+mC40SiRO4rJCQEoaGhVu2rUqnMA9r79u2L/fv3Y9myZfj444+ter+4wszDDz+MZcuWWV2uPWwOVwaDAWvWrMHWrVtRXFxsXkNM9PPPP0tWOSJqHlIt2mw0Gs3zxrhLyxWnYyDyToIg2PWjKT09HWVlZVi8eDGOHTsGmUyGxMREzJgxA2FhYZLUzeZwNWvWLKxZswZ33nknkpKSeGszkReolKhbULzQiQPJ3UH96RiIyDO9+OKLGDNmDNq2bYvKykpkZGRg+/bt2Lhxo83HOnDgAEaPHo2AgAD069cPgiBgyZIleOONN7Bp0yZJeuNsvpJmZGTg3//+N+644w6HCyci9yD1uoLu0moFsFuQyBtcuHABU6ZMQWFhIcLCwtCzZ09s3LgRo0aNsvlYc+bMwdixYy0Wb9br9Xj00Ucxe/Zs7Ny50+H62nwlrd/nSUTewRvXFRSx5YrI861evVqyYx04cMAiWAGAn58fnnvuOfTt21eSMmye0OaZZ57BsmXLIAiCJBUgItczj7lSe2/LFcMVEQFAaGgo8vPzG2w/e/YsQkJCJCnDqivphAkTLJ7//PPP+PHHH9G9e/cG4yq++eYbSSpGRM1Hqnmu3G2OK4DdgkRkaeLEiXjkkUewaNEiDBo0CDKZDLt378azzz6LSZMmSVKGVVfSa0fP33PPPZIUTkTuoULibkF3mYYBYMsVEVlatGgRZDIZpk6dar67WalU4oknnsBbb70lSRlWhav09HRJCiMi91ThxS1XYtATBAE6nc5t7mIkItdQqVRYtmwZ0tLScOrUKQiCgE6dOlnM2O4om6+kp0+fhl6vR+fOnS22nzhxAkqlEvHx8VLVjYiaiTig3dFJRN1xzJVcLoefnx/0ej20Wi3DFREBME2C3qNHD6cc2+ZwNX36dMyYMaNBuNq7dy8++eQTzrJO5IGkGnPljt2CgCns6fV6aDQaBAUFubo6RORiW7dubXIy9E8//dTh49t8t2BWVhb+8pe/NNg+YMAAHDp0yOEKEVHzM7dcOTjmyh27BQGOuyKiq+bPn4+UlBRs3boVJSUlKC0ttXhIweafqTKZDJWVlQ22l5eXw2AwSFIpImo+OoMRdTrTLzepWq7cLVxxrisiEq1YsQJr1qzBlClTnFaGzS1XQ4YMQVpamkWQMhgMSEtLw+DBgyWtHBE5n9glCADBDsxzJQ4YB9wvXHE6BiISabVaDBo0yKll2HwlfeeddzB06FB07doVQ4YMAQDs2rULFRUVXLSZyAOJXYKBKgX8FDb/3jKr3yrkruGKLVdE9Oijj2LdunV45ZVXnFaGzeEqMTERhw8fxocffojff/8dAQEBmDp1KmbOnInw8HBn1JGInEjqwewqlcrtFnRntyCRb0tNTTX/2Wg0YuXKldiyZQt69uzZ4A7iJUuWOFyeXVfTuLg4vPnmmw4XTkSuJ9UEou46mB1gyxWRr8vKyrJ43rt3bwDAkSNHLLZL9cPQsZ+qROTxnNFy5W445orIt23btq1Zy7N/gAUReYWKWu9d+kbElisiaowgCBAEQfLjMlwR+Tix5SrUC5e+EXHMFRHVt3r1aiQlJcHf3x/+/v5ISkrCJ598Itnx2S1I5OOudgt639I3IrFOBoMBBoMBCoXCxTUiIld55ZVX8N577+Gpp57CwIEDAQC//fYb5syZgzNnzmDhwoUOl8FwReTjrs7O7r1jrvz8/CCXy2E0GqHRaCRdoJWIPMvy5cuxatUqTJo0ybxt7Nix6NmzJ5566ilJwpVk3YIvvvgiZsyYIdXhiKiZSDWgXewWdMcxVwDHXRGRicFgQN++fRts79OnD/R6fSPvsJ1k4aqgoABnzpyR6nBE1EwqNdIOaHfHliuA466IyOShhx7C8uXLG2xfuXIlJk+eLEkZknULfvbZZ1IdioiakdRTMbh7yxWnYyCi1atXY9OmTRgwYAAAYM+ePTh79iymTp1qMeGovROKcswVkY+rkGBAuyAIbt9yxW5BIgJME4fefPPNAIBTp04BAKKiohAVFWUxqagjE4raHK7ef//9RrfLZDL4+/ujU6dOGDp0KO/GIfIQleYZ2u3/raXX681zxbhruGK3IBEBzTOhqM1X0/feew8XL15ETU0NWrZsCUEQUFZWhsDAQAQHB6O4uBgdOnTAtm3b0LZtW2fUmYgkJEW3oNjVJt6V547YckVEzcXmq+Cbb76JW265BSdOnMClS5dw+fJl/Pnnn+jfvz+WLVuG/Px8xMTEYM6cOc6oLxFJTJyhPdSBbkF37xIEOOaKiJqPzT9VX375ZXz99dfo2LGjeVunTp2waNEi3HvvvcjNzcU777yDe++9V9KKEpH0tHojNHojAGnClbsOZgfYLUhEzcfmlqvCwsJG54HQ6/UoKioCAMTFxaGystLx2hGRU4njrQAgWIJuQU9ouWK4IiJnszlcDR8+HH//+9+RlZVl3paVlYUnnngCt912GwAgOzsbCQkJ0tWSiJxCHG8VpFJAIbf/zhhPaLlityARNRebw9Xq1asRHh6OPn36QK1WQ61Wo2/fvggPD8fq1asBAMHBwVi8eLHklSUiaUm1rqC7z84OXK2bXq+H0Wh0cW2IyJvZ3A8QExODzZs3448//sCff/4JQRDQrVs3dO3a1bzP8OHDJa0kETmHFNMwAJ7RLejn5weZTGaek8vf39/VVSIiNzJy5Ejk5uYiNzfX4WPZfEXdsWMHhg0bhm7duqFbt24OV4CIXKfCR9YVBExz8alUKmg0GoYrImpg/PjxuHTpkiTHsvmKOmrUKMTExODBBx/EQw89hKSkJEkqQkTN72rLlfd3CwKm+mk0Go67IqIGZs6cKdmxbA5X58+fR0ZGBtavX4933nkHSUlJeOihh/Dggw+iTZs2Nh1r+fLlWL58uXnB5+7du+PVV1/FmDFjbK0WEdlBqnUFPSlcARzUTuTrtm7diq1bt6K4uLjBGMxPP/3U4ePbPKA9MjISM2fOxC+//IJTp05h4sSJWLt2LeLj4813C1qrTZs2eOutt3DgwAEcOHAAt912G8aNG4ejR4/aWi0isoMYrkID7G+5MhqN0OlMLWAMV0Tk7ubPn4+UlBRs3boVJSUlKC0ttXhIwaGfqwkJCXj++efRq1cvvPLKK9ixY4dN77/77rstnr/xxhtYvnw59uzZg+7duztSNSKyghQD2sVpGGQyGZRKx7oXnY3hiohWrFiBNWvWYMqUKU4rw+5FwH755Rf84x//QGxsLB588EF0794d33//vd0VMRgMyMjIQHV1NQYOHNjoPhqNBhUVFeYHJyolckxFnbRL3ziyinxzYLgiIq1Wi0GDBjm1DJvD1YsvvoiEhATcdtttyMvLw9KlS1FUVIR//etfdo2Vys7ORnBwMNRqNR5//HFs2LABiYmJje6blpaGsLAw86Op/YjIOuZuQQlmZ3f3LkGA4YqIgEcffRTr1q1zahk2X1G3b9+OuXPnYuLEiYiMjHS4Al27dsWhQ4dQVlaGr7/+GtOmTcOOHTsaDU4vvPACUlNTzc8LCgoYsIgcIMUkop4UrjhLOxHV1dVh5cqV2LJlC3r27NlgOMOSJUscLsPmcPXrr786XGh9KpUKnTp1AgD07dsX+/fvx7Jly/Dxxx832FecEV5UUVEhaV2IfI25WzDAu9cVFLHliogOHz6M3r17AwCOHDli8ZpUQxvsvqLm5OQgPz+/wSKoY8eOdahCgiDwwkfUTCpqHZ/nypNarsQ6arVaCILg9mPEiEh627Ztc3oZNoer3Nxc3HPPPcjOzjYvJQFcTXsGg8HqY7344osYM2YM2rZti8rKSmRkZGD79u3YuHGjrdUiIjtcHXPlG+FKbF0TBAE6nc4jWtuIyPPYPKB91qxZSEhIwIULFxAYGIijR49i586d6Nu3L7Zv327TsS5cuIApU6aga9euGDFiBPbu3YuNGzdi1KhRtlaLiGwkCIK5W9CRqRg8KVzJ5XKOuyLyQfn5+TbtX1BQ4FB5Noer3377DQsWLEBUVBTkcjnkcjkGDx6MtLQ0PP300zYda/Xq1Thz5gw0Gg2Ki4uxZcsWBiuiZqLRG6EzmFqeHZlEVBwa4AnhCuC4KyJfdMstt+Cxxx7Dvn37mtynvLwcq1atQlJSEr755huHyrP556rBYEBwcDAA02zt58+fR9euXdG+fXscP37cocoQUfMRx1vJZUCQSmH3cTyp5Qow1bOyspLhisiHHDt2DG+++SZuv/12KJVK9O3bF3FxcfD390dpaSlycnJw9OhR9O3bF++++67Dy/DZHK6SkpJw+PBhdOjQAf3798c777wDlUqFlStXokOHDg5VhoiaT0W9aRjsHdhd/wYUTxm/xJYrIt8THh6ORYsWYeHChfjhhx+wa9cunDlzBrW1tYiMjMTkyZMxevRoJCUlSVKezeHq5ZdfRnV1NQBg4cKFuOuuuzBkyBBERETgyy+/lKRSROR8UkzDoNfrzTe1eFLLFcBwReSL/P39MWHCBEyYMMGp5dh8VR09erT5zx06dEBOTg4uX76Mli1b8rZmIg9inoZB7fidgkqlEnK53atpNSsOaCciZ3No4WZReHi4FIchomZknobBRyYQFbHlioiczTN+ahKR5K5Ow+Abc1yJGK6IyNkYroh8lK9NICpiuCLyPGlpabjlllsQEhKCVq1aYfz48W49QwHDFZGPurr0jW9MICqqvwQOEXmGHTt24Mknn8SePXuwefNm6PV6pKSkmG+wczeSjLkiIs9zdcyVb7ZcGY1G6HQ6KJX2f34iah7XLouXnp6OVq1aITMzE0OHDrXpWNOnT8eMGTNsfp8t2HJF5KPMUzE40HLlabOzA6YlcMRAxa5BIs9UXl4OwL4b6iorK5GSkoLOnTvjzTffdHipm8YwXBH5KLFb0JExV3V1dQA8K1wBnI6ByF1UVlaioqLC/LDm36QgCEhNTcXgwYPtmvTz66+/RkFBAWbOnIn//Oc/iI+Px5gxY/DVV19Bp9PZ8zEaYLgi8lFSTMUghit/f39J6tRcOKidyD0kJiYiLCzM/EhLS7vhe2bOnInDhw9j/fr1dpcbERGBWbNmISsrC/v27UOnTp0wZcoUxMXFYc6cOThx4oTdxwY45orIZzk6FYMgCB7ZLQgwXBG5i5ycHLRu3dr8/EbXkqeeegrfffcddu7ciTZt2jhcfmFhITZt2oRNmzZBoVDgjjvuwNGjR5GYmIh33nkHc+bMseu4DFdEPsrRqRh0Oh2MRiMAhisisk9ISAhCQ0NvuJ8gCHjqqaewYcMGbN++HQkJCXaXqdPp8N133yE9PR2bNm1Cz549MWfOHEyePBkhISEAgIyMDDzxxBMMV0RkG0enYhC7BFUqlccsfSMSuzEZrog8w5NPPol169bhv//9L0JCQlBUVAQACAsLQ0BAgE3Hio2NhdFoxKRJk7Bv3z707t27wT6jR49GixYt7K4vwxWRD9IbjKjWGgDYPxWDJ07DIBLrLAZEInJvy5cvBwDceuutFtvT09Mxffp0m441a9YsPPPMMwgMDLTYLggCzp49i3bt2qFly5Y4ffq03fX1rJ+bRCSJKo3e/GdHW648bTA7cLXODFdEnkEQhEYftgYrAJg3bx6qqqoabL98+bJD3Y31MVwR+aCKWlO4ClAqoFTYdxnw5JYrdgsS+S5BEBrdXlVVJdmPRXYLEvkg8wSiDkzDIAYTT2y5qj+g3Wg0etyYMSKyXWpqKgBAJpPh1VdftegWNBgM2Lt3b6Pjr+zBcEXkgxydhgHw3AlEAdMgfJlMZp5OwhMDIhHZJisrC4Cp5So7O9s8mTBguib06tULc+fOlaQshisiH3R1Ggbfm0AUMP1yVavVqKurQ11dnUd+BiKyzbZt2wAADz/8MJYtW2bVFBD2Yrgi8kFXp2HwrUWb6/P39zeHKyLyHenp6U4vg+GKyAdVmJe+cbxb0FNbfTionch3pKam4vXXX0dQUJB57FVTlixZ4nB5DFdEPqi8xrRsTQs7w5Ver4fBYJony5NbrgBOx0DkC7KyssyLMotjrxojk8kkKY/hisgHlV/pFmwR6NgEogqFAn5+nnkZ4USiRL5DHG917Z+dhfcfE/mgsivhKszOlitP7xIE2HJF5Ktqa2tRU1Njfp6Xl4elS5di06ZNkpXBcEXkg8pqHAtXnj6YHWC4IvJV48aNw9q1awEAZWVl6NevHxYvXoxx48aZl9lxFMMVkQ8qZ8sVB7QT+aiDBw9iyJAhAICvvvoKMTExyMvLw9q1a/H+++9LUgbDFZEPujrmSnWDPRvnybOzi+rP0t7UchhE5H1qamoQEhICANi0aRMmTJgAuVyOAQMGIC8vT5IyGK6IfJCjA9o9eXZ2kThLO8DWKyJf0qlTJ3z77bc4e/YsfvrpJ6SkpAAAiouLJZtYlOGKyMcYjQLKrkzF4MvdguIs7QDHXRH5kldffRVz585FfHw8+vfvj4EDBwIwtWIlJydLUoZn3kNNRHar0uphvNILZm+4qq2tBQAEBARIVS2X4CztRL7nvvvuw+DBg1FYWIhevXqZt48YMQL33HOPJGUwXBH5mPIrdwqq/eTwVyrsOoY3tFwBHNRO5KtiYmIQExNjsa1fv36SHZ/hisjHODreSqfTmWdn9/RwxW5BIt+0detWbN26FcXFxTAajRavffrppw4f36VjrtLS0nDLLbcgJCQErVq1wvjx43H8+HFXVonI60k1DYNSqYRCYV/Ll7vgXFdEvmf+/PlISUnB1q1bUVJSgtLSUouHFFzacrVjxw48+eSTuOWWW6DX6/HSSy8hJSUFOTk5CAoKcmXViLyWOIFoiwD7pmEQg4inj7cCGK6IfNGKFSuwZs0aTJkyxWlluDRcbdy40eJ5eno6WrVqhczMTAwdOtRFtSLybmW1V+4UtLNbUBzM7uldgsDVgMhwReQ7tFotBg0a5NQy3GoqhvLycgBAeHi4i2tC5L04O/tV4mcQAyMReb9HH30U69atc2oZbjOgXRAEpKamYvDgwUhKSmp0H41GY3FXT2VlZXNVj8hrlJu7BX17Ggbg6mfQ6XTQ6/Xw83ObSyIROUldXR1WrlyJLVu2oGfPnlAqLa+FS5YscbgMt7mSzJw5E4cPH8bu3bub3CctLQ3z589vxloReR9HF232ppYrPz8/KJVK6HQ61NXVITg42NVVIiInO3z4MHr37g0AOHLkiMVr4qoNjnKLcPXUU0/hu+++w86dO9GmTZsm93vhhReQmppqfl5QUIDExMTmqCKR15Bq6RtvaLkCTCFRp9OhtraW4YrIB2zbts3pZbh0zJUgCJg5cya++eYb/Pzzz0hISLju/mq1GqGhoeaHuPAiEVnv6oB2++4W9KYB7cDVkMhxV0S+Y9euXXjooYcwaNAgFBQUAAA+//zz6/ae2cKl4erJJ5/Ev/71L6xbtw4hISEoKipCUVERL3JETlReqwdgX7egwWCATmdq+fKWcMXpGIh8y9dff43Ro0cjICAABw8eNI/lrqysxJtvvilJGS4NV8uXL0d5eTluvfVWxMbGmh9ffvmlK6tF5NXKryzabM+AdjGAKBSKBoNAPRVbroh8y8KFC7FixQqsWrXK4jo2aNAgHDx4UJIyXDrmShAEVxZP5JNKa+wfc+Vt460AhisiX3P8+PFG59IMDQ1FWVmZJGW41TxXRORctVoDanWmdQHDg2wfc+VNdwqKGK6IfEtsbCxOnjzZYPvu3bvRoUMHScpguCLyIZevdAkqFTIEq21vuPa2wewAZ2kn8jV///vfMWvWLOzduxcymQznz5/HF198gblz5+If//iHJGW4xVQMRNQ8SqtN4So8SGXXfC41NTUAgMDAQEnr5UpiUNTr9dDpdF4zloyIGvfcc8+hvLwcw4cPR11dHYYOHQq1Wo25c+di5syZkpTBcEXkQy5fCVctHZyGwZvGXImD88W5rhiuiLzfG2+8gZdeegk5OTkwGo1ITEyUdJ47dgsS+ZDL9Vqu7CG2XHlTuALYNUjkS/Lz8yEIAgIDA9G3b1/069fPHKzy8/MlKYPhisiHOBquxJYrb+oWBDionciXJCQk4OLFiw22X7p06YaTmVuL4YrIh5TW2B+utFotDAbTnYbe2nLFcEXk/QRBaHTMaVVVlWQ363DMFZEPueRAy5UYPNRqNeRy7/pdJl5QGa6IvJe4NrFMJsMrr7xi0QJvMBiwd+9e84LOjmK4IvIhpQ6EK2+8U1AkfiaGKyLvlZWVBcDUcpWdnQ2V6up1UKVSoVevXpg7d64kZTFcEfkQR+4W9MY7BUViuBIDJBF5n23btgEAHn74YSxbtgyhoaFOK4vhisiHiOEqwoFuQW8OV7W1tTAajV7X7UlEV6Wnpzu9DIYrIh8iDmhvyW5BC2q1GgqFAgaDAbW1tQgKCnJ1lYjIibZu3YqtW7eiuLgYRqPR4rVPP/3U4ePz5xmRjzAaBfOizY4MaPfGlivg6udi1yCRd5s/fz5SUlKwdetWlJSUoLS01OIhBbZcEfmIijodDEYBgH1jrrx1AlFRYGAgqqqqGK6IvNyKFSuwZs0aTJkyxWllsOWKyEeI461C1H5Q+dn2T19cdw/w7nAFsOWKyNtptVoMGjTIqWUwXBH5CPOdgg6Mt1IqlV679h7DFZFvePTRR7Fu3TqnlsFuQSIfUVJ15U7BYPvDlTcP9Ga4IvINdXV1WLlyJbZs2YKePXs2+MG4ZMkSh8tguCLyERerNACAqGC1ze+trq4GwHBFRJ7v8OHD5pnYjxw54pQyGK6IfERJ5ZVwFWJ/uPLGaRhE4mfTarXQ6/Xw8+PlkcgbiZOJOhOvHkQ+Qmy5irSj5coXugXF8WQ6nQ61tbUICQlxdZWISCKpqal4/fXXERQUZF5jsDEymQyLFy92uDyGKyIfcZEtVzcUGBiI8vJy1NTUMFwReZGsrCzzHc/iGoONkclkkpTHcEXkI+wNV0aj0TyBqDe3XAFXw5UYJonIPezcuRPvvvsuMjMzUVhYiA0bNmD8+PFWv79+V2BzdAtyKgYiH2FvuKqtrYUgCJDL5VCrbW/18iRieOSgdiL3Ul1djV69euHDDz90dVWswpYrIh8gCAJK7LxbsP54K6mazN2VGK6qqqpcXBMiqm/MmDEYM2aMq6thNYYrIh9QqdFDozctTmpry5UvTMMgCg4OBgB2CxI1k8rKSlRUVJifq9Vqr2ghZ7cgkQ8QuwRD1H7wVypseq+vDGYHLLsFjUaji2tD5P0SExMRFhZmfqSlpbm6SpJgyxWRD3DkTkFfmIZBpFar4efnB71ej+rqat4xSORkOTk5aN26tfm5N7RaAQxXRD5BDFeRdoQrcfyRL4QrwNQ1WFZWxnBF1AxCQkIQGhrq6mpIjt2CRD7A3HJl42B2o9Fo7hYUxyN5Ow5qJyJHseWKyAdcqKwDYHu3YE1NDQRBgEKhgL+/vzOq5nY4qJ3I/VRVVeHkyZPm56dPn8ahQ4cQHh6Odu3aubBmjWO4IvIBReWmcBXXwraAJLbeBAcHe/00DCKx5Yrhish9HDhwAMOHDzc/F5ewmTZtGtasWeOiWjWN4YrIBxReCVcxYQE2va9+uPIV7BYkcj+33norBEFwdTWsxjFXRD5AbLmKDbO/5cpXiOFKo9FAr9e7uDZE5IkYroi8nCAI5nAVE8pwdSNKpdI8vqyystLFtSEiT8RwReTlLldroTWYJsSMZriyijgFA8MVEdnDpeFq586duPvuuxEXFweZTIZvv/3WldUh8krieKvIYDVUftb/k9doNNDpdAAYroiIbOHScOVpq1wTeSJ7x1uJwSIwMBByuW81cjNcEZEjXHq3oKetck3kiQorxDsF7QtX3jh78o2In7n+grJERNbyqKkYNBoNNBqN+Tl/VRLdWFF5LQDbW67EYOGL4UrsBtVoNNBqtVCpVC6uERF5Eo9q609LS7NYPTsxMdHVVSJye1fnuGK4spafnx8CAwMB8EccEdnOo8LVCy+8gPLycvMjJyfH1VUicnvnSk0tV61bWD+BqCAIPt0tCHDcFRHZz6O6BdVqNdTqq2ujcTwE0Y2dvVwDAGgbHmj1e2pqamAwGKBQKMwtOL4mJCQEFy5cYLgiIpt5VMsVEdlGozeg6MqA9nY2hCvxh0tISIjPrCl4LQ5qJyJ7ubTlytNWuSbyNAWltRAEIFClQESQ9YOyfXm8lUj87OXl5RAEwWdDJhHZzqXhytNWuSbyNPlil2DLQJvCAcOV6Y5BhUIBg8GAqqoq8xgsIqIbcWm48rRVrok8zdkrg9ltGW8FAGVlZQB8O1zJZDKEhYXh8uXLKC8vZ7giIqtxzBWRFxMHs9sy3qqurg51dXWQyWRo0aKFk2rmGcLCwgCYugaJiKzFcEXkxfIviXcKWj8Ng9hqFRISAoVC4YxqeQwxXInnhIjIGgxXRF7sbKntLVdikPD1Vivg6jkQB7UTEVmD4YrISwmCgDMl1QCA9hEMV/aoP6i9urra1dUhIg/BcEXkpYoq6lCtNcBPLkP7iCCr38dwdZU4qB0ASktLXVwbIvIUDFdEXupUsamlpV1EIJQK6/6pV1dXQ6fTQS6X8+64K8LDwwEAly9fdnFNiMhTMFwReamTxaZlWzpFBVv9HjFAhIWFQS7n5QFguCIi2/HqSeSlTl00tVx1bGV7uIqIiHBKnTyRGK6qqqqg1WpdXBsi8gQMV0Re6tTFKgBARxtari5dugSA4ao+pVJp7iJl6xURWYPhishLnSwWw5V1g9nr6urMd8S1bNnSafXyROwaJCJbMFwReaGSKg2KKzWQyYAu0dYNTK8/3kqpVDqzeh5HDFdiyx4R0fUwXBF5oZzzpoWXEyKCEKS2bgnRkpISAFeDBF0VGRkJwDRNhU6nc3FtiMjdMVwReaGcQlO4uinO+oWXL168CACIiopySp08mb+/v3nclXieiIiawnBF5IXElqvuVoarqqoq1NTUQC6XczB7E8TQyXBFRDfCcEXkhY6eLwcAJMZaF66Ki4sBmLoE/fys60b0Na1atQLAcEVEN8ZwReRlKut0yL2ypmCilS1XYmAQAwQ1FB4eDrlcjtraWlRWVrq6OkTkxhiuiLzMobNlEASgbXgAWoX433B/vV5vHszOcNU0hUJhHtheVFTk4toQkTtjuCLyMpl5pgWG+7Szbq6qCxcuwGg0IigoiOsJ3kBsbCwAoLCw0MU1ISJ3xnBF5GXM4SreuikVzp8/DwCIi4tzWp28RUxMDGQyGcrLy80TrhIRXYvhisiLGIwCsvLLAFjXcqXX682D2RmubkylUpm7Btl6RURNYbgi8iKHzpahSqNHqL8fusbcuIuvqKjI3CUYGmr9nFi+TOwaLCgocHFNiMhdMVwReZGdf5ru+hvcORIKueyG++fn5wMA2rRp49R6eZO4uDjI5XJUVFSgrKzM1dUhIjfEcEXkRXaeMIWroZ1vPMt6dXU1Ll26BJlMhnbt2jm7al5DqVSaW6/EcEpEVB/DFZGXKK3W4vezZQCAoV1uHK7y8vIAmKZf8Pe/8ZQNdJUYRgsKCqDX611cGyJyNwxXRF5i49EiGAXTrOxxLQKuu69OpzOHq/bt2zdH9bxKZGQkgoODodfrzeeRiEjEcEXkJb4/bJpS4a5esTfcNy8vD3q9HiEhIZw41E4dO3YEAOTm5sJoNLq4NkTkThiuiLxAcUUdfjt1CQBwV4/rT6mg1+uRm5sLAOjUqRNkshsPfKeG2rRpA39/f9TV1eHs2bOurg4RuRGGKyIvsH7fWRgFoE/7lmgXEXjdfU+dOgWNRoPAwEDObeUAuVyOTp06AQCOHz/OsVdEZMZwReThdAYj1u0zjfuZOvD646dqa2tx6tQpAEBiYiLkcl4CHNG+fXsEBQVBo9HgxIkTrq4OEbkJXlmJPNyGgwW4UKFBZLAKtyfFXHff33//HQaDAeHh4ebpBMh+crkciYmJAEwtguXl5S6uERG5A4YrIg+m0RuwbKupxeTxYR2h9lM0ue/p06dx8eJFyOVy9OrVq7mq6PViYmIQFxcHQRBw8OBBGAwGV1eJiFyM4YrIg3207RQKymoRHarGQwOa7hIsKSnB0aNHAQA33XQTgoODm6uKPqFHjx5Qq9WoqqpCVlYWBEFwdZWIyIUYrog81OFzZfho+0kAwCt3JcJf2XirVVlZGQ4cOABBENC6dWt06NChOavpE1QqFfr27Qu5XI7CwkIcOXKEAYvIhzFcEXmgwvJa/P3zTOgMAm7vHoM7ezQ+furixYv47bffoNPpEB4ezu5AJwoPD0fv3r0BAGfOnMGhQ4fYRUjko1werj766CMkJCTA398fffr0wa5du1xdJSK3drqkGn/9+DcUltehQ1QQ3rm/Z4O5qgwGA/744w/s2bMHer0eERER6N+/PxSKpsdkkeNat26Nm2++GTKZDOfOncPu3btRUVHh6moReQ1PyQwuDVdffvklZs+ejZdeeglZWVkYMmQIxowZw8VQiRqhNxjxxd483PX+Lpy9XIv2EYFYO6MfQv2V5n0MBgPy8vKwbds289QA7du3x4ABA+Dn5+eqqvuU1q1bo3///lCpVKioqMCOHTtw6NAhVFZWurpqRB7NkzKDTHDhwID+/fvj5ptvxvLly83bbrrpJowfPx5paWk3fP+5c+fQtm1bnD17Fm3atHFmVYlcQhAEnLpYjZ+OFuHfB84i71INAKBfQjg+nJSMVqH+qK2tRWlpKS5cuIALFy5Ap9MBAAICAtC9e3dOueAidXV1yMnJQUFBgXlby5YtERMTg8jISISGhnKeMfJZ9nx/O5oZmpPLfspqtVpkZmbi+eeft9iekpKCX3/91UW1Mqms1aCotKrBdsFofQ5tLLIKEBod5NpUvG1sc1NZWGjk9aZre/06WHOMxqtR7331j9HYuWjyM1t3fpo8DzYe15afFo2VaWziDN3o/FgcQwCq6nSo0uhRVadHRZ0O58vqcLa0GicvVOFSZS0EowEQjAhVyTH5ljikdFXhRHYmsqqqoNVqLY4XGBiIhIQEtG/fnt2ALuTv74+bb74ZCQkJOHXqFIqKilBaWorS0lIApjmygoODERAQgICAAPj7+0OpVMLPz8/8kMvlkMlkkMlkFn+2Zskie5Y14lJI1BiFQgGVSuXSOrhzZmiMy8JVSUkJDAYDoqOjLbZHR0ejqKio0fdoNBpoNBrzc2c1s3+16wj+ufK/Tjk2ka2Ucjk6RQdhQEIE+saHwV9Zh6LC8+bXZTIZQkJCEBUVhejoaISHh/NL0o20bNkSffv2RV1dHYqKinDhwgWUlpZCp9OhoqKCY7LI7YljCZ2hsrLS4t+AWq2GWq1usJ89mcGVXD4I49ovAUEQmvxiSEtLw/z5851eJ4VcDpWy8VPTWM2u+zV2zWeRNbF3k8eQNXaMpne19sDmzfWObVfdGhzG+i/1BnteOUhj//ubqpv1x7btV3mT51hu/d+AJo8hu/a5DP5KOQKUCgSq/BCgUiAqWI1Wof6ICwtAp9hQBKrV8PPzM/+CFFs7AgMD2b3kIfz9/REfH4/4+HgApqWIKioqUFdXh9raWtTV1UGv11s8jEYjBEEwP+o/t4Ytoz44dQQ1xZk/1sQVDkSvvfYa5s2bZ3VdrpcZXMll4SoyMhIKhaJB4iwuLm6QTEUvvPACUlNTzc8LCgoa/I+RwtSRyZg6Mlny4xIRicSATOTLcnJy0Lp1a/PzxlqtAPsygyu57OeuSqVCnz59sHnzZovtmzdvxqBBgxp9j1qtRmhoqPkREhLSHFUlIiIiJwgJCbH4Xm8qXNmTGVzJpd2CqampmDJlCvr27YuBAwdi5cqVyM/Px+OPP+7KahEREZGb8aTM4NJwNXHiRFy6dAkLFixAYWEhkpKS8MMPP6B9+6bXSCMiIiLf40mZwaXzXDmK81wRERF5Hm///uYtRkREREQSYrgiIiIikhDDFREREZGEGK6IiIiIJMRwRURERCQhhisiIiIiCTFcEREREUmI4YqIiIhIQgxXRERERBJiuCIiIiKSkEvXFnSU0WgEABQWFrq4JkRERGQt8Xtb/B73Nh4dri5cuAAA6Nevn4trQkRERLa6cOEC2rVr5+pqSM6jF27W6/XIyspCdHQ05HJpezgrKyuRmJiInJwchISESHpsb8NzZT2eK+vxXFmP58o2PF/Wc9a5MhqNuHDhApKTk+Hn59HtPI3y6HDlTBUVFQgLC0N5eTlCQ0NdXR23xnNlPZ4r6/FcWY/nyjY8X9bjubIPB7QTERERSYjhioiIiEhCDFdNUKvVeO2116BWq11dFbfHc2U9nivr8VxZj+fKNjxf1uO5sg/HXBERERFJiC1XRERERBJiuCIiIiKSEMMVERERkYR8Ply98cYbGDRoEAIDA9GiRYtG95HJZA0eK1assNgnOzsbw4YNQ0BAAFq3bo0FCxbA24azWXOu8vPzcffddyMoKAiRkZF4+umnodVqLfbxhXPVmPj4+AZ/j55//nmLfaw5f77io48+QkJCAvz9/dGnTx/s2rXL1VVyuXnz5jX4OxQTE2N+XRAEzJs3D3FxcQgICMCtt96Ko0ePurDGzWfnzp24++67ERcXB5lMhm+//dbidWvOjUajwVNPPYXIyEgEBQVh7NixOHfuXDN+iuZxo3M1ffr0Bn/PBgwYYLGPr5wre/l8uNJqtbj//vvxxBNPXHe/9PR0FBYWmh/Tpk0zv1ZRUYFRo0YhLi4O+/fvxwcffIBFixZhyZIlzq5+s7rRuTIYDLjzzjtRXV2N3bt3IyMjA19//TWeeeYZ8z6+cq6asmDBAou/Ry+//LL5NWvOn6/48ssvMXv2bLz00kvIysrCkCFDMGbMGOTn57u6ai7XvXt3i79D2dnZ5tfeeecdLFmyBB9++CH279+PmJgYjBo1CpWVlS6scfOorq5Gr1698OGHHzb6ujXnZvbs2diwYQMyMjKwe/duVFVV4a677oLBYGiuj9EsbnSuAOD222+3+Hv2ww8/WLzuK+fKbgIJgiAI6enpQlhYWKOvARA2bNjQ5Hs/+ugjISwsTKirqzNvS0tLE+Li4gSj0ShxTV2vqXP1ww8/CHK5XCgoKDBvW79+vaBWq4Xy8nJBEHzvXNXXvn174b333mvydWvOn6/o16+f8Pjjj1ts69atm/D888+7qEbu4bXXXhN69erV6GtGo1GIiYkR3nrrLfO2uro6ISwsTFixYkUz1dA9XHvNtubclJWVCUqlUsjIyDDvU1BQIMjlcmHjxo3NVvfm1tj327Rp04Rx48Y1+R5fPVe28PmWK2vNnDkTkZGRuOWWW7BixQqLlbx/++03DBs2zGIekNGjR+P8+fM4c+aMC2rrGr/99huSkpIQFxdn3jZ69GhoNBpkZmaa9/Hlc/X2228jIiICvXv3xhtvvGHR5WfN+fMFWq0WmZmZSElJsdiekpKCX3/91UW1ch8nTpxAXFwcEhIS8MADDyA3NxcAcPr0aRQVFVmcN7VajWHDhvn8ebPm3GRmZkKn01nsExcXh6SkJJ88f9u3b0erVq3QpUsXPPbYYyguLja/xnN1Y963WqITvP766xgxYgQCAgKwdetWPPPMMygpKTF36RQVFSE+Pt7iPdHR0ebXEhISmrvKLlFUVGT+3KKWLVtCpVKhqKjIvI+vnqtZs2bh5ptvRsuWLbFv3z688MILOH36ND755BMA1p0/X1BSUgKDwdDgXERHR/vUeWhM//79sXbtWnTp0gUXLlzAwoULMWjQIBw9etR8bho7b3l5ea6ortuw5twUFRVBpVKhZcuWDfbxtb93Y8aMwf3334/27dvj9OnTeOWVV3DbbbchMzMTarWa58oKXtly1digz2sfBw4csPp4L7/8MgYOHIjevXvjmWeewYIFC/Duu+9a7COTySyeC1cGaF+73d1Ifa4a+7yCIFhs99Rz1Rhbzt+cOXMwbNgw9OzZE48++ihWrFiB1atX49KlS+bjWXP+fEVjf0988TzUN2bMGNx7773o0aMHRo4cif/9738AgM8++8y8D89b0+w5N754/iZOnIg777wTSUlJuPvuu/Hjjz/izz//NP99a4ovnqumeGXL1cyZM/HAAw9cd59rW09sMWDAAFRUVODChQuIjo5GTExMg7QuNqFe+0vJ3Uh5rmJiYrB3716LbaWlpdDpdObz4MnnqjGOnD/x7puTJ08iIiLCqvPnCyIjI6FQKBr9e+JL58EaQUFB6NGjB06cOIHx48cDMLXAxMbGmvfheYP5jsrrnZuYmBhotVqUlpZatMgUFxdj0KBBzVthNxMbG4v27dvjxIkTAHiurOGVLVeRkZHo1q3bdR/+/v52Hz8rKwv+/v7m6QgGDhyInTt3Woyf2bRpE+Li4hwKcc1BynM1cOBAHDlyBIWFheZtmzZtglqtRp8+fcz7eOq5aowj5y8rKwsAzBd7a86fL1CpVOjTpw82b95ssX3z5s28cF9Do9Hg2LFjiI2NRUJCAmJiYizOm1arxY4dO3z+vFlzbvr06QOlUmmxT2FhIY4cOeLz5+/SpUs4e/as+VrFc2UFlw2ldxN5eXlCVlaWMH/+fCE4OFjIysoSsrKyhMrKSkEQBOG7774TVq5cKWRnZwsnT54UVq1aJYSGhgpPP/20+RhlZWVCdHS0MGnSJCE7O1v45ptvhNDQUGHRokWu+lhOcaNzpdfrhaSkJGHEiBHCwYMHhS1btght2rQRZs6caT6Gr5yra/3666/CkiVLhKysLCE3N1f48ssvhbi4OGHs2LHmfaw5f74iIyNDUCqVwurVq4WcnBxh9uzZQlBQkHDmzBlXV82lnnnmGWH79u1Cbm6usGfPHuGuu+4SQkJCzOflrbfeEsLCwoRvvvlGyM7OFiZNmiTExsYKFRUVLq6581VWVpqvSQDM/97y8vIEQbDu3Dz++ONCmzZthC1btggHDx4UbrvtNqFXr16CXq931cdyiuudq8rKSuGZZ54Rfv31V+H06dPCtm3bhIEDBwqtW7f2yXNlL58PV9OmTRMANHhs27ZNEARB+PHHH4XevXsLwcHBQmBgoJCUlCQsXbpU0Ol0Fsc5fPiwMGTIEEGtVgsxMTHCvHnzvG5qgRudK0EwBbA777xTCAgIEMLDw4WZM2daTLsgCL5xrq6VmZkp9O/fXwgLCxP8/f2Frl27Cq+99ppQXV1tsZ81589X/N///Z/Qvn17QaVSCTfffLOwY8cOV1fJ5SZOnCjExsYKSqVSiIuLEyZMmCAcPXrU/LrRaBRee+01ISYmRlCr1cLQoUOF7OxsF9a4+Wzbtq3R69O0adMEQbDu3NTW1gozZ84UwsPDhYCAAOGuu+4S8vPzXfBpnOt656qmpkZISUkRoqKiBKVSKbRr106YNm1ag/PgK+fKXjJB8IGpsYmIiIiaiVeOuSIiIiJyFYYrIiIiIgkxXBERERFJiOGKiIiISEIMV0REREQSYrgiIiIikhDDFREREZGEGK6IiIiIJMRwRURERCQhhisicsj27dshk8lQVlbmkvJ//vlndOvWDUaj0bxt5cqVaNu2LeRyOZYuXXrd93///fdITk62eD8RkSO4/A0RWe3WW29F7969LQKLVqvF5cuXER0dDZlM1ux16tu3L2bNmoUpU6YAACoqKhAZGYklS5bg3nvvRVhYGAIDA697jJtvvhmpqal46KGHmqPKROTl2HJFRA5RqVSIiYlxSbD69ddfceLECdx///3mbfn5+dDpdLjzzjsRGxt7w2AFAA8//DA++OADZ1aViHwIwxURWWX69OnYsWMHli1bBplMBplMhjNnzjToFlyzZg1atGiB77//Hl27dkVgYCDuu+8+VFdX47PPPkN8fDxatmyJp556CgaDwXx8rVaL5557Dq1bt0ZQUBD69++P7du3X7dOGRkZSElJgb+/v7nsHj16AAA6dOhgruPvv/+O4cOHIyQkBKGhoejTpw8OHDhgPs7YsWOxb98+5ObmSnvSiMgn+bm6AkTkGZYtW4Y///wTSUlJWLBgAQAgKioKZ86cabBvTU0N3n//fWRkZKCyshITJkzAhAkT0KJFC/zwww/Izc3Fvffei8GDB2PixIkATK1HZ86cQUZGBuLi4rBhwwbcfvvtyM7ORufOnRut086dOzFp0iTz84kTJ6Jt27YYOXIk9u3bh7Zt2yIqKgp33XUXkpOTsXz5cigUChw6dAhKpdL8vvbt26NVq1bYtWsXOnToIOFZIyJfxHBFRFYJCwuDSqVCYGAgYmJirruvTqfD8uXL0bFjRwDAfffdh88//xwXLlxAcHAwEhMTMXz4cGzbtg0TJ07EqVOnsH79epw7dw5xcXEAgLlz52Ljxo1IT0/Hm2++2Wg5Z86cMe8PAAEBAYiIiABgCn5iPfPz8/Hss8+iW7duANBoWGvdunWjQZGIyFYMV0QkucDAQHOwAoDo6GjEx8cjODjYYltxcTEA4ODBgxAEAV26dLE4jkajMYelxtTW1pq7BK8nNTUVjz76KD7//HOMHDkS999/v0X9AFMwq6mpserzERFdD8MVEUmufpcbAMhkska3idMfGI1GKBQKZGZmQqFQWOxXP5BdKzIyEqWlpTesz7x58/Dggw/if//7H3788Ue89tpryMjIwD333GPe5/Lly4iKirrhsYiIboThioisplKpLAahSyU5ORkGgwHFxcUYMmSITe/Lycmxat8uXbqgS5cumDNnDiZNmoT09HRzuKqrq8OpU6eQnJxsV/2JiOrj3YJEZLX4+Hjs3bsXZ86cQUlJiWQTb3bp0gWTJ0/G1KlT8c033+D06dPYv38/3n77bfzwww9Nvm/06NHYvXv3dY9dW1uLmTNnYvv27cjLy8Mvv/yC/fv346abbjLvs2fPHqjVagwcOFCSz0NEvo3hioisNnfuXCgUCiQmJiIqKgr5+fmSHTs9PR1Tp07FM888g65du2Ls2LHYu3cv2rZt2+R7HnroIeTk5OD48eNN7qNQKHDp0iVMnToVXbp0wV//+leMGTMG8+fPN++zfv16TJ482ao5sYiIboQztBORR3vuuedQXl6Ojz/+2K73X7x4Ed26dcOBAweQkJAgce2IyBex5YqIPNpLL72E9u3b2z0W7PTp0/joo48YrIhIMmy5IiIiIpIQW66IiIiIJMRwRURERCQhhisiIiIiCTFcEREREUmI4YqIiIhIQgxXRERERBJiuCIiIiKSEMMVERERkYQYroiIiIgk9P8B0RxESOUwrigAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot Average Charge\n",
    "fig, ax = plt.subplots()\n",
    "ax2 = ax.twinx()\n",
    "\n",
    "ax.plot(times / FS2AU, avgCharge, label='avg. charge')\n",
    "ax2.plot(times / FS2AU,\n",
    "         fluence * pulseEnvelope(times) * 1 / AU2FS * 1 / ((1e6 * AU2M)**2),\n",
    "         label='pulse',\n",
    "         c='black',\n",
    "         alpha=0.3)\n",
    "\n",
    "ax.set_xlabel('time (fs)')\n",
    "ax.set_ylabel('avg. charge')\n",
    "ax2.set_ylabel(r'intensity (n. photons / $\\mathrm{\\mu m}^2$ fs)')\n",
    "ax.legend()\n",
    "ax2.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:22.503391Z",
     "start_time": "2022-05-31T09:13:22.161539Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmUAAAHACAYAAADjth/nAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAACazUlEQVR4nOzdeXxU9b3/8deZfZKZ7CSENQiIgLIIBcENWzesxb1ea8EFbdFWRYRW5ScKVbFVUVuva5HN9brgUqyCCu64sAgaZIewZCF7JpPMds7vj5MZErJNkklmknye95ELmTlzzndSSd75fj/n81U0TdMQQgghhBBRZYj2AIQQQgghhIQyIYQQQoiYIKFMCCGEECIGSCgTQgghhIgBEsqEEEIIIWKAhDIhhBBCiBggoUwIIYQQIgZIKBNCCCGEiAESyoQQQgghYoCEMiGEEEKIGCChTAghhBDN+uyzz/jNb35Dr169UBSFt99+u0Wvr66u5tprr+Wkk07CZDJx8cUX1zsmNzeX3/3udwwZMgSDwcDMmTMjMvbOQkKZEEIIIZpVWVnJyJEjefLJJ1v1+kAggN1u59Zbb+Xss89u8BiPx0OPHj2YO3cuI0eObMtwOyVTtAcghBBCiNg3efJkJk+e3OjzXq+X//f//h8vvfQSpaWlnHjiifz9739n0qRJAMTHx/P0008D8OWXX1JaWlrvHFlZWTzxxBMAvPDCCxF/D7FOQpkQQggh2uy6665j3759vPrqq/Tq1YuVK1dy/vnns3XrVgYPHhzt4XUKsnwphBBCiDbZvXs3r7zyCq+//jqnn346AwcOZPbs2Zx22mksWbIk2sPrNGSmTAghhBBtsnHjRjRN4/jjj6/zuMfjITU1NUqj6nwklAkhhBCiTVRVxWg0smHDBoxGY53nHA5HlEbV+UgoE0IIIUSbjB49mkAgQEFBAaeffnq0h9NpSSgTQgghRLNcLhe7du0Kfb537142b95MSkoKxx9/PFdffTXTpk3j0UcfZfTo0RQWFvLJJ59w0kknccEFFwCQnZ2N1+uluLiYiooKNm/eDMCoUaNC5w0+5nK5OHLkCJs3b8ZisTBs2LCOeqtRo2iapkV7EEIIIYSIbevWreOss86q9/g111zD0qVL8fl83H///SxfvpxDhw6RmprKhAkTmD9/PieddBKgt7zYv39/vXPUjiKKotR7vn///uzbty9ybyZGSSgTQgghhIgB0hJDCCGEECIGSCgTQgghhIgB3a7Q3+/3s2nTJjIyMjAYJJMKIYQQnYGqquTn5zN69GhMpq4ZX7rmu2rCpk2bGDduXLSHIYQQQohW+Pbbb/nFL34R7WG0i24XyjIyMgD9f9TMzMwoj0YIIYQQ4cjNzWXcuHGhn+NdUbcLZcEly8zMTPr06RPl0QghhBCiJbpy6VHXfWdCCCGEEJ2IhDIhhBBCiBggoUwIIYQQIgZIKBNCCCGEiAESyoQQQgghYoCEMiGEEEKIGCChTAghhBAiBkgoE0IIIYSIARLKhBBCCCFiQFRD2WeffcZvfvMbevXqhaIovP32282+5tNPP2XMmDHYbDaOO+44nnnmmfYfqBBCCCFEO4tqKKusrGTkyJE8+eSTYR2/d+9eLrjgAk4//XQ2bdrE3Xffza233sqbb77ZziMVQgghhGhfUd37cvLkyUyePDns45955hn69evH448/DsDQoUP5/vvveeSRR7jsssvaaZRCCCGEEO2vU21I/vXXX3PuuefWeey8885j8eLF+Hw+zGZzlEYGAVUjt6yqwecURan/WCPnaeBQlAaObui4Rs8b5jkbv35DxzXy+jDP2egYGnqoBe+1I7/W4Z6zJedt6OVGg9Lo11vUF+3vBUII0VqdKpTl5eWRkZFR57GMjAz8fj+FhYVkZmbWe43H48Hj8YQ+r6ioaJexFVV6OO3va9vl3KJ7MyjgtJlJsJtIsJnpmWAjKy2ewekOxh+XSlZqnIQ2wOv18s0331BaWkpKSgrjxo2TcCaE6FQ6VSiD+jMhmqY1+HjQwoULmT9/fruPC8Bqql+ipzV0YAMPag0fidbgsQ3TGji4oWMbOqeIXaoGZVU+yqp8QBU/HS6v83z/1Diu/EVfrh7Xn8S47htCfvjhB0pLSwEoLi5m69atnHzyydEdlBBCtECnCmU9e/YkLy+vzmMFBQWYTCZSU1MbfM1dd93FrFmzQp8fOnSIYcOGRXxs6U4b2+8Pvz6us2gw6DUS6hoOgOEFxcbO25Kw2pLjGjpvhwbgFrxXj1+lotpHWZWf8iofB0ur2FdYydZDZWzOKWV/kZt/fLCd5z/bw10XDOWKMX263cxZeXk5eXl5KIrCiBEj+OGHHzh06BCDBw/G6XRGe3hCCBGWThXKJkyYwHvvvVfnsdWrVzN27NhGlymsVitWqzX0eXl5eYPHiYY1WKPVop/33SsctJeMBFuDj7u9flZtyeW5z/aws8DFX97YwqacEv520YmYjN2nDeHevXsByMzMpF+/fuTn55OXl0dOTg7Dhw+P8uiEECI8Uf2u7XK52Lx5M5s3bwb0b6ybN28mJycH0Ge5pk2bFjp+xowZ7N+/n1mzZrFt2zZeeOEFFi9ezOzZs6MxfCGiLs5i4oqxffnvbacz57whKAq88u0B7npra4OzeV2RpmmhGfT+/fsD0K9fPwAOHz4ctXEJIURLRTWUff/994wePZrRo0cDMGvWLEaPHs28efMAyM3NDQU0gAEDBvD++++zbt06Ro0axd/+9jf++c9/SjsM0e2ZjAb+dNYgnvrdyRgUeH3DQZZ+tS/aw+oQJSUleL1ezGYzKSkpAKSlpWE0GqmurpbZcSFEpxHV5ctJkyY1+dv80qVL6z125plnsnHjxnYclRCd1+STMrnnwmHMfy+bhe//zOmD0xiU3rVrqo4cOQJAeno6BoP+e6bRaCQtLY38/HwKCgpISEiI5hCFECIs3afoRIhu4tqJWZw1pAfegMr9q7ZFezjtrri4GKDezT5paWl1nhdCiFgnoUyILkZRFOb9Zjhmo8K67Uf4andhtIfUblRVpaSkBCC0dBkU/FxCmRCis5BQJkQXNCAtnv/5hV7s/u/P90Z5NO2nvLycQCCA2WzG4XDUeS4hIQGDwYDP58PlckVphEIIET4JZUJ0UdefNgBFgU9+LmD3ka4ZSoLNYpOTk+u1bzEYDCQlJQGEZtOEECKWSSgToosakBbPL4ekA/DGhoNRHk37CN5ZmZiY2ODzwVAmd2AKIToDCWVCdGGXnNwbgPd+ONwl+5YFw1Zjd1cGH5dQJoToDCSUCdGF/eqEDOItRg6WVLExpzTaw4koTdOaDWXBLZYklAkhOgMJZUJ0YXaLkbOHZQDw0bb8KI8mstxuN4FAAKPRSHx8fIPHOJ1OFEXB6/Xi8Xg6eIRCCNEyEsqE6OImDekBwGc7jkR5JJEVnP0KBq+GGI1G4uLi6hwvhBCxSkKZEF3caYP0UPbT4XIKXV1ntijY5iK4RNkYqSsTQnQWEsqE6OJ6OK0My9SDyRc7u04j2WAoa2zpMijYv6yysrLdxySEEG0hoUyIbuDUQfoWRN/u6zrd7YMh69imsceSUCaE6CwklAnRDYzpnwzApi50B2a4M2XB56WrvxAi1kkoE6IbGN1PD2Xb88qp9PijPJq283q9+Hw+IPxQVl1dTSAQaPexCSFEa0koE6IbyEiw0SvRhqrBloNl0R5OmwVnvex2O0ajscljLRYLZrMZkCVMIURsk1AmRDcRnC3bdKDz7wMZbj1ZUHC2TEKZECKWSSgTops4qY++P2T24c7fGiLcerIgCWVCiM5AQpkQ3cSQnno/r+15FVEeSdu53W4g/FAmd2AKIToDCWVCdBMn1ISyPYWVePydu+C9qqoK0GvKwhHs6i+hTAgRyySUCdFN9EywkWg3E1A1dhd07nASnCkLhq3mBI8LhjkhhIhFEsqE6CYURTm6hJnfeevKAoFAaHPxcENZcEatqqoKTdPabWxCCNEWEsqE6EaCS5g/53beurLgbJfJZAq1umiOzWZDURQ0TQsFOiGEiDUSyoToRgal6wXvewo77/JlS+vJQJ8lDB4fXPoUQohYI6FMiG4kK1W/W3FfJw5lLa0nC6q9hCmEELFIQpkQ3UgwlO0vdqOqnbO2qq2hTGbKhBCxSkKZEN1IryQbZqOC16+SW14d7eG0SmuWL0HuwBRCxD4JZUJ0Iyajgb4pejjprEuYMlMmhOiqJJQJ0c2E6sqKOmcok5kyIUQ47rvvPhRFqfPRs2fPaA+rSaZoD0AI0bE6c7G/qqpUV+vLri0NZVLoL0T3M3z4cD766KPQ50ajMYqjaZ6EMiG6mf6p+ozR/qLOt4wX7DFmMBiwWCwtem0wlAUCAbxeb4tfL4TofEwmU8zPjtUmy5dCdDO9k/RwklvW+Qr9g7NkVqsVRVFa9NraQS54HiFE17Zz50569erFgAED+J//+R/27NkT7SE1SUKZEN1Mr5pQdri08y3jBcOUzWZr1euDr5NQJkTnVVFRQXl5eeijsV06xo8fz/Lly/nwww95/vnnycvLY+LEiRQVFXXwiMMnoUyIbiY4U1ZU6aXaF4jyaFpGQpkQYtiwYSQmJoY+Fi5c2OBxkydP5rLLLuOkk07i7LPPZtWqVQAsW7asI4fbIlJTJkQ3k2A3EW8xUukNcLi0iuN6OKI9pLC1tsg/SEKZEJ1fdnY2vXv3Dn1utVrDel18fDwnnXQSO3fubK+htZnMlAnRzSiKUmsJs3OFk+Cdk22dKZM7MIXovJxOJwkJCaGPcEOZx+Nh27ZtZGZmtvMIW09CmRDdUGetK2vr8mVwhq2xGhQhRNcxe/ZsPv30U/bu3cs333zD5ZdfTnl5Oddcc020h9YoWb4UohsKhrJD3SyUyUyZEN3HwYMHueqqqygsLKRHjx6ccsoprF+/nv79+0d7aI2SUCZEN9Q7SQ8n3W2mTGrKhOg+Xn311WgPocVk+VKIbqhnoj5TlteJNiX3+XwEAvrdom0NZV6vF1VVIzY2IYSIBAllQnRD6U69MLagvPPUVgVnt8xmc6u3SrFYLBgMhjrnE0KIWCGhTIhuKD2hJpRVdJ5g0tZ2GEGyhCmEiFUSyoTohtKdejApcfvw+jvHMl5b68mCJJQJIWKVhDIhuqHkODNmo753ZKGrcyxhSigTQnR1EsqE6IYURaGHI7iEKaFMCCFigYQyIbqpHgl6OCnoJHdgBkNUuN27GxN8vTSQFULEGgllQnRToTswO8lMWTBERWqmTEKZECLWSCgTopvqrKFMZsqEEF2VhDIhuqkeNaHsSCdpixHpUCY1ZUKIWCOhTIhuKtgWozM0kPX7/aFu/pEKZV6vF03T2jw2IYSIFAllQnRTwZmywkpvlEfSvOAsmclkanU3/yCLxYKiKHXOK4QQsUBCmRDdVEq8GYDiytgPJpFaugS9HYjFYqlzXiGEiAUSyoToplLi9YBTUumL8kiaF8lQVvs8EsqEELFEQpkQ3VRKnD5b5PL48fgDUR5N0yLVoyxI2mIIIWKRhDIhuqkEuwmTQa+tivXZMpkpE0J0BxLKhOimFEUhOV6fLSuK8bqySDWODZJQJoSIRRLKhOjGgkuYxTF+B2Z7zZRJrzIhRCyRUCZEN5YS371DmcyUCSFiSdRD2VNPPcWAAQOw2WyMGTOGzz//vMnjX3rpJUaOHElcXByZmZlcd911FBUVddBohehaumsok0J/IUQsimooe+2115g5cyZz585l06ZNnH766UyePJmcnJwGj//iiy+YNm0a06dP56effuL111/nu+++44YbbujgkQvRNQRDWUk3C2UyUyaEiEVRDWWLFi1i+vTp3HDDDQwdOpTHH3+cvn378vTTTzd4/Pr168nKyuLWW29lwIABnHbaafzxj3/k+++/7+CRC9E1HC30j91Q5vP5UFUViHwoq31uIYSItqiFMq/Xy4YNGzj33HPrPH7uuefy1VdfNfiaiRMncvDgQd5//300TSM/P5833niDX//6141ex+PxUF5eHvqoqKiI6PsQojNL7QTLl8HZLLPZjMEQmW9Ztc8ls2VCiHD4fD4OHDjA9u3bKS4ubpdrRC2UFRYWEggEyMjIqPN4RkYGeXl5Db5m4sSJvPTSS1x55ZVYLBZ69uxJUlIS//rXvxq9zsKFC0lMTAx9DBs2LKLvQ4jOLLkThLJIN44NkiVMIURzXC4Xzz77LJMmTSIxMZGsrCyGDRtGjx496N+/PzfeeCPfffddxK4X9UL/4MbAQZqm1XssKDs7m1tvvZV58+axYcMGPvjgA/bu3cuMGTMaPf9dd91FWVlZ6CM7Ozui4xeiM+tMM2XtFcqkLYYQoiGPPfYYWVlZPP/88/zyl7/krbfeYvPmzWzfvp2vv/6ae++9F7/fzznnnMP555/Pzp0723xNUwTG3SppaWkYjcZ6s2IFBQX1Zs+CFi5cyKmnnsqcOXMAGDFiBPHx8Zx++uncf//9ZGZm1nuN1Wqt8828vLw8gu9CiM6tM9x9GenGsUHB7wteb+y+dyFE9Hz11VesXbuWk046qcHnx40bx/XXX88zzzzD4sWL+fTTTxk8eHCbrhm1UGaxWBgzZgxr1qzhkksuCT2+Zs0aLrroogZf43a7MZnqDtloNAL6DJsQomWSa5rHllb5mpyljqb2nimT5UshRENef/31sI6zWq3cfPPNEblm1EIZwKxZs5g6dSpjx45lwoQJPPfcc+Tk5ISWI++66y4OHTrE8uXLAfjNb37DjTfeyNNPP815551Hbm4uM2fOZNy4cfTq1Sti49I0Db/fTyAQ25s0xxKj0YjJZIrJH+qicYl2MwABVaPSG8Bhjeq3hAYFZ7IsFktEzxs8n8yUCSGaU1VVhaZpxMXFAbB//35WrlzJ0KFDOe+88yJ2nah+B77yyispKipiwYIF5ObmcuKJJ/L+++/Tv39/AHJzc+v0LLv22mupqKjgySef5I477iApKYlf/vKX/P3vf4/YmLxeL7m5ubjd7oids7sINvSN9A9P0X5sZgMWowFvQKWsytetQpnMlAkhwnXRRRdx6aWXMmPGDEpLSxk/fjxms5nCwkIWLVrETTfdFJHrRP078M0339zotN/SpUvrPXbLLbdwyy23tMtYVFVl7969GI1GevXqhcVikZmfMGiahtfr5ciRI+zdu5fBgwdHrHWBaF+KopBgN1Po8lDq9tI7yR7tIdXTXsuXMlMmhAjXxo0beeyxxwB44403yMjIYNOmTbz55pvMmzev64SyWOL1elFVlb59+4amKEV47HY7ZrOZ/fv34/V6I16ULdpPot1EoctDWZUv2kNpkMyUCSGize1243Q6AVi9ejWXXnopBoOBU045hf3790fsOjKd0QCZ5Wkd+bp1Tkk1xf7l3SyUyUyZECJcgwYN4u233+bAgQN8+OGHocb3BQUFJCQkROw68lNUiG4uWOwfizNlqqri8+njkpkyIUS0zJs3j9mzZ5OVlcX48eOZMGECoM+ajR49OmLXkeVLIbq5WA5lwUCmKApmszmi5w6GPE3T8Pl8ET+/EKJz27JlCyeeeCIGg4HLL7+c0047jdzcXEaOHBk65le/+lWdtl5tJTNlolElJSVMnTo1tEXV1KlTKS0tjfawRITFciirve9lpG+6MRgMob6HMlsmhDjW6NGjKSwsBOC4447DbDYzevToOqU648aN44QTTojYNSWUiUb97ne/Y/PmzXzwwQd88MEHbN68malTp0Z7WCLCEmpCWak79kJZe9WTBUlXfyFEY5KSkti7dy8A+/btQ1XVdr+mLF92IZWVldx000289dZbOJ1OZs+ezXvvvceoUaN4/PHHW3Subdu28cEHH7B+/XrGjx8PwPPPP8+ECRPYvn07Q4YMaYd3IKIhKYZnyoJhKdLtMIIsFguVlZUyUyaEqOeyyy7jzDPPJDMzE0VRGDt2bGgXoWPt2bMnIteUUNYMTdOo8kWns7/dbGzRks2cOXNYu3YtK1eupGfPntx9991s2LCBUaNGATBjxgxefPHFJs+RnZ1Nv379+Prrr0lMTAwFMoBTTjmFxMREvvrqKwllXUgsL1/KTJkQIlqee+45Lr30Unbt2sWtt97KjTfeGGqL0V4klDWjyhdg2LwPo3Lt7AXnEWcJ738il8vF4sWLWb58Oeeccw4Ay5Yto0+fPqFjFixYwOzZs5s8T3C7qry8PNLT0+s9n56eXm8TedG5BUNZLLbEaO9QFjyvzJQJIRpy/vnnA7BhwwZuu+02CWUiPLt378br9YZu0wVISUmpM6OVnp7eYNBqTEOzdLG6abVovcS42J0pC4YlmSkTQkTTkiVLOuQ6EsqaYTcbyV4Quc1GW3rtcGma1uwxLVm+7NmzJ/n5+fWeP3LkCBkZGWGPS8S+4ExZaQyGMpkpE0J0JxLKmqEoSthLiNE0aNAgzGYz69evp1+/foDe0mLHjh2ceeaZQMuWLydMmEBZWRnffvst48aNA+Cbb76hrKyMiRMntuM7ER0tqdbypapqGAyxMxPa3oX+MlMmhIglsZ82RFgcDgfTp09nzpw5pKamkpGRwdy5c+v0U2nJ8uXQoUM5//zzufHGG3n22WcB+MMf/sCFF14oRf5dTLAlhqqBy+snwRY7TVRlpkwI0Z1IKOtCHn74YVwuF1OmTMHpdHLHHXdQVlbW6vO99NJL3HrrraE9vqZMmcKTTz4ZqeGKGGEzG7GaDHj8KmVuX7cKZTJTJoSIJRLKuhCHw8GKFStYsWJF6LFVq1a1+nwpKSnN1qCJriHBbuZIhYeKan+0h1JHR/QpC15HbmIRQgRVVVVRXFxM79696zz+008/MXz48Ha7rnT0F0LgtOq/n1VUx06xv8/nC3XQbu/ly+D+l0II8cYbb3D88cdzwQUXMGLECL755pvQc+29q42EMiEETpseylye2JkpC86SmUymOrWRkWQwGEIbkcsSphAC4P7772fjxo388MMPvPDCC1x//fW8/PLLQHidDtpCli+7uHXr1kV7CKITcNbUkcXS8mV715MFWSwWfD4fHo8Hh8PRrtcSQsQ+n89Hjx49ABg7diyfffZZqLN/e5c4yEyZEAJHDC5fdlQok2J/IURt6enpbNmyJfR5amoqa9asYdu2bXUebw8SyoQQoeXLihhavmzvbv5B0hZDiO5h4cKFKIrCzJkzmzxuxYoV9dpHWSwWXnnlFT799NN2HKEsXwohAEcwlHXD5cvgTJmEMiG6ru+++47nnnuOESNGNHts7T2ja6uursZsNvOf//wndBNS0JQpUyIyTgllQohQTZkrBkNZe7XDCKrdFkMI0fW4XC6uvvpqnn/+ee6///5WneODDz5g6tSpFBUV1XtOURQCgUBbhwnI8qUQgthsiSE1ZUKIhlRUVFBeXh76aG6W+09/+hO//vWvOfvss1t9zT//+c/89re/JTc3F1VV63xEKpCBhDIhBLHdEqOjasoklAnROQwbNozExMTQx8KFCxs99tVXX2Xjxo1NHhOOgoICZs2aRUZGRpvO0xxZvhRChGrKymNo+bKjC/0llAnROWRnZ9fptN9YicOBAwe47bbbWL16NTabrU3XvPzyy1m3bh0DBw5s03maI6FMNOqBBx5g1apVbN68GYvFQmlpabSHJNpJLNeUSSgTQtTmdDpJSEho9rgNGzZQUFDAmDFjQo8FAgE+++wznnzySTweD0ajMaxrPvnkk1xxxRV8/vnnnHTSSaGm00G33npry95EIySUiUZ5vV6uuOIKJkyYwOLFi6M9HNGOQn3KPLFXUyaF/kKI1vjVr37F1q1b6zx23XXXccIJJ/DXv/417EAG8PLLL/Phhx9it9tZt25dnSayiqJIKBP1VVZWctNNN/HWW2/hdDqZPXs27733HqNGjeLxxx9v8fnmz58PwNKlSyM7UBFzEmKsJYaqqvj9+lg6aqZMVVV8Pl+934CFEJ2T0+nkxBNPrPNYfHw8qamp9R5vzv/7f/+PBQsWcOedd7bbtm8goax5mgY+d3SubY6DFmzpMGfOHNauXcvKlSvp2bMnd999Nxs2bGDUqFEAzJgxgxdffLHJc2RnZ9OvX7+2jFp0QsGaMle1H03T2n0rkeYEZ60URWn3kGQ0GjEajQQCAbxer4QyIUQ9Xq+XK6+8sl0DGUgoa57PDQ/2is617z4MlviwDnW5XCxevJjly5dzzjnnALBs2bI6TfAWLFjA7NmzmzxPr15Req8iqoI1ZX5Vo9qnYreEP63fHjqqnizIYrFQVVWF1+slPj68f3NCiM6ntftBX3PNNbz22mvcfffdkR3QMSSUdRG7d+/G6/UyYcKE0GMpKSkMGTIk9Hl6enq9rSOEAIgzG1EUfWK4wuOLeijrqDsvg2qHMiGEOFYgEOAf//gHH374ISNGjKg3o75o0aKIXEdCWXPMcfqMVbSuHSZN05o9RpYvRWMMBgWH1URFtZ+Kaj/pzuiOp6OK/IOk2F8I0ZStW7cyevRoAH788cc6z0Wy3ENCWXMUJewlxGgaNGgQZrOZ9evXh0JVSUkJO3bs4MwzzwRk+VI0zVkTymKhLUY0li9rX1cIIWpbu3Zth1xHQlkX4XA4mD59OnPmzCE1NZWMjAzmzp1bpyixpcuXOTk5FBcXk5OTQyAQYPPmzYAeAB0OR6Tfgogyp80MZdUxcQdmR4cy2WpJCNEaR44c4aOPPuKqq66KyPkklHUhDz/8MC6XiylTpuB0OrnjjjsoKytr9fnmzZvHsmXLQp8Hp27Xrl3LpEmT2jpcEWNCd2DGQK+yaNSUgYQyIUTDFixY0ODju3fv5q233pJQJupzOBysWLGCFStWhB5btWpVq8+3dOlS6VHWjThjaKslWb4UQsSSlStX1vk8EAhw4MABysvL+dvf/hax60goE0IAtbr6x1Aok0J/IUQs2LRpU73H/H4/M2fOJDs7O2LXkVAmhACO9iqrqI7+8mW0ZsqCy6ZCCNEck8nEzJkzGTFiROTOGbEziZjU2kZ5ovtxWPXeZG5vIMojkeVLIUTnsH//fgYMGBCx80koE0IAEF+zfFnpid7ypcfjYdWqVaxcuRK73U7v3r05+eST2/26wVDm8/liYpspIURs+ec//1nvsby8PJYsWcJvfvObOs+3ZXNyCWVCCADiLdENZZqm8eabb7Jt2zb8fj8VFRW88847eL1eTjnllHa9du0ZOZ/P12EzdEKIzuGxxx5r8HGbzcaaNWtYs2YNoDeSlVAmhGiz0ExZlJYvd+zYwY4dO9A0jTPOOIPi4mJ8Ph8ffPABaWlpDBo0qN2uHdz43Ofz4fF4JJQJIerYu3dvh1ynfbc7F0J0GvE1NWXRmin76quvABg1ahR9+vThlFNOYcyYMQC8/fbbuN3udr2+1JUJIaJNQpkQAqi1fBmFmbKSkhL279+PoiihO5msVivnn38+aWlpuFwu3n///XYdg4QyIUS0SSgTQgAQF8WZsm3btgHQv39/bDYboIcks9nMJZdcgqIo/Pjjj2zfvr3dxiBbLQkhok1CmRACONo81h2FULZjxw4Ahg4dWm+Lpd69ezNhwgQA/vOf/1BdXd0uY5CZMiFEtEkoEw3at28f06dPZ8CAAdjtdgYOHMi9994rP7C6sDhLcO/Ljg1lfr+fgwcPAjBw4MAGu/lPmjSJ5ORkKioq+Oijj9plHBLKhBDHuvvuu/n222877HoSykSDfv75Z1RV5dlnn+Wnn37iscce45lnnuHuu++O9tBEOwnNlHkDaJrWYdfNy8vD7/cTFxdHampqg41jLRYLU6ZMAeD7779n3759ER+HhDIhxLFyc3O58MILyczM5A9/+AOrVq1q150/JJR1IZWVlUybNg2Hw0FmZiaPPvookyZNYubMmS0+1/nnn8+SJUs499xzOe6445gyZQqzZ8/mrbfeivzARUwI1pT5VQ2PX+2w6+bk5ADQt29fFEVptJv/gAEDQo1k3333XXy+yG4HJaFMCHGsJUuWkJ+fz//93/+RlJTEHXfcQVpaGpdeeilLly6lsLAwoteTUNYMTdNw+9xR+WjpbMWcOXNYu3YtK1euZPXq1axbt44NGzaEnp8xYwYOh6PJj+APyIaUlZWRkpLS6q+liG3Buy+hY7daCi5d9u3bF6BeTVlt55xzDg6Hg+LiYj799NOIjkNCmRCiIYqicPrpp/OPf/yDn3/+mW+//ZZTTjmF559/nt69e3PGGWfwyCOPcOjQoTZfS5rHNqPKX8X4l8dH5drf/O4b4sxxYR3rcrlYvHgxy5cv55xzzgFg2bJl9OnTJ3TMggULmD17dpPn6dWrV4OP7969m3/96188+uijYY5edDZGg4LNbKDap1Lp8ZMS3zENVPPz8wHIzMwEmt730m638+tf/5rXXnuNL7/8kuHDh4de11YSyoQQ4Rg6dChDhw7lL3/5C0eOHOHdd9/l3XffBWj2Z2xzJJR1Ebt378br9YbuUgNISUlhyJAhoc/T09NJT09v8bkPHz7M+eefzxVXXMENN9wQkfGK2OSwmqj2ean0dkyxv8/no7i4GCD032ZDhf61DR06lGHDhpGdnc3bb7/N9OnTI9KBP3iO9qwXEUJ0LT169GD69OlMnz49IueTUNYMu8nON7/7JmrXDlc4S50zZszgxRdfbPKY7Oxs+vXrF/r88OHDnHXWWUyYMIHnnnsu7PGIzkm/A9PbYb3Kjhw5gqZpxMXF4XA4CAQCBAL60mlTQeuCCy5g37595Ofn8/bbb3P55ZdjMLStGiN4veAYjEZjm84nhBAtJaGsGYqihL2EGE2DBg3CbDazfv36UKgqKSlhx44dnHnmmUDLly8PHTrEWWedxZgxY1iyZEmbf+iJ2Bfa/9LTMTVlwaXL9PT0OkX+BoMBk6nxb08Oh4Mrr7ySZcuWkZ2dzX//+18mT57cpv9GzWYziqKgaRo+n09CmRCiw0ko6yIcDgfTp09nzpw5pKamkpGRwdy5c+v8kGrJ8uXhw4eZNGkS/fr145FHHuHIkSOh53r27Bnx8YvYEG/p2K7+BQUFAGRkZABN15Mdq3///lx00UWsXLmS7777Dq/Xy4UXXojZbG71eCwWCx6PB6/XG9pZQAghOoqEsi7k4YcfxuVyMWXKFJxOJ3fccQdlZWWtOtfq1avZtWsXu3btqnOzAIS3VCo6p9BMWQfdfRmsJ0tNTQVaFsoARo4cCegblv/www/k5+dz8cUXt/oXh2Aok7oyIUQ0RH096qmnnmLAgAHYbDbGjBnD559/3uTxHo+HuXPn0r9/f6xWKwMHDuSFF17ooNHGNofDwYoVK6isrCQvL485c+a0+lzXXnstmqY1+CG6rvgO3v8yGMqCrVaCYaixIv+GjBw5kquvvpq4uDjy8vJ49tlneeuttzh48GCL/3uVOzCFEA2pqqrC7XaHPt+/fz+PP/44q1evjuh1ojpT9tprrzFz5kyeeuopTj31VJ599lkmT55cr9i8tt/+9rfk5+ezePFiBg0aREFBAX5/x+/VJ0RXFOxV1hF3X2qaRklJCQDJyclAy2fKggYNGsSMGTP48MMP+emnn9iyZQtbtmzB6XSSlZVFeno6KSkp2O12rFZraImz9i8bqqpSXFxMfn4+O3fuxO12EwgEUFW1zX+q6tFmvA0FxcbCY6QeF6Kj2e12Lr744mgPI2IuuugiLr30UmbMmEFpaSnjx4/HbDZTWFjIokWLuOmmmyJynaiGskWLFjF9+vRQm4XHH3+cDz/8kKeffpqFCxfWO/6DDz7g008/Zc+ePaHfrLOysjpyyEJ0aUcL/ds/lFVUVOD3+1EUhaSkJKD1oQwgISGBK664ggkTJvDdd9/x008/UVFRwdatW8M+R35+PqWlpWzevJm0tLQWj0EIoXM4HNEeQkRt3LiRxx57DIA33niDjIwMNm3axJtvvsm8efM6fyjzer1s2LCBO++8s87j5557Ll999VWDr3n33XcZO3Ys//jHP1ixYgXx8fFMmTKFv/3tb9jtDbePOLY+pKKiInJvohNYt25dtIcgOpGjy5ftX1MWnCVLTEwM3enYVDf/cPXp04c+ffpw4YUXcuDAAQ4cOEBRURElJSV4PB6qq6tDs+uKooT+NBqNeL1eVFUlMTGRzMxMDAYDRqMxIn/WFrxuc49F8nEhOlJbbriJRW63G6fTCeg115deeikGg4FTTjmF/fv3R+w6UQtlhYWFBAKB0F1XQRkZGeTl5TX4mj179vDFF19gs9lYuXIlhYWF3HzzzRQXFzdaV7Zw4ULmz58f8fEL0RXFWTpupuzYejJo20zZscxmM8cddxzHHXdc2K/Zu3cvP/74I7169WLMmDFtHoMQomsYNGgQb7/9Npdccgkffvght99+O6DfQZ6QkBCx67Sq0D8/P5+pU6fSq1cvTCYTRqOxzkdLHPtbnaZpjf6mp6oqiqLw0ksvMW7cOC644AIWLVrE0qVLqaqqavA1d911F2VlZaGP7OzsFo1PiO7EUbN82RF7XzYVylpS6B9JUugvhGjIvHnzmD17NllZWYwfPz60e87q1asZPXp0xK7Tqpmya6+9lpycHO655x4yMzNbNV2elpaG0WisNytWUFBQb/YsKDMzk969e5OYmBh6bOjQoWiaxsGDBxk8eHC911it1jrf4MvLy1s8ViG6i7iaPmWuDpgpKy0tBY4W+UNkZ8paQ0KZEKIhl19+Oaeddhq5ubmhVjwAv/rVr7jkkksidp1WhbIvvviCzz//nFGjRrX6whaLhTFjxrBmzZo6b2jNmjVcdNFFDb7m1FNP5fXXX8flcoWKCHfs2IHBYKjXS0sI0XJHZ8raP5QFf0GqPfUvoUwIEat69uxZrwfiuHHjInqNVoWyvn37RuTW61mzZjF16lTGjh0b2lsxJyeHGTNmAPrS46FDh1i+fDkAv/vd7/jb3/7Gddddx/z58yksLGTOnDlcf/31jRb6CyHCF1cTylwdUOh/bCjTNE1CmRAiZn388cd8/PHHFBQU1GlzA0SsX2qrQtnjjz/OnXfeybPPPtumlhRXXnklRUVFLFiwgNzcXE488UTef/99+vfvD0Bubi45OTmh4x0OB2vWrOGWW25h7NixpKam8tvf/pb777+/1WMQQhzlqLn7sr1nyjRNqxfKfD5f6PlohzJVVfH5fF3uDjIhROvMnz+fBQsWMHbs2FaXbYWjVaHsyiuvxO12M3DgQOLi4up94woW8Ibj5ptv5uabb27wuaVLl9Z77IQTTmDNmjUtGq8QIjwddfdlsDErELrNPNgOI7gxeDQEb1YKBAJ4vV4JZUIIAJ555hmWLl3K1KlT2/U6rZ4pE13flClT2Lx5MwUFBSQnJ3P22Wfz97//nV69ekV7aKKdOELNY9t3+TI4SxYfH4/JpF8z2kuXQRaLhaqqKrxeL/Hx8VEdixAiNni9XiZOnNju12lVKLvmmmsiPQ4Rg8466yzuvvtuMjMzOXToELNnz+byyy9vtLmv6PyCd19W+QIEVA2joX1mrJoq8o9WO4yg2qFMCCEAbrjhBl5++WXuueeedr1Oq5vHBgIB3n77bbZt24aiKAwbNowpU6a0uE+ZiJzKykpuuukm3nrrLZxOJ7Nnz+a9995j1KhRrZrdDDbHA+jfvz933nknF198sdTadGHB5UvQg1lw5izSgjtrxNKdl0FS7C+EOFZ1dTXPPfccH330ESNGjKj3M3DRokURuU6rvuPu2rWLCy64gEOHDjFkyBA0TWPHjh307duXVatWMXDgwIgMLhZomobWSGPa9qbY7S2qrZkzZw5r165l5cqV9OzZk7vvvpsNGzaEWpfMmDGDF198sclzNLYZfHFxMS+99BITJ06UQNaF2cwGFAU0TS/2b69QFovtMIIklAkhjrVly5bQz9Iff/yxznORrIFt1XfcW2+9lYEDB7J+/fpQN+6ioiJ+//vfc+utt7Jq1aqIDTDatKoqtp8cne1WhmzcgBIXF9axLpeLxYsXs3z5cs455xwAli1bVqd/24IFC5g9e3aT5zm2Xuyvf/0rTz75JG63m1NOOYX//Oc/LXwXojNRFAW72YjbG6CqHbv6B0NZsMgfjhb6R3v5Mnh9CWVCiKC1a9d2yHVaFco+/fTTOoEMIDU1lYceeohTTz01YoMT4du9ezderze09QPo29cMGTIk9Hl6ejrp6ektOu+cOXOYPn06+/fvZ/78+UybNo3//Oc/sulxFxZn0UNZe261JDNlQojOprS0lMWLF9cp27r++uvr7DLUVq0KZVarNVQTUpvL5Yr6N9RIU+x2hmzcELVrhyucZr6tWb5MS0sjLS2N448/nqFDh9K3b1/Wr19fJ/yJrkWvK/O2a68yCWVCiM7k+++/57zzzsNutzNu3Dg0TWPRokU88MADrF69mpNPPjki12lVKLvwwgv5wx/+wOLFi0NbDHzzzTfMmDGDKVOmRGRgsUJRlLCXEKNp0KBBmM1m1q9fHwpVJSUl7NixgzPPPBNo3fJlbcHgF1xmEl1T8A7M9pwpq6ysBKjTciL435WEMiFErLn99tuZMmUKzz//fKiNj9/v54YbbmDmzJl89tlnEblOq0LZP//5T6655homTJgQKvr2+/1MmTKFJ554IiIDEy3jcDiYPn06c+bMITU1lYyMDObOnYvBYAgd05Lly2+//ZZvv/2W0047jeTkZPbs2cO8efMYOHCgzJJ1cfZ2DmWBQICqmptngnvYQuzNlMkvH0KIoO+//75OIAMwmUz85S9/YezYsRG7TqtCWVJSEu+88w47d+7k559/RtM0hg0bxqBBgyI2MNFyDz/8MC6XiylTpuB0OrnjjjsoKytr1bnsdjtvvfUW9957L5WVlWRmZnL++efz6quvRr0QW7SvUK+ydgplwVkyRVHq7FkbS33KQGbKhBBHJSQkkJOTwwknnFDn8QMHDtS5Yamt2nS/++DBgxk8eHCkxiLayOFwsGLFClasWBF6rLV3wp500kl88sknkRqa6ETsZv3bQnvNlAVDWVxcXGgm1+/3hzb4jZWZMp/Ph6ZpclOLEIIrr7yS6dOn88gjjzBx4kQUReGLL75gzpw5XHXVVRG7TtihbNasWfztb38jPj6eWbNmNXlspJqoCSE63tGasvYp9A+GsoaWLoN7T0ZT7VDo9XqjPnMnhIi+Rx55BEVRmDZtGn6//r3RbDZz00038dBDD0XsOmGHsk2bNuHz+UJ/F0J0TR21fFm7yD9W6slAX1Y1m834fD4JZUJ0Yk8//TRPP/00+/btA2D48OHMmzePyZMnt/hcFouFJ554goULF7J79240TWPQoEHERfhGwLBDWe3GaR3VRE203bp166I9BNHJhAr9fe0TylwuFxC7oQz0ujafz4fH44lovYgQouP06dOHhx56KFTvvmzZMi666CI2bdrE8OHDW3SunJwc+vbtS1xcHCeddFK95xraCac1DM0fUt/111/fYJ+yyspKrr/++jYPSggRPR01U1Z7+TJW2mEESbG/EJ3fb37zGy644AKOP/54jj/+eB544AEcDgfr169v8bkGDBjAkSNH6j1eVFTEgAEDIjFcoJWhbNmyZaFb2murqqpi+fLlbR6UECJ6gpuSt3dNWUMzZbGyVCihTIjYVVFRQXl5eegjnPY1gUCAV199lcrKyla1dWrsph+Xy4XNZmvx+RrTorsvy8vL9Q26NY2Kioo6AwkEArz//vst3sZHCBFb7Ob27VMW6zVlIPtfChHLhg0bVufze++9l/vuu6/BY7du3cqECROorq7G4XCwcuXKeq9vSvDGRkVRuOeee+rUkAUCAb755pvQRuWR0KJQlpSUpHe4VxSOP/74es8risL8+fMjNjghRMeLt7ZvKGuopixWly+lgawQsSc7O5vevXuHPm9qhn3IkCFs3ryZ0tJS3nzzTa655ho+/fTTsINZ8MZGTdPYunVrne9RFouFkSNHNrtTTku0KJStXbsWTdP45S9/yZtvvllnQ3KLxUL//v2b3KZHCBH77B20fNlQSwxZvhRCNMfpdNbZN7cpFoslVOg/duxYvvvuO5544gmeffbZsF4fvLHxuuuu44knngj7uq3VolAW3ENx79699O3bt84WPkKIriHO3H6F/pqmyfKlECJqNE1r1Qz4kiVLKC0t5dFHH2Xbtm0oisKwYcO4/vrrSUxMjNj4WtXRv3///gC43W5ycnLqfeMaMWJE20cmhIiK9tyQvKqqKtS5X5YvhRDt6e6772by5Mn07duXiooKXn31VdatW8cHH3zQ4nN9//33nHfeedjtdsaNG4emaSxatIgHHniA1atXc/LJJ0dkzK0KZUeOHOG6667jv//9b4PPBwLtU4siosPj8TB+/Hh++OEHNm3aFNGiRhF72nND8uAsmdVqrbOxb6wtX8pMmRCdX35+PlOnTiU3N5fExERGjBjBBx98wDnnnNPic91+++1MmTKlzqbkfr+fG264gZkzZ/LZZ59FZMytCmUzZ86kpKSE9evXc9ZZZ7Fy5Ury8/O5//77efTRRyMyMBE7/vKXv9CrVy9++OGHaA9FdIBgS4yqdmge63a7gbqzZKqqhrYtibWZMq/XK/tfCtFJLV68OGLn+v777+sEMgCTycRf/vIXxo4dG7HrtKoo7JNPPuGxxx7jF7/4BQaDgf79+/P73/+ef/zjHyxcuDBigxMtU1lZybRp03A4HGRmZvLoo48yadIkZs6c2epz/ve//2X16tU88sgjkRuoiGntufdlsL+h3W4PPRacjQpubxQLgqFM07RQYBRCdF8JCQnk5OTUe/zAgQMR3fWjVaGssrIy1I8sJSUl1OX2pJNOYuPGjREbXCzQNA2fJxCVD03TWjTWOXPmsHbtWlauXMnq1atZt24dGzZsCD0/Y8YMHA5Hkx+1/6PLz8/nxhtvZMWKFRHf30vEruDyZbVPRVVb9t9gc4IzZbX/e4q1In8Ag8EQ+o1Y6sqEEFdeeSXTp0/ntdde48CBAxw8eJBXX32VG264gauuuipi12nV8uWQIUPYvn07WVlZjBo1imeffZasrCyeeeYZMjMzIza4WOD3qjx326dRufYfnjgTc03PqOa4XC4WL17M8uXLQ+vly5Yto0+fPqFjFixY0Gw/lWBLE03TuPbaa5kxYwZjx44Nbegqur7gTBnoS5jx1lZ9m8BfWk3JmzsxpdhIumgQikFpcKYs1or8gywWC36/X+rKhBA88sgjKIrCtGnTQrPnZrOZm266iYceeihi12l1TVlubi6gd9I977zzeOmll7BYLCxdujRigxPh2717N16vt872ESkpKQwZMiT0eXp6etg7LvzrX/+ivLycu+66K+JjFbHNZjoaytze1oeyslV78ewsxQNYsxKJG53e5ExZrBT5B1mtVtxut4QyIQQWi4UnnniChQsXsnv3bjRNY9CgQRFfRWrVd9urr7469PfRo0ezb98+fv75Z/r160daWlrEBhcLTBYDf3jizKhdO1zhLHXOmDGDF198scljsrOz6devH5988gnr16+v94Ny7NixXH311SxbtizssYnOxWBQsJuNVPkCre5VpvlVqrcXhz53by0kbnR6p5spA7kDUwhxVFxcHCeddFK7nb91vwIfIy4uLmI9OmKNoihhLyFG06BBgzCbzaxfv55+/foBUFJSwo4dO0JNf1uyfPnPf/6T+++/P/T44cOHOe+883jttdcYP358O70LESvirXooq2xlsb/3sAvNq4Y+9+wuRVO1TlNTBtKrTAhR18cff8zHH39MQUFBqN9i0AsvvBCRa4QdyoKbcoZj0aJFrRqMaD2Hw8H06dOZM2cOqampZGRkMHfu3Dq7LrRk+TIY7GqfH2DgwIF16tRE19TWXmX+Aj18WQcm4s2pQPME8BdVNXn3ZSwuX4LMlAkhYP78+SxYsICxY8eSmZnZbm1ywg5lwU05myP9fKLn4YcfxuVyMWXKFJxOJ3fccQdlZWXRHpbohOLMNb3KWhnKfDWhzJwRj+ZV8R6owHe4ssGZslhfvpSZMiHEM888w9KlS5k6dWq7XifsUBbclFPELofDwYoVK1ixYkXosVWrVkXk3FlZWS1u0SE6L3sbe5X58/XwZUqPQ/PXhLJcV5MzZbEaymSmTAjh9XqZOHFiu19HdhQXQtQTbIvR2q7+/uJqAExpNswZ+qyYt8DdqWrKZPlSCBF0ww038PLLL7f7dVpV6H/WWWc1uUz5ySeftHpAQojoa+um5IEyPcgYE61ofn2GtbqoMlQc2xlqymT5UojurXYtvaqqPPfcc3z00UeMGDGi3u4jkaqlb1UoO3ZDap/Px+bNm/nxxx+55pprIjEuESHr1q2L9hBEJ2Sv2f+yNaFMrfaj1bzOmHg0aLlKysGg7xdXexujWJ0pk+VLIbq3Y2vpg9nnxx9/rPN4JGvpWxXKHnvssQYfv++++3C5XG0akBAi+uLMNcuXragpC5TpM0uKzYTBYkRJtoEC1T4PWBueJYPYDWXBDdNrb0QshOj6olFLH9Gast///vcR69UhhIietrTEOLp0qYcaxWTAmGjFo/iAhuvJzGZzzN25bTKZQi1lZLZMCBGkaVq73fgW0VD29ddfY7PZInlKIUQUtKWmLFB+tJ4syJRioxo9lHWGerIgKfYXQgQtXryYE088EZvNhs1m48QTT+Tf//53RK/Rqvn4Sy+9tM7nmqaRm5vL999/zz333BORgQkhoid092Vrasoq9fBljD9aCGtManimLFZ7lAVZLBaqqqqk2F+Ibu6ee+7hscce45ZbbgntMf31119z++23s2/fvjo74LRFq0JZYmJinc8NBgNDhgxhwYIFnHvuuREZmBAieuJqCv1bs81SwK2HL0PtUJZopVppfKYslkMZyEyZEN3d008/zfPPP89VV10VemzKlCmMGDGCW265JbqhbMmSJRG5uBAiNkVipswQd/TbizHRiofGa8piNZTJ8qUQAiAQCDB27Nh6j48ZMwa/v3VNthvSppqy77//nhUrVvDiiy+yYcOGSI1JCBFlbSn0V936N6i6M2WW0PJlZ6opk15lQgjQb2R8+umn6z3+3HPPcfXVV0fsOq2aKTt48CBXXXUVX375JUlJSQCUlpYyceJEXnnlFfr27RuxAYroycrKYv/+/XUe++tf/8pDDz0UpRGJjhJcvnS3oqN/4zNlelirHco6Q00ZyEyZEEIv9F+9ejWnnHIKAOvXr+fAgQNMmzatTqPZtjSSbVUou/766/H5fGzbto0hQ4YAsH37dq6//nqmT5/O6tWrWz0gEVsWLFjAjTfeGPrc4XBEcTSioxxdvmz5tLwarCmLq1tT5q2ZKbNZjt6hHeszZbJ8KYQAvWHsySefDMDu3bsB6NGjBz169KjTTLatrX1aFco+//xzvvrqq1AgAxgyZAj/+te/OPXUU9s0INF6lZWV3HTTTbz11ls4nU5mz57Ne++9x6hRo3j88cdbdU6n00nPnj0jO1AR89q2fFn/7ktDnAmPogc8S8AYejzWa8pkpkwIAR3XSLZVNWX9+vXD5/PVe9zv99O7d+82DyqWaJqGr7o6Kh8tbU43Z84c1q5dy8qVK1m9ejXr1q2rU+s3Y8YMHA5Hkx85OTl1zvn3v/+d1NRURo0axQMPPCA/nLqJ1hb6a6p2tKas1kyZoih4Dfq5zL6j33Y6y/Kl1JQJITpCq2bK/vGPf3DLLbfwv//7v4wZMwZFUfj++++57bbbeOSRRyI9xqjyezz885rLo3LtW5e9gTnMZrwul4vFixezfPlyzjnnHACWLVtGnz59QscsWLCA2bNnN3meXr16hf5+2223cfLJJ5OcnMy3337LXXfdxd69eyPeLE/Enjhz6/a+1Kr9UPO7RO2aMgCvpv8iZ/Ycnd6P9ZkyWb4UQnSkVoWya6+9Frfbzfjx40P7wQX3hrv++uu5/vrrQ8cWFxdHZqSiSbt378br9Yaa2gGkpKTUWWJOT08nPT097HPefvvtob+PGDGC5ORkLr/88tDsmei6gsuXVb4AqqphMIRXJxGomSVTrEYU09EZMZ/PRwAVAHO1EnosOBscqzVlwbDo9/tRVTW07ZIQQrSHVoWy1tYndUYmq5Vbl70RtWuHK5ylzhkzZvDiiy82eUx2djb9+vVr8LngHSe7du2SUNbFBZcvAar9gdDdmM0J3XlZq54MoLq6GgBFA2OV/t9qcPap9h6TsSa4J6emaXi9XtlGTgjRrloVyq655ppIjyNmKYoS9hJiNA0aNAiz2cz69etDoaqkpIQdO3Zw5plnAi1fvjzWpk2bAMjMzIzQqEWsspuPhjK3twWhrKqmnsxe9/hgKLNgQnPpx8R6PVmQxWLB4/Hg8XgklAkh2lWrQhno3W3ffvtttm3bhqIoDBs2jClTpmA0Gpt/sYg4h8PB9OnTmTNnDqmpqWRkZDB37tw6MxAtWb78+uuvWb9+PWeddRaJiYl899133H777UyZMqXRmTTRdRgMCjazgWqf2qJif626JpRZ634fCIUyzUTApc+QBUNZrC5dBlmtVjwej9SVCSHaXatC2a5du7jgggs4dOgQQ4YMQdM0duzYQd++fVm1ahUDBw6M9DhFGB5++GFcLhdTpkzB6XRyxx13UFZW1qpzWa1WXnvtNebPn4/H46F///7ceOON/OUvf4nwqEWsireYqPZ5W1Tsr3r0Y5XGQhlm1IrOFcqkLYYQojFnn302e/bsYc+ePRE5X6tC2a233srAgQNZv349KSkpABQVFfH73/+eW2+9lVWrVkVkcKJlHA4HK1asYMWKFaHHWvu/xcknn8z69esjNTTRCdktRqgEdwsayGo1oezYmbKqqioArJqJQIVed9bZQpm0xRBCHOviiy+mqKgoYudrVSj79NNP6wQygNTUVB566CFpHitEFxHXigayzc+UmVDdPrSAFvPd/IOC45NQJoQ41p///OeInq9VocxqtVJRUVHvcZfLFfNFu0KI8NgtLe9VpoVCWeOF/mj6XZqdpdBfQpkQAuDjjz/m448/pqCgAFVV6zz3wgsvROQarQplF154IX/4wx9YvHgx48aNA+Cbb75hxowZTJkyJSIDE5Gxbt26aA9BdFJx5uBMWQuWL70NL18GQ5nNbAUfBCq8nWb5UkKZEGL+/PksWLCAsWPHkpmZ2eY9LhvTqlD2z3/+k2uuuYYJEyZgNuv9iHw+HxdddBFPPPFEi8711FNP8fDDD5Obm8vw4cN5/PHHOf3005t93ZdffsmZZ57JiSeeyObNm1vzNoQQTWjNVkvNLV/aLDZwQ8AloUwI0Xk888wzLF26lKlTp7brdVoVypKSknjnnXfYtWsX2dnZAAwbNoxBgwa16DyvvfYaM2fO5KmnnuLUU0/l2WefZfLkyU02MAUoKytj2rRp/OpXvyI/P781b0EI0YzWbEreWKF/MJRZ7TYoBbXCKzVlQohOw+v1MnHixHa/TqvbaC9evJiLL76YK664giuuuIKLL764xXsiLlq0iOnTp3PDDTcwdOhQHn/8cfr27cvTTz/d5Ov++Mc/8rvf/a7OlkJCiMiKq7XVUriamymzx9kB8JV58Pn0uzA7SyiTlhhCdF833HADL7/8crtfp1UzZffccw+PPfYYt9xySygYff3119x+++3s27eP+++/v9lzeL1eNmzYwJ133lnn8XPPPZevvvqq0dctWbKE3bt38+KLL4Z1HSFE68SFCv1bXlPWaCiLjwOgqsQFVjAYDKESiFgVDGWqquLz+WJ+vEKIyKuurua5557jo48+YsSIEfW+DyxatCgi12lVKHv66ad5/vnnueqqq0KPTZkyhREjRnDLLbeEFZYKCwsJBAJkZGTUeTwjI4O8vLwGX7Nz507uvPNOPv/889BG6M0Jbo8S1NBdo0KI+iK5fBnsU2ZPiAPceijrGft3XoIeHE0mE36/H4/HI6FMiG5oy5YtjBo1CoAff/yxznORLPpvVSgLBAKMHTu23uNjxozB7w//t2qo/2Y0TWvwDQYCAX73u98xf/58jj/++LDPv3DhQubPn9+iMQkhjt592bJCf/3fv2JpeKYsLjEecOMudUPP2F+6DLJaraFQ5nA4oj0cIUQHW7t2bYdcp1U1Zb///e8brPt67rnnuPrqq8M6R1paGkajsd6sWEFBQb3ZM9BnuL7//nv+/Oc/YzKZMJlMLFiwgB9++AGTycQnn3zS4HXuuusuysrKQh/BGxOEEE2L1EyZpmlHly+T9EDjKa8EOlcoAyn2F0K0r1ZvSL548WJWr17NKaecAsD69es5cOAA06ZNY9asWaHjGltntVgsjBkzhjVr1nDJJZeEHl+zZg0XXXRRveMTEhLYunVrnceeeuopPvnkE9544w0GDBjQ4HWsVmudb/zl5eXhv0nBqlWrWLBgAVu2bCE+Pp4zzjiDt956K9rDEh0groXNYzVVQ/PqDRVr15R5vV40TQMgPsWJB6gqcwP2mAxlmqbi8Rbg9RzB6y3E6yuiuHgjuXn5bPt5HWXlVgKBKgIBN2qgGk3zo2mB0J+q5j/mMRXQgmcPXgTt2MdqHtcf0eq9Jvg1FKIzsFrSGD/+/WgPo01ycnKa7ARxrEOHDtG7d+82XbNVoezHH3/k5JNPBmD37t0A9OjRgx49etRZa21unXXWrFlMnTqVsWPHMmHCBJ577jlycnKYMWMGoM9yHTp0iOXLl2MwGDjxxBPrvD49PR2bzVbvcREZb775JjfeeCMPPvggv/zlL9E0rV4wFl3X0bsvwytJ0GqFt9ozZcFZMoPBgDVJv/vSU1mN5rdGvaYsEPBQUbGVsrINlJdvxV21F7d7H6paXee4vHwvBfnBIv/Yr4MTItoUxdj8QTHuF7/4BVOmTOHGG28MNco/VllZGf/3f//HE088wR//+EduueWWNl2zVaEsUmurV155JUVFRSxYsIDc3FxOPPFE3n//ffr37w9Abm4uOTk5EblWd1BZWclNN93EW2+9hdPpZPbs2bz33nuMGjWKxx9/vEXn8vv93HbbbTz88MNMnz499PiQIUMiPGoRq1q6fBlcusQAmI5WRoQax9psGOLMYFTwBryo1f6ozJSpqpcjR9aQX/A+RUWfoqpV9Y5RFCMWcxoWSxpmSwrV1Qb8/irSUjM4bsAgjMY4jEY7BqMdg2JCUUwoirHWn0YUxYxiMKJgABSo+SVVofYvq0rwgrWeU+o9V/c1QsQ+RWn1QlzM2LZtGw8++CDnn38+ZrOZsWPH0qtXL2w2GyUlJWRnZ/PTTz8xduxYHn74YSZPntzma0b9q3bzzTdz8803N/jc0qVLm3ztfffdx3333Rf5QdWiaRqaT23+wHagmA0tuqtjzpw5rF27lpUrV9KzZ0/uvvtuNmzYELpjZMaMGbz44otNniPYuHfjxo0cOnQIg8HA6NGjycvLY9SoUTzyyCMMHz68LW9LdBLxNcuX4Rb6h3qUWUx1/rsN1ZPZ7SiKgtFhxhvwoVYFOjSUBQIeDh5cSs6BpXi9BaHHLZY0EhPHkJg4mvi4QcTFDcBm64PBcPTbY0ryftTAFnr27MmAAb/osDELIaInJSWFRx55hPvvv5/333+fzz//nH379lFVVUVaWhpXX3015513XkRX66IeymKd5lM5PK/xvmntqdeCifXuYmuMy+Vi8eLFLF++nHPOOQeAZcuW0adPn9AxCxYsYPbs2U1fs1cvAPbs2QPowXfRokVkZWXx6KOPcuaZZ7Jjxw5SUlJa85ZEJ9LambJG97202fTnHRa8AT9Kta/DQllR0Wf8vP0eqqsPAmCx9CAz8zLS0yfjdAxv9pcfKfQXovuy2WxceumlXHrppe1+LQllXcTu3bvxer11djlISUmps9yYnp5Oenp6WOdTVX12cO7cuVx22WWA3ri3T58+vP766/zxj3+M4OhFLIprYShrrJt/sEdZMJTpM2VeLNXtP1OmaQF27X6YnJznAbBae3LcgNvp2XMKBkP4tWESyoQQHUFCWTMUs4FeC9p/v6vGrh2ucO7MasnyZWZmJqDvaRpktVo57rjjpM6vmzi6IXmYhf5hzpQp8WZ8AT8mt79dC/1V1cNPP91BwZH/AtCn91QGDfoLRmNci88loUwI0REklDVDUZSwlxCjadCgQZjNZtavXx+6hbekpIQdO3Zw5plnAi1bvhwzZgxWq5Xt27dz2mmnAeDz+di3b1/oRgzRtYWWL32BRps616Y2s8VSMJQFbHo7iPYs9Ne0AFt/vJXCwo9QFAvDhz1MRsaFrT5fcJyBQAC/3x/2jiJCCNES8p2li3A4HEyfPp05c+aQmppKRkYGc+fOxWA4OtvWkuXLhIQEZsyYwb333kvfvn3p378/Dz/8MABXXHFFu7wHEVuCfco0DTx+FZu56V9OtGA3/2ZCmd+mhzuTX4no9iS17dj5AIWFH2EwWBg54t+kpJzapvMZjUaMRiOBQACPxyOhTIhOYOHChbz11lv8/PPP2O12Jk6cyN///veY7iIg31m6kIcffhiXy8WUKVNwOp3ccccdlJWVtel8JpOJqVOnUlVVxfjx4/nkk09ITk6O4KhFrLLXCmGVHn8YoSy85Uu/Ra9XNHlbtaFIs/Ly3uXgwWUADBv2aJsDWZDVasXtduPxeIiPj4/IOYUQ7efTTz/lT3/6E7/4xS/w+/3MnTuXc889l+zs7Jj9NyyhrAtxOBysWLGCFStWhB5btWpVq89nNpt55JFHeOSRRyIxPNHJGA0KVpMBj1/F7Q2Q2szxjRX6HxvKfGa9/tHsi/wsWXX1YbbvmAfAgKxbyEi/IGLnrh3KhBCx74MPPqjz+ZIlS0hPT2fDhg2cccYZLTrXtddey/XXX9/i17VU+/yqKoToEo529W/+DszmZsrsdr2bv9+sz5QZvZEPZdt3zMfvryAhYRRZWX+O6Lml2F+Izi24ctSalk4VFRWce+65DB48mAcffJBDhw5FeniAhDIhRBNasv9luDNlXoN+nMVvRPNHrjFzUdFnNYX9JoYN/Xud5q+RIKFMiNhQUVFBeXl56COcf5OapjFr1ixOO+20VjV7ffPNNzl06BB//vOfef3118nKymLy5Mm88cYb+Hy+1ryNBkko6+LWrVvX4i2WhAg62kC2+bYYwb0vDZam+5T5lAAoChaTmYArMt/MNC3Azl0PAtCnzzTi4wdF5Ly1SSgTIjYMGzaMxMTE0MfChQubfc2f//xntmzZwiuvvNLq66ampnLbbbexadMmvv32WwYNGsTUqVPp1asXt99+Ozt37mz1uYOkpkwI0aijvcrCX75UrHW/rdSbKfN5MdiM2IwWVJcXktreFqPgyIdUVu7EZEpgQFbbNgRujIQyIWJDdnY2vXv3Dn3eXGudW265hXfffZfPPvuszi43rZWbm8vq1atZvXo1RqORCy64gJ9++olhw4bxj3/8g9tvv73V55ZQJoRoVPAOzJYtXx6dgFdVNRRigqGsuroag92ExWSJyEyZpmns2/e/APTtex1mc0Kbz9mQ4Dd+r9fbLucXQoTH6XSSkND8v3NN07jllltYuXIl69atY8CAAa2+ps/n491332XJkiWsXr2aESNGcPvtt3P11VfjdDoBePXVV7npppsklAkh2kdrZsoMtWbKas8qHRvKrAYLakXbA05x8We4XD9jNDro22dam8/XmGAoC878CSFi25/+9Cdefvll3nnnHZxOJ3l5eQAkJiaGbjwKV2ZmJqqqctVVV/Htt98yatSoesecd955JCUltWnMEsqEEI06WujffE1ZQ4X+wQBjMpkwmUyoqorP50OxmbCokakpO3joZQB6ZV6O2ZzU5vM1RpYvhehcnn76aQAmTZpU5/ElS5Zw7bXXtuhct912G3fccQdxcXW3adM0jQMHDtCvXz+Sk5PZu3dvW4Yshf5CiMbV3mqpOQ21xDi2HUYw0JjizFiM5jbPlFVX51JY+AkAvXv/rk3nak5om6hAIKJ3Wwkh2oemaQ1+tDSQAdx33324XK56jxcXF7dpWfRYEsqEEI2Kb8HyZWimzFI/lNVeugSwJejdtAOutoWyw4f/D1BJShpPfPzANp2rOUajEbPZDMhsmRDdjaZpDT7ucrlC398iQZYvhRCNsofZp0wLaFDTc6z28uWx7TCCYSYuKR7yQW3D8qWmaeTmvQVA717/0+rztITVasXn81FdXY3D4eiQawohomfWrFkAKIrCvHnz6ixfBgIBvvnmmwbry1pLQplo0Lp16zjrrLMafO7bb7/lF7/4RQePSERDnCW8uy+Dm5FDw8uX9WbKEts+U1Ze/gPV1QcxGuPo0eOcVp+nJWw2Gy6XS4r9hegmNm3aBOi/BG7duhWLxRJ6zmKxMHLkSGbPnh2x60koEw2aOHEiubm5dR675557+Oijjxg7dmyURiU62tG7L5su9FeDoc2ooJiOVkUcG8qCM2X25HjA26aZsvyC/wCQlvYrjMaW3UnVWse+DyFE17Z27VoArrvuOp544omwWnG0hdSUdSGVlZVMmzYNh8NBZmYmjz76KJMmTWLmzJktPpfFYqFnz56hj9TUVN59912uv/56FCXyexaK2GQPe6as6X0vj50pi0vWl/5Ut79VWy1pmkpBwX8ByEj/dYtf31rHvg8hRPewZMmSdg9kIDNlzdI0LWp3WpnN5hYFoDlz5rB27VpWrlxJz549ufvuu9mwYUNovXvGjBm8+OKLTZ4jOzubfv361Xv83XffpbCwsFV3rYjOK9wNycPd9/JoTZkDDMWgQqDShymxZV39y8o24vHkYTI5SU09o0WvbQvpVSZE9zFr1iz+9re/ER8fH6ota8yiRYsick0JZc3w+Xw8+OCDUbn23XffXWf9uikul4vFixezfPlyzjlHr69ZtmxZnS0lFixY0Ozad69evRp8fPHixZx33nn07ds3zNGLrsBuDrPQv4UzZTa7jUC83jxWrfBCC0NZYZG+pJCaOgmDoe3bNIVLZsqE6D42bdoUmpQJ1pY1JJKrRxLKuojdu3fj9XqZMGFC6LGUlBSGDBkS+jw9PZ309PQWn/vgwYN8+OGH/N///V9Exio6j/AL/Zve9/LYPmVWqxWPQ+9T1poGskVF6wBIS234ZpT2IjVlQnQfwXqyY//eniSUNcNsNnP33XdH7drhaqyHSm2tXb5csmQJqampTJkyJezxiK4h7EL/MJYvNU2rsw+mz2mB3Ep9U/IWqK7OxeX6GVA6dOkSZPlSiO6qqqoKTdNCLTH279/PypUrGTZsGOeee27EriOhrBmKooS9hBhNgwYNwmw2s379+lCoKikpYceOHZx55plA65YvNU1jyZIlTJs2rUUhUXQNbS30r92nzOv1hn55sFqtuB36f08tnSkrKvoUgMSEUZjNyS16bVsd29Vf/k0I0T1cdNFFXHrppcyYMYPS0lLGjRuHxWKhsLCQRYsWcdNNN0XkOhLKugiHw8H06dOZM2cOqampZGRkMHfuXAyGozfYtmb58pNPPmHv3r1Mnz490kMWnUBcmM1jG+rmD3VnymovXSqKgsGp/7LT0q2WgkuXqamTWvS6SAh29ff5fHg8HgllQnQTGzdu5LHHHgPgjTfeoGfPnmzatIk333yTefPmSSgT9T388MO4XC6mTJmC0+nkjjvuoKysrE3nXLx4MRMnTmTo0KERGqXoTI7WlPnRNK3RgtZwCv2Dfw8uARpbMVOmaQFKStcDdPjSZZB09Rei+3G73TidTgBWr17NpZdeisFg4JRTTmH//v0Ru470KetCHA4HK1asoLKykry8PObMmdPmc7788st8+eWXERid6IyCy5eqBp4m+ompNR39a9eU1d64u+FQ1vKZMpfrZ/z+CoxGBw7HsBa8k8iRYn8hup9Bgwbx9ttvc+DAAT788MNQHVlBQUFE+5dJKBNCNCrOXGsfyyaWMDVv/X0vaxfDW63WendiGloxU1ZS8g0ASUljMRiiM9Evxf5CdD/z5s1j9uzZZGVlMX78+FCng9WrVzN69OiIXUeWL4UQjTIZDViMBrwBFbcvQGNl9cG9Lxva99JisWA0Guv1LDMGa8pacPdlSakeypKTxrfofUSS9CoTovu5/PLLOe2008jNzWXkyJGhx3/1q19xySWXROw6Esq6uHXr1kV7CKKTi7Ma8brVJttiNFTof+zMWGMzZarbjxZQUYxNT9xrWoDS0m8BSE4+pTVvJSIklAnRPQW3Haxt3LhxEb2GhDIhRJPizEZK8TV5B2ZDhf6122E09LkhzqwXUKigunwYm+nq73Jtx+8vj2o9GUhNmRDd1ccff8zHH39MQUEBqlq3xvaFF16IyDUklAkhmhROr7KGmsc2usVSzeeKQcEQb0at8BEII5QFly6TksZErZ4MZKZMiO5o/vz5LFiwgLFjx5KZmRnRrZVqk1AmhGhSsFdZk4X+zYQyVVXxevXaseDyJeh3YOqhrPm6svKyzQAkJo5p2RuIMCn0F6L7eeaZZ1i6dClTp05t1+vI3ZdCiCaFM1OmeesvX9YOZcGly2Dz1aBQXVlF83dglpX/AOid/KMpGCprt/wQQnRtXq+XiRMntvt1JJQJIZpUu4FsY9QGNiRvqHFscOkvKHgHZnMzZV5vIdXVBwCFhIQRLXsDEWYwGEJbrwXDphCia7vhhht4+eWX2/06snwphGhSaFNyX8MzZZpfhYC+p2VDM2V2u73RUGYINpBtpldZcJYsPn4QJpOzpW8h4ux2O16vl6qqqog2jhRCxKbq6mqee+45PvroI0aMGFFvi7VFixZF5DoSyoQQTbKbm97/MjhLBg23xKi9fFm7ngxqbbXUTFf/8rJNACQkjGzyuI5it9spKyuTmTIhuoktW7YwatQoAH788cc6z0Wy6F9CmWjUjh07mDNnDl9++SVer5eTTjqJ+++/n7POOivaQxMdKK6ZmrJQkb/ZgGI8+s2pdguMRmfKwmwgGyv1ZEHH9l4TQnRta9eu7ZDrSE2ZaNSvf/1r/H4/n3zyCRs2bGDUqFFceOGF5OXlRXtoogOFli8bqSlrqB0GNFxT1uhMWRPLl5oWoLx8CwAJMRLKju29JoTo+j7//HN+//vfM3HiRA4dOgTAihUr+OKLLyJ2DQllXUhlZSXTpk3D4XCQmZnJo48+yqRJk5g5c2aLz1VYWMiuXbu48847GTFiBIMHD+ahhx7C7Xbz008/RX7wIjaV5zLxyGsMUHKbmCmrvxk5hFfof7SmrPGZssrK3QQCLozGOOLjB7fufURYMFxKKBOie3jzzTc577zzsNvtbNy4MdQ8uqKiggcffDBi15FQ1gxN0wgE3FH50DStRWOdM2cOa9euZeXKlaxevZp169axYcOG0PMzZszA4XA0+ZGTkwNAamoqQ4cOZfny5VRWVuL3+3n22WfJyMhgzJjo9okSHcRXDS+cx5l7H+NNy71oVWUNHhbq5m9pPJQd280/yOisaYlR6ddvGGhARYVev+FwDItq09jaZPlSiO7l/vvv55lnnuH555+vU+Q/ceJENm7cGLHrxMZ3uBimqlWs+/SkqFx70plbMRrjwjrW5XKxePFili9fzjnnnAPAsmXL6NOnT+iYBQsWMHv27CbP06tXL0AvXFyzZg0XXXQRTqcTg8FARkYGH3zwAUlJSa17Q6Jz2fYulO4HIEVxMbLov8Bp9Q5raPlS07RQYLFaraHfKo9dvjTEm8GoQEAjUO7FlFI3tAFUuLYB4HRGb2ulY8nypRDdy/bt2znjjDPqPZ6QkEBpaWnEriOhrIvYvXs3Xq+XCRMmhB5LSUlhyJAhoc/T09NJT08P63yapnHzzTeTnp7O559/jt1u59///jcXXngh3333HZmZmRF/DyLGbH0dAK85AYuvnGGVXzd4WEP7Xvr9fgKBmscNBjRNQ1GUUH+vIEVRMCZaCRRXEyj3NBzKKvTlcqdjeNvfU4QEQ5mqqng8nlCXfyFE15SZmcmuXbvIysqq8/gXX3zBcccdF7HrSChrhsFgZ9KZW6N27XCFs9Q5Y8YMXnzxxSaPyc7Opl+/fnzyySf85z//oaSkJNSH6amnnmLNmjUsW7aMO++8M+yxiU5IDUDOegC2jbybkd/fyfHVP0LAD8a63zYamikLziApihIKZzabrcFbx40JFj2UldWvK9M0DVcMzpQZDIZQrVx1dbWEMiG6uD/+8Y/cdtttvPDCCyiKwuHDh/n666+ZPXs28+bNi9h1JJQ1Q1GUsJcQo2nQoEGYzWbWr19Pv379ACgpKWHHjh2ceeaZQMuWL91uN6D/8KnNYDCgqg3X/oguJG8reMrBmkDpwClUfncv8YoHindDjyF1Dg21xGikR1ljd14GBTciD5R56j1XXX0Qv78cRTETHz+o7e8rgoLvraqqisTExGgPRwjRjv7yl79QVlbGWWedRXV1NWeccQZWq5XZs2fz5z//OWLXkVDWRTgcDqZPn86cOXNITU0lIyODuXPn1glVLVm+nDBhAsnJyVxzzTXMmzcPu93O888/z969e/n1r3/dXm9DxIoD3+p/9h1PnN3Odq0vJyu79LB2TChTm9n3Mhjw4+Ia/uXGmFiz1VJ5/ZmyClc2APHxgzEYLPWejya73U5paanUlQnRTTzwwAPMnTuX7OxsVFVl2LBhOByOiF5D7r7sQh5++GHOOOMMpkyZwtlnn81pp53W6jsl09LS+OCDD3C5XPzyl79k7NixfPHFF7zzzjuMHBkbXdVFO8qv6VjdaxR2s5Ftav+6j9eiNbB8GU43/yBjQuMzZRUVeiiLpaXLILkDU4juIycnB03TiIuLY+zYsYwbNy4UyIJdCyJBZsq6EIfDwYoVK1ixYkXosVWrVrX6fGPHjuXDDz+MxNBEZ1Og13GRPpQ4i5Fdmr6sTfEeil98iaJnnyX1D38gZervQ6Fs748b6NG/it5DhrZwpqzxUOYKhjJH7IUyuQNTiO5jwIAB5Obm1lttKioqYsCAAaHa2baSmTIhRF2aViuUDSPOYiJH078RBfL2kP/QQ/iPHCH/gQfwFRSg1jSP3Ze9iTcfnIe7vKxOHVmzM2VhLF86nbFz52WQNJAVovsI3kF+LJfLVa//YlvITJkQoq6yA+CtAIMZUgdh98KOeDf3JaTw259yUfxHa7ty3vuSqj2JJGHGr3rxVVex/evPqbI5gTCXL4MzZeVeNFVDMejf+LzeIjwefUsvh2NIg6+NJgllQnR9s2bNAvSb/u655546M/6BQIBvvvkmtFF5JEgo6+LWrVsX7SGIziY4S5Z2PBjN5Lh+oqT3+7ypOMgo8zOp5jANhU83WBlrqwKLGZ+mz3Tt3fgd8SdPBMBkMuH36zNpjYYyhwUUQNVQXT6MCXroCzaNtdv7YzI52+WttkXtmjJVVevdqSyE6Pw2bdoE6DNlW7durdNr0WKxMHLkyGa7GrSEhDIhRF2hpcsTAHj55xdRFL0PXu+avegTLriA/V/tpEq1YVL0mrGkPsdzeMcucnftoM+wk+uc0mq1YjTW3YYpSDEqGJwW1HIvgXJPKJRVurYD4HAMjejbixSbzRZqEVNdXd1ozZwQovNau3YtANdddx1PPPFEqG9ne5Ff7RrQ0j0nhU6+bl1EyV79z9RBeAIe1h7QvymNqlDoX6A/lXjRFEqS9M3BTYr+bcRg74fRZKLaVUFFmb5PZrAGo7FZsqCGiv1dlTsAcMQfH4E31T6CQSx4M4MQomtasmRJuwcykJmyOoKbjLrd7mZ/iIj6gj+Yam/WKjqh4ppQljyALUe2UOmrRAkk8KsDCpZAEV4T2E85hQrnD2hqFSYluNzoIK3fAPL37KSitAQ4GtSb+/dkSrDggzpd/SsrdwIQ74jtUOZyuSSUCdENfPzxx3z88ccUFBTUa6L+wgsvROQaEspqMRqNJCUlUVCgTwfExcU1eLeFqEvTNNxuNwUFBSQlJTW6TCU6iZJ9+p/JWWwu2AyALTCI+KICoIi8JEiszqUiMQstUITJoPcwKytVyRzal/w9O6k6JqS0dKZM09RQKJOZMiFEtM2fP58FCxYwduxYMjMz2y0bRD2UPfXUUzz88MPk5uYyfPhwHn/8cU4//fQGj33rrbd4+umn2bx5Mx6Ph+HDh3Pfffdx3nnnRWw8PXv2BAgFMxG+pKSk0NdPdFIBH5Qd1P+enMUPu14CwMEgAhXFAOSlKFQd2IrX5MDgOYChZvnSr4LN2QMAj1ef8Qr+NtlcvVWwLYa/JpRVVx8iEHCjKBbs9v4RfIORJaFMiO7hmWeeYenSpUydOrVdrxPVUPbaa68xc+ZMnnrqKU499VSeffZZJk+eHNoU+1ifffYZ55xzDg8++CBJSUksWbKE3/zmN3zzzTeMHj06ImNSFIXMzEzS09Px+XwROWd3YDabZYasKyg7AFoATDa0+HR+OPIDAMnGwRhc+gbleclQdjAPIwOx+IsAfbbUDxhMKQD4AiooCqqqoihK8zNlyXqfn0CJHspCS5fxx2EwRP13x0ZJKBOie/B6vUycOLHdrxPV73aLFi1i+vTp3HDDDQA8/vjjfPjhhzz99NMsXLiw3vGPP/54nc8ffPBB3nnnHd57772IhbIgo9EoIUN0P7WWLnOr8in1lGIymEg3HofNpc9+5SYrOHPL6AmY1QoAAoo+I6YoSWhAsNrC7/djNpubnSkz1YQyf7HedNbl0ov842N46RIklAnRXdxwww28/PLL3HPPPe16naiFMq/Xy4YNG7jzzjvrPH7uuefy1VdfhXUOVVWpqKggJSWlPYYoRPdTK5TtK9f/3s/ZD6ffSrxbb5J6JBEMRXpAM6LPbKk193H7fU4wGkFRCAQCobqLZpcvU/RQplZ40XwBKkN3Xg6O1DtrF8H35fF4CAQC8oucEF1UdXU1zz33HB999BEjRoyod0PbokWLInKdqIWywsJCAoEAGRkZdR7PyMggLy8vrHM8+uijVFZW8tvf/rbRYzweDx7P0dvsKyoqWjdgIbqD0J2XWewr2wdA/4T+xFeaSKiu1A9xKvQo0791GA01M2QWI6BSWRbAnpyGCwj4/ZhMJqxWKyZT099qDHEmFIsRzRvAX+rBFVq+jO2ZMrPZHGqQ63a7cTpjr8mtEKLttmzZEurc/+OPP7bbdaJerHHsHQyN7S91rFdeeYX77ruPd955p94GobUtXLiQ+fPnt3mcQnQLoZmyAaGZsqzELAzuAA6fvrTojldx5qegaSqKsSacxVkAHxXFHuJS0wC9ST9AfHx8s5dVFAVjshV/vht/kQu3e1fNa2N7pgz02bLy8nIJZUJ0YcEmsu0taqEsLS0No9FYb1asoKCg3uzZsV577TWmT5/O66+/ztlnn93ksXfddVdo7yqAQ4cOMWzYsNYPXIiurPby5f6NAAxIGED57lL9cRMkGlXivYmguTEb9LsmTXH6t5KK4moSspKhoppAzY0y4Xa6N6XY8Oe7cRXvQVW9GAw27Pa+EXtr7aV2KBNCdB2zZs3ib3/7G/Hx8XVyxLEUReHRRx+NyDWjFsosFgtjxoxhzZo1XHLJJaHH16xZw0UXXdTo61555RWuv/56XnnlFX796183ex2r1YrVag19Xl5e3raBC9FVaVqdULZ/635AX77MqaxpkxFnoJffS5w3EU0twmzQ/20ZTSooEPCpGGzxUFGN36vPrIUzUwZHi/1d5TvAos+SKUrsbzoixf5CdE2bNm0KdWEI7oHZkEj2LIvq8uWsWbOYOnUqY8eOZcKECTz33HPk5OQwY8YMQJ/lOnToEMuXLwf0QDZt2jSeeOIJTjnllNAsm91uJzExMWrvQ4guoaoEPPovLdXODHIrcwF9+bK0Uq+h8MdZ6OUxYNRMBNRKTCZ9pkzRfMQ5LbjLvWgGPVx5a24MCDeUBYv93Z5dYIn9Iv8gCWVCxK7PPvuMhx9+mA0bNpCbm8vKlSu5+OKLw3pt7SXLjlq+jOqvoVdeeSWPP/44CxYsYNSoUXz22We8//779O+vN4vMzc0lJycndPyzzz6L3+/nT3/6E5mZmaGP2267LVpvQYiuI7jnpTOT/VUFaGgkWBJItiYTV6Y3jq2Os5NRre//pprKQjNl+KuJq2kAqyr6XUnVbhfQguXLmpkyt6qPI5a3V6otGDorKyujPBIhxLEqKysZOXIkTz75ZLSHEpaoF/rffPPN3HzzzQ0+t3Tp0jqfr1u3rv0HJER31UA7jKyELBRFwVamN4kttztI9pg4AviUEsyKXtSv+tzEJ6ZSeMCFTzWgqiqeSheapoU/U5asB7xqk75sGsvbK9VWO5SFe6OSEKJjTJ48mcmTJ0d7GGGL/YINIUTHqNUOY3+5HoyyErMAsJTooazI5sTp1UsFfFoJpppCf63aFZop8/pVfD4fqteH6vdhsVjCurwpxYam+PHa9bKEznDnJegzgQaDHkSrq6ujPRwhuryKigrKy8tDH7XbXnV2EsqEELraM2U1PcqyErIAMBQXApBrTcbiSwLAHygPLV+q7griazYV93i8+Hw+FC0APm/YlzfYTPhSj4AhgNHgwGrNbPNb6giKooSWaF0uV5RHI0TXN2zYMBITE0MfDe0A1FlFfflSCBEjavcoO/gOcHSmTCkqRAP2mVMZ7te71iuBqlBLDM1VQnzNTJnHW62HMjWAWt2y4nd/zyMAxBkGdKplQIfDgcvlorKykh49ekR7OEJ0adnZ2fTu3Tv0ee0OC52dhDIhhK4mlGlJ/dn3k/73/gn90VQVtVAPSzm2dKoCGgDGgDc0UxYoLyGuZqbM6/fg9XpBDaBVVbVoCN6kwwDYfP3a+m46lBT7C9FxnE4nCQkJ0R5Gu5BQJoQAvxfK9F5kxXFJVPgqUFDo5+xHoLgY/H5UFPZYM6l0V6NpGmafiknRZ8f8pYWhmTK/6sXr9aIEWj5T5ok7BIDF1buZI2OLw+EAZPlSCNE2EsqEEFB2ANDAZGdfQA8WvRy9sJlsVOXvBqDc5qDEkIArEA9aNQaN0PKlWlxAfKIVDQ0VH16vF5Oq4q9sWUipNu6HAJiLekb07bU3mSkTIja5XC527doV+nzv3r1s3ryZlJQU+vWLvRl5CWVCiKN3XqYMYF/50U7+AP78AgBK45IwaQoezYmm6cuZptDyZTE2q4amBAioAfx+P2ZVxVdRFvYQVNVDdUCfrTMeTutU7SWCocztdqOqKgaD3EMlRCz4/vvvOeuss0KfB7dLuuaaa+q13YoFEsqEEEcbxx7TowzAX5APQIUzGYeqhyRNK0NBCc2U4a+GshIs8eD1ejEYDJgNCpXFRWEPodK9F1Ax+OIwupyo5V6MiZ2jgNdms2EymfD7/bjd7tByphAiuiZNmoSmadEeRtjk1zkhRN07L4OhrObOS1++HsrcCck4ND2UKRwJ9SgD0PzV+IuKsDj1UGY2mbGYTLjLSvF7w2uLUenaAYDV0wcFBd+Rlt0kEG3B2TKpKxNCtJaEMiFEgz3Kjl2+rE5MDc2UmSgKFflrqh9UP4HiYszxml5PZrQQF2cHoKLoSFhDqHTrdR92svTrFnauvSSl2F8I0VYSyoQQoVDmS+zLwQq9rmtAwgAA/DUzZd6kVJw1ocyo1tr3Er/+/4uKMdr0UGZUTKT2SAegvDDMUFa5E4A460B9LPmdK5Q5nU4AysvLozwSIURnJaFMiO5O00KF/odscfg1PzajjYz4DOBoTVkgNS20fGlUXaF6MkUJ6M+XFKNYAnpNmWYktad+B2VF2KFMnylzJJ8AgC+3c93JGAxlFRUVUR6JEKKzklAmRHdXWQi+SkBhPz5AX7o0KPq3B1/N8qWSlo5D1V+iqdWYDTYAVKMeyvxFxQQUDz6fDzQT6b30XmPhzJSpqoeqKv2uz8TMYfp1cys7VYFusJmly+XqVOMWQsQOCWVCdHfBOy8TerOvUm/eGqwnU6uqUGuW48wZGaGaMn/AjyXYDsNUM1NWVITbW4GmaRg1Cz0yg6GsoNkhuN370LQAJpOTuMz+YFTQPAECJZ1no2G73Y7RaERVVelXJoRoFQllQnR3tYr895bpAS1452Wwnkyx24lLScShKWiaRnUALDUzZT5TTU1ZcTEut96XzKLE40zT94AMZ/nSVanfeRkfPxiDyYg5Xd/g25fbeYrmFUWRJUwhRJtIKBOiuwuGspT6PcqCS5fm9HQS7OaamTIvAU0JLV9WWfQlz0BREeWuUv14LZ6EmkL/cO6+DBb5x8cN0l+fqbeXkLoyIUR3IqFMiO6u+Gjj2P013fyPbRxrysjAgQETCpqqByWLSS/0r7DqS4z+4mLKKkr154jHnpAK6DVlmqo2OYRgkX+843ig84cyuQNTCNEaEsqE6O5qZspcCb0orCoEajWOzcsDwNQzA6tfL173ogclq0mvKSu26J8HiotDYSTOmIjB6ABFIeDz4S5verul+jNles8v7+HOs3wJR4v9ZaZMCNEaEsqE6O5qCv33WfSZr1RbKk6LPuMTbBxrzuiJ2auHsmr0oGQz6scXmPUAUu31hgrc7aYEPG4NR3IK0HRdWe07L+MdgwGw9HGAAoESD4GK8HYEiAXBmbLKykoCgUCURyOE6GwklAnRnXndUJELwF5FL9g/Lum40NP+/KMzZYYqPWR4VD2UWWtqynKNxShxcRSZ9K10LRYLJsWGu8wTKvZv6g7M4J2XRqMDq0XvjWawmTDVFPt7czrPrJPNZsNisaBpmixhCiFaTEKZEN1ZkV7LhT2FPVV6cAp28gfw5ek1ZeaePVEra+6yrAlltpqWGLn+PEwpKZQYjfqpbA4UFNzlXhJSg6Gs8Zmy4J2XjvjBKIoSetzaT18K9B7oXOEmKSkJgLKyppdshRDiWBLKhOjOCvVARNrgUDuMAYlHQ5k/WFOWnkF1hX6XJcFCf0Xf27JAK8KQkkxRTShLiNfDVGWZ5+gdmE2EslCRf/zgOo9b+ulLgZ1ppgwgMTERkFAmhGg5CWVCdGfBmbIGQpnm8+Ev1Av/zT0zcJfqd1kaakKZUdGXF10GN/7EeEoM+ixXUkIyAO5yb1jLl6Ei/8ZC2YEKtEDn6ZAfDGWlpaXRHYgQotORUCZEd1YzU+ZLGUhORQ5wNJT5Cwv1fTFNJoypqbhqQplZrURBQamZKXMZ3bjiTFQY9G8nqalpALjLvCSkhbF86doGgMMxpM7jph5xKDYjmk/tVE1kg8uXFRUVqM20AhFCiNoklAnRnRXqs1SHHKn4VT92k52e8fpG4sF2GOb0dBSDgcqaUGYNVIYax4IeyvJNKj6zBbPZHAol7nIvCWlNL1/6/ZVUVelh0OE4oc5zikHBOkCfdfLs7jxLgXa7XYr9hRCtIqFMiO5KVUPLl3st+p2TWQlZoY3Ig1ssmTIy8HsDeNx+NM2HWfOFtljyKtUEFJVcxYPPYsZms+FM1HuMVda6+7KqohxfdXW9IVTWFPlbLD2wWFLrPW8dmARA9a6SSL3rDiFLmEKI1pBQJkR3VXEYfG4wmNgbcANHm8ZCrVDWMyO0dKkGe5TV9DHzGvQi/MNaFV6zBbvdTlKyHkiqKnyYbXFY7PoyZ0N1ZRWhpcsT6j0HYBucpF9nXzmav/MsBQZnC0tKOleYFEJEl4QyIbqrmqVLkgewp2Z7pQbbYWT0DC1dqiZ9tstq15vCBhR9WfGQ34XHasFms5HaIwkU0FQNT6WfpIxeAJTk5dYbgsu1HWg8lJnS4zA4zWg+Fc/+zrMUmJKif30klAkhWkJCmRDd1ZGf9T97DGFHib6MODj56B2QocaxGRm4SvRQhlmfUTNZ9NChKGWoXpV8tRKf1YrVaiUhMQG7wwyAu9xDUqYeykpzD9Ubgsulj6GxUKYoCraaJUzPzs4TcJKT9TtQKysrqW5g2VYIIRoioUyI7ir/RwD86UPZXbobgCHJR++A9B7SQ5S5dy8qivVgYTRV1fypL1+alDICrgCFxiqMCQkYDAacDgdxiXpj2coyL8k9gzNlh+tcXtO0UChzOoY2OkzbED0AVmUXteHNdiyz2RzaB1Nmy4QQ4ZJQJkR3lf8TAPsTM/CqXuJMcfR29g497TsYDGW9qSjSQ5nJWBPKDPEA2LRS/C4/JYFq7HY7hkAAq89HfIK+L6a7zEtycKbsmFBWXX2IQMCFopiJixtAY2wnpIBBwV9Qhe+IOxLvvEMElzCLi4ujPBIhRGchoUyI7kgNQIE+S7XdrN95OTh5cOjOS9XtJlCkz0xZ+vShoqgmjNVssWRCL96P10oJuAJoqoLdbsdWXY1aXExcUs1MWamHpOBMWW7dmrLgLFl8/EAMBkujQzXYTVgH1tw88FPnmS2TUCaEaCkJZUJ0R8V7wV8F5jh2eEuBukuXvpqlS4PTiTExkfKamTJ8emG/CT10KYEirD4DFlW/89LursJfVIwzWX++oqQ6NFNWUXQEn+dofVVzd17WZh+uN6TtjKGsrKwMn88X5dEIIToDCWVCdEc19WSkD2V7qV7kf3zy8aGnvQcPAmDu0wdN1UI1ZX53KQA29EL+kqoiEhUz8eZ4jEYj9qoqAsVFOFL0PmauomrszgSs8fpyZ2nNzQMAFRVbAXA6T2x2uPbhqaCA70BFp1nCtNvtxMfHo2kaRUWdJ0wKIaJHQpkQ3VEolA1jR7Eeyoak1J4p0+u/LH164y73ovo1UFSqy0sBcCp6KCt2F5OoGXDa9cJ/e1XNTFmqHsoqiqtRFCVU7F+ae7SurLxcD2UJzpOaHa7RacF2fM2emhsb30cz1vTooTfPPXKk8W2mhBAiSEKZEN3RoY0AHOkxmIKqAhSUOu0wfMGZst59QkuXcQ4vmqYSUIwk1XzrKHKXkaiqJNmS9GOq3PgLCnCmHA1lmqaF6sqKa9piVHvy8HoLAANO5/Cwhhw3JgMA98Z8NLVzbFCenq5vMyWhTAgRDgllQnQ3mgaHNgCwJU5fVhyUPIh4c3zoEN+hYCjrHSryt8bpf6rWVMyKQoWnEq/fRYrqI8mSBIDdXYUvLxdHTU2Z36tSXekjtXdfAIoO6vtcVtTMkjniB2M02sMatn1oKordRKDM22l6lqWmpqIoCpWVlVRWVkZ7OEKIGCehTIjupmg3VJeCycYWv164PyJtRJ1DvPv2AWDp34+yI3oYM1v1Wi5rzYblBd4SFMXPcaZqHKq+36XD5cKfl4/JbCSupi1GRVE1af31lheFOfp5y8u3AOBMqHvdpihmA/Gj9Zkn11eHmzk6NphMplDBv8yWCSGaI6FMiO7m0Pf6n5kj2Vqk9yob0eNoONL8frz79G2XLMcNpCRPD2MGgz7TExevB6NDNXdiDjJ5sQdqWmRUVuKrKeZ31FrC7NEvC4CigwcI+P2UV4RfT1abY2IvUKB6ewm+/M4x8xRcwszLy2vmSCFEdyehTIju5qAeyvy9TubHQr3gv/ZMme/QITSfD8Vqxdwrk9J8PZRpqr73ZLI9hWq/h1xvBZjjGWCLR0EhgB+rx4M/Lx9N047WlRVVk9AjHYvdjhrwU3z4IBUV+nUTEloWykxpdmzDUvXzfl5/26ZYlJmZCUBhYSFerzfKoxFCxDIJZUJ0Nwe+AWBnah+q/FXEm+MZkHi0o75nzx4ALFlZoCihUOZxFQKQHJ9KgauIcgVS0tKJN+jbCVWaK1EAzeMhUFpKQpoeysqPVKEoCml9swDI3f81Pl8JBoMVh+PoHZ/hcp6u7zrg3lSAvzj295WMj48nISEBTdNktkwI0SQJZUJ0J+5iyNOXDr81KQCcnH4yRoMxdIh3z14ALMcNwF3uxecJoChQUaQHimRrAvmuIopVlcw+fSlH77bvMleiJesBzZ+XR1JGHAAlNaGuR/8sAIoLvwYgIWEUBoO1xW/BmpWIdVASBDTKP9rf4tdHQ69e+t2nucfsaiCEELVJKBOiO9n3BaBBjxP4plivJxufOb7OIZ49+ubk1gEDKK2pJ4tPMVBZqm8XZAkolFaXc8QfILPvAEpqQlmlqRJvqt6vzJebR3JNKAvOtKX102fjqnzZACQl/aLVbyPxvCxAny3z5cV+bVlwCfPIkSN4PJ4oj0YIEasklAnRnez9DABf1mlsyNfbYtQLZdv0PSmtQ06g8JC+16XDqQcre0IixTWF/CXGOGxp/ShEb+paaa6kPEUv+PcdPEBSTz2UuUo8+DwBeh43CADFpteCtSWUWfo6sZ+YChqUvL0r5vuWORwOkpKS0DSNgzU94IQQ4lgSyoToLjQNdn8MwNa0/rj9bhKtiXW2V9K8Xqp37gTANnwYhQf1UGaN02ejkjJ6kVvTw6zImUKlvSdH0Avvy83l5KXqS6LeffuxOyxY4/XNzksL3PTIGoA9CcwOL2AkMWF0m95O4q+PQzEb8O4rx70hv03n6gj9+/cHICcnJ8ojEULEKgllQnQXR7ZD8R4wWvhE08PWab1Pw6Ac/Tbg2b0bfD4MCQmYe/em8EAFAIpSCoDZlkK1pxqDwUixI5lcJSO0fFluKWe3Q59R8+7Xa71qL2EaTWYyR+oNak30xWQ62qy2NUzJNhLO0YNO6ft78ZfG9rJgr169MBqNuFwu2QtTCNEgCWVCdBc//wcALesMPjn8BQC/6verOodUZ+v1XrahQ1EDGsWH9Rmy6gq9QL3Kq898xSVloCoGdrjiAQUrVXgMHn6063VnweazKb30prLBGbfE/npo8xVnROQtOU7tjbmPA63KT/GrP6MFYncZ02Qy0adPHwB2794d5dEIIWKRhDIhuott7wGwc8B4DlQcwGKwcGqvU+scUvWD3mnfNmwYxYcrUQMa1jgTJbn7cXt9+P1GDIqBxEw9XOwq0vtuZVCExWhmf5L+uS83F9XrpUc/vfD/SE4FquoH+z4A8rLViLwlxaiQetUJKFYj3n3llK3eF5HztpeBAwcCkJ+fT3l5eZRHI4SINRLKhOgOCrZB7mYwmHjPoC/zndr7VOLMcXUOc3/3HQBxY8dweGcpAGl9zZQV5FNYUYlVsdPTkYa5h9764khBAQDpFJFlS6MsDtQ4G6gqvpwc0vvXhLL9FZSVbULDjb/aSH52BeVHCiLy1kypdpIv0zdTd316ENc3sdt2Ij4+PtQeY2dN7Z4QQgRJKBOiO9j0IgC+wefx7oGPALhk0CV1DvEfOYJ3715QFOLGHA1ljqQKKj1e/FYbVATISu6NuebOSm+53lC2N3kMNMSBolDZKwkAz86dpPZyYDAqVFf6OHxwtX6dskzQFPb9sDFiby9uRA+cv+oHQOk7u3Bvid19JgcP1gPk4cOHKSnpHBurCyE6hoQyIbo6rxt+eAWAtf1GUFxdTKotldP6nFbnsOAsmXXIEAwJiaFQpvoOkVdWQWpmbzLVJOItcSQNSEJBw+rV97/sQx4nBfTz7O+pN6Kt/uknjGYDqb0dgMqRwlUAJMafDsDezd9H9G0mnN2PuDEZoELxKz9T+X1s3pGZkJBA3759Afjpp5/QtNitgxNCdCwJZUJ0dZtWgLsILakfi4v12anLjr8Ms8Fc57CKT9YCEH/KKRw5UEF1pQ+TxcDP2d9S7feT2XsgWfZeYIB+x6eSrLgxomI1GUilmJGl+uzUxmQ9qFX9qDen7X18Eva03QS0fIxGBwOH/w6AfZs34nFHrvGroigkXzaY+F/01PuXvbGDsg/2xmQPsxNOOAGj0UhJSQn7am6KEEIICWVCdGVeN3z5TwC+GHER2cXbsBltXD306jqHqV4vrnXrAHCeey57NusBKyXLwE9b9OL/8cNPw2I0Y06PI8FpZZijCoDE1B4YgBPytmM2mNnaQ3+8+qef0AIB+g5LITFL31qpR49z6XncUFJ698Xv87Ljmy8j+nYVg0LSpYNwnKnfiFCx7iCF/96Kv6gqotdpK5vNxrBhwwDIzs6moqIiyiMSQsQCCWVCdGVfPgHlB/El9uUfZZsBuHLIlaTYUuocVvnFF6guF6b0dGwjR7B74xF8fi8HCr/E5/OS3qMHAyz6kptlgN6XrL9Jv3uw0tEX7ClYVB/Dnf05mAaBOCtqRQXV2dn06K+S0H89AA7zb1AUhWGnnwXAT+s+jvhbVhSFpMkDSLlqCIrZgGdPGfmPb6R87QFUbyDi12utrKwsevTogaqqfPfdd3i93mgPSQgRZRLKhOiqcrfAF48B8Ozws9hXvp8UWwp/HPnHeoeWvPQyAAm//jWHd5VTeLiMgyXbKS/ejdVkZNK5F+DZpS9L2gYmUVlZibla70m2y+OA3icDcKopmYBRYd/x+t2Zrs8/J6/gVQxGH1XF/Tm4Vd8DctgZv8RgNHLo55/I3bm9Xd5+3Mh00m87GetxiWg+lfIP95H3j++o+OwggUpfu1yzpU4++WTi4uKorKzk22+/xeeLjXEJIaJDQpkQXVFlIbx+DQQ8fDboNJ7L/RSAO8fdidPirHNo1Y8/Ufnll6AoJP/uKr5+72f25P9EUl8o2L2dvilJDBt+Kv7CKjApWAclsWXLFtA0jqjxfHXAgz9rEgBnFentKNb10gNc6ef/ZX/OvwEo3n4OO9bn4/cFcKamMfQ0fbbsq9dfardid3OanbQbTyL5t8djTLGhunyUvb+X3IXfUPTKz7i3HkGt9rfLtcNhsVgYN24cZrOZkpISvvrqK6qqYmupVQjRcaIeyp566ikGDBiAzWZjzJgxfP75500e/+mnnzJmzBhsNhvHHXcczzzzTAeNVIhOoiIPVlwMxXv4Mq0vs7Q8NDQuG3wZkwdMrnOo5veT/+CDACRceCEbd5Sy7pNP8fqrsNvz6eOMJ7P/AJwlepCzDUlBNcG3334LQL6pJxUePxus+ubix+//juMSsvhycADVbOTIyJ8JBFw4HMPRXKfhLveS/cVhAMZfcgVGk4l9P2xk+1eftduXQ1EU4k/OoOcdY0i+bDDmzHjwa1T9cITil37m8N/WU/DMD5S+v5eqHwvxFVZ16M4ATqeTiRMnYrVaKS8v59NPP+Xw4cMddn0hurqW5oxoimooe+2115g5cyZz585l06ZNnH766UyePLnRDXv37t3LBRdcwOmnn86mTZu4++67ufXWW3nzzTc7eORCxCBNg5/fh2fPpDr/R/43PZObnUY8qpdJfScx95S5xxyuUbDoMdwbNlBsNrNlxDhe/t9VBFQ/J4xORzn4M1azifFnX4b7O729hOOUTL799ltKSkqIj4/nhOEnAvDiTgukDUFRfVxtz8IVp7DtCqgap4IKQ46/lzHnDwDgm3f2UF5YRXJmb8ZdfAUAq597kvw9u9r1y6MYDcT/oifpt44m/U+jcJzeG1OaHQIa3n3luD47SNGL28h/5HsOzfuSvEXfU7j0J0re2knZ6n241h+m6sdCqneW4D1Qga/ATaDcg+rxt/kOz4SEBE477TQSExPx+Xxs2LCBr776iiNHjkjLDCHaoKU5I9oULYr/4sePH8/JJ5/M008/HXps6NChXHzxxSxcuLDe8X/9619599132bZtW+ixGTNm8MMPP/D111+Hdc2DBw/St29fDhw4ENqHTohOzXUEdq1B+/4F9uVvYnV8HK8nJpNv0P9pXzzoYu455R4sRkvoJVX5+ey6/372/mcVhX4/1Rdfx47cOLxuleMGpJNs+Y7SwwfpM3A4k3r/D74DLixZCRSdaeH1119HVVV+85vfYMscxK//+QUGBT4/axe9v5qH25nOoyckMsFWgqKA810j/XvdSNrM23n78c3k7SknIc3G5BkjSO5p462F88j5cQtmq42zrvsDw8/8FQaDscO+fL7CKrz7yvHmlOth60gV+FuxDZQCGA0oRgXFqBz9u8mgP6cAiqIfqtS8oPbjCqiqyp7CA+w7cgAVFVCIs9pJT0glzZlCgt2B2WhuZAAxSIn2AES0GOLMpE0bFtFztubnd0tzRrRFLZR5vV7i4uJ4/fXXueSSo53Fb7vtNjZv3synn35a7zVnnHEGo0eP5oknngg9tnLlSn7729/idrsxm5v/ZtVeoayipIjv31tR73FVA02r+w1e/5JroNV5sO4Ltdp/0Wr+D/TvctqxB6EGX9/gabR6jx99XgNNq3W5Y89T64Va8Dz1T6bVeUypOWfD163zuFL31Zp29Dmt5lqhd1zvXFpjb6v+1xMANfSWGv+P/ujXtrEvQd0XN/a10Dj2J5JW9//VfX1DA1JqDVbTQFNBU1FVP4GAD3/Ag6r68KoBvApUKQb8NT/0VU3FqJlIs6fhMDnxe31Uu9143FVUlpbiLnOhYkQ1WfHZEwADqU47GWkWTFo5Rsw47Sn0dGTh8/ipNHso7unlUJ5eM3bccQM466yzUBSFR1dv55s9+QxILOFPSa9RlliNx6qHqt27jZz2qAEFBU+vVFzjTuHnql/i8dgAjaQML/EpleRvX42rSF+ys9jjSerTF0dqGpa4OMw2GwajCcVgQFEMNX82+D9eM5p4Ue2nNDB5jZg9ZsweI0afEZPPiNGv/90QUDCoBv3PgP7eIq3K52F/6SFyywvwa3XvGI0z27CbbNjMVmwmKyaDEZPBiNFgxGQwYVD0pGdQFBRFQQn+PcxxKq364kr6EnWpFpUhC8+N6Dlb+vO7NTkj2kzRunBhYSGBQICMjIw6j2dkZJCXl9fga/Ly8ho83u/3U1hYSGZmZr3XeDwePB5P6PP26geUve49vsx5qF3OLUQdhpqPJn4HMaBnuSPAERX9X3pCzUdPsAJWm4GEBAMpKUYSE40YjfV/sAb/tRiBHkCPWv/ENmzUf/M8M1X/ACioOdrqCeA/6ONFYyI/XKhw3RqV+MNFWN9exS/Mn7H9+Cs50mM0pflWSvOtaNoVmOyb8Fd/i7eqkoKdP1PQSbaGNComjIoJg2LEgEH/M/iBvoG7ouiVIkrN1JgS/ExRaj1e/zFDnEZVtZuyKjcuTxXeQPRuSuhOJGK2nclo4UEiG8qCKioqKC8vD31utVqxWq31jmtNzoi2qIWyoGN/K9M0rcnf1Bo6vqHHgxYuXMj8+fPbOMowKZFZcmnJb96Nzfi0x2/vzWvhNYOHNzlX2/Q5m/6lvvEnw50frn/+ln9dFQXQWve/R2hCscGXN35OoxEMBgWDQcFoBJNRwWIxYLEoWK0G4uJMmExHX+/3Qt0f90eX1hRFwWgwYDSZsFisGI31S1E1jBRXJ7GvJIWd+b0YVbCHMepP/NtQwleZKq9eYyFpr4HeuQppJRWk5f6bHnkZlCWOwB3fH58lCZ85i4BtAH4tl4B6hIBaAqobTasGAmjoM4UQ6V5jbVssCGgaAc0HtF87C4sdUuzgD6h4/D58ARVfIIAvEEBVNQKahqqpqGrNPG1w9lvTjk60tvF9NkZK3kRDzCZ7u5072Hg56N577+W+++5r9PiW5oxoilooS0tLw2g01kurBQUF9VJtUM+ePRs83mQykZqa2uBr7rrrLmbNmhX6/NChQ/X+B42E8Zdcy/hLro34eYXo7I6P9gCEEF1KdnY2vXv3Dn3e0CwZtC5nRFvU7r60WCyMGTOGNWvW1Hl8zZo1TJw4scHXTJgwod7xq1evZuzYsY3Wk1mtVhISEkIfTqezweOEEEIIEfucTmedn+uNhbLW5Ixoi2pLjFmzZvHvf/+bF154gW3btnH77beTk5PDjBkzAH2Wa9q0aaHjZ8yYwf79+5k1axbbtm3jhRdeYPHixcyePTtab0EIIYQQMaq5nBFrolpTduWVV1JUVMSCBQvIzc3lxBNP5P3336d///4A5Obm1uklMmDAAN5//31uv/12/vd//5devXrxz3/+k8suuyxab0EIIYQQMaq5nBFrotqnLBqkT5kQQgjR+XSHn99R32ZJCCGEEEJIKBNCCCGEiAkSyoQQQgghYoCEMiGEEEKIGCChTAghhBAiBkgoE0IIIYSIARLKhBBCCCFigIQyIYQQQogYIKFMCCGEECIGSCgTQgghhIgBUd37MhpUVQX0fTWFEEII0TkEf24Hf453Rd0ulOXn5wMwbty4KI9ECCGEEC2Vn59Pv379oj2MdtHtNiT3+/1s2rSJjIwMDIbIrt5WVFQwbNgwsrOzcTqdET13VyNfq/DJ1yp88rVqGfl6hU++VuFrr6+Vqqrk5+czevRoTKauOafU7UJZeyovLycxMZGysjISEhKiPZyYJl+r8MnXKnzytWoZ+XqFT75W4ZOvVetJob8QQgghRAyQUCaEEEIIEQMklEWQ1Wrl3nvvxWq1RnsoMU++VuGTr1X45GvVMvL1Cp98rcInX6vWk5oyIYQQQogYIDNlQgghhBAxQEKZEEIIIUQMkFAmhBBCCBEDJJS1wgMPPMDEiROJi4v7/+3da1BUBRsH8P/Kyy43XZFFlkUCbFrJIG4V4kR4RRkulUpIGtiMzdgMKoo5NtqAdNNySNOULIfMD6wf0qZRdNTiouEFYasVmkRhIY1LIAiBsgjP+8E808ICGwK7eZ7fzH7gnOccnvOfw/JwdveAiRMnmqyRSCT9HtnZ2UY1Op0OERERsLe3h4eHBzIzM/GovcXPnKxqa2sRGxsLR0dHKBQKrFmzBgaDwahGDFmZ4u3t3e882rRpk1GNOfmJxd69e+Hj4wM7OzuEhITg7Nmzlm7J4jIyMvqdQ0qlUlhPRMjIyIBKpYK9vT1mzZqF8vJyC3Y8doqKihAbGwuVSgWJRIJvv/3WaL052XR1dWH16tVQKBRwdHREXFwcbty4MYZHMTaGymrFihX9zrMZM2YY1Yglq4fBQ9kwGAwGxMfH48033xy0LicnB3V1dcIjOTlZWNfW1ob58+dDpVKhpKQEu3fvxo4dO5CVlTXa7Y+pobLq6elBdHQ0Ojo6cO7cOWg0GnzzzTdIS0sTasSS1UAyMzONzqMtW7YI68zJTywOHz6M1NRUbN68GVqtFuHh4YiKikJtba2lW7O4p556yugc0ul0wrqPPvoIWVlZ2LNnD0pKSqBUKjF//ny0t7dbsOOx0dHRgYCAAOzZs8fkenOySU1NxdGjR6HRaHDu3Dn89ddfiImJQU9Pz1gdxpgYKisAWLhwodF5lpeXZ7ReLFk9FGLDlpOTQ3K53OQ6AHT06NEBt927dy/J5XK6e/eusOzDDz8klUpFvb29I9yp5Q2UVV5eHo0bN45u3rwpLMvNzSWZTEa3b98mIvFl9U9eXl70ySefDLjenPzE4rnnnqNVq1YZLfP19aVNmzZZqCPrkJ6eTgEBASbX9fb2klKppG3btgnL7t69S3K5nLKzs8eoQ+vQ9znbnGxaW1vJ1taWNBqNUHPz5k0aN24cnTx5csx6H2umfr8lJyfTiy++OOA2Ys3q3+IrZaMoJSUFCoUCzz77LLKzs43+s/358+cRERFhdB+XBQsW4I8//oBer7dAt5Zx/vx5+Pn5QaVSCcsWLFiArq4ulJaWCjVizmr79u1wcXFBYGAg3n//faOXJs3JTwwMBgNKS0sRGRlptDwyMhLFxcUW6sp6VFZWQqVSwcfHB0uXLkVVVRUAoLq6GvX19Ua5yWQyREREiD43c7IpLS1Fd3e3UY1KpYKfn58o8ysoKMDkyZOhVqvxxhtvoLGxUVjHWZnn0fyPnlbg3Xffxdy5c2Fvb4/vv/8eaWlpaGpqEl56qq+vh7e3t9E2bm5uwjofH5+xbtki6uvrheN+wNnZGVKpFPX19UKNWLNau3YtgoOD4ezsjEuXLuHtt99GdXU1vvzySwDm5ScGTU1N6Onp6ZeFm5ubqHIwJTQ0FF9//TXUajUaGhrw3nvvYebMmSgvLxeyMZVbTU2NJdq1GuZkU19fD6lUCmdn5341YjvvoqKiEB8fDy8vL1RXV+Odd97BnDlzUFpaCplMxlmZia+U/c3Um2H7Pi5fvmz2/rZs2YKwsDAEBgYiLS0NmZmZ+Pjjj41qJBKJ0df09xvX+y63NiOdlanjJSKj5f/VrEz5N/mtW7cOERERePrpp7Fy5UpkZ2fjwIEDaG5uFvZnTn5iYeo8EWMO/xQVFYXFixfD398f8+bNw/HjxwEABw8eFGo4t4ENJxsx5peQkIDo6Gj4+fkhNjYWJ06cwNWrV4XzbSBizGowfKXsbykpKVi6dOmgNX2v1vwbM2bMQFtbGxoaGuDm5galUtnvr4MHl3r7/mVmbUYyK6VSiYsXLxota2lpQXd3t5DDfzkrUx4mvwefZrp27RpcXFzMyk8MFAoFbGxsTJ4nYsrBHI6OjvD390dlZSVeeuklAPev+Li7uws1nBuET6gOlo1SqYTBYEBLS4vRFaDGxkbMnDlzbBu2Mu7u7vDy8kJlZSUAzspcfKXsbwqFAr6+voM+7Ozshr1/rVYLOzs74bYQYWFhKCoqMnp/0KlTp6BSqR5q+BsLI5lVWFgYrly5grq6OmHZqVOnIJPJEBISItT8V7My5WHy02q1ACD8kjAnPzGQSqUICQnB6dOnjZafPn2an/D76Orqwq+//gp3d3f4+PhAqVQa5WYwGFBYWCj63MzJJiQkBLa2tkY1dXV1uHLliujza25uxu+//y48V3FWZrLYRwz+w2pqakir1dLWrVvJycmJtFotabVaam9vJyKi7777jvbv3086nY6uXbtGX3zxBU2YMIHWrFkj7KO1tZXc3NwoMTGRdDodHTlyhCZMmEA7duyw1GGNiqGyunfvHvn5+dHcuXOprKyMzpw5Q1OmTKGUlBRhH2LJqq/i4mLKysoirVZLVVVVdPjwYVKpVBQXFyfUmJOfWGg0GrK1taUDBw5QRUUFpaamkqOjI+n1eku3ZlFpaWlUUFBAVVVVdOHCBYqJiaHx48cLuWzbto3kcjkdOXKEdDodJSYmkru7O7W1tVm489HX3t4uPCcBEH7eampqiMi8bFatWkVTpkyhM2fOUFlZGc2ZM4cCAgLo3r17ljqsUTFYVu3t7ZSWlkbFxcVUXV1N+fn5FBYWRh4eHqLM6mHwUDYMycnJBKDfIz8/n4iITpw4QYGBgeTk5EQODg7k5+dHO3fupO7ubqP9/PLLLxQeHk4ymYyUSiVlZGQ8crd4GCorovuDW3R0NNnb29OkSZMoJSXF6PYXROLIqq/S0lIKDQ0luVxOdnZ2NG3aNEpPT6eOjg6jOnPyE4vPPvuMvLy8SCqVUnBwMBUWFlq6JYtLSEggd3d3srW1JZVKRYsWLaLy8nJhfW9vL6Wnp5NSqSSZTEYvvPAC6XQ6C3Y8dvLz800+PyUnJxORedncuXOHUlJSaNKkSWRvb08xMTFUW1trgaMZXYNl1dnZSZGRkeTq6kq2trb02GOPUXJycr8cxJLVw5AQieC26IwxxhhjVo7fU8YYY4wxZgV4KGOMMcYYswI8lDHGGGOMWQEeyhhjjDHGrAAPZYwxxhhjVoCHMsYYY4wxK8BDGWOMMcaYFeChjDHGGGPMCvBQxhhjjDFmBXgoY4yNuYKCAkgkErS2tlrk+//www/w9fVFb2+vsGz//v3w9PTEuHHjsHPnzkG3P3bsGIKCgoy2Z4yxh8X/ZokxNqpmzZqFwMBAo0HHYDDg1q1bcHNzg0QiGfOennnmGaxduxavvfYaAKCtrQ0KhQJZWVlYvHgx5HI5HBwcBt1HcHAw1q9fj+XLl49Fy4wxEeArZYyxMSeVSqFUKi0ykBUXF6OyshLx8fHCstraWnR3dyM6Ohru7u5DDmQA8Prrr2P37t2j2SpjTGR4KGOMjZoVK1agsLAQu3btgkQigUQigV6v7/fy5VdffYWJEyfi2LFjmDZtGhwcHLBkyRJ0dHTg4MGD8Pb2hrOzM1avXo2enh5h/waDARs3boSHhwccHR0RGhqKgoKCQXvSaDSIjIyEnZ2d8L39/f0BAFOnThV6/PnnnzF79myMHz8eEyZMQEhICC5fvizsJy4uDpcuXUJVVdXIhsYYE63/WboBxtija9euXbh69Sr8/PyQmZkJAHB1dYVer+9X29nZiU8//RQajQbt7e1YtGgRFi1ahIkTJyIvLw9VVVVYvHgxnn/+eSQkJAC4f7VKr9dDo9FApVLh6NGjWLhwIXQ6HZ544gmTPRUVFSExMVH4OiEhAZ6enpg3bx4uXboET09PuLq6IiYmBkFBQdi3bx9sbGzw008/wdbWVtjOy8sLkydPxtmzZzF16tQRTI0xJlY8lDHGRo1cLodUKoWDgwOUSuWgtd3d3di3bx8ef/xxAMCSJUtw6NAhNDQ0wMnJCdOnT8fs2bORn5+PhIQEXL9+Hbm5ubhx4wZUKhUAYMOGDTh58iRycnLwwQcfmPw+er1eqAcAe3t7uLi4ALg/MD7os7a2Fm+99RZ8fX0BwOSQ5+HhYXLAZIyx4eChjDFmFRwcHISBDADc3Nzg7e0NJycno2WNjY0AgLKyMhAR1Gq10X66urqEIcuUO3fuCC9dDmb9+vVYuXIlDh06hHnz5iE+Pt6oP+D+QNfZ2WnW8THG2FB4KGOMWYV/vjQIABKJxOSyB7eh6O3thY2NDUpLS2FjY2NU989Bri+FQoGWlpYh+8nIyMCrr76K48eP48SJE0hPT4dGo8HLL78s1Ny6dQuurq5D7osxxszBQxljbFRJpVKjN+ePlKCgIPT09KCxsRHh4eH/aruKigqzatVqNdRqNdatW4fExETk5OQIQ9ndu3dx/fp1BAUFDat/xhjriz99yRgbVd7e3rh48SL0ej2amppG7IararUay5YtQ1JSEo4cOYLq6mqUlJRg+/btyMvLG3C7BQsW4Ny5c4Pu+86dO0hJSUFBQQFqamrw448/oqSkBE8++aRQc+HCBchkMoSFhY3I8TDGGA9ljLFRtWHDBtjY2GD69OlwdXVFbW3tiO07JycHSUlJSEtLw7Rp0xAXF4eLFy/C09NzwG2WL1+OiooK/PbbbwPW2NjYoLm5GUlJSVCr1XjllVcQFRWFrVu3CjW5ublYtmyZWfc0Y4wxc/Ad/RljorNx40bcvn0bn3/++bC2//PPP+Hr64vLly/Dx8dnhLtjjIkVXyljjInO5s2b4eXlNez3ulVXV2Pv3r08kDHGRhRfKWOMMcYYswJ8pYwxxhhjzArwUMYYY4wxZgV4KGOMMcYYswI8lDHGGGOMWQEeyhhjjDHGrAAPZYwxxhhjVoCHMsYYY4wxK8BDGWOMMcaYFeChjDHGGGPMCvwfqNkokh0WLeEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot charge population\n",
    "fig, ax = plt.subplots()\n",
    "ax2 = ax.twinx()\n",
    "for q in range(chargePops.shape[0]):\n",
    "    ax.plot(times / FS2AU, chargePops[q, :], label=f'q={q}')\n",
    "\n",
    "ax2.plot(times / FS2AU,\n",
    "         fluence * pulseEnvelope(times) * 1 / AU2FS * 1 / ((1e6 * AU2M)**2),\n",
    "         label='pulse',\n",
    "         c='black',\n",
    "         alpha=0.3)\n",
    "ax.legend()\n",
    "ax.set_xlabel('time (fs)')\n",
    "ax.set_ylabel('population')\n",
    "ax2.set_ylabel(r'intensity (n. photons / $\\mathrm{\\mu m}^2$ fs)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:22.521270Z",
     "start_time": "2022-05-31T09:13:22.504524Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Singly core ionized configuration has electronic occupation numbers [1 2 4]\n",
      "The probability for PP is 0.40088334372339834.\n",
      "The probability for PA is 0.4820026678699649.\n",
      "Ratio PP/PA is 0.8317035785195053\n"
     ]
    }
   ],
   "source": [
    "# ratio of PP vs PA processes:\n",
    "# i.e. either A or P starting from configuration 1,2,6\n",
    "# determine the singly core ionized configuration\n",
    "cfgP = [\n",
    "    p['fconfig'] for p in processData\n",
    "    if np.all(p['iconfig'] == np.array(configurations[startCfgIdx]))\n",
    "    and p['type'] == 'P'\n",
    "][0]\n",
    "print(\n",
    "    f\"Singly core ionized configuration has electronic occupation numbers {cfgP}\"\n",
    ")\n",
    "# index of singly core ionized configuration\n",
    "idxP = np.argmax(np.all(configurations == cfgP, axis=1))\n",
    "# find the total photoionization c.s. and the Auger-Meitner decay rate for configuration idxP\n",
    "sigmaPP = np.sum([\n",
    "    p['pcs'] for p in processData\n",
    "    if np.all(p['iconfig'] == configurations[idxP]) and p['type'] == 'P'\n",
    "])\n",
    "GammaPA = np.sum([\n",
    "    p['rate'] for p in processData\n",
    "    if np.all(p['iconfig'] == configurations[idxP]) and p['type'] == 'A'\n",
    "])\n",
    "# probability for PP is\n",
    "# int dt sigma * F(t) p_P(t)\n",
    "# probability for PA is\n",
    "# int dt Gamma * p_P(t)\n",
    "dt = times[1] - times[0]\n",
    "probPP = np.sum(sigmaPP * fluence * pulseEnvelope(times) * p[idxP, :]) * dt\n",
    "probPA = np.sum(GammaPA * p[idxP, :]) * dt\n",
    "\n",
    "print(\n",
    "    f\"The probability for PP is {probPP}.\\nThe probability for PA is {probPA}.\\nRatio PP/PA is {probPP/probPA}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[{'type': 'P',\n",
       "  'i': '45',\n",
       "  'f': '44',\n",
       "  'pcs': 0.008532284158464506,\n",
       "  'energy': '213',\n",
       "  'iconfig': array([2, 2, 4]),\n",
       "  'fconfig': array([1, 2, 4])},\n",
       " {'type': 'V',\n",
       "  'i': '45',\n",
       "  'f': '42',\n",
       "  'pcs': 0.0003873479622629958,\n",
       "  'energy': '721',\n",
       "  'iconfig': array([2, 2, 4]),\n",
       "  'fconfig': array([2, 1, 4])},\n",
       " {'type': 'V',\n",
       "  'i': '45',\n",
       "  'f': '36',\n",
       "  'pcs': 0.0001015465060570277,\n",
       "  'energy': '736',\n",
       "  'iconfig': array([2, 2, 4]),\n",
       "  'fconfig': array([2, 2, 3])}]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[p for p in processData if np.all(p['iconfig'] == np.array(configurations[startCfgIdx]))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[{'type': 'A',\n",
       "  'i': '44',\n",
       "  'f': '39',\n",
       "  'rate': 0.0009139013893710296,\n",
       "  'energy': '513',\n",
       "  'iconfig': array([1, 2, 4]),\n",
       "  'fconfig': array([2, 0, 4])},\n",
       " {'type': 'A',\n",
       "  'i': '44',\n",
       "  'f': '33',\n",
       "  'rate': 0.001850087889632101,\n",
       "  'energy': '528',\n",
       "  'iconfig': array([1, 2, 4]),\n",
       "  'fconfig': array([2, 1, 3])},\n",
       " {'type': 'A',\n",
       "  'i': '44',\n",
       "  'f': '27',\n",
       "  'rate': 0.002553547336205528,\n",
       "  'energy': '542',\n",
       "  'iconfig': array([1, 2, 4]),\n",
       "  'fconfig': array([2, 2, 2])}]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[p for p in processData if np.all(p['iconfig'] == np.array(configurations[idxP])) and p['type']=='A']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "42004.27766762723"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fluence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### below is the solution for exercise 3b, 3c, 3d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:47.852629Z",
     "start_time": "2022-05-31T09:13:22.522784Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fluences = np.linspace(1e11, 1e13, 21) / (1e-6 * M2AU)**2\n",
    "finalAvgCharges = []\n",
    "probPPs = []\n",
    "probPAs = []\n",
    "fwhmPulse = 25 * FS2AU\n",
    "print(f\"pulse duration = {fwhmPulse*AU2FS: 5.2f} fs \")\n",
    "for fluence in fluences:\n",
    "    sigmaPulse = fwhmPulse / (2. * np.sqrt(2. * np.log(2)))\n",
    "    # pulse Envelope: simple Gaussian with stddev 25fs and t0 = 0fs\n",
    "    pulseEnvelope = lambda t: gaussian(t, sigma=sigmaPulse, x0=0.0 * FS2AU)\n",
    "    res = solveREQ(pIni, processData, configurations, fluence, pulseEnvelope)\n",
    "    # p is population, times is timepoints\n",
    "    p = res['y']\n",
    "    times = res['t']\n",
    "    # compute the average charge for the last time step\n",
    "    finalAvgCharge = nuclearCharge - np.dot(Nelectron, p[:, -1])\n",
    "    # computer the PP and PA probs\n",
    "    dt = times[1] - times[0]\n",
    "    probPP = np.sum(sigmaPP * fluence * pulseEnvelope(times) * p[idxP, :]) * dt\n",
    "    probPA = np.sum(GammaPA * p[idxP, :]) * dt\n",
    "    print(\n",
    "        f\"fluence = {fluence *  (1e-6 * M2AU)**2:10.2e} charge = {finalAvgCharge:5.2f} probPP = {probPP:5.2f} probPA = {probPA:5.2f} \"\n",
    "    )\n",
    "    finalAvgCharges.append(finalAvgCharge)\n",
    "    probPPs.append(probPP)\n",
    "    probPAs.append(probPA)\n",
    "\n",
    "finalAvgCharges = np.array(finalAvgCharges)\n",
    "probPPs = np.array(probPPs)\n",
    "probPAs = np.array(probPAs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:47.976163Z",
     "start_time": "2022-05-31T09:13:47.853982Z"
    }
   },
   "outputs": [],
   "source": [
    "# plot average charge as a function of fluence\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(fluences * (1e-6 * M2AU)**2, finalAvgCharges)\n",
    "ax.set_xlabel(r'fluence (photons per $\\mathrm{\\mu m}^2$)')\n",
    "ax.set_ylabel('final avg. Charge')\n",
    "plt.ylim((0, nuclearCharge))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:13:48.232011Z",
     "start_time": "2022-05-31T09:13:47.978189Z"
    }
   },
   "outputs": [],
   "source": [
    "# plot probability of PP / PA as a function of fluence\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(fluences * (1e-6 * M2AU)**2, probPPs, label='PP')\n",
    "ax.plot(fluences * (1e-6 * M2AU)**2, probPAs, label='PA')\n",
    "\n",
    "ax.set_xlabel(r'fluence (photons per $\\mathrm{\\mu m}^2$)')\n",
    "ax.set_ylabel('probability')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:14:10.479858Z",
     "start_time": "2022-05-31T09:13:48.233187Z"
    }
   },
   "outputs": [],
   "source": [
    "finalAvgCharges = []\n",
    "probPPs = []\n",
    "probPAs = []\n",
    "pulseDurations = np.array([\n",
    "    0.2, 0.3, 0.5, 0.7, 1.0, 1.5, 2.0, 3.0, 5.0, 7.0, 10.0, 15.0, 20.0, 25.0,\n",
    "    30.\n",
    "]) * FS2AU\n",
    "fluence = 2e12 / (1e-6 * M2AU)**2\n",
    "print(f\"fluence = {fluence *  (1e-6 * M2AU)**2:10.2e}\")\n",
    "for fwhmPulse in pulseDurations:\n",
    "    sigmaPulse = fwhmPulse / (2. * np.sqrt(2. * np.log(2)))\n",
    "    # pulse Envelope: simple Gaussian with stddev sigmaPulse and t0 = 0fs\n",
    "    pulseEnvelope = lambda t: gaussian(t, sigma=sigmaPulse, x0=0.0 * FS2AU)\n",
    "    res = solveREQ(pIni, processData, configurations, fluence, pulseEnvelope)\n",
    "    # p is population, times is timepoints\n",
    "    p = res['y']\n",
    "    times = res['t']\n",
    "    finalAvgCharge = nuclearCharge - np.dot(Nelectron, p[:, -1])\n",
    "    # computer the PP and PA probs\n",
    "    dt = times[1] - times[0]\n",
    "    probPP = np.sum(sigmaPP * fluence * pulseEnvelope(times) * p[idxP, :]) * dt\n",
    "    probPA = np.sum(GammaPA * p[idxP, :]) * dt\n",
    "\n",
    "    print(\n",
    "        f\"pulse duration = {fwhmPulse*AU2FS: 5.2f} fs  charge = {finalAvgCharge:5.2f} probPP = {probPP:5.2f} probPA = {probPA:5.2f}\"\n",
    "    )\n",
    "\n",
    "    finalAvgCharges.append(finalAvgCharge)\n",
    "    probPPs.append(probPP)\n",
    "    probPAs.append(probPA)\n",
    "finalAvgCharges = np.array(finalAvgCharges)\n",
    "probPPs = np.array(probPPs)\n",
    "probPAs = np.array(probPAs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:14:10.584887Z",
     "start_time": "2022-05-31T09:14:10.481136Z"
    }
   },
   "outputs": [],
   "source": [
    "# plot Average Charge as a function of pulse duration\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(pulseDurations * AU2FS, finalAvgCharges)\n",
    "ax.set_xlabel('pulse duration FWHM (fs)')\n",
    "ax.set_ylabel('final avg. charge')\n",
    "plt.ylim((0, nuclearCharge))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-05-31T09:14:10.712138Z",
     "start_time": "2022-05-31T09:14:10.586075Z"
    }
   },
   "outputs": [],
   "source": [
    "# plot probability of PP / PA as a function of pulse duration\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(pulseDurations * AU2FS, probPPs, label='PP')\n",
    "ax.plot(pulseDurations * AU2FS, probPAs, label='PA')\n",
    "ax.set_xlabel('pulse duration FWHM (fs)')\n",
    "ax.set_ylabel('probability')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}