stephentu/cartpole_noncollocated.ipynb

## cartpole_noncollocated.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import jax.numpy as jnp\n",
    "import matplotlib.pylab as plt\n",
    "import jax\n",
    "import scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = 9.81\n",
    "cartpole_params = (5.0, 1.0, 1.0, 0.0, 0.0)\n",
    "def cartpole_dynamics(q, q_d, params):\n",
    "    x, theta = q\n",
    "    x_d, theta_d = q_d\n",
    "    m_c, m_p, r, k1, k0 = params\n",
    "    \n",
    "    M = jnp.array([\n",
    "        [m_p + m_c, m_p * r * jnp.cos(theta)],\n",
    "        [m_p * r * jnp.cos(theta), m_p * (r ** 2)]\n",
    "    ])\n",
    "    \n",
    "    v = jnp.array([\n",
    "        -m_p * r * jnp.sin(theta) * (theta_d ** 2) + k1 * x_d + k0 * x,\n",
    "        m_p * g * r * jnp.sin(theta)\n",
    "    ])\n",
    "    \n",
    "    return M, v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def angdiff(th1, th2):\n",
    "    d = th1 - th2\n",
    "    d = jnp.mod(d+jnp.pi, 2*jnp.pi) - jnp.pi\n",
    "    return d\n",
    "\n",
    "def compute_feedback(q, q_d, thetades, thetades_d, thetades_dd, cartpole_params, pd_params):\n",
    "    k0, k1 = pd_params\n",
    "    M, v = cartpole_dynamics(q, q_d, cartpole_params)\n",
    "    theta_dd = thetades_dd -  k1 * (q_d[1] - thetades_d) - k0 * angdiff(q[1], thetades)\n",
    "    f_x = -M[0, 0] * v[1] / M[0, 1] + theta_dd * ( - M[0, 0] * M[1, 1] / M[0, 1] + M[0, 1] ) + v[0]\n",
    "    \n",
    "    return f_x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_pd_params(pd_params):\n",
    "    k0, k1 = pd_params\n",
    "    roots = np.roots([1, k1, k0])\n",
    "    return np.all(np.real(roots) < 0) \n",
    "pd_params = (1, 1)\n",
    "assert check_pd_params(pd_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def thetades(t):\n",
    "    return 0.0\n",
    "thetades_d = lambda _: 0.0\n",
    "thetades_dd = lambda _: 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jax.jit\n",
    "def closed_loop(t, state):\n",
    "    q, q_d = state[:2], state[2:]\n",
    "    f_x = compute_feedback(q, q_d, thetades(t), thetades_d(t), thetades_dd(t), cartpole_params, pd_params)\n",
    "    M, v = cartpole_dynamics(q, q_d, cartpole_params)\n",
    "    q_dd = jax.scipy.linalg.solve(M, jnp.array([f_x, 0.0]) - v, sym_pos=True)\n",
    "    return jnp.hstack((q_d, q_dd))\n",
    "\n",
    "@jax.jit\n",
    "def compute_inputs(t, state):\n",
    "    q, q_d = state[:2], state[2:]\n",
    "    return compute_feedback(q, q_d, thetades(t), thetades_d(t), thetades_dd(t), cartpole_params, pd_params)\n",
    "\n",
    "@jax.jit\n",
    "def compute_accelerations(t, state):\n",
    "    return closed_loop(t, state)[2:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/stephentu/anaconda3/lib/python3.7/site-packages/jax/lib/xla_bridge.py:122: UserWarning: No GPU/TPU found, falling back to CPU.\n",
      "  warnings.warn('No GPU/TPU found, falling back to CPU.')\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "DeviceArray([ 0.       ,  0.       , 13.427354 ,  1.0000004], dtype=float32)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "state_0 = jnp.array([0.0, -1, 0.0, 0.0])\n",
    "closed_loop(0.0, state_0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_end = 10\n",
    "t_eval = np.linspace(0, t_end, num=100)\n",
    "res = scipy.integrate.solve_ivp(closed_loop, (0, t_end), state_0, t_eval=t_eval)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = res.y.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 4)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "traj.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a25f27350>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xUVf7/8dchlNBDJ5QQqjF0iIIdRXcRUVDsqNgWdXXtCNiwrC64rmW/YsGy6uLSAgqoawFBl1VAWgolEEILCQQkDULqnN8fM/jLsmGBTCY3c+f9fDzyyNw7d3I/N3fy5nDm3nOMtRYREXGXWk4XICIiVU/hLiLiQgp3EREXUriLiLiQwl1ExIVqO10AQMuWLW10dLTTZYiIBJU1a9YcsNa2qui5GhHu0dHRrF692ukyRESCijFm5/GeU7eMiIgLKdxFRFxI4S4i4kIKdxERF1K4i4i4kMJdRMSFFO4iIi6kcBcRcUBpmYdpS1NJ2J0TkJ9fI25iEhEJJVv25fPo3AQS03M5VFRK344RVb4PhbuISDUpLfPwzg9pvL54K43CazPtxgFc1icyIPtSuIuIVIOUvfmMj/e21i/rHclzI3vSolG9gO1P4S4iEkDlW+uNw2vz5pgBDO8dmNZ6eQp3EZEASdnr7VtP2lM9rfXyFO4iIlWstMzD299v4/UlW2kSXqfaWuvlKdxFRKrQ5r15jJ+b6G2t94nkuSuqr7VensJdRKQKlJR5eMfh1np5CncRET9t3pvHo3MTSN6Tx4g+kTzrUGu9PIW7iEgllZR5eHvZNv76nbe1/taYAVzqYGu9PIW7iEglbMrMY3y8t7V+ed92PHtFT5o3rOt0Wb864dgyxpgPjDFZxpjkcuv+bIzZbIxJNMZ8aoyJKPfcJGNMqjEmxRjz20AVLiLihJIyD39dspUr3ljO3txC3hozgP+7oX+NCnY4uYHDPgSGHbPuW6CXtbYPsAWYBGCMiQWuB3r6XvOmMSasyqoVEXHQxow8Rk37N698u4VhvSL55qELakw3zLFO2C1jrf3BGBN9zLpvyi2uAK72PR4JzLLWFgHbjTGpwJnAT1VSrYiIA0rKPLy5dBv/991WIhrU4e2bBjCsV80M9aOqos/9dmC273F7vGF/VLpvnYhIUNqY4b0SZmNmHlf0bcczNaxv/Xj8CndjzBNAKfDJ0VUVbGaP89pxwDiAqKgof8oQEalyJb7x1t/4LtXXWh/IsF5tnS7rpFU63I0xY4ERwFBr7dEATwc6ltusA5BR0euttdOB6QBxcXEV/gMgIuKE8q31kf3a8czlPWkWBK318ioV7saYYcAE4AJrbUG5pxYC/zDGvAK0A7oDq/yuUkSkGhSXenhz2dHWet2ga62Xd8JwN8bMBIYALY0x6cBkvFfH1AO+NcYArLDW3m2t3WCMmQNsxNtdc6+1tixQxYuIVJUNGbmMn5vIxsw8RvVrx+QgbK2XZ/5/j4pz4uLi7OrVq50uQ0RCUHGpt2992tJUmjWsywujevGbnsHRWjfGrLHWxlX0nO5QFZGQlbwnl0fnJrB5bz5X9m/P5MtjiWgQvK318hTuIhJyiks9vPHdVqYt20bzhnWZfvPAoGmtnyyFu4iElKT0XMbHe1vrVw1oz9Mj3NNaL0/hLiIhoai0jL8u2crb36fRomFd3rsljotj2zhdVsAo3EXE9RJ25zA+PoEt+w5x9cAOPHVZLE0b1HG6rIBSuIuIaxWWlPH6kq288/02WjcO54Nb47goxr2t9fIU7iLiSut2ZTM+PpHUrENcG9eBJy6LpWl9d7fWy1O4i4irFJaU8eriLbz7QxptmoTz4W1nMOS01k6XVe0U7iLiGmt2ZvNYfALb9h/mhjM7Mmn46TQJD53WenkKdxEJeoUlZfzlmxTeW76ddk3r8/HtZ3J+j1ZOl+UohbuIBLXVOw4yPj6R7QcOM2ZQFJOGn06jeoo2/QZEJCgVFJfy569T+PDHHbSPqM8/7hzE2d1aOl1WjaFwF5GgsyLtFybMS2TnLwWMPasTjw2LoaFa6/9Bvw0RCRqHi0qZ+tVmPv5pJ51aNGDWuMEM7tLC6bJqJIW7iASFf6ceYMK8RPbkHOH2czrz6G970KCuIux49JsRkRotv7CEF7/czMxVu+jSsiFz7zqLuOjmTpdV4yncRaTGWpaSxePzk9ibV8i487vw8CU9CK8T5nRZQUHhLiI1Tm5BCc9/sZH4Nel0a92IefecTf+oZk6XFVQU7iJSoyzeuI/HP03il8PF3HthV+4f2p16tdVaP1UKdxGpEbIPF/PMog0sWJ9BTNvGvD/2DHp3aOp0WUFL4S4ijvtnUiZPLUgm90gJD13cg3uGdKVu7VpOlxXUFO4i4pj9+UVMXpjMl0l76d2+KTPuHERM2yZOl+UKJwx3Y8wHwAggy1rby7euOTAbiAZ2ANdaa7ONMQZ4HRgOFAC3WmvXBqZ0EQlW1loWrM/gmUUbKCgu47FhpzHuvC7UDlNrvaqczG/yQ2DYMesmAkustd2BJb5lgEuB7r6vccBbVVOmiLjF3txC7vxoNQ/OXk/nlg358v5z+f2Qbgr2KnbClru19gdjTPQxq0cCQ3yPPwKWARN86z+21lpghTEmwhgTaa3NrKqCRSQ4WWuZ/fNuXvhiEyUeD0+NiOXWs6MJq2WcLs2VKtvn3uZoYFtrM40xR6c5aQ/sLrddum/df4W7MWYc3tY9UVFRlSxDRILB7oMFTJqfxPLUAwzu0pypo/vQqUVDp8tytar+QLWif4JtRRtaa6cD0wHi4uIq3EZEgpvHY/n4px1M/SqFsFqGF67sxQ1nRFFLrfWAq2y47zva3WKMiQSyfOvTgY7ltusAZPhToIgEp237DzEhPpHVO7M5v0cr/nRVb9pH1He6rJBR2XBfCIwFpvi+Lyi3/j5jzCxgEJCr/naR0FJa5uHdf23n1cVbqF8njJev6cvoAe3xXkwn1eVkLoWciffD05bGmHRgMt5Qn2OMuQPYBVzj2/xLvJdBpuK9FPK2ANQsIjXUxow8JsxLJGlPLr/t2YbnR/WideNwp8sKSSdztcwNx3lqaAXbWuBef4sSkeBSVFrGtO9SeXPZNiIa1OHNMQMY3jvS6bJCmu5QFRG/rNuVzWPxiWzNOsSV/dvz9IhYmjWs63RZIU/hLiKVcqS4jL98k8L7/95O2ybh/O3WM7gwpvWJXyjVQuEuIqfsx20HmDgviV0HC7hxUBSTLo2hcXgdp8uSchTuInLS8gpL+NOXm5i5ajedWjRg5u8Gc1ZXTVBdEyncReSkLN64jyc+S2J/fhHjzu/CQxf3oH5dTaJRUyncReR/OnCoiGcXbWRRgncSjek3x9G3Y4TTZckJKNxFpEJHh+V9dtEGDheV8fAlPbj7Ak2iESwU7iLyX/bkHOHJT5NYmrKfAVERTB3dh+5tGjtdlpwChbuI/MrjscxYuZOp/9yMx8LTI2IZq2F5g5LCXUQASM06xKT5ify8I5vzurfkxSt707F5A6fLkkpSuIuEuJIyD9N/SOP1xVupX1cDfbmFwl0khCWm5zBhXhKbMvO4rHckk6+I1UBfLqFwFwlBR4rLeHXxFt77VxotG9XjnZsH8tuebZ0uS6qQwl0kxPyYeoCJ871DB9xwZkcmXno6Tetr6AC3UbiLhIjcghJe/HITs1fvJlpDB7iewl0kBHyVnMlTCzZw8HAxd1/QlQcv7k54HQ0d4GYKdxEX25dXyNMLkvl6wz5iI5vwt1vPoFf7pk6XJdVA4S7iQtZaZv28mxe/3ERxqYcJw2K487zO1AnT0AGhQuEu4jLbDxxm0vxEVqQdZFDn5kwZ3YfOLRs6XZZUM4W7iEuUlHl4919pvLZ4K/Vq12LKVb257oyOuhkpRCncRVwgKT2XCfMS2ZiZx7CebXluZE9aN9HNSKFM4S4SxI69GentmwYyrJduRhI/w90Y8xBwJ2CBJOA2IBKYBTQH1gI3W2uL/axTRI6xfOsBHv/06M1IUUy8NEY3I8mvKv3RuTGmPXA/EGet7QWEAdcDU4FXrbXdgWzgjqooVES8cgqKGT83gZveX0lYLcOscYP501W9FezyH/ztlqkN1DfGlAANgEzgIuBG3/MfAc8Ab/m5H5GQZ63li6RMnlm4gZyCEn4/pCv3D9XNSFKxSoe7tXaPMeZlYBdwBPgGWAPkWGtLfZulA+0rer0xZhwwDiAqKqqyZYiEhMzcIzz1WTKLN2XRp0NTPr59ELHtmjhdltRglQ53Y0wzYCTQGcgB5gKXVrCprej11trpwHSAuLi4CrcRCXUej+WTlTuZ+lUKZR7Lk5edzm3ndNbMSHJC/nTLXAxst9buBzDGzAfOBiKMMbV9rfcOQIb/ZYqEntSsfCbOS2L1Ts2MJKfOn3DfBQw2xjTA2y0zFFgNLAWuxnvFzFhggb9FioSS4lIPby3bxrSlqTSop5mRpHL86XNfaYyJx3u5YymwDm83yxfALGPMH33r3q+KQkVCwZqd2Uyan8iWfYe4om87nr48lpaN6jldlgQhv66WsdZOBiYfszoNONOfnysSag4VlfLy1yl89NMOIpuE88GtcVwU08bpsiSI6Q5VEYctTcniyU+Tycg9wi2DOzF+WAyN6ulPU/yjd5CIQ345VMTzn2/ks/UZdGvdiPi7z2Jgp+ZOlyUuoXAXqWbWWhasz+C5zzeSX1jCA0O78/sLu1Kvtm5GkqqjcBepRunZBTz5WTLLUvbTr2MEU0f34bS2jZ0uS1xI4S5SDco8lr//tIOXvk4BYPLlsdxyVrRuRpKAUbiLBNjWfflMmJfI2l05nN+jFS9e2YsOzXQzkgSWwl0kQIpLPbz9/Tbe+C6VhvXCePW6vozqp5uRpHoo3EUCYP3uHCbEJ5KyL5/L+7Zjsm5GkmqmcBepQkeKy3jl2xTeX76d1o3Dee+WOC6O1c1IUv0U7iJV5MfUA0yc750Z6cZB3pmRmoRrAg1xhsJdxE95hSX86ctNzFy1m+gWDZg1bjCDu7RwuiwJcQp3ET8s3riPJz5LYn9+EXed34WHLumhmZGkRlC4i1TCL4eKeHbRRhYmZBDTtjHv3hJHnw4RTpcl8iuFu8gpsNbyeWImkxduIL+whIcu7sE9Q7pSt3al55oXCQiFu8hJysor5MnPkvlm4z76dmjKS1cP1tABUmMp3EVOwFrLvLV7eG7RBopKPUy6NIY7zu1M7TC11qXmUriL/A8ZOUd4/NMklqXs54zoZkwd3YcurRo5XZbICSncRSpgrWXWz7t54YtNlHksky+PZexZ0dTSQF8SJBTuIsdIzy5g4rwklqce4KwuLZg6ug9RLTTQlwQXhbuIj7WWT1bu4k9fbgLgj6N6ceOZUWqtS1BSuIvgba1PmJfIv1N/4ZxuLZhyVR86NldrXYKXX+FujIkA3gN6ARa4HUgBZgPRwA7gWmtttl9VigSItZZ/rNrFi194W+svXOltrWtYXgl2/l7L9TrwlbU2BugLbAImAkustd2BJb5lkRpnT84RbvlgFU98mky/qAi+evB8xgzqpGAXV6h0y90Y0wQ4H7gVwFpbDBQbY0YCQ3ybfQQsAyb4U6RIVbLWMmf1bp7/fBMea3l+VC9uGqTWuriLP90yXYD9wN+MMX2BNcADQBtrbSaAtTbTGNO6ohcbY8YB4wCioqL8KEPk5O3NLWTi/ESWpexncJfm/PnqvupbF1fyJ9xrAwOAP1hrVxpjXucUumCstdOB6QBxcXHWjzpETshay4L1GTy9IJniMg/P+Cao1pUw4lb+hHs6kG6tXelbjscb7vuMMZG+VnskkOVvkSL++OVQEU98msxXG/YysFMzXr6mL51bNnS6LJGAqnS4W2v3GmN2G2NOs9amAEOBjb6vscAU3/cFVVKpSCV8s2Evj3+aRN6RUiZeGsPvzutCmFrrEgL8vc79D8Anxpi6QBpwG94rcOYYY+4AdgHX+LkPkVOWV1jCsws3Mm9tOrGRTfjkzn4awVFCil/hbq1dD8RV8NRQf36uiD9+3HaA8XMTycw9wn0XduP+od013rqEHN2hKq5RWFLGn79O4f3l2+ncsiHx95zNgKhmTpcl4giFu7jChoxcHpy1nq1Zh7h5cCcmDY+hQV29vSV06d0vQa3MY5n+QxqvfJtCswZ1+ej2M7mgRyunyxJxnMJdglZ6dgEPz0lg1faDDO/dlhdG9aZZw7pOlyVSIyjcJegcvSHpqc+SscBfrunLVQPaa/gAkXIU7hJUco+U8ORnySxKyOCM6Ga8cm0/DR8gUgGFuwSNlWm/8PCcBPblFfLob3pwz5BuuiFJ5DgU7lLjlZR5eG3xFt5cto3oFg2Zd8/Z9O0Y4XRZIjWawl1qtB0HDvPA7PUk7M7h2rgOTL68Jw3r6W0rciL6K5EayVrL/LV7eHpBMmG1DG+OGcDw3pFOlyUSNBTuUuPkF3o/NF2wPoMzOzfntev60S6ivtNliQQVhbvUKOt353D/zHXea9gv6cG9F+pDU5HKULhLjeDxWKb/K42Xv06hTZNw5tx1FnHRzZ0uSyRoKdzFcfvzi3h4znr+tfUAl/Zqy5Sr+tC0QR2nyxIJagp3cdTyrQd4cPZ68gtLeOHKXtx4piaqFqkKCndxRGmZh9cWb2XaslS6tWrEJ3cO0mQaIlVI4S7VLjP3CPfPXMfPO7K5Lq4jk6+I1fC8IlVMf1FSrZZuzuLhOespLvXw2nX9GNW/vdMlibiSwl2qRUmZh5e/SeGd79M4PbIJ027sT5dWjZwuS8S1FO4ScJm5R7jvH+tYszObGwdF8fSIWMLrhDldloirKdwloJalZPHQbG83zOvX92NkP3XDiFQHhbsExNGrYd5YmkpM28ZMGzOAruqGEak2foe7MSYMWA3ssdaOMMZ0BmYBzYG1wM3W2mJ/9yPBIyu/kPtnrmNF2kGuP6Mjz1zRU90wItWsVhX8jAeATeWWpwKvWmu7A9nAHVWwDwkSK9J+4bK/Lmf97hxevqYvU0b3UbCLOMCvcDfGdAAuA97zLRvgIiDet8lHwCh/9iHBweOxvLVsGze+u4LG9Wqz4N5zuXpgB6fLEglZ/nbLvAY8Bhy9tbAFkGOtLfUtpwMVfoJmjBkHjAOIioryswxxUu6REh6Zk8DiTfu4rHckU0b3pnG4xoYRcVKlw90YMwLIstauMcYMObq6gk1tRa+31k4HpgPExcVVuI3UfBsycrlnxloyco7w9IhYbjsnWmPDiNQA/rTczwGuMMYMB8KBJnhb8hHGmNq+1nsHIMP/MqUmil+TzhOfJhHRoA6z7xrMwE4aolekpqh0n7u1dpK1toO1Nhq4HvjOWjsGWApc7dtsLLDA7yqlRikqLWPS/CQenZvAgKhmfHH/eQp2kRomENe5TwBmGWP+CKwD3g/APsQhe3KO8PsZa0hIz+WeIV155JIe1A6riouuRKQqVUm4W2uXAct8j9OAM6vi50rN8q+t+7l/5jpKyyzv3DyQ3/Zs63RJInIcukNVTshay5vLtvGXb1Lo1roRb980UIN+idRwCnf5n/ILvZc5frNxHyP6RDJ1dB8a1tPbRqSm01+pHNfWffnc9fc17DxYwFMjYrldlzmKBA2Fu1Ton0mZPDo3gfp1w/jkzkEM7tLC6ZJE5BQo3OU/lHksf/46hbe/30b/qAjeGjOQtk3DnS5LRE6Rwl1+lX24mD/MXMfy1AOMGRTF05fHUq+2Bv0SCUYKdwEgeU8ud89YQ1ZeES+N7sO1Z3R0uiQR8YPCXZi/Np1J85No3rAuc+4+i34dI5wuSUT8pHAPYSVlHl74YhMf/riDQZ2bM23MAFo2qud0WSJSBRTuIWp/fhH3/mMtq7Yf5PZzOjNpeAx1NIyAiGso3ENQwu4c7p6xhuyCYl69ri9X9tekGiJuo3APMXNW7+bJz5Jp1age8XefTa/2TZ0uSUQCQOEeIopLPTz3+QZmrNjFud1a8tcb+tO8YV2nyxKRAFG4h4CsvELu+WQta3Zmc/cFXRn/29MIq6VhBETcTOHucqt3HOSeT9ZyqLCUN27sz4g+7ZwuSUSqgcLdpay1zFi5i+cWbaBdRH1m3DGI09o2PvELRcQVFO4uVFhSxtMLkpmzOp0LT2vFa9f3p2n9Ok6XJSLVSOHuMntyjnD339eQtCeX+y/qxoMX96CW+tdFQo7C3UV+TD3AfTPXUVLq4d1b4rgkto3TJYmIQxTuLmCtZfoPaUz9ajNdWjXinZsH0lXT4ImENIV7kDtUVMpj8Ql8mbSX4b3b8tLVfWmkafBEQl6lU8AY0xH4GGgLeIDp1trXjTHNgdlANLADuNZam+1/qXKs1KxD3D1jDWn7D/H48Bh+d14XTYMnIgD4M1JUKfCItfZ0YDBwrzEmFpgILLHWdgeW+Jaliv0zKZORbywn+3AxM+4YxLjzuyrYReRXlW65W2szgUzf43xjzCagPTASGOLb7CNgGTDBryrlV6VlHl76OoXpP6TRr2MEb900gMim9Z0uS0RqmCrpnDXGRAP9gZVAG1/wY63NNMa0Ps5rxgHjAKKioqqiDNfLyivkvpnrWLX9IDcNjuKpEZoGT0Qq5ne4G2MaAfOAB621eSfbNWCtnQ5MB4iLi7P+1uF2K9N+4b6Z68gvLNEwvSJyQn6FuzGmDt5g/8RaO9+3ep8xJtLXao8EsvwtMpR5PJZ3fkjj5W9S6NS8gYYREJGT4s/VMgZ4H9hkrX2l3FMLgbHAFN/3BX5VGMJyCop5ZE4CSzZncVnvSKaM7k3jcA0jICIn5k/L/RzgZiDJGLPet+5xvKE+xxhzB7ALuMa/EkPT+t053PvJWrLyC3nm8ljGnh2tq2FE5KT5c7XMcuB4aTO0sj831FlreX/5dqb8czNtmoQz566z6B/VzOmyRCTI6FbGGiSnoJjx8Yl8u3Efl8S24eWr+9K0gbphROTUKdxriFXbD/LArHUcOFTEUyNiuf0cdcOISOUp3B1W5rG88V0qry/ZQlTzBsy/5xx6d9Ck1SLiH4W7g9KzC3ho9np+3pHNqH7t+OOVvTXol4hUCSWJQxYlZPD4p0lYC69e15dR/dqrG0ZEqozCvZrlHilh8oJkPlufwYCoCF67rj9RLRo4XZaIuIzCvRr9mHqAR+YmkJVfxIMXd+e+C7tRO8yfgTlFRCqmcK8GBcWlvPRVCh/+uIMuLRsy/56z6dsxwumyRMTFFO4Btmr7QcbHJ7DzlwJuPTuaCcNiqF9XIzmKSGAp3APkUFEpL3+dwkc/7aBDs/rMGjeYwV1aOF2WiIQIhXsALE3J4slPk8nIPcItgzvx2LAYGuoSRxGpRkqcKpSVV8jzX2xiUUIG3Vo3Iv7usxjYqbnTZYlICFK4V4Eyj+XvP+3gL99soajUwwNDu/P7C7tqliQRcYzC3U8r0n7h2UUb2ZSZx3ndW/LcyF50btnQ6bJEJMQp3Ctp98ECpvxzM18kZdI+oj7TbhzA8N5tdZepiNQICvdTlH24mGlLU/n4p53UqgUPXdyDced30eWNIlKjKNxP0qGiUj76cQdvf7+Nw0WljB7QgYd/04PIpvWdLk1E5L8o3E/gUFEpH/+0g3d/SCO7oIShMa15bFiMJqkWkRpN4X4cWXmFfPjjDmas2EleYSkXntaKBy/uoWEDRCQoKNzLsdaSkJ7LjBU7Wbg+gxKPh2E923LXBV3pp1AXkSCicMf7IekXSZnMXLWLDRl5NKgbxrVndODOc7sQrcsaRSQIhWy45xQUsyxlP4sSMvh+y35KPZaYto15flQvRvVrR+NwTUwtIsErYOFujBkGvA6EAe9Za6cEal8no6i0jKT0XFZuP8j3KftZvfMgHgttm4Rz+7mdGdmvHbGRTXSduoi4QkDC3RgTBkwDLgHSgZ+NMQuttRsDsb/yDheVsi+vkMzcQrbuy2dL1iFS9uaTtCeX4lIPALGRTbj3wm5cGNOafh0iqFVLgS4i7hKolvuZQKq1Ng3AGDMLGAlUabgvS8ni+c83Uljioai0jIJi71d5TevXoUebRtx6djRxnZoxsFMzWjSqV5VliIjUOIEK9/bA7nLL6cCg8hsYY8YB4wCioqIqtZMm9esQ07YJ9erUIrxOGPXrhNGqcT3aNKlHm8bhdGvdiFaN66mrRURCTqDCvaI0tf+xYO10YDpAXFycrWD7ExoQ1YwBY5pV5qUiIq4WqNmZ04GO5ZY7ABkB2peIiBwjUOH+M9DdGNPZGFMXuB5YGKB9iYjIMQLSLWOtLTXG3Ad8jfdSyA+stRsCsS8REflvAbvO3Vr7JfBloH6+iIgcX6C6ZURExEEKdxERF1K4i4i4kMJdRMSFjLWVun+oaoswZj+ws5IvbwkcqMJygkUoHncoHjOE5nGH4jHDqR93J2ttq4qeqBHh7g9jzGprbZzTdVS3UDzuUDxmCM3jDsVjhqo9bnXLiIi4kMJdRMSF3BDu050uwCGheNyheMwQmscdiscMVXjcQd/nLiIi/80NLXcRETmGwl1ExIWCOtyNMcOMMSnGmFRjzESn6wkEY0xHY8xSY8wmY8wGY8wDvvXNjTHfGmO2+r67ctYSY0yYMWadMeZz33JnY8xK33HP9g0p7RrGmAhjTLwxZrPvnJ8VCufaGPOQ7/2dbIyZaYwJd+O5NsZ8YIzJMsYkl1tX4fk1Xn/15VuiMWbAqewraMO93CTclwKxwA3GmFhnqwqIUuARa+3pwGDgXt9xTgSWWGu7A0t8y270ALCp3PJU4FXfcWcDdzhSVeC8DnxlrY0B+uI9dlefa2NMe+B+IM5a2wvvMOHX485z/SEw7Jh1xzu/lwLdfV/jgLdOZUdBG+6Um4TbWlsMHJ2E21WstZnW2rW+x/l4/9jb4z3Wj3ybfQSMcqbCwDHGdAAuA97zLRvgIiDet4mrjtsY0wQ4H3gfwFpbbK3NIQTONd7hx+sbY2oDDYBMXHiurbU/AAePWX288zsS+Nh6rQAijDGRJ7uvYA73iibhbu9QLdXCGBMN9AdWAm2stZng/QcAaO1cZQHzGvAY4PEttwByrLWlvmW3nfMuwH7gb76uqPeMMQ1x+bm21u4BXgZ24Q31XGAN7j7X5R3v/PqVccEc7kREnNsAAAGvSURBVCechNtNjDGNgHnAg9baPKfrCTRjzAggy1q7pvzqCjZ10zmvDQwA3rLW9gcO47IumIr4+phHAp2BdkBDvF0Sx3LTuT4Zfr3fgzncQ2YSbmNMHbzB/om1dr5v9b6j/0Xzfc9yqr4AOQe4whizA2+X20V4W/IRvv+6g/vOeTqQbq1d6VuOxxv2bj/XFwPbrbX7rbUlwHzgbNx9rss73vn1K+OCOdxDYhJuXz/z+8Ama+0r5Z5aCIz1PR4LLKju2gLJWjvJWtvBWhuN99x+Z60dAywFrvZt5qrjttbuBXYbY07zrRoKbMTl5xpvd8xgY0wD3/v96HG79lwf43jndyFwi++qmcFA7tHum5NirQ3aL2A4sAXYBjzhdD0BOsZz8f5XLBFY7/sajrf/eQmw1fe9udO1BvB3MAT43Pe4C7AKSAXmAvWcrq+Kj7UfsNp3vj8DmoXCuQaeBTYDycDfgXpuPNfATLyfK5TgbZnfcbzzi7dbZpov35LwXk100vvS8AMiIi4UzN0yIiJyHAp3EREXUriLiLiQwl1ExIUU7iIiLqRwFxFxIYW7iIgL/T8ht9I/OrMo2QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(traj[:, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a26938550>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD5CAYAAAAk7Y4VAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhU933v8fd3NNqFkJDYhRCLbDBeWGSMjfEG1GsMXhIntRuc2iXN0qTLTZtc52nv7XIft07TNDduUy52TBwvcbxig41tvGODAWNL7IvYhDYkARIS2mZ+948ZXEwkJJiRzmjm83qeeeacmaP5fY8P/ujoO2cx5xwiIhL/fF4XICIi/UOBLyKSIBT4IiIJQoEvIpIgFPgiIglCgS8ikiD80fgQM7sB+HcgCVjqnHvwtPf/Ergf6AQOA3/snNt/ps/Mz893RUVF0ShPRCRhbNy4sc45N7Sr9yIOfDNLAh4G5gMVwHozW+6c23rKYpuAEudci5l9C/gX4K4zfW5RUREbNmyItDwRkYRiZt3uTEejpTMT2O2cK3fOtQNPAwtOXcA597ZzriU8uxYoiMK4IiJyFqIR+KOBg6fMV4Rf6859wKtRGFdERM5CNHr41sVrXV6vwczuAUqAq7t5fzGwGKCwsDAKpYmIyEnR2MOvAMacMl8AVJ6+kJnNAx4AbnXOtXX1Qc65Jc65EudcydChXX7nICIi5ygagb8eKDazcWaWAnwVWH7qAmY2DfgvQmFfG4UxRUTkLEUc+M65TuC7wCpgG/CMc26Lmf29md0aXuwhIAv4nZl9ambLu/k4ERHpI1E5Dt85txJYedprf3vK9LxojCMiIucuKoEvcrY6AkEOHTnBvvpmqo61cqI9wImOAB2BIJkpfjJT/WSn+ynIzWDskAxyMpIx6+r4ABHpLQW+9Ivmtk7WlteHHw1srWokEOz9zXey0/xMGTWYiwsGc8mYHC4bN4S8rNQ+rFgk/ijwpc8Ego615fU890kFr22upqU9QIrfx/TCHL551XjG5WdSlJ/JqJx0MlOSSE9Jwu/z0dLeyfG2To62dFBx5AT765spr2tm86FjPLpmLx2B0C+KKaOymVM8lPkXDGPamFx8Pv0FIHImFqu3OCwpKXG6tMLA1BkIsvyzSn7x1m7K65oZlOrnlktG8qWLRzF9bC5pyUnn/NltnQG2VDby4e463ttVx6YDR+gIOIZnp3L9lBHceskoZozNVftHEpaZbXTOlXT5ngJfosU5x8ulVfz09R3sq29h0ohBfOuaCVw/ZUREIX8mTa0dvLW9llfLqnlnZy2tHUHG5mVw+7QC7iwpYHROep+MKxKrFPjS5w7Ut/DAi2W8v6uOySOz+fN5xcyfPLxf2yzH2zp5bXM1z22s4KPyenwG100axt2XjeWq84aSpJaPJAAFvvSZYNDx6Jq9/OT1Hfh9Pn5w/fncM2us5+F6sKGF364/yNPrD1J3vI2xeRnce0URXy4ZQ1aqvrqS+KXAlz5x7EQHf/XMZ7y5rYZ5k4fzDwunMHJwbLVQ2juDvL61ml+t2cfG/UcYlOrnqzPHcN+V4xkxOM3r8kSiToEvUbe1spFvPbGRQ0dO8MDNk7n3iqKY/6L004NHeeSDvawsq8JnsGDqaP706vFMHDbI69JEokaBL1H13s7DfPPxjWSn+/mPu6czY+wQr0s6KwcbWlj6fjm/3XCQts4gN0wZwXeunciFowd7XZpIxBT4EjWvlFbyF7/9lInDBrHsjy9l2KCB2xZpaG7nV2v28tiH+2hq7eTa84fyvbnFTCvM9bo0kXOmwJeoeGLdfn784mZKxuaydNGlDE5P9rqkqGhs7eDxj/az9P1yjrR0MKc4n+/PLaakaGD95SICCnyJgmfWH+SvnyvluknDePgPp5Oe0jfH1Xupua2T36zdz5L3yqlvbmf2xDz+fN55XKrglwFEgS8ReW1zFd9+4hNmT8znkUWXkuKPxm0UYteJ9gBPrNvPL98tp+54G1dMyOP7c4u5bHye16WJ9EiBL+dsze46vvGr9UwZnc0T919GRkriHMN+Mvj/671yDje1MWv8EL4/9zwun6Dgl9ilwJdzsqO6idv/Yw0FuRn89puzyMlI8bokT7R2BHhy3QH+8909HG5qY2bREL43t5jZE/Ni/lBUSTwKfDlrR1vaufUXazjREeDl716pk5QIBf/THx/gl++WU93YyvTCHL573USuPX+Ygl9ixpkCP76bsXJOOgNB/uypTVQdO8Ev75mhsA9LS07i3tnjePevr+EfF15ITWMbf/zYBm7++QesKK06q+v7i3hBgS+/56FVO3h/Vx3/sOBCZozVMemnS/Uncc+ssbzzg2t46M6Lae0I8J0nP2HeT9/lqY8P0NYZ8LpEkS6ppSNf8PqWahY/vpF7ZhXyjwsv8rqcASEQdLy2uZpfvruHskPHGDYolXtnF3H3zLEMzoiPcxVk4FAPX3qltrGV63/2HqNz03n+W7Pj/vDLaHPOsWZ3Pf/13h7e31VHRkoSXykZwzdmFzE2L9Pr8iRBnCnwE+cYOzmjYNDxV7/7jBMdAX521zSF/TkwM64szufK4ny2Vjay9INynli3n2Uf7WPupGHce8U4HdkjnlLgCwCPrtnL+7vq+D+3XcTEYVlelzPgXTAqm59+ZSp/c8Mknli7nyfWHeDNbeuYMDSTuy8byx0zCuLm0hQycKilI+ysaeKWn3/A1ecPZckfzdAeaB9o7QiworSK36zbz6YDR0lL9nHzRaO469IxXFqke/BK9KilI90KBB1/81wpWWl+Hrz9IgVPH0lLTuKOGQXcMaOAzYeO8eTHB1j+aSXPfVLB+PxM7phRwMJpo3UPXulTUdnDN7MbgH8HkoClzrkHT3s/Ffg1MAOoB+5yzu0702dqD79/LPtwH3+3fAs/u2sqC6eN9rqchNLS3smK0iqe2XCQ9fuOAHDZuCEsnDaa66eMYEhmYp7ZLJHp06N0zCwJ2AnMByqA9cDXnHNbT1nm28DFzrk/NbOvArc55+460+cq8Pte5dETzP/pu8woGsKyb1yqvXsPHahv4aVPD/HCpkOU1zWT5DOumJDHTReNZO7kYQP6vgPSv/q6pTMT2O2cKw8P9jSwANh6yjILgP8Vnn4W+IWZmYvVLxASgHOOv31pM0EH/7TwQoW9xwrzMvizucV897qJbK1qZEVpFSvKqvjR82WYwdQxOcybPJyrzxvKlFHZ2l5xpDMQ5ERHgLbOYOjRESDJZ31yKG80An80cPCU+Qrgsu6Wcc51mtkxIA+oi8L4cg5WbanhzW21/PjmyYwZkuF1ORJmZkwZNZgpowbzg+vPZ0dNE69vqeGNrTU8tGoHD63awdBBqcyZmM/lE/K4fEIeBbnafv0pGHQ0tXXSeKKDxtYOmlo7w48OjreFpo+3ddLcFnpuaQvQ3N5JS3sg/OjkRHsg9OgI0NnFJTmmFebwwrdnR732aAR+V7sap69Bb5bBzBYDiwEKCwsjr0y61NoR4J9WbuX84YO494oir8uRbpgZk0ZkM2lENt+bW0xtUyvv7azjnR21vLPzMM9vOgRAQW46M8bmMmNsLtMLczlv+CCdR9ELwaCjqbWThpZ2jrS0c6S5nSMtHRxtaedoSwdHT4Tmj4Wnj50ITTe1ddJTb8LvMzJT/WSl+slMTSIz1U9mip8hmSlkpCSRkZJEqj/0nJacRFqyj/Tk0GupyT6GZqX2yTpHI/ArgDGnzBcAld0sU2FmfmAw0HD6BznnlgBLINTDj0Jt0oVH1+zlYMMJnrj/MvxJCoaBYtigNO6cUcCdMwoIBh07a5tYu6eedXsb+GhPPS99GvrfLiXJx6SRg7hw9GAmjRjE+cMHcf6IQXF/eeu2zgBHmjtoaG6nobmd+uY2jnw+HQr1+uOh54ZwuHd3wTufQU5GCjkZyeSkJzM0K5WJQ7MYnJ7M4PRksk8+0vwMSktm0CnPWal+Uv2+mGy7RSPw1wPFZjYOOAR8FfjD05ZZDiwCPgLuBN5S/94btY2tPPzWbuZfMJzZE/O9LkfOkc/333v/984eh3OOQ0dPsOnAUTYfOkbZoWO88lklT67r/PxnhmSmMC4/k6K8TMYMSacgN4OC3HRGZKcxPDstZm5b6ZzjeFtnaI86vFd9pKWDIy3tHG1p/3z6SHM7DS0dNDS3caQ51E7pihnkZqSQm5FMXmYq4/IzmTE2lyGZKeRmpISew9O5GcnkZKQwKNWPzxd7gR2piAM/3JP/LrCK0GGZjzrntpjZ3wMbnHPLgUeAx81sN6E9+69GOq6cm4dW7aA9EOSBmyZ7XYpEkZmFAzyDL10yCggFZ01jG9urG9lVc5zyumb21h3ng92HqW1q+722xKBUP3lZofDLy0wh+5Q914zUJNKTQ4/UZB9+n4/kJB9JPsMIhSpA0EHQOQJBR0cgSEfA0d4ZpLUjQGtngNaOICfaO2luD9DS1snxtgDH20Jh3Xgi1AdvbO0846WmM1KSvhDURXkZDAnXfLL2k+8PyUwhJyOFpDgM73OhM20TSFnFMW59+AMWzxnPjxT4Ca2tM0DV0VYqjpygprGV2qY2ahpbP2+HNDS3n/KFZAfRvNT/yR52Rkqo/ZEVboOcbI9kp/s/b50MTg/tceeebK9kJJPqj42/RGKVzrQVAP5l1XZy0pP57nUTvS5FPJbqT6IoP5Oi/J4P/XPO0RbeS29pD4T33EOHEAaD4HCf/7WQ5DPMQs/JST6SfT5S/D7Skn2hLyT9vrhslQwUCvwE8eGeOt7fVcePb57MoDRdtEt6z8zCR5IkkaMjQAc0HaKRAJxz/GTVDkYOTuOeWWO9LkdEPKLATwCrt9XyyYGjfG9uMWnJ6n+KJCoFfpwLBh0/eX0HRXkZ3DmjwOtyRMRDCvw490pZFdurm/jLPzifZJ1kJZLQlABxLBh0/OKtXZw3PItbLhrpdTki4jEFfhx7fWsNO2uO851rJ+pQOBFR4Mcr5xwPv72borwMbtbevYigwI9b7+48TNmhY3zrmgm6QJqIAAr8uPXw27sZNTiN26bpyBwRCVHgx6F15fWs33eEb149QddFF5HPKQ3i0C/f3UN+Vgp3XTqm54VFJGEo8OPM7tom3t5xmD+aVaSzakXkCxT4ceaRD/aS6vdxzyzdIlJEvkiBH0fqj7fx3CeHuH16AXl9dE9MERm4FPhx5PG1+2nvDHLfleO8LkVEYpACP060dgR4/KP9zJ00jInDsrwuR0RikAI/Try46RD1ze3cN0d79yLSNQV+HHDO8diH+7hgZDaXj8/zuhwRiVEK/Diwft8Rtlc38fXLx2Kmi6SJSNcU+HHg1x/tIzvNz4Kpo70uRURimAJ/gKttbOW1zdV8pWQM6Sk60UpEuqfAH+Ce/PgAnUGnm5OLSI8U+ANYRyDIk+sOcPV5QynKz/S6HBGJcREFvpkNMbM3zGxX+Dm3i2WmmtlHZrbFzErN7K5IxpT/9vqWGmqb2vj65dq7F5GeRbqH/0NgtXOuGFgdnj9dC/B159wU4AbgZ2aWE+G4Ajy+dh8Fuelcc/4wr0sRkQEg0sBfACwLTy8DFp6+gHNup3NuV3i6EqgFhkY4bsLbW9fM2vIGvjazkCTdr1ZEeiHSwB/unKsCCD+fcVfTzGYCKcCeCMdNeE+vP0CSz/jyDN3RSkR6x9/TAmb2JjCii7ceOJuBzGwk8DiwyDkX7GaZxcBigMJCXd63O+2dQZ7bWMF1k4YxLDvN63JEZIDoMfCdc/O6e8/MasxspHOuKhzotd0slw2sAH7snFt7hrGWAEsASkpKXE+1JarV22qoO97OH87UL0UR6b1IWzrLgUXh6UXAS6cvYGYpwAvAr51zv4twPAGeWn+QUYPTuOo8fRUiIr0XaeA/CMw3s13A/PA8ZlZiZkvDy3wFuAq418w+DT+mRjhuwjrY0ML7uw7z5ZIx+rJWRM5Kjy2dM3HO1QNzu3h9A3B/ePo3wG8iGUf+2zMbDgLwFd2gXETOks60HUACQcezGyu4qngoo3PSvS5HRAYYBf4A8uGeOqqOtfKVEu3di8jZU+APIM9urGBwejJzJ+vMWhE5ewr8AaKxtYPXNldz6yWjSEvWZZBF5Owp8AeIFaVVtHUGuVNn1orIOVLgDxDPbqzgvOFZXFww2OtSRGSAUuAPAHsOH2fj/iPcOaNA96wVkXOmwB8AnttYQZLPWKh71opIBBT4MS4YdLyw6RBXFefrQmkiEhEFfoxbu7eeqmOt3D5dX9aKSGQU+DHuxU2HyEr1M2/ycK9LEZEBToEfw1o7ArxaVs0NF44gPUXH3otIZBT4MWz1tlqa2jq5bZq+rBWRyCnwY9gLmyoYkZ3GrPF5XpciInFAgR+jGprbeWfHYRZMHaXr3otIVCjwY9SK0ko6g46FaueISJQo8GPUC5sOMWnEICaPzPa6FBGJEwr8GHSwoYVPDhxlgc6sFZEoUuDHoOWfVQLwpUtGelyJiMQTBX4MevmzSmaMzaUgN8PrUkQkjijwY8zOmia2Vzdx6yWjvC5FROKMAj/GLP+0Ep/BTRepnSMi0aXAjyHOOZZ/VsnsifkMHZTqdTkiEmcU+DHk04NHOdDQwpfUzhGRPqDAjyHLP6skxe/jhgtHeF2KiMShiALfzIaY2Rtmtiv8nHuGZbPN7JCZ/SKSMeNVIOhYUVrFtecPJTst2etyRCQORbqH/0NgtXOuGFgdnu/OPwDvRjhe3Fq/r4HapjZuuVjtHBHpG5EG/gJgWXh6GbCwq4XMbAYwHHg9wvHi1iullaQnJzF38jCvSxGROBVp4A93zlUBhJ9/L63MzAf8K/CDCMeKW52BIK9trua6ycPISPF7XY6IxKke08XM3gS6+hbxgV6O8W1gpXPuoNmZL/NrZouBxQCFhYW9/PiBb93eBuqOt3OLjr0XkT7UY+A75+Z1956Z1ZjZSOdclZmNBGq7WOxyYI6ZfRvIAlLM7Lhz7vf6/c65JcASgJKSEtfblRjoXimtJDMliWsnqZ0jIn0n0v7BcmAR8GD4+aXTF3DO3X1y2szuBUq6CvtE1RFu58y7YDhpybpvrYj0nUh7+A8C881sFzA/PI+ZlZjZ0kiLSwQf7qnnSEsHN6udIyJ9LKI9fOdcPTC3i9c3APd38fpjwGORjBlvVpRWMijVz1XnDfW6FBGJczrT1kPtnUFWbalhvto5ItIPFPgeWrOnjmMnOrj5YrVzRKTvKfA9tLK0ikGpfq4szve6FBFJAAp8j3QEgry+NdTOSfWrnSMifU+B75E1u0PtHN3oRET6iwLfIyvLQu2cOeepnSMi/UOB74GT7Zx5aueISD9S4Hvgwz31HG1RO0dE+pcC3wMrS6vISvUzR0fniEg/UuD3s45AkFVbq5k3eZhOthKRfqXA72dry0PtnBvVzhGRfqbA72cry6rJTEnial07R0T6mQK/H3UGgqzaUs11k3XtHBHpfwr8fvTx3gYamtu56cKubiAmItK3FPj9aOXmKtKTk7jmfN3ZSkT6nwK/nwSCjtc213DtpKGkp6idIyL9T4HfTzbsa6DueBs3Xqijc0TEGwr8fvLq5mpS/T7dqFxEPKPA7wfBoOPVzVVcfd5QslIjvW+8iMi5UeD3g00Hj1DT2KZr54iIpxT4/WBlWTUpST6um6x2joh4R4Hfx5xzvFpWxZzifLLTkr0uR0QSmAK/j31WcYzKY626do6IeE6B38deLavC7zPmTx7udSkikuAU+H3IOcfKzVXMnpjP4Ay1c0TEWxEFvpkNMbM3zGxX+Dm3m+UKzex1M9tmZlvNrCiScQeKLZWNHGw4wU0X6do5IuK9SPfwfwisds4VA6vD8135NfCQc24yMBOojXDcAWFlWRVJPmP+BQp8EfFepIG/AFgWnl4GLDx9ATO7APA7594AcM4dd861RDhuzHPOsbKsisvH5zEkM8XrckREIg784c65KoDwc1cHmp8HHDWz581sk5k9ZGZxf/Ww7dVN7Ktv4Ua1c0QkRvR4nr+ZvQl0lVoPnMUYc4BpwAHgt8C9wCNdjLUYWAxQWFjYy4+PTSvLqvAZ/IHaOSISI3oMfOfcvO7eM7MaMxvpnKsys5F03ZuvADY558rDP/MiMIsuAt85twRYAlBSUuJ6twqxxznHirIqLhuXx9BBqV6XIyICRN7SWQ4sCk8vAl7qYpn1QK6ZnbyJ63XA1gjHjWk7apooP9zMTRfrZCsRiR2RBv6DwHwz2wXMD89jZiVmthTAORcA/gew2szKAAP+X4TjxrSVpaF2zg1T1M4RkdgR0bV6nXP1wNwuXt8A3H/K/BvAxZGMNVCcbOfMHDdE7RwRiSk60zbKdtYcZ8/hZm7WtXNEJMYo8KNsRVkVZnD9hWrniEhsUeBH2cqyKi4tGsKwQWlelyIi8gUK/CjaVdPE7trjaueISExS4EfRK6Whds6NaueISAxS4EeJc45XSiuZWTSEYdlq54hI7FHgR8n26ib2HG7mlktGeV2KiEiXFPhR8kppJT61c0QkhinwoyDUzqniign55GfpZCsRiU0K/CjYUtnI/voWbtG1c0Qkhinwo+Dl0kr8PuN6XTtHRGKYAj9CzjlWlIZuVJ6rO1uJSAxT4Efo04NHqThyQu0cEYl5CvwIvfxZFSlJPv5A7RwRiXEK/AgEgo6XSyu5dtJQBqcne12OiMgZKfAjsLa8nsNNbSyYOtrrUkREeqTAj8CLmw6RlernuknDvC5FRKRHCvxz1NoR4LXN1Vw/ZQRpyUlelyMi0iMF/jl6Z0ctTW2dLJiqa+eIyMCgwD9HL31aSX5WKldMyPO6FBGRXlHgn4PG1g5Wb6/llotH4k/Sf0IRGRiUVufgtc3VtHcG1c4RkQFFgX8Onv+kgqK8DKaOyfG6FBGRXlPgn6WDDS2sLW/gjukFmJnX5YiI9JoC/yy9sOkQALdN18lWIjKwRBT4ZjbEzN4ws13h59xulvsXM9tiZtvM7Oc2QHeNnXM890kFl4/PoyA3w+tyRETOSqR7+D8EVjvnioHV4fkvMLMrgNnAxcCFwKXA1RGO64kN+4+wv76FO2YUeF2KiMhZizTwFwDLwtPLgIVdLOOANCAFSAWSgZoIx/XEcxsryEhJ0n1rRWRAijTwhzvnqgDCz793URnn3EfA20BV+LHKObctwnH7XWtHgBWlVdx44UgyU/1elyMictZ6TC4zexPoapf2gd4MYGYTgcnAyT7IG2Z2lXPuvS6WXQwsBigsLOzNx/ebVVuqaWrr5I4Z+rJWRAamHgPfOTevu/fMrMbMRjrnqsxsJFDbxWK3AWudc8fDP/MqMAv4vcB3zi0BlgCUlJS43q1C/3hmw0FG56Qza5wupSAiA1OkLZ3lwKLw9CLgpS6WOQBcbWZ+M0sm9IXtgGrp7KtrZs3uer42cww+34A8wEhEJOLAfxCYb2a7gPnhecysxMyWhpd5FtgDlAGfAZ85516OcNx+9dT6AyT5jC+XjPG6FBGRcxbRt4/OuXpgbhevbwDuD08HgG9GMo6X2juDPLuhgnmThzE8O83rckREzpnOtO3B61urqW9u52szY+tLZBGRs6XA78GT6w4wOiedq4qHel2KiEhEFPhnsLeumQ/36MtaEYkPCvwzePpjfVkrIvFDgd+N5rZOnvr4ANdPGa4va0UkLijwu/HcJxU0tnZy35XjvC5FRCQqFPhdCAYdv1qzj6ljcphe2OUVn0VEBhwFfhfe2l7L3rpm7rtynO5qJSJxQ4HfhaUflDNqcJougywicUWBf5otlcdYW97AvbOL8CfpP4+IxA8l2mmWvr+XjJQk7rpUZ9aKSHxR4J+i/PBxXvr0EHdfVsjg9GSvyxERiSoF/il+vnoXqf4kvnn1BK9LERGJOgV+2O7aJl76rJJFVxSRn5XqdTkiIlGnwA/72Zu7yEhOYvFV470uRUSkTyjwge3Vjawoq+Ibs8cxJDPF63JERPqEAh/419d3kpXi5/45uoyCiMSvhA/8t7fX8sbWGv70mgnkZGjvXkTiV0IH/on2AH+7fDMTh2XxJ3PUuxeR+BbRPW0Hul+8vYuDDSd46k9mkeJP6N99IpIAEjbldtc2seS9cm6fPprLJ+R5XY6ISJ9LyMDvCAT50fNlZKT4+Z83Tfa6HBGRfpGQLZ1/WrGN9fuO8G93XaKTrEQkYSTcHv4z6w/y2If7uO/Kcdw2rcDrckRE+k1CBf7G/Uf48YubmVOcz49unOR1OSIi/SqiwDezL5vZFjMLmlnJGZa7wcx2mNluM/thJGOeqzW767h/2XpG5qTxf782Tde6F5GEE2nqbQZuB97rbgEzSwIeBm4ELgC+ZmYXRDhurznn+I93dvNHj6wjLyuVx74xUydYiUhCiuhLW+fcNqCn+77OBHY758rDyz4NLAC2RjJ2T2qbWvlgVx0vbDrE+7vquOXikfzzHReTmZqQ31OLiPTLUTqjgYOnzFcAl/XVYIeOnuC+x9azvboJgCGZKfzdly7g3iuKdENyEUloPQa+mb0JdHU37weccy/1YoyuUtZ1M9ZiYDFAYeG53WJw+KBURuWks2DqaOYU53PByGx8PgW9iEiPge+cmxfhGBXAmFPmC4DKbsZaAiwBKCkp6fKXQk/8ST4evffSc/lREZG41h+HqqwHis1snJmlAF8FlvfDuCIicopID8u8zcwqgMuBFWa2Kvz6KDNbCeCc6wS+C6wCtgHPOOe2RFa2iIicrUiP0nkBeKGL1yuBm06ZXwmsjGQsERGJjM4+EhFJEAp8EZEEocAXEUkQCnwRkQShwBcRSRDm3Dmd39TnzOwwsD+Cj8gH6qJUzkCRiOsMibneibjOkJjrfbbrPNY5N7SrN2I28CNlZhucc91esjkeJeI6Q2KudyKuMyTmekdzndXSERFJEAp8EZEEEc+Bv8TrAjyQiOsMibneibjOkJjrHbV1jtsevoiIfFE87+GLiMgp4i7wY+GG6f3BzMaY2dtmti18I/nvh18fYmZvmNmu8HOu17VGm5klmdkmM3slPD/OzNaF1/m34ctwxxUzyzGzZ81se3ibXx7v29rM/iL8b3uzmT1lZmnxuK3N7FEzqzWzzae81uW2tZCfh/20U1UAAALySURBVPOt1Mymn81YcRX4Xt8wvZ91An/lnJsMzAK+E17XHwKrnXPFwOrwfLz5PqFLbZ/0z8C/hdf5CHCfJ1X1rX8HXnPOTQIuIbT+cbutzWw08D2gxDl3IZBE6F4a8bitHwNuOO217rbtjUBx+LEY+M+zGSiuAp9TbpjunGsHTt4wPe4456qcc5+Ep5sIBcBoQuu7LLzYMmChNxX2DTMrAG4GlobnDbgOeDa8SDyuczZwFfAIgHOu3Tl3lDjf1oQu355uZn4gA6giDre1c+49oOG0l7vbtguAX7uQtUCOmY3s7VjxFvhd3TB9tEe19BszKwKmAeuA4c65Kgj9UgCGeVdZn/gZ8NdAMDyfBxwN32gH4nObjwcOA78Kt7KWmlkmcbytnXOHgJ8ABwgF/TFgI/G/rU/qbttGlHHxFvi9vmF6vDCzLOA54M+dc41e19OXzOwWoNY5t/HUl7tYNN62uR+YDvync24a0EwctW+6Eu5ZLwDGAaOATELtjNPF27buSUT/3uMt8Ht9w/R4YGbJhML+Cefc8+GXa07+iRd+rvWqvj4wG7jVzPYRatddR2iPPyf8Zz/E5zavACqcc+vC888S+gUQz9t6HrDXOXfYOdcBPA9cQfxv65O627YRZVy8BX7C3DA93Lt+BNjmnPvpKW8tBxaFpxcBL/V3bX3FOfcj51yBc66I0LZ9yzl3N/A2cGd4sbhaZwDnXDVw0MzOD780F9hKHG9rQq2cWWaWEf63fnKd43pbn6K7bbsc+Hr4aJ1ZwLGTrZ9ecc7F1YPQvXR3AnuAB7yupw/X80pCf8qVAp+GHzcR6mmvBnaFn4d4XWsfrf81wCvh6fHAx8Bu4HdAqtf19cH6TgU2hLf3i0BuvG9r4H8D24HNwONAajxua+ApQt9TdBDag7+vu21LqKXzcDjfyggdxdTrsXSmrYhIgoi3lo6IiHRDgS8ikiAU+CIiCUKBLyKSIBT4IiIJQoEvIpIgFPgiIglCgS8ikiD+P5tV9Juc5/bTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot([angdiff(th, thetades(0.0)) for th in traj[:, 1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a25dadf90>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAfo0lEQVR4nO3deXhc9X3v8fd3Fo1Wy9ply7LlDe8GgwxmX8xOAmQhl4QkkKX03pubpclNkzTP0970tmnz3LZp2pubPNQkJWkakgINDoSwGJtAAINsY7zvmyxbi2XZ1jqa0e/+MWMwRMayNaMzZ+bzep55Zs7Mked7fOzPHH3nd87PnHOIiIj/BLwuQEREzo0CXETEpxTgIiI+pQAXEfEpBbiIiE+FxvLNKisrXUNDw1i+pYiI761Zs6bDOVf17ufHNMAbGhpoamoay7cUEfE9M9s33PNqoYiI+JQCXETEpxTgIiI+pQAXEfEpBbiIiE8pwEVEfEoBLiLiU2M6DvxcHe2JcuBoLy1dfTQf7WNqZRFL59R4XZaIiKdGFOBm9ifAZwEHbAA+BUwAHgbKgbXAJ5xz0XQU+YWH1/Hijo53PPftDyzgY5dMTsfbiYj4whkD3MzqgC8Ac51zfWb2S+Bu4Fbgu865h83sh8BngB+ko8j/evV0PrFkCnVlBVSX5POnj6znm7/aQFEkyB0X1KXjLUVEMt5Ie+AhoMDMQkAhcAi4Dngk+fpDwJ2pLy/h8hmV3DivlnkTS6kqifCDj1/EJVPL+fIv1/PMpsPpelsRkYx2xgB3zh0E/g7YTyK4jwFrgC7nXCy5WjMw7KGwmd1vZk1m1tTe3p6SovPDQZbdu5j5daV86RdvcKR7ICV/roiIn5wxwM2sDLgDmApMBIqAW4ZZddjJNZ1zDzjnGp1zjVVVf3AxrXNWHAnx93edT99gnAde3J2yP1dExC9G0kK5HtjjnGt3zg0CjwGXAeOTLRWASUBLmmo8rRnVxdx+/kR++so+HYWLSM4ZSYDvB5aYWaGZGbAU2AysBD6cXOde4PH0lPjePn/dDPoG4/zLi3u8eHsREc+MpAe+msSXlWtJDCEMAA8AXwO+bGY7gQrgwTTWeVozqkt4/8KJ/OSVvXT2pGUUo4hIRhrRKBTn3F8452Y75+Y75z7hnBtwzu12zl3snJvhnLvLOedZD+MLS08ehasXLiK5IytOpZ9RXcL7Fk7kJy/vpXsgduYfEBHJAlkR4AD3XdZATzTOUxsOeV2KiMiYyJoAv3DyeBoqCnls7UGvSxERGRNZE+BmxgcvnMQru4/QfLTX63JERNIuawIc4AOLEieD/mqdjsJFJPtlVYDXlxdyydRyHl17EOeGPTFURCRrZFWAA3zookns6ehh3YEur0sREUmrrAvwW+bXkh8O8OiaZq9LERFJq6wL8JL8MDfNq+XX61sYiMW9LkdEJG2yLsAB7rygjuP9MV7eecTrUkRE0iYrA/yyGRUU5QV5ZrMmexCR7JWVAR4JBblmdjXPbm4lPqTRKCKSnbIywAFumldLR3eUtfuPel2KiEhaZG2AXzurirxggKc3qo0iItkpawO8JD/MZTMqeGZzq07qEZGslLUBDok2yv7OXrYePuF1KSIiKZfVAX79nBrM4OlNaqOISPbJ6gCvKolw0eQynt7U6nUpIiIpl9UBDok2ypZDxznQqUvMikh2yfoAv35uDQArt7V5XImISGplfYBPrSyioaKQ57cqwEUku2R9gANcO7uaV3YdoS+qi1uJSPbIiQC/bnY1A7EhXtnd4XUpIiIpkxMBfvHUcgrzgmqjiEhWyYkAj4SCXD6jkpVb23VWpohkjZwIcIBrZ1VzsKuPHW3dXpciIpISuRPgs6sAWKk2iohkiZwJ8AmlBcyZME59cBHJGjkT4JC4xGzTvqMc6xv0uhQRkVHLqQC/bnY18SHHSzs0nFBE/C+nAvyC+vGMyw/xwna1UUTE/3IqwEPBAFfOrOKF7RpOKCL+l1MBDnD1eVW0Hh9gW6smeRARf8u9AJ+VGE64alu7x5WIiIxOzgV4zbh8ZteW8IICXER8bkQBbmbjzewRM9tqZlvM7FIzKzezZ81sR/K+LN3Fpso1s6pp2tdJ90DM61JERM7ZSI/Avwf81jk3Gzgf2AJ8HVjhnJsJrEgu+8LV51UxGHe8vFPDCUXEv84Y4GY2DrgKeBDAORd1znUBdwAPJVd7CLgzXUWm2kVTyijKC/LCdrVRRMS/RnIEPg1oB35sZuvMbJmZFQE1zrlDAMn76uF+2MzuN7MmM2tqb8+MwMwLBbh8RiWrtmk4oYj410gCPARcCPzAObcI6OEs2iXOuQecc43OucaqqqpzLDP1rp5VxcGuPna16+qEIuJPIwnwZqDZObc6ufwIiUBvNbMJAMl7X53eePV5Gk4oIv52xgB3zh0GDpjZrORTS4HNwHLg3uRz9wKPp6XCNJlUVsiM6mL1wUXEt0IjXO/zwM/MLA/YDXyKRPj/0sw+A+wH7kpPielzzXlV/OTVffRF4xTkBb0uR0TkrIwowJ1zbwCNw7y0NLXljK2rZ1Wx7KU9vLr7CNfOHvY7WBGRjJVzZ2KeanFDOQXhIKu2+ap9LyIC5HiA54eDXDq9glXqg4uID+V0gANcM6uKfUd62dvR43UpIiJnJecD/O3hhGqjiIi/5HyAT6koYmplkYYTiojv5HyAQ+Io/JXdR+gfjHtdiojIiCnASQwn7B8cYvWeTq9LEREZMQU4cOm0CiKhgPrgIuIrCnBOGU6o66KIiI8owJOunVXNno4e9mg4oYj4hAI86dpZiVPpV25VG0VE/EEBnjS5opDpVUWsVB9cRHxCAX6Ka2dVs3p3J71RTXYsIplPAX6Ka2dXE40P8fLOI16XIiJyRgrwUzQ2JCY7fl5tFBHxAQX4KSKhYGKy461tmuxYRDKeAvxdrptdTcuxfra3arJjEclsCvB3uSY5nPB5DScUkQynAH+X2tJ85k4Yx4otrV6XIiLynhTgw7h+TjVr9x+lsyfqdSkiIqelAB/G0jk1DDlN8iAimU0BPowFdaVUlURYsUUBLiKZSwE+jEDAWDq7mhe2txONDXldjojIsBTgp7F0Tg3dAzFe0yQPIpKhFOCnccWMSiKhAM9pNIqIZCgF+GkU5CXOylyxtVVnZYpIRlKAv4elc6o50NnHzjadlSkimUcB/h6um504K/NZtVFEJAMpwN/DhNICFtSV8uxmBbiIZB4F+BncNK+Gdfu7aD3e73UpIiLvoAA/gxvn1QLoKFxEMo4C/AxmVhcztbKIZxTgIpJhFOBnYGbcOLeGV3Z1cLx/0OtyRETeogAfgRvn1TIYd6zUNcJFJIMowEdgUf14qkoiPLNJbRQRyRwjDnAzC5rZOjN7Irk81cxWm9kOM/uFmeWlr0xvBQLGDXNrWLWtjf7BuNfliIgAZ3cE/kVgyynL3wG+65ybCRwFPpPKwjLNjXNr6InGeXlXh9eliIgAIwxwM5sE3AYsSy4bcB3wSHKVh4A701FgprhseiUlkRC/3XjY61JERICRH4H/I/CnwMmLY1cAXc65WHK5Gagb7gfN7H4zazKzpvb29lEV66W8UIClc6p5ZnMrg3FdI1xEvHfGADez9wFtzrk1pz49zKrDXrLPOfeAc67ROddYVVV1jmVmhlsXTKCrd5BXdx/xuhQRkREdgV8O3G5me4GHSbRO/hEYb2ah5DqTgJa0VJhBrjqviqK8IL/ZcMjrUkREzhzgzrlvOOcmOecagLuB551z9wArgQ8nV7sXeDxtVWaI/HCQpXNqeHpTKzG1UUTEY6MZB/414MtmtpNET/zB1JSU2W5dMIHOniiv7tZUayLirdCZV3mbc24VsCr5eDdwcepLymzXzEq0UZ7ccIgrZlZ6XY6I5DCdiXmW8sNBrptTwzObDquNIiKeUoCfg1vn13KkJ6oZ60XEUwrwc3DNrGoKwkGe0GgUEfGQAvwcFOQFuX5uDU9tOKSTekTEMwrwc3TH+RM52jvISzt1bRQR8YYC/BxddV4VpQVhlr+R9ecviUiGUoCfo7xQgFsX1PL0psP0RXWJWREZewrwUbj9/Dp6o3Ge26KJHkRk7CnAR+HiqeXUjIuwfL3aKCIy9hTgoxAMGO9fOJFV29o41qsJj0VkbCnAR+mOC+oYjDue2qgx4SIythTgozS/bhzTKov41RsHvS5FRHKMAnyUzIw7F9Xx6u5ODnT2el2OiOQQBXgKfGBRYja5/1yno3ARGTsK8BSoLy9kybRyHlvbjHPDziwnIpJyCvAU+fBF9ew90suafUe9LkVEcoQCPEVumV9LYV6QR9c2e12KiOQIBXiKFEVC3Dy/lifWH6J/UKfWi0j6KcBT6MMXTuLEQIynNx32uhQRyQEK8BRaMq2CuvEFPLJGbRQRST8FeAoFAsaHLqzjpZ0dNB/VmHARSS8FeIp9ZHE9AL98/YDHlYhItlOAp9ikskKumlnFL5uaNWu9iKSVAjwNPnpxPYeP9/PC9navSxGRLKYAT4Olc2qoLI7w89fURhGR9FGAp0E4GOCuxkk8v7WVw8f6vS5HRLKUAjxN7l5cz5CD/2jSUbiIpIcCPE2mVBRx+YwKHn79APEhXeBKRFJPAZ5G91wyhYNdfazc2uZ1KSKShRTgaXTD3Bpqx+Xz0Ct7vS5FRLKQAjyNwsEA91wymRd3dLC7vdvrckQkyyjA0+zuiycTDho/fXWf16WISJZRgKdZVUmE2xZM4JGmZnoGYl6XIyJZRAE+Bj55WQMnBmKaM1NEUuqMAW5m9Wa20sy2mNkmM/ti8vlyM3vWzHYk78vSX64/Laofz/y6cfzklb2aM1NEUmYkR+Ax4CvOuTnAEuBzZjYX+Dqwwjk3E1iRXJZhmBn3XTaV7a3dvLijw+tyRCRLnDHAnXOHnHNrk49PAFuAOuAO4KHkag8Bd6aryGzw/vMnUFUSYdlLe7wuRUSyxFn1wM2sAVgErAZqnHOHIBHyQHWqi8smkVCQ+y5r4Hfb29l2+ITX5YhIFhhxgJtZMfAo8CXn3PGz+Ln7zazJzJra23P78qr3XDKZgnCQZS/u9roUEckCIwpwMwuTCO+fOeceSz7damYTkq9PAIY9X9w594BzrtE511hVVZWKmn1rfGEedzVO4vE3Wmg7oasUisjojGQUigEPAlucc/9wykvLgXuTj+8FHk99ednn05dPZXBoiJ+8rBN7RGR0RnIEfjnwCeA6M3sjebsV+FvgBjPbAdyQXJYzaKgs4sa5Nfzb6n30RnVij4icu9CZVnDOvQTYaV5emtpycsP9V03n6U2t/Pvq/Xz2ymlelyMiPqUzMT1w0ZQyLp1Wwb+8uJuBWNzrckTEpxTgHvnctTNoPT7Ao2t0er2InBsFuEcun1HB+ZNK+eELu4jFh7wuR0R8SAHuETPjc9fOYH9nL0+8ecjrckTEhxTgHrp+Tg3n1RTz/ZU7GdK8mSJylhTgHgoEEkfhO9q6eWrjYa/LERGfUYB77H0LJzKzupjvPrdds9eLyFlRgHssGDC+dP157Gzr5tfrW7wuR0R8RAGeAW6ZX8vs2hK+t2KHRqSIyIgpwDNAIGB8+Ybz2NPRw2Oadk1ERkgBniFumFvDwkml/NOKHURjOgoXkTNTgGcIs8RRePPRPn7+2n6vyxERH1CAZ5Crz6vi0mkVfG/FDo73D3pdjohkOAV4BjEz/uzWOXT2RPnhql1elyMiGU4BnmEWTCrlzgsm8uBLe2jp6vO6HBHJYArwDPSVG2fhHPzDs9u9LkVEMpgCPAPVlxdy3+UNPLq2mU0tx7wuR0QylAI8Q33umhmUFebxreWbcU6n2IvIH1KAZ6jSwjBfvWkWr+3tZLlOsReRYSjAM9hHGutZOKmUv35yC90DmgBZRN5JAZ7BggHjW7fPo+3EAP/8/A6vyxGRDKMAz3CLJpfx4Ysm8aOX9rCzrdvrckQkgyjAfeBrN8+mIBzkm/+5QTP3iMhbFOA+UFUS4Zu3zWH1nk5+0XTA63JEJEMowH3iI431LJlWzrd/s4W24/1elyMiGUAB7hNmxt98cCEDsSH+Yvkmr8sRkQygAPeRqZVFfOn6mTy18TBPbTjkdTki4jEFuM/80ZXTmF83jm/+aiNtJ9RKEcllCnCfCQcDfPcjF9AzEOPrj27QafYiOUwB7kMza0r4+i2zeX5rGw+/rlEpIrlKAe5T917awBUzKvnfT2xmb0eP1+WIiAcU4D4VCBj/566FhALG53++jv7BuNclicgYU4D72ITSAv7+Ixew4eAx/vrJLV6XIyJjTAHuczfMreH+q6bx01f36bKzIjkm5HUBMnpfvWkWa/cd5RuPvsm8ieOYXlXsdUmSYWLxIboHYm/degbi9A/G6YvG6RuME40NMRhP3GJDjnjyduoYJwMCZphBKGCEggFCASMcDBAOBcgLBoiEA0RCASKhIPnhAAXhIAV5wbfu84IBzMyrv4asY6MZhmZmNwPfA4LAMufc377X+o2Nja6pqemc309O79CxPm77p5cYXxjmP//75ZQWhL0uSdJsaMjR3j1A89E+Wo/303q8n8PH++k4EeVIzwCdPVE6e6Ic6x3kRIZcTz4UMAryghTlhSiMJO6LIkGKIyGKIyGKIiGK80OUJJeL88MUR0KU5Ifevs8PMS4/TCSUOx8GZrbGOdf4B8+fa4CbWRDYDtwANAOvAx91zm0+3c8owNNr9e4jfPzB1SyZVsGP71tMKKgOmd/FhxwHj/axq72bXe3d7D3Sw74jvezv7KWlq4/B+Dv//4aDRmVxhMriCBXFeZQV5lFaEGZ8YZiS/DAlyXAsyAtSmBeiIBwkEk4cPeeFAoSTR9XBoBFMHm0bhsPhHLhkTfEhR2xoiFjcvXXkPhAbIhpL3PcPxukfTNz3Dcbpjcbpi8bojcaTt8RvAT3RGD0DMboH4vQMJB6fGIgRjQ2d8e8mHLRkyIcoiYTfCv6TIV8cCVOc/HAoOuUD4u3Hb78WzvD/K6cL8NG0UC4Gdjrndiff4GHgDuC0AS7pdcm0Cv7qzvl87dEN/NWTW/hft8/zuiQ5C129UTa1HGdTyzG2Hj7BtsMn2NnWzcApYTYuP0RDZREL6kq5dcEEJo4vYNL4AmpL86kZl09ZYTgrjkqjsWTLpz/G8f7Btx6fGBhM3sc40Z8M/P4YJ5LrHD7ez462s/sgAMgLBiiMBClMtnqKIokPt8K8xHJ+8nF+6O3lSChAfvjtx5FQgEg4+NaHYST09odiXihAdUkk5R8UownwOuDUs0iagUvevZKZ3Q/cDzB58uRRvJ2MxH9ZPJkdrd0se2kP06uK+MSlDV6XJMPoH4yz4eAx3tjfxRvNXaw/0EXz0b63Xq8uiTCrtoRPLJnCzJpiplcVM62qmPKiPA+rHjt5oQDlobxRb+9ALJ440j8Z+NGT3wGcvCVe6x2M05v8TaBv8O3fFDq6o/Sd8l1Bf/J2Lpflf+7LVzGjumRU2/Nuownw4T7m/2CznHMPAA9AooUyiveTEfrGrXPYe6SHP1++ifKiCLctnOB1STnveP8gr+/pZPWeTpr2drLh4LG32h914wu4oH4891wyhfl145g3sTRngjrdIqEgkVAwpX+fzjmi8SH6B4cYSLaKovHE/dttpMQXw9FkeykaG6J6XH7KajhpNAHeDNSfsjwJ0Di2DBAMGP/80Qv55I9W86VfrKMkP8RV51V5XVZO6R+M07T3KC/ubOflnUfY1HKMIZf4VX3hpFI+fcVUGqeUc0H9eKpKIl6XK2fBzN76YMDjwQKj+RIzROJLzKXAQRJfYn7MOXfai1XrS8yxdaxvkLsfeJW9HT3822cv4aIpZV6XlLWcc+zu6GHVtnZWbWvjtT2dDMSGCAeNRfVlLJlewaXTKlg0eTz54aDX5YrPpPxLTOdczMz+B/A0iWGEP3qv8JaxV1oQ5qFPL+YjP3yF+370Gv/66cVcNKXc67KyxmB8iNf2dPLcllae39rGviO9AEyvKuKeS6Zw5cxKLp5aTlFEp1tIeoxqHPjZ0hG4N1q6+rhn2Wpaj/fzo/sWs2Rahdcl+VbPQIwXtrfzzKbDrNjaxon+GHmhAJdNr2Dp7GqumVVNfXmh12VKlkn5OPBzoQD3Ttvxfj62bDXNR3v5l082cuVM9cRH6ljfICu2tPLUxsP8bns7A7EhygrDXD+nhhvm1nDFzEoK83SULemjABc6ugf4+LLV7Grv5tsfWMBdjfVn/qEcdbx/kOc2t/Lkm4f43Y52BuOO2nH53Dy/lpvn19I4pUwnSsmYSceJPOIzlcURfvHHl/K5n63lq4+8yd4jPXzlhlkEAv4/8SMVeqMxVmxp49frW1i1rZ1ofIiJpfnce2kDtyyYwKL68fq7koyiAM8xpQVhfvypxfz54xv5/spd7Ono4TsfWkhJfm5eOyUaG+J329t5fH0Lz21upW8wTnVJhI8vmcJtCxXaktkU4DkoHAzw7Q8sYFplMX/7261sbnmJ//uxC5lfV+p1aWMiPuRYvecIy99o4amNhznWN0hZYZgPXFjH+xdO5OKp5QQV2uIDCvAcZWb80VXTuGDyeD7/7+v44P97mT+7dTafvLQhK484nXOsbz7G8jdaeHJDC63HByjKC3LjvFpuP38iV8yszPgLGom8m77EFDp7ovzP/1jP81vbWNxQxt98cEHKr9ngBeccWw+f4Ik3W3jizUPsO9JLXjDA1bOquOOCiSydXUNBnk6qkcynUSjynpxzPLKmmb96cgt90Tj/7Zrp/PHV03w3PM45x7bWE/xmw2GefLOFXe09BAPGZdMreP/5E7lpXq2ulS6+owCXEenoHuAvf72Z5etbqCqJ8IWlM7l7cX1GtxeGhhzrm7t4ZnMrT286zO72HgIGF08t530LJ3LL/FoqinW9EfEvBbiclaa9nXznt1t5fe9RJpcX8pkrpvKhiyZRnCGnhfcMxPj9zg5WbmtjxZY22k4MEAoYl0wr55b5E7hpXq0uEiVZQwEuZ805x6pt7XxvxQ7eONBFSSTEXY31fPDCOuZNHDemEwfEhxwbDx7jpZ0d/H5nB017jxKND1EcCXHlzEpunFfDdbNqKC1Ue0SyjwJcRmXd/qP8+Pd7+c2GQ8SGHFMri3jfwglcObOK8+tLE5fWTKGjPVE2thxj3f4uXt/bybr9XXQn53WcXVvClTMruXZWNY0N5eSFMre9I5IKCnBJiaM9UX676TBPvNnCK7uOMOQgPxzgwsllzJ0wjvNqS5hZXczE8QWUF+W9Z++8ZyD21kS8Bzp72dXew+72brYcOsHBrsTsNGYwq6aExoYyFjeUc9n0SrVGJOcowCXlunqjvLank1d2H+H1vZ3saH3n/I0AZYVhCsJBQsEAoaAxGB+iNzmZbf/gO9fNCwaYWlnEjJpiFtSVsqCulPkTS9UWkZyna6FIyo0vzOPGebXcOK8WSPSp93f2srOtm9bj/XR0D9DRPUD/4BCx+BCDQ47Iyclj80KUFeZRWxqhpiSfSWWF1JUV6AxIkbOgAJeUCQaMqZVFTK0s8roUkZygb39ERHxKAS4i4lMKcBERn1KAi4j4lAJcRMSnFOAiIj6lABcR8SkFuIiIT43pqfRm1g7sO8cfrwQ6UliOX+TidufiNkNubre2eWSmOOeq3v3kmAb4aJhZ03DXAsh2ubjdubjNkJvbrW0eHbVQRER8SgEuIuJTfgrwB7wuwCO5uN25uM2Qm9utbR4F3/TARUTknfx0BC4iIqdQgIuI+JQvAtzMbjazbWa208y+7nU96WBm9Wa20sy2mNkmM/ti8vlyM3vWzHYk78u8rjXVzCxoZuvM7Ink8lQzW53c5l+YWZ7XNaaamY03s0fMbGtyn1+a7fvazP4k+W97o5n93Mzys3Ffm9mPzKzNzDae8tyw+9YS/imZbW+a2YVn814ZH+BmFgS+D9wCzAU+amZzva0qLWLAV5xzc4AlwOeS2/l1YIVzbiawIrmcbb4IbDll+TvAd5PbfBT4jCdVpdf3gN8652YD55PY/qzd12ZWB3wBaHTOzQeCwN1k577+V+Dmdz13un17CzAzebsf+MHZvFHGBzhwMbDTObfbORcFHgbu8LimlHPOHXLOrU0+PkHiP3QdiW19KLnaQ8Cd3lSYHmY2CbgNWJZcNuA64JHkKtm4zeOAq4AHAZxzUedcF1m+r0lM4VhgZiGgEDhEFu5r59zvgM53PX26fXsH8BOX8Cow3swmjPS9/BDgdcCBU5abk89lLTNrABYBq4Ea59whSIQ8UO1dZWnxj8CfAienqK8AupxzseRyNu7vaUA78ONk62iZmRWRxfvaOXcQ+DtgP4ngPgasIfv39Umn27ejyjc/BPhw05Rn7dhHMysGHgW+5Jw77nU96WRm7wPanHNrTn16mFWzbX+HgAuBHzjnFgE9ZFG7ZDjJnu8dwFRgIlBEon3wbtm2r89kVP/e/RDgzUD9KcuTgBaPakkrMwuTCO+fOeceSz7devJXquR9m1f1pcHlwO1mtpdEa+w6Ekfk45O/ZkN27u9moNk5tzq5/AiJQM/mfX09sMc51+6cGwQeAy4j+/f1Safbt6PKNz8E+OvAzOS31XkkvvhY7nFNKZfs/T4IbHHO/cMpLy0H7k0+vhd4fKxrSxfn3Decc5Occw0k9uvzzrl7gJXAh5OrZdU2AzjnDgMHzGxW8qmlwGayeF+TaJ0sMbPC5L/1k9uc1fv6FKfbt8uBTyZHoywBjp1stYyIcy7jb8CtwHZgF/BNr+tJ0zZeQeJXpzeBN5K3W0n0hFcAO5L35V7XmqbtvwZ4Ivl4GvAasBP4DyDidX1p2N4LgKbk/v4VUJbt+xr4FrAV2Aj8FIhk474Gfk6izz9I4gj7M6fbtyRaKN9PZtsGEqN0RvxeOpVeRMSn/NBCERGRYSjARUR8SgEuIuJTCnAREZ9SgIuI+JQCXETEpxTgIiI+9f8BpZmLGAEhJ14AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot([compute_inputs(t, state) for t, state in zip(t_eval, traj)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a26b2b190>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxddZ3/8dfn3put2dMkbbO0TXe6UhpKKTvoQBmg8hgFigooyoyKOqMzijo4/pyf89AZ/SECIxZUQGdAdgpUURYFhEJT2qYt3ULpkqZL2iZpmn35/v64txBCAm1zkpN7z/v5eNxH7jn3cL+fw0nf5+T7PYs55xARkcQX8rsAEREZGgp8EZGAUOCLiASEAl9EJCAU+CIiARHxu4D+5Ofnu/Hjx/tdhohIXFm1atUB51xBX58N28AfP348FRUVfpchIhJXzGxHf5+pS0dEJCAU+CIiAaHAFxEJCAW+iEhAKPBFRALCk8A3s1+Z2X4zW9/P52ZmPzOzKjOrNLNTvGhXRESOnVdH+PcAF33A54uAybHXDcDPPWpXRESOkSeB75x7ETj0AYssBu5zUSuAHDMb40XbvdU3t3Prs1tZV90wGF8vIhK3hqoPvxjY1WO6OjbvPczsBjOrMLOK2traE2ooFDJueXYLf968/8QqFRFJUEMV+NbHvPc9ecU5t9Q5V+6cKy8o6PPK4A+VlZrEhPx0KnfrCF9EpKehCvxqoLTHdAlQM1iNzSzOZr0CX0TkPYYq8JcB18TO1lkANDjn9gxWY7OKs9nT0MqBI22D1YSISNzx5OZpZnY/cC6Qb2bVwL8BSQDOuTuB5cDFQBXQDHzGi3b7M7M4G4B1uxs4b2rhYDYlIhI3PAl859ySD/ncAV/yoq1jMbM4C4D11Qp8EZGjEvJK20wN3IqIvE9CBj5o4FZEpLeEDfzZJRq4FRHpKWEDv+fArYiIJHDgzyh6d+BWREQSOPA1cCsi8l4JG/gAs0o0cCsiclRiB76uuBUReUdCB74GbkVE3pXQgT+jKAszdG98ERESPPAzU5Moy0/XEb6ICAke+ACzi7N1hC8iQgACf1ZJDnsPt7K/sdXvUkREfJXwgT+7JDpwq9MzRSToEj7wp4+JDtxWqltHRAIu4QM/PSXCpIIM9eOLSOAlfOBD9Ipb3WJBRIIuEIE/uzib2sY29h3WwK2IBFcgAn9WbOBW/fgiEmSBCPzpY7IJGayrrve7FBER3wQi8NOSw0wZlal+fBEJtEAEPkTvnLl+dwPOOb9LERHxhSeBb2YXmdlmM6sys5v6+Hysmb1gZqvNrNLMLvai3eMxqySbA0fa2dOggVsRCaYBB76ZhYE7gEXAdGCJmU3vtdi/Ag865+YCVwH/PdB2j9esYg3cikiweXGEPx+ocs5tc861Aw8Ai3st44Cs2PtsoMaDdo/LSWOyiISMSg3cikhAeRH4xcCuHtPVsXk9fQ/4lJlVA8uBL/f1RWZ2g5lVmFlFbW2tB6W9KzUpOnCrWyWLSFB5EfjWx7zeI6NLgHuccyXAxcBvzOx9bTvnljrnyp1z5QUFBR6U9l5zSrOprNbArYgEkxeBXw2U9pgu4f1dNtcDDwI4514FUoF8D9o+LrOKc2ho6WDnoeahblpExHdeBP5KYLKZlZlZMtFB2WW9ltkJXABgZicRDXxv+2yOwdFbJa/VwK2IBNCAA9851wncCDwDbCR6Ns4GM/u+mV0WW+zrwOfNbC1wP3Cd86FfZeroTFIiISp3aeBWRIIn4sWXOOeWEx2M7Tnvuz3evwmc4UVbA5EUDjG9KEtX3IpIIAXmStujZseuuO3q1sCtiARL8AK/JIfm9i7eqj3idykiIkMqcIE/pzQ2cKt+fBEJmMAF/oT8DNKTw7oAS0QCJ3CBHwoZM4uzdWqmiARO4AIfYE5pDhtrDtPe2e13KSIiQyaQgT+7JJv2rm627Gv0uxQRkSETzMAvzgFgre6cKSIBEsjAL81LI3dEEpW71I8vIsERyMA3M2aX5OgIX0QCJZCBD9GB2y37Gmlq6/S7FBGRIRHYwJ9bmkO3g/U6H19EAiKwgf/urZLVrSMiwRDYwB+ZkUJpXhprdIsFEQmIwAY+wJySHNbqTB0RCYhAB/7JpTnsrm9hf2Or36WIiAy6wAc+oPPxRSQQAh34M4qyCYdMA7ciEgiBDvy05DBTR2Vq4FZEAiHQgQ/RC7DW7qqnW488FJEEF/jAn1uaw+HWTrYfbPK7FBGRQeVJ4JvZRWa22cyqzOymfpa5wszeNLMNZva/XrTrhTmlunOmiATDgAPfzMLAHcAiYDqwxMym91pmMvAt4Azn3AzgHwfarlcmFWYwIjms8/FFJOF5cYQ/H6hyzm1zzrUDDwCLey3zeeAO51wdgHNuvwfteiIcMmaXZLN6Z53fpYiIDCovAr8Y2NVjujo2r6cpwBQz+6uZrTCzi/r6IjO7wcwqzKyitrbWg9KOzdyxuby55zCtHV1D1qaIyFDzIvCtj3m9T3mJAJOBc4ElwN1mlvO+/8i5pc65cudceUFBgQelHZu5pTl0dDk21KhbR0QSlxeBXw2U9pguAWr6WOYJ51yHc+5tYDPRHcCwcPLY6L5n9U4N3IpI4vIi8FcCk82szMySgauAZb2WeRw4D8DM8ol28WzzoG1PFGamUpKbpsAXkYQ24MB3znUCNwLPABuBB51zG8zs+2Z2WWyxZ4CDZvYm8ALwL865gwNt20tzx+Zq4FZEElrEiy9xzi0Hlvea990e7x3wtdhrWJpbmsOTa2vY29DK6OxUv8sREfFc4K+0PWpurB9/zS4d5YtIYlLgx0wvyiI5HFI/vogkLAV+TEokzIziLN5QP76IJCgFfg9zS3OprG6go6vb71JERDynwO9h7tgc2jq72bSn0e9SREQ8p8Dv4ZRxuQCs1sCtiCQgBX4PRdmpFGamaOBWRBKSAr8HM2Pu2BxW7dARvogkHgV+L/PG5bLzUDO1jW1+lyIi4ikFfi/zxuUB6ChfRBKOAr+XmcVZJEdCrNpxyO9SREQ8pcDvJSUSZnZxNhU6wheRBKPA78O8cbms392gJ2CJSEJR4Pdh3rhcOroc63brCVgikjgU+H04egGWBm5FJJEo8PuQn5FCWX46FdsV+CKSOBT4/ThlbC5v7Kwj+uwWEZH4p8DvR/n4XA41tfP2gSa/SxER8YQCvx/l6scXkQSjwO/HxIIMslIjCnwRSRgK/H6EQsa8cbm6AEtEEoYngW9mF5nZZjOrMrObPmC5j5uZM7NyL9odbOXj86jaf4RDTe1+lyIiMmADDnwzCwN3AIuA6cASM5vex3KZwFeA1wba5lA5rSx6I7WV23VfHRGJf14c4c8Hqpxz25xz7cADwOI+lvt34D+BVg/aHBKzSrJJjoR4/W0FvojEPy8CvxjY1WO6OjbvHWY2Fyh1zj3lQXtDJiUSZm5pjgJfRBKCF4Fvfcx752olMwsBtwBf/9AvMrvBzCrMrKK2ttaD0gbutLI8NtQ0cKSt0+9SREQGxIvArwZKe0yXADU9pjOBmcCfzWw7sABY1tfArXNuqXOu3DlXXlBQ4EFpA3dqWR7dTufji0j88yLwVwKTzazMzJKBq4BlRz90zjU45/Kdc+Odc+OBFcBlzrkKD9oedKeMzSUcMl5/+6DfpYiIDMiAA9851wncCDwDbAQedM5tMLPvm9llA/1+v6WnRJhZnM3Kt3WELyLxLeLFlzjnlgPLe837bj/LnutFm0Np/vhc7n1lB60dXaQmhf0uR0TkhOhK22Mwv2wk7V3drN1V73cpIiInTIF/DE4dH72Rmi7AEpF4psA/Bjkjkpk6KpPXdD6+iMQxBf4xml+Wx6oddXR0dftdiojICVHgH6PTJ46kub2Lymo92FxE4pMC/xgtmDASgBXbdD6+iMQnBf4xyktPZtroTF59S4EvIvFJgX8cTp84kpXbD9HW2eV3KSIix02BfxxOnzCSts5u1uzU+fgiEn8U+MfhtLKRmMGr6scXkTikwD8O2SOSmFmUzSvqxxeROKTAP06nTxzJmp31tHaoH19E4osC/zidPiF6Xx3dH19E4o0C/zidWpZHOGQ6PVNE4o4C/zhlpESYXZLNK28d8LsUEZHjosA/AQsnjqSyWs+5FZH4osA/AWdMyqez27FC3ToiEkcU+Cdg3rhc0pLCvLS11u9SRESOmQL/BKREwiyYkMdLW9WPLyLxQ4F/gs6aXMC2A01U1zX7XYqIyDFR4J+gsybnA/CyjvJFJE4o8E/QpMIMRmelqltHROKGJ4FvZheZ2WYzqzKzm/r4/Gtm9qaZVZrZc2Y2zot2/WRmnDU5n5erDtDV7fwuR0TkQw048M0sDNwBLAKmA0vMbHqvxVYD5c652cDDwH8OtN3h4KwpBTS0dLButx57KCLDnxdH+POBKufcNudcO/AAsLjnAs65F5xzR0c3VwAlHrTruzMn5WMGL+v0TBGJA14EfjGwq8d0dWxef64Hft/XB2Z2g5lVmFlFbe3wD9G89GRmFGXxovrxRSQOeBH41se8Pju1zexTQDnwX3197pxb6pwrd86VFxQUeFDa4DtrcgFv7KijsbXD71JERD6QF4FfDZT2mC4BanovZGYfAb4DXOaca/Og3WHhnCkFdHY7/lqlo3wRGd68CPyVwGQzKzOzZOAqYFnPBcxsLvALomG/34M2h41543LJTI3w/KaEWi0RSUADDnznXCdwI/AMsBF40Dm3wcy+b2aXxRb7LyADeMjM1pjZsn6+Lu4khUOcPaWAFzbX0q3TM0VkGIt48SXOueXA8l7zvtvj/Ue8aGe4umBaIU9X7mFDzWFmlWT7XY6ISJ90pa0HzplSgBnq1hGRYU2B74GRGSmcXJrD85v2+V2KiEi/FPgeOX9qIWurG6htTJgTkEQkwSjwPXLetEIA/rxZ3ToiMjwp8D0yoyiLUVkpvKDAF5FhSoHvETPjvKmFvLTlAB1d3X6XIyLyPgp8D50/rZDGtk5e23bI71JERN5Hge+hs6cUkJYU5pkNe/0uRUTkfRT4HkpNCnPu1AKe2bBXV92KyLCjwPfYhTNGs7+xjdW76v0uRUTkPRT4HjtvWiFJYVO3jogMOwp8j2WnJbFwYj7PbNiLc+rWEZHhQ4E/CC6cMZodB5vZtLfR71JERN6hwB8EH50+CjP4w3p164jI8KHAHwQFmSmUj8tVP76IDCsK/EFy4YzRbNrbyPYDTX6XIiICePQAFHm/RbPG8H+f3shTlTXceP5kv8uRYayr23GwqY3axjb2N7ZR39xOQ3MHh1s7aWrrpLWji5aOLjq6HF3dju7YyQCRkBEOhUgKG6lJYVKSQqQlhclIiZCVmkRmaoTsEUnkpCWTm55EXnoyKZGwz2srflLgD5LinDTKx+XyxJoavnTeJMzM75LEZwePtLF5XyNb9x1h6/5GdhxsZtehZnbXt9DR1fcZXSmREGnJYVIjYZIiRtiMkBmO6I6iq9vR0dVNa0cXrZ3dtHd+8H2cslIj5GemUJiZwqisVAozUxidncaY7FTGZKdSlJNGQUYKoZB+XxORAn8QLT65iJuf2MCmvY2cNCbL73JkCLV3dlNZXc/r2w9RuauByup6ahpa3/k8MzVCWX46M4qzuWjmGIpyUinISKEwK4W89BSy06JH6Enh4+t17ezqpqm9iyNtnRxu6aC+uYP65nYONbdz6Eg7B460ceBIO/sbW1m9s559h1tp67WTSAobY7LTKM5JoyQ3jZLcEZTkplGaN4LSvDRGZaZqhxCnFPiD6OJZY/jek2+ybG2NAj/BOefYsu8If968nxe31rJqRx2tHdEgHT9yBPPG5/GZ4mymjclkyqhMCjNTBuWvvkg4RHZaiOy0JIpz0o6p7vrmDvY0tLKnoYWahlZq6lvYXdfC7voWXtxay77D732oT3IkRElOdAcwNvYqzTu6QxhBVmqS5+sl3lDgD6KRGSmcOSmfZWtq+MaFU9Wtk2C6uh2rdtSxfN0e/rhh7ztH8FNHZXLVqWNZMCGPU8fnMTIjxedK+2dm5KYnk5uezPSivg9KWju6qKlvYVddC7sONbOrrpnqQy3sONTEml31NLR0vGf5nBFJlOZGdwQleWnRHULsr4Ti3DSNI/jIk8A3s4uAW4EwcLdz7oe9Pk8B7gPmAQeBK51z271oe7i7bE4RX39oLW/srGPeuDy/yxEPbNp7mIcrqlm2tob9jW0kR0KcM6WAr1wwmXOmFjAm+8OPrONJalKYCQUZTCjI6PPzhpaO6I7gUDM7e7ze3HOYP7659z3jE2ZQmJnyTjdRcU50J1Cck0ZR7JWRouPQwTLg/7NmFgbuAD4KVAMrzWyZc+7NHotdD9Q55yaZ2VXAj4ArB9p2PPibGaNIeSzEsjU1Cvw41tTWyWOrd/PAyp2s332YSMg4b1ohl8wewwUnjQp0SGWnJZFdnM3M4uz3fdbV7dh3uDX2l0HLO4PU1XXNrNpRx9OVe+jsdWfZzNQIRdlpjI4NJBdmpTI6K5VRWSkUZqZSmJXCyPRkIsc5viHeHOHPB6qcc9sAzOwBYDHQM/AXA9+LvX8YuN3MzAXgZjOZqUlccFIhT6/bw82XTNcvaZzZfqCJe17ZziOrqmls6+SkMVn826XTuWxO0bDuqhkuwiF758j9tD4+P7pD2F3fQk19CzX10bGEPQ2t7G1oZUPNYQ42tdE7Kcwgb0Qy+Rkp5Gcmk5ce3QmMjHVP5Y5IJndEUvS01BHJZKclkZ4cDny3qheBXwzs6jFdDe/btu8s45zrNLMGYCRwoOdCZnYDcAPA2LFjPShteLhsTjHL1+3l5aoDnDu10O9y5BhUVtdz51/e4vfr9xIJGRfPGsM1p4/jlLG5gQ8NL/XcIfSno6ubA0fa2NvQ+s61Cvsb26JnHMV+VtbVc+hIO41tnR/YVmZqJPpKSSIjJUJGaoT0lAjpyWFGJEdITwmTmhQmLSlMWnKYlEgoeo1DJERyJERyOERS7GckbERCodj1EEYkHDttNhT9aRYdIzGD3r8xDqI7MQcOR7eDbhe9xqK7O1prQab3BxReBH5fv/29j9yPZRmcc0uBpQDl5eUJc/R/3rQCckck8VBFtQJ/mHtjZx23/GkLL209QGZqhH84ZyKfWTiewqxUv0sLrKRwiDHZacc0NtLW2UV9cwd1ze3UNXXQ0NJOQ0sHDS0dHG7p5HBrB4dbOmhs7eRIWyf7G1s5UttJc3tX7NXJcHh20dyxOTz2xTM8/14vAr8aKO0xXQLU9LNMtZlFgGwgMA9+TYmE+djcYn67YgeHmtrJS0/2uyTpZf3uBn7yx828sLmWvPRkblo0jU+eNpZMnWIYV1IiYUZlhRl1gjto5xztXd20tnfT0tFFW2cXrR3dtHV20R67sK2ts5uOrm46Yxe9dTtHZ5ejs/vdK6G7uh3OHT2Sdz2+P9oddZSZEYr9BRAOWWzayM8YnIzwIvBXApPNrAzYDVwFXN1rmWXAtcCrwMeB54PQf9/TlaeW8uu/bufx1bv57JllfpcjMdV1zfzkj1t4bPVuckYk8c2LpnHN6eNID/AgbJCZGSmRMCmRMNkk3s5+wL/VsT75G4FniJ6W+Svn3AYz+z5Q4ZxbBvwS+I2ZVRE9sr9qoO3Gm2mjs5hTks2DFbv4zBnj1Q/ss6a2Tm5/oYpfvvw2Bnzh3Il84dyJumhIEponhzHOueXA8l7zvtvjfSvwCS/aimdXnFrKdx5bT2V1A3NKc/wuJ5CcczxVuYcfPL2RvYdbuXxuMf984dRjuipVJN7pHMEhdOmcIlKTQvyuYteHLyyee/tAE5+8+zW+fP9qRmYk88gXTueWK09W2EtgqKNyCGWlJnHxrDE8uaaGm/92OmnJusR8KHR0dbP0xW3c+txWUiIh/v1jM7l6/ljCugGYBIyO8IfYleWlNLZ18lRl7xOZZDCs393Apbe9zH89s5mPnFTIc187h08vGKewl0DSEf4Qm1+Wx5RRGdzzynY+Pq9Eg7eDpKOrmzteqOL256vITU/mrmvK+ej0UX6XJeIrHeEPMTPjuoVlbKg5zMrtdX6Xk5Cq9jdy+X//lZ8+u5VL5xTxp386W2EvggLfF5fPLSY7LYl7Xnnb71ISinOO36zYwSW3vUxNfSt3fmoet1x5MjkjdKGbCKhLxxdpyWGuml/K3S+9ze76Fp0l4oFDTe184+G1PLtxP2dPKeDHn5hNYaZuhyDSk47wfXLN6eNxznHfq9v9LiXuvbbtIBff+hIvbjnAdy+Zzj3XnaqwF+mDAt8nxTlpXDhjNA+8vouW9i6/y4lL3d2O257bypK7VpCaFOLRLy7ks2eW6XmrIv1Q4PvoM2eU0dDSwSNvVPtdStw51NTOdfes5Cd/2sKlc4p46itn9fkADhF5lwLfR6eOz+Xk0hzu/MtbdHR1+11O3Fi9s45LfvYSK946yH9cPoufXnlyoJ84JXKsFPg+MjO+csEkqutaeGz1br/LGfacc/zPazu44hevEgoZj3xhIVefNlbXMogcIwW+z86bWsjM4izueKGKTh3l96u1o4tvPlLJdx5bz8KJ+Tz15TOZVaIuHJHjocD3mZnx5fMns+NgM0/qdgt9qqlv4YpfvMqDFdV85fxJ/Oq6U3VuvcgJUOAPAx89aRTTRmdy+/NVdA2H56sNI6+/fYhLb3uZbbVNLP30PL72N1N1HxyRE6TAHwZCoehR/lu1TTy9bo/f5QwLR6+avfquFWSnJfH4l87gb2aM9rsskbimwB8mFs0czbTRmfzkj5tp6wz2efntnd18+7H13Pz4es6eUsDjN57BpMIMv8sSiXsK/GEiFDK+ffFJ7DjYzH2v7PC7HN/UNrZx9V0ruP/1nXzx3IncdU25Hjso4hEF/jBy9pQCzp1awM+e38qhpna/yxly63c3cNntL7O+poHblszlGxdNU3+9iIcU+MPMdy4+ieb2Lm59dovfpQypJ9bs5u9+/gohMx7+h4VcOqfI75JEEo4Cf5iZPCqTJfNL+e1rO6naf8TvcgZdV7fjh7/fxFcfWMOckhyeuPEM3SJBZJAo8Iehf/rIFEYkhfnesg04l7inaTY0d/DZe1Zy51/e4pOnjeW3nzuN/IwUv8sSSVgDCnwzyzOzP5nZ1tjP3D6WOdnMXjWzDWZWaWZXDqTNIBiZkcI3F03j5aoD/G7lLr/LGRRb9jWy+I6XeeWtA/zg8pn84PJZJEd0/CEymAb6L+wm4Dnn3GTgudh0b83ANc65GcBFwE/NLGeA7Sa8q+ePZcGEPH7w9Eb2NLT4XY6n/rB+D5ff8Vea2ru4//ML+ORp4/wuSSQQBhr4i4F7Y+/vBT7WewHn3Bbn3NbY+xpgP1AwwHYTXihk/OjvZtPZ7fj2o+sSomuns6ubH/5+E//w2zeYPCqTJ288k/LxeX6XJRIYAw38Uc65PQCxn4UftLCZzQeSgbf6+fwGM6sws4ra2toBlhb/xo1M518unMoLm2t59I34vpvmwSNtXPvr17nzL29x9Wlj+d3fL2B0tp5KJTKUPvQm4mb2LNDXNe3fOZ6GzGwM8BvgWudcn7eFdM4tBZYClJeXx/8hrQeuXTie36/fw81PrGd2STaTR2X6XdJxW7XjEDf+72oONrXznx+fzRXlpX6XJBJIH3qE75z7iHNuZh+vJ4B9sSA/Guj7+/oOM8sCngb+1Tm3wssVSHThkHHbklMYkRzm73+zisbWDr9LOmbOOe56cRtX/mIFSeEQj35hocJexEcD7dJZBlwbe38t8ETvBcwsGXgMuM8599AA2wuk0dmp3H71Kew41Mw/P7Q2LvrzDzW18/n7VvGD5Ru54KRCnvzymTq/XsRnAw38HwIfNbOtwEdj05hZuZndHVvmCuBs4DozWxN7nTzAdgNnwYSRfGvRNJ7ZsI87Xqjyu5wP9ErVARbd+iIvbqnl5kumc+en5pGdpvvhiPhtQA8Cdc4dBC7oY34F8LnY+98Cvx1IOxJ1/ZllrNvdwI//uIXsEcl8esHwOp2xtaOLW57dwtIXtzEhP51fXnuqjupFhhE9+TmOmBk//sQcmto6ufnx9aRGQnximPSJr91Vzz8/tJat+4+wZP5Ybr7kJEYk69dLZDjRv8g4kxQOcfvVp/D5+yr45iOVJEdCLD652Ld6Wtq7uO35rfzixW0UZKRw72fnc84UXWYhMhwp8ONQalKYpZ8u57pfv85XH1jD9gPNfPn8SYSG+FbCz765j39btoHd9S18Yl4J/3rJdPXViwxjCvw4lZYc5t7Pzufbj67jlme3sGnvYX78iTmkpwz+Jt2yr5Ef/X4Tz23az5RRGfzuhgWcNmHkoLcrIgOjwI9jqUlhfnLFHKYXZfEfyzeydf8R/uPyWcwvG5zbFeyub+GWP23h0TeqSU+O8K1F0/jsmWUkhXXTM5F4oMCPc2bG586awNTRmdz0yDqu+MWr/N0pJXzr4mme3Wp4XXUDv3x5G09V7iEUMq4/s4wvnjuJ3PRkT75fRIaGDdeLeMrLy11FRYXfZcSV5vZObn++irte2kZSOMTlc4v59OnjmDY667i/q7axjT9s2MsTq3dTsaOOjJQIV5SX8rmzyijKSRuE6kXEC2a2yjlX3udnCvzEU7X/CHf+5S2eXFtDW2c3c8fmsHDiSE4dn8fJpTlkpyVh9u4Ab1e348CRNiqrG1izq46V2+uo2H6IbgcTC9JZMn8sV5xaqoeJi8QBBX5A1TW18/Cqap6qrGF9zWG6uqPbOilsZKclk54S5nBLB/UtHRz9NQiHjGmjM7lgWiF/O7uIKaMy3rNzEJHhTYEvNLd3smZnPetrGqhr7qC+uYOmtk6y0iLkpacwMj2Z6UVZzCzKJi057He5InKCPijwNWgbECOSIyyclM/CSfl+lyIiPtH5dCIiAaHAFxEJCAW+iEhAKPBFRAJCgS8iEhAKfBGRgFDgi4gEhAJfRCQghu2VtmZWC+wYwFfkAwc8KideBHGdIZjrHcR1hmCu9/Gu8zjnXJ+PnRu2gT9QZlbR3+XFiSqI6wzBXO8grjMEc729XGd16YiIBIQCX0QkIBI58Jf6XYAPgrjOEMz1DuI6QzDX27N1Ttg+fBEReZxuUpoAAAOUSURBVK9EPsIXEZEeFPgiIgGRcIFvZheZ2WYzqzKzm/yuZ7CYWamZvWBmG81sg5l9NTY/z8z+ZGZbYz9z/a7Va2YWNrPVZvZUbLrMzF6LrfPvzCzZ7xq9ZmY5ZvawmW2KbfPTE31bm9k/xX6315vZ/WaWmojb2sx+ZWb7zWx9j3l9bluL+lks3yrN7JTjaSuhAt/MwsAdwCJgOrDEzKb7W9Wg6QS+7pw7CVgAfCm2rjcBzznnJgPPxaYTzVeBjT2mfwTcElvnOuB6X6oaXLcCf3DOTQPmEF3/hN3WZlYMfAUod87NBMLAVSTmtr4HuKjXvP627SJgcux1A/Dz42kooQIfmA9UOee2OefagQeAxT7XNCicc3ucc2/E3jcSDYBiout7b2yxe4GP+VPh4DCzEuBvgbtj0wacDzwcWyQR1zkLOBv4JYBzrt05V0+Cb2uij2BNM7MIMALYQwJua+fci8ChXrP727aLgftc1Aogx8zGHGtbiRb4xcCuHtPVsXkJzczGA3OB14BRzrk9EN0pAIX+VTYofgp8A+iOTY8E6p1znbHpRNzmE4Ba4Nexrqy7zSydBN7WzrndwI+BnUSDvgFYReJv66P627YDyrhEC3zrY15Cn3dqZhnAI8A/OucO+13PYDKzS4D9zrlVPWf3sWiibfMIcArwc+fcXKCJBOq+6Uusz3oxUAYUAelEuzN6S7Rt/WEG9PueaIFfDZT2mC4BanyqZdCZWRLRsP8f59yjsdn7jv6JF/u536/6BsEZwGVmtp1od935RI/4c2J/9kNibvNqoNo591ps+mGiO4BE3tYfAd52ztU65zqAR4GFJP62Pqq/bTugjEu0wF8JTI6N5CcTHeRZ5nNNgyLWd/1LYKNz7v/1+GgZcG3s/bXAE0Nd22Bxzn3LOVfinBtPdNs+75z7JPAC8PHYYgm1zgDOub3ALjObGpt1AfAmCbytiXblLDCzEbHf9aPrnNDbuof+tu0y4JrY2ToLgIajXT/HxDmXUC/gYmAL8BbwHb/rGcT1PJPon3KVwJrY62KifdrPAVtjP/P8rnWQ1v9c4KnY+wnA60AV8BCQ4nd9g7C+JwMVse39OJCb6Nsa+D/AJmA98BsgJRG3NXA/0XGKDqJH8Nf3t22JduncEcu3dUTPYjrmtnRrBRGRgEi0Lh0REemHAl9EJCAU+CIiAaHAFxEJCAW+iEhAKPBFRAJCgS8iEhD/H6d/f0ue/uVOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot([compute_accelerations(t, state)[1] for t, state in zip(t_eval, traj)])"
   ]
  },
  {
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}