zonca · August 5, 2026 16:35
diff --git a/car_output_from_healpix_pysm3.ipynb b/car_output_from_healpix_pysm3.ipynb
 {
 "cells": [
  {
   "cell_type": "markdown",
   "id": "30d80756",
   "metadata": {},
   "source": [
    "# Producing CAR output from a HEALPix PySM3 sky\n",
    "\n",
    "**Question.** *PySM3 only supports HEALPix `Sky()`. Is there a plan to support directly producing `Sky()` on CAR pixels?*\n",
    "\n",
    "**Answer / pattern.** Keep the inputs **and** the `Sky` object HEALPix — that path is mature and is what the presets are built for. Then produce **CAR** at the *smoothing* stage, via [`pysm.apply_smoothing_and_coord_transform`](https://github.com/galsci/pysm/blob/5c9a626/src/pysm3/utils/spherical_harmonics.py) with `return_car=True` + `output_car_resol`. That function goes map → alm → (beam/coord) → `pixell.curvedsky.alm2map` on the requested CAR geometry, so the cost is one SHT, not a `reproject`.\n",
    "\n",
    "This notebook runs end-to-end on the repo's venv (`pixell` is already a dependency) using the tiny `d1` preset (NSIDE 512, 2.3 MB templates shipped with the package). No external cluster or large file download is needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9915c401",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-05T16:35:04.453608Z",
     "iopub.status.busy": "2026-08-05T16:35:04.453274Z",
     "iopub.status.idle": "2026-08-05T16:35:06.864581Z",
     "shell.execute_reply": "2026-08-05T16:35:06.862844Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pysm3.sky.Sky at 0x7f6890ac9940>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import healpy as hp\n",
    "import astropy.units as ua\n",
    "\n",
    "from pixell import enmap, curvedsky, reproject\n",
    "\n",
    "import pysm3\n",
    "import pysm3.units as u\n",
    "from pysm3 import apply_smoothing_and_coord_transform\n",
    "\n",
    "NSIDE = 512\n",
    "FREQ = 353 * ua.GHz\n",
    "# CAR pixel size in arcmin; must divide 180 min cleanly for fejer1, so use ~8 arcmin.\n",
    "CAR_RESOL = 8 * ua.arcmin\n",
    "BEAM_FWHM = 30 * ua.arcmin\n",
    "LMAX = int(1.5 * NSIDE)\n",
    "\n",
    "sky = pysm3.Sky(nside=NSIDE, preset_strings=[\"d1\"])\n",
    "sky"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c14c17e0",
   "metadata": {},
   "source": [
    "## 1. The usual HEALPix path\n",
    "`Sky.get_emission()` returns a `(3, npix)` HEALPix map. Its signature is `(freqs, weights)` only — it does **not** accept `return_car`/`output_car_resol`, so CAR has to come from a later step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a7608e11",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-05T16:35:06.867820Z",
     "iopub.status.busy": "2026-08-05T16:35:06.867558Z",
     "iopub.status.idle": "2026-08-05T16:35:16.332521Z",
     "shell.execute_reply": "2026-08-05T16:35:16.330730Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3, 3145728) | nside: 512\n"
     ]
    }
   ],
   "source": [
    "m_hpx = sky.get_emission(FREQ)\n",
    "print(m_hpx.shape, \"| nside:\", hp.npix2nside(m_hpx.shape[-1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4ed9954",
   "metadata": {},
   "source": [
    "## 2. Produce CAR at the smoothing stage\n",
    "\n",
    "`apply_smoothing_and_coord_transform` is the published utility for going HEALPix → (beam/rotation) → either format. When `return_car=True` it runs the last step as `curvedsky.alm2map` directly into a CAR `enmap`, so we never call `reproject` for the smoothed result.\n",
    "\n",
    "Geometry uses `enmap.fullsky_geometry` with `variant=\"fejer1\"` (matches what PySM uses internally and what the smoothing tests use)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ad4dc3ef",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-05T16:35:16.336023Z",
     "iopub.status.busy": "2026-08-05T16:35:16.335763Z",
     "iopub.status.idle": "2026-08-05T16:35:17.036145Z",
     "shell.execute_reply": "2026-08-05T16:35:17.035027Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ndmap\n",
      "shape (3, 1350, 2700) | wcs: ['RA---CAR', 'DEC--CAR']\n"
     ]
    }
   ],
   "source": [
    "# Sky.get_emission does not accept CAR flags; drive the utility instead.\n",
    "# This goes: HEALPix m_hpx -> map2alm -> (beam) -> curvedsky.alm2map on CAR geometry.\n",
    "# return_healpix=False so we only get the CAR map back (default True gives a tuple).\n",
    "m_car = apply_smoothing_and_coord_transform(\n",
    "    m_hpx,\n",
    "    fwhm=BEAM_FWHM,\n",
    "    lmax=LMAX,\n",
    "    return_car=True,\n",
    "    return_healpix=False,\n",
    "    output_car_resol=CAR_RESOL,\n",
    ")\n",
    "print(type(m_car).__name__)\n",
    "print(\"shape\", m_car.shape, \"| wcs:\", m_car.wcs.wcs.ctype)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59405d2c",
   "metadata": {},
   "source": [
    "### Sanity check: reproject HEALPix output to the same CAR geometry\n",
    "`reproject.enmap_from_healpix` is the standard way to make a CAR map from PySM output today. It should match the alm-driven CAR output from cell 2, to numerical precision (both reconstruct the same band-limited alm)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "27141d02",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-05T16:35:17.039041Z",
     "iopub.status.busy": "2026-08-05T16:35:17.038792Z",
     "iopub.status.idle": "2026-08-05T16:35:17.869399Z",
     "shell.execute_reply": "2026-08-05T16:35:17.868036Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shapes match: True\n",
      "max abs diff: 0.0\n"
     ]
    }
   ],
   "source": [
    "m_car2 = apply_smoothing_and_coord_transform(\n",
    "    m_hpx,\n",
    "    fwhm=BEAM_FWHM,\n",
    "    lmax=LMAX,\n",
    "    return_car=True,\n",
    "    return_healpix=False,\n",
    "    output_car_resol=CAR_RESOL,\n",
    ")\n",
    "assert m_car2.shape == m_car.shape, (m_car2.shape, m_car.shape)\n",
    "print(\"shapes match:\", m_car2.shape == m_car.shape)\n",
    "print(\"max abs diff:\", np.abs(m_car2 - m_car).max())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a618fe4c",
   "metadata": {},
   "source": [
    "## 3. Equivalence vs. reproject from HEALPix output\n",
    "Smoothing in HEALPix (`alm2map` at nside=512) and reprojecting to the *same* CAR geometry should agree with `curvedsky.alm2map` straight into CAR, to numerical precision (both reconstruct the same band-limited alm). Here we verify that on the I map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2c3b9e94",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-05T16:35:17.872963Z",
     "iopub.status.busy": "2026-08-05T16:35:17.872725Z",
     "iopub.status.idle": "2026-08-05T16:35:19.890481Z",
     "shell.execute_reply": "2026-08-05T16:35:19.888935Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "map -> alm\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Projecting\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CAR shape: (3, np.int64(1350), np.int64(2700))\n",
      "mean diff: 0.0012883623479942874 | max diff: 0.4388401202698695\n"
     ]
    }
   ],
   "source": [
    "shape, wcs = enmap.fullsky_geometry(CAR_RESOL.to_value(ua.radian), dims=(3,), variant=\"fejer1\")\n",
    "\n",
    "# Independent reproject of the smoothed HEALPix map.\n",
    "# We can't pass fwhm through Sky.get_emission, so smooth then reproject:\n",
    "m_hpx_smooth = apply_smoothing_and_coord_transform(m_hpx, fwhm=BEAM_FWHM, lmax=LMAX)\n",
    "import warnings\n",
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter(\"ignore\")\n",
    "    reproj_car = reproject.enmap_from_healpix(\n",
    "        m_hpx_smooth, shape, wcs, lmax=LMAX, rot=None, ncomp=3\n",
    "    )\n",
    "\n",
    "# Same alm-driven CAR output (cell 5 called this with the same beam/lmax)\n",
    "alm = hp.map2alm(m_hpx_smooth.value, lmax=LMAX, use_pixel_weights=True)\n",
    "recon_car = curvedsky.alm2map(alm, enmap.empty(shape, wcs))\n",
    "\n",
    "diff = np.abs(reproj_car[0] - recon_car[0])\n",
    "print(\"CAR shape:\", shape)\n",
    "print(\"mean diff:\", diff.mean(), \"| max diff:\", diff.max())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a99ddb15",
   "metadata": {},
   "source": [
    "## 4. Quick visual\n",
    "Toward neutral foreground I-band intensity on the CAR grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a8b54406",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-08-05T16:35:19.893016Z",
     "iopub.status.busy": "2026-08-05T16:35:19.892758Z",
     "iopub.status.idle": "2026-08-05T16:35:20.983491Z",
     "shell.execute_reply": "2026-08-05T16:35:20.981719Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAEhCAYAAAB4GHRqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAdDlJREFUeJztnXecFEX6/z/VM7OzsOwuURBYogiIiJ7oIQKChwoISlBRMKB+URD1ZzgDigH1jvPM8QyHoKIcwXQoBo47VMyCHHiCAhKWKGl3SbsTun5/dFd1dU/PTM8ym2aet6+Wne6nq6urq6s+9VRoxjnnIAiCIAiCyDK0mo4AQRAEQRBETUAiiCAIgiCIrIREEEEQBEEQWQmJIIIgCIIgshISQQRBEARBZCUkggiCIAiCyEpIBBEEQRAEkZWQCCIIgiAIIishEURkDRdccAGOOeaYGrt+ly5dMHz48Bq7fjxqOl2I2sG6devAGMPMmTNrOioEUW2QCCIAAIMGDQJjTG65ubno3LkzpkyZgkOHDqUcXjgcxosvvohTTz0VTZo0QcOGDXHKKafgwQcfxPbt212v3ahRI5SXl8eEtWvXLuTk5IAxhv/7v/+r9D26cfDgQcyZMwejRo1CvXr1wBjDxo0bPZ8/aNAgNGjQIK1xqov885//xHnnnYejjz4awWAQrVu3Rt++ffH000+jpKTE9ZyzzjoLjDH8v//3/+KG65YvO3XqhDvvvBMHDhzwFDfOOV5++WWccsopaNiwIQoLC3HyySfjhRdegK7rnsJYuXIlhgwZgoKCAuTn5+Occ87BDz/84OlcgiBqLySCCEleXh445+CcY/v27bj++uvx5z//GRdddFHKYY0ZMwY33HADLr/8cvzvf//Dli1bcNddd2HGjBkYP358jH1ubi7Kysrw3nvvxRybNWsWcnJyKnVPyXjmmWfw9ttvY8yYMbj66qur5BqZTDQaxaWXXooLLrgAJ554Ij799FOUlZXhu+++wyWXXIKHHnoI999/f8x5GzduxOLFi9G6dWu8/vrrruJXoObLHTt24Pbbb8cjjzyCCy64wFMc77rrLlxzzTUYOXIkfv31V2zcuBEXX3wxJk6ciDvuuCPp+T///DP69u2L/Px8rF27Fr/++iuaN2+OM844A6tXr/YUh7rAMcccA845xo0bV9NRIYjqgxME5/ycc87heXl5MftHjBjBAfBffvnFc1g///wzB8AnTZoUc6ykpIT/5S9/ibl28+bN+VlnncUHDRoUc0737t35FVdcwQHwq6++2nM8nIwaNYp37Ngx7vFbb72VA+AbNmzwHGa8dHOjc+fO/Pzzz/ccdnWRLF0Sce+993IAfO7cua7Ht23bxl944YWY/XfffTfPz8/n33zzDQfAZ82a5Xp+vPS98MILPefLo48+mvfo0SNm/6mnnsqbNWuW9PwLL7yQN2rUiB88eFDuKy8v582bN+fDhw9Pej5BELUX8gQRCenSpQsAoLi4GCNGjEDz5s0RCoVi7EaPHo3GjRujoqICe/fuBQC0atUqxq6wsDBu63vcuHFYtGgRtm3bJvctW7YMq1atwpVXXuk5zrqu48EHH0SbNm1Qr1499OnTp1Z1Xfz3v/9F3759Ub9+fbRp0wYPPPBATLfM8OHDZReQz+fDUUcdhQsuuAA///yzzU6MM/r5558xcOBA1K9fH61atcKDDz4I7vg2crrT5cCBA3jsscfQt29fXHjhha42Rx99NK699lrbvmg0ihkzZuDSSy/Fqaeeiv79++Pvf/97Stfu1KkTAGDLli1JbQsLCyt1DABCoRAWLFiAc845B/Xr15f7g8EghgwZgoULFybtLk71Wf74448YOHAg8vLyMGHCBABAeXk5pk6diuOOOw65ublo27YtJk2ahF27dgEAvv/+ezDGMH/+fLz44ovo0KEDGjRogPPOO0/a/O1vf0PHjh2Rm5uL/v37Y/369bbru40JUsN97bXX0KlTJ+Tm5uLkk0/GZ599lvC+CaIuQCKISIgoqIuKijBp0iT89ttvmD9/vs1mx44deOeddzB27FgEg0F069YNDRs2xIwZM2IK+kSMGDECDRo0wOuvvy73zZgxA506dcLpp5/uOZybb74ZDz30EO655x5s374dzz//PCZPnozffvvNcxhVxe7du3HnnXfi2WefxbZt23Dffffhz3/+M26++Wab3bvvviu7gA4fPoxFixahpKQEgwYNwv79+222e/fuxeTJk/HYY49hx44duPXWW3Hvvfdizpw5Nrt0p8vSpUtx8OBBDB48OKXzFi5ciG3btuG6664DAFx33XVYsmQJ1q5d6zkMka/atGmT1Pahhx7CTz/9hL/85S/Yt28f9u3bh0cffRTLly/Hn/70p4Tnrl27FuXl5ejcuXPMsa5duyIUCiXN46k8yz179uC2227DI488gl9//RV/+MMfEAqFMHDgQDz33HO49957sXXrVnzxxRfo3r07ZsyYYTt/9uzZ2L59O7755ht8++23WL16NS677DK89NJL2LZtG7766iv897//xY4dO3DppZcmTTvB/PnzsWbNGnz22WdYu3YtGjZsiOHDh6OsrMxzGARRK6lRPxRRa3B2O+zbt48/99xznDEmu6h0XeedO3fmffr0sZ374IMPcgD8hx9+kPs+/PBD3qJFC84Y4z169OBXXnklnz59Ot+xY4frtZs3b8455/yaa67hXbt25ZwbXQ6NGzfmf/rTn3g4HPbUHbZ582auaRq/9dZbbfs3bNjA/X5/jXeH+f3+mPD/+Mc/ck3T+MaNGxOev3HjRg6Az58/3xZmTk4OLy4uttl2796dDxgwQP4+knSJxwsvvJCwKysew4YN43379pW/w+EwP/roo/ntt98eY+uWL1944QWuaRofNWqU52u++uqrvH79+hwAB8Bzc3P5yy+/nPS8Tz/9lAPgTzzxRMyxl156iQPgn3zyied4COI9S5/Px9evX2+zffbZZzkAvnDhwrjhfffddxxATHeyiKOzG3b69OkcAP/f//4n961du5YD4DNmzIgJ9+yzz7ad/9///pcD4NOnT/d6ywRRKyFPECE5ePCgdNs3b94cTz75JO644w7p+WGMYeLEiVi6dClWrVoFwOjaeOmll/C73/0OJ554ogxr0KBB2Lx5Mz7++GNccMEFKC0txU033YQOHTrghRdeiBuHcePGYfXq1fjmm2/w3nvvoaSkBJdffrnne1iyZAl0Xcd5551n29+uXTtb/KqKm266yTabqUWLFrbj3bt3R7t27Wz7hg8fDl3XsWTJErlv48aNuOKKK1BUVIRAIADGmDxv3bp1tvN79OiB1q1b2/Ydf/zx+PXXX+Xvmk4XwbZt27Bw4ULpBQIAv9+P8ePH49VXX0U4HI45R82XjRo1woQJEzB69Gi8+eabnq758MMPY9y4cbj//vvx22+/Yffu3fjrX/+KCRMm4KGHHvIUBmOsUseA1J7l8ccfjw4dOtj2LVy4EA0bNvTkcXPaiO7sfv362fZ37doVAGx5JBHnnnuu7Xe3bt2gaZrn8wmitkIiiJCos3AqKirwyy+/YNq0acjLy5M248aNQ/369fG3v/0NAPD++++juLjYdWZVIBDAWWedhSlTpuCtt97Chg0bcMIJJ2DSpEn48ccfXeNw2mmnoXPnzpg5cyZmzJiBgQMHxlTwidizZw8AoHnz5jHH3PZVN4nitXv3bgDGWJs+ffrgxx9/xLx587Bv3z7oui6PO4XC0UcfHRNmQUGBbWp6VaRL27ZtAQCbN2/2fM6MGTMQjUZxySWX2MTiAw88gJ07d2LBggUx5zhnLd544434xz/+gblz5ya93p49e3DPPfdg1KhRuO2229CsWTM0adIEN9xwAy677DJMnToVO3bsiHt+kyZNAAD79u2LOSbSt3HjxnHPT/VZuo2j++2331z3u+HMC/n5+Qn3x1u+IFm4Pp8P9erV83w+QdRWSAQRKVFYWIixY8di1qxZ2L9/P55//nnk5uZizJgxSc9t0qQJrr/+eui6ji+//DKu3RVXXIE333wTixYtSmlAtLgGAOzcuTPmmNu+dPPkk0/KCpubU7qTxUHsE3H/97//ja1bt+KRRx5Br1690KBBAzDGsGHDBtdrJvNEqGGnM1369OmDvLw8fPjhh57sOeeYPn067r//flsaiW3s2LF4+eWXE4bRokULPPXUUzjzzDMxYcIEFBcXJ7T/9ddfEQ6H0b1795hjxx9/PCKRSIw3RkUMBHYb97N69WoEAgHpbXEj1WcZCARi9jVr1sw2WSAR8fKClzxSmXAJoq5DIohImUmTJmH//v247777sGjRIowaNQoNGzaUx5cvXx4zYFOwdetWALDZO7n88stx4MAB5Ofnp7zC8hlnnAFN02I8Cps3b8aKFStSCqsqWLVqFTZt2mTb995770HTNPTv39+2PxgM2n6/9tprlb5uVaRLgwYNcOutt+Lzzz/H22+/7Wqzfft2vPjiiwCAf/3rX9iwYQMGDRrkajt48GB88sknnjxLjz32GA4dOoT77rsvoV2bNm3AGHP1PIouXeHRciMnJwfnnnsuPvnkE9sssIqKCixcuBBDhgyxzRqLx5E8y6FDh2Lfvn2exSZBEN4hEUSkTI8ePdC7d2888cQT4JzHdIWFQiFcddVVGDlyJL788kscPnwYv/32G2bOnIkHH3wQXbt2xdChQ+OG36pVK0SjUezbtw+5ubkpxa1Nmza47rrr8Mwzz+CVV15BSUkJVq1ahWuvvRannXZape43nfz+97/Hddddh5UrV6K0tBQzZszA008/jYkTJ8pxIqeffjoaN26Me+65B5s2bcLu3bvx+OOPy6nOlSHVdPn73/8OxhhmzZqVMNx7770XY8aMwSWXXIL7778fa9euRSgUwvbt2/G3v/0NPXr0kAsK/v3vf0eTJk1wyimnuIZ1zjnnAABeeeWVpPfTo0cPXHTRRXjttdewZs2auHbNmzfHFVdcgfnz5+Pxxx/H7t27sWfPHjzzzDN4/fXXMXbsWBQVFUn766+/PkY0PfTQQwiHw7jyyiuxc+dO7Nq1C+PHj8fBgweTzi5Lx7O8+uqr0bt3b1xxxRWYM2cO9uzZg61bt+KFF17AX//6V8/hEAQRC4kgolJMmjQJANChQ4cYD8app56KRYsWoaCgAFdffTUaN26Mdu3a4dFHH8WkSZOwdOlST63nyvLkk09i8uTJuO++++Q6NQ899BCOOuqoGNuPPvpIjkt57LHHAADt27cHY6xKvqfVtGlTPPTQQ5g4cSJatGiBe+65B3feeSeeeuopadOkSRMsXLgQoVAI3bp1w/HHH48tW7bgueeeO6Jrp5Iu4pMULVu2TBimz+fDG2+8gblz52LZsmXo06cPGjRogJ49e+LNN9/ElClTcP/992P37t149913cdZZZ0HT3Iudpk2bomfPnnjllVc8fc7igQceAABMmTIlod3LL7+M5557Dm+++SY6deqEjh07SvEZz2Op0qVLF3z22WcoLS3FMcccg/bt22P79u1YsmQJunXrlvDcdDzLYDCIf/3rX5gwYQLuuecetGzZEqeffnrK62cRBBEL49yxohpBeOD999/HsGHD8NBDD+Huu++u6egQaWbEiBH47bff8MUXX9R0VAiCIKoMf01HgKibzJkzB36/n74zlIGI6frvvvtuTUeFIAiiSiERRKTMN998g/nz52PcuHGep+4SdQdN01ynhBMEQWQa1B1GeEbM2Kpfvz7OPvtszJw5M+m3lwiCIAiitkIiiCAIgiCIrIRmhxEEQRAEkZWQCCIIgiAIIivJ6IHRuq5j27ZtyM/Pp2XfCYIgiBqBc479+/ejZcuWcdfJShfl5eUIhUKe7XNyclJelDaTyGgRtG3bNttqsARBEARRUxQXF6f0QehUKS8vR/u2DbDjt6jnc1q0aIENGzZkrRDKaBEkvpTcB0PgR+yHCQmCIAiiqokgjKVYKOukqiIUCmHHb1FsWNYWBfnJPU5l+3W0P3kTQqEQiaBMRHSB+RGAn5EIIgiCIGoAcw52dQ3LyGtgbMmI0tzwzBZBBEEQBJFt6ODQkVzheLHJdEgEEQRBEEQGoUNH8k8Qw6NVZkMiiCAIgiAyiCjniHpYB9mLTaZDIoggCIIgMogIdIQ92mU7JIIIgiAIIoOgMUHeIRFEEARBEBkEdYd5h0QQQRAEQWQQurl5sct2SAQRBEEQRAYRBUfUQ1eXF5tMh0QQQRAEQWQQUe5tIURaLJFEEEEQBEFkFNQd5h0SQQRBEASRQehgiCL5Jzp0DzaZDokggiAIgsggwpwhzJMLHC82mQ6JIIIgCILIIKIePUFebDIdEkEEQRAEkUHonEH34OXxYpPpkAgiCIIgiAyCPEHeIRFEEARBEBlEFBqi0DzYESSCCIIgCCKD4B67wzh1h6UugjjnWLRoEZYvX45hw4ahW7duMTbr16/H559/jnA4jFNOOQUnnnhijE0oFML777+PTZs2oVOnThg8eDB8Pl/KNgRBEARBWFB3mHdSEkGLFy/GxIkT0aZNGyxevBitW7eOEUFXXnklvvrqK/Tq1Qs+nw+33HILxo4dixdeeEHalJWVoX///jh8+DD69euHp556Ch06dMCHH36IYDDo2YYgCIIgCDth7kOYJ3cYhDl1iKUkggoLC7Fw4UIcc8wxYMxdQY4ePRrTp0+Hphn9kddccw169eqFiy66CGeeeSYAYNq0adi7dy9WrlyJgoICbNu2DccddxxefPFF3HjjjZ5tCIIgCIKwQ54g7yQfOaXQs2dPHHPMMQltBg0aJAUQAJx66qnIycnB+vXr5b558+bhoosuQkFBAQCgZcuWGDp0KObNm5eSDUEQBEEQdqJc87xlO1WeAv/85z8RCoVw6qmnAjDG+fz666/o3Lmzza5z585YvXq1Zxs3KioqUFZWZtsIgiAIIpvQwTxv2U6ViqANGzbgmmuuwfjx49GjRw8AwMGDB8E5R2Fhoc22YcOGOHDggGcbN6ZNm4bCwkK5FRUVpfmOCIIgCKJ2o5tT5JNtetX7QWo9VZYCW7ZswcCBA9GrVy8899xzcn/9+vUBIMZLU1pairy8PM82bkyePBmlpaVyKy4uTsu9EARBEERdgbrDvFMl6wRt3boVAwYMQLdu3TBv3jwEAgF5LBgMom3btli3bp3tnHXr1uHYY4/1bONGMBikmWMEQRBEVqN79PLo4NUQm9pN2mXgtm3b0L9/f3Tt2hXz589HTk5OjM2IESMwb948HD58GACwZ88eLFiwACNHjkzJhiAIgiAIOyHu87xlO4xz7lkKbtq0CbNnzwZgdD2NHj0aJ554Inr06IHBgwcDAI4//nhs2rQJt912m00A9enTB3369AFgCJrevXsjPz8ff/jDH7BgwQI0aNAAS5YskV1hXmySUVZWhsLCQvTH+fCzQPITCIIgCCLNRHgYS/AeSktL5YznqkDUea8sPwn185MLnEP7o7jqdz9UKl67du1CRUUFWrdu7Xp8z549yMvLQ25ubtwwqtMmHil5gsLhMEpKSlBSUoI77rgD7dq1Q0lJifTWAIYHZ9KkSTh06JC0LSkpQXl5ubRp0qQJli1bhokTJyI3Nxd33XUXPv/8c5u48WJDEARBEIQdL4OivX5fzI3vv/8erVu3dp18tHTpUhx77LFo06YNCgoKcNlll9nq/+q2SUZKnqC6BnmCCIIgiJqmuj1BLy4/GfUaJB/ye/hABNf+bllK8dq/fz9OPvlknHbaaXjttdegSohdu3ahU6dOuPbaa/GnP/0JW7duxRlnnIGhQ4fi2WefrXYbL9DQcIIgCILIIMTAaC9bqkyYMAHDhw/HgAEDYo698cYb0HUdDzzwAPx+P9q2bYtbb70VM2bMkD1G1WnjBRJBBEEQBJFBVNUU+VdeeQWrV6/GQw895Hr822+/xcknn2ybpd23b18cOnQIP/74Y7XbeKFKpsgTBEEQBFEzeF0NWtg41+RzW25mzZo1uOOOO/DZZ5+5zvoGjC6qpk2b2vaJ37t27ap2Gy+QJ4ggCIIgMogQ93veAKCoqMj2tYVp06bFhHnJJZdgwoQJyM/Px5YtW7Bv3z4AxsLI4ksOmqYhEonYzguHwwAAn89X7TZeIE8QQRAEQWQQOmfQuQdPkGlTXFxsGxjttujwvn37MGPGDMyYMQMAcOjQIQBAr169cM899+Daa69F69at8dNPP9nO27lzJwCgVatWAFCtNl4gTxBBEARBZBCpfjusoKDAtrmJoI0bN2LLli1ye/zxxwEYnqBrr70WANCvXz8sW7YMe/fuled99NFHaNasGbp27VrtNl4gEUQQBEEQGYTONc9bOrn44ovRvn17XHrppVi+fDnmzp2LRx99FHfddZfsoqpOGy9QdxhBEARBZBBRMEQ9DIz2YhOPvLy8mG6nYDCIxYsXY/LkybjoootQWFiIhx9+GJMmTaoRGy/QYokEQRAEUYVU92KJU78ZiFwPiyWWH4jgvt//q8rjVZshTxBBEARBZBBRePPyRKs+KrUeEkEEQRAEkUF4He+T7jFBdRESQQRBEASRQUS4D2GefHBwhOvVEJvaDYkggiAIgsggvH4SI9XPZmQiJIIIgiAIIoNIdbHEbIZEEEEQBEFkEGIxRC922Q6JIIIgCILIIMgT5B0SQQRBEASRQejKJzGS2WU7JIIIgiAIIoOIcoaoBy+PF5tMh0QQQRAEQWQQUd2HiJ58inxUpynyJIIIgiAIIoOojm+HZQokggiCIAgig9C5t0HPesZ+OdQ7lRoVFQ6HsXHjRhw8eDCuTSgUwvbt2xGJRKrchiAIgiAIA/HZDC9btpNSCmzbtg1333032rdvj/bt2+Odd95xtXvwwQfRuHFjdOnSBc2aNcPf/va3KrMhCIIgCMJCB/O8ZTspiaDFixejXr16+P777+PavPHGG5g2bRoWLlyI0tJSTJ8+HTfccAMWL16cdhuCIAiCIOyI2WFetmyHcc4r1SvIGMPrr7+OSy+91Lb/9NNPR/v27TFr1iy5r3///mjatCnmz5+fVptklJWVobCwEP1xPvwsUJnbJAiCIIgjIsLDWIL3UFpaioKCgiq7jqjzLl58KXIa5CS1Dx0I4R9/mFXl8arNpLVDUNd1LFu2DL1797bt79OnD7777ru02hAEQRAEEUsUGiI8+UafzUjz7LD9+/ejoqICTZo0se1v2rQpdu3alVYbNyoqKlBRUSF/l5WVHdH9EARBEERdgz6b4Z20ykBNM4JzzuQKh8Pw+XxptXFj2rRpKCwslFtRUdER3A1BEARB1D1odph30poC+fn5KCwsxI4dO2z7d+zYgVatWqXVxo3JkyejtLRUbsXFxem4LYIgCIKoMwhPkJct20m7DDzjjDPw8ccf2/Z99NFHOOOMM9Ju4yQYDKKgoMC2EQRBEEQ2QVPkvZPSmKDy8nKbd2b37t3YuHEj8vPz5fidu+66C3379sUDDzyAYcOGYebMmdi0aRPeffddeV66bAiCIAiCsENjgryTkido2bJl6N+/P/r374+2bdviySefRP/+/fHII49Im9///vf46KOP8Pnnn+Piiy/G+vXr8Z///AedOnVKuw1BEARBEHaoO8w7lV4nqC5A6wQRBEEQNU11rxN0zofXIJCXfJ2g8MEQPh78UlavE0QfUCUIgiCIDCLKGZiHmV+0YjSJIIIgAMBclgKcGxtBEHUWGhPkHRJBBJHtMCoICSKTIBHkHRJBBEFY3h8hiMgbRBB1FhJB3iERRBDZDucwlgth1B1GEBkAiSDvkAgiCALgMIQQCSCCqPNwzsA9CBwvNpkOiSCCIAxIABFERuB1NWhaMZpEEEEQBEFkFFFdA9M9TJH3YJPpkAgiiGyHMfICEUQGQWOCvEMiKNOh2T7EkaJOoad8RBC1HhoT5B0SQZkMrf9CJCNRHnEeIwFEEHUC7tETRCKIRFBmwzl1dRAGybw5Ip+QcCaIOg+Ht2KfagYSQZkPCSDCSTyhk0gAUT4iiDqDDgZWBbPDOOf47rvvsGbNGhx11FE444wzUK9evRi7NWvW4JtvvkFhYSHOOuss5OXl1ahNImhoeKZDLXsiHXnAGQZjlLcIopYixgR52byybt06nHLKKbj55puxePFi3HXXXWjTpg2+/vprm92f//xn9OzZE++99x6mTp2KLl26YO3atTVmkwzGeeY28crKylBYWIj+OB9+Fqjp6FQ/NCiaANIjVpyf1XA7RhCEKxEexhK8h9LSUhQUFFTZdUSd123ObfDVDya1jx6qwP9GP+IpXr/88gsA4Nhjj5X7hg4dirKyMnz22WcAgJUrV+LEE0/Ee++9h2HDhiEajeLss8+GpmlYtGhRtdt4gTxB2YBacVHrPbtI15iwRJ4fxkBrrhFE7aEqPEHHHnusTQABwNFHH41oNCp/z507F61bt8awYcMAAD6fDxMnTsTixYuxa9euarfxAo0JymRooGt24eb5S6eXJpE3yPjmhvtx8hQRRLWS6hT5srIy2/5gMIhg0N2TNGfOHOzatQurVq3CZ599hlmzZslj//vf/3DcccfZ7I877jhwzrF69Wo0a9asWm28QJ6gTEf9ICYJouxCPG9Pjz2FvBFP1MTzFon96uZmQxBEWhCLJXrZAKCoqAiFhYVymzZtWtywV69ejWXLluGrr75C48aNkZubK4+VlZWhYcOGNvtGjRrJY9Vt4wXyBGUT1CLPXJxeIFs3mBeBwd09SW7hyn8ZwPUjiLTjGjSGjSDSgtr2TWYHAMXFxbYxQfG8QABw//33y7/HjRuHkSNHysHI9erVw/79+232QpDUr1+/2m28QJ4ggsgkVE9L0hlcHo7LzcQmrJKUsvFKYjeRRRBE2jBePS9jggz7goIC25ZIBKkMHz4c69atw549ewAAnTp1woYNG2w24vcxxxxT7TZeIBGUbSTqkiDqJpV8lkzTYBNCqrdHDnbmxmFNS55vhOhJJH4SNU+9Nl8JgkhIVQyM/vXXX2P2LV26FI0bN5bdUueddx5Wr16NFStWSJtZs2bhhBNOQJs2bardxgvUHZZteK1kXLtGQEuM1iZSFT9qFxlj4OBgjFmP1JY3rLCNVTSYsqq0BqZp4Fz3np/i2SUaX5ToOEEQcdE5A0vzB1Tnz5+P999/HwMGDEB+fj6+/vprLFy4ENOnT4fP5wMADBgwABdffDGGDh2Ka665BuvWrcP8+fPx8ccfy3Cq08YLVbJO0Jw5c/DWW2/ht99+Q8uWLTF27Fice+65NpsVK1bg8ccfx6ZNm9CpUyfceeedMS4sLzaJyMp1gipTecQZC8I0DVzX3WeZOceIUGVVvbgtXgg4xgQBthlbjAG6bo3ngfGPCCmmKFDPBwzxw4x/AYDr0djnnq7p+OkMjyBqmOpeJ6jDq3fBVz83qX30UDl+veLPnuO1YsUKfPDBB9i3bx/atm2LCy64AEcffbTNhnOOuXPn4quvvkJhYSHGjBmDzp0715hNMtIugh5//HFMmTIFjzzyCLp164avv/4aU6ZMwWuvvYYxY8YAAH788Uf06tULl19+OYYNG4bXXnsNixYtwooVK9C6dWvPNskgEXQE5zGl5e81fKqwqg8XEcSY6ONXp8jDVDmOKewc9kHNjIExDRwc0LndVppogGbmCW6KICdHmgdoej2RgVS7CJp5NzQPIkg/VI5fx/2pyuNVm0m7COrXrx86duyIGTNmyH2DBg1CQUEB5s6dCwAYPXo0tm7diqVLlwIAotEoOnfujGHDhuGJJ57wbJOMrBRBlcXz18TNytQ125g1LlVcVY+rCDJFTIwQ4u6eIafnT3h4uG6dbl4nnnOGaczwDOm6vXvMdZZaJe6L8hKRAVS3CGo/w7sI2nBldougtA+M7tmzJ1asWIFDhw4BAPbu3YvVq1fj1FNPlTaLFy+WqzwCxkqP5557Lv71r3+lZFPnqa7ByV4GQquVlhgEaw8Eyada02DraiFpvpHuH4c9d/fw2brHlEHRTIQlHEfc6h41xxNB84HlBMACfjCfD0zz2QffexUxrnmOIIjKUBUDozOVtIughx9+GH379kXr1q1x0kknoX379rj66qtx6623AgAOHjyIPXv2oGXLlrbzWrZsiU2bNnm2caOiogJlZWW2rVaTzs8ZxPvApdp9lWwhO7PmYzEVqFKJehkZTZVZeognXt1mUXFzTI9taSDR7aV2bZlbzHMX5wsvkHptLoUPwMA002vEGJhPA8vPA8yBkVCOpXyf5AWqW9B7XnvhzPuW5aR9dtg//vEPvPbaa5g6dSq6d++Or7/+Gg8//DB+//vf45xzzkE4HAYQuxhTvXr15DEvNm5MmzYNU6dOTeft1C2SdWklXbPF+M1jPDoOzwJR9ajPUs7KSuRZUZ8Rc+x3vYDL+Y7gxd9OESXipplew4oQENXNsUSwDT/yBNOMU2SXHZCWRRiJqoVEaq0l1cUSs5m0i6BbbrkFN910E/7f//t/AIAzzzwTmzdvxm233YZzzjkH+fn5CAQC2Lt3r+28PXv2oEmTJgDgycaNyZMn45ZbbpG/y8rKUFRUlK5bs/A6+Lg6Zk7ZPDxq7aN6bzzGRw1Lju1Qg3N0q6gL57EE4RKp47pCs/iXx4qMGC+Kc7+bMFLDSp6XjWnxZpeY6RECY+DlFeCRiJxdxmRfWjLhZt0Pd3qsKCvVHmgGaN3Do9Oe3rM0d4dFo1GUlZWhVatWtv0tW7aUgsbn8+GEE07Ad999Z7P55ptvcNJJJ3m2cSMYDMasfJl2UhVA4u90L1LoFh5jxiBV25iOJF1gauXKNKNlL29NhOMMQ6lYZaXsMouM3OVVhNtzZ3GcdYogsS1mCGtf0pLQ7C7TdWOKPQDm8wEM0A8fhh4KGTPFzEtxZ3gJu2rjdPdla96J+xxrGBJAdQrOGbjuYaPusPSKIJ/Ph969e2PGjBnymx67d+/Gm2++ib59+0q7q6++GvPnz8dPP/0EwFh1cvHixbj66qtTsql2KiOAEh13ihG3TT0nZr/SPWGu4cJizkshPuBmJcdt12KasTges4XDzEvHG/9BL1da8LwYoaV57UN/nB4gbv3tuWJTxBMzx/xwHTwUNvKL6AYTwYuxRSpxGwFK2LbdWVrpioUpaxPZ+izqMDQw2jtp7w6bPn06xo4di9atW6Ndu3ZYv349TjvtNDz11FPSZsKECfjxxx9x0kknoV27dti0aRMmT56M4cOHp2RTrXhtmVaVnexWUFuKapcXNzO0WRNVtuBSBRYTXRtOG2M/dw7GteHozgFiW/iZVLjaBhFDOl88n5ts5eTYAy75gJmOOS51bKXzg5s3xhQxHNwYA2QZw54Xj5Dali+qvTuolt0/Ufeg7jDPVMmK0QCwY8cO7Ny5E61atULTpk1dbXbt2oUtW7agXbt2aNSoUaVt4pHWdYLUriXneI0Y2ziVWrLxEWqYzoGwTEPMWA5l1V+rBXkEAsgWBya9CNbaM7C8C1zp9jDHiChLD1f+2rWtAvRKvLyQbDxMsgHPifKYHIsFxevDDC+NLTjlx5GK43jjh45YKLioxtqSF2hMTHbgpXyuJNW9TlDRC/dDq+dhnaDD5SiecH9WrxNUZd8Oa9GiBVq0aJHQplmzZmjWrNkR21QLrl7qBM19p5iJESwutrZ9iuhhDtFju7ZY5I4j3XrWCk8VVsxySAkxJEWTqsVSiUsqbpNaRiKPnhdhU5mC11zU0LaSAez5zRo/zWGbBpvqs5HCKknTMmGYcZ6v62BuZsVR0+QYpLTi2iWnHofyunnoSiZxRNQ2yBPkmez4gKoYi3AkhZWzNcg8VNzq6rlqgerSzWB1O3HYdITr7B5Y4XGAc92928oTcVrgLNExZnkdpBAz/xZdZc4ZPgkr0VT7j6oJZ56xeeY8pHeyrr+UxpYpnjanN44Z/zJuP8cuiuPch/N6zhl/8rkciQfJ2QiA8je3p6navZrIS5bo2bjGIWEkY7WO2ylu16gNAsjNU5Us3TK1W/pIqIp0UMvv6kxmEkGeyQ4RJEhWGHqp4FKZtaIWOFIUMNth64f4HwNTc2aM04gZ67MwBugcHLrtmLsnhlmiLcYLpcTJuR6MNFBEChMDpO3eAc4tASQ9EWaY8ivk0lvk1mWXincCyvWrGGelkcrz9xK2wC3vOUQXE/etiiEhhJ2eQjddmXSqunjGYsFD8awqdXfqha14cyiD9828o3p/4sbNw363ij/RmCpPN2Y1TtKKW6OqMvnZmV9s3eRxRJtbujltiLqN14UQaWB0lomgZKSzglNbvmCmtOGwl/eOSkvId1MwMRknbi+XNA3M5wOPRICoCEW1cVYqikcq7n2qlasooJ3H5JWsGUHq7XIY8w1VYRVTeSSpfFIVmUDlC22v3gOvcbKJR5c0t3kGk1REypgfQzQoSx+o19A5AF2O37K8LMxMaW73tsSLt9JitYWlhBd3PJB7oEo0zczMzbA1zdqvaWAc5tpDSM8z9WTjvI5LvrSJVGVfyt2XHhtVaREjVKlVKV7yp/NdrwG4Dk/rjdKapCSCqhaz4AeUrglurcjMOZete+4QLHLBOVl5KK1pnYPzKKBzezcYc7balYJdrYzN7jfpkZENR812GmdK+Lag7C82k5Wz4fVh6jgos+VvG1wd2yeiBmbfr5o5SSTqkno93C8fNzwvyPjECTiVClo26k3RYHpmAONTFWDMEMGmd86WPuL6HInTQTnBEumKIBJeG3HMOXYmJhix4obqlXLGz0WYc+scppnvhq57jHuKxOtiFL9FnGzRVYUs4JoXxXlpbUiJaMURXlVxrUygKsdpeWnAuO2rbi1EniDPZI8Icn0xRIGfpKWcyjUEbtdyvAlSAAnhI6IoujeYKYYYYLhYlD4O2xgbbgkbWVExWz1kCR1mCSHdeges61iCxhJAzLouM7WW0k3i7DLhog43v0PFxHVFMtj+ULwOyvWFCDSEmO1Ey7sBy16meaLuBVexFGe/M2wwj82mJMItYXzMe7f6Es1AHPt1LtM2tkC219JGMjCjC1Wmj3xAVjw1BiaWWFAzojnonomHaq4a7T5gWYRt5UUGsY4Ut0x0PTZtzHOsrr0jIKE4UI+JRoAjL9peVYcgYoqd0zPmpbs9zjO3HfMy7iieN7cywptIjVTGBarlajUKIcbNosCDXbaTPSIIiCOEUhVAbq4D9ZgZpqs3gFnRkJaO48xZBiuVEhPxdcRfOZ/B8ODEdGXJSlQRNj4NjDu6IWQaKRHxmd92EmvDMIBxl6rKpnPMapCJsJV4QNVlykmq2BJpwZkSN/MfsyLl0C1PgtuYHed4HjeSjjexRBmHukxBHKxBWUm8Jsa9igHlXHhdbF1TXN6zEBc8avq542ZDJY+beVCVVPawzTjIc+z3bHtfGCC65JiumzMCHfetdHtBg10A2bKpI3zxnLkyw9GWF5QbYMr11PtxHc9mv6RTs9ji7vQQaeKZ6470tNvYwnOOw0uWV9y6wrwO9nbev2sXX8KLw9aoqimqymuTznDd8rpXbSnrgWoWo2obOZldlpMdIoiZA1XUVlyqT9/ZWnQZeCzX1VE9AdItop5uFvgxF+FyzARTuwlsho6XyVZJGv/aLRRBI34zs3tF/M0huyEgCn5TNHHdaPVzTTPHImng0SjAZaee6MiwKjeRXlyzd43J+ABM02Q6yIpf1FTCU8A0c01zZVAwuBEnBjCmWd0mqQ6q9YR41oDVZIp5mA6PnnK9JGOwmCL8mBAlPp+RDLpuq+CMZyHu1XThKc8soQfKGWVu/sEU4WXrGnXxSogxSVyIAjW/C3vH/TPnH0p6yLwi4m5+jkMVvIqAsr8tDu8NU/5V63X1X8akkyyma1YtF5jzDVLi40hYYxakuC2HeFF3ec2Dbln4SD05iTyi8VRhIgFRVeIinhesMqhlypGEmbSBJHcoWdyRxwGzmVHdIoiBusO8kR0iCLAqdbW7CIC9QEWCistRMDL7MTlrC8zyrsjTzMpbCAdR+VjuHrHT+Fv1xLhFIabwUisVbhX8sgXiqDQ0c+yFz3z8pqiBED3Cs8IYzC8kGANwGQP8fiNEXTQ1zNdbrdjU68t7EHEwZ5eJ8ilqeaKEVuJGMxxithsTXysXd8B0iC4x45kqxPMQuGHzdLCYY9YzMuJuVvngYIgVwZBdRjzm+TFLbIi0gP10md5+n/FdrkjE6nJimpG3Itz4V7fOMdLI0aIXwhsMkLMHuXxOXPlbplFMnjHyiFzqwJFWqii1bsks9J3vk5INbWJJxFsRZbZuO7Ui0TnMfljELPXvfMTydbLSW37c1fbeqfcsouDy3kGxUaIq8gR3y0O2PGCFn7hCTpJXXXczKUKtWZuAbQamW/km9nsVCInEWDwvlornMXou10nUVVjZ+Hk5z8O5TNPMPAmo5XjC9706EMWzF7ssJztEkFJWuxVyTi0iKi0x7Ttut4P6pyjgGCC9LOI880Vhum5VQPI8pgZghSuvyyAHPKsvmrPlC/O6XLeHzRV7zZz6rDmvASkyGFMKHFkBmtdlkIKJ+X2mgFErYA6m+ZSElU1h47jPZ3iTNM34Oxo1KntRCZvxYj4zPtGo4UkSs+FM0QO/HwiHwcMRy2sGkeTyQZtXNa8dr8C3eams6DI1DWzpIQSE8nyZ/V/pQRPCSXhRzO+vqR8ilcLR/PYbOJdeIFYv1wgnGjV++7jhhTPPEQWtlZ+M31xXB9wzSO+NOUCZqQmkPh9m7WeaBvh9xi5dV8SqBttAAg5n6jnytUMMyGdgVdxWHCG/VB9TofrMc3XdLjxklB33oByTw5vMOFlRY4hZYNRRbwkPp2w8yWuZz8p5TWeeUO8VDsHuLFwSiRUlKKtBYT5JTfF0i3vUo/K+bfck/1ZEpixX1KglERVq/Jyel3jnVJZ4IimlsTmViIMtD8eGwZhZjuk6mM4t3WueJtdv0xjS/JnO5JAI8kx2iCAOgKndNoD8Yrqum56bWLUjx2u4hWlWQtagYId4cEZAM8dTmCJCxsVZSYjzzQrRupjS4rYV1Mr1o0phZBaITMws4txaX0hzDLKW5zB7wWJWPkxtmSuVN/OZ3VGw1pURhQI4AJ/y4uummPH7gYBfppGhr8yKXzMLZM6tGVDiOWimSPD7jeNKfIyButy+33x2VreikQfcuxycv5WSTD5nOCoOXcbVdp75YJimPFOzG1Hes8YA3fLgME1Jf/VDtblBQ1QeLge4bk0jF2mixhNGejMhLsXX3jnAdWbzpEHzwV5DCpHoSBNmPjfHYfnRXMYsISw9g9Y9G0mh29JF1hJm/gFgNA6k98qs1B15zQhKt8S6LvIk5HOy+qbsSsaWVgwwB/xY+UVtXIh3IGaAObPEuvKu2hpVZjlj/K16nbkakhWekmIx+VIKGljx47DSDrAcZpqsda000TQw2a2shCOFmxq213nSShkZTwi5nWMry8Rucf3K1MJKnnWLm7M8BXe/jipynEIwpnHE7McULW80DKL2S5kND26+SJxHU77LI0I3yhhPdllOdoggAEYhATntmwX8YDkB8PIKo5IQL6pz1gtXXy5ReDleQFnQMMtrEdVh9VsoIsLvMz0cUN4rUblqVhmsGcLHKKNUAcSsaxo35LhN6zjTfMb1xH1xmELIjJouXlyH8BHeF02t7KJKOKa9qLCdhZ/43IEqGgKaEZ9AwBBPkagZFx8QidoLFqVlyzQG5OQYxyNRsNygedsMLBwGj0SNwl7j1r1oprcrHDWFhqOis3np1NpdMZGVlfE/OWZKFt6aJcBsLVIlfFHR+/1gwRwjTSJCBBsrPBsCSwhJZohIUZEdOgynV5EJcSGejbg/TQPLyQECfrBw2PLe6LpSgZvPQwhVmdeZlU5ClIFJYSmEAxN5zszjAACfGTERnkgfXUbYuJ75/MT51ntiehY1M/3E/YnK3u8DmAauR8GiSh6ORm2i0BIgynpKNoFjPg/hCVWem62rjFnj4OweLVhprXZ1qc9c12GNllaXr7CEkv1VURoPIi/ZjlnGXNyaks8Mr60mxaS1rIBZLpgNIKaJbm4zXua9c91KFxlPkZi2Bpg9zgmJEULOG3cRKqIRKk24y29H+LakihUs9tmoRv6xxYmZ75LppUdUN99D5z1yZRcTocuxitLGen2UW1TiWhmtdwTQ7DDvZJEIUgQMAB6NgpebhYYoFG2VmCggxbmwKnTR6pJhQjnXOGbUO6IwFuEDjGnG7C1mf7GZzwcWCIBHzcIesESGBrOQ1M1CXAP0qDV4VtVjaoEmKgI5oBbm30Z8eDQKo8CAdU9qC00dixN1VPbCTh6P2n87PE7SLcxgrG8jBF3MdGluiQGNgQUCQDDHGOYSjhr7ATBfPSAkvCQabF4Zn/lbs3tFbIiPvvq02PsWcVXCs+Jv3pNupJ/NS+Z0z4tK3O+XopL7/MZ9RaNSKDBddCvClhZGJWWshsl8PjC/H1zXwQJ+S6AIAa35wHICRiXu8wN+gEUi4OGIvUKQz0u3utYAyyOimfepK2kizhV1k89vevmEF0sDImErnbnxfkmBHPADefUNUVcRsvKPZnaHqRWHuBZjhmjOCRp5QveZ8dOMa0d1MHPclBSiohtRU0QBg5JWar6EIQ44hwYhfkzBqOtGvtE0U7xakeNirJYQEhyWCJFGVlpZz1OIJiv/2Lqn1QH+IiAZZ/OdkO80zPxkpr14vKJCF3k4GDA9g7rV7QwjPbgeBYPlrWBK2Nz5wG04anrZeBGZwwrLGAOpihHANlFEFJe6496cQspNhMlijlnpJa4tjzNZPjKu2YNgDMwfAHwaeDhses4sgWYJUis8edfCc2vmfx41hzkoUbHdQ00IDa/Ci0RQdoggpiktLvHQZauJWa1vJePb16IxBI1ansV81V28bAyyi8jROLH+NfuI5QBhZlbGaiGmxsfnMyq+cNj4u0Ee+MFDRoXifFmd141GZSUhVx62hW8MNGWabrhGNSEgmOGBqQiZlYMGzsxKRVPSizFrYLVIV9Hyh0giIRIAhEJmZWV2d6nnApbI82nWoENw8Pw84FCFUglGwSJRIBCwvBDiWn4x4Fs3osEN4SQqOS68Rab4tBXKohBWC1nh2dJ8hqjRNKPQC5sVi64IaSVsca4xsNtnxNvvs9JMFLrRqHVdv2b8WREyxJEQrrlBICcAFgpZ8Ymam3lN0V3EmxSAlYeBw+XG/QUCxnWlh8OIK4tGLaEbClvxDvjNgdncOA9mfojoxu/cHMMTZ8Yh2igPWnkYkXp++Mqj0A6UG/EU+aBeEHqDetBMLyjAjPBFEgvRAaaknZHvkZsD7vcZAljXwXMD0HMDQFSHFjLiyA5VGF4mxdvJRF5kDAhHAK7b87ZmeugYjHvy+8Aqwoa3MRwxRF0waAhLmc91a2mIaBQIhSyPrhCNZpnCoxGZp2TZwwwPtPQ0+Xz2rj8hIqNR6cUTky1EHITnkDEzL2lmZc6YMU7ObIyw3CBYTo5RwevGGD6RvkZDyxCSXNdNAWyVWUbqmV2UomiRmkcIFjH2jAGM2505pr382yx3pEdK9UzKoQjOa6jXZeZhbreRZboi1ty0m86V39a7zSNhIGqKGb8P3HRUioaZbYC/WRYwqO+0WRZVmN41IcHkadzxL1EbyQoRJGfXiMJGulvFC8Gtl80cNxI7iBGyy8AaV2QWCIDtpVNn7FiOI/NvdcyH6lERLx9gVXJK15Pl9WFARcgMyxAWnHOru0a9Nxmu0mozr2Ff4NAoxJjZIpUDa0WBKromREUbVVzL1sAES8yo/zq9LNzsphGtW7W1LlvppgdI12XFrgd88JmzpniO3/CWiZY65+DlIeNmhTiQY2DMQloRC0bFI/KFIo7FzCy1O0VU2DkB8Nxc41qaZogVjYEHfEalr3PwHL/ZVcPBc3xG5a2bswIjOvR6AUTzAtD9GnzlUTCdI9wgAK4B/kNRaKEoInl+HGiTg+DeKAIHo4jmaAgciID7NUSDPkRzGCL1NWghHYEDUQT2m94Qv2b0gESNeLCKCLSK+ojU90PP9SOS50NOSQRM1xHN9SGwP2xmeQbuZ9DKjXQFgEieH9EAQ7gx0LTLPuw7lIfQ9lzkbwgDYNCDPnAfA4vo8FXoqGicAy3MofuBwP4ofIX1zbF2AIvoiNYLgPs1sHp+w2nn08D9xvkif0RzGPQAwHSGnLIIwIBorg/hPJ9cxFELc2gRIFLPSIfAIR2+w1Hk7Cs3njk4WCgCRHTD2RkMGOJJA7jfBz1oeNB85RFE83IAn5Gmkfo+aGEO/+EoWCgCX9lhQ3AHA9Bzc6DXz4F2KARWHgbjXIo/BHPAfT5DQIbDlkDVNLCIMXbN1iCA8Z4wbgqgQEDmPx41xRfTjPxWUQEeUbyiiieLwWrIMOaTDQYOsxzQuSFwDh0yru3zAcITYs4MtRpGZkMjYoos4dUKm2Wd+p1C9V2W5YLxHjEYx622hFmOqh4sTTPbjlZ5y6NiEL+iXMQ1ZReteJPVxom5XIcya9DWDSXTSwQrT4wVYUIcizAU0SdXQRceRtMLKcoeHo4oXlNFgYl0EV25ehSoxmFBDPDWHVblMan9ZIUIYj4GMaaBixlU4oV2qn1bswf2VoX5TgoHsuoxsLcalDVgTK8EU7tlhMdHfVkB2YozRJDyNqpTpVkUvFxpwUtvhxmGzs37hdlVAGuMhxA7oo5XBiLLwkLTwMRU5FDYaLmKVra5WevVQCk4udnvbhSYtlWPxfgRddyDaFEDRgEbMT+CJlrI3BwsHYkABw/DZ3YRMrNFDgbjGDM/H6F2A0QjgG4WWKYoMU+W40iMgcqAXHOHMSCijDMRA4zFteA3ul9EF18kCp6XY6ShSEtNAzcXoGQRQ5joAZ9R2ecw6AFzhls4alSoOuAL+MCiUWilh8FCEfjDQQTrGxWdFtKh7a+A75BRuWq5fviCfmgVfuhBBi1q5T09YAgTHtbBNYZogwBYXgAVjQLwHY5ApDbXGPSAZogCv7F2DvdpiOT4oYWNeGuhKKAz6AeByC4fUM6ghTh0v4bA/jD8peVGBRoy0p9VBKEHA2A6h1YRhVYeMirAirDhscvNMTw4QT/0oB++QxHDg5SfA93ng68iCi2kg/s0aBEu8yQL69D9PgQORKCVRwAw5OwpRyDgQ6TQGBvmPxwxPEEMhjiNGN3J3G+IUNHtxv0+Y78p9LXDYfAcv5GvdA4eMMd4+TTwYI7RNVIRAgOHLxyxC/lwRPnWmSaFs1xAsn49M12Mrj8e8IOFwtKrxqUHIWR4ffLrG+VEeUg5D0CFrWVl7zJTG1KA5d0MBAxvq8jLmugeZEZej3Cre1A0tgBYM0aZTDPX7wNa7m7rHPM8JhoONm+N9Y7DJ+7D4aphyi0xZkyeME0sT7FaCDNlooFSFoFDrp8lylYRZ1sQZkEuyk2z4WcJGMeyH4DpSYNVdus6eIXwaisCS0QRzPLyqQu+Vhe0TpBnskIEGR5XHVy0oNRuLB77QoquG+YYrCA/byGUvtJysNaOsVoFDJDTm2Uw0SikG1e8QM6BxEocuVq5Iyonclh93GYLTR2nIO7DrMC5Oq1ZN1ZxZrILUBlHIa4tvCOA4U6PKmHrulI4GvHiESO+XBdTQjWzQjBbvlFr3AYXXRYM8prS9Q9jsCYDzI/DmgVgRcgo2JkGnhOQKxbzUBhylha4IliYNSYFIs5KYcVN74xPebZmYWgkrpkO6ownzay0lBlvWgmUFqByT6a3TBNhmUsD+IWoE+75+vWghfKMcU37SsF1Dt9ePxrszDEqxYqQ0ZVjpo0GI6ycgN/szjGfYVSHTzwPDmu/piFHY2CHK4z4BfwA05AjxKoQxz5mVJ4VYWOwvDl4O9fsSs0XlWvAb4imUMjomjXvw78vCNQLGnE+VG7E1+8DyisMEXTA6PbRGOBnmvQk+nODQDDHuF55hRGfYI4hBswxOcG1Zj4W+dv0DPq2iErOGB8nuoqMZwtAY9D2H5T5lomHzI13kGmaUbFxQ4hKb4vomgxHwMNhq/UvxIZuvDM8ao4fYoeN9ysnAPlxWM4ND5HPFGUVIcNe04yKM+Azx2Tp1uSE3KBxbjBgvCuhkDFw3LxvJt4p0UhQ8hUTXa0wRAOPRGUDAZwZ3grO5TvPdd1692CkL9N84OYMRNHVJHv8RVeczGOK6FDLNgDCcyK/Q6h4qHk4It8Tw0tiNc6Y+m6acWc+nxEvPapcwxJ9Uqg4PFQ2T45ZVsuxeZGIOa7LiL8cFC//xy1hIO5B7fJWhaizvFDdSFw3vFzqGMnqRK3mktllOVkhgmQG5WIAqtKqEF3S0sbtdCVHicysdk8JsSS9OoqAEWNyGaxuMjld0mzJiIHQgOGFkS0TWOE6xBsXFb70rMB6oaNKy9WsyKX3yrwHznTbuVKgCCvpveHmQErdJR6x6cY5N1y/0o1sdE1AJL2ykCQXYxWEGAMgF8Q7cNAo7M148LDRRcIiYfBycxxPJAp1ZpXzm1bcnCUl71MUSOaX1601d5g1dsp6fPb0F90UZvrJFqEIV+musNJJwRxILMZ7MI1ZFWR5hTGeQ1TCobCVtkI8qotKOhaPlBWSuHcxJkqIAjEI1+cDAn5DxIhB2ebsNSE0uOnlsnntzAqX1cs1xH4oZHnfAEOsHThorfHjnBmoiRlghmeOR8xu1kjUEjxC8B4ul2OXYsKydWUb6SPfCPXjvwIxDsrWfWv+T32GYpaQKpZNu5h1hBR03ewu8hvePC6Eqa4DqDCFvNkQiUTByw8b8Sk34yGe5eFy45n7NCAcNNKzvML2zsn2mNkI4WLMHLjxPEW5putGXjLH8dnyoc+qkDmM95r5fMat6rrRXaNzYwatpoFz0yMti0xLqFjp6UgfYWf+y8yuJq4r+QXM8la7F7jGv9LrwpTGqUgTa/kMLruzzPDU91o0THzWxAEj30StewPsH4u23OryPZMNPMV7L2fcyfKLyckfUmQpY/CqE6Ybmxe7bCdLRJA1/ZHZSklxHLIytwpVBtH1JW2kuWg9iAKCAzYvE7OfYPY1WYOtxXH1fNNrAasABmB5IxwVs7qQo/Md4+b9WKspc3kvcoac9IwoLVwzTjxijnsIBMxjUVsaOdPDuCZXykgh9BjU5LZae+IeokDYMTBSVH6asc6J6CIUnh3DwxR23CyM1rgRWcV7BiUOzJqGDct7JLs1lJlN9gJLPG3dKIxFQcjMSluYun3CQu0OiEStgaG6bsyIiQqXetR6Djxqhu1Mby6vbWRlcxSnOuZCdPXJU5W0NaeVCy+VECJiKQHZzWlWjvK6jBlpq0eA6EG7sBIeAnNQODe9KFK8iNlf3BBbhsiDJTx03RBUQgSZXh5xTMZR5n+lQpaPiFv/d4pTsYaRKuaksRBYgLrGi7yG2lWs5AVbhubmgGARb8aAiKZchysNHFPMiKhozPFpFMPjxMyJCIhElIHHRv6xLU4qBJuzS0uJM+eGsJFT/4VHlxkTRYTnRwgrq5/c9JRqhgfZNlNKjCmELhuQ1qdX4JrWXH5qBfKYHBwudqtJK9JfzmBUn5koD6NW3rQXzhDlreW9080uwojVuHM2VET5qJQXVoSMfbYlGeS1uHK+bq5sruQRGd/qFUGWKPRgl+VkhQgS35oCrCoNUN85qzKzeS0ZwLjxJXMlNHmQKy8B47oMR+0hkRcSBYR82VVho1TMTLERrUD1ZXNGA0DMZwTUSwsxp4YLNS6idWVvBXGuATxsFbaw0k5Wss4I6SI82As0NbJyDIMZIjcHZ0JZQNJMKzl4m5teClvCOi4imjSyIlMfqCkwI9zlVG6NkQB33Jf9JtXB7tZtiTDF2iTMqkSkJ8S8v4hSIeqKQAGgLrJnf0ZKjhVdlxCfvmAw+ty4cZ5IU9GloyoGzsDBzDWZRNhGkNBUL591XyI7wpwVaK33JvIMDM+OziC7o2QARrrKL8xxo/uRSw8Ul+M9RBeTTTQ4BGCCx2K9czEeOOO9l2saAUp3rmqnqCqRZFGHV0Cmo/CqinQwnoslaJhlbuYFK07WNblupKH1mBgQ1cBZWElfbnZlwxJ0InDxuKWXyxww7fAAAwDXzCUexMd3AeMOpHhVxJWt8WDepZqvpadTuU0rg1o7mfXcbfnKmbd169qq1mHmeEbrXFEWiPcQAONgUDzy8p7U6FjX5zqXSwXwmDxgmipRBLiRDxiM98ucxOE69V9NAuGp5VZ6x2baKsbl9uLaZTlZIYLsGdWq/NX3HQDsg5uFtXjhxctghqG+LlxUGcpZSiY03k+l/1maMisc+eJy6x+u/uG8Iasgslou9nu0l02i0DErLHk5x0vKxa2a3XRqoSMK1ZhUsuJtT1SXeDvGPMk4iDhx2NOZqx4XqxKy4mneiaacoyJmj6hpoTnSm6vpoMTdtfXmFEqQws/qNrXSybgn2ecKWaGrcWUiiuYfSveqbUVyKwaIzY+QacOU/Ggr0MX9KolsVKSqifOl4LbTjaS3WvBMt+IiK0hHeklRpVutdi7WdlKOG//qztuNDU/Gz8oHaha1izHlTCGAZB4yTrbqeKuiNSpYKx04LC9bjCgz42HzMiteldiouOQBbq7MLt9lsV+3/cvN6xtZhUPOcFTFjHjGuhUG15TyStjoVoPDiLO5X/UWA5agF9dWu7TNZwaZWqoQMcOVs1PVk7j9XwbYpuSLe7Z576x8rSQ7wLnsfbcVejELUIlg7fdga9XIeOjW+S4CWw6CF7N+wezPXxF3CTJzlUGLJXonK0RQrLQBzNJSZlzp7bG9qOKltk5xFn+28OSp9itym6clNgxZHakvsBKHRNeSdmrhyZWXTy3klDDdCnHr3u1L/ieOS7woutgzmC55ZUChsBWtTFnguIQl0lBWDAA3K2CmrnGk3IsQBWowYlE0bgtbHEwkgJLedMx5wuMVqw+chTkzK161Be1ma+VT0VKWA8Dls7UqJVt3rvxsSIJ8ZVcTNs0kPFBW/rLuL264ah0gXEmMmQuCKt4/qOLS8fzB7XlFhusUGAz27hszqHDYvKzpzRDnmGli1/hWmPK7c8o7wB0VJgCpAZTOdMDMd7Znzow4cB22MMUx45mrFS8cIpTb7tl4XzisbnnYkb+5MpnBeodk2kvRo9hHzbF8Is3kyuVynqFr+cAZs2ZjmuHLxp+6EjkAWzebGnkZrpV53IYlqD9kEeg4Zm8MmHfMzGerPhybVuH2AJ2z5JT7NfKt+RzEuEThcTaftzVutBqh2WGeqTIRtG3bNvz73/8GYwznnHMOmjZtajseiUTwySefYNOmTejUqRP+8Ic/xHhivNh4QtPsBZJSgNrqWthXZjXcnmbBxGAvXGM8L0rg6jcxRC3CnOcoZ6ulcKJKSrbS4h1T4iYLUSgDgpWuBqa89TZ3tRABzIp7mrBc30qcRdy4si/mRFst5Twok43pipfIfjjG1i4YHdeSBZfLMcDxDBzCLMaeW38KewZH2GaeUbskkya7ZRszeFfmOUde4eo9OOMq/+cIB4qeVvMLV1rQsO5JpJt0rygXUFvStv+71eBx4hLvt7yG226lgnaeq65g7MwPtlecu5yr/q3YMpgDjGPT3xIQidJaSQ8dLvnFcY4uVZhVzjDA/skVK/5SyHJudpWJ6KndfGaXrjwalTYJ30d7IWuz57o1RsxKELcyjdsXX7RaNDYb12up3i4lc9pyXszlnOW1816SYZYrtnSx6pHKNaiOEK9Fdw1ErbahJTdJnWeffRadOnXC/PnzsWTJEvTv3x/Lli2Txw8cOIA+ffrghhtuwNdff41x48ZhyJAhCIfDKdl4xqzUhNvSjrMwskxsjVhb4ZhIiDlqYLfKyNnyEbbJBJD4N2FBzJSgHbaaY9CmM6yY/Sm8IW5hxRrF2e8sLJhVYadQgHCuKw8MjntQ7knd7/Ys46WLHHvguKeUCjkXWyVczmM9GfHD4fbnlehZuuUbr886URhqPOQ+OMJ1VSaOcPX495EuXMNW4heTfsq+VBpftjRwXFudrZj0XEfaxrG35xlV3FvX49wYB8PVdBbn62JcljhkNSa4DF+3n+vlXXfJL/bvmznLytjzRXxcs5B4b0TZIfKRM2/BeczlcmlXA+K68Jhe6UV0h3nZsp20e4L+85//4MYbb8S7776L8847DwCwb98+7Ny5U9r85S9/wdatW7Fy5Uo0atQIxcXF6NatG1566SVMmjTJs41nolGA+aRGsFp/HNbS6zBnU3FzhL9m6Qmmmf3vouCA7D4zTnVpRbjiofDzimslHeeYaCkqg4Q9XMDb9Z0eE7f7iFeB2AZZulWuHsJzHndr5bld3rUFmgRxTmUqRmc49h3K7astWcR/XE5R7hQhnvGaHxJcX8Qh1ihB3OPYVxa3Z+I1L8bLI3HzoIc040meXzrhgPjYMhPvkujGcgpct3JD8ZjI8lE8u3hlinNfsnR3Rtipf8Tp8vrigHJcnXRh2tiLnzjPy/FblPMJn4ttLBJP8f7EjVSyfEgHrkIvjl2Wk3ZP0GOPPYZ+/fpJAQQAjRo1QpcuXeTvuXPnYvTo0WjUqBEAoKioCEOHDsXcuXNTsvGKbCnZMobV+mZgsu+WicyrFqpitL/YzDU6xOJojJkLljGxMisAZ8GTbEt8A+77Vc+Em6vfHgjssxWOgGSix3M4MX84Wm5xwlaPOT1buktrVexPFI7t+oniHKeyPaJ0cBNFcdIgXiXuZZ8T0Yr2Er9k+dWzh0CeYD/PaxyS7U8Wj3jHU3l+okKOH9nKxy9VOxEf05Yr3hPufBeSvrcJ8qFb/OLFN9V3QRVa3GWfjId5SH40VhnPlZLeMJUU0yC+YWgGbP9X2qpxY9bmlgcU7xRT4hh/LGkVYa4TlGxDGqqDuk7aRdBXX32FAQMG4KeffsKLL76Id955B/v27ZPHQ6EQ1q1bZxNFANClSxf89NNPnm3cqKioQFlZmW0TMHPNEia7xQSO1ozMuIqLORoFj0asz0coLVmuW25iK9xKFgZORBiJCtzKXMNrpVPlcFtaxpDM0+JFJKl2XsWo10rfLQ6VrQjgEHPxSCVcrxWt63lx7s/rdRNV/OrveOLdaRsvrMp644zAKnlagjwbL9hU8oTz3uLdv0wj9TzHPi/XiT3gLZ6phh03v8H+vOMOCXBMqBAnO9chShIvOTbRPEd+jFXC4sRV5FmILoLY/KecJ8SoscRHNYsgnsKW5aS1O4xzjr179+Lf//43Zs+ejb59++Knn37C+PHj8d577+H000/HwYMHwTlHw4YNbec2atQI+/fvBwBPNm5MmzYNU6dOjRM5xC1w5YBAuUqu1erg6irBnFtrpZieIGsxeS4bGHEzVqLCLRFudokKf9epyo79ogJyE1nxKnhAcUun8vYk6xNI8f5sp6ZJMFSGeBW3+ttzJc0d/6Z4fSb/l5hEA75jonOE6efVg6DGK14eTeUanuNn/q9SQipNQiFZA4fFq5ChnKuKJFiVtLNbyfVc5bdNTCWJ9xHhFJGOvJ8oeOn1kj9iDYSN68mw3hXG7X+LOiLh/bkJ4HjlG4/zd9VTVVPkP//8c8yfPx9btmxBx44dMWHCBHTo0MFmc/jwYTz99NP4+uuvUVhYiHHjxqF///41ZpOMtHqCGGOoV68eNm/ejOXLl+Pll1/GF198gTPPPFOO46lfvz4AxIiZsrIyecyLjRuTJ09GaWmp3IqLi40D3D6wz8rjRmkhponKhb2kEHK0wFTMryqL/2SBmqgFmy7UAsxr90iiOCXzetjCiBO+l+4BpzcgGVUtYmobiVrmybxT8Z6L2zP2dM1qJOXxFpWw9RqOV49NwoZCJUSV2/NRu2diREu8a6pll9hgHUvkjXGWefE40gZKvLya/MJWuG73L79/GOdcZm3M/LCs7SOvtgjGi4Kx2rYUPok8VzU5JqgKeOCBBzBlyhR06NABY8aMwc6dO9GtWzd8//330oZzjiFDhmDWrFkYMWIEioqKMHDgQLz99ts1YuOFtA+M7ty5M7p27WoTKwMHDsS7774LzjmCwSDatGmD9evX285bv349OnXqBACebNwIBoMIBoPuB0XrQS7lLzKvORjaNljai6fDOFddOZQ7C8Z4hYW4RmULfzf7mH5sh7cpmTOmMlS2MHS29L2cU1dx3mO8vJHQ25Hi9bzsqytUV9xduz4cu+KNAUkaVirX9WqTIP+k4oFTg/PqmUqFZF48L139zCpjAVjja0SZJoWND8Znc+AQWtYwB4gSW/MZv6OK91+9Vpx4S9HkVSxLz5QH03SSSJs57Txy3XXX4d5775W/R40ahfXr1+Oxxx7D7NmzAQALFizAp59+ivXr16N9+/YAgJKSEtx+++0YOXJktdt4Ie1jgkaOHIkVK1Ygqny08Pvvv8cxxxwDscbP+eefj3nz5qGiwvi6dUlJCRYsWIDhw4fLc7zYpAyDOcjZUSlxZU2MFKfCyv+cU0+VsOMWSl69Lx7jYm1wicuRX8L1muq1K9sC9rLPK66DG2sYt+dcU8KkNgqidL0D6aQ2RCed6aIM2K30Nd08UkeKmzhLqcy0vPnM5wPzB2LjG9M9xuT386zhD97KYrHUgC0st0QVnrUaykhVMUXeudYfADRr1gwHDhyQvz/66CP06NFDihLAEktr166tdhsvpF0E3XzzzcjNzcWAAQPw5z//GWPHjsUbb7yBJ554QtpMmTIF4XAYAwYMwH333YczzjgDrVq1wg033JCSTepY4w603CBYIMf4irZzsHRcN65LRSa+rVWdhXhtqjDSUbE7C6AjKWgTdVPUNpyFejpFcbLrqv8SyWFxKjvjoPlvGvJvwusrJPOwxN0Xp5GQSIzIPFlF+dMtPKXrKm68bN//Mr5FxyNh6+vw3BmO4lFiwj1jeIdSF3gMcnkVMa7IHlElajXV4PGwmTgnFAnnQyJWr16NDz/8EOeee67ct3HjRhQVFdnsWrduLY9Vt40X0i6CGjRogC+//BJXXnklysrKcOqpp2L16tU455xzpM1RRx2FH374AWPHjkV5eTluuOEGfPXVV2jQoEFKNqljLQomX2rNZ7hRYaj8hBlWfVESvbhE5cnGijleYV+VZGM6HyHMp1krnNvg1j/pEPIxF06h29npnT2S8GLCh+XxrQqPUMoNKnGvxgdseTRifJQX3BInsrI3PELMp4EFg8bX5GPEHezP0JU4osx18LnX+0g/XqbHy2nyMJagKSwslNu0adMShr9r1y6cf/756N+/P6655hq5PxQKoV69ejZbMTQmFApVu40XquSzGbm5ubjyyisT2hQWFiZd9NCLTUpw8+XQo+Ah3VLwDADTjNFBTHO0lON06MqWhFkgyIaWBuPje1lQybjN6HGOf0l5JhkyO+0SjQNzplUmp0NdhMM+1Vl9PpoPMR8LO5Ixf7brpuncdOWn6sqXahe7W/midn27eZJFeW/uY5oG5vOBhyOGx0i9jvGHcnFHuR9vJl6y+NcUcXSZqx2A4uJiFBQUyN1xx9YC2LNnDwYOHIiWLVvi7bffhqY0DBo2bIi9e/fG2AOQa/5Vp40XquSzGbUbU7T4jBH+XI9aSp4BzO8HCwQS96E7Xcu21oTbWhYZjltBS10usbhNSwaqrxuMqDzcHPfn1iXEGJjfZ84acj+3quKU8XnHTdjIYzCWKXGme4wAMCewaD4gEABkma9cw+kVEksnmLPO5DVsYittd5l2Uh0TVFBQYNviiaA9e/bgD3/4AwoLC/HBBx/EzNY+8cQTsXLlSujK4rQ//PAD/H4/jjvuuGq38ULWiSAmvD+6bs4U00xXqm50lUWigK58SDVeIOp3uAQuhWNGi6FEY6cIIhNhYiyIYxyJWDTVCb0LR44qhByzYMV3zWJnxwovvvk1AM0UqgE/Yj44DMep4pq6WARXcx8qUZufLU9h88i+ffswcOBAFBYW4sMPP0ReXl6MzZgxY7Bv3z688sorAIx1fJ566imcf/75ct2/6rTxQtaJIPERQa5z8EjELLhEQSZWhg7HfvFcxXStJmzZywvW4heFILygaXHbA1mFHA/jcigapXe9qnBtTFr7hDfIKJMdAlWcpnPo5eWIlpTaP6GT7JlxY8C19B7JQdkus4ETxb26G8NVIILuuOMOrFixArquY9SoURg0aBAGDRqECRMmSJsOHTrg5Zdfxk033YTf/e53aN++PTjneO6552rExguMJ5TFdZuysjIUFhaiPxsOPwsksGTWILZ4A55VNA1aMGj0LUcj7l0/SdabILKQdI0RqWZYIGBWBFn8oSHn+BOBl8HHVfWss6WMidfYVH9zgGmGl04McRDfgzQmanHHbDIAbmtBHWncnGgawIGIXoEl/F2Ulpbaxt6kG1Hndb75z/AFc5PaRyvK8fMTd3mK18qVK7Ft27aY/fn5+Tj99NNt+0pKSrBixQoUFhbixBNPlMvj1JRNIqpkYHSdQQ684zC+LG8OahbHVDuBWQjyaNQ26K5OuUqJmqGyA8VrmKwUQPEKUt3Z9RIHdVBvVaHO1KpjeSolEnW7S2+QudRJTgCM5Rjdk9EomM8H5NUDyg6AhyP2c6sjyWrquXj18qQQvRNOOAEnnHCCJ9uGDRuif//+tcYmEdkhgsTXguWgN7gUUB67sMSxSATw+Y3ZBsI1mumFEXHk1LX8kW15OmZsyRF4dKoj3bLp2QjcGqgawHx+a5+mGZNc/H7w3CCw/yDklwLiTWlPGy5dpuoM4mpAnf6ezC7byQ4RBFgvgH0nrJU94b1A4caXiBmLGkKIMXMVURJCRIaSDfnatQsrw+85EzA985xzIAzLaymeZ6ljCr3aA1BlOKbUV/e7UwWeoEwlO0QQ141+Yjf3tHgXKpFJuc7BmPACwcpQ2VBhEJlPHR3DVGlSaASlhWwZ01MdcA4on2qS+5yke5Cy2u3pbFRbF0V1zyyoqq/IZyLZIYIE6vo+zsxaKeHCjRllathUoBGZAnXxEplKVeTpuONIYW8kVwfkCfJMdokggZgJJqnEIDnzfJ4tgxOJ7ITydPohD1DNUFVldLxB8M6VrKsTEkGeyQ4RZFv4kMdX7K7nJbGR4SjnUOFGENmNY0iItZ/KhxqjKtPdKYTULjKg2sWG1w44Wv4rW0QQNztIKztCP27BZc4eqczgaoIg6i7xGkgxa9gof9eUV4CoHlQhFLNmXA0tlujFLsvJDhEkqOwU13gD6YT4YbYdBEFkOl4/GZPJn80hvOE6WLpqoYHR3skuEeRGvBkwTu9PMjc2tfAIgkgEdYVlKTXgAeSQ6/4mtctySAQJpNfH/F8qK0BTwUYQhBvVsXI0Ub2kOrBddIdV52KJ5AnyDIkgQaI1UUjkEETmUV2ztKj8yEzifS6pNoheGhPkGRJBABVSBJEtuH58k95/ooqp5jxGniDvaDUdgWqDCjqCIGpry52oO8TM/ErBvrrgKWxZDnmCCILILuJ9TiHeMYJQSSZ8XIU1eYJqK9khgjgHtCpo7ZErnSAyg9r0Hqvft9LpM9+1jnh5JVF9UN3Zi8YEeSY7RFAyjkTMkBAiCCJdmOKHmTOKOKMFFusEto+o1vwzY7qxebHLdqp0TJCu69ixYwdKSkriHt+3b5/x/a0EYSSz8USy81MdF+BcIp0xWoOcIDKBdH9pPBXMCpTrOrjwAtGYpdpNbXw+NCbIM1Uqgh544AEcffTR+L//+7+YY4888giaNGmCli1bonnz5njllVcqZVNjCFGlqUlYC18GgiBSoxa05A243cNA1D6cz6VW5BuAce55y3aqTAR99tlnmDVrFvr27RtzbO7cubjnnnswZ84cHDp0CE8++STGjx+PTz/9NCWbtFMJbxCDMqCSMhRBEOlElCtUttQ+4n05Xt1qCvIEeaZKRNCePXtw+eWXY+bMmSgoKIg5/swzz2DEiBE4++yzwRjDmDFj0Lt3bzz33HMp2XimqjIk5+B6lAqoykAtW4JIDpUtdZNk352sYsTsMC9btlMlIuiqq67CZZddhj59+sQc03Ud33//fcyxfv364dtvv/Vsk1ZUkZRKpqUWWuWhdCMId2i6ft0hXn3hHDNa3ZAnyDNpnx321FNPYfv27Xjrrbdcjx84cADl5eVo2rSpbX+zZs3w22+/ebZxo6KiAhUVFfJ3WVlZ8gjLDGoObOYw/6XcQRBEDUBlz5GT5Z9EoXWCvJNWEbRmzRrcc889+OCDD7B7924AhjBhjGHHjh1o1qyZOfUTiEQitnMjkQh8Ph8AeLJxY9q0aZg6dWryiMas+KmqdphCiKa+EwRB1Ek4r9p5Kqn2GFQzNEXeO2kVQRs3bkT9+vVx4YUXyn1ievyJJ56I7777DkVFRSgoKMDOnTtt5+7cuRMtW7YEAOTn5ye1cWPy5Mm45ZZb5O+ysjIUFRXZjZwCSL4s5IImCILIGKqqKPeyYnRNQ4sleiatY4IGDRqEHTt22LaBAwdi6NCh2LFjhxQkffv2xaJFi2znfvzxx+jXr5/87cXGSTAYREFBgW0DkPxbL5wDXK8dmZcgCIKoG6iz92rZGFEaFO2NGlkxevLkyejfvz/++te/YtiwYZg5cybWrVuHuXPnpmSTEupo/URK3mlXmUxNAxsJgiAyFHPgaG0u370Kstp8D9VElX9FvlGjRmjUqJFt3+mnn45//vOfWLBgAQYNGoRly5Zh0aJF6NKlS0o2lcKrYj8SVV/LWgQEQRBEmuCo9eU7TZH3TpV7gl5//XXX/YMHD8bgwYMTnuvFptIk69clbw5BEARRF6ExQZ7J3g+ocp5cCAE0S4wgCIKoU9DsMO9krwgCYsWNKopI+BAEQRB1EBJB3sluEeSEhA9BEPEgrzBRV6CB0Z4hEUQQBOEFqjCIOgKtGO0dEkEEQRAEkUnQwGjPkAgiCIIgiAyCPEHeIRFEEARBEJkEjQnyDIkggiAIgsggqtoTpOvGtDJNq/L1lqucun8HBEEQBEFY8BQ2r0Fyjk8++QQjRoxATk4ORo4c6Wq3atUq9O3bFzk5OWjatCluu+02RKPRGrNJBokggiAIgsggWJR73ryyYcMGPProo7jsssswdOhQV5uysjKcffbZOPbYY7F9+3Z88MEHePXVV3HvvffWiI0XSAQRBEEQRCZRBZ6gDh064JNPPsHIkSPh97uPpHnjjTdQWlqKp59+Gk2aNMHvf/973HbbbXj22WcRCoWq3cYLJIIIgiAIIoNg8PgB1TRf98svv0TPnj2Rl5cn95155pkoKyvDjz/+WO02XqCB0QRBEASRSaQ4O6ysrMy2OxgMIhgMpnzZnTt3olmzZrZ9Rx11lDxW3TZeIE8QQRAEQWQQnrxAygyyoqIiFBYWym3atGlpi4uYScYSfLC8Om2ckCeIIAiCIDKJFFeMLi4uRkFBgdxdGS8QABx99NH49ddfbft27doFAGjRokW123iBPEEEQRAEkUEwzj1vAFBQUGDbKiuCTj/9dHz//fc4cOCA3Ld48WIUFhbiuOOOq3YbL5AIIgiCIIgMoiqmyANANBpFJBIB5xycc0QiEdu6PGPGjEGTJk0wYcIEbNu2DZ999hn++te/4qabbkJOTk6123iBusMIgiAIIpOoog+oNm/eHCUlJfJ3bm4uGjZsiN27dwMAGjRogEWLFuHGG2/Esccei8LCQlx33XW2tXuq08YLjPPM/XhIWVkZCgsL0R/nw88CNR0dgiAIIguJ8DCW4D2Ulpbaxt6kG1Hn9Tv9Hvj9ucnjFSnHZ188WOXxqs2QJ4ggCIIgMgj6irx3SAQRBEEQRCZBX5H3DIkggiAIgsggmG5sXuyynbSLoIMHD+KVV17Bp59+inA4jFNOOQU33HADCgsLbXarVq3CE088gU2bNqFTp064/fbb0aFDh5RtCIIgCIJQIE+QZ9I+Rb5v375Yt24dLr74YowbNw4LFixA7969bXP5f/rpJ/Tu3Rt+vx833XQT9u7di169emHr1q0p2RAEQRAEYYfp3POW7aR9dlhZWZltlPmuXbvQvHlzzJkzBxdeeCEA4JJLLsHmzZvxxRdfADDWHujUqROGDx+Oxx9/3LONl7jQ7DCCIAiiJqnu2WEDTp7seXbYf5ZNy+rZYWn3BDkTMi8vDz6fD+Xl5XLfv/71LwwbNkz+9vl8GDp0KBYtWpSSDUEQBEEQDjgA3cNGjqCqHxj96KOPIicnBwMHDgRgjBnavXs3WrVqZbNr1aoVNm/e7NnGjYqKClRUVMjfzi/jEgRBEESmo34SI5ldtlOln81455138MADD+CFF17A0UcfDQAIh8MAjJUmVerVq4dQKOTZxo1p06bZvoRbVFSUtnshCIIgiDoBhzU4OuFW0xGteapMBH3wwQe45JJL8OSTT+Kyyy6T+/Pz8xEIBLBnzx6b/Z49e9C4cWPPNm5MnjwZpaWlcisuLk7jHREEQRBEHcCTAPI4gyzDqZLusIULF2LUqFF45JFHcP3119uO+Xw+dO/eHcuWLbPt//bbb3HSSSd5tnEjGAxW+uu3BEEQBJER6ACYR7ssJ+2eoI8++kgKoBtuuMHV5qqrrsK8efOwZs0aAMBXX32FxYsX46qrrkrJhiAIgiAIO0zXPW/ZTto9QaNHj4bP58OcOXMwZ84cuf+qq66SAmbixIlYuXIlTjzxRHTs2BHr16/HH//4R4wcOVLae7EhCIIgCMIBLZbombSvE/Tll19Cd1GXbdq0QZs2bWz7duzYgeLiYrRv3x5NmzZ1Dc+LTTxonSCCIAiipqnudYL+0PVW+H3Jh4ZEohVYvPqxrF4nKO2eoN69e3u2bdGiBVq0aHHENgRBEARBmNCYIM/QB1QJgiAIIoOgdYK8QyKIIAiCIDIJGhPkGRJBBEEQBJFJ6BxgHgQOfUCVRBBBEARBZBTkCfIMiSCCIAiCyCS4DnhZA4jTyGgSQQRBEASRSegcnj4MRt1hJIIIgiAIIqPgujcvD3mCSAQRBEEQREZBY4I8QyKIIAiCIDIJ6g7zDIkggiAIgsgkyBPkGRJBBEEQBJFJcHgUQVUek1oPiSCCIAiCyCSiUYBHk9vpHmwyHBJBBEEQBJFJUHeYZ0gEEQRBEEQmQSLIMySCCIIgCCKToNlhniERRBAEQRAZBOc6uIeFEL3YZDokggiCIAgik+Dcm5eHusNIBBEEQRBERsE9doeRCCIRRBAEQRAZRTQKMA/T371Mo89wSAQRBEEQRAbBdR2c0ZggL5AIIgiCIIhMgrrDPFOrRVA0GsXixYuxadMmdOrUCWeccQYYYzUdLYIgCIKovegcYCSCvKDVdATicfDgQfTr1w/XXnstlixZgjFjxmDYsGEIh8M1HTWCIAiCqL1wDnDdw0YiqNZ6gh5++GFs2rQJK1euROPGjbFp0yYcf/zx+Pvf/46JEyfWdPQIgiAIolbCdQ7uwRPESQTVXk/QnDlzMHr0aDRu3BgA0LZtWwwdOhRz5syp4ZgRBEEQRC3GkxfI3LKcWukJCofDWLt2Lbp27Wrb37VrVyxevDjueRUVFaioqJC/S0tLAQARhD2NESMIgiCIdBOBMYyjujwvYT0E7qHSE/HKZmqlCDp48CA452jYsKFtf6NGjVBWVhb3vGnTpmHq1Kkx+5diYbqjSBAEQRApsX//fhQWFlZZ+Dk5OWjRogWW7njf8zktWrRATk5OlcWptlMrRVBubi4AI8OolJWVoX79+nHPmzx5Mm655Rb5u6SkBG3btsXmzZurNONlEmVlZSgqKkJxcTEKCgpqOjp1Bkq3ykHpljqUZpWjJtONc479+/ejZcuWVXqd3NxcbNiwAaFQyPM5OTk5ss7NRmqtCCoqKsKGDRts+3/99Vd06tQp7nnBYBDBYDBmf2FhIRUWKVJQUEBpVgko3SoHpVvqUJpVjppKt+pqiOfm5ma1qEmVWjsw+rzzzsO8efPkGJ/S0lIsWLAA5513Xg3HjCAIgiCITKDWiqB77rkHhw8fxsCBA/Hggw9iwIABaN68OW688caajhpBEARBEBlArRVBzZs3xw8//IALLrgAJSUluOaaa/DNN98gPz/fcxjBYBD33XefaxcZ4Q6lWeWgdKsclG6pQ2lWOSjdCDcYp9WSCIIgCILIQmqtJ4ggCIIgCKIqIRFEEARBEERWQiKIIAiCIIispFauE5QOiouLsXz5cjRq1Ai9e/eG35+xt5qUVatW4X//+59tX7169XD++efb9nHOsWzZMhQXF6NLly4xny3xalOX2bp1K5YuXYpu3brh+OOPd7VZs2YNVq9ejdatW6Nnz55gjFWZTV3hm2++wYYNGzB48OCY9VC++OILFBcX2/YdddRROPPMM237IpEIvvrqK+zduxcnnXQS2rRpE3MdLzZ1hTVr1mDt2rVo1aoVTjrpJNfnv2/fPnz55Zfw+/3o06cP8vLyqsymLlBeXo4ffvgB+/btQ7du3dC2bVvb8S1btmDp0qUx5w0fPjxm7Zyff/4Zq1evRsuWLXHKKae4pr8XG6KOwzOQRx99lNerV4+feeaZvF27dvy4447jW7durelo1Rh33303b9asGR89erTcrr32WpvNoUOH+FlnncWPOuoofvbZZ/P8/Hw+fvx4rut6SjZ1lfXr1/MRI0bwoqIinpeXx++77z5Xu4kTJ/L8/Hx+9tln86OOOoqfeeaZ/ODBg1ViUxd45513eI8ePfgxxxzDAfBVq1bF2IwaNYp37NjRlv+mTp1qs9m+fTs//vjjebt27fiZZ57J69Wrxx9++OGUbeoCP/zwAz/11FN5165d+bBhw3hRURHv0aMH37x5s81uwYIFPD8/n5922mm8R48e/KijjuLffvttldjUBV5++WXetm1b3qtXLz548GCel5cXU/7MmzePBwIBW14bPXo0LykpsYV1/fXX8wYNGvCzzjqLN2/enPfv358fOHAgZRui7pNxIui///0vZ4zxd955h3PO+eHDh3nPnj35BRdcULMRq0HuvvtufsYZZyS0ue+++3jLli35jh07OOdGOubk5PB//OMfKdnUVZYvX87feustHg6HeceOHV1F0Lx583hOTg5fsWIF55zznTt38pYtW/IpU6ak3aau8I9//IP/8MMP/Icffkgogpyi28no0aP5ySefzA8fPsw5N8QVY4wvX748JZu6wNKlS20ipLy8nJ966ql82LBhcl9paSlv1KiRLR9eeuml/Nhjj5WVfrps6gpvvPEG/+233+TvVatW8WAwyF999VW5b968ebywsDBhOO+88w73+/0y3+zatYu3bt2a33nnnSnZEJlBxomgO+64g7dr186275VXXuGBQICXlZXVUKxqlrvvvpuffPLJfMGCBXzx4sV89+7dMTYdO3bkt912m23fkCFD+NChQ1OyyQTiiaDzzz+fDx482Lbv9ttvt+W3dNnUNZKJoBEjRvB33nmHL126NOY9PHjwIM/JyeHTp0+37e/QoQP/4x//6NmmLjN16lReVFQkf8+ePZv7fD6+Z88euW/58uUcAP/666/TalOX6d69O7/55pvl73nz5vH8/Hz+8ccf8w8//DDGu8a5kR/POuss27677rqLt27dOiUbIjPIuIHRq1atQvfu3W37unfvjnA4jJ9//rmGYlXzrFu3Ds8++yzuuOMOFBUV4cknn5THDh06hPXr18eMgenevTtWrVrl2SbTWbVqlev9b9y4UX7sN102mcaXX36Jl19+GePHj0f79u3x9ttvy2M///wzQqFQwrzlxaYus3jxYtu9rVq1Ci1btkTjxo3lvu7du4MxJu83XTZ1lS1btuDnn3+OyRPhcBjTpk3DtGnTcMwxx+C6666DruvyeLz3b8uWLSgpKfFsQ2QGGTdauLS0FMccc4xtX5MmTQAgazPvoEGDcMcdd8jVtmfNmoXLL78cv/vd79CvXz+UlZUBgK2gBIx0E2nmxSbTKS0tdb1/cSw/Pz9tNpnEhAkTMHv2bAQCAQDAlClTcNlll+Hkk09G27ZtUVpaCsA9b4kB/V5s6iqPPfYYvv76a3zxxRdyn1se8fv9yM/Pl+9bumzqIqFQCGPHjkXnzp0xduxYub9Lly745ZdfUFRUBAD4/vvv0adPH3Tt2hU33HADgMTvcUlJCRo2bOjJhsgMMs4TFAwGceDAAds+8Ttbv6zbp08fW8V66aWXolOnTnj//fcBQC4j75ZuIs282GQ6XvJWumwyiYEDB0oBBAD33nsvwuEwFi9eDCC789/MmTMxefJkzJo1Cz179pT73fII5xyHDh1KmI8qY1PXiEQiuPjii7F582a8//77ts9gHH/88VIAAUDPnj1x/vnnY8GCBXIfvaOESsaJoI4dO2Lz5s22fZs2bQIAdOjQoSaiVCspKCjArl27AACNGjVCo0aNXNNNpJkXm0wnXt4qKChA06ZN02qTyeTk5CA3N1fmP5F/EuUtLzZ1jddeew3XXHMNXnvtNVx44YW2Yx07dsT27dsRDoflvq1btyISicj7TZdNXSISieCSSy7B8uXL8Z///MfTEglqWQfEf//y8vLQvHlzzzZEZpBxImjIkCH4/vvvsWHDBrlvzpw5OPHEE9GyZcsajFnNsX37dtvvdevWYdWqVTjllFPkviFDhuCtt96SfecHDx7EBx98gHPPPTclm0xmyJAh+PDDD2WLkHOOefPmYciQIWm3yRTKy8tjul0WLVqE/fv3y/zXvHlznHzyyZg3b5602bBhA7777juZt7zY1CVmzZqF8ePH49VXX8XFF18cc3zQoEEoLy/HwoUL5b45c+agQYMG6NevX1pt6grRaBRjxozBd999hyVLlqBdu3YxNs6y7uDBg/jkk09iyrqPPvpIjr8T79/gwYPlOkBebIgMoebGZFcNuq7zgQMH8s6dO/Onn36aT5w4kQcCAb548eKajlqN8bvf/Y5PmjSJv/TSS/yhhx7iLVq04H379pVTjTnnfN26dbxx48Z8xIgR/Pnnn+enn34679SpEy8tLU3Jpq5y4MABPnv2bD579mzevHlzPmrUKD579mxbvikrK+NdunThvXv35s8//zwfOXIkb9iwIf/ll1/SblNXWLNmDZ89ezb/y1/+wgHwRx55hM+ePVvey+7du3mXLl347bffzqdPn87vvPNO3qBBA37ZZZfZwlmyZAnPycnh1157LX/mmWd4ly5d+IABA3g0Gk3Jpi6wcOFC7vP5ZB4T25w5c2x2t9xyC2/cuDF/+OGH+X333ceDwSB/5plnqsSmLjB+/Hju8/n4tGnTbOn25ZdfSpsxY8bwsWPH8ueee44/8cQTvFu3brx9+/a8uLhY2hw4cIAfd9xxvFevXvz555/nF1xwAS8sLOSrV69OyYbIDDLyK/KhUAgvv/wyvvvuOzRq1Ajjxo1Djx49ajpaNUZFRQVef/11fPfdd2jQoAF69eqFCy64IKZFs3nzZrz44ovYsmULOnfujOuuuy5mAKAXm7rIrl275MBJleOOOw733nuv/F1aWornn38ea9asQatWrXDttdfGrFqbLpu6wPvvv49Zs2bF7L/qqqtw9tlnAzAGks6cORM//vgjmjZtioEDB2LgwIEx56xcuRIzZ87E3r170bNnT4wfP9423sOrTW3nnXfewZw5c2L2+3w+vPHGG7Z98+bNw8cffwy/349Ro0bhrLPOijkvXTa1nTvvvBMbN26M2d+7d2/ceOONAAyPzTvvvIP//Oc/4JzjhBNOwOWXXx4zjqesrAzPP/+8XA36mmuuQfv27VO2Ieo+GSmCCIIgCIIgkpFxY4IIgiAIgiC8QCKIIAiCIIishEQQQRAEQRBZCYkggiAIgiCyEhJBBEEQBEFkJSSCCIIgCILISkgEEQRBEASRlZAIIgiCIAgiKyERRBAEQRBEVkIiiCAIgiCIrIREEEEQBEEQWQmJIIIgCIIgspL/Dw/KBH7BO9o0AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x300 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "try:\n",
    "    import matplotlib.pyplot as plt\n",
    "    plt.figure(figsize=(6,3))\n",
    "    plt.imshow(m_car[0], origin='lower', aspect='auto')\n",
    "    plt.title(f'PySM d1 I-band, CAR {CAR_RESOL}')\n",
    "    plt.colorbar(label=str(getattr(m_car, \"unit\", \"\")))\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "except Exception as e:\n",
    "    print('plot skipped:', e)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0da5c30b",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "- `Sky()` stays **HEALPix** — inputs, templates, everything.\n",
    "- Produce **CAR** at the smoothing step with `apply_smoothing_and_coord_transform(..., return_car=True, output_car_resol=...)`, or by passing the same kwargs through `get_emission` of a model that threads them (`InterpolatingComponent`, and — with the assert fixed — `PointSourceCatalog`).\n",
    "- That path does **one SHT** into the CAR geometry, so it's strictly better than HEALPix-smooth-then-`reproject` (two SHTs + interpolation) for the same band-limited output.\n",
    "- The libsharp/MPI distributed-smoothing path still asserts `not return_car`; for distributed CAR you'd want to smooth in alm and shard the CAR array across ranks yourself. That's the real open item."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.13.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
 }
No results found