mooring.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geometric analysis of weighted moring system stiffness\n",
    "\n",
    "Algebraics for a post at [insignificancegalore.net](https:/insignificancegalore.net/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "import sympy.geometry as geom\n",
    "from sympy import plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Construct geometry\n",
    "\n",
    "The sliders from geogebra are replaced with symbols, and we construct the points $A$, ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "a,b,c,d, v = sympy.symbols('a b c d v')\n",
    "A = geom.Point(-a, 0)\n",
    "F, G  = [geom.Point(0, dy).rotate(v, pt=A) for dy in (b/2, -b/2)]\n",
    "C, D = [geom.Point(c, dy) for dy in (d/2, -d/2)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyse mooring line lengths\n",
    "\n",
    "This is easy here, as sympy works directly with the algebraic expressions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "direct = F.distance(C)+G.distance(D)\n",
    "crossed = F.distance(D)+G.distance(C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sqrt{\\left(a \\sin{\\left(v \\right)} - \\frac{b \\cos{\\left(v \\right)}}{2} + \\frac{d}{2}\\right)^{2} + \\left(a \\cos{\\left(v \\right)} - a + \\frac{b \\sin{\\left(v \\right)}}{2} - c\\right)^{2}} + \\sqrt{\\left(a \\sin{\\left(v \\right)} + \\frac{b \\cos{\\left(v \\right)}}{2} - \\frac{d}{2}\\right)^{2} + \\left(a \\cos{\\left(v \\right)} - a - \\frac{b \\sin{\\left(v \\right)}}{2} - c\\right)^{2}}$"
      ],
      "text/plain": [
       "sqrt((a*sin(v) - b*cos(v)/2 + d/2)**2 + (a*cos(v) - a + b*sin(v)/2 - c)**2) + sqrt((a*sin(v) + b*cos(v)/2 - d/2)**2 + (a*cos(v) - a - b*sin(v)/2 - c)**2)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "direct"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the lengths to verify our construction. For this we have to choose values for the free parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de3zN9R/A8dfZxszdmPvdyGXmNqVicp8klUso5paKSr9SEqKQS0K5D0nulzARijHlfg8hmbvFLu67b5/fHx8W2tjsnPM9l/fz8dgD55yd8/Y93+19vp/P+/P+mJRSCCGEELbGxegAhBBCiLRIghJCCGGTJEEJIYSwSZKghBBC2CRJUEIIIWySJCghhBA2SRKUEEIImyQJSgghhE1yMzoAISzFZDIVBp4FigOxwBFgr1IqxdDAhBAZYpJOEsLRmEymRsAngCdwALgC5AAqARWA5cDXSqkbhgUphHgkSVDC4ZhMpq+ASUqpc2nc5wa8ALgqpX60enBCiAyTBCUclslkclVKJRsdhxDi8UiRhHBkp00mU5DJZGpiMplMRgcjhMgcSVDCkT0BbAT6opPVZJPJVN/gmIQQGWTJIT4ZOxQ24+rVq/Tr148FCxaQnCyjfkJYQZZHLeQKSji00NBQ+vTpQ+3atYmLi2Pp0qVGhySEyCC5ghIOq1y5ctSsWZMOHTrw4osvkitXLqNDEsKZZPkKShKUcFg3btwgb968RochhLOSIT4hHjRixAiio6PTTU4hISGsWbPGylEJITJLWh0Jh1O9enVat25Njhw5qF27Nl5eXsTFxXHy5EkOHjxI06ZN+fTTT40OUwjxCDLEJxxOly5dmDdvHmPHjqVw4cKEh4fj4eFBlSpV8Pf3x8PDw+gQhXAGWR7ikyso4XD27dvH2bNnWbBgAZs3b77vvtjYWElQQtgJSVDC4bz11lsEBAQQFhaGn59f6u1KKUwmE2FhYQZGJ4TIKBniEw7r7bffZtq0aUaHIYSzkjJzIYQQNsl2y8zXrIG9ey317EIIIWzR7t3w88/meS6LXUGVKoXKkwf274ccOSzyEkLYrICAANavX290GEJYVWws1Kql/zx71oavoJa/txXXY4cZNsxSryCE+fXo0YPChQvj4+OTettHH31E5cqV8fX15eWXX+batWuPfJ7IyEhLhimETRoyBNxPHGL5/343y/NZbg6qeHF1Pt6LSld3s3m7O/XqWeZlhDCnrVu3kjt3brp27cqRI0cA+OWXX2jcuDFubm4MGDAAgDFjxjz0efz8/NgrY9zCiWzfDk2ejeNUAT+K57oO58/b7hUUQUGUiv6DcbmH0a2bvuQTwtb5+/vj6el5323NmzfHzU2vyKhXrx4XLlwwIjQhbFZMDHTrBt/kHUzxq0dh5kyzPK/lEtQLL0CvXvS5PZaCJ7bx2WcWeyUhrOa7776jZcuWad4XFBSEn58ffn5+REREWDkyIYwzeDAUOxnKGzfHw9tvQ0CAWZ7XsmXmN29CjRpERJmocOMgG7bn4emnLfVyQpjHmTNneOGFF1KH+O4aOXIke/fuZcWKFTxqB3kZ4hPO4rff4AX/G4Tl8aVgkWxw8CDorW1seIgPIE8emDuXQjdPs6xMf7p105eCQtibuXPnsmbNGhYsWPDI5CSEs7h9G3r0gFm538fz9nmYN+9ucjILy2+30aABpo8+osXZILz/WsuQIRZ/RSHMav369YwZM4bVq1eTM2dOo8MRwmZ8+ilU/TuY9rfmYBo4EHNXw1mnk0R8PNSty/W/r+Ade4RVvxfi2Wct9bJCPL5OnTqxZcsWIiMjKVKkCJ9//jmjRo0iPj6eggULArpQYvr06Q99HhniE44uNBTaPXeFUx6+5K1cHHbuhOzZ732IHbU6OnQIVbcu67K14f3iSzl4yIR8GBWOShKUcGS3bkGd2orJF1+mceWLuM6dA/esHbzDxueg7lWjBqbhw3k+Zjm1/l7K4MFWe2UhhBBm9Omn8MqtH2gWE4xr545pJSezsG6z2ORk6NWL24tWUzn+EIt+K0n9+pZ6eSGMI1dQwlGFhkLX585yIlt1cjxdC0JCwNU1rYfa0RUU6P/E4MHkdItnUY7u9OiWIlV9QghhJ27fhp7dU1iSoxvu7grmzk0vOZmFdRMUQIUKmMaPp37cRpqdmsaXX1o9AiGEEI9h4EBoffpb6sVtwfTNN1C2rEVfz5j9oJSCVq1I/HULNTnI9JBKNGhgqTCEsD4Z4hOOJjQU3n7uTw651ibb880hOBgevibQzob47jKZYPZs3PJ4sMC1C290T5KhPiGEsFG3b0Pv7oksce+KW/48uteeFRasG5OgAIoVwzR9GjXjd9Pu1Gip6hNCCBs1cCB0Oj2S6vH7MM2YDkWKWOV1jUtQAB06QKdOfO7yOaET9vO7ebYQEUIIYSahobB90h66Fv0FunSBtm2t9trGzEHdKzqaFJ/qHI0sSqdSv7P7sIcs4BV2T+aghCO4fRvq+sTy08ValPO6hcvRI5A/f0a/3U7noO7l6YnL3O/x9oymW9gQGeoTQggbMXo09D4zkAqJJ3CZ+31mkpNZGJ+gAJo1w+OVlnzAePZPCGXbNqMDEkII5xYaCttHhPA+38C770LTplaPwfghvrtu3yaleg0unkvmpfJ/8NvBPDLUJ+yWDPEJe3b7Njzrc52fz1enaHkPXA4e4DF+ITvAEN9duXLhMv8HSqpzDIz4QHbgFUIIgwwaBP3OvE8xdQmX+fMeJzmZhe0kKIBnnsE0YADtrs3ixNdr2L7d6ICEEMK5/P47nP1mFd35HtOgT+HJJw2LxXaG+O6KjyfZ70mi/rzMi+WOsPlwITw8zByZEBYmQ3zCHsXEQGOfK3z5T3caVrmC645tD+7xlBkONMR3l7s7rgvmUcglmv6n3mLYUIslUCGEEPcYMljxyeneNEzapPd4evzkZBa2l6AAfH1xGTGcdvzIpXEL2bXL6ICEEMKxbd8OURN+4CWCcR010mJ7PGWG7Q3x3ZWcTNKz/uw8mIMhpeay7nBJcuQwU2RCWJgM8Ql7EhsLLaudY/WZ6uR8piZuoenu8ZQZDjjEd5erK27z51KPnQz8uwdffC5DfUIIYQmjRqbwVWR3cuZIwW3e9xbd4ykzbDdBAXh74zbxa5rzK9fHTGPPHqMDEkIIx7JzJ1z/cjJ1b4bg9u0EKFfO6JBS2XaCAnjzTRKbtGAc/Zn58Uni440OSDiyHj16ULhwYXzuGX9ftmwZ1apVw8XFRYbthEOJi4MvOh9nDANIbNEKevY0OqT72H6CMpnINnc2brnc6bYlkC+/SDI6IuHAunXrxvr16++7zcfHhxUrVuDv729QVEJYxhefJTHsdFdccuck2xzr7PGUGbafoABKlCBb0FSeYQdJo75i/36jAxKOyt/fH09Pz/tuq1KlCk888YRBEQlhGXv2gOmrMeQr4kH2mVOhWDGjQ/oP+0hQAB07ktCmPUPVUEZ1PERCgtEBCXG/oKAg/Pz88PPzIyIiwuhwhEhXfDyMfXUfwxhGuQYl4NVXjQ4pTfaToEwmss+aSkqBggw52YWxw2UyStiW3r17s3fvXvbu3YuXl5fR4QiRrlFD4xh2uivJBQuTfcZko8NJl/0kKIBChcgxbxa+HCb7yKH88YfRAQkhhH3Zvx/yjh1MNf4kx/zZ8MCQti2xrwQF0KoVca/34kP1FRPb/U5iotEBCSGEfUhIgMntQ3lfjSe++1sQEGB0SA9lfwkKyDF1PJfrPE/fk/2YMPyW0eEIB9KpUyeefvppTpw4QcmSJZk9ezYrV66kZMmS7Nixg1atWtGiRQujwxTisXw97CYjLnQjtmh53L/9yuhwHsl2Wx09ytatpDR8jpkub/LsoWm20DZKiFTS6kjYmj/+gD0136C7mo3L77/Bs89a+iUduNXRo/j7E/f2B7yZMp2Zr6wjSZZHCSFEmhITYcYr62nLcuLf+9gaycks7DdBATnHj+B6yaoMONmTKcOjjQ5HCCFs0qShkQw51Y3EqjXwGPu50eFkmF0nKHLkIO+qeRQ2RVB0xDscP250QEIIYVuOHlGUGf02hVyi8VrwDbi7Gx1Shtl3ggJMdWoT+9FQXk1ZxPzWS0hONjoiIYSwDUlJsLjNItqq5cR9+gXUqGF0SJliv0US90pKIvKJZzGF/c2PQw/Te1hxq720EGmRIglhC6Z+eoFOo6qTXLEKhY79Zu1tNJy4SOJebm4UXPsDuV1iKTO8F3+flL2jhBDO7fgxRcXRPfFwTaDgmrk2s8dTZjhGggJMlZ8gdtgYWqSs46cXZ5KSYnREQghhjKQkWNt6Os3UL8SPHIepUkWjQ3osDpOgAPIP6kt4tSb0PvEB878IMzocIYQwxPyhJ6l3aQWXfFuQ7+O3jA7nsTlUgsLFhaI/z0G5ulHmi56c/lsqJoQQzuXv40lUHh2Ib9I+iq2dbXN7PGWGYyUowFS6FHHjp1GCi/zacjyWqwERQgjbkpICIS2/ol7KDhImTMVUsoTRIWWJwyUogELvdMS1hg+Bfw/mx6HS8lwI4RyWDTpItzNDCav7KgX7djQ6nCxzjDLztF78SgTXSvpwIbko+U/sppS3/SxOE/ZPysyFtZ09EcetKnUpmi0Kz4uHMRUqaHRIUmaeHlNhL+Inz6R6yh9sa/G5DPUJIRyWUrCzxWdUU0dImjHbFpKTWThsggIo2vtFjtbrQfuwMfw8eLvR4QghhEWsGfAb7c+O42j9NynSraXR4ZiNww7x3ZVy7QaXi9YgNtEVj+MHKVYxt9EhCScgQ3zCWi4ev8mfPu0pnus6Vc7/iktem/kdJ0N8j+KSPy/Js7+nbEoYR1p+JEN9QgiHoRQcbv4hTZJ/IX/QV7aUnMzC4RMUQMnXGrK/4Qc0OzWd0E/WGR2OEEKYRehHawg4P5P9TT6mxKv1jQ7H7Bx+iO+u5NtxnPHyI1dcNC5/HqFwZU+jQxIOTIb4hKX9cyQSF18f4vIXoeTF3bh42FylsgzxZZRrrhyY5s0jm4pj94sjjA5HCCEem0pRhLV4i/wqmpTv59licjILp0lQAOXb1uJQi0944eQEdv1vkdHhCCHEY9n13gKeufQjO1oOp+yLvkaHYzFOlaAA/Fd9wB+56vHEN324eviC0eEIG9OjRw8KFy6Mj49P6m3R0dE0a9aMihUr0qxZM65evWpghMLZRR04R+KM2WzL15JnV/Y3OhyLcroE5ZbDjeyL55FNJXCheQ9kXw5xr27durF+/fr7bhs9ejRNmjTh5MmTNGnShNGjRxsUnXB6KSmEB3SjZtJevBZPxs3d/vZ4ygynS1AAlV/wZnOrr6n+z6/8+c5Uo8MRNsTf3x9Pz/sLaIKDgwkMDAQgMDCQVatWGRGaEPzR61t8rmwmtM0EKgWUNzoci3OaKr4HJcQrDpZoRfXoLSTu3E/eJysbHZKwEWfOnOGFF17gyJEjAOTPn59r166l3l+gQIE0h/mCgoIICgoCICIigrNnz1onYOEUrm3/kxzP1mZX3uY8ExFMtuw2v42GVPE9ruzuJtwXzCZGeRDVqgskJhodkrBzvXv3Zu/evezduxcvLy+jwxGOJCGBa627cJM8FFwx0x6Sk1k4bYICqNGiGBtenkG5yL383X2k0eEIG1WkSBHCw8MBCA8Pp3DhwgZHJJzNycDhlI3ez8b2Qfg0KWJ0OFbj1AkK4JWF7Vid93XKLhjBrc27jQ5H2KAXX3yRuXPnAjB37lzatGljcETCmdzauJPyi79kQ7FutJ3/stHhWJXTzkHda8+v1yjavDouBfJS4vxuyJXL6JCEQTp16sSWLVuIjIykSJEifP7557z00kt06NCBc+fOUbp0aZYtW/afQooHSScJYRa3b3O5eE1ibyQSFfIHdRrlNTqizMjyOKQkqDumtg+h8fI+5H6pGSVXTjI6HGHnJEEJczjX5l1Orz7EsY7DeWtRQ6PDySwpkjCXbnMbsyNfS0qumkxc8AajwxFCOLmYZWspvXoyfxd4isDv7C45mYUkqDty5gTvZaM4SlXiX+sOUVFGhySEcFYRESQF9uQPqlN1xQg8PIwOyBiSoO7RoFkO1rw6H4/bkUS2fxvZPEoIYXVKEflKb9xjr7L+tfk8/ZxjNoLNCElQD+g7qxaTvL6g0OZlJMxZYHQ4QggnEx80l0K/r2K850j6znDcRrAZIQnqAblzQ835H/E7z5L8dl84d87okIQQzuL0adR777GFhjy99H9OX1AsCSoNTZq7svbVeSQlpHDj5UBpKCuEsLzkZG681JXkZAjtPpfnmjh2I9iMkASVjk9mlGNYgW/Ju38LSV9NMDocIYSDSxw1jrx//M7gfJP54JsyRodjEyRBpSNfPmg6vxsreQkGfQqHDxsdkhDCUR04gGnoEJbRjpYLu5Anj9EB2QZJUA/R8nkTIa8GEZWcn2t9P4X4eKNDEkI4mrg4bnfswfqU5vz+2nSat3CORrAZIQnqET6f6kW/Aj+Q+7f1JA8cbHQ4QggHk/TxQHL9dZBFBd/liykFjQ7HpkiCegRPT+g0pwUz6YXLxK8hJMTokIQQjuLXX3GbNJFveZfOc1uQL5/RAdkWSVAZ0KYN7G4/jpOqIgmdAyGNzeqEECJToqNJeK0bf1KFQ53G0KqV0QHZHklQGfTV1Fz0zb8Alyv/kPKWdJkQQmSBUqS88SamyAje81zAV5OdtJfRI0iCyqBCheDduX4MYxguS5fAwoVGhySEsFfz5uGyYjmD1XD6zqrFI3ZvcVqSoDKhdWs48dInbDM9S/JbfeDsWaNDEkLYm9OnSe7zDjuy+ZPwbn9edq49CDNFElQmmEwweZor7+SdR1ysQnXpCsnJRoclhLAXycmo17sQG2eiT+4f+HSIdIt4GElQmVSkCPSfUo4+yZMw/bYVxo0zOiQhhL0YN44LF028nTyJT6aVwcvL6IBsmySox9C5M1x9oSsrXNqhhgyB/fuNDkkIYev27EENHsyO8yW4/VIXOnQwOiDbJ1u+P6ZLl+DZKtHsjqtOoacrYfp5rd71UAhky3fxgFu3ULVrc/lsHM/kPMT2YwUoWtTooCxOtnw3SvHiMPQbT9onzCdx5z7o39/okIQQtuqDD+Dvv3k1YR7DJztFcjILSVBZEBgIOVs2Ylpyb5g2DX76yeiQhAV98803+Pj4UK1aNSZOnGh0OMJerFwJM2cyzmUA+Vo3pHNnowOyHzLEl0Xnz0OtqvHs5Ckq5LiI6fBh5OOR4zly5AgdO3Zk9+7dZM+enYCAAKZNm0bFihXTfLwM8QkALl1C+fryV3wZGrju4OCf2Sle3OigrEaG+IxWqhSMGu/Oi7cWknz9FnTrJhscOqBjx45Rr149cubMiZubGw0bNmTlypVGhyVsWUoKdO9O0o0Y2txawNiJTpWczEISlBn06gX1uldleP7xsGEDTJ5sdEjCzHx8fNi6dStRUVHExMTw888/c/78+fseExQUhJ+fH35+fkRERBgUqbAZM2aQ+OdJhuWdQPmWlQkMNDog+yNDfGZy9iz4VFP86vEiT938FdOePVC9utFhCTOaPXs2U6ZMIXfu3FStWhUPDw8mTEh7t2UZ4nNyBw+innqKrXlf4HX35WzfYaJUKaODsjoZ4rMVZcrA2K9MvBg5m1j3/HqxVFyc0WEJM+rZsyf79+9n69ateHp6pjv/JJxcTAx07kyMR0HaRs5g6DCnTE5mIQnKjN58E3waFeb1xO/hyBH45BOjQxJmdOXKFQDOnTvHihUr6NSpk8ERCZv04Ydw7Bgd43+gdrNC9OxpdED2y83oAByJiwvMmgXVqwewtOzHtN8RimntWmSjF8fQtm1boqKiyJYtG1OmTKFAgQJGhyRszapVMH06S8v0Z0tUU47M1D08xeOROSgLmDwZ+r8bx6WST+IZ/w8cOgTFihkdlrAimYNyQhcvgq8vEbnLUuLcDiZNz86bbxodlKFkDsoW9ekDT/nnoOW1xaTcvKVX9ErpuRCOKyUFunYlJTaOFlELadA4O717Gx2U/ZMEZQEuLjB7NhxOrspU7wnw668wfrzRYQkhLGXcOAgJ4Zty3/AXTzBrlgztmYMkKAvx9oYvv4R3j/Tm/JOvwKefwr59RoclhDC3PXtgwQLO1H6FD/7syZgxUK6c0UE5BpmDsqDkZHjtNdi1Lpq/c9XANbeH3pojd26jQxMWJnNQTuLmTahVi6SEFNqW3sN1t4KEhOhRFCFzUDbN1RW++AL+SfBkcNn5qL//hnffNTosIYS59OmDOn2agcXn8uv+gsyeLcnJnORQWlilSnqob/SOhhxuPQi+/x4WLzY6LCFEVs2bB/Pnc7D1Z4zb1YAxY6BCBaODciwyxGcFycng7w8njiZxyduf7Of+hl27ZKDagckQn4M7eRJq1SLepw7FjoXgU8OVLVvk6ukBMsRnD1xdYc4cuB3vRp98C1ClSutWSImJRocmhMishATo2BGVPTtveMwnLtGV776T5GQJckitpFIlGDkSZoeU47dnBsDOnTBkiNFhCSEya+BA2L+fza9/x7wtpRg9WlftCvOTIT4rSk6GBg3g+HE43+pNcs0PgvXroUULo0MTZiZDfA5q3Tp4/nluBfah5Kop+PoiQ3vpkyE+e3J3qC82FrpdnYjy8YEuXSA83OjQhBCPEh4OkyahfH3p8s84EhKQoT0Lk0NrZU88oYf6rtz0YE2XJXDrlk5SyclGhyaESM/dRY2hoax6dTGrNngwapQM7VmaJCgD9Oun6yMCR1fl2vBJsGkTjBljdFhCiPSMGAGbN3N1xBS6j61CgwaypNEaJEEZ4N6hvq5beqA6doTPPoNt24wOTQjxoM2b4fPPUV260GVjIH5+yIJcK5FDbJC7Q30/rTGxpNEMvSVvp04QFWV0aEKIu65c0UtCKlViwTNTWfuzidatQTZTtg6p4jNQcjI0bAhHj8JfC/bgNfw98PKC4GBphWznpIrPAaSkQMuWsHUrl1fvolI7X2rUkKq9TJAqPnvm6qqrgOLioMe0uqiOneCnn+Drr40OTQgxZgz88gtq4jcEfu1LUpIempfkZD1yqA12t1ffmjUwL/+78Mor8MknsH270aEJ4bx+/10vpH/1Vb5zfYMNG5BeewaQIT4bkJwMzz0HR47A0W3XKP5CbV3md+AAFCpkdHjiMcgQnx2LjISaNSFHDs4H76fa03mpU0cX28rVU6bIEJ8juFvVV6wYDB6XH7V0mZ6cla3ihbCulBS9uWhiImrJUnr+Ly8pKbIg1ygWu4KqVq2a8vDwsMhzW0NERAReXl5Wfc0rV+D8eV3QVyjlzj9KlICiRTP9XEbEb062GP/ly5eJjIwEwMPDg3LlymFKp5jl8OHDVK9e3ZrhmZ0tvgeZ8Vjxh4fDpUtQujQReHHuHJQurWuXrM3ej/++ffuOKqV8svQkSimLfNWpU0fZMyPiT05WqmFDpfLkUersmRSl2rdXytVVqd9+y/RzyfE3rwsXLqiyZcuqmJgYpZRS7du3V3PmzEn38Tlz5rRSZJZja+9BZmU6/k2blHJxUapTJ3U6LEXlzq1Ukyb659II9n78gb0qi3lELlptiIuLHkpISYGevUyombOgbFl49VWIiDA6PKeXlJREbGwsSUlJxMTEULx4caNDEuZy6ZJeh/jEE6RMD6JnL31lLAtyjSWH3saULw/jxsHGjTBjUV5YtgwKFID33pN+fQYqUaIE/fv3p3Tp0hQrVox8+fLRvHnz+x4TFBSEn58ffn5+JCUlGRSpyLTERP0h8PZtWL6c6fNzk5SkV3uUKWN0cE4uq5dg6X3NmDHDklePFmdk/CkpSjVtqlSuXEqFhSmlgoKUAqWGDMnwc8jxN6/o6GjVqFEjdeXKFZWQkKDatGmj5s2bl+7jS5cubcXoLMPW3oPMynD8/fvrn68FC9TJk0rlzKlUixb659BI9n78gd4qi3lEysxt1LlzUL061KoFIZsULm/01KV+a9ZAq1ZGh+d0li1bxvr165k9ezYAP/zwAzt37mTq1KlpPl7KzO3EypV67WGfPiR/O4WGDfVyjyNHoGRJo4Oze1Jm7qhKl4YJEyA0FKZMNcGUKXptRpcucPq00eE5ndKlS7Nz505iYmJQSrFp0yaqVKlidFgiK06dgm7doG5dGD+eiRN1v+ZJkyQ52Qq5grJhSsELL+hmyocOQUXXMKhTB8qV050mcuQwOkSnMnToUJYsWYKbmxu1atVi1qxZuLu7p/lYuYKycbGx8PTTeqjiwAGOxZShVi0ICNAXVdIK0yyyfhSzMj4ItAeOAimA3wP3p1q3bp2qVKmSqlChgho1alTq7WFhYerJJ59U3t7eqkOHDio+Pt6cQ6CPFBUVpZo2baq8vb1V06ZNVXR09H8eExISomrUqJH65e7urlauXKmUUiowMFCVLVs29b4DBw6YPcaLF5XKl0+pLl2USkpSSv30kx4v79UrQ/ErpZSLi0tqjK1bt0693R6O/4EDB1S9evVU1apVVfXq1dXixYtT77PG8U9LeufzXXFxcapAgQKqQoUK6sknn1SnT59Ove/LL79UFSpUUJUqVVLr16+3SrwPelT8X3/9tapSpYqqXr26aty4sTpz5kzqfemdS9b0qPjnzJmjChUqlBrnzJkzU+/7/vvvlbe3t1qcJ4/+OVq7ViUmKlW3rlIFCyr1zz/Gx//++++nxl6xYkWVL1++1Pts4fh3795deXl5qWrVqqV5f0pKinr33XcV8DfwB1Bb/ZszAoGTd74C1aNyzKMe8NBvhirAE8CW9BJUUlKSKl++vDp16pSKj49Xvr6+6ujRo0opvZZk0aJFSiml3nzzTTV16lTzHMEM+uijj1JPkFGjRqmPP/74oY+PiopSBQoUULdv31ZK6V+Qy5Yts3icS5bod2rs2Ds3DBqkFKglAQEZij9Xrlxp3m4Px//EiRPqr7/+UkopdfHiRVW0aFF19epVpZT1jv+9HnY+3zVlyhRVqFAhpZRSixYtUh06dFBKKXX06FHl6+ur4uLiVFhYmCpfvrxKSkqyufhDQkJSz/GpU6emxq9U+ueStWQk/jlz5qi+ffv+53ujoqJUuXLl1K3x4ztG1P0AACAASURBVFVi1apqXP78Kjo6Wo0YoX++li61jfjv9e2336ru3bun/tvo46+UUqGhoWrfvn3pJqi1a9eqgIAAdecKqh6wS+l84QmE3fmzwJ2/F1CWWgellDqmlDrxsMfs3r0bb29vypcvT/bs2enYsSPBwcEopQgJCaFdu3YABAYGsmrVqqyEk2nBwcEEBgZm+PWXL19Oy5YtyZkzpzXCS9W+Pbz8MgwerLfm4PPPoWlT2mzYQM/atYHMHz97Of6VKlWi4p3Nd4oXL07hwoWJMHBNWHrn872Cg4MpdKeHYrt27di0aRNKKYKDg+nYsSPu7u6UK1cOb29vdu/ebXPxN2rUKPUcr1evHhcuXLBqjA+TkfjTs2HDBt709SXXgAG4lSrFyXbtmDFjJ59/rqvM27e3cPBkPv5FixbRqVMnyweWCf7+/nh6eqZ7f3BwMF27dr17kbITyG8ymYoBLYBflVLRSqmrwK9AwMNey+JFEhcvXqRUqVKp/y5ZsiQXL14kKiqK/Pnz4+bmdt/t1nT58mWKFSsGQLFixbhy5cpDH7948eL/nCyDBg3C19eX//3vf8THx1skTpMJpk+HfPmga1dITHGFhQvZ7OqK13vvQUTEQ+OPi4vDz8+PevXqpSYBezz+u3fvJiEhgQr3tJS2xvG/V3rn84OPyZYtGwBubm7ky5ePqKioDH2vpWU2htmzZ9OyZcvUf6d1LllTRuP/8ccf8fX1pV27dpw/fx6Aa8eO0SckBEqVgoULKVK8LN9+WwdPT12DZEvxA5w9e5bTp0/TuHHj1NuMPv4Z8eD/EbgAlLjzdT6N29Pl9qgXM5lMG4G0msENUko98qOLSqMIw2QypXu7uTVt2pR//vnnP7ePHDkyU88THh7O4cOHadGiRepto0aNomjRoiQkJNC7d2/GjBnDZ599luWY73Vv/DlyNGX//omUKDGFmTNLMsHdnYAzZ/THvw0b0n2Oc+fOUbx4ccLCwmjcuDHVq1cnb968/3mcrR//Ll26MHfuXFzuLO23xvF/UEbOW6PP+YfJTAzz589n7969hIaGpt6W1rlUwYp7UGQk/tatW9OpUyfc3d2ZPn06gYGBhKxfT+t583CPi9NVEJ6ebNniT3h4YVavhoIFbSf+uxYvXky7du1wdXVNvc3o458Raf0f0UVzaf1HH1pM98gEpZRqmrGw0layZMnUTzAAFy5coHjx4hQqVIhr166RlJSEm5tb6u3mtnHjxnTvK1KkCOHh4RQrVozw8HAKFy6c7mOXLl3Kyy+/nPrJGEj99O/u7k737t0ZN26c+QK/48H4u3SBxYv7UrIkhJcowbW+fcnfrx+3+/RJN/67x7V8+fI899xzHDhwgLZt29rN8b9x4watWrVixIgR1KtXL/V2axz/B6V3Pj/4mHPnzgG6PdL169fx9PTM0PdaWkZj2LhxIyNHjiQ0NPS+SsW0ziVr/oLMSPwF78k2b7zxBgMGDID//Y+SZ84ws0kT3vD1Zfdu+O23Z/D3D6N16/I2Ff9dixcvZsoDl3ZGH/+MePD/CJQELqGvmJ574PYtD32yh01QZfSLhxRJJCYmqnLlyqmwsLDUScEjR44opZRq167dfZP0U6ZMyex8XZb079//vkn6jz76KN3HPvXUUyokJOS+2y5duqSU0lUr/fr1UwMGDLBcsHdERytVvLhSVasq9f77n+j4339fKVCLW7ZM4/HRKi4uTimlVEREhPL29k6dlLWH4x8fH68aN26sJkyY8J/7jDj+Dzuf75o8efJ9RRLt27dXSil15MiR+4okypUrZ/UiiYzEv3//flW+fPnU4pS7HnYuWUtG4r97Xiil1IoVK9QX5csrBSrm3XdV2bJl1cWL0apu3UTl6npRnT6dduWrkfErpdTx48dVmTJlVMo97Sxs4fjfdfr06XSLJNasWfNgkcRu9W+RxGl0gUSBO3/3VBas4nv5TlaMBy4DG+7cXrzlPb8s165dqypWrKjKly+vRowYkXr7qVOnVN26dVWFChVUu3btUg++tURGRqrGjRsrb29v1bhxYxUVFaWUUmrPnj2qZ8+eqY87ffq0Kl68uEp+oK1xo0aNlI+Pj6pWrZp67bXX1M2bN60S97p1+p17550Y1bhxY/VEhQpqX4ECKsXdXandu++Lf9u2bcrHx0f5+voqHx8fNWvWrNTnsYfjP2/ePOXm5nZfqf/dcnKjjn9a5/OQIUNUcHCwUkqp2NjY1DLzunXrqlOnTqV+74gRI1T58uVVpUqV1M8//2yVeDMbf5MmTVThwoX/U878sHPJluL/5JNPVNWqVZWvr696s3ZtlZw9u+4dlpioZs+erfLl+065uESq/v1t8/grpdTQoUP/84HLVo5/x44dVdGiRZWbm5sqUaKEmjVrlpo2bZqaNm2aUkp/YOzTp48CTgGH7714AXqgy8//BrqrR+QYWahrp/r21RV9I0ZA/froXUD9/CApCfbufaw9pIT5yEJdG3Dlil7Y7uamfyYKFmTzZmjcWP/8TJ5sdIAOL8sTrJKg7NTNm1Cjhv77oUOQJw9w8CA88wzUrg0hIZA9u6ExOjNJUAZLTNQlr8HBuutKzZpcvw6+vuDuDgcOQK5cRgfp8KQXn7PKkwfmzoUzZ6B//zs31qypN5Tatg369TMyPCGMoxS88w6sWAHz5umfC+D99+HCBfjhB0lO9kISlB1r0EAnp6Ag+PnnOzd27Agff6wXTs2YYWh8Qhhi0iT9Q/Hhh9C2LQCrVsH338PAgXBPIaiwcTLEZ+fi4/XUU2Sk3iKgYEH0xoYdO0JYGIwfDw0bGh2m05EhPoNs2ADPPw8vvgg//gguLly5Aj4+UKIE7NolI99WJEN8zs7dXY9iREXB22/r0Q1cXWHmTIiJ0XvdnDxpdJhCWN7x43rRevXq+ofCxQWloHdvuH5d3yTJyb5IgnIANWvq9nzLlsHixXduzJ9fb25oMuk9O6KjDY1RCIuKjobWrfUnttWrIXduQM83BQfDyJH6KkrYFxnicxBJSeDvD8eOweHD92y49vvv0KQJPPusHv64pxOGsBwZ4rOixERo0UIXB23Zovd5As6e1VV7NWroPdXu6RgkrEOG+ITm5qY/Lfr56SGNlJQ7d9Svr4f7Nm+GPn3ujAEK4SCUgnff1ef3rFmpySklBYYP13d//70kJ2sZMGAAU6dOTf23yWQaZjKZPnzc55ME5UC8vaFdO1i37oFFiF27wqBB+gd4/HjD4hPC7CZN0ldOAwboRpV3TJwIs2fru8tbr9We0+vYsSNLliy596YOwLLHfT4Z4nMwSumh+E2bYN8+qFr1zh0pKXoC+ccfdTfnNm0MjdPRyRCfFfz4o97EqX17WLQI7nS5P3xYjyQEBOjyctm+3bqqVKnCpk2bKFGiRE1gqlLq2cd9LrmCcjAmk75Qyp0bXn8dEhLu3OHiolf2+vlB586wf7+hcQqRJdu36xO8Xj09hncnOcXHw2uvQYEC+udAkpP1tWvXjuXLlwO8Cix+xMMfShKUAypaVP9wHjgAw4bdc0fOnLrCqUkTaNUKTp82KkQhHt/Jk3qdU8mS+nz28Ei9a/BgfQU1ezZ4eRkYoxPr2LEji3U5cTtgeVaeSxKUg2rTBnr2hDFjdCFfqqJFYfRo/VGzRQswcPt0e3LixAlq1qyZ+pU3b14mTpxodFjOJyICWrbUl0br1kGhQql3bd4MX38Nb72lP38JY1SrVo2bN28CXFRKhWfluWQOyoHdvKnXSKWk6Iay922iu20bNG2qFzWGhKSuGxGPlpycTIkSJdi1axdlypRJ8zEyB2UBMTG6FfmhQzob3dOz6No1XVKeI4c0grUhUmYu0pcnjy49P3dON8q8z7PP6lW9+/ZBhw56LYnIkE2bNlGhQoV0k5OwgORkPee0ezcsXPifhnp9+8KlSzB/viQnRyIJysE9+yx88gnMmaOL9+7Tpg1Mm6aHSt54Q9ZIZdDixYvp1KnTf24PCgrCz88PPz8/ImTo1HyU0hkoIgImTICXX77v7sWLdc767DN48kmDYhSW8agdDbPwJWxEfLxSnTopVaqUUhcvpvGAYcP0Fr0DB1o9NnsTHx+vChYsqP7555+HPq5OnTpWisgJfPppuufnuXNK5c+vVL16SiUmGhCbeJgs5xG5gnIC2bPD0KG643m3bvd0mbjrs890+4lRo/TKRpGudevWUbt2bYoUKWJ0KM5h3Dj48kt4803dUO8eKSn6fE5M1I1g3dyMCVFYjiQoJ/HEE3p1/a+/wjffPHCnyQRTp8JLL+mNDpcuNSRGe7Bo0aI0h/eEBXz3HXz0kZ4jnTLlP4uavv1WJ6fx43UXFeF4pIrPiSilh+/XrdNzzXe3jE8VG6tXOR44oHslSa3ufWJiYihVqhRhYWHky5fvoY+VKr4sWrFCd4ho2hR++uk/+2QcPAhPPaUrzleulAW5NirL74okKCcTEaHLcQsWhD177lvjqF27pn8pHDmiF0E2b25InPZOElQWbNyoPxzVqaMv+R8oy4uJ0Xddvw5//HHfUihhW6TMXGSOl5fuDHP0qO6v+R/588Mvv0DlyrrKb8sWK0conNrOnfrK6YknYO3aNGvG+/fXexP+8IMkJ0cnCcoJtWih10VNmgQ//5zGAzw99SfXChX0Zofbtlk9RuGEdu3SJ2eDBnrvsgIF/vOQ1av1yogPP9QX+sKxyRCfk4qL02tGLl/WvcsKF07jQf/8Aw0bQni4TlhPPWX1OO2VDPFl0u7d0KyZviTasgVKlfrPQ8LD9fB0yZL6Qsvd3fphikyRIT7xeHLk0Isbb9/WDWXT/JxStKhug1S4sP5kKx3QhSXs2aPnOh+SnO6WlN++rXfWkOTkHCRBOTEfH12iO22aLtlNU4kSOknlz68/4R4+bNUYhYPbs0efV56eur9eGskJ9KbQISG6kUTlylaOURhGEpSTe+MNvcHhxx/r6vI0lS6tfzvkzKnXpRw8aNUYhYPau/ff5LRliz7P0rB/v97VvVs3vZ5cOA+ZgxJERuo1Ublz696x6TY2DwuDRo10fe/atbrRn0iTzEE9wp49esvbvHl1ckqn8e6tW1C7tl6id+iQzmXCbsgclMi6QoV0F+iTJ3UjiXSVLw+//abnpJo31+XoQmTW5s1624wGDfTfH9IV/p134NQpfX5KcnI+kqAEoC+MPv1Ud5dZ/LBNmkuX1kmqYkVdgv7jj1aLUTiAVat0+4cyZXT7orJl033oggUwd67eJbdhQ+uFKGyHDPGJVImJ+hfB0aN6mqlcuYc8+OpVvdp/1y69v3z37laL0x7IEF8a5s6FHj3Az08vwCtYMN2HnjoFtWrpoefNm6URrJ2SIT5hPtmy6dJzgM6dH7GHYYECem1Ukyb6l85/OtAKcY+JE3WVQ6NGsGnTQ5NTQgJ06gSurvoqSpKT85IEJe5TtiwEBemFkKNHP+LBuXLpRp6vvKJbU3z5pWx6KO6nlD4v/vc/fZ6sXfuQKhxt7Fg4c0ZfmKdT2CechCQo8R+vvgqDBuntdx5ZB+HuDkuW6B1Pp07VQ33x8VaJU9i4uDjo0kUvtuvXT58nj1hhu24dDBmir+DbtrVSnMJmyRyUSFNMjG6FdOWKno8qXvwR36AUDB+ud0Zs0EBvl+DEnTydfg7qyhW9t8v27foK6pNPHrknxsWLULOmPtd27kyj076wNzIHJSwjZ05Ytky3luncGZKSHvENJpPemXfRIt1X7amn4Ngxq8QqbMyff+r3f/9+fRINHPjI5JSUpM+z2Fi9X6YkJwGSoMRDVKmi2yCFhsIXX2Twmzp21Asvb92Cp5/WhRTCefzyi37f4+Jg61Zo1y5D3/bFF/rh06frnTaEAElQ4hG6dtXFVyNG6H3kMqRePX0VVaqUXvMyfbolQxS2YupUeP55XWmzaxfUrZuhb9u0SZ9f3bvD669bNkRhX2QOSjzS7dt6PioyUs9HFSuWwW+8cUNfUUVF6X0Svv3WacZunGoOKiZGt3z46y+9/GDhQsiTJ0PfGh4OPXvqqr09e9Lcn1DYL5mDEpaXK5eeF7h5E157DZKTM/iNefPqHeYCAnTN8FNP6a1QheM4cUJfMc+Zo9fErVqV4eSUlKTXO23ZoqeqJDmJB0mCEhlSrZruTLN5sy7WyzA3N/j8c10/HB6uuwgsWGCxOC3p2rVrtGvXjsqVK1OlShV27NhhdEjGWrxYv5+XLun39/PP9eraDBo6VM9vTpumzy8hHiQJSmRYt256TmrOHL37RqYEBOjxwdq19URDr156aMiO9OvXj4CAAI4fP86hQ4eoUqWK0SEZIy4O+vTRlz++vnqfloCATD3Fzz/r6vOePSEw0EJxCrsnc1AiU27d0iM5p0/rKuKSJTP5BElJ+qPzqFH6Y/OyZXaxA92NGzeoUaMGYWFhmB5RMg0OPAd1/LjuCrF+vd4bbORI3SMrE86d0332SpWCHTucZlrSGckclLCu3Ll1z8/YWGjfXvdNyxQ3N/1Lbf16uHxZt62YMkXv6W3DwsLC8PLyonv37tSqVYtevXpx+/Zto8OynuRk+PprvZL22jUIDtY9iTKZnBISoEMH3edx2TJJTuLhJEGJTKtcWQ/z7dwJ/fs/5pM0b66H/GrU0BVgzz2nq8BsVFJSEvv37+ftt9/mwIED5MqVi9EPNCsMCgrCz88PPz8/IiIiDIrUAv76C/z99ZsdEAArV8KLLz7WU33yia5AnzNH79gixEMppSz1JRzcBx8oBUotXJiFJ0lJUeq775TKn1+pHDmUGjNGqcREs8VoLuHh4apMmTKp/966dat6/vnn0318nTp1rBCVhSUlKTV+vH5fChRQat48/X49pkWLlCpfXqlBg8wYo7BlWc4jcgUlHtvo0VC/vq53OHr0MZ/EZNIrNP/8Uy/qHTBAly0fOmTWWLOqaNGilCpVihMnTgCwadMmqlatanBUFrRvn65e+OADaNoUjhzRxS0ZmH9Ly5EjuiCiaFHdEUuIjJAiCZEl4eF6wjtfPr3QMm/eLDyZUnqH3r59IXt2PT81eDDkz2+2eLPi4MGD9OrVi4SEBMqXL8+cOXMoUKBAmo+12yKJy5d1K/vvvgMvL5g0SU82PmZiArh+XTeVuHlT571HNh4WjiLLRRIyxCeyLDRUKVdXpdq2zdII0L8iI5X69FOlTCalChZUavJkpRISzPDE1mN3Q3xxcUqNHatUnjxKubkp9eGHSl27luWnTU5W6sUX9VNu3WqGOIU9kSE+YTx/fz3c9+OPZmq7V7CgrvTbv//fIorq1WHNGtkQ0dyU0ptO+vjAxx/rN/PIERg3Tl8WZ9Ho0bqZyLhxehcWITJDEpQwiw8/1Mtj3n3XjA3Ma9bUHWpXr9a/SFu31otDt20z0ws4MaV0j6HOnXVFnpub7gaxZo3Z2on/8oseoe3cGd57zyxPKZyMJChhFiaT3jKhShU9dXTqlBmfuHVr/an+2291x9r69XVZ+q+/yhVVZimlE1GDBtCokU5Ss2fDH39kuhvEw4SF6asnHx8ICsrSFJZwYpKghNnkzq3Xb5pM0KaNnhQ3m2zZ9OVZcDBMmAAnT+q1VE89pW+z8YW+hktJ0euX/Pz0lhjnzsHkyTqT9OiR6QW3D3Prln7/Dx7ULylNYMXjkgQlzKp8eViyRG+mGxhogbyRKxe8/77+xRoUpLfyeOkl3Z5gyhRdMib+FRUFEyfq7ddfeUVvgTJ7Nvz9t66WNHMrh5QU3a/xzz91L9kKFcz69MLJSIISZte0qZ4UX7lSb0RnEe7u8MYberuH+fP1L9533tE1zL166Zp3Zx3+U0q3CX/tNShRQk8ORkb++8mhRw9dxm8Bw4fr933cOH2BK0RWyDooYRFK6SuoefP0FkFt2ljhRffuhRkz9IZ5MTHQtq0e0mrf3uof5Q1ZB3X+PKxYoXe2/esvXYXXtatO5NWrW/zlV67UF2ldu8L338u8k8j6OihJUMJiYmOhYUP9oX3nTivu+XP9uk5SS5boKwmAOnX0MGD79lCunMVDsFqCCgvTmWHZMt3kzt9fj7O98Qa0awc5c1o+BuDwYXj6af0eh4ZCjhxWeVlh2yRBCdt24YK+iClRQpcdFyxo5QDOnoXly3Wy2rNHVwBeu6YLBVq00L9VLdBS22IJKiZGZ/uff9ZfoD8B1Kqlk2+7dlbvwhoRobtUXbqkL2KlU4S4QxKUsH07dujpkNKldZKy0PTHo50+rQNYulSXmF2/rneAbdVK18fXq6erAgsXzvJLmS1BXb4Mu3frJLR6tf57rVo6/oYNdaJ94QXw9s76az2GuDg953jxov4M8OSThoQhbJMkKGEfFi7USSowUG+1YPj8xM2b8Ntveg/78HBdcpacrIcCr1zRHSz8/XVW9fbWX5norJDpBHX9uq6sO3FCJ6PISN24bs8efX+NGrow5Lnn9PqlBg0Mr99WSvePXbhQH75XXzU0HGF7svxT7maOKIR4lM6d9dKlYcN0o4KBAw0OKE8effXx/PP630FBOiEcOwZbt/57hXX7tm7VHh+vhwdLlNDNa7Nn1wu/ihYFT0+dLDw89O2urjoBhobqHYQTE/WE3K1bOvklJOj7//hDX3q4uurv37dPv17Fivoqydtbd86oWxdq17bafFJGDRumk9PIkZKchGXIFZSwmns/cS9dqqdMbFp8vL6quXtlc+XKv0nr7Fk93FaokM64KSn6Kz4ebt3C78wZ9vr66mqBpCT9ZTLpooZq1fTtly/r6sLKlaFYMf1VsaL+cnc3+n//UPPm6Wq97t31sirDr4iFLZIhPmFf7s5Z7NunR9fq1TM6oixQSq+/un5dXxHFx+uvpCT8evVi7+zZusdd9uz66ipnTihQQF+92fFv9K1b9XtYvz6sX2/gnKKwdZKghP2JiNCJ6dYtXRldtqzREZmf3e4H9Qh//aXfuyJFYPt2nW+FSEeWE5R0khBW5+UFa9fqEbEPP4ToaKMjEhkRGakLHt3c9PsnyUlYmiQoYYjKlfX60jVrdJV0TIzREYmHiY2F/v11s4rgYN1zUQhLkwQlDFO/vi6Y2LlTN3lITDQ6IpGWpCRdpffDD7q45emnjY5IOAtJUMJQbdvq1nFr1+ruPM7a39VWKaXfl59+0rtzvPii0REJZyLroITh3npLV1wPG6abOIwda3RE4q4BA3Tj12HDoE8fo6MRzkYSlLAJn32mlxl99ZWuEPvwQ6Mj+q+yZcuSJ08eXF1dcXNzc8gqvXt99ZX+6ttXvz9CWJskKGETTCa9o/uVK3oyvnBh6NLF6Kj+a/PmzRQqVMjoMCzu++/h44/13NO339r1si1hxyRBCZvh6qr3HoyO1pvmFiigK/yEda1dq/d8bNZMF0a4yEy1MIicesKmuLvr8nN/f11AcXdHCVtgMplo3rw5derUISgoyOhwLGLtWr07fOfOeu9D6RIhjCSdJIRNunpVt9M5elTvyBsQYHREcOnSJYoXL86VK1do1qwZkyZNwt/fP/X+oKCg1MQVERHB2bNnjQr1sWzYoKv0qleHjRt1T1whskBaHQnHFR0NTZr8uxVS8+ZGR/SvYcOGkTt3bvr375/m/fbW6mjjRj2cWqUKbNqkG7QLkUXS6kg4Lk9P/YuzcmVo00b/3Si3b9/m5s2bqX//5Zdf8PHxMS4gMwoJgdatdVP2jRslOQnbIQlK2LSCBfUvzYoV9fBTSIgxcVy+fJn69etTo0YNnnzySVq1akWALYw7ZtGWLTo5eXvr41ywoNERCfEvGeITdiEiAho3hlOn9ER+o0ZGR/Rw9jDE98sveuNIFxd9TM2w070Q95IhPuEcvLz03EiTJrpvX3Cw0RHZt+XL9ZxTcrJu2CvJSdgiSVDCbhQuDHPm6E1oX3kFZs0yOiL7NHu2XoBbt64e4itSxOiIhEibJChhVwoV0ldSzZvrJqZffikNZjPj66//XYT7yy9SSi5smyQoYXdy5dJl56+/DoMGQb9+evNDkT6lYORI3UaqfXt9/HLlMjoqIR5OWh0Ju5QtG8ydq4f9xo/XRRRz50rng7TExkK3bnohbv/+MHq0bislhK2TBCXslosLjBun51AGDAAPD/1vWcfzr/BwvYZs716dmD76SBq/CvshQ3zCrplMuuv2kiV6t1c/Pzh0yOiobMOBA/Dkk7pd1IoV+jhJchL2RBKUcAgdOuiFpgkJekvyhQuNjshYq1ZB/fr679u2wUsvGRuPEI9DEpRwGPXqwb59+irqtdfgf/+DxESjo7Ku5GQ9lNehA/j4wJ49ULOm0VEJ8XgkQQmHUqSILkN/7z2YOFGXU1+5YnRU1nHunO6wMXCgLobYsgWKFjU6KiEenyQo4XCyZYNvvoF582DXLr1myqgeftayfDnUqKHnnX74QZeUe3gYHZUQWSMJSjis11+H7duhZEndIumtt+DGDaOjMq/bt/XC2/btoVIlOHgQunSRYgjhGCRBCYdWq5au7vvwQ5g5E6pVg3XrjI7KPPbuhdq14bvv9LDe77/rNlBCOApJUMLh5cyp10dt3w5588Lzz0NgoN4Q0R5FRkLfvnpb9sREXb345Zd6aFMIRyIJSjiNp56C/fth8GBdhl61KixbZj+9/BIT4dtv9d5YM2ZAixZ6zqlxY6MjE8IyJEEJp+LuDsOH6/Lr4sVh8mRd+bZjh9GRPdyGDboIol8/3YX80CGYNAny5TM6MiEsRxKUcEo1a+oKv06d4NgxeOYZvZj18GGjI7vfgQN636aAAH0FFRysk1W1akZHJoTlSYISTitbNl3Zd+qUvqravFlv4dG6tf67UUN/Sum1XAEBuggiMRHGjoUjR/S291KhJ5yFbPkuxB1RUTB1qh46i4iAl1/W8zsdO+p9qDLjcbZ8j47Wc2MhIbBypV50/N57OolKA1xhh7L8UUoSlBAPiIuDBQtg/Xq9ANbNDXr0u6VkdAAAAfZJREFU0D3+AgIy1p0hownq8mV9tbRoEdy8CaGhulVT3746MebIYYb/kBDGkAQlhCX98YfuSHHkiL6ycXNTwGGKFj3BqFHtqVIFvL3/27UhrQQVEwMnT8Lx43puad06OH1aV+VFRkLXrtC2rfTOEw5DEpQQ1qCU7tIwaNDvHD+ezI0b+alatQbXrumhuXr19Bqr7Nn1uquFC/147bW9xMbqbg+RkXDpkt7F1tVVXxkppee76tfXSUnmloSDyfIZLRsWCpEBJhN4eV0gPn4os2YNYty4IXz11WqOH9dXQTduwLVrcOsWxMdDUhJcvaqTUYECUKoUNG0KpUtDlSrwxBOyq60QjyIJSogMev/99xk7diw3b97ExSWFatXuL/cOCgoiKCgIgCJFYvnuO4MCFcJBSJm5EBmwZs0aChcuTJ06ddJ9TO/evdm7dy979+7l6NGjVoxOCMckCUqIDNi2bRurV6+mbNmydOzYkZCQEF5//XWjwxLCoUmRhBCZtGXLFsaNG8eaNWuMDkUIW5blIgm5ghJCCGGT5ApKCCGEJcgVlBBCCMckCUoIIYRNkgQlhBDCJkmCEkIIYZMkQQkhhLBJlmx1JK0vhRBCPDa5ghJCCGGTJEEJIYSwSZKghBBC2CRJUEIIIWySJCghhBA2SRKUEEIIm/R/3EfDeI0gd8MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "values = {'a': 5., 'b': 1.5, 'c': 2., 'd': 2.1}\n",
    "h=plotting.plot(direct.subs(values), crossed.subs(values),(v, -1,1), show=False)\n",
    "h[0].line_color='b'\n",
    "h[1].line_color='r'\n",
    "h.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare stiffness\n",
    "\n",
    "Now we can finally compare the second derivatives of the summed lengths to get a measure for the relative stiffness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "kd = direct.diff(v,2).subs('v',0).simplify()\n",
    "kc = crossed.diff(v,2).subs('v',0).simplify()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2.50805786205887$"
      ],
      "text/plain": [
       "2.50805786205887"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(kd/kc).subs(values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For reference, this is the algebraic form of the relative strength:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\left(4 c^{2} + \\left(b + d\\right)^{2}\\right)^{\\frac{3}{2}} \\left(\\left(4 c^{2} + \\left(b - d\\right)^{2}\\right) \\left(4 a^{2} + 4 a c + b^{2} - b \\left(b - d\\right)\\right) - 4 \\left(a \\left(b - d\\right) + b c\\right)^{2}\\right)}{\\left(4 c^{2} + \\left(b - d\\right)^{2}\\right)^{\\frac{3}{2}} \\left(\\left(4 c^{2} + \\left(b + d\\right)^{2}\\right) \\left(4 a^{2} + 4 a c + b^{2} - b \\left(b + d\\right)\\right) - 4 \\left(a \\left(b + d\\right) + b c\\right)^{2}\\right)}$"
      ],
      "text/plain": [
       "(4*c**2 + (b + d)**2)**(3/2)*((4*c**2 + (b - d)**2)*(4*a**2 + 4*a*c + b**2 - b*(b - d)) - 4*(a*(b - d) + b*c)**2)/((4*c**2 + (b - d)**2)**(3/2)*((4*c**2 + (b + d)**2)*(4*a**2 + 4*a*c + b**2 - b*(b + d)) - 4*(a*(b + d) + b*c)**2))"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(kd/kc).simplify()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}