{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import torch\n",
    "from torchvision.transforms.functional import rotate\n",
    "from pytomography.utils import rotate_detector_z\n",
    "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Last Time"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finished implementing operation $g=Hf$ and $\\hat{f} = H^T g$ for the simple case:"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* $g = Hf$ (where $f$ is `obj` and $g$ is `image`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = torch.linspace(-1,1,128)\n",
    "xv, yv, zv = torch.meshgrid(x,x,x, indexing='ij')\n",
    "obj = (xv**2 + 0.9*zv**2 < 0.5) * (torch.abs(yv)<0.8)\n",
    "obj = obj.to(torch.float).unsqueeze(dim=0)\n",
    "\n",
    "angles = np.arange(0,360.,3)\n",
    "image = torch.zeros((1,len(angles),128,128))\n",
    "for i,angle in enumerate(angles):\n",
    "    object_i = rotate_detector_z(obj,angle)\n",
    "    image[:,i] = object_i.sum(axis=1)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* $\\hat{f} = H^T g$ (where $\\hat{f}$ is `obj_bp`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "obj_bp = torch.zeros([1, 128, 128, 128])\n",
    "for i,angle in enumerate(angles):\n",
    "    bp_single_angle = torch.ones([1, 128, 128, 128]) * image[:,i].unsqueeze(dim=1)\n",
    "    bp_single_angle = rotate_detector_z(bp_single_angle, angle, negative=True)\n",
    "    obj_bp += bp_single_angle"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This Time"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start making $H$ more representative of true image modeling. In this tutorial, we'll consider attenuation modeling; in the next tutorial, we'll consider PSF modeling"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"../images/emission.png\" alt= “” width=\"400\">"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose the green is atteunating material with linear attenuation coefficient $\\mu(x,y,z)$ at the energy of the emissions. The probability of the emission reaching the detector is the probability of it not being attenuated\n",
    "\n",
    "$$p(x,y,z,\\theta) = e^{-\\int_l \\mu(\\vec{l}) \\cdot d\\vec{l}}$$\n",
    "\n",
    "where the path $l$ forms a perpendicular line from voxel $(x,y,z)$ to the detector at angle $\\theta$\n",
    "\n",
    "* Thought: If an emission at (2,3,1) was going to yield 5 counts per second under no atteunation at angle $20^{\\circ}$, and $p(2,3,1,20^{\\circ}) = 0.2$, then it would instead receive 1 count per second under attenuation.\n",
    "\n",
    "* Thought: Our $g=Hf$ above is implemented using $H = \\sum_{\\theta} P(\\theta) \\otimes \\hat{\\theta}$. Maybe we can adjust for the probabilities before projecting to get $H = \\sum_{\\theta} P(\\theta) A(\\theta) \\otimes \\hat{\\theta}$ where $A(\\theta)$ is related to $p(x,y,z,\\theta)$. In other words, $A(\\theta)$ is used to adjust the voxel value of each voxel in $f$ before forward projecting\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evaluating $p(x,y,z,\\theta)$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First lets make a CT object, which will be a small cylinder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = torch.linspace(-1,1,128)\n",
    "xv, yv, zv = torch.meshgrid(x,x,x, indexing='ij')\n",
    "mu = (xv**2 + 0.9*zv**2 < 0.3) * (torch.abs(yv)<0.6)\n",
    "mu = mu.to(torch.float).unsqueeze(dim=0) * 0.1 #cm^-1"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the concept of linear attenuation coefficient, we need to specify the voxel dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = 0.3 #cm"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute $p(x,y,z,\\theta)$, we simply need to rotate the detector to angle $\\theta$ and then compute the integral\n",
    "\n",
    "* Example: $10^{\\circ}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu10 = rotate_detector_z(mu, angle=10)\n",
    "\n",
    "def rev_cumsum(x: torch.Tensor):\n",
    "    return torch.cumsum(x.flip(dims=(1,)), dim=1).flip(dims=(1,)) - x/2\n",
    "\n",
    "p = torch.exp(-rev_cumsum(mu10 * dx))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reverse cumulative sum is the opposite of the cumulative sum in that it starts at the RHS of the array as opposed to the LHS (required because the detector is at the RHS of the array).\n",
    "\n",
    "Then we simply multiply the object by this tensor $p$ and project to get our attenuation-corrected projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "object_10 = rotate_detector_z(obj,10) * p\n",
    "projection_10 = object_10.sum(axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7f24c35a0d00>"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABD9UlEQVR4nO3dfXScdZ3//9dk7pLmPilNGmghsl0LgoBUaoWzq0u+FmURlqoLp65d5NBdbJXSXYHuWlxUrLArcopI1bOLeha84RxBYb+ypxYoP9ZSSguuCJQqtS1tk97kZnLTTCaZ6/cHX2fm/Zlkrk4zaa60z8c5PSefua6ZXLnmpp/5vK7P+xPyPM8TAABAgJRN9gEAAAC46KAAAIDAoYMCAAAChw4KAAAIHDooAAAgcOigAACAwKGDAgAAAocOCgAACJzIZB/AsUin09q3b5+qq6sVCoUm+3AAAMBR8DxPvb29amlpUVlZ4TGSKdlB2bdvn2bNmjXZhwEAAI7Bnj17dNpppxXcZ0p2UKqrqyVJl+gjiig6yUcDAACOxrBSek7/N/P/eCFTsoPyx1gnoqgiITooAABMCf9v9b+juTyDi2QBAEDg0EEBAACBQwcFAAAEDh0UAAAQOHRQAABA4NBBAQAAgUMHBQAABA4dFAAAEDh0UAAAQODQQQEAAIFDBwUAAAQOHRQAABA4dFAAAEDg0EEBAACBQwcFAAAEDh0UAAAQOEV3UJ599lldccUVamlpUSgU0mOPPZbZlkqldOutt+rcc89VZWWlWlpa9KlPfUr79u0zj9HZ2anFixerpqZGdXV1uv7669XX1zfuPwYAAJwYiu6g9Pf367zzztP999+ft21gYEDbtm3T6tWrtW3bNv30pz/V9u3b9dGPftTst3jxYv32t7/V+vXr9cQTT+jZZ5/V0qVLj/2vAAAAJ5SQ53neMd85FNKjjz6qq666asx9tmzZoosuuki7du3S7Nmz9dprr+nss8/Wli1bNG/ePEnSk08+qY985CN666231NLS4vt7E4mEamtr9QFdqUgoeqyHDwAAjqNhL6Vn9DP19PSopqam4L4Tfg1KT0+PQqGQ6urqJEmbNm1SXV1dpnMiSW1tbSorK9PmzZtHfYxkMqlEImH+AQCAE9eEdlAGBwd166236tprr830lNrb2zVjxgyzXyQSUUNDg9rb20d9nDVr1qi2tjbzb9asWRN52AAAYJJNWAcllUrpE5/4hDzP0wMPPDCux1q1apV6enoy//bs2VOiowQAAEEUmYgH/WPnZNeuXXrqqadMztTc3KwDBw6Y/YeHh9XZ2anm5uZRHy8ejysej0/EoQIAgAAq+QjKHzsnO3bs0C9/+Us1Njaa7QsWLFB3d7e2bt2aue2pp55SOp3W/PnzS304AABgCip6BKWvr0+/+93vMu2dO3fq5ZdfVkNDg2bOnKmPfexj2rZtm5544gmNjIxkritpaGhQLBbTWWedpcsuu0w33HCD1q1bp1QqpeXLl+uaa645qhk8AADgxFf0NONnnnlGH/zgB/NuX7Jkif7lX/5Fra2to97v6aef1gc+8AFJbxdqW758uR5//HGVlZVp0aJFWrt2raqqqo7qGJhmDADA1FPMNONx1UGZLHRQAACYegJVBwUAAKBYdFAAAEDg0EEBAACBQwcFAAAEDh0UAAAQOHRQAABA4NBBAQAAgUMHBQAABA4dFAAAEDh0UAAAQODQQQEAAIFDBwUAAAQOHRQAABA4dFAAAEDg0EEBAACBQwcFAAAEDh0UAAAQOHRQAABA4NBBAQAAgROZ7AMAkBWKRE07fEqjaXuNtaY9XD/NtJP19v7J2nDOzyGzLVVtf/dQjW0PV4/YG6qHTbOyZtC0T6nqM+1ZVd2Zn8+Ydths+9Py/abdGLH3dR0erjLtNwZnmvYfBux52tNXl/n5YJ+9b3+i3D54r/0YjPSGTTuWsLtHe2073uM57ex5i3el7GN3DZh26HCPaY8ctOfJG7b3B04mjKAAAIDAoYMCAAACh4gHKLG8mKaxPvOzN73ObMuPaGK2XRt22jamGXJimpQT06Sq09lGzZDZNq3aRjQt1f2mnRvRSNIZFU5MU9Fu2nNitv3OaPb31ZfZv7N4nbZZbdtdaRudbE9lz+OOoWaz7Y0jtv2HI2PHQ5J0qLfStHt6nYgoYZ/vaE5kFE3Yj9hYb4Vpx3sanPZsu3+Xfc7yIqJD3aY9crgr8zPxEKY6RlAAAEDg0EEBAACBQ8QD+AiFbcwSnu7MrGmoM+3hRhtnDNZl44ZknRPZ1BQ7syZt2l6NHcavqE6a9sya7OyYUyvtjJF3TDtk2meWHzDtufF9pj0nah+7scxGH/kfJ8fv48WNkN4Xz/3Z/p2HK3eZ9o5U3LRfT7aY9u8HZ5j2mwPTTXtvv51ZdTCRnTU00Gsf+4gTB0V67XfEmBMJRZ0ZRvGEExF119n7d2ePPXLYiYM6u0175JAzY2jEmbUFTDJGUAAAQOAwgoKTTt6ISKO9UNHLuahVkkacC1kHnAtZB51RkSFnVCR3FMQdIUnVuCMittZIhXMha5PPhaynT7MXj/5JeUfm53fGbO2ROVH7DXtG2NYLmcwRkYnkjvw02kEOvS9uR5IOVL5p2jtq7Oth+5CtyfK7wabMz7sG7GvLvQD3sHMBbr9zAW7IHVFJ2O+U7ghL7ghMzBltKXdHW7rsSFE4r0ZLl2mPHLavLUZcMNEYQQEAAIFDBwUAAATOiTFmCzgKxjjuRa0Ndph9qNgIp4haJF6tjXDKq+yFp9NrbMn3UyttnfXWSnvBZ26EIxWOcfIjHLeN0bjnbYZ9OWhO9PemvSPnOdhe7sQ/05pMe2elc8FttX3xHKqyv3uwxuZRqR77EZ6qzn7nHKp2o0YnDnJe1+XO6z7mRJsR530j96LbnAiI+AelwAgKAAAIHDooAAAgcIh4MDWFbN86Mt3OlpAzEye3pPxQgx0mL7qcfJEzccpzapM0OrNwTquytUlOd1f9zSsn70Q6URsJEeMcf4UiIDf+2e48fzvKbeTzxjRbhn9Xpa2581afrblyuMrGk4M5dVdS1U7848wA8ouA4s77Il5MBOTMABo+5CxV4Nn3DDAaRlAAAEDg0EEBAACBQ8SDqWEckY4kDeWWm3cjHb9ZOeOIdCQb67Q4s3LcSMctN+9GOnN8Ix0Eifv8pNU3xp5vG/EKf2dMe6GC23NfTYPOtlTex73P99O83+VMX0rHNBa/En9EPjgajKAAAIDAoYMCAAACh4gHweAX4fgVV6uzw815qwbX5hawsg+dV1jNJ8KJOxHOdGdmjhvj5BZXy4tw4naWztxor2k3MyvnhOI+n815Rd92mPbrzutjR7md5fP7SrvScm7ht31u0bdq+55JukXfEmMXfZPyZ/0ka+zKzPGc91y8bnxF34iAIDGCAgAAAogOCgAACBwiHkwON9JpsLNw8iKd+mOPdCQb6/hGOtWFI52GKrssfd7MnEo7PJ0b6xQf6eBkkv/8946631jSBb5zujOAnBBFSc+2i571U3C7M+PH+V15s37SdofhTlv4jcjn5MAICgAACBw6KAAAIHCIeHD85MQ6eZHOdLfQmhPp1BeOdIZqxo50JFtszTfSqSkc6eStn+NEOnPKnfVzcmKdOVEbBzWHnQMFcriRz4gSY+w5fnmRj9POi3wKFJUL+RScy4t8HL6F3oh8TgpFj6A8++yzuuKKK9TS0qJQKKTHHnvMbPc8T7fffrtmzpypiooKtbW1accOO3Wus7NTixcvVk1Njerq6nT99derr69whUUAAHDyKLqD0t/fr/POO0/333//qNvvvvturV27VuvWrdPmzZtVWVmphQsXanAwW3h58eLF+u1vf6v169friSee0LPPPqulS5ce+18BAABOKCHP8zz/3ca4cyikRx99VFdddZWkt0dPWlpa9A//8A/6x3/8R0lST0+Pmpqa9L3vfU/XXHONXnvtNZ199tnasmWL5s2bJ0l68skn9ZGPfERvvfWWWlpafH9vIpFQbW2tPqArFQlFfffHcVLEzJzhRjt0nTcrp764CGfImZkz7MQ46Zxia7HqIbOt0S20VuUWWrPr5bzDKbY2N77PtN0Y59QpEuMk0kdM+2B6xLTbh22xrX3D9vndl7LtvUPZdsegPQedSftY3ckKeyxHyk37yKB9nw87bQ0637Vi2ec/Wpkymxpr7PP9znr7fP5F/Wum/VeVb5l2VZk9tiDZO5Kd9bPDma72etJ+tr456BR562807X199v6He53Ytde+Z8tyCr1FEvb5iDnJVMyZnBRL2PdrvMu+9mLd9j0bOeyMuOcUeiP+CbZhL6Vn9DP19PSopqam4L4lvUh2586dam9vV1tbW+a22tpazZ8/X5s2bZIkbdq0SXV1dZnOiSS1tbWprKxMmzdvLuXhAACAKaqkF8m2t799MWBTU5O5vampKbOtvb1dM2bYnnskElFDQ0NmH1cymVQymb1kK5GYuAvFAADA5JsSs3jWrFmjO+64Y7IPA64ii62N5Kyf40Y6Q0XOyikm0pFsrOMX6cyeZuczuJHOnJgzS+c4RjpuDNOZtn9n+4iNSvbmxDBuBLMnaYf09w/ak9o+YP+OQ/12iL+v30Ydw86Qf7gv+5xGe22RsFCRo+5+szqcGmRSKHebPa7DYft3/E/kFNP+/6LvMu0vxp2DnebED5U2fjilNhs/nFXfYbb9n7rfmvZfVtrXVkXIro9TLPvas6/LEZ+ZNSN5J7Gww0479ywMj7PIWyjvwoOYs90+h7mfHsz4OXGUNOJpbn57IauODvum7OjoyGxrbm7WgQP2TTk8PKzOzs7MPq5Vq1app6cn82/Pnj2lPGwAABAwJR1BaW1tVXNzszZs2KDzzz9f0ttxzObNm3XjjTdKkhYsWKDu7m5t3bpVF154oSTpqaeeUjqd1vz580d93Hg8rnh8fN8sUALOiEm4vtZuLzBiIknJ+uxz6I6YFCpNL9k6JtIoIybVzrda50LY3FomfiMmcypsB9sdMZkb6zbt8YyY9KUHTftQ2l7Q6Y6ItA9PN+09KTsKsjfpXKg6WJe9rzMicnjAPj+9fc6ISJ+9EDXcZz8uIs6oSNQ+BcdV3jfunLY7LpD3BXrY7lHmFADxjtjXqtdnX6vDMXue9h7MPmd7p9lVuZ+tOtO076u1V4ueU7/ftC+te9W0/9J5rcYLTBLIe106r9vx8gqUzh9yttlxvqOpk+K3PT5my1kgWmHnCR/psnWMGFEJrqI7KH19ffrd736Xae/cuVMvv/yyGhoaNHv2bK1YsUJf+cpXNGfOHLW2tmr16tVqaWnJzPQ566yzdNlll+mGG27QunXrlEqltHz5cl1zzTVHNYMHAACc+IruoLz44ov64Ac/mGmvXLlSkrRkyRJ973vf0y233KL+/n4tXbpU3d3duuSSS/Tkk0+qvDz77eyhhx7S8uXLdemll6qsrEyLFi3S2rVrS/DnAACAE8G46qBMFuqgTI6wcxFsqNEOX7uRzlCdHYbNrW2SdC+CdVcYdiOeWp9Ix6c8/anV2WFdN9I507kI9p1xO8w+J9pt2rMjxUU6uTGOf4RjYzO/CGf/oN2//Yg9toP92Xoz441wyiYxwjmeCl1wO9p2z8kUvEj2IzXtVHRPx3wuuK2y0eQMnwjoQ3WvmPbl07Kv82io8PfP3cP2sXek6kx7e3Kmaf/eqZuye8C+//f2Zl+LnX22zs1Qwn4WlPXakxbtsZ8H0bw6KbYdL1A3JdZtPwvCnfaieO+wff+PuBfRYkJNWh0UAACAUqCDAgAAAmdK1EHB5MiLdJx2MZGOZGMd30inxiaP+bN0Ckc6MwvM1Cl1pOM3E+fgSDZa2TvslBN3apO8NeQMmx+pM203wnFn4iScGCfVl80YcuuSSFLcKUd+skQ4fgrNCHp7h8L7eyPZHULO9JWyIXvO0057yJkx9Fa/zYgOJuwSEa9229IMv8yJgBY68c9lFc57Iu913a3xSBeoo9LpnKMhZxZOyi9Xy1Pou7V9bHfd5LBzVYM764fIJzgYQQEAAIFDBwUAAAQOEQ8yfCOdvBWInUjHLb5WoFy9b6STV6reiXSqC0c6ZzgrEOfGOnPibqn6btM+LeKUdHcinE4nwukYsbFKoRhnIiMcaZQYpzf7HJS51bJwTIqJgArFP1J+BBRK2feMl7Tt5IB9fnf32fdgR0/29fNq7djxjyR9qO43pt1WYV/XbuST9o1djl6n0x5y2qm8/5rcqVPuI5YV2OYT+fi0iXwmDyMoAAAgcOigAACAwCHiARx9aRsndTsrBruRTvuILZZWaGaOW1jtwKCNzcYb6UR6nZk5xDrHXW6s4zMBSHJmTpU5e6SdRyhzvlO6q8gkc+KMdhUugpXPRj5zY8kx9gOOD0ZQAABA4NBBAQAAgUPEcxIL19WZdtGF2NxZO7Vjz9qRpFRN7s926Hqk2mYRcWfWTl3VEdMuVIhNyi/GljtzpyVi7+sOk29P2bdF+4izHo4T4bjr47w1aNu5M3MO5UU4di0eN8IJ+RVXI8IJtPEWfQs5s348v1k/Q9n9B50ZP39wXlvtPTYC+k2tXXvnXGfWz3lVu037HXH7HstVqGjbaNu7ne3JvJk4PrN6DL/v3c6sHp/CbeG081nV3e3z+CgVRlAAAEDg0EEBAACBQ8RzksmNdULTbdGwvEin3s4gKaYQm2Qjnbfb2aHSvEinpnCkc2p1j2m7kc6cig7bdoqx5cY6lSE7deJw2g75urNyxhPpSDbWKTbSiRLpnNDyCrm52907uGv75M36yd0mZ5v9uB90truzfkJ5+ZTDTkArGPkUq9tp588nKiLy8fy+h9vPuaILuRH5TBhGUAAAQODQQQEAAIFDxHOCy5up05iNI4qepeNGOnkRjtsee6aOX6TjN0vnHRUHTduNdJrDvaadG+u4kc6+4TrTLmWkI9lYxzfSodDaSa3oyKdAoTd3dlpe5BPyiXxChQu9lYWc35DzsnfjnhGf78J+s366nXbhyMdvzSCf7+Wez9o9BWb9EPeUFiMoAAAgcOigAACAwCHiwQkhkbazY6LOTJ1BL7uge7GzdPYlC6+fU0zxNd9Ix13xHic138jHzXFyFJrh8/Z2OdvtfwdHnO1+a/scGsq+L94x7ZDZ1hDpL3hfYDSMoAAAgMBhBOUEU+iiWEkaacx+y8m7KLbevhzci2JTfnVOqguXr49VZ0cx/C6KPaPysGm7F8W2RLtMuzHcZ9qVoSHTzh01KXbEpP2I/UMP9NsRlN5+Z8Xh/qhp546aMGKC8ZgqIypl7pFNs80zyu37u1hdzkW15t3uFa6Rkl/exf3LC/+3WKhOCjVSSosRFAAAEDh0UAAAQOAQ8WBKGEjbgdXOYRuzpJ1y1t0hGy+159Q6KTbSOdRvL4J1I50ht7ZJr31b5cY6RDoopUKRT6G4RypB5OOUG9nvZd9H/UM2PnYvLG+ZZpeuqArbSBaQGEEBAAABRAcFAAAEDhHPFOc7a6feDq3mztwZqi1u1s6QTULyZ+3UjD1rR5LqqwcyP+eXsrezcprjdgi4IWzrKFSHbWHuqBPpHBy2Mc3+obrMz8VGOj39tq5JMZGORKyD4yc38ilmho80/sgn9x3puxKyw4183M8DP7l754dFRax8LMn9S0Npn1k9OSeKlY9LixEUAAAQOHRQAABA4BDxTDHhWhtP+EY6DfZq+txYZ7DOJ9IZRyE2SaqrGjDt3FjHHcL9k2kdpl1sIbYDI87BAwiMyphdf3hGuX0/n15hVyuvDbtl4QrLXWC429k2lJc2Of/t+aykXEwht0JF3KT8lZBHenqEsTGCAgAAAocOCgAACBwiHgRC34gtfnYwZPOllLO+RrkT8bizdt4aarDtwbrMz6WetRPpY9YOgmc86/ZI45vV47duj1vI7VDSvgebK+wsPwq5nZwYQQEAAIFDBwUAAAQOEQ8mxZF01LS7hu0Q74iztk7S2d8tzNaRsrObciMdycY6pY50GH3GVHA8I5+8dXucX5ZK2fktfYP2PXe43L5HG6fZQo0Rv4PFCYERFAAAEDh0UAAAQOAQ8QRcXmG26XZ2iluYLVU/dmE2SUrWZvukvoXZanwKs9UULszWUj32ejvuWjv1ETuEWxe2j+U3awdAgDm10KLREdOuKrfvbzfS8ZvVU1Zg7Z9upz3etXpC3tiF2/yWIPJdq4fCbQYjKAAAIHDooAAAgMAh4sGESTtrXOTO3Cn1rJ39Sds+cMTmV7kzd3oGnFk7/czawclnImf1hJyUxAvZMOOIbBTNrB6MhhEUAAAQOHRQAABA4BDxBEzerJ3GetP2m7WTrBt71o5kZ+rkzdqpLjxrJ1pdeNZOc1WvaTdV2HZdNLtCR42znHp1eND+riIjnf0D9o850G/PU26sQyE2IF9JIx8n4kk734Xdxx5WYX0+23MjH/dzyNXlFo1zd/CKm9Vjv+cX919q2LMHw6weq+QjKCMjI1q9erVaW1tVUVGhM888U1/+8pfl5TwRnufp9ttv18yZM1VRUaG2tjbt2LGj1IcCAACmqJJ3UO666y498MAD+uY3v6nXXntNd911l+6++27dd999mX3uvvturV27VuvWrdPmzZtVWVmphQsXanBwsMAjAwCAk0XJI55f/epXuvLKK3X55ZdLks444wz98Ic/1AsvvCDp7dGTe++9V1/4whd05ZVXSpJ+8IMfqKmpSY899piuueaaUh8SAACYYko+gvL+979fGzZs0BtvvCFJ+vWvf63nnntOH/7whyVJO3fuVHt7u9ra2jL3qa2t1fz587Vp06ZRHzOZTCqRSJh/AADgxFXyEZTbbrtNiURCc+fOVTgc1sjIiO68804tXrxYktTe3i5JampqMvdramrKbHOtWbNGd9xxR6kPFSU2OGJfTj1Dtt7ISE5dlGTY7tsXthf7RkO2FHaHc0Wve1HswYGxL4qV7IWxeRfF9nJRLOAaz0Wzbo0UV95Fs07NpFTKbh+OOBf/O3VSeo9kPz+qypMFfzemjpKPoPzkJz/RQw89pIcffljbtm3T97//ff3bv/2bvv/97x/zY65atUo9PT2Zf3v27CnhEQMAgKAp+QjK5z//ed12222Za0nOPfdc7dq1S2vWrNGSJUvU3NwsSero6NDMmTMz9+vo6ND5558/6mPG43HF4/FRtwEAgBNPyTsoAwMDKitzhszDYaXTb48Htra2qrm5WRs2bMh0SBKJhDZv3qwbb7yx1IcTeL51Txpsyfa8uicFViuWpKECKxbn1T2pKVz3pL762OueSLb2SbF1TwCcwJw8KRS1+VEk7qykHrPt3NWQiy2D3+W08+qi+K12bA7dCSXyaqoUxmrHVsk7KFdccYXuvPNOzZ49W+9617v00ksv6Z577tGnP/1pSVIoFNKKFSv0la98RXPmzFFra6tWr16tlpYWXXXVVaU+HAAAMAWVvINy3333afXq1frMZz6jAwcOqKWlRX/3d3+n22+/PbPPLbfcov7+fi1dulTd3d265JJL9OSTT6q8vLzUhwMAAKagkndQqqurde+99+ree+8dc59QKKQvfelL+tKXvlTqXw8AAE4ALBYIAAAChw4KAAAIHFYzxjErpjCbZIuzlbowW3f/2IXZJFucjdWKgeIFunBbMpr5Obdom5RfuC3tuzoxgoIRFAAAEDh0UAAAQOAQ8UyCcHW2eppbmC1db6ML38JsdT6F2Zx6Z7nF2fwKs9VW20Jr4ynMJtnibG5htkPDzoEDOHkVWbgtntOudiKd+gpbYNKvcJtbCq3owm2GOwZQ+L/ckGf/bvfe4bQ99pHeXp3IGEEBAACBQwcFAAAEDh0UAAAQOHRQAABA4NBBAQAAgcMsHhw1tzBb19A0006l7eLgSaeYUm5xtuNZmE2yxdnC9iJ/AMegUOE2n4kypS/clsy+3wcHo2Zbz4BdhNYt3IbgYgQFAAAEDh0UAAAQOHRQAABA4NBBAQAAgUMHBQAABA4dFAAAEDh0UAAAQODQQQEAAIFDBwUAAAQOHRQAABA4lLo/DsLV1aYdOqUx83O63pZwH2qwZZmTtRGnbfuUQ/ahlapx2tW2HvVIdbbEfLR6yGyrrT5i2i1VCdOeXt5n2o2xAdOuCdv7V4cHTTtelsr8fMA9UAAYi1NXPxS1tfQj8eHMz/GcnyWp2ilt31jRX9Sv7nHaKc8t059d4iOUt83huWMC7n/B9vM/5mx17527uMhIb2/h3z0FMYICAAAChw4KAAAIHCIejMldvTiRsisIp53hykKrF0t2BWN39eKOIzarclcvdlckHeo/+tWLJVYwBiZabgrjudtKvbqxsz13dePclY0lVjeeyhhBAQAAgUMHBQAABA4dFAAAEDh0UAAAQODQQQEAAIFDBwUAAAQOHRQAABA4dFAAAEDg0EEBAACBQwcFAAAEDh0UAAAQOHRQAABA4NBBAQAAgUMHBQAABA4dFAAAEDh0UAAAQODQQQEAAIFDBwUAAAQOHRQAABA4dFAAAEDg0EEBAACBMyEdlL179+qTn/ykGhsbVVFRoXPPPVcvvvhiZrvnebr99ts1c+ZMVVRUqK2tTTt27JiIQwEAAFNQyTsoXV1duvjiixWNRvWLX/xCr776qr7+9a+rvr4+s8/dd9+ttWvXat26ddq8ebMqKyu1cOFCDQ4OlvpwAADAFBQp9QPeddddmjVrlh588MHMba2trZmfPc/Tvffeqy984Qu68sorJUk/+MEP1NTUpMcee0zXXHNNqQ8JAABMMSUfQfn5z3+uefPm6eMf/7hmzJihCy64QN/97ncz23fu3Kn29na1tbVlbqutrdX8+fO1adOmUR8zmUwqkUiYfwAA4MRV8g7Km2++qQceeEBz5szRf//3f+vGG2/U5z73OX3/+9+XJLW3t0uSmpqazP2ampoy21xr1qxRbW1t5t+sWbNKfdgAACBASt5BSafTes973qOvfvWruuCCC7R06VLdcMMNWrdu3TE/5qpVq9TT05P5t2fPnhIeMQAACJqSd1Bmzpyps88+29x21llnaffu3ZKk5uZmSVJHR4fZp6OjI7PNFY/HVVNTY/4BAIATV8kvkr344ou1fft2c9sbb7yh008/XdLbF8w2Nzdrw4YNOv/88yVJiURCmzdv1o033ljqw8E4lIeHTbsmesS0G2MDdnvYbq8O21lZ8bJU5udo2UjB3532QgW3dzvbhzy7fTjvpZ3ti4eTBR8awDEwb0nn7es5X4W9sG2no/YNnI7ZthdP2ztU2M+PaDz7WRWP28+t6nL7hm+s6Fch7f3VBbfj+Cl5B+Xmm2/W+9//fn31q1/VJz7xCb3wwgv6zne+o+985zuSpFAopBUrVugrX/mK5syZo9bWVq1evVotLS266qqrSn04AABgCip5B+W9732vHn30Ua1atUpf+tKX1NraqnvvvVeLFy/O7HPLLbeov79fS5cuVXd3ty655BI9+eSTKi8vL/XhAACAKSjkeZ7nv1uwJBIJ1dbW6gO6UpFQdLIPp2jh6uwQYuiURrMtXV9l2kMNttOWrLV9ymSdHTsdckYnh2ptO1WdfbpHauxQaLR6yLRrq21k01DhRDoxO3RaFbHtmoiNeCrC2cePhuwQbceQva5o/4BtHxyw56W7v8K0h/piph3qtecp0kfEA0ykCY143HbURj6hnHY4Zj9bYjH7OVflRD5p52B7eu1nS6rXfraEE/azJdqbvX+sx2xSrNe24932uOM99thinfYzs6yrz7S9g4dNe6TX+QVTwLCX0jP6mXp6enyvJ2UtHgAAEDh0UAAAQODQQQEAAIFDBwUAAAQOHRQAABA4dFAAAEDg0EEBAACBQwcFAAAEDh0UAAAQOHRQAABA4JR8LR6cuNzVjeud1YzrnNWOC61uHA05ZfZ9Vjf20+20h5y2Xd3Y9sspfQ8UL2/B8Zz28Vy9WLIrGLurF9c7S3REQvax9/UVLreOycMICgAACBw6KAAAIHDooAAAgMChgwIAAAKHDgoAAAgcZvFMgpHe3szP4TLbR3R7jNGQe6m8y30K/fqcuY9n75ty9uzJe2R7pf3giL1/T6zCtKuig7Ydzs6tqQjbeTbR0Phm8QA4gThThLyU/Vwb9uxnj5ezfzpt7zvitNOy7Z5e+7mV6o2Zdjhhf1e0194/lvNBGes1mxTvtjOG4j3O7MUuO+OorKvPtL3DXaad+3/HyYARFAAAEDh0UAAAQOAQ8eCYuYXbamO2MFuhwm25Rduk/MJtZSEbJxWr22nnBkrD7sveqSoVdqu8AShYmE2yb6OJLswWi9l2VXn2Tds4rd9scwuztfdVC1MDIygAACBw6KAAAIDAoYMCAAAChw4KAAAIHDooAAAgcJjFM8lGemw5NOfi97y2r5BP4bbQmA35FW7rUmHFFG6rcqbKULgNOIkVW5gtbbfnFmrzfCYAdvVOM+2iC7Ml7OPlFmeL9xRXmC3caQuv5RVm63HLZZ5cGEEBAACBQwcFAAAEDhEPSqaYwm25Rduk0hduc4d5cwdK3TpseYXbeincBhRTmE2yxdmOZ2E2yRZnozDbiYMRFAAAEDh0UAAAQODQQQEAAIFDBwUAAAQOF8kGjG9dlJB75Zofn7oohvPYocJ1UbqLPZQiuBfNNkVLVw/AfaS8i2admgvq46JZnPjGc1GsJKUj2Qthi70oNjKOi2Ile2Gse1Fsd59P3ZPeIuueOO3c2ifxbp+6J119pk3dk8IYQQEAAIFDBwUAAAQOEQ8mjFu7pKIsGxLVR+wQrduuCw+YdnmouFL4I854dTqvrH9Wj1MzZchp59VJIfLBCaCUkY4keTmxjudEPG6kU1Flo4/qCtt2I5zmCpuruEtl7B6oF048jKAAAIDAoYMCAAACh4gn4Ma72nHeMG7OU+6F/PqnhVc7dpONbp9Hy3UkHTXtrmilaRdbCh9AgDmJTyplP7n6QnZmjZ9C5ezdWTtDCZ9ZO4nCs3aidsHhvBWLYzkrFvvO2jnUadrM2imMERQAABA4dFAAAEDgEPFgUuTO6JEmflZPMXwLuTGrB1PQRM7akaR0PKc9zUawzNrBsWAEBQAABA4dFAAAEDhEPAiEKmeWzikRO8TbGLZXw1dOYMTjKjbyCfXafn+Zu4gRcBwc10hHMrGOG+k019r382mV9l11ZuVB0651ZvG9eeQU4eTDCAoAAAicCe+gfO1rX1MoFNKKFSsytw0ODmrZsmVqbGxUVVWVFi1apI6Ojok+FAAAMEVMaMSzZcsWffvb39a73/1uc/vNN9+s//qv/9Ijjzyi2tpaLV++XFdffbX+53/+ZyIP54SQV7gtZMdt3cJthcsfuU+/T3815FO4zdnc7Xf3AgbK7ZF3hqtM2y3c5s7qARAc/UNx0z7gvJ8jZTai9Zu1s7+vxrRzi7MN9foUZustrjBbeffYhdkkKdaZjbPyCrMd7jJtCrMVZ8JGUPr6+rR48WJ997vfVX199sXV09Ojf//3f9c999yjv/iLv9CFF16oBx98UL/61a/0/PPPT9ThAACAKWTCOijLli3T5Zdfrra2NnP71q1blUqlzO1z587V7NmztWnTplEfK5lMKpFImH8AAODENSERz49+9CNt27ZNW7ZsydvW3t6uWCymuro6c3tTU5Pa29tHfbw1a9bojjvumIhDnfJGurtN22+tnliBLqlXVmTkU+RaPV06emXOWhuKj77fH4XDdn93FtBIzt+Szl+gqCh+s3pScoeUs7+bGT2YSOalfRxn7UhSeWX2neAWYquvsIUWW6bZd1GxkU5Xr7PeTk6sE074RDrOGzhvrZ2ET6TTPfZ6O3mRjvP5jOKUfARlz549uummm/TQQw+pvLy8JI+5atUq9fT0ZP7t2bOnJI8LAACCqeQdlK1bt+rAgQN6z3veo0gkokgkoo0bN2rt2rWKRCJqamrS0NCQup2eZUdHh5qbm0d9zHg8rpqaGvMPAACcuEoe8Vx66aX6zW9+Y2677rrrNHfuXN16662aNWuWotGoNmzYoEWLFkmStm/frt27d2vBggWlPhycIKaV2SHghoi9Wt6vkFs45AxPj0PaGTt3RogLRj5RirihhAoVYyt5pFNReH2d3GJsbiG2M6YdNu0GZ62tPww2CnCVvINSXV2tc845x9xWWVmpxsbGzO3XX3+9Vq5cqYaGBtXU1Oizn/2sFixYoPe9732lPhwAADAFTUqp+2984xsqKyvTokWLlEwmtXDhQn3rW9+ajEMBAAABdFw6KM8884xpl5eX6/7779f9999/PH79SaWYWT2Fi7hJfoXc8ifDODeExp7V013CyGVUTuTTHD5+BZISzolJ5bTdRIfIB8UoZn2dkkc61YXX12mpzLaLjXT+0G/b+YXYKky7UDE2d9ZO1KcQmztrJ97lM2vn8NjF2Ji1U1qsxQMAAAKH1YxxQqgps6ufNpbZugvVZdlvRWXjHL1JO1cf+tVVyf0Cl3LGrRhRQSHjWZF4okdM3l2/z7TPrXwr83NL1NYD2ZEcfYYmUAgjKAAAIHDooAAAgMAh4jnBFbpo1l0JOe+iWd/lh4+9FH7S2dLt95t8YpmwCpfGLwtlh6cby5zfHvH77aXjriLlG/kknMhnWDiB+UY6brvAhbB5kU6scKRTXmUr+LiRzjl1dimS3EhHsrHOm8kZZtubR04x7d0DDabtd1FsMmHf0HkrFCdCOT+bTXmrFeddFNttV1LOuyi2017gSzn744cRFAAAEDh0UAAAQOAQ8ZxkcocjfVc+zru3z+KPoWIiH/vSKzbyGY8WZ6VjN/JJH8+aKU6byOfkUnSk43xiF5qpkxfpuKsPFxnpnFe127TdmTq5sc6OI01mmxvp7O2tNe3xRDpvt7M/x/zqnLiRTteg/V1upHOo07SJdI4fRlAAAEDg0EEBAACBQ8RzEsub4VNmh03zIp+8WT3OVJk8hfq/Y8/wkaQhZ3O3+8g+s3qKKcb2Tmeo+p1ROxTeFO4wbbds/qnO/WfG7HD2W+XZ9v4KO7TdPq3atA9WVpl2otLGaslKp8R3n32WIm6hNyKgSTXuWTlRJ8KJOtvjzuy1adn4orzAasOSdG79ftP+UJ1dhX5u7JBpv5mqM223+NrvB7MRT7GzdIZ6jz3SkWysE+8Z5yydTmbpBAUjKAAAIHDooAAAgMAh4gEcVWV2uLksL44aVCEjUb8Cd1nuOj7FrOsj5c/6ceVGPsQ9x4d5Co9jpCNJ8ZxYp9hIp63CTn/pGJngFccBH4ygAACAwKGDAgAAAoeIBxkjztXrxRdyc2b15KUVZQW22RtS7qweZ+9OlU6Z3KHsbtOaHbEzbarK7MyahjIb+ZwStgWuTs1Z68ctbrUvVm/ab1U4Ba0q6kzbnfVzuKrStBN9zqyfquyzlDfjxy36ZtMCjKGYmTluhONF3PVx7PZiIhxJaqq1scw5OTHOwrpXzLbLKmzkEw3Z99juYfs63lHErB1J+kN/Y+Znd9ZOZ+8003Zn7ZQlfGbtOMXXChVj8521c7jPtPNm7ThtTB5GUAAAQODQQQEAAIFDxIMx+UY+TuE238jHKGbdHikVsr99yHnsTp+JM8UUbsvXbVp+kY/bnp4TAbnxT3PEFn1zI6A9sUbT3ltuI6FiCr/15hV9s1NEwn324yDSa0/qyRIBjbe4Wm6MkxfhxApHODFnfZwZBSIcSfqQE+NcPi37enIjHPfjfvewfWw30tmenGnabqRTqBhbZ59PpNNrT5pvpONMX8tbX6crex6LLsRGpBNYjKAAAIDAoYMCAAACh4gHR22ky8YRxc7y8UIFIp+QX1/Zbk85W/MiH59HK/ibfOKgcKjbtE8NV4++4/+TG/kUin8kqdknApoVtZHOnnjhCCh3FpDfDKDevBlAJ2YEVMoIR/KJccYZ4Vxa96pp/+U0+8qOh5xKbgU+0veOFI50dgwVnqXjRjp7e+1rMTfWGUr4RDo99v1czCwdSYoVmKmTF+kctufM/RxDcDGCAgAAAocOCgAACBwiHhw9zw6zukOlESemcSOfeIFWvsJ9Z8/5Xe4yM0MhO+6eO8gbGteMnlHEuk3TL/LJ5TcD6AznHZqIHjLtzriNhNor7DL2e4ezkc++lI1/9iRtPLR/0BbXah+wf8ehfhsJ9fU7kVCvPee5heGiThwUciazjFeh2KboyKbI9W9ilTa2OaU2WwjsrPoOs+3/1P3WtP+y8oBpVxSKQSVJbqQzNjfSeX2ozrTdSGfHkSbTdiOdfW7xtbyZOtnn34103KKAvpFOT+FIJ95VYKZOZ7fZlhfpeCV+8WHCMIICAAAChxEUHDvnm8iwU0/AfXHlfqeKhdyvvPbbd/6ivsVdRDtcoFT+YZ9HcoV9L5p1v5HZog3FjKj4qSmrcNp2uzvionjuV1Pna6p2m1YifcS0D6btt9b2YfuNed+wHZFxR2j2DmXbHYP2HHQm7WN1J+3flThiR2eODNqRg2GnrUHnRORcqBqttJdUN9bYiyjfWW9HMf6i/jXT/qvKt0zbHeUaH78Rk+LkjprsSNkRD3fE5E2fi2DdEZPDvXYEbcgZMcstV++OmLh1TNwRk5jvRbB2lMq9EDZ31MT9HGLEZOpiBAUAAAQOHRQAABA4RDwoHTfyOeyEKTlDr5FOGwdEptt2rNsZTq63w8nJOqf0fa/tayeddiqRvX+q1r7s91c7F9TWOFFGlR3q3lVlh8LfrLRD5XPK7YWrc3IuZJ0Tnbj4Z7z84qMz8z4t3HoSbvsP4z+oQChlpDM++bVMnBgn+afZnwdtpLOr375u3+obu46JJCWdWiahXvsCiBWoZRJ3Ix23NL0b4XTZCCfS5UQ4h2xsQ4xzcmAEBQAABA4dFAAAEDhEPDh+coZh/Wb8+L8w89dONr+qUOl8Z1vK+W22wsL4yub7CTsza5rDVWPsiZNR+0ifaedHOk4tk5xYp9SRTrSIWiZEOigFRlAAAEDg0EEBAACBQ8SDyeFX5K3MVmrLe6HmFXIrHPkU7ov7RD7O73IjH7/Vj8tUzHA0kc/JzI10Xk/ZWV5upOOuOJwb6+zrL1yaPtlbZKTjU2wttzx9XqTT7RPpOOXpiXQgMYICAAACiA4KAAAIHCIeBIMb+Ry0q/bqkA1WItPtDIVIo1Pord4OZ8frshFQst6JcJyVdt2ib7lF3iRpqMbef2+NHSo/XG2LzO2rzg6176ycbrb93inytqPcrgszN77ftOdEbQRABDS15M/Ksc/f68k5pu1GODv77evHjXFy18sZdCOchFNozW9WjhPpxBM2yox32TXEc2OcSNeAvfNhZ1aO834mwsFoGEEBAACBQwcFAAAEDhEPpgY3AnIjH2f3gi/sMncKkF3XR6G8KULuA5iWO+tn0NnbWZHIOZTCM4DCocJD32WykcEMIp9AOeAT6ewYajLtvFk5A42mXSjSkWys40Y6xRRak0aJdHr8ZubkxDpEOigBRlAAAEDg0EEBAACBQ8SDqcln1k/ILfzWmJ31E+mqM9tiDXaYPF5vZ+0M1tkIaKjGZ9ZPr71/sjr7Nttb68z4qbK/e2+NHcLfVWmH+H9XaSOBHeUdpv3OmDvrJzvsTvwzMfJjnOwMsu1DZ5ptvxu0z587K2evE+EcStjnbLDPmZnT48zMyXkt+hVWizkRTrnfejmdhYurjRzOxjjeiH0s4FiUfARlzZo1eu9736vq6mrNmDFDV111lbZv3272GRwc1LJly9TY2KiqqiotWrRIHR0dYzwiAAA42ZS8g7Jx40YtW7ZMzz//vNavX69UKqUPfehD6u/P9r5vvvlmPf7443rkkUe0ceNG7du3T1dffXWpDwUAAExRIc/zCk8jGKeDBw9qxowZ2rhxo/7sz/5MPT09OuWUU/Twww/rYx/7mCTp9ddf11lnnaVNmzbpfe97n+9jJhIJ1dbW6gO6UpFQdCIPHyegUNhGNuFGW/TNc4q+jThF34aKjYByRu2dpVWUqrFRlVdji19VVNs5QY3Vdph9VlW3aZ8+zc6W+JOcCKhQ/CMRAY2lUIQjSduHZpp2boyza8C+tvb01Zm2OwvnSG+5aRc9EyeR+3NxEU7YKa4Wcmbi5EY4EjEOjs2wl9Iz+pl6enpU40Targm/SLanp0eS1NDw9ht169atSqVSamtry+wzd+5czZ49W5s2bZrowwEAAFPAhF4km06ntWLFCl188cU655xzJEnt7e2KxWKqq6sz+zY1Nam9vX3Ux0kmk0omk5l2IpEYdT8AAHBimNAOyrJly/TKK6/oueeeG9fjrFmzRnfccUeJjgonO3doevjAQbuD03YjoWnT7cyaioY6+3iNTiSUuw6QEwcla5wZQNXuuj+2vbfaxjBv1dgI4eXqU037lJrWzM+nVvaYbe+YZmc+nZm3DtA+054TTZp2Y5mNJ6aKw2kbk+1I2ZkxrydbTPv3g+eY9psD7sybWtM+mDPz5kjeejg2ko44M8Aqne9eboSTVzzNjW1y18M57EQ27qybQ7aE4AiRDQJmwiKe5cuX64knntDTTz+t0047LXN7c3OzhoaG1N3dbfbv6OhQc3PzqI+1atUq9fT0ZP7t2bNnog4bAAAEQMlHUDzP02c/+1k9+uijeuaZZ9Ta2mq2X3jhhYpGo9qwYYMWLVokSdq+fbt2796tBQsWjPqY8Xhc8Xh81G3ARMsbcemwIw1y2mUR+y25Muei22nT6+xj+VyAm6x1Rlxq7XeKoWr7vkg5Kyvvzhlx2V1jR35+XW1HCqb7XIB7RoX9xv2nFTaSnROz7XdGs9/m68vs31lqXWk7WrA9lT2PO4bsF583jtgRkT8csefFvZD1kHMh64BzIaucUZFozqjINJ9aJPEen3Lybi0S90LWQ92mPZJzYevIcErAVFbyDsqyZcv08MMP62c/+5mqq6sz15XU1taqoqJCtbW1uv7667Vy5Uo1NDSopqZGn/3sZ7VgwYKjmsEDAABOfCXvoDzwwAOSpA984APm9gcffFB/+7d/K0n6xje+obKyMi1atEjJZFILFy7Ut771rVIfCgAAmKImvA7KRKAOCk5UISceCp9i4wev0V6Q6UZEyXp7/9yIKFlr67O4NVmGnJIEw9XORZPVtkZLZY2t0XJKla0XkhsRnTHNiYfKbU2Wxoi9r+vwsL04+I1BW3vkDwNjxzQH++x9+xNORNNrv6dFem2sFvO7cLVATBPvsjFLXkRz2F64PHLQniePmAYnmEDVQQEAACgWHRQAABA4rGYMBIg7pD+83yle6LRtaCOV+7SPp3bzs/2oeV6zSvzbbFQSzWm3uLtOIiqNAEePERQAABA4dFAAAEDg0EEBAACBQwcFAAAEDh0UAAAQOHRQAABA4NBBAQAAgUMHBQAABA4dFAAAEDh0UAAAQODQQQEAAIFDBwUAAAQOHRQAABA4dFAAAEDg0EEBAACBQwcFAAAEDh0UAAAQOHRQAABA4NBBAQAAgUMHBQAABA4dFAAAEDh0UAAAQODQQQEAAIFDBwUAAAQOHRQAABA4dFAAAEDg0EEBAACBQwcFAAAEDh0UAAAQOHRQAABA4NBBAQAAgUMHBQAABA4dFAAAEDh0UAAAQODQQQEAAIFDBwUAAAQOHRQAABA4dFAAAEDg0EEBAACBQwcFAAAEDh0UAAAQOHRQAABA4NBBAQAAgUMHBQAABM6kdlDuv/9+nXHGGSovL9f8+fP1wgsvTObhAACAgJi0DsqPf/xjrVy5Ul/84he1bds2nXfeeVq4cKEOHDgwWYcEAAACYtI6KPfcc49uuOEGXXfddTr77LO1bt06TZs2Tf/xH/8xWYcEAAACIjIZv3RoaEhbt27VqlWrMreVlZWpra1NmzZtyts/mUwqmUxm2j09PZKkYaUkb+KPFwAAjN+wUpIkz/P/z3tSOiiHDh3SyMiImpqazO1NTU16/fXX8/Zfs2aN7rjjjrzbn9P/nbBjBAAAE6O3t1e1tbUF95mUDkqxVq1apZUrV2ba3d3dOv3007V7927fPxBZiURCs2bN0p49e1RTUzPZhzMlcM6ODeeteJyzY8N5K95knjPP89Tb26uWlhbffSelgzJ9+nSFw2F1dHSY2zs6OtTc3Jy3fzweVzwez7u9traWF+QxqKmp4bwViXN2bDhvxeOcHRvOW/Em65wd7cDCpFwkG4vFdOGFF2rDhg2Z29LptDZs2KAFCxZMxiEBAIAAmbSIZ+XKlVqyZInmzZuniy66SPfee6/6+/t13XXXTdYhAQCAgJi0Dspf//Vf6+DBg7r99tvV3t6u888/X08++WTehbOjicfj+uIXvzhq7IOxcd6Kxzk7Npy34nHOjg3nrXhT5ZyFvKOZ6wMAAHAcsRYPAAAIHDooAAAgcOigAACAwKGDAgAAAmdKdlDuv/9+nXHGGSovL9f8+fP1wgsvTPYhBcaaNWv03ve+V9XV1ZoxY4auuuoqbd++3ewzODioZcuWqbGxUVVVVVq0aFFe0byT2de+9jWFQiGtWLEicxvnbHR79+7VJz/5STU2NqqiokLnnnuuXnzxxcx2z/N0++23a+bMmaqoqFBbW5t27NgxiUc8uUZGRrR69Wq1traqoqJCZ555pr785S+bdUk4Z9Kzzz6rK664Qi0tLQqFQnrsscfM9qM5R52dnVq8eLFqampUV1en66+/Xn19fcfxrzi+Cp2zVCqlW2+9Veeee64qKyvV0tKiT33qU9q3b595jKCdsynXQfnxj3+slStX6otf/KK2bdum8847TwsXLtSBAwcm+9ACYePGjVq2bJmef/55rV+/XqlUSh/60IfU39+f2efmm2/W448/rkceeUQbN27Uvn37dPXVV0/iUQfHli1b9O1vf1vvfve7ze2cs3xdXV26+OKLFY1G9Ytf/EKvvvqqvv71r6u+vj6zz9133621a9dq3bp12rx5syorK7Vw4UINDg5O4pFPnrvuuksPPPCAvvnNb+q1117TXXfdpbvvvlv33XdfZh/OmdTf36/zzjtP999//6jbj+YcLV68WL/97W+1fv16PfHEE3r22We1dOnS4/UnHHeFztnAwIC2bdum1atXa9u2bfrpT3+q7du366Mf/ajZL3DnzJtiLrroIm/ZsmWZ9sjIiNfS0uKtWbNmEo8quA4cOOBJ8jZu3Oh5nud1d3d70WjUe+SRRzL7vPbaa54kb9OmTZN1mIHQ29vrzZkzx1u/fr3353/+595NN93keR7nbCy33nqrd8kll4y5PZ1Oe83Nzd6//uu/Zm7r7u724vG498Mf/vB4HGLgXH755d6nP/1pc9vVV1/tLV682PM8ztloJHmPPvpopn005+jVV1/1JHlbtmzJ7POLX/zCC4VC3t69e4/bsU8W95yN5oUXXvAkebt27fI8L5jnbEqNoAwNDWnr1q1qa2vL3FZWVqa2tjZt2rRpEo8suHp6eiRJDQ0NkqStW7cqlUqZczh37lzNnj37pD+Hy5Yt0+WXX27OjcQ5G8vPf/5zzZs3Tx//+Mc1Y8YMXXDBBfrud7+b2b5z5061t7eb81ZbW6v58+eftOft/e9/vzZs2KA33nhDkvTrX/9azz33nD784Q9L4pwdjaM5R5s2bVJdXZ3mzZuX2aetrU1lZWXavHnzcT/mIOrp6VEoFFJdXZ2kYJ6zKbGa8R8dOnRIIyMjedVmm5qa9Prrr0/SUQVXOp3WihUrdPHFF+ucc86RJLW3tysWi2VelH/U1NSk9vb2STjKYPjRj36kbdu2acuWLXnbOGeje/PNN/XAAw9o5cqV+qd/+idt2bJFn/vc5xSLxbRkyZLMuRnt/XqynrfbbrtNiURCc+fOVTgc1sjIiO68804tXrxYkjhnR+FozlF7e7tmzJhhtkciETU0NHAe9fY1dbfeequuvfbazGKBQTxnU6qDguIsW7ZMr7zyip577rnJPpRA27Nnj2666SatX79e5eXlk304U0Y6nda8efP01a9+VZJ0wQUX6JVXXtG6deu0ZMmSST66YPrJT36ihx56SA8//LDe9a536eWXX9aKFSvU0tLCOcNxkUql9IlPfEKe5+mBBx6Y7MMpaEpFPNOnT1c4HM6bPdHR0aHm5uZJOqpgWr58uZ544gk9/fTTOu200zK3Nzc3a2hoSN3d3Wb/k/kcbt26VQcOHNB73vMeRSIRRSIRbdy4UWvXrlUkElFTUxPnbBQzZ87U2WefbW4766yztHv3bknKnBver1mf//znddttt+maa67Rueeeq7/5m7/RzTffrDVr1kjinB2NozlHzc3NeRMnhoeH1dnZeVKfxz92Tnbt2qX169dnRk+kYJ6zKdVBicViuvDCC7Vhw4bMbel0Whs2bNCCBQsm8ciCw/M8LV++XI8++qieeuoptba2mu0XXnihotGoOYfbt2/X7t27T9pzeOmll+o3v/mNXn755cy/efPmafHixZmfOWf5Lr744rwp7G+88YZOP/10SVJra6uam5vNeUskEtq8efNJe94GBgZUVmY/dsPhsNLptCTO2dE4mnO0YMECdXd3a+vWrZl9nnrqKaXTac2fP/+4H3MQ/LFzsmPHDv3yl79UY2Oj2R7IczYpl+aOw49+9CMvHo973/ve97xXX33VW7p0qVdXV+e1t7dP9qEFwo033ujV1tZ6zzzzjLd///7Mv4GBgcw+f//3f+/Nnj3be+qpp7wXX3zRW7BggbdgwYJJPOrgyZ3F43mcs9G88MILXiQS8e68805vx44d3kMPPeRNmzbN+8///M/MPl/72te8uro672c/+5n3v//7v96VV17ptba2ekeOHJnEI588S5Ys8U499VTviSee8Hbu3On99Kc/9aZPn+7dcsstmX04Z2/PqHvppZe8l156yZPk3XPPPd5LL72UmXFyNOfosssu8y644AJv8+bN3nPPPefNmTPHu/baayfrT5pwhc7Z0NCQ99GPftQ77bTTvJdfftn835BMJjOPEbRzNuU6KJ7neffdd583e/ZsLxaLeRdddJH3/PPPT/YhBYakUf89+OCDmX2OHDnifeYzn/Hq6+u9adOmeX/1V3/l7d+/f/IOOoDcDgrnbHSPP/64d84553jxeNybO3eu953vfMdsT6fT3urVq72mpiYvHo97l156qbd9+/ZJOtrJl0gkvJtuusmbPXu2V15e7r3jHe/w/vmf/9n8J8E587ynn3561M+xJUuWeJ53dOfo8OHD3rXXXutVVVV5NTU13nXXXef19vZOwl9zfBQ6Zzt37hzz/4ann3468xhBO2chz8spYQgAABAAU+oaFAAAcHKggwIAAAKHDgoAAAgcOigAACBw6KAAAIDAoYMCAAAChw4KAAAIHDooAAAgcOigAACAwKGDAgAAAocOCgAACBw6KAAAIHD+f9PeAzv2KL+gAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.pcolormesh(projection_10[0].T)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the points in the center are now darker because they're attenuated. To get the full image (collection of projections), we modify our $g=Hf$ loop above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = torch.linspace(-1,1,128)\n",
    "xv, yv, zv = torch.meshgrid(x,x,x, indexing='ij')\n",
    "obj = (xv**2 + 0.9*zv**2 < 0.5) * (torch.abs(yv)<0.8)\n",
    "obj = obj.to(torch.float).unsqueeze(dim=0)\n",
    "mu = (xv**2 + 0.9*zv**2 < 0.3) * (torch.abs(yv)<0.6)\n",
    "mu = mu.to(torch.float).unsqueeze(dim=0) * 0.1 #cm^-1\n",
    "\n",
    "angles = np.arange(0,360.,3)\n",
    "image = torch.zeros((1,len(angles),128,128))\n",
    "for i,angle in enumerate(angles):\n",
    "    mu_i = rotate_detector_z(mu, angle)\n",
    "    p_i = torch.exp(-rev_cumsum(mu_i * dx))\n",
    "    object_i = rotate_detector_z(obj,angle)\n",
    "    image[:,i] = (object_i*p_i).sum(axis=1)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot multiple projections:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAz8AAADFCAYAAAB6tEyhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACEb0lEQVR4nO39e5As2VXfj37X3jszq/pxzpkzj3M0aEYSSEhCBoGFNJ5r2T8eA4P4BcaYGzekq2uwL4HiEoa4tnA4LNtYhrCv/AM/JcsofoR/VviBhR3Xwj9jW/wkwQ9dOaQRCAFGL5AYpBGjc0bzOKcf9cjce6/7x35kVnV1d3Wf7q6q7vWJU6eqsjKrsnPnzp1rr7W+i5iZIQiCIAiCIAiCcM5Ri94BQRAEQRAEQRCEs0CMH0EQBEEQBEEQLgRi/AiCIAiCIAiCcCEQ40cQBEEQBEEQhAuBGD+CIAiCIAiCIFwIxPgRBEEQBEEQBOFCIMaPIAiCIAiCIAgXArPoHTgO3ns8+eST2NzcBBEtencuFMyM7e1t3H///VDqzmxnacfFIG24+kgbng+kHVcfacPzgbTj6nOUNlxJ4+fJJ5/EAw88sOjduNA88cQTeP7zn39H3yHtuFikDVcfacPzgbTj6iNteD6Qdlx95mnDlTR+Njc3AQCvxXfDoFjw3lwsLBp8GP81t8GdIO24GKQNVx9pw/OBtOPqI214PpB2XH2O0oYrafwkN6JBAUNyYp0pHJ5OwpUr7bggpA1XH2nD84G04+ojbXg+kHZcfY7QhiJ4IAiCIAiCIAjChUCMH0EQBEEQBEEQLgRi/AiCIAiCIAiCcCEQ40cQBEEQBEEQhAuBGD+CIAiCIAiCIFwIxPgRBEEQBEEQBOFCIMaPIAiCIAiCIAgXAjF+BEEQBEEQBEG4EIjxIwiCIAiCIAjChUCMH0EQBEEQBEEQLgRi/AiCIAiCIAiCcCEQ40cQBEEQBEEQhAuBGD+CIAiCIAiCIFwIxPgRBEEQBEEQBOFCIMaPIAiCIAiCIAgXAjF+BEEQBEEQBEG4EIjxIwiCIAiCIAjChUCMH0EQBEEQBEEQLgRm0TtwnvnlJ3/7SOs/ev8rT2lPLiaq6oF6FejyJlCW4F4Bt9mDLzVcpdBsaPiS4Ao60/3SDaPY8VC1R7HdQNUONG5AW7vAaAwe13C3b5/pPp0l5r57gc0N8HoP42vraNY17JqC13GF1Bzc2Yim3s9atsjtaOp51nadz8ptj2LXwexamJu3ge0d+O0d+OEQywy/9hsXvQsrA334txa9C4cyqy/Wm1NzosvUz466zhy/N6sv2qe+gvOKqnpQlzeBXgW+tA6uCvhSo9ks4EsFVxKa9bOdF9cNQ9UMXTOKbQtVu/C4PQDqGhiN4W9vw49HZ7pf552LfI8qxs8Jc9STab9tz9NJdpaQ1qCqgrpyGShLoFfCXQoGj680mg0DVxBcRWjWCb4AmjXAl4Av0A6YJw0DZhdQDaAbgu0RVKPh+hpq7KGbCmathBpZUG1h1teAcQ0/HMIPBqe0U2eL6vehNjfA1++B26hg1w2G9xZo+gS7BoyvAnzOfdHkAT0EXBkMvmJNo4fL0Os9qM0N4MkbS20ADZ7XW/QurAzri96BAzioL47uBlz/YvdFA6zEZMS8kNbQV64AVRkmAi+twfcM7GYJVyj4SqHeUPAFgvGzFsZDX4bHaaHq8CgGFMbGmmF7KhhDY4+ib4IRNLJQa2tQdQ2Ma7hbt8DOnd6OnWPkHjUgxs8JcScn1EHft+on2FmiL18G9SqgKvfOaFUKriA0G2FmKxs9BWCz8cPwax5ceNBJGUG1AjUEagisCMoCrga8IegaYK2gGoIeK3hN0HUBNXbQRoHGFmrcB90uwKPxSs96kdZQd18Fb66juasPu2HQ9BXGlwi2H262mk2G63twFY4/M2Y+A0f/DNh//cRB23XXmbVsnnUAwNwy8CYM9L4guEJB2RKmUNClgR5fBe3swm9vL+XgPrp6zu+IT5BlNX70lSugjfV9+2KzCbj+wX2xy1H722n06+l9OWz/aKygh2rfvlg0d0GVJWhnF+7WrRM9/mdJjn7oVcClTXBlwFURjJ4yRD+4iuCLMBnoSsAbwK6HZ18Cbs2DCz6xcZEZoCaMi3qgoGqAdZgYVJbi2MhQDYELCpODdQFT6BAhMbbQRoNH43MdIXHSyD3qJGL83CEnfUId9P2repKdJmQKqF4F6vfCxb0swD0Du1lFb49Cs6bgKoIrkL09vogX+AKw6x5sGFww9HoDUzgUhcN6VR97v6xXGDcGw0EJ32j4WoGNgqopXNRNNII0QTcEVTNYE1TD0I0JN8NjB6orqKIAjcZQdQ333G2wbU7wCJ4+eZb56iW4jQr1lQLNukKzRqgvBcPH9gC36UB9i7JvcWX9fMy4JqxX2BlWqB3BFwrKqngeEnSj4Q3BlAqqvgSqSuiyWMrQm/GVRe+BcCeofh90913g9d6+fbG55ME9D+pbXL26C6P8onf7RLFe4bln1+GM2bcvqmYNutCgqoQaj1fOA5THxcuXQnhbWcBf6sFVOkwIbmj4glBvUAz97k4CAnaN4UsGGwbWHVThYAqHjbXxHZ0P1ivsDCrYRsM3GlZrkKUwNjbBAGKFMCHYEFiHiUHVMLxRUHUBPXZQRoPqBqYswMMR/Gi8cuPiWXFW96irdn8q03iCIAiCIAiCIFwIxPNzB5y2RT3r91bNuj5VSEFd2gD1eiGOeb0H3zPwpYZdM2FGq1TBnV8QfAnYPsAF4ErA9hlcANzzQOmhjEfVa1AWFv3C4u7+LtZMA6OOFn60Vfew25QY6QLeE6zxsFqD2cDpMMtFnuANAAbYAMoQiAFdE3zNgDPgQkGNFcgxqDDAqIByLoTADUcAL/+MbPL6YHMj5xY06yrk+cRwt/AcZpqLvsVGf4wHLt1Cqeyid/9EqL3BM8M1AMBzfQOGge3HeScGmj4BUGACzEYFjZB6ppYw58CuLXoPhOOSPbDR67NfX0xen6Jvcaka4e7+4Nz1xZ1+hQbYty+aQbg10kC4fgFL1xf3Q62txVC3HrC5Di5NyO/ZKEJ+T0mwayqMgWshEiJ5ftgArgwhjygYKD2KfgNjHMoinA/rRX3kcdF6jYEtsNuUcJ5QawNrPBoGuFZwcTzUmkA+hNyxAZCW1wTyDGUIXCgYZqiRBjkPKgqosoB77vZKjIlnxSLuT4HV8QCJ8XMMzvqkmvXbq3KCnRaq6oH6PdDmRqvktl7CVxquVLDrOsQvFwTbb2OZXR857M33GFx6oOegy+DW3+yP0SsarBc1rlYDbBYjbOgxLpvDBz7HCjuuwpfpMnrGYmALeCaMGoOxNhh7AoyCNwzLGqoOF3o2gNLIBpE2ALEG1wSlCeQYqtCgQoM4XOypLFciFj0ZPrzeizdbGs1auNmyvRju1g9hNmU0fC5VI9xd7uJatYWeakMZHCtoOvrgNmu7eb7rpLb7wvBueA7B8ummy41Uzii3I8REBQW7vtw3Xc3GtDSXsArMnoSY3Re7kxB39YYz++JhHKWvdtc97nbzMPIFbo4vwTNhq9/DDrBvXzTrOm9nNjdyiMwy9cVZqCoKp/QqoCzh16ss9mPXwtjoSqCJRo/ttaFuroc21K3vAeOhS4d+r0avsOgVDe7uD7Cma2wWI1wvt+Y6/o4VbtSXsN300DMWHoSRKTBqDLwnOKXhGSCtgsHDBJWMHwa4ThOFCroI78kVYK2gmUEjAyoMtOcQArfCubEnxaLvUVfh/lSMnyOyyJNKaKFeBapKoChCEmcRYpldqeDLNobbFTHHx2DyuQg5PmwYpvAwhUNpHEpj0TcNetpi3YyxocPjWnEbao8ma4sHYddXAIB1Mw7LmFAaC2aAmWALDweAWYE1x/1J2bxRbQ4AmKALBpgAVvBlGIgVM6goAObgGah6S32hJ62DYVoW8KUJ7WLCQOZN8MClNkDhURiHUjv0iwbrZoxNPcJVs5O/T4HhjyHHd5LbzfNd0+s83Wxg6AqMCoPCODij4WKOmbft8XAFwjGqFag00GUJqipQXS+N+AEbMX5WEaqqti8W6sC+qA3nvrhhxjP74mEct88dhaP+xrN2AzumwtAVKLU7sC96Q7kvchkmm6iqgCU2fpK4ATrjoi/bcTGPhUWb95ofJo6JOhwLGA8dx8XK2Dwuruk6j4uXzQDrajzXuLjjqjwB1DclPBOYgXFhACbYRgUlBA77l77SG4DS64JAHMZFVSqAAVXoMB4yh/sBZihgqcfF02YZ7lFXwQAS4+cILMNJlbioHqCcyLmxDvQq+M0efBXqE9iujPVaULBxMdQtXeRdv5PI2XfQhUOvX6MqLPpFg7t7A6yZGptmjOvVFjb0CFf0APcXz+GKGqBHFnrqYu9AuOEuYdv1sa7G8CDc1n30dQPPhB1dojQOngmN0bBGwzFAtYLzCtwQWIeLvNKAMgAxQRuCLxjEBmqsoKK7P8x0FVDeg4YGbmf+m5KzgrSG2tgANtbg1yu4jQJNP4QgToe7oe+g+xabvTHu6g9wtRrgvnIL14rbuN88h7v17qL/nGPzjFvHV8pNAOFG4LleiBvbHWkEc0bB9dMNHKHpKwAGIEBvrAUjl3lpVI14YzmMMGF+SOtwvYx9Mam77dcXe/0698V7qp3cF7+++tKi/5Q74n+Mnw8XjaWn+60W36y+2KwHTyxgoHer4GEAQEuqwqg3NoLoT1mC19fAPQNfFbDrrdJpCnXzBWWPTxsJwfA9Dy49qGCU/QZF4VAVFpd6I2wUNdZMjXuqHVw2Q2zoEa6b29jUQ1zXW3vGRADY5QK3/Bq2XR8Na6zpGrd1H5YVetpiR5dwrDDWHkMGvNF5MpA1gZv9IiPi8iIYTKo2UIUOk4PGgMoC/Jy7cAIIy3R/Ciy/AXRk4+dDH/oQfuZnfgYf//jH8eUvfxnvfe978Wf/7J/NnzMz3vrWt+Lnfu7ncOvWLfzJP/kn8bM/+7N4yUtektd59tln8WM/9mP4z//5P0Mphe///u/HP/2n/xQbGxsn8kedBst2YiWOc4I9x1/BF/B72MJzqDHCN+Bh3EdflT9nZvwBPoU/wuOwqHEF9+Bl+Cas0WZe59lnn8WP/MiPnGkbktbB8KminHVZgMuQG+PLcOH0BeVnVwBeT3p90uwmFwxtQqhbYeLslrZYMzXWTY2+rrGmamyqUXimETZVjU1y6E3pfY6YsetHgAYcFNZUDacVPBPWTA3PBM+Ewrgsz+oLDeY44wgA3OYAtbNfHJcT4AmABpcG8MCzu1/EH259GFv1UxjzYOnaMMyUlqF9Sh1nm2nPrCMMg4yHNh4902DNNFjX6biPsa7GuF+vZr7BgBkjHmNTjTDQFXZ1hZ5pUBcaA+PBRmHwhc/h5vt+DeMnvwS7s4WX/Om/gGtXXwFfqNDWZYHPNb+FJ/iTE+1Yoq23c6bXUyMx9dOMPvM4tv7rh1D/4R/B3drGvf/v/wfWXvWK/Pky9EV0+6KhA/tiVdjcFzf0OPfFK6rB2onVADhbBsxYU21fXDMNhnP2RS517oufx6fwJf795eiLkVTbDhPtHMfF6fEwP6P1cqUoiIJBBUMZB2PacXHNNHlc3NDjPC5u6iE2aYSrqtkzJgLAVzzDsQI0sOlG8KzgtMK6CSqqngkDY8OkYOFgGfDQYXxmwKPjBeLgGQfCuKiLkCPrSw1igNmAywLPDb6Ix7cew5a9iTEPl25cvGgsswF0ZLW33d1dvPKVr8Q73/nOmZ//9E//NN7+9rfjXe96Fx577DGsr6/j0UcfxWjUuiHf+MY34pOf/CTe//7345d+6ZfwoQ99CG9605uO/1cIR8LBYgOX8TJ808zPv4DP4gl8Di/DH8er8W1Q0PgEPgzH7YzXD//wD0sbLhDnG2wW9+Lr1l8783Npw9XA1zWqa/fj+nf8uZmfP37rY/jC4H/sbUe07SjX08XC4xrFg8/D1R/43pmfS19cDebpi1+0n5G+uMSkcfHl1Z+Y+bn0RSFxZM/P6173Orzuda+b+Rkz45/8k3+Cv/W3/ha+93vDQPCv/tW/wrVr1/CLv/iLeP3rX49Pf/rTeN/73odf//Vfxzd/8zcDAN7xjnfgu7/7u/EP/sE/wP33338Hf87psKxen+NyDz0P9+B54c2Ut5qZ8UV8Di/Cy3Afhbb4Y/wafAj/GV/Bk7gH1wEAH/jAB868DalqvT4oC6BIHoUQv+4KygVM21mtyZhmb9IMV4hp1tqjZ9qZzp5usK5Trs8Im3qETTXEJTXGJjlsKo27VCt51bDFDte4pMaAB5xS2NDB0PdM2NI9eA5eoKGx8RgTSHsg5vwoAJ5TTHN83QDgEOPsihiCwT7UpWCNe69+Le4rHwBGY2D7/wApldty4W1IKs9E+tLkNnLJ+2aCRy7l+qiY79Mv4vGPOQaX1BCbqsY9eg1qxVT5PTye80Psco1NPcK272HdjNEvGoydgSo8nPXof8PLUH7VK2BGAN4bjosrQky7KzS+cOvj+JpLr8F92/cD7HM7Po0vAwA++9nPnun1VJfLF/KzaDZe/WJsvPrFAICvvB1QMVEcWKK+WJjcF12a/d+nL1bG5r64prt9UeEu1V/ZvnhJDTHQFbZ9Dz3dzN0XfaGA2Be/uv9NuG/8/KXoi0DH61OVeVzkqjMuxlwflzx9ptPundxXb8KYSDqMi2XhZo6Lm3rUjos0CuOi0tigEgW1t5MNW4x4CB/HxU09govnzZbuwTHBs0LPNGAmjLUBFwTmcLw9h5l5b4Lnh7iT/8MhnB2cVFJj+Fuhce9dL8W9618N3N4GRr+6XOPiKbLM96jL6v050Zyfxx9/HDdu3MAjjzySl12+fBkPPfQQPvKRj+D1r389PvKRj+DKlSv5xAKARx55BEopPPbYY/i+7/u+Pd87Ho8xHo/z+62trZPc7QM5yZPq5f/rj+xZ9uk3/ewdf+9J5v8MsYsaI1zFtbzMUIFLfBW38Uy+QFy+fPlIbQjcWTuqfh/q8qVg9CQVm8qEmOaOpPWePJ8y5vkUyCo2qnBQhUdZWmjl4Txh2IT4Y0WMLd1HTzfYsn30dY1KWfxBsY11FVz+d+ttrFFw3T9h78Utt4Yv11ewZXvYdRVuDjex05QY2QK3Bz00jYFrFLxTUNqDFFD0LJyNEp8NgZSC4xDrrIqY/2Pai7+vAW2Ci18bBTY6xKHrIIZA/T4U9+GHw4W3ob58KVeQd+tFzjGwUVI35BcwXM9D9xx6vQaX+iPcU+3mHIPr5hZ6ZLHtS3zR7uKy0nt+Z5l52nk85dbgoXDd3GqXVyF0YqvXw4gJDoBbU7ksvaso5xs0egdjt4urd78UWhH89jaMAy7xVWzhOQDAxz72sTO9nq6vjQ9fCcDOTu/wlVaYjY39E6qrqsnHaZF9kbSG2twEbazDbfTavrhOB/bFK71h7ov3F7dw3dxCQQ6fqfu4Tw9wj14t4+e2d3jGV+iRzX3xqepSTtQ/rC8WQ4Nm+CzGbhf3XHoxdL22FH0RQGjffg9YXwP3Qhi4WytiMdOQ+5oMn27uq+0HZbeseFowdM9CGw/2hK3bfdxywXD4Uv8uFIVFWTjc1R/E/J8G13u3ccmM8LzyFq7oAa7oXWzSCAMuMeI1fMV+FXZ9hYEv8VyzjqErMfQFnqvXUDuNkTNh3PUhJFwphtIeNQOsNXwTwr3ZhHGxLQvR5v9kVVSjQB5Q46iK6hzwdBgXMQjHatHj4mlw0kbP9D3qSdyfAstpAJ2o8XPjxg0AwLVr1yaWX7t2LX9248YN3HfffZM7YQyuXr2a15nmbW97G37yJ3/yJHd1Lu70xJpl7By2zp2cbCdxgtUIg3qJamJ5iV7+DADuvffeic8Pa0Pg+O2oql6IWy9MSOjMKjYqVqdOUtYUjZyOx0enmS7OscRwwegYOgIpYKg9SDGIGE+rdRABihhae2jloYlRGgtD4X2pHQx5jJzBrWEfo6bAYFDBNQrcKKjdULVa2SDYBkL8zvCaAVjF4TMFgILKTjJ0svJbnvGiqHjD0CVlWVZVGqQkIioMCCUwHC68DaksouKQzjOPfo9HjoGCYUqLqrCotEVfN1hTNRpv8D9GD6AgBw2PnmqyslPNSQraw0Hted7vMwD7rp84aLvuOrOWdd+XZFFRg4YNGtYY+AoNa3hW4W80DarCwjoVFKYKBjUp0TeqMZVAzUHooeptgsoxqCzBw+FEO968efNMr6e9Yr78q61b5bG+f1Xo3bW/yEipfT5Oi+yLVJa5L3LZ6YvFwX1xzTS5L95sLuNZu4GCLNZUjT8kC02MW671fh+1vx21787Tr6f3JdGjBh4ExwojLnI/XFP13H3RFYRR7Itl/3IQPliCvqiqXlQ8NUHxtDQdxVNq29pMKZ52oiFYM2AJaBT8jg7DCYdxSKXxh0uMGRgD2NJ3hXG0ZGDNQhcea2uhPERVWNxVDWFZwXqFxmlYDqI+tdNB3IcB63RWQHVegT3BOwX2AHsCu7AeURi3fdyfVg01jY8EV0ZVVASlTLAOQ2oZrj9UmKyKuuhxcdlYxfvTk2Ql1N7e8pa34M1vfnN+v7W1hQceeOBUf/O4hs88J9Q825+Uxb1MHKcdVVm2stbR8EkXeV+oNnnX0MRFvRv2Bp1UBgByQWKTFAWZVGJ4R9H4AZzyIAKIGKTCdooYSpXQyoOIoRXDOoXGaoyHBbzVwK6GagjUEFTdJn9280C5FfUK/1E0gNJKnYt9CHuLg09O9EyeIJWTPVVSUNAGZEqosgyj1ClxWBvuJ2+dw2ym5K2N9lne2iiHHVdhzAXG3kARoyCHQgWFPUUeO265PQprqsY9ZhtjLlBn46dEwxojX8Aoh55uUGqHsfawhrPULhBvuJLsdVQ+8oUBSj4xqd07uZ6uFfVc65FdzcT4eTnoOFSmmfs43QmH9sWuvHV5sLx1ty/2dOiLA1+i8A4FOVSKoDhMQChm7LgeBn65DVwFRk81cAghViMf6q41rDH2Zu6+GI5ZPJ+LAnC08L6YxkUURWzjVtY6t3PZKfUwNTay4RBFMGgLvKaQsvQ+QZ3XaIAwdgFuTHAFY7vW2DUOpnTYqnow2ocxzFM2eJwPv8NMwcDhzuv4ABPgAXYE8gS4uC86CiDMmBxME4NhnGxD3KiKt7bagHoVFPuFjounwSLuUbvbHucedZkMoBM1fq5fD27Dmzdv4nnPe15efvPmTXzjN35jXuepp56a2M5ai2effTZvP01VVaiqauZny8SdGj6zvuuoJ9idnlxJuabGGBX6eXmNETZxJb//yle+MrHdYW0IHK8dqarC7GUyfArd3lSXCm76Aq+RvT+cZraAeNEMF1cCwAjGDogQnSpgAohVXM4gH9bxAEgxLCkQQi6Pswq+UfAjA2oIeqBy0dKw42h/F1PL8ntqV+saRhT3PeX/xBkvYsCVKQZaxYt97MKFAYpQi6Jym4BdTBsmZSkUGlx0vT7dQZhnKkvV3sAToWEGNKDYwxPlWh4FuXwDAwAeCgoHq4/NWmee7Y7zXSlPTJNHE/fTs4KLj6T6VyqbVd8aG+qM+HjDlc5fbwlm4xIAYIwB1qp1UB1uqGuMsI7w2bVr1870ero+r/Hjz7fxc9BxqLTFelGjxmKvpxMe2GLKA3tAXyyVzeeqYxW8rvH8dVBQ5LChRxj5Ajsu/P5Z9rN5tlMUDJ9koDVeB68PFBqvMPYGljUM+cP7YtnpizREr7qy0L5IpmjHxSpFQrS5ld7QzHExj4kmTk4kQ6E7Rk3lANPU++46eleBSg7hgoWCbzScU9AmREiEWj4xXI1bgwfxNTh6ejiqmcYxmnz3ffzJzpjIHP4eJJXCODGoy2TIhfwfAEBhQLoAvF/ouHjSLGpyfvq7VnmS/kSDd1/0ohfh+vXr+OAHP5iXbW1t4bHHHsPDDz8MAHj44Ydx69YtfPzjH8/r/Mqv/Aq893jooYdOcneEY9DHOkr08Czai7jlBlt4Fpdxd152+/ZtacMlRdrwfFCtXUVZbuLZW3+Ql6V2vIS7AACvec1r5Hq6xEhfPB/kvrj9eF4mfXG1kL4odDmy52dnZwef+9zn8vvHH38cv/Vbv4WrV6/iwQcfxF/+y38Zf/fv/l285CUvwYte9CL8xE/8BO6///5cC+jlL385vuu7vgs//MM/jHe9611omgY/+qM/ite//vVLo6RxHKv6JC3q6e89ae+PZYsh2pj1IXaxzbdQoESP1vAgvxiP49NY4w30sY7P45Oo0Me9uB+xKg0eeeSRU29DMkVw7c+qX2C6Mc3oPHfimZNpH1w6sR4AAB88LkxtxBlT8AfFf8ErQ5zzdeBjbg0x7NgAVgGWoIcqeH6GUzPdacass3jmLNre1cLmGnmeM9X/IQ4zeeQBZ8fYHT8FXQdVqQHvYAvPoSgcev0eHtx9KR73nzrTNiStc26WLw18pdu6EklhyEwqS6WQNwXG2Bk4RSjIYxwPYKMcCnZQxGigYVmj8QqeT3Te5kS4r9zCugqxFQ4qPqKXigkeKs6gK/SLBiNbQDdD1F9+GhiFv6fefha7z/4RKqyhKi/jq776tXj8c7+Kzesb6KPC79FvoOJ+Vmt86UtfeqbX0zUzZ+HAc14OqHsc7KDBzpdut+9vPof6D76MEQ8Wdj1N4ae5L5Zq7r6YztHk7QGFQs4OwQuk4bGuxhjoEgNfYtcuV1SGIo9CeTQ+5P+kvjd24Xan4ej58RoedGhf5HEfG+UVfNVXvxaf//1fxfrzriy0L6pelcfFfWvddZTduuMi6zB+0HFEG6e9Qg5x3FNRMY7gPMGZoBiXB7yOpwcI3h6g9QIhjsl5bJ7y+qSwCFZxLEQbDZHq4RHH3N+6xmj8NIom9M8wLt5GUagwLm4vx73NnbBs96fA0SKUliX07cjGz2/8xm/gW7/1W/P7FOf4gz/4g3j3u9+Nv/bX/hp2d3fxpje9Cbdu3cJrX/tavO9970Ov18bq/9t/+2/xoz/6o/j2b//2XETq7W9/+wn8OXfOUU+s0zqpZv3GSbkYt/AsfhMfyu9/H78DAHgeXoBX4NV4AV4KB4dP4+OwaHAF9+Ab8Vpo0rAcLhA/93M/h7/xN/7GqbVhuJEuQkxzSuosOjHN3cG8e6Ev2gtke/EM12EGYgwbALQGEGJyZaifRvFzjgYRohEUL8CegKEGGgVlg9FDM3LAqfPbxyKlKZmu8RMLvTFh56kv4VMffVde/bM3PwAAuH/tZfiGtf8JL3LfADeoz7QNZxY2LabysWYUNjXkUXsNRSGnwJMPuQVgNKxgSUORD5/HmP3amyyCkEjqTdPL/JRpedB26bN5vzstW4/F/xQ8/JRDffr3Q+ibQ880eOYPb+LLb/03+bOnPvCf8BSAu17+amz+8dfj/pd+K3w9xu9+8b/BuhGu0L34Rn4tNFr1u7O8nm4U+6ucdaFzbvx0j8OXP/dlfOBHfjm//513fATAYq+nMwubztEXS+VyWGli+vwFAAWPNRUKQT/tN47cz/Zb57Dt5umfpbLhOq+KiYmSsTfwIDRex6R8Dc+UQ9/264tXX/JqvOQ1oS9iZ4jfvbG4vpilrYuibd8oLBOEUlIo+N5xEYwwVt1hROrU6QE9JHBDoILhmIAihIVTksVP4W7xdVb9ASaNnmTw+LivnWXUMYaYujk/yIO4Lgjbu0/idw8YF184+Do4t9h7m7PmrO5RV80AOrLx8y3f8i1g3v+OjojwUz/1U/ipn/qpfde5evUqfv7nf/6oP710nMVJdRpcpfvwCP6v+35ORPgavAJfg1fsu84yt+HkzFYUFmBqxQcIYA8gqZqquIzCutk4irk/HBM8QQhJmI2CcmEgUXNOhE+wXzLplMF0kKfo8n0vxp/+7v8FuvYw2zVo7KDGDbC1A4zr0Ia0Gm04cgZGhYHHsIIhDxtzDRQxrHLxtYdnhYYV6mjheqZsFKl4wLo3b+mzLt31591uv3WAcOO1roPHx6Ob36Nyzk8TPVYNK1hWqH24Ydr8hhfghf/6/wMeGtBIoXpWQQ8BMwTwDIOI8KIXfydetvknQbsj+CdvwA+HsNyeeMvSjl3I3eEd1grxvFc9D3/hY39hz/IvPRS866t0Pa2jYdCwgornrWJG4w08WfiopuapNfDHLhgV+/Wz7ntgb59N68xadtTtAMAxwZLG2IfrRVB7o+zpsdn4CY9ZfdE8Z2CGgB4CvWcADENf/JoXPIKvvftPrVRfTGT1UaA1gDpNNtF8B3TfWeMSRSEELjhO/niwVpNjWvLgcPs+iRV0c3yS54cYQVzCo324g8fFK/d+Df70d/8vKHYs1NiGcXF7EOrhrdi4eBKc5T3qquUArYTa2zKyCMPnKB6gZbCsj00szJdD3oqOwlup4Aq1V8YzuvQnvqY7kwREA6Zj4DBlz0+ybyY+B7L3JwggENAQ9IigGprf8DkoibQzo9X1FnVnuqDC30dI6jbBzR+kPQFVaigPMHMr8ckeGI6ilXcGHFTYtJxUGeoWNrWsANcaFdYrGOWzEQRgwhBS5GG9Rh1vuBaNAuOucoC+bnLI27wQMUrtcpFFtmGmVplW8U3Fgqe+NNC1Ccd4ND6+R/EOuWTm8/wsav/OirmPwyI4qLBpOTvkrTAueLiPQAp96+sGz9VrS9MfFThcM72JIW+UJyNsnJioXTKEWknmib5oQhiXMuGYORv6oi8UaFF9MbVrVD7tRkO4fL1NstbtuAgVJwNjCHcqvwBMGkPd1tvzJ83RtKoBMFCggsFOBXtGte4a6o5zXU/QlOcnf86t0ZM/7+wc6zb8DTEkPHymWtlrjyAGlMLWqwp8luPiCTNvZNKiJuaPYgAt+h51+QLnF8iyn1hH/f1lrvp7EKpXxRCqCqiKaPgY+ELDm1DDwGfXfifHR8dH9uS0D+Rnbq+78X1m6n07MHB4YwnkOobPvIPedDrQ9EAyYyDaM0CpblgfZSnkoN6kYx0PA1RFULipKqheBdDZdHHV26emiJmOPW/DbKrCQlGor6SIs4EDtGEsinx+nX+LGEaFWkvpZuegx7zrzbPd9LI1U+NSESqeF+Sg4KFjiF4I1Qvvg2R3qBllyIeH8jDaQRsPMj7Gzbfnc1eulksNLg2oLEK7LoigaHf447yzzMcg9cUc8hbrvuS+WMzui0b5fG7q2CcLcuH8VTaHnmryUAjn9BU9wKVihDVTn0ifu5NH6FOuvYbEyZPp64qaut6k649WPvdFv09f9KVaTF8klcdFFCaPiyG8WLftO1X2AaozltDk2LPHC7Tvb8+5jxwMINVQUJOzkz846VnaOyM4PU5P7LPiyb8jjvXdfCYXlfnCea7yuMjdcbEsz3RcPElW5X5u0ffH87J6Z4AgCIIgCIIgCMIxkLC3yKp4fY7Kol2Lx4GqVNS0aBXeSh0TOrszmNTW9onqbtnrk14rnpj9mqyn084uTcxuzZihIhtms8KsFiZjmeeZGUvhdYd8HoXpZi5HdPNTTBpN1a19QSBWUECYkey4+JVn+PHph+hMFFTs1JvIXrlOQUVtGIVxIcykMysLIM/KTpOWaWJoCknZQWdi//mblCPU3f6oHLQvALBuQm2fy3oQKs6TgudQad6DYu2RONPMHoYcnKIY2udRKIfCODijwcn7YwkqzWZaQJlOuE30iNJYA8dRbLpDLps5CzsuPgLqVJnvOGyc+n7MYt/w026Ryxl90ajojVQOhkJx0+Qh0Qgqb6FSDod+CI9NNcKGHmFgiiwikEjez5NSZjzsu1T0WhUqdAw3I3+ou28pT0gRg5nhiXNfTMVOvW29Zc4iXNfqs++LKv1WKvZdplDwVlTGlZOeEK9nRBkAe8fBKWZucxjdkG0b8tyoCLF0bJCjJzp6QnEgmwy2Y1AOtCAE0SFSbY7uxO9lwQPEmj8hJBwM6JJiqJ2CKsNtLnkOxw84s3FxESzDPeq84W+LvD8V42dFWbXkstOGHKBUTO9RbWIkp0RK1eb5ZGMnXbAnwuDa9yHdJ7yh6MonS3uVrA4ygI4ibjD9OieDTuUAnROcJ4ydgeZw0+KZYDiE1pgcOsY590cR5+TkOsXzx8d+NznH4SjfdXeFUMQ0jswubteVus45BwgCCNbrnGztllCy+yQ472pv5w3PIf/FUjg/i2hoeCZ4oj1S144JmuJ5Hs9hz4SBLQ75paP1r6N8l1HhmmE7fSoZYym/J10vkuCB56D+5ryC8+erLxJ3RA7ihGAW/Ul/Ku01dibMkelJu1nNNh3B5sPvejtjRnF63POTy7viB3nsi++VQx4LyYX3x5LsPscsg+GTWPZ7VDF+sLpen2U/uY6DKstOZfKilbfuJnR2qlenHB/WYaYLXe8PpddTHp7pXB9MxRnPgBxBHVfdrfMb+zkgJmbEppdTa191kzzDTFeYjUTcVhUaoe4Bh+PIwePg6/0r0p8ElCXJ9UQ7tfkFyNXkdeFigvWkp+cgD9D0OoY87JTX57jenVnM811GeWzoMS6bAa7oAYDgmTrI86NiXkLy/Fgfks2bwoVZ5oJBTUeO2FIQjCgVVKOBqgTVdWjbOZ0wJ8nl+HdedOY7Dgvy/GQZ5E5fLGjC8zPdF3X2+qS8GZ9zfPbz/ABAj2pcNgOMvUFtDAa2nGsfT7qvdidKZqnD7flNnrymcHykvugMt33RtH3RRQGSs+6L1KvanMqqI3SQyj6UFEs/dMbFbu7orMiGGYZPl/1yU+faXxcNIEren2TdoDOY7R0T03iHtFoeq2lyLIwTnVEktB0X0+RTkQb/KAjEHMbFqozeIQZWyPMzzz3qst2fzsuivD9i/MzJsp5Y59EAOhK+vVgqhYliakkthjW19X2SB4g6HiBCVDqifeU+cyLnfjNNs7w/c3h99lN4o1jrIP0N1Hl9HjxA3hNqaDADRBpatbV9iBhWqZmeHxXrAdVOTyg55e/tSN4eaX+Osd0aNdmb41gBNJ/nJ8+oM8GxwgGVA1YW8fysDsyA63hFpoucHub5CXLSCo4JI7f/LcWsPjZPv5t3u2T4lLr7t1Be33b7nVe5X9bWgDl4h7w/wh3+shI9I7HKQ54EzLXsqH3dNYxiNHUey9LhzZfXQ8a4aVQTvsxTZ6NpLxHHH5gx9iUvEKVxL3l9kuR15/m8F1VeVZb5/lSMnzlYVsNnXlYl7yfFNaOMcc2VmYhb91NeH2/QenrU5EzXHoWbrmdnOt45XuX3mwUjF8PdDnOxHxD+NjHD1QklyMs73p08I4Z2/TxZRmiV32KBNx1zfjwDvjTBHxJlrwkAvAdO2fMT8rOKtpq86VSTz1XlGejI6galJUCrdoY5K791FJlmYZSH9WrimHumXC8oMe/N1TzbdZdlpTc9xIYeYVOHqV+NEKDe9fwU5OBIhZwf5dBwlPP2QQnO6HBMXKH25hqYpGKk4BsFLovQrmW5EM/Ppl6d2dLTZKmPQ+6Lndw7gz25d6kvGu1RaJcVCFvvpGsV32Z4fhSF9xt6hJEv0BiFp2gzS15P96Hpejxp2WHMu11QgfR7tpu1jSeGigaPVj6EvHGYBCuMQ10E6XkuGFwQvO2qvp1tX1SxeHQeF8t2XORCtR6f6KHi6RzXaU/PjM/uKM9nH8gh5P4Qt2UopsPogNZC63qF0ss49kGhLVQeyggFAy0NAdzm/4RoiFRDSMEXCuAwLlJZAp7DPGhZnnpExEmw6l6feQygRdyjnq8gV0EQBEEQBEEQhH248J6fVdFOv1OW3ftDpmgLm1YluOio2aTY9akaBnmmK+b9eANkpbdc74dz/HMOe8OMWa80KzjtGbJ0fOGq7oY8ObtG6b+0nDueHe7McnEMWYhhyuTbZ1+ElV3ZTo/5UoFY5/hmACH/ZzgC2ztIWDqMmO/Dpe7k+0y2FUyoK1IYh0I5aMUojQ0zycrnhOU0Aw1gKuyt9RBNFy48awx5bJpRzvfpUYOGg78HFLw9iYI0HClAAWNv4BWhVBZeh5C4kXXggtBYPZlrUCDMOKdcg0IBRcw1qBZT6+eK3p1rPX/OR5Z5j8NCSLl3xVTuXX5w7otVYVFqh1LHZ2VRKRvqUu3z0NHjoxHqAV3RAzSsMfYGG0W9kP6owNlTla4PQEfwIIa9pde1D2GzBRMap+F1CIWroUNEtPFgo0K9nyIofHobPArOnl1fTOMiyrIdFwudx0VXUBgXC0zkdbX1cTrjTkoN3TP2pR9rf3e/8O/jwgTApEENbeRDDHujXAwVnXUI7DkKFwGkKI9/Km1LcSxU7XbE1KrAccj5Ccs1uDAh5wcI9xueT3dcvEMuyv3pIjjnQ9TFYJnjKk+NaCDke0yOccydcZc8wNFvThSLpKl4JZ/K92GausJT+x1kQ/LmXLk206Fvc+T95M+mX+e458mcn1VXfvOeULOGgQesARFDs4fzPoe9mU4oXDcMDmjDV1L8fu31RHhcyh2a+M05l+3Z17hOd10Pwoapc15Eoqv45mOORMr7afd3Rt6PV3DnIdegg6gwrQbOx/Mvn4+d8zOdzyCoeD4HoyduCwUdNZ59x7AY2OLQvnbcPpqMmullhjwMq4lrQ3rOD7SvrVfguKxxGswE5wl1fL2ycFBBy+IAaTxK82MpDyhetmZFoe3J9+HJz460LzZsxCW3x/WgvB+gzRdMuT7xM+XaZTkHqJv7s+Lj4nFZ5pC3xDLeo4rxcwircGIBy3lynSatDGZcwMGu4bgcGq1xw8EIok6cczaGCABR63UBJmbCyMbEzaNeVGdd4Luf8dRynrx4HyR4kP72VaWxGkqFwZB1iAlXTsFpnyuuOx/+wJwP1DF8Jj0/4UbmLFHE8Ka9C3BRzABoDaAu4UaS9gogxH13TPDnTWZXjJ+VwHsFx36PkeBA2YgP5297AXOsQg5bFD7oinx4UBYjOUumJ0y6eUfdSQYgGEtJ4IGZ0ERDyHmCdXqlRQ9y8j918kZn5P1QWqfz0cSQlSb/piIWDsprnUk0gHxD4W5z1qRfep/HvklDaFrsJ016Tggf+M66F4hVuT9dRi608XPRXIrLHPqmelVIRoyFMtsibrqT0JmS6CknduYEz5DT2CZ66u7rrqEzmfQJAKymLZX45DtvjsrUiMKdwWQivKC7fsfbwx03PqJRl4URGCAV6+pFF38yAlWpQ30HICTkRllP1avghx7sTumu1HTEKYqO2EGnqCIKD1V4lIWFVsHIKbXLBk+RwlY6icshnKWj9tYRHJiu3XFWGOWxYca4VITippf0CAoeOlnLMFDMWfQghcA5UqjCVChKbbMR1DNhmfM0mWjdSVQPYh8KvjTQtQGKxVy6r6j5pK75nI8s8x6HhZD6YrlPXzSTfbFnLArl2vC3GPo2Hf6WxA4KsijI5rC3S3qEhg1GbLBhxrmG1VnRvWbMMnzS627Ym+0YPzqGwDmvoBXHfmjgbOiL3gCqIz6izqgvktadcbHI42IIA49Fv4tU9Du0cxvujU6pB+wZ9yYEEfIPTj7vuaTe6SW2K3zAE08T9X6YORs8edxjamv3xSgPSs+dSVCKYX4qyl634eFhXKRU8DQKAilmuJ3lDXu7SJz1/ek5H6LujFWzqi+S9yfNLFN0heeraFR+Cx/GB+/1+gAd1z/RXq9P/O5Q2PSIO5cMnP28PwfMfu0JeZue8fLtslWVE26shlcMpXy4+SAf2id5gYhhfVfeeq/nB8CEB+Uguevu+27tj8Pkc2fVB/FM2CjqKFNN+Uar/bydXU6EWfRJ+d18AwbKBRZXOtxmCvH8rAb53FPxfFST56hjBUUeqXRpUjD0rKA7F6DsIYrbDWwx0c8O6m/z9s/9tsu1v6Lxk/YnrZv6cTLIUl/k1A9dG/bn3OqGoE6MB7O8PWpiPm3S69NdEJ8nQt+Afce1Q+EwjgJAjplM34epcLWOsTMxBqZIj/R+P8+Pw4Xz/qwKy3Z/er5iLQRBEARBEARBEPbhwnp+Dgt5WzWvzypDWgOxejWKIqiGVUGtKNdVyCpvIQwhV7Du1vnpuPvDF2NyBgzI00I8NSs2c78sQTXH9LAcNEs2NcvGNDVb1V2OznP8jNCGxSXFuxTA7U2o+QNmqEIjFNlhoCxA3gPjMXAK6jZcGiCq8uVQjCmFKTIe2gQvT/K+EYVaP/vV7Zh+n2Z2VZLDO2SSduZ3zFh22HbdZWumwboZY0OPsalG6FENB4WGNRT5MFse3XMhKZzhyYewPt+pRq9CroXRDoUPx2y0n8pUDP3kUodjHZX8zppNNWddjNWcPJ+buY/DAkh90XXDojp9MSm9aRPOSaPdhNpi9xGU3WJNHwqeFYVQCyjRoxqbaoSRLrBuxqi9yZ6fRLffdT2qR+mfs7abDoXNHqFO+Fv3e2Z9FupdcygfQwwigjYePvfF0AfJtiHXZ9EXg8pbHBfLzrhYtuNiVnozRwxx6zAr6mHPZ3nB8f4W1SRv3N6aPzk3acLzlGK8O0/dcTOO/WkYyPsav9sXFFXdCLoMtdeIFVSlYyh4PK7sQVqfXjj4HXAe71EP8/6cZeibeH4EQRAEQRAEQbgQXFjPz3ll2eIqT5yU7+IQqj6nnB8fSviwCp+zRjtLlJIhCeAodU2dmTAi7PEIEeP4Xp/Ovib2qL2lv2VWvg+mXndjnTt5P8qtXt6PdwRLCt7Hyt9McCrOMnvKXqCU+5Oek8eHOjO+AHJewrTgQcqfoc6B32/ZUdeZxiHkDIT96Uhfd/bLsZqd98MqSw27c6T4Jjk/KwAjn3fTuT7dnJ+CHHw8fxW5PTk/WeK9M5c6duHWYt7+dCf9MF0LCpVkt2niO5KXJym8pdfMHeU3H3KArIt91K2W61K5IASQvT8qCgMArZeko/AWF+9Ve5v6s/eovR0T8mE8DXlWnL87/DDNHvfiuN5dPivnp5sDm3Nhz3nezyp6fZYNMX4uGMuo+EZlt5BbEZWKNHyhQ4E+Ewr1uVjEjXXr5uf0PFHYtOPqn3oPmlR+A/Zx+x9V5OAQ9qi9ddz56CzK0txdtTegVb3xkzZdkDYlgDkUdyvjyOcBXxqEOq/cFuJjBuzoZP84AL7UoCIUWXXdIrRRXYo1Qxcexnho7WG0A1G4YSEKkta6I1WbjB4AORwHwGTYG7CnfsdZsGbGuGRG2NQjbOohNDxAgIeCiw2kyeebroIcnCLAo1V7U611UGmbk69V4YPKVDNDZcoQfKGCYpFZlNrbfFYN68PXWWXmPQ6LIPfFpI5ZYKbSm9YepXaoYoHTnm5QKjep9JYV32wIgSOPglwMg+OocuixqYcYcYFLZoSxN+ibs1HQmhA8iMaYnzJ6sjGnaGIyItfYIg7KoSAopWCdDtdjptAXTZC2VzHE2Nkz6IukQL0qj4tB5c3EwqaqHRcNohJqZxxMxk9n/JsZCoe9r7vPp3U5TUqQ3Pkviyt0jBwGAE97C3x3JjkpGXW+HQ9b4yeMi95QPiiq1GGiFACVZVZ949F4qQymi6ZE3OWs7lHPz1SjIAiCIAiCIAjCAYjxIwiCIAiCIAjChUDC3maw6vGU5znvp1vMbEIdppMfk1VfqH1P2d0f8kqANgwgxEW3OT/ECEpv7mRd//uljXSXH1bnJ9U6SLk/K1XVOtbRAJAL1bEOz5oJOua8KGJo5XN8v4s1f1IonGfKSmnhczUR4pLWSRxUQ2R6vcO2y6E1M+aN3IyTxUHBpX3r1EgJoTgq1/oJdUeOdjiXmhXLRbuIpFDLNmxUwbMPz1CIZZThQFDo5rIRCgrPauqUT/2icToXFQXm61v7LTtoO6ANhU1hs9arvF0376db28flmj+U839yzo9X4Tp1RmG0d0rOAwXC2KEB9sjj3L51fjqh2BN/ahxX99T5uUOUC5cF1UyFgac3++S7zsqNzbk+3To/bnJcPK/I/enJcCGNn4scTwksX95PK+kZKlhzqYPEaKngCwVXhpwHNoDXs3N+QJ2cH2CfKte8VwYbnVhoIFxQTzjfJ3HkvB9u18nVrmOyD3Nw27b3F2ElXShQvIHhUoPDqB7yqrwPstfDzih4Un9bocGlauV1Ozk/vmBwydBFlNTVHkp5aMU5n4eiwTOd87NH1hZ7pXC7ldsPKnA6z7LDijIa5bGha1zSQ2zoEQqk3I+QAxHuOgDFHgoKJVk0FBJgPCmYmP9TaQsPQk8TBsrB6XADliR2g5xtyDGgJkq8l4AvFVQTJXYXwGUqDl8J5z/nZ97jsAhSX/Rlpy8WbV/MMtfGodAOhXLoaQujQv6PUQ6GwvKQ7xMeSeZaE0On5yx/7bCphrikhxhrg37RZCMksV/R4OMWG54lcw0A0JO5gAnrJwVIkjGUhBBcFFTwPlx/nGJ4S+GYNZ3jaE+/L6qyABVFOy4WOo+LXFAeF7vX2pTj2h3z8tiHqdczcn+6z6dt93GcnAwpPzwxgccdkQOOYghdoR/mTq4PtflARJM5P4qiIE7JAMK4yEaBS4CZQVUBsAdVFZRnYNwsRd7PRb8/PSsk7E0QBEEQBEEQhAvBhfT8HMSquxTPM9Rx5VP0jviOBwVAlHzpbIPujBe1nqGOt2daPpUsQZ2C92c/uev82fTrqVC3mSFwSe5zBfAuTUUC5AleMbwOjaCVh+7MFKewt2mZ664XKGG9CuEsZ1BZU9FeFSsPmpD57RJmmFWWCm4lhFup666s8DLMPJ4Eq3JOXli4o4LWOQ9dZ1k+d8HwUdFwGhclr9XUieuZstz1aZK8xEBQhUy/nZ5nhb0Brew1M9D46Jn1IeyNY3iu9ysS9uZiCLFqxzvqjIN5rOs2UXyfw9+mohFOOuQtoRzANhSa5c7vdEPfpkPcJsZCIKu7TYS/+fDdefk5D31bdZYh9E08P4IgCIIgCIIgXAjE+BEEQRAEQRAE4UIgxo8gCIIgCIIgCBcCMX4EQRAEQRAEQbgQXDjBA5ERDCyT3DWVJVCWUc7TRBnPINXaynlSkLvWITfSFwDUlMRnlNdlCp9lOc9k4iep66n6B91EXrU3n/3E2CN13f2ssyi/TjKeSJKgYRlFqWtE2WtPKcE8ynqmJN8iHBSFeIyZQQgyqhif7N/mCwWUCt5QlrjmAuCCwSbI6xoTJK4L40JtHxXkcrt1fVIdnSR7DbSJzCmpuitxCwTBAwB7pG1PgzXTYN2MsalHuKIH0OTjCebh4GMxjfB3eXh4EIqone4UofAGUEBBHl45eFYotQPHZGxjHJzRYOPBRZDaVVHK1hVBXtfXQcZ8EWyo3lzrnXep63mPwyJIfdF15eZN2xe1YRjjUJggbV3qIHVdKoeCPCplJySuCwrvNU3KXAfZ6yB2oMljXYV+MTAl1swpXkgjXTGU7jUhCTYAyCIo+RoRP2NFWSJfu1gTSAFK6SCCQBreK7BXcIZzX0zH8rT7IvUqoGrHRV/qPC66ojMumjgmxjIPDISxT7XjIDDj9awyD51L55mUOFKAT2NzR+46l3hAEMYJC5Fr2iVBB0Yred0dD5PQQWjPcJyCRjZBlfEgsQYVGuRLUBXGRsX1iY+LwvE4i/vTC2f8XBSWQU3jROkWctNRNaZT64ZiCRtQVL7pKNlQvuiHAqfd2gbByOD4OaAaCuo5DjPVje6UrpLNxDP2WT6ldtNVe5su9LbsSmHOtXU2iACvCE6FmxbyCpoYNtX8icUeFDEa16q9pWWzjJ+uslO6KUqvD1qWbpC6ilWz1tHkJ27qZhU1BWJh085n3UKn4bsnb8QmCk3mu4HVRpSWlhvvFbz3nQKnXYMhFDgNKoahwGnqG44Jes7rYu30nj653/OsdfbbDmj76rTSGzBZ78vNmBjp1vYB2vo/zIB1un1eAbW3fKOflNNUW/cm19GJ41wqotytJxc2nv5S7C18elIwoCxisdmO0mr3x6aV3dKy6dccFd6SYZTGwc7r83YdOk9qxIu+Rz0fI60gCIIgCIIgCMIhiPEjCIIgCIIgCMKFQIwfQRAEQRAEQRAuBGL8CIIgCIIgCIJwIRDjRxAEQRAEQRCEC4EYP4IgCIIgCIIgXAjE+BEEQRAEQRAE4UIgxo8gCIIgCIIgCBcCMX4EQRAEQRAEQbgQiPEjCIIgCIIgCMKFQIwfQRAEQRAEQRAuBGbROyCcDp9+088uehdOFgJYAawBJsBrAArwZvIzIHwOFZ/jZ+E7OC/L39l9D8Azg4igAJA/+T+DCSCe/M3wQbtPmP48/R1xHfJxdR1f63b5MqO1h1IehXEgArTy0MRQxNDKg4hh4h+hlYei8AcbFZapeJBU3CZhfWhgD4Ln6QN7in9PZx8mliP8XZ7b9x4KTfxcIf498e9QxFCIz2rJG3FOWKbVlhqlPJTqnHedPpXPTzA0OJzPcdl+5/wsSu1Ofsc7pH6TX8d989xeB7wOz/kaET9jRbDxJFVOh2UMUOc660nBu+U+kfO4p+JroB0P4zgITI5z02NedzzM33tal1EKYzYbDmN5OrzMyKcWd8c6ysvIx7GzMz56imOgCn8q02Qb8vyn60rw6Tf9LF7+v/7IonfjRFj0PeqFM34evf+V+OUnf3vRu7FwHr3/lYvehQzXNVDXIBWuXKox8YafoGsACEaL1gTvOxdMhXChVPGC5zsXP91ey1khvqG8bRoAwuKzuUJO/AxPPueLevc1x9ediz75djm58F7Z+NwwdA3ohqFqhmo8VO1AtQXqGhjX4HENXzc4aVTjgdpDFQRlwz75BqCGQAXBWwVrFZQiMBOIGFopaN0aO5oYNh6kZAwpYjQuLKPuDdqU8cNM8HssytNh11bYVj1UqsG9egseBAcFDwXH4bVjgmMFxwoNGzSs0fjwPPYGDSvUXoeH0xg7g9ppWKvhLYGtAjXhWFITjqduAFV7qMaD69O9sdyPHT+aaz1azO6dGfMeh0WQ+qJuNHyBTn9MfZFgrQYRMHYm96ue1lDkMfYGUIBiRpNmVgAUHM5yUDi/QeFuW8e7z11fYdv1sGsrDGyBoS1O/28Fg6YmSbLhw+FaA7QTI+m98yobPI0Pf6P3BOvCtcQ5Be/DsSJLuS+mx2n3RR6NgXENEIEUQdUuDGwE+CYOeGAoS2CPMNGiQpN4xPEhGgF5TATCmMftGJrm2iagdqw61bkkD6h06NIPMSaMoD1jIeIY6NsxMI2LyrbLVRO+m9zkuKhjm6nagRoH1GFM5NH4VMZF4Xicxf3pck9rCIIgCIIgCIIgnBBi/AiCIAiCIAiCcCEQ40cQBEEQBEEQhAuBGD+CIAiCIAiCIFwILpzgwWGcJzWN8wYnRRuNrBSTFWMIE+oxeRm1r7PaG9L7uG5X/YZDEq1HEls4wf2Pv5HV3jrPnMQXuq/TKiqq1kypvVFU9OEkhLDkKB1U3UxUfdOKYaKqG0XFt6TglJKw1ZTIQVKn6mLVpJJT9/Pp94ctm34+bDsgqdB5zEqX1eShyEOB4uuobhefzcQyv1cFcEXp5MkLywhNqg2m81B3luVzl/zM8x6I5/felHkoYlTa7umL6bPussP62UHbtQp1vOfzafVHG0UOAMAplQUQtA/93akgxOJ9UgNQWAXdjjQeojPWTT8DmBjz0vuJ8bDzfFpCB14Hpbes0gpEMYO27TiOZcQIAhUd4Z88Pqr4Pm4PCqIPRGE5xzHzvAuvrCqLVnoDTsH4+Tt/5+/gJ3/yJyeWvfSlL8VnPvMZAMBoNMKP//iP4z3veQ/G4zEeffRR/PN//s9x7dq1k94V4Zh8nj+Jx/HpiWVr2MT/hR7N73/8x38c//E//kdpwyXm81/8IB5/4lcnlq2ry3g1vg1A6ItvectbpC8uMU//6vvwzK/9HxPL1vr34OGv+X8CkDZcFfa7pr5G+uLKIH1x9fmDxz+AP/zDX5lYtq4u40+W35Pfy73NxeBUPD+veMUr8IEPfKD9EdP+zF/5K38F/+W//Bf8h//wH3D58mX86I/+KP7cn/tz+O///b+fxq4Ix2Qdl/DH8afze5qakn7f+94nbbgCrPfuxatf9EbQ9iBIp9YWvLULAHjLW96C97///dKOS051z3V87aP/L5gho/eMRTH0wCBMaUobrg4HXVOlHVeD6up1fO3rpC+uMutr9+GbX/LnoXbHwKgGtgfAuP1c7m0uBqdi/BhjcP369T3Lb9++jX/xL/4Ffv7nfx7f9m1hxutf/st/iZe//OX46Ec/ij/xJ/7EaezOHi56rZ95NNQJhIp6e5bbGNzz9/7e3zuxNuTxGKgbgFTwzNcOpABFBNXESm0EOAsobt3iKeQrhcMh1TYgADEMgNAJayOaWeytW+vnNOsa7Ffn58C6BlO1fXJ9g26dnybWMmhCjQ9qGKr2IMdQUOhxD/AW7BncONTxAP7rf/2vT6wvUuNAtYcuFVTDUIagTKi9wA2BNcE1GqxDqIlSBEuA052CpzHcpFv0FEAufJpC4TwTjPJtCItXucZPN2wt0S1+OB3SNr3eYdsZ8thxJSrXR081aKCh4dvaPlDwHGr+eBDqVOMnPixrWB/q+livMXIGjddoXKj346wKNX5qyvWbiBSqYhNFzehTDc0NuL594m04D7d5vloY5z3cZN7j0GXWNdXG7zmNvqhqhioQ+mIDKBP6orMKzio0FM47rTxGjrHmNcZg9JWGJkbjPRrVxi+WZOGhAPZwFJ4VFDSF83/b97Hl+thxJYZNkYuIAvP1rf2WHbQd0BY/TmGz1qu83aw6P54JLl5rXHzNMcTKOgXvVajz40Jf1M3Z90VfN+CmAalQvIeqAkQEUgjXd/IAFFzZVvzkFAqmMRFSncLbumNhN+oaaGsA5dDqzr6cxpiYrg9hjGtr/LTL4uvu2If2da7vE8dA8gDZ+Fmqfdep8ZPGRfIMYoUe1sJ3eMA3Q/h6DMvBAjrJe5vjcNHvT8+KUzF+fv/3fx/3338/er0eHn74YbztbW/Dgw8+iI9//ONomgaPPPJIXvdlL3sZHnzwQXzkIx/Z9+Qaj8cYj1vTfGtr6zR2O7PqeT8nEU85wA4+xL8EDY3LuIoX4+vRozVs4xYA4Fu+5VvyuvO0IXAy7djN89nvdY5lpr3vc97PtBGEznuEG21ydKK5NNwpHrff8ul8n4k3HMcnj4mB6qDBaVA/i1/93D+DYoUr+l68WL0S6XbmzPsiMbT20NrDRKOnUB3jZ0bOD4CcFwS0NzoT1dw1TVR1PwsU9p4Ymhhuqn01PDQYDWJeAtocCkU+/G0xX4H22f3xc0/jt3/hJ6GUwZX1B/CS69+GXjwxlvZ6KlI6e5h1TTUIhUAX0Y45rw6cz8f83Dm/NRi6+z4VIp6Z/xbWK7RDcQYZM928nzRJ4lVr8ABtLmD3GuFYZSPIEYM5rE/EsHG3vd/bIevby9cXU15rfsT5wsl8V0yMe91L5Z7L5inl/YR8H8B3at9yTPrpGjfcMXbYU8z9aQ2hnPeTTsn0d6W8IIS/e78JmMH4Gfyfn/zHUNC4UlzDi/mPoYJZ+L3NvMj96clw4kPUQw89hHe/+9143/veh5/92Z/F448/jj/1p/4Utre3cePGDZRliStXrkxsc+3aNdy4cWPf73zb296Gy5cv58cDDzxw0rstdLiMq3gFXo1vwmvxMnwThtjFb+D/hOUGdfQPH7UNAWnHs+by+vPxxx74M/jm5//f8HWX/ycM/TZ+ffzL2XsnfXH56d//Ajz/u16Pl3zHD+Or//j3Yzh6Dr/xmf8N1gc1DmnD1WDfa6r0xZWhf/8L8OC3SV9cZS5tPoA/9qLvw6u++v+Or7v7EQztVhgT5d7mwnHinp/Xve51+fU3fMM34KGHHsILXvAC/Pt//+/R7/eP9Z1vectb8OY3vzm/39rakpPrmMwT8nYPPW/i/SW+ig/jv+ImvnRHv71fO3LdgMdjUJz6VlWRzXJVq+yt0SUAUBv2puJEj4rhbn5K+S3ODmVvCgGsqH0dZ7YoeoSAKUW2E+RIIW/ozHIhuvY7oXDpsxSKoRuGsghhLjVDNx6qcbhv7UUgY4FRjQ3SuFR+G/5/g/8vvoInj/137NeGqnagxkPVHrrR8AYh5MYC3hLIEXyjYInDBJ2PykqaoBRD+RA+Q1GFKs1KEzEaz9njk0Lhup6fNIvbJYW7UOfA77fsqOsMbIUtZVEpi23Xx6YewnEIe/OsYghceOSQNx/C38beYOwNaq8xcgVGzmDsDGqn0VgN3yigUSAXQm02X/hymCHQf4Zh9HXc9/L78N9//R/g5q3J5PmTaMN5uOXnk3E772Fv8x6HxH7X1FPti5ahG4I3ISQo9UU0Ct4oOK1QOx283MQYuTAlP/bxtkABjTfQKpz7Bbs2xAweDh5BEU1h2/Vxy61hy/YwsBWGNnzXPP1p1rJ5twPaa0GhwkmXrgXdkDcghsd2QmWTt8d5Be+pDXvjEBqY+6Jt+2JvEyiGZ9AX2YNH43gMCFSXYUhUgG50eEEMZWPYgKLJCIj4fcStB4g6Xp/8GvuEvk1FI5wkZNt9yz8wK+St48npjn8T4W4p/M2342EKBc/jom3Hxfs2vgbKWFBtAVPh8sZ34kPP/TvcxBfv6G866XvUixz6Ns896klw6lLXV65cwdd+7dfic5/7HL7jO74DdV3j1q1bE9b1zZs3Z+YIJaqqQlVVp72rwj4UVGKdNzHEDi7jbgDArVu3cOnSpbzOYW0ISDsumoJKrGEDQwTBA+mLq0dh+lirrmJobwGQNlxV2muq9MVVRfri6lOoKo6Jcm9z0Tj1yOydnR18/vOfx/Oe9zy86lWvQlEU+OAHP5g//+xnP4svfvGLePjhh097Vy4EpxFPadligB2U6GETVwAAv/Zrv5Y/P9M2JEzm+qSaPwbwJrz2Zu9n+XPDsc5AeOaC4YtQd8CbzrpFePCd9JD94qo7OUbczU+aisnu1mroxnOn116Hx7xYbjDALkqEi/Qi+qLSDKM9CuPyo9TxYRxKY1Eai55p0DMNKm3Dw1j0iwZrpn30tEXPWPRNM/GodLt9v/Md+y2bZ7tK27xsGg0favQAMRE5Lu/UR9HkJ+r8dOuraBU+0+rwBDPrxhjUz6HU6wAW04bzkPveOX3cKe01dXF9EYR83k2flxO1pxDO6/w65q0lUj5QNy8o9ZduP079rNuvu/2uu2ze7br9vqftxPWhXzSojN3T/ytt43WmvfZ0r0dGeyh9uOt/mfpiGgv25P0kMSDqLJ8ab7pj0DQn5fVhlcZUnhxjTVxWxHE5PqdaQBNj99SYnz/fp+7fPLRj4hLc2xyBZcmbWWVO3PPzV//qX8X3fM/34AUveAGefPJJvPWtb4XWGm94wxtw+fJl/NAP/RDe/OY34+rVq7h06RJ+7Md+DA8//PCZKWkIh/N7/Nu4F/ejhzWMMcQf4FMgEK7jQahoL//Nv/k38fznP1/acIn57JPvx339r0HP9zFqnsLnRr8OAuE+PB+fw+/iz//5Py99ccm58Sv/O6488HVQ7i4MnruFP/r0+0BQuL7xtfi9Zz4kbbgi7HdNlb64Otz4lf8dd9/3dVg30hdXld/7wi/jvvUXo491jMdfwee2PyL3NheUEzd+vvSlL+ENb3gDnnnmGdx777147Wtfi49+9KO49957AQD/+B//Yyil8P3f//0TRaTOmsNiKlddUeNOGGOI/4HH0KBGiQpXcDdejW9DSVWWZn300UdPrA3ZuSx1DSKgcSCiIHVtQ0Avk4JqCEwhIJk1AN0qvgSFthjbHOOcu7k1bWBzCGburkOYPcPFhuERZGKPrPrW/f09Xzz5vCeMfSr+eSLfJ71Osc02xjfXHGKbbZC5Vo3HaHwbv/3sL6J2Q5TUi+34rXm2+W1vext6vd6JtCPVFlQ7qCJKXReAajoSu6aV2DXa5fh9ZgIjHOfpCfWZld4pqKZ1K7hPK71Ny+bOkraeZ9l0Zfm0bGALlKpCXzfY9j30VJNVrnx0FXq0uUguvm68nthf6xUsK1ing9y11WCrAEv52LmtW/ji+/4N3GgXRbmBu9YewMMv/AGUI5x4G87Dti/nW/GEc+aWjbmPQ2Tfa+op9kXdMHynL7IJfRGWwEnu2mtoFzyQ1it4NdmnHCso8iGXjRVAPsi4s8qeoBGX2PY93HZr2LUVBrbY06dmKTEepS/ut93EcmpzR/yMqf+83dRnSeY6SWEzI0vOh76I3Bf/8Hf/Dez4bPoij8dZ6jqMjxTGxcqDKSS06jq4bJRtIx7ymIbJ13nMm2JivKR9PssLcKy+7Yswnk54TWfk+bTPNLFszzop99W3Y+Fk+Qdux8Xa53FxPLqN33n6P6K2Q5Sqh7vUvWFMPKV7mzvhPN6jHuaxOqt8H+AUjJ/3vOc9B37e6/Xwzne+E+985ztP+qdPnFU7uU7KFfr1dPgsxz/8h/8QP/dzP3civ3ccpt3crKILPLnydfc1z3T3Z0lPNXU1Txf7GPLGTKD6CDvXCWvrDhRZ0ro70HQGl+nxhzH1ZtqYO2QQeuXzvw80roFxDb69BT8cgp3LF/lF9cXCOGgVkqxL7bKwQaFcGw6mUohN53VH1lZFg8nGJOYzlbnOIUGtEEP7mQ83jOAs9KvB8LNC3lS44Uwy37OSuR/4n38gCh54FAOP6uYQemcEO3oOwPJeT08iNOw8sd81ddF9sXvumRnhb5qCVHuiK9feJchlt9uvmaPXQTou3f6UrgvA5ASGZ4JVqn0drxvMBEUanj2cD4IsyhPG42LP7zzwP/8Aes8yigEvZV+cVe7hQKnrjugPgMlxCzMMpanP5zaCCDnEPIeRd8dFAJRkxRlR6po7r+MYHMUN8raq/X4QshDSQfv1DS/8/jA5MK6B29vg8RhuZ2dinUXf25xXli1U79QFD5aZi6aocZZW9VHxozEUEUgRqI4DpwJUHV08iuEsB/cOcZjJjDdY2fjhTlxzXDUovtHExT6pu6VlSWWOpy/ud3I/PW1P8dRn016fqdmtbk2DNKs18dojFDVN3p6Gg7pN7YLaU23DcawboK7hR+PgYTstrIWqLbjUcV8IvkCrMtVQVpmqGwOlQiHT5AFSYDitsiGUanZo5aF8x+iJ63a9PtafbVGZdKNVKYvb5RoKsrikR3BQaNgE1TeEgqdhWVvkNCi9GdRR4W3kDEbWoG4MrNNZXYoayscuePfiDGZtgSY+FsAtvzbXerSY3Tsz5j0OCyH2RVV7qIJykU4fzyWyBO70RQAgYtROQ4FR6/a2oNcpcuqgoOHRsIFmDo3MHru+wm23htt2DTu2wk5TnulkBND2ye71AZis8ZMKm2aFt2j8NF6BmeA8wToN72mP0lvqi7o5u77IzrXjYt20aqi1iWNVuNayCvsGIng/OQmYPTcxB2ha7W1iErD744S9am/cfna8P4ja68KM8a8d+2hiHOyOfd0ip7mwaScKglynneK4qPYZFzmOi8JycNb3p1KKThAEQRAEQRCEC8GF9vzMw6qEvi2bS/G06aqeJVe/j+UPJhRuYghc9vSodvt25osPdPF7hJj5I8U6T82WcWcmbUIFrrMohbkRd2K3uY3d3hP2dnaRJSdKYVzw7miPIobcaOVhyOf6Pt1Z3G4YDrA3B6f2eiL8zIP2hKPNu2yatE533entQvibB6Bj/ZPJmC9FHorDrLnG5N9klA/bz6H0tkpI2NtqoNL5132AodGGdYZwtkkVw6zwFsM/gaD2phBDWQ/pa8fto8mjM73MkM/XDGAyLyg/dMdbTMHb40Egp4PnJ0UAKMLoGMdyGZhQe6O9r9tQ77hBN/RtRsTDvuPivFDMRToo5I0xmePjuf3Qh5A3MKBcuyztm+rEibNqP7sorMr96TIixs854KIZPgAmDJtu/s+0FDaIJ/J/EPN7utKeHA2g9L0TUtSWoDhcwJXF4QbQfoPFVOLo9GrTuT6Ufr+b38PtsiMLMCwRSoX8Hq0YpbHhZku1Ny/pRgZAfh9e+xlhbyqve6Z/A824KURIBtfEUMzxPUGD0Uxt082PCHkX58sJL8bPaqBj+Gkrw945P7OEdTSGonGTt+3caSpq++tZ5vvk3wfDpJzBTtibjbUAPBNsEiOJEyZZDCWun8LimiMWsF0qKE4CJqNGdcad6Ym/qdf7Gj3HDXOjVsraawCG2/Gz88xMIccnLc9GUihoTqEaLTxRKFqew/w64yMAdjGc74KxCgbQMt6jXnjj56Lk/Sxzvg8AsG3A43DloqoK8c0KUHVMPCUEZRsCmIICT1ruQ6HxcCGN4jieU1xziIPOeT4qfAehnf2iOH10x7Nce/6oqa87St5P53U31hlJ1WZK0UbnOHQPakJsM8Y1UNfg0RhsT/mGZFwDPdeqTNX7q0w1VoM1oDmoQKW8nkK5fT0/ACZmp7tJy934fkWTB33Wsmnm3a67LN1s3bZrKMihiMHsDgqOqZPnYybyfSbzfnTI93Ea48bAdZXeooqfijkGumGoxgONCwIW48XEqt9y63Otp857zs+cx2EhxL6oGg/dqNAXG4qqi6FPOhP64rgxuc/VTkMRY+xN6GfMaFhDx0zypmPROlJBGpg9bqV8Hxfzfe6wP85inu263p/k+Un5gNPqkDbm+VhWcD4+mGBdUF+cVl3Mj/ps+2IaF6kOijtEBGocVLzL140CCHBlmCFjouztyeqnyRDiyXwfjuNgnpjrGBLdoW86F+g4EANo9iYX5SbtCB5MjIEx/yfl+eTXts3zgQ9eofaaGfJg23HRteNi0xkXx2cwLt4hF+X+dBFceONnHlbBsj6IZTd8Er6uQ05mXQOKgqOmsgAZgAiqUdlAcQ3lWa6kcJoumlkhJn5v+zpNFaH1phByTNoeWc8I61iYrZuwOYsDBoY9hs/08u5FH5MDwIQhxO1FPyR2hoROanxOwg1qNjZe4M/oRrluQHUDVei9idYNoArkROvGamjl4WNbOa8AFbwkKkzytQnThD2hMAAmDJ9ucrX1auYN12HMu11XLndgS2y5PnqqQY8arKtxFDtQWezAIRhCPi63XsOzyjLXjhWsC8fEd8QOstBBSt5Nhm0UsOD6KPKDJ8e26y3kd5eNpT4OqS/WHqoIfTGdS75BENMo2r5otENDOsiuexXPTx1UClU4fxV5OKgYEufhmKARbrx3XA9bro8t28velcQsCfr9ZOn3Y97tkvpjqd2e307v0yOFvAHh+pMdD0xo7JTwSNMRHrF85n3R1zVoXIe9VQSqi+AcIROu+6SgmzDbxxRkpaHaMZGBPN7l8Gp0XqdQ62OEtB0UCREKjDO4K5w3a/30w12v0Iwxr/uaupOBUfI6Gz82GKl5XGzacRF1ncdFv6Br6FGZxwBa5nvUebw+i7hHFeNnTpb15FpGd+KZ0sn7mc75ybk/iieqXO/J80kzY8STBkzntS/ClV4lyc1pZg0cncEh38t3w9g6Yc7d0LZuzk/eP263PcbE6dKhVMgNKLQLnh8wCt1KXZcxhAUIM7r5tXJ5meWzl7me+Buolf1NNU80MTwjh7z5eMMYZtJnhLyRP5ehGlP3wMISQxRVFbNndeo87YS9pdC3iTyfGCKnidHTi3H55Rw68vkakcJiZ9UEmwh7I269P+ch/DSFg0+Md50Qt87rmeFt3ZA4zPh8etk+45EvONT30Qw2nZVmhbx1vyN6e0J4d5gNYxXHXR/2S8X9SxHP52FMPI8s8/3pOejpgiAIgiAIgiAIhyPGD+Z3uS2bFbts+3MS+LoG102Iy22aEN9cu+i69jGGN8Vdo9X7d1ENJrnBJ9zlNBFCFh6T01h7FGemYM0xefP4f9tBs1PdkLdZy/Ozm8wFUbGCta6De58aF0LemgZcN+DR+Ezc+9ykNnMT7ZTDEJqomGcJrtEh74f3FiJMTHt0ptWbul6fmcpOd/CY97tqp7HjKty2a7jlQt0XF7fvhr3lWWe0IUXWhxAjxyH0yDUhz2BWvk9qY6pjjkEd2nYR3HZrcz3OO8t8DFJfpLrTF5tOX7R7+6KL52N46Ogxab0mDm3dKg/K5/mIy5zvs2vLE+1fR/muFAY77Q2edU05aB9SX5zIvevkkqTcyrPuizwa53GRxjaPi7pO4c6dWm9pXJzKn2m/rH0+dEya2Ikj7K+OKm8HeH1m/casnNc0Lk+MhTHPZ2L8b9ow8NRO4fi4dlxMbbZi9X3muUdd1fvBRaVlSNjbirKqJ/ppkcLesuJbUntTk2pvrNucH1acwwBal38b+pZD4QCkeYJQSJUm1dYOClvqfMazF7e/NSPmOcVkp/htnKNEcq0YlbY55KZUbkLxLYXgdMPeihT25lLY2/7zN4r8gZ+f2N8xJf17WNibiX+nYQ/rV1i27wAk7G21SCIj4eGOFPaWQj5nKb2dZB887LvStaKKoXeOCV65bMBZ0u1rr2ApGHLKMRwpWKWwOy5PZF+XAaZggADYP9xtRtjbXOFu05937RsVDR/DQeLapMEsfFF3oo+nJiVT0evwIQBHSImgHhRChENVgVCCIk12CgCWKz1j2e9RxfiJzKuqsUwn1zysithBl5SgT0qFJEUARATV6JyroyzAFFSM8oW9I3WdLogck2SomwCa82umrtrpIs/d12F7jjKdHoAyIQl2Yr1D/6j5Pt9vti2p3UwogNWdKtZxNrJN6mzA47Px+gCxzaoaVJfZS6eKqJRUTCZaextmWGunURoL4jZ5uSt00CXH5sc4ftuZmd4PH9WpppOl52XWdtPLdm2FHd2gIIerZvdQz49lPeH5aXyYefdRCS+JHVCTZjJTG8eK8kmlyM1KPDt9btv+fCue8xj8uY/DAkh9UdUWXOro/dFBrr9Iwgd7+2IyCKzXWa7dM8FT6/npCh6AFAa+xI7rYddWqF17O5H6iT9BzfPDv0vBAlBeZ+nt/fr9Hu8PKIsdeEt7hUeajvDIAvqir2soFY2HugbF16o2Qf3UckjeorCv2cgxM4yZOL+Xh7+pz+9I/CD9ZjR6DjJ89kZkUFtxovu+o3A6/VA2Kr1llTeEcbHujout2EESAFoVsYPjsAz3qPMaPou8PxXjp8OqGECrcGLdCX40DlKeikDjKqi+JeNHAUwKqg5XWtcgX7x952xmNelC53TBZwKBc9HRZByFjYCuokBrKMVlJnyPL8Jvq3mjHWa59md8fqDy24SaDbc3xTVDNS6GuFhg3ACNDRf40Rjgs5kW86Ma3GuCKlHtoAoFXaqwj7a9ieCG4KyCizK7a2UNhfYGBBxvnBDrb7CCj7PP+beYQmjOnJbncQURDlJ7SwxsiS3VgyGXpYBdDr0JBppjFeSuk5xuNN6sD7K6zqoJw6cNT2pDPJNhy3UT2nVB7LhqYb+9TMx3HBbjps19sbag2kAV4XoZDJ84IWH29kXrg/fDsoJJYWCsg+S1N/Bk4RFU3zwpeAZuuTVsNT0MbDnRHxchQuKZoJiipxjZS9RVhJwQPkBH+Y1DnZ/UF/U+fbE7yXSmfZF9HhepsQCpMOdXOygFqFq1E4MlwNEQcgXyxB/QMXowNRF4UHPNO8lHISzcJ29PMnzSxxOhd1NfmAydPPa1hlEY/2hiPKcY8tYNEdZRijyMOT6Pi9QdF+v6TMfFk2RVZK+X3eOTEOPnmCzCADrKSbWqhg8AgH3r/elF2WsiqDpIXgOAj/UNQrGzOKtUYKK4Ikf3flZ964S6MSG40zvPOcwtDQ6E1gWfPjMh3p0Lhm8IZjDHqDB9nZ/+KIUbdIyfHOYGgLo5IJY7uSAeuunENCevzzjU9jnTC3xsMyqLvTPOdbjZUhbwloBGwRuFxmgY8ih0CANLYW9JsSmFixX5dQx7UUCJduY2zTLv8R7hYO/NUbbbbx0AGLoCQ1dg11dYV/PfDDGHvCHfqCytq2xs7yhx3dZvstmoXeTAvWXnlHg++3vfM2W+47Bz6vsxk9QXG5v7om7C9crX7fk13Rf5iAbLrq8wcBWGrkCp7aFe0m4fmu6D3fXvZDsgKEIW5KGUDRMQMUfJeg2D8KyIg7x9uho7YODK2X2x7vTFZoF9MY2L4zodmJD7owDd6NjnFFxDYAqeINUEg4RLzA5364ZlzwiBOwq+ANyaD6qrhQcq35nQa70/IdSNJwyepGiKaOiAY+mJ5DECgz1NjJcTkRB2MtdHdXJ90HTGxQVfP++Uo0zQA5B71AOQyGxBEARBEARBEC4E4vm5AxYd/ibMJtX9Cc9B8AAqhASwTt6edhkUt16eGCudA6PRWY4YTlDERHYOs12zCp8eWHdgv9m2iT+i9RAdwZmwkvS0hYo1QoxywRvU8fYUneRrzwTDMScBPBFqk2Zxp5dNh8cdtF36bN7v3rvMB49kqgEEgiaPghy8UijIw1Lwcg1xB9KBS04WFhFWijL1v/xw4fxVNvdBHev6JCptj9zP9lvnsO3m6Z+lsijIoVAejW9FUcZk4EFoYo0wSxpeBcGDxp1cbtKyMiF+cAfMGrPYxLo+yetTMEjnxNvoweE4rMaQNjDgKXh8fBvqzQ7ZAxRC1EMIXz7lFpPuuJKc5T3qqoS7JcT4meKocZVn4V5ctZPqJGDngLoBmibFnwUXP6VETw+mEJsdrqgUQ+CQY4+7rv2UAIpU9E0xWAejJ4e/gUGqDXNLRg+lkDeEZRYACg5hIxRyNfRwckSYqaZzFKjj1rftQzeATsmcKaEzhbw1UeJ6PF6IFDLHhNKca1CqsI9FkNZNuQaIyf3OKoxsgX7RYEPXKJXNISuVChZleB8Knmp4jLkAfFCKWyY8h+TvWWFv3XwlIBhGtdcY2QLjxoCtaqV1s6BFSuANoTYhbv30q8kfxk4zX9jbeVd7m+84LCjsDaEvYlyDChP6YhJGiX2Rzf59UU3N3EyfvwDgEc73gS/Rnzv58WxQ5FEoj4oaGJ3kutsT0pDD2JtsQO3YckZf3BtSpSx3ckgW0xfZuRDSqOKYWJWASuNhsBB0nQatkG+Zx8SyExZ+nMvnPtu4flB184aBvgeMhy58OwPIBO4UMmWfwt+oTbFNYW+d3B5GDI2L43cqBk08NS42MQc2hbw1PBkKHstmLFIoZpGchQF01HvURYe8AWL8zOQ4iWWndYIdx/BZhhPrJGAb9PhDnC+B6jBTTkndRjFck2aG4oUeCDPvRQwdVsgX4dYAYrCKhk72/ERDJ3qBQGFZvuDGdQiA0g2cVfCNgiMDaoJVpeqOBPas2bHpZbM8QPFZj7uDb3uRVxN1DByoSWo2Tc718aMx2C7A+HEOXNdtrsFYxX2l7CELcfQEjrkGtdMYNgXWTINKW5ho6ITZZB9ncMOAVZCDcwoqlvUOMtIHx2/PWmee7Y7zXQNXYVfVWFfjCalrABNS14o8hk2B2gWlt5RjQA3lY5Tbu/ZQ4zZufdGD98DO6a0658bP3MdhQYTJoxpoSqjaZvXFlMjvD+iLqt+ep8kQmpa6Tvk+ngl9XZ9pP5tnu5QrWEQ3QeM1lPZhHR9OTk2MxuuZfXE6905N5d4tsi9m4QMgXPcVQE0BRQBIQzUcoh66E4IIY4lHzAFS2Cu0A+wxcCY+o8l1WAG+ZLg1Dy4YVHjonoU2Hlr5LCaRPD8cn0lFASIOeTw5rycJK6S8H+ruBLXqbt1xMSq85THRxnExKrxR3YS2iuPieeC496fAyU/Sr/L9qRg/J8hJGkDH9fYsy4l1UvB4DCgFoqD8BqKg/FaHKawwywUgqa8RgCYMXMFhRGEGKXlvsuGD6AUKRlAyekj5HAKXvECKGEoxtPIgYmjFsC4UpxwrhrcaDgAXUR617lg1Bxg4SZ1nOvk03QRPCB0kRZukZtORXE3F27L06gIMn0SacUZZgJpQYNHXUzPODcGZMOM8bgxGusDAFrjWa8NVKgoz0OG9hY4hcM0JSueeNB6EHddDjxpo8lEaONRB8aBYt4hRe5NnmvcUVOxKmccEa2p8Ls63aHab+eqg8JJ55k6aeY7Dos0jroP6IsoCqvHQjQp9saGg9rZPX6y9QV83OcRNdwyhVNNnx/XgQVjTyysZrMDoqQYOhIIcRr4IwggqSuAzY+iK2X2xe92t276Ixi28L7JtwGMK4+K4DpOBlQUoCAKpSmVvj2vaCcF2TAT8WrJg2gi0rLIW36enPcZR+u41H5RP1x20cTClQ69qYHQwfLwPCnrMgPOpTh5Fg6f1/MCnfSDAB8ON0ntq9y2NiTQhbx2EKHTdETroFvtORU0XPC6eNMdVfluGe9RlQYyffbiTkytxnJNs1U+ok8bXNRSFCyfV4UIPICq/AcqqrOqWa/4AEzNevqQQ2hZzgEiHq6kyHqQApX327Oho/ChiaB1msDQxSmNhKLwvdYiLHzmDW70+Rk2BQVnBNQrcKKhdnWcPs3HTDb0D2hyk5JkiQI1UO/B2ZVaT0pCdNnw6Cm9NAx7XC69fkGacqS6hCr13xrlT74cbBevaGWfrNe4qB6iUxZoeh3wDePRUk0NU1lT4+3SU3J1+3u8zAPuunzhou+46s5Z133tQMN7ibLlXQeq6IIdn63WMXJhptk7NriliOzVF0mxzNGwXzWBe48ecb+NnnuNw+Qz24yCS+iLVZeyHsfZWUnvbpy+OXIF1XeMuM8CaqsMEBFn0qEFBFpoYG3qEDT0CcPT+dtS+O0+/Tkz3xR41Ue1NYcQFGtbwrHDb9vGVenPfvkgT4cbtjXX2JCy4L6ZxkXpVUH2ri05IuAIojIu6jNYKoZ0cpKD05wuG70fLg4NnhZIqW5KU5o4KG2IurWFwycCahS481tbG6BUNqsLirmqYZfwbp2FZwflwXiVDyDqdjR/nFdgTvFNgH0PiSIFdHCotAY7yeKqmJwSTBHmuhebCuDjujIuj81nX507vURdxf7pMk/PnPDhBEARBEARBEAQhIJ6fA7jTolLTVvIsS/skPT3LZFWfJH48glKpuFuIdVaVCy7+OgUuE1wZhA84zXIBSLWAfEreqQDSYRa012uglUdlLLQK3p3L5QildujpBleLAfq6RqUs7im2sa7GWFM17tbbWKMwk/SEvRu33Bq+XF/Blu1h11W4OdzEThOSaG8PemgaA9coeKeilwnQhQsFBgcF0BCoVu2sVkfogNLMY8O5crXOs5Ax5K2xoXL1ksxuTYTb1B66nAy3ycVONcHWBmPtMTYGQ1dg4EtsmBG+vvcEemRRkMV1PcRltbzhbrN42nk85dbgoTBigxv2Cp5sroS/MYbZ2NrsLaiY8wvagopoXAjdWIL2HTXzDRnqyuL39TSZ5zgs3PNT17kvThQebhjK7N8XB7bAmg598WuLG7hubgXvDxzu0wPco1drzvS2d3jGVxj4MvfFL9eXD+6Ltu2LqXYM1W3I2zL0RT8egcY9kFJhHCCCAqBrndUBQk4sJgWB0HqB2BG49EDhodYaaOPBnvJ4BQaKvkVRWJSFw139ATaKGmumwfXebVwyIzyvvIUreoArehebNMKAS4y4xFfsZsgN8yWea9YxdCWGvsBz9Vr0MBpsj3twnjC2Bs5H7+OgALwOnic7fW3sRkWEMVF3rpeT0RBhXPTj0YJaaLmZde85fY96nu9Pxfg5hJOsqntaIW3LdlKdBn44DBfwfh+oLJRWoMaFpEurQQ2DScNVKeSN4KLxwBRyaHwBOKXgHYGdR60YxjislQ368YJ+T28Hl8wIG3qMr6qewxU9wKYa4kHzHC4ri02lcZdaQ3KaNvwV7HCNL1nClq9wy6/hieZu7Lgebts+nhxfxsBW2GlKPDNcx8gajGuD0aAMykJDBWUJVAfFOF0D1ABm2Ob4FLtxAB57mN0muPZHFrQ7AkZj+NtbS6Vk425vQSOMs3q9iksNmmGIuXBDxBhABTfSGFHIjni6t57la28UV3Dd3MJVPcCDZh1qxZzUl5XHVT3EDUf4zdGDeLK5gqfqS3h6vI5boz5GowJupIGhhh4QzAjQQ8AMQ3sXQw+920DvjICdAeztraUozrc7qA5fCWhVE88p8x6HRcLO5b6od9by8qKf+tLsvnhr1AcQwjefLK8AAK6bW/j6agd3qbWV7IubfojfHK/lSYinxpt4dry2b180nb5odmzui7yzC7ckfREA/PY2qK6DAIKrQE0BrQlkPbzVYBXCqF0FsCJ4C7hoBDkXJgrZK3jLcMrAGw9TOly6PETPWPSLGvf1d7FmxrhkRviq6hY29AhX9AAPmGdwSY3xfMPYoBIFGYRx0aLhEb7sbuK2N9jyVZ4k3HE9/JEOE4UDW8F6hWFTorEa1mrYWsOPDKgmKBvGxDAOhuujjgaQGTL0mGF2XSj0PXbQO+No+DTAzgB+OFx4eOJpk+795B716IjxMwcnaQAJx4fHY7BS4fa4DhLYqvHR20PQTZKs7uTboBPrDAJbCipwUHCNAhFjZE1QewMwcgV0fL3jejmufMsPgvKRd2h4UsJ2xIwtv45t7mHb97Hjeth2Pey4CqM4uziwBUbWoG40nFNgF/KDdMr3iLk+WeQgq7uhrVqdZI/HLifepqrVy2L4AECuMF8WucK8qj10o1vVN7e3wvywKTAyBXZthW3Xw5rqY9OP8LQbLPovOhYDZmz7PrZdLyjB2QrDpsDYmslK8m5SWUrPUJZalpstV6+WB+7Ck/piUmAsdRAh6SowTvXFsTW5Lw5chW0V+uK2fw6OV7cvbnX64sgV8/fFpu2LPB4vTV8EYo7leBzGgujtoXHM/yGCqj2CHBzgSuTICG8BFfNkPQWTxTVxvYZRN6GfEzEGtoCikPu67YLEu4bHtu4BHtj2u2hor3dlmzW2fBXGRdfL4+JuZ1wcxXHRuTAes1Mx/7GTB9lgjwqmrrvqbr4Vo0iqp+Px8o2Lp8gy36Muo+EDiPEzN8t6ci3riXUasHNB5hMAjaugAFeboHZDgKpUFhEIYW+Uw96y+EFDWb7aWw1LQFPorMQ1sCGRWYEx8GW86Htscw/aezQUlMe6OFC8wPdz7YsdV2HoSgxsiYEtUDvTzm41GtykxNqO0tdE2FusXZBEDpK09bSE55Je4GfW/JlWfZtRZ2RgC/R0iW3fw6YfYddXeNLNqCK7Iuz6Cts+GMK7blY9kb3KUntq+4wXH2KTsas16y/sU/PngL7YVWDccRXWdOiLt3yBW4v+Y+6AQacvphvv/ftim0i/6No+hzFR+wfItX+ICKrhLK6jG8qy0T5OCKYoCY8gfsFE8AjjVBL+CcZPnBA0VVb/23Z9QAPP+nrPmAgAu1zkcXHb9/K4uGvbcXFs47jY6OCp6ho+U+I/SXI8S48nWeupcg+o61DuYQnHxYvGMt+fivFzBJbJAFrmk+o0YdvA7TTQWoOsDe7+xkGVQeaTKgVlg9qNK0PMMAiTIXAuKt0owFnCiJClq4kYo8Jg5MJs1yVTYqBDiMsNXNl3vzwIu77CbbuGG+NL2LUVdmyFZ0Zr4WaiKTAalq0i3CAowulhNHzqNtRN1UAx4JjfE8IuVO2gxhZqexQ8AeMafmtnaeOZ2Tn4nR3onY2sQlSsh7AI22/1v91QwSEU7d7utaFET1WXoMFoWOMLzT0AMLMa/Dyc5HbzfNf0OiHcbQPPjtewPaowGpYhxGaoQqjjsBNmM/Q5zAY7A/DOLvzO4oplTkM74vlZNdg58M5uyJUkgtkpEIb+2X1x1GuV7J6uNgAg98V5OW6fOwpH/Y1n7Ubui88ND+6LKfTU7Fio3XHui8t6Q+12dqAaC+pVIK1ATQlqgmSabzSoUmFS0BJcGUPgXFSTdmGcZKXAlsAlo6YCzgYVPE0eY2cwKgwMeQxdgUumREEO276HG82VffcrjYtP1ZdCTqyt8Ox4DbtNiXFjsDsqYWsDNzSgJtTK00MFVUfv23SoWx3GRTNw0GMPsxMn2EYNsLMLHo3D4xzJWs/LSYfA3SnLfo8q03iCIAiCIAiCIFwIxPNzRJbJ+3OR4dE4eHp6Taj9QwRVuxzephoFjq5631Cb/5ND4AjchO1CrHOY1aytySFwu7bK7v6bzeHaTY5VdOtXGLgSI2dQW4OxNaitDr9jQ2w5uVmFLTs1fSZyP2LtgsYBTZNDoZbV65PYr+bPfnVGGqtzzZ9dW2Fb9yZmmx2rXGTxKMzabp7vOqntkhdw2BRorD68tk8zWdtnmWabyZ7ubL5wOkzU/Gn8gX3RW8p9ccdW6OtmT188jKP01e66x91uHka+yH2xdvrAvtjNJVmW2j6H4cejEBLeC8XAQ/HTeO0ghNwmIjAF5U2gmw/b5v+gQQhFBgDSGBemDYFzMSycGLftGnaod+h+TY+LQ1vkcdE2YVzM9ZXc7Dp3ORQxFzX1eUwMiqehzh2Pxks/Lp42y3CPuuxeH0CMn2OxSPfiKpxUZ4Efh/AvrQjU64GaEpoI5Awoq9wQXBnc/WQJvoz5QBZwLrj+2Wl4ApwneKuwTYyx1cEIAmPkCjyn1vAE7pp737bqHnaj1PX2sMoqNhhpoCGoRsEMggLdXrd+yPuYULHZDa59jGrw9k64wA9X4wLvt3fCgAzA7IawNrumshqFGyGEKAJohgYpwOuZeh3bdvkVteah9gbPjfrYGvfQDA14aGCGKii8DQAzAIohwww8zK4NqkW7I/jtnaByuEQUO2L8rCLpPFIA9Hp7wzqrL9qeQQNgB8BzVR+K+Nz1xZ1hdWBfLHZd7ovY3lnKvjgLPx4B2wrUNKBeL6qiGpAzYE2gJho+NBn+5iwmwuKcUkBBcJ4wVIzGaoytDmNiDAt/rl47fIci1msMbIHdpsTtYQ91Y2CthhsZoFZQQw1qAF0TzDAonuomtIeOZR6KXR9CwRsPPWiC4un2MBimo9FSqfAtmkUaQKtyjyrGzx1w1ifYqpxUZwZ7+K0dqLoB9XsgraCaAtQzUdZTQzkApKCqeIGnIPdJKdHTIsc6+0JhrBjOKVin4fn4N3rWh8Th8aiAbzS4VlCjENOsGoIetXk+Oia6F7uc60qYgYUeR4GDKGmNOuT5rFI886ybrmJNI2bkwqb7CZ686XpCXTn7nT1FrFfYGVbgoQGNVGj3IWBG6WbLhzj2aPhgSW+2zGqKfQno9MXNDSQfzqy+6PoKjNAXt/o9DJpy1tetLNarbPjs2xeXfBLiIPxgAKqbMC4SgXoVVFPAEIEq3ZHAjspvJcXxEPAlQI4AUvCOwZbQEOAKBecUtOI7Oh+sVxiOyiBw0Gigm/uaIh9SmYcmtIceh3HRDIK3R48d1PY4eOR2d8HDEfxouVT4loGznqRftftTMX7ukNM+wVbthDprkgACdnag6yYkfFYlimYdXBXwpQY1BXyl4AqCcgqujDKvDlHuVcGXCCIIjjAuPGoCdjH/zNa+1Cor15jdSXGDMJsVDZ4xo9hxwZ0/dtDbo+DOH4/Bt7dX2p3vh0NwXUOXJXSzjrJQUNZA1UFa18ZBjwsNZwnjWuOprQocatbueQb2LjvsM2D/9RMHbdddZ9ayedYBALOtoYeEYqsVOahuuyxyQM9uBZGD7e3jH/BTpLq16D0Q7gQ/HIKeeQ40Xt+3L3qj4BqCs4Rn6kt31G9Oou/O06+n9+Ww/aOxyiIHM/vicwPQ9m7oiytk+CTSuJhFEHoVtN2Eqgy4KqBsCVcqkNNwFcEXBOWC18ebYAB5E6IlnCP4QodxcfvOPYDMADVhXNSDIG5gBm2YW7Eb6tuphlHueKixh649zHYNGjdhXNwKY6K7ffsEjtb5Ru5RZyPGzwlx0ifYqp5Qi8Tdvg3a0aCqgrIOVJZQvRLkevClhq80lDNwBcFVSfUtGkEl4AsCdjSAU1C0YsDsdma0dkP8cr64N+HirkZByhpbO0HRbTiEH6z+dDs7B//Ms1B1jaLQUE0FVRv4soBqCNaG4292T+n4Lwnkw8xmudXOMpfP1e0s8zPPLvXNVu9ZmV1dddytW1Dj8b590Rch/OlC98Vnnls5j88s/HgEjEegHQ1tHagqQWWJolmD7hkoV8IVCr5SIKdCIfCS8sSgLwE9TrpYJ38uqDqpm3aKeu+0Rb2LnU5R761BVjp1t24tVS7kKiD3qJOI2psgCIIgCIIgCBcC8fycMF1r+KgW9qpb0ssAOwceDOAHA6iqF9z9g03osgT3CuhB8AK5SsEMNXxJKLfONolbN4xiJyguFdthZovGDWhrFxiNweP6XLrz/XAIPxzCADCbGzEHaB1mXcMOVRtSlZqjWzePpt7PWrbI7WjqedZ2nc/KbZ+Tqs3N2yuTVL325dUMvRQmOagv6npqTnSZ+tlR15nj92b1RfvUV3CeYOdgn3kGAKCqHtRgE7pXQQ3a8HA9LODLEBZu1s92Xlw3bV27YttmhVN1O3p7RmP429srG/q9TMg9akCMn1PkPJ0oq0hy+WPKkFDxUSxkr/Zy0QKJ7FNfAeLNhYFchOyid+AI0Id/a9G7IJwgs/pif6F7tFhWqS8eFz8ewT+114hYlvEw4XHxxsaz5iLfo0rYmyAIgiAIgiAIFwIxfgRBEARBEARBuBCI8SMIgiAIgiAIwoVAjB9BEARBEARBEC4EYvwIgiAIgiAIgnAhEONHEARBEARBEIQLgRg/giAIgiAIgiBcCMT4EQRBEARBEAThQiDGjyAIgiAIgiAIFwIxfgRBEARBEARBuBCI8SMIgiAIgiAIwoVAjB9BEARBEARBEC4EYvwIgiAIgiAIgnAhEONHEARBEARBEIQLgRg/giAIgiAIgiBcCMT4EQRBEARBEAThQrBQ4+ed73wnXvjCF6LX6+Ghhx7Cxz72sUXujnAMpA1XH2nD84G04+ojbXg+kHZcfaQNzzcLM35+4Rd+AW9+85vx1re+Fb/5m7+JV77ylXj00Ufx1FNPLWqXhCMibbj6SBueD6QdVx9pw/OBtOPqI214/iFm5kX88EMPPYRXv/rV+Gf/7J8BALz3eOCBB/BjP/Zj+Ot//a9PrDsejzEej/P727dv48EHH8Rr8d0wKM50vy86Fg0+jP+KW7du4Tu/8zvnbkNA2nFZkDZcfbptePnyZbmerijSF1cf6YvnA+mLq890XzwQXgDj8Zi11vze9753YvkP/MAP8J/5M39mz/pvfetbGYA8lujx+c9//khtKO24fA9pw9V/PPHEE3I9PQcP6Yur/5C+eD4e0hdX//HEE0/MbKcuBgvg6aefhnMO165dm1h+7do1fOYzn9mz/lve8ha8+c1vzu9v3bqFF7zgBfjiF794uHV3gdna2sIDDzyAJ554ApcuXTqR72RmbG9vA8CR2hCQdjwO0obng5Nux9SG999/P27cuCHX0zNA+uL5QPri6iN9cfU5zTa8//77D113IcbPUamqClVV7Vl++fLlEzto55lLly6d6HG6fPkynnzyySNvJ+14fKQNzwcn2Y7HHVSlDe8M6YvnA+mLq4/0xdXnNNpwHhYieHDPPfdAa42bN29OLL958yauX7++iF0Sjoi04eojbXg+kHZcfaQNzwfSjquPtOHFYCHGT1mWeNWrXoUPfvCDeZn3Hh/84Afx8MMPL2KXhCMibbj6SBueD6QdVx9pw/OBtOPqI214QTg0K+iUeM973sNVVfG73/1u/tSnPsVvetOb+MqVK3zjxo1Dtx2NRvzWt76VR6PRGezp6nLax+lO2vAs9u88IG14PljmdpQ2nI9lbsOz2L/zwjK3o7ThfCxzG57F/p0HFn2MFmb8MDO/4x3v4AcffJDLsuTXvOY1/NGPfnSRuyMcA2nD1Ufa8Hwg7bj6SBueD6QdVx9pw/PNwur8CIIgCIIgCIIgnCULyfkRBEEQBEEQBEE4a8T4EQRBEARBEAThQiDGjyAIgiAIgiAIFwIxfgRBEARBEARBuBCspPHzzne+Ey984QvR6/Xw0EMP4WMf+9iid+nM+NCHPoTv+Z7vwf333w8iwi/+4i9OfM7M+Nt/+2/jec97Hvr9Ph555BH8/u///sQ6zz77LN74xjfi0qVLuHLlCn7oh34IOzs7Z/hXXOw2BM5HO0obrn4bAtKO56EdpQ1Xvw2Bi92O0obng1Vpx5Uzfn7hF34Bb37zm/HWt74Vv/mbv4lXvvKVePTRR/HUU08tetfOhN3dXbzyla/EO9/5zpmf//RP/zTe/va3413vehcee+wxrK+v49FHH8VoNMrrvPGNb8QnP/lJvP/978cv/dIv4UMf+hDe9KY3ndWfcOHbEFj9dpQ2XP02BKQdgdVvR2nD1W9DQNpR2vB8sDLtuDCR7WPymte8hv/SX/pL+b1zju+//35+29vetsC9WgwA+L3vfW9+773n69ev88/8zM/kZbdu3eKqqvjf/bt/x8zMn/rUpxgA//qv/3pe57/9t//GRMR/9Ed/dCb7LW04ySq2o7ThJKvYhszSjtOsYjtKG06yim3ILO3YRdrwfLDM7bhSnp+6rvHxj38cjzzySF6mlMIjjzyCj3zkIwvcs+Xg8ccfx40bNyaOz+XLl/HQQw/l4/ORj3wEV65cwTd/8zfndR555BEopfDYY4+d+j5KGx7OsrejtOHhLHsbAtKO87Ds7ShteDjL3oaAtONhSBueD5apHVfK+Hn66afhnMO1a9cmll+7dg03btxY0F4tD+kYHHR8bty4gfvuu2/ic2MMrl69eibHUNrwcJa9HaUND2fZ2xCQdpyHZW9HacPDWfY2BKQdD0Pa8HywTO24UsaPIAiCIAiCIAjCcVkp4+eee+6B1ho3b96cWH7z5k1cv359QXu1PKRjcNDxuX79+p7kO2stnn322TM5htKGh7Ps7ShteDjL3oaAtOM8LHs7ShsezrK3ISDteBjShueDZWrHlTJ+yrLEq171Knzwgx/My7z3+OAHP4iHH354gXu2HLzoRS/C9evXJ47P1tYWHnvssXx8Hn74Ydy6dQsf//jH8zq/8iu/Au89HnrooVPfR2nDw1n2dpQ2PJxlb0NA2nEelr0dpQ0PZ9nbEJB2PAxpw/PBUrXjiUknnBHvec97uKoqfve7382f+tSn+E1vehNfuXKFb9y4sehdOxO2t7f5E5/4BH/iE59gAPyP/tE/4k984hP8hS98gZmZ//7f//t85coV/k//6T/x7/zO7/D3fu/38ote9CIeDof5O77ru76Lv+mbvokfe+wx/vCHP8wveclL+A1veMOZ/Q0XvQ2ZV78dpQ1Xvw2ZpR2ZV78dpQ1Xvw2ZpR2lDc8Hq9KOK2f8MDO/4x3v4AcffJDLsuTXvOY1/NGPfnTRu3Rm/Oqv/ioD2PP4wR/8QWYOUoI/8RM/wdeuXeOqqvjbv/3b+bOf/ezEdzzzzDP8hje8gTc2NvjSpUv8F//iX+Tt7e0z/Tsuchsyn492lDZc/TZklnY8D+0obbj6bch8sdtR2vB8sCrtSMzMJ+dHEgRBEARBEARBWE5WKudHEARBEARBEAThuIjxIwiCIAiCIAjChUCMH0EQBEEQBEEQLgRi/AiCIAiCIAiCcCEQ40cQBEEQBEEQhAuBGD+CIAiCIAiCIFwIxPgRBEEQBEEQBOFCIMaPIAiCIAiCIAgXAjF+BEEQBEEQBEG4EIjxIwiCIAiCIAjChUCMH0EQBEEQBEEQLgT/fw8P7i4XyKqXAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x200 with 7 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "iis = [0,10,20,30,40,50,60]\n",
    "fig, ax = plt.subplots(1,7,figsize=(10,2))\n",
    "[a.pcolormesh(image[0,i].T) for (a, i) in zip(ax, iis)]\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we also need to implement it in back projection. Since \n",
    "\n",
    "$$H = \\sum_{\\theta} P(\\theta) A(\\theta) \\otimes \\hat{\\theta}$$\n",
    "\n",
    "it follows that \n",
    "\n",
    "$$H^T = \\sum_{\\theta} A^T(\\theta)P^T(\\theta)  \\otimes \\hat{\\theta}^T$$\n",
    "\n",
    "$A$ is effectively a diagonal matrix, so it doesn't matter. In this case, the important thing is that we apply $P^T$ first to go back to object space, then attenuate correct after."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "obj_bp = torch.zeros([1, 128, 128, 128])\n",
    "for i,angle in enumerate(angles):\n",
    "    obj_bp_i = torch.ones([1, 128, 128, 128]) * image[:,i].unsqueeze(dim=1)\n",
    "    mu_i = rotate_detector_z(mu, angle)\n",
    "    p_i = torch.exp(-rev_cumsum(mu_i * dx))\n",
    "    obj_bp_i = obj_bp_i * p_i\n",
    "    obj_bp_i = rotate_detector_z(obj_bp_i, angle, negative=True)\n",
    "    obj_bp += obj_bp_i"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "torch",
   "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.9.16"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
