fig2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figure 2\n",
    "\n",
    "> Evolution of the Sun's mass $M(t)$ and radius $R(t)$ (both in solid red) and Earth's semi-major axis $a(t)$ for splitting intervals $dt_\\textrm{split}$ of $400~\\textrm{yr}$ (solid yellow) and $4000~\\textrm{yr}$ (dotted blue).  The simulation starts approximately $4~\\textrm{Myr}$ before the TRGB phase.  Earth's orbital radius starts at $1~\\textrm{au}$.\n",
    "\n",
    "For more information and additional explanation of the code in this notebook, see the official REBOUNDx effect example documentation for [\"Parameter Interpolation\"](https://github.com/dtamayo/reboundx/blob/master/ipython_examples/ParameterInterpolation.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load\n",
    "\n",
    "We first load in the extracted data, exported from MESA log files, as done in [mesa2txt.ipynb](https://github.com/sabaronett/REBOUNDxPaper/blob/master/mesa2txt.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "data = np.loadtxt('m.txt') # return (N, 2) array\n",
    "mtimes = data[:, 0]        # return only 1st col\n",
    "masses = data[:, 1]        # return only 2nd col\n",
    "\n",
    "data = np.loadtxt('r.txt')\n",
    "rtimes = data[:, 0]\n",
    "Rsuns = data[:, 1]         # data in Rsun units\n",
    "\n",
    "# convert Rsun to AU\n",
    "radii = np.zeros(Rsuns.size)\n",
    "for i, r in enumerate(Rsuns):\n",
    "    radii[i] = r * 0.00465047"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "Next we set up the Sun-Earth system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import rebound\n",
    "import reboundx\n",
    "\n",
    "M0 = 0.8645388227818771 # initial mass of star\n",
    "R0 = 0.3833838293200158 # initial radius of star\n",
    "\n",
    "def makesim():\n",
    "    sim = rebound.Simulation()\n",
    "    sim.G = 4*np.pi**2               # use units of AU, yrs and solar masses\n",
    "    sim.add(m=M0, r=R0, hash='Star')\n",
    "    sim.add(a=1., hash='Earth')\n",
    "    sim.collision = 'direct'         # check if RGB Sun engulfs Earth\n",
    "    sim.integrator = 'whfast'\n",
    "    sim.dt = 0.1*sim.particles[1].P\n",
    "    sim.move_to_com()\n",
    "    return sim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can create an `Interpolator` object for each parameter set and pass the corresponding arrays as arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = makesim()\n",
    "ps = sim.particles\n",
    "rebx = reboundx.Extras(sim)\n",
    "starmass = reboundx.Interpolator(rebx, mtimes, masses, 'spline')\n",
    "starradius = reboundx.Interpolator(rebx, rtimes, radii, 'spline')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Runs\n",
    "\n",
    "We integrate twice for 4 Myr each, updating the primary's mass and radius interpolated at the time between outputs, i.e., $dt_\\textrm{split}$. We interpolate and update 1000 times ($dt_\\textrm{split} = 4000$ yr) in our first run, and then for 10,000 times ($dt_\\textrm{split} = 400$ yr) in our second."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 18.3 s, sys: 1.36 ms, total: 18.3 s\n",
      "Wall time: 18.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "Nout = 1000\n",
    "a1 = np.zeros(Nout)\n",
    "ts1 = np.linspace(0., 4.e6, Nout)\n",
    "T0 = 1.23895e10 # Sun's age at simulation start\n",
    "\n",
    "for i, time in enumerate(ts1):\n",
    "    sim.integrate(time)\n",
    "    ps[0].m = starmass.interpolate(rebx, t=T0+sim.t)\n",
    "    ps[0].r = starradius.interpolate(rebx, t=T0+sim.t)\n",
    "    sim.move_to_com() # lost mass had momentum, so need to move back to COM frame\n",
    "    a1[i] = sim.particles[1].a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 17.8 s, sys: 20 ms, total: 17.9 s\n",
      "Wall time: 17.9 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# reset sim\n",
    "sim = makesim()\n",
    "ps = sim.particles\n",
    "rebx = reboundx.Extras(sim)\n",
    "starmass = reboundx.Interpolator(rebx, mtimes, masses, 'spline')\n",
    "starradius = reboundx.Interpolator(rebx, rtimes, radii, 'spline')\n",
    "\n",
    "Nout = 10000\n",
    "mass = np.zeros(Nout)\n",
    "radius = np.zeros(Nout)\n",
    "a2 = np.zeros(Nout)\n",
    "ts2 = np.linspace(0., 4.e6, Nout)\n",
    "\n",
    "for i, time in enumerate(ts2):\n",
    "    sim.integrate(time)\n",
    "    ps[0].m = starmass.interpolate(rebx, t=T0+sim.t)\n",
    "    ps[0].r = starradius.interpolate(rebx, t=T0+sim.t)\n",
    "    sim.move_to_com() # lost mass had momentum, so need to move back to COM frame\n",
    "    mass[i] = sim.particles[0].m\n",
    "    radius[i] = sim.particles[0].r\n",
    "    a2[i] = sim.particles[1].a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot\n",
    "\n",
    "We plot the Sun's mass over time, and a comparison of the Sun's radius and Earth's semi-major axis over time, adjacent to one another."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAFACAYAAACm+Ov/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hUVfrA8e+ZmUx6I5UWAoQeJEBApAdQkJWiq7uKgAgu6wrYd3XXte+Krqyuu+hPsaEriiwWRKlSRCnSew2EEiAJ6b3MzPn9MSEmIaQMSSbl/TzPPMnce+697wxh3jn3NKW1RgghhCjP4OwAhBBCNEySIIQQQlRIEoQQQogKSYIQQghRIUkQQgghKiQJQgghRIVMzg6gtgQGBurw8HCHjs3JycHT07N2A6oFElfNNdTYJK6akbhq5lri2rVrV7LWOqjCnVrrJvHo27evdtSGDRscPrYuSVw111Bjk7hqRuKqmWuJC9ipr/K52mRqEI7K3bUL33ffI2HTjxj9/Ysffpj8/TG2aFGyzWA2OztUIYSoV80+QVjT0zHFnyMzNhZrRsZVyxk8PH5JIC387QnEz//KbZcfvr4oo7EeX4kQQtSuZp8gvEeOJMVoZPjw4WiLBWtGBta0NKxpaVjS0rCmpWNNSy37PDWNwtiTWNLT0bm5FZ9YKYy+viUJwxQQgDEwAFNAIKaAFhgDAjAFBtq3BwRi8PRAKVW/L14IISrR7BNEacpkwhQQgCkgoNrH2PLzsaan/5JUUtNKfremFyeVlFQKTp3CumMH1vT0iq/t5oapRQuMxUnDFBiAZ1Y2qWfPYQq0JxFTYID9tpefnyQTIUSda/QJQik1DhgXERHhlOsb3NwwhIbiEhparfK6qMieRFKSsaSkYklJxpqSgiU5pfj3VIouXiTv4AE8U1JJXLXqypOYTJhatLDXQIKDMQUFVfwzoAXK1Oj/iYUQTtLoPz201suB5dHR0b9zdizVoVxccAkJxiUkuMqyG9evZ3Dv3liTk7EUJxFr6i/JxJKcTFFiInkHD2JNSYHyM/MaDBgDWuASVEHyCA7CFBQsiUQIcVXyqdCQGQyY/O2N366dOlVa1F4zScWSlGR/XLpU8rMoKYmiS0nkHTpUcSJRyt4+EhSES0goLi1DMYW2xCU0BFNoKC4tW2IKCZGeXEI0M9eUIJRSYYAPcFhrbaudkIQj7DWTEFxCQiotVyaRXE4iSZewXEqiKDGRogsXyN29G1sFPbqMAQG4hIbiazKR8NPm4gTSEpeW9ltspuBglItLXb1EIUQ9cyhBKKVcgPeBU0AK8Del1ONa69jaDE7UvuomEltuLkUJiVgSLlJ0MYGihItYEhIpSkjAdDKWjK+/xpadXe7kyt4u0rIlLqGhuLRuXfxoZf/ZqjVGr4Y3ClUIUTFHaxAPAf8rvv+PUupz4B3g1toKTDiXwcMD1w7tce3Q/op9GzduZPjw4Vizs7EkJFB0MQFLYsIvieRiAgWxsWRv2oTOzy9zrNHPr1TiKJtAzK1bY2iA0xgI0Vw5miD6AW8ppdyBQ1rrDkop+Z/dzBi9vDBGROB6lR5kWmusqakUxcdTdP48hefPU3T+PEXnL9gTyA8/oAsKyp6zfAJp0xpzWDvM7cJwadVKGtOFqEeO/m/bAYzSWn8DdFBKhQA51T1YKTUGeAMwAu9prV8utz8M+AjwKy7zpNZ6hYOxCidRSpWMK3Hv1euK/VprrCkpxUmjXAI5cYLsjRvRhYW/HGAy4dK6lT1hhIXZk0ZYGOawdri0aS2N6ELUMkcTxBvAAqVUP+xtEMOAx6pzoFLKCLwJ3AjEAzuUUt9orQ+XKvZXYInW+v+UUt2BFUC4g7GKBkpdbrMIDKw4gdhsWC4lU3T2DIVnz1F49iyFZ89QdOYsGbt3Y8sp9Z3EYMAlNBQ/H28ubthYNoG0a4fB1bUeX5kQTYNDCUJrXQTcq5RqA3gD/65BL6b+QKzW+hSAUmoxMAEonSA09t5RAL7ABUfiFI2bMhhKxox49OtXZp/WGmtaGoVnzlB09iyFZ85SeO4cOQcPkLV6ddkR6waDvY2jQ3tcw9tj7tABc/twXDt0wBgQIKPShbiKKhOEUioZ2AfsAfYW/zyitbZpreMduGZr4Fyp5/HA9eXKPAesUUrNATyBUQ5cRzRhSin7aPIWLaB375LtJzZupNfw4VgzMuy1jjNnKIyLozDuFAWn4sj9eXuZhnODj489WbTvgLl9e3sS6dABc9u2KLllJZo5pcsPmipfQKlc4CywGegB9AQUcBB7wtirtX6r2hdU6g5gtNb6vuLnU4D+Wus5pco8WhzbP5VSN2DvUhtZvpailJoJzAQICQnpu3jx4uqGUUZ2djZeXl4OHVuXJK6aqzI2mw1DWhqmhESMiQnFPxMxJSRgLDX2QxsMWAMDsYaGYGnZCkurllhatcISGgoOjPVoqO+ZxFUzTTGumJiYXVrr6Ir2VecWUydgLjAUeBL4CugC9AaigIlAtRME9hpD21LP23DlLaQZwBgArfVWpZQbEAgklS6ktV4ALACIjo7Ww4cPr0EYv7jcbbOhkbhq7lpis2ZnUxh32l7biIuj8FQchadOUrBuHVgs9kIGA+awMFw7ReDaqROuERGYIyJwDQ+vtMbRUN8ziatmmltcVSYIrfV5YKpSqg/wKvAI8IjW+lPgUweuuQPopJRqD5wH7gQmlStzFhgJLFRKdQPcgEsOXEuIajN6eeHeMxL3npFltuvCQgrPnKEgNpaCE7EUnDhBQWwsWevWg624UmsyYQ5vV5I0XDt1wrVTJ8xhYbIuiGi0qt1IrbXeDYwsnj31Q6XUPmCO1jq1JhfUWluUUrOB1di7sH6gtT6klHoB+9J332DvEfWuUuoR7A3W03RV98KEqCPKbC75wOfmX7bbCgoojIuzJ4wTsRTExpJ/8BBZq1aXzHel3N1x69wZb18f0hIScevWFdfOnTG4uzvp1QhRfdVKEEopL6A7EIm9HSIJuAv76OlNNb1o8ZiGFeW2PVPq98PAoJqeV4j6ZHB1xa1rV9y6di2z3ZabS8GpOAqOHSP/2FEKjhzFbcdOEjb9WHygAXO7dvZk0bUbbt3s5zAFVbxuvBDOUp1eTHGAP3AI2A8cwN4OcUBrffU1OoVopgweHrhH9sA9skfJto0bNjCwU2cKjh4h/+gx8o8eIW/ffjJXrCwpYwwMxK1bN9wie+DesydukZG4BFc9LbwQdaU6NYh22NsKUrDXHC4C5yU5CFEDSmFu0xpzm9Z4j/ql17Y1I4P8Y8coOHqU/CNHyT90iJTNm0vaNkzBwbj17Il7z0jcekTiFtkDk7+/s16FaGaqkyBaYO+tFIW959JEoFtx99d92Lu5PlR3IQrRdBl9ffHs3x/P/v1Lttny8sg/coT8gwfJO3CQ/AMHyF63rmS/S9u29lpGZE/cr7PXNKRNQ9SF6vRiSgc2Fj8AUEqZsbdHXO7qKoSoJQZ3dzz69MGjT5+SbdasLPIPHfolaezbT9bK4uVoTSbcunXDvXcUHlFRuPfujUvLlk6KXjQljk61UQjsLn4IIeqY0dsbzwED8BwwoGSbJTWVvH37yNuzl7w9e0hf8j/SPv4vAKbQUNyjovDobU8Ybl27yshwUWMyd7IQjZSpRQu8Y2LwjokB7KsF5h87Tt6ePfbH3r1krbLXMpSrK249I/Ho1w/Pfv1wj4rC4OHhzPBFI9DoE0TxuIxxEVdZk0CI5kK5uPzSe2rKZACKEhPtNYy9e8ndtYuUBe+S8n9vg8mEe2QkXqEhZBtNePTpLYs1iSs0+gRRvKrd8ujo6N85OxYhGhqXkBBcxozGZ8xoAKzZOeTt2U3u9u3kbt+Bx5q1nFu1GoxG3Hr0wKNfNB79+uERHY2xAc45JOpXdcZBnANWYh/YtlZrXe2FgYQQDYvRyxOvIUPwGjIEgB9WraKvtze5O3aQu30HqR//l9T3PwCjEfdevfAcOBDPgQNxv66nrObXDFXnX7w/MBaYgn36i73Yk8V3WuvjdRmcEKJuaTc3vAYNwmuQfeICW14eeXv3krN1GzlbtpD85pskz5+PwcsLj+uvx3PgDXgNGoRLu3ayjkYzUJ1urhexT7f9vlLKhH1W17HA18XdXS9Pm7FBa11w9TMJIRo6g7s7njfcgOcNN8Cjj2BJSyP355/J2byFnC1byF63jkTApVUrPAcNxHPQYDwHD5LbUU1UjeqMWmsLsL748bhSKhz4FTAH+7iIebUcnxDCiUz+/viMGYPPmDForSk6e9aeKDZvJnPlKtL/txRMJjyio/EaNgyv4cNwbd/e2WGLWnJNNxW11qexry/9Zq1EI4RosJRSmNu1w9yuHf533YW2WMjbu5fsH34ge+NGkl55haRXXsGlXZg9WQwbhke/fhhk/EWjJa1OQgiHqOKag0d0NMGPPUZh/HmyN9mTRfriz0n7+L8YPDzwHDQI7xtH4TV8OEYfn6pPLBoMSRBCiFphbtOaFpMm0WLSJGx5eeRs22avXazfQNbateDigueAAXjfOArvkSMxBQQ4O2RRBUkQQohaZ3B3LxnlrZ95hvz9+8lcu5asNWtJeOZZEp59Do++ffG+6Ua8R43CpVUrZ4csKtDoE4SMpBaiYVMGA+5RUbhHRRH8+OMUHD9O1pq1ZK1ZQ+JLc0l8aS5uPXviM3YsPmNvxiUkxNkhi2LVXVEuBLgJ6AX4AenYp/peq7VOqLvwqiYjqYVoPJRSuHXpgluXLgTNmU1BXBxZa78na9UqeyP3P/6BR79+uHeKwNKrl6x94WSGynYqpboppZYCh7EPlHMBEop/TgEOKaWWKqW613mkQogmx7V9ewJn/o72X35BhxUrCHzgASxJSfgs+pQTQ4Zy7v4/kLH8W2w5MoGDM1RVg1gIvArcDRRprW2ldxYPlJuAfSDdDXURoBCieXDt0J6gObMJnD2Lzf/9L50SEslcsYLsjRtR7u743HQjvrfehkf/fihDpd9tRS2pNEFora+//LtS6hIQVG5/IfC/4ocQQlwzpRSWsDBCpk4l+PHHyNu1i4zl35K5YgUZy77BpU0bfCdOxHfiRMxtWjs73CatyjSslLpOKeUCuF1l/9laj0oIIbA3cHv060fLF56n04+baPXqq5jD2pL85pucHDWKM/dMI2PZMmy5uc4OtUmqTiP1Cuw1B4NS6jNgL8VrUWNPML51F54QQtgZ3N3xHXcLvuNuoejCBTKWLSP9q6+58MSTGF54EZ/x4/C/8y7cunR2dqhNRpU1CK11G6A1UAT8CHQAngdOAueAT+oyQCGEKM+lVSsC//AHOq5eRbtP/ov3qFFkfPElcRMmcHrS3WQsX46tQOYOvVbV6uaqtU5WSvXUWp+8vE3Z5/p111pL3U4I4RRKqV+m+3jyCTK++pq0zxdz4Y9/wvjSXHxvuxX/3/4Wc1iYs0NtlKrdFaB0cih+riU5CCEaCpO/PwHT76XjypWEffA+HtHRpC78iJM3jebs739PzpYtaK2dHWajUtU4iAeVUq5VlHFVSj1Yu2EJIYRjlMGA58CBtPnPv4lYv57A2bPJP3iIs9NnEDdhIulffCG3n6qpqhpEKBCrlHpHKTVJKdVXKdW5+OddSql3gBNAcN2HWjGl1Dil1IKMjAxnhSCEaKBcQoIJmj2LiA3rafnSSwBcfOqvxI4YyaX5b2JJSXFyhA1bpQlCa/0XoA/2JDAD+9rUB7H3bJoOHAV6a63/WsdxVhbjcq31TF9f6UwlhKiYwWzG77Zbab/sa8I+/AD3yEiS588nNmYEF59+hsIzZ5wdYoNUnSVHL2FfKU5WixNCNGpKqZIlVQtOnSL1o4/J+Oor0r/4Ap+bbyZg5kzpJltKdQbKJSul1iml5imlJiuluiulZJy7EKJRc+3QgZbPP0fH79fS4t5pZG/YQNyECZx7YBZ5+/Y5O7wGoTof9B7Yx0H4A7OBHUC2Umq7UmqBUuqBugxQCCHqkktwMCF//CMR69cROHs2ubt2cfq3d3Lm3nvJ3bHD2eE5VXUSRCdgOzAU+8R93kBf4F/Yp/2eWGfRCSFEPTH6+dkbtNetI/hPf6IgNpYzU6Zydvp08vbudXZ4TlGdkdTntdZTgd8CDwCbAC+t9ada6z9prW+q6yCFEKK+GL08CZh+LxFr1xL8xBPkHznK6Tvv4tz9f8B09pyzw6tXNRkot1trPRJ4BfhQKbVIKdXCkYsqpcYopY4ppWKVUk9epcxvlFKHlVKHlFKfOnIdIYRwlMHNjYB7pxHx/VqCHn6Y3N27CXjpJeIffIiCEyecHV69qFaCUEp5KaX6K6WmA8OBJOAuILKmF1RKGYE3gZuB7sBd5RccUkp1Av4MDNJa9wAerul1hBCiNhg8PQm8//dEfL+W7F+NJWfzZk5NmMiFv/6VosQkZ4dXp6rTiykOiAdeB/phn6TvGcBfa73JgWv2B2K11qeK15NYjH3RodJ+B7yptU4D0Fo37X8FIUSDZ/TxIWfcOHuvpylTyFj2DSfHjOHSv//TZFe8q04Noh2QBaRgrzlcBM5rrR0dutwa+yywl8UXbyutM9BZKbVZKbVNKTXGwWsJIUStMvn7E/LnJ+m44ju8hg8j+a23iB09hrTPl6AtFmeHV6tUVZNXKaX8gKjiR+/in92AXIrXhdBaP1TtCyp1BzBaa31f8fMpQH+t9ZxSZb7FPr34b4A22KcZj9Rap5c710xgJkBISEjfxYsXVzeMMrKzs/Hy8nLo2LokcdVcQ41N4qqZxhSXy6lTeH3xJeaTJ7G0bEnWbbdRGNkDlHJqXNUVExOzS2sdXeFOrXWNH4AZ+xQcM4D/1PDYG4DVpZ7/GfhzuTJvA9NKPV8H9KvsvH379tWO2rBhg8PH1iWJq+YaamwSV800trhsNpvOWL1an7jpJn24S1d9+p5pOu/QIafHVR3ATn2Vz1WHRkRrrQu1vVfT+7rUN/9q2gF0Ukq1V0qZgTuBb8qV+RqIAVBKBWK/5XTKkViFEKKuKaXwuekmOi5fTshTT1Fw9Chxv76dC088SdHFi84Oz2H1PmWG1tqCfUT2auAIsERrfUgp9YJSanxxsdVAilLqMLAB+KPWWqZdFEI0aMpspsWUyXRcs5qA+2aQuXIlJ8fcTNJrr2PNznZ2eDVWrRXlapNS6gPgFiBJa13STVZr/UypYsOAEYAGUrXWjjUuCCGEExh9fAh+7DH877yTpH+9QcqCBaR+9BEGNzd7gdLtE5d/L/+zzLbLPyou4zFoEAwfXuuvo94TBLAQmA98XNHO4kbxt4AxWuuzSimnrTUhhBDXwqV1a1q/+g9aTJ1K5rfL0RYrlOkYZP9dX95W8rN0kXL7Su28fFxyYFDtB48TEoTWepNSKrySIpOAL7XWZ4vLyxgIIUSj5t4zEveeNR5XXG3HN26sk/M2xGm7OwP+SqmNSqldSqmpzg5ICCGaI2fcYqqKCftssSMBd2CrUmqb1vp4+YKlx0G4u7vTpUsXhy5os9kwGBperpS4aq6hxiZx1YzEVTPXGFefq+65Wv/XunwA4cDBq+x7Eniu1PP3gTuqOqeMg6g/DTUurRtubBJXzUhcNdOgxkHUsWXAEKWUSSnlAVyPvTusEEKIeuSMbq6fYZ8RNlApFQ88C7gAaK3f1lofUUqtAvYDNuA9rfXB+o5TCCGaO2f0YrqrGmVexb56nRBCCCdpiI3UtaawsJCTJ0+Sm5tbaTlvb2927dpVT1FVn8RVc7Udm4eHBx07dsRsNtfaOYVoLBp9glBKjQPGRUREXLHv5MmT+Pn50aVLlwbZ80A0bDabjcTERE6ePEm3bt2cHY4Q9a7Rf2pqrZdrrWf6+vpesS83N5eQkBBJDsIhBoOBkJCQKmugQjRVTf6TU5KDuBby9yOaM/nrF0IIUSFJEEIIISokCaKBy8vLY9iwYVitVuLj4/n8888Bew+toUOHYqnhGrgzZszgu+++AyhzPmewWq307t2bW265pWTbqlWr6NKlCxEREbz88stlyle2T4iGymK1oTOWwPnfwvnbuXD0Hqa8uYA/fvQ6xN9W/LiVRz98nUnzF3DhyCSIHw/nxvHxiif4zb/fZfn6h+DcWDh3Mwd3T+bX/3qXZxa9BGdvgrM30tH7rTqJXRJEA/fBBx9w2223YTQaWbduHbt37wbAbDYzcuTIGn/A7927l169egGUOZ8zvPHGG2V6B1mtVmbNmsXKlSs5fPgwn332GYcPH65yX01orbHZbLX2GoSoTH5hEb9//xNe+Pp7dOYSyPqCvIx1/HiuNTvPe0L2V8WPr9lz0ZMt8a3Jz/wBspdDzrecTUlh+4VWJKSeg5yVkLOKrIz97EpoxbEkA+Suhdzv8XY5VifxN+lurg3F0qVLmTdvHnl5eXh7e/PVV18RFFS9+dsXLVrEp59+yk8//cSjjz6Kn58fq1ev5quvvmLixIn8+c9/5u67777q8cePH2f69OlkZGQwffp0EhISaNOmTYXna9++fW295CrFx8fz3Xff8dRTT/Haa68BsH37diIiIujQoQMAd955J8uWLaN79+6V7rvs6aefJjAwkIceegiAp556ipCQEMaPH8/NN99MTEwMW7du5euvv6Zdu3b19lpFM2XL4cj+R/jx9Fg8XGK4b2hnWge0o2UgfHxnAe4ugdBqacniP/NuNZBfZCG09ZtgVoBiiodiVB9NmN908J0BKLoHwv9CFN5uvSB4FSjFyT2n6FsHL6HRJwit9XJgeXR09O8qLXhUVbrbYV11lUViYmK4/fbbAXj++edZsmQJs2bNqvK4wsJCTp06RXh4OOHh4fTr14958+YRGWmfV95qtbJjx46rHm+xWJg8eTLz58+nf//+PPDAA3Tt2hWAwYMHX3G+mkpLS8Pf37/MtiFDhpCVlXVF2Xnz5jFq1KiS5w8//DD/+Mc/ypQ9f/48bdu2LXnepk0bfv755yr3XTZjxgxuu+02HnroIWw2G4sXL2b79u1kZWVx7NgxPvzwQ956q26q4kKUYbkI8ePo7bWL+TEnad/pSVqH2yeR8ACGBl55SN8Khtq084byX2V8vaBfueOzLBtrI+oryC2merBw4UL69+9Pr169eOutt3C7vOwgcOzYMR544AEee+wxLly4UOa45ORk/Pz8ypQtPaW50WjEbDZX+IEM8OWXX9KtWzf69+8PQI8ePYiKirrq+aoybdo0wP5NHeCRRx65osyPP/7I3r17r3iUTg7ffvstwcHB9O1b9juP1lcmW1X87aqyfZeFh4cTEBDAnj17WLNmDb179yYgIACAdu3aMWDAgGq/ViEcFXdhD4d3/Rryd4FLB24a+iadwm90dlgOafQ1iGqrxjf9uvDxxx+zfft21q9fj5eXF0OHDqVHjx6APQHMnTuXwYMH07FjRx5//HEWLlxYMq2Du7s7+fn5AKSkpODr64uLi0uZ8xcUFJRJOKXt37+/zIfwrl27GF68bm358507d44XXngBX19f+vfvz8svv8zo0aM5duwYS5cuJT8/H09PTxISErBYLKxatYqjR48yb948Hn/88ZJrVKcGsXnzZr755htWrFhBfn4+mZmZTJ48mVmzZnHu3LmSY+Lj42nVqhVgrzFcbV9p9913HwsXLiQhIYHp06eXbPf09KzwPRKiNp04vYY7PkzFpB7gq1/70jbiYzDVzXKg9UFqEHXswIEDDBw4EC8vL7744gu2bNlCz549Adi0aRM2m41///vfPPfcc3h5eXHo0KGSY/39/bFareTn5xMXF3fFB2JKSgpBQUElH/IjR47k/PnzJfsDAgI4eNA+Ee6uXbv47LPPSmoQ5c939OhRzGYzDz74IJ6entx1113MnTuX4OBgUlJS2L17N3369GHPnj1ERUURGBjI5MmTyyQHqF4NYu7cucTHx3P69GkWL17MiBEj+OSTT+jXrx8nTpwgLi6OwsJCFi9ezPjx4wEq3VfarbfeyqpVq9ixYwejR4+u+T+YEI5Kf5/w3PF0b3GK60KzaRGxpFEnB2hONQgnueeee5gwYQJLly5l7NixdOjQoeTbbKdOnejVqxcDBw5k8ODBfPDBB2XuswPcdNNN/PTTTwwYMIDk5GQiIyNZsGABAwcOZMOGDYwdOxawzxsUGxtLixYtSo6dMmUKY8eOJSoqii5duuDn51fSa6hr165lznfjjTfStm1bZs+eTUhICHfccQcAGRkZBAUF8cknnzBy5Ei+++47brvtNjZv3lzSG6q2mEwm5s+fz+jRo7FarUyfPr2ktlXZvtLMZjMxMTH4+flhNBprNT4hKmKxFGBL+gvmzNdwMcCCX+fg3vKJpvH3d7WVhBrLAxgHLIiIiLhipaSdO3fWfHmlevbmm2/q0aNH61tuuUUvXbr0iv27d+/WkydPrvDYW2+9VR89elRrrfWBAwf0I4884nAcf/rTn/Rjjz2mZ82ape+88049e/Zsff/99+vvv/9ea631jBkztMVi0dOnT9dWq1UvW7ZMT506VR8+fNjha9YFq9Wqe/XqpY8fP15r56yNv6OmuBJZXWoscaWkn9OT//Mf/fi7D2nbYZPWae80iLhqgkpWlHP6B3xtPSpacrQxJIjqeP/997XFYimzraCgQH/00Ud1cr2rJaSG7tChQ7p9+/b60UcfrdXzSoKof40irtwt+ujP0brLX5bq3s9+pi9c3NQw4qqhyhKE3GJqBEo3tl5mNpuZOnVqnVzvv//9b52ct651796dU6dOOTsM0dRpDen/B4kP08WniP8bvZwu3V+gZWCYsyOrdfXeSK2U+kAplaSUqnQZUaVUP6WUVSl1e33FJoQQlSmwZPDgB/9h854FQBH4P0TMkHdp1QSTAzinF9NCYExlBZRSRuAVYHV9BCSEEFXKWcfp9CV8c6Ijf976IJbQTyHkX6Bcqj62kXLGmtSblFLhVRSbA3wB9KvzgIQQojK6EC79FVLncV93xYX8ztw78k5MftUfZNpYNbg2CKVUa+BWYARVJAil1ExgJkBISAgbN24ss9/b27tughTNTvm/rZrKzs6+5nPUBYmrcql5Z1kdd5bXhvwHD5MiNu1ORrQZzOljFzl97KKzwytRZwEJ+fIAACAASURBVO/X1Vqv6/IBhAMHr7Lvf8CA4t8XArdX55xNtRdTbm6uHjp0qLZYLPrcuXN68eLFWmt7L6YhQ4booqKiGp1v+vTp+ttvv9Va6zLncwaLxaKjoqL0r371q5JtK1eu1J07d9YdO3bUc+fOLVP+avsqO6Y2SC+m+uf0uKz52pb4rL513jzd7olv9UuLHtI6Z7Pz47qKuurF1BBHUkcDi5VSp4HbgbeUUhOdG5LzyHTflU/3LdOAi9qmc7fC6T6o1Od5ccCbTOp5mlkTngGPgc4Ord41uAShtW6vtQ7XWocDS4EHtNZfOzmsa7J06VIGDBhAr169GDx4MJcuXar2sYsWLWLChAkl03MvXbqUqKgo4uLimDhxIosWLar0+OPHjzN48GB69uzJ66+/fsV036XPV58uT/d93333lWwrPaW32WwumdK7sn2VHVPa008/zRtvvFHy/KmnnuLf//433bp144EHHqBPnz5l5noSzU9WThovLnmLl79cAIWHwSWCHr0/5KW7Z+Hj2aLqEzRBzujm+hmwFeiilIpXSs1QSt2vlLq/rq8d/uR3hD/5XZltMxbuIPzJ7/j+cGLJtk9/Pkv4k9/x5y/3l2xLzMwn/Mnv6P/372t83ZiYGLZt28a+ffu48cYbWbJkSbWOKz3d9+XpuZctW8bevXtp3749kZGR1Zru+7XXXuPAgQOcOHHiium+S5+vptLS0q7YNmTIEKKioq54fP992fft8nTfBsMvf4IVTel9eW6pq+2r7JjSZsyYwUcffQRQMhX4oEGDOHbsGFOnTmXPnj2yRkRzpTVkfMq5g6P5cE8bPjg8gQSXp6H9fvAY5uzonMoZvZjuqkHZaXUYSr1ZuHAhn3/+OQUFBSQkJPDSSy+V7Dt27BhvvPEG7u7uPPbYY2Um0KvJdN8VNchXNN23u7v7Vc9XlWnTprFw4UKefvppXnzxRR555BEWLlxYpsyPP/5Y5XlKT/ddumHNfju0rKqm+67smNJKTwWemJhYMhW4TAPevCVd2kNw7sOQt4nuvvDsoNX07T6F0A4vODu0BqHB9WKqS6df/tUV296fdmVHqUnXhzHp+rIDX0J83Co8vioy3bddXU33Xd1pwKHiqcBlGvBmypLMs0uX8On+Vnw59jw9QwIh6BXu6TINVIO78+40jf6dUEqNU0otyMjIcHYoFZLpvut2uu/qTgMOMhW4AGy5kPwSnOqIa9FObNrAzuz7ocNx8JsuyaGcRl+D0NVdctRJZLrv6nN0uu/qTAMOMhV4c1ZkKWTxpqV0M79NdKD9NuisASn8NqY1HdtU/IVCILO5OptM911/HJ0KXMZB1L9ai8tm1Trjc/3RNzN0uye+1eNffU1bT/bVOvt758ZVy2S67yaaIKpDpvu+dtcyFbgkiPp3rXEVFBbpuLglWp/sofURdPZ+N33ba2/oFT8v0TarpeoT1FFcdUWm+27GZLrvaydTgTcT2srJ018wdVEhbsZs1kw4gtHcFs+Qp/jioWlgcHV2hI2KtMgIIRo9myUH0t6EU50Jy5sEuhCDMpLotQA6xoL/7yU5OKDR1yCUUuOAcREREc4ORQhRzwryE3h1+bd8f0KxavwjuJmKcHHtyOIpNlq3nYrBKEnhWjT6GoTWernWeqavr6+zQxFC1JeCg5BwP+bT7dl2OofTmaFsunQntFoKHY7RNvw+SQ61oNHXIIQQzYQu5Nufv+G/2y/x9tAn8HfLQgHPDjuEW4soenb8CCoYRS8cJwlCCNGwFZ2F9Hcg/T2+3Hk/P1/oz+LYcfxhsB/4z6Gfa1dnR9hkNfpbTI3BO++8Q2hoKL169aJjx458/PHHDp8rLy+PYcOGYbVaAfvUEpen/C4sLGTo0KFYLJZaiVsIZ1EUsXbX10z6zzuc3D0YUl4CaxK/772Hv92UyrRfvQmhb4IkhzolCaIe7N+/n+eee459+/bx2Wef8eijjzp8rtLrQwC1skaEEA2BzWpD5+2EhDkMDL6d9fvWseV8Gz4/Php8JkHYT1w/cBmTR0zB3d2v6hOKayYJoh4cOHCgZNbU9u3bl0zG54jL60MADq8RIUSDUhTPR2vfZ+jLi9i6bTqkz8fFkMndkbE8PyqJWROfh1aLwGOQtDHUs0bfBlHdbq4JL71EwZGjtXpt125dCf3LX6osdzlBaK2ZP38+f//73x26Xun1IeCXNR3mzZtHZGQkYF95rbI1IoRoCJLS4vEu+gb3vM8g7yeSL00mPutO1sbHMLD7CHaeiCR68H1EOjvQZq7R1yAaejfXc+fOkZWVxdixYwkODuann35i2rRpDp2r/PoQUPkaEUI0KNYUSH+X5xb9jetf2cOanUsg7ydQbtzVx8Cnd1t4etI/IORfZFtkXFNDUO81CKXUB8AtQJLW+oovCEqpu4Enip9mA3/QWu+71utW55t+Xdi/fz9Dhw5l/fr1pKWlERkZydatW2nbtm3J+gtjxowhJCSE1atX8/jjj/PAAw8wb948UlNTmTJlCuPHj2fbtm28/fbbJetDgGNrRAhRn84mHGflnq1MDP8fIWoVYCXMfTwuhh6cLxoKLX8HXuNpZfSm4lU8hDM5owaxEBhTyf44YJjW+jrgRWBBfQRVVw4cOEDv3r0B+/oOkyZN4rvvviuz/sKoUaPYs2dPyfTZubm5eHh4sG/fPiZOnMgjjzyCyWQqsz4EXLmmA1y5RoQQ9cne0LwbLj0LcVH8/cuPmPtDICuPKkCB5038ZtBEdj01mAcmvAC+d4PxytUQRcNQ7wlCa70JSK1k/xat9eXFjrcBbeolsDpSOkEAjBs3jhUrVnDjjTcyZ84cZs+ezfnz5zlw4ADXXXcdmZmZJUtm7tu3r2Rhm8vbLq8PAWXXdNiyZQtAmTUihKgX1kzI+oq5/3uDfi8u4dDu30DKC1Cwj4kR2xjfOY5uEZOh0yVouxqvkHvx9gxydtSiGhp6I/UMYKWzg7gW5XsUDR06lD179vDEE09gtVoJCwsjODiYrl27Mm/ePEwmE1272vt2x8bG0rlzZ5KTkwkNDQVg9uzZvPbaa4waNQovLy+2b99e5vyffvopc+fOrZ8XJ5onrTl0egc/HNnPfV0XYS78CbCQmTWblPwItiUOJLLjKPCayM2dY7hZJslrtJSuYNH3Or+oUuHAtxW1QZQqEwO8BQzWWqdcpcxMYCZASEhI38WLF5fZ7+3tXWZN5qbigw8+4J577rliVbTLS27W1TTgzdWuXbuuudE/OzsbLy+vWoqo9lQnLq01OUVJRPjuxs+8B3/zLsYtf5Hj6e34dPSfuSH0EBlFPdiXPILUgp74uYWj1LWt2NeY3y9nuJa4YmJidmmtoyva1yAThFLqOuAr4Gat9fHqnDM6Olrv3LmzzLZdu3Y1yQQh6ldt/B1t3LiR4cOH105AteiqcRXFQ+4GLNkbiVk4lAvZ/uy96068zXkAvH3oXk7n9uWeG1rSrf0IMNbuwLVG93452bXEpZS6aoJocLeYlFJhwJfAlOomByHENdAais6QnrqVl79PIzkzi/dGPAnYPyCC3LqRU2TmtOUuera5DjxGcH+XSBm01gw4o5vrZ8BwIFApFQ88C7gAaK3fBp4BAoC3ihtmLVfLbkIIB+hCEi/tYvmeo+Rn7Gf4ybvBcgEPq4mvjnxOgbUdKQWtCGgRBR4xvDttGP5+vTEY73J25KKe1XuC0FpX+lemtb4PuK+2rmez2TAYGv14QOEkNpvN2SFcs5SMeHae2EOg6QB9/VdB/g4uJbfmbxveINy7O7OvuwAGf8yeA3llTBIdWkXiF34ajPau0gHODV84UYO7xVSbPDw8SExMJCQkRJKEqDGbzUZiYiIeHh7ODqXaMrKT2HdqNx08D9DGdSvk72D1oR78ZescxrdPo++wHwHoGmLmrsg4WpoL0OGHUa5dQBmY2NbJL0A0KI0+QVQ2F1PHjh05efIkFy5cqP/ARJPg4eFBx44dnR1GhVIykzh85gB9A/fiYdsO+Tt5+YexfHZ8DH/tt5n7eiwDoG+oF4Nan+G6sDBosxzcBmAyBTI3wt64qdy6OfmViIaq0ScIrfVyYHl0dPTvyu8zm81061b1H39T7JlQlxpqXNCwY3OUzWrhdOJR4i4eY2TYHsjfBwX7uferh9mf0pklN39I/5BDAESHnOB4Zg9atBgALX8NbtF0MXdh0fVSgxY11+gThBBNSU5uMntPHaAwP46Y1tugYB+W3COM/uRjLDYzB+/+J54u9qlWokOOY3Zxw+oxEUIfBbdoft2lO79W8t9a1A75SxKivmkN1iRizx/kx+Pn6ehziqEtN0HBYU4lenH3t2/Q2S+fmInvAmA2wMBWxzCZvMj0fBLPwG7geh3PdOkE1zggTYjKSIIQoq5oKxSdgcLjnErbwepPDnJ31w1E+mwAWxrbjt7M89tm8euOexjquwGAjn6FRIeeoVuwBYL/Ba69wK0nH3WVvkSi/kmCEOIa2KxW0rNO08IlDgpPQOEJ/vBNOEcu+fDNLQ/jY84A4FjqHD4/0YeunquJ7JYGBl/6tDVzV3YcA9t3gzbfgbk7Hi5hLI2U9gLRMEiCEKIqupDktFMcPHcaf5cL9Ao8AoUnOJWUxpgvHqaNZyLrb/tDSfG41P9wOjOI05lBXBfqCS6dGNDSQpe2yQzq/CC0WQDGULorxVyZCUY0YA4lCKXUC1fbp7V+xvFwHIqlWkuOCnE1NqsNZUtAFZ2GolN8vS+dH+MUk7pupm/AZrDEs/7ECP60+WFu7XCI14e+BkBLsytFVhMFNndsboMxuHYCcyfmjm+Jp2cI4aG7wWxf68A/biO3NrHeVaLpc7QGUX44TSgwDPsEe/Wqsm6uQgCgreTnXeDw+dPk5yUxsM0ZKIyDolNMWHILR1JD2XrHNALcMgH4OXYWXxy/mZ4+q+jrexZQdAooZGDrM3RuGQSBL4K5E+7mThx6tj0e7v7APSWX6y3NBaKJcChBaK3vLb9NKTUGkMlaRL3SNhupWYl4qAu4q3hae6xny86tfLbfh+tDjzK58zdgucCFjGBu+2oBbb3S+fH2x0qOL7KOptBq5nxuJwL8FLh0YGJUK3p1SKV/hzkQ+jq4hNFbmfl04JXXbzxjrIWoudpsg1gDfF6L5xPNnM1agMGWCJYL6KKLfLIzi6TMPB7ptxWD9RwUneP3q+9izdn+vBPzN0a320YnH9gXO4Llxx9FFZ1mcoezALT0MRAZeIF2fvngPwdc2oNLe96d2hY/nw54evwyVfz1reF6Z71oIRoQR9sgOpTb5AFMAs5dc0Si6bMVgDUBLBc5FH+RnWcy6RF4keiQY2C5yP6LimkrphLhe5YlN9unnVbAaz8sIq2gDVPClhHskQ5AkHsyvuZs8lQEeLbgwiUj13cdwj/9sugSOgbC7gdTa9wNbnxbwaD61rIcshBX5WgNIhbQ2P/fAuQCeyh9I1Y0K5aiLEw6GaxJ5Ocn8vX+LHLyc5gRdRgsiWBNYs7qEfwYH8HCUc8QFWRf6mPj/jt4dfc9/D7yB6LdFgLgTStS831IcmkBplbFj5bc0zsBrXIwhc4Dnzbg0pYXOrbh7y4+wG8BOH5yI8P7DKet9FkQ4po52gbRYDpqSy+mulFYmEd6dhIBbhkYiz/4NxzPYsdZG729f4b418GSyNpTwTy4/h5GtNnBm8NfAUBZTTy58mtMyot7w+7HoOyrFuYUDCa9wIvE3EAw5YCpJb3DAri76DR9OvSC0HfB1JK2qiXbOwUR4DMajL/0PXi4zZVxyjhiIepOox8HIb2YqqA12LLBmkJCehLHE1IJ9Uinc4tLYE3hfFoOr/wUjq85gxcHfQHWFLClMHTxfBJyA9l8+7209roEwMaDv+ejo+N4pv8KyP4GAE91HXkWN9ILfMEUBqZgXI3B3N3zFN5uLhQF/BNX1yAwhfC3O1rg4hpMgPcPYLT/6Q0Mh4EDyoZsAoI96/E9EkJUyNE2CB/gOexdWwP55VYTWuuwWolMlKUt5OencjE9GZPOpK1vJljTKSpM473tJvIKC3h0wD6wJIM1hcfWjeDH+I7MH/Yq/UP2AfDtoYn8bcd93NttPc9eb5/nx5YVwjfH36e1ZxIvRu8uuVyQezoWmwu5qge4azAFE9O1JQEBFwgw9YbW08AYQnRYEAevD8DT7WYw/Lnk+L/ffeVLaCUf+kI0Ko7WIN4C2gAvAJ8Ak4E/Al9UdaBS6gPgFiBJax1ZwX4FvAGMxd62MU1rvbt8uUZHayyWbJIzkykqSqetTw5Y08CWxme7rSRlFTGj93G8jMlgS+etHZ35+lgEWQVvMK7992DLZtPZ65m5/mlGtf2Z90a+CIBRK/6xaRkaAw91noPRYF8BLSt/MEm5fqTkeYJyB2Mg7QNcGdTmLOFBoeD3ABgDCGoRyOvjcgjyDoV228EYAMYAvnnMG2UwYO97YDe8tX2t2I0bN4L3cADMxQ8hRNPjaIK4CeimtU5RSlm11suUUjuB5cDrVRy7EJgPfHyV/TcDnYof1wP/hzN7HWoNOhdsmVgKMzh4IZX8wmwGtM0GaybYMlmyz0RssmLydccJ80oAWyZLD4fx3r6+3BaxlZmRX4EtkxOpbbn5m/l09jvDmomzSi7x1ub3OJfdmgmhz+HlcxGAlKz7OJ7emoQsZb9FhCLAw0aYdzIBnibwGAlGPwwGf2Zffw43FzOWoPkYXe0f8M/92o/njf4E+KwBs/2r+8gIGDm87MtzA24NvvJly3L0QghHE4QByCj+PVsp5QdcBKpsKdZab1JKhVdSZALwsdZaA9uUUn5KqZZa64sOxlq5jE84cmk/S98/xOQeexjQ8ijYMll9Mpw/bprMiDY7+NfQeQAUFrkycdEXuBnzOTrl9pJTLNv/NzZfjGJIwALCWu0FIDv3Fo6m3kJ8phlsaQD4uRUS5J5OC/cicB8IBj8w+jOpVwa5VgueLZ8Ab28w+jNtpA9dA85x46B/gPeHYPCmb1cDmwZf+RIeu/XKbXI7RwhxrRxNEPuwtz+sA34E3gSygeO1EFNryo6niC/edkWCUErNBGYChISE2G991FBnn8WcyezAtyfCGRi4nAF+mwEwaTcyCz3JKPTEql2x2jzQeBIZcBYPUxGJeYPReGCxeRLTNo/uAYexcTOH03+FVXvSzteHV4dewsc8jJ8Sx2DVnmhMvDrM/hI3xvUqiaFbcV/8QyfLxuauA9m7N774LWg4srOzHXqv60NDjU3iqhmJq2bqLC6tdY0fQAegY/HvQcB72EdRd6/m8eHAwavs+w4YXOr5OqBvVefs27evdkj2ev3Vd8/qr7cs0afPfKt1ziat8/bqvOxYnZqZoAuL8h07by3YsGGD065dmYYal9YNNzaJq2Ykrpq5lriAnfoqn6uOjoM4Ver3S8B9jqWnCsVTdjLANsCFWjx/WZ4x+Hkoht8wvMxmt+KHEEI0Vw1mwFsp3wBTld0AIEPXVfuDEEKIq6r3gXJKqc+w95YMVErFA88CLgBa67eBFdi7uMZi7+Z6xcyxQggh6l69JwitdaVTghffE5tVWRkhhBB1r0a3mJRSf1NKDSwezCaEEKIJq2kbRA7wCnBRKbVIKXW3Usqp62cppcYppRZkZGRUXVgIIUS11ShBaK3naq2HAF2xNyaPBg4qpbYopf6qlOpTF0FWEdNyrfVMX1/f+r60EEI0aQ71YtJap2utP9daTwVaAQ9jb2h+Ryl1Xin129oMUgghRP275kbq4kbl7cWPZ5VSwYB8nRdCiEau1nsxaa2TgKTaPq8QQoj61RAHygkhhGgAGn2CkF5MQghRN64pQSilDEqplrUVjCOkF5MQQtQNhxJE8RoNnwL52KfEQCk1Xin1t9oMTgghhPM4WoN4G/uCQe2AwuJtWwHp3iqEEE2Eo72YRgKttNZFSikN9mm/i7u4CiGEaAIcrUFkAIGlNyilwqhg1TchhBCNk6MJ4j3gC6VUDGBQSt0AfIT91pMQQogmwNFbTK9gb6B+E/sUGx8A7wBv1FJc1aaUGgeMi4iIqO9LCyFEk+boXExaa/0vrXV3rbWn1rpb8XNd2wFWIxbp5iqEEHXA0W6uTyql+pXb1l8p9adqHj9GKXVMKRWrlHqygv1hSqkNSqk9Sqn9SqmxjsQphBDCcY62QTwEHC637TD2WV0rpZQyYr81dTPQHbhLKdW9XLG/Aku01r2BO4G3HIxTCCGEgxxNEGagqNy2QsCtGsf2B2K11qe01oXAYmBCuTIa8Cn+3Re44GCcQgghHORoI/Uu4AHgX6W23Q/srsaxrYFzpZ7HA9eXK/McsEYpNQfwBEY5GKcQQtQba0YGefv3U3DsGHn79lNw/DiFZ84A0Pq1f+IztnHdLVeOtCsrpXoAa7GPezgJRAAhwI1a6/K3nsofewcwWmt9X/HzKUB/rfWcUmUeLY7tn8VdaN8HIrXWtnLnmgnMBAgJCem7ePHiGr8WgOzsbLy8vBw6ti5JXDXXUGOTuGqm0cRlseC2cxfmI4cxJSRiOncOZfvlY6qgSxdwNeO6/wD5vXuT8fuZ9RNXDcTExOzSWkdXtM+hGoTW+pBSqjMwDmgDfAl8q7XOrsbh8UDbUs/bcOUtpBnAmOJrbVVKuWEfmFdmnQmt9QJgAUB0dLQePnx4zV8MsHHjRhw9ti5JXDXXUGOTuGqmIcc1qFMnMlevIWfrVvL27MGWbf/YMwYF4jNpEp4Db8CtRw9cQkJKjjs1bjzegYH0rqPXVFfvl8MLBhUng88cOHQH0Ekp1R44j70RelK5MmexT+exUCnVDXvbxiVHYxVCCEdprSk4fpysNWtp8fXXxJ4/D4AxMBDvUaPwHjMaj+hojFV+g6/3UQDXzKEEUfzh/ncgCijzrmitwyo7VmttUUrNBlYDRuCD4hrJC8BOrfU3wGPAu0qpR7C/q9OcMcZCCNE82ZPCCbLWfU/WylUUnDgBSqE7dCDgvhl4jRyJe1QUSqnqnVApGuNHmKM1iE+xtz08BuTW9GCt9QpgRbltz5T6/TAwyMHYhBCixrTNRt6+fWStXkPmmtVYLlwEpXCPiiL02WfwGjGSzUcO08uRWznVTSQNjKMJogcwqHyjsRBCNCbaaiX/wAEyV63+JSm4uODZL5rAmTPxHjkSU1DQLwccqbQPThUXu/Z465ujCWIT0Bt7d1enkrmYhBA1oW028vbuI2v1ajJXrsSSlAQmE56DBhL04IN4jxiB0cen6hPVRDOrQZwGViulvgQSSu8ofauoPmitlwPLo6Ojf1ef1xVCNB7aaiVv334yV60ka81aLAkJYDTiOWAAQQ8/jPfIERjrej63ZtQG4Qksxz6Ta9sqygohRL3TNhu5P/9M5po1ZH+/DsulS+DigtegQfg8+gieQ4Zg8vevn2CUaj4JQmt9b20HIoQQ10oXFZGz7Wey1qwha906rKmpKHd3vAYPxnvUSLyGD6/7mkJFGucdJsfHQQAopbyxD2Arefla61PXGpQQQlSXLS+PnK1byVq9hqwNG7BlZqI8POxJYfRNeI8cicGtOtPE1bHmUoMonn11EdALe9u84pc2emPthCaEEBWzpKWR89NPZK1bT/YPP6Dz8jD4+uIdMxzv0aPxHDQIg6urs8MsoRppFcLRGsRbwAYgBogDwoG5wJbaCUsIIcoqunjRnhDWrydn+3awWDC2aIHvhPF4j7oRj/79MJjNzg7z6ppLDQJ7zeFGrXWRUkpprTOUUn8EDgKf1F54VZNurkI0TVprCo4cIWv9BrI3bCD/0CEAzOHhBNw7De9Ro3Dr2RNlcHTVgnqkFLoRDoRwNEHkY+/BVAQkK6XCgDQgoLYCqy7p5ipE02HLzyf355/xXvQpsc89b++OqhTu111H0KOP4j1yBK4dOzo7zJprZuMgfgR+AywElgIrgQJgfe2EJYRoLoouXiT7h01kb9pEzpYt6Px83M1m3IcNw2vOHLyGD8MUUO/fPWtfc7nFpLX+Tamnf8F+a8kb+Kg2ghJCNF26qIi8vXvJ3vQj2Zs2UXDsGAAurVrhd9uteMXEsDM/n+433ujkSGtRc6pBKKUe11rPAyiej+mT4u2PAq/VXnhCiKagMP48OVs2k7NlKzlbtmDLzASTCY+oKIIffwyv4cMxd+z4y+yoGzc6Nd460fgqEA7fYnoGmFfB9r8iCUKIZs+Wk0Purl3kbN1G9saNFMbFAWAKDcV7xAi8RsTgOXBgNdZQaCKaQw1CKTWi+FejUiqGsuMDOwBZtRVYDWKSXkxCOJm22cg/fIScrVvI2byFvF270EVFJTOj+v3mN3gNGVy2ltDcNIM2iPeLf7oBH5TaroFEYM4VR9Qx6cUkRP3TWlN05gzZW7aQu3UruTt2Yk1PB8C1Uyf8J0/Gc+BAPPr2weDh4eRoGwBF008QWuv2AEqpj7XWU+smJCFEQ1R0/jw523eQu2MHOdu22tdOAEytWuIVE4PngOvxuOEGXIKDnRxpw9OsRlKXTw7Ft5ssWusfq3O8UmoM8Ab2aTne01q/XEGZ3wDPYa+d7NNal1+3WghRh4ouXCDn5+3k7d9Hzk+bKTp3DuD/27v36LiqeoHj399kZpI0j6avJG2TJi19QIttCpUiBWlRniq4rqBFrKBgfdwKLNSriHIvKN6rXl1epWt5y2NdRKS+sXCreFGCKIotAoUWKKUtzbS0KaVNSZpkXr/7xzkp0+lMMpNk5swkv89as3LOnD1nfme3M7/Z+5yzN77qaioWL2bMNddQuWQJgWnTRm+3UTZGeguij4g8BnxZVf8iIl8EbgCiIrJaVb8xwGtLgNXAuUAI2CAi69xpRvvKzAJuxJm17qCI2E8SY3JIVYm0tXFkwwaO/H0DRzZuJLJ7NwBSXk7F6acz7ooPU/GOMyideQJSYkOuZaVIE+hgr2I6Gfibu/wJYCnQCfwF6DdBAKcBDcKiqQAAF3dJREFU2/pGfRWRtcAlQOJcfp8AVqvqQQBVbR9knMaYFDQapfeVVyhvbSX0wG/ofvppovv2AVAybhxjFi1i/EdXMGbxYkpnzbKEMCyKrwUhOohmj4j0DasxHfi9qp7gPv+mqlYN8NpLgQtU9Rp3fQWwWFVXJZR5ANgKLMHphvo3Vf1din2tBFYC1NXVnbp27dqsjwWgs7OTygK83M7iyl6hxuZ1XHL4MMEdOwi4D//OV/H19gIQq6khMnMm4ZknEJ49h9jkes9/8XpdX+kMNq5x3/o2Ggxw6PrrcxDV0Opr2bJlT6nqolTbBtuC+DNwOzAZ+DWAiJwAvJ7Ba1P9z0vOUn5gFk7LpAF4XEROVtVDx7xIdQ2wBmDRokW6dOnSzI8gQWtrK4N9bS5ZXNkr1NjyGZeGw/S89BLdzzxL97POo+/8AX4/ZSeeSPkHPkD5gvlsCoc589JLC+4cwkj7d9y55g6kNEhLjo4pV/U12ARxFfA5YD/wbfe5E3FOPA8kxLHTlDYAe1KU+ZuqRoAdIvISTsLYMMh4jRmRVJXo3r1OInATQs/mzWg4DIC/ro7ylhbGXX455S0LKJs795jJc2KtrQWXHEas4uthGvRVTAdwxmBKfO5/M3z5BmCWiEwHdgPLgeQrlB4ALgf+R0QmArMBm6nOjHrx3l56nnvumIQQbXdO0UlpKWXz5jHuiisoX7CA8pYFBOrrPY7YAJ532Q1WxglCRG5S1dvc5VvTlVPVm/vbj6pGRWQV8DDO+YW7VXWzu8+NqrrO3XaeiGwBYsAX3KRkzKih4TA9L79M79aX6XluEz1bt9L97CaIRAAITJvGmMWLnWSwYAFlc2YjhTxhzmg3wi9zbUhYbkxbKgOquh5Yn/TczQnLinPp7A1DeR9jikW8t5fel16iZ8sWejZvpmfzFnpefvloMvCNGUPprFmMX7GCMaeeQnlLy8gYAnu0GOl3UqvqpxOWP5abcIwZ+eLd3fS8+KKTBNyE0LttG8RiAJSMHUvZvLlMuOpKyubOpXTOiQSnNSL+wZ4yNF4bVXdSi8hc4CxgPPAG8HjijW75ZIP1mUIW6+yi98UXKP/DH9nz29/Ss2ULva9sh3gcgJLx4ymbN4/KpUspmzeXsrnzCEydYieOR6KR3IIAEOd/7V3AlThXGu0BpgJTRORe4OM6mBsrhsAG6zNeU1W6n3kGX0UFnY+2cuTJJympGUt456v0vPgixONUA12TJlE2bx5V557nJIN58/DX1VkyGA2K9N842xbESpx7E05X1aOXnIrI24H7gU8CPxy26IzxmKqikQgSCNC79WWCzU0Qj3Nk41PE3jhA7/YddD3xBD3PPXf0NSXjxoHPR9mc2Uz81Ccpmz+fpzs6eOcll3h4JMZzI70FAawArk1MDgCqukFErscZP8kShCk6GokQ2b2b3p07Ce/ciYhQOmcOe2+5lej+/QQm19P78rajQ1fHjxw5+tpgczPjr7wSf20tFWcuoWzOnOP2Hx+JM6SZzImMigQxF3gszbbHgHuHFo4xuRXr7KJny2bCO5xEEN6xw/kbCkE0elx5/6RJlJ14Ir3btzPpumvp3fYKEghQefY7CTQ0UDpzJr7ycg+OxBQbLcI75bJNECWqmnLWOFV9U0R8wxCTMUMW6+hwprn0ByibexJEoxy46y4O3HU38c5OwLmxLNjcTOmcOVSdfz7B5maCzU0Em5uJ7t/PkSf/TvV7LsI/frzHR2OKnsiouJM6kGKq0aHsb8jsKqbRS2MxIq+9Rnj7dnq3bye8fQfjnn6arTd9hdiBt+6r9E+ZDJEo0f37qTr33dRcdhnBGScQmDIZ8aX+TeMfN46y2bPzdShmpBslJ6nbOXaq0VTb88quYhrZVJXYoUNON9DOV92/bz3UHZEUnPsHmDiRymVLKZ0+g+D06cQOd3B4/Xoiu/fQcMstVJ2zzMOjMaPaSD8HoarNOYrDjHIaixHv6qKkuprYoUMcuOtuOlsfJdwWQnt63iro9xNsaCDY1ETFGWcQnDGd0hkzCM6YgX/cOFpbW1mQNKplzfvfn9+DMSZZcTYg8t8lZEaveFcX4VCISFsb4V1tREJthNtCRHbtIrxnD0QiBJuaiL7+OvHubirOOIOKM88iMLmeYFMTwaYmAlOnIoGA14diTPZGegvCmP5oPE60vZ3wrl1E2kKEQ21E2tyEEAodc14AwFdVRbCxkdKTTqLqvHPxVVTS9de/EmhuovaGG1JeLmpMMRIR8nwP8bCwBGGyEuvswh8K8eYjjzi//tva3koEoRDqDi4HgM9HYPJkAo2NVJ2zjEDjNIKNDQQaGglOa3TOGSSZ+KlP5vFojMmX4uxjKvoEYVcxDR9VJXbgAJE9e4jsec396z5ec9bjHR1MwBlnBRJaAbNnU/Wucwg0NBJobCDY2Ehg8mTrDjKmj7Ug8s+uYsqeqh5t8nY98QQdv/wlPS++RGT37mOuCgLwVVQQmDKFwJQpjFnYgn/yZF7uOEzL+ecTbGygpKbGo6MwpoiMkstcTYGLdXYS3buXyL59RPfuI9q+j8jefUT37XOe27ePWEcHZSedhIbD9G7dSsmECZQvbKHy7LOdZDB1ytGk4KuqOm4wuedbWyl/28keHaExxWk03Ek9LETkApz5q0uAO1X1P9KUuxT4OfB2Vd2YxxALjsbjxA4eJLJ3L9F97UT37U344nef27v3mDGC+pSMG4e/vp5AbS3l8+fjq6ig+9lnAZh829epft/78NlMZMbkzii5k3rIRKQEWA2ci9OVvUFE1iXPJyEiVcC1wJP5jjHfNOrc5Rvd1/drfy+Rfe2M3bSJnXfcSdT95X/MCWCAkhL8kyYRqKujdOZMKs5cQqCuHn9dHYH6Ovx1dfhra/GVlnpzYMYYh3UxZew0YJuqbgcQkbXAJUDyhENfA74FfD6/4Q0fVSV++DDR9nai+/cTcf8ebQG4v/qjr79+dAKZPlJair+6GmlupnzhQgJ1tfjr6vHX1xGoq3OWJ05ASkq8OThjTHbsJHVGpgJtCeshYHFiARFZCDSq6kMiUnAJom/4h2j7fucLf//+o0kg+a+Gw8e93ldZ6X7R11M6c6b7a7/+rV/9dXWU1NTw2GOP8baku4KNMUWoOBsQSL5v3hCRy4DzVfUad30FcJqqftZd9wF/BK5S1Z0i0gp8PtU5CBFZiTOJEXV1daeuXbt2UDF1dnZSWVkJ8TjS1UVJRwe+hEfJoQ58hzvwHepwth0+jKQYGjpeVka8pob42GpiY52/8bE1xMZWO89XjyU+thotK8surgJTqHFB4cZmcWVnpMVV84Pb8XV28saNX8pBVEOrr2XLlj2lqotSbfOiBRECGhPWG3CmLu1TBZwMtLpXz9QD60Tk4uQkoaprgDUAixYt0qWD+LXd+fjjvLJ6NWN6ep2unuR+fsBXXY1/0iT89XX433YygdpaZz3p73DPC9Da2spgjinXCjUuKNzYLK7sjLS4dt1/PzER5ufomHJVX14kiA3ALBGZDuwGlgMf7tuoqh3AxL71/loQw0ECQbS0jIr5C/DXTkr9xZ/hL35jjEnLzkEMTFWjIrIKeBjnMte7VXWziNwKbFTVdfmMp+L0xRy6/jpaCvDXijFmZBBGx5Sjw0JV1wPrk567OU3ZpfmIyRhjcqZIL3O1KUKNMSYPivFO6qJPECLyPhFZ09HR4XUoxhiTmrUgvKGqD6rqyrEpho42xpiCUXwNiOJPEMYYU/CkOE9SW4Iwxphcsy4mY4wxaVkLwhhjzHGKswFhCcIYY/LCWhD5Z5e5GmMKndhJam/YZa7GmMJXnH1MRZ8gjDGmOFgLwhhjTDK7zNUYY0w6+Z6cbThYgjDGmFwTKcYepuJPEHYVkzHG5EbRJwi7iskYU/DsMldjjDEpFec5am8ShIhcICIvicg2EflSiu03iMgWEdkkIn8QkSYv4jTGmGFjLYiBiUgJsBq4EJgLXC4ic5OKPQ0sUtX5wC+Ab+U3SmOMGT5il7lm7DRgm6puV9UwsBa4JLGAqj6qqkfc1b8BDXmO0RhjhlcRtiAk39fmisilwAWqeo27vgJYrKqr0pS/Hdirql9PsW0lsBKgrq7u1LVr1w4qps7OTiorKwf12lyyuLJXqLFZXNkZaXGNvfMu/Lt2ceDWW3IQ1dDqa9myZU+p6qKUG1U1rw/gMuDOhPUVwA/SlP0ITguidKD9nnrqqTpYjz766KBfm0sWV/YKNTaLKzsjLa7QDZ/TbeedP7zBJBhKfQEbNc33qn9QKWdoQkBjwnoDsCe5kIi8G7gJOFtVe/MUmzHG5IQW4Z1yXpyD2ADMEpHpIhIElgPrEguIyELgv4GLVbXdgxiNMWb42EnqzKhqFFgFPAy8APxMVTeLyK0icrFb7NtAJfBzEXlGRNal2Z0xxhQ8KfFBNOZ1GFnzoosJVV0PrE967uaE5XfnPShjjMkRKS0j3lt8PeV2J7UxxuSYlJWiPT1eh5G1ok8QNlifMabQ+awF4Q21wfqMMQVOykohEkFjxXUeougThDHGFDpfWRlA0XUzWYIwxpgck1InQRRbN5MlCGOMybGSKmcYjFiRnSu1BGGMMTnmr60FILp/v8eRZKfoE4RdxWSMKXRHE0R7cSUIT26UG06q+iDw4KJFiz7hdSzGGJOKv64egEioLW0ZVSXy6qtE9uwhHg5TUlVFSU0NgSlT8JWX5yvUYxR9gjDGmEJXUllBsKmJ7uefT7m9+7nneO2mr9C7dWvq10+aSLChkeC0Rqrf+z4qzzozl+EeZQnCGGPyoLylhc7WVjQcRoLBo89HDxxg11Ufwze2mrqbv0rZ7NlIMEjszTeJvXGQyO7dhENtRNpCdD72JzoefIjGNWuoPHNJzmO2BGGMMXlQdeEFdPzmN7z5aCvV55939PmD991H/MgRmn/+M0pnzOh3H/GuLnZ++ApC117L1P/8NlXnnJPTmIv+JLUxxhSDyjPPJDBtGq+vXo3G4wDEjxzh4H0/ofKccwZMDgC+igqm3XkHwaYmQp/5Z3Z9/GrefOQRiERyEnPRJwi7iskYUwzE76f2+uvo3bqVgz/+MQCHfvVrYh0dTLj64xnvxz9pEtN/upbaL3ye3m3bCK36LJUP/CYnMRd9grCxmIwxxaLqwgupPPts2r/zXbo3beKNu++mvKWFMaecktV+JBhkwtVXM/MPj9B4xxq6c3Q+ougThDHGFAsRof5rt1JSU8POD36IyJ49TPzsqsHvLxCg8qyziE2ePIxRvsWTBCEiF4jISyKyTUS+lGJ7qYj81N3+pIg05z9KY4wZfoHaWpp/ch/jVqxgyre+SeWS3F+NNFh5v4pJREqA1cC5QAjYICLrVHVLQrGrgYOqOlNElgPfBD6U71iNMSYXAlOnUn/Tl70OY0BetCBOA7ap6nZVDQNrgUuSylwC3OMu/wJ4l0iRzvptjDFFyosEMRVIvN885D6XsoyqRoEOYEJeojPGGAN4c6NcqpaADqIMIrISWOmu9ojI5n7edyxOokllIvD6IF870PahvNbiyn57f7FZXBZXocY10PZcxtWUdouq5vUBvAN4OGH9RuDGpDIPA+9wl/3ugcsA+10z2O3Axhzu2+LKU1wDxWZxWVyFGtcwxD3ouPp7eNHFtAGYJSLTRSQILAfWJZVZB1zpLl8K/FHdI+3Hg0Pcnqt9W1zZbbe4sttucWW3vVDjGmh7LuNKSwb+3h1+InIR8D2gBLhbVW8TkVtxsuA6ESkD7gUWAm8Ay1V1ew7j2aiqi3K1/8GyuLJXqLFZXNmxuLKTq7g8GaxPVdcD65OeuzlhuQe4LI8hrcnje2XD4speocZmcWXH4spOTuLypAVhjDGm8NlQG8YYY1IaVQmiUIf4yCCuq0Rkv4g84z6uyVNcd4tIu4iknAZLHN93494kItmNOJa7uJaKSEdCfd2cqtwwx9QoIo+KyAsisllErktRJu/1lWFcea8v933LROTvIvKsG9stKcrk/TOZYVxefSZLRORpEXkoxbbhr6vBXv5UbA+cE+KvADOAIPAsMDepzGeAH7rLy4GfFkhcVwG3e1Bn7wROAZ5Ps/0i4Lc4962cDjxZIHEtBR7Kc11NBk5xl6uArSn+HfNeXxnGlff6ct9XgEp3OQA8CZyeVMaLz2QmcXn1mbwB+Emqf69c1NVoakEU6hAfmcTlCVX9E85VZOlcAvxIHX8DakQkN8NKZhdX3qnqa6r6D3f5TeAFjh8hIO/1lWFcnnDrodNdDbiP5JOief9MZhhX3olIA/Ae4M40RYa9rkZTgijUIT4yiQvgA263xC9EpDHHMWUq09i98A63i+C3IjIvn2/sNu0X4vzyTORpffUTF3hUX26XyTNAO/B/qpq2zvL4mcwkLsj/Z/J7wL8A8TTbh72uRlOCGLYhPoZZJu/5INCsqvOBR3jrV4LXvKivTPwDaFLVBcAPgAfy9cYiUgn8ErheVQ8nb07xkrzU1wBxeVZfqhpT1RagAThNRE5OKuJJnWUQV14/kyLyXqBdVZ/qr1iK54ZUV6MpQYSAxCzfAOxJV0ZE/DhjmOS6K2PAuFT1gKr2uqt3AKfmOKZMZVKneaeqh/u6CNS55yYgIhNz/b4iEsD5Er5PVX+Voogn9TVQXF7VV1IMh4BW4IKkTV58JgeMy4PP5BLgYhHZidMNfY6I/DipzLDX1WhKELka4iPncSX1U1+M049cCNYBH3Wvzjkd6FDV17wOSkTq+/peReQ0nP/nB3L8ngLcBbygqt9NUyzv9ZVJXF7Ul/tek0Skxl0uB94NvJhULO+fyUziyvdnUlVvVNUGVW3G+Y74o6p+JKnYsNeVJ3dSe0FVoyKyCmcgwL4hPjZLwhAfOB+ke0VkG+4QHwUS17UicjEQdeO6KtdxAYjI/ThXuEwUkRDwrzgn7FDVH+LcDX8RsA04AnysQOK6FPi0iESBbpyhWnKd6JcAK4Dn3L5rgC8D0xLi8qK+MonLi/oC5wqre8SZRMwH/ExVH/L6M5lhXJ58JpPluq7sTmpjjDEpjaYuJmOMMVmwBGGMMSYlSxDGGGNSsgRhjDEmJUsQxhhjUrIEYQwgIl8WkXRj3BgzKlmCMKOCiHQmPOIi0p2wfoWqfkNV8zJksxvPh0XkJymeXyoiKiK/Snp+gft8a75iNGbU3ChnRjdVrexbdocruEZVH/EuIi4iadrdBPuBM0Rkgqr23dF8Jc5Q3VkTEb87eJsxWbEWhDGAiPxb39g2ItLs/lr/mIi0ichBEfmUiLzdHb3zkIjcnvT6j4szKc9BEXlYRJr6eS8fcC7wuzRFwjgD5i13y5cAHwTuS9jHahH5TtJ+HxSR693lnSLyRRHZBHS5Y/MYkxVLEMaktxiYBXwIZ6jlm3DG5ZkHfFBEzgYQkffjDF/xT8Ak4HHg/n72exqwXVVf76fMj4CPusvnA5s5dmC/e4DL3WSDO7jeu5Le93Kc+QNqrAVhBsMShDHpfU1Ve1T190AXcL+qtqvqbpwksNAt90ng31X1BfeL+BtASz+tiPeQvnsJAFV9AhgvInNwEsWPkrb/HWe8/3e5Ty0HWlV1X0Kx76tqm6p2Z3rAxiSyBGFMeolftt0p1vvOazQB/+V2PR3CGShNSD8ZUH/nHxLdC6wClgG/TrH9HqBvRM+PuOUTtWHMEFi/pDFD1wbcpqr3DVRQROpxRgv9Rwb7vRdn5NcfqeoROX72yB8Dz4vIAuAkjp/ox0biNENiLQhjhu6HwI3iTtUpImNF5LI0ZS8CfpfJcNqqugM4G+fcR6rtIZz5RO4FfmldSWa4WYIwZohU9dfAN4G1InIYeB64ME3xTLuX+vb9Z1Xtb9a5e4C3cXz3kjFDZvNBGJMn7qWme4ETVLVjmPb5TpyupmZVTTeZvTGDYi0IY/JnPPDVYUwOAeA64E5LDiYXrAVhTBESkZOAjcCzwAWqetjjkMwIZAnCGGNMStbFZIwxJiVLEMYYY1KyBGGMMSYlSxDGGGNSsgRhjDEmJUsQxhhjUvp/HSHfUyszkY8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as mticker\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, sharex=True, figsize=(6, 5),\n",
    "                               gridspec_kw={'height_ratios': [1, 3]})\n",
    "fig.subplots_adjust(hspace=0)\n",
    "\n",
    "ax1.set_ylabel(\"$M(t)$ / $M_{\\odot}$\", fontsize='large')\n",
    "ax1.yaxis.set_minor_locator(mticker.AutoMinorLocator())\n",
    "ax1.plot(ts2/1e6,mass, color='tab:red')\n",
    "ax1.grid()\n",
    "\n",
    "ax2.set_xlabel('Time / Myr', fontsize='large')\n",
    "ax2.xaxis.set_minor_locator(mticker.AutoMinorLocator())\n",
    "ax2.set_ylabel('Distance / au', fontsize='large')\n",
    "ax2.set_ylim([0, 1.6])\n",
    "ax2.yaxis.set_minor_locator(mticker.AutoMinorLocator())\n",
    "\n",
    "ax2.plot(ts2/1e6, a2, lw=2, color='gold',\n",
    "         label='$a_\\oplus(t), dt_{split}=400$ yr')\n",
    "ax2.plot(ts1/1e6, a1, ':', lw=2, color='tab:blue',\n",
    "         label='$a_\\oplus(t), dt_{split}=4000$ yr')\n",
    "ax2.plot(ts2/1e6, radius, color='tab:red', label='$R_{sun}(t)$')\n",
    "ax2.legend(fontsize='medium', loc='upper left', labelspacing=0.1,\n",
    "           framealpha=1.0)\n",
    "ax2.grid()\n",
    "\n",
    "plt.savefig('img/fig2.eps', bbox_inches='tight', pad_inches=0.01)\n",
    "plt.savefig('img/fig2.pdf', bbox_inches='tight', pad_inches=0.01)"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "5143b235c2103413b967559f03b877662679c3138437357f05bf8079dc5edae9"
  },
  "kernelspec": {
   "display_name": "REBOUNDx (3.1.0)",
   "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.9.5"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}