{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **SpaceX Falcon9 Launch Trajectory (Booster Recovery)**\n",
    "(c) 2023 Devakumar Thammisetty\n",
    "\n",
    "MPOPT is an open-source Multi-phase Optimal Control Problem (OCP) solver based on pseudo-spectral collocation with customized adaptive grid refinement techniques.\n",
    "\n",
    "https://mpopt.readthedocs.io/\n",
    "\n",
    "Download this notebook: [falcon9_to_orbit.ipynb](https://github.com/mpopt/mpopt/blob/docs/docs/source/notebooks/falcon9_to_orbit.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Install mpopt from pypi using the following. Disable after first usage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import mpopt (Contains main solver modules)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip install mpopt\n",
    "from mpopt import mp\n",
    "import numpy as np\n",
    "import casadi as ca"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining OCP\n",
    "\n",
    "OCP: (Multi stage launch vehicle OCP : Delta-III rocket)\n",
    "http://dx.doi.org/10.13140/RG.2.2.19519.79528\n",
    "\n",
    "   \\begin{aligned}\n",
    "&\\!\\min_{x, u}        & \\qquad & J = -x_6(t_f) + \\int_{t_0}^{t_f}0 dt\\\\\n",
    "&\\text{subject to} &      & \\\\\n",
    "&\\qquad \\text{Define}&      & \\mathbf{r} = [x_0, x_1, x_2]; \\ \\mathbf{v} = [x_3, x_4, x_5]; \\ m = x_6; \\ \\mathbf{u} = [u_0, u_1, u_2];\\\\\n",
    "&&      & \\mu=3.986012e14;\\ \\omega_e = [0, 0, 7.29211585e-5]^T;\\ R_e = 6378145;\\\\\n",
    "&&      & \\rho_0=1.225;\\ H_0=7200;\\ C_d=0.5;\\ A^{\\text{ref}}=4\\pi;\\\\\n",
    "&&      & T^0=4854100;\\ T^1=2968600;\\ T^2=1083100; T^3=110094;\\\\\n",
    "&&      & \\dot{m}^0=1723.273;\\ \\dot{m}^1=1044.682;\\ \\dot{m}^2=366.092;\\ \\dot{m}^3=24.029;\\\\\n",
    "&\\quad \\textit{Dynamics:} &      & \\dot{\\mathbf{r}} =  \\mathbf{v}\\\\\n",
    "&                  &      & \\dot{\\mathbf{v}} =  -\\dfrac{\\mu}{|{\\mathbf{r}|}^3}\\mathbf{r} + \\dfrac{T^p}{m}\\mathbf{u} - \\dfrac{C_dA^{\\text{ref}}\\rho_0|{\\mathbf{v}-(\\omega_e\\times \\mathbf{r})}|(\\mathbf{v}-(\\omega_e\\times \\mathbf{r}))e^{\\frac{|{\\mathbf{r}}|-R_e}{H_0}}}{2m}\\\\\n",
    "& &      & \\dot{m} = -\\dot{m}^p \\qquad \\qquad p = 0, 1, 2, 3\\\\\n",
    "& \\quad \\textit{Path constraints: } &      & |{\\mathbf{u}}| = 1 \\\\\n",
    "&                  &      & |{\\mathbf{r}}| \\geq R_e \\\\\n",
    "& \\quad \\textit{Terminal constraints: } & & t_0^0 =0;\\ t_0^1 = 75.2;\\ t_0^2=150.4;\\ t_0^3 = 261   \\\\\n",
    "&                  & & t_f^0 =75.2;\\ t_f^1 = 150.4;\\ t_f^2=261 \\\\\n",
    "&&      & \\text{Let}\\ \\mathbf{r_f} = [x_0(t_f^3), x_1(t_f^3), x_2(t_f^3)]; \\ \\mathbf{v_f} = [x_3(t_f^3), x_4(t_f^3), x_5(t_f^3)];\\\\\n",
    "&                  & & a_f(\\mathbf{r_f}, \\mathbf{v_f}) = 24361140; e_f(\\mathbf{r_f}, \\mathbf{v_f}) = 0.7308; \\\\\n",
    "&                  & & i_f(\\mathbf{r_f}, \\mathbf{v_f}) = \\frac{28.5\\pi}{180}; \\Omega_f(\\mathbf{r_f}, \\mathbf{v_f}) = \\frac{269.8\\pi}{180}; \\\\\n",
    "&                  & & \\omega_f(\\mathbf{r_f}, \\mathbf{v_f}) = \\frac{130.5\\pi}{180};  \\\\\n",
    "&                  & & \\text{Where}\\ a_f, e_f, i_f, \\Omega_f, \\omega_f\\ \\text{are computed as follows}\\\\\n",
    "&                  & & \\mathbf{h} = \\mathbf{r_f}\\times \\mathbf{v_f};\\ \\mathbf{n} = [0, 0, 1]^T\\times \\mathbf{v_f};\\\\\n",
    "&                  & & \\mathbf{e} = \\frac{1}{\\mu}\\mathbf{v_f}\\times \\mathbf{h} - \\dfrac{\\mathbf{r_f}}{|{\\mathbf{r_f}}|};\\\\\n",
    "&                  & &e_f = |{\\mathbf{e}}|;\\ a_f = -\\dfrac{\\mu}{|{\\mathbf{v_f}}|^2-\\frac{2\\mu}{|{\\mathbf{r_f}}|}};\\ i_f = \\cos^{-1}\\left(\\dfrac{\\mathbf{h}.[0, 0, 1]^T}{|{\\mathbf{h}}|}\\right);\\\\\n",
    "&                  & & \\text{if }\\mathbf{n}[1]\\text{>0:}\\ \\Omega_f = \\cos^{-1}\\left(\\dfrac{\\mathbf{n}.[0, 0, 1]^T}{|{\\mathbf{n}}|}\\right)\\ \\text{else:}\\ \\Omega_f = (2\\pi - \\Omega_f);\\\\\n",
    "&                  & & \\text{if e[2]>0:}\\ \\omega_f = \\cos^{-1}\\left(\\dfrac{\\mathbf{n}.\\mathbf{e}}{|{\\mathbf{n}}|e_f}\\right)\\ \\text{else:}\\ \\omega_f = (2\\pi - \\omega_f);\\\\\n",
    "& \\quad \\textit{Event constraints: } &      &  \\mathbf{r}(t_f^p) - \\mathbf{r}(t_0^{p+1}) = 0; \\ \\mathbf{v}(t_f^{p}) - \\mathbf{v}(t_0^{p+1}) = 0; \\quad p = 0, 1, 2\\\\\n",
    "&  &      & m(t_f^p) - m(t_0^{p+1}) = [13680, 6840, 8830]\\quad p = 0, 1, 2\\\\\n",
    "& \\quad \\textit{Initial conditions: } &      &  \\mathbf{r}(t_0^0) = [5605222.973, 0, 3043387.761]; \\mathbf{v}(t_0^{0}) = [0, 408.74, 0];\\\\\n",
    "&  &      & m(t_0^0) = 301454\n",
    "  \\end{aligned}\n",
    "\n",
    "\n",
    "We first create an OCP object and then polulate the object with dynamics, path_constraints, terminal_constraints and objective (running_costs, terminal_costs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "ocp = mp.OCP(n_states=7, n_controls=4, n_phases=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize parameters\n",
    "Re = 6378145.0  # m\n",
    "omegaE = 7.29211585e-5\n",
    "rho0 = 1.225\n",
    "rhoH = 7200.0\n",
    "Sa = 4 * np.pi\n",
    "Cd = 0.5\n",
    "muE = 3.986012e14\n",
    "g0 = 9.80665\n",
    "\n",
    "# Variable initialization\n",
    "lat0 = 28.5 * np.pi / 180.0\n",
    "r0 = np.array([Re * np.cos(lat0), 0.0, Re * np.sin(lat0)])\n",
    "# v0     = omegaE*np.array([-r0[1], r0[0], 0.0])\n",
    "v0 = omegaE * np.array([0.1, 0.1, 0.1])\n",
    "m0 = 431.6e3 + 107.5e3\n",
    "mf = 107.5e3 - 103.5e3\n",
    "mdryBooster = 431.6e3 - 409.5e3\n",
    "mdrySecond = mf\n",
    "x0 = np.array([r0[0], r0[1], r0[2], v0[0], v0[1], v0[2], m0])\n",
    "q_max = 80 * 1e3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step-1 : Define dynamics\n",
    "# Thrust(N) and mass flow rate(kg/s) in each stage\n",
    "Thrust = [9 * 934.0e3, 934.0e3, 934.0e3]\n",
    "\n",
    "\n",
    "def stage_dynamics(x, u, t, param=0, T=0.0):\n",
    "    r = x[:3]\n",
    "    v = x[3:6]\n",
    "    m = x[6]\n",
    "    r_mag = ca.sqrt(r[0] * r[0] + r[1] * r[1] + r[2] * r[2])\n",
    "    v_rel = v  # ca.vertcat(v[0] + r[1] * omegaE, v[1] - r[0] * omegaE, v[2])\n",
    "    v_rel_mag = ca.sqrt(v_rel[0] * v_rel[0] + v_rel[1] * v_rel[1] + v_rel[2] * v_rel[2])\n",
    "    h = r_mag - Re\n",
    "    rho = rho0 * ca.exp(-h / rhoH)\n",
    "    D = -rho / (2 * m) * Sa * Cd * v_rel_mag * v_rel\n",
    "    g = -muE / (r_mag * r_mag * r_mag) * r\n",
    "\n",
    "    xdot = [\n",
    "        x[3],\n",
    "        x[4],\n",
    "        x[5],\n",
    "        T * u[3] / m * u[0] + param * D[0] + g[0],\n",
    "        T * u[3] / m * u[1] + param * D[1] + g[1],\n",
    "        T * u[3] / m * u[2] + param * D[2] + g[2],\n",
    "        -T * u[3] / (340.0 * g0),\n",
    "    ]\n",
    "    return xdot\n",
    "\n",
    "\n",
    "def get_dynamics(param):\n",
    "    dynamics0 = lambda x, u, t: stage_dynamics(x, u, t, param=param, T=Thrust[0])\n",
    "    dynamics1 = lambda x, u, t: stage_dynamics(x, u, t, param=param, T=Thrust[1])\n",
    "    dynamics2 = lambda x, u, t: stage_dynamics(x, u, t, param=param, T=Thrust[2])\n",
    "\n",
    "    return [dynamics0, dynamics1, dynamics2]\n",
    "\n",
    "\n",
    "ocp.dynamics = get_dynamics(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step-2: Add path constraints\n",
    "def path_constraints0(x, u, t):\n",
    "    return [\n",
    "        u[0] * u[0] + u[1] * u[1] + u[2] * u[2] - 1,\n",
    "        -u[0] * u[0] - u[1] * u[1] - u[2] * u[2] + 1,\n",
    "        -ca.sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]) / Re + 1,\n",
    "    ]\n",
    "\n",
    "\n",
    "def path_constraints2(x, u, t, dynP=0, gs=0):\n",
    "    r_mag = ca.sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2])\n",
    "    h = r_mag - Re\n",
    "    rho = rho0 * ca.exp(-h / rhoH)\n",
    "    v_sq = x[3] * x[3] + x[4] * x[4] + x[5] * x[5]\n",
    "    r_rf = ca.vertcat(x[0] - x0[0], x[1] - x0[1], x[2] - x0[2])\n",
    "    r_rf_mag = ca.sqrt(r_rf[0] * r_rf[0] + r_rf[1] * r_rf[1] + r_rf[2] * r_rf[2])\n",
    "    rf_mag = np.sqrt(x0[0] * x0[0] + x0[1] * x0[1] + x0[2] * x0[2])\n",
    "    glide_slope_factor = np.cos(80.0 * np.pi / 180.0)\n",
    "\n",
    "    return [\n",
    "        dynP * 0.5 * rho * v_sq / q_max - 1.0,\n",
    "        u[0] * u[0] + u[1] * u[1] + u[2] * u[2] - 1,\n",
    "        -u[0] * u[0] - u[1] * u[1] - u[2] * u[2] + 1,\n",
    "        -ca.sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]) / Re + 1,\n",
    "        gs\n",
    "        * (\n",
    "            r_rf_mag * glide_slope_factor\n",
    "            - (r_rf[0] * x0[0] + r_rf[1] * x0[1] + r_rf[2] * x0[2]) / rf_mag\n",
    "        ),\n",
    "    ]\n",
    "\n",
    "\n",
    "ocp.path_constraints = [path_constraints0, path_constraints0, path_constraints2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step-3: Add terminal cost and constraints\n",
    "def terminal_cost1(xf, tf, x0, t0):\n",
    "    return -xf[6] / m0\n",
    "\n",
    "\n",
    "ocp.terminal_costs[1] = terminal_cost1\n",
    "\n",
    "\n",
    "def terminal_constraints1(x, t, x0, t0):\n",
    "    h = ca.vertcat(\n",
    "        x[1] * x[5] - x[4] * x[2], x[3] * x[2] - x[0] * x[5], x[0] * x[4] - x[1] * x[3]\n",
    "    )\n",
    "\n",
    "    n = ca.vertcat(-h[1], h[0], 0)\n",
    "    r = ca.sqrt(x[0] * x[0] + x[1] * x[1] + x[2] * x[2])\n",
    "\n",
    "    e = ca.vertcat(\n",
    "        1 / muE * (x[4] * h[2] - x[5] * h[1]) - x[0] / r,\n",
    "        1 / muE * (x[5] * h[0] - x[3] * h[2]) - x[1] / r,\n",
    "        1 / muE * (x[3] * h[1] - x[4] * h[0]) - x[2] / r,\n",
    "    )\n",
    "\n",
    "    e_mag = ca.sqrt(e[0] * e[0] + e[1] * e[1] + e[2] * e[2])\n",
    "    h_sq = h[0] * h[0] + h[1] * h[1] + h[2] * h[2]\n",
    "    v_mag = ca.sqrt(x[3] * x[3] + x[4] * x[4] + x[5] * x[5])\n",
    "\n",
    "    a = -muE / (v_mag * v_mag - 2.0 * muE / r)\n",
    "    i = ca.acos(h[2] / ca.sqrt(h_sq))\n",
    "    n_mag = ca.sqrt(n[0] * n[0] + n[1] * n[1])\n",
    "\n",
    "    node_asc = ca.acos(n[0] / n_mag)\n",
    "    # if n[1] < -1e-12:\n",
    "    node_asc = 2 * np.pi - node_asc\n",
    "\n",
    "    argP = ca.acos((n[0] * e[0] + n[1] * e[1]) / (n_mag * e_mag))\n",
    "    # if e[2] < 0:\n",
    "    #    argP = 2*np.pi - argP\n",
    "\n",
    "    a_req = 6593145.0  # 24361140.0\n",
    "    e_req = 0.0076\n",
    "    i_req = 28.5 * np.pi / 180.0\n",
    "    node_asc_req = 269.8 * np.pi / 180.0\n",
    "    argP_req = 130.5 * np.pi / 180.0\n",
    "\n",
    "    return [\n",
    "        (a - a_req) / (Re),\n",
    "        e_mag - e_req,\n",
    "        i - i_req,\n",
    "        node_asc - node_asc_req,\n",
    "        argP - argP_req,\n",
    "    ]\n",
    "\n",
    "\n",
    "def terminal_constraints2(x, t, x_0, t_0):\n",
    "    return [\n",
    "        (x[0] - x0[0]) / Re,\n",
    "        (x[1] - x0[1]) / Re,\n",
    "        (x[2] - x0[2]) / Re,\n",
    "        (x[3] - x0[3]) / np.sqrt(muE / Re),\n",
    "        (x[4] - x0[4]) / np.sqrt(muE / Re),\n",
    "        (x[5] - x0[5]) / np.sqrt(muE / Re),\n",
    "    ]\n",
    "\n",
    "\n",
    "ocp.terminal_constraints[1] = terminal_constraints1\n",
    "ocp.terminal_constraints[2] = terminal_constraints2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scale the variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ocp.scale_x = np.array(\n",
    "    [\n",
    "        1 / Re,\n",
    "        1 / Re,\n",
    "        1 / Re,\n",
    "        1 / np.sqrt(muE / Re),\n",
    "        1 / np.sqrt(muE / Re),\n",
    "        1 / np.sqrt(muE / Re),\n",
    "        1 / m0,\n",
    "    ]\n",
    ")\n",
    "ocp.scale_t = np.sqrt(muE / Re) / Re"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initial state and Final guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial guess estimation\n",
    "# User defined functions\n",
    "def ae_to_rv(a, e, i, node, argP, th):\n",
    "    p = a * (1.0 - e * e)\n",
    "    r = p / (1.0 + e * np.cos(th))\n",
    "\n",
    "    r_vec = np.array([r * np.cos(th), r * np.sin(th), 0.0])\n",
    "    v_vec = np.sqrt(muE / p) * np.array([-np.sin(th), e + np.cos(th), 0.0])\n",
    "\n",
    "    cn, sn = np.cos(node), np.sin(node)\n",
    "    cp, sp = np.cos(argP), np.sin(argP)\n",
    "    ci, si = np.cos(i), np.sin(i)\n",
    "\n",
    "    R = np.array(\n",
    "        [\n",
    "            [cn * cp - sn * sp * ci, -cn * sp - sn * cp * ci, sn * si],\n",
    "            [sn * cp + cn * sp * ci, -sn * sp + cn * cp * ci, -cn * si],\n",
    "            [sp * si, cp * si, ci],\n",
    "        ]\n",
    "    )\n",
    "\n",
    "    r_i = np.dot(R, r_vec)\n",
    "    v_i = np.dot(R, v_vec)\n",
    "\n",
    "    return r_i, v_i\n",
    "\n",
    "\n",
    "# Target conditions\n",
    "a_req = 6593145.0  # 24361140.0\n",
    "e_req = 0.0076\n",
    "i_req = 28.5 * np.pi / 180.0\n",
    "node_asc_req = 269.8 * np.pi / 180.0\n",
    "argP_req = 130.5 * np.pi / 180.0\n",
    "th = 0.0\n",
    "rf, vf = ae_to_rv(a_req, e_req, i_req, node_asc_req, argP_req, th)\n",
    "xf = np.array([rf[0], rf[1], rf[2], vf[0], vf[1], vf[2], mf])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Intial guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Timings\n",
    "# Timings\n",
    "t0, t1, t2, t3 = 0.0, 131.4, 453.4, 569.7\n",
    "\n",
    "# Interpolate to get starting values for intermediate phases\n",
    "x1 = x0 + (xf - x0) / (t2 - t0) * (t1 - t0)\n",
    "\n",
    "# Update the state discontinuity values across phases\n",
    "x0f = np.copy(x1)\n",
    "x0f[-1] = x0[-1] - (9 * 934e3 / (340.0 * g0) * t1)\n",
    "mFirst_leftout = 409.5e3 - (9 * 934e3 / (340.0 * g0) * t1)\n",
    "x1[-1] = x0f[-1] - (mdryBooster + mFirst_leftout)\n",
    "\n",
    "# Step-8b: Initial guess for the states, controls and phase start and final\n",
    "#             times\n",
    "ocp.x00 = np.array([x0, x1, x0f])\n",
    "ocp.xf0 = np.array([x0f, xf, x0])\n",
    "ocp.u00 = np.array([[1, 0, 0, 1.0], [1, 0, 0, 1], [0, 1, 0, 1]])\n",
    "ocp.uf0 = np.array([[0, 1, 0, 1.0], [0, 1, 0, 1], [1, 0, 0, 0.5]])\n",
    "ocp.t00 = np.array([[t0], [t1], [t1]])\n",
    "ocp.tf0 = np.array([[t1], [t2], [t3]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Box constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step-8c: Bounds for states\n",
    "rmin, rmax = -2 * Re, 2 * Re\n",
    "vmin, vmax = -10000.0, 10000.0\n",
    "\n",
    "lbx0 = [rmin, rmin, rmin, vmin, vmin, vmin, x0f[-1]]\n",
    "lbx1 = [rmin, rmin, rmin, vmin, vmin, vmin, xf[-1]]\n",
    "lbx2 = [rmin, rmin, rmin, vmin, vmin, vmin, mdryBooster]\n",
    "ubx0 = [rmax, rmax, rmax, vmax, vmax, vmax, x0[-1]]\n",
    "ubx1 = [rmax, rmax, rmax, vmax, vmax, vmax, 107.5e3]\n",
    "ubx2 = [rmax, rmax, rmax, vmax, vmax, vmax, x0f[-1] - 107.5e3]\n",
    "\n",
    "ocp.lbx = np.array([lbx0, lbx1, lbx2])\n",
    "ocp.ubx = np.array([ubx0, ubx1, ubx2])\n",
    "\n",
    "# Bounds for control inputs\n",
    "ocp.lbu = np.array(\n",
    "    [[-1.0, -1.0, -1.0, 1.0], [-1.0, -1.0, -1.0, 1.0], [-1.0, -1.0, -1.0, 0.38]]\n",
    ")\n",
    "ocp.ubu = np.array([[1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0]])\n",
    "\n",
    "# Bounds for phase start and final times\n",
    "ocp.lbt0 = np.array([[t0], [t1], [t1]])\n",
    "ocp.ubt0 = np.array([[t0], [t1], [t1]])\n",
    "ocp.lbtf = np.array([[t1], [t2 - 50], [t3 - 100]])\n",
    "ocp.ubtf = np.array([[t1], [t2 + 50], [t3 + 100]])\n",
    "\n",
    "# Event constraint bounds on states : State continuity/disc.\n",
    "lbe0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -(mdryBooster + mFirst_leftout)]\n",
    "lbe1 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -107.5e3]\n",
    "ocp.lbe = np.array([lbe0, lbe1])\n",
    "ocp.ube = np.array([lbe0, lbe1])\n",
    "\n",
    "ocp.phase_links = [(0, 1), (0, 2)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "ocp.validate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve with atmospheric drag disabled "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "******************************************************************************\n",
      "This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
      " Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
      "         For more information visit http://projects.coin-or.org/Ipopt\n",
      "******************************************************************************\n",
      "\n",
      "Total number of variables............................:      956\n",
      "                     variables with only lower bounds:        0\n",
      "                variables with lower and upper bounds:      956\n",
      "                     variables with only upper bounds:        0\n",
      "Total number of equality constraints.................:      746\n",
      "Total number of inequality constraints...............:      641\n",
      "        inequality constraints with only lower bounds:        0\n",
      "   inequality constraints with lower and upper bounds:      300\n",
      "        inequality constraints with only upper bounds:      341\n",
      "\n",
      "\n",
      "Number of Iterations....: 153\n",
      "\n",
      "                                   (scaled)                 (unscaled)\n",
      "Objective...............:  -2.9204480860707361e-02   -2.9204480860707361e-02\n",
      "Dual infeasibility......:   1.9040474700928261e-07    1.9040474700928261e-07\n",
      "Constraint violation....:   6.7175836255209499e-06    6.7175836255209499e-06\n",
      "Complementarity.........:   2.5059035596800655e-09    2.5059035596800655e-09\n",
      "Overall NLP error.......:   6.7175836255209499e-06    6.7175836255209499e-06\n",
      "\n",
      "\n",
      "Number of objective function evaluations             = 203\n",
      "Number of objective gradient evaluations             = 154\n",
      "Number of equality constraint evaluations            = 203\n",
      "Number of inequality constraint evaluations          = 203\n",
      "Number of equality constraint Jacobian evaluations   = 154\n",
      "Number of inequality constraint Jacobian evaluations = 154\n",
      "Number of Lagrangian Hessian evaluations             = 153\n",
      "Total CPU secs in IPOPT (w/o function evaluations)   =      2.930\n",
      "Total CPU secs in NLP function evaluations           =      0.117\n",
      "\n",
      "EXIT: Solved To Acceptable Level.\n",
      "      solver  :   t_proc      (avg)   t_wall      (avg)    n_eval\n",
      "       nlp_f  | 744.00us (  3.67us) 720.91us (  3.55us)       203\n",
      "       nlp_g  |  28.51ms (140.46us)  28.24ms (139.09us)       203\n",
      "    nlp_grad  | 384.00us (384.00us) 383.90us (383.90us)         1\n",
      "  nlp_grad_f  |   1.51ms (  9.72us)   1.41ms (  9.13us)       155\n",
      "  nlp_hess_l  |  34.36ms (224.61us)  34.18ms (223.38us)       153\n",
      "   nlp_jac_g  |  41.84ms (269.92us)  41.74ms (269.30us)       155\n",
      "       total  |   3.07 s (  3.07 s)   3.06 s (  3.06 s)         1\n"
     ]
    }
   ],
   "source": [
    "ocp.dynamics = get_dynamics(0)\n",
    "mpo = mp.mpopt(ocp, 5, 6)\n",
    "sol = mpo.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve the problem with drag enabled now with revised initial guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of variables............................:      956\n",
      "                     variables with only lower bounds:        0\n",
      "                variables with lower and upper bounds:      956\n",
      "                     variables with only upper bounds:        0\n",
      "Total number of equality constraints.................:      746\n",
      "Total number of inequality constraints...............:      641\n",
      "        inequality constraints with only lower bounds:        0\n",
      "   inequality constraints with lower and upper bounds:      300\n",
      "        inequality constraints with only upper bounds:      341\n",
      "\n",
      "\n",
      "Number of Iterations....: 61\n",
      "\n",
      "                                   (scaled)                 (unscaled)\n",
      "Objective...............:  -3.1423421485119993e-02   -3.1423421485119993e-02\n",
      "Dual infeasibility......:   6.7634340985309009e-12    6.7634340985309009e-12\n",
      "Constraint violation....:   1.9999893324795792e-08    1.9999893324795792e-08\n",
      "Complementarity.........:   9.1006695117999755e-10    9.1006695117999755e-10\n",
      "Overall NLP error.......:   1.9999893324795792e-08    1.9999893324795792e-08\n",
      "\n",
      "\n",
      "Number of objective function evaluations             = 172\n",
      "Number of objective gradient evaluations             = 63\n",
      "Number of equality constraint evaluations            = 172\n",
      "Number of inequality constraint evaluations          = 172\n",
      "Number of equality constraint Jacobian evaluations   = 63\n",
      "Number of inequality constraint Jacobian evaluations = 63\n",
      "Number of Lagrangian Hessian evaluations             = 62\n",
      "Total CPU secs in IPOPT (w/o function evaluations)   =      1.441\n",
      "Total CPU secs in NLP function evaluations           =      0.079\n",
      "\n",
      "EXIT: Solved To Acceptable Level.\n",
      "      solver  :   t_proc      (avg)   t_wall      (avg)    n_eval\n",
      "       nlp_f  | 539.00us (  3.13us) 523.99us (  3.05us)       172\n",
      "       nlp_g  |  25.06ms (145.73us)  24.79ms (144.13us)       172\n",
      "    nlp_grad  | 305.00us (305.00us) 303.70us (303.70us)         1\n",
      "  nlp_grad_f  | 579.00us (  9.05us) 660.80us ( 10.33us)        64\n",
      "  nlp_hess_l  |  25.98ms (419.10us)  25.94ms (418.32us)        62\n",
      "   nlp_jac_g  |  23.84ms (372.56us)  23.79ms (371.73us)        64\n",
      "       total  |   1.55 s (  1.55 s)   1.54 s (  1.54 s)         1\n"
     ]
    }
   ],
   "source": [
    "ocp.dynamics = get_dynamics(1)\n",
    "ocp.path_constraints[2] = lambda x, u, t: path_constraints2(x, u, t, dynP=1, gs=0)\n",
    "ocp.validate()\n",
    "\n",
    "mpo._ocp = ocp\n",
    "sol = mpo.solve(\n",
    "    sol, reinitialize_nlp=True, nlp_solver_options={\"ipopt.acceptable_tol\": 1e-6}\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Retrive data from the solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Payload:\n",
      "Mass to orbit :  16940.3665 kg(MPOPT) vs 17310kg(Beigler paper)\n",
      "Time to orbit (SECO) :  454.6864 s(MPOPT) vs 453.4s(Beigler paper)\n",
      "\n",
      "Booster:\n",
      " Mass at landing :  22100.0 kg(MPOPT) vs 22100kg(Beigler paper)\n",
      "Time of landing :  579.5091 s(MPOPT) vs 569.7s(Beigler paper)\n"
     ]
    }
   ],
   "source": [
    "post = mpo.process_results(sol, plot=False, scaling=False)\n",
    "mp.post_process._INTERPOLATION_NODES_PER_SEG = 200\n",
    "x, u, t, _ = post.get_data(phases=ocp.phase_links[0], interpolate=False)\n",
    "print(\"Payload:\\nMass to orbit : \", round(x[-1][-1], 4), \"kg(MPOPT) vs 17310kg(Beigler paper)\")\n",
    "print(\"Time to orbit (SECO) : \", round(t[-1][0], 4), \"s(MPOPT) vs 453.4s(Beigler paper)\")\n",
    "x, u, t, _ = post.get_data(phases=ocp.phase_links[1], interpolate=False)\n",
    "print(\"\\nBooster:\\n Mass at landing : \", round(x[-1][-1], 4), \"kg(MPOPT) vs 22100kg(Beigler paper)\")\n",
    "print(\"Time of landing : \", round(t[-1][0], 4), \"s(MPOPT) vs 569.7s(Beigler paper)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot trajectory of launch to orbit (Second stage) and Booster return to Ship."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Text(0, 0.5, 'Velocity, km/s')]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtkAAAHgCAYAAABw0HFmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeVyVZf7/8dfFOayCiEKIogJu4AKouGWamrZoqWWb1WQ1jWUu35mWaR+bX8u02CypTWPbmFZj2lSWaVmJYO77BpoICCoKKAiynO36/QE6mmYC53Cfc/g8H/GQc3Pu+357oMsP17kWpbVGCCGEEEII4Tw+RgcQQgghhBDC20iRLYQQQgghhJNJkS2EEEIIIYSTSZEthBBCCCGEk0mRLYQQQgghhJNJkS2EEEIIIYSTmY0O0BDh4eE6JiamzuedOnWKZs2aOT9QI5H8xvHk7CD5jfbz/Js3by7SWkcYGKnR1afd9rbvu6eR/Mbx5Ozgffnr3GZrrT32o0+fPro+Vq5cWa/z3IXkN44nZ9da8hvt5/mBTdoN2tLG/KhPu+1t33dPI/mN48nZtfa+/HVts2W4iBBCCCGEEE4mRbYQQgghhBBOJkW2EEIIIYQQTubREx8vxGq1kp+fT1VV1S8+JzQ0lIyMjEZM5VxNIX9AQADR0dH4+vo2UiohhFF+rd1uCm2eO5M2W4j68boiOz8/n5CQEGJiYlBKXfA5ZWVlhISENHIy5/H2/FpriouLyc/PJzY2thGTCSGM8Gvttre3ee5O2mwh6sfrhotUVVXRqlWrXyywhftTStGqVauLvhshhHAPSqk/KKV2K6V2KaU+VkoF1PUa0m57NmmzhbgwryuyAWmovYB8D4Vwf0qptsB0IEVr3QMwAbfX81rOjCYamXz/hDifVxbZRnvxxRfp3r07iYmJJCcns379+kbPkJqayvXXX39Jzy0pKeHNN990caKatxSnT59OUlISiYmJbNmyxeX3FEK4nBkIVEqZgSDgsMF56kXa7QvLzMzkqquuwt/fn5kzZ7r8fkJ4E68bk220tWvX8tVXX7Flyxb8/f0pKirCYrEYHeuiTjfWDz30kEvvs2zZMn766Se2bdvGnj17mDx5siH/kAkhnENrfUgpNRM4CFQC32qtvzU4Vp1Ju/3LWrZsyauvvsqKFStceh8hvJEU2QB5GyAnHWIGQ7t+DbrUkSNHCA8Px9/fH4Dw8PAzX9u8eTMPP/ww5eXlhIeH8+9//5uoqCj279/Pgw8+SGFhISaTiUWLFhEXF8cf//hHli1bhlKKZ555httuu43U1FSeffZZIiMj2bVrF3369GHBggUopVi+fDm///3vCQoK4oorrrhgvt27d3PvvfdisVhwOBx8+umnPPvss2RlZZGcnMzIkSOZMWMGY8eO5cSJE1itVl544QXGjh0LwPPPP8+CBQuIiIigXbt29OnTh0cffZSsrCymTJlCYWEhQUFBvP3228THx59z7y+++IK7774bpRQDBgygpKSEI0eOEBUV1aDXXAhhDKVUGDAWiAVKgEVKqbu01gt+9rxJwCSAyMhIUlNTz7lOaGgoZWVlv3gfu91+3td9Dm/GnLcWW7uBONr0adDf48CBA7Ro0QKLxYLFYsHf3x9/f3/KysrYunUrTz31FKdOnaJly5a89dZbtG7dmqysLP7whz9QVFSEyWRi3rx5xMbG8uyzz7JixQqUUjz22GOMHz+eVatW8corr9CqVSv27NlDcnIy77zzDkopVqxYwRNPPEFQUBADBgzAZrOd93fNyMhg8uTJWK1WHA4H8+fP54UXXiArK4vExESGDRvGE088wYQJEygpKcFqtfLss88yevRoAF555RUWLlxIeHg4bdu2pVevXkyfPp0DBw7wyCOPUFxcTGBgILNmzaJLly7n3DswMJDk5GS++eYbqqurL/p9qqqqOu976w7Ky8vdMtel8OTsIPkN32K3IR8X2p53z549v7pN5smTJ//34OB6rZ+P1Pq5sJo/D67/1fMvpqysTCclJenOnTvryZMn69TUVK211haLRQ8cOFAfO3ZMa631f/7zH33vvfdqrbXu16+f/u9//6u11rqyslKfOnVKL168WI8YMULbbDZdUFCg27Vrpw8fPqxXrlypmzdvrvPy8rTdbtcDBgzQ6enpurKyUkdHR+t9+/Zph8Ohb7nlFj169Ojz8k2dOlUvWLBAa611dXW1rqio0NnZ2bp79+5nnmO1WnVpaanWWuvCwkLdsWNH7XA49IYNG3RSUpKurKzUJ0+e1J06ddKvvfaa1lrr4cOH63379mmttV63bp0eNmzYefcePXq0Tk9PP/P6Dx8+XG/cuPEXX8tL+V42Nm/bItbTeFt+PHxbdeAW4N2zHt8NvHmxc+rTbp/TZmvtce320qVLPbbd1rrm9Z8xY8aZ836JO7bZWnt2u+HJ2bX2vvx1bbOlJzsnHewW0PaaP3PSG9SbHRwczObNm0lPT2flypXcdtttvPzyy6SkpLBr1y5GjhwJ1PTMREVFUVZWxqFDh7jxxhuBmrVGAVavXs2ECRMwmUxERkZy5ZVXsnHjRpo3b06fPn2Ijo4GIDk5mZycHIKDg4mNjaVz584A3HXXXcydO/e8fAMHDuTFF18kPz+fm2666czzz6a15qmnniItLQ0fHx8OHTrE0aNH+fHHHxk7diwBAQEEBARwww03ADW/6a1Zs4ZbbrnlzDWqq6vr/RoKITzGQWCAUiqImuEiVwGbXH5XD2u3zWYz/fr1k3ZbiCZGiuyYwWDyq2moTX41jxvIZDIxdOhQhg4dSs+ePZk3bx59+vShe/furF279pznXuytt1/i5+d3zr1sNtsln3vHHXfQv39/li5dyqhRo/jXv/5FXFzcOc/58MMPKSwsZPPmzfj6+hITE3PRpZkcDgctWrRg27ZtF71327ZtycvLIykpCahZG7dt27aXnF0I4V601uuVUouBLYAN2AqcXyU6mwe226eHEJ6+l6e020KI+pPVRdr1g4lLYPjTNX82cEz23r17+emnn8483rZtGx06dKBr164UFhaeaaytViu7d+8mJCSE6OhoPv/8c6CmJ6GiooLBgwezcOFC7HY7hYWFpKWl0a/fL2eLj48nJyeHrKwsAD7++OMLPu/AgQPExcUxffp0xo4dy44dOwgJCTnnH43S0lIuu+wyfH19WblyJbm5uQAMGjSIL7/8kqqqKsrLy/nqq68AaN68ObGxsSxatAio6VHZvn37efceM2YMH3zwAVpr1q1bR2hoqIzHFsLDaa1naK3jtdY9tNa/0Vq7vjtU2u1Ga7eFEPUnPdlQ00A3sJE+rby8nGnTplFSUoLZbKZTp07MnTsXPz8/Fi9ezPTp0yktLcVms/H73/+e7t27M3/+fB544AH+9Kc/4evry6JFi7jxxhtZu3YtSUlJKKV49dVXad26NZmZmRe8b0BAAHPnzmX06NEEBQUxePDgC/a2fPLJJ8yfPx9fX19at27NU089RcuWLRk0aBA9evTguuuu4/HHH+eGG26gZ8+epKSknJnA2LdvX8aMGUNiYiKRkZH07NmT0NBQoKYXZfLkybzwwgtYrVZuv/32Mz3Wp40aNYqvv/6apKQkgoODef/9953ymgshmiAPard/iSe02wUFBfTp04eysjJ8fHz4+9//zp49e2jevLlTXnshvFpdBnC724dTJj56ICPzl5WVaa21PnXqlO7Tp4/evHlzna9xqfndcRKNt03i8DTelh8Pn/hYnw+nTHz0MEbnb2i77cltttae3W54cnatvS9/Xdts6ckWdTJp0iT27NlDVVUVEydOpHfv3kZHcm95G2D7R4CCpAnQrh8Oh6Yqex2+eT9ijhuMat/f6JRCCC8m7ba4GGvOOqxZaQR1Geq0d4dEDSmyRZ189NFHRkdwb2etuV5ts2OePwYfR82mFtZN87lPz6DCYuNDv5dQ2KhaaeZ36k8UhyXTLiyQhKjmDAnKplv1dgI7D5UGTwjRYNJui1+UtwGf+WPxs1mwr3kd0z1fyr87TiRFthDOkrcBPW8M2laNTfnyX9tgbvWxolTNl83YuKdtPgG+Pvjn2vDBgQ92Jkbl8bH/ALIKyynKTOdB35fwxUZ12ky+6vUWvS6/mriIYGP/bkIIIbxPTjrKbsGkHGiHtcHLYYpzSZEthBMcP2Uh8/vP6WerxowDH22l42XN4IQv2mFBAT4mP0ZcN77mhHnvg92Cj8mPkaNuZmS7vgBUr9yI3yobCgdgJWvDch5Z609KhzBuTWnH9UlRBPnJ/7ZCCCGcIGYwNuWLj7ZictJymOJ/5F9rIRqgwmLj3fRs/pV2gC6WCD4O8MUHGyaTH/3GTQWmnjcmG6hZdqx2WMnZvQb+na6EH18HuwWzyY/7bvsNzY+14ZNNefzx0x28tCyDx7uXMqZFNs1k/JwQQoiGaNePP7V4iUTbTu649U75N8XJpMgWoj7yNvDT+mW8nBnO9+UxXN0tkseuuR9/y6Dzi+cLNVq/tPzY6fV/a68R3q4fDybAA0Pi2Jhzgu9XfMm4HY/iiw3r6plUTviM5l0GufbvKoQQwis5HJqvjrfDv/cAaNfD6DheRzajcQGTyURycjJJSUn07t2bNWvWOPX6M2fOdOr1TktNTXV61gvJycmhf//+dOrUidtuuw2LxeLyezrTyX0/YnnvemJ3/p05tudYdpM/c+9OoXNkSE2RPPiRhvUGXOAaSin6xbbkyYQiAnzsmJUD5bDy3ofzmf3DT1RYLn33OCHE+Vzdbr/00ktOvd5pjdVuz549m06dOqGUoqioyOX3E43jUEkl5dU2uraWdc9dQYrsWoWzZjvtWoGBgWzbto3t27fzl7/8hSeffNJp1wZ4/fXX63yO3W7/1ec0VmM9Y8YM/vCHP7B//37CwsJ49913XX5PZ9lfYufDhR/i47BiVg78lZ2E6kbcJS1mMMrkB8qEj68flW0HMvPbfVz1+irWpn6NTnu9ZoUTIZoAT2q361Nku1O7PWjQIL777js6dOjg8nuJxpNZULP5UXxUiMFJvJMU2bWK5sxxyXVPnjxJWFgYULPxz2OPPUaPHj3o2bMnCxcuvOjxI0eOMGTIEJKTk+nRowfp6ek88cQTVFZWkpyczJ133gnAggUL6NevH8nJyTzwwANnGubg4GAeeeQRkpKSzmwLfNobb7xBt27dSExM5PbbbycnJ4e33nqLv/3tbyQnJ5Oens6XX35J//796dWrFyNGjODo0aMAFBYWMnLkSLp37879999Phw4dzvRs/FKW07TWrFq1iptvvhmAiRMnntma2N19sDaHv6yvYpupBz7mmkJXNfZEkbO2k/aZ+CVPTprIogcHMsA3i+SVE3H88AKOeTdIoS2aBE9pt2fMmOHR7TZAr169iImJcf6LLQyVeeQkAF0ipch2ibrsXONuH87c8XFP1/hfPe9S+fj46KSkJN21a1fdvHlzvWnTJq211osXL9YjRozQNptNFxQU6Hbt2unDhw//4vGZM2fqF154QWuttc1mO5O7WbNm/8u9Z4++/vrrtcVi0VprPXnyZD1v3jyttdaAXrhw4QUzRkVF6aqqKq211idOnNBaaz1jxgz92muvnXnO8ePHtcPh0Fpr/fbbb+uHH35Ya631lClT9EsvvaS11nrZsmUa0IWFhRfNclphYaGOjY098/jgwYO6e/fuF8zoLruH2e0O/eclu3WHx7/SY2Yu0yWnLFofXK912syaP92ALfU1bZ/RQusZzbX1Ty30zrfu1bbU187L5227b3ka2fHRuTs+ekq7ffLkSY9ut89+/Tt06KALCwsvmM9d2uyf8+R2w9XZH1qwWQ9+5QeXXd+TX3utZcfHBimcNfucnpCM+AQAwqdMIWLa1Hpf9/TbjgBr167l7rvvZteuXaxevZoJEyZgMpmIjIzkyiuvZOPGjb94vG/fvtx3331YrVbGjRtHcnLyeff6/vvv2bx5M3371iwBV1lZyWWXXQbUjDEcP378BTMmJiZy5513Mm7cOMaNG3fB5+Tn53Pbbbdx5MgRLBYLsbGxAKxevZrPPvsMgGuvvfZMj8/FsniqKqudRz7ZztKdR7h3UAyDg48RGuQLQb8wcdEgprghkD4Tbbeg8aHz4S/gyGc40l7Bp9dd565sIoQH88R2u6ys7Jx7Sbst3EVmwUniW0svtqs06eEiEdOmkpCZQUJmBsCZzxvSUP/cwIEDKSoqorCwsM7nDhkyhLS0NNq2bcs999zDBx98cN5ztNZMnDiRbdu2sW3bNvbu3ctzzz0HQEBAACaT6YLXXrp0KVOmTGHLli307dsXm+38iXPTpk1j6tSp7Ny5k3/9619UVVVdNO/FspzWqlUrSktLz9wvPz+ftm3bXsKr0fiqrHYemL+ZpTuP8MzoBGbc0B2f0zvLuJvaYSRq+NOY+9yFr3JgwoGyWdCb3kfPGyNDSIRXkHa78dtt4Z0qLXayi05Jke1CTbrIbgyZmZnY7XZatWrF4MGDWbhwIXa7ncLCQtLS0ujXr98vHs/NzSUyMpLf/e533H///WzZsgUAX19frFYrAFdddRWLFy/m2LFjABw/fpzc3NyLZnI4HOTl5TFs2DBeeeUVSktLKS8vJyQk5Jwel9LS0jMF8Lx5884cHzRoEJ988gkA3377LSdOnLjkLEophgwZwuLFi89cd+zYsfV7cV3IkrOWL+c8StlPP/LK+J7cPzjO6Ei/rnZVEpV8Bz5mfzSKmv80Dls1JzNWGp1QCI8g7fbFswjvkFlwEoeGbm1CjY7itZr0cJGzhU+Z4rRrnZ7gAjW9BPPmzcNkMnHjjTeydu1akpKSUErx6quv0rp16188Pm/ePF577TV8fX0JDg4+0yNyzz33kJiYSO/evfnwww954YUXuPrqq3E4HPj6+jJnzpyLzgC32+3cddddlJaWorVm+vTptGjRghtuuIGbb76ZL774glmzZvHcc89xyy23EBYWxvDhw8nOzgZqVgeZMGEC8+fPZ+DAgbRu3ZqQkBDCw8MvKcuf//xn7r//fp555hl69erFb3/7W6e99s5gy10H88Zyo8PKjYF+mFsPAtobHevSne7V3v4ReutHOOw2LNrElDVBXNdhN5g2n7cJjhCeyJPa7UmTJnl0u/3GG2/w6quvUlBQQGJiIqNGjeKdd95x2usvGt/uwzWTHru3keX7XKYuA7jd7cOZEx89idH5q6qqtNVq1VprvWbNGp2UlFSn8y81vxGTaBwOh1725qPa+qeaSYT6ubCaCY61PG4SR+0EzfwdK/Ujr8/VFX8K17YZLbTj+Ui3mbRZFx73+v+MTHx07sRHT+EO+RvSbrtzm30pPLndcGX2Jz7doROf++bMRFlX8OTXXmuZ+CgMcPDgQW699VYcDgd+fn68/fbbRkdymtk/7GflwTZcFegL2gaNvUSfs9XuLNkWePn4JtRKGyYc2GzVHN/5HZdJb7YQTYI3t9uifvYcOUm3qOYod51r5AWkyBZ11rlzZ7Zu3Wp0DKf7fOshXl+xj5t6DcM88HLIXe1VwyrMcUOwr3oFh7Zj0yamrmnGqNBsJl4eI42sEF7OW9ttUT+23HUMKfiAlt2HAwOMjuO1XFZkK6XeA64Hjmmte9QeWwh0rX1KC6BEa52slIoBMoC9tV9bp7V+0FXZhPi5XYdKefzTHfSPbcnL4xNRZh9o39/oWM7Vrh/bk56nd8tTVEX2p9mPfjz35R7Sfyri1ZsTaRXsb3RCIYQQrpa3AZ/5Y/k/Hwvs+xzyYr2mM8nduLIn+9/AbODM+kVa69tOf66Ueh0oPev5WVrr8xeCrgettfTMebiaoU+No6TCwuQPNxMW5MecO3vjZ/beRXdOhsbD4KG0AN7rrPn3mhz+8nUm1/0jnb/emswVncONjiiaKGm3PVtjttmigXLSwWbBrBxoba15LEW2S7ismtBapwHHL/Q1VdOS3gp87Oz7BgQEUFxcLP/DezCtNcXFxQQEBLj8Xg6H5v/+s42C0irevKs34U2oN1cpxb2DYvl8yiBCAsz85r31fPDJJ9hWzZQ1tUWjknbbszVmmy2cIGYwduWLTft4/rwjN2fUmOzBwFGt9U9nHYtVSm0FTgLPaK3T63Ph6Oho8vPzL7qJQFVVlUc3Bk0hf0BAANHR0S7PMjf9AKv2FfLCuB70bh/m8vu5o25tmvPVtMH8e+En3LJ7OmqPDUf6a/hM/FJ6N0Sj+LV2uym0ee7Mndps4QTt+vFCq5fpWLGNuyfcJe28CxlVZE/g3F7sI0B7rXWxUqoP8LlSqrvW+uTPT1RKTQImAURGRpKamlrnm5eXlxMcHFyv4O6gqeR39YYIpw7voWzPNu5o1YO2lUGkpmb/6jnl5eX1+plzFxfLf43jR/yVDR8c2KwWVn3+PvQ45VZv4Xvz69+U+fr6ntn++0JSU1Pp1atXIyZyLskv3InWms+LoxnVsy+062l0HK/W6EW2UsoM3AT0OX1Ma10NVNd+vlkplQV0ATb9/Hyt9VxgLkBKSooeOnRonTOkpqZSn/PcheRvuOoDa9GpM7jGbMNU/QWq06X12rpD9oa4aP68IJi3GG23YMfMrMNdiGoVyp/HdHebYTRe/foLIUQjOFRSSWmlVTahaQRGzPAaAWRqrfNPH1BKRSilTLWfxwGdgQMGZBNNxKoVn2HWNsw4UPbaiR9N3emdIoc/jfneLxk+8nq+3V3AiL+u4r9b8mW8rBBCeAHZ6bHxuHIJv4+BoUC4UiofmKG1fhe4nfMnPA4B/p9Sygo4gAe11hecNClEQ/24v4i3cqIYFuALeMGGM85Uu3mNCZjaAa7p3prHP93Bw59s54tth3mlXxWtT2z0qvXDhRCiKdl9+CQ+CuJbS5Htai4rsrXWE37h+D0XOPYp8KmrsghxWqXFzlOf7US1TMZx8xeQv0YKxovoHBnCogcvZ/7aHL755ktCc5/HoWwosx9KJkYKIYTH2XWolI4RwQT6mYyO4vVkx0fRpPz9+33kFlfw0e/64x8bDrEDjY7k9kw+insGxTK+shK/1bUTI20WNv7wOd1u6UVokK/REYUQQlwCrTU78ku4sstlRkdpErx31w0hfmbXoVLeSc/mtpR2XN5RNl2pq5CuwzCZ/dHKhF358lpmBJe//D3PLdlN1pYf0Gmvy/raQgjhxg6XVlFUbiGpXajRUZoE6ckWTYLDoXnqs52EBfnx1KgEo+N4ptMTI3PS8Y8ZzAumeP6VlsXuDStos/lFHMqGXfmyJ+lJ2vhV0qr7VZg6eNnW9EII4cF25JUAkBjdwuAkTYMU2aJJ+OG7L7niyFf0Hz5Whjc0RO3ESIBuwD9u70XlDyvwT6sZRqIdFrpvfR6FxmfLGzBxiYzbFkIIN7E9vxRfkyIhKsToKE2CDBcRXu9U1hquWPNbHvFdxJC1v5UhDU4W2HkoPmZ/UCZMJhNmpTErB8pukaURhRDCjezILyG+dXP8zTLpsTFIkS28XsbyufhpCyak8HOJ2mEkDH8aNep1VG3BLUsjCiGE+3A4NDvzS0mMlvHYjUWGiwivlr8jlZ7HvuTMzuA+Zin8XOGsYSREdqv5RUaWRhRCCLeRXXyKsmobSTIeu9FIkS282rofvmCcclBTYyvodYcUfq52dsEthBDCLezIr530KCuLNBoZLiK81voDxXx0tD3ax7dm+II5AJLuMDqWEEII0ei255US6GuiU0Sw0VGaDOnJFl5Ja80ryzM5FNITxwTZ2VEIIUTTtiO/hB5tm2M2Sf9qY5EiW3ilDenLGXDoM5KHXI9/7AjZ2VEIIUSTZbU72H34JHcN6GB0lCZFimzhdey560n64W5SfG34bFwC3VpLD7YQQogma9/RMqptDllZpJHJewbC6+xZuxSzttVryb7CWbNdmEwIIYRofDvySwHZ6bGxSZEtvIrV7mBOdmvsyoyux1rNRXPm/OpzpBAXQgjhSbbkniAsyJeYVkFGR2lSZLiI8CpfbDvM8tIO7Bg1n75qd50mO56oOgHArK2zOFJ+hBPVNY8DzYFENYuiY4uO2Kw2IufMIWLaVJf9HYTwFEqprsDCsw7FAX/SWv/doEhCiAvYfPAEvduHoc5sGiEagxTZwmvYHZo3V+6nW1RzUgZfAeraXz3H6rCy4YWHafnRd2eOjZjwJgCpIy9j1TVRHLIdIi0/jWp7NQCfALO3zmZ85/FEBUdROGu2FN2iSdJa7wWSAZRSJuAQ8JmhoYQQ5yipsHCg8BTje0cbHaXJkSJbeI2vdhzmQNEp/nln71/9bd3usPPZ/s/4145/UdChgPYvxXFNzDWMvOOfxOzaQqA5kARgcu3zj77xBsff/OeZ86+aMIcS5pA+JoGeSzKImDZVim3R1F0FZGmtc40OIoT4n60Hazah6d0+zOAkTY8U2cIrOByaOSv30yUymGu6t77oczOPZ/LnNX9mV/EukiKSeHbAswxuOxilFBn8k0Bz4HnnRE6fTuT06aSmphL54GSab/qOBXsWsHDvQj5cAu/vep8BMoxENG23Ax9f6AtKqUnAJIDIyEhSU1PrdOHy8vI6n+NOJL+xPDm/M7L/9ycLPgpOZu8gNa9xh4t48msPDc8vRbbwCt/sLmDf0XLemNALH58LNyJaaz7K/IjXN71OqH8oLw9+mVGxo87p9Q6fMuWS7uf3/meMmfNvxtQ+HnDzqwBk//VlYh9+okF/FyE8jVLKDxgDPHmhr2ut5wJzAVJSUvTQoUPrdP3U1FTqeo47kfzG8uT8zsg+96d1JERZuWbEpS8C4Cye/NpDw/PL6iLC42mtmb1yP3HhzRjdM+q8rxfOmo3FbuHJ1U/y8oaXGdhmIP8d819Gx40+b1jJpfREh0+ZQsS0qSRkZpxXlFfNnUdGfIKsQCKamuuALVrro0YHEUL8j92h2Z5XIkNFDCJFtvB4a7OK2X34JJOGxGG6QC920Zw5TP5uMksPLGVar2nMHj6bsID6NzhnF+KnP0/IzADg4RdjuOvpIPaM61nv6wvhgSbwC0NFhBDG2VtQximLnT4dpMg2ghTZwuO9nX6A8GA/xvVqe97Xyi3lAGw5uoWXrniJSYmTnL6E0dm92VaY+8sAACAASURBVB+P/pjOYZ35/crfsyJ3hVPvI4Q7Uko1A0YC/zU6ixDiXAe3p/KQ6QsG+GYZHaVJkiJbeLSfjpaxcm8hdw+MIcDXdOZ44azZZMQnkJfYF4CPXqqi0+g/umQYx+ne7PApU2gR0IK3r36bHuE9eGzVY6Tlpzn9fkK4E631Ka11K611qdFZhBBnydvA8A3387DvIiI/vxXyNhidqMmRIlt4tHfSswnw9eGuAR3OOR425UHemD2U25/yA2qGcyRkZrh09Y/T1w7xC+GtkW/RJawLj656lB2FO2SMthBCiMaVk45JWzHjQNktkJNudKImR4ps4bGO703nsu1z+EN8CS2b+Z3ztdc3vc7qQ6t5ZsAzhmRr5tuMN0e8SauAVkz5fsolbdcuhBBCOEtpZH8s2owDE5j8anZAFo1KimzhmfI2ELJwPP/n8wm/y/79OW+DLd63mAUZC7gr4S5u6XLLJS/L52zhgeG8NfIt7NoOQKWt0pAcQgghmp4Nts7caXmKQ70fholLoF0/oyM1OVJkC49kzUpDOayYlQMfu/XM22AZxRm8tP4lBrUZxCMpjwCXtiyfKxTOmk1Fv2t5988nAMjp0VuW9xNCCNEoNuYcZ5dPPBHXPSkFtkGkyBYeabU1Hqs241D/exvslPUUj6U9RlhAGH8Z/BfMPsbutXR6Le3Ty/vd+qSZ3Uuel10hhRBCuNz67OMktQs9Z1EA0bikyBYe6R/7wvhjs+dRw54+8zbYS+tfIq8sj5cHv9ygdbBdZUDUAF7d+Cq5J3ONjiKEEMKLnaq2setQKf1iWxodpUmTIlt4nJ35pWzLK6HPoGtQQx6Bdv1YkbuCJVlLeCDxAfq27mt0xPOET5nC84Oex+xj5qn0p7A5bEZHEkII4aW2HizB7tD0i21ldJQmTYps4XHmr8shyM/ETX2iASitLuXFdS+S0DKBSYmTDE53YRHTptK6WWv+NPBP7Cjawds73zY6khBCCC+1IbsYH4Xs9GgwKbKFRympsPDFtsOM69WW5gG+ALy68VVKq0vP9BS7s2tjruW62OuYu2MuB0oPGB1HCCGEF1qffZwebUMJ9nfvfxO9nRTZwqMs3pxPtc3Bb2o3n1lzaA1LspZwb4976dqyq8HpLs3jfR8nyBzE/1v7/3Boh9FxhBBCeJFqm52teSX0i5Hx2EaTIlt4DIdDM39dLn1jwkiIao7VbuUvG/5CTPMYHkx60Oh4l6xVYCseSXmEzUc388X+L4yOI4QQwovsyC/FYnPIpEc3IEW28BhrsorJLa44s4V66nOTyTmZwx/7/hE/k9+vnO1exnUaR+/LejNz00xKqkqMjiOEEMJLbMg+DkBf6ck2nBTZwmMs3JRHaKAv13RvTVFlEdGLfmRI9BAGR3veVrE+yodnBjxDubWcf27/p9FxhBBCeIn12cfpGhlCWDPP6nzyRlJkC49QUmHhm90F3NirLQG+Jt7Y8gYAj6U8ZnCy+usc1plbutzCwr0LOVAikyCFEEI0jM3uYHPOcfrGyqoi7kCKbOERvth2GIvNwa27lpERn8CE3y4CoLL/dR69VflDyQ8RZA5i5qaZRkcRQgjh4fYcOckpi13Wx3YTsraL8AgLN+bRo21zuk17lKn9c9hybAvv/vnEmS3LPVXLgJZMSpzE65tfZ82hNVze9nKjIwkhhPBQ6w/UjMeWlUXcg/RkC7e361Ape46c5LaUdmw7to1V+au4r8d9RsdymjsS7iA6OJq/bv4rx2bNMjqOEEIID7Umq4i48Ga0Dg0wOopAimzhARZuzMPf7MMNSW34x5Z/0DKgJXfE30H4lClGR3MKP5MfDyU/xN4Teyme86bRcYQQQnggq93BhuzjXN5Jhoq4CymyhVurstr5fNshru3Rmj0lm9h0dBOTEicR5BtExLSpRsdzmlGxo4gLjQPA7rAbnEYIIYSn2ZFfwimLnUEdw42OImpJkS3c2je7CyirsnFrn2je3PYmUc2iuKXLLUbHcqrCWbPZ160HLz+xD4B93Xp49GROIYQQje/H/cUoBQPipCfbXUiRLdzaZ1sP0bZFIOZm2Wwv3M59Pe7zuI1nfk3EtKkkZGbQNWM3AI+8FEunPTu8qqdeCCGEa63JKqJbVHNZH9uNSJEt3NaJvavpnvUOkzsW8+6ud2gV0IpxncYZHctlfFTN/455ZXksPbDU4DRCCCE8RaXFzpbcEgZ1kqEi7sRlRbZS6j2l1DGl1K6zjj2nlDqklNpW+zHqrK89qZTar5Taq5S6xlW5hIfI20DIwpv4g2kRidm/Z+2Rtdzd/W4CzN49Y7rVlIfoEtaF93a9h0M7jI4jhBDCA2zOPYHF7mBgRxkq4k5c2ZP9b+DaCxz/m9Y6ufbjawClVDfgdqB77TlvKqVMLswm3F1OOsphxawcvBccSIiPH7d2udXoVC532bRp/LbHb8kuzeaHgz8YHUcIIYQH+DGrCLOPkvWx3YzLimytdRpw/BKfPhb4j9a6WmudDewH+rkqm3B/h1r0waLN7N3Tgu+bBXJH+2sI9gs2OlajuDrmaqKDo3l357torY2OI4QQws2t2V9Er/YtaOYvewy6EyPGZE9VSu2oHU4SVnusLZB31nPya4+JJmphQRvutD6FY0cQfj5m7uj/qNGRGo3Zx8y9Pe5lV/Eu1hesNzqOEEIIN1ZaaWXnoVIGytJ9bqexf+X5J/A8oGv/fB2o09Z9SqlJwCSAyMhIUlNT6xyivLy8Xue5C2/Pr7Xm4zWVtGrZHoDegSnsWLejkdJdXGO99mE6jOam5ryW9hrTIqc57bre/rPj7jw9vxDC/aw/UIxDwyAZj+12GrXI1lofPf25Uupt4Kvah4eAdmc9Nbr22IWuMReYC5CSkqKHDh1a5xypqanU5zx34e35d7zwGh98/N6Zxw8/sxpYTfiUKYYva9eYr33erjz+tvlvtE5sTXzLeKdc09t/dtydp+cXQrifNVnFBPqa6NU+7NefLBpVow4XUUpFnfXwRuD0yiNLgNuVUv5KqVigM7ChMbMJ97G453WMu/k1vh4WAkBCZgYJmRmGF9iNbXzn8QSaA1mwZ4HRUYQQQriptH2FDIhriZ9ZVmV2N65cwu9jYC3QVSmVr5T6LfCqUmqnUmoHMAz4A4DWejfwCbAHWA5M0VrL3tJNkM3uYOmOIyTFH2LUyjKj4xgq1D+UMR3H8HX21xRXFhsdRwghhJs5ujuNa098xPjLDhsdRVyAy4aLaK0nXODwuxd5/ovAi67KIzzD+uzjFJ+y0DYwFYCWD002NpDB7ky4k4V7F7Jo3yIeTHrQ6DhCCCHcRd4GWn16Cw+bLfhs/QIS20A7WZjNnch7C8Kt5P31Hyz7/FFeeXobAMff/CcZ8QkUzpptcDJjxIbGMqjtIBbuXYjVbjU6jhBCCHeRk45yWDArB8puhZx0oxOJn5EFFYXbsNkdvBo1hMuePsUxUvnwxUoSMjOMjmW43yT8hge/e5DlOcu5oeMNRscRQgjhBmztB2HVZvyUDZPJD2IGGx1J/Iz0ZAu3se7AcY5XnOK4WsuI9iOMjuM2Lm9zObGhsXyY8aFsTiOEEAKAzfbO3Gl5iqwev4eJS2SoiBuSIlu4jaU7j9AsbDeV9nLGdxlP+JQpRkdyC0op7oi/g93Fu9ldvNvoOEIIIdzAqn2F7FBdibr+KSmw3ZQU2cIt2OwOvtldQFjrrbQLaUff1n2b3JJ9F3N93PUEmgNZtG+R0VGEEEK4gbSfCundIYyQAF+jo4hfIEW2cAvrDhznhPUQpXovN3W+CR8lP5pnC/YLZlTsKJZlL6PM0rSXNhRCiKausKyaXYdOcmWXCKOjiIuQSka4haU7D9Os5WZMysTYjmONjuOWbu5yM5W2SpYeWGp0FCGEEAZavb8QQIpsNydFtjCc1e5g2a7D+Idt54q2VxARJI3GhXRv1Z2Elgks2rdIJkAKIUQTlraviFbN/OgW1dzoKOIipMgWhtuQfZwbdizAwgmu73i90XHcllKKm7vczL4T+9hZtNPoOEIIIQxgd2hW7StkSJcIfHyU0XHERUiRLQz37e4C7ti1mWa+wQyNHmp0HLc2Om40QeYgmQAphBBN1La8Exw/ZWF4/GVGRxG/QopsYSitNd9k5AFwdYeRBJgDDE7k3pr5NmNU3CiWZy8n/++vGx1HCJRSLZRSi5VSmUqpDKXUQKMzCeHNvs84htlHMUTGY7s9KbKFYQpnzSYzoRvv//thACb8dlGT3kL9Ut3U6Saq7FWUvfWO0VGEAPgHsFxrHQ8kAbJNqxAu9H3GMfrGtCQ0UJbuc3eyrbowTMS0qaxoE8kHR9/gb28U0TVjtyzddwl6hPcgLjQO2Gd0FNHEKaVCgSHAPQBaawtgMTKTEN4s73gFe4+W8czoBKOjiEsgFY0wTt4Ghu2eQkHQCQB88jcZHMj9ne79f/mJmgI7Iz5Bev+FkWKBQuB9pdRWpdQ7SqlmRocSwlv9kHkMgKsSIg1OIi6F9GQLw5zY/QPfN/PDrhQqsQJy0mVr2F8RMW0qEdOmcqziGMW9r2TFR5OZ3nu60bFE02UGegPTtNbrlVL/AJ4Anj37SUqpScAkgMjISFJTU+t0k/Ly8jqf404kv7E8Of/Psy/aVEXrIEXuro3kGhfrknnyaw8Nzy9FtjBMmrUrS4Ob0dViIT6xGmIGGx3JY1wWdBnFwJKsJUxJnoLJx2R0JNE05QP5Wuv1tY8XU1Nkn0NrPReYC5CSkqKHDh1ap5ukpqZS13PcieQ3lifnPzv7qWob+1as4O6BMQwd2s3YYJfIk197aHh+GS4iDPPeYRO7Avy4rvVAmLhEerHrqOTOazhacZQNBRuMjiKaKK11AZCnlOpae+gqYI+BkYTwWuk/FWGxO2SoiAeRIlsYoqi8moyTqwG4ZugLUmDXQ++nXiHEN4QlWUuMjiKatmnAh0qpHUAy8JLBeYTwSj9kHiUkwExKTJjRUcQlkiJbGOL7jKOYQ3YS1zye6JBoo+N4JH+TP9fFXsd3ud9Rbik3Oo5oorTW27TWKVrrRK31OK31CaMzCeFtHA7ND5mFDO16Gb4mKd08hXynhCG+3L0LU2A+YztdZ3QUjzam0xiq7FWsyF1hdBQhhBAusj2/hKLyaobHywY0nkSKbNHoLHbNluJVAFwTe43BaTxbYngi0cHRLMteZnQUIYQQLvLN7qOYfRTDu8p4bE8iRbZodBnH7dBsOx2C42kb3NboOB5NKcV1sdexvmA9RZVFRscRQgjhZFprlu86wsCOrQgNkl0ePYkU2aLRrS86iinwEOM6jzI6ilcYHTcah3bwTc43RkcRQgjhZPuOlpNTXME13VsbHUXUkRTZolFprdlTvQ2AUXEyVMQZOrboSJewLjJkRAghvND2td/ykPkLRrU4aHQUUUdSZItGtf9YOZaAnbQJ6Eqb4DZGx/Ea18Vex/bC7eSX5RsdRQghhJM0L81kzPYHecS8iJaLb4E82RfBk0iRLRrVf3duxxR4iOs7Xmt0FK9yXWzNKi3Lc5YbnEQIIYSz+BbuxKxtmHCA3QI56UZHEnUgRbZoVN/mfAfA+HgZj+1MbYPbkhyRzNIDS42OIoQQwknSbQlYMaOVCUx+EDPY6EiiDqTIFo2mtNLKEdsmgh3RMlTEBUbFjWJ/yX72ndhndBQhhBBO8FlJR55u/iJq+NMwcYnsjuxhpMgWjWb1j59iCsilvynK6Che6eoOV2NSJpkAKYQQXuBYWRX7SxzEJA+DwY9Ige2BpMgWjSNvA2XbHwcFDx5aJpM3XKBVYCv6R/Xnm5xv0FobHUcIIUQDrNhzFA2ydJ8HkyJbNApHdjqpQb5EW610rq6SyRsuMrLDSPLK8th7Yq/RUYQQQjTA8l0FRAYpukQGGx1F1JMU2aJR7PCPZ0NgAMMqqsDHLJM3XGR4++H4KB++zfnW6ChCCCHqqbi8mjVZxfRtbUYpZXQcUU9SZItG8d/sTViVYnDUELYnPS9jy1ykZUBL+kb2ZUXuChkyIoQQHmr57gLsDk2/1iajo4gGkCJbuF7eBiqK5xNmt9Nvr/SwutrIDiPJOZlDVkmW0VGEEELUw1fbjxAX0Yx2IVKmeTL57gmXK9v3Az8G+TO0ohKT3UKLkl1GR/JqV3W4CoViRe4Ko6MIIYSoo2NlVazPLub6xDYyVMTDSZEtXO4TayDX/agZVlENJj9KWvQwOpJXCw8Mp3dkb77NlXcNhBDC0yzbWYBDww2Jstytp5MiW7jc58dzuGW1pl/f6TBxCSdD442O5PVGdhhZszHNzOeNjiKEEKIOvtx+mPjWIXSODDE6imggKbKFS2mtya3aBECzYU/IhMdGMqL9CADs73xkcBIhhBCX6nBJJZtyT3C99GJ7BSmyhcsUzppNZkI3Fr5aBEBGfAIZ8Qk0+/Irg5N5v8hmkSRHJBsdQwghRB1sWr2ch0xfMD7isNFRhBOYjQ4gvFfEtKn8Oa6KVUXz+OQvNhIyMwBITU01NpiXK5w1m6I5c3iq9nFGfAIA4VOmQE8ZDy+EEG4pbwPXbJrEaF8bpi+WQIslRicSDSQ92cKlthatxcfa1ugYTUrEtKkkZGYQuul7ANYt/iMJmRlETJtqcDIhhBC/pGTPD5i0DRMOsFtkZ2QvIEW2cJnjFSWc1PuJC0qp6UUVjapNcBsAVuatNDiJEEKIX/N9VResmNHKBCY/2RnZC8hwEeEy/9n9HSgHI2KuJOKOkUbHaZJyxvdj27GtFFUWER4YbnQcIYQQF6C1ZtXeQpoFjeDa7q0haULNQgFZqUZHEw0gPdnCZb7LXoW2B3Fb4hVGR2myYh95Eo0mLT/N6ChCCCF+wYGtK3nl1DNcXbkctn1sdBzhJFJkC5dwaAe5p9bTzRJC+ImdRsdpsrqGdaVNszasPChDRoQQwl3lbv4GX2z4yHhsr+KyIlsp9Z5S6phSatdZx15TSmUqpXYopT5TSrWoPR6jlKpUSm2r/XjLVblE49i8/T9YfCr5TdkemDcG8jYYHalJUkoxtN1Q1h5ZS4W1wug4QgghfsZmdzC/oD0OH1+Q8dhexZU92f8Grv3ZsRVAD611IrAPePKsr2VprZNrPx50YS7RCL7J+AKlNYMrK+W3coMNaz+Mans1a4+sNTqKEEKIn1m9v4iVp2LYMnQeDH8aJi6Rjdu8hMuKbK11GnD8Z8e+1Vrbah+uA6JddX9hrFXWE/SothKqlfxWbrA+kX0I8QuRISNCCOGGPtt6iNBAX3oNuhoGPyIFthcxcnWR+4CFZz2OVUptBU4Cz2itL9j1qZSaBEwCiIyMrNfGJuXl5R69IYq75y+zl1GgjxBkGUh2TDAlLXpwMqvizCxpd89/MZ6avatvV77P/p5+Lfp5ZP7TPPX1P83T8wshnKu82sY3uwsY3zsaf7PJ6DjCyQwpspVSTwM24MPaQ0eA9lrrYqVUH+BzpVR3rfXJn5+rtZ4LzAVISUnRQ4cOrfP9U1NTqc957sLd83+69yvIh25xNxM3+vrzvu7u+S/GU7NX51Tz6KpH8f/uc4b+4z2j49Sbp77+p3l6fiGEcy3fVUCV1cFNvWXTNm/U6KuLKKXuAa4H7tRaawCtdbXWurj2881AFtClsbMJ51i2Pw1tD+CGeHnLy11c0fYKfH18af+NjMsWQgh3sWvdCp4M+Zre6iejowgXaNQiWyl1LfBHYIzWuuKs4xFKKVPt53FAZ+BAY2YTzqG1ZufxjejKTvTp0MroOKJWM99m9Iuq+aWn9ndbIYQQBirOSOfxY3/kd9aPUB+MlVW4vJArl/D7GFgLdFVK5SulfgvMBkKAFT9bqm8IsEMptQ1YDDyotT5+wQsLt5Z7MpcKRxHtApNkfJmbKJw1m4z4BP5v6ioAMhO6kRGfQOGs2QYnE0KIpmvfhmWyNraXc9mYbK31hAscfvcXnvsp8KmrsojG833OagAGt73c4CTitIhpU4mYNpWjp45yvM9Q1i5+jPt63Gd0LCGEaLK01sw/0o7eyhezsskqXF5KdnwUTvVtdhoOSxjXdO1hdBTxM5HNIgFki3UhhDDY+uzjfF3SnrVXvCdrY3sxKbKF09gcNvaVboOqLiS1a2F0HHEBW67qxLZj28j7+0yjowghRJP1yaY8QvzN9B9ynayN7cWkyBZOs6toFzYq6disF74m+dFyR9U33Ihd2yl/64Ijt4QQQrjYySorX+88wg3JbQj0k7lL3kwqIeE032Wno7VieIyMx3ZXMf4xtPCXdxmEEMIoX24/TJXVwW0p7YyOIlzMyB0fhZdZeXA1jqq2DO8Sa3QUcQGFs2YTNWdOzU5OQEZ8AgDhU6YQMW2qccGEEKIJ+WRTPl0jQ0iMDjU6inAxKbKFU1RYK8ir2Iup+kq6t5GGwx1FTJvK7p49qGhfQeyoR7CmLyQxItHoWEII0WTsLShje14Jz17fDaWU0XGEi8lwEeEUW49tRWOna4tkTD7ScLizQW0HAbLKiBBCNLaFG/PwNSlu7CXbqDcFUmQLp1iZuw6tfRgeIzOk3V2ofyirr2kjRbZwCqVUjlJqZ+0GY5uMziOEu6q22flsaz4ju0XSspmf0XFEI5AiWzjFmkPrcVRGM6RTtNFRxCVQ908g43gGhRWFRkcR3mGY1jpZa51idBAh3NU3u49yosLKrTLhscmQIls0WIW1gvyKfZiqO5EQ1dzoOOISDIkeAkD6IdnGVwghGsOCdbm0bxnEkM4RRkcRjUQmPooG21a4DY2dzqFJMh7bQ3Ru0ZnIoEjS8tO4qfNNRscRnk0D3yqlNPAvrfXcnz9BKTUJmAQQGRlJampqnW5QXl5e53PcieQ3ljvkP1TmYEN2Jbd28SUtbdUln+cO2RuiqeeXIls0WNrBmvHYQzvIeGxPoZRiSPQQlh5YisVuwc8k4wNFvV2htT6klLoMWKGUytRanzPgv7bwnguQkpKihw4dWqcbpKamUtdz3InkN5Y75H9uyW78TAd5/NYraRXsf8nnuUP2hmjq+WW4iGiw1fk147EHdZTZ0p5kSPQQKmwVbD662egowoNprQ/V/nkM+AyQ37aFOEvVgbW02DyLyZ2K61RgC88nRbZokAprBQdP7YXqOHq2lfWxPUn/qP74+fjJKiOi3pRSzZRSIac/B64GdhmbSgg3krcB84JxTGUh0/MfgbwNRicSjUiKbNEg2wu3o7HTMTgJP7P8OHmSQHMgKa1TWHN4jdFRhOeKBFYrpbYDG4ClWuvlBmcSwn3kpKMcFszKgY/DCjky2bwpkapINMiazM9QGm4JCzQ6iqiHy9tczoHSAxScKjA6ivBAWusDWuuk2o/uWusXjc4khDv5KTAZizbjwIQy+UHMYKMjiUYkRbaov7wNbP/pC3pUV3PnnoflbTAPNKhNze6PPx760eAkQgjhfeZmh3OffhbrlU/CxCXQTqYsNCVSZIt6qzqwkp3+vqRUVcvbYB6qY4uORAZF8uNhKbKFEMKZSiusfLnjMDHJw/Af9pgU2E2QFNmi3naFtcGmFElVFnkbzEMppRjUdhDrDq/D5rAZHUcIIbzGfzYepMrq4DcDOhgdRRhEimxRbxvtlQBUtZogb4N5sMvbXE6ZtYydRTuNjiKEEF7BZnfwwdpc+se2pFsb2Qm5qZIiW9Rbet5G7NURBF3xhBTYHmxA1AB8lI+MyxZCCCf5LuMoh0oquXdQrNFRhIGkyBb14tAO9pbswlHZgT4xYUbHEQ0Q6h9Kz/CespSfEEI4yXs/5hAdFsjIbpFGRxEGkiJb1Et2aTYWXU6kXwLNA3yNjiMaaFCbQewq2sWJqhNGRxFCCI+2+3ApG7KPM3FgDCYfZXQcYSApskW9bCzYAkDfqF4GJxHOMKjtIDSadUfWGR1FCCE82vs/5hDkZ+LWvu2MjiIMJkW2qJdVuRu4eZWJ4XHdjY4inKB7q+6E+ofKuGwhhGiAovJqlmw7zPje0YQGyru8TZ0U2aJedhZt59Y11aTEtDQ6inACk4+JAVEDWHN4DcdmzTI6jhBCeKSP1h/EYndwz6AYo6MINyBFtqizosoiSm1HALiseYDBaYSzDGoziMLKQornvGl0FCGE8DgWm4P563K5sksEHSOCjY4j3IAU2aJOCmfNprDXYD75S83GJRnxCWTEJ1A4a7bByURDXd7mcqMjCCGEx/pqx2EKy6q5V3qxRS2z0QGEZ4mYNpXZyWV8uv8TFr1SRUJmhtGRhBMUzprN8Tlz+KT2cUZ8AgDhU6YQMW2qccGEEMIDaK2Zm3aArpEhXNklwug4wk1IT7aosw0FW7BXtTU6hnCiiGlTScjM4OsF9wPQYedmEjIzpMAWQohLkPZTEZkFZfxuSBxKybJ9ooYU2aJOqmxVHKr4CZMllpYPPWR0HOFkp4eMbDq6yeAkQgjhIfI2cOSrFxkRnMuYpDZGpxFuRIaLiDrZU7wHjZ24kB5E3veA0XGEk/WO7M3Mwb6YD69lSPQQo+MIIYR7y9uA4983cLPNws0mX8xHUqBdP6NTCTchPdmiTjYe2QbAgDayCY038jf5k3NLf9mUpolRSr2qlGqulPJVSn2vlCpUSt1ldC4h3F5OOtpuwawcmLQNctKNTiTciBTZok7W5G/BYWnB4I6xRkcRLjKwzUD2l+ynsKLQ6Cii8VyttT4JXA/kAJ2AxwxNJIQHKAjri0WbsWNCmfwgZrDRkYQbkSJb1Mnekt3Yq9qT3K6F0VGEiwyIGgAgvdlNy+mhg6OBRVrrUiPDCOEp3sxqyd32p6m44nGYuESGiohzSJEtLllRZRGn7EW0MnciJEC2i/VW8S3jaeHfQorspuUrpVQm0Af4XikVAVQZnEkIt3b8lIVPNuURmzyMkBGPS4EtznPJRXbteL2Wpz9cGUq4p+3HdgCQGNHT4CTClXyUD/2j+rPu8Dq01kbHES6klGoDoLV+ONjmiAAAIABJREFUArgcSNFaW4EKYKyR2YRwd/PW5FBldTBpSJzRUYSb+tUiWyn1gFKqANgBbK79kPW9mqBVuZvR2odhsTLp0dsNiBrAscpjZJdmGx1FuNY7Sql1SqmXgURAAWitT2mtC4yNJoT7Kquy8v6P2VzdLZJOl4UYHUe4qUtZwu9RoIfWusjVYYR723J0O47q1gyMjTI6inCx0+Oy1x5ZS1wL6aXxVlrrUUqpAGAocCMwUyl1EFj+/9m77/Aoq7SP498zkwYJnZCEGjoBaRISqlJERUVdXRArViwQu2vdV13XvrZNYgFFVFDELigKKKGToSOQ0AMBkhB6Gkkmc94/ZkBkKWmTM+X+XNdcyfNknsyPMU7unLnPOcAvWutdJvMJ4ak+XbqTo8fsJAxpbzqK8GDlaRfZhvOtQ+HHHNrB7sLNBNmjad6gluk4ws2a12lOizotWLZX+rJ9ndb6mNb6F631A1rrWOARnAMwSUopm+F4QnicwhI7Hy3awaCO4XRtXs90HOHByjOS/SSwRCmVChQfP6m1vt9tqYTHyTiSQRlFtKsTI1vG+ok+UX34ecfPlDpKCbTIRFd/oJSqCxwBprlu+WYTCeF5Pk/dxcGCEhKGtDMdRXi48oxkfwD8Dizjz57sle4MJTzPkt2rAIiL6mE4iagpfZv2paC0gPX715uOItzsTHNvtNYlZpMJ4VmOlZYxYcF2+rVtRK9WsgaEOLvyjGQHaq0fdnsS4dEWZa5ElwUzqE0X01FEDYmLjEOhWLZ3GT2byGRXHydzb4Qoh69WZLIvr5i3R8uAkzi38oxkz1JKjVVKRckSfv4r7eAGHMda0K15A9NRRA2pF1yPzo06y3rZ/kHm3ghxDiV2B+/P305sqwb0bdPIdBzhBcozkn296+OTJ53TwDmXHFBKTcK5Te8+rfV5rnMNgS+BaJzb947SWh9Szkbfd4DLcL7Y36q1XlW+f4ZwpyJ7EQdLd1LPMozQ4PL8yAhf0bdpXyavn0xBaQGhgaGm4wj3kbk3QpzDd6t3s+dwES/+7TyZmyTKpTwj2W211q1PvgGdy/n9JwOXnnLuCeA3rXV74DfXMcBwoL3rNhZ4r5yPIdws7UAaKAedGkqriL/pE9UHu7azIluWxvdxMvdGiLMoLXOQPG8b3ZrX48IO4abjCC9RniL7w5MPlFKhwE/l+eZa6wXAwVNOXwV84vr8E+Dqk85/qp2WAfWVUrIgsweYn+H8Xdu/ufTl+pseTXoQbA1madZS01GEewVqrR/WWn+stf7k+M10KCE8RcrcmVxx5Av+2T1fRrFFuZWnyN6jlHoXQCnVAJgDTKnCY0ZorbNcn2cDEa7PmwGZJ91vt+ucMMyWtQZHaX0GtJFNSfxNsDWYXhG9ZL1s3ydzb4Q4g5KMpQxcegePBH5F7PxbIVOWjxflc84GW631P5VSryml3gd6Aa9orb+pjgfXWmullK7INUqpsTjbSYiIiCAlJaXCj5ufn1+p6zxFTeffcmgDHGvO3rSVZKdX/S94b37+vTk7VC5/eGE4S44s4Ye5P1AvwOzGC/74/NeQSs+9EcLXrV/8M920HatyQFkJZCyEFnGmYwkvcMYiWyl1zUmHqcA/ARuglVLXaK2/reRj5iilorTWWa52kH2u83uAFifdr7nr3F9orScAEwBiY2P1oEGDKhwgJSWFylznKWoy/5HiIxTvPEBkrQsYMnhwtXxPb37+vTk7VC5/xIEIfpj5A6q1YlDbil1b3fzx+a8hbbXWjpNPuLZbF8KvFZWUkbQ9kvdUIFZlR1mDIHqg6VjCS5ytXWTESbcrgNVA4EnHlfUjMMb1+Rjgh5PO36Kc+gBHTmorEYas2+B80+KC2rKVur/q2LAjDYIbSMuIb6v03BshfNmUZTv5vSCarcOnooY8DWN+lFFsUW5nHMnWWt9W1W+ulPoCGAQ0VkrtBp4FXgGmK6XuAHYCo1x3/xnn8n1bcS7hV+XHF1WUaWPzghegQRjjt78DmUPlxcUPWZSF+Kh4lmUtQ2stk3580x6l1Lta6/tcc29+AiaaDiWESQXFdt6bv42B7RvTJT4eGGY6kvAybl30WGt9/Rm+NPQ099XAOHfmERWUsZCNgRaaldppaJc+NH8WHxXPLxm/sOPoDtrUkzZdX+POuTdCeKvJSzI4WFDCw8M6mI4ivFR5VhcR/ip6IBuDg+hUUiJ9aH4uPjIeAFuWzKr3JUqpa47fcM696YOzNVCfMi9HCL9ypKiUD+ZvY2inJvRsKTsdi8qRIluc0ZEmHdkdGICydpM+ND/XvE5zokKjsGVLke1j3DX3RgivNmHBNo4es/OQjGKLKqhwu4hS6iogW2ud6oY8woMs3/MHALVa3iQFtp9TShEXGUfK7hQc2oFFyd/nvqA65t4I4WsOpC3EungKCR0Gcl4zs8uWCu9Wmd+U8cAzSqlZ1R1GeJZ5O1cBMLi17PQonH3ZR4qPsPnQZtNRhBDCPTJt1Jl+Lfer6Ty09xHZeEZUSYWLbK31U1rrEVrr4e4IJDzHun3rcZQ0oG/rlqajCA/QO7I3AKlZ8iaW+CullFUptVopNdN0FiGq4vDSTwlwFBOgHFjKSp0T/oWopHMW2Uqp2kqpfyqlJrqO2yulpFfPD2QVbaU20YQFu3URGuElIkMjia4bLX3Z4nQeANJMhxCiSjJthG6chsK53SmWAJnwL6qkPCPZHwPFQF/X8R7g325LJDzCkWNHKFa5tKjd3nQU4UHiIuNYkb2CUkep6SjCjZRSVyml4st53+bA5ZyyoY0Q3mb36tkoXYZSoFDQ8waZjySqpDxFdlut9WtAKYDWuhCQ3Sh83OLMNQB0a9LFcBLhSeKi4ii0F7LxwEbTUYR7VWTuzdvAPwDHue4ohKfSWpO0IxK7CkArKwSEQPcbTMcSXq48fQAlSqlauN49UUq1xTmyLXzYgp2rAZn0KP7qeF+2LctG9/DuhtMId9FaP1We+7laB/dprVcqpQad5X5jgbEAERERpKSkVChPfn5+ha/xJJLfrPLkX5FtZ1pWFO3bPsvgoHQO1z+Po9sKYdvZr3M3f3juPVlV85enyH4W+AVooZSaCvQHbq30IwqvsOHARhwlDYlrKZMexZ8ahjSkQ4MOpGancle3u0zHEdVEKbUSmAR8rrU+VIFL+wNXKqUuA0KAukqpKVrrm06+k9Z6AjABIDY2Vg8aNKhC+VJSUqjoNZ5E8pt1rvwldgfPv72Adk0CGXP7fQRYPWeJUl9/7j1dVfOf8ydJaz0HuAZnYf0FEKu1Tqn0Iwqv4Jz02IpaQVbTUYSHiYuMY82+NWT9923TUUT1uQ5oCixXSk1TSl2ilDpnW6DW+kmtdXOtdTQwGvj91AJbCE/3yy8/MPzQ57wad8yjCmzh/c7406SUOv/4DWgFZAF7gZauc8JHHT52WCY9ijOKj4qnuKyYw+9+YDqKqCZa661a66eBDsDnOEe1dyqlnldKNTSbTgj3ObJ5EcNWjOWRwK/olTJG1sUW1eps7SJvuD6GALHAWpwTHrsBK/hztRHhYxZnrgWgR0RXw0mEJ+oV0Ut2fPRBSqluwG3AZcA3wFRgAPA70ONc17ve4UxxX0Ihqt+y339kqLZjVQ4oK3Guiy0riohqcsbflFrrwVrrwThHsM/XWsdqrXsBPXEu4yd81ELXpMcLo2XSo/ir3MQkdneLY9pLJQCkdYohrVMMuYlJhpOJqnD1ZL8FLAe6aa3v11qnaq3fALabTSeEe6RlHWXCrii0JRCUFaxBsi62qFblmfjYUWv9x/EDrfV6pVSMGzMJw9Yf2OCa9NjcdBThYcITxhOeMJ63V77NJTd+QKs/VlI7sLbpWKLqRmqt/1JMK6Vaa613aK2vMRVKCHfRWvOvGRvZFtKFY9d/R2D2MmeBLaPYohqV5z3fdUqpD5VSg1y3icA6dwcT5mQXbac2LWXSozij+CjnPiWr9q0ynERUk6/LeU4In/DrhhyWbj/Aw8M6UKd9fxj4iBTYotqVZyT7NuBenNvmAiwA3nNbImFUXnEexWof7UIvNB1FeLAeTXrwxsAAamfZGNBsgOk4opKUUp2ALkA9pdTJI9Z1cc7HEcLnHCst46Wf0+gQEcYNcbJMrXCfcxbZWutjOHv13nJ/HGHakt3OzqCu4dIRJM6sVkAttl0bS2F2qukoomo6AlcA9YERJ53PA2QhdOGTPlq0g10HC5lyR7ws2Sfc6pxFtlJqB67dHk+mtW7jlkTCqEU7nZ1AF0afczEB4efiouJ4b817HCk+Qr3geqbjiErQWv8A/KCU6qu1Xmo6jxDutvtQIYm/b+HizhEMaN/YdBzh48rTLhJ70uchwEhA1k31URv2p6HtofRr1dp0FOHh4iPjeZd3WZGzgqEth5qOIypBKfUPrfVrwA1KqetP/brW+n4DsYRwm3/N2IhC8eyVXUxHEX6gPDs+Hjjptkdr/TZweQ1kEwbE/mojRLekVlB5/v4S/qxr467UCqiFLUs2b/Biaa6PK4CVp7kJ4TN+T89h9sYc7h/anmb1a5mOI/xAedpFTt7d0YJzZFsqMB9UYi/hb4sPs2zIINNRhBcItAZyfpPzsWVLke2ttNYzXB8/MZ1FCHcq3rGULd9MYUTDbtwxYLjpOMJPlKdYfuOkz+3ADmCUe+IIk1L3pNEY6NK4k+kowkvERcXx1sq32F+0n8a1pL/RWyml5uBcK/uw67gBME1rfYnZZEJUXd0j6Vg++z/uKCtFOb7GmtVTlusTNaI802rvOL77o9Z6mNZ6LFDi7mCi5uQmJpHWKYbGl4wG4JYHXpFd/ES5xEc618tenr3ccBJRReHHC2wArfUhoInBPEJUG+u+P1BlpQQoB1ZHqXPrdCFqQHmKbNmkwMeFJ4wnJj2NZ1+4CoCWa9cSk55GeMJ4w8mEp+vUsBN1AuuQmiVL+Xm5MqXUiQWDlVKtOM2qUkJ4G601nx3sSCkBaNk6XdSwM7aLyCYF/iezYCsAocFBhpMIb2G1WImNjJW+bO/3NLBIKTUfUMBAYKzZSEJU3cx1WXx9qC0XDpjIiHrbZOt0UaPO1pMtmxT4EYfDQb7eya8XtEW2oREVER8Vz7zMeezN30vTsKam44hK0Fr/4prk3sd16kGt9X6TmYSoqkMFJTz34wZa17Nw2WVXgUWZjiT8zBmLbNmkwL/8kbMDLMfYebWszigqJi7SOSpky7ZxdburDacRVdAPuOCk45mmgghRHf79UxpHikp5sG8IVimwhQFnaxeRTQr8yLwdawGIa9bVcBLhbdrVb0fDkIbYsqTI9lZKqVeA3sBU16kHlFL9tNZPGYwlRKUt3JLLN6t2M35wO1oEZ5mOI/zU2dpFTt6kQPi4NTkb0FpxUdtupqMIL6OUondkb2zZNrTWKCUjRl7oMqCH1toBoJT6BFgNSJEtvE5hiZ2nvvuDNo1DGT+kHcsWS5EtzDhbu8gM16eFWuuvTv6aUmqkW1OJGrcjbzPWsiaEh9U1HUV4objIOH7N+JVdebtoVbeV6TiicuoDB12f1zMZRIiqeGvOZjIPFvHl2D6EBFpNxxF+rDxL+D1ZznPCix2276RRYGvTMYSXio9yrpctS/l5rZeB1Uqpya5R7JXAi4YzCVFh63Yf5qNFO7ghviXxbRqZjiP83Nl6sofjfAuxmVLqvyd9qS7OnR+Fj9h9JBeH9RBt6nYwHUV4qZZ1WtKkdhOWZy9nVEfZENbbaK2/UEql4OzLBnhca51tMJIQFVZid/CPr9cRXieYJ4bLzsXCvLP1ZO/FOZpxpevjcXnAQ+4MJWrW3G1rADg/qovhJMJbKaWIj4xn8d7F0pftRVzL9p1st+tjU6VUU631qprOJESlZNpYNuc7auc04bFbbqBuSKDpREKctSd7LbBWKTVFay0j1z4sbdtsAC6tbTiI8GpxUXHM2D6DrYe30r5Be9NxRPm8cZavaWBITQURotIybTgmj6CfvYQ+IYEEhfUHIkynEuKs7SJ/4NpW95RRKQVorbUsQ+ELMm3ofd/RJCSQ1j/cDvV/lN2wRKWcvF62FNneQWs92HQGIaqqdNsCLGUlBCgHGjtkLJTfY8IjnK1d5IoaSyHMyVjI5qAAYkpKUWUl8uIkKq1pWFOahzXHlmXjxpgbTccRFaCUqg08DLTUWo9VSrUHOmqtZUMa4fGmZLdgtA4g2FKGxRrk3DpdCA9wxtVFtNY7T3cDWgD/qLmIwp3ymvYmIzCADsWlIC9Ooorio+JZnrOcMkeZ6SiiYj4GSnDu+giwB/i3uThClE/q9gP8a20Yn3VIxDLkaRgj78YKz1GeJfxQSvVUSr2ulMoAXgDS3ZpK1Jh5JUGUKUWdiIvlxUlUWe/I3uSV5LHp0CbTUUTFtHXt8FsKoLUuxNkaKITHKii28+jXa2nRoDY3/v3vMPAR+R0mPMrZerI7ANe7bvuBLwElPXy+ZWnmHwB06PcotOhoOI3wdif6srNsdG7U2XAaUQElSqla/DkPpy1QbDaSEGf37582svtQEV+O7Uto8Nm6X4Uw42wj2ek4Z5ZfobUeoLVOBOQ9YB+z8WAauiyEuObtTEcRPiC8djht6rUhNVs2pfEyzwG/AC2UUlOB35C2QOHBZv2RxRe2TO6+oC1xrRuajiPEaZ2tyL4GyALmKaUmKqWGIm8f+pysom3UpgWBVtl6VlSP3pG9WZWzilJHqeko4hyUUslKqf5a69k4X/NvBb4AYrXWKSazCXEmWUeKeOLbP+jWvB4PD5NN1ITnOtvEx++11qOBTsA84EGgiVLqPaXUxTUVULhPqd1OEZlE1WpjOorwIfFR8RTaC9mwf4PpKOLcNgP/cc23eRzYq7WeqbXebzaWEKdX5tA89OUaSsscvDO6J0EB5ZpaJoQR5/zp1FoXaK0/11qPAJoDq3G+GAsvt3zPFrCUEtMoxnQU4UN6Rzh35rZl2wwnEeeitX5Ha90XuBA4AExSSqUrpZ51zcsRwqN8sGAby7Yf5Lkru9C6cajpOEKcVYX+BNRaH9JaT9BaD3VXIFFz5mc4t1Pv11z2FRLVp35IfTo26IgtS4psb+FaovVVrXVPnJPdrwbSDMcS4i+2rPiNwrmvMb79QUb2am46jhDnVOPTcZVSHXGuVHJcG+D/gPrAXUCu6/xTWuufazieX1mXuxGtrQxu09V0FOFj4qLimL5pOsVlxQRbg03HEeeglAoAhgOjgaFACs7JkEJ4hMJtS2kxczQPBtixZv2A2t1RlusTHq/Gm5m01pu01j201j2AXkAh8J3ry28d/5oU2O63K38rgWWR1AkJMR1F+Jj4yHiKy4pZl7vOdBRxFkqpYUqpScBunIMcP+FcM3u01voHs+mEcNJaM3fW1wRoOwE4UGWlzt2JhfBwpmcMDAW2uXaSFDXsqGMn4cEy6VFUv/MjzseiLKRmyVJ+Hu5JYAkQo7W+0jX/psB0KCFO9tmynUze0xxtDQRlld2JhdcwvXr7aJzLRR03Xil1C7ACeERrfchMLN+3KXcPWPNoX182oBHVr05QHbo06sLy7OWmo4iz0FoPMZ1BiLNZm3mYF2ZuZGDHAVgHz4Bdi5wFtrSKCC9grMhWSgUBV+IcSQF4D+eW7dr18Q3g9tNcNxYYCxAREUFKSkqFHzs/P79S13mK6sj/c7Zzp8cGhbVq/Lnw5uffm7NDzeaPLInkt6O/8evvvxJsqZ6+bHn+hfAfRwpLGff5KprUCeGNkd2xhAZBq3jTsYQoN5Mj2cOBVVrrHIDjHwGUUhOBmae7SGs9AZgAEBsbqwcNGlThB05JSaEy13mK6sg/5YdlUAz3XHItzes1qp5g5eTNz783Z4eazR+0J4g5c+cQ1jGM/s36V8v3lOdfCP+gd6Uy97tpROW14qmxY2gQGmQ6khAVZrLIvp6TWkWUUlFa6yzX4d+A9UZS+YntR7eg7I1qvMAW/qNHkx4EWAKwZduqrcgWQviBTBv2ySO4qqyUq4ICCVD9AGkPEd7HyMRHpVQoMAz49qTTryml/lBKrQMGAw+ZyOYvDpbuoEFAtOkYwofVDqxNt8bdZL1sIUSFZM3/GGtZMQHKgVXbZSUR4bWMjGS7Zq83OuXczSay+KPcgiOUWXOJriN7Cgn3iouKY8K6CeSV5FEnqI7pOEIID5e9fgENt05HKecELWUJkJVEhNcyvYSfMOC3bWsB6BHR2XAS4eviIuNwaAfprz9vOooQwsMVFNuZNfMrrDhQgEJBzxtkJRHhtaTI9kOpe5wbhAxu09NwEuHruod3J9gaTNhnP5mOItxAKRWilLIppdYqpTYopeSvKVEpelcqcyc8zua8QJQ1yLkedkAIdL/BdDQhKs30OtnCgE2HNjFyvoVuN7c0HUX4uCBrED2a9AAWm44i3KMYGKK1zldKBQKLlFKztNbLTAcTXsQ10fHyslIuDw7EOvw1KDog62ELrydFth/KObadkUtKsFjkjQzhPrmJSexPTuYR13FapxgAGo8bR3jCeHPBRLXRWmsg33UY6Lppc4mEN9q2fBatykoJUA60tjsL7IGPnPtCITycVFl+pqi0mGK1x3QM4QfCE8YTk55GycJpAGT+8g4x6WlSYPsYpZRVKbUG2AfM0Vqnms4kvMeGvUd4Zk0DylQAWlmdrSIy0VH4CBnJ9iPHRxa/ch3LyKKoCZ0bdWYrYMu2cXH0xabjiGqmtS4Deiil6gPfKaXO01r/ZZ+Dqu7U6+07ZUr+0zt0zMG/lh5DWdqxrMvztCjcyOH653F0WyFsq77H8+bn35uzg+SXItuPhCeM5+NWB/ns0FdMf9lOTHqa6UjCDwRaAkkdHo0tW9bL9mVa68NKqXnApZyymVhVd+r19p0yJf//Kii2M+qDpZToEr6+px8xUZdV6/c/mTc//96cHSS/tIv4k0wbjozJBDv0iWMhasQd17HjyA72Fe4znURUI6VUuGsEG6VULZybjKWbTSU8XZlD88C01aRlHSXphvOJiaprOpIQbiFFtj/JWMjmICvtS0tofF6+7KIlakxclHOFABnN9jlRwDzXTr3LcfZkzzScSXi4f/+0kblp+3juyi4M7tTEdBwh3EaKbD/iaNmf9KAgOhaXEt6jVCaXiBrTsUFH6gbVlS3WfYzWep3WuqfWupvW+jyt9b9MZxIeLNOG7dOnWbtkNrf1j+aWvtGmEwnhVtKT7UfWBDYhz2ohrH5/uPp+WX9U1BirxUpsRKyMZAvhrzJt2D8ewfllJUwLCcTarZ/pREK4nYxk+5F521fz7BQ7LbrfIwW2qHFxUXHsyd/D7rzdpqMIIWrYVtssKCshQDkIxI511yLTkYRwOymy/UjWtm/okgmXFKw/952FqGbxkfEALM9ebjiJEKImLc84yNOr68ta2MLvSJHtL1ZMprjAuUdE/V8fhxWTzeYRfqdt/bY0DGkoLSNC+JH07KPcMXk5ufW7c+yG71BDnoYxP8q7qcIvSE+2H8i4+RaKli/nPoIASJvWFKa9Sq3evxP92aeG0wl/oZQiLjIOW5YNrTVKKdORhBButHVfPjd9aKNWkJVPbo+jXsPa0GGA6VhC1BgZyfYD0Z99StiEexj1pPNvqpjRe4mZ8rgU2KLGxUXFsa9oHxlHM0xHEUK4Ucb+Am6YuAzQTL2zDy0a1jYdSYgaJ0W2n/g5pP2fB1e8A7G3Gssi/FdcpPMtYunLFsJ3ZR4s5IaJyygtczD1zj60axJmOpIQRkiR7SeWZ/0BgKVnDymwhTEt67QkonYEqVmppqMIIdwgd+NCfnr3MdoWb2TKnfF0jKxjOpIQxkhPtp/YdmQTo+YH0/GLL0xHEX5MKUV8VDwLdy/EoR1YlPydL4SvyN24kDrTr+FObWdsQBCWsnhAJjgK/yW/4fxEfukm/r6kADJlZQdhVu/I3hwqPsTWw1tNRxFCVJOdBwr45ttpBGg7AcqBxVEKGQtNxxLCKCmy/UD2prkcCzjqPPjkSim0hVHH+7Jli3UhfMOWnDxGvr+UJWWdsAQEgbKCrIUthBTZvi43MYlDVyUw/WU7AGlTGpA2bAy5iUmGkwl/1TSsKS3qtCA1W/qyhfB2G/Ye4boJy3BoePruW7HcOgNkLWwhAOnJ9nnhCeN5PTKDWYW/Mv1lOzE3HZIXP2FcXGQcszNmU+Yow2qxmo4jhKiElTsPcdvHNsKCA5h6Vx9aNw4F4uT3ixAuMpLtB1YfyyPQHuo8kAJbeIC4yDjySvNIP5huOooQohLmbMzhhonLaBgaxPR7+roKbCHEyWQk2w/klmynTkB7Go+LlwJbeIS4KOfPYWp2Kl0adzGcRghRbpk21iycwfvrG9KpWRyTxsTSKCzYdCohPJIU2T7uQGEedmsO0WEXED5mvOk4QgDQuFZj2tRrgy3bxu3n3W46jhCiHPSuVOyTR3BeWSmfBweih/9AiBTYQpyRtIv4uLlb16CU5vzI80xHEeIv4iLjWJWzitKyUtNRhBDnUGwv4+eZX6HKSglQDoKwE7JnielYQng0KbJ93NI96wAY2ran4SRC/FV8VDxF9iLWH1hvOooQ4iwOFzu4fsIyPspshrYGopUVJUv0CXFO0i7i49IObICyMM5r0tJ0FCH+IjYiFoXClmWjZxP5I1AIj5NpI3vdHL5fHkpaWSfeuGE0gfX7OzeZiR4oc3yEOAcpsn3cvuIt1Le2xWKRNy2EZ6kfUp+ODTtiy7Zxd/e7TccRQpxsxWQcPz1CuKOM91Uge6/5kjZdo4AoKa6FKCepvHxYdt4hSq05tK0XYzqKEKcVFxnHmn1rOGY/ZjqKEMKleMdSymY+jHLYsSpNEKW0yV9tOpYQXkeKbB82a/MKRi2yE99M3ooXnikuMo4SRwlrc9eajiKEyLSx7+eXmPP5O2jtQCnQAMoi/ddCVIK0i/iwdVt/ZuwiTdhDsqOe8Ey9InphVVZs2TYzc8cxAAAgAElEQVTio+JNxxHCbzl2puL4ZAQNy0oZpixYrIHgsKMsFja3G0tHaRERosKkyPZVmTb0gZlAAC2+GwN1ZadH4XnCgsLo0qgLtiwbyBsuQhiRc/QYKd98wbWu5fmsSqHOvxHqtYDogWRtK6Sj6ZBCeCFpF/FBuYlJpA0bw93/df4NlTalAWnDxpCbmGQ4mRD/q3dkb9bvX8/ed94yHUUIvzPrjywueXsB3x9qA9agP5fn634DDHxEBmeEqAIpsn1QeMJ4Aqa/zqgnnUV2zE2HiJnzCeEJsuOj8DxxUXHYtZ0j700wHUUIv3H0WCkPf7mGe6euomXD2rx4/+0E3DYDNeRpGCPvfApRHaRdxEfNzDtpFz15wRQerGeTngRYAgC76ShC+IV1S2ez5Lcf2FXUnvuHXkbCkHYEWi1AnPyuEKIaSZHto1Zkr0VrC6Fj75IXTeGxchOT2J+czOeu47ROzuUmG48bJ++8CFHNco4eY8rXX3PfzofoouyMDQnC0qkfWOVNbSHcQf7P8lE7jqYTVNaUlg8/bDqKEGcUnjCemPQ05k1LAKDp2mXEpKdJgS1ENSootpM8bytD35gPOxYSrOxYcWBxlDp3bxRCuIUU2T7I4XCQp7cTVau96ShClEtcpPPdlpU5Kw0nEcJ3FJWU8eHC7Vzw2jxe/3UT8a0bcsOoG7AEBIOygjVI1r8Wwo2kXcQHLc3cBNYizmvU1XQUIcqlW3g3Xh8YiDXbxpCWQ0zHEcKr5WxYwMalPzNpdzMWHmvDgHaNeWhYB3q1auC8Q/0fnSPY0QOlnVAIN5Ii2wfN3bYCgAujzzecRIjyCbIGkTEyngPZNtNRhPBKWmuWbDvAgt9/4sE9jzIQO/0tgWy95nM6x52y0VMLmeAoRE2QdhEftCpnDdoRzJA2MpItvEd8VDxbDm3hQNEB01GE8BoFxXY+W7aTi99awI0fphKWtYxgZSdAOQjCTufitaYjCuG3ZCTbB2UWpVFHtSEkMMh0FCHKrXdkbwCW5yzn0uhLDacRwoNl2ji08Xe+O9Sat9Lrk3fMTtdm9fjPyO6MaNgQy9TvoKxEeq6FMEyKbB+Tm3+UEstuOte51nQUISqkS6MuhAaGsjxLimwhjssYdTXRCQPJXbCfeo+/yKrFv9Az5VbqOEq5ngAOtXqLQUMv5/yW9VFKAc2deyNIz7UQxhkrspVSGUAeUAbYtdaxSqmGwJdANJABjNJaHzKV0Rv9mL4UpTT9W/Q2HUWICgmwBNArohc26csWAoDSX/9J0bpNHP19Pvu/iOSxsmZ0Ll5HbGApAcqBVZXxSPscOD6h8TjpuRbCI5juyR6ste6htY51HT8B/Ka1bg/85joWFbAwcwV/X+BgRMc+pqMIUWFxkXFkHM0gpyDHdBQhzFoxmf9s+gyA/tHNAagbOYfYYRdjdS3Bp6QdRAiP5mntIlcBg1yffwKkAI+bCuON9hxaxuOLHTQ/uhXqyUiG8C7H18u2ZdsY0XaE4TRCmJFx8y0ULV/O1dQCYPrLdgDufXMz8E9+v/IChlwRh2p9gYxYC+HBTI5ka2C2UmqlUmqs61yE1jrL9Xk2EGEmmneyZywln+3Og0+uhEx52114l44NO1I3qK60jAi/Fv3Zpxx44cITxzGj9wIQkvoz/00axINd0rjfvpP8iM4A5CYmGckphDg7kyPZA7TWe5RSTYA5Sqn0k7+otdZKKX3qRa6CfCxAREQEKSkpFX7g/Pz8Sl3nKU6XP3TGTMJ++olJruO0KQ1gyhjyL7+cghFX1HjGs/Hm59+bs4N35I+2RjN/x3zmlcxzTeT6kzfkPxtvzy9qSKaNdXt/ZTAhzuP+D8K06bSu15p3h77L1LSpvLHiDe6YfQfvX/Q++5OTCU8YbzazEOJ/GCuytdZ7XB/3KaW+A+KAHKVUlNY6SykVBew7zXUTgAkAsbGxetCgQRV+7JSUFCpznac4bf5Bg3i+fwhf5/3A9JftxNx0yDnD3APfSvTm59+bs4N35M9Kz+Kl1JdosXIJ7R59+i9f84b8Z+Pt+YX75SYmsT85+c8CG0hLmE6t3s7J7Eopbup8Ey3rtuThlIe59ZdbecVUWCHEWRlpF1FKhSql6hz/HLgYWA/8CIxx3W0M8IOJfN5qSeFBAspcL8weWmALcS59o/oCUPrhFMNJxLkopVoopeYppTYqpTYopR4wncnbhSeMp9Y3b3LdE1YAYm46RMycT4j+7NO/3C/mu3V89u8CXnliMwBpnWJI6xQjrSNCeBBTPdkRwCKl1FrABvyktf4FeAUYppTaAlzkOhbllF2cTh1rDI3HjZMCW3itVnVbERkaaTqGKB878IjWujPQBxinlOpsOJPX+zh7H/p4q9QZBkzCE8YTk57G3tnOojo5eSgdNq6XthEhPIiRIltrvV1r3d1166K1ftF1/oDWeqjWur3W+iKt9UET+bzRtgPZOAJy6dSgq7zICq+Vm5hEekxn/vvMbkBG5zyd1jpLa73K9XkekAY0M5vK+y3csxDK6tDwvnvPOWAytOVQAObvns97a9+riXhCiHIyvU62qCbfblwIwJDoeMNJhKi846NzO35+AwD7wunEpKfJH45eQCkVDfQEUs0m8W7HSkvIta+jWVBPIu6/v1zXNBp3H1e3u5oP1n3A/Mz5bk4ohCgvT1snW1TS4j02tCOQER2lyBbeLz4qnn3A0qyldA3vajqOOAelVBjwDfCg1vroab5epVWhvH1Vlork/33fZrAW0U61Kf+/uWtXBhaXsDJoJY/Ne4zHmz5Oo4BGlc57Kn96/j2NN2cHyS9Fto/YVbieMNWW0OBg01GEqLJGtRox7aJwVmYtY2y3see+QBijlArEWWBP1Vp/e7r7VHVVKG9flaUi+d/9ag5aW3nmytuJrNPg3BecpHNeZ0bOGMmPpT8yacgkrBZrJdL+L396/j2NN2cHyS/tIj4g8/ABSix76Fivh+koQlSbY7dexep9q9n7zlumo4gzUM6FzD8C0rTWb5rO4wu6zZpLmG5f4QIboHmd5jwV/xSr9q1i0vpJ575ACOFWUmT7gG82LmTUIjtDovuYjiJEtekb1Re7w86R9yaYjiLOrD9wMzBEKbXGdbvMdChvtXLl11y75CiDwlpX+ntc0eYKLo2+lHfXvMuG/RuqMZ0QoqKkyPYBa7b/xMhFmr/VCzQdRYhq0zOiJ4EW+Zn2ZFrrRVprpbXuprXu4br9bDqXV8q0sW7BowDct/UzyLRV6tsopXimzzM0qtWIJxY+wTH7sepMKYSoACmyvV2mjaJC58oidb8YVekXZiE8SW5iEhnnnc/UF4uAP5fyC50x03AyIapfbmISacPG0OeDWgAUfFKftGFjKr10Zb3gerzQ/wUyjmbIsn5CGCRFthc7/sL83JvOTQvSpjSo0guzEJ7i+FJ+i796xHm8eiEx6WkUjLjCcDIhql94wnhCv3mb0afs8liVpSv7Nu3LNe2vYfKGydI2IoQhUmR7sfCE8Sz7792MetK5SEx1vDAL4UmOb7GemiVLLwvfNiFrL45z7PJYUY/EPkLjkMb8c8k/KS0rlQEYIWqYFNle7pcj+0FX7wuzEJ6iU8NO/HhhCEv3LjUdRQi3WrgnBVVWn0b33Vdtr+N1g+ryTJ9n2HJoCx+u/5D9ycnV8n2FEOUjRbaX2573ByGO1jQeN04KbOFzrBYrOdcPZmnWUrTWpuMI4RaHCvM54PiD6Fq9aXJ/QrV+78EtBzM8ejgT1skqPULUNCmyvdiBwjyKLBm0rdNNWkSEz+oT1Yd9hfvYcXSH6ShCuMUHK35CWUr5W4fh1f69cxOTuO3uGXz+onOVkeOTiKV1RAj3kyLbi329fiGjFpVyQUvZSl34rr5NnX3Z0jIifNXsjNlQFsaN3QdX+/c+Pol4x89vAGD75gli0tNkYEaIGiBFthdL3fw9Ixdprm8YZjqKEG7Tok4Lmoc1Z1nWMtNRhKh2hwrzuXC+jVYh8QQFBLjtcYa3do6SJ61OIvNoptseRwjxJymyvVWmjV5zFwDQ4MvRsj628Gl9mvZhefZyynSZ6ShCVKup85IZtdjOjRHt3fo4SilCxo4hwBLAc0ufkzkOQtQAKbK9UOiMmaQNG8OQZbI+tvAPfaP6UlBawM7inaajCFF9Mm1kZHwIwN+XvuD2wZLWDz/Bw7EPY8u28e2Wb936WEIIKbK9kvWCduy4M//EccyNB2R9bOHT4iLjUCjSj6WbjiJEtTi+mdgdic4Wka01NFhybftr6R3Zm/+s+A85BTlufSwh/J0U2V4mNzGJWo+/Q+sP/+zDTpvaiNzvpV1E+K76IfXp3Kgzm45tMh1FiGoRnjCeeW/eXuObiVmUhef6Pkepo5R/p/5b2kaEcCMpsr1MeMJ4il59gAcfdR437loko9jCL/SJ6kNGcQYFpQWmowhRLb45sJmAshDnQQ1uJtaybkvG9xhPSmYKv+78tUYeU/geaVE9NymyvZDN0pD+rn7s8DenyyY0wi/0bdoXBw6WZy83HUWIKjtcVECOfTVNQ/oa2Uzsps430aVRF15OfZnDxw7X6GMLz1eeAvrUHUSl6P5f7lsvSLjNxpwU7l2kOfa3C6XAFn6jR5MeBKpAlu5dyqAWg0zHEaJKJq38FWUp5qr2wwm/ofo3oTmXAEsAz/d7ntEzR/Pa8td4aeBLNZ5BeCa7w87+5GRWjmhH5tFMcgpzyC/N55j9GIGWQIKtwTSs1ZDhwLKsZZzX6DzCgsLYn5xMeMJ4chOT5N11FymyvU2mDUvJHCCA7qHfQ+btUmgLvxBsDaZdcDuW7F1iOooQVTZrx69QVpubegwxlqFjw47c2e1O3l/7Ppe2vpQLml9gLIswa9fRXczdNZdle5exNnctk4HH5j8GQP3g+oQGhlIroBaDZ+cw7LdDJ66rN/g2MoHfLwpnCLB+/3qsrmJbSJHtVXITk9ifnMzdrv9sm6Y0gCljaDxunPxAC78QUyuGbw99y578PTQLa2Y6jhCVcqAwjwEpS7ANHUztwGCjWe7qehdzMubwwrIX+O7K7wgLks3N/MXRkqPM2DaDb7d8y+ZDmxm5sIyHFv05EXb6y3YAGo+77s8a4yrnB4d2sCmmC0W3XEmtT39kyNxcAKwDRgIw65kx9Hnm7Zr7x3go6cn2IuEJ49k64dEan40uhKeIqRUDwOI9iw0nEaLyJs19m1GL7dzR9DzTUQiyBvF8/+fJKcjh7VVSFPmy4z3TmXmZPL/0eS766iJesb1CkCWIx2If48635xKTnkZMehrAic9PV2NYlLN8PP+pV4lJT3POKzhJ9Nc2snv0Y9f01zhSfOR/MvgLKbK9zPTcPaCdkx5rcja6EJ4gIiCCqNAoKbKF98q0kb73MwCuXvR/HrFbb/fw7twYcyNfbvpSJhb7sP3Jyfzf4v9jxHcjmLFtBpe1vowvr/iSL674glu63ELTsKYV+n4nF9bHC/HjBXrg0hlM/fBa/hO/m8u+vYwv0r/AoR3/M1nS10mR7WU2HrERWNrSyGx0IUxTStG/WX9Ss1MpdZSajiNEhRzfgObRt62As+XPU3brTeiZQKu6rXh60dPkleR5RCZRPQpLC3ln1TsA/LT9J0Z3Gs2sa2bxXL/n6Nyo82mvOXVk+nROHeE++Zp2Ddrx7wH/5omoJ4hpFMNLqS9x79x7q/Cv8E5SZHuRTbl7KbHupKW1s7SICL/Vv2l/CkoLWLtvrekoQlRIeMJ4Jjx/qUe2/NUOrM2LA14kpzCHV2yv+N2Ioy/SWrP0hQfZ2bUXF9/wPgBTXizkqls+gY++POu1lfmZPH7NycV206CmvLSpO9NftvNgwgIA0jrFkNYpxi/+kJMi24tMWTsbgH71uxhOIoQ58VHxWJWVJXuX+MWLtPAdDoeDRr8upoHd9ba8h7X8dQ/vzp1d7+THbT+ajiKqaH/Rfu6fdz9jW/7Gs6935diCz4Gz91lXl1O/d5OEBGLS07As+R6AO55tAIu/9Yg/Lt1NimwvsnjvIkbODyC2fgvTUYQwpk5QHbqHd2fRnkUy2ia8yuwln3DtkjyGNujtkS1/uYlJXHT9uydWlTg+4hg6Y6bhZKK8tNbM2jGLq3+4mqV7l/Jo7KN8cfkX9GzS03Q0OjbsCEBoYCj3zr2XzKOZhhO5nxTZXuJYaQkH7asZueQY9fM2m44jhFH9m/Un7WCa6RhClF+mjRWrXwDgwY0TCb/aswpscI5AxqSnEbTMWVQnJw+lU9pGCkZcYTiZKI+C0gIeX/g4/1jwD1rVacX0EdMZ02UMVotzDkB5+qzdrfG4cXww7APKdBn3z7ufInuR6UhuJUW2l/h58WTKrCUAdF/7T4+YkS6ECbmJSQwc9db/jLZJ64jwVKEzZpI2bAzXJjvXxN47paHHTHg8nbb12wIwf/d8vt3yreE04lxyE5PYdHATo2eO5teMXxnfYzyfDP+ENvXa/OV+ntCeEZ4wnjb12vDqBa+y7fA2XrG9YjqSW0mR7QVyE5OIuefPomLz5409+gVaCHcKTxhPx7QN3P1cOFAzPYZC+JtG4+4jPjKeV5e/yv7S/abjiDPQWrM/OZkbf76RgtICPrz4Q+7ufjcBFs/ea7Bf037c2fVOvt3yLSmZKabjuI0U2V4gPGE8tz8WwSsPOgDocMN+j5mRLoQJFmWhb9O+gHPnMSE8mfWCdgA8/KjzOOZGz38Nb5KQwAv9X8CqrHx24DPKHGWmI4lTlJSV8OySZwHo2aQn00dMp3dkb8Opyu/e7vfStl5bXk59mcLSQtNx3EKKbC+wfPdW8gMO0KThCADWdn/B4ybMCFHTBjQbwFcDFOkH001HEeLMMm002PAcAFflu/pPb/3JK17Do8KieCr+KbYXb+fjDR+bjiNOsvOt19jWpTvX3/EVAA/fv4jcngO96h3uQGsgz/R5hr0Fe5n4x0TTcdxCimwvMHmNcxLK3+LvpPG4cRyt18lwIiHM69u0L18NtMruj8JjHd985tinDQC48P1azvPfe8+cmivaXEHP2j1JXpNM2gGZbOwJ0g+mc2eLOdzyTBg7f3kL8N62udjIWC5vczlTNk4htzDXdJxqJ0W2F1iRu5AAe1P6terkdf8DCeEujWs1JqZhDIv3SpEtPFN4wnhaz/qQcY8pwLM2nykvpRTXNbyOBsENeHLhkxSXFZuO5Nd+2/kbt8y6Ba01nwz/hEujLzUdqcrGdR+H3WFnwroJpqNUOymyPdz2gzkUqC10qd/PdBQhPE6/pv1Yu28t+SX5pqMIcVrv7skhN8C5hJqnbT5TXqHWUP7V/19sO7KN/676r+k4fik3MYmpaVN5KOUh2jdoz7Qrpp3YEt0TluarihZ1W/C39n/j6y1fk5WfZTpOtZIi28N9tPInlNJc13m46ShCeJz+zfpj13ZSs1NNRxHitL7aMg3s9ah/7z1eWWAfN6DZAK7reB2fbvyU1Cz5/60m7UtMZH9yMq/YXmFwi8F8dPFHNK7V+MTXvemdkTO5q+tdaK2ZmjbVdJRqJUW2h6v35VSUvT6Xd4w1HUUIj9MjvAe1A2pLX7bwSPO2/8Hwxetpr/oT9cADpuNU2cO9Hia6bjTPLH6GoyVHvWqSnbeyazsHkt8F4LqO1/HmoDcJCQgxnKr6RYVFcXGri/lmyzc+9c6kFNke7MiWFEYs2ku3oHZYLPKfSohTBVoDiY+KZ8neJWitTccR4i8+XPwGIxdp7qwbYTpKtagdWJuXBrxEbmEuL6e+zP7kZNORfNqet9+k2b0JJ46vvXUqmzuf57N/3Nzc+WbyS/P5but3pqNUG6ncPFWmDdsPNwMwLmuO7PAoxBn0b9qfPfl7yDiaYTqKECfs2/QbW+3LAbhk48s+8xreNbwrd3e7m5nbZ5qO4tN2vfU6R9//32XtGo8b5xPtIafTNbwr5zc5n6lpU9mXmGg6TrWQItsDHV/2qfmHYQDUn1RHdngU4gz6NXNOCl6yd4nhJEI45SYmceCq8Ux+1blRki/t0pubmMSg0f89sQNxWqcY0jrF+MS/zRPkJiaxr3Af97dezM3PhPL7W38W1N64RF9Fje40mj35e060yHg7KbI9UHjCeJr99AG3Pe78z+ONyz4JUVNa1GlBq7qtpC9beIxG4+5jzGONef4h39ulNzxhPDHpadRKnQXAO4kX0ilto0/82zzB/uRkbv3lVrIKsnjvovfoUqsL4P0riJTXkJZDqBdcz3SMaiNFtoeamJVLwfE+bC9d9kmImtK/aX+WZy+XNXyFR3jP9hNFAYc5P+pGwDd36Y2uFw3A4r2L+X7r92bD+IgdR3YAcLj4MBMvnnhii3RfbhE5WW5iEtu79GDicwcA33iXRIpsD/Xzjl+gLNTrl30Soib0b9afY2XHWJWzynQUIfgs7RNUWT3GXfK4T+/S22jcffSK6MXrK15nf9F+ry6GTMpNTCKtUwzH4i8D4KPnDxE48LoTz6c/FNjw57sk1qU/ALD8mye9vkWmxotspVQLpdQ8pdRGpdQGpdQDrvPPKaX2KKXWuG6X1XQ2T3GoMJ8c+ypahfTxiWWfhHC32IhYAi2B0jIijPthYyqXLd5Iv/C/ERoc7NUFwrk0SUjg2b7PUmwvltVGqiB79IXc+WxDxr/QFPDeLdKrS4cGHQB8YnKtiZFsO/CI1roz0AcYp5Tq7PraW1rrHq7bzwayeYQPVvyEspRwbcfLTUcRwivUDqzN+RHnyxbrwrgPU//DyEWaf3XoYTpKjWhdrzX3dL+H2Ttnm47ildblrmPsnLHUCarDJ5d+YjqOx9gzagBpB9NOtNB4qxovsrXWWVrrVa7P84A0oFlN5/BkAZ9OhLI63Nh9sOkoQniNAU0HsPXwVrILsk1HEX5qzapv2aXSAGjy1Y0+s2zf2eQmJjFg1Juy2kglrMtdx91z7qZ+cH0mXzqZ5nWa+80Ex3Pp9sSLKBSzdswyHaVKAkw+uFIqGugJpAL9gfFKqVuAFThHuw+d5pqxwFiAiIgIUlJSKvy4+fn5lbquRhxYx+WL9rAkrjtLFi067V08On85eHN+b84Ovp0/qCQIgI/nfUzfsL41mKr8vP35P5VSahJwBbBPa32e6Twm5SYmEZyczDTXcdqUBjBljLNo6uq7T014wnjCE8azYf8GLAP+zlcfX8/I277w21aH8voj948TBfbHl35MZGgk4D/91+fSpHYTekf2ZtaOWdzb/V6UUqYjVYqxIlspFQZ8AzyotT6qlHoPeAHQro9vALefep3WegIwASA2NlYPGjSowo+dkpJCZa5zu0wbP678DxGE8tThRcS2ffi0kx49Nn85eXN+b84Ovp1fa81HX33E/rD9Hvtv9Pbn/zQmA0nAp4ZzGJd340iuC/uIywvzuCUxkJibDv25MpQP/WF1Jl0adyEN+GrzV4w0HcbDrd+//kSBPemSSScKbPFXw1sP5/mlz7Px4Ea6NOpiOk6lGFldRCkViLPAnqq1/hZAa52jtS7TWjuAiYBfLalxfAOa9hNDAQiVDWiEqBClFAOaD2Dp3qWUOkpNx/ELWusFwEHTOTzBU78l4VBlXBH7gvOEny29evx3lbSNnFluYhLr969n7Oyx1A2uy6RLJhEVFmU6lsca1moYAZYAZm333pYRE6uLKOAjIE1r/eZJ50/+SfsbsL6ms5kUnjAe9eVrXPeEFZANaISojAuaX0BeaR6rc1abjiL8SMbBfcT88gNNLHH0j7/e2SLiRwU2/Ln82uHfJwEw+/N7/HqFjNPZn5x8osD++JKPpcA+h3rB9ejXtB+zd85Ga206TqWYaBfpD9wM/KGUWuM69xRwvVKqB852kQzgbgPZjFoy8XN0b1ffkZ+NgghRHfpG9SXQEkjK7hTiouT/H09R1bk0nt7LPnHzRzyzuJToIdHOnF3P+0uLiKfnP5fK5J/0xyQa7W9E86Dm7glVAZ7w/O8q3kVvIEgHMbbeWDat2MQmNp3zOk/IXhVVzd+iqAULChbw2ezPaBncsvqClVNV89d4ka21XgScroPdb5fsA3DsXEbfuWuY3LMljceNkAJbiEqoHVibuKg4FuxewD96/8N0HOFS1bk0ntzLvmfjbLbtcG6CNGZHElxw8f+8fnty/vKoaP7Me+6gQchMZpTMYOpFUwmwGF1jwejzn5uYxP7kZCJcx++9kAM8Xe5dHP3tZ+dU3Y91Z9r0aRwOP8wt599SfcHKqar5ZcdHT5BpY/20awEYdySd8KulwBaisgY1H8TOozu9fn1V4flyE5M4es0DfPyqA3CuKCJzaaDFg4/yZPyTbDywkSkbp5iOY9SRmy5l7HONeeDfLQDZaKaiGoQ0IDYilt92/WY6SqVIkW3Y8QmPgR/XBaDbhNryIi1EFVzQ/AIA5mfON5zE9ymlvgCWAh2VUruVUneYzlSTCm4cxdf9A/jgfudkP5lL86eLW13MoBaDSF6TTObRTNNxjNh5dCd3zbmLQEsgH138kek4XmtIyyFsP7Kd7Ye3m45SYVJkGxaeMJ6Wkx9k7D+cHTQxo/cSM+VxeZEWopKahjWlQ4MOzN8tRba7aa2v11pHaa0DtdbNtdZ+VUnMefwm/r7Yzqg2tzpPyFyaE5RSPBP/DAGWAJ5f+rzXTlyrrL35e7lz9p2UOcqYePFEWtRtIRvNVNKQlkMAvHI0W4psw3ITk9h169tMeM35ApQ2rSlpN70qI9lCVMGFzS9k9b7VHCk+YjqK8FHpv7/BgIXOEdp+KxJpfPPVUmCfIiI0god6PURqdirfb/3edJwas69wH3fOvpOC0gImXDyBtvXbArLRTGVFhkbSrXE3KbJFxYUnjOfBRyNOHMvbjUJU3YUtLqRMl7F4z2LTUYQPyn3xSfR9H544TpvSgP2ffS+DI6fx9w5/p1dEL15f8To733rNdBy3O3jsIHfNvosDRQd4/6L36dSwk8PWAH4AABzXSURBVOlIPmFIyyFsOLCBrPws01EqRIpsw5bYprE38MCfJy59RUZDhKiiro270jCkISm7U0xHET4mNzGJ/Z/976hs45uvlsGR07AoC8/2fZYiexGFH3xsOo5bHSk+wt1z7mZv/l6ShibRLbyb6Ug+Y2jLoYD3tYxIkW1QbmISDW55/sQOWYC0ighRDSzKwsBmA1m0Z5Hs/iiqVfjVcSQ+YP/LuZgpjxP+9MuGEnm+1vVac1uX2wBYkb3CcBr3KCgt4L6597Ht8DbeHvw2vSN7m47kU6LrRdOufjvm7pprOkqFSJFtUOCdt3H9P0KZKDPThah2Q1oOIa8kj+XZy01HET7i+GpQCe/8dd3n3KX5hhJ5vtzEJNI6xTDshvcACB10s09tt56bmESRvYhxv41jw4EN/OfC/9C/WX/TsXzSRa0uYvW+1RwoOnDuO3sIKbINemvJ15RZi+nd8UHnCZmZLkS16de0H7UCajF3p3eNfAjPFXrPXdz6WDgPP+o8bty1SAZGzuH4dusx6WkAjHoyANs3T/jMc7Y/OZkH5z3IqpxVvDzw5RMrYYjqd1HLi3BoBymZKaajlJsU2YY4HA6Cpn6A1d6E6wbf61zaRwpsIapNSEAIFza/kN92/UbOf/9rOo7wAV89fAeFAYe4pvldAIS/OV1etytoYLOBvLvmXXIKckxHqbLjrWhL9i7h+X7PM7z1cMOJfFuHBh1oFtbMq1pGpMg25IcFH3LNkkOMqB+PxWLxmb/qhfAkF7W6iIPHDnLw3fdMRxFebseGn4mbs5LIsmjuuPhBGRiphMbjxvFk3JPYHXbeWPGG6TiVdrwFZmtn58TG6S/b6TTiKZ9pgfFUSikuankRy7KWkVeSZzpOuUiRbUKmjfkbXwfgifRJkGkzHEgI3zSw2UCCrcGmYwhvl2lj0rwEAP6/vfsOj6pM3zj+fZLQm5RQBARRutIRUSyoKALr4gqIq0gT1wJSVMq67lpXcRGEBAuCKyAKroUVUBGF37LWQAQBCSggGHoQFAKh5v39MSdspEibmTMnuT/XNVfmnGl3ksmTZ855z3lf2ppK3IaF2jByGhL79aVqyar0vrA3H6z9gC83fel3pNNStu89vDGhM12Ghcbma6r06Lmm2jUczD7I/PXz/Y5yUtRkR1nOgTN/GhP64/xxcmlNoy4SARlJyay7sCmTn9gNQFqdunnqgCuJjpya/cfkggDse7WUavYZ6nVBLyoXr8zTXz3NweyDJ35ADHHOMTxlOO+uepe7Gt7ld5x8p0FiA8oVKReYU/mpyY6yxH59efovl3LL0HhAZxQRiZScA65Wvx/aa7Rv/hva2iSnrPCdvenxYHn6P2CAanY4FE4ozIPNHmT1L6uZtnKa33FOyZhFY3h9xevcXu927ml4j6ZKj7I4i+Pqc67m0w2fsvfgXr/jnJCa7CjbvGsHqQcWU83VC63QGUVEIuqKKlcA8NG6j3xOIkE0dWB39iRs55Zz7gytUM0Oi6vOuYoWlVowdvFYduzd4XeckzJ+6XjGLx1Pp1qdeKDZA5iZPmz54OpzribrYBafb/zc7ygnpCY7yv4ydxwWt4+eLQbpwBmRKChRsAQL29Xgwx8+5FD2Ib/jSIB8njKVVnOXUiO7Lj11sGNYmRlDmg9hz4E9jF081u84JzQlbQqjvx5Nu3Pb8ZcWf8HM/I6UbzWr2IySBUsGYsiImuwo2rEnk3P+/QYl3QV0rHexPgGLREnF/gPIyMogZbMOMpaTc3DtF7yw6G8ATNg0H9JTVLPDrGbpmnSp3YV/ffcvVm5f6Xec43r3+3d5OuVpWldtzROtniA+Lt7vSPlagbgCXFn1Sualz4v5GX3VZEfR8FmP0/nzfTxQvbXfUUTylcurXE6JAiWYuWam31EkADKSkvm+bS/+PCr0LzJjkg5Qj5R7G91LiYIlGL5gOM45v+McZfba2TzyxSO0rNSSEVeMoEBcAb8jCaEhI0GY0VdNdpRkrplPys4ZANz46V912j6RKCoUX4hrq1/Lx+s+Zs+BPX7HkRi3evsmbh5SkBfuC535Qgc7Rk6pQqXo26gvCzYviKlJRjKSkpm/fj5D5w+lUWIjnmv9HAXjC/odSzw5M/p+si62h4yoyY6CjKRk0tv9ibH/CH1KT3tNW0VEoq19jfbsObgnUFPySvRlrplPqTfeId4VYtAV3kyhOtgxojrV6kTN0jUZsWBEzJwxYtvYsQycN5BaZWqRfHUyRQsU9TuS5FI4oTCtKrfi4x8/jukZfdVkR0HBPr247cGzeHxgNqCtIiJ+aFqhKZWKVdKQETm+9BRenNkDgJE/baJ6ybI62DEKEuISGNp8KBt3b2TitxP9jsM3Gd8AcE7Jc3jpmpcoUbCEz4nkWNqf2z7mZ/RVkx0FQ+e8yP6ETH5Xs39ohbaKiERdnMXRoUYHPt/4OVt2b/E7jsSYtd1uJ61Nd9o/H5ohtOKE4qS16e5zqvzjokoX0aZaGyYsm8Dm3Zt92dObM116wcu6AvD3ISvY1Kil9jrHqMuqXBbzH4DUZEfYxp3bqfTuREq4C7jlqnu0VUTERzfWvJFD7hDvrHrH7ygSS9JTyFoQOoBq6P2hYX3a4xh9g5oO4lD2IUaljmLb2Oif1m/X7e2469Hy3PdEFUDTpceyjKRkVtdvyIRHQ+dYj9UZfdVkR9ijMx+i8+f7ebTm9QD6YxXxUdUSVbnk7Et4+7u3Azeds0RIegqbB3U+vPhshjcxivY4Rl2VElXocUEP3v/h/ai/dvqudPp81AfDeLnNy1F/fTk1OTP67vnPawCsmvVMTH4gUpMdQd8tmcHX++YD0GbeYJ1RRCQGdK7VmS17tvDphk/9jiI+y0hKJq1Nd3Ys/d9BbQcmlQndNl31OtoykpK55pbnefOp0AfgaG2d3JC5gd6ze7Pv0D5eavMS1UtV13TpAdG4fGMAZq2Z5XOSY1OTHSEZSckc6jKYScNDBzvqjCIiseGKqleQWCSRN1e+6XcU8Vlix4tY33vXr9bVvfUnDRPxSc7WyTXvjwAgbcaTEd06mZGUzKbMTfSe3ZvMA5mMazOO2mVqH84isS/O4lj9h2Z8sekLNu/e7Heco6jJjpC0takAfHR3FqDxfSKxokBcAW6seSOfbviU9F3pfscRn+Rsxa4y4dcHTmXE9dYwEZ+1O7cdAGO+HkPm/syIvc62sWPp/VFvdu7bycttXqZe2XoRey2JnEbDnsI5x1vfveV3lKOoyY6A/SkTSJz1JQB37PImvtD4PpGYcXPtm4mPi2fy8sl+RxGfZF5VA4C+DxoA5S7Moly3jiQ+9JSfsQQwM7J7duanvT8xfun4iLzG2pHDAdi+dzsvtnmR+uXqR+R1JPKqlKjCZVUu4+3v3465adbVZIdZxpPDWH37iMPLP75WNrRe4/tEYkb5ouVpf257pq+azs97f/Y7jkRZ5pr5PPzpQABGewc6Jo58Uw12DKk/5DFuOO8GJi2fxPpd68P2vDmn6csa9yoA/3zsFwpcdrOGcgbczbVvZlvWNub+ONfvKL+iJjuMMpKS2TZ5+lHry3XrqGEiIjGmR/0eZB3MYtrKaX5HkSjaMmYM6e3+xEOjQv/+4l89C9CGkFh0X+P7SIhLYGTqyLA9Z4E+t/Hn4XUOL+s0fXnDpWdfSuXilZmSNsXvKL+iJjuMyt7QjJH9D/1qXd3XhmjriEgMOr/0+VxW+TJeX/E6WQez/I4jUbB1TBL3VNhIl2EJ/OcuHS8T6yoUq0CvC3oxZ90cFm5eeMbPt3nMaDY1asnfh6w4vC4Wz60spy4+Lp5u9bqxaOuisLxXwkVNdpis7XY7K6/ryaDR8b9an/FF5A7aEJEzc8eFd7B973amrpjqdxSJtPQUfnr+eb7bO4vzCrXlzhu937mOl4lp3et3p2Kxijyz4BkOZR868QN+w+SW++kyLIEfP3wOgHL33qut2HnITTVvokzhMoxbMs7vKIepyQ4Hb7awPw4pzFMDQ6fsK3dhlraOiMS4JhWa0KpyKyYsm8Cu/btO/AAJpvQU3nkzNOHMtbv38tbFfyCu2sWagTcAiiQUYWCTgaRtT+O91e+d9vOkbkll0vJJdK3dleuqXwfoNH15TeGEwvSo34MvNn3BkowlfscB1GSfufQU9k+8AYDCbj9Dzu8DhA6iUfEWiX39Gvfjl32/MGn5JL+jSASs7XY7aW26U/fl0IQzd4xJ4Pu2vchISlaTFRDXn3s9DRIbMGbRGHYf2H3Kj8922aQ8cT/li5ZnULNBAJpsJo/qUrsLpQuVZmTqSJxzfsdRk30mcs6zuvq10gC88kw2cfdOoEjz5mqwRQKiXtl6tKnWhknfTmLrnq0an5mXeHsZbx5agKcHhPYyagx28JgZQ5oPYVvWNiYsnXDKj0/dncqVc7bQv0l/iiQUAbQVO68qVqAYfRv3JXVLKh+t+8jvOGqyz9Qnz/agy7AE4H/Fu/pkbRETCZKBTQZyMPsgIxaOYNvYsX7HkTBY26UjH07rBECjfft4rsGA0A0agx1IDRIb0KFGByZ+O5ENmRtO+nF7D+5lxs8zAOhQo0Ok4kkMuanmTdQqXYuRC0f6flC7muzTlZ7CtrFjGbf1Dc46WCm0TsVbJJCqlqxK7wt788EPH/gdRcJgbcdryVqykmrjiwEwbFQcm3uN0V7GgOvfpD9xFseo1FEndf+MpGR+uKAxLz2RAcDKuvV1NpF8ID4unqEXDWXj7o2MXeTvRhM12acjPYXNr/0egArZ+3n1ogE6gEYkwDKSkmndNYk3nzoIhE7rpX/GATXnb2StSAfgmQGh36f2MuYNFYtVpNcFvZi9djZfb/n6hPePu+MWev61FI8/1hjQObHzk+YVm9OlVhcmLZ/E4q2LfavlarJPUc5BNDsmhiYwGPMPx/7O9/ucSkTORGK/vtRdkcaOuaHxnv+e1F3/jANobeeOpPV78/Dy4OdCQ/nWLrlYG0HyiB4X9KBC0QoMXzCcbJf9m/d94ZsX2HtwL78v/fsopZNYMqjZICoVq8TDnz3s2zBANdmnwjuI5u4HytBtSOh82DqIRiTvuOTsSwCYkjaFzzZ85nMaORVrbruNrKUrj1pfpE5Vqr959Ey8EkxFEoowsOlAlv+0nBmrZxz3fmt+XsNb371F51qdqVCggs4mkg8VK1CMRy99lLU71/qWQU32Scp4chgZk0Ofhn9J+IUHzmoXukHjsEXylLPu/hPnn3U+Q/87lPW71vsdR07CphUfs29hKl2GJfD+PfsOr6+b1IXq0/0/w4CEV7tz29EgsQGjvx7NngN7jjkU4NnUZymSUIS7G90N6Gwi+VFGUjKlWvf0dRigmuwTyEhKJuPJYWybPJ1tk0JDRF5/+hANh83QQTQieVCl/gN4rvVzHHKH6De332mdl1eiIyMpmYUL3qLPf/sB8Mi2ndx/yV8BKNKgNrR51M94EiFmxuDmg8nIymD80vFHDQX4ctOXzF8/nz4N+lCmcBmfUorfcoYBnv3NlwD87R8XUnP50qh+4FKT/Vu8M4gsLTARgDsHG6CDaETyumolqzHiihH88MsP9Jvbz/fTQMkxePW5WLeH+cezoVX1Xi5K2m3DKdK8uYaI5HENExvSvkb7oyaROpR9iBELRnB2sbO5te6tPqWTWFKqUCkA0ranMW3ltKi+tprs48h4chgb7u8CQKVXQr+gcc+EZg/KsJ7agi2Sx11y9iU82epJFm5eyIB5A9hzYI/fkYTQ1uv9a79g9Ds3A/Dn+x1Fu+8AtAEkP8lISqZ7n38z+YnQnqacoQCfP34fK3esZEDTARSKL+RzSokVZe+9hxYVWzBuybio1vKYa7LNrK2ZrTSzVWY2NJqvnTNOZ+sTQ9k2eTo7lxQ59h1LVo5iKhHxS/sa7Xns0sf4ctOX9PiwB1t2bwHQqf1yiWrN9rZer27bi2tfKAzA35819kwMzbqrY2Tyj5yhAHOnhnb975v/BtWWpvJI3RU0KNeAttXb+pxQYkn5fv3o16Qf2/du540Vb0TtdWOqyTazeGAscD1QD7jFzOqF+3WKzZh5+J/k4X+WXvGe/Pb9DEh8F4DB94e2XNftujH09dZtOpOISD7T8fyOJF2VxLqd6+gyswtz1s05agxofm26o1qznxzGnGk3AdBtSBxpfUJbMHO2XmuugvypZ/2eAAxPGc64JePYmrWVB5o/gJn5nExiTcPEhrSo1IJpK6exZcyYqLxmTDXZwEXAKufcGufcfmAqEN4TXKanUHzWrNA/yYWvsm3sWKa9O5in3u0KQLOH3ufhUaEfyzPPhv5IM/Z1Dj22xywVcZF86PIqlzOl3RQqFqvIoP8bBEDKppTD5+nNx1OxR69mT55OlfHFAZg8PJu6L4dmc8zZeq2NH/lT0QJF2f7Ha1j20zImLJvADefdQOPyjf2OJTGqa+2ubNq9ie3PvxCV10uIyqucvMpAeq7l9UCLsD17egqT/9WJZoSK89XfDCcZaDBsBg04euxW3dt2kGE9SXzoKTg7WQ22SD5W6rUPeWTsksPLJVp3ZyXw6XVn08q/WH6LeM1Oeb0jJQgdF/PJXbu5+sVi1L1tB3R/j4zpKarLwiUPj6FH6kg2ZG5gcPPBfseRGHZl1SspX6Q8sDEqr2fOuai80Mkws05AW+fcHd5yN6CFc65vrvvcCdwJUKFChaZTp049qecuNmMmxWfNOuH9av1xG6vOv4Psx6aTNbw/O0vVOY3vJLIyMzMpXry43zFOW5DzBzk7KH+4VLjrbr67tjG1Plp01G2Z7duz+3cdjvm4I/O3bt061TnXLGJBI+xkara3/pTr9m/V7APXtGB7px5nFj6KYuV9e7qU3z9Bzg6xlf94NSWiNds5FzMXoCUwO9fyMGDY8e7ftGlTd0p+/Mq5x8q55bXruOW16zj3t5Khr49XcG7BP0PXf/zKOefc1jFJp/bcUTRv3jy/I5yRIOcPcnbnlD9clteu85vLx3NkfmChi4Hae7qXU63Z7lTr9pE1+7FEt/WJoSf/+BgRK+/b06X8/glydudiN3+0anasjcleANQ0s3PNrCDQFXgvbM9e9aLQuOocHUaHvnZ/D5r1+NWBMxrfJyLHoymaD4tuze4xMzR8T0QkAGJqTLZz7qCZ9QVmA/HAK865b8P6IlUvIrN9e6pXr+411plqrEXklBxZK/Jr0x31mq3x1yISBtGq2THVZAM4594H3o/ka+z+XQcSr7wSUGMtImcuP9eRaNdsEZEzFa2aHWvDRUREREREAk9NtoiIiIhImKnJFhEREREJMzXZIiIiIiJhpiZbRERERCTM1GSLiIiIiISZmmwRERERkTBTky0iIiIiEmZqskVEREREwkxNtoiIiIhImKnJFhEREREJMzXZIiIiIiJhpiZbRERERCTM1GSLiIiIiISZmmwRERERkTAz55zfGU6bmWUA607joeWAbWGOE03K758gZwfl99uR+as55xL9CuOH06zbee33HjTK758gZ4e8l/+Uanagm+zTZWYLnXPN/M5xupTfP0HODsrvt6Dn90vQf27K768g5w9ydlB+DRcREREREQkzNdkiIiIiImGWX5vscX4HOEPK758gZwfl91vQ8/sl6D835fdXkPMHOTvk8/z5cky2iIiIiEgk5dct2SIiIiIiEZPvmmwza2tmK81slZkN9TvPsZjZK2a21cyW5VpXxszmmNn33tfS3nozszHe97PEzJr4lxzMrKqZzTOz5Wb2rZn1D1j+wmaWYmbfePkf9dafa2ZfeTmnmVlBb30hb3mVd3t1P/N7meLNbJGZzfSWg5R9rZktNbPFZrbQWxeI946X6Swze8vMVphZmpm1DFL+WKSaHVmq2Vbdz/xepsDWbAh23Y50zc5XTbaZxQNjgeuBesAtZlbP31TH9CrQ9oh1Q4FPnHM1gU+8ZQh9LzW9y53AC1HKeDwHgfudc/WAi4F7vZ9xUPLvA65yzjUEGgFtzexiYDgwyjl3PrAD6O3dvzeww1s/yruf3/oDabmWg5QdoLVzrlGu0yYF5b0DMBr40DlXB2hI6PcQpPwxRTU7KlSz/Rf0mg3BrduRrdnOuXxzAVoCs3MtDwOG+Z3rOFmrA8tyLa8EKnnXKwErvesvAbcc636xcAH+DbQJYn6gKPA10ILQyegTjnwfAbOBlt71BO9+5mPmKl5RuAqYCVhQsns51gLljlgXiPcOUAr44cifYVDyx+JFNduX70M1O7qZA12zvSyBrNvRqNn5aks2UBlIz7W83lsXBBWcc5u865uBCt71mP2evF1ZjYGvCFB+b9fdYmArMAdYDfzsnDvo3SV3xsP5vdt/AcpGN/GvPAcMBrK95bIEJzuAAz4ys1Qzu9NbF5T3zrlABvBPb9fveDMrRnDyx6Ig/4wC93tXzfZF0Gs2BLduR7xm57cmO09woY9QMX1aGDMrDrwNDHDO7cx9W6znd84dcs41IrSF4SKgjs+RToqZdQC2OudS/c5yBlo555oQ2i13r5ldnvvGGH/vJABNgBecc42B3fxvNyMQ8/klQoLwe1fNjr48UrMhuHU74jU7vzXZG4CquZareOuCYIuZVQLwvm711sfc92RmBQgV6ynOuXe81YHJn8M59zMwj9DuurPMLMG7KXfGw/m920sBP0U5ao5LgRvMbC0wldDux9EEIzsAzrkN3tetwLuE/mEG5b2zHljvnPvKW36LUAEPSv5YFOSfUWB+76rZqtlnIsB1O+I1O7812QuAmt6RuwWBrsB7Pmc6We8B3b3r3QmNm8tZf7t31OvFwC+5dnNEnZkZMAFIc86NzHVTUPInmtlZ3vUihMYmphEq3J28ux2ZP+f76gTM9T75Rp1zbphzropzrjqh9/Zc59ytBCA7gJkVM7MSOdeBa4FlBOS945zbDKSbWW1v1dXAcgKSP0apZkeYarZq9pkIct2OSs32Y7C5nxegHfAdoTFbD/md5zgZ3wA2AQcIfdLqTWjc1SfA98DHQBnvvkbo6PvVwFKgmc/ZWxHatbIEWOxd2gUofwNgkZd/GfBXb30NIAVYBfwLKOStL+wtr/Jur+H3+8fLdSUwM0jZvZzfeJdvc/4+g/Le8TI1AhZ675/pQOkg5Y/Fi2p2xLOrZsfGeyhwNTtX1sDW7UjXbM34KCIiIiISZvltuIiIiIiISMSpyRYRERERCTM12SIiIiIiYaYmW0REREQkzNRki4iIiIiEmZpsyXPMrKyZLfYum81sg3c908ye9zufiIj8muq25EU6hZ/kaWb2CJDpnBvhdxYRETkx1W3JK7QlW/INM7vSzGZ61x8xs4lm9l8zW2dmfzCzZ8xsqZl96E0zjJk1NbP/mFmqmc3OmWr1N17jilxbYxblzIQlIiKnTnVbgkxNtuRn5wFXATcArwHznHMXAllAe69gJwGdnHNNgVeAJ0/wnA8A9zrnGgGXec8lIiLhobotgZHgdwARH33gnDtgZkuBeOBDb/1SoDpQG7gAmGNmePfZdILn/AwYaWZTgHecc+sjEVxEJJ9S3ZbAUJMt+dk+AOdctpkdcP87QCGb0N+GAd8651qe7BM65542s1lAO+AzM7vOObci3MFFRPIp1W0JDA0XETm+lUCimbUEMLMCZlbfu97XzPoe+QAzO885t9Q5NxxYANSJamIRkfxNdVtihppskeNwzu0HOgHDzewbYDFwiXdzHeCnYzxsgJktM7MlwAHgg6iEFRER1W2JKTqFn8hp8I52/4NX0EVEJMapbku0qckWEREREQkzDRcREREREQkzNdkiIiIiImGmJltEREREJMzUZIuIiIiIhJmabBERERGRMFOTLSIiIiISZmqyRURERETC7P8B+gXDaqS4IIAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mp.plt.rcParams[\"figure.figsize\"] = (12,8)\n",
    "x0, u0, t0, _ = post.get_data(phases=ocp.phase_links[0], interpolate=True)\n",
    "x1, u1, t1, _ = post.get_data(phases=ocp.phase_links[1], interpolate=True)\n",
    "x0o, u0o, t0o, _ = post.get_data(phases=ocp.phase_links[0], interpolate=False)\n",
    "x1o, u1o, t1o, _ = post.get_data(phases=ocp.phase_links[1], interpolate=False)\n",
    "\n",
    "r0 = 1e-3 * (np.sqrt(x0[:, 0] ** 2 + x0[:, 1] ** 2 + x0[:, 2] ** 2) - Re)\n",
    "v0 = 1e-3 * np.sqrt(x0[:, 3] ** 2 + x0[:, 4] ** 2 + x0[:, 5] ** 2)\n",
    "y0 = np.column_stack((r0, v0))\n",
    "fig0, axs0 = post.plot_single_variable(y0, t0, [[0], [1]], axis=0)\n",
    "\n",
    "r0o = 1e-3 * (np.sqrt(x0o[:, 0] ** 2 + x0o[:, 1] ** 2 + x0o[:, 2] ** 2) - Re)\n",
    "v0o = 1e-3 * np.sqrt(x0o[:, 3] ** 2 + x0o[:, 4] ** 2 + x0o[:, 5] ** 2)\n",
    "y0o = np.column_stack((r0o, v0o))\n",
    "fig0o, axs0o = post.plot_single_variable(\n",
    "    y0o,\n",
    "    t0o,\n",
    "    [[0], [1]],\n",
    "    axis=0,\n",
    "    fig=fig0,\n",
    "    axs=axs0,\n",
    "    tics=[\".\"] * 15,\n",
    "    name=\"Second stage\",\n",
    ")\n",
    "\n",
    "r1 = 1e-3 * (np.sqrt(x1[:, 0] ** 2 + x1[:, 1] ** 2 + x1[:, 2] ** 2) - Re)\n",
    "v1 = 1e-3 * np.sqrt(x1[:, 3] ** 2 + x1[:, 4] ** 2 + x1[:, 5] ** 2)\n",
    "y1 = np.column_stack((r1, v1))\n",
    "fig1, axs1 = post.plot_single_variable(y1, t1, [[0], [1]], axis=0, fig=fig0o, axs=axs0o)\n",
    "\n",
    "r1o = 1e-3 * (np.sqrt(x1o[:, 0] ** 2 + x1o[:, 1] ** 2 + x1o[:, 2] ** 2) - Re)\n",
    "v1o = 1e-3 * np.sqrt(x1o[:, 3] ** 2 + x1o[:, 4] ** 2 + x1o[:, 5] ** 2)\n",
    "y1o = np.column_stack((r1o, v1o))\n",
    "fig1o, axs1o = post.plot_single_variable(\n",
    "    y1o,\n",
    "    t1o,\n",
    "    [[0], [1]],\n",
    "    axis=0,\n",
    "    fig=fig1,\n",
    "    axs=axs1,\n",
    "    tics=[\"+\"] * 15,\n",
    "    name=\"Booster stage\",\n",
    ")\n",
    "\n",
    "axs1o[0].set(ylabel=\"Altitude, km\")\n",
    "axs1o[1].set(ylabel=\"Velocity, km/s\")"
   ]
  }
 ],
 "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
}