PS-N1-Solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2a8e87d6-2e37-4118-bc61-b506adc84a0f",
   "metadata": {},
   "source": [
    "# Problem Set N1 Solutions\n",
    "\n",
    "_PHY 323 Spring 2024, Prof. Ian Beatty_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e50c094-80d1-4252-8ba2-0304ea6758ec",
   "metadata": {},
   "source": [
    "## Problem 1: Burn Rate"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4df91f5b",
   "metadata": {},
   "source": [
    "### a) Plot the rocket's acceleration vs. time function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc3193ef",
   "metadata": {},
   "source": [
    "The function we're given is:\n",
    "\n",
    "$$ a(t) = A\\,t^2\\,e^{-C\\,t^2} $$\n",
    "\n",
    "I could _sketch_ this qualitatively by hand, but to make a more careful plot, I'll need to calculate the acceleration numerically for each of many different time values and plot the $(t,a)$ pairs. I can't do this without having numerical values for the parameters $A$ and $C$.\n",
    "\n",
    "So, I'll transform my equation into a dimensionless form, rescaling time and distance to absorb those pesky parameters.\n",
    "\n",
    "If I define $\\underset{\\sim}{t} = c^{1/2}\\,t$, the exponent becomes parameter-free and necessarily dimensionless. (Note: I'll use under-tildes to signify the dimensionless, scaled versions of variables.) I'll have to replace all occurrences of $t$ with $\\underset{\\sim}{t}/c^{1/2}$ in the equation, and $dt$ with $d\\underset{\\sim}{t}/c^{1/2}$ in the derivative. (Note that the second derivative involves $dt^2$.)\n",
    "\n",
    "$$ \\frac{d^2 x}{dt^2} = A\\,t^2\\,e^{-C\\,t^2}\n",
    "    \\ \\Rightarrow \\ \n",
    "    C\\, \\frac{d^2 x}{d\\underset{\\sim}{t}^2} = A\\, \\frac{\\underset{\\sim}{t}^2}{C}\\, e^{-\\underset{\\sim}{t}^2}\n",
    "    \\ \\Rightarrow \\ \\frac{C^2}{A}\\, \\frac{d^2 x}{d\\underset{\\sim}{t}^2} = \\underset{\\sim}{t}^2\\, e^{-\\underset{\\sim}{t}^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53523ad2",
   "metadata": {},
   "source": [
    "If I now define a scaled position variable $\\underset{\\sim}{x} \\equiv C^2\\,x / A$, it absorbs the remaining combination of parameters to leave:\n",
    "\n",
    "$$ \\frac{d^2 \\underset{\\sim}{x}}{d\\underset{\\sim}{t}^2} = \\underset{\\sim}{t}^2\\, e^{-\\underset{\\sim}{t}^2} $$\n",
    "\n",
    "_This_ is a parameterless, dimensionless equation that I can deal with numerically! Here's a plot of the dimensionless acceleration $d^2 \\underset{\\sim}{x} / d\\underset{\\sim}{t}^2$ vs. $\\underset{\\sim}{t}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d39c5244",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2AAAAFDCAYAAABcNW0VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB2vklEQVR4nO3dd1wT9/8H8FcSSNh7qaCAoogiKArurdRV92wVt3XUWtta7XBUW0dd/bZaR+uoVVu3tcO9xclQHLgVRUER2Ssk9/sDyc8IKpDAMV7PhzxMPrncve4+Ocg7ufucRBAEAURERERERFTspGIHICIiIiIiqihYgBEREREREZUQFmBEREREREQlhAUYERERERFRCWEBRkREREREVEJYgBEREREREZUQFmBEREREREQlhAUYERERERFRCWEBRkREREREVEJYgBGVMTNnzoREItFqc3V1xdChQ8UJpEf5rZu+HT16FBKJBEePHi3W5ZQnEokEM2fOFDtGHkOHDoWrq6vYMSq0V18b69atg0Qiwb1790TLpC8l8bpv3bo1WrduXazL0Cfuc0T6wQKMiIhKrUePHmHmzJkIDw8XOwpRhcB9jqj4GYgdgIh0d/36dUilZf/zlK+++gpTp04VOwaVIo8ePcKsWbPg6uoKX19frcdWr14NtVotTjDK1+DBgzFgwAAoFAqxo+gsPT0dBgYV720S9zmi4lfxfrMQlUPl4c0OABgYGFTINzwVSUZGBuRyuV4+MDA0NNRDotIpOzsbarUacrlc7CiFIpPJIJPJxI6hF0ZGRmJHKHXK8z5HVJLK/kfmROXYyZMn0ahRIxgZGaF69epYuXJlvtO9eg5Y7nkYJ0+exMSJE2Fvbw8rKyuMGTMGWVlZSEhIwJAhQ2BtbQ1ra2tMmTIFgiBozVOtVmPp0qWoU6cOjIyM4OjoiDFjxuD58+d5lt21a1ecPHkS/v7+MDIygru7O3777Tet6ZRKJWbNmgUPDw8YGRnB1tYWzZs3x4EDBzTT5HcOWHZ2NmbPno3q1atDoVDA1dUVX3zxBTIzM4uU43XOnj2Ld955B5aWljAxMUGrVq1w6tQprWmSk5MxadIkuLq6QqFQwMHBAR06dEBoaKhmmps3b6J3795wcnKCkZERnJ2dMWDAACQmJr5x+SdOnEDfvn1RtWpVKBQKuLi44OOPP0Z6enqeaSMjI9GvXz/Y29vD2NgYtWrVwpdffqk1TXR0NEaMGIHKlStDoVDAzc0NY8eORVZWlmaahIQETJo0CS4uLlAoFKhRowbmz59foE+4o6OjMXz4cDg6OkKhUKBOnTpYs2aN1jS559v98ccf+Oqrr1ClShWYmJggKSkJ8fHx+PTTT+Ht7Q0zMzNYWFigU6dOuHjxotbzGzVqBAAYNmwYJBIJJBIJ1q1bByD/81FSU1PxySefaNapVq1aWLhwYZ7Xt0QiwYQJE7Br1y7UrVtXsw579+5943rHxsbCwMAAs2bNyvPY9evXIZFI8NNPPxVqG9+7dw8SiQQLFy7E0qVLNa/1q1evAgB+/PFH1KlTByYmJrC2tkbDhg2xadMmzfNfd15OfvvTgQMH0Lx5c1hZWcHMzAy1atXCF1988cZ1BoDMzEx8/PHHsLe3h7m5Od599108fPgwz3T5nQOWu28ePXoUDRs2hLGxMby9vTXnYe7YsQPe3t4wMjKCn58fwsLC8sw3MjISffr0gY2NDYyMjNCwYUP89ddf+S771KlTmDx5Muzt7WFqaoqePXvi6dOnWtNeuHABgYGBsLOzg7GxMdzc3DB8+HCtafI7BywsLAydOnWChYUFzMzM0K5dO5w5c6bIOfKTmZmJGTNmoEaNGprfBVOmTMnzO68gffm2186rCrvPvfzaXbZsGdzd3WFiYoKOHTviwYMHEAQBs2fPhrOzM4yNjdG9e3fEx8fnWe5///2HFi1awNTUFObm5ujSpQuuXLny1m1FVFbxo2aiUioiIgIdO3aEvb09Zs6ciezsbMyYMQOOjo4FnseHH34IJycnzJo1C2fOnMGqVatgZWWF4OBgVK1aFd999x3+/fdffP/996hbty6GDBmiee6YMWOwbt06DBs2DBMnTsTdu3fx008/ISwsDKdOndL6JPTWrVvo06cPRowYgaCgIKxZswZDhw6Fn58f6tSpAyDnzeDcuXMxcuRI+Pv7IykpCRcuXEBoaCg6dOjw2nUYOXIk1q9fjz59+uCTTz7B2bNnMXfuXFy7dg07d+7UmrYgOfJz+PBhdOrUCX5+fpgxYwakUinWrl2Ltm3b4sSJE/D39wcAfPDBB9i2bRsmTJgALy8vPHv2DCdPnsS1a9fQoEEDZGVlITAwEJmZmZptHx0djb///hsJCQmwtLR8bYatW7ciLS0NY8eOha2tLc6dO4cff/wRDx8+xNatWzXTXbp0CS1atIChoSFGjx4NV1dX3L59G3v27MG3334LIOcQIn9/fyQkJGD06NHw9PREdHQ0tm3bhrS0NMjlcqSlpaFVq1aIjo7GmDFjULVqVQQHB2PatGl4/Pgxli5d+tqssbGxaNy4saaIsbe3x3///YcRI0YgKSkJkyZN0pp+9uzZkMvl+PTTT5GZmQm5XI6rV69i165d6Nu3L9zc3BAbG4uVK1eiVatWuHr1KipXrozatWvjm2++wfTp0zF69Gi0aNECANC0adN8cwmCgHfffRdHjhzBiBEj4Ovri3379uGzzz5DdHQ0lixZojX9yZMnsWPHDowbNw7m5ub43//+h969eyMqKgq2trb5LsPR0RGtWrXCli1bMGPGDK3H/vzzT8hkMvTt2xcACr2N165di4yMDIwePRoKhQI2NjZYvXo1Jk6ciD59+uCjjz5CRkYGLl26hLNnz2LQoEGv7aP8XLlyBV27dkW9evXwzTffQKFQ4NatW3k+aMjPyJEj8fvvv2PQoEFo2rQpDh8+jC5duhR42bdu3cKgQYMwZswYvP/++1i4cCG6deuGFStW4IsvvsC4ceMAAHPnzkW/fv20Dqu+cuUKmjVrhipVqmDq1KkwNTXFli1b0KNHD2zfvh09e/bUWtaHH34Ia2trzJgxA/fu3cPSpUsxYcIE/PnnnwCAJ0+eaH63Tp06FVZWVrh37x527Njx1u3XokULWFhYYMqUKTA0NMTKlSvRunVrHDt2DAEBAYXKkR+1Wo13330XJ0+exOjRo1G7dm1ERERgyZIluHHjBnbt2qXJ8ra+LMprp7D7XK6NGzciKysLH374IeLj47FgwQL069cPbdu2xdGjR/H555/j1q1b+PHHH/Hpp59qfVizYcMGBAUFITAwEPPnz0daWhp+/vlnNG/eHGFhYRz0g8ongYhKpR49eghGRkbC/fv3NW1Xr14VZDKZ8OquW61aNSEoKEhzf+3atQIAITAwUFCr1Zr2Jk2aCBKJRPjggw80bdnZ2YKzs7PQqlUrTduJEycEAMLGjRu1lrN379487dWqVRMACMePH9e0PXnyRFAoFMInn3yiafPx8RG6dOnyxnWeMWOG1rqFh4cLAISRI0dqTffpp58KAITDhw8XOseRI0cEAMKRI0cEQRAEtVoteHh45NlWaWlpgpubm9ChQwdNm6WlpTB+/PjX5g8LCxMACFu3bn3jeuYnLS0tT9vcuXMFiUSi9Rpo2bKlYG5urtWWux65hgwZIkilUuH8+fN55pk73ezZswVTU1Phxo0bWo9PnTpVkMlkQlRUlKYNgDBjxgzN/REjRgiVKlUS4uLitJ47YMAAwdLSUrMuudva3d09z/plZGQIKpVKq+3u3buCQqEQvvnmG03b+fPnBQDC2rVr86xLUFCQUK1aNc39Xbt2CQCEOXPmaE3Xp08fQSKRCLdu3dJaJ7lcrtV28eJFAYDw448/5lnWy1auXCkAECIiIrTavby8hLZt22ruF3Qb3717VwAgWFhYCE+ePNGatnv37kKdOnXemOfV7ZDr1f1pyZIlAgDh6dOnb5zfq3L3w3Hjxmm1Dxo0KM9rI/d3z927dzVtuftmcHCwpm3fvn0CAMHY2FjrtZy7bXP3T0EQhHbt2gne3t5CRkaGpk2tVgtNmzYVPDw88iy7ffv2WvvDxx9/LMhkMiEhIUEQBEHYuXOnACDf/eNlr65bjx49BLlcLty+fVvT9ujRI8Hc3Fxo2bJloXMIgiC0atVK63fvhg0bBKlUKpw4cUIry4oVKwQAwqlTpwRBKFhfFuS1k5/C7HO5r117e3ut9Zo2bZoAQPDx8RGUSqWmfeDAgYJcLtf0ZXJysmBlZSWMGjVKazkxMTGCpaVlnnai8oKHIBKVQiqVCvv27UOPHj1QtWpVTXvt2rURGBhY4PmMGDFC6xCkgIAACIKAESNGaNpkMhkaNmyIO3fuaNq2bt0KS0tLdOjQAXFxcZofPz8/mJmZ4ciRI1rL8fLy0nxSCgD29vaoVauW1jytrKxw5coV3Lx5s8D5//33XwDA5MmTtdo/+eQTAMA///xT6ByvCg8Px82bNzFo0CA8e/ZMs66pqalo164djh8/rjlczMrKCmfPnsWjR4/ynVfuN1z79u1DWlpagdcTAIyNjTW3U1NTERcXh6ZNm0IQBM0hWU+fPsXx48cxfPhwrdcFAE0/q9Vq7Nq1C926dUPDhg3zLCd3uq1bt6JFixawtrbW6uP27dtDpVLh+PHj+eYUBAHbt29Ht27dIAiC1nMDAwORmJiodUgmAAQFBWmtH5Bz3mLuNxwqlQrPnj3THEb16vML6t9//4VMJsPEiRO12j/55BMIgoD//vtPq719+/aoXr265n69evVgYWHxxtcLAPTq1QsGBgZa32RcvnwZV69eRf/+/TVthd3GvXv3hr29vVablZUVHj58iPPnzxdsI7yBlZUVAGD37t2FGkghdz98dbu++k3nm3h5eaFJkyaa+7nfFrVt21brtZzbntsH8fHxOHz4MPr164fk5GTNNnz27BkCAwNx8+ZNREdHay1r9OjRWr/3WrRoAZVKhfv37wP4/+3w999/Q6lUFii/SqXC/v370aNHD7i7u2vaK1WqhEGDBuHkyZNISkoqVI78bN26FbVr14anp6fWa6Zt27YAoPndW5C+1Odr52369u2r9Q1/bj++//77Wuf1BgQEICsrS9NnBw4cQEJCAgYOHKi1vjKZDAEBAXn+1hCVFyzAiEqhp0+fIj09HR4eHnkeq1WrVoHn8+qb9Nw/kC4uLnnaXz636+bNm0hMTISDgwPs7e21flJSUvDkyZM3LgcArK2tteb5zTffICEhATVr1oS3tzc+++wzXLp06Y3579+/D6lUiho1ami1Ozk5wcrKKs8bmYLkeFVuQRgUFJRnXX/55RdkZmZqzt9asGABLl++DBcXF/j7+2PmzJlab9bd3NwwefJk/PLLL7Czs0NgYCCWLVv21vO/ACAqKgpDhw6FjY0NzMzMYG9vj1atWgGA5vm5y6pbt+5r5/P06VMkJSW9cZrc9d67d2+edW7fvj0A5Onjl+efkJCAVatW5XnusGHD8n2um5tbnvmo1WosWbIEHh4eUCgUsLOzg729PS5dulSg7ZWf+/fvo3LlyjA3N9dqr127tubxlxXl9QIAdnZ2aNeuHbZs2aJp+/PPP2FgYIBevXpp2gq7jfPbTp9//jnMzMzg7+8PDw8PjB8/vkCHDOanf//+aNasGUaOHAlHR0cMGDAAW7ZseWsxlrsfvlysAsX3uwiApg9u3boFQRDw9ddf59mOuYeAvu33kbW1tdY8W7Vqhd69e2PWrFmws7ND9+7dsXbt2jznWL3s6dOnSEtLy3eda9euDbVajQcPHhQqR35u3ryJK1eu5FnXmjVraq1rQfpSn6+dtylq/+b+/m3btm2edd6/f/9rfw8RlXU8B4yoHHvdaGT5tQsvDVKgVqvh4OCAjRs35vv8Vz+lf91yXp5ny5Ytcfv2bezevRv79+/HL7/8giVLlmDFihUYOXLkG9ejoBdnLkiOV+W+Yfn+++/zDLmcy8zMDADQr18/tGjRAjt37sT+/fvx/fffY/78+dixYwc6deoEAFi0aBGGDh2qWc+JEydi7ty5OHPmDJydnfOdv0qlQocOHRAfH4/PP/8cnp6eMDU1RXR0NIYOHVoswz6r1Wp06NABU6ZMyffx3Dd8+T0PyPlkOygoKN9p6tWrp3X/1W+/AOC7777D119/jeHDh2P27NmwsbGBVCrFpEmTSmyY66K8XnINGDAAw4YNQ3h4OHx9fbFlyxa0a9cOdnZ2mmkKu43z2061a9fG9evX8ffff2Pv3r3Yvn07li9fjunTp2sGAnnd/qFSqfLM//jx4zhy5Aj++ecf7N27F3/++Sfatm2L/fv3F+vohYX5XQT8fx/kvhY+/fTT1377/+oHNG+bp0QiwbZt23DmzBns2bMH+/btw/Dhw7Fo0SKcOXNGs7/rqqi/j7y9vbF48eJ8H88taArSlwV57eiLrv27YcMGODk55ZmOo+JSecVXNlEplDu6XX6H612/fr3Yl1+9enUcPHgQzZo1y/dNYVHZ2Nhg2LBhGDZsGFJSUtCyZUvMnDnztQVYtWrVoFarcfPmTc23GEDOIBAJCQmoVq2azplyP9W3sLDQfDPxJpUqVcK4ceMwbtw4PHnyBA0aNMC3336rKcAAwNvbG97e3vjqq68QHByMZs2aYcWKFZgzZ06+84yIiMCNGzewfv16rYFQXh4hEoDm0KfLly+/Np+9vT0sLCzeOA2Qs94pKSkFWudX529ubg6VSlXo575s27ZtaNOmDX799Vet9oSEBK0ipqDFN5Dzejl48CCSk5O1vgWLjIzUPK4vPXr0wJgxYzSHId64cQPTpk3Tmqao2/hVpqam6N+/P/r374+srCz06tUL3377LaZNmwYjIyNYW1sjISEhz/PyO9RNKpWiXbt2aNeuHRYvXozvvvsOX375JY4cOfLanLn74e3bt7W+ASqJ30W5r3lDQ0Odt+OrGjdujMaNG+Pbb7/Fpk2b8N577+GPP/7I9/eRvb09TExM8l3nyMhISKXSPN/2FEX16tVx8eJFtGvX7q2v/YL05dteO/kpzD6nq9zfvw4ODnrvX6LSjIcgEpVCMpkMgYGB2LVrF6KiojTt165dw759+4p9+f369YNKpcLs2bPzPJadnZ3vm723efbsmdZ9MzMz1KhR442H/XTu3BkA8owWl/vpcGFGYXsdPz8/VK9eHQsXLkRKSkqex3OHjVapVHkOjXNwcEDlypU165CUlITs7Gytaby9vSGVSt+4nrmfEr/8ybggCPjhhx+0prO3t0fLli2xZs0ardfFy8+VSqXo0aMH9uzZgwsXLuRZVu50/fr1w+nTp/N9PSUkJORZj5ez9u7dG9u3b8+3yCvIMNu583n1m4CtW7fmOZ/H1NRUk+ltOnfuDJVKpTUMPAAsWbIEEolEq0jWlZWVFQIDA7Flyxb88ccfkMvl6NGjh9Y0Rd3GL3t1v5HL5fDy8oIgCJrzl6pXr47ExEStQ3ofP36cZ5TQ/Ib/zv3W902vz9zt9r///U+r/U0jZeqLg4MDWrdujZUrV+Lx48d5Hi/o6+1lz58/z/Pae9t2kMlk6NixI3bv3q01xH5sbCw2bdqE5s2bw8LCotBZXtWvXz9ER0dj9erVeR5LT09HamoqgIL1ZUFeO/kpzD6nq8DAQFhYWOC7777LN1NR+peoLOA3YESl1KxZs7B37160aNEC48aNQ3Z2tuaaLm87d0pXrVq1wpgxYzB37lyEh4ejY8eOMDQ0xM2bN7F161b88MMP6NOnT6Hm6eXlhdatW8PPzw82Nja4cOGCZkj31/Hx8UFQUBBWrVqFhIQEtGrVCufOncP69evRo0cPtGnTRtdVhVQqxS+//IJOnTqhTp06GDZsGKpUqYLo6GgcOXIEFhYW2LNnD5KTk+Hs7Iw+ffrAx8cHZmZmOHjwIM6fP49FixYByBnOfsKECejbty9q1qyJ7OxsbNiwQVO0vI6npyeqV6+OTz/9FNHR0bCwsMD27dvzPVfkf//7H5o3b44GDRpg9OjRcHNzw7179/DPP/8gPDwcQM7hffv370erVq00Q1k/fvwYW7duxcmTJ2FlZYXPPvsMf/31F7p27aoZqj81NRURERHYtm0b7t27p/VN1MvmzZuHI0eOICAgAKNGjYKXlxfi4+MRGhqKgwcP5vvm8FVdu3bFN998g2HDhqFp06aIiIjAxo0btQY4AHKKCysrK6xYsQLm5uYwNTVFQEBAvudLdevWDW3atMGXX36Je/fuwcfHB/v378fu3bsxadKkPOcw6ap///54//33sXz5cgQGBmoGRsilyzbO1bFjRzg5OaFZs2ZwdHTEtWvX8NNPP6FLly6ab/kGDBiAzz//HD179sTEiRM1w3jXrFlTa0CTb775BsePH0eXLl1QrVo1PHnyBMuXL4ezszOaN2/+2gy+vr4YOHAgli9fjsTERDRt2hSHDh3CrVu3ir7xCmHZsmVo3rw5vL29MWrUKLi7uyM2NhanT5/Gw4cPta4dVxDr16/H8uXL0bNnT1SvXh3JyclYvXo1LCwsNB/65GfOnDmaa2+NGzcOBgYGWLlyJTIzM7FgwQJdVxMAMHjwYGzZsgUffPABjhw5gmbNmkGlUiEyMhJbtmzBvn370LBhwwL1ZUFeO/kpzD6nKwsLC/z8888YPHgwGjRogAEDBsDe3h5RUVH4559/0KxZszwfqBCVCyU55CIRFc6xY8cEPz8/QS6XC+7u7sKKFSvyDC0tCK8fhv7VYZZzn/vq0MVBQUGCqalpnuWvWrVK8PPzE4yNjQVzc3PB29tbmDJlivDo0SOtZec3vPyrwyvPmTNH8Pf3F6ysrARjY2PB09NT+Pbbb4WsrKw8+V6mVCqFWbNmCW5uboKhoaHg4uIiTJs2TWtI6sLkeHUY+lxhYWFCr169BFtbW0GhUAjVqlUT+vXrJxw6dEgQBEHIzMwUPvvsM8HHx0cwNzcXTE1NBR8fH2H58uWaedy5c0cYPny4UL16dcHIyEiwsbER2rRpIxw8eDBPrlddvXpVaN++vWBmZibY2dkJo0aN0gyL/upw0JcvXxZ69uwpWFlZCUZGRkKtWrWEr7/+Wmua+/fvC0OGDBHs7e0FhUIhuLu7C+PHjxcyMzM10yQnJwvTpk0TatSoIcjlcsHOzk5o2rSpsHDhQq1+wSvDcQuCIMTGxgrjx48XXFxcBENDQ8HJyUlo166dsGrVqjzbOr9h+TMyMoRPPvlEqFSpkmBsbCw0a9ZMOH36dJ7+EgRB2L17t+Dl5SUYGBhobY/8hl9PTk4WPv74Y6Fy5cqCoaGh4OHhIXz//fdaw4HnrlN+lxR4dV96k6SkJMHY2FgAIPz+++/5TlOQbZw7lPf333+f5/krV64UWrZsqXldVq9eXfjss8+ExMREren2798v1K1bV5DL5UKtWrWE33//Pc/+dOjQIaF79+5C5cqVBblcLlSuXFkYOHBgnmHy85Oeni5MnDhRsLW1FUxNTYVu3boJDx48KPAw9Pntm/n1weu2xe3bt4UhQ4YITk5OgqGhoVClShWha9euwrZt2/Is+9Xfe6/u86GhocLAgQOFqlWrCgqFQnBwcBC6du0qXLhwIU++V1/3oaGhQmBgoGBmZiaYmJgIbdq00RpevzA5BCHv7ydBEISsrCxh/vz5Qp06dQSFQiFYW1sLfn5+wqxZszT9XpC+LOhrJz8F3ede11+v2/fftG0CAwMFS0tLwcjISKhevbowdOjQPH1CVF5IBKEAZxsTERERERGRzngOGBERERERUQlhAUZERERERFRCWIARERERERGVEBZgREREREREJYQFGBERERERUQnhdcBeoVar8ejRI5ibm5fo1eCJiIiIiKh0EQQBycnJqFy5MqRS/Xx3xQLsFY8ePYKLi4vYMYiIiIiIqJR48OABnJ2d9TIvFmCvyL06/N27d2FjYyNymopHqVRi//796NixIwwNDcWOUyGxD8THPhAf+0B87APxsQ/Exz4QX3x8PNzc3DQ1gj6wAHtF7mGH5ubmsLCwEDlNxaNUKmFiYgILCwv+ohEJ+0B87APxsQ/Exz4QH/tAfOwD8SmVSgDQ66lJHISDiIiIiIiohLAAIyIiIiIiKiEswIiIiIiIiEoICzAiIiIiIqISwgKMiIiIiIiohLAAIyIiIiIiKiEswIhIQxAEJKUrkZAJZChVYschIiIiKnd4HTCiCkStFnDlURLO34tHbHIG4pKz8Cw1E3EpmXiWkoVnKVnIUqkBGGBG6CGYymWwMZPD1lQBuxf/25rJ4WJjggA3G7jZmer1uhhERERE5R0LMKJyLupZGk7eisOpW3E4dTsOCWnKtz5HKhGgFiRIzVIhNT4dD+LT853OwVyBxu62L35YkBERERG9DQswonImW6XG4cgnOHL9KU7dikNUfJrW42YKAwS42cDVzhS2ZnLYmeV8u2VnpoCtmQIWcgkO7d+Hlu06IilT/eIbspxvx+Jf3L72OAlhUQl4kpyJvy4+wl8XHwEAHC1yCrL2tR0RWMcJcgMe5UxERET0MhZgROVESmY2/jgXhbWn7iE64f+/sTKQStCgqjWa1bBDcw9b+DhbwUD2+sJIqVRCIgHMjQxgY24IVzvTfKfLUKoQFpWAM3ee4cydZwiLSkBsUiZ2hz/C7vBHsDWVo29DFwzyr4qqtiZ6X18iIiKisogFGFEZ9yghHeuC72Hz2SgkZ2YDAGxM5XjXpzJa1rSDv5stzBT639WNDGVoUt0WTarbAsgpyEKjnuPkzThsD32I2KRMrDh2GyuO3UYLDzu8F1AN7Ws7vLH4IyIiIirvWIARlVERDxOx+sQd/BPxGCq1AABwtzfFyObu6NWgCowMZSWax8hQhqbV7dC0uh0md6iJQ5FPsPFsFI7feIoTN+Nw4mYcHC0U6N+oKoY0qQY7M0WJ5iMiIiIqDViAEZUx9+JSMf2vKzh+46mmrYm7LUa1dEPrmg6QSsUfBMNAJkVgHScE1nFC1LM0bD4fhS3nHyA2KRP/O3QTa07exQet3DGiuTuM5SVbKBIRERGJiQUYURmRla3G6hN38L9DN5GZrYaBVIKu9SphZAt31K1iKXa816pqa4LP3/HEpPYe2H8lFquO30FEdCIW7r+BDWfuY3KHmujj5wJZKSgciYiIiIobCzCiMiDkfjy+2HEZ12OTAQDNa9hhTo+6rx0gozRSGMjQzacyunhXwt8Rj7FgbyQePk/H59sj8OvJu5jayRNtajlwGHsiIiIq11iAEZViSRlKLNgbiY1noyAIOYNrfN21Nnr4VimzhYpUKsG7PpURWMcRG07fx4+Hb+FGbAqGr7uAJu62mNbZE/WcrcSOSURERFQsWIARlUKCIOC/yzGY+dcVPEnOBAD09XPGF51rw9pULnI6/VAYyDCyhTv6+rlg+dFbWBt8D6fvPMO7P53C0Kau+PwdT54fRkREROUOCzCiUiY9S4WpOy5hd3jOxY3d7Ezxbc+6aFrdTuRkxcPSxBDTOtfG4CbVsHDfdewKf4R1wfdw7MZTLOxbD37VbMSOSERERKQ3vCAPUSny8Hkaev8cjN3hj2AglWBi2xr476MW5bb4epmztQmWDqiPdcMawcnCCHfjUtF3xWnM+y8SmdkqseMRERER6QULMKJS4vTtnMPvrj5Ogq2pHBtHBmByx1olfj0vsbWu5YB9k1qiV4MqUAvAimO30e3Hk4h4mCh2NCIiIiKdsQAjEpkgCFgffA/v/3oW8alZqFvFAns+bI4Ad1uxo4nG0sQQi/v5YtVgP9iZyXEjNgU9lp/CkgM3oFSpxY5HREREVGQswIhElJmtwufbL2HGX1egUgvo4VsZ2z5oispWxmJHKxU61nHC/o9boUu9SlCpBfxw6CZ6Lj+FB/FpYkcjIiIiKhIWYEQiiU3KQP+VZ7DlwkNIJcCXnWtjSX/fCnfI4dvYmMqxbFAD/DiwPqxMDHE5OgldfzyJ4zeeih2NiIiIqNBYgBGJ4HJ0Irr9eBLhDxJgaWyI9cP9Maqle5m9tldJ6OZTGf991AI+LlZITFciaO05LDtyC4IgiB2NiIiIqMBYgBGVsIsPEjBo9Rk8Sc5ETUcz/DWhGVp42Isdq0yoZGmMLWMaY0AjFwgC8P2+6xj7eyhSMrPFjkZERERUICzAiEpQyP3neP+Xs0jKyIZfNWtsH9sU1WxNxY5VpigMZJjXux7m9vKGXCbF3isx6LHsFG4/TRE7GhEREdFbsQAjKiHn78VjyK9nkZyZDX83G/w23B/mRoZixyqzBvpXxZ9jGsPJwgi3nqSg+0+nsP9KjNixiIiIiN6IBRhRCTh9+xmC1pxDapYKTavbYt2wRjBVGIgdq8yrX9Uaez5sDn83G6RkZmP0hhAsOXCD54URERFRqcUCjKiYnbwZh2HrziEtS4UWHnZYM7QRTOQsvvTF3lyBjSMDMKyZKwDgh0M3MWXbJV4vjIiIiEolvRdgSUlJOH78OP73v//pe9ZEZc6xG08xYv15ZCjVaFPLHquHNOQw88XAUCbFjG51MK+XN2RSCbaGPMTI9ReQysE5iIiIqJTR6WP4qKgohIeHa/3cv38fgiDA1NQUEydO1FdOojLn0LVYjP09FFkqNdrXdsSy9+pDYcDiqzgN8K8Ke3MFxm8KxbEbTzFw9RmsGdoIdmYKsaMRERERASjiN2Bt27aFra0tXF1dERQUhH379sHe3h5RUVH49ddfcf/+fSQnJ+s7K1GZEXwrDh/8HoIslRrv1HHC8vcasPgqIe1qO2LzqMawMZXj0sNE9P45GPfiUsWORURERASgiAXYyZMn8cEHH+DBgwd4/vw5Tp06hZUrV0IikcDf3x8uLi76zklUZlyPScaYDSFQqgR0quuEHwfVh9yAp1uWpPpVc4b4r2pjgvvP0tD752CEP0gQOxYRERFR0Qqws2fP4sSJExg/fjxu3Lih70xEZVZMYgaGrj2nGWp+SX9fGMpYfInBzc4U28c2hXcVSzxLzcLAVWdwODJW7FhERERUwRXpnWH9+vVx/Phx9OvXD4GBgRg/fjyePHmi72xEZUpyhhLD1p3H48QMVLc3xarBfhxwQ2T25gr8MboxWta0R7pShVG/hWBbyEOxYxEREVEFptNH84MGDcKVK1dgbW2NOnXqQK1WQ6VS6SsbUZmhVKkxbmMorj1Ogp2ZAuuG+cPKRC52LAJgqjDAr0EN0buBM1RqAZ9tu4jN56LEjkVEREQVlM7HRpmYmGDOnDk4e/Ysunbtinbt2mHhwoVIT0/XRz6iUk8QBHyxIwInbsbBRC7D2qGN4GJjInYseomhTIqFfethaFNXCAIwbUcEfjt9T+xYREREVAHp7eQUd3d37N69G7///jvWrl0Ld3d3fc2aqFT74dBNbA15CKkEWDaoAbydLcWORPmQSCSY0c0Lo1vm/G6avvsKfjlxR+RUREREVNEUqQCbPn06QkJC8n2sY8eOuHjxIj7//HOdghGVBVsvPMDSgzcBAHN6eKONp4PIiehNJBIJpnXyxIQ2NQAAc/65huVHb4mcioiIiCqSIhVgDx8+RKdOneDs7IyxY8fiv//+Q1ZWluZxAwMDTJo0SV8ZiUql4zeeYtqOCADA+DbVMSigqsiJqCAkEgk+DayFj9vXBAAs2HsdPxy8CUEQRE5GREREFUGRCrA1a9YgJiYGmzdvhrm5OSZNmgQ7Ozv07t0bv/32G+Lj4/Wdk6hUuReXivEbQ5GtFtDDtzI+7VhL7EhUSB+198CUd3L6bcnBG1i4/zqLMCIiIip2RT4HTCqVokWLFliwYAGuX7+Os2fPIiAgACtXrkTlypXRsmVLLFy4ENHR0frMSyS69CwVPvg9BMmZ2WhYzRoL+vhAIpGIHYuKYFzrGviqS20AwLIjtzH3v0gWYURERFSsdBqEIzk5WXO7du3amDJlCk6dOoUHDx4gKCgIJ06cwObNm3UKuGzZMri6usLIyAgBAQE4d+7ca6fdsWMHGjZsCCsrK5iamsLX1xcbNmzQaflELxMEAV/ujEBkTDLszBRY9l4DyA14oeWybGQLd3zTvQ4AYNXxO5jHIoyIiIiKkU7vHFu0aIGYmJg87fb29hgxYgR2796NTz/9tMjz//PPPzF58mTMmDEDoaGh8PHxQWBg4Gsv+mxjY4Mvv/wSp0+fxqVLlzBs2DAMGzYM+/btK3IGopf9fjYKO8KiIZNK8NOg+nC0MBI7EunBkCau+LZnXQDAyuN38MOhmyInIiIiovJKpwKsfv36CAgIQGRkpFZ7eHg4OnfurFMwAFi8eDFGjRqFYcOGwcvLCytWrICJiQnWrFmT7/StW7dGz549Ubt2bVSvXh0fffQR6tWrh5MnT+qchSgs6jm+2XMFAPD5O7XQ2N1W5ESkT+8FVMP0rl4AgKUHb2LFsdsiJyIiIqLyyECXJ69duxYzZsxA8+bNsWvXLjg4OOCrr77C9u3bdS7AsrKyEBISgmnTpmnapFIp2rdvj9OnT7/1+YIg4PDhw7h+/Trmz5//2ukyMzORmZmpuZ+UlAQAUCqVUCqVOqwBFUXuNi9t2/5ZahbG/h4CpUpARy8HDG3sUuoy6ktp7YOSMDjAGakZWVh08Bbm/RcJuRQY3LjkR7esyH1QWrAPxMc+EB/7QHzsA/EVx7aXCHo42eG7777D7NmzoVKp0K5dO8yaNQv+/v46zfPRo0eoUqUKgoOD0aRJE037lClTcOzYMZw9ezbf5yUmJqJKlSrIzMyETCbD8uXLMXz48NcuZ+bMmZg1a1ae9k2bNsHExESndaDyQS0AP1+T4kaiFA5GAj7xVsFIp48uqLT7J0qK/dE5BwgMrK5CYweeE0ZERFQRpaWlYdCgQUhMTISFhYVe5qnT28jY2Fh89913WL16Nby8vBAZGYmhQ4fqXHzpwtzcHOHh4UhJScGhQ4cwefJkuLu7o3Xr1vlOP23aNEyePFlzPykpCS4uLmjTpg1sbXmIWUlTKpU4cOAAOnToAENDQ7HjAAAWH7yJG4l3YWwoxbqRjeHhaCZ2pGJVGvugpHUSBMzdewNrg+/jjzsyNGrgjW71KpXY8tkH4mMfiI99ID72gfjYB+J79uyZ3uepUwHm5uaGWrVqYevWrejSpQv27t2L/v37IyoqCp999plOwezs7CCTyRAbG6vVHhsbCycnp9c+TyqVokaNGgAAX19fXLt2DXPnzn1tAaZQKKBQKPK0Gxoa8oUuotKy/Q9ejcXPx+4CAOb1rgcvZ2uRE5Wc0tIHYpnerQ6yVAI2no3CZ9svw0Qhxzt1X/+7pzhU9D4oDdgH4mMfiI99ID72gXiKY7vrNAjHmjVrEBYWhi5dugAA3nnnHRw5cgRLlizB+PHjdQoml8vh5+eHQ4cOadrUajUOHTqkdUji26jVaq1zvIgK6v6zVHy8JRwAENSkGrr7VhE3EJUoiUSC2d3roleDKlCpBXy4ORRHr+c/AisRERFRQelUgA0YMCBPW4MGDRAcHIzDhw/rMmsAwOTJk7F69WqsX78e165dw9ixY5Gamophw4YBAIYMGaI1SMfcuXNx4MAB3LlzB9euXcOiRYuwYcMGvP/++zpnoYpFqVLjw81hSM7IRoOqVviyi5fYkUgEUqkEC3rXQxfvSlCqBIzZEIJzd+PFjkVERERlWLEMJeDq6org4GCd59O/f388ffoU06dPR0xMDHx9fbF37144OjoCAKKioiCV/n8NmZqainHjxuHhw4cwNjaGp6cnfv/9d/Tv31/nLFSxLD14A5ceJsLS2JAXW67gDGRSLB3giwylCocin2DE+vPYMqYJalfSz4m4REREVLEU21hu1tb6OVdmwoQJmDBhQr6PHT16VOv+nDlzMGfOHL0slyquc3fjsfxozjWg5vbyRiVLY5ETkdgMZVL8NKgBBv96FhfuP0fQmnPYPrYpXGw4UioREREVDj/WJ3pJUoYSH/8ZDkEAejdwRmfvkhv5jko3Y7kMvwY1Qi1HczxJzsTgX88iLoXnlxIREVHhFKoAS09PR3R0dJ72K1eu6C0QkZhm7r6C6IR0uNgYY+a7PO+LtFmaGOK3Ef6oYmWMe8/SMHTtOSRn8OKYREREVHAFLsC2bdsGDw8PdOnSBfXq1dO6EPLgwYOLJRxRSdpz8RF2hEVDKgGW9POFuRGHe6W8HC2MsGGEP2xN5bgcnYQxG0KQma0SOxYRERGVEQUuwObMmYOQkBCEh4dj7dq1GDFiBDZt2gQAEASh2AISlYTHien4cmcEAGBCmxpo6GojciIqzdztzbBumD9M5TIE336Gj/8Mh0rN34NERET0dgUehEOpVGpGH/Tz88Px48fRs2dP3Lp1CxKJpNgCEhU3tVrAJ1suIikjGz7OlviwnYfYkagM8Ha2xOohDTF07Xn8GxEDK5PL+LZHXf4+JCIiojcq8DdgDg4OuHTpkua+jY0NDhw4gGvXrmm1E5U1v568i+Dbz2BsKMOS/r4wlHFsGiqYpjXssHSALyQSYNPZKCw5eFPsSERERFTKFfid5oYNG+Dg4KDVJpfLsXnzZhw7dkzvwYhKwtVHSfh+33UAwPRuXnC3NxM5EZU1nb0rYXb3ugCA/x26ic3nokRORERERKVZgQswZ2dnODk55ftYs2bN9BaIqKRkKFWY9GcYslRqdPByxIBGLmJHojLq/cbVMPHFoatf7bqMw5GxIiciIiKi0qpIx1pdvXoVX3/9NRISEvQch6jkLD5wAzdiU2BnpsC8Xt48d4d08nF7D/Txc4ZKLWD8xjBcfJAgdiQiIiIqhYpUgM2dOxeXL1+GlZVVnscyMjIQGRmpay6iYnXxQQJ+OXEHADC/tzdszRQiJ6KyTiKRYG4vb7TwsEO6UoXh687j/rNUsWMRERFRKVOkAuzMmTOYOHFivo8ZGRlh1KhRmDt3rk7BiIpLZrYKn227CLUA9PCtjHa1HcWOROWEoUyKn9/3Q53KFniWmoWha88jPjVL7FhERERUihSpAHv48CFq1Kjx2sc/+OAD/PXXX0UORVSclh25jRuxKbA1lWN6tzpix6FyxkxhgLVDG6GKlTHuxqVixPrzSM/ihZqJiIgoR5EKMBsbGzx+/Pi1j/v7++PWrVtFDkVUXK49TsLyIzmvzW+614WNqVzkRFQeOVgYYf3wRrA0NkRYVAI++iOMF2omIiIiAEUswFq2bIl169a9fqZSKTIyMoqaiahYZKvUmLLtErLVAgLrOKKzd/6jehLpQw0Hc/wS1BByAyn2X43FrD1XIAgswoiIiCq6IhVgn376KVavXo1Vq1bl+/jp06fh7u6uUzAifVt94i4iohNhYWSA2d3rctRDKnaNXG2wtH/OhZp/O30fq47fETsSERERiaxIBZifnx+WL1+OcePGoUOHDti1axeioqIQHx+P3bt34/PPP8egQYP0nZWoyG4/TcGSgzcAAF939YKDhZHIiaii6OxdCV918QIAzP0vEv9cev3h20RERFT+FakAA4BRo0bh6NGjSExMRK9eveDm5gZ7e3v07NkT3t7e+Pjjj/WZk6jI1GoBn2+7hKxsNVrWtEcfP2exI1EFM6K5G4Y2dQUAfLwlHCH348UNRERERKIx0OXJzZs3x7lz5xAZGYnQ0FCkpaWhbt26aNy4sb7yEenst9P3cOH+c5jKZfiuJw89JHF83dULD5+n4+C1WIz6LQQ7xjaFq52p2LGIiIiohOlUgOXy9PSEp6enPmZFpFcP4tOwYN91AMDUTp5wtjYRORFVVDKpBP8b6IsBq87g0sNEDFt3HjvGNoU1R+IkIiKqUIp8CCJRaScIAqbtiEBalgr+bjZ4L6Ca2JGogjORG+CXoIaaa4SN3nABGUpeI4yIiKgiYQFG5daO0GicvBUHhYEU83vXg1TKQw9JfA7mRlg3rBHMjQxw/t5zfLbtEtS8RhgREVGFoXMB9vDhQ6jV6jy3icSUkJaFb/+9BgCY2M4DbjzXhkoRD0dzrHzfDwZSCfZcfIRFB66LHYmIiIhKiM4FmJeXF+7du5fnNpGY5u+NRHxqFjwczDCqBa9JR6VP0xp2mNe7HgBg2ZHb+ONclMiJiIiIqCToXIAJgpDvbSKxhNyPx+ZzDwAAc3rUhdyAR9pS6dTHzxkftfMAAHy56zKO33gqciIiIiIqbnxnSuWKUqXGlzsvAwD6+jkjwN1W5EREbzapvQd61a8ClVrA+I2huB6TLHYkIiIiKkYswKhcWXvqLiJjkmFtYohpnWuLHYforSQSCeb29kaAmw2SM7MxfN15PE3OFDsWERERFRMWYFRuRCekY8mBmwCAaZ1rw4bXV6IyQmEgw8rBfnC3M0V0QjrGbAxDFkenJyIiKpdYgFG5MWP3FaQrVfB3tUGfBs5ixyEqFCsTOdYOawQbUzkiopPw200pVByenoiIqNxhAUblwv4rMTh4LRYGUgnm9KzLa35RmVTN1hSrBvvBUCZBxHMpFuy7IXYkIiIi0jMWYFTmpWZmY+ZfVwAAo1q6o6ajuciJiIquoasN5veqCwBYE3wfG07fEzcQERER6ZXOBdgXX3wBGxubPLeJSsoPh27iUWIGnK2NMbGth9hxiHTWrV4ldHHJOQlsxl9XcCTyiciJiIiISF90LsCmTZsGKyurPLeJSsK1x0n49eRdAMA33evAWC4TORGRfnSoIqBX/cpQC8CETaG4+ihJ7EhERESkBzwEkcostVrAV7suQ6UW0KmuE9p6OoodiUhvJBJg9rteaOJui9QsFYavO4+YxAyxYxEREZGOWIBRmbUjLBoh95/DRC7D9G5eYsch0ju5gRQr3vdDdXtTxCRlYMT680jNzBY7FhEREemABRiVSYnpSsz77xoA4KN2HqhkaSxyIqLiYWliiLVD/WFrKseVR0n4cHMYh6cnIiIqw1iAUZm09OANxKVkwd3eFMOauYkdh6hYVbU1weqghlAYSHE48glm/31V7EhERERURDoVYMnJyfrKQVRgkTFJ+O30fQDAzG51IDfg5whU/jWoao0l/X0BAOuC72HNi8FniIiIqGzR6Z1rixYtEBMTU6BpY2JioFardVkcEQRBwIzdV6BSC3injhNa1rQXOxJRiensXQlTO3kCAGb/cxX7rxTs9y8RERGVHjoVYPXr10dAQAAiIyO12sPDw9G5c2ettkOHDmH+/Pm6LI4Iey49xtm78VAYSPFV19pixyEqcWNaumOgf1UIAvDRH+GIeJgodiQiIiIqBJ0KsLVr12Lo0KFo3rw5Tp48iRs3bqBfv37w8/ODTKZ9PaYePXpgz549OoWlii01Mxvf/pNz7sv4NjXgbG0iciKikieRSPBN9zpo4WGHdKUKw9efR3RCutixiIiIqIB0Pnlm1qxZmDx5Mjp06IC6desiOTkZp0+fzlNsXb16FS4uLroujiqwHw/fQmxSJqramGB0S3ex4xCJxlAmxfL3GsDTyRxPkzMxfO15JGUoxY5FREREBaBTARYbG4uPPvoIc+bMgZeXFwwNDTF06FD4+/vnmXbnzp0YNGiQLoujCuz20xT8evIOAGBGNy8YGcre8gyi8s3cyBC/Dm0Ee3MFrscmY/zGUChVPM+WiIiotNOpAHNzc8Px48exdetWhISEYPv27Rg9ejS+//77PNMmJSXByclJl8VRBSUIAmb+dQVKlYC2ng5oV9tR7EhEpUIVK2OsCWoEY0MZTtyMw5c7IyAIvEYYERFRaaZTAbZmzRqEhYWhS5cuAIB33nkHR44cwZIlSzB+/HitaQcPHoz169frsjiqoPZfjcWJm3GQy6SY3tVL7DhEpYq3syV+GlQfUgmw5cJD/HT4ltiRiIiI6A10KsAGDBiQp61BgwYIDg7G4cOHtdoDAgKwfPlyXRZHFVCGUoVv9uQMvDG6pTtc7UxFTkRU+rSr7YhZ79YBACw6cAM7wx6KnIiIiIhep1iuYOvq6org4ODimDVVMD8fvY3ohHRUtjTCuDbVxY5DVGoNbuKqGZxmyrZLCL4dJ3IiIiIiyk+xFGAAYG1trZf5LFu2DK6urjAyMkJAQADOnTv32mlXr16NFi1awNraGtbW1mjfvv0bp6fS7eHzNKw4dhsA8FVXL5jIDURORFS6TX3HE128K0GpEjBmQwhuxiaLHYmIiIheUWwFmD78+eefmDx5MmbMmIHQ0FD4+PggMDAQT548yXf6o0ePYuDAgThy5AhOnz4NFxcXdOzYEdHR0SWcnPThu3+vITNbjSbutuhUlwO4EL2NVCrBon4+8KtmjeSMbAxdex5PkjLEjkVEREQvKdUF2OLFizFq1CgMGzYMXl5eWLFiBUxMTLBmzZp8p9+4cSPGjRsHX19feHp64pdffoFarcahQ4dKODnp6uzdePwbEQOpBJjxrhckEonYkYjKBCNDGVYPaQg3O1NEJ6Rj+PrzSM3MFjsWERERvVBqj+nKyspCSEgIpk2bpmmTSqVo3749Tp8+XaB5pKWlQalUwsbG5rXTZGZmIjMzU3M/KSkJAKBUKqFU8sKmJU2pVEIlALP/uQYAGNjIBdVtjdkXJSh3W3Obi0fXPjCXS7D6/frou+osLkcnYfzGEPw8yBcGslL9mVupwv1AfOwD8bEPxMc+EF9xbHuJUEovGvPo0SNUqVIFwcHBaNKkiaZ9ypQpOHbsGM6ePfvWeYwbNw779u3DlStXYGRklO80M2fOxKxZs/K0b9q0CSYmJkVfASqykzESbL0rg4lMwFf1VTA1FDsRUdl0NxlYdkUGpSBBM0c1+rqpwS+TiYiICi4tLQ2DBg1CYmIiLCws9DLPUvsNmK7mzZuHP/74A0ePHn1t8QUA06ZNw+TJkzX3k5KS4OLigjZt2sDW1rYkotJLnialYdqSEwCAzzrVRt+AqiInqniUSiUOHDiADh06wNCQ1a8Y9NkH1a/EYuKfF3EqVooA75oY28pdTynLN+4H4mMfiI99ID72gfiePXum93nqXIAdOnQIhw4dwpMnT6BWq7Uee925WgVhZ2cHmUyG2NhYrfbY2Fg4Ob15QIaFCxdi3rx5OHjwIOrVq/fGaRUKBRQKRZ52Q0NDvtBF8PPx+0jLlqCWoxkGN3HjIVMi4j4gPn30QTdfZzxLVWLmnqtYfPAWKlmZoG9DFz0lLP+4H4iPfSA+9oH42AfiKY7trtO721mzZqFjx444dOgQ4uLi8Pz5c60fXcjlcvj5+WkNoJE7oMbLhyS+asGCBZg9ezb27t2Lhg0b6pSBStb1mGRsOp9zAdmvOnuy+CLSk6HN3DDmxTdfU3dE4Oj1/EeSJSIiouKn0zdgK1aswLp16zB48GB95dEyefJkBAUFoWHDhvD398fSpUuRmpqKYcOGAQCGDBmCKlWqYO7cuQCA+fPnY/r06di0aRNcXV0RExMDADAzM4OZmVmxZCT9EAQBs/ZcgUotwMdGjcburx84hYgK7/NAT8QmZmBX+COM2xiKP0c3gbezpdixiIiIKhydvmLIyspC06ZN9ZUlj/79+2PhwoWYPn06fH19ER4ejr1798LR0REAEBUVhcePH2um//nnn5GVlYU+ffqgUqVKmp+FCxcWW0bSj31XYhB8+xkUBlJ0r6Z++xOIqFCkUgkW9PFBsxq2SMtSYdi6c4h6liZ2LCIiogpHpwJs5MiR2LRpk76y5GvChAm4f/8+MjMzcfbsWQQEBGgeO3r0KNatW6e5f+/ePQiCkOdn5syZxZqRdJOhVGHOi2HnRzZ3he3rx0whIh3IDaRY8b4faleyQFxKFoLWnkN8apbYsYiIiCqUQh+C+PKIgWq1GqtWrdIMdvHqSWqLFy/WPSGVe7+cuIOHz9NRydIIo1u44ujBG2JHIiq3zI0MsW5YI/RaHoy7cakYvu48No9qDGO5TOxoREREFUKhC7CwsDCt+76+vgCAy5cva7VLeLEZKoDHielYduQ2AGBqJ0+YyMvtlRGISg1HCyOsH+6PPiuCEf4gAR9uDsWK9/048A0REVEJKPS73SNHjmhuR0VFwdnZGVKp9h9tQRDw4MED3dNRuTfvv0ikK1VoWM0a7/pURnZ2ttiRiCqEGg5m+GVIQ7z3y1kcvPYEX+68jHm9vfnhGRERUTHT6eNONzc3xMXF5WmPj4+Hm5ubLrOmCuDCvXjsDn8EiQSY+W4dvvEjKmENXW3w48D6kEqAPy88wIJ918WOREREVO7pVIAJgpBve0pKCoyMOJICvZ5aLWDWnqsAgP4NXVC3CofDJhJDxzpOmNvLGwDw89Hb+OXEHZETERERlW9FOuEmdyAOiUSC6dOnw8TERPOYSqXC2bNnNeeGEeVnW+hDREQnwlxhgE861hI7DlGF1r9RVcSnKjF/byTm/HMN1iZy9PZzFjsWERFRuVSkAix3IA5BEBAREQG5XK55TC6Xw8fHB59++ql+ElK5k5yhxIK9OYc6TWznAXtzhciJiOiDVu6IT83E6hN3MWX7JViZGKJdbUexYxEREZU7RSrAcgfiGDZsGH744QdYWFjoNRSVbz8duYW4lEy42ZkiqKmr2HGICDlHNEzrVBvPUrOwIzQa4zaG4veRAWjkaiN2NCIionJFp3PA1q5dy+KLCuVuXCrWnLwLAPi6a23IDTjsNVFpIZVKML93PbTzdEBmthrD153HtcdJYsciIiIqV4p0IebZs2fD1NRU66LM+eGFmOlV3/5zFUqVgFY17dGmloPYcYjoFYYyKX4a1ABD1pzF+XvPEbTmHLaPbQoXG5O3P5mIiIjeqkgXYlYqlZrbr8MhxelVx248xcFrT2AgleDrrl58jRCVUsZyGX4Z0gj9V51GZEwy3v/1LLaOaQIHC45uS0REpCudLsT88m2iN1Gq1Jj9d86w80FNXVHDwUzkRET0JpYmhlg/3B99VgTj/rM0vP/rWfw5ugmsTeVvfzIRERG9Fk/AoRLx+5n7uPUkBTamckxs5yF2HCIqAEcLI/w+IgAO5grciE3BkDXnkJShFDsWERFRmaZTATZkyBCsWbMGt2/f1lceKofiU7Ow5MANAMCnHWvB0thQ5EREVFDVbE2xcWQAbEzliIhOxIh155GWlS12LCIiojJLpwJMLpdj3rx58PDwgIuLC95//3388ssvuHnzpr7yUTmw+MB1JGVko3YlC/Rv5CJ2HCIqJA9Hc/w23B/mRgY4f+85xmwIQWa2SuxYREREZZJOBdgvv/yCGzdu4MGDB1iwYAHMzMywaNEieHp6wtnZWV8ZqQy7+igJm85GAQBmdPOCTMqBN4jKorpVLLFuWCOYyGU4cTMOEzaFQalSix2LiIiozNHLOWDW1tawtbWFtbU1rKysYGBgAHt7e33MmsowQRAwc88VqAWgi3clNHa3FTsSEenAr5oNVg9pCLmBFAeuxuLTrRehUgtixyIiIipTdCrAvvjiCzRt2hS2traYOnUqMjIyMHXqVMTExLxxiHqqGP6+9Bjn7sbDyFCKaZ09xY5DRHrQrIYdfn6vAQykEuwOf4SvdkVAEFiEERERFVShh6F/2bx582Bvb48ZM2agV69eqFmzpr5yURmXlpWN7/69BgAY26oGnK15EVei8qJdbUcs6e+Lj/4Iw+ZzD2BsaICvu9bmtf2IiIgKQKcCLCwsDMeOHcPRo0exaNEiyOVytGrVCq1bt0br1q1ZkFVgPx+9jceJGahiZYwxrdzFjkNEetbNpzLSs1SYsv0S1py6C5kU+KIzizAiIqK30akA8/HxgY+PDyZOnAgAuHjxIpYsWYLx48dDrVZDpeIoWRVR1LM0rDx+BwDwddfaMDKUiZyIiIpDv0YuyFSp8fWuy1h94i6kEgmmdvJkEUZERPQGOhVggiAgLCwMR48exdGjR3Hy5EkkJSWhXr16aNWqlb4yUhnz7b9XkZWtRrMatgis4yR2HCIqRoMbVwMEAV/vvoKVx+9AIpHg83dqsQgjIiJ6DZ0KMBsbG6SkpMDHxwetWrXCqFGj0KJFC1hZWekpHpU1J24+xb4rsZBJJZjRrQ7fhBFVAIObuEIAMH33Faw4dhtSCfBZIIswIiKi/OhUgP3+++9o0aIFLCws9JWHyjClSo1Ze64CAIY0qYaajuYiJyKikjKkiSvUagEz91zF8qO3IZVI8EnHmizCiIiIXqFTAdalSxd95aBy4LfT93HrSQpsTOWY1J4DsBBVNEObuUEAMGvPVfx05BYkEmByBxZhREREL9PLhZiJ4lIysfTADQDAlMBasDQ2FDkREYlhWDM3fN3VCwDw4+FbWHrwpsiJiIiIShcWYKQX3++9juTMbHhXsUTfhi5ixyEiEY1o7oavutQGAPxw6CYW77/OizUTERG9wAKMdHbxQQK2hDwAAMx81wsyKQ83IqroRrZw1xRh/zt8C9/9e41FGBEREViAkY5yTrq/AkEAetWvAr9qNmJHIqJSYmQLd8x6tw4AYPWJu/hq12Wo1SzCiIioYitUAZaeno7o6Og87VeuXNFbICpbtoU8RFhUAkzlMnzeyVPsOERUygQ1dcWC3vUgkQAbz0bh020Xka1Six2LiIhINAUuwLZt2wYPDw906dIF9erVw9mzZzWPDR48uFjCUemWkJaFeXsjAQCT2teEo4WRyImIqDTq18gFS/v7QiaVYEdoND76IxxZ2SzCiIioYipwATZnzhyEhIQgPDwca9euxYgRI7Bp0yYA4HH9FdT3+64jPjULNR3NMLSZq9hxiKgU6+5bBcvfawC5TIp/Ih7jg99DkKFUiR2LiIioxBW4AFMqlXB0dAQA+Pn54fjx41i5ciW++eYbXuOlArr4IAGbzkUBAL7pXheGMp5OSERvFljHCauDGkJhIMXhyCcYsf480rKyxY5FRERUogr8rtnBwQGXLl3S3LexscGBAwdw7do1rXYq/1RqAV/vvgxBAHrWr4LG7rZiRyKiMqJVTXusH+4PU7kMp249w5BfzyExXSl2LCIiohJT4AJsw4YNcHBw0GqTy+XYvHkzjh07pvdgVHr9cT4Klx4mwlxhgGmdOfAGERVOY3dbbBgZAAsjA1y4/xz9V55GbFKG2LGIiIhKRIELMGdnZzg5OeX7WLNmzfQWiEq3+NQsLNh7HQAwuWNNOJhz4A0iKrwGVa3xx+gmsDdXIDImGb2WB+P20xSxYxERERU7nrhDhTL/v0gkpitRu5IFBjeuJnYcIirDvCpbYMfYpnCzM0V0Qjr6/ByM8AcJYsciIiIqViVWgMXExECt5rDDZVlo1HP8eeEBAGB29zow4MAbRKQjFxsTbPugCeo5W+J5mhIDV53BketPxI5FRERUbHR6B920aVMkJSUVaNpDhw5h/vz5uiyORKRSC/h612UAQB8/ZzR0tRE5ERGVF7ZmCmwe1Rgta9ojXanCqPUXsD3kodixiIiIioVOBdiZM2eQkZH3xOmkpCR8/vnnWm09evTAnj17dFkciWjj2fu48igJFkYGmNqJA28QkX6ZKgzwy5CG6Fm/CrLVAj7ZehGrjt8WOxYREZHeFakA69OnD+bNmweJRIInT/IeKpKamoqFCxdqtV29ehUuLi5FS0miepqcie/35Qy88VlgLdiZKURORETlkdxAikV9fTCqhRsA4Lt/IzH776tQqwWRkxEREemPQVGeVLVqVfz9998QBAE+Pj6wtbWFj48PfHx84Ovri+vXr6NSpUpaz9m5cycGDRqkl9BUsub+ew3JGdmoW8UCgwI48AYRFR+pVIIvu3jBwdwI3/57Db+evIsH8WlYOsAXJvIi/ckiIiIqVYr012zx4sUAcq4DdurUKTx69AhhYWEIDw/Hzp07oVarsWDBAq3nJCUlvXYYeyq9Tt6Mw46waEgkwOzudSGTSsSOREQVwKiW7nCwUOCzbZew/2os+q08jV+GNIKTJS99QUREZZtOHyempqbC0NAQANC9e/c3Tjt48GCsX78eAQEBuiySSlB6lgpf7IwAAAQ1cUX9qtYiJyKiiqS7bxU4Wxtj9G8huBydhO7LTuLXoEaoW8VS7GhERERFptMgHLnFV0EEBARg+fLluiyOStgPh24iKj4NlSyN8GlgLbHjEFEF5FfNBrvGN4OHgxlikzLRd8Vp7L8SI3YsIiKiIuOFnChfVx8lYfWJOwByDj00U/DcCyISh4uNCbaPa4oWHnZIV6ow5vcQrDp+G4LAwTmIiKjsKXQB1qRJE4wbNw6rVq3CuXPn8h2GXp+WLVsGV1dXGBkZISAgAOfOnXvttFeuXEHv3r3h6uoKiUSCpUuXFmu28kqlFjBtxyWo1AI6ezuhvZej2JGIqIKzMDLE2qGNMLhxNQhCzgiJ03ZEQKlSix2NiIioUApdgHXp0gVxcXFYtGgRmjZtCnNzc3h5eWHQoEFYsGAB9u/fn+/Q9EXx559/YvLkyZgxYwZCQ0Ph4+ODwMDA184/LS0N7u7umDdvHgf80MH64Hu4+DAR5kYGmNmtjthxiIgAAAYyKb7pXgczunlBKgH+OP8Ag389i7iUTLGjERERFVihC7CvvvoKW7ZswfXr1xEcHAx7e3vUrVsXEokE69atwzvvvINKlSqhcuXKOodbvHgxRo0ahWHDhsHLywsrVqyAiYkJ1qxZk+/0jRo1wvfff48BAwZAoeC1qooiOiEdC/fnXPNrWqfacLDgiGNEVHpIJBIMa+aGX4IawlQuw5k78ej240mEP0gQOxoREVGB6HRiz9ixY7Fs2TL07NlT0/bvv/9i9OjRCAoK0ilYVlYWQkJCMG3aNE2bVCpF+/btcfr0aZ3m/bLMzExkZv7/p6dJSUkAAKVSCaVSqbfllAWCIODLHZeQlqVCw2pW6O3rVOLbIHd5FW3blybsA/GxD96uRXUbbBsTgPGbw3EnLg19VwRjVrfa6OvnrJf5sw/Exz4QH/tAfOwD8RXHttepALt27Rp8fX212jp37ozly5djyZIluswacXFxUKlUcHTUPv/I0dERkZGROs37ZXPnzsWsWbPytB85cgQmJiZ6W05ZEBYnwdGbMsgkAjpYxWHv3v9Ey3LgwAHRlk052AfiYx+83Rg34HeVFBHPpfhi11X8HXwZvd3UMNDTEFPsA/GxD8THPhAf+0A8aWlpep+nTgVYo0aNsH79esycOVOr3dvb+42DZZQm06ZNw+TJkzX3k5KS4OLigjZt2sDW1lbEZCUrMV2J2f87BSAL41vXwPC21UXJoVQqceDAAXTo0KFQlzkg/WEfiI99UDg91AJWnriLJYduIfiJFKlya/w4wAeVdLhoM/tAfOwD8bEPxMc+EN+zZ8/0Pk+dCrDFixejbdu2uH//Pj7++GPUrVsXWVlZWLRoEezs7HQKZmdnB5lMhtjYWK322NhYvQ6woVAo8j1fzNDQsEK90BftuYa4lCxUtzfF+HYeMDSQiZqnom3/0oh9ID72QcFNbF8LPlVtMHFzGC4+TESvFWfw06AGaOyu2wdp7APxsQ/Exz4QH/tAPMWx3XU6SMPPzw9nz57FgwcP4OvrC2NjY5ibm+PXX3/F3LlzdQoml8vh5+eHQ4cOadrUajUOHTqEJk2a6DRv0nbmzjNsPvcAADCvdz0oRC6+iIiKolVNe+yZ0By1K1kgLiUL7/1yFquO34ZazeuFERFR6aHz1XU9PT1x8OBBREVFITw8HFKpFH5+fqhUqZLO4SZPnoygoCA0bNgQ/v7+WLp0KVJTUzFs2DAAwJAhQ1ClShVNsZeVlYWrV69qbkdHRyM8PBxmZmaoUaOGznnKo5TMbHy27SIAYFBAVTRytRE5ERFR0VW1NcGOsU0xbccl7Ap/hO/+jcSpW8+wqJ8P7Mw4Oi4REYmvSAXY9OnT0b17d/j5+WnaqlatiqpVq+otGAD0798fT58+xfTp0xETEwNfX1/s3btXMzBHVFQUpNL//xLv0aNHqF+/vub+woULsXDhQrRq1QpHjx7Va7by4rt/r+FBfDqqWBljWidPseMQEenMWC7Dkv6+8Hezxaw9V3DsxlN0+uEElvb3RbMauh0eT0REpKsiFWAPHz5Ep06dIJfL0a1bN7z77rto164d5HK5vvNhwoQJmDBhQr6PvVpUubq6QhB4qElBHbvxFJvORgEAFvb1gbkRjy0movJBIpFgUEBV+FWzxoebQ3EjNgXv/3oW41pXx6T2NWEo09MwiURERIVUpL9Aa9asQUxMDDZv3gxzc3NMmjQJdnZ26N27N3777TfEx8frOyfpWWKaElNeHHo4rJkrmlSvOCM+ElHFUcvJHLvHN8dA/6oQBGDZkdvov/I0Hj7X/7DCREREBVHkjwClUilatGiBBQsW4Pr16zh79iwCAgKwcuVKVK5cGS1btsTChQsRHR2tz7ykJzP3XEFsUibc7UwxJZCHHhJR+WUsl2FuL28sG9QA5kYGCI1KQOcfTuC/iMdiRyMiogpIb8dg1K5dG1OmTMGpU6fw4MEDBAUF4cSJE9i8ebO+FkF6svfyY+wMi4ZUAizs5wNjOUc9JKLyr0u9Svh3Ygv4ulghKSMbYzeG4tOtF5GUoRQ7GhERVSDFchC8vb09RowYgd27d+PTTz8tjkVQEcWlZOLLnZcBAB+0qo4GVa1FTkREVHJcbEyw9YMmGNu6OiQSYFvIQwQuOY7jN56KHY2IiCoInQqw5ORkfeWgEiAIAr7cGYFnqVnwdDLHR+09xI5ERFTiDGVSfP6OJ7aOaQJXWxM8TszAkDXn8MXOCKRkZosdj4iIyjmdCrAWLVogJiamQNPGxMRArVbrsjjS0a7waOy7EgtDmQSL+/nygstEVKE1dLXBvx+1wNCmrgCATWej8M7S4zh9+5m4wYiIqFzTqQCrX78+AgICEBkZqdUeHh6Ozp07a7UdOnQI8+fP12VxpIPHiemYvvsKAOCjdh7wqmwhciIiIvGZyA0w89062DQqAFWsjPHweToGrj6Db/6JRJZK7HRERFQe6VSArV27FkOHDkXz5s1x8uRJ3LhxA/369YOfnx9kMu1vV3r06IE9e/boFJaKRhAEfL49AskZ2fBxtsQHraqLHYmIqFRpWt0O+z5uiYH+VQEAG85EYf5FGU7x2zAiItKzIl2I+WWzZs2CQqFAhw4doFKp0K5dO5w+fRr+/v5a0129ehUuLi66Lo6K4NeTd3H8xlMoDKRY1M8XBrwAKRFRHmYKA8zt5Y136jphyraLiE3KxNB1IejuWxlfdfGCvblC7IhERFQO6PROPDY2Fh999BHmzJkDLy8vGBoaYujQoXmKLwDYuXMnBg0apMviqAjCop5j3n85h4h+1aU2ajiYiZyIiKh0a1XTHv992AwtndSQSIDd4Y/QdtFRbDx7H2q1IHY8IiIq43QqwNzc3HD8+HFs3boVISEh2L59O0aPHo3vv/8+z7RJSUlwcnLSZXFUSIlpSkzYFIZstYAu3pXwfuNqYkciIioTzI0M0NtNje1jAlC3igWSM7Lx5c7L6L0iGNceJ4kdj4iIyjCdCrA1a9YgLCwMXbp0AQC88847OHLkCJYsWYLx48drTTt48GCsX79el8VRIQiCgE+3XUR0Qjqq2Zpgbm9vSCQSsWMREZUp3lUssXt8c8zo5gUzhQHCohLQ9ceT+O7fa0jlkPVERFQEOhVgAwYMyNPWoEEDBAcH4/Dhw1rtAQEBWL58uS6Lo0JYc+oeDlyNhVwmxbJBDWBhZCh2JCKiMkkmlWBYMzccnNwKnb2doFILWHX8DtotOobtIQ95WCIRERVKsYzG4OrqiuDg4OKYNRVA+IMEzPvvGgDg6661UbeKpciJiIjKPidLIyx/zw9rhzaCi40xYpIy8MnWi+ix/BTO3Y0XOx4REZURxTYcnrW1dXHNmt4gMU2J8RtDoVTxvC8iouLQxtMBBz5uhc/f8YSZwgCXHiai38rTGLcxBFHP0sSOR0REpRzHIy9HXj7vq6oNz/siIiouRoYyjG1dHUc/a41BAVUhlQD/RsSg/eJjmPvvNSRlKMWOSEREpVShCrD09HRER0fnab9y5YreAlHR8bwvIqKSZWemwHc9vfHvRy3QwsMOWSo1Vh6/gzbfH8XaU3eRoVSJHZGIiEqZAhdg27Ztg4eHB7p06YJ69erh7NmzmscGDx5cLOGo4F4+7+urrrXh7czzvoiISoqnkwV+G+6PtUMbobq9KZ6lZmHWnqtoszDn+mFZ2WqxIxIRUSlR4AJszpw5CAkJQXh4ONauXYsRI0Zg06ZNAHIOfSPxPEnO0Jz31dnbCYN53hcRUYmTSCRo4+mAvZNa4tuedVHJ0giPEzPw5c7LaLvoKLZceIBsFQsxIqKKzqCgEyqVSjg6OgIA/Pz8cPz4cfTs2RO3bt3ieUYiylCqMOq3EEQnpMPNzhTzetdjfxARichQJsV7AdXQu4Ez/jgXhWVHb+Ph83RM2XYJPx+9jY/aeaCbT2XIpPxdTURUERX4GzAHBwdcunRJc9/GxgYHDhzAtWvXtNqp5KjVAj7ZchEXHyTAysQQa4Y24nlfRESlhJGhDEObueH4Z23wRWdP2JjKcTcuFZP+DMc7S49jR+hDKPmNGBFRhVPgAmzDhg1wcHDQapPL5di8eTOOHTum92D0dosOXMc/EY9hKJNgxft+cLMzFTsSERG9wlguw+iW1XF8Sht8FlgLFkYGuPkkBZO3XESrBUfw68m7SM3MFjsmERGVkAIfgujs7AwAaNKkCerXrw9fX1/4+vqiXr16aNasWbEFpPxtvfAAy47cBgDM7VUPjd1tRU5ERERvYqYwwPg2NfB+42r4/cx9rD11D48SMzD776v436GbCGpSDUOausLOTCF2VCIiKkYFLsBydenSBZcuXcKiRYtw+/ZtSCQSeHh4aAqy3J9Xvy0j/Tl9+xm+2BkBAJjQpgb6+DmLnIiIiArK0tgQ49vUwIjmbtgRGo1Vx2/j3rM0/O/wLaw8fgf9GrpgZAs3VLPlUQ1EROVRoQuwr776SnP73Llz6NGjB+rWrQuJRIJ169YhMjISEokEjo6OePTokV7DEnDnaQo++D0ESpWALvUqYXKHmmJHIiKiIjAylGFQQFX0b+SCfVdisOLYbVx6mIgNZ+7j97P30aaWAwY3qYZWHvaQcsAOIqJyo9AF2MvGjh2LZcuWoWfPnpq2f//9F6NHj0ZQUJDO4Ujb89QsDF93HonpSvi6WGFRXx/+USYiKuNkUgk6e1dCp7pOOHMnHiuO3caxG09xOPIJDkc+QVUbE7zfuCr6+rnA2lQudlwiItJRgQfhyM+1a9fg6+ur1da5c2csX74cwcHBusyaXpGZrcKYDSG49ywNVayMsXpIQxgZysSORUREeiKRSNCkui3WD/fH4U9aYURzN1gYGSAqPg3f/RuJxnMP4dOtF3HpYYLYUYmISAc6FWCNGjXC+vXr87R7e3vj3LlzusyaXpKtUuOTLRdx7l48zBUGWDusEezNeZI2EVF55W5vhq+7euHsF+0xv7c36lS2QGa2GttCHuLdn06h648nsPbUXcSnZokdlYiICkmnQxAXL16Mtm3b4v79+/j4449Rt25dZGVlYdGiRbCzs9NXxgotW6XGR3+G459LOcPNL3uvAWo6mosdi4iISoCxXIb+jaqiX0MXhD1IwIbT9/HPpce4HJ2Ey9FX8e0/19DG0wF9/JzRppYD5AY6fa5KREQlQKcCzM/PD2fPnsWECRPg6+sLQ0NDqNVqGBgY4Ndff9VXxgpLqVJj0h/hmmt9LX/PDy1r2osdi4iISphEIkGDqtZoUNUaX3f1wp6Lj7At5CEiohNx4GosDlyNhbWJIbr7VkHvBs6oW8UCEgnPESYiKo10KsAAwNPTEwcPHkRUVBTCw8MhlUrh5+eHSpUq6SNfhaVUqTFxcxj+uxwDuUyKn99vgHa1HcWORUREIrMxlSOoqSuCmrriRmwytoc8xI6waDxNzsS64HtYF3wPbnam6OzthM7eleBVicUYEVFponMBlqtq1aqoWrWqvmZXoWVlq/Hh5lDsuxILuUyKFYMboK0niy8iItJW09Ec0zrXxmeBtXDiVhy2hzzE/quxuBuXimVHbmPZkdtwtTVBJ+9K6OJdCXUqsxgjIhKb3gow0o+sbDUmbArF/quxkBtIsXKwH9rU4kWtiYjo9QxkUrSp5YA2tRyQnKHE4cgn+DfiMY5ef4p7z9Lw89Hb+PnobVS1MUEnbye0r+2I+i5WMJDxnDEiopLGAqwUycpWY9zGUBy8llN8rRrsh9YsvoiIqBDMjXLOBevuWwUpmdk5xdilxzhy/Qmi4tOw8tgdrDx2B1YmhmhV0x5tPR3QqqY9rEx4jTEiopLAAqyUyMxWYfzGUBy89gQKAylWD2nIATeIiEgnZgoDvOtTGe/6VEZqZjaOXH+CfVdicez6EySkKbE7/BF2hz+CVAI0qGqNNp4536J5OplDKuWhikRExYEFWCkQk5iBsRtDEBaVAIWBFL8ENUQLDxZfRESkP6YKA3StVxld61VGtkqN0KgEHI58giORT3A9NhkX7j/HhfvP8f2+67AxlaOJuy2aVLdF0+q2cLMz5bljRER6wgJMZOfuxmPcxlDEpWTCwsgAK973Q9MavIYaEREVHwOZFP5uNvB3s8HUTp6ITkjXFGNn7jxDfGoW/ol4jH8iHgMAKlkavSjG7BDgZgNna2MWZERERcQCTCSCIGBd8D18+881ZKsFeDqZY+VgP1SzNRU7GhERVTBVrIwxuHE1DG5cDUqVGpceJiD41jOcuh2H0PsJeJyYgR2h0dgRGg0AcDBXwK+aNfyqWaNBNWvUqWwBhYFM5LUgIiobWICJID1LhS93RmBHWM4fsnd9KmNeb2+YyNkdREQkLkOZFH7VbOBXzQYftvNAhlKFkPvPEXw7DsG3n+FydCKeJGfiv8sx+O9yDABAbiBFvSqWaFDNGvWcLeFdxRJVbUz4LRkRUT74jr+EPYhPw5gNIbj6OAkyqQTTOnliRHM3/pEiIqJSychQhmY17NDsxeHxGUoVIqITceHec4Tcf47QqOeIT83SnEOWy8LIAHWr5BRjdV78X83GhIN7EFGFxwKsBJ24+RQfbg5DQpoStqZy/DioPppW5/leRERUdhgZytDI1QaNXG0A5BxSf+9ZmqYYuxydiMjHyUjKyEbw7WcIvv1M81wzhQFqOpqhlpMFajmaoaaTOWo5msPWTCHW6hARlTgWYCXgSVIG5u+9ju2hDwEAPs6W+Pl9P1S2MhY5GRERkW4kEgnc7EzhZmeKPn7OAHKua3kjNhmXoxMREZ2Iy4+ScO1xElIysxEalYDQqAStediZyVHT0RweDmZwtzeDm50pqloroBZEWCEiomLGAqwYZShV+PXkXSw/cgupWSoAwKCAqpje1QtGhjxZmYiIyie5gRR1q1iibhVLDHjRplSpcftpCm7EpuBGTDIiY5JxIzYZUfFpiEvJQlyK9rdlAGAokeHnO8Go7pBTlFWzNYGLtQlcbExQydIIBjJpya8cEZGOWIAVA0EQsO9KLL799yoexKcDAHxdrDDz3TrwdbESNxwREZEIDGVSeDpZwNPJAvD5//bUzGzcepKC6zHJuP00BXfiUnE3LhX3n6VCqQJuPEnBjScpeeYnk0pQydLoRUFmDBdrE1S2MkYlSyM4WRqhkqUxjOX8sJOISh8WYHoWGZOEb/Zc1XyK52ihwNROnujuU4UnHhMREb3CVGEAHxcr+LzyAWV6RiY27d6LanUbIep5Ju7EpeBBfDoePE/Dw+fpyMpW4+HzdDx8no7Td/Kft5WJIZwsjFDZyhiOFkZwMFfA/qUfB3MF7MwUPCqFiEoUCzA9UKkFnL79DDtCH2JXeDTUQs7hF2NauuODVtVhquBmJiIiKgwDmRR2RkCrmvYwNDTUekytFvA0JRMP4tPw4HlaTmEWn4aYpAw8SkjH48QMpGWpkJCmREKaEpExyW9cloWRAezMFbAxkcPGVA5bs5z/rU1ybytgZWwIKxNDWBobwtzIEDJ+qEpERVTqK4Nly5bh+++/R0xMDHx8fPDjjz/C39//tdNv3boVX3/9Ne7duwcPDw/Mnz8fnTt31nsuQRAQ/iABu8Mf4e9LjxGXkql5rLO3E6Z1qg0XGxO9L5eIiKiik0olcLQwgqOFERq+GI3xZYIgICkjGzGJGXiUmI6YxAw8TsxAXEomniRl4mlKJuKSM/E0ORNZKjWSMrKRlJGNO0gt0PIlEsBcYQArEzksjXOKMgtjA5grDGFmZABzIwOYGxnC3MgAFkYGMFMYwlQhg6nCIOdHnnPbkOewEVVIpboA+/PPPzF58mSsWLECAQEBWLp0KQIDA3H9+nU4ODjkmT44OBgDBw7E3Llz0bVrV2zatAk9evRAaGgo6tatq3MeQRBw60kK/rr4CLvDHyEqPk3zmJWJIbp4V0JvP2c0qGqt87KIiIioaCQSiaYwquVk/trpBEFAUno2niRnIC4lC/GpWYhPy0J8ShbiUzMRn6ZEfGomnqVkITFdicR0JdKyVBAEaIo2XchlUpgoZDCVG8BYLoOx4YufF7dN5DIYvbitMJDCyFAGI0MpFAba/8sNcm7LDaSQy3Lv5/yf22ao+ZHw2qNEIivVBdjixYsxatQoDBs2DACwYsUK/PPPP1izZg2mTp2aZ/offvgB77zzDj777DMAwOzZs3HgwAH89NNPWLFiRaGW/felx0gW4hGdkHM4Q+5P7miGAGBsKEPHOo7o7lsZzWvYQ27AT7KIiIjKColEAksTQ1iaGMLDsWDPycpWvyjGcoqy3MMcUzKzkZyhRHJGNpIzs3P+z1AiJSPndmpWNtKyVEjJzEZWtjpnXio1stLUSEhTFuNa5mUok2gVZAZSKQxetBlIJZBJc27LpEByggybY87D0EAGmfTFtFIJZDJJzv+SnOllUgmk0pw26cttEglkUkAmySn8Xm6XSpDzv/Sl2xJoppMgp03yoi13mtz7Lz+ec0To/z9fArxolwA5/7TaJXgxn5ynad2XvJjni4de3M5peHka7cf/fx4ve/m+BJI87fk9/nJbtjIbMWnAzScpMDQwyPP4/6fMf5n5edPDbyvO31a661LbS946d3HEP097+0SFVGoLsKysLISEhGDatGmaNqlUivbt2+P06dP5Puf06dOYPHmyVltgYCB27dr12uVkZmYiM/P/Dx9MSkoCAHz91zVIFXkPITSQStDSww7d6jmhrac9TOQvNqGgglKpyjM9FY5SqdT6n0oe+0B87APxsQ/EV1r7QALAykgKKyMjwNqoSPNQqtRIy1JpCrK0LBUylCqkK1VIz3r5fzXSs1TIyFYhQ6lGZrYamUoVMrPVyMhWIVOpRka2Glm5P6qc/zNfup2lUkMQXl2+AKVKBaAg71skuJ38vEjrSfpigLkXg8UOUWGpMytQARYXFweVSgVHR+2PpBwdHREZGZnvc2JiYvKdPiYm5rXLmTt3LmbNmpWnvYaFGvYWalgrAGuFAGv5i/8VgKH0MfDwMY4+LMKKUYEcOHBA7AgVHvtAfOwD8bEPxFcR+8DoxY/WCQ0SAIYvfgpBEAA1AJUaUAk5P9kvbmcL2u0qAVALklfu5/zkd1slAAL+v10tAGpItO7nPq75/0Wel/8XXvlf/WJavPS/kM9zAMn/P+9F+8vT5tadeR57cePlx/HKNHhpupfb8tx+0+P5PPameb3akN91yN96bXLhjXcL89RCK2y2kqTLolVS/QcvtQVYSZk2bZrWt2ZJSUlwcXHBprGtYGtrK2KyikmpVOLAgQPo0KFDnlGvqGSwD8THPhAf+0B87APxsQ/Exz4Q37Nnz1BpsX7nWWoLMDs7O8hkMsTGxmq1x8bGwsnJKd/nODk5FWp6AFAoFFAoFHnaDQ0N+UIXEbe/+NgH4mMfiI99ID72gfjYB+JjH4inOLZ7qR01Qi6Xw8/PD4cOHdK0qdVqHDp0CE2aNMn3OU2aNNGaHsg5dOF10xMREREREZWkUvsNGABMnjwZQUFBaNiwIfz9/bF06VKkpqZqRkUcMmQIqlSpgrlz5wIAPvroI7Rq1QqLFi1Cly5d8Mcff+DChQtYtWqVmKtBREREREQEoJQXYP3798fTp08xffp0xMTEwNfXF3v37tUMtBEVFQWp9P+/xGvatCk2bdqEr776Cl988QU8PDywa9cuvVwDjIiIiIiISFelugADgAkTJmDChAn5Pnb06NE8bX379kXfvn2LORUREREREVHhldpzwIiIiIiIiMobFmBEREREREQlpNQfgljShBdX1EtOTuZwnyJQKpVIS0tDUlISt79I2AfiYx+Ij30gPvaB+NgH4mMfiC85ORnA/9cI+sAC7BXPnj0DALi5uYmchIiIiIiISoNnz57B0tJSL/NiAfYKGxsbADkjLOprI1PBJSUlwcXFBQ8ePICFhYXYcSok9oH42AfiYx+Ij30gPvaB+NgH4ktMTETVqlU1NYI+sAB7Re6w9paWlnyhi8jCwoLbX2TsA/GxD8THPhAf+0B87APxsQ/E9/Klr3Sel97mRERERERERG/EAoyIiIiIiKiEsAB7hUKhwIwZM6BQKMSOUiFx+4uPfSA+9oH42AfiYx+Ij30gPvaB+IqjDySCPsdUJCIiIiIiotfiN2BEREREREQlhAUYERERERFRCWEBRkREREREVEJYgBEREREREZUQFmBEREREREQlpEIWYMuWLYOrqyuMjIwQEBCAc+fOvXH6rVu3wtPTE0ZGRvD29sa///5bQknLp8Js/3Xr1kEikWj9GBkZlWDa8uf48ePo1q0bKleuDIlEgl27dr31OUePHkWDBg2gUChQo0YNrFu3rthzlmeF7YOjR4/m2Q8kEgliYmJKJnA5M3fuXDRq1Ajm5uZwcHBAjx49cP369bc+j38L9KcofcC/B/r1888/o169erCwsICFhQWaNGmC//77743P4T6gX4XtA+4DxWvevHmQSCSYNGnSG6fTx35Q4QqwP//8E5MnT8aMGTMQGhoKHx8fBAYG4smTJ/lOHxwcjIEDB2LEiBEICwtDjx490KNHD1y+fLmEk5cPhd3+AGBhYYHHjx9rfu7fv1+Cicuf1NRU+Pj4YNmyZQWa/u7du+jSpQvatGmD8PBwTJo0CSNHjsS+ffuKOWn5Vdg+yHX9+nWtfcHBwaGYEpZvx44dw/jx43HmzBkcOHAASqUSHTt2RGpq6mufw78F+lWUPgD490CfnJ2dMW/ePISEhODChQto27YtunfvjitXruQ7PfcB/StsHwDcB4rL+fPnsXLlStSrV++N0+ltPxAqGH9/f2H8+PGa+yqVSqhcubIwd+7cfKfv16+f0KVLF622gIAAYcyYMcWas7wq7PZfu3atYGlpWULpKh4Aws6dO984zZQpU4Q6depotfXv318IDAwsxmQVR0H64MiRIwIA4fnz5yWSqaJ58uSJAEA4duzYa6fh34LiVZA+4N+D4mdtbS388ssv+T7GfaBkvKkPuA8Uj+TkZMHDw0M4cOCA0KpVK+Gjjz567bT62g8q1DdgWVlZCAkJQfv27TVtUqkU7du3x+nTp/N9zunTp7WmB4DAwMDXTk+vV5TtDwApKSmoVq0aXFxc3vrJEOkf94HSw9fXF5UqVUKHDh1w6tQpseOUG4mJiQAAGxub107D/aB4FaQPAP49KC4qlQp//PEHUlNT0aRJk3yn4T5QvArSBwD3geIwfvx4dOnSJc/rOz/62g8qVAEWFxcHlUoFR0dHrXZHR8fXnksRExNTqOnp9Yqy/WvVqoU1a9Zg9+7d+P3336FWq9G0aVM8fPiwJCITXr8PJCUlIT09XaRUFUulSpWwYsUKbN++Hdu3b4eLiwtat26N0NBQsaOVeWq1GpMmTUKzZs1Qt27d107HvwXFp6B9wL8H+hcREQEzMzMoFAp88MEH2LlzJ7y8vPKdlvtA8ShMH3Af0L8//vgDoaGhmDt3boGm19d+YFCoqYlKWJMmTbQ+CWratClq166NlStXYvbs2SImIyo5tWrVQq1atTT3mzZtitu3b2PJkiXYsGGDiMnKvvHjx+Py5cs4efKk2FEqrIL2Af8e6F+tWrUQHh6OxMREbNu2DUFBQTh27NhrCwDSv8L0AfcB/Xrw4AE++ugjHDhwoMQHM6lQBZidnR1kMhliY2O12mNjY+Hk5JTvc5ycnAo1Pb1eUbb/qwwNDVG/fn3cunWrOCJSPl63D1hYWMDY2FikVOTv78+iQUcTJkzA33//jePHj8PZ2fmN0/JvQfEoTB+8in8PdCeXy1GjRg0AgJ+fH86fP48ffvgBK1euzDMt94HiUZg+eBX3Ad2EhITgyZMnaNCggaZNpVLh+PHj+Omnn5CZmQmZTKb1HH3tBxXqEES5XA4/Pz8cOnRI06ZWq3Ho0KHXHm/bpEkTrekB4MCBA288PpfyV5Tt/yqVSoWIiAhUqlSpuGLSK7gPlE7h4eHcD4pIEARMmDABO3fuxOHDh+Hm5vbW53A/0K+i9MGr+PdA/9RqNTIzM/N9jPtAyXhTH7yK+4Bu2rVrh4iICISHh2t+GjZsiPfeew/h4eF5ii9Aj/tB4ccKKdv++OMPQaFQCOvWrROuXr0qjB49WrCyshJiYmIEQRCEwYMHC1OnTtVMf+rUKcHAwEBYuHChcO3aNWHGjBmCoaGhEBERIdYqlGmF3f6zZs0S9u3bJ9y+fVsICQkRBgwYIBgZGQlXrlwRaxXKvOTkZCEsLEwICwsTAAiLFy8WwsLChPv37wuCIAhTp04VBg8erJn+zp07gomJifDZZ58J165dE5YtWybIZDJh7969Yq1CmVfYPliyZImwa9cu4ebNm0JERITw0UcfCVKpVDh48KBYq1CmjR07VrC0tBSOHj0qPH78WPOTlpammYZ/C4pXUfqAfw/0a+rUqcKxY8eEu3fvCpcuXRKmTp0qSCQSYf/+/YIgcB8oCYXtA+4Dxe/VURCLaz+ocAWYIAjCjz/+KFStWlWQy+WCv7+/cObMGc1jrVq1EoKCgrSm37Jli1CzZk1BLpcLderUEf75558STly+FGb7T5o0STOto6Oj0LlzZyE0NFSE1OVH7pDmr/7kbvegoCChVatWeZ7j6+sryOVywd3dXVi7dm2J5y5PCtsH8+fPF6pXry4YGRkJNjY2QuvWrYXDhw+LE74cyG/bA9B6XfNvQfEqSh/w74F+DR8+XKhWrZogl8sFe3t7oV27dpo3/oLAfaAkFLYPuA8Uv1cLsOLaDySCIAiF+86MiIiIiIiIiqJCnQNGREREREQkJhZgREREREREJYQFGBERERERUQlhAUZERERERFRCWIARERERERGVEBZgREREREREJYQFGBERERERUQlhAUZERERERFRCWIARERERERGVEBZgRERERdCzZ09YW1ujT58+eR4bOHAgTpw4gdatW8PLywv16tXD1q1bRUhJRESljUQQBEHsEERERGXN0aNHkZycjPXr12Pbtm2adrVajcaNG2P37t2IjY2Fr68vYmJi4Ofnhxs3bsDU1FTE1EREJDZ+A0ZERFQErVu3hrm5eZ72kJAQNGjQAJUqVYKvry8AwMnJCXZ2doiPjy/hlEREVNqwACMiojLn448/Rq9evcSOka/9+/ejY8eOWm0hISFQqVRwcXHJ9zmleX2IiEi/WIAREVGZc+7cOTRs2FDsGPk6dOgQ2rVrp7kfHx+PIUOGYNWqVa99TmleHyIi0i8WYEREVGZkZWXB0NAQwcHB+PLLLyGRSNC4cePXTv/06VNERka+8ScrKyvf50ZFRSEoKAiOjo4wNjaGj48PTp48+cZ8ycnJUKvVsLS0BABkZmaiR48emDp1Kpo2barz+hARUdlnIHYAIiKigjIwMMCpU6cQEBCA8PBwODo6wsjI6LXTL1q0CPPnz3/jPK9duwZPT0+ttvv37yMgIAAtW7bEX3/9BRsbGxw9ehQWFhZvnNfhw4fRpk0bAIAgCBg6dCjatm2LwYMH62V9iIio7OMoiEREVKbs2rULI0eORFxcXLEto3PnzpDL5di1a9drp2nfvj0uXryI1NRU2NjYYOvWrdiwYQMGDx6MJk2a4OTJk2jZsiXq1aunec6GDRvg7e2tNZ+SWB8iIio9+A0YERGVKWFhYfDx8SnQtFOnTi30N2D379/Hf//9h7CwsDc+7+DBg3naPvzwQ/z4448AgObNm0OtVr81Y2HWh4iIyj4WYEREVKaEh4cXuGD55JNPMHTo0DdO4+7unmf+crlcM4R8YVy4cKHQzynM+hARUdnHQTiIiKhMiYiIKFJxVFCGhobIzs5GWlpasS3jZcW9PkREVLrwGzAiIipT1Go1rl+/jkePHsHU1FQz4mB+ijIIR0BAACwtLTF27FhMnToVgiDg+PHjaNeuHTw8PPS2HrkKsz5ERFT28RswIiIqU+bMmYN169ahSpUqmDNnzhunnTdvHgRBeOPPqyMg2traYs+ePbh58yYaNWqE5s2b46+//oKDg4Po60NERGUfR0EkIiIiIiIqIfwGjIiIiIiIqISwACMiIiIiIiohHISDiIjKhDNnzmDDhg0wNzdHr1694OXlhV27dqFt27aoXLmy2PGIiIgKhAUYERGVCevWrUNQUBBiYmKwcOFCREVFoVu3bsU2OAYREVFx4CAcREREREREJYTngBEREREREZUQFmBEREREREQlhAUYERERERFRCWEBRkREREREVEJYgBEREREREZUQFmBEREREREQlhAUYERERERFRCWEBRkREREREVEJYgBEREREREZUQFmBEREREREQl5P8Aa/Lmi2eSergAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Make the plot…\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "def A(T):\n",
    "    return T**2 * np.exp(-T**2)\n",
    "T_values = np.linspace(0, 4, 101)\n",
    "A_values = A(T_values)\n",
    "fig, ax = plt.subplots(figsize=(10, 3))\n",
    "ax.plot(T_values, A_values)\n",
    "ax.set_xlim(0, 4)\n",
    "ax.set_xlabel(r'$\\underset{\\sim}{t} \\equiv c^{1/2}\\,t$')\n",
    "ax.set_ylabel(r'$d^2 \\underset{\\sim}{x} / d\\underset{\\sim}{t}^2$ with $\\underset{\\sim}{x} = C^2\\,x / A$')\n",
    "ax.grid(True)\n",
    "ax.set_title('dimensionless acceleration versus dimensionless time')\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a9c19d3",
   "metadata": {},
   "source": [
    "I've also plotted the function in the accompanying spreadsheet, to illustrate how to do it that way."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e597d17d",
   "metadata": {},
   "source": [
    "### b) Find an expression for the rocket's speed vs. time."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "790acb2b",
   "metadata": {},
   "source": [
    "This is a 1D problem and I can define the rocket's direction of motion to be the positive $x$ direction, so speed vs. time is eqivalent to velocity vs. time.\n",
    "\n",
    "I can get dimensionless velocity vs. time by integrating the above expression for dimensionless acceleration vs. time, noting that the rocket's initial velocity is zero:\n",
    "\n",
    "$$ \\underset{\\sim}{v}(\\underset{\\sim}{t}) = \\cancel{\\underset{\\sim}{v}(0)} + \\int_0^{\\underset{\\sim}{t}} \\underset{\\sim}{a}(\\underset{\\sim}{t})\\, d\\underset{\\sim}{t}' \n",
    "\\ \\Rightarrow \\ \\boxed{ \\underset{\\sim}{v}(\\underset{\\sim}{t}) = \\int_0^{\\underset{\\sim}{t}} \\underset{\\sim}{t}'^2\\, e^{-\\underset{\\sim}{t}'^2}\\,d\\underset{\\sim}{t}' } $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8112423",
   "metadata": {},
   "source": [
    "Unfortunately, this integral is not analytically solvable. This is as far as I can go to \"find an expression for\" the velocity vs. time."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6507767a",
   "metadata": {},
   "source": [
    "### c) Find the rocket's asymptotic $(t \\to \\infty)$ speed."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5a31df9",
   "metadata": {},
   "source": [
    "However, this specific definite integral from 0 to infinity _can_ be evaluated analytically (or looked up or solved with Wolfram Alpha), and the value turns out to be:\n",
    "\n",
    "$$ \\underset{\\sim}{v}(\\infty) = \\int_0^{\\infty} \\underset{\\sim}{t}'^2\\, e^{-\\underset{\\sim}{t}'^2}\\,d\\underset{\\sim}{t}' = \\frac{\\sqrt{\\pi}}{4} $$\n",
    "\n",
    "(If you didn't figure that out, it's okay. You could also proceed with part (d), and then read off an approximate value of the asymptotic speed from the graph.)\n",
    "\n",
    "To get an actual value for the speed (in terms of the unspecified parameters A and C), I'll have to transform back to SI coordinates, using $dx = A\\,dX / c^2$ and $dt = d\\underset{\\sim}{t} / c^{1/2}$:\n",
    "\n",
    "$$ \\begin{split}\n",
    "    v(t) &\\equiv \\frac{dx}{dt} = \\frac{A}{C^2}\\,C^{1/2}\\,\\frac{d\\underset{\\sim}{x}}{d\\underset{\\sim}{t}} = \\frac{A}{C^{3/2}}\\,\\underset{\\sim}{v}(\\underset{\\sim}{t})\n",
    "    \\\\ \\\\\n",
    "    \\Rightarrow \\ v(\\infty) &= \\frac{A}{C^{3/2}}\\,\\underset{\\sim}{v}(\\infty) \\ \\Rightarrow \\ \\boxed{ v(\\infty) = \\frac{A}{C^{3/2}}\\,\\frac{\\sqrt{\\pi}}{4} }\n",
    "\\end{split}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "387acb6f",
   "metadata": {},
   "source": [
    "### d) Make a plot of the rocket's speed vs. time."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c328166",
   "metadata": {},
   "source": [
    "This is just a matter of evaluating the numerical integral for $\\underset{\\sim}{v}(\\underset{\\sim}{t})$ and plotting the result. I'll do that in the accompanying spreadsheet, using the velocity-update portion of the Euler-Cromer update recipe:\n",
    "\n",
    "$$ \\underset{\\sim}{v}(\\underset{\\sim}{t} + \\Delta \\underset{\\sim}{t}) = \\underset{\\sim}{v}(\\underset{\\sim}{t}) + \\underset{\\sim}{a}(\\underset{\\sim}{t})\\,\\Delta \\underset{\\sim}{t} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d00d5dc",
   "metadata": {},
   "source": [
    "### e) Figure out how to determine the time required for the rocket to reach a specified distance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89e8f3aa",
   "metadata": {},
   "source": [
    "Now that I've got numerical values for velocity vs. time, I can numerically integrate that to get position vs. time. I've added a column for that to the spreadsheet, using the position-update portion of the Euler-Cromer recipe:\n",
    "\n",
    "$$ \\underset{\\sim}{x}(\\underset{\\sim}{t} + \\Delta \\underset{\\sim}{t}) = \\underset{\\sim}{x}(\\underset{\\sim}{t}) + \\underset{\\sim}{v}(\\underset{\\sim}{t} + \\Delta \\underset{\\sim}{t})\\,\\Delta \\underset{\\sim}{t} $$\n",
    "\n",
    "To find the time required to reach any specific distance $x_1$, I can (given numerical values for the parameters $A$ and $C$), calculate the corresponding scaled distance $\\underset{\\sim}{x}_1 = C^2\\,x_1 / A$, and then scan down the spreadsheet to find the scaled time $\\underset{\\sim}{t}_1$ at which the rocket's scaled position passes that value. (I can find this to any desired precision by decreasing the step size $\\Delta \\underset{\\sim}{t}$ in the numerical integration.) Then, I can calculate the actual time $t_1 = \\underset{\\sim}{t}_1 / C^{1/2}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b977959",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "628bb487-d42e-4b74-9549-cb1ab030722c",
   "metadata": {},
   "source": [
    "## Problem 2: The Precession of Mercury"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0831c4fd",
   "metadata": {},
   "source": [
    "### a) Spreadsheet Euler-Cromer Solution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d67d22c",
   "metadata": {},
   "source": [
    "I could attack this problem in either Cartesian or plane-polar coordinates. I'll use Cartesian, because Excel doesn't support polar plots (without some really hacky workarounds), so I'll eventually need $\\{x(t), y(t)\\}$ to make a trajectory plot.\n",
    "\n",
    "The only relevant force acting on Mercury is gravity with Einstein's general relativity correctioon, so the net force is\n",
    "\n",
    "$$ \\vec{F} = \\frac{G\\,M\\,m}{r^2} \\left( 1 + \\frac{\\alpha}{r^2} \\right) \\left(-\\hat{r}\\right) $$\n",
    "\n",
    "where $G$ the universal gravitation constant, $M$ is the mass of the Sun, $m$ is the mass of Mercury, and $\\alpha$ is the general relativity parameter.\n",
    "\n",
    "I want to rewrite this in Cartesian coordinates. To do that, I'll note that\n",
    "\n",
    "$$ \\hat{r} = \\frac{\\vec{r}}{r} \\quad \\Rightarrow \\quad \\vec{F} = \\frac{G\\,M\\,m}{r^3} \\left( 1 + \\frac{\\alpha}{r^2} \\right) \\left(-\\vec{r}\\right) $$\n",
    "\n",
    "and $\\vec{r} = x\\,\\hat{\\imath} + y\\,\\hat{\\jmath}$, so that\n",
    "\n",
    "$$ \\begin{align*}\n",
    "    F_x &= -\\frac{G\\,M\\,m}{r^3} \\left( 1 + \\frac{\\alpha}{r^2} \\right) \\, x\n",
    "    \\\\ \\\\\n",
    "    F_y &= -\\frac{G\\,M\\,m}{r^3} \\left( 1 + \\frac{\\alpha}{r^2} \\right) \\, y\n",
    "\\end{align*} $$\n",
    "\n",
    "with $r = \\sqrt{x^2 + y^2}$. Newton's second law therefore yields:\n",
    "\n",
    "$$ \\begin{align*}\n",
    "    a_x &= -\\frac{G\\,M}{r^3} \\left( 1 + \\frac{\\alpha}{r^2} \\right) \\, x\n",
    "    \\\\ \\\\\n",
    "    a_y &= -\\frac{G\\,M}{r^3} \\left( 1 + \\frac{\\alpha}{r^2} \\right) \\, y\n",
    "\\end{align*} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05b64f45",
   "metadata": {},
   "source": [
    "Breaking these two second-order ODEs into four first-order ODEs, I get:\n",
    "\n",
    "$$ \\begin{align*}\n",
    "    \\frac{dv_x}{dt} &= -\\frac{G\\,M}{r^3} \\left( 1 + \\frac{\\alpha}{r^2} \\right) \\, x\n",
    "    \\qquad\n",
    "    \\frac{dx}{dt} &= v_x\n",
    "    \\\\ \\\\\n",
    "    \\frac{dv_y}{dt} &= -\\frac{G\\,M}{r^3} \\left( 1 + \\frac{\\alpha}{r^2} \\right) \\, y\n",
    "    \\qquad\n",
    "    \\frac{dy}{dt} &= v_y\n",
    "\\end{align*} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a2f287c",
   "metadata": {},
   "source": [
    "The usual Euler-Cromer update procedure is therefore:\n",
    "\n",
    "$$ \\begin{align*}\n",
    "    v_x(t + \\Delta t) &= v_x(t) - \\frac{G\\,M}{r(t)^3} \\left( 1 + \\frac{\\alpha}{r(t)^2} \\right) \\, x(t)\\, \\Delta t\n",
    "    \\quad & \\quad\n",
    "    x(t + \\Delta t) &= x(t) + v_x(t + \\Delta t)\\, \\Delta t\n",
    "    \\\\ \\\\\n",
    "    v_y(t + \\Delta t) &= v_y(t) - \\frac{G\\,M}{r(t)^3} \\left( 1 + \\frac{\\alpha}{r(t)^2} \\right) \\, y(t)\\, \\Delta t\n",
    "    \\quad & \\quad\n",
    "    y(t + \\Delta t) &= y(t) + v_y(t + \\Delta t)\\, \\Delta t\n",
    "\\end{align*} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8134b0fa",
   "metadata": {},
   "source": [
    "where $r(t)$ is just a shorthand for $\\sqrt{x(t)^2 + y(t)^2}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "886051f9",
   "metadata": {},
   "source": [
    "The second sheet (tab) of the accompanying spreadsheet file, `PS-N1-Solutions-2.xlsx`, has an implementation of these equations. To simplify the formulas, I've included columns for the acceleration components and for $r$. I also included a column for $\\theta$ just to make it easy to see when each orbit is completed.\n",
    "\n",
    "(Note: I used an Excel trick of giving a name such as `dT` to specific cells, so that I can then call them by that name in formulas rather thn using the hard-to-read cell references like `$L$2`. It's helpful, but not essential.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f591cc25",
   "metadata": {},
   "source": [
    "### b) Numerical Investigation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83a4fe62",
   "metadata": {},
   "source": [
    "By playing around with different values of $\\alpha$, I discovered that:\n",
    "\n",
    "- For $\\alpha = 0$, the orbit is an apparently-perfect ellipse, as expected. The eccentricity also looks credible, since according to NASA, Mercury's orbit has an eccentricity of about 0.2, which implies a difference between perihelion and aphelion of about 41%.\n",
    "\n",
    "- For small values of $\\alpha$, say around 0.003, the orbit still looks elliptical but precesses slightly with each revolution.\n",
    "\n",
    "- Increasing $\\alpha$ increases the precession rate.\n",
    "\n",
    "- When $\\alpha$ gets to around 0.03, the orbit ceases to look very elliptical at all, but the details depend strongly on the time-step chosen. For a step size of $\\Delta t = 0.001$ yr, the planet escaped the sun entirely whenever $\\alpha \\gtrsim 0.033$. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b663224c-398f-4e2a-b81c-47bb124135fe",
   "metadata": {},
   "source": [
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}