From a294a51807121a2dc194bd2db6e2f53447ef1e0b Mon Sep 17 00:00:00 2001 From: Roger Labbe Date: Sat, 31 Jan 2015 08:03:39 -0800 Subject: [PATCH] Started section on Runge Kutta. Mostly a placeholder for now. There is also some experimental code for quaternions, which I will be needing to introduce later. --- 07_Kalman_Filter_Math.ipynb | 185 +++++++++++++++++++++++++++++++++++- experiements/RungeKutta.py | 42 +++++++- experiements/__init__.py | 0 experiements/quaternion.py | 76 +++++++++++++++ 4 files changed, 296 insertions(+), 7 deletions(-) create mode 100644 experiements/__init__.py create mode 100644 experiements/quaternion.py diff --git a/07_Kalman_Filter_Math.ipynb b/07_Kalman_Filter_Math.ipynb index f5fde71..2fb5132 100644 --- a/07_Kalman_Filter_Math.ipynb +++ b/07_Kalman_Filter_Math.ipynb @@ -1,7 +1,7 @@ { "metadata": { "name": "", - "signature": "sha256:27a5ed2ddd67c228610a6ba471f1ec1e67b9a62114cc2be9574fab84b11cb41b" + "signature": "sha256:5f83a9ad789c1e39b3baf7784f5a29653a668d597037e30e671943b708366e53" }, "nbformat": 3, "nbformat_minor": 0, @@ -259,7 +259,7 @@ "output_type": "pyout", "prompt_number": 1, "text": [ - "" + "" ] } ], @@ -737,7 +737,7 @@ "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAGaCAYAAACSU9UtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4TGf/BvD7TPZ1JJHEmgWR1hoEEVuQqCUE9aqqkiiq\n6Gvti7e/WspbLdVWaemiiaoilqJCEBIapCR2QZGEIGMJCYns8/z+0EyNmUwSSSSZ3J/rysU85znn\nfM85M5O5c55zRhJCCBAREREREZWBrLILICIiIiKi6o/BgoiIiIiIyozBgoiIiIiIyozBgoiIiIiI\nyozBgoiIiIiIyozBgoiIiIiIyozBgohKxcXFBa6urpVdRpWgbV+EhIRAJpNh7dq1lVRV9RYVFQWZ\nTIYFCxaotfv4+EAmq5m/subPnw+ZTIbDhw9XdimlVpOPG1FNxFc7kR6RyWRF/hJPSkqCu7s7ZDIZ\npkyZgrJ8hY0kSS88b3Xg4uKi2pfafqZPnw7g6X4oal883y6TyapMIAsMDNTYJgsLC7zyyiuYMmUK\nUlJSKrtEjf2na1+XRGHgez6w6KPCD/OHDh0qsk9hWHkZ+6M83i+q0uuHiIpmWNkFEFH50vZL/MyZ\nM+jbty/u3LmDxYsXY9asWZVQWfUzdepU1KpVS6Pdy8sLAHDw4MFSLa+qBbJBgwbBw8MDAHD37l2E\nh4djxYoV2LhxI/7880+4uLhUboHP+Pnnn5GVlVXm5VS1Y1BRShrEqtP+qE61EtVUDBZEei4yMhKD\nBg1CdnY2QkJC8Pbbb1d2SdWCJEmYOnUqnJyciuxT3f+COmjQIIwaNUr1OD8/H6+99hoiIyOxaNEi\n/Pjjj5VYnbqGDRuWy3LKcqauOqkp20lEVQuHQhHpsdDQUPTp0wdCCPz+++8aoSIvLw8rV65Ev379\n4OzsDFNTU9ja2sLX1xdhYWElXs+zw0xiY2PRp08fyOVy2Nra4l//+hdu3boFALh8+TKGDh2K2rVr\nw9zcHD179sS5c+c0lvfXX39h9uzZ8PT0hL29PUxNTeHi4oJx48YhOTlZo3/huPygoCAkJSVh+PDh\nqF27NszMzNC+fftSbUtplOR6k8LagKfD0Z4dfhQUFKTW9+rVqxg7dqzqWDg4OGDIkCE4deqUxnIL\nh7KsXbsWu3fvRrdu3WBtbQ1bW9sX3h5DQ0NMmDABAHD8+HG17ZTJZMjNzcX8+fPh5uYGExMTTJs2\n7YVqB4A7d+7gnXfegaOjI8zNzdGmTRv8/PPPRdama6z+/v37MXDgQDg6OsLU1BQNGzaEv78/du3a\nBeDp0K8xY8YAABYsWKB2DEpy3cL27dsxcuRING3aFJaWlrC0tES7du2wfPlyKJVKjf6FQ80OHTqE\nLVu2oEOHDrCwsICdnR3efPNN3L59W+t64uLi0KdPH1hZWUEul8PPzw8xMTHF1leenj3W//3vf+Hi\n4gJTU1O4ublh0aJFyMvL0zrfxo0b0a5dO5ibm8PR0RGjRo0qcjtL875TUa8fIqoYPGNBpKdWrFiB\nqVOnwt7eHmFhYWjXrp1Gn9TUVEydOhWdO3fGa6+9Bnt7e9y+fRu///47BgwYgNWrV2P8+PElXueJ\nEyewZMkS+Pr64t1330VMTAy2bt2Kc+fOYfPmzejatSvatWuHMWPG4MKFC9izZw98fX2RkJAACwsL\n1XK2bduG7777Dj179kSXLl1gbGyM8+fP46effsLvv/+OuLg41K9fX2P9169fR8eOHdG4cWOMHj0a\nqamp2LRpEwICAhAREQEfH59S7cOS/NW3uOEZrq6umDdvHhYsWAC5XK72YbxwGBLwdFhVQEAAcnNz\n4e/vDzc3N9y8eRPbtm3Dnj17sGPHDvTu3Vtj+Zs3b8bevXvh7++PSZMm4c6dO6XYQk2FH5S1bVfh\nh7S+ffuidu3aqlBV2trv378Pb29vJCYmonPnzujWrRtu376N9957D76+vkXWpq2mefPmYeHChbC0\ntMSgQYPg5OSElJQUxMTE4KeffoK/vz8GDx6M9PR07NixAz4+PmrPA2dn52L3yZw5c2BgYIBOnTqh\nfv36SE9Px4EDBzBt2jQcP34c69ev1zrft99+i507dyIgIAA9evRATEwMNm3ahDNnzuD06dMwNjZW\n9T169Ch8fX2Rl5eHIUOGwM3NDWfOnIGPjw969uxZbI3lbejQoTh58iSGDh0KIyMj/Pbbb5g7dy7i\n4uLw22+/qfX98ssvMWPGDNSqVQujRo2CjY0N9u7di86dO0Mul2ssuzTvOxX9+iGiciaISG9IkiRk\nMpmYPXu2kCRJNGnSRFy7dq3I/jk5OeLWrVsa7enp6aJFixbC1tZWZGVlqU1zdnYWrq6uam3BwcFC\nkiQhSZLYtm2b2rQ+ffoISZKEXC4XX3zxhdq08ePHC0mSxPLly9Xab926JXJzczXq2rdvnzAwMBAT\nJkxQa4+MjFSt/+OPP1abtnfvXiFJkujXr18Re0GTs7OzkCRJTJ06VcybN0/t59NPPy3Rvli7dq1a\nuyRJGn0LpaWlCTs7O1G7dm1x8eJFtWkXL14UVlZWol69eiInJ0fVPm/ePCFJkjAwMBB79+4t8bYJ\nIcTo0aO11pibmyt8fHyEJEli3LhxGvujdevWIjU1tcy1jxs3TkiSJP7973+r9T916pQwNjYWkiSJ\nBQsWqE3r3r27kMlkam2Fx9bV1VXcvHlTYzufbSs8Ls8vtyQSEhI02pRKpWo/xsTEqE0rbJfL5eL8\n+fNq00aMGCEkSRKhoaFqy3J3dxcymUzj9fPNN9+oXteHDh0qUb3du3cXkiTp7F/4/Hl+fxQea3d3\nd5GWlqZqz8rKEh06dBCSJIkNGzao2hMTE4WRkZGwsbERiYmJats0bNgwVe3PepH3nfJ+/RBRxWCw\nINIjhR+uJUkSxsbG4sqVKy+8rGXLlglJksThw4fV2nV9mO7Ro4fGctatW6cKOc87fPiwkCRJjBkz\npsR1tWzZUjRq1EitrTBYuLq6CqVSqTGPk5OTsLe3L/E6Cj9cafuxsbFR61ceweLrr78WkiSJFStW\naJ0+bdo0IUmS2L17t6qt8IPhkCFDSrxdhQo/+A4aNEgVmCZOnCgaNWokJEkSjo6OIikpSWN/7Ny5\ns8y15+bmCnNzc2FlZaX2wbXQmDFjShws/P39hSRJYsuWLcVuc1mCRVHi4uKEJEli4cKFau2F+/ej\njz7SmKfwufrBBx+o2qKjo4UkSaJLly4a/ZVKpXBzcys2KDyrPILFL7/8ojHPvn37hCRJws/PT9W2\naNEiIUmS+L//+z+N/omJicLAwEDjuOlS1PtOeb9+iKhicCgUkR7q06cPwsPDMWLECOzduxc2NjZF\n9r1w4QKWLl2Kw4cPQ6FQIDs7W216UeOktWnTpo1GW506dQAArVq10phWt25dAMDNmzc1pv3yyy8I\nCQnBmTNnkJaWhoKCAtU0ExMTrev38PDQOlymYcOG+PPPP0u2EX+TJAmJiYk6L94uL0eOHAHw9O5d\n8+fP15h++fJlAMDFixfRt29ftWkdOnR44fXu2LEDO3bsAACYmprC2dkZ77//PmbPnq06NoUkSdK6\nrtLWfunSJWRlZcHb21vrMJlu3bohODi4RPXHxMRAkiSNfVLeUlNTsXTpUuzevRsJCQl48uSJ2vTC\na4ie5+npqdHWoEEDAMDDhw9VbSdPngQAdO/eXaO/JEno3Lkzrl69+sL1l5YkSVpr6dq1KwDg9OnT\nqjZdtbu4uKBhw4a4ceOGxrTyfN8py+uHiMoXgwWRnpEkCTt37sSwYcOwfft2+Pj4ICIiAvb29hp9\nY2Ji0LNnTyiVSvTq1QuDBg2CtbU1ZDIZTp06hR07diAnJ6fE69b2QdHQ0LDYac9fEDpt2jQsX74c\n9erVQ9++fVG/fn2YmZkBAIKDg7V+UAGg9dawhevRdpFtVZGamgoAWLNmTZF9JElCZmamRnthcHsR\nISEhaneFKo6jo6NGW2lrT09PL3JZutq1SUtLg7W1NczNzUs8T2mlpaWhffv2SEpKQseOHREYGAhb\nW1sYGhri4cOHWL58eZGvEW3Px8Ln/LNBuTz3CQDVxc66nvOF04q6IF7bOk1NTWFtba2qFyhZ7c+/\nXsv7facsrx8iKl8MFkR6yNDQEJs3b8bo0aPx66+/olu3boiIiNC44HnRokXIzs5GVFQUunXrpjZt\n8eLFqr9mv0x3797F119/jZYtW+Lo0aNqF3UDKPJC2eqsMHSdPHlS7YLUkqjse/uXtvbC/kVdZF6a\ni89r1aqFBw8eIDMzU+N5Ul5+/PFHJCUlYf78+Zg7d67atGPHjmH58uVlXkd57pNnl1f4gVub+/fv\nAyg6jCsUCo1b/GZnZ+PRo0eoXbu2xrru3LmDli1blqj28n7fKcvrh4jKF283S6SnDAwMsG7dOowd\nOxaXL19Gt27dkJSUpNbn6tWrsLOz0/jlDkDnt/ZWpISEBAgh0Lt3b40Pizdv3kRCQkKl1FVWkiSp\n/ZX6Wd7e3gBQolufVjWlrf3VV1+FmZkZzp49i7S0NI3ppXnederUCUII7Nmzp9i+BgYGAFDkMShK\n4RCk119/XWNaeb1GCu/YFhUVpTFNqVQiOjq6VMsr/HCta76jR48CAFq3bq0xTQihddsKj/GzQx51\n1Z6YmKj19tAv8r6jr68fIn3DYEGkxyRJwvfff48pU6YgMTERXbt2xV9//aWa7urqitTUVI3vkliz\nZg327dv3sstV1QQAf/zxh9pQjoyMDIwbN67UHwyrCjs7O9y7d09jLDkABAUFwcbGBgsXLtT6vQVC\nCERHRxf5HQKVqbS1GxoaYuTIkcjIyNA4A3Dq1Cn88ssvJV73+++/DwD44IMPtF6n8+y1D4V/Zb9+\n/XqJlw/883yMjIzUqHXx4sWlWlZRvL294e7ujqNHj2Lbtm1q01atWoVr166V6szU22+/DUNDQ/z4\n4484c+aMxvQff/wRZ8+eRdOmTdGlSxety1i4cKHakKesrCz83//9HwCofX/EW2+9BSMjI3zzzTdI\nTExUtSuVSsyePVvrcKwXed/R19cPkb7hUCiiGuDLL7+Eubk5Fi9ejO7du2P//v1o0aIFpk6dir17\n96JLly4YNmwYrK2tERsbiyNHjmDo0KHYsmXLS6/V0dERw4cPx8aNG+Hh4QE/Pz+kp6dj//79MDc3\nh4eHh9rFoyUlSvlNxKXtX5zevXvj119/RZ8+fdC1a1eYmJjAw8MD/v7+sLGxwdatWzFo0CB4e3uj\nZ8+eaNasGYyMjJCcnIw///wTycnJSEtLg5GRUbnWVVYvUvsnn3yCAwcOYOXKlTh58iS6du0KhUKB\n0NBQ9OvXDzt37tS6ruePiZ+fHz766CMsXLgQzZo1Q0BAAJycnHD37l3ExMSgSZMmqu9c8Pb2hoWF\nBTZu3AgjIyM4OTlBkiSMGjVK5wX6o0aNwtKlSzF16lRERkaiSZMmuHLlCsLCwvD6669j48aN5bIf\n16xZAz8/PwwbNgxDhgxBkyZNcPbsWURERKhuxlBSjRo1wsqVKzFx4kR06NAB/v7+cHd3R15eHmJi\nYnDkyBHY2tpi/fr1RQaWV199Fc2bN8fQoUNhYGCA7du3IzExEYMGDcLw4cNV/ZydnfHpp59ixowZ\naNu2LYYNG6b6Hov09HS0atUKZ8+eVVv2i7zv6Ovrh0jvFHfbqJUrV4pWrVoJa2trYW1tLTp16iTC\nwsKK7B8ZGSkGDhwo6tatK8zNzUWrVq3ETz/9VB53sCKiYmi7Z/yzPvnkEyFJkqhdu7aIi4sTQgix\na9cu4eXlJaysrISNjY147bXXxB9//CFCQkKETCbTuG2qi4uLxm0fC/tqu5Vn4e01g4KCNKYlJiZq\nvU3tkydPxIcffiiaNGkiTE1NhZOTk5g8ebJITU0VPj4+Gtuoax1CCK3z6OLi4iJkMpm4fv16sf2K\n2hfP77d79+6JUaNGibp166puwfl8vTdu3BBTpkwR7u7uwszMTFhZWQl3d3fx5ptvio0bN6rdSnf+\n/Pla11MSgYGBpZq3cH/oUprahRBCoVCIMWPGCHt7e2FmZibatGkj1q5dK6KiorTeBlXXMQwPDxf9\n+vUTdnZ2wtjYWDRs2FAMGDBA4/ai+/fvF126dBFWVlal+m6I+Ph4MXDgQOHg4CAsLCyEp6enWLNm\njUhKStL6vCvcv9qWXfic1/ZcjYuLE3369BFWVlbCyspK+Pn5iZiYGNWxLuntZgvFxMSIESNGCCcn\nJ2FiYiLMzc1Fs2bNxLRp07R+74cQT283K5PJRE5OjpgzZ45wcXERJiYmonHjxmLhwoUiLy9P63wb\nNmwQbdu2FaampsLBwUG8/fbbIiUlpcjjVtr3nfJ+/RBRxZCE0P1nuZ07d8LExARubm5QKpUICQnB\nkiVLcOLECa1jMxcvXoysrCz07dsXdevWRXh4ON5//338/PPPePPNNyssIBEREVHZuLi4IDk5udoO\nOSSiylVssNDGzs4On376KcaNG1ei/m+88QYKCgoqZVgFERERlQyDBRGVRaku3i4oKMDGjRuRnZ2t\n9W4ORUlPT4etrW2piyMiIiIiouqhRBdvnzt3Dp06dUJOTg7MzMwQGhoKd3f3Eq1g165dOHjwoOrW\nds969o4TREREVLkKBzHw9zMRFUfbF9+WaChUXl4ekpOTkZ6ejs2bN2PFihWIjIyEp6enzvmOHDmC\nfv36YcmSJXj33Xc1pvONi4iIiIio+nnhYPE8Pz8/NGjQAMHBwUX2iY6ORv/+/bFw4UL8+9//1tqH\nwYKIiIiIqPrRFixe6HssCgoKtH7pTaHDhw/D398fH3/8cZGhoiTFUfUWGxtb7Fkt0l88/jUbjz/x\nOVCz8fjrp+JOChQbLGbPng1/f380aNAAjx8/xq+//opDhw6pvqxnzpw5OHHiBCIiIgAAUVFR6N+/\nPyZPnow333wTCoUCAGBgYAB7e/uybg8REREREVVBxQaLO3fuYOTIkVAoFJDL5WjdujXCw8Ph5+cH\nAFAoFEhISFD1X7t2LbKzs7F06VIsXbpU1e7i4qLWj4iIiIiI9EexwULXdRTapgcHBxc7DxERERER\n6ZdSfY8FERERERGRNgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERU\nZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwW\nRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWRERE\nRERUZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERU\nZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwW\nRERERERUZgwWRERERERUZgwWRERERERUZgwWRERERERUZgwWREREVCP4z/kVQZ/tUD32mRqCf3+9\np0zLDPx0Owb8d0NZSyPSC4aVXQARERHRyyBJEqRnHm9fNBxGBiX7G2vU6ST0nL4W97f/B7bWZqr2\nFf/uCyHKuVCiaorBgoiIiKqN3LwCGBsZlMuyalmalnoe8VyKsDI3KZdaiPQBgwURERFVGp+pIXjV\nuTaMDQ2wbv9ZAMDYfm3x2bu+kCQJLsO/QlBfD1xXpOO36Evo7dkYm+YNxdHzyZjzwwHE/nUbNpam\nGOjtjs/e9VV90H+SnYeJX4Vh6+GLsDA1wpTXOwIAxHPrbtnIASv+3Q/A09AyPyQKvx44B8WDDNS3\nt8bU1ztioLc7ek5fCwCwH7wUABD4mgd+mhWAwE+3I/VRFn7/5E0AQE5uPmZ9H4GNB88jPTMHHk3q\n4PMJfujc0gnAP2c+Ij4fhTk/HMD5pLto5myP72f4o41b3Qrf30QViddYEBERUaVaH3EOABDzzVh8\nN90f3++Kw1dbYlTTv9gcg2Yu9oj7bjw+GdsT5xLu4LX//IJBXdxx9scJ2PbxGzh9TYExS3aq5pm5\nah8i4hKw7eNhOLBsFE5dUeDw2etqQ6GeDo36p2X0p9uxbv9ZfDnpNVz6eTLWzh4EWyszNHSQY+uC\nYQCA+JBJUGydieXv9/lnGc8s9D/f7Udo1AUEzwrA6R/eRUtXB/SZtR6KBxlq2/zfHw9gybu+OPnd\neNhZm+Gt/20rr91JVGl0BotvvvkGrVu3hlwuh1wuh7e3N3bv3q1zgefOnUP37t1hbm6OBg0aYOHC\nheVaMBEREemXenZWWP5+XzRtaId/+TTHB8O98cXmf4KFT2sXzHzDG43q2aBxfVss3XgUb/Rojmn/\n6oTG9W3R4dX6+HZqf2w9HI/76U+QkZWLn/acwtIJfvDzbIzmrg4InhUA2bMJ4DlXbqZiU+R5rPlg\nIAZ3fRUudWqhS0snvOXXCjKZBBurp9dVONSygIONherMiBBCdY1FZlYuVv8ehyXv+qFvRze4O9XG\n6un+cLSxwDfbj6utb+GYHuju4QJ3p9qYO6o7Lt24j9v3H5fzniV6uXQOhWrYsCGWLFkCNzc3KJVK\nhISEYNCgQThx4gRat26t0f/Ro0fw8/ODj48PYmNjcfHiRQQFBcHCwgLTp0+vsI0gIiKi6kmSJHg1\na6DW5vVqA3z0UyQeP8mBJEnwdFcfIhT3Vwqu3X6ATZEXVG1CCEiShGu3HsDU2BC5+QXo1LyharqF\nmTFaNnIoso5TVxSQSRJ6tHF54W25dvsh8vIL0LnFP+uVySR0at4Q8Un31fq2auSo+n9dO0sAwN20\nTNSrbfXC6yeqbDqDxcCBA9UeL1q0CKtWrcLx48e1Bov169cjOzsba9euhYmJCZo1a4ZLly7hiy++\nYLAgIiKqAZRKJZKT85CfXxtKpRIyWfGjrou7q5KFqbF6fwiM82+LaUM7afStV9sKl2/c12gvyXoq\nihACMpn62RIjw38uQJf+PpOiVPL2UlS9lfgai4KCAmzcuBHZ2dno1q2b1j7Hjh1D165dYWLyzx0S\nevfujdu3b+P69etlr5aIiIiqLKVSiX37cuDlZYTBg52wb18OlEqlznmEEPjz4k21tpiLN1G/tnWR\nd1xq61YX5xPvoVE9G40fU2NDNK5vCyNDAxy7kKyaJzMrF+cT7xZZh0eTOlAKgYMnE7VON/47CBTo\n2J7G9WxgbGiA6HM3VG0FBUoci7+JZs61i5yPSF8UGyzOnTsHS0tLmJqaYvz48QgNDYW7u7vWvgqF\nAo6OjmpthY8VCkU5lEtERERVVXJyHoKCTKBQyKBQyBAUZILk5Lxi57ud+hhTV4bj8o372HIoHp9v\nOoZpQ70AaN7eFQBmvdkZxy/dwntf7sKpKym4eusBdh37CxO+2AUAsDQzxjt922DW9xGIiEvAhcS7\nGLNkJ5TPLUsIAfH3faKaNrTDMJ/mGPv579h2+CISUx7ij7PX8cvfd6pydpRDkiTsOvYX7qVlIjMr\nV6MuCzNjvBfgiVnfR2DPn1dw8fo9vPdlGO6lZWJiQPvS7UyiaqjY282+8sorOHv2LNLT07F582YM\nHz4ckZGR8PT01Ogr6bgoqjixsbEvPC9VXTyuNRuPf83G41/z5OfXBuCk1paSkoJ797QPTQKAjIzH\n6N26DlIUCrSf8B0kSBjYoSG6uhoiNjYWebm5SE5ORmyskdp8q9/tiFXhl9F13xkolQL1bc3Ro2Ud\n1fNuRMfauHErBQEf/gpTI0O80cUZrZzkuJ96X9UnI+Mx7t2VVI+n9HaCicjCe1/sRFpmLhzkphjR\nzRWxNk9DxHg/N/xn9V6MXboT/T0bYO6w1khNTUX6k1zVMoa1s8WdO/YYuWgrMrLz4F7fGl8GtUNy\nwiUkJwCXr6VCAnDq9CnIzZ8O8br94AkkAPHx8VA+ulXWw1Bl8D1A/7i5uemcLgltfwrQwc/PDw0a\nNEBwcLDGtNGjRyM1NRW7du1StZ04cQIdO3ZEYmIinJ2d1fqnp6er/i+Xy0tTBlUDsbGxWgMo1Qw8\n/jUbj3/NVDgUKijo6RCm4OAc9O5tovM6ix7T1qKlqwO+/nffl1UmvQR8D9BPxX12L/X3WBQUFBQ5\nXrJTp074448/kJOTo2rbv38/6tevrxEqiIiISL/IZDL07m2CmJg8/PbbjWJDBaA+HImIqjedr/bZ\ns2cjOjoaSUlJOHfuHObMmYNDhw5h5MiRAIA5c+bA19dX1X/EiBEwNzdHYGAgLly4gG3btuGzzz7j\nHaGIiIj01JOsXKRlZKsey2QyODubwNDwfonuCPX8l9QRUfWl8xqLO3fuYOTIkVAoFJDL5WjdujXC\nw8Ph5+cH4OkF2QkJCar+1tbW2L9/PyZNmgRPT0/Y2tpi5syZmDZtWsVuBREREb00WTl52PPnVazd\newa/H7sMSZKQuXsOTE2MUPA4HU8O7oZkU7f4BQGI/HJ0BVdLRC+LzmCh7TqK4qa3aNEChw4dKltV\nREREVKXk5OYjLOYKQqMuYNexv5CZ/c/dngxkEgwNnp6dyNy3E+k/LYeZ32Cgm08lVUtElaHYu0IR\nERERjVq8HaFR/3zTtYmRAXLyCgAAw3u2gOHf3/MgcrIAAFKu5u1YiUi/MVgQERFRsTq8Wg93Hmag\na0sn/BJxFkmKdFW4eNuvVWWXR0RVQKnvCkVEREQ1z4xh3vj1/17HpqgLSFKkw72hHXLyCmBnbYYe\nbVwquzwiqgIYLIiIiKhYt+8/hs+0EFy5+QAeTepggHdTAMDgrq/A6O9hUERUszFYEBERkU7Ph4qI\nz9/G4TM3AAD/6t68kqsjoqqC11gQERFRkbSFCju5OYb1aIaWjRzQs61rZZdIRFUEgwURERFpVVSo\nAJ5ec0FE9CwOhSIiIiINukIFEZE2DBZERESkhqGCiF4EgwURERGpMFQQ0YtisCAiIiIADBVEVDYM\nFkRERMRQQURlxmBBRERUwzFUEFF5YLAgIiKqwRgqiKi8MFgQERHVUAwVRFSeGCyIiIhqIIYKIipv\nDBZEREQWcpx/AAAgAElEQVQ1DEMFEVUEBgsiIqIahKGCiCoKgwUREVENwVBBRBWJwYKIiKgGYKgg\noorGYEFERKTnGCqI6GVgsCAiItJjDBVE9LIwWBAREekphgoiepkYLIiIiPQQQwURvWwMFkRERHqG\noYKIKgODBRERkR5hqCCiysJgQUREpCcYKoioMjFYEBER6QGGCiKqbAwWRERE1RxDBRFVBQwWRERE\n1RhDBRFVFQwWRERE1RRDBRFVJQwWRERE1RBDBRFVNQwWRERE1QxDBRFVRQwWRERE1QhDBRFVVQwW\nRERE1QRDBRFVZQwWRERE1QBDBRFVdQwWREREVRxDRc0UGBiIAQMGVHYZRCXGYEFERFSFMVTov6io\nKMhkMjx48ECtfcWKFVi/fn2Fr3/+/Plo2bJlha+H9J9hZRdARERE2jFU1CxCCLXHVlZWlVQJ0Yvh\nGQsiIqIqiKGiZHx8fDBp0iT897//hb29PRwdHfHBBx+ofUjPzc3FrFmz0LBhQ1hYWKBDhw7Yt2+f\narqXlxc+++wz1eORI0dCJpPhzp07AIAnT57AxMQER48eLbKO+Ph49O/fH9bW1nB0dMSIESNU8wPA\nuXPn0KtXL8jlclhZWcHDwwNRUVFISkpCz549AQD29vaQyWQYM2YMAM2hUD4+Ppg4cSJmzJgBOzs7\nODg44Ouvv0Z2djYmTJiAWrVqwdnZGRs2bFCrbfbs2XjllVdgbm4OV1dXzJo1Czk5OQCAkJAQfPzx\nx7hw4QJkMhlkMhl+/vlnAEB6ejrGjx8PR0dHWFtbw8fHB3FxcaU7QFSjMFgQERFVMQwVpbN+/XoY\nGxvj2LFjWLlyJb766its2rRJNT0oKAh//PEHNmzYgAsXLmD06NEYMGAAzp49CwDo0aMHoqKiVP0P\nHToEe3t7VdvRo0dhZGSEDh06aF1/SkoKunXrhlatWuHEiRM4cOAAMjIyEBAQoOozYsQI1K9fHydO\nnMCZM2ewYMECmJqawsnJCVu3bgXwNJwoFAosX74cACBJEiRJ0thWuVyO48ePY/bs2Zg6dSoCAgLQ\nvHlznDx5EqNHj8aYMWPUQo2lpSWCg4Nx6dIlfPvtt9i4cSP+97//AQCGDx+OGTNmwN3dHQqFAgqF\nAsOGDYMQAv3790dKSgrCwsJw+vRpdOvWDT179oRCoXjBI0V6T+jwySefCE9PT2FtbS3s7e3FgAED\nxPnz53XNIoQQIiwsTHTs2FFYWVmJ2rVri4CAAPHXX39p9EtLS1P9kP45ceJEZZdAlYjHv2bj8X9x\nt+49Em4jvxbwmS88xq4W99MyK7ukUklb/5240a+diP9s7ktZX/fu3YW3t7dam5+fnxg7dqwQQoir\nV68KmUwmbty4odYnICBATJw4UQghxJ49e4SlpaUoKCgQV65cEdbW1uKjjz4S7777rhBCiA8//FD4\n+fkVWcNHH30kevXqpdb24MEDIUmS6rVgbW0t1q5dq3X+yMhIIUmSSE1NVWsfPXq08Pf317mt9vb2\nIiAgQPU4Ly9PGBsbi61btxZZ76pVq0STJk1Uj+fNmydatGih1ufAgQPC0tJSZGVlqbV7eHiIJUuW\nFLnsQnwP0E/FfXbXecbi0KFDmDx5Mo4dO4aDBw/C0NAQvr6+ePjwYZHzXL16FYMGDYKPjw9Onz6N\niIgIZGdno1+/fuUeioiIiPRJTTlTUVBQgJ9++kljyM6LkCQJrVq1UmurW7cu7t69CwA4efIkhBBo\n1qwZrKysVD+7d+9GQkICAKBLly7IycnB8ePHERUVha5du6JXr16qMxZRUVHw8fEpsoa4uDgcPnxY\nbflOTk6QJAnXrl0DAEyfPh1jx45Fr1698Mknn+Dy5cvlsq0ODg5qF14bGhrCxsZGtf0AsGXLFnTp\n0gV169aFlZUVpk+fjuTkZJ3riouLw5MnT2Bvb6+2XRcuXFDtN6Ln6bx4Ozw8XO3xunXrIJfLcfTo\nUfTv31/rPKdPn4ZSqcTixYtVp+9mzZqFXr164cGDB7C1tS2n0omIiPSHvoQK8fe1DUqlOZRKJWQy\n9b9hRkVF4bfffsPo0aPRtm3bclmnkZGR2mNJkp6pQwlJkhAbG6vRz8zMDMDToULt2rVDZGQk4uPj\n0aNHD3h5eeHGjRu4du0aYmNjsWTJkiLXL4SAv78/Pv/8c41pDg4OAIB58+bhrbfewp49e7B3714s\nWLAAq1evRlBQUJm3VVubUqkEAMTExODNN9/E/Pnz0adPH9SqVQs7duzAzJkzda5HqVTC0dER0dHR\nGtOsra1LVTPVHKW6K9SjR4+gVCphY2NTZJ/OnTvD0tISP/zwA9555x08efIEISEh6NChA0MFERGR\nFvoSKpRKJa5dLYADgE2hcnh75KB3bxPIZDJcu3YN3377LTp16oRly5YhMzMTycnJePToER49eoT0\n9HTV/x89eoTHjx+rPhxLkoTAwEC4uLiUuqY2bdpACIGUlBSdZx18fHxw8OBBXL58GVOnToWJiQk6\nduyIRYsW6by+AgDatm2L0NBQODk5wdCw6I9WTZo0wfvvv4/3338fEydOxI8//oigoCAYGxsDeHom\np7wdOXIE9evXx4cffqhqS0pKUutjbGysse527drhzp07kCQJrq6u5V4X6adSBYspU6agTZs26NSp\nU5F96tati927d2PQoEGYNGkSlEol2rRpgz179pS5WCIiIn2jL6ECAJKT87B9uwHG1wMyMyQEBZkg\nJiYPxsYP0L17d/j5+eH8+fOIj4+HlZUVrK2tYW1tDSsrK8jlcjRs2FDVZmlpCQMDg2LXKYTQuE1r\nYTsANG3aFG+99RYCAwOxbNkytGnTBg8ePEBUVBQaN26MwYMHA3gaLD7//HNYWlqqzqT4+Phg0aJF\n6NGjh87AMGnSJPzwww944403MGvWLNSuXRsJCQnYvHkzli1bBkNDQ8yYMQPDhg2Ds7Mz7ty5g+jo\naHh5eQEAnJ2dIUkSdu3aBX9/f5ibm8PCwqJE26pt25/l7u6OW7du4ddff4WXlxf27t2LjRs3qvVx\ndXXF9evXcerUKdUx8PX1RefOnREQEIAlS5aoLu4ODw+Hn58funTponO9VDOVOFhMnz4dR48eRXR0\ntMYdCp6VkJCAQYMGISgoCCNGjMCjR48wd+5cDBs2DAcPHixy3tjY2NJXT1Uej2vNxuNfs/H4F+9e\nejYmrI7BjfuZaFrPGktHtkTilXgkVnZhLyg/v7ZGW0pKCgwN72PLli2IiorC+fPn8eqrr6JTp04a\nH9azsrKQlZWldkej4mRkZODevXtqz7fU1FSkp6er2iZNmgRzc3NMmTIFd+/ehbW1NVq0aIFx48ap\n+hQOi2rZsqXqlqp16tRBQUEB3Nzcin0+r169Gt988w38/PyQk5ODOnXqwMvLC+fPnwcAXLt2DSNG\njMD9+/chl8vRtWtXjBw5UrXc8ePH4z//+Q/Gjh2L/v37Y+7cuRrboW1bs7OzcevWLbW2vLw8XL9+\nHbGxsahTpw5GjhyJyZMnIycnB15eXnjnnXewZMkS1TwuLi7o1KkTfHx88PjxY8ybNw/9+/fHokWL\nsGrVKowePRoPHz6Era0tPDw80LZt2xK9vvkeoH/c3Nx0TpdEcVEXwLRp0xAaGorIyEg0bdpUZ99Z\ns2YhIiJC7T7Ht27dQsOGDREdHQ1vb29Ve3p6uur/crm8uDKomomNjYWnp2dll0GVhMe/ZuPxL54+\nnakopFQqcWrhajgc/wk/pIyF9yejVUOhnhUbG4utW7fC3t4e06ZN0/kHS6qe+B6gn4r77F7sGYsp\nU6Zg8+bNJQoVwNNTcs+/gRQ+LhwrSUREVJPpY6gAnv6+b9zEAI+PA28MS8erWkIFAHh6esLT01Pt\nQwoRVX86bzc7adIkhISEqL6MpfCLUzIzM1V95syZA19fX9XjgQMH4uTJk1i4cCGuXLmCkydPIigo\nCE5OTmjXrl3FbQkREVE1oK+holDh2QeZ7InWUPEsuVzOsxVEekTnK37VqlXIyMhAr169UK9ePdXP\nsmXLVH0UCoXa/Yy7dOmCTZs2YceOHWjbti369u0LU1NThIeHq8YvEhER1UT6HiqIqGbTORSqJEOX\ngoODNdqGDh2KoUOHvnhVREREeoahgoj0ne5zlERERFRmDBVEVBMwWBAREVUghgoiqikYLIiIiCoI\nQwUR1SQMFkRERBWAoYKIahoGCyIionLGUEFENRGDBRERUTliqCCimorBgoiIqJwwVBBRTcZgQURE\nVA4YKoiopmOwICIiKiOGCiIiBgsiIqIyYaggInqKwYKIiOgFMVQQEf2DwYKIiOgFMFQQEaljsCAi\nIiolhgoiIk0MFkRERKXAUEFEpB2DBRERUQkxVBARFY3BgoiIqAQYKoiIdGOwICIiKgZDBRFR8Rgs\niIiIdGCoICIqGQYLIiKiIjBUEBGVHIMFERGRFgwVRESlw2BBRET0HIYKIqLSY7AgIiJ6BkMFEdGL\nYbAgIiL6G0MFEdGLY7AgIiICQwURUVkxWBARUY3HUEFEVHYMFkREVKMxVBARlQ8GCyIiqrEYKoiI\nyg+DBRER1UgMFURE5YvBgoiIahyGCiKi8sdgQURENQpDBRFRxWCwICKiGoOhgoio4jBYEBFRjcBQ\nQURUsRgsiIhI7zFUEBFVPAYLIiLSawwVREQvB4MFERHpLYYKIqKXh8GCiIj0EkNFzZR/5zaS/dsj\n9+qlIvvkXolHsn975N9VlOu6U7+Yj3sLppXrMomqE8PKLoCIiKi8MVTUXAb2dVDvl72QWclf/sol\n6eWvk6gK4RkLIiLSKwwV1ZvIyyvT/JJMBoNatpAMDMqpolIQ4uWvk6gK4RkLIiLSGwwV1c/d2eNh\n2LARZCamyDwYBkPHerCdNg9pa5Yj58JpSCYmMG3dHrXGTYeBjR0AIDfpKtK+X4bcK/GAEDCs0wC1\nxk+HaStP5N+5jZR3AuD41ToYN3kFAJAVexRpPyxDwV0FjN2awaLf62o1ZO7/HQ+/W4oGWw6r2rLP\nxuLef99DvQ0RMLCSo+BxOtK+/Qw58WegfJwGgzr1YT34bVj4DXh5O4uoimOwICIivcBQUX09idwD\ny75D4LD0RygfP8LdWeNg8dpg1Bo3DSI/H+lrv8H9hTPg+EUIAODBkg9h1NgdjpN+hiQzQF7SVUjG\nJlqXnX9PgfuLZsKy7xBY+v8LeQlXkPbDF6UftpSbCyO3ZrAaFgSZuQWyT/2JBys/gYFDHZi2bl/G\nPUCkH3QOhVq8eDHat28PuVwOBwcHDBw4EBcuXCjRgr/66iu88sorMDU1Rb169TBnzpxyKZiIiOh5\nDBXVm2Gd+qj1zhQY1XdG9okjMHJ1R63AyTBq4AJjlyawm74AuX9deHqGAk/DgqlHBxjVd4Zh3QYw\n6+QDk1daal12xu6tMHSsC5t3Z8KovjPMu/rCst/rpR62ZGBnD+shI2Hs6gZDx3qw7DMY5t498OTQ\n3jJvP5G+0HnG4tChQ5g8eTLat28PpVKJuXPnwtfXF/Hx8bCxsSlyvunTpyMsLAyff/45WrZsifT0\ndKSkpJR78URERAwV1Z2kGrIEALlXLyLnwkncHNrtuW4S8lNuwtitGawGv4UHXy9C5oEwmLZuD7PO\nPWHUwEXr0vOTE2Hsrh46jIsIIbqIggI83hyCJ3/sR8GDexB5eRB5eTBt5VnqZRHpK53BIjw8XO3x\nunXrIJfLcfToUfTv31/rPJcvX8bKlStx7tw5uLu7q9pbt25dDuUSERH9g6FCP0imZs88EjBr3xW1\n3pmi0U9WyxYAIB8xHuY+fZEdewTZJ2OQ/usPsJk8B5Z+A7UtHUAxZydkkuYZjIJ8tYePt63D4+2/\nota7M2Hk0gQyUzOkrf0GyrQHxW8gUQ1RqrtCPXr0CEqlUufZih07dqBRo0bYvXs3GjVqBFdXVwQG\nBuLevXtlLpaIiKgQQ4V+Mm78CvKuX4OBQx0Y1m2g9iMz++f4GtVrCKuBw2E//ytY9A5A5t7tWpdn\n2NAFuZfPq7XlXjqn9lgmt4HIyYbySeY/fRL+UuuTE38Gph27waJHXxi7usGgTn3k37zOW8wSPaNU\nwWLKlClo06YNOnXqVGSfhIQEXL9+HaGhofj555+xbt06XLp0CQMGDIDgbdiIiKgcMFToE6F2tsDS\nfxiUTzKQ+ukc5Fw+j/yUm08vlF7xPyiznkDk5uDht58h+1wc8u/cRs6l88iNPw0jp8Zal27Z73Xk\n30nBw++XIe9mEp5ERyAjfJtaH2P3FpBMzZC+diXybifjyZEDyAjbotbHqL4zck4fR078aeQlJyFt\n1RLk373NW8wSPaPEd4WaPn06jh49iujoaEg60rlSqUROTg7WrVuHJk2aAHg6hMrd3R2xsbFo3177\nnRNiY2NLWTpVBzyuNRuPf81WUcf/Xno2JqyOwY37mWhazxpLR7ZE4pV4JFbI2qi0LG7dhuXf/y/J\nc8DmcQby793DtWf6GoyaCsuI7cj8cCKk/DwUyG2R27gZEs4+PdNgfSMJxp/+F7KMdCjNLJDj3goZ\n7bojITYWsof3URtAfHw88tMyAADGw8Yhf+8WPA7bgrx6zsjq1h/W24Jx9uxZKGvdBACYDBqNvH3b\n8HjvduQ6N0V2l9dgvS0Ep0+dhjC3gOTeFtaXLiDv/yZDGBohu403pGbtkHNfgRt/126dmgpZViau\n870PAH8H6CM3Nzed0yVRgtMI06ZNQ2hoKCIjI9G0aVOdfefNm4fFixcjNzdX1SaEgLGxMTZu3IjX\nX//n3tHp6emq/8vllfANmVShYmNj4enJi9pqKh7/mq2ijj/PVFR96b9+j0frv0dG9/549T8LKrsc\nqiT8HaCfivvsXuxQqClTpmDTpk04ePBgsaECALp06YL8/HwkJCSo2hISElBQUABnZ+eS1k1ERKSG\noYKIqGrTGSwmTZqEkJAQrF+/HnK5HAqFAgqFApmZ/1zcNGfOHPj6+qoe+/r6om3bthgzZgxOnz6N\nU6dOYcyYMfDy8mJyJSKiF8JQQURU9ekMFqtWrUJGRgZ69eqFevXqqX6WLVum6qNQKNTOTkiShF27\ndsHBwQHdunVDnz594OTkhB07dlTcVhARkd5iqCAiqh50XrytVCqLXUBwcLBGW506dRAaGvriVRER\nEYGhgoioOinV7WaJiIheFoYKIqLqhcGCiIiqHIYKIqLqh8GCiIiqFIYKIqLqicGCiIiqDIYKIqLq\ni8GCiIiqBIYKIqLqjcGCiIgqHUMFEVH1x2BBRESViqGCiEg/MFgQEVGlYaggItIfDBZERFQpGCqI\niPQLgwUREb10DBVERPqHwYKIiF4qhgoiIv3EYEFERC8NQwURkf5isCAiopeCoYKISL8xWBARUYVj\nqCAi0n8MFkREVKEYKoiIagYGCyIiqjAMFURENQeDBRERVQiGCiKimoXBgoiIyt299GyGCiKiGsaw\nsgsgIiL9cvv+Y0xYHYMb9zMZKoiIahCesSAionJTOPyJoYKIqOZhsCAionLx7DUVTetZM1QQEdUw\nHApFRERl9vyF2ktHtmSoICKqYXjGgoiIykTb3Z9qWRhXdllERPSSMVgQEdEL4y1liYioEIMFERG9\nEIYKIiJ6FoMFERGVGkMFERE9j8GCiIhKhaGCiIi0YbAgIqISY6ggIqKiMFgQEVGJMFQQEZEuDBZE\nRFQshgoiIioOgwUREenEUEFERCXBYEFEREViqCAiopJisCAiIq0YKoiIqDQYLIiISANDBRERlRaD\nBRERqWGoICKiF8FgQUREKgwVRET0ohgsiIgIAEMFERGVDYMFERExVBARUZkxWBAR1XAMFUREVB4Y\nLIiIajCGCiIiKi86g8XixYvRvn17yOVyODg4YODAgbhw4UKJF37lyhVYWVnBysqqzIUSEVH5Yqgg\nIqLypDNYHDp0CJMnT8axY8dw8OBBGBoawtfXFw8fPix2wbm5uRg+fDi6d+8OSZLKrWAiIio7hgoi\nIipvhromhoeHqz1et24d5HI5jh49iv79++tc8KxZs+Dh4YFu3brh0KFDZa+UiIjKBUMFERFVhFJd\nY/Ho0SMolUrY2Njo7BcWFoawsDCsWLECQogyFUhEROWHoYKIiCqKzjMWz5syZQratGmDTp06Fdnn\n9u3bGD9+PLZv3w5zc/6yIiKqKhgqiIioIpU4WEyfPh1Hjx5FdHS0zmsm3n77bbz33nto3759qQqJ\njY0tVX+qHnhcazYe/6rjXno2JqyOwY37mWhazxpLR7ZE4pV4JFbgOnn8ayaLW7dh+ff/+Ryo2Xj8\n9Y+bm5vO6ZIowViladOmITQ0FJGRkWjatKnOvjKZDAYGBqrHQggolUoYGBhg1apVGDt2rGpaenq6\n6v9yuby4MqiaiY2NhaenZ2WXQZWEx7/qqIwzFTz+NVf6r9/j0frvkdG9P179z4LKLocqCd8D9FNx\nn92LPWMxZcoUbN68uUShAgDOnz+v9nj79u343//+hxMnTqBevXolqZmIiMoJhz8REdHLojNYTJo0\nCb/88gu2b98OuVwOhUIBALCysoKFhQUAYM6cOThx4gQiIiIAAM2aNVNbxvHjxyGTyTTaiYioYjFU\nEBHRy6TzrlCrVq1CRkYGevXqhXr16ql+li1bpuqjUCiQkJCgcyX8HgsiopeLoYKIiF42nWcslEpl\nsQsIDg7WOT0wMBCBgYGlKoqIiF4cQwUREVWGUn2PBRERVW0MFUREVFkYLIiI9ARDBRERVSYGCyIi\nPcBQQURElY3BgoiommOoICKiqoDBgoioGmOoICKiqoLBgoiommKoICKiqoTBgoioGmKoICKiqobB\ngoiommGoICKiqojBgoioGmGoICKiqorBgoiommCoICKiqozBgoioGmCoICKiqo7BgoioimOoICKi\n6oDBgoioCmOoICKi6oLBgoioimKoICKi6oTBgoioCmKooGpLJlP/l4hqDMPKLoCIiNQxVFB1Zt61\nN/KTk3C/hWdll0JELxn/nEBEVIUwVFB1Z1TfCXYfLEJBbcfKLoWKEBUVBZlMhgcPHlR2KaRnGCyI\niKoIhgrSB0qlEtev5yA/vzaUSmVll/PS5efnV3YJJSaEqOwSSM8wWBARVQEMFaQPlEol9u3LgZeX\nEQYPdsK+fTklDhc+Pj6YOHEiZsyYATs7Ozg4OODrr79GdnY2JkyYgFq1asHZ2RkbNmxQzXPr1i0M\nHz4ctra2sLW1hb+/P65evaqafu3aNQQEBKBu3bqwtLREu3btEBYWprbebdu2oVWrVjA3N4ednR18\nfHxw9+5dAMD8+fPRsmVLtf4hISGwsrJSPS7sExISgsaNG8PU1BRPnjxBeno6xo8fD0dHR1hbW8PH\nxwdxcXEaywkPD8crr7wCCwsLBAQE4NGjR9i0aROaNm2KWrVqITAwEDk5OWo1LFmyBE2aNIG5uTla\ntWqF9evXq6YlJSVBJpNh27Zt8PPzg4WFBZo3b46IiAjV9J49ewIA7O3tIZPJMGbMmBIdI6LiMFgQ\nEVUyhgrSF8nJeQgKMoFCIYNCIUNQkAmSk/NKPP/69eshl8tx/PhxzJ49G1OnTkVAQACaN2+OkydP\nYvTo0RgzZgzu3r2LJ0+eoEePHjA3N8fhw4cRExODunXrwtfXF1lZWQCAzMxM9O/fHxERETh79ixe\nf/11DBkyBJcvXwYAKBQKDB8+HEFBQbh06RIOHz6MUaNGlXq7ExMTsXHjRmzduhVnz56FsbEx+vfv\nj5SUFISFheH06dPo1q0bevbsCYVCoZovJycHX3zxBTZs2IADBw4gNjYWQ4YMwfr167Ft2zZs374d\nO3fuxKpVq1TzfPjhhwgODsa3336LixcvYs6cOXj33Xexe/dutZo+/PBDTJ06FWfPnkX79u0xfPhw\nZGZmwsnJCVu3bgUAxMfHQ6FQYPny5aXeZiKtRCVKS0tT/ZD+OXHiRGWXQJWIx79kbt17JNxGfi3g\nM194jF0t7qdlVnZJ5YLHv2ZKSsoWdeoUCEAIQIg6dQpEUlJ2iebt3r278Pb2Vmuzt7cXAQEBqsd5\neXnC2NhYbNmyRaxZs0a4ubmp9c/Pzxd2dnYiNDS0yPV4eXmJRYsWCSGEiIuLE5IkievXr2vtO2/e\nPNGiRQu1tuDgYGFpaanWx8jISNy9e1fVduDAAWFpaSmysrLU5vXw8BBLlixRLUeSJPHXX3+pps+c\nOVMYGBiI1NRUVVtgYKDw9/cXQgiRkZEhzMzMRHR0tNpyp0yZIvr16yeEECIxMVFIkiS+//571fRb\nt24JSZLEkSNHhBBCREZGCkmS1NZT3vgeoJ+K++zOu0IREVUSnqkgfdOwoRGCg3MQFGQCAAgOzkHD\nhiYlmleSJLRq1UqtzcHBQW0okqGhIWxsbHD37l2cP38eiYmJasOSACArKwsJCQkAnp6xWLBgAcLC\nwpCSkoK8vDxkZ2ejdevWAAAPDw/4+vqiRYsW6N27N3x9fTF06FDUrl27VNvdoEED2Nvbqx7HxcXh\nyZMnam0AkJ2draoNAExMTODm5qa2vXXq1IGtra1aW3x8PICnZxiys7Px2muvQZIkVZ+8vDy4urqq\nrevZfVm3bl0AUA3xIqooDBZERJWAoYL0kUwmQ+/eJoiJyUNKSgo6dHCCrBTfZ2FkZKT2WJIkrW1K\npRJCCHh4eGDTpk0ayyn8YD5z5kzs3bsXy5Ytg5ubG8zMzDBq1Cjk5uaq6t23bx9iYmKwb98+rFmz\nBnPmzMGhQ4fQqlUryGQyjQuc8/I0h3ZZWFioPVYqlXB0dER0dLRGX2tra9X/DQ3VP4bp2t7C5QLA\nrl274OTkpNbv+fmefVwYQmrixfT0cjFYEBG9ZAwVpM9kMhmcnU1w7959yGQuFbIOSZLQtm1bbNiw\nAXZ2dpDL5Vr7HTlyBKNHj8bgwYMBPD1jcPXqVbi7u6v18/LygpeXF+bOnYvmzZsjNDQUrVq1gr29\nPQKkBOUAABjxSURBVO7cuaPW9/Tp08XW165dO9y5cweSJGmcSSiLZs2awcTEBElJSfDx8Xnh5Rgb\nGwMACgoKyqkyoqd48TYR0UvEUEGknRBC4+zA84+f9dZbb8HR0REBAQE4fPgwEhMTcfjwYcycOVN1\nZ6imTZti27ZtOHXqFM6dO4eRI0eq3WEpJiYGixYtQmxsLG7cuIEdO3YgOTkZzZo1A/D0TlUPHjzA\nJ598gmvXrmHNmjWqC5918fX1RefOnREQEIDw8HAkJibi2LFjmDdvntazGCVlZWWFmTNnYubMmQgO\nDsbVq1dx+vRprF69Gj/88EOJl+Ps7AxJkrBr1y7cu3cPmZmZL1wT0bMYLIiIXhKGCqKiSZKkdt1A\nYVtRzMzMcPjwYTRq1Aj/+te/8OqrryIwMBBpaWmwsbEBAHzxxRdwcHBA165d0b9/f3h7e6Nr166q\nZdSqVQtHjx6Fv78/mjZtig8++ABz587FiBEjAAD/3969B0V1330c/+xKVUBYLwkoEKP2UaONGhEV\nvEStiPVSpO2o9dI8aOstahBqYknaSZQ2jkGJdlpN7FSghlq0atraxMQEdKVgEzRUpSbBYowTuzak\nBiIpprr7/JHxPCHcWZZld9+vmZ3h/PZ3zvkezvktfDgXhgwZol27dmn37t0aMWKEXn/9dT3++OO1\n6qqvbkl66aWX9PWvf13Lli3Tfffdp/nz56usrEzh4eENbl9D34MvtqWlpempp57S1q1bjXtDDh8+\nrAEDBjTr+yZJ4eHh2rhxo5544gn17t1ba9eubbQ/0FwmR2N/DnCxyspK4+uGTmPCcxUXFysqKsrd\nZcBN2P+1+VqoYP+DY8C3sf+9U1O/u3PGAgBczNdCBQDANxEsAMCFCBUAAF9BsAAAFyFUAAB8CcEC\nAFyAUAEA8DUECwBoY4QKAIAvIlgAQBsiVAAAfBXBAgDaCKECAODLCBYA0AYIFQAAX0ewAAAnESoA\nACBYAIBTCBUAAHyOYAEArUSoAADg/xEsAKAVCBUAANRGsACAFiJUAABQV5PBYvPmzRo9erQsFotC\nQkIUHx+v0tLSRuc5fvy45syZo7CwMAUGBmrEiBHKzMxss6IBwF0IFQAA1K/JYHHixAmtWbNGRUVF\nysvLk5+fn2JjY3X9+vUG5ykqKtKIESN08OBBlZaWatWqVVq+fLn27dvXpsUDQHsiVAAA0DC/pjoc\nPXq01vTevXtlsVhUWFioWbNm1TtPampqremVK1cqPz9fBw8e1IIFC5woFwDcg1ABAEDjWnyPRVVV\nlex2u3r06NGi+SorK9WzZ8+Wrg4A3I5QAQBA05o8Y/FlSUlJGjlypGJiYpo9z5EjR5SXl6fCwsKW\nrg4A3IpQAQBA85gcDoejuZ1TUlK0f/9+FRQUqF+/fs2a5y9/+YtmzpypZ555RitWrKj1XmVlpfF1\nWVlZc8sAgHbxYWWNVj53Su9XVGtQWLB+uXysugd2dndZAAC4xcCBA42vLRZLnfebfcYiOTlZ+/fv\nV35+frNDRUFBgWbNmqW0tLQ6oeLLoqKimlsKPERxcTH71Yd5+v6/WvGJFiZn6f2Kas5UtIKn7384\nj2PAt7H/vdMXTwrUp1nBIikpSQcOHFB+fr4GDRrUrBVbrVbNnj1bmzZt0iOPPNKseQCgI+DyJwAA\nWq7Jm7dXr16trKws5eTkyGKxyGazyWazqbq62uiTmpqq2NhYY/r48eOaMWOGVq1apQULFhjzfPjh\nh67ZCgBoI4QKAABap8lgsWvXLt24cUNTp05VWFiY8dq2bZvRx2azqby83JjOzs5WTU2N0tPT1adP\nH2OesWPHumYrAKANECoAAGi9Ji+FstvtTS7ky/9VOzMzk/+0DcCjECoAAHBOi/+PBQB4G0IFAADO\nI1gA8GmECgAA2gbBAoDPIlQAANB2CBYAfBKhAgCAtkWwAOBzCBUAALQ9ggUAn0KoAADANQgWAHwG\noQIAANchWADwCYQKAABci2ABwOsRKgAAcD2CBQCvRqgAAKB9ECwAeC1CBQAA7YdgAcArESoAAGhf\nBAsAXodQAQBA+yNYAPAqhAoAANyDYAHAaxAqAABwH4IFAK9AqAAAwL0IFgA8HqECAAD3I1gA8GiE\nCgAAOgaCBQCPRagAAKDjIFgA8EiECgAAOhaCBQCPQ6gAAKDjIVgA8CiECgAAOiaCBQCPQagAAKDj\nIlgA8AiECgAAOjaCBYAOj1ABAEDHR7AA0KERKgAA8AwECwAdFqECAADPQbAA0CERKgAA8CwECwAd\nDqECAADPQ7AA0KEQKgAA8EwECwAdBqECAADPRbAA0CEQKgAA8GwECwBuR6gAAMDzESwAuBWhAgAA\n70CwAOA2hAoAALwHwQKAWxAqAADwLgQLAO2OUAEAgPchWABoV4QKAAC8E8ECQLshVAAA4L0IFgDa\nBaECAADvRrAA4HKECgAAvF+jwWLz5s0aPXq0LBaLQkJCFB8fr9LS0iYXeu7cOU2aNEkBAQGKiIhQ\nWlpamxUMwLMQKgAA8A2NBosTJ05ozZo1KioqUl5envz8/BQbG6vr1683OE9VVZWmTZumPn36qLi4\nWDt27FB6eroyMjLavHgAHdu/rlcTKgAA8BF+jb159OjRWtN79+6VxWJRYWGhZs2aVe88OTk5qqmp\nUXZ2trp06aKhQ4fq7bffVkZGhlJSUtqucgAdkt1u15Ur/9WtW3fpxYK3CRUAAPiIFt1jUVVVJbvd\nrh49ejTYp6ioSBMnTlSXLl2Mtri4OF29elWXL19ufaUAOjy73a5XX72p6Oiv6Fvf6qu7Hf+j36Qm\n6Piz/0uoAADAy7UoWCQlJWnkyJGKiYlpsI/NZlNoaGittjvTNputFSUC8BRXrvxXS5Z0kc1mls1m\n1sMrgvTg4Ptk6dbV3aUBAAAXa/RSqC9KSUlRYWGhCgoKZDKZGuzX2HuNKS4ubtV86NjYr77l1q27\nJPWt1fbPf/5TH35Y4Z6C4FaMf3AM+Db2v/cZOHBgo+83K1gkJydr//79ys/PV79+/Rrt27t37zpn\nJq5du2a815CoqKjmlAIPUlxczH71MXa7XZmZN7VkyeeXQmZm3tSYMX1lNvdzb2Fod4x/cAz4Nva/\nd6qsrGz0/SYvhUpKSlJubq7y8vI0aNCgJlcYExOjkydP6ubNm0bbsWPHFB4ernvvvbcZJQPwVGaz\nWXFxXXTq1H91+PD7iovrIrOZf5cDAIAvaPQn/urVq5WVlaWcnBxZLBbZbDbZbDZVV1cbfVJTUxUb\nG2tML1y4UAEBAUpMTFRpaakOHTqkLVu28EQowEeYzWbde28X+flVECoAAPAhjf7U37Vrl27cuKGp\nU6cqLCzMeG3bts3oY7PZVF5ebkwHBwfr2LFjunr1qqKiorR27VqtX79eycnJrtsKAAAAAG7V6D0W\ndru9yQVkZmbWabv//vt14sSJ1lcFAAAAwKNwnQIAAAAApxEsAAAAADiNYAEAAADAaQQLAAAAAE4j\nWAAAAABwGsECAAAAgNMIFgAAAACcRrAAAAAA4DSCBQAAAACnESwAdEhms1mHDh1ydxnNNnnyZD3y\nyCPuLgMAALfxc3cBAFAfm82m7t27u7uMOrKysrR27Vp98skntdpffPFFfeUrX3H5+hMTE/XRRx/p\nT3/6k8vXBQBASxAsAHQon332mTp37qyQkJA2WU576YghCACA9sSlUABcZvLkyVq1apWSkpLUs2dP\n9ezZU4899pgcDofRp1+/ftq4caOWLl2qHj166Hvf+56kupdCnTt3TrGxsQoICFCvXr20ZMkSVVVV\nGe8nJibqm9/8prZs2aKIiAj17du3wboKCws1adIkBQYGKiIiQg8//HCtMxBWq1XR0dEKCgpS9+7d\nNXbsWJWWlur48eNaunSpqqurZTabZTabtWnTJmNb165dW2u70tLSlJiYqODgYPXt21f79+/X9evX\nNW/ePAUFBWnw4MHKy8sz5rHb7fr+97+vAQMGKCAgQIMGDVJ6errx/Xrqqaf0m9/8Rn/+85+N9Vut\nVknSBx98oO9+97vG93n27Nm6ePFiq/YbAACtQbAA4FI5OTmSpFOnTun555/X7t27tX379lp9MjIy\nNHToUJ0+fVpPP/10nWVUV1dr+vTpCg4O1ptvvqnDhw+rsLBQS5curdXvxIkTOn/+vF599VW9/vrr\n9dZz7tw5TZ8+XQkJCTp79qwOHTqkkpISY1m3bt3SnDlz9OCDD+rs2bN64403lJycrE6dOmn8+PHa\nvn27AgICZLPZZLPZtH79ekmSyWSSyWSqta7t27crOjpab731lubNm6fExEQtWLBA8fHx+tvf/qaJ\nEydq0aJFunnzpqTPg0VERIQOHDigt99+Wz/72c/09NNPKzMzU5L06KOPat68eZo2bZqx/piYGH36\n6aeaMmWKAgICZLVaderUKfXp00exsbH6z3/+09JdBgBAq3ApFACXCgsL044dOyRJgwYN0rvvvquM\njAwlJycbfSZPnmz8gl6f3/72t/r000+1d+9eBQYGSpJ2796tKVOmqLy8XAMGDJAk+fv7a8+ePY3e\n65Cenq758+cb6//qV7+qnTt3KjIyUhUVFTKbzaqsrNTs2bPVv39/o+47goODZTKZmnWp1je+8Q2t\nXLlSkrRx40ZlZGTovvvu0+LFiyVJP/nJT7Rnzx6VlpYqMjJSfn5+2rhxozF/3759dfr0ae3bt09L\nly5VYGCgunbtWudSsb1790qS9uzZY7Q999xzCg0N1ZEjRzR37twmawUAwFmcsQDgMiaTSdHR0bXa\noqOj9cEHH+jGjRtGn6ioqEaXc+HCBY0YMcIIFZIUExMjs9msv//970bb/fff3+QN1KdPn9YLL7yg\noKAg4zVhwgSZTCb94x//UM+ePZWYmKjp06dr9uzZevbZZ3XlypWWbrpMJpOGDx9uTAcGBiogIEDD\nhg0z2u6Eg3/9619G23PPPaeoqCiFhIQoKChI27dvb3L9p0+f1qVLl2ptU/fu3fXxxx+rvLy8xbUD\nANAanLEA4FJfvJ+iIV8MDC1dzhcvPwoICGjWcpYtW1brjMkdYWFhkj7/y/+6det09OhR/fGPf9QT\nTzyhF198UXFxcU0u/4u+HHJMJlOttju12+12SVJubq6Sk5O1bds2jRs3TsHBwfrFL36hw4cP11nO\nF9ntdj3wwAPKzc2tU0OPHj1aVDMAAK1FsADgMg6HQ3/9619rtZ06dUrh4eHq1q1bs5czdOhQZWZm\n6saNG8Z8hYWFstvtGjJkSItqioyM1Pnz543LpxoyfPhwDR8+XI899phmzpyp7OxsxcXFqXPnzrp9\n+3aL1tlcBQUFGjt2rB5++GGj7eLFi7WCROfOnXXr1q1a840aNUq/+93v1KtXL1ksFpfUBgBAU7gU\nCoBLXb16VevWrdM777yj3//+99q6dWu9Zwsas2jRIgUEBOihhx7S+fPnZbVatWLFCn3nO99pMiB8\n2YYNG/TGG29o1apVeuutt3Tx4kUdOXLEuBfi0qVL+tGPfqSioiJdvnxZ+fn5Onv2rL72ta9J+vxp\nTzU1NXrttddUUVFh3BztcDiadXamMYMHD9aZM2d09OhRlZWVKS0tTVartdZy+/fvr/Pnz+vdd99V\nRUWFbt26pUWLFik0NFRz5syR1WrVpUuXZLVatX79ep4MBQBoNwQLAC5jMpm0ePFi3b59W9HR0Vq+\nfLl+8IMfaN26dS1ajr+/v1555RVVVVVpzJgxSkhI0Pjx42vdrFzfU5nqM2zYMFmtVr333nuaPHmy\nHnjgAT3++OPq3bu3pM8vyyorK9PcuXM1ePBgJSYmavHixdqwYYMkady4cVq5cqUWLFigkJAQpaen\nt2j9jVmxYoXmzZunhQsXasyYMXr//ff1wx/+sNZyly1bpiFDhigqKkqhoaEqLCyUv7+/rFarBgwY\noLlz52rIkCFKTEzUxx9/zKVQAIB2Y3I4+yc2J1RWVhpfc/re+xQXFzd5Uy68V3FxsR599FENGzZM\nP//5z91dDtoZ4x8cA76N/e+dmvrdnTMWAFymLS4PAgAAnoFgAcBl2uLyIAAA4Bl4KhQAl8nPz3d3\nCQAAoJ1wxgIAAACA0wgWAAAAAJxGsAAAAADgNIIFAAAAAKcRLAAAAAA4jWABAAAAwGkECwAAAABO\nI1gAAAAAcBrBAgAAAIDTCBYAAAAAnEawAAAAAOA0ggUAAAAApxEsAAAAADiNYAEAAADAaQQLAAAA\nAE4jWAAAAABwWpPBwmq1Kj4+XhERETKbzcrOzm5yoS+99JKio6MVHBysu+++WwkJCSorK2uTggEA\nAAB0PE0Gi+rqag0fPlw7duyQv7+/TCZTo/0vXryohIQETZ48WSUlJXrttddUU1OjmTNntlnRAAAA\nADoWv6Y6zJgxQzNmzJAkJSYmNrnAkpIS2e12bd682QghGzZs0NSpU/Xvf/9bPXv2dK5iAAAAAB1O\nm99jMX78eHXr1k2/+tWvdPv2bX3yySfKysrSmDFjCBUAAACAlzI5HA5HczsHBQXpl7/8pR566KFG\n+xUWFiohIUHXr1+X3W7XyJEj9fLLL+vuu++u1a+ysrJ1VQMAAABwG4vFUqetzc9YlJeXKyEhQUuW\nLFFxcbGOHz+uoKAgzZs3Ty3IMAAAAAA8SJP3WLTU888/r3vuuUdbtmwx2l544QXdc889Kioq0rhx\n49p6lQAAAADcrM2DhcPhkNlc+0TInWm73V6rvb5TKAAAAAA8T7MeN1tSUmI87eny5csqKSnRlStX\nJEmpqamKjY01+sfHx+vMmTNKS0tTWVmZzpw5oyVLlqhv374aNWqU67YEAAAAgNs0GSzefPNNRUZG\nKjIyUjU1NXryyScVGRmpJ598UpJks9lUXl5u9J8wYYJyc3P1hz/8QZGRkZoxY4a6du2qo0ePyt/f\n33VbAgAAAMBtWvRUKAAAAACoT5s/FQq+Y+fOnerfv7/8/f0VFRWlgoKCBvu+9957MpvNdV6vvvpq\nO1aMtmK1WhUfH6+IiAiZzWZlZ2c3Oc+5c+c0adIkBQQEKCIiQmlpae1QKVyhpfuf8e89Nm/erNGj\nR8tisSgkJETx8fEqLS1tcj7Gv/dozTHAZ4DvIFigVXJzc7Vu3Tr9+Mc/VklJicaNG6cZM2YY9940\n5JVXXpHNZjNeU6ZMaaeK0Zaqq6s1fPhw7dixQ/7+/jKZTI32r6qq0rRp09SnTx8VFxdrx44dSk9P\nV0ZGRjtVjLbU0v1/B+Pf8504cUJr1qxRUVGR8vLy5Ofnp9jYWF2/fr3BeRj/3qU1x8AdfAb4AAfQ\nCmPGjHEsX768VtvAgQMdqamp9fa/dOmSw2QyOYqLi9ujPLSjbt26ObKzsxvts3PnTofFYnHU1NQY\nbT/96U8d4eHhri4PLtac/c/49143btxwdOrUyXHkyJEG+zD+vVtzjgE+A3wHZyzQYp999pnOnDmj\nuLi4Wu1xcXEqLCxsdN5vf/vbCg0N1YQJE3Tw4EFXlokOpKioSBMnTlSXLl2Mtri4OF29elWXL192\nY2VoT4x/71NVVSW73a4ePXo02Ifx792acwzcwWeA9yNYoMUqKip0+/ZthYaG1moPCQmRzWard56g\noCBt27ZNBw4c0Msvv6ypU6dq/vz5ysnJaY+S4WY2m63O8XJnuqFjBt6D8e+9kpKSNHLkSMXExDTY\nh/Hv3ZpzDPAZ4Dva/B/kAfXp1auXkpOTjenIyEh99NFHeuaZZ7Ro0SI3Vob20Nxr8OGdGP/eKSUl\nRYWFhSooKGh0jDP+vVdzjwE+A3wHZyzQYnfddZc6deqka9eu1Wq/du2a+vTp0+zljB49WmVlZW1d\nHjqg3r171/nL5J3jp3fv3u4oCW7G+PdsycnJys3NVV5envr169doX8a/d2rJMVAfPgO8E8ECLda5\nc2eNGjWqzmPijh07pnHjxjV7OSUlJQoLC2vr8tABxcTE6OTJk7p586bRduzYMYWHh+vee+91Y2Vw\nF8a/50pKSjJ+oRw0aFCT/Rn/3qelx0B9+AzwTgQLtEpKSoqysrL061//WhcuXFBSUpJsNptWrlwp\nSUpNTVVsbKzRPzs7W/v27dOFCxf0zjvvaOvWrdq5c6fWrl3rrk2AE6qrq1VSUqKSkhLZ7XZdvnxZ\nJSUlxuOGv7z/Fy5cqICAACUmJqq0tFSHDh3Sli1blJKS4q5NgBNauv8Z/95j9erVysrKUk5OjiwW\ni/HY0OrqaqMP49+7teYY4DPAh7j7sVTwXDt37nT069fP0aVLF0dUVJTj5MmTxnuJiYmO/v37G9PZ\n2dmOoUOHOgIDAx3BwcGO0aNHO3JyctxRNtpAfn6+w2QyOUwmk8NsNhtfL1myxOFw1N3/DofDce7c\nOceDDz7o6Nq1qyMsLMyxadMmd5SONtDS/c/49x5f3ud3Xhs3bjT6MP69W2uOAT4DfIfJ4XA43B1u\nAAAAAHg2LoUCAAAA4DSCBQAAAACnESwAAAAAOI1gAQAAAMBpBAsAAAAATiNYAAAAAHAawQIAAACA\n0wgWAAAAAJz2fxs+cVrkrHxBAAAAAElFTkSuQmCC\n", "text": [ - "" + "" ] } ], @@ -1337,6 +1337,185 @@ "> For reference, the Appendix **Symbols and Notations** lists the symbology used by the major authors in the field." ] }, + { + "cell_type": "heading", + "level": 2, + "metadata": {}, + "source": [ + "Numeric Integration of Differential Equations" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "** author's note: this is just notes to a section. If you need to know this in depth, \n", + "*Computational Physics in Python * by Dr. Eric Ayars is excellent, and available here.\n", + "http://phys.csuchico.edu/ayars/312/Handouts/comp-phys-python.pdf **\n", + "\n", + "So far in this book we have been working with systems that can be expressed with simple linear differential equations such as\n", + "\n", + "$$v = \\dot{x} = \\frac{dx}{dt}$$\n", + "\n", + "which we can integrate into a closed form solution, in this case $x(t) =vt + x_0$. This equation is then put into the system matrix $\\mathbf{F}$, which allows the Kalman filter equations to predict the system state in the future. For example, our constant velocity filters use\n", + "\n", + "$$\\mathbf{F} = \\begin{bmatrix}\n", + "1 & \\Delta t \\\\ 0 & 1\\end{bmatrix}$$.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Kalman filter predict equation is $\\mathbf{x}^- = \\mathbf{Fx} + \\mathbf{Bu}$. Hence the prediction is\n", + "\n", + "$$\\mathbf{x}^- = \\begin{bmatrix}\n", + "1 & \\Delta t \\\\ 0 & 1\\end{bmatrix}\\begin{bmatrix}\n", + "x\\\\ \\dot{x}\\end{bmatrix}\n", + "$$\n", + "\n", + "which multiplies out to \n", + "$$\\begin{aligned}x^- &= x + v\\Delta t \\\\\n", + "\\dot{x}^- &= \\dot{x}\\end{aligned}$$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This works for linear ordinary differential equations (ODEs), but does not work (well) for nonlinear equations. For example, consider trying to predict the position of a rapidly turning car. Cars turn by pivoting the front wheels, which cause the car to pivot around the rear axle. Therefore the path will be continuously varying and a linear prediction will necessarily produce an incorrect value. If the change in the system is small enough relative to $\\Delta t$ this can often produce adequate results, but that will rarely be the case with the nonlinear Kalman filters we will be studying in subsequent chapters. Another problem is that even trivial systems produce differential equations for which finding closed form solutions is difficult or impossible. \n", + "\n", + "For these reasons we need to know how to numerically integrate differential equations. This can be a vast topic, and scipy provides integration routines such as `scipy.integrate.ode`. These routines are robust, but " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "** material about Euler here**\n" + ] + }, + { + "cell_type": "heading", + "level": 3, + "metadata": {}, + "source": [ + "Runge Kutta Methods" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Runge Kutta integration is the workhorse of numerical integration. As mentioned earlier there are a vast number of methods and literature on the subject. In practice, using the Runge Kutta algorithm that I present here will solve most any problem you will face. It offers a very good balance of speed, precision, and stability, and it is used in vast amounts of scientific software. \n", + "\n", + "Let's just dive in. We start with some differential equation\n", + "\n", + "$$\n", + "\\ddot{y} = \\frac{d}{dt}\\dot{y}$$.\n", + "\n", + "We can substitute the derivative of y with a function f, like so\n", + "\n", + "$$\\ddot{y} = \\frac{d}{dt}f(y,t)$$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$t(t+\\Delta t) = y(t) + \\frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) + O(\\Delta t^4)$$\n", + "\n", + "$$\\begin{aligned}\n", + "k_1 &= f(y,t)\\Delta t \\\\\n", + "k_2 &= f(y+\\frac{1}{2}k_1, t+\\frac{1}{2}\\Delta t)\\Delta t \\\\\n", + "k_3 &= f(y+\\frac{1}{2}k_2, t+\\frac{1}{2}\\Delta t)\\Delta t \\\\\n", + "k_4 &= f(y+k_3, t+\\Delta t)\\Delta t\n", + "\\end{aligned}\n", + "$$\n" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "def runge_kutta4(y, x, dx, f):\n", + " \"\"\"computes 4th order Runge-Kutta for dy/dx.\n", + " y is the initial value for y\n", + " x is the initial value for x\n", + " dx is the difference in x (e.g. the time step)\n", + " f is a callable function (y, x) that you supply to compute dy/dx for\n", + " the specified values.\n", + " \"\"\"\n", + " \n", + " k1 = dx * f(y, x)\n", + " k2 = dx * f(y + 0.5*k1, x + 0.5*dx)\n", + " k3 = dx * f(y + 0.5*k2, x + 0.5*dx)\n", + " k4 = dx * f(y + k3, x + dx)\n", + " \n", + " return y + (k1 + 2*k2 + 2*k3 + k4) / 6." + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 15 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's use this for a simple example. Let\n", + "\n", + "$$\\dot{y} = t\\sqrt{y(t)}$$,\n", + "\n", + "with the initial values\n", + "\n", + "$$\\begin{aligned}t_0 &= 0\\\\y_0 &= y(t_0) = 1\\end{aligned}$$" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "import math\n", + "t = 0.\n", + "y = 1.\n", + "dt = .1\n", + "\n", + "ys = []\n", + "ts = []\n", + "\n", + "def func(y,t):\n", + " return t*math.sqrt(y)\n", + "\n", + "while t <= 10:\n", + " y = runge_kutta4(y, t, dt, func)\n", + " t += dt\n", + "\n", + " ys.append(y)\n", + " ts.append(t)\n", + "\n", + "exact = [(t**2+4)**2 / 16. for t in ts]\n", + "plt.plot(ts, ys)\n", + "plt.plot(ts, exact)\n", + "plt.show()\n", + "\n", + "error = np.array(exact) - np.array(ys)\n", + "print(\"max error {})" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "display_data", + "png": "iVBORw0KGgoAAAANSUhEUgAAAyIAAAGNCAYAAADtmMVMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XtcVHXi//E3V5GLqMQwCl5QMfCGXCRzW8pEt1rLbNtM\nt7LNza3MRCs2u8GW2VquEV7askzTdd3addvabuR207UL4D3EK2agMyYiiALCzPz+sOYX3zIuIoeZ\neT0fDx8PPedzZt6nDyjvznzO8XI4HA4BAAAAQBvyNjoAAAAAAM9DEQEAAADQ5igiAAAAANocRQQA\nAABAm6OIAAAAAGhzFBEAAAAAbY4iAgAAAKDNNVpE6uvr9eCDD6pPnz7q2LGj+vTpo0ceeUQ2m63B\nuKysLEVGRiowMFAjR45UYWFhg/21tbWaPn26wsPDFRwcrHHjxqm0tLR1zwYAAACAS2i0iMydO1fP\nP/+8Fi5cqF27dunZZ5/VkiVL9OSTTzrHzJs3TwsWLNCiRYuUl5cnk8mk0aNHq6qqyjkmPT1da9eu\n1Zo1a7R+/XpVVlZq7Nixstvt5+fMAAAAALRbXo09Wf3qq6/WBRdcoJdfftm5bfLkyTp27JjefPNN\nORwOde/eXffcc49mz54tSaqpqZHJZNL8+fM1depUVVRUyGQyafny5Zo4caIkqaSkRL169dI777yj\nMWPGnMdTBAAAANDeNHpF5Morr9QHH3ygXbt2SZIKCwv14Ycf6pe//KUkqbi4WFartUGZCAgIUGpq\nqjZu3ChJKigoUF1dXYMxUVFRiouLc44BAAAA4Dl8Gxtw1113qaSkRHFxcfL19VV9fb0efvhh3XHH\nHZIki8UiSYqIiGhwnMlk0qFDh5xjfHx8FBYW1mBMRESErFar888VFRXndjYAAAAADBEaGtqs8Y0W\nkZycHL388stas2aNBg4cqM2bN2vGjBnq3bu3brvttp881svLq1lhAAAAAHiGRj+a9cQTT+jBBx/U\nDTfcoIEDB+qmm27SrFmznIvVzWazJDW4svHdn7/bZzabZbPZVFZW1mCMxWJxjgEAAADgORq9IuJw\nOOTt3bCveHt767s17tHR0TKbzcrNzVVSUpKkM4vVN2zYoPnz50uSkpKS5Ofnp9zc3AaL1YuKijRi\nxIgffd/mXtqB68rPz1dycrLRMdDGmHfPw5x7HubcMzHvnuVcllY0WkSuvfZa/elPf1J0dLQGDBig\nzZs365lnntHkyZMlnfn4VXp6uubOnavY2FjFxMRozpw5CgkJ0aRJkySdKRVTpkxRRkaGTCaTunbt\nqlmzZik+Pl5paWktDg8AAADANTVaRJ555hl16tRJ06ZNk9VqVbdu3TR16lQ9+uijzjEZGRmqrq7W\ntGnTVF5eruHDhys3N1dBQUHOMdnZ2fL19dWECRNUXV2ttLQ0rVq1inUkAAAAgAdq9Dkiben7l3b4\naJbn4BKuZ2LePQ9z7nmYc8/EvHuWc/n5vdHF6gAAAADQ2igiAAAAANocRQQAAABAm6OIAAAAAGhz\nFBEAAAAAbY4iAgAAAKDNUUQAAAAAtDmKCAAAAIA2RxEBAAAA0OYoIgAAAADaHEUEAAAAQJujiAAA\nAABocxQRAAAAAG2OIgIAAACgzVFEAAAAALQ5iggAAACANkcRAQAAANDmKCIAAAAA2hxFBAAAAECb\no4gAAAAAaHMUEQAAAABtjiICAAAAoM1RRAAAAAC0OYoIAAAAgDZHEQEAAADQ5igiAAAAANocRQQA\nAABAm6OIAAAAAGhzFBEAAAAAbY4iAgAAAKDNUUQAAAAAtDmKCAAAAIA2RxEBAAAA0OYoIgAAAADa\nXKNFpHfv3vL29v7Br7Fjx0qSHA6HsrKyFBkZqcDAQI0cOVKFhYUNXqO2tlbTp09XeHi4goODNW7c\nOJWWlp6fMwIAAADQ7jVaRAoKCmSxWJy/Nm3aJC8vL02YMEGS9NRTT2nBggVatGiR8vLyZDKZNHr0\naFVVVTlfIz09XWvXrtWaNWu0fv16VVZWauzYsbLb7efvzAAAAAC0W40WkbCwMJlMJuevt956S6Gh\nobrhhhvkcDiUnZ2t2bNna/z48Ro4cKBWrFihEydOaPXq1ZKkiooKLVu2TPPnz9eoUaOUkJCglStX\natu2bVq3bt15P0EAAAAA7U+z1og4HA699NJLuummm9ShQwcVFxfLarVqzJgxzjEBAQFKTU3Vxo0b\nJZ25olJXV9dgTFRUlOLi4pxjAAAA4J4qyo5pX94mo2OgHWpWEXn//fd14MAB3X777ZIki8UiSYqI\niGgwzmQyOfdZLBb5+PgoLCyswZiIiAhZrdYWBwcAAED7t+vF59X7iSn6YtFCo6OgnfFtzuClS5cq\nJSVFgwcPbnSsl5dXi0NJUn5+/jkdD9fCfHsm5t3zMOeehzn3TN/N+8ljxzXs839JkirNkXw9uKGY\nmJgWH9vkInLkyBG98cYbWrJkiXOb2WyWJFmtVkVFRTm3W61W5z6z2SybzaaysrIGV0UsFotSU1PP\n+n7JyclNPwu4tPz8fObbAzHvnoc59zzMuWf6/rznZz2oDvY6bet/idKuv87gZDgfKioqWnxskz+a\ntXz5cgUEBGjixInObdHR0TKbzcrNzXVuq6mp0YYNGzRixAhJUlJSkvz8/BqMKSkpUVFRkXMMAAAA\n3Evprr0avPU91Xt5K/y2e4yOg3aoSVdEHA6HXnzxRd14440KDAx0bvfy8lJ6errmzp2r2NhYxcTE\naM6cOQoJCdGkSZMkSaGhoZoyZYoyMjJkMpnUtWtXzZo1S/Hx8UpLSzs/ZwUAAABDHXkpR2aHXZuH\nXKHk2JZ/fAfuq0lF5KOPPtK+ffuct+T9voyMDFVXV2vatGkqLy/X8OHDlZubq6CgIOeY7Oxs+fr6\nasKECaqurlZaWppWrVp1zutIAAAA0P7s37RVQ3avV623n3r9bprRcdBONamIjBw5Ujab7az7MzMz\nlZmZedb9/v7+ysnJUU5OTvMTAgAAwKWcWn7mZ74dKeOU0jOqkdHwVM26fS8AAADwU8qL9iju4CZV\n+QWq/+/uMDoO2jGKCAAAAFqF3W5X1EdrJUm7fn6DOl8Q1sgR8GQUEQAAALSKHe+8p/5H96q8QycN\nunWK0XHQzlFEAAAAcM5sdfXq9NpfJEnFYyYrsFOwwYnQ3lFEAAAAcM62/HOtehw/KEvgBYr/zSSj\n48AFUEQAAABwTmpOnlL3N5ZKknaMGC//gACDE8EVUEQAAABwTra9slymU0d1oGtvdUlJNjoOXARF\nBAAAAC12/GiZLvzgr5Kkmkn3yNubHy/RNHylAAAAoMV2L31OwXWntLNHggamjTQ6DlwIRQQAAAAt\nYin+SoPz/i1JCrpthsFp4GooIgAAAGiRQy88K397vbYOuFzRCfFGx4GLoYgAAACg2fZv2qr4nR/q\ntLevut9+j9Fx4IIoIgAAAGgWu92u6mXPSJK2p1wrc3QvgxPBFVFEAAAA0Cw7P/hIsSVbdcIvSBfe\nfqfRceCiKCIAAABoMltdvTquXihJ2jPqNwoN62pwIrgqiggAAACabOvatep17ICOBF2gIZNvNToO\nXBhFBAAAAE1SffKkIv/9vCTp8DW3q0PHjgYngiujiAAAAKBJtr/0oi44dUzFF/RV/PW/MjoOXBxF\nBAAAAI06WmpR3Cd/kyTV35wuHx8fgxPB1VFEAAAA0KgDz2crsL5WO/perNhLLzE6DtwARQQAAAA/\nqXjLDsVvy1W9l48u+P0so+PATVBEAAAAcFZ2u12nlj4tbzm0LfkaRfbvZ3QkuAmKCAAAAM7qy/fe\nV2zpNlX6B6v/HXcbHQduhCICAACAH1V3+rQ6rTnz8MK9o2/h4YVoVRQRAAAA/Kitf/2roipKdSik\nm4bePNnoOHAzFBEAAAD8wIny4+rz7suSpPIJd8svwN/gRHA3FBEAAAD8wM4XFiu09oR2dR+kQVdd\nYXQcuCGKCAAAABoo3b1XQz77lySpw5R75e3Nj4xofXxVAQAAoIGjf5kvP4dNWwaPUd+koUbHgZui\niAAAAMDpy3UfatD+z3XKN0C975ppdBy4MYoIAAAAJJ25XW/wqgWSpJ2jblZYN7PBieDOKCIAAACQ\nJG15ZaWijpfocIhZQ397m9Fx4OYoIgAAANDxI2WKee/M7XqPTZwh/4AAgxPB3VFEAAAAoD3PZSuk\n7qR29kjQoCvGGB0HHqBJReTw4cOaPHmyTCaTOnbsqIEDB+qTTz5pMCYrK0uRkZEKDAzUyJEjVVhY\n2GB/bW2tpk+frvDwcAUHB2vcuHEqLS1tvTMBAABAixRv2aH4zW+p3stbne64n9v1ok00+lV2/Phx\n/exnP5OXl5fefvttFRUVadGiRTKZTM4x8+bN04IFC7Ro0SLl5eXJZDJp9OjRqqqqco5JT0/X2rVr\ntWbNGq1fv16VlZUaO3as7Hb7+TkzAAAANMput6v6hafkI4e2JV2tngPjjI4ED+Hb2ICnnnpKkZGR\nWr58uXNbr169nL93OBzKzs7W7NmzNX78eEnSihUrZDKZtHr1ak2dOlUVFRVatmyZli9frlGjRkmS\nVq5cqV69emndunUaM4bLfwAAAEbY/uZbGnRouyo6hOjCO2cYHQcepNErIq+//rpSUlI0YcIERURE\nKCEhQYsXL3buLy4ultVqbVAmAgIClJqaqo0bN0qSCgoKVFdX12BMVFSU4uLinGMAAADQtmpOnlL4\na4skSfuvnKJOYV0MTgRP0mgR2b9/v5YsWaJ+/fopNzdXM2bM0AMPPOAsIxaLRZIUERHR4DiTyeTc\nZ7FY5OPjo7CwsAZjIiIiZLVaW+VEAAAA0Dzblj6viKojOtill4b+ZpLRceBhGv1olt1uV0pKip54\n4glJUnx8vPbs2aPFixdr2rRpP3msl5dXi4Pl5+e3+Fi4HubbMzHvnoc59zzMeftVdeSoUj75myRp\n7+gbVbZ1a6u9NvPuOWJiYlp8bKNFpHv37howYECDbbGxsTp48KAkyWw+88RNq9WqqKgo5xir1erc\nZzabZbPZVFZW1uCqiMViUWpq6o++b3JycjNPBa4qPz+f+fZAzLvnYc49D3Pevm2ddac62Ou0NfYy\njZ50Y6u9LvPuWSoqKlp8bKMfzfrZz36moqKiBtt2796t3r17S5Kio6NlNpuVm5vr3F9TU6MNGzZo\nxIgRkqSkpCT5+fk1GFNSUqKioiLnGAAAALSNHbnrNGj/ZzrlG6Ae0zOMjgMP1egVkZkzZ2rEiBGa\nO3eubrjhBm3evFkLFy7Uk08+KenMx6/S09M1d+5cxcbGKiYmRnPmzFFISIgmTTrzWcPQ0FBNmTJF\nGRkZMplM6tq1q2bNmqX4+HilpaWd3zMEAACA0+maGoWuXCBJKkqbrGGR3QxOBE/VaBFJTk7W66+/\nrgcffFCPP/64evXqpTlz5ujOO+90jsnIyFB1dbWmTZum8vJyDR8+XLm5uQoKCnKOyc7Olq+vryZM\nmKDq6mqlpaVp1apV57SOBAAAAM2z5cWlSjpxWCWdozT0ttuMjgMP1mgRkaSrrrpKV1111U+OyczM\nVGZm5ln3+/v7KycnRzk5Oc1LCAAAgFZh/eprDfxwlSTp1OQM+fn7G5wInqzRNSIAAABwD4cXzlOA\n7bS29f+54kb+3Og48HAUEQAAAA9Q+N+PNHjv/1Tt00GRd7NAHcajiAAAALi5uprTCn5lviSp8PKb\nZeoZ1cgRwPlHEQEAAHBzm5ctVWRFqUo7dVfClNuNjgNIoogAAAC4NUvxVxr431ckSVWT75dfAAvU\n0T5QRAAAANyYNWfumQXqF6ZqwKjLjI4DOFFEAAAA3NS2t97RoOIvdNK3o6LuecDoOEADFBEAAAA3\ndKqySuF/PfME9V1X/k7hPEEd7QxFBAAAwA3teG6hTKeOqviCvkq85Raj4wA/QBEBAABwM19tL1T8\nZ/+UXV7y+v2D8vHzNToS8AMUEQAAADdis9lUvfgJ+Tps2pI4Vn2HJRodCfhRFBEAAAA3suXvryrG\nUqhjAZ0VO32W0XGAs6KIAAAAuInjR8sU/e/nJEkl19+tkC6dDU4EnB1FBAAAwE3sffYphdaeUFFU\nvOKvG290HOAnUUQAAADcQNHHG5SwPVenvX0VevdD8vbmxzy0b3yFAgAAuLja6moFvTRXkrT9spsU\nFRtjcCKgcRQRAAAAF7d1ySJ1rzysrzv3UMLUO42OAzQJRQQAAMCFfbW9UPEb/i5JOv37h+UX4G9w\nIqBpKCIAAAAuymazqWbhY/J12LQpYaz6X5xidCSgySgiAAAALmrzK6+o35FdOtqxq2Jn3Gd0HKBZ\nKCIAAAAu6JuDJer/1lJJkmXSvQrpHGpwIqB5KCIAAAAu6NCzcxVUX60dfS9W/NVXGR0HaDaKCAAA\ngIvZ+uZbGrTvU5307ajImQ8bHQdoEYoIAACACzlRflzdVv9ZkrR77FRdENXd4ERAy1BEAAAAXMiu\nBU8qrLpceyJilXDzzUbHAVqMIgIAAOAiCv/7kRK256rW209B6X+Uj4+P0ZGAFqOIAAAAuIBTlVXq\nvGyuJGnH5beoR1x/gxMB54YiAgAA4AK+fPZpRZz8RsUX9FXC1KlGxwHOGUUEAACgndu14VMlFryh\nei8f+dydJT9/f6MjAeeMIgIAANCO1Zw8pcAXHpckbf35jYoeOsjgREDroIgAAAC0Y9sWP6vulYf1\ndeceGnrX3UbHAVoNRQQAAKCd2pe3SUM3viabvFR/Z6b8AwKMjgS0GooIAABAO1RXc1reS/4oHzm0\n5aLr1O+iJKMjAa2KIgIAANAObV6So57lB3U4xKzB98w0Og7Q6hotIllZWfL29m7wq3v37j8YExkZ\nqcDAQI0cOVKFhYUN9tfW1mr69OkKDw9XcHCwxo0bp9LS0tY9EwAAADexL2+T4tf/TXZ56eTvH1XH\noCCjIwGtrklXRGJjY2WxWJy/tm/f7tw3b948LViwQIsWLVJeXp5MJpNGjx6tqqoq55j09HStXbtW\na9as0fr161VZWamxY8fKbre3/hkBAAC4sNrqavkszpKvw64tKdfqwksuNjoScF74NmWQj4+PTCbT\nD7Y7HA5lZ2dr9uzZGj9+vCRpxYoVMplMWr16taZOnaqKigotW7ZMy5cv16hRoyRJK1euVK9evbRu\n3TqNGTOmFU8HAADAtW1dmK2k41+rtFN3DZ5xn9FxgPOmSVdE9u/fr8jISPXp00cTJ05UcXGxJKm4\nuFhWq7VBmQgICFBqaqo2btwoSSooKFBdXV2DMVFRUYqLi3OOAQAAgLT7sy+cd8mqvStLAUGBRkcC\nzptGi8jw4cO1YsUKvffee1q6dKksFotGjBihY8eOyWKxSJIiIiIaHGMymZz7LBaLfHx8FBYW1mBM\nRESErFZra50HAACAS6s+eVIdl2SduUvWiBsUM3yY0ZGA86rRj2ZdccUVzt8PGjRIF198saKjo7Vi\nxQpddNFFZz3Oy8vrnILl5+ef0/FwLcy3Z2LePQ9z7nmY86Y7+fo/dEnlYX0VGiX7JZe69H87V86O\n5omJiWnxsU1aI/J9gYGBGjhwoPbu3atrr71WkmS1WhUVFeUcY7VaZTabJUlms1k2m01lZWUNropY\nLBalpqae9X2Sk5ObGw0uKj8/n/n2QMy752HOPQ9z3nRF6zcqfkeu6r285bjncQ1PGmp0pBZj3j1L\nRUVFi49t9nNEampqtHPnTnXr1k3R0dEym83Kzc1tsH/Dhg0aMWKEJCkpKUl+fn4NxpSUlKioqMg5\nBgAAwFOdqqxSyPOPyVsObf35JPV14RICNEejV0Tuu+8+XXPNNerRo4eOHDmixx9/XNXV1Zo8ebKk\nM7fmnTt3rmJjYxUTE6M5c+YoJCREkyZNkiSFhoZqypQpysjIkMlkUteuXTVr1izFx8crLS3t/J4d\nAABAO/flM39SYpVVB8KilXDXdKPjAG2m0SJSWlqqiRMn6ujRowoPD9fFF1+szz77TD169JAkZWRk\nqLq6WtOmTVN5ebmGDx+u3NxcBX3vwTvZ2dny9fXVhAkTVF1drbS0NK1ateqc15EAAAC4su3vrVPi\n5rd02ttXPjMek1+Av9GRgDbTaBH529/+1uiLZGZmKjMz86z7/f39lZOTo5ycnOalAwAAcFPlR76R\n+eUnJEk7xkzRsCGDDE4EtK1mrxEBAADAubHb7To4L1Nda45rV/dBSpzyO6MjAW2OIgIAANDGtvz9\nNQ3a96mq/ALVNeMJ+fg1+0amgMujiAAAALQhy74D6vfPZyVJ+36dLnPvngYnAoxBEQEAAGgjtrp6\nHX/6IQXVV2tbzCUaev2vjI4EGIYiAgAA0EYKXnheMZZClXXsoug/ZMnbmx/F4Ln46gcAAGgDxZu3\nasi6lyVJR377kDpfEGZwIsBYFBEAAIDzrPrkSXk9+4j8HDZtSrpag8aMMjoSYDiKCAAAwHm2Y/6f\n1OP41yoJjdSgWX8wOg7QLlBEAAAAzqNtb72jxM3/0WlvX9lnPqmOQUFGRwLaBYoIAADAeXK05JB6\nrJgrSdpx5VRFDx1scCKg/aCIAAAAnAc2m03WJ2er0+kqFfZMUuJttxkdCWhXKCIAAADnQcELzyu2\ndJvKA0IVOfsJ+fj4GB0JaFcoIgAAAK1sX8FmxecukyQd/u3D6totwuBEQPtDEQEAAGhFJytPyD/7\nIfk6bNo0bLwG/yLN6EhAu0QRAQAAaEVFTz2u7icO66uuvTVk1v1GxwHaLYoIAABAK9n06j80dMf7\nqvXxk++9T6pDx45GRwLaLYoIAABAKyjdtVcxr86XJO0cP0M9B8YanAho3ygiAAAA56i2ulo18+5X\nYH2ttsZeqoRJE42OBLR7FBEAAIBztO3pJ9X72AEdCumm/g/8Ud7e/IgFNIbvEgAAgHOw5d9vKnHT\nmzrt7au6mU8quHOo0ZEAl0ARAQAAaKHD+w4oetWfJEk7xt6lPonxBicCXAdFBAAAoAXqak6r6k8Z\nCq47pe19L1bSrZONjgS4FIoIAABAC2x+5in1+WaPrMEm9XnwCdaFAM3EdwwAAEAzbXv7XSV9/k/V\ne/no1PS56hTWxehIgMuhiAAAADTDob371evlxyVJW6+4Xf0uSjI4EeCaKCIAAABNVHPylKqfvM+5\nLiT59tuNjgS4LIoIAABAE+14ao56lxXrcIhZfR95knUhwDnguwcAAKAJCta8qoSt76jWx0919z6t\nEJ4XApwTiggAAEAjDmzbodjX/ixJ2vmrWYoeOsjgRIDro4gAAAD8hBPHK+QzP0MBttPaPHiMkibd\naHQkwC1QRAAAAM7Cbrdr35yH1L3ysA507a1BD2QaHQlwGxQRAACAsyh48UUN3vs/VfkFKmD2fAUE\nBRodCXAbFBEAAIAfsfPj9RryzvOSpAO3PKTImL4GJwLcS7OKyJNPnrlN3fTp0xtsz8rKUmRkpAID\nAzVy5EgVFhY22F9bW6vp06crPDxcwcHBGjdunEpLS889PQAAwHlg/eprmZY8LF+HXQWX3Kj4q68y\nOhLgdppcRD777DMtXbpUQ4YMkZeXl3P7vHnztGDBAi1atEh5eXkymUwaPXq0qqqqnGPS09O1du1a\nrVmzRuvXr1dlZaXGjh0ru93eumcDAABwjmqrq1U5Z5Y611ZqZ89EJc241+hIgFtqUhGpqKjQTTfd\npJdfflldunRxbnc4HMrOztbs2bM1fvx4DRw4UCtWrNCJEye0evVq57HLli3T/PnzNWrUKCUkJGjl\nypXatm2b1q1bd37OCgAAoAXsdrt2zM1Sn2/2yhIcoR6ZT8vHz9foWIBbalIRmTp1qn7961/r0ksv\nlcPhcG4vLi6W1WrVmDFjnNsCAgKUmpqqjRs3SpIKCgpUV1fXYExUVJTi4uKcYwAAANqDTa+s1NDt\nuarx8dfp++YrNKyr0ZEAt9VoxV+6dKn279/vvMLx/Y9lWSwWSVJERESDY0wmkw4dOuQc4+Pjo7Cw\nsAZjIiIiZLVazy09AABAK9n96Rca9O+FZ34/MUMJPLQQOK9+sojs2rVLDz30kDZs2CAfHx9JZz6O\n9f2rImfz/cICAADQnh0ttajLs3+Qn8OmTSnjNez6XxkdCXB7P1lEPv30Ux09elQDBw50brPZbFq/\nfr2ef/557dixQ5JktVoVFRXlHGO1WmU2myVJZrNZNptNZWVlDa6KWCwWpaamnvW98/PzW3ZGcEnM\nt2di3j0Pc+55XGHO60/XqdPL2bqw5rgKI+Jkv3yMS+Ruz/jv5zliYmJafOxPFpHx48crJSXF+WeH\nw6Hf/va36t+/vx588EHFxMTIbDYrNzdXSUlJkqSamhpt2LBB8+fPlyQlJSXJz89Pubm5mjhxoiSp\npKRERUVFGjFixFnfOzk5ucUnBdeSn5/PfHsg5t3zMOeexxXm3G63a3PmbF34zW59EximqD9mq4vZ\nZHQsl+YK847WU1FR0eJjf7KIhIaGKjQ0tMG2wMBAdenSRQMGDJB05ta8c+fOVWxsrGJiYjRnzhyF\nhIRo0qRJzteYMmWKMjIyZDKZ1LVrV82aNUvx8fFKS0trcXAAAIBzVbBsmRK/XZx+6r4/y0wJAdpM\ns+9H5+Xl1WD9R0ZGhqqrqzVt2jSVl5dr+PDhys3NVVBQkHNMdna2fH19NWHCBFVXVystLU2rVq1i\nHQkAADDMl+s+VPx/lkiS9tz0oIYmxhucCPAszS4iH3744Q+2ZWZmKjMz86zH+Pv7KycnRzk5Oc19\nOwAAgFZXunuvIp9/RD5yqCD1N0oZP87oSIDHafKT1QEAANzBieMVqn8iXSF1J7Wjz3AlzZhpdCTA\nI1FEAACAx7DV1as4635FVZTq68491S/rKecjCgC0LYoIAADwGAXPzNfAA3mq9A9Wh0eeVVCnEKMj\nAR6LIgIAADzCpr/9XUkb/656L29Zps1Vt769jY4EeDSKCAAAcHtF6/+nAa8+LUnafu0MxV36c4MT\nAaCIAAAAt1a6e68iFj4gP4dNm1LGK3nyLUZHAiCKCAAAcGMVZcdkmzNDnU5X6cveKUq8/0GjIwH4\nFkUEAACn6ehjAAAgAElEQVS4pbqa0yp9dKYiKw/pq6691fex+fLxa/Yj1ACcJxQRAADgdux2u7bO\nfVQXlm7TsYDOCs5ayB2ygHaGIgIAANxOwQvPK2Hbe6rx8VflvX+WqWeU0ZEA/B8UEQAA4Fa2vvmW\nhr67VJK055aH1XdYosGJAPwYiggAAHAbuz/9QjHL/yhvOVQw+ncaOu5qoyMBOAuKCAAAcAulu/bq\nguz71MFWp01Dr1LynXcaHQnAT6CIAAAAl1duOSI9frdCa0/oy97DlPBQlry9+TEHaM/4DgUAAC6t\n+sRJffPI3TJXWbX/gr7q99gC+fr5GR0LQCMoIgAAwGXZ6uq1O/Ne9f1mj6zBJnV9fJECOwUbHQtA\nE1BEAACAS7Lb7dr05GMatP9zVfoHy/bQQoV1MxsdC0ATUUQAAIBLyn/+L0rc9KZqvf105J6n1COu\nv9GRADQDRQQAALicgjWvKum9M88K2X3zQ7rwkosNTgSguSgiAADApex4930N+vs8SdKmsXdr6Phx\nBicC0BIUEQAA4DJ2f/aFopc+LF+HXQWXTNSw300xOhKAFqKIAAAAl3DwyyKFL7hXAbbT2jzkF0qe\ndZ/RkQCcA4oIAABo944cLJH/E3er0+kq7ehzkYY+8jgPLARcHN/BAACgXasoO6bqR+9U+Kky7TEP\n0IWP88BCwB1QRAAAQLt1qrJKltl3Kep4iQ526anuTyxSQFCg0bEAtAKKCAAAaJdqq6tV/ODd6ndk\nl6xB4Qp5/Dl1CutidCwArYQiAgAA2p36ujoVPTRLsSVbdSygsxyZS3RBVHejYwFoRRQRAADQrths\nNm3NnK1B+z9TpX+wTj64SJH9+xkdC0Aro4gAAIB2w263a9MTf9TQwv/qlG8Hld2brd5DBhodC8B5\nQBEBAADtRv6Cp5W46U3VevupZNrT6ndRktGRAJwnFBEAANAufLFksZI2rFG9l7f2/u4xxY38udGR\nAJxHFBEAAGC4vJeWKSn3RdnlpcLfPKghV11hdCQA5xlFBAAAGCr/lZVKfHOhJGnrdTOVcP2vDE4E\noC1QRAAAgGEKVq9RwtoFkqRNV09X8i03G5wIQFtptIgsXrxY8fHxCg0NVWhoqEaMGKG33367wZis\nrCxFRkYqMDBQI0eOVGFhYYP9tbW1mj59usLDwxUcHKxx48aptLS0dc8EAAC4lM2v/VNDXn1KklRw\n5R0aNuU2gxMBaEuNFpEePXroqaee0ubNm1VQUKDLL79c1157rbZu3SpJmjdvnhYsWKBFixYpLy9P\nJpNJo0ePVlVVlfM10tPTtXbtWq1Zs0br169XZWWlxo4dK7vdfv7ODAAAtFtbXn9DA//6hLzlUEHa\nbUr5/e+NjgSgjTVaRK655hr94he/UJ8+fdSvXz/NmTNHISEh+uKLL+RwOJSdna3Zs2dr/PjxGjhw\noFasWKETJ05o9erVkqSKigotW7ZM8+fP16hRo5SQkKCVK1dq27ZtWrdu3Xk/QQAA0L5s+8/bilvx\nR/nIoYJLb1bK3dONjgTAAM1aI2Kz2bRmzRrV1NQoNTVVxcXFslqtGjNmjHNMQECAUlNTtXHjRklS\nQUGB6urqGoyJiopSXFyccwwAAPAM2995T/1felS+DrsKfnajkmekGx0JgEF8mzJo+/btuvjii1Vb\nW6uOHTvq1Vdf1YUXXugsEhEREQ3Gm0wmHTp0SJJksVjk4+OjsLCwBmMiIiJktVpb4xwAAIAL2Pb2\nu+q/9GH5OWwquOg6Jd97v7y9uW8O4KmaVERiY2O1bds2VVRU6LXXXtONN96oDz/88CeP8fLyOqdg\n+fn553Q8XAvz7ZmYd8/DnHue7+a8bMs2XfrWEvk5bFo/6AoFjbpCmzZtMjgdzhe+1z1HTExMi49t\nUhHx8/NTnz59JEkJCQnKy8vT4sWL9eijj0qSrFaroqKinOOtVqvMZrMkyWw2y2azqaysrMFVEYvF\notTU1LO+Z3JycvPPBi4pPz+f+fZAzLvnYc49z3dzvvXNt3TZW0vk67Cp4OLrlXr/bK6EuDG+1z1L\nRUVFi49t0d8CNptNdrtd0dHRMpvNys3Nde6rqanRhg0bNGLECElSUlKS/Pz8GowpKSlRUVGRcwwA\nAHBPW/79pmKXPXqmhIy4QcmUEADfavSKyAMPPKCxY8cqKirKeTesjz/+WO+++66kM7fmnTt3rmJj\nYxUTE+O8q9akSZMkSaGhoZoyZYoyMjJkMpnUtWtXzZo1S/Hx8UpLSzu/ZwcAAAxTlr9Zl733lzML\n0y+5UcmzWBMC4P9rtIhYrVbddNNNslgsCg0NVXx8vN59912NHj1akpSRkaHq6mpNmzZN5eXlGj58\nuHJzcxUUFOR8jezsbPn6+mrChAmqrq5WWlqaVq1adc7rSAAAQPu0+Z//+v8l5OeTlDzzXkoIgAa8\nHA6Hw+gQ3/n+Z8xCQ0MNTIK2xGdJPRPz7nmYc89RsPpvGvLq02ceVph6k5LTZ1JCPAjf657lXH5+\nb9JidQAAgKbIe2mZEt9cKEn6cNj1Spt1r8GJALRXFBEAANAqvli8WEnvvyhJ2vTLu9Q5IcHgRADa\nM66TAgCAc2K32/XFgvlKev9F2eWlzdfdq2G33250LADtHFdEAABAi9ntduXPe0JJn69VvZe3vpw4\nW8k3XG90LAAugCICAABaxFZXr81PPKqkLe+ozstHu277oxKv/qXRsQC4CIoIAABottM1Nfry0fuV\nsHuDan38tP/3cxU/hueDAWg6iggAAGiWU5VV2v/QPRry9WZV+QXq8D1Pa9DPRxgdC4CLoYgAAIAm\nO360TEdm36m4b/aoPCBUJx5YqNihg42OBcAFUUQAAECTfHOwRKcevVN9j5fIGmyS49HFiu7fz+hY\nAFwURQQAADSqpGiPfB6/S1Enj+pgl54KnvMXhUd2MzoWABdGEQEAAD9pX94mhc6fqS61ldobESvz\n3MUKDetqdCwALo4iAgAAzmrHu++r99JH1NFWq8KeSer3xLPqGBJkdCwAboAiAgAAflT+qr9q8D8X\nyNdh15aBaRqS+YT8/P2NjgXATVBEAABAA3a7Xfk52Ur6aKUkqSD1N0pOnyVvb2+DkwFwJxQRAADg\nVF9Xpy2PP6Kkbe/JJi9tu26mUm652ehYANwQRQQAAEg686DCvY/OVMKBfNX4+GvflD8q+aorjI4F\nwE1RRAAAgI6WHFJF5nQNLNuvig4hKpv5Zw0ZPszoWADcGEUEAAAPV7xluwKfSlfvU8d0OMQsr4cX\nKuZCHlQI4PyiiAAA4MF2vPu+er34iALra7XHPEDmx3LU2RRmdCwAHoAiAgCAh8p7aZni31wkHzm0\nNW6kBmXOlX9AgNGxAHgIiggAAB7GVlevTfPmKDH/35KkgstuUfI9M7g9L4A2RREBAMCDnDheoeKs\n+5R4IF+nvX2188Y/KOWG642OBcADUUQAAPAQpXv2qf6JdA08XqKKDiH6Zvo8JV5ysdGxAHgoiggA\nAB6g8MNP1O25h9TpdJUOdumpgIef1YV9exsdC4AHo4gAAODG7Ha7Cpav0JA3F8nXYdeO6BT1++N8\nBXUKMToaAA9HEQEAwE3V1ZzW1j9lKXHLO5KkgktuVNLM++Tj42NwMgCgiAAA4JaOHbbK+ti9Sjj8\npWp9/LRr4h+Ucv2vjI4FAE4UEQAA3MzezwsU8myG+p86pqMdu6pq1tNKGJZodCwAaIAiAgCAGylY\n9VcNWJstf3u99pjjdMEjC9Q30mx0LAD4AYoIAABu4HRNjbY9+ZgStp5ZD7Ip6WoNvf9h+QX4G5wM\nAH4cRQQAABf3zcESlc+5VwlHdqvW2087b7hPw268wehYAPCTKCIAALiwnR+vV/hzj6hvTYWOBF2g\nU7OeVlLSUKNjAUCjKCIAALggm82mgueWaOi6l+Ujh4oih6h75nx1M4UbHQ0AmoQiAgCAizl+tEwH\n5zygpAP5kqSCSyYqccZM+fr5GZwMAJrOu7EBTz75pIYNG6bQ0FCZTCZdc801+vLLL38wLisrS5GR\nkQoMDNTIkSNVWFjYYH9tba2mT5+u8PBwBQcHa9y4cSotLW29MwEAwAPs/bxA1TNu1MAD+ar0D1bh\nnU8r5b4MSggAl9NoEfn44491991369NPP9UHH3wgX19fpaWlqby83Dlm3rx5WrBggRYtWqS8vDyZ\nTCaNHj1aVVVVzjHp6elau3at1qxZo/Xr16uyslJjx46V3W4/P2cGAIAbsdvtyntpmaL+9HuZTh7V\nvvAY1T71Vw3+RZrR0QCgRRr9aNa7777b4M8rV65UaGioNm7cqF/+8pdyOBzKzs7W7NmzNX78eEnS\nihUrZDKZtHr1ak2dOlUVFRVatmyZli9frlGjRjlfp1evXlq3bp3GjBlzHk4NAAD3UFlWrv1PPqLE\nvf+TJG1KHqeh9z3IrXkBuLRGr4j8X5WVlbLb7erSpYskqbi4WFartUGZCAgIUGpqqjZu3ChJKigo\nUF1dXYMxUVFRiouLc44BAAA/tPuzL1R1zwQN3vs/nfTtqO2TszTs4SxKCACX1+zF6jNmzFBCQoIu\nvvhiSZLFYpEkRURENBhnMpl06NAh5xgfHx+FhYU1GBMRESGr1dqi4AAAuLMzd8V6TvH/fVm+Drv2\nh8co6A9/0tB+fYyOBgCtollFZNasWdq4caM2bNggLy+vRsc3ZczZ5Ofnt/hYuB7m2zMx756HOW+a\n6uOV6vTPl5V0eIck6ZMhVynwqqtVcfyYDuUfMzhd8zDnnol59xwxMTEtPrbJRWTmzJl69dVX9eGH\nH6p3797O7WazWZJktVoVFRXl3G61Wp37zGazbDabysrKGlwVsVgsSk1N/dH3S05ObtaJwHXl5+cz\n3x6Iefc8zHnT7Mj9r3ove1xdaipU3qGTLFMe0cgxrrkgnTn3TMy7Z6moqGjxsU1aIzJjxgz9/e9/\n1wcffKD+/fs32BcdHS2z2azc3FzntpqaGm3YsEEjRoyQJCUlJcnPz6/BmJKSEhUVFTnHAADgyWpO\nnlLenCzFLblPXWoqVBQ5RFrwdw1y0RICAI1p9IrItGnTtGrVKr3++usKDQ11rgkJCQlRUFCQvLy8\nlJ6errlz5yo2NlYxMTGaM2eOQkJCNGnSJElSaGiopkyZooyMDJlMJnXt2lWzZs1SfHy80tL4CxYA\n4NmKt2yXV/ZDSjz+teq9fLR11K1K+v0d8vHjucMA3Fejf8M999xz8vLyct529ztZWVl69NFHJUkZ\nGRmqrq7WtGnTVF5eruHDhys3N1dBQUHO8dnZ2fL19dWECRNUXV2ttLQ0rVq16pzWkQAA4MpsNps2\nLV2qwe+9KD+HTaWhkaq7Z45SkoYaHQ0AzrtGi0hTHziYmZmpzMzMs+739/dXTk6OcnJymp4OAAA3\ndeRgiY7Oe0iJpdskSZsSr9age/+gjt/7n3gA4M645gsAQBuy2+3a8o9/qs9rz+rCupMqDwiV5daH\nNOyK0UZHA4A2RREBAKCNHDts1ddPZyl+/2eSpB3RKer5h8c1yGwyOBkAtD2KCAAAbWDL62+o59/m\na1DtCVX5BWrfdfdo6IRfy9u7STewBAC3QxEBAOA8On60TPufflzxuz6WJO3skSDT/Y8psWdUI0cC\ngHujiAAAcJ5s+8/b6rbqacXXHNcp3w7adfVdSrz5Jq6CAIAoIgAAtLpjh6366pk5GrJ7gyRpd7cB\nCr33cSX362NwMgBoPygiAAC0ku/uiBX9jxwNOV2lU74dVHTV75U4+Rb5+PgYHQ8A2hWKCAAArcD6\n1dc6suCPiv+qQNKZtSAX3JulYb17GpwMANoniggAAOfAVlevTStXKvbtFzSgvkaV/sHaf910Jdxw\nPWtBAOAnUEQAAGih4s1bVb9kjhK/2StJ2t5vhHrc+4iSupkNTgYA7R9FBACAZjpVWaUdi5/R0M//\nJR85dDSwqywT79XQq68yOhoAuAyKCAAAzbDt7XcVvurPSjp1VDZ5aVPKeA24e6biO4UYHQ0AXApF\nBACAJrAUfyXronkatO9TSdL+C/rK565HNCwx3uBkAOCaKCIAAPyE2upqbV36ggZ9/FeF2+p0yjdA\nRVf8TomTJ8vHj39GAaCl+BsUAICz2PHu++r81wVKOmGRJG2NvVQ9787QsKjuBicDANdHEQEA4P+w\n7Dsg6+I/adD+zyVJX3fuqZrf3q/ESy8xOBkAuA+KCAAA3zpVWaUdS5/ToP+9pnB7nU76dtSuMbdq\n6K23ys/f3+h4AOBWKCIAAI9ns9m05bV/qse//6Kk6nJJ0tbYy9RjeoaGRXYzOB0AuCeKCADAo+3+\n9At5LZuvod/skSTtD4+RY8p9ShyeYnAyAHBvFBEAgEeyHDioQ39ZoPiijyVJZR27qOTaOxR//a/k\n4+NjcDoAcH8UEQCAR6ksK1fRi3/RoM//pXB7nWp9/LTjZzdo0O/uUGKnYKPjAYDHoIgAADzC6Zoa\nbX3lFfVbt1JJp6sknVkH0v2OmUrp3dPgdADgeSgiAAC3Zrfbte31NxS+9i9KrLJKknZ1H6wOU2Yq\nMSnB4HQA4LkoIgAAt2S327Xzg4/lv2axBh/dJ0n6unMPVd04XQPGjJK3t7fBCQHAs1FEAABuZ/en\nX6h+5ULFHtoh6cxC9K9/OUXxE26Qr5+fwekAABJFBADgRg5s26HKZQs18MAXkqRK/2DtGTlRg2+Z\nrKSgIIPTAQC+jyICAHB5B3fu0tHlz2nwrk/UQw6d8u2gnRf/SnG33a6ULp2NjgcA+BEUEQCAyyrd\ntVfW5Us0eOdHipRDp719tSNprPredodSukUYHQ8A8BMoIgAAl1O6e68sy/+iwYUfyCyH6rx8tCX+\nF+px6+81jFvxAoBLoIgAAFxGSdEeWVe+oEGFH8jssKvey0ebh4xR1G/voIAAgIuhiAAA2r3iLTt0\nfPULGrR7g7rJoXovb20e8gtFTb5DyX17Gx0PANACFBEAQLu157M8Va95UQMPfKGekk57++rLIaMV\nefNUCggAuDiKCACgXbHb7Sr6aL0ca19WbMlWSVKNj78Kk36p6Ft+p+So7gYnBAC0BooIAKBdqK+r\n0/Y3/qNOb6/ShWX7JUlVfoHaddE1irnptxpmNhmcEADQmrwbG/DJJ5/ommuuUVRUlLy9vbVixYof\njMnKylJkZKQCAwM1cuRIFRYWNthfW1ur6dOnKzw8XMHBwRo3bpxKS0tb7ywAAC6r+sRJ5b20TEd+\nO1ZDVj6m3mX7dSygswrSbpP3C28p5b4/qAslBADcTqNF5OTJkxoyZIieffZZdezYUV5eXg32z5s3\nTwsWLNCiRYuUl5cnk8mk0aNHq6qqyjkmPT1da9eu1Zo1a7R+/XpVVlZq7NixstvtrX9GAACXcLTU\noi8WzFfNlCuU+OZCRVQd0aFO3bTl+vvU6eW3lXL3dIXwMEIAcFuNfjTryiuv1JVXXilJuvXWWxvs\nczgcys7O1uzZszV+/HhJ0ooVK2QymbR69WpNnTpVFRUVWrZsmZYvX65Ro0ZJklauXKlevXpp3bp1\nGjNmTCufEgCgPdu/aauO/+MVDdj5sZIcNknSPlN/VV89WYOuGKMefnxqGAA8wTn9bV9cXCyr1dqg\nTAQEBCg1NVUbN27U1KlTVVBQoLq6ugZjoqKiFBcXp40bN1JEAMAD2OrqVbZ5q3a+lK3+h79UL0k2\neWl7v5+p4/ibFHNxiry9G71IDwBwI+dURCwWiyQpIiKiwXaTyaRDhw45x/j4+CgsLKzBmIiICFmt\n1nN5ewBAO3f8aJl2v/Z3Rf3v30qrOiLp2wXoiVeq5403a2h0L4MTAgCMct6uf//ftSTNlZ+f30pJ\n4AqYb8/EvLuv48UH1fGLj5Sw7zMl2eskSYeCI7RzaJqCU4bLPzBAJWXfqKTsG4OT4nzj+9wzMe+e\nIyYmpsXHnlMRMZvNkiSr1aqoqCjndqvV6txnNptls9lUVlbW4KqIxWJRamrqWV87OTn5XKLBheTn\n5zPfHoh5dz81J0+p8K23FLzuH0o4sluSZJeXvuydLK8rJ6i6SyeNSUkxOCXaEt/nnol59ywVFRUt\nPvacikh0dLTMZrNyc3OVlJQkSaqpqdGGDRs0f/58SVJSUpL8/PyUm5uriRMnSpJKSkpUVFSkESNG\nnMvbAwDagYNfFsn6+qvqt/V9xZ8+c8fESv9g7Um8QlHXT9SQfn0k8X9IAQANNVpETp48qT179kg6\n87Tbr776Slu2bFFYWJh69Oih9PR0zZ07V7GxsYqJidGcOXMUEhKiSZMmSZJCQ0M1ZcoUZWRkyGQy\nqWvXrpo1a5bi4+OVlpZ2fs8OAHBefHf1I/CD1xVjKVTkt9v3h/dT5WXjNeDacUoJCjI0IwCgfWu0\niOTl5enyyy+XdGbdR2ZmpjIzM3Xrrbdq2bJlysjIUHV1taZNm6by8nINHz5cubm5CvreP0DZ2dny\n9fXVhAkTVF1drbS0NK1ateqc15EAANqO3W7XgS3bVfbWWvXb/oHz6sdJ347aNfhyhV39a8Ukxhuc\nEgDgKhotIpdddlmjDx78rpycjb+/v3JycpSTk9P8hAAAQx0/UqY9b/xLXTe+pd7HDui7+1ztD++n\nikuv1YBrxmlYp2BDMwIAXA9PjQIA/EDd6dMq+u+Hqv/gP4rb+6kSv33wYEWHEO0bMkoX/PI6xQwd\nbHBKAIAro4gAACSd+ehVccEWHXvvDUVv/1ADaislnXnwYGGvJNkuH6cBvxitYQEBBicFALgDiggA\neDhL8Vc6+J83ZM5/T70rStX72+0lnaNkTf6F+lx7nQZHdTcyIgDADVFEAMADlR22aP9bb6lT3jr1\nsxYp/Nvt5QGhKh58ubpecY2iE4aol7e3oTkBAO6LIgIAHqKi7Jj2vPOOOn6aq5jS7UqUQ5JU7dNB\nu/sNl9/lYxV3+aUa5udncFIAgCegiACAGzt+tEz73ntP/p9/oP5fb1HCt4vO67x8tDM6WfafXaHY\nMaOVGMIzPwAAbYsiAgBuptxyRPtycxWQ94FiSrZqqOPMLdjrvby1M2qoakeMUcwvrlB8WBeDkwIA\nPBlFBADcwOF9B1Ty31wFb/5EfQ8XKuHbj13Ve/loZ89E1aRcrr5jxmiQKbyRVwIAoG1QRADABdnt\ndn21dbu++ei/umD7evU6dkCmb/fVefmoqFei6i8apb5jRmtQWFdDswIA8GMoIgDgImpOntLe9f9T\n9ecfq8euz9TzVJl6frvvpG9H7e2bIu+LLlPfy0dqSOdQQ7MCANAYiggAtGPfHCzRVx99KP9NG9T3\n6y2Ks5127ivr2EVfxY5Q4M8uV8wlI5TIgwYBAC6EIgIA7cjpmhrt3fi5qr7YIFPRF+px/KC+/8Gq\n/eH9VD5whMJ+PlK9hw6WycfHsKwAAJwLiggAGMhut+v/tXfnsVGd9xrHvzPjWT32eB3vxDZhdVkc\niHPrcFNaJai0KmrVQpS2KUorRVXTFEKkViVUolJiWrVKGwhOUFShSFUVqnulpEtUJRK0wReSm1xs\nEjAYEhYbsAe8je2xZzvn3D/sODgQoAmeMfbzkV7N+J33HJ7hyJ73N2e7cOIDOg804X73INUd7zLP\niI2/Ppzh4dRttSTvWMFtK7/EnLLiNKYVERG5eVSIiIikWE9nF+1N/4N5+C3KTh2iZLiHksteP5Nf\nRc/8u8j+jxXMvvNOlnhcacsqIiIyWVSIiIhMsnBPL2cPvkns8NsUvH+IWX3t5Fz2ep87m/bKWqza\nem77z3uYXVbM7LSlFRERSQ0VIiIiN9lATx9n3nqL6OG3yX+/mVk9Z6gZu68HQNTh4lTZIkZq7qTw\nrruZtWgBBTrXQ0REZhgVIiIin1H3uQt0vPW/JI8eouD0u5T3tU8oPOL2DM4UzWNw3h1k19ZRfddy\nFukKVyIiMsOpEBER+TcYiSQdx9robjmEre1dis++R/FQiNzLxsTtGZwNzmVg7h1k1dZRVbeMBZm+\ntGUWERGZilSIiIhcw2BfP+3/10zkSAu+D96joquNisQwFZeNiWR4aS9dwPCcpWQvXUbV8lrme71p\nyywiInIrUCEiIjImEY3TfvQofe+2YDt5lMLzxykLn2f+x8ZdzCygs7wGY+5i8muXMWvRQj7ndKYl\ns4iIyK1KhYiIzEiJeJxzrSfobX0P82QrgXNtlPecodJMUHnZuJjdybmCKvqrPodrwRLK7riDklnl\nEy63KyIiIv8+FSIiMu3FRkY4d+QY/W3HMD84TuD8Ccq7TzPLTDDrY2PPB8q4VD4f6/YachcvZVZN\nDXN1Hw8REZGbToWIiEwrvZ0hulqPMfR+GxlnT5B34X1KwueotMwrxnZmFXOpZA7J6gX4532O8qWL\nmJWbc0VxIiIiIjefChERuSUNDwxx4Xgb4fdPYp45ge/CKYounSY3NkDgY2MNbHTkVNBbPBujch6Z\n8xZSvngR5fm5lKclvYiIiKgQEZEpbag/TOjkB/R/8D7G2Q/wdJ6i4NJZiiKXqLraeKePzvxKBkuq\nsVXNI2f+AsprFlCZ6Ztw7oeIiIiklwoREUk7wzC41HGe3lOniZw5BedP4wu1U9DTQcFI71ULiLg9\ng65AGf3B20hU3I539lwKFywgOKucgN2e6rcgIiIi/yYVIiKSEoZh0NcZovv0GSIdZxk80UrLS7sJ\n9JyjaKCTQiNB4VWWi9mddOWUES6oIFFWjbtyNnlz51FyexVVumSuiIjILUuFiIjcNNHIMBfbOwi3\ndxC70IHZdR73pfMEei9QONhFrhGfcAfyy/V6cujOLSNSWIFVXoXntmoKbr+d4Kxyqp36UyUiIjLd\n6NNdRG7YyGCES+0dDJw/R7TzAubFTjJ6OsnsC5EX7iIv2k8ZUPYJy/e7s+kOlBDJL6XPl0PewkVk\nV1ZRNLuKwpzAVfeIiIiIyPSkQkREgNG9GT0XOhnsCjHS1UniYhe2nhDu3hCZAxfJG+omOz50zUIj\nab40/QcAAAyISURBVHNwyV9If6CYaEEJZmEZrrIKsmfdRmF1Jfk5AfLHxr7zzjvULl+eoncnIiIi\nU40KEZFpLh6N0tt1kcGuEMOXLpHouYjZcwl7fzeecDeZg93kRHrJjg9RDBRfa132DLozCxkIBBnJ\nLcIsKMFRVEpmWRm5s2ZRUFpCuTNDl8QVERGR61IhInILGhmMEL54iaHuHkZ6u4n39WL09mAL9+IY\n6MUz2EvmcB/Zw/1kx4cohOse9pSwOej15THgz2ckUEgirwhbQRGuohL8xSXklJeSEyykwuFIxVsU\nERGRaU6FiEgamaZJLDLCYF8fkd4+on39xMJ9JMP9mOF+bIP9OIb6cUXCuIfD+EfCZEcH8BjxGyou\nAJI2O/3eHAYy8xj255EI5GPmFGDPK8RdGMRfXExOSQnZhfmUOhyUTvabFhEREUGFiMhnlkwkiAwM\nMtI/wMjgILFwmPjgAImhIcyhAayhQWyRQezDAzgjg7iig3iiQ/hiQ/hjQ7jNBPkwfu7EjYjZnYQ9\nAYZ8AaKZOcT9ORjZ+ZCThzOvAHd+AZmFhWQXFZKdn0eRw0HRZP0HiIiIiHwKKS1EGhsb+c1vfkNX\nVxc1NTX8/ve/Z8WKFamMIIJhGMQiI0QjEWLDI8QiQ8QjwyQiEZKRIZIjI1gjEczhYRgZwjYyjD0a\nwREdJiM2gisWwRMfxhMfxhcfxmvE8AP+T5knZncy6M5i2ONnxJtN3JtFMjOAmZUD2Tlk5OTizMnD\nm5uLLy+XQLAQrz8Tn27aJyIiIrewlBUie/bsYePGjTz33HOsWLGCnTt3snr1alpbW6moqEhVDJmC\nLNMiHo0Sj8ZIRGPEYzGSsdHnRixGIhrFiMUwYlHMWGysRbHiMaxoFOJRSMSwxaLYEzEcsSiORBRH\nIoZzrLmSUVyJKJ5kDK8Rww24b1J+AxvDTh8jLi8jbj8xdyZxTyZJbxamz4/l82PzZ2P3Z+MMBHBl\nB/AEAmTmBPDn5uLJ9KqoEBERkRknZYXI008/zUMPPcQPfvADALZv384//vEPnnvuORoaGlIV45Zn\nmiaWYWIaJqZlYCQMDNPANC5riSSGaWImk1iGiZFMYhpJLMPATBoYRhIraYz2ffwxkcAyDKxkEstI\njj4mE5AcfY4x1sae25IJbEYCm2GMPU9iN5LYjASOZAK7kcBhJMdagozLmtNM4jQS3GEmAPCOtVQY\ncbiJZriJOz3EnB7iLh8Jl5eky4Ph8WG6fVgeL3j92Lw+HJl+HD4fGZl+nP5MPFnZ+ALZ+LKz8fh9\n5Njt5KQou4iIiMh0kJJCJB6Pc+jQIX76059O6F+1ahUHDhy4oXU0/9d/w9tvXPmCZX34ZMKD7Yp+\nC9uEZazRPuvDMdboMhbYxl4bXY/50bLWh+PMj5azRpezWebYmNHH0THWRw1zfBwW2C1zfIydj8ba\nLXNC+7DPMTbeMfZ+7GNtupzkk7A5iDtcxB1OEg4XSYcTw+Ek4XSRzHBhZLgxnC7MDBem043pdGG5\nPeB0g8sDbjd2t2e0eb1keLw4vB6cXh9OrxeXz4fb58Xr9+PyefA7HJ/6UCoRERER+exSMo/t7u7G\nMAyKiiaeLhsMBunq6rrqMuFweMLP1ffdC/fdO2kZpyJz7NEAEukMkkLOsTZZ4maS+NDQJP4LcqPm\nzJlzxe+5TG/a5jOPtvnMpO0uN0oHpouIiIiISMqlpBApKCjA4XAQCoUm9IdCIUpKSlIRQURERERE\nppCUHJrlcrlYtmwZr732Gt/85jfH+19//XXWrl07/nMgEEhFHBERERERSbOUneu8adMmHnzwQerq\n6qivr+f555+nq6uLH/7wh6mKICIiIiIiU0TKCpF169bR09PDk08+SWdnJ4sWLeLVV1/VPURERERE\nRGYgm2WNX+dWREREREQkJabUVbMaGxupqqrC6/WyfPlympqa0h1JJsm2bdu48847CQQCBINB1qxZ\nw9GjR9MdS1Jo27Zt2O12Hn300XRHkUnW2dnJ+vXrCQaDeL1eampqeOONq9wXSqaFZDLJ5s2bqa6u\nxuv1Ul1dzS9+8QsMw0h3NLlJ3njjDdasWUN5eTl2u50XX3zxijFbt26lrKwMn8/HF7/4RVpbW9OQ\nVG6Wa23zZDLJz372M5YsWYLf76e0tJTvfOc7dHR0XHe9U6YQ2bNnDxs3bmTLli20tLRQX1/P6tWr\nb+hNyK3nX//6Fz/+8Y85ePAge/fuJSMjg3vvvZe+vr50R5MUePPNN3nhhRdYvHgxNpvt+gvILau/\nv5+7774bm83Gq6++yvHjx3n22WcJBoPpjiaTpKGhgV27drFjxw7a2tp45plnaGxsZNu2bemOJjdJ\nJBJh8eLFPPPMM3i93iv+jv/617/m6aef5tlnn+Xtt98mGAxy3333MaT7eN2yrrXNI5EIzc3NbNmy\nhebmZl555RU6Ojr48pe/fN0vIKbMoVl33XUXS5cuZdeuXeN9c+fO5Vvf+hYNDQ1pTCapEIlECAQC\nvPLKK3z1q19NdxyZROFwmGXLlvGHP/yBrVu3smjRIrZv357uWDJJNm/ezP79+9m/f3+6o0iKfO1r\nX6OgoIDdu3eP961fv56+vj7+8pe/pDGZTIasrCx27tzJ9773PQAsy6K0tJSf/OQn/PznPwcgGo0S\nDAb57W9/y8MPP5zOuHITfHybX82xY8eoqanhvffeo6am5hPHTYk9IvF4nEOHDrFq1aoJ/atWreLA\ngQNpSiWpNDAwgGma5ObmpjuKTLKHH36YtWvX8oUvfIEp8j2ITKKXX36Zuro67r//foqKiqitrWXn\nzp3pjiWTaPXq1ezdu5e2tjYAWltb2bdvH1/5ylfSnExS4fTp04RCoQlzOo/Hwz333KM53QwSDocB\nrjuvS9lVs66lu7sbwzAoKiqa0B8MBunq6kpTKkmlDRs2UFtby+c///l0R5FJ9MILL3Dq1Cn+9Kc/\nAeiwrBng1KlTNDY2smnTJjZv3kxzc/P4eUGPPPJImtPJZPjRj37EuXPnWLBgARkZGSSTSbZs2aLL\n9c8QH87brjanu3DhQjoiSYrF43Eef/xx1qxZQ2lp6TXHTolCRGa2TZs2ceDAAZqamjQxncba2tp4\n4oknaGpqwuFwAKO78LVXZHozTZO6ujqeeuopAJYsWcLJkyfZuXOnCpFpavv27ezevZuXXnqJmpoa\nmpub2bBhA5WVlXz/+99PdzxJI33GT3/JZJLvfve7DAwM8Le//e2646dEIVJQUIDD4SAUCk3oD4VC\nlJSUpCmVpMJjjz3Gn//8Z/bt20dlZWW648gkOnjwIN3d3ROOFTUMg/3797Nr1y4ikQhOpzONCWUy\nlJaWsnDhwgl98+fPp729PU2JZLI99dRTbNmyhXXr1gFQU1PD2bNn2bZtmwqRGaC4uBgYncOVl5eP\n94dCofHXZHpKJpM88MADHD16lH/+8583dLj9lDhHxOVysWzZMl577bUJ/a+//jr19fVpSiWTbcOG\nDezZs4e9e/cyd+7cdMeRSfaNb3yDI0eOcPjwYQ4fPkxLSwvLly/ngQceoKWlRUXINHX33Xdz/Pjx\nCX0nTpzQFw/TmGVZ2O0Tpxd2u117P2eIqqoqiouLJ8zpotEoTU1NmtNNY4lEgvvvv58jR46wb9++\nG74y4pTYIwKjh+c8+OCD1NXVUV9fz/PPP09XV5eOKZ2mHnnkEf74xz/y8ssvEwgExo8pzcrKIjMz\nM83pZDIEAgECgcCEPp/PR25u7hXfmMv08dhjj1FfX09DQwPr1q2jubmZHTt26FKu09jXv/51fvWr\nX1FVVcXChQtpbm7md7/7HevXr093NLlJIpEIJ0+eBEYPvzx79iwtLS3k5+dTUVHBxo0baWhoYP78\n+cyZM4cnn3ySrKwsvv3tb6c5uXxa19rmpaWlrF27lnfeeYe//vWvWJY1Pq/LycnB4/F88oqtKaSx\nsdGqrKy03G63tXz5cmv//v3pjiSTxGazWXa73bLZbBPaL3/5y3RHkxRauXKl9eijj6Y7hkyyv//9\n79aSJUssj8djzZs3z9qxY0e6I8kkGhoash5//HGrsrLS8nq9VnV1tfXEE09YsVgs3dHkJtm3b9/4\n5/bln+UPPfTQ+JitW7daJSUllsfjsVauXGkdPXo0jYnls7rWNj9z5swnzutefPHFa653ytxHRERE\nREREZo4pcY6IiIiIiIjMLCpEREREREQk5VSIiIiIiIhIyqkQERERERGRlFMhIiIiIiIiKadCRERE\nREREUk6FiIiIiIiIpJwKERERERERSbn/BxKZXADfmv0TAAAAAElFTkSuQmCC\n", + "text": [ + "" + ] + } + ], + "prompt_number": 27 + }, { "cell_type": "heading", "level": 2, diff --git a/experiements/RungeKutta.py b/experiements/RungeKutta.py index 0a2fa8c..b995d5b 100644 --- a/experiements/RungeKutta.py +++ b/experiements/RungeKutta.py @@ -76,11 +76,8 @@ def rk4(y, x, dx, f): return y + (k1 + 2*k2 + 2*k3 + k4) / 6. -def rk2 (y,x,dx,f): - """computes the 2nd order Runge-kutta for dy/dx""" - - + def fx(x,t): return fx.vel @@ -126,7 +123,44 @@ def ball_scipy(y0, vel, omega, dt): ys.append(solver.integrate(t)) +def RK4(f): + return lambda t, y, dt: ( + lambda dy1: ( + lambda dy2: ( + lambda dy3: ( + lambda dy4: (dy1 + 2*dy2 + 2*dy3 + dy4)/6 + )( dt * f( t + dt , y + dy3 ) ) + )( dt * f( t + dt/2, y + dy2/2 ) ) + )( dt * f( t + dt/2, y + dy1/2 ) ) + )( dt * f( t , y ) ) + +def theory(t): return (t**2 + 4)**2 /16 + +from math import sqrt +dy = RK4(lambda t, y: t*sqrt(y)) + +t, y, dt = 0., 1., .1 +while t <= 10: + if abs(round(t) - t) < 1e-5: + print("y(%2.1f)\t= %4.6f \t error: %4.6g" % (t, y, abs(y - theory(t)))) + + t, y = t + dt, y + dy(t, y, dt) + +t = 0. +y=1. + +def test(y, t): + return t*sqrt(y) + +while t <= 10: + if abs(round(t) - t) < 1e-5: + print("y(%2.1f)\t= %4.6f \t error: %4.6g" % (t, y, abs(y - theory(t)))) + + y = rk4(y, t, dt, test) + t += dt + if __name__ == "__main__": + 1/0 dt = 1./30 y0 = 15. diff --git a/experiements/__init__.py b/experiements/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/experiements/quaternion.py b/experiements/quaternion.py new file mode 100644 index 0000000..b7f4929 --- /dev/null +++ b/experiements/quaternion.py @@ -0,0 +1,76 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Jan 29 18:36:36 2015 + +@author: rlabbe +""" + +import numpy as np +from math import sin, cos, atan2, asin + +''' +psi = yaw +phi = roll +theta = pitch +''' + + +def e2q(vector): + + roll = vector[0] + pitch = vector[1] + heading = vector[2] + sinhdg = sin(heading/2) + coshdg = cos(heading/2) + + sinroll = sin(roll/2) + cosroll = cos(roll/2) + + sinpitch = sin(pitch/2) + cospitch = cos(pitch/2) + + q0 = cosroll*cospitch*coshdg + sinroll*sinpitch*sinhdg + q1 = sinroll*cospitch*coshdg - cosroll*sinpitch*sinhdg + q2 = cosroll*sinpitch*coshdg + sinroll*cospitch*sinhdg + q3 = cosroll*cospitch*sinhdg - sinroll*sinpitch*coshdg + + return np.array([q0, q1, q2, q3]) + + +def q2e(q): + roll = atan2(2*(q[0]*q[1] + q[2]*q[3]), 1 - 2*(q[1]**2 + q[2]**2)) + pitch = asin(2*(q[0]*q[2] - q[3]*q[1])) + hdg = atan2(2*(q[0]*q[3] + q[1]*q[2]), 1-2*(q[2]**2 + q[3]**2)) + + return np.array([roll, pitch, hdg]) + + +def e2d(e): + return np.degrees(e) +def e2r(e): + return np.radians(e) + +def add(q1,q2): + return np.multiply(q1,q2) + + + +def add2(q1, q2): + Q1 = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3] + Q2 = q1[0]*q2[1] + q1[1]*q2[0] + q1[2]*q2[3] - q1[3]*q2[2] + Q3 = q1[0]*q2[2] + q1[2]*q2[0] + q1[3]*q2[1] - q1[1]*q2[3] + Q4 = q1[0]*q2[3] + q1[3]*q2[0] + q1[1]*q2[2] - q1[2]*q2[1] + + return np.array([Q1, Q2, Q3, Q4]) + + +e = e2r([10, 0, 0]) +q = e2q(e) +print(q) +print(e2d(q2e(q))) +q2 = add2(q,q) +print(q2) +e2 = q2e(q2) +print(e2d(e2)) + + \ No newline at end of file