{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "38bf69d0",
      "metadata": {},
      "source": [
        "# Spontaneous emission of TLS: TCP socket\n",
        "\n",
        "Here, we introduce the two-level system (TLS) spontaneous emission tutorial using `maxwelllink.Molecule` and a TCP socket."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f71d059c",
      "metadata": {},
      "source": [
        "## 1. Setting up the socket communication layer\n",
        "\n",
        "Using the TCP socket requires to set the `hostname` and `port number`. In local machines, we can use the helper function `get_available_host_port()` from **MaxwellLink** to obtain these two pieces of information. Then, we initialize a `SocketHub` instance to provide the socket communication in **MaxwellLink**."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "2c8428dc",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Using MPI version 4.1, 1 processes\n",
            "SocketHub listening on 127.0.0.1:50692\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "import maxwelllink as mxl\n",
        "from maxwelllink import sockets as mxs\n",
        "\n",
        "try:\n",
        "    import meep as mp\n",
        "except ImportError as exc:  \n",
        "    raise RuntimeError(\n",
        "        \"Meep is required for this tutorial.\"\n",
        "        \"Install via conda: conda install -c conda-forge pymeep=*=mpi_mpich_*\"\n",
        "    ) from exc\n",
        "\n",
        "host, port = mxs.get_available_host_port()\n",
        "hub = mxl.SocketHub(host=host, port=port, timeout=10.0, latency=1e-5)\n",
        "\n",
        "print(f\"SocketHub listening on {host}:{port}\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3f43f478",
      "metadata": {},
      "source": [
        "## 2. Bind Molecule and EM solver to the SocketHub\n",
        "\n",
        "Then, we create a `Molecule` instance to define the information of this molecule in the EM simulation environment, including the `center`, `size`, `sigma` (width of the molecular polarization distribution), and `dimensions`. \n",
        "\n",
        "We also need to setup the EM solver (**MEEP**) using `mxl.MeepSimulation`. This class is a wrapper of the `meep.Simulation` object with extended parameters for **MaxwellLink**."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "a5f84aa1",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[Init Molecule] Under socket mode, registered molecule with ID 0\n",
            "\n",
            "\n",
            " ######### MaxwellLink Units Helper #########\n",
            " MEEP uses its own units system, which is based on the speed of light in vacuum (c=1), \n",
            " the permittivity of free space (epsilon_0=1), and the permeability of free space (mu_0=1). \n",
            " To couple MEEP with molecular dynamics, we set [c] = [epsilon_0] = [mu_0] = [hbar] = 1. \n",
            " By further defining the time unit as 1.00E-01 fs, we can fix the units system of MEEP (mu).\n",
            "\n",
            " Given the simulation resolution = 10,\n",
            " - FDTD dt = 5.00E-02 mu (0.5/resolution) = 5.00E-03 fs\n",
            " - FDTD dx = 1.00E-01 mu (1.0/resolution) = 3.00E+00 nm\n",
            " - Time [t]: 1 mu = 1.00E-01 fs = 4.13E+00 a.u.\n",
            " - Length [x]: 1 mu = 3.00E+01 nm\n",
            " - Angular frequency [omega = 2 pi * f]: 1 mu = 6.5851E+00 eV = 2.4200E-01 a.u.\n",
            " - Electric field [E]: 1 mu = 6.65E+07 V/m = 1.29E-04 a.u.\n",
            " Hope this helps!\n",
            " ############################################\n",
            "\n",
            "\n"
          ]
        }
      ],
      "source": [
        "molecule = mxl.Molecule(\n",
        "    hub=hub,\n",
        "    center=mp.Vector3(0, 0, 0),\n",
        "    size=mp.Vector3(1, 1, 1),\n",
        "    sigma=0.1,\n",
        "    dimensions=2,\n",
        ")\n",
        "\n",
        "sim = mxl.MeepSimulation(\n",
        "    hub=hub,\n",
        "    molecules=[molecule],\n",
        "    cell_size=mp.Vector3(8, 8, 0),\n",
        "    boundary_layers=[mp.PML(3.0)],\n",
        "    resolution=10,\n",
        "    # fix a units system \n",
        "    time_units_fs=0.1,\n",
        ")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "20cf3a66",
      "metadata": {},
      "source": [
        "## 3. Python way to lunch `mxl_driver` on a separate terminal\n",
        "\n",
        "Generally, using the Socket Interface requires to launch the EM simulation in one terminal and then start the molecular driver simulation in a separate terminal. To avoid openning a second terminal, below we introduce a python helper function `launch_tls_driver(...)`, which will launch `mxl_driver` from Python (so we can stay within this notebook to finish this tutorial). \n",
        "\n",
        "Here, we set the TLS starting at the initial excited-state population of 1e-4.\n",
        "\n",
        "Immediately after launching this driver in the background, we run the simulation using `sim.run(...)`. This function is a wrapper of the `meep.Simulation.run(...)` function, which can accept user-defined step functions."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "d2c4f2ad",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Launching TLS driver via subprocess...\n",
            "If you prefer to run it manually, execute:\n",
            "  /Users/taoli/miniforge3/envs/mxl/bin/mxl_driver --model tls --address 127.0.0.1 --port 50692 --param \"omega=0.242, mu12=187, orientation=2, pe_initial=1e-4\"\n",
            "-----------\n",
            "Initializing structure...\n",
            "time for choose_chunkdivision = 0.000397 s\n",
            "Working in 2D dimensions.\n",
            "Computational cell is 8 x 8 x 0 with resolution 10\n",
            "time for set_epsilon = 0.003466 s\n",
            "-----------\n",
            "[SocketHub][initialization] Time step in atomic units: 0.20670686667500004\n",
            "[initialization] Assigned a molecular ID: 0\n",
            "init TLSModel with dt = 0.206707 a.u., molecule ID = 0\n",
            "[initialization] Finished initialization for molecular ID: 0\n",
            " CONNECTED: mol 0 <- 127.0.0.1:50693\n",
            "Meep progress: 284.35/400.0 = 71.1% done in 4.0s, 1.6s to go\n",
            "on time step 5708 (time=285.4), 0.000700828 s/step\n",
            "run 0 finished at t = 400.0 (8000 timesteps)\n",
            "[SocketHub] DISCONNECTED: mol 0 from 127.0.0.1:50693\n",
            "Received STOP, exiting\n"
          ]
        }
      ],
      "source": [
        "import shlex\n",
        "import shutil\n",
        "import subprocess\n",
        "import time\n",
        "\n",
        "\n",
        "def launch_tls_driver(host: str, port: int, sleep_time: float = 0.5):\n",
        "    executable = shutil.which('mxl_driver')\n",
        "    if executable is None:\n",
        "        raise RuntimeError('mxl_driver executable not found in PATH.')\n",
        "    cmd = (\n",
        "        f\"{executable} --model tls --address {host} --port {port}\"\n",
        "        f' --param \"omega=0.242, mu12=187, orientation=2, pe_initial=1e-4\"'\n",
        "    )\n",
        "    print('Launching TLS driver via subprocess...')\n",
        "    print('If you prefer to run it manually, execute:')\n",
        "    print('  ' + cmd)\n",
        "    argv = shlex.split(cmd)\n",
        "    proc = subprocess.Popen(argv)\n",
        "    time.sleep(sleep_time)\n",
        "    return proc\n",
        "\n",
        "launch_tls_driver(host, port)\n",
        "\n",
        "sim.run(until=400)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "cdabe7d1",
      "metadata": {},
      "source": [
        "## 4. Retrieve molecular simulation data\n",
        "\n",
        "After the simulation, we can retrieve molecular simulation data from `molecule.additional_data_history`, a Python list which stores the molecular information sent from the driver code at each step of the simulation."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "9120e20d",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Collected 8001 samples.\n"
          ]
        }
      ],
      "source": [
        "population = np.array([entry[\"Pe\"] for entry in molecule.additional_data_history])\n",
        "time_au = np.array([entry[\"time_au\"] for entry in molecule.additional_data_history])\n",
        "\n",
        "print(f\"Collected {population.size} samples.\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ae529503",
      "metadata": {},
      "source": [
        "## 5. Compare with the Analytical Golden-Rule Decay\n",
        "\n",
        "Finally, we can compare this numerical simulation with analytical golden-rule rate calculations, with the 2D  spontaneus emission in vacuum as:\n",
        "\n",
        "$$\n",
        "k_{\\rm {FGR}} = \\frac{|\\mu_{12}|^2 \\omega^2}{2\\hbar\\epsilon_0 c^2}\n",
        "$$\n",
        "\n",
        "The corresponding TLS excited-state population decay dynamics obey:\n",
        "\n",
        "$$\n",
        "P_{\\rm e}^{\\rm {QM}}(t) = P_{\\rm e}(0) e^{-k_{\\rm{FGR}} t}\n",
        "$$\n",
        "\n",
        "When the EM field is described entirely classically, more than half a century ago, Jaynes and collaborators (https://ieeexplore.ieee.org/document/1443594) calculated the semiclassical spontaneous emission rate. More recently, we have also reproduced this semiclassical excited-state population decay (https://doi.org/10.1103/PhysRevA.97.032105):\n",
        "\n",
        "$$\n",
        "P_{\\rm e}^{\\rm {sc}}(t) = \\frac{e^{-k_{\\rm{FGR}} t}}{e^{-k_{\\rm{FGR}} t} + \\frac{1 - P_{\\rm e}(0) }{P_{\\rm e}(0) }}\n",
        "$$\n",
        "\n",
        "When $P_{\\rm e}(0)\\rightarrow 0$, the semiclassical decay dynamics exactly agree with the quantum correspondance $P_{\\rm e}^{\\rm {QM}}(t)$.\n",
        "\n",
        "As shown below, using $P_{\\rm e}(0)= 10^{-4}$, our semiclassical simulation exactly reproduces the quantum golden-rule decay."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "ebe9c2ed",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "std_dev=2.012e-03, max_abs_diff=7.928e-03\n"
          ]
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAGGCAYAAACNCg6xAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAcyxJREFUeJzt3Xl8TNf7wPHPrNkTCVkJScTS2CVF0ga1JCil1qoqVS3dNNVWVbV0+draoqqqWmppv6i1aqvYgoo9dqUlxJKIIAmyz9zfH/mZb9MEMyQmy/N+veYlc+655z5nLubJueeeq1IURUEIIYQQQtyT2toBCCGEEEKUFZI4CSGEEEKYSRInIYQQQggzSeIkhBBCCGEmSZyEEEIIIcwkiZMQQgghhJkkcRJCCCGEMJMkTkIIIYQQZtJaO4CKzmg0cunSJZycnFCpVNYORwghhKhwFEXhxo0b+Pj4oFbffUxJEicru3TpEr6+vtYOQwghhKjwzp8/T7Vq1e5aRxInK3NycgLyT5azs7OVoxFCCCEqnvT0dHx9fU3fyXcjiZOV3b485+zsLImTEEIIYUXmTJmRyeFCCCGEEGaSxEkIIYQQwkySOAkhhBBCmEnmOAkhhKjwDAYDubm51g5DlBCdTodGoymWtiRxEkIIUWEpikJSUhKpqanWDkWUsEqVKuHl5fXAayZK4iSEEKLCup00eXh4YG9vLwsRl0OKopCRkUFycjIA3t7eD9SeJE5CCCEqJIPBYEqaKleubO1wRAmys7MDIDk5GQ8Pjwe6bGf1yeEzZszA398fW1tbgoOD2b59+13rx8TEEBwcjK2tLQEBAcycObNQnWXLlhEUFISNjQ1BQUGsWLHC4uMuX76cyMhIqlSpgkql4uDBg4XayM7O5o033qBKlSo4ODjw1FNPceHCBcs+ACGEEFZxe06Tvb29lSMRD8Pt8/ygc9msmjgtXryYqKgoPvjgA+Li4ggPD6djx44kJCQUWT8+Pp5OnToRHh5OXFwco0aNYtiwYSxbtsxUJzY2lj59+tC/f38OHTpE//796d27N7t377bouLdu3eKxxx5jwoQJd4w/KiqKFStWsGjRInbs2MHNmzfp3LkzBoOhGD4dIYQQD4NcnqsYiu08K1bUrFkzZejQoQXK6tatq4wcObLI+iNGjFDq1q1boGzIkCFKixYtTO979+6tdOjQoUCdyMhI5Zlnnrmv48bHxyuAEhcXV6A8NTVV0el0yqJFi0xlFy9eVNRqtbJ+/foi4y9KWlqaAihpaWlm7yOEEOLBZWZmKsePH1cyMzOtHYp4CO52vi35LrbaiFNOTg779+8nIiKiQHlERAQ7d+4scp/Y2NhC9SMjI9m3b59p6O1OdW63eT/HLcr+/fvJzc0t0I6Pjw/169e3qJ2SsunEZdIy5dZaIYSoiFQqFStXrizx4/j5+TF16tRS087DYLXEKSUlBYPBgKenZ4FyT09PkpKSitwnKSmpyPp5eXmkpKTctc7tNu/nuHeKRa/X4+rqalE72dnZpKenF3gVt/PXMhgybxdfT3qf00nXi719IYQQ1pWcnMyQIUOoXr06NjY2eHl5ERkZSWxsLACJiYl07NjRylEWNnfuXCpVqlSofO/evbz88ssPP6D7YPW76v59zVFRlLtehyyq/r/LzWnT0uOa617tjB8/no8//viBj3M3yelZ/Nd2Es2Uo6xfqabm0EklejwhhBAPV48ePcjNzWXevHkEBARw+fJlNm3axLVr1wDw8vKycoSWcXd3t3YIZrPaiFOVKlXQaDSFRmeSk5MLjQbd5uXlVWR9rVZrupX0TnVut3k/x71TLDk5OVy/XnBE517tvP/++6SlpZle58+fN/uY5gr2c8M1bAAATyTNIefyn8V+DCGEENaRmprKjh07mDhxIk888QQ1atSgWbNmvP/++zz55JNAwUt1Z8+eRaVS8csvvxAeHo6dnR2PPvoop06dYu/evYSEhODo6EiHDh24cuWK6TitW7cmKiqqwLG7devGwIED7xjb5MmTadCgAQ4ODvj6+vLqq69y8+ZNALZu3coLL7xAWloaKpUKlUrF2LFjgcKX6hISEujatSuOjo44OzvTu3dvLl++bNo+duxYGjduzIIFC/Dz88PFxYVnnnmGGzdu3P8HayarJU56vZ7g4GCio6MLlEdHRxMWFlbkPqGhoYXqb9iwgZCQEHQ63V3r3G7zfo5blODgYHQ6XYF2EhMTOXr06F3bsbGxwdnZucCrJAS0HcxOVRNsyOXWkldIy8gukeMIIUR5oigKGTl5D/11++qJORwdHXF0dGTlypVkZ5v/f/uYMWMYPXo0Bw4cQKvV0rdvX0aMGMFXX33F9u3bOX36NB999NH9fGwmarWaadOmcfToUebNm8fmzZsZMWIEAGFhYUydOhVnZ2cSExNJTEzknXfeKdSGoih069aNa9euERMTQ3R0NKdPn6ZPnz4F6p0+fZqVK1eyevVqVq9eTUxMzF3vhC8uVr1UN3z4cPr3709ISAihoaHMmjWLhIQEhg4dCuSPzly8eJH58+cDMHToUKZPn87w4cN56aWXiI2NZfbs2SxcuNDU5ptvvknLli2ZOHEiXbt25ddff2Xjxo3s2LHD7OMCXLt2jYSEBC5dugTAyZMngfyRJi8vL1xcXHjxxRd5++23qVy5Mm5ubrzzzjs0aNCAdu3alfhndy8ajZrd9T6k4ZFncU05wNhxI2jw9Lv0CK5m7dCEEKLUysw1EPTR7w/9uMc/icReb95XslarZe7cubz00kvMnDmTpk2b0qpVK5555hkaNmx4x/3eeecdIiMjgfzvyr59+7Jp0yYee+wxAF588UXmzp37QP345wiVv78/n376Ka+88gozZsxAr9fj4uKCSqW666XEjRs3cvjwYeLj4/H19QVgwYIF1KtXj7179/Loo48CYDQamTt3Lk5OTgD079+fTZs28Z///OeB+nAvVl3HqU+fPkydOpVPPvmExo0bs23bNtauXUuNGjWA/BGcf66t5O/vz9q1a9m6dSuNGzfm008/Zdq0afTo0cNUJywsjEWLFvHjjz/SsGFD5s6dy+LFi2nevLnZxwVYtWoVTZo0MQ17PvPMMzRp0qTAgptTpkyhW7du9O7dm8ceewx7e3t+++23YnuQ4IN6onkwE/L6AvCuZiE//b4dg9H832qEEEKUTj169ODSpUusWrWKyMhItm7dStOmTe+a+Pwzqbo9paRBgwYFym4/luR+bdmyhfbt21O1alWcnJx4/vnnuXr1Krdu3TK7jRMnTuDr62tKmgCCgoKoVKkSJ06cMJX5+fmZkibIf5TKg8ZvDqtPDn/11Vd59dVXi9xW1F+AVq1aceDAgbu22bNnT3r27HnfxwUYOHDgXa/jAtja2vL111/z9ddf37WetTT2rcTpp4cTv+UA/rcO8WrmLE4kdqJ+VRdrhyaEEKWSnU7D8U8irXJcS9na2tK+fXvat2/PRx99xODBgxkzZswdv7tuT2mB/90g9e8yo9Foeq9WqwtdQrzbqtvnzp2jU6dODB06lE8//RQ3Nzd27NjBiy++aNFq3Xe6yerf5f+Mvaj4S4rVH7kiSlaPkOr4vzCHY7bBjMvrxx9/p1g7JCGEKLVUKhX2eu1DfxXHXd1BQUEWjezci7u7O4mJiab3BoOBo0eP3rH+vn37yMvL48svv6RFixbUrl3bNN3lNr1ef8+nawQFBZGQkFDg5qnjx4+TlpbGI488cp+9KT6SOFUEVQKJfewH4hVv/jh91drRCCGEeABXr16lTZs2/PTTT6a5QEuWLGHSpEl07dq12I7Tpk0b1qxZw5o1a/jzzz959dVXSU1NvWP9mjVrkpeXx9dff82ZM2dYsGBBoefJ+vn5cfPmTTZt2kRKSgoZGRmF2mnXrh0NGzakX79+HDhwgD179vD888/TqlUrQkJCiq1/90sSpwqiVe38NTK2nbrCu98uYd7Os9YNSAghxH1xdHSkefPmTJkyhZYtW1K/fn0+/PBDXnrpJaZPn15sxxk0aBADBgwwJS3+/v488cQTd6zfuHFjJk+ezMSJE6lfvz4///wz48ePL1AnLCyMoUOH0qdPH9zd3Zk0qfA6g7eXUnB1daVly5a0a9eOgIAAFi9eXGx9exAqxZJ7IEWxS09Px8XFhbS0tBJbmuC2DlO30S3lO17WrOHV3DcZMjSKJtVd772jEEKUQ1lZWcTHx+Pv74+tra21wxEl7G7n25LvYhlxqkA+6Vofezs71CqFT3U/8vveE/feSQghhBAmkjhVIM383Xj+vW+46VwTd1UajY9NsGjRNSGEEKKik8SpotHaoOs+A4OiooMxhguxS6wdkRBCCFFmSOJUAdn4tWBDpfyl6x2j32XKr7Fcu5Vj5aiEEEKI0k8SpwrqVugIThqr4aqkUmvfWD757Zi1QxJCCCFKPUmcKqgnm/oxu8q7XFMc2WJozNojiWTk5Fk7LCGEEKJUk8SpgrLTa5g0bCCuo06yyzmSHIPC7vhr1g5LCCGEKNUkcargVDaOtKxdBYBdJxJA7rITQggh7kgSJ0HLWu60Uh9iUFxPNi/5hr1nZeRJCCFE/iNSpk6d+kBtbN26FZVKddfHtVji7NmzqFQqDh48WCztWUoSJ0HrOh48qovHU3WdpsfG8fbsdSSnZ1k7LCGEEPewc+dONBoNHTp0sHYoALRu3ZqoqKgCZWFhYSQmJuLi4mKdoIqZJE4CO70Gl4gRHFUCqKS6xRi+Z2XcBWuHJYQQ4h7mzJnDG2+8wY4dO0hISLB2OEXS6/V4eXmhUqmsHUqxkMRJAND/sVrUGbIAg0pHW00cysH/WjskIYQQd3Hr1i1++eUXXnnlFTp37szcuXNN225fHtu0aRMhISHY29sTFhbGyZMnTXVOnz5N165d8fT0xNHRkUcffZSNGzfe8XiDBg2ic+fOBcry8vLw8vJizpw5DBw4kJiYGL766itUKhUqlYqzZ88Weanujz/+oFWrVtjb2+Pq6kpkZCTXr18HYP369Tz++ONUqlSJypUr07lzZ06fPl08H1oxkMRJmOh86pPa/B0Anr02g8wrZ60bkBBCiDtavHgxderUoU6dOjz33HP8+OOPhR6j9cEHH/Dll1+yb98+tFotgwYNMm27efMmnTp1YuPGjcTFxREZGUmXLl3uOHI1ePBg1q9fT2Jioqls7dq13Lx5k969e/PVV18RGhrKSy+9RGJiIomJifj6+hZq5+DBg7Rt25Z69eoRGxvLjh076NKlCwaDAchPCIcPH87evXvZtGkTarWap59+GqPRWBwf2wPTWjsAUbq4tX+bo3uWUd94ipNzBhAdMosXwwOx02usHZoQQjw8ObfuvE2lAZ2tmXXVoLO7e129g+XxAbNnz+a5554DoEOHDty8eZNNmzbRrl07U53//Oc/tGrVCoCRI0fy5JNPkpWVha2tLY0aNaJRo0amup999hkrVqxg1apVvP7664WOFxYWRp06dViwYAEjRowA4Mcff6RXr144Ojrmd0Wvx97eHi8vrzvGPWnSJEJCQpgxY4aprF69eqafe/ToUaifHh4eHD9+nPr165v9+ZQUSZxEASqNjt2NxhFwoB/HbzowPfoY6dlGRnV6xNqhCSHEwzPO587bakVAv3885/PzQMjNKLpujcfhhTX/ez+1AWRcLVhnbJrF4Z08eZI9e/awfPlyALRaLX369GHOnDkFEqeGDRuafvb29gYgOTmZ6tWrc+vWLT7++GNWr17NpUuXyMvLIzMz865zpQYPHsysWbMYMWIEycnJrFmzhk2bNlkU+8GDB+nVq9cdt58+fZoPP/yQXbt2kZKSYhppSkhIkMRJlE5d27bkzbMziE7M/y1pyb7zvNehLhp1+ZjYJ4QQZd3s2bPJy8ujatWqpjJFUdDpdKa5QgA6nc708+3J2bcTkXfffZfff/+dL774gsDAQOzs7OjZsyc5OXd+dunzzz/PyJEjiY2NJTY2Fj8/P8LDwy2K3c7O7q7bu3Tpgq+vL99//z0+Pj4YjUbq169/17geJkmcRCFVHG34/s2e5BqMNP00musZORw5f5XGNapYOzQhhHg4Rl268zbVv6YuvPv3Xer+aypx1JH7j+n/5eXlMX/+fL788ksiIiIKbOvRowc///yzWSMz27dvZ+DAgTz99NNA/pyns2fP3nWfypUr061bN3788UdiY2N54YUXCmzX6/WmuUp30rBhQzZt2sTHH39caNvVq1c5ceIE3333nSkh27Fjxz378jBJ4iTuSKdR09FPQ/vTE8nbUA9e+traIQkhxMNhybyjkqp7B6tXr+b69eu8+OKLhdZG6tmzJ7Nnz2bKlCn3bCcwMJDly5fTpUsXVCoVH374oVkTsAcPHkznzp0xGAwMGDCgwDY/Pz92797N2bNncXR0xM3NrdD+77//Pg0aNODVV19l6NCh6PV6tmzZQq9evXBzc6Ny5crMmjULb29vEhISGDly5D1jepjkrjpxV12qJNJes5+mFxYwa/58tvyZbO2QhBCiQps9ezbt2rUrckHJHj16cPDgQQ4cOHDPdqZMmYKrqythYWF06dKFyMhImjZtes/92rVrh7e3N5GRkfj4FJwL9s4776DRaAgKCsLd3b3I+VK1a9dmw4YNHDp0iGbNmhEaGsqvv/6KVqtFrVazaNEi9u/fT/369Xnrrbf4/PPP7xnTw6RS/n3vonio0tPTcXFxIS0tDWdnZ2uHU0h6Vi6/j+tFL/UWLihV6GqYxO8jO1PF0cbaoQkhxAPJysoiPj4ef39/bG1t772DACAjIwMfHx/mzJlD9+7drR2O2e52vi35LpYRJ3FXzrY6tJ0mcFHlSTVVCqNUP7LmcOK9dxRCCFGuGI1GLl26xIcffoiLiwtPPfWUtUOyCkmcxD093aIuVV9YgBE1PTTbuXlgqbVDEkII8ZAlJCRQtWpVfvnlF+bMmYNWWzGnSUviJMxTvTnXm+YviPbslSlkXpVn2QkhREXi5+eHoiicP3+etm3bWjscq5HESZjNrdNoTqoCuKnYMXD6Gp77YTfJN7KsHZYQQgjx0EjiJMym0toQG/IVHXPGszuzGjv+TmHqxr+sHZYQQgjx0EjiJCzSp/1j9A2vh7dL/h0J6w5dwGCUGzOFEGWX3FxeMRTXeZbESVjETq/hgyeD2P5ua16y3czPxvc4dFrmOwkhyp7bjyPJyLjDc+ZEuXL7PP/zMTT3o2JOiRcPTJt3i1e1K3HNS+HQhvfJrDEPO73m3jsKIUQpodFoqFSpEsnJ+Qv72tvbm57nJsoPRVHIyMggOTmZSpUqodE82HeVJE7i/tg6czLsS5rFDKTRld94bexYnEN6M757A2tHJoQQZvPy8gIwJU+i/KpUqZLpfD8ISZzEfWva8inm/9GDgYaljNf9QKe9Nfn7cT8CPZysHZoQQphFpVLh7e2Nh4cHubm51g5HlBCdTvfAI023SeIk7pteqybi1Skkzz+FR9phpuq+YcPRFgS2qWvt0IQQwiIajabYvlhF+SaTw8UD8ansjMfABeRoHQlRn6Lyvq+sHZIQQghRYiRxEg/O1Y9b7T8nT1GTkJ7Hrwcvsv2vK9aOSgghhCh2cqlOFAvX5s/y5j49v563g0UHAfjh+RDaBXlaNzAhhBCiGMmIkyg2z3R4Altd/l8pG3KYs+OMlSMSQgghipeMOIliE1qzMgc+bM+Fvw5hWDyAXxLakJYZgovdgy02JoQQQpQWMuIkipW9Xkvtm/t5RH2e9zU/cXD3FmuHJIQQQhQbSZxE8Wv2EqdcW6NXGfDf8ioDv41m85+XrR2VEEII8cAkcRLFT6XCtue3nFfcqa66Qp9LExm++CAZOXnWjkwIIYR4IJI4iRJRvaoPp1tNJw8tHTV76ZqzhpiTskSBEEKIsk0SJ1FiWrfpgLbDfwAYpf2ZEwdirByREEII8WAkcRIlq/kQrlWP4LASwLp4hVOXb3D+Woa1oxJCCCHui0pRFMXaQVRk6enpuLi4kJaWhrOzs7XDKRGGzHRaT4nlfHr+HCcbrZpVrz9OHS95GLAQQgjrs+S7WEacRInT2DkztO3/HvxbJe8yP2yXxTGFEEKUPZI4iYeiX/MabIwKY8Mj69ls8w4Xj+8kz2C0dlhCCCGERSRxEg9NoGclAnUp2KhymWj4kriT8dYOSQghhLCIJE7i4VGpUD89g6s6H3zVV8hcPIg+M/9g15mr1o5MCCGEMIskTuLhsnMl5ckfyFJ0tFQdpPn5Obz9yyG5bCeEEKJMkMRJPHR1Gj9GXMOPAIjSLiMwfRd74q9ZOSohhBDi3iRxElYR2mMYBA9ErVKYovuGmEN/WTskIYQQ4p601g5AVGAdJnLt/J+MvtCM3cdv0qROEt4utjTyrWTtyIQQQogiyYiTsB6dLY4vrWWvQyuu3sph6E/76fHtTk5dvmHtyIQQQogiSeIkrEqv0/BuZB20ahUeXOdxDvDjH2etHZYQQghRJLlUJ6yud4gvT/rmoP0xAiUzjUHHPDB0q49GrbJ2aEIIIUQBVh9xmjFjBv7+/tja2hIcHMz27dvvWj8mJobg4GBsbW0JCAhg5syZheosW7aMoKAgbGxsCAoKYsWKFRYfV1EUxo4di4+PD3Z2drRu3Zpjx44VqJOUlET//v3x8vLCwcGBpk2bsnTp0vv4FISDuz+6ak2wVeUyPvdz9v8ZT1auwdphCSGEEAXcd+KUk5PDhQsXSEhIKPCyxOLFi4mKiuKDDz4gLi6O8PBwOnbseMd24uPj6dSpE+Hh4cTFxTFq1CiGDRvGsmXLTHViY2Pp06cP/fv359ChQ/Tv35/evXuze/dui447adIkJk+ezPTp09m7dy9eXl60b9+eGzf+N/+mf//+nDx5klWrVnHkyBG6d+9Onz59iIuLs+hzEIBajbr7LK7pvKmhTubWwoE0GLOOJfvOWzsyIYQQ4n8UC506dUp5/PHHFbVaXeClUqkUtVptUVvNmjVThg4dWqCsbt26ysiRI4usP2LECKVu3boFyoYMGaK0aNHC9L53795Khw4dCtSJjIxUnnnmGbOPazQaFS8vL2XChAmm7VlZWYqLi4syc+ZMU5mDg4Myf/78Au24ubkpP/zwwx37/G9paWkKoKSlpZm9T3kWf2SnkvFRFUUZ46x8/UF/JfjTDUpWbp61wxJCCFGOWfJdbPGI08CBA1Gr1axevZr9+/dz4MABDhw4QFxcHAcOHDC7nZycHPbv309ERESB8oiICHbu3FnkPrGxsYXqR0ZGsm/fPnJzc+9a53ab5hw3Pj6epKSkAnVsbGxo1apVgdgef/xxFi9ezLVr1zAajSxatIjs7Gxat25t9ucgCvKrH8rVtl8A8Lr2V5plbGPDsctWjkoIIYTIZ/Hk8IMHD7J//37q1q37QAdOSUnBYDDg6elZoNzT05OkpKQi90lKSiqyfl5eHikpKXh7e9+xzu02zTnu7T+LqnPu3DnT+8WLF9OnTx8qV66MVqvF3t6eFStWULNmzTv2Ozs7m+zsbNP79PT0O9atqKq1HABZf8HOrxmsXcuUvU/SpZGPtcMSQgghLJ/jFBQUREpKSrEFoFIVvHNKUZRCZfeq/+9yc9osjjqjR4/m+vXrbNy4kX379jF8+HB69erFkSNH7hj/+PHjcXFxMb18fX3vWLdCazuW1ND36Z87iu1/X6XlpC2MX3fCdL6FEEIIa7A4cZo4cSIjRoxg69atXL16lfT09AIvc1WpUgWNRlNodCk5ObnQSM9tXl5eRdbXarVUrlz5rnVut2nOcb28vADuWuf06dNMnz6dOXPm0LZtWxo1asSYMWMICQnhm2++uWO/33//fdLS0kyv8+dl8nORNFoqRY6kaWA1ABKuZfBdzBm2nrpi5cCEEEJUZBYnTu3atWPXrl20bdsWDw8PXF1dcXV1pVKlSri6uprdjl6vJzg4mOjo6ALl0dHRhIWFFblPaGhoofobNmwgJCQEnU531zq32zTnuP7+/nh5eRWok5OTQ0xMjKlORkYGAGp1wY9Qo9FgNBrv2G8bGxucnZ0LvMSdjXu6Ad0b+/CG/jdGav/Loj2W3bkphBBCFCeL5zht2bKl2A4+fPhw+vfvT0hICKGhocyaNYuEhASGDh0K5I/OXLx4kfnz5wMwdOhQpk+fzvDhw3nppZeIjY1l9uzZLFy40NTmm2++ScuWLZk4cSJdu3bl119/ZePGjezYscPs46pUKqKiohg3bhy1atWiVq1ajBs3Dnt7e5599lkA6tatS2BgIEOGDOGLL76gcuXKrFy5kujoaFavXl1sn1FF5+tmz+SwHPhzIajh7ZM1SM1oSCV7vbVDE0IIURGV6P19Zvjmm2+UGjVqKHq9XmnatKkSExNj2jZgwAClVatWBepv3bpVadKkiaLX6xU/Pz/l22+/LdTmkiVLlDp16ig6nU6pW7eusmzZMouOqyj5SxKMGTNG8fLyUmxsbJSWLVsqR44cKVDn1KlTSvfu3RUPDw/F3t5eadiwYaHlCe5FliMw08aPFWWMs5L5UWVlyMQflM9WH1Ou38q2dlRCCCHKAUu+i1WKYvls29TUVGbPns2JEydQqVQEBQUxaNAgXFxcij+zK+fS09NxcXEhLS1NLtvdjdFAwjddqX51OxeVyjyV/RlhjR7h675NrB2ZEEKIMs6S72KL5zjt27ePmjVrMmXKFK5du0ZKSgqTJ0+mZs2aFq3jJIRF1BqqDJjHZV01qqquMlM/hU1HEki+kWXtyIQQQlQgFo84hYeHExgYyPfff49Wmz9FKi8vj8GDB3PmzBm2bdtWIoGWVzLiZKErp+CHdpCdxlJDS1LbT2VwyzuvmyWEEELcS4mPOL333numpAlAq9UyYsQI9u3bZ3m0QljCvTb0+pE8tQ37jbWYvPEvZm07zdGLadaOTAghRAVgceLk7Oxc5EN4z58/j5OTU7EEJcRdBbbl1tD9LKM9GTkGxq39kz7fxZKcLpfthBBClCyLE6c+ffrw4osvsnjxYs6fP8+FCxdYtGgRgwcPpm/fviURoxCFuHj4MqJDHWx1aipxA6/cBBbtlcVEhRBClCyL13H64osvUKlUPP/88+Tl5QGg0+l45ZVXmDBhQrEHKMSdDA4PYNAjChlzupJ6K5PX9rnzRpvAuz6yRwghhHgQ97UcAeSvnH369GkURSEwMBB7e/vijq1CkMnhDyjzOsbv26K+dpp9xtpsaT6b+tXd6VDfSxIoIYQQZinRyeG32dvb06BBAxo2bChJk7AeO1fUzy4mQ+1IiPoUAbtG8crP+/lu2xlrRyaEEKIcMutSXffu3Zk7dy7Ozs507979rnWXL19eLIEJYbYqtUjr/D36Vf3oodnOX8aq/LDdlhcf90enue/fDYQQQohCzEqcXFxcTJc9nJ2d5RKIKHW8m3bCkDsR1r3LCN1i4jO82fxnfSLreVk7NCGEEOXIfc9xEsVD5jgVszXvwN7vOWGszrtuX/P8YzV5vFYVfCrZWTsyIYQQpVSJznFq06YNqampRR60TZs2ljYnRPHqMIFrTV/n2ZxRHE26xYhlh+n57U4ycwzWjkwIIUQ5YHHitHXrVnJycgqVZ2VlsX379mIJSoj7ptHi9tR/6NyigakoMS2DXw9etGJQQgghyguz13E6fPiw6efjx4+TlJRkem8wGFi/fj1Vq1Yt3uiEuE+fdqvPyI512bd8Mqrjq5i25zOeaVbd2mEJIYQo48xOnBo3boxKpUKlUhV5Sc7Ozo6vv/66WIMT4kE45FwlPP4r1JqbXLn8OePXelCvaiW6NPSWGxyEEELcF7MTp/j4eBRFISAggD179uDu7m7aptfr8fDwQKPRlEiQQtwXJ0/Uvedh+KkXPTQ7mLrzc4bl9SQtI4f+oX7Wjk4IIUQZZHbiVKNGDQCMRmOJBSNEsQtsx5WW4/Ha9h5R2uWcN3rw3TY7+jWvgVoto05CCCEsY/Gz6m47fvw4CQkJhSaKP/XUUw8clBDFyavNUBTDZVR/TGaC/nueT6vMH6cbEF7L/d47CyGEEP9gceJ05swZnn76aY4cOYJKpeL2MlC354wYDHLbtyh9VG0/hNSz6I4t5zvdFHourEF40wb0CqlGXS9ZP0sIIYR5LF6O4M0338Tf35/Lly9jb2/PsWPH2LZtGyEhIWzdurUEQhSiGKjV0O1bMrweZWZeF05lODB7RzzP/bCHjJw8a0cnhBCijLA4cYqNjeWTTz7B3d0dtVqNWq3m8ccfZ/z48QwbNqwkYhSieOhssX9pHfZt36W6mwMAKTezWXXwkpUDE0IIUVZYnDgZDAYcHR0BqFKlCpcu5X/p1KhRg5MnTxZvdEIUN42O19vUYtuIJ/iovS8jtItYuvtva0clhBCijLB4jlP9+vU5fPgwAQEBNG/enEmTJqHX65k1axYBAQElEaMQxU9ReO7cKPTa7fgnJ/LkVDUNq7sxpks9bHWyrIYQQoiiWTziNHr0aNOSBJ999hnnzp0jPDyctWvXMm3atGIPUIgSoVKhbzOSPJWOjpq99EmZzsI9CczYetrakQkhhCjFVMrt2+IewLVr13B1dZXVmO+DJU9kFsUvff9SHH8bjBqFz3N7s9iuNztHtkWvtfh3CiGEEGWUJd/FxfLt4ObmJkmTKJOcg3ui7jgJgHd1v/BE5gY2HE+6x15CCCEqKrPmOHXv3t3sBpcvX37fwQhhFc1fhptJsP1Lxmt/4OVFzry75FFGdKjDC4/5Wzs6IYQQpYhZiZOLi0tJxyGEdbX5kIyrF7h17HeSFDcycw2MW3uCJxt44+Fsa+3ohBBClBLFMsdJ3D+Z41SKGHL5O/4MMUl65uyI52JqJlHtahHVrra1IxNCCFGCLPkuvu9n1QlR7mh0BAbWITAQqjjqWbB4Edt3p9KlkQ+2Og1VK9lZO0IhhBBWZnHi5O/vf9eJ4GfOnHmggIQoDTo6/kWkzQSSs13o+aVCmrYyK197jEe8ZVRQCCEqMosTp6ioqALvc3NziYuLY/369bz77rvFFZcQVqX3fIR0e0+qZ5znJ/04eud8xNSNp/iuf4i1QxNCCGFFFidOb775ZpHl33zzDfv27XvggIQoFZw8cX5pDYYfIqh96yJz9ZN47vgozl8LwtfN3trRCSGEsJJiW+WvY8eOLFu2rLiaE8L6XGugGfAr2LnRWH2a77STGTBrG5/8dpyEqxnWjk4IIYQVFFvitHTpUtzc3IqrOSFKB4+68NxS8rQOPKY5xnu3PmfeH3/z8oJ9GI1yQ6oQQlQ0Fl+qa9KkSYHJ4YqikJSUxJUrV5gxY0axBidEqVA1GM2zC8lb0BNHOxvUuQp/Jt1gx98ptKztbu3ohBBCPEQWJ07dunUr8F6tVuPu7k7r1q2pW7duccUlRKmiCmiF9qVoHvNqQL/VfzJ351nmx56TxEkIISoYWQDTymQBzLLn7+SbtJ+8hXbqA+zUNqOmhxM/D26Ok63O2qEJIYS4DyW+AKbBYGDFihWcOHEClUrFI488QteuXdFqZT1NUf4Fujvwk8d/eSx9LVNye/DVhR7M3hEvK4wLIUQFYHGmc/ToUbp27UpSUhJ16tQB4NSpU7i7u7Nq1SoaNGhQ7EEKUaqoVDRr/hhEr+Ut3TJy0DFvZw+GtqqJrU5j7eiEEEKUIIsv1bVo0QIPDw/mzZuHq6srANevX2fgwIEkJycTGxtbIoGWV3KprgzbPhk2fQzAJ7n9uVBnII/6udG3eXUcbWT0VQghygpLvostTpzs7OzYt28f9erVK1B+9OhRHn30UTIzMy2PuAKTxKmM2zIeYiYAMDr3BX4ytKdDPS9m9g+2cmBCCCHMZcl3scXrONWpU4fLly8XKk9OTiYwMNDS5oQo21qPJDd0GACf6X6kl2Yr648lcebKTevGJYQQokRYnDiNGzeOYcOGsXTpUi5cuMCFCxdYunQpUVFRTJw4kfT0dNNLiHJPpUIX8Qk0fwXUOmr4+ACwYNc5KwcmhBCiJFh8qU6t/l+udXshzNtN/PO9SqXCYDAUV5zlllyqKycUBS4fY1u6J8/P2QOAq72OhtUq8f3zIei1xbZIvxBCiGJWossRbNmy5b4DE6LcUqnAqz7hngrBNVy5nHCS2lkX2HyqKYv3JtA/1M/aEQohhCgGFidOrVq1Kok4hCgXVCoV83r5opo9CNvMZF7NeYPvttnRt1l1tBoZdRJCiLLuvu6ZTk1NZfbs2aYFMIOCghg0aBAuLi7FHZ8QZY6jmw8EhsORX/hGP43X02DwfEcaVHXh5ZYBssK4EEKUYRb/Crxv3z5q1qzJlClTuHbtGikpKUyePJmaNWty4MCBkohRiLJFrYGnZ0LDPmgxMl03Dbu/VvP15r8ZteKotaMTQgjxACyeHB4eHk5gYCDff/+96REreXl5DB48mDNnzrBt27YSCbS8ksnh5ZjRQO7yoeiO/oIBNW/kvM46pQXb3n0CXzd7a0cnhBDi/5X4AphxcXHUrVu3QPnx48cJCQkhIyPD8ogrMEmcyjmjAX59HQ79FwNqBue8TUBYd1583J/KjnpstPKIFiGEsLYSXQDT2dmZhISEQuXnz5/HycnJ0uaEKN/UGug6HRr3I9OlFoeMNZm9I56wCZvp9s1OsnJlyQ4hhChLLE6c+vTpw4svvsjixYs5f/48Fy5cYNGiRQwePJi+ffuWRIxClG1qDTz1NQ5Df6dhnZqm4hOJ6fx3d+FfQoQQQpReFt9V98UXX6BSqXj++efJy8sDQKfT8corrzBhwoRiD1CIckGtQWXnypwBj/Jn0g2uxHzHiiPXmL3Djv6hNdDJUgVCCFEmWDzH6baMjAxOnz6NoigEBgZiby+TXe+HzHGqgM7+AXM7YUTFB7mD2OX6FO6ONnzeqyE1KjtYOzohhKhwSnSO02329vZUqlQJNzc3SZqEsET1UGj2MmoUxutm88T1pew5e433lx+xdmRCCCHuweLEKS8vjw8//BAXFxf8/PyoUaMGLi4ujB49mtzc3JKIUYjyRa2GjpMwhr0JwEe6BbymWcnO01fZf+66lYMTQghxNxbPcXr99ddZsWIFkyZNIjQ0FIDY2FjGjh1LSkoKM2fOLPYghSh3VCrU7T8GvQNsHce7ul+wU2UzfLEd7YO86NjAm+AartaOUgghxL9YPOK0cOFC5s6dy5AhQ2jYsCENGzZkyJAhzJkzh4ULF1ocwIwZM/D398fW1pbg4GC2b99+1/oxMTEEBwdja2tLQEBAkYnasmXLCAoKwsbGhqCgIFasWGHxcRVFYezYsfj4+GBnZ0fr1q05duxYoXZiY2Np06YNDg4OVKpUidatW5OZmWnhpyAqJJUKWr8HEZ8B8Lr2V6qn7uaHHfE8P3s3SWlZVg5QCCHEv1mcONna2uLn51eo3M/PD71eb1FbixcvJioqig8++IC4uDjCw8Pp2LFjketEAcTHx9OpUyfCw8OJi4tj1KhRDBs2jGXLlpnqxMbG0qdPH/r378+hQ4fo378/vXv3Zvfu3RYdd9KkSUyePJnp06ezd+9evLy8aN++PTdu3ChwrA4dOhAREcGePXvYu3cvr7/+Omq13CElLBD2Bjz5JUcDXuSGTzgAt3IMTNv8l5UDE0IIUYhioY8//ljp27evkpWVZSrLyspS+vXrp4wdO9aitpo1a6YMHTq0QFndunWVkSNHFll/xIgRSt26dQuUDRkyRGnRooXpfe/evZUOHToUqBMZGak888wzZh/XaDQqXl5eyoQJE0zbs7KyFBcXF2XmzJmmsubNmyujR482p6t3lJaWpgBKWlraA7Ujyo/dZ64qdd9bqgS896syb2e8svvMVWuHJIQQ5Zol38UWD43ExcWxevVqqlWrRrt27WjXrh3VqlXjt99+49ChQ3Tv3t30upucnBz2799PREREgfKIiAh27txZ5D6xsbGF6kdGRrJv3z7TxPQ71bndpjnHjY+PJykpqUAdGxsbWrVqZaqTnJzM7t278fDwICwsDE9PT1q1asWOHTvu2u/s7GzS09MLvIT4p2bV7Fjm8hUzdF/xn1/j6P1dLMv2X7B2WEIIIbiPyeGVKlWiR48eBcp8fX0tPnBKSgoGgwFPT88C5Z6eniQlJRW5T1JSUpH18/LySElJwdvb+451brdpznFv/1lUnXPnzgFw5swZAMaOHcsXX3xB48aNmT9/Pm3btuXo0aPUqlWryD6MHz+ejz/+uOgPRQiAxMPUzT1BkCaHH5nEy7nD+WLDSbo08kGvlcvAQghhTRYnTj/++GOxBqBSqQq8VxSlUNm96v+73Jw2H7SO0WgEYMiQIbzwwgsANGnShE2bNjFnzhzGjx9fZPzvv/8+w4cPN71PT0+/r8RTlGPVm6Puvwxl4bOE5RxniXY8/dLeZWXcRUJrVqZqJTvU6jv/GxFCCFFyrPbra5UqVdBoNIVGl5KTkwuN9Nzm5eVVZH2tVkvlypXvWud2m+Yc18vLC+Cudby9vQEICgoqUOeRRx654+R2yL/k5+zsXOAlRCH+LVEN/A3sK/OIcpol+o+Zumwz4ZO2MOSn/aZfGIQQQjxcVkuc9Ho9wcHBREdHFyiPjo4mLCysyH1CQ0ML1d+wYQMhISHodLq71rndpjnH9ff3x8vLq0CdnJwcYmJiTHX8/Pzw8fHh5MmTBdo5deoUNWrUMOszEOKufJrAoN8xOlejpjqRZTZjCVRdIPr4ZTadSLZ2dEIIUTGV6DT1e1i0aJGi0+mU2bNnK8ePH1eioqIUBwcH5ezZs4qiKMrIkSOV/v37m+qfOXNGsbe3V9566y3l+PHjyuzZsxWdTqcsXbrUVOePP/5QNBqNMmHCBOXEiRPKhAkTFK1Wq+zatcvs4yqKokyYMEFxcXFRli9frhw5ckTp27ev4u3traSnp5vqTJkyRXF2dlaWLFmi/PXXX8ro0aMVW1tb5e+//zb7M5C76sQ9pZ5XjF8/quSN81VmLF6l1HhvtfLktG1Kbp5ByckzWDs6IYQo8yz5LrZq4qQoivLNN98oNWrUUPR6vdK0aVMlJibGtG3AgAFKq1atCtTfunWr0qRJE0Wv1yt+fn7Kt99+W6jNJUuWKHXq1FF0Op1St25dZdmyZRYdV1HylyQYM2aM4uXlpdjY2CgtW7ZUjhw5Uqid8ePHK9WqVVPs7e2V0NBQZfv27Rb1XxInYZZbVxXl4gEl5UaWUnf0OqXGe6uVGu+tVlqM26jEX7lp7eiEEKJMs+S7WKUo9z9ZIisrC1tb2+Ia/KqQLHkisxAAU6JPsWvzKtxU6awzNqdVbXfmDWpm7bCEEKLMsuS72OI5TkajkU8//ZSqVavi6Ohoui3/ww8/ZPbs2fcXsRDCbFFNNfzsOIUZ+mkM0Gwg5tQVdp+5au2whBCiQrA4cfrss8+YO3cukyZNKvCIlQYNGvDDDz8Ua3BCiMJUrn5oG/VGhcLHurm8p13IM7N20viTDczbedba4QkhRLlmceI0f/58Zs2aRb9+/dBoNKbyhg0b8ueffxZrcEKIIqg18ORkaDMagFe0vzFZ9y23MjL5ZPVxzl29ZeUAhRCi/LI4cbp48SKBgYGFyo1Go+mxJ0KIEqZSQct3oesMFLWWpzV/8Ivjl9gbb/HdtjPWjk4IIcotixOnevXqsX379kLlS5YsoUmTJsUSlBDCTE36oXp2MegdaZJ3iMHaNfx3dwLhkzbzzKxY0jLklxkhhChOFj9yZcyYMfTv35+LFy9iNBpZvnw5J0+eZP78+axevbokYhRC3E1gOxi4BmK/4cTNQXDiGuevZXL+WiaTo0/ycdf61o5QCCHKjftajuD3339n3Lhx7N+/H6PRSNOmTfnoo4+IiIgoiRjLNVmOQBSnrFwDaw4ncvTCNXbu+oMz6hpsfrs1vm721g5NCCFKLUu+ix9oHSfx4CRxEsVOUWD9SHJ3zyYq5xV+JxRnOx1jugTRtXFVa0cnhBClTomu4xQQEMDVq4XXjElNTSUgIMDS5oQQxc2YB2kX0JHLN/ppvKxaybVb2YxcdoSktCxrRyeEEGWaxYnT2bNnMRgMhcqzs7O5ePFisQQlhHgAGh30ng/NXwFghG4x0x1mk5ebzeTok8ggsxBC3D+zJ4evWrXK9PPvv/+Oi4uL6b3BYGDTpk34+fkVa3BCiPuk1kDHCVC5JqwbQWfDZirrkhi6L4pf9l0guIYr8wY1w9HG4vtDhBCiQjN7jpNanT84pVKpCv3GqtPp8PPz48svv6Rz587FH2U5JnOcRIn7ayMsGQg5Nzhq9KNLzmcoqHm5ZQCjOj1i7eiEEMLqSmSOk9FoxGg0Ur16dZKTk03vjUYj2dnZnDx5UpImIUqjWu3gxd9RKlXH8Pg7jOgQBMBPu86RfCOLrNzCl96FEEIUTe6qszIZcRIPTW4W6GxRFIXOX+/g/KVLpOOIWgVjn6rH86F+1o5QCCGswpLv4vua4HDr1i1iYmJISEggJyenwLZhw4bdT5NCiJKmswXyL7dPbOeKx+LnWGh4gil5PflszQmeqOMh6z0JIcQ9WJw4xcXF0alTJzIyMrh16xZubm6kpKRgb2+Ph4eHJE5ClAH1b+wEVSpvalfQzPEqg1IH8vFvx3i5ZU0CPRxxc9BbO0QhhCiVLF6O4K233qJLly5cu3YNOzs7du3axblz5wgODuaLL74oiRiFEMWt+cvw1HRQ6wjN2sZS/cccP3Gc3t/FEjl1G5fTZb0nIYQoisWJ08GDB3n77bfRaDRoNBqys7Px9fVl0qRJjBo1qiRiFEKUhKb9YcAqsK9CPfU5VtmMJkT1J1duZDN14ylrRyeEEKWSxYmTTqdDpVIB4OnpSUJCAgAuLi6mn4UQZUSNMHh5C3g1oIoqnV/sxvOY+giL957nmy1/88ve8+TkGa0dpRBClBoWz3Fq0qQJ+/bto3bt2jzxxBN89NFHpKSksGDBAho0aFASMQohSlKl6jDod1j5Kuprp7HXhmL8+yaf/34SgH3nrjGpZyMrBymEEKWDxcsR7Nu3jxs3bvDEE09w5coVBgwYwI4dOwgMDGTOnDk0bty4hEItn2Q5AlFqKApkpXI51463Fh/k/LVbpF6/xg3sWR8VTl0v+fsphCifLPkulnWcrEwSJ1FqbfucK9vn8OzNKLJda9OwmgvtHvGkW5Oq1o5MCCGKVYmsHH5bmzZtSE1NLfKgbdq0sbQ5IURplHML4n7CPfcSK/Uf8UhqDKsPJxK1+CCxp69aOzohhLAaixOnrVu3Flr0EiArK4vt27cXS1BCCCvTO8DgzeAXjoMqi+/0U5jgshw1RkavPCITxoUQFZbZk8MPHz5s+vn48eMkJSWZ3hsMBtavX0/VqjKEL0S54VAZ+q+A6DGw6xueyV6Kn90pXrnyKrVHr6NqJTu+fa4pDatVsnakQgjx0Jg9x0mtVpuWIShqFzs7O77++msGDRpUvBGWczLHSZQJR5bCqjcgN4PTRm865EwkFy0B7g5siGqJVmPx4LUQQpQaJfKsuvj4eBRFISAggD179uDu7m7aptfr8fDwQKPR3H/UQojSq0FP8AiCxf3Q1H2Jad7NGLH0MGeu3GJmzGmaVHelYTUXnGx11o5UCCFKlNxVZ2Uy4iTKlJwM0Oc/CPj7bWeYv24rlxU3ctAR6OHIytcew9Hmvp4dLoQQVlOid9XNmzePNWvWmN6PGDGCSpUqERYWxrlz5yyPVghRdvx/0gTwQhNnVjhMZLH+U7y4yt/JN5m8QR7VIoQo3yxOnMaNG4ednR0AsbGxTJ8+nUmTJlGlShXeeuutYg9QCFE6aVPPUEWbSRP132xzGcNj6iPM+SOefj/sYvy6E2TmGKwdohBCFDuLL9XZ29vz559/Ur16dd577z0SExOZP38+x44do3Xr1ly5cqWkYi2X5FKdKNOuxcPi/nD5CEZUTMt7mml53TGipm8zX8Z3b2jtCIUQ4p5K9FKdo6MjV6/mL4C3YcMG2rVrB4CtrS2ZmZn3Ea4Qosxy84fB0dD0edQoRGmXs9Z1MlVIY/He85xMukFGTp61oxRCiGJjceLUvn17Bg8ezODBgzl16hRPPvkkAMeOHcPPz6+44xNClHY6O3jqa3h6FujsqZt5gK/cf8WoQOTUbQR99DujVx4pchkTIYQoayxOnL755htCQ0O5cuUKy5Yto3LlygDs37+fvn37FnuAQogyolEfeGkLBLYnoN8UPJxsTJt+2pXAxhPJVgxOCCGKhyxHYGUyx0mUV2mZufx9+QZXNk7l/b8fwdbFg2ebVad6ZXueauRjWlBXCCGsrUQWwCxKgwYNWLt2Lb6+vg/SjBCiHHKx0xF89Te48BWNbSvzWvrrfBmdBcC1Wzm88Ji/lSMUQgjLPdBzEs6ePUtubm5xxSKEKG+qhUDlWnhxlV9sPuUjl7WoMTI5+hTXbhV+WLgQQpR28oApIUTJ8awHL2+BBr3QYGRQ9k+scJyEfVYyTT+Npu6H6/hh+xlrRymEEGZ7oMQpPDzctBimEEIUycYJun8PXWeAzoFGeYdZb/M+rdSHyMo18tmaExw6n2rtKIUQwiwyOdzKZHK4qFBS/oZlg1ASD3OozXy+Ou3FlpNXqOnuQLfGVant5URkPS9rRymEqGAs+S42K3FatWqV2Qd/6qmnzK4rJHESFVBeNsRvg1rtuZSaScSUbeRl3yKL/OULpvRpxNNNqlk5SCFERVLsiZNaXfCKnkqlKrCY3T9vKzYY5PlUlpDESVR0cXH7CPitOz/qn2Vq6uNUcbRh4/BW6LVq7PUPdOOvEEKYpdgfuWI0Gk2vDRs20LhxY9atW0dqaippaWmsXbuWpk2bsn79+mLpgBCi4miSvAIXYxpRWd+ywOFr8m5epfEn0dQf8zuztp22dnhCCFGAxXOc6tevz8yZM3n88ccLlG/fvp2XX36ZEydOFGuA5Z2MOIkKz2iEXTNg41gw5pKsuPJu7svEGBsBMH9QM1rWdrdujEKIcq1EH/J7+vRpXFxcCpW7uLhw9uxZS5sTQlR0ajWEvQ6DN0KV2niorjNPP5GFVZdgSzYDf9zDsIVx/LTrnDzvTghhdRYnTo8++ihRUVEkJiaaypKSknj77bdp1qxZsQYnhKhAfBrDkG3QfCgAoVdX8FalbRgVWHXoEqNXHuW/exKsG6MQosKzeOblnDlzePrpp6lRowbVq1cHICEhgdq1a7Ny5crijk8IUZHo7KDjRKgdCXu+J7LNWM5sTyA+5RZ7zl7jP2tOUNnBhsqOekJquMrz7oQQD919reOkKArR0dH8+eefKIpCUFAQ7dq1k//E7oPMcRLi3oy52Wz+sj//SYsgXvEGoG+z6ozv3sDKkQkhyoNiX47gTrKysrCxsZGE6QFI4iSEGbZOgK3jycSGCYbnmJfbBlAxrW8TnmrkY+3ohBBlXIlODjcajXz66adUrVoVR0dH4uPjAfjwww+ZPXv2/UUshBB30+Q58G+FHdl8rJnNVq+v8SGFYQvjeOKLrby75BBZubKGnBCi5FmcOH322WfMnTuXSZMmodfrTeUNGjTghx9+KNbghBACAJdq0H8ldJgAWlv8Unex0e49+mi2EJ9ykyX7LzB21TFrRymEqAAsTpzmz5/PrFmz6NevHxqNxlTesGFD/vzzz2INTgghTNRqaPEKDN0Bvs2xVzKZqPuepYHRqFSwaO95en67k7GrjpGcnmXtaIUQ5ZTFidPFixcJDAwsVG40GsnNzS2WoIQQ4o6q1IIX1kHEf8DOjZDub/Jm21oA7Dt3nbk7zzJ4/j5yDUYrByqEKI8sXo6gXr16bN++nRo1ahQoX7JkCU2aNCm2wIQQ4o7UmvxFM0NeAL0Db7RR0KhU+JxezLRzNTh8Ad5dcogWAZUJq1mF6pXtrR2xEKKcsDhxGjNmDP379+fixYsYjUaWL1/OyZMnmT9/PqtXry6JGIUQomh6BwA0ahVvBF6B7V/wlJ09H2b2ZdHBJ1h58BJOtlpWvBpGoIeTlYMVQpQHFl+q69KlC4sXL2bt2rWoVCo++ugjTpw4wW+//Ub79u1LIkYhhLg3R0/wbYYu7xYTdD+w1OFzgmyvcSMrj1d+OsDF1EzSs2Q6gRDiwTzQOk7iwck6TkIUI6MBds+ETZ9AXhaK1o5pxl5My2iPAQ0atYovejXk6SbVrB2pEKIUKdF1nAICArh69Wqh8tTUVAICAixtjhkzZuDv74+trS3BwcFs3779rvVjYmIIDg7G1taWgIAAZs6cWajOsmXLCAoKwsbGhqCgIFasWGHxcRVFYezYsfj4+GBnZ0fr1q05dqzo250VRaFjx46oVCp57IwQ1qTWQOhr8MpO8AtHlZfJm8b5/GT7OaBgMCp8sOIoZ1NuWTtSIUQZZXHidPbsWQyGwgvNZWdnc/HiRYvaWrx4MVFRUXzwwQfExcURHh5Ox44dSUgo+kGe8fHxdOrUifDwcOLi4hg1ahTDhg1j2bJlpjqxsbH06dOH/v37c+jQIfr370/v3r3ZvXu3RcedNGkSkydPZvr06ezduxcvLy/at2/PjRs3CsU1depUWT1diNKkck0Y8Bs8NR1sXQjuMoQjYyNp7u9GRo6BdpNjqPvhOt5ffgSjUQbdhRDmM/tS3apVqwDo1q0b8+bNw8XFxbTNYDCwadMmoqOjOXnypNkHb968OU2bNuXbb781lT3yyCN069aN8ePHF6r/3nvvsWrVKk6cOGEqGzp0KIcOHSI2NhaAPn36kJ6ezrp160x1OnTogKurKwsXLjTruIqi4OPjQ1RUFO+99x6Qnxh6enoyceJEhgwZYtrv0KFDdO7cmb179+Lt7c2KFSvo1q2b2Z+BXKoTooTdugr2bqBScTE1k/HTpnM9y8Afxvzn3L3XoS6vtK5p5SCFENZkyXex2XfV3U4GVCoVAwYMKLBNp9Ph5+fHl19+aXaQOTk57N+/n5EjRxYoj4iIYOfOnUXuExsbS0RERIGyyMhIZs+eTW5uLjqdjtjYWN56661CdaZOnWr2cePj40lKSipwLBsbG1q1asXOnTtNiVNGRgZ9+/Zl+vTpeHl5md13IcRD5FDZ9GNVm2ym2X2P2niZ4x5P8mzCU0xc/yfLDlygups947s3wNPZ1orBCiFKO7Mv1RmNRoxGI9WrVyc5Odn03mg0kp2dzcmTJ+ncubPZB05JScFgMODp6Vmg3NPTk6SkpCL3SUpKKrJ+Xl4eKSkpd61zu01zjnv7z3vF9tZbbxEWFkbXrl3N6jPkj1ylp6cXeAkhHhK1BnW9roCKoOQ1bHd4j6fV2/k7+Qab/0zmpfn7yM6TZ94JIe7M4jlO8fHxVKlSpdgC+PfcIEVR7jpfqKj6/y43p80HrbNq1So2b95sGsky1/jx43FxcTG9fH19LdpfCPEAbJyg0+fw4gZwfwQnQypT9N+yy3sKTWyTOHwhjae+/oPB8/YRc+qKtaMVQpRCZl2qmzZtGi+//DK2trZMmzbtrnWHDRtm1oGrVKmCRqMpNLqUnJxcaKTnNi8vryLra7VaKleufNc6t9s057i3L7slJSXh7e1dZJ3Nmzdz+vRpKlWqVKCdHj16EB4eztatW4vsw/vvv8/w4cNN79PT0yV5EuJh820GQ7ZB7NcQ8zle1/exTH2IlqrPOXkZTl6+wbZTV1g0pAVNq7taO1ohRCliVuI0ZcoU+vXrh62tLVOmTLljPZVKZXbipNfrCQ4OJjo6mqefftpUHh0dfcdLX6Ghofz2228FyjZs2EBISAg6nc5UJzo6usA8pw0bNhAWFmb2cf39/fHy8iI6Otr0GJmcnBxiYmKYOHEiACNHjmTw4MEFYmnQoAFTpkyhS5cud+y3jY0NNjY2d/9whBAlT6uH8Lehfk9Y9x5qG0c+qdeJfWevsyf+GvvOXWfQ3L1EBHkS4O7Ii4/7o9NYPEgvhChvFCtatGiRotPplNmzZyvHjx9XoqKiFAcHB+Xs2bOKoijKyJEjlf79+5vqnzlzRrG3t1feeust5fjx48rs2bMVnU6nLF261FTnjz/+UDQajTJhwgTlxIkTyoQJExStVqvs2rXL7OMqiqJMmDBBcXFxUZYvX64cOXJE6du3r+Lt7a2kp6ffsT+AsmLFCos+g7S0NAVQ0tLSLNpPCFHMcrNMP964HK9s+qSj8tjIH5Ua761Wary3Whm1/LAVgxNClCRLvostflZdcerTpw9Xr17lk08+ITExkfr167N27VrTA4QTExMLrK3k7+/P2rVreeutt/jmm2/w8fFh2rRp9OjRw1QnLCyMRYsWMXr0aD788ENq1qzJ4sWLad68udnHBRgxYgSZmZm8+uqrXL9+nebNm7NhwwacnOR5V0KUS9r/jQQ7xoyljeEPHrfdx4Yqz/PW+XB+3p3Ajaw8vCvZ8myz6tSo7GDFYIUQ1mLxI1d69uxJSEhIodv5P//8c/bs2cOSJUuKNcDyTtZxEqIUunIS1rwNZ/OfKHDVzo+otGfYbmwIgLuTDb+9/jheLrJ0gRDlgSXfxRYnTu7u7mzevJkGDRoUKD9y5Ajt2rXj8uXLlkdcgUniJEQppShw+BfY8AHcyr/D7rBDGJ/mPsfe9ErUdHcgsp4Xj3g707mhtzw9QIgyrESfVXfz5k30en2hcp1OJ2sSCSHKD5UKGvWB1/dBi9dAraXhrZ380PAkTrZaTl+5xYytp3ljYRzfbz9j7WiFEA+JxYlT/fr1Wbx4caHyRYsWERQUVCxBCSFEqWFXCTqMy39wcMM+uLQfwdwXmhER5EmkvxZQmLDuT0avPMKU6FMkp2dZO2IhRAmy+FLdqlWr6NGjB88++yxt2rQBYNOmTSxcuJAlS5ZY9Jw2IZfqhCizjEaUH9pyNjWX16/35ZjiB4BfZXtWvvYYlewLj8wLIUqnEp3jBLBmzRrGjRvHwYMHsbOzo2HDhowZM4ZWrVrdd9AVlSROQpRRSUdhdnvIzUBBxQ7nToy92Z3TGXbU9XIitGZlmlZ3pUsjH2tHKoS4hxJPnETxkcRJiDIs7SJEfwRHlwJg0DszOasLP+S0J5v8EaePn6rHgDA/KwYphLiXEp0c/uGHH2IwFH4IZlpaGn379rW0OSGEKLtcqkLP2fDCOvBqgCYnnXfVP7PTaSQ9axoB+Pi3Y7z28wE+XHmUi6mZVg5YCPGgLE6c5s+fz2OPPcbp06dNZVu3bqVBgwacPXu2OGMTQoiyoUYYvBwDXWeAkzeV3b35/MVOdG9SFaMCa44ksmDXOZ6ZFcvVm9nWjlYI8QAsTpwOHz6Mn58fjRs35vvvv+fdd98lIiKCgQMHsmPHjpKIUQghSj+1Bpr0gzf2Q88fUak1jOvegFHtfFlWdRHBzqmcv5ZJl6930H/2br7fdgaZKSFE2XPfc5w++OADxo8fj1arZd26dbRt27a4Y6sQZI6TEOXclvEQMwGjWs98QwSTs58iHUcA3mxbi7fa17ZygEKIEp8c/vXXX/Pee+/x9NNPs3//fjQaDf/9739p1KjRfQddUUniJEQ5d/kY/P4BnNkCQKbWha2eA3nzdFNy0FHN1Y4qjjZ80rUeDatVsm6sQlRQJTo5vGPHjnz88cfMnz+fn3/+mbi4OFq2bEmLFi2YNGnSfQcthBDlkmc96L8C+i0F97rY5aXR8eJX7HEeSTf1Di5cz+Tg+VT6z97DmSs3MRjl8p0QpZnFI07t27dn3rx5+PgUXJtkzZo1DB48mMTExGINsLyTESchKhBDHsQtgK0T4GYSl2t0Zk/Tz/l++xkOX0gzVetQz4upzzTGVqexYrBCVBxWW8cpJSWFKlWqFFdzFYIkTkJUQDkZsHsm1OsGbgFcuZHNGzNWkJt6if1KHQDaB3kyIrIOznY6PJ1trRuvEOVciV6qA9i+fTvPPfccoaGhXLx4EYAFCxbw559/3k9zQghRsejtIXw4uAUA4O5kw08B0Syz+Zg9/t9TT3OR6OOXaT9lG6HjN/Hz7nNWDlgIcZvFidOyZcuIjIzEzs6OuLg4srPz1yS5ceMG48aNK/YAhRCi3DMa0do6gUqDR+IWVutGMNVmJlW5glGBD1ceZdqmv1i8N4HL8hBhIazK4kt1TZo04a233uL555/HycmJQ4cOERAQwMGDB+nQoQNJSUklFWu5JJfqhBAmKX/B5k/h+K8AKGodu1w7E3WxLZdxA6Cyg54lQ0MJcHe0ZqRClCsleqnu5MmTtGzZslC5s7MzqampljYnhBDitiq1oPd8eGkz+LdCZcwl9OoKxgccop6PMx5ONly9lUOfWbt4+5dDTN5wkoycPGtHLUSForV0B29vb/7++2/8/PwKlO/YsYOAgIDiiksIISquqsEwYBXEb4ddM2jTfQxtbJxIuZnN8OmLOJaqY9mB/GkScedTmT3gUfTa+5qyKoSwkMWJ05AhQ3jzzTeZM2cOKpWKS5cuERsbyzvvvMNHH31UEjEKIUTF5B+e//p/VRz0/OD8A+T8xR6Pnrx7sRXb/0qh2biN2Os0DAjzY0irmlYMWIjyz+LEacSIEaSlpfHEE0+QlZVFy5YtsbGx4Z133uH1118viRiFEEIAZFxFr9WAMZPHkxaw3WY5M4lkVkYHLuHI+HV/kpFjoGMDLzydbHF10Fs7YiHKnftexykjI4Pjx49jNBoJCgrC0VEmKt4PmRwuhLCIosCp9bDlP5B0BIBcnRN7PXrxyukWpP3/c/BsdWq+7RfME3U9rBmtEGWC1RbAFJaTxEkIcV+MRvhzNWwdD8nHAVhd93PeOlKNPKOCooCNVs2gx/2x12noEVwNn0p2Vg5aiNJJEqcyRBInIcQDMRrhxKr8V/cfyFNAAb76/nuWnLU3LWPg4WTDL0NC8aviYN14hSiFJHEqQyRxEkIUu5wMlKkNMGSmsdOpA1/ndGZvqhMOeg1eLrY84u3Mf55ugIudztqRClEqlPgjV4QQQpRit66gcq+DVsmlZfpv/JL9KjMdf8Aj9wKnr9xi9eFEnp+9m7iE65xITEd+fxbCfDLiZGUy4iSEKDFn/4Btn8OZLQAoqLlYtQOvXurA4cz/PZD9yYbeTO3TGJ1GfpcWFZNcqitDJHESQpS4C/tg2xdwah0Ap3ps5OX1N0m+kU1GjgGAR/1cCajiSLCfK72Cq6FSqawZsRAPlSROZYgkTkKIhybxMMTHQNgbpqIzSz/iq4MKq/MexYAGgCEtAxjRoS6KoqCVUShRAUjiVIZI4iSEsJrUBPiqMSgGrul9iKnch5Hxjcgmf+FMO52GUU8+Qv8WNawbpxAlTBKnMkQSJyGE1WReh93f5b8yr+UX6Vz5NrMt8/IiTItpDgzzw9PZltCalWnsW8mKAQtRMiRxKkMkcRJCWF1OBsT9BLFf549CAUadPYtqTmLUQTdTNbUKJvduTLcmVa0VqRAlQhKnMkQSJyFEqWHIg+Mr4Y+pcPUMyvBjLDiYxqYTySgZV9l2wQCoqOygx8lWy9in6tG6jjzSRZR9kjiVIZI4CSFKHUWBq6ehSqDpvfJ9GxKv3eCL9LasNoaSgw6dRsXrT9Sikr2OlrXd8ZdVyUUZJYlTGSKJkxCi1Lt+Fr5pAXmZAOTaubPerjNjLjXnGvn/b9lo1cwe8CiP16pyl4aEKJ0kcSpDJHESQpQJGddg/1zY8z3cuARAnkpPjO0T/Kx5ms0p+f9/uTnocXPQM+7pBjTzd7tLg0KUHpI4lSGSOAkhyhRDLhxbCbu+gUtxAOT2WsAbB6qy/liSqZqNVs2QVjWx02mIrOdJgLujlQIW4t4kcSpDJHESQpRJigLnd8PhX6DT5ygqNUcvpuNydC77/oxnXFIzUnAB8teD+v75EB71d0WnVqNWy6rkonSRxKkMkcRJCFFuGHJhSn24mYRBpWWvfUt+UXdk+RUfID9Z8nCyYWqfxoQFylwoUXpI4lSGSOIkhCg3DLlwdFn+PKiL+0zF521q8fXN1qwyhJGFDXqtmq6NfNBq1PQKqUbT6q5WDFoISZzKFEmchBDl0qU42PMDHF0KeVkA3Gw0iLdv9uP3Y5dN1XQaFZ91q09tTyequdrj7mRjrYhFBSaJUxkiiZMQolzLuAZxC2DvbHjmZwwe9Vm0N4Gb5w6Rd/k4X12sSw46APQaNRN7NuDpJtWsHLSoaCRxKkMkcRJCVAhGI6jV/3u/8jU4+BMZWhd+yQ1nuaodh7PyVyEP9HBEp1Hz4uP+9AyWJEqUPEu+i7UPKSYhhBAV2T+TJoDKNcHJB/sblxioWs1AVnOuShO+vBbG78mPko2ed5Yc4siFVOpXdaGaqz2hNStbJ3Yh/kFGnKxMRpyEEBWWIQ/+js5fWPOvDaAYAUhzrcf0wB/4fnt8ger9mlfnvY510WvU2Oo0VghYlFdyqa4MkcRJCCGAtIsQ9xMcmA8hA6Hlu6w7ksiqA/HUuxrNt8n1uIUdACoVdG9Sjf88XV8SKFEsJHEqQyRxEkKIfzAawJADuvwkiWMrYckA8jR2rDc246fscHYb66KgprqbPb5udgS6O/J2ZB2cbXVWDV2UXTLHSQghRNmk1oDarmCZW020107TmRg662PIdPDlx1uh/Pf6Y/xxzZ0//r7K9r9T6PtodWz1GjrV96KyoyxrIEqGjDhZmYw4CSHEPSgKnN8DB3+Go8sh50Z+MSqWPr6GL/Zkcjk921Tdy9mW8T0aUN3NHh8XO+z0cjlP3J1cqitDJHESQggL5GTAn6vz50PlZsLgaJLTs/h++xlqJixhy1VXNtz0RyH/Lj4XOx1fPdOY1nU8rBy4KM0kcSpDJHESQoj7lJcN2v+/JJd5Hb6oDYYcrum8WJkXxn+zWvC3kr8OVBVHGxxtNAxtVZNnmlW3YtCiNJLEqQyRxEkIIYpB2kXY8h84vsp0KQ/gol1tfkx/lFWGMJLJfyZeMz83qjjpaVrdlRce80ejVlkralFKSOJUhkjiJIQQxSgnA06tg8NL8teIMuYBkNTiI+bzJN/GnOaf33rN/N1oW9cDZzsdXRr54Ggj90xVRJI4lSGSOAkhRAm5dRWOr4AjS6HXXHDy4kRiOtd2/UTl85uYdqUJ0bkNyP3/G8wD3B0YEVkHZzsdQd7OVLLXWzd+8dBI4lSGSOIkhBAP2bynID4GgEyNEwcdw5mX3pTozDoYyL8Dz8VOx/juDWju74aDjVYW2iznJHEqQyRxEkKIh+zSQTiyJP9187Kp+KbGhS2axxidM5C0rDxTuV6j5s12tXilVU3UMh+qXJLEqQyRxEkIIazEaIBzf8CxFXD8V8i4Cv6tyO63gi9+P8l/dydQI/c0fyrVMaLGTqfBTq+hTV0PPuoSJCuVlyOSOJUhkjgJIUQpYMiDs9tBowe/x/LLbiShfFmXLJsqLMsKZkVOcw4otVBQ4+agx9fNHg8nG97rUIdADyfrxi8eiCXfxeqHFNMdzZgxA39/f2xtbQkODmb79u13rR8TE0NwcDC2trYEBAQwc+bMQnWWLVtGUFAQNjY2BAUFsWLFCouPqygKY8eOxcfHBzs7O1q3bs2xY8dM269du8Ybb7xBnTp1sLe3p3r16gwbNoy0tLT7/CSEEEJYjUYLNZ/4X9IEkHwClY0zdtlXeE61nmU2H3PC9R0+d/iZ2pkHOXr+KtHHL9P56x0MWbCPNxfF8fuxJOv1QTwUVk2cFi9eTFRUFB988AFxcXGEh4fTsWNHEhISiqwfHx9Pp06dCA8PJy4ujlGjRjFs2DCWLVtmqhMbG0ufPn3o378/hw4don///vTu3Zvdu3dbdNxJkyYxefJkpk+fzt69e/Hy8qJ9+/bcuJG/PsilS5e4dOkSX3zxBUeOHGHu3LmsX7+eF198sYQ+LSGEEA9VzSfg3b+g72Jo+AzonbDNTKKXYQ2L9J+xLDyRsJqVyco18vuxy/x68BJDFuxn8Ly9fLr6OHN2xJOdZ7B2L0Qxs+qluubNm9O0aVO+/fZbU9kjjzxCt27dGD9+fKH67733HqtWreLEiROmsqFDh3Lo0CFiY2MB6NOnD+np6axbt85Up0OHDri6urJw4UKzjqsoCj4+PkRFRfHee+8BkJ2djaenJxMnTmTIkCFF9mfJkiU899xz3Lp1C63WvLVA5FKdEEKUEblZcHpz/iNf/oqG1/dgtKlEzF9XsI37EbvEPXyfEsQWQ2MysAWgtqcjXRtXxUarJrKeF75u9lbuhChKmbhUl5OTw/79+4mIiChQHhERwc6dO4vcJzY2tlD9yMhI9u3bR25u7l3r3G7TnOPGx8eTlJRUoI6NjQ2tWrW6Y2yA6QM3N2kSQghRhuhsoW4n6DYD3j4Jdq6o1SqeqONB6M2NNE7byDe6aRy2G8p6zxk8b/cHly8n8fnvJ/lszQk6TN3G1I2nWLgngd1nrlq7N+I+We0bPiUlBYPBgKenZ4FyT09PkpKKvkaclJRUZP28vDxSUlLw9va+Y53bbZpz3Nt/FlXn3LlzRcZ29epVPv300zuORt2WnZ1Ndvb/nuKdnp5+1/pCCCFKIfW/xh0ix8GJVXDiN7TX46mbtoNP2MFYWw3HnVrwjvo9/rx8k6kb/zLt8mQDb56o64FOo6J9kCf2evmluyyw+llSqQquiaEoSqGye9X/d7k5bRZXHchPfp588kmCgoIYM2bMHWMHGD9+PB9//PFd6wghhChjfB/Nf7X/BJKPw4nf4MRq1JePUL+qK2t6t2TR3gR2/JVCePLPrEypyrojRtYcSQSgaiU7XmldE0cbLQ2ruRDg7mjlDok7sVriVKVKFTQaTaHRpeTk5EIjPbd5eXkVWV+r1VK5cuW71rndpjnH9fLyAvJHnry9ve8a240bN+jQoQOOjo6sWLECne7u63q8//77DB8+3PQ+PT0dX1/fu+4jhBCijFCpwLNe/qv1SLh2BvKy0ahV9Gteg36BefD1bJ7Vw02NM3E2zVid3Zg1qXUZvTITALUK+jWvQWjNyjjYaHk8sIo8iLgUsdocJ71eT3BwMNHR0QXKo6OjCQsLK3Kf0NDQQvU3bNhASEiIKWG5U53bbZpzXH9/f7y8vArUycnJISYmpkBs6enpREREoNfrWbVqFba2tvfst42NDc7OzgVeQgghyim3APB45H/vjQao3xNsXXA0pBOesZGJhi84ZDeUX52/oJ/HGYwKLNh1jld/PsCAOXt4ctp2ftl3nlWHLnE5Pct6fRGAlS/VDR8+nP79+xMSEkJoaCizZs0iISGBoUOHAvmjMxcvXmT+/PlA/h1006dPZ/jw4bz00kvExsYye/Zs091yAG+++SYtW7Zk4sSJdO3alV9//ZWNGzeyY8cOs4+rUqmIiopi3Lhx1KpVi1q1ajFu3Djs7e159tlngfyRpoiICDIyMvjpp59IT083zVdyd3dHo5HnGgkhhPgX99rQczYYciFhF5xaDyfXobl2mkY5B2j09NtEaJrx/bYzOGWcJ+faObYlBTJiaf5SOHqNmmebV6eRrwueTra0CKgsj4F5yKyaOPXp04erV6/yySefkJiYSP369Vm7di01atQAIDExscDaSv7+/qxdu5a33nqLb775Bh8fH6ZNm0aPHj1MdcLCwli0aBGjR4/mww8/pGbNmixevJjmzZubfVyAESNGkJmZyauvvsr169dp3rw5GzZswMkpf3XY/fv3m9aGCgwMLNCv+Ph4/Pz8iv3zEkIIUU5odOAfnv+K/A+k/AUn10HNJ2ild6BVbXfYHA3bPifH3o44bSN2q5uyOLUOc3caTc3Ur+pMnxBf7PVaHq9VBU/ne1/5EA9GHrliZbKOkxBCiCLFfA57ZsGt5ALFF7W+xNmE8MnNbiRn/29erV6rpk+IL7U9Hanqaker2h4yN8pM8qy6MkQSJyGEEHdkNELSYfh7Y/7r/B5QDGDnxtVXjvHttrMcvZSGX/oBdl61I0H53w1MdTyd6NzQGzu9hvZBntSo7GDFjpRukjiVIZI4CSGEMFtmKpzZCpnXIeSF/DJFQZlSD1X6Ra5ovTlu14RVN+qwKasuqeRPL9GoVXRu6E3VSnb4VXbgqcY+2OpkLu5tkjiVIZI4CSGEeCCZ12HRc3B+FxjzTMUKKs7b1GKTTVs+Tg4vsIu3iy3htapgq9PQpZEPj/q5PeyoSxVJnMoQSZyEEEIUi+wbcG5n/ojUma35C3ECtHiV3bXfYc2RRFR5WXgc/5F1GXU5pvih/P+qRI18K1HFQU8tTycGPeaHRwWbZC6JUxkiiZMQQogScSMJzsTkryPl3TC/7MxWmN8VgEytCyftm7L0WgB/GIKIV7wAFXqtmkB3R+z0Gp5uUpVeIdWw0Zbvy3qSOJUhkjgJIYR4aBJ2wx9TIX475NwosCnDxp0pdq/zfVKtAuU2WjW2Og31qzozpGVNmlSvhE6jLldzpCRxKkMkcRJCCPHQGXLh4gE4swXO7si/W8+QjTJkGwdyqpOakYNybCWGY6vYmlOHXcb/jUgB6DQqujTyoU+ILzY6DbU8HHGwsfrjb++bJE5liCROQgghrC43Cy7shRqPgfr/n8a2Yigc+t+TOW7oqrA1uzY7DY8USqScbLQ827w6tT2dcLHT0bK2O3qt1Z7qZjFJnMoQSZyEEEKUSuf3wt/R+SNSF/aCIafA5jerLyc2USHXYESTkcI1nDD+/2RzbxdbIut5odeqaRHgRuvaHqX60TCSOJUhkjgJIYQo9XIz85OnszvyX7kZ8PJWAIxGhdQZEdhdPcpf+kfYmVuLmKyaxBkDycIGAL/K9lRztaeyo56+zarT3N8Nlar0JFKSOJUhkjgJIYQoc4zG/13SMxrhyzqFHg1jQMMlu9qsymjA59ndCmxzstGiVqsIr1WF/i1q4Olsi6uDHhc7HdYgiVMZIomTEEKIMs9oyF83KmEXJMTCuVi4cQmA3MBI1tafQp5BYX/CdQLjxnPCUJU4YyCnFR/TWlI6jYqnGlUlrGZltBoVYTWr4O5k81DCl8SpDJHESQghRLmjKJB2Pj+Rsq8MgW3zy1PPw9T6pmqZagfiDDU5qNRiT14AccZapOEIgF6jpu0jHrjY6fB1s6dXSDU8nEpmYU5JnMoQSZyEEEJUGOmXYPd3cGEfXDqQP1fqH2JcezLL4SWu3swhPukqgapL/Kn4YkCDVq3Cu5It127mcOCj9sW6KKcl38Vld9EFIYQQQpQtzj7Q/uP8nw15+Zf3Luw1vVo98SSt6rcA4OTuddRZN4pctS2nNIFsz/QnLjWQbcYGJKdn4+tmb5UuSOIkhBBCiIdPo81/FIx3Q3j0xfyyf1wEq2N/C2xc0GWnUc94lHraowAkDNiDt4v1nqUniZMQQgghSod/LlHQoCfU6w4pp/43KnXtDNX9ahes95BJ4iSEEEKI0kmtBo+6+a+m/a0dDQBlZz10IYQQQggrk8RJCCGEEMJMkjgJIYQQQphJEichhBBCCDNJ4iSEEEIIYSZJnIQQQgghzCSJkxBCCCGEmSRxEkIIIYQwkyROQgghhBBmksRJCCGEEMJMkjgJIYQQQphJEichhBBCCDPJQ36tTFEUANLT060ciRBCCFEx3f4Ovv2dfDeSOFnZjRs3APD19bVyJEIIIUTFduPGDVxcXO5aR6WYk16JEmM0Grl06RJOTk6oVKpiazc9PR1fX1/Onz+Ps7NzsbVbGkjfyp7y2i+QvpVV0reyqaT6pigKN27cwMfHB7X67rOYZMTJytRqNdWqVSux9p2dncvdP5zbpG9lT3ntF0jfyirpW9lUEn2710jTbTI5XAghhBDCTJI4CSGEEEKYSRKncsrGxoYxY8ZgY2Nj7VCKnfSt7Cmv/QLpW1klfSubSkPfZHK4EEIIIYSZZMRJCCGEEMJMkjgJIYQQQphJEichhBBCCDNJ4lQOzZgxA39/f2xtbQkODmb79u3WDsliY8eORaVSFXh5eXmZtiuKwtixY/Hx8cHOzo7WrVtz7NgxK0Z8Z9u2baNLly74+PigUqlYuXJlge3m9CU7O5s33niDKlWq4ODgwFNPPcWFCxceYi+Kdq++DRw4sNB5bNGiRYE6pbFv48eP59FHH8XJyQkPDw+6devGyZMnC9Qpq+fNnL6V1fP27bff0rBhQ9MaP6Ghoaxbt860vayeM7h338rqOfu38ePHo1KpiIqKMpWVtvMmiVM5s3jxYqKiovjggw+Ii4sjPDycjh07kpCQYO3QLFavXj0SExNNryNHjpi2TZo0icmTJzN9+nT27t2Ll5cX7du3Nz3CpjS5desWjRo1Yvr06UVuN6cvUVFRrFixgkWLFrFjxw5u3rxJ586dMRgMD6sbRbpX3wA6dOhQ4DyuXbu2wPbS2LeYmBhee+01du3aRXR0NHl5eURERHDr1i1TnbJ63szpG5TN81atWjUmTJjAvn372LdvH23atKFr166mL9myes7g3n2DsnnO/mnv3r3MmjWLhg0bFigvdedNEeVKs2bNlKFDhxYoq1u3rjJy5EgrRXR/xowZozRq1KjIbUajUfHy8lImTJhgKsvKylJcXFyUmTNnPqQI7w+grFixwvTenL6kpqYqOp1OWbRokanOxYsXFbVaraxfv/6hxX4v/+6boijKgAEDlK5du95xn7LSt+TkZAVQYmJiFEUpX+ft331TlPJz3hRFUVxdXZUffvihXJ2z2273TVHK/jm7ceOGUqtWLSU6Olpp1aqV8uabbyqKUjr/rcmIUzmSk5PD/v37iYiIKFAeERHBzp07rRTV/fvrr7/w8fHB39+fZ555hjNnzgAQHx9PUlJSgX7a2NjQqlWrMtdPc/qyf/9+cnNzC9Tx8fGhfv36ZaK/W7duxcPDg9q1a/PSSy+RnJxs2lZW+paWlgaAm5sbUL7O27/7dltZP28Gg4FFixZx69YtQkNDy9U5+3ffbivL5+y1117jySefpF27dgXKS+N5k2fVlSMpKSkYDAY8PT0LlHt6epKUlGSlqO5P8+bNmT9/PrVr1+by5ct89tlnhIWFcezYMVNfiurnuXPnrBHufTOnL0lJSej1elxdXQvVKe3ntWPHjvTq1YsaNWoQHx/Phx9+SJs2bdi/fz82NjZlom+KojB8+HAef/xx6tevD5Sf81ZU36Bsn7cjR44QGhpKVlYWjo6OrFixgqCgINMXaFk+Z3fqG5Ttc7Zo0SIOHDjA3r17C20rjf/WJHEqh1QqVYH3iqIUKivtOnbsaPq5QYMGhIaGUrNmTebNm2ea8Fge+nnb/fSlLPS3T58+pp/r169PSEgINWrUYM2aNXTv3v2O+5Wmvr3++uscPnyYHTt2FNpW1s/bnfpWls9bnTp1OHjwIKmpqSxbtowBAwYQExNj2l6Wz9md+hYUFFRmz9n58+d588032bBhA7a2tnesV5rOm1yqK0eqVKmCRqMplGEnJycXytbLGgcHBxo0aMBff/1luruuPPTTnL54eXmRk5PD9evX71inrPD29qZGjRr89ddfQOnv2xtvvMGqVavYsmUL1apVM5WXh/N2p74VpSydN71eT2BgICEhIYwfP55GjRrx1VdflYtzdqe+FaWsnLP9+/eTnJxMcHAwWq0WrVZLTEwM06ZNQ6vVmmIrTedNEqdyRK/XExwcTHR0dIHy6OhowsLCrBRV8cjOzubEiRN4e3vj7++Pl5dXgX7m5OQQExNT5vppTl+Cg4PR6XQF6iQmJnL06NEy19+rV69y/vx5vL29gdLbN0VReP3111m+fDmbN2/G39+/wPayfN7u1beilJXzVhRFUcjOzi7T5+xObvetKGXlnLVt25YjR45w8OBB0yskJIR+/fpx8OBBAgICSt95K/bp5sKqFi1apOh0OmX27NnK8ePHlaioKMXBwUE5e/astUOzyNtvv61s3bpVOXPmjLJr1y6lc+fOipOTk6kfEyZMUFxcXJTly5crR44cUfr27at4e3sr6enpVo68sBs3bihxcXFKXFycAiiTJ09W4uLilHPnzimKYl5fhg4dqlSrVk3ZuHGjcuDAAaVNmzZKo0aNlLy8PGt1S1GUu/ftxo0byttvv63s3LlTiY+PV7Zs2aKEhoYqVatWLfV9e+WVVxQXFxdl69atSmJioumVkZFhqlNWz9u9+laWz9v777+vbNu2TYmPj1cOHz6sjBo1SlGr1cqGDRsURSm750xR7t63snzOivLPu+oUpfSdN0mcyqFvvvlGqVGjhqLX65WmTZsWuM24rOjTp4/i7e2t6HQ6xcfHR+nevbty7Ngx03aj0aiMGTNG8fLyUmxsbJSWLVsqR44csWLEd7ZlyxYFKPQaMGCAoijm9SUzM1N5/fXXFTc3N8XOzk7p3LmzkpCQYIXeFHS3vmVkZCgRERGKu7u7otPplOrVqysDBgwoFHdp7FtRfQKUH3/80VSnrJ63e/WtLJ+3QYMGmf7vc3d3V9q2bWtKmhSl7J4zRbl738ryOSvKvxOn0nbeVIqiKMU/jiWEEEIIUf7IHCchhBBCCDNJ4iSEEEIIYSZJnIQQQgghzCSJkxBCCCGEmSRxEkIIIYQwkyROQgghhBBmksRJCCGEEMJMkjgJIYQQQphJEichxEOxdetWVCoVqampVjn+5s2bqVu3Lkaj0SrHt7bWrVsTFRV13/sfOXKEatWqcevWreILSogySBInIUSxK+pLOiwsjMTERFxcXKwS04gRI/jggw9Qq/P/25s7dy4qlYpHHnmkUN1ffvkFlUqFn5/fQ46y5CxfvpxPP/3U9N7Pz4+pU6eavX+DBg1o1qwZU6ZMKYHohCg7JHESQjwUer0eLy8vVCrVQz/2zp07+euvv+jVq1eBcgcHB5KTk4mNjS1QPmfOHKpXr/4wQyxxbm5uODk5PVAbL7zwAt9++y0Gg6GYohKi7JHESQhRrAYOHEhMTAxfffUVKpUKlUrF2bNnC12qmzt3LpUqVWL16tXUqVMHe3t7evbsya1bt5g3bx5+fn64urryxhtvFPiizsnJYcSIEVStWhUHBweaN2/O1q1b7xrTokWLiIiIwNbWtkC5Vqvl2WefZc6cOaayCxcusHXrVp599tlC7fz2228EBwdja2tLQEAAH3/8MXl5eabtaWlpvPzyy3h4eODs7EybNm04dOiQafvYsWNp3Lgx3333Hb6+vtjb29OrV6+7Xr68/Tn908qVKwskoLfbXbBgAX5+fri4uPDMM89w48YNU51/jgK2bt2ac+fO8dZbb5nOEcC5c+fo0qULrq6uODg4UK9ePdauXWtqIzIykqtXrxITE3PHeIUo7yRxEkIUq6+++orQ0FBeeuklEhMTSUxMxNfXt8i6GRkZTJs2jUWLFrF+/Xq2bt1K9+7dWbt2LWvXrmXBggXMmjWLpUuXmvZ54YUX+OOPP1i0aBGHDx+mV69edOjQgb/++uuOMW3bto2QkJAit7344ossXryYjIwMID9R6dChA56engXq/f777zz33HMMGzaM48eP89133zF37lz+85//AKAoCk8++SRJSUmsXbuW/fv307RpU9q2bcu1a9dM7fz999/88ssv/Pbbb6xfv56DBw/y2muvmffh3sXp06dZuXIlq1evZvXq1cTExDBhwoQi6y5fvpxq1arxySefmM4RwGuvvUZ2djbbtm3jyJEjTJw4EUdHR9N+er2eRo0asX379geOV4iyShInIUSxcnFxQa/XY29vj5eXF15eXmg0miLr5ubm8u2339KkSRNatmxJz5492bFjB7NnzyYoKIjOnTvzxBNPsGXLFiA/OVi4cCFLliwhPDycmjVr8s477/D444/z448/3jGms2fP4uPjU+S2xo0bU7NmTZYuXYqiKMydO5dBgwYVqvef//yHkSNHMmDAAAICAmjfvj2ffvop3333HQBbtmzhyJEjLFmyhJCQEGrVqsUXX3xBpUqVCiR+WVlZzJs3j8aNG9OyZUu+/vprFi1aRFJSktmfcVGMRiNz586lfv36hIeH079/fzZt2lRkXTc3NzQaDU5OTqZzBJCQkMBjjz1GgwYNCAgIoHPnzrRs2bLAvlWrVuXs2bMPFKsQZZnW2gEIISoue3t7atasaXrv6emJn59fgVEOT09PkpOTAThw4ACKolC7du0C7WRnZ1O5cuU7HiczM7PQZbp/GjRoED/++CPVq1fn5s2bdOrUienTpxeos3//fvbu3WsaYQIwGAxkZWWRkZHB/v37uXnzZqE4MjMzOX36tOl99erVqVatmul9aGgoRqORkydPmhKY++Hn51dgDpO3t7fpczPXsGHDeOWVV9iwYQPt2rWjR48eNGzYsEAdOzs70+icEBWRJE5CCKvR6XQF3qtUqiLLbi8hYDQa0Wg07N+/v9Ao1j+TrX+rUqUK169fv+P2fv36MWLECMaOHcvzzz+PVlv4v0aj0cjHH39M9+7dC22ztbXFaDTi7e1d5Hyrf89R+qfb84vuNGlerVajKEqBstzc3EL17va5mWvw4MFERkayZs0aNmzYwPjx4/nyyy954403THWuXbtWINkVoqKRxEkIUez0en2J3HnVpEkTDAYDycnJhIeHW7Tf8ePH77jdzc2Np556il9++YWZM2cWWadp06acPHmSwMDAO25PSkpCq9XedRmDhIQELl26ZLp0GBsbi1qtLjSKdpu7uzs3btzg1q1bODg4AHDw4ME7tm+uO50jX19fhg4dytChQ3n//ff5/vvvCyROR48epWfPng98fCHKKpnjJIQodn5+fuzevZuzZ8+SkpJSbItO1q5dm379+vH888+zfPly4uPj2bt3LxMnTixw99e/RUZGsmPHjru2PXfuXFJSUqhbt26R2z/66CPmz5/P2LFjOXbsGCdOnGDx4sWMHj0agHbt2hEaGkq3bt34/fffOXv2LDt37mT06NHs27fP1I6trS0DBgzg0KFDbN++nWHDhtG7d+87XqZr3rw59vb2jBo1ir///pv//ve/zJ079x6f1L35+fmxbds2Ll68SEpKCgBRUVH8/vvvxMfHc+DAATZv3lxgnauzZ89y8eJF2rVr98DHF6KsksRJCFHs3nnnHTQaDUFBQbi7u5OQkFBsbf/44488//zzvP3229SpU4ennnqK3bt33/HOPYDnnnuO48ePc/LkyTvWsbOzu+s8qcjISFavXk10dDSPPvooLVq0YPLkydSoUQPIvzS2du1aWrZsyaBBg6hduzbPPPMMZ8+eLXCHXmBgIN27d6dTp05ERERQv359ZsyYccfjurm58dNPP7F27VoaNGjAwoULGTt27F0+IfN88sknnD17lpo1a+Lu7g7kz9l67bXXeOSRR+jQoQN16tQpENvChQuJiIgw9VmIikil/PviuRBClEMjRowgLS3NdBecNYwdO5aVK1cWy6W2hy07O5tatWqxcOFCHnvsMWuHI4TVyIiTEKJC+OCDD6hRo4asen2fzp07xwcffCBJk6jwZHK4EKJCcHFxYdSoUdYOo8yqXbv2HSewC1GRyKU6IYQQQggzyaU6IYQQQggzSeIkhBBCCGEmSZyEEEIIIcwkiZMQQgghhJkkcRJCCCGEMJMkTkIIIYQQZpLESQghhBDCTJI4CSGEEEKYSRInIYQQQggz/R8U7s/k5jln5wAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 600x400 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "time_fs = time_au * 0.02418884254\n",
        "time_meep = time_fs / 0.1\n",
        "initial = population[0]\n",
        "dipole_moment = 0.1 # meep units of mu12\n",
        "frequency = 1.0 # meep units of omega\n",
        "\n",
        "# analytical golden-rule decay rate\n",
        "gamma = dipole_moment**2 * frequency**2 / 2.0\n",
        "\n",
        "# simple exponential decay reference\n",
        "reference = initial * np.exp(-time_meep * gamma)\n",
        "\n",
        "# below is a more accurate reference which works under any initial population\n",
        "#reference = np.exp(-time_meep * gamma) / (\n",
        "#    np.exp(-time_meep * gamma) + (1.0 - initial) / initial\n",
        "#)\n",
        "\n",
        "std_rel = np.std(population - reference) / initial\n",
        "max_rel = np.max(np.abs(population - reference)) / initial\n",
        "print(f\"std_dev={std_rel:.3e}, max_abs_diff={max_rel:.3e}\")\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "plt.figure(figsize=(6, 4))\n",
        "plt.plot(time_meep, population, label=\"Simulation\")\n",
        "plt.plot(time_meep, reference, label=\"Analytical\", linestyle=\"--\")\n",
        "plt.xlabel(\"time (Meep units)\")\n",
        "plt.ylabel(\"excited-state population\")\n",
        "plt.legend()\n",
        "plt.tight_layout()\n",
        "plt.show()\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "mxl-install",
      "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.11"
    },
    "title": "Socket TLS Workflow"
  },
  "nbformat": 4,
  "nbformat_minor": 5
}