docs/source/notebooks/multi_stage_launch_vehicle_ascent.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Multi-stage Launch Vehicle (Delta III) Ascent to Orbit**\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: [multi_stage_launch_vehicle_ascent.ipynb](https://github.com/mpopt/mpopt/blob/docs/docs/source/notebooks/multi_stage_launch_vehicle_ascent.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:\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=3, n_phases=4)"
   ]
  },
  {
   "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",
    "m0 = 301454.0\n",
    "mf = 4164.0\n",
    "mdrySrb = 19290.0 - 17010.0\n",
    "mdryFirst = 104380.0 - 95550.0\n",
    "mdrySecond = 19300.0 - 16820.0\n",
    "x0 = np.array([r0[0], r0[1], r0[2], v0[0], v0[1], v0[2], m0])"
   ]
  },
  {
   "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 = [6 * 628500.0 + 1083100.0, 3 * 628500.0 + 1083100.0, 1083100.0, 110094.0]\n",
    "mdot = [\n",
    "    (6 * 17010.0) / (75.2) + (95550.0) / (261.0),\n",
    "    (3 * 17010.0) / (75.2) + (95550.0) / (261.0),\n",
    "    (95550.0) / (261.0),\n",
    "    16820.0 / 700.0,\n",
    "]\n",
    "\n",
    "\n",
    "def dynamics(x, u, t, param=0, T=0.0, mdot=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 = 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 / m * u[0] + param * D[0] + g[0],\n",
    "        T / m * u[1] + param * D[1] + g[1],\n",
    "        T / m * u[2] + param * D[2] + g[2],\n",
    "        -mdot,\n",
    "    ]\n",
    "    return xdot\n",
    "\n",
    "\n",
    "def get_dynamics(param):\n",
    "    dynamics0 = lambda x, u, t: dynamics(\n",
    "        x, u, t, param=param, T=Thrust[0], mdot=mdot[0]\n",
    "    )\n",
    "    dynamics1 = lambda x, u, t: dynamics(\n",
    "        x, u, t, param=param, T=Thrust[1], mdot=mdot[1]\n",
    "    )\n",
    "    dynamics2 = lambda x, u, t: dynamics(\n",
    "        x, u, t, param=param, T=Thrust[2], mdot=mdot[2]\n",
    "    )\n",
    "    dynamics3 = lambda x, u, t: dynamics(\n",
    "        x, u, t, param=param, T=Thrust[3], mdot=mdot[3]\n",
    "    )\n",
    "\n",
    "    return [dynamics0, dynamics1, dynamics2, dynamics3]\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",
    "ocp.path_constraints = [path_constraints0] * ocp.n_phases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step-3: Add terminal cost and constraints\n",
    "def terminal_cost3(xf, tf, x0, t0):\n",
    "    return -xf[-1] / m0\n",
    "\n",
    "\n",
    "ocp.terminal_costs[3] = terminal_cost3\n",
    "\n",
    "\n",
    "def terminal_constraints3(x, t, x0, t0):\n",
    "    # https://space.stackexchange.com/questions/1904/how-to-programmatically-calculate-orbital-elements-using-position-velocity-vecto\n",
    "    # http://control.asu.edu/Classes/MAE462/462Lecture06.pdf\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 = 24361140.0\n",
    "    e_req = 0.7308\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",
    "ocp.terminal_constraints[3] = terminal_constraints3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scale the variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ocp.scale_x = [\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",
    "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": [
    "# Intial guess \n",
    "# Initial guess estimation\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 = 24361140.0\n",
    "e_req = 0.7308\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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Intial guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Timings\n",
    "t0, t1, t2, t3, t4 = 0.0, 75.2, 150.4, 261.0, 924.0\n",
    "# Interpolate to get starting values for intermediate phases\n",
    "xf = np.array([rf[0], rf[1], rf[2], vf[0], vf[1], vf[2], mf + mdrySecond])\n",
    "x1 = x0 + (xf - x0) / (t4 - t0) * (t1 - t0)\n",
    "x2 = x0 + (xf - x0) / (t4 - t0) * (t2 - t0)\n",
    "x3 = x0 + (xf - x0) / (t4 - t0) * (t3 - t0)\n",
    "\n",
    "# Update the state discontinuity values across phases\n",
    "x0f = np.copy(x1)\n",
    "x0f[-1] = x0[-1] - (6 * 17010.0 + 95550.0 / t3 * t1)\n",
    "x1[-1] = x0f[-1] - 6 * mdrySrb\n",
    "\n",
    "x1f = np.copy(x2)\n",
    "x1f[-1] = x1[-1] - (3 * 17010.0 + 95550.0 / t3 * (t2 - t1))\n",
    "x2[-1] = x1f[-1] - 3 * mdrySrb\n",
    "\n",
    "x2f = np.copy(x3)\n",
    "x2f[-1] = x2[-1] - (95550.0 / t3 * (t3 - t2))\n",
    "x3[-1] = x2f[-1] - mdryFirst\n",
    "\n",
    "# Step-8b: Initial guess for the states, controls and phase start and final\n",
    "#             times\n",
    "ocp.x00 = np.array([x0, x1, x2, x3])\n",
    "ocp.xf0 = np.array([x0f, x1f, x2f, xf])\n",
    "ocp.u00 = np.array([[1, 0, 0], [1, 0, 0], [0, 1, 0], [0, 1, 0]])\n",
    "ocp.uf0 = np.array([[0, 1, 0], [0, 1, 0], [0, 1, 0], [0, 1, 0]])\n",
    "ocp.t00 = np.array([[t0], [t1], [t2], [t3]])\n",
    "ocp.tf0 = np.array([[t1], [t2], [t3], [t4]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Box constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmin, rmax = -2 * Re, 2 * Re\n",
    "vmin, vmax = -10000.0, 10000.0\n",
    "rvmin = [rmin, rmin, rmin, vmin, vmin, vmin]\n",
    "rvmax = [rmax, rmax, rmax, vmax, vmax, vmax]\n",
    "lbx0 = rvmin + [x0f[-1]]\n",
    "lbx1 = rvmin + [x1f[-1]]\n",
    "lbx2 = rvmin + [x2f[-1]]\n",
    "lbx3 = rvmin + [xf[-1]]\n",
    "ubx0 = rvmax + [x0[-1]]\n",
    "ubx1 = rvmax + [x1[-1]]\n",
    "ubx2 = rvmax + [x2[-1]]\n",
    "ubx3 = rvmax + [x3[-1]]\n",
    "ocp.lbx = np.array([lbx0, lbx1, lbx2, lbx3])\n",
    "ocp.ubx = np.array([ubx0, ubx1, ubx2, ubx3])\n",
    "\n",
    "# Bounds for control inputs\n",
    "umin = [-1, -1, -1]\n",
    "umax = [1, 1, 1]\n",
    "ocp.lbu = np.array([umin] * ocp.n_phases)\n",
    "ocp.ubu = np.array([umax] * ocp.n_phases)\n",
    "\n",
    "# Bounds for phase start and final times\n",
    "ocp.lbt0 = np.array([[t0], [t1], [t2], [t3]])\n",
    "ocp.ubt0 = np.array([[t0], [t1], [t2], [t3]])\n",
    "ocp.lbtf = np.array([[t1], [t2], [t3], [t4 - 100]])\n",
    "ocp.ubtf = np.array([[t1], [t2], [t3], [t4 + 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, -6 * mdrySrb]\n",
    "lbe1 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -3 * mdrySrb]\n",
    "lbe2 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -mdryFirst]\n",
    "ocp.lbe = np.array([lbe0, lbe1, lbe2])\n",
    "ocp.ube = np.array([lbe0, lbe1, lbe2])"
   ]
  },
  {
   "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............................:      474\n",
      "                     variables with only lower bounds:        0\n",
      "                variables with lower and upper bounds:      474\n",
      "                     variables with only upper bounds:        0\n",
      "Total number of equality constraints.................:      374\n",
      "Total number of inequality constraints...............:      276\n",
      "        inequality constraints with only lower bounds:        0\n",
      "   inequality constraints with lower and upper bounds:      132\n",
      "        inequality constraints with only upper bounds:      144\n",
      "\n",
      "\n",
      "Number of Iterations....: 320\n",
      "\n",
      "                                   (scaled)                 (unscaled)\n",
      "Objective...............:  -2.7222444839176720e-02   -2.7222444839176720e-02\n",
      "Dual infeasibility......:   1.5901114178071890e-09    1.5901114178071890e-09\n",
      "Constraint violation....:   1.9999240413043351e-08    1.9999240413043351e-08\n",
      "Complementarity.........:   2.5890735722630404e-09    2.5890735722630404e-09\n",
      "Overall NLP error.......:   1.9999240413043351e-08    1.9999240413043351e-08\n",
      "\n",
      "\n",
      "Number of objective function evaluations             = 806\n",
      "Number of objective gradient evaluations             = 312\n",
      "Number of equality constraint evaluations            = 807\n",
      "Number of inequality constraint evaluations          = 807\n",
      "Number of equality constraint Jacobian evaluations   = 323\n",
      "Number of inequality constraint Jacobian evaluations = 323\n",
      "Number of Lagrangian Hessian evaluations             = 321\n",
      "Total CPU secs in IPOPT (w/o function evaluations)   =      4.019\n",
      "Total CPU secs in NLP function evaluations           =      0.182\n",
      "\n",
      "EXIT: Solved To Acceptable Level.\n",
      "      solver  :   t_proc      (avg)   t_wall      (avg)    n_eval\n",
      "       nlp_f  |   2.18ms (  2.71us)   2.15ms (  2.66us)       806\n",
      "       nlp_g  |  62.39ms ( 77.31us)  61.96ms ( 76.78us)       807\n",
      "    nlp_grad  | 161.00us (161.00us) 159.72us (159.72us)         1\n",
      "  nlp_grad_f  |   1.81ms (  5.78us)   1.76ms (  5.63us)       313\n",
      "  nlp_hess_l  |  37.45ms (117.03us)  37.34ms (116.69us)       320\n",
      "   nlp_jac_g  |  51.28ms (158.28us)  51.17ms (157.94us)       324\n",
      "       total  |   4.24 s (  4.24 s)   4.22 s (  4.22 s)         1\n"
     ]
    }
   ],
   "source": [
    "ocp.dynamics = get_dynamics(0)\n",
    "mpo = mp.mpopt(ocp, 1, 11)\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............................:      474\n",
      "                     variables with only lower bounds:        0\n",
      "                variables with lower and upper bounds:      474\n",
      "                     variables with only upper bounds:        0\n",
      "Total number of equality constraints.................:      374\n",
      "Total number of inequality constraints...............:      276\n",
      "        inequality constraints with only lower bounds:        0\n",
      "   inequality constraints with lower and upper bounds:      132\n",
      "        inequality constraints with only upper bounds:      144\n",
      "\n",
      "\n",
      "Number of Iterations....: 119\n",
      "\n",
      "                                   (scaled)                 (unscaled)\n",
      "Objective...............:  -2.4977981075384650e-02   -2.4977981075384650e-02\n",
      "Dual infeasibility......:   1.2976399535433416e-08    1.2976399535433416e-08\n",
      "Constraint violation....:   6.8977004524795049e-09    6.8977004524795049e-09\n",
      "Complementarity.........:   2.5143327708305730e-09    2.5143327708305730e-09\n",
      "Overall NLP error.......:   6.8977004524795049e-09    1.2976399535433416e-08\n",
      "\n",
      "\n",
      "Number of objective function evaluations             = 262\n",
      "Number of objective gradient evaluations             = 120\n",
      "Number of equality constraint evaluations            = 262\n",
      "Number of inequality constraint evaluations          = 262\n",
      "Number of equality constraint Jacobian evaluations   = 120\n",
      "Number of inequality constraint Jacobian evaluations = 120\n",
      "Number of Lagrangian Hessian evaluations             = 119\n",
      "Total CPU secs in IPOPT (w/o function evaluations)   =      1.941\n",
      "Total CPU secs in NLP function evaluations           =      0.106\n",
      "\n",
      "EXIT: Optimal Solution Found.\n",
      "      solver  :   t_proc      (avg)   t_wall      (avg)    n_eval\n",
      "       nlp_f  | 801.00us (  3.06us) 779.38us (  2.97us)       262\n",
      "       nlp_g  |  26.59ms (101.48us)  26.88ms (102.59us)       262\n",
      "    nlp_grad  | 183.00us (183.00us) 248.72us (248.72us)         1\n",
      "  nlp_grad_f  | 786.00us (  6.50us) 762.01us (  6.30us)       121\n",
      "  nlp_hess_l  |  36.99ms (310.85us)  37.56ms (315.60us)       119\n",
      "   nlp_jac_g  |  31.55ms (260.72us)  32.24ms (266.46us)       121\n",
      "       total  |   2.07 s (  2.07 s)   2.10 s (  2.10 s)         1\n"
     ]
    }
   ],
   "source": [
    "ocp.dynamics = get_dynamics(1)\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": [
      "Final time :  924.1393 (MPOPT) vs 924.1413s(PSOPT)\n",
      "Final mass :  7529.7123 (MPOPT) vs 7529.7123kg (GPOPS-II)\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(interpolate=True)\n",
    "print(\"Final time : \", round(t[-1], 4), \"(MPOPT) vs 924.1413s(PSOPT)\")\n",
    "print(\"Final mass : \", round(-sol[\"f\"].full()[0, 0] * m0, 4), \"(MPOPT) vs 7529.7123kg (GPOPS-II)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "plot steering data (Thrust vector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3iUVfrw8e+Zkt4mvVdCKAkEEiAISFBpKqDYQUV3QUXc1dVtbnFX3dd13dWfrhVQQFBBsSCuiEqTGlIgJLSQAilACum9zXn/mBBCpISQSUJyPtc1V/K0ee55rjA3pwspJYqiKIpypTQ9HYCiKIpybVIJRFEURekUlUAURVGUTlEJRFEURekUlUAURVGUTtH1dADdydXVVQYGBnbq2urqamxtbbs2oGuQeg4m6jmYqOdg0tefQ1JS0hkppVv7/f0qgQQGBpKYmNipa7dt20ZsbGzXBnQNUs/BRD0HE/UcTPr6cxBCZF9ov6rCUhRFUTpFJRBFURSlU1QCURRFUTpFJRBFURSlU1QCURRFUTqlRxOIEGKZEKJQCHHwIseFEOK/QogMIUSKEGJkm2PzhBDpLa953Re1oiiKAj1fAlkBTLvE8elAaMvrEeBdACGEM/A3YAwwGvibEMJg1kgVpTNy42HHq6afitLH9Og4ECnldiFE4CVOmQWslKY55+OEEE5CCC8gFvhRSlkCIIT4EVMiWm3eiJW+QkqJ8UQ8jZl7aPSKocklgsb6ZpoajDQ2NNPU0IyxWbZ5GU0/jabtonRJYs1xjEYQwvSeQggQpm0hBJTnIpKWg7EJjSYV7dhGtO7B6PQatDoNWr0GrU6g1WvR6TXoLDRYWOnQW2nRW2pN76EovVhvH0joA+S22c5r2Xex/T8jhHgEU+kFDw8Ptm3b1qlAqqqqOn1tX9LTz8Gh/ChOZQcpcwqn3CEMYyM01ZlezQ3nXsZGed52cyMYzx5rAmOzBCmAYUANsPeKYylMOd6Bs+ae+3VjPXCkY28uQKMzvbR60LS8tBagswStpWjzu+l19nenqrTWZ1ThOOiKP9eV6Om/h96ivz6H3p5ArpqUcgmwBCA6Olp2drRoXx9p2lFd9hxy42HX61CZDyMeBI8hcGIHBE6gwW0k1WX1VJXWU1PRQE1FA7WVDdTmn+JkWg0ZzcOpkU7USonReJH3F2BprcPSRoe1jR5LR9PvltY69JY6dAUJ6HM2o6MOvaYBXfh09MNmoLMw/e9fq9eg1WrQaEWblwaNxvT7zl07mDhpIpqWUoKUEglgBIlESiA3Efnx3dDchFFrRfOdH9PsPpzmJiNNjUaaG400N5l+NjUaaWpopqGumYa6JhpbfjbUNdNY20RDfTMNtU3UVTdSfaaR+pqmiz5aS+GNrdYKW20ZtmEabL19sHWyNL0cLbEzWGLjYIHQXH0JR/27MOmvz6G3J5CTgF+bbd+WfScxVWO13b+t26JSOi43vjUxADSk76Kq2Y3qrR9Q2WSg2hhC1eGjVBmLqWp2pqq5gAa5/Wdvo9EJbCzqsTHaYaMtxUWTjU3QEGyGTsDaXo+1vQVWtnpTkrDVY2GpvfQXZG4pfPi0qUiitYDY34Hfz6b6uSiNTqDVnmtCFLTcS9vmpJAx8PDac5/fb3SH3/9yjM1G6qqbqKtqpK66gdqqRuqqGqk9uJXq9BSqm52oNrpQkl5DzcEcpPH8lUe1Og32LlY4uFhh72pt+ulihYOrNQ6uVljZ6lUVmnJZvT2BrAeeEEKswdRgXi6lPC2E+B54qU3D+RTg2Z4KUmmRG48xawfVLuMo14dRkXaY8h2fUdHoSnnzASqaPKiXw1tO/mvrZdaaUuw0xThq8/GxOITdwOHYjZzc+j9mawcLLKy0iLwE+PChc1/6t68HP78LhnJZfqNh3nqzfLn/7D5meG+NVoONgwU2DhZAm0n8AgfCh78994zmrcfoM4raygaqy+qpLm+gqqSOiuI6KotrqThTR0F2BfXV55doLKy0OHnYYPC0bflpg5OHDY7u1uj0WhQFejiBCCFWYypJuAoh8jD1rNIDSCnfAzYANwMZmCqqH245ViKEeBFIaHmrF842qCtm0K4UYTRKKlL3UnY4lXKLwZQ3ulOee5qK7FwqmsIxUgPsB0DDLdhrC3HUFuBhlY69tgA7bQl22mLsNEXYakvQiiZTBb80mr70brkb/Dx/HkdXf+mb6cu9R13gGWkAW0dTMr6YhtomKloSSsWZWsqLaikrqCEvrZS0vfmt5wkB9i5WOHnY4uJtS1ml5ExeFQZPG7Q6zfl/K33t2So/09O9sO67zHEJLLrIsWXAMnPE1a+1+QJo9oymLGUvJV+8TGmDJ6XNdZzSjeTo2q00NwOEAA3oLfJwtK3BRZdNkGUcjrpCHKJvxHHQcOzWzUZjrAONFhCmFmytBUx7GTJ+uGAbyCW/ePril35X68QzsrDW4eprj6uv/c+ONdQ1UV5YS2l+NaUFNZTl11CaX0NeWgnGJsmncfFoNAKDq8Clajsu2uO4WH6P65wXsBk8RlWF9WG9vQpLMYd2JYqmjJ0U24yl9IyR0m1rKWnworQ5g4rmKlNjME8CRhy0RdhaljBwoAbDqS9w1ubgqC/A6sZfI4ImwIdvnqs6ifmt6UvMad159zovSUQ/dH5cKjH0ShZWOtz87XHzPz+5NDcb+fF/PxHiPYQzJ6soTknhVF0Yx4zjoRL4bw02DrtwD3TAPcAe9wAH3APtsbaz6JkPonQ5lUD6urPJwtoFUj+lseAEZ6qcKGoIpqipnKKGIEqaIpDUAqBhOk6607jqjhM6SGIYEIhh99MYNNnodLAv4u+MHDkSPtx5LlkETbh49VL7/w2rJNFnaLUarBwFoaM8CB3lAVHF8OFM6hr1lBhDKBrxD4oqnCk8UcGJ1DPQ0o5v72LVmkw8AhxwD3RAb6naVa5FKoH0VbnxyORPKEnYQX5dMIWNAyhovJOSJn9kS1cha00ZbrpMAm3jcdOfwDlyJI7py01VTloLmLXe9IU/9M3WxFCRWdPxZKH0Ly1/F1YnduAdOAHvNn8LDbVNFOVUUpBdQeGJSopyKsjcVwiARiNw9bPDa4ATXiGOeIY4XrK9Ruk9VALpS3LjqTgYR26xF3lJxzhZP5Fa40wALEUlHvp0gmwTcNNn4q7PxFZfaaqfPtsuMe55GDfr0okhc9vP9ynKWRf5u7Cw1uETZsAn7NyMQ7VVDRQcryA/s5zTmeUc3H6SA5tN44Md3Kzxbkkm3qFOOHnYqLaUXkglkGtZS/VUmcM4Mo80kbnrKEWN4QDYaMLws9iPr0UKXhZHcNTmc96/v4Dr4KbnTb9fKGEoiplZ21kQGOFKYIQrAM1NRopyKjmdUc7pzDJOHCzmaJypB5itowU+gwz4hjnjO8iAvbNVT4autFAJ5BrVfGIvWUv+yaGqSZxsqAHAXd/MWLsPCbRKwqA/jUCaej+NuB88h0Pqp1B6AiLuhsnPn3szlTCUXkCr0+AZ7IhnsCMj8EdKSXlhLSePlZKXVkru4RKO7S0AwNHNGt9BBnwHOeMT5mRqmFddiLudSiDXmObjezmy5ShJqa5U1f0ae00hY+w+IWyMG/Zpy881bE97FWqLz//H1L7Xk6L0YkIInDxMAxiHTvBBGiUlp6vJO1pK3tESjiUUcGjHKQDcPDX4V39BgEUiHlb/QfPQ1yqJdAOVQK4RUkqOb9rDjnW5VDX74WmRxkSXdfjrk9DodBCzHmJuV/8DU/osoRG4+Njh4mPH8Bv9MDYbKcyuJO9oKTl7UthXOZMkbsdSVOL/US4Bk/zwG+LSMlpfMQeVQHq73Hhqjuxha0okJzKacdHVcIPhb/haHkJEPwiOU1T7hdIvabTnqryiIwqpW34fuTWDyGmIJvtMLOkrjoAAd397AsJdCBzmipu/vWqM70IqgfRmufHkLf4jP5Ysot5Yy3Wjqxl+8s/nutkOn6MShqIA+I3G6uHVhJ7YQWjgBKTPKIpyK8k+WEzOoWISNpwg4dsT2BksCYp0IzjSDe8Bjmi0Pb2m3rVNJZBe7MjmQ2w98yxO2tPMdH0Rl9AHYPI6VU2lKBfSpguxANNgxQAHRt0SRG1VAydSzpCVfIbDO0+RujUPS1sdQRGuBEW64TfEGb2FGsx4pVQC6Y1y49m3IY09+wPws0xhmtMrWOiN55KGShyKckWs7SwYfJ03g6/zpqGuidwjJWQlF3E85QxH4/LRWWjwH+JCSJQbgRGuWFipr8aOUE+pA5KyS/lfZgP2QaVEBZh56fXceFLffpM9ZQ8Tar2bG+/1RdvwjCpxKEoXsbDSETLCnZAR7jQ3Gzl1rIys5KLWl06vISDChQFRHgRGuKBTJZOLUgnkMpKyS3npHyupQcc9mY3cM8qf2SN9zZNIcuPJ+Hwt28vmEWgZz02Or6NpeBYmPNP191IUBa1Wg99gZ/wGO3P9PQM5nVlGemIhmfsKydxXhM5SS9AwV0Kj3fEf4oJWr9pM2lIJ5DK+2JfHPYe/J6L4OKkuQawqmMraxFDuivbr2kSSG0/x+wvZXPginvo0phr+z9Q99+wstoqimJXQCLxDDXiHGphwdygn08vISCwkc38h6QkFWFjrCB7uysAYT3ys09Hk7Oz3NQMqgVyGAP46dj7TsuO5I2Mbr+x6j62+I1hSM5Mv9uXx8fyYLkkiDd//g41nnkQvapnm9Cq6kOsg9tl+/cepKD1Fo9XgN8gZv0HOXH/fQPKOlJKRWEBmchFH4/Kx0xYz0OoEYXaf4Dz/3Z4Ot8eoBHIZs0f6sjYpj69DJvBdYAx3H9vCXelbiTiTyb9HzSUuK/TqE8gXC4g/HExZsy+zDH/FVl+ukoei9BJarYaAcBcCwl2Y2NDM8U9XkJZYwv7qmeyrvh23dwvRBthSM7Kh3w1aVAnkMqICDKxeEMPqTQmMDI/g4KkQfrd5GL/f8yEv7VxM3iBb3gZigl06l0hy4ylI2s+BmpcJt/kOX8uDYOutkoei9EI6Cy2hsRGEZs+kpsGa9PqJpFk8RP5+IysO7MJ/qDNhYzwJGu7aL9aOVwmkA6ICDFSGWBA7xh+ApJG+JBwajf0HLxOw9D9sGJbBm2HXd6o6y5i1g60VC7HVlDDWbpVp57C7u/ojKIrSVVrWPbE5sYPhgRMY7jeajV9vxdHoR9reArJTD2Fpo2PgKA8Gj/P+2UqOfYlKIJ0QFWAgKsDAO7q/YvXPv7IwZR2VFrZXXp2VG8/RQ1DcFMRUp39joan7+Uy5iqL0Pu3GY1k5CsbGDmDMrBBOHivl6O7THN51mtSfTuLmb8+QcV6EjvLA0kbfg0F3PZVArsKYME8eGvsgdtsX85ukNRyfPIK3t3awOis3nsbldxKf/yoe+nRCxg2CyOdV1ZWiXMM0GtHa+D7hnkbSEwo4tPMUP60+xq7PMwgZ6c6Q8V54DXDqE3Ny9WgCEUJMA94AtMD7UsqX2x3/P2BSy6YN4C6ldGo51gykthzLkVLO7J6oz4kKMLDi0QkkjvOn+d+/w/X/XuD5Sb/hTTvHy1dnndhBSuVkqo0uTDG8jnCarZKHovQhVrZ6ImJ9CZ/oQ1FOJUd2neZYfD5pe/Nx8rBh8HVehMV4XtPL9/ZYAhFCaIG3gclAHpAghFgvpTx89hwp5W/anP8rYESbt6iVUkZ2V7wXc7Y6a0XJ74n859M8k/gJfxu3gLis4ksmkEbv8SRXDyDAMglv60w13kNR+ighROu8XNfdOYDMfYUc2XWaPV9lsvfrLIJHuBE+0Qfv0GuvVNKTJZDRQIaUMgtACLEGmAUcvsj59wF/66bYrljE9dF8sPF2Hk9ayy058cQEj7/k+YezvakzVhN1vT2MXq9KH4rSD+gttAyK8WJQjBel+dUc2nmKo7tPk5FUiMHLlvCIRsIcE7AcOO6a+E4QUsqeubEQdwLTpJTzW7YfAMZIKZ+4wLkBQBzgK6VsbtnXBCQDTcDLUsp1F7nPI8AjAB4eHlFr1qzpVLxVVVXY2dld8pyMkiYC33sTj4IcSv/2V4zOzhc8Txol6f+T6G0h6MZra2qEjjyH/kA9BxP1HEyu5jkYmyTlOVCeVkt1uRU6UUeo9U70I0Jp9gvq4kg7Z9KkSUlSyuj2+6+VRvR7gc/PJo8WAVLKk0KIYGCLECJVSpnZ/kIp5RJgCUB0dLSMjY3tVADbtm3jctfGAg2jB5E1YyZBmzbj987bFzzv6Le7OFxTz5QZVgTGXtepeHpKR55Df6Ceg4l6DiZd8hx2vErhxjUcrJlMes0EmnZZ4hFkR/hEHwZEuffKcSU9+d/fk4Bfm23fln0Xci+wuu0OKeXJlp9ZwDbObx/pMRa+vrg+vpCqLVv45J3PScouPf+E3HhSNx7GoMslYOdsyI3vmUAVReldAifgbp3LDU7vMc/7CcZPtqC+ponNK46w8k+72bs+i+ry+p6O8jw9WQJJAEKFEEGYEse9wJz2JwkhBgEGYE+bfQagRkpZL4RwBcYBr3RL1B2QHTuDyndX4LTiHR44Zc2qBde1NqgXJiVS2DiI6+2XIIwNpsWhroG6TkVRzKxlgCIndmDVMkBx2GxJ3tFSUrbmkfjdCfZ9n82AaHeG3+CHe4BDT0fccwlEStkkhHgC+B5TN95lUspDQogXgEQp5fqWU+8F1sjzG2sGA4uFEEZMpaiX2/be6mlxeVXsDb+FZ+NXEZsRR1xWWGsCOXg6HJ2oIsxmh2lZWtX7SlGUs9oNUBRCtE43X1ZYQ+rWPI7sPs2xvQV4DXBk2CQ/giNde2xp3h5tA5FSbgA2tNv3XLvtv1/gut1AhFmDuwoxwS686R/JkYzt3Je2CRvfRQDUVTeSftjIwEg7LIY83e+nglYUpeOc3G2YcM9ARs8M5uju06RszeX7pQexc7ZkWKwfQ4JOY5nfvVPMXyuN6NeUqAADHy8YS5pPLa5v/AWPvZshdC5pcfk0NRoJvzkK/K7v6TAVRbkGWVrrGH6jHxGTfDmRcoaULbns/jKDBFFDuE02wx2WYvuLld2SRK6tPqTXkKgAA/c9Nhvr6CiKFy/GWFfHkT2ncQ90wM2v706upihK99BoBMGRbtz29Ejunn6UQMt9JFfPYOXpN9iy5gSl+dXmj8Hsd+jHhBC4/erXNBUVkb3iS4rzqggb49nTYSmK0se4jYxmius7zHX7FUNst3Is25VPnt/Ld++lkh+3B3a8apYen6oKy8xsx4zGZswYUjceAU930B7F1GNZURSli7T04HI8sYOJgRMY5RhJytZcDm7NJisZvC0smGR4BKf5S7q0akuVQLpB9cRoCpyG4dWUSujm+ziasKmnQ1IUpa/xGw0TngG/0dg4WBAzK4QHZx5gnP0KapsdsaLENGygC6kE0g1OVlVTZ+WCW24SOmMTpYe39HRIiqL0Axah44h0/IH73J7CSt/Q5cMGVBVWN6g2RqOV9TjlHaTcxwbDzBt6OiRFUfqDlqotcWKHWbr3qgRiZs1NRs7k6nENkYgUHUWnwxkUfWNPh6UoSn/RbnBiV1JVWGaWc7iE+uomoqdH4vPU7+F4DjV79/Z0WIqiKFdNJRAzS4/Px8pOj98QZxxvm4XWzZXiJUt7OixFUZSrphKIGTXUNXH8wBkGRLmj1WrQWFri/OCDVO/ezYcfbvz5TL2KoijXEJVAzOh4chFNjUYGjj43eDB73DSqdVbUr1zO3PfjVBJRFOWapRKIueTGc+z7vdg7CjyDz027HFdQz/+CxzHuZCrupfnEZRX3YJCKoiidpxKIOeTGU7PsfnJP2zLQ+BUiL6H1UEywC9+FXU+DVsfdGduICXbpwUAVRVE6TyUQczixg/Tq0Ui0DLTcdt7oz6gAA+88MZmiCVO5ITeJYfranotTURTlKqgEYg6BEzhWez2uuuM4WxX+bPRnVICBSc/9BoCtL76u2kEURbkmqQRiBmWW4RQ2DmBghIVpicoLDOJJabRmi99I3LZvZOHbm1QSURTlmqMSiBmkJxSAgNC7777oCNC4rGI+GzAJi+Ympqf9pBrTFUW55qgE0sWklByLL8BnoBN2BsuLnhcT7EKhwZPd3hHMyNzFmVNFqhSiKMo1RSWQLlaUU0lZQQ0DR1164aioAAMfz4+h9K4HsW6qR7NmpRoXoijKNUUlkC6WtjcfrU5DyEi3y54bFWDAdlAY2/xGMCNzJzYVpbz83ZFuiFJRFOXq9WgCEUJME0KkCSEyhBB/vMDxh4QQRUKI5JbX/DbH5gkh0lte87o38gszNhtJTyggcJgLljb6Dl0TE+zCx4OnopVG5qRtIuFEKS9vUElEUZTe74oSiBBCI4RwuPyZHXovLfA2MB0YAtwnhBhygVM/lVJGtrzeb7nWGfgbMAYYDfxNCGHoiriuRs72OGorGxkYXNnha6ICDBg9vfkuMIZp2XvxripiXfJJM0apKIrSNS6bQIQQnwghHIQQtsBB4LAQ4nddcO/RQIaUMktK2QCsAWZ18NqpwI9SyhIpZSnwIzCtC2LqvNx4jq3fjKWoJGDPXVe0gP1tkT6sDruJeq2ex1K/pqiyXrWFKIrS63VkQakhUsoKIcRc4Dvgj0AS8O+rvLcPkNtmOw9TiaK9O4QQ1wPHgN9IKXMvcq3PhW4ihHgEeATAw8ODbdu2dSrYqqqqS17rnbWO47VTGWS9BdFcQ9aWleQE1HTovWNsYIunE6sGTeXRg+uJOZnK6k16KkMsOhWrOV3uOfQX6jmYqOdg0l+fQ0cSiF4IoQduA96SUjYKIaSZ4zrrG2C1lLJeCPEo8CFwRevBSimXAEsAoqOjZWxsbKcC2bZtG5e69mi1nibqCbPZgUZnSfANDxJ8BauA2QeV8kBlM1NyEng09Wts//QAUYMumBN71OWeQ3+hnoOJeg4m/fU5dKQNZDFwArAFtgshAoCKLrj3ScCvzbZvy75WUspiKWV9y+b7QFRHr+1uaek2OBgEHlPvuujo80uJCjCw6pFxlDzyFK515ZS++h9VjaUoSq922QQipfyvlNJHSnmzNMkGJnXBvROAUCFEkBDCArgXWN/2BCGEV5vNmcDZ7knfA1OEEIaWxvMpLft6RMWZWvKOljJofCDi+mc6vf5wVICBwZMn8HXoRHx2bOS1F5epJKIoSq912SosIYQH8BLgLaWc3tJTaizwwdXcWErZJIR4AtMXvxZYJqU8JIR4AUiUUq4Hfi2EmAk0ASXAQy3XlgghXsSUhABekFKWXE08l/LZ0c9YW7CWfYn7cLB0INojmkj3yNbjh3eeQggYfJ3XJd6lY+Kyilk+aBrD89NYlPApSQdiiQqIuvyFiqIo3awjbSArgOXAn1u2jwGfcpUJBEBKuQHY0G7fc21+fxZ49iLXLgOWXW0Ml7M2bS1blr5DULWOfalfcsq1jnedJbNCZzEzZCYRLsM4svs0ARGu2Bmsrvp+McEuvGlpwavR9/Hatv8SufzfvGPzPGNC3YkK6PGeyoqiKK06kkBcpZSfCSGehdaSQ7OZ4+o1NuVsotyuAadKPZHpjoxId6LKqomMtB95MuBbFvj8iZoKe6wj6rrkfmenOInLKuZ0iIaApf+h9PXXmDviNj6eH6OSiKIovUZHEki1EMIFkABCiBig3KxR9SI3+d/EC8G7ORRciWWDBp8ia0JO2jIsw4GITDhh+x16qxj+ePwtloS9f17VVmdFBRiICjDwNnAgOI7ZGdvJtfckLitUJRBFUXqNjiSQpzE1bocIIXYBbsCdZo2qF7kr7C4APkv+jLFDx1LdWM26jHVY1sCIY24E553GWP0FA9Mt+cZ9KYk+I37WRtJZMcEuPDhiFr5VhTyxfy1N+SOAAVf9voqiKF3hsglESrlPCDERCAMEkCalbDR7ZL3IXWF34XbajdjoWABmhMwgsSARTUkzlVWunNC/x4hjBkry01k+Yg+LHbQsnbL0qpNIVICBlY+MI2F8ACx5HquXnmO/pRW7DQOICXZRpRFFUXrURROIEGL2RQ4NFEIgpfzSTDH1epHukQyyGcqHJ3fg7bANW580Clydsc3w4tZdnuyKLCGxILFLq7Oax77P0TkPov3Lb9kTPZc3AyJVm4iiKD3qUiWQGZc4JoF+m0AADm7Po6kJbrhvOi7VjiQbvHhq9xuMTXRkYqILGq983uf9LqvO0jo5sfPxF3B56Vn+GL+KxfWVxGUOUAlEUZQec9EEIqV8uDsDuZY01jeTsjUP/6HOuIyMBMYSCbzuPoSEqDga1iVT+k0c+4+Ws3jQ4i6pzgKIjgjg4YmP8cyeD3n8wFecWVPLu00LGB3mpRKJoijdriOz8boIIf4rhNgnhEgSQrzR0iur30r9eg+1lY1ER9eetz/SPZIFUY9hcedIMnyqGZHuyMBjViQWJHbJfaMCDCx/7Hqq/vpP8mfOwXX7Rrz+9hS//b9v1Ih1RVG6XUfmwloDFAF3YOp9VYRpIGG/1JCxl/3bivG33I/X5tsvOG37KK/RJI2oIsu7mpFpjrgfqr/AO3VOVICBRTeGcXj6fbw45iG8q87w2o//IXf5SpJOlPD21gyVTBRF6RYd6cbrJaV8sc32P4QQ95groN4u5Yd06ozejLZdDc0NcGLHz+a+inSPZMnUpSQMi8dqQybpX33HRgc78vwau7SL75v+w3jC4MeT+z9j5Jr3SN38I2uH3c6bzt6qgV1RFLPrSAL5QQhxL/BZy/ad9ODEhT2pqrSOpCPeBFvF42GZBVoLCJxwwXMj3SOJdI+kaUgjH/79N6Ss+oxNo4pY7N41bSJtR6wPemoaR1d+QuDaZby15TW+CR7P+l0G4rJcVXdfRVHM5lLdeCsx9bYSwFPARy2HNEAV8FuzR9fL7P4oHtlsZNxsH9D92ZQ8LjPzrk6vxzhrCJWL04lNcuH764q6vIsvAPPm8liVB/elbGBm5g6q/5nAF6GxvBc6gZkxA5g90lclEkVRutSlemHZd2cgvVpuPE6JX7MrYwbRdl/isOvLK1rzY1TgdXw85kNu2unMjQluREwN6/IQowIMvLfoJuKyRvDd4SO4rl3OQ4e/Y1bmDj7LvIGH42IYFORBfZORe0b5M2eMf5fHoChK/9KRKixa1twIBVqnm5VSbjtiJlQAACAASURBVDdXUL1KbjwNy+7iQMG/cNSeYqTtWmhuvmDbx8VEukfy+m2L2Ru4iZpVuzjw3ipS7jvEKL+YLimJnHW2RJIU7MLcfB1BBVk8eHgjjx5cz5y0H9kYOIb1weP5U14529IKeXRiiCqVKIrSaR1ZD2Q+8CSmVf+SgRhgD1e4tOw168BqdpbNodrowmznP6MXjaC1umjbx8VEukcSOSWSTTpf9i9exckVmSwd/QFLpnbNGJG2zraPfLHPl7+5BxNSmMVtGduZnf4Tt2XsYLN/FJ9XTWJuepFqbFcUpdM6UgJ5EhgFxEkpJwkhBmFaYKqfkLjps7DXnsHTIg18omDay51edfCEWxUJQ8sYc8jA8NRmEoYndHkCgXOlkTtG+vLFPl9ejg/EvaqY2RnbmZq9lynZCezxDme7cxVxwyOJCXFViURRlCvSkQRSJ6WsE0IghLCUUh4VQnR9JX5vNXwOEftvwdjcaOp1dRXJAyDaI5rFwYuxq6lk6HF7inbsJ9kz2SxJBM4lknBvR/6yLpV3h9/O6rCbuC1rBzcf34P94udIM/jx7sBYFv51PlED3M0Sh6IofU9HEkieEMIJWAf8KIQoBbLNG1Yv4jcaHvqWE1tWEnzDg1eVPMBUlbV0ylK+CV5Pzqc/4bc9m79XL2LkhKnMDJlptkQyZ4w/YZ72xGUVY7CxoLQmmv8VPkTpl18xM2MHz+xdRf2c/5Ey83b2Dr2ekSMHqhKJoiiX1JHp3G9v+fXvQoitgCOw0axR9TZ+o8kJqCH4KpPHWZHukSQWJPLV8GIm17kxdr8jP+jX88WxL5g3dB5PRz/dJfdp77xuv0BSdilzD07gu6CxjCpO58mKA9h/tJwYsYI4n2FonlpA5C2TEEKYJR5FUa5tlxoH4iClrBBCOLfZndry0w4oMWtkfVy0RzRaCz1booqYHufBTYnu/DiqkOWHluNn79e6kJU5tR2MGBM8nrisYj7+ajfTs3YzNTseq98u4viSgVRMv409QaMYPdjH7DEpinLtuFQJ5BPgViCJcwMK2/4MNnt0fVhrVVbmN3wjv2BanAeTE0xJ5KMjH3VLAoGfl0redHRjRcQM1kZMZ2VgOcb1n2P5xitE663YFDgan5mxENstoSmK0stdaiDhrcJUdzFRSpljjpsLIaYBbwBa4H0p5cvtjj8NzAeaME3i+AspZXbLsWbOlYhypJQzzRGjOZ2d7iSzLJONJDMtzoMp8e7saMwnudB8DesXc36JxIXIAANvh45lw5ofuCVrFzdn7ET36nZy9m6hePIsdruGEjPAXbWVKEo/dck2ECmlFEJ8C0R09Y2FEFrgbWAykAckCCHWSykPtzltPxAtpawRQiwEXgHOTuRYK6Xs3m9YM3kq6inmFc1j49h8bkx0Z1KiK6u1/0bO/S0jPEZ0ayztSyQxIa686RHMqy6BrBwxi7/XJuC3fy9WO3cyxMaZTweMh+ceJyo8oFvjVBSl53VkOvd9QohRZrj3aCBDSpklpWzANG38rLYnSCm3SilrWjbjMA1m7HMi3SP5y5i/0GilYWNMIdkeNbjvLWfFP39Hwok9PRrb2VLJ01PCePuJKVjdeSs/Pb+Uf45+gCJrRx5OWY/lnFmcfv559u3Yr6aTV5R+REgpL32CEEeBAZi67lbT0gYipRx2VTcW4k5gmpRyfsv2A8AYKeUTFzn/LSBfSvmPlu0mTCPjm4CXpZTrLnLdI8AjAB4eHlFr1qzpVLxVVVXY2dl16tqOOl5/nA1lG0irTWNolj0jjznRbK2j5roAgvxHEGQZZNb7d0RVVRX5jda8klBHkxFCK07ybPke3A4komlqIsl9IN+GjOfGWyMZ4KLv6XDNpjv+Hq4F6jmY9PXnMGnSpCQpZXT7/R1JIBesmzjbFtFZV5JAhBD3A09gao+pb9nnI6U8KYQIBrYAN0opMy91z+joaJmY2LnVAbdt20ZsbGynrr0SyYXJLPhhAY3GRtzLrBi33wm7Gh1ZfrXMW/gCFva2JBYkdtm6Ilfq7HNIyi5tbSuJCjCwZH0Sx1d8xM1Zu3Gtq6DGzYugBQ9zPHoScQV1fW5a+e76e+jt1HMw6evPQQhxwQTSkXEgZxut3WkzmWIXOAn4tdn2bdl3HiHETcCfaZM8WuI62fIzSwixDRgBXDKBXAvO9s5KLEjkVNUp1tl9SUS6PeHHHdj+/CukBpZzMKAMjaW+y9Za74z2bSVRw4N5begUvgidxPX5B3mqYj8FL70EulepDBjFMwOv59Wnbu1TSURR+ruOTKY4E3gV8AYKgQDgCDD0Ku+dAIQKIYIwJY57gTnt7j0CWIyppFLYZr8BqJFS1gshXIFxmBrY+4SzvbOSC5P5JvMbDgyuJM+/kSk5AxiW1szA49YcCaziG68vSXTvudJIW+f34JpAeICBlcs3UPnJx9yctYcZmbvIy/6Rj6fNZtCU64kKdL78myqK0qt1ZCqTFzHNwLtJSjlCCDEJuP9qbyylbBJCPIFpdUMtsExKeUgI8QKQKKVcD/wb06DFtS2joc921x0MLBZCGDF1BHi5Xe+tPqFtaSTaw1R6/MOnCxmaZsuIY440ZySR5LWDzwOX8c973kYI0aPVW+1LJUNvGMvcLMGKylu49cQepmXsxv+1BDKW+1J491zWOQ0iv7pJrU+iKNeojiSQRillsRBCI4TQSCm3CiFe74qbSyk3ABva7Xuuze83XeS63Ziha3FvdLY0cta/7nmXxIJETudkkrNtNyEnbRlw0pZNqf/giEcxGZ6VLHbU8PvRv6e8obxHSydtSyWnyiJ4aPcNxOYkcXvmDvwW/4v7rRxZHzyOl7JiyCmuxt5a3+faShSlL+tIAikTQtgB24GPhRCFmHpjKT2gtXrLI5kFpZtIHlRBYKEd0SV+DE63ZWi6HZU2Tfxw8G3y3GpZ7PouM8JmMdh5cI8klNZFrrJL+WJfHj8Ej2VTcAyRp49yW+Z2fnF4A3PSfmTTkVGsDpnAmwYPtUaJolwjOpJAZgG1wG+AuZgmU3zBnEEpl3eh6q0nvnkUr1M6fAutGJBjw6ATdjRqjRQm/MQ65+/Jd67nLccGxviOZZTXqG5NJm1LIwYbC577WpDgOZjA8lPcnrmDqdl7ufn4HhI8B3PUp5aRj96uJnFUlF6uIwnkUeDTll5PH5o5HuUKtK/eemvGYhILEnG0cOTVuFcwFAp8iqzwLLZi5DHT/+ibNEaKDJn85HyIr1yX8eBNv6JSVndLMmnbRhLmac97P2VSWOFI7YzxLNiaypT0ndxyfA+Or/+Z499/hPO8B3G8+WaEhYVZ41IUpXM6kkDsgR+EECXAp8BaKWWBecNSOqNtQgk1hPJN5jesy1hHk7EcfQN4lFjiWWKFR4klw9MdEemC7L0fccaxgd0uHzF3yiIcA/1Irjho9oQSFWBg6YPnupVPHupJXNZI9L5/wuvATopXrOD0H5+l8NVXcZ47F6d77kFnUNVaitKbdGQcyPPA80KIYZjmofpJCJF3sQZupXc4m0xmhMwgsSCRhNMJ7LbYTY5nLQAWjRq8Sm1wK9HjUWzJ4Exbkt/5ECOSMvtGEp0/4baJDzJ21HSy5CmSCpNaq8rM4bweXKF34njHHVTv2k3JihUUvf4GZ959D8dZszg9+Tb2NNmrxnZF6QU6UgI5qxDIB4oBte7pNeJsIpkfMZ/XEl9jc85mIlwjGGAYgKOFI68kvMJ+YwXWRgtusZrAsdR4XMssCDhpxaGVn3No5efUWjZTYKhnq/MqYkNmMqFpPKklB83aZVgIgd34cdiNH0d9RgYlH66k9Kt1WH/2GfYeg3glbCK//9ODajyJovSgjgwkfBy4G3AD1gIL+uKYi/7g6einf7baYagh9LyG+AVN22k0VmIh9Lw69AVSDuwgK2UP7qWWBObrqDi8lf9+v5PTDlXkO9WxznU5z9/5f0T5jSa5MNksScVywAC8XnyBb0bN4sTyVdyStZu/b19MVeb3lP/6MRxuvhmh77vzbilKb9WREogf8JSUMtncwSjdr31DfNueXZHukdj7erFC/yONxlIc6i2Z03QrjdVF6NOOMizDAZEh2Lb3Rfb5eROvO8Zpp1pWuC7hrRlLWkfTd1VSiRoezGvhU/kidBI3nE7miaJ4Tv3hjxS+8QYuDz2M0513oLGxudpHoihKB3WkDeTZ7ghE6R3aJ5T23YXLDpfhNMSJBT8sgPomPMtteMjxNnKPHCT0hA2DjLYAbE16mSNhg/iyZjOnnKpZbKe96sGN50+Xcj1D/J2o+uknipe+T8FLL3HmnXcw3H8/hrlzVIO7onSDK2kDUfqptkll2+FtP0sqZ0saj363AIdSgVepNTdowzi17wAxdY6AI9VWTWzf9y4FTvWscfuARyb9hl2nd1FYW8jsAbM7vIRv++lS7GNjsY+NpWbfforff58zb71F8Qcf4HTXnbg89BB6b29zPBJFUVAJROmkC5VUFk8/P6nsz9/H775chPMZDR6llriVWBB4ylTFlLFrBcJQh3Cu571j/yKnNBtHW6dOl05sRo7A5p23qU9Pp/iDZZR+sprST1bjeMstuMz/JZahoV310RVFaaESiNJl2ieVEZ4j+c+d77YObnwl/hUsq424l1jiWqLHo8QKvyJTQmmO38oxx3rTeJTJjzN+9C1Y2thecQyWoaF4v/xP3H79K0o+/JDSz9ZS/vXX2E2ahMuC+diMHNlln1dR+ruLJhAhRCXQdrUp0bJ9dkVCBzPHpvQB7Qc3nk0mL+19iSZZgmW9qXTiUWqFe0nLeJSMlRx47yNcAwLxHTQUn0FD8Bk0FDuDc4cb5fXe3ng8+ywujz1G6SefULrqI7LnzMV65EhcFszHbuJEhKYjKzorinIxF00gUkr77gxE6fvaJ5PlB5dTWFvIqJGjWH10NY3GCqyMel4K+SNW+XWcPHqY1K0/sH/jNwBYuzpzyPok+U61fOS+lNdvW3zZ6i6dwYDbokW4PPwwZV98SfHyZeQtfBzL0AFUzLqP3YEjGRPqoQYlKkondKgKSwgxHJjQsrldSplivpCU/iDSPZI3bnijdfsG/xsuWLJobmqi8EQmJ48eZlf8BryPWxCcaw3AlsT/R37kdfgNHYbf0AgcXC8+vlVjY4PzA/djuPceKr77jrx3FmP5nxcZau3EJ2Gx8PyviApTDe6KciU6MpDwSWAB8GXLro+FEEuklG+aNTKlX2nffnKWVqfDa0AYXgPC0I0OYsH3C7CpkHiX2DBdN5TMfQkc+mkzAI4envgNiWhNKPbOrj97P6HX4zhzJh/ZDWb7qnXccWwr85PX0TBnM0UPPYjh/rmqC7CidFBHSiC/BMZIKasBhBD/AvYAKoEo3SrSPZKlU8/v6SWNRs7kZpN7KIXcw6lkxO/h4NYfATB4eePbklD8hw7D1ulcYogJceVN33ASPYcQXpbNP2r3c+bttyletgynO+/E5aF56H18euqjKso1oSMJRADNbbabW/YpSrdrX1IRGg1uAUG4BQQx8uZZGI3NFGWfaE0oabt3kLr5ewCcvX1bSyeDh0Tw/F12/JC1mynB1xEx7PFzXYBXr6Z09Wocb70Vl0cWYBkc3FMfV1F6tY4kkOXAXiHEVy3btwEfmC8kRek8jUaLR1AIHkEhRN96O0ZjM4XHs0wJ5VAKh3ds5cCPplWUy+yb0DjXsib9G4LtX2FU6HWtXYCLl6+gbK2pC7D95Mm4PPII1uFDe/jTKUrvcskEIoTQAHHANmB8y+6HpZT7zRyXonQJjUaLZ0goniGhjJp5h6lR/ngmX2z+gJrUJELybBmUrWF74ksc8Q9sKaEMw/fJX+O68DFKVq6k9ONPqPzhB2zHj8f1sUexiTbftPaKci25ZAKRUhqFEG9LKUcA+7opJkUxG61Oh1doGOMd5/Ch3Vaamxpxr7DhUcPdNJ0oImXTRvZ9tx6EwD0gGL+hEfi+/m8cDhyi8qOPyb7/AayjorAYG4OcOFEtu6v0ax2pwtoshLgD+FJKKS979hUQQkwD3gC0wPtSypfbHbcEVgJRmNYhuUdKeaLl2LOYGvibgV9LKb/vytiUvu1C83kBNDU2kp+RRu6hVHIPpZD8w7ckNTai0WrxmBSDu9BhG5+Ew9tvc3zzFlwffQT7yZMRWm0PfyJF6X4dXRP9aaBJCFFHF41EF0JogbeByUAekCCEWN9urZFfAqVSygFCiHuBfwH3CCGGAPcCQwFvYJMQYqCUshlF6aALdR3W6fX4Dg7Hd3A4Y++8j8aGek6lHSH3UAo5qQdIyTiCNFihcQ7FpaEe55eex/PNNwh+6JcYZs1U67cr/UpHpnM314j00UCGlDILQAixBpgFtE0gs4C/t/z+OfCWMNUZzALWSCnrgeNCiIyW99tjpliVfkpvYUlARCQBEZFwL9TXVJN35CC7v/uW5rJi0vSCNGDXpx/g+vH7BIwcTdjcB3ANHqCqt5Q+ryMDCTdLKW+83L5O8AFy22znAWMudo6UskkIUQ64tOyPa3ftBTvtCyEeAR4B8PDwYNu2bZ0KtqqqqtPX9iXqOZgYIkdjZ2eHT001lSdzqD54gPJTeew9tI+9f9qHXqvDPiCIem8H8pyrCXEdSpBlUE+H3eXU34NJf30Ol5pM0QqwAVyFEAbOjf1w4CJf1r2RlHIJsAQgOjpaxsbGdup9tm3bRmev7UvUczC52HPI3/QjaatWcPJ0HsX19TRkaTEA2TYHsYqMYUT0JPzDh2Pj6NTtMZuD+nsw6a/P4VIlkEeBpzC1MSRxLoFUAG91wb1PYlou9yzfln0XOidPCKEDHDE1pnfkWkXpdp43TcbzpsnUHjrE9n/8AcfDORQ52pDmbc3JhP2c2p0IgKt/IP5Dh+EfMRzfwRFYqqV4lWvQpWbjfQN4QwjxKzPNe5UAhAohgjB9+d8LzGl3znpgHqa2jTuBLVJKKYRYD3wihHgNU4ILBeLNEKOidIr10KHsnPcYew/+hdvjKph2oByj9gyaW2+mYshATuYcb+0yrNFq8RwQZmprGTYCrwED0aheXco1oCON6G8KIa4DAtueL6VceTU3bmnTeAL4HlM33mVSykNCiBeARCnlekwj3le1NJKXYEoytJz3GaYG9yZgkeqBpfQ2tw0ex+dJj/H6mAzWDHDjzYZStN9uwPDtdwTNno3ji69xprqC7NT9ZKcks+eL1ez5/BMsrG3wDx9GQMQIAoZF4uTprRrklV6pI43oq4AQIJlzc2JJTOMzroqUcgOwod2+59r8XgdccLFsKeX/A/7f1cagKOYSFWDgo/vvIS6rmJhgF4YFGGg89WuK33+fsrWfU/bFFzjOmsmYRx9lwn3zqK2sIOdgSmtCyUgw9RNxcHNvLZ34hw/H2l6t5ab0Dh0ZBxINDOnqQYSK0h9EBRjOW6xK7+2N53PP4fLooxS//wFln31G+bqvTRM3PvooYWPHEzZ2PFJKygpOk52STHbKfo7F7SJ1yw8gBB5BIa0JxTtsCDq9vgc/odKfdSSBHAQ8gdNmjkVR+g29hweef/4Tro8soHjZckrXrKH8m29wmD4d14WPYTlgAAZPbwye3kROuRljczP5memtpZPE/31F/Nefo7OwxHfwUFNCGT4SV78ADhQd6NCyv4pytTqSQFyBw0KIeKD+7E4p5UyzRaUo/YTOzQ2PP/wel/m/pGTFCko+/oSKDRuwnzoV14WPYRUWBoBGq8V74CC8Bw5i7B330VBbQ+7hVFMJJTWZnz5aBh8tw8LBjmP2ReS51PChx1LenHn5ZX8VpbM6kkD+bu4gFKW/07m44P7MMzj/4hemGYBXfUTlxo3Y3XQjrgsXYj30/KnkLaxtCIkaQ0iUaextZfEZslP2s2XHF3iklxN40rTs79bkf1E66gYChkXiNzgCvZVVt382pe/qSC+sn4QQHsColl3xUspC84alKP2TzmDA/ckncXnoIUpWfUTJypWc2LQZu9hYXB9fiPWwYRe8zt7FlfBJk2ka6saC7xdgVwa+xTbcIENI+fE79m34Go1Wh3fYoNbeXR7BA9BoVHdhpfM60gvrbuDfmNYEEcCbQojfSSk/N3NsitJvaR0dcXtiEc7zHqT0448pWb6CE3ffY1qTZNHj2IwYccHrLrTsb1NDAyePHm5tP9n16Sp2fboKK1s7/Fq7C4/AycOzmz+lcq3rSBXWn4FRZ0sdQgg3YBOmyQ0VRTEjrb09ro89huH+Byhbs5riZcvJvm8OtteNxXXRImyion52TftZhnUWFgQMiyRgWCTMhZqKcnJSTW0n2SnJpO/dDYCjhycB4abz/MKHY21nrnlUlb6iIwlE067KqhjQmCkeRVEuQGtni8v8+RjmzKF09RqKly0je+792MTE4LbocWxGjbr8m7SwcXBk0LiJDBo3ESklpadPtiaTo7u3k7J5IwiBZ/AAAoaNICAiEq+Bg1V3YeVnOpJANgohvgdWt2zfA3xnvpAURbkYjY0NLr/8BYY591G65lOKP/iA7AcexGb0aFwXLcJ2zOgrej8hBM7evjh7+zJi6q0Ym5s5nXGM7JT9ZKcmE//15+z96jN0lpb4Dg5vHX/i6hegRscrHWpE/50QYjbn1kRfIqX8yrxhKYpyKRpra1wefgjDvfdQ9tlnnHn/fXLmzcMmOpqSux5kt30gMSGu5w1i7ND7arX4hA3GJ2ww1901h/qaGvKOpLYOaPxp1QcA2DoZ8I+IpFZnSdWwCOycXczxMZVe7lLTuQ8APKSUu6SUXwJftuwfL4QIkVJmdleQiqJcmMbaGud583C65x7KPlvL6cVLsPrDr3F3CeJfQ6fyhz89QFSgc6ff39LmAt2FU03JJDtlPzXlZSzeuhEXX38CIiJp9nfkuGMpo/xi1PiTfuBSJZDXgWcvsL+85dgMs0SkKMoV01hZ4fzgA3zqGcXhDz7i7rQtPL/9PcryfqLqT89gO+66LqlysndxJTz2JsJjb0IajXz35ee4WerITtlP8qYNGBubMArJasOXHB87nTEx0/EMCVWzC/dRl0ogHlLK1PY7pZSpQohAs0WkKEqnjQ7z4o2BE/ghYAzT8hJ4JHcnufPnYz18OK5PLMJ2/Pgua7sQGg02ru6Mio1l1IzZLN2/mM+3fIDnGUt8iq3J+nYzWd9uxtLGFr+hw1raT9Tswn3JpRLIpZZMs+7qQBRFuXpRAQY+nh/TMgPw9Qz2eo7yL7/izJLF5C54BKthw3Bb9Di211/f5V/io3zGsMTjfQrcKjmsqePt617HscDY2iCfkbAHOH92Yb+hw7BxcOzSOJTuc6kEkiiEWCClXNp2pxBiPqYVChVF6YXazwBsuPcenGbfTtm6dRQvXkLuo49hFR6O66LHsYuN7bJEEukeydIp5w9iJBjCxk645OzC7oHBrd2FfcKGoLOw6JJ4FPO7VAJ5CvhKCDGXcwkjGrAAbjd3YIqidB1hYYHh7rtxuv12yr/+mjPvLSZv4eNYDRliSiQ33NAliaT9IMbW+wtxwdmFzw5oTPrfVyR8/Tk6vQU+Z2cXHjYCN/9AhEYNO+utLrWkbQFwnRBiEhDesvtbKeWWbolMUZQuJ/R6nO68E8dZsyhf/w1nFi8mb9ETWA4ejOvjC7G/8cZu+cJuO7twzB330lBXS97hg63VXds/Xg4fL8fawdHUfhI+HL/wYTh5eKn2k16kI+NAtgJbuyEWRVG6idDrcbpjNo6zZlL+zf848967nPzVr7EMC8P18cexn3xTt/7P38LKmuCRowgeaRpRX1VSTHZqMjmpyeQcPMCxPTsAsHd1w3/ocPwjhuM/dJgaf9LDOjISXVGUPkrodDjdfhuOM26l4ttvOfPue5x88kksBw40lUimTOmRKiQ7ZxeGTryRoRNvbJ1uJedgCjkHk8lM2suhnzYB4Ozti1/4cALCh1PmDgeqDquFtLqRSiCKoiB0OhxnzcLh1lup2LCBM++8y8mnfoNl6ABcFy7EfupURA+N5Wg73UrklJuRRiOF2cfJPXiAnIMHOPzTZg788C0SSYlDI3FuH3PflIVcHzNTrX9iZj2SQIQQzsCnQCBwArhbSlna7pxI4F3AAWgG/p+U8tOWYyuAiZgGNQI8JKVM7o7YFaUvE1otjjNm4HDzzVR8t5Ez777LyaefwSLkHVwXLsRh+rQeSyStMWo0eASF4BEUQvSM2TQ3NbH0h1fZtft/eBZbMvC4DcnvriRlySd4hQ7EP3w4/kOH4zUwDK1OTQjZlXqqBPJHYLOU8mUhxB9btv/Q7pwa4EEpZboQwhtIEkJ8L6Usazmu1iRRFDMRWi2Ot96Cw/Rp/7+9e4+rqkwXOP57gA1bEJCbaKKEoqQpYaAiHSfF0bx0zNNY081Lps500aYZa+g0nWk6TWPljFliJy+Fml29NmbNmIhpXgpT8YKmGSooKHgDDQV5zx97Q3gB5bqB/Xw/n/3Za6291trPel3y7PWutZ5F/r/+Re5bb3Fk8mRyZ84k8NHf4jNkiMMTSSlXNzdiY+5g7olF7CjJx1pi4a/t/4gls4DDO7ezafFHbFz0AW4eHrSJ6EK7rrcQ2i2KoBvD9IFaNeSoBHIX0Nc+PA/bw6ouSSDGmO/LDR8RkWNAEHAKpVS9EFdXfIYMwXvQIPL//W9yE2dy5OlnyE20JRKaN3d0iEAF96DYFZ4tIHP3Tg7Zu7zWvZ/EOsDq1ZyQLt1o1zWSkC7dCAxpp5cMV5GjEkiwMeaofTgbCK5sZhHpie3+k/IFHP8qIv8DrAYSjDHn6yRSpRTi4oLPoEF4DxxI/pdfkjvzLY78MYGAoCBOncnH9z/vRNwce0q1ontQrF7NCe8RS3iPWADOnjrJoV1pHNphSyild8hbvX0Iuelm2nbpSkiXbnoPynUQY0zdrFjkS+Bqz8h8DphnjGlRbt6Txpir1p0WkdbYjlBGG2M2lZuWy89hBQAAGz5JREFUjS2pzAJ+MMa8WMHyE4AJAMHBwdEffvhhtbanoKCA5g3k15YjaTvYOH07lJTgkZZGs0//iceRIxQHBXF28CB2dYoh/TTc5O9KuF/j6B46f+Y0BUcOk380k/ysw1zIt51adXX3oHnrELzbtKX5DSF4BrSsMKE09f2hX79+W4wxMZdPr7MEUhkR2Qv0NcYcLU0QxpiIq8zngy15vFzR+Q4R6QtMNsbcea3vjYmJMampqdWKOSUlhb59+1Zr2aZE28FG28EmZc0aoo0hd0Yihbt3k+3lz0ed+rOufU/mT7itys8jaQjO5B4nM30nmbt3cHj3Dk5l2zpL3Jt50uamLoR07spPra1873GEHq17EtUyqsnvDyJy1QTiqGPOT4HRwBT7+/LLZxARd2ApMP/y5CEire3JR4DhwM66D1kpdQURvPv1o3m/fnz05gc0e/9dntz6Cfft/ZJDHg9wa8IEpJHVtvIJDKJLn3506dMPsN3UeLgsoezkx622H6FFriW85/8JGbF34uYWyMXiIqe7ystRCWQK8LGIPAIcBO4FEJEY4LfGmHH2ab8AAkRkjH250st1F4pIECDANuC39Ry/UqocEaHjXYN58FgLbsnazQN7V3HT+4nsT1lC4IQJ+N59Ny6NLJGUau4fQOfbbqfzbbcDMGvjDFasfY+WeR60PmHlh3+uAuCHFYtpFd6JGyI60+amLtzQsTPWJtytBQ5KIMaYPKD/VaanAuPsw+8B71WwfHydBqiUqrLoUD8Wju/NpgOdCAkbS9vDu8idkUj2C38h9+1ZBIwfR4sRIxptIinVs8N/MPuHeRy64TQWl3Mk9n6drJRdtHATjuzdzbefLuabZZ8AEBDSjjYRXWwJJaILvi2Dm1QtL4ecA3GUq50DKSoqIjMzk8LCwkqXLSwsxKp3tV7SDlarlZCQECwW5zpsBz0HUupa7WCM4ezXG8hNTOSnrVtxCw4mYPx4WtwzAhcPj/oLtJZtO7btkkuGy7dDUWEhR/d/z5G9u8n6Pp0je9O58NM5wPYs+RsiOtMm4mbaRHQm6Mb2uDr46rXr0dDOgTQYmZmZeHt7c+ONN1b6yyA/Px9vb+96jKxhKm0HYwx5eXlkZmYSFhbm6LBUAyUiNP+P2/C6LY5zGzdyPHEmOS+9RN6sWQSMG0eLe+/BpRH+MKvokmEAi9VKu66RtOsaCUBJyUXyDh8ia2+6Lans3c2+zRsAcPPwoHWHTrTqGEHr8E6c8CtmR+H3jaael9MnkMLCwmsmD3UlESEgIIDjx487OhTVCIgIXnFxePbuzbnNm8mdkUjOyy+TO3sWgePG0eLXv26UieR6uLi4EhQaRlBoGFEDhwCQfyKXI3vTydq7myN797BlxTJKLhYDcNZazLt+HxLfazjR3fsRHNYBi0fDbBunTyCAJo9q0nZTVSUieMXG4hUby9nN35CbmEjO36aQO3sOAY88gt+v78XF09PRYdY5b/9AInr3IaJ3HwCKL1xgzurXSdm8nMBT7gSd8mD/si/Yv+wLxMWFoHZhtO4YQeuOEbQK74R/6zYN4iZHTSBKKYfw6tUTr149OffttxxPnMmxV14hb84cAsY+jN/99ztFIinl5u5ObPRA5uZ9wp6Ss1hcLMyInUbAKQvZ+/dydN9e0tevYfuqlQB4eHnRqkMnW0Lp0JFWHTqx78LBq5ZyqdO46+VbVJ3LyMhgw4YNPPDAA1VaLikpidTUVGbMmHHJdGMMTz75JCtXrsTT05OkpCRuvfXW2gxZKQA8e/QgNOldzm3ZQm5iIsdem0renLn4j30Y/wcewMXLy9Eh1ouK6nmFx/QCbOdSTmRlcnT/XrL3fc/RfXvYvORjjCkB4Jz1Isd9z/OV3wLG9J9IXPQgrF7NrzjhX5s0gVTDloMn2XQgj9j2AQ3mTtuMjAzef//9qyaQ4uJi3Kp4pcfnn3/Ovn372LdvH5s3b+bRRx9l8+bNtRWuUlfwjI6m3TvvcO67reQmJnL87//gxNx38H/4YfwefBDX5k0/kVR2ct7FxZXAtqEEtg2lW7+BgO2Kr5yMH1i+fgE/7txMwGl3QnM82bonia0k0SwogD3umRzzKSSp7Sxm/OesWk0imkCqaMvBkzw4ZxMXiktwd3Nh4bjYGieR+fPnM3XqVESEyMhIFixYQEZGBmPHjiU3N5egoCDeffdd2rVrx5gxY/Dx8SE1NZXs7GxeffVVRowYQUJCAunp6URFRTF69Gj8/PxYsmQJBQUFXLx4kaVLlzJ27FgOHDiAp6cns2bNIjIyssKYli9fzqhRoxARYmNjOXXqFEePHm3S9X5Uw+B5a3fazZ3DT9u2cTxxJsenTePEO+/g//AY/B56CFfdB8tYrFZCbrqZ2/zvI8m6mqKSIjyL3fnf8GdonlfCpq2rCDxoITSrGZmtjpKak6oJxJE2HcjjQnEJJQaKikvYdCCvRglk165dvPTSS2zYsIHAwEBOnDgBwMSJExk9ejSjR4/mnXfeYdKkSSxbtgyAo0ePsn79evbs2cOwYcMYMWIEU6ZMYerUqaxYsQKwdU199913pKWl4e/vz8SJE+nevTvLli0jOTmZUaNGsW1bxc/gysrKom3btmXjISEhZGVlERFxRckypepEs6go2s2exU9paRxPTOT469PJezcJ/9Gj8B85Ele9rL5MRd1fHrd1Yvy/x+N27iLFzVyJCb7iVo4a0QRSRbHtA3B3c6GouASLmwux7QNqtL7k5GTuueceAgMDAfD39wdg48aNLFmyBICRI0fyzDPPlC0zfPhwXFxc6NKlCzk5ORWue8CAAWXrW79+PYsXLwYgPj6evLw8zpw5U6PYlaoPzSIjaff22/y0Ywe5iTPJfeNNTiTNw3/UKPxHjcTVx8fRITYIV+v+quw5KbVBE0gVRYf6sXBcrEPPgXiUu4O3skoCXjU4+dimTRsOHz5cNp6ZmUmbNm2qvT6laqpZt260/b+3+GnXLnJnvkXujBmcmDcP/5EP4TdyJNvP0ODOTTYElZ1XqSnHX0jcCEWH+vF4v/Ba2Unj4+P55JNPyMvLAyjrwoqLi6P02SULFy6kT58+la7H29ub/Pz8Cj/v06cPCxcuBGzlJwIDA/Gp5JfbsGHDmD9/PsYYNm3ahK+vL61bt67StilVF5rdfDNtE2cQtnQJXrG9yJ35Ft/3jefLJ54lafk3PDhnE1sOnnR0mE5Bj0Ac7Oabb+a5557j9ttvx9XVle7du5OUlMSbb77Jww8/zGuvvVZ2Er0ykZGRuLq6cssttzBmzBj8/C5Nbi+88AJjx44lMjIST09P5s2bV+n6hgwZwsqVKwkPD8fT0/Oa369UfbN27kzIm29yft8+vn5pGkO/Wcvg/etZHRrD9q5eRIf+wtEhNnlOX0wxPT2dzp07X3NZrYVlc3k7XG/7NTVaTNGmobTDloMneWr6ZwxLT2bAwW+xUILv4MEETBiPtR4u/Ggo7VBXtJiiUqrJig71Y9qTQ9l0IJZiX0Or1cs49cGHnPnsM5r360fgbybQLKrhFydsbDSBKKWahOhQv5/PS976NIETJnBi4UJOzptPxn3349mrFwETxuMVF6d13GqJnkRXSjVJrr6+BD32GOHJq2mZ8Ecu/Pgjhx8ZR8Y993Jm1SpMSYmjQ2z0NIEopZo0Fy8vAsaMocOXq2j14l+4eOYMWRMncWDYME4vX44pLnZ0iI2WJhCllFNwcXfH79576bDyM26YOhURF478MYEf7hjEyQ8+oOQaTyVVV9IEopRyKuLmhu+dQwlbvoyQmTNxCwwk+y8vsj++P8dnzqT4pN5Dcr00gTQRpdV4qyopKYknnnjiiul79uyhd+/eeHh4MHXq1NoIUakGRVxc8I7vR+iHH9Bu/jyadetG7htvsj++P9n/+xIXMjMdHWKD55AEIiL+IrJKRPbZ3696S7eIXBSRbfbXp+Wmh4nIZhHZLyIfiYh7/UUPHP4G1v3d9t5AVJZAiqvRx+vv788bb7zB5MmTaxqaUg2aiODVsydt3/4/2v/zU3wGDeLkxx/zw8A7yPr97/lp5y5Hh9hgOeoIJAFYbYzpCKy2j1/NT8aYKPtrWLnprwDTjDHhwEngkboNt5zD38C8YZD8V9t7LSSR+fPnExkZyS233MLIkSMBW0KIj48nMjKS/v37c+jQIQDGjBnDpEmTiIuLo3379ixatAiAhIQE1q1bR1RUFNOmTSMpKYlhw4YRHx9P//79OXHiBMOHDycyMpLY2FjS0tIqjally5b06NEDi8VS4+1TqrHw6NiRG/72MuFfriJg7MMUfLWOjBEjODh6DAVffVVp7Tln5KgEchdQWktjHjD8ehcU2wXc8cCi6ixfYxnr4OIFMBdt7xnrarS60nLuycnJbN++nenTpwM/l3NPS0vjwQcfZNKkSWXLlJZzX7FiBQkJttw7ZcoU+vTpw7Zt23jqqacA+O6771i0aBFr167lz3/+M927dyctLY2XX36ZUaNG1ShupZoyS3AwLSdPJnxNMi2ffpoLGRkcnvAbfhx2F6eWLWPL/mMkrtnv9DW3HHUjYbAx5qh9OBsIrmA+q4ikAsXAFGPMMiAAOGWMKe2XyQQqLBMrIhOACQDBwcGkpKRc8rmvr2+lRQhLXbx4kfz8fFxaRuPpaoGLgKuFcy2jKbmO5SuycuVK7rrrLjw8PMjPz8disZCfn8+GDRuYN28e+fn5DB8+nKeffpr8/HyKioq44447OHv2LG3btiUnJ4f8/HzOnTtHcXFx2bYUFhbSt2/fsvV99dVXLFiwgPz8fHr06EFubi5ZWVkUFhZy4cKFCtvg/PnzZeso3w6lCgsLr2hTZ1BQUOCU2305p2iHDu3h+T9h/fZbPFd9yfmEZzln9eVAhz7Mbt+LJ+L8aGX5qem3w1XUWQIRkS+BVlf56LnyI8YYIyIVHReGGmOyRKQ9kCwiO4DTVYnDGDMLmAW2WliX16tJT0+/rhpXZTWgIvrC6H/ajjxu7INX255VCecKVqsVd3f3K2IQEby9vbFYLBQVFV0y3qJFi7L5jTF4e3vj6emJm5tb2XSr1XrJfC4uLjRv3rxsvHR9FX1/KQ8PDzw8PMo+v7wWltVqpXv37jVqg8aoqdc+ul5O1Q6//CUmIYEPZ3yEy0fvMW7XCh7Yu4q8goE06xdDnzvvdHSE9a7OurCMMb80xnS9yms5kCMirQHs78cqWEeW/f0AkAJ0B/KAFiJSmvxCgKy62o6ratsT+vzB9l5DDbWcu1LqSiJCp2F38Od+j/Nk39+xuU03QtZ9TsCfX+DwY49zdtNmpzpP4qgurE+B0cAU+/vyy2ewX5l1zhhzXkQCgduAV+1HLGuAEcCHFS3fWDTUcu7Z2dnExMRw5swZXFxceP3119m9e7fWEFJO7+eHynWkV/t76disiG1TXsFt40YOJSfjERGB/6hR+Nw5FJdyD39rihxSzl1EAoCPgXbAQeBeY8wJEYkBfmuMGSciccDbQAm2I6XXjTFz7cu3x5Y8/IGtwEPGmPPX+l4t515zWs7dxqm6biqh7WCTkpLCL2JjObNiBSfmzef8vn24+vvjd999+N1/H25BQY4OsUYaVDl3Y0we0P8q01OBcfbhDUC3CpY/ANS8/0gppWqJi9VKixEj8P3Vrzi3aRMn5s0nd+ZM8mbPxmfIEPweeohm3bo6OsxapeXclVKqFokIXr1749W7N+d//JGT7y3k1NKlnF6+HGtkJH7334/PkMFNontLS5kopVQd8QgLo9Xzf6Lj2hSC//QnSgoKOPrss+y/vS85r73W6MulaAJRSqk65urtjf9DD9L+sxW0S0rCs2dPTiTN44cBAzn0m99QsHZto3w+iXZhKaVUPRERvGJ74RXbi6KcHE599DEnP/mYw7/5LZaQEM4OuotNEXHERIb9/HTFBkyPQJRSygEswcEETZpIx+Rk2kz7B4X+QbjPSST2j6P5dvwTbP1ncoO/p0QTSBNR2+XcFy5cSGRkJN26dSMuLo7t27fXRphKqcuIxYLP4MGsffRFHu0/mS9CexJ9dDfWpx/nwOAh5M2dS7H9RuOGRhNINWw7to05O+aw7dg2R4dSprbLuYeFhbF27Vp27NjB888/z4QJE2oaolKqErHtA8j2v4FZUXfzyJ0vcH7y87gGBHDstans69uPzN89RcHXXzeocyV6DqSKth3bxvh/j+fCxQu4u7oze+BsolpG1Wid8+fPZ+rUqYgIkZGRLFiwgIyMDMaOHUtubm7Znejt2rVjzJgx+Pj4kJqaSnZ2Nq+++iojRowgISGB9PR0oqKiGD16NH5+fixZsoSCggIuXrzI0qVLGTt2LAcOHMDT05NZs2YRGRlZYUxxcXFlw7GxsWQ28qtFlGrofr7DPY/Y9gFEhfrBuAc4v38/pz5ZxOnly8n/4gssISG0GPErfP/rbizBLR0asx6BVFFqTioXLl6ghBKKSopIzUm99kKVaAzl3OfOncvgwYNrtJ1KqWuLDvXj8X7hl5xA9wgPJ/jZBMK/WssNf5+KJSSE469PZ398PIcfe5z8NWswxcVsOXiy3kvM6xFIFcUEx+Du6k5RSREWFwsxwVfc3V8lycnJ3HPPPQQGBgK2JwECbNy4kSVLlgAwcuRInnnmmbJlhg8fjouLC126dCEnJ6fCdQ8YMKBsfevXr2fx4sWArYBjXl4eZ86cuWZ8a9asYe7cuaxfv756G6iUqhUu7u74Dh2K79ChXDh4kFOLFnNq6VIKkpMp8Q/k85aRfBHSgzf9glk4LrZeruLSBFJFUS2jmD1wNqk5qcQEx9S4+6o6PMrdwVrZVRpeXl41+p60tDTGjRvH559/TkBAQI3WpZSqPe6hobT8w+8JmjSR/JQUts9awPCdKfxqTzK7AsI4aD1I99+NxKWGfwOuRbuwqiGqZRTjuo2rleTRUMu5Hzp0iLvvvpsFCxbQqVOnKm2TUqp+iMWCz4ABeE6dzvih/8O7Nw/F98JZOi+Yzvd9fsGR/36OrZ+lkJi8r066tvQIxMEaajn3F198kby8PB577DEA3NzcuLySsVKqYYgO9WPGxIFsOhCNf9h/E3riIKeXLuHkipVYlywhwiuQZ297mL/94b9qtWvLIeXcHUXLudeclnO30TLmNtoONg21Hd76YiffzV9En8xtTOk1iseHRPJ4v/Aqr6dBlXNXSilV93p2bsP0Dr1IDu2Bxc2F2Pa1ey5TE4hSSjVRl99bUttXZmkCwXYlkz6qteqcqftTqcYqOtSvzi7pdfqrsKxWK3l5efrHsIqMMeTl5WG1Wh0dilLKQZz+CCQkJITMzEyOHz9e6XyFhYX6x5JL28FqtRISEuLgiJRSjuL0CcRisRAWFnbN+VJSUujevXs9RNSwaTsopUo5fReWUkqp6tEEopRSqlo0gSillKoWp7oTXUSOAweruXggkFuL4TRW2g422g422g42Tb0dQo0xQZdPdKoEUhMiknq1W/mdjbaDjbaDjbaDjbO2g3ZhKaWUqhZNIEoppapFE8j1m+XoABoIbQcbbQcbbQcbp2wHPQeilFKqWvQIRCmlVLVoAlFKKVUtmkCug4gMEpG9IrJfRBIcHU9dEZG2IrJGRHaLyC4RedI+3V9EVonIPvu7n326iMgb9nZJE5FbHbsFtUtEXEVkq4issI+Hichm+/Z+JCLu9uke9vH99s9vdGTctUlEWojIIhHZIyLpItLbGfcHEXnK/n9ip4h8ICJWZ9wfLqcJ5BpExBVIBAYDXYD7RaSLY6OqM8XAH4wxXYBY4HH7tiYAq40xHYHV9nGwtUlH+2sC8Fb9h1ynngTSy42/AkwzxoQDJ4FH7NMfAU7ap0+zz9dUTAe+MMbcBNyCrT2can8QkTbAJCDGGNMVcAXuwzn3h0sZY/RVyQvoDfyr3PizwLOOjquetn05MADYC7S2T2sN7LUPvw3cX27+svka+wsIwfbHMR5YAQi2O43dLt8vgH8Bve3Dbvb5xNHbUAtt4Av8ePm2ONv+ALQBDgP+9n/fFcAdzrY/XO2lRyDXVrrzlMq0T2vS7Ifd3YHNQLAx5qj9o2wg2D7clNvmdeAZoMQ+HgCcMsYU28fLb2tZO9g/P22fv7ELA44D79q78uaIiBdOtj8YY7KAqcAh4Ci2f98tON/+cAVNIOoKItIcWAz8zhhzpvxnxvazqklf+y0idwLHjDFbHB2Lg7kBtwJvGWO6A2f5ubsKcJr9wQ+4C1tCvQHwAgY5NKgGQhPItWUBbcuNh9inNUkiYsGWPBYaY5bYJ+eISGv7562BY/bpTbVtbgOGiUgG8CG2bqzpQAsRKX0IW/ltLWsH++e+QF59BlxHMoFMY8xm+/gibAnF2faHXwI/GmOOG2OKgCXY9hFn2x+uoAnk2r4FOtqvuHDHdvLsUwfHVCdERIC5QLox5h/lPvoUGG0fHo3t3Ejp9FH2q29igdPlujYaLWPMs8aYEGPMjdj+vZONMQ8Ca4AR9tkub4fS9hlhn7/R/yo3xmQDh0Ukwj6pP7AbJ9sfsHVdxYqIp/3/SGk7ONX+cFWOPgnTGF7AEOB74AfgOUfHU4fb+R/YuiPSgG321xBs/bergX3Al4C/fX7BdoXaD8AObFepOHw7arlN+gIr7MPtgW+A/cAngId9utU+vt/+eXtHx12L2x8FpNr3iWWAnzPuD8BfgD3ATmAB4OGM+8PlLy1lopRSqlq0C0sppVS1aAJRSilVLZpAlFJKVYsmEKWUUtWiCUQppVS1aAJRqgZEJEBEttlf2SKSZR8uEJGZjo5Pqbqkl/EqVUtE5AWgwBgz1dGxKFUf9AhEqTogIn3LPUfkBRGZJyLrROSgiNwtIq+KyA4R+cJePgYRiRaRtSKyRUT+VVoupJLvuL3c0c9WEfGuj21TqpQmEKXqRwdsNbWGAe8Ba4wx3YCfgKH2JPImMMIYEw28A/z1GuucDDxujIkC+tjXpVS9cbv2LEqpWvC5MaZIRHZgeyDRF/bpO4AbgQigK7DKVm4JV2ylwyvzNfAPEVkILDHGZNZF4EpVRBOIUvXjPIAxpkREiszPJx9LsP0/FGCXMab39a7QGDNFRD7DVq/saxG5wxizp7YDV6oi2oWlVMOwFwgSkd5gK6svIjfbh58QkScuX0BEOhhjdhhjXsFWNfqmeo1YOT1NIEo1AMaYC9hKf78iItuxVUKOs398E1d/nsTvRGSniKQBRcDn9RKsUnZ6Ga9SDZz9aq677UlGqQZDE4hSSqlq0S4spZRS1aIJRCmlVLVoAlFKKVUtmkCUUkpViyYQpZRS1aIJRCmlVLX8P4YY03ISLFo+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "figu, axsu = post.plot_u()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "plot mass of the vehicle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxU9b3/8dcnkw3CDiEiQQIIoQiKSC24EbVeLV1svdhqa6W2vdy6FNvrbav96a16a3+19ddevYuW1q3V4m2rveJSvFYN4oIWUFllEQHDDrIFzP75/XFOwhACTEImJzPzfj4e88icbeYzXya8c77ne84xd0dERAQgK+oCRESk81AoiIhIE4WCiIg0USiIiEgThYKIiDTJjrqAY9GvXz8vKSlp07b79u2joKCgfQtKQWqHgNohoHYIpHs7LFiwYLu7F7a0LKVDoaSkhPnz57dp2/LycsrKytq3oBSkdgioHQJqh0C6t4OZrTvcMnUfiYhIE4WCiIg0USiIiEiTlD6mICLpqba2loqKCqqqqiJ5/549e7J8+fJI3rs95efnU1xcTE5OTsLbKBREpNOpqKige/fulJSUYGYd/v579+6le/fuHf6+7cnd2bFjBxUVFQwZMiTh7dR9JCKdTlVVFX379o0kENKFmdG3b99W720lLRTMLN/M3jSzd8xsqZndFs4fYmZvmNlqM/tvM8sN5+eF06vD5SXJqk1EOj8FwrFrSxsmc0+hGjjP3U8BxgIXmdkE4E7gl+5+IrAT+Ea4/jeAneH8X4brJcXKLXuZubyaqtr6ZL2FiEhKSlooeKAynMwJHw6cB/wpnP8w8Pnw+cXhNOHy8y1JfypU7NzPc+vq+NvaD5Px8iIiKSupB5rNLAYsAE4E/hN4D9jl7nXhKhXAwPD5QOADAHevM7PdQF9ge7PXnAZMAygqKqK8vLzVddXWO9nmPPLXhdRvyGv19umksrKyTW2YbtQOgc7SDj179mTv3r2RvX99fX0k7z958mQ2b97Mj3/8YyZPnnzI8muuuYbZs2dTWFjIG2+80TT/5ptv5rHHHmP69OlMnz79oG2qqqpa9W+a1FBw93pgrJn1Av4MjGyH15wBzAAYP368t/VU9JELZ/PeR/mUlU061pJSWrqfzp8otUOgs7TD8uXLIx39E9Xoo1gsxsyZMxk/fnyLy//hH/6B7373u1x55ZUH1Xf33XfTu3dv8vLyDqk7Pz+fU089NeEaOmRIqrvvMrOXgIlALzPLDvcWioEN4WobgEFAhZllAz2BHcmqaUxhjJnvVlKxcz/Fvbsm621E5Bjd9tRSlm3c066vOer4Hvzosycddvm6deuYMmUKEyZM4LXXXuPjH/84V111FT/60Y/YunUrjz76KKeffjpvvvkm119/PVVVVXTp0oUHH3yQ0tJSli5dylVXXUVNTQ0NDQ08/vjjHH/88Xzxi1+koqKC+vp6brnlFr70pS+1qu5zzjmHtWvXHuOnP7Jkjj4qDPcQMLMuwAXAcuAlYEq42lTgyfD5rHCacPmLnsQbSI/pFwOgfMW2ZL2FiKSw1atXc8MNN/Duu+/y7rvv8vvf/55XXnmFu+66i5/85CcAjBw5krlz5/LWW29x++2388Mf/hCA++67j+uvv563336b+fPnU1xczOzZszn++ON55513WLJkCRdddFGUH++wkrmnMAB4ODyukAX8wd2fNrNlwGNm9mPgLeD+cP37gd+Z2WrgQ+CyJNbGgAKjuHcXylds44oJg5P5ViJyDI70F30yDRkyhDFjxgBw0kkncf7552NmjBkzpumv9d27dzN16lRWrVqFmVFbWwvAxIkTueOOO6ioqOCSSy5h+PDhjBkzhhtuuIEf/OAHfOYzn+Hss8+O5HMdTTJHHy1y91Pd/WR3H+3ut4fz17j76e5+ortf6u7V4fyqcPrEcPmaZNUGwfjdstJCXntvOzV1Dcl8KxFJQXl5BwahZGVlNU1nZWVRVxeMlbnllls499xzWbJkCU899VTTiWJf/vKXmTVrFl26dGHy5Mm8+OKLjBgxgoULFzJmzBhuvvlmbr/99iO+/wcffMDYsWMZO3Ys9913X5I+5aEy+jIXZSP688i89cxf+yFnnNgv6nJEJMXs3r2bgQODAZQPPfRQ0/w1a9YwdOhQpk+fzvr161m0aBEjR46kT58+XHHFFfTq1Yvf/OY3R3ztQYMG8fbbbyez/BZl9GUuJg7rS24si/KVOq4gIq33/e9/n5tuuolTTz21ae8B4A9/+AOjR49m7NixLFmyhCuvvJLFixdz+umnM3bsWG677TZuvvnmVr/f5ZdfzsSJE1mxYgXFxcXcf//9R9+olTJ6T6EgL5uPD+lN+Yqt/HDyx6IuR0Q6icGDB7NkyZKm6fi9gJKSkqZlEydOZOXKlU3LfvzjHwNw4403cuONNx70mhdeeCEXXnjhMdU1c+bMY9o+ERm9pwBBF9LKLZVs3PVR1KWISIbr06cPX/va15g1a1artvve977HI4880i73lc7oPQWAstJC7nh2OXNWbuPy00+IuhwRCbl7xl0U74knnmjTdj//+c/5+c9/fsj8tozqz/g9hRP7d+P4nvmUr9gadSkiEsrPz2fHjh1t+k9NAo33U8jPz2/Vdhm/p2BmTCrtz1PvbKSmroHc7IzPSZHIFRcXU1FRwbZt0QwCqaqqavV/pp1R453XWiPjQwGCLqSZb65n4fqdTBjaN+pyRDJeTk5Oq+4W1t7Ky8tbdb2gdKI/i4EzT+xHdpbpkhcikvEUCkC3vGzGl/TWcQURyXgKhVBZaX/e3byXzbtbdz9TEZF0olAIlZUWAjBnpfYWRCRzKRRCpUXdOa5HPnN0yQsRyWAKhZCZMWlEIXNXbaeuXldNFZHMpFCIU1ZayN6qOhau3xV1KSIikVAoxDlzeOPQVB1XEJHMpFCI0yM/h3GDe+u4gohkLIVCM5NGFLJ04x627tXQVBHJPAqFZpqGpursZhHJQAqFZkYN6EFh9zzdjU1EMpJCoZnGoamvaGiqiGQghUILykoL2f1RLe9UaGiqiGQWhUILzj6xkCxDV00VkYyjUGhBz645jDuht0JBRDKOQuEwJo0oZPGG3WyvrI66FBGRDpO0UDCzQWb2kpktM7OlZnZ9OP9WM9tgZm+Hj8lx29xkZqvNbIWZXZis2hJRVtofgJc1CklEMkgyb8dZB9zg7gvNrDuwwMyeD5f90t3vil/ZzEYBlwEnAccDfzWzEe5en8QaD+uk43vQr1su5Su2ccm41t3jVEQkVSVtT8HdN7n7wvD5XmA5MPAIm1wMPObu1e7+PrAaOD1Z9R1NVpZxzohCXl61jfoGj6oMEZEOlcw9hSZmVgKcCrwBnAlcZ2ZXAvMJ9iZ2EgTGvLjNKmghRMxsGjANoKioiPLy8jbVVFlZedRt+9fVsWt/LQ/NepFhvWJtep/OLpF2yARqh4DaIZDJ7ZD0UDCzbsDjwHfcfY+Z3Qv8K+Dhz/8HfD3R13P3GcAMgPHjx3tZWVmb6iovL+do256yr4YZi59nT8EgyspGtOl9OrtE2iETqB0CaodAJrdDUkcfmVkOQSA86u5PALj7Fnevd/cG4Ncc6CLaAAyK27w4nBeZ3gW5nDKoly55ISIZI5mjjwy4H1ju7r+Imz8gbrUvAEvC57OAy8wsz8yGAMOBN5NVX6LKRvRnUcUudmhoqohkgGTuKZwJfBU4r9nw05+Z2WIzWwScC3wXwN2XAn8AlgGzgWujGnkUb1JpIe7wyurtUZciIpJ0STum4O6vANbComePsM0dwB3JqqktTh7Ykz4FwdDUi8ceafCUiEjq0xnNR5GVZZwzvB8vr9xGg4amikiaUygkoKy0Pzv21bB4w+6oSxERSSqFQgLOHt4PM3TvZhFJewqFBPTtlsfJA3tSvmJr1KWIiCSVQiFBk0r78/YHu9i1vybqUkREkkahkKCy0kIaHF5epaGpIpK+FAoJOqW4F7265jBHN94RkTSmUEhQLMs4e3ghczQ0VUTSmEKhFcpGFLK9spplm/ZEXYqISFIoFFrhnBGFABqFJCJpS6HQCoXd8xg9sIfOVxCRtKVQaKWyEf1ZuH4Xuz+qjboUEZF2p1BopbLSQuobnFc0NFVE0pBCoZXGDupFj/xsHVcQkbSkUGil7FhW09BUdw1NFZH0olBog0mlhWzdW83yTXujLkVEpF0pFNqgrHFo6kp1IYlIelEotEH/HvmMGtCDcl3yQkTSjEKhjSaVFrJw3U72VGloqoikD4VCG5WNKKSuwXlttYamikj6UCi00bjBvemel60uJBFJKwqFNsqJZXHW8H6Ur9DQVBFJHwqFYzBpRCGb91Sxcktl1KWIiLQLhcIxmFSqq6aKSHpJWiiY2SAze8nMlpnZUjO7Ppzfx8yeN7NV4c/e4Xwzs3vMbLWZLTKzccmqrb0M6NmFkcd113EFEUkbydxTqANucPdRwATgWjMbBdwIvODuw4EXwmmATwHDw8c04N4k1tZuJpUWMn/dh1RW10VdiojIMUtaKLj7JndfGD7fCywHBgIXAw+Hqz0MfD58fjHwWw/MA3qZ2YBk1ddeJo0opLZeQ1NFJD1kd8SbmFkJcCrwBlDk7pvCRZuBovD5QOCDuM0qwnmb4uZhZtMI9iQoKiqivLy8TTVVVla2edt4dQ1Ofgx+X/4OudvePebX62jt1Q6pTu0QUDsEMrkdkh4KZtYNeBz4jrvvMbOmZe7uZtaq8ZzuPgOYATB+/HgvKytrU13l5eW0ddvmzqmYz9KNe5g0aRLxny8VtGc7pDK1Q0DtEMjkdjhq95GZXWpm3cPnN5vZE4keBDazHIJAeNTdnwhnb2nsFgp/Ng7d2QAMitu8OJzX6U0qLWTDro9YvVVDU0UktSVyTOEWd99rZmcBnwTuJ4GDwBb8yXw/sNzdfxG3aBYwNXw+FXgybv6V4SikCcDuuG6mTq2stD+A7t0sIikvkVCoD39+Gpjh7s8AuQlsdybwVeA8M3s7fEwGfgpcYGarCELmp+H6zwJrgNXAr4FrEv8Y0RrYqwvD+3fT0FQRSXmJHFPYYGa/Ai4A7jSzPBIIE3d/BThcB/v5LazvwLUJ1NMplZUW8vBr69hXXUdBXoccvxcRaXeJ7Cl8EXgOuNDddwF9gO8ltaoUNGlEf2rqG3j9vR1RlyIi0maJ/MW/n6Dff5+ZnQDkAKk39jLJPj6kN11zYzquICIp7aj9HGb2beBHwBagIZztwMlJrCvl5GXHOGNYX8pXbsXdU25oqogIJHZM4Xqg1N3VL3IUk0r789flW1mzfR/DCrtFXY6ISKslckzhA2B3sgtJB2UjGq+aqi4kEUlNiewprAHKzewZoLpxZrNzDwQY1KcrQwsLmLNyG984a0jU5YiItFoiewrrgecJzk3oHveQFpSN6M+8NTv4qKb+6CuLiHQyR91TcPfboOkaRri7ruVwBGWlhTzw6vvMW7ODc0f2j7ocEZFWSeTaR6PN7C1gKbDUzBaY2UnJLy01nT6kD/k5Wbobm4ikpES6j2YA/+Tug919MHADwWUopAX5OTEmDu2r8xVEJCUlEgoF7v5S44S7lwMFSasoDZSV9mftjv2s3b4v6lJERFolkVBYY2a3mFlJ+LiZYESSHEZZaePQVHUhiUhqSSQUvg4UAk8Q3BuhH3BVMotKdYP7FlDStyvl6kISkRSTyHkKn3T36fEzzOxS4I/JKSk9TBjal+eWbo66DBGRVklkT+GmBOdJnPycGPUNrbrTqIhI5A67p2BmnwImAwPN7J64RT2AumQXJiIiHe9I3UcbgfnA54AFcfP3At9NZlEiIhKNw4aCu78DvGNmv3f32g6sSUREIpLITXYUCCIiGSKRA80iIpIhFAoiItIkkdtxjgC+BwyOX9/dz0tiXSIiEoFETl77I3AfwUXwdJOABOXEjKq6BvbX1NE1N5FmFhGJXiLdR3Xufq+7v+nuCxofSa8sxX1qzABq6hp47M0Poi5FRCRhiYTCU2Z2jZkNMLM+jY+kV5bixp3Qm08M6cNv5q6hpq4h6nJERBKSSChMJTim8BrBSWwLCE5qOyIze8DMtprZkrh5t5rZBjN7O3xMjlt2k5mtNrMVZnZh6z9K53N12TA27q5i1jsboy5FRCQhiZynMKSFx9AEXvsh4KIW5v/S3ceGj2cBzGwUcBlwUrjNf5lZLPGP0TlNGlHIxwb04L4579Gg6yCJSAo4bCiY2Xnhz0taehzthd39ZeDDBOu4GHjM3avd/X1gNXB6gtt2WmbG1WXDWL21kr8u3xJ1OSIiR3WkYTGTgBeBz7awzAnur9AW15nZlQRdUDe4+05gIDAvbp2KcN4hzGwaMA2gqKiI8vLyNhVRWVnZ5m1bo6DBKexi/HTWW+RszcfMkv6erdFR7dDZqR0CaodARreDuyftAZQAS+Kmi4AYwR7KHcAD4fz/AK6IW+9+YMrRXv+0007ztnrppZfavG1r/e71tT74B0/76+9t77D3TFRHtkNnpnYIqB0C6d4OwHw/zP+rHXpGs7tvcfd6d28gOO+hsYtoAzAobtXicF5amHJaMf265XFv+XtRlyIickQdGgpmNiBu8gtA48ikWcBlZpZnZkOA4cCbHVlbMuXnxPj6WSXMWbmNpRt3R12OiMhhJS0UzGwm8DpQamYVZvYN4GdmttjMFgHnEt6Xwd2XAn8AlgGzgWvdPa3Onr5iwmC652Vz35w1UZciInJYRw0FM7vUzLqHz282syfMbNzRtnP3y919gLvnuHuxu9/v7l919zHufrK7f87dN8Wtf4e7D3P3Unf/y7F9rM6nR34OX5kwmGcWbWTdjn1RlyMi0qJE9hRucfe9ZnYW8EmCg8D3Jres9PT1M0vIjmUx42XtLYhI55RIKDR243wamOHuzwC5ySspffXvkc+U04r544IKtu6tirocEZFDJBIKG8zsV8CXgGfNLC/B7aQF084eSl19Aw++ujbqUkREDpHIf+5fBJ4DLnT3XUAfgmshSRuU9Ctg8pgBPPL6OvZU6U6nItK5JBIKA4Bn3H2VmZUBl5JGw0Wj8K1Jw9hbXcej89ZHXYqIyEESCYXHgXozOxGYQXCS2e+TWlWaGz2wJ+eMKOT+V96nqjatRt6KSIpLJBQa3L0OuAT4d3f/HsHegxyDqycNY3tlNY8vrIi6FBGRJomEQq2ZXQ5cCTwdzstJXkmZYcLQPowd1ItfzVlDXb1uwiMinUMioXAVMBG4w93fDy9D8bvklpX+Gi+rvf7D/fxlyeaoyxERARK7yc4yd5/u7jPD6ffd/c7kl5b+LvhYEcMKC7i3/L3Gq8OKiEQqkctcDDezP5nZMjNb0/joiOLSXVaW8a1Jw1i2aQ8vr9oedTkiIgl1Hz1IcFmLOoKL2P0WeCSZRWWSi8cOZEDPfO4tXx11KSIiCYVCF3d/ATB3X+futxJc8kLaQW52Ft88eyjz1nzIwvU7oy5HRDJcIqFQbWZZwCozu87MvgB0S3JdGeWyjw+iV9cc7tNNeEQkYomEwvVAV2A6cBrwVWBqMovKNAV52UydWML/LtvC6q17oy5HRDJYIqOP/ubule5e4e5Xufsl7j6vI4rLJFPPKKFLTkw34RGRSGUfboGZzTrShu7+ufYvJ3P1KcjlstMH8bvX1/FPF4zg+F5doi5JRDLQYUOB4IS1D4CZwBuAdUhFGeybZw/ld6+v4zdz3+dfPjsq6nJEJAMdqfvoOOCHwGjgbuACYLu7z3H3OR1RXKYZ2KsLF48dyMw317NzX03U5YhIBjpsKLh7vbvPdvepwARgNVBuZtd1WHUZ6FuThvJRbT0PvbY26lJEJAMd8UCzmeWZ2SUEJ6tdC9wD/LkjCstUw4u6c8GoIh5+fS37quuiLkdEMsxhQ8HMfgu8DowDbnP3j7v7v7r7hg6rLkNdXTaMXftreexvH0RdiohkmCPtKVwBDCc4T+E1M9sTPvaa2Z6OKS8zjTuhN58Y0offzF1DTZ0uqy0iHedIxxSy3L17+OgR9+ju7j06sshMdHXZMDbtruLJt7VjJiIdJ5EzmiUCk0YU8rEBPbhvzns0NOiy2iLSMZIWCmb2gJltNbMlcfP6mNnzZrYq/Nk7nG9mdo+ZrTazRWY2Lll1pYrGm/C8t20fzy/fEnU5IpIhkrmn8BBwUbN5NwIvuPtw4IVwGuBTBMcvhgPTCC7VnfEmjz6OE/p05b90Ex4R6SBJCwV3fxn4sNnsi4GHw+cPA5+Pm/9bD8wDepnZgGTVliqyY1lMO2co73ywi3lrmjeliEj7O9JlLpKhyN03hc83A0Xh84EEl9RoVBHO20QzZjaNYG+CoqIiysvL21RIZWVlm7ftSP3rnR65xh1//hv/PD6/3V8/Vdoh2dQOAbVDIJPboaNDoYm7u5m1uk/E3WcAMwDGjx/vZWVlbXr/8vJy2rptR/tWbDU/m72CfsNPZfTAnu362qnUDsmkdgioHQKZ3A4dPfpoS2O3UPhzazh/AzAobr3icJ4AV0wYTPe8bO6bo5vwiEhydXQozOLADXqmAk/Gzb8yHIU0Adgd182U8Xrk5/CVCYN5dvEm1m7fF3U5IpLGkjkkdSbBZTJKzazCzL4B/BS4wMxWAZ8MpwGeBdYQXHTv18A1yaorVX39zBKyY1nMmKub8IhI8iTtmIK7X36YRee3sK4TXHBPDqN/j3ymnFbMn+ZX8J3zh9O/R/sfdBYR0RnNKWTa2UOpa2jggVfXRl2KiKQphUIKKelXwOQxA3h03jr2VNVGXY6IpCGFQor51qRh7K2u45F566IuRUTSkEIhxYwe2JNzRhTywCtrqaqtj7ocEUkzCoUUdPWkYWyvrOZPCyqiLkVE0oxCIQVNGNqHsYN6MePlNdTV6yY8ItJ+FAopqPGy2us/3M+zSzZHXY6IpBGFQoq64GNFDOzVhaff2Rh1KSKSRhQKKSory+hTkEud7somIu1IoSAiIk0UCiIi0kShkOJ27a+hXl1IItJOFAop7NzSQhau38VX73+DrXuroi5HRNKAQiGFffeCEfxsysksXL+TyXfP5ZVV26MuSURSnEIhhZkZXxw/iFnXnUWvrrl89YE3+MX/rlB3koi0mUIhDYwo6s6s687k78cVc8+Lq/nyr+exZY+6k0Sk9RQKaaJrbjZ3XXoKd116CosqdjP57rm8vHJb1GWJSIpRKKSZKacVM+u6M+nbLZepD77JXc+t0PWRRCRhCoU0NLyoO09eexZfPG0Q//HSar786zfYvFvdSSJydAqFNNUlN8adU07ml186hSUbdzP5nrmUr9gadVki0skpFNLcF04tZtZ1Z9G/ex5fe/Bv3Dn7XXUnichhKRQywIn9u/E/157J5acP4t7y97j81/PYtPujqMsSkU5IoZAh8nNi/N9LTubuy8aybOMeJt89l5feVXeSiBxMoZBhLh47kKe+fRbH9ezCVQ/9jf9eUUOtupNEJBRJKJjZWjNbbGZvm9n8cF4fM3vezFaFP3tHUVsmGFrYjT9fcwZf+cQJ/OX9Wr70q9fZsEvdSSIS7Z7Cue4+1t3Hh9M3Ai+4+3DghXBakiQ/J8YdXxjD1afksXJLJZ++Zy5/XbYl6rJEJGKdqfvoYuDh8PnDwOcjrCVjfGJANk99+ywG9urCN387nzueWabuJJEMZu4df/E0M3sf2Ak48Ct3n2Fmu9y9V7jcgJ2N0822nQZMAygqKjrtsccea1MNlZWVdOvWra0fIW00tkNNvfPYihpeXF/H0J5ZXDM2j35dOtPfDMml70NA7RBI93Y499xzF8T10hwkqlAY6O4bzKw/8DzwbWBWfAiY2U53P+JxhfHjx/v8+fPbVEN5eTllZWVt2jadNG+HZxZt4sbHF2EGd116Cn930nHRFdeB9H0IqB0C6d4OZnbYUIjkT0F33xD+3Ar8GTgd2GJmAwDCnxovGYFPnzyAp6efxQl9uzLtdwu4/all1NSpO0kkU3R4KJhZgZl1b3wO/B2wBJgFTA1Xmwo82dG1SWBw3wIev/oMvnZGCQ+8+j6X3vcaH3y4P+qyRKQDRLGnUAS8YmbvAG8Cz7j7bOCnwAVmtgr4ZDgtEcnLjnHr507ivivGsWb7PibfM5fZSzZHXZaIJFl2R7+hu68BTmlh/g7g/I6uR47sotEDGDWgJ9+euZBvPbKAr51Rwk2TR5KXHYu6NBFJgswZXiJtdkLfrvzxW2fw9TOH8NBra5ly7+us36HuJJF0pFCQhORmZ/Evnx3Fr756Gut27OPT98zl2cWboi5LRNqZQkFa5cKTjuOZ6WcztH83rnl0IT96cgnVdfVRlyUi7UShIK02qE9X/viPE/nmWUN4+PV1/P29r7F2+76oyxKRdqBQkDbJzc7i5s+M4tdXjueDDz/iM//+Ck8v2hh1WSJyjBQKckwuGFXEM9PPYnhRN677/VvcOmtp1CWJyDFQKMgxK+7dlT/840T+flwxD722lo26DLdIylIoSLvIiWUxcVhfAOobOv56WiLSPhQKIiLSRKEgIiJNFAoiItKkw699JOnvgVffZ0DPfArysumWl03X3GwK8mIU5GZTkBc+z8umIDebWJZFXa6IxFEoSLsZ0q+AHvnZPPjq2oS3yc/JiguObApyYweCo4UQOWSd+Pl5MXJjWQQ37hORtlAoSLs5bXBvFt16IfUNzv6aOvZV17Ovpo591eHz6rpwup79NXVUVofLauoPWmfX/hoqdtaxv6a+aZ1EBzRlZ1mz4GghYMLnXXNjdAvXeW9rHV3W7DiwTbh919yYQkYyikJB2l0sy+ien0P3/Jx2eT13p7qugcrqOvZXB0FxIFQOBE9jiOyvrqOyWfDsqKwJtwvWaelucv+2cN4h88yga87BAdM1t7Fb7ECoNIVIXjbd4tY5KKDCbrTsmA7lSeelUJBOz8zIz4mRnxODdrqXem19QxAwNUGIvPz6m4wcfUpc4NSzPwyU+IBpDJWte6vYV30ghPbVJH5RwNzsrENCpfF5ECaxMFyCQGl8fvA62XTNC6bzstVlJu1HoSAZKSeWRc+uWfTsGuzNbOgV48wT+7X59RoanI9q65t1hwXdZY0B0xgq8d1ojYPgapUAAAgtSURBVCG0p6qOzburDtq+LsE+s1iWxQVG7KBjNPEBc/CyA91o8QHzUZ3T0OBkaQBAxlIoiLSDrMZjGXnt8yvl7tTUNzQ7FhN/bOZw8w8837Dro4PCqar20C6zFv31WbqEXWbxXWFdmw7sH3yAv3GdQw/+Hwih3Gx1maUKhYJIJ2Rm5GXHyMuO0acgt11es66+gf219YcNkX01dSxatpLjige3uMezo7KG9R/uD47fhF1vnuAAgNxYVhAquYcGR6Ijy+Kfd8nRAIBkUSiIZIjsWBY9Yln0OMIAgIEfvU9Z2YiEXs/dqaptiBtF1uzAf3X9oSPMmg0K2La3OlinJli/pj6xvRkzWg6YZiHSeOD/iMGTl03XHA0AaKRQEJE2MTO65MbokhujsHteu7xmTV3DIQf14/dsmgYBNK4T7rE0htCm3VUHrbO/FQMA8nOymoLDa6soWv7agVCJC5sWR581C56uubGUHQCgUBCRTiM3O4vc7Fx6dW2fLrPGc2YODFc+eEhz44H/+BDaX13Huo1byM3OYvf+GjYcwzkzB40wa36MpqXRZ0dYp0tOrEMGACgURCRtxZ8zU9SK7crLyykrm3DI/MZzZpr2XmrqWgyV+OMxzYc076jcf9CxnOoWzplpSeM5M42jyb7yiRP45tlDW/GpEqNQEBFJUPw5M33b+ZyZppFkzQ7wNw0EOOjYTD39urVPl11zCgURkQg1P2cmap3ucLuZXWRmK8xstZndGHU9IiKZpFOFgpnFgP8EPgWMAi43s1HRViUikjk6VSgApwOr3X2Nu9cAjwEXR1yTiEjG6GzHFAYCH8RNVwCfiF/BzKYB0wCKioooLy9v0xtVVla2edt0onYIqB0CaodAJrdDZwuFo3L3GcAMgPHjx3tZWVmbXicYcta2bdOJ2iGgdgioHQKZ3A6drftoAzAobro4nCciIh2gs4XC34DhZjbEzHKBy4BZEdckIpIxOlX3kbvXmdl1wHNADHjA3ZdGXJaISMYwT/Tat52QmW0D1rVx837A9nYsJ1WpHQJqh4DaIZDu7TDY3QtbWpDSoXAszGy+u4+Puo6oqR0CaoeA2iGQye3Q2Y4piIhIhBQKIiLSJJNDYUbUBXQSaoeA2iGgdghkbDtk7DEFERE5VCbvKYiISDMKBRERaZKRoZBJ92wws0Fm9pKZLTOzpWZ2fTi/j5k9b2arwp+9w/lmZveEbbPIzMZF+wnaj5nFzOwtM3s6nB5iZm+En/W/w7PoMbO8cHp1uLwkyrrbk5n1MrM/mdm7ZrbczCZm6Hfhu+HvwxIzm2lm+Zn4fWhJxoVCBt6zoQ64wd1HAROAa8PPeyPwgrsPB14IpyFol+HhYxpwb8eXnDTXA8vjpu8EfunuJwI7gW+E878B7Azn/zJcL13cDcx295HAKQTtkVHfBTMbCEwHxrv7aIKrJ1xGZn4fDuXuGfUAJgLPxU3fBNwUdV0d+PmfBC4AVgADwnkDgBXh818Bl8et37ReKj8ILq74AnAe8DRgBGesZjf/XhBcZmVi+Dw7XM+i/gzt0AY9gfebf5YM/C40XqK/T/jv+zRwYaZ9Hw73yLg9BVq+Z8PAiGrpUOFu76nAG0CRu28KF20GisLn6do+/wZ8H2gIp/sCu9y9LpyO/5xNbRAu3x2un+qGANuAB8NutN+YWQEZ9l1w9w3AXcB6YBPBv+8CMu/70KJMDIWMZGbdgMeB77j7nvhlHvwJlLZjk83sM8BWd18QdS0RywbGAfe6+6nAPg50FQHp/10ACI+ZXEwQkscDBcBFkRbViWRiKGTcPRvMLIcgEB519yfC2VvMbEC4fACwNZyfju1zJvA5M1tLcIvX8wj61nuZWeOVguM/Z1MbhMt7Ajs6suAkqQAq3P2NcPpPBCGRSd8FgE8C77v7NnevBZ4g+I5k2vehRZkYChl1zwYzM+B+YLm7/yJu0Sxgavh8KsGxhsb5V4YjTyYAu+O6FlKSu9/k7sXuXkLw7/2iu38FeAmYEq7WvA0a22ZKuH7K//Xs7puBD8ysNJx1PrCMDPouhNYDE8ysa/j70dgOGfV9OKyoD2pE8QAmAyuB94D/E3U9Sf6sZxF0BywC3g4fkwn6RF8AVgF/BfqE6xvB6Kz3gMUEIzQi/xzt2B5lwNPh86HAm8Bq4I9AXjg/P5xeHS4fGnXd7fj5xwLzw+/D/wC9M/G7ANwGvAssAX4H5GXi96Glhy5zISIiTTKx+0hERA5DoSAiIk0UCiIi0kShICIiTRQKIiLSRKEg0oyZ9TWzt8PHZjPbED6vNLP/iro+kWTSkFSRIzCzW4FKd78r6lpEOoL2FEQSZGZlcfdiuNXMHjazuWa2zswuMbOfmdliM5sdXloEMzvNzOaY2QIze67xchJHeI9JcXspb5lZ9474bCKNFAoibTeM4DpKnwMeAV5y9zHAR8Cnw2D4d2CKu58GPADccZTX/GfgWncfC5wdvpZIh8k++ioichh/cfdaM1tMcKOW2eH8xUAJUAqMBp4PLrFDjOBSzUfyKvALM3sUeMLdK5JRuMjhKBRE2q4awN0bzKzWDxygayD43TJgqbtPTPQF3f2nZvYMwfWpXjWzC9393fYuXORw1H0kkjwrgEIzmwjBJczN7KTw+XVmdl3zDcxsmLsvdvc7Ca7oO7JDK5aMp1AQSRJ3ryG41PKdZvYOwRVqzwgXj6Tla/J/J7yZ/CKgFvhLhxQrEtKQVJEIhKOYLgmDQ6TTUCiIiEgTdR+JiEgThYKIiDRRKIiISBOFgoiINFEoiIhIE4WCiIg0+f+IQp75WsBAiQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "# Plot mass\n",
    "figm, axsm = post.plot_single_variable(\n",
    "    x * 1e-3,\n",
    "    t,\n",
    "    [[-1]],\n",
    "    axis=0,\n",
    "    fig=None,\n",
    "    axs=None,\n",
    "    tics=[\"-\"] * 15,\n",
    "    name=\"mass\",\n",
    "    ylabel=\"Mass in tons\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot height and velocity profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAAFzCAYAAAD8AIVCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeXiU1f3+8feZmSTsYQ9bIIDsO4RFEI2KiAiouOKKWrFq1br8tGq/1draWqu1VOtWFxBEEVxYFBTQsENYJOw7gYQk7ARCSDLL+f2R0FIrmIRMnpnJ/bquXJl5Msv9GJx85sw5n2OstYiIiIiISOm5nA4gIiIiIhKuVEyLiIiIiJSRimkRERERkTJSMS0iIiIiUkYqpkVEREREykjFtIiIiIhIGXmcDnAu6tevbxMSEkp9vxMnTlC9evXyDxSidL6RTecbvlatWnXQWtvA6RwVRa/ZJaPzjWw63/B1ptfssC6mExISWLlyZanvl5ycTFJSUvkHClE638im8w1fxpjdTmeoSHrNLhmdb2TT+YavM71ma5qHiIiIiEgZqZgWERERESkjFdMiIiIiImWkYlpEREREpIxUTIuIiIiIlJGKaRERERGRMlIxLSIiIiJSRiqmRURERETKSMW0iIiIiEgZqZgWERERESkjFdMiIiIiUjmkp8DCV4q+lxNPuT2SSCVU6AtwMLeA/ccLyCv0/ft41Sg3DWrGUL9GDFWyV0HaQkgYCPF9HEwrIiJSiaWn4B83HFfAi3FHwx3Ty+XvsoppkRLy+QOs3H2ElWmHWbc3h/V7j7H36Mmz3qen2cqkmD8RhY+AiWL5hR/Qofcg6tWIqaDUIiIiAnBiazIx/kIMAay/EJO2UMW0SLD5/AG+33KAr9Zm8v2WA+Sc9ALQsn51ejSvzfWJzWhYswoNa8ZQPcaDMUX3yyv0ceB4Ac03LCYqzYebADbgZcncadz6DbRvVJOkdg0Z1rUxnZrUwpy6o4iIiJQ7ay2vbmvIY9aD2+UvGplOGFguj61iWuQnHMwt4JOUPUxavofMnHzqVItiUIc4BnVoSP/z6hNbNapkD9Toahj/PvgLcXuiuWr4DVQ/nsCibQd5d+FO3pq/g5b1q3NNj6bc2DueuFpVgntiIiIildD4JWm8m9aAnhe+y9Ca28t16qWKaZHT5Jz08q8FO3l/8S7yCv1ccF59fje8E4M6NMTjLsN63fg+RXOy0hZiEgbSLr4P7YAHLj6PIycKmb0hmxmpmfxtzlb+MW8bgzvFMbp/S/q0rFvepyYiIlIpbco6xp9mbeaS9g254opEKOdPg1VMiwD+gGXckjT+MW8bOSe9DOvamF8PasN5DWue+4PH9/nJd791qkczqk9zRvVpTtrBE0xK2cOnK9P5el02vRPq8MDF53FR2waaAiKOM8a8DwwD9ltrOxcfqwtMBhKANOAGa+0RpzKKiPyUfK+fhz7+gdiqUfz1uq5B+Zuq1nhS6W3KOsbINxbzh5kb6doslpkPXsDrN/csn0K6hBLqV+fpoR1Y+ptLeW54RzKOnGT0Byu4+o0lpOw6XGE5RM5gHDDkR8d+A8yz1rYB5hVfFxEJKS98tYlt+3N55fpuQVv8r5FpqbT8Acs/v9/OP+ZtI7ZqFK+N6sGwro0dHQmuGu1m9ICW3Ny3BZ+vzuDVuVu54e2lDO4Yx5NXtKd1gxqOZZPKy1q7wBiT8KPDVwFJxZfHA8nAkxUWSkTkZ8zZuI8Jy3Zzz8CWXNi2QdCeR8W0VEo5WxfxzVdTSd7fgiu7JfHc8E7UqR7tdKx/i/a4uKlPc67q3pT3Fu3kzeQdDH51Abf1a8Fjg9tSs0oJF0CKBE+ctTar+HI2EOdkGBGR0x3YuJBNUyZxbYMePH75jz9YK18qpqXS2bJyLi1mjmKk9TGyajTu/v0xIVRIn65qtJtfXdKGm/o059U5Wxm/NI3Z67N5bkQnLu8Up/nUEhKstdYYY3/qZ8aYMcAYgLi4OJKTk0v9+Lm5uWW6X7jS+UY2nW/wVT+6iS5rfsf91ge5U1n7VYBjse2D9nwqpqVS+WZDNuunTeFhtw+PCYD1wu5F0Lyv09HOqn6NGF64pgvX9WrGU5+v45cTVzGoQxzPX9XJ6WhSee0zxjS21mYZYxoD+3/qRtbad4B3ABITE21SUlKpnyg5OZmy3C9c6Xwjm843+L771yw8tvjvPH561j0BA4OXQQsQpdKYsGw3901cxf56vXF7osG4oRybtleEHs3rMOPBC3h6aHsWbz/I5a8uYNFeL9b+5KCgSDBNB+4ovnwHMM3BLCIiAMxen83rOxthXVEV9nc+aCPTxph44EOK5tFZ4B1r7dgztVMyRZ9XjwWGAnnAaGvt6mDlk8rl73O38ve527i0fUOeu3kIZl9XSFtYrk3bK0qU28WYC1tzRefGPPZpKu+uO8yeCav408gu1Nc25RIExpiPKVpsWN8YkwE8C7wIfGqMuRvYDdzgXEIREdhzKI//NzWVVk17Y4dNg/QlFfJ3PpjTPHzAY9ba1caYmsAqY8wcYDRF7ZReNMb8hqJ2Sk8CVwBtir/6Am8Wfxc5J2PnbuPvc7dxXa9mvDiyS9HmK2fo/RxO4utW4+Mx/Xh6/Fy+2HKAIX9fwF+u7cqlHbQOTMqXtXbUGX50aYUGERE5gwKfnwcmrcYAr9/ck+i61SDh/Ap57qBN87DWZp0aWbbWHgc2AU0paqc0vvhm44Griy9fBXxoiywDahfPwxMps39+v51X527lul7NeOnarmXbxTCEuV2GK1pGMePBC2hQswp3j1/JB59Mxjf/ZUhPcTpe5ZCeAgtf0X9vEREHvfDVJtbtzeHl67sRX7dahT53hSxALO5P2gNYzpnbKTUF0k+7W0bxsazTjmlleBlU1vP9Js3Lx5sLOb+Jm6H1D7NgwXynowVFbm4uWZtX8Uhny0p2ctOm5zGbffjm/4W13f4Q1BXMTgilf881j26i85rfEYUP6/KQGoH/vUVEQt3MtZl8uLSon/TgTo0q/PmDXkwbY2oAnwG/ttYeO72V19naKZ2JVoaXXmU83xN12/Hx7NUM7dKIf9zUI+JGpE93+u93cPRqAvN8uAjgC/io4d1Hz6RfOhuwnIXSv+f06YtwWx8uEwAb/BXjIiLy37buO84TU9fSs3ltnhjizGBGUCsMY0wURYX0R9baz4sP7zs1feNH7ZT2AvGn3b1Z8TGRUtl2xM8jn64hsUUd/nZD94gupP9HwkBcnhiscePHw29W1eL5GRvx+gNOJ4tIk/e3wGs82DDsDCMiEu5yTnq5d8IqqkV7ePPWXkQ59Pc+mN08DPAesMla+7fTfnSqndKL/Hc7penAr4wxn1C08DDntOkgIiWSdvAEY1fn07R2dd65PZEqUW6nI1Ws+D5wx3RM2kJc8QPourYm7y/exYbMHP55S091+yhHmUdP8ubOetTv9jqjm2SEZWcYEZFwFQhYHp28hvTDeXw8ph9xtao4liWY0zwGALcB64wxa4qPPc2Z2yl9TVFbvO0Utca7M4jZJAKdKPDxiw9XAvDB6N7UDdFdDYOuuFNJFPBcAnRtFstTn69jxGuLeOu2XnRtVtvphBHhg8W7ABg0eDjUqdjFLiIild3YeduYt3k/z1/Vid4JdR3NErRi2lq7CDjTXsf/007JFu068UCw8khks9by5Gdr2Xkgl8cTq5BQv7rTkULGyJ7NaBtXk3snrOK6t5byp+KdFKXsjud7+SQlnaFdGtNMhbSISIWas3EfY+cVtby9rV8Lp+NoB0SJDOOWpDFzbRaPX96OjvUq2dSOEujcNJbpvxpAYos6PD4llWenrdc86nMweUU6xwt83DOwpdNRREQqlR0Hcnlk8hq6NI3lj1d35vTGFk5RMS1hb9Xuw7zw1SYu6xjHfRe1djpOyKpXI4YP7+rDLy5oyfilu7nl3eUcyi1wOlbY8foDvL9oF31b1tWUGRGRCnQ838uYD1cS7XHx1m29QmZdlIppCWvH8r089PEamtSuysvXdwuJd6ihzON28dthHfn7jd1JTT/KiNcXsyEzx+lYYeXz1Rlk5uTzS71xExGpMIGA5dFPU0k7lMfrN/egae2qTkf6NxXTEtaenbaB7GP5jL2pO7FVo5yOEzau7tGUKb88H3/Act2bS/lqrRrnlITXH+D177fTtVksSe0aOB1HRKTSePnbLczZuI9nhnagf+v6Tsf5LyqmJWzNXJvJFz/s5VcXn0eP5nWcjhN2ujarzfQHB9CxSS0emLSal7/ZQiBQqj2UKp1pazJJP3yShy5po09BREQqQnoK6yc/y7L5sxjVpzl3DkhwOtH/qJDtxEXKW3ZOPs98sZ5u8bX51SXnOR0nbDWsWYVJ9/Tld19u4PXvt7M5+xiv3tidmlU0yv9jPn+Af36/nY6Na3Fph4ZOxxERiXzpKQTGDae9r5CPY6Jw9zo/JAcyNDItYedUG7xCX4C/39jdsR2PIkWMx82L13bh9yM68f2WA4x8YwlpB084HSvkzJszkyuOTOLZHrkh+WIuIhJpcjZ9R8BfiMcEiDY+PHsWOx3pJ6kKkbAzPTWT+VsP8MSQdrRUP+lyYYzhjv4JTLi7DwdzCxjx+iIWbD3gdKyQkb9zKRctu5vHoqbQZ8GdkJ7idCQRkYh2osDHc2vr4sWDNW6MO7pop9kQpGJawsqRE4U8P2Mj3eNrc/v5CU7HiTj9W9dn+q8uoEntqoz+IIV3F+6kaD+lym31ghl4rA83AYy/ENIWOh1JRCRiBQKWhz9Zw/TDzdh82UTMJc/AHdOLdvgNQZozLWHlha83kXPSy8SRXXC79FF7MMTXrcZn9/XnsU9T+eNXm9iYdYw/XdMlZPp5VrT9x/L5567GJLqiAB+E8OiIiEgk+Ms3m5m7aR+/H9GJHv0TgMudjnRWKqYlbCzefpCpqzK4P6k1HRrXcjpORKse4+GNW3ry2nfbeXXuVnYcOME7t/UirlYVp6NVuD/P2swK33kcvGkKTY6uLCqkQ3R0REQk3E1ctpu35+/k1n7Nuf1857cKLwlN85CwUOgL8H/T1tOiXjUeurSN03EqBZfL8PCgNrx9Wy+27zvO8NcW8cOeI07HqlDLdx7iix/2MubCVjTpchEMfEyFtIhIkHy3eR+/m7aeS9o35LnhncJmsbeKaQkLHy5NY+eBEzw7vGOlnW7glMs7NeLz+wcQE+XixreXMXVVhtORKoTXH+DZ6RtoWrsqD1ys9osiIsG0LiOHX036gY5NavHaqB54wqhTV/gklUrrwPECxs7dRlK7BlzSPs7pOJVSu0Y1mf7ABSQm1OHxKak8P2MjPn/A6VhB9cb3O9icfZxnh3ekarTewImIBEvGkTzuGr+COtWieX90b6rHhNcsZBXTEvJe+XYLJ71+/m9YR6ejVGp1qkfz4V19uHNAAu8v3sWd41ZwNK/Q6VhBsSnrGK9/v40R3ZowuFMjp+OIiESsnJNe7vxgBfleP+Pu7E3DmuG3NkfFtIS0dRk5TF6Zzuj+CbRuUMPpOJWex+3i2eGdeOm6rizfeZin//EeB2b9OaL6Lhf6Ajw+JZXYqlE8N6KT03FERCJWgc/PvRNWknboBO/clkibuJpORyqT8BpHl0rF7lnO+imTuKhqGx4aNNjpOHKaGxLj6Wq3kPDV7/As8+Ff+Xfco2dExOK8v8zezIbMY7x9Wy/qVo92Oo6ISEQKBCxPTF3Lsp2HGXtTd85vXc/pSGWmkWkJTekpBMaP4PpjH/Iuz1PrwA9OJ5IfaZ+fSozx4zEBrK+QJfO+DPsNXuZs3Md7i3Yxun8Cl2t6h4hIUFhr+cNXG5m2JpMnhrTjqu5NnY50TlRMS0gK7FoI/kI8JoDb+rTjXChKGIhxR2ONm4Aripe3NOCeD1dx5ER4zqPedfAEj09JpXPTWjw1tL3TcUREItYbyTv4YHEad1/Qkvsuau10nHOmYlpC0iJfewqthwBujHacC03xfeCO6ZhLniHqzhkMu/Jq5m/dz5CxC1iy46DT6UolJ8/L3eNW4DLwxs29iPGoe4eISDB8nLKHv36zhWt6NOWZoR3Cppf02WjOtIScQl+A366sRmKtF3i593FoqR3nQlZ8H4jvgwHuag59WtbloU9+4JZ3l3PfRa155LK2RIV4r9ACn5/7J60i/UgeH/2iH83rVXM6kohIRJq9PotnvlhHUrsGvHRdV1yu8C+kQSPTEoI+WbGHPYfzGDHsalwXase5cNK5aSwzH7yAm3rH80byDq7+52I2ZOY4HeuMvP4AD076gcXbD/HiyK70aVnX6UgiIhFpyY6DPPTxGrrH1+aNW3qG/EBLaUTOmUhEyPf6+ef32+nTsi4XtW3gdBwpg2rRHv48sitv3dqLfccKuOr1xbzy7RYKfH6no/0Xnz/AY5+m8u3GfTw3vCPX9mrmdCQRkYi0fm8OYz5cRUL9arw/ujfVoiNrYkRknY2EvU9XprPvWAGv3tg9IuZRVWZDOjeiX6u6PD9zI699t53Z67N58dou9Grh8OhvegreHQv46+YGTE+rz5ND2jN6QEtnM4mIRKjsEwEeez+F2KpRfHhXX2pXi7yWoyqmJWQU+Py8mbyD3gl1OL9V+PablP+oXS2av93QneHdmvD05+u49s2lXNOjKb+5oj1xtRzY5So9BTt+OC5fIY9YDz0veo8hSeG/klxEJBTtPXqSv67Ix7qj+PDuPjSKDb/dDUtC0zwkZExZmUFWTj4PXdpGo9IR5uJ2DZn76EU8cHFrvlqbxcUvJzN27jaO53srNEf66m/w+wpxEyDG+BlSY1uFPr+ISGWx/1g+t/xrGXk+y4d39YnoXYw1Mi0hodAX4M3kHfRoXpsLzqvvdBwJguoxHv7f5e25MbE5f/p6E6/O3coHS3Zx74Wtuf38FlSPCd7LUV6hj1e+3coPy2swKdqD2/hxqeWiiEj5S0/hxNZk/rgqlv25LXisVxU6N411OlVQqZiWkPD56gz2Hj3JC9d01qh0hGterxpv3daLdRk5/G3OFv4yezNvJG/npt7x3H5+AvF1y681nT9gmZ66l5e/2creoye5td8g/F36YjKXFhXS6hQjIlJ+iqfSVfEV8hfrYfvQSRzKj/y+/SqmxXH+gOWt+Tvo0jRWHTwqkS7NYvngzj78sOcI7y3axfuL03h30S4uOK8+w7s24fJOjYitFgXpKUU7YJai+M0t8PHFD3sZvySN7ftz6dSkFq/c0I1+p+bit+4fxDMTEamcCrbPx3NqKp3LTxfvWpLp5XSsoFMxLY6bszGbtEN5/PPmnhqVroR6NK/D6zfXISvnJJOW7+HLNXt54rO1PPPlOm5slMWzR57CE/CCJxpun/6Tj+HzB9i2P5c16Uf5bvN+Fm47QL43QOemtXhtVA+u7NI4YjYHEBEJRScL/fxxfV1+az3EnD6Vbkee09GCTsW0OMpay1vzd9K8bjWGdG7kdBxxUOPYqjw2uB2PXtaWtRk5fL0uiybrZ2P8XowJ4PMW8vr7H/CJewSNNy7GZQxef4BDuYXsO5aPL2ABaBJbhRsT47mqR1N6xNfWGzQRkSAr8PkZM2ElizIbMeSycQyM2vyfTxN3JDsdL+hUTIujVqQdYU36UZ6/qhNujRwKYIyhW3xtusXXhi63Ysd/SsBfCC4P1c9Lovl+FzExHqyFGjEezmtQg0axVWgbV5OuzWJpWb+6CmgRkQpS4PPzwEerWbjtIC9d15WBifHAlU7HqlAqpsVR7yzYQZ1qUVzfK97pKBKK4vtg7piOSVuIK2Eg98T3oU1yMklJfZ1OJiJS6RX4/Nw/cTXzNu/nD1d35obEyvm3XMW0OGb7/uPM3bSfhy9tQ9XoyF/tK2UU30ddN0REQkyBz899E1fz3eb9/PHqztzar4XTkRyjYloc8+7CXcR4XNx+fuX9H1BERCTc5Hv93DdxFd9vOcCfrunCzX2bOx3JUSqmxRFHThTyxQ97GdmzKfVqxDgdR0REREog3+vnlxNXkaxC+t9UTIsjJq9Mp8AX4I7+CU5HERERkRLI9/q5d8Iq5m89wIsju3BTHxXSoGJaHODzB5iwdDf9WtWlfaNaTscRERGRn5Hv9TNmwioWbjvAX67two29VUif4nI6gFQ+czftZ+/Rk4zWqLSIiEjIyyv08YvxK4sK6ZFdVUj/iEampcKNX5JGk9gqDOoQ53QUEREROYtj+V7u+mAFq/cc4a/XdeO6Xs2cjhRyNDItFWpL9nGW7jzEree3wOPWPz8REZFQdfhEIbf8azlr0o/y2qieKqTPQCPTUqE+XJpGtMfFTfqISEREJGTtP5bPre8tJ+1QHu/c3otL2uvT5DNRMS0V5kSBj2lrMhnetQl1q0c7HUdERER+QsaRPG59dzn7jxcw7s7e9G9d3+lIIU3FtFSYmWszyS3wcXPfyrndqIiISKjbdfAEt/xrGccLfEy4uy+9WtRxOlLIUzEtFWbS8j20jatBz+b6H1NERCTUbMk+zi3vLidgLR/f04/OTWOdjhQWtAJMKsT6vTmkZuQwqk9zjDFOxxGRcmKMecQYs8EYs94Y87ExporTmUSk9FbtPswNby/FZWDyGBXSpaFiWirEJyv2EONxcU2Ppk5HEZFyYoxpCjwEJFprOwNu4CZnU4lIaX23eR+3vLucutWj+ey+/rSJq+l0pLCiaR4SdHmFPr78IZMruzSmdjUtPBSJMB6gqjHGC1QDMh3OIyKl8NmqDJ74bC0dGtdk3J19qF8jxulIYUfFtATdzNQscgt8jOqrdngikcRau9cY8zKwBzgJfGut/fb02xhjxgBjAOLi4khOTi718+Tm5pbpfuFK5xvZQul8Z+3yMnlLIR3qunigg4/1K5eW+3OE0vkGi4ppCbpPV6bTukF1ErUiWCSiGGPqAFcBLYGjwBRjzK3W2omnbmOtfQd4ByAxMdEmJSWV+nmSk5Mpy/3Clc43soXC+VpreXHWZiZv2cmVXRrztxu7EeNxB+W5QuF8g01zpiWodh08wcrdR7iuV7wWHopEnkHALmvtAWutF/gc6O9wJhE5C9/uZcx68wlWLJzNrf2a849RPYJWSFcWGpmWoPpsVQYugxYeikSmPUA/Y0w1iqZ5XAqsdDaSiJzJyR1LcU28isEBL5dVjcLTsz/GpYGuc6WRaQmaQMDy+eoMBrZpQKNYdcsSiTTW2uXAVGA1sI6ivynvOBpKRH7S/mP5fPrZJNwBLx4TIMr6MLsXOR0rIqiYlqBZuvMQmTn5XNurmdNRRCRIrLXPWmvbW2s7W2tvs9YWOJ1JRP7b1n3HueaNJczObYPxRINxgzsaEgY6HS0iaJqHBEd6Cke++ZQLqrRgcMchTqcRERGplJbsOMi9E1ZRJcrN02PuwG37QdrCokI6vo/T8SKCimkpf+kp2PHDGeItZLAriujsvvofVkREpIJ98UMGT0xdS0K96nxwZ2+a1akG9NHf5HKmaR5S/tIWYn2FRXOy8BW9AxYREZEKYa3ltXnbeGRyKokt6jL1vv7FhbQEg0ampfwlDMRLFG68uDUnS0REpMJ4/QF++8V6Jq9M55oeTfnLtV2J9mjsNJhUTEu52xfblfsKn+LRNvu5YNA1+jhJRESkAuTkefnVx6tZuO0gD15yHo9e1lZ7PFQAFdNS7makZrI60JbGw+6BBjWcjiMiIhLxdh7I5RfjV5J+JI+XruvKDYnxTkeqNFRMS7mbkZpJpya1aK1CWkREJOgWbTvI/R+twuN28dEv+tGnZV2nI1UqQZtEY4x53xiz3xiz/rRjzxlj9hpj1hR/DT3tZ08ZY7YbY7YYYy4PVi4JrrSDJ0jNyGFEtyZORxEREYl4E5amcccHKTSOrcq0BwaokHZAMEemxwGvAx/+6Pir1tqXTz9gjOkI3AR0ApoAc40xba21/iDmkyCYkZoJwDAV0yIiIkHj9Qd4fsZGJizbzSXtGzL2pu7UrBLldKxKKWjFtLV2gTEmoYQ3vwr4pHjnrF3GmO1AH2BpkOJJEFhrmZ6aSe+EOjStXdXpOCIiIhEpJ8/L/ZNWsXj7IcZc2Ionh7TH7dJCQ6c40SvlV8aYtcXTQOoUH2sKpJ92m4ziYxJGNmcfZ9v+XE3xEBERCZIdB3K55o3FpOw6zF+v68rTQzuokHZYRS9AfBP4A2CLv78C3FWaBzDGjAHGAMTFxZGcnFzqELm5uWW6X7iqqPOdsqUQl4HY47tITk4L+vOdiX6/ka2yna+IyClzN+7jkclriPK4mHRPP3onaH50KKjQYtpau+/UZWPMv4CZxVf3Aqf3cGlWfOynHuMd4B2AxMREm5SUVOocycnJlOV+4aoiztdayzPLvueCNrUZMdjZvtL6/Ua2yna+IiKBgGXsvG2MnbeNzk1r8datvbSjYQip0GkexpjGp129BjjV6WM6cJMxJsYY0xJoA6RUZDY5N6v3HGHv0ZOa4iEiIlKOjuV7uefDlYydt41rezZj6i+1NXioCdrItDHmYyAJqG+MyQCeBZKMMd0pmuaRBtwLYK3dYIz5FNgI+IAH1MkjvMxIzSLa4+LyTnFORxEREYkI2/YdZ8yEVaQfzuP5qzpxW78W2tEwBAWzm8eonzj83llu/wLwQrDySPAEApbZ67O5qG0DteUREREpB7PWZfH4lFSqRnuYdI82Ygll2gFRzllqxlGyj+Xz/zq1czqKiIhIWPMHLK98u4U3knfQPb42b93ai0axVZyOJWehYlrO2ez12XhchkEdNMVDRESkrA7mFvDI5DUs3HaQUX2a89yIjsR43E7Hkp+hYlrOibWW2RuyOb91PWKraYqHiIhIWaTsOsyDH6/mSJ6XF0d24aY+zZ2OJCWkYlrOyebs4+w+lMe9F7Z2OoqIiEjYCexezor50/nrlgZUrd2ND+7vQ8cmtZyOJaWgYlrOyaz12RgDl3XUFA8REZHSOL5tMTGTrqFXwMuk6Ci8I7+kugrpsOPEduISQb5Zn03vhLo0qBnjdBQREZGwsSb9KB9NnoQr4MVjAkTho3rmUqdjSRmomJYy23kgly37jjOkUyOno4iIiIQFay3jFu/i+reWsNrVCZcnGowb446GhIFOx5My0DQPKUhx2xsAACAASURBVLPZG7IBGNJZxbSIiMjPycnz8s81Bazct5FL2zfkpRsG4TqUCGkLiwrp+D5OR5QyUDEtZTZ7fTbdmsXSpHZVp6OIiIiEtJRdh/n1Jz+w75ifp65ozz0DW+FyGajWR0V0mNM0DymTvUdPsjYjhyGdGzsdRUREJGT5/AFenbOVm95ZSpTHxTP9qnDvRa2LCmmJCCqmpUy+Wa8pHiIiImez9+hJRv1rGWPnbePq7k356qGBtIrVJiyRRtM8pEzmbNxH27gatKxf3ekoIiIiIWfWuiye/Gwt/oDl1Ru7cU2PZk5HkiBRMS2llnPSy4q0w9xzYSuno4iIiISUk4V+np+5kY9T9tAtvjb/uKk7Lepp4CmSqZiWUluw9QC+gOXS9g2djiIiIhIyUtOP8sina9h18AT3JbXm0cvaEuXWjNpIp2JaSm3H6u94rOpiepj6QF+n44iIiDjK6w/w+nfbef377cTVjOGju/vS/7z6TseSCqJiWkrFv3s59+5+hGh8uCd8AXdMV0sfERGptLbvz+WxT9eQmpHDyB5NeXZEJ2KrRjkdSyqQimkplazUOTSyPtwmAP7CokbzKqZFRKSSCQQsHy5N48+zNlMt2s2bt/Tkii5qF1sZqZiWUvkuvx3X48Ft/Nr6VEREKqWsnJP8vylrWbT9IBe3a8Bfru1Kw1pVnI4lDlExLaUycW9DdjT4C7/vdkRbn4qISKVireXLNXt5dtoGfAHLn67pwqg+8RijDVgqMxXTUmLph/PYui+XG65MgoFqiyciIpVHdk4+z3yxjnmb99OrRR1eub4bCdprQVAxLaUwb9M+AAZ1iHM4iYiISMWw1jJlZQZ/+GojXn+A/xvWkdH9E3BrO3AppmJaSmze5v20alBd78RFRKRS2Hv0JL/5bC0Ltx2kT8u6vHRtV/0NlP+hYlpKJH/Ze4zZPY5jLa8AkpyOIyIiEjSBgGVSyh7+/PUmLPD8VZ24tW8LXBqNlp+gYlp+3spxxMx+lAsMsHstrGwCiaOdTiUiIlLustcvIPmbz/n8YAu6tz6fF0d2Jb5uNadjSQhTMS0/b9M0AIwBe+q6imkREYkgXn+A6TO/ZOjqe7kOH9dVjcY9eDpGhbT8DG0YLz8r0H4EUFRIG4AOVzkZR0REpFyt2n2EYf9YxI6U2UQbHx4TwGO9mN2LnI4mYUDFtPysbSaeb3yJHKndBYaN1ai0SJgyxrxkjKlljIkyxswzxhwwxtzqdC4Rp+Sc9PLbL9dx3VtLOJbvJenykbg9MWDcoI3JpIQ0zUPOLj2FVl/fTGu3F3duNMR1dDqRiJTdYGvtE8aYa4A0YCSwAJjoaCqRCmat5et12Tw3YwOHcgu4s39LHh3clhoxHmg5HdIWamMyKTEV03J2aQtxBby4TQD83qIXGL24iISrU6/5VwJTrLU52rlNKpv0w3k8O30D323eT+emtXj/jt50aRb7nxvE99HfOSkVFdNyVvvq9qaW9RBj/Lj0kZdIuJtpjNkMnATuM8Y0APIdziRSIfK9ft6ev5M3krfjdhl+e2UHRvdPwOPWjFc5NyUupo0xtU6/vbX2cFASSUiZdTSe6YVP8+5FBdTtdInerYuEIWNME2ttprX2N8aYl4Aca63fGJMHaEWxRLy5G/fx/MyN7Dmcx5VdG/PbKzvQOLaq07EkQvxsMW2MuRf4PUWjF7b4sAVaBTGXhIh5m/dztH4P6g5JcjqKiJTdu8aYukAyMBtYBGCtPQGccDCXSFDtPnSC38/YyHeb93Newxp89Iu+DDivvtOxJMKUZGT6caCztfZgsMNIaDme72XZzkPcOaCl01FE5BxYa4caY6pQtH3pNcDLxpg9FBXWs621e5zMJ1LeThb6eTN5O28t2EmUy/DM0A6MHpBAlKZ0SBCUpJjeAeQFO4iEnnXL5nAPXzKs3g1AB6fjiMg5sNbmU1w8AxhjWgJXAK8bYxpZazWHS8KetZZZ67N54atN7D16kqu7N+HpoR1oWKuK09EkgpWkmH4KWGKMWQ4UnDporX0oaKnEeekpJM4fTZ8oL+4506HJdM2XFokQxWtgcoBPir9ynU0kcu7WZhzlDzM3siLtCO0b1WTymH70bVXP6VhSCZSkmH4b+A5YBwSCG0dCRWDXQlzWi4cA+AvVEk8kApxpDYy1tsxrYIwxtYF3gc7Fj3mXtXbpuWYVKamsnJP8dfYWPv9hL/VrRPPnkV24ITEet0ttH6VilKSYjrLWPhr0JBJSNsd0o6X14FJLPJFIEow1MGMpmnd9nTEmGqhWjo8tckZ5hT7enr+TtxfsIGDhvqTW3J/UmppVopyOJpVMSYrpWcaYMcAM/nuah1rjRbCp+5uw3v9bJlxaSMx5F2lUWiQylOsaGGNMLHAhMBrAWlsIFJbX44v8lEDA8vkPe/nrN5vZd6yAYV0b8+SQ9sTX1fs4cUZJiulRxd+fOu2YWuNFMGst32zIpkOb/sRc3NvpOCJSfsp7DUxL4ADwgTGmG7AKeLi45R4AxYMxYwDi4uJITk4u9ZPk5uaW6X7hSuf7v2rlbCb26DpW0ZHXMlqRfjxAq1gXz/StQps6x9ixNoUdFRP3nOn3G3lKUky3ttb+11zp4hZLEqHW7z3G3qMn+fWgNk5HEZHyVd5rYDxAT+BBa+1yY8xY4DfA/526gbX2HeAdgMTERJuUlFTqJ0lOTqYs9wtXOt8fSU8hMO5ZrL+QJtbD91We5/GbRjC8axNcYTgvWr/fyFOSYvpd4K5TV4wx1YHpwKXBCiXOmrU+C7fLMKhDnNNRRKR8lfcamAwgw1q7vPj6VIqKaZFysX3/cdZ98QnDfYV4TIAY4+e1/ifwdG/qdDSRfytJ9/K9xpg3AIwxdYA5wMSgphLHWGuZvT6b81vVo071aKfjiEj5mmWMGWOMaWyMqXvqq6wPZq3NBtKNMe2KD10KbCyXpFKpZR49yZNT1zL41QVMPdQS647CGjcuTzSeVhc6HU/kv/zsyLS19v+MMS8ZY94CegEvWms/C340ccK2/bnsPHiCOy/QrociESgYa2AeBD4q7uSxE7jzHB5LKrkDxwt4a/4OJi7bjbUwun9LHrh4EFFH+ha1aE0YqAXxEnLOWEwbY0aednU5RXPgUgBrjBlprf082OGk4s1en40xcHlHTfEQiUDlvgbGWrsGSDynVFLpHcwt4O35O5iwbDeFvgDX9GjGrwe1+U+Hjhp9VERLyDrbyPTwH13/AYgqPm4BFdMRaNea7/ljvQ00zGkAtfTCJRJhtAZGQsqh3AImbynkvnnfU+Dzc3X3pjx4aRta1q/udDSREjtjMW2t1Ud1lUzW+vn86dgzxBgfjJ8Md2gLcZEIs9cY84a19v7iNTBfAf9yOpRUPodPFPLOgp18uDSNk4V+rurehAcvbUPrBjWcjiZSaiXp5iGVxJ7V39IAHy5tIS4SkbQGRpy2/1g+7y3excSlu8nz+hnetQn9ah7h5mE9nI4mUmYqpuXfvjjSku7Gg8f4QVuIi0QMrYERp+05lMdbC3YwdVUGPn+AK7s24aFLzqNNXM2I39BDIp+KaQEgOyefT7Iak9jvLa6rl6YV0yKRRWtgxBGbso7xZvIOZq7NxONycW2vZtx7YSsSNCdaIkipi2ljzFVA9mlN+iUCfLUuC4AeAy4HzVkTiShaAyMVJj0F0hayuUo3/rohlnmb91M92s0vBrbi7gtaEldLGyhL5CnLyHRfoIsxxmOtvaK8A4kzZqRm0rFxLS3+EBGRMvHtXob58CrwF9LCevC5fsejlw3hjvMTiK0W5XQ8kaApdTFtrX06GEHEOemH81iTfpQnh7R3OoqIiISZY/leJqek418wgV+ctu33vy7KJzqpjdPxRILuZ4tpY0w14DGgubX2HmNMG6CdtXZm0NNJhZixNhOAYV0bO5xERETCxe5DJ/hgcRpTVqZzotDPrU17Yo5OxQa8uNzRRLe+yOmIIhWiJCPTHwCrgPOLr+8FpgAqpiPEjNQsejSv/Z+dpkSkUtAaGCktay1Ldxxi3JI05mzah8dlGN61CXdd0JLOTWMhvYu2/ZZKpyTFdGtr7Y3GmFEA1to8Y4wJci6pINv3H2dT1jF+N6yj01FEpOJpDYyUyNG8QqauymDS8j3sPHiCOtWiuD+pNbefn/Dfiwrjte23VD4lKaYLjTFVKWqfhDGmNVAQ1FRSYWakZmGMpniIVEZaAyNnY61lTfpRPlq+hxmpmRT4AvRsXpu/3dCNoV0aUyXK7XREkZBQkmL6WWA2EG+M+QgYAIwOZiipGHbPcmqsmMBtTXvSUO2KRCKeMWYV8D4wyVp7xOk8EppyC3zMSM1k4rLdbMg8RvVoN9f1asYtfVvQsUktp+OJhJyfLaattXOMMauBfoABHrbWHgx6Mgmu9BTs+BHc6SuEw1MgvbM+mhOJfDcCdwIrjDErKVoT86211jobS5wWCFg2pMxlz+pvmJjVnKXe1rRvVJM/XN2Zq7s3oWYVtbYTOZMzFtPGmJ4/OpRV/L25Maa5tXZ18GJJ0KUtBH9RCyNrvUXXVUyLRDRr7XbgGWPM/wHDKBql9htjPgDGWmsPOxpQKlz64Tw+W53B5hVzeTX/WTrg4zJPFDtHTKJd4kC0RErk551tZPqV4u9VgEQglaKR6a7ASv7T3UPCkG1xAYV4iMKH2x1dtPJaRCKeMaYrRaPTQ4HPgI+AC4DvgO4ORpMKklvg45v12UxdlcHSnYcwBl5ssJWYAh8uAoCP9vmpYAY5HVUkLJyxmLbWXgxgjPkc6GmtXVd8vTPwXIWkk6D5gbb8seBpnu92lM4DrtSotEglUDxn+ijwHvAba+2pxeTLjTEDnEsmQZWegnfHAlaZjkzIaMS8zfvI9wZoUa8aj13WlpG9mtH0eAMY/wn4C0EDLCKlUpIFiO1OFdIA1tr1xpgOQcwkFWBGaibr3e1pfvUg0Fw4kcriemvtztMPGGNaWmt3WWtHOhVKgmDlOAIbvyQ9pg2NN4/HFfDSDQ9ve57lhsQkhndrQmKLOv+ZxlG7D9wxXT2iRcqgJMX0WmPMu8DE4uu3AGt/7k7GmPcpmpO331rbufhYXWAykACkATdYa48U960eS9HHjnnAaM3JDh5/wPLV2iwubteAWiqkRSqTqcCP18NMBXo5kEWCZeU47MyHMUBz+z3WGFzG4jJ+3r0wH/dFnX/6fuoRLVImrhLc5k5gA/Bw8dfG4mM/Zxww5EfHfgPMs9a2AeYVXwe4AmhT/DUGeLMEjy9ltHzXIfYfL2B4tyZORxGRCmCMaW+MuRaINcaMPO1rNEXrYiSSbJoGFC1ywlA0+mzcuNzRuFtd6Gg0kUhUktZ4+cCrxV8lZq1dYIxJ+NHhq4Ck4svjgWTgyeLjHxa3Z1pmjKltjGlsrc1Cyt2M1CyqRbu5pH1Dp6OISMVoR9EnhbWB4acdPw7c40giCZpAXBfMju+wFBfU/R+CKrU0fUMkSH62mDbG7KJ498PTWWtbleH54k4rkLOBuOLLTYH0026XUXxMxXQ58/oDzFqfxaAOcVSLLsksHxEJd9baacA0Y8z51tqlTueRIEpPwS5/m4AFl3FhBjwEl/3e6VQiEa0k1VTiaZerANcDdc/1ia211hhT6o0CjDFjKJoKQlxcHMnJyaV+7tzc3DLdL1ydfr5rD/g4muelpftQxP43qMy/38qgsp1veTDGPGGtfQm42Rgz6sc/t9Y+5EAsCYa0hRh/IW4D1piiEWkRCaqSTPM49KNDfy9ur/S7MjzfvlPTN4wxjYH9xcf3AvGn3a5Z8bGfyvMO8A5AYmKiTUpKKnWI5ORkynK/cHX6+U7/dA01q+zj/pEXE+NxOxssSCrz77cyqGznW042FX9f6WgKCbpAc+0hIFLRSjLN4/SV3y6KRqrLOj9gOnAH8GLx92mnHf+VMeYToC+Qo/nS5S/f6+fbDfu4onOjiC2kReR/WWtnFH8f73QWCa4faMMLBU/zh+5H6dRfewiIVISSFMWvnHbZB+wCbvi5OxljPqZosWF9Y0wG8CxFRfSnxpi7gd2nPc7XFLXF205Ra7ySdAuRUkrecoDcAh8juquLh0hlZIyZQ1Gv6aPF1+sAn1hrL3c2mZSXr9dls97VnuZXaQ8BkYpSkmL67p9q8v9zd7LW/s+8vGKX/sRtLfBACbLIOVi/fA6PV13C+dH1gAZOxxGRitfgVCENUNznX219IkQgYJm1LosL29anpgppkQpTkj7TU0t4TELYyR1LeWDPo9xnP8Ez4WpIT3E6kohUPL8xpvmpK8aYFvxEtyYJT5tWzOXq3MncHr//528sIuXmjCPTxpj2QCeKm/yf9qNaqMl/2Nm1ajZt8eEmAP7Coi1jNZdOpLJ5BlhkjJlPUQvigRR3R5Iwl55Cm9m38KjHi3vJdGgzXa/xIhXkbNM81OQ/gnx5pBWPGg9u48dohbdIpWStnV28qLxf8aFfW2sPOplJyod3xwJcAS8eowETkYp2xmJaTf4jxwmv5YM9DWjc5XXubJqhXbBEKrf+wOl7Ss90KoiUn+W2A73w4DJ+XBowEalQZ5vmoSb/EWL1Ph9ev6XngMshvrbTcUTEIcaYF4HewEfFhx42xvS31j7tYCwpBx/sbsi4qN/zzoUnoaUGTEQq0tmmeajJf4RYnuWned1qdG0W63QUEXHWUKC7tTYAYIwZD/wAqJgOYwdzC5i/9QB3D7wY14UdnI4jUumcbZrHjOKLedbaKaf/zBhzfVBTSbk5mFvAxsN+7ktqjDHG6Tgi4rzawOHiy3qHHQFmpGbiC1hG9mjmdBSRSqkkrfGeKuExCUGz1mcTsDC8mzZqERH+DPxgjBlXPCq9CnjB4UxyDqy1TFmZQacmtWjXqKbTcUQqpbPNmb6Coo8Emxpj/nHaj2pRtBOihIEZqZk0qWFoF6cXWZHKzlr7sTEmmaJ50wBPWmuzHYwk52htRg4bs47xh6s7Ox1FpNI625zpTIpGLUYUfz/lOPBIMENJ+cjKOcmKtMNc3TpKUzxEKrHidninyyj+3sQY08Rau7qiM0n5WJL8NQ9Fz2dkgzpAC6fjiFRKZ5sznQqkGmMmWms1Eh2GvlqbhbXQt3FJdo0XkQj2yll+ZoFLKiqIlJ+8HUsZvf1hol0+3B9Pgzu0UYuIE842zWMdxdvM/mhU0wDWWts1uNHkXM1IzaRz01o0qu53OoqIOMhae7HTGaT8bV0+i87a2VbEcWcbshxWYSmk3GWtn8+ArI9o228IUNfpOCISAowx1YBHgebW2jHGmDZAO2utNm4JQxOz4/mjdrYVcdzZpnns/qnjxpgLgFHAA8EKJecoPYX6n1/Pox4vrtRprOnyeyDJ6VQi4rwPKFoD07/4+l5gCtoFMeysy8hh6v4mXHLRuwytsV0724o4qESTaY0xPYCbgeuBXcDnwQwl5yhtIa6AF7cJgN9L7aPrnU4kIqGhtbX2xlO72lpr84xWJ4elcUvSqBbtZkDSYKga5XQckUrtbHOm21I0Aj0KOAhMBozm3oW+jNhe1LMeYowflzuao7XVMklEACg0xlTlP+thWgMFzkaS0jpwvIAZqZnc2DueWBXSIo4728j0ZmAhMMxaux3AGKOWeGFgyr4mLPQ+zYeXFFKj3cUc25HndCQRCQ3PAbOBeGPMR8AAYLSTgaT0Ji3fQ6E/wOgBCU5HERHOXkyPBG4CvjfGzAY+oaiTh4Qway0z1mYSl3A+NQb1Kzq4I9nRTCLiLGPMP4FJ1tpvjTGrgH4UvZ4/bK096Gw6KY1CX4CJy3dzUdsGtG5Qw+k4IsJZthO31n5prb0JaA98D/waaGiMedMYM7iiAkrpbMo6zs4DJxjWrbHTUUQkdGwFXjbGpAFPApnW2pkqpMPP0vlfc33epzzc7ojTUUSk2BmL6VOstSestZOstcOBZsAPFL0YSwiasTYTt8twRWcV0yJSxFo71lp7PnARcAh43xiz2RjzbPH6GAkDds9y+i28k8eiptDj+zsgPcXpSCJCCYrp01lrj1hr37HWXhqsQFJ21lpmpGYy4Lz61K0e7XQcEQkx1trd1tq/WGt7ULS4/Gpgk8OxpIQy1szBbYs2aTGnNmkREceVqpiW0JaakUPGkZMM76pRaRH5X8YYjzFmePHiw1nAForWx0gYGJfZDJ/xYI0btEmLSMgoUZ9pCQ8zUjOJdrsY3KmR01FEJIQYYy6jaCR6KJBC0YLyMdbaE44GkxLbkJnDe2kN6NDvLa6rl6ZNWkRCiIrpCBEIWL5am8WFbRuo76iI/NhTwCTgMWutVq6Fobfm76RGjIfLLh+hTVpEQoymeUSIlbuPkH0sn+Hq4iEiP2KtvcRa+26wCmljjNsY84MxRtuSB0HawRN8tTaTW/o112CJSAhSMR0h1iz5loeipzO45h6no4hI5fMwWsgYNG8v2InH7eLuAS2djiIiP0HFdATw7V7G7dse5GHXp1T9+Bq1SxKRCmOMaQZcCbzrdJZIlHn0JJ+tyuC6Xs1oWKuK03FE5CeomI4AGavn4Clul4TaJYlIxfo78AQQcDpIJHrtu+1YLPcntXY6ioicgRYgRoCvj7fmTjy4jR+jdkkiUkGMMcOA/dbaVcaYpDPcZgwwBiAuLo7k5ORSP09ubm6Z7heuTp2vf99G6q5bw+iGXdieWo3tTgcLksr6+60sKsP5qpgOc4W+AG/vqk9hwqv8+rx9apckIhVpADDCGDMUqALUMsZMtNbeeuoG1tp3gHcAEhMTbVJSUqmfJDk5mbLcL1wlJyeT1LoahfOf4yK3F/fxaZjWMyL2tb1S/n51vhFFxXSYW7z9IDknvXTpdxl0iHM6johUItbapyhqu0fxyPTjpxfSUnaHN8yjVsCLxwTA7y2avhehxbRIuNOc6TA3c20WNat4GNimgdNRRESknIzPiseLdjsUCQcamQ5jXn+AuZv2cVmHOKI9el8kIs6x1iYDyQ7HiAi7j/n5x9Y61OvxT25vnK7peyIhTsV0GFu64xA5J70M6aztw0VEIoG1lslbComtGsVVw66BatqkRSTUaTgzjM1an021aDcXttUUDxGRSJC89QAbDwV46JI2xKqQFgkLKqbDlD9gmbMxm4vbN6RKlNvpOCIico58/gB/+moTDasZbu3Xwuk4IlJCmuYRpjatmMsNJ6fQv/FVTkcREZFzlZ7ChoUzqXmgDoO7ddE6GJEwomI6HKWn0G72LbT3eHEvng7nTdfiFBGRcJWegh0/gk6+Aj6OiWJ9ld8Dg5xOJSIlpLe+YSiwayHGFvUfNdo+XEQkvKUtxPoK8BAg2viok7PB6UQiUgoqpsPQtmrd8FoPAdR/VEQk3O2p1YsC68GPC+OO5mjtzk5HEpFS0DSPMPT5/qas9j3DhEu9VGlzkaZ4iIiEKWstTyyPIcr8jrcGnqR62ySO7chzOpaIlIKK6TBjrWXW+mwSWvenyiUqokVEwtn01EyW7TzMH68eQfVTHTx2JDuaSURKR9M8wszGrGPsOZzHFdqoRUQkrB3P9/LCV5vo0jSWUX2aOx1HRMpII9NhZvb6bFwGBneMczqKiIicg5dmb+FAbgFv39YLt8s4HUdEykgj02Fm9vps+rSsS70aMU5HERGRMtqwfA41V/yD33Y9To/mdZyOIyLnQCPTYSRjbTKXHZpE5/ZXOh1FRETKqGDXUlrPuplHo3y4d0yH9FZaSC4SxjQyHS7SU4j78gYe9UzhitVjID3F6UQiIlIGS+dNw2N9eNBeASKRQMV0uEhbiCtwaqMWr158RUTC0NqMo7y2M46AKwqM9goQiQSa5hEmsuv0JtZ6iDF+XHrxFREJOycL/Tz6aSrHa3Sh8IYvic5aWvRarikeImFNxXSYmH64KbMLn+b9pAJqd7xEL74iIuEkPYVFX0+l1oHGPHvXrdQ4rwGc19/pVCJSDlRMh4P0FKotn0irBp2pffkjTqcREZHSSE/BP244F/sKSaoSRVSVAUADp1OJSDlRMR3q0lOwHwzlFr+XUS4PpHfXqLSISBg5sTWZGH8hHhPA4ita86LXcZGIoQWIoW7xWAh4MQZc1ld0XUREwoK1lrHb4/BaD9a4MVrzIhJxNDId6o5n/fui+dF1EREJYekprPh+GivT6jN3wL8YHrtDCw5FIpCK6RCX1+R8qmaswpriYrrH7U5HEhGRn1M8T7qnr5BPqkQR1XUGNL/a6VQiEgSa5hHK0lOIXvUufsDiggG/hsTRTqcSEZGfcXLFRIw/H48JEIUPs3uR05FEJEhUTIeytIWYQCEeA8YYqFLL6UQiIvIz/N/+jpi14zEWLGBcHs2TFolgKqZDWE5cXwqtBz9atCIiEhZWjsO1ZCzGgjFgMNDjZs2TFolgmjMdwpbvPMx+/0CGdmlM3f536MVYRCTUbZoGFBXSluJPFbvd7GwmEQkqjUyHqvQULlr+C0Z5vqfOts+cTiMiIiXgbdAZTk3vAOh8nQZCRCKciukQdWJrMu6Al//f3r2HS1XXexx/f+eyIUVJvOwQUS6ShBaCOwWFI2UF4YUe9ViahWWHbj5aeTpH63R7OnWs09X0VFQq5QXLG4QIIrkRERhAEESQ+3bYgBqJshXYc/meP2Zt2xHCZpg9a83M5/U88zBrze372789X777t37rt+LksVxrYZF/ERGJtE0tCXIEhTQGxw0MNyAR6XQqpiNqXmYgGQqL/KP50iIiFWF6S38yVlfI3Ymuyt0iNSCUOdNmtgnYCeSArLs3mFkP4F6gD7AJuMzdXwkjviiYtPk4phz2XW4Z/jr01SL/IiJRl8nl+d2mY7E+P+Pa/tt0gRaRGhHmyPT73P10d28Itm8AZrv7AGB2sF2TtrfsYf767fQ5fRT2L9crGYuIVICVC2dxZeZ+hvU7GkYqd4vUiiit5jEOGBXcnwQ0Av8ZVjBhenjFVvIOF7zn+LBDVy3b1wAAFZhJREFUERGRjkinGDTrE5yWyBB/Yir0napiWqRGhFVMO/ComTnwa3efCNS7+9bg8W1A/b5eaGYTgAkA9fX1NDY2HvSHt7S0FPW6cknNe4YbD3uOXak0jd0P/eSVqLe31NTe6lZr7ZXKkN84l5hnSFge2k4aVzEtUhPCKqZHuHuzmR0HzDKz1e0fdHcPCu1/EhTeEwEaGhp81KhRB/3hjY2NFPO6ctiyYg7/2/rfdLEssRVTYPyhj25Eub2dQe2tbrXWXqkMq7sOpq8niFmOmE4aF6kpocyZdvfm4N+XgAeBM4EXzawnQPDvS2HEFrYNi2eQJEuMdqMbIiISafe/dDzjc/9F5twbSzIIIiKVo+zFtJkdbmZHtN0HPgQ8C0wFxgdPGw9MKXdsYXN37n7xRHKxJGhJPBGRipDPOzOe3Ua3k8+my/u+qkJapMaEMc2jHnjQzNo+/253n2Fmi4A/mtnVQBNwWQixhWppegfTd5zIJef9jvO6rtGySiIiFWBx0ys079jFV0efEnYoIhKCshfT7r4BGLyP/duB88odT5Q8+HQzXRIxzhz5Ieh6YdjhiIhIByyZN5Nr6xoZfWQPoFfY4YhImUVpabya1prNM235Fj44qJ4juibDDkdERDogs2kBV627lrpYlvg9pTlpXEQqiy4nHhHLnnqUy/fcx/jeNXnepYhIRWpaMpOkZ4nrpHGRmqViOgrSKU5//JN8JfknGuZcBelU2BGJiEgHTNnRj6wlcJ00LlKzNM0jAlqef5yueS32LyJSSVr2ZPnNpmPoNvBmPnviFp00LlKjVExHwKw3BjCGBHHLYRrZEBGpCNOXb2V3Jk/DiNFwUo+wwxGRkKiYDpm784s1PVjU4ya+P3SHRjZERCrE5EUvcPJx3Rh64lFhhyIiIVIxHbLFTa+w4a+v87lLPwQNvcMOR0REOuCFZxoZ1nw3A4d/mOC6CSJSo1RMh+zeRWkOr4tz/rt7hh2KiIh0RDpFz4cu4yuJDPFlU2BwLx1RFKlhWs0jRDt3Z3h4+VYuHHw8h3fR3zUiIpUgu+EJzAsnjVsuo+XwRGqciukQPbx8K7syOS57r6Z3iIhUigW5d5HxBHkthyciaJpHqJbNf5RvdF/GEI4Bzgo7HBER6YCb1/bgqC7f4ZcjdkFfnTQuUutUTIekaVkj3/rbjXSxLPb7+3QJWhGRCrByy6ukNv6Nr48dQ+xf+oUdjohEgKZ5hGTNwukkyRLTJWhFRCrGpKc28bZknMu0+pKIBFRMh6BlT5Y7tvQmH0uC5tyJiFSE7S17eGjZFi45oxfdD0uGHY6IRISmeYTgoaXNzNvTj42X3s0pu57RhVpEpCKZWW/g90A94MBEd/95uFF1nscfe5jP+Ayu6HdF2KGISISomC4zd+fOBU2cevyRvPOMEWAfCDskEZFiZYHr3f1pMzsCWGJms9z9ubADK7Vs0wIuWPZZksks8alT4Sid5yIiBZrmUWZPv/AKq7ft5MphJ+mqWSJS0dx9q7s/HdzfCawCeoUbVedYNX86Cc8S13kuIrIXjUyXUzrF1ml/ZESXkxh3+uiwoxERKRkz6wMMARbutX8CMAGgvr6exsbGg37vlpaWol5XKnl37l7Tg1ssQYwsbnGe+dvhvNZJMYXd3nJTe6tbLbRXxXS5pFP4pAsZk2lldCxJ8sWzdIhQRKqCmXUD7ge+5O6vtX/M3ScCEwEaGhp81KhRB/3+jY2NFPO6Upm+YiuPvfEGS0ZP4pzEKugzkqGdmL/Dbm+5qb3VrRbaq2K6XDbNxbOtJCyPky0cIlQxLSIVzsySFArpu9z9gbDjKTV359bH19HvmMMZdu65EBsbdkgiEjGaM10m+RNH0EqCHDFMS+GJSBWwwokfvwNWuftPwo6nMyx96lHOffEPfH3wTuIxneciIv9MI9Nl8sTuvty852t8b8irvGv4WI1Ki0g1OAf4BLDCzJYF+77m7tNDjKlk/IWFnDrrSt6TzBJfOBUGagUPEflnKqbL5M4FTbxw+Gn0v/g8SOiAgIhUPnd/Eqja4dq1qRn08ywJa7eCh4ppEdmLqroyaNr+OrNXv8QVZ51EnQppEZHIy+WdWza8g6wlcF2pVkT2QyPTnS2dYt0j99MQ78mVZ50XdjQiItIBDy1tZurfTuCj7Vbw0Ki0iOyLiunOlE7hky7i3MweRiaT1L16NhypZCwiEmW7Mzl++tgaTut1JMPPHaEVPERkvzTnoDNtmotn95CwPMm25fBERCTSpj38EBe9NpnvN+wiphU8ROQANDLdidqWw0uQJaH5diIikffXVXM5f+nnqEtmic+eCidoBQ8R2T+NTHeiObv7csWer7Fm0LUwXglZRCTq5s+eQpIscdqt4CEish8ame5Et8/bRPMR72bAJe+HuP5uERGJsoUbtnN7cy/GvC0JntUKHiLSISqmO8nm5Y2ctv5uLmwYTVKFtIhIpO3ZOJ/lk+/i2CMGkfvYFJLNT2kFDxHpEBXTnSGd4rgHL+MriQzxlVOgobcSsohIVKVTxP4wjk/lMnw6UUc8+WcYeX3YUYlIhdCQaSfYvXYOsXyGhOWxXEZz7kREImzb8llYrpCz43nlbBE5OCqmO8Gjb5xMBl01S0Qk6vZsmM+iZcvJWVw5W0SKomkeJZbLOz9c2Z2FR9/E94bs0Jw7EZGoWnwHiWlfYaznIZ7Eho6HwZcrZ4vIQVExXWKLn5zBRa89yMizPgIjrw47HBER2Zd0ivy0LxPzPGYUVu/ofoIKaRE5aCqmSymd4vS/jOeMZIb4k1Ohv9aWFhGJotyj3yRGoZB2wEDTO0SkKJozXULb591B0veQII9psX8RkWhKp7D0/EIVTVBIHz1Agx8iUhQV06WSTtF99WSMID/HEhrlEBGJoMz6J8hDYXpHm2FfCCscEalwKqZLZPcj3yDuOcyCUY4TGjTKISISQQu2Oe7xwsCHxeCcL0HDVSFHJSKVSnOmSyGdom7Lgn/cl90dTiwiIvKW9mycz3tX/5CY5bFYAsb+WIW0iBwSjUyXwJ51c8j7XocMh3wytHhERGTfVsx7mIRniePgDru2hx2SiFQ4FdMlMH1nfzLU4ZgOGYqIRNTO3Rlu3fgOcpYEXaBFREpE0zwO0e4N82leOot7enyBTw89UhdpERGJonSKpbMe5NU3juOFS+5hwK5lytciUhIqpg9FOkXiznF8Lp8h9lod9PmzErOISNSkU+QnXcjZmVaGdU1SVz8Nel8fdlQiUiU0zeMQZNY/AfkMCcsTy2e0rrSISBRtmotnW0lYniRZ5WoRKSkV04fgkZb+ZDyBa+6diEhkPVc3mFZPkCOOKVeLSIlpmkeRWtY+xcbFj3J3j89z9dDumnsnIhJBrdk8181L0rfLd7jl7NeJ9z9XuVpESkrFdDHSKbrc/RG+6JorLSISWekUS/7yEEe8fCwf/eQV1L2rPuyIRKQKqZguwo7n/kK3YK40bXOlVUyLiERHcNLhezOtTO6apK7bOYCKaREpPc2ZLsJv0seTQXOlRUSiKrvhCZ10KCJloZHpg/T46pe4dd3R9DnrV/zrMZs0V1pEJILuffkkLvYEXSxHTIMeItKJVEwfhN0b5vP8fXdyYY/3MO7Cf4OEBvZFRCIlnaJpyUweeLobrwz8Gdf03aZBDxHpVCqmOyqdIv6HcXwmn8H8PuJbhyg5i4hESTBPule2lbvqEuTPmQr9Px52VCJS5TS02kFNS2ZiwUmHcV2gRUQkct6cJ02eLpbjsC1PhR2SiNQAFdMdsOONVr777NFkTScdiohEkbtz64Z30OoJ8ro4i4iUkaZ5HECuaSGzH5zMq7tPovnie+n/+lLNvxMRiZJ0ijkzH2DO+np6nvkrLtPJ4SJSRiqm9yeYfzcul2FcXZLEsdNg6PVhRyUiIm3SKbK3X8CIXIazuyZJnvFnOPGSsKMSkRqiaR778czcaVgumCftWqdURCRqFjVOhSBPJ8liTU+GHZKI1BgV0/uSTrHy3m9xz8oWcrEkbpp/JyISKekUCyZ9nftXvUFeeVpEQqRpHnsLDhmeksvwnWQSG/M/WOsOzb8TEYmIbNMC8ndcREM+w9AuSWJjb8J2v6I8LSKhUDHdJp0is/4JlixfQUPb1A7LFgrpkZonLSISBa+ueZKm+7/JqflW4uY42UIhrTwtIiGJXDFtZmOAnwNx4LfuflOnfVg6BZvm4m/rQf6RG4jlWhniMYjFcUyHDEVEoiKdYsuc2zh63X0M8iwxc7CY8rSIhC5SxbSZxYFbgQ8Cm4FFZjbV3Z8r2YcsvoOhS26BjcfhmxfhuSx5DPM8cXMsZsTO+AR0761DhiIiYUuneOcz3yfXuIh3eB4zMAOIQb9RMOpG5WkRCVWkimngTGCdu28AMLPJwDigNMX04jvwaddxBOA714JDzCgU07E4jhOL18HgK5ScRUTClk6Rve18euZbgUIR7YV7kOiiQlpEIiFqxXQvIN1uezNwVvsnmNkEYAJAfX09jY2NHX7z9zxzO0cBFmy7UZjOEUuy9uTPkMzsZMfbT+O19W/A+o6/b9S1tLQc1M+p0qm91a3W2lvTNs0l5plgJLrALA5njIfBl6uQFpFIiFoxfUDuPhGYCNDQ0OCjRo3q+Iu7fQqmXYdTKKjbknJ88OWcUsVJubGxkYP6OVU4tbe61Vp7a1qfkcTiSTzXWhgEsTic/xNouCrkwERE/i5qxXQz0Lvd9gnBvtIIEvDOObdwZK9T4JzrNLIhIhJVvc+Eqx5my/Qf06vX8RqNFpFIiloxvQgYYGZ9KRTRHwOuKOknNFzF0y19NLIlIlIJep/J2lM+Ty/lbBGJqEgV0+6eNbNrgJkUlsa7zd1XhhyWiIiIiMg+RaqYBnD36cD0sOMQERERETmQWNgBiIiIiIhUKhXTIiIiIiJFUjEtIiIiIlIkFdMiIlI0MxtjZs+b2TozuyHseEREyk3FtIiIFMXM4sCtwIeBQcDlZjYo3KhERMpLxbSIiBTrTGCdu29w91ZgMjAu5JhERMoqckvjiYhIxegFpNttbwbOav8EM5sATACor6+nsbHxoD+kpaWlqNdVKrW3uqm91UfFtIiIdBp3nwhMBGhoaPBirj7b2NhYU1etVXurm9pbfTTNQ0REitUM9G63fUKwT0SkZpi7hx1D0czsZaCpiJceA/y1xOFEmdpb3dTeynWSux8bdhDFMrMEsAY4j0IRvQi4wt1XvsXzlbM7Ru2tbmpv5dpnzq7oaR7F/idkZovdvaHU8USV2lvd1F4Ji7tnzewaYCYQB257q0I6eL5ydgeovdVN7a0+FV1Mi4hIuNx9OjA97DhERMKiOdMiIiIiIkWq1WJ6YtgBlJnaW93UXql2tdbnam91U3urTEWfgCgiIiIiEqZaHZkWERERETlkNVdMm9kYM3vezNaZ2Q1hx3OozKy3mT1uZs+Z2Uozuy7Y38PMZpnZ2uDfo4L9ZmY3B+1fbmZDw21BccwsbmZLzWxasN3XzBYG7brXzOqC/V2C7XXB433CjLsYZvZ2M7vPzFab2SozG17N/WtmXw5+l581s3vMrGs196/sn3J25X+nQTm7mvtXObvGimkziwO3Ah8GBgGXm9mgcKM6ZFngencfBAwDvhi06QZgtrsPAGYH21Bo+4DgNgH4ZflDLonrgFXttn8A/NTdTwZeAa4O9l8NvBLs/2nwvErzc2CGuw8EBlNod1X2r5n1Aq4FGtz9NArLrX2M6u5feQvK2UCFf6fbUc6uwv5Vzg64e83cgOHAzHbbNwI3hh1Xids4Bfgg8DzQM9jXE3g+uP9r4PJ2z3/zeZVyo3CVtdnA+4FpgFFYED6xdz9TWP92eHA/ETzPwm7DQbS1O7Bx75irtX+BXkAa6BH01zRgdLX2r24H/H1Qzq7w73QQs3J2lfavcnbhVlMj0/y909tsDvZVheBwyRBgIVDv7luDh7YB9cH9avgZ/Az4DyAfbB8N7HD3bLDdvk1vtjd4/NXg+ZWiL/AycHtwiPS3ZnY4Vdq/7t4M/Ah4AdhKob+WUL39K/tX0b/PB6KcXZXfaeXsGszZtVZMVy0z6wbcD3zJ3V9r/5gX/gSsimVbzOwC4CV3XxJ2LGWSAIYCv3T3IcDr/P3wIFB1/XsUMI7Cf0jHA4cDY0INSqQTKGdXLeXsGszZtVZMNwO9222fEOyraGaWpJCU73L3B4LdL5pZz+DxnsBLwf5K/xmcA1xkZpuAyRQOG/4ceLuZtV3Rs32b3mxv8Hh3YHs5Az5Em4HN7r4w2L6PQqKu1v79ALDR3V929wzwAIU+r9b+lf2r9N/nfVLOVs6mevpXOZvaK6YXAQOCs0zrKEySnxpyTIfEzAz4HbDK3X/S7qGpwPjg/ngK8/La9n8yOIN4GPBqu0NPkefuN7r7Ce7eh0L//cXdPw48DlwaPG3v9rb9HC4Nnl8xIwLuvg1Im9kpwa7zgOeo0v6lcKhwmJkdFvxut7W3KvtXDkg5u8K/08rZytlUUf++pbAnbZf7BowF1gDrga+HHU8J2jOCwuGi5cCy4DaWwhyk2cBa4DGgR/B8o3B2/HpgBYUzcENvR5FtHwVMC+73A1LAOuBPQJdgf9dge13weL+w4y6inacDi4M+fgg4qpr7F/gOsBp4FvgD0KWa+1e3A/4+KGdX+He6XduVs6uwf5WzXVdAFBEREREpVq1N8xARERERKRkV0yIiIiIiRVIxLSIiIiJSJBXTIiIiIiJFUjEtIiIiIlIkFdNSdczsaDNbFty2mVlzcL/FzP4v7PhEROQfKW9LJdPSeFLVzOzbQIu7/yjsWERE5MCUt6XSaGRaaoaZjTKzacH9b5vZJDOba2ZNZnaxmf3QzFaY2Yzgcr+Y2RlmNsfMlpjZzLbLwe7nM85tN7qy1MyOKEfbRESqkfK2VAIV01LL+gPvBy4C7gQed/d3A7uA84PE/AvgUnc/A7gN+N4B3vPfgS+6++nAyOC9RESkNJS3JXISYQcgEqJH3D1jZiuAODAj2L8C6AOcApwGzDIzgudsPcB7zgN+YmZ3AQ+4++bOCFxEpEYpb0vkqJiWWrYHwN3zZpbxv59AkKfw3TBgpbsP7+gbuvtNZvYwMBaYZ2aj3X11qQMXEalRytsSOZrmIfLWngeONbPhAGaWNLNTg/vXmNk1e7/AzPq7+wp3/wGwCBhY1ohFRGqb8raUnYppkbfg7q3ApcAPzOwZYBlwdvDwQGD7Pl72JTN71syWAxngkbIEKyIiytsSCi2NJ1KE4Ozyi4PELSIiEae8LZ1FxbSIiIiISJE0zUNEREREpEgqpkVEREREiqRiWkRERESkSCqmRURERESKpGJaRERERKRIKqZFRERERIqkYlpEREREpEj/D62Zn5i/zJesAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mp.plt.rcParams[\"figure.figsize\"] = (12,6)\n",
    "# Compute and plot altitude, velocity\n",
    "r = 1e-3 * (np.sqrt(x[:, 0] ** 2 + x[:, 1] ** 2 + x[:, 2] ** 2) - Re)\n",
    "v = 1e-3 * np.sqrt(x[:, 3] ** 2 + x[:, 4] ** 2 + x[:, 5] ** 2)\n",
    "y = np.column_stack((r, v))\n",
    "fig, axs = post.plot_single_variable(y, t, [[0], [1]], axis=0)\n",
    "\n",
    "x, u, t, _ = post.get_data(interpolate=False)\n",
    "r = 1e-3 * (np.sqrt(x[:, 0] ** 2 + x[:, 1] ** 2 + x[:, 2] ** 2) - Re)\n",
    "v = 1e-3 * np.sqrt(x[:, 3] ** 2 + x[:, 4] ** 2 + x[:, 5] ** 2)\n",
    "y = np.column_stack((r, v))\n",
    "fig, axs = post.plot_single_variable(\n",
    "    y, t, [[0], [1]], axis=0, fig=fig, axs=axs, tics=[\".\"] * 15\n",
    ")\n",
    "axs[0].set(ylabel=\"Altitude, km\", xlabel=\"Time, s\")\n",
    "axs[1].set(ylabel=\"Velocity, km/s\", xlabel=\"Time, s\")\n",
    "\n",
    "# Plot mass at the collocation nodes\n",
    "figm, axsm = post.plot_single_variable(\n",
    "    x * 1e-3, t, [[-1]], axis=0, fig=figm, axs=axsm, tics=[\".\"] * 15\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}