Skip to content

Instantly share code, notes, and snippets.

@stephentu

stephentu/cartpole_task_space.ipynb

Created April 24, 2020 02:13
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save stephentu/71fc82928bf2183e2bc7e341f173f218 to your computer and use it in GitHub Desktop.
Save stephentu/71fc82928bf2183e2bc7e341f173f218 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
cartpole_task_space.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 194,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from jax.config import config \n",
"config.update(\"jax_debug_nans\", True)\n",
"import jax.numpy as jnp\n",
"import matplotlib.pylab as plt\n",
"import jax\n",
"import scipy"
]
},
{
"cell_type": "code",
"execution_count": 195,
"metadata": {},
"outputs": [],
"source": [
"g = 9.81\n",
"cartpole_params = (5.0, 1.0, 1.0, 0.0, 0.0, 0.5)\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, kp = 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) + kp * theta_d\n",
" ])\n",
" \n",
" return M, v\n",
"\n",
"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 task_space_fn(q):\n",
" x, theta = q\n",
" gamma_1 = 0.1\n",
" gamma_2 = 0.5\n",
" #return gamma_1 * x**2 + gamma_2 * angdiff(theta, jnp.pi)**2 \n",
" #return x\n",
" return angdiff(theta, jnp.pi)**2\n",
" #return theta\n",
" #return x \n",
" \n",
"grad_task_space_fn = jax.grad(task_space_fn)"
]
},
{
"cell_type": "code",
"execution_count": 180,
"metadata": {},
"outputs": [],
"source": [
"def check_feedback_linearization(q, q_d, cartpole_params, x_dd_des):\n",
" M, v = cartpole_dynamics(q, q_d, cartpole_params)\n",
" f_x = v[0] - M[0, 1] * v[1]/M[1, 1] + (M[0, 0] - M[0, 1]**2 / M[1, 1]) * x_dd_des\n",
" q_dd = jax.scipy.linalg.solve(M, jnp.array([f_x, 0.0]) - v, sym_pos=True)\n",
" assert np.allclose(q_dd[0], x_dd_des)\n",
" \n",
"check_feedback_linearization(np.array([-0.1, 0.2]), np.array([-0.3, 0.4]), \n",
" cartpole_params, -8.383)"
]
},
{
"cell_type": "code",
"execution_count": 204,
"metadata": {},
"outputs": [],
"source": [
"def compute_feedback(q, q_d, cartpole_params, pd_params):\n",
" k0, k1 = pd_params\n",
" #print(\"q\", q, \"q_d\", q_d)\n",
" M, v = cartpole_dynamics(q, q_d, cartpole_params)\n",
"\n",
" H = grad_task_space_fn(q)\n",
" hbar = H[0] - H[1] * M[0, 1]/M[1, 1]\n",
" if np.abs(hbar) < 1e-6:\n",
" raise Exception(\"rank condition\")\n",
" \n",
" Hdotqdot = jnp.dot(jax.jvp(grad_task_space_fn, (q,), (q_d,))[1], q_d)\n",
" \n",
" #y_dd_des = - k1 * jnp.dot(H, q_d) - k0 * angdiff(task_space_fn(q), np.pi)\n",
" y_dd_des = - k1 * jnp.dot(H, q_d) - k0 * task_space_fn(q)\n",
"\n",
" x_dd_des = 1/hbar*( -Hdotqdot + H[1]*v[1]/M[1, 1] + y_dd_des)\n",
"\n",
" f_x = v[0] - M[0, 1] * v[1]/M[1, 1] + (M[0, 0] - M[0, 1]**2 / M[1, 1]) * x_dd_des\n",
" print(\"f_x\", f_x)\n",
" \n",
" return f_x\n",
"\n",
"def compute_explicit_feedback(q, q_d, cartpole_params, pd_params):\n",
" k0, k1 = pd_params\n",
" M, v = cartpole_dynamics(q, q_d, cartpole_params)\n",
" theta_dd = - k1 * q_d[1] - k0 * angdiff(q[1], np.pi)\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": 197,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"f_x 14.037316\n",
"14.037316\n",
"21.628494\n"
]
}
],
"source": [
"state_0 = jnp.array([0.0, -0.1, 0.1, 0.2])\n",
"print(compute_feedback(state_0[:2], state_0[2:], cartpole_params, pd_params))\n",
"print(compute_explicit_feedback(state_0[:2], state_0[2:], cartpole_params, pd_params))"
]
},
{
"cell_type": "code",
"execution_count": 198,
"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": 205,
"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, 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))"
]
},
{
"cell_type": "code",
"execution_count": 206,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"f_x -13.563093\n"
]
},
{
"data": {
"text/plain": [
"DeviceArray([ 0. , 0. , -2.512715 , 1.5207963], dtype=float32)"
]
},
"execution_count": 206,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"state_0 = jnp.array([3, 0.1, 0, 0])\n",
"closed_loop(0.0, state_0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 207,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"f_x -13.563093\n",
"f_x -13.56308\n",
"f_x -13.562799\n",
"f_x -13.562653\n",
"f_x -13.561922\n",
"f_x -13.561795\n",
"f_x -13.561633\n",
"f_x -13.561633\n",
"f_x -13.558741\n",
"f_x -13.557388\n",
"f_x -13.55077\n",
"f_x -13.54965\n",
"f_x -13.548281\n",
"f_x -13.548283\n",
"f_x -13.524956\n",
"f_x -13.52234\n",
"f_x -13.520065\n",
"f_x -13.523757\n",
"f_x -13.532532\n",
"f_x -13.534466\n",
"f_x -13.726027\n",
"f_x -14.551561\n",
"f_x -18.700415\n",
"f_x -18.783491\n",
"f_x -20.521353\n",
"f_x -22.113558\n",
"f_x -27.404533\n",
"f_x -31.638523\n",
"f_x -56.364258\n",
"f_x -58.58406\n",
"f_x -68.77334\n",
"f_x -76.734085\n",
"f_x -106.34645\n",
"f_x -132.06326\n",
"f_x 9536.273\n",
"f_x 833.31824\n",
"f_x 323.60867\n",
"f_x 284.4687\n",
"f_x -87.36459\n",
"f_x -93.85395\n",
"f_x -140.65024\n",
"f_x -152.9139\n",
"f_x -172.3255\n",
"f_x -173.32945\n",
"f_x -254.6078\n",
"f_x -336.50647\n",
"f_x 587.0704\n",
"f_x 388.47525\n",
"f_x 270.9444\n",
"f_x 271.26132\n",
"f_x -204.63098\n",
"f_x -225.28287\n",
"f_x -442.34717\n",
"f_x -533.6478\n",
"f_x -718.88385\n",
"f_x -718.2861\n",
"f_x -1900.407\n",
"f_x -12715.804\n",
"f_x 458.7968\n",
"f_x 384.0638\n",
"f_x 320.15262\n",
"f_x 321.80304\n",
"f_x -1132.7705\n",
"f_x -1604.6472\n",
"f_x 1479.1552\n",
"f_x 1096.1578\n",
"f_x 829.557\n",
"f_x 831.95374\n",
"f_x -917.5696\n",
"f_x -1067.2863\n",
"f_x -5769.298\n",
"f_x -27544.826\n",
"f_x 7429.9185\n",
"f_x 7471.8467\n",
"f_x -771.6475\n",
"f_x -801.5243\n",
"f_x -993.74536\n",
"f_x -1038.073\n",
"f_x -1099.3712\n",
"f_x -1099.3407\n",
"f_x -1838.836\n",
"f_x -2791.2942\n",
"f_x 1753.0596\n",
"f_x 1354.977\n",
"f_x 1056.3066\n",
"f_x 1057.9884\n",
"f_x -1482.4808\n",
"f_x -1798.8511\n",
"f_x 26925.408\n",
"f_x 6975.4463\n",
"f_x 3623.6333\n",
"f_x 3629.0232\n",
"f_x -1244.9625\n",
"f_x -1333.6576\n",
"f_x -2071.0945\n",
"f_x -2297.4783\n",
"f_x -2661.0461\n",
"f_x -2660.7686\n",
"f_x -5685.493\n",
"f_x -13319.967\n",
"f_x 2329.6152\n",
"f_x 1924.7859\n",
"f_x 1581.5221\n",
"f_x 1582.248\n",
"f_x -4179.991\n",
"f_x -5863.4707\n",
"f_x 5782.369\n",
"f_x 4268.8643\n",
"f_x 3216.6284\n",
"f_x 3217.583\n",
"f_x -3456.6067\n",
"f_x -4067.0203\n",
"f_x -34746.645\n",
"f_x 100992.19\n",
"f_x 17164.68\n",
"f_x 17171.492\n",
"f_x -2942.8594\n",
"f_x -3107.8508\n",
"f_x -4318.317\n",
"f_x -4639.8164\n",
"f_x -5116.064\n",
"f_x -5116.018\n",
"f_x -10070.526\n",
"f_x -19600.58\n",
"f_x 5251.703\n",
"f_x 4283.5435\n",
"f_x 3481.206\n",
"f_x 3481.6362\n",
"f_x -7840.793\n",
"f_x -10698.05\n",
"f_x 13013.398\n",
"f_x 9330.519\n",
"f_x 6891.8784\n",
"f_x 6892.4688\n",
"f_x -6572.7896\n",
"f_x -7666.335\n",
"f_x -45592.82\n",
"f_x -381203.34\n",
"f_x 46471.445\n",
"f_x 46479.12\n",
"f_x -5534.732\n",
"f_x -5771.0747\n",
"f_x -7337.4956\n",
"f_x -7709.695\n",
"f_x -8231.637\n",
"f_x -8231.515\n",
"f_x -15291.708\n",
"f_x -26820.473\n",
"f_x 9682.386\n",
"f_x 7793.955\n",
"f_x 6265.9536\n",
"f_x 6266.232\n",
"f_x -12296.951\n",
"f_x -16338.892\n",
"f_x 25394.418\n",
"f_x 17458.832\n",
"f_x 12553.882\n",
"f_x 12554.441\n",
"f_x -10394.012\n",
"f_x -11967.428\n",
"f_x -49243.848\n",
"f_x -110422.91\n",
"f_x 199601.14\n",
"f_x 199601.14\n",
"f_x -8919.92\n",
"f_x -9309.162\n",
"f_x -11908.324\n",
"f_x -12530.351\n",
"f_x -13405.779\n",
"f_x -13405.78\n",
"f_x -25442.336\n",
"f_x -46222.74\n",
"f_x 14986.485\n",
"f_x 12128.595\n",
"f_x 9793.813\n",
"f_x 9793.983\n",
"f_x -20412.98\n",
"f_x -27646.414\n",
"f_x 35821.41\n",
"f_x 25434.625\n",
"f_x 18667.178\n",
"f_x 18667.797\n",
"f_x -17303.416\n",
"f_x -20249.688\n",
"f_x -136186.0\n",
"f_x 7458742.0\n",
"f_x 105483.21\n",
"f_x 105483.21\n",
"f_x -14659.359\n",
"f_x -15378.769\n",
"f_x -20379.67\n",
"f_x -21630.059\n",
"f_x -23427.174\n",
"f_x -23427.174\n",
"f_x -46586.71\n",
"f_x -92186.43\n",
"f_x 23668.152\n",
"f_x 19344.648\n",
"f_x 15747.721\n",
"f_x 15747.721\n",
"f_x -36849.383\n",
"f_x -51656.06\n",
"f_x 51184.508\n",
"f_x 37801.082\n",
"f_x 28487.348\n",
"f_x 28488.791\n",
"f_x -31256.676\n",
"f_x -37530.992\n",
"f_x 10339798.0\n",
"f_x 206114.44\n",
"f_x 92622.72\n",
"f_x 92622.72\n",
"f_x -26459.332\n",
"f_x -28290.771\n",
"f_x -43259.258\n",
"f_x -47750.312\n",
"f_x -54873.55\n",
"f_x -54873.55\n",
"f_x -132619.48\n",
"f_x -455675.12\n",
"f_x 40766.71\n",
"f_x 34150.44\n",
"f_x 28390.973\n",
"f_x 28390.973\n",
"f_x -94131.6\n",
"f_x -146553.86\n",
"f_x 82075.33\n",
"f_x 64260.05\n",
"f_x 50535.793\n",
"f_x 50535.793\n",
"f_x -79466.31\n",
"f_x -102415.51\n",
"f_x 230543.11\n",
"f_x 146075.89\n",
"f_x 100200.17\n",
"f_x 100200.17\n",
"f_x -68088.23\n",
"f_x -77419.54\n",
"f_x -245619.19\n",
"f_x -400200.8\n",
"f_x -1876090.5\n",
"f_x -1876090.5\n",
"f_x -61194.273\n",
"f_x -64934.023\n",
"f_x -93536.38\n",
"f_x -101472.59\n",
"f_x -113522.75\n",
"f_x -113522.75\n",
"f_x -258139.23\n",
"f_x -710996.6\n",
"f_x 91492.18\n",
"f_x 76191.08\n",
"f_x 63028.15\n",
"f_x 63028.137\n",
"f_x -192336.84\n",
"f_x -294521.56\n",
"f_x 177706.61\n",
"f_x 138278.88\n",
"f_x 108266.28\n",
"f_x 108266.28\n",
"f_x -164046.84\n",
"f_x -210959.98\n",
"f_x 490370.38\n",
"f_x 308132.47\n",
"f_x 210435.14\n",
"f_x 210435.14\n",
"f_x -141144.42\n",
"f_x -160670.78\n",
"f_x -521639.88\n",
"f_x -869330.5\n",
"f_x -5174654.0\n",
"f_x -5174654.0\n",
"f_x -126019.03\n",
"f_x -133374.73\n",
"f_x -188089.25\n",
"f_x -202888.83\n",
"f_x -225065.67\n",
"f_x -225065.67\n",
"f_x -506140.47\n",
"f_x -1349241.0\n",
"f_x 184168.81\n",
"f_x 153243.39\n",
"f_x 126629.695\n",
"f_x 126629.695\n",
"f_x -382240.22\n",
"f_x -587620.1\n",
"f_x 348862.75\n",
"f_x 271783.22\n",
"f_x 212986.48\n",
"f_x 212986.48\n",
"f_x -328036.53\n",
"f_x -425334.06\n",
"f_x 880594.9\n",
"f_x 569913.5\n",
"f_x 395218.25\n",
"f_x 395218.25\n",
"f_x -283958.34\n",
"f_x -326512.2\n",
"f_x -1314534.5\n",
"f_x -2844356.5\n",
"f_x 6220662.5\n",
"f_x 6220662.5\n",
"f_x -245726.55\n",
"f_x -257665.95\n",
"f_x -339323.44\n",
"f_x -359735.66\n",
"f_x -388850.53\n",
"f_x -388850.53\n",
"f_x -942183.6\n",
"f_x -3275236.2\n",
"f_x 287636.44\n",
"f_x 241123.92\n",
"f_x 200455.22\n",
"f_x 200455.22\n",
"f_x -678391.0\n",
"f_x -1079720.5\n",
"f_x 549845.6\n",
"f_x 433631.0\n",
"f_x 342904.34\n",
"f_x 342904.34\n",
"f_x -580939.94\n",
"f_x -772556.44\n",
"f_x 1197660.4\n",
"f_x 825022.2\n",
"f_x 593374.9\n",
"f_x 593374.9\n",
"f_x -506140.47\n",
"f_x -596347.44\n",
"f_x -5477166.5\n",
"f_x 12123803.0\n",
"f_x 2408831.8\n",
"f_x 2408831.8\n",
"f_x -432867.16\n",
"f_x -458649.4\n",
"f_x -655451.75\n",
"f_x -709202.5\n",
"f_x -791037.75\n",
"f_x -791037.75\n",
"f_x -1764219.2\n",
"f_x -4621594.0\n",
"f_x 654136.6\n",
"f_x 543466.7\n",
"f_x 448865.8\n",
"f_x 448865.8\n",
"f_x -1355749.0\n",
"f_x -2100408.5\n",
"f_x 1190054.0\n",
"f_x 930134.7\n",
"f_x 731607.56\n",
"f_x 731607.56\n",
"f_x -1172017.5\n",
"f_x -1545895.0\n",
"f_x 2610157.5\n",
"f_x 1765704.2\n",
"f_x 1256557.0\n",
"f_x 1256557.0\n",
"f_x -1024615.8\n",
"f_x -1204669.0\n",
"f_x -9271141.0\n",
"f_x 45456052.0\n",
"f_x 5491441.0\n",
"f_x 5491441.0\n",
"f_x -874741.44\n",
"f_x -922092.1\n",
"f_x -1278653.5\n",
"f_x -1372297.5\n",
"f_x -1512617.2\n",
"f_x -1512617.2\n",
"f_x -3182522.8\n",
"f_x -7059105.0\n",
"f_x 1379940.5\n",
"f_x 1137095.2\n",
"f_x 933222.8\n",
"f_x 933222.8\n",
"f_x -2547845.0\n",
"f_x -3886654.5\n",
"f_x 2419198.5\n",
"f_x 1877768.5\n",
"f_x 1466328.9\n",
"f_x 1466328.9\n",
"f_x -2225119.5\n",
"f_x -2918184.8\n",
"f_x 5335089.0\n",
"f_x 3549615.2\n",
"f_x 2505457.0\n",
"f_x 2505457.0\n",
"f_x -1954360.5\n",
"f_x -2288530.5\n",
"f_x -15345048.0\n"
]
},
{
"ename": "Exception",
"evalue": "rank condition",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mException\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-207-b481877f231c>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mt_end\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m20\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mt_eval\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_end\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnum\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mintegrate\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve_ivp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mclosed_loop\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_end\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstate_0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_eval\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mt_eval\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/scipy/integrate/_ivp/ivp.py\u001b[0m in \u001b[0;36msolve_ivp\u001b[0;34m(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, **options)\u001b[0m\n\u001b[1;32m 500\u001b[0m \u001b[0mstatus\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 501\u001b[0m \u001b[0;32mwhile\u001b[0m \u001b[0mstatus\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 502\u001b[0;31m \u001b[0mmessage\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msolver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstep\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 503\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 504\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0msolver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstatus\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'finished'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/scipy/integrate/_ivp/base.py\u001b[0m in \u001b[0;36mstep\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 180\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 181\u001b[0m \u001b[0mt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 182\u001b[0;31m \u001b[0msuccess\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmessage\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_step_impl\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 183\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 184\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0msuccess\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/scipy/integrate/_ivp/rk.py\u001b[0m in \u001b[0;36m_step_impl\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 141\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 142\u001b[0m y_new, f_new, error = rk_step(self.fun, t, y, self.f, h, self.A,\n\u001b[0;32m--> 143\u001b[0;31m self.B, self.C, self.E, self.K)\n\u001b[0m\u001b[1;32m 144\u001b[0m \u001b[0mscale\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0matol\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmaximum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my_new\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mrtol\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 145\u001b[0m \u001b[0merror_norm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnorm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merror\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mscale\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/scipy/integrate/_ivp/rk.py\u001b[0m in \u001b[0;36mrk_step\u001b[0;34m(fun, t, y, f, h, A, B, C, E, K)\u001b[0m\n\u001b[1;32m 68\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mzip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mC\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 69\u001b[0m \u001b[0mdy\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mK\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0ms\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mh\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 70\u001b[0;31m \u001b[0mK\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0ms\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mc\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mh\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mdy\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 71\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 72\u001b[0m \u001b[0my_new\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0my\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mh\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mK\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mB\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/scipy/integrate/_ivp/base.py\u001b[0m in \u001b[0;36mfun\u001b[0;34m(t, y)\u001b[0m\n\u001b[1;32m 137\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mfun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 138\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnfev\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 139\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfun_single\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 140\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 141\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfun\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfun\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/scipy/integrate/_ivp/base.py\u001b[0m in \u001b[0;36mfun_wrapped\u001b[0;34m(t, y)\u001b[0m\n\u001b[1;32m 19\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 20\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mfun_wrapped\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 21\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 22\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 23\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mfun_wrapped\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m<ipython-input-205-4172c020afcf>\u001b[0m in \u001b[0;36mclosed_loop\u001b[0;34m(t, state)\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mclosed_loop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mq\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq_d\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mstate\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstate\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mf_x\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcompute_feedback\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mq\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq_d\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcartpole_params\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpd_params\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 5\u001b[0m \u001b[0mM\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcartpole_dynamics\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mq\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq_d\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcartpole_params\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0mq_dd\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mjax\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mM\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mjnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mf_x\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msym_pos\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m<ipython-input-204-ff03e39c954d>\u001b[0m in \u001b[0;36mcompute_feedback\u001b[0;34m(q, q_d, cartpole_params, pd_params)\u001b[0m\n\u001b[1;32m 7\u001b[0m \u001b[0mhbar\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mH\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mH\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mM\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mM\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhbar\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0;36m1e-6\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"rank condition\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 11\u001b[0m \u001b[0mHdotqdot\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mjnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mjax\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mjvp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgrad_task_space_fn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mq\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mq_d\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq_d\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mException\u001b[0m: rank condition"
]
}
],
"source": [
"t_end = 20\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": 187,
"metadata": {},
"outputs": [],
"source": [
"traj = res.y.T"
]
},
{
"cell_type": "code",
"execution_count": 188,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
" message: 'The solver successfully reached the end of the integration interval.'\n",
" nfev: 392\n",
" njev: 0\n",
" nlu: 0\n",
" sol: None\n",
" status: 0\n",
" success: True\n",
" t: array([ 0. , 0.2020202 , 0.4040404 , 0.60606061, 0.80808081,\n",
" 1.01010101, 1.21212121, 1.41414141, 1.61616162, 1.81818182,\n",
" 2.02020202, 2.22222222, 2.42424242, 2.62626263, 2.82828283,\n",
" 3.03030303, 3.23232323, 3.43434343, 3.63636364, 3.83838384,\n",
" 4.04040404, 4.24242424, 4.44444444, 4.64646465, 4.84848485,\n",
" 5.05050505, 5.25252525, 5.45454545, 5.65656566, 5.85858586,\n",
" 6.06060606, 6.26262626, 6.46464646, 6.66666667, 6.86868687,\n",
" 7.07070707, 7.27272727, 7.47474747, 7.67676768, 7.87878788,\n",
" 8.08080808, 8.28282828, 8.48484848, 8.68686869, 8.88888889,\n",
" 9.09090909, 9.29292929, 9.49494949, 9.6969697 , 9.8989899 ,\n",
" 10.1010101 , 10.3030303 , 10.50505051, 10.70707071, 10.90909091,\n",
" 11.11111111, 11.31313131, 11.51515152, 11.71717172, 11.91919192,\n",
" 12.12121212, 12.32323232, 12.52525253, 12.72727273, 12.92929293,\n",
" 13.13131313, 13.33333333, 13.53535354, 13.73737374, 13.93939394,\n",
" 14.14141414, 14.34343434, 14.54545455, 14.74747475, 14.94949495,\n",
" 15.15151515, 15.35353535, 15.55555556, 15.75757576, 15.95959596,\n",
" 16.16161616, 16.36363636, 16.56565657, 16.76767677, 16.96969697,\n",
" 17.17171717, 17.37373737, 17.57575758, 17.77777778, 17.97979798,\n",
" 18.18181818, 18.38383838, 18.58585859, 18.78787879, 18.98989899,\n",
" 19.19191919, 19.39393939, 19.5959596 , 19.7979798 , 20. ])\n",
" t_events: None\n",
" y: array([[ 3.00000000e+00, 2.94289564e+00, 2.78785310e+00,\n",
" 2.55850312e+00, 2.27687667e+00, 1.96291516e+00,\n",
" 1.63409302e+00, 1.30523666e+00, 9.88467633e-01,\n",
" 6.93236946e-01, 4.26480505e-01, 1.92812066e-01,\n",
" -5.21327703e-03, -1.66825428e-01, -2.92730530e-01,\n",
" -3.84823958e-01, -4.45899551e-01, -4.79380989e-01,\n",
" -4.89082597e-01, -4.78998922e-01, -4.53117392e-01,\n",
" -4.15283908e-01, -3.69076135e-01, -3.17728315e-01,\n",
" -2.64072055e-01, -2.10504819e-01, -1.58984029e-01,\n",
" -1.11032079e-01, -6.77626414e-02, -2.99122951e-02,\n",
" 2.11578594e-03, 2.82076735e-02, 4.84876284e-02,\n",
" 6.32724837e-02, 7.30225570e-02, 7.82997565e-02,\n",
" 7.97282394e-02, 7.79589397e-02, 7.36422896e-02,\n",
" 6.74031022e-02, 5.98227645e-02, 5.14259816e-02,\n",
" 4.26711300e-02, 3.39462522e-02, 2.55671584e-02,\n",
" 1.77791644e-02, 1.07609976e-02, 4.63027520e-03,\n",
" -5.49504453e-04, -4.76159801e-03, -8.02780575e-03,\n",
" -1.04009265e-02, -1.19569125e-02, -1.27878889e-02,\n",
" -1.29959562e-02, -1.26872943e-02, -1.19678742e-02,\n",
" -1.09392829e-02, -9.69596064e-03, -8.32299987e-03,\n",
" -6.89464056e-03, -5.47365483e-03, -4.11100327e-03,\n",
" -2.84621315e-03, -1.70795463e-03, -7.15009954e-04,\n",
" 1.22640466e-04, 8.02548346e-04, 1.32852778e-03,\n",
" 1.70937387e-03, 1.95759286e-03, 2.08832134e-03,\n",
" 2.11822102e-03, 2.06464271e-03, 1.94483007e-03,\n",
" 1.77530419e-03, 1.57141171e-03, 1.34693899e-03,\n",
" 1.11392372e-03, 8.82507320e-04, 6.60920316e-04,\n",
" 4.55525977e-04, 2.70926087e-04, 1.10115997e-04,\n",
" -2.53345367e-05, -1.35074207e-04, -2.19767100e-04,\n",
" -2.80872771e-04, -3.20457795e-04, -3.40999749e-04,\n",
" -3.45225150e-04, -3.35964638e-04, -3.16023829e-04,\n",
" -2.88092971e-04, -2.54660687e-04, -2.17965318e-04,\n",
" -1.79955128e-04, -1.42270099e-04, -1.06239122e-04,\n",
" -7.28861887e-05],\n",
" [ 1.00000001e-01, 1.34350562e-01, 2.05200369e-01,\n",
" 2.60945888e-01, 2.62730256e-01, 1.97022797e-01,\n",
" 7.87991490e-02, -5.71595741e-02, -1.70086040e-01,\n",
" -2.29498465e-01, -2.24683268e-01, -1.68083819e-01,\n",
" -8.75293989e-02, -1.60182521e-02, 2.16793430e-02,\n",
" 1.66780349e-02, -2.25021472e-02, -7.49653919e-02,\n",
" -1.16628080e-01, -1.29381121e-01, -1.07842577e-01,\n",
" -5.95055614e-02, -1.64079135e-03, 4.66005148e-02,\n",
" 7.08554235e-02, 6.68512945e-02, 4.04292095e-02,\n",
" 4.99690262e-03, -2.45697805e-02, -3.72159335e-02,\n",
" -2.95955459e-02, -6.30027699e-03, 2.23953633e-02,\n",
" 4.51299486e-02, 5.36817111e-02, 4.57271499e-02,\n",
" 2.52196386e-02, 3.40446096e-04, -1.99238804e-02,\n",
" -2.91322714e-02, -2.55070499e-02, -1.20632246e-02,\n",
" 4.84234364e-03, 1.83167890e-02, 2.34051985e-02,\n",
" 1.88555840e-02, 7.06572375e-03, -6.99664069e-03,\n",
" -1.80138990e-02, -2.22054250e-02, -1.87016577e-02,\n",
" -9.44439487e-03, 1.63778313e-03, 1.03934001e-02,\n",
" 1.38860701e-02, 1.14678952e-02, 4.66288519e-03,\n",
" -3.47898807e-03, -9.75949186e-03, -1.19447335e-02,\n",
" -9.58096226e-03, -3.89439788e-03, 2.72638321e-03,\n",
" 7.80065101e-03, 9.61950297e-03, 7.86111216e-03,\n",
" 3.50905836e-03, -1.56657995e-03, -5.43793786e-03,\n",
" -6.79330617e-03, -5.39862291e-03, -2.03444893e-03,\n",
" 1.84455757e-03, 4.74888826e-03, 5.68095246e-03,\n",
" 4.47882504e-03, 1.76927414e-03, -1.31020139e-03,\n",
" -3.60337282e-03, -4.34770113e-03, -3.42898555e-03,\n",
" -1.34741949e-03, 1.00961839e-03, 2.74584556e-03,\n",
" 3.27914411e-03, 2.52952692e-03, 8.95328553e-04,\n",
" -9.32276127e-04, -2.25983325e-03, -2.64454771e-03,\n",
" -2.03234723e-03, -7.41819598e-04, 6.87698870e-04,\n",
" 1.71811472e-03, 2.01210318e-03, 1.53197929e-03,\n",
" 5.32592656e-04, -5.66956127e-04, -1.34886390e-03,\n",
" -1.55648511e-03],\n",
" [ 0.00000000e+00, -5.45041834e-01, -9.70308885e-01,\n",
" -1.28205400e+00, -1.48955475e+00, -1.60424275e+00,\n",
" -1.63887286e+00, -1.60680609e+00, -1.52143586e+00,\n",
" -1.39564347e+00, -1.24146645e+00, -1.06976797e+00,\n",
" -8.90079935e-01, -7.10480686e-01, -5.37567617e-01,\n",
" -3.76486172e-01, -2.31004514e-01, -1.03627050e-01,\n",
" 4.26627402e-03, 9.22673686e-02, 1.60772140e-01,\n",
" 2.10825729e-01, 2.43958864e-01, 2.62046295e-01,\n",
" 2.67175309e-01, 2.61525416e-01, 2.47278484e-01,\n",
" 2.26530001e-01, 2.01234504e-01, 1.73154941e-01,\n",
" 1.43834533e-01, 1.14580510e-01, 8.64576155e-02,\n",
" 6.02950676e-02, 3.66975254e-02, 1.60650764e-02,\n",
" -1.38506854e-03, -1.55920906e-02, -2.66261759e-02,\n",
" -3.46613065e-02, -3.99500121e-02, -4.28001565e-02,\n",
" -4.35525929e-02, -4.25634267e-02, -4.01877338e-02,\n",
" -3.67664577e-02, -3.26170320e-02, -2.80253636e-02,\n",
" -2.32414729e-02, -1.84767382e-02, -1.39030773e-02,\n",
" -9.65403231e-03, -5.82666839e-03, -2.48484200e-03,\n",
" 3.37229800e-04, 2.63064328e-03, 4.40767458e-03,\n",
" 5.69732192e-03, 6.54122715e-03, 6.98988886e-03,\n",
" 7.09902620e-03, 6.92677693e-03, 6.53092299e-03,\n",
" 5.96695504e-03, 5.28640586e-03, 4.53564106e-03,\n",
" 3.75517731e-03, 2.97917325e-03, 2.23540217e-03,\n",
" 1.54535481e-03, 9.24619852e-04, 3.83376603e-04,\n",
" -7.29797077e-05, -4.43170243e-04, -7.29327016e-04,\n",
" -9.36275223e-04, -1.07089182e-03, -1.14144098e-03,\n",
" -1.15704966e-03, -1.12719186e-03, -1.06128050e-03,\n",
" -9.68343588e-04, -8.56741545e-04, -7.34003123e-04,\n",
" -6.06683532e-04, -4.80310084e-04, -3.59364471e-04,\n",
" -2.47307094e-04, -1.46639820e-04, -5.89858064e-05,\n",
" 1.48062115e-05, 7.45556254e-05, 1.20629001e-04,\n",
" 1.53832636e-04, 1.75296997e-04, 1.86379057e-04,\n",
" 1.88570395e-04, 1.83415450e-04, 1.72449305e-04,\n",
" 1.57137647e-04],\n",
" [ 0.00000000e+00, 3.02924253e-01, 3.53810889e-01,\n",
" 1.64929950e-01, -1.59327509e-01, -4.77408687e-01,\n",
" -6.62482878e-01, -6.47845762e-01, -4.44093888e-01,\n",
" -1.34119406e-01, 1.69861891e-01, 3.66518593e-01,\n",
" 4.01772099e-01, 2.85062439e-01, 8.07033947e-02,\n",
" -1.22722651e-01, -2.46963689e-01, -2.51605839e-01,\n",
" -1.45035252e-01, 2.29443838e-02, 1.84151994e-01,\n",
" 2.78882658e-01, 2.77405678e-01, 1.87559359e-01,\n",
" 4.92383636e-02, -8.40393260e-02, -1.65505303e-01,\n",
" -1.72390828e-01, -1.10706335e-01, -1.17459573e-02,\n",
" 8.32389769e-02, 1.38182465e-01, 1.36047802e-01,\n",
" 8.19695329e-02, 6.62565168e-04, -7.58808531e-02,\n",
" -1.19911564e-01, -1.18578229e-01, -7.65421097e-02,\n",
" -1.31597775e-02, 4.64924266e-02, 8.09710789e-02,\n",
" 8.04309875e-02, 4.86572349e-02, 7.60265622e-04,\n",
" -4.37795849e-02, -6.85370832e-02, -6.60824817e-02,\n",
" -3.96811762e-02, -1.17704844e-03, 3.42383682e-02,\n",
" 5.39089324e-02, 5.22129449e-02, 3.18840906e-02,\n",
" 2.23413015e-03, -2.49000037e-02, -3.97500142e-02,\n",
" -3.80874159e-02, -2.21272886e-02, 8.12874064e-04,\n",
" 2.15625352e-02, 3.26056060e-02, 3.07933744e-02,\n",
" 1.79453643e-02, -1.65863908e-04, -1.64159084e-02,\n",
" -2.50044358e-02, -2.35751748e-02, -1.36074369e-02,\n",
" 3.56983930e-04, 1.27902397e-02, 1.92213050e-02,\n",
" 1.78898615e-02, 9.99717698e-03, -8.87720720e-04,\n",
" -1.04789069e-02, -1.53356352e-02, -1.41481826e-02,\n",
" -7.89619442e-03, 6.13055589e-04, 8.05167121e-03,\n",
" 1.17690522e-02, 1.07875654e-02, 5.90190512e-03,\n",
" -6.84654900e-04, -6.39023221e-03, -9.17624766e-03,\n",
" -8.31473074e-03, -4.45089264e-03, 6.86633090e-04,\n",
" 5.09427687e-03, 7.20537527e-03, 6.47864565e-03,\n",
" 3.43734879e-03, -5.61722521e-04, -3.96289968e-03,\n",
" -5.56283527e-03, -4.95661768e-03, -2.57072627e-03,\n",
" 5.31657862e-04]])"
]
},
"execution_count": 188,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"res"
]
},
{
"cell_type": "code",
"execution_count": 189,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x1a25a806d0>]"
]
},
"execution_count": 189,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3RcZ3nv8e8zo/tdsi6WZdnyRXbs3JwgjEMKpVxKEiBh9QRITgtJgJpDaYEeulpoeyht1zmLdh3u4UBNCCEpDZdAgwuBlgRYoQU7lnOxYzuOFceOZcuWbMu27peZ5/wxo0RWJFv27NFo9vw+a409e8/r/T7b2+un7Xe/s7e5OyIiEn6RTBcgIiJzQ4EvIpIjFPgiIjlCgS8ikiMU+CIiOSIv0wXMpLa21ltaWjJdhohIVtm+fftxd6+b7rN5G/gtLS20t7dnugwRkaxiZgdn+kxDOiIiOUKBLyKSIxT4IiI5QoEvIpIjFPgiIjki5cA3syIze8zMnjKzXWb2t9O0KTSz75hZh5ltNbOWVPsVEZELE8QZ/gjwene/ElgHXGdmG6a0eR/Q6+4rgc8B/xBAvyIicgFSDnxP6E8u5idfU++5fBPwzeT7B4A3mJml2vd0xmNx/s9Dezh8aigdmxcRyVqBjOGbWdTMngS6gZ+5+9YpTZqAQwDuPg6cBhZMs52NZtZuZu09PT0XVUtn7xD3P/YCv/+1LXT3DV/UNkREwiiQwHf3mLuvAxYD683ssilNpjubf9mTV9x9k7u3uXtbXd203ww+r5baUu65Yz3dfSO8+67H6B0YvajtiIiETaCzdNz9FPBL4LopH3UCzQBmlgdUAieD7HuyVyyt5q73tPH8iQFu+8Zj9A2PpasrEZGsEcQsnTozq0q+LwbeCDwzpdlm4Lbk+5uBn3uan6346pW1fPUPrmbXkTPc+fOOdHYlIpIVgjjDbwR+YWY7gG0kxvB/ZGZ/Z2Y3Jtt8HVhgZh3A/wQ+HkC/5/X6Sxr4ndV1bH7qCPG4nt0rIrkt5btluvsO4Kpp1n9y0vth4B2p9nUxblzXxMN7utl24CSvWv6y68QiIjkj9N+0feOaeorzo/zwqSOZLkVEJKNCH/glBXn87qUNPLSzi9HxeKbLERHJmNAHPsCNVy7i1OAY/9lxcXP7RUTCICcC/zWtdVSV5LP5SQ3riEjuyonAL8iLcMPljfzH7mMMjo5nuhwRkYzIicCHxLDO4GiMh/d0Z7oUEZGMyJnAX99SQ21ZIY/sOZbpUkREMiJnAj8SMV61rIb2A72ZLkVEJCNyJvABXtlSzeFTQ7p1sojkpJwK/LaWGgDaD6Ttvm0iIvNWTgX+msYKygrzeOx5Bb6I5J6cCvxoxLh6aTXbdIYvIjkopwIfYH1LNc8e6+fUoB6MIiK5JecC/6VxfM3WEZHcknOBv665ivyose2ghnVEJLfkXOAX5Ue5vKmSbbpwKyI5JucCH+CVy2rYefg0w2OxTJciIjJncjPwl9YwFnOePHQq06WIiMyZnAz8tpZqAA3riEhOSTnwzazZzH5hZnvMbJeZfWSaNq8zs9Nm9mTy9cnptjVXqkoKWFlfxhM6wxeRHJLyQ8yBceBj7v64mZUD283sZ+6+e0q7X7n7WwPoLxCXN1Xy6+eOZ7oMEZE5k/IZvrt3ufvjyfd9wB6gKdXtptuliyo4dmaEnr6RTJciIjInAh3DN7MW4Cpg6zQfX2NmT5nZT8zs0iD7vRhrF1UAsOvI6QxXIiIyNwILfDMrA74PfNTdz0z5+HFgqbtfCXwJeHCGbWw0s3Yza+/pSe8Dxy9dVAnAriNTSxURCadAAt/M8kmE/bfc/QdTP3f3M+7en3z/EJBvZrXTtNvk7m3u3lZXVxdEaTOqLM6nuaaY3Qp8EckRQczSMeDrwB53/+wMbRYm22Fm65P9nki171Rd2ljJ0xrSEZEcEcQsnWuBdwM7zezJ5Lq/BJYAuPtXgZuBD5rZODAE3OLuHkDfKbmsqYKf7jrKmeExKoryM12OiEhapRz47v6fgJ2nzZ3Anan2FbSJcfw9R87wquULMlyNiEh65eQ3bSdc+uJMHY3ji0j45XTg11cUUVdeqHF8EckJOR34kDjL10wdEckFCvxFFezr7tetkkUk9HI+8C9bVEks7uw92pfpUkRE0irnA1/fuBWRXJHzgd9cU0x5UZ7uqSMioZfzgW9mrG2sYHeXzvBFJNxyPvABVi8sZ9+xfubBl39FRNJGgU8i8PtHxjl8aijTpYiIpI0CH1jdUA7As8c0U0dEwkuBD7QmA3/v0f4MVyIikj4KfBL3xm+sLGLvUV24FZHwUuAnrV5Yzt5jOsMXkfBS4Cetbijnue5+xmPxTJciIpIWCvykVQ3ljMbiHDgxkOlSRETSQoGftHqhLtyKSLgp8JNW1pcRMdirqZkiElIK/KSi/CgtC0p5VnfNFJGQUuBPsqqhXGf4IhJaKQe+mTWb2S/MbI+Z7TKzj0zTxszsi2bWYWY7zOzqVPtNh9ULyzlwYkAPQxGRUAriDH8c+Ji7rwE2AB8ys7VT2lwPtCZfG4GvBNBv4FYvLMcdOrp14VZEwiflwHf3Lnd/PPm+D9gDNE1pdhNwrydsAarMrDHVvoO26sVbLGhYR0TCJ9AxfDNrAa4Ctk75qAk4NGm5k5f/UMDMNppZu5m19/T0BFnarLQsKKEgL6JxfBEJpcAC38zKgO8DH3X3qTelsWn+yMtuPu/um9y9zd3b6urqgipt1vKiEVbWlfGMzvBFJIQCCXwzyycR9t9y9x9M06QTaJ60vBg4EkTfQVvVUEaHzvBFJISCmKVjwNeBPe7+2RmabQbek5ytswE47e5dqfadDq0N5Rw5PUzf8FimSxERCVReANu4Fng3sNPMnkyu+0tgCYC7fxV4CLgB6AAGgTsC6DctWuvLAHiuZ4B1zVUZrkZEJDgpB767/yfTj9FPbuPAh1Ltay60Tnr6lQJfRMJE37Sdorm6mIK8iObii0joKPCnyItGWF5byj5duBWRkFHgT2NVQzn7dIYvIiGjwJ9Ga30Znb1DDIyMZ7oUEZHAKPCn0dowMVNHZ/kiEh4K/GlMzNTZp4eai0iIKPCnsbSmhPyoaRxfREJFgT+NxEydMs3UEZFQUeDPYGVDmc7wRSRUFPgzWFVfzqHeQYZG9fQrEQkHBf4MWhvKcNdMHREJDwX+DCZuoravW+P4IhIOCvwZLF1QSl7ENDVTREJDgT+DgrwIy2pLdeFWREJDgX8OrQ2amiki4aHAP4eV9eW8cHKQ4THN1BGR7KfAP4fW+jLiDvt7BjJdiohIyhT457Bq4p46mqkjIiGgwD+HltoSohHT069EJBQCCXwzu9vMus3s6Rk+f52ZnTazJ5OvTwbRb7oV5kVZuqCEZ3XhVkRCIOWHmCfdA9wJ3HuONr9y97cG1N+caa3XPXVEJBwCOcN390eBk0Fsa75prS/n4IlBRsY1U0dEsttcjuFfY2ZPmdlPzOzS6RqY2UYzazez9p6enjksbWatDWXE4s6B44OZLkVEJCVzFfiPA0vd/UrgS8CD0zVy903u3ububXV1dXNU2rm11idm6mgcX0Sy3ZwEvrufcff+5PuHgHwzq52LvlO1vK6UiKFxfBHJenMS+Ga20Mws+X59st8Tc9F3qoryoyypKaFDc/FFJMsFMkvHzO4HXgfUmlkn8DdAPoC7fxW4GfigmY0DQ8At7u5B9D0XWhvKdddMEcl6gQS+u996ns/vJDFtMyu11pfxi2e6GYvFyY/qu2oikp2UXrPQ2lDGeNw5cFz31BGR7KXAn4WJmTq6cCsi2UyBPwsr6sowQ+P4IpLVFPizUFwQpbm6RHfNFJGspsCfpdb6Mp3hi0hWU+DPUmtDOfuP9zMWi2e6FBGRi6LAn6VVDWWMxZyDJzRTR0SykwJ/liaefrX3qIZ1RCQ7KfBnaWV9GRGDvbqJmohkKQX+LBXlR2lZUMqzRxX4IpKdFPgXYFVDuW6TLCJZS4F/AVYtLOfAiQGGx/T0KxHJPgr8C7CqoYy4w3M9unArItlHgX8BVjfo6Vcikr0U+BegpbaU/KhpaqaIZCUF/gXIj0ZYUVemM3wRyUoK/Au0qqGcvZqaKSJZSIF/gVYvLOfwqSH6R8YzXYqIyAVR4F+g1voyAPZpWEdEskwggW9md5tZt5k9PcPnZmZfNLMOM9thZlcH0W8mrF6omToikp2COsO/B7juHJ9fD7QmXxuBrwTU75xrri6hKD+imToiknUCCXx3fxQ4eY4mNwH3esIWoMrMGoPoe65FIqZbLIhIVpqrMfwm4NCk5c7kuqykwBeRbDRXgW/TrPOXNTLbaGbtZtbe09MzB2VdnNUN5XT3jXByYDTTpYiIzNpcBX4n0DxpeTFwZGojd9/k7m3u3lZXVzdHpV24NY0VAOzpOpPhSkREZm+uAn8z8J7kbJ0NwGl375qjvgO3pjExU0eBLyLZJC+IjZjZ/cDrgFoz6wT+BsgHcPevAg8BNwAdwCBwRxD9ZsqCskLqywvZrcAXkSwSSOC7+63n+dyBDwXR13yxprGCPV26cCsi2UPftL1Iaxor6OjuY3Q8nulSRERmRYF/kdY0ljMWczq69QUsEckOCvyLtFYzdUQkyyjwL9Ky2lIK8iIKfBHJGgr8i5QXjbC6oZw9RxX4IpIdFPgpWJucqZOYhCQiMr8p8FOwprGckwOjdPeNZLoUEZHzUuCnYOIWC/oClohkAwV+Ci7RTB0RySIK/BRUFufTVFWsb9yKSFZQ4KcocYsFneGLyPynwE/R2sZy9vf0MzwWy3QpIiLnpMBP0aVNlcRdF25FZP5T4KfoisWVAOzsPJ3hSkREzk2Bn6KFFUXUlhWyQ4EvIvOcAj9FZsYViyvZefhUpksRETknBX4ALm+qpKO7n4GR8UyXIiIyIwV+AK5YrAu3IjL/KfADcHlT4sKtxvFFZD5T4AegvqKIhRVF7OzUOL6IzF+BBL6ZXWdme82sw8w+Ps3nt5tZj5k9mXy9P4h+55PLF1ey47DO8EVk/ko58M0sCnwZuB5YC9xqZmunafodd1+XfN2Var/zzRVNlezvGaBveCzTpYiITCuIM/z1QIe773f3UeDbwE0BbDerXJ78AtbTh3XhVkTmpyACvwk4NGm5M7luqv9mZjvM7AEza55uQ2a20czazay9p6cngNLmzhWLqwA0H19E5q0gAt+mWTf1mX//BrS4+xXAw8A3p9uQu29y9zZ3b6urqwugtLlTU1rA4upizdQRkXkriMDvBCafsS8Gjkxu4O4n3H3iOYBfA14RQL/zTuIbtwp8EZmfggj8bUCrmS0zswLgFmDz5AZm1jhp8UZgTwD9zjuXN1Vx8MQgpwZHM12KiMjLpBz47j4O/DHw7ySC/LvuvsvM/s7Mbkw2+7CZ7TKzp4APA7en2u98dGVz4sLtE4c0ji8i809eEBtx94eAh6as++Sk958APhFEX/PZuuYqohGj/cBJfmd1fabLERE5i75pG6CSgjwuW1RB+4HeTJciIvIyCvyAvWJpDU8eOsXoeDzTpYiInEWBH7C2lmpGxuPsOqLZOiIyvyjwA9a2tBpAwzoiMu8o8ANWX1HEkpoS2g+ezHQpIiJnUeCnQVtLNe0HenGf+oVjEZHMUeCnQdvSGk4MjHLgxGCmSxEReZECPw3aWhLj+NsOaFhHROYPBX4arKwro7I4n+26cCsi84gCPw0iEeMVS6vZpgu3IjKPKPDTpK2lmv09A5wc0I3URGR+UOCnyfqWGgC27j+R4UpERBIU+GlyZXMVZYV5PLrveKZLEREBFPhpkx+NcM2KBTz6bI/m44vIvKDAT6PXttZy+NSQ5uOLyLygwE+j17Qmnsv7q33Z9UB2EQknBX4atdSWsqSmhEef1Tj++cTiztBojL7hMUbH4xoGE0mDQJ54JTN7TWstDz5xmLFYnPxo7v58HRmP8UxXH88cPcOerj6ePz5Ad98IPX3D9A6OEYufHfARg6L8KNUlBdSUJl6NlUU0VhazqKqIpQtKWbqghPryQswsQ3slkl0U+Gn2mtY6vrX1BZ544RTrl9Vkupw5deD4AA/vOcav9h1n6/MnGB5LPBSmpCDKiroymqqKWNdcRU1pPoV5UQryIkTNGI3FGRmLMTga4+TgKL0DoxzvH2XXkTMc7x85q4/i/CjLaktZXlfKiroyVtYnXstqSynKj2Zit0XmrUAC38yuA74ARIG73P3TUz4vBO4FXgGcAN7l7geC6Hu+u2bFAqIR41f7enIi8PuGx/jxji4e2N5J+8HErSVW1JVyyyuX8KplNaxdVEFzdQmRyMWdlY+Mx+g6NcwLJwc5eGKAAycG2d/Tz47O0/x4ZxcTI0ERg+aaElYmfwisqC9jRV0py2vLqC4tCGp3RbJKyoFvZlHgy8CbgE5gm5ltdvfdk5q9D+h195VmdgvwD8C7Uu07G1QW57OuuYpH9x3nY7+7OtPlpM2J/hHu/q/nuffXB+kbGWdFXSl/cd0lvO3KRhZXlwTWT2FelJbaUlpqS4G6sz4bHovx/PEBOrr72dfdz3Pd/XR09/OrfccZjb30yMnqknyWLihlWW1iWKi5uoTmmhKaqotpKC8kL4eH3iTcgjjDXw90uPt+ADP7NnATMDnwbwI+lXz/AHCnmZnnyJW517bW8flHnqV3YDR0Z5d9w2Pc+fMO7v3NQYbHY1x/2ULe/5rlXNVcNedj60X5UdY0VrCmseKs9eOxOIdPDfFcTz/PdQ+w//gAB44PsGX/CR588jCT/xVGDOrLi1hYWUR9eSH1FYXUlhVSU1pAdUniVVGcR3lRPuVFeZQURCnKi170/1jOJx53RmPxxGv8pddYLM5I8vfEsjMaizE67ozF4i++RmPOeCzOeMwZi8eJxZyxuBOLxxmPO7GYMx534p54xeKJPhPL4D7pPbx4Md1f/OUcDCb+Vsxs0vvEejNLvrdJ6yYtW2IjZonjcna7l/6+p9vG1D4n6jl78extzNDsnGb6K5gp2XymPzFl9aKqYm57dcsFVDI7QQR+E3Bo0nIn8KqZ2rj7uJmdBhYAZ01fMbONwEaAJUuWBFDa/PCaVbV87uFneXRfDzeta8p0OYGIx50fPHGYT//kGU4MjPD2dU186HdWsLK+PNOlvUxeNJK8yFvK6y85+7OR8RiHe4c41DvE4d4huk4PceTUMMfODHPwxCDtB3tndT+k4vwohfkRCqIRCvIi5EcjRCNGXvIHwUT4TISmO8TcicVfeo3FEiH8UmD7yy5mByUvYkSnvsyIJH9PhKwRiSSCMWKTAnQiUGHGH+ruk6LNX8qzifXuifBzT76fbj0TwZlYjk9uM+UHj0/pc2I7TFo+q74ZFmYM5HOwmX5EXNjqs37oXLG4at4G/nT1T/1bm00b3H0TsAmgra0tNGf/6xZX0VBRyI93dIUi8A8cH+DPvvcU7Qd7Wddcxd23t3HF4qpMl3VRCvOiLK8rY3ld2YxtxmNxTg2N0TswysmBUfqGx+kbGaNveJzB0cTF5cGR8RfPwkfG42eF+OTgmghNM8iLRJIBC9FIhPyokRc18qOJHxwvvs976QfJ5B8oBXkvX5cfNQqT7/OSy/mRl7Y1EfSa2ZSbggj8TqB50vJi4MgMbTrNLA+oBHLm3sGRiHHD5Y18a+sL9A2PUV6Un+mSLoq78+1th/j7H+0mL2L8481XcPPVi9M2nDFf5EUj1JYlhnZEslkQV6e2Aa1mtszMCoBbgM1T2mwGbku+vxn4ea6M30946xWNjI7HeWRPd6ZLuSinB8f4w3u384kf7GRdcxU//ehreWdbc+jDXiRMUj7DT47J/zHw7ySmZd7t7rvM7O+AdnffDHwduM/MOkic2d+Sar/Z5qrmahori/jRji7eflV2DevsPdrHxvvaOXJqiL9+yxree+0yBb1IFgpkHr67PwQ8NGXdJye9HwbeEURf2WpiWOe+3xzkzPAYFVkyrPPQzi7+7HtPUVqYx/1/uIG2lvB/l0AkrDTheA695YpGRmNxHt59LNOlnJe78/9+2cEffetxVi8s50d/8lsKe5Esp8CfQ1c1V9FUVcyPd3RlupRzisWd//XDp/nHn+7lbVcu4tsbN9BQUZTpskQkRQr8OWRm3HD5Qh7d18PpobFMlzOt4bEYH7hvO/+85QU+8NvL+cK71lGYp3vSiISBAn+OveWKRYzFnJ/snH9n+f0j49z+jcd45JljfOpta/nE9Wt0cVYkRBT4c+zKxZVcsrCcb/7m4Ly65/vpwTH+4K6tbDvQy+feuY7br12W6ZJEJGAK/DlmZtz+6hb2dJ3hsefnx3fPTvSPcMvXtrD7yBm+/N+vzrppoyIyOwr8DLhpXRNVJfl8478OZLoUevpGuPVrW3j+eD933dbGdZctzHRJIpImCvwMKC6Icuv6JfzH7qN09mbuAefdfcPc+rUtHDo5xN23v5LXrqo7/x8SkaylwM+QP9iwFDPjvi0HM9J/95lhbt20hSOnhvjGHa/k1StqM1KHiMwdBX6GNFUV8+ZLG/j2Y4cYHB2f076Pnh7mXZu20HV6mHvuWM+G5QvmtH8RyQwFfgbd/uplnB4a47vbDp2/cUC6Tg9xy6bf0NM3wr3vXZ8Tj10UkQQFfga9sqWaDctr+Pwj++idxUM2UnX41BDv+qctHO8f5ZvvXa9bJYjkGAV+BpkZn7rxUvqGx/nMz/amta/njw/wjq/8mt7BUe5733pesbQ6rf2JyPyjwM+wSxZW8O4NS/mXrS+w68jptPTxzNEzvOOrv2F4PM79f7iBq5Yo7EVykQJ/HvjTN66iqqSAT23eFfi3b7cf7OVd/7SFvIjx3Q9cw2VNlYFuX0SyhwJ/HqgsyefP37yabQd6+V57Z2Db/benjnDr17ZQVZLP9/7HNaysn/m5rSISfgr8eeKdbc28alkNf/3g0/z6ueMpbcvd+dIj+/iT+59g3eIq/vWPrqW5piSgSkUkWynw54lIxNj07jaWLijhA/du55mjZy5qOycHRvngPz/OZ372LL93VRP3vX89NaUFAVcrItlIgT+PVJbkc89711NSGOX2u7dx+NTQBf35X+7t5s2ff5RHnjnGX95wCZ9555W6l72IvCilwDezGjP7mZntS/4+7fQPM4uZ2ZPJ1+ZU+gy7pqpi7rljPQMj41z/+Uf5l60vEI+f+0Lu3qN9fPj+J7j9G9uoKSnghx/6LTa+dgVmupe9iLzEUpkVYmb/CJx090+b2ceBanf/i2na9bv7BV0xbGtr8/b29ouuLdt1dPfz1w/uZMv+k1y1pIoP/vYKLm2qZFFl4lGDnb1D7Dpyhu8/3snPdh+jpCDKHde28Cevb6UoX2f1IrnKzLa7e9u0n6UY+HuB17l7l5k1Ar9099XTtFPgXwR351+fOMz//vEeTiS/iVtelAdA33Di/juVxfnccW0Lt13TQrXG6kVy3rkCPy/FbTe4exdAMvTrZ2hXZGbtwDjwaXd/MMV+c4KZ8XtXL+b6yxrZdeQ0e472sTd5MXdtYyVrGstZ01ihM3oRmZXzBr6ZPQxM91SMv7qAfpa4+xEzWw783Mx2uvtz0/S1EdgIsGTJkgvYfLgVF0Rpa6nRvW9EJCXnDXx3f+NMn5nZMTNrnDSk0z3DNo4kf99vZr8ErgJeFvjuvgnYBIkhnVntgYiIzEqq0zI3A7cl398G/HBqAzOrNrPC5Pta4Fpgd4r9iojIBUo18D8NvMnM9gFvSi5jZm1mdleyzRqg3cyeAn5BYgxfgS8iMsdSumjr7ieAN0yzvh14f/L9r4HLU+lHRERSp2/aiojkCAW+iEiOUOCLiOQIBb6ISI5I6dYK6WRmPcDBFDZRC6R2Y/nsk4v7DLm537m4z5Cb+32h+7zU3eum+2DeBn6qzKx9pvtJhFUu7jPk5n7n4j5Dbu53kPusIR0RkRyhwBcRyRFhDvxNmS4gA3JxnyE39zsX9xlyc78D2+fQjuGLiMjZwnyGLyIikyjwRURyROgC38yuM7O9ZtaRfM5uKJlZs5n9wsz2mNkuM/tIcv2sHiyfzcwsamZPmNmPksvLzGxrcp+/Y2ahe9ajmVWZ2QNm9kzymF8T9mNtZn+a/Lf9tJndb2ZFYTzWZna3mXWb2dOT1k17bC3hi8l822FmV19IX6EKfDOLAl8GrgfWArea2drMVpU248DH3H0NsAH4UHJfPw484u6twCPJ5bD5CLBn0vI/AJ9L7nMv8L6MVJVeXwB+6u6XAFeS2P/QHmszawI+DLS5+2VAFLiFcB7re4Drpqyb6dheD7QmXxuBr1xIR6EKfGA90OHu+919FPg2cFOGa0oLd+9y98eT7/tIBEATif39ZrLZN4G3Z6bC9DCzxcBbgLuSywa8Hngg2SSM+1wBvBb4OoC7j7r7KUJ+rEncvr3YzPKAEqCLEB5rd38UODll9UzH9ibgXk/YAlQlnzY4K2EL/Cbg0KTlzuS6UDOzFhKPjdzKlAfLAzM9WD5bfR74cyCeXF4AnHL38eRyGI/5cqAH+EZyKOsuMyslxMfa3Q8D/xd4gUTQnwa2E/5jPWGmY5tSxoUt8G2adaGed2pmZcD3gY+6+5lM15NOZvZWoNvdt09ePU3TsB3zPOBq4CvufhUwQIiGb6aTHLO+CVgGLAJKSQxnTBW2Y30+Kf17D1vgdwLNk5YXA0cyVEvamVk+ibD/lrv/ILn62MR/8c71YPksdS1wo5kdIDFc93oSZ/xVyf/2QziPeSfQ6e5bk8sPkPgBEOZj/UbgeXfvcfcx4AfAqwn/sZ4w07FNKePCFvjbgNbklfwCEhd5Nme4prRIjl1/Hdjj7p+d9NF5Hyyfrdz9E+6+2N1bSBzbn7v775N4VvLNyWah2mcAdz8KHDKz1clVbwB2E+JjTWIoZ4OZlST/rU/sc6iP9SQzHdvNwHuSs3U2AKcnhn5mxd1D9QJuAJ4FngP+KtP1pHE/f4vEf+V2AE8mXzeQGNN+BNiX/L0m07Wmaf9fB/wo+X458BjQAXwPKMx0fWnY33VAe/J4PwhUh/1YA38LPAM8DdwHFIbxWAP3k7hOMUbiDP59Mx1bEkM6X07m204Ss8HZGEIAAAA4SURBVJhm3ZdurSAikiPCNqQjIiIzUOCLiOQIBb6ISI5Q4IuI5AgFvohIjlDgi4jkCAW+iEiO+P9PhI4ABEzN2gAAAABJRU5ErkJggg==\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": 193,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x1a27187ed0>]"
]
},
"execution_count": 193,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhbV53/8ffRYkuyJct7vGZfmzRrt5SWtrQ0oUDZ2w7QFiiFeWDYZoAyhfkxzAzrDDCdYYBSoGUvFNqGNqW0pUugS/ZmcxY7ifd9kRdZtiWd3x+SXMeRba2xc/V9PU+eWNK1z1Wu8/Hx95x7jtJaI4QQwvhMs30CQgghzg0JfCGEyBAS+EIIkSEk8IUQIkNI4AshRIawzPYJTKWoqEgvWLBgtk9DCCHOK3v27OnSWhdHe23OBv6CBQvYvXv3bJ+GEEKcV5RS9VO9JiUdIYTIEBL4QgiRISTwhRAiQ0jgCyFEhpDAF0KIDCGBL4QQGUICXwghMkRGB34gqDnc4uHnL53m1ca+2T4dIYRIqzl741W6/cfjR/j1zkYGR/wAXLqogN/cedksn5UQQqRPRga+1prf7GxkSUkuH7h8AQ/va6auc3C2T0sIIdIqI0s6PUOjDIz4eevacm5cV8GKeS7aPD6CQdn9SwhhXBkZ+PU9XgDmFzoAqHDbGAtougZHZvO0hBAirTIy8Bu6zwz8crcdgOa+4Vk7JyGESLeMDPz6bi9KQWX+mYHf0uebzdMSQoi0yszA7xlinsuGzWoGJga+9PCFEMaVkYHf0O2lusAx/thls5CTZZaSjhDC0DIy8Ot7vOP1ewClFOVuu/TwhRCGlnGBPzTip3NghPmFOWc8X+620+qRGr4QwrgyLvAbJk3JjJAevhDC6DIu8OsjUzILzuzhV7htdA+N4hsLzMZpCSFE2mVc4Df0DAFQHaWHDzJTRwhhXBkX+PXdXtwOK3l26xnPy1x8IYTRZVzgN/R4mV/gOOv58jzp4QshjC3jAr++20v1pBk6AKV52SgFLR4JfCGEMWVU4I8FgjT3DUft4WdbzBTnZksPXwhhWBkV+M29wwSC+qwB24jQ1Eyp4QshjCmjAj+yLPKCKCUdgAqZiy+EMLCMCvyG7tCUzMk3XUWUu2009w2jtWyEIoQwnowK/PpuLzariRJndtTXy/LsjPiD9AyNnuMzE0KI9MuswO8JrZKplIr6uszFF0IYWUoCXym1RSl1TClVq5S6K8rrn1FKHVFKHVBKPaOUmp+KduMVWhY5ev0eQjV8kKmZQghjSjrwlVJm4HvAVmAVcItSatWkw/YBm7TWFwIPAd9Mtt14aa1p7D1zHfzJyt02QG6+EkIYUyp6+BcDtVrrk1rrUeA3wI0TD9BaP6u19oYfvgxUpqDduPT7/HhHA+OhHk1BThbZFpMEvhDCkFIR+BVA44THTeHnpvIh4IloLyil7lRK7VZK7e7s7EzBqb2mLbzWfalr6sBXSoWnZkoNXwhhPKkI/GgjoFHnNSql3gdsAr4V7XWt9b1a601a603FxcUpOLXXtIbr8mV5Uwc+QJnbJjV8IYQhpSLwm4CqCY8rgZbJBymlrgXuBt6qtR5JQbtxae8P9drnzRD4bkcWHu/YuTglIYQ4p1IR+LuApUqphUqpLOBmYNvEA5RS64EfEgr7jhS0GbfI9oUlzukDP89uxTMsgS+EMJ6kA19r7Qc+DjwJ1AC/1VofVkp9RSn11vBh3wJygd8ppfYrpbZN8eXSps3joyg3myzL9G85z26l3zcmd9sKIQzHkoovorXeDmyf9Ny/TPj42lS0k4y2ft+M9XsAl83KWEAzPBbAkZWSfx4hhJgTMuZO2zaPb9oZOhGRnbCkrCOEMJqMCfxWT4w9fHuoV98/7E/3KQkhxDmVEYE/PBrAMzw24wwdkB6+EMK4MiLw2yJTMuMo6fRL4AshDCYjAj/Wm64gNGgL0sMXQhhPRgR+rDddwYQevk8CXwhhLBkR+JGbrmIJfKctNGgrPXwhhNFkROC3eXy4bJaY5tVbzCZysy0S+EIIw8mIwA9NybTHfHye3SrTMoUQhpMRgd/e76M0hnJOhNMmPXwhhPFkROC3enyUxTAlMyKyno4QQhiJ4QN/LBCka3AkpgHbCJfdKvPwhRCGY/jA7xgYQevYZuhE5EngCyEMyPCB3xa+6SquHr5N1sQXQhiP4QM/Mgc/lrtsI/LsVoZGA4wFguk6LSGEOOcMH/iRzctjWUcnIi+8YuaAT6ZmCiGMIyMC32Y1jS+ZEAuXrJgphDAgwwd+a3/opiulVMyfIytmCiGMyPCB3+7xxVXOAenhCyGMyfCB3+rxxTVDB2TFTCGEMRk68P2BIB0Dse1lO5HseiWEMCJDB/6x9gHGApqVZc64Pk82QRFCGJGhA39vQx8AG6rz4/o8m9VEltkkK2YKIQzF0IG/r76XotxsKvNjXxoZQCmFyy4rZgohjMXQgb+3oZcN1e64pmRGuGTFTCGEwRg28LsHRzjd7WXD/PjKOREumyygJoQwFsMG/r4E6/cRsmKmEMJoDBv4ext6sZgUF1bmJfT5eXZZMVMIYSyGDvxV5S5sVnNCny+DtkIIozFk4PsDQV5t9CRczoHINod+tNYpPDMhhJg9hgz8o20DDI8FWF/tTvhruGxWAkHN0GgghWcmhBCzJyWBr5TaopQ6ppSqVUrdFeX1K5VSe5VSfqXUu1LR5nT2NfQCiQ/YgqyYKYQwnqQDXyllBr4HbAVWAbcopVZNOqwBuB34VbLtxWJvQx/FzvhvuJpI1tMRQhiNJQVf42KgVmt9EkAp9RvgRuBI5ACt9enwa+dkz8BkbriKcEkPXwhhMKko6VQAjRMeN4WfmxVdgyPUd3uTKueA9PCFEMaTih5+tG50QlNblFJ3AncCVFdXJ3Qyjiwz33/vBlaVuxL6/IjIipn9sq+tEMIgUtHDbwKqJjyuBFoS+UJa63u11pu01puKi4sTOhlHloWta8qYX5iT0OdHSA9fCGE0qQj8XcBSpdRCpVQWcDOwLQVfd1bl2kK//EjgCyGMIunA11r7gY8DTwI1wG+11oeVUl9RSr0VQCl1kVKqCXg38EOl1OFk2003s0nhtFlk0FYIYRipqOGjtd4ObJ/03L9M+HgXoVLPeUUWUBNCGIkh77RNFZdN1sQXQhiHBP40ZMVMIYSRSOBPw2W3yL62QgjDkMCfhssmPXwhhHFI4E/DabMyNCI9fCGEMUjgTyPXZmFw1E8waLw18Uf9QR7Z14zHK7/BCJEpUjIt06hcNgtaw9CoH2d4qQWjeGhPE//88EGcNgsfff1ibt+8gJxs+XYQwsikhz+N3HAADhhwPZ0nD7dR4bZzycJCvvXkMV7/redo7PHO9mkJIdJIAn8akeUVBg1Wx+/3jfFiXRc3XFjGfbdt4rcfuYyuwREe3d8826cmhEgjCfxpRMo4Awa7+erZox2MBTTXX1AKwMULC1hTkcezxzpn+cyEEOkkgT8No5Z0/ny4naLcbNZXvbZnwNUrStjX0Evv0OgsnpkQIp0k8KfhMmBJxzcW4LljHVy3qhST6bWtDK5eXkxQwwsnpJcvhFFJ4E8jUsM3Ug//xbouhkYD4+WciAsr3RTkZPHcOSjrHGnp5zc7G9DaeNNdhZjLZB7eNCIlnUEDBf6Th9pxZlvYvLjojOfNJsXrlxXz/PFOAkGN2ZT4fsDTqe0Y4JYfvYxneIzT3V4+v2V5UnsPCyFiJz38aeRkWVDKOIO2gaDm6Zp2rl5RQpbl7Et/9YoSeoZGebWpLy3tt3l83PrjnVjNJt6+voIfPF/HPc/UpqUtIcTZpIc/DZNJkZtlYcAgNfw99b10D41y/QXzor5+5dIiTAqeO9qR9Cbwk3mGx7j9pzvxDI/x4EcuY1WZC7NJ8Z2nj2PPMnHnlYtT2p4Q4mzSw5+B02YxTElnX0MvAJcvKYz6utuRxYbq/LRMz/zq4zXUdgzyg/dvZHVFHiaT4hvvvJAb1pTxtSeO0tQrN30JkW4S+DPItVkMM2hb2zFIsTMbtyNrymOuXlHCwWYPHQO+lLU74g+w/WArN66r4Iqlr21ObzYp7tq6Aq3h4b1y05cQ6SaBPwOnzWqYaZl1nYMsLs6Z9pirlocC+W+1XSlrd8fxLgZG/Lx5bdlZr1UVOLhsUSG/39sks3aESDMJ/BnkZlsMMWirtaa2Y5AlJbnTHrdinoucLDP7GlI3cPvYgRby7FZet6Qo6uvv2ljJ6W4ve+p7U9amEOJsEvgzyLUZY9C2c3CEfp+fxcXTB77ZpFhb5U5Z4PvGAjx1pJ0tF8zDao7+7bZl9TwcWWYe2tOUkjaFENFJ4M/AZZBB27qOIYAZe/gA66rc1LT24xsLJN3u88c7GRoNcMOFZ5dzInKyLbxpTRmPH2hleDT5NoUQ0UngzyBU0jn/A7+2cxCIPfD9Qc3hFk/S7T5+oJV8h5XLFkefGRTxzg2VDIz4+fORtqTbFEJEJ4E/A6fNyvBYAH8gONunkpS6jkFysszMc9lmPHZdtRsg6bKObyzA0zXtbFldNmU5J+KShQVU5tulrCNEGkngz2B8eYXzvI5f1znI4pLcmJYxKHHaqHDb2deYXOA/e7QD72iAN09TzokwmRTv2FDJX2u76OhP3ZTQyfY29HKkpT9tX1+IuUwCfwbpXEAtGNR8809HefxAa9qnJNZ2DM44YDvRumo3+5Ps4W8/1EZhThaXLCyI6fgb1pShNTxd05FUu9EMjvi5++GDvOP/XuSG/9nBF/5wUJaCFhlHAn8GrjQG/qnuIf7vuTo+9qu9fOTne9LWsx0c8dPq8cVUv49YX+WmuW844RuwgkHNjhOdXLW8BMsM5ZyIZaW5VBc4eCrFdfyX6rq5/jsv8KudDXz4ioV8YPNCfru7kav/6zkeO9CS0raEmMsk8GeQmx3a9SodJZ2DTaFB0fdfOp/nj3dy7befZ9fpnpS3czI8YBtPD399uI6faC//SGs/fd4xXrd0+sHaiZRSXLuylL/VdTOUon/v9n4fH7x/F9kWEw99dDN337CKf3nLKrZ/4goWFObwmd++Ov7vI4TRSeDPwDm+CUrqb7460OTBZjXx/96yiic+eQVmk+LXrzSkvJ3ajsgMnenvsp3ogvI8LCbF/gTr+DtOhO7UvXyKm62mct2qUkb9QXakaCOWb/7pGIGg5v4PXMzG+a8tCLd8npN7b92IzWLi878/QDCY3pJaY4+Xe545wSP7mvEMn/838onzk6yWOYN01vAPNveFgtVsYlFxLhuq8zmUgqmQk9V1DmIxKeYXxh74NquZlWWuhGfq/K22i+WlTkqcM88KmuiiBfm4HVb+fCQ0uycZB5s8/H5vEx99/WKqCx1nvV7itPGlN6/isw8d4Jev1PP+yxYk1V40h5o9/OD5OrYfbCXyM8ViUly2uJBPXbvsjB9CQqSbBP4MnGna1zYQ1Bxq7uemi6rGn7ugIo9nj3UwPBrAnmVOWVu1HYNUFzpmnBo52boqN3/Y2xT3hii+sQA7T/fwvkvmx3uqWMwmrllewl+OduAPBGOu/0+mteYrjx2mKDeLj1099dLL79pYybZXW/j6E0e5ZmUpFW57Qu1Fs/1gKx/71V5ysyx8+IpF3LZ5AW39Pp483Ma2/S3c8qOXuefmdUn/YItGa01N6wA7T3Wz63QvdZ2DrCxzsXF+PhcvLGBZqTPlbYq5T0o6M3DaQjX8VAd+Xecgw2MBLqzMG39udbmLoIaattROG6zrHGJJHPX7iPXVboZGA+MloVjtqe9l1B/kiqXxlXMirltVSp93jN1JrK3zxKE2dp3u5R/fuHz8GkajlOKrb1+DBr70yKGE25tsT30Pn3pwPxuq8/nbF67hC29aSbnbzobqfL6wdSWPf+IKLih38fe/3MsDL55OWbsALX3D3P7TXbzpnh18+Y9H2NfQS4nLxo4TXXzxkUO88Tsv8MH7d8V9XeOhtaZ7cIRXG/vYdbqH7sGRtLUlYpeSHr5Sagvw34AZuE9r/fVJr2cDPwM2At3ATVrr06loO91sVhNmk0p5Df9AeMD2jMCvCH18uNmTsg1IxgJBTncN8cZVpTMfPMm6qtDA7d6GXpbPi71HuONEFxaT4uIYp2NOdsWyYrLMJp4+0s6li2If9I3wB4J87YkaVsxz8p5NVTMeX1Xg4FPXLuWr24/yt9quuMcdJjvZOcgdD+ymwm3nR7duwhXlB05BTha/uuNS/uHX+/h/2w7T6x3lU9cuS6pdrTUP7mrk3x+vIRDU3P2mlWxdM4/KfMf46409wzx+sJX/e7aWLd99gfddOp/PvHFZ1HOMl28swJOH2/jd7ib21PcyPGlpjnyHlVXlLt60powb1pRNu0x3ogZH/HT0++gbHsPjHcNqNlHqyqbEacNlt2T8dppJB75Sygx8D7gOaAJ2KaW2aa2PTDjsQ0Cv1nqJUupm4BvATcm2fS4opdKyCcrBpj4cWWYWFr3W8y7Ls1GQk8Wh5tT18Ou7vfiDOq4ZOhELi3Ioys3m5ZPd3HJxdcyf97faLjZU55OTndi3V262hc1LCnmqpp27b1gZ93/S7YfaaOwZ5ke3boq5FHXrZQt44MV6vvZEDds+9jpMCe7p2+cd5faf7sKkFPd/4CIKcqYONXuWmR+8bwN3/eEg3336BLnZFu64YlFC7frGAvzT717lsQOtXLqogG++c+1Z4xZKKaoLHfz9VYt596ZKvv3UcX720mmePNzGV9+xhquXlyTUdke/j+8/X8dDe5oY8PmpKrBz88VVVBc4qHDbsVpMnOwcoq5zkFdOdnP3w4f48rbDvH5ZCbdcXMVVy0sS3kO5zzvK88c72Xmqhz31vRxrH2CqW1ryHVbWVrlZW+lmw/x8NlS7p/3tbzpaazoGRjjaNkB99xAN3V6aeocZGvUzPBpgxB/EbjXjsltw2a2U5dmoyndQXeCgutBBWZ49bftGTycVPfyLgVqt9UkApdRvgBuBiYF/I/Dl8McPAf+rlFL6PFkAPR3r6Rxo9rC6PO+Mi66U4oJyV0oHbuviWENnMqUUmxcX8mJdN1rrmIK3d2iUQy0ePp1kb/W6VaXc/fAhaloHWFXuivnztNbc+0Idi4pzeMOK2APMZjXzj29cxmd++yp/PNDCjesq4j5nrTV3P3KIVs8wD37kspgGyS1mE99454UMjwb498dryM22cHMcP1wBeoZG+fDPdrOnvpfPbVnOR69cPOMPrKLcbL769jW8Z1MVn/3dq3zgp7t418ZKvrB1BYW52TG12zkwwg+er+MXL9fjD2recmEZN11UzSULC85q/+rlob+11hxu6efR/c08vK+Fp2vaKc+z8Z6Lqtiyeh7LS53Tfp9pranrHOT54108faSdnad7CAQ1zmwL66rdXH/BPBYW5ZDnsJJntzLqD9IxMEJHv4/j7QPsb+zj+eOdaA0mBavKXWyozmf5PCcr5jmZX5hDbraFbIsJpRTDowF6vKO0eXzUdgxwvH2Q4+0DHGnpp3vCjXvZFhOV+XZcdit2qxmnzcLwWICWPh81rQO09/vwT5gJlmUxjf9AnOeyUerKxmW3kmUxkWU2UeqycXUc37+xSkXgVwCNEx43AZdMdYzW2q+U8gCFwBm7bCil7gTuBKiuju+bPp2cNmtKl0geCwQ50tLP+y49e1BzdUUe9+04yYg/QLYl+YHbk52hVTIXzbDxyVQ2Ly5k26st1HYMsjSGgb6XTnajdfzTMSfburqML287zO/3NrGqfFXMn/fyyR4ONffztXesibuX/rZ1Ffxoxym+9eQxtqyeF/e//7ZXW3j8QCufvX55XCU5s0nxnZvWMTTq5wsPHyTbauLt6ytj+ty6zkE+dP8uWjw+vvd3G6ZdlTSadVVuHvvE67jnmRP84PmTbD/Yyu2bF/DhKxaRP8VvJzWt/fzkr6d4dH8L/mCQt6+v5BNvWBLTDzilFKsr8lhdkcfntqzg6SPt/GpnA999+gTfffoEFW47Vy0vprrAQYkrtDtb18AIzX3D1Hd7eamum7bwDYpLS3L5yJWLuG5VKRdWumPuMQ+O+NnX0MuuUz3sPN3D7/c0MTRplVaLSWExK3xjZ66hZbOaWFKSyxtWlrCqzMWKMhcLi3Iozs2e9vvNHwjS1u+jocdLfbeX011DnOoaotXjo6a1n67BESbODF5f7Z6zgR/tXU7uucdyDFrre4F7ATZt2jRnev/OFG+CcqJ9kBF/8Iz6fcTq8jzGApoT7YPjNf1kNPR4KczJSvhX10hwv1jXHVPg7zjRhTPbwtoo7y0eBTlZXLuylEf2NXPX1hUxzzD60Y6TFOVm8fb18ffQTSbFF7au4Naf7OTnL9XHVV5p9QzzpUcOsaHazUeujL8sk2Ux8f33buQD9+/k0w++Sl3HEJ+5btm0IfLIvmb++eGD2Kxmfv3hS9g4P7Exk2yLmc9ev4K3r6/gv5+p5fvP1/Gzl+q5dFEhK+Y5WTbPicc7Sk3bAIeaPRxo8mC3mnnPRZV88PKFLEqgXAhgNZvYuqaMrWvKaO/38ezRDp6u6eDR/S1Rb3QsdWWzaX4Br1taxOuWFFFVcPZU21jkZlu4Ymnx+HabwaCmuW+Y4+0DNPZ4GRoNMDTixx/UuB1WCnOyKMrNZklJLpX5joRKMRazicp8B5X5DjZHmTTmDwTxjgUY9QcZ9QcxpWmsIRWB3wRMHBmrBCbfrx45pkkpZQHygNTfUpomuTZLSvd4Pdgcmtu+Jkqgr64IlS8ONXtSEviNPd6E/2NAaECzMt/Oi3Vd3LZ5wbTHBoOa5451cNniwoSnU070ro2VPHGojb8c7eD6C+bNePyJ9gH+crSDz1y3DJs1sd+OrlxWzBVLi7jnmRO8+cJy5uXNfB9BMKj53EMH8Ac1337PuoTfuz3LzM8+eAlfeuQQ//tsLcfaB/jOTevGF/CLGBrx869/PMxvdzdx0YJ87rllPWV5yU8nXVLi5H9uWc8/XLOEe184yf7GPp491kEg3PV02SysLHNx19YV3HxRVUoHXUtdNm6+uJqbL65Gax0afB0Yoc87SlFuNvPybCn5jTcak0lRVeBI6v9JsixmE64U/J+ZsZ0UfI1dwFKl1EKgGbgZ+LtJx2wDbgNeAt4F/OV8qd9D6G7bus7UlXQONHlwZltYEOVX4OoCB06bJWV1/IYeL2vDs20StXlxIX861DbjfPz9TX20enx8bsvypNqLeP2yYoqd2Ty0pymmwL9vxylsVlPUUlk8/vWtF3DDPX/lsw+9ygMfuHjG0tD3n69jx4ku/uPtq1lQlFjpLCLLYuLr71zDyjIn//Z4DZd+9RkuX1I4PrD51JF2dpzoZMQf5ONXL+FT1y5NyQ/XiZaVOvnPd68FQoPBp7qGyAsPPJ6LWS6hiRLWhH8rFVNLOvDDNfmPA08Smpb5E631YaXUV4DdWuttwI+Bnyulagn17G9Ott1zKTc7tbN0DoZ779GCZHzgNgUzdfyBIM19w7wlyubh8bh8SRG/3d3EkZZ+1kxTqnniYCtWs+INK+OfAhqNxWziHesr+PFfT9E1OELRNIOJjT1eHt7XzHsuqpx2ZkwsFhXncvcNK/niI4f4+cv10/5m8+ThNr715DFuXFfO38U52DoVpRS3X76QNZVuHtrTyHPHOnnycDsQmsn1nk1VvG19Rcqm7k4ncse1MIaUzMPXWm8Htk967l8mfOwD3p2KtmZDKgdtR/wBalr7+eDlC6c8ZnV5Hj9/uT6pO00BWj0+AkFNdZK/ql4Wngv/Yl3XlIGvtWb7wTauWFqckjndEe/aWMkPXzjJI/uap62pf/1PRzGZ4GNXL0lJu++9pJqna9r56vYaLl9SFHWWU01rP59+cD9rq9x8450Xprz3u3F+Phvn56O15kTHIGOBIKvKXBk/l1wkTu60jYHTZmHUH2TEn/x+q409w4wF9LQ3Mq2uyGPEH6QuPMMm8ba8AEnXJktcNpaU5PJiXfeUxxxs9tDcN8zW1TOXXuKxtNTJ2io3v9vdNOWeAbtO9/D4gVY++vrFKallQ6iX/c13Xogjy8zHfrn3rO0ej7cPcMcDu3HaLPzo/RsTHjOI9VyWlTq5oDxPwl4kRQI/BuO7XqWgrNMenlI2XTBNHLhNRkM48JPt4UOojr/zVA+j/uhbPW4/2IbFpLgugTt6Z/LujZUcax/g2WNnb4wSDGq+8scjzHPZuDOB2THTKXHZ+M5N62gf8PHm//krn3lwP08daeeD9+/ijd95gT7vKD+6dRMlMWwbKcRcIIEfA2cKV8xs84QCf7rZHwuLcrFbzUkP3Db0eLGYVEp6vZsXFzE8FuDVprNXzwyVc1rZvKQoLbfLv2NDBSvmOfnkb/aftf7LH/Y1c7DZw+e3LseRlfq1AK9aXsLzn72aj1y5mMcOtvLhn+1mf2Mfn752GTs+fw0XViY3IC7EuSSBH4NU7msbuWlkus3EzSbF4pKcpEs6DT1eKvNTcwv3pYsKMCl47NWzd4g63NJPQ4+XN6W4nBPhyLJw322byLaYuOOBXfR5R9Fa80xNO19/4ihrq9zcuDb+efexyrNbuWvrCp79p6v44fs38uJd1/DJa5cmPTgsxLkmyyPHIJUrZrZ5fOTZrTMuf7ygMGd8gbVEJTsHfyK3I4ubLqrmF680cMsl1ayY99rMjScOtWI2Kd4Yw9TJRFXmO/jh+zdyy72vcPtPd+EbC3C0bYAKt53/eNvqhNe+iUeF257S5ZOFONekhx+D10o6yd9t2+rxURbDzTyLinJo6vVOWTOPRUMKAx/gc9cvx2Wz8KVHDo0PoJ7uGuLBXU1cuqgg7T3ejfML+No71rC/sQ9/UPNf717Lc5+9KiU3qAmRCaSHH4NUlnTa+32UxjDIt6Aoh6AOhXYiC58N+Mbo9Y6lZMA2Ij8ni7u2ruDzvz/I7/c2s3F+Pjff+xJBrfnSm2Nf7yYZ79xYyaWLCylz2c5Jr14II5HAj0FKB237fVwQw+qPC8N3bJ7uGkoo8Bt7hoHUzNCZ6N0bq3hwVyNf216DxawYC2h+9eFLzijxpJuUVYRIjJR0YpBrS00PfywQpGtwJKYefiTwT3UlNnCbytFRqvwAAA95SURBVCmZE5lMin9722p6vaOzEvZCiMRJDz8G2RYzWRZT0j38joERtJ5+SmaE25GF22HlVHdigZ+qm66iuaA8j1/ccQmVbkfUzcGFEHOTBH6MUrFEcpsnVGaJJfAh1Ms/leDUzIYeL3n20CYQ6bB5cXLr3Qshzj0p6cTIabMkXdJp84Q2cp5uDv5ECwtzOJ1gD7+hx5vyco4Q4vwmgR+jXFvy2xy2jS+rEFvgLyjKodXjY3g0/jV8QnPwZXBTCPEaCfwYpWKJ5DbPMNkWU8xllvGZOnH28gNBTVPv8Kxu6CCEmHsk8GPktFnpT7aG3z/CvDg2kZg4NTMe7f0+RgNBKekIIc4ggR8jZ3YqavjDMdfvgfHdk+KdqdOYpimZQojzmwR+jJwpquHHOkMHQmWkYmd23DN10jUHXwhxfpPAj5HLbmXAN0YwmNhWvFpr2j0jcQU+JDZTp6HHi0lBudyRKoSYQAI/Rnl2K0FNwlsd9gyNMhoIxlXSAVhQ5OBUlzeuz6nrHKS6wIE1xZtbCyHOb5IIMYrMrOkfTmzgNpZ18KNZWJRL1+BIXDd91XUktv6OEMLYJPBjFAl8T6KBH8NOV9EsLArV4U/H2Mv3B4Kc6hpicbEEvhDiTBL4MYoEfp83yR5+3IEfCu5YZ+o09Q4zGgiyWHr4QohJJPBjFNmrNdEefrvHh0lBcW52XJ83P7w4WawzdSJ7vkoPXwgxmQR+jMZ7+MOjCX1+q8dHsTMbS5wDqTarmfI8G6e6Bmc+mNCALcASCXwhxCQS+DFyO5Ks4ff74h6wjVhS6uRYe2yBX9sxSFFuNnmO9KySKYQ4f0ngx8hmDa2Jn8ygbbz1+4hVZS5qOwZi2t+2rnOQJSU5CbUjhDA2Cfw45NmteJIYtE20h7+yzMlYQI+Xa6aitaa2Y1Dq90KIqCTw4+C2WxPq4Q+N+Bnw+ZmXl9idr6vKQlsI1rT2T3tc5+AI/T6/zMEXQkQlgR+HvAQD/7UpmfHN0IlYWJRDlsU0Y+DXdYRm8kjgCyGikcCPQ57dmtA8/I7+0E5Xpc7ESjoWs4nlpU5qWgemPa62U6ZkCiGmJoEfhzxHYj38zsFQ4Bc7E+vhQ6iOX9Paj9ZTL95W1zGII8sc845aQojMklTgK6UKlFJPKaVOhP/On+K4Pyml+pRSjyXT3mzLs1sTWkunI1zSKUmwhw+wssxF99AonQMjUx5T1xkasI11gxUhRGZJtod/F/CM1nop8Ez4cTTfAt6fZFuzLs9uZWDEjz8w8/TIiToHR8gym3DZLQm3vTI8cHtkmjp+Xceg1O+FEFNKNvBvBB4If/wA8LZoB2mtnwGmL0CfB9yRFTPj3Ailc2CEYmd2Uj3vlfMiM3Wi/zMOjfhp8fhYXCxz8IUQ0SUb+KVa61aA8N8lyXwxpdSdSqndSqndnZ2dSZ5a6kXuXu3zxre8QiTwk227wm2fcqbOyU6ZoSOEmN6MNQal1NPAvCgv3Z3qk9Fa3wvcC7Bp06bEtpZKo0SXSO4cGKEqBdsNRgZuo6ntDPX8ZYaOEGIqMwa+1vraqV5TSrUrpcq01q1KqTKgI6VnN8fk2RNbMbNzYIQN86OOZ8dlZZmLZ4914hsLYLOaz3ittmMQs0kxv1BKOkKI6JIt6WwDbgt/fBvwaJJfb05LpIc/FgjSPTRKSZIlHQgFfiCoORFlIbXnj3dyQbmLLIvMtBVCRJdsOnwduE4pdQK4LvwYpdQmpdR9kYOUUjuA3wFvUEo1KaWuT7LdWZFI4HcPhur9ydbw4bWZOpPLOrUdgxxq7ueta8uTbkMIYVyJzxMEtNbdwBuiPL8buGPC4yuSaWeuGA/8OO62jcybj3fjk2jmFzhwZJnPmpr56P5mTAoJfCHEtOT3/zhkWUw4ssz0xdHD7xgI33SV4EqZE5lMipVlLv5a2zV+L4DWmkf3t7B5cVFK2hBCGJcEfpziXUBtvIefgpIOwB2vW0htxyC/eLkegL0NfTT0eLlxnfTuhRDTk8CPU6KBX5SblZL2t6yexxVLi/ivPx+nc2CER/c3k20xsWV1tJmzQgjxGgn8OMW7CUrHwAhuh5Vsi3nmg2OglOLLb70Anz/Avz12hMcOtHLtylKcNtnSUAgxPQn8OCXSw0/FgO1Ei4tz+fAVi9j2ags9Q6NSzhFCxEQCP05uh5W+4diXVugcHKHEldrAB/j4NUsoz7ORZ7dy1fKkVrQQQmSIpKZlZqJ4e/gdAz42Vid/l+1kjiwLP/nARQz4/HKzlRAiJhL4ccqzW/GNBaMubzCZ1jolC6dNZUV4BU0hhIiFdA3jlOcIzbaJZSOUwRE/vrFgUhufCCFEqkjgxyme5RU6UjwHXwghkiGBH6d4Aj/VN10JIUQyJPDjFNn1qi+GufiRwE/FSplCCJEsCfw4SUlHCHG+ksCPk9sRX0kny2wa/yEhhBCzSQI/TpElDGJZMTMVm5cLIUSqSODHyWxSOG2WmKZldgz4KJJyjhBijpDAT4DbYaXPO/PyCulYR0cIIRIlgZ+AWJdX6ErTOjpCCJEICfwExBL4/vDm5dLDF0LMFRL4CXDbs2YctO0eGkVrpIcvhJgzJPAT4LJbZxy0TeXm5UIIkQoS+AmIlHS01lMe0+oJbV5eKhuLCyHmCAn8BBQ7sxkLaHqnWV6hvnsIgOoCx7k6LSGEmJYEfgIq8+0ANPZ4pzymoceL02YZvzNXCCFmmwR+AqryQ732xt7pA39+oUPushVCzBkS+AmoKoj08IenPKah2yvlHCHEnCKBnwCnzYrbYZ2yhx8Iahp7vVQX5JzjMxNCiKlJ4CeoKt8xZQ2/rd/HWEBLD18IMadI4CeoqsBOU2/0kk5khs78Qgl8IcTcIYGfoKp8B829wwSDZ8/Fb+gO9fylhy+EmEsk8BNUWeBgNBCkfcB31msNPV4sJkVZntx0JYSYO5IKfKVUgVLqKaXUifDf+VGOWaeUekkpdVgpdUApdVMybc4VVeG5+NHKOvU9Xirz7VjM8vNUCDF3JJtIdwHPaK2XAs+EH0/mBW7VWl8AbAG+q5RyJ9nurKsKl2uiDdw29njHXxdCiLki2cC/EXgg/PEDwNsmH6C1Pq61PhH+uAXoAIqTbHfWVbinnotf3+2VAVshxJyTbOCXaq1bAcJ/l0x3sFLqYiALqJvi9TuVUruVUrs7OzuTPLX0slnNlLqyz5qL7/GO4RkeY77MwRdCzDGWmQ5QSj0NzIvy0t3xNKSUKgN+DtymtQ5GO0ZrfS9wL8CmTZumXopyjog2F78h/FhKOkKIuWbGwNdaXzvVa0qpdqVUmda6NRzoHVMc5wIeB76otX454bOdY6oKHOw81XPGc/U9MgdfCDE3JVvS2QbcFv74NuDRyQcopbKAh4Gfaa1/l2R7c0pVvp1WzzBjgdd+Yanvlh6+EGJuSjbwvw5cp5Q6AVwXfoxSapNS6r7wMe8BrgRuV0rtD/9Zl2S7c0JlgYOghpa+1wZuG3u8FOVmkZs94y9PQghxTiWVSlrrbuANUZ7fDdwR/vgXwC+SaWeuGl8muWeY+YWhQdp6WSVTCDFHyZ1BSRhfJnnCTJ2GHgl8IcTcJIGfhLI8OxaTGp+pM+oP0uIZprpQpmQKIeYeCfwkmE2KcredxvDyCk29XrSWRdOEEHOTBH6SqgrsNPZ40Vrzl6OhWakyJVMIMRfJVJIkVeU7eOJQGx+8fxfPHutkbZWbNRV5s31aQghxFgn8JFUVOPAMj/HKqR6+eMNKbt+8QFbJFELMSRL4SbpxXTl93lFu27yAynwp5Qgh5i4J/CRV5ju4+4ZVs30aQggxI6k9CCFEhpDAF0KIDCGBL4QQGUICXwghMoQEvhBCZAgJfCGEyBAS+EIIkSEk8IUQIkMorefmXuFKqU6gPokvUQR0peh0zheZ+J4hM993Jr5nyMz3He97nq+1Lo72wpwN/GQppXZrrTfN9nmcS5n4niEz33cmvmfIzPedyvcsJR0hhMgQEvhCCJEhjBz49872CcyCTHzPkJnvOxPfM2Tm+07ZezZsDV8IIcSZjNzDF0IIMYEEvhBCZAjDBb5SaotS6phSqlYpdddsn0+6KKWqlFLPKqVqlFKHlVKfDD9foJR6Sil1Ivx3/myfa6oppcxKqX1KqcfCjxcqpV4Jv+cHlVJZs32OqaaUciulHlJKHQ1f88uMfq2VUp8Of28fUkr9WillM+K1Vkr9RCnVoZQ6NOG5qNdWhdwTzrcDSqkN8bRlqMBXSpmB7wFbgVXALUopo25H5Qf+UWu9ErgU+Fj4vd4FPKO1Xgo8E35sNJ8EaiY8/gbwnfB77gU+NCtnlV7/DfxJa70CWEvo/Rv2WiulKoBPAJu01qsBM3AzxrzW9wNbJj031bXdCiwN/7kT+H48DRkq8IGLgVqt9Umt9SjwG+DGWT6ntNBat2qt94Y/HiAUABWE3u8D4cMeAN42O2eYHkqpSuAG4L7wYwVcAzwUPsSI79kFXAn8GEBrPaq17sPg15rQFqx2pZQFcACtGPBaa61fAHomPT3Vtb0R+JkOeRlwK6XKYm3LaIFfATROeNwUfs7QlFILgPXAK0Cp1roVQj8UgJLZO7O0+C7wOSAYflwI9Gmt/eHHRrzmi4BO4KfhUtZ9SqkcDHyttdbNwH8CDYSC3gPswfjXOmKqa5tUxhkt8FWU5ww971QplQv8HviU1rp/ts8nnZRSbwY6tNZ7Jj4d5VCjXXMLsAH4vtZ6PTCEgco30YRr1jcCC4FyIIdQOWMyo13rmST1/W60wG8CqiY8rgRaZulc0k4pZSUU9r/UWv8h/HR75Fe88N8ds3V+aXA58Fal1GlC5bprCPX43eFf+8GY17wJaNJavxJ+/BChHwBGvtbXAqe01p1a6zHgD8BmjH+tI6a6tkllnNECfxewNDySn0VokGfbLJ9TWoRr1z8GarTW357w0jbgtvDHtwGPnutzSxet9Re01pVa6wWEru1ftNbvBZ4F3hU+zFDvGUBr3QY0KqWWh596A3AEA19rQqWcS5VSjvD3euQ9G/paTzDVtd0G3BqerXMp4ImUfmKitTbUH+BNwHGgDrh7ts8nje/zdYR+lTsA7A//eROhmvYzwInw3wWzfa5pev9XAY+FP14E7ARqgd8B2bN9fml4v+uA3eHr/QiQb/RrDfwrcBQ4BPwcyDbitQZ+TWicYoxQD/5DU11bQiWd74Xz7SChWUwxtyVLKwghRIYwWklHCCHEFCTwhRAiQ0jgCyFEhpDAF0KIDCGBL4QQGUICXwghMoQEvhBCZIj/D1jNiR8SGfdUAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(traj[:, 1])"
]
},
{
"cell_type": "code",
"execution_count": 192,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x1a27080a10>]"
]
},
"execution_count": 192,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXgcV53u8e+vW6tleZW8y7uzOIsdUEIWYLISxzc3BoYlgUvCFsNcAszcGYbkhoGBwDNhwgMXLhm4TgjbBMIyhJjEZGVNQkiUxPESZ5EdO5ZlWYu1WbvUv/tHl+yO07Jkd0stV72f59GjruqjOqdU9tulU6fqmLsjIiLhF8t1A0REZGwo8EVEIkKBLyISEQp8EZGIUOCLiEREXq4bcCRlZWW+cOHCXDdDROS48fTTTze6e3m698Z14C9cuJCqqqpcN0NE5LhhZruGei8rXTpmdoeZ1ZvZliHeNzP7lplVm9kmM3tDNuoVEZGRy1Yf/g+AVUd4/zJgWfC1FvhOluoVEZERykrgu/ufgP1HKLIG+JEnPQFMMbPZ2ahbRERGZqxG6cwFdqcs1wTrXsfM1ppZlZlVNTQ0jEnjRESiYKwC39KsS/sQH3df5+6V7l5ZXp72QrOIiByDsQr8GqAiZXkeUDtGdYuICGMX+OuBq4PROmcDre6+d4zqFhERsjQO38x+CpwPlJlZDfAFIB/A3b8LbABWA9VAJ/ChbNR7PHB3tu1t5/HtjSycXsJ5S8soLojnulkiEkFZCXx3v2qY9x34RDbqOl509Q7wtQdf5P4tdexp6Tq4vjAvxnlLy/i785dw5sJpOWyhiETNuL7T9njV3TfAtT+q4rHtjVx88kw+ddFS3rKsnB0NHTzywj5+u7mOK9c9wQ2XncRH3rwIs3TXtEVEskuBn2U9/QN87MdP89j2Rm551wre9cZ5B9+bM6WYNy8r4x8uOYHP/OI5vnzfNp55tZl/f9cKJhbqUIjI6NLTMrOofyDBJ+58hj++1MC/veO014R9qklF+Xz3f7yRGy47ifu31PHBO56ku29gjFsrIlGjwM+inz75Kg9vq+eLV5zClWfNP2JZM+Njf7OEb111BlW7mvlfP99IIqH5hUVk9Cjws6Sjp59vPvIyb1o0javPWTDin7v89DncuPpkNmyu4ysbto1iC0Uk6tRxnCW3//kVGg/0ctvVJx31RdiPvmURe1q6+N6jr7Bg+gSuPmfh6DRSRCJNZ/hZ0High3V/2s6qU2ZxxvypR/3zZsa/XL6cC0+awU33Ps+WPa2j0EoRiToFfhZ8+3fVdPcn+MyqE495G/GY8bV3r2DqhAI+fdezdPXqIq6IZJcCP0O793dy51938Z7KCpaUT8xoW9NKCvjGe1eyo7GDm+57PkstFBFJUuBn6BdP19CfcD554dKsbO+8pWWsfetifvLXV7l/S11WtikiAgr8jLg792zcw3lLypgzpThr2/3HS07k1LmTuPHuzTR39GZtuyISbQr8DDy7u4VdTZ2sWTknq9styItxy7tW0NrVx033qmtHRLJDgZ+Be57dQ2FejFWnzsr6tk+ePYm/O38Jv3p2D394sT7r2xeR6FHgH6O+gQS/2bSXi5fPpLQof1TquO7CpSwpL+HGu7fQ0dM/KnWISHQo8I/Roy83sr+jl3esTDs1b1YU5sX56t+eTm1rF7c88OKo1SMi0aDAP0Z3P7uHKRPyeesJozvvbuXCaXzg7AX86C872VTTMqp1iUi4ZSXwzWyVmb1oZtVmdn2a9z9oZg1mtjH4+mg26s2VAz39PPh8HZefPpuCvNH/zPynS09k+sRCbrx7CwN6wJqIHKOM08rM4sCtwGXAcuAqM1uepujP3H1l8HV7pvXm0iPb9tHdl2DNKHbnpJpUlM/nL1/O5j2t/OcTu8akThEJn2ycnp4FVLv7DnfvBe4C1mRhu+PWoy83MmVCPm84hufmHKvLT5/NW5aVccsDL7KvrXvM6hWR8MhG4M8Fdqcs1wTrDve3ZrbJzH5pZhVDbczM1ppZlZlVNTQ0ZKF52eXuPL69iXMWTyceG7upCc2Mm9acSu9Agi9pbL6IHINsBH661Du8o/k3wEJ3Px14GPjhUBtz93XuXunuleXlo3tB9FjsaupkT0sX5y4tG/O6F5aV8Inzl3Lfpr388aXx92EoIuNbNgK/Bkg9Y58H1KYWcPcmd+8JFm8D3piFenPise2NAJy3ZHpO6v/4+YtZXFbC5+/ZomkRReSoZCPwnwKWmdkiMysArgTWpxYws9kpi1cAx+3UTo9XNzF7chGLykpyUn9hXpwvv/1UdjV18h+/r85JG0Tk+JRx4Lt7P3Ad8ADJIP+5u281sy+Z2RVBsU+Z2VYzew74FPDBTOvNhUTCeXx7I+cuKTvqWa2y6dylZbzjjLl854/bqa4/kLN2iMjxJSuDyN19g7uf4O5L3P0rwbrPu/v64PUN7n6Ku69w9wvc/YVs1DvWnt/bRnNnH+ctzU13Tqr/vfpkivPjfO7Xm3HX2HwRGZ7utD0Kjw/23+fggu3hyksLuf6yk3lix35+UVWT6+aIyHFAgX8UHqtuYkl5CTMnFeW6KQBceWYFZy2axpfve556jc0XkWEo8Eeotz/Bk6/sHxdn94NiMePmd55Gd3+Cf/3N1lw3R0TGOQX+CG3c3UJX3wDnLhk/gQ+wuHwin75oGRs21/HAVk2JKCJDU+CP0JOvNAFwzuLcX7A93Nq3LuakWaX8y6+30NrZl+vmiMg4pcAfoU01rSwuK2HyhNGZ7CQT+fHklIj7O3r5/PotuW6OiIxTCvwR2lTTyunzJue6GUM6bd5kPnnhMu7ZWMu9m2qH/wERiRwF/gjUt3VT19bNafOm5LopR/SJC5awomIKn/v1Fj1RU0ReR4E/Apv3tAKM6zN8gLx4jK+/ZwXdfQN89r826YYsEXkNBf4IPFfTSszglDmTct2UYS0pn8gNl53MH15s4PY/v5Lr5ojIOKLAH4HNNS0sm1HKhIK8XDdlRK4+ZwGrTpnFzfe/wFM79+e6OSIyTijwh+HubN7TymnjvDsnlZnx7+8+nYqpxVz3k2doPNAz/A+JSOgp8Iext7WbxgO9477//nCTivK59f1voKWzj0/f9awmPxcRBf5wNtW0AHDa3OMr8AFOmTOZm9acymPVTXzxN1t1EVck4o6PTukc2lTTSl7MOHn2+L9gm857zqzg5fp2bvvzK1RMncC1b12c6yaJSI4o8IexeU8rJ84qpSg/nuumHLMbLjuZPS1dfGXDNuZOLWb1abOH/yERCR116RyBu4/7O2xHIhYzvv6elbxxwVT+/mcbeay6MddNEpEcyMoZvpmtAr4JxIHb3f3mw94vBH5EcvLyJuC97r4zG3WPplf3d9La1cdpc8f3HbYjUZQf57arK3nfbU/woR88xf/7wBu54MQZuW5WTiUSzoHefjp7Bujo7aerd4DegQS9/Qn6B5yEO07ygz9mlvyKJZ9dlB+PkRczivJjFObFKcyPUZQfpygvTn7ccjoFpshQMg58M4sDtwKXADXAU2a23t2fTyn2EaDZ3Zea2ZXAV4H3Zlr3aNtUc3zcYTtS00oK+Om1Z/OBO/7Kx370NN9+3xm87ZRZuW7WqGjt6mNXUwe793dR09zJ3tZu6tu7qW/roamjl5bOXlq7+hiNwUvxmFGcH6e4IE5xfpwJBcnXEwriFOfnMaHg8HVxilLKF+XHKcqPURR8kBTmxSnMi1EQfA1+4OTH7eAHz3j8gHF33CHhTsLBSS4Prhv8ME04kPr+4M+SLAvJ90g5Vv6aeo6tfUP9ymyIBUtZMEtb5LDt27BlhmqHYaPyoMZsnOGfBVS7+w4AM7sLWAOkBv4a4F+D178Evm1m5uN82MiW2lYK4jFOmFma66ZkzdSSAu786Nlcc8eT/M87n+Gmt5/KVWfNz3Wzjllnbz8v1LWzbW8b2/a28fK+A2xv6HjdvQelhXmUTypkRmkhp8yZxNQJBUyZkM+konxKCvMoKUwGbUFejMJ4jLx4jJil/qd1BhIwkHD6E8m/AHr6E/QOJOjpG6C7P/jeN0BX3wBdvQm6+vrp7B2gs3eArt4BOnv7aTrQS3dfyrq+gawMmY3HjLyYEY8ZcTNiwevBfYhZMkTMkuFz+AfEYMAmXx8evofeT7iTSBwK42SYHwrt1yzLMSubWEjV5y7O+nazEfhzgd0pyzXAm4Yq4+79ZtYKTAde15lsZmuBtQDz5+c2iF7ed4DF5SUU5IXrUsfk4nx+/JGz+MRPnuWGX23mud0tfHHNKRTmje8L0+7OjsYOqnbuZ+PuFp59tYWX9rUfDJfSwjxOmFXKhSeVs6R8IgvLSqiYOoF504qZVDT+HmsNyX3qHUjQ3ZdIflj0DtDdP3BwubtvgN6DHywJ+gYSh7qdEk5ff4K+hDOQSC4PDDgDQSgPDJ5ZO0FIHwrxQ/UfOsNMfhAMvg4+HCy5FDOI2aEPjFgs2cVlJD9sCD5Q4jEOrregG8yMgx88ZqnvJ18TlB1cN/izg68HG5XavsE2Hnx9lH/gDHWq6Sm/HR/iL4rUN4b6XHvNzx7hvHaod4pHaZBINgI/3a/68P0YSZnkSvd1wDqAysrKnJ4nvLSvnZUVx3//fTqlRfl8/4Nn8vWHXuTW329nW107//fKM5g/fUKum3ZQIuG8VN/OX7Y38cSOJqp2NtPU0QskP7RWVEzhkuUzOXXuZJbPnsS8qcXjsmvjSMws6LKJM7l4fH4oSXhkI/BrgIqU5XnA4Q9kHyxTY2Z5wGRgXD/kpbO3n5rmLt5TWTF84eNUPGZ85tKTOG3uFP7pF89xyTf+yCcvXMq1b12cs7P9V5s6ebS6kce2N/LE9qaDAV8xrZjzT5zBmQunUrlwGkvKS467cBfJtWwE/lPAMjNbBOwBrgTed1iZ9cA1wF+AdwG/G+/999X1BwA4YebEHLdk9K06dRYrKiZz073P87UHX+JXz+7hny89kUuWz0r+uT6KGtp7+MuOJh6vbuTR6kZqmrsAmDWpiL85sZxzl5Rx9uJpzJs6fv7yEDleZRz4QZ/8dcADJIdl3uHuW83sS0CVu68Hvgf82MyqSZ7ZX5lpvaPtpX3JwF86IzwXbI9k9uRi/uP9b+T3L9bzxfVb+fh/PsOC6RP48HmLePsZc7PW3bC3tYuqnc08tXM/f9nexMvBB2tpUR7nLJ7OtW9ZzHlLy3QGLzIKbDyfaFdWVnpVVVVO6v63327j+4/u5PkvXUpePFwXbYczkHAe2FrHbX/ewbOvtpAXM85cOI2LTp5xsDuldJiLoD39A+xp7mJHQwfP723j+do2NtW0UNuanIlrQkGcyoXTOGfxdM5ZMp1T50yK3O9ZZDSY2dPuXpnuPT1aYQiDI3SiGELxmLH6tNmsPm02G3e38MDWOn63rZ4v37ftYJlZk4qYNbno4JhygAM9/RzoSQ49rGvrPjhSwQwWTi/hjAVTuXbBVCoXTOOk2aXkR/B3K5JLCvwhhHmEztFYWTGFlRVT+Oyqk9jT0sWWPa1U1x9ge/0BGjt66ezpp7mzDwMmFuZRPrGQE2aWMn/aBCqmTmBhWQknzSqlpFD/1ERyTf8L04jCCJ1jMXdKMXOnFHPpKbluiYgcC/1NnUaURuiISHQo8NOI2ggdEYkGBX4aL9e3UxCPsXAc3XUqIpIpBX4aUR6hIyLhpURL46V97Sydof57EQkXBf5hBkfohOmRyCIioMB/HY3QEZGwUuAfRiN0RCSsFPiH0QgdEQkrBf5hXmnoYMH0CRqhIyKho1Q7zK6mThZML8l1M0REsk6BnyKRcHbt71B3joiEkgI/xb72brr7Eiwo0xm+iIRPRoFvZtPM7CEzezn4PnWIcgNmtjH4Wp9JnaNpZ2MngM7wRSSUMj3Dvx54xN2XAY8Ey+l0ufvK4OuKDOscNbuaOoDkZB0iImGTaeCvAX4YvP4h8PYMt5dTO5s6yY8bc6YU57opIiJZl2ngz3T3vQDB9xlDlCsysyoze8LMjvihYGZrg7JVDQ0NGTbv6Oxq6qBi2gTiMU2eLSLhM+yMV2b2MDArzVs3HkU989291swWA78zs83uvj1dQXdfB6yD5CTmR1FHxnY2dao7R0RCa9jAd/eLh3rPzPaZ2Wx332tms4H6IbZRG3zfYWZ/AM4A0gZ+rrg7u5o6OHvxtFw3RURkVGTapbMeuCZ4fQ1wz+EFzGyqmRUGr8uA84DnM6w36xrae+jsHdAZvoiEVqaBfzNwiZm9DFwSLGNmlWZ2e1DmZKDKzJ4Dfg/c7O7jLvB3NgVDMjUGX0RCatgunSNx9ybgojTrq4CPBq8fB07LpJ6xsPPgkEyNwReRcNKdtoFdTR3kxYy5GpIpIiGlwA/sbOpk3tRiPSVTREJL6RbY2dihp2SKSKgp8Bkcktmp/nsRCTUFPtDU0cuBnn6d4YtIqCnwOfTQtEUakikiIabA59BjkReoS0dEQkyBT/IMP2Ywb6oCX0TCS4FPckjmnCnFFOTp1yEi4aWEA2qaO5k/TWf3IhJuCnygprlLd9iKSOhFPvC7+waob+9R/72IhF7kA7+2pQuAeVN1hi8i4Rb5wN+jwBeRiIh84Nc0JwN/rgJfREJOgd/cSTxmzJpUlOumiIiMqowC38zebWZbzSxhZpVHKLfKzF40s2ozuz6TOrOtprmL2ZOL9FhkEQm9TFNuC/BO4E9DFTCzOHArcBmwHLjKzJZnWG/W7GnuUv+9iERCRoHv7tvc/cVhip0FVLv7DnfvBe4C1mRSbzYlx+BrSKaIhN9Y9GPMBXanLNcE69Iys7VmVmVmVQ0NDaPasN7+BPvau3WGLyKRMOwk5mb2MDArzVs3uvs9I6jD0qzzoQq7+zpgHUBlZeWQ5bJhb2sX7hqSKSLRMGzgu/vFGdZRA1SkLM8DajPcZlZoSKaIRMlYdOk8BSwzs0VmVgBcCawfg3qHVdOcfA5+hR6rICIRkOmwzHeYWQ1wDnCfmT0QrJ9jZhsA3L0fuA54ANgG/Nzdt2bW7Oyoae4iZjBrssbgi0j4DdulcyTufjdwd5r1tcDqlOUNwIZM6hoNe5q7mD25mHyNwReRCIh00umxyCISJREP/E6N0BGRyIhs4PcNJKhr0xh8EYmOyAZ+XWs3CdfE5SISHZEN/N3BkEyNwReRqIhs4A/edKUuHRGJikgHvhnMnqzAF5FoiGzg72nuYmZpEQV5kf0ViEjERDbtapo71X8vIpES2cCva+tmjm66EpEIiWTguzt7W7uZrWfoiEiERDLwmzp66e1PKPBFJFIiGfh1rd2ARuiISLREMvBrW5Jj8HWGLyJREsnA3zt4hj9FgS8i0RHZwM+PG2UlhbluiojImIlo4Hcxc1IRsVi6+dVFRMIp0ykO321mW80sYWaVRyi308w2m9lGM6vKpM5s2NvSzRxdsBWRiMn0DH8L8E7gTyMoe4G7r3T3IT8Yxsreti7134tI5GQ6p+02ALPjp2skkXDqWrs1cbmIRM5Y9eE78KCZPW1ma49U0MzWmlmVmVU1NDRkvSGNHT30Dbi6dEQkcoY9wzezh4FZad660d3vGWE957l7rZnNAB4ysxfcPW03kLuvA9YBVFZW+gi3P2KHbrrSGb6IRMuwge/uF2daibvXBt/rzexu4CxG1u+fdbUtustWRKJp1Lt0zKzEzEoHXwNvI3mxNyf2tgZ32eqirYhETKbDMt9hZjXAOcB9ZvZAsH6OmW0Iis0EHjWz54Angfvc/f5M6s1EXWs3BfEY00sKctUEEZGcyHSUzt3A3WnW1wKrg9c7gBWZ1JNNtcEIneNpZJGISDZE7k7butYuXbAVkUiKXODXtmimKxGJpkgF/kDC2demm65EJJoiFfhNB3roTzhzFPgiEkGRCvxazXQlIhEWqcDfG8x0pS4dEYmiaAV+cIavi7YiEkURC/wuCvNiTJ2Qn+umiIiMuUgFfm1rN7N105WIRFSkAr+utVsXbEUksiIY+LpgKyLRFJnATySc+vZuZkxS4ItINEUm8Pd39tI34MyaVJjrpoiI5ERkAn9wpiuNwReRqIpM4O9rSwb+THXpiEhERSjwewAFvohEV6YzXt1iZi+Y2SYzu9vMpgxRbpWZvWhm1WZ2fSZ1Hqu6tm7MoLxUffgiEk2ZnuE/BJzq7qcDLwE3HF7AzOLArcBlwHLgKjNbnmG9R21fazdlEwvJj0fmjxoRkdfIKP3c/UF37w8WnwDmpSl2FlDt7jvcvRe4C1iTSb3Hoq6tm1nqzhGRCMvm6e6Hgd+mWT8X2J2yXBOsS8vM1ppZlZlVNTQ0ZK1x+9q6makhmSISYcMGvpk9bGZb0nytSSlzI9AP3JluE2nW+VD1ufs6d69098ry8vKR7MOIJANfZ/giEl15wxVw94uP9L6ZXQNcDlzk7umCvAaoSFmeB9QeTSMz1d03QHNnn7p0RCTSMh2lswr4LHCFu3cOUewpYJmZLTKzAuBKYH0m9R6thnYNyRQRybQP/9tAKfCQmW00s+8CmNkcM9sAEFzUvQ54ANgG/Nzdt2ZY71GpG7zpSnfZikiEDdulcyTuvnSI9bXA6pTlDcCGTOrKxMHHKugMX0QiLBKD0gcfq6DAF5Eoi0zgF+bFmFSc0R80IiLHtUgEfl1bD7M0taGIRFwkAn9fq8bgi4hEI/Db9VgFEZHQB767U9eqxyqIiIQ+8Fu7+ujpT6hLR0QiL/SBP3jTlaY2FJGoC33ga6YrEZGk8Ae+7rIVEQEiEPiDXTozdNFWRCIuEoE/raSAwrx4rpsiIpJToQ/8+rZuZmjichGR8Ad+XVu3RuiIiBCBwN/X1sPMUgW+iEioA79/IEHjgR5NfCIiQoYToJjZLcB/B3qB7cCH3L0lTbmdQDswAPS7e2Um9Y5U44Fe3FEfvogImZ/hPwSc6u6nAy8BNxyh7AXuvnKswh6gvj2Y2lBj8EVEMgt8d38wmLMW4AlgXuZNyp5Dd9nqDF9EJJt9+B8GfjvEew48aGZPm9naLNZ5RINTG+oMX0RkBH34ZvYwMCvNWze6+z1BmRuBfuDOITZznrvXmtkM4CEze8Hd/zREfWuBtQDz588fwS4Mrb6tGzOYXlKQ0XZERMJg2MB394uP9L6ZXQNcDlzk7j7ENmqD7/VmdjdwFpA28N19HbAOoLKyMu32Rqq+vYeyiYXkxUM9GElEZEQySkIzWwV8FrjC3TuHKFNiZqWDr4G3AVsyqXek9rVp4hMRkUGZnvp+Gygl2U2z0cy+C2Bmc8xsQ1BmJvComT0HPAnc5+73Z1jviOimKxGRQzIah+/uS4dYXwusDl7vAFZkUs+xqm/vZkXF5FxULSIy7oS2c7tvIEFTRy8zdIYvIgKEOPAbD/TgriGZIiKDQhv4uulKROS1Qhz4wUxX6tIREQFCHPj17TrDFxFJFd7Ab+smZjB9ogJfRARCHPj72ropLy0kHrNcN0VEZFwIceD3qP9eRCRFaAO/vr1H/fciIinCG/ht3czQGHwRkYNCGfi9/cm7bPUcHRGRQ0IZ+A0HkkMyZ6hLR0TkoFAGfv3Bma4U+CIig0IZ+IOPVdAoHRGRQ0IZ+PXtmstWRORwoQz8fW3dxGOmuWxFRFKEMvDr23oon1hITHfZiogclHHgm9lNZrYpmOLwQTObM0S5a8zs5eDrmkzrPZJ9uulKROR1snGGf4u7n+7uK4F7gc8fXsDMpgFfAN4EnAV8wcymZqHutHTTlYjI62Uc+O7elrJYAniaYpcCD7n7fndvBh4CVmVa91D2tXUzo1Rn+CIiqTKaxHyQmX0FuBpoBS5IU2QusDtluSZYl25ba4G1APPnzz/qtrg7F5w4g8qFo/YHhIjIcWlEZ/hm9rCZbUnztQbA3W909wrgTuC6dJtIsy7dXwK4+zp3r3T3yvLy8pHuR2pb+fp7V/KOM+Yd9c+KiITZiM7w3f3iEW7vJ8B9JPvrU9UA56cszwP+MMJtiohIFmRjlM6ylMUrgBfSFHsAeJuZTQ0u1r4tWCciImMkG334N5vZiUAC2AV8HMDMKoGPu/tH3X2/md0EPBX8zJfcfX8W6hYRkREy97Rd6eNCZWWlV1VV5boZIiLHDTN72t0r070XyjttRUTk9RT4IiIRocAXEYkIBb6ISESM64u2ZtZAcuTPsSgDGrPYnONBFPcZornfUdxniOZ+H+0+L3D3tHetjuvAz4SZVQ11pTqsorjPEM39juI+QzT3O5v7rC4dEZGIUOCLiEREmAN/Xa4bkANR3GeI5n5HcZ8hmvudtX0ObR++iIi8VpjP8EVEJIUCX0QkIkIX+Ga2ysxeNLNqM7s+1+0ZLWZWYWa/N7NtZrbVzD4drJ9mZg8Fk8U/NJpzB+eKmcXN7FkzuzdYXmRmfw32+WdmVpDrNmabmU0xs1+a2QvBMT8n7MfazP4h+Le9xcx+amZFYTzWZnaHmdWb2ZaUdWmPrSV9K8i3TWb2hqOpK1SBb2Zx4FbgMmA5cJWZLc9tq0ZNP/CP7n4ycDbwiWBfrwcecfdlwCPBcth8GtiWsvxV4BvBPjcDH8lJq0bXN4H73f0kYAXJ/Q/tsTazucCngEp3PxWIA1cSzmP9A14/x/dQx/YyYFnwtRb4ztFUFKrAB84Cqt19h7v3AncBa3LcplHh7nvd/ZngdTvJAJhLcn9/GBT7IfD23LRwdJjZPOC/AbcHywZcCPwyKBLGfZ4EvBX4HoC797p7CyE/1iTn6yg2szxgArCXEB5rd/8TcPj8IEMd2zXAjzzpCWCKmc0eaV1hC/wRT5YeJma2EDgD+Csw0933QvJDAZiRu5aNiv8D/DPJCXcApgMt7t4fLIfxmC8GGoDvB11Zt5tZCSE+1u6+B/ga8CrJoG8Fnib8x3rQUMc2o4wLW+CPeLL0sDCzicB/AX/v7m25bs9oMrPLgXp3fzp1dZqiYTvmecAbgO+4+xlAByHqvkkn6LNeAywC5gAlJLszDhe2Yz2cjP69hy3wa4CKlOV5QG2O2jLqzCyfZNjf6e6/ClbvG/wTL/hen6v2jYLzgCvMbCfJ7roLSZ7xTwn+7IdwHvMaoMbd/xos/5LkB0CYj/XFwCvu3uDufcCvgDGN488AAAElSURBVHMJ/7EeNNSxzSjjwhb4TwHLgiv5BSQv8qzPcZtGRdB3/T1gm7t/PeWt9cA1wetrgHvGum2jxd1vcPd57r6Q5LH9nbu/H/g98K6gWKj2GcDd64DdwdzRABcBzxPiY02yK+dsM5sQ/Fsf3OdQH+sUQx3b9cDVwWids4HWwa6fEXH3UH0Bq4GXgO3Ajbluzyju55tJ/im3CdgYfK0m2af9CPBy8H1arts6Svt/PnBv8Hox8CRQDfwCKMx1+0Zhf1cCVcHx/jUwNezHGvgi8AKwBfgxUBjGYw38lOR1ij6SZ/AfGerYkuzSuTXIt80kRzGNuC49WkFEJCLC1qUjIiJDUOCLiESEAl9EJCIU+CIiEaHAFxGJCAW+iEhEKPBFRCLi/wNGHA6c4s0seQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot([compute_accelerations(t, state)[0] 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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment