CE_UQ_Propagation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5db64b5b",
   "metadata": {},
   "source": [
    "# Uncertainty propagation\n",
    "\n",
    "Using the force on the sail example seen in class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "66f46415",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n",
      "Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4922d8ce",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp8AAAEnCAYAAADvmSiNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABLGElEQVR4nO3de1xUdf4/8NeZOyAzyh0EEc1bXjDBG2ZqJq662nW1db+plZWlldrVdTcv22+prVy7qbWZbmVm28WyzGJXvGspoal4FwQVREABuQzMzOf3xzCjI6AMzsyBmdfz0XkwnDmX99Hx3Xs+n/P5HEkIIUBERERE5AEKuQMgIiIiIt/B4pOIiIiIPIbFJxERERF5DItPIiIiIvIYFp9ERERE5DEsPomIiIjIY1h8EhEREZHHsPgkIiIiIo9h8UlEREREHsPiswX54osvIEkS1qxZU+e9+Ph4SJKEH3/8sc57HTt2RJ8+fQAAmzZtgiRJ2LRpk8viys7OhiRJWLlypcuO2RzOW9/xd+zYgfnz5+PixYtuOee6deswduxYhIeHQ6PRICgoCMOHD8eqVatQU1Nj306SJPuiVCrRpk0bxMfH47HHHsOuXbsavJb6lsTERLdcC1FzxVzq2fN6MpdOmTLFIb9ptVp06dIF8+bNQ1VVlX27+fPnO2zn7++P6OhojBw5Em+//TbKysque+wrl++++86l1+HtWHy2IEOHDoUkSUhLS3NYX1xcjP379yMgIKDOe6dPn8bJkycxbNgwAECfPn2wc+dOewKlhkVGRmLnzp0YM2aMfd2OHTuwYMEClydMIQQefPBBjBs3DhaLBYsWLcJ///tf/Pvf/0Z8fDyeeOIJLFmyxGGf++67Dzt37sS2bdvw2WefYdKkSdi1axcGDhyIp59+ut7zPPnkk9i5c6fD4un/0RHJjbnUszyZSwHAz8/Pnt/Wrl2L/v37Y+HChZg8eXKdbTds2ICdO3diw4YNeP3119GuXTs8//zz6N69O/bt23fNY1+53HrrrS6/Dm+mkjsAaryQkBD06NGjzjftzZs3Q6VS4eGHH66TMG2/2xKmXq/HgAEDPBJvS6fVaj32Z/Xaa69h5cqVWLBgAV566SWH98aOHYvnn38ex48fd1gfHh7uEN/IkSMxc+ZMPProo3jrrbfQtWtXPP744w77tGvXjn//5POYSz3Lk7kUABQKhcP5Ro0ahezsbHz++edYtGgR2rZta38vISEBISEh9t/vv/9+zJgxA0OGDMG4ceNw9OhRaLXaBo9NTcOWzxZm2LBhOHLkCPLy8uzrNm3ahL59+2L06NFIT0936C7YtGkTlEolBg8ebP/96q6iKVOmoFWrVjh+/DhGjx6NVq1aISYmBs888wyMRqPD+c+ePYvx48cjMDAQBoMBEyZMQH5+/nXj3rdvHyRJwvLly+u898MPP0CSJHz77bf2dceOHcPEiRMRFhYGrVaLbt264d13323Un9G2bdswfPhwBAYGwt/fH0lJSfj+++/rbHfmzBk8+uijiImJgUajQVRUFO677z6cO3cOQN2uovnz5+O5554DAMTFxdm7WzZt2oSHH34YQUFBqKioqHOe22+/Hd27d28w3pqaGrz66qvo2rUr/vrXv9a7TURERKO+WSuVSrzzzjsICQnBa6+9dt3tiXwVc+n1tbRcei22gvHUqVPX3TY+Ph5z585FTk5Ovbdm0I1j8dnC2L51X5nw0tLSMGTIEAwaNAiSJGHr1q0O7/Xp0wcGg+Gax62pqcG4ceMwfPhwfPPNN3jooYfwz3/+E6+++qp9m8rKStxxxx346aefkJKSgv/85z+IiIjAhAkTrht3fHw8brnlFqxYsaLOeytXrkRYWBhGjx4NAMjMzETfvn1x4MABvPHGG/juu+8wZswYPPXUU1iwYME1z7N582bcfvvtKCkpwfLly7F69WoEBgZi7NixDknkzJkz6Nu3L77++mvMnj0bP/zwAxYvXgyDwYALFy7Ue+ypU6fiySefBAB89dVX9u6WPn364Omnn8aFCxfw6aefOuyTmZmJtLQ0TJ8+vcGY9+zZg+LiYtx5552QJOma19cYfn5+uOOOO5CVlYXTp087vGexWGAymRwWIcQNn5OopWEu9b5cei22nqPQ0NBGbT9u3DgAwJYtW+q8d3UONZvNTYrJpwlqUYqLi4VCoRCPPvqoEEKIwsJCIUmS2LBhgxBCiH79+olnn31WCCFETk6OACCef/55+/5paWkCgEhLS7Ovmzx5sgAgPv/8c4dzjR49WnTp0sX++9KlSwUA8c033zhs98gjjwgAYsWKFdeM/a233hIAxJEjRxyuR6vVimeeeca+buTIkSI6OlqUlJQ47D9jxgyh0+lEcXGxEEKIrKysOucdMGCACAsLE2VlZfZ1JpNJ9OjRQ0RHRwuLxSKEEOKhhx4SarVaZGZmNhhvfcd/7bXXBACRlZVVZ/shQ4aI3r17O6x7/PHHhV6vd4jnap999pkAIJYtW9bgNlcDIKZPn97g+y+88IIAIH7++WeHa6lvSU1NbfR5ibwFc6n35VIhrH8HAQEBoqamRtTU1Ijz58+LN998U0iSJPr27Wvfbt68eQKAOH/+fL3HqaysFADEqFGjHI5dXw4dNGjQNWOiutjy2cLYRjXbvq1v3rwZSqUSgwYNAgAMGTLEfm/S1fcoXYskSRg7dqzDul69ejl0UaSlpSEwMND+jdBm4sSJjYr9T3/6E7RarcMAl9WrV8NoNOLBBx8EAFRVVeF///sf7r77bvj7+zt8uxw9ejSqqqrqHc0NAOXl5fj5559x3333oVWrVvb1SqUSDzzwAE6fPo0jR44AsHZPDRs2DN26dWtU7I3x9NNPY+/evdi+fTsAoLS0FB9//DEmT57sEI8niAZaM59++mns3r3bYenfv79HYyNqDphLvTeXlpeXQ61WQ61WIzQ0FDNnzsSoUaPw9ddfNzqGhnKon59fnRxa3y0QdG0sPlugYcOG4ejRozh79izS0tKQkJBg/wc5ZMgQZGRkoKSkBGlpaVCpVI26V9Df3x86nc5hnVardZiaoqioCOHh4XX2jYiIaFTcQUFBGDduHD766CN7N8XKlSvRr18/+308RUVFMJlMePvtt+3Jw7bYupIKCwvrPf6FCxcghEBkZGSd96KiouzHB4Dz588jOjq6UXE31p133on27dvb76dauXIlysvLr9tN1K5dOwBAVlaWy2Kx/Y/Odt020dHRSExMdFgCAwNddl6iloS51Ltyqc2VBeJvv/2Gixcv4vvvv3cYaHQ9DeVQhUJRJ4d26dKl0cclKxafLdCV9ypt2rQJQ4YMsb9nS45btmyx3zzvqla34OBg+w3kV2rMTfI2Dz74IM6cOYPU1FRkZmZi9+7d9m/qgLU1QqlUYsqUKXW+XdoWW+K8Wps2baBQKBwGENicPXsWAOyjGkNDQ+vcD3mjFAoFpk+fji+++AJ5eXlYsmQJhg8fft3ElJiYiKCgIHzzzTcuuf+ysrIS//3vf9GxY0eX/0+ByJswl3pXLr1yf1th2LNnT+j1eqdjsA3aGjp0qNP70vWx+GyBbrvtNiiVSnzxxRc4ePCgwz8Og8GA3r1749///jeys7Mb1U3UWMOGDUNZWZnDSEoAdW4Mv5bk5GS0bdsWK1aswIoVK6DT6fDHP/7R/r6/vz+GDRuGjIwM9OrVq843zMTERAQHB9d77ICAAPTv3x9fffUVKisr7estFgs++eQTREdHo3PnzgCsU2+kpaXZu44ayzblxpXHv9LUqVOh0Wjwpz/9CUeOHMGMGTOue0y1Wo0XXngBhw8fxt/+9rd6tykoKLB3QV2L2WzGjBkzUFRUhBdeeOG62xP5MuZS78qlrrJv3z78/e9/R/v27TF+/HiPndeXcJ7PFkiv16NPnz5Yu3YtFAqF/R4lmyFDhmDx4sUAGnePUmNNmjQJ//znPzFp0iT8v//3/9CpUyesX7++3ieBNESpVGLSpElYtGgR9Ho97rnnnjqjR998803ceuutGDx4MB5//HG0b98eZWVlOH78ONatW4eNGzc2ePyUlBSMGDECw4YNw7PPPguNRoMlS5bgwIEDWL16tX00+cKFC/HDDz/gtttuw5///Gf07NkTFy9exIYNGzB79mx07dq13uP37NnTHuPkyZOhVqvRpUsXe9d169atMWnSJCxduhSxsbF17v1qyHPPPYdDhw5h3rx5+OWXXzBx4kTExMSgpKQEW7Zswfvvv48FCxY4/F2fO3cOu3btghACZWVlOHDgAD766CPs27cPs2bNwiOPPNKocxP5KuZS78ulzkpPT4fBYEBNTQ3Onj2L//3vf/j4448RFhaGdevWQaPRuOW8Pk/O0U7UdM8//7wAIBITE+u8t3btWgFAaDQaUV5e7vBeQyM0AwIC6hzHNhrwSqdPnxb33nuvaNWqlQgMDBT33nuv2LFjR6NGaNocPXr0uiOts7KyxEMPPSTatm0r1Gq1CA0NFUlJSeLll1922Ka+827dulXcfvvtIiAgQPj5+YkBAwaIdevW1TlHbm6ueOihh0RERIRQq9UiKipKjB8/Xpw7d+6ax58zZ46IiooSCoWizp+lEEJs2rRJABCvvPJKo/48rvTNN9+IMWPGiNDQUKFSqUSbNm3EsGHDxLJly4TRaLRvhytGWioUCqHX60XPnj3Fo48+Knbu3FnnuLZree2115yOicibMZd6Vy5t6O/gara/E9ui1WpFZGSkSE5OFm+++aYoLS1t8rHp+iQhOMkfkSs988wzWLp0KXJzcxvs1iIiomtjLvVe7HYncpFdu3bh6NGjWLJkCR577DEmSyKiJmAu9X5s+SRyEUmS4O/vj9GjR2PFihUen9uTiMgbMJd6PxafREREROQxnGqJiIiIiDyGxScREREReQyLTyIiIiLymBYx2t1iseDs2bMIDAy0T2xLRORKonay/qioKCgU3ve9nHmUiNytsXm0RRSfZ8+eRUxMjNxhEJEPyM3NRXR0tNxhuBzzKBF5yvXyaIsoPm2P28rNzYVer5c5GiLyRqWlpYiJibHnG2/DPEpE7tbYPNoiik9bF5Fer2fSJCK38tYuaeZRIvKU6+VR77uxiYiIiIiaLRafREREROQxLD6JiIiIyGOcLj63bNmCsWPHIioqCpIkYe3atdfdZ/PmzUhISIBOp0OHDh2wbNmypsRKROQ1lixZgri4OOh0OiQkJGDr1q0Nbrtp0yZIklRnOXz4sAcjJiJyDaeLz/LycsTHx+Odd95p1PZZWVkYPXo0Bg8ejIyMDPz5z3/GU089hS+//NLpYImIvMGaNWswc+ZMzJ07FxkZGRg8eDBGjRqFnJyca+535MgR5OXl2ZdOnTp5KGIiItdxerT7qFGjMGrUqEZvv2zZMrRr1w6LFy8GAHTr1g179uzB66+/jnvvvdfZ0xMRtXiLFi3Cww8/jKlTpwIAFi9ejB9//BFLly5FSkpKg/uFhYWhdevWHoqSiMg93D7V0s6dO5GcnOywbuTIkVi+fDlqamqgVqvdHQJRi2CxCJiFgEUIWCyAWQiYLQKi9qdZCAgB64La17A+UUII6zHqfQ+X34PDe1dse8VrZ49zJXH1itp9HX5vzD5Xb1XvNtc+jkalQL+4oLo7yqy6uhrp6el48cUXHdYnJydjx44d19z3lltuQVVVFW6++Wb85S9/wbBhw9wWZ7nRhEN5pTCaLBh0U4jbzkNEvsftxWd+fj7Cw8Md1oWHh8NkMqGwsBCRkZF19jEajTAajfbfS0tL3R0meTEhBCqqzSg3mlBZY0ZljRlVNRZUVptRZTKjqvZnZbWl9j3rUm2yoMYsUGO2wGSxoNpU/+sak0CNxYIa8+XXJrO1YLSIyz8tAtbXtYWktbC8XGSSa4Xrtfj5z3fIHUYdhYWFMJvN9ebF/Pz8eveJjIzE+++/j4SEBBiNRnz88ccYPnw4Nm3ahNtuu63efW40j+47fRET//Uz4kICkPbsUKf2JSK6Fo9MMn/1ZKO2lpCGJiFNSUnBggUL3B4XtRxCCJRWmlBcUY3iciOKy2tQXG5EUXk1SiprcKnKhLIqE8qqanDJaHt9+XdvqO0kCVBIEqTa1xIk1P5n/12y/27dDlf8Dvt2l9+XajeSrnEc+362c9ZznLrBNmpVvTmg/u2u3qae/a5aFRSgqS+yZqO+vNhQTuzSpQu6dOli/33gwIHIzc3F66+/3mDxeaN5NKaNPwDgzIVKWCwCCoV3Tr5PRJ7n9uIzIiKizrf5goICqFQqBAcH17vPnDlzMHv2bPvvtsc1kXeqrDYjr6QS+SVVyCupQn5plcPvBWVGXCivhukGK0iFBPipldDZFwX8NEroVEr4aZTQ1v7UqRTQqZXQqhRQqxRQKySoldbXKoUEjUoBtdLxtVqpgEopQVP7WqmQrIskQaGA/bUk1V2vkKyLfRsFoKz9XWH/6b1P3vE1ISEhUCqV9ebFq1tDr2XAgAH45JNPGnz/RvNopEEHpUJCtdmCc2VViDT4NXpfIqJrcXvxOXDgQKxbt85h3U8//YTExMQG7/fUarXQarXuDo08qKrGjKzCcmQXliOryPozu6gC2YXlKCgzXv8AtVppVWgToEZQgBZB/tafrf3VCNSp0Eqrgl6nRiudCoE6FQJ16tp1KrTSqeCnVrKAI9lpNBokJCQgNTUVd999t319amoq7rzzzkYfJyMjo97blmxuNI+qlApEtdYht7gSucWVLD6JyGWcLj4vXbqE48eP23/PysrC3r17ERQUhHbt2mHOnDk4c+YMPvroIwDAtGnT8M4772D27Nl45JFHsHPnTixfvhyrV6923VVQs5JXUonMs6U4nF+GzLxSHM4rRVZh+TW7vgM0SkQYdIg0+NX+1Nl/hgXqENxKgzb+GujUSs9dCJGbzJ49Gw888AASExMxcOBAvP/++8jJycG0adMAoE4eXbx4Mdq3b4/u3bujuroan3zyCb788ku3T1kX08YfucWVOH2holkO3iKilsnp4nPPnj0OIyxt3TqTJ0/GypUrkZeX5zBXXVxcHNavX49Zs2bh3XffRVRUFN566y1Os+QlqmrMOHi2BL+euoiM3Av49dRF5JdW1butwU+N9iEBiAv2t/4MCUD7YOti8OesB+Q7JkyYgKKiIixcuBB5eXno0aMH1q9fj9jYWACok0erq6vx7LPP4syZM/Dz80P37t3x/fffY/To0W6N03rfZxFyiyvdeh4i8i2SuHoelGaotLQUBoMBJSUl0Ov1cofj0ywWgYNnS7H1+HlsO1aIPdkXUG22OGyjVEjoFNYK3SL16BoRiK6RenSLCERooJbd3tRseXueacr1vbPxGF7/6SjuS4jG63+Id3OERNTSNTbPeGS0O7VsZVU12Hi4AD9lnsOO44W4UFHj8H5IKy36tGuNW9q1wS3tWqNXtAH+Gn60iFq6mCDriPfc4gqZIyEib8IKgepVVlWDHw7k44f9edh+vMihdbOVVoUBHYIxuFMIbu0Ugg4hAWzRJPJC0W2sg4xOX2C3OxG5DotPsjNbBHacKMQX6afx48F8VNVcLjg7hAbgd90jcHvXMMTHtIZaqZAxUiLyBNtcn3kllagxW/jvnohcgsUn4WJFNVb/kotPdp3CmYuXWzg6hgbgzt5tMapHBDqFB8oYIRHJITRQC61KAaPJgryLVWgX7C93SETkBVh8+rDjBZfw4fYsfPXraXsrp8FPjXHxUbg3IRrx0QZ2pxP5MEmSEN3GDyfOlyP3QgWLTyJyCRafPuh4wSW89b9jWPfbWdjmOugWqceDg9pjXHwU59IkIruYIH9r8clBR0TkIiw+fUh2YTkW//covt131j7h+x3dwjF1cBz6xwWxlZOI6rDd95l7gcUnEbkGi08fcMlowtsbj+HDbVmoMVurzhE3h2PmHZ3QPcogc3RE1JzFBFlHvHOieSJyFRafXkwIgbV7zyBl/WH789Nv6xyK55K7oGc0i04iur5otnwSkYux+PRSucUVeP6L37DzZBEAIDbYHy/9/mbc3jWM3etE1Gi2bnfO9UlErsLi08sIIfD5nlz87btDuGQ0wU+txIzbb8LUwXHQqjiQiIicY+t2P19mRFWNmQMSieiGsfj0IufLjHjhy9+w8XABACAxtg1e/0M82ocEyBwZEbVUBj81ArUqlBlNOH2hAjeFcc5fIroxLD69xO7sYkxf9SsKyozQKBV4Jrkzpg7uAKWCXexE1HSSJCE6yB+H8kqRW1zJ4pOIbhiLzxZOCIEPt2fj7+sPwWwR6BTWCu9M7IMuEfwfBBG5RkwbP2vxyUFHROQCLD5bsBqzBX/5+gDW7MkFANzZOwp/v7snArT8ayUi14mtfbLRqSIWn0R041iltFAlFTV4fFU6dpwogkIC/vr7mzElqT1HshORy8UGW+8bP1VULnMkROQNWHy2QAVlVZi0/Bcczi9DgEaJdyb2wbCuYXKHRUReqn1t8ZlVyOKTiG4ci88W5vSFCvzfBz8ju6gCYYFarHywH26O0ssdFhF5MVu3e25xJcwWwYGMRHRDFHIHQI2XU1SB8ct2IruoAtFt/PDFtCQWnkTkdlGt/aBRKlBttiCvhJPNE9GNYfHZQpy9WIk//msXzpZUoWNoAL6YloR2ta0RRETupFRI9snmsws56IiIbgyLzxagoKwKf/rgZ5y5WIm4kACsfnQAIgw6ucMiIh9iu+8zm4OOiOgGsfhs5sqqajBp+S/IKixH29Z+WDW1P8ICWXgSkWdxxDsRuQqLz2asxmzBE6t+xeH8MoQGavHpI/0R1dpP7rCIyAfFhVhv88litzsR3SAWn82UEAJ/XXsAW48Vwk+txIeT+9pbHoiIPI0tn0TkKiw+m6nl27Lw2e5cKCTg7T/egp7RBrlDIiIfZrvn81RxBSwWIXM0RNSSsfhshnaeKELKD4cBAH8ZczPuuDlc5oiIyNdFtdZBpZBQbbIgr7RK7nCIqAVj8dnM5JdU4cnVv8JsEbj7lrZ4cFB7uUMiIoJKqUC7oNpnvPNJR0R0A1h8NiM1Zgumf/orCi9Vo2tEIP5+d08+q52Img3bk46yizjoiIiajsVnM/L2/44h/dQFBOpUeO+BBPhplHKHRERkF8u5PonIBVh8NhN7sovxTtpxAMDf7+7Jke1E1Oy0t7V8studiG4Ai89moKyqBjPX7IVFAPfc0hZj46PkDomIqI72IbbpltjtTkRN16Tic8mSJYiLi4NOp0NCQgK2bt16ze1XrVqF+Ph4+Pv7IzIyEg8++CCKioqaFLA3evm7Qzh9oRIxQX5YcGd3ucMhIg9wNo/abN++HSqVCr1793ZvgPW48hGbnG6JiJrK6eJzzZo1mDlzJubOnYuMjAwMHjwYo0aNQk5OTr3bb9u2DZMmTcLDDz+MgwcP4j//+Q92796NqVOn3nDw3mDH8UKs2ZMLAHjjD70RqFPLHBERuZuzedSmpKQEkyZNwvDhwz0UqaPoNn5QKyUYOd0SEd0Ap4vPRYsW4eGHH8bUqVPRrVs3LF68GDExMVi6dGm92+/atQvt27fHU089hbi4ONx666147LHHsGfPnhsOvqWrqjFjztf7AQAPDIhFv7ggmSMiIk9wNo/aPPbYY5g4cSIGDhzooUgdqZQK+/3oJwouyRIDEbV8ThWf1dXVSE9PR3JyssP65ORk7Nixo959kpKScPr0aaxfvx5CCJw7dw5ffPEFxowZ0/SovcTi/x7DqaIKROh1eP53XeQOh4g8oCl5FABWrFiBEydOYN68ee4O8Zo6htYWn+dZfBJR06ic2biwsBBmsxnh4Y5P3AkPD0d+fn69+yQlJWHVqlWYMGECqqqqYDKZMG7cOLz99tsNnsdoNMJoNNp/Ly0tdSbMFuFwfin+tfUkAOBvd/VgdzuRj2hKHj127BhefPFFbN26FSpV49K2u/Jox9BWAM6x+CSiJmvSgKOrJz4XQjQ4GXpmZiaeeuopvPTSS0hPT8eGDRuQlZWFadOmNXj8lJQUGAwG+xITE9OUMJstIQQWfJsJs0VgZPdwjODjM4l8TmPzqNlsxsSJE7FgwQJ07ty50cd3Vx61Fp/AcXa7E1ETOVV8hoSEQKlU1vl2XlBQUOdbvE1KSgoGDRqE5557Dr169cLIkSOxZMkSfPjhh8jLy6t3nzlz5qCkpMS+5ObmOhNms/fjwXPYebIIGpUCfxlzs9zhEJEHOZtHy8rKsGfPHsyYMQMqlQoqlQoLFy7Evn37oFKpsHHjxnrP46482jHMWnyeOM+5PomoaZwqPjUaDRISEpCamuqwPjU1FUlJSfXuU1FRAYXC8TRKpfXJPULUP1WHVquFXq93WLxFVY0Z/299JgDgsds6IKb2WclE5BuczaN6vR779+/H3r177cu0adPQpUsX7N27F/3796/3PO7Ko7Z7Ps+XGVFSWeOSYxKRb3Hqnk8AmD17Nh544AEkJiZi4MCBeP/995GTk2PvRp8zZw7OnDmDjz76CAAwduxYPPLII1i6dClGjhyJvLw8zJw5E/369UNUlO9Npr58WxZyiysRodfh8aEd5Q6HiGTgTB5VKBTo0aOHw/5hYWHQ6XR11ntCoE6NcL0W50qNOHn+Em5p18bjMRBRy+Z08TlhwgQUFRVh4cKFyMvLQ48ePbB+/XrExsYCAPLy8hzmqpsyZQrKysrwzjvv4JlnnkHr1q1x++2349VXX3XdVbQQF8qrsWzTCQDAi6O6wl/j9B8/EXkBZ/Noc9MxtBXOlRpx4nw5i08icpokGur7bkZKS0thMBhQUlLSorvgU344hPc2n8TNkXp89+StUCjqH6RFRJ7nLXmmIa68vr+uPYCPd53C40M74oXfdXVRhETU0jU2z/DZ7h5yrrQKK7dnAwCeG9mFhScRtVj2uT454p2ImoDFp4e8vfEYjCYLEmPbYGiXULnDISJqsssj3ll8EpHzWHx6QG5xBT77xTrNyXMjuzQ4JyoRUUtgm+vzVFEFaswWmaMhopaGxacHLN18AiaLwOBOIejfIVjucIiIbkikQQd/jRImi0BOcYXc4RBRC8Pi080KSqvwxZ7TAIAnb+8kczRERDdOkiQ+6YiImozFp5t9sC0L1WbrvZ794oLkDoeIyCXsg4543ycROYnFpxuVVNRg1a5TAIAnhnFCeSLyHvaWz3MsPonIOSw+3ejfO7NRXm1G14hADOsSJnc4REQu0zkiEABwtKBM5kiIqKVh8ekmVTVmrNyRDQB4fGhHjnAnIq/SJdxafB47dwlmS7N/VgkRNSMsPt1k3b6zKC6vRpRBhzE9I+UOh4jIpdoF+UOnVsBosuBUUbnc4RBRC8Li0w2EEPZWzwcGtodKyT9mIvIuCoWEzrWtn0fy2fVORI3HqsgN9py6gINnS6FVKXB/3xi5wyEicgt78XmOxScRNR6LTzewtXrefUtbtAnQyBsMEZGbdLUNOmLxSUROYPHpYnklldhwIB8AMDmpvbzBEBG5ka3l8zC73YnICSw+XWz1zzkwWwT6xwWhW6Re7nCIiNzG1vKZXViOqhqzzNEQUUvB4tOFzBaB/6RbH6X5fwNiZY6GiMi9QgO1aO2vhkXwMZtE1HgsPl1oy7HzyCupQmt/NZK7h8sdDhGRW0nS5RHvvO+TiBqLxacLfb47FwBwV++20KqUMkdDROR+tq53TrdERI3F4tNFii4Z8d9D5wAAEzi9EhH5CE63RETOYvHpIl9nnEGNWSA+2sCBRkTkM7rYpltiyycRNRKLTxcQQmBNbZf7eLZ6EpEPsbV8ni2pQklljczREFFLwOLTBQ6eLcWxgkvQqRUYGx8ldzhERB5j8FMjyqADwPs+iahxWHy6wDd7zwAAhncLh16nljkaIiLPujnKeqvRwbMlMkdCRC0Bi88bZLYIfLvvLADgTrZ6EpEPujnKAMDaC0REdD0sPm/QL1nFOFdqhF6nwpAuoXKHQ0Tkcd3tLZ8sPono+lh83iBbq+eoHpGc25OIfJKt+Dx2rgxGEx+zSUTXxuLzBlSbLFi/Pw8AcGdvdrkTkW9q29oPBj81TBaBo/l8zCYRXRuLzxuw5eh5lFTWICxQi/4dguUOh4hIFpIkXdH1zkFHRHRtLD5vwHe/Wbvcx/SKhFIhyRwNEZF8eN8nETUWi88mqjZZ8L/DBQCA3/eKlDkaIiJ5dbePeGfLJxFdG4vPJtp5sghlVSaEBmpxS0wbucMhIpJVj7bWls9DeWUwW4TM0RBRc9ak4nPJkiWIi4uDTqdDQkICtm7des3tjUYj5s6di9jYWGi1WnTs2BEffvhhkwJuLjYcyAcAjLg5HAp2uRORk5zJo9u2bcOgQYMQHBwMPz8/dO3aFf/85z89GO31xYW0gp9aicoaM7IKy+UOh4iaMZWzO6xZswYzZ87EkiVLMGjQILz33nsYNWoUMjMz0a5du3r3GT9+PM6dO4fly5fjpptuQkFBAUwm0w0HLxezRSA18xwAYGT3CJmjIaKWxtk8GhAQgBkzZqBXr14ICAjAtm3b8NhjjyEgIACPPvqoDFdQl1IhoWtkIDJyLuLg2RLcFNZK7pCIqJmShBBO9Y/0798fffr0wdKlS+3runXrhrvuugspKSl1tt+wYQPuv/9+nDx5EkFBQU0KsrS0FAaDASUlJdDr9U06hivtyS7Gfct2IlCnQvpfRkCj4t0LRC2dJ/OMs3m0Pvfccw8CAgLw8ccfN2p7T1zfX9buxye7cvDYbR0wZ3Q3t5yDiJqvxuYZp6qm6upqpKenIzk52WF9cnIyduzYUe8+3377LRITE/GPf/wDbdu2RefOnfHss8+isrLSmVM3Kz8etHa5D+8axsKTiJzSlDx6tYyMDOzYsQNDhgxxR4hNZht0dICDjojoGpzqdi8sLITZbEZ4eLjD+vDwcOTn59e7z8mTJ7Ft2zbodDp8/fXXKCwsxBNPPIHi4uIG7/s0Go0wGo3230tLm8/UHUII/HiQXe5E1DRNyaM20dHROH/+PEwmE+bPn4+pU6c2uK0cebRnW2vx+dvpElgsgvfDE1G9mtRsJ0mOCUUIUWedjcVigSRJWLVqFfr164fRo0dj0aJFWLlyZYOtnykpKTAYDPYlJiamKWG6xZFzZcgproBWpeCz3ImoyZzJozZbt27Fnj17sGzZMixevBirV69ucFs58miXiEBoVQqUVZmQVcRBR0RUP6eKz5CQECiVyjrfzgsKCup8i7eJjIxE27ZtYTAY7Ou6desGIQROnz5d7z5z5sxBSUmJfcnNzXUmTLdKO3weAJDUMRj+GqfHaxGRj2tKHrWJi4tDz5498cgjj2DWrFmYP39+g9vKkUfVSoW99XNf7kW3n4+IWianik+NRoOEhASkpqY6rE9NTUVSUlK9+wwaNAhnz57FpUuXn/d79OhRKBQKREdH17uPVquFXq93WJqLtCPWieWHdQ2TORIiaomakkfrI4Rw6Fa/mlx5ND6mNQBgL4tPImqA093us2fPxgcffIAPP/wQhw4dwqxZs5CTk4Np06YBsH7bnjRpkn37iRMnIjg4GA8++CAyMzOxZcsWPPfcc3jooYfg5+fnuivxgNKqGqSfugAAGNqZxScRNY2zefTdd9/FunXrcOzYMRw7dgwrVqzA66+/jv/7v/+T6xIaZCs+2fJJRA1xut94woQJKCoqwsKFC5GXl4cePXpg/fr1iI2NBQDk5eUhJyfHvn2rVq2QmpqKJ598EomJiQgODsb48ePx8ssvu+4qPGT7sUKYLQIdQgPQLthf7nCIqIVyNo9aLBbMmTMHWVlZUKlU6NixI1555RU89thjcl1Cg3pHtwYAZOaVwmgyQ6tSyhsQETU7Ts/zKYfmMs/n81/sw+d7TuOhQXF4aezNssVBRK7XXPKMu3jq+oQQSHj5vygur8ba6YPQu7YllIi8n1vm+fRlQghsPmodbDSUo9yJiOolSRLio62DjvbmXJA5GiJqjlh8NtKhvDKcKzXCT61Ev7imPamJiMgX2O/7PM3J5omoLhafjWQb5Z7UMRg6Ne9hIiJqCAcdEdG1sPhsJHa5ExE1jm3Q0cnCcpRU1MgbDBE1Oyw+G6Gi2oSM2nuXBndi8UlEdC1tAjSIrZ0RZO/pi/IGQ0TNDovPRtiTfQE1ZoG2rf3sCZWIiBp2S23Xu21uZCIiGxafjbD9RCEAYGDH4Os+e5mIiIDE9taBmXuyi2WOhIiaGxafjbDzRBEAYNBNwTJHQkTUMvStLT4zci6ixmyRORoiak5YfF5HSUUN9p+xThcysEOIzNEQEbUMncJaweCnRmWNGZlnS+UOh4iaERaf17ErqwhCAB1CAxBh0MkdDhFRi6BQSEiMbQMA2M2udyK6AovP67B3uXdkqycRkTMu3/fJQUdEdBmLz+vYftw62CipI+/3JCJyRt/21pbPPaeKIYSQORoiai5YfF5DQVkVjhVcgiQBAzqw+CQickbPaAM0KgUKL1Uju6hC7nCIqJlg8XkNu05a71PqFqFHmwCNzNEQEbUsWpUS8dEGALzvk4guY/F5DbuzrMmSrZ5ERE3D+T6J6GosPq/B9k3ddt8SERE5x5Y/d3PQERHVYvHZgJLKGhw5VwYASGDxSUTUJAmxQVBIQFZhOfJLquQOh4iaARafDfg15wKEANoH+yMskPN7EhE1hcFPjR5trfd97qh9VDER+TYWnw2w3Z9ku1+JiIiaJql2nuQdtfMmE5FvY/HZANv9Sbzfk4joxtjmSd55oojzfRIRi8/6VJss2Jd7EQBbPomIblRi+zZQKyWcuViJU5zvk8jnsfisx4GzJTCaLAgK0KBDSIDc4RARtWj+GhVuaWftRWLXOxGx+KyH/X7P2DaQJEnmaIiIWj5b1zsHHRERi896XL7fk13uRESuYBt0tPNEESwW3vdJ5MtYfF5FCIH0U9bik/N7EhG5Ru+Y1vBTK1FUXo2jBWVyh0NEMmLxeZXc4koUl1dDo1SgR5RB7nCIiLyCRqVA3zhrb9L247zvk8iXsfi8SkautdXz5ig9NCr+8RARucrgm6xd75uPnpc5EiKSE6urq+zLLQFg7SIiIiLXGdY1FACw62QRKqvNMkdDRHJh8XmVfacvAgDiY9jlTkTkSh1DW6Ftaz9UmyzYeZKj3ol8FYvPK9SYLThwxtbyycFGRESuJEkShnaxtn6mHWbXO5GvYvF5hSP5ZTCaLNDrVGgf7C93OEREXmdYlzAAQNqRAj5qk8hHNan4XLJkCeLi4qDT6ZCQkICtW7c2ar/t27dDpVKhd+/eTTmt2+2tfaRmfExrTi5PRG7lTB796quvMGLECISGhkKv12PgwIH48ccfPRit6yTdFAyNUoHTFypx4ny53OEQkQycLj7XrFmDmTNnYu7cucjIyMDgwYMxatQo5OTkXHO/kpISTJo0CcOHD29ysO5me547BxsRkTs5m0e3bNmCESNGYP369UhPT8ewYcMwduxYZGRkeDjyG+evUaFf7ZRLm44UyBwNEclBEk72e/Tv3x99+vTB0qVL7eu6deuGu+66CykpKQ3ud//996NTp05QKpVYu3Yt9u7d2+hzlpaWwmAwoKSkBHq93plwnZL8z804eu4Slk9OxPBu4W47DxE1P57KM0DT8+iVunfvjgkTJuCll15q1PaevL7r+WDrSbz8/SHcelMIPpnaX9ZYiMh1GptnnGr5rK6uRnp6OpKTkx3WJycnY8eOHQ3ut2LFCpw4cQLz5s1z5nQeVVZVg2MFlwAAvaJbyxsMEXmtpubRK1ksFpSVlSEoqGU+Anho7X2fv2QV45LRJHM0RORpKmc2LiwshNlsRni4Y6tgeHg48vPz693n2LFjePHFF7F161aoVI07ndFohNFotP9eWlrqTJhNsv9MCYQA2rb2Q2ig1u3nIyLf1JQ8erU33ngD5eXlGD9+fIPbyJFHG6tjaADaB/sju6gCm4+cx5hekXKHREQe1KQBR1cPxhFC1DtAx2w2Y+LEiViwYAE6d+7c6OOnpKTAYDDYl5iYmKaE6RROLk9EntTYPHq11atXY/78+VizZg3CwsIa3E6OPNpYkiRhZPcIAMCPBxtXcBOR93Cq+AwJCYFSqazz7bygoKDOt3gAKCsrw549ezBjxgyoVCqoVCosXLgQ+/btg0qlwsaNG+s9z5w5c1BSUmJfcnNznQmzSQ6ctRafPaM5uTwRuY+zefRKa9aswcMPP4zPP/8cd9xxxzW3lSOPOiO5tvhMO1yAapNF5miIyJOcKj41Gg0SEhKQmprqsD41NRVJSUl1ttfr9di/fz/27t1rX6ZNm4YuXbpg79696N+//hvNtVot9Hq9w+JumWetXVLdo+S9EZ+IvJuzedRm9erVmDJlCj799FOMGTPmuueRI48645aY1ggL1KLMaMKOE3zaEZEvceqeTwCYPXs2HnjgASQmJmLgwIF4//33kZOTg2nTpgGwfts+c+YMPvroIygUCvTo0cNh/7CwMOh0ujrr5XTJaEJWoXW+ue5RbPkkIvdyJo8C1sJz0qRJePPNNzFgwAB7q6mfnx8MhpaZsxQKCSNuDseqn3Pw48Fz9kFIROT9nC4+J0yYgKKiIixcuBB5eXno0aMH1q9fj9jYWABAXl7edef8bG4O5VlbPSMNOgQFaGSOhoi8nbN59L333oPJZML06dMxffp0+/rJkydj5cqVng7fZUZ2j8Cqn3OQmnkOL9/VA0oFH+5B5AucnudTDu6en27l9izMX5eJO7qF4YPJfV1+fCJq/prTPJju0Byvr9pkQcLLqSirMuGLaQOR2L5lTh1FRFZumefTWx2svd/zZna5ExF5jEalwPCu1u72DQc46p3IV7D4xOXik4ONiIg863c9rHN8fvdbHiyWZt8RR0Qu4PPFp9FkxtFzZQBYfBIRedrQLqEI1KmQX1qFX7KL5Q6HiDzA54vPY+cuwWQRMPip0ba1n9zhEBH5FJ1aiVE9rHN+frP3rMzREJEn+HzxebB2cvnuUfpGPV2EiIhca1x8WwDADwfyOOE8kQ9g8cn7PYmIZDWwYzBCWmlxsaIGW4+dlzscInIzFp/24pMj3YmI5KBUSBgbbx14xK53Iu/n08Wn2SLsE8yz5ZOISD7j4qMAAKmZ51BuNMkcDRG5k08XnznFFaioNkOrUqBDaCu5wyEi8lm9Y1ojNtgflTVm/MA5P4m8mk8Xn0fyrVMsdQ4P5GPdiIhkJEkS/pAQDQD4fHeuzNEQkTv5dPFpm9+zc3igzJEQEdF9CTFQSMAv2cU4cf6S3OEQkZv4dPFpa/nsGsHik4hIbhEGHYZ0DgUAfL6HrZ9E3sq3i09byyeLTyKiZmFC33YAgC/Tz6DGzDk/ibyRzxafRpMZWYXlAIAu7HYnImoWhncLQ0grDQovGZF2uEDucIjIDXy2+DxRUA5z7WM1w/VaucMhIiIAaqUC9/SxDjz6jAOPiLySzxaftsFGXcID+VhNIqJm5P6+MQCAtCMFyCmqkDkaInI1ny0+D9umWYrg/J5ERM1Jh9BWuK1zKIQAPtqZLXc4RORiPlt82ls+I/hkIyKi5ubBpPYAgDV7cvnEIyIv47PFp22aJQ42IiJqfoZ0DkX7YH+UVZnwdcYZucMhIhfyyeKzrKoGZy5WAmDxSUTUHCkUEiYNbA8AWLkjG0IIeQMiIpfxyeLz6DnrkzMi9DoY/NUyR0NERPW5LzEaARoljhdcwtZjhXKHQ0Qu4pPFp/2Z7pxcnoio2dLr1PhDonXk+7LNJ2SOhohcxSeLT9tgIz5Wk4ioeXvktg5QKSTsOFGEjJwLcodDRC7gk8XnsQJr8XlTGKdZIiJqztq29sOdvdsCAJZsYusnkTfwyeLzRIH1sZosPomImr/Hh3aAJAGpmefsPVdE1HL5XPF5yWhCfmkVAKBjKItPIqLm7qawQIy8OQIAsIytn0Qtns8VnyfPW0e6hwZqYfDjSHciopbgiWEdAQBr957Bido8TkQtk88Vn7ak1TE0QOZIiIiosXpFt8Yd3cJgEcCi1KNyh0NEN8Dnis/jBbbik13uREQtyTPJXSBJwPe/5eHAmRK5wyGiJvK54tM22IjFJxFRy9ItUo9x8VEAgDd+OiJzNETUVL5XfNq63TnSnYioxZl1R2coFRLSjpzH7uxiucMhoiZoUvG5ZMkSxMXFQafTISEhAVu3bm1w26+++gojRoxAaGgo9Ho9Bg4ciB9//LHJAd8Ik9mC7CJbyyfv+SQi+TiTR/Py8jBx4kR06dIFCoUCM2fO9FygzUz7kACMr33q0cvfZcJi4TPfiVoap4vPNWvWYObMmZg7dy4yMjIwePBgjBo1Cjk5OfVuv2XLFowYMQLr169Heno6hg0bhrFjxyIjI+OGg3dW7oVK1JgF/NRKRBn8PH5+IiLA+TxqNBoRGhqKuXPnIj4+3sPRNj+zRnRCK60K+06X4KuMM3KHQ0ROkoQQTn1t7N+/P/r06YOlS5fa13Xr1g133XUXUlJSGnWM7t27Y8KECXjppZcatX1paSkMBgNKSkqg1+udCdfBfzPPYepHe9A9So/vnxrc5OMQkfdxVZ5pjBvJo0OHDkXv3r2xePFip87pyevzhPe3nMDf1x9GaKAWac8ORSutSu6QiHxeY/OMUy2f1dXVSE9PR3JyssP65ORk7Nixo1HHsFgsKCsrQ1BQkDOndonj5znSnYjk5Yo8SsCUpDjEhQTgfJkR72w8Lnc4ROQEp4rPwsJCmM1mhIeHO6wPDw9Hfn5+o47xxhtvoLy8HOPHj29wG6PRiNLSUofFFU5wmiUikpkr8mhjuCuPNhcalQJ/GdMNALB820kcL+BjN4laiiYNOJIkyeF3IUSddfVZvXo15s+fjzVr1iAsLKzB7VJSUmAwGOxLTExMU8Ks4/JIdw42IiJ5NTWPNpa78mhzcnvXMNzeNQw1ZoEXv9zPwUdELYRTxWdISAiUSmWdb+cFBQV1vsVfbc2aNXj44Yfx+eef44477rjmtnPmzEFJSYl9yc3NdSbMegkhcOI85/gkInndSB51hjvyaHMjSRL+dlcPBGiU2HPqAlb9Uv+ALSJqXpwqPjUaDRISEpCamuqwPjU1FUlJSQ3ut3r1akyZMgWffvopxowZc93zaLVa6PV6h+VGFZVXo6SyBpIExIWw5ZOI5NHUPOosd+TR5qhtaz88N7ILAODVHw4jr6RS5oiI6Hqc7nafPXs2PvjgA3z44Yc4dOgQZs2ahZycHEybNg2A9dv2pEmT7NuvXr0akyZNwhtvvIEBAwYgPz8f+fn5KCnx7KPRbPd7Rrfxg06t9Oi5iYiu5GweBYC9e/di7969uHTpEs6fP4+9e/ciMzNTjvCbnQcGtsct7VrjktGEF7/cDycncSEiD3N6booJEyagqKgICxcuRF5eHnr06IH169cjNjYWgHUy5CvnqnvvvfdgMpkwffp0TJ8+3b5+8uTJWLly5Y1fQSOdKqoAAMSFsMudiOTlbB4FgFtuucX+Oj09HZ9++iliY2ORnZ3tydCbJaVCwj/u7YUxb2/D5qPn8dHOU5ic1F7usIioAU7P8ykHV8xP9+qGw1i66QQmDYzFwjt7uDhCImrpvG0ezKt5+/UBwMrtWZi/LhNalQLfPXkrOoUHyh0SkU9xyzyfLdmp2sdqxgbzfk8iIm80Oak9hnQOhdFkwVOf7UVVjVnukIioHj5TfGYXWrvd2wf7yxwJERG5gyRJeO0PvRAUoMGhvFIsWHdQ7pCIqB4+UXwKIewtn+050p2IyGuFBeqweEJvSBKw+pdcrNnN6ZeImhufKD7PXzKivNoMhWQd7U5ERN7rts6heGZEZwDAX785iN9OX5Q3ICJy4BPFp22ke1RrP2hVnGaJiMjbPTH0JtzRLQzVJgse/Sgd+SVVcodERLV8ovjMLrR2uXNyeSIi36BQSHhjfG90DA1AfmkVHlq5G5eMJrnDIiL4SvFpH+nOwUZERL7C4KfGygf7IaSVBpl5pXjy019hMlvkDovI5/lI8Wkb6c6WTyIiXxIT5I8PJveFTq1A2pHzeP6L32CxNPvprYm8mk8Un5zjk4jId/WOaY23/9gHSoWErzLO4K/fHOAjOIlk5PXFpxDCPsdnXAi73YmIfNGIm8OxaHw8JAlY9XMO/vbdIRagRDLx+uKzqLwal4wmSBIQ3YbFJxGRr7qzd1u8ek8vAMCH27Mwd+0BmNkFT+RxXl982rrcowx+0Kk5zRIRkS8b3zcGr9zTE5IEfPpzDmZ/vhc1HIRE5FFeX3xm1Xa5c6Q7EREBwP392uHN+2+BSiHhm71n8fC/96CsqkbusIh8htcXn3ysJhERXW1cfBTeeyABfmolthw9j/uW7sTpCxVyh0XkE7y++Lw8zRJbPomI6LLh3cLx+WMDERaoxZFzZbjr3e3YeaJI7rCIvJ7XF5+2ls92QWz5JCIiRz2jDVg7fRC6RepReKkaf/pgF5ZuOsGR8ERu5PXFZ04x7/kkIqKGRbX2w1ePJ+GePm1hEcCrGw5jyordOFfK58ETuYNXF59lVTW4WGG9iTwmiMUnERHVz0+jxBt/iMff7+4JjUqBzUfPI/mfW/DN3jNsBSVyMa8uPk9fqAQAtPFXo5VWJXM0RETUnEmShIn92+H7J29Fz7YGlFTW4OnP9mLGpxkovGSUOzwir+HVxWdubZc7J5cnIqLG6hQeiK+eSMKsOzpDpZDw/f48DHttE5Zvy+KcoEQu4N3FZ23LZ0yQn8yREBFRS6JWKvD0HZ3w9ROD0LOtAWVGE/72XSZGvbkVW4+dlzs8ohbNu4vP2pbPGLZ8EhFRE9hGw79yT08EBWhwvOASHlj+C/70wS7syS6WOzyiFsmri0/bhMHRHGxERERNpFRIuL9fO6Q9MxRTktpDrZSw/XgR7lu2Ew8s/xnppy7IHSJRi+LlxWdtt3sbdrsTEdGNMfirMX9cd2x8Ziju7xsDlULC1mOFuHfpDtyzZDvW7TvLe0KJGsFrh4ALIS53u7Plk4iIXCQmyB+v3NsL04fdhLc3HsPXGWfwa85F/JqTgUiDDn/q3w73JkQj0sCGD6L6eG3L54WKGpRXmwEAbVszARARkWvFBPnjH/fFY/uLt+Pp4Z0Q0kqDvJIqvP7TUSS9shEPLP8Z3+w9g8ra/xcRkZXXtnzaWj3DArXQqZUyR0NERN4qLFCHWSM644lhHfHdvjys2ZOLX7KKsfVYIbYeK4SfWolhXUMxsnsEbu8ahkCdWu6QiWTlvcXnBXa5ExGR52hVStybEI17E6KRU1SBL389jS9/PY3TFyqxfn8+1u/Ph0apQL+4INzaKQS33hSCmyP1UCgkuUMn8ijvLT6LOdiIiIjk0S7YH7NGdMbMOzrhwJlSbDiYhx8O5OPk+XJsO16IbccLAQDBARok3RSCgR2C0Se2NTqFBULJYpS8nNcWn6fZ8klERDKTJAk9ow3oGW3AcyO74njBJWw5eh7bjhdi18kiFJVXY92+s1i37ywAoJVWhfgYA/q0a4Ne0a3RLTIQbVv7QZJYkJL38Nri0/Z0o2i2fBIRUTNxU1gr3BTWCg/dGodqkwUZORew/Xgh9py6gH25F3HJaML240XYfrzIvk+gToVuEXp0jQxE5/BAdAgJQPuQAETodeyypxapScXnkiVL8NprryEvLw/du3fH4sWLMXjw4Aa337x5M2bPno2DBw8iKioKzz//PKZNm9bkoBvjNJ9uRETNWEvIo+ReGpUC/TsEo3+HYACA2SJw9FwZfs25gF9PXcTBsyU4cf4SyqpM+CW7GL9c9UQlrUqB2GB/tA8OQFxIANq28UOEXodIgx8iDDoEB2hYnFKz5HTxuWbNGsycORNLlizBoEGD8N5772HUqFHIzMxEu3bt6myflZWF0aNH45FHHsEnn3yC7du344knnkBoaCjuvfdel1zE1SwWcXmCeXa7E1Ez0xLyKHmeUiGhW6Qe3SL1+FP/WABAtcmCk4WXcDivDIfyS3Hs3CVkF5Yjp7gCRpMFR89dwtFzl+o9nlopIVyvQ6RBh3C9tRgNCtAiqJUGQf4aBAVoENxKgzb+GrTxV0Ol9NrZF6mZkYQQwpkd+vfvjz59+mDp0qX2dd26dcNdd92FlJSUOtu/8MIL+Pbbb3Ho0CH7umnTpmHfvn3YuXNno85ZWloKg8GAkpIS6PX6626fX1KFASn/g1Ih4cjffsd/UER0Xc7mmRvREvIoNW8mswVnLlYiq7Ac2YXlyC6qQF5JJfJLqpBXUoXzl4xw7v/ugJ9aiUCdCq10KgTq1AjUqhCosy6ttGr7a51aCT+10vpTo4BOpYROc8U6tRI6tQI6tRJalYL3q/qQxuYZp1o+q6urkZ6ejhdffNFhfXJyMnbs2FHvPjt37kRycrLDupEjR2L58uWoqamBWu36+c5sg40iDToWnkTUrLSUPErNm0qpQGxwAGKDA4Audd+vMVtQUGZEfkkl8kqqkF9ShQsV1Sguv7wUlVfjQnk1LlbWQAigssaMyhozCsqMro1VIUGtVECttP1UQK2SoFZcfq1SKKBRKqC6chulBIVCglKSoJBgf61UWNcrJFjfu2K9JElQKhzXW7e1rldI1teSBEiwDgizvYYk1a4DJFyx/op1uHI/2/ortsfV79VzHOCq89vXX1ZfvS7h6o3q2+bq40jXfL++c119nn5xQdCoXFtLOVV8FhYWwmw2Izw83GF9eHg48vPz690nPz+/3u1NJhMKCwsRGRlZZx+j0Qij8fKHv7S01JkwL8/xyfs9iaiZaSl5lFo2tVKBtq39GvWEP5PZgrIqk3Ux1thfX7ritXWpwSWjCVU1ZlTWWFBVY7a+rjajymRGZbUFxtoC1mS53OxqsgiYLGZU1rjzisldfv3rCASpNC49ZpMGHF1dSQshrtmsXt/29a23SUlJwYIFC5oSGgBAo1Ti5kg9ukQENvkYRETu1NzzKPkOlVKBNgEatAlwXYFRY7YWp0aTBSazQI3ZUrs0/NpktqD6itc1ZgvMFgGzsH7era8FLBYBi7AO0LJcsV7UrrtyvUVYx4Fc3u/y8QQACEDAuq+48jVwxW0LV64T9vdsvwNX7XvVcSDsR2nwOA7quV/i6jX13VIhrtrq6m3q3+fqbepupHTDbRNOFZ8hISFQKpV1vp0XFBTU+VZuExERUe/2KpUKwcHB9e4zZ84czJ492/57aWkpYmJiGh3nmF6RGNOrbksAEZHcWkoeJboRtq5zNgFRfZzqxNdoNEhISEBqaqrD+tTUVCQlJdW7z8CBA+ts/9NPPyExMbHB+5S0Wi30er3DQkTkDZhHicjXOX0H6ezZs/HBBx/gww8/xKFDhzBr1izk5OTY55ubM2cOJk2aZN9+2rRpOHXqFGbPno1Dhw7hww8/xPLly/Hss8+67iqIiFoQ5lEi8mVO3/M5YcIEFBUVYeHChcjLy0OPHj2wfv16xMZa5yTLy8tDTk6Offu4uDisX78es2bNwrvvvouoqCi89dZbnJuOiHwW8ygR+TKn5/mUA+enIyJ38/Y84+3XR0Tya2ye4SSYREREROQxLD6JiIiIyGNYfBIRERGRxzRpknlPs92Wyid0EJG72PJLC7gNvkmYR4nI3RqbR1tE8VlWVgYAnCCZiNyurKwMBoNB7jBcjnmUiDzlenm0RYx2t1gsOHv2LAIDA6/5+Lkr2Z7mkZub65MjO3n9vH5ev3PXL4RAWVkZoqKioFB43x1JzKPO4/Xz+nn97smjLaLlU6FQIDo6ukn7+vqTPXj9vH5ef+Ov3xtbPG2YR5uO18/r5/W7No9639d7IiIiImq2WHwSERERkcd4bfGp1Woxb948aLVauUORBa+f18/r993rdxVf/3Pk9fP6ef3uuf4WMeCIiIiIiLyD17Z8EhEREVHzw+KTiIiIiDyGxScREREReQyLTyIiIiLyGK8sPpcsWYK4uDjodDokJCRg69atcofkEfPnz4ckSQ5LRESE3GG5zZYtWzB27FhERUVBkiSsXbvW4X0hBObPn4+oqCj4+flh6NChOHjwoDzBusH1rn/KlCl1Pg8DBgyQJ1g3SElJQd++fREYGIiwsDDcddddOHLkiMM23v4ZcCdfzaMAcylz6VqH9705l8qVR72u+FyzZg1mzpyJuXPnIiMjA4MHD8aoUaOQk5Mjd2ge0b17d+Tl5dmX/fv3yx2S25SXlyM+Ph7vvPNOve//4x//wKJFi/DOO+9g9+7diIiIwIgRI+zPuG7prnf9APC73/3O4fOwfv16D0boXps3b8b06dOxa9cupKamwmQyITk5GeXl5fZtvP0z4C6+nkcB5tIrefu/I1/OpbLlUeFl+vXrJ6ZNm+awrmvXruLFF1+UKSLPmTdvnoiPj5c7DFkAEF9//bX9d4vFIiIiIsQrr7xiX1dVVSUMBoNYtmyZDBG619XXL4QQkydPFnfeeacs8cihoKBAABCbN28WQvjeZ8CVfDmPCsFcylz6tcM6X8qlnsqjXtXyWV1djfT0dCQnJzusT05Oxo4dO2SKyrOOHTuGqKgoxMXF4f7778fJkyflDkkWWVlZyM/Pd/gsaLVaDBkyxGc+CwCwadMmhIWFoXPnznjkkUdQUFAgd0huU1JSAgAICgoCwM9AUzGPWjGXWvHfkZWv5FJP5VGvKj4LCwthNpsRHh7usD48PBz5+fkyReU5/fv3x0cffYQff/wR//rXv5Cfn4+kpCQUFRXJHZrH2f6+ffWzAACjRo3CqlWrsHHjRrzxxhvYvXs3br/9dhiNRrlDczkhBGbPno1bb70VPXr0AMDPQFP5eh4FmEuvxH9HvpNLPZlHVU0Ps/mSJMnhdyFEnXXeaNSoUfbXPXv2xMCBA9GxY0f8+9//xuzZs2WMTD6++lkAgAkTJthf9+jRA4mJiYiNjcX333+Pe+65R8bIXG/GjBn47bffsG3btjrv+fJn4Eb48p8bc2ldvvx58JVc6sk86lUtnyEhIVAqlXWq8YKCgjpVuy8ICAhAz549cezYMblD8TjbyFR+Fi6LjIxEbGys130ennzySXz77bdIS0tDdHS0fT0/A03DPFoXcyn/HV3JG3Opp/OoVxWfGo0GCQkJSE1NdVifmpqKpKQkmaKSj9FoxKFDhxAZGSl3KB4XFxeHiIgIh89CdXU1Nm/e7JOfBQAoKipCbm6u13wehBCYMWMGvvrqK2zcuBFxcXEO7/Mz0DTMo3Uxl/Lf0ZW8KZfKlkebPFSpmfrss8+EWq0Wy5cvF5mZmWLmzJkiICBAZGdnyx2a2z3zzDNi06ZN4uTJk2LXrl3i97//vQgMDPTaay8rKxMZGRkiIyNDABCLFi0SGRkZ4tSpU0IIIV555RVhMBjEV199Jfbv3y/++Mc/isjISFFaWipz5K5xresvKysTzzzzjNixY4fIysoSaWlpYuDAgaJt27Zec/2PP/64MBgMYtOmTSIvL8++VFRU2Lfx9s+Au/hyHhWCuZS51HdyqVx51OuKTyGEePfdd0VsbKzQaDSiT58+9ikDvN2ECRNEZGSkUKvVIioqStxzzz3i4MGDcoflNmlpaQJAnWXy5MlCCOsUEfPmzRMRERFCq9WK2267Tezfv1/eoF3oWtdfUVEhkpOTRWhoqFCr1aJdu3Zi8uTJIicnR+6wXaa+awcgVqxYYd/G2z8D7uSreVQI5lLmUt/JpXLlUan25EREREREbudV93wSERERUfPG4pOIiIiIPIbFJxERERF5DItPIiIiIvIYFp9ERERE5DEsPomIiIjIY1h8EhEREZHHsPgkIiIiIo9h8UlEREREHsPik4iIiIg8hsUnEREREXkMi08iIiIi8pj/D06DhjtV0hC9AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x300 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = 7 # m^2, sail area\n",
    "rho = 1.2 # air density kg/m^3\n",
    "v0 = 2 # m/s mean wind speed \n",
    "beta = 1/v0 * np.sqrt(2/(rho*A))\n",
    "from scipy.stats import expon\n",
    "velocity_dist = expon(scale=v0) # input velocity is exponentially distributed\n",
    "\n",
    "xvals = np.linspace(0.001, 20, 200)\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(8, 3))\n",
    "ax[0].plot(xvals, velocity_dist.cdf(xvals))\n",
    "ax[0].set_title('Wind velocity CDF')\n",
    "ax[1].plot(xvals, velocity_dist.pdf(xvals))\n",
    "ax[1].set_title('Wind velocity PDF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0d0bbbf",
   "metadata": {},
   "source": [
    "In class we saw that the CDF and PDF of wind force can be computed analytically as:\n",
    "$$ F_{Y}(y)=1-e^{-\\beta \\sqrt{y}} $$\n",
    "$$ f_{Y}(y)=\\frac{\\beta}{2 \\sqrt{y}} e^{-\\beta \\sqrt{y}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "80751939",
   "metadata": {},
   "outputs": [],
   "source": [
    "# those were computed in class\n",
    "def force_cdf(y):\n",
    "    cdf_vals = np.where(y>0, 1.-np.exp(-beta*np.sqrt(y)), 0.)\n",
    "    return cdf_vals\n",
    "\n",
    "def force_pdf(y):\n",
    "    pdf_vals = np.where(y>0, beta/(2*np.sqrt(y))*np.exp(-beta*np.sqrt(y)), 0.)\n",
    "    return pdf_vals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "7ce635c3",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqAAAAEnCAYAAAB2YfE4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMaklEQVR4nO3deVxU5f4H8M/swy6LssgimgtFmkIamJlZkJYtt37SctNSuxGWCnlLs65Lt7hW18wKrZS8lan3upTduAaVW0mluORCZgmiLCKgbMIMzJzfH8OMDDMoA7Mww+f9ep0XM895zjnPo/X45dmOSBAEAUREREREdiJ2dAGIiIiIqGdhAEpEREREdsUAlIiIiIjsigEoEREREdkVA1AiIiIisisGoERERERkVwxAiYiIiMiuGIASERERkV0xACUiIiIiu2IA6sI2bdoEkUiEjRs3mpwbNmwYRCIRvv76a5NzAwYMwIgRIwAAO3fuhEgkws6dO61WrsLCQohEIqxdu/aqeb/99lvExsbCw8MDIpEIn3/+udXKYQ/nzp3DvHnzcP3118PT0xNKpRIDBw7E7NmzcfLkSUO+RYsWQSQSGQ53d3eEhoYiMTER77zzDmpra03u/fjjjxtd0/r473//a89qErk0tqWOo6+j/hCLxfD398fEiRORm5trlLd1PolEAl9fXwwbNgxPPfUUfvzxx6veu/URGxtrryr2WFJHF4Bs59Zbb4VIJMKOHTuQlJRkSK+qqsKRI0fg4eGBHTt2IDEx0XDu7NmzOHXqFNLS0gAAI0aMQG5uLq699lq7l18QBEyePBmDBg3Ctm3b4OHhgcGDB9u9HJ31888/4+6774YgCHjmmWcQFxcHuVyOEydO4NNPP8XIkSNx4cIFo2u2b98OHx8fqNVqlJSU4Ntvv8Xzzz+PN954A19++SWGDRtmlN/NzQ3fffedybOHDBli07oR9SRsSx3v2WefxSOPPAKNRoNjx45h8eLFGDduHHJzczF8+HBDvgcffBDPPfccBEFATU0Njh49io8//hgffPABZs2ahbfffrvde7fm6elp8zr1dAxAXVhAQACio6NNfuPetWsXpFIppk+fjh07dhid038fN24cAMDb2xs33XSTXcrbVklJCaqqqnD//fdj/PjxVrlnQ0MDlEolRCKRVe7XnpqaGtx7771QKpXYu3cvQkNDDeduvfVWPPXUU9i0aZPJdTExMQgICDB8f+ihh/DMM89g7NixuOeee/Dbb79BoVAYzovFYof9/RD1FGxLTdmrLdULDw83/PmNHj0a11xzDcaPH4+MjAx8+OGHhnyBgYFGf86JiYmYM2cO/vKXv2DFihUYMmQInn766XbvTfbDIXgXN27cOJw4cQKlpaWGtJ07d+LGG2/ExIkTkZeXZzS8u3PnTkgkEowZM8bwve2w0eOPPw5PT0/8/vvvmDhxIjw9PREWFobnnnsOKpXK6PklJSWYPHkyvLy84OPjg6SkJJSVlV213IsWLTIEbS+88AJEIhH69etnOP/9999j/Pjx8PLygru7O+Lj4/HVV18Z3WPt2rUQiUTIzs7GtGnT0Lt3b7i7uxvK+NlnnyEuLg6enp7w9PTEDTfcgDVr1hjd45tvvsH48ePh7e0Nd3d3jB49Gt9+++1Vy//hhx+irKwMr7/+ulHw2dqDDz541fsAuiG+BQsWoKioyOwQIBHZHttSx7Sl7dEHjKdPn75qXolEgnfffRcBAQF44403Ov1Msi4GoC5O/9t360Zvx44dGDt2LEaPHg2RSIQ9e/YYnRsxYgR8fHyueN+mpibcc889GD9+PL744gtMmzYNb731FpYuXWrI09DQgNtvvx3Z2dlIT0/Hf/7zHwQFBRkNYbVnxowZ2LJlCwDd8Ehubi62bt0KQNfrcNttt6G6uhpr1qzB+vXr4eXlhUmTJpkN0KZNmwaZTIZPPvkEmzZtgkwmw9/+9jc8+uijCAkJwdq1a7F161ZMnTrVqDH79NNPkZCQAG9vb/zrX//Cv//9b/j5+SExMfGqDWd2djYkEgkmTZp01bp2xD333AMA2L17t8m55uZmo0Oj0VjlmUR0GdtSx7Sl7fn9998BAL179+5Qfjc3N9x+++0oKCjA2bNnjc5ptVqTdlQQhE6ViywgkEurqqoSxGKx8Je//EUQBEGoqKgQRCKRsH37dkEQBGHkyJHC3LlzBUEQhKKiIgGA8Pzzzxuu37FjhwBA2LFjhyFt6tSpAgDh3//+t9GzJk6cKAwePNjwfeXKlQIA4YsvvjDK9+STTwoAhI8++uiKZS8oKBAACG+88YZR+k033ST06dNHqK2tNaQ1NzcL0dHRQmhoqKDVagVBEISPPvpIACBMmTLF6PpTp04JEolEePTRR9t9dn19veDn5ydMmjTJKF2j0QjDhg0TRo4cecWyDxkyRAgKCrpintYWLlwoABDOnz9v9nxDQ4MAQJgwYYIhTf/30PYYPXp0h59LRB3DttQxbam+7EuXLhWampqExsZGIS8vT7jxxhsFAMJXX31lyAtAmDlzZrv3euGFFwQAwk8//WR0b3NHTk7OFctFXcceUBenXwWo/619165dkEgkGD16NABg7NixhrlKbecsXYlIJDLp3Rs6dKjRb707duyAl5eXofdOr+1kb0vU19fjp59+woMPPmg0SVwikeCxxx7D2bNnceLECaNrHnjgAaPvOTk50Gg0mDlzZrvP2bt3L6qqqjB16lSj34q1Wi3uvPNO7Nu3D/X19Z2uh6WEdn4bd3Nzw759+4yOtkNfRNR1bEsd25a+8MILkMlkUCqViImJQVFREd5//31MnDixw3Vurx2dPXu2STs6atSoDt+XOoeLkHqAcePGYdmyZSgpKcGOHTsQExNjaHDGjh2Lf/7zn6iursaOHTsglUpx8803X/We7u7uUCqVRmkKhQKNjY2G75WVlQgMDDS5NigoqNN1uXDhAgRBQHBwsMm5kJAQw3Nba5v3/PnzANDu3ExAt30ScOV5mlVVVfDw8DB7Ljw8HCdPnkR9fX27eSyh/8dIX0c9sVjM7UKI7IRtqf3bUr3Zs2fjz3/+M8RiMXr16oXIyEiLF0C1146GhoayHXUABqA9gL7R3LlzJ3bu3Gn0G6O+gdy9e7dhQr21tp/w9/fHzz//bJLekYnz7fH19YVYLDZaCKBXUlICAEaryAGYNFL6OUNnz55FWFiY2efo7/HOO++0uzrS3D8IeomJicjOzsaXX36Jhx56qN18HbVt2zYAuhX0ROQYbEvt35bqdTVIbGhowDfffIMBAwZcMWAm++EQfA9wyy23QCKRYNOmTTh27JhREOPj44MbbrgB//rXv1BYWNihIaOOGjduHGpraw3Bk95nn33W6Xt6eHhg1KhR2LJlCxoaGgzpWq0Wn376KUJDQzFo0KAr3iMhIQESiQQrV65sN8/o0aPRq1cvHD9+HLGxsWYPuVze7vXTp09HUFAQnn/+eRQXF5vNo18YcDWHDx/Ga6+9hn79+mHy5MkduoaIrI9tqTF7tKXWoNFo8Mwzz6CyshIvvPCCTZ9FHcce0B7A29sbI0aMwOeffw6xWGyYs6Q3duxYLF++HEDH5ix11JQpU/DWW29hypQpePXVVzFw4EBkZWWZfWOIJdLT03HHHXdg3LhxmDt3LuRyOTIyMnD06FGsX7/+qsMy/fr1w4svvohXXnkFDQ0NePjhh+Hj44Pjx4+joqICixcvhqenJ9555x1MnToVVVVVePDBB9GnTx+cP38ehw8fxvnz56/Y6Pr4+OCLL77A3XffjeHDhxttRH/y5El8+umnOHz4MP70pz8ZXZeXlwcfHx80NTUZNqL/5JNP0KdPH3z55Zc2b6iJqH1sS43Zoy211Llz5/Djjz9CEATU1tYaNqI/fPgwUlNT8eSTT1rtWdQ1DEB7iHHjxmHfvn0YPnw4vL29jc6NHTsWb731FuRyOeLj4632THd3d3z33XeYPXs25s2bB5FIhISEBGzYsKFLzxk7diy+++47LFy4EI8//ji0Wi2GDRuGbdu24e677+7QPZYsWYKBAwfinXfewaOPPgqpVIqBAwdi1qxZhjx//vOfER4ejtdffx1PPfUUamtr0adPH9xwww14/PHHr/qMkSNH4siRI3jrrbfw73//G0uXLoVGo0FYWBjGjx+Pd9991+SaO++8E4BuDpifnx+uv/56LF26FE888QS8vLw69gdERDbDttSYPdpSS2zatAmbNm2CWCyGp6cnIiIiEBcXh1WrVnGz+W5GJLS3LIyIiIiIyAY4B5SIiIiI7IoBKBERERHZFQNQIiIiIrIrBqBEREREZFcMQImIrCwjIwORkZGG1wbu2bOn3bylpaV45JFHMHjwYIjFYsyZM8dsvs2bN+Paa6+FQqHAtddei61bt9qo9EREtscAlIjIijZu3Ig5c+ZgwYIFOHjwIMaMGYMJEyagqKjIbH6VSoXevXtjwYIFGDZsmNk8ubm5SEpKwmOPPYbDhw/jsccew+TJk/HTTz/ZsipERDbjFNswabValJSUwMvLy+J3vxIRdYR+4+qQkBCIxZ3/3XzUqFEYMWKE0ebaUVFRuO+++5Cenn7Fa2+99VbccMMNhs3M9ZKSklBTU4P//e9/hrQ777wTvr6+WL9+fYfKxXaUiGzNknbUKTaiLykpafc9s0RE1nTmzJlOvytarVYjLy8P8+bNM0pPSEjA3r17O12m3NxcpKamGqUlJiaaBKqtqVQqqFQqw/fi4mJce+21nS4DEVFHdaQddYoAVP8GmDNnzpi8eYKIyBpqamoQFhbWpTdOVVRUQKPRIDAw0Cg9MDAQZWVlnb5vWVmZxfdMT0/H4sWLTdLZjhKRrVjSjjpFAKofLvL29mbDSUQ2ZY3h6bb3EAShy/e19J7z589HWlqa4bv+Hwa2o0Rkax1p75wiACUicgYBAQGQSCQmPZPl5eUmPZiWCAoKsvieCoUCCoWi088kIrIlroInIrISuVyOmJgY5OTkGKXn5OQgPj6+0/eNi4szuWd2dnaX7klE5EgWB6C7d+/GpEmTEBISApFIhM8///yq1+zatQsxMTFQKpXo378/Vq1a1ZmyEhF1e2lpaVi9ejUyMzORn5+P1NRUFBUVITk5GYBuaHzKlClG1xw6dAiHDh1CXV0dzp8/j0OHDuH48eOG87Nnz0Z2djaWLl2KX3/9FUuXLsU333zT7p6hRETdncVD8PX19Rg2bBieeOIJPPDAA1fNX1BQgIkTJ+LJJ5/Ep59+ih9++AEpKSno3bt3h64nInImSUlJqKysxJIlS1BaWoro6GhkZWUhIiICgG7j+bZ7gg4fPtzwOS8vD5999hkiIiJQWFgIAIiPj8eGDRvw0ksv4eWXX8aAAQOwceNGjBo1ym71IiKypi7tAyoSibB161bcd9997eZ54YUXsG3bNuTn5xvSkpOTcfjwYeTm5nboOTU1NfDx8UF1dTUnzxORTbh6O+Pq9SMix7OknbH5HNDc3FwkJCQYpSUmJmL//v1oamqy9eOJiIiIqJux+Sr49vava25uRkVFBYKDg02uabuBck1Nja2LSUTdgCAI0GgFNGsFaAXdT41GgKZ1estPjVYLjRZo1mrNnBMgFYswqr+/o6vktGoam/DLmWooZGLc2M/P0cUhIhdjl22YzO1fZy5dr70NlImoY4SW4E3VrIW65VA1a1p+ag3p+jS1Rve9SaNFk0ZAs0aLZq0AtUaL5pbvTdqWnxoBTS3pTdqWn/rrjL7r7mF0v2atIU/bgLFZq4XWii8GDvJW4scXx1vvhj1Mwfl6/HnNT+jbyw0/zLvN0cUhIhdj8wC0vf3rpFIp/P3N9060t4EykbMSBF0g1qDWoKGp5Wj53Njqc+vvl9Qtn9vkbWwyDiaNg0qN7qdGi87P7u6epGIRxGIRpGIRJC2H4bNIBIlEBKlYDLEIkIrFCPCSO7rITk3fP6B1tf+QiKhbsHkAGhcXhy+//NIoLTs7G7GxsZDJZGav4QbK1F00a7SoaWxGdUMTahqaUK9qRl3LofusQZ2qCfUqTau0y+frVRrUNjahXq2BxprdexaSikVQSMWQtxwKqaTlZ0ua5PJPqUQEqUQMmbjlp0QMWUtwJ5OIIGvJI5OIIW3JI2+5Rio2Pq+/7vL3y3kk7QWThu9iw2exyDpvKKKOE7f8eTP+JCJbsDgAraurw++//274XlBQgEOHDsHPzw/h4eGYP38+iouL8fHHHwPQrXh/9913kZaWhieffBK5ublYs2YN1q9fb71aEF2FulmLqno1KutVqKxT48IlNWoamlBt9mg2nKtTNVu9LBKxCO4yCZRyCdxkukP3Waz7LpfATSaFm1zc5nzLZ5nEEDjqA8m2waRCJoZCIoFCJjYEe0SWYA8oEdmSxQHo/v37MW7cOMN3/VD51KlTsXbtWpM97iIjI5GVlYXU1FS89957CAkJwYoVK7gHKHVZbWMTztWocK6mEedrVaisV6OyTtUSaBp/rm3sWiDpqZDCWymFp1IKD4UUngopPOS6754KKTwUEngopPBS6M4b8rT89FRIWwJLXcBI1N3pe0Ad2HFPRC7M4gD01ltvxZW2Dl27dq1J2tixY3HgwAFLH0U9lCAIOF+nwtkLDSirbkRZdSPO1TbiXHWjIeA8V9OIerXGovtKxSL4esjh7yGHr7scPm4y3eGu++mt/97m8FZKIZUwaKSeRd8D2oWtoomI2mWXVfBErQmCgMp6Nc5UXcLZCw04e6EBZy7oP19C8YUGqJq1HbqXl1KKQG8l+ngp4O+pgL+HHH4ecvh76gJNf08F/DzkCPBQwNtNynmERB1kmAPq4HIQkWtiAEo209ikQUFFPU6dr8ep83U4VXH559WGxMUi3TY6Ib3cEOitbDkUCPJRoo+XsuWnAh4K/idMZAtizgElIhviv97UZc0aLQor65FfWotfy2rwa2ktTpyrRfHFhnZX0IpaAsxQXzeE+roj1NcNYS0/Q33dEeSj5FxJIgfSjxZoOQmUiGyAAShZRNWsQX5pLQ4VXcCxkhrkl9Xgt3N1ULczZO6tlKJ/b0/07+2BAb090T/AA/17eyLC3x1KmcTOpSeijtJPVmEHKBHZAgNQapcgCCisvIRDZy7g8JlqHDxzEfklNVBrTINNd7kEg4O8MCTIG1HBXhgc6IUBfTzh7yHnvEsiJ8Q5oERkSwxAyUAQBPxxvh4/nqpsOapQUacyyefrLsOwsF4YGtoL1wbrgs5wP3eIudckkcu4vA0TQ1Aisj4GoD1ceU0jdpwox56TFWYDTrlUjOgQb9wQ5othYT4YHuaLMD839moSuThuRE9EtsQAtIfRaAUcPnsRO34tx3e/luNYSY3ReblUjJhwX9zU3x839ffDsLBenKtJ1ANd3gfUseUgItfEALQHaNZo8VNBFf77Swm+PnYOVfVqwzmRCBga2gu3DuqN+AH+DDiJCADfBU9EtsUA1EVptYIh6Nx+tAyVrYJOb6UUtwzqjXGD+2Ds4N4I8FQ4sKRE1B1xDigR2RIDUBdTcrEB/9l/Fv/efwbFFxsM6X4ectwZHYS7rg/GqEg/vlqSiK6IG9ETkS0xAHUBzRotvsk/h/U/n8Huk+cNQ2beSikmRAfj7mHBiOvvz6CTiDpOPwfUsaUgIhfFANSJ1TQ2YePPZ7B2b6FRb2f8AH8k3RiGxOuCOJ+TiDql9RxQQRC48wURWRUDUCdUVt2ID3afwsZ9RahXawDohtgfHhmGybFhiPD3cHAJicjZiVsFnIJweVU8EZE1MAB1IudqGrFy5x/47Ociw6svB/bxxPSbI3Hf8L7s7SQiq2n9XgmtIEAMRqBEZD0MQJ1AVb0a7373O9b9dBqqlsAzNsIXz9x2DcYO6s2hMSKyOlGrgJPzQInI2hiAdmPqZi0+zi3E29+eRG1jMwAgJsIXqbcPwuhr/Bl4EpHNiFqtWeRKeCKyNgag3dSOX8vxyn+P41RFPQDg2mBvzJswBGMGBjDwJCKbazsHlIjImhiAdjOVdSos/vI4th0uAQAEeMoxN2Ew/i82DBIxA08iso+2c0CJiKyJAWg3IQgCth0uwaJtx3DhUhPEImDa6EjMvn0gvJQyRxePiHoYozmgjD+JyMoYgHYD1Zea8MLmX7D9WBkAYEiQF15/cCiGhvZybMGIqMcSsQeUiGyIAaiD7S+swuwNh1B8sQEyiQizbhuI5FsHQMa3FhGRA7WeA6pl/ElEVsYA1EEEQcD7u0/hja9PQKMV0M/fHe88PALXh/o4umhEREZzQLkPExFZGwNQB2hs0uCFzb/gi0O6hUb3D++LV+6LhqeCfx1E1D2IjHpAGYESkXUx4rGzczWN+MvH+3H4bDUkYhEWTboWj8X1c3SxiIiMcBU8EdkSJxra0anzdbj/vR9w+Gw1ernL8Mn0kQw+iVxQRkYGIiMjoVQqERMTgz179lwx/65duxATEwOlUon+/ftj1apVJnmWL1+OwYMHw83NDWFhYUhNTUVjY6OtqtCmB9RmjyGiHooBqJ0cLa7G/63KRUl1I/oHeOCLmaMRPyDA0cUiIivbuHEj5syZgwULFuDgwYMYM2YMJkyYgKKiIrP5CwoKMHHiRIwZMwYHDx7Eiy++iFmzZmHz5s2GPOvWrcO8efOwcOFC5OfnY82aNdi4cSPmz59v07roY1CBk0CJyMo4BG8H+wqrMO2jfahVNSO6rzfWPjESAZ4KRxeLiGxg2bJlmD59OmbMmAFA13P59ddfY+XKlUhPTzfJv2rVKoSHh2P58uUAgKioKOzfvx9vvvkmHnjgAQBAbm4uRo8ejUceeQQA0K9fPzz88MP4+eefbVoXsUgEjSBwH1Aisjr2gNpY3ukLmJr5M2pVzRgZ6Yf1T97E4JPIRanVauTl5SEhIcEoPSEhAXv37jV7TW5urkn+xMRE7N+/H01NTQCAm2++GXl5eYaA89SpU8jKysJdd93VbllUKhVqamqMDkvp54FyDigRWRt7QG3oaHE1Hv/oZ1xSazBmYAA+nBILpUzi6GIRkY1UVFRAo9EgMDDQKD0wMBBlZWVmrykrKzObv7m5GRUVFQgODsZDDz2E8+fP4+abb4YgCGhubsbTTz+NefPmtVuW9PR0LF68uEv10c0DFTgHlIisjj2gNvJ7eR2mZP6M2sZmjOznhw8eY/BJ1FO0XsAD6Pb9bZt2tfyt03fu3IlXX30VGRkZOHDgALZs2YL//ve/eOWVV9q95/z581FdXW04zpw5Y3k92pSHiMha2ANqA5V1Kkxbuw9V9WoMDfXBmsdj4SZn8Enk6gICAiCRSEx6O8vLy016OfWCgoLM5pdKpfD39wcAvPzyy3jssccM80qvv/561NfX4y9/+QsWLFgAsdi0L0GhUECh6Np0H/3bkBh/EpG1sQfUyhqbNHjqkzwUVV1CuJ87Pnr8RngpZY4uFhHZgVwuR0xMDHJycozSc3JyEB8fb/aauLg4k/zZ2dmIjY2FTKZrOy5dumQSZEokEgiCYNPeSc4BJSJbYQBqRYIgYN7mX7D/9AV4KaXIfPxG+HPBEVGPkpaWhtWrVyMzMxP5+flITU1FUVERkpOTAeiGxqdMmWLIn5ycjNOnTyMtLQ35+fnIzMzEmjVrMHfuXEOeSZMmYeXKldiwYQMKCgqQk5ODl19+Gffccw8kEtuNrrAHlIhshUPwVvTJj6fx+aESSMUirHw0Btf08XR0kYjIzpKSklBZWYklS5agtLQU0dHRyMrKQkREBACgtLTUaE/QyMhIZGVlITU1Fe+99x5CQkKwYsUKwxZMAPDSSy9BJBLhpZdeQnFxMXr37o1Jkybh1VdftW1l2ANKRDYiEpxgdnlNTQ18fHxQXV0Nb29vRxfHrCNnq/HAyr1Qa7R46a4ozBjT39FFIiILOEM70xWdqd+wxdmobmjCN2lj+Qs1EV2VJe0Mh+CtoLqhCSmf5UGt0SLh2kBMvznS0UUiIuoy/RxQJ+inICInwwDUCl7+/CjOVDUg1NcNbzw47IrbrRAROQvDHFAHl4OIXA8D0C76+lgZth0ugVgEvPPwcPi4c8U7EbkGEeeAEpGNMADtgouX1Fiw9SgA4KmxAzA83NfBJSIish79aI5W6+CCEJHLYQDaBYu/PI6KOhWu6eOJ2eMHOro4RERWxX1AichWGIB20p6T57H1YDHEIuCNB4fyNZtE5HLEnM9ORDbCALQT1M1aLNp2DAAwJa4fh96JyCXpw0/2gBKRtXUqAM3IyEBkZCSUSiViYmKwZ8+eK+Zft24dhg0bBnd3dwQHB+OJJ55AZWVlpwrcHXycW4g/ztfD30OO1DsGObo4REQ2YZgDyviTiKzM4gB048aNmDNnDhYsWICDBw9izJgxmDBhgtGbPVr7/vvvMWXKFEyfPh3Hjh3Df/7zH+zbtw8zZszocuEdoaJOheXfnAQAPH/nYPi4cdU7Ebkm/evn2QNKRNZmcQC6bNkyTJ8+HTNmzEBUVBSWL1+OsLAwrFy50mz+H3/8Ef369cOsWbMQGRmJm2++GU899RT279/f5cI7wns7fkedqhnRfb3xfzFhji4OEZHN8F3wRGQrFgWgarUaeXl5SEhIMEpPSEjA3r17zV4THx+Ps2fPIisrC4Ig4Ny5c9i0aRPuuuuuzpfaQYovNmDdj7qe3nl3RkEs5gR9InJd+haOb0IiImuzKACtqKiARqNBYGCgUXpgYCDKysrMXhMfH49169YhKSkJcrkcQUFB6NWrF9555512n6NSqVBTU2N0dAcrvjkJtUaLuP7+GH2Nv6OLQ0RkU2LOASUiG+nUIqS2r5oUBKHd108eP34cs2bNwt/+9jfk5eVh+/btKCgoQHJycrv3T09Ph4+Pj+EIC3P8UHdBRT02HTgLAJibOJiv2yQilyfiu+CJyEYsCkADAgIgkUhMejvLy8tNekX10tPTMXr0aPz1r3/F0KFDkZiYiIyMDGRmZqK0tNTsNfPnz0d1dbXhOHPmjCXFtIkPdv8BjVbAuMG9ERPBbZeIyPWxB5SIbMWiAFQulyMmJgY5OTlG6Tk5OYiPjzd7zaVLlyAWGz9GItFt2t7eb9UKhQLe3t5GhyOV1zZic14xACBl3DUOLQsRkb2wB5SIbMXiIfi0tDSsXr0amZmZyM/PR2pqKoqKigxD6vPnz8eUKVMM+SdNmoQtW7Zg5cqVOHXqFH744QfMmjULI0eOREhIiPVqYkMf/VAItUaLmAhf3NjPz9HFISKyC/aAEpGtSC29ICkpCZWVlViyZAlKS0sRHR2NrKwsREREAABKS0uN9gR9/PHHUVtbi3fffRfPPfccevXqhdtuuw1Lly61Xi1sqLaxCZ/+eBoA8NQt/R1cGiIi+9HPdRfACJSIrMviABQAUlJSkJKSYvbc2rVrTdKeffZZPPvss515lMNtyjuL2sZm9O/tgdujzM9zJSJyRfqd5tgDSkTWxnfBX4EgCFj3k64394n4ftz3k4h6FJEhAGUESkTWxQD0CvYVXsDv5XVwk0lw7/C+ji4OEZFdXX4TEgNQIrIuBqBXsO4n3dzPe28IgbeS73wnop5FxFdxEpGNMABtR1W9Gv87otvv9NFREQ4uDRGR/XEOKBHZCgPQdmw9WAy1Rovr+/rg+lAfRxeHiMju9LPeOQeUiKyNAWg7th3SbTz/f7GhDi4JEZFjiDkET0Q2wgDUjIKKehw+Ww2JWISJ1wc7ujhERA7BRUhEZCsMQM3YdqgEAHDzNQEI8FQ4uDRERI4h4hxQIrIRBqBtCIKALw7rht/vvcE5XhVKRGQL3AeUiGyFAWgbx0pqcOp8PRRSMRKuC3J0cYiIHMYwBO/gchCR62EA2sZ/fykFANweFQhPRafeVEpE5BI4B5SIbIUBaBs5x3V7fyZGs/eTiHo2DsETka0wAG3l1Pk6/HG+HlKxCLcO7u3o4hAROZT+TUharYMLQkQuhwFoK9/knwMA3NTfn6/eJKIeT/8mJPZ/EpG1MQBtJee4LgC949pAB5eEiMjx9HNAOQRPRNbGALTFxUtq5J2+AAAYH9XHwaUhImeWkZGByMhIKJVKxMTEYM+ePVfMv2vXLsTExECpVKJ///5YtWqVSZ6LFy9i5syZCA4OhlKpRFRUFLKysmxVBQCXX8XJRUhEZG0MQFvs/aMSWgG4po8nQn3dHV0cInJSGzduxJw5c7BgwQIcPHgQY8aMwYQJE1BUVGQ2f0FBASZOnIgxY8bg4MGDePHFFzFr1ixs3rzZkEetVuOOO+5AYWEhNm3ahBMnTuDDDz9E3759bVoXwxxQxp9EZGXcZ6jFnpMVAHRvPyIi6qxly5Zh+vTpmDFjBgBg+fLl+Prrr7Fy5Uqkp6eb5F+1ahXCw8OxfPlyAEBUVBT279+PN998Ew888AAAIDMzE1VVVdi7dy9kMt389IiICJvXRdoyCVTDCJSIrIw9oC2+//08AOCWQQxAiahz1Go18vLykJCQYJSekJCAvXv3mr0mNzfXJH9iYiL279+PpqYmAMC2bdsQFxeHmTNnIjAwENHR0Xjttdeg0WjaLYtKpUJNTY3RYSmJhAEoEdkGA1AApyvrcaaqATKJCKMi/R1dHCJyUhUVFdBoNAgMNF7IGBgYiLKyMrPXlJWVmc3f3NyMigrdyMypU6ewadMmaDQaZGVl4aWXXsI///lPvPrqq+2WJT09HT4+PoYjLCzM4vroe0CbGYASkZUxAMXl4ffh4b7w4NuPiKiL9HMn9QRBMEm7Wv7W6VqtFn369MEHH3yAmJgYPPTQQ1iwYAFWrlzZ7j3nz5+P6upqw3HmzBmL6yExDMFzI1Aisi5GWwB+PFUJgPM/iahrAgICIJFITHo7y8vLTXo59YKCgszml0ql8PfXjcgEBwdDJpNBIpEY8kRFRaGsrAxqtRpyudzkvgqFAgqFokv1YQ8oEdlKj+8BFQQB+wqrAAA39vNzcGmIyJnJ5XLExMQgJyfHKD0nJwfx8fFmr4mLizPJn52djdjYWMOCo9GjR+P333+HtlVP5G+//Ybg4GCzwae1SMS6fyI0GgagRGRdPT4APXuhAedqVJBJRLghrJeji0NETi4tLQ2rV69GZmYm8vPzkZqaiqKiIiQnJwPQDY1PmTLFkD85ORmnT59GWloa8vPzkZmZiTVr1mDu3LmGPE8//TQqKysxe/Zs/Pbbb/jqq6/w2muvYebMmTatC3tAichWevwQvL73M7qvD9zkkqvkJiK6sqSkJFRWVmLJkiUoLS1FdHQ0srKyDNsmlZaWGu0JGhkZiaysLKSmpuK9995DSEgIVqxYYdiCCQDCwsKQnZ2N1NRUDB06FH379sXs2bPxwgsv2LQuEm7DREQ2wgC0JQAdyeF3IrKSlJQUpKSkmD23du1ak7SxY8fiwIEDV7xnXFwcfvzxR2sUr8PYA0pEttLjh+D3FepevxnLAJSIyMjlfUC5Cp6IrKtHB6AX6tX4vbwOABAT4evg0hARdS/6HtAmLkIiIivr0QHokeJqAEBkgAf8PGy3kpSIyBkZVsFzCJ6IrIwBKIDr+/o4uCRERN0P54ASka306AD0l7MXAQBDQxmAEhG1xTchEZGt9OgA9MhZ9oASEbWHPaBEZCs9NgA9X6tCSXUjRCLgOgagREQmuA8oEdlKjw1Aj7bM/xzQ2xOeih6/HSoRkQn2gBKRrfTYAPSXluH3oez9JCIySyrhu+CJyDZ6bAB6rEQXgHL4nYjIPPaAEpGt9NgA9LdztQCAIUFeDi4JEVH3xFXwRGQrPTIAbWzS4HTVJQDAoEAGoERE5kgl7AElItvokQHo7+V1EATAz0OOAE++AYmIyBy+CYmIbKVHBqAnynTD7wP7eEIkEjm4NERE3RPngBKRrfTIAPS3cl0AOpjzP4mI2sV9QInIVnpmANrSA8r5n0RE7WMPKBHZSs8MQM/VAWAASkR0JVwFT0S20qkANCMjA5GRkVAqlYiJicGePXuumF+lUmHBggWIiIiAQqHAgAEDkJmZ2akCd1VtYxOKLzYAAAYFejqkDEREzkDasgipmRvRE5GVWfwOyo0bN2LOnDnIyMjA6NGj8f7772PChAk4fvw4wsPDzV4zefJknDt3DmvWrME111yD8vJyNDc3d7nwnfHH+XoAQG8vBXq5cwU8EVF7OAeUiGzF4gB02bJlmD59OmbMmAEAWL58Ob7++musXLkS6enpJvm3b9+OXbt24dSpU/Dz8wMA9OvXr2ul7oLTlboANDLAw2FlICJyBvp9QBmAEpG1WTQEr1arkZeXh4SEBKP0hIQE7N271+w127ZtQ2xsLF5//XX07dsXgwYNwty5c9HQ0NDuc1QqFWpqaowOayms0G1A38/f3Wr3JCJyRRIuQiIiG7GoB7SiogIajQaBgYFG6YGBgSgrKzN7zalTp/D9999DqVRi69atqKioQEpKCqqqqtqdB5qeno7FixdbUrQO0/eARvizB5SI6EqkHIInIhvp1CKktpu3C4LQ7obuWq0WIpEI69atw8iRIzFx4kQsW7YMa9eubbcXdP78+aiurjYcZ86c6UwxzSpsCUD7MQAlIrqiyz2gXAVPRNZlUQ9oQEAAJBKJSW9neXm5Sa+oXnBwMPr27QsfHx9DWlRUFARBwNmzZzFw4ECTaxQKBRQKhSVF67DTlboh+AgOwRMRXRFXwRORrVjUAyqXyxETE4OcnByj9JycHMTHx5u9ZvTo0SgpKUFdXZ0h7bfffoNYLEZoaGgnitx5NY1NqKxXA2AASkR0NXKp7p8IdTN7QInIuiwegk9LS8Pq1auRmZmJ/Px8pKamoqioCMnJyQB0w+dTpkwx5H/kkUfg7++PJ554AsePH8fu3bvx17/+FdOmTYObm5v1atIBRS29nwGecngpZXZ9NhGRs5G1rIJXaxiAEpF1WbwNU1JSEiorK7FkyRKUlpYiOjoaWVlZiIiIAACUlpaiqKjIkN/T0xM5OTl49tlnERsbC39/f0yePBl///vfrVeLDirkAiQiog4z9IBqtFec609EZCmLA1AASElJQUpKitlza9euNUkbMmSIybC9I3D+JxFRxykkEgCAIOi2YtL3iBIRdVWPehd8YQVXwBMRdZS+BxTgPFAisq4eFYAWVbEHlIiooxiAEpGt9KgAtKRat+9o3172XfxEROSMJGKRYS9QLkQiImvqMQGoRiugrLoRABDCAJSIqEMMK+HZA0pEVtRjAtCKOhWaNALEIqCPl202uSciAoCMjAxERkZCqVQiJiYGe/bsuWL+Xbt2ISYmBkqlEv3798eqVavazbthwwaIRCLcd999Vi61eXLJ5ZXwRETW0mMC0JKLuuH3IG8lpJIeU20isrONGzdizpw5WLBgAQ4ePIgxY8ZgwoQJRtvTtVZQUICJEydizJgxOHjwIF588UXMmjULmzdvNsl7+vRpzJ07F2PGjLF1NQzkUt1KePaAEpE19ZhIrOQih9+JyPaWLVuG6dOnY8aMGYiKisLy5csRFhaGlStXms2/atUqhIeHY/ny5YiKisKMGTMwbdo0vPnmm0b5NBoNHn30USxevBj9+/e3R1UAAAq+DYmIbKAHBaC6HtBgBqBEZCNqtRp5eXlISEgwSk9ISMDevXvNXpObm2uSPzExEfv370dTU5MhbcmSJejduzemT59u/YJfQevN6ImIrKVTG9E7o+KWADSkl9LBJSEiV1VRUQGNRoPAwECj9MDAQJSVlZm9pqyszGz+5uZmVFRUIDg4GD/88APWrFmDQ4cOdbgsKpUKKpXK8L2mpqbjFWnFMAeUPaBEZEU9rgeUWzARka21fWXl1V5jaS6/Pr22thZ//vOf8eGHHyIgIKDDZUhPT4ePj4/hCAsLs6AGl8mkXAVPRNbXY3pAy2p0c0CDvNkDSkS2ERAQAIlEYtLbWV5ebtLLqRcUFGQ2v1Qqhb+/P44dO4bCwkJMmjTJcF6r1QWDUqkUJ06cwIABA0zuO3/+fKSlpRm+19TUdCoI5Sp4IrKFHhOAltfohqL6MAAlIhuRy+WIiYlBTk4O7r//fkN6Tk4O7r33XrPXxMXF4csvvzRKy87ORmxsLGQyGYYMGYIjR44YnX/ppZdQW1uLt99+u92gUqFQQKHo+pZzci5CIiIb6BEBqFYroKKuJQDlHqBEZENpaWl47LHHEBsbi7i4OHzwwQcoKipCcnIyAF3PZHFxMT7++GMAQHJyMt59912kpaXhySefRG5uLtasWYP169cDAJRKJaKjo42e0atXLwAwSbcFbsNERLbQIwLQqktqNGt1c6oCPBmAEpHtJCUlobKyEkuWLEFpaSmio6ORlZWFiIgIAEBpaanRnqCRkZHIyspCamoq3nvvPYSEhGDFihV44IEHHFUFIxyCJyJb6BEBqH743c9DbhhOIiKylZSUFKSkpJg9t3btWpO0sWPH4sCBAx2+v7l72Ip+H1BVk8ZuzyQi19cjorHzHH4nIuoUpUw3BN/IIXgisqIeEYCWt6yA780AlIjIIm5y3T8TDWr2gBKR9fSMALRW3wPKFfBERJZwa+kBbeAQPBFZUY8IQM/rA1Bv9oASEVnCTa5bKsAeUCKyph4RgJbX6obgOQeUiMgy7AElIlvoGQFoyyp4zgElIrKMm4xzQInI+npEAFpVrwbAPUCJiCzlrh+CZw8oEVlRzwhAL+kCUD8PuYNLQkTkXJTyliF49oASkRW5fADarNGiuqEJAODrzgCUiMgSnANKRLbg8gFodUMTBN1bONHLXebYwhARORl39oASkQ24fAB6oWX43VsphUzi8tUlIrIqJXtAicgGXD4iq6rXDb/7cwESEZHFOARPRLbQAwJQXQ+oL4ffiYgsph+Cv6RqdnBJiMiVuHwAeoEr4ImIOs1TqduGqV6tgUYrOLg0ROQqXD4AvdwDygCUiMhSngqp4XO9mr2gRGQdLh+AXqhnDygRUWcppGLIJCIAQF0jA1Aisg6XD0D1m9D7MgAlIrKYSCSCl1I3h76O80CJyEpcPgA19IByCJ6IqFP0w/C17AElIitx+QC06pJuGyZuQk9E1Dn6AJQ9oERkLS4fgNboX8PJIXgiok7Rr4SvbWxycEmIyFX0mADUSym9Sk4iIjLHS98DyiF4IrISlw9A9XOW9JPoiYjIMl5KzgElIuty6QC0sUkDtUYLQPcueCIislyvlkWc+l1FiIi6yqUD0JqW+UoiEeAhZwBKRNQZ/i1z6KvqGIASkXW4dACqHy7yVEghFoscXBoiIuekX8TJHlAishaXDkD1C5C8Of+TiKjTDD2g9QxAicg6XDoAvbwAicPvRESd5ccAlIisrEcEoOwBJSLqPH9PXQBaWadycEmIyFV0KgDNyMhAZGQklEolYmJisGfPng5d98MPP0AqleKGG27ozGMtpl+ExB5QIqLO8/NQAABqGpvR1LKzCBFRV1gcgG7cuBFz5szBggULcPDgQYwZMwYTJkxAUVHRFa+rrq7GlClTMH78+E4X1lL6t3Z4u7EHlIios3zcZBC1rOO8wIVIRGQFFgegy5Ytw/Tp0zFjxgxERUVh+fLlCAsLw8qVK6943VNPPYVHHnkEcXFxnS6spTgHlIio6yRiEXzdOQ+UiKzHogBUrVYjLy8PCQkJRukJCQnYu3dvu9d99NFH+OOPP7Bw4cIOPUelUqGmpsbo6Ay+hpOIyDr8uBcoEVmRRQFoRUUFNBoNAgMDjdIDAwNRVlZm9pqTJ09i3rx5WLduHaTSjgWC6enp8PHxMRxhYWGWFNOAi5CIiKxDH4BWsgeUiKygU4uQRCLjTd0FQTBJAwCNRoNHHnkEixcvxqBBgzp8//nz56O6utpwnDlzpjPFRA3fA09EDmDpQs1du3YhJiYGSqUS/fv3x6pVq4zOf/jhhxgzZgx8fX3h6+uL22+/HT///LMtq2DCj0PwRGRFFgWgAQEBkEgkJr2d5eXlJr2iAFBbW4v9+/fjmWeegVQqhVQqxZIlS3D48GFIpVJ89913Zp+jUCjg7e1tdHQGV8ETkb1ZulCzoKAAEydOxJgxY3Dw4EG8+OKLmDVrFjZv3mzIs3PnTjz88MPYsWMHcnNzER4ejoSEBBQXF9urWgjw0gWg52u5FRMRdZ1FAahcLkdMTAxycnKM0nNychAfH2+S39vbG0eOHMGhQ4cMR3JyMgYPHoxDhw5h1KhRXSv9VRiG4LkKnojsxNKFmqtWrUJ4eDiWL1+OqKgozJgxA9OmTcObb75pyLNu3TqkpKTghhtuwJAhQ/Dhhx9Cq9Xi22+/tVe10LeXOwCg+GKD3Z5JRK7L4q7BtLQ0PPbYY4iNjUVcXBw++OADFBUVITk5GYBu+Ly4uBgff/wxxGIxoqOjja7v06cPlEqlSbot1Kv074KX2PxZRET6hZrz5s0zSr/SQs3c3FyThZ2JiYlYs2YNmpqaIJOZ/gJ96dIlNDU1wc/Pz3qFv4pQXzcAwNkLl+z2TCJyXRYHoElJSaisrMSSJUtQWlqK6OhoZGVlISIiAgBQWlp61T1B7aWhSQMAcJNxCJ6IbK8zCzXLysrM5m9ubkZFRQWCg4NNrpk3bx769u2L22+/vd2yqFQqqFSXh8s7u5uIXt+WALT4AntAiajrOrUIKSUlBYWFhVCpVMjLy8Mtt9xiOLd27Vrs3Lmz3WsXLVqEQ4cOdeaxFmtQtwSgcvaAEpH9dHSh5pXym0sHgNdffx3r16/Hli1boFQq272ntXYT0dP3gJbVNELdzLchEVHXuOy74AVBaNUDygCUiGzP0oWaABAUFGQ2v1Qqhb+/v1H6m2++iddeew3Z2dkYOnToFctird1E9AI8FJBLxdAKQFl1Y5fuRUTksgFok0aARqvrRWAPKBHZg6ULNQEgLi7OJH92djZiY2ON5n++8cYbeOWVV7B9+3bExsZetSzW2k1ETywWIbRXyzzQi5wHSkRd47IBqH74HWAPKBHZT1paGlavXo3MzEzk5+cjNTXVZKHmlClTDPmTk5Nx+vRppKWlIT8/H5mZmVizZg3mzp1ryPP666/jpZdeQmZmJvr164eysjKUlZWhrq7OrnXra1iIxHmgRNQ1Lrs6Rz/8LhWLIJe6bJxNRN2MpQs1IyMjkZWVhdTUVLz33nsICQnBihUr8MADDxjyZGRkQK1W48EHHzR61sKFC7Fo0SK71AsAQn11WzGdrWIPKBF1jcsHoOz9JCJ7S0lJQUpKitlza9euNUkbO3YsDhw40O79CgsLrVSyrhnQ2wMA8Ns5+/a8EpHrcdmuwUtq3R6gSs7/JCKyisFBXgCA387VOrgkROTsXDYAbWzpAXVnAEpEZBWDA3UBaGFlvaGNJSLqDJcNQC+pOQRPRGRNvb0U8HWXQSsAv5dzGJ6IOs9lA1BuQk9EZF0ikQiDAjkMT0Rd57oBKBchERFZnX4e6PGSrr3ak4h6NtcNQDkET0RkdcNCewEADhRdcGxBiMipuW4A2sQheCIia4vt5wsAOFpcw4VIRNRprh+AsgeUiMhqwv3cEeCpgFqjxZHiakcXh4iclOsGoGpuw0REZG0ikQixEbpe0P2FHIYnos5x+QCUG9ETEVnXyEg/AMCek+cdXBIiclYuG4Be4hA8EZFN3DakDwDg54Iq1DY2Obg0ROSMXDYAbeQQPBGRTfQL8ED/AA80awXsOVnh6OIQkRNy2QCUi5CIiGxnXEsv6NfHyhxcEiJyRi4bgBpexSmXOrgkRESu566hwQCA7GPnUKdqdnBpiMjZuGwAyh5QIiLbGR7WC5EBHmho0mD7UfaCEpFlXDYAbTRsRO+yVSQichiRSIQ/De8LAPjsp9MOLg0RORuXjc4MQ/AyDsETEdlC0sgwyCViHCi6iLzT3BOUiDrOZQNQwz6gMpetIhGRQ/XxUuLeG0IAABk7fndwaYjImbhsdKbWaAEACinngBIR2UryrQMgEYvw7a/l2PsHt2Qioo5x2QC0qSUAlUtFDi4JEZHrGtDbE4+OCgcAvJaVD61WcHCJiMgZuG4A2twSgErYA0pEZEuzxw+El0KKo8U1yPyhwNHFISIn4LIBqH4IXsYeUCIim/L3VGDexCEAgNe3n8CvZTUOLhERdXcuGYAKgoAmjW4YSC5xySoSEXUrj4wMx/ghfaDWaJH8SR6q6tWOLhIRdWMuGZ3pez8BQCZ1ySoSEXUrIpEISx8cir693FBYeQlPfrzfsBsJEVFbLhmd6Xs/AfaAEhHZS4CnAmufuBFeSinyTl/AlMyfUNPY5OhiEVE35JLRmbq5VQ8oA1AiIrsZGOhlCEL3FV7A/63MRUFFvaOLRUTdjEtGZ/otmCRiESRiLkIiIrKnmAg/rH/yJgR4KnDiXC0mvfM9vjxcAkHgFk1EpOOSAajasAWTS1aPiKjbi+7rg6xZN2NkPz/UqZrx7PqDmLZ2H4oqLzm6aETUDbhkhGbYgknC3k8iIkfp463EZ0+OwuzxAyGTiLDjxHmMX7YT87f8gjNVDESJejLXDED1PaB8DScRkUNJJWKk3jEI/5t9C8YMDECTRsD6n89g3Js7kbIuD9+frODbk4h6IKmjC2ALhtdwsgeUiKhbuKaPJz6ZPgr7Cqvw9jcn8f3vFcg6UoasI2Xo28sNidcFIfG6QMT28+PcfaIewCUD0Ms9oC7ZwUtE5LRu7OeHT2eMwvGSGmzYV4StB4pRfLEBmT8UIPOHAvi6yzAq0h839ffDTQP8MaiPF8QMSIlcjmsGoIY5oAxAiYi6o2tDvLHk3mi8ODEKu387j+3HyvDN8XO4cKkJ24+VYfuxMgCAh1yCa0O8cV2ID64L8cbgIC/0C/CAt1Lm4BoQUVe4ZgDKHlAicqCMjAy88cYbKC0txXXXXYfly5djzJgx7ebftWsX0tLScOzYMYSEhOD5559HcnKyUZ7Nmzfj5Zdfxh9//IEBAwbg1Vdfxf3332/rqticUiZBwnVBSLguCE0aLY4UVyP3j0r8eKoS+wsvoF6twb7CC9hXeMHougBPOSIDPBAZ4IFwP3cE+bgh2EeJQG8lgn2U8FC45D9vRC7DJf8P1b8JiT2gRGRvGzduxJw5c5CRkYHRo0fj/fffx4QJE3D8+HGEh4eb5C8oKMDEiRPx5JNP4tNPP8UPP/yAlJQU9O7dGw888AAAIDc3F0lJSXjllVdw//33Y+vWrZg8eTK+//57jBo1yt5VtBmZRIwR4b4YEe6LmeOuQbNGi1MV9ThaXI1jJTU4WlyNP87Xo6JOhYo6NSrq1CaBqZ6XUopAbyX83OXw9ZDBz0MOX3e54aevhww+bjJ4KKTwVEjhpZDBQyGBlP9uENmFSHCCnYFramrg4+OD6upqeHt7XzV/1pFSpKw7gJGRfvj3U3F2KCEROTtL25n2jBo1CiNGjMDKlSsNaVFRUbjvvvuQnp5ukv+FF17Atm3bkJ+fb0hLTk7G4cOHkZubCwBISkpCTU0N/ve//xny3HnnnfD19cX69es7VC5r1a87qG1sQmHFJZyqqENhxSUUX7yE0upGlLUctarmTt9bKRPDUyGDp0ICT6UUHnIp3OQSKKUSKGViKGUSKGUSKGTilrTW6ZfTpBIRZBIxZBIRpGKx4bNMIm51ruWzWHdOIhZBJOJ8V3JelrQzLtkDyo3oicgR1Go18vLyMG/ePKP0hIQE7N271+w1ubm5SEhIMEpLTEzEmjVr0NTUBJlMhtzcXKSmpprkWb58uVXL7yy8lDJcH+qD60N9zJ6vbWzCuZpGlNeocOFSE6ouqXGhXo2qejUuXFLr0upVqG1sRl1jM+pUzVC1/LvR2KRFY5MKFXX2rNFl8pagVCoWQSoRQywSQSIGJCIRxC1v9zN8NqTBJE0s0r0NUHe98TUSsf4zIBaLIIIIIhEgFsHwWSRq+QlALNKf1wXH+s+iVp/Rcq1YZHweIpHRfcUt94bRfds8t/UzWuJxo7BcXw7jr2h5osk1bc/D5Lzx/Yyuaefel/NZVhbje7RX3nbueZV6mCuXUZrZfGbS2uQcGelnkymNrhmAajgHlIjsr6KiAhqNBoGBgUbpgYGBKCsrM3tNWVmZ2fzNzc2oqKhAcHBwu3nauycAqFQqqFQqw/eamhpLq+O0vJQyeClluKaPV4evUTdrUa/SBaNGR2MzGpo0UDVpWoJTDRqbdZ9Vza3SDN91n5s0uqNZK6BZI0Ct0aJZo738WStAY2b/U7VGC7XGmn8aRF1z4OU74CeVW/2+nQpALZlgv2XLFqxcuRKHDh2CSqXCddddh0WLFiExMbFLBb8Sb6UUUcHeCPdzt9kziIja07b3QRCEKw6tmsvfNt3Se6anp2Px4sUdLnNPJ5eKIZfK4eth/X9o26PVCmjS6oJSXcAqoFmrRVOzLkjVCrogVaMVDJ91P2GUphEEaLVtzrdKM3w2SoMhTRAAAS0/Dd8BbavP+nStILR8b0lrOacVYHIfbavvuli71T1aPQNt7qttdV+0pOnpP+tytP6ONnkFM9fovwtmr2l7vrV2n9tOulEJrlpm4/KYy4N2y9x+Wa+aZppktu4SG00LsTgAtXSC/e7du3HHHXfgtddeQ69evfDRRx9h0qRJ+OmnnzB8+HCrVKKtO6ODcWd0sE3uTUTUnoCAAEgkEpOeyfLycpMeTL2goCCz+aVSKfz9/a+Yp717AsD8+fORlpZm+F5TU4OwsDCL6kO2JRaLoBBLwAX71BNZPEa9bNkyTJ8+HTNmzEBUVBSWL1+OsLAwown3rS1fvhzPP/88brzxRgwcOBCvvfYaBg4ciC+//LLLhSci6k7kcjliYmKQk5NjlJ6Tk4P4+Hiz18TFxZnkz87ORmxsLGQy2RXztHdPAFAoFPD29jY6iIi6C4sCUP0E+7YT5q80wb4trVaL2tpa+Pn5tZtHpVKhpqbG6CAicgZpaWlYvXo1MjMzkZ+fj9TUVBQVFRn29Zw/fz6mTJliyJ+cnIzTp08jLS0N+fn5yMzMxJo1azB37lxDntmzZyM7OxtLly7Fr7/+iqVLl+Kbb77BnDlz7F09IiKrsKjjvzMT7Nv65z//ifr6ekyePLndPJy7RETOKikpCZWVlViyZAlKS0sRHR2NrKwsREREAABKS0tRVFRkyB8ZGYmsrCykpqbivffeQ0hICFasWGHYAxQA4uPjsWHDBrz00kt4+eWXMWDAAGzcuNGl9gAlop7Fon1AS0pK0LdvX+zduxdxcZf313z11VfxySef4Ndff73i9evXr8eMGTPwxRdf4Pbbb283n7nVm2FhYS6xfx0RdU+utE+mOa5ePyJyPJvtA9qZCfZ6GzduxPTp0/Gf//znisEnoJu7pFAoLCkaERERETkJi+aAdmaCPaDr+Xz88cfx2Wef4a677upcSYmIiIjIJVi8+UNaWhoee+wxxMbGIi4uDh988IHJBPvi4mJ8/PHHAHTB55QpU/D222/jpptuMvSeurm5wcfH/FssiIiIiMh1WRyAWjrB/v3330dzczNmzpyJmTNnGtKnTp2KtWvXduiZ+mmqXA1PRLaib18smBbvVNiOEpGtWdKOWrQIyVHOnj3LDZSJyC7OnDmD0NBQRxfD6tiOEpG9dKQddYoAVKvVoqSkBF5eXld89ZyeftX8mTNnXGK1pyvVh3XpvlypPp2piyAIqK2tRUhICMRii9/R0e2xHXWd+rhSXQDXqo8r1QWwvD6WtKNO8QIwsVjcqR4JV3v7hyvVh3XpvlypPpbWxZXnpbMd1XGl+rhSXQDXqo8r1QWwrD4dbUdd79d8IiIiIurWGIASERERkV25ZACqUCiwcOFCl9nM3pXqw7p0X65UH1eqi6O42p+hK9XHleoCuFZ9XKkugG3r4xSLkIiIiIjIdbhkDygRERERdV8MQImIiIjIrhiAEhEREZFdMQAlIiIiIrtyyQA0IyMDkZGRUCqViImJwZ49exxdJBO7d+/GpEmTEBISApFIhM8//9zovCAIWLRoEUJCQuDm5oZbb70Vx44dM8qjUqnw7LPPIiAgAB4eHrjnnntw9uxZO9ZCJz09HTfeeCO8vLzQp08f3HfffThx4oRRHmepz8qVKzF06FDDprtxcXH43//+53T1MCc9PR0ikQhz5swxpDlTfRYtWgSRSGR0BAUFGc47U12cAdtR+2I72v3qYQ7bUSvWRXAxGzZsEGQymfDhhx8Kx48fF2bPni14eHgIp0+fdnTRjGRlZQkLFiwQNm/eLAAQtm7danT+H//4h+Dl5SVs3rxZOHLkiJCUlCQEBwcLNTU1hjzJyclC3759hZycHOHAgQPCuHHjhGHDhgnNzc12rUtiYqLw0UcfCUePHhUOHTok3HXXXUJ4eLhQV1fndPXZtm2b8NVXXwknTpwQTpw4Ibz44ouCTCYTjh496lT1aOvnn38W+vXrJwwdOlSYPXu2Id2Z6rNw4ULhuuuuE0pLSw1HeXm5U9alu2M7yna0K9iOdt/6dKd21OUC0JEjRwrJyclGaUOGDBHmzZvnoBJdXduGU6vVCkFBQcI//vEPQ1pjY6Pg4+MjrFq1ShAEQbh48aIgk8mEDRs2GPIUFxcLYrFY2L59u93Kbk55ebkAQNi1a5cgCM5fH19fX2H16tVOW4/a2lph4MCBQk5OjjB27FhDw+ls9Vm4cKEwbNgws+ecrS7dHdtRx/83wXa0e9WD7aj16+JSQ/BqtRp5eXlISEgwSk9ISMDevXsdVCrLFRQUoKyszKgeCoUCY8eONdQjLy8PTU1NRnlCQkIQHR3t8LpWV1cDAPz8/AA4b300Gg02bNiA+vp6xMXFOW09Zs6cibvuugu33367Uboz1ufkyZMICQlBZGQkHnroIZw6dQqAc9alu2I72j3+m2A72r3qwXbU+nWRWqEu3UZFRQU0Gg0CAwON0gMDA1FWVuagUllOX1Zz9Th9+rQhj1wuh6+vr0keR9ZVEASkpaXh5ptvRnR0NADnq8+RI0cQFxeHxsZGeHp6YuvWrbj22msN/3M5Sz0AYMOGDThw4AD27dtncs7Z/l5GjRqFjz/+GIMGDcK5c+fw97//HfHx8Th27JjT1aU7Yzvq+LqyHe0+9QDYjtqqLi4VgOqJRCKj74IgmKQ5g87Uw9F1feaZZ/DLL7/g+++/NznnLPUZPHgwDh06hIsXL2Lz5s2YOnUqdu3aZTjvLPU4c+YMZs+ejezsbCiVynbzOUt9JkyYYPh8/fXXIy4uDgMGDMC//vUv3HTTTQCcpy7OgO0o29GuYDvaPevTndpRlxqCDwgIgEQiMYnCy8vLTSL67ky/Iu1K9QgKCoJarcaFCxfazWNvzz77LLZt24YdO3YgNDTUkO5s9ZHL5bjmmmsQGxuL9PR0DBs2DG+//bbT1SMvLw/l5eWIiYmBVCqFVCrFrl27sGLFCkilUkN5nKU+bXl4eOD666/HyZMnne7vpjtjO8p21BrYjnbP+rTlyHbUpQJQuVyOmJgY5OTkGKXn5OQgPj7eQaWyXGRkJIKCgozqoVarsWvXLkM9YmJiIJPJjPKUlpbi6NGjdq+rIAh45plnsGXLFnz33XeIjIw0Ou9s9WlLEASoVCqnq8f48eNx5MgRHDp0yHDExsbi0UcfxaFDh9C/f3+nqk9bKpUK+fn5CA4Odrq/m+6M7SjbUVtgO9o96tOWQ9tRi5YsOQH99iFr1qwRjh8/LsyZM0fw8PAQCgsLHV00I7W1tcLBgweFgwcPCgCEZcuWCQcPHjRsc/KPf/xD8PHxEbZs2SIcOXJEePjhh81uhRAaGip88803woEDB4TbbrvNIds6PP3004KPj4+wc+dOo60dLl26ZMjjLPWZP3++sHv3bqGgoED45ZdfhBdffFEQi8VCdna2U9WjPa1XbwqCc9XnueeeE3bu3CmcOnVK+PHHH4W7775b8PLyMvy/7Ux16e7YjrId7Qq2o923Pt2pHXW5AFQQBOG9994TIiIiBLlcLowYMcKwjUV3smPHDgGAyTF16lRBEHTbISxcuFAICgoSFAqFcMsttwhHjhwxukdDQ4PwzDPPCH5+foKbm5tw9913C0VFRXavi7l6ABA++ugjQx5nqc+0adMM/+307t1bGD9+vKHRdKZ6tKdtw+lM9dHvRyeTyYSQkBDhT3/6k3Ds2DHDeWeqizNgO2pfbEe7Xz3aw3bUOnURCYIgWNZnSkRERETUeS41B5SIiIiIuj8GoERERERkVwxAiYiIiMiuGIASERERkV0xACUiIiIiu2IASkRERER2xQCUiIiIiOyKASgRERER2RUDUCIiIiKyKwagRERERGRXDECJiIiIyK4YgBIRERGRXf0/J9CfqqbBAhoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x300 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "yvals = np.arange(1, 500, step=1)\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(8, 3))\n",
    "ax[0].plot(yvals, force_cdf(yvals))\n",
    "ax[0].set_title('Wind force CDF')\n",
    "ax[1].plot(yvals, force_pdf(yvals))\n",
    "ax[1].set_title('Wind force PDF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "0be7842b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.00000000000000\n"
     ]
    }
   ],
   "source": [
    "from sympy import Symbol, sqrt, exp, integrate\n",
    "y = Symbol('y')\n",
    "integral_check = integrate(beta/(2*sqrt(y))*exp(-beta*sqrt(y)), (y, 0, np.inf))\n",
    "print(integral_check)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "dc7da718",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected wind force = 33.6 N\n",
      "Variance of wind force = 5644.8 N^2\n"
     ]
    }
   ],
   "source": [
    "from sympy import Symbol, sqrt, exp, integrate\n",
    "y = Symbol('y')\n",
    "expect_force = integrate(y*beta/(2*sqrt(y))*exp(-beta*sqrt(y)), (y, 0, np.inf))\n",
    "print('Expected wind force = {:.1f} N'.format(expect_force))\n",
    "\n",
    "var_force = integrate(y**2*beta/(2*sqrt(y))*exp(-beta*sqrt(y)), (y, 0, np.inf)) - expect_force**2\n",
    "print('Variance of wind force = {:.1f} N^2'.format(var_force))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92f53828",
   "metadata": {},
   "source": [
    "## Linearization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "c1e0ed39",
   "metadata": {},
   "outputs": [],
   "source": [
    "def func_h(x):\n",
    "    return 1/2 * rho * A * x ** 2\n",
    "\n",
    "def deriv_h(x, order=1):\n",
    "    if order == 1:\n",
    "        return rho * A * x\n",
    "    elif order == 2:\n",
    "        return rho * A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "24c108a7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.0 2 4.0 4\n",
      "Expected wind force by 1st order linearization = 16.8 N\n",
      "Expected wind force by 2nd order linearization = 33.6 N\n"
     ]
    }
   ],
   "source": [
    "mean_x, var_x = velocity_dist.mean(), velocity_dist.var() # note mean is actually v0 and v0**2\n",
    "print(mean_x, v0, var_x, v0**2)\n",
    "\n",
    "first_order_approx_mean = func_h(mean_x)\n",
    "print('Expected wind force by 1st order linearization = {:.1f} N'.format(first_order_approx_mean))\n",
    "\n",
    "second_order_approx_mean = func_h(mean_x) + 1/2 * var_x * deriv_h(mean_x, order=2)\n",
    "print('Expected wind force by 2nd order linearization = {:.1f} N'.format(second_order_approx_mean))\n",
    "\n",
    "# Comment on accuracy of 2nd order..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "5a847e92",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variance of wind force by 1st order linearization = 1129.0 N^2\n"
     ]
    }
   ],
   "source": [
    "first_order_approx_var = var_x * deriv_h(mean_x, order=1) ** 2\n",
    "print('Variance of wind force by 1st order linearization = {:.1f} N^2'.format(first_order_approx_var))\n",
    "\n",
    "# Better approximations would include effect of 3rd and 4th order moments"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0981f6de",
   "metadata": {},
   "source": [
    "## Monte Carlo sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "ba588157",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxQAAAEiCAYAAABgP5QIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABFwElEQVR4nO3deVxU9f4/8NewIwgKKouCoGnSNTUhDYjINMyFrmVXUq9LqMVFBUFvYmaoqYR6jdygFNTcsNzS5CqYiqZ4VcSuqVldUVwgBBMQkfXz+8PfnC/DDMsMwza8no/Heeh8zuec8/mM43zmfT7LkQkhBIiIiIiIiDSg19QFICIiIiKilosBBRERERERaYwBBRERERERaYwBBRERERERaYwBBRERERERaYwBBRERERERaYwBBRERERERaYwBBRERERERaYwBBRERERERaYwBRQPZvXs3ZDIZdu3apbSvb9++kMlkOHLkiNK+7t27o3///gCAEydOQCaT4cSJE1or182bNyGTybB58+Za8/7www9wc3ODmZkZZDIZ9u/fr7VyNIY//vgDYWFheP7552Fubg4TExP06NEDwcHB+O2336R8CxcuhEwmk7Y2bdqgS5cuGDp0KNasWYOCggKlc0+ePFnhmMrb999/r7U6yMump6eHGzduKO0vLCyEhYUFZDIZJk+erPF7oEs2b94MmUyGmzdvSmmTJ0+Gk5NTk5WJqDZsM5qOvI7yTU9PD9bW1hg+fDhSUlIU8lbOp6+vj/bt26Nv37744IMPcPbs2VrPXXlzc3PTaj3k51XVFgDA4sWLpTyVvx/lDh48CF9fX9jY2MDIyAhWVlYYPHgwtm/fjtLSUq2WtblwcnJSeL/U+byTIoOmLoCuevXVVyGTyXD8+HH4+flJ6Q8ePMDly5dhZmaG48ePY+jQodK+O3fu4MaNGwgNDQUA9O/fHykpKXjuuecavfxCCIwZMwY9e/bEgQMHYGZmhmeffbbRy6Gpc+fOYeTIkRBCYMaMGXB3d4eRkRGuX7+Obdu2YcCAAfjzzz8Vjjl8+DAsLS1RUlKCe/fu4YcffsCHH36IFStW4ODBg+jbt69CflNTUxw7dkzp2r169dJ6fczNzbFp0yZ8+umnCunffvstSktLYWhoqHSMJu+BrlqwYAGCg4ObuhhE1WKb0fRmzpyJcePGoby8HFeuXMGiRYswaNAgpKSk4IUXXpDyvfPOO5g9ezaEEMjPz8fPP/+Mr7/+Gl999RWCgoLwxRdfVHvuyszNzbVeh7Zt2+Lbb7/FmjVr0LZtWyldCIHNmzfDwsIC+fn5CscIIeDv74/Nmzdj+PDhWLVqFRwcHJCXl4fjx48jMDAQOTk5reI71M7ODikpKejevXtTF6XlEdRgnn/+efHss88qpO3du1cYGhqKoKAgMWDAAIV9X3/9tQAgDh482GBlSk9PFwDEpk2basx3584dAUBERkZq7dqPHz8WFRUVWjtfdfLy8oStra1wcHAQt2/fVpnn22+/lf4eHh4uAIj79+8r5bt06ZKwtLQUjo6O4smTJ1L6pEmThJmZmfYLX4W8bFOnThUODg6ivLxcYf/LL78sxo4dK8zMzMSkSZOkdHXfA12yadMmAUCkp6c3dVGI1MI2Q1FjtRnyOq5YsUIh/YcffpC+f+UAiOnTpyudo6ysTPj7+wsAYv369bWeuyEAEH//+9+Fqamp+OqrrxT2HT16VAAQ06ZNU/p+jIyMFADEokWLVJ43MzNTnDp1qiGL3mS6du2q0HaS5jjkqQENGjQI169fR2ZmppR24sQJvPjiixg+fDhSU1MVhtOcOHEC+vr68PLykl5X7b6ePHkyzM3N8fvvv2P48OEwNzeHg4MDZs+ejeLiYoXr37t3D2PGjEHbtm1haWkJPz8/ZGVl1VruhQsXokuXLgCAuXPnQiaTKQwX+fHHHzF48GC0bdsWbdq0gYeHBw4dOqRwDvmwk8TERPj7+6Njx45o06aNVMYdO3bA3d0d5ubmMDc3R79+/RAbG6twjqNHj2Lw4MGwsLBAmzZt4OnpiR9++KHW8m/YsAFZWVlYvny5VI+q3nnnnVrPAzwdajB//nxkZGSoHIrQWPz9/XH79m0kJSVJab/++it+/PFH+Pv7K+XXxnvw+PFjzJkzB87OzjAxMYGVlRXc3Nywc+dOKc+FCxfw7rvvwsnJCaampnBycsLYsWNx69YthXPJPw/Hjh3DtGnTYG1tDQsLC0ycOBGFhYXIysrCmDFj0K5dO9jZ2WHOnDkKXezybujly5dj6dKlcHR0hImJCdzc3Or0mVA15Ekmk2HGjBnYunUrXFxc0KZNG/Tt21flkLXvvvsOffr0gbGxMbp164YvvvhCGo5GpC1sM5qmzajOSy+9BABK32eq6OvrY+3atejQoQNWrFih8TXry9LSEm+99Rbi4uIU0uPi4uDp6YmePXsqpJeWliIyMhK9evXCggULVJ7T1tYWL7/8co3XPXbsGF599VVYW1vD1NQUjo6OGD16NB4/fizlWbRoEQYOHAgrKytYWFigf//+iI2NhRBC4VxOTk4YOXIkvv/+e7zwwgswNTWFi4uL9N28efNmuLi4wMzMDAMGDMCFCxcUjpd/5q9cuYLBgwfDzMwMHTt2xIwZMxTKo4qqIU/y7/orV65g7NixsLS0hI2NDfz9/ZGXl6dw/MOHDzFlyhRYWVnB3NwcI0aMwI0bNyCTybBw4cIar93SMaBoQIMGDQIAhS/348ePw9vbG56enpDJZDh16pTCvv79+8PS0rLG85aWluLNN9/E4MGD8d1338Hf3x+ff/45IiMjpTxFRUUYMmQIEhMTERERgW+//Ra2trYKXenVmTp1Kvbu3QvgaTdtSkoK9u3bBwBITk7Ga6+9hry8PMTGxmLnzp1o27YtfH19Vf7g9vf3h6GhIbZu3Yrdu3fD0NAQn3zyCcaPHw97e3ts3rwZ+/btw6RJkxS+tLdt2wYfHx9YWFhgy5Yt+Oabb2BlZYWhQ4fW2kAkJiZCX18fvr6+tda1Lt58800AwMmTJ5X2lZWVKWzl5eVauWZVPXr0gJeXl0IjERcXBycnJwwePFgpvzbeg9DQUERHRyMoKAiHDx/G1q1b8be//Q25ublSnps3b+LZZ59FVFQUjhw5gsjISGRmZuLFF19ETk6O0jmnTp0KS0tLxMfH4+OPP8aOHTswbdo0jBgxAn379sXu3bsxadIk/Otf/8KaNWuUjl+7di0OHz6MqKgobNu2DXp6ehg2bJjSOOe6OnToENauXYvFixdjz549sLKywltvvaUwX+Xw4cN4++23YW1tjV27dmH58uXYuXMntmzZotE1iarDNqNp2ozq/P777wCAjh071im/qakphgwZgvT0dNy5c0dhX0VFhVJ7UfWHtLZMmTIFZ8+exbVr1wA8/ZG7d+9eTJkyRSnvhQsX8ODBA/z1r3/V+AbJzZs3MWLECBgZGSEuLg6HDx/GZ599BjMzM5SUlCjk++CDD/DNN99g7969ePvttzFz5kylobwA8NNPP2HevHmYO3cu9u7dC0tLS7z99tsIDw/Hxo0bsWzZMmzfvh15eXkYOXIkioqKFI4vLS3F8OHDMXjwYOzfvx8zZszAl19+WafPc3VGjx6Nnj17Ys+ePQgLC8OOHTsQEhIi7a+oqICvry927NiBuXPnYt++fRg4cCDeeOMNja/ZojR1F4kue/DggdDT0xPvv/++EEKInJwcIZPJxOHDh4UQQgwYMEDMmTNHCCFERkaGACA+/PBD6fjjx48LAOL48eNS2qRJkwQA8c033yhca/jw4Qpd5dHR0QKA+O677xTyybs7a+u+rq6b9qWXXhKdOnUSBQUFUlpZWZno3bu36NKli9Q9LR92MnHiRIXjb9y4IfT19cX48eOrvXZhYaGwsrISvr6+Cunl5eWib9++St3+VfXq1UvY2trWmKeymoY8CSFEUVGRACCGDRsmpcn/Hapunp6edb6uumXbtGmTMDY2Frm5uaKsrEzY2dmJhQsXCiGE0pAndd8DVXr37i1GjRql1jFlZWXi0aNHwszMTHzxxRdSuvzzMHPmTIX8o0aNEgDEqlWrFNL79esn+vfvL72Wfx7t7e1FUVGRlJ6fny+srKzEkCFDlK5VuUt/0qRJomvXrgrXACBsbGxEfn6+lJaVlSX09PRERESElPbiiy8KBwcHUVxcLKUVFBQIa2trwa9Q0ia2GU3TZsjLHhkZKUpLS8WTJ09EamqqePHFFwUAcejQISkvqhnyJDd37lwBQPznP/9ROLeqLSkpqcZyqUtetoqKCuHs7Cx9VtatWyfMzc1FQUGBWLFihcL3Y3x8vAAgYmJiNL7u7t27BQBx6dKlOh9TXl4uSktLxeLFi4W1tbXC0LauXbsKU1NTcefOHSnt0qVLAoCws7MThYWFUvr+/fsFAHHgwAEpTf6Zr9wGCSHE0qVLBQDx448/KlyrctupaoifvB1evny5wvkCAwOFiYmJVPZDhw4JACI6OlohX0REhAAgwsPD6/z+tETsoWhA8tUf5HebkpOToa+vD09PTwCAt7c3jh8/DgDSn/I7VDWRyWRKd5779OmjcLfm+PHjaNu2rXR3Xa7qpDB1FBYW4j//+Q/eeecdhclk+vr6mDBhAu7cuYPr168rHDN69GiF10lJSSgvL8f06dOrvc6ZM2fw4MEDTJo0SeFuTkVFBd544w2cP38ehYWFGtdDXaKau0impqY4f/68wla1C76qqnep1OnR+Nvf/gYjIyNs374dCQkJyMrKqnY1D20YMGAA/v3vfyMsLAwnTpxQugMEAI8ePcLcuXPxzDPPwMDAAAYGBjA3N0dhYaF0d6yykSNHKrx2cXEBAIwYMUIpXdUwg7fffhsmJibSa/mdzpMnT2rUOzRo0CCFiYs2Njbo1KmTdO3CwkJcuHABo0aNgpGRkZTP3Nxcaz1gRHJsM5q2zZg7dy4MDQ1hYmICV1dXZGRk4Msvv8Tw4cPrXOfq2ovg4GCl9mLgwIE1nkvTHg35Sk9bt25FWVkZYmNjMWbMmAaZBA4A/fr1g5GREd5//31s2bJF5YqEwNNhUUOGDIGlpSX09fWl3qfc3FxkZ2crnbNz587Sa3lb8eqrr6JNmzZK6arai/Hjxyu8ln+W5f931FX1/0afPn3w5MkTqezJyckAgDFjxijkGzt2rEbXa2kYUDSwQYMG4ddff8W9e/dw/PhxuLq6Sv+pvb29kZaWJq2kYGBgUOs4RQBo06aNwo8qADA2NsaTJ0+k17m5ubCxsVE61tbWVuO6/PnnnxBCwM7OTmmfvb29dN3Kqua9f/8+AFQ7rh94utQp8HSMv6GhocIWGRkJIQQePHhQ7fGOjo64f/++1oIO+ReVvI5yenp6cHNzU9hqW9VE3p0v31QNV6qOmZkZ/Pz8EBcXh9jYWAwZMgRdu3ZVmVcb78Hq1asxd+5c7N+/H4MGDYKVlRVGjRqlsNzsuHHjsHbtWkydOhVHjhzBuXPncP78eXTs2FFlAGJlZaXwWv4jXVV65c+znKrPr62tLUpKSvDo0SO162htba2UZmxsLJVd/plX9X9JVRpRfbHNaPw2Q07+oz81NRX/+9//kJmZiffff7/W4yqrrr3o0qWLUntR+WZGVTdv3lSqi/wHa1289957uH//PpYtW4aLFy+qHO4EPG0rACA9Pb3O566qe/fuOHr0KDp16oTp06eje/fu6N69u8JqV+fOnYOPjw+Ap3P8Tp8+jfPnz2P+/PkAoNReqNNWAFBqLwwMDJS+3+Wf5aqfubqqej5jY2MA/1f23NxcGBgYKJWxtbQVXDa2gQ0aNAirVq3CiRMncOLECYU7HfKG4OTJk9LEO23dQbC2tsa5c+eU0usywa467du3h56ensKEQbl79+4BADp06KCQXnVMpnws6p07d+Dg4KDyOvJzrFmzRpoUV1VN/0GHDh2KxMREHDx4EO+++261+erqwIEDAJ7eGamvhQsXYsaMGdLrmhoUVfz9/bFx40b897//xfbt26vNp433wMzMDIsWLcKiRYvwxx9/SL0Vvr6++OWXX5CXl4fvv/8e4eHhCAsLk44rLi6uU+OtCVWf36ysLBgZGTXI3bf27dtDJpNJP1hqKwtRfbHNaPw2Q07+o19TRUVFOHr0KLp3715jAFQX9vb2OH/+vEKaOsvwOjg4YMiQIVi0aBGeffZZeHh4qMzn5uYGKysrfPfdd4iIiNB4HoWXlxe8vLxQXl6OCxcuYM2aNZg1axZsbGzw7rvvIj4+HoaGhvj+++8VgtuGelZJWVkZcnNzFYIA+WdZ1Y0kbbC2tkZZWRkePHigEFS0lraCPRQN7JVXXoG+vj52796NK1euKPwotbS0RL9+/bBlyxbcvHmzTl3XdTVo0CAUFBRIP4blduzYofE5zczMMHDgQOzdu1fhbkJFRQW2bduGLl26KK0gUZWPjw/09fURHR1dbR5PT0+0a9cOV69eVbqjI98qDz+pasqUKbC1tcWHH36Iu3fvqswjn0BYm59++gnLli2Dk5OTUjemJpycnNTq0ajK3d0d/v7+eOutt/DWW29Vm0+b7wHwtDGePHkyxo4di+vXr+Px48eQyWQQQkh3aeQ2btzYYJPT9+7dq3AnqqCgAAcPHoSXlxf09fW1fj0zMzO4ublh//79CpMLHz16pNUHGBLJsc1Q1BhthjaUl5djxowZyM3Nxdy5c+t9PiMjI7V6NFSZPXs2fH19q129CQAMDQ0xd+5c/PLLLyonRwNAdnY2Tp8+Xadr6uvrY+DAgVi3bh0A4OLFiwCeBooGBgYK39NFRUXYunVrXaujtqo33eSfZW3cHFTF29sbAJQWG4iPj2+Q6zU37KFoYPKl0fbv3w89PT1pLKyct7c3oqKiANRtLGxdTZw4EZ9//jkmTpyIpUuXokePHkhISFD5pFV1RERE4PXXX8egQYMwZ84cGBkZYf369fj555+xc+fOWu9uODk54aOPPsKnn36KoqIiaQm2q1evIicnB4sWLYK5uTnWrFmDSZMm4cGDB3jnnXfQqVMn3L9/Hz/99BPu379fY+NiaWmJ7777DiNHjsQLL7yg8FC33377Ddu2bcNPP/2Et99+W+G41NRUWFpaorS0VHqw3datW9GpUyccPHiwwRukuqptngag+XtQ2cCBAzFy5Ej06dMH7du3x7Vr17B161a4u7tLY1hfeeUVrFixAh06dICTkxOSk5MRGxuLdu3aaau6CvT19fH6668jNDQUFRUViIyMRH5+PhYtWtQg1wOePl12xIgRGDp0KIKDg1FeXo4VK1bA3Ny8wXpiqPVim6GoMdoMdf3xxx84e/YshBAoKCiQHmz3008/ISQkBNOmTdPaterDx8dHGmZUk3/+85+4du0awsPDce7cOYwbN056sN3Jkyfx1VdfYdGiRUqfRbmYmBgcO3YMI0aMgKOjI548eSKtSDhkyBAAT+fJrVq1CuPGjcP777+P3NxcrFy5UumGlLYYGRnhX//6Fx49eoQXX3wRZ86cwZIlSzBs2LA6DRPUxBtvvAFPT0/Mnj0b+fn5cHV1RUpKCr7++msAT4dJ6zIGFI1g0KBBOH/+PF544QVYWFgo7PP29sbnn38OIyOjarskNdGmTRscO3YMwcHBCAsLg0wmg4+PD+Lj4+t1HW9vbxw7dgzh4eGYPHkyKioq0LdvXxw4cEBpwm11Fi9ejB49emDNmjUYP348DAwM0KNHDwQFBUl5/v73v8PR0RHLly/HBx98gIKCAnTq1An9+vWr00TkAQMG4PLly/j888/xzTffIDIyEuXl5XBwcMDgwYOxdu1apWPkS7sZGxvDysoKzz//PCIjI/Hee++pfWeoOdDkPajstddew4EDB/D555/j8ePH6Ny5MyZOnCiNeQWe3vEJDg7Ghx9+iLKyMnh6eiIpKUlpkrW2zJgxA0+ePEFQUBCys7Pxl7/8BYcOHaq2odOGN954A3v27MEnn3wCPz8/2NraIjAwEPfu3WvQu2vUerHNUNQYbYY6du/ejd27d0NPTw/m5ubo2rUr3N3dERMTU+2Qq+ZMJpNh06ZNeOutt/DVV19h1qxZ+PPPP9G2bVv069dPager069fPyQmJiI8PBxZWVkwNzdH7969ceDAASmgee211xAXF4fIyEj4+vqic+fOmDZtGjp16lTt/I76kA+vCgoKwpIlS2Bqaopp06Y16DNC9PT0cPDgQcyePRufffYZSkpK4OnpiW3btuGll15qsBttzYVM1HXZACKiJnLz5k04OztjxYoVmDNnTlMXB6WlpdIqJImJiU1dHCIi+v8mT56M3bt3a7RQR0PYsWMHxo8fj9OnT2v1JkBzwx4KIqJaTJkyBa+//jrs7OyQlZWFmJgYXLt2TWEVEyIiat127tyJu3fv4vnnn4eenh7Onj2LFStW4JVXXtHpYAJgQEFEVKuCggLMmTMH9+/fh6GhIfr374+EhARpfDAREVHbtm0RHx+PJUuWoLCwEHZ2dpg8eTKWLFnS1EVrcBzyREREREREGtPtKedERERERNSgGFAQEREREZHGGFAQEREREZHGWsSk7IqKCty7dw9t27bV+LHwREQE6WFc9vb2WnnQ0vr167FixQpkZmbiL3/5C6KiouDl5aUyb2ZmJmbPno3U1FT89ttvCAoKkh7SVtmePXuwYMEC/O9//0P37t2xdOnSGp8MXxXbDCIi7ahzmyFagNu3bwsA3Lhx48ZNS9vt27fr/d0cHx8vDA0NxYYNG8TVq1dFcHCwMDMzE7du3VKZPz09XQQFBYktW7aIfv36ieDgYKU8Z86cEfr6+mLZsmXi2rVrYtmyZcLAwECcPXu2zuVim8GNGzdu2t1qazNaxCpPeXl5aNeuHW7fvq301FAiIqq7/Px8ODg44OHDh7C0tKzXuQYOHIj+/fsjOjpaSnNxccGoUaMQERFR47Gvvvoq+vXrp9RD4efnh/z8fPz73/+W0t544w20b98eO3furFO52GYQEWlHXduMFjHkSd5lbWFhwcaBiEgL6jsUqKSkBKmpqQgLC1NI9/HxwZkzZzQ+b0pKCkJCQhTShg4dqnJolFxxcTGKi4ul1wUFBQDYZhARaUttbQYnZRMRkdpycnJQXl4OGxsbhXQbGxtkZWVpfN6srCy1zxkREQFLS0tpc3Bw0Pj6RESkPgYURESksap3rYQQ9e79UPec8+bNQ15enrTdvn27XtcnIiL1tIghT0RE1Lx06NAB+vr6Sj0H2dnZSj0M6rC1tVX7nMbGxjA2Ntb4mkREVD/soSAiIrUZGRnB1dUVSUlJCulJSUnw8PDQ+Lzu7u5K50xMTKzXOYmIqGGpHVCcPHkSvr6+sLe3h0wmw/79+2s9Jjk5Ga6urjAxMUG3bt0QExOjSVmJiKgZCQ0NxcaNGxEXF4dr164hJCQEGRkZCAgIAPB0KNLEiRMVjrl06RIuXbqER48e4f79+7h06RKuXr0q7Q8ODkZiYiIiIyPxyy+/IDIyEkePHsWsWbMas2pERKQGtYc8FRYWom/fvnjvvfcwevToWvOnp6dj+PDhmDZtGrZt24bTp08jMDAQHTt2rNPxRETUPPn5+SE3NxeLFy9GZmYmevfujYSEBHTt2hXA0wfZZWRkKBzzwgsvSH9PTU3Fjh070LVrV9y8eRMA4OHhgfj4eHz88cdYsGABunfvjl27dmHgwIGNVi8iIlJPvZ5DIZPJsG/fPowaNaraPHPnzsWBAwdw7do1KS0gIAA//fQTUlJS6nSd/Px8WFpaIi8vj0sAEhHVQ2v4Pm0NdSQiagx1/T5t8EnZKSkp8PHxUUgbOnQoYmNjUVpaCkNDQ6Vjqq4pnp+f39DFJKJW5u7DIvxZWNLUxdBYezMjdG5n2tTFICIiaviAoro1xcvKypCTkwM7OzulYyIiIrBo0aKGLhoRNYLm+MM9t7AEAVtTUVRa3tRF0ZipoT6OzvZmUNEAnMIO4eZnI5q6GERELUajLBurak1xVely8+bNQ2hoqPRa/thvIlJPU/+Yb84/3E0N9bHFfwCszYyauihq+z37EWbtuoQ/C0sYUBARUZNr8ICiujXFDQwMYG1trfIYrilOrZm2goDm8mO+uf5w55AhIiIi7WjwgMLd3R0HDx5USEtMTISbm5vK+RNELZU2AgFtBwHN4cc8f7gTERHpNrUDikePHuH333+XXqenp+PSpUuwsrKCo6Mj5s2bh7t37+Lrr78G8HRFp7Vr1yI0NBTTpk1DSkoKYmNjsXPnTu3VgkjL1A0OtBkIaDMI4I95IiIiamhqBxQXLlzAoEGDpNfyuQ6TJk3C5s2bldYdd3Z2RkJCAkJCQrBu3TrY29tj9erVfAYFNZrGCg60FQgwCCAiIqKWRO2A4tVXX0VNj67YvHmzUpq3tzcuXryo7qWI6qSmgKExgwMGAkRERNQaNcoqT0TaoCpwqEvAwOCAiIiIqOEwoKBmpbrehpoCh9oCBgYHRERERA2HAQU1marBQ229DdUFDgwYiIiIiJoOAwpqcOoMVaqpt4GBAxEREVHzw4CCtK5yAKHuUCUGDUREREQtCwMKqpe6DFviUCUiIiIi3cWAgtQmDyLqOmyJgQMRERGR7mJAQXVSXRDBYUtERERErRsDCqpWXYIIBg9ERERErRsDCpLUNJmaQQQRERERqcKAohWrbTUmBhFEREREVBsGFK1MTROqOZmaiIiIiNTFgKIV+fluHv4Wk1LthGoGEERERESkLgYUOkzVkCYAHMZERERERFrDgEJHVe2NAJ72SHwb4I7enS2bsGREREREpEsYUOiYuw+L8Hv2I6XeCIBDmoiIiIhI+xhQ6IjKgURRaTl7I4iIiIioUTCgaMFUrdgkn2j9TCdz9kYQERERUYNjQNECVe2NAMBAgoiIiIiaBAOKFqbyZGs+eI6IiIiImhoDihZC1WRr9kYQERERUVNjQNECVO2V4GRrIiIiImouGFA0c/JgAmCvBBERERE1PwwomrHKwQR7JYiIiIioOWJA0QxVnS/BYIKIiIiImiu9pi4AKfr5bh6G/CsZk+LOAWAwQUTN2/r16+Hs7AwTExO4urri1KlTNeZPTk6Gq6srTExM0K1bN8TExCjliYqKwrPPPgtTU1M4ODggJCQET548aagqEBFRPbGHohnhfAkiakl27dqFWbNmYf369fD09MSXX36JYcOG4erVq3B0dFTKn56ejuHDh2PatGnYtm0bTp8+jcDAQHTs2BGjR48GAGzfvh1hYWGIi4uDh4cHfv31V0yePBkA8Pnnnzdm9YiIqI4YUDQDHOJERC3RqlWrMGXKFEydOhXA056FI0eOIDo6GhEREUr5Y2Ji4OjoiKioKACAi4sLLly4gJUrV0oBRUpKCjw9PTFu3DgAgJOTE8aOHYtz5841TqWIiEhtHPLUxDjEiYhaopKSEqSmpsLHx0ch3cfHB2fOnFF5TEpKilL+oUOH4sKFCygtLQUAvPzyy0hNTZUCiBs3biAhIQEjRoxogFoQEZE2sIeiCXGIExG1VDk5OSgvL4eNjY1Cuo2NDbKyslQek5WVpTJ/WVkZcnJyYGdnh3fffRf379/Hyy+/DCEEysrK8I9//ANhYWHVlqW4uBjFxcXS6/z8/HrUjIiI1MUeiiZSdUlY754dGUwQUYsjk8kUXgshlNJqy185/cSJE1i6dCnWr1+PixcvYu/evfj+++/x6aefVnvOiIgIWFpaSpuDg4Om1SEiIg2wh6IJ8PkSRNTSdejQAfr6+kq9EdnZ2Uq9EHK2trYq8xsYGMDa2hoAsGDBAkyYMEGal/H888+jsLAQ77//PubPnw89PeX7YPPmzUNoaKj0Oj8/n0EFEVEj0qiHQt1lArdv346+ffuiTZs2sLOzw3vvvYfc3FyNCtzSMZggIl1gZGQEV1dXJCUlKaQnJSXBw8ND5THu7u5K+RMTE+Hm5gZDQ0MAwOPHj5WCBn19fQghpN6MqoyNjWFhYaGwERFR41E7oJAvEzh//nykpaXBy8sLw4YNQ0ZGhsr8P/74IyZOnIgpU6bgypUr+Pbbb3H+/Hnp7lNrwmCCiHRJaGgoNm7ciLi4OFy7dg0hISHIyMhAQEAAgKc9BxMnTpTyBwQE4NatWwgNDcW1a9cQFxeH2NhYzJkzR8rj6+uL6OhoxMfHIz09HUlJSViwYAHefPNN6OvrN3odiYiodmoPeVJ3mcCzZ8/CyckJQUFBAABnZ2d88MEHWL58eT2L3rIwmCAiXePn54fc3FwsXrwYmZmZ6N27NxISEtC1a1cAQGZmpsLNJmdnZyQkJCAkJATr1q2Dvb09Vq9eLS0ZCwAff/wxZDIZPv74Y9y9excdO3aEr68vli5d2uj1IyKiulEroJAvE1h1tY2algn08PDA/PnzkZCQgGHDhiE7Oxu7d++ucQlAXVux4+7DIgYTRKSTAgMDERgYqHLf5s2bldK8vb1x8eLFas9nYGCA8PBwhIeHa6uIRETUwNQa8qTJMoEeHh7Yvn07/Pz8YGRkBFtbW7Rr1w5r1qyp9jq6tmLH79mPUFRajpgJrgwmiIiIiEinaDQpW51lAq9evYqgoCB88sknSE1NxeHDh5Geni6NsVVl3rx5yMvLk7bbt29rUsxm4ee7eQjYmgpTQ30808m8qYtDRERERKRVag150mSZwIiICHh6euKf//wnAKBPnz4wMzODl5cXlixZAjs7O6VjjI2NYWxsrE7RmqWq8yb4nAkiIiIi0jVq9VBoskxgdUsAAqh2CUBdwEnYRERERNQaqD3kSd1lAn19fbF3715ER0fjxo0bOH36NIKCgjBgwADY29trrybNCCdhExEREVFrofayseouEzh58mQUFBRg7dq1mD17Ntq1a4fXXnsNkZGR2qtFMyOfhL3FfwCDCSIiIiLSaWoHFID6ywTOnDkTM2fO1ORSLc7dh0WchE1ERERErYZGAQWpdvdhEc6nP5B6JzgJm4iIiIh0HQMKLbn7sAhD/pWMotJy9k4QERERUavBgEJL/iwsQVFpOaL8+uFFZyv2ThARERFRq6DRg+1IWW5hCQDgmU7mDCaIiIiIqNVgQKEFlSditzczauriEBERERE1Gg550oLKy8Syd4KIiIiIWhP2UNQTl4klIiIiotaMPRT1xN4JIiIiImrN2ENRD+ydICIiIqLWjj0U9cDeCSIiIiJq7dhDoSH2ThARERERsYdCY+ydICIiIiJiD4VG2DtBRERERPQUeyg08GdhCXsniIiIiIjAHgqN5BaWAACs+VRsIiIiImrlGFCoqfJwp/YMKIiIiIioleOQJzVxMjYRERER0f9hD4UaOBmbiIiIiEgReyjUwN4JIiIiIiJF7KGoI/ZOEBEREREpYw9FHbF3goiIiIhIGXso6oC9E0REREREqrGHog7YO0FEREREpBp7KGrB3gkiIiIiouqxh6IWfxaWsHeCiIiIiKga7KGoRW5hCQDAmk/FJiIiIiJSwoCiBpWHO7VnQEFEREREpIQBRQ3kw51iJrhyuBMRkQrr16+Hs7MzTExM4OrqilOnTtWYPzk5Ga6urjAxMUG3bt0QExOjlOfhw4eYPn067OzsYGJiAhcXFyQkJDRUFYiIqJ4YUNQBhzsRESnbtWsXZs2ahfnz5yMtLQ1eXl4YNmwYMjIyVOZPT0/H8OHD4eXlhbS0NHz00UcICgrCnj17pDwlJSV4/fXXcfPmTezevRvXr1/Hhg0b0Llz58aqFhERqYmTsmsgnz9BRETKVq1ahSlTpmDq1KkAgKioKBw5cgTR0dGIiIhQyh8TEwNHR0dERUUBAFxcXHDhwgWsXLkSo0ePBgDExcXhwYMHOHPmDAwNDQEAXbt2bZwKERGRRthDUQ3OnyAiql5JSQlSU1Ph4+OjkO7j44MzZ86oPCYlJUUp/9ChQ3HhwgWUlpYCAA4cOAB3d3dMnz4dNjY26N27N5YtW4by8vKGqQgREdUbeyiqwYfZERFVLycnB+Xl5bCxsVFIt7GxQVZWlspjsrKyVOYvKytDTk4O7OzscOPGDRw7dgzjx49HQkICfvvtN0yfPh1lZWX45JNPVJ63uLgYxcXF0uv8/Px61o6IiNShUQ+FupPwiouLMX/+fHTt2hXGxsbo3r074uLiNCpwY+DD7IiI6kYmkym8FkIopdWWv3J6RUUFOnXqhK+++gqurq549913MX/+fERHR1d7zoiICFhaWkqbg4ODptUhIiINqN1DIZ+Et379enh6euLLL7/EsGHDcPXqVTg6Oqo8ZsyYMfjjjz8QGxuLZ555BtnZ2SgrK6t34RsKH2ZHRFSzDh06QF9fX6k3Ijs7W6kXQs7W1lZlfgMDA1hbWwMA7OzsYGhoCH19fSmPi4sLsrKyUFJSAiMj5SGo8+bNQ2hoqPQ6Pz+fQQURUSNSu4ei8iQ8FxcXREVFwcHBodq7R4cPH0ZycjISEhIwZMgQODk5YcCAAfDw8Kh34RsaV3ciIlLNyMgIrq6uSEpKUkhPSkqq9vvd3d1dKX9iYiLc3NykCdienp74/fffUVFRIeX59ddfYWdnpzKYAABjY2NYWFgobERE1HjUCig0mYR34MABuLm5Yfny5ejcuTN69uyJOXPmoKioqNrrFBcXIz8/X2EjIqLmJTQ0FBs3bkRcXByuXbuGkJAQZGRkICAgAMDTnoOJEydK+QMCAnDr1i2Ehobi2rVriIuLQ2xsLObMmSPl+cc//oHc3FwEBwfj119/xaFDh7Bs2TJMnz690etHRER1o9aQJ00m4d24cQM//vgjTExMsG/fPuTk5CAwMBAPHjyodh5FREQEFi1apE7RtIrLxRIR1c7Pzw+5ublYvHgxMjMz0bt3byQkJEjLvGZmZio8k8LZ2RkJCQkICQnBunXrYG9vj9WrV0tLxgKAg4MDEhMTERISgj59+qBz584IDg7G3LlzG71+RERUNxqt8qTOJLyKigrIZDJs374dlpaWAJ4Om3rnnXewbt06mJoqz1FoyvGwXC6WiKjuAgMDERgYqHLf5s2bldK8vb1x8eLFGs/p7u6Os2fPaqN4RETUCNQKKDSZhGdnZ4fOnTtLwQTwdIKdEAJ37txBjx49lI4xNjaGsbGxOkXTGk7IJiIiIiKqO7XmUGgyCc/T0xP37t3Do0ePpLRff/0Venp66NKliwZFbljy4U6ckE1EREREVDu1V3lSdxLeuHHjYG1tjffeew9Xr17FyZMn8c9//hP+/v4qhzs1JQ53IiIiIiJSj9pzKNSdhGdubo6kpCTMnDkTbm5usLa2xpgxY7BkyRLt1UJLONyJiIiIiEg9Gk3KVncSXq9evZSGSTVnHO5ERERERFQ3ag95IiIiIiIikmNAUQmfP0FEREREpB4GFP8fJ2QTEZGcU9ihpi4CEVGLodEcCl3ECdlEREREROpjD0UVnJBNRERERFR3DCiIiIiIiEhjDCiIiIiIiEhjDCj+P67wRERERESkPgYU4ApPRERERESa4ipP4ApPRERERESaYg9FJVzhiYiIiIhIPQwowPkTRERERESaavUBBedPEBERERFprtXPoeD8CSIiIiIizbX6Hgo5zp8gIiIiIlIfAwoiIiIiItIYAwoiIiIiItIYAwoiIiIiItJYqw8ouGQsEREREZHmWnVAwSVjiYiIiIjqp1UvG8slY4mIiIiI6qdV91DIcclYIiIiIiLNMKAgIiIiIiKNMaAgIiJSwSnsUFMXgYioRWBAQUREREREGmNAQUREREREGmNAQUREREREGmvVAQUfakdEREREVD+tNqDgQ+2IiIiIiOqv1QYU8ofaxUxw5UPtiIg0tH79ejg7O8PExASurq44depUjfmTk5Ph6uoKExMTdOvWDTExMdXmjY+Ph0wmw6hRo7RcaiIi0qZWG1DI8aF2RESa2bVrF2bNmoX58+cjLS0NXl5eGDZsGDIyMlTmT09Px/Dhw+Hl5YW0tDR89NFHCAoKwp49e5Ty3rp1C3PmzIGXl1dDV4OIiOqp1QYUnD9BRFQ/q1atwpQpUzB16lS4uLggKioKDg4OiI6OVpk/JiYGjo6OiIqKgouLC6ZOnQp/f3+sXLlSIV95eTnGjx+PRYsWoVu3bo1RFSIiqgeNAgp1u7jlTp8+DQMDA/Tr10+Ty2oN508QEdVPSUkJUlNT4ePjo5Du4+ODM2fOqDwmJSVFKf/QoUNx4cIFlJaWSmmLFy9Gx44dMWXKlDqVpbi4GPn5+QobERE1HrUDCnW7uOXy8vIwceJEDB48WOPCagvnTxAR1U9OTg7Ky8thY2OjkG5jY4OsrCyVx2RlZanMX1ZWhpycHABPbzzFxsZiw4YNdS5LREQELC0tpc3BwUHN2hARUX2oHVCo28Ut98EHH2DcuHFwd3fXuLDaxvkTRET1I5PJFF4LIZTSassvTy8oKMDf//53bNiwAR06dKhzGebNm4e8vDxpu337tho1ICKi+jJQJ7O8izssLEwhvaYubgDYtGkT/ve//2Hbtm1YsmRJrdcpLi5GcXGx9Jrd10REzUuHDh2gr6+v1BuRnZ2t1AshZ2trqzK/gYEBrK2tceXKFdy8eRO+vr7S/oqKCgCAgYEBrl+/ju7duyud19jYGMbGxvWtEhERaUitHgpNurh/++03hIWFYfv27TAwqFv8wu5rIqLmzcjICK6urkhKSlJIT0pKgoeHh8pj3N3dlfInJibCzc0NhoaG6NWrFy5fvoxLly5J25tvvolBgwbh0qVLbAuIiJoptXoo5OraxV1eXo5x48Zh0aJF6NmzZ53PP2/ePISGhkqv8/Pz2ZAQETUzoaGhmDBhAtzc3ODu7o6vvvoKGRkZCAgIAPD0u/zu3bv4+uuvAQABAQFYu3YtQkNDMW3aNKSkpCA2NhY7d+4EAJiYmKB3794K12jXrh0AKKUTEVHzoVZAoW4Xd0FBAS5cuIC0tDTMmDEDwNPuayEEDAwMkJiYiNdee03pOHZfExE1f35+fsjNzcXixYuRmZmJ3r17IyEhAV27dgUAZGZmKizY4ezsjISEBISEhGDdunWwt7fH6tWrMXr06KaqAhERaYFaAUXlLu633npLSk9KSsJf//pXpfwWFha4fPmyQtr69etx7Ngx7N69G87OzhoWm4iImoPAwEAEBgaq3Ld582alNG9vb1y8eLHO51d1DiIial7UHvKkThe3np6eUjd1p06dVHZrExERERFRy6N2QKFuFzcREREREekujSZlq9vFXdnChQuxcOFCTS5LRERERETNjNoPttMFuYUlTV0EIiIiIiKd0OoCirsPixCwNRWmhvpozydlExERERHVi0ZDnlqyPwtLUFRaji3+A9C5nWlTF4eIiIiIqEVrdT0UctbsnSAiIiIiqrdWG1AQEREREVH9MaAgIiIiIiKNtbqAgis8ERERERFpT6sKKLjCExERERGRdrWqVZ64whMRERERkXa1qh4KOa7wRERERESkHa0yoCAiIqoLp7BDTV0EIqJmjwEFERERERFpjAEFERFRDdhLQURUMwYURERERESkMQYURERERESkMQYURERERESkMQYURERERESkMQYURERERESkMQYURERERESkMQYURERERESkMQYURERERESksVYVUOQWljR1EYiIqAXiw+2IiKrXagKKuw+LELA1FaaG+mhvZtTUxSEiIiIi0gkGTV2AxvJnYQmKSsuxxX8AOrczberiEBERERHphFbTQyFnzd4JIiIiIiKtaXUBBRERERERaQ8DCiIiIiIi0hgDCiIiIiIi0hgDCiIiIiIi0hgDCiIiIiIi0hgDCiIi0tj69evh7OwMExMTuLq64tSpUzXmT05OhqurK0xMTNCtWzfExMQo7N+wYQO8vLzQvn17tG/fHkOGDMG5c+casgpERFRPDCiIiEgju3btwqxZszB//nykpaXBy8sLw4YNQ0ZGhsr86enpGD58OLy8vJCWloaPPvoIQUFB2LNnj5TnxIkTGDt2LI4fP46UlBQ4OjrCx8cHd+/ebaxqVYtPyyYiUk2jgEKdO1J79+7F66+/jo4dO8LCwgLu7u44cuSIxgUmIqLmYdWqVZgyZQqmTp0KFxcXREVFwcHBAdHR0Srzx8TEwNHREVFRUXBxccHUqVPh7++PlStXSnm2b9+OwMBA9OvXD7169cKGDRtQUVGBH374obGqRUREalI7oFD3jtTJkyfx+uuvIyEhAampqRg0aBB8fX2RlpZW78ITEVHTKCkpQWpqKnx8fBTSfXx8cObMGZXHpKSkKOUfOnQoLly4gNLSUpXHPH78GKWlpbCystJOwYmISOsM1D2g8h0pAIiKisKRI0cQHR2NiIgIpfxRUVEKr5ctW4bvvvsOBw8exAsvvKBZqYmIqEnl5OSgvLwcNjY2Cuk2NjbIyspSeUxWVpbK/GVlZcjJyYGdnZ3SMWFhYejcuTOGDBlSbVmKi4tRXFwsvc7Pz1enKkREVE9q9VBockeqqoqKChQUFNR4t6m4uBj5+fkKGxERNT8ymUzhtRBCKa22/KrSAWD58uXYuXMn9u7dCxMTk2rPGRERAUtLS2lzcHBQpwpERFRPagUUmtyRqupf//oXCgsLMWbMmGrzsHEgImreOnToAH19faXv/uzsbKU2Qs7W1lZlfgMDA1hbWyukr1y5EsuWLUNiYiL69OlTY1nmzZuHvLw8abt9+7YGNSIiIk1pNClb3TtScjt37sTChQuxa9cudOrUqdp8bByIiJo3IyMjuLq6IikpSSE9KSkJHh4eKo9xd3dXyp+YmAg3NzcYGhpKaStWrMCnn36Kw4cPw83NrdayGBsbw8LCQmFrKFzpiYhImVpzKDS5IyW3a9cuTJkyBd9++22NY2GBp42DsbGxOkUjIqJGFhoaigkTJsDNzQ3u7u746quvkJGRgYCAAABPbw7dvXsXX3/9NQAgICAAa9euRWhoKKZNm4aUlBTExsZi586d0jmXL1+OBQsWYMeOHXBycpLaG3Nzc5ibmzd+JYmIqFZq9VBockcKeNozMXnyZOzYsQMjRozQrKRERNSs+Pn5ISoqCosXL0a/fv1w8uRJJCQkoGvXrgCAzMxMhRUAnZ2dkZCQgBMnTqBfv3749NNPsXr1aowePVrKs379epSUlOCdd96BnZ2dtFVeWpaIiJoXtVd5UveO1M6dOzFx4kR88cUXeOmll6S7TaamprC0tNRiVYiIqLEFBgYiMDBQ5b7NmzcrpXl7e+PixYvVnu/mzZtaKhkRETUWtQMKPz8/5ObmYvHixcjMzETv3r1rvCP15ZdfoqysDNOnT8f06dOl9EmTJqlsbIiIiIiIqOVQO6AA1LsjdeLECU0uoXW5hSVNXQQiItIRTmGHcPMzDuElIgI0XOWppbn7sAgBW1NhaqiP9mZGTV0cIiIiIiKdoVEPRUvzZ2EJikrLscV/ADq3M23q4hARERER6YxW0UMhZ83eCSIiIiIirWpVAQUREVF98eF2RESKGFAQEREREZHGGFAQEREREZHGGFAQERFpiMOfiIgYUBAREWmEwQQR0VMMKIiIiIiISGMMKIiIiOqBPRVE1NoxoCAiIiIiIo0xoCAiIiIiIo0xoCAiIiIiIo0xoCAiIqonzqMgotasVQQUuYUlTV0EIiIiIiKdpPMBxd2HRQjYmgpTQ320NzNq6uIQEZGOY28FEbU2Oh9Q/FlYgqLScsRMcEXndqZNXRwiItJRDCSIqLXS+YBCzpq9E0REREREWtdqAgoiIiIiItI+BhRERERaUnXYU+XXHBJFRLqKAQUREREREWmMAQUREZGWsTeCiFoTBhRERESNiMEGEekaBhREREQNTB5EVP2TiEgXMKAgIiJqALUFDQwqiEhXMKAgIiJqQHUJHBhcEFFLxoCCiIiIiIg0xoCCiIioGVA1v4I9F0TUEjCgICIiaiJOYYdUBg0MJIioJWFAQURE1AIwyCCi5sqgqQtAREREyrjELBG1FOyhICIiaiGqzq+oHHRU3qrmJSJqSBoFFOvXr4ezszNMTEzg6uqKU6dO1Zg/OTkZrq6uMDExQbdu3RATE6NRYYmIqHlpiPZgz549eO6552BsbIznnnsO+/bta6jiNzuaLDFbn2VpGXQQkTaoHVDs2rULs2bNwvz585GWlgYvLy8MGzYMGRkZKvOnp6dj+PDh8PLyQlpaGj766CMEBQVhz5499S48ERE1nYZoD1JSUuDn54cJEybgp59+woQJEzBmzBj85z//aaxqtQh1DQSq68WofA5VE8Ory6ft8hGRblA7oFi1ahWmTJmCqVOnwsXFBVFRUXBwcEB0dLTK/DExMXB0dERUVBRcXFwwdepU+Pv7Y+XKlfUuPBERNZ2GaA+ioqLw+uuvY968eejVqxfmzZuHwYMHIyoqqpFqpRs0eUp3XdKqW5WqtvxVj9PW0rgMXIiaB7UmZZeUlCA1NRVhYWEK6T4+Pjhz5ozKY1JSUuDj46OQNnToUMTGxqK0tBSGhoZKxxQXF6O4uFh6nZeXBwDIz89Xp7gAgEcF+agofoxHBfnIz5epfTwRUXNTn+81+feoEKJeZWio9iAlJQUhISFKeWoKKLTZZgBARfFjjY7TBY4h39Ypvbp8lfepyuMY8i1+XjQUFcWPFfbXdD65nxcNBQD0Dj8iva4ofoz8/HwpTZ5e+XVNaVXPJ39d1321/b026uStKX9N56m6r2qddJG67ytVr85thlDD3bt3BQBx+vRphfSlS5eKnj17qjymR48eYunSpQppp0+fFgDEvXv3VB4THh4uAHDjxo0btwbabt++rc7Xf6O1B4aGhmL79u0KebZv3y6MjIyqLQvbDG7cuHFr2K22NkOjZWNlMsU7YkIIpbTa8qtKl5s3bx5CQ0Ol1xUVFXjw4AGsra1rvI4q+fn5cHBwwO3bt2FhYaHWsS0Z6816twast/r1FkKgoKAA9vb2WilLQ7QH6p6TbUb9sd6sd2vAejdcm6FWQNGhQwfo6+sjKytLIT07Oxs2NjYqj7G1tVWZ38DAANbW1iqPMTY2hrGxsUJau3bt1CmqEgsLi1b14ZFjvVsX1rt10bTelpaW9b52Q7UH1eWp7pwA2wxtYr1bF9a7dWnINkOtSdlGRkZwdXVFUlKSQnpSUhI8PDxUHuPu7q6UPzExEW5ubirnTxARUfPXUO1BdXmqOycRETU9tVd5Cg0NxcaNGxEXF4dr164hJCQEGRkZCAgIAPC063nixIlS/oCAANy6dQuhoaG4du0a4uLiEBsbizlz5mivFkRE1Ogaoj0IDg5GYmIiIiMj8csvvyAyMhJHjx7FrFmzGrt6RERUR2rPofDz80Nubi4WL16MzMxM9O7dGwkJCejatSsAIDMzU2ENcmdnZyQkJCAkJATr1q2Dvb09Vq9ejdGjR2uvFjUwNjZGeHi4Une4rmO9We/WgPVu2no3RHvg4eGB+Ph4fPzxx1iwYAG6d++OXbt2YeDAgY1Sp+by3jY21pv1bg1Y74art0yIeq4dSERERERErZbaQ56IiIiIiIjkGFAQEREREZHGGFAQEREREZHGGFAQEREREZHGdDqgWL9+PZydnWFiYgJXV1ecOnWqqYtULydPnoSvry/s7e0hk8mwf/9+hf1CCCxcuBD29vYwNTXFq6++iitXrijkKS4uxsyZM9GhQweYmZnhzTffxJ07dxqxFuqLiIjAiy++iLZt26JTp04YNWoUrl+/rpBHF+seHR2NPn36SA+icXd3x7///W9pvy7WuaqIiAjIZDKFJUN1td4LFy6ETCZT2GxtbaX9ulrv5oRthm58pthmsM1gm9EE9RY6Kj4+XhgaGooNGzaIq1eviuDgYGFmZiZu3brV1EXTWEJCgpg/f77Ys2ePACD27dunsP+zzz4Tbdu2FXv27BGXL18Wfn5+ws7OTuTn50t5AgICROfOnUVSUpK4ePGiGDRokOjbt68oKytr5NrU3dChQ8WmTZvEzz//LC5duiRGjBghHB0dxaNHj6Q8ulj3AwcOiEOHDonr16+L69evi48++kgYGhqKn3/+WQihm3Wu7Ny5c8LJyUn06dNHBAcHS+m6Wu/w8HDxl7/8RWRmZkpbdna2tF9X691csM3Qnc8U2wy2GWwzGr/eOhtQDBgwQAQEBCik9erVS4SFhTVRibSrauNQUVEhbG1txWeffSalPXnyRFhaWoqYmBghhBAPHz4UhoaGIj4+Xspz9+5doaenJw4fPtxoZa+v7OxsAUAkJycLIVpX3du3by82btyo83UuKCgQPXr0EElJScLb21tqHHS53uHh4aJv374q9+lyvZsLthm6+5lim8E2Qxfr3dzaDJ0c8lRSUoLU1FT4+PgopPv4+ODMmTNNVKqGlZ6ejqysLIU6Gxsbw9vbW6pzamoqSktLFfLY29ujd+/eLep9ycvLAwBYWVkBaB11Ly8vR3x8PAoLC+Hu7q7zdZ4+fTpGjBiBIUOGKKTrer1/++032Nvbw9nZGe+++y5u3LgBQPfr3dTYZjylq58pthlsM3S13s2pzVD7SdktQU5ODsrLy2FjY6OQbmNjg6ysrCYqVcOS10tVnW/duiXlMTIyQvv27ZXytJT3RQiB0NBQvPzyy+jduzcA3a775cuX4e7ujidPnsDc3Bz79u3Dc889J/1n18U6x8fH4+LFizh//rzSPl3+tx44cCC+/vpr9OzZE3/88QeWLFkCDw8PXLlyRafr3Rywzfg/uvaZYpvBNgPQzXo3tzZDJwMKOZlMpvBaCKGUpms0qXNLel9mzJiB//73v/jxxx+V9uli3Z999llcunQJDx8+xJ49ezBp0iQkJydL+3Wtzrdv30ZwcDASExNhYmJSbT5dqzcADBs2TPr7888/D3d3d3Tv3h1btmzBSy+9BEA3692csM3Qvc8U2wy2GYDu1Rtofm2GTg556tChA/T19ZUirOzsbKVoTVfIZ/bXVGdbW1uUlJTgzz//rDZPczZz5kwcOHAAx48fR5cuXaR0Xa67kZERnnnmGbi5uSEiIgJ9+/bFF198obN1Tk1NRXZ2NlxdXWFgYAADAwMkJydj9erVMDAwkMqta/VWxczMDM8//zx+++03nf33bi7YZvwfXfpMsc1gm8E2o/HqrZMBhZGREVxdXZGUlKSQnpSUBA8PjyYqVcNydnaGra2tQp1LSkqQnJws1dnV1RWGhoYKeTIzM/Hzzz836/dFCIEZM2Zg7969OHbsGJydnRX263LdqxJCoLi4WGfrPHjwYFy+fBmXLl2SNjc3N4wfPx6XLl1Ct27ddLLeqhQXF+PatWuws7PT2X/v5oJtxlO68plim/F/2GawzZBr8HqrPY27hZAvARgbGyuuXr0qZs2aJczMzMTNmzebumgaKygoEGlpaSItLU0AEKtWrRJpaWnSsoafffaZsLS0FHv37hWXL18WY8eOVblEWJcuXcTRo0fFxYsXxWuvvdbsl0b7xz/+ISwtLcWJEycUlkd7/PixlEcX6z5v3jxx8uRJkZ6eLv773/+Kjz76SOjp6YnExEQhhG7WWZXKK3YIobv1nj17tjhx4oS4ceOGOHv2rBg5cqRo27at9J2lq/VuLthm6M5nim0G2wy2GY1fb50NKIQQYt26daJr167CyMhI9O/fX1oyrqU6fvy4AKC0TZo0SQjxdJmw8PBwYWtrK4yNjcUrr7wiLl++rHCOoqIiMWPGDGFlZSVMTU3FyJEjRUZGRhPUpu5U1RmA2LRpk5RHF+vu7+8vfX47duwoBg8eLDUMQuhmnVWp2jjoar3la4QbGhoKe3t78fbbb4srV65I+3W13s0J2wzd+EyxzWCbwTaj8estE0II9fs1iIiIiIiIdHQOBRERERERNQ4GFEREREREpDEGFEREREREpDEGFEREREREpDEGFEREREREpDEGFEREREREpDEGFEREREREpDEGFEREREREpDEGFEREREREpDEGFEREREREpDEGFEREREREpDEGFEREREREpLH/B406cEY6Vd7OAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x300 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# sample from input RV\n",
    "rvs_velocity = velocity_dist.rvs(10000)\n",
    "# propagate RV through nonlinear function\n",
    "#rvs_force = func_h(rvs_velocity) # broadcasting is faster\n",
    "rvs_force = []\n",
    "for x in rvs_velocity:\n",
    "    y = func_h(x)\n",
    "    rvs_force.append(y)\n",
    "\n",
    "# Plot \n",
    "bins = np.arange(1, 500, step=1)\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(8, 3))\n",
    "ax[0].hist(rvs_force, cumulative=True, histtype='step', density=True, bins=bins)\n",
    "ax[0].set_title('Wind force CDF - MC sampling')\n",
    "ax[1].hist(rvs_force, density=True, bins=bins)\n",
    "ax[1].set_title('Wind force PDF - MC sampling')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "52f74f8a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected wind force by MC sampling with large nb of samples = 32.6 N\n"
     ]
    }
   ],
   "source": [
    "print('Expected wind force by MC sampling with large nb of samples = {:.1f} N'.format(np.mean(rvs_force)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "2073e02f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variance of wind force by MC sampling with large nb of samples = 5171.5 N^2\n"
     ]
    }
   ],
   "source": [
    "print('Variance of wind force by MC sampling with large nb of samples = {:.1f} N^2'.format(np.var(rvs_force)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c076c0a1",
   "metadata": {},
   "outputs": [],
   "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}