examples/3dq8_20M.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "53cc2cde",
   "metadata": {},
   "source": [
    "# Example usage of the `3dq8_20M` fitting model\n",
    "\n",
    "`3dq8_20M` models ampitudes and phases of the ringdown from a quasi-circular, non-precessing black-hole binary. It is calibrated up to mass ratio 8 and at a starting time $t_{\\rm start}=20M$ from the peak of the $(2,2)$ strain.\n",
    "\n",
    "The model natively fits for the effective cartesian coordinates $\\tilde A^{x,y}$ defined by\n",
    "$$\n",
    "\\tilde A_{\\rm mode}^x+i\\tilde A^y_{\\rm mode}=\\tilde A_{\\rm mode}e^{i\\beta\\tilde\\phi_{\\rm mode}}\n",
    "$$\n",
    "\n",
    "where $\\beta=1$ if we are computing predictions for a mode corresponding to spherical-harmonic indices $(l,m)$ with even $m$, while $\\beta=2$ for odd $m$.\n",
    "\n",
    "The amplitudes $\\tilde A$ are the relative amplitudes $\\tilde A_{\\rm mode}=A_{\\rm mode}/A_{220}$, with the only exception that $\\tilde A_{220}\\equiv A_{220}$.\n",
    "\n",
    "The relative phases are defined as $\\tilde\\phi_{\\rm mode}=\\phi_{\\rm mode}-\\frac{m}{2}\\phi_{220}$.\n",
    "\n",
    "They are extracted from $\\tilde A_{\\rm mode}^{x,y}$ as\n",
    "$$\n",
    "\\tilde A_{\\rm mode}={\\rm abs}(\\tilde A_{\\rm mode}^x+i\\tilde A_{\\rm mode}^y)\n",
    "$$\n",
    "$$\n",
    "\\tilde\\phi_{\\rm mode}=\\frac{{\\rm angle}(\\tilde A_{\\rm mode}^x+i\\tilde A_{\\rm mode}^y)}{\\beta}~{\\rm rad}\n",
    "$$\n",
    "meaning that $\\tilde\\phi_{\\rm mode}$ has support in the interval $[-\\pi,\\pi]$.\n",
    "\n",
    "1-$\\sigma$ uncertainties are propagated as\n",
    "$$\n",
    "\\sigma_{\\tilde A}=\\frac{\\sqrt{\\left(\\sigma^x\\tilde A^x\\right)^2+\\left(\\sigma^y\\tilde A^y\\right)^2}}{\\tilde A}\n",
    "$$\n",
    "$$\n",
    "\\sigma_{\\tilde \\phi}=\\frac{\\sqrt{\\left(\\sigma^y\\tilde A^x\\right)^2+\\left(\\sigma^x\\tilde A^y\\right)^2}}{\\beta\\tilde A^2}~{\\rm rad}\n",
    "$$\n",
    "\n",
    "__Note:__ the model does not return predictions for mixing modes, i.e. modes $(l',m,n)$ contributing to a spherical harmonic $(l,m)$ with $l'\\neq l$. Rather, it is assumed that they are be derived from the non-mixing modes $(l',m,n)$ and from the corresponding spherical-spheroidal mixing coefficients $\\mu_{mll'n}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "73e3c2fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('..')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "46f904d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "38e14dd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import postmerger as pm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be3f220e",
   "metadata": {},
   "source": [
    "## Load fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a185475c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3dq8_20M models ampitudes and phases of the ringdown from a quasi-circular, non-precessing black-hole binary.\n",
       "It is calibrated up to mass ratio 8 and at a starting time 20M from the peak of the (2,2) strain."
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fitname = '3dq8_20M'\n",
    "fit = pm.load_fit(fitname)\n",
    "fit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "936a88cb",
   "metadata": {},
   "source": [
    "## Read documentation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "74861d09",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on AmplitudeFit3dq8 in module postmerger._fits object:\n",
      "\n",
      "class AmplitudeFit3dq8(builtins.object)\n",
      " |  AmplitudeFit3dq8(fit_dict)\n",
      " |  \n",
      " |  Methods defined here:\n",
      " |  \n",
      " |  __init__(self, fit_dict)\n",
      " |      Initialize self.  See help(type(self)) for accurate signature.\n",
      " |  \n",
      " |  __print__(self)\n",
      " |  \n",
      " |  __repr__(self)\n",
      " |      Return repr(self).\n",
      " |  \n",
      " |  predict_amp(self, mass_ratio, chi1z, chi2z, lm, mode, return_std=False, start_time=None)\n",
      " |      Predict the value of abs(A^x+iA^y) corresponding to the query points (mass_ratio,chi1z,chi2z).\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      mass_ratio : array_like of shape (n_samples,) or float\n",
      " |          Mass ratio of the query points.\n",
      " |          \n",
      " |      chi1z : array_like of shape (n_samples,) or float\n",
      " |          Projection along the z axis of the primary spin.\n",
      " |      \n",
      " |      chi2z : array_like of shape (n_samples,) or float\n",
      " |          Projection along the z axis of the secondary spin.\n",
      " |      \n",
      " |      lm : tuple_like object\n",
      " |          Ordered couple (l,m) specifying the angular number l and azimuthal number m of the harmonic.\n",
      " |      \n",
      " |      mode : tuple_like object\n",
      " |          Ordered tuple specifying the queried quasi-normal mode.\n",
      " |          For linear modes, ordered triple (l,m,n).\n",
      " |          For quadratic modes, ordered couple of the form ((l1,m1,n1),(l2,m2,n2)).\n",
      " |      \n",
      " |      return_std : bool. Default=False.\n",
      " |          Whether or not to return the standard deviation of the predictive distribution at the query points.\n",
      " |      \n",
      " |      start_time : float or None. Default=None.\n",
      " |          Start time of the ringdown waveform, relative to the peak of |h_{22}|.\n",
      " |          If None (default), assume that start_time=self.t0.\n",
      " |  \n",
      " |  predict_phase(self, mass_ratio, chi1z, chi2z, lm, mode, return_std=False, start_time=None)\n",
      " |      Predict the value of angle(A^x+iA^y)/beta corresponding to the query points (mass_ratio,chi1z,chi2z).\n",
      " |      beta is a correction factor introduced in CITE-THE-PAPER to help fitting of A^x and A^y.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      mass_ratio : array_like of shape (n_samples,) or float\n",
      " |          Mass ratio of the query points.\n",
      " |      \n",
      " |      chi1z : array_like of shape (n_samples,) or float\n",
      " |          Projection along the z axis of the primary spin.\n",
      " |      \n",
      " |      chi2z : array_like of shape (n_samples,) or float\n",
      " |          Projection along the z axis of the secondary spin.\n",
      " |      \n",
      " |      lm : tuple_like object\n",
      " |          Ordered couple (l,m) specifying the angular number l and azimuthal number m of the harmonic.\n",
      " |      \n",
      " |      mode : tuple_like object\n",
      " |          Ordered tuple specifying the queried quasi-normal mode.\n",
      " |          For linear modes, ordered triple (l,m,n).\n",
      " |          For quadratic modes, ordered couple of the form ((l1,m1,n1),(l2,m2,n2)).\n",
      " |      \n",
      " |      return_std : bool. Default=False.\n",
      " |          Whether or not to return the standard deviation of the predictive distribution at the query points.\n",
      " |      \n",
      " |      start_time : float or None. Default=None.\n",
      " |          Start time of the ringdown waveform, relative to the peak of |h_{22}|.\n",
      " |          If None (default), assume that start_time=self.t0.\n",
      " |  \n",
      " |  predict_xy_amp(self, mass_ratio, chi1z, chi2z, lm, mode, return_std=False)\n",
      " |      Predict the values of A^x and A^y corresponding to the query points (mass_ratio,chi1z,chi2z).\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      mass_ratio : array_like of shape (n_samples,) or float\n",
      " |          Mass ratio of the query points.\n",
      " |      \n",
      " |      chi1z : array_like of shape (n_samples,) or float\n",
      " |          Projection along the z axis of the primary spin.\n",
      " |      \n",
      " |      chi2z : array_like of shape (n_samples,) or float\n",
      " |          Projection along the z axis of the secondary spin.\n",
      " |      \n",
      " |      lm : tuple_like object\n",
      " |          Ordered couple (l,m) specifying the angular number l and azimuthal number m of the harmonic.\n",
      " |      \n",
      " |      mode : tuple_like object\n",
      " |          Ordered tuple specifying the queried quasi-normal mode.\n",
      " |          For linear modes, ordered triple (l,m,n).\n",
      " |          For quadratic modes, ordered couple of the form ((l1,m1,n1),(l2,m2,n2)).\n",
      " |      \n",
      " |      return_std : bool. Default=False.\n",
      " |          Whether or not to return the standard deviation of the predictive distribution at the query points.\n",
      " |  \n",
      " |  sample_amp(self, mass_ratio, chi1z, chi2z, lm, mode, n_samples=1, start_time=None)\n",
      " |      Draw samples of abs(A^x+iA^y) from the predictive distribution corresponding to the query points (mass_ratio,chi1z,chi2z).\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      mass_ratio : array_like of shape (n_samples_X,) or float\n",
      " |          Mass ratio of the query points.\n",
      " |      \n",
      " |      chi1z : array_like of shape (n_samples_X,) or float\n",
      " |          Projection along the z axis of the primary spin.\n",
      " |      \n",
      " |      chi2z : array_like of shape (n_samples_X,) or float\n",
      " |          Projection along the z axis of the secondary spin.\n",
      " |      \n",
      " |      lm : tuple-like object\n",
      " |          Ordered couple (l,m) specifying the angular number l and azimuthal number m of the harmonic.\n",
      " |      \n",
      " |      mode : tuple-like object\n",
      " |          Ordered tuple specifying the queried quasi-normal mode.\n",
      " |          For linear modes, ordered triple (l,m,n).\n",
      " |          For quadratic modes, ordered couple of the form ((l1,m1,n1),(l2,m2,n2)).\n",
      " |      \n",
      " |      n_samples : int. Default=1.\n",
      " |          Number of samples to be drawn per query point.\n",
      " |      \n",
      " |      start_time : float or None. Default=None.\n",
      " |          Start time of the ringdown waveform, relative to the peak of |h_{22}|.\n",
      " |          If None (default), assume that start_time=self.t0.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      amp : array-like of shape (n_samples_X, n_samples)\n",
      " |  \n",
      " |  sample_phase(self, mass_ratio, chi1z, chi2z, lm, mode, n_samples=1, start_time=None)\n",
      " |      Draw samples of angle(A^x+iA^y)/beta from the predictive distribution corresponding to the query points (mass_ratio,chi1z,chi2z).\n",
      " |      beta is a correction factor introduced in CITE-THE-PAPER to help fitting of A^x and A^y.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      mass_ratio : array_like of shape (n_samples,) or float\n",
      " |          Mass ratio of the query points.\n",
      " |      \n",
      " |      chi1z : array_like of shape (n_samples,) or float\n",
      " |          Projection along the z axis of the primary spin.\n",
      " |      \n",
      " |      chi2z : array_like of shape (n_samples,) or float\n",
      " |          Projection along the z axis of the secondary spin.\n",
      " |      \n",
      " |      lm : tuple-like object\n",
      " |          Ordered couple (l,m) specifying the angular number l and azimuthal number m of the harmonic.\n",
      " |      \n",
      " |      mode : tuple-like object\n",
      " |          Ordered tuple specifying the queried quasi-normal mode.\n",
      " |          For linear modes, ordered triple (l,m,n).\n",
      " |          For quadratic modes, ordered couple of the form ((l1,m1,n1),(l2,m2,n2)).\n",
      " |      \n",
      " |      n_samples : int. Default=1.\n",
      " |          Number of samples to be drawn per query point.\n",
      " |      \n",
      " |      start_time : float or None. Default=None\n",
      " |          Start time of the ringdown waveform, relative to the peak of |h_{22}|.\n",
      " |          If None (default), assume that start_time=self.t0.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      phase : array-like of shape (n_samples_X, n_samples)\n",
      " |  \n",
      " |  sample_xy_amp(self, mass_ratio, chi1z, chi2z, lm, mode, n_samples=1)\n",
      " |      Draw samples of of A^x and A^y from the predictive distribution corresponding to the query points (mass_ratio,chi1z,chi2z).\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      mass_ratio : array_like of shape (n_samples_X,) or float\n",
      " |          Mass ratio of the query points.\n",
      " |      \n",
      " |      chi1z : array_like of shape (n_samples_X,) or float\n",
      " |          Projection along the z axis of the primary spin.\n",
      " |      \n",
      " |      chi2z : array_like of shape (n_samples_X,) or float\n",
      " |          Projection along the z axis of the secondary spin.\n",
      " |      \n",
      " |      lm : tuple-like object\n",
      " |          Ordered couple (l,m) specifying the angular number l and azimuthal number m of the harmonic.\n",
      " |      \n",
      " |      mode : tuple-like object\n",
      " |          Ordered tuple specifying the queried quasi-normal mode.\n",
      " |          For linear modes, ordered triple (l,m,n).\n",
      " |          For quadratic modes, ordered couple of the form ((l1,m1,n1),(l2,m2,n2)).\n",
      " |      \n",
      " |      n_samples : int. Default=1.\n",
      " |          Number of samples to be drawn per query point.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      amp_x : array-like of shape (n_samples_X, n_samples)\n",
      " |      amp_y : array-like of shape (n_samples_X, n_samples)\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors defined here:\n",
      " |  \n",
      " |  __dict__\n",
      " |      dictionary for instance variables (if defined)\n",
      " |  \n",
      " |  __weakref__\n",
      " |      list of weak references to the object (if defined)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97c17f07",
   "metadata": {},
   "source": [
    "## Predict relative amplitudes and phases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "cb512da4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{(2, 2): [(2, 2, 0), (2, 2, 1)],\n",
       " (2, 1): [(2, 1, 0), (2, 1, 1)],\n",
       " (3, 3): [(3, 3, 0), (3, 3, 1)],\n",
       " (3, 2): [(3, 2, 0), (3, 2, 1)],\n",
       " (4, 4): [(4, 4, 0), (4, 4, 1), ((2, 2, 0), (2, 2, 0))],\n",
       " (4, 3): [(4, 3, 0), (4, 3, 1)],\n",
       " (5, 5): [(5, 5, 0), (5, 5, 1), ((2, 2, 0), (3, 3, 0))]}"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## inspect available modes\n",
    "## syntax: {(l,m): lm_modes}\n",
    "fit.modes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4b2a36f6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "relative amplitude 330:\t  0.115(+/-)0.018\n",
      "relative phase 330:\t  1.256(+/-)0.024\n"
     ]
    }
   ],
   "source": [
    "## binary parameters\n",
    "mass_ratio = 1.5 # defined to be >=1\n",
    "chi1z = 0.1\n",
    "chi2z = 0.\n",
    "\n",
    "## spherical-harmonic indices\n",
    "lm = (3,3)\n",
    "\n",
    "## mode indices\n",
    "mode = (3,3,0)\n",
    "\n",
    "## predict\n",
    "amp, sigma_amp = fit.predict_amp(mass_ratio,chi1z,chi2z,lm,mode,return_std=True)\n",
    "phi, sigma_phi = fit.predict_phase(mass_ratio,chi1z,chi2z,lm,mode,return_std=True)\n",
    "\n",
    "## results\n",
    "print('relative amplitude 330:\\t  %.3f(+/-)%.3f'%(amp,sigma_amp))\n",
    "print('relative phase 330:\\t  %.3f(+/-)%.3f'%(phi,sigma_phi))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a7f2dd6",
   "metadata": {},
   "source": [
    "You can convert relative amplitudes to absolute ones by simple multiplication to $A_{220}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fa3baa8d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Absolute amplitude 220:\t  0.173(+/-)0.001\n",
      "Absolute amplitude 330:\t  0.020(+/-)0.018\n"
     ]
    }
   ],
   "source": [
    "amp, sigma_amp = fit.predict_amp(mass_ratio,chi1z,chi2z,lm,mode,return_std=True)\n",
    "amp_220, sigma_amp_220 = fit.predict_amp(mass_ratio,chi1z,chi2z,(2,2),(2,2,0),return_std=True)\n",
    "amp_abs = amp*amp_220\n",
    "sigma_amp_abs = np.sqrt(sigma_amp**2+sigma_amp_220**2)\n",
    "print('Absolute amplitude 220:\\t  %.3f(+/-)%.3f'%(amp_220,sigma_amp_220))\n",
    "print('Absolute amplitude 330:\\t  %.3f(+/-)%.3f'%(amp_abs,sigma_amp_abs))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f4e1a9f",
   "metadata": {},
   "source": [
    "You can also predict $\\tilde A^{x,y}$ directly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a4152c3f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "relative amplitude x:\t  -0.093(+/-)0.022\n",
      "relative amplitude y:\t  0.068(+/-)0.005\n"
     ]
    }
   ],
   "source": [
    "## predict x and y\n",
    "amp_x, amp_y, sigma_x, sigma_y = fit.predict_xy_amp(mass_ratio,chi1z,chi2z,lm,mode,return_std=True)\n",
    "\n",
    "## results\n",
    "print('relative amplitude x:\\t  %.3f(+/-)%.3f'%(amp_x,sigma_x))\n",
    "print('relative amplitude y:\\t  %.3f(+/-)%.3f'%(amp_y,sigma_y))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1eb7734",
   "metadata": {},
   "source": [
    "By default, predicitions are returned at a starting time $20M$ after the peak. You can retrodict (or postdict) the predictions to a different starting time using the `start_time` argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "71996b18",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "relative amplitude 330:\t  0.117(+/-)0.019\n",
      "relative phase 330:\t  -0.566(+/-)0.024\n"
     ]
    }
   ],
   "source": [
    "## start time in units of the initial total mass M\n",
    "start_time = 10\n",
    "\n",
    "## predict\n",
    "amp, sigma_amp = fit.predict_amp(mass_ratio,chi1z,chi2z,lm,mode,return_std=True,start_time=start_time)\n",
    "phi, sigma_phi = fit.predict_phase(mass_ratio,chi1z,chi2z,lm,mode,return_std=True,start_time=start_time)\n",
    "\n",
    "## results\n",
    "print('relative amplitude 330:\\t  %.3f(+/-)%.3f'%(amp,sigma_amp))\n",
    "print('relative phase 330:\\t  %.3f(+/-)%.3f'%(phi,sigma_phi))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84f01c1f",
   "metadata": {},
   "source": [
    "## Batch evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "2a8b39a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "## binary parameters\n",
    "mass_ratio = np.linspace(1,8,100)\n",
    "chi1z = 0.\n",
    "chi2z = 0.\n",
    "\n",
    "## spherical-harmonic indices\n",
    "lm = (3,3)\n",
    "\n",
    "## mode indices\n",
    "mode = (3,3,0)\n",
    "\n",
    "## predict\n",
    "amp, sigma_amp = fit.predict_amp(mass_ratio,chi1z,chi2z,lm,mode,return_std=True)\n",
    "phi, sigma_phi = fit.predict_phase(mass_ratio,chi1z,chi2z,lm,mode,return_std=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "dc18a8b3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEJCAYAAAC3yAEAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABGuElEQVR4nO3dd3zV9fX48de5K3uQBWFPWQqoEVx1oChSLY5vnbXaRdVaR2ur32qto7a23w47HKVqf7V1tVrrKI66cNQqoMgQZQkSVkjInvfez/n98bkJAZJwAze5N7nn+XjkkdzPPPdCPifvLaqKMcYYEyueeAdgjDGmf7HEYowxJqYssRhjjIkpSyzGGGNiyhKLMcaYmPLFO4B4Kygo0JEjR8Y7DGOM6VOWLFlSrqqFHe1L+sQycuRIFi9eHO8wjDGmTxGRjZ3ts6owY4wxMWWJxRhjTExZYjHGGBNTSZtYROQMEZlfXV0d71CMMaZfSdrEoqrPquq8nJyceIdijDH9StImFmOMMT3DEosxxpiYssRijDEmppJ2gKSInAGcMXbs2HiHYowxMeE4SjAMwbASDLV+37WtOag0tShNQaUpCKOLvIwe5I95HEmbWFT1WeDZkpKSb8Q7FmOMaS8UjiSDkNIciiSJkNISdhNDY+SrKag0t34PKiEHpIPrqYLi7vN6wOOBphYlI4AlFmOM6WscR2kJQXNIaQkqLaHI66BDQzPUNztukmgrSSihMHhaM4TgZgXA0V3JoTVBeEXweiAzVfAIiHSUWvam6vTE2wUssRhjTLeouqWJ5kgpoTmkNAehqcWhvllpaFIampWGSKmiJaSI7F6SUNxShFd2JQmvx00Q2aniHh9lgkhElliMMQZw1E0Q7dshGlsc6pqU+ialvlmpb3JLF60lBwAkUtWk4NsjSWSkQFaq9OkksT+SNrH09cb7O+64g9LSUu69914AKisrKSoqora2ltTU1DhHZ0xiCYW1rbqpMag0NjvUNUFto0Ntk1Lf6NAYdOubditZ6K4Shc8r+DyQkyZ4PMmVKLoraRPL/jTen3DCCXttO/fcc7niiitoaGhgzpw5e+2/9NJLufTSSykvL+d//ud/9tp/+eWXc95557Fp0yaGDRsWdfzLly/nxBNPbHu9dOlSxo8fb0nFJJ3WkkZDpK2isUWpaVBqGhxqGt0SR0vQrY7adY5bDeXzugnD74O0QPKVLHpKwiUWEZkN/AbwAver6p177J8L3A44QAi4RlXfiuzbANQCYSCkqiW9GHqvWr58Oddee23b66VLlzJ16lQA3n77bW666SYcx22cu+qqqzjnnHPYsGEDRxxxBJMnT6ahoYG7776bwsLCtm319fXcfvvtzJ49Oy7vyZiOqCrNIWhocmhoVuqbHaoalOp6h5oGN3G09nhyj3cbvn1e8HuFND9kpljSaFVXU8nO8m20aApjB/VMjU1CJRYR8QJ3A7OAUmCRiDyjqh+1O+wV4BlVVRGZAvwNmNBu/4mqWt4T8b3++uud7ktPT+9yf0FBQZf7u1NaaWlpYd26dRxyyCFt2z788EOmTZtGRUUFV1xxBS+88ALFxcUEg8HdFjI7/vjjeeKJJ3j33Xe58cYbmT9/ftu20tJSzjrrLEssptc5qjQ2t7ZjKLWNDjvrHKrrHaobHcJh2kocTqQtw+8Fv0/IzRA8SZg0VJX62iqqq8qpqSxn/CEz8Hg8vP/OSyx971Vqqircr+oKGupquPvx9xERHrjrBl7911855pQLOOrQB3sktoRKLMB0YK2qrgcQkceAuUBbYlHVunbHZ9DWES95fPTRRwwZMoT09HTA/Q/2+uuvc+GFF7JgwQLOOussiouLAfD7/Rx11FF7XWPatGls2rRpt21VVVWoJt3HaXpJa8mjrtGhvkmpbnCTx846h9pGB8dxk4dGGsb9PjdxZKUK3iRq06irrWLrpnVU7SyjemcZVZU7qNpZxnlfvYGsnDyef/KPPP7AT6mpqiAcDrWd95cXN5KdW8DKpf/hhX88QM6AArJy8snKyaN46GhCoSB+f4BT5l7KoUeeTM7AnmtfTrTEMgRo/7QrBWbseZCInAX8FCgCPt9ulwIviYgCf1DV+R3dRETmAfMAhg8fHpvIe9Hy5cspKytj3bp1DB48mDvuuIONGzcycuRIPvjgg7ak8uqrr3LbbbeRnZ3NM888s9s1Fi5cyIQJE9p+Puqoo1i+fDlPP/10r78f07+0lj5qG5XaJoeKmjAVdUpVvdPW1tH650vACwGf9OsGcVWltmYnKSnppKSmsXXTOv7z2j+pLN9GZUUZVTu3U1mxnet+/P8YfdBU/vPqP7n7J9/a7Rqp6ZmcMvcrZOXkUVQ8nMOPmU3ugEKyBxSQk1tA9oACUtMyAbjomzdz8eW3dBrPxKlHMRGobkiecSwdDhrda4PqU8BTInIcbnvLyZFdx6jqFhEpAv4tIh+r6hsdnD8fmA9QUlLS5/5EX758OaeeeiqnnXYa9fX13HDDDYwePZo77riDk08+mTVr1gAwc+ZMZs6cSUnJrqamhQsXcsIJJ5CZmcldd90F7Koee/jhh3n11Vc56aST4vG2TB+j6o7XqG1UahrClNc6lNc6VNU7hCMjwFurrQI+t60jK7V/TU/Y1FhPxY4tVJRtpnjoGAoHDaN042r+eu+t7CzfSmX5VnaWbyMUbOEHP3+cGcefztbSdTx0982kpWcxIH8gufkDGTFmMj6fOwJ+2vSZ3PiLv5ObV0RuXhE5AwpISU1vu+fhR5/K4Uef2mlMHk/8P+NESyylQPvGhqHAls4OVtU3RGSMiBSoarmqbolsLxORp3Cr1vZKLH3d8uXL+frXv84TTzzRtu3b3/42ABUVFZxwwglcdtllDB48mFAotNu5rUmk1YYNG9p+vuiii5gxYwbXXXcdAwYM6Nk3YfqUYDjSy6pBqagNs6PaoaLOIRhWd2C4uo3lAZ+Qndo/Sh8tzU1UlG1mx/ZN7Ni2ieGjJzFu0uHs2LaJ2797DhVlW6irqWw7ft51v+TzX7wMVNn06cfkFQxk0rRjyCssJi9/EMNHTwTgkJITePz1MlLTMjq8b1HxcIqK+15NSnuJllgWAeNEZBSwGTgfuLD9ASIyFlgXabw/DAgAFSKSAXhUtTby8ynAbb0bfu9Yvnw5EydO7HBffn4+9913HxdeeCEigsfj4Zprron62pdccgnz58/n+uuvj1G0pq9palGqGhyq6sOUVTmU1YSpbdC29g+PQIofMlIEbwL8dby/GuprKdu6kbKtn1G2dSODhoyi5JjZNDc1MO+syVTtLNvt+LMvvpZxkw4nIyuXgcUjmTztWPKLBpNfOJj8osEMH+X+Tg4dOZ67H1/S6X39/gB+f6BH31u8SaI11orIHOAu3O7GD6rqHSJyGYCq3ici1wNfBoJAI/A9VX1LREYDT0Uu4wMeUdU79nW/kpISbd9rKtFVVlYycOBA6uvr8ftjP3mcSS5NLW7bR2VdmK1VDmVVYRpbdiURvxdS/ELA1/emGAmHQpRt+4ztmz9l2+YNZGTl8LlZ7liyeWdNZvuWDbsdf8Ls87n21gcA+MP/fYfcvCIKBw2lYOAwCgcOJb9oCIGU/jNOrLrBYewgH0dN2L/3JCJLOhvSkXCJpbe0G3n/jdY2CWP6s2BIqax32FkXZlulw/aqMA3NuxrTA15I9Qs+b99JIsFgC9u3bGDLZ2sJtjRxzElnA3DrNWfy4Xuv7dZr6uDDPscd974AwOMP3InX52Pg4JEMHDySokHDyMkr6jPvOxYssfSgvlZiMSYajrqN6pV1DtsqQ2ytdBvVW0siAZ9bEvH3gSSiquzcsZUd2zcx4RC3k+hf7r2Ft15+krItG9oGAhcOGsb9T38MwD/+8mvq66opHjKKgUNGMWjoaPIKivF6vXF7H4lCVQk7UN2gTBjSM4kl0dpYjDH7obU0UlETpnSnw/aqEMEwEJnrKjUg5GUm9ujzcDiMx+NBRFj81gu8+fITlG74hNKNq2lqqMPr9fG3N8rx+fxkZGYzdsJhHHfKuQweNobBw8ZSPGx027XOvvjaLu7U/7Qmi5AD4bC7Lkso7O5rnVm5dUZlBFJ9Ql6mh8Kcnkm0lliM6YOaWtzeWdurHUorQuysdf9qV3Ub1tMDgs+buEmkunIHHy9/j8/Wf8Rn6z9i0/pVlG5czR/+sZL8wmI2rFvB8sULGTpqPCd9/ksMHTmeoSMPQiIjEs6++Dtxfge9Q3VXkgiF3eQRdnbNQgC7xmOk+oT0VCEj00tGKmSmekgLCKl+ISXyFYi0l/X0TAWWWIzpAxqaHSpqHbZWhtlUHqKmYdekimkBYUBGYpZGGupr2bh2BZ+uWc6GtSs447wrGDZqAu+/82/uutWd/7Vw0DCGj5rI1Okz297TOV/+Lv9zyXVxjLznOaptCSMUxi1hsitptCaMNL+QkSpkpnrJTHNnIkgLeEjxC6mBSNLwkVBdvJM2sfT1afNN/9bYolTUhNlSGeaz8hC1je5jxgOkpyRmtVZlxTZEPOTmFbFh7Qp+ev0FbCtd37Y/M3sA04+dw7BREzjsqFn87P5XGT5qIumZ2XtdK9HeW3e1Vk0Fw7uWGVbdPWkIkJEq5GR4yUoTstM8ZLSVMtzqy1R/3xwTZI331nhvEkAwpJTXOGytCrOxLERVvTt0vTWRJFp331AoyAf/fZm1q95n7ar3WffJUirLt3HRN2/m3K9eT011BffeeTWjDprC6HFTGDHuYAqKhiTUezgQ7RNHMKQEHaB94lC3h11WuoesNCEn8j0t4CaOtICQ4k+sf9Pusl5hXbDEYuLBcdzG9m1VYTaWhSmrDrftS0+wh05DXQ1rP36f1SsXkzOgkFlfuIRwKMT5MwcSCrYwZMR4xkyYxpjx05g2YybDR0+Kd8gx4ThuSSMYVoKhyHrz7RNHQMhO95CdJuRmeMhK85CeIm2JI5HbuGLBeoUZkwAamh22VzlsKg/xWXmIUMitEkkLkJBTv//57h+y+K3n2fTpx22zXh898yxmfeESvD4fP7//NYqHjel0apK+oDV5tISUltDujeIi7txm+dlecjM85KQLGSlu8khP6f+J40BYYumDmpubueKKK3j55ZfZuXMnY8eO5Sc/+QmnnXZavEMz7TiOUlHrsKUyzKfbQlRGZpP1e9zqLV9a/B9MzU0NfLJiER8v+y+rlr1D9c4d/OqhtwGoqaqgcNBwjjnpHA46uISDJpWQmb1rDrlRB02JV9jdopFG8tbk4bBrtlsRyE53u93mZXrITveQkeI2lqcGEi/Z9xVJm1j2p/H+n+82sLOu56aazsv0cOaM9H0eFwqFGDZsGAsXLmT48OEsWLCAc889l+XLlzNy5Mgei8/sW1NQKasKs6EsUiqJ1HClB4S8BOi5VVdTSUZWLiLC3/70cx774x1to9OHj57IxClHEQ6F8Pp8fPvGe+Iaa3epuoljz9KHqttIXpDtbSt9ZKYKGSlCWoolj56QtIllf9a831nnUJDVc5PulddGl7QyMjK45ZZb2l6ffvrpjBo1iiVLljBy5EiCwSC33HILDz/8MJs3b26b4XjKlCl8+OGHPRF60tLICPetlSHWbwuzLdJW4vO4kzTGu7qkprqClR+8zYolb7Di/TfZuG4lv31kEcNHT2TsxEM560vXMHHKUUyYMmO30kgi6yiBtPayyk73MDTXS0GWh5x0D5lpHjJSBb9VW/WqpE0s/cn27dtZvXo1kydPBuCmm25i4cKFvPnmm+Tl5TF37lyys7P5xS9+EedI+wdH3alSNleEWbstRHW9+wdBqp+4l0oa6mtxwiEyswfw4aLX+NG3z0BVCaSkMWHKDC6c90MyMnMAOOzIWRx25Ky4xRqNUFhpDrorT7aXlSYMzfdRmOOJNKC7CSSZVppMZJZY+rhgMMhFF13EJZdcwoQJE6itreW3v/0ty5YtY9gwd2mbc845h8cff5zRo90pL95++21uuummtjmWrrrqKg4//HCOOOIIJk+eTENDA3fffTeFhYVt2+rr67n99tuZPXt23N5rPDmOUl7r8NmOEGu3hWhscRuzM+I8VUooFGT1ikV8uOg1PnzvNVavXMSF3/wh/3PJdYyZcCgXzvshhxx+HGMnHZ7QU7W3lkKag25jemspJMUn5Gd7GZjtZUCWm0Cy0uJfEjRdS9rE0h8GSDqOw8UXX0wgEOD3v/89AG+88QajR49m3LhxbcdVVlYyaNAgwF0I7IorruCFF16guLiYYDBIa3fr1kXA3n33XW688Ubmz5/ftq20tJSzzjorqRJLazLZUBZi7dYQzZFldTNThfzM+K1D0thQR1p6JqFQkK+ePpbqynI8Hg9jJh7GWRdfy2FHuguqZmblcu5XE29dndYk0hR0G9VbJ8bMShdGFLmlkNwML9lpbrfdeLdLme5L2sSyP20siURV+drXvsb27dtZsGBB29osO3bs2G31R1Xlqaee4uqrrwZgwYIFnHXWWRQXFwPg9/s56qijdltJctq0aWzatGm3+1VVVZEMY54cdQcqbtwRYs2WSDLxQFaKkBmnZXVbmptY8cGbLHn7Rd5/599k5eTx8wdew+fzc/bF36GoeARTSo5PyDYSVbcE0tSyqyQCblvImHwvRTmtXXk9+H2WQPqLhEssIjIb+A3uQl/3q+qde+yfi7vOvQOEgGtU9a1ozu1PLr/8clatWsXLL79MWlpa2/aDDz6Y999/n6VLlzJ+/HhuvfVWRITzzjsPgC1btrQllVdffZXbbruN7Oxsfvvb37ZdY+HChUyYMKHt56OOOorly5fz9NNP9+I77D2qys46h407wqzeHKSxRfFESibxSiatHrv/Jzz50K9oaW4kkJLKwYd9jpJjdnUrP/Oiq+MY3d7CjtLY4raJtM6om5nilkQG5XrJzfSQm+GxxvR+LqESi4h4gbuBWUApsEhEnlHVj9od9grwTGRp4inA34AJUZ7bL2zcuJE//OEPpKSktFVxAfzhD3/goosu4sYbb2TOnDmoKqeeeupuJZrBgwfTurDZzJkzmTlzJiUl7uDZhQsXcsIJJ5CZmcldd90F7Koee/jhh3n11Vc56aSTevfN9qDaRoeNZSE+3hyktrFdNVcP9vzrjOM4rPtkKe+98S8Wv/08P7rrn+TmFVE8bAyzvnAJhx99CgcfdhwpqWn7vlgvaa3SamzRthl3fV4YmONlcJ6X/Cw3kaT6LYkkm4RKLMB0YK2qrgcQkceAuUBbclDVunbHZ7BrEtB9nnug8jI9UXcJ3t/rR2PEiBFdVkvdeOON3HjjjR3umzNnDieccAKXXXYZgwcPbuuKDLuSSKv21WMXXXQRM2bM4Lrrrtutqq2vaQ4qpeUhVm0OUlbttE0EGK8G+O1bNvLUX3/Nu288x84dW/F4PEyYchTVlTvIzSvi+FPP4/hTz+v1uDqiqjQFoaFl1/+9rDRh3GAfxQN85GV6yEyzcSEm8RLLEKB95X4pMGPPg0TkLOCnQBHw+e6cGzl/HjAPYPjw4VEHF83gxUSXn5/Pfffdx4UXXoiI4PF4uOaaa6I695JLLmH+/Plcf33iNQh3xXGU7VUOq7cG+XR7CEchzQ/5cUgmTY31vP/Ov8kZUMjkQ48B4LUFj3DokbOYftwcSo4+lezcgl6NqTOOozQF3WTSOjNvfpaH8UP8FEVGqqcGLImYvSXUJJQi8kXgVFX9euT1xcB0Vf12J8cfB9ysqid399xWNgll/1Xb6LB+e4hVm9x2E5/Xrerq7bEODXU1LHr7ed559WmWvPMSLc2NHD/7PL5z64MABFua8QdSejWmjqgqjS1ub61WRTlehhW4jex5WdY2YnbpS5NQlgLD2r0eCmzp7GBVfUNExohIQXfPNf1TKKxsqQzz0WdBtlaGESAzrffbTYLBlrZxIzdeMZv1n3zIgIJBnHzGlzl65plMmnp027HxSiqqbiN7Q/OuRFKQ7WXycEsk5sAkWmJZBIwTkVHAZuB84ML2B4jIWGBdpPH+MCAAVABV+zrX9F+1jQ7rtoZYWRqkJaQEvPR6u0lTYz2L336BN//9BKs+/A9//OfHpKSmcfEVt5KWnsX4g6fj8cS3l1korNQ3u+NHEMhJ93DICD+DB3jJz/YQsC6/JgYSKrGoakhErgRexO0y/KCqrhSRyyL77wPOAb4sIkGgEThP3fq8Ds+NyxsxvcJxlG1VDis3tVBaHo5Mcy5k9XIX4Y3rVvLEn3/Be2/8i6bGegYUDOJzs75Ic3MDKalpcZ02xa3ecqu4EAj4hJFFPoYXuAMR01Pim+hM/5RQbSy9qd3I+2+0dr81fUNzUNlQFmTZhiC1TUrA5yaU3iqdqCofL3+XrOwBDB05njUfLeHWq8/k6Jln8rlTvsikacfg9Xp7JZaO7FkqKcz2Mmqgl+LIOBLrtWViwVaQ7II13vcdtY0On2wO8tGmIGHH7Sbcm2Mktm5ax2vPP8rCFx5j2+ZPOfXMr3LF//7OXe8jFIzbXFyt40nqmxUF/F4YUehjRKGPolyvjSMxPaIvNd4bsxtVd76ulZ+5XYVFIDsOkxDecd25vPfmvxARphxxAud9/X858vgvAO4Swr2dVFp7cLVOhpmVJkwdGWBIvjtlvMdm+TVxZInFJCRVZWulw9JPW9hWFcbnFQZk9s7gO8dxWLb4dd57cwFfv/bneDweJk49kglTZnD8qedTMHBIj8fQEVW3iqs5CAgUZHmYOtLP4Hwf2Wk2WaNJHJZYTEJxVNlcEWbJuhZ21jmk+HpvIGPZ1s945bm/8Mpzf2HHtk1kZOVyxnlXUDx0NGdf/J0ev39HHFUa2iWT4lwvY4t9FA/wkhHnecyM6UzSJpb+MG1+f+I4yqbyEO+vD1JZ75AW6N2R8Ss/eJsbLz8VgKnTZ3LJlT9mxnGnE0hJ7ZX7t7dnyWTwAC/jin0U5/lIs5Hupg+wxntrvI8rR5UtFWHeW9tCVZ1Degq90gV2a+l6Xnr6T+QVFHPGeVcQCgV56q93cfyp51FUHP00P7GiqjS0KE2RbsHFkWQyJM9n06aYhGSN9ybhqLpjUBatbWZHjUNGAAqyezahhEMhFr39PM8/+UeWvvsKHq+XU8/8KgA+n58vXvq9Hr3/ntomdYyMfC/M8VAy1s+QPK+NLzF9miUW0+t21jksXttMaXmY1AAU9FKV190/vZJXnvsL+UVDuOAbNzHrC5eQXzS4x++7p2BIqW1SHIXcDA8zRvgZVuAjK82SiekfLLGYXtPQ7Pby+mRzCJ9XyM/quYSiqqxesYh/PXEfF877IYOGjGL22V9n+nGf54hjTsPr693/+mFHqWtyBy2m+IVDRgQYVeRlQKbHenOZfidpE4s13veesKOs3hJk0doWHAcGZEiPjbMIBlt4+5V/8Nzj97DmoyWkZ2RzzMyzGTRkFAdN7rA6uMe0TqfS0AJeD4wq8jGu2M/AXBtnYvo3a7y3xvseVVYd5q1VzVTVO+SkSY+uax4MtnDFF6dRtnUjQ0YcxOnnXs6Jcy4kLT2zx+7ZYRyRqi5Vt91o0lA/Qwt8pNgIeNOPWOO96XVNQeX99S18XBokze8O5usJWzet4703/8XcC6/C7w9w+nmXM2zkeKbNOLlXZxJWdau6mkNuVdeUkQFGD/SRk24DF03yscRiYkrVHY/y1qoWWkJKXg9Ve635aAn/+Muveee1f+L1+Tn6pLMpHDiUuRd0ua5bzLVESicAQ/O8TBzqZ9AAb69POWNMIrHEYmKmKai8t7qZNVtDZKUKmZmxLzFs37KB399xBcsWLyQjM4dzvvxdPn/uZeQVFMf8Xp1pXzpJDwglYwKMGugj00bCGwMkcWKxxvvY2rozzOsrm2gOKgUx7u2lqlRWbCOvoJis7DwqK8r4ylU/4ZQzv0p6RlbM7rMv7dtOhhZ4mRQpnfT2UsfGJDprvLfG+wMSCitLP21h2cYgGSkS0ylHHMfh3Tee4/EHfko4FOQ3D7+Hx+NBVXt17ZWGZqUxCCk+YfIwP2OKbcyJMX2q8V5EZgO/wV0F8n5VvXOP/RcB10de1gGXq+qHkX0bgFogDIQ6e9MmNmobHV5f0cSOaocBmRKzv9xVlUVvPc8jf7iNT9csp3joGM796vdp/SOoN5JK2FFqGyHkKEXZHo6a4Gdovs/aToyJQkIlFhHxAncDs4BSYJGIPKOqH7U77FPgeFWtFJHTgPnAjHb7T1TV8l4LOklt2Rnm1eWNqMZ+oOO7bzzHT79/PoOGjubqH83n+FPO67UBjc1Bt7rL64FxxT7GD/GTZ4MYjemWhEoswHRgraquBxCRx4C5QFtiUdX/tDv+v8DQXo0wyTmqrPgsyOI1LWSmSswmSFy76gMqyjYz4/jTOeKY0/jObQ9yzEln4/P5Y3L9rrSv7koPCNPHBRgzyG8zCRuznxItsQwBNrV7XcrupZE9fQ14vt1rBV4SEQX+oKrzOzpJROYB8wCGD+/9mWz7qmBIefvjZtZtCzEgIzarOJZt/Yy/3HMzb7z0d4aNmsgRn5uD1+fj+FPPi0HEXXMcpabJbScqyvZwzKQAQ/KsMd6YA5VoiaWj3+gOexeIyIm4ieXYdpuPUdUtIlIE/FtEPlbVN/a6oJtw5oPbeH/gYfd/9U0Ory5vorzGiUmvr4a6Gp748y945rHfI+Lh3K9cz1lfuqZXBjWGwkpNo/vPPmqgj0nD/BRkWXWXMbGSaImlFBjW7vVQYMueB4nIFOB+4DRVrWjdrqpbIt/LROQp3Kq1vRKL6Z6ddQ7/XtpIcwjyYzSCfvVHi3nyoV9ywmkX8KXLb6FwYM/XaDYH3fEnPi8cMiLA+CE29sSYnpBoiWURME5ERgGbgfOBC9sfICLDgX8AF6vq6nbbMwCPqtZGfj4FuK3XIu+nyqrDvPhBI14P5KYf2EN4/SdLWfPR+5x61leZNn0m9/79QwYP79lxRO5EkNDQoqSnCEeODzB6oN/m7TKmByVUYlHVkIhcCbyI2934QVVdKSKXRfbfB9wM5AP3RKouWrsVDwSeimzzAY+o6gtxeBv9xqbyIK8sayYtcGDjUxob6nhk/u089/g9DCgo5oTTLiAlNa1Hk0rr8r5NQRiQ4WHGQQGGFVh3YWN6Q9IOkGw38v4ba9asiXc4CWfd1iALP2omK1UO6K/7RW89z30/v4by7aXMPvvrXHzFrWRm5cYu0D2oKrWNSksYBuZ4mDY6QPEALx5rPzEmproaIJm0iaWVjbzf25otQd5Y1UzuAU5zv2PbJr559sEMHXkQl9/wOyZOOTKGUe7OUaWm0W2YH1bgZerIAIXZ1iBvTE/pUyPvTXx9sjnIW6uayU3f/6SyfvWHjD5oKoWDhnHr755lwpQj8fsDMY7U5ahS0+COkB9V5GPKSD/5Wd4euZcxJjrWJca0aUsqGfuXVBob6vj9T77FtRcfzbLFCwE45PDjeiSpOKpU1SuVdcrIIi9nzUjnxENSLakYkwCsxGIAt03lrVXN+z3wcc1HS/jlD7/Cts3rOefL32Xi1KN6IEq3DaW6wZ3La/QgH1NHBsjNsL+PjEkkSZtYbNr8XTaUuQ31uen7l1QWPDGf+3/1PQYUDOLH9zzPwYd9LuYxqrqj5IMhZdRAH9NGBhjQA+u9GGMOXNImFlV9Fni2pKTkG/GOJZ42V4R4bXkz2QfQUJ+SmkbJMbO56of3kZk9IKbxtS2qFYRhBV4OGxOw6i5jEpz1CkviXmE7asIsWNJIWkBI7WaX4g1rV7Bt86ccefwZAD2yRkpDszsWZWCuhyPGplCUY728jEkUMe0VJu5v9hjcQYoVwDpN9uzUB1U3OLz0QRMpPrqdVN557WnuuvUb5OYVcfjRp+L3B2L6wG8JKTUNSna6MGtqCsMKfJZQjOlDupVYRGQm8GvcNVGqgAHASBG5VlVfjX14pic0NDu89EEjCKSnRN9Ooao8+sc7ePyBn3LQ5CO44WePxrTHV9hRqhsUvxeOnZjC6EE2Ut6Yvqi7JZbbcRfS2tm6QUTygGeBPpVYkrXxPhhSXl7WRGMQctOjf2iHw2HuuvXrvPHi35j5+S9x+fW/IZCSGpOYtHVwo6McMtzPwSMC3S5FGWMSR3cTiwA1e2yrpePp7hNaMjbeO47y5qpmKmqcbs9S7PV6ycsfxMWX38I5l1wXs6qpxha3cX54oZfp41LIOcCJLo0x8dfdxHIP8K6IvAFUAnm466HcFeO4TIypKu+vb+HT7SEKsqJPCju2baKhvoYRYyZz6VU/iVlCCTvuAMf0VOGUaSkMzbd2FGP6i24lFlX9q4g8i7uqY2vj/Y9UtbongjOxs2ZriKUbguRnRr9I1+bP1nDzlaeTkprO7x5djNd74N18VZXayHiUqaMCHDLcf0DzkRljEk+3e4VFkshLPRBLr0qmNpay6jBvr2pmQLpEvezu+tUfcstVcwHlxv97PCZJJRhSqhqUQQM8HD0+zQY4GtNPdes3W0RGishDIvKAiExot/03sQ+tZ6nqs6o6LycnJ96h9Ki6JoeXP2wiLRD9AMjVKxdz0+Wn4Q+k8JM//JvR46cdUAwamderrlk5dmIKpx1qScWY/qy7v933A48CjwB/ifzVDzAlVgGJyGwR+URE1orIDR3sv0hElkW+/iMiU6M9N9kEQ8qry5oIO5CeEn110z8e+hVZOXncOf9lho446IBjKK9VinI9nH1kOuOH+PFEWWoyxvRN3a0K86nq8wAicgLwNxEpBGIyQFJEvMDdwCygFFgkIs+o6kftDvsUOF5VK0XkNGA+MCPKc5OGqvLOJ82U1zoUdLMH2LW33k9dTRX5RYMP6P61TW4j/bETUxg32GeLbRmTJLpbYvFF1pNHVeuBM4E5QIfD+vfDdGCtqq5X1RbgMWBu+wNU9T+qWhl5+V9gaLTnJpOPS4Os2RoiPzP6NpVbrzmT+rpqUlLTDyiphB2lok7JSRfOnBEppVhSMSZpdLfEchWQCdQDqGpQRM4FvhijeIYAm9q9LsXtgdaZrwHPd/dcEZkHzAMYPnz4/saasMqqw7yzuoUBGdH1ANv82RpuueoL+AOpNNTVkJG5/+1OreNSDhsdYMpIf9SdBYwx/Ud3uxu/38E2B3g8RvF09BTqsJpNRE7ETSzHdvdcVZ2PW4VGSUlJv5rnrKHZ4ZVlTaQHopsCv3z7Zn70bbep7NbfPUvhoGH7dV9Vt8eX3yvMOTyN4gE2A7ExyepAeoVNbLc9Vr3CSoH2T7ahwJYO4piC25FgrqpWdOfc/izsKAtXNhMMaVSN9TXVFfzoqi9QV1PFj+7653431IcdpaJWGZTr5cwZllSMSXYH0ivsoR7oFbYIGCcio0QkAJwPPNP+ABEZDvwDuFhVV3fn3P7uww1BtuwMkxPlHGB11ZWActMv/86YCYfu1z2bg8rOOmXa6ACzpqZ2a1JLY0z/lFC9wlQ1JCJXAi8CXuBBVV0pIpdF9t8H3Iw76v+eSPtBSFVLOjs3FnH1BZsrQnzwaQt5UYysb13lYPDwsfz2kUX7PfixttEh7MCsqSkML/Tv1zWMMf1Ptxb6EpG3gFMjPcIQET9uCeYUVc3umRB7RruR999Ys2ZNvMM5IPVNDv98txG/F1ID+y6tPHT3zdRUV3D59b/dr6TS2p6SniLMmppma84bk4S6Wuiru0+E1l5hgNsrDDgX6HMzBPeXkfduu4o7CDKapPLyM3/myYd+icfjxePpfkJw1G1PKcrxcnpJuiUVY8xeEq1XmOmmFRuDbKt0yI9ixuJVy/7LvT+7mqnTT+Sb1/2q27MJhx23PWX8EB9HHpRii3AZYzrU3RUkhwCHA4e1fqnq0K7PMj2lrDrMkvUt5EYxXqV8+2buvP4CCgYN43s/fgivr3vNa8GwO9/XoaMDTBtlAx6NMZ3r8ukiImezeyIpiOxS4J/AL3syONO5pqDy2vIm0gNEVXLYsmkNHq+XG//vcbJy8rp1r+agUtvoTs1y0GBbN8UY07UuG+9FxMFdc+UxYAWwLPJ9NTBVVct6I8ie0Jcb71Xd8SobykLkdWOW4Jbmpm4vJ9wUVOqblZMOsZ5fxphdDqTx/k4gCxgO/FtV31HV2lgHGA99ufH+07IQ67aFGJARRWP9sw/x3N/uRVW7nVQaW5SGZuWUqamWVIwxUesysajqD4CpQDqwQkR+LCLpvRKZ6VBdk8Pbq5rJTtt3u8q6jz/gvp9fw6I3F+A4Trfu09CsNAdh9qFpDMnv9npwxpgkts96FFX9RFVn4c7L9RXgE9xSjOlljipvr2pGVUjxd51U6moq+dn/fons3AK+c9uD3Rqv0tiitIRg9mGpDLLpWYwx3RR1Bb2qPgqMB54EUoAHRKSg67NMLK3eEqS0IkzOPsqMqsrv7riC8u2lfP8nfyFnQGHU92hqURpblFMOTaUox5KKMab7ukwsInJE+9eqWqeq1+D2FMvDbcjvk0TkDBGZX11dHe9QolLT4PDuJy3kpu+7Cmz1ysW8u/BZLr7iViYc0tWqA7trDir1Lcqp09IYlGtJxRizf/bVK2wNkAo8jdu9+DVVDbfb/xVV/VNPB9mTSkpKdPHixfEOo0uOKi990ERZjUNulBNMrvloCWMmHBr16PpgSKluVGZNTWVYgbWpGGO6tt+9wlR1HHAasA34GVAmIn8RkXNEJL2vJ5W+Yv22EJt3hslJ6/q45qYGVi37LwDjJh0edVIJhd25v46flGJJxRhzwKJpvF+hqj9W1cOBacB7wLeAbSLyjIh8VUS6N+LORK2h2eGd1c3kRFEF9qff/oAfXHYK27dsiPr6jqNU1iszxgUYU2xdio0xB65bMwiq6iZV/Z2qzgRG4jbknwF8tQdiS3qqyrtrWnAcCPi6TipL/vMizz/5R8447woGDh4Z9fUr6pRDRviZPNySijEmNrpd7yHun81jcNdEeRt4SLsz936CaDfyPt6hdGpLZZj120IU7GOCyZqqcn7348sZMWYSX7rslqivX1mvjCryUTImYNO0GGNiprtLE88ElgK/AC7HnStsaWR7n5LoI+9DYeWdj5vJSOm6CkxVuefOq6mt3sm1tzwQ9ej66gZlQKaHYyel4PFYUjHGxE53Syy3Ayeq6s7WDZH2lWeBV2MRkIjMBn6Duwrk/ap65x77JwB/wp0U80ZV/UW7fRuAWiBMZGXJWMQUD59sDlLTqBRk7Tv3TzhkOhOnHsmog6JbIbqh2cHnEU6akrrPKjZjjOmu7iYWAWr22FYb2X7ARMQL3A3MAkqBRSLyjKp+1O6wnbgLjp3ZyWVOVNXyWMQTL/VNDkvWtZCTtu+PVUQ486Kro752MKQ0tsDpJalkptoiXcaY2Ovuk+Ue4F0R+bWI3CwidwHvRLbHwnRgraquV9UW3FmV57Y/QFXLVHUREIzRPRPO++tbcBT8XZQmVJXf3DaP/7z6z6iv60SWFD56QgqFNqreGNNDutsr7K/ATOB5YA2wADgpsj0WhgCb2r0ujWyLOkTgJRFZIiLzOjtIROaJyGIRWbxjx479DLVnlNeEWbM1RO4+Zi5+7flHePVfD1O+vTTqa1fWu6s/HjTYxqoYY3pOt58wqloNvNQDsUDHVWrd6XF2jKpuEZEi4N8i8rGqvrHXBVXnA/PBHXm/f6HGnqqyZF0LAS9drtBYsWMr9//q+0yaejSnn3dFVNeublTyMj3MOCjFeoAZY3pUd3uFnSwi74nIfBGZIyKfisiHInJijOIpBYa1ez0U2BLtyaq6JfK9DHgKt2qtz9he7bB5Z5isLtpWVJX7fnY1wZYmvn3TvVGNrm8OKqhy4sGp+G2demNMD+tuieUO4AIgG3gBOAKow32IvxaDeBYB40RkFLAZOB+4MJoTRSQD8KhqbeTnU4DbYhBTr1BVlqxtJsVHlyWKj5b+h/fe/BeXfvsOBg/f9xgcR5WaRuWkKSlkp1tjvTGm53U3sTSq6joAEVmtqp9Ffg7FIhhVDYnIlcCLuN2NH1TVlSJyWWT/fSIyCFiMm9wcEbkGmAQUAE9FHso+4BFVfSEWcfWGLZVhtlU5+xwMOWna0fzwV09y6IyTo7puZZ0yYaiPEYXWrmKM6R3dfdoEREQiI+0vPoDrdEpVF+B2Cmi/7b52P2/DrSLbUw3uapdRSaSR944qi9a0kB7ourRSU11Bdk4+JcfMjuq6dU0O2ekepo+1dhVjTO/pbt3IGUQa2FV1A4CIBIDrYxtWz0ukkfdbKsLsrHNIT+m6CuzrX5jAssWvR3XNUFhpDsGJB6d02W3ZGGNirVslDVWt6GBbC/DfmEWUZFSVDz5tIdXfeWklFApyz51XkTOggIMmH9HhMXtes6pBmT42QF6WjVcxxvSuROsVlnR21DjsqHHI6KK08tzj97Dp01XMu+5XpKZl7POaNY0wMNfDpGE2Y7ExpvclWq+wpLPisyB+b+ellYqyLTx6/0844nNzOOLY0/Z5vWBIcVT53MQ0m1zSGBMXCdUrrDclQuN9dYPDhrIQ+ZmdJ4Cl772COg5fv/bn+7xeaxXYsROta7ExJn66+/QJyK4/rXukV1hvSYTG+1Wbgng9XfcEO+n0i5n/1EoGDRm1z+tVN8KQPC/jbMoWY0wc7U+vMKDv9wqLt4Zmh483B8nuZJR9OBRi3ccfAJCbV7TP6wXDiuO4E0x2NR2MMcb0tO5OQlmx52qRqtqiqtYrrJvWbwuhCt5O2kFeevpPfOeSY1nz0ZJ9XktVqa5XjhgbsCowY0zcdavORESGAIfjLrJ1GHCYqnY0WNF0wVFlZWmQzNSOk0p9XTWP/vHHTD70WMZOPGyf16trUgqyPUwYar3AjDHx12ViEZGz2T2RFER2KfBP3KWJ+6R4Nt5vr3JoaFLyO1kd8ok//4LqynJu/vWd+xwxH3aUlhAcOzGl09KPMcb0pn3VmzwBzAPWAjcDxwK5wA7gClX9dY9G14Pi2Xj/yeYgvk7GLW7fspFnH7ubE067gLETD93ntaoalCkjbSCkMSZx7Cux3AlkAcOBf6vqO6pa2/Nh9V+NLcqGshBZnVSDbVi7gvSMbL50+S1RXSs9IEwZYVVgxpjE0WViUdUf4E7smA6sEJEfi0h6r0TWT20oC6LQ6eDFGcd9ngee+YTCgV03XakqdU3K0RMCNheYMSah7LMLkap+oqqzgK8BXwE+wS3F9GkicoaIzK+uru61e6oqKz8LkhHoOBGsWvZfVBV/IGWf16pqgJFFPobm25gVY0xiibpvqqo+CowHngRSgAdEpKDrsxJXPNpYdtQ41DYqqR0klk9WLOKGb5zEi089uM/rhMLuipDTxwVsOnxjTMLp7jiWOlW9BrenWB6wItYBichsEflERNaKyA0d7J8gIu+ISLOIXNedc+Nt/bZQp1Vgj8y/nezcAo479dx9Xqe6QTl0TICsNBuzYoxJPPv1ZFLVZap6DPC/sQxGRLzA3cBpuKtCXiAik/Y4bCdwFfCL/Tg3bsKOsnZbkKzUvfet/OBtlr77CmdffC3pGV3XMja2KBmpwkQbs2KMSVAH9Cevqv4pVoFETAfWqur6yDovjwFz97hnmaouAoLdPTeeyqocgmHwefcusTwy/3YG5A9kzv/M6/IarQ32R40P4O/gOsYYkwgSrS5lCLCp3evSyLaYnisi80RksYgs3rFjx34F2l3rtwfpKBdU7Sxj+5YNnPPl75KS2nWHu+pGGFbgtQZ7Y0xCS7QnVEd/hmsH2w7oXFWdD8wHKCkpifb6+y0UVtZvD3U4hUtuXhH3PbkcdZwurxF2lLCjTB9n69cbYxJbopVYSoFh7V4PBbb0wrk9ant1mJCzdzVYTXUFwZZmfD7/PrsYV9UrBw/3k5uRaP9kxhizu0R7Si0CxonIqMh0/OcDz/TCuT1q/bYQvg4+6T///odcef7hhENdr5PWElL8PjhkRKCHIjTGmNhJqKowVQ2JyJXAi4AXeFBVV4rIZZH994nIIGAx7vLIjohcA0xS1ZqOzo3LG2knGFY+7WAKl53lW3n9+UeZ9YVL8fq6/meoaXRH2Kf6rQrMGJP4EiqxAKjqAmDBHtvua/fzNtxqrqjO7UxvzW68vSpM2Nl73ZXnHr8XJxxi7oXf7vL8hmYlO00YO8i6Fxtj+oZEqwrrNb018n5j2d7VYA31tbzwj/s56sQzKR46uqsYqW9Wjhyf0mE3ZWOMSURJm1h6g6PKhh1hMlJ2TwrvvfEc9XXVzL2g69JKbaMyOM/LkDybEt8Y03ckXFVYf1JZ5xAMKlmpu+fv42efT/GwsRw0uaTTc1WVljAcMdbmAzPG9C2WWHrQ1p1htIOcICKMP/iILs+tboDRA30UZFtpxRjTtyRtVVhvTJv/aVmItD16ct3/q+/xl3t+1OV5rYMhDx1t3YuNMX1P0iaWnm68b2xRymsd0trlhtrqnbz4zwepq63q8tzqBmXiUD856Un7z2OM6cPsydVDyqrDoOzWPvLagkdoaW5i9llf6/S8UFgRgSkjrXuxMaZvssTSQzaWhfC1ax5RVV546gHGHzydUQdN6fS86gblkBEB0lPsn8YY0zfZ06sHOI7yWfnu3YxXffgOmzeu5tSzvtrpeaGw4vXCJFtrxRjThyVtYunJxvuddQ7BsO42qDF7QAGnnTOPo2ee1el51Y3KlBGBDpcuNsaYviJpE0tPNt5vrQzvtW3oiIO47Pu/Ji09s8NzQmHF64EJQ6y0Yozp25I2sfSkLZVhUtuNEFq/+kM+Wf4eqp0v/VLdoEy10ooxph+wxBJjqsqO6vBuCeLvf/o/7vjeF3HCe5dkwJ0B2euF8da2YozpByyxxFhdkxIM75rNuK62ikVvLeBzs77Y6fT4Na2lFZsW3xjTDyRtYumpxvvq+t2XGH7ntacJtjRz/OzzOzw+7CgeD4wbbLPrGGP6h6RNLLFovF+2oYXKut0TSXmtQ/tyx+vPP8bgYWMZN+nwDq9R3aBMGOK3cSvGmH4j4Z5mIjJbRD4RkbUickMH+0VEfhvZv0xEDmu3b4OILBeRpSKyuKdjLat2+Kx892WFt1aGSY00lTTU1bBh7XKOn31+hzMUO46iCpOGWduKMab/SKj6FxHxAncDs4BSYJGIPKOqH7U77DRgXORrBnBv5HurE1W1vDfidVRZuy3IlBF+RARHlfKaMJmRZYjTM7P503NrCQVbOjy/phHGDPKRlZZw+d0YY/Zboj3RpgNrVXW9qrYAjwFz9zhmLvCQuv4L5IpIcW8H2qq82qGuye1GXNeoey1DHEhJJT0ze6/zVJWQo0webqUVY0z/kmiJZQiwqd3r0si2aI9R4CURWSIi8zq7iYjME5HFIrJ4x44dBxRwS9hd1x6gqmFXe0tlxTauPL+EDxe91uF5dU3u6pD5WbbeijGmf0m0xNJRf9s9RxV2dcwxqnoYbnXZt0TkuI5uoqrzVbVEVUsKCwv3P1ogzQ/rt7vtLBU1uxLLuwv/xaZPV5EzoOPrN4fgkBFWWjHG9D+JllhKgWHtXg8FtkR7jKq2fi8DnsKtWutR6SnClp1hWkIaabh38947rz9N8dAxjBgzea9zmoNKZqpQnGulFWNM/5NoiWURME5ERolIADgfeGaPY54BvhzpHXYkUK2qW0UkQ0SyAEQkAzgFWNHTAQtucWl7VZjy2jApfqirqWT54oUcNXNuh73B6pqVQ0b48XhsQKQxpv9JqF5hqhoSkSuBFwEv8KCqrhSRyyL77wMWAHOAtUAD8JXI6QOBpyIPch/wiKq+0BtxewU+2hTEiTTcv/fWAsLhEEedsGe/A3dApACjiqwazBjTPyVUYgFQ1QW4yaP9tvva/azAtzo4bz0wNdr7iMgZwBljx47d/2AjMlKFbVW75gEbWDyCWXMvZezEw/Y6tqYRDhrss8kmjTH9VqJVhfWaWE6b7/cKwZC29SqYfOixXPmDu/F4dv94VZWwo4y3qfGNMf1YwpVY+qqcdA8+L2z+bA1er49BQ0btdUxDs1KY7SEvM2nzuTEmCdgTLkZS/ILXI/ztgTu57ivHE+5givymIBw83N9hg74xxvQXllhiKBwKsfg/L1JyzKl4vbt3JQ477gqRQ/KtkGiM6d+SNrH0xLT5q5a9Q11NJdM/9/m99tU0wthiHwGflVaMMf1b0iaWnljz/r03/4XPH+DQI0/ea1/YUcYWW6O9Mab/S9rE0hMWv/0iU0tOIC09c7ftTUElK00oyLaP2xjT/1mFfwzd+ceXqa3eudf2+mbliLEBPNZob4xJApZYYig7J5/snPzdtrnjOWFUkX3UxpjkkLR1M7FuvP/jL6/j7Vf+sdf21unxM1KT9qM2xiSZpH3axbLxvqpiO//6+32Ubli9177mEEy0kfbGmCRi9TMH4D+vPctbr7/I8kWvoKpM/9yc3faHworfC8V5Nj2+MSZ5WGI5AE8/di+frFjCIYcfx0WX3cyog6bstr+2CcYV+/B7rdHeGJM8kjaxxGJ24+//+AFC3gHkZKZ0uN+xsSvGmCRkbSwH0MaSX1iMz9dx4mgOKukpNnbFGJN8Eu6pJyKzReQTEVkrIjd0sF9E5LeR/ctE5LBoz+1N9c3KhKF+G7tijEk6CZVYRMQL3A2cBkwCLhCRSXscdhowLvI1D7i3G+f2ClVFFUba2BVjTBJKqMQCTAfWqup6VW0BHgP2XN93LvCQuv4L5IpIcZTn9orGFiU/20N2mpVWjDHJJ9ESyxBgU7vXpZFt0RwTzbkAiMg8EVksIot37NhxwEHvqaEFJg61dVeMMckp0RJLR09ijfKYaM51N6rOV9USVS0pLCzsZohdq29yyM3w2BQuxpiklWhPv1JgWLvXQ4EtUR4TiOLcHuU4SkMQTpqags/GrhhjklSilVgWAeNEZJSIBIDzgWf2OOYZ4MuR3mFHAtWqujXKc3tUVYMyYYiPohwbaW+MSV4JVWJR1ZCIXAm8CHiBB1V1pYhcFtl/H7AAmAOsBRqAr3R1bm/F3hx0p285fHTHgyWNMSZZJFRiAVDVBbjJo/22+9r9rMC3oj23M7EYed9eXZNy7KQUUgNWBWaMSW6JVhXWa2I5u3HrmivD8hMuTxtjTK9L2sQSS81ByM3wWGnFGGOwxBITDS3KCOtebIwxgCWWmFCFQbnWE8wYYyCJE0usliZWAIH8rKT9KI0xZjdJ+zSMVeN9U9BNKil+a18xxhhI4sQSK40tyohCa18xxphWllgOlMJAa18xxpg2llgOUFqKWPuKMca0k7RPxFg13g/K9eK3CSeNMaZN0iaWWDTeezzCiEKrBjPGmPaSNrHEQmGWh8F51nBvjDHt2VPxAEwdFYh3CMYYk3CsxGKMMSamkjaxxKrx3hhjzO6SNrHEctp8Y4wxuyRMYhGRPBH5t4isiXwf0Mlxs0XkExFZKyI3tNt+i4hsFpGlka85vRe9McaYVgmTWIAbgFdUdRzwSuT1bkTEC9wNnAZMAi4QkUntDvm1qk6LfEW1kqQxxpjYSqTEMhf4c+TnPwNndnDMdGCtqq5X1Rbgsch5xhhjEkQiJZaBqroVIPK9qINjhgCb2r0ujWxrdaWILBORBzurSgMQkXkislhEFu/YsSMWsRtjjIno1cQiIi+LyIoOvqItdXQ0d4pGvt8LjAGmAVuBX3Z2EVWdr6olqlpSWFjYnbdgjDFmH3p1gKSqntzZPhHZLiLFqrpVRIqBsg4OKwWGtXs9FNgSufb2dtf6I/BcNDEtWbKkXEQ2drI7B+iqP3IBUL6PW+zrGtEcE4trwL7jjdV9YnGMfbb7d0wsYu3NWOyz7blYevqzHdHpWaqaEF/A/wE3RH6+Afh5B8f4gPXAKCAAfAhMjuwrbnfctcBjMYhp/j72Lz7Qa0R5nwO+RjTxxvA+sXjP9tn23Puxz9Y+25h9th19JdKULncCfxORrwGfAV8EEJHBwP2qOkdVQyJyJfAi4AUeVNWVkfN/LiLTcKvGNgDfjEFMz/bSNfZ1TCyuEY1Y3acvxduXYo3mmFjEGu117LPdP/3ts92LRLKS2Q8islhVS+IdR7T6Urx9KVboW/H2pVihb8Xbl2KFnos3kXqF9UXz4x1AN/WlePtSrNC34u1LsULfircvxQo9FK+VWIwxxsSUlViMMcbElCUWY4wxMWWJZT9ERvaXiciKeMeyLyIyTEReE5FVIrJSRK6Od0xdEZFUEXlPRD6MxHtrvGPaFxHxisgHIhLV2Kl4EpENIrI8MlHr4njH0xURyRWRJ0Tk48j/36PiHVNnRGR8uwlwl4pIjYhcE++4OiMi10Z+v1aIyKMikhrT61sbS/eJyHFAHfCQqh4c73i6EhlsWqyq74tIFrAEOFNVP4pzaB0SEQEyVLVORPzAW8DVqvrfOIfWKRH5DlACZKvq6fGOpysisgEoUdV9DeKLOxH5M/Cmqt4vIgEgXVWr4hzWPkUmy90MzFDVzgZfx42IDMH9vZqkqo0i8jdggar+v1jdw0os+0FV3wB2xjuOaKjqVlV9P/JzLbCK3edXSyjqqou89Ee+EvavHxEZCnweuD/esfQnIpINHAc8AKCqLX0hqUScBKxLxKTSjg9IExEfkE5kBpNYscSSRERkJHAo8G6cQ+lSpGppKe60Pv9W1USO9y7g+4AT5ziipcBLIrJERObFO5gujAZ2AH+KVDPeLyIZ8Q4qSucDj8Y7iM6o6mbgF7gD0bcC1ar6UizvYYklSYhIJvAkcI2q1sQ7nq6oalhVp+HOBTddRBKyulFETgfKVHVJvGPphmNU9TDcNY2+FanWTUQ+4DDgXlU9FKingzWaEk2kyu4LwN/jHUtnIjO/z8WdGmswkCEiX4rlPSyxJIFIW8WTwMOq+o94xxOtSNXH68Ds+EbSqWOAL0TaLR4DZorIX+MbUtdUtXXS1jLgKdw1jhJRKVDarrT6BG6iSXSnAe9ru0lxE9DJwKequkNVg8A/gKNjeQNLLP1cpDH8AWCVqv4q3vHsi4gUikhu5Oc03F+Cj+MaVCdU9X9VdaiqjsSt/nhVVWP6l18siUhGpAMHkWqlU4CE7NmoqtuATSIyPrLpJCAhO5zs4QISuBos4jPgSBFJjzwfTsJte40ZSyz7QUQeBd4BxotIaWTizER1DHAx7l/TrV0h58Q7qC4UA6+JyDJgEW4bS8J34+0jBgJviciHwHvAv1T1hTjH1JVvAw9H/i9MA34S33C6JiLpwCzcEkDCipQCnwDeB5bj5oGYTu1i3Y2NMcbElJVYjDHGxJQlFmOMMTFlicUYY0xMWWIxxhgTU5ZYjDHGxJQlFmOSkIg8LyKXxDsO0z9Zd2Nj+jkRuQUYm8iDN03/YiUWY/qwyOy0xiQUSyzG7ENkcazvicgyEakXkQdEZGCkOqlWRF6OTOzXevzfRWSbiFSLyBsiMrndvjki8lHkvM0icl1ke4GIPCciVSKyU0TeFJEOfz9FREXkWyKyBlgT2fYbEdkUWWBqiYh8LrJ9NvAD4DwRqYuMukdEXheRr0d+9ojITSKyUdwF7B4SkZwe+jhNErDEYkx0zsGdruMg4AzgedwHdgHu79FV7Y59HhgHFOFOm/Fwu30PAN9U1SzgYODVyPbv4k68WIg79coP6HodmjOBGcCkyOtFuNOe5AGPAH8XkdTIlC0/AR5X1UxVndrBtS6NfJ2IO119JvD7Lu5tTJcssRgTnd+p6vbIWhZvAu+q6geq2ow7S/ChrQeq6oOqWhvZdwswtV0JIAhMEpFsVa1sXYQtsr0YGKGqQVV9U7tuAP2pqu5U1cbIPf+qqhWqGlLVXwIpwPguzm/vIuBXqro+ssja/wLnWzWb2V+WWIyJTvtp0Bs7eJ0JbYuU3Ski60SkBtgQOaYg8v0cYA6wUUQWyq513P8PWIu7CNd6EdnX2iOb2r8Qke+Kuy58tYhUATnt7rkvg4H2qx1uxF0PZWCU5xuzG0ssxsTWhbiLKJ2M+3AfGdkuAKq6SFXn4laT/RP4W2R7rap+V1VH41a1fUdETuriPm2lmUh7yvXAucAAVc0Fqlvvyb6Xdt4CjGj3ejgQYvfkaUzULLEYE1tZQDNQgbuWeNtU7yISEJGLRCQnssBSDRCO7DtdRMZG1sdo3R7uxj1DuEv5+kTkZiC73f7twMjOOgPgrh9yrYiMiqw02tomE4ry/sbsxhKLMbH1EG5V0mbchan+u8f+i4ENkWqyy4DWsSXjgJeBOty1fu5R1dejvOeLuB0GVkfu3cTuVWWty+RWiMj77O1B4C/AG8CnkfO/HeW9jdmLDZA0xhgTU1ZiMcYYE1OWWIwxxsSUJRZjjDExZYnFGGNMTFliMcYYE1OWWIwxxsSUJRZjjDExZYnFGGNMTP1/HefQzcnGBocAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1,1,figsize=(6,4))\n",
    "axs.plot(mass_ratio,amp,color='k',linestyle='--',label=r'$\\mu_{\\rm GPR}$')\n",
    "axs.fill_between(mass_ratio,amp-2*sigma_amp,amp+2*sigma_amp,alpha=0.5,\\\n",
    "                 color='cornflowerblue',label=r'$2\\sigma_{\\rm GPR}$')\n",
    "axs.set_xlabel('mass ratio',fontsize=12)\n",
    "axs.set_ylabel(r'$A_{330}/A_{220}$',fontsize=13)\n",
    "axs.minorticks_on()\n",
    "axs.legend(fontsize=12,frameon=False)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}