{
 "cells": [
  {
   "cell_type": "raw",
   "id": "c2b00ed3-9e2f-40f5-b9fa-890401747f5c",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/restructuredtext",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    ".. _cookbook-rot-broaden:\n",
    "\n",
    "Rotational Broadening\n",
    "=====================\n",
    "\n",
    "This :ref:`cookbook` recipe demonstrates how to synthesize a rotationally broadened irradiance spectrum, using the MSG Python interface. At first glance, this might seem to be a trivial task: simply convolve a non-rotating irradiance spectrum with a rotational broadening function. However, as discussed in chapter 17 of :ads_citet:`gray:2005`, such an approach rests on a number of assumptions:\n",
    "\n",
    "* that line profiles do not change shape across the stellar surface;\n",
    "* that limb darkening depends only weakly on wavelength;\n",
    "* that Doppler shifts arising from the rotational velocity fields can be represented as convolutions.\n",
    "\n",
    "While these assumptions are often reasonable, there are situations where they can lead to incorrect results. A better strategy is to use equation :math:numref:`eq:irrad`. \n",
    "\n",
    "Preparation\n",
    "-----------\n",
    "\n",
    "Download :download:`sg-HeI-6678.h5 <../downloads/sg-HeI-6678.h5>`\\ [#specgrid-file]_ and place it in your working directory. This :f-schema:`specgrid` file spans the HeI :math:`\\lambda6678\\,\\angstrom` line that will be the focus of our irradiance calculations.\n",
    "\n",
    "Initialization\n",
    "--------------\n",
    "\n",
    "To start, import modules and set other parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2eb4244c",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Import standard modules\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import astropy.constants as con\n",
    "import scipy.io as si\n",
    "\n",
    "# Import pymsg\n",
    "\n",
    "import pymsg\n",
    "\n",
    "# Set constants\n",
    "\n",
    "M_SUN = con.M_sun.cgs.value\n",
    "R_SUN = con.R_sun.cgs.value\n",
    "G = con.G.cgs.value\n",
    "C = con.c.cgs.value\n",
    "PC = con.pc.cgs.value\n",
    "\n",
    "# Set plot parameters\n",
    "\n",
    "plt.rcParams.update({'font.size': 12})"
   ]
  },
  {
   "cell_type": "raw",
   "id": "b815eec6-c5c8-43b4-bfb2-a166674fed69",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/restructuredtext",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Define a Synthesis Function\n",
    "---------------------------\n",
    "\n",
    "Now let's create a function that synthesizes an irradiance spectrum for a rotating star. Neglecting any oblateness and gravity darkening due to the rotation, the disk integral in equation :math:numref:`eq:irrad` can expressed in polar coordinates as\n",
    "\n",
    ".. math::\n",
    "   :label: eq:irrad-disk-int\n",
    "\n",
    "   \\irrad(\\lambda) = \\frac{R^{2}}{d^{2}}\n",
    "   \\int_{0}^{2\\pi} \\int_{0}^{1} \\intsy(\\lambda; \\mu; \\vx; \\vproj) \\, \\rho \\, \\diff{\\rho} \\, \\diff\\phi, \n",
    "\n",
    "where :math:`R` is the stellar radius, :math:`\\rho` the radial coordinate, and :math:`\\phi` the angular coordinate, oriented so that the rotation axis lies along :math:`\\cos(\\phi) = 0`. The angle parameter :math:`\\mu` depends on disk position via\n",
    "\n",
    ".. math::\n",
    "\n",
    "   \\mu = \\cos\\left[ \\sin^{-1}(\\rho) \\right].\n",
    "\n",
    "Likewise, the projected velocity :math:`\\vproj` (appearing here as an additional argument of the specific intensity function) depends on position via\n",
    "\n",
    ".. math::\n",
    "\n",
    "  \\vproj = v \\sin(i) \\, \\rho \\, \\sin(\\phi),\n",
    "\n",
    "where :math:`v \\sin(i)` is the projected equatorial rotation velocity of the star.\n",
    "\n",
    "The :py:func:`synthesize_irradiance` function, defined below, implements equation :math:numref:`eq:irrad-disk-int` using a simple numerical quadrature scheme. It takes the following arguments:\n",
    "\n",
    "* ``specgrid`` --- :py:class:`pymsg.SpecGrid` object\n",
    "* ``lam`` --- wavelength abscissa (Å)\n",
    "* ``M`` --- stellar mass (:math:`\\Msun`)\n",
    "* ``R`` --- stellar mass (:math:`\\Rsun`)\n",
    "* ``Teff`` --- stellar effective temperature (K)\n",
    "* ``d`` --- distance to star (pc)\n",
    "* ``vsini`` --- projected equatorial velocity (km/s)\n",
    "* ``nrho`` --- number of :math:`\\rho` (radial) points in quadrature\n",
    "* ``nphi`` --- number of :math:`\\phi` (angular) points in quadrature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "25ee1683-a67f-483d-8278-b021f2033714",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def synthesize_irradiance(specgrid, lam, M, R, Teff, d, vsini, nrho=128, nphi=128):\n",
    "\n",
    "    # Set up a grid of rho and phi coordinates for the quadrature\n",
    "\n",
    "    drho = 1/nrho\n",
    "    dphi = 2*np.pi/nphi\n",
    "\n",
    "    rho, phi = np.meshgrid(np.linspace(drho/2, 1-drho/2, nrho), \n",
    "                           np.linspace(dphi/2, 2*np.pi-dphi/2, nphi))\n",
    "\n",
    "    # Evaluate mu (angle), dA_proj (projected area) and \n",
    "    # v_proj (projected velocity)\n",
    "\n",
    "    mu = np.cos(np.arcsin(rho))\n",
    "    dA_proj = ((rho + drho/2)**2 - (rho - drho/2)**2)*dphi/2\n",
    "    v_proj = vsini*1E5*rho*np.sin(phi)\n",
    "\n",
    "    # Evaluate dOmega (solid angle) and z (redshift) \n",
    "\n",
    "    dOmega = (R*R_SUN)**2/(d*PC)**2*dA_proj\n",
    "    z = v_proj/C\n",
    "\n",
    "    # Set up (uniform) photospheric parameters\n",
    "\n",
    "    x = {'Teff': np.full(nrho*nphi, Teff), \n",
    "         'log(g)': np.full(nrho*nphi, np.log10(G*M*M_SUN/(R*R_SUN)**2))}\n",
    "\n",
    "    # Synthesize the irradiance. The .flatten() calls are because\n",
    "    # MSG expects 1-D arrays\n",
    "\n",
    "    F_obs = specgrid.irradiance(x, mu.flatten(), dOmega.flatten(), z.flatten(), lam)\n",
    "\n",
    "    # Return the irradiance\n",
    "\n",
    "    return F_obs"
   ]
  },
  {
   "cell_type": "raw",
   "id": "4b4697b4-3500-4344-b904-366393dbe6d8",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/restructuredtext",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "The function uses the following algorithm:\n",
    "\n",
    "* divide the disk into a grid of ``nrho`` by ``nphi`` quadrature elements; \n",
    "* evaluate angle parameters ``mu``, projected areas ``dA_proj`` and projected velocities ``v_proj`` of the elements;\n",
    "* evaluate corresponding solid angles ``dOmega`` and redshifts ``z`` (as viewed from the observer at distance ``d``);\n",
    "* pass ``mu``, ``dOmega`` and ``z``, together with photospheric parameters ``x`` and wavelength abscissa ``lam``, to the :py:meth:`pymsg.SpecGrid.irradiance` method, which sums the contribution of each element toward the irradiance.\n",
    "\n",
    "Evaluate the Irradiance\n",
    "-----------------------\n",
    "\n",
    "Using the :py:func:`synthesize_irradiance` function, evaluate the irradiance spectrum spanning the HeI :math:`\\lambda6678\\,\\angstrom` line profile for a B2 main-sequence star:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "54df2467-6dee-4cbe-b4fe-d50c400d1331",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Create the SpecGrid object\n",
    "\n",
    "specgrid = pymsg.SpecGrid('sg-HeI-6678.h5')\n",
    "\n",
    "# Set up parameters (mass, radius and effective temperature from Wikipedia,\n",
    "# https://en.wikipedia.org/wiki/B-type_main-sequence_star; the choices of \n",
    "# vsini and distance are arbitrary)\n",
    "\n",
    "M = 7.30\n",
    "R = 4.06\n",
    "Teff = 20600.\n",
    "vsini = 100\n",
    "d = 10\n",
    "\n",
    "# Set up the wavelength abscissa\n",
    "\n",
    "lam = np.linspace(6670, 6690, 201)\n",
    "lam_c = 0.5*(lam[:-1] + lam[1:])\n",
    "\n",
    "# Evaluate the irradiance spectrum\n",
    "\n",
    "F_obs = synthesize_irradiance(specgrid, lam, M, R, Teff, d, vsini)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "48829726-45f8-48c9-809c-23a80a773426",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/restructuredtext",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "Plot the Irradiance\n",
    "-------------------\n",
    "\n",
    "As a final step, plot the irradiance spectrum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5e460a32-bdb0-4c8d-a3b5-1fe918c8d013",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAHKCAYAAABPM0zGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABqV0lEQVR4nO3deVxU5f4H8M9szMCw7yCbu+KaXrW0XDM107R91VxKy9Ly1zXNyiULs3u76rW6ppa23kozy1a30puYZZLiloqiAoqAMMAwAzPz/P4YZhTZBAbOnOHzfr14KWfO8j2Mh/n4PM95jkIIIUBEREREklBKXQARERFRc8YwRkRERCQhhjEiIiIiCTGMEREREUmIYYyIiIhIQgxjRERERBJiGCMiIiKSEMMYERERkYQYxoiIiIgkxDBWT4WFhZg1axZuueUWhIWFQaFQYP78+ZLU8sMPP6Bfv37w9vZGQEAARo0ahUOHDklSCxEREdUNw1g95ebm4p133oHZbMaYMWMkq2PTpk0YMWIEwsPDsWHDBvznP//B8ePHcdNNN+HkyZOS1UVERETXRsFnU9aP48emUCiQk5ODsLAwzJs3r8lbxzp06ACtVouUlBQoFAoAQHp6Otq1a4e77roLH330UZPWQ0RERHXDlrF6UigUzvBTm08//RQ33HAD9Ho9fH19MWzYMOzfv7/BNeTm5uLYsWMYMWJEhVri4+PRuXNnfPnll7BarQ0+DhERETUehrFG9uqrr+L+++9HYmIiPvvsM3zwwQcoLCzETTfdhMOHDzdo36WlpQAArVZb6TWtVguj0ciuSiIiIjenlroAT3b27FnMmzcPTz75JJYvX+5cPnToULRt2xYLFizAp59+Wu/9R0REIDg4GL/88kuF5fn5+UhNTQVgbz0jIiIi98WWsUb0ww8/wGKxYNy4cbBYLM4vnU6HAQMG4KeffnKuu3btWmfXZ21f+fn5AAClUolp06Zh27ZtePnll5GdnY0TJ07goYcegtFodK5DRERE7ostY43owoULAIBevXpV+fqVQalfv35YtWrVNe3Xx8fH+feXXnoJRUVFWLRoEV566SUAwMiRIzFhwgSsXr0aLVq0qG/5RERE1AQYxhpRaGgoAGD9+vWIj4+vcd22bduibdu2dT6GWq3GG2+8gYULF+LUqVMIDQ1FVFQUhg0bhpYtWyImJqZetRMREVHTYBhrRMOGDYNarcbJkydx5513NuqxfH190aVLFwDAH3/8gW3btuGf//xnox6TiIiIGo5hrAG+++47FBcXo7CwEABw+PBhrF+/HgBw6623IiEhAQsXLsTcuXORlpaG4cOHIygoCBcuXMDevXuh1+uxYMGCBtXw008/4bfffkPXrl0hhMDevXvx2muvYfjw4XjyyScbfI5ERETUuDjpawMkJCQgPT29ytdOnTqFhIQEAPZZ8pctW4Z9+/bBbDYjMjISvXr1wtSpUzFkyJAG1bB7927MnDkTR44cgdlsRtu2bfHII49g+vTp0Gg0Ddo3ERERNT7ZhLG9e/fixRdfxO7duyGEQK9evbBo0SL069fvmrbPzs7GrFmzsHnzZhiNRnTr1g2LFi1qcBgiIiIiaghZzHvw22+/oX///igpKcEHH3yADz74ACaTCUOGDEFycnKt25vNZgwZMgTbtm3DsmXLsGnTJkRERGD48OH4+eefm+AMiIiIiKomi5ax4cOHIyUlBWlpac5pHQoLC9GqVSu0a9eu0qSnV3vrrbcwbdo07N69GzfccAMAwGKxoFu3bvD19cWvv/7a6OdAREREVBVZDOD/5ZdfMHLkyArza/n5+aF///744osvkJWVhaioqGq337hxI9q3b+8MYoB9SoiHHnoIzz//PDIyMq5pPi6bzYbMzEz4+fld83MpiYiISFpCCBQWFiI6OtotJ0OXRRgrLS2t9vmLAHDw4MEaw1hqaipuuummSsu7du0KADh06NA1hbHMzEzExsZea9lERETkRs6ePeuW82/KIowlJiZiz549sNlszkRrsVic3Yu1PX8xNzcXwcHBlZY7llW3vdlshtlsdn7v6NE9e/Ys/P39634iRERE1OQMBgNiY2Ph5+cndSlVkkUYe+qppzBp0iQ8+eSTmDt3Lmw2GxYsWOCcVuJamhxr6las7rWkpKQq5wHz9/dnGCMiIpIZdx1i5H4dp1WYOHEiFi9ejA8++AAxMTGIi4vD4cOH8eyzzwJArV2MISEhVbZ+5eXlAUCVrWYAMGfOHBQUFDi/zp4928AzISIiIqpIFmEMAJ577jnk5OTg4MGDOH36NHbv3o1Lly5Br9ejZ8+eNW7bpUsXHDx4sNJyx7LOnTtXuZ1Wq3W2grE1jIiIiBqDbMIYYA9HnTt3Rnx8PM6cOYNPP/0Ujz76KLy9vWvcbuzYsTh69GiFKSwsFgs+/PBD9OnTB9HR0Y1dOhEREVGVZDHPWGpqKjZs2IC//e1v0Gq1+PPPP7F48WIkJCRgx44d8PX1da47adIkrFu3DidPnkR8fDwA+0D8nj17wmAwYPHixQgPD8dbb72Fr7/+Glu3bsWAAQOuqQ6DwYCAgAAUFBSwlYyIiEgm3P3zWxYD+L28vLB9+3YsX74cRUVFiIuLw9SpUzF79mzo9foK61qtVlitVlyZMbVaLbZt24ZZs2bhqaeegtFoRPfu3fHdd99dcxAjIiIiagyyaBlzF+6erImIiKgyd//8ltWYMSIiIiJPwzBGREREJCGGMSIiIiIJMYwRERERSYhhjIiIiEhCDGNEREREEmIYq4dzeUYUGMtgtXFWECIiImoYWUz66m6GL9sFpdYHAKD3UsFPp4GfTg0/nRr+3pqK31f6++Xv/XQa+GnVUCrd8ynyRERE1PgYxurBS62EpfzvxaVWFJdacd5Q//35aisHt8uhTeMMef5XvGazCZgtNpgtNpjKrCi12OCrUyPUV4twPy2iA72hYsgjIiJye5yBvw6unMFX66NHoclS/lXm/NNw1TJDSflr5rIK6xtMFpRabI1Wq0alQGywDxJC9IgL9kFCiA/iQ/VICNEjJsgbGhV7qImIqHlw9xn42TJWT1q1ClpfFUJ9tfXeh9lirRjYrgxxpsvL7H+//H2R2QKFwl6DTqOEVq2ERqVEocmCnCIzsg1mlFptSLtYjLSLxZWOq1Iq0CLQG/EhPogPsQe2+BA9EkJ8EBvsA51G1ZAfDREREdUBw5iEXBHoqmK1CWQVlOBMrhGnc41Izy3G6dxipOcacTq3GKYyG87kGXEmz4hdxytvHxWgqxDSHKEtPkQPXy3/yRAREbkSuynrwN2bOa+FEAIXC804XR7M7EHNaA9uOcUoNFtq3D7AW4PoQG9EB+gQHeiNqEAdogO87X8P0CEyQMcuUCIicivu/vnNZo5mRqFQINxfh3B/HXq3DK7wmhACl4xlzpCWnmt0tqal5xqRV1yKgpIyFJSU4UhW1XcsKBRAuJ8WUQHeaBHojQh/HSL8tQj31yLCT4dwfy3C/XXw06qhUPAGAyIiIoYxclIoFAjWeyFY74UecUGVXi80lSEz34TMghJk5ZuQmV9y+e/lf5ZabbhgMOOCwYyUs/nVHkunUSLCX4dwP3s4C/fTOr93/Bnmp0WAt4ahjYiIPBrDGF0zP50G7SM1aB/pV+XrNptAbnEpsgpK7KEtvwQXCk3INpiRXWjCBYMZ2QYTDCYLTGU2Z8tbTdRKBUJ8vRCi1yLUT4tQXy+E+tr/dCwL0XshzE+LYL0Xu0iJiEh2GMbIZZRKBcLKW7S6xlS/XkmpFdmFJmQXmnHBYA9rFwpNuFj+Z7bBvtxgssBiE86WNmTVXkOgjwYhekdgKw9t5X8PuSLIhfpq4eOlYqsbERFJjmGMmpy3l6r8Lk19jeuZLVbkFZcip7AUOUXm8q9S5Jb/Pbe4FBcL7cvyis2wCSDfWIZ8YxlOVjGlx9V0GmV5SNMizNnSdkUrnN7L2fIW5OPFJyUQEVGjYBgjt6VVqxAV4I2oAO9a17XZBC4ZS5FbXIqcQjNyHH8WmZFbVIrcYjMuXhHkTGU2mMpsOHepBOculdS6f5VSgSAfr4rdpFe0uIWV/z0yQIcQPYMbERFdO4Yx8ghKpQIh5a1c7SKqHtPmIISAsdTqbGlzBDb7n5eXOVrf8ssfCu9YBhTWuH+NSoEIfx2iyqf/iA/2QUKoHonR/mgX7segRkREFTCMUbOjUCig16qh16pr7SoFgFKLDZeM9i5RR8tbbnHlIHex0IyLRWaUWcUVLW6XKuwr0EeD3gnBGNElEkMTIzmJLhERcdLXunD3SeNIemVWG7ILzThfUIKsAhMyLpXgdK4RaReLcDCjAMZSq3NdnUaJIR0iMKpbFAa2D+djqIiIGom7f34zjNWBu7+Z5N7KrDakZhRgx7GL+PrPTJzKuXyTga9WjQHtwjCkYzgGtQ9HkN5LwkqJiDyLu39+M4zVgbu/mSQfQgikZhjw9YFMfP1nJrIKTM7XlArgb/HBGJoYgZFdoxAdWPsNDEREVD13//xmGKsDd38zSZ5sNoE/z+Vj25FsbD1yAUfPX75BQKEAeicEY8x1LXBr5ygE+GgkrJSISJ7c/fObYawO3P3NJM9w7pIR249m45sDWfj1VJ5zuZdKicEdwjHuhnjc0DqEE9YSEV0jd//8ZhirA3d/M8nzZOSX4KuUTGxKyajQYtYh0g8T+iXg9u4tOPCfiKgW7v75zTBWB+7+ZpJnO5JlwEe/pmPDvgyUlNnvygzWe+HBPnGYdGNLBPpw0D8RUVXc/fObYawO3P3NpOahwFiGT38/g3W705GRb396gK9WjYk3tsTkm1rCX8dxZUQ1EULAbLGhyGxBkcli/9NsQUmZFeYyG0qtNpjLrDBbbDBbbCi12GC2WCv+vXw9q03AJgSsNgGrDc6/X152xd8FoIB9LKhSoYBSASiguPy98qrvFfZ5ER1/KmBfrlLav9QqBdRKBdQqpf1PpfLysvLlKuVV6zhfd6xrX0ejsq+jUSngpVLCS62EpvxPL9Xlv2tUCvufSqWsJrB2989vhrE6cPc3k5oXi9WGHw9fwL+3n8CRLAMAINTXC38f1h539YyFSka/KInqq6TUivMGE/KK7RMwXzKWIq+4DHnFZuQVl+GSsRQFJWUVQlex2QKLjR99DaVRKa4IacorQpyi0rIrA51Oo4ROo7J/qZXQalTQqq9YplFCp77i7+V/aq9aplEpr7lWd//8ZhirA3d/M6l5stkEvks9j3/+eAxp5XOXdWkRgMV3dkGn6ACJqyOqPyEELhnLcCqnCGfyjMjMNyGroARZ+SZkFphwvqAEl4xlDTqG3ksFvVYNX50a3uVhwEulhFajhFZtDwBe6st/dyx3BAxHK5WzxUqhgFKpgEqJysvKb7qxCQFRfn42AQhhX2YTAkIAAgK28lY25/fOdezXvMUmYLXZUGa1t7pZrDaUlbfElZW32Nlfs8Fita9vqfB3+zaOP602gVKr/e+lVhvKLPY/Sy32YzhaAd2JSqmA7ooQp3WGuMrBDmUl+NfDfd3285thrA4YxsidlVpseD/5NJZtPY5CswVqpQKPD2yNJwe3gVbNQf7k3gpKynAoowAHMwpw7EIh0i4W41ROMQpKag9bPl4qhPpqEaT3QrCPBsF6LYL1mvLvvRDgrYGfTgO9VgU/nf1RaL5aNfReall1tUnNEfQuhzQbyiwCpVYrSi32wFZWHuTMVwS6svL1S63icvdvmRUmiw2mMmv5V/nfy5eZHcss9tfNznVt9ardZjbi7NJ73Pbzm2GsDhjGSA4uFprx0qZUfJd6HgDQuYU/3nygxzU9h5OoKdhsAn9lF2L3iVzsS7+EgxkFOJNnrHb9FoHeiAv2QXSgN6ICdIgK1CE6wBtRgTpEBXjDX6fmVC/NhGO8n/mKoOYMclcFvCvXyc834Lnbe7jt5zfDWB0wjJGcfHswC3M3HsQlYxn8tGq8fndXDO8cJXVZ1ExlFZRg+9Fs7D6Ziz0nc5FbXFppndhgb3SODkBilD9ah/uiZageCSF6eHuxZZcaxt0/vxnG6sDd30yiq2Xml+CpT/ZjX/olAMDfh7XHEwNbsxWBmsT5AhO+PZiFbw5mOf8NOnhrVOjVMhh9WgajW0wgOrfw5/Qs1Gjc/fObYawO3P3NJKpKmdWGpG+P4t1fTgEA7u8dh5dv7wR1He5EIrpWxlILvkrJxIY/zuG305cDmEIB9IwLwk1tw9C3TQi6xQTCS81/g9Q03P3zWy11AUTUuDQqJV4alYj4EB/M//oQPtl7BgUlpVh+33UMZOQymfklWL3rFD7fdxaFJotz+d/igzCyaxRGdI5CZIBOwgqJ3JdsfhPv378fY8aMQXR0NHx8fNChQwcsXLgQRmP1gz4dduzYgaFDhyI8PBy+vr7o2rUrli9fDqvV2gSVE7mH8X0T8J+HesJLpcS3B8/juQ0HYXOzW9VJfs7kGjFr/Z/ov2QH3v3lFApNFsSH+GDOiA5InjMY6x/viwn9WjKIEdVAFi1jhw8fRt++fdG+fXssXboUoaGh2LlzJxYuXIh9+/Zh06ZN1W67detWDBs2DP3798eqVaug1+vx1VdfYcaMGTh58iSWLVvWhGdCJK1hnSLx7weuwxMf/YENf5yDXqvCgtGdOIaM6izfWIoV209gXfJplFntob5v6xA81r8V+rcN45QRRHUgizFjL7zwAl555RWcOHECrVu3di6fMmUK3nnnHeTl5SEoKKjKbR966CGsX78eubm50Osv39o/bNgw7NmzBwUFBddch7v3ORNdq00pGXj60xQIAbxxTzfc0SNG6pJIJoQQ2JSSiQVfH3JOuHpT21A8M7QdesRV/XuYSGru/vkti5Yxjcb+rL2AgIqziQcGBkKpVMLLq/o7cDQaDby8vODt7V1pW52OzebUPN3evQXO5Brxzy1/Yd6mQ+jdMhgxQT5Sl0Vu7oLBhDlfHMT2o9kAgHYRvpg7MhED2oVJXBmRvMlizNj48eMRGBiIxx9/HGlpaSgsLMTmzZuxcuVKTJs2rUKL19WmTp2K0tJSTJ8+HZmZmcjPz8cHH3yAjRs3YtasWU14FkTu5fGBrdEzPgiFZgtmfvan2z3qhNzL7hM5GLl8F7YfzYaXSolnb2mHb6bfxCBG5AKy6KYEgKNHj2Ls2LE4evSoc9n06dOxdOnSWse77N69G3fffTcyMzMBACqVCklJSfj73/9e43Zmsxlms9n5vcFgQGxsrNs2cxLV1ZlcI0Ys24niUivmjUrEhH4tpS6J3IwQAit3pmHJ90dhE0CHSD/8+/7r0DbCT+rSiK6Zu3dTyqJl7PTp0xg1ahRCQkKwfv16/Pzzz1iyZAnWrl2LyZMn17jtvn37MHbsWPTs2RNff/01tm/fjjlz5uCFF17Ayy+/XOO2SUlJCAgIcH7Fxsa68rSIJBcX4oPZIzoAAFbvOsXWMarAZhNYuPkwFn9nD2J39ojBxif6MYgRuZgsWsbuu+8+7NixA2lpaRW6JN977z1MnDgRP/30EwYMGFDlttdffz2MRiP2798PleryIzXmzZuHRYsW4fjx42jVqlWV27JljJoDU5kVNyRtwyVjGd55uCdu6RQpdUnkBqw2gTlfHMBnv58DALx0WyIm9EvgnbckS2wZc4GUlBQkJiZWGhvWq1cvAEBqamqN2/bs2bNCEHNsa7PZcOTIkWq31Wq18Pf3r/BF5Gl0GhXu7RUHAHg/OV3iashdvLz5MD77/RyUCuAfd3fDxBtbMogRNRJZhLHo6GgcOnQIRUVFFZYnJycDAGJiqr8tPzo6Gr///nulCV6vZVui5uLBPnFQKoD/ncjBieyi2jcgj/bVn5lYu/s0AGD5/dfhrp78PUnUmGQRxp5++mnk5ORg6NCh+Oyzz7B9+3a8+uqrmDlzJhITEzFixAgAwKRJk6BWq5Gefvl/98888wxSU1MxatQobNq0CVu2bMHs2bOxZMkS3HzzzejWrZtUp0XkNmKDfTC4QwQA4MM9bB1rzk5kF2L2hgMAgCcHtcFtXaMlrojI88kijI0ePRrbtm2Dv78/ZsyYgdtuuw3r1q3DlClTsHPnTuc8Y1arFVarFVcOg3vqqaewYcMGFBYWYvLkyRg7diw2b96MefPm4csvv5TojIjcz/i+8QCA9fvOwVTGR4U1RxarDdM+2g9jqRV9W4fgmaHtpC6JqFmQxQB+d+HuAwCJGsJmE+i7eDvOG0z4aHIf9GsTKnVJ1MQ2H8jEkx/vR5CPBj8+MwBhflqpSyJyCXf//JZFyxgRNT6lUoG+rUMAALtP5khcDTU1IQRW7ToFAHj4hgQGMaImxDBGRE7Xl4ex5JO5EldCTW1f+iX8eTYfXmolHr4+XupyiJoVhjEicnK0jP15rgBFZovE1VBTWl3eKja2ewu2ihE1MYYxInKKCfJBXLAPrDaB307lSV0ONZH03GL8cPg8AGDyTXwkFlFTYxgjogpuaFXeVZnGrsrm4tPfzkIIYEC7MD7qiEgCDGNEVEHfNhzE39zsOm5/r8dcxznFiKTAMEZEFThaxg5lGlBgLJO4Gmpsl4pLkZpZAADo15rTmRBJgWGMiCoI99ehdZgeQgB7TrGr0tPtPpkLIYB2Eb4I99dJXQ5Rs8QwRkSV3FB+V+VeDuL3eP87cREAcGObMIkrIWq+GMaIqJKuMYEAgCNZBmkLoUb3vxP28WI3tg2RuBKi5othjIgq6Rhpf1zIkSwD+MQ0z3Um14izeSVQKxXo05JhjEgqDGNEVEnbCF8oFcAlYxmyC81Sl0ONZFd5F2WPuCDotWqJqyFqvhjGiKgSnUaFVmG+ANhV6cl+Ke+i5EPhiaTFMEZEVeoQaZ/880hWocSVUGOw2QR2lz+D9Ma2DGNEUmIYI6IqdYyyjxs7ep4tY57o7CUj8o1l8FIr0TUmQOpyiJo1hjEiqlLHKHvL2FG2jHmkw5n2kN0uwhcaFT8KiKTEK5CIqtSh/I7KkxeLYLZYJa6GXM0xFjCxvAWUiKTDMEZEVYoK0CHAWwOLTeBEdpHU5ZCLHWYYI3IbDGNEVCWFQuEcxM+uSs/juDGjI8MYkeQYxoioWo4Pak5v4VnyjaXIyC8BAHSMZhgjkhrDGBFVy9kydp4tY57E0UUZG+wNf51G4mqIiGGMiKrF6S08k7OLMpKtYkTugGGMiKrVLsLeMpZTVIq84lKJqyFXcUxrkcguSiK3wDBGRNXy9lIhKkAHADiVUyxxNeQqvJOSyL0wjBFRjRJC9ACA0wxjHqHUYsOJbN5JSeROGMaIqEYJoeVhLJdhzBOcyC5CmVXAT6dGTJC31OUQERjGiKgWLUN9ALCb0lNcOfO+QqGQuBoiAhjGiKgWzm5Ktox5hL8u2LsoHdOWEJH0GMaIqEaObsr0HCOEEBJXQw3laOFsFeYrcSVE5MAwRkQ1igv2gUIBFJotyOX0FrLnCGMty0M2EUmPYYyIaqTTqBAdYB/ozTsq5c1qE0jPMwJgGCNyJwxjRFSrBA7i9wiZ+SUotdjgpVIiOpB3UhK5C4YxIqoVB/F7BkeYjg/xgUrJOymJ3AXDGBHVytGldTrHKHEl1BCOMJ3ALkoit9IoYcxms8Fo5C9tIk/BljHPkHax/E5KhjEit+KSMGYymbB27VrcfffdiI6OhpeXF/z8/ODj44O//e1vmDVrFv78809XHIqIJOAYM3Y6p5jTW8gY76Qkck/qhmxcUlKCJUuWYNmyZSgoKECHDh0wZMgQhIeHQ6fTIS8vD2lpaVi1ahX++c9/om/fvliyZAluuOEGV9VPRE0gNtgHSgVQXGrFxSIzwv10UpdE9cBuSiL31KAw1rZtW+j1erzwwgt48MEHERERUeV6Qgjs2LED7733HgYNGoQVK1Zg8uTJdTrW/v37sWDBAuzduxf5+fmIi4vDAw88gGeffRY+Pj61br9p0ya88cYb2L9/P6xWKxISEjBjxgw89thjdaqDqDnSqlWIDvTGuUslOJ1jZBiToVKLDWfLp7VgNyWRe2lQGFu4cCHGjx8PlUpV43oKhQKDBw/G4MGDsWDBApw5c6ZOxzl8+DD69u2L9u3bY+nSpQgNDcXOnTuxcOFC7Nu3D5s2bapx+8WLF2Pu3LmYOnUq5syZA41Gg6NHj6K0lBNYEl2rlqH68jBWjN4tg6Uuh+roTJ4RNgHovVQI89NKXQ4RXaFBYWzixIl13qZVq1Zo1apVnbb5+OOPYTKZsGHDBrRu3RoAMHjwYGRlZeGdd97BpUuXEBQUVOW2+/btw9y5c5GUlIRZs2Y5lw8ZMqTOtRM1Zwkheuw6nsNB/DLlmLA3IVTPB4QTuRmX3U1ps9nw/vvvu2p3FWg0GgBAQEBAheWBgYFQKpXw8vKqdtsVK1ZAq9XiqaeeapTaiJqLmCD7JKEZ+SUSV0L1wcH7RO7LZWGsrKwMEyZMcNXuKhg/fjwCAwPx+OOPIy0tDYWFhdi8eTNWrlyJadOmQa+v/pfLzp070bFjR2zYsAHt27eHSqVCTEwMZs+ezW5KojqICbKPzTx3iWFMjtJyOK0FkbuqUzflwoULq32trKyswcVUJyEhAcnJyRg7dqyzmxIApk+fjqVLl9a4bUZGBi5evIjp06fj5ZdfRmJiIrZt24bFixfj7Nmz+Oijj6rd1mw2w2w2O783GAwNPhciuWrhaBljGJOlK7spici91CmMLVq0CGPHjoW/v3+l16xWq8uKutrp06cxatQoREREYP369QgLC8Ovv/6KRYsWoaioCGvWrKl2W5vNhsLCQnzyySe47777AACDBg1CcXExli5digULFqBNmzZVbpuUlIQFCxY0yjkRyU2L8mcZXig02Z9vqOYDPOSE3ZRE7qtOYaxLly4YN24cRo4cWek1x8SvjWH27NkwGAxISUlxdkn2798foaGhmDhxIsaNG4cBAwZUuW1ISAjOnz+PYcOGVVg+YsQILF26FH/88Ue1YWzOnDmYOXOm83uDwYDY2FgXnRWRvIT6ekGrVsJssSGroATxIfxQlwtTmRXnDSYA4PtG5Ibq9F/bRx99tNoWMI1Gg3nz5rmkqKulpKQgMTGx0tiwXr16AQBSU1Or3bZr165VLnfMIq5UVv8j0Gq18Pf3r/BF1FwpFAp2VcpUZvlNF3ovFYJ8NBJXQ0RXq1MYmzp1KkaPHl3layqVqtHCWHR0NA4dOoSioqIKy5OTkwEAMTEx1W575513AgC+++67Csu//fZbKJVKZ6AjotpxEL88Od6vmCAfTmtB5IYaNM9YU3n66acxZswYDB06FM888wxCQ0OxZ88eJCUlITExESNGjAAATJo0CevWrcPJkycRHx8PAJgwYQJWrlyJJ554Ajk5OUhMTMTWrVvx5ptv4oknnnCuR0S1c4wbO8fpLWTlchjzlrgSIqpKg0fgNlZr2JVGjx6Nbdu2wd/fHzNmzMBtt92GdevWYcqUKdi5c6dznjGr1Qqr1VrhQcYajQZbtmzBfffdh1dffRW33norNm7ciMWLF2PZsmWNXjuRJ4lhN6UsnbtkfwwSwxiRe1KIK5NLPfj4+MBoNLqqHrdmMBgQEBCAgoICjh+jZmlTSgZm/DcFfVoG49MpN0hdDl2j6Z/sx1d/ZmLurR3xaP+6PQGFyBO4++d3g1vGGpjliEhGHN2UnIVfXtgyRuTeGhzGOBiUqPlw3E2ZVWCCxWqTuBq6VlcO4Cci98NZG4nomoX76aBRKWC1CVwoNNe+AUnOVGZFdvl7xZYxIvfEMEZE10ylVCAqgIP45eTKOcYCOccYkVvimDEiqhNH64pjHBK5N0cXZYsgbw4rIXJTDQ5j/fr1c0UdRCQTzkH8bBmTBY4XI3J/DQ5jW7dudUUdRCQTLZwtYwxjcsA7KYncH8eMEVGdOFpYOL2FPHD2fSL3xzBGRHXCucbk5XLLGLspidxVk4Wxv/76q6kORUSN6MpHItlsvIHH3bFljMj91etB4YcOHYLFYkG3bt0qvWYymXDgwAGkpKQ4vw4ePAij0Qir1drggolIWhH+OigUQKnVhjxjKUJ9tVKXRNWoOMcYW8aI3FWdw9g//vEPLFmyBEFBQbj33ntx4403Vghex48fh81mgxACCoUCQUFB6N27d5XBjYjkx0utRKivFhcLzThfYGIYc2OOOcZ8vFQI4hxjRG6rzmEsKSkJ77//Pjp16oRWrVo5561RKpVo27Yt7rzzTnTt2hUtW7bEgw8+iC+++AL9+/d3eeFEJJ2oAB0uFpqRVWBC5xYBUpdD1biyi5JzjBG5rzqHMZ1Oh6KiIhiN9kGhQgjcdtttWLduHYKCgpzrFRQUuK5KInIrUQE6HDhXgPMFHMTvzjjHGJE81DmMLVu2DFOmTAEAPPLII7DZbPjggw/Qo0cPvP7667jrrrtcXiQRuRfHI5EyC0wSV0I1ySoPy1EBOokrIaKa1Pluyrvuugvp6ek4evQo3n33XaxduxY7d+5EYGAg7r33XgwZMgSHDh1qjFqJyE1Eln+4n2cYc2uO9yc6kHdSErmzek1t4evri7CwMOf3/fr1w759+7B8+XKkpKTguuuuwzPPPMMxCkQeytHSksVuSrd23mAPYxH+bBkjcmcum2dMqVRi2rRpOHbsGB5++GGsW7cOALBr1y5XHYKI3ESkP1vG5CCr/P1hNyWRe3P5pK+hoaFYs2YNkpOT0aNHD7z44ovo378/jhw54upDEZFEHGPGsgpMEIITv7orR1iOZBgjcmuNNgN/7969sXfvXrz99ts4fPgwunfv3liHIqImFhFgn1vMbLEh31gmcTVUlUJTGYrMFgCXWzKJyD016uOQFAoFpkyZgr/++gsTJ05szEMRURPSqlUI9fUCcLkrjNzLhfLxYn46NfTaej1shYiaSJM8mzI4OBhvv/12UxyKiJpIJAfxuzWOFyOSjyZ7UDgReZZI/8vjxsj9OMaL8U5KIvfHMEZE9RLFucbc2nm2jBHJhsvCmM1mw/vvv++q3RGRm7vcTckw5o6yDI47KTnhK5G7c1kYKysrw4QJE1y1OyJyc9GB5S1jBo4Zc0cXHNNasJuSyO3V6RabhQsXVvtaWRlvbydqTjhmzL1xAD+RfNQpjC1atAhjx46Fv79/pdesVqvLiiIi93flmDEhBB9/5mbOGzjhK5Fc1CmMdenSBePGjcPIkSMrvWYymbB27VpX1UVEbs7xIW8stcJQYkGAj0biisjBVGZFXnEpAHZTEslBncaMPfroo9W2gGk0GsybN88lRRGR+9NpVAgqD2BZHDfmVrINZgCAVq1EIEMykdurU8vY1KlTq31NpVIxjBE1M5EB3rhkLENWgQkdIisPXyBpOCbijQrQsfuYSAY4zxgR1Vs05xpzS47xYpzwlUgeGhzG2BpG1HxFcK4xt8QJX4nkpcFh7PXXX3dFHUQkQ47B4RcYxtyKIxxzwlcieWhwGBNCuKIOIpIhRxhzdIuRe7jgmNbCXytxJUR0LRocxjg4lKj5cnRTXmAYcytsGSOSFw7gJ6J6Y8uYe+KYMSJ5kU0Y279/P8aMGYPo6Gj4+PigQ4cOWLhwIYxGY53288ILL0ChUKBz586NVClR8+EIY/nGMpjK+BQOd2Cx2nCxyD7PGGffJ5KHOs0zVpWmGDN2+PBh9O3bF+3bt8fSpUsRGhqKnTt3YuHChdi3bx82bdp0TftJSUnBP/7xD0RERDRyxUTNg7+3GjqNEqYyGy4YTIgP0UtdUrOXU1QKq01ApVQg1JdjxojkoMFhrF+/fq6oo0Yff/wxTCYTNmzYgNatWwMABg8ejKysLLzzzju4dOkSgoKCatyHxWLBhAkTMGXKFPz555/Iyclp9LqJPJ1CoUCkvw6nc404X8Aw5g4cE75G+GmhUnJML5EcNLibcuvWra6oo0Yajf1xHgEBARWWBwYGQqlUwsvLq9Z9LF68GHl5eXjllVcapUai5iqC48bciuNmigh2URLJhsvGjG3fvh2ff/658/sLFy7g1ltvRWRkJMaNGweTqf6/qMePH4/AwEA8/vjjSEtLQ2FhITZv3oyVK1di2rRp0Otr/t/44cOHsWjRIrz99tvw9fWtdx1EVFkk76h0K1kcvE8kOy4LYy+99BIOHz7s/H7WrFnYtWsX+vbti/Xr1zdoctiEhAQkJycjNTUVrVu3hr+/P0aNGoXx48dj2bJlNW5rs9kwceJE3HHHHbj11lvrdFyz2QyDwVDhi4gqct5RWWCWuBICLt9JGenPaS2I5MJlYeyvv/5Cjx49ANjHZ23cuBGvvfYavvjiCyxcuBCffPJJvfd9+vRpjBo1CiEhIVi/fj1+/vlnLFmyBGvXrsXkyZNr3PaNN97A8ePHsXTp0jofNykpCQEBAc6v2NjYep4BkedydFOyZcw9OLqLIwM4eJ9ILho8gN/BYDAgMDAQALBv3z4UFxdj9OjRAIDevXtj/vz59d737NmzYTAYkJKS4uyS7N+/P0JDQzFx4kSMGzcOAwYMqLTdmTNn8NJLL2Hx4sXw8vJCfn4+AHtYtNlsyM/Ph1arhbd31f+DnDNnDmbOnFnhHBnIiCpydFNyzJh74ISvRPLjspax8PBwHD9+HIB9UH98fDxiYmIAAIWFhc5B+PWRkpKCxMTESmPDevXqBQBITU2tcru0tDSUlJRgxowZCAoKcn798ssvOHLkCIKCgjBnzpxqj6vVauHv71/hi4gqcg7g5/Mp3QInfCWSH5e1jA0fPhzPP/88Dh06hLVr12L8+PHO144ePYqEhIR67zs6OhqpqakoKiqqMAA/OTkZAJyh72rdu3fHjh07Ki1/+umnUVBQgPfee6/abYno2jhaxrILTbDZBJScTkEyQojL3ZT+DGNEcuGyMPbqq6/izJkzWLVqFXr37o0XXnjB+drHH3+Mvn371nvfTz/9NMaMGYOhQ4fimWeeQWhoKPbs2YOkpCQkJiZixIgRAIBJkyZh3bp1OHnyJOLj4xEYGIiBAwdW2l9gYCAsFkuVrxFR3YT7aaFQAGVWgTxjKScaldAlYxlKLTYAQDgfEk4kGy4LY6Ghofj++++rfG3Hjh3Q6er/v7TRo0dj27ZtWLx4MWbMmIGCggLExsZiypQpmDNnjnOeMavVCqvV2iRPBSAiO41KiRC9FjlFZpwvMDGMScgx4Wuorxe0apXE1RDRtXJZGKuJK8ZaDRo0CIMGDapxnbVr12Lt2rW17uunn35qcD1EdFlkwOUw1rlFQO0bUKNwTvjKLkoiWZHNg8KJyH1FchZ+t8AJX4nkiWGMiBqMc425B+eErwxjRLLCMEZEDRbJ6S3cwuXZ9xnGiOSkzmHMZDLh4MGDMBqNlV775ZdfXFIUEclLBCd+dQuXZ9/nhK9EclKnMJacnIzY2FgMHDgQYWFhWLx4cYXXHVNMEFHzEsluSrfAMWNE8lSnMPZ///d/+Oc//4nc3Fzs27cPX3zxBSZOnAibzT6vDaeUIGqenI9EYjelpC4U8G5KIjmqUxg7fPgwxo0bBwDo0KEDfv75Z2RnZ+Ouu+5CaWlpoxRIRO7P8eFvMFlQUmqVuJrmqdBUhkKzBQAH8BPJTZ3CmL+/PzIyMpzfe3t748svv4ROp8Pw4cOdLWRE1Lz469Tw1tgnGeW4MWk4uoj9dGr4aptkCkkicpE6hbGbb74Z7733XoVlarUaH330EVq3bo2SkhKXFkdE8qBQKNhVKbHzBWYAHC9GJEd1CmP/+c9/MHPmzErLFQoFVq1ahdOnT7uqLiKSmYjyZyFyEL80HI9C4ngxIvmpU1u2l5eX8zmQVYmLi2twQUQkT1Hl0ymwm1Ia53knJZFsNXjS13nz5rmiDiKSuQhO/CopzjFGJF8NDmOvv/66K+ogIpmLZDelpDj7PpF8NTiMcW4xIgKumGuMYUwSnPCVSL4aHMYUCoUr6iAimXM+LJzdlJK44OymZBgjkhs+KJyIXMIRArILzbDZ2GLelExlVuQW2yfeZjclkfwwjBGRS4T5aqFUABabQE6xWepympVsg/3nrVUrEeijkbgaIqorjhkjIpdQq5QI9S0fxF/AMNaUHOP0ogJ0HDpCJEMNDmP9+vVzRR1E5AE4iF8anPCVSN4aHMa2bt3qijqIyAM45xpjGGtSnPCVSN44ZoyIXCaSd1RKwjGtBSd8JZKnOj0OqTZffvklPvroI6Snp8NkqvjLWKFQ4M8//3Tl4YjIzbCbUhrOaS3KJ94lInlxWRh7/fXX8dxzzyEsLAxt2rSBXq931a6JSCacc40xjDUptowRyZvLwthbb72FiRMnYuXKlVCpVK7aLRHJSCSfTykJjhkjkjeXjRnLzc3FAw88wCBG1IxFBti7ydhN2XQsVhsuFtmnEuHs+0Ty5LIw1q9fPxw5csRVuyMiGXJ0UxaaLDCWWiSupnnIKSqF1SagUiqc87wRkby4LIwtXboUb775Jr766iuUlpa6ardEJCN+Og30XvbWcXZVNg3nHGN+WqiUnPCVSI5cNmasTZs2uPnmmzF27FgoFAr4+PhUeF2hUKCgoMBVhyMiNxURoEPaxWKcN5jQKsxX6nI8nuNmiQh2URLJlsvC2KxZs7BixQp0794dHTt2hJeXl6t2TUQyEulvD2O8o7JpZHHwPpHsuSyMrV27Fs899xySkpJctUsikiHHIPIsdlM2CUd3cKQ/p7UgkiuXjRmzWq0YOnSoq3ZHRDLVItAeCjIulUhcSfPguHPVcScrEcmPy8LYLbfcgj179rhqd0QkU84wls8w1hQ44SuR/Lmsm/LFF1/EvffeC71ej5EjRyI4OLjSOlUtIyLP0iKILWNN6XI3JceMEcmVy8JYt27dAAAzZ87EzJkzq1zHarW66nBE5KaubBkTQkCh4HQLjcVmE86pLaIDGcaI5MplYeyll17iL10iQnR5GDOWWpFvLEOQnndWN5bsQjPKrPYJX9kyRiRfLgtj8+fPd9WuiEjGdBoVQn21yCkyIyO/hGGsEWXkGwHYuyjVKpcNASaiJuayq7esrAzFxcVVvlZcXIyysjJXHYqI3Jxj3Ng5jhtrVI6fr+PnTUTy5LIwNnnyZEyePLnK1x577DE8/vjjDdr//v37MWbMGERHR8PHxwcdOnTAwoULYTQaa9zuiy++wP333482bdrA29sbCQkJePDBB3H8+PEG1UNE1YvhHZVNwvHzdfy8iUieXBbGfvrpJ4wePbrK10aNGoVt27bVe9+HDx9G3759cfr0aSxduhSbN2/Gfffdh4ULF+L++++vcdvXXnsNRqMRc+fOxffff49FixZh//796NGjBw4dOlTvmoioeryjsmk4WsZi2DJGJGsuGzN24cIFREVFVflaZGQkzp8/X+99f/zxxzCZTNiwYQNat24NABg8eDCysrLwzjvv4NKlSwgKCqpy26+//hrh4eEVlg0ePBgJCQn417/+hdWrV9e7LiKqmuOOynOXam65pobJYDclkUdwWctYYGAgTpw4UeVrJ06cgJ+fX733rdFoAAABAQGVjqlUKmt8DubVQQwAoqOjERMTg7Nnz9a7JiKqHid+bRqOn2+LQB+JKyGihnBZGBs0aBCSkpKQl5dXYXleXh4WL16MwYMH13vf48ePR2BgIB5//HGkpaWhsLAQmzdvxsqVKzFt2jTo9fo67S8tLQ3p6eno1KlTvWsiouo5uykZxhqNEMLZ8siWMSJ5c+nUFr169ULbtm1x7733okWLFjh37hw+//xzlJWVYcGCBfXed0JCApKTkzF27FhnNyUATJ8+HUuXLq3TviwWCyZNmgRfX18888wzNa5rNpthNpud3xsMhjodi6i5coSDfGMZis0W6LUu+1VD5fKKS2EqswHghK9Ecuey35Dt27fHrl27MHPmTKxatQpWqxUqlQoDBgzAG2+8gfbt29d736dPn8aoUaMQERGB9evXIywsDL/++isWLVqEoqIirFmz5pr2I4TApEmTsGvXLmzYsAGxsbE1rp+UlNSgEEnUXPnrNPDTqVFosiAjvwTtIuo/TIGq5mh1DPfTQqtWSVwNETWES/+72q1bN2zbtg0lJSW4dOkSgoODodM1/H9ss2fPhsFgQEpKirNLsn///ggNDcXEiRMxbtw4DBgwoMZ9CCEwefJkfPjhh1i3bh1uv/32Wo87Z86cCo92MhgMtQY4IrJrEeiNo+cLkXGJYawxcPA+kedolCmbvb29ER0d7ZIgBgApKSlITEysNDasV69eAIDU1NQat3cEsffeew+rV6/GQw89dE3H1Wq18Pf3r/BFRNfGMd3COY4baxTOCV85xxiR7Mni+RnR0dE4dOgQioqKKixPTk4GAMTExFS7rRACjz76KN577z2sXLkSEyZMaNRaicjOeUcl5xprFM4JX4N4JyWR3MkijD399NPIycnB0KFD8dlnn2H79u149dVXMXPmTCQmJmLEiBEAgEmTJkGtViM9Pd257fTp07FmzRpMmDABXbp0wZ49e5xf+/fvl+qUiDwe76hsXHwUEpHnkMUtTqNHj8a2bduwePFizJgxAwUFBYiNjcWUKVMwZ84c5zxjVqsVVqsVQgjntl9//TUA4N1338W7775bYb/x8fE4ffp0k50HUXPiaLHJ4MSvjYKPQiLyHLIIY4B9HrNBgwbVuM7atWuxdu3aCssYtoikcXkWfraMNQbOMUbkOVzWTalUKqFSqar8UqvVCA0NxfDhw7Fjxw5XHZKI3FhcsL1lLLvQDFOZVeJqPIvBVIZCkwUAB/ATeQKXhbGXXnoJ8fHxCA4Oxvjx4zFr1iw8/PDDCA4ORlxcHB566CGcO3cOQ4cOxZYtW1x1WCJyU4E+9rnGAOBMHrsqXclxU0SQj4YT6hJ5AJddxcHBwYiMjMTBgwcrTEFRVFSEoUOHokWLFkhJScHQoUPxyiuvYOjQoa46NBG5IYVCgYQQPQ5mFOB0TjHnGnMhzjFG5Flc1jK2fPlyPPvss5XmAvP19cWzzz6Lt956C2q1GlOnTsUff/zhqsMSkRuLD7F3VabnsmXMlRwtjTF8QDiRR3BZGDt37hw0Gk2Vr6nVapw/fx4AEBUVhbKyMlcdlojcmDOM5RVLXIlnSc+1/zwTQvW1rElEcuCyMNa+fXssW7YMFoulwnKLxYJly5Y5n02ZlZWFsLAwVx2WiNxYfIg9LLBlzLVOl/88E0LYMkbkCVw2ZmzhwoW488470aZNG4wZMwYRERG4cOECvvzyS2RkZGDDhg0AgC1btuCGG25w1WGJyI0llIex07lsGXMlR8uYI+wSkby5LIzdfvvt2Lx5M1566SX8+9//hhACCoUCf/vb37By5UoMGzYMALB69WpXHZKI3Jyj5SbjUglKLTZ4qWXx0A+3Vma1OeduSwhlyxiRJ3DpPdHDhw/H8OHDYTQacenSJQQFBcHHh78siJqrMD8tvDUqlJRZkZFfgpYc49RgmfklsNgEtGolIvx0UpdDRC7gkv+mlpSUoEWLFs5HD/n4+KBFixYMYkTNnEKhuOKOSnZVuoJjvFh8iA+USoXE1RCRK7gkjHl7e6OkpKTStBZERJzewrU4XozI87hsAMeQIUOwdetWV+2OiDxEPAfxu9TpHN5JSeRpXDZm7Pnnn8edd94JnU6HO+64A1FRUVAoKjahBwcHu+pwRCQTjpaxM2wZcwm2jBF5HpeFsZ49ewIA5s+fjwULFlS5jtXKhwUTNTec3sK1HD/HBIYxIo/hsjD20ksvVWoJIyKKC7a3jJ3NK4HVJqDioPN6s9oEzubZp7WIZzclkcdwWRibP3++q3ZFRB4kOtAbGpUCpVYbzhtMaBHIh1vXV1ZBCUqtNmhUCkTz50jkMVw+A2NBQQF++OEHfPTRR7h06ZKrd09EMqNSKhBb3jqWnsOuyoZw3JEaG+zDFkYiD+LSMPbyyy8jOjoaI0aMwLhx43Dq1CkA9jstFy9e7MpDEZGMOMY3neK4sQbheDEiz+SyMPbWW29hwYIFmDRpEr755hsIIZyv3Xbbbfjmm29cdSgikhnHzPtpFxnGGiLd+YBwhjEiT+KyMWMrVqzAzJkzsWTJkkp3TbZt2xbHjx931aGISGZah/kCAE5eLJK4Enk7Xd7Ny2dSEnkWl7WMpaWlOR8GfjU/Pz/k5+e76lBEJDOtw+wtOQxjDeNoGXPcoUpEnsFlYSwgIAAXLlyo8rXTp08jPDzcVYciIplpHW5vGTt3qQSmMs43WB9Wm3COuXO0NBKRZ3Dp45CWLFmC4uLLY0IUCgUsFgvefvvtalvNiMjzhei9EOCtgRCc/LW+zl0yotRig1at5LQWRB7GZWPGFi5ciF69eiExMRFjx46FQqHAihUrsH//fpw5cwafffaZqw5FRDKjUCjQKkyP/WfycTK7GB0i/aUuSXYcXbytwnw5rQWRh3FZy1ibNm3wyy+/oGPHjnjrrbcghMD777+P0NBQ7Nq1C3Fxca46FBHJEAfxN8yJbPvPzTH+jog8h8taxgAgMTER33//PcxmM3JzcxEUFARvbzanExHDWEOdzOZ4MSJP5dIw5qDVahEdHd0YuyYimeIdlQ3j+Lk5boYgIs/RoG7Kzp07Y+PGjde8flZWFqZPn87Z+ImaIUeISLtYXGFSaKqdEAInysNYG7aMEXmcBoWxe+65B+PGjUNcXBzmzJmDH374ARcvXnT+oi0pKUFqaipWr16NUaNGIT4+Hvv27cPo0aNdUjwRyUdcsA/USgWMpVacN5ikLkdW8opLkW8sg0Jx+WkGROQ5GhTGXnrpJfz111+4//77sXr1aowYMQKRkZHQaDTw9vaGr68vunXrhsceewwGgwH//e9/8csvvyAxMdFV9RORTGhUSsSF2CcrdYx/omtzsvwxUi0CveHtpZK4GiJytQaPGYuKisJrr72GRYsW4ddff0VycjIyMzNRUlKC0NBQdOjQAQMHDkRMTIwr6iUiGWsd5ou0i8U4ebEIN7YNlboc2XCMF2vD8WJEHsllA/g1Gg1uvPFG3Hjjja7aJRF5mNZhvtiCCxzEX0eXp7VgGCPyRC6bZ4yIqDateEdlvTjvpGQYI/JIDGNE1GQcYcLR0kPX5nIY4+B9Ik/EMEZETaZthD2MXTCYUWAsk7gaeTCVWXHuUgkAjhkj8lQMY0TUZPx1GkQH6AAAf2UXSlyNPNjnZQMCfTQI1ntJXQ4RNQKGMSJqUu0i/QAAf11gGLsWJ64YL6ZQ8AHhRJ5INmFs//79GDNmDKKjo+Hj44MOHTpg4cKFMBqNtW6bnZ2NRx55BKGhofDx8cENN9yAbdu2NUHVRHS1dhHlYew8w9i1cPycHD83IvI89QpjeXl5zhBktVrx/vvv46OPPmq0R5wcPnwYffv2xenTp7F06VJs3rwZ9913HxYuXIj777+/xm3NZjOGDBmCbdu2YdmyZdi0aRMiIiIwfPhw/Pzzz41SLxFVzxEqjrFl7JocLQ9j7SM4XozIU9VrnrFbbrkFq1evRvfu3TFnzhx899138PLywm+//YalS5e6uETg448/hslkwoYNG9C6dWsAwODBg5GVlYV33nkHly5dQlBQUJXbrlmzBqmpqdi9ezduuOEGAMCgQYPQrVs3zJo1C7/++qvL6yWi6rV3tIxd4B2V18LRnds+0l/iSoiosdSrZez48ePo3r07AODDDz/Ed999h23btuHTTz91ZW1OGo0GABAQEFBheWBgIJRKJby8qh/UunHjRrRv394ZxABArVbjoYcewt69e5GRkdEoNRNR1dqE+0KhsD9vMafILHU5bs1YasGZPHsvRDu2jBF5rHqFMS8vLxQXF+O3335DREQEYmJi4Ofnh+Lixnne3Pjx4xEYGIjHH38caWlpKCwsxObNm7Fy5UpMmzYNen31c++kpqaia9eulZY7lh06dKhRaiaiqnl7qRAXbH9GJceN1czRehjqq0WIr1biaoiosdSrm/L+++/HoEGDUFhYiKlTpwKwD7BPSEhwZW1OCQkJSE5OxtixY53dlAAwffr0WrtFc3NzERwcXGm5Y1lubm6125rNZpjNl//nbjAY6lg5EVWlXYQf0nONOHahEH3b8BmV1XGE1Q6RHLxP5MnqFcaWLVuGH3/8ERqNBoMHDwYAKBSKRhkvBgCnT5/GqFGjEBERgfXr1yMsLAy//vorFi1ahKKiIqxZs6bG7Wu6Hbym15KSkrBgwYJ6101EVWsf4Ycthy9w3FgtHDc58E5KIs9WrzCmUCgwbNiwCst69uzpkoKqMnv2bBgMBqSkpDi7JPv374/Q0FBMnDgR48aNw4ABA6rcNiQkpMrWr7y8PACostXMYc6cOZg5c6bze4PBgNjY2IacChGBc41dq2OOOykjOV6MyJPVK4wBwE8//YT//ve/yMrKQlRUFO69914MGjTIlbU5paSkIDExsdLYsF69egGwjwurLox16dIFBw8erLTcsaxz587VHler1UKr5TgNIldzDEb/63whhBCczLQax3gnJVGzUK8B/P/4xz/wwAMPIDQ0FCNHjkRYWBgeeughvP76666uDwAQHR2NQ4cOoaioYpdGcnIyACAmJqbabceOHYujR49WmMLCYrHgww8/RJ8+fRAdHd0oNRNR9VqF+kKtVKDQbEFWgUnqctxSXnEpLhbax6y25TMpiTyaQtRjptbo6Gjs2LED7du3dy47duwYBg4ciKysLJcWCABfffUVxowZgz59+uCZZ55BaGgo9uzZg6SkJMTFxWH//v3w8vLCpEmTsG7dOpw8eRLx8fEA7IPwe/bsCYPBgMWLFyM8PBxvvfUWvv76a2zdurXaFrWqGAwGBAQEoKCgAP7+/J8qUUMMfeNnHM8uwnsTemFQ+3Cpy3E7ySdzcf+qPYgN9sauWYOlLodI1tz987teLWMajaZSa1R0dHSN8301xOjRo7Ft2zb4+/tjxowZuO2227Bu3TpMmTIFO3fudB7XarXCarVWeBKAVqvFtm3bMGjQIDz11FMYNWoUsrKy8N1339UpiBGRa7UvHzd2NIvjxqrinOw1wv0+OIjIteo0ZiwzMxMA8Pe//x133XUXnn/+ecTExODs2bN47bXXMGvWrEYpErDPml/bmLS1a9di7dq1lZZHRERg3bp1jVQZEdVHYrQ/Nh/IwuEsThlTlaMcvE/UbNQpjMXExEChUDhbnn744YcKr3///feYNm2a66ojIo+VGGVv8TnCMFalY+ftPxdOa0Hk+eoUxmw2W2PVQUTNTGK0PYylXSxCSakV3l4qiStyHzabcLaMdYxiNyWRp6v31BYFBQX4/vvvkZmZiRYtWmD48OFuOSiOiNxTuJ8Oob5a5BSZcexCIbrHBkpdkts4nVsMY6kVWrUSrUKrf9wbEXmGeg3g3717N1q1aoVly5Zh7969WLZsGVq2bIndu3e7uj4i8mCO1rHDmeyqvJJjHF2HKH+oVfX6NU1EMlKvlrEZM2bg7bffxj333ONc9vnnn2P69On4/fffq90uPz8fqampOHjwIA4cOIC33367PocnIg+RGOWPnX9dxOGsAqlLcSuHysNpIrsoiZqFeoWxEydO4M4776yw7I477sBjjz3m/D47OxtbtmzBwYMHnV8ZGRmIi4tD165d0a1bt4ZVTkSy1zHKPjidLWMVOcJYp2iGMaLmoF5hrEuXLnj77bfx5JNPOpetXLkSnTp1cn5/0003QafT4frrr0dKSgratWuHffv2ISwsrOFVE5FHcISNo+cLYbUJqJR8LBJwOZwmMowRNQv1CmP/+c9/MGrUKLzxxhuIjY3F2bNnoVAosGnTJuc6hYWFOHz4MFQqFQoKCvDss8+if//+WLNmDfr27euyEyAi+WoZ6gudRgljqRXpucVoFcY5tbINJuQUmaFUAB35TEqiZqFeYSwxMdH5vMfMzExER0ejT58+0Gg0znUcE8QCQEBAAFatWoXt27fjkUcewYgRI5CUlAQfH5+GnwERyZZKqUD7SH/8eTYfh7MMDGMADpUP3m8Zqud0H0TNRL1v09FoNPj3v/+Nvn374sYbb6wQxKozePBg/Pnnn9BoNOjevXt9D01EHsQxSJ3jxuwOO8eLBUhcCRE1lQbdM/35559XaAEDgIMHD9a4jbe3N/7xj3/gk08+acihichDOMZFcSZ+u0OZ9jtLOV6MqPlw6QQ2VqsV3bt3xx9//FHruj179nTloYlIphwtYwczDM5HrTVnh3knJVGz4/LZBPnLlIjqIjHKHyqlAjlFZpw3mKQuR1KFpjKczjUC4BxjRM0Jp3YmIkl5e6nQNtw+cP/AueY9+euRLPvzKCP9dQjx1UpcDRE1lQaHsYsXL7qiDiJqxrrFBAIADjbzMHbgXD4AoHMLtooRNScNDmOjR49GUFAQBg4ciKeffhoKhQIZGRmuqI2ImokuMfY7B/8sDyPNlaNl0BFOiah5qNc8Yw7ffPMN/vjjD/zxxx/Yt28fdu7cCQAYM2YM9Ho9unfvjp49e6JHjx7o0aNHhRn6iYgcupaHsYMZBRBCQKFonjPxO1rGusYGSloHETWtBoWxESNGYMSIEc7v8/LynMHMEdKWL1/u/OVqtVobXDAReZ72kX7QqBTIN5bh3KUSxAY3vwmhC4yXB+93bcE5xoiakwaFsasFBwfj5ptvxs033+xcZjAYsG/fPuzfv9+VhyIiD6JVq9Axyh8HzhXgwLmCZhnGDmTkAwDign0QpPeSthgialKNfjelv78/Bg0ahJkzZzb2oYhIxrqUtwYdaKbjxhzjxRxdtkTUfHBqCyJyC44Q0lynt/jzbD4ADt4nao4YxojILXRpEQgASM0ogM3W/CaPdt5JycH7RM0OwxgRuYV2Eb7QqpUoNFtwKrdY6nKa1AWDCecNJigVnGOMqDliGCMit6BWKZ3PY3R02TUXjvNtG+4HHy+X3ldFRDLAMEZEbuO6uCAAwB9nLklcSdPi4H2i5o1hjIjcRs94exjbl54vbSFN7E9O9krUrDGMEZHbcISxY+cNKDJbJK6maVhtAiln8gEA1zGMETVLDGNE5DYi/HVoEegNm2g+48aOZxei0GyBj5cKHSL9pC6HiCTAMEZEbqVHeevYH+nNY9zY76ft53ldXCDUKv5KJmqOeOUTkVvpERcIANjXTAbxO0Jnz/KbF4io+WEYIyK30vOKlrHmMPnr744wlhAscSVEJBWGMSJyKx2j/KHTKGEwWZCWUyR1OY0qu9CEM3lGKBT2bkoiap4YxojIrWhUSnQtfz7jPg8fN+boomwf4Qd/nUbiaohIKgxjROR2ejgmf/Xw+cYcYdNx0wIRNU8MY0Tkdhzjxn5Lz5O4ksblGC/2N4YxomaNYYyI3E6vhCAoFEDaxWJkF5qkLqdRmMqsSM2wPwapJ8MYUbPGMEZEbifQxwvtI+wToO495ZmtYwfOFaDMKhDqq0VcsI/U5RCRhGQRxh555BEoFIpqv/bs2VPj9jt27MDQoUMRHh4OX19fdO3aFcuXL4fVam2iMyCiurq+VQgA4Nc0zwxje9JyAQB9WgZDoVBIXA0RSUktdQHX4sUXX8TUqVMrLR81ahS0Wi169epV7bZbt27FsGHD0L9/f6xatQp6vR5fffUVZsyYgZMnT2LZsmWNWToR1dP1rYKxdvdp/HoqV+pSGoUjjF3fivOLETV3sghjrVu3RuvWrSss+/nnn5GTk4MXXngBKpWq2m3Xrl0LjUaDzZs3Q6/XAwBuvvlmHDt2DGvXrmUYI3JTvVvaW8b+ulCE3CIzQny1ElfkOmaL1Xkn5Q2tQySuhoikJotuyqqsWbMGCoUCEydOrHE9jUYDLy8veHt7V1geGBgInU7XmCUSUQME673QLsIXgOeNG0s5kw+zxYZQXy1ah/lKXQ4RSUyWYaygoADr16/HkCFD0LJlyxrXnTp1KkpLSzF9+nRkZmYiPz8fH3zwATZu3IhZs2Y1UcVEVB/OcWMeFsb2lI+Du74Vx4sRkUy6Ka/2ySefoKSkBJMmTap13T59+mD79u24++678eabbwIAVCoVkpKS8H//9381bms2m2E2m53fGwyGhhVORHXSp2UI3k9Od46v8hTJaTkALodNImreZBnG1qxZg5CQEIwdO7bWdfft24exY8eiT58+WLlyJfR6PbZv344XXngBJpMJL774YrXbJiUlYcGCBa4snYjqoHdL++D2YxcKkW8sRaCPl8QVNZypzIo/zuQD4HgxIrKTXRg7cOAAfv/9d8yYMQNabe0DeqdNm4aIiAhs3LjROdB/0KBBUCqVmD9/Ph588EG0atWqym3nzJmDmTNnOr83GAyIjY11zYkQUa3C/LRoE+6LE9lF2JOWi+Gdo6QuqcH2n8lHqcWGcD8tWoXqpS6HiNyA7MaMrVmzBgAwefLka1o/JSUFPXv2rHTHZa9evWCz2XDkyJFqt9VqtfD396/wRURN68Y2oQCAXcdzJK7ENS5PaRHC8WJEBEBmYcxsNuPDDz9E79690blz52vaJjo6Gr///nulCV6Tk5MBADExMS6vk4hcp387exjbefwihBASV9Nwu0/aQyW7KInIQVZh7Msvv0ReXl61rWKTJk2CWq1Genq6c9kzzzyD1NRUjBo1Cps2bcKWLVswe/ZsLFmyBDfffDO6devWVOUTUT30aRkCjUqBs3klSM81Sl1OgxhMZc7xYo4WPyIiWYWxNWvWQK/X47777qvydavVCqvVWuF/z0899RQ2bNiAwsJCTJ48GWPHjsXmzZsxb948fPnll01UORHVl16rdj5Ie9fxixJX0zC7T+TAahNoHaZHLJ9HSUTlFMIT2v2biMFgQEBAAAoKCjh+jKgJvfXTCSz5/hhu7hiB1eP/JnU59TbniwP4ZO9ZTOiXgHmjOkldDlGz4e6f37JqGSOi5ql/2zAAQPLJHJRZbRJXUz9CCPx8zN6yN6BdmMTVEJE7YRgjIreXGOWPYL0Xikut2F8+5kpuTmQXIbPABK1aycleiagChjEicntKpcI54H3nX/IcN/Zzed19WoVAp1HVsjYRNScMY0QkC/3Lu/Z++itb4krqxxHG2EVJRFdjGCMiWRjYPgwKBZCaYUBWQYnU5dRJSanV+bBzhjEiuhrDGBHJQqivFj3i7FNcbD0ir9axX07koNRiQ4tAb7QO4yOQiKgihjEiko2bO0YAALYeviBxJXXz4+HzAIChiRF8BBIRVcIwRkSyMTTRHsaST+aiyGyRuJprY7HanC15t3SKkLgaInJHDGNEJButw/RoGapHqdUmm7sq96VfQl5xKQK8NeidECx1OUTkhhjGiEg2FAoFbu4YDkA+XZU/ltc5pGM41Cr+yiWiyvibgYhkZWhiJABg+7FsWNx8Nn4hhHO82LBOkRJXQ0TuimGMiGSlR1wggnw0yDeWYW/5dBHu6khWIc7mlUCnUTof6UREdDWGMSKSFbVK6Wxl+vpApsTV1MzRKnZT2zB4e3HWfSKqGsMYEcnO6G7RAIBvD55HqcU9uyqFENh8IAsAuyiJqGYMY0QkO31ahSDcT4uCkjLsOu6ed1UeyjTgRHYRtGolp7QgohoxjBGR7KiUCozsGgUA+OpP9+yq/HJ/BgD7RLX+Oo3E1RCRO2MYIyJZcnRVbjl8ASWlVomrqchqE86QOOa6FhJXQ0TujmGMiGSpe2wg4oJ9YCy1YusR95pzLPlkLrILzQj00fDB4ERUK4YxIpIlhUKBUd3sXZWOLkF38WWKvZ6RXaLgpeavWSKqGX9LEJFs3dEjBgCw41g2zheYJK7GzlRmxfep9ikt2EVJRNeCYYyIZKt1mC96twyGTQCf/35W6nIAAJsPZKHIbEFMkDd6xgVJXQ4RyQDDGBHJ2v29YwEA//3tLKw2IXE1wId70gEA9/eOg1KpkLgaIpIDhjEikrURnaMQ4K1BRn6J5HOOpWYUIOVsPjQqBe75W6yktRCRfDCMEZGs6TQqjC0fm/XfvdJ2VX706xkA9hn3w/y0ktZCRPLBMEZEsnd/7zgAwNYjF3DBIM1A/kJTGTaV30X50PXxktRARPLEMEZEstc+0g+9EoJgsQm8+8spSWr4cn8GjKVWtAn3RZ+WwZLUQETyxDBGRB5hSv/WAICP95yBwVTWpMe2WG1YtcseAh/sEweFggP3iejaMYwRkUcY3CEcbcN9UWi2OO9obCpf/ZmJM3lGBOu9cG8vDtwnorphGCMij6BUKjB1gL117N3/nYaprGmeV2m1CazYcQIAMOnGlvDxUjfJcYnIczCMEZHHGN09GtEBOuQUmbF+37kmOea3B7OQdrEYAd4ajLuBA/eJqO4YxojIY2hUSky+qRUA4M0dJ1BS2ritYzabwIrt9laxCf0S4KfTNOrxiMgzMYwRkUd5oE8cWgR6I6vA1Oh3Vn6ZkoFjFwrhp1VjQt+WjXosIvJcDGNE5FF0GhX+Pqw9AODtn04ip8jcKMcpKbViyffHAABPDGqDAB+2ihFR/TCMEZHHGd0tGp1b+KPIbMHybccb5RirdqXhvMGEFoHemNAvoVGOQUTNA8MYEXkcpVKB52/tCMD+iKKD5wpcuv8LBhPe/ukkAGD2iA7QaVQu3T8RNS8MY0Tkkfq2DsWtXSJhtQlM/+9+FJstLtmvzSYwd2MqSsqs6BEXiNu6Rrlkv0TUfMkijD3yyCNQKBTVfu3Zs6fWfWzatAkDBgyAv78/9Ho9OnXqhHfeeacJqiciqbw6tguiAnQ4lVOMeV8dcsk+V+w4ga1HLsBLpcTC2ztztn0iajCFEEJIXURtTp48iYsXL1ZaPmrUKGi1WqSnp0Olqr6bYPHixZg7dy6mTp2KUaNGQaPR4OjRoxBC4Mknn7zmOgwGAwICAlBQUAB/f/96nQsRNa09abl4YNUe2ASw9N7uGHNdi3rva8fRbExc9xuEAJbc2RX3cLZ9Illw989vWUwV3bp1a7Ru3brCsp9//hk5OTl44YUXagxi+/btw9y5c5GUlIRZs2Y5lw8ZMqTR6iUi93F9qxA8Obgtlm87juc2HEB8iA+uiwuq836Onjdg+n/3Qwj78ycZxIjIVWTRTVmVNWvWQKFQYOLEiTWut2LFCmi1Wjz11FNNVBkRuZsZQ9ri5o7hMFtsePT933HukrFO25+7ZMT4d/ei0GRBr4QgzBvVqZEqJaLmSJZhrKCgAOvXr8eQIUPQsmXNEy3u3LkTHTt2xIYNG9C+fXuoVCrExMRg9uzZKC0tbaKKiUhKKqUCy+67Dh2j/JFTVIpH3vsNf10ovKZtc4rMGP/uXlwwmNE23Berxv0NXmpZ/uokIjcly98on3zyCUpKSjBp0qRa183IyMDx48cxffp0TJ8+HVu3bsUjjzyCf/zjH5gwYUKN25rNZhgMhgpfRCRPeq0aa8b/DeF+WpzILsLI5buwbOtxmC1VPzLJahP46Nd0DPnnzzh5sRjRATq8P6k3An28mrhyIvJ0shjAf7VevXrh1KlTyMjIgFarrXFdLy8vlJWV4ZNPPsF9993nXP7MM89g6dKlOH78ONq0aVPltvPnz8eCBQsqLXfXAYBEVLsLBhPmbjyIrUeyAQChvl64v3cchnWKhLeXCiWlVmw7ko3NBzJxPLsIANAxyh//vv86tAn3lbJ0Iqondx/AL7swduDAAXTr1g0zZszA0qVLa10/KioK58+fR15eHoKCLg/a/fHHHzFs2DB8+umnuOeee6rc1mw2w2y+/CgVg8GA2NhYt30ziejaCCHw1Z+ZWPzdUWQVmKpdz0+rxv/d0g4PXR8PtUqWHQlEBPcPY7K4m/JKa9asAQBMnjz5mtbv2rUrzp8/X2m5I4MqldX/gtVqtbW2vBGR/CgUCtzevQVGdonCj4cv4P3k0ziRXQSzxQYIoHfLYAzvHIlbEiP5zEkianSyCmNmsxkffvghevfujc6dO1/TNnfeeSd+/PFHfPfdd3jggQecy7/99lsolUr06tWrscolIjenVilxa5co3NqFs+gTkXRk1e7+5ZdfIi8vr9pWsUmTJkGtViM9Pd25bMKECejRoweeeOIJLF++HFu3bsXs2bPx5ptv4oknnkB8fHxTlU9ERERUiaxaxtasWQO9Xl9hIP6VrFYrrFYrrhwGp9FosGXLFjz//PN49dVXkZeXh5YtW2Lx4sWYOXNmU5VOREREVCXZDeCXkrsPACQiIqLK3P3zW1bdlERERESehmGMiIiISEIMY0REREQSYhgjIiIikhDDGBEREZGEGMaIiIiIJMQwRkRERCQhhjEiIiIiCTGMEREREUmIYYyIiIhIQgxjRERERBKS1YPCpeZ4jKfBYJC4EiIiIrpWjs9td30cN8NYHeTm5gIAYmNjJa6EiIiI6io3NxcBAQFSl1EJw1gdBAcHAwDOnDnjlm9mYzEYDIiNjcXZs2fd8mn3jYXnzfNuDnjePO/moKCgAHFxcc7PcXfDMFYHSqV9iF1AQECz+kfs4O/vz/NuRnjezQvPu3lpruft+Bx3N+5ZFREREVEzwTBGREREJCGGsTrQarWYN28etFqt1KU0KZ43z7s54HnzvJsDnrd7nrdCuOt9nkRERETNAFvGiIiIiCTEMEZEREQkIY8JY//73/9w6623IigoCN7e3mjbti1efvnlCuuUlZXhjTfeQJcuXeDt7Y3AwED07dsXu3fvdq6zdu1aKBSKar8WL15cYZ/Z2dl45JFHEBoaCh8fH9xwww3Ytm3bNdedlpaGO+64A4GBgfD19cXQoUPxxx9/uPV5f/HFF7j//vvRpk0beHt7IyEhAQ8++CCOHz9+TTXPnz+/ymPodDq3Pu+a1j1//vw11S3H93vgwIE1rlvbubvT+w0A58+fx5NPPolWrVrB29sb8fHxmDRpEs6cOVPpuJ5yfdflvD3p+q7LeXvS9V2X8/a06zsrKwuPPPIIwsPDodPp0LVrV6xZs6bK40p9fV/JI+YZ+/jjj/Hwww/jnnvuwfvvvw9fX1+cPHkSmZmZznWsVivGjh2L//3vf5g1axb69u2L4uJi7Nu3D8XFxc71Ro4cieTk5ErHeOmll7BlyxaMHTvWucxsNmPIkCHIz8/HsmXLEB4ejjfffBPDhw/H1q1bMWDAgBrrvnjxIm666SYEBQXh3XffhU6nQ1JSEgYOHIjffvsN7du3d8vzfu211xAZGYm5c+eiVatWOHv2LF599VX06NEDe/bsQadOnWqs2+H777+vMHnutc7/ItV5O7z33nvo0KFDhWUhISG11i3X9/utt96q9Agwo9GI4cOHo2fPnoiMjKz13AH3eL/NZjP69++PS5cuYcGCBUhMTMSxY8cwb948/PDDDzhy5Aj8/Pyc63rK9V2X8/ak67su5+3gCdd3Xc7bk67vgoIC3HjjjSgtLcWSJUsQFRWFTz75BJMnT0ZBQQFmzpxZ4Wck5fVdiZC5c+fOCb1eLx5//PEa1/vXv/4llEqlSE5OrvMxioqKhK+vr7jxxhsrLH/zzTcFALF7927nsrKyMpGYmCh69+5d637//ve/C41GI06fPu1cVlBQIEJDQ8U999xT47ZSnveFCxcqrZuRkSE0Go2YNGlSrfudN2+eACAuXrxY55qkPO/33ntPABC//fZbnfcphHzf76qsXbtWABCrV6+udV13er+3bNlSZd0ff/yxACC++OIL5zJPur7rct6edH3X5bw96fquy3lXRa7Xd1JSkgAgfv/99wrLb7nlFqHX68WlS5ecy6S8vqsi+zA2f/58AaDCD6QqCQkJYuDAgfU6xpo1awQAsXbt2grLb775ZtG+fftK67/66qsCgDh37lyN+23Tpo0YNmxYpeWPPfaY8Pb2FmVlZdVuK+V5V6dly5billtuqXW9hly8Up53Q39Ze9L7fdNNNwlfX19RWFhY67ru9H7/9NNPAoD4/PPPKyz/9ttvBQDx7bffOpd50vVdl/Oujhyv77qctydd3w19v+V6fd92220iIiKi0vLly5cLAOKTTz5xLpPy+q6K7MeM7dy5E8HBwTh69Ci6d+8OtVqN8PBwTJ061dn0evbsWZw+fRpdunTB888/j4iICKjVanTq1Anr1q2r9Rhr1qyBv78/7r777grLU1NT0bVr10rrO5YdOnSo2n2WlJTg5MmT1W5fUlKCtLQ0tzzvqqSlpSE9Pf2auzAAoEuXLlCpVIiIiMC4ceOqHLNzNXc479tuuw0qlQrBwcG44447kJqaWus+Pen9Pn78OHbt2oX77rsPvr6+te7XwR3e7379+qFnz56YP38+fvvtNxQVFeGPP/7A888/jx49euDmm292rutJ13ddzrsqcr2+63PennB9N+T9lvP1XVpaWuU8Yo5lBw4ccC6T8vquUp2imxtq37690Ol0ws/PT7z66qtix44dYsmSJcLb21v069dP2Gw2kZycLAAIf39/kZiYKD777DPxww8/iLvuuksAEO+88061+z9y5IgAIKZMmVLpNY1GU+Xy3bt3CwDi448/rna/GRkZAoBISkqq9JqjKfnK5lN3Ou+rlZWViYEDBwp/f39x5syZWtd///33xSuvvCK+/fZbsX37drF48WIRHBwsIiIiav3fiJTn/d1334m5c+eKr7/+Wvz8889ixYoVIiYmRuj1epGSklJj3Z70fj/33HMCwDV3hbrb+20wGMSoUaMEAOfXwIEDRW5uboX1PO36vtbzvprcr+9rPW9Pu77r+37L+fp++umnhVKpFOnp6RWO8/DDDwsA4rHHHnMuk/L6rorsw1jbtm2r/KEsXbpUABBbtmwRv/zyiwAgvLy8KjSH2mw20aNHDxETE1Pt/p999tlqm641Go2YOnVqpeWON/PKJtGrOd7MxYsXV3rN8WbWdDFIed5XstlsYty4cUKlUokvv/yyxnVr8uuvvwqlUimmT59e43ruct4Op06dEr6+vmL06NE1rucp73dZWZmIjIwUnTp1qnG92kj1fpeWlooRI0aI2NhYsWrVKrFz506xbt060bZtW9GjRw+Rn5/vXNeTru+6nPeV5H591/e8HeR6fdf3vOV+fR8+fFhotVpx4403itTUVJGTkyNWrFghvLy8BIAK17OU13dVZB/Grr/+egFA/PHHHxWWHzt2TAAQr732mjh69KgAILp27Vpp+zlz5ggAVQ5aLS0tFeHh4aJbt25VHjsyMlLcfffdlZZv3rxZABA//PBDtXUbjUahUCjE3//+90qvrVixQgAQx44dq3Z7Kc/bwWaziYkTJwqlUik++OCDGte9Fh06dKh14KQ7nPfVhg8fLsLDw2tcxxPebyGE2LRpkwAg/vWvf9W6bm2keL/ffvvtKkPnyZMnBQAxf/585zJPur7rct4OnnB91+e8rybH67u+5y3361sI+7i42NhYZ2tgbGys+Pe//y0AiJdfftm5npTXd1VkP2asqj5bABDlT3lSKpVo3bo1fHx8al3vaps3b0Z2djYmT55c5bZdunTBwYMHKy13LOvcuXO1dXt7e6NNmzbVbu/t7Y1WrVpVu72U5+3YfvLkyXjvvfewevVqPPTQQ9Wue62EELXeDi31eVe3z9rqlvv77bBmzRp4eXnh4YcfrnXd2kjxfqekpEClUqFHjx4V1mvVqhVCQkIqjA/ypOu7Luft2N4Tru+6nnd965b7++0g9+sbAEaMGIH09HT89ddfOHz4ME6dOuWcmqR///7O9aS8vqs7GVn74YcfBADxyiuvVFj+xhtvCABi165dQggh7r//fqHRaMSpU6ec69hsNtG9e3fRunXrKvc9cuRIodPpRF5eXpWvv/XWWwKA2LNnj3NZWVmZ6NSpk+jTp0+ttc+aNUt4eXlVGIdhMBhEWFiYuPfee2vcVsrzttlsYtKkSUKhUNQ4DqkukpOThVKpFE8//XSN60l53lVJS0sTvr6+YsyYMbWuK9f32yErK0uo1ep63bZ9Nane7wULFlS6ZoW4/D/xK+vxpOu7LuftSdd3Xc67KnK9vutz3p5wfVfFbDaLPn36iO7du1dYLuX1XRXZhzEhhBg1apTQarXi5ZdfFlu2bBFJSUlCp9OJ2267zbnOiRMnRGBgoGjfvr345JNPxDfffCPGjh0rFApFpdt/hbD3CatUKvHAAw9Ue1yTySQ6deokYmNjxUcffSS2bNkixo4dK9Rqtfjpp58qrDt48GChUqkqLMvOzhZRUVGiS5cuYuPGjeLbb78V/fv3F35+fuLIkSNue95PPvmkACAmTpwokpOTK3xd3dxc1Xl37dpVLFmyRHz99ddiy5Yt4pVXXhGBgYEiOjpaZGZmuu15DxkyRCxYsEBs3LhRbNu2TSxdulRER0cLPz8/cfDgwVrPW67vt8PixYsFAPHjjz9Wu467v99nzpwRgYGBokWLFuLtt98W27dvF6tXrxatWrUSer1eHD161LmuJ13fdTlvT7q+63LennR91+W8HTzh+hbC/u93/fr1YseOHWLNmjWiW7duIiQkRKSmplZYT+rr+2oeEcaMRqN47rnnRGxsrFCr1SIuLk7MmTNHmEymCusdPHhQjBw5Uvj5+QmdTieuv/568fXXX1e5z1deeUUAENu3b6/x2OfPnxfjxo0TwcHBzn1u2bKl0noDBgwQVTVEnjhxQowZM0b4+/sLHx8fMWTIELFv3z63Pu/4+Hhnf/zVX/Hx8bWe93333SfatGkj9Hq90Gg0Ij4+XkydOvWaLlwpz/vpp58WiYmJws/PT6jVahEdHS0eeuihKscGeNL77dCuXTuRkJAgbDZbtevI4f0+fvy4ePjhh0VCQoLQarUiLi5O3HvvveLQoUOV1vWk6/taz9vTru9rPW9Pu77r8u9cCM+5vm+//XYRFRUlNBqNiIyMFI888ki185hJeX1fTSFEeacrERERETU52Q/gJyIiIpIzhjEiIiIiCTGMEREREUmIYYyIiIhIQgxjRERERBJiGCMiIiKSEMMYERERkYQYxoiIiIgkxDBGREREJCGGMSIiIiIJMYwRERERSYhhjIioBuvXr0d0dDSCgoKwZMkSqcshIg/EB4UTEVWjsLAQ8fHxWLFiBfz8/DBx4kRs374dXbp0kbo0IvIgaqkLICJyVyUlJRBCoGvXrvD19UVgYCAsFovUZRGRh2E3JRF5FJPJBLVaDT8/Pzz11FM1rrtw4UIkJibCZrNV+Xp4eDjGjBmDLl26oGXLlrj33ntx3XXXOV9fs2YNWrRogeLiYpeeAxE1LwxjRORRFAoFtm/fjt69e2PFihX466+/qlwvMzMTS5YswcKFC6FUVv+r8MqgduONN1Z4bfz48dDr9RxLRkQNwjBGRB5Fq9Wif//+mDVrFgBg3759Va63bNkyBAYG4o477qh2X7/88gs++OADjBw5EgCQkpJS4XW1Wo0pU6Zg2bJlMBqNrjkBImp2GMaIyCN16NABQOUABQClpaVYs2YNHnjggWpbxaxWK5588knExMRg3bp1UKlUVe7rwQcfhMFgwH//+19Xlk9EzQjDGBF5pNdeew1A1WHs119/RW5uLgYNGlTt9v/5z3+QkpKC119/HSEhIWjbtm2V+4qMjESHDh3wzTffuKp0ImpmGMaIyOP8+OOPePvttxEUFIT9+/dXej05ORkA0KNHjyq3v3jxIl588UUMGDAA9957LwCga9euOH78eJXdkT169MAvv/ziwjMgouaEYYyIPEp+fj4mTpyI22+/HU888QQuXryIzMzMCutkZmZCoVAgNDS0yn3Mnj0bBoMBy5cvdy7r2rUrbDYbDhw4UGn98PBwZGdnc9oLIqoXhjEi8ijTpk1DWVkZVq1a5ZyG4uruxZKSEmg0GqhUqkrb7927F++99x4eeughxMXFIT8/H/n5+WjVqlWV+wIAnU4HIQRMJpPLz4eIPB/DGBF5jPXr1+Pjjz/GmjVrEBYW5uyGvDpAhYaGorS0tNL8YDabDdOmTYMQAuvWrUNQUJDz64EHHqhyXwCQl5cHrVYLX1/fRjkvIvJsnIGfiDzC+fPnMXXqVDz22GO47bbbAAAtW7ZEYGBgpXFjjjstT548ia5duzqXr1q1Cr///jsWLFiA/v37VzrGnXfeWWUYS0tLQ2JiogvPhoiaE4YxIvIIjz76KIKCgvDGG29UWH7ddddVClADBw4EAOzZs8cZxvLy8jB37lz07dsXL774IhQKRaVjdOvWDb/++itsNptzSgybzYa9e/di0qRJrj8pImoW2E1JRLK3evVqfPfdd/jggw+g1+srvHbdddfh5MmTKCwsdC6LjY3FTTfdhE2bNjmXPf/88ygoKMDKlSurDGKAPYwZjcYKs/r/9NNPKCgowIMPPujisyKi5kIhhBBSF0FE1NQ2bNiAe++9F+np6WjRokW99/Pwww8jLS2NU1sQUb0xjBFRsySEQN++fdGzZ0+sWLGiXvs4efIkOnbsiO3bt1d6biUR0bViNyURNUsKhQKrVq1CdHR0hYeB18WZM2ewYsUKBjEiahC2jBERERFJiC1jRERERBJiGCMiIiKSEMMYERERkYQYxoiIiIgkxDBGREREJCGGMSIiIiIJMYwRERERSYhhjIiIiEhCDGNEREREEmIYIyIiIpLQ/wPHEk6fVq56twAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the irradiance spectrum\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "plt.plot(lam_c, F_obs)\n",
    "plt.xlim(6670, 6690)\n",
    "\n",
    "plt.xlabel(r'$\\lambda ({\\AA})$')\n",
    "plt.ylabel(r'$F^{\\mathrm{obs}}_{\\lambda}\\ ({\\rm erg\\,cm^{-2}\\,s^{-1}}\\,\\AA^{-1})$');"
   ]
  },
  {
   "cell_type": "raw",
   "id": "4467201c-0aa1-4c3d-8e95-ebcba7edd994",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/restructuredtext",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "This figure shows the characteristically rounded shape due to rotation broadening. Note that the line center is around :math:`6680\\,\\angstrom`, rather than :math:`6678\\,\\angstrom`, because the :file:`sg-HeI-6678.h5` grid adopts vacuum rather than air wavelengths.\n",
    "\n",
    ".. rubric:: Footnotes\n",
    "\n",
    ".. [#specgrid-file] Generated by applying the :command:`subset_specgrid` tool to the high-resolution :mad-star:`BSTAR2006 grid <msg-grids#BSTAR2006_Grids>`."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}