RV_Fitting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notebook used to fit a curve to radial velocity data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Inital Imports\n",
    "\"\"\"\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "from matplotlib.ticker import (MultipleLocator, FormatStrFormatter, AutoMinorLocator)\n",
    "import matplotlib.ticker as ticker\n",
    "import mpld3\n",
    "# mpld3.enable_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Read in the data and convert to numpy\n",
    "\"\"\"\n",
    "\n",
    "data = pd.read_csv('wasp2b_radialvelocity.csv')\n",
    "time_wasp = data['Obs_Date'].to_numpy()\n",
    "time_wasp = time_wasp - 2450000\n",
    "vel_wasp = data['Velocity'].to_numpy()\n",
    "vel_err_wasp = data['Error'].to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Define Constants\n",
    "\"\"\"\n",
    "\n",
    "mass_sun = 1.989e30\n",
    "radius_sun = 6.960e8\n",
    "G = 6.67e-11\n",
    "au = 1.496e11\n",
    "alpha = np.power(G/(4*np.power(np.pi, 2)), (1/3))\n",
    "\n",
    "wasp_mass = 0.79 * mass_sun\n",
    "wasp_mass_err_pos = 0.15 * mass_sun\n",
    "wasp_mass_err_neg = 0.04 * mass_sun\n",
    "\n",
    "incl = 84.4\n",
    "incl_err = 0.4\n",
    "\n",
    "transit_time = 2453991.5146 - 2450000\n",
    "\n",
    "guess_freq = 3\n",
    "guess_amplitude = (np.max(vel_wasp) - np.mean(vel_wasp)) \n",
    "guess_phase = 0\n",
    "guess_offset = np.mean(vel_wasp)\n",
    "guess = [guess_freq, guess_amplitude, guess_phase, guess_offset]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_sin(x, freq, amplitude, phase, offset):\n",
    "    \n",
    "    return np.sin((x * freq) + phase) * amplitude + offset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func_fitter(time, rad_vel, err, guess):\n",
    "\n",
    "    # Initial guesses\n",
    "    p0=[*guess]\n",
    "\n",
    "    # now do the fit\n",
    "    fit, pcov = curve_fit(my_sin, time, rad_vel, p0=p0, sigma=err, bounds=([0,-np.inf,0,-np.inf],[np.inf,np.inf,360,np.inf]))\n",
    "\n",
    "    return fit, pcov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chisquared(fit, time):\n",
    "    \"\"\"\n",
    "    Function to return chi squared and reduced chi squared \n",
    "    \"\"\"\n",
    "    phase = fit[0]\n",
    "    chisq = np.sum(np.power((vel_wasp - my_sin(time, *fit)), 2)/np.power(vel_err_wasp, 2))\n",
    "    rchisq = chisq/(len(vel_wasp)-4)\n",
    "\n",
    "    \n",
    "    return chisq, rchisq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {},
   "outputs": [],
   "source": [
    "def radvelparams(fit, pcov):\n",
    "    \"\"\"\n",
    "    Calculate and return Planetary Properties and their errors\n",
    "    # params - [0]=freq, [1]=period_sec, [2]=star_vel, \n",
    "               [3]=semi_m_a, [4]=min_planet_mass, [5]=true_planet_mass\n",
    "    \n",
    "\n",
    "    # errors - [0]=freq_err, [1]=period_sec_err, [2]=star_vel_err, \n",
    "               [3]=semi_m_a_err_pos, [4]=semi_m_a_err_neg, [5]=min_planet_mass_err, \n",
    "               [6]=true_planet_mass_err\n",
    "    \"\"\"\n",
    "    \n",
    "    perr = np.sqrt(np.diag(pcov))\n",
    "    \n",
    "    freq = fit[0] / (2 * np.pi)\n",
    "    freq_err = perr[0]\n",
    "    freq_err1 = (perr[0] / fit[0]) * freq\n",
    "    \n",
    "    star_vel = fit[1] * 1e3\n",
    "    star_vel_err = perr[1] * 1e3\n",
    "    \n",
    "    period_day = 1 / freq\n",
    "    period_sec = period_day * 24 * 3600\n",
    "    period_day_err = period_day * (perr[0] / fit[0])\n",
    "    period_sec_err = period_day_err * 24 * 3600\n",
    "\n",
    "    semi_m_a = np.power((G * np.power(period_sec, 2) * wasp_mass) / (4 * np.power(np.pi, 2)), 1/3)  \n",
    "    semi_m_a_err1 = np.power(period_sec_err, 2) * np.power((2/3) * alpha * np.power(wasp_mass, (1/3)) * np.power(period_sec, (-1/3)), 2)\n",
    "    semi_m_a_err2 = np.power(wasp_mass_err_pos, 2) * np.power((1/3) * alpha * np.power(wasp_mass, (-2/3)) * np.power(period_sec, (2/3)), 2)\n",
    "    semi_m_a_err3 = np.power(wasp_mass_err_neg, 2) * np.power((1/3) * alpha * np.power(wasp_mass, (-2/3)) * np.power(period_sec, (2/3)), 2)\n",
    "    semi_m_a_err_pos = np.sqrt(semi_m_a_err1 + semi_m_a_err2)\n",
    "    semi_m_a_err_neg = np.sqrt(semi_m_a_err1 + semi_m_a_err3)\n",
    "\n",
    "    min_planet_mass = (wasp_mass * star_vel * period_sec) / (2 * np.pi * semi_m_a)\n",
    "    min_planet_mass_err_1 = np.power((wasp_mass_err_pos), 2) * np.power((star_vel * period_sec) / (2 * np.pi * semi_m_a), 2)\n",
    "    min_planet_mass_err_2 = np.power((perr[1]), 2) * np.power((wasp_mass * period_sec) / (2 * np.pi * semi_m_a), 2)\n",
    "    min_planet_mass_err_3 = np.power(period_sec_err, 2) * np.power((wasp_mass * star_vel) / (2 * np.pi * semi_m_a), 2)\n",
    "    min_planet_mass_err_4 = np.power(semi_m_a_err_pos, 2) * np.power((wasp_mass * star_vel * period_sec) / (2 * np.pi * np.power(semi_m_a, 2)), 2)\n",
    "    min_planet_mass_err = np.sqrt(min_planet_mass_err_1 + min_planet_mass_err_2 + min_planet_mass_err_3 + min_planet_mass_err_4)            \n",
    "    \n",
    "    true_planet_mass = min_planet_mass / np.sin(incl * (np.pi/180))\n",
    "    true_planet_mass_err = np.sqrt(np.power(incl_err, 2) * np.power((min_planet_mass * (1/np.tan(incl*(np.pi/180)))*(1/np.sin(incl*(np.pi/180)))), 2))\n",
    "    \n",
    "    \n",
    "    params = np.array([freq, period_sec, star_vel, semi_m_a, min_planet_mass, true_planet_mass])\n",
    "    errors = np.array([freq_err, period_sec_err, star_vel, semi_m_a_err_pos, semi_m_a_err_neg, min_planet_mass_err, true_planet_mass_err])\n",
    "    \n",
    "    \n",
    "    return params, errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Plotting the Radial Velocity Curve Fit and the Photometric Curve Fit\n",
    "\"\"\"\n",
    "\n",
    "def plotting(xdata1, xdata2, ydata1, ydata2):\n",
    "    plt.rcParams[\"font.family\"] = \"serif\"\n",
    "    fig, (ax1, ax2) = plt.subplots(2, sharex=False)\n",
    "\n",
    "    ax1.xaxis.set_minor_locator(AutoMinorLocator())\n",
    "    ax1.yaxis.set_minor_locator(AutoMinorLocator())\n",
    "    ax1.tick_params(which='major', direction='in', length=6, width=2, colors='black', top=True, right=True)\n",
    "    ax1.tick_params(which='minor', direction='in', length=4, top=True, right=True)\n",
    "\n",
    "    ax2.xaxis.set_minor_locator(AutoMinorLocator())\n",
    "    ax2.yaxis.set_minor_locator(AutoMinorLocator())\n",
    "    ax2.tick_params(which='major', direction='in', length=6, width=2, colors='black', top=True, right=True)\n",
    "    ax2.tick_params(which='minor', direction='in', length=4, top=True, right=True)\n",
    "\n",
    "    fig.suptitle('WASP-2: Radial Velocity Curve Fit')\n",
    "    ax1.errorbar(time_wasp, vel_wasp, yerr=vel_err_wasp, fmt='o', ms=7.5, label='Radial Velocity Data', ecolor='k', marker='x', mew=1.5)\n",
    "    ax1.plot(xdata1, ydata1, label='Fitted Curve', color='darkorange')\n",
    "    ax1.set_xlabel(r\"Barycentric JD$ - 2450000.0$\")\n",
    "    ax1.set_ylabel(r\"Radial Velocity \" +'\\n' r\"$[kms^{-1}]$\")\n",
    "\n",
    "    ax2.errorbar(photophase, vel_wasp, yerr=vel_err_wasp, fmt='o', ms=7.5, label='Radial Velocity Data', ecolor='k', marker='x', mew=1.5)\n",
    "    ax2.plot(xdata2, ydata2, label='Fitted Curve', color='darkorange')\n",
    "    ax2.set_xlabel(\"Photometric Phase\")\n",
    "    ax2.set_ylabel(r\"Radial Velocity \" + '\\n' r\"$[kms^{-1}]$\")\n",
    "    fig.tight_layout()\n",
    "    # plt.savefig('rad_vel_fit_phase.png')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Plotting the Residuals of the Radial Velocity Curve Fit\n",
    "\"\"\"\n",
    "\n",
    "def plotting_residuals(xdata1, xdata2, ydata1, ydata2):\n",
    "    fig, (ax1, ax2) = plt.subplots(2, sharex=True)\n",
    "\n",
    "    ax1.xaxis.set_minor_locator(AutoMinorLocator())\n",
    "    ax1.yaxis.set_minor_locator(AutoMinorLocator())\n",
    "    ax1.tick_params(which='major', direction='in', length=6, width=2, colors='black', top=True, right=True)\n",
    "    ax1.tick_params(which='minor', direction='in', length=4, top=True, right=True)\n",
    "\n",
    "    ax2.xaxis.set_minor_locator(AutoMinorLocator())\n",
    "    ax2.yaxis.set_minor_locator(AutoMinorLocator())\n",
    "    ax2.tick_params(which='major', direction='in', length=6, width=2, colors='black', top=True, right=True)\n",
    "    ax2.tick_params(which='minor', direction='in', length=4, top=True, right=True)\n",
    "\n",
    "    fig.suptitle('WASP-2: Radial Velocity Curve Fit Residuals')\n",
    "    ax1.errorbar(xdata1, ydata1, yerr=vel_err_wasp, fmt='o', ms=7.5, ecolor='k', marker='x', mew=1.5)\n",
    "    ax1.axhline(y=0, color='darkorange')\n",
    "    ax1.set_ylabel(\"Radial Velocity\" +'\\n' \"Residuals\")\n",
    "\n",
    "    ax2.errorbar(xdata2, ydata2, yerr=new_y_err, fmt='o', ms=7.5, ecolor='k', marker='x', mew=1.5)\n",
    "    ax2.axhline(y=0, color='darkorange')\n",
    "    ax2.set_xlabel(r\"Barycentric JD$ - 2450000.0$\")\n",
    "    ax2.set_ylabel(\"Radial Velocity\" +'\\n' \"Residuals\")\n",
    "\n",
    "    fig.tight_layout()\n",
    "#     plt.savefig('RV_residuals.png')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Fitting the data and the photophase plot\n",
    "\"\"\"\n",
    "\n",
    "fit1, pcov1 = func_fitter(time_wasp, vel_wasp, vel_err_wasp, guess)\n",
    "fitx1 = np.linspace(time_wasp[0]-0.5, time_wasp[-1]+0.75, 10000) \n",
    "fity1 = my_sin(fitx1, *fit)\n",
    "\n",
    "params1, errors1 = radvelparams(fit, pcov)\n",
    "\n",
    "photophase = ((time_wasp - transit_time/2) % (params1[1]/(24*3600))) - (params1[1]/(24*3600))/2\n",
    "\n",
    "fit2, pcov2 = func_fitter(photophase, vel_wasp, vel_err_wasp, guess)\n",
    "fitx2 = np.linspace(-1.5, 1.5, 10000) \n",
    "fity2 = my_sin(fitx2, *fit2)\n",
    "\n",
    "params2, errors2 = radvelparams(fit2, pcov2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEdCAYAAABZtfMGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABkW0lEQVR4nO2dd5wUVfLAv7XkHCQqYckCgmEBUTxZcwJzOs/Aeup5nunuTHfqipjTmX7m08V4xvMI5sCiIqKAiOQMEiRIzmHr90f1LLPDzGzPbM/s7PK+n898Zrr79evqMK/61atXJaqKw+FwOByZRlZ5C+BwOBwORzScgnI4HA5HRuIUlMPhcDgyEqegHA6Hw5GROAXlcDgcjozEKSiHw+FwZCROQVVAROQjEVER+UFEDvLWdRSRnSJS31vuJiKTRGSWiBwWtu+/ReTDKHWKiAwWke9FZJSIfCsiT4tIU2/76yKyVkQWikihiPwoIh+LSHaUuuqIyB0iMlpEvvLKH5LA+X0hIltFZKa37yQRGS8ihyZxrap4dWhIVhE5QkQ+9bFvH+/YC2JsH+PVO0lEuoSt/5uILBeRUXHq/q93jrmJnlNEPdkiMq0sdUSp8woR+SbsOXhTRI4L8hhJytVDRL7zrvk4776GPjO8MheIyL/LW1ZHQKiq+1TAD7AYOC9s+Y+AAieHrfsTkBu2XBWYDWwHGkfUdwnwDVDdW64H/Az0CitTCNzt/RbgNWBMFNlygbFADW/5cmApUCuB81sAXBa2/DywoAzXS4HsMNkb+NwvN95xgfnALVHWf+rzHHP9yFFKPQ3Dfg8CCstQ113AF0A9bzkLGAKsTdWznKB82d697BixvtD7rgLUDVs/FBhc3nK7T3If14OquIzGGs8Q/YCPI9YdCnwXtnws8DSwEjgjor4+wLequh1AVTdgjdXqaAdX+/e/BRwe6rWFsQy4S1W3ecv/AVoCB/g4r1gMB9qKSJMy1AGY7Kq6rqz1eLwDnBe+QkRygIkB1V8qqro2iHpEpAPwT+Aq7/6jqkXAHcCSII6RQm4CUNVdqrqxvIVxBINTUBWXQqB/2HJdYETEuhqqujVs+UxMWbxNRKMKLATOFZHuoRWq+raqzosjQzXsbXZn+EpVnamq4WbEmt73qjh1lUZVYCOwDkBE9vPMZF97prbB4YVFJE9EpojI5yJyWdj6piEzUdi6M7x6RnkmyX4JyPU2cJCIdA5bdx6mvBGRi73jjRaRN6Io83CZQ2W/EpH3RKR52LZDwsxZ34rIVd76L0LmSxHpD9ziyVMoIk96568iMtIrf5GILBWRd6OIcAbwi6rODF/pvYz09fZ/Jdw0KSLDIsynIdPljSIyQkTWi0i+iGz3zMJdxMzPU0Vkhoi0EJFm3n5feed2egLXHxEpVNXvPRNgsUlWRK4DTgQGedfjj4nU68gAyrsL5z7JfYDOmHJohpk9HgC6AzswZdURuDOsfDXgTe93H69ck7DtDbBeWREwCrgCqB1xzEJ2m/iqYz22D3zIOii8HNAKWAG0irPPAjwTH2a2eRw4J2x7J+CksOVRwDHe7+7AZqCdt/xnSpr4svHaXW/5QjyTp7dtUdi2XEoxLQJzgdvDlj/1vvthSrmpt/wQ8O+Ic8z1fh/hXZNQ2duAL8LuzUqgv7fcBpgaVk/4uQ0izMTn3fdVwOFh60YAEuU8niGKyTbGvcmNdvyw7S95v48CTgH+C/wzrMy/gE6h6wUM8X63BH4Lry/i2Nne8cZhz2NhxPmWuF84E1+F/rgeVAVFVWdhprT+3mc0MA1YizV2R2J/3hDHYQoFVf0ea0TOCqtvnar2xxqUOcDDwBQR2S/i0BeJSCE2XrUIuCienCKyD3AdcFXY6s3ADO87HjeKyDfAGsw8OCJs2y/AMd4bdyHQFcjxtp0NjFXV+d7yW6Uc5yegwDvWUKC1iDQrZZ9wis18Yo4c33vrBwEjVHWlt/wG8AcRkSh1XAKMDCtbABwtIm2AAcAGVR0NoKqLsBeIUlHVHcCbwMWefAcDk9VrvVPIMO/4o1T1A+CVMBmqYC8Ps73n6zjgRa/8MuzZ+n0p9f9BVXNVNTdF8jsyAKegKjahcah+wDdeo/MVprD6UnL86Uzg8pCZCOsBnRtZoaqOVtXLsbf0NZijRTiveg1DH1W9QlVXi8gtYeanE0MFRaQG8CpwjaouDDvGalU9UlWjjm+F8ZCqHoEpzSMwZ4sQt2BK+BivkfoYqO1ta0mYOdHHcYYDX6vqEWENXu045SN5G+jumUfP9ZbBeopHh13zJ4HlwD5R6miF9ZJCrAxbH7kNVR2TgHyvYObbGtgLxasxys0BIl9IkiVyjO8DYB9PgZ+A9ZrAzg3glbDr1AGo4/dATklVXqqWtwCOMlEIXAv8rKrrvXWjsbfPOaq6BUBEqmPmuuKxFa+hGCMizVV1uYhcC3ymqtPBBt5F5AvMvBQXVb0fuD98nYhkAS8Dz6nqNyLSABsTW5HoSarqBBF5AbhWRP7PU8R9gK9C54iZskIsw0ygIVmiKYTQtpCJ9GNvuVqssnHkmygic4DzgQNUdbK36Rdgnqr+Jex4TVQ12ljcL0DTsOXQ78VRtoV6Qj+pOTGUJt/3IrIcG2PqrKozYhT9L3C/iHT2euihY9UD3gVOUxvT3A7U8LY1LO34ngw7ROQtrKdYB/irt+kX7/vsUO9RRGri2iYHrgdV0RkNdMNMbeHr+mAmvBDHA19H7Ps9NuYRMvM1Bv4UMj95jdIJWI8sGf4P+AH4QkTqAocDJ3t17+MNiDdOoL6HMEVyirc8B+gtIlkiUgfrYYV4FzhMRNp7yxfEqfc3zCwammN1YuyicXkHM2WG91qHAqeISCMAsblSI/bctbjsybLbS/ES4EvPnDcSqCciR3r1tAeeiaGcNuD1/jxHi1BD/yo27hNz/pdnEr0DeNq7ZyGF/SgwRXc73Mxnt0fmybHqi8IrmBLPCvVqVXWpJ1O4qfhZrNccBBuA2mJz814PqE5HuijvQTD3KdsH6y2cErYsWKMbchjojymricB+YeVu9cr9io1nHIx5+H2L9czGAzeGlX8da8gXAgWlyHQcNpAd+RnkbW+NmayiOklg83C2AjOBv4WtfxnrUdwNtMAcIyZ5so3yzvMCr2weMBX4ErjeO/53QFvvW73zrAKcDswDPgHuCSvbx6t/K/BOKed8oLdf94j1F3p1fYn10jp76//r1TsJyIko+5W3vXlYPTmevKO9ug4Iu1YhefcDGnn3egzweNj+rYFtQDMfz9TlXn2F2Hy22zClEtp+hHdvPscUafjxXwk7rwui1D0N64mFr2uGvVR87X3+GUOuHmH3bhxwUZTtJe4XcBg23vl9NHncJ7M/4t1Eh8NRifFMcS+p6pnlLYvD4Rdn4nM4KjEiEjIb/h7z5nM4KgxOQTkclZsOmNv2SZjZ0OGoMDgTn8PhcDgyEteDcjgcDkdG4hSUw+FwODISp6AcDofDkZE4BeVwOByOjMQpKIfD4XBkJE5BORwOhyMjcQrK4XA4HBmJU1AOh8PhyEicgnI4HA5HRuIUlMPhcDgyEqegHA6Hw5GROAXlcDgcjozEKSiHw+FwZCROQTkcDocjI6la3gJkOk2aNNHs7Oyo29atW0eDBg3i7h9EmdK2T5gwAYCcnJyUHieIMkHJ6qdMEHX4kdc9B4mXSedz4KeMu7bJlQnqfCZMmLBKVZvusaG8c85n+icnJ0djcfnll8fcFmSZ8O2jZizXoqKiEtsBBXTUjOWBHSdVZUKypvo4QdXhR97yeA6iEdS1rWzPgZ8y7tomVyao8wHGa5T215n4KhCFM1cwqOAHhoychkYkmmx0zOUMKviBwpkrykk6h8ORbgpnrtijLQjfVtFxCqoC0b9zU/L6ZVMwZkGxklJVGh1zOfV7nUZev2z6d96zl+xwOCofe8MLqxuDKoV169ZxxRVXFC8PHDiQgQMHFv8ujSDKhLaLCPkDugFQMGYB77z7LgLU73Ua68cPI/++5wB7cHO7NEv6OOkok47jVCRZ/ZQJQtagjuOubeqO47dM+AsrUNw2+H1hLc//2IgRIxgxYkT4qqiDVBKre+gwevXqpePHjy9vMUpQOGMFg4b+ULzc87dhbCx8gS9mFTFk5DQKxixgaF7vqEoqLSz6EuZ/BA3aQ/dLoFrt4k0iAhDTLJEy1i+CqS+D7oJuF0HDDr52S6m8u7bDtFfht+nQ5mhodxJ4x0uGwGRd/A3MHQ71WsEBeVC9Xtnqi0LCsm5cBlMLYOcW2P8C2Kdr4DLFI6q8RTth+huwcjLs1w86nl6m+5cMqlr8n8/rlw3Yy2uPVcMYfmVX5IA8qNkwrTIlg4hMUNVee6x3Cio+maigioqKGPDEaKb9uhmAQTWHk1/nee6q9ywFC1qR1y+b/AHdiv9UaeXrf8D390NWVfsDNz0Izvkcau2TfllCLP4G3j8Ftm+wBqRKTTj1PWh3YvnJtH0DvHsCLBu7+1od+Gc45qm0N3IlGHcvfHPrbpka7w/nfAl1W5afTMsn2LXa+htIlsl2yn+g05nlJ9POrfD+AFj0xe5r1fVCOOllkzGNhCspgLxaw8iv+xLCLqifDeeOggbZaZUpUWIpKDcGVcFQVe4aOZ1pv26mW9V5AAzdeiq/W/G0Kacum8pPOf38kimnHpfDNRvh9BGwejqMPA+0KP3yAGxYAsNOgzot4I+z4fKF0KgzjDgH1swuH5lU4aNL4NfvraG9djPk/B1+egYmPlY+MgHMfMeU0/4XwNXr4JwvYMMvMPwMa4DLg82r4P2BUL0uDJoGVy6DZjnwwQXWcykvPr/SlNNxL8B1W+CwwTD9NRh7V/nJFKJBe/jLb/D7b2HbWvjfQNi5rbylSgqnoCoQxW9K3y4gr9YwRp5VjbzDswFYUqUN3Wot5/YNVyBbVqZfuM0rYPTfoXUuHPsMVK0BHQbAUY/bH3nqK+mXCaDweti5Gc4YaWa9eq3gjBH2lvv5VaYs0s2cYTDnfTjiXtj/fKhSDfo/BO0HmoJYtyD9Mm1dC19eA817wYkFZpZtc7Q1wMvGwaSn0i8TwNe3wJaVcNowM+vVbganD4Pq9eGzK8rnxWfhF2YuPvRW6HmZ9aAOy4euf4Bx95jJNk2E957OrTqSs6p8QMGvPRjy2WK0ZV84+TVYNQV+eDBtMgVJqQpKRE5NhyCO0hk9a6XZmut/Tn7Hn5EDL4ewjtK0Lc0Zvalz+bzFfXcP7NgIxzwDWVV2r+95ObQ8DL75Z/rf4pZ+B7PehT7/hEaddq+v1wr63Q2LPoeFn6VXpqJd8NWN0OQAyPnr7vUiZt5DYeyQ9MoEMP5he8k49hmoUn33+v3Ph7bHm0zbN6RXplVTYMqLcPB10OzA3etrNzWFvmwczH4/vTKp2otYg/amoEKIQO6jptjH3J4mUcLGnzr+xgMNn+XhP51T0tO33cnQ+Wz44QHrjVYw/PSgbhORB0XkgJRL44hLbpdmDM1dTX71x+DwwQz5YAYFYxZQb9l46i2zcbKv6l6A/vxva2zSxeZV8PMLZoPfZ/+S2yQL+g2BTctgWpp7UT88ADUblVQEIXpeAXX3g+/vS69Ms/8La+eYSahKtZLb6reGHlfA9FfT24vavsF6SJ3OhBYlhwEKZ61ED78Ttq6GSc8Ur1fV1Lswj38YqtaGQ/+557auF9pLx3d3p7cXvPBTWPkT9L0dqtUqua12Uzj4Wpj9HqyamnJRil9YD2tN/o5bkLbHIm2OIn9At2IlNXrWSjj8TtixqXzNx0niR0FdAgwGjhSRp0QkGH/MCkLIzTz0iXCNTC+q5C69B5rnMOTnVsWeO5OH5jN5aL49lEs6MWTdReiPz5ReX1D89LR5V/W+Kfr2NseY6Wj8I+lrTFbPMlPaQVfb+EUkVWtAr7/DL4WwfGJ6ZAIY/5A1rB1Pj7699412jSY9nT6Zfn7Rxip631hidfE8mwn10NZHw6T/g6JdxW/uKZ1ns2EJTH/deuC1Gu+5PasK9L4ZVk6CJV+nRoZo/PAw1N0Xul4Qffsh10HVWvDjkykXJbdLM4bm9Sa/3Xhk8zLoczOwezpKsSfvPt3sefvp2YwZixoxYkSJdpUYbualhvoBDgBqAH8EfgRGAk8C55a2b2X4xAt1lHYWf6P6MDrqk9e07c0jdfDwKSXCHhUVFeng4VO07c0jddRjJ6kW7Uq9TEW7VJ9vq/rOcfHL/fyS6sOo/vJ16mVSVR19s+ojVVQ3LotdZssa1cdqqn52VXpkWj7JrsGEJ+KX+98Zqk81Vd25LfUyFRWpFnRXfb1vlE27n6fBL/9Pix5Ci+aM3L0u4vkLlLF327VaMzd2me2bVJ9soPrBH1IjQyRr5ppMY++OX+6jS1SfqKe6bUNaxNLX+6oWdLN7GYv5n5js0/+THpkShDKEOnoNmAb09JTSAFW9xlNcjnQypQCq1SE39zR7c4rw1it+czp2M7k7P7L5SKnml0JYvxC658Uv1/kcqFbXxhRSTdEuM5O1O8m892JRs6HJNeN12LE59XJNexmyqsV++w7R83JzDJgzLPUyLZ8Av02F7oP22BR6nvL6ZVMwrSq3rv8TVz73SXHPPWXeoqp2rVr1h4btY5erVttMfbPehS2rg5cjkmmvAALdLo5frscVZjad+VbqZVo9E5Z9Z/cv3r1oe6y5nP/8QuplChA/CmomcJCqXqeqswFEpDrQKKWSOUqyY5M98J3Phep1ye3SLGrjICLk5g6Emo1NoaWaqUOhRoPYJqsQ1evagPvMt+1cUsmiz2Hj0qiN7h4c8EfYti71ymDXDjNZdRhY+pywtsdD3Vbmtpxqpr4MVWpAl/Oibg4pqeO6NeON7QP5JOsE8g5tUaycNBVjUcu+sykA3S8pvewBl8KubTbuk0q0yBRUm6NtrDAe+x5m88fScf+mvWLjvF0vjF9OsmzS9aIvYcPi1MsVEH4U1HBV3QAgIgeJyL9UdbvXiyo3RORREblHRP4hIm+KSAtv/SARmSAihd5nnogMirL/CSLyhojcIiL/E5HOaT+JRJj3oXnJdS/l7Q2gak3odJZFA9i5NXUy7dgCs96zxi1ywDgaXf9gLt/zPkydTADTXjPniPYDSi/b6nc2pjDrndTKtOhzc1zp5qPRzaoCXc6BBR+b8kwVu3bAjP9Ah9PiRhsYPXMln00LU0LrbP6dpmosatprNo7T+ezSyzY72Mb0Zr4d3PGjseRbWDffn9IUgS7nwy+jYdOvqZNJ1a5V2+P9TaQOvYSkWpkHiB8FVRwTRlUnAeU043IPNqnqrap6HzY2FvL5nAmcoaq5qpoL/Ay8G2X/x4CHVPV+YDRwY5QyTJgwARHZ45N25o2wXtF+R/gr3/kcU2gLPkmdTL+MMoXjd0b/fr+DWk3NJJMqinbC/A9sTlHVGqWXlyxT5gs+gu0bUyfXnGFm4sw+wV/5zudYKKS5w1Mn09JvLTpDl3NjFlFVCmftVj7dqv9CwYza3DliKkNG7A6xE1iQYlWYO8zMs35CLImYVeGXL2FzCuf/zfmfmWc7nOavfJdzALUXuFSxYhJsWBT3/pWgcRdoemDqlXkcorWl8drTmApKRK4TkfnAX71eyHwRmU2GBJhV1dvCFrOAjd76saq6CKzHB8xS1Wgtz3Ig9K9qCkxInbRlpGiX9TranWyTAv3QOtcUWiqVwbwRUK0OtMr1Vz6riimz+R+kbsxn6bewdY2Z0vzS+Rzrac4bmRqZVK3u7OP9KU2AlodCvdYW3SFVzB1hc56yj4+6OdRDGvrtQgYd3pa8w7OZtr013arOY+i3C23C+OEBj0WtmAQbl9gLhl+6nGsmuNn/DUaGaMwbaf+pGvX9ld+nG+zTHWalUBnMGwkItD/F/z6dz7H/yPpfUiZWkMRUUKr6uKq2A65Q1faq2k5VO6nq9ekSTkQ+EZFJUT6nhpVpCBwPPBSlimuA/4tR/XVAvog8AhwORJ2xmZOTE8u7MX0sHWtvuok0ulWq2bjQ3OGpcS1Vhbkjzbzgt9EFM9vs2GTmq1Qwd6S96baN3uhGZb9+UKdl6pR5Mo2uZNm1WvhJ6ibIzhthLxcxeirF82z6ZXPHwO7kD+xG3sG1mbZzt+PCkZ2bBGtRmDcCa3RP9r9Pkx4WvipVpqs1s2HNzMTuH5gyWPw1bFqeGrnmjbAXmdoJBIUO9bbmpHmCs0e0tjRee1qqiU9VS7zCiUj00dQUoKonqOpBUT7DPVkaAE8Dl6pqCTceEWkO1FDVhZH1ek4eHwDnq+rfgUeA/6T8hJJl3kjrOfk1D4XodCZsXw+LC4OXacUk2Lg4MaUJXs+uUeqcEuaNMO8vv2+6YMqgw6lmDk2FMk+m0QUzJ+3anppoF6tnwZpZccfpiufZhPeQau1uDI/r1jz4iPlzk2h0Rez+/VKYGmU+15v72MHHmGY4HU8D1CL7B83GZfDrD4n//xp1gsZdU2ctCJh4Jr43vO/5nolvnmfyS+MMwtiISBPgKeBGVZ0vImdFFLkKeDZin9YiUgWb19UYCBmtlwE1Uyxy8hQ3utHnssWk9dE22Dzvg9TIlKh5ATxFe6L9aYOOo7ZmDqyekfifFqyh3rERFn8VrExgjUGijS7AvodDjYbWK0yFTFDqtQp5i4bMfQXfLiCv9Tzy6nzIZ9OWR02WlzQbl8Hy8cnfv6IdqVHm80aaua5Bu8T2a3qgRStJhTKY7zka+XEEiqT9gNQp84CJ14Ma7H3/yzPxtfdMfrfF2SedfIrNzXpdRAqBQaENIlID6KWq30Ts8xZwmOeV+FfgORH5B3ADptAyj7Xz4Ldpyf1pq9Uyt9h5I4OP4JBsowv2B9my0t4AA5Up9KabxLVqc7R5PwbdmCT7pgtmpm13ko3ZBa3M542weIA+0jAUK6fQ/KcT2pFf62nyelQpkd257DJ5L1KJmtLAzLQ1Gu7u7QTF1rUWqSKZ+ydiz3oqeuZzR0C9NmbeTJQOnjJf8GmwMqWAeGNQs7yfL4lIPQAR2UdV0xhDJzaqeoiq9gx566nqwLBt21R1j1d7VT08pLRU9TlVHaSq96nq+VGUWWYQanST+dOC/UHWzbeeRVCUpdEF60FJVvA9u7kjknvTBS969zHBK/OyNLpg92/zimCV+da1Njbi8+07fCwqf0A3pO3RSLWa5Lf8uGTMt7IyL9ToJhEDIKuqKfN5ASvzBR+bZ2gyPRVITc9851brKXYYmFzusFDPvAKY+fxGkjjO+91fRDKlB7V3MHeEeQTFm1Efj3beuEeQyqCsjW6txvYnCfIPUpY33RDtB9gcn9UzAxOLeSOgftvkGl0IU+YBXqsFH1tmYZ/3b4+xKE+Zy7wR5J/SNZjszTu2wMLPk290YXfPfNn3ZZMlnHkjoVYTaNk3uf1T0TMPTe9IVmmmSpmnAD8Kapyq/hfA+96rUvCWa7DYbetg8ejkFQFAfc8MMD9IBVXGRheg3Smw4keL+BAECz4p25suhCnzgBqTUKPbfkDyjW6txrBvv2DHoeaO8BrdQ33vskfkEk+Zy5pZwThKhBrdsrxgZJ8IUmW31aGsFO20sZ52J5dMIZMI1WpDm2NNpqB65nO96R2tc5OvI1Vmdp/4DRbrR0FFxmUpx9zd6adBgwY8//zzxZ+BA9MYzD3U6JblTwvmyLDkm2CiEuzYYuaFsjS6IZkguKgS80aU7U0XTJk37RlcAxdEowt2rVdOCmbuStFOm5RclkYX7AUDgrtWic6pi0atxjYWFZS1YOlYm1NXlpceCDOzB5DIMDSnru3x1jNLllDPPOgxO58MHDiwRLsKRG2c/CioWSLysxcO6CcscGxMRKSxj0/DhM9ob2RuAI0ueB5OO4OJKvHLKEutUdZGt8kBNt4QRG8liDfdEO0HwJIx1jCVlXkjLHpEWRpd2H2t5wegzJOZyByN+q3NSy2I+5fsnLpotB9g+ZqCUOZzR9icukSnd+whk6fMg+gFr5wMG34pu9IM9cxT4eEbIH7mQb0AnIONRZ2rqv8uZZelwHgsMkOsz+QyyLx3EGSj27KvRZUIojEJqtENeTgt/Kzs8QKDanTBZNJdZVfmoUY3kegRsWi8v2VwDaK3Emp0E5nIHIuQMi9rJPGVPyU3py6qTCFlHkDDm8ycumjUawVNDwru/5fM9I5odBgYXM88RfhJ+S7AEUB/4Agpfdr49LDIE1E/wG8ByF65WTrWspiW9U0JTMG1O9nMaUW7kq8nyEYXzN1152abk1EWgmx0W/SxeIFlNX2EGt0g7l9ImS/6ouwhouaNDKbRheCU+Vyv0W2X4ETmaDTuAg07lP3+Fc+pC+D+gSmDpQEo87kjoGUfqNO87DKFns0gx6cDxo+J71HgEGA+0MtbjsdhPur0U2bvJijzQoj2Ayxc0rLvkq8jFD2iLE4b4bQ+ylJ6l7UxmRvQmy6YMm9/so3TFO0sm0xBvemC3b+dW8uW46ssE5mj0dJT5mXtGcwLsNEtoczLkNYldE5BPevtB5jH3IIyRJXYtBx+/T44mYp75mW8fyN/D98ODkSkSPwoqJWqepWq/ktV/wysjVdYVUu11/gpkymUmxdfUOaFEO1ONPfSsjyMyYbsiUXVmtD2uLLNPQrFSQuq0QVPma+xXmyyJBMnLR6tjjTTapnvH8H06sAG2dufUjZlvulX8yQLqtEFO79d28qmzOeVcXpHJC162bNQlnGo4ukdQd2/AHrm29ZbDMSdWxLaLUgvvsgdfbeYInKz37KZSrl48a2ZHeybLliYpFZHlq23Mne4jWcF1eiCneOGRbDq5yRlKkP0iFi0Pd56r8kqg41LyzaRORpVa1hvuizKfK4XPSKoRhfClPm3ScoUCrl0avxyidDqSAuAm+z927rGcjkFKZNkmefjgo8tD1cyzB1ujkVNewYnV4eBXs/8i+T2X/CxRaVI8FoF6cU324sgXqoXn4i8HfZ5B7gsIakdRiiQakefuWf80n6ApfdeNz/xfTcstvTgQctU1rlHc4fZPK9kokfEokZ9670mq8yL49yl4P5tXGKm1kTZusaiGQTZ6IKNR2ZVS75nMHeYpSIvy5y6SKpUL5syn+9NZA76WnUYCNvW2lhUouzYAgs/NZmCjB5f1p753OHBeBrHwK8X3+/x58W3XlXP9T7nAJ8HJOfexdzh5sJbv22w9YbMKMk0JsU9lYD/tHVbQoveySmDLb/Z/K6gZQIbHF89HdbOTXzfOcPMtr9Pt2BlancSIMk1JvM/TE2jW72eTRhNxsNwxybLNBx0owueMl+anDKfO9ysBC36BCtT22NNeSbz/1v0hTe9I+D7VxZlvmuHmR3bDyi7p3EM/PSgUNXpqvquqs4UkSviFL0HiiONw+4st3svu7YnVn7zKnvDCvrtG6BRR2jUJbkGbu5waNjRBlaDpv0AWDbOYs4lwvwPbeA56F5dSCZIfJ7I9o3WmKSi0a3T3JwJkrl/c4ZD7eb2MhA07QeYSXrNnMT2W/i5mZdS8YJRrMwTVJy7ttuYWioa3er1bHpGUv+/Yd7LQP9gZQLr2SWjzJd8Yz3CVNw/j3jpNlaHpdkIT7dxX6x9VDVkO3rJWy6jT2UF578nw0eXJLbPvJGpa3TBHsbFhYmF2t++wVJqdzgt+EYXvJ6dJh5VYs4wSzTYPCd4mRp2sLw5ifbsFn5mA/Spun/tB5gnVyJJ8EKNboeBNhaSCpkgcXflucN3j40GTe1m5qSSqDJY/LVFXElVo9t+gDn1rJntfx8tsucw+yTr8QRNssp87nCoUiNmRuYgiPe0Xh2WZiM83cY1PupNQStWAanf1m5iHHfX3NxccnNzd6+YO8xSfTc7ODUydRhoDVYiSdTmf2z7dEzRn7bZQZY3Z87//O+zY4vNv+lwamoaXbDGZPHoxKJKzHnfIkXv2y9FMoXMtAk0Jr+MspeMVDW6DT1zZiJJKIt2hjW61VIjV/sB5qySSLzHOe973qXHpkam0LyqucP977P0O9i8PHX/v5AyT0QmVfu/tjnGQlSliHjpNt4I/RaRniJyjIi0wl/m2UoTULZMbuZdzreJqH5tzlvXmuLoeEZqeipgDWedljAjgQTCM98081CqGl0RS4+94CP/ymD+h5bGoFNknsoA6XKueSjN/q+/8ju22J+205mpa3Sb9rTe3cw3/e8z4z/WU2l7XOllk6XT2Tbh2q8yWPSlBSsNpSBPiUzeszHzbX/li3bCrHdMsaWq0W3Qzl4+ZyR4/6rWDNYVP4zCmSvQTmebE1RYz05VKZwZw+y+bBysX5D0/QvMzVxEbgQeBy4GOgH3+zh+pelBlcnNfL8joO6+/huT2f8181DXPyQnrB+yqsD+51sD70cZbF1r4zD7n5+ygVAAul5gvTS/ymDGG6Y02xydOpma51iK7BlvlF4WzMS1fQPsf0HqZBKB/f9gDbwfZbBji13TjmeWLbhoaXS9AFD/De+MN6B6fc+8lCL22R+aHQLTX/dXftGXNg66/+9TJxPY/3v5eFg9q/SyRTth1tvmoh7UnMgwCmeuYFDBDwxZdASqAtPtWQ8lqRxU8EN0JTXjDTPvdTwjqeMG6WZeV1WPAqap6ijAT2rIf/gXtRKTVQU6n2vKwE8k8Rlv2NtxKgayw9k/AWUw5/3UK02A5r1MGfhpTNKlNEXsWi0a5U8ZzPgP1GlRtjQIfggpg5lvlV52/oemNLumUGmChRhqnuNPme/cCrPftx5OKpUmJKYMZvzHU5oBTUSPRZfzAfH3rC8aZUozRfevf+emlnTyh1UMkcHo9NfRoqISGZT7d25acqeindYrbT8gJUozHD8KKtQCaMRyTFR1StISVTb2/70pg9LeLNf/Ym9w+1+QOvNeiFDPYNqrpZed+rJ57zXvlVqZQsrgl8LSg1fOetuUZqrfdME7hha/WcZk8ypTml3OS63ShN3KYNprpZed+rL1NFsflVqZwO7f8gnwWylpJeb8D7avT73SBHuJQWB6Kddq+waLiNApxT1NgHr7QZujYMbrpScMnPZySpWmiJA/oJspqZU5DFmcS+6Vt5fMoBzZHi341MbE0vD/86OgdonIx8AAEXkX8BXgSkRGi0h97/eVInK9iATmgiIij4rIPSLyDxF5U0RaeOuri8hQERksIs+IyO0x9s8WkRe9/Z8TkbpByVaCFr0tkvFPT8efZzD5Wfs+4NKUiFECEeh+qTkArJoau9yqKVamx+WpV5oA3T2Px9C1iIYqTHrK5okFPU8lGo272NjbT8/ED7Q75SVTmj0uT71MYM/Jiok2gB6LdQvMi63HZalXmmC9laxqMOnp+OUmPW3zxFJpng1Rd18L8/XzC/GnfEx/3ZRUzz+lXiaA7oNg7VwKv/4EjdIuqCqFk2fbmFj3S1KqNIuVVN99KdhyGgsbHx5bOYG1ZXVaBBspJQbx3MyfAlDVO4B/AcOB51T1Xp91N1TV9SKSA1wONAJeKKO84WxS1VtV9T7gR3bPuToDaKSqg4G/AH8Tkf2i7P8sdj73AVOAqGGZJkyYgIjs8fGNCBz0F8vjsuSbqEWqZRXB5BfshjfI9l93WehxmT30Pz4Ru8ykp83OnA6lCTaA3OFU+Ok5GzuJxpIxdi0P+kt6lCbAIddaKvhY+ZiKdpkCa50LTbqnR6ZuF5vjw8THY5eZ/Jxdo3Q1unWaW49l6tDYJu2Vk2HJ13Dgn1PnfRnJIddZzL9YzhKq9qw3OyShLMNlovO5FMoxDPrQzGnhSqp4/OeNWRRuOcCuVToId2GPFZtv7TybDtLj8qRc3qO1pfHa03hPSLbXA7kZ+ElVH1bVzxKQZYeIVMWcKx7wFF1g/15VvS1sMQvY6P1eDoQmCtfH8lOVmI8lItWAo4BQvuMxQEBhp6PQ9QLLxzQuum4/pc0y82g6+NqUibAHtZvYYPu0V6OPr2xYYr2Crn+wsunikOss6vrUgujbx90NNfdJj3koRMczoG4r+P6B6L3g6a+bR1M671/1unDAH2H2u9GjXWz5zRrdTmdZcsF0cch15l056Zno27+720LrpOulByy2YuP9YfzD0U1qc4dbLMhDrk3bS0/h3HUceeiR5NUaRsGYBcVKqqioiCtenWAmtnqf0r9jE9ina0plCSnEgjELyMtpaDJ9v2IPxQlYG1alWtpeeuIpqAu96OVvAX8Wkf8TkURyPzwB/AQMAEK+2QmZ0UTkEy8OYOTn1LAyDYHjgYcAVLUQmCgirwBvAi+rauTreBNgi+6++uuBqBFQc3JyUNU9PglRrTYc+k8LrBgRYblWlZ1c3GmhxX5Lh8kjnEP/aQOeY6JYQcfeaX/mvrftuS2VtM61iZtj77RIyeEsGmVzn/rcktK5F3tQpZpdh6Vj9pyrtWMLjB1srsOpmpwbi143QFZ1+PqWPbeNu89MVofdkV6ZmudYhPPv74PNK0tuWz7BTFY5f7WMrulCBPrmW46uyHHXXTtgzG3QqHPqHYE8Qp5zd604hdsbv0tesx9NSY2YxoAnv+Gzacs5rsUa8ms8gRwxJKWylFBO/bLJP/tw8nM2kVd7ZAnFCdjY4tSh0PNKG0dL8ngJtaexdohSQT3gf8C8BPapC9TyfncECvzu67P+BsAbQLuwdddipjswBfwDcHLEftUwb0Txlg8BJkY7Rk5OjgbCji2qz2ervtBedetaW1dUpB9c20J3PoTq0u+COU6ijPq76sOozvto97r5n9q6UX8rH5mWfa/6sKh+NGj3uq3rVP/dQfWFdqrbN6Vfpl07VAu6qT67n+qm5bvXf3mdXasFn6dfJlXVb++04898Z/e6xWNUH8lS/eSy8pFp1TTVR6qoDjtTtajI1m3fpDq0h+ozLVW3rEm/TEVFqq8fqvpUE9V1i3av/ybfrt+s99MoSpEOHj5F2948Ugf/+1Xd9SA6+KW3tO3NI7XtzSP1pIc/0V2PVFcdfm7KZRk1Y7nJMXyKFoXu1bqFWvRoTR38yF3a9uaROmrGctWd21Rf76v6f41VN/4auBzAeI3S/sYbg2rmfXcXkceBBdg4UiKv1NvU672o6hxVzUtg37h48f6eAm5U1fkiEpqx2RpY5h2zCDP51fT2aS0iVVR1BzAKCPlz9wNSm1ayak04+TVYvxDeO8HikH1xNSe3+ZU3ZrdJn+07kn5DLBr4yHPN42vqKzDiLIsM0O/u8pGpRW/oe7u9rX3+F1j4hV2z9QvhxJetR5pusqrCSa/BllXw7vHmyTT6JhsDOuhqaHtM+mUC6022PBQ+uhgmP2/eov871aKY9H+4fGTapysc+YBNY/jkUuv5vj/AnG5OeBFqNky/TCL27OzaBu8eZxPivx0M3w2x8bxOp6dRlDDPudmNGJI1hHFzd4euOnTrR0jtJnDM/6VcltwuzRia17ukQ0T9NsjRT5BfdDtDu31Obq2pMOx0S3Z67LPBJJb0SagHsecGkU+B2kAX4GXgBVWd6btikX8DZ2Jef0uBycBkVX2yrEJ79U8EqrJ7fGmDqg70vPmexBwfamHjUNeo6i4R+Ra4SVW/EZFsIB+YB7QB/qaqGyOP06tXLx0/fnwQIhuz34ePLymOhffuvP14amoHRhWODu4YibJhsTVqK3605WYHw+nDoV6r8pOpaBeMvgEmPmbL1evBCQXQOYWRI/yw4BMYeb4FyQRzNjn2GVNg5cWm5TDsDFjmJVhs3BXOGGFz6soLVTMdj7sXUKhaC457DrpdVH4ygaUcGX62jfmCuUqfODQ1Me5KQVW5c8RUhn67EIBuVeZyaLUpFGw9jbxD6pF/zu8Sc8gKVjj44UH45laLgl+lOhz1OBx4ZUoOJyITVHWPuSzxFNRkLDDse6qaYEhuEJGZwAGqusPzojsQ6KmqfiJRZAyBKyiALavh13HQsCO5Z5hbcmFhYbDHSJSiXbuTzu17WPk2uOGsmQ1r50CLQ9M7bhGPbevtWtVvE3xKjWTRIgs/s2s77Ht46kItJcq6+TZ20aI31G5aevl0sGOTedTWaRls8r8EKSoqYsCT3zBtmb2sDjqwJvmHwV0/NaVg7C/xXb3TxfpfzIGk+SHmWp4iYimoeK3Qn1U1icxaxXyHmQRXqOoSYAmQYLjqSkqtxqkN8ZIMWVWg1e/KW4o9adTJPplEjfo2tyaTkCx7scg0GrQLNplkEFSrYzmQyhFV5U+vTWDasg1U27ScGut/YSi9kLrZ3D6wK2RVoWDMAvp3bkpulwAzWCdK/dbp9QKNIKaCKqNyAngeGC0iLwLjMPOej3g/mUUoWGyIgQMHpiftu8PhqJSo5zn32bQVHNetGc9deBIiUuxNB3D7KV3LXzmlkBEjRkQG3o4aLDamia+seLmjXgG24Jn3gJqqWo6G8cRJiYkvjFCqjXI38TkcjrQQcjOPNOGFFFfBmAUMzetdaZVTNJIx8cWqqKWqLvNRdLHa5NzwfWskejyHw+GoTIQ85/p3blpifCnk3VeZe06JElNBicjFMTYNBM7xUfckEblOVYtjsaiqn0joDofDUamJpYBExCmnMOL1oC4Bvoqy3q8rVXPgWC9U0kQsqsQkVX0nMREdDofDsTcST0Fdq6p7hLoWkQP8VKyq53rla2Ax+HoAhwJ7vYIqnLlij+49mA169KyV7g3K4XA4iJ/yvVg5iUhXETlSRI4ErvNTsYjc7dWzTVUnqurLxIgYvjdRnMEyVgTjWBksHQ6HYy+jVCcJEXkQiyaxLzAL8Btadz8RuUBV3/DqaQq8jUURrzAE7WZenMHScycNqai4GSwdDoejEhGYm7mI3Kmqd4jIzar6gIj8XVUfKU0ALznhJ1ivSYEC4GZVTW3Mu4BJhZt5uDtpvWVW94aWvTJj5rjD4XCkmbK4mYeCVDXy8jvllHKgVzCniB+xhIFvADuB01V1TkJSV1JC7qQABd50aKecHA6HoyR+UlruEJGBwHhgA1BaXL6XvXovxZRTNrAGuFBEzk5eVIfD4XDsTZTag1LV/NBvEfmMUhSUqn4BfBG2T1WgGxZNoi/wbrLCVhaimfhCPSnXi3I4HA4jXj6oJt73kWEefAdimXJj4qXBKEZVd6rqZFV9VVVviFZmbyJcOeX1y6bxwlE0Xjiq2HEiapplh8Ph2AuJ14N6FTgJeBwbTwq91vcopc6uXqqOWAgxPDYykaC9+EbPWrk7vfKAbozy3E12j0llQARjh8PhSCFBevH9TlW/DlvuFy/SuYi09SHfLlVd7KNcuZMKL77wibrhwWLdRF2Hw7E3krQXX7hy8ugBxFRQqrowcfH2LlwcLofD4SideMFiVwNrgWpY/L2VQDNgHfBsOoRzOBwOx95LvB7U1ar6hohcA7ykqptEpC5wfXpEi4+IPApsBjZizhvXq+qv3gTh54EFWMDapap6V8S+guWqmoU5inTAMghvSt8ZOBwOhyMe8WLxveH9bB5quFV1I1A/HYL5YJOq3qqq92FOHLd6688AGqnqYGyi8N9EZL+IfbOAeap6l6reCWwCrkyT3A6Hw+HwgZ9IEl1F5EZgNtAZ6JRakfyhqreFLWZhPSmA5UAT73d9YCmwOmLfXUB4MsXw/UswYcKEqPOSnCu4w+FwJEaiczz9KKjLgVuAXGAacFnCUiWJiHyCmekiyVfV4V6ZhsDxwFkAqlooIhO9kEvNgJdVdUucY2QD7YFrg5Xe4XA4HGXBjxffauCm0LKIHAL8lkqhwo59QrztItIAeBq41JMTEbkWqK6qF4tIFjBORKao6odR9m8F3AecFyvbb05ODkG7mTscDsfeSCzLU6yeVamx+ESklYg8ICIvichLmANCueNFungKuFFV54vIWd6m1sAyAFUtwkx+Nb19WotIFe93B0w5/UlVV4ft73A4HI4MwI+J737gPeBELI7erpRK5J9PMflf97TvBkzOR4AnReQOoBawCBjm7fMWcJOIjMfS2S8Bhnv7z/b2dzgcDkcG4EdBTVLV90Wks6p+JiJ9Uy6VD1T1kBjrfwXOibHt8LDFSM8+h8PhcGQQftJt5Hjhi5qKyIVUsIy4DofD4aiY+OlBPQbUBZ4BHgKeTKVAmUbQwWIdDodjbyewYLF77CBynqq+VQbZKhSpCBYbTniwWIfD4dgbiRUsNl4+qOoicp6IHBW27gxgcGpEdDgcDodjN/FMfK8B1YHGItIGCwW0FhiUerEcDofDsbcTT0EtVdXrAURkDHCzqn6VFqkcDofDsdcTz4svPH7dRyHlJCLHpFYkh8PhcDji96BO8NJrAPQRkVAU877AF6kVy5EIoTAhFSGAbUWSFSqWvE7W1FGR5K1IspZGPAW1HUtDATAqbP2O1ImTeTg3c4fD4QgWv27mqGrUD9A7xvqcWPtUxk9OTo7GYvjw4TG3BVmmtO2A2q1M7XGCKBOUrH7KBFGHH3ndc5B4mXQ+B37KuGubXJmgzgcYr1Ha33gJC3+IsX6Cfz1ZuYl4A0hZGT91+CGI4wRVJh3HqUiy+injnoPUlXHXNrkyqZbVT6gjh8PhcDjSjlNQDofD4chIEg51tLchIiuBhTE2NwDWlVJFEGX81OGHII4TVJnSSNd1yxRZ/ZRxz0Hqyrhrm1yZoGRtq6pNI1c6BeVwOByOjMSZ+BwOh8ORkTgF5XA4HI6MxE8+KEeKEZEsYAQwDgvQ2wG4FNgfuA6YBnQHblfVRSJSHXgeWAA0x+Im3iUijYCngZ+A+t7nOlXdFXG8XCzP11pv1Qeq+lC65fXq+h/QMKz6c1R1ZZTj3QtsBNoCL6rqd+Ug6xJgtld1dWC7quZGHC+XJK9tErLWAB4HfgFaAt+o6pteXQcBfwHmA82AG1R1Z5TjJXVdg5Q3Hc9twNf2f6TwmU2BvCl9blNKtMlR7pPeD9aTvS1seRjwB2AScLC3biAwzPt9XtjvLGANlsK+FXBlWD0fAWdFOV4ukFve8nrLg30c73zgae93Y2AWUKUcZD0vrJ484MIgr20Ssl4PPO79roo587QEBJgCtPC2PQL8McjrGrC8KX9ug5I1Hc9sCuRN6XObyo8z8WUAqlqkqncDiEhV7A87E+gELPKKzQNCgXqXA0283/WBpcBqVV2sqs+GVZ2FvcFF4yIRuUFEhohI6/KQ11tuLSL/FJFbReSiGIc8BRjrHXs1sBV7e0yrrFoyUec5QKzEnUld2yRkLV6v1jtaAfQD2gO1VPVXr9wY7BpGkvR1DVLedDy3AV5bSPEzG7S8qX5uU4lTUBmEiJwAjARGqup44BssOC9AH6C2iFRV1UJgooi8ArwJvKyqWyLq6g1sAD6NcqhpwF2q+jD2sH7mmRTKQ97nVPVeVb0HGBDjD9/MO5cQ67116ZY1VNdRwLeqGi0uZZmvrV9Zw9eLSEOgM6ZU/V6vMl/XgOQNryulz21AsqblmQ1Q3lBdKX1uU4Ebg8ogVPUT4BMReUVErgIuAq4VkY7Yw75UVXeKyLVAdVW92HuIxonIFFX9EEBEegBXY135PeYRqOqKsN9TvQe6NbHne6VMXlX9PqzKUcDRwKsRh1oB1Atbru+tS6usYdX9Cbu+0Y5T5mubgKxvAvVF5J/e+h+wt2i/16vM1zUgeYH0PLdByJquZzYoecNI6XObCspdQzpARLqJSLgJZj5mpmmpqrep6uPAYnb/CVoDy8BMAZhZqqZXV2/sIbwcUBEZ4K2vIyJNvd+3iEhj73djbOB0ebrlFZF6IpIfVk8nYE6kvMAHwGFh8tYEpqZT1rD62gPrVHVV2LpArm0SsjbCcrXdCzwD1AC+wkw/W0SkhVeuH3YNA7uuAcub8uc2KFnT8cwGKW9YfSl7blOJm6ibAYhIB+AhYCJQDegKXAtchj2U07DZ2PeqaqjheRIbCK+FvZ1dgw04TwUmA7uAKsAXqjpYRP4C9FDVK0XkfMxGPg3oBrypqh+Ug7zVgNc9mQVz9LhaVTdHyJsF3AdsBtoAL6h/L75AZFXPo0xEHsU8sqaEHSOQa5uErD2AfwGjsbf111T1Z6+ug7xrvBAbpL/Be9MO5LoGKa+ItCLFz22AstYkxc9skPKG1Zey5zaVOAXlcDgcjozEmfgcDofDkZE4BeVwOByOjMQpKIfD4XBkJE5BORwOhyMjcQrK4XA4HBmJU1AOh8PhyEicgnI4HA5HRuIUlMPhAEBEWovIKBGZLiJTReS6KGWqiMiPIjIybN0CEflZRCaJyPiw9SeKyEwRmSMit0TUE3VbMvskWsZRcXATdR0VHhHpAzyIhWf5FGgB7AP8QVW3l6ds4YjI9ar6WCllxgOHakQupBhlj8aiC3QGbgF+peR1qOUVvU9V1/qoryUWSmeiiNQDJgCnq+q0sDJ/A3oB9VU1FI5oAdArIoxOFSzFxHFYSJ4fgN+r6rRY27Bo3QntEyFbqWUcFQvXg3JUeLzgnYVYpObBqnolljLjxHIVbE+u91Gmtx/lBKCqXwL/A6ap6vAo1+Fmb/lLsYjXpdW3TFUner83ANOxUD4AeCGJTgH+7UO8PsAcVZ3nvSS8CZxWyrZk9vF7TEcFxEUzd1Q6vDho+wArRaQulj7gK6AL8Iaqfi4ieVistGewXDs5WNDPLljE6J3Ay8A9QE+8TKRAT1U92zvOEOw/tAvYoKoPhtX7KBYzrSswAFOWDUVkMDAD692EH/9A4G7gCRHJVdUFXlDSPY6bCKr6kYjcARwLfOx3PxHJBg7GMrqGeAy4iZJRugEU+FREFEtF8Tym2H4JK7MYONT7HWtbMvuE46eMowLhFJSjMtHHG3c4GEsfPgF7xh/1lFJj4BPgc1UtEJFLgB9U9U4R6YWZh37C0mNnYYE3twB9VfVkABH5o/d9grf+eG+5UEQ+Dat3kqo+ICJPAcep6tsi8qCqDg4JG3l8VR3vmdDC6y9x3CRZiAUsRUQ+x0ygkdyqqsO8MnWB94DrVXW9t24AsEJVJ4ilBw+nn6ouFZFmWB6hGVgg1UhC4wmxtiWzTzh+yjgqEE5BOSoT36vq/QAi8g8sGvQ/gVwROQzYATSN2Gc6gFoyOLzB/z9g6QpeAU7AS6fglXvR+9kTSxYXGoj/JaLuWd73SvbsccQ8fhg9Yxw3GdqyO5fRsfEKikg1TDm9rqr/DdvUDzhVRE7GUkfUF5HXVPVCVV3q1b1CRN7HTG1jsNQlIVph2YnBejbRtsVaH2+fcPyUcVQg3BiUo7KyDNgXS0+wr6reBTwSpVzkG/aTWF6iA1R1Mtaj6hDaKCKXikh1b/0KVb3fU4oF2CB/rHoBdolxcCnliHPchBCR4zCF8rmPsgK8CExX1X+Fb1PVf6hqK1XNBs4HvlTVC8VyCtXz9q8DHI+lKvkB6CQi7Ty5zweGe9XF2pbMPuH4KeOoQLgelKPC45nnjgSqi8it2HPdA8gHioCzReQhYDXQQETOwlJwtwWuFpEHVHUlgKrOEpEleOM1qvqpiPQVkfuArcBv3gD8pyLSx1u/E1MCt3gKoS1wqYgM9eTqISIfYEnsHgaqiEiTyOOLJahrC1ypqrfEOG40NMp1uA2ojeUSOlpVd/q4lP2w8befRWSSt+6fWjKbcCTNgfdNt1EVG+P72JPnasykWgV4SVWnetd0Z6xtSe7zIXCZZ2aMWsZRMXFu5g6Hh4jUUNVt3rjRNWoZdTMSETnWG1e7HuigqteUt0wOR9C4HpTDsZtHRGQt8HUmKyePv4jIMVj205vLWxiHIxW4HpTD4XA4MhLnJOFwOByOjMQpKIfD4XBkJE5BORwOhyMjcQrK4XA4HBmJU1AOh8PhyEicgnI4HA5HRuIUlMPhcDgyEqegHA6Hw5GROAXlcDgcjozEKSiHw+FwZCROQTkcDocjI3HBYkuhSZMmmp2dvcf6devW0aBBg7j7TpgwAYCcnJxSj+OnvkTLBl1nIueTyPHL89wryj2qbOfjt6x75vaOZ27ChAmrVDUymSioqvvE+eTk5Gg0Lr/88qjrw8Hy9JRazm99iZYNus5EzieR46fq3EfNWK5FRUV7bCsqKtJRM5arasW5R37LVZTz8Vu2oj1zfqgo9yid5wOM1yjtrzPxOSolq6o1Y1DBDwwZOQ0Ni9ivqgwZOY1BBT9QOHNFOUrocDhKwykoR4WncOaKEkoIYJ8dKxh0eFsKxiwoVlIh5VQwZgF5/bLp33lPi4LD4cgc3BhUkgwcOLDc6vNbNhV1JkI65CycuYJBBT+Q1y+b/AHd8FKPM3DAQH70yhSMWcDoka9TNUuZ1eAI8vbfTH67CciMnzj3QNi0HVg6FmruQ+GSKvTv3g7JKvnupqp0OKL87lF53p9Eyu4Nz1wQdQZ97ETKpuJ8avrUJIleI5ewsBQ6deqkRx11VPHywIEDfV/kUGNZWa5xJp6PqjJk+E8UjF1CXpeN5Lf5Fn6bypC5B1Gw7ijyag5DgaFbTwMgr9Yw8uu8gHcqJSjcfgiD1g0hr86H5Lf5BmnYARp2RJsexJCpbSiYsI6heb3J7dIsvSeZAJl4j8pCZTsfqMDnVLQT1syGVT/Dysmwagqsnc2mZdN4dzJc8h//5zNixAhGjBhRvPzCCy/MUdVOkeWcgiqFXr166fjx48tbDEc4m1fCoi9h6bewdAy6fBJDNl5KwZbTuKjGcNZvr8YwPYm87MXc3msbd01rS8E0e8XLy6lP/nH7ISgU7YKiHbBjE2z9Dd38G0O+g4LZjchrMZn8+q/CurkMWX8JBVtOM8XVaSrSJhdaHwUt+0LVGuV7LRyOVLFzKyz7DpZ8A4u/tv/bjo22TapAo87QuAs0aAf7HQGdzkz6UCIyQVV7Ra53Jj5H5qNFsHwCzPsQ5n8Iv/4AKFStDS0PRQ69hfwWfWFSYwrGnwpAXr9sbj/lJO76YDoF0xZQb5m9ZBRM6AU1G5cwB4YQIL+bwshpFIwBOrwDTYoo+HaRmQVbV0GWbYTv7oKxd0LVmrDfkdDxdOhwKtTbL62XxeEInI1LYd4HMG8kLPwcdm4GBJocAN0uhpaHQtOe0Hh/e/5TjFNQjsxEFVb8CDPehJlvwYZFgNgf5PDBkH0iNDsYqlTbXX76NGC1twxDRk5j6LcLvfGpkwGKnSSA6EpKhPwB3QCKy+0e3zrHCm1dC4u/gl9GwfwP4Iur7NOiN3Q6C7r+Aeq1StmlcTgCZdNy+4/NeAOWjbN19dpA90HQ7kTrHdVsVC6iOQXlyCw2LoUpL8G0V2HNLMiqCtknQL+7oN3JULvJHrtEeufBbuUy6PC2JRRRuPLp37lpcuNJNRtCx1Ptk/svWD0D5vwP5rwPX98CX/8D2hwD3S+GjmdA9bpJXAiHI4Xs3Aqz3rX/2aLPzUrR9EA44l7oMBD26U7Ugdo0U6qCEpFTVXV4OoRx7KUU7YKFn8Hk52DuCNBd0DoXet1gdu1a+8TcNVI5hRQQmBISoveQYimn8PqKzYJjbFu0HhcisE9X+xz6D1gzx/70016Bjy6GanWh+yVw0F+sTBQKZ66gf+eme9StqoyetTKjnTIcFYw1s+Gn52BqAWxdDfWzoc8tsP8F0KR7eUu3B6U6SYjI90Ah8IqqTkmHUJlEWbz4HKWwfQP8/CJMfBzWL4BaTeGAPOhxOTTq6KuKWG7m4YrGr+ddLGUXuW4PJRW1siJYMgYmPw+z3oZd261XddDV9oaaVSVw+R2OqGiRjd9OfNx6S1lVbdz0wD+bs0859JQC8+ITka7AQmAQ0B34WFVHxN2pEuG8+FLAxqXw45Pw07Owba3ZuA++xv40VaonXF1QPZCUKYvNK+Dnf9v5bvgFGnaAXjdC90vQKjWCU4oORzi7tsP0N2D8Q/DbNKjbCg78ExzwR6jbsrylK0EsLz4/CuoAYDZwIXA1sASYD3ytqm+nQNaMwimoAFkzB8bdC9NfMzNep7Og19/N8SFDSKm5rWgnzBkGPzxgnoh1WkDO39AeVzDks8V7jKE55eRIiu0bzIw38THYuMS87nrfBJ3P3e1UlGGURUFNAuoBI4H/U9XZ3vohqpqfAlkzCqegAmDtXPjubhubqVLd3uByrreexN6IqnkAjrvPTC41GqIHX8+QladS8N1SgFKVkxu3cuzBjk3w41Pww4Ow9Tcz3/W5GdoenxEOD/EoyzyomcBlqrohrLLqQPn4He6W4VFgM7AROBC4XlV/FZFBwDVASN42wBBVHRqx/wnAJcBkoC9wk6rOSo/0ewnr5ptimvqyvbkdcq29ydVpUd6SlS8i0OZo+/w63nqVYwdbo8IJpe7uxq0cJdixBSY/C9/fb+bk7BNtKkYGWSaSxY+CGh5STiJyEHCxqv4NUwLlySZVvQ1ARG4GbsVkmgmcoaqLvG3DgHej7P8YcIGq/igifwVuBC5Ph+CVnk3L4ds7YMqLNuP8oL+Yp1CG2b0zgha90FPfY8jbhRT8uJm8WsOgam0KxhwHRbvIP7XHHr2k/p2bktcvu8R8LsAFwt3b2LUdJr8A399r47ptjobDh8B+/cpbssDwo6CK7TCqOklELkyhPL4JKSePLKwnhaqODa30FOosVd0YpYrlQOhf3BSYEO04EyZMiGpmcSGiorBzq3kKjbsHdm6Bnn+CPv9wERbiUNzr+XGz9YgOvgTG3A4zNlMw9jRYO5v8i05HPK8/iD6ZOPTbjVvtBajadIyvbjC38f2OgJNft6kZGU6iz2VMBSUi1wHXAw09s5kAO4EPkhcvMUTkE6B5lE35oblZItIQOB44K0q5a4AhMaq/DnjKM/XlAC+WWeC9FVWb9PfVTeYu3n4g9H/I4nQ54jJ61so9Fct5X5K/4HN4ZzQF0/vS/9kLyB1wJbTZPd0hfsQLp5wqLSsnQ+HfYNEXFm7ojA+g3UkZP8aULDEVlKo+DjwuIueo6jtplClchrgGeRFpADwNXKqqqyO2NQdqqOrCKPtVxxRtX1VdLCIDgf8AfSLL5uTk4Jwk4vDreCj8qwWUbNIDzv4M2h5b3lJVGHK7NGNoXu89HB4k+1jy/34U/b98j9x538E7b5viP/KBmBN+HZWYzSusZ/3zv6FGQzjqCTjwyoz1yotFLMtTrJeqUk18kcpJRM5T1beSES5IRKQJNo50o6ouEZGzVPW9sCJXAc9G7NMaWArUABoDK71Ny4DURz6sTGxYAt/80yIm1G4Gxz0PB1xaPAHV4Z9YzgySVYXcY8+FnafCxCfMdPpyD+h5BXrYHQz5cpX/iBeOisnObZ7Z/G4zmx98DfTNh1qNy1uytBDPxPeGql4gIvOxnPNgZr76QLkrKOBTTP7XvT/jBuA9ABGpAfRS1Tsi9nkL89b7xnOMeE5EZmJegFelTfKKzI7N8MND5sqqO6H3zXDoP6FG/fKWrPJStSb0ucmibIwdgv70LEPGZlGw6STyDmtN/qn+AuE6KhCqMPu/8NWN5g3bfgD0f3ivM5vH60EN9r7/papPhlaKyJ9TKpFPVPWQONu2AadEWX942O/ngOdSI10lRItg+usWCHXjEuh8jpmbGrQrb8n2Hmo3hWOeZHSDPAreXWbJF5d/isx8ALqcG0wgXEf5s3yCjTMt/srSXJz1KWQfV95SlQt+JurWAbJUdYOI7KOqv6VHtMzATdTFYsoV/tWiHzTvBbmPQqsjyluqvZrCmSvoX2My8tUNsPInaHkYHPUo2qKPm6hbUdm4FL651eYN1mpiEfx7/NFi51VyyhJJ4n3gVVX9r4icCXRT1btTJGfGsVcHi123AL662YKd1t0Pfnef5TqSrPKWzBGiaBdMHWoN2+blFpX6d/dB/TblLZnDLzu2wIRHbKLtru1wyPXQ91ao0aC8JUsZQQaLvUVV7w9bvlVV7wlS2Exmr+xBbVsP398HEx41ZdT7Juh9I1SrU96SOWKxfQN8/4A1dAA5f7fJ0S4XVeaiagk5v77FEnJ2OhOOfHCvDAFWllBHkcl4YifncVRsinZZssAxt5lba7eLLIGZyw6b+VSvB0fcDT2vsHHCcfdYJI9+d1tmVOddmVksGwej/grLxlpm6JNegdb9y1uqjMOPgpolIj8Dc4F2wJPxCouIH//HIlVd66OcI10s+tLGmVZOhn37wRkjLYW5o2JRvw2c8rrFPRz1V/j0MkttkvsvC4XjKF/W/wLf/MMcjuq0gBNegm4XuxeIGJRq4gMQkf2BA4CfVXVmKWW3YnON4vm4VlHVCmEkr/QmvtUzLQLE3OGWXfPIB8xDz7koV3xUYebb8PXNsH4hdDgVjnwIGncub8n2PrZvtDQr4x+2+9LrBos0Xr1eeUuWESRt4hObTHEENleogYjM0vhabbqqHlxKnT+WdlxHitm8CsbeaVGQq9aCI+6zFBhV3XzlSoMI7H8edDwNJjxmQUVf7m7Be/eiyZ7lSsiJZcxtsOlX2P/3nhNL2/KWrELgx0niMaA6MAfoBGxT1evjlK+pqltLqbPUMplCpetB7dxqJp9x99hbXc8rLDR/beeWXOnZtBy+zffC5TSAw+6wtN9JZDF2+GDh5zD6757Z/HAzs1aCFBipoCxOEivDvfZEZHC8wn4UT0VRTgDr1q3jiiuuKF6usG7mxeaeW7yArqeYx9A+3cpbMke6qNMcjnsODrraGs5R18Okp+DIh6HDQGfWDYrfplsEiHkf2ET2AW9D57Pd9Q0j0s0ciOpT76cH9aCq3hS2/C8vH1SpiMjNqvqAn7KZSqXoQS360uLmLRtn6Z/7P+ICuu7tqML8j0xRrZ5hDhT9H4FmB5W3ZBWX9Ytg7BAz6VWrA31vt9h5VWuUt2QZT1l6ULO9tO8LKMWLT0TeDl8EDgIqtIKq0CwbZxM4F30BdVvB8S9C90ucx5DD3ubbnwxtj4PJz1mCyVcPsXh//e52ySUTYfMKy4r80zO2fPA1Fp+ytksaWVb8RDN/QUS+AbpTuhffelW9LLQgIs8EIKMjUVb+bIOyc4dDraZw1GOWPNA5QDgiqVINDr7aIoR8d7eNT858yxwpet3gGtl4bFtnXnkTHrWx3QPyrNfkongEhi838xI7iFyhqs/H2NZOVeeLSBNVXSUijSPzNFU0KpSJb8UkGHcfzHrHoov3vgkOvtZFE3D4Z80cc6SY8aZ5dx74Z4siUida3tC9lM2r4MfHTZlvWwedz4V+Q/a6SONBknAsPhFZDayNXA3UV9W40SREZLiqnpqkrBlFhVBQi78xF+L5H0H1+vZG3OsGqNmovCVzVFR+m2GenjPegCo1rAfe+0aou295SxYoubm5ABQWFpZeeOMy6zH99KzlZup0Jhx6KzSPO6vG4YNkxqCuVtU3olR0gZ/jJSJcJpOxXnxFu2D+h5abacnXZso74l446KpKHWTSkSb22R9OftVMVt/fa72FSU/B/ufDIX8ttVEunLlijyzBYBlVK1y09d+mWcLIqUOhaKfNZTr0H84DtgwE5sUHICI9gabATGBJKRN1XQ8qlWxda/HyJj0F6+aZ80PvG6HHZVCtdnlL56isrJ1nmV2nvAQ7NkLrXFNU7U/Zw+mmcOYKBhX8QF6/7BKJE1W1OKni0LzeGaGkYvagtMjcxCc+AYs+t/Hbbhdbgs6G7dMuZ2WnLJEkbgROBhYBQ4FrgJtL2y0JGR2xUIXl4+HnF2Haq7BzM+x3hM1I73iGDXQ7HKmkYXs4+nE4/E6b6PvjEzDsNKjXBg641BwEPOeA/p2bktcvu0R2X9id8TevXzb9O2eo88XGpTD1FZjyb1g719LMHHEv9Lgcajcpb+n2Ovy4mddV1aO8OU2jROSo0nfhH2UVzIH9Waa9ZqaF1dPtLW7/C2yipbN7O8qDmg2h9w0WFmv2+/DzCzB2sIXNyj4euuchHQaUyO4bIqScMi0dfbWsIpj1LkwpgAUfW+9pvyOg3z02zuReAMsNPwoq1H/XiOWYqOqUpCXyiYg8CmwGNmJxAq9X1V9FpDrwPDZvqzmwVFXvirJ/NnA7FsIpG/i7qm5Mtdylsv4XmPO+fRZ/ZX+WffvBcS9Al3Pc+JIjM8iqas9jl3MsseWUAjP/fXA+VK2NtD+F/E7nQFGnYiUVtHIq0zjXjs2w4GNuPXg6hzX/DUZ8bb2lPrdYepJGe+TOc5QDfhTULhH5GKgtIn2AiX4qFpHRwEBVXS8iVwI1gadVdXvy4pZgk6re5h3rZuBWzPx4BtBIVQeJSBbwm4i8pKpLIvZ/FshX1e9FJGS2vD0g2fyzc5vlhFn0Bcz/2Ex5APt0Nw+hrhe66NOOzKZBNvS7Ew7Lt5eqWe/A7Pdg5juw+UpggJXb8puZqwNQUAmPc6laxIyFn9uY0sLPYOcW+jSryuilTTnlby9Dm2PcJPYMI56b+VOq+hfv9/FAT+AnVf3MV8UiP6nqgSKSg/VoRgLZqnpJMKKXONY/MPf3f4hILnCPqvYTkYbAGKCXqm4JK18N63nVVFUVkUOAf6vqIVHqjnqBEp0/5u1kZrvl42H5BFj2HSz5xlxWJQta9IGOp9u4klNKjgqM7trJkLe+oGDyTvIafQVb11Cw5TTy6n9OfvdFyL59LXBq85ykMjWHK6KQkoKwca7D9iP/kPXI8vHw6w/m6brRe0dt2AHangCdz+KYS+5kl4o/N3NHmYnTe07YSSLbiwSxABiqqp8mKMsOEakKXAw8oKpvi0hC7nAi8glmposkX1WHe2UaAscDZwGoaqGITBSRV4BmwMvhysmjCbAlzBtxvVc2WOaOhFWTYc1sWDsH1syysCgAUgWadLfB1zbHWDZNZ75zVAJUlSEfzjTl1C+b/AEnW5y6/46jYOaxMK2Q/Nk37e5I1WsDjfe3ia4N2kGdllCnBYXL69O/czOkWi2bNJxVDYq2ozu2MHr2KvIPE9hQ1UyIKybBtnUULGpLXoNC8uf8C5lTtLv+/Y6ANsdC22PsGB67dEi6L48jAeIpqAtVdY03VvNnEWkCjFDVT3zW/QTwE2bau8Vbl1BIA1U9Id52EWkAPA1cGopYISLXAtVV9WLPxDdORKao6odhu64CaomIeEqqPrAi2jFycnJI2s38+/th6Rib3NiwI7QfCE0PhBa97Nu5hTsqIaNnrdzTIaJBW/IHtYGR0ygYk0v/C64gt85MWPEjrJlp5rcpBebCDhRuP4RB64aQV+tZ8uu8UKzMVGHIpssp2HIaQxvkk19tItS6nILZpwENyWs8lvyui5F9bofmvSwrtIuCkTHEsdhFXR9TQanqGu97ATBYROoBr4rIM6pa6kQAVX1FRP4L7FLVLSLSERhb6hn4xFOYjwE3quoSETlLVd8DWgPLPBmKRGQ5piQRkdaY08QOERkF9Aa+B/oBHwQlWzED3rRoDkmYMByOikpul2YMzeu9hwODiJA/oBv9Ozf1xoY6QYcBu3dUhW1rYdMy+m9cRl7hOgqmnwat+pN/wGIo2sGQKW0oWLUPeZ030L/vDVCnGYyrBhM3WB1dfw8DurnUFpWEeGNQzVR1hYh0B64ALgSmAM9FizARo45qqrojMGlL1j0RU7ChWH8bVHWgiLTAIq5PAWphvaNrVHWXiHwL3KSq33g9w3xgHtAG+Fs0L76Mm6jrcOwlhI8z1Vtm/8ENLXtFHXOKtt2Pt2BCoY4cKSOZibqviUhtoAvwMnB4KZHMIw/4b+BMEdkELAUmA5NVNWa6jkSI5tDgrf8VOCfGtsPDfi8ALg1CFofDETyhHhdAwRhbF9UhIjTOFbYOyLj5Vo7EiaegWgD3Ae8l6Rr+O6C5Z07bD5ur1DOJehwOh6MEUce5oMQE4d2mREdFJZ6C+rOqjilD3d8BjYAV3hykJcCH8XfJPDI2WKzDUcmJdCUHSvSO/I1zOTKRQIPFJoOI9AP+DbwIjMPMe+tScrAU4sagHI70U+o8p4CiUrgxqMygLCnfk+U14BXvGFcBPUWkpqp2SOExHQ5HJcCZ8ByQhIISkZaqusxH0cWqekfEvjUSPZ7D4dj78O+q7qjMxFRQInJxjE0DieElF8EkEblOVR8PrVDVbQnK53A49lJiKSARccppLyErzrZLgHZRPo191t0cuFJElorISBG5R0T8KDaHw+FIGYUzV0SNaKCqFM6MGlDGUU7EM/Fdq6pTI1eKyAF+KlbVc73yNYDuQA/gUOCdJOR0OByOMhMZBT1EJmb7dcQPdVSsnESkK5byHeAi4PLSKhaRu1X1Ns+sNxGYKCIVLpa9czN3OCoPkdl+Q/2oCpHttxIRmJu5iDyIRZPYF5gFdI0VxSFivwLgs1BYJBFpCrytqn4y8mYMzs3c4ahclBZCyUWfSD9lcTPfoqqneSnfHxCRv/s85p+AT0RkDvaiUoAlBXQ4HI5yI14IJaecMgs/Cqq6993Iy++UE6+wl4dpIvAj8BfgDWAncLqqzimDrA6Hw+HYi4jnxRdih4gMBMYDG4DS4vK97NV7KaacsoE1wIUicnbyojocDkfZiTTx1Vs2noIxCxgyclpymbIdKaPUHpSq5od+i8hnlKKgVPUL4IuwfaoC3bBgsX2Bd5MV1uFwOMpCZAilLx95CICzzz7bRUHPQOJN1G2iqqtE5MiITXG9+ERkYrgTharuxEu1AbwarUwm47z4HI7KQ2QIpVGP2HoXQim9lNmLT0Q+UtWTRORHbDwp9ErRI5q3Rdh+W4DZcWQToIGqtolTJmNwXnwOR+WicOaK4hBK4cFiVZXRs1Y65VQOJOzFp6oneT+vVdWvwyrqV8qx9vchzy4fZRwOhyNwXAilioOfMaivI1b1AGLmiVLVhWUVyuFwOByOeGNQq4G1QDUs/t5KoBmwDng2HcLFQ0QeBTYDGzEHjOtV9VcRqQ48DyzA4gEuVdW7IvYVLBXILMzjsAOWoHFT+s7A4XA4HPGI52Z+taq2Bx4EmqlqNqagnkqHYD7YpKq3qup92BjZrd76M4BGqjoYm4f1Ny/lfDhZwDxVvUtV7wQ2AVemSW6Hw+Fw+CDeGNQb3s/moZ6Fqm4UkfppkawUVPW2sMUsrCcFsBxo4v2uDywFVkfsuwsIz1UVvn8JJkyYENXl1M2XcDgcjsRI1H3fTySJriJyI+aZ1xnolIRcSSEin2BmukjyVXW4V6YhcDxwFoCqForIRC+iRTPgZVXdEucY2UB74NpgpXc4HA5HWfCjoC4HbgFygWnAZakUKBxVPSHedhFpADwNXKqqq7111wLVVfViEckCxonIFFX9MMr+rYD7gPNiJVPMycnBuZk7HA5H2YkzrSnq+lJDHanqalW9SVVPUdUbgbZlkjAgRKQJNh52o6rOF5GzvE2tgWUAqlqEmfxqevu0DqX8EJEOmHL6k6quDtvf4XA4HBlAqT0or5dxDbvzQfUEYk7UTSOfYvK/7mnfDcB7wCPAkyJyB1ALWAQM8/Z5C7hJRMYDXwFLgOHe/rO9/R0Oh8ORAfgx8d2PNdwnYnH0MmKSbaxQSar6KxA1tbyqHh62GOnZ53A4HI4Mwk8080mq+j7mlv0Z1iNxOBwOhyOl+OlB5YhIW6CpiFwIHAXcVco+lQYXLNbhcDiCJciU74dic4S2Ag8Br3o9qr0CFyzW4ai8hAeLdZQfSad8V9VxYYtnish5gUrmcDgcDkcUYo5BiUh1ETlPRI4KW3cGMDgdgjkcDodj7yZeD+o1oDrQWETaYLHq1gKDUi+Ww+FwOPZ24imopap6PYCIjAFuVtWv0iKVw+FwOPZ64rmZhwdY/SiknETkmNSKVHkQkYSDI2Yyle18oPKdkzufzKeynVMqzydeD+oEEanr/e4TFsW8L/BFSqTJQJybucPhcASLXzdzVDXqBxiFpaSI/HwRa5/K+MnJydFoDB8+POr6cAC1S1w6fupLtGzQdSZyPokcvzzPvaLco8p2Pn7Lumdu73jmgPEapf2NZ+K7SVXvjPwANyWjMSsbEdo/rfX5LZuKOhOhPOUsz/NJpOzeej6Jlg26zvI+96CPnUjZinI+EGcMSlV/iLF+QoIyORwOh8ORMH5i8TkcDofDkXZKDXW0tyMiK4GFUTY1ANYFeKhE6vNbNhV1JkJ5ylme55NI2b31fBIt6xf3zKW/XCLEqrOtqjaNXOkUlMPhcDgyEmficzgcDkdG4hSUw+FwODISp6AcDofDkZH4SVjoSAARyQIux5I6Hq2qU2KUWwAs8BaXqOof0iJggiRwPhcCBwO7gLmq+lz6pEwMEWkM3A/MAzoB/1TV5VHKLSCD75GIHAucCazAJkreGbG9JvAwsAQ7z/tVdVbaBfWJj/MZhAWt3uqtelFVX02rkAkgIi2Au4EDVbV3lO0V7f6Udj6DCPj+OAUVPAcC44DNpZQbqqqDUy9OmSn1fESkFXADcLCqqoj8ICJfqursdAmZIPcCn6vq2yIyEGskLopSLmPvkYjUBp4FuqvqNhF5T0SOUdXwMGTXA4tU9UER6QG8CPyuHMQtFZ/nA3C+qi5Iv4RJcQQwDDgoxvbrqSD3x6O084GA748z8QWMqv6oqpN8FD1SRG4SkbtE5PBUy5UsPs/nBGCC7nYJHQuclFLBysYpmIwAY7zlaGTyPToMWKiq27zlaOdRfJ6q+jNwYFhMzUzDz/kAXC0iN4hIvtcTzlhU9V1gQ5wiFen++DkfCPj+uB5UEojIJ0DzKJvyVXW4z2puUdXvvTfHiSIyQFXnBCelfwI4n2aUfHDXe+vKjXjnREl51wONRKSqqu6MKJsx9ygKfq55rDLrUytaUvg5n9HAB6q6UkROBt4BKnJ2hYp0f/wQ+P1xCioJVPWEAOr43vveLCKTgH5AuTR+AZzPCqBj2HJ9yulcQsQ7JxFZAdTDEnDWB9ZEUU4ZdY+iEDqHEPW9dYmWyRRKlVVV54ctfgkMF5EqqrorDfKlgop0f0olFffHmfjShIhkeZmJEZFjROTEsM0dgbnlI1lyhJ8P8AmQI7uTwhwGfFQ+kvniA0xGMKXzAVS4ezQWaCsiNbzlfsAHItI4zExUfJ7eGMdPqpqpb+elno+I3CcioZfqTsD8iqacKvD9iUqq74+LJBEwItII+Avwd+BV4A1V/U5EDgJeVdUe3sM4GJgA7It5iN1XTiLHxc/5eOUuBHphXnyzKoAX3wNYCKsOmClveUW7RyJyHHA2sBLYoap3isiDwGpVvV9EamEOIMswBXtvhnuJlXY+1wEHAPOBHsDjqvpd+UkcHxHpD1wMnAg8AzwC3EnFvT+lnU/g98cpKIfD4XBkJM7E53A4HI6MxCkoh8PhcGQkTkE5HA6HIyNxCsrhcDgcGYlTUA6Hw+HISJyCcjgiEJE+IlIoIt+KyGAReUpEnvTmSV3lBZFNpL7TRSQ7NdIWH6OhF6wzXpljROThBOq8QkTWisjTIjJERN4VkYEi0tqLlTe4rHI7HPFwbuYORxS8xreuqt7gLY8GHlDVD0VkgapmJ1DXUCzwbGEKRA0dI9s7Rm4p5UQT+NN7yniAqk4RkXbAOFVt5inD7EwNpuuoHLhQRw5HKXiz45sAq7xVWSJyIzYxeZOqXuqVOw04DvgFmwD8d6ArFv15kIj09SY0Rit3ADb58Wvsf3kQ8BBwlHeca1R1vDdr/zFgNtAKGK6qnwBXANmeYv0YOAP4PfACcCgWoqkIOATIFZEq2ETLVUBNoLGqXlXKpWgZdg0AOnoTaQ8H7lPVD0SkPfAv4FtssuYjqjpJRPoBlwIzgRwsLcMG71xWAA2ASZmcPsNRDqiq+7iP+0R8sCgS33rfTwIXhG3bAtTzfk8B9gEaAYuB6t76m7HIAABDgVzvd7xyg4G7vN/XY1EtwJTNk97v+7D8VQC1sGgYVYFsoDDiHLZ6x6uCKbziMsCfgKfDyv4xxnVYADyN5QN7EujirR8EvO797gWM8H63Ag7xfh8CvOP9fgy4HRtW6AnU8WR43tsuwAygWXnfe/fJnI/rQTkcsflWo5uwVqhqKAr1KizgZ1Ms5Mt2b/0crMcQScdSyoXi/a0N+72G3UFFewK/icgt3vLPQKy0BstVdY33e1LEOFhPwgLfquqLMeoAU2TRElWG9g9dA4AdwPkichIW/LSpt/4e4Fbgeyzu3k2eDC3DzmUK0IIKHDDVESxOQTkciRNtDGcO0FhEqnvKpxMwydu2Cxv+6Qj8GqecH34CflXVJ7BKLwJ+A2pgvRBE5CC1HF7xxpp+wpJR4u1zpao+m4AcxKj/Fiw6/D0i0hkzLwL0VdXrvYDCb2L5wn4Ctqnq/Z4MZ7A7g7HD4RSUwxGJiPQCjgSqi8hZqvpe2LbLgAYichawGmgLXKqq+SLyF+AJEVmMmdP+6u32OXAZZt76IxZ8t0Q5rzE/EughIt8CA7E8VZ2xbL89RaQPZuJ7UERuw8Zt5qnqLhFZBmwVkX8B07xzaCAif1PVf3lyXIVFDD8Jy976iDdmVRXrvUReh0u9Y1whIk+olwtLRJqHydcRuNCr9xjgPeA+Lyp59bD12SLyKNY72gIUAuu8cxnsld2iqu/7vlGOSo/z4nM4HA5HRuLmQTkcDocjI3EKyuFwOBwZiVNQDofD4chInIJyOBwOR0biFJTD4XA4MhKnoBwOh8ORkTgF5XA4HI6M5P8B1piDzz5REOYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Plot the fitted curve for raw data and photometric phase plot\n",
    "\"\"\"\n",
    "plotting(fitx1, fitx2, fity1, fity2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The reduced chi squared value for the original fit is: 1.5433422548104576\n",
      "The reduced chi squared value for the photophase fit is: 4.4633400448587635\n"
     ]
    }
   ],
   "source": [
    "chi1, rchi1 = chisquared(fit1, time_wasp)\n",
    "chi2, rchi2 = chisquared(fit2, photophase)\n",
    "print(f\"The reduced chi squared value for the original fit is: {rchi1}\")\n",
    "print(f\"The reduced chi squared value for the photophase fit is: {rchi2}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEdCAYAAABZtfMGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6qUlEQVR4nO3de5hcRZ3G8e/LXSSJRBJEkASRZAFB3AQEwppBVASCFxbFVdBEAUEusi4giywgKiKCgigoIEFAFhFvJKjIJRMgopBoFI0mIISLsiQIJCDIzXf/qOpw0unu6enu6enO/D7P08/0OV1d9evTZ071qVOnSrYJIYQQOs0agx1ACCGEUElUUCGEEDpSVFAhhBA6UlRQIYQQOlJUUCGEEDpSVFAhhBA6UlRQDZL0U0mWdKekHfK610l6QdLwvLyNpPmSFknapfDeiyX9pEKeknSqpDskzZL0C0nnSxqVX/+OpCck3S+pV9JvJP1M0tgKeb1c0imSZku6Jaf/1358vpsk/UPSwvze+ZLmSnpTA9tqzZyHS7FK2k3Sz+t470657MVVXp+T850vaXxh/SclPSJpVo28f5A/Y09/P1NZPmMlLWgmjwp5HirptsJ+cJWkt7WyjAbj2k7SL/M2/1X+XkuPP+U0H5B0cT/yPE7SYknLcj6/kPRnSf/dgnh/Uun7rbRPNlFGzX00NMF2PBp8AA8BBxSWPwoY2Luw7mNAT2F5LeBu4DlgZFl+HwZuA9bJy8OAu4CJhTS9wOfycwFXAHMqxNYD3A6sm5cPAf4KvKwfn28xcHBh+UJgcRPby8DYQuwj6nxfT61ygfuAEyqs/3mdn7Gnnjj6yOcVhedTgd4m8voscBMwLC+vAZwGPDFQ+3I/4xubv8vXla3vzX/XBDYorL8UOLWPPE8Fbiss75rLaOq7AYYDqvH6in2yyXJq7qPxaOwRZ1DNmU3aMUsmAT8rW/cm4JeF5bcC5wNLgfeU5bcT8AvbzwHYfpJ0sHqsUuFO/xnfBXYtnbUVPAx81vazefl/gU2A19fxuaq5FhgjaaMm8gBS7LaXNZtP9j3ggOIKSROAX7co/z7ZfqIV+UjaEjgR+Hj+/rH9T+AU4C+tKGMAHQ9g+0XbTzWTke1fAH8DJjaZz/L8fxK6UFRQzekFJheWNwBmlK1b1/Y/Csv7kSqLqyk7qAL3A++TtG1phe2rbd9bI4a1Sb8CXyiutL3QdrEZcb3899EaefVlLeApYBmApE1zM9mtuant1GJiSdMk/V7SjZIOLqwfVWomKqx7T85nVm6SnNSPuK4GdpA0rrDuAFLljaQP5fJmS7qyQmVejLmU9hZJ35e0ceG1fy00Z/1C0sfz+ptKTUWSJgMn5Hh6JZ2XP78lzczpD5L0V0nXVAjhPcCDthcWV+aD7M75/ZcVmyYl/bis+bTUdHmcpBmSlks6WdJzSs3C45Wan/8g6U+SXiVpdH7fLfmzvbsf2x9JvbbvyE2AK5q7JH0CeAcwNW+Pj/Yj27VIrRSlMvaUdHvOZ4akV+f160u6On+/t0n6cl5/nKT/K+6XNfbJbcri3iVvm95Cmrr3UaUm2l/mfWOmpH/px+cOJYN9CtfND2AcqXIYTWr2+CKwLfA8qbJ6HfCZQvq1gavy851yuo0Kr48gnZX9E5gFHAqsX1ZmLy818a1DOmO7ro5YpxbTAZsBS4DNarxnMbmJj9Rscy7w3sLrWwF7FZZnAXvk59sCTwNb5OXDWbmJbyz5uJuXDyQ3eebXHii81kMfzSfAn4H/KSz/PP+dRKqUR+XlLwEXl33Gnvx8t7xNSmlPAm4qfDdLgcl5eXPgD4V8ip9tKoUmvvy9PwrsWlg3gwpNT8AFVGiyrfLd9FQqv/D6Jfn57sA+wA+AEwtpvgxsVdpewGn5+Saks5exVcoem8v7FWl/7C37vCt9XzTWxLcX8C1grby8BfAkMD4vHwHcmJ9/HLigsJ/eWals+t4ny+Mu/x7r2kdJ//t/46Xm9U8AU/v6TuOx6iPOoJpgexGpKW1yfswGFgBPkA52byb985a8jVShYPsO0kHk3wv5LbM9mXRAuQc4C/i9pE3Lij4o/7K7DXgAOKhWnJJeSfon+Xhh9dPAn/LfWo6TdBvwOKl5cEbhtQeBPfIv7l5ga2BCfm1/4Hbb9+Xl7/ZRzm+B6bmsS4HXSBrdx3uKVjTzKXXkuCOvnwrMsL00L18JfFCSKuTxYWBmIe104C2SNgemAE/ang1g+wHSD4g+2X4euAr4UI7vjcDvnI9eA+jHufxZtq8DLivEsCbpQH133r/eRqoQsP0wad/6jz7y/6DtHts9LYp3u3x2dC9pH7jEdqll4APAXL90Znklad/bhNQE/m+Sdrb9Iiu3YBT1d58sV+8++iKp4vuQpPWBr+d4Qz9FBdW80nWoSaRfgAZuIf2T7MzK15/2Aw4pNRORzoDeV56h7dm2DyH9Sn+c1NGi6PJ8YNjJ9qG2H5N0QqH56R2lhJLWBS4HjrJ9f6GMx2y/2XbF61sFX7K9G6nS3I3U2aLkBFIlvEc+SP0MWD+/tgmF5sQ6yrkWuNX2boUD3vo10pe7GthWqXn0fXkZ0pniWwrb/DzgEeCVFfLYjHSWVLK0sL78NWzP6Ud8l5Gab9cl/aC4vEq6e4DyHySNKr/Gdx3wylyB70k6a4L02QAuK2ynLYGX11tQiyqpu3I+rwcWAl8tvLYZsE0hvh+SmsQ3tn0VcAZwjqS7ST80KunvPlmurn3U9jO89KP1PuCbpM4aoZ/WGuwAVgO9wNGkf67led1s0q/Pe/LOiqR1SM11K9qt84FijqSNbT8i6WjgBtt/hHThXdJNpOalmmyfQfonXUHSGsC3gW/avk3SCFKzw5L+fkjb8yRdBBwt6Wu5It4JuKX0GUlNWSUPk5pAS7FUqhBKr5WaSH+Wl9eulrZGfL+WdA/wfuD1tn+XX3oQuNf2EYXyNrJd6Vrcg8CownLp+UMVXiudCf3WqRNDX/HdIekR0jWmcbb/VCXpD4AzJI3LZ+ilsoYB1wDvcrqm+Rywbn7tFX2Vn2N4XtJ3SQfwlwP/mV96MP/dv3T2KGk9Bun4YPtpSZ8E7pT0b7ZvzTHOtb1PKZ2kDYHlSp12vmv7ivyd3CjpT7ZnlWXd1z65YptmryikrXsfza89YvvA/D93KXA21SvOUEWcQTVvNrANqamtuG4nUhNeyduBW8veewfpmkepmW8k8LFS81M+KO1JOiNrxNeAO4GbJG1A6rq7d877lflC78h+5Pcl0j9p6SBxD7CjpDUkvZx0hlVyDbCLpNfm5Q/UyPdvpGbR0j1W76ietKbvkZoyi2etlwL75IMZSvdKzVj1rSvS7q2Xeil+GLg5N+fNBIZJenPO57Wk6x6VKqcnyb+slTpalA70l5Ou+1S9/ys3P50CnJ+/s9IB7yvA7/1Sh5v7eKlH5t7V8qvgMlIlvkbpDML2X3NMxabib5DOmlvhSWB9pXvzvlPPG2zPBW4kfZ+QOha9SdIYWFFhzCYdw47kpX3yLlKT35oVsu1rn7wf2Eipw8gapP+9kv7so5sCF+XPsQyYXyWe0JfBvgi2OjxIv8z2KSyLtEOXOgxMJlVWvwY2LaT7dE73f6TrGW8k/SP+gnRmNhc4rpD+O6R/kvuB6X3E9DZSO3j5Y2p+/TWkJquKnSRI9+H8g9TU8snC+m+Tzig+B7yK1DFifo5tVv6cH8hppwF/AG4Gjsnl/xIYk/86f841gXcD9wLXA58vpN0p5/8P4Ht9fOY35PdtW7b+wJzXzaRfwOPy+h/kfOcDE8rS3pJf37iQz4Qc7+yc1+sL26oU76bAhvm7ngOcW3j/a4BngdF17FOH5Px6SfeznUSqVEqv75a/mxtJFWmx/MsKn+sDFfJeQDoTK64bTTqA35ofJ1aJa7vCd/cr4KAKr6/0fQG7kK533lElnuPyfrOMwr1rpAryhbwNNiD9yJuTl28Gds7pds7fwc2k/5kvFPL9v5z3R/vYJzfPr58BLMrb4rOk/7fz8mvvpo59lHR2+i3SdbxbSJX/5oN9nOrGh/KGDyEMsNwUd4nt/QY7lhC6QTTxhTDAJJWaDf+D1JsvhFCHqKBCGHhbkpp79iI1G4YQ6hBNfCGEEDpSnEGFEELoSFFBhRBC6EhRQYUQQuhIUUGFEELoSFFBhRBC6EhRQYUQQuhIUUGFEELoSFFBhRBC6EhRQYUQQuhIUUGFEELoSFFBhRBC6EhRQYUQQuhIXV1BSXqrpPMlnSrplCpp3ifpz5KmtDu+EEIIjVur7ySdSdL6pGmpt7X9bJ5aew/bNxXSbEGaNfbBwYozhBBCY7q2giJNI32/7Wfz8hxgH9LUzwDYvg+4r9rZVYmkinOOvO51r2PEiBF9BrJs2bI+09WTZt68eQBMmDCh6bxaGVc96eqNvd1x1ZtmMLZ9vek6Nf7Yd5Ju3vbQ3vhLZVXwqO1Rq6wd7DnnG32QZif9UWH5YOCKKml7gSk18nKlxyGHHOJ61JOunjSlcluRVyvjqiddvbG3O6560wzGtq83XafGH/tO0s3b3m5v/NWOtcBcVzg2d/M1qCXAsMLy8Lyu3yZMmFCt4gohhNAipWPrIYccUtextpsrqNuBMZLWzcuTgOskjZQ0fBDjCiGE0AJdew3K9tOSDge+Kmkp8DvbN0k6E3gMOEOSgE8DY4ADJD1v+/p6y9h3331blq7evFpVXr3pWplXvdodV6du+3rTdWr8se/0z1DY9vWmqzcvRVMWTJw40XPnzm3ovT09PQD09vY2HUeqT+nK5sVujh0i/sHUzbFDxN+iGObZnli+vmvPoFpp2bJlHHrooSuW991335b+8qhXt+7g0N2xQ8Q/mLo5doj4GzFjxgxmzJhRXFWx21+cQdE5Z1AhhDAUVTuD6uZOEiGEEFZjUUGFEELoSFFBhRBC6EjRSYLO6SQRQghDQXSS6IfoJBFCCIMnOkmEEELoKk1VUJLe2apAQgghhKJmr0GdJGk34DLbv29FQIMhrkGFEEL7tOUalKStgfuBqcC2wM9sz6j5pg4U16BCCGHwDNQ1qDWBF4FngV2Bj0k6T9L7msw3hNBhehcuqTgsjm16FzY0000INTVbQV0BLAC2B95ne4rto4DXNx3Zaiz+0UO36V24hKnT7+S0mQtW2ndtc9rMBUydfmfsuwNkKB8vmq2gFgI72P6E7bsBJK0DbNh0ZKup+EcP3WjyuFFMmzSW6XMWc9rMBUzu6WFyTw+nzVzA9DmLmTZpLJPHrTpjd2hOK48X3VjRNdtJ4lrbTwJI2gH4kO1PAkc1G1g7tbOTRPEfHeDms48A4C3/9fX4Rw8dSxInT9kGgOlzFjNszO4rnk+bNJaTp2yzYtqG0DqtOl6UKrry76pU0U2fs5hLp+1Iz/jRdcXVu3AJk8eNWuU7t83sRUv7zKfeThLNVlBbFgKbL+nAJvMbFCNGjODCCy9sS1nxjx661cr7bloX++zAatXxopU/jFtR2ZWfBFx00UXLKqVrqIKS9AngGOAVkqYCAl4Armskv6Em/tFDN6n2axngocefGYSIhpZWHC9a+cO4vLIr5TsQzb0NVVC2zwXOlfRe299rSSQNkPRWYD9gSQrLnyl7fT3gLOAvwFbAGbYXtT3QELpU8dfy/+yzNZ+97o/pAPfwXP4x/DXcsAAOvXweFx40IX5cDYDSj4NVOL22+79sXHderfphXF7ZXXPNNQA8ucnElv/QbulYfJIOsP3dlmVYu6z1gd8B29p+VtL3gfNt31RIcwLwT9tnStouv/5v5XmNHy1/c/9Vy+iZPLnPOObPnw/ADjvsUHfsNpz2l92ZvnQiU/gpADPZi2mj5nLyprOI//PQKYr76jYve4QFz2zMFH6KgevYa8W62Hdbr3f5WKb++b1M3WguAqY/uvLxAmD6a7/H7iMW151n8fsEmvreGsmrd/bsiut3v4DW3Qcl6cr89z5J9+bHfcD5jeTXoF2A+20/m5fnAPuUpdkHuB3A9l3AGyQNb1+Iqyp+qdNGzeWN/I6PcAXTRs1l+tKJnPaX3bFTut7lYwcz1BCQ4ORNZ/G2EYtY8MzGbMHiFZXTtFFzmTnushX77uwnxw52uKuVycMWM3WjuVz66ESmPzqRaRvN5SNcQfGU4pblW1DvOUb5D+Mp/HSlY04narSTxKn575dtn1daKenwpiOq32jgycLy8ryunjTLi4mGbT6BnvMbG0nimAt6AOj9Qm9d6WcvXML0+anJZPJWOzL10okrmk+47o9MnwNv3ucwbrn7UabP71/PmhAGgoALbQ69fC43LID7GJubcvZGOoWTbSYvWkrP+PLfh6EZAnr+tIRLL70TgO/dBdfwnzy5yUSm7joGIab/Aibve3ifx4gVHRiWpmtEN5+dRqab9l+Hp+a+cfv3q2mumN+wh9OxczoT+8yn54AqGV5QOX2j16BK13EukTTM9pOSXmn7gkbya9ASYFhheXhe1980bdUzfjSXTttxRbty8WLjTWcfwagRY7nl7rHR5Tx0nM02XL/ieknxI2qA9PzLaC6duiO9i5Zw6S/SutJ1HoDJ40fVte1nL1q6UoeIU/btBVhxX9T0OYuZPK6+vIq99VJl9yUA9t9//5U6TrTiOlSz3cyvAC4HfgBMlrSN7c81HVV9bgfGSFo3N/NNAs6XNBJ4wfZyUq/CXYBb8zWo3+b1g6q4E0SX89DpigekFb+W80X22EcH3uTxo5h999JV1vfnh0Hxh3Hx+yp1eKi3coLWVnZ9abaC+pXtHwDY/kEePLZP+abe39t+odGCbT+dmxS/Kmkp8DvbN0k6E3gMOAM4FzhL0knA64CPNlreQIku56GTlf9aPnnK3sBLXYohKqmB1MofB9UqjP6eAbeysutLsxXUK/tYXkHSVcA5wDuAtwF3AYc1U7jtG4AbytYdX3j+DHBEX/nEdBshVFb+a7l0QCqe+bfygBRe0s6mtP5qtrJr13QbhwBHA38GtgDOs31xlbTH2j5L0h2kkc+Pt316w4W30GBOt1HpF9JA3E8QQqOaHdYmNKbVwxN1smrTbTR1BmX7Ikm3kkYvv8v2whrJt5B0EDDf9guSXtZM2auDTv6FFEJJq5qGQv+0symtUzU75buA3YDJwG6qfSS9iTTqwxckTWm27NXBKs0npK6lJ0/ZZkXvvtmLVr04GkIYGnrGj674A3Wo/DhotonvHGAd4B7SUELP2j6mzvfubPuXDRfeQoPZxBfNJyGEoW5AmviApbY/Xyjk1AoFX1LlvdsDqwQ0GAazk0Q0n4QQhpp2TbdRnmmlYYT+SbpXqtxBTZbdMu2cbiOEEIa6AZ1uo+BuSfOBxeRefBXSHGP7qfKVkh5tsuwQQgirsVb04rsN2JYqvfhKlZOkDUhnTaWxe94MvLWZ8kMIIay+mu5JZ/uPtq+xvVDSoTWSXgCsR+qS/jDwRLNlhxBCWH01OqPuY6xawYh0DaraxZy7bH9F0jr5zKvqqBPtFiNJhBBC+wx0J4kjbV9ZvlLSB2q8Z7ykYcAoSbsBu5PGyxt00UkihBDaZ0A7SRQrJ0nbk64rLQT+t8bbrgV2AK4ELs6PEEIIoaKmOklIOg7YG3gAuBQ4CvhUpbS2f1xY/FdJmzRTdgghhNVbs93MN7C9u6RP2Z4lafdqCSV9qGzVvsB7myw/hBDCaqrZCmrN/Ndly5VMA3rz882BhueCarXoJBFCCO3Truk2PgO8CVifNJX6r6tNoSFpK9t3F5Y/bvv8hgtvof6OxVccP684Fl+MnxdCCP1XbSy+hu6DkvR1ANunAF8mdYD4Zh/zOz0rafP8eD3Q00jZufyRki6UdIKkb0nauEq610n6kaRrGi2rXGmOltNmLqBYuZemzpg6/U56Fy5pVXEhhDBkNdrEN1bSBaQhji61/fM63jMbuI90v9Ry4KIGywY4HbjR9tWS9gXOovLYfm8CfgK8vYmyVjJ53KgVU2HAS22bxXmdJo8bVfX9IYQQ6tNQE5+kDW0/LmksMBXYCJhh+/oa79mz1uv9LP9BYFfbD0oaCdxje2SVtFOBKbb3r5FfxY1QbdvELLghhNB/NY6NrZtuw/bj+e9i4NR8A+7lki6w/doq71mpcpJ0mO1vVCtD0vVApaa7k4HRwJN5eTmwoaS1bLel40VpRkuA6XPSuqicQgihtRod6mi07SWStgUOBQ4Efg+cVCFtaViktYGRwFJSBbMMqFpB2d6zRvlLgGE53+HA481UThMmTKDRCQtDCCHUp1qrVLUf9o1eg7pC0vrAeODbpOa2VUYyz460faWko4BLbP89j2x+TINlA1wH7AI8CEzKy0haA9jM9gNN5N2nSk18pTOpOIsKIYTWaHQ081cBXwc2tX1sjcqpOCzSxrb/ntc9ReXJDet1IvA2SScB+wHH5vXbkysrAEnvIt0Q/C+Sjm+ivBWKldO0SWMZef8sRt4/a0XHifLefSGEEBrT6BnU4bbn9PM9W+ehke4GxgFbNVg2th8DDqmwfj6wXWH5x8CPy9M1Y/aipSsqp5OnbMOss9P6l65JLWbyuFFxL1QIITSp0U4S/a2cIFUoJ5Duf1oAHNxI2QOhPyNJ9IwfzaXTdlxxo25JqeNEVE4hhFBbW0aSaIaknW3/clAKL9PfkSSKiiNJhBBC6L9qI0k0OxZfeSGb2H64bN0Rtr8u6ZKy5NsDqwQUQgghQOPdzMtHJi+pNEL506W3kabkKKk08kMIIYQANH4G9WHglgrrVxnNwfb0/PRo26Wba5E0r8GyQwghDAGNVlBH2/5D+co8CGw1x0n6CWmqjXOAc4EvNlh+S8V0GyGE0D5t6yQhaWvSlO8AB9lepft3TneS7c9Jug14D6mS+5+mCm+R6CQRQgiDZ0A6SUg6kzSaxKuBRcDWNZKvK2k8sMT2UklP10gbQghhiGt0JImSZ2y/C7jG9geB79RI+0/gZuAcSVOACU2WHUIIYTXWbDfzdfLfDSWtRY1KJ09ueEph1cwmyw4hhLAaa7aCej5PGDiXNP3Fd6slzHNHnQ08RRov7yHbv2iy/JaIThIhhNA+bR9JQtII4Dnbz1R5/SLgCmB30oy459j+eEsKb1J0kgghhMHT0k4Skjay/aikN5e9dBAVBnHN/mR7tqQ32X5O0kONlB1CCGFoaLSJ73JgL9K9TL8hjRIBhZHEK9he0s7Aevl+qbENlh1CCGEIaHQ0873y06Nt31paL2lSjbd9EbiYNAbf24EZNdKGEEIY4prqJFGsnLLtgJWm4pC0DunG3CW2d83r3kO6DhUjSYQQwhAzoJ0kJD0GPAGsTRp/bykwGlhme5OytFeTuqOPBL4FfAxYBpxq+1f9LnwARCeJEEIYPK0eSeJI21dKOgq4xPbfJW0AHFMh7V9tH5ODmAOcYLvSQLN1kzQSOAO4lzQz74m2HylLs2OO5zek0S7usH1RM+WGEEJon0avQV2Zn25s++953VOShldI/ljh+U9LlZOkPWzf1Ej5pObBG21fne/DOotVp+/YBDjX9h2S1gaWSPqh7UfLM5s3b95Ks+OWDNZkjiGEsDqqdJytpdkbdbeWdBxwNzCOdDZTbs98dgWwU6ES2xlotILaB/h8fj4H+HZ5AtvXlq16AXi+wfJCCCG0WbMV1CHACUAPsAA4uEKa54C/5+ezCutrVhaSrgc2rvDSyaTrXaW5pZaTh1qy/UKV7I4ETre9rNKLEyZMoNFrUCGEEOpTrVWq2plVs734HgOOLxTyr8DfypIdb/vOCgHVHCzW9p7VXpO0BBhG6qgxHHi8WuUk6QPAy21/rlZ5IYQQOkuz021sBhzFS/NBbQ+s1BOjUuWU1zczo+51wC7Ag8CkvIykNYDNbD+Qlw8GNsjzUG0HPGt7URPlhhBCaJNmp9s4A/glqbnuf0k95trhROBtkk4C9gOOzeu356XK6l2kwWnfLakXuJI0b1UIIYQu0Ow1qPm2fyhpnO0b8lBGAy43La4y5p/t+eThlmz/mCo3f5WLG3VDCKF92jKauaT/JXWSOAqYD3zE9lsaznCQxI26IYQweKrdqNtsE985wAbABaSmtvOazC+EEEIAmu/FVxyqaD9JBzQZTwghhAA0eAYlaR1JB0javbDuPcCprQoshBDC0NboGdQV5AFgJW0OHEa6J2lqa8Jqr+gkEUII7TPQo5mfUzYA7H83OwDsYIpOEiGEMHha3Umi6gCwDeYXQgghrKTRJr6BGAA2hBBCWKHRCqqhAWBDCCGEejVaQTU0AGynik4SIYTQPm0ZSWJ1EZ0kQghh8AzUSBIhhBDCgIgKKoQQQkeKCiqEEEJHana6jdVCdJIIIYT2iU4S/dApnSQkAdCN30k3xw4R/2Dq5tgh4m9RDBU7SXTlGZSkkaTZfO8FtgJOtP1IWZrRwHTgNmA0aezAo2z/s5WxRO+9EEIYGN16Dep04EbbZwA/As6qkGYt4Ee2v2D7P4HdgF36U0jZKWhT6erNq1Xl1ZuulXnVq91xdeq2rzddp8Yf+07/DIVtX2+6evPq1gpqH+D2/HxOXl6J7b/avgggD8u0AXB/pczmzZuHpFUesaP3P12r8hoK277edJ0af+w7/TMUtn1f6UrH1ne+850rHWur6dgmPknXAxtXeOlkUpPdk3l5ObChpLVsv1Ahn/cDhwNn2n5ooOINIYTQWh17BmV7T9s7VHhcCywBhuWkw4HHK1VOOZ+rgN2B/5C0d6U0EyZMwPYqjxBCCK1TOrYecsghdR1ru7IXn6RvADfbvlrSvsD7bB8kaQ1gM9sPSJoMPGP7jvyebwN32v5ahfyWUrn5bwSwrI6Q6klXb171GIy42h3/UNj29abr1Phj3+mfobDt601XnmaM7VHlibq1ghoJfJFUqWwJnGD7EUk7AJfb3k7Sm4Bjgd+QzrY2BQ63/fcq2YYQQuggXVlBhRBCWP117DWoEEIIQ1tUUCGEEDpSVFAhhBA6UlRQIYQQOlJUUCGEEDpSVFAhhBA6UlRQIYQQOlJUUCGEEDpSxw4W204bbbSRx44du8r6ZcuWMWJExYke+52unjTz5s0D0tiAzebVyrjqSVdv7O2Oq940g7Ht603XqfHHvpN087aHzoh/3rx5j1Ya6qjiIKlD7TFhwgRXcsghh1Rc30i6etIATl9J83m1Mq560tUbe7vjqjfNYGz7etN1avyx7yTdvO3tzogfmOsKx+Zo4gshhNCRooIKIYTQkeIaFKk99NBDD12xvO+++6541KOedPXm1ary6k3Xyrzq1e64OnXb15uuE+Pv6enhb3/7W0vLbFVeq/u2b3Ve9WplmRtvvPFKx1zS9BuriNHMgYkTJ3ru3LmDHcaKqY+78Tvp5tgh4u+vnp4eAHp7e5vOK7b94OqE+CXNsz2xfH2cQXWQbt3Bobtjh4h/MHVz7BDxD6S4BhVCCF2ip6dnxdnrUBAVVAghhI7UFU18kt4K7AcsIfXX/0zZ6+sBZwF/AbYCzrC9SNKOwDGkad/HA3fYvqidsYcQQmhMx1dQktYHvgFsa/tZSd+XtIftmwrJjgEesH2mpO2AbwH/BmwCnGv7DklrA0sk/dD2o+3+HCGEEPqnG5r4dgHut/1sXp4D7FOWZh/gdgDbdwFvkDTc9rW27yikewF4vryAefPmIWmVRzsMtTblEMLQVek4W+tY2w0V1GjgycLy8ryuv2mOBE63vazlEYYQQmi5bqiglgDDCsvD87q600j6APBy21+pVMCECRMqjtEXQgihdSodZ2sda/usoCS9s6UR9t/twBhJ6+blScB1kkZKGp7XXUdqCiRfg/qt7eV5+WBgtO3PSdpO0rg2xx9CCKEB9ZxBnSTpTEmvH/BoKrD9NHA48FVJnwN+lztInAB8PCc7l1SJnQT8F/BRAEnvAs4G3i2pF7gSeHV7P0EIIXSWbrn2XU8vvg8D9wNTJR0O/Mz2jIENa2W2bwBuKFt3fOH5M8ARFd73Y6qM8RRC6J/ehUuYPG7UKhe1bTN70VJ6xpdf9g2hOfWcQa0JvAg8C+wKfEzSeZLeN6CRhRA6Ru/CJUydfienzVyw0jUD25w2cwFTp99J78LyS8MhNKeeM6grSB0QZgLvs303gKTTgKsHMLYQQoeYPG4U0yaNZfqcxUCa4Q7gtJkLmD5nMdMmjWXyuFUnRA2hGfVUUAuBg22v6MYtaR1gwwGLKoTQUSRx8pRtAJg+ZzHDxuy+4vm0SWM5eco2bbt3MAwd9VRQ15YqJ0k7AB+y/UngqIEMrJ2qzQcVQnjJypVUWheVU2jEjBkzmDFjpa4MFfsK1FNBbVl6Ynu+pAObjK3jjBgxggsvvHCwwwghhFWsjp1Tyk8CLrrooooDKFTtJCHpE5LuA/5T0r2S7pN0N10wfl8IofVKHSKmz1nMsIfnMuzhuUyfs3iVjhOhdYZ655SqlY3tc4FzJb3X9vfaGFMIocMUK6dpk8Zy89lfAmD//fdf0XEimvpab6h3TunzbKi8cpJ0gO3vDlxIIYROM3vR0pU6RMw6O60vdpyYPG5UVzY3dbKh3jmlagUl6UrbH8jNfKWKW6Rx7qKCasLq2KYcVm8940dz6bQdV9lvSwfQqJwGzlDunFLrRt1T898v235tfmwBnDTwYa2+hnqbcuhePeNHVzwgSorKKQyIWtegFuWnl0gaZvtJSa+0fUGbYmubdnYzH+ptyiGE/invnAIvnUl1wllUaUy/3t7eut/Tym7mVwCXAz8AJkvaxvbn+npTvmfq97ZfqKOMQdXObuZDvU05hFC/1bVzStPdzAt+ZfsHAPlv1f6kkq6StLOkU4GvA1/rT9BDRamSmjZpLE9uMpEnN5kYlVMIYRXlnVNE6ghQOn5Mn7OY2YuW1pVX78IlFW8HsN2xlxXqqaBe2cdy0VzbvwT2BiYDDzQaWAghDHWlzinlP15LP3IvnbZjXdf/uvXadz1NfIsk3QX8GdgCOK9G2i0kHQTMt/2CpJe1IsjVTae3KYcQOke1Cqg/nVO69dp3PfdBXSTpVuD1wF22F9ZIfhNwEPBJSVPojinl22p1bVMOIXSubr323WcFpRT1bsAbgBGSFrnKuCb5GtUP8uJ9kh5tWaSribjhMYQwGLrxfqp6mvi+AqwD3ANMBLYDjikmkHRJlfdun98TsrjhMYQQ6lNPBbXU9udLC7mHXrl/krqilzuowbjaqt3TbbSiTTmEwVC65yV0n0669t3K+6DK3zi8QppjbD9VvrJbmvhiuo0Qwuqs065913sfVD0V1N2S5gOLqdKLr1Q5SdqAdNZU6g7yZuCt/Yg7hBBCi3Xrte96e/HdBmxL3734LgB+TerxdwPwRCuCDCF0lv4MaxMGX7de+66rG7jtP9q+xvZCSYfWSHqX7a8A82xfBMxtSZQhhBCa0srBfts1KkWt6TYeY9UzoNJ0G9Uu2IyXNAwYJWk3YHfgjGaDlPRWYD9gCWDbnyl7fT3gLOAvwFbAGaXBbvMU9W8EXgT+bPubzcYTQghDVWlUilJzYUnxOle9I1z0pVYT35G2ryxfKekDNd5zLbADcCVwcX40RdL6wDeAbW0/K+n7kvawfVMh2THAA7bPlLQd8C3g3yRtBhwLvNG2Jd0p6WbbdzcbVwghDEXtHJVCVe65XTmRtD2p48NC4C/VbtSt8L5NbD/cVIDSHsCJtvfIy58ENrP9yUKaW3OaW/PycmAz4L3ArrY/mtd/FbjH9leLZYwfLX9z/1XL7pk8uZnQ6zJ//nwAdthhhwEvK4Rmxf46uFq1/ZvNx4bT/rI705dOZAo/BWAmezFt1FxO3nQW1ToD9s6eXXH97hcwz/Yq98zWM5LEcaTBXx8ALgWOAj5VJe2HylbtS6okmjEaeLKwvDyvqydNPe8NIYTQDxKcvOksAKYv3Qugz8qpEfV0M9/A9u6SPmV7lqTda6SdBvTm55sDrZgLagkwrLA8PK+rJ80S4HVl6+8pL2DY5hPoOX9w+nMcc0EPAL1f6B2U8kPoj9hfB1ertn9L8rFh5gJYujgtj9sfppxMrRqq54AqL1xQ+T31VFBrlsIpW67k0OL1HUkfryP/vtwOjJG0ru1ngUnA+ZJGAi/YXg5cB+wC3JqvQf3W9nJJ1wNHSVJultyF2qOxhxBCR2pk5tqB0q5RKeqpoF6U9DNgfUk7ke5zquZZSZvn58OBHuD8ZgK0/bSkw4GvSloK/M72TZLOBB4j9RI8FzhL0kmkM6aP5vc+JOks4CuSXgQujg4SIYTQuHaOSlGrm/nXbR9h+xRJbycN/Ppb2zfUyG82cB+pO/py4KKmIwRymTeUrTu+8PwZ4Igq772CNG19CKEBvQuXrHKDJ6QD1exFSzvyBs8wcNo5KkWtM6ixki4gDXF0qe2f15HfYbavbzqqIaQTTtdDqKad97yE7tDOUSlqVVAH2n5c0ljgcEkbATNqVUDlr0k6zPY3WhJpCKHtunUm1jCw2jUjQ9UKyvbj+e9i4NQ8QsTlki6w/dqyoEqjTqwNjASWkrpzLyPdZBtC6ELdOhNrWD3UugY12vYSSdsChwIHAr8HTqqQ/EjbV0o6CrjE9t/zyObHDETQrdbu+aBC6CbdOBNr6GytmA/qijzM0Hjg26QRGSqOZF4YEmlj23/P656SVGnuqI4T80GFEEL7tGI+qFcBXwC+b/u5OsvdOo88cTcwjjRwawihi3XSTKyhNbqlc1at6TYOt/2//aicAA4hjdlX+ntwM8GFEAZX+T0vI++fxcj7Z63oOHHazAUVp10IoRVqdZKY09/MbD8GrLg/SdLOwN8aC23o6aQ7xUOA7p2JNawe6hlJok+SjrD9dUmXlL20PbDKCLUhhO7QrTOxhtVDXTPqFknapMLqp0svkzpUlB7zG44shNARWjkTa+ifds1c26lqdTMvnzqjZJUpNGxPz0+Ptr1iegtJ85qOsA2im3kIodOszqN4tKKb+YeBWyqsH1njPcdJ+glpqo1zSIO4frFmpB0gupmHEDrN6jyKRyu6mR9t+w/lKyW9vsZ7nrP9yzyC+BuAo+uMN4QQQkGM4lG7F9+KyknS1qRu4wAHkbqRV7KupPHAEttLJT1dJV0IIYQ+dMMoHgPZ67ieKd/PJI0m8WpgEbB1jeT/BG4G/kPSFGBCK4IMIYQw9NTTzfwZ2+/KU75/UdJ/VUto+xTglMKqmU1HGEIIQ9RQH8Wjnm7m6+S/G0paixpnRZLGSvq+pG9Lep+kXVsSZQghDDExikd9FdTzkvYF5gJPArWGPvo08FXSrLo/Io2AHkIIoZ/KR/EQ6UbTk6dss6KSmr1o6WCHOaD6bOKzfXLpuaQbqF1B/cn2bElvsv2cpIdaEeRAi/ugQgidZnUexaPp+6AkbWT7UUlvLnupVi++7fP4e+vl7uhj6w958MR9UCGETtSumWvbrRX3QV0O7EW62fY3pLNLgO1qvOeLwMWkMfjeDsyokTaEEEKoqtZ9UHvlp0fbvrW0XtKk8rSS1gHeQ7r/ade87j3A6XTBSBIhhPrESPuhneq5BnVr2artgPKpOK4g9fYbKWlz4GPAMmBqC2IMIYQwBNW6BvUY8ASwNmn8vaXAaFLF842y5H+1fUx+3xzgBNuVxvHrF0kjgTOAe0mz855o+5EK6Q4E3gi8CPzZ9jeVripeRrq5eA1gS9IkjH9vNq5W6l24ZJWLoJC6mM5etLSr25lDCKEZtbqZH2n7tcCZwGjbY0kV1NcrpH2s8PynpcpJ0h5Nxnc6cKPtM0jd1s8qTyBpM+BY4FjbxwMHS9qK9Nnutf1Z258B/g4c1mQ8LVUarbj8fobS/Q9Tp985JIbUDyGESmpdg7oyP924dNZh+ylJwysk31PSBvn5ToU0OwM3NRHfPsDn8/M5pDmmVikbmOeXjvC3A3vZ/iorj2qxBvBUpULmzZtX8Y7sgb4JbnUerTiEEMr1d+SLeoY62lrSccDdwDhSU1u550hnKACzCuuf7ytzSdcDG1d46WTSGVtpfqnl5NEsbL9QSFdMU0q3UruYpLHAa+mw0dVjtOIQQqiungrqEOAEoAdYABxcIc3xtu8sXympz8Fibe9Z7TVJS4BhpGthw4HHyyongCXA6wrLw4F7CnlsBnwBOMD2s5XKmTBhAnPnzu0r1AHRDaMVhxBCK1Rrlap2rOtzqCPbj9k+3vY+to8DxlRIs0rllNc3O6PudcAu+fmkvIykNXJvQYDrgQl66RPuAvw0p9uSVDl9zPZjkv69yXhCCCG0ST3TbWwGHMVL80FtD0wcyKAKTgS+KGkcqRfesYUYLge2s/1QniDxK5JeBC62fbek9UgzAv8FuDbXX3cD329T7HUZ6qMVhxBCNfU08Z1BOqi/A7iG1JW7LWw/RoVhlWzPpzCihe0rSPdiFdP8A9h0gENsSvloxTef/SUA9t9//xUdJ6KSCiEMVfVUUPNt/1DSONs35LH2QguUj1Y86+y0vthxopsHhAwhtM5QHMWjngpqgqQxwKh8Q+zuwGcHNqyhYXUerTiEEJpVz3xQ5wAbABcA+wHnDWRAQ03P+NEVm/C6fbTiEEJoVj1j8f2qsLifpAMGMJ5BEfNBhRBC+9Q7H5Rq9EsvjlA+K697D3C67a1bG+7gmjhxogfrPqiinp4eYGi2NYcQhi5J82yv0ju81hlU+Qjlh5FumJ06EAGGEEIIRbUqqPIRyj/VihHKQwghhHrU6iQxUCOUhxBCCH2qdQY1UCOUhxBCCH2qVUE1PEJ5CCGE0KxaFVTDI5R3m+hmHkII7dN0N/OhJLqZhxDC4KnWzbyekSRCCCGEtosKKoQQQkeKCiqEEEJHqmc089Ams2fPHuwQGlYa8LZbr2lG/IOnm2OHiH8gxRlUCCGEjhQVVA1l3SCbSldvXq0qr950rcyrXu2Oq1O3fb3pOjX+2Hf6Zyhs+3rT1ZtXVFC8dB9U6VHaeLGj9z9dq/IaCtu+3nSdGn/sO/0zFLZ9vem+9rWvrXTMpcp9UHENChgxYgQXXnjhYIcRQghDwpgxY1Y65l500UXLKqWLM6gQQggdKUaSACQtBe6v8NIIoGLN3kC6evOqx2DE1e74h8K2rzddp8Yf+07/DIVtX2+68jRjbI8qTxQVVAghhI4UTXwhhBA6UlRQIYQQOlJUUCGEEDpSdDNvIUlrADOAXwHrAFsCHwH+BfgEsADYFvgf2w9IWge4EFgMbAz81fZnJW0InA/8FhieH5+w/WJZeT3AOcATedV1tr802PHnvH4EvKKQ/XttL61Q3unAU8AY4Fu2f9kBsf8FuDtnvQ7wnO2esvJ6GNxtvy5wLvAgsAlwm+2rcl47AEcA9wGjgWNtv1ChvMHa9hVj76L9vta2/xFt3O8HIP627/s12Y5Hix6kM9KTCss/Bj4IzAfemNftC/w4Pz+g8HwN4HFgU2Az4LBCPj8F/r1CeT1AT6fFn5dPraO89wPn5+cjgUXAmh0Q+wGFfKYBB3bgtj8GODc/X4vUC3UTQMDvgVfl184GPtph275a7N2y31eMfzD2+wGIv+37fq1HNPG1kO1/2v4cgKS1SP9wC4GtgAdysnuBPfLzR4CN8vPhwF+Bx2w/ZPsbhazXIP3aquQgScdKOk3Sazoh/rz8GkknSvq0pIOqFLkPcHsu+zHgH6RfeoMau+3vFrJ+L1BcLhrMbb9ivdPZ0RJgEvBa4GW2/y+nm0PazuUGc9tXjL2L9vtq2x7avN+3Ov7B2PdriQpqAEjaE5gJzLQ9F7gN2Dm/vBOwvqS1bPcCv5Z0GXAV8G3bz5TltSPwJPDzCkUtAD5r+yzSjnRDPt3vhPi/aft0258HplT5Zx2dP1vJ8rxusGMv5bU78Avbz1coalC3fXG9pFcA40gVbb3bdNC2fY3Yi3l17H7fR/yDst+3MP5SXm3f9yuJa1ADwPb1wPWSLpP0ceAg4GhJryPtmH+1/YKko4F1bH8of8G/kvR72z8BkLQdcCTpNHuVG9ZsLyk8/0Pe2V5D5ZuO2xq/7TsKWc4C3gJcXlbUEmBYYXl4XjeosRey+xhp+1cqZ7C3/VXAcEkn5vV3kn4V17tNB3PbV4sd6Ir9vmr8g7Xftyr+grbv+5XEGVQLSdpGUrE55T5Sk8smtk+yfS7wEC/tsK8BHoZ0mk5qdlov57UjaQc5BLCkKXn9yyWNys9PkDQyPx9Juqj5yGDHL2mYpJML+WwF3FMeP3AdsEsh/vWAPwxm7IX8Xgsss/1oYV0nbfsNgZ/aPh24AFgXuIXUlPOMpFfldJNI27mTtn212Ltlv68Y/2Ds962Mv5BfW/f9mp+twg+U0CBJWwJfAn4NrA1sDRwNHEzaYRaQhvg43XbpIHIe6aL2y0i/pI4iXTD+A/A74EVgTeAm26dKOgLYzvZhkt5Pas9eAGwDXGX7ug6If23gO/kziNTx40jbT5fFvwbwBeBpYHPgIjfek6wlsTv3GJP0FVLvqt8Xyuikbb8d8GVgNunX+BW278p57UD6Hu4nXYQ/Nv9y7pRtXzF2SZvRHft9tfjXo837fSvjL+TX1n2/5meLCiqEEEIniia+EEIIHSkqqBBCCB0pKqgQQggdKSqoEEIIHSkqqBBCCB0pKqgQQggdKSqoEEIIHSkqqBCGKEmvkTRL0h8l/UHSJyqkWVPSbyTNLKxbLOkuSfMlzS2sf4ekhZLukXRCWT4VX2vkPf1NE7pX3Kgbuo6knYAzSUOs/Bx4FfBK4IO2nxvM2IokHWP7nD7SzAXe5LI5j6qkfQtphIBxwAnA/7HydnhZTvoF20/Ukd8mpOFwfi1pGDAPeLftBYU0nwQmAsNtl4YdWgxMLBsKZ03StBFvIw2rcyfwH7YXVHuNNOJ2v95TFlufaUJ3izOo0HXygJy9pNGWT7V9GGnqjHcMamCrOqaONDvWUzkB2L4Z+BGwwPa1FbbDp/LyzUqjVveV38O2f52fPwn8kTQ8DwB56KF9gIvrCG8n4B7b9+YfCVcB7+rjtUbeU2+ZYTUQo5mHrpfHNnslsFTSBqQpAG4BxgNX2r5R0jTS+GcXkObLmUAayHM8adTnF4BvA58HtifPJgpsb3v/XM5ppP+ZF4EnbZ9ZyPcrpHHPtgamkCrLV0g6FfgT6eymWP4bgM8BX5XUY3txHmh0lXL7w/ZPJZ0CvBX4Wb3vkzQWeCNpVtaSc4DjWXnkbQADP5dk0vQSF5IqtgcLaR4C3pSfV3utkfcU1ZMmdLGooEI32ylfd3gjaZrweaR9+iu5UhoJXA/caHu6pA8Dd9r+jKSJpOah35KmuF6DNHjmM8DOtvcGkPTR/HfPvP7teblX0s8L+c63/UVJXwfeZvtqSWfaPrUUbHn5tufmJrRi/iuV26D7SYOQIulGUhNouU/b/nFOswHwfeAY28vzuinAEtvzlKb4Lppk+6+SRpPmAvoTaXDUcqXrB9Vea+Q9RfWkCV0sKqjQze6wfQaApP8mjeh8ItAjaRfgeWBU2Xv+COA0oRv54v8HSVMOXAbsSZ4iIaf7Vn66PWnCt9KF+AfL8l6U/y5l1TOOquUXbF+l3EaM4aX5id5aK6GktUmV03ds/6Dw0iTgnZL2Jk0HMVzSFbYPtP3XnPcSST8kNbXNIU1hUrIZaZZiSGc2lV6rtr7We4rqSRO6WFyDCquLh4FXk6YYeLXtzwJnV0hX/gv7PNL8Q6+3/TvSGdWWpRclfUTSOnn9Ettn5EpxOukif7V8AV5U8sY+0lGj3H6R9DZShXJjHWkFfAv4o+0vF1+z/d+2N7M9Fng/cLPtA5XmBRqW3/9y4O2kKUvuBLaStEWO+/3AtTm7aq818p6ietKELhZnUKHr5Oa5NwPrSPo0aT/eDjgZ+Cewv6QvAY8BIyT9O2la7THAkZK+aHspgO1Fkv5Cvl5j++eSdpb0BeAfwN/yBfifS9opr3+BVAmckCuEMcBHJF2a49pO0nWkienOAtaUtFF5+UqTzI0BDrN9QpVyK3GF7XASsD5pPqC32H6hjk05iXT97S5J8/O6E73yrMLlNgZ+mOo21iJd4/tZjudIUpPqmsAltv+Qt+kL1V5r8D0/AQ7OzYwV04TVQ3QzD0OWpHVtP5uvGx3lNLNuR5L01nxd7RhgS9tHDXZMIQy0OIMKQ9nZkp4Abu3kyik7QtIepBlMPzXYwYTQDnEGFUIIoSNFJ4kQQggdKSqoEEIIHSkqqBBCCB0pKqgQQggdKSqoEEIIHSkqqBBCCB0pKqgQQggd6f8BUkHhD8m+BwQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "residuals = vel_wasp - my_sin(time_wasp, *fit1)\n",
    "new_x = np.delete(time_wasp, [0])\n",
    "new_y = np.delete(residuals, [0])\n",
    "new_y_err = np.delete(vel_err_wasp, [0])\n",
    "\n",
    "plotting_residuals(time_wasp, new_x, residuals, new_y)"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}