{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9d805639",
   "metadata": {},
   "source": [
    "# Tutorial 4: Galacticus calling *your* Python functions\n",
    "\n",
    "The tutorials so far called from Python *into* Galacticus. Some Galacticus\n",
    "methods invert the relationship: they take a **function** as an argument —\n",
    "an integrand, a cross-section — and call it, possibly thousands of times,\n",
    "from inside the Fortran numerics. The library interface bridges this with\n",
    "ctypes callbacks: you pass a plain Python callable, and Galacticus's\n",
    "quadrature engines drive it directly.\n",
    "\n",
    "Two showcases:\n",
    "\n",
    "1. volume integrals of arbitrary Python integrands over computational\n",
    "   domains, with exact analytic checks;\n",
    "2. integrating a Python-defined cross-section against a black-body\n",
    "   radiation field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c128e14f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-05T19:39:57.697387Z",
     "iopub.status.busy": "2026-07-05T19:39:57.697226Z",
     "iopub.status.idle": "2026-07-05T19:39:57.889777Z",
     "shell.execute_reply": "2026-07-05T19:39:57.889289Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Galacticus library interface loaded.\n"
     ]
    }
   ],
   "source": [
    "import os, sys\n",
    "\n",
    "# Locate the Galacticus library interface. Two supported layouts:\n",
    "#  * a Galacticus source tree built with\n",
    "#      make GALACTICUS_BUILD_OPTION=lib libgalacticus.so\n",
    "#    (galacticus.py at the tree root, the library under galacticus/lib/);\n",
    "#  * an unpacked binary distribution (the `galacticus/` folder from\n",
    "#    libgalacticus.tar.bz2, with python/ and lib/ inside it).\n",
    "# Set GALACTICUS_LIBRARY_PATH to the directory CONTAINING the `galacticus/`\n",
    "# folder if the auto-detection below does not fit your setup.\n",
    "root = os.environ.get('GALACTICUS_LIBRARY_PATH',\n",
    "                      os.path.abspath(os.path.join(os.getcwd(), os.pardir)))\n",
    "os.chdir(root)   # galacticus.py loads galacticus/lib/libgalacticus.so relative to here\n",
    "for candidate in (root, os.path.join(root, 'galacticus', 'python')):\n",
    "    if os.path.exists(os.path.join(candidate, 'galacticus.py')):\n",
    "        sys.path.insert(0, candidate)\n",
    "        break\n",
    "else:\n",
    "    raise RuntimeError(f\"galacticus.py not found under {root} - build the library \"\n",
    "                       \"(make GALACTICUS_BUILD_OPTION=lib libgalacticus.so) or set \"\n",
    "                       \"GALACTICUS_LIBRARY_PATH\")\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import galacticus\n",
    "print(\"Galacticus library interface loaded.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3e9e8a1a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-05T19:39:57.891947Z",
     "iopub.status.busy": "2026-07-05T19:39:57.891585Z",
     "iopub.status.idle": "2026-07-05T19:39:57.893658Z",
     "shell.execute_reply": "2026-07-05T19:39:57.893299Z"
    }
   },
   "outputs": [],
   "source": [
    "plt.rcParams.update({'figure.figsize': (7.0, 4.5), 'font.size': 11,\n",
    "                     'axes.grid': True, 'grid.alpha': 0.3})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6f8abcf",
   "metadata": {},
   "source": [
    "## Volume integrators\n",
    "\n",
    "`computationalDomainVolumeIntegrator` performs quadrature over a spatial\n",
    "domain in its native coordinate system. Your integrand receives the\n",
    "position as a 3-element numpy array of **Cartesian** coordinates\n",
    "(conversion from the domain's own coordinates happens on the Fortran side)\n",
    "and returns a float.\n",
    "\n",
    "First, exact checks on the unit cube:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "dbd86b9f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-05T19:39:57.894482Z",
     "iopub.status.busy": "2026-07-05T19:39:57.894208Z",
     "iopub.status.idle": "2026-07-05T19:39:59.131316Z",
     "shell.execute_reply": "2026-07-05T19:39:59.130486Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "volume                     = 1.000000   (exact: 1)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "integral of 1              = 1.000000   (exact: 1)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "integral of x              = 0.500000   (exact: 0.5)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "integral of x^2+y^2+z^2    = 1.000000   (exact: 1)\n"
     ]
    }
   ],
   "source": [
    "boundaries = np.array([[0.0, 1.0], [0.0, 1.0], [0.0, 1.0]])   # (axis, lower/upper)\n",
    "integratorCube = galacticus.computationalDomainVolumeIntegratorCartesian3D(boundaries)\n",
    "\n",
    "print(f\"volume                     = {integratorCube.volume():.6f}   (exact: 1)\")\n",
    "print(f\"integral of 1              = {integratorCube.integrate(lambda pos: 1.0):.6f}   (exact: 1)\")\n",
    "print(f\"integral of x              = {integratorCube.integrate(lambda pos: pos[0]):.6f}   (exact: 0.5)\")\n",
    "print(f\"integral of x^2+y^2+z^2    = {integratorCube.integrate(lambda pos: float(np.sum(pos**2))):.6f}   (exact: 1)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc246c9f",
   "metadata": {},
   "source": [
    "Now something with real content: the mass of a\n",
    "[Plummer sphere](https://en.wikipedia.org/wiki/Plummer_model)\n",
    "$$\\rho(r) = \\frac{3M}{4\\pi a^3}\\left(1+\\frac{r^2}{a^2}\\right)^{-5/2}$$\n",
    "inside a spherical shell, where the enclosed-mass profile\n",
    "$M(<\\!R) = M R^3/(R^2+a^2)^{3/2}$ gives us the exact answer. The density\n",
    "is *defined in Python*; the spherical-coordinates quadrature runs in\n",
    "Fortran."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6c4ef7d1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-05T19:39:59.133524Z",
     "iopub.status.busy": "2026-07-05T19:39:59.133433Z",
     "iopub.status.idle": "2026-07-05T19:39:59.782166Z",
     "shell.execute_reply": "2026-07-05T19:39:59.781095Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shell mass, Galacticus quadrature of Python density: 9.776424e+11 Msun\n",
      "shell mass, analytic:                                9.776424e+11 Msun\n",
      "fractional error: 2.22e-16\n"
     ]
    }
   ],
   "source": [
    "massTotal, radiusScale = 1.0e12, 0.05        # Msun, Mpc\n",
    "def densityPlummer(pos):\n",
    "    r2 = float(np.sum(pos**2))\n",
    "    return 3.0*massTotal/(4.0*np.pi*radiusScale**3)*(1.0+r2/radiusScale**2)**(-2.5)\n",
    "\n",
    "def massEnclosed(R):\n",
    "    return massTotal*R**3/(R**2+radiusScale**2)**1.5\n",
    "\n",
    "radiusInner, radiusOuter = 0.01, 0.5         # Mpc\n",
    "integratorShell = galacticus.computationalDomainVolumeIntegratorSpherical([radiusInner, radiusOuter])\n",
    "massNumerical = integratorShell.integrate(densityPlummer)\n",
    "massExact     = massEnclosed(radiusOuter)-massEnclosed(radiusInner)\n",
    "print(f\"shell mass, Galacticus quadrature of Python density: {massNumerical:.6e} Msun\")\n",
    "print(f\"shell mass, analytic:                                {massExact:.6e} Msun\")\n",
    "print(f\"fractional error: {abs(massNumerical/massExact-1.0):.2e}\")\n",
    "assert np.isclose(massNumerical, massExact, rtol=1.0e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8a8b927",
   "metadata": {},
   "source": [
    "## A Python cross-section in a Galacticus radiation field\n",
    "\n",
    "`radiationField.integrateOverCrossSection(wavelengthRange, crossSection,\n",
    "node)` integrates a photon cross-section against the field's spectrum —\n",
    "the operation at the heart of photo-ionization and photo-heating rates.\n",
    "The cross-section is any Python callable $\\sigma(\\lambda)$ (wavelengths\n",
    "in Angstroms).\n",
    "\n",
    "Two exact properties make good tests: a zero cross-section integrates to\n",
    "exactly zero, and the integral is linear in $\\sigma$. (The `node` argument\n",
    "lets fields depend on a galaxy's environment; the black-body field ignores\n",
    "it, so we pass a null handle.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "dc4c10dc",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-05T19:39:59.784362Z",
     "iopub.status.busy": "2026-07-05T19:39:59.784277Z",
     "iopub.status.idle": "2026-07-05T19:39:59.787657Z",
     "shell.execute_reply": "2026-07-05T19:39:59.786908Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sigma = 0     ->  0.0   (exactly zero)\n",
      "sigma doubled ->  ratio = 2.000000000000   (exactly 2)\n"
     ]
    }
   ],
   "source": [
    "import ctypes\n",
    "nodeNull  = ctypes.c_void_p(0)\n",
    "radiation = galacticus.radiationFieldBlackBody(3.0e4)          # 30,000 K\n",
    "uvRange   = np.array([100.0, 912.0])                           # ionizing UV, Angstroms\n",
    "\n",
    "rateZero = radiation.integrateOverCrossSection(uvRange, lambda w: 0.0, nodeNull)\n",
    "rate1    = radiation.integrateOverCrossSection(uvRange, lambda w: 1.0e-18, nodeNull)\n",
    "rate2    = radiation.integrateOverCrossSection(uvRange, lambda w: 2.0e-18, nodeNull)\n",
    "print(f\"sigma = 0     ->  {rateZero}   (exactly zero)\")\n",
    "print(f\"sigma doubled ->  ratio = {rate2/rate1:.12f}   (exactly 2)\")\n",
    "assert rateZero == 0.0 and np.isclose(rate2/rate1, 2.0, rtol=1.0e-9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8b57cfe5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-05T19:39:59.789892Z",
     "iopub.status.busy": "2026-07-05T19:39:59.789812Z",
     "iopub.status.idle": "2026-07-05T19:39:59.964747Z",
     "shell.execute_reply": "2026-07-05T19:39:59.964224Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAG5CAYAAADoNcIGAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAfbVJREFUeJzt3XlclNX+B/DPDMuA7KuAoiypiIjmvovknmuuuSJqomR1LTe8ZVYuZYu/ktTKMk27JWrlbnrVct9FRVJR1AAFZZhB2WfO7w8ukyOLA8wwMHzer5evmmfOnPOdh2GeL2d7JEIIASIiIiKq0aTGDoCIiIiIKo9JHREREZEJYFJHREREZAKY1BERERGZACZ1RERERCaASR0RERGRCWBSR0RERGQCmNQRERERmQAmdUREREQmgEkdYd26dZBIJDh06JCxQ6kwHx8fhISEmGx7JQkLC4NEIjFqDAAQEhICHx8fY4dhdBKJBGFhYQapOyYmBi1atIC1tXWN/10tjSl8D5F+vfvuu5BIJEhMTDR2KDUGkzoTdejQIUgkEq1/NjY2CA4OxgcffICcnJwqiePChQt49913+UtZQStWrMC6deuMHYbBvfvuu/jll1+MHUa1dO3aNbz88stwcHDAypUrsWHDBjRt2tSgbRYlWE/+s7OzQ7t27bBy5UqoVKoK1Xvo0CG8++67yMjI0G/AVGOtW7cOK1asMHYYpkOQSTp48KAAIIYPHy42bNggNmzYID7++GPRpk0bAUD07dtXU/a7774TAMTBgwf1Hoch635STk6OyM3NNWgbxmivYcOGonv37iU+l5eXJ7Kzsw0ew7Pk5uaKnJycStUBQEycOFE/ARmJod7DmjVrBABx9uxZvdddmqLf2+nTp4sNGzaI9evXi6VLl4rGjRsLAGLatGkVqnfhwoUCgLh161apbRr6u4Kql+7du4uGDRuW+Fx+fr7Izs4WarW6aoOqwcyNlEtSFWnRogXGjRunefzaa6+hbdu22LNnD06fPo22bdsaMTr9kclkJt1eSSwsLGBhYWHsMGBpaWnsEPRKqVTC3t7e2GFo3Lt3DwDg7Oxc5W136tRJ6/tj+vTpCAwMxFdffYX33nsP7u7uVR5TTVXdPlc1gbm5OczNmaaUB4dfaxkLCwv07NkTAHDjxg2t54QQWLFiBRo3bgyZTAZfX198+umnJdZz8uRJDBgwAM7OzrCyskJAQADef/995OXlacqEhYVh0qRJAIAePXpohnGenHeUk5ODRYsWISAgAFZWVnB2dsbAgQNx5syZcr2v0ua47dmzBz169IC9vT2sra3RsmVLREdHQwihVa5ofppSqcTMmTPh6ekJmUyGVq1aYe/evc9sr2juR2n/ioZQ1Wo1lixZgpCQEHh6esLS0hL16tXDxIkTcefOHU19iYmJkEgkuH37Ng4fPqxVV9FQdmlz6uLj4zF69GjUrVsXMpkMfn5+eOutt6BUKrXKFQ2xHTx4UOefe0lKmlNXdOzevXsYP348XFxcYG1tjW7dumn9bIumCQDA999/r/U+n3T+/HkMHz4c7u7usLS0hJ+fH+bNm4esrKxi8Zw6dQqhoaGwsbGBk5MTRo4ciTt37pT4GSn6PB46dAghISGwt7dHixYtAACZmZl4++230aFDB7i5ucHS0hI+Pj549dVXkZ6ervP5KYlarcbnn3+umSdnb2+P0NBQ/P7775oyRZ+BhQsXAgB8fX0hkUieOX8xPj4ekZGRCAoKgoODA6ytrdG8eXN8/PHHFR42LeLg4ICOHTtCCIGEhATExMRAIpHgyy+/LLH8oEGDYGVlhQcPHiAkJASLFi3Sei8SiQTvvvuu1mv0/T0E/PP7ee3aNbzzzjto2LAhZDIZmjZtio0bN5brHPz5558YPHgw3NzcIJPJ0KBBA4wZMwYJCQmaMkWftdjYWLz44otwcnKCg4OD5nldf0dzc3Px/vvvIzAwEDY2NrC3t0eTJk0QHh6O7OxsTbkTJ05g4MCB8PLygkwmg6enJ3r06KHzlIarV6/i5Zdfhre3N2QyGdzd3dGpUyd88803WuWEEPj666/Rrl072NjYwMbGBp06dSq1nWedK4lEgsOHD+P27dtav/tF8ypLm1OXlJSEKVOmoF69erC0tET9+vXxyiuvICUlRatc0ffLunXrsGHDBgQHB8PKygr16tVDVFRUpX8fqiOmwLXQtWvXAABubm5ax6OioqBUKjFp0iTY2tpi/fr1ePPNN+Hl5YXRo0dryu3ZsweDBg2Cvb09ZsyYAQ8PD+zatQvvvPMOjh07hp07d0IqlWLatGmQyWT46quvEBUVpZkH5O/vDwBQqVTo378/Dh48iP79++PVV1/FvXv3sGrVKnTp0gW7d+9Gjx49Kvw+165di6lTp6JBgwaYPXs2bG1tERMTg1dffRUXL17EV199Vew1ffr0gaOjI+bPn4+srCysWLECgwYNwvXr19GgQYNS23rppZfw3HPPFTv+/vvv49q1a6hbty4AIC8vDx9++CFeeuklvPjii3BwcEBsbCy+/fZbHDhwALGxsXB2doabmxs2bNiAf/3rX3B1dcWCBQs0dT79c3vShQsX0K1bNxQUFGDGjBnw8/PDkSNH8Mknn+DAgQM4evQo6tSpo/UaXX/u5fX48WN07doVrVu3xvvvv4/79+/js88+Q79+/XDz5k3Y2dmhadOm2LBhA8aPH4+uXbvilVdeKVbPnj17MGTIEHh7e2PmzJmoW7cuLl68iE8//RRHjx7FwYMHNX/Nnz59GiEhIbC0tMTrr7+O+vXr48CBA+jevTseP35cYpxnzpxBTEwMwsPDMWbMGGRmZgIovHB89dVXeOmllzBq1ChYWVnh1KlTWLNmDY4cOYLTp09XuKc0LCwMGzZsQOfOnbFkyRI8evQI33zzDfr06YP169dj3Lhxms/A1q1bsW3bNnz22WdwdXWFra1tmXUfOnQIBw8exIABA+Dr64ucnBzs2rULs2fPxs2bN0tNwHQhhMD169cBFH4O27RpAw8PD3zzzTeYMWOGVtmkpCTs2rULI0eO1HyGnZ2dtd4LAAQHB2u9Tt/fQ0+aOHEiJBIJXnvtNUilUnz55ZcYN24c/P390aFDh2e+/2+++QbTpk2Dm5sbpkyZAl9fX9y7dw979uzB5cuXNd9tAHD37l10794dQ4cOxdKlSzU9ruX5HX311VfxzTffYOzYsXjttdcAALdu3cKOHTvw+PFjWFtb49q1a3jhhRfg7u6OGTNmwMvLCw8ePMDZs2dx/PhxDBkypMz39PDhQ/To0QNqtRrTpk2Dr68v5HI5Ll26hMOHD2PKlCmaspMmTcL69esxePBgjB07FgCwdetWDB06FKtWrUJERES5ztWGDRuwePFiPHjwAJ999pnmtWXNGU1KSkLbtm2RmpqKKVOmoEWLFrh48SK+/vprzQhU0fdtkTVr1mgSQTc3N2zduhVLly6Fvb095s2bV+b5qXGMOfZLhlM0p27+/PkiLS1NpKWliStXroi5c+cKAMLX11czD6poLktwcLDW3KhHjx4JFxcX0bFjR82xgoIC4ePjI6ytrcX169e12pw0aZIAIDZs2KA5VtY8mbVr1woAYurUqVrH//rrLyGTyUSjRo2ESqXS6f0+PfcsIyND2NraCk9PT5GWlqY5np+fL3r16iUAiD///FNzfOLEiQKAeOWVV7TqPX78uOY8ltVeSRYtWiQAiDlz5miOqdVq8fjx42Jlf//9dwFAfPTRRzq3UxTzk7p27SokEok4cuRIibG8//77mmPl+bmXpaQ5Md27dxcAxJIlS7SO//jjjwKAWLNmjdZxlDIfLTs7W3h4eIh27doVm7cXExMjAIh169ZpjnXq1EmYmZmJ2NhYrbIzZ84UAIqdSwACgNi9e3extnNzc0VeXl6x419//bUAIH7++Wed3sPTDhw4IACIfv36iYKCAs3x1NRU4e7uLhwdHUVmZqbmeFnz0Ery6NGjEo+PGTNGmJmZiZSUlGfWUfTZWLVqlUhLSxOpqani/Pnzmt/xTp06acpGRUUJAOLMmTNadbz//vvFfvd1mVNniO+honb79eun9Z1y584dYWFhIV5++eVnnpO///5byGQy4evrq/WdUuTJehs2bKg5f08rz++ok5OT1vznkvzf//2fACBOnDjxzPdQkl9//VUAEP/5z3/KLPfLL78IAOLTTz8t9tzAgQOFvb29UCqVQojynauy5tSV9HkZP368ACA2btyoVfb7778XAMTkyZM1x4qugx4eHiI9PV2r/aZNmwpPT88y33NNxKTORBV9mEv616NHD60vwqIv06+++qpYPQMGDBAuLi6ax6dOnSox+RFCiNu3bwsAYsiQIcXqLimp69+/vwAgkpOTiz1X9MV84cIFnd7v08nPzz//XGJSIYQQhw8fFgDEG2+8oTlWlCBdu3atWHlbW1sxbNiwMtt7WtH7Hj16dKmTfFUqlZDL5Zqk28HBoVztPJ3UpaamCgCid+/exco+fvxY2NjYiJYtWxaLUZefe1lKS+qkUmmxhRwPHjwQAMSbb76pdby0hGj79u0CgPjyyy8156noX2pqqqhTp47mgnz//n0BQAwYMKBYPcnJyaUmdS1atHjme8zPz9f8rBISEsr1Hp42Y8YMAUAcO3as2HNFF/ZffvlFc6y8Sd2TcnJyxMOHD0VaWppYv369ACC2b9/+zNcVfTae/mdmZiaGDh0q7t+/rymbmJgopFKp1uIJtVotfHx8ROPGjbXq1SWpM8T3UFG7+/btK1Y+KChItG7duuwTIoT47LPPBADx3XffPbNsw4YNhbOzs1bSLkT5f0d9fX2Ft7d3md+DRcnMvHnzRFZW1jNje1rR9+GECROEXC4vtdywYcOEtbW1SEpKKva7WPQH+t69e4UQ5TtX5UnqVCqVsLOzE02aNClWVq1WC39/f+Hk5KT5zi26DkZFRRUr/+qrrwoAWn9AmQIOv5q4sLAwjB07FhKJBNbW1mjUqFGpw3d+fn7Fjrm4uODhw4eaxzdv3gQANG/evFjZBg0awN7eXmtuSVlu3rwJFxcXeHp6FnuuqP6EhAS0aNEC2dnZUCgUWmVsbW1LHYoqK84n636aLufgWfbu3YupU6eie/fumnlrT/rll1/w0Ucf4ezZs8Xm/lRmrlZZ77lOnTrw9/c32HsuiZeXF6ysrIrVC0Dnuq9evQoAmDFjRrHhvSL3798H8M/7DwgIKFbG09NTa07Tkxo3blxq+19//TW+/PJLXL58GQUFBVrPVfRnVdHPpq6ysrLw/vvv4z//+U+JWwmVJ+558+bhhRdegEQiga2tLZo0aQJHR0etMg0bNkTfvn2xadMmfPLJJ7CxscHvv/+OxMRELF++vNzxG/J7qLS6b9++/cy4iqattGrV6pllgcJpJmZmZlrHyvs7+n//938YP348WrZsiQYNGqBr167o06cPRowYofndGj16NH788UcsW7YMn332Gdq1a4du3bph9OjRCAoKAlA41SUtLU2rPUtLSzg7O6Nbt24IDw/Ht99+i02bNqFVq1bo0qULhg8fjo4dO2rKX716FdnZ2ahXr16p77nod7G850pXaWlpyMzM1LyvJ0kkEjRr1gy//fYb5HK51sKi0n7uQOF30bOmNNQkTOpMnL+/v2ZhxLM8/QVUEYbaDPenn37SLLoosnDhwmKTrCurtHMgnlpYUZqiCf2NGzfGL7/8UmyV7K+//oqhQ4eiTZs2+PTTT9GgQQNYW1sDKPxyVqvVlXsDFaCPn3t569X1fBadj8WLF6Ndu3YllnFycip/cE94eo5hkf/7v//DG2+8gZ49e+LLL7/UTEIvKChAv379jPKz0sXYsWPx66+/YsqUKejWrRtcXV1hbm6Os2fPYt68eeWKu1mzZjp9f0RERGDXrl346aefEB4ejq+//hqWlpYV2ozZkN9Dlf39Lo/SPlflMXDgQCQmJmLv3r04dOgQDh06hI0bN2LRokU4fvy4ZgHP7t27ce7cOezduxdHjhzBZ599hiVLlmD58uV48803cffuXfj6+mrV3b17d82ChLVr12L27NnYvXs3jhw5gm+//RaffvopZs6cic8//xxA4e+ig4MDYmJiSo23WbNmlX7PhqCP76KagkkdlUvRROArV64Ue+7u3btQKBRak4XLSvL8/f0RHx+P+/fvF5vYevnyZa32+vTpo7UyECj5r6+S4hwwYECZdevL7du38eKLL8LOzg67du0q1qMBFK7wtLKywuHDh7W+9B8/fgy5XF6sfHmS5KLzUdLPJjs7Gzdv3ixxMUd1VtSLZmVl9czkouj9x8fHF3suJSWlWE/vs3z//ffw8fHB3r17tSbcF/UeVtSTn8327dtrPVfZz6ZCocCvv/6KcePGFVsIVLTAwRD69+8Pb29vfP311xgwYAB+/fVXDBs2TLMYooi+/ugr7/eQPhR9Fs+fP19scYeuKvI76ujoiFGjRmHUqFEAgNWrV2P69OmIjo7W+qO2VatWmp4xuVyOTp06ISoqCjNnzoSHh0ex78+n/xgKCAhAQEAA/vWvfyE7Oxv9+/fHF198gVmzZsHHxweNGzdGfHw8nn/+eU0vV2nKc67K85lwc3ODnZ1diedPCIErV67Aycmp0n/o1WTc0oTK5fnnn4ePjw82bNhQbMjivffeAwAMGzZMc6yoW7ukIZ+XXnoJQOEK0SfduHEDmzZtQqNGjTRfCJ6enujZs6fWv7KSul69esHW1hYrV67USpZUKhUWL15cLM7Kksvl6NevHzIzM7Fz5040bNiwxHJmZmaQSCTFekvef//9EntQbG1tdR4uc3NzQ9euXbF3716cOnVK67lPPvkEjx490ut71qfS3mefPn1Qt25dLF++XLN68EkFBQWa17m7u6Njx47YvXs3Ll26pFVu2bJl5Y6p6K/7J38uQgjN57yiij73S5Ys0ar7wYMHiI6OhqOjI1544YUK1V2UfD7d+5CZmVmubWrKy8zMDFOmTMGJEycwe/Zs5Ofnl7iSuazvg/Io7/eQPowYMQIymQzvv/9+ifHr0gNant9RlUpV4h96rVu3BvDPFIYHDx4UK+Pk5AQ/Pz/k5eUhMzNT80fRk/+K6klPTy8Wu7W1NQIDA7XamTBhAgBgzpw5JfZuFQ29AuU7V7a2tpDL5Tr1mEmlUgwZMgTx8fHFegw3btyIhIQEvPTSS9Xi9onGwp46KhczMzOsWrUKgwYNQtu2bREREQF3d3fs3r0bu3btQp8+fTBmzBhN+bZt20IqlWLx4sWQy+WwsbGBr68v2rdvjwkTJuCHH35AdHQ07ty5gz59+mi2NBFCYM2aNRX+5XRwcMCKFSswdepUtGnTBuHh4bCxsUFMTAyOHj2KqVOnokuXLvo6LZg+fTquXr2KCRMm4MqVK8X+kuzUqRP8/PwwYsQIxMTEoHv37ggLC4MQAnv37kVcXFyxXg0A6NChA9auXYu3334bTZs2hVQqxcCBA2FjY1NiHJ9//jm6deuG0NBQTJ8+XbNdwqZNm9CiRQvMmjVLb+9Znzp06ID9+/fjww8/RIMGDSCRSDB69GjUqVMHGzZswODBg9G0aVNMmjQJAQEByMzMREJCArZu3Yply5Zphvk+++wzhISEoGvXrpgxY4ZmS5Pz58/D1dW1XJ+nESNGYO7cuejTpw+GDx+OrKwsbNu2rdg8yPIKDQ3F+PHjsWHDBvTo0QNDhw7VbGmSmpqK9evXV3iOj52dHfr27YuNGzdCJpOhffv2SElJwdq1a4v1huvblClT8P7772P9+vVo1KhRidsRFW0bMnfuXIwdOxZWVlYICgoqcY5UWcr7PaQP9erVw+eff46IiAg0a9YMkyZNgq+vL1JTU7Fnzx689dZbGDx48DPr0fV3NDMzE56enhg4cCBatmwJT09PJCcn4+uvv4a5ublmS5EPPvgAe/bs0WxhY25ujsOHD2PXrl0YMGDAM3vV1q9fj08//RRDhgyBv78/6tSpg7Nnz+Kbb75BixYt0LJlSwCFSfLUqVPx9ddf4+LFixgyZAg8PDyQnJyMs2fPYteuXcjPzy/3uerQoQN27NiBV199FZ06dYKZmRlCQ0NL3dh6yZIl2L9/P15++WUcPHgQzZs312xp4u3trfmjvdYy0gINMrCiVT9PLo8vTVkrVEvaNkOIwq0++vfvLxwdHYWlpaVo3LixeO+990q8dda6detE06ZNhYWFRbEVgtnZ2WLhwoWicePGwtLSUjg6OooBAwaIU6dOlev9lrZKdNeuXaJ79+7C1tZWyGQyERwcLL744otiK1JLe5+l1f30saItPEr79+QqsLVr14qgoCBhZWUl3NzcxJgxY8Tdu3dLbOf+/fvipZdeEk5OTkIikWitBCst5ri4ODFy5Ejh6uoqLCwsRMOGDcWsWbNERkaGVrmK/NxLUtrq19JWtD39GRBCiGvXrolevXoJOzs7zTl70tWrV8XEiRNF/fr1hYWFhXB1dRWtW7cW8+fPF3fu3NEqe/z4cRESEiKsra2Fo6OjGDFihLhz545wdnYW/fr1e2YsRVQqlfjwww9Fo0aNhEwmE15eXmL69OkiPT29xNeVVVdJda9YsUI0b95cyGQyYWtrK3r06KFZPfik8q5+ffjwoZg2bZqoV6+ekMlkokmTJuKjjz4S+/fv13lFYtFn48ltQXQxdOjQErfmedKHH34ofH19hbm5uQAgFi5cqNWmIb6HyjqHZX1WS3LgwAHRt29f4eTkJCwtLUWDBg3E2LFjRUJCgqbMs1bH6/I7mpubK+bPny/at28vXF1dhaWlpahfv74YPny4OHnypKbcwYMHxahRozRbvNjb24vg4GDx4Ycf6rQa9vz58yIsLEw0atRI2NraChsbGxEQECAWLFigtQ1IkU2bNomQkBDh4OAgLC0thbe3t+jXr1+J27focq4eP34swsPDhbu7u5BKpVqfgdJ+bnfv3hWTJ08Wnp6ewtzcXHh5eYmpU6cW20mh6DpY0me+MqvKqzOJECY2S5CI6ClpaWlwd3dHREQEVq1aZexwTNbo0aOxbds2/P3332Vukk1EhsE5dURkUp68fVKRDz74AEDhHD0yjKSkJGzduhXDhg1jQkdkJJxTR0Qmo6CgAN7e3nj55ZcRGBiIx48fY9++ffj999/Ro0cPDBo0yNghmpyTJ08iPj4eq1atglqtNr3bLhHVIBx+JSKTIYTA1KlT8ccffyA5ORkFBQXw8fHBiBEjEBUVpdkTkPQnLCwM69evR8OGDfHee+9h/Pjxxg6JqNZiUkdERERkAjinjoiIiMgEMKkjIiIiMgFcKFFOarUaycnJsLOzq9W7VhMREVHVEEIgMzMTXl5eWrctfFqtTOq++eYbfP7558jPz8fzzz+PtWvX6jyBOjk5Gd7e3gaOkIiIiEjb3bt3Ub9+/VKfr3ULJeLi4jBgwACcO3cOjo6OiIyMhJeXFxYsWKDT6xUKBRwdHXH37l3Y29sbJEa1Wo20tDS4ubmVmZETERGRcVXFNVupVMLb2xsZGRlwcHAotVyN6am7ceMGPv74Y5w4cQKXL19GQEAALl++XKxcfHw8Zs6ciWPHjsHOzg4TJkzABx98AEtLSwDAlStX0LZtWzg6OgIo3Iw0KipK56SuaMjV3t7eoEldTk4O7O3tmdQRERFVY1V5zX7WtK8akzFcuXIFO3fuxHPPPYfAwMASy8jlcoSGhiIvLw9bt27FkiVL8NVXX2ndxDw4OBhHjx5FUlISVCoVYmJicOfOnap6G0REREQGUWN66gYOHIjBgwcDKNzs8syZM8XKrF69GkqlEtu2bYOzszOAwh3mZ8yYgaioKHh5eaFJkyZYtmwZBg8eDAsLC4SGhsLcvMacBiIiIqIS1ZieOl26NHfv3o2ePXtqEjoAGDlyJNRqNfbt26c5Nm7cOJw5cwbHjx9Hy5YtERAQYJCYiYiIiKqKSXVRxcfHIzw8XOuYo6MjPD09ER8frzl2//591K1bFwqFAh999BFmz55dap25ubnIzc3VPFYqlQAKx9DVarWe3wE0dQshDFY/ERER6UdVXLN1rdukkjq5XK5ZAPEkJycnpKenax6PGzcOycnJyM3NxbRp0zBy5MhS61y6dCkWLVpU7HhaWhpycnL0EvfT1Go1FAoFhBBcKEFERFSNVcU1OzMzU6dyJpXU6er333/Xuez8+fO1FloULSt2c3Mz6OpXiUTCLU2IiIiquaq4ZltZWelUzqSSOicnJygUimLH5XK51jy78pDJZJDJZIiOjkZ0dDRUKhWAwjl+hky4JBKJwdsgIiKiyjP0NVvXek0qYwgICNCaOwcUbhackpJS6cUQkZGRiIuLw+nTpytVDxEREZkGlVrgxM2H2BefjhM3H0KlNu79HEyqp65fv35YsmQJMjIyNHPrNm/eDKlUit69e1eq7qd76oiIiKj22nM5BYu2xyFFUTS//hY8HaywcGAg+gZ5GiWmGnObsKysLOzatQtAYYKVkJCATz/9FADQvXt3uLm5QS6Xo1mzZmjcuDGioqKQlJSEWbNmYezYsVi5cqVe4lAqlXBwcIBCoTDonLrU1FS4u7tz+JWIiKia2XM5BdN/OIenE6ii+z2sGtdKr4mdrrlHjUnqEhMT4evrW+JzBw8eREhICADg6tWrxW4TtnjxYs1twirqyZ66a9euMakjIiKqhVRqgS4f/veJHjptEgAeDlY4MjcUZtKyb+ulK5NL6qoL9tQRERHVXscTHuLlr088s9yPUzugo7+LXtrUNfcwqTl1RERERPomhEBC2iMcS3iImLN/6/Sa1EzD7GVbFiZ1OuJCCSIiotpBCIG76dk4lvAAx28+xLGEh0jLzH32C5/gbqfb3nL6xOHXcuLwKxERkem5p8jB8ZsPcOxGYRKXlJGt9bzMXIo2Pk5o7+uM74/dRvrjvGILJQDjzqljTx0RERHVOg8f5eLEzfTC3riEh7j54LHW8+ZSCVp6O6KTvws6+rvi+QaOsLIwAwA0rmuH6T+cgwTQSuyKUriFAwP1ltCVB5M6IiIiMnmK7HycuvVPEhd/T/t+qhIJ0LyeAzr6u6CTvyvaNHSCjazkNKlvkCdWjWv11D51hT10xtynjkmdjjinjoiIqObIyivAmUQ5jiU8xPGEB7iUpMDTN3wI8LBDBz8XdPJ3QXtfFzjUsdC5/r5BnugV6IGTNx/gxt9peK6+G9r7uRqlh64I59SVE+fUERERVT+5BSqcv5OhSeIu3M1Avko7xfF1tflfT5wLOvi5wNVWVul2q+KazTl1REREZLIKVGrEJilwPOEhjic8xOnEdOQWqLXKeDlYodNzrv+bF+cCTwdrI0VbNZjUERERUbWnVgtcvafE8YTC1amnbqXjUW6BVhlXW5kmgevk74IGznUgkRhvOLSqManTEefUERERVZ0nN/w9duMhTtx6iIysfK0yDtYW6ODnjE7+hb1xz7nb1qok7mmcU1dOnFNHRESkf7ps+GtjaYZ2voVJXEd/FzT1tDfqwgSAc+qIiIiIdN7wt6Nf4V5xwfUdYGHGzo7SMKkjIiKiKlGZDX/p2ZjUERERkUHoc8NfejaeOSIiItILQ2/4S2VjUqcjrn4lIiJTp1ILnLqVjtTMHLjbWaGdr3OZCxGMteEvlYyrX8uJq1+JiMgU7bmcUuxepp5P3cuUG/4Wx9WvREREVG3suZyC6T+cw9O9PPcUOYj44RyGt6qH9Kx8bvhbzTGpIyIiqsVUaoFF2+OKJXQANMdiziVpjnHD3+qLSR0REVEtdupWutaQa2nGtW+A0e0aVIsNf6lkTOqIiIhqoRupmdgRm4L/nLqjU/m2vs4Iqudg4KioMpjUERER1RK3HjzGjovJ2Hkppdiecc/ibmdloKhIX5jUERERmbA7D7Ow41IydlxMQVyKUnPcXCpB10au6N/cEx/v+wupytwS59VJAHg4FG5vQtUbkzodcZ86IiKqKf6WZ2FnbAp2XkpB7N8KzXEzqQSdn3PFgOae6N2sLhzrWAIA7KzMMf2Hc5AAWold0cy5hQMDOY+uBuA+deXEfeqIiKg6Ss7Ixq5LKdgRm4ILdzM0x6USoKO/CwYEe6FPMw8421iW+Hpd9qmj4rhPHREREVXafWWOJpE7e1uuOS6RAO19nTEg2At9gzx0uotD3yBP9Ar0KNcdJah6YVJHRERUg6Rl5mL35cJE7nRiOorG2yQSoG1DZ7wY7Il+zT0qtLDBTCpBR38XPUdMVYVJHRERUTX38FEu9ly5hx0XU3Dy1kOon5g41aqBIwYEe6F/c094OHCFam3GpI6IiKgakj/Ow94r97DzUgqOJTyE6olMroW3IwY090T/YE/UczTte6uS7pjUERERVROKrHzsjbuHnbEpOHrjAQqeSOSC6tljQLAXXmzuCW/nOkaMkqqrWpnU7d69G1FRURBCwNzcHF988QU6duxo7LCIiKgWUubkY3/cfeyITcGf19OQr/onkWvqaY8BwZ54sbknfFxtjBgl1QS1MqmbPHky9u/fj8DAQOzYsQOvvfYaTp8+beywiIiolniUW4ADVwsTucN/pSFPpdY816SuHV4M9sSLwZ7wd7M1YpRU09SYpO7GjRv4+OOPceLECVy+fBkBAQG4fPlysXLx8fGYOXMmjh07Bjs7O0yYMAEffPABLC3/2ZdHKpVCoSjcjFGhUKBevXpV9j6IiKh2ysorwIGrqdgZm4KDf6Uit+CfRM7fzaZwaDXYE43r2hkxSqrJakxSd+XKFezcuRPt27eHWq2GWq0uVkYulyM0NBSNGjXC1q1bkZSUhFmzZiErKwsrV67UlPvhhx8waNAgWFtbQwiBo0ePVuVbISKiWiI7T4VDf6ViR2wKDsTfR07+P9cuH5c6GBDshQEtPNGkrh0kEu4HR5VTY5K6gQMHYvDgwQCAsLAwnDlzpliZ1atXQ6lUYtu2bXB2LrxHXUFBAWbMmIGoqCh4eXmhoKAAy5Ytw/79+9GiRQt8/fXXmDJlCvbt21el74eIiExTTr4Kh6+lYWdsCvZfvY+svH9uL+ntbK1Z7NDMy56JHOlVjUnqdLn1xu7du9GzZ09NQgcAI0eOREREBPbt24ewsDBcuHAB6enpaNGiBQBg7NixeP3110utMzc3F7m5uZrHSmXhzZBL6y3UB7VaDSGEweonIiL9yi1Q4ciNh5pE7lHuP4mcl6MVXmxeuNiheb1/EjkhBHinzpqvKq7ZutZdY5I6XcTHxyM8PFzrmKOjIzw9PREfHw8AqF+/Pm7cuIG7d+/C29sbe/bsQbNmzUqtc+nSpVi0aFGx42lpacjJySnhFZWnVquhUCgghOC9X4mIqql8lRqn72TiwHU5Dt/IwKMneuTcbC3wQiMn9GzsjGYedf6XyOUiLS3NeAGTQVTFNTszM1OnciaV1Mnlcjg6OhY77uTkhPT0dACAh4cHPvnkE/Tt2xfm5uawsbHB2rVrS61z/vz5mDVrluaxUqmEt7c33NzcyrypbmWo1WpIJBK4ubkxqSMiqkYKVGocv/kQOy/dw94r96HIztc8524nQ/8gD/QP9kQrb0dIec/UWqEqrtlWVrrdKcSkkjpdTZo0CZMmTdKprEwmg0wmQ3R0NKKjo6FSFf4lJpVKDZpwSSQSg7dBRETPplILnLz5ENtjU7D3yj2kP87TPOdqa4l+QZ4YEOyJNj7OMGMiVysZ+pqta70mldQ5OTlptip5klwu15pnVxGRkZGIjIyEUqmEg4NDpeoiIqLqTaUWOJ2Yjp2xKdh9OQUPHv2TyDnbWKJvkAcGNPdEez8XJnJUbZhUUhcQEKCZO1dEoVAgJSUFAQEBlar76Z46IiIyLWq1wLk7cuyITcGuSylIzfxnkZxjHQv0beaBF4M90dHPBeZmHEWh6sekkrp+/fphyZIlyMjI0Myt27x5M6RSKXr37l2putlTR0RkeoQQOH83Azv/l8ilKP5ZAGdnZY4+zTwwINgTnZ9zhQUTOarmakxSl5WVhV27dgEAbt++DaVSiZiYGABA9+7d4ebmhoiICHzxxRcYMmQIoqKikJSUhNmzZyMiIgJeXl7GDJ+IiKoJIQQuJSmwIzYFO2NTkJSRrXnOVmaO3oF18WKwJ7o0coXM3MyIkRKVT41J6lJTUzFixAitY0WPDx48iJCQEDg5OeHAgQOYOXMmhgwZAjs7O0yZMgWLFy+udPscfiUiqrmEELiSrMTOS4WJ3J30LM1zdSzN0LNpXQwI9kS3xm6wsmAiRzWTRHDnw3IpGn5VKBQG3dIkNTUV7u7uXP1KRFRBQgjE38vEztgU7LyUglsPHmues7YwQ2hTdwwM9kRIE3cmclRhVXHN1jX3qDE9dURERLq4fj8T22NTsDM2GQlp/yRyMnMpQgPc8WKwJ0ID3FHHkpdAMi38ROuIw69ERNVXQtoj7IxNwY7YZFy7/0hz3NJMiu5N3DAg2BMvNK0LWxkve2S6OPxaThx+JSIyLJVa4NStdKRm5sDdzgrtfEve1DfxwWPsvJSCHbEpuJqi1By3MJOgWyM3vBjsiZ6BdWFvZVGV4VMtw+FXIiKiEuy5nIJF2+O0thbxdLDCwoGB6BvkibvpWYWrVi8l43LSP4mcuVSCLo1c8WJzT/QO9IBDHSZyVPswqdMRh1+JiAxrz+UUTP/hHJ4ePkpR5CDih3No6FIHtx/+s2rVTCpBJ38XDAguTOScbCyrNmCiaobDr+XE4VciIv1TqQW6fPhfrR66kkgAdPR3wYvBnujbzAMutrKqCZCoFBx+JSIiesKpW+nPTOgA4MtxrdAvyLMKIiKqedgNRERERnfy1kOdyuUVqA0cCVHNxZ46IiIyCiEEDl9Lw6pDCTh5K12n17jbWRk4KqKai0mdjrhQgohIPwpUauy6fA+rDiVotiIxlwIWZmbIzi/5O1YCwMOhcHsTIioZkzodRUZGIjIyUjNZkYiIyicnX4Ut5/7GV3/c1KxirWNphjHtGmByV19cvJuB6T+cAwCtFbBFO9QtHBhY4n51RFSISR0RERlUZk4+fjhxB2uP3MKDR7kAAKc6Fgjr5IsJHRtqtiLxdLDGqnGtiu1T5/HEPnVEVDomdUREZBBpmbn47ugtbDhxG5k5BQAALwcrTO3mh1FtvUu892rfIE/0CvTQ6Y4SRKSNSR0REenV3fQsfPXHTfx85i5y/7da9Tl3W0R098egFl6wNC974wUzqQQd/V2qIlQik8KkTkdcKEFEVLb4e0qsOpSAHbEpUKkLZ8W18HbEjBB/9GpaF1L2thEZFO8oUU68owQRkbbTielYdSgB/41P1Rzr2sgVM0KeQwc/Z0gkTObIdPGOEkREVKMJIXDwr1R8eTABZ27LAQASCdC/uSemd/dHUD3uEkBU1ZjUERGRzgpUauyITcHqwwmIv5cJALA0k2JY63p4pZs/fF1tjBwhUe3FpI6IiJ4pJ1+FzWfuYs0fN/G3PBsAYGNphrEdGmJyF1/UteedHoiMjUkdERGVSpGdjx9O3MZ3R2/hwaM8AICzjSXCO/tgfAcfONSxMHKERFSESR0RERWTqszB2qO3sPHEHTzKLdxjrp6jNV7p5oeRbbxhbWlm5AiJ6GlM6oiISCPxwWN89edNxJz9G3n/22OucV1bTA/xx4BgL1iYcUU+UXXFpE5H3KeOiEzZ5SQFVh9OwK5LKfjfFnNo1cARM0KeQ2iAO/eYI6oBuE9dOXGfOiIyFUIInLxVuMfc4WtpmuMhTdwwI+Q5tPVx4h5zRM/AfeqIiMho1GqB/VfvY9XhBJy/kwEAkEqAAcFeiOjuj0Avw/zBSkSGxaSOiKiWyFep8duFZKw+nIDrqY8AAJbmUoxoXR+vdPNDQxfuMUdUkzGpIyIycdl5Kvx0+g6+/vMWkjIK95izk5ljXMeGmNTZB+523GOOyBQwqSMiMlGKrHx8fzwR644lIv1x4R5zrraWCO/ii3EdGsLeinvMEZkSJnVERCbmniIHa4/cxKaTd/A4r3DFvrezNaZ188fw1vVhZcE95ohMUa1L6lJTU9G7d2/N4wcPHsDNzQ3nz583YlRERJV3M+0R1hy+ia3n/0a+qnBjgwAPO0wP8ceLzT1hzj3miExahZO6QYMGVeh1K1asgJ+fX0WbrTR3d3dcuHBB83js2LEIDg42WjxERJV16W8FVh2+gd2X76Fok6p2Ps6YHuKPkCZu3JaEqJao8J9tO3bsQFJSEjIzM3X6p1QqsXPnTmRkZFSovRs3biAiIgItW7aEubk5goKCSiwXHx+PXr16wcbGBh4eHpgzZw7y8vJKLJuZmYnffvsN48aNq1BMRETGIoTAsRsPMO6bkxi48gh2XSpM6F4IcEdMREf8HNERPQLcmdAR1SKVGn5dtWoV2rVrp1PZgoICWFpaVritK1euYOfOnWjfvj3UajXUanWxMnK5HKGhoWjUqBG2bt2KpKQkzJo1C1lZWVi5cmWx8jExMejYsSPq1atX4biIiKqSWi2wL65wj7mLdzMAAGZSCQa18MK07n4I8OAec0S1VYWTurlz55YrGTIzM8PcuXPh6elZofYGDhyIwYMHAwDCwsJw5syZYmVWr14NpVKJbdu2wdnZGUBhMjljxgxERUXBy8tLq/z69esxZcqUCsVDRFSV8grU+OVCElYfTsDNtMcAAJm5FKPaemNqVz94O9cxcoREZGw18jZhRUnd5cuXtY5369YNzs7O+OWXXzTHMjIy4OzsjG+//RZhYWGa43fu3EFwcDCSk5NRp47uX4a8TRgRVaXHuQX4z+m7+ObPm0hR5AAA7KzMMbGjD8I6+8DVVmbkCIlqN5O7Tdjq1asxZswYgyU5uoqPj0d4eLjWMUdHR3h6eiI+Pl7r+IYNGzBs2LBnJnS5ubnIzc3VPFYqlQBQ6hCwPqjVagghDFY/EVV/8qw8fH/sNtYfv42M7HwAgJudDJM7++Dldt6w+98ec/yeIDKuqrhm61q3XpK6WbNm4c0338SwYcMwefJkdO/eXR/VlptcLoejo2Ox405OTkhPT9c6tmHDBqxZs+aZdS5duhSLFi0qdjwtLQ05OTkVjrUsarUaCoUCQgj21BHVMvcz8/Djufv45dID5BQUfpHXd5BhXJu66NfUBTJzKbKVcmQrjRwoEQGommt2ZmamTuX0ktQlJydj48aN+O6779CjRw/4+fkhPDwcEydOrLaLEJ7uuSvN/PnzMWvWLM1jpVIJb29vuLm5GXT4VSKRwM3NjUkdUS2RkPoIq/+4iV8vJKNAXTgrppmXPSK6+aFvkAfMpFzFSlQdVcU128pKt1v56SWpc3R0RGRkJCIjIxEbG4tvv/0WK1aswMKFC9G7d2+Eh4dj0KBBsLAw7C1pnJycoFAoih2Xy+WahRPlJZPJIJPJEB0djejoaKhUhbuzS6VSgyZcEonE4G0QkfFduJuBVYduYF/cfc0ecx38nDE95Dl0a+TKLUmIagBDX7N1rVfvrQcHB2PFihW4cOECOnfujN27d2PEiBGoV68eFi5ciOzsbH03qREQEFCsB06hUCAlJQUBAQGVqjsyMhJxcXE4ffp0peohIhJC4M/raRjz9QkMiT6KvVcKE7pegXWxdUYn/OeVjujemJsGE1H56PU2YUII7NmzB2vXrsWOHTvg6OiI2bNnY+jQodi1axe++OILXL58GVu2bNFnsxr9+vXDkiVLkJGRoZlbt3nzZkilUq1bg1XE0z11RETlpVIL7Ll8D6sO38DlpMJJceZSCQa3rIeI7n5oVNfOyBESUU2mly1NEhIS8O2332L9+vVITk5Gr169MHXqVAwePBjm5v/kjb/99htGjx6NrKyscreRlZWFXbt2AShMsBISEvDpp58CALp37w43NzfI5XI0a9YMjRs3RlRUlGbz4bFjx5a4+XBFcEsTIiqv3AIVtp1Lwpo/buLWg8I95qwspBjdtgGmdPVFfSfuMUdUU5ncliaNGjVCvXr1MGnSJEyePBkNGzYssVxAQADat29foTZSU1MxYsQIrWNFjw8ePIiQkBA4OTnhwIEDmDlzJoYMGQI7OztMmTIFixcvrlCbT2JPHRGV16PcAmw6eRvf/HkLqZmFWyM5WFtgYicfhHXygbNNxe+yQ0T0NL301O3YsQP9+/evFb1K7KkjIqBwKPXUrXSkZubA3c4K7XydNStUHz7Kxbpjifj+WCKUOQUAgLr2Mkzt6ofR7RrAVqbXmS9EZEQm1VOXl5eHvLw8JCYmws/Pr7LVERFVe3sup2DR9jjNHR4AwNPBCq/2eA7XUx/hP6fvICe/cI85P1cbTOvuhyHP14PM3MxYIRNRLVDppM7S0hJjxozBnj17TDqp4/ArEQGFCd30H87h6SGOFEUOFvzyz60Lm9dzwIwQf/Ruxj3miKhq6GUMICAgAHfu3NFHVdVW0T58RV2gRFT7qNQCi7bHFUvonmRpLsU349uga2PuMUdEVUsvg79Lly7FBx98gDNnzuijOiKiaunUrXStIdeS5BWoYWEuZUJHRFVOLz11c+bMwcOHD9G+fXu4uLigbt26Wl9oEokEFy9e1EdTRERGk5qp2/2edS1HRKRPeknqWrdujTZt2uijqmqLc+qIarfsPBX2Xr6nU1l3O93u00hEpE962dKkNuGWJkS1z+nEdMyJidVsHFwaCQAPByscmRvKxRFEtUR12tKEGQMRUSmy81R4b3scRq45jlsPHsPD3gqRPfwhQWEC96SixwsHBjKhIyKj0NsOmBs2bMCaNWtw7do15OQUn0+iVCr11RQRkcGdupWOOTEXkfiw8LaGI9vUx4IXA+FgbYHm9RyK7VPn4WCFhQMD0TfI01ghE1Etp5ek7ocffsDUqVMRFhaGY8eOITw8HCqVCtu3b4ejoyMmTJigj2aMinPqiGqH7DwVPtobj3XHEiEE4GFvhaXDmqNHE3dNmb5BnugV6FHqHSWIiIxBL3Pqnn/+eQwfPhzz5s2DhYUFzpw5g1atWiEzMxO9e/fGiBEjMGvWLH3Ea3ScU0dkuk7dSsfsmIu4/UTv3L8HBMLeysLIkRFRdWVyc+quX7+Ozp07w8zMDGZmZpqhVjs7O8ydOxeff/65PpohIjKIrLwCvPvbFYz66jhuP8yCp4MV1k1qi4+Gt2BCR0Q1hl6GXx0cHJCbmwsAqFevHuLi4hASEgIAUKlUePjwoT6aISLSu5M3H2LOllhN79yoNt5YMKApkzkiqnH0ktS1adMGsbGx6NOnDwYNGoRFixZBrVbDwsICy5YtQ4cOHfTRDBGR3mTlFeCjPX9h3bFEAICngxWWDQtG98Zuxg2MiKiC9JLUzZ8/H7dv3wYAvPfee7h9+zbeeOMNqNVqtG3bFmvWrNFHM0REenHi5kPMiYnFnfTC3rnRbb0R9SJ754ioZqt0UpeXl4fk5GS0bdsWAODo6Ihff/0Vubm5yM3NNdhigqrG1a9ENd/TvXNe/+ud68beOSIyAXpZ/WplZYU9e/Zo5tGZMq5+JaqZjic8xNwt//TOvdzOG/P7s3eOiCqnOq1+1cvwa0BAAO7cuaOPqoiI9OpxbgE+2hOP748XThFh7xwRmSq9JHVLly7F66+/jsDAQLRp00YfVRIRVdrxhIeYs+Ui7qZnAwBebtcAUf0DYMfeOSIyQXpJ6ubMmYOHDx+iffv2cHFxQd26dSGR/LOzukQiwcWLF/XRFBHRMz3OLcCHe+Kx/n+9c/UcrbFsWHN0bcTeOSIyXXpJ6lq3bs0eOiKqFo4lPMDcLbHsnSOiWkcvSd26dev0UQ0RUYU9zi3Ast3x2HDin965D4cFo0sjVyNHRkRUNfSyTCM8PBy3bt0q8bnbt28jPDxcH80QEZXo2I0H6LPiD01CN7Z9A+z9VzcmdERUq+glqVu3bh3S0tJKfO7Bgwf4/vvv9dGMUUVHRyMwMFCzHx8RGd/j3AL8+5dLGPPNSfwtz0Y9R2tsnNIei4c2h61MLwMRREQ1ht6+9Z5cGPGk69evw8XFRV/NGE1kZCQiIyM1e8UQkXEdu/EAc7bE4m954dy5se0bYH7/pkzmiKjWqvC336pVq7Bq1SoAhQndmDFjYG1trVUmJycHiYmJGDFiROWiJCL6n0e5BVi66yo2nizcG7OeozU+Gh6Mzs9xqJWIarcKJ3VeXl5o3bo1AODy5cto0qQJ3Ny0twuwtLRE06ZNMXny5MpFSUQE4OiNB5gTE4ukjMLeuXEdGmBeP/bOEREBlUjqBg8ejMGDB2sev/POO/D19dVLUERET3qUW4Alu65i0/965+o7WeOjYcHoxN45IiINvfx5+9133+mjGiKiYo5cL9x3rqh3bnyHhpjXLwA27J0jItLCb0UiqpYyc/KxdHe8du/c8GB08mfvHBFRSSqc1AUHB2PTpk0ICgpC8+bNS139ClS/24RlZ2fjtddew+HDh2FhYYFhw4bhvffeM3ZYRPQ/f15Pw7wtlzS9cxM6NsTcvuydIyIqS4W/IVu3bg0bGxvN/5eV1FU3b731Fjw8PHDt2jUAwP37940cEREBhb1zS3bF48dThb1z3s6Fd4Vg7xwR0bNJhBDC2EHo4saNG/j4449x4sQJXL58GQEBAbh8+XKxcvHx8Zg5cyaOHTsGOzs7TJgwAR988AEsLS0BAI8ePYK/vz/u3LkDmUxW7jiK9qlTKBSwt7ev9PsqiVqtRmpqKtzd3SGV6mV/aKJq7+neuYkdG2IOe+eIqJqrimu2rrmH3r4tHzx4gM8++wwnT55ESkoKPD090aFDB7zxxhtwda38X9lXrlzBzp070b59e6jVaqjV6mJl5HI5QkND0ahRI2zduhVJSUmYNWsWsrKysHLlSgBAQkIC3N3d8dZbb+HYsWNwdXXF8uXLERwcXOkYiaj8CnvnruLHU3cBAA2c6+DDYcHo6F/zNy0nIqpKekkpT548iUaNGmHlypVwcHBA9+7d4eDggC+++AL+/v44efJkpdsYOHAg7t69i5iYGLRq1arEMqtXr4ZSqcS2bdvQp08fhIeH46OPPsLq1auRnJwMACgoKMDly5fRs2dPnD17Fm+++SaGDBlS6fiIqPz+uJaGPp/9oUnoJnZsiD1vdGVCR0RUAXoZfm3Tpg2srKywa9curW5BhUKBfv36IT8/H6dPn65sMxphYWE4c+ZMseHXbt26wdnZGb/88ovmWEZGBpydnfHtt98iLCwMaWlp8PHxwePHjzVl3N3dERcXV2KPYm5uLnJzczWPlUolvL29IZfLDTr8mpaWBjc3Nw6/kklS5uRj6a54/HTmbwBAA2drfPhSc7T3YzJHRDVLVVyzlUolnJycqmb49cqVK9i8eXOxhhwcHDBv3jyMGjVKH808U3x8PMLDw7WOOTo6wtPTE/Hx8QAANzc3dOjQAUeOHEGXLl1w5swZyGSyUu9Pu3TpUixatKjY8bS0NOTk5Oj/TaDwA6JQKCCEYFJHJudEogJL9t9G6qN8AMCIlm6Y0bkerC1USE1NNXJ0RETlUxXX7MzMTJ3K6SWpe+6555CRkVHicwqFAn5+fvpo5pnkcjkcHR2LHXdyckJ6errm8erVqxEeHg6FQoE6depg8+bNpa7enT9/PmbNmqV5XNRT5+bmZtCeOolEwp46MinK/61s/Zm9c0RkQqrimm1lZaVTOb0kdcuXL0dkZCS8vb3RvXt3zfFDhw7h3Xff1SxSqC4aNWqEP//8U6eyMpkMMpkM0dHRiI6OhkqlAgBIpVKDJlwSicTgbRBVlUN/pWL+1ktIURT2bod18sGcvk1Qx5IrW4mo5jP0NVvXeiv8jfr0hsMKhQKhoaFwcHCAm5sb0tLSoFAo4OTkhLlz56Jfv34VbUpnRePNT5PL5XB2dq5U3ZGRkYiMjNQsKyaiZ1Nk52PxzjhN71xDlzr4aFgwe+eIiAygUpsPP5nUtW7dWi8BVUZAQIBm7lwRhUKBlJQUBAQEVKrup3vqiKhsB/9Kxfwtl3BPmQOJpLB3bnYf9s4RERlKhb9d161bBwDIy8vDjh070LJlyyqbO1eafv36YcmSJcjIyNDMrdu8eTOkUil69+5dqbrZU0ekG0V2Pj7YEYfNZwt753xc6uCj4S3QzrdyveVERFS2Sv/JbGlpiTFjxmDPnj0GTeqysrKwa9cuAMDt27ehVCoRExMDAOjevTvc3NwQERGBL774AkOGDEFUVBSSkpIwe/ZsREREwMvLy2CxEVGhg/GFc+eKeucmdfLF7D5NYG1pZuzQiIhMnl7GQQICAnDnzh19VFWq1NRUjBgxQutY0eODBw8iJCQETk5OOHDgAGbOnIkhQ4bAzs4OU6ZMweLFiyvdPodfiUqnyM7H+zviEPNE79zyES3Q1oe9c0REVUUvmw/v3r0br7/+OjZt2oQ2bdroI65qi/d+JdJ2MD4V87bG4r4yFxIJEN7ZF2/1Zu8cEdUOJnfv1zlz5uDhw4do3749XFxcULduXa1FFBKJBBcvXtRHU0RUTSiy8vHejjhsOVfYO+fraoOPhgezd46IyEj0ktS1bt3a5HvoOPxK9I//xt/H/K2XNL1zkzv74k32zhERGZVehl9rEw6/Um1WUu/c8uHBaMPeOSKqpUxu+JWITB9754iIqjcmdTri8CvVVoqsfCzacQVbzyUBAPxcbbB8RDBaN2TvHBFRdaK34dcNGzZgzZo1uHbtGnJycoo9r1Qq9dGM0XH4lWqT/XH3EbXtElIzC3vnpnQp7J2zsmDvHBERYILDrz/88AOmTp2KsLAwHDt2DOHh4VCpVNi+fTscHR0xYcIEfTRDRFVEkZWPRduvYOt59s4REdUUekkpP/nkE7z99tuIjo4GAMyYMQPfffcdbt26BTc3N9ja2uqjGSKqAvvj7qPXZ4ex9XwSpBLglW5+2PV6VyZ0RETVnF6SuuvXr6Nz584wMzODmZmZZqjVzs4Oc+fOxeeff66PZojIgDKy8vCvny5gyvozSM3MhZ+bDTZHdEJU/6YcbiUiqgH0ktQ5ODggNzcXAFCvXj3ExcVpnlOpVHj48KE+mjGq6OhoBAYGom3btsYOhUjvfo+7j16f/YFt/+udm9bND7te64rWDZ2MHRoREelIL3Pq2rRpg9jYWPTp0weDBg3CokWLoFarYWFhgWXLlqFDhw76aMaoIiMjERkZqZmsSGQKMrLy8O5vV/DLhWQAgL+bDZaPaIFWDZjMERHVNHpJ6ubPn4/bt28DAN577z3cvn0bb7zxBtRqNdq2bYs1a9booxkiKieVWuDUrXSkZubA3c4K7XydYSYtvIXfviv3ELXtMh48yoVUAkzt5od/9WzMoVYiohpKL0ldhw4dNL1xjo6O+PXXX5Gbm4vc3FyDbftBRGXbczkFi7bHIUXxzxZDng5WeLNXY/x54wF+faJ37uMRLfA8e+eIiGo0g20+LJPJIJPJDFU9EZVhz+UUTP/hHJ7ehDJFkYO3YmIBgL1zREQmpsILJZYsWYKUlJRyv+b+/fsVbdKouFCCagqVWmDR9rhiCd2TzKUS/DytI+b348pWIiJTUeGk7u2338bff/+tc3mVSoW3334bSUlJFW3SqCIjIxEXF4fTp08bOxSiMp26la415FqSArVAvkovN5MhIqJqosLDr0IIvPnmm3B0dNS5PBEZXmpm2QldecsREVHNUOGkrlu3bpBIJMjMzCzXa+zs7CraJBHpwN3OSq/liIioZqhwUnfo0CE9hkFE+tK6oRPqWJohK09V4vMSAB4OhdubEBGR6dDLHSWIqHrILVDhjZ/Ol5nQAcDCgYGa/eqIiMg0MKkjMhGZOfkI+/Y0dl26BwszCaZ09YWng/YQq4eDFVaNa4W+QZ5GipKIiAzFYPvUEVHVSc3MQdi3pxGXooStzBxrxrdG5+dcMb9f01LvKEFERKaFSZ2OoqOjER0dDZWq5GEtImNJfPAYE749hTvpWXC1tcS6Se0QVK/w/sRmUgk6+rsYOUIiIqoKEsG9RspFqVTCwcEBCoXCYLdAU6vVSE1Nhbu7O6RSjpBT6S4nKRD23Sk8eJSHBs51sD68HXxcbYwdFhFRrVEV12xdcw/21BHVUEdvPMAr68/gcZ4KgZ72WBfeltuUEBHVYkzqiGqgHbHJ+NdPF5CvEujk74I141vDzsrC2GEREZERMakjqmG+P5aId7dfgRBA/+Ye+GxUS8jMef9WIqLajkkdUQ0hhMAn+65h5cEbAIAJHRti4cBmXM1KREQAmNQR1QgFKjUWbLuMn87cBQC82asxXg19DhIJEzoiIipUK5M6Hx8f1KlTB5aWlgCATZs2ITAw0MhREZUsJ1+FmT+ex+9x9yGVAIuHNsfL7RoYOywiIqpmDJrU5efn4+rVq7h06RLGjh1ryKbKbdeuXfDx8TF2GERlUmTlY8r60zidKIeluRRfvPw8+jTzMHZYRERUDeltQ5WEhAT88ssv+OCDDzBq1Cg0a9YMNjY2aNmyJSZMmFDp+m/cuIGIiAi0bNkS5ubmCAoKKrFcfHw8evXqBRsbG3h4eGDOnDnIy8urdPtEVe2eIgcj1xzH6UQ57KzMsSG8HRM6IiIqVYV76i5cuIDo6GhcvHgRcXFxyM7OBgBYWloiMDAQ7dq1w+TJkxEcHIzg4OBKB3rlyhXs3LkT7du3h1qthlqtLlZGLpcjNDQUjRo1wtatW5GUlIRZs2YhKysLK1eu1Co7dOhQqNVq9O/fH++99x4sLLgdBFUfN1IfYeK3p5CUkQ13Oxm+D2+Hpp6G2eyaiIhMQ4WTupEjRyI5ORl9+/ZFaGgo1q9fj5ycHKxfvx4DBgzQZ4wAgIEDB2Lw4MEAgLCwMJw5c6ZYmdWrV0OpVGLbtm1wdnYGABQUFGDGjBmIioqCl5cXAODIkSOoX78+Hj16hPHjx+Ojjz7CggUL9B4zUUVcuJuBSd+dgjwrH36uNvg+vB28nesYOywiIqrmKpzUFQ23Dhw4EAAQFRWFuXPnYujQoRg9ejRWrFgBFxf93XNSl1tv7N69Gz179tQkdEBh8hkREYF9+/YhLCwMAFC/fn0AgK2tLSZPnow1a9aUWmdubi5yc3M1j5VKJQCU2luoD2q1GkIIg9VP1dfha2mYsfE8svNVCK7vgLUTWsPFVsbPAhFRNVUV12xd665wUrd37148//zzmsf29vZYtWoVJkyYgGnTpiEgIACffvopxo8fX9Emyi0+Ph7h4eFaxxwdHeHp6Yn4+HgAwOPHj6FSqWBvb4+CggJs2bKlzOHhpUuXYtGiRcWOp6WlIScnR79v4H/UajUUCgWEELz3ay2y++pDfPB7IlRqoH0Deywd4AtVlgKpWcaOjIiISlMV1+zMzEydylU4qevZs2eJxzt27Ihz585h+fLliIiIwMaNG7F69eoqWWkql8vh6OhY7LiTkxPS09MBAPfv38dLL70EtVqNgoICdOrUqcyh1/nz52PWrFmax0qlEt7e3nBzcyvzprqVoVarIZFI4ObmxqSulvjmyC0s2ZsIABjUwhMfDQuGpTl/9kRE1V1VXLOtrHS7r7dBtjQxNzfH/PnzMWrUKEybNg1BQUF49OiRIZoqNz8/P1y4cEHn8jKZDDKZDNHR0YiOjoZKpQJQOBxsyIRLIpEYvA0yPrVa4MM98Vjzx00AwOQuvljQvymkvEsEEVGNYehrtq71VjipW79+vU7lxo8fjx49elS0mXJxcnKCQqEodlwul2vNs6uIyMhIREZGQqlUwsHBoVJ1EQFAvkqNuVtisfVcEgBgXr8ATOvmx7tEEBFRhVQ4qStadFCk6EIkhCh2DChcSGFoAQEBmrlzRRQKBVJSUhAQEGDw9ol0lZVXgMiN53DwrzSYSSVY9lJzjGjjbeywiIioBqtwP6FcLtf8O336NBo2bIh///vfuHjxIu7du4eLFy9iwYIFaNiwIU6ePKnPmEvVr18/7N+/HxkZGZpjmzdvhlQqRe/evStVd3R0NAIDA9G2bdtKRkm1nfxxHsZ8fRIH/0qDlYUUX09ozYSOiIgqTSKe7FqroF69euGFF17AvHnzij23dOlS7N+/HwcOHKhUG1lZWdi1axeAwgQrISEBn376KQCge/fucHNzg1wuR7NmzdC4cWNERUVpNh8eO3Zssc2HK6po+FWhUBh0oURqairc3d05p87EJGVkY8Lak0hIewwHawt8G9YWrRs6GTssIiKqoKq4Zuuae+hlocSxY8cwZ86cEp9r3bo1Pvjgg0q3kZqaihEjRmgdK3p88OBBhISEwMnJCQcOHMDMmTMxZMgQ2NnZYcqUKVi8eHGl2396oQRReV27n4kJa0/hnjIHng5WWB/eDo3q2hk7LCIiMhF66anz9fXFCy+8gG+++abYc+Hh4fjvf/+LxMTEyjZTLbCnjiriTGI6wtedhjKnAI3cbfF9eDt4OVobOywiIqokk+upi4qKwrRp05CQkIAhQ4bA3d0dqamp2LZtG/74448y79hAZOr2x91H5KZzyC1Qo3VDJ6yd2AaOdSyNHRYREZkYvSR1U6dOhaenJxYvXozZs2ejoKAA5ubmaNWqFX799VfNrcRqMg6/UkX8fPou5m+7BJVaIDTAHdFjWsHa0szYYRERkQnSy/Drk9RqNdLS0kz2bggcfiVdCCHw5aEELN/7FwBgeOv6WPpSc1iY8edJRGRKTG749UlSqRR169bVd7VENYZaLfDejjisO5YIAJge4o85fZpwU2EiIjIog9wmjKi2yitQ483NF7H9YjIA4O0BgZjcxdfIURERUW3ApE5HnFNHz/IotwARG87iyI0HMJdK8MnIFhjcsp6xwyIiolpC73PqTB3n1FFJHjzKxaTvTuNSkgJ1LM2welxrdGvsZuywiIjIwExuTp1arcaFCxdw8uRJpKSkIDs7Gy4uLmjSpAm6dOkCNzde3Mh03U3Pwvi1J5H4MAvONpb4LqwtWng7GjssIiKqZSqV1CUkJCA6OhobN25EWloazMzM4OjoCJlMhoyMDGRlZUEikaBr166YOnUqXn75ZfY8kUmJS1Zi4nenkJaZi3qO1tgwuR383GyNHRYREdVCFc6wXnnlFTRr1gwXL17EokWLcOHCBeTk5CAtLQ1///03Hj16hNTUVOzYsQMtWrTAnDlzEBgYiGPHjukz/ioTHR2NwMBAtG3b1tihUDVxPOEhRq05jrTMXAR42GHrjE5M6IiIyGgqPKfulVdewdy5c+Hv769TeZVKhU2bNsHMzAxjxoypSJPVAufUEQDsvpSC1/9zAXkqNdr5OuPrCW3gYG1h7LCIiKiKmcScuq+++qpc5c3MzDB+/PiKNkdUbfxw4jbe/vUyhAD6NKuL/xv9PKwseJcIIiIyLm5pQqQjIQT+78B1rNh/HQDwcrsG+GBIEMyk3FSYiIiMz6Bje/n5+YiNjcXGjRsN2QyRwanUAv/+5bImoXvthUZYMpQJHRERVR9666lLSEjApUuXcPnyZc1/r1+/joKCAkgkEowdO1ZfTRFVqZx8Fd74zwXsuXIPEgnw3qBmGN/Rx9hhERERaalwUnfhwgVER0fj4sWLiIuLQ3Z2NgDA0tISgYGBaNeuHSZPnozg4GAEBwfrLWBj4R0laidlTj5eWX8GJ26mw9JMihWjW6J/c09jh0VERFRMhVe/Nm7cGMnJyejbty+ee+45rF+/Hjk5OVi/fj0GDBig7zirDa5+rT1SlTmY+N1pXE1RwlZmjq8mtEYnf1djh0VERNVIdVr9WuHWExIS8OOPPyImJgbLli1DfHw8Ro0ahaFDh2L8+PF4+PBhRasmMrpbDx5j2OpjuJqihKutDP95pQMTOiIiqtYqnNTt3bsXnTp10jy2t7fHqlWr8Mcff+DixYsICAjAhg0b9BIkUVW69LcCw1cdw930bDR0qYMt0zsiqJ6DscMiIiIqU4WTup49e8LFxaXY8Y4dO+LcuXOYNWsWIiIi0LdvXyQmJlYmRqIq8+f1NIz+6jgePs5DUD17xER0QkMXG2OHRURE9EwGGfw1NzfH/PnzERsbC5VKhaCgIEM0Q6RXv11MRvi603icp0Ln51zw49QOcLOTGTssIiIinRh082F/f3/8/vvv3KeOqr3vjt7Cou1xAIAXgz3x6cgWkJnzLhFERFRzVLinbtWqVcjNzdWpbNEedZcuXcJ///vfijZJpHdCCHy0J16T0E3s2BBfjH6eCR0REdU4FU7q1q1bh4YNG+Jf//oXjh07hvz8/BLLJScnY+3atejZsyc6deoEuVxe4WCNKTo6GoGBgWjbtq2xQ6EKUqkFjic8xK8XknA84SFy81WYuyUWXx5KAADM7tME7w5qBinvEkFERDVQhfepA4Bt27bh//7v//Dnn3/CwsICjRs3hpubG2QyGTIyMnDr1i2kpqbC2dkZEydOxFtvvQUPDw99xl/luE9dzbTncgoWbY9DiiJHc0xmLkVugRpSCbBkaHOMbtfAiBESEVFNVJ32qavUnLqhQ4di6NChSExMxIEDB3D69GmkpKQgJycHDRs2RO/evdG5c2eEhITAwsKiMk0RVdieyymY/sM5PP3XS26BGgDwSjc/JnRERFTj6WWhRIsWLbBkyRKsXr1aH9UR6Y1KLbBoe1yxhO5Jv15Ixuw+ATDjsCsREdVgeknqMjMz8eeff+LBgwfIyMiAl5cXXnzxRQQGBuqjeqIKO3UrXWvItSQpihycupWOjv7F910kIiKqKfS2pcnmzZthbm4OT09PPHjwAPPmzUOXLl2wYcMGNGjAoS0yjtTMshO68pYjIiKqrvQ2o2/QoEFIT09HYmIiFAoFfv/9d+Tn56NVq1a4efOmvpohKhd3Oyu9liMiIqqu9JLUSaVSTJ8+HTY2hbdTMjMzQ2hoKI4cOYLg4GD8+9//1kczejd16lRIJJxHZcra+TrDTlZ6h7QEgKeDFdr5OlddUERERAagl6Subt26+Pvvv4tXLpXitddew969e/XRjF7t37+/1L31yHT8eT0NmbkFJT5XlM4vHBjIRRJERFTj6SWp69evHxYuXFjiMGt2djays7Mr3caNGzcQERGBli1bwtzcvNT7ycbHx6NXr16wsbGBh4cH5syZg7y8PK0yjx8/xoIFC/DJJ59UOi6qvm4/fIzXfjwPAOjynCs8HbSHWD0crLBqXCv0DfI0RnhERER6pZeFEsuWLcO5c+cQHByMiRMn4oUXXoCHhweuX7+OqKgodOzYsdJtXLlyBTt37kT79u2hVquhVquLlZHL5QgNDUWjRo2wdetWJCUlYdasWcjKysLKlSs15ebPn4/XX38dLi5c7WiqHucW4JX1Z6HMKcDzDRyxNqwNzKVSnLqVjtTMHLjbFQ65soeOiIhMhV6SOldXV5w8eRLLli3DDz/8gFWrVkEikUAIgUaNGmHVqlWVbmPgwIEYPHgwACAsLAxnzpwpVmb16tVQKpXYtm0bnJ0L50gVFBRgxowZiIqKgpeXF44ePYqEhAR8/vnnlY6JqichBOZsicVf9zPhZifD6nGtNfdy5bYlRERkqvS2pYmFhQXefvttvP3227h58ybu3r0LR0dHNG/eXC+3zdCljt27d6Nnz56ahA4ARo4ciYiICOzbtw9hYWH4888/cf78efj4+GjK+Pj44MSJEyXewiw3Nxe5ubmax0qlEgBK7S3UB7VaDSGEweo3dV/9cRM7Y1NgYSZB9Jjn4WZryXNJREQGURXXbF3r1ltS9yQ/Pz/4+fkZouoyxcfHIzw8XOuYo6MjPD09ER8fDwCYN28e5s2bp3leIpEgMTGx1DqXLl2KRYsWFTuelpaGnBzD7G2mVquhUCgghOC9X8vp5G0lPtp7HQDwr+7eaGCdj9TUVCNHRUREpqoqrtmZmZk6lTNIUmcscrkcjo6OxY47OTkhPT29QnXOnz8fs2bN0jxWKpXw9vaGm5tbmTfVrQy1Wg2JRAI3NzcmdeVwNz0LC/fEQi2AkW3qY9oLgdyyhoiIDKoqrtlWVrrtpWpSSV1FCFHWXUEBmUwGmUyG6OhoREdHQ6VSASgcDjZkwiWRSAzehinJzlMhYuN5ZGTno4W3I94bHAQzMzNjh0VERLWAoa/ZutZbpRnD/fv3DVq/k5MTFApFseNyuVxrnl1FREZGIi4uDqdPn65UPaR/QgjM3RKLqylKuNpaYvW4VrCyYEJHRES1i0F66vLy8nDlyhXExsZq/QMMm9gFBARo5s4VUSgUSElJQUBAQKXqfrqnjqqPtUdu4beLyTCXShA9phU8HayNHRIREVGV03tP3bRp0+Dg4IBevXohIiICdnZ2mDVrFs6fP2/wnrp+/fph//79yMjI0BzbvHkzpFIpevfuXam62VNXPR278QBLdl0FALw9IBDt/bhlCRER1U56T+p+//13pKamajb+3bJlC5ydneHl5VWperOyshATE4OYmBjcvn0bSqVS8zgtLQ0ANEnkkCFDsG/fPnz33XeYPXs2IiIiKt0+VT9/y7MQuekc1AIY3ro+JnRsaOyQiIiIjEYinrVSoJzS0tLg5uameXz+/HlMnToVoaGh+OCDD2BpaVmhehMTE+Hr61vicwcPHkRISAgA4OrVq5g5cyaOHTsGOzs7TJgwAYsXL65wu0WeHH69du0aFAqFQVe/pqamwt3dnQslSpGTr8KwVcdwJVmJ4PoO+HlaR86jIyKiKlcV12ylUgkHB4dn5h56T+pKolKpsGTJEmzZsgUXLlwwdHMGpeuJrQwmdWUTQmDWzxex7XwSXGwssX1mF3g5ch4dERFVveqU1FW49aZNm+LatWs6lTUzM8Pbb7+NH3/8saLNEWl8dzQR284nwUwqwcoxrZjQERERoRJJ3V9//aW1IEEIgSVLlmgde1rTpk0r2pzRRUdHIzAwEG3btjV2KLXa8YSHWPy/hREL+jflvVyJiIj+R2/9hGq1WnPfV1PE1a/Gl5SRjVc3nYNKLTD0+XqY1NnH2CERERFVG3od/K2C6XlUS+XkqxCx4SwePs5DMy97LH2pOW8BRkRE9IRKJXW16aLK4VfjEUJgwbbLuJSkgFMdC6wZ35orXYmIiJ5SqTtKREREoE2bNggKCkJAQAAkEonJJnqRkZGIjIzUrEChqrPhxG1sOfc3pBIgekwr1HeqY+yQiIiIqp0KJ3ULFixAbGws9u3bh6+//lpzfPDgwWjZsiWaN2+O4OBgNG/eHAEBAdyagyrk1K10vLc9DgAQ1b8pOj3nauSIiIiIqqcKJ3Xvv/++5v+VSqXWPV4vXryIzz//HI8fPwYAWFlZISsrq/LRUq2SosjGjI1nUaAWGNTCC5O7lLz5NBEREVVy+LWIvb09unTpgi5dumgdT0hIwIULF3Dp0iV9NEO1SE6+ChE/nMODR3lo6mmPD4cFm+zQPhERkT7oJakrjb+/P/z9/TFs2DBDNlMlnrxNGBmWEALv/HoZF+9mwLGOBb4a3xrWllwYQUREVBadJ7qlpaVhwoQJ6NWrlyHjqba4T13V2XjyDn4+U7gw4ouXn4e3MxdGEBERPYvOSd3o0aOxZcuWGn1XCKr+ziSmY9H2KwCAuX0D0LWRm5EjIiIiqhl0TurOnTuHHj16YMWKFQYMh2qz+8ocTN94DvkqgReDPfFKNz9jh0RERFRj6Dyn7uWXX8bq1avh5OSE+vXrw97eHra2trCzs9P5vw4ODvD29oaFhYUh3xPVQLkFKkT8cBZpmbkI8LDD8uFcGEFERFQeOid1n332GRwcHLBz504kJyfjxo0byM/PL/M1JV2Ura2tMWfOHLzzzjvlj9aIuFDCsN79LQ7n72TA3soca8a3Rh1Lg67hISIiMjkSUYkbtubn5yMzMxOPHj3S6b9yuRxHjx7FzZs3sW/fPrzwwgv6fC9VouiOEgqFAvb29gZpQ61WIzU1Fe7u7rVi0+ZNJ+8gatslSCTAd2FtEdLE3dghERER6aQqrtm65h6V6g6xsLCAs7MznJ2ddX5NYmIi/Pz8cPz48RqZ1JF+nb0tx8LfLgMA3urdhAkdERFRBVVpN9Dhw4fRrFkzyGQyDB48uCqbpmooVZmD6T+cRb5KoH9zD8wI8Td2SERERDVWlSZ17u7umDhxIk6dOoXmzZtXZdNUzeQVqDFj4zmkZuaicV1bLB/eggsjiIiIKqFKZ6M3bdoUX375ZVU2SdXU+zvicOa2HHZW5lgzvg1sZFwYQUREVBmmPwufqp2fT9/FhhO3IZEA/ze6JXxdbYwdEhERUY3HpI6q1IW7Gfj3L4ULI2b1bIzQgLpGjoiIiMg0MKnTUXR0NAIDA9G2bVtjh1JjpWXmImLDWeSp1OgdWBeRPZ4zdkhEREQmg0mdjiIjIxEXF4fTp08bO5QaKV+lRuTGc7inzIG/mw0+GdkCUikXRhAREekLkzqqEot3XsWpxHTYyczx1YQ2sLPireKIiIj0iUkdGVzM2b+x7lgiAOCzUS3h72Zr3ICIiIhMEJM6MqhLfysQte0SAOCNno3QM5ALI4iIiAyBSR0ZzINHuZi24QzyCtTo2dQdr4U2MnZIREREJotJHRlE0cKIZEUO/Fxt8OmollwYQUREZEC1bht/tVqNzp07Izs7GyqVCgEBAVi7di3s7e2NHZpJWborHidvpcPG0gxfTWgNey6MICIiMqha11MnlUqxd+9eXLhwAZcuXYK3tzeWL19u7LBMyrbzf+Pbo7cAAJ+Oaonn3O2MHBEREZHpqzFJ3Y0bNxAREYGWLVvC3NwcQUFBJZaLj49Hr169YGNjAw8PD8yZMwd5eXlaZYp65dRqNR4/fswbyevR5SQF5m0pXBjxWuhz6NPMw8gRERER1Q41Zvj1ypUr2LlzJ9q3bw+1Wg21Wl2sjFwuR2hoKBo1aoStW7ciKSkJs2bNQlZWFlauXKlVtmfPnrhw4QKCgoLwySefVNXbMGnpj/MwbcNZ5Bao0aOJG97o2djYIREREdUaNaanbuDAgbh79y5iYmLQqlWrEsusXr0aSqUS27ZtQ58+fRAeHo6PPvoIq1evRnJyslbZ/fv34/79+2jTpg2+/PLLqngLJq1Apcarm84hKSMbPi51sGL081wYQUREVIVqTE+dVPrs/HP37t3o2bMnnJ2dNcdGjhyJiIgI7Nu3D2FhYVrlzczMEBYWhtGjR2POnDkl1pmbm4vc3FzNY6VSCQCl9hbqg1qthhDCYPUbwrLd8TiW8BB1LM2welwr2MnMalT8REREFVEV12xd664xSZ0u4uPjER4ernXM0dERnp6eiI+PBwA8ePAAAODq6gohBGJiYkqdnwcAS5cuxaJFi4odT0tLQ05Ojh6j/4darYZCoYAQQqdk1tj2xafjmyOFCyP+3ashHCXZSE3NNnJUREREhlcV1+zMzEydyplUUieXy+Ho6FjsuJOTE9LT0wEA9+/fx4QJE5Cfnw8hBJo1a4bPP/+81Drnz5+PWbNmaR4rlUp4e3vDzc3NYNugqNVqSCQSuLm5Vfuk7mqKEksO3AYATO/uh9Gdmxg5IiIioqpTFddsKysrncqZVFKni2bNmuHs2bM6l5fJZJDJZIiOjkZ0dDRUKhWAwuFgQyZcEonE4G1UlvxxHiI2nkNOvhrdG7vhrT4BnEdHRES1jqGv2brWW30zhgpwcnKCQqEodlwul2vNs6uIyMhIxMXF4fTp05Wqx1So1AKv/ec87qZno4FzHXw++nmYMaEjIiIyGpNK6gICAjRz54ooFAqkpKQgICDASFGZpuV7/8Kf1x/A2qLwjhEOdXjHCCIiImMyqaSuX79+2L9/PzIyMjTHNm/eDKlUit69e1eq7ujoaAQGBqJt27aVjLLm2xGbjNWHEwAAy0cEI8CDt1gjIiIythqT1GVlZSEmJgYxMTG4ffs2lEql5nFaWhoAICIiAnZ2dhgyZAj27duH7777DrNnz0ZERAS8vLwq1T6HXwvF31Ni9uZYAMC07n4YEFy580pERET6UWMWSqSmpmLEiBFax4oeHzx4ECEhIXBycsKBAwcwc+ZMDBkyBHZ2dpgyZQoWL15c6fafXihRG2Vk5eGV9WeRna9C10aumNOHQ9pERETVhUQIIYwdRE2iVCrh4OAAhUJh0C1NUlNT4e7uXm1Wv6rUAuHrTuPwtTTUd7LG9le7wMnG0thhERERGVVVXLN1zT2qR8ZA1d6nv/+Fw9fSYGUhxVfj2zChIyIiqmaY1OmoNi+U2H0pBdEHCxdGfDgsGIFeXBhBRERU3TCp01FtXShx7X4m3tx8EQAwtasvBresZ+SIiIiIqCRM6qhUiux8TNtwFll5KnTyd8HcvlwYQUREVF0xqaMSqdUCb/znPG49eIx6jtZYOaYVzM34cSEiIqqueJXWUW2bU7di/zUc/CsNMnMp1oxvDWcujCAiIqrWmNTpqDbNqdt75R4+/+8NAMDSl5ojqJ6DkSMiIiKiZ2FSR1pupGbizZ8LF0ZM6uyDl1rVN3JEREREpAsmdaShzMnHKxvO4lFuATr4OSOqf1Njh0REREQ6YlKnI1OfU6dWC8z66QJupj2Gl4MVVo5pBQsujCAiIqoxeNXWkanPqfv8v9ex/2oqLM2lWD2+NVxtZcYOiYiIiMqBSR3h97j7WLH/OgBg8ZAgBNd3NG5AREREVG5M6mq5hLRHmPXTBQDAxI4NMaKNt3EDIiIiogphUleLZebk45X1Z5CZW4B2Ps7494BAY4dEREREFcSkrpZSqwXe/PkiEtIew8PeCtFjuTCCiIioJuNVXEemtvr1y0M3sC/uPizNChdGuNlxYQQREVFNxqROR6a0+vW/8ffxye/XAAAfDAlCS29H4wZERERElcakrpa59eAxXv/PBQgBjOvQACPbcmEEERGRKWBSV4s8yi0oXBiRU4DWDZ3wzoBmxg6JiIiI9IRJXS0hhMDszRdxPfUR3O1kWDW2FSzN+eMnIiIyFbyq1xKrDidg9+V7sDCTYNW41nC3tzJ2SERERKRHTOpqgUN/pWL53r8AAO8NDkLrhk5GjoiIiIj0jUmdibv98DFe+/E8hABebtcAL7drYOyQiIiIyACY1OmoJu5Tl5VXgGkbzkKZU4DnGzji3UG8YwQREZGpYlKno5q2T50QArNjYhF/LxNudjKsHtcaMnMzY4dFREREBsKkzkR99cdN7IxNgblUglVjW6EuF0YQERGZNCZ1JujP62n4cE88AGDhoGZo4+Ns5IiIiIjI0JjUmZi76VmY+eN5qAUwsk19jGvPhRFERES1AZM6E5Kdp8IrG84iIysfLbwd8d7gIEgkEmOHRURERFWASZ2JEEJg3tZYXE1RwtXWEqvHtYKVBRdGEBER1Ra1Lql7+PAh+vXrhyZNmqB58+YIDw9Hbm6uscOqtLVHbuHXC8kwl0oQPaYVPB2sjR0SERERVaFal9RJJBLMnz8ff/31Fy5evIjs7GysXLnS2GFVyrEbD7Bk11UAwNsDAtHez8XIEREREVFVqzFJ3Y0bNxAREYGWLVvC3NwcQUFBJZaLj49Hr169YGNjAw8PD8yZMwd5eXma552dndGtWzcAgFQqRZs2bXDnzp0qeQ+G8Lc8C5GbzkEtgGGt6mNCx4bGDomIiIiMwNzYAejqypUr2LlzJ9q3bw+1Wg21Wl2sjFwuR2hoKBo1aoStW7ciKSkJs2bNQlZWVom9cTk5OVi3bh2WL19eFW9B73LyVZi24SzkWfloXs8Bi4dyYQQREVFtVWOSuoEDB2Lw4MEAgLCwMJw5c6ZYmdWrV0OpVGLbtm1wdi7cm62goAAzZsxAVFQUvLy8NGXVajUmTpyIHj16oG/fvlXzJvRICIH5Wy/hSrISLjaWWD2+NRdGEBER1WI1JqmTSp89Urx792707NlTk9ABwMiRIxEREYF9+/YhLCxMczwyMhJSqRQrVqwos87c3FythRRKpRIASu0t1Ae1Wg0hRJn1rzuWiG3nk2AmleCLl1vC015msHiIiIioZLpcs/XRhi5qTFKni/j4eISHh2sdc3R0hKenJ+Lj4zXH5syZg7t372Lbtm3PTBaXLl2KRYsWFTuelpaGnJwc/QT+FLVaDYVCASFEifGdvZuJxTuvAQBe61oPfrYqpKamGiQWIiIiKt2zrtn6kJmZqVM5k0rq5HI5HB0dix13cnJCeno6gMK5ecuXL0dAQADatm0LAOjVq1ep8+rmz5+PWbNmaR4rlUp4e3vDzc0N9vb2+n8TKPyASCQSuLm5FfuAJGVk4509l6ASwJCWXni1N+fRERERGUtZ12x9sbLS7f7tJpXU6aJZs2YQQuhcXiaTQSaTITo6GtHR0VCpVAAKh4MN9cMDCrdeebqNnHwVZmw8j4eP89DMyx5LXwqGmRnn0RERERlTSddsfdK13hqzpYkunJycoFAoih2Xy+Va8+wqIjIyEnFxcTh9+nSl6qkoIQQWbLuMS0kKONWxwOpxrWFtyYSOiIiICplUUhcQEKA1dw4AFAoFUlJSEBAQUKm6o6OjERgYqBmyrWobTtzGlnN/QyoBVo5pBW/nOkaJg4iIiKonk0rq+vXrh/379yMjI0NzbPPmzZBKpejdu3el6jZmT92pW+l4b3scAGB+v6bo/JxrlcdARERE1VuNmVOXlZWFXbt2AQBu374NpVKJmJgYAED37t3h5uaGiIgIfPHFFxgyZAiioqKQlJSE2bNnIyIiQmuPupokRZGNGRvPokAtMKiFF6Z09TV2SERERFQN1ZikLjU1FSNGjNA6VvT44MGDCAkJgZOTEw4cOICZM2diyJAhsLOzw5QpU7B48eJKt//0QglDUakFTt58iBt/p6NhhgSf7L+OB4/y0NTTHh8OC+ZKVyIiIiqRRJRnKShBqVTCwcEBCoVC71ua7LmcgkXb45Ci0N7/ro6lGfa+0Y3z6IiIiKoZtVqN1NRUuLu7G2z1q665h0nNqavJ9lxOwfQfzhVL6AAgK0+FK8nFV/USERERFWFSpyNDrn5VqQUWbY9DaV2mEgCLtsdBpWanKhEREZWMSZ2ODLn69dSt9BJ76IoIACmKHJy6la73tomIiMg0MKmrBlIzdbuHrK7liIiIqPZhUqcjQw6/utvpdk83XcsRERFR7cOkTkeGHH5t5+sMTwcrlLZZiQSAp4MV2vlW7lZnREREZLqY1FUDZlIJFg4MBIBiiV3R44UDA2Em5R51REREVDImddVE3yBPrBrXCh4O2kOsHg5WWDWuFfoGeRopMiIiIqoJaswdJWqDvkGe6BXogZM3H+DG32l4rr4b2vu5soeOiIiInolJnY6q6jZhZlIJOvi5wM9WBXd3F0iZ0BEREZEOOPyqI0MulCAiIiKqLCZ1RERERCaASR0RERGRCWBSR0RERGQCmNTpyJB3lCAiIiKqLCZ1OuJCCSIiIqrOmNQRERERmQAmdUREREQmgJsPl5MQAgCgVCoN1oZarUZmZiasrKwglTLvJiIiqq6q4ppdlHMU5SClYVJXTpmZmQAAb29vI0dCREREtUlmZiYcHBxKfV4inpX2kRa1Wo3k5GTY2dmhXbt2Blk4oVQq4e3tjbt378Le3l7v9VPt07ZtWy7yMZDaem5r8vuuzrFXl9iqOo6qaM9QbVTFNVsIgczMTHh5eZXZG8ieunKSSqWoX78+AMDMzMygSZe9vT2TOtILQ39Wa7Paem5r8vuuzrFXl9iqOo6qaK+mX7PL6qErwglblRAZGWnsEIh0ws+q4dTWc1uT33d1jr26xFbVcVRFe9Xl3BoSh1+rIaVSCQcHBygUimrxFxsRERGVrDpds9lTVw3JZDIsXLgQMpnM2KEQERFRGarTNZs9dUREREQmgD11RERERCaASR0RERGRCWBSR0RERGQCmNQRERERmQAmdTXY1KlTIZFIjB0GERERlcLHxweBgYFo2bIlWrZsibi4OIO1xTtK1FD79+9Hfn6+scMgIiKiZ9i1axd8fHwM3g576qrIjRs3EBERgZYtW8Lc3BxBQUEllouPj0evXr1gY2MDDw8PzJkzB3l5eVplHj9+jAULFuCTTz6pitCJiIhqFX1es6sSe+qqyJUrV7Bz5060b98earUaarW6WBm5XI7Q0FA0atQIW7duRVJSEmbNmoWsrCysXLlSU27+/Pl4/fXX4eLiUpVvgYiIqFbQ5zUbAIYOHQq1Wo3+/fvjvffeg4WFhWECF1QlVCqV5v8nTpwomjVrVqzMkiVLhI2NjXj48KHm2Jo1a4SZmZlISkoSQghx5MgR0b9/f83z/BESERHpl76u2UIIcffuXSGEEJmZmWLIkCHigw8+MFjcHH6tIlLps0/17t270bNnTzg7O2uOjRw5Emq1Gvv27QMA/Pnnnzh//jx8fHw04/M+Pj64d++eQeImIiKqbfR1zQaA+vXrAwBsbW0xefJknDhxQv8B/w+TumokPj4eAQEBWsccHR3h6emJ+Ph4AMC8efOQnJyMxMREJCYmAgASExPh4eFR1eESERHVWrpcsx8/fgylUgkAKCgowJYtWxAcHGywmDinrhqRy+VwdHQsdtzJyQnp6elVHxARERGVSJdr9v379/HSSy9BrVajoKAAnTp1woIFCwwWE5O6Gk4IYewQiIiIqAR+fn64cOFClbXH4ddqxMnJCQqFothxuVyuNWZPRERExlUdr9lM6qqRgIAAzTh8EYVCgZSUlGLj9kRERGQ81fGazaSuGunXrx/279+PjIwMzbHNmzdDKpWid+/exguMiIiItFTHa7ZEcFJWlcjKysKuXbsAANHR0UhISMCnn34KAOjevTvc3Nwgl8vRrFkzNG7cGFFRUZqNDMeOHVtsI0MiIiIyjJp6zWZSV0USExPh6+tb4nMHDx5ESEgIAODq1auYOXMmjh07Bjs7O0yYMAGLFy+GpaVlFUZLRERUe9XUazaTOiIiIiITwDl1RERERCaASR0RERGRCWBSR0RERGQCmNQRERERmQAmdUREREQmgEkdERERkQlgUkdERERkApjUEREREZkAJnVEREREJoBJHREREZEJYFJHREREZAKY1BHVMGFhYQgKCqp0GUO1XZXxFMnIyIBEIsG6desqXdcvv/yCL7/8svJB1WDV6RycO3cOEonkmf9Onz5d7LUhISGa5z/44APN8XfffRe2trbFyr/55puQSqVYu3YtAGDevHma1w8YMMBwb5JIT8yNHQARUXXyyy+/4MyZM5gxY4axQzGa6nQOfHx8cPz4cc3j77//HmvWrMGhQ4dgaWmpOf7888+X+PrOnTvj448/hre3d5ntzJ07F5999hlWr16NyZMnAwAiIyMxZMiQanEeiHTBpI6IyMQJIZCXlweZTFbj2nd2dkaHDh00j7/88kv4+vqiW7duOr3e0dFR6/Ul+fe//42PPvoIX375JV555RXNcW9vb3h7e8Pe3r7ccRMZA4dfiWqo3bt3IygoCFZWVmjdujVOnDhRatnjx49j0KBB8PLygo2NDVq2bIkNGzaUWrZ3796wt7eHnZ0d2rdvj99//73Esmq1GlOmTIGrqyvOnDmjl5i3bt2Kli1bwsrKCl5eXpg1axZycnK0ynz99dfw8fFBnTp18MILL+DGjRtaz2/fvh0SiQTXr1/XOi6Xy2FtbV3q0GJYWBi+//57XLlyRTPsFhYWpnVuQkNDYWNjAwcHB4wZMwapqalarw8KCsL+/fsRHBwMa2trdO/eHYmJiUhPT8fIkSNhb28Pf39//PTTT8XaDgoK0ukc6RrHrl270KJFC8hkMmzfvl3z2rI+C2Wdg5CQkGLDkBcuXIBEIsGhQ4d0br+s2J/l0qVLaN68uc7ln+Xdd9/F4sWL8cUXX2D69Ol6q5fIGJjUEdVAKSkpmDFjBmbPno2ff/4ZMpkMffr0KfXiePv2bXTu3BnffPMNtm/fjmHDhmHy5Mn4/vvvtcodPXoUISEhyM3NxTfffIMtW7Zg8ODBuHPnTrE6CwoKMHbsWOzcuROHDh1CmzZtKh3zb7/9huHDhyMwMBC//PIL5syZg9WrV2PcuHGaMjt27MArr7yCHj16YNu2bXjhhRcwYsQIrbb69++PevXq4dtvv9U6vmnTJgDAmDFjSozx7bffRv/+/eHn54fjx4/j+PHjePvttwEUJiMhISFwcHDATz/9hK+++gqnT5/G4MGDteq4d+8e3nzzTSxYsAAbN25EQkICxo4di1GjRqF58+bYsmULWrdujXHjxuH27dvlPke6xpGcnIzXXnsN//rXv7Bnzx60bNkSwLM/C2Wdg/IoqX1dYy9NQUEBrl69qrekbvHixVi0aBE+++wzvPrqq3qpk8ioBBHVKBMnThQAxIEDBzTHMjIyhJ2dnZg3b56mTLNmzUp8vVqtFvn5+eKVV14RHTt21HquU6dOIjAwUBQUFJTadrNmzUROTo4YNGiQaNCggbh27ZpeYhZCiOeff75YTGvWrBEARGxsrBBCiPbt24uuXbtqlXn77bcFAPHdd99pjv373/8WXl5eWu+lVatWYsyYMc+MtaRz161bN9GpUyehVqs1x65cuSIkEonYuXOn5rUSiURcvnxZU+aLL74QAMTcuXM1x+RyuTAzMxMrVqwo9znSNQ4A4sSJE2W+19I+C6Wdg+7du4sXX3xR69j58+cFAHHw4MFi7+Xp9nWJvSxXrlwRAMTPP//8zLKlxSuEEAsXLhQABAAxZcqUCtdDVN2wp46oBnJwcEBoaKjW4549e+LkyZMllpfL5XjttdfQsGFDWFhYwMLCAl999RWuXbumKZOVlYUTJ05g4sSJMDMzK7Xt7OxsDBgwAFevXsWff/6JRo0aaZ4TQqCgoEDzT6VS6Rzzo0ePcOHCBQwfPlyrvVGjRgEAjhw5ApVKhbNnz2Lo0KFaZZ5+DQBMnjwZKSkp2LNnDwAgNjYW586d00yCL4+srCwcPXoUI0aMgEql0ry/xo0bw9vbW2vlpZeXF5o1a6Z53LhxYwBAz549NcccHR3h7u6Ou3fvarXzrHNUnjhcXFzQvn37Yu9Fl8+CPjzdfnliL82lS5cAQC89ddbW1ujWrRs2bdqEo0ePVro+ouqASR1RDeTm5lbsWN26dZGSklJi+bCwMPz444946623sG/fPpw+fRrh4eFac9XkcjnUajW8vLzKbDstLQ2HDx/Giy++iAYNGmg9d/jwYU2iYGFhgRdeeEHnmDMyMiCEQN26dbXKODg4QCaTIT09HWlpaSgoKIC7u3uxep7m4+ODXr16aban+Pbbb+Hr64sePXqU+f5KIpfLoVKp8K9//Uvr/VlYWODOnTtayZmjo6PWa4tWaJZ0/Om5gs86R+WJo6RzAuj2WdCHp9svT+yluXTpEqysrLT+kKgoqVSK3377DY0bN8aAAQM0CSNRTcbVr0Q1UFpaWrFj9+/fh6enZ7HjOTk52LFjBz799FPMnDlTc1ytVmuVc3R0hFQqRXJycpltN2jQAO+++y5Gjx4NV1dXLFiwQPNc69attXpc7OzsdI7Z0dEREomk2LxAhUKB3NxcODs7w83NDebm5sXK3L9/v8RYp06dijFjxiApKQkbN27Ea6+9BolEUub7K0lRbFFRURgyZEix511dXctdZ0l0PUe6xFHS+9T1s1AaKysr5OXlaR2Ty+Ulln26fX2cw9jYWDRt2rTMnuTycHBwwN69e9G5c2f06dMHR48eha+vr17qJjIG9tQR1UAKhQL//e9/tR7v37+/xOG23NxcqNVqrT29MjMz8dtvv2mVs7GxQceOHbF+/XqtYdOSDB8+HN9//z3eeecdrFixQnPczs4Obdq00fxr0qSJzjHb2tqiZcuWiImJ0Wrr559/BgB06dIFZmZmaNWqFbZt26ZV5unXFBk8eDCcnJwwZswYpKena61kLU1JPWhF5+bq1ata76/on4+PzzPr1cWzzlFl49D1s1DSOQCA+vXr46+//oIQQnNs3759Or03fZzDS5cuITg4WKf2dOXu7o7ff/8dEokEvXr1wr179/RaP1FVYk8dUQ3k7OyMyZMnY9GiRXB0dMSyZcsghMAbb7xRrKyDgwPatm2LZcuWaXq6li1bBgcHh2I9XsuWLUNoaCh69uyJGTNmwMnJCefOnYOrqyvCw8O1yo4dOxbZ2dmYNm0arK2tMW3atErH/O6772LIkCEYN24cxo0bh7/++gtRUVEYNmyYZh7VggULMHjwYEyaNAmjR4/G2bNnS92excLCAhMnTsTy5cvRp0+fZ25ACwBNmzbFt99+ix9//BGNGjWCq6srfHx8sHz5coSGhmLUqFEYPXo0nJyc8Pfff+P333/HpEmTEBIS8sy6n0WXc1SZOHT9LJR2DoYPH461a9di5syZGDJkCI4dO1ZqQl2SysSemZmJ27dv63U7kyI+Pj7Yu3cvunXrhr59++Lw4cNwcHDQeztEBmfMVRpEVH5FKxN37NghmjZtKiwtLcXzzz8vjh49WqxMkevXr4vQ0FBRp04d4e3tLZYvXy4WLlwobGxsitV/9OhR0aNHD1GnTh1hZ2cnOnToIPbv319ivUIUru6USqXi+++/r1TMRWJiYkRwcLCwtLQUHh4e4o033hDZ2dlaZVavXi28vb2FlZWV6N69uzh58mSx1a9Fjh07JgCIn376qdT4nqRQKMTo0aOFi4uLACAmTpyoee706dOif//+wsHBQVhbW4tGjRqJiIgIcffu3VLPz8GDBwUAcfr0aa3jDRs2FJGRkRU6RxWJo4gun4WyzsFHH30kvL29hY2NjRgxYoTYv39/iatfS2v/WbGXpujnuHfv3jLLPams1a8lffaPHTsm6tSpI7p06SKysrKeWQ9RdSMR4ol+dCIiE/POO+/gyy+/RFJSktHuqKCLsLAwnDlzBpcvXzZ2KCYjJCQENjY2+PXXX2FmZlbu+ZRqtRpqtRovvPAC7OzssGPHDgNFSqQfnFNHRCbpr7/+wm+//YYvvvgCERER1TqhI8PZtWsXLCwssHjx4nK/NioqChYWFvjjjz8MEBmR/rGnjohMUkhICE6cOIG+ffti48aNsLGxMXZIZWJPnf799ddfyMzMBADUq1evxNXhZUlKStJsJ+Po6IjnnntO7zES6ROTOiIiIiITwOFXIiIiIhPApI6IiIjIBDCpIyIiIjIBTOqIiIiITACTOiIiIiITwKSOiIiIyAQwqSMiIiIyAUzqiIiIiEwAkzoiIiIiE8CkjoiIiMgE/D/x455L35IQUgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 700x450 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rate(1e5 K)/rate(1e4 K) = 1.209e+07\n"
     ]
    }
   ],
   "source": [
    "# A hydrogen-like toy cross-section — zero below the ionization edge,\n",
    "# falling as (lambda/912A)^3 above threshold energy — integrated against\n",
    "# black bodies of increasing temperature. Hotter stars are dramatically\n",
    "# more ionizing; the steep rise mirrors the Wien tail crossing 13.6 eV.\n",
    "def crossSectionHydrogenLike(wavelength):\n",
    "    return 6.3e-18*(wavelength/912.0)**3 if wavelength <= 912.0 else 0.0\n",
    "\n",
    "temperatures = np.array([1.0e4, 1.5e4, 2.0e4, 3.0e4, 4.5e4, 7.0e4, 1.0e5])\n",
    "rates = np.array([galacticus.radiationFieldBlackBody(T)\n",
    "                  .integrateOverCrossSection(np.array([50.0, 912.0]),\n",
    "                                             crossSectionHydrogenLike, nodeNull)\n",
    "                  for T in temperatures])\n",
    "plt.loglog(temperatures, rates, 'o-')\n",
    "plt.xlabel('black-body temperature $T$ [K]')\n",
    "plt.ylabel('$\\\\int F_\\\\lambda\\\\, \\\\sigma(\\\\lambda)\\\\, \\\\mathrm{d}\\\\lambda$  [arbitrary]')\n",
    "plt.title('Photo-ionization integral of a Python cross-section')\n",
    "plt.show()\n",
    "print(f\"rate(1e5 K)/rate(1e4 K) = {rates[-1]/rates[0]:.3e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfd728aa",
   "metadata": {},
   "source": [
    "## Notes on the callback contract\n",
    "\n",
    "* The Python callable is invoked synchronously from inside the Fortran\n",
    "  quadrature; keep it fast (vectorization happens per-call, not across\n",
    "  calls).\n",
    "* Exceptions raised inside a callback cannot propagate through the\n",
    "  Fortran stack — ctypes prints the traceback and returns garbage.\n",
    "  Validate inputs *before* integrating.\n",
    "* The callback slot is per-method global state: don't call the same\n",
    "  method concurrently from multiple threads.\n",
    "\n",
    "Combined with the earlier tutorials, this closes the loop: Galacticus\n",
    "objects in your Python analysis, and your Python physics inside\n",
    "Galacticus's numerics."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.14.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}