{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "93a49da0",
   "metadata": {},
   "source": [
    "# Differential cross-section uncertainty propagation with KDUQ\n",
    "\n",
    "KDUQ is an uncertainty-quantified global optical potential. This notebook shows how to sample that posterior information, run the solver repeatedly, and summarize the resulting uncertainty in differential cross sections.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9fbc70e4-6ac7-4a99-8a84-3197e65a8666",
   "metadata": {
    "editable": true,
    "scrolled": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import stuff for nice plotting\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from tqdm import tqdm\n",
    "\n",
    "import jitr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a13cc786-67c3-416b-a632-1524e9496f70",
   "metadata": {},
   "outputs": [],
   "source": [
    "from jitr.optical_potentials import kduq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "943628b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "#  elastic reaction\n",
    "target = (54, 26)\n",
    "proton = (1, 1)\n",
    "neutron = (1, 0)\n",
    "projectile = proton\n",
    "\n",
    "reaction = jitr.reactions.ElasticReaction(\n",
    "    target=target,\n",
    "    projectile=projectile,\n",
    ")\n",
    "\n",
    "# energy\n",
    "Elab = 35\n",
    "kinematics = reaction.kinematics(Elab)\n",
    "\n",
    "# for plotting differential xs\n",
    "angles = np.linspace(0.1, np.pi, 180)\n",
    "\n",
    "# Lagrange Mesh\n",
    "core_solver = jitr.rmatrix.Solver(40)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24060cde",
   "metadata": {},
   "source": [
    "## Set up the solver for differential cross sections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "59075875",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get kinematics and parameters for this experiment\n",
    "\n",
    "a = jitr.utils.interaction_range(target[0]) * kinematics.k + np.pi * 2\n",
    "N = jitr.utils.suggested_basis_size(a)\n",
    "assert N < core_solver.kernel.quadrature.nbasis\n",
    "channel_radius_fm = a / kinematics.k\n",
    "\n",
    "# build solver\n",
    "solver = jitr.xs.elastic.DifferentialWorkspace.build_from_system(\n",
    "    reaction=reaction,\n",
    "    kinematics=kinematics,\n",
    "    channel_radius_fm=channel_radius_fm,\n",
    "    solver=core_solver,\n",
    "    lmax=50,\n",
    "    angles=angles,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87220ced-3aa9-4dff-bc60-411c7e945fb4",
   "metadata": {},
   "source": [
    "## Sample the KDUQ posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f6c3b321-1561-4c30-a5c6-dfb24cfeb66b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function get_samples in module jitr.optical_potentials.kduq:\n",
      "\n",
      "get_samples(projectile: tuple[int, int], posterior: str = 'federal') -> numpy.ndarray\n",
      "    Get the posterior samples for the given projectile (neutron or\n",
      "    proton) from the KDUQ Federal or Democratic posteriors.\n",
      "\n",
      "    See [Pruitt, et al., 2023]\n",
      "    (https://journals.aps.org/prc/pdf/10.1103/PhysRevC.107.014602) for\n",
      "    details on the KDUQ posteriors.\n",
      "\n",
      "    :param projectile: tuple A tuple representing the projectile, with format (Ap,\n",
      "                       Zp), where Ap is the mass number and Zp is the atomic number.\n",
      "                       Must be either (1, 0) for neutron or (1, 1) for proton.\n",
      "    :param posterior: str Which KDUQ posterior to return samples from. Must be\n",
      "                      either \"federal\" or \"democratic\". Defaults to \"federal\".\n",
      "    :returns: An array of shape (NUM_POSTERIOR_SAMPLES, num_params) containing the\n",
      "              posterior samples for the given projectile, where num_params is the\n",
      "              number of parameters in the Koning-Delaroche potential, and the\n",
      "              parameters are ordered according to the order set by the\n",
      "              get_param_names function.\n",
      "    :rtype: np.ndarray\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(kduq.get_samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ffa64ad7-c663-40e6-a6c4-b340122f705b",
   "metadata": {},
   "outputs": [],
   "source": [
    "kduq_samples = kduq.get_samples(projectile)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08fda017-a059-4da4-8dfd-dbd0e7c85055",
   "metadata": {},
   "source": [
    "### Set up the potential\n",
    "\n",
    "Let's take a look at the functionality built into jitr:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f1f685f7-17b7-43bd-9de6-c28fb10e65fa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on KDUQ in module jitr.optical_potentials.kduq object:\n",
      "\n",
      "class KDUQ(jitr.optical_potentials.omp.SingleChannelOpticalModel)\n",
      " |  KDUQ(projectile: tuple)\n",
      " |\n",
      " |  Koning-Delaroche Uncertainty Quantification (KDUQ) optical\n",
      " |  potential model.\n",
      " |\n",
      " |  Method resolution order:\n",
      " |      KDUQ\n",
      " |      jitr.optical_potentials.omp.SingleChannelOpticalModel\n",
      " |      builtins.object\n",
      " |\n",
      " |  Methods defined here:\n",
      " |\n",
      " |  __init__(self, projectile: tuple)\n",
      " |      Initialize self.  See help(type(self)) for accurate signature.\n",
      " |\n",
      " |  evaluate(self, rgrid: float | numpy.ndarray, reaction: jitr.reactions.reaction.Reaction, kinematics: jitr.utils.kinematics.ChannelKinematics, *params: float) -> tuple[complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]], complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]], float | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]]\n",
      " |      Evaluate the KDUQ central, spin-orbit, and Coulomb terms.\n",
      " |\n",
      " |  ----------------------------------------------------------------------\n",
      " |  Methods inherited from jitr.optical_potentials.omp.SingleChannelOpticalModel:\n",
      " |\n",
      " |  __call__(self, rgrid: 'ArrayOrScalar', reaction: 'reaction.Reaction', kinematics: 'kinematics.ChannelKinematics', *params: 'float') -> 'tuple[PotentialArray, PotentialArray, PotentialArray | ArrayOrScalar]'\n",
      " |      Call self as a function.\n",
      " |\n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors inherited from jitr.optical_potentials.omp.SingleChannelOpticalModel:\n",
      " |\n",
      " |  __dict__\n",
      " |      dictionary for instance variables\n",
      " |\n",
      " |  __weakref__\n",
      " |      list of weak references to the object\n",
      "\n"
     ]
    }
   ],
   "source": [
    "omp = kduq.KDUQ(projectile)\n",
    "help(omp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b4b6f520-5169-4349-b426-333edd436676",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function central in module jitr.optical_potentials.kduq:\n",
      "\n",
      "central(r: float | numpy.ndarray, Vv: float, Rv: float, av: float, Wv: float, Rwv: float, awv: float, Wd: float, Rd: float, ad: float) -> complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]]\n",
      "    Koning-Delaroche central terms at a given energy.\n",
      "\n",
      "    This matches Eq. (7) in Koning and Delaroche (2003).\n",
      "\n",
      "    :param r: float or np.ndarray The radius at which to evaluate the potential.\n",
      "    :param Vv: float The real central depth.\n",
      "    :param Rv: float The real central radius parameter.\n",
      "    :param av: float The real central diffuseness parameter.\n",
      "    :param Wv: float The imaginary volume depth.\n",
      "    :param Rwv: float The imaginary volume radius parameter.\n",
      "    :param awv: float The imaginary volume diffuseness parameter.\n",
      "    :param Wd: float The imaginary surface depth.\n",
      "    :param Rd: float The imaginary surface radius parameter.\n",
      "    :param ad: float The imaginary surface diffuseness parameter.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(jitr.optical_potentials.kduq.central)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "020770dd-b6c7-40f6-9447-0eaf46bc089f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function spin_orbit in module jitr.optical_potentials.kduq:\n",
      "\n",
      "spin_orbit(r: float | numpy.ndarray, Vso: float, Rso: float, aso: float, Wso: float, Rwso: float, awso: float) -> complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]]\n",
      "    Koning-Delaroche spin-orbit terms at a given energy.\n",
      "\n",
      "    This matches Eq. (7) in Koning and Delaroche (2003).\n",
      "\n",
      "    :param r: float or np.ndarray The radius at which to evaluate the potential.\n",
      "    :param Vso: float The real spin-orbit depth.\n",
      "    :param Rso: float The real spin-orbit radius parameter.\n",
      "    :param aso: float The real spin-orbit diffuseness parameter.\n",
      "    :param Wso: float The imaginary spin-orbit depth.\n",
      "    :param Rwso: float The imaginary spin-orbit radius parameter.\n",
      "    :param awso: float The imaginary spin-orbit diffuseness parameter.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(jitr.optical_potentials.kduq.spin_orbit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7d53236a-9334-4286-ab9d-805ef7a8d1ea",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function coulomb_charged_sphere in module jitr.optical_potentials.potential_forms:\n",
      "\n",
      "coulomb_charged_sphere(r: 'ArrayOrScalar', zz: 'float', r_c: 'float') -> 'ArrayOrScalar'\n",
      "    Return the Coulomb potential of a uniformly charged sphere.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(jitr.optical_potentials.potential_forms.coulomb_charged_sphere)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1304a608-88b1-49d3-a64d-1457e3bd24b7",
   "metadata": {},
   "source": [
    "## Run the calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db00328f-fa98-4c1a-a4fe-c49931baf401",
   "metadata": {},
   "source": [
    "This should demonstrate how the `KDUQ` class (`omp`) is built to interface with `solver`, which is an instance of `jitr.xs.elastic.DifferentialWorkspace`. This is how calculations are done in jitr."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "2398265d-eae5-467c-969b-41441cfa7190",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 95%|██████████████████████████████████████████████████████████████████████████▍   | 397/416 [00:11<00:00, 229.31it/s]/home/kyle/umich/jitr/src/jitr/optical_potentials/kduq.py:467: RuntimeWarning: overflow encountered in exp\n",
      "  d2 = d2_0 + d2_A / (1 + np.exp((A - d2_A3) / d2_A2))\n",
      "100%|███████████████████████████████████████████████████████████████████████████████| 416/416 [00:11<00:00, 36.63it/s]\n"
     ]
    }
   ],
   "source": [
    "kduq_xs = np.zeros((len(angles), kduq.NUM_POSTERIOR_SAMPLES))\n",
    "rgrid = solver.radial_grid()\n",
    "\n",
    "for j, sample in enumerate(tqdm(kduq_samples)):\n",
    "    central_term, spin_orbit_term, coulomb_term = omp(\n",
    "        rgrid, reaction, kinematics, *sample\n",
    "    )\n",
    "    xs = solver.xs(central_term, spin_orbit_term, coulomb_term)\n",
    "    kduq_xs[:, j] = xs.dsdo / solver.rutherford"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "164af22c-5f75-476d-a15d-4c81ff2690d6",
   "metadata": {},
   "source": [
    "## Compute intervals and plot results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "82fae416-346a-4593-86f2-8ff7152df82c",
   "metadata": {},
   "outputs": [],
   "source": [
    "kduq_pred_post = np.percentile(kduq_xs, [16, 84], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "07abb330-463d-47d7-bb1e-826db06dc9a0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$\\\\frac{d \\\\sigma}{d\\\\Omega} / \\\\frac{d \\\\sigma_{R}}{d\\\\Omega} $')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiQAAAF3CAYAAACVAVC0AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAe7NJREFUeJzt3Xl4W/WVP/639l2y5d2xkzgLgewQSDAdIBQKTRlKusHwo09oC8zQSVsoHcqk05n2285gWsrSb8sAbaGhX6C0tAU67AEaKCUECMmQEEjiYGLHa7xpl+6VdH9/HF1LsiVbkrVd+byeR09irVfrPfd8zud8VJIkSWCMMcYYKyF1qTeAMcYYY4wDEsYYY4yVHAckjDHGGCs5DkgYY4wxVnIckDDGGGOs5DggYYwxxljJcUDCGGOMsZLjgIQxxhhjJact9QaUu2g0ir6+PthsNqhUqlJvDmOMMaYYkiTB4/GgubkZavX0ORAOSGbQ19eH1tbWUm8GY4wxplg9PT1oaWmZ9jockMzAZrMBoBfTbreXeGsYY4wx5XC73WhtbZ3Yl06HA5IZyMM0drudAxLGGGMsB5mUPCiqqPWee+7B6tWrJ4KD9vZ2PPvss2mvv337dqhUqqST0Wgs4hYzxhhjLBOKypC0tLTg1ltvxdKlSyFJEh588EFceuml2Lt3L1asWJHyNna7HYcOHZr4mwtTGWOMsfKjqIDkkksuSfr7v/7rv3DPPffgjTfeSBuQqFQqNDY2FmPzGGOMMZYjRQ3ZJIpEInj00Ufh8/nQ3t6e9nperxcLFixAa2srLr30Urz33nvT3m8oFILb7U46McYYY6ywFBeQ7N+/H1arFQaDAddddx0ef/xxLF++POV1ly1bhgceeABPPvkkHnroIUSjUZx11lk4fvx42vvv6OiAw+GYOPGUX8YYY6zwVJIkSaXeiGwIgoDu7m64XC784Q9/wK9+9Su88soraYOSRKIo4pRTTsEVV1yBH/7whymvEwqFEAqFJv6Wpyy5XC6eZcMYY4xlwe12w+FwZLQPVVQNCQDo9XosWbIEALBu3Tq89dZb+OlPf4r77rtvxtvqdDqceuqp6OzsTHsdg8EAg8GQt+1ljDHG2MwUN2QzWTQaTcpoTCcSiWD//v1oamoq8FYxxhhjLBuKypBs27YNmzZtwvz58+HxePDII49g586deP755wEAW7Zswbx589DR0QEA+MEPfoAzzzwTS5Yswfj4OG677TYcO3YM11xzTSmfBmOMMcYmUVRAMjQ0hC1btqC/vx8OhwOrV6/G888/j0984hMAgO7u7qTFe8bGxnDttddiYGAA1dXVWLduHV5//fWM6k0YY4wxVjyKK2ottmwKcliZikaB4eOAaxgIBYGmhUBNc6m3ijHGKl5FF7UyljXPKHD4bWDwGOBsBhAFquoBDX/8GWOsXCi+qJWxGY0PARERcNQBzkbAPQKMnyj1VjHGGEvAAQmrbEIQGOkHjLFUoUZLp6FjNJTDGGOsLHBAwiqbewTwuwGTNX6erYayJq7h0m0XY4yxJDyIziqXJAEneiggkY0PAXoj4BoBBrqA6vrSbR9jjLEJnCFhlcvvAfqP0v+Dvvi/fjegUtGwTShQuu1jjDE2gTMkrHL5xgGDGWiZT7UkA11AYxtlSCSJpgL7XIDBVOotZYyxOY8zJKxyecYAjS71ZSoVABXgHS/mFjHGGEuDAxJWmSIR6j+inyb7YTQD44N0XcYYYyXFAQmrTEEfEPLTkE06RisQ8AIBT/G2izHGWEockLDKFPQCYQHQ6dNfR6cHRCF5Fg5jjLGS4ICEVSafG1Bl8PHW6oFx7kfCGGOlxrNsWOWRJGBsCICKZtcAU/+VabTUm0QIAXpDUTeTMcZYHAckrPKEAsBoH9WHTB6OGe1P/jsaBcQADfFwQMIYYyXDAQmrPEEfoNYALcsAdWzYRghSMOJsoj4kiU70cIM0xhgrMQ5IWOUJeAGtlqb1TqY3Tg1IDCa6DWOMsZLholZWeTyjgDaL4RedkRukMcZYiXFAwipLNEp1I7psAhIDDfOExcJtF2OMsWlxQMIqixii3iLZBiRikOtIGGOshDggYZVFCFJQok2zhk0qOj1lRwQOSBhjrFQ4IGGVRQgCUoT6i2SLMySMMVYyPMuGVRYhSNN4RSH5fLk+ZKArdfbEdQKobQEaFxZ8ExljjE3FAQmrLD4XYLIDZnvy+UKQCleNlqnTfgHq7jo+SP+qVMXZVsYYYxM4IGGVJegFapoBe03y+UKQsiNV9akDEpOVursKQepLwhhjrKi4hoRVDlGggCKbGTYynQEI+LiOhDHGSoQDElY5JmbY6LO/rVpDU395pg1jjJUEBySscghBKl7NZspvonAECPrzu02MMcYywgEJqxxikP7NtShVowZ843nbHMYYY5njgIRVDr93djNkNDrA66KZNowxxoqKZ9mwyuE+Afg9wPjQ1MuEWPZk8CPAZEszrKMCvKM07KPLoQ6FMcZYzhSVIbnnnnuwevVq2O122O12tLe349lnn532No899hhOPvlkGI1GrFq1Cs8880yRtpYVVSQCeMapqNXvnnpynQBcQ1S8OtRN/UomX0ejpVk2YqjUz4YxxuYcRWVIWlpacOutt2Lp0qWQJAkPPvggLr30UuzduxcrVqyYcv3XX38dV1xxBTo6OvD3f//3eOSRR7B582a88847WLlyZQmeASsYMURZkJZlU/uIhAKUEVnxd0B1PfDhu0AkPLVXSe8RICLS/Zhtxdt2xhhjUEmSsgfMnU4nbrvtNlx99dVTLrv88svh8/nw1FNPTZx35plnYu3atbj33nszun+32w2HwwGXywW73T7zDVhpeMaA158AFqygLEiikT6gZh5w0mn091APcGQPUN2YPHTT1wmE/MBpnwDqWoq26YwxVqmy2YcqasgmUSQSwaOPPgqfz4f29vaU19m1axcuuOCCpPMuuugi7Nq1qxibyIopLFBjtMnBSDQCSFGgtjl+Xk0znVwpak2iEg/ZMMZYCShqyAYA9u/fj/b2dgSDQVitVjz++ONYvnx5yusODAygoaEh6byGhgYMDAykvf9QKIRQKL5Dcrvd+dlwVlhiiIZbJvO5AYsjeXhGowEa24DRARq6SVwZWKMBAt7Cby9jjLEkisuQLFu2DPv27cPu3bvx1a9+FVdddRUOHjyYt/vv6OiAw+GYOLW2tubtvlkBCSHKbkwW8NAqvpNn1VirAJOFFtxLpNXRbRhjjBWV4gISvV6PJUuWYN26dejo6MCaNWvw05/+NOV1GxsbMTg4mHTe4OAgGhsb097/tm3b4HK5Jk49PT153X5WIH4PNTZLJAQBvQGoqpt6fa0OcNSmCEj08Y6vjDHGikZxAclk0Wg0aYglUXt7O1566aWk83bs2JG25gQADAbDxLRi+cQUwDsGSKBgQj6NDwFWJ2BO8x7aaqiIVb5+WKSak6Cf60gYY6zIFFVDsm3bNmzatAnz58+Hx+PBI488gp07d+L5558HAGzZsgXz5s1DR0cHAOD666/Hueeei9tvvx0XX3wxHn30Ubz99tv4xS9+UcqnwfItEgYCsT4iA13x810ngMVr0ndvNdtpqKfnULwRmhii4EYIAiZr4bedMcYYAIUFJENDQ9iyZQv6+/vhcDiwevVqPP/88/jEJz4BAOju7oZaHU/6nHXWWXjkkUfw3e9+F9/5znewdOlSPPHEE9yDpNKIAgUl6lixKkDt33UGwDhNUGGyAvULKJixOSmYMdkAvZEzJIwxVmSK70NSaNyHRAG848CbzwKQqDEaQM3Qgl5g5dlUvJpOzyGg5wMqfO3rpKyJGAIWrQaaFhVj6xljrGLNiT4kjE0QQ0BkUg8SIQAYzIDRPP1trdWUTYlG4+dptDz1lzHGiowDEqZ8opAcUACUIbHXzLz6r8VOQUvIHz9Pq+eAhDHGiowDEqZ8YWFqD5JohBqizURvpDoTMRg/T6ujDEsknN/tZIwxlhYHJEz5Al7qsCqTu68ap6kdSWRz0qwamVYPhMOUeWGMMVYUHJAw5Qv6ktu/h+T6kQwDEqNlag1JRKTMC2OMsaJQ1LRfxqaIRADXMHVqNVlppoxnFGhbHe8tMhOjmW6j0QHhEcDvBvyuuZEhiUap70okQn/bncnBHWOMFQn/8jBlCwu0mq+8EzXbabpvTVPm92EwUwGsa5hqSsx2qiGp9AxJJAwcex/46ABlhFQqYNl6oHVZqbeMMTYHcUDClC0sUDBS3UAzZarqaRqwYYbpvon0RqojGR+iwKSqnu63kgOSsAh0vw/0HwXq59Nr4PcAAx/Ra2mtKvUWMsbmGK4hYcomhOJFrAD9X6sDDKbM70Oloh1wYgCi1iRPBa40fUdpeKuqgYIRADDbaLZR39Gp06gZY6zAOCBhyjY5iyEEaZaMzpjd/UyeIqzRAQFf6usqnSgAxw8BZge110/kqKMW+mODqW/LGGMFwgEJU7bJAUlYoJ1spgWtMr2JMivySgoaLWVIKjFT4B4GRvtTr4KsM1B2qPv9+GvBGGNFwAEJUzYhmNwyXhQo2zFTh9bJjGZaWE8OcLS62NRfMX/bWg4kCRjqplk16jRff4sDGDpGawQxxliRcEDClG1yD5KISLUQ2dKbAJM9vsqvRkfBSKUVtvrdwIme6V8jnQHwjNGsI8YYKxKeZcOUzeuiI/lorI+GEIgXaWZDowEctcAHR2koQ+7PUWkBydggBSSOeppVlIoQpNew9wjQuJCyRYwxVmAckDDlioSp54gQoELWaJRO2Ra0yqrqAClCWQQACHgqqzlaJAwcPxx/feTnOVlYpPV9hrqpYVx1Q/G2kTE2Z3FAwpRLFGioZd5JFDyMDwE183LLkAA0VGGtBpqX0N8neiorQ+JzUw1JXSsV7MrPczIhSDNt7E5gbIADEsZYUXANCVOuiAhEw/EhhbAI6A25ByR6IxXIyrNLVCpaF6dS+N3JPVtmYnbQEE8lvQaMsbLFAQlTrnBsFow8yyYs0EyZdLNHZqIzUEM1eWaNRkdFs5VibCC7YM1oAYJ+wOcq3DYxxlgMByRMueThFHmKbySSurdGpnQGqq9InPpbKQGJHFhkugIyEHtdVTz9lzFWFByQMOWa3CNEo8muZfxkWh3VpMgBiTz1NxLO/T7Lhc9FQy/ZrPEDUAAzNhhfDZgxxgqEAxKmXGIIQCw7Eo3STJtc60eA2Jo2jqnN0YTQrDe15Dyj9PyybhhnoZlMAU9htosxxmI4IGHKFfRTUasQpB3mbApaZRYHEPDSfUbCNBtF6TNtImFgdABQa+l5CUHK/Mj/T3UC6F8pSqsAp5sizBhjecLTfplyeUaAkX4aigj6qLHZ5MXisqU30c5bHg7yu5QfkPg9wNBHgEpDzd5kA13T3260n/51j9AU6Pr5BdtExhjjgIQpUzQKSACaFlFW4/hhmmEz266iOgNQ0ww4m2h4Y/i48pujBX302tTMo7/Hh+i8xrbU1xeCFIw4myjj5KgDfJ54B1fGGCsADkiYMoVFGk4wWeP9QyxVs79fvYEW2tNoKbjRGaiORMn8bnoecjChj80kmim40BvppNUDI71UGMsBCWOsQLiGhClTRKST3OQrGgZsVbO/X52RApGJmTZaZU/9lSQqaNXPYvaRWg1AUvbrwBgrexyQMGUSBSCc0HVUpcrP0btOTxmBpOZo/tnfb6mIISrMne1rozVQLQljjBUIByRMmSKxIRu5S6tKM/uCVoACG5MtNqUYFPCE/PF28koT9FHR72wDEoOJ1sKZ3PuFMcbyhAMSpkyJhaaR2Ho2Wn1+7ttkje94tTq6f6UWtgZjwZQcuOVKbwLEIA/bMMYKhgMSpkyJhaaRMA215CNDAsS6vcYyIppYczSlTv31uwBVHr7m2ljXWg5IGGMFwrNsmDL5vdQbQwxRIzO7k4KSfNAZaCE6MUTDQkGFNkeTJGC4j56LkLBir5z96etMfTv58oGu5GnU40NAyzKgdl5htpcxNqcpKkPS0dGBM844AzabDfX19di8eTMOHTo07W22b98OlUqVdDIaeeqi4oV8gLWKFtMzmACDJV7gOls6A00hNpjpX41WmUM2QpDqX+y19DrJJ1nieYkneQE+oyX5fJuTghLGGCsARWVIXnnlFWzduhVnnHEGwuEwvvOd7+DCCy/EwYMHYbGkX8XUbrcnBS6qbNfzYOUnLALVjbSTlKK0080XnQGwVtMO2WCiTIkSMyRBHxXp1rXGpu4m8LuBqvrUtxOClB2pqk8uhjWYgXCIG6QxxgpCUQHJc889l/T39u3bUV9fjz179uCcc85JezuVSoXGxsZCbx4rlkg4uQdJOAyYbfm7f60+3ovEYAKgis+6URJ5dtDkYCRXBhM1Rwv6OCBhjOWdooZsJnO5XAAAp9M57fW8Xi8WLFiA1tZWXHrppXjvvffSXjcUCsHtdiedWJmJhOk0MUQj5XcHqdHQrBK5cFarBQIKLOb0u2c/uyaRWhOrqVHga8EYK3uKDUii0ShuuOEGfOxjH8PKlSvTXm/ZsmV44IEH8OSTT+Khhx5CNBrFWWedhePHj6e8fkdHBxwOx8SptbW1UE+B5SosAtEIrV4ry1dBq8xkSW6OFlLYTliSAO94/mYeydQaypIwxlieKTYg2bp1Kw4cOIBHH3102uu1t7djy5YtWLt2Lc4991z86U9/Ql1dHe67776U19+2bRtcLtfEqaenpxCbz2YjLMa6tGooU6LOU1O0RCYr3TdAwzdiCIhE8vsYhSQK1BAt36+L3gh4xpTbKI4xVrYUVUMi+9rXvoannnoKr776KlpaWrK6rU6nw6mnnorOztRTHg0GAwyGPP+Is/yKhONdWkOBWM1HnjMkifen0dEwRVgANLNYE6aYhAAFUYmzavJBb6TaFCEYq69hjLH8UFRAIkkSvv71r+Pxxx/Hzp070daWZvn0aUQiEezfvx+f+tSnCrCFrCgiIuAbp6P/gIeGV/I9ZCPf99ggDQ/JAYlSdsKhAK09k6qGRF6TJt0UXiEYv3xybU40QkM2Ib9yXgvGmCIoKiDZunUrHnnkETz55JOw2WwYGBgAADgcDphM9OO4ZcsWzJs3Dx0dHQCAH/zgBzjzzDOxZMkSjI+P47bbbsOxY8dwzTXXlOx5sFkKi4DfQ4Wn3nHAXpO/HiQynT5WhzFG9+13KasXSchPRa3TFfv60xRsy7UzchA2mc9FzeLsNbPfTsYYi1FUQHLPPfcAADZu3Jh0/q9//Wt86UtfAgB0d3dDnTDNcWxsDNdeey0GBgZQXV2NdevW4fXXX8fy5cuLtdks38IiUDcfqGsBRgfy24NEpjMAzkZAZ6R+JMPHldWLxD0CNC2iXi2TjQ9RMNK8JPVt5T4kjW2pAxq9kQJCxhjLI0UFJFIGhXQ7d+5M+vvOO+/EnXfeWaAtYiUhBOIZkUiYGnblm1ZPtSOJq9sqJUMSjVLhqa5AvUL0RmB8AAAH9Yyx/FHsLBs2h4UCCbUREqAvQBGyWh3rTCoHJCrl9N8IBYCgN/8zbGQ6I+AeAwQFNotjjJUtDkiY8ojB5JoRjS79dWfDbIsP02h1VOSqBEKgMFN+ZXojBTwhf2HunzE2J3FAwpQlsUurJNFaLfmeYSPTm2h6MUBBj1K6tYYCNBMpXy3jJ9NoY1N/AzNflzHGMsQBCVOWsBgLSDTxwKRQGZLEoSCNlhaWS6wpKVd+DxAtcOOysAgEvIV9DMbYnKKoolbGEAnTLBCtPt4vo1AZEq2eMg2hAPXfCPlpR6wtUACUL2ODAKT46zOZEKTnMd3lif+mMz4EtJyU82YyxlgiDkiYsoRF2hH6XBQoVDcWLkOiMwDBANB7GFCpASkSqykpwKyefAmLgH+cAreBrumvO9Plo/3pL9MZANeJSYscMsZY7viXhClLRAQsVUBNM+AZBaxVNHxTCDoDUFVPtRhGC3D8/fIfshGCgM9Dr0ttmmUVxodoxlBjGwVYiR1b9UbAaKVgxNmUvrHa8cOx2Tx+wJLn9vSMsTmJAxKmLJEwDZnojbFAwVq4x9LqAKMJEMXY42nLvzmaEKQ1bBy16YMJvZGeh84AjA8C85bS9SWJAg3vOP1fb0x/HwYT1dQIAQ5IGGN5wUWtTFkSMxSRMGAs4PCJSkW9SCKxx5SiysiQRATa9pm4RwBrNdCyDKifDzQsAJaeRtkVz8j0t1Wp6PUP8Uwbxlh+cEDClCUsAIjtbCWpcL02ZEZLPAiRQEMU5UwIUKfWmYghynC0LqPFCWX2GmDRaurBIs7Q+EylSb8eDmOMZYkDEqYswqSmaIWaYSOb0oukzHfAPldmRaa+caBhIdWJTGavAZqX0sKC09Hp6fEyWNKBMcZmwgEJUxa5bbwkAVAVboaNLDHg0WjKe1G5aJR6g8z0moghKgyua009tKNWU11JNErDMuloDZSRmSmTwhhjGeCAhCmL3DY+GqEAodAZEq0egIoCILlDaSRS2MfMlRCkBQBn6pPicwHNi6lWJJ2qOip09Y2nv45OT+vZcB0JYywPeJYNU45ImFax9YxSlkRnLE6GxDceW7slAJisVMeiMRX2cXMhBIGRXpoN5B5JXd8RjdJzaWybvvBVo6Xr7H6aXvPJ1w2LQNhF6/sIAQDVeX0qjLG5hwMSphxycam1Oj7DphgZErM9/vhCgP41lGFAIgZpuEVew8acYjqudxyomw846ma+v6p6oK6Fnq/RknyZOzYLR2fgDAljLC84IGHKIfcgsTrpKF9nKHyXUK2e+myo1PR4I73l24tECFJtiBiiYKSqPsV1AsCC5Zm1vzeagfoFwImeqffld9NjSFEKchhjbJa4hoQpR1ik+g2NtvA9SGQaDQUiYTHee6NcAxK/Z/ohLCFIjc4yyY7IquvpOaebSaMzUnCSyVRjxhibBgckTDkiIh2Rq9WxYZMirSljsMSbo0Uj5dkcTZIoMJhuCCvoo6yGKYvuttZq6lMS9KW+XGegjAwP2zDGZokDEqYckXC8uFKKUo+QYjBZgHBs+ms0Up7TXEUhvgpyOiE/UNUQrzHJhMEE2OvSN0CTAxKBAxLG2OxwQMKUY/LQQSZ1EPmgM4DatCI29bcMd75iML4+TSrRCAVz0031TWe6YRu1moZryvE1YYwpCgckTDnEhLbxwPTZgHxKfBy1Nv3wRSkJQRpKShekBf00UybVzJuZzDRso1ZzhoQxNms8y4YpR9BL7cw1Ghp+KFqGRE87Y7+bsgQ+F2UFshn6KDQxRK+NXIDrnrQ43tgg0LYK0Oew9o/RTLN3ut+Pz7ZJfAy/e+rjMcZYlsroF5WxGQR8dCTuHaPAoFgBiZwhiYSpIVvQW34zbfxemmWTWOvhd8dPQV/qacCZqm6gJmjy/SU+hhACXCM804YxNiucIWHKERaA5iWxQspg4bu0yrR62pn73BQQGcyUIdAbi/P4mQj5aP2ZqnqgrzO5D4kQpNV7LY7c799ko54kVfVUR5P4GKEABTxCsDhTsRljFYkzJEwZIpFYJ1JNLFOhLe6QjVZHU381mvLrRSIvqpeupibkB4zW7Kb7Tma2USAW8k+9TKePzbQJ5n7/jLE5jwMSpgyRMGUl5KZoOn3xajhUKppiHAlTx1ZJKq9eJGIoNsMmXUASAOw1s3u9tDoqbg2mCEjUGkCKcGErY2xWOCBhyhARqReInCEpVg8SmclCjysrpwzJxCq/aQKSaCS32TWTOWriDeKmUHGGhDE2KxyQMGWQh0nkDEmxF7czmKkZGwBAVX4Zkmgk9bo+ciHu5MXxcmGyUYYoMTCTaXVUY8MYYznigIQpQ1iko3O1JtaltcgFpVp9vEusRlNevUjEUPq1ZoQgBW95CUisdD+p6ki0BiDgTr8djDE2Aw5ImDLIGRIxRKdizbCRafVUJxEWKRvh9xb38acT8NLrIwTjDdLk//vGKbuTS/+RyXR66vTqGU1+DCFIgUjQV55t9RljiqCogKSjowNnnHEGbDYb6uvrsXnzZhw6dGjG2z322GM4+eSTYTQasWrVKjzzzDNF2FqWV5Ew9f8Y6ALcw8WbYSPT6QHEVvsd7gO8o+WTDXAPA6MD9NoMdNF5QR/9v68TsDnz91iOWmCoO/kxBrqAkV5gtJ/rSBhjOVNUQPLKK69g69ateOONN7Bjxw6IoogLL7wQPl/69Pnrr7+OK664AldffTX27t2LzZs3Y/PmzThw4EARt5zNWkQENHqgYSH1vih2QKLVU2FoJAw0L6Zpx+VQRxKJAKJI/Vka2+gE0NBKYxtQ1wpYZ9F/ZDKTle4zGo0/RmMb9UAx2TggYYzlTFGN0Z577rmkv7dv3476+nrs2bMH55xzTsrb/PSnP8UnP/lJ3HTTTQCAH/7wh9ixYwd+/vOf49577y34NrM8CQVoZo1GS0MQxR6y0emp6ZfPTTtln2v6qbbFIhe0mm3xuhqtjv4vF7Ma8tiszGgBLHYaJtIbk2t5tHoOSBhjOVNUhmQyl8sFAHA606ekd+3ahQsuuCDpvIsuugi7du1Kef1QKAS32510YmUg4CtNUzSZRhvvRaKJZUfEMpj6KwQpKEkVGAkBQGfMT0GrTG+kJmupakXUGmovzxhjOVBsQBKNRnHDDTfgYx/7GFauXJn2egMDA2hoaEg6r6GhAQMDAymv39HRAYfDMXFqbW3N63azHIV88S6pmhIEJAANSURj02ijkWl6chSRGKJZR2rN1MuEIGUz8v1a2Zyp+7DoDJQ5YoyxHCg2INm6dSsOHDiARx99NK/3u23bNrhcrolTT09PXu+f5UCSaKqpnCHRalP33Cg0uYYEoCnA5ZIhQZriWjFE3VXzzWxL/frrDPEmbYwxliVF1ZDIvva1r+Gpp57Cq6++ipaWlmmv29jYiMHBwaTzBgcH0djYmPL6BoMBBkMepkiy/ImEY1N9Y03RzLbSbIdclwFQQFIOGZKQj4aw0snncE3ifRotU4t6tXqaCSUES19bwxhTHEVlSCRJwte+9jU8/vjjePnll9HW1jbjbdrb2/HSSy8lnbdjxw60t7cXajNZvoVFmkmi1tCQSbG7tMq0+oSARJ26QVix+dMsqhcWaaimEKvvGi20crA4qYBVp4/1J+E1bRhj2StYhkSj0SASieT1Prdu3YpHHnkETz75JGw220QdiMPhgMlEO6ktW7Zg3rx56OjoAABcf/31OPfcc3H77bfj4osvxqOPPoq3334bv/jFL/K6bayAImFA8NOOzucCWk4qzXbIC/r1Hga844UZDslGJAKMDwKuYWqAJguLtH3VDfmdYSPT6qgfyeE9U3uxuE7wTBvGWE4KliGRCtA06p577oHL5cLGjRvR1NQ0cfrd7343cZ3u7m709/dP/H3WWWfhkUcewS9+8QusWbMGf/jDH/DEE09MWwjLykxEpCJKtZrqOHQlGlLTGQCdCdCbKRgRAqVtjiYEKVizOel1kU8A1XE46gpX/Otsii/al3gyWrmwlTGWk1llSCRJgkpe32OSdOfP9vFmsnPnzinnfeELX8AXvvCFvG8PK5JImHb8jjqqJSnFDBuAhkaMFmqfjioKlCLh0m2P3ELf2Zg8y8Yf65VSVV+4xzZaALuT3pPE77paTT1KGGMsSzkHJJ/73Ofg9/tx5ZVXYv/+/ejs7MQf//jHfG4bYyQSps6gAABV8ZuiybQ6wGik7dEbKUMiCqUNSNJN+VWrAFMBClplRgtgiPUjSWyOpjNSbY1YBk3jGGOKkvOQTV1dHZ599lncd999+Kd/+qe0wcgbb7yByy67DGvWrMG6detw/fXXo7e3l4tKWebkBe0iEepFUqoAQKUCTA7aHrk5Wiln2qRbyC4SpsBAX8DiX4OZgpLJ9SJaPW0X15EwxrKUc4bk7bffxve//304nc6U02QlScLDDz+Mn/3sZ7jllluwatUqeDwePPXUUzjvvPMwPj4+m+1mc0lYpExANJaZKFWGBKDZJd4xCkiikdL23Ah6U2dHRIGKTgtR0CpTqwFHDTA+lHy+VhdfCdhiL9zjM8YqTtYByR/+8Ad8/vOfx9tvv43Ozk6sWLECd999N7q7u/HQQw8lXbejowOvvfYaqqqqAFBW5Rvf+AY2bdqEL37xi3l5AmwOEAKxpmiR0nVplZmsFIjIUnUsLZa0U35DVGCqL3Dxr70GGO5NPk+uJ+Gpv4yxLGUdkFx++eX47//+b4yPj+O0007DZz/72bQFo9FodCIYSbR06VLcf//9WW8sm6PEULwHibqEQzYA1UUkFleXKiCJRKhWI9VrIQqFLWiVpcvAqNSAn9e0YYxlJ+saEkmSsHv3bmg0Gjz44IPYsGEDPvroo5TXbWlpSTnr5Sc/+QnWrFmT7UOzucrnpqLW8RO0syvADK6Myd1Ixwap10epZpSIIRo68rlo2CTxBBSnm63BTAHZcG/y4wc8wFjqtaIYYyydrDMkJ510Eh544IGJv9966y388z//M5555pmk66lUKtxzzz34/Oc/j9WrV2P16tUTNSSLFy/G0qVLZ7/1rPJFIkAwttKvGCpsoWYm5AyJd4yGJbwl6rkhBCkY0ZuSszTRKHWyLWT9iMxgogDNPZzcol4M0bbJ3WIZYywDWQckTqcTBw8exPLlywEAZ5xxBnp7e1Ned/HixdizZw+ee+45fPDBB6iursYvf/lLnHrqqfjRj340uy1nc0MkDBiMVMyqM5S+UFKrj3dAFQVAilCAUuysjRiihmh1k1ajDvmB0b7itNfX6qhBms4IVNUlbJsADPcAoQAHJIyxjGUdkNx11134zGc+g3POOQcrV67Eu+++i+bm5inXkyQJH374IRYtWoRPfepT+NSnPpV0+c0335z7VrO5IxKmLIlaQ1NsS9WlVabV0yyfiam/Av2/2D030k35FQXqllqs18lWDYxOGp7R6ih7xDNtGGNZyDogWb9+Pfbt24fnn38eBw8exOmnn44rr7wy5XWvu+46HD58GE1NTRPDNvLJ4XDMeuPZHBCJ9SDRxIpaUy17X0waDWVrAh76N+SjoKTYAUnQR/U0k8mZE3WR1s00WgBM6qCsUtHQ0eTF9xhjbBoZ/7ofPXoUixcvBgCYTCZs3rwZmzdvTnt9lUqFF154AQBwyy234K233kJvby/+/Oc/48UXX0RbWxs6Oztnt/Ws8kXCdFKp6VQOQwBGC+AZi2dKRAEodmmL15U6CxIWaTpusRjMsezVpGBRkqgYmTHGMpRxQHLdddehs7MTjY2NWWc7fv/732Pfvn0Tf7/wwgt4+OGHc95oNoeEEzqhlnrKr8xgph2wWk073mJP/Y1EKDOT6rWQosUpaJUZTNTvRAwlByQ6A61CzBhjGco4r7tjxw50dXXhkksuwdDQEHp7e/Gf//mfcDqdWLJkybS3NRqNOHjw4MTfF154IQ4cOJD7VrO5IxKmf6MR2gGXskurTG9A0jBFsQMSMRRbQ2fSMFE0SkFSMQpaZToDBUBTWsjrAL8r/v4xxtgMsh6QzyXbcf/99+Pyyy/Hxo0bsXbtWuzfv78gqwGzChQJUy1CKABYHaWvIQFoJxwWaCcsCoCQpsC0UIQg9T8x25IDASFIWSSdMf1t802lAixVwEh/8tTfaJRWHRaC1N2WMcZmkHXlWy7ZjhUrVmDPnj04++yz8dFHH2HBggV49tlns99aNveEAsBwHw1FlEsNiVZPfTb6jgKj/dSTpJjEEOA6AQx1AwNd8VPvERrO0RcxIAFoJs1ob/K2hEUg4KP3jzHGMpD14ebkbMeBAwcyynbo9XpcdtlluOyyy3LaUDZHhQWgYSF1RjWYyyNDotVT/w29maa9FntlW3k6bWNb8vmuYeqRokmx4F4h6U2Acx69JvJvwfgQ4BnlVX8ZYxnL+tddznY88cQT2L9/P375y1+mbA8vSdLUGzOWLSFIQwFqFfXXKIehPp2edsIaDQ1HRMTYwn9FCgS8Y4DBMjUTogJgrSrONiQymOh1UKni26Q30vaE/MXfHsaYIuV0uJmY7Th+/Dg+9alPQa/X45JLLsGnP/1pnH/++YhGo/neVjYXhUPxlX4TaxRKSaujk9zKPhjrRaIpUjGpZzTNDBsJMBZxho1Mb6K6GjGUHCTJQ1uMMZaBWXdP+vWvf42BgQH89re/hc1mww033IDa2lp87nOfw29+8xuMjo7mYzvZXBQJA+FYf4toBDCVSUACUHAkr9USEYs30yYs0k5+ciM2SQKgKs1aPxoNDSFNmWmjp+LbSKT428QYU5y8tHNUq9U4++yz8eMf/xiHDh3C7t27sX79etx3331obm7GOeecg5/85Cdp17xhLKVIONadVUPDAaVuG5/IaKFARKONzQQqUkAihqhYdPKU37AA6HTFL2iVWaqmBmVaXXw2EmOMzaAg/aVPOeUU3Hzzzfjb3/6Gnp4eXHXVVfjrX/+K3/72t4V4OFapwrG28WotDduUQ0GrTG+MZSViipUhEYKA4J8akIgCoDWULiBJNVQkD2txC3nGWAYK9guv0WgQiURQV1eHq6++GldffXWhHopVqokhGw01/CqHKb+yxCETlap4GRK598nk4t5wCDDZSvca6U1TW8irYp1seeovYywDOWVI/vCHP8x4HZ5lw2YtLAJj/UD/0ViGpIwCEq2eptn2HqF+IAFvcR43FKDXZaAL6OuMn44fpoCkVAwmyoZ0v0/b4x6h7Rzq4Zk2jLGM5BSQXH755bjvvvvwox/9CDt27EAkRdGa3JvkjTfewGWXXYY1a9Zg3bp1uP7669Hb24v29vbZbTmrfNEwdR01WCgAKKsMiYEKOQ1mwFpNXUmLwTtOGQijBTDb4yejlbanVHQGwO6kzJE5YTusVbTNjDE2g5wCEkmSsHv3bmg0Gjz44IPYsGEDurq6plzn4Ycfxg033IDrrrsOL774Ih577DEsXrwY5513Ho4ePZqXJ8AqWCQMWBzUfMxgLK8aEq2edrbWKmoIFhaoXXohSRIQ9FJgVlUfPznq6DUqVf0IQENIVfWUpamqpxWHtTqgppmmRXMbAMbYDHL6hT/ppJPwwAMPTPz91ltvYevWrXjmmWeSrtfR0YHXXnsNVVVVAIC6ujp84xvfwKZNm/DFL34x961mc4O80m8klikppwyJVkdBSVikrIAQoKCkkEGBvG7O5KErefpxKab8JjLZqAg5kc5AAYkQLE2PFMaYYuSUIXE6nUnr2Zxxxhkpp/RGo9GJYCTR0qVLcf/99+fy0GwuEROaounKLEOiVlMAEAlTgBAWaXsLSQwC4eDUwCws0I6/lBkSIB4QJdaP6fT0uvDUX8bYDHL6hb/rrrvwmc98Bueccw5WrlyJd999F83NzVOu19LSgp07d2Ljxo1J5//kJz/BzTffnLL2hLEJoUB85kY5rhhrstCaLdpYQCKEgEL2bhNCgChODczEIGUnJjdLKzaDibYhcQq0WkMLIwo804YxNr2sApKjR49i8eLFWL9+Pfbt24fnn38eBw8exOmnn44rr7wy6boqlQr33HMPPv/5z2P16tVYvXo1PB4PnnrqKSxevBhLly7N6xNhFUgM0c43LJRXl1aZwZw8RFGMDImUohZDFIDaEha0yvRG6oWSago0T/1ljM0gq4DkuuuuQ2dnJxobGyeCjLPPPhurV6+G3Z78gyhJElQqFfbs2YPnnnsOH3zwAaqrq/HLX/4Sp556Kn70ox/l9YmwCiNJFIjIGZJyWccmUVJzMlXhhyUCXurtMVk0QrNsSk2ro/dp8vo1vKYNYywDWQUkO3bsAADccssteOutt9Db24s///nPeOmll7Bw4UJ0dnYmXf+6667DkSNHkgIYr9cLl8uFm2++OX/PglWesAj4PZR18LvLq228TKsHAm4atvG5Cj+91TUSfy3Gh+I1I34XDZeUA4sD6D4YD9bGh2ibNRqaaaMuSHNoxlgFyOnX4fe//z0ef/xx3HLLLXj++efxzDPP4GMf+1jSdVQqFV544QV0dXXhkksuwdDQEHp7e/Gf//mfcDqdWLJkSdaP++qrr+KSSy5Bc3MzVCoVnnjiiWmvv3PnTqhUqimngYGBrB+bFVkkTLMzQn4aqiinGTYyvQGISrT6rhgEvGOFe6ywCHhG6HUB6LXxu+mxoSp9QavMZEluhOZ3U1DpGy/8kBZjTNFyKmo1Go04ePAgli9fDgC48MILsW3btrTX//3vf499+/ZN/P3CCy/g4YcfzvpxfT4f1qxZg6985Sv47Gc/m/HtDh06lDSkVF9fn/VjsyKLiIDZBthqgLGB8urSKtPqAWcjNW6LhKl9uzwFN9+EIBWMzlsKjPQBjW0UhAS8senGZZIh0ZuA2lZ6DQIeoHkJDSmNDlAdSblkchhjZSengOT+++/H5Zdfjo0bN2Lt2rXYv3//RGfWVLINYNLZtGkTNm3alPXt6uvrU04/ZmUsEqaTSlV+C+vJtPrYlN/YtNtAbIipEAFJKEDFoqZJxatiiPp7lHqGjUxvpFWHRTF+njzThhfZY4xNI6df+RUrVmDPnj144oknsH//fixYsAD/9m//lvb62QYw+bZ27VqEQiGsXLkS3//+96cMLyUKhUIIheKpZbe7SC3BWbKwSIWt8mJt5Thko9HEshQemnYrhuhUiCnK8s58yqJ6AmBuzP/j5cpgigVnKdb24V4kjLFp5Fxhptfrcdlll+GHP/whvvnNb6Kuri7p8sTF9eQA5uyzz8ZHH32EBQsW4Nlnn819qzPU1NSEe++9F3/84x/xxz/+Ea2trdi4cSPeeeedtLfp6OiAw+GYOLW2thZ8O1kKcq1ENEIBSTkO2QA0qyQcpmJNKVq4Ool0M2wi4eS1Y0pNXmdn8tRfrQ7wVtBMGyEEBHnRQMbyqWB58OiktSvkAOayyy4r1ENOsWzZMixbtmzi77POOgtHjx7FnXfeif/3//5fytts27YNN95448Tfbrebg5JSSGwbX64ZEoB6kcjBkySl7sGRD75pZhqVS0GrzFJFmZvE90xroBlJkjQ1y6MUQT8w2g+MDVJRMQC0nATUtfLsIcbyoAwH5gtr/fr1eO2119JebjAYYDCU4RTTuSYs0I5L7kFSrj/4eiOAWDZQpaapyvkWFmkHODkgkYO1cgtIjJapQYfeQLNvQgHlrWkjhIChY8DgMRqeM1joNY+EgQ/eBEYHgcWry+99YExh5lxAsm/fPjQ1NZV6M9hMQoH41F9bdam3Jj25VboQ66LqHc3/YwhBWuXXZIvXYQjBWF8PTfnMsJHpjfFskby9kkTBmqCwgMQ1DHTuBdwjgNkBOOrjwZZGS/VC770GQAIWr6XAizGWE0UFJF6vN6n5WldXF/bt2wen04n58+dj27Zt6O3txW9+8xsAtOZOW1sbVqxYgWAwiF/96ld4+eWX8cILL5TqKbBMeceAkX7a2Tdn37OmaLR6wOcBQp2UAdCb8t8ALBQAhnspIJF3hqP91IytaVH5zLCRGUw0XCOKwEBX/Hz3MA172GtKt23ZGO4F9r1MDelsTuoD4xmZej2bEzj4OiXKTjqtfIcXGStzigpI3n77bZx33nkTf8u1HldddRW2b9+O/v5+dHd3T1wuCAK+9a1vobe3F2azGatXr8aLL76YdB+sDEkSnRoXxjIDZdg2XqbVA84GCkRCAdpeMZTffhtikHbitS2UcRjtB5xN9BjOpvKrydCbAGs1cOI40Bqv4YJOT43SlMA7Dny4n4bLHLXU9yUV+f1oPQUY7AKMRmDhqvIdYmSsjCkqINm4cWPS7J3Jtm/fnvT3t7/9bXz7298u8FaxvItGKCAxWoFIpHxn2AC0k9WbaOjEbKOCRyGY34Ak4KXHSKxR0BuBgAawlNEMG5lGQwHUiZ7kbbY4lLGmjRACjr1Pje6cTRREzVQfYrbRqb+LhnYaFxZlUxmrJBzGs/ITFoForGATUnmnwLU6CkoiIm2vEMh/vw0lzbCR2WtoOnQinZFqgoQybiEvSUDvYWAsloXKht5Iw2o9H1DtCWMsKxyQsPIjd2lVq2k4opwDEoAyORPTlMX8BiTpZthEIzSrp9wKWmVGC6DVJJ+nN9JrEyrj/h2uYZpN46inrFe2rFWU1Tv2Hg3hMcYyxgEJKz9hkX7UgVhTtDIfWTRZKBABKIDK57BEKEA1JJMDEjFEMzrKNUOiN1FGJJKQJdFoKZAq14AkEgZ6Y0Xzsxlyq6oH3KNA7xHKuDDGMsIBCSs/kTBNoY1GAXUZd2mV6QzxHY9GRwWR+doRCQEK0CbPpBFj6+ekG8opNaOZaiqmdK5VlW+H0+FeGqpx1M183emo1RSUDHTRQoiMsYxwQMLKj5xtkNvGl/uQTWJQoNXTsES+OrYKwdTBTVigeoVync2hM9BCgJMDEp0B8BSgV8tsBf2UHTHZ85OR0xvp1PNB6nV9GGNTlHkunM1JkTAwfBxwnQBq55X/kI28k/W5KXga7aXMRj6aZHnH6bUIxwIcuVZl4ENg4YrZ33+hqFRAVR3150jMOAS8NHslLJZXoDl4DOg+SNua2GtEfr37OlPfbuL96Jr6fCQJGB+iepRFq/K/zYzlSyQ2lBry00FViZpRlvkvPZuTRAHQm2mHbrSXX5+NyXQGmtIaFimYEsL5K2z1jNCMFXkBPSFIRa4mGxWOljObk9b6SVz8T2+imphQoHwCEp8bOH4IqG6culKzOxacpFvAUH4/jJbU9TwaLfDRfqBhQXlO0WZzmyhQq4LeI/RbE/QB85cDJ68vyeZwQDJX+FxU36CEtt1CAKiuB6ACLLZSb83MdAZqBCaGYi3GbfkJSIQQHbnUzIvvKIUg0H+UdvblWtAqM5hoOx11yUHliR764SuXHfRAFwCJut6m4ndTTUgqQpBuX1Wf/v04+Dd6zpYyzmixuWdsCDj8FmUtLVX0O+Mdz212WZ5wQFLpolFKRx8/RB+02nm0Oqm5jHf0YoiKWfPd8bRQNBrKBMi1AmpNfjqSikF6DUyT3quwCNiMNIulnBnM8bV+JhffyqvllppnjPqOVDUU7jGMNuD4YfrelUsQxuYuSaIA+ej/Av0fAktOjWdbS5yNLtOKOJYXkTDQtR848FcaF9ToKDD58H/jY9/lKBSgnTwk2m4lMCf0ItEZ8hOQyAsMTh7aCAtTO7eWI4OJXotUha3uFGvClMJAF+D3Fnb4y2imRRdP9BTuMRjLRDRKwXHnXupjZK8paUZkMg5IKtlIH/De63SEbXFQVqRmHhXaDfeWeutSi0RoJ6zRKqMpmsxgoe6yQP5m2qRrrBUWaQhHUz4/JClptNQ0bnJnVr0R8I2XPij2jAHd79N3o9DMVRSQKKF1PqtMkQjN+jp2kIaYrVWl3qIpOCCpVEIQ6DtK00ITj6TVGpra2He0PPtBRMK0Y1drKLVY7j1IZHpDPN0pZwVmW0fid6d+/hGRajOUwFY9NUOiN1GdUKkbpJ3oodlRxRgWNJrp8zDUPfN1Gcu3SIRmkfV8QDVdZVoQzwFJpRrqBrxjqY/+rFW0sxs8VvTNmlEk1qVVpaaTUjIkOgNtrxSlbQ6LtNPNlSRRFmFyQzSAghQl1NYAsR++SX1UtDoKUkpZR+IdpyxhutkzhSAvOOgZK95jMhaN0lB9XyfVSpXxbwcXtVYiv4fScvKaMKnG66Ox9F1dS3kVuIZF2lmEAoAUUVZAEonQgnLjQ1QzMJu1TMQQvQ7yeygL+SkbU8Y/Kkn0Jtrmkf7kISb3WGkbhg0eA0b7aXhtunoW+bLxodSXy1mw8aH0NT1hMX4/Y4O0IrC1quQFhGwOkCSqGTn8FmXG/e54fVuqz657mALnEuGApBIN9QAnupNnDqQqshwboixKOQUkkTDtqMRQrE5CQQGJVksZnqCPghHXcO7L0Af9tBPTT1oPxu+hmqByXVRvMqOZgkr3cHIQFRapxqnlpOJvUyhAfReQEBDMVISc7nK5Diboizevm+4+1Brg2AH6XNgVMuzGlGuoG3jvNfrchYXkz2iqz27ASwerJcIBSaURBWBsAJi/go7C+jopLZ2qj4LJRsFLXWv5HK1FROrwKafSlZIh0eoAex0wOgA0tsUaZnlj6/HkMDIa9FHhWf385PPHBilAKfcZNjKdAXA2UYCVeORlqaKATRRSD0sVknuEsjVtq6inS7rvB0BHj3430Lwk9eVyH5LGtvTvyeTv4PBxYLiHAxJWWGNDwHt/AxraUhewpvrsjg+VdGYj15BUGs8oEPBklvUw22kHWk5rbchRu7zTLde1WlKx2OPr8OiN9IXPddjG70o9HU8MAY5a5bwuKhX9GE4pbDUCXlfx60gkiTIzWn3pgnCbk+pXvOOleXxW+fwemtrrHSvL2TTpKORXjWVsbACAKrO55QYT7fi8ZVRkFxYBqGIBiUKGJWRmG9WRAJQZkNuKZysapaN4Q4quuuFYQKIkZvvUNLBGCwi+4s+08XvotbVUFfdxExnMlBkq16n3TNkkiRaKHOkFbKWrB8kFBySVJOinlFs2fRW0OiqyS7WibCkIsaZokbByCjdlOkN8iEmlotc0lx1uyE+ZlVTPX6Uu2yl7aelNAFQUaCUKhylAKCb3SGzhwxIPeVmrgRPHaR0dxvLJMwaM9tFnrFyG4jPEAUklcY/EF/qaLCymDjqMVmDkePF3DOmEAtQ2XopObTde7nQGGgqQX2etPrdGWEEfZVcmP/9oBFBplJc5MsQWSpxc9Kkz0Dh3sYLhSITqN8ohoDNZKfAcGyz1lrBKc6JbWXVmCTggqRSSRFGxNqFBVzRCsxvGBmglx+Ge5BkbAH1oPePU86IcCEHldWmV6WLTcRPrSLxjUzMDM5GHeSYf3YgC3afSfmgMpnhNTSKdAQi4i9egzzdOdRvmInRmzYTJRgFSqTvWssrhHacaKYUN1ch4lk2lCPqA8RO04Jr8wz/SC7Qso5PRTN1Z+z+kTn3yTk0I0g+iZ2zqjI5ii4RpQblolIoglRqQ+D3xSvWAl46EJy9rP50TvdRLbPIO3OeiGSlKC0jUauqBMPhR8udOZ6DnFPAApiJkLUYHgJAPiDoBIVbTEhZpW9J11ZW/H9NdnvhvKukeQ6enIVb3COBszO65MJbKUDfNCjM7cvvsCiHKTpcIBySVwueiKVzWavrb76ZpliedHq+yNttpJ/C/O6mltyqWINMbgIGPgIUrKDtRKmERcI/SDl0F5fQgkekN9Br3HaXXX5JoZxvMIiARBcA1RFNig5NmP7lHgJM3KGeGTSJbNbD/leSaGpUKGB+k16rQO2QhRI0APeNUu5Io6KPvznRmuny0f/rL0z2G6wQw2M0BCZs97zgtnBoKAuGEz1o2n13PKFCdZgp8EXBAUilcw9TvwdlEP756I3Dy+uQpX1od0LqMdgKiSJeND1FmQiXR0Xwpp4iFRcow2JxUeKi0DAkA2GvpS12/gP4ePp5dYWvQRzvN1pOnBodabUl/LGbFYAJqmoHaFvpbCNIPYU0LfQbnLS1sAZ5nhIYwF6xIDugGuqimZLo+JEEf9WpIRX4ezqb0mavpHqO6ARjqAhacUl4NCpnyDB2j38zGRfR3Lp9do7Wkv7sckFQCIUQBib2GPliuE8C8JfEf/0R6I9C0mCJpnYH+Dgs0RFDqgESub9EZKEOiyIDECfSrk7/g2fSbCPriXWon0xqUV9AqM5hph6tSJRfr2qrjnW2NKaY558tIP23D5MfQ6qavy5G/HzMNk013H9M9ht5I69u4hjkgYbnzjNE08pp58c9ZLp9dvSGeOS8BBeZ+2RQ+FwUTRgtFvVodHaGnS+07m6igLpAws0alLn1hayQ2EygaoZk2ShuyAeiLnfiFNpjo6DzTwsXxE7SGz2RhkWoOCrnTLiSjheqbJjeKk88rZHM+nzs2Hb5q9vcVik2tHz5OsxlG++l7NJt223ojcPwDLm5luZEkyo6ExdS9ixSEA5JK4B2jjIJaQ8MFzubpMx1GMxWwJjZEM5gB10j2M0LyKbFLq06fvBibUuiM9D7IOyiTlXa2mWRJRIFmSqX6URGCdN9K/cHRaKi+aXKBnRw0z7SWzGy4h6lYejZ9bUJ+ymQEvJSJXLACWLqOhpqczRSY5Np51Wynxf5cJ3LfPjZ3uUcoO2JXWMPEFDggUbpolKZ5GcyxHboE1M6beTy+ppluIx+xGs00A6GUbeTDIm23Eru0ygwmChzEWM8NuWPudCvKynwu2qmlC0iMFmUOY8ls1fEp0Yl0RgoaCiESpp29YRazeOSAsnkxsLydarPmLaGgfvEaYP0mKh6XopQ5yTZbIgewQ92lPSBgyhMWaaHIaFR5s+9S4IBE6fxuav9utFJ2xFGXWWtxs42OWOWeFzoD7UQnz+woJnm9EyV2aZXJdTmJa7cYLPFF8abjc9GRfKqhtnCIduhKZrQgZcdWg4k+x7mu+zMd9wjVZ2TTvThR0E/Zx7ZVdEqVedTpKVg5eQPgqKegZKb3ejKzg5rEZRK4MiYbPEbZuXRF2QrDAYnS+d00zUujpZ1g/fzMp4U6G5OPWFWq0rayDvrpeUgKjvbVahqmSexKarZRsDFd19ZolBrY6dI8b0nKrpdJOTJapgZr8vlBf2EWmxs4Rv/mMp1dXrF55dlAw4KZr2+xU8akdh71AMomKNHqKPDh9W1YprzjQNd+OrAsZbuGPFLUs3j11Vdx2223Yc+ePejv78fjjz+OzZs3T3ubnTt34sYbb8R7772H1tZWfPe738WXvvSlomxvUQz30Y4s4AWq6ihDkilrFaWLgz5aIt07TkMlrctKswaCd5TW94iGaXxeqaxVwKE3k+sixgcB9ylUf5BKwEtHO0II0En0fsiiEdo5KnUYSyZnvboPxoelBrpoZzw+BDQuBGqa8vd4PjfQd4QyL4mvZ6KwSFmJyTUsUpRqqlb9HdC8JPPvg9EMLFpD2cZDuyljEgmnfozEbQDoPfZ7gKa23DM6bG6QJKD7fWDgQ5pNOT409Try52q6zz4Q/w4ClGUv4UxLRWVIfD4f1qxZg7vvvjuj63d1deHiiy/Geeedh3379uGGG27ANddcg+eff77AW1okQojSw/ZaCizmnUTp40yZbNQHIRSgwjpHHRXvTdd1slAkiY6crVUlnws/a0YLZTPM9vjJXktFkelqBHwuyqrI71/ibXUGCmSUOsNGplJRwKHRxdeTMVroOVY3UDAqhKa/j2yMDdIPtb0m+fVMPMkmny+EgIXLqR9Mto3o9EZgyalA6ynJ36WZtqGqnobsRmZossaYe4QCe0t19p9t+TT5O2i20+9WJivFF4iiMiSbNm3Cpk2bMr7+vffei7a2Ntx+++0AgFNOOQWvvfYa7rzzTlx00UWF2szi8bupGHDeSfEMSTbUaupJ0nuEfgyjUZrlEfQVv4YjLNJU36oGqolRckCiNwJWJwV48pG1pYreK/dI6vdp/AQdTavV9L4mjgm7RygYUdpig6lYqynorKqnI7Oqenq97DW0I/aNA/qG2T+OEKIiUa0+/hip+N30Q5z4evs91ICubXXuQ4cmayzLJ1EX5LrW6cf55fdcZ6CDjMaFyh22ZIUlScBQDw1zNi/J7HOVihBM/g7KShiQKCpDkq1du3bhggsuSDrvoosuwq5du9LeJhQKwe12J53KllwIGPBSwaM1h6JHWzWlzyNh2hlKUnadRfMlLMbrWTQK7UEi0xsBnS65jkSnp+Ck5wNadTaRaxgYH0jfGEvOHClsKfGUjBb6wZtcX5HNbKRMjA1QAJhts7FolIYO552UWXH4dOxOypJAFZ91NROznfqapErBMwZQE7TRvooc1qvogGRgYAANDclHWw0NDXC73QgEUlf0d3R0wOFwTJxaW1uLsanZkyRaLEyjowCipjm3vh3WKjoyl6f7qjXTF18WSmRSQKL0DInWMHUn5KgDjr1H07RlYRE4foSCFGOaqanRSHIKVsmMFsq+TS5sBSirMDow+wZhokBZCYM5+66T3lFaKbUuT9/72nnA0tMAzzB9Z2eiUgF6M2V3JgeujEkSNeQLhysjYzpJRQckudi2bRtcLtfEqaenp9SblJo8K0Gtph/4XI/mNFqaQSBP/9WbqLApkx/PfAqLsUXPJOUHJBpt6p2uRgtABRw7CAR89BoPHgPGppm2F43STkqpDdEm0xuodinVFF+TlbIDs51tMzZAgUW2nVkjYfpeNS2i7cwHtRqYfwqt2ePKsNeKtYoyRYXqzcKU68RxWowxXXG8wlV0QNLY2IjBwcGk8wYHB2G322Eypa6RMBgMsNvtSaey5HdTI7NImDIcszmCtjvj7cr1RppGHCzysE0k1tQtEqasj9KnsVmrqEBxMouDpoS+/wbw7ivA8UPTT9sTAvSepMueKJG9BgineG00WlpTabRv6mWZEgWgv4sCuGyLUT2jVN+T75V3jWZg0SoKUDPJ/siB6/AsXgdWecaGKMNqMCu3T9MMKjogaW9vx0svvZR03o4dO9De3l6iLcojue27FAVqm2dXX2C2UfGfKFAaUAgWv44kEqaMQVisjJ2v2Z46y6RSUQCp0VJK3mCZPpgUAjTrqJIKHM02CjxSvT52J+2Icx02HO2n7Ei29VSRMNX8NC0qTHauupFmGGWa9bA5KdNTyr5ArHz43MBH+yljqvQGidNQ1GGo1+tFZ2d8TnVXVxf27dsHp9OJ+fPnY9u2bejt7cVvfvMbAMB1112Hn//85/j2t7+Nr3zlK3j55Zfx+9//Hk8//XSpnkJ+RCKUtvOO09FmLsWsiYwWGms/0U07R89IbAgny1k7sxEWqVgr6M+sCVW5M5gpqBvpT67tCYtAeDz5PCFh+EIu6pSLGkf7gVPaK6OgVWayUkASFuh5Tg62RvtoFkFblkV7oQBw7ABl+Nwj8Sm3qR5DJvchGR+iGWdVeZjhk4paTcvA9x2NrdWUMCQ0+T2XjcZ6DFnKNEvLikOSgOOHaXi3ujH+OZE/u+mk+1zJUn0/3MMlHQ5SVEDy9ttv47zzzpv4+8YbbwQAXHXVVdi+fTv6+/vR3d09cXlbWxuefvppfPOb38RPf/pTtLS04Fe/+pXyp/z6XPFmaC0nzX7Zco2Wahj6j9LfoQCNdzcunPWmZizgpS9IRKiMegmjhQpb3SdSP5+ZFpOTLxdDgLXCqukNZsBsjQe+4ckzUNRA93v0+cumO+2JHqDvQyoe9rvjwyMpHyNBWAQQBVqWFnZBR0cttZj/YHfqmqHJn4loFOj+AGhYmF1/IVZZxk9QdkSjS16hXZbpb8lkqb4fAe/sVq6eJUUFJBs3boQ0TbHl9u3bU95m7969BdyqEvC7qRZBb6SjrnwcPVfVUce/ulZq4hX00g9ituPwufKN03OBpOyCVpneQCl6vzW54Livc2rfi0TjQ/T+Ni+hHwyTtTKGsBKp1bTEQc8h6mqaKnsx1E1t1FuXZXaffg8w+BGwcGV8CEzus9DYlj5D0tdJGYuFq7PrcpwLlYpWB3aPUFAmB1uJ73miaGyxPtcJmq3D5p5IhLqxOuoA56Quxtn8lqSS6vsxPkTD9yVS0TUkFUmSKDsC0I5qtsM1MqM13h9Cb6QsSbHqSCSJInP56LSEX4i8sjlTF7ZmKuSvvIJWmbV6+myEtQr4cF9mtSTRKGX3gr7si7vluqmGhcUJvm3VtAPwjMy8sq9aTZmRQV4FeM4aG6Cp8IUOlssEByRKE/QDXhcVs9prZz9cIzNZ43UPEyv/FikgiYQpWlepKSiqlPS02UbvU66EAE1dVfqMo1Tk9vrpFqAz24GxE0D3oZl3xkPd1Hckl/oP3zhlLWbbBC0bjQspWPWMznxdq5NqSTKdMswqh98DHHufWjFU4m9AChyQKI3fRUGDSjX72TWJdHo6KpXvW4omF1sWUlikgESKDdcouUtrIqOFnku2S9HLwiLtuCqRyUrB1nRZOGsV0HuICnvTcY9Q91uTLftGUaEADX025WnYM1N6IwVBQmDmDq46PX0veg9zlmQuEQXgyB7AO1KxPUdS4YBEadwj1LNDb8zfcI3MXhP/gVRrUhdQFUJYiA9taLSVkyExmOPDX9mSJACqyhyuASjwdDZN/9podZQ16z1CC+VN3iH7XMDhtyjgy2WFUs8IsGBFaYI+ZxPVa7lPzHxdm5NWds1XW31W3qJRoOd94OheoKqxsmbYzYADEiUJi1R0FJUolZuv4RqZ0ULj1tEo7Uhn2zEzU2EREGON0bQV0BRNptNTJiCX1ZOFIL0H2cwyURpn0/SzXwAKugM+mplyZA9Now4FaOd88A1qu5/LUE3IT5+1pkWl+cFXq2l6u0o98+dDZ6Aaq8Fjxe+gzIpvuJc+1yZH5fwWZmhuPVul845Tr45ohGZw5PuH1GihFXflQkK/l1ZNzVcb7XQiIk33DfqoLqaSjghsNTT7Qw4e5eGpdDshIRjryTJKTcIqtCMjgPhrEvSlXmE0LNK056p6ClwGjwH9H8aazkWpjbalKn1QI7/GqV7rkT5q6V7KJlP2GsDZDBx6c/rAVQjS0NLAhxTEZLuqN1MOv4eWlgDogGa6YDXT35LpLk/8F6Df+9nUvc0SByRKMn6Cpno5m/M/XAPQzk8FamdutFA7eSFQ+IBEFGh620gf0HpyYR+r2EwWYHw4tk5PTNBH0+2m09cJzP/7ygrOJjPb6Qiw+/30Q1OTXyspSv1GJCl+m5ley8k1KEKQZnLVz8992/NBpaIC14/2x6dgpqPVx46cD1EgU6zp+Kx4IhHgowPUHl6eyjtd/RSQ2W9JNt8PzyhQnWYacRFwQKIUkTB9cKobqSAu38M1AP1A1s+nlLizidZcCQVyG5/PRligYaKqusI8r1Iy2eioVq4nGeiiHel0vQMCHqovyHZxOKXR6YH6BfS5q06xfkwmr1XQF+tfk4IQpO+Msym5D8lILzDvpMJ/rjNhq6beKftfoYxNKvLzaD2ZpoCOD+V/vR0liEYpYyaGKKOm1VNgFhbp99FoVu4QhyTR5/1ED7DkNMqCp/rsJirE98NoLWkfKIW+e3OQd5yalZms+Z1dM5klNm6pN9IXvhgzbYL+WDGroXJm2MhMVtrpBH2x11RH/6b7kZFrd6xVc6NleHU9MHI89euRyWsVFmZe5yfxPsQQTaOsa5nddueLSkV1LD0fUPZnui7FZnus50oX9aUoZFfZUgn4KCD3e+ILbobDdGAkBGg42eeiadpaPdXgRMJANNY/yeqkz5S9tvCZ3Xw60UOZwqqG5OG76T7/hfh+6A30mpYIByRK4R6hQjxLVWGGa2QGM/1IRqP0gfeMAU0z32xWgl464lGpKmeGjUyloh+ZY+9lfhshlpWqpAX10jHb41Oji3F06xkDqhvKazq1rRpYsJz6qcy0bIKjlo5qxwYqp3urKNDv20gf0NtJM+4sjnhdkUoNaLVU3xaN0ufE6qSdrSTR90Slpu9N/4fAe6/Rkhp18+m1Ndspe1KuRvqArv3J3XvnKA5IlCASplRtJELpuUIeORst9AUXg4AudrReyBby0SgVUqnVgEpTeRkSgH5cocq8j0Q0Omc6M8Jsox/igLfwQyiRMJ3q55dfDUZjGwUaIf/0QYlGS9/Lvk7lZQEmkySazt3XSfVxWh397jjqKGhMJSxQ5kSnn3rwotNTljXko9/K44cBSFQbZ3HQb6cp9nkzmEr/GZAkyox8dIB+98phCLHEOCBRAu84pSkNJppdU0h6I31hhSCltn3jlC41FagfRliMz5LQ6iovQwJQAGk0Z9aKPxKm1zrbFuhKpdECjhpqj17oH2TvOD1WMbuyZspaRXVDfUeBuhmO5u01wHAP7czmpVmnpNz5XFQDMdRNmY/aeZQRCfpmPxytUsfX+pIk+t55xumgDirKpFTV0+ttq6brFrtuIixSNqjvCGCwcDASwwGJEoz0AUEPFf7ZC/xjqlJROrvvKA0PCbGx20IFJBGRvpySRGPClZgh0Rvph29scObrhgL0+s+F+hGZrYbqIgopGqUj5wWnlG/hY10rTWUO+qZviKdW07Bt/9HCZ0zzLeinIGTwIzrocdQVdmhSFWsuKL+ekgR0H6QhL7+bLjfFFsC0OuPXLVTmKRIBxgepZujY+0DbqsptfpiDMv1msgl+DzWDgopW4y1GBkHu8yA3Scul02imRCG+DLZWVxkr/abiqKMj2pkEfXTUW647zUIw22Lp+lD27d8zJa+QnW5GQjmwVlGxbd/RmXdSZjtwopt6kyxaU/7Tw4UgBVtH91EW0NlUmmFJlYo+Y/IquZEwfecGu4G+D+N1K1UNNHnAUZufXkCiQLNeho7RukR+Nx2kcDCSZA796pURIURFjguWz3x0MNoPHDtAmZFifYGN5lgRbYAe31/AFvIRERjrByTEFtcrs7H9fDHbAZ+bjtDCIv0gTRYJA+FQ5RQrZkpeh6br3eSC7bBIn8NUr5V8OUA1CNNdPtBFwzWnX1j+hcL184GP3qPXQq4lSXweiQG7KACde6kvUQl7R8xobJC67B4/RPUvRitlKMZSXDef73m6g5vpHiMapVk+Ax9StrammaZl187LvTDW56ZOw72d1OfJbMdETVmhnkfifWTzGJ7Rkg4fcUBSCmKQmhzpDMDCFemvJwSp+lprAJoWF++DYjDTkUEkTEeVriEAywvzWGGRvvhCALA6CvMY5cBsp5T8sYP0PqaqEXGPAI2LCj8sV27Uavp8D3Unvy7y2i3p6mlmulwIxmoS1LRDcRZ6ulgeWBy0vs7hN+PPS34ecsF5orHBWLM0Z/ll1cIiZXs+/F8KnuoWUD1Hquchy9d7PpvHiEboN9DmpODkwF/pO9nYRlkTuTVCJnxu4NBblB1KLKYuxvPI5THCQuquyUVSZp/gOUJnoDRg90H6kKZrBjY2SIFL/XygeXHxsgd6I30ZhSB9Mb1jhZuWGRYoEBn20pFypdJoaCim5wPKAkweOojGVleef3Jl9peYia2aZlZYHPFhG787nlpPx+9Of7nc/VSrpe+PUpruNSygujF5LSP5eVTVT9052Zx03ZF+oL61NNubimeMPutyM8eZnocsX+95Ph6juoFOkkSzwAY+TG5G5qihwn+DmYbSVSqqEQl6Y31TgjRMG/RSrcjkYbViPY9sH4MDkjlGb6Ro+60PKI23ZO3UD6vPTdmRoI+6MhZz7Fuloh1Ef+zLN36CvmCF+EEXglRlr1IrewpjJqobqTOp3z2106bfDZjLvMahkMx2CkYCnvzWkYghwGyl+iulsNgpuOg5NHNfCo2WhhI+fJc6HZd6SEoI0U647ygNP9bMK7/MTbZUKvrtM9vowMw1DLyzg+p9DGb6vOoM9D74PZQFOtFDRfrVDTSkVu41PmWiQgfsFaC2BahfSOnWkb7ky+Qak5E++pFpWFj8o2ZT7Mun1ccj/kII+Skg0erosSqZTk9DdEIgedVWSQL8LtoJlXqHUipqNQ2pBDOYGp0NvxuonV/aRfRyUdca38HNxOqko/fjR0q3GnBYBIZ6gL0vUv2LVlcZwchkGi0NuVbV0/Nz1NHvlhgCxoYo02l1UhbU2UQHIZVaF1cA/EqVit4ALF5NO+Sj+2JHFCLNzz9+iFKdai3Qsqw0R80Gc7yFsBgqXAv5oB9Qq+JHGZVOTgN7RunvaJSG5ew1VEA3l9mcsdk2aVbvzZYYAkx2oGG+8o5QzTbKpnlTVX5OolZThqlr/8yLseVbNEoHTh/sBva9RLUSta2xZoBzgEYbWx7CSd9ra1Vl9lIqkgoLXxXG2URFjCO9dFQx1E0BihiKde5z0IJbpYiwDWYKmsRQfAw13yJhmmUjxR6v0jMkANUOtZ4C9B6maZtQUTCyaM2cbxsNiyM2bOMGdHko7PWMAvOX0eurRHWt9BlJN4sikd5I36fuDyg4KcZnKSxSN9TeTtoJ22qp1wtnBFiO+JNTSlodsPRUGosMi1QQZXVSLYFaDZxyZum6ShpMlLEQQxQoZHKkli1RoIWz5PUo5sqRRfMiYPlZlP2qbgAWr1VWc6tCUatpKDOTjrYzEQXKisxbqrzsiMxkoeFanyuzoRhrNQ39HTuY/6GvyQJemj1z/BD9RlXVz81ibJZXnCEpNYsDWHwqcPhtapqjVlO2YNn60vajUKvpKEs++nGP0I98PoOGsEBHsUEfBWWV2hQtFYsdsCynHY1Sd5iF4KilYYDBY/E+C+nIl40PTb1sbABoXa78NYFq5lFfEu84Pc90NUZhkb5Llirgo/1U97V0Xf6HTgI+YPg4Fdy6TwCOesrg+N3Tvx9AfAXbmZ5Hru85P8bsH8M9XNKMIgck5cBiB07eALQuozFJrb481jaw19APm7WafuyEQH4DElGgIzohVNgVjMsZByPJzDagpoX6cMifiZmGLCZfLgRp6mLLEuUPH5gsVAg9+BFlJcIz1NeEfHRA89F7FDycvCE/TdOEIBWtfvi/VKcir0ybKps1U6OuoG/m55Hte86PkZ/HCHipD0uJcEBSLkyWwq0XkyuDmbI0okDZi1Agv0dcYoiq0EWhPAIwVh4a5lP33lAgPqMhlfEh+rFtnrTA3IkeoGkRFRpWgpp5sU6htvTBRV9nck+K5qX0+nS+Qx2h61pzC34jYSpa7e2kvkmOOuCU9tT3le79kMl9Lxrb0h/1T34e/BjFfYzxoZLW8in88IEVlNzwR46k5RRfvoQFACpaNycf60WwymCvoR/aYA6F1H43NatqWFA52SeTBWhbTdnETKf1qlRUn6TRUXv5A3+jovmAL7PbR6M0jfXw28CRdygwsdfS+1IprysrO5whYekZTPTjHh4E1PbMqv2zEfDSEJVaPTdm2LDM6PSUEfjoQHZ9NUSBCkDbVlXetNO6VmDgI3p+2WQTrVV09Nv/ITVPa1lChcO181IHF3Ljr6FjwOhgvD+MRluYwnbGEnBAwtLT6mjGjyjQjBvveH6LMIOxKYJqzdzoQcIyVzuPshyeUTrSn0k0SsM8TYtoZkqlMVmAlqVUw2G2Z1cbozdSsbBOD2j0NDNm8Bg1i3M20eWSFE/h+9w0Y8ZRN3dmvrGywAEJm57dCURiAYkYpLqPfHQTjUTizdZ0Bs6QsGQmK800e/cVyqRN11dDkigYqaqnqdSVOv20rpXqOdzDuTdLlFughwK03sxoPwAVDZ+O9AKNi2lZg0rrsMoUgT91bHpGC2VKdHoasgkF8hOQiCHqQQLQ/Vf6OjYsezVN1KOl+2DsM5jiMyJJwGgfZQ0WrqzsWiSdnhYJ/OCt2U/BN5iSXyshSPdpr+FghJUMf/LY9PSxOhIh1j5eCADIwxTdsAAIfkBS0Tg3/wiyVOYtoZ3l8UMUdCQupR7w0FH94lMpcJkLM7WqG2mJgcEumn0jC4v0OqUrPBeC8eukuzzx31TK4TECXqqjcQ1Twa5GQ0PIKjX9q5TnUa6PIYRokkGJqCSpVKsx5e7uu+/GbbfdhoGBAaxZswY/+9nPsH79+pTX3b59O7785S8nnWcwGBAMZjZjxO12w+FwwOVywW6fg900hSDw+pM0rux3A+2XAC0nzf5+x4aAvz5GWZJTzgRO2TD7+2SVKRql2oYDfwM8I5iYmWW2USfWJadWXhHrdNwjwBtP0c6mkhdjlKTYOlpBOoDRaGjozhx7r4UAnR8F/SsG6ODJaImvw8Wy4xkFTjodWHFW3u4ym32o4g5Lf/e73+HGG2/Evffeiw0bNuCuu+7CRRddhEOHDqG+PvW4qt1ux6FDhyb+VvG0tczpDPRjr4oVuXnyVGkvhmgaIVSV0y+CFYZaTUMV1ioKjINeClJqmmmIQenNz7Jlr6EgrPsgLWQHUMBmtEzfkyLoo54TqQhBqieRi1xTKcZj9H9I76fOCEhhKry1OOi3x2SlxzfG+jWFBZoVJEkUnLlO0NTm3iNAVUP6WphKea0K8RhGa0k7ZisuILnjjjtw7bXXTmQ97r33Xjz99NN44IEH8K//+q8pb6NSqdDY2FjMzawcKhV9YANHqKNs0EcFqbMtHAwL8ZoArh9hmbDXKHehvHxrbKOdSVigYSytjnYq6XZeeiNdd6aMynT3UcjHkFvfB7xAbTNQN5+awFkcsYxHioPIyTU01fXUBKxtFXW2dceyafaa5MdT+mtVyMfQG0qaXVJUQCIIAvbs2YNt27ZNnKdWq3HBBRdg165daW/n9XqxYMECRKNRnHbaabjllluwYsWKlNcNhUIIhUITf7vdee69oUS2ajoS0Rlo2EYIzH410YCXxoAjIs+wYSxbZhsFJR8doKNapZIkYGyQ2pVX1QGL11DGNNfiZL2B1sVyNtFspBPHaV0jTzS2AKCidnlzjqJyncPDw4hEImhoSO5L0NDQgIGBgZS3WbZsGR544AE8+eSTeOihhxCNRnHWWWfh+PHjKa/f0dEBh8MxcWptbc3781AcvYl6hWj1VIkfCsz+PuUeJCo19yBhLBf18+noP93iaeUuEqaF+sxWqiM7eQP1n8nHTCmNhvrXLD2N7tdeSwXQPtfs75sVjKICkly0t7djy5YtWLt2Lc4991z86U9/Ql1dHe67776U19+2bRtcLtfEqaenp8hbXIYMZvqRCAtUTCjMMiCJRGj8UpIo7coZEsaypzcC85cDahUF+EoihihAqGmiVYmr6wtTC6RSUVO4ZadT+/1oGBgfzM9BFcs7ReWvamtrodFoMDg4mHT+4OBgxjUiOp0Op556Kjo7O1NebjAYYDDwEXsSuYW8vIrqbFvIh4XY+jgSBSOVPFOAsUKqqqNmcAPH4mtOlbuAh35LmhcDLScXp4ZMq6Mp5NUNVPDa8wHV4FTV028aKwuKypDo9XqsW7cOL7300sR50WgUL730Etrb2zO6j0gkgv3796OpqalQm1l51GpaCl4MUvDgHZ/d/Ykhqh2RJCpYq9TOmowVQ8NCmnXjPpHd2j/FJkk0S08MUhO7BSuLX9ButlGdytqPU23c8HEg5C/uNrC0FJUhAYAbb7wRV111FU4//XSsX78ed911F3w+38Ssmy1btmDevHno6OgAAPzgBz/AmWeeiSVLlmB8fBy33XYbjh07hmuuuaaUT0N5bFXAwdcAkw0w2amBTq4/JmKsTbUQAhoX5XUzGZtzNBqaWTI2CBx+K/U0ejl70pc6Mzxx+UBX+mmfYZFmrqTLkE73GJEwbV/TYqrpcJZ41mNVHc3g6e0EPthFw8hyL5tSv1alfAzPaEkbDCouILn88stx4sQJ/Md//AcGBgawdu1aPPfccxOFrt3d3VAnjEWOjY3h2muvxcDAAKqrq7Fu3Tq8/vrrWL58eamegjIZzNSQyOak8WohMIuAJEhDNRodFbQxxmbHZKHmgvv+QjNJDObky90j9K85TWMqIUjfa6Ml/RDqTPeR7nK/m3Z+p5xJgZPZNvPzKQadHlhwCrUzOLKHsjeOuviqxqV4rQr9GGExviZZwEf/Gs1UWwNVydcVU2Sn1mKa851aZT43cOCvgNVJRWEnr6fGVLno/gDoOQRAApa3Z7aaK2NsepJE36vug0B1U3KfjvEhCgyal6S+rbzSb2Nb+h1gXyft3KZrxJX4GGKIzjOa6bz6+eU77dY9Chx7j5qraXQ0jFPM1yqfjxGNUNHuSB8FV7XzYkN5Ej03vYmCQmtVbAalOqHPi4rerzx2Pq7oTq2sRIxm+mKIQQAq6iOSK5+LvgRqNRe0MpYvKhUVboYCtNZNbUtpCjajEVprJhoG6hdQ8aqlzA/m7E6a7dP9PtC1v7xrcVKJRul3NeCh99xgokClcSHV/+mNlP0wmOILppYhDkhYZjRaisbHT9CH2z2a2/1EIpSKVKkoNajjgISxvNFoaRgiHKIj5Jp5xWutL0l01D98nIY+mmOzWpTS2t9oBhatpqaPR/bQ61fOzdQkCQj6gdEBmrlocdD2W6umH+opY2X6SrOyZK0Ghnvp36A3tyXQhQClcgGK2GezhDpjbCq9kWaxRMJUPO4swoxCIQiM9VOjtsVrKTujxO+2Vge0LqPeJT2HY69fc3llFEJ+yoZ4xqjNvrMJcDZQEKjE1zwBByQsc3K7eL0RGHdTpiPrgCS2cqdaQxkXXuiQsfwzWYHFpwLHDtBBhLpAP/WRMBVPRiPA0tMpK2KyFOaxisleQ11ejx2k4a+qhtJnHCKxpm4aHQUhS06lg0OjeebbKgQHJCxzBjMdKUSjND4c8tNc/mzIHVqjkfKptmesEpkswKK1gEYPdL1LBwL5EhYBzwj966iND89U0gGG3gC0raTfvP6jtGZQsafESlJ8xk3QS4FIyzKqealAHJCwzBnNgArAiR76cvg92d9HyE9z3bmglbHC0xuorsDuBD7YDXz4v/Eix0RCkP4dH0r/vQyL1Etk/ASgUQP1C4HmReVdZzFbWh2wYDkdPHXuBU50Uy8mYObXSp52m4p8Wap1iOTbekfpfbFV06ltNdCwoLyGj/KsQj9FrCA0WvoxGzhGAcVoHzD/5Ozuw3WC6kh00yybzRjLH42GdmTWaqD3MPDRe1SwKRc+qlTxJllBX3ImRZLobyFIM+uMZjpCn39KrO26QgpWZ0OtptfPbAeO7gM+fDfe52WmrNNMy2zIl0ejdLAW9NNrXDuP+rZYqygYMloqK/uUBgckLDvVjdR/xFod+6HKomNrNEqFsA0LqYaEZ9gwVjwWO01tbVgIDB6jqbmCn1bc1ugoSKltoSPwSBhAbOqrxUE7SUc9Dc/YnHNzuQdbNbDiY/T79/4bVD9XMy/1gVWmfUiqGwHfOJ1XVU+LDDqb6Pd1Dr7GHJCw7MhHBnojfaFC/swDEiFIAQlUvKgeY6Ugr37rqKVOnd4xyooEfLTDNVlp9pveRN9rnZGCEYN5bmRDZqLTx2fhDByjmUXuMHWxNlkzy2IIQRq2ltu0Ny2iwMTiUPwsmdnigIRlJ6mwNdZTJNPCViHWslijpR+5OXgEwFjZMFkqY0ZMKdhrKFPkWUhD1yN9VFun0dKBlhCkIthEfg8FgHoDUNMUX9PHxMtnyDggYdkxmunoSQjQEdNMY6SJhAAgRWm82lZTuG1kjLFCU6moWNjupEVC3SOA30X/6oxUlBqJ1eZIUaoDmX8KBSMmG2ecUuCAhGVHo6Ujg+FeCkxcJyhbksmXSwjSl1iK8pEZY6xyGM2xfiCt8Q6qIX9s9lKsFsdaza0OZsABCcuerZqaBVmr6Igg4M1srQqvi76sWj0dLTDGWKVRqXg4LEecM2LZM9koKxIJUzASyKAfSSQMeIapqFWr44CEMcZYEg5IWPbkIqyeD4ATx2n64EwCXir6Gj5OxV6V2kiJMcZYTnivwLKn01OjIKiob8HgMVrMa7pZM34PBSJGK1WWM8YYYwk4Q8JyY6+lQlabkxr7BL3TX989TEWwOgMP1zDGGJuCAxKWG7ONirc0GsDvBXzTTP+NRGhWjkbDAQljjLGUOCBhuTFaKbAIBQApMn0dSdAHeMYAqKixmqFylstmjDGWHxyQsNzoDbRWQ9BHAcZYX6wtfAoBDw3pSFHAVsUNgRhjjE3BewaWO0cdtYI3mClDkm65bc8YdWeNiIClqqibyBhjTBk4IGG5s1XTjJtoGIhKtJ6DJCVfJxoFRvtpdV+tnhaQYowxxibhgITlzuIADCbAM06Nz45/MHVtG7+bAhUxRAWt5gw6ujLGGJtzOCBhuVOpaFEpKUIrXIaCwNhQ8nUGj9FQjkoFNC3m+hHGGGMp8d6BzY7NSbUkRgtQvwAY6aV6EYBqR0Z6gZpm6ltid5Z2WxljjJUtDkjY7JhtQM08mkVjsQOeUeD4YVrlcuhYrJg1DNTN4+EaxhhjaXHreDY7KhXQ1Ear/6o1lC05fog6swa8gK2GhmwaF9F1GWOMsRQ4Q8Jmz15DgYdnjKYA17YCQT+gUgPRCGB20IwcxhhjLA0OSNjsmazAkrU0bCOGqHC1qh6wVgOeEWDBch6uYYwxNi0OSFh+1LZQUev4IPUeiUaBsX46r2UpD9cwxhibliIDkrvvvhsLFy6E0WjEhg0b8Oabb057/cceewwnn3wyjEYjVq1ahWeeeaZIWzqHqNVAy0k026b3MNDzPq13M/8UQMOlSowxxqanuIDkd7/7HW688UZ873vfwzvvvIM1a9bgoosuwtDQUMrrv/7667jiiitw9dVXY+/evdi8eTM2b96MAwcOFHnL5wCTBVjxd8AZnwTWfwpYfhZg5IX0GGOMzUwlSZN7fZe3DRs24IwzzsDPf/5zAEA0GkVrayu+/vWv41//9V+nXP/yyy+Hz+fDU089NXHemWeeibVr1+Lee++d8fHcbjccDgdcLhfsdq6DYIwxxjKVzT5UURkSQRCwZ88eXHDBBRPnqdVqXHDBBdi1a1fK2+zatSvp+gBw0UUXpb1+KBSC2+1OOjHGGGOssBQVkAwPDyMSiaChoSHp/IaGBgwMDKS8zcDAQFbX7+jogMPhmDi1trbmZ+MZY4wxlpaiApJi2LZtG1wu18Spp6en1JvEGGOMVTxFTX+ora2FRqPB4OBg0vmDg4NobGxMeZvGxsasrm8wGGAwGPKzwYwxxhjLiKIyJHq9HuvWrcNLL700cV40GsVLL72E9vb2lLdpb29Puj4A7NixI+31GWOMMVZ8isqQAMCNN96Iq666CqeffjrWr1+Pu+66Cz6fD1/+8pcBAFu2bMG8efPQ0dEBALj++utx7rnn4vbbb8fFF1+MRx99FG+//TZ+8YtflPJpMMYYYyyB4gKSyy+/HCdOnMB//Md/YGBgAGvXrsVzzz03Ubja3d0NtTqe+DnrrLPwyCOP4Lvf/S6+853vYOnSpXjiiSewcuXKUj0FxhhjjE2iuD4kxcZ9SBhjjLHcVGwfEsYYY4xVJg5IGGOMMVZyHJAwxhhjrOQUV9RabHKJDbeQZ4wxxrIj7zszKVflgGQGHo8HALiFPGOMMZYjj8cDh8Mx7XV4ls0MotEo+vr6YLPZoFKpplzudrvR2tqKnp6eOTELZ649X2DuPee59nwBfs5z4TnPtecLlMdzliQJHo8Hzc3NSS05UuEMyQzUajVaWlpmvJ7dbp8zH3Jg7j1fYO4957n2fAF+znPBXHu+QOmf80yZERkXtTLGGGOs5DggYYwxxljJcUAySwaDAd/73vfmzArBc+35AnPvOc+15wvwc54L5trzBZT3nLmolTHGGGMlxxkSxhhjjJUcBySMMcYYKzkOSBhjjDFWchyQMMYYY6zkOCCZhbvvvhsLFy6E0WjEhg0b8Oabb5Z6k/Kio6MDZ5xxBmw2G+rr67F582YcOnQo6TobN26ESqVKOl133XUl2uLZ+/73vz/l+Zx88skTlweDQWzduhU1NTWwWq343Oc+h8HBwRJu8ewtXLhwynNWqVTYunUrAOW/x6+++iouueQSNDc3Q6VS4Yknnki6XJIk/Md//AeamppgMplwwQUX4MiRI0nXGR0dxZVXXgm73Y6qqipcffXV8Hq9RXwW2ZnuOYuiiJtvvhmrVq2CxWJBc3MztmzZgr6+vqT7SPW5uPXWW4v8TDI30/v8pS99acrz+eQnP5l0HSW9zzM931TfaZVKhdtuu23iOuX6HnNAkqPf/e53uPHGG/G9730P77zzDtasWYOLLroIQ0NDpd60WXvllVewdetWvPHGG9ixYwdEUcSFF14In8+XdL1rr70W/f39E6cf//jHJdri/FixYkXS83nttdcmLvvmN7+J//mf/8Fjjz2GV155BX19ffjsZz9bwq2dvbfeeivp+e7YsQMA8IUvfGHiOkp+j30+H9asWYO777475eU//vGP8X//7//Fvffei927d8NiseCiiy5CMBicuM6VV16J9957Dzt27MBTTz2FV199Ff/4j/9YrKeQtemes9/vxzvvvIN///d/xzvvvIM//elPOHToED796U9Pue4PfvCDpPf961//ejE2Pyczvc8A8MlPfjLp+fz2t79NulxJ7/NMzzfxefb39+OBBx6ASqXC5z73uaTrleV7LLGcrF+/Xtq6devE35FIRGpubpY6OjpKuFWFMTQ0JAGQXnnllYnzzj33XOn6668v3Ubl2fe+9z1pzZo1KS8bHx+XdDqd9Nhjj02c9/7770sApF27dhVpCwvv+uuvlxYvXixFo1FJkirrPQYgPf744xN/R6NRqbGxUbrtttsmzhsfH5cMBoP029/+VpIkSTp48KAEQHrrrbcmrvPss89KKpVK6u3tLdq252ryc07lzTfflABIx44dmzhvwYIF0p133lnYjSuQVM/5qquuki699NK0t1Hy+5zJe3zppZdKH//4x5POK9f3mDMkORAEAXv27MEFF1wwcZ5arcYFF1yAXbt2lXDLCsPlcgEAnE5n0vkPP/wwamtrsXLlSmzbtg1+v78Um5c3R44cQXNzMxYtWoQrr7wS3d3dAIA9e/ZAFMWk9/vkk0/G/PnzK+b9FgQBDz30EL7yla8kLSJZae+xrKurCwMDA0nvqcPhwIYNGybe0127dqGqqgqnn376xHUuuOACqNVq7N69u+jbXAgulwsqlQpVVVVJ5996662oqanBqaeeittuuw3hcLg0G5gnO3fuRH19PZYtW4avfvWrGBkZmbiskt/nwcFBPP3007j66qunXFaO7zEvrpeD4eFhRCIRNDQ0JJ3f0NCADz74oERbVRjRaBQ33HADPvaxj2HlypUT5/9//9//hwULFqC5uRnvvvsubr75Zhw6dAh/+tOfSri1uduwYQO2b9+OZcuWob+/H//n//wfnH322Thw4AAGBgag1+un/Gg3NDRgYGCgNBucZ0888QTGx8fxpS99aeK8SnuPE8nvW6rvsHzZwMAA6uvrky7XarVwOp0V8b4Hg0HcfPPNuOKKK5IWXvvGN76B0047DU6nE6+//jq2bduG/v5+3HHHHSXc2tx98pOfxGc/+1m0tbXh6NGj+M53voNNmzZh165d0Gg0Ff0+P/jgg7DZbFOGl8v1PeaAhE1r69atOHDgQFI9BYCk8dVVq1ahqakJ559/Po4ePYrFixcXezNnbdOmTRP/X716NTZs2IAFCxbg97//PUwmUwm3rDjuv/9+bNq0Cc3NzRPnVdp7zOJEUcRll10GSZJwzz33JF124403Tvx/9erV0Ov1+Kd/+id0dHQopgV5on/4h3+Y+P+qVauwevVqLF68GDt37sT5559fwi0rvAceeABXXnkljEZj0vnl+h7zkE0OamtrodFopsyyGBwcRGNjY4m2Kv++9rWv4amnnsJf/vIXtLS0THvdDRs2AAA6OzuLsWkFV1VVhZNOOgmdnZ1obGyEIAgYHx9Puk6lvN/Hjh3Diy++iGuuuWba61XSeyy/b9N9hxsbG6cUqYfDYYyOjir6fZeDkWPHjmHHjh0zLku/YcMGhMNhfPTRR8XZwAJbtGgRamtrJz7Hlfo+//Wvf8WhQ4dm/F4D5fMec0CSA71ej3Xr1uGll16aOC8ajeKll15Ce3t7CbcsPyRJwte+9jU8/vjjePnll9HW1jbjbfbt2wcAaGpqKvDWFYfX68XRo0fR1NSEdevWQafTJb3fhw4dQnd3d0W837/+9a9RX1+Piy++eNrrVdJ73NbWhsbGxqT31O12Y/fu3RPvaXt7O8bHx7Fnz56J67z88suIRqMTwZnSyMHIkSNH8OKLL6KmpmbG2+zbtw9qtXrKsIZSHT9+HCMjIxOf40p8nwHKeq5btw5r1qyZ8bpl8x6XuqpWqR599FHJYDBI27dvlw4ePCj94z/+o1RVVSUNDAyUetNm7atf/arkcDiknTt3Sv39/RMnv98vSZIkdXZ2Sj/4wQ+kt99+W+rq6pKefPJJadGiRdI555xT4i3P3be+9S1p586dUldXl/S3v/1NuuCCC6Ta2lppaGhIkiRJuu6666T58+dLL7/8svT2229L7e3tUnt7e4m3evYikYg0f/586eabb046vxLeY4/HI+3du1fau3evBEC64447pL17907MKLn11lulqqoq6cknn5Teffdd6dJLL5Xa2tqkQCAwcR+f/OQnpVNPPVXavXu39Nprr0lLly6VrrjiilI9pRlN95wFQZA+/elPSy0tLdK+ffuSvtuhUEiSJEl6/fXXpTvvvFPat2+fdPToUemhhx6S6urqpC1btpT4maU33XP2eDzSv/zLv0i7du2Surq6pBdffFE67bTTpKVLl0rBYHDiPpT0Ps/0uZYkSXK5XJLZbJbuueeeKbcv5/eYA5JZ+NnPfibNnz9f0uv10vr166U33nij1JuUFwBSnn79619LkiRJ3d3d0jnnnCM5nU7JYDBIS5YskW666SbJ5XKVdsNn4fLLL5eampokvV4vzZs3T7r88sulzs7OicsDgYD0z//8z1J1dbVkNpulz3zmM1J/f38Jtzg/nn/+eQmAdOjQoaTzK+E9/stf/pLyc3zVVVdJkkRTf//93/9damhokAwGg3T++edPeR1GRkakK664QrJarZLdbpe+/OUvSx6PpwTPJjPTPeeurq603+2//OUvkiRJ0p49e6QNGzZIDodDMhqN0imnnCLdcsstSTvvcjPdc/b7/dKFF14o1dXVSTqdTlqwYIF07bXXTjlwVNL7PNPnWpIk6b777pNMJpM0Pj4+5fbl/B6rJEmSCpqCYYwxxhibAdeQMMYYY6zkOCBhjDHGWMlxQMIYY4yxkuOAhDHGGGMlxwEJY4wxxkqOAxLGGGOMlRwHJIwxxhgrOQ5IGGOMMVZyHJAwxhhjrOQ4IGGMFYUkSbjjjjvQ1tYGs9mMzZs3w+Vypbzuxo0boVKpoFKpJhb1S2fjxo244YYb8rqtX/rSlyYe/4knnsjrfTPGUuOAhDFWFDfddBPuuecePPjgg/jrX/+KPXv24Pvf/37a61977bXo7+/HypUri7eRMT/96U/R399f9MdlbC7jgIQxVnC7d+/GHXfcgd/97nc455xzsG7dOlx77bV45pln0t7GbDajsbERWq22iFtKHA4HGhsbi/64jM1lHJAwxgruJz/5Cc4//3ycdtppE+c1NDRgeHg4q/vx+XzYsmULrFYrmpqacPvtt0+5TjQaRUdHB9ra2mAymbBmzRr84Q9/mLjc4/HgyiuvhMViQVNTE+68886CDPswxrLDAQljrKBCoRCefvppfOYzn0k6PxgMwuFwZHVfN910E1555RU8+eSTeOGFF7Bz50688847Sdfp6OjAb37zG9x7771477338M1vfhNf/OIX8corrwAAbrzxRvztb3/Dn//8Z+zYsQN//etfp9wHY6z4ip8LZYzNKe+88w4CgQC+9a1v4dvf/vbE+aIo4rzzzsv4frxeL+6//3489NBDOP/88wEADz74IFpaWiauEwqFcMstt+DFF19Ee3s7AGDRokV47bXXcN999+G0007Dgw8+iEceeWTiPn7961+jubk5H0+VMTYLHJAwxgrq8OHDsFgsU2bLXHzxxfjYxz6W8f0cPXoUgiBgw4YNE+c5nU4sW7Zs4u/Ozk74/X584hOfSLqtIAg49dRT8eGHH0IURaxfv37iMofDkXQfjLHS4ICEMVZQbrcbtbW1WLJkycR5x44dw5EjR/C5z30ur4/l9XoBAE8//TTmzZuXdJnBYMDo6GheH48xlj9cQ8IYK6ja2lq4XC5IkjRx3n/913/hU5/6FJYvX57x/SxevBg6nQ67d++eOG9sbAyHDx+e+Hv58uUwGAzo7u7GkiVLkk6tra1YtGgRdDod3nrrrYnbuFyupPtgjJUGZ0gYYwX18Y9/HMFgELfeeiv+4R/+AQ8//DD+53/+B2+++WZW92O1WnH11VfjpptuQk1NDerr6/Fv//ZvUKvjx1U2mw3/8i//gm9+85uIRqP4u7/7O7hcLvztb3+D3W7HVVddhauuugo33XQTnE4n6uvr8b3vfQ9qtRoqlSrfT50xlgUOSBhjBdXQ0IDt27fjpptuwg9/+EN8/OMfx2uvvYbW1tas7+u2226D1+vFJZdcApvNhm9961tTur3+8Ic/RF1dHTo6OvDhhx+iqqoKp512Gr7zne8AAO644w5cd911+Pu//3vY7XZ8+9vfRk9PD4xGY16eL2MsNyopMY/KGGNlYOPGjVi7di3uuuuugj+Wz+fDvHnzcPvtt+Pqq69OukylUuHxxx/H5s2bC74djM11XEPCGCtL//3f/w2r1Yr9+/fn9X737t2L3/72tzh69CjeeecdXHnllQCASy+9dOI61113HaxWa14flzE2Pc6QMMbKTm9vLwKBAABg/vz50Ov1ebvvvXv34pprrsGhQ4eg1+uxbt063HHHHVi1atXEdYaGhuB2uwEATU1NsFgseXt8xlhqHJAwxhhjrOR4yIYxxhhjJccBCWOMMcZKjgMSxhhjjJUcBySMMcYYKzkOSBhjjDFWchyQMMYYY6zkOCBhjDHGWMlxQMIYY4yxkuOAhDHGGGMlxwEJY4wxxkqOAxLGGGOMldz/D1/FQD255QI1AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(6, 4))\n",
    "ax.fill_between(\n",
    "    angles * 180 / np.pi,\n",
    "    kduq_pred_post[0],  # / solver.rutherford,\n",
    "    kduq_pred_post[1],  # / solver.rutherford,\n",
    "    color=\"#ff4500\",\n",
    "    hatch=\"|-|-\",\n",
    "    alpha=0.2,\n",
    ")\n",
    "plt.xlabel(r\"$\\theta$ [deg]\")\n",
    "plt.ylabel(r\"$\\frac{d \\sigma}{d\\Omega} / \\frac{d \\sigma_{R}}{d\\Omega} $\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "d58d9b30-fad9-48ae-b640-dd4841a164f8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.163423\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_CHECK_OUTPUT\n",
    "print(f\"{kduq_pred_post[1][9]:1.6f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc74a182-825d-463c-ae69-4315af0aeadc",
   "metadata": {},
   "source": [
    "In fact, `KDUQ` is just one definition of a global optical potential. `jitr` also has others built in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "9ec83fe1-0baa-41db-83ca-5eee93eb8825",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on class WLH in module jitr.optical_potentials.wlh:\n",
      "\n",
      "class WLH(jitr.optical_potentials.omp.SingleChannelOpticalModel)\n",
      " |  WLH(projectile: tuple)\n",
      " |\n",
      " |  The Whitehead-Lim-Holt global optical potential for nucleon-nucleus\n",
      " |  scattering.\n",
      " |\n",
      " |  Method resolution order:\n",
      " |      WLH\n",
      " |      jitr.optical_potentials.omp.SingleChannelOpticalModel\n",
      " |      builtins.object\n",
      " |\n",
      " |  Methods defined here:\n",
      " |\n",
      " |  __init__(self, projectile: tuple)\n",
      " |      Initialize self.  See help(type(self)) for accurate signature.\n",
      " |\n",
      " |  evaluate(self, rgrid: float | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]], reaction: jitr.reactions.reaction.Reaction, kinematics: jitr.utils.kinematics.ChannelKinematics, *params: float) -> tuple[complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]], complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]], float | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]]\n",
      " |      Evaluate the central, spin-orbit, and Coulomb terms of the WLH\n",
      " |      potential on the given radial grid for the specified reaction and\n",
      " |      kinematics, using the provided potential parameters.\n",
      " |\n",
      " |\n",
      " |      :param reaction: Reaction The reaction for which to calculate the\n",
      " |                       parameters.\n",
      " |      :param kinematics: ChannelKinematics The kinematics of the reaction\n",
      " |                         channel.s\n",
      " |      :param uv0, uv1, ..., aso1: float Parameters of the WLH potential. See\n",
      " |                                  [Whitehead et al., 2021]\n",
      " |                                  (https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.127.182502)\n",
      " |                                  for details.\n",
      " |      :returns: U_central: np.ndarray Central term of the optical potential\n",
      " |                evaluated on the radial grid; U_spin_orbit: np.ndarray\n",
      " |                Spin-orbit term of the optical potential evaluated on the\n",
      " |                radial grid; U_coulomb: np.ndarray Coulomb term of the optical\n",
      " |                potential evaluated on the radial grid\n",
      " |      :rtype: tuple[U_central, U_spin_orbit, U_coulomb]\n",
      " |\n",
      " |  ----------------------------------------------------------------------\n",
      " |  Methods inherited from jitr.optical_potentials.omp.SingleChannelOpticalModel:\n",
      " |\n",
      " |  __call__(self, rgrid: 'ArrayOrScalar', reaction: 'reaction.Reaction', kinematics: 'kinematics.ChannelKinematics', *params: 'float') -> 'tuple[PotentialArray, PotentialArray, PotentialArray | ArrayOrScalar]'\n",
      " |      Call self as a function.\n",
      " |\n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors inherited from jitr.optical_potentials.omp.SingleChannelOpticalModel:\n",
      " |\n",
      " |  __dict__\n",
      " |      dictionary for instance variables\n",
      " |\n",
      " |  __weakref__\n",
      " |      list of weak references to the object\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(jitr.optical_potentials.wlh.WLH)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "cbc59d17-2ff7-4f09-9d53-895d3493acf6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on class CHUQ in module jitr.optical_potentials.chuq:\n",
      "\n",
      "class CHUQ(jitr.optical_potentials.omp.SingleChannelOpticalModel)\n",
      " |  Chapel-Hill Uncertainty Quantification (CHUQ) optical\n",
      " |  potential model.\n",
      " |\n",
      " |  Note that CH89 is Lane consistent, so the same parameters can be\n",
      " |  used for both neutron and proton projectiles.\n",
      " |\n",
      " |  Method resolution order:\n",
      " |      CHUQ\n",
      " |      jitr.optical_potentials.omp.SingleChannelOpticalModel\n",
      " |      builtins.object\n",
      " |\n",
      " |  Methods defined here:\n",
      " |\n",
      " |  __init__(self)\n",
      " |      Initialize self.  See help(type(self)) for accurate signature.\n",
      " |\n",
      " |  evaluate(self, rgrid: float | numpy.ndarray, reaction: jitr.reactions.reaction.Reaction, kinematics: jitr.utils.kinematics.ChannelKinematics, *params: float) -> tuple[complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]], complex | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.complex128]], float | numpy.ndarray[tuple[typing.Any, ...], numpy.dtype[numpy.float64]]]\n",
      " |      Evaluate the CHUQ central, spin-orbit, and Coulomb terms.\n",
      " |\n",
      " |  ----------------------------------------------------------------------\n",
      " |  Methods inherited from jitr.optical_potentials.omp.SingleChannelOpticalModel:\n",
      " |\n",
      " |  __call__(self, rgrid: 'ArrayOrScalar', reaction: 'reaction.Reaction', kinematics: 'kinematics.ChannelKinematics', *params: 'float') -> 'tuple[PotentialArray, PotentialArray, PotentialArray | ArrayOrScalar]'\n",
      " |      Call self as a function.\n",
      " |\n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors inherited from jitr.optical_potentials.omp.SingleChannelOpticalModel:\n",
      " |\n",
      " |  __dict__\n",
      " |      dictionary for instance variables\n",
      " |\n",
      " |  __weakref__\n",
      " |      list of weak references to the object\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(jitr.optical_potentials.chuq.CHUQ)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a477e6f-b8bd-41bb-9b63-f1e0c07f6536",
   "metadata": {},
   "source": [
    "In fact, these are all derived from `jitr.optical_potentials.SingleChannelOpticalModel`, which sets the interface between the model for the effective interaction and the solver. You can come up with your own optical model and derive your own class from `jitr.optical_potentials.SingleChannelOpticalModel` to plug it into jitr!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "7c602e29-abdd-4662-8419-ca3a99bc701a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on class SingleChannelOpticalModel in module jitr.optical_potentials.omp:\n",
      "\n",
      "class SingleChannelOpticalModel(builtins.object)\n",
      " |  SingleChannelOpticalModel(params: 'list[str]') -> 'None'\n",
      " |\n",
      " |  Base class for local single-channel optical potentials.\n",
      " |\n",
      " |  Methods defined here:\n",
      " |\n",
      " |  __call__(self, rgrid: 'ArrayOrScalar', reaction: 'reaction.Reaction', kinematics: 'kinematics.ChannelKinematics', *params: 'float') -> 'tuple[PotentialArray, PotentialArray, PotentialArray | ArrayOrScalar]'\n",
      " |      Call self as a function.\n",
      " |\n",
      " |  __init__(self, params: 'list[str]') -> 'None'\n",
      " |      Initialize self.  See help(type(self)) for accurate signature.\n",
      " |\n",
      " |  evaluate(self, rgrid: 'ArrayOrScalar', reaction: 'reaction.Reaction', kinematics: 'kinematics.ChannelKinematics', *params: 'float') -> 'tuple[PotentialArray, PotentialArray, PotentialArray | ArrayOrScalar]'\n",
      " |      Evaluate central, spin-orbit, and Coulomb terms on ``rgrid``.\n",
      " |\n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors defined here:\n",
      " |\n",
      " |  __dict__\n",
      " |      dictionary for instance variables\n",
      " |\n",
      " |  __weakref__\n",
      " |      list of weak references to the object\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(jitr.optical_potentials.SingleChannelOpticalModel)"
   ]
  }
 ],
 "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}