Skip to content

Instantly share code, notes, and snippets.

@moorepants
Last active June 5, 2020 14:19
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save moorepants/22415db65f0e2b5932079cd9b4860b44 to your computer and use it in GitHub Desktop.
Notebooks for lecture 18 Fall 2019.
Display the source blob
Display the rendered blob
Raw
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 169,
"metadata": {},
"outputs": [],
"source": [
"import sympy as sm\n",
"import sympy.physics.mechanics as me\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from pydy.system import System\n",
"from pydy.viz.shapes import Sphere\n",
"from pydy.viz import Scene, VisualizationFrame\n",
"import pythreejs as pjs\n",
"me.init_vprinting()"
]
},
{
"cell_type": "code",
"execution_count": 170,
"metadata": {},
"outputs": [],
"source": [
"q1, q2, q3, u1, u2, u3 = me.dynamicsymbols('q1:4, u1:4')"
]
},
{
"cell_type": "code",
"execution_count": 171,
"metadata": {},
"outputs": [],
"source": [
"kc, cc, mu, vs, m, g = sm.symbols('k_c, c_c, mu, v_s, m, g')"
]
},
{
"cell_type": "code",
"execution_count": 172,
"metadata": {},
"outputs": [],
"source": [
"N = me.ReferenceFrame('N')"
]
},
{
"cell_type": "code",
"execution_count": 173,
"metadata": {},
"outputs": [],
"source": [
"O = me.Point('O')\n",
"P = O.locatenew('P', q1 * N.x + q2 * N.y + q3 * N.z)"
]
},
{
"cell_type": "code",
"execution_count": 174,
"metadata": {},
"outputs": [],
"source": [
"P.set_vel(N, u1*N.x + u2*N.y + u3*N.z)"
]
},
{
"cell_type": "code",
"execution_count": 175,
"metadata": {},
"outputs": [],
"source": [
"kdes = [u1 - q1.diff(), u2 - q2.diff(), u3 - q3.diff()]"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"z = P.pos_from(O).dot(N.z)\n",
"penetration = (sm.Abs(z) - z)/2\n",
"vertical_force = (kc*penetration**3 - z)*(1-cc*P.vel(N).dot(N.z))\n",
"friction_force_x = -mu*vertical_force*((2/(1+sm.exp(-P.vel(N).dot(N.x)/vs))) - 1)\n",
"friction_force_y = -mu*vertical_force*((2/(1+sm.exp(-P.vel(N).dot(N.y)/vs))) - 1)\n",
"contact_force = friction_force_x*N.x + friction_force_y*N.y + vertical_force*N.z"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"loads = [(P, -m*g*N.z + contact_force)]\n",
"bodies = [me.Particle('P', P, m)]"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"kane = me.KanesMethod(N, (q1, q2, q3), (u1, u2, u3), kd_eqs=kdes)\n",
"fr, frstar = kane.kanes_equations(bodies, loads=loads)"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAhgAAACWCAYAAACcofGRAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae2dW9LVttKGgeJ6FyE3+/aHGZAwApIZJGQEZM+AFCOgkhkkGUGAGZCMAMgMSG73TQi1J5D/ffSp/cleB9s62LKXVOXlkw59UnerJXndvHHjxgMd73QcS6/++eefr4+9qPnZzZs3vxd8rwX7rzXD2WBLo4D4fEc1/LxFGU3DvJVOoUDTDynUa2WhgJehv3R5X8ddHU+khz7qvLmUokdV9r0QvncMadHj5u3gxU+6JnOY/ghvtnDtGX+jORdb4FYajHRo8ftHHe90/Vlaba30JVCg6YdL4HJZHNE5agGdg83E2eD+pY4vud9aStSjDOYZ6IXpoW6+4kHoYHyvhjbnUIRYidFPdX9PeGwu6hLi0a6nU0C8/pUOrqM5GdPJdpE5m364SLaXQprIhSVmAL61my2eY/WoyjknK8RZ/QznwjkYt8IXW772SD1rzsWWuRgHuxfyD5IBvOmWGgUOKND0wwFJ2oNICkjf/EfHd0Fxoqebn44voUfDCEZAr21dSnm4uXhB3SIX22JdTmjh/Z+Shbb2JidVd1BX0w87YGKlKHjZeizwHlUK4lywsurRvUQwmP96IQ9s817kXGlo+a8oIN6zwOqJDmShpUaBkAJNP4TUaNdZKCDngin533Q8l/75PUulK1eSW49u3sEQk5nr+UJHGLJamU2t+TUooM7xSu22qZI1iF9pm7XrB8H3LUel5Ns1WKL7HR3RAxLpmx90MD1yX/W83guxcurRzTsYYirz7j+JKJvcIrQXoawIDxzNp+rw9yqCqYGyHgWq1Q+SUUbAzOcfLJRbj1yX07K3GW6BeCLWyNgX4qdb2JhYVy3Fs+jRTTsYYiieP4akLe6rRSxXhsN73+yGahGtlXmxdvM16wfBRtT1mY69zN2vze6o9qUvmFb/RfyYFMlQPqIef+vg+1HDxPcwdpFy6dFNOxjiJEbkVxFj09trdyGRdSGBw0noebg/uy4oGzSlKVClfvByiUF7JN3VIq8DKRB9cL6ypnN1igc/qLF7yjM6VeX59UH5Q74BL/cvdOwpJevRzToYEgY8yBa92JM458PFOvqowsjXZKupJgpUrh9+Fq0YGO1iYWABvpeIPo7Vye6J7yU3UwYl5P2P8jIVixHmA1uf7dBZTNajt0WYrab/CPCPYmrbObIQB31nogPi2LlpiBo7FTAJVpT3NzoYnbR0eRSoUj9ILhkYMVcffqjp8rhTGcbSGX+IN9gSnD8ciJNJedEtu3cOc+jRzUYwxGD2Hjfn4mQ3yPvCOxc/SuhYlIbHznwjW7RqTb8IsAeCG2eopcujQK36AQPGfzy1ad36ZPK5QPpKOmNKFKM+6MtAlKRHN+lgSAAYBSAEu9kaVEY2stY6XCFNZ6zZgLNllZR9Pveq2vZbKwVq1Q/e2UV3jYXrayXtanBh9HWw44NpCY7suzbk9Flkok2tXnM6SY9udYrEjEaLYFwLwhJXNjVCWx99g+GzJWCY1AYjRCkhYCTa0rYBTqLabjLVqh9wLH5fKnoh+WcQxhoB+xNLDOgdtW9GQ7f1J48H//fxNU6A7nHSuHfTTP7+c92DL8+ItMZOYfDHZdCsTa2KCMiq6ButR7fqYGA0WH/RwowixBJJtB7OGdPJSW+vTlX+ApvBORtAdSycpy+Ee3NQZlMvvoDoTrQsxRDXqh+YtiHyVzyJhhjbpYxyaXww+MNFsaH+5z1/1skfHyI7TEPF/rsyg1aiJfT7NoAVMZSi9egmp0iEMN5qzYYNpuw90am/UyfEu601MYrBSZidvIJmJNSci9nUSysgmr9SDazoj+Kdys7SD/Bax2sdfN/ARvtpSAxKq15wwegvZbSmGOU/RGtG6kw1Y5RrTThm4XQ4EaqQjkQ27P7h4N0snFQPg1Ybsc8qu+PM0Xp0cw6GOiqdlGPx6IXaZgEQSoj2LzYJf1NetYcRnbEQvDFRDBawshOhpQIUEE/G+tITNTvp40cheL5vztIPMiqMhol6oFPMUIXV5rh20zZqJzZ0PxeGxYzyXMDm5A/4GdINXr2xeuAffVwHeumB7lPXuCADNs1mzVzyOVqPbs7BEJdtVFNkpDGUIgRcx0svvGx7THYuVBfK9emwrS3cC24WQDGPuwXja1GuWcrC4/hWOC7uxNYuAymyq7KT+xJGQ7Qg5D23n6ToBxzRcKSckx2E7BeRJ+istjiWNMo3UmTjFKG9HCALLqkNeDSMYLBWgCk1HAsiX0wNpSScl5hBSUqbNZeN0qMgtGUHY5HOioDrIASH8JZSPjULVwebOi4d+75o4ZwL7n2H7/JUdmEyMlw/MgamU1Rjmdr7eRSI6EusV+Bz2nOSORjG+0llvWyTt2QEYxZMkwA/kgk66/HSRvkIJNke2UewcDYZ5DmHwmoX7zonVLjDv9TdbY5PqtdkyZq61LPJ7Vw9ukkH467nsiF9qUxfFG/vSBCCfKNrIjAspsIQV8sHr2ihk8kM12eTKRWVrRavswjs6KXnH6vY50SgjNdz+eemSAKZyU1JjNVcmFJgWNoop8B6tqx4wgJPBnlMyRK57iIzkg0iNZ0TqnsiD+ECUN3OTsYn6rr4FPQJ61uTaXJ7cs4MGb3yfieAPwmr03M80IcIUfj8xDUCReo89Kvb9luYAqxJgPa9efFA+Ao3n1S9ycyUSnCcSo1iu/Yl8zhrf3UP1KZo2SnO4PnBpcqCD1NVjCjeq1y3Fsa/+1PP/28jvBGoZxO8QC9M5Ynx+uPZWg9f9sLuoiNygEJlaoMdCmZ0DktOf7LItC7gCF7o5WgmXJCVTra8jGCUndzoPodRVnWLJPjRyQIyLvi/02H96aHeP0qExHg926AmthtVXLjjvBJVRr6QfxamO557OdBtlmR9a3JlizoYggoifDgCHc9ZOT4lfeozHatnSvmWJ4ICEtSeUxhRxVpFMDRzFAUj2Z4TlRNwrwyYamMHjpN5PeOe0C/Kc0rC4PFfCBhFYO0cDF3bGpmegVVeU8A4JdDjiero5dGzGhNKE/0wNc3WD6INihMj+5xGdM+AB96guPkuAnP60buJfP2qYrVB0RJGGfyKJtERXjjDqWsWczqnSedo3hwDmH6h+nmFXFSdBKfpAP4LxTlGeoZO+Er3Dgnd5+j76Io5etTRbWkHA2J03icQCHmYSEd+w/2EVD3TJ+DQsixLgTkygyyWdF5xJvhUdOhQ03lRnqNJ/QX4bBEbI3sbbVlZHKRhH3NGUm06Raw6uEcJkbf2BC/m8G9OXsMdvUSyRaU/iVYYGdrm0/ipBmy2YnbQZPjxvC5ulDOAGlYxlGn3TnzA0ZzjbI7WGWY4cm3O6pFX6z8Sb5F1+jGDlSHNnPPl+c+sQY6+P7tvTXIwPCKEyOckRkgOSQr5OrpRQlBR17mDZ+cuXWdV3R/PZbJ3OWC3uuae1Taeo+EXFodRd/XeLVgKX+iavemjU0Uq+8+g3EXdikbOO5+ANEZijoIn71HZSpUllcew4yC4kbLOLk3ht+XV+YPym8F7rPteXbrnGxDDZ3rU+3MtHBSMzskkWIvJ7slGj79AcUKzqWmWfvCV4mjBcxcZsoZEZ54Zre1xzJn+TkIWT6ZU+bKKVc8x3fBA79/pnWU7eR7rW6VkI4Rb12fl8yTwZ15YnWP4BVUY34JHcZdqm7qSbOiRlumj7Ogbyig2J3wWLs4c7ftH2uHRXD3qqpnkYPiONjV8ewK+ztD2RlfKTOdmixGdOXvKBHsUXGqbRZAHScLG/O49vQ9D2wf5zj1Q2XFNca6C9m42BTLIEg4Bay2iZd3KSoZQIiitTpHoGUaEZ70+pjLDUV8vZK78B0llisnuQWPrP4CW0Oy9aEiE6bXwj+6bseh43qbqWdZfFNUNpWSjNNyxfMlRLhdvB7A4fRI+k/zijKMDkGNkYXbfp1yuNMnByNQYjgSj86FypXPbXDTe2PD9sHk3ChAhp+Qdlm33K1BAvHIOFU2Lv7awDCV+EKJXXkJ+udcHnB05DkhCXjpo1oS8qkKO36dUrPw2b+qyG92Cssf6E32JFfQn2/BwoJhSF8IFoBS9JCIxDP+ea9DxGjxFhzFdYpFVlDJTIUyRoItwNBj09Bw1a1TvxnhjWe1scLjoij1c+uzhRgbBF5rW/iXekyQqiYvqhkYk49vVXUW/HkbgdI5EABo6wBb5Bo+7WYSUvu/6Vq/SkZtbI+9zvgbxnqIQkRB0Dlt/MSUsZkxftbPmJMyF1IXRc6E68Z2RtvHRoc8zHafC8ikkQk56bY1Uhoxml63A2LFo8SAJd5wwl3RNHhtFY+SeXb3p/ULDoSOB0/GWXKrjWx0ooC7p/qluCNM+FzzDsl2+yi7QD3P4Z3mn8hCFGypkU6KOdqIZMtkl3U/hTZffX/TqHL5c4t7jwafvcaSQE+iDLGwuLYCLyQ68rj31bKqAZXr9oG+LZql9H5pY35pMk0UcDCFHZzVnIgTOwrD2D3lDYoV57fovf+EUgD1c6LxGmwuhVq4ZKTRGhQi+eds4m+ZUuoYxeDqQB1PG7nmmnylyZU3ROZPD1FbZ4Mx0Ri9qQ9/QwdoMB6NdixZu9Axd9O5YtKGHk8pBU9Zf2HM+iNZTCLr/QQe43Vd+44VuV0lT+9LodM4A+rn6AX50kYqAZkbHrvoZvOnKcBHU+WnvxbI3nQPrm2WdTuoHqZbF4Lq10riYbJbQRddYJFx5mUI3YFddknwyQEcPdPJ89eYqcpyh7x/0Cav/1Pn2qReZn4M0iRBkF14Uwmy1w0tknugv3U+Z9zSlaV4m9RZNXrHQnuHxm54xUuxW5xYFYB+VY/zMoUSps3edrVQ4HyUTymKOosDwYvCzJy/vfMqYPmCjI0L54TwpI+onYeN6jyIZJmj5s+piZEKi82OMqZ9nv+g4lVz7C9G/B4PahLZz+hJy06NHr8LDm7n6AfkY8hvamm7CEFuayhvLH56BqzMG4YsFr2nfjITRKXy2ICjJTYVw58bFbIvVmwxsoQoYeAx1AE2dGzzE9v07qneOHgWOG0s5GBgURqihInUA6NkUp8Ll9T+zkQwLx1wfgzumngsvgxInbE+nRcl9o+NcR9DrtKS26BSkyYpCvMYJZoqBRbimjK9qyfCrOs3JOlUbMB9zKHr5VQ84Hdtt1HvmacCHtx6pzLBeU6S9ukveCIYDHXCqPQ87DtgQ7lNFeD5LP6hudFMv6dkpnTSJN73Krm8YkKzmYAgnNz15DU73XxvAtam0AC7Gp6ppIzr0dID6i1tioOcugpGr7/t6kBHam5WWcjDw/MORwCwgB5lN6ZsQDF5v4vYghLUJqBOAlNB3ux0Sqplb1GTEwuZTy+Plu1Hs1AIZ82FMgdvkHGfHzZ3HtIESUnmMbqgciMRx/yKiziVl95ngG0YXxkA2uhnvx/LPeZ/CG8pOWWM2B56UvE7GvZFKqScsu6RshO3mxgVnjMXSYZ8J26v12g3kDbiMfd/60lw9Wj6CIeUGcHd0ZBE+EY31GtCQcPAm0wYFd5N0FtDWMczoTMJD/PlJMvYa2dX1rLKTGjifiQgE00cPdEbB0XfGoh7KcjZRJyF/FATrAKALX/6brUBjyqit2UmwgvcXam9WP1f+kvohhTesOXoKXkvR8BTRBQMGme3SpyI1p4qefb4GXoVwwQHPYq/OEiz/S+AeDuRy9P0oPQp6S0QwAI7tqXjwuRJKn7nZlnZCASkKDCodBHlhHQG7KFLXZzz05ImRPTomW2YPQui+ziIn4YxsT55GmAKE73sxNJhSfak80B4exKQi+iGRN2awiOYOjUAMjlFl1K+IogzX/UTVtXahErioThxbdNHcyNlq5BDMOIzoTmDnjyg7/mbq+9F6tLiDIQTpWMP5Pz1KStRZU7gxCZktFfYdEIE+lXqOgfIf+6rgqbL2HN6yXsPdS4ZiPxzkPG9vGKzuSWeVYWqBUf+3ul7NIEwCdmeZoLlQYjtlbPSoOv3g5Qm4cFhXkSfRFSPE7iLnwPr7D7rfmvPJtGEpXKiXBK82kcS/1AjnGJ7RerS4gzEGeeT7d5STkK0eboyEP2sx0QGPm88cFx9tqw3C6pNH2Mof6xzkoBHKIlpRCHYM3CrGIAfyW61DdE+lea36gajMOee8GMtMR6iB57q2bZ7049goUTFYxyoujItbx+D7/hgol/I+Wo9u1cGw1b0gnhpG36yQ+I5mxt55mZtFJjPgog30IGRYdKdKZrBbdXkoUKt+YFHtjxh4GbCl9RYf1aI/4OR0SXAwYNhaKonLwVbkrREnJ7ypevSWgPmXB8jOOeErUpc6BSE9OsY3RRrYSKXQQQcORk9pnANfAjM577l6NvDORmlLK/INkGbfINaqHwQXOotwNrtjFk1q+xMdN4fHokBkaqwULtKNTwFR9Tedcc2rGD36byuOg/E/f2Nne1f7mdGAIV87rDXBxyhm8cSoTYd9BIrpLdZYuA5dCBhCnSwujp3HLwRWq3YhClSpHySP7Ny4J9lnWrOluiiA41d6PUNdGI9DE6NH/2vV4mBsNbmRuDoq0yQtZaSAaMrnq50zoDNb63I5AkxX2IJfDABbJksl5GIzK8FLEeGC661ZP2DEfr5g3lSHutdxDEhS1/9Uh1siQEl6dLMOhgSBxXuMTje3SCmR4UWLq6MR4WCOkz/EYrRFh7NtSrqMS6qLsCO8sjURd4NrtypcbRPVINKR5NAE5ZuyiGPX5kvVrB8EG3L5IZDTzdN7ywiID/cEP9GLZksCRgbyGa1HN+tgeDqwIhujtErYP+DFni4ZWeG4PRZd2S74WAoxV8f7XPXZrg4Wutm1HrudKXz8B0eEueqUhLJ4pbpS60mBoZVdnwI16wf6FNug21TJinLibQeDnifSF206tc+LZD3KFkIEnG1d7I/eHIElIH8L9p8E+8XOnYkGLowlGtj0g0jitvESJsY7DxP3Qz7zKWrnpXp6fq370PiH5aOvVTcOC84ghh9nooND75BDFqzihNDZo/bmqx7W5YD3JuVZcLeUkQJenqvUD4KNvohx+zLsCxnRb1WNUEA8gP4vRf/oUfpIE5t8naJHraxo6r5ktHUHwwwKq6QvcsQqhh51MI5JNh1KdDr5vQy958NYPeOsZ0W/N6L6GWni5PxBW7p+pusoh1HlcZbfqrxt39VtS5dKAclD1fpB8OFk0L+iHOpL5WsOvL2u+Vy0zz6YygHfmnWk6FHrc6LrzdtrIpGjbSHxSgjROTFSl2pUMMq5Ep0Np9NFF3wnZNok1zTJMTj56/Iv1BbvwOU5F3OTF2wU9qO5ZVv+fVKgdv0g+Loo3j45UC9Woj0D0uZcDFiUU49u3sHwtMH4vRNhGAVfzEhA+GJMGekTwWDrG1MDb0QDFmfGJmjJDhLqdkn1lXQu2HeeHJ4UvDgmzhFSfSiOlhoFjAIXqR8M+XZuFJhKgdx6dBcOhgwKofUnIqKbe59KzK3nA2/hkDVq441z1joXojPOxQvB30YkCxF8K81cqn7YCn8anFVRIKsevVUVagnASImw++BXORo4GS1dEAXEc7d4VDKwRcfogji1HqpNP6xH+9byNihQQo/uxsGAhWZgRKik7yhsQxyiodyVERavmR7inzdPLlyNplQruCsKNP2wK3Y2ZDJSoJQexcEgzM48/oeM8K5WlZQI8633PcFWg6PWhkWf3SwqE49Zd4HD1BZ11ipwlcHV9ENlDGngrE6BAnqUdZBuFyAOxj0d7MC4q2MvqS3y2wsnJ+Aho9H4PYFOLUtHgSYvHSnaRQwFZJTtrxT4d9yX3kjHVFVNmYx6lF2I+BQ3djVFAkIwnrOI1Rb7QYgdJ98h6OB8+6KlRoFRCjT9MEqilmGEApIh/uPovfTPDzqIoBL93+zav5J6dFcOhhjP2ot7IljUR5pUtqWNUUC8xpFsTsbG+LYGuE0/rEH13bYZfjWZAQ5rwTabSunR3TgYUh58sY8vQBb9ZsNmJWjHgIvnfEeDP49y0asdo9pQi6RA0w+RhGvFDiggfcOi8nAQ+5kybT5iXkKP7uI7GFIeLPYr/bXJA0FrD6qiAI7ln5IFPoW++c5eFWU3DkzTDxtnYMXge9l6LBD3stA8qx7dSwSD+a/2kaWKO2Jp0ORUsHDPPrZWurlW/7Yo0PTDtvi1CWjlXDAl/5uO59I/u/iCdG49unkHQ0xmaoT5rzBkpduWLo0C6hx8bK1NlVwa48/gW7t+EHzfcpxBob0qRAHR/Y6O6MWZ0jcs8mR6hM8ivC4E5uLV5tSjm3cwRH3m3fk75rb1bHFRrLJBHM2n6vD3qoSuAbU0BarVD5JRRsDM5yf/F8/SRN1De95m5FggjozxZ40MdveSsujRTTsYYiieP4akLe7bi1gn4uG9b/t4XGJtrfiWKVCzfhBsRF2f6djL3P0mRUX6gvVav4gfkyIZykfU428dfOthmHbzLalcenTTDoa4i5f1q4ixm69TDiW23UdRAIeT0DOLf1u6XApUqR+8XGLQHkl3tcjrQD5Fn+xbPs/VKR7w79P8G/XoVJXnF9+9CPkGvNy/0LGnlKxHN+tgSBjwIFv0Yk/inA8X6+ijCiNfk62mmihQuX5gxxsDo10sDCzAdxzD3GmsTnZP8HXOKYMS8v5HeZmKxQjzP0if7dBZTNajt0WYrSa+oPZRTG1bEhfioO9MdEAcOzcNUWOnAibBivL+Rgejk5YujwJV6gfJJQMj5urDDzVdHncqw1g64w/xBlsy+rkD5UW37N45zKFHiWD8y/Pazv62+hN7j5tzsRCbvHPxo4TO/rmU+Ua2aNWafhFgDwQ3zlBLl0eBWvUDBuwVBu3yWFI9xs8F4VfSGVOiGNUjkwnAGD36b2sbB+N//sbO9q7aswSAUQBCsJutQdUS+xqw4QppOmPNBpwtq6Ts87lX1bbfWilQq37wzi66ayxcXytpV4MLo6+DHR9MS3Bk37URRCba1Oo1p2P06H+t+FanSMxotAiGcXKZs02N0NpH32T4bBkoJrTCCFFKCBiZH23bACfQbEdZatUPOBa/LxW9kPwzCGONwHvPW8L6d9S+GQ3/uO6Tx4P/+/gaJ0D3OGncu2kmf/+57sGXZ0RaY6cw+CMzaNamVkUEZFX0jdajW3UwMBqsv2hhRhFiiSRaD+eM6eSkt1enKn+BzeCcDaA6Fs7TF8K9OSizqRdfQHQnWpZiiGvVD0zbEPkrnkRDjO1SRrk0Phj84aLYUP/z/nv101+97DANxQewYhKDVqIl9Ps2gL2iYLQeZYpkiwlvtWbDtkWazoWZTv2dOiHeba2JUQxOwuzkFTQjoeZczKZeWgHR/JVqYEV/FO9UdpZ+gNc6Xuvg+wY22k9DYlBa9YILRn8pozXFKP8hWjNSZ6oZo1xrwjELp8OJUIV0JLJh9w8H72bhpHoYtNqIfVbZHWeO1qObczDUUemkHItHL9Q2C4BQQrR/sUn4m/KqPYzojIXgjYlisICVnQgtFaCAeDLWl56o2UkfPwrB831zln6QUWE0TNQDnWKGKqw2x7WbtlE7saH7uTAsZpTnAjYnf8DPkG7w6o3VA//o4zrQSw90n7rGBRmwaTZr5pLP0Xp0cw6GuGyjmiIjjaEUIeA6XnrhZdtjsnOhulCuT4dtbeFecLMAinncLRhfi3LNUhYex7fCcXEntnYZSJFdlZ3clzAaogUh77n9JEU/4IiGI+Wc7CBkv4g8QWe1xbGkUb6RIhunCO3lAFlwSW3Ao2EEg7UCTKnhWBD5YmooJeG8xAxKUtqsuWyUHgWhLTsYi3RWBFwHITiEt5TyqVm4OtjUcenY90UL51xw7zt8l6eyC5OR4fqRMTCdohrL1N7Po0BEX2K9Ap/TnpPMwTDeTyrrZZu8JSMYs2CaBPiRTNBZj5c2ykcgyfbIPoKFs8kgzzkUVrt41zmhwh3+pe5uc3xSvSZL1tSlnk1u5+rRTToYdz2XDelLZfqieHtHghDkG10TgWEhHoa4Wj54RQudTGa4PptMqahstXidRWBHLz3/WMU+JwJlvJ7LPzdFEshMbkpirObClALD0kY5BdazZcUTFngyyGNKlsh1F5mRbBCp6ZxQ3RN5CBeA6nZ2Mj5R18WnoE9Y35pMk9uTc2bI6JX3OwH8SVidnuOBPkSIwucnrhEoUuehX92238IUYE0CtO/NiwfCV7j5pOpNZqZUguNUahTbtS+Zx1n7q3ugNkXLTnEGzw8uVRZ8mKpiRPFe5bq1MP7dn3r+fxvhjUA9m+AFemEqT4zXH8/WeviyF3YXHZEDFCpTG+xQMKNzWHL6k0WmdQFH8EIvRzPhgqx0suVlBKPs5Eb3OYyyqlskwY9OFpBxwf+dDutPD/X+USIkxuvZBjWx3ajiwh3nlagy8oX8szDd8dzLgW6zJOtbkytb1MEQVBDhwxHoeM7K8SnpU5/pWD1Tyrc8ERSQoPacwogq1iqCoZmjKBjJ9pyonIB7ZcBUGztwnMzrGfeEflGeUxIGj/9CwCgCa+dg6NrWyPQMrPKaAsYpgR5PVEcvj57VmFCa6IepabZ+EG1QnBjZ5zSiewY88AbFzXcRmNOP3k3k61cVqw2KljDK4Fc0iY7wwhlOXbOY0zlNOkfz5hjA9AvVzyvkouokOE0H8F8ozjHSM3TCV7p3SOg+R99HV8zRo45uSzsYEKPzPoFAyMNEOvIb7iek6pk+AYeWZVkKzJEZZLGk84ozwaeiQ4eazovyHE3qL8Bni9gY2dtoy8riIA37mDOSatMpYtXBPUqIvLUneDGHf3PyGu7oJZItKv1JtMLI0Dafxk81YLMVs4Mmw4/ndXGjnAHUsIqhTLt34gOO5hxnc7TOMMORa3NWj7xa/5F4i6zTjxmsDGnmnC/Pf2YNcvT92X0LBwPAmEs/qVQ9IoTI5yRGSA5JCvk6ulFCUFHXuYNn5y5dZ1XdH89lsnc5YLe65p7VNp6j4RcWh1F39fA0Rd8AACAASURBVN4tWApf6Jq96aNTRSr7z6DcRd2KRs47n4A0cj1HwZP3qGylypLKY9hxENxIWWeXpvDb8ur8QfnN4D3Wfa8u3fMNiOEzPer9uRYOCkbnZBKsxWT3ZKPHX6CfoNnUNEs/+EpxtOC5iwxZQ6Izz4zW9jjmTH8nndSxvEyVL+ogqZ5juuGBXr3TO5fn3M9Y3yolGyHcuj4rn+fgP/XO6hzDLyhvfAsexV2qbepKsqFHWqaPsqNvKKPYnPBZuDhztO8faYdHc/Qodh+f4gYOBp0XQBlR0aEOku9oU8O3B+X9AzO0vdGV3tG52WJ0tO1TlU19ngn2qc318qltR+TeQ91I2Jjfvaf3YWh7mO3svcqOa4qzNSz30vClRcNZz17regsj6I5QghcZTekHOASstYiWdSsr+tGf7ujoFImeYUR41utjKjMc9fVC5sp/kFSmmOweNLb+A2gJzd4jlzojm9F9MxYdz9sU+XJNq56iuqGUbJSGO5YvOcrl4u0AFqdPwmeSX+w5OgA5Rt/O7vuUS0zoIXyKH3AwlkoYE0bnQ+VK57a5aLyx4fshfG4UIEJOyTss2+7XowBerYvMeEPY47N/ZhEdBJTQ9DDsFwu9k5mJhclLB82akFdVyAEdRpPy00H/soyixdDgHetP9CVW0J9sw8OBYkpdCGeglT7fVQNz5MDxGjwn6BKLGqCUkTemSNBFOBoMenqOmiGqd2O8sax2NlkHl9WShxsZBF8XuZ5Co9UAPtNwSVxUNzQiGd+u7ir69TACp3MkAtDQAbbIN3jcRchS+r7rW71KR25ujbzP+RrEe4pCRELQOWz9xZSwmDF91c6akzB7r0tKDKWNc2GdAVkwnhv6hKeZS2TkzGF57X3sGTkxmZlSBzKaXbYCRc6ixYOkvkBUyyVdk8dG0Ri5Z1dver8PdDd0JHA63pJLdXyrAwXUJd0/1Q1h2ueCZ1i2y1fZBfphDv8s71QeonBDhWxK1NFONMOZ6JLup/Cmy+8venUOXy5x7/Hg0/c4UsgJ9EEWNpcWwMVkB17Xnno2VcCiZw/6tmiW2vehifWtyTRZxMEQcnRWcyZC4CwMa/+QNyRWmNeubVTnFIA9XOi8RpsLoVa8GdYG2IgQBWd/TGQNW6fmHjlAXnKlKXJlbdE5k8PUVtngzHQGuHeJvqGDtRkORruWEXC08o7AsWhDDyeVw2mDxvacD6L1FILuf9ABbveVP5cDp+qi0tS+NDqdM2h9rn5wsmh1BDQzOtornDbHJ+UZ401Xhougzk97L5a96RxY3yzrdFI/SLUsBtetlcbFZNMcw+uWK7nyMoWu6vSk5JMBOnrA9GwHrfLn6PsHfaJr4MTF7RPPcz8GaRJGpQsvCmm22uElMk/0F0Qg00gypRkapJEiaa+9YqE9w+M3PWOk2K3OTWvhIkqjnBlVwz8ElemQzsiJ96FRh86zhVlljiWUxRxFAUzAmj15eedTxvQBGx0Ryg/nSRlRPwkb13sUyTDhnP+suhiZkKAXNKR+nv2i41Ry7SsfW9mILi2W1Ca0ndOXcJp69BgBdq5+QD6G/Ia2ppswxJam8sbyh2fg6oxB+GLBa9q3fmV0Cp8tCEpyUyHcuXEx22L1JgNbqAIGHkMdQFOdXj3Sbmzfv6O65uhR1/RSDgajBOY0Q0XqANCzKU6Fy+t/ZiMZFo65PgZ3TD2XXEY07BYjnqODDBCC/EwHMpOUfF3UMVlRCE6cYEarLMI1ZZwER1hYdVrULnwcXoP/MYcizGOj4mO7jXrPPA348NYjtT2s1xRpr+6SN4LhQAecas/DjgM2hPtUEZ7P0g+q+0DO9OyUTprEmxPAMSBZzcEQTuFOAkB84OF0U2r+ehOnBXAxPlVNG9EBvdb1d/UXt8RAz10EI1ff9/UgG5P1qAnSLbsofMbzPzeimtO8KX0Tgjlla8mLAEwyuLUAvAQcXpB/VlsYQ+NzStMmIxY2n1oXXv6YIzC1rrn5MKYGtysrugxH2JPr9EoIoxsqByJE3L+YXNF1xiVlF0dzLu4mNz0aXoOfdJXCm4OySZCkF3Yy7uUjvbarGpaUjRDm3LjgjLFYOuwzYXu1XruBvAGXse9bX5qrR902VYOnyFnKEeDu6HBeVWojIhrrNagmDKmnVrto+Q0KbnH6eOcCRfGExnTvFn0mNmwdw4zOpOrEn5/U/msdRaIYI0AwIuHTx4wwUXD0nVRnhzoJ+aMgWAcAXfjy32wFGlNGbc1OghW8v1B7s/q58pfUDym8YVHzU/BaioaniC4Y6Gdslz4VqTlV9OzzNfAqhAsOeBZ7dZZg+V8C93DgmqPvR+lR0FtiigTg2J6KB58rYTCYm21pPxR4J1SQFdtJBI9TDetDT54Y2aNj8pW8gxC6r7PISf0EvCdPI0wBwve9GBpMqb5UHmgPD2JSEf2QyBszWERzh0YgBseoMjLI9K/hup+outYuVAIX1Ylji3M/N3K2Gjm8k4VzAez8EWXH30x9P1qPFncwhCAdazj/l8oM6jRDlFpXKz+DAr4DMgI6ldhe2S0cVP5jXxU8VTZ8fs/Kqr7YDwc5z1vlMTizksrwqWhG/d/qejWDMAvonWSG5kIl5Tso1ekHL0/AhcO6ijyJrhghdhc5B9bf82XYrTmfRDhL4UK9JHi1iST+pQ7ExvCM1qPFHYwxyCPfM9pFyFYPN0bCn7WY6IDHzZRC8dG22iCsPnmErfyxzkEOGqEsohWFYMcxWcUY5EB+q3WI7qk0r1U/EJU555wXY5npCDXwXNe2zZN+HBslKgbrWMWFcUGHsiFh9qBkDO4Nv4/Wo1t1MGx1L4h3o+UNMzAKdN/RzNg7LzOqoh0WEm2gByHDc1u2doh5Q0kUqFU/sKj2Rwy8DNjSeouPatEfcHK6JDgYMGwtlcTlYCvy1oiTE95UPXpLwPzLA2TnnPAVqUudgpAeHeObIg1spFLooAMHo6c0zoEvgZmc91w9c9+hVHXYNxqIPvFNjKdz65mR30ZpSyvyGSC2rCUoUKt+EFzoLMLZ7I5ZNKntT3TcHB6LApGpsVK4mD5S/U1nXPMqRo/+24rjYPzP39jZ3tV+ZjRgyNcOa03wMYpZKxFNsPU48I8dDaUSoU4WF7dQZykK111vlfpB8viDyMb6IqY1W6qLAjh+pdcz1IXxODQxevS/Vi0OxlaTG4mrozJN0lJGCoimfL7aRRt0ZmtdcqTBjwqY77Upi7vBtVu0pXaIahDpSG5PdSMXm1kJLlhbykuBmvUDRozvvbRUCQW8zmFAkrr+pxKMsoGRpEc362BIEFi8x+h0c4uUsrG+QEXqaEQ4mOPkD7EYbdHhbJuSLpPS5yptiy6Zh7ZrKmWqh735hCcJJUenwEFpyiKaitsuWLN+EGzI5YdATrdN7I1DLz7cEwpEL5otCXgZyGe0Ht2sg+HpwIrsg3+NDGjULudTgJEVjttjCRjbBR9LIebqeEQU4Bf1Duc5n+sZH5hiB8BbHSkJZfFKcCc5KikAtLJVUKBm/UCfYht0mypZUVREfwZURFWfSF+06dQ+L5L1KFsIEXCUOvujN0dgCcjfgv0nwX6xc2eigQtjiQa2vkEkcdt4CRPjnYeJ+yGf+Rtn56V6en6t+zC6EJbPfq02MQTA8Ieu6fDPdB3FT5VnXQ54b1KeBXdLGSng5blK/SDY6IsYty+R/Yxot6omUkA8gP4vRf/oUfrEpjaVLUWPWlnR9ObtTWF9HNgnevxSSBHSbyPWgEaix0HkgQ6l5yzcOZUw8D1lpzKlvzfCP4t+oXaAifaJZsQmvG4MSg+H2Mpauc1ToFr9gIxK5umLyHxLC1MAvaYm+X7QYoOphVFMaS6LHt28gyHheCVBYdsqo2Dm8S8x5VRQdDaiWs5A+07ItMmBs5KL0OJhltGDYCV6wajwUS7YWj3bpkDt+kHwNUd4JRET7T+q6eZcDOifU49u3sHwtMH4vRNhCLPjbFxEEr4YU6YSmCJh6xtTA29EAxZnxiZoyQ4S6nZJ9RVzLqyN1LPgxclyjpBXHKlVtvL7ocBF6of9sK9hshQFcuvRXTgYMiiEGl0oVIzorUNYijFrtAPeajdr1MYb56x1LkQbnIsXgr+NSBYi+FaauVT9sBX+NDirokBWPXqrKtQSgJESYVfCr3I03P73hKpa0Y1RQDxnV0r3D4IbA7+BuwAFmn5YgMitiU1ToIQe3Y2DAWelRNzIW4TK8aGmTQvLGeC3GJ04iY54zfQQ/7x5buHqyfLtxeVQoOmHy+F1w3QeBUrpURwMwuzM43+YB1KduaVEmG+97wlWJ5ArQiX67GZRmXjMugscpraoM4NMiZ6svcFh221q+mG3rM2OGPpFx+4j4uAp4uXUo6yDdJ8Z2Px3MLJLVauwUeACKSAlwy4sooBR3x+5QJI1lI9QwMvRX3rFWjj+DoAPWLFbY5NJ+OBws5X1s00isALQopn7FpFodnNXUyQr0LI12SiweQpIITCleK85F5tn5aoISI5+FADvJUc/6GBETFR80xEA4cGi8R+FGx+jbGkmBZqDMZNgLXujwJ4o4EcbfDm1+q3Ie6L7jnEJd/FhlDc/5aa+wXd6+O8YF+XbMe+yo3Y7e42twkaBRoFNUEAKs/t2yCYAbkBWTQEftQhhZFphL9vGccD/VJ/hS8h7wSnkVZHr5mAUIWurtFFgExQgfN2+HbIJVm0LSO+8PhbUu1iALafio3Cyby19si1urAdtmyJZj/at5UaB1SggZclCLMLXbVHnalzYZ8OSLdb0/KaD/4fazZeVhQvfWmpTJTPEtjkYM4jVsjYK7IgCzCfzp3BVrvCXkfqWY0f03gwqonvS9kzJFIs8mR7hcwGvN4P4NEBxyJ8Kr3vTsl92ruZgXDb/G/YXSAFvuFGQVS5aE3yMgPl4WpY/wbtAFieh7J3OHDsnkC/+JZlo2S6SaEMUw74dtQucSiLRHIyS1G11NwrUSQFGYb9KWVb30TUZI6Zt+KvoXczd18n+cagkGyxk/EX8mLTNVPmIevytg39iHia+h7GnhONEhI1F0i2doUBzMM4Qp71qFNgbBbwBqDJ64RU2Bu2RH0XvjfxJ+Ig+OF9Z07k6xQP+lZl/aR6dqvL84rsX4ZQb8HL/QseekuEzSpc9IR2DS3MwYqjWyjQKbJcCfADpox+h1oYF/+RIZGU3CwMzE7jEgtyxOtmeySfkp4zWyfsf5WWNAqN8/h/oM+986HIfyeODjH6zD4zKYXG7XNWt5kaBRoEKKcDWwer28csgEVpnrj78UFOF5LsskGRM/xBvkBecv7MfY/OO4aU4h7+IHjhefAG3uqnGWqSUCMa/dHzw51rganA0CjQKZKaAN+KMRGtc2Y8Be9WUdWam56nuuar5SvIzJYqRp8X6a2GxJyn7tNVVtZv+/begx6e4QQTjfzpYhMO5pUaBRoH9UsCUYVURDEaBIjkRjLMj5Fxs8YaSEP57Xyej7jtybsxo5GqqeD1L4EJkQu1AI9YcsC7j4hOOsGjC+hKmgdpup75E/Fe3+BTOwei/aneNAo0Ce6UAypD1F7WFdFkH8PsScHmDzH9kfO0NJ44N95ubmhnDRe/B7XMdRB7A70dw1jkm8UdmOGXNwbim3ltdQuOWTlCgLfI8QZj2+DQFpLgmrSw/XUN7E0MB0Z0wNaP92ISxQSnWllgXwpz2EgkjOVxI2jldGGUdbEFkoSLfgqjZgJzFRXjy/g85FTgFTIsxDRWbiHqxFdWiYLH17KkczlpKf9wTLY7i0hyMo2RpD09RQAqG0RAjoRYWPEWkQs9Fc0L4bmHZ3CY83+Dd5OgFZXS81sH3DWw6YW7TZ/OrXhQ0cC01bYMzE65BwWCGbec0yqq6aBrDhSiN4fZQkNj1bKBUD3LzUQdRsJauKOD6hGS4Zid0VV41B2NV8m+ycf5jgK2OLRWggJQVUQoMOkb3WHqih5M+fjQobCOtyY6CjAojewwKxiXaOA3gGN66EbHaiQ3dD+s7ee9pCl3DtsDvTVAom1EO6sx+OQUX+Kd8RGRwmh7ofmxL6hicyECLYFxTyaKBjSbXNOldNQejR4793kjJYLiepmCo8izyeitFNXkUnNLelsqm0FdliRS81IEhYG/9KefiBkZD739V3rm8NAcjhneM0MJRv26zJf6zIgam2QB42kE/l0RD8OpFMMjDc8+LHEb5hupK7nse5O40BRcyKx9rW3AsiHyx1iQl4Yi10fo1BU1uN7d+5xqFsle3y1bfat8ZBVBULUSamaneWLgdFDICOHFj/93AtsE/dcxZcHfXg21K0d+ePwkeG52VjGDMguk8xKNvoTPGFmP5Kbkxwpwt+Xt2TvA/Gu90jxNUYzqLi2B/KtidjOiMU4rjlPLdBsenxDpqpGMUTKIpzihlrW9F1bPnQi2CsWfuZsQNpUJ16lRLGoOMGOynKhSbsGGbnBn/KchZVKQbwU8ppDxuisS3ObHIrGzI1WIyJTxY4Mk0CIaX6aKecyGadpEh8uq9M8o6V5fO4SI84Df/6eKS7ok8dItZ/eO5J+NTi2L0KWd9q/+03bVtqqVkQB2aUSjK2ea86Zz8Q+SXeodh4B3K1eZFbdTKYqw3KECvFAiZ/6WD5yyuLDWSVPVnE/AVb1s4M00AvpaGK/7t+cFZZenoRAAIWb6HhpbJv2PU/396PtfIWjU1neEFI9ipPHGjdeX/MBMJZLVrQ3REDhixMar/XrQ0o6Pb6GR9JLqCyILgEOJmRtnJjXDNYZQjQZtdrIcLMi74v9Nh/Qn98Wh2rf0Cxmv4X30S7uhX1oshX/CWhenoB6aNOr7rPiWhSzZBjxQkY8veji3Yyp2mgO/UzN+66QQv6Ai57bz4Uu/o/Mxrs7/8te5NqdERWOSnk+sIzgHRPYqdxX2f8GKFBC60XyR5GkGP70QLdkswd809DhbKc0rC4PFfCBhFYO0cDF2jWO7oPQqhS8prChinBEXxZJiny1zXBfI0Z7EtcjUriTaUwcgyJQM/GN3DGxQ3csucvsm0LuclXz+FejyZV0tcbrUN/M7Y6Jq+itEpYZTjAJxR6hguFBc+0bw51rynD69my9Kx+ko+E01MB/BfKM4x0jN0wle6d8pV97n6fvX0KEnrc3U3B+McdSLeSWhxBFDEnSOAgOs5tbHlD4X9hhslFDVhy9AQmjf8qZ6HXjbP1xRkYJ07+lWRyQlngk9FO+fCl8LwYAhGk+gKfLaIjZG9jbasLA5SSE8MpjOSatMpYn+PEnKOoRWs9Awv5siDkyvhCk2nJpQ0yRaV/kR50Ym2icalGjCTddfIkj+CHefswEHLgFNJNIYy7do6hctEQI7WOVL205H3q76WfNIv6McMVob4/Q5wGfs+fWE1OQaXmlNzMALueMH8LXg05ZIRrxNan/lnnQnrd4pc9T7w75yBI79vC6M4NGaW140afTlOBwYyeNddql68cjMM3XNd0Onu6j0RgWHiYzwY5XOJTtThFGb0uETTTeUx9NCih/MEmEIwPii/GbzHw7p0//mRZ5QPV4DjoDCqPZkK0vdkmydeoDihWcmEzMFzFxmyhkRnnhmt7XHMGZkkoaTPplQZo3LV8c/ZRkZeCm83SjiVrZRshHDr+qx8noLt3HOrcwy/oA7jW/Ao/lLtU1+0/jjSMjqQaOVQRtGL4bNZff9IO+3RCAWagxEQyCvOqeH4oOTVpToKzgGdBe85TAg2Rvxj8NA5AXrWG1XrPUp9mJdi5O8ZYB4Ok+pzUyrD54KNyAoryMNoyTBb1L3HK5puahSHoOeUzQXEw4ARgU7woFMkAV96tFaZ4QgWHHp5hnCozOL0HcIQee+MuGhxME10pj5oCT3eqxwRpm4q70yZIq8yyBhTBmcdhFTAS8lGabhT8U4tn4O3AxicPgmfSX5xxtELyDGyMLvvU+5EGnWQT5Tb/ePbu8dwWQQZJZOGRupY9IFnYeTDFdQPSj2cJgiNpq1NmGMkrN7UM52IDpo1YfBUIccxWhy0pfw2b+reSVEMHaZjDho0ZSrqZBseDhRT6kK4A5gLPSCiNAz/nmvqo395MhIVFvb0QCkzFcIUCbKHo3FygZzejfEmbILrEKbhu8XuPdzIIPhCU0LrBtticORoqCQuqhsakaqljYcROJ0j4aC9+kEH4FgMdTP6lfyxfX9uP7yC5kJ+b10InkuhiRFGiIeKH+EeCvwxLxtBR8kN8zJ9gWJnLQfvXWfReckETnSmrClQ5CxaPEjCl8iLS7omj42iURTPrt70fh/obuhI4HS8JZfq4H8moHOXdP9UN4RonwueYdkuX2UXyMHHGTDZzpwe7mfKI5+hQnayrUeuvGiGM9El3U/hTZffX/TqHL5c4t7jwe4st8NLbSLjyMLm0gK4WP+H17WnoQ52OnQItGiWo+8P2xo2c7H3zcHIy3rnHUtoUf4u6dqmSzrP2b9HUQ8diVNeNpERK48i7EU4rloq/ovhTZkGOQcg0xk4AV0Sje7oYG2G67x2LdwdHbwjcCza0OvsKgdNoZ89v6+yPcOs+x90gNt95R/yRI8XTVMdgNHpnAHUhrMZicHrg1v4YTKHo2HljY5dgRm86cpwEdT5ae/FsjedA+ubZRqy2m9fjJCmNC4mm+YYjoCz/GsvU+iqUAezbgU90MmzQab8qX0fmlRLD8NzrfPttRreY7sIt5Qtijn8UiAKebimAuHn2VDged6tHdC1JZQe38+go2B010gY3iJtiw5sLYVmjIptdMQ0UDhPyoj6SYi43qNIhok1Ej+rLkYmJOiPMaZ+nv2i41Ry7SsfW9kWdeLUJrTF+DsnU+ff9IyoC1+SPCYTOE09euj+XJqrBFGcQ35DW3gFj5BJS1N5Y/nDM45LZwzCFwteu/7o2wMeUvjs6sk2fkO4c+NizqnVWytFGHgMdQCwnhs8zO776gd3PAFqp4cHc/lTczAy01zGAKehcxwkhAh1d09zPk+4gpnHPB+uJ7DnGLtFDZ5rOPgBZuHCFEPKp4aDGvuXqv/o4skgF535mEMRZOlGxcd2xPSeeeXAh7ceqe1hvaZIe3WXvBEMoTN1tikPOw7YEO5z5SzyMMmYq24c5V7Ss6PyqUyTeNOr7PoGJ2oSTNdF8l0Jp2E/fOBrB65NpQVwMT5VTRvRAYPf9Xf1F7fzRs+dHs7Y940eNv24KXlZAthbSzRyKW1IcE05OZS9IDMitWmSNUlB5zo2Ep4DE17+mCMwp745eTGm1qFdOdF3OMKeXJ9XQh9UIBx9wCvuX0yu6DpjDvpe13b+6plez8Jd+JozQjQnd0rhzUHZ3MDNrM/JuJePmUVPZl9SNkIgcuOCM8Zi6bDPhO3Veo2zbPJvU3MfwCUAOKbvmz4y5z2orl1CgRbByCQHMnZ05qc6fxJ0QBwLPk7kPOdMTUVVE8AUVZ5CquMn4cfHwopEMUYAY0TC109x4lAMd7jXkZKok5A/IxDWAaAw+PJfqHj0aDzFlBmv9TCHYAXvL9RejKOAImRqJXdK4c0bAUO/ISIzm+45EREM9GG2S5+K1EQ1twZehXDBCK+uyyKYANzDwVWOvv/Qw9I5LxGw7bpIczDysZepEAzUY28E8PZZnb7q1EY+9Lqa6Jg4Tgch9C5HgQvREeM4eRphCgiqE8WwNeUA7eFBTMI4ZP9QUyJvzGCxjmNoBGJwjCqjPgtdhut+oupau1AJXLxOw7mfFTlbkxaCGYcR5wKn/CtwkKw6HZKp77sIhpf/NVGttu3mYGRijYQMRWnKMlOt9VUjPFnIyqj/W12vZhDqo0x5iKC5WmEXUWxIli+Vun3/8LE8xOMteHmi3+CwriJPoitGiN1Fzvj4e74MuzXnE96WwoV6SZvRceJfaoTzCuPTv9BkM/Q4jUa5N7fKVd1q3isF1HHZAbOKMdgrTafgBc2h/ZS8J/LY4jwzFieyLf6YqMwqMMkgMypnpPtG14xy2eqJYUqhs4ovnwrj4tYxJMrf8kQp1KJoTfSCyMi5nSmFWt9OtXw6lw7GyAYPfnOdajukbpA2CqxPASnGvwUF6wxip1myIyGYUNTA9bXgWnRK0dOD9ntJcBT9rHivsUw3JXHxdfO/S4vyJxNpslcjerDlHce02c0BdUUbnPSX9KEWwRgQp902CuycAuyQQQFUk6SImK4hasDumEWT2mZR9s3hsSgQmRorhYs3pizybs7FNa+I6BDJbYPya5ocXDUH44Ak7UGjwK4pwHSEzdVXg6gUNTs37smYEVFtqS4K4PiVXs9QF8bj0DClt5kFr+PolMnRHIwydG21NgpUSQEZchalMeqqZookIBRG7Ofgvl2uTAEfvWhrrgI+WERHj9o6tIAuxy6bg3GMKu1Zo8C+KcDc8cGfvq2NspwfFPaHQIGvDdJFty8+3BMBiF7U6IyuyRto8kryytReS2co0ByMM8RprxoF9kgBb8hRjijK2hLGjG3QbapkRc6I/ix8ZYcECzvbOgPPC9GF9UvQpk0ZeZqcOzUH4xx12rtGgf1S4IlQc1/QrAlFPypkAd1LP4KuCbxLgoW1Ot+LH21hZ5/rOOWp28X7Ne74rjkYO2ZuQ61R4BQFvOHgQ1JMl1SV/IgZJ4ORYksLU8BHL3Au2hqDgPY+esG0UYteBHQ5d9kcjHPUae8aBfZNAaYj+LR9ddMROBk6NvclzT2Ii+jOH5q1L1QGzPROFwuQ+VYL04stTaBAczAmEKllaRTYIwWkKJlbZ6rEbV3dI44Np0aBTBTAuXjRHK951GwOxjx6tdyNAruigBQmc+y/aoTWnIxdcbYhk4sC6hu7+SO8XDSZWk9zMKZSquVrFNgpBeRk2J988fnjlhoFGgU8BeRc8EEt/mCQNUEtzaRA6GC8FzH/GRxtVDOToC17o8AWKSAFynqM+16hbhGFBnOjQFYKqC+wyBjn+1HWindUmWh04DcIvc5v4A99ICIhZf23cwAAACpJREFUoGPpdymettjnGGXas0aBRoFGgUaBRoELpoAcDDd9dIwE8h1++H+Ur1H9K93QXwAAAABJRU5ErkJggg==\n",
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- \\mu \\left(-1 + \\frac{2}{1 + e^{- \\frac{u_{1}}{v_{s}}}}\\right) \\left(- c_{c} u_{3} + 1\\right) \\left(k_{c} \\left(- \\frac{q_{3}}{2} + \\frac{\\left|{q_{3}}\\right|}{2}\\right)^{3} - q_{3}\\right)\\\\- \\mu \\left(-1 + \\frac{2}{1 + e^{- \\frac{u_{2}}{v_{s}}}}\\right) \\left(- c_{c} u_{3} + 1\\right) \\left(k_{c} \\left(- \\frac{q_{3}}{2} + \\frac{\\left|{q_{3}}\\right|}{2}\\right)^{3} - q_{3}\\right)\\\\- g m + \\left(- c_{c} u_{3} + 1\\right) \\left(k_{c} \\left(- \\frac{q_{3}}{2} + \\frac{\\left|{q_{3}}\\right|}{2}\\right)^{3} - q_{3}\\right)\\end{matrix}\\right]$"
],
"text/plain": [
"⎡ ⎛ 3 ⎞⎤\n",
"⎢ ⎛ 2 ⎞ ⎜ ⎛ q₃ │q₃│⎞ ⎟⎥\n",
"⎢-μ⋅⎜-1 + ─────────⎟⋅(-c_c⋅u₃ + 1)⋅⎜k_c⋅⎜- ── + ────⎟ - q₃⎟⎥\n",
"⎢ ⎜ -u₁ ⎟ ⎝ ⎝ 2 2 ⎠ ⎠⎥\n",
"⎢ ⎜ ────⎟ ⎥\n",
"⎢ ⎜ vₛ ⎟ ⎥\n",
"⎢ ⎝ 1 + ℯ ⎠ ⎥\n",
"⎢ ⎥\n",
"⎢ ⎛ 3 ⎞⎥\n",
"⎢ ⎛ 2 ⎞ ⎜ ⎛ q₃ │q₃│⎞ ⎟⎥\n",
"⎢-μ⋅⎜-1 + ─────────⎟⋅(-c_c⋅u₃ + 1)⋅⎜k_c⋅⎜- ── + ────⎟ - q₃⎟⎥\n",
"⎢ ⎜ -u₂ ⎟ ⎝ ⎝ 2 2 ⎠ ⎠⎥\n",
"⎢ ⎜ ────⎟ ⎥\n",
"⎢ ⎜ vₛ ⎟ ⎥\n",
"⎢ ⎝ 1 + ℯ ⎠ ⎥\n",
"⎢ ⎥\n",
"⎢ ⎛ 3 ⎞ ⎥\n",
"⎢ ⎜ ⎛ q₃ │q₃│⎞ ⎟ ⎥\n",
"⎢ -g⋅m + (-c_c⋅u₃ + 1)⋅⎜k_c⋅⎜- ── + ────⎟ - q₃⎟ ⎥\n",
"⎣ ⎝ ⎝ 2 2 ⎠ ⎠ ⎦"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fr"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAEcAAABLCAYAAAAicppkAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGAklEQVR4Ae2c320cNxCH7ww9B4oNBHnNqQM7qSBKB05SQaQO5BIMuwO7g0TuwE4Fit2BlNcgQCQhDSjftyJXPN7uelfS8nSXHWC05PDv/HY45HHuNJ/NZk/hj3ATvbu6uvqxqWAbZPP5/BQ9Fk26oPd8Jyl4S9rKKZ2lmS1Mv0Kn3Uyv78g/VzaHo+Xsgda2g6HOnYQ1CcyxlvOos+b/vHCjweEt78IXcJvPvNPr3WhwMP1LtD+HR3EHqUO+E8rragxAe2ONvdGWMxYosd8iloNP2GfAH+AF/CIMXm2XpN06T7CA19Rz5/wZ/ifI3yD/QLoiyt12D2Ct5dQ2VQF/QtmfJL9B7nK7F3JCV/CCTmdjMH2/sl/ofeCjOA55FbbwCN5P5IJ3EfOhvWAhngl2XmZ7iu+mA304btXP6MsqWMMJA0pazmNGrt+4eQugJ8hrKyGvXOAqoh/bxl3JU3vuhLXMtH3Vzj/OAfbFDKJey4qOneTvg3qezX5B2U+0OfMZ+lBBlUhJy5VeXj/qv7my5/TjKV76Cc7rf5vLBAXZISw59mCyg1GXFUo5qdpczUdG/gbWf9SyUP8Ceb38Yjkyl5Qd7iayqMPTKEufoc3KGGmdmKZuPc9HZEqRlqAl5aSy71Ihb1yZ1lrJg9XFKvajNV5GAU/rXyJr6j+pNixZEhyXwpJPSJZa7g/0KZ9Q9ow6LgeVj6SV5CAI2B9WoP5BBqbiW1ERcIKCWkIOQqU0ICyBRj39R5QdUp5a1pIjpm/7sH6U+wE6tSqKbkdFwGFqvn2XQlQ4zlZ5dLJR5lNn6+cmzzT6pJSqcxJlRzIFgv4Mfhzyv6aV75qOzmy0cw6gLDnb0nkA0roetEO+60ss3r7UsiquWDagS28wbTU4+KAFrM/yOtT0cfBLvYDa6VVrQyvh29zB4gl5sBZbbTmD0cgaTOBkgKTZCZwUjSw9gZMBkmYncFI0svQETgZImp3ASdHI0hsNDge6KaiXvdA6yyHPq4lzOF5X1GX3kdj4EzIATUG9+7CEoX0UsRx8g/cpXmUu4Cmol15oAcjagnqODfvJ3Ctan3XUIp1jTFNeRx9GtxysxpvGNKhnlKBIUI+xK2AYr3LY5I+Zi/E3r1U/S73AoVMvizYuqMectYL0Dtq76Y/o45XwZ3c4wfkC/js8eawSHbll9kI7bx3aKtbveJecX7I3xaFifZWpKPaDYvbjy6ov5pFpncryvhFdX+6bgNRD0ve1gfM1ZeIxE5x/4a/Ck8doJAh5vMnBVDYNvRh7igDUQb0IDnWbwLT+SlCPNvk2L4hSFeO6Tq78/QuJeMxKnpAfQlBPH/QiAVoMWqkIOK5xZqDZry2oxxwE5gPApJtBKzAWFAGHcQRnbUE9gDE46BY++D7ZdTj6tyyY2FoCe+imP6rOWc4h5Bu/jRHK63NOKctxUsUJi/HFu5xOSD+XSXtCb9upKLqhnZvkVqY8m+nrPPzVhIXELb2WNSW2GhxA+LJJ6b6yrV5WfUFoqzeB04YM8gmcCZwOBDqKJsuZwOlAoKNospwJnA4EOoo22nL4ODAF9dpebvgYMAX1OgDKb/vaqg6Wb/SyGqztwAZFPnjiG7xTmYJ6XhjlDDDrDur5MwBDNF5d9A7qjb6ssBovnNKgXrFf6jG2gFS/BeWFeUWq81662yHfSr2WlVsmPWxiUE/FU4ftzyC9T+5NRe6QmU19N5suO+Tx7S4tR+T39ku9OF4Y633MNz3TefayHBrcBxUP6qWTDtZv7Oz7VN6VHt3nJIOvLagHMDpk3cJLrKUp6ppM8yZZBBwmt2DItQX1AOQ1bKx/j7nkgcUbNLJUEXAYU3DWFtRLdDZMsw9A+r9eVMQhNzm/MWVorqXq1OsAHmlfkpG9g7axKas3jlKWw5hlCeWNTXmu8RnJk7r536Kg61lyt+qax1hl/gT7kGXkPw55Ams5zwJwJLtpq8EBBHem3rtTDtXWLqtc0dvkJ3A6UEvBOWVtXmXc+0NaxxgPtghdV3RmsrXO+hy/jhG/OJ0rcuv1mnf0QPOee9zyG+k/YlpTdqC7p/AAAAAASUVORK5CYII=\n",
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- m \\dot{u}_{1}\\\\- m \\dot{u}_{2}\\\\- m \\dot{u}_{3}\\end{matrix}\\right]$"
],
"text/plain": [
"⎡-m⋅u₁̇⎤\n",
"⎢ ⎥\n",
"⎢-m⋅u₂̇⎥\n",
"⎢ ⎥\n",
"⎣-m⋅u₃̇⎦"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"frstar"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"sys = System(kane)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAKIAAAAVCAYAAADFPTXWAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGW0lEQVRoBe2bi3GUNxCAfR4XYEgHpANiKgA6IEkFwR2YEjJOB9AB43SAUwETd2A6MFwHzvcpkizpf59vGHw+zQhJu9qnVivp97G6vb09aMtqtXoL7Bfqv+A/tPj9eO+BpR4gpp5Dc0q9pn4grtYlj1UZiEw+BvkP9RL4u3Livr/3wDY8QIy9gs8F9SUxdpV4toFoJnzPhFWasCtt3GQ64IT6FRt/3hXbHpodrMV714E18NQN5TB1YuvifGlgOzHE6DX1dbTvcieMerhGeCx7+ubSBmJG7HDHu8qnHbbvQZr2qAIx3k9cqH1G/MHC9VEFIr4PR7PH9A+2Do9enaNNPEBmOYfupqD1lZ1fQAX83l1kvYGJAeSz3+Id9jTe9wJgwT++2HI2jLyfAvPSfA7PrdyPY+ZV52fU9PVBOywvqJ+R9RfzvCb8TtWXwn0oZv0YLyrwU56f3J6UhMDPGL8A/msJn+pD5z3OB6xvh2t1TjTgXA99tunnPf1+V2B0kCpQX5UakmFlH5yGqsCbBKfvfWuQJs3bpIWvAf8p0dJXvh8+XbBeHYfg0OhUaYPutC6O/AxO4W+HaJfC4eUC0QTf6J+zxIN+0kP5rwq4gfotjTdpoddfBkzlG2HiWvjUONHQ6qNKN8byvJji0YeHzuAGdadnVhhEcATt83JC2QfXMQiYwZsDs5x/n748qa7mccknwhbLK/nRNwgCX1vqNoPQLJeCXX9Vm5Rx2kxVYKgDFVPvFmdpX1nUapMyToG/yGfQqWeyw88tVYCLo1Y2LNEX2kpXxgc6zl2r054NMQOnMp3AGJp/XziyvlFzNpQfRV030gE69ZdntVD31bOlh38Z4Oqas160YWiDuaEre1veY2PlUpVXBZzjCK829BivqGeeD71+y1k94s2Slawpni0eem027s7KjBjSL8DejAi8Exgt422N1YGqU6tMxdhMVu3MuTKlo2q4PKrjci6PJfOQEQKgpQHeyS7OoXQWu6UdGyd5tDmAIl/lVVl5jE+LgzZdXVq+VWC2dFNjdaLmo/2IQSgQXnIB9aLsPcOLdi7x0uqOm/UgYX71mIF3vuRmpuOdk4huL+7q1cLGOYGN+nvU+MjRzr/pX9NeOe5jsAUb1LXPXy6s8nNBljD9G+CMXfR1njCvo7wvPXRZ3pb5/tSq1fhMXSo703zmeQ0x4eW/rOSMCNB5vZfdiOtkKOERl1M0PMw84TiiNbNVl9xEM9ZCM5RN1CHLGuNR4qCp7l+Mq2OMcXXXYbwNG8xwLd8ktz2uc9aCJt/NShum+lHn6miPvLLPGC/OYtB4iuTMpR4U16c6ORkbO+W6V7qU+se5VVwcApxbfKZ3MiXRrRPDZ4/UR2jIMrRmhJelAOYcU2+ppuahEuiZ46KEQl+HWKoMBnwOvyqTolfKNp3PNfAL9gzZMEcec9TboPMKUBazk5u3sgGQJ0CCmbVzZoTXlK/KjJ/9pRxK+nR0BR+TQrZ3jh2Bw/80ma90wP0U1GZ7+b8Db5AqZ+pT0Vfm5HKUexMdBJ8i5NzKVDOGxSPEn/ak8hudP9LAtlWY8RoeKmoA9R5BcY7Bo7zPtB4D0nSOnjn8oNN5BlhZXCRt0pY/C8SoDTPluXDqmoIrsRfe991N+a/Rxcyd9ZwpS94hwGm9duRrEfTJPtfohnG+Is3lzTyD6yLyDesuTKFNUYayU8Lw22Pa8M3UniFM8/EMOgRZCVvSh968PfjyLnkxz51TXYBLfNtnrtll8MW7lF/LP43hM8uGbclLcsfaKVngDd6NHiRTvMf0SrjII687Y4Nx8BoArhNnhwDb8rQFLBibrnMal46dlHd4w8f03rtjoDHN58LYjOauT7st44rOIL9izpzuXBu2JW+OTlOyzOIf5zDqmTPFu4ekAzKwXKNUPHbb0yDhetv2aDb1lgx7iUaA3gtM5QaSQSavThqPgXUDrlPAadQZ7ZMiUA1Af9Xba9wYv46AacCkDVuWN6rRlCzwbnz93OubMeZTvMdoG5xrfAI/77rq4uPGDT1UnLOukBDko9k+xUdE9dpr59x3DP/q+2DJD1zKfN6X/OZnRh19KYMf5Ffy3lb/e8qbkgVef236bfW7+k3/R339olB9Oah+oe0sC5Ht/c3HgveOvsu10/Zl74HZHiCmPCV90Hjqdv7Pyn/FUX1UBp5i7QAAAABJRU5ErkJggg==\n",
"text/latex": [
"$\\displaystyle \\left\\{c_{c}, g, k_{c}, m, \\mu, v_{s}\\right\\}$"
],
"text/plain": [
"{c_c, g, k_c, m, μ, vₛ}"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sys.constants_symbols"
]
},
{
"cell_type": "code",
"execution_count": 134,
"metadata": {},
"outputs": [],
"source": [
"sys.constants = {cc: 0.1, #0.85,\n",
" g: 9.81,\n",
" kc: 5.0e7,\n",
" m: 1.0,\n",
" mu: 0.2, #1.0,\n",
" vs: 0.01}"
]
},
{
"cell_type": "code",
"execution_count": 176,
"metadata": {},
"outputs": [],
"source": [
"sys.times = np.linspace(0.0, 5.0, num=5000)"
]
},
{
"cell_type": "code",
"execution_count": 177,
"metadata": {},
"outputs": [],
"source": [
"sys.initial_conditions = {q1: 0.0,\n",
" q2: 0.0,\n",
" q3: 5.0,\n",
" u1: 5.0,\n",
" u2: 10.0,\n",
" u3: 0.0}"
]
},
{
"cell_type": "code",
"execution_count": 178,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<string>:2: RuntimeWarning: overflow encountered in exp\n"
]
}
],
"source": [
"x = sys.integrate()"
]
},
{
"cell_type": "code",
"execution_count": 179,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3xUVfrH8c9JT0ggkISEkAQCofcu0quCdEUXFVFAdAW7a1l3V3ftu65lf2LBioJKERABRREQQYr0GkgIJaGkEUJ6ppzfHzdUExKSmdxM8rxfr7wmmbkz8wzKl5Nzzz2P0lojhBDC9biZXYAQQojykQAXQggXJQEuhBAuSgJcCCFclAS4EEK4KI/KfLPg4GDduHHjynxLIYRwedu2bUvTWodceX+lBnjjxo3ZunVrZb6lEEK4PKXUseLulykUIYRwURLgQgjhoiTAhRDCRVXqHHhxLBYLSUlJ5Ofnm12K0/n4+BAREYGnp6fZpQghqgHTAzwpKYmAgAAaN26MUsrscpxGa016ejpJSUlER0ebXY4QohowfQolPz+foKCgah3eAEopgoKCasRvGkKIymF6gAPVPrzPqymfUwhROUyfQhFCCJdlt0FhNhRkX3574fusotsc6DgB6jVx6NtLgAshxHn55yA9DtLiIO0QZKeUEMxFt9a8sr92ZA8JcCGEqBCt4dwJI6DPB/X577NOXTxOuYN/ffDyB29/47Z2xMXvvWqBd8Dlj5f0vacfuDl+xloCvBzy8vK48cYbWb16Ne7u7iQlJbFhwwbGjh3L4MGDWb16NR4e8kcrhKks+XDmcDFBHQ+WnIvHedeB4GbQZIBxG9zc+KrbGDy8TCu/LCRlyuGTTz5h3LhxuLu7A/Dzzz+zf/9+brvtNgYNGsS8efO44447TK5SiBpAa8hJu3wUff77s8eB8y0jFQRGGsHcqNfFoA5qZoyyXXSBgQR4kbi4OO69916ys7MZOHAgixYtIj4+vthj586dy5dffgnA+vXreeyxxwgMDGTlypXMnj2bZ555RgJcCEeyWSHjaPFBnX/24nEevhAcAxFdoePtF4O6XlPw8jOtfGepUgH+z+/2sf/kOYe+Zuvw2jw3ss1Vj7HZbNx1113MnDmTzp078+CDD9KmTfHPKSwsJCEhgfPb4vbu3Ztu3brx+uuv07ZtW2w2G7///rtDP4MQNY7dDid3QPwqiP8JTu4Eu+Xi4/5hRji3vbloyiPGuK0d4ZS55qqqSgW4WZYsWULr1q3p3LkzAK1atSIwMJCEhAReeuklMjMzWbhwIQBpaWkEBgZe9vyDBw/SokULANzd3fHy8iIrK4uAgIDK/SBCuLKcdDj8M8T9ZNzmpgMKGnaBntMhpOXFsPapY3a1VUKVCvDSRsrOsmPHDjp27Hjh5127djFkyBCaNGnCxx9/zC233HLhMV9f38uupkxPT6dOnTqX7W9SUFCAj49P5RQvhKuy2+DEdmOEHb/K+B4NfsEQMxhihkDTgVAryOxKq6xSA1wpFQl8DoQBdmCW1vptpdTzwL1AatGhf9Var3BWoc4UFBREbGwsAJs3b+bzzz/n8ccfL/bYunXrYrPZyM/Px8fHhyNHjhAeHn7h8fT0dEJCQmTDKiGKk516ySh7NeSdAeUGDbtC/2eg2WBo0KlGTYNURFlG4Fbgca31dqVUALBNKfVT0WNvaq1fd155lWPixIncdNNNtGvXjuHDhxMUFERMTEyJxw8dOpT169czePBgWrZsSVpaGm3btmXWrFmcPHmS4cOHV2L1QlRhdhskbb04yj65w7i/Vgg0v8EYaTcdCH71zK3TRZUa4FrrU8Cpou+zlFIHgIbOLqwyBQcHs3nzZgASExNZu3Ytbm5upKen8+yzz7Jjxw5eeeUVnnnmGQBmzJjBG2+8weDBg/H392fLli0XXmvcuHG88sorpnwOIaqE7BQjrM+PsvPPGqPsiG4w4G/GKDusg4yyHeCa5sCVUo2BTsBmoBcwQyl1F7AVY5SeUcxzpgHTAKKioipYrvPt2rWL9u3bA8bUyvvvv/+HYzp16sSAAQOw2WwX1oKDsUJlzJgxF05oClEj2KyQ9PvFUfapXcb9/qHQ8iZjlN2kv4yynUBprUs/ClBK+QO/AC9prRcppUKBNIyV8i8ADbTWk6/2Gl27dtVXNjU+cOAArVq1Kk/tLqmmfV5RTWWdvjjKTlgD+ZnGpeeR3Y3AbjYEQtvJKNtBlFLbtNZdr7y/TCNwpZQn8A0wV2u9CEBrnXzJ4x8CyxxUqxCiqrFZIHHLxVH26T3G/f5h0Gpk0Sh7APgGXv11hEOVZRWKAj4GDmit37jk/gZF8+MAY4G9zilRCGGKcycvGWWvhYJzxig76joY9FzRKLuty16GXh2UZQTeC5gI7FFK7Sy676/ABKVUR4wplKPAfU6pUAhRebKSYccXsG8xJBeNyQLCofVoI7Cb9JeLaKqQsqxCWQ8U90+sS675FkJcQWs4vgk2vw+xy8BuhaieMPifRmjXby2j7CqqSl2JKYSoRDYL7P8WNs6Ek9vBJxB63A9d7jb2GRFVngS4EDVN3lnYPhs2f2A0NqjXFG76L3SYYDQpEC5DAlyImuJMAmx6H3bMMRoaNO5jBHezG2S5n4uSABeiOtMajm80pklil4Obh7EFa88HoEEHs6sTFST/7JZDXl4e/fr1w2azAZCUlMS8efMoLCykb9++WK1WkysUNZ7NArsXwIcD4NNhcGwD9HkMHtkD4z6Q8K4mJMDLobiWatu3b8fLy+tCSzUhTJGXAevfhLfaw6KpUJAFN70Bj+6HQf+A2g3MrlA4kEyhFJGWasKlpR+GTe/BzrlgyYXovjDyLWNPbZnfrraqVoB///TFS3QdJawdDHv1qodISzXhkrQ2pkY2vgsHVxjz2+3GG/PbYe3Mrk5UgqoV4CYpqaXakiVLWL58OSkpKUyfPp2hQ4dKSzVhPmuhcaXkppnGzn++9aDvE9BtKgSEmV2dqERVK8BLGSk7S0kt1caMGcOYMWPIyMjgiSeeYOjQodJSTZgn9wxs+xS2fAhZp4z+kCPegva3VcuO66J0VSvATVJaS7UXX3yR6dOnA9JSTZggLR42vwc7vzTmt5v0h1H/B00Hyfx2DScBTskt1bTWPP300wwbNuzC9ApISzVRCbSGo78a89uHfgB3T2h3K1z3Zwhra3Z1ooqQAKfklmr/+9//WLVqFZmZmcTHx3P//fcD0lJNOJG1EPYtgo3vGCf0/YKg35PQdQoEhJpdnahiJMCvcGlLtYceeoiHHnroD8dISzXhcLlnYOsnxvx29mkIbgEj3zbmtz19za5OVFFlbqnmCNJSreZ9XlGKtDjY9C7s/AqseUZXm54zjE7tMr8tilSopZoQwoG0hiPrjP1J4laCuxe0vxWuewBCi7/+QIjiSIALUVmsBbD3GyO4k/eCXzD0exq6TQH/+mZXJ1yQBLgQzqY17PoaVj0H2ckQ0tJYBtjuVvCU6wVE+UmAC+FM6Ydh2aNw5BeI6AZj3jPmt6VFmXAACXAhnMFmgd/+D355zZjjvum/0GWynJgUDiUBLoSjJW2F7x425rlbjYRh/4ba4aU/T4hrJAEuhKMUZMHPL8CWWRDQAG6bC61GmF2VqMYkwIVwhNgVsOIJOHcSut8LA/8OPrXNrkpUczIhVw7SUk1ckHUa5k2EryeATx2Y8hMM/4+Et6gUEuDlIC3VBHY7/P4xvNMdDq002pXdtw4iu5ldmahBJMCLxMXF0b9/f7p27cqTTz5JTExMicfOnTuX0aNHAxdbqi1cuJCOHTsyZswY5s6dW1llCzOkxBqNgpc/BuEd4IGN0OdxY8dAISpRlZoDf23La8SeiXXoa7as15Knuj911WOkpZooE0s+/Ppfo2mwt7+xprvDBFnTLUxTpQLcLCW1VDtw4ABvv/02aWlpDBo0iD//+c/SUq2mOroevnsE0uOMHQJveBlqBZtdlajhSg1wpVQk8DkQBtiBWVrrt5VS9YB5QGPgKHCr1jqjIsWUNlJ2lpJaqrVq1Yr3338fu93OvffeCyAt1WqavAz46R+w/XMIbAR3LoKYQWZXJQRQtjlwK/C41roVcB0wXSnVGnga+Flr3Qz4uehnl1RcS7Xze4IvXbqU3r17M2iQ8Zf20pZqgLRUq660Njaeeqcb7JgLvR6GBzZJeIsqpdQA11qf0lpvL/o+CzgANARGA7OLDpsNjHFWkc42ceJEtm7dSrt27Vi0aNGFlmoAo0aN4rfffrvsxOT5lmrAZS3VfvvtN9asWSMt1Vzd2ePw5a2wcDLUiYBpa2HIv6RxsKhyrmkOXCnVGOgEbAZCtdanwAh5pVSx+2EqpaYB0wCioqIqUqvTlNRSbe3atSxatIiCgoLLQllaqlVTNits+QBWvwgouPFV6D4N3NxLfaoQZihzgCul/IFvgEe01udUGc+8a61nAbPA6MhTniIr06Ut1fr370///v3/cIy0VKuGTu2CpQ/BqZ3Q7AZj86nASLOrEuKqyhTgSilPjPCeq7VeVHR3slKqQdHouwGQ4qwiK9OIESMYMaL0/SsmT578h/u8vLy46667nFGWcJbCHFj7itH93S8IbvkU2oyVpYHCJZRlFYoCPgYOaK3fuOShpcAk4NWi22+dUqEQzhK/ytir++xx6DwJhvwTfOuaXZUQZVaWEXgvYCKwRym1s+i+v2IE93yl1BTgODDeOSUK4WDZqbDyGdizAIKbwz3fQ6Prza5KiGtWaoBrrdcDJf0+6ZA1VVpryjqn7sq0rvKnAKo3rWHnXFj5rDF10u9p6PMYeHibXZkQ5WL6lZg+Pj6kp6cTFBRUrUNca016erpc4GOWtHhY9ggc/RWiesLItyFETjYL12Z6gEdERJCUlERqaqrZpTidj48PERERZpdRs1gL4be34Zf/gIcPjHjLmO+W1maiGjA9wD09PYmOjja7DFEdJW4xWpul7IfWY2DYaxAQZnZVQjiM6QEuhMPlZ8LP/zL2667dECZ8DS2GmV2VEA4nAS6qlwPLjNZmWaehx/0w8Fnwll0hRfUkAS6qh3MnYcVfIHYZhLYzGgpHdDG7KiGcSgJcuDa7HbZ+DKv+CXYLDP4n9Jwu3XFEjSABLlxX8n7jJGXSFmjSH0a8CfWamF2VEJVGAly4Hks+rPsPbHjL6AQ/9gOjS041vo5AiOJIgAvXcmSd0drszGGjH+XQl6BWkNlVCWEKCXDhGnLPwI9/h51zoG40TFwCTQeYXZUQppIAF1Wb1rBnIfzwtNGfsvej0O8p8PQ1uzIhTCcBLqqujKOw7DE4/DM07AJ3fQthbc2uSogqQwJcVD02K2x6F9a8bLQzG/Zv6DZVWpsJcQUJcFG1JO+DxffD6d3QfBjc9LrRWFgI8QcS4KLqOLwG5k00ur/f+jm0GiVLA4W4CglwUTXsng9LHjA65Ny5EGqHm12REFWebIoszKU1bHgbFt0LUdfBPSskvIUoIxmBC/PY7bDyr7D5PaMT/NgPpL2ZENdAAlyYw5IPS+6HfYvhugeMKyqlS44Q10QCXFS+vLPw9R1wbD0MfRGuf9DsioRwSRLgonKdOwlzboa0OBj3EbQfb3ZFQrgsCXBReVIOwJxbjJZndy40toAVQpSbBLioHMc2wle3GZ3h71kBDdqbXZEQLk/OGgnn278UPh8NterDlJ8kvIVwEAlw4VxbPoT5d0GDDjDlR6jbyOyKhKg2ZApFOIfW8PO/YP0b0GI43PyxcYm8EMJhJMCF49kssPRB2PUVdLkHhr8O7vK/mhCOJn+rhGMVZMH8ScYe3gOehb5/kQ2phHCSUufAlVKfKKVSlFJ7L7nveaXUCaXUzqKv4c4tU7iE7BT4bAQkrIVR70C/JyW8hXCispzE/Ay4sZj739Radyz6WuHYsoTLST8MHw+BtEMw4SvoPNHsioSo9kqdQtFar1NKNXZ+KcJlJW2DL4uuqJy0DCK6mFuPEDVERZYRzlBK7S6aYqlb0kFKqWlKqa1Kqa2pqakVeDtRJR36EWaPAO8AY423hLcQlaa8Af4e0BToCJwC/lvSgVrrWVrrrlrrriEhIeV8O1Elbf8CvvoTBDczwjuoqdkVCVGjlCvAtdbJWmub1toOfAh0d2xZokrTGta+BktnGPuZ3L0C/OubXZUQNU65lhEqpRporU8V/TgW2Hu140U1YrPCisdh22fQ4XYY9T9w9zS7KiFqpFIDXCn1FdAfCFZKJQHPAf2VUh0BDRwF7nNijaKqKMyFhZPh0PfQ53EY+HdZJiiEicqyCmVCMXd/7IRaRFWWk27sJpi01biysvu9ZlckRI0nV2KK0mUcNZowZCbBbV9Aq5FmVySEQAJclObULpg7HqwFcNe3Rud4IUSVINvJipIdXg2fDgd3L2MrWAlvIaoUCXBRvF3zjJF33cbGGu+QFmZXJIS4ggS4uJzWsP5NWDwNonoa7c9qNzC7KiFEMWQOXFxkt8EPT8OWWdD2FhjzLnh4m12VEKIEEuDCYMmHRffCgaXQcwYMeQHc5Bc0IaoyCXABeRnw1e1w/De44WXoOd3sioQQZSABXtOdTYS5t8CZBLjlE2h7s9kVCSHKSAK8JkveB3NugcJsuPMbiO5rdkVCiGsgAV5THfkVvr4dvGrB5B8gtI3ZFQkhrpGcpaqJ9iyEOeOgdrixxlvCWwiXJAFek2gNv74B30yBiG7GyDsw0uyqhBDlJFMoNYXNCiuegG2fyhpvIaoJCfCaoCAbFt4DcT9C70dh4D9kjbcQ1YAEeHWXlWx0jD+9B0a8CV0nm12REMJBJMCrs5RYY0Oq3HSY8DU0v8HsioQQDiQBXl0dXW8sE3T3hnuWQ3gnsysSQjiYTIRWR7sXwBdjwT8Mpq6S8BaimpIAr07OLxNcNBUiusOUlVC3kdlVCSGcRKZQqgubFVY8Dts+k2WCQtQQEuDVwWXLBB+DgX+XZYJC1AAS4K4u6zR8eassExSiBpIAd2UpscZWsLlnYMI8aD7U7IqEEJVIAtxVHfkVvr4DPH1kmaAQNZRMlLqi88sEA2SZoBA1mQS4K9Eafv2vsUwwsoexTDAwyuyqhBAmkSkUV3HpMsF242H0TFkmKEQNV+oIXCn1iVIqRSm195L76imlflJKxRXd1nVumTVcQTZ8PcEI7z6Pw9hZEt5CiDJNoXwG3HjFfU8DP2utmwE/F/0snCHrNHw2HOJ/hhFvwSDZClYIYSg1CbTW64AzV9w9Gphd9P1sYIyD6xJgLBP8aDCkxRu7CXa9x+yKhBBVSHmHcqFa61MARbf1SzpQKTVNKbVVKbU1NTW1nG9XAx35FT4eCrZCuGeFrPEWQvyB038X11rP0lp31Vp3DQkJcfbbVQ+751+xTLCj2RUJIaqg8gZ4slKqAUDRbYrjSqrBtIZ1r8OieyHqOlkmKIS4qvIG+FJgUtH3k4BvHVNODWazwncPw+oXjGWCd34DvrK4RwhRslLXgSulvgL6A8FKqSTgOeBVYL5SagpwHBjvzCKrvYIsWHAPxP9kLBMc8DdZaSKEKFWpAa61nlDCQ4McXEvNlHXa6FuZvM9YJigrTYQQZSRXYpop5UBR0+EzRU2HZaWJEKLsJMDNcmQdfH1n0W6CK2SliRDimslEqxl2z4cvxkHtBrJMUAhRbhLglenKZYKTf5BlgkKIcpMplMr0499g4zuym6AQwiEkwCvL7vlGeHe7F4b/B5QyuyIhhIuTKZTKkHLAuEgn6nq48VUJbyGEQ0iAO1tBFsybCF7+MP5TcJdfeoQQjiFp4kxaw7cz4EwCTFpqbE4lhBAOIiNwZ9r0HuxfYjRhaNzb7GqEENWMBLizHN8EP/0dWtwEvR42uxohRDUkAe4M2amw4G6oEwlj3pWTlkIIp5A5cEez2+CbyZCXAVN+At9AsysSQlRTEuCOtuYlY5+TUe9Ag/ZmVyMqidYam11jtWsKbXasNo3FZi/60lht9ivuN26tdjs2u9nVi8rQKSqQYH/HXrwnAe5IB3+AX/8LnSZC54lmV+OS8gptZBdYsdovhp3VXnRr01jt5wNRYyk6xmqzY7EX3V4SmMbzLn+8+Ne65HGbvuz7wqKQtViN97NcEcLW82Fst6O12X96oir77J5u9G9RYvvgcpEAd5SMo7B4GoS1M660FGWitWb/qXOsPZjKLwdT2XY8A5vdOUno4abwcFd4urnh4a7wcHfD0824Lf5+RYCnB57ubngW3e/l7oaHm8LTwzjG092t6P6i5xUd63n+NS/92c0NLw+FR9H7eBU9110pOU1SA0QF+Tn8NSXAHcGSD/PvAg3c+gV4+ppdUZWWmWthXVwqvxwyvlKzCgBo3aA29/ZpQnigzxUhZ4Se51VC19P9YjB6FoWsx/ngLXqukpQU1YwEuCN8/ySc2gV/+grqRZtdTZVjt2v2nsw0RtmHUtlxPAO7hto+HvRpHkL/5iH0ax5C/do+ZpcqhEuRAK+onV/C9tnQ6xFoOdzsaqqMjJxC1sWlsvZgKusOpZKeUwhA+4g6TB8QQ/8WIXSICMTDXVayClFeEuAVcXovLHsUGveBgX83uxpT2e2aPSeMUfbaQynsSjyLXUNdP0/6Fo2w+zYPcfhZeCFqMgnw8srPhPkTwScQbv64Rm5SdX6U/UvR1Eh6TiFKQfuGdZgxsBkDWoTQPiIQdzeZexbCGWpe6jiC1vDtdMg4Bncvg4BQsyuqFOfnstfEFj/K7t8ihL7NQgiSUbYQlUICvDw2vgMHvoOhL0Kj682uxqkuHWWvi0slLVtG2UJUFRLg1+rYb/DTc9BqJPScYXY1DldgtbHtWAYb4tPYEJ/O7iQZZQtRVUmAX4usZGOTqrqNjJ6W1WBdsd1uXEizPj6NDfFp/H70DPkWO+5uio6RgTLKFqIKkwAvK5sVFk6G/HNw5yLwqWN2ReWiteb4mdwLgb3xcDoZuRYAmof6M6F7FL1jgukeXY8AH0+TqxVCXI0EeFmtfgGOrYcx70FYW7OruSZp2QX8djidDXFprI9P48TZPAAa1PFhUKtQescEc33TILmQRggXIwFeFrErYMNb0HkSdLzd7GpKlVNgZcvRMxcCO/Z0FmBc+dizaRD392tCr5hgooNryeXlQriwCgW4UuookAXYAKvWuqsjiqpSziTA4vuhQQcY9m+zqymWxWZnd9JZ1selsyE+jR2JGVhsGi8PN7o2qstfbmhB75hg2jasI/PYQlQjjhiBD9BapzngdaoeS56xSZUCbv0cPKvGFIPWmriUbNbHGfPYm4+cIbvAilLQNrwOU3o3oXdMMF0b18XH093scoUQTuISUyi/nfiNQxmH8PXwxcfDh3o+9Wjo35AG/g3w9XDizn8r/gKn98CEeVC3sfPepwxOns0rWtqXxobD6Rd28IsOrsXojuH0jgmmZ9MgAv28TK1TCFF5KhrgGvhRKaWBD7TWs648QCk1DZgGEBUVVa43WXV8FQsOLSj2sca1G9MmuA0dQjrQO7w3kbUjy/Uef7BjDuz4Avo8Di1udMxrXoPMXAsbE9IvhHZCWg4Awf5e9IoJplfTYK6PCSKiruP3GBZCuAalK9BGRCkVrrU+qZSqD/wEPKi1XlfS8V27dtVbt2695vex2q0U2ArIs+aRZ8kjPT+dE9knOJ51nAPpB9iXto+UvBTACPRBUYMYHTOa6Drl3Nr11G74eAhEdoeJS8DN+dMQ+RYb249lXFjet+dEJnYNfl7uXNckiOubBtG7WTAtQgPkxKMQNYxSaltx5xgrFOBXvMHzQLbW+vWSjilvgJdF4rlE1p1Yxy+Jv7D59Gbs2k6HkA7c0eoOhjQagodbGX/ZyDsLs/qDNR/u+xX8Q5xSr82u2X/y8gtoCqx2PNwUnaICub5pML2bBdMhIhAvD9lyVYiazOEBrpSqBbhprbOKvv8J+JfW+oeSnuPMAL9Uam4qyxKWsShuEUfPHaWhf0PubnM3Nze7GU/3q1ycojXMuxMO/QB3L4eo6xxaV+KZXH45lMqG+DR+O5xOZp5xAU3LsABjWiQmiO7RQfh7u8SpCSFEJXFGgDcBFhf96AF8qbV+6WrPqawAP8+u7axJXMMnez9hd+puIgMiebTLowyOGlz8NMSGt+Gnf8ANL0PP6Q6rw2KzM3NNPP+3Oh6bXRNex4fezYLpVXTisX5A1VjdIoSompw+hVIWlR3g52mtWX9iPW9se4P4s/F0rt+Z53o+R5PAJhcPOroBZo+EViNg/GyH7XMSn5LFY/N3sTspkzEdw3loUDO5gEYIcU1qdICfZ7VbWRK/hLe2v0WuJZdp7acxpe0UPHPT4YO+4B0A964Bn9oVfi+7XfPpb0f59w+x+Hm589LYdgxv18ABn0IIUdOUFOA1arLVw82DW5rfQv/I/ry25TVm7pzJ6mOr+E9qBo0KsowVJw4I76SMXJ5YsItNCWcY1LI+r9zcTqZJhBAOVyOXNwT7BvOffv/hrQFvcfJsAre6nWZZ72kQ2rpCr6u1ZtH2JIa99St7kjL5983t+WhSVwlvIYRT1KgR+JUGZaTS5thRnmrajmeOLmKrp+LZHs9efaVKCc7mFvLs4r0s33OKro3q8uZtHYmsJxfZCCGcp+YG+Knd8N3DhEVez8fjFzJzzwd8tOcjjp07xhv936CuT90yv9Svcak8sWAX6dmF/OWGFtzfr6lsGiWEcLoaOYVC7hljvbdvPbjlUzw8fXi488O82udVdqfu5vblt3P47OFSXybfYuP5pfuY+PEWAnw8WTK9F9MHxEh4CyEqRc0LcLsNFt0LWafgti8uu9LypiY38emNn5JnzWPSD5PYnbq7xJfZeyKTEf+3ns9+O8rd1zdm2YO9advQNbv0CCFcU80L8LWvQvwqGPYaRPxx+/L2Ie35YvgXBHgGMPXHqWw8ufGyx212zcw18Yx9dwPn8izMntyd50e1kW1bhRCVrmYFeOwKWPdv6HQndLmnxMMiAyL5fNjnRAREMP3n6aw6tgowLoX/06yN/GflQYa0DmXlI33p19w5e6UIIURpas6FPGnx8OEACGoK9/xQpuYMmQWZTP95OvvS9nFzxDN8/YsxRfLPUW0Y17mhXE0phM9otkkAABLlSURBVKgUJV3IUzNG4AXZxklLN49r6qxTx7sOL/d8Gx8dxdfHXyaq4TG+f7gPN3eJkPAWQpiu+ge41rB0BqQdhFs+gcCyN5X4Ye8pxr6znfT4SdT3bkyy7wck5u10YrFCCFF21T/AN86EfYth0HPQdECZnnI2t5CHv97B/XO2E1bHh++mD2HxuNlE14nm4TUPsyt1l5OLFkKI0lXvAD+yztgettUo6PVwmZ6yOjaZoW+uY/nuUzwyuBlLpveiRVgAdbzr8MGQDwjxC2HGzzM4mnnUubULIUQpqm+AZybBgnuMk5Zj3i11e9hz+RaeXLiLyZ9tpa6fF0um9+KRwc3xdL/4RxTkG8T7g9/HTblx/6r7SctLc/anEEKIElXPALfkw/y7wFoAt801tom9ivVxadz45joWbkvigf5NWfpgrxIvyomqHcXMQTM5k3+GB1Y9QI4lxxmfwOlsdhsFtgKyC7PJteRSmauRhBCOUf32QtEalj8OJ7bBbXMgpHmJh+YUWHnl+wPM2XScJiG1+ObP19MpqvQ9UNoGt+X1fq/z0OqHeObXZ3hrwFu4qar1b6HWmsSsRA5lHCL+bDwJZxM4lXOKtLw00vPTybPmXXa8m3LDz8OPQO9AGvo3JNw/nKjaUbSu15rWQa0J9Ak06ZMIIUpS/QJ8y4ewcw70ewpajSzxsI2H03nqm90kZuQytXc0T9zQ4pqupuwb0Ze/dPsLr255lXd2vMNDnR9yRPUVkpiVyLqkdWxL3sb25O2k56dfeKyhf0Mi/CNoH9KeYN9g/L388XTzxNPNE5u2kWPJIceSw5m8M5zMOcn6E+tJjU+98PwI/wh6hvekZ3hPejToQW2viu+bLoSomOoV4Ed+hR+ehhbDod/TxR6SlW/hle9j+XLzcRoF+TFvWk+6R9cr19vd3vJ24jLi+HDPhzSr24xh0cMqUn25JGQmsCJhBasTVxOXEQdAeK1weob3pHNoZ1rXa010nWj8PK99a9vMgkxiz8SyP30/21O2s+LIChYcWoCH8uD6htczLHoYAyIHUMuzlqM/lhCiDKrPlZhnj8Os/uAXDFNXFdtZZ3VsMs8u3kvyuXym9I7msSEt8PWq2B4mFpuFqT9OZV/6PmYPm02boDYVer2yyLHksPLoShbHLWZn6k7clBud6ndiQOQABkYOJLJ2pFPe12K3sCd1D2sS1/DD0R84nXMaXw9fbmpyE39q8Sda1GvhlPcVoqar3j0xC3Phkxsg4xjcuxqCYy57+ExOIf/6bh9Ldp6kRWgAr93Sno6RjpvTTc9LZ8LyCdi1nQUjF1zTXuLXIiU3hbkH5rLg4AKyLFlE14lmbMxYRjYdSbBvsFPesyR2bWdnyk6WxC9hxZEVFNgK6BLahantptIrvJdcqSqEA1XfANcavpkKe7+B2+dD86GXPKRZtvsUzy/dx7l8Cw/0j2H6gBi8PBx/wnF/+n4mrphIt7BuvDv4XYee1EzKSmLW7ll8l/Addm1nUNQgJraeSMeQjlUiKDMLMlkct5gvY7/kVM4p2gW34/4O99OnYZ8qUZ8Qrq76BviGt42LdQb9A/o8fuHu5HP5PLt4L6sOJNMhog6v3dKelmHOPfE2/+B8Xtj0AjM6zuC+DvdV+PXS8tKYtXsWCw4twF25MyZmDJNaT3LaFElFWWwWvj38LR/u/pCTOSfpXL8zT3V/itZBFes1KkRNVz0DPH4VzB1vXGk5/jNQCq0187cm8uLyAxRa7TwxtAWTe0dXSpccrTXPrH+GFQkrmDV0Ftc1uK5cr1NgK+CTvZ/w6d5PKbQVMrbZWO5vfz+htUIdXLFzWGwWFscvZubOmWTkZzA6ZjQPd3640qd5ciw5HM08SkJmAsm5yaTnpZOen05WYRZWuxWbtmHXdnw9fPH39KeWZy1C/UIJ9w8n3D+cJnWaEOQbVKk1C1Gc6hfg6YeN7WHrRMKUH8GrFgmp2Ty7eC8bE9LpEV2P125uT+Pgyl0hkWvJZcLyCZwtOMuCkQuo71f/mp6/Lmkdr2x+haTsJIY0GsJDnR6icZ3GzinWybIKs5i1exZzDszB192Xx7s+zrhm45wyrWKz24jNiGVnyk52puxkV+ouTuWcuuyYWp61CPIJIsArAE83T9zd3FEo8qx55FhyOFd4jjP5Zy57Tn2/+rSq14r2Ie3pHtadtsFt8XCrXou3RNVXvQK8IAs+GgLZp2HaWgoCInl/bQIz18Tj7enGM8Na8adukbiZ1Jvy8NnDTFg+gbbBbflwyIe4u5W+0iUlN4WXNr3E6sTVNK7dmL/2+Cs9w3tWQrXOd+zcMZ7/7Xm2Jm+lR1gPnuv5nEOmgXIsOaw/sZ51SetYf2L9hfANqxVGh5AOtKzXkuja0UTXiSbcPxwfj9K3ES60FXI65zRJ2UnEZcQReyaWA+kHSMhMQKOp5VmLHmE9GNxoMP0j+xPgdfWrfIVwhOoT4HY7zJ8IB7+HiYvYRDv+ungPCak5jOoQzt9GtKJ+QNn2+3amxXGL+cdv/+CRzo8wpd2UEo/TWrP8yHJe3vwyFpuF+zrcx6TWk/B096zEap3Pru18E/cNb2x9A6vdylPdn+LmZjdf82jcarey+dRmlh5eyurjq8m35VPbqzZ9IvrQp2EfuoR2IaxWmMPrz8jPYMvpLWw+tZlfkn4hJTcFTzdPeoX3YkyzMfSL6OeUkXmuJZcj545wOvs0qXmppOalklWYRaGtkAJbATa7DW8Pb7zdvfHz8CPIN4hQv1Dq+9WnUe1GMgVUTVSfAP/l37DmJXIGvMBzKf1YuC2JyHq+vDimXZVqb6a15olfnmD18dXMGT6HNsF/XB+enpfOi5teZNXxVXQI6cBLvV+iUe1GJlRbeU7nnObvG/7OplObGNJoCM/1fI463qU3g87Iz2DBoQXMi51HSl4KAV4BDGs8jGHRw+hYv2OlTmvYtZ3dqbtZeXQlK4+uJDUvlfp+9RnXbBzjm4+/5mmz8/Kt+exN28uOlB3sSt1F/Nl4TmSfuOwYN+VGgFcA3m7eeLl74e7mToGtgHxrPrmWXArthZcdX8+nHs0Cm9EqqBVdQrvQqX6nMv15l5fWmgJbAQW2Auzajre7UadMO1WMUwJcKXUj8DbgDnyktX71asdXOMBjV8DXEzgWMYoxJ+4kq8DGtL5NeHBgswpfkOMMmQWZ3Lz0Znw8fJg/Yv5lV0NuOrWJp9c9zbnCczzY6UHuan1XmaZaqgO7tjN732z+t/1/BPsF82qfV+kS2qXYY49kHmH2vtksS1hGga2A68OvZ3zz8fSN6IuXu1clV/5HVruVdUnrWHBoARtObMDDzYNRTUcxue1komqX3jwk8Vwia5PWsjZxLdtTtmO1WwFoUqcJLeq2oGlgU5oGNiXcP5wQ3xDq+dQr8f8TrTVZliySc5JJzk3mSOYR4jLiiMuI41DGIQrthSgUzes2p19kPwZGDaR1vdbX9FuQ1poT2Sc4fPYwcWfjOJJ55ML7peallri5m6ebJ0G+QdT3rU+IXwhRtaOICYyhaR3j85Vleqs0VruV7MJs8qx55FpzsdgteLp54uHmgbe7N4HegQ55HzM4PMCVUu7AIWAIkAT8DkzQWu8v6TkVCvDUg9hnDeQI4QzP+ittG4Xy8th2tAir2nOQv5/+nSkrpzCu2Tiev/55bHYbs3bP4r1d7xFdJ5rX+71Os7rNzC7TFPvS9vHkuidJyk5iWvtp3Nf+vgsjtWPnjvHBrg9YfmQ5Xm5ejGw6kjta3UHTwKYmV12yxKxEZu+bzeK4xVi1lRsa3cDU9lNpXvfyDdVSclNYlrCM7w5/R/zZeABiAmMuTAF1rN/R4aPkAlsBe9P2si15GxtPbmR7ynbs2k6oXyjDoocxJmZMsX+2WmsOnDnA76d/Z3vydnak7CCjIOPC46F+oTSo1YAQvxBCfEOo7V0bb3dvfNx9UEphsVkosBWQY80hPS+d1NxUknOTScxKxGK3AODh5kGboDZ0rt+ZLqFd6BbWrcStH87/A3Io49CFjdpOZZ/idO5p0vLSsGv7Vf8czk8zNfRvSHQd4/xITGAMbYLalGm7iezCbA5lHOJkzkky8jPIyM/gXOG5y1Y1ebt74+/pj5+nH3W96xJWK4zQWqE0qt0IXw/fUt+jOM4I8J7A81rrG4p+fgZAa/1KSc8pb4AXZJ8h951+2PIymaBe4Z5hfUw9SXmt3t7+Nh/t+Yi/9fgbq46vYtOpTYxsMpK/Xfe3cu1RUp3kWHJ4efPLLD28lE71O/Fol0dZFLeI7w5/h6ebJxNaTuDutndTz6d8+9WYIS0vjc/3f8682HnkWnMZGDmQqe2mcjr3NN/EfcPGkxuxazsdQjpwY+Mb6RfZj8iAyl3bn5Gfwbqkdaw6tor1J9Zj1VbaBbdjdNPRDIwayMGMg6xNXMuaxDWk5KYAEBkQSef6nelQvwPNApvRNLBpuU/iWu1WjmcdJz4jnr3pe9mRvIO96Xux2q14u3vTo0EP+kf2p1d4L3Itufye/DtbT29la/LWCyerFYqG/g1pGNCQML8wwmqFUdenLr4evvh5+OHp5olFW7DarRRYC8goyDCWkualk5iVyJFzRy78xuCu3GletzntQ9rTIaQD7UPa467cOXjmIAczDl64vXJKy125X7aqyQ038m355FhyKLAVXHbszEEz6RvRt1x/Xs4I8FuAG7XWU4t+ngj00FrPKOk55Q3w398cT8ezPzMz6k1uv/W2KnGS8lpY7BYmfT+JPWl78HLz4tnrnmVszFi5SvESyxOW88KmF8ix5ODl5sVtLW9jctvJlb523JEyCzKZe2Aucw7MIaswCzBGrKOajmJU01FVZnloel46yxKWsSR+yYXfCAB8PXzpFd6L/pH96Rnes9xz+2WVb81nR8oOfkn6hbWJa/8QlqF+oXQL60an+p1oWa8lMYExFRoAaa1JzUsl9kwsu1J3sSt1F3tS95Brzb3sOIWiUe1GNK/bnBb1WtCibgsia0deWJJa0lXXFruFjPwMknOSOZ17ms71O5f7pLIzAnw8cMMVAd5da/3gFcdNA6YBREVFdTl27Ng1v9fhuFjyj22hzeC7ylVrVZCck8yHez5kfPPxsulTCRKzElmesJzRTUfTwL+B2eU4THZhNssTlhNaK5Q+DftU2XMdWmtiz8Sy4eQGmtdtTo8GPfB29zatlsNnD7P59GZ8PXzpFtaNCP8Ipw96bHYbhzMPsyd1D3bstKjbosL/UDiCS0+hCCFETVZSgFdkx6XfgWZKqWillBfwJ2BpBV5PCCHENSj34kyttVUpNQNYibGM8BOt9T6HVSaEEOKqKrS6Xmu9AljhoFqEEEJcg6rViVcIIUSZSYALIYSLkgAXQggXJQEuhBAuSgJcCCFcVKVuJ6uUSgWu/VJMQzCQ5sByXIF85ppBPnPNUJHP3Ehr/Yf9sis1wCtCKbW1uCuRqjP5zDWDfOaawRmfWaZQhBDCRUmACyGEi3KlAJ9ldgEmkM9cM8hnrhkc/pldZg5cCCHE5VxpBC6EEOISEuBCCOGiXCLAlVI3KqUOKqXilVJPm12PsymlPlFKpSil9ppdS2VQSkUqpdYopQ4opfYppR42uyZnU0r5KKW2KKV2FX3mf5pdU2VRSrkrpXYopZaZXUtlUEodVUrtUUrtVEo5tKNNlZ8DV0q5A4eAIUASRiOJCVrr/aYW5kRKqb5ANvC51rqt2fU4m1KqAdBAa71dKRUAbAPGVPP/xgqopbXOVkp5AuuBh7XWm0wuzemUUo8BXYHaWusRZtfjbEqpo0BXrbXDL1xyhRF4dyBea52gtS4EvgZGm1yTU2mt1wFnzK6jsmitT2mttxd9nwUcABqaW5VzaUN20Y+eRV9VezTlAEqpCOAm4COza6kOXCHAGwKJl/ycRDX/y12TKaUaA52AzeZW4nxFUwk7gRTgJ611tf/MwFvAk4Dd7EIqkQZ+VEptK2ry7jCuEODFtaGu9iOVmkgp5Q98AzyitT5ndj3OprW2aa07AhFAd6VUtZ4uU0qNAFK01tvMrqWS9dJadwaGAdOLpkgdwhUCPAmIvOTnCOCkSbUIJymaB/4GmKu1XmR2PZVJa30WWAvcaHIpztYLGFU0J/w1MFApNcfckpxPa32y6DYFWIwxLewQrhDgvwPNlFLRSikv4E/AUpNrEg5UdELvY+CA1voNs+upDEqpEKVUYNH3vsBgINbcqpxLa/2M1jpCa90Y4+/xaq31nSaX5VRKqVpFJ+ZRStUChgIOW11W5QNca20FZgArMU5uzdda7zO3KudSSn0FbARaKKWSlFJTzK7JyXoBEzFGZDuLvoabXZSTNQDWKKV2YwxSftJa14hldTVMKLBeKbUL2AIs11r/4KgXr/LLCIUQQhSvyo/AhRBCFE8CXAghXJQEuBBCuCgJcCGEcFES4EII4aIkwIUQwkVJgAshhIv6f+n58Pk0zFZYAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"axes = plt.plot(sys.times, x[:, :3])\n",
"leg = plt.legend([sm.latex(q, mode='inline') for q in (q1, q2, q3)])"
]
},
{
"cell_type": "code",
"execution_count": 180,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"<string>:2: RuntimeWarning: overflow encountered in exp\n"
]
}
],
"source": [
"viz_frame = VisualizationFrame('ball', N, P, Sphere(0.2))\n",
"scene = Scene(N, O, viz_frame, system=sys)\n",
"scene.create_static_html(overwrite=True, silent=True)\n",
"matrices = viz_frame.evaluate_transformation_matrix(x, list(sys.constants.values()))"
]
},
{
"cell_type": "code",
"execution_count": 181,
"metadata": {},
"outputs": [],
"source": [
"ground_mesh = pjs.Mesh(pjs.PlaneBufferGeometry(40.0, 40.0),\n",
" pjs.MeshStandardMaterial(color='blue', side='DoubleSide'),\n",
" name=\"ground\")\n",
"\n",
"ball_mesh = pjs.Mesh(pjs.SphereBufferGeometry(0.2),\n",
" pjs.MeshStandardMaterial(color='red'),\n",
" name=\"ball\")"
]
},
{
"cell_type": "code",
"execution_count": 182,
"metadata": {},
"outputs": [],
"source": [
"ball_mesh.matrixAutoUpdate = False\n",
"ball_mesh.matrix = matrices[0]\n",
"ground_mesh.position = [20.0, 20.0, -0.1]"
]
},
{
"cell_type": "code",
"execution_count": 183,
"metadata": {},
"outputs": [],
"source": [
"view_width = 800\n",
"view_height = 600\n",
"camera = pjs.PerspectiveCamera(position=[15, 15, 15], aspect=view_width/view_height)\n",
"key_light = pjs.DirectionalLight(position=[0, 100, 100])\n",
"ambient_light = pjs.AmbientLight()\n",
"scene_pjs = pjs.Scene(children=[ground_mesh, ball_mesh, camera, key_light, ambient_light])\n",
"controller = pjs.OrbitControls(controlling=camera)\n",
"renderer = pjs.Renderer(camera=camera, scene=scene_pjs, controls=[controller], width=view_width, height=view_height)"
]
},
{
"cell_type": "code",
"execution_count": 184,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "d8540ac45bc247b7a6abf1551a58f5c9",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Renderer(camera=PerspectiveCamera(aspect=1.3333333333333333, position=(15.0, 15.0, 15.0), projectionMatrix=(1.…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"renderer"
]
},
{
"cell_type": "code",
"execution_count": 185,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "0492fa862e094f24b9d707fb554b086e",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"AnimationAction(clip=AnimationClip(duration=5.0, tracks=(VectorKeyframeTrack(name='scene/ball.matrix', times=a…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"track = pjs.VectorKeyframeTrack(name='scene/ball.matrix',\n",
" times=list(sys.times),\n",
" values=matrices)\n",
"clip = pjs.AnimationClip(tracks=[track], duration=sys.times[-1])\n",
"action = pjs.AnimationAction(pjs.AnimationMixer(scene_pjs), clip, scene_pjs)\n",
"action"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 186,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAALIAAAAVCAYAAADik7Q+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAG7klEQVRoBe2bi3FUNxRAMeMCgHRAOiCmAqADklQQ6MCUkCEdQAeM0wFOBUzcgekAcAfOObIkJL2fdv2y41lWM7Kkq/vX1ZX0bB9dX1/fa8vR0dErYL9Q/2X+fTt/GB88sGsPEJNPkPmaekl9T1xelToclYEM8gMm/6GeA39TIh76Bw/cBQ8Qo8/R44z6jBi9SDq1gWwmfgfCUULYlzZuUh1wQv2KjT/vi20/mh2s5TvXkTX01hDK/dSJrYv7uYHtxRCjr6gvon3ne2HUDo0geC6pj3cock6U1wpvD7m0gZwn9rjjXevjHtu3umkE8EuYPiYR3Nkk90MFMgvi/cpyyMg3fuj96UmW76O9RLvE+6ECGceGq4XXjF06eQ9kmQDu9OY/3sbJZLa30H0paP3K8b/s2HisGYB+drF4vL2O990A2OBHtSCR9yPofTS8XevojJlfnb1Tpq8/Hs+Wp9RPyPoLPK85v1P1pXAf2lsHDPyU5yfTh7S5AD9l8BT4rxm40Ik2iC9P6xNgPrIu1Z12sYDvPdYPCL69KjrmXE99vu3nXdfte4HRvVSB+qrXERlW9pnTIBV4meD0vW9O0iS8bVr4umE+Jlr6yvfDtws+quMUHBqdKm3QndbFlZ/BLfzVFO2mcHi5QDTBN/rnNPGgn/RQ/vMCbqB/S+NtWuj1lwFT+UaYcy28ZwydeoFa8+wZJ5m0+riyjbE6nfXwaXGgc3MA/q5TNpiJ4EjaJyVC2Wdu4BBgBn8O7BL/Nn15Uo2GByWfCNtYXsmPvkEU+NpS1wxis2zaLPqr2uSM02asAksdqJj6fXE27SuLWm1yxmnjbOwz5cuPWtnQoxc02pn8kDJ5G2+VD3r4Jhx1omZb6d/T8WYNne7LNAsr+8ypzCCwSpw1+8j6Rs3ZWN4Udd1KB+jUX57Z+DX1TbzgX24Qdc1ZN9owtUFNCJW9iWdPq1yq8qqAdRzhVULo4Rn1HSSvHlpkZnn09Xs+lSJfs3Slaw/fEgd6faZ+pzloGYT0TzuakYEPAqtkumZfHaguSpUpGZtJB0dnj2zpqBouj+q476HfFAcZIYBaOuCD7CQOZbDYLe3cOMmjzQEU+Spv44waadPmqDbjnB7tHLLT1a3Vqwrslm5prE3UfDU5ZhAKhOdcwH1oeM/yoZJLvLRrVNeDDvzqMQjvrsdBFnjz2zeH7cNHvVpYQTbejfp71PlI1M6/6fuB/8LxGNUKNqjrmL9cWOXngixh+jfAGbvoVxmhr6O8zyN0Wd4WfE8UXfpoCx5Tev3UmgXvMm60pfJTwgfPa5gJN/9mL2dkgOKNPhbi3CBDCo9z+YiAh5kv7GBaM2t1yU80cy00U9lMHbKsOR7lHDTV/ZNxyjSBF+PqrsZ4DRvMsC3fJLfKcODlrEk/3y1LG5b6UefqahJ5ZZ8x3igLgm885GxO33UZPbGn9APfUzBnTvEoAz7AlFXGTWVLyT/iVnF1H2Bv8TOJuysXdyfVRQi/8Ul9hIYsR2tGepYJ6ESaa1qPhqkS6MFxUUOhr0MsVQbt5KfemQ69UrYb/KYKfsGeKRt65IGj3gatV5iymB2rDBcnzXxJP0+NnJnhteSr4FPolZn9FfmmT38X8DGpZHt77ABfG9La2vcNFU6ZTnpIAn3WSzpgfgpsTyv1e8O8Qa7MpU+FX8HJ5Tj3FjoIfo2Qt1ZQzVgWj0D/tC6V3+j8kQa2rcKMr+Chom6C0SM04hh8yvtE6zEkzeDo7OEHnc4zQMviImuTtvxZTMza0CnPhVPXFJyJvfCx76bKf4EunhxZz05Z8g4bhNZrk+vzRSD0yT7X6AvjfMXr5C2vd1Ev+W1KL43BeRb1CnEjTP2aoo7KSwnLb89XDc70EKYKCxWsEKRpvGkLvefG5JePkh947rzqAVDOt31wzW6TXxw25dfyT2P4dNmwlrwkd65dksW8wZ+vAHO82rkl3i1+O74tvfwijxw3jA3myWsQc4M4vQ+wLY9awAZjj4t8jEjHTswZpuHj8TK646DxmMmFsRnVrJN2a54rOpP8Cpyebq8Na8nr0WlJlqfIhx5GIzhLvEdIKtBt6WVmYLrGqXhtaE+zNDfatlcLU3/JcJRoBui9xqPEQDRI5TU4RmJghuOP+aowp1GntA+LQDeA/a+AUePm+FXM+waLNqwsb1arJVnMmzj086hv5pgv8Z6jde629AV/Y+QEfr4VtMXHoQllqohzVU1CkK8W9ik+wqrXdotz2zH8q+/DJT/mUub1vug3XzP67JcK5if5lbzX6u9S3pIs5vXXtt/Wb+W3Jd3W8nfJJ9rrF6Hqy0/1HyJMhsLO8P7qY8t719jj5Abx8PPggR15gJj0lPdB6K1h8D97/wH/L9OrhgpBSwAAAABJRU5ErkJggg==\n",
"text/latex": [
"$\\displaystyle \\left\\{c_{c}, g, k_{c}, m, \\mu, t, v_{s}\\right\\}$"
],
"text/plain": [
"{c_c, g, k_c, m, μ, t, vₛ}"
]
},
"execution_count": 186,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fr.free_symbols"
]
},
{
"cell_type": "code",
"execution_count": 187,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAIAAAAAVCAYAAACOleY7AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAEeElEQVRoBe2ajVEVMRCAeY4FOFiB0IE/FYgdYAvQAZbgaAfagVoCWoFiB2gFCh3g9z2SmOTu3rscpwzwdiYk2Z/sZnfzd7zFxcXFVg2LxeIA3BPKCfT3NX3TvzkeIJaPsfaQckp5TzzPc+sXeQLA/ADiF8pn8K9yxk37ZnuA2O4xg0+U58T2e5xNnQCu/HcwLCLDpr49HiAJ3jGbp8TX3X0J92Ij1LvUPyrcpnt7POD27y6foE6ARNg07oYHNglwN+I8OMtNAgy65m4Q7rdOk4vEDjLxWeF54jPRy+N3LhefqWeB8CJxXO8lp4z9Ng4caD/pPwLvuTYrXJfu/6R3u3AWDtyKBYLPBN/+CZe3ofmUOKPsRDxtZeheytB+QzmieOOU9iDSWmrl5QeWOnNZcI4PqmsneN+9x320sTjkm3Uj47yd83Gom+etnDYCo+cM72h/w+uCQsVfv6VAQ9inGNzHOUNsg3e1Sz+IOGsgJQ1tHZDood8cDOTcZZbjhDGKpASnk4tx6Rt49VvcMdLcWtrITtFtEOpFUdi8zoaJepv9jZ4TyjLRQvwuVwxIvxSlSdQGK0QBXToWnEnxJgwmz7Id+p2Mq+X7+oyRVk8Y/yjn68NFOjRXz1USoFm3+ijJd7RNRldHwkX7hmp4p+id5G90uWi1+SitEjpxex/aAQx0vepcLU50r29i4DWwkOnjG8I5LsXxc+dE5w7ZeaUEiLa06IZXZyYf0F7pl6ijr27RW8sju9bf8LgDfIqyKQFESKB0AgbO7V+GeiUOrvAgM3ikRANW1Yzh1lqsZm2gnA3JQZsrAZp1R5uwweO0SNxIW1cjN0kvcvGI7l0Y6gU68bonNgO/Apq9Q1B/JXwJY/quHIW4zRqkL5TXKO7QI9+I2tVey78A901Z9ByEm7PdueEqug3iK+Y+5YXSrLfB376oCpvqBOh1YpiIgUjJofPpu9o6Tz/431L83rwLnxe2Agwa5YLidrQKioSDX31PKRG/2+LkBr3aFHXYNtlG6YbP4PvPtPRsDfL/bM7r/K3+DH5n7a1RCRAEnlM/Y4JHFtpxoE6AMwU6Yw9+t8QEIWg6WKe4dQ3B8j+SmU55TaztYMOHIcE+fINexZt1Y5OLwvuK30kKaNDdrDdT1OvvjN5tYli6B0B1gOLMzel5G77iPKFvcIozn747hodPehpWY5gY6YKX06a2GW/tHQCef6U3fwFpR+c8nks340zxdye+fTvANoOPAc/idD4TsHP67grWEXSC/Y8RUdXPglyFvlJXx6yDWfWy8j23de5X2vsW2q7k4hihL8yie6K/Ly3I/t7P2jZ9zoxxoLwGt/61kJfCQxzwi/ohxR3gSV+Q4VGPfLMA46lLp2vXDn1fNF/R3TmLwc+ml7EEL7zOR50J6nnPPWcUjfZ3MEobz5OBNjCyKKC8mKWtrIdupsvj1m7CpK9KNe+qPnK9x8IqmTlo16VX269ZtwvDIzp9r9Cm4hdBWiiQqW5hbvGb3wTqkBsMxNLjKf7zrvObwD83Hza8rBnJTgAAAABJRU5ErkJggg==\n",
"text/latex": [
"$\\displaystyle \\left\\{q_{3}, u_{1}, u_{2}, u_{3}\\right\\}$"
],
"text/plain": [
"{q₃, u₁, u₂, u₃}"
]
},
"execution_count": 187,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"me.find_dynamicsymbols(fr)"
]
},
{
"cell_type": "code",
"execution_count": 188,
"metadata": {},
"outputs": [],
"source": [
"eval_force = sm.lambdify((q3, u1, u2, u3, cc, g, kc, m, mu, vs), fr)"
]
},
{
"cell_type": "code",
"execution_count": 196,
"metadata": {},
"outputs": [],
"source": [
"fr_vals = eval_force(*x[:, 2:].T, sys.constants[cc], sys.constants[g], sys.constants[kc], sys.constants[m], sys.constants[mu], sys.constants[vs])"
]
},
{
"cell_type": "code",
"execution_count": 203,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7f7d7b498748>]"
]
},
"execution_count": 203,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAUMUlEQVR4nO3dbYxc1X3H8d9v1zyFxAmIBTleE9PWigpUScrWQUJqq5AGJyExL4pkqgS/oLKEiESUShHkTRVVlnhRRSlNQLKSCKOkoZaSCisNTZEDjZAIZp2QGGMcVuFpseM1IWAT8NPuvy/mzHpsZp/su3Pn3PP9SKO5c+bemf/d3fnN2XPPnXFECABQhoG6CwAA9A6hDwAFIfQBoCCEPgAUhNAHgIIsqbuAuVx00UWxcuXKussAgKzs2LHj1YgYOrW970N/5cqVGh0drbsMAMiK7Re7tTO8AwAFIfQBoCCEPgAUhNAHgIIQ+gBQEEIfAApC6ANAQQj9PvbS79/Sz35zoO4yADRI35+cVbK/+ddHFCG9cNen6y4FQEPQ0+9jfL8NgKoR+gBQEEIfAApC6ANAQQh9ACgIoQ8ABSH0AaAghD4AFITQB4CCEPoAUBBCHwAKQugDQEEIfQAoCKEPAAUh9AGgIIQ+ABSE0AeAghD6AFAQQh8ACkLoA0BBCH0AKAihDwAFIfQBoCCEPgAUhNAHgILMO/RtD9r+pe0fpdsX2n7Y9nPp+oKOde+0PWZ7j+3rOtqvsr0z3Xe3bVe7OwCA2Sykp3+7pN0dt++QtC0iVknalm7L9uWS1km6QtIaSffYHkzb3Ctpg6RV6bLmjKoHACzIvELf9rCkT0v6VkfzWkmb0/JmSTd0tD8QEUci4nlJY5JW214maWlEPB4RIen+jm0AAD0w357+1yV9WdJUR9slEbFPktL1xal9uaSXO9YbT23L0/Kp7e9ge4PtUdujBw4cmGeJAIC5zBn6tq+XNBERO+b5mN3G6WOW9nc2RmyKiJGIGBkaGprn0wIA5rJkHutcI+mztj8l6VxJS21/V9J+28siYl8auplI649LWtGx/bCkval9uEs7AKBH5uzpR8SdETEcESvVOkD704j4nKStktan1dZLejAtb5W0zvY5ti9T64Dt9jQEdMj21WnWzs0d2wAAemA+Pf2Z3CVpi+1bJL0k6UZJiohdtrdIekbScUm3RcRk2uZWSfdJOk/SQ+kCAOiRBYV+RDwq6dG0/HtJ186w3kZJG7u0j0q6cqFFAgCqwRm5AFAQQh8ACkLoA0BBCH0AKAihDwAFIfQBoCCEPgAUhNAHgIIQ+gBQEEIfAApC6ANAQQh9ACgIoQ8ABSH0AaAghH4GWt8jDwBnjtAHgIIQ+gBQEEI/A4zuAKgKoQ8ABSH0M0BHH0BVCH0AKAihnwGmbAKoCqEPAAUh9AGgIIR+BhjcAVAVQh8ACkLoZ4DjuACqQugDQEEI/QwEo/oAKkLoA0BBCH0AKAihnwEO5AKoCqEPAAWZM/Rtn2t7u+1f2d5l+6up/ULbD9t+Ll1f0LHNnbbHbO+xfV1H+1W2d6b77rbtxdktAEA38+npH5H0sYj4kKQPS1pj+2pJd0jaFhGrJG1Lt2X7cknrJF0haY2ke2wPpse6V9IGSavSZU2F+wIAmMOcoR8tb6abZ6VLSForaXNq3yzphrS8VtIDEXEkIp6XNCZpte1lkpZGxOPR+tjI+zu2wSwY0wdQlXmN6dsetP2UpAlJD0fEE5IuiYh9kpSuL06rL5f0csfm46lteVo+tb3b822wPWp79MCBAwvZHwDALOYV+hExGREfljSsVq/9yllW7zZOH7O0d3u+TRExEhEjQ0ND8ykRADAPC5q9ExGvS3pUrbH4/WnIRul6Iq02LmlFx2bDkvam9uEu7ZgDZ+QCqMp8Zu8M2X5fWj5P0sclPStpq6T1abX1kh5My1slrbN9ju3L1Dpguz0NAR2yfXWatXNzxzYAgB5YMo91lknanGbgDEjaEhE/sv24pC22b5H0kqQbJSkidtneIukZSccl3RYRk+mxbpV0n6TzJD2ULpgDB3IBVGXO0I+IX0v6SJf230u6doZtNkra2KV9VNJsxwMAAIuIM3IBoCCEfgYY3QFQFUIfAApC6GcgOJILoCKEPgAUhNDPAP18AFUh9AGgIIQ+ABSE0M8Ax3EBVIXQB4CCEPo5oKcPoCKEPgAUhNDPAJ+nD6AqhD4AFITQB4CCEPoZYMomgKoQ+gBQEEI/A3T0AVSF0AeAghD6AFAQQj8DfIkKgKoQ+gBQEEI/A/TzAVSF0AeAghD6GWBIH0BVCH0AKAihDwAFIfQzwEcrA6gKoQ8ABSH0c0BHH0BFCH0AKAihnwE6+gCqQugDQEHmDH3bK2w/Ynu37V22b0/tF9p+2PZz6fqCjm3utD1me4/t6zrar7K9M913t20vzm4BALqZT0//uKR/iog/l3S1pNtsXy7pDknbImKVpG3pttJ96yRdIWmNpHtsD6bHulfSBkmr0mVNhfvSWJyRO7Od42/o6VfeqLsMIBtzhn5E7IuIX6TlQ5J2S1ouaa2kzWm1zZJuSMtrJT0QEUci4nlJY5JW214maWlEPB6tzwq+v2Mb4LR85huP6fp/f6zuMoBsLGhM3/ZKSR+R9ISkSyJin9R6Y5B0cVptuaSXOzYbT23L0/Kp7d2eZ4PtUdujBw4cWEiJjcTJWQCqMu/Qt/1uST+Q9MWIODjbql3aYpb2dzZGbIqIkYgYGRoamm+JAIA5zCv0bZ+lVuB/LyJ+mJr3pyEbpeuJ1D4uaUXH5sOS9qb24S7tmANj+gCqMp/ZO5b0bUm7I+JrHXdtlbQ+La+X9GBH+zrb59i+TK0DttvTENAh21enx7y5YxsAQA8smcc610j6vKSdtp9KbV+RdJekLbZvkfSSpBslKSJ22d4i6Rm1Zv7cFhGTabtbJd0n6TxJD6ULAKBH5gz9iHhM3cfjJenaGbbZKGljl/ZRSVcupEBwRi6A6nBGLgAUhNDPQHAkF0BFCH0AKAihDwAFIfQzwOgOgKoQ+gBQEEIfAApC6ANAQQj9DDCmD6AqhD4AFITQB4CCEPoZ4EtUAFSF0AeAghD6GeBALoCqEPoAUBBCPwN09AFUhdAHgIIQ+gBQEEI/A3yJCoCqEPoAUBBCPwP08wFUhdAHgIIQ+gBQEEI/AxzHBVAVQh8ACkLoZ4GuPoBqEPpAH3r6lTf0pf98SlNTvOGjWoR+BhjTL8+G+0f1w1++on0HD9ddChqG0Af60DlnDUqSDh+brLkSNA2hD/Qhp2v+y0PVCP0M8Lovj9upz28fFSP0gT7klPocx0XVCP0M8C9+eRjewWKZM/Rtf8f2hO2nO9outP2w7efS9QUd991pe8z2HtvXdbRfZXtnuu9u+8Q/sABO1n51BMM7qNh8evr3SVpzStsdkrZFxCpJ29Jt2b5c0jpJV6Rt7rE9mLa5V9IGSavS5dTHxAx44ZfHqa9PTx9VmzP0I+Jnkl47pXmtpM1pebOkGzraH4iIIxHxvKQxSattL5O0NCIej9Y3gtzfsQ2AU0z39Al9VOx0x/QviYh9kpSuL07tyyW93LHeeGpbnpZPbe/K9gbbo7ZHDxw4cJolAvmbIvVRsaoP5HYbp49Z2ruKiE0RMRIRI0NDQ5UVlyte9+U5MXuHXz6qdbqhvz8N2ShdT6T2cUkrOtYblrQ3tQ93aQfQRbuXxJRNVO10Q3+rpPVpeb2kBzva19k+x/Zlah2w3Z6GgA7ZvjrN2rm5YxvMgc5eedpj+pOkPiq2ZK4VbH9f0t9Kusj2uKR/lnSXpC22b5H0kqQbJSkidtneIukZSccl3RYR7Q8PuVWtmUDnSXooXQB0ceJALqGPas0Z+hFx0wx3XTvD+hslbezSPirpygVVB0lM2SxRe8omPX1UjTNygT7U7unnkPljE4f0H0+8VHcZmKc5e/oAeu/Egdz+T/1P/dtjOjo5pX/46KV1l4J5oKefgQxe97Vr3Nh3RlM2j05OSWrg76ChCH00QtPyZiDD2TvHM6q1ZIQ+GqFpcZPjp2zm9AZVMkIfjZDDMMhCtM/IzSlIj6VhHvQ3Qh+N0LDMz+pAbltOb1AlI/QzkNHrvjZNO5fhxJTNfPbr2GQ+tZaM0EcjZJSN89I+OSunzjM9/TwQ+hloWi92MTQt9JXh7J3JPvkl5PQzqwOhj0Zo2htjlmP6fTC88/Qrb+hPv/Jj/d9v+B6OmRD6GcjodV+bpnXuchzT74ee/pMvtL7k76e799dcSf8i9NEITTsbdHpMP6NZkJN9VGx7yiveidBHIzQr8js+Tz+jN7N+mKaf0Y+rNoR+Bvg7nlv0QeBUKcfP0z/eRz19zIzQRyM070Bu+4zcmgtZgH7I/Lr+Ct4+OqnvPfGipjI4uMRHK2cgp95eXZr2I8rxQG7JPf1vPjKmbzwypqXnnqXPfOj9dZczK3r6yFbnm2E+0bgwOYV+TrVW7Q9vHZUkvf72sZormRuhn4FyX0qz68yYpgVOe/ZJDsMFbcf7YJ4+5kboI1udEdOwzJ8+OSunHOVM2DwQ+sjWycM7zQqc9peo5NDT76fppe2/iV5P03dGX4BA6Gcgg7+jWjS6p5/R1yW2866fvjnL4uSsmRD6yFZnHmaQjQtyYnin/3esn44/tN94cnizrAuhnwX+gLvpHNJp2vDOiZOz6q1jPvqpp89xhbkR+hnI6QSdXjp59k59dSyOfL4u0X10/KE9g6iuc1ty+H0R+hnI4Q+pbk07gS2nk7Pa4+f90dNv9ZDqGhY7mkEPjdDvY87wizR6qclj+u396Yfe85z66A2q/cZz7HhNoX+c0McZGHS7B9X/f0h1OGlMv/68qVT7P5eMMr8vTs5qd5CO1fSaIfRxRgYG8pm2V4eTevoNO5DbHp7IY/ZO67ofam339Hv9BtR+vsOEPs7EdE+/D3pQ/ajJ8/QnM5p66D466Nyuodf/HbfH8g8fm+zp854OQr+PMaY/u86DtzmE40JMh34Gv/t++jtth/3RHo/pH0k9/CPH6OnjDAykV1M//Nvc75r2EzpxklHNhSxAP4R+bT39FPqHj9PTxxnopx5UP2ry8E67h5/D7366c9IHtbaHQns9JDod+gzvvJPtNbb32B6zfUevnz8n/TQroh+dHPTN+hm1e/o5nX/QD6HfrqHX8+WPpB7+YYZ3TmZ7UNI3JX1S0uWSbrJ9eS9ryEl79g7DO911jnf3w4lBVZrKaPZOO2iP9MHQRnt45UiPZ9Hk1NPv9dclrpY0FhG/lSTbD0haK+mZHtfRN6amQpMRmorQ1FTrRT451bq0Q23f64c1cfBwzZX2n4lDR6aXX3j1j7rgXWfXWE21Jg629u13b/T/7/7tFHR79r+p/fOsdSHvZQuZjvv8q29JksZfe0v7Dx6efp72Y5y43VnLyfedWuN8tj14+Lgkae8bb1f6+xp6zznTH2hXFffy30fbfy9pTUT8Y7r9eUkfjYgvnLLeBkkbJOnSSy+96sUXX1zwc9206ed6/tU/djxmx+Of/Fyz1DvDcscjdLZHtHo9Uym4pyI01W5L4T45Fa310jIAzOTZf1mjc88aPK1tbe+IiJFT23vd0++WsO9IvojYJGmTJI2MjJxWMv7Vygt06YXvSk/Q/czN2Q4EntS7mHGbOKl9wJbdml8/OGANDFiDtgas6eV2+0Bab3qdAWvA1uBA63GWDFhLBgd09pIBHcvg8zzq8q6zB2VZbx3t/3+rF2r4gvP08h/eqruMOZ09OKA/GTpfz/7u0II+x34hHdiF9HX/Yvi92vXKwekhv/bztB/jxO2T7zhxv7uvf8p2J3cKrcuXLdWOF1+rdKhxyUC1vXyp96E/LmlFx+1hSXsX44m+9IkPLsbDApjBVR+4sO4Spl3x/vfW8rx/dvG7a3nehej17J0nJa2yfZntsyWtk7S1xzUAQLF62tOPiOO2vyDpJ5IGJX0nInb1sgYAKFmvh3cUET+W9ONePy8AgDNyAaAohD4AFITQB4CCEPoAUBBCHwAK0tOPYTgdtg9IWvjnMLRcJOnVCsvJAftchtL2ubT9lc58nz8QEUOnNvZ96J8J26PdPnuiydjnMpS2z6Xtr7R4+8zwDgAUhNAHgII0PfQ31V1ADdjnMpS2z6Xtr7RI+9zoMX0AwMma3tMHAHQg9AGgII0MfdtrbO+xPWb7jrrr6QXb37E9YfvpumvpBdsrbD9ie7ftXbZvr7umxWb7XNvbbf8q7fNX666pV2wP2v6l7R/VXUsv2H7B9k7bT9kerfSxmzamb3tQ0m8k/Z1a39T1pKSbIqLRX75u+68lvSnp/oi4su56FpvtZZKWRcQvbL9H0g5JNzT59+zW9/idHxFv2j5L0mOSbo+In9dc2qKz/SVJI5KWRsT1ddez2Gy/IGkkIio/Ia2JPf3VksYi4rcRcVTSA5LW1lzToouIn0l6re46eiUi9kXEL9LyIUm7JS2vt6rFFS1vpptnpUuzem1d2B6W9GlJ36q7liZoYugvl/Ryx+1xNTwMSmd7paSPSHqi3koWXxrmeErShKSHI6Lx+yzp65K+LGmq7kJ6KCT9r+0dtjdU+cBNDP1uXx/f+N5QqWy/W9IPJH0xIg7WXc9ii4jJiPiwpGFJq203eijP9vWSJiJiR9219Ng1EfGXkj4p6bY0fFuJJob+uKQVHbeHJe2tqRYsojSu/QNJ34uIH9ZdTy9FxOuSHpW0puZSFts1kj6bxrgfkPQx29+tt6TFFxF70/WEpP9Sa9i6Ek0M/SclrbJ9me2zJa2TtLXmmlCxdFDz25J2R8TX6q6nF2wP2X5fWj5P0sclPVtvVYsrIu6MiOGIWKnWa/mnEfG5mstaVLbPT5MTZPt8SZ+QVNmsvMaFfkQcl/QFST9R6+DelojYVW9Vi8/29yU9LumDtsdt31J3TYvsGkmfV6vn91S6fKruohbZMkmP2P61Wp2bhyOiiCmMhblE0mO2fyVpu6T/joj/qerBGzdlEwAws8b19AEAMyP0AaAghD4AFITQB4CCEPoAUBBCHwAKQugDQEH+H6fb+Z2cBixJAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(sys.times, fr_vals.squeeze().T[:, 2])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.7"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
@renefritze
Copy link

Since I'm struggling with an Animation setup in pythreejs myself, I wanted to replicate your excellent demo. Put the notebook in a repo, added a conda description and now you can try it out on mybinder.org

@moorepants
Copy link
Author

Note that the development version of pydy works like this: https://pydy.readthedocs.io/en/latest/examples/multidof-holonomic.html

I'll try to do a new release soon.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment