norheim_final_project.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 15.083 Writeup "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction and overview of methods\n",
    "\n",
    "### Motivation\n",
    "Spacecraft design tends to be a laborious and manual procedure, where engineers iterate on evaluating the perfomance while having a consistent design, i.e satisfying so-called budgets (a.k.a. constraints).\n",
    "\n",
    "One interesting part of the problem is that although one could start from scratch to custom design all parts of a spacecraft(e.g. a satellite), companies generally try to save on development efforts by picking components from a database - either internal to the company, or from larger market availability. \n",
    "\n",
    "Building a spacecraft based on commercially-off-the-shelf components is especially attractive within the new paradigm of CubeSats: a standard for cubesat component sizing and interfaces - turning the spacecraft problem into picking the right \"LEGO\" bricks, and reducing the engineering part to simple assembly and integration(saving a lot on development time). \n",
    "\n",
    "Yet it still leaves the frustration of picking the right component to the engineer; normally the procedure will be quite manual and the engineer will iterate through both picking components, and making sure all constraints are satisfied. It can often be very much a whack-a-mole process, leading either to the engineers frustration, or indifference to looking for better solutions once the design converges on something that works - given the feel that any design changes will just be the same as starting from scratch and having to suffer through yet more frustrations.\n",
    "\n",
    "\n",
    "### The MINLP optimization approach\n",
    "Here we present the optimization approach; the components selection process translate easily into a Mixed Integer formulation. The challenge is that most of the physical sizing constraints are highly non-linear, and oftentimes non-linear in more than one variable. \n",
    "\n",
    "With further careful consideration, one can notice though, that many of the relationships are power laws: products of powers, or sum of products of powers. This is very promising, as this becomes a linear problem in log space, as long as any factors in front of the powers are positive. \n",
    "\n",
    "Unfortunately there are still sometimes models that cannot be captured by this reformulation. These are oftentimes univariate expressions, and therefore the approach we take is to use PiecewiseLinear in JuMP to model them.\n",
    "\n",
    "Finally, we might sometimes encounter disjoint constraints, where the choice is between two component with different physical models. This can also be introduced with a simple modelling technique.\n",
    "\n",
    "### The optimization problem\n",
    "\n",
    "It is fairly accepted that mass can be used as a proxy for cost when it comes to spacecraft. We also choose it over cost, as many times the cost models or cost(in terms of \\$) information is proprietary. Therefore we will use mass as the objective function in the optimization problem.\n",
    "\n",
    "Here we will focus on the optimization of an Earth-Observation satellite with a minimum resolution required, while minimizing the satellite's mass. We will also introduce a minimum lifetime constraint to ensure the missions lasts for a minimum amount of time.\n",
    "\n",
    "**Component selection MIP formulation**\n",
    "\n",
    "For the toy problem(which is much less a toy problem now, and getting closer to covering a lot of the real constraints) we will have a catalog for four type of components:\n",
    "\n",
    "1. Antenna (antenna diameter, maximum power, and mass) \n",
    "1. Battery (energy storage, mass)\n",
    "1. Payload (optical diameter, mass)\n",
    "1. Solar panel (solar panel efficiency, surface mass density)\n",
    "\n",
    "Obviously there are a lot of other satellites components that we can find in catalogs; this can easily be extended for the propulsion and magnetorquer models that are currently being optimized over a continuous domain. We keep it restricted for now just out of time constraints(and not because the framework is limitted).\n",
    "\n",
    "For the MIP formulation, we will create binary variables for each component in the catalog, and set a constraint that only one can be picked for each type of catalog.\n",
    "\n",
    "Next we introduce all the models in their natural form. We have included constraints from 7 different but highly interconnected domains:\n",
    "\n",
    "1. Orbits\n",
    "1. Power\n",
    "1. Payload performance \n",
    "1. Communication\n",
    "1. Satellite lifetime\n",
    "1. Momentum management\n",
    "1. Mass models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models useds\n",
    "### 1. Orbits\n",
    "Given a certain altitude h, the semi-major axis is defined as $a = R_e + h$, where $R_e = 6.378km$ is the radius of the Earth. Then the period T is defined as:\n",
    "\n",
    "$T = 2\\pi\\sqrt{\\frac{a^3}{\\mu}}$\n",
    "\n",
    "$\\mu = 3.98\\cdot10^{14}m^3/s^2$\n",
    "\n",
    "The ground station communication time $g$ is given by the relationship:\n",
    "\n",
    "$g = \\frac{1}{\\pi}arccos\\left(\\frac{R_e}{R_e+h}\\right)$\n",
    "\n",
    "This can be derived from equation (5-24) in SMAD(see figure 5-13 for geometry):\n",
    "\n",
    "$\\lambda_0 = arccos\\left(\\frac{R_e}{R_e+h}\\right) $ \n",
    "\n",
    "This is half of the angle that is actually visible from the ground station. If we normalize by $2\\pi$ we get the fraction of the circle, since we assumed circular orbit:\n",
    "\n",
    "$g = \\frac{1}{2\\pi}2 arccos\\left(\\frac{R_e}{R_e+h}\\right)$\n",
    "\n",
    "The daylight fraction time is given by a similar expression:\n",
    "\n",
    "$d = \\frac{\\pi+2\\lambda}{2\\pi} = g + 0.5$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Power"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a solar panel of efficiency $\\eta_s$, and area $A$, the solar panel charging power will be:\n",
    "\n",
    "$ P_t = d A \\eta_s Q $\n",
    "\n",
    "Where $Q=1367W/m^2$ is the solar constant, and $d$ is still the daylight fraction as already defined.\n",
    "\n",
    "Battery\n",
    "\n",
    "$ E_b = \\frac{1}{d} P_t \\cdot T $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Payload\n",
    "The payload is assumed to be an Earth-imaging camera, and we assume that its resolution, which we have been given a requirement on, is diffraction limitted. This gives the following relationship:\n",
    "\n",
    "$ X_r = h \\frac{\\lambda}{D} $\n",
    "\n",
    "Where $h$ is the satellite's altitude, $\\lambda$ the imaging wavelength and $D$ the aperture of the camera.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Communications"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will assume that every orbit we need to downlink all the data that has been gathered by the payload. Assuming the payload is always recording a strip underneath it with length $N\\cdot X_R$ where $X_R$ is the resolution of the payload. The total data $D$ recorded then is just the amount of pixels needed to represent a full sweep of the earth, i.e. $\\frac{2\\pi R_e}{X_R} $ time the width $N$, times the amount of bits per pixel $B$ :\n",
    "\n",
    "$D = \\frac{2\\pi R_e}{X_R}\\cdot B\\cdot N$\n",
    "\n",
    "All this data needs to be downloaded during a ground station pass, so the bit rate $b$ then becomes:\n",
    "\n",
    "$b = \\frac{D}{T} $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Link budget\n",
    "This equation describes the signal to noise ratio of a communication link. Normally the SNR, here through the measure energy per bit over noise energy, needs to be above a given threshold in order to be decoded by the ground station.\n",
    "\n",
    "$ E_{b}/N_{o} = \\frac{P_{Tx} \\cdot G_{Rx}}{L_{other} \\cdot k \\cdot T_{sys} \\cdot  b}\\eta \\left( \\frac{D_{Tx}}{4r} \\right)^2$\n",
    "\n",
    "Since $b$ has been previously defined we can easily replace it's expression into this equation. This has been done for the optimization formulation.\n",
    "\n",
    "Where $r = \\sqrt{h^2 + 2R_eh}$ is the distance from the satellite to the furthest visible point on the surface of the Earth."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Lifetime of the satellite"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Total lifetime $L_t$ is the sum of two parts, the natural decay time due to drag $L_n$, and the time we can keep the satellite from decaying by countering the drag with propulsion, $L_p$.\n",
    "\n",
    "$ L_n = \\frac{H \\cdot m_t}{2\\pi \\cdot C_D \\cdot A \\cdot \\rho \\cdot a^2} $ Eq. (6-28) from SMAD \n",
    "\n",
    "Where $H(h)$ and $\\rho(h)$ are given by tables at the end of SMAD, which have been digitized and rendrered here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApYAAAJfCAYAAADFFZZ2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzde5xdVX338c83hICgQAxYUUpCENBqqzapjtCCULVYtdqKT2tREY3UC9710UfRolW0rVcQq4iKCi1U8FJaL1gBUTFqxqq1KgghQS0IhMgdQpjf88feo+MwlzOTM3POzHzer9d5bWbttfb+7cWZ4cfae6+VqkKSJEnaVot6HYAkSZLmBxNLSZIkdYWJpSRJkrrCxFKSJEldYWIpSZKkrjCxlCRJUlcs7nUAauy+++61YsWKXochSZI0qcHBweuqao/R5SaWfWLFihWsW7eu12FIkiRNKsnGscq9FS5JkqSuMLGUJElSV5hYSpIkqStMLCVJktQVJpaSJEnqChNLSZIkdYWJpSRJkrrCxFKSJGkeGNy4mZMvuIzBjZt7FoMTpEuSJM1xgxs3c+Spa9mydYglixdxxpoBVi1fOutxOGIpSZI0x61dv4ktW4cYKrhz6xBr12/qSRzzOrFMsiFJjfO5epw2Byb5XJLrk9ya5PtJXpZkuwnO88QkFya5IcnNSb6Z5KiZuzJJkqRfG1i5jCWLF7FdYPvFixhYuawncSyEW+E3AO8Zo/zm0QVJngycA9wOnAVcDzwJeDdwEPC0MdocC5wEbAJOB7YARwCnJfndqnpVdy5DkiRpbKuWL+WMNQOsXb+JgZXLenIbHCBV1ZMTz4YkGwCqakUHdXcBLgN2BQ6qqnVt+Y7A+cCjgKdX1Zkj2qwAfgzcAqyqqg1t+VLg28C+wIFV9Y3Jzr969epat25dx9cmSZLUK0kGq2r16PJ5fSt8io4A9gDOHE4qAarqduC49scXjGrzHGAH4H3DSWXbZjNwQvvj82cqYEmSpH6yEG6F75DkGcDeNCOL3wcuqqq7RtU7rN1+YYxjXATcChyYZIequqODNp8fVUeSJGleWwiJ5X2BT4wquyLJ0VX1lRFlB7TbS0cfoKq2JrkCeDCwEvhRB22uSnILsFeSnarq1m25CEmSpH4332+FfxT4Y5rkcmfgd4EPAiuAzyd56Ii6u7bbG8Y51nD5btNos+tYO5Mck2RdknXXXnvteNcgSZI0J8zrxLKq3lRV51fVL6rq1qr6QVU9H3gXcA/g+CkcLsOH7VabqjqlqlZX1eo99thjCoeVJEnqP/M6sZzAB9rtwSPKJhxdBHYZVW8qbW6cUnSSJElz0EJNLK9ptzuPKLuk3e4/unKSxcA+wFZgfYdt9myP/zOfr5QkSQvBQk0sH9VuRyaJ57fbw8eofzCwE3DxiDfCJ2vz+FF1JEmS5rV5m1gmeXCSe49Rvhx4X/vj6SN2nQ1cB/xVktUj6u8IvKX98Z9GHe6jwB3Ase1k6cNtlgKva3/8AJIkSQvAfJ5u6GnAa5NcAFwB3ESzEs4TgB2BzwHvGK5cVTcmeR5NgnlhkjNplnT8M5pphc6mWeaREW2uSPJq4ERgXZKz+PWSjnsB7+xk1R1JkqT5YD4nlhfQJIQPp7n1vTPwS+BrNPNafqJGrWdZVZ9JcgjweuCpNAnoZcArgBNH12/bnNQuHfkq4Fk0o8A/BI6rqo/NzKVJkiT1n3mbWLaTn39l0op3b/d14E+n2OZc4NypnkuSJGk+mbfPWEqSJGl2mVhKkiSpK0wsJUmS1BUmlpIkSeoKE0tJkiR1hYmlJEmSusLEUpIkSV1hYilJkqSuMLGUJElSV5hYSpIkqStMLCVJktQVJpaSJEnqChNLSZIkdYWJpSRJkrrCxFKSJEldYWIpSZKkrjCxlCRJUleYWEqSJKkrTCwlSZLUFSaWkiRJ6goTS0mSJHWFiaUkSZK6wsRSkiRJXWFiKUmSpK4wsZQkSVJXmFhKkiSpK0wsJUmS1BUmlpIkSeoKE0tJkiR1xYJKLJM8M0m1nzXj1HlikguT3JDk5iTfTHLUJMc9Ksm32vo3tO2fODNXIUmS1J8WTGKZ5LeBk4CbJ6hzLHAu8BDgdOBDwP2A05K8Y5w27wBOA/Zs658O/C5wbns8SZIWtMGNmzn5gssY3Li516Fohi3udQCzIUmAjwKbgE8BrxqjzgrgHcD1wOqq2tCWvxn4NvDKJOdU1TdGtDkQeCVwOfAHVbW5Lf9HYBB4R5J/Hz6WJEkLzeDGzRx56lq2bB1iyeJFnLFmgFXLl/Y6LM2QhTJi+RLgMOBo4JZx6jwH2AF438hEsE0WT2h/fP6oNsM/v3U4qWzbbABObo939DbGLknSnLV2/Sa2bB1iqODOrUOsXb+p1yFpBs37xDLJg4C3A++tqosmqHpYu/3CGPs+P6rOtrSRJGnBGFi5jCWLF7FdYPvFixhYuazXIWkGzetb4UkWA58ArgReN0n1A9rtpaN3VNVVSW4B9kqyU1XdmmRn4P7AzVV11RjH+0m73X960UuSNPetWr6UM9YMsHb9JgZWLvM2+Dw3rxNL4I3Aw4E/rKrbJqm7a7u9YZz9NwA7t/Vu7bA+wG6dhSpJ0vy0avlSE8oFYt7eCk/yCJpRyneOfOFmWw7ZbmuK7catn+SYJOuSrLv22munH5kkSVIfmJeJ5Yhb4JcCb+iw2fAI467j7N+l3d7YYf3JRjSpqlOqanVVrd5jjz06DFOSJKk/zcvEErgnzbONDwJuHzEpegF/29b5UFv2nvbnS9rt3Z6JTLInzW3wn1XVrQBVdQvwc+Ce7f7R9mu3d3tmU5IkaT6ar89Y3gF8eJx9v0/z3OXXaJLJ4dvk5wMHAYePKBv2+BF1RjofeGbb5qMdtpEkSZqXUjXVRwbntiTH04xaPq+qTh1Rvg/wI5p5LleNmCB9Kc0E6fsCB44xQfrXufsE6StoJkjfGXhgJxOkr169utatW7fN1ydJkjTTkgxW1erR5fN1xHLKquqKJK8GTgTWJTkL2AIcAezFGC8BVdXFSd4FvAL4fpKzgSXAXwL3Bl7sqjuSJGmhMLEcoapOSrKBZsnHZ9E8g/pD4Liq+tg4bV6Z5PvAscAxwBDwHeAfq+rfZyVwSZKkPrDgboX3K2+FS5KkuWK8W+Hz9a1wSZIkzTITS0mSJHWFiaUkSZK6wsRSkiRJXWFiKUmSpK4wsZQkSVJXmFhKkiSpK3o+QXqSAI8BHgscDOwN7A7cBlwDfJdmve1/q6qf9ypOSZIkTaxniWWSnYCXAH9Dk0ym3XU7TUJ5D2AlzRrdTwXem+RcmqUVL579iCVJkjSRntwKT3I08BPgBJqRyTfRjFjuVlU7VdVeVbWMJvH9HeA5wDnA44GvJjkryd69iF2SJElj69WI5YeBzwBvq6pvj1epmvUmf9x+TkuyC3AU8Frg2cCbZz5USZIkdaJXieXqqvrOVBtV1Y3ASUk+BKzoelSSJEmatp4kltNJKke1v51mFFOSJEl9wumGJEmS1BU9n25opCSLgPsDewHbj1Wnqi6a1aAkSZLUkb5JLJO8GngVzRyWE9luFsKRJEnSFPVFYpnkeOCNwCbgY8DPga29jEmSJElT0xeJJfBcYD2wqqpu6HUwkiRJmrp+eXlnGc2SjSaVkiRJc1S/JJaXAUt7HYQkSZKmr18Sy/cDT0xy314HIkmSpOnpi2csq+oDSfYHvp7kzcB3gDFvi1fVlbManCRJkjrSF4ll63s0639/ZII6RX/FLEmSpFZfJGlJ1gAfpJli6ELgf3G6IUmSpDmlLxJL4JXANcCBVXVFr4ORJEnS1PXLyzsrgLNNKiVJkuaufkksf844a4NLkiRpbuiXxPLjwJ8muVevA5EkSdL09EtieQLwLeA/kzzaBFOSJGnu6ZfE8g7gz4E/AL4M/DLJXWN8pvSmeJK/T/LlJD9NcluS65P8V5K/TbJsnDYHJvlcW/fWJN9P8rIk201wnicmuTDJDUluTvLNJEdNJVZJkqS5rl/eCv8qzRyV3fZymsnWv0Tz1vnOwABwPHBMkoGq+ulw5SRPBs4BbgfOAq4HngS8GzgIeNroEyQ5FjgJ2AScDmwBjgBOS/K7VfWqGbguSZKkvpOqmcjn+kOSHavq9jHK3wq8DvinqnphW7YLzZrluwIHVdW64WMA5wOPAp5eVWeOOM4K4MfALcCqqtrQli8Fvg3sSzOF0jcmi3X16tW1bt26aV+rJEnSbEkyWFWrR5f3y63wjiTZYSr1x0oqW//abvcbUXYEsAdw5nBSOeIYx7U/vmDUcZ4D7AC8bzipbNtspnluFOD5U4lZkiRpruqLxDLJazqoswT4TJdO+aR2+/0RZYe12y+MUf8i4FbgwFHJ7URtPj+qjiRJ0rzWL89YnpDkyqr6l7F2ti/OnA08bjoHT/Iq4J40t7lXA39Ik1S+fUS1A9rtpaPbV9XWJFcADwZWAj/qoM1VSW4B9kqyU1XdOp3YJUmS5op+SSwvBj6a5OqqumDkjjapPAt4Is1LMtPxKuC3Rvz8BeDZVXXtiLJd2+0N4xxjuHy3KbbZua1nYilJfWRw42bWrt/EwMplrFq+tNfhSPNCX9wKp7k1vR74VJKHDBcmCc3k6X8BfKiqXjqdg1fVfasqwH3bY60E/ivJ70/hMBk+XLfaJDkmybok66699tqxqkiSZsDgxs0ceepa3nneJRx56loGN27udUjSvNAXiWVV/RI4HLgN+HyS+7e7PgI8HfhEVf1NF87zi6r6NM0t9WU0Seuw4VHHXe/WsLHLqHpTaXPjOPGcUlWrq2r1HnvsMWHskqTuWbt+E1u2DjFUcOfWIdau39TrkKR5oS8SS4CquhJ4PE0y9sUkpwJHAZ8Eju7yuTYCPwQenGT3tviSdrv/6PpJFgP7AFtpRlbpoM2eNLfBf+bzlZLUXwZWLmPJ4kVsF9h+8SIGVo65ZoakKeqbxBKgqr4HPJVmGqCjgc8Cf11VQzNwuvu127va7fnt9vAx6h4M7ARcXFV3jCifqM3jR9WRJPWJVcuXcsaaAV7xuAM4Y82Az1hKXdKTCdKTPGuSKk+jWSHndTTLPf5KVX18zBZ3P8cDgV9W1dWjyhcBf9ce++KqOqgt3wW4nGbEtNMJ0veheUPcCdIlSdKCMd4E6b16K/w0xn6hJW358EsvHxhVVvzmc5ETORz4xyQX0SSMm2jeDD+E5uWdq4HnDVeuqhuTPI9mWqMLk5xJs6Tjn9FMK3Q2zdvpjGhzRZJXAycC65Kcxa+XdNwLeGcnSaUkSdJ80KvEsqvPTI7jP4FTaNb4fijNNEG30Mw5+QngxKq6fmSDqvpMkkOA19Pckt+RZpnHV7T175YMV9VJSTbQTGn0LJrHC34IHFdVH5uZS5MkSeo/83qt8LnEW+GSJGmumBdrhUuSJKl/mVhKkiSpK3qSWCb5XpInT7PtfZK8N8lruh2XJEmSpq9XI5abgU8n+VGS17TT9owryQ5JHp/kX4ANNBOn/88sxClJkqQO9eSt8Kp6dJIjgLcAbwNOSHI1sA64iibx3JFm2cUHAr8HbA/cSTNV0Rur6poehC5JkqRx9Gq6IarqbODsJI8FngscCjxpjKp3Ad8FzgE+XFXXzl6UkiRJ6lTPEsthVfUl4EsASQ4A9qYZqbwNuAb4n6q6sXcRSpL6yeDGzaxdv4mBlctcilHqMz1PLEeqqkuAS3odhySpPw1u3MyRp65ly9Yhlixe5DrfUp9xuiFJ0pyxdv0mtmwdYqjgzq1DrF2/qdchSRrBxFKSNGcMrFzGksWL2C6w/eJFDKxc1uuQJI3QV7fCJUmayKrlSzljzYDPWEp9ysRSkjSnrFq+1IRS6lPeCpckSVJXmFhKkiSpK0wsJUmS1BUmlpIkSeqKvnh5J8k9gYPbz97A7vx65Z3vAhdU1Q97F6EkSZIm09PEMskA8ELgCGAHIONUrSSXAP8EnFZVN81SiJIkSepQTxLLJPsD7wCeAAwBXwW+DnwbuBq4HrgHzZrhDwQeBRwGvBd4Q5K/BT5YVUOzH70kSZLG0qsRyx/Q3OZ+LXB6VV01Qd0LgQ8kCfBY4G+A9wG7AW+b4TglSZLUoV4llq8F3l9Vt3faoKoKOA84L8lDgT1nKjhJkiRNXU8Sy6p61za2/x7wvS6FI0mSpC5wuiFJ0qwY3LiZky+4jMGNm3sdiqQZ0hfTDUmS5rfBjZs58tS1bNk6xJLFizhjzYDrfUvzUN+MWCY5JMm/J7kmyZ1J7hrjs7XXcUqSpm7t+k1s2TrEUMGdW4dYu35Tr0OSNAP6YsQyyROAzwDbAVcClwAmkZI0TwysXMaSxYu4c+sQ2y9exMDKZb0OSdIM6IvEEjgeuBN4QlWd1+NYJEldtmr5Us5YM8Da9ZsYWLnM2+DSPNUvieVDgDNNKiVp/lq1fKkJpTTP9cszljfTrLYjSZKkOapfEssv0yzb2DVJliVZk+TTSS5LcluSG5J8Lclzk4x57UkOTPK5JNcnuTXJ95O8LMl2E5zriUkubI9/c5JvJjmqm9cjSZLU7/olsXwNsG+S49qlG7vhacCHgEcC3wTeA5xDc9v9VOBfR58ryZOBi4CDgU8DJwNLgHcDZ451kiTHAue2xz29Pef9gNOSvKNL1yJJktT30qyUOMsnTT4yRvEK4BBgI/Bd4Jdj1Kmqem6H5zgM2Bn4j6oaGlF+X+BbwG8DR1TVOW35LsBlwK7AQVW1ri3fETifZkT16VV15ohjrQB+DNwCrKqqDW35UuDbwL7AgVX1jcniXb16da1bt66TS5MkSeqpJINVtXp0ea9e3nn2BPtWtJ+xFNBRYllV549TfnWSDwBvBR5NM4oJcASwB/Dx4aSyrX97kuNobte/gN8cuXwOsAPw98NJZdtmc5ITgA8DzwcmTSwlSZLmul4llvv06LzD7my3I+fKPKzdfmGM+hcBtwIHJtmhqu7ooM3nR9WRJEma13qSWFbVxl6cFyDJYuBZ7Y8jE8ID2u2lo9tU1dYkVwAPBlYCP+qgzVVJbgH2SrJTVd3ajfglSZL6Vb+8vDOb3k7zos3nquqLI8p3bbc3jNNuuHy3abTZdZz9kiRJ80ZPEssk9+/CMfacRpuXAK+keeHmmVNt3m6n8rbThG2SHJNkXZJ111577RTDkSRJ6i+9GrG8LMm7p5pgpvHkJP8FPG+KbV8EvBf4IXBoVY2ekH2y0cVdRtWbSpsbx9pZVadU1eqqWr3HHnuMG7skSdJc0KvE8h+BY4ANST6f5Ogk+41VMck9kxyW5O+BnwKfAm5vtx1J8jLgfcAPaJLKq8eodkm73X+M9otpXjjaCqzvsM2eNNMd/cznKyVJ0kLQk8Syqt5I8+LLR4E/opmw/MdJfpnkx0m+keS/klwJbAa+BLwa+AVwZFU9qqp+0Mm5kryGZoLz79IkldeMU3V4eqLDx9h3MLATcPGIN8Ina/P4UXUkaU4a3LiZky+4jMGNm3sdiqQ+15MJ0n8jgGZi8r8GHgMcBPzWiN1bgP8GLgTOqaq1Uzz2G4A3A4PA48a4/T06jstpbl93OkH6PjRviDtBuqR5aXDjZo48dS1btg6xZPEizlgzwKrlS3sdlqQe67cJ0n+lqm4EPtB+SLI9sAy4rarGe9t6Uu1a3W8G7gK+CrxkjNUiN1TVacNxJHkecDZwYZIzgeuBP6MZXT0bOGtU7FckeTVwIrAuyVk0yfARwF7AOztJKiWpX61dv4ktW4cYKrhz6xBr128ysZQ0rp4nlqNV1Z3AWM9ATtXwJOzbAS8bp85XgNNGnPszSQ4BXg88FdiRZpnHVwAn1hjDu1V1UpINwKto5sdcRPOC0HFV9bEuXIck9czAymUsWbyIO7cOsf3iRQysXNbrkCT1sZ7fClfDW+GS+tXgxs2sXb+JgZXLHK2UBPThrfAke0+nXVVd2e1YJEnjW7V8qQmlpI708lb4BqY22Tht/b67fS9JkqTeJmlXcvfEcjeaycZ7tpa4JEmSpqdniWVVrRhdluR44A1Vtc/dGkiSJKmv9WrlnfH4JpEkSdIc1W+JpSRJkuYoE0tJkiR1hYmlJEmSusLEUpIkSV1hYilJkqSu6OXKO3dNY19VlROkS5Ik9aFeJmmZpTaStCC5xrek2dbLCdK9DS9JM2Rw42aOPHUtW7YOsWTxIs5YM2ByKWnGmdxJ0jy0dv0mtmwdYqjgzq1DrF2/qdchSVoATCwlaR4aWLmMJYsXsV1g+8WLGFi5rNchSVoAfBFGkuahVcuXcsaaAZ+xlDSrTCwlaZ5atXypCaWkWeWtcEmSJHWFiaUkSZK6wsRSkiRJXWFiKUmSpK4wsZQkSVJXmFhKkiSpK0wsJUmS1BUmlpLURwY3bubkCy5jcOPmXociSVPmBOmS1CcGN27myFPXsmXrEEsWL+KMNQNOcC5pTnHEUpL6xNr1m9iydYihgju3DrF2/aZehyRJU2JiKUl9YmDlMpYsXsR2ge0XL2Jg5bJehyRJU+KtcEnqE6uWL+WMNQOsXb+JgZXLvA0uac6ZtyOWSY5IclKSrya5MUklOX2SNgcm+VyS65PcmuT7SV6WZLsJ2jwxyYVJbkhyc5JvJjmq+1ckaSFYtXwpLzr0ASaVkuak+TxieRzwUOBm4GfAAyeqnOTJwDnA7cBZwPXAk4B3AwcBTxujzbHAScAm4HRgC3AEcFqS362qV3XrYiRJkvrdvB2xBF4O7A/sArxgoopJdgE+BNwFPLqqnltVrwYeBnwDOCLJX41qswJ4B00CurqqXlRVLwd+D7gceGWSR3X1iiRJkvrYvE0sq+qCqvpJVVUH1Y8A9gDOrKp1I45xO83IJ9w9OX0OsAPwvqraMKLNZuCE9sfnTzN8SZKkOWfeJpZTdFi7/cIY+y4CbgUOTLJDh20+P6qOJEnSvGdi2Tig3V46ekdVbQWuoHkedWWHba4CbgH2SrJTd0OV1I9cMUeS5vfLO1Oxa7u9YZz9w+W7TbHNzm29W7cpOkl9zRVzJKnhiGVn0m47eV6z4zZJjkmyLsm6a6+9dtrBSeotV8yRpIaJZWN41HHXcfbvMqreVNrcON5Jq+qUqlpdVav32GOPjgKV1H9cMUeSGt4Kb1wCrKaZnmhw5I4ki4F9gK3A+lFtdm/bfGNUmz1pboP/rKq8DS7Nc66YI0kNRywb57fbw8fYdzCwE3BxVd3RYZvHj6ojaZ5zxRxJMrEcdjZwHfBXSVYPFybZEXhL++M/jWrzUeAO4Nh2svThNkuB17U/fmCG4pXUJb7NLUndM29vhSd5CvCU9sf7tttHJTmt/efrhpdcrKobkzyPJsG8MMmZNCvq/BnNtEJn0yzz+CtVdUWSVwMnAuuSnMWvl3TcC3hnVf3GLXJJ/cW3uSWpu+ZtYkmzHONRo8pW8uu5KDcCv1rLu6o+k+QQ4PXAU4EdgcuAVwAnjrWCT1WdlGRDe5xn0YwA/xA4rqo+1tWrkdR1Y73NbWIpSdM3bxPLqjoeOH6Kbb4O/OkU25wLnDuVNpL6w/Db3HduHfJtbknqgnmbWErSZHybW5K6y8RS0pw0uHFzVxLCVcuXmlBKUpeYWEqac3zpRpL6k9MNSZpV3ZjexyUUJak/OWIpadZ0a6TRl24kqT+ZWErqSDeeaezW9D6+dCNJ/cnEUtKk+nGk0ZduJKn/mFhK85wjjZKk2WJiKfWpbiSEjjRKkmaTiaXUh7qVEDrSKEmaTSaWUhd1a9LubiWEjjRKkmaTiaXUJd2ctLtbCaEjjZKk2WRiKdFfL7hAdxNCRxolSbPFxFILXj++4AImhJKkucfEUnNaP400ettZkrTQmVhqzurHkUZHGSVJC5mJpeYsRxolSeovJpbqiW7cwnakUZKk/mJiqVnXrVvYjjRKktRfTCw167o9LY8JpSRJ/WFRrwPQwjN8C3u70JVpeSRJUn9wxFKzzlvYkiTNTyaWmpJurYXtLWxJkuYfE0t1rJtrYUuSpPnHZyzVsbFeupEkSRpmYqmO+dKNJEmaiLfC1TFfupEkSRMxsdSU+NKNJEkaj7fCJUmS1BUmll2QZK8kH0nyv0nuSLIhyXuSOLQnSZIWDG+Fb6Mk+wIXA/cBPgv8GHgE8FLg8CQHVZWvT0uSpHnPEctt936apPIlVfWUqnptVR0GvBs4AHhrT6NrDW7czMkXXMbgxs29DkWSJM1TqapexzBnJVkJXA5sAPatqqER++4FXAUEuE9V3TLRsVavXl3r1q2bkTid2FySJHVTksGqWj263BHLbXNYuz1vZFIJUFU3AV8HdgIGZjuwkZzYXJIkzQYTy21zQLu9dJz9P2m3+89CLONyYnNJkjQbfHln2+zabm8YZ/9w+W6zEMu4nNhckiTNBhPLmZV2O+aDrEmOAY4B2HvvvWc0ECc2lyRJM81b4dtmeERy13H27zKq3m+oqlOqanVVrd5jjz26HpwkSdJsMrHcNpe02/Geodyv3Y73DKYkSdK8YWK5bS5ot49L8ht92U43dBBwG7B2tgOTJEmabSaW26CqLgfOA1YALxq1+03AzsDHJ5vDUpIkaT7w5Z1t90KaJR1PTPLHwI+ARwKH0twCf30PY5MkSZo1jlhuo3bUcjVwGk1C+UpgX+BE4FGuEy5JkhYKRyy7oKp+Chzd6zgkSZJ6ybXC+0SSa4GNvY5jlu0OXNfrIPqcfTQx+2dy9tHE7J+J2T+TW6h9tLyq7jZXoomleibJurEWsNev2UcTs38mZx9NzP6ZmP0zOfvoN/mMpSRJkrrCxFKSJEldYWKpXjql1wHMAfbRxOyfydlHE7N/Jmb/TM4+GsFnLCVJktQVjlhKkiSpK0wsJUmS1BUmlpqyJMuSrEny6SSXJbktyQ1JvpbkuUnG/F4lOTDJ55Jcn+TWJN9P8rIk201wricmubA9/s1JvpnkqJm7upmT5JlJqv2sGafOlK83yVFJvtXWv6Ft/8SZuYruS/JHSc5JclWSO9rteUn+dIy6C+47lOQJbX/8rP1dW5/kk0keNU79edVHSY5IclKSrya5sf39OX2SNrPSB/3yuzeVPkqyX5LXJDk/yU+TbEnyiySfTXLoJOeZ0vUm2a7t992f3UMAACAASURBVO+3393r238vB27rNU/VdL5Ho9p/eMTf7weMU2fK15vkHknelOSSJLcnuSbJvyZ50HSusy9UlR8/U/oAzwcK+F/gDOBtwEeAX7blZ9M+vzuizZOBrcDNwIeBfwR+3Nb/5DjnObbdfx1wMvBu4Kdt2Tt63Q9T7LPfbvvnpjb+Nd24XuAd7f6ftvVPBja1Zcf2+ro76Jfj2livBT4KnEDzIPy3gX9Y6N8h4O9HxH8q8Pb292sLMAQ8Y773EfDdNpabgB+1/3z6BPVnpQ/66XdvKn0EnNnu/x/ggzR/vz/V9lkBL+nG9QIBPtnu/3H77+HD7b+XrcCT+7WPxmj7pBFtC3hAN64X2AH4Wtvm2+3v+z8DdwK3AI+czT7qWl/3OgA/c+8DHNb+oi0aVX5f4Mr2l+SpI8p3Aa4B7gBWjyjfEbi4rf9Xo461Ari9/cO1YkT5UuCyts2jet0XHfZXgP8ELm//2NwtsZzO9QIHtuWXAUtHHWtTe7wVM3VdXeiXp7Xxfwm41xj7t1/I36H29+ku4GrgPqP2HdrGv36+91F7rfu1v0ePZuKkaVb6oN9+96bYR88GHj5G+SE0/8NyB7Dntl4v8PS2zdeBHUeU/0F7jmsY4/e+H/poVLs92t/BM4ELGT+xnPL1Av+vbfNJRvz3lOZ/joaT/0XTud5efnoegJ/59QFe1/5CnDSi7Dlt2cfGqH9Yu+8ro8rf3Ja/aYw24x6vHz/AS2lGlw4GjmfsxHLK1wt8vC0/eow24x6vHz40j+Gsp/m/8j06qL/gvkPAI9sYPzvO/huBmxZSH02WEMxWH/Tz795kfTRJ2/MYNTAw3esFLmrLDx2jzbjH67c+Aj5Nk1guY+LEckrXS5PgbmzL95nK8fr94zOW6rY72+3WEWWHtdsvjFH/IuBW4MAkO3TY5vOj6vSt9jmZtwPvraqLJqg6neudy310ILAP8Dlgc/sc4WuSvHScZwcX4nfoJzQjSI9IsvvIHUkOBu5FMxI+bCH20Wiz1Qfzrd+GjfX3G6Z4vW3/HkjT31/tpE0/SvJs4CnA86tq0wT1pnO9+wJ7A5dW1RUdtpkTTCzVNUkWA89qfxz5B+iAdnvp6DZVtRW4AlgMrOywzVU0I117JdlpG8OeMW1/fILm8YDXTVJ9StebZGfg/sDN7f7RftJu959G6LPhD9rtL4DvAP9Ok4C/B7g4yVeS7DGi/oL7DlXV9cBrgN8CfpjklCRvS/KvNCNLXwL+ZkSTBddHY5jxPpgHv3tjSrIc+GOa5OiiEeXTud4HANvRPKoxOkkdr01fafvjvTSjmp+ZpPp0rnfc790EbeYEE0t109uBhwCfq6ovjijftd3eME674fLdptFm13H294M3Ag8Hnl1Vt01Sd6rXO50+7Sf3abfPB+4BPIZmBO4hwBdpHhv45Ij6C/I7VFXvAf6CJhl6HvBammdTfwqcVlXXjKi+IPtolNnog7n+u3c37YjbGTQvkxxfVZtH7J7JPu3LPkozs8nHaF68eUkHTRZcH03ExFJdkeQlwCtp3oZ75lSbt9ua4TazJskjaEYp31lV3+jGIdvtVK+3L/uH5v/uobmuI6rqy1V1c1X9D/DnwM+AQ8abUmcM8+47BJDk/9K8BX4aza2znYFVNM+nnpHkH6ZyuHY7r/poimazD+ZEn7VTMH0COAg4i+bt7+mYT9+rl9O8zPS8UUn2dC2o3z0TS22zJC+iuWXwQ5oHja8fVWWyUY9dRtWbSpsbpxDqrBhxC/xS4A0dNpvq9U5Wf7L/G+614T/W66vqeyN3tKO7wyPej2i3C+o7BJDk0TTTj/xbVb2iqtZX1a1V9R2a5PvnwCuTDN/WXXB9NIbZ6IO5/rv3K21SeTrNKPi/0kxfNTqRmc71TuffQ19Ish/wVuCjVfW5DpvN5Peu7/poMiaW2iZJXga8D/gBTVJ59RjVLmm3d3tWpE3C9qF5WHx9h232pBm5+VlV3Tr96GfMPWnifhBw+4hJdQv427bOh9qy97Q/T+l6q+oWmsTinu3+0fZrt+M9v9Nrw9f7y3H2Dyee9xhVf6F8hwCGJ56+YPSONuZv0fwNf3hbvBD7aLQZ74N58LsH/Ko//gX4K5q5E/96rOcDp3m9l9FMlbWyPU8nbfrFg2keCTh65N/u9u/3IW2dn7RlT2l/ns71jvu9m6DNnGBiqWlL8hqaiXK/S5NUXjNO1fPb7eFj7DsY2Am4uKru6LDN40fV6Td30EyMO9bnv9o6X2t/Hr5NPp3rnct9dBHNf+D3S7JkjP0Pabcb2u1C+w5B8x83aObRG8tw+ZZ2uxD7aLTZ6oM53W/t79zZNCOVHweeWVV3TdBkStfb9u/FNP39R5206SMbGP/v9/DAySfbnzfAtK/3cpoXO/dPsk+HbeaGXs935Gdufmhu8RawDrj3JHV3oVlZZSqTFu9Dn0/cPM1+O56x57Gc8vXSZ5M0T6MvTm/jf8uo8sfSzPv5S2C3hfodAv5PG+PVwP1H7Xt820e3AcsWSh/R2QTpM94H/fy710Ef7QD8R1vnVDqYgHs610tnE4bv0o99NEG7C9m2CdJ3GdVmXk6QnvYipI6lWT/3NJqh/5MY+xmQDVV12og2T6H5P+TbaVYwuB74M5opF84G/k+N+jImeTFwIs0frrNoRmaOAPaieSnmVd28rtmQ5Hia2+HPq6pTR+2b8vUmeSfwCpqXXc4GlgB/STOZ74ur6n0zdjHbKMl9aP4IP4Bm7rdvActpnh8smltznxxRf0F9h9o3U79I88b8Tfx6ouYH0dwmD/CyqnrviDbzro/aaxq+5Xhf4E9obmUPzxd43cgYZ6sP+ul3byp9lOSjNKvvXAe8n7FfDrmwqi4cdY4pXW+S0Dy3eQTNS53ntnX/kibRf2pVfXa61zxVU/0ejXOMC2luh+9XVZeN2jfl623fxj+fJnFfB3yZZm7Lp9F8Bw+rqm9O43J7q9eZrZ+59+HXo24TfS4co91BtBNi04y0/DfN23fbTXCuJwFfofkP6y0066ke1es+6ELf3W2t8OleL3BUW++Wtt1XgCf2+lo77I97A++imV9wC81/2D8LDIxTf0F9h4DtgZcBa2leHtlKM/Lx78DjFkIfdfD3ZkOv+qBffvem0kf8etRtos/x3bhemmmyXt72/23tv4/PAQf2cx9NcIzhvrvbiOV0r5fmOfI30cxbeQfNiPsngd+Z7T7q1scRS0mSJHWFL+9IkiSpK0wsJUmS1BUmlpIkSeoKE0tJkiR1hYmlJEmSusLEUpIkSV1hYilJkqSuMLGUJElSV5hYStIMSrIiSSU5rcvHfXOS25P89kyfq8N4dkpydZJPzPa5JfUPE0tJmmPaZPJVwClV9dNexwNQVbcCbwOOTPKIXscjqTdMLCVp7nkDsAPwj70OZJQPAr8E3tLrQCT1homlJM0hSXYFjgS+3C+jlcOq6nbgLOAxSfbrdTySZp+JpSTNkvYZyDOTXNc+H7kuyROneJinAzvRJHCdnndRkhPb5y8/lWTHtvz4tuzRSZ6eZDDJrUn+N8m7kuzQ1jssyYVJbkyyOcknkiwb53RnAgGeM8XrkjQPmFhK0uxYDnwLWAF8giYxfAjw2SSHTuE4j2m3X+ukcptE/ivwYuBk4Ih2ZHGkFwMfBi4B/gnYBLwc+GCSPwc+D1wPnAL8CHgGcPo4p/wWcCfw2A6vR9I8srjXAUjSAvFo4PiqetNwQZJ/Br4AvBq4oMPj/CFwE3DpZBWT3Bv4LHAQ8Nqq+vtxqj4GWFVVP2rb7QB8B3gm8CTgcVX1lXbfIuCLwOFJHlZV3x15oKq6Lcn/AA9Pcq+quqnD65I0DzhiKUmzYyOjXmqpqi8CVwIdvUWdZAnwW8DVVVWT1F0OfB14JPDMCZJKgBOHk8o2rjtoRlQXAf8xnFS2+4b49WjlQ8c53tVt2/tPfEWS5hsTS0maHd+tqrvGKP8psLTDYww/17h5knoHAN8A7gc8vqrOmKT+ujHK/rfdDo6x7+ftdq9xjnd9u919kvNKmmdMLCVpdvxynPKtdP63+LZ2u+Mk9fYH9gTW09zSnswN48Q12b7txznePdrtbePslzRPmVhK0hxRVb8EtvDrkcvxnAu8DngY8OUksz1yOBzfNbN8Xkk9ZmIpSXPLfwN7JtllokpV9TaaN7sfDlyQ5LdmI7jWATRvlv9sFs8pqQ+YWErS3HIhzd/uSV/4qar3AC8AHgx8Jcn9ZjY0SLIPzQtGF072gpGk+cfEsouS7JnkY0mubSc//mGSQ3odl6R55Zx2+yedVK6qD9BMVr4fcFGSvWcqsNbj2u05E9aSNC+ZWHZJkt1opvYI8ATgQTSTDvuMkaSuqapvAP8FHJlkuw7bnEYzqflymuRy5cxFyFHAtZhYSgtSvFPRHUlOAA6pqoN6HYuk+S3J04F/Bv6iqj7d63iGJfk94HvAG6rqLZPVlzT/zNkRyyRHJDkpyVfb9WsryXhLjA232SvJR9p1cO9IsiHJe5J0OofcRJ4CfDPJWUmuSfLdJMcmSReOLUkjnQl8Ezi+z/7GvJnmhZ139joQSb0xZxNL4DjgWJrpNH4+SV2S7Esz0e/RNGvZvptmjreXAt9IMtn0HZNZCbywPeafAO8F3g68aBuPK0m/oX0p5hjg0zSToPdckp1obtE/s6qcv1JaoObsrfAkh9L8n/FlwCE06+yeUVXPGKf+F2keKn9JVZ00ovxdNFNyfLCqnj+i/C3A6ycJ49CqurCtvwVYV1UHjjjGCcCfV9WDpn6FkiRJc8ucTSxHSvJoJkgs2wfVLwc2APu2a90O77sXcBXNSzf3qapb2vLdmXw5siur6ta2/kbgS1W1ZsSxnwl8oKp2nuwadt9991qxYsVk1SRJknpucHDwuqraY3T54l4E0wOHtdvzRiaVAFV1U5Kv04xmDgBfbsuvA66bwjm+TjMp8Ej7Axs7abxixQrWrRtruV5JkqT+0g6o3c1cfsZyKoYTvkvH2f+Tdrv/Npzj3cBAktcneUCSpwEvAU7ehmNKkiTNGQslsdy13d4wzv7h8t2me4Kq+jbNm+H/B/gB8FbgDcD7x2uT5Jgk65Ksu/baa6d7akmSpL6wUG6FT2Z4uo5teuC0qv4D+I8p1D8FOAVg9erVc/9hV0mStKAtlBHL4RHJXcfZv8uoepIkSZqihZJYXtJux3uGcr92O94zmJIkSZrEQkksL2i3j0vyG9fcTjd0EHAbsHa2A5MkSZovFkRiWVWXA+cBK7j7SjhvAnYGPj48h6UkSZKmbs6+vJPkKTRvYQPct90+Kslp7T9fV1WvGtHkhcDFwIlJ/hj4EfBI4FCaW+CTrbIzpw1u3Mza9ZsYWLmMVcu7sTS6JEnSb5qziSXNGuFHjSpb2X6gmZj8V4llVV2eZDXwZuBw4E9pVtw5EXhTVV0/4xH3yODGzRx56lq2bB1iyeJFnLFmwORSkiR13ZxNLKvqeOD4Kbb5KXD0TMTTz9au38SWrUMMFdy5dYi16zeZWEqSpK5bEM9YLnQDK5exZPEitgtsv3gRAyuX9TokSZI0D83ZEUt1btXypZyxZsBnLCVJ0owysVwgVi1fakIpSZJmlLfCJUmS1BUmlpIkSeoKE0tJkiR1hYmlJEmSusLEUpIkSV1hYilJkqSuMLGUJElSV5hYSpIkqStMLCVJktQVJpaSJEnqChNLSZIkdYWJpSRJkrrCxFKSJEldYWIpSZKkrjCxlCRJUleYWEqSJKkrTCwlSZLUFSaWkiRJ6goTS0mSJHWFiaUkSZK6wsRSkiRJXWFiKUmSpK4wsZQkSVJXmFhKkiSpK0wsJUmS1BUmlpIkSeoKE0tJkiR1hYmlJEmSusLEUpIkSV1hYtklSbZL8ndJrkhye7t9S5LFvY5NkiRpNpj0dM9rgBcBRwH/Dfwe8DHgDuDvehiXJEnSrDCx7J4DgXOr6tz25w1J/g14ZA9jkiRJmjVz9lZ4kiOSnJTkq0luTFJJTp+kzV5JPpLkf5PckWRDkvckWdqFkL4GHJrkge25fgc4DPhcF44tSZLU9+byiOVxwEOBm4GfAQ+cqHKSfYGLgfsAnwV+DDwCeClweJKDqmrTNsTz98C9gB8muYumb99aVe/fhmNKkiTNGXN2xBJ4ObA/sAvwgg7qv58mqXxJVT2lql5bVYcB7wYOAN46snL74k1N8nn0iCZ/CTwL+Gvg99t/fmGS527rhUqSJM0Fqapex7DN2gTvAuCMqnrGGPtXApcDG4B9q2poxL57AVcBAe5TVbe05bsDu09y6iur6ta2/k+Bd1TVe0cc+zjg2VX1gMmuYfXq1bVu3brJqkmSJPVcksGqWj26fC7fCp+Kw9rteSOTSoCquinJ14HHAQPAl9vy64DrpnCOnYC7RpXdxdweFZYkSerYQkl6Dmi3l46z/yftdv9tOMe5wGuTPCHJiiR/DrwC+PQ2HFOSJGnOWCgjlru22xvG2T9cvts2nOPFNPNVDj/LeRXwIeDN4zVIcgxwDMDee++9DaeWJEnqvYWSWE4m7XbaD5xW1U3Ay9pPp21OAU6B5hnL6Z5bkiSpHyyUW+HDI5K7jrN/l1H1JEmSNEULJbG8pN2O9wzlfu12vGcwJUmSNImFklhe0G4fl+Q3rrmdbugg4DZg7WwHJkmSNF8siMSyqi4HzgNWAC8atftNwM7Ax4fnsJQkSdLUzdmXd5I8BXhK++N92+2jkpzW/vN1VfWqEU1eSLOk44lJ/hj4EfBI4FCaW+Cvn/GgJUmS5rE5m1gCDwOOGlW2sv0AbAR+lVhW1eVJVtNM/3M48Kc0UwKdCLypqq6f8YglSZLmsTmbWFbV8cDxU2zzU+DomYhHkiRpoZtSYpkkwGOAxwIHA3vTrKd9G3AN8F3gfODfqurn3Q1VkiRJ/ayjxDLJTsBLgL+hSSaHJxS/nSahvAfNLeh9gacC701yLvDOqrq420FLkiSp/0z6VniSo2nW0j6BZmTyTTQjlrtV1U5VtVdVLaNJUn8HeA5wDvB44KtJzkrieoWSJEnzXCcjlh8GPgO8raq+PV6lqirgx+3ntCS70Lxc81rg2UywZrYkSZLmvk4Sy9VV9Z2pHriqbgROSvIhmvkjJUmSNI9NmlhOJ6kc1f52mlFMSZIkzWMLYuUdSZIkzbyuJZZJ7p1k124dT5IkSXNLx4llkvsneV+SLyb5hyTL2vKHJfk+cC1wfZKLkjxwpgJWbw1u3MzJF1zG4MbNvQ5FkiT1mU7nsbw3sBa4f1v0WOAxSf4E+A9gGfA94H7AHwL/meQhVfXL7oesXhncuJkjT13Llq1DLFm8iDPWDLBq+dJehyVJkvpEpyOWx9IklSfQrNH9xnZ7GnATsH9V/X5V3Rd4G02C+eKuR6ueWrt+E1u2DjFUcOfWIdau39TrkCRJUh/pNLF8CvDtqjquqr5fVW8BvgYcDrymqq4cUfc4YD3wxO6Gql4bWLmMJYsXsV1g+8WLGFi5rNchSZKkPtLpWuHLgdNHla0DDgJ+Y8nGqqokXwH+YtvDUz9ZtXwpZ6wZYO36TQysXOZtcEmS9Bs6TSzvAdwyquwGgKq6doz6vwB23oa41KdWLV9qQilJksbU6a3w64D7jCq7BbhmnPrLAF/ckSRJWkA6TSwvBX5nZEFVvaOq9hyn/r7Az7YlMEmSJM0tnSaWg8CqJEsmq5jkt4CDga9vS2CSJEmaWzpKLKvqNVW1Q1Vt6aD6bsCrgQ9sU2SSJEmaUzp9eadjVXUJcEm3jytJkqT+1rW1wiVJkrSwTWvEMslewMtpVt/ZC9h+jGpVVftuQ2ySJEmaQ6acWCZ5NPA5YEdgK82clVvHqrpNkUmSJGlOmc6I5T8A2wHPAv65qoa6G5IkSZLmoukklr8L/EtVjV7iUZIkSQvYdF7e2Qxc3+1AJEmSNLdNJ7H8d+CQbgciSZKkuW06ieXrgF2TnJxk524HJEmSpLlpys9YVtV1SQ4Hvgk8K8mlwA1jV60/3tYAJUmSNDdMZ7qhBwMXAEvbooePU7WmG5QkSZLmnuncCn8XsAx4I7Ac2L6qFo3x2a6rkWpeGdy4mZMvuIzBjZt7HYokSeqS6Uw39CjgU1X1lm4Ho4VhcONmjjx1LVu2DrFk8SLOWDPAquVLJ28oSZL62nRGLLcAG7ochxaQtes3sWXrEEMFd24dYu36Tb0OSZIkdcF0EssLgUd0OY6+luTgJP+W5OdJKsmzx6n3wiRXJLk9yWCSP5rlUOeEgZXLWLJ4EdsFtl+8iIGVy3odkiRJ6oLpJJb/F/idJK9NslDWA78n8APgpcBtY1VI8pfAe4ETaF5ouhj4fJK9ZyvIuWLV8qWcsWaAVzzuAG+DS5I0j6Rqai9vJ/kIsA9wMM0t8e8y/nRDz93WAPtNkpuBY6vqtFHl3wS+X1XPG1H2E+Dsqvp/kx139erVtW7dum6HK0mS1HVJBqtq9ejy6by88+wR/7xP+xlLAdNKLJMcQbO6z8OAhwL3As6oqmdM0GYv4M3A4TRvrV8FfAZ4U1XN6KvHSZYAq4B3jNp1HnDgTJ5bkiSpX0wnsRwvkeym42gSypuBnwEPnKhykn1pbj3fB/gs8GOa50BfChye5KCqmsk3RHYHtgN+Mar8F8BjZvC8kiRJfWM6K+9snIlARnk5TUJ5Gc3I5QWT1H8/TVL5kqo6abgwybvaY70VeP6I8rcAr5/kmIdW1YVTjHv0cwUZo0ySJGle6iixTPIe4FPAV2uqD2VOQ1X9KpGc7P2gJCuBx9E873nyqN1/CxwDPDPJK6vqlrb8PcDpk4Rx5RRCvg64C7jvqPL7cPdRTEmSpHmp0xHLY4EXA5uS/BtNkvmfVbVlxiLr3GHt9ryqGhq5o6puSvJ1msRzAPhyW34dTTLYFVW1Jckg8FjgkyN2PRY4p1vnkSRJ6medTjd0P+AFwCDwDOBc4LokZyb5yyT3mqkAO3BAu710nP0/abf7T/cESe6Z5GFJHkbTZ3u3P4+cSuhdwLOTrEnyoCTvpem3D0z3vJIkSXNJR4llVV1TVadU1eNpbu8+E/gi/7+9O4+Xo6rzPv75hn0xgAFcUJZoyKiMgFwhgguLBuTREQVndBABBcYVXJgRGecx4VHcRtnEZ8igZBBQBBREUaMYdiPcC+joIFtMBAEhJKwJCSG/+eOclqbpvr3c6v37fr3qVblVp6rOOem++eXUWWB/4DvA/ZJ+JOn9krZoX3ar2iTvq015VH580wk8YwS4KW8bALPzn08oJYiI84GPkQYe3Qy8Fth/vD6pko6SNCpp9IEHHphA9szMzMy6r5XBO48A5wHnSVoP2Bd4O/AWUqD5H5KuI70uv7hDg33GU+qk2XLf0DyIp+5k8BHxDdJAokbvOweYA2key1bzZ2ZmZtYLWll5568iYmVE/DAiDgeeB+xDevW7HXASsDD3PWynUovkJjXOT65IZ2ZmZmZtMKHAslxErImI+RHx0YjYmjSP5JdIr47b6da8r9WHclre1+qDaX1sbPEyTp9/B2OL2zoHvpmZmTWglQnSGxIRo8AocHy7npGVpiaaKWlS+cjwPKhoD9L63gvanA/rsLHFyzj4zAWsWr2Gddee5HXHzczMuqywFstuiYg7SUsnbgt8uOL0bGAj4OyyOSxtQCxY+CCrVq9hTcCTq9ewYGE7F1cyMzOzeppusZS0sIFka4BHgFuA70dEU3M5SjoAOCD/WJp0/DWS5uY/L4mIY8su+RBpScdTJe2Tn7sbsBfpFXi9VXasD82YOoV1157Ek6vXsM7ak5gxdUq3s2RmZjbU1OxCOpIWkQLSF+ZDq4EHgSk8HajeQxo0szFpNPZlwAER8VSDz5hFWjWnlsURsW3FNS8mTf+zX87LvcDFwOyIWNrIc7tpZGQkRkdHu52NvjO2eBkLFj7IjKlT/BrczMysQySNRcTIs463EFhOBn4OPAF8GlgQEWskTQJeA5wIrEdadeb5pOUT9wP+OSK+NqFSDDAHlmZmZtYvagWWrfSx/Dxpap99IuK60mCZPCr8WlJAuSnw+Yi4HXgn8Gfg4JZzb2ZmZmY9r5XA8u3ADyNidbWTef3wS4F35J+Xk9bobnlJRTMzMzPrfa0EllOAdeukWSenK7mPNk5tZGZmZmbd10pguRA4MM8R+Sy5D+aBwB/LDr8A6PkBNGZmZmbWulYCyznAVsCvJR0saVtJG+T9e4Bfk0aMnwEgScCewM0F5dnMzMzMelDTr6cj4hRJ04EPAGdXSSJgTkSckn/eEvgOaSS5mZmZmQ2olvo9RsSHJJ0HHAbsRBol/ghwE2mVm6vK0v6FNC2RmZmZmQ2wVlbeWTciVkXENcA146TbNiIWTSRzZp3iidbNzMwmrpUWy/OAg8ZLkFfB+SUwtZVMmXXS2OJlHHzmAlatXsO6a0/i3CNmOLg0MzNrQSuDd94h6ZRaJyW9gBRUbtVyrsw6aMHCB1m1eg1rAp5cvYYFCx/sdpbMzMz6UiuB5deBj0g6tvKEpC1JQeV2wHsnmDezjpgxdQrrrj2JtQTrrD2JGVOn1L/IzMzMnqWVV+HHkKYT+qKkuyPiuwCSpgC/IK2wc3hEnF9cNs3aZ5dtNuPcI2a4j6WZmdkEtTLdUEg6mBREzpV0H2k0+DxgB+ADEVFtGiKznrXLNps5oDQzM5ugVqcbWinpbcC1wA9Iq/HsDHw8IuYUmD8zMzMz6xOt9LEEICKWAvsCy0lzWX66bFJ0MzMzMxsydVssJX2rTpLFwLrA9Iq0ERHvn0jmzMzMzKx/NPIq/LAG71WZLgAHlmZmZmZDopHAcru258LMzMzM+l7dwDIiFnciI2aDwEtDmpnZMGtpVLiZPZuXhjQzs2FXd1S4pA0m+pAi7mHW67w0pJmZDbtGphv6o6RjJK3X7M0l7SjpEuBZyz+aDRovDWlmZsOukVfh84CvD3EqUQAAHIZJREFUAZ+VdD7wPWBBRKyolljSVNL8lu8FdgXuAr5STHbNepeXhjQzs2GniKifSBoBTgT2yYeeAm4B7gWWAesDU4DpwOaAgL8ApwAnRcTKwnM+YEZGRmJ0dLTb2TAzMzOrS9JYRIxUHm9o8E5EjAIzJU0jzU25D2m1nb+tSPoA8H3gIuCiiHhyQrk2MzMzs77R1KjwiLgdOA5A0obAVqSWyhXA/RFxb+E5NDMzM7O+0PJ0QxGxHLg9b2ZmZmY25BoZFW5mZmZmVpcDS7MeNbZ4GafPv4Oxxcu6nRUzM7OGeOUdsx7kVXzMzKwfucXSrAd5FR8zM+tHDizNepBX8TEzs37kwLJBkl4v6YeS/iwpJB1WJc2nJd0g6RFJD0i6VNIOXciu9bnSKj6fmDndr8HNzKxvuI9l4zYGfgecnbdq9gS+AdxAWn3oBOAXkl4eEUs7kUkbHLtss5kDSjMz6ysOLBsUEZcBlwFImlsjzb7lP0s6BHgY2AO4tM1ZNDMzM+uqpl6FK3mTpC9LWiDpHkmrJD0s6XZJF0j6oKStJpoxSQdJOk3S1fnVckg6p841L5L0rZyvlZIWSTpZUreafZ5DqmPPF2NmZmYDr6EWy7x849HAPwFbk17zAjwB3A9sAEwFXgIcCJwi6VLgqxFxXYt5+wywI/AYcDfwN3Xy+BLgOmBL4BLgD8CuwDHAfpL2iIhOD609BbgZ+FWHn2tmZmbWcXVbLCUdTlq28UTSmuCzgTcBm0bEhhHxooiYQgpSXw68D7gIeDNwtaTzJW3dQt4+DmwPTAY+2ED6b5CCyqMj4oCIOC4i9gZOAqYDn68o1+dyK+h4254t5Lt0/68BrwUOjIinWr2PmZmZWb9opMXym8DFwBci4oZaiSIiSK2EfwDmSpoMHAocBxxGGsjSsIiYX/qzpPGSImkqMBNYBJxecfqzwFHAIZI+GRGP5+MnA+O+Wgf+1ESWy/NzEvAuYK+IWNjKPcyKMrZ4GQsWPsiMqVM8GMjMzNqqkcByJCJubPbGEfEIcJqk/wS2bfb6Ju2d9/MiYk1FPh6VdC0p8JwBXJ6PLwGWFJ0RSaeQgso9I+IPRd/frBlewcfMzDqp7qvwVoLKiuuf6ECANT3vb6tx/va8377VB0jaWNJOknYi1dvW+eety9KcDhwOvBtYJun5edu41eeaTYRX8DEzs04alAnSN8n7h2ucLx3fdALPGAFuytsGpL6mN/HMV/wfIo0Evxy4t2w7ttoNJR0laVTS6AMPPDCBrJlV5xV8zMysk1qax1LSRsBOwPrAnRGxqMhMtUGpk2a0eoOIuKLsPrXSjN8Z9Nnp5wBzAEZGRlrOm1ktpRV83MfSzMw6oenAUtKupOl8tiw79jCp9W4MuBEYi4jbq9+hLUotkpvUOD+5Ip3Z0PAKPmZm1imttFieBDwP+D5wC2n+ylcBbwD2IrcKSno0Iiby6rkZt+Z9rT6U0/K+Vh9MMzMzM5ugVgLLVwI/joiDyg/mASo7A7uQ+iO+auLZa1hpaqKZkiaVjwyX9BzSkoorgAUdzJOZmZnZUGklsHwM+H3lwYh4DLg6bx0VEXdKmkeaUujDwGllp2cDGwFnlM1haWZmZmYFayWwvBx4WdEZqSTpAOCA/OPz8/41kubmPy+JiPLR1h8iLel4qqR9SK/pdyO9nr8N+Nd259lskHmidTMzq6duYCnpXcCNEVHqn/h5YFTSjhHxmzbmbSfSyj3lpuYNYDFl0/jkVssR0vQ/+wH7k6b6ORWYHRFL25hXs4HmidbNzKwRjbRYngeEpMeBm0mjvn8O/FTSoRExrx0Zi4hZwKwmr7mLNEG5mRWo2kTrDizNzKxSI4HlsaRBOa8CdgdeW3buJ5IWAvNIUw2NAb+LiKeKzqiZdU9povUnV6/xROtmZlaTIhqfl1vShqRX1LuQAs1Xkfpbrs3Tk4+vAn5LmsvyQ4XmdoCNjIzE6Ohot7NhVpP7WJqZWYmksYgYedbxZgLLGjdeD9iRZwabOwBrR8RaE7r5EHFgaWZmZv2iVmDZ0pKO5SJiJXB93koPWwf424ne28zMzMz6x4QDy2oi4knSIB8zMzMzGxKT6iWQtMFEH1LEPczMzMyst9UNLIE/Sjom96VsiqQdJV1C2XyTZjbcxhYv4/T5dzC2eFm3s2JmZgVr5FX4POBrwGclnQ98D1gQESuqJZY0FdgXeC+wK3AX8JVismtm/cwTrZuZDba6gWVEvFfSqcCJwFF5e0rSLaSVbZYB6wNTgOnA5oCAv5CWUTwpD/AxsyHnidbNzAZbQ4N3ImIUmClpGvB+YB/SfJaVI78fAL4PXARclAfxmJkBnmjdzGzQNTUqPCJuB46Dv06WvhWppXIFcH9E3Ft4Ds1sYOyyzWace8QMT7RuZjagWp5uKCKWS3oC+FNE3FNgnsxsgO2yzWYOKM3MBlQjo8KfQcm/SnoQWATcJekRSRdL2r/wHJqZmZlZX2g6sASOB/4fsBlwK2ki9MeAvwMulfQjSZsUl0UzMzMz6wetBJZHkkaC7xwRL4+IV0fEC4HdgAuA/YGfSVq3wHyamT2D58M0M+s9rfSx3AqYExG/KT8YETcA75J0NXAa8AngixPPopnZM3k+TDOz3tRKi+Vy4NFaJyPidOBa0gTpZmaFqzYfppmZdV8rgeVNwMw6aa4Gtmvh3mZmdZXmw1xLeD5MM7Me0sqr8M8AV0o6MSKOr5FmC2BJ69kyM6vN82GamfWmpgPLiLhG0heA4yW9jrSO+E8i4gkASW8G/hE4udCcmpmV8XyYZma9p6UJ0iPiM3keyxOAC4E1kpYA6wKbApeSpiQyMzMzsyHRSh9LACLiJGB74F+AX5ECys0AAW8Flkm6QdIZkj4gadciMmxmZmZmvanlJR0B8trgXwW+Kmkt4BXALmXbK/MeIIC1JvI8M7OijS1e5r6aZmYFmVBgWS4ingJ+m7ezACRNAl4OjACvKupZZmZF8HyYZmbFavlVeCMiYk1E/C4i5kbE0e18lplZszwfpplZsdoaWJqZ9TLPh2lmVqzCXoWbmfUbz4dpZlYsB5ZmNtQ8H6aZWXH8KtzMzMzMCuHA0sysAGOLl3H6/DsYW7ys21kxM+savwo3M5sgT1tkZpa4xdLMbII8bZGZWeLA0sxsgjxtkZlZ4sCyAZJeL+mHkv4sKSQdVif98Tnd1zuURTProtK0RZ+YOd2vwc1sqLmPZWM2Bn4HnJ23miTNAI4kLW1pZkPC0xaZmbnFsiERcVlEHB8RFwJraqWTtAlwLvB+wENDzaxpHl1uZv2sJwNLSQdJOk3S1ZIeya+Vz6lzzYskfUvSPZJWSlok6WRJnWxCmANcGBG/7OAzzWxAlEaXf3XerRx85gIHl2bWd3r1VfhngB2Bx4C7gb8ZL7GklwDXAVsClwB/AHYFjgH2k7RHRLR1mKakI4GXAoe08zlmNriqjS7363Uz6yc92WIJfBzYHpgMfLCB9N8gBZVHR8QBEXFcROwNnARMBz5fnljS53Ir6Hjbno1mVtJ04ETg4IhY1eh1ZmblPLrczPqdIqLbeRhXDvDmA+dGxHuqnJ8K3AksAl4SEWvKzj0HuBcQsGVEPJ6Pbw5sXufRf4qI5VWe9xjwkYiYW3bsMOAs4KmypGsBQeqTuVFErBzvYSMjIzE6OlonS2Y26MYWL2PBwgeZMXWKWyvNrGdJGouIkcrjvfoqvBl75/288qASICIelXQtMBOYAVyejy8BlhSYh4uByqjwLOB2UkumWzHNrCFFji53kGpmnTYIgeX0vL+txvnbSYHl9uTAslmSNib1n4TUfWBrSTsBSyPiTxHxEPBQxTWP5/O/G+e+RwFHAWy99datZM3MrCovM2lm3dCrfSybsUneP1zjfOn4phN4xghwU942AGbnP58wgXsSEXMiYiQiRrbYYouJ3MrM7Bm8zKSZdcMgtFjWo7xvuTNpRFxRdp9Gr9mz1eeZmU1UaSDQk6vXeCCQmXXMIASWpRbJTWqcn1yRzsxs4JWWmSyij6X7appZowYhsLw177evcX5a3tfqg2lmNpCKGAjkvppm1oxB6GM5P+9nSnpGefJ0Q3sAK4AFnc6YmVm/c19NM2tG3weWEXEnMA/YFvhwxenZwEbA2aU5LM3MrHGetN3MmtGTE6RLOgA4IP/4fGBfYCFwdT62JCKOLUtfuaTjLcBuwF6kV+C7t3tJx4nyBOlm1quK6mPpvppmg6PfJkjfCTi04tjUvAEsBv4aWEbEnZJGSNP/7AfsT1px51RgdkQsbXuOzcwGlPtqmlmjejKwjIhZwKwmr7kLOLwd+TEzs4mp1lfTgaXZ4On7PpZmZtb7iuyrObZ4GafPv4OxxcsKzKGZFaEnWyzNzGywFDWvpl+pm/U2B5ZmZtYRRfTVLPKVugcTmRXPgaWZmfWNopaqdMunWXs4sDQzs75R1Ct1DyYyaw8HlmZm1leKeKVeVMunmT2TA0szMxs6RbV8mtkzObA0M7OhVETLp5k9k+exNDMzM7NCOLA0MzMzs0I4sDQzMzOzQjiwNDMzM7NCOLA0MzOztvL67sPDo8LNzMysbbzK0XBxi6WZmZm1TbVVjmxwObA0MzOztimtcrSW8CpHQ8Cvws3MzKxtvMrRcHFgaWZmZm3lVY6Gh1+Fm5mZmVkhHFiamZmZWSEcWJqZmZlZIRxYmpmZmQ2AXpiI3oN3zMzMzPpcr0xE7xZLMzMzsz7XKxPRO7A0MzMz63O9MhG9X4WbmZmZ9blemYjegaWZmZnZAOiFiej9KtzMzMzMCuHA0szMzMwK4cDSzMzMzArhwNLMzMzMCuHA0szMzMwK4cDSzMzMzArhwNLMzMzMCqGI6HYeDJD0ALC42/nosM2BJd3ORI9zHY3P9VOf62h8rp/xuX7qG9Y62iYitqg86MDSukbSaESMdDsfvcx1ND7XT32uo/G5fsbn+qnPdfRMfhVuZmZmZoVwYGlmZmZmhXBgad00p9sZ6AOuo/G5fupzHY3P9TM+1099rqMy7mNpZmZmZoVwi6WZmZmZFcKBpZmZmZkVwoGlNU3SFElHSPqBpDskrZD0sKRrJL1fUtXPlaTdJV0maamk5ZJ+K+ljktYa51lvkXRFvv9jkn4t6dD2la59JB0iKfJ2RI00TZdX0qGSrs/pH87Xv6U9pSiepNdJukjSvZJW5v08SftXSTt0nyFJ/yfXx935u7ZQ0gWSXlMj/UDVkaSDJJ0m6WpJj+Tvzzl1rulIHfTKd6+ZOpI0TdKnJP1S0l2SVkn6i6RLJO1V5zlNlVfSWrnef5s/u0vz38vuEy1zs1r5HFVc/82y398vrZGm6fJK2kDSbEm3SnpC0v2SvifpZa2UsydEhDdvTW3AB4AA7gHOBb4AfAt4KB+/kNx/t+yatwGrgceAbwJfAf6Q019Q4zkfyeeXAKcDJwF35WP/3u16aLLOXpzr59Gc/yOKKC/w7/n8XTn96cCD+dhHul3uBurlMzmvDwBnASeSOsLfAHx52D9DwJfK8n8m8MX8/VoFrAHeM+h1BNyc8/IocEv+8znjpO9IHfTSd6+ZOgK+m8//HjiD9Pv7+7nOAji6iPICAi7I5/+Q/x6+mf9eVgNv69U6qnLtW8uuDeClRZQXWA+4Jl9zQ/6+nwc8CTwO7NbJOiqsrrudAW/9twF75y/apIrjzwf+lL8kB5YdnwzcD6wERsqOrw9cl9O/q+Je2wJP5F9c25Yd3wy4I1/zmm7XRYP1JeAXwJ35l82zAstWygvsno/fAWxWca8H8/22bVe5CqiXd+b8/xx4TpXz6wzzZyh/n54C7gO2rDi3V87/wkGvo1zWafl7tCfjB00dqYNe++41WUeHATtXOf4G0n9YVgIvmGh5gXfna64F1i87/ur8jPup8r3vhTqquG6L/B38LnAFtQPLpssLfDpfcwFl/56S/nNUCv4ntVLebm5dz4C3wdqA4/MX4rSyY+/Lx/6rSvq987krK46fkI/PrnJNzfv14gYcQ2pdej0wi+qBZdPlBc7Oxw+vck3N+/XCRuqGs5D0v/ItGkg/dJ8hYLecx0tqnH8EeHSY6qheQNCpOujl7169Oqpz7TwqGgZaLS9wVT6+V5Vrat6v1+oI+AEpsJzC+IFlU+UlBbiL8/Htmrlfr2/uY2lFezLvV5cd2zvvf1ol/VXAcmB3Ses1eM1PKtL0rNxP5ovAKRFx1ThJWylvP9fR7sB2wGXAstyP8FOSjqnRd3AYP0O3k1qQdpW0efkJSa8HnkNqCS8Zxjqq1Kk6GLR6K6n2+xuaLG+u391J9X11I9f0IkmHAQcAH4iIB8dJ10p5XwJsDdwWEX9s8Jq+4MDSCiNpbeC9+cfyX0DT8/62ymsiYjXwR2BtYGqD19xLaul6kaQNJ5jttsn18W1S94Dj6yRvqrySNgK2Ah7L5yvdnvfbt5D1Tnh13v8FuBH4ESkAPxm4TtKVkrYoSz90n6GIWAp8Cnge8D+S5kj6gqTvkVqWfg78U9klQ1dHVbS9Dgbgu1eVpG2AfUjB0VVlx1sp70uBtUhdNSqD1FrX9JRcH6eQWjUvrpO8lfLW/NyNc01fcGBpRfoisANwWUT8rOz4Jnn/cI3rSsc3beGaTWqc7wX/F9gZOCwiVtRJ22x5W6nTXrJl3n8A2AB4I6kFbgfgZ6RuAxeUpR/Kz1BEnAy8gxQMHQkcR+qbehcwNyLuL0s+lHVUoRN10O/fvWfJLW7nkgaTzIqIZWWn21mnPVlHSjOb/Bdp4M3RDVwydHU0HgeWVghJRwOfJI2GO6TZy/M+2nxNx0jaldRK+dWI+FURt8z7Zsvbk/VD+t89pHIdFBGXR8RjEfF74O3A3cAbak2pU8XAfYYAJP0LaRT4XNKrs42AXUj9U8+V9OVmbpf3A1VHTepkHfRFneUpmL4N7AGcTxr93YpB+lx9nDSY6ciKILtVQ/Xdc2BpEybpw6RXBv9D6mi8tCJJvVaPyRXpmrnmkSay2hFlr8BvA/6twcuaLW+99PX+N9xtpV/WCyPiN+UncutuqcV717wfqs8QgKQ9SdOP/DAiPhERCyNieUTcSAq+/wx8UlLpte7Q1VEVnaiDfv/u/VUOKs8htYJ/jzR9VWUg00p5W/l76AmSpgGfB86KiMsavKydn7ueq6N6HFjahEj6GPB14HekoPK+Ksluzftn9RXJQdh2pM7iCxu85gWklpu7I2J567lvm41J+X4Z8ETZpLoBfDan+c987OT8c1PljYjHSYHFxvl8pWl5X6v/TreVyvtQjfOlwHODivTD8hkCKE08Pb/yRM7z9aTf4Tvnw8NYR5XaXgcD8N0D/lof3wHeRZo78R+r9Q9ssbx3kKbKmpqf08g1veIVpC4Bh5f/7s6/v9+Q09yejx2Qf26lvDU/d+Nc0xccWFrLJH2KNFHuzaSg8v4aSX+Z9/tVOfd6YEPguohY2eA1b65I02tWkibGrbbdlNNck38uvSZvpbz9XEdXkf6BnyZp3Srnd8j7RXk/bJ8hSP+4QZpHr5rS8VV5P4x1VKlTddDX9Za/cxeSWirPBg6JiKfGuaSp8ub6vY5U369r5Joesojav79LDScX5J8XQcvlvZM0sHN7Sds1eE1/6PZ8R976cyO94g1gFHhunbSTSSurNDNp8Xb0+MTNLdbbLKrPY9l0eemxSZpbqItzcv4/V3H8TaR5Px8CNh3WzxDw9zmP9wFbVZx7c66jFcCUYakjGpsgve110MvfvQbqaD3gxznNmTQwAXcr5aWxCcMn92IdjXPdFUxsgvTJFdcM5ATpyoUwa5jS+rlzSU3/p1G9D8iiiJhbds0BpP8hP0FawWAp8HekKRcuBP4+Kj6Mkj4KnEr6xXU+qWXmIOBFpEExxxZZrk6QNIv0OvzIiDiz4lzT5ZX0VeATpMEuFwLrAv9Amsz3oxHx9bYVZoIkbUn6JfxS0txv1wPbkPoPBunV3AVl6YfqM5RHpv6MNGL+UZ6eqPllpNfkAj4WEaeUXTNwdZTLVHrl+HxgX9Kr7NJ8gUvK89ipOuil714zdSTpLNLqO0uAb1B9cMgVEXFFxTOaKq8kkfptHkQa1HlpTvsPpED/wIi4pNUyN6vZz1GNe1xBeh0+LSLuqDjXdHnzaPxfkgL3UeBy0tyW7yR9BveOiF+3UNzu6nZk663/Np5udRtvu6LKdXuQJ8QmtbT8N2n03VrjPOutwJWkf1gfJ62nemi366CAunvWWuGtlhc4NKd7PF93JfCWbpe1wfp4LvA10vyCq0j/sF8CzKiRfqg+Q8A6wMeABaTBI6tJLR8/AmYOQx018PtmUbfqoFe+e83UEU+3uo23zSqivKRpsj6e639F/vu4DNi9l+tonHuU6u5ZLZatlpfUj3w2ad7KlaQW9wuAl3e6jora3GJpZmZmZoXw4B0zMzMzK4QDSzMzMzMrhANLMzMzMyuEA0szMzMzK4QDSzMzMzMrhANLMzMzMyuEA0szMzMzK4QDSzOzNpK0raSQNLfg+54g6QlJL273sxrMz4aS7pP07U4/28x6hwNLM7M+k4PJY4E5EXFXt/MDEBHLgS8AB0vatdv5MbPucGBpZtZ//g1YD/hKtzNS4QzgIeBz3c6ImXWHA0szsz4iaRPgYODyXmmtLImIJ4DzgTdKmtbt/JhZ5zmwNDPrkNwH8ruSluT+kaOS3tLkbd4NbEgK4Bp97iRJp+b+l9+XtH4+Pisf21PSuyWNSVou6R5JX5O0Xk63t6QrJD0iaZmkb0uaUuNx3wUEvK/JcpnZAHBgaWbWGdsA1wPbAt8mBYY7AJdI2quJ+7wx769pJHEOIr8HfBQ4HTgotyyW+yjwTeBW4P8DDwIfB86Q9HbgJ8BSYA5wC/Ae4Jwaj7weeBJ4U4PlMbMBsna3M2BmNiT2BGZFxOzSAUnnAT8F/hmY3+B9Xgs8CtxWL6Gk5wKXAHsAx0XEl2okfSOwS0Tckq9bD7gROAR4KzAzIq7M5yYBPwP2k7RTRNxcfqOIWCHp98DOkp4TEY82WC4zGwBusTQz64zFVAxqiYifAX8CGhpFLWld4HnAfRERddJuA1wL7AYcMk5QCXBqKajM+VpJalGdBPy4FFTmc2t4urVyxxr3uy9fu9X4JTKzQePA0sysM26OiKeqHL8L2KzBe5T6NS6rk2468CvghcCbI+LcOulHqxy7J+/Hqpz7c96/qMb9lub95nWea2YDxoGlmVlnPFTj+Goa/128Iu/Xr5Nue+AFwELSK+16Hq6Rr3rn1qlxvw3yfkWN82Y2oBxYmpn1iYh4CFjF0y2XtVwKHA/sBFwuqdMth6X83d/h55pZlzmwNDPrL/8NvEDS5PESRcQXSCO7dwbmS3peJzKXTSeNLL+7g880sx7gwNLMrL9cQfrdXXfAT0ScDHwQeAVwpaQXtjdrIGk70gCjK+oNMDKzwePA0sysv1yU9/s2kjgi/oM0Wfk04CpJW7crY9nMvL9o3FRmNpAcWJqZ9ZGI+BVwE3CwpLUavGYuaVLzbUjB5dT25ZBDgQdwYGk2lOQ3FWZm/UXSu4HzgHdExA+6nZ8SSa8EfgP8W0R8rl56Mxs8DizNzPqMJJHmqdwA2KlX+jJKuhjYBdg+IjzVkNkQ8qtwM7M+kwPJo4AfkCZB7zpJG5Je0R/ioNJseLnF0szMzMwK4RZLMzMzMyuEA0szMzMzK4QDSzMzMzMrhANLMzMzMyuEA0szMzMzK4QDSzMzMzMrhANLMzMzMyvE/wLBJG0j6IhT2AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams.update({'font.size': 20})\n",
    "dragtable = pd.read_excel(\"SMAD_tables.xlsx\", sheet_name=\"Drag\")\n",
    "H = dragtable[\"Scale Height\"][2:].values\n",
    "h = dragtable[\"Altitude [m]\"].values[2:].astype(float)\n",
    "rho = dragtable[\"Unnamed: 3\"].values[2:].astype(float)\n",
    "f, (ax1, ax2) = plt.subplots(2,1, figsize=(10,10))\n",
    "ax1.plot(h, H, '.')\n",
    "ax1.set_xlabel(\"h (km)\")\n",
    "ax1.set_ylabel(\"H(h) (km)\")\n",
    "ax2.semilogy(h, rho, '.')\n",
    "ax2.set_xlabel(\"h (km)\")\n",
    "ax2.set_ylabel(r'$\\rho(h)$ (kg/m3)')\n",
    "plt.savefig('NLRMSISE00.eps')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I believe that the atmospheric models are based on the NRLMSISE-00 standard by NOAA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the propulsion lifetime, it is the total momentum available from propulsion, $l=m_P \\cdot I_{sp} \\cdot G$ divided by the drag force:\n",
    "\n",
    "$F_D = 0.5 \\rho \\cdot C_D \\cdot A \\cdot v^2 = 0.5 \\rho \\cdot C_D \\cdot A \\cdot \\frac{\\mu}{a}$\n",
    "\n",
    "\n",
    "$L_p = \\frac{l}{F_D} = \\frac{m_P \\cdot I_{sp}\\cdot G \\cdot a}{0.5 C_D \\cdot A_D \\cdot \\rho \\cdot \\mu} $ \n",
    "\n",
    "Turns out we actually need to add a constant to make $L_p$ realistic:\n",
    "\n",
    "$L_p = \\kappa \\frac{m_P \\cdot I_{sp}\\cdot G \\cdot a}{0.5 C_D \\cdot A_D \\cdot \\rho \\cdot \\mu} $ where $\\kappa$ is a constant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n",
      "The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAAFPCAYAAACRa82TAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeVwV9f7H8deXHRFRBENRBDfcQyUzyz2XMs1scym1LJf6qWW2l1ravbZo4m0zyzRbzMy8eTWXXNIsV0DNBRdSwjRUEJBFtu/vjwECBeTggTkcPs/H4zzozMyZeUPRh+/Md1Faa4QQQghhDgezAwghhBBVmRRiIYQQwkRSiIUQQggTSSEWQgghTCSFWAghhDCRFGIhhBDCRE5mB7AWHx8fHRgYaHYMIYQQVdjevXvPa619LfmM3RTiwMBA9uzZY3YMIYQQVZhS6pSln5Fb00IIIYSJpBALIYQQJpJCLIQQQphICrEQQghhIinEQgghhImkEAshhBAmkkIshBBCmEgKsRBCiCppdfRq+izvQ9vFbemzvA+ro1ebksNuJvSwxOXLl4mPjyc5OZns7Gyz4whR6Tg6OuLp6Ym3tzeurq5mxxHCYqujVzP91+mkZ6cDcCblDNN/nQ5A/0b9KzRLlSvEly9fJiYmhlq1ahEYGIizszNKKbNjCVFpaK3JzMwkKSmJmJgYAgICpBiLSicsPCy/COdJz04nLDyswgtxlbs1HR8fT61atfDx8cHFxUWKsBAWUkrh4uKCj48PtWrVIj4+3uxIQljsbMpZi7aXpypXiJOTk6lRo4bZMYSwCzVq1CA5OdnsGEJYzM/Dz6Lt5anKFeLs7GycnZ3NjiGEXXB2dpZ+FqJSmtR+Em6OboW2uTm6Man9pArPUuUKMSC3o4WwEvldEpVRVk4WcalxPHfTc9T1qItCUdejLtM7T6/w58NQBTtrCSGEqNpWnVjFnL1z+E/P/7D+vvVmx6maLWIhhBBVU2ZOJvP3z6dV7VZ0q9/N7DiAFGIhhBBVyKoTqzh96TRPhDxhM49WpBALIYSoEjKzM5m/bz5tfNrQxb+L2XHySSGu4nr27IlSim3btpkdRQghylViRiJNazVl/I3jbaY1DNJZq0rTWhMeHo6DgwPt2rUzO44QQpQrH3cf3uv1ntkxriIt4irs2LFjJCYmEhwcTPXq1c2OI4QQ5Wb32d38mfyn2TGKJIW4Ctu7dy8AHTp0MDmJEEKUn4zsDF765SVe+eUVs6MUSW5Nl5OVEad5e10Uf11Mo15Nd57tG8ygdv5mxypkz549AISGhpqcRAghys/3x77nbMpZXuv8mtlRiiSFuBysjDjNiysOkJZpTP13+mIaL644AGBTxVgKsRDC3mVkZ7DgwALa1WnHLXVvMTtOkaQQX+HB+b9dte2utnV5+JZA0jKyGfXZrqv239ehPveHNiA+JYPxX+wlIuYiGdk5hY5Jy8zm3z8e5utdMVd9/vEujbi95Q2cOHeJl3ILdkETejbltqY+HPwrkVb1vK7ju/tHTk4OERERODo6EhISYpVzCiGErfnu2Hf8nfo3M2+baVM9pQuSZ8Tl4MoinCcu6XIFJyne0aNHSU5Opnnz5nh4eJgdRwghykVCegKd6nbiZr+bzY5SLGkRX+GbscXfunB3cSxxv7eHC9+MvYVbZ23i9MW0q/bXq+le4ucb+1Yvcb+1WsMgt6WFEFXDEyFPkKNzbLY1DNIiLhfP9g3G3dmx0DZ3Z0ee7RtsUqKr5RXi0vSYjo2NRSnFf//7X/r164eHhweNGzdm8+bN5R1TCCHKJD0rnYi4CAAclG2XOttOV0kNaufPvwe3wb+mOwrwr+nOvwe3samOWnlDl0rTIo6MjARgzpw5TJkyhX379tG6dWsmT55crhmFEKKsvj36LSN+HMGR+CNmR7kmuTVdTga187epwluQpR219u3bR40aNfjmm2/w8/MD4L777uPFF18s76hCCGGxtKw0Pj3wKR39OtLcu7nZca5JCnEVdPjwYVJSUvDw8GDChAlFHuPj48OsWbMAo0Xcv3///CIMcPz4cZo0aVIheYUQwhLLopZxIf0Cs0Nmmx2lVKQQV0F5t6VTUlL49NNPizymb9+++f+8b9++qwp2RESEDHsSQtic1MxUFv6+kJvr3kyHGyrHrIHyjLgKGjFiBFrrEl9r164FjGJ94sSJqxaFkEIshLBFxy8eJ1tn82TIk2ZHKTVpEYsS7d+/H6BQ0b1w4QKxsbFSiIUQNqetb1s23LcBdyd3s6OUmrSIRYn27dtH06ZNC63OFBERgbOzMy1btjQxmRBCFPZH4h9k52RXqiIMUojFNYwbN44jRwp3/7/99tvJyMjAxcXFpFRCCFFYamYqI38cyb92/svsKBaTQiyEEKLS++rIVyRcTmBQk0FmR7GYFGIhhBCVWkpmCosOLqKLfxfa+LYxO47FpBALIYSo1L46/BWJlxMZf+N4s6OUiRRiIYQQlZbWmi2xW+hav2ulbA2DDF8SQghRiSmlWNxvMUkZSWZHKTObaRErpforpdYrpWKVUmlKqWil1LdKqeLXBRRCCFFlpWamkpqZipODE95u3mbHKTObKMRKqTeB/wHtgbVAGBAO3A1sV0o9ZGI8IYQQNmjxocXcueJOEi8nmh3luph+a1op5QdMAf4G2mqt4wrs6wFsAl4HvjAnoRBCCFuTlJHEkoNLuMnvJrxcvcyOc11soUXcECPHzoJFGEBrvRlIBnzNCCaEEMI2fXHoC5IzkxkfUjl7ShdkC4X4GJABdFRK+RTcoZTqCngCP5kRTAghhO1JvJzIkkNL6BXQq1KsN3wtpt+a1lrHK6WeB+YAh5RSK4ELQGNgILABGGtiRCGEEDZkw6kNXMq8VGnHDV/J9EIMoLWeq5Q6CSwEHi+w6ziw6Mpb1nmUUmOAMQABAQHlHVMIIYQNuLfpvbT1bUuzWs3MjmIVtnBrGqXUc8ByYBFGS9gD6ABEA18qpd4q6nNa64+11qFa61BfX3mMLIQQ9i49Kx2llN0UYbCBQqyU6g68CfygtZ6stY7WWqdqrcOBe4DTwDNKqUZm5hTCXk2fPh2lVKGXn5+f2bGEuMrF9Iv0/a4vP5z4wewoVmULt6bvyv26+codWutUpdQujILcDqOFLISwsuDgYLZs2ZL/3tHR0bwwQhRj8aHFJKQn0Kp2K7OjWJUtFGLX3K/F3VvO255RAVmEqJKcnJykFSxsWkJ6Al8e/pJ+gf1oXLOx2XGsyvRb08C23K9jlFL+BXcope4AbgXSgV8rOlhV0LNnT5RSbNu27doH27Hly5czYcIEunTpQo0aNVBK8dBD5TOhW1muFRsby6OPPkq9evVwdXUlMDCQp556ioSEBKtkio6Oxt/fn6CgIIYMGUJ0tNx8Erbls4OfkZ6Vzrgbx5kdxepsoUW8HGOc8O3AYaXU98BZoAXGbWsFvKC1vmBeRPuktSY8PBwHBwfatWtndhxTzZw5k3379lG9enXq16/PkSNHbOZaJ06coHPnzsTFxXH33XfTvHlzdu3aRVhYGGvXrmX79u3Url27zHluvvlmFi1aRPPmzYmLi2PmzJl07tyZgwcPXtd5hbCW1MxUlkUt446gO2hU0w67C2mtTX8BzsBTwA4gCcgC4jDmn+5TmnN06NBBl8ahQ4dKdVxVEBUVpQHdokULs6OYbtOmTfro0aM6JydHb968WQN6+PDhNnGtPn36aEDPmzev0Pann35aA3rs2LGFtr/88ssaKPG1efPmYq+XnJysfX199ezZs0v1/cjvlKgIJy6e0LHJsWbHuCZgj7awBtrCrWm01pla67la605a6xpaayetdR2t9V1a6/Vm57NXe/fuBaBDhw4Vet3Dhw+jlGLw4MGFtu/cuRN/f38aNmxIZGRkhWbq0aMHTZs2RSllU9eKjo5m/fr1BAYG8uSTTxba99prr+Hh4cGSJUtISUnJ3/7UU09x+PDhEl8dO3Ys9prVq1enVatWHDt2rOzfpBBWkqNzAGjk1Qj/6v7XOLpysoVb08Ike/bsASA0NLRCrxseHg5Q6Hb44sWLGTt2LDfffDPLly9HxoUbNm3aBECfPn1wcCj8d7Onpye33nor69evZ8eOHfTq1QsAHx8ffHx8rjpXaaWnp3PkyBF69OhR9uBCWMm7e9/lTMoZ3ur6Fg7KJtqOVieFuLzsXwYbX4fEWPCqD72mQtsHzE5ViNmFuH379mRnZzNlyhTmzp3L+PHjCQsLw9nZucTPz507l4sXL5b6eiEhIQwaNOi6Ml/LoUOHyMjIIDg4GHd3d6udNyoqCoBmzYqevKBp06asX7+eo0eP5hdiS02ZMoUBAwYQEBBAXFwcM2bMICUlhZEjR5Y5txDWcD7tPEuPLKVPYB+7LcIghbh87F8GqyZCZprxPvFP4z3YTDHOyckhIiICR0dHQkJCKvTaeYU4KCiIfv36sXXrVhYsWMBjjz1Wqs/PnTuXU6dOlfp6I0eOLPdCfOedd3Lq1CkiIiKs+vNMTDTWWfXyKnqZt7ztlvxhcqXY2FiGDh3K+fPn8fX1pVOnTuzYsYOGDRuW+ZxCWMPC3xeSmZPJ2Lb2vdyAFOIrfdb/6m2tBkHHxyEjFb68/+r9IcOg3XBIuQDLRkDsbsi+XPiYzDTYMBX2Lr76853/D4LvgPPHYNVTV+/vOgUa94Az+6Fu27J9X1c4evQoycnJtGrVCg8PD6ucs7QiIyNxcXFh4MCBnDhxgo4dO5a6CAOcPHmy/MJVMkbfEK7r2fbSpUutFUcIqzmXeo5lUcu4q9FdBNSw77UEpBCXhyuLcJ7ks+BtGwPRzbotHR0dnd96q1OnDtWqVWPXrl2sWrWKAQMGVGgWa1q4cCGpqak0amTdoRV5Ld68lvGVkpKSCh0nhL1YdHARWTlZdt8aBinEV3tkdfH7XKqVvN+jtrH/3dbG7egredUv+fM+TUveb6XWMPxTiCu6x3Tebel77rmHpUuXsm7dOgYOHMhLL71E//79r+qQVBRbfEbcs2fPcjlvcHAwYNzBKEpez+biniELUVmNbjOaG31vpEGNBmZHKXdSiMtDr6mFnxEDOLsb221E3tCl0rSI854dzpkzh6+++ooDBw4QEBDAe++9R58+fSy6bl4hHj58OC4uLgwYMCD/meQXX3zBiBEjrnkOW3xGXF7yei6vX7+enJycQn+oJCcns337dtzd3enUqZNZEYUoF95u3vQJtOz/L5WV/XZDM1PbB2DAPPBqACjj64B5lbajVkREBADvvfces2bNYv/+/bRt25Zhw4aRlvbPHxujRo1CKcWiRYuKPVdRQ5dmzJgBwLRp08jIuPaU4idPnrRosHxJeWxd48aN6dOnDydPnuT9998vtG/atGmkpKQwYsSICn/OL0R5OZtylpE/juRIfPnNbmdrpEVcXto+YDOF90qHDx8mJSUFDw8PJkyYUOQxPj4+zJo1CzA6Vzk6OrJmzZr8W6VvvvkmTZo04ciRI/lFNSfHGHjv5FT8f1YRERF4eXkRFBSUv+3222+ne/fubNmyhQ8//JBJkyZZ5fu0xMqVK1m5ciUAZ8+eBeC3335j1KhRgPHzeOedd0y51gcffEDnzp2ZOHEiGzdupEWLFuzcuZPNmzfTrFkz3njjDavkEsIWfHLgE/af24+ni6fZUSqOpVNx2epLprgsvcWLF19zCsS+ffvmHz9s2DA9aNCgQueIiYnRgA4PD8/fFhISoj09PXV8fHyR142NjdWA7tat21X7tm/frgHt6+urk5KSrPONWmDatGkl/jwaNmxo6rViYmL0qFGjtJ+fn3Z2dtYBAQF64sSJ+sKFC1bLVVbyOyWs5cylM7rd5+309F+nmx2lzCjDFJdK5w5/sIRSygXwAdK01tZZ/uU6hYaG6rwOSCU5fPgwLVq0qIBE9qNVq1Y8+OCDTJ36zzPuFStWMHz4cC5cuEC1atW4ePEitWvX5plnnuGtt94yMa2oaPI7Jaxlxm8zWHF8BWvuWUPd6nXNjlMmSqm9WmuLhqOU6ta0UsoTGAL0BrpSYO1gpVQWsB/YBKzQWu+0JICwbenp6URFReXfds4zd+5chgwZQrVq1QDYtm0bzs7OTJ482YyYQohK7sylM6w4voLBTQZX2iJcViUW4tz1gV8FhgN5vUEuAlFAPOAO1AZCgA7AFKVUJPCO1vrr8gotKs6BAwcA+Prrr+nVqxd16tRh5syZHD9+nGXLluUfN2DAANLT082KKYSo5LzdvXk29Fl6BpTPUEBbVmyvaaXUaxgFdzSwDRgFNNVae2utW2qtb9Nad9BaBwJeQE/gLYzW8pdKqR1KKesNfBWmiIyMpGnTprz22msMHTqUdu3akZyczO7du/Hz8zM7nhDCTrg6ujKsxTD8PKre/1dKahE/C3wEvK21PlPSSbTWqcAWYItS6iXgLmA6MAjjtrWopCIjI2nTpg1DhgxhyJAhZscRQtih9yLew7+6P/c0vcfsKKYoaRxxE6315GsV4SvldhxbpbXuAHx8ffGE2SIjI2nbVm5sCCHKR2xyLJ8e+LRKjRu+UrGFWGv91/WeXGt99nrPIcyjtebAgQNSiIUQ5ebj/R/joBwY3Wa02VFMY5UJPZRStYAMrXWKNc4nbINSKn9RASGEsLY/k/7khxM/MKT5EOpUq2N2HNOUeopLpVQvpdRbuUU3b1sdpdTPwHkgXik1pzxCCiGEsD/z98/HycGJ0a2rbmsYLGsRTwBaa62fK7DtHaALcAzwBCYppXZorZcVdQIhhBAiT++GvWlRuwW+1XyvfbAds6QQ3wj8nPdGKeUO3Ads0Fr3zZ304wAwDpBCLIQQokTdGnQzO4JNsGT1pTpAwQ5cNwNuwCIArXUy8D8g2FrhhBBC2J+TiSeZFz6P5Ixks6PYBEsK8WWMmbTydMGYpH5rgW1JgLcVcgkhhLBT8/fPZ8mhJWRkX3vZ06rAkkL8B8bsWXnuBY5prU8X2NYAo+OWEEIIcZU/Ev9gzR9rGNJ8CLXda5sdxyZYUogXA22UUjuVUtuANsBXVxzTHmNaTCGEEOIqH+37CFdHV0a1GmV2FJthSWetD4FOwIOAAlYBb+btVEp1BFoAstiDEEKIq0RfjObHP35kVOtR0houoNSFWGudCQxTSo0z3uorn7JHA+2Ak9aLJ4QQwl4opbi94e3SGr5CqQuxUmohcEBr/W5R+7XW55Hnw0IIIYoR5BXEnO4y79OVLHlGPAxjCJMQQghhkaVHlhKTFGN2DJtkSSE+iRRiIYQQFjqacJQ3dr7ByuMrzY5ikywpxF8BdxSca1oIIYS4lo/2fUR15+qMbDXS7Cg2yZJC/G9gD7BZKXWXUuqGcsokhKhA06dPRylV6OXn52d2LGEnouKj2HBqA8NbDMfL1cvsODbJkuFL6blfFfBfMHrAFUFrra2yvKIQomIEBwezZcuW/PeOjo7mhRF25cN9H+Lp7MnDLR82O4rNsqRgbsOY0lIIYWecnJykFSysLkfn4OvuyyOtH5HWcAlKfWtaa91da92jNK/yDCysq2fPniil2LZtm9lRTHPhwgU++eQT7rnnHpo0aYK7uzteXl7cdtttfPrpp+Tk5Fj1esuXL2fChAl06dKFGjVqoJTioYceKvEzsbGxPProo9SrVw9XV1cCAwN56qmnSEhIsEqm6Oho/P39CQoKYsiQIURHR1vlvKJqc1AOvNzpZR5v+7jZUWya3EKuwrTWhIeH4+DgQLt27cyOY5pvv/2W8ePHU7duXXr06EFAQAB///03K1as4LHHHuPHH3/k22+/Le5RjMVmzpzJvn37qF69OvXr1+fIkSMlHn/ixAk6d+5MXFwcd999N82bN2fXrl2EhYWxdu1atm/fTu3aZZ+l6Oabb2bRokU0b96cuLg4Zs6cSefOnTl48OB1nVdUbScuniDxciLtb2hvdhTbp7W2i1eHDh10aRw6dKhUx1UFUVFRGtAtWrQwO4qpNm7cqH/44QednZ1daPuZM2d0gwYNNKCXL19utett2rRJHz16VOfk5OjNmzdrQA8fPrzY4/v06aMBPW/evELbn376aQ3osWPHFtr+8ssva4zHSMW+Nm/eXOz1kpOTta+vr549e3apvh/5nRJF+b+f/k93/qqzTslIMTtKhQL2aAvrlyW9pgFQStVVSj2klHpeKTW1iNer1vkTQZS3vXv3AtChQ4cKve7hw4dRSjF48OBC23fu3Im/vz8NGzYkMjKywvL07NmTAQMG4OBQ+NfBz8+PcePGARTqyHS9evToQdOmTUvVwo6Ojmb9+vUEBgby5JNPFtr32muv4eHhwZIlS0hJScnf/tRTT3H48OESXx07diz2mtWrV6dVq1YcO3as7N+kqNIOnj/IltgtjGg5gmrO1cyOY/MsujWtlHoNeOGKzyn+6cSV988zyhJGKdUFeArojLGucTxwAJirtV5TlnOK4u3ZsweA0NDQCr1ueHg4QKHb4YsXL2bs2LHcfPPNLF++HF9f3wrNVBxnZ2fA6Mxkhk2bNgHQp0+fq/5Q8PT05NZbb2X9+vXs2LGDXr16AeDj44OPj0+Zr5mens6RI0fo0UO6e4iy+WDfB3i5ejG8xXCzo1QKlsw1PRx4FdgEvA98BywC1gPdgdHAt8D8sgRRSr2CUcDPA/8DzgA+GAtJdAekEFuZ2YW4ffv2ZGdnM2XKFObOncv48eMJCwvLL37FmTt3LhcvXiz19UJCQhg0aJDFObOysvj8888B6NevX4nHHjp0iIyMDIKDg3F3d7f4WsWJijJWFW3WrFmR+5s2bcr69es5evRofiG21JQpUxgwYAABAQHExcUxY8YMUlJSGDlSJl8Qljtw7gBbY7cysd1EqrtUNztOpWDJn/njgVign9Y6K/e22kmt9VJgqVLqe2A1ZVgGUSl1P0YR/gkYrK9Y2UkpVfL/mW3Q6ujVhIWHcTblLH4efkxqP4n+jfqbHStfTk4OERERODo6EhISUqHXzivEQUFB9OvXj61bt7JgwQIee+yxUn1+7ty5nDp1qtTXGzlyZJkK8QsvvMDvv//OnXfeSd++fUs89s477+TUqVNERERY9eeZmJgIgJdX0UM/8rZb8ofJlWJjYxk6dCjnz5/H19eXTp06sWPHDho2bFjmc4qq60zKGQJrBDKsxTCzo1QalhTiNsDXWuusAtvyR/1rrdcppdYBz2KsVVwqSikHjHWNU4FhVxbh3HNnWpDTdKujVzP91+mkZxtzoJxJOcP0X6cD2EwxPnr0KMnJybRq1QoPD48KvXZkZCQuLi4MHDiQEydO0LFjx1IXYYCTJ0+WX7hc8+bNY/bs2TRv3pwlS5aU+/XKyugbUuzkOqWydOlSa8URgj6Bfbi94e04KIu7IFVZlhRiZ+BCgfdpwJV/pv8OjLMwQ2cgCFgOJCil+gOtMWby2qW1/s3C812XR9Y+ctW2voF9GdJ8CGlZaTzx0xNX7b+7yd0MajKIhPQEJm+ZzP5z+8nIySh0THp2OnP2zmH50eVXfX5kq5F0b9CdPxL/4PXfXr9q/5i2Y7il3i0ciT9Cc+/m1/Hd/cOs29LR0dH5rbc6depQrVo1du3axapVqxgwYECFZinO+++/z6RJk2jZsiUbN27E29v7mp9ZuHAhqampNGrUyKpZ8lq8eS3jKyUlJRU6Tggz7Tizg5tuuAlHB5mZzRKWFOIzQN0C72OAtlcc4w9kYZmbcr/+DYRjtLzzKaW2Avdprc9ZeF7TXFmE85xLPUeAZ0AFpylaXiGu6B7Tebel77nnHpYuXcq6desYOHAgL730Ev3797+qQ1JRyvMZ8dy5c3n66adp3bo1GzdupE6d0i041rNnz1LnsURwcDBg3MEoSl7P5uKeIQtRUSLjInl8/eO82ulVHgh+wOw4lUtpxzlhtFh3F3j/PpAJPAx4AP0xWsk/WTJ+CmMxCY1RwI8BvYDqQCtgbe6+LcV8dgzGQhR7AgICSjXGqyLGPPb+trduvaj1Va/e3/Yu92uX1m233aYB/euvv17z2HPnzmlAz5kzR4eGhmpXV1fdtGlTvW7dOouv++KLL141LrdTp04a0IsXLy7VORo2bHjNcbIFXyNHjizVeWfNmqUBHRISos+dO2fx91YW1xpHfPz4cQ3owMDAq8Y5JyUlaQ8PD+3u7q4vXbpUEXGLJOOIhdZaP77ucd11adcqN274SpTzOOL/Aa2UUkG572cBiRg9p5OAHzCGL71iwTnhn+fMCqPlu1FrfUlrfRC4B6ODWDel1C1XflBr/bHWOlRrHWorw10AJrWfhJujW6Ftbo5uTGo/yaREhVnaUSsiIgKA9957j1mzZrF//37atm3LsGHDSEtLyz9u1KhRKKVYtGhRsecqaujSjBnGaLdp06aRkVH03YSCTp48adF/5CXlKZjhhRdeoEOHDmzcuPG6hv9YU+PGjenTpw8nT57k/fffL7Rv2rRppKSkMGLEiAp/zi9EQRFxEfx25jceafWIjBsug1LfmtZaL8Iounnv/1RK3QQ8AzQGTgIfaK0PWJghb7LcaK31viuumZbbAWw00BGo0OfFZZXXIctWe00fPnyYlJQUPDw8mDBhQpHH+Pj4MGvWLMDoXOXo6MiaNWvyb5W++eabNGnShCNHjuQX1bw5mUsacxsREYGXlxdBQUH5226//Xa6d+/Oli1b+PDDD5k0qWL/YFm8eDFTp07F0dGRLl26MG/evKuOCQwMZNSoUVa53sqVK1m50lgg/ezZswD89ttv+ef38fHhnXfeyT/+gw8+oHPnzkycOJGNGzfSokULdu7cyebNm2nWrBlvvPGGVXIJUVbvR76Pt5u33JIuK0ub0NZ+AYMxbiHuLmb/27n7XyjpPDLFZektXrz4mrdz+/btm3/8sGHD9KBBgwqdIyYmRgM6PDw8f1tISIj29PTU8fHxRV43NjZWA7pbt25X7du+fbsGtK+vr05KSrLON1pK06ZNu+bPo6jM5XW9hg0bXvWZmJgYPWrUKO3n56ednZ11QECAnjhxor5w4YLVcpWV/E5VbRfTL+q7VtylF/2+yOwoNoEy3JpWWpu7sqFSygejI1gKUEdrnXHF/h+Bfp7jZB0AACAASURBVMBQbYxZLlJoaKjO64BUksOHD9OiRYvrC13FtGrVigcffJCpU6fmb1uxYgXDhw/nwoULVKtWjYsXL1K7dm2eeeYZ3nrrLRPTioomv1MiKycLrTXOjpVuygerU0rt1VpbNBylLHNND1BKLVVK7VNKHS+wvYVS6jmllL8l59Nanwe+wRgKNbXgPqVUb6AvxrPotZZmFdcvPT2dqKioq5YCnDt3LkOGDKFaNeN50LZt23B2dmby5MlmxBRCmOD0pdOkZqbi5OAkRfg6WDLFpcJ4Rpy3cGoaUHAuvwTgXxidrt60MMdk4GbgZaVUV2AX0BCjs1Y28LjWuuxTB4kyO3DAeOT/9ddf06tXL+rUqcPMmTM5fvw4y5Ytyz9uwIABpKenmxVTCGGCV355heSMZL4dYL1lQqsiS1rET2AMVfoMY0GGdwru1FqfBbZjDGOyiNY6DqMQvws0ACYCPTGmzOyitf7W0nMK64iMjKRp06a89tprDB06lHbt2pGcnMzu3bvx8/MzO54QwiS7zuxiz997uKfpPVKEr5MlE3qMBvZhtE61Uqqoh8vHMG4lW0xrHY/RMpZ7mzYkMjKSNm3aMGTIEIYMGWJ2HCGEDdBa837k+9Rxr8N9ze4zO06lZ0mLOBjYrEvu3RUH2M6AXnHdIiMjadv2ygnUhBBV2c6zOwmPC2d0m9G4OrqaHafSs6QQZwFu1zjGH7hU9jjClmitOXDggBRiIUQhW2O3UqdaHe5tdq/ZUeyCJbemDwHdlVKqqFaxUsoN47luhLXCCXMppfIXFRBCiDzP3fQcj7Z+VFrDVmJJi3gJ0Bx4N3fpwnxKKUdgDlCPArNvCSGEsB9aa86nnQfAx902poG1B5YU4vnAeowezX8CQwGUUsuBUxjLH/6gtf7S2iGFEEKYb/tf2+m7vC8RcXLj05pKXYi11tnAXcDrgAvQDGPM8GCgGjADuL8cMgohhDCZ1poPIj+gtnttWtdubXYcu2LJM2K01lnAdKXUaxiFuDbGrFdHcgu1EEIIO/TL6V84cP4AU2+ZKrNoWZlFhThPbmetKCtnEUIIYYPyWsP1POoxqPEgs+PYnVLfms6dW3q8UsqzPAMJIYSwLVEJURyKP8SYtmOkNVwOLOms1RJ4D/hLKbVAKWXR6hK2xOwVp4SwF/K7VDU0927OyrtXMrDJQLOj2CVLCnF94FXgHMZ0lzuVUnuUUo8rpTzKJV05cHR0JDMz0+wYQtiFzMxMHB0dzY4hylF6lrGYS5BXEM4O0houD5b0mv5ba/0vrXUj4A5gJdAW+AijlfyBUiqknHJajaenp0xSIYSVJCUl4ekpT6vsldaaET+O4M1dli6oJyxh8XrEAFrrdVrrezFWSnoVOA+MBfYqpXYopUblzrRlc7y9vUlISOD8+fNkZGTIrTUhLKS1JiMjg/Pnz5OQkIC3t7fZkUQ52fTnJg7HH6a5d3Ozo9g1db2FKHed4oEYz4/9czdrIB54Q2s997ouUEqhoaF6z549pTr28uXLxMfHk5ycTHa2jLoSwlKOjo54enri7e2Nq6tMc2iPcnQOD6x6gPTsdFbevRInhzINsrFpKyNO8/a6KP66mEa9mu482zeYQe38r/3BEiil9mqtLepDVeafrFLKH3gM43mxP5AD/AAsBNpjzLQ1WylVW2v9almvUx5cXV2pW7cudevWNTuKEELYpE0xm4hKiOJft/3LbovwiysOkJZpNMZOX0zjxRUHAK67GFvKolvTynCnUuq/wB/ANMAZ+BfQSGs9SGv9g9Z6OtAU2ItRqIUQQlQinx/6nMAagdwRdIfZUcrF2+ui8otwnrTMbN5eV/FTZJT6zxyl1CsYLeAGGFNbbgU+AFbkzrhViNY6WSm1CphunahCCCEqSliPMM6knLHL1jDAXxfTLNpeniz5Cb8OJGEU3w+11odK8Zm9wOdlCSaEEKLi5fUbquVWi1putUxOU348XJ24dPmqNiT1arpXeBZLbk2PB/y11hNKWYTRWq/RWj9StmhCCCEq2rpT6xi2ehjnUs+ZHaXcLNv9J5cuZ+HooAptd3d25Nm+wRWex5JxxPO11inlGUYIIYR5snOy+SjyI1KzUvF2s89hab+eOM9L3x+gS1Mf3rq3Df413VGAf013/j24TYV31ILr6DUthBDCvqw/tZ4TiSd4u+vbODrY34xpJ85dYtySvQT5ePDesPZ4uTtzb4cGZscqvkWcu8jD3WU5qVKqjlIqTCn1fNmjCSGEqCjZOdl8uO9DmtRsQp/APmbHKRdODormfjVYOOomvNxtZ7rOkm5NJwDfK6UOK6WeV0oFlXQipZSrUuoOpdTXwElgJHDQelGFEEKUlw0xG/gj8Q/G3TgOB1WmSRdtVmZ2DlprGtb24JuxnWjgXc3YsX8ZvNsaptc0vu5fZkq+EmfWUkrdB8wEmmHMlnUW2AOcwSjUbkBtoDnGvNPOQCawCJiqtY4rx+yFWDKzlhBCiMIysjNYe3ItdzW6y64KsdaaZ5btQwNzHrgRYzJIjKK7aiJkFhiu5OwOA+ZB2wfKfD2rz6yltV4OLFdK9caYmKMHMKCIQ7OBSOA74FOttf12txNCCDvk4ujCwMb2t8zhe5uOsyLiNJN7N/unCANsfL1wEQbj/cbXr6sQl0WpOmtprTcAGwCUUsFAAEZLOA2IAw5qrWVJIyGEqGSycrIY99M4hgQP4faGt5sdx6r+t/8vZm84yj3t/JnQs0nhnYmxRX+ouO3lyOJe01rrKKDi5wATQghhdWv+WMPOMzsZGjzU7ChWFR6TwORl+7gpsBaz7m1TuDX8V0TxH/SqX/7hrmA/DwKEEEJYJCsni/n75tPcuzk9A3qaHceq0jOzae7nyfyHQ3F1KjAU6/Re+PxucK8FTles1uvsDr2mVmxQpBALIUSV9b/o/xGTHMP4G8cXbjFWYjk5Rgfkzo19+O+Tt+Lt4VL4gI0zwK0mjP0ZBv4HvHKXT/BqcN0dtcpKJvQQQogqKDMnk/n75tPCuwU9GvQwO45VZGXn8Njne7i1sQ+Pd21U9B8X938GGSnGLeiaAaYU3itJIRZCiCrISTnxfMfn8XD2sIvWsNaa6asOsiXqHHe09iu8889d8Nt7cM9845a0u20tZiGFWAghqiClFN0bdDc7htV8tv0kX+yIYWy3Rjx4U8A/O2J2whf3gocPpCcaz4FtjDwjFkKIKmbViVXMC59HZnam2VGs4qdDfzNj9SH6tfLj+b7N/9lx6jf4YjBUrwOPrAFPv+JPYiJpEQshRBWxOno1YeFhnEk5g7ODM428GnFX47vMjnXdzialc2P9mrz7YAgOeUsbnvrNaAnXqAsjV0GNeuaGLEGZCrFSqjnQAqiutV5i3UhCCCGsbXX0aqb/Op307HTA6Kz12m+voZSif6P+JqcrG601Sike6tSQITc1wMmxwE1eV0/wawP3LzKKsQ2z6Na0UipEKbUHYzGH5RhzSuft66aUSlVKFTUFphBCCBOFhYflF+E86dnphIWHmZTo+qRmZDF0wQ42RxlLGuQX4YSToDX4tYZH19p8EQYLCrFSqhmwBQgGwoAfrzhkKxAP3GetcEIIIazjbMpZi7bbspwczVNLI9n1R3z+uGEA/tgKH9wCuxYY7ytJb3BLWsTTABego9Z6MrC74E5tLOP0G3DT9YZSSj2slNK5r8eu93xCCFHV+XkU3VGpuO227M21R1h/6G9e6d+SXi1uMDZG/wxfPgA1G0KrQeYGtJAlhbgXsEJrfbiEY2KA63oirpRqAPwHuHQ95xFCCAGXsy/z3dHvmNhuIm6Ohad0dHN0Y1L7SSYlK5ulu2KYvzWahzs15JFbA42NJzbDVw+Ad5DRMat6HVMzWsqSQlwTuNayFA4YreYyUcao8s+AC8BHZT2PEEIIY43hSZsnMf236QTUCGB65+nU9aiLQlHXoy7TO0+vdB21DpxOpFszX6YNaGlMRJJyHpYOB+/GuUXY1+yIFrOk13Qc0OQax7QC/ix7HCYCPYHuuV+FEEKUQWZ2JpO3TGb76e1Mu2UabX3b0ta3baUrvFeaOag1l7Ny/umc5eED934CDW4Gj9rmhisjS1rEm4ABuesRX0UpdRPG7et1ZQmilGoBzALCtNZby3IOIYQQxtCkZ35+hp9jf+bVTq9yX7PK3Yc2PiWDhz/dyfG4ZJRSuDk7wrENELXWOKD5nZW2CINlhfjfQBawVSk1ntxnwUqpVrnvVwHJwDuWhlBKOQFLMJ4xv2Tp54UQQvzj9/O/sy12Gy92fJEHgs1f1OB6XM7KZszne9j1RzxJ6VnGxqPrYekw2PYO5OSYG9AKSn1rWmsdpZS6F/gaeC93swL25369CAzWWseUIcdUoB1wm9Y6rbQfUkqNAcYABAQEXONoIYSoGtrVaccP9/xAA88GZke5Llprnl++nz2nEnhvWDvaB9QyWsHLHoY6LWDYMnCo/DM1W/QdaK3XAkHAZGAZ8BOwAngWaKK13mRpAKVUR4xW8Gyt9W8W5vlYax2qtQ719a18D+iFEMJasnOyeeWXV1h/cj1ApS/CAPM2Hmdl5F9M6dOMu9rWg6gf4ZuH4IZWMOK/UM3b7IhWYfEUl1rrixgTelz3dCwFbkkfBV693vMJIURVlJ2TzdRfp/LDiR8I8goyO45VZGbnsO3YOQa39+fJHrn9hP/Yakxb+fD34F7T3IBWpIx5OEy6uFI1gYRSHh6mtX6quJ2hoaF6z5491gkmhBCVRI7OYfqv0/n++Pf8X8j/MfbGsWZHspr0zGyUAleywMnVmLoyMxVcPMyOViyl1F6tdagln7G4RayU8sFY8KE+4FzUMVrrz0t5usvAp8Xsa4/x3PgXIApj1i4hhBC5tNbM2DGD749/z/gbx9tFEf4zPpVZa4/wr0Ft8KrmDId+gPWvwMgfoFagTRfhsip1IVZKuQJzgEcpftIOBWigVIU4t2NWkVNYKqWmYxTixVrrT0qbUwghqpLqztV5vM3jjL9xvNlRrltiWiaPLNrNueTLXEi5jNcfq2H5o+DfAdzt43lwUSxpEb8DjAcOA98ApzGGMwkhhKhAWmvOp53Ht5ovkztMBjBmmarEMrNz+L+vwjl5PoUlo2+mUdwGWD4a6t8EDy03ljW0U5YU4gcwhirdpLXOLKc8QgghSqC15u09b7Mmeg3LBiyjTrXKNa9yUbTWTPvhINuOneet+9pyi44winCDjjD8W7suwmDZ8CUPYENFFWGt9XSttZLb0kIIYdBa8+7ed1lyaAl9A/vi624fwzbjUzLYciSO8d0b80BoA6jfETo+DsPtuyWcx5IW8UHA9ldYFkIIO6S15j8R/+Gzg5/xYPCDvNDxhUp/OzpP7equ/G9iF2qe3goZDcCtBtzxptmxKowlLeJ3gHuUUs3KK4wQQoiirTy+kgUHFnBfs/t46eaX7KII/346kTdWHyIrOwfv49/j8PX98PNbZseqcJZMcfmtUqousE0p9QEQDiQWc6ws2iCEEFbUN7AvyRnJPNTyIRxU5Z/W8WxiOqMX78ZRKSb67MHzx4kQ1AW6PW92tApn6TjiWhjPiqde4zjHssURQghR0KoTq+jRoAfVXaozotUIs+NYRcrlLEYv3s2l9Cw29IzF88cpENQVhi4Fl2pmx6twlowjfhGYBlzAGL70FzJ8SQghys2i3xcxe+9sngx5knE3jjM7jlVk52gmLY3k8JkkFg1rSb11k6BRdxj6NTi7mx3PFJa0iMcA0UAHrXWRt6SFEEJYx5JDS5i9dzb9AvvxWJsi5z2qlKLOJvPL8XNMG9CKrm0Cod4aqFGvyhZhsKyzlh/wgxRhIYQoX18d/oq3dr9F74a9+XeXf+PkYPFsxDarZb0a/NYnlpEpC425o2s3rtJFGCxrEUcD9rPchRBC2KCUzBQW/r6QXgG9eLPrm3ZThH85dp7YhFSGOPxErY1ToElvyMkCxyKXLKhSLPk3/CEwTSnlp7U+W16BhBCiKvNw9uCLO7+gtlttnB3so0gdj0tm/Jd7GeO+GdI+hKZ94cElUoRzWVKIVwHdgV+VUq8Deyl++FLM9UcTQoiq4/tj33Mk/ggvdHwBPw8/s+NYzYVLl3lk0W6GO2xgQtoCaNYPHvjcWNZQAJYV4j8wVlZSFL90IbnH2Me9FCGEqACrTqxi2q/TuKXeLWTmZOLiWNwCd5VLemY2Y5bsJS7pMg/2aQdnBsC9n0oRvoIlBfNzjCIrhBDCSlZHr+aV7a/QsW5HwnqE2U0RBtgSdY5zMUd4d9gdBLWpC3oY2MGMYNZmycxao8oxhxBCVDlrT67lpV9eosMNHfhPz//g5uRmdiSr6pf8HX3dp6G8WgF1pQgXQ24hCyGESdwc3ejoZ7SE3Z3sZwjP6v1naPPnEgJ2v4FqMRD825sdyaZJIRZCiAp2Pu08Pu4+dG/QnW71u9nFAg559pyM58C3M+jv+CW65SDUvZ9I7+hrKLYQK6UWYjwTfklr/Xfu+9LQWuvRVkknhBB25uc/f+aZn59hTvc5dK3f1a6K8KkLKXz8+WI+dvySjOC7cbn3U3CU9t61lPQTGoVRiN8E/s59XxoakEIshBBX2Ba7jae3PE2zWs0IqRNidhyrSkzN5NFFu7mgWxLX+z/U6TRMinAplfRTCsr9evqK90IIISz06+lfeWrzUzSp2YT5vedTw6WG2ZGsat93s3CK9+Oj0YOp06i22XEqlWILsdb6VEnvhRBClM6fSX8ycfNEgryC+Lj3x3i5epkdybp+fpuuJ96hRcij+DZ63Ow0lU6pF31QSk1VSnW9xjFdlFLXWqtYCCGqlPqe9Xkm9BkW9FlATTf7mrL/4NcvweaZ0HYIvve+Y3acSsmS1ZemY0xxWZKuGGsWCyFElRf+dzhR8VEopRjafCi13GqZHcmqji97mVZR7xNe6w4Y9AE4OJodqVKypBCXhhOQY+VzCiFEpRMZF8n4n8bzxs430Nr+JiX8/dQ5Eg5u5Ce33rQc97kU4etg7S5tHYDzVj6nEEJUKvvP7WfcT+PwrebLO93esashSmjNmfhEHv1iH9XdpvHNE91xc7WfaTnNUGIhVkptumLTKKVU9yIOdQQaAA2Br60TTQghKp+D5w8ybsM4vN28+bTPp9SpVsfsSNajNWx8nbS9G8jJeI4Px/fAt4b9zAhmlmu1iLsX+GcNBOa+rpQDXAC+AZ62Qi4hhKiUPjv4GTVca7Cw70Ju8LjB7DjWozX8NB22z6VByEgWtLuVYD9Ps1PZhRILsdY6/xmyUioHmK61fr3cUwkhRCX1xm1vkJCeYFdrChtFeBpsDyOz3SicB7xLOwdrdzGquor9SSqlJiqlOhbY9AiwsvwjCSFE5XIs4RhP/PQEiZcTcXV0ta8iDLDtHdgexudZvVlUcyJIEbaqkn6ac4F+Bd4vBAaVbxwhhKhcTlw8wWPrHyMqPoqky0lmxykXu6t1Y172YLY0fo5HuzQyO47dKenWdDrgWuC9HXX7E0KI6xedGM3odaNxUA582vdTGtRoYHYk69EajvyPo7W68eiqBPxrP8LyYe1xdJBSYG0ltYj/APoqpQr2NrC/wXBCCFEGp5JO8di6xwD4tO+nBHoFmhvImrSGtS/ANw+x/Mv5uLk4snDUTVR3lUUcykNJP9X5GLen/yowBm66Umr6Nc6ptdbyb0sIYdcclSM3VLuBGbfOoJGXHd2u1Rp+fA52fQy3/B8PdXiCgZezqFdThimVl5IWfZinlIoD+gP1gB5ADHCyYqIJIYTtuZB2gVputajvWZ+v+n9lX5N15OTAj8/C7k840fRRGvWeQYB0zCp31xq+tBRYCvnDlz6T4UtCiKrqr0t/8ei6R+kV0Itnb3rWvoowwN+/w97F7K73MPcf6MXHIXH0aWVnPcBtkCW3kF8DtpRTDiGEsGlnU84yet1okjKS6N+ov9lxykfdtqzvupwxay/xQGgDere0owlJbFipC7HW+rXyDCKEELbq75S/Gb1uNBcvX2RBnwW0rN3S7EjWk5MDa56BwNvY5dGDJzekcksjH2YOamN/LX4bVWwhVkoF5P7jaa11doH316S1jrnuZEIIYQOyc7J5YuMTXEi/wPze82nt09rsSNaTkwOrJkDEF2S51ebJnV408K7GRw91wMVJng1XlJJaxCcxhiu1AI4WeH8t+hrnFUKISsPRwZFnOjyDu7M7N/reaHYc68nJhh8mQOSX0O15nLq/SFjQBep5ueNVzdnsdFVKSQXzc4yimnjFeyGEsHvx6fHs/XsvvRv2prN/Z7PjWMf+ZbDxdUiMBWd3yEwlq+sL7GrwOJ2VonNjH7MTVkklDV8aVdJ7a1FK1QbuwRgm1QbwBzKAA8BnGD21c8rj2kIIUZSE9AQeW/8YscmxtK/Tntrutc2OdP32L4NVEyEzzXifmYp2cOLbaGde2rCTdU91pdkNspqSGWzhFvL9wIfAGWAzxljlG4DBwCfAHUqp+7XW0hoXQpS7xMuJjNkwhpikGN7r9Z59FGEwWsJ5RTiXysmiS8yHTOixRoqwiWyhEB8FBgKrC7Z8lVIvAbuAezGK8nfmxBNCVBWJlxN5fP3jnLh4gv/0/A+d6nYyO5J1XL4EiX8Wucvf4QJP925WwYFEQSX1ml5YxnNqrfVoCw7eVMz2s0qpj4A3gO5IIRZClLONMRs5dvEYYT3CuNX/VrPjXL+MFNj9CWwPK/YQXcMfBxmmZKqSWsSjynhODZS6EF9DZu7XLCudTwghijW46WBCbwgloEapR2vatlO/wYap0LgXh5yaE3RkAe4qI393mnbh9yYTucnEiKLkQhxUYSmKoJRyAkbkvl1rZhYhhP1KzUzl+W3PM/7G8bSs3bJyF+GsyxD+OWSmwq2ToEkvGPMz1Avh8Vmb6JCpec5pGfXUBf7StXkr6wH2HmrK9oFmB6/aSuo1faoigxRhFtAaWKO1XlfUAUqpMcAYgICASvzLI4QwRWpmKk9sfILIuEgGNh5YeWfMysqAyC9g6zuQdBqa9IbOE0EpqBdCbEIqpy+mcZrb+CHjtkIfVRfTijmpqCi20FnrKkqpicAzwBHg4eKO01p/DHwMEBoaKr2qhRCllpaVxoRNE4iIi2BWl1n0btjb7EhlE70F/jsBEmOgfkcY9AEEdQOluHDpMu9vPsEXO4pvV8nyhuazuUKslHoSCAMOAb201vEmRxJC2Jn0rHQmbprI7rO7eeO2N7gj6A6zI1kmOwsyksG9FlTzgep1YMC70LiX0QrOdeB0Iot+/YP7OzSgRT1P3vwxirTM7Pz97s6OPNs32IzvQBRgU4VYKfUU8C7wO0YRjjM5khDCDjkoB9wc3Zhx6wwGNB5gdpzSy8mG37+DLbPAvz3c+wn4tYbHNwJwOSubr3bGkJaZzRPdm9CtmS8/P9uDBt7VAKjp7sLb66L462Ia9Wq682zfYAa18zfzOxLYUCFWSj2P8Vw4EuittT5vciQhhJ3JyM4gLSsNL1cv5vWcV3lWF8rJgUPfGwX4/FG4oTW0Gpy/OztH89/I08zZcJTYhDR6Nq+D1hqlVH4RBhjUzl8Krw2yiUKslHoVeB3YC/SR29FCCGvLzM5kys9TOJtyli/v/BJnx0q0sMEvc2DTDPBtDvcvhhYDwcFYHSkiJoEXVxzgyNlkWvvX4N+D23BbE5/K80eGML8QK6VGYhThbGAbMLGI/4BOaq0XVXA0IUQltzp6NWHhYZxNOYuroyvp2em8dPNLtl+EtYaoH41nv/VDod3DUCsQWt0DDo4AZGTl4OLkgIerExnZObw3rB13tq6Lg4MU4MrG9ELMP+OVHYGnijnmZ2BRhaQRQtiF1dGrmf7rdNKz0wFIz07HycGJGi41TE5WAq3h2AbY/AaciYQ290P9T8DzBmhzHwBHzibx9too3FwceX9Ye5rd4MlPT3eTAlyJmb7ys9Z6utZaXePV3eycQojKJSw8LL8I58nKySIsvPjpHk31xzb4tDd8dT+kxcPd78Ogj/J3/xmfyuRvIrkjbBu7T8bTxt+LvLVwpAhXbrbQIhZCCKs7m3LWou2m0doYcnRmHySdgbvmQshwcHLJP2Tt72eY8HUEDkoxtmtjxndrjFc1G7+9LkpNCrEQwq5k5WTxyYFPcHF04XL25av2+3n4mZCqCKd+hc3/MopuyFDo+LjxcnIFIDk9kwuXMgj08SA00JuhHQN4onsT/LzcTA4urE0KsRDCbpxMPMnLv7zM/vP7udHnRqISogrdnnZzdGNS+0kmJgT+3GUU4OjNUP2G/M5XeQU4PTObL3ac4v3Nxwn08WDF+M74VHfl9btbmxhalCcpxEKISk9rzTdR3zB7z2xcHF14u9vb9AvsV6jXtJ+HH5PaT6J/o/7mBV3zHOyab8yG1ecNCH0UXIxxvtk5mhXhscz96RinL6bRpakPz/YNlmFIVYAUYiFEpXcp8xIL9i+gww0deP3W16lTrQ4A/Rv1N7fwgvHst1YQuNWARt2hRl3oOAZcPAodtnzvnzz/3QHa1vfirfvacmsTH1Piioqn8nrdVXahoaF6z549ZscQQlSgn//8mc7+nXF2cObMpTP4efjZTgvy74PGLegj/4NeU6HLM1cdsjP6AqmZ2fQIrkN6ZjY/Hz1Hn5Y32M73ICymlNqrtQ615DPSIhZCVDqJlxOZuWMma0+u5dVOr/JA8APUrV7X7FiGuCPw8yw4+D241oBuL0Do6EKHHPwrkbfXRbEl6hwdGtaiR3Ad3Jwd6dvKRjqSiQolhVgIYfMKPuut5VaLzOxMYxnDdhMY3HTwtU9Qkda9aHTI6vosdHoCqnnn74q5kMrsDVH8N/IvvNydefGO5ozsHGheVmETpBALIWzalTNkxafHo1A8GfIkY9qOMTkdEB8NW2dD9+ehZgD0n2O0hD1qX3XooTOJrDt4lie6N2Zst8Z4uctYOKlmbgAAHOJJREFUYCGFWAhh44qaIUuj+e7Yd4y9caxJqYCEU7D1bYj8ChydoWlvoxB7B+UfkpSeyYKt0Xi4OjGuW2P6tvJj23M98fV0NS+3sDlSiIUQNuvMpTOcSTlT5D7TZsjSGn58DvZ8BsrBmITjtqfB85/nu+mZ2Sz57RTvbznOxdRMhnZsAIBSSoqwuIoUYiGEzcnMyeSLQ1/w4b4PUSg0V4/uqPAZslLjjee9SkFOFnQYCbdNBq/C6/v+fPQcL3y3nzOJ6XRt5stzfYNp7e9VsVlFpSKFWAhhU6Lio3hh2wscv3icHg160NGv41W3pyt0hqzkv2H7XNizEB5ZA/4djOfABYYYaa35//buO76KKm3g+O/JTQ8plCQQqgQJFqoNsACC4IorNkBXxV11dde1vPjKvqjgIitrW1asa3v3ZdeyiEjRVQwWEBE7KqASpUSFBBLS20275/1j5oaQOzcNknsJz5fP/Qw5c87MuXdO7pOZOeeMu9pDVLiLztFhJMdFsnDaUEan6lhg1TQNxEqpoBIfYZ09Pnb2Y4ztPRaAzpGd23+GrLL9VgD+9DmorbLmg46xJgqpH4Q/2pHHA29t49ikTjw0dShDeiWw4sbROhZYNZtO6KGUCiiP8bDs+2V8nP0xC8csREQwxgQ2kNVWw6IhULoXhky3hiJ1TT0oy9Y9RTyYnsH673PpER/JzHMGMu3k3gGqsAoWOqGHUiro1R8T3DWqK5GuSHaX7uaU7qdQWl1KbHhsYIJwRQF8/bI1/aQrDH7xACQOgsSBPlmXfPoTs5dvISE6jLvOO46rRvUlMszV/nVWHYIGYqVUu2k4Jnh/xX4Apg6cytyRcwMTgN1F8PFT8NETUFkEPUdA71Ph+AsOypZT7KasqpZjusUwJi2Rm8YN4Pox/YmL1LHA6tBoIFZKtRunMcEAG/ZsaP8gXO2Gjx6HjY+BuxAGnQ9jZ0P3wQdlK6qo5un3d/CPD3dxSr8uPH/tafSIj+L2SWntW1/VYWkgVkq1OWMMb//4dnCMCTbG6mwV4oJN/4Lep8G4OyBl+EHZ3NW1/HNjJk+u20FRRTUXDE3hvyf6XqZW6lBpIFZKtalPsj9h0ReL2Jq3lVAJpcbU+ORpszHBm5fCu/OhaDfE9YS+oyA3A65dA2FRcMN6iEpwLPr8Rz9y3+ptjE1LZNakNE5I0bHAqm1oIFZKtYmiyiL+uP6PbMzaSPeY7swfPZ9QCWX+x/PbZ0zw5qXw+i1QXWH9XLwbtrwC3QZZk3PE9zwoCBtjeGvrXqIjQhkzMJErRvZhcK94Rvb3nTNaqcNJA7FS6rAqry4nOiya2PBYPMbD7SffzmWDLiPCZU3tGBIS0j5jgt+dfyAI11dd5jMb1ofb9/PgW9v4encRE45LZszARKLDQzUIq3ahgVgp1Sr1hyF1j+nONSdew/bC7azJXMPrF71OfEQ8z5zzjE8nrMn9J7fdZBz7voWtr1r/L9rtnKde+jdZRdy/ehsf/LCflPhIHrp0CBeP6NU2dVPKDw3ESqkWazgMKbssmwWfLEAQpg6cisd4ANqnJ3T+Lti6DLYuh5xvrQcxDDof4ntB0c+++eMPBNrv95WwdU8RcyYfx5UjdSywCgwNxEqpFvM3DKlbVDfmjpp7+HdYv9NVfC8YfTOccp3V8/nTZ+HjJ6DPKDjvr3D8FOiUBJuXUrPqZkLr1bMmJJJX435DyQc7ue7M/kwZ2pMJxyUTq2OBVQBpIFZKtcie0j1+hyF5J+g4rBp2uir62XoMYek+GH83jPoDjPw9JBw8veTK2tPZUH0d/8USUiSPLNOVB6um8Z/tg7ixVzUAISGiQVgFnAZipVSzbM7dzD+/+Sfv/PSO3zxtMgzp7budO1199ZIViBt0vPJ6KD2DPVWjWcbog9KT4yJ0Mg4VVDQQK6UalZGfwYJPFvBlzpfEhsVy9QlXkxydzKIvFh3+YUi1NbD7U8h4EzolW5egS/xM9uEnPb+sitVbs8kqdAjeQE5x5aHVUanDTAOxUke5hr2fbx1xK+N6jyPfnU+v2F7ER8RT4C5g9qmzuXDAhcSExQCQEJHQumFIDe/3jr8bIhPgm+XwfTpU5ENIGAz7lZW/GZ2uit3VrPlmH69/ncWG7fup9RjiIkMpdvtOHpKSENWqz0mptqKPQVTqKNaw9zOAS1yEhoQyPGk4z058FuDwPZaw4f1esGa4Sh4MeT/AsZMg7ReQejZExtWV8el05YrENeUxZMg0qms9nLrgHQrKq+nVOYpfDk3hl0NSyNhbzJ0rtlJRXVtXLirMxX0XD+bC4c6Xs5U6VPoYRKVUizj1fq41tYQRxk3Db6pLO+QgXFkKP34I/5npe7+3ugKK98Dt28Hl+5Xk1Onqr9XT+So9kXVDIMwVwpzJx3NMYgzDeyfU1fX4lDhEhIfSM8gqrCAlIYpZk9I0CKugo4FYqaOQMYYt+7f4fdhCZW0lQxOHNr4Rp0vMQ6ZZ62prIGuTlR6XAj+kw7Jr/G+rOMsxCIP/TlfRJZVU1tQSEerikpOcJ+G4cHhPDbwq6GkgVuooYYwhoyCD1btWk56Zzp7SPSREJFBYWeiTt8nez05Dil67Gba/Zz3TN3MDVBbDxHutDlep42HGKlj5B2vO54bq3e/1eAxb9hSxLiOXS07q6bfTVUWVFYSVOtJpIFaqA3HqeDW5/2SKKou48s0rySzOJFRCGZkykhuH3UhVbRUPfPpAy3s/O83jXOOGzf+Gzv3gxIuh/zg45ixrXVQC9B8LE/7keL+3dswc3vpqD+sycln/fS55ZVWIQP/EGFISotjjEIy105XqKDQQK9VBOE07OWfDHMCa33lE8ghmnDCDCX0m0Dmyc125qKyveWTnCvaGQHcP3Nr3XOfezz9/ap3p7v7MuRczAAK3fu23jk73exd5LmNY1SjmrvqKztHhnHVsN8amJXHmsd3o2imCWo/hjuVbfDpdzdKxwKqD0F7TSnUQY18eS547zye9R0wP1ly6xrmQv17MZ8+FmCQoyIQxs6z05y+GHe9C1wHWPd3qct/txfeGmVsdd/VZZj7XLv7McUhRz4Qo/vHrUxiQ1AlXiG/HsJVf7tFOV+qIoL2mleog/F1i9tpXto9NOZv4Yt8XzDplFhGuCMcgDPjtkAU4X2KuroD0O63/h8fC6Jus4HzeQ9Z435iufocUhY6/G4/HsHN/KZt+KiSn2M1NZx8LWJ2unIIwQFZhBWndY/1WUztdqY4saAKxiPQC5gPnAl2BbGAlcI8xpiCQdVOqPTldYp63cR5ZpVlkFmeyad8mdpdaHZ5iwmKYljaNgZ0HkhydzL7yfT7bq+t4VVkCe7dYTyjK+c56+b3EDNywHpJOONCbuWtq3SqnS8wLa6bz7XvdyV6+pi7gdo0J54YxqYS5Qrj/4sFc8dwnZBf5PixC7/eqo1lQBGIRSQU2AknAKmAbcCpwK3CuiJxujHH+c1+pINfU2W197ho3Cz9f6DO2113r5sXvXsRjPIxIHsGvjvsVI5JHkNY5jdAQ69d4ZuIo5u1agbvepd1Ij4dbowdYP+xYC0uvsv4fHgtJx0FYDFSX+VYkvjf0OHj4UrG7mm+zitmZW8aCN7+lzGFIUWhOKVNP7s3wPgmM6JNA/26dCLHr0z+xE/9z7iC936tUA0ERiIEnsYLwLcaYx7yJIvI3YCawAPhdgOqmFNCygFq/zLwNc3Eb62k/2WXZzNswh7yKPFITUsksziQ1IZWRPUayr2wf5yw7B4Nzv418dx5fX/U1UpEPhT/Bvu2QsRZ6nwq9T2Xy50vBFPNI5wT2hrroXlPLrQWFTM55C84D+p4OVyyzAnBcTxDhs9ee5sQv5hAlVXX7qTDhvJlwDdnv/cDO3DKuH9OfQd3jWLsth1uXfNXo+631GO67eLDf9d7Ly3q/V6kDAt5ZS0T6AzuATCDVGPuJ4ta6WKxL1AIkGWMc/nS3HI7OWo++churitLJDRUSawxT4idxy9S/tUm5jrqvI6GOrSnzxs43uPuDOVRx4B5nOKHMP/Peg4Kxu8ZNgbuA/Mp8CtwF3LH2NgprncfBek1Pm86ckXPwGA/PfPQXXsxYQqFDh6UetYY12Xm+naTGzIZxd2DmJSAOQdwgyLxCKmtqydhbQm5JJftLK8ktqeTv63YwvuZ9/hi69MCjAmum8ZrnDGuf8ZHcd/FgxqYlkVtSyXfZxaQmdWLqUxvJKvS9xNwzIYoPZ5/d6PtVqiNrTWetYAjE1wHPAs8YY25wWJ8OTAQmGGPe9bedQw3Ej75yG8+XpuMOCalLi/R4uKpT41/SrSnXUfd1JNSxqTLuGjdl1WVU1FTUvVzi4g+rf0uep8Rne2G4mDbwUmaPsoYJnfvKBPY43Kf1YQyLq+PpW15M12PORi54xEpf0IM3wmFety4+dfzT/gLOSptBRWQyJZE9KArvTkFoEsnJKQzunUD2vAH0INdnV1l0I2XeDrIKKxh9/3tN1w3rL9+t90wiJsL5otnKL/c4XmLWeZzV0e5I7TXtvTn0vZ/1P2AF4oGA30B8qFYVpeMOCzkozR0Swisl6XR6/R5Sug7g3NFXALD4jT/j8Vgn7q+UpON2+S/nlZoyjDEnTcFdWe63zKqidK4uzefVtY/R0In9Tm+0jmds+4ARg87k5327WPPJ4rr1/va1tCSd6Nfurks7c+jFpPUdRkbml6zfvJylfsqtKE4nYtVdVkK9P+LOG3UdvZP7s7w4HXeob7mVxencAmz46nU27VoHxuD990rJ2477WlKSTvkrv8VjDDMveISoyBiWr3+SzdkfY4whvXyTY7kXStPZ9sJkaqmlxngIExdPXbkagJca+exvAa5fMpkva3MOWt/bRJCPGxzmW642NZRuXgF2IJ6xN4+Iyjy6eDx0rq3ltqRu5Ib6/pr1qKklPD+KDEmivLA759jpD3aaxTX5DzOPfJ9LzCeVxXDyJ2OoqvVeNCoDdnHlyFoG907g/qqp3Bf2HNH1LjOXm3Dur57Go0BibATPzTiZbrERJMZG0DUmnPEL3/c7WYa/IAx6iVmpwykYAnG8vSzys96bntBwhYhcD1wP0KdPn0OqRG6o86T2ha4QHs5fxmlZneoC8ZM5S6jwnq00+FJvWM5r4r6NjDlpCkWlBRT6KZMbKmTl7DqonNfUgoxG6/jOl88zYtCZbN2+gUUFyw+s9LOvIlcIjxSsOJBtSxhpfYfx4ZbXeLRgpd9y+13C44Wv+aT33XEivZP7k+dyruN+O/3dzS+xzDQYZ+pnXyWuEF4s/xiA31YUERUZwwffr+SdsGwrg8PlW4AKEXa6dxJmDKEYojwH1pX5KeP9bIcXRDO0UggzIbg8IYR6XBhXN5Z3ySQ7zDmgdis+v+7nb6uvZV9hIW6JpFIiuS7vYR5ODPE5u70yv5Y7ov9EuCuEUxO71gXi/T3HM39fLveXPsfksqy6MuUmnDuqp/Lg1CGEh4YQExFKTLiLmIhQEmMjAPg87hxmF+NzmfmLOGvrYa4QJhyffFD9Z01Ka3XnKR1SpNThEQyBuCneb06fa+jGmGeAZ8C6NH0oO0msMeSE+X5Jd6nxcOfxd5EYf+ALZ8Fxd+Mx1hfXX75dQH6obyDxlvPqmWgN/egc240uNR7HMok1hj7dB/LXgXf4rDsm5UTef/tyv3WcOuE2AE47YSILQw4c1gXf3uu3fnNOOHBGfGLqKAAmjfw1fXYex5+33uNcx1rDXSceONMX+/AMP+6suvVOfzAk1liH59cT/syEvdutshKCAHdumsl+x8/Dw6JRixFx0Tk2CYDZF/yLG4uLEQnhhrVTHD+PpBrDE+etx4gLxIXL5TponVMZb/1m/GYZVbUeBEHEanyhrhBinxzMom4en4A6o8DD9Nv+VJd27203HbTdeffmckfuP3iqS6e6s9vf5Zeyueoa/vPHM33q8eClQzl9e57fgPpII4HPCqpVvFZ1Rl1aVJiL+xoJqnpmq1TgBUMg9p7xxvtZH9cgX5uYEj/J8d7hJXGTmDTqVwflPWfk9Lr/Z+z+vNnlAMLDI7gkznlfU+InERMd61iuqToe03MQAF0Skpk46vK69dt2f+a3TP334dUzsS89E/vy7U8fO5a7MG4S40+7xLF+ABc28t4A+qYMoG/KgIPKXLTducyFcZMYcvwpB+VNTupOcpI1LnbKJv/7Sk1JdKyfv8/QW7+unSIcyw1Mu407t93L3+sF1N/nl9Jr0BzC/JzRAwybfD0bVtTwj9J6UzoygzMuut5vmdYEVGh9UNUzW6UCSztr1RPsPXePhH0dCXVs7b4+e+1pem96iCSznxzpxs8jZnHKBT5N1kdrpmfUKR2VOjIdqb2mU4HtND58KQRIbOvhS0oppdShaE0g9n9NrZ0YY3YAa4B+wB8arL4HiAH+1VgQVkoppY5UwXCPGOBGrCkuHxWR8cB3wGnAOKxhTXc1UlYppZQ6YgX8jBjqzopPBhZjBeD/BlKBR4FROs+0UkqpjipYzogxxvwM/CbQ9VBKKaXaU1CcESullFJHKw3ESimlVABpIFZKKaUCSAOxUkopFUABn9DjcBGRXODHQNdDBUw3YH+gK6GClrYP1ZjD2T76GmOc59j1o8MEYnV0E5HPWzqbjTp6aPtQjQl0+9BL00oppVQAaSBWSimlAkgDseoongl0BVRQ0/ahGhPQ9qH3iJVSSqkA0jNipZRSKoA0ECullFIBpIFYBZyIZIqI8fPa66fMaBF5U0TyRaRcRDaLyH+JiKuR/ZwvIutEpEhESkXkExG5uu3emWouEblURB4TkQ9EpNg+9i80UaZd2oCIXC0in9r5i+zy57f2vaqWa0n7EJF+jXyfGBFZ0sh+WnSsRcRlt7nNIlJht8U3RWR0S95f0Dx9SR31ioBFDumlDRNEZArwKuAGXgbygV8CDwOnA1MdytwEPAbkAS8AVcClwGIRGWyMuf3wvA3VSnOAoVjHezcwqLHM7dUGROSvWI9l3Q08C4QDlwGvi8jNxpjHW/NmVYu1qH3YvgZWOqRvdcrc0mMtIgIswWpDGcDjQBdgOrBeRC4xxqxqRj3BGKMvfQX0BWQCmc3MGwfkAJXAyfXSI4GNgAEua1CmH9YXdh7Qr156Z2C7XWZUoD+Ho/kFjAOOBQQYax+TFwLZBoDRdvp2oHODbeXZ2+t3KO9bX23SPvrZ6xe3YPstPtbA5XaZD4HIeumn2G0zB4htzv710rQ60lwKJAJLjDGfexONMW6sv5oBft+gzDVABPC4MSazXpkC4C/2j79rqwqrphlj1hpjfjD2N1kT2qsNeH9eYOfzlskEnrC3p89QbwctbB+t0Zpj7W1jc+y25y3zGdZVmkSsttokDcQqWESIyJUicqeI3Coi4/zc6zvbXr7lsG49UA6MFpGIZpZZ3SCPCn7t1Qa03RzZUkTkBvs75QYRGdJI3hYda7ttjcZqax80p0xj9B6xChbdgecbpO0Skd8YY96vl5ZmL79vuAFjTI2I7AJOAPoD3zWjTLaIlAG9RCTaGFN+KG9CtYs2bwMiEgP0BEqNMdkOdfjBXg48hPeh2tY59quOiKwDrjbG/FQvrTXHegDgAnYaY2qaWcYvPSNWweD/gPFYwTgGGAw8jXV/ZrWIDK2XN95eFvnZljc9oRVl4v2sV8GlPdpAa/ahgkM58GfgJKw+AJ2BMcBarPvL79rB16st21Oz2ocGYhVwxph7jDHvGWP2GWPKjTFbjTG/A/4GRAHzWrA58W62jcuo4NWebUDbTJAxxuQYY+42xmwyxhTar/XAROATrLPZ61qz6RbkbVF70kCsgtlT9vKsemlNnb3GNcjXkjLFLaqdCpT2aANN5W/qjEgFGfsS8nP2jy35TnE61q1pg35pIFbBLMde1r+MlGEvfe69iEgocAxQA+xsZpke9vZ36/3hI0abtwFjTBmwB+hkr2/oWHvpc89ZBbVce1n3ndLKY70dqAX6222uOWX80kCsgtkoe1n/C/U9e3muQ/6zgGhgozGmspllftEgjwp+7dUGtN10PCPt5c4G6S061nbb2ojV1s5sTplGBXqgtr6O7hdW79YuDul9sXoeGuDOeulxWH/VtmQyh2PQCT2OmBfNm9CjzdsAOqFHUL6a0T5OA8Id0s+2j5kBRh/qsaZ5E3rENec96WMQVUCJyDxgNlaPxl1ACZAKTMb6Yn0TuMgYU1WvzIXAMqxfjiVY0xtegDVEZRkwzTRo2CJyM/Ao1i/VyxyY3rAXsNDoFJcBZR/TC+0fuwOTsM5avGM099c/Ru3VBkRkIXAb1rSHy7CmPZwOdAV0ist20pL2YQ9ROgFYh3XcAIZwYEzvXGPMvQ77aNGxtqe4XIrVhrYBr9t5p2N9d+kUl/o6Ml5Ywwr+bTfkQqAa62znbWAG9jOzHcqdjhWkC4AKYAswE3A1sq9fAu9jBfsy4DOsMYUB/xyO9hdWz3jTyCszUG0AuNrOV2aXex84P9Cf2dH0akn7AK4F/oM1dW4p1tnpT1h/fJ15OI811lwcM+22V2G3xTdpcMbd1EvPiJVSSqkA0s5aSimlVABpIFZKKaUCSAOxUkopFUAaiJVSSqkA0kCslFJKBZAGYqWUUiqANBArpZRSAaSBWCmllAogDcRKBQkR6SciRkQWH+btzhcRt4j0but9NbM+0SKyV0Seb+99KxWMNBAr1YHZwfd24BljzM+Brg+AsR41eB9whYicGuj6KBVoGoiV6tjmAhHAQ4GuSANPY80t7jP5vlJHGw3ESnVQIhIPXAG8Gyxnw17GGDfWJPwTROTYpvIr1ZFpIFYqCNn3cJeIyH77/u7nInJ+CzdzOdaDy19uwX5DRORR+/7xchGJtNPn2WljReRyEflCRMpFJEtE/iYiEXa+s0VknYgUi0iBiDwvIl397G4JIMA1LXxfSnUoGoiVCj59gU+xHkr+PFYgPRFYJSLjWrCdCfZyQ3My20F3KXAz8ARwqX3mWt/NwP8CGcDfsZ7tOxN4WkQuAlZjPRv4GeA74ErgBT+7/BTrsZfnNPP9KNUhhQa6AkopH2OBecaYe7wJIvIS8BYwC1jbzO2cgfVM1e+byigiXYBVWM/4nW2MecBP1gnAScaY7+xyEcAm4CqsZ/1ONMa8b68LAdKBc0VkmDHmq/obMsZUiMg3wHARiTXGlDTzfSnVoegZsVLB50cadGIyxqRjPdy8Wb2MRSQcSAb2miYeOi4ifYEPgdOAqxoJwgCPeoOwXa9KrDP2EOANbxC213k4cDY81M/29tplezb+jpTquDQQKxV8vjLG1Dqk/wx0buY2vPdlC5rIlwZ8BKQAvzDGvNhE/s8d0rLs5RcO6/bYy15+tpdvL7s1sV+lOiwNxEoFn0I/6TU0/3e2wl5GNpFvINAD2Il1ibkpRX7q1dS6MD/bi7KXFX7WK9XhaSBWqgMyxhQCVRw4M/bndeBOYBjwroi095mpt3457bxfpYKGBmKlOq4tQA8RiWsskzHmPqyez8OBtSKS3B6Vs6Vh9bze3Y77VCqoaCBWquNah/U73mQHL2PMIuD3wAnA+yKS0rZVAxE5BqtD2bqmOpQp1ZFpIFaq43rVXk5qTmZjzFNYk2scC6wXkT5tVTHbRHv5aqO5lOrgNBAr1UEZYz4CvsR6uIKrmWUWY03C0RcrGPdvuxpyNZCLBmJ1lBO9IqRUxyUilwMvARcbY1YEuj5eIjIE+BqYa4zRBz+oo5oGYqU6MBERrHHCUcCwYLkXKyIrgZOAgcYYHbqkjmp6aVqpDswOvNcDK7Am7Qg4EYnGumR+lQZhpfSMWCmllAooPSNWSimlAkgDsVJKKRVAGoiVUkqpANJArJRSSgWQBmKllFIqgDQQK6WUUgH0/6POq+mwFa4GAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 540x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = 7.240\n",
    "CD = 2.2\n",
    "A = 0.05261451671609104\n",
    "R = 6378e3\n",
    "bc1 = m/(CD*A)\n",
    "ai = R+h*1000\n",
    "year = 3600*24*365\n",
    "mu = 3.98e14\n",
    "Isp = 536\n",
    "g = 9.81\n",
    "ar = 2*np.pi*rho*ai**2/bc1\n",
    "Ln = H*1e3*m/(2*np.pi*rho*ai**2*CD*A)/(3600*24*365) \n",
    "Lp = lambda mp: mp*Isp*g*ai/(0.5*CD*A*rho*mu)/year\n",
    "ub = 800\n",
    "#plt.rc('text', usetex=True)\n",
    "plt.figure(figsize=(7.5,5))\n",
    "plt.plot(h, Ln, '--o',\n",
    "         h, 1e-5*Lp(4.3), '--o', \n",
    "         h, 2e-5*Lp(4.3), '--o')\n",
    "plt.legend([r\"$L_n$\", \"$L_p, \\kappa=1\\cdot10^{-5}$\", \"$L_p, \\kappa=2\\cdot10^{-5}$\"])\n",
    "plt.xlabel(\"h (km)\")\n",
    "plt.ylabel(\"liftime (years)\")\n",
    "plt.savefig('lifetime_P.eps')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The constraint is to satisfy a minimum lifetime $L_{min}$:\n",
    "\n",
    "$L_n + L_p \\leq L_{min}$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Momentum budget and disjoin constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the constraints spacecraft have to manage is accumulated torque due to disturbance forces acting on the spacecraft. Although there are several (solar radiation pressue, aerodynamic...) here we focus on one specific type: gravity gradient. The expression for gravity gradient is given in table 11-9A in SMAD:\n",
    "\n",
    "$ T_g = \\frac{3\\mu}{a^3}c$\n",
    "\n",
    "Where $c = |I_z-I_y|sin(2\\theta)$. We choose to just set this as a constant instead of going deeper into the modelling of where $c$ comes from.\n",
    "\n",
    "The gravity gradient can build up over time, and we investigate two ways to manage it:\n",
    "\n",
    "1. Magnetorquers: they create a dipole moment in the spacrcraft that generates a torque on the satellite in the presence of Earth's magnetic field\n",
    "2. Propulsions: small burns with attitude reaction thrusters; this requires propellant mass\n",
    "\n",
    "Both can be tied town to the mass of their system, and the optimizer should choose the lightest one. Below are the models used:\n",
    "\n",
    "For the thrusters, the amount of momentum that the system will have to dissipate is given by $h_p$(the letter h is also often associated with momentum, not to be confused with satellite altitude). This is modelled through the equations given in Table 11-13:\n",
    "\n",
    "$h_p = T \\cdot T_g \\frac{1}/{4}\\frac{sqrt(2)}{2}\\cdot 3 \\cdot 365\\cdot L_{min} $\n",
    "\n",
    "For the magnetorquers, we have applied a model that is not correct from the literature, but that makes the choice between thrusters and magnetorquer interesting. The assumption is that the magnetorquer can counter the torque(through the dipole moment it creates); so:\n",
    "\n",
    "$T_m = T_g$ \n",
    "\n",
    "Both key metrics from the systems: $h_p$ and $T_g$ will then be used in the mass models.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. Mass models\n",
    "\n",
    "Mass models for any system that's not in the component catalog uses a simple linear relationship that was tuned to generate meaningful masses, i.e.:\n",
    "\n",
    "$ m_i = \\rho_i \\cdot f_i $ \n",
    "\n",
    "Where $m_i$ is the mass of certain component not in the catalog(except for the solar panel), $\\rho_i$ is the density, and $f_i$ is the performance. $f_i$ for momentum management is both $h_p$ and $T_g$. $f_i$ for the surface is the area, and this time the constant $\\rho_i$ is given by the component selection in the catalog. \n",
    "\n",
    "Finally, we have included the structural mass as a percentage of the total system(20% in this case)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some notes on reformulation\n",
    "Most constraints in the modelling are equality constraints; they can however be turned into inequality constraints if one know on which side of the equality the pressure lies. This is what effectively allows sums of powers to be put in logspace through inequalities. Sometimes, however there is no such pressure, and it was realized later that the expression\n",
    "\n",
    "$a = h + R$, which was reformulated as $h + R \\leq a$:\n",
    "\n",
    "1. $e^{R - a} \\leq exp_{Ra}$\n",
    "1. $e^{h - a} \\leq exp_{ha}$\n",
    "1. $exp_{Ra} + exp_{ha} \\leq 1$\n",
    "\n",
    "Where $R$, $h$, and $a$ are in logspace. This assumes a certain pressure on $a$ which turned out not to be true, due to it's relationship with the orbit period $T$ and the comms link budget, where $T$ ideally wants to be larger to satisfy the link budget constraint. Therefore, this constraint had to be encoded through the Piecewiselinear formulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Catalog components\n",
    "A mini-catalog was generated with two options for each component except the battery which has 5. This leads to 2x2x2x5 = 30 possible combinations; which is manageable; but not scalable when including perhaps 100s of components. The components where found from an online database of spacecraft components [satsearch.co](https://satsearch.co), and some of the values were tweaked to play around with the problem. \n",
    "\n",
    "For simplicity of the writeup consult data.jl to see the catlog components."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this is a python notebook (the project was originally started in python, when using GpKit, so some legacy code like extracting the excel tables, was kept in python), we will just display some of the results from the code, and show that the results are non trivial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we change the resultion requirement and lifetime we get very different designs. Here we show the results from two runs:\n",
    "\n",
    "1. 15m resolution 5 years\n",
    "2. 5m resolution 3 years"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. 15m resolution, 5 years\n",
    "![results1.PNG](https://www.dropbox.com/s/1260hniqdfjy3bh/results1.PNG?raw=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Payload 1 is chosen(lower aperture than payload 2, but much lighter), and the orbit is brough to 1200km. This is the altitude that satisifies the resolution constraint given the camera aperture. Choosing this altitude comes at some significant propellant mass penalty, but seems to outweight the choice of having a heavier payload with better performance. \n",
    "\n",
    "The optimal mass comes to ~5.24kg."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. 10m resolution, 3 years\n",
    "![results2.PNG](https://www.dropbox.com/s/umcno6gn2a0tbkq/results2.PNG?raw=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time Payload 2 is chosen(larger aperture than payload 1, but much heavier), to meet the higher perfomance metrics required by the resolution. This also allows the orbit to be higher (1500km), giving us a longer lifetime and at less propellant cost.\n",
    "\n",
    "The optimal mass comes to ~4.40kg, which is lighter than the first option."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "It took less than a minute (53s from executing the command to results), of which 21s where spent by the backend solver(there is a significant overhead to compile the problem, unfortunately). \n",
    "\n",
    "This gives us the opportunity to quickly compare two different options in a couple of minutes, versus the several weeks it might take to do it manually. \n",
    "\n",
    "### Limitations and future considerations\n",
    "\n",
    "* The solver seems to sometime run into errors for no known reasons. \n",
    "* Further investigation should be given to how the solving time increases with the amount of catalog components\n",
    "* Make the magnetorquer and propulsion system integer variables(choices) from a catalog."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mdo",
   "language": "python",
   "name": "mdo"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}