Playing+with+Ellipses (5).ipynb

{"nbformat_minor": 2, "cells": [{"source": "# Playing with Ellipses\n\n$$\n    \\newcommand{\\vect}[1]{\\mathbf{#1}}\n    \\newcommand{\\hypot}[2]{\\sqrt{#1^2 + #2^2}}\n$$\n\nAn ellipsis centered at the origin can be parameterized as\n$$\n    \\vect{p} = \\vect{a} \\sin\\theta + \\vect{b} \\cos\\theta\n$$\nWhere $\\vect{a}$ and $\\vect{b}$ are orthogonal to each other.\n\nThus, on the $y$ axis,\n\\begin{align\\*}\n    p_y &= a_y \\sin\\theta + b_y \\cos\\theta \\\\\\\\\n        &= \\hypot{a_y}{b_y} \\left( \n            \\frac{a_y}{\\hypot{a_y}{b_y}} \\sin\\theta + \\frac{b_y}{\\hypot{a_y}{b_y}} \\cos\\theta\n        \\right)\n\\end{align\\*}\n\nPicking $\\alpha$ such that $\\cos\\alpha = \\frac{a_y}{\\hypot{a_y}{b_y}}$ we then have\n\\begin{align\\*}\n    p_y &= \\hypot{a_y}{b_y} ( \\cos\\alpha \\sin\\theta + \\sin\\alpha \\cos\\theta ) \\\\\\\\\n        &= \\hypot{a_y}{b_y} \\sin(\\alpha + \\theta)\n\\end{align\\*}\n\nSolutions for the above are\n$$\n    \\theta = \\arcsin\\left( \\frac{p_y}{\\hypot{a_y}{b_y}} \\right) - \\alpha + 2k\\pi\n    \\quad \\text{for $k$ in $\\mathbb{Z}$}\n$$\nand\n$$\n    \\theta = \\pi - \\arcsin\\left( \\frac{p_y}{\\hypot{a_y}{b_y}} \\right) - \\alpha + 2k\\pi \n    \\quad \\text{for $k$ in $\\mathbb{Z}$}\n$$", "cell_type": "markdown", "metadata": {}}, {"execution_count": 54, "cell_type": "code", "source": "%pylab inline\nimport ipywidgets", "outputs": [{"output_type": "stream", "name": "stdout", "text": "Populating the interactive namespace from numpy and matplotlib\n"}], "metadata": {"collapsed": false}}, {"execution_count": 24, "cell_type": "code", "source": "#         double ay = a.y;\n#         double by = b.y;\n\n#         if ((ay*ay + by*by) < 0.0) {\n#             return {0.0, 0.0};\n#         }\n\n#         double alpha = acos(ay / sqrt(ay*ay + by*by));\n#         double theta = asin((y-circle_center.y)/sqrt(ay*ay + by*by));\n#         double theta1 = theta - alpha;\n#         double theta2 = M_PI - theta - alpha;\n\n#         double z1 = circle_center.z + radius*(sin(theta1)*a.z + cos(theta1)*b.z);\n#         double z2 = circle_center.z + radius*(sin(theta2)*a.z + cos(theta2)*b.z);\n        \ndef find_z(ay, by, cy, y, az, bz, cz):\n    alpha = arccos(ay / sqrt(ay*ay + by*by))\n    theta = arcsin(  (y - cy) / sqrt(ay*ay + by*by)  )\n    theta1 = theta - alpha\n    theta2 = pi - theta - alpha\n    \n    z1 = cz + az*sin(theta1) + bz*cos(theta1)\n    z2 = cz + az*sin(theta2) + bz*cos(theta2)\n    return z1, z2\n", "outputs": [], "metadata": {"collapsed": true}}, {"execution_count": 27, "cell_type": "code", "source": "def rotation_matrix(axis, theta):\n    \"\"\"\n    Return the rotation matrix associated with counterclockwise rotation about\n    the given axis by theta radians.\n    \"\"\"\n    axis = asarray(axis)\n    axis = axis/sqrt(dot(axis, axis))\n    a = cos(theta/2.0)\n    b, c, d = -axis*sin(theta/2.0)\n    aa, bb, cc, dd = a*a, b*b, c*c, d*d\n    bc, ad, ac, ab, bd, cd = b*c, a*d, a*c, a*b, b*d, c*d\n    return np.array([[aa+bb-cc-dd, 2*(bc+ad), 2*(bd-ac)],\n                     [2*(bc-ad), aa+cc-bb-dd, 2*(cd+ab)],\n                     [2*(bd+ac), 2*(cd-ab), aa+dd-bb-cc]])", "outputs": [], "metadata": {"collapsed": true}}, {"execution_count": 47, "cell_type": "code", "source": "rx, ry, rz, rd = 0.14, 0.21, 0.08, 2.08\n\n# @ipywidgets.interact(rx=(0, 1, 0.01), ry=(0, 1, 0.01), rz=(0, 1, 0.01), rd=(0, 2*pi, 0.01))\n# def go(rx, ry, rz, rd):\n\nax, ay, az = [0, 3, 0] @ rotation_matrix([rx, ry, rz], rd)\nbx, by, bz = [0, 0, 6] @ rotation_matrix([rx, ry, rz], rd)\ncx, cy, cz = 0, -10, 10\n\ntheta = arange(0, 2*pi, 0.001)\ny = cy + ay*sin(theta) + by*cos(theta)\nz = cz + az*sin(theta) + bz*cos(theta)\nplot(y, z, 'r')", "outputs": [{"execution_count": 47, "output_type": "execute_result", "data": {"text/plain": "[<matplotlib.lines.Line2D at 0x7fd3c5e16278>]"}, "metadata": {}}, {"output_type": "display_data", "data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VVXWx/HvgoCADqgIYoVBSoCAgIg0IQrqWF5BwV5H\n1FFHUWEUEQRERUfGrjjqgB0bVuxlCIoUFSlCABkdy1hQx4ptBPb7xwoaEVLuPeckOfl9nicPyU1y\n9z4Pedbdd52117YQAiIikg41KnoCIiISHQV1EZEUUVAXEUkRBXURkRRRUBcRSREFdRGRFCk1qJvZ\nJDNbaWaLij02zswWmtl8M3vazJrEO00RESkLK61O3cx6AauAO0IIHYoe2yyEsKro8zOAtiGEU+Oe\nrIiIlKzUlXoIYSbwxXqPrSr25abA2ojnJSIiGcjJ9BfN7GLgWOBLYI/IZiQiIhnL+EZpCGFUCGFH\n4G7gjOimJCIimcp4pV7MFOBJYOyGvmlmai4jIpKBEIKV93fKulK3og//wqxFse8NAJaWMrHUfowZ\nM6bC56Dr07Xp+tL3kalSV+pmNgXIBxqa2XvAGGB/M2sNrAHeBU7JeAYiIhKZUoN6COHIDTx8awxz\nERGRLGlHaZby8/MregqxSvP1pfnaQNdXXZW6+SjrAcxC3GOIiKSNmRFivFEqIiJVgIK6iEiKKKiL\niKSIgrqISIooqIuIpIiCuohIiiioi4ikiIK6iEiKKKiLiKSIgrqISIooqIuIpIiCuohIiiioi4ik\niIK6iEiKKKiLiKSIgrqISIooqIuIpIiCuohIiiioi4ikiIK6iEiKKKiLiKSIgrqISIooqIuIpIiC\nuohIiiioi4ikSE5FTyBrIcBf/wqFhbDttrDbbtCrFzRqVNEzExFJXKkrdTObZGYrzWxRsccuN7Ol\nZrbAzB40s/rxTrMEH34II0ZA795Qrx7ccgu0bAl9+sA118D771fY1EREklaW9MutwD7rPfYs0C6E\n0BFYAYyIemJltu22sO++sHw5jB4NTz4JH38M55wDCxZAp06w555w993w/fcVNk0RkSRYCKH0HzJr\nCkwLIXTYwPcGAANDCMds5HdDWcbIyn//Cz17whFHeGA3++V7P/4I06bBpEnwyitw2GEweDB07vzr\nnxMRqUTMjBBCuYNUFDdKTwCeiuB5MtewIcyYAVOnwtlnw5o1v3xvk01g0CB46ilfuW+zjX+9665w\n++3www8VN28RkYhltVI3s5FA5xDCwBJ+N/6V+jpffAGHHOKB/J57oP5GUv1r18LTT8O118L8+XDy\nyXDKKbDddsnMU0SkFJmu1DOufjGz44H9gD1L+9mxY8f+/Hl+fj75+fmZDluyLbbwFfmQIdC9Ozz0\nELRu/dufq1ED9tvPP5Yvh+uvh/btYe+9f/ldpWZEJEEFBQUUFBRk/TxlXak3w1fq7Yu+/gNwBdA7\nhPDfUn43uZX6OiF4FczIkR6wDzus9N/56iu47Tb/+fr14S9/8VV/TtWv+hSRqifTlXqpQd3MpgD5\nQENgJTAGOB+oDawL6HNCCKdt5PeTD+rrzJ/vgXmffeDKKz0tU5q1a72CZsIEePddz9EPHgybbRb/\nfEVEisQW1LNVoUEdfAU+eDD8+9+eZ2/Vquy/O3euB/cZM+BPf4IzzoCtt45vriIiRSqy+qVya9AA\nHngATjgBevSAm2/29ExZ7LabV9TMmuVlk7m5HtzffDPeOYuIZCj9K/XiCgvhqKOgaVPPuZe3lcAn\nn8ANN8DEibD77r7BqXv3eOYqItWaVupl0bYtzJnjFTEdO3pZY3k0bgwXXgjvvOO7VI880tsRPPdc\n2Vf/IiIxql4r9eKmT4fjjoP+/eGyy2DTTcv/HKtXe57+0kv9RuqoUXDAAV4yKSKSBa3Uy2uPPWDh\nQt+w1LEjzJxZ/ufIyYFjjoHFi2H4cBgzxp/r3nt/vatVRCQh1XelXtwjj8Bpp8Hhh8PFF3u3x0yE\n4JufLrkEPv3Uu0cefTTUqhXtfEUk9bRSz8aAAbBoEXz0ka+0Z83K7HnMfJfqzJleZXP33dCihd9Y\nVY8ZEUmAVurre/BBOP10X2GPGwd162b3fHPm+Mp93jwYOtR7zGgjk4iUQiv1qAwc6Kv2d9/19rxz\n52b3fN26eevfp56CV1+F5s09yH/9dTTzFREpRkF9Qxo1gvvv9/LF/v29VcC332b3nDvvDPfdBy++\nCMuWwU47KbiLSOQU1Ety6KFe2fLZZ5CX5/Xo2crNhTvv9Ly7gruIRExBvTRbbeVBeOJEOOkkOP54\n+Pzz7J+3dWsFdxGJnIJ6We27L7zxhrflzcvz9EwUN4AV3EUkQqp+ycTs2d75sWVLX8FHeWLS8uVe\nK//003DWWd4ZcmMnOIlIaqn6JUndu3uv9k6dvK79ppu8D3sUtHIXkSxopZ6txYvhxBOhTh0P7hs6\nPi8bWrmLVEtaqVeUvDx4+WU46CDo2dPLIH/8Mbrn39DKffx4WLUqujFEJDUU1KNQsyaceSYsWOBN\nwjp08C6QUSoe3Bcv9vYDV14J338f7TgiUqUpqEdp++3hoYfg8su9re/xx3uNe5Rat4YpU7xmfuZM\nD+433BDtuwMRqbIU1OPQv7+fstSwIbRrB7feGv0hGu3b+wvIY4/5QdmtWsGkSfDTT9GOIyJVim6U\nxu311/1c03r14O9/hzZt4hln1iy44AJ47z0YO9bbCNesGc9YIhI73SitrDp39k6NhxwCvXvD6NHx\ntOHt0QNeeMErcCZO9Lz+gw9GV2opIlWCVupJ+uADv6G6cCHceCP06xfPOCF4CeSoUR7UL7oI9t/f\n+72LSJWQ6UpdQb0iPP6492zv1csrWBo3jmecEODRRz0ts+mmXu/et6+Cu0gVoPRLVXLAAbBkCWyz\njde5T5wYz5mmZn6q08KFvnHpz3/2s1lfein6sUSkUtBKvaItXuzno373nadkdt01vrFWr4a77vIN\nUq1be1omzvFEJGNaqVdVeXkwYwYMGQIHHginnhpNa98Nycnx2vnly30Ff/DBvhN2yZJ4xhORxCmo\nVwZmcOyxsHSplyG2beu17XFVrtSu7Welvvmm5/X33NM3S73zTjzjiUhiFNQrk803h+uvhyee8Jr2\n3r39vNS41K0Lw4bBihXQrBnssos3DFu5Mr4xRSRWpQZ1M5tkZivNbFGxxwaZ2WIzW2NmneOdYjW0\nyy7es/3YY73scejQeFvv1q/vefbi7xRGjYIvv4xvTBGJRVlW6rcC+6z32BvAQcCMyGckrkYNOPlk\nz3d/+aUH2nvvjb7dQHGNG8PVV/su2A8/9NYDl1/uN3FFpEooNaiHEGYCX6z32PIQwgpABc9xa9QI\nJk+G++6DSy+FvfbyG51xatrUx5wxA+bO9eB+003qKyNSBSinXlX07Anz5nmNe8+eMHJk/CvoNm28\n1cDDD8PUqf5u4Z571HpApBIrU526mTUFpoUQOqz3+HRgWAjh9RJ+N4wZM+bnr/Pz88nPz894woKn\nRoYN87z71Vd7V8gkdon+858wYoS3+b3kEthvP+1OFYlIQUEBBQUFP3994YUXxtcmINugrs1HMXnh\nBa9W2XFHuOaa6I/S25AQ4JFH/Ebqllv6KUy77x7/uCLVTNybj4yN58+1VKsofft6C4C99/aUzHnn\nxX/MnZlvWFq0yM9mPeYYbxa2YEG844pImZSlpHEKMAtoZWbvmdkfzWyAmb0PdAMeN7On4p6obESt\nWl7y+MYb8NFHkJvree+43x3VrOkblpYvhz/8AfbdF444wmveRaTCqPdL2rz8sneArF8frrvO+6on\nYdUqTwFddZX3jh892huWiUhG1PtFXM+e8NprfvJRv37eUyaJTUSbbeYVOcuX+ylPeXke2OPcNCUi\nv6GgnkY1a3pjsMJC+N//PCUzeXIypYgNG8IVV/gGpnffhZYt4dprdTC2SEKUfqkO5s3zlMzatd5b\nJsl2u4sW+Q3cZcv8kI7DD/fdsiJSIp18JCVbuxbuuMPrzA84wEsRGzVKbvzp02H4cO/p/te/+s5Y\nEdko5dSlZDVqeC/1Zcs8/922ra/aV69OZvw99vCWA+ef7ycw7bWXv4MQkUgpqFc3DRp4hcr06fDQ\nQ94RMqnj7cxg0CBvUnbwwf6O4cgj4e23kxlfpBpQUK+u8vJ8R+rIkR5YjzoKPvggmbFr1fIbuStW\neH+ZXXf1Kp1PPklmfJEUU1Cvzszg0EO9j3qzZl7Tfskl8MMPyYy/2WZwwQU+vpmnhC66KP5dsSIp\npqAuHlwvuQRefdVr3Nu29c6MSd3gbtzYNy7NnesBvlUrP4RbrX5Fyk3VL/Jbzz8PZ50FW2/twTYv\nL9nxX3/dK2XeecerdAYNUjdIqXZU0ijRWr3az0kdN85TNOPGeVfGJD33nAf3nBw/gUktm6UaUUmj\nRCsnxzcsLV3qaZjcXLjhhuRKIMHLHl97Dc4+G044wfu3L16c3PgiVZCCupSsYUMP5s8/76cfderk\nh2UkpUYN7/64bBnssw/suae3/P3ww+TmIFKFKKhL2XTo4MH8wgth8GAYOBD+/e/kxq9dG848E958\n019o2rf3hmHffJPcHESqAAV1KTsz3zRUWAidO3t9+ahR8O23yc1h8829zcDrr/uLSqtWnvtPMi0k\nUokpqEv51a3rm5YWLPDAmpsLd9+dXAkkQNOmcOed8MQT8MADvnJ/7LFk5yBSCan6RbL38sueGqld\n29vsdumS7PghwFNPwTnnwFZbwd/+lmwnSpEYqPpFKk7PnvDKK34D8//+zytVPv44ufHNvDJm4UI4\n+mjo399vriaZ8xepJBTUJRo1angwX7bMb2Tm5cGECckejpGTAyed5DdTc3P9HcOwYfD558nNQaSC\nKahLtBo08GA+axbMmAHt2iXbcgC87cGYMV7TvmoVtG7tpzHp9CWpBpRTl3g99xwMHeq57quugo4d\nk59DYaHvTF282NsOHHaYTl+SSk9tAqTyWr0a/vEPGDvWe6hffDE0aZL8PAoK4C9/8YA+YQL06ZP8\nHETKSDdKpfLKyYFTTvF8+xZbeL790kuTa/G7Tn6+39A96yw47ji/obpsWbJzEImZgrokZ/PNfYU8\nZ463+c3NhfvvTzbfXqOGHwqybBnsvrt/nHoqrFyZ3BxEYqSgLslr0cKP0rv1Vl+x7767N+5KUp06\nnopZtsw3U7Vt6z3lv/8+2XmIRExBXSrOHnt4MP/jH+HAAz0lktSReus0bAhXXulpmQULvFLm7rth\n7dpk5yESEQV1qVg1a3qDsOXLYbvtvHHYuHHw3XfJzmOnnbzdwJQpcPXV0K0bzJyZ7BxEIqCgLpXD\n737n5Ybz5sGSJb/0k0l6xdyrlx+rd+aZnnsfNAjeeivZOYhkodSgbmaTzGylmS0q9tgWZvasmS03\ns2fMrEG805Rqo1kzuO++X1bMPXrA7NnJzqFGDTjqKM+3d+oEXbt6/v3LL5Odh0gGyrJSvxXYZ73H\nzgOeDyG0Bv4JjIh6YlLNrVsxn3YaHHKIr5rfey/ZOdSr590olyyBr7/2fPv11+tAbKnUSg3qIYSZ\nwBfrPdwfuL3o89uBARHPS8RXzMce6/n2li191XzBBb71P0lNmsDNN/vu2Ecf9Ta/06apza9USpnm\n1BuHEFYChBA+BhpHNyWR9Wy6qZ+4tK5/e+vWXg65Zk2y8+jQAZ591qtlhg+Hfv18TiKVSE5Ez1Pi\nkmXs2LE/f56fn0++ToWXTOywA9x1l6dlhg6Fa67x3un9+iU3h3VtfvfeG265Bf7wB9h/f7joIth2\n2+TmIalTUFBAQUFB1s9Tpt4vZtYUmBZC6FD09VIgP4Sw0syaANNDCG028rvq/SLRCwEefNBXzLm5\nvlO1bdvk5/HVV161849/eMXMsGH+zkIkS3H3frGij3UeA44v+vw44NHyDiySFTMvNyws9JV6fn7F\nbPdv0MDPTH3ttV9KMe+4Q5uXpMKUpaRxCjALaGVm75nZH4HLgL3MbDnQt+hrkeRtsgmcfbaXH9ap\n4/3bx49Pfrv/73/vpZj33QcTJ/pxejNmJDsHEdR6V9LmX/+C887zbf/jx3spZNK900Pw4H7eeV6x\nc/nlXr0jUg5qvSsC3ixs6lTfvHTddb5xKOkVsxkcfri/e+jWDbp391y7Ni9JAhTUJZ169fKdqMOG\neaOwAQP87NIk1anjN3LXbV7KzYWbbkq+FFOqFQV1Sa8aNeCII3zF3KOHfwwZAp99luw8tt7ayx+f\nesrfQXTqBP/8Z7JzkGpDQV3Sr04dOPdcWLrUq1LatPESyKRPXurUyY/UGzMGTjwRDjrI7wGIREhB\nXaqPRo28d8vMmf7Rpo3f0EzyRr4ZDBzopZi77eY593PP9fSMSAQU1KX6ad3ae7hMnuyVKd27w6xZ\nyc6hTh2vjnnjDU8HtW7tG5iUb5csqaRRqre1a731wMiRvmq+7DI/MCNp8+b5gdirVnnL4T59kp+D\nVCoqaRTJRPFOkB07ekpk6FD473+Tnccuu8CLL8KIEV6tM3AgvP12snOQVFBQF4FfeqcvXuw3UHNz\nPTWT5M5UMzj0UL+h27mz19iPGAHffJPcHKTKU1AXKa5JE9/m/9JLXufeunXyvVzq1vUXmEWL4KOP\nfA6TJ6ufjJSJcuoiJZk5E845x1fsEybAXnslP4dXX/V8+w8/eL59992Tn4MkLtOcuoK6SGlCgIce\n8mqV5s29K2PHjsnP4d57fYdq9+4+h2bNkp2DJEo3SkXiUry2/MAD/WCM445L9sxUs192x7Zr5zdW\nR42Cb79Nbg5SJSioi5RVrVrw5z97D5mmTX2H6PDhyTbqqlcPRo+GhQv9aL/cXLjnHp2XKj9TUBcp\nr/r1Ydw43zj0xRfQqhVcdRX8+GNyc9h+e7j7bk/JTJgAvXvD/PnJjS+VloK6SKa23RZuvtn7uUyf\n/suqOckqlZ49/UbqscfCvvvCKack37BMKhUFdZFstW0Ljz0Gt93mK/auXZPtwlizJpx0kte3b7KJ\nz+e662D16uTmIJWGql9EohQCPPCAbxpq3dqrVNq3T3YOS5b4IdgrV8I118CeeyY7vkRC1S8ilcG6\nXaGFhbDPPtC3L5xwAvznP8nNoV07eO45uPBCGDwYDjkE3n03ufGlQimoi8Rhk018tbxihe9S3Xln\nOP98+OqrZMY3g4MP9heXDh287cDYsfDdd8mMLxVGQV0kTg0a+AHYCxf6lv+kK2Xq1oULLvDKmKVL\nvYf8Aw+oBDLFlFMXSdLixb5iX7TI0yNHH+03OpMyY4Yf6bfllnDttcnn+6XMlFMXqQry8rxS5q67\n/NzSjh1h2rTkVs59+njv9kMO8Xz/6afD558nM7YkQkFdpCL06uWdIMeP90qZ3XeHl19OZuycHDjt\nNE/HhOApmRtv1KlLKaH0i0hFW7PGV+6jR/sN1fHjfUWflIULPSXz1Vd+hmuvXsmNLRul9ItIVVWz\npjcIW74c9tjD68r/+MfkGobtvLPvij3vPG8aduyx8PHHyYwtkVNQF6ks6tSBs8/2MsjttvOGYcOG\nJXO0nhkcfriXQDZp4jdQr75au1KrIAV1kcqmQQO4+GKvlPn+e9+ZesklybTZ/d3v/Bi/F1+Exx/3\nF5YXX4x/XIlMVkHdzM40szeKPoZENSkRAbbZxo/WmzPHO0K2bOk3NH/6Kf6x27TxXamjR3vZ5VFH\nwYcfxj+uZC3joG5m7YDBQBegI3CAmTWPamIiUqRFC2+xO20aPPywN+y67774u0GaeeljYSHsuKPv\nTL3iimReVCRj2azU2wBzQwg/hhDWAC8CB0czLRH5jV12gWefhb//Hf72N+8G+dxz8Y+72WZw6aVe\ncvnss15bX1AQ/7iSkYxLGs0sF3gE6A78CDwPvBpCOHO9n1NJo0jUQoCpU2HkSF9FX3YZdOmSzLgP\nP+w3dHv08BeX7baLf9xqKNOSxpxMBwwhLDOzvwLPAauA+cAGdy+MHTv258/z8/PJz8/PdFgRgV9S\nIwMGwOTJ0L8/dOsGF13k6Zk4xz34YD+ndfx4L4c891w46yyoXTu+cauBgoICCiJ4BxTZ5iMzuwR4\nP4Tw9/Ue10pdJG7ffw833OCVK/vu6x0Zf//7+MddscK7Uf77375xqW/f+MesJipk85GZNSr6d0fg\nIGBKNs8nIhmqWxf+8hf41788mHfp4odkf/RRvOO2bAlPPOHpnxNP9F7y778f75hSomzr1B80s8XA\no8BpIYSvI5iTiGSqfn1fpS9f7oE+Lw+GD493A5OZp38KC70UslMnz7WrSqZCqPeLSJr95z++kWnq\nVE+TnHWWbzCK04oV/i7h44+9rr5nz3jHSyn1fhGR39p+ey+BnDPHV+8tWvghHT/8EN+YLVvCM8/A\nqFGejhk8GD77LL7x5FcU1EWqgxYtvBPk88/7QRktW3o/97hSJOvOal261N8ZtGsHkybFv2FKlH4R\nqZbmzvUa93fegXHjvJlXjRjXePPnwymneC/3G2/03alSIqVfRKTsdtvNV+033eTH2sV9AlOnTjBr\nlrf17dfPK3VWrYpnrGpOQV2kOuvbF2bP9pup55/vu0SnT49nrJo14U9/8u6Tn37qm6QeekiHYEdM\n6RcRcWvWeKOw0aOhWTPfndq9e3zjzZgBp57qdfXXXQfN1Q+wOKVfRCQ7NWvCkUf6zc3DDvM8+377\nwauvxjNenz6wYIGfz9q1q/eM//HHeMaqRhTUReTXatWCk06CN9+EAw6Agw6CAw/0ABy12rX9GL3X\nXvObtx076lCOLCmoi8iGbbIJnHaatx7o29d7ygwc6DnxqDVrBo8+6qv1I4/0F5Uvvoh+nGpAQV1E\nSlanju9Gfestz7H37eupmWXLoh1nXQfIJUt8Bd+unef4dU+uXHSjVETKZ9Uqv7F55ZW+eh892jc3\nRW32bDj5ZNhhBz/Wr1mz6MeoxHSjVESSsdlmMGKEr9xbtPA+7iec4O13o9S9O8yb5zdSu3Txo/RW\nr452jBRSUBeRzNSv76v0FSu8x0yXLl6HHmXr3dq1/QVk9mx48kmvkpk3L7rnTyEFdRHJzhZbeKuB\nN9+ELbf0CpbTT4cPP4xujJYtfQfsWWd5meXQodqRuhEK6iISjYYN/YDqpUu9ciYvz4NwVAd1mHmb\ngSVLvD98u3Z+QIf8ioK6iESrcWPPfy9Z4oG4XTsYMsR7u0dhq63g9tu96+OQIXD00WrtW4yCuojE\nY5ttvHd7YaHnxjt08Lr3996L5vn79YNFi/xFpH17uP9+lT+ikkYRSconn3gZ5C23wKBBfgM0qjLF\n2bP9MI7cXD+Ae5ttonneCqSSRhGp3Bo39gOqly/3FMouu3ggfuut7J+7e3d4/XXv/LjzznDbbdV2\n1a6VuohUjM8/h2uu8ZX1AQf4oR0tW2b/vPPne9381lvDzTfDjjtm/5wVQCt1EalattwSLrzQe8s0\nb+693I85Jvv2A506wSuvQO/e/m5g4sRqdYyeVuoiUjl8/TVcfz1cfbXfBB050itnsrF0qa/a69WD\nyZOhadNo5poArdRFpGqrX99PX3rrLc+L9+0LAwZ4S95MtWkDM2fCXnv5jtdbb019rl0rdRGpnL77\nzlfXEyZ4j5kRIzzQW7kXr27RIt+81LSpn83apEm0842YVuoiki716nm7gX/9y4PxkCF+YPbDD2eW\nI+/QwXPteXneymDq1OjnXAlopS4iVcPatX6QxqWXet+X4cP9QI1atcr/XHPmwHHHeUrm+uu9f00l\no5W6iKRbjRp+tN7cuXDttXDHHV4Cef31nqopj27dvPRxXQOymTPjmXMF0EpdRKquV17xlfvs2XDG\nGXDKKd5YrDyeeAJOPNHbBo8aBTk58cy1nCpkpW5mZ5vZYjNbZGZ3m1ntbJ5PRKRcunb1HPsLL3jV\nTMuW8Oc/e4/3stp/f9+N+vLLkJ8P774b23STkHFQN7NtgTOAziGEDkAOcHhUExMRKbN27bxSprDQ\n8+M9eng55Esvla2EcZtt4Jln4MADYddd/YWiiso2p14T2NTMcoB6QIRd8UVEyqlJE7j4Yl9t77OP\n95bp2hXuvbf0o/Bq1IBzz4XHH/c+8CNGwJo1ycw7Qlnl1M1sCHAJ8B3wbAjhmA38jHLqIlIx1q71\nIH3FFfDOO553P+EEv0Fakk8/hSOO8Jr4KVOgUaNEpltc4jl1M9sc6A80BbYFNjOzIzN9PhGRyNWo\n4SmVGTPgwQdhwQLYaScP7CWdddqoETz9tPeO6drV2w1UEdnc5u0HvB1C+BzAzB4CegBT1v/BsWPH\n/vx5fn4++fn5WQwrIpKBLl3grrt8FT55svd0b9zYD+449FCoW/fXP5+T462C27b1G6j33w99+sQ2\nvYKCAgoKCrJ+nozTL2bWFZgE7Ar8CNwKvBpCuGG9n1P6RUQqnzVrfDU+caKXRh51lO9c7dTpt60I\nXnjB0zH33OOtChKQafol25z6GLzi5SdgPnBiCOGn9X5GQV1EKre33/ZzT++801fsxxzjQX6HHX75\nmRdfhIEDoaAg++6RZVAhQb1MAyioi0hVEYLXq995px+wAXDBBd4KuHNnP1Fp2jQvf4yZgrqISJS+\n+caD+3vvwfPPew38//7nPWMmTYp9eAV1EZE4hQA//PDbG6oxUVAXEUkRdWkUEREFdRGRNFFQFxFJ\nEQV1EZEUUVAXEUkRBXURkRRRUBcRSREFdRGRFFFQFxFJEQV1EZEUUVAXEUkRBXURkRRRUBcRSREF\ndRGRFFFQFxFJEQV1EZEUUVAXEUkRBXURkRRRUBcRSREFdRGRFFFQFxFJEQV1EZEUUVAXEUkRBXUR\nkRRRUBcRSZGMg7qZtTKz+Wb2etG/X5nZkCgnJyIi5ZNxUA8hvBlC6BRC6AzsAnwLPBzZzKqIgoKC\nip5CrNJ8fWm+NtD1VVdRpV/6AW+FEN6P6PmqjLT/YaX5+tJ8baDrq66iCuqHAfdE9FwiIpKhrIO6\nmdUCDgQeyH46IiKSDQshZPcEZgcCp4UQ/rCR72c3gIhINRVCsPL+Tk4E4x5BCamXTCYlIiKZyWql\nbmb1gHeWPPGKAAAD6UlEQVSB5iGEbyKblYiIZCTr9IuIiFQesewoNbNBZrbYzNaYWef1vtfBzGYV\nfX+hmdWOYw5xKun6ir6/o5l9Y2ZDK2J+2drY9ZlZPzN7rej/7VUz26Mi55mpUv4+R5jZCjNbamZ7\nV9Qco2JmO5vZ7KINgq+YWZeKnlPUzOyMov+vN8zssoqeT9TMbJiZrTWzLcvy81Hk1DfkDeAg4Kbi\nD5pZTeBO4KgQwmIz2wL4KaY5xGmD11fMFcCTyU0nchu7vk+BA0IIH5tZO+AZYPukJxeBjf19tgEO\nBdrg1/W8mbUMVfvt7OXAmBDCs2a2LzABqJIvxhtiZvnA/wHtQwirzWyrCp5SpMxse2AvPM1dJrEE\n9RDC8qIJrX+TdG9gYQhhcdHPfRHH+HEr4fows/7A2/gO2yppY9cXQlhY7PMlZlbHzGqFEKrUC3MJ\n/3/9gXtDCKuBd8xsBdAVmJvwFKO0FmhQ9PnmwAcVOJc4nApcVvR/RgjhswqeT9SuAs4BHivrLyTd\n0KsVgJk9XfQ2/pyEx4+VmW0KnAtcCKS66sfMBgGvV7WAXortgOK7oj8oeqwqOxv4m5m9h6/aR1Tw\nfKLWCuhtZnPMbHqa0ktF5eLvhxDeKM/vZbxSN7PngK2LPwQEYGQIYVoJ4/UEugA/AC+Y2WshhOmZ\nziMuGV7fWOCqEMJ3RYvAShvYM7y+db/bDrgUf1tYKWVzfVVNSdeKt/A4M4TwSNEL8WQq8f/bhpRw\nfaPwmLJFCKGbme0K3A80T36WmSnl2s7n1/9XZYonGQf1EEImfxj/AV5cl3YxsyeBzkClC+oZXt9u\nwEAzuxzYAlhjZt+HECZGO7vsZXh963J8DwHHhBDeiXRSEcrw+j4Adij29fZUgXRFSddqZneGEM4s\n+rmpZjYpuZlFo5TrOwX/eySE8GrRDcWGIYT/JjbBLGzs2swsD2gGLCxKE24PzDOzriGET0p6ziTS\nL8VfXZ4B2hflYnOAPkBhAnOI08/XF0LoHUJoHkJoDlwNjK+MAb2cfr4+M2sAPA4MDyHMqbgpRar4\n3+djwOFmVtvMfg+0AF6pmGlF5gMz6wNgZn2BNyt4PlF7BNgTvB04UKuqBPSShBAWhxCaFMWT3+ML\n4k6lBXSIr6RxgJm9D3QDHjezp4om+iVwJfAa8DrwWgjhqTjmEKeNXV9alHB9pwM7AaOL9dKvctUG\nJfx9FuJv3wvx6qXTqnjlC8BJwBVmNh+4GDi5gucTtVuB5mb2BjAFOLaC5xOXQBnTL9p8JCKSIjrO\nTkQkRRTURURSREFdRCRFFNRFRFJEQV1EJEUU1EVEUkRBXUQkRRTURURS5P8BQ9ybkgJkrcQAAAAA\nSUVORK5CYII=\n", "text/plain": "<matplotlib.figure.Figure at 0x7fd3c5e3f0b8>"}, "metadata": {}}], "metadata": {"collapsed": false}}, {"execution_count": 48, "cell_type": "code", "source": "zb = find_z(ay, by, cy, y, az, bz, cz)[0]\nze = find_z(ay, by, cy, y, az, bz, cz)[1]\nplot(y, zb, 'r', y, ze, 'g')", "outputs": [{"execution_count": 48, "output_type": "execute_result", "data": {"text/plain": "[<matplotlib.lines.Line2D at 0x7fd3c6032a90>,\n <matplotlib.lines.Line2D at 0x7fd3c6032da0>]"}, "metadata": {}}, {"output_type": "display_data", "data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xmcj+X+x/HXNcaSJcTILpNdlCKRNKVNWSMcJyqhHadO\nncphGhT1E9rPKUsoZM1Stg6TFCLZJVvWElmyNMPMXL8/rmEmMcN3vvc9M995Px8Pj+7rO9/7/nzu\nZh6fueb6Xvd1GWstIiISGsIyOwEREQkeFXURkRCioi4iEkJU1EVEQoiKuohICFFRFxEJIekWdWPM\nCGPMPmPMmlSv9TPGrDbGfG+MmWOMKeltmiIiciFMevPUjTGNgGPAGGtt7eTXClprjyUfPwXUsNY+\n5nWyIiKStnR76tbaxcChs147lqpZAEgKcl4iIhKA8EBPNMYMADoDh4FbgpaRiIgELOAPSq21/7bW\nlgc+Bp4KXkoiIhKogHvqqYwDPgdeOtcXjTFaXEZEJADWWnOx51xoT90k/3MNYyql+lorYGM6iYXs\nv+jo6EzPQfene9P9hd6/QKXbUzfGjAOigGLGmJ1ANHCPMaYqkAjsAB4NOAMREQmadIu6tbbjOV4e\n5UEuIiKSQXqiNIOioqIyOwVPhfL9hfK9ge4vp0r34aMMBzDGeh1DRCTUGGOwHn5QKiIi2YCKuohI\nCFFRFxEJISrqIiIhREVdRCSEqKiLiIQQFXURkRCioi4iEkJU1EVEQoiKuohICFFRFxEJISrqIiIh\nREVdRCSEqKiLiIQQFXURkRCioi4iEkJU1EVEQoiKuohICFFRFxEJISrqIiIhREVdRCSEqKiLiIQQ\nFXURkRCioi4iEkJU1EVEQki2L+pJNonCgwpjYgwmxtBxSke2HNyS2WmJiGSKdIu6MWaEMWafMWZN\nqtdeM8ZsNMasMsZMMcZc6m2a57d383f8Hv87TSrcQt5ceRm/bjyV36qMiTFUfKMi3+39LrNSExHx\nnbHWpv0GYxoBx4Ax1traya/dBiyw1iYZYwYB1lr7wnnOt+nFyDBj3H+T4xyNP0qfhX14Y9kbZ95y\nU/mbmNJuChEFIrzNRUQkCIwxWGvNRZ93IQXXGFMBmHm6qJ/1tVZAG2ttp/Oc631R37QJqlWDsDBI\nSEgp8sDhuMM8NfspPlrz0ZnXRrcaTeerO3ubk4hIBmRmUZ8BTLDWjjvPud4XdYDVq+Gaa9xxfDzk\nyfOXtyzZtYQbR96IxeXTsFxDZnSYQbH8xbzPT0TkImRKUTfG9Aautda2SeNcf4o6wE8/QcWK7njn\nTihX7pxvO5l4kl5zevHeivfOvLb04aXUL1vfhyRFRNLne1E3xjwIdANutdbGp3GujY6OPtOOiooi\nKirqYvO8cCdOQIEC7nj6dGjRIs23z9kyh6YfNz3THnbnMHrU74ExF/3/UkQkYLGxscTGxp5px8TE\neFrUr8AV9VrJ7buA14HG1trf0jnXv556ag0awNKlcN99MHFium/fdWQXN426iR1HdgDQomoLJrSZ\nwCW5L/E6UxGRv/Csp26MGQdEAcWAfUA08CKQBzhd0Jdaax8/z/mZU9QBPvgAund3xwcPQtGi6Z5y\nKvEUj8x6hFGrRgFQulBplndbTulCpb3MVETkTzwdfsmITC3qALt2Qfny7njaNGjV6oJPff+793lk\n1iNn2qseWcXVJa8OdoYiIn+hop4Wa6F+fVi+HG68Eb766k/THtPzxbYvuH3s7Wfan3f8nKaVm6Zx\nhohIxqioX4iJE6F9e3e8caOb234RNuzfQM13a55pv930bZ64/olgZigiAqioX7gjR6BIEXf82GPw\n7rsXfYl9x/bRcGRDth3aBsCzDZ/l1dte1YwZEQkaFfWL1a8fnJ5quXcvlCp10Zc4ceoELSe05Itt\nXwBwf+37GdVyFOFh4cHMVERyIBX1QOzenfKA0quvwnPPBXSZhKQEus7oyujVowFoVqUZk++bTN7w\nvMHKVERyGBX1jHj4YRg50h0fPgyFCwd0GWstz8x7hqFLhwJwY7kbmXv/XArkKRCsTEUkh1BRz6j1\n6+Gqq9zx0KHQq1fAl7LW0n9Rf6Jj3fDOVSWu4quHvqJIviLByFREcgAV9WCwFtq1g8mTXfvAASiW\nscW+3lj6Br3mul8QZQqVYeUjKylRoERGMxWREKeiHkwbN0KNGu54wADo3TvDlxy9ajQPTn8QgMJ5\nC7PhiQ16SlVEzktF3QsPPQQffuiOf/4ZSpbM8CUnrZ9Eu8ntALgk/BI2PbmJcoXPvZqkiORcKupe\n2bYNrrzSHT//PAwcGJTLTv9hOq0+cUsW5DK52NZzG+ULlw/KtUUk+1NR91qPHvDWW+54x46U9WQy\naPbm2dw97m4AwsPC2dZjm3ruIqKi7ovUi4N16gRjxgTt0qmLe55cedjaYytlLy0btOuLSPaiou6n\n/v2hb193vGwZXH990C792Y+f0Wx8MwDy5srLlh5bVNxFciAVdb+lXkOmRg23R2p48JYHmPXjLJqP\nbw5AvvB8bHlqC2UuLRO064tI1qainlmmTIG2bd3x+PHQoUNQLz9z00xaTHBb8hXIXYDNT22mVKGL\nX6dGRLIXFfXMlJgI9erB99+79qFDKb34IJmxaQYtJ7QEoHj+4mx4fAMRBSKCGkNEso5Ai3qYF8nk\nOLlywcqVsGKFaxctmrICZJC0qNoCG22ZdN8kDpw4QInBJSg/tDwH/zgY1Dgikr2pp+6Fbt1g+HB3\nvHkzVKoU9BAfrfmITtM6AVD5ssos77acwvkCW4hMRLIeDb9kNXv3QpnkDzbvuQdmzryoLfQu1IiV\nI+g6sysAdUvX5csHvyR/7vxBjyMi/lJRz6reess9uAQwa5Yr8B5IvXDY7ZG3M6vjLPLkyuNJLBHx\nnop6VhYXB1Wrws6dru3BB6mn9fuy35klf++rcR/j24wnV1guT2KJiHf0QWlWli+fW1pg6VLXLloU\nnn7ak1B9b+5LUt8knr7haSZtmER4/3C6z+xOjv/FKpJDqKeeGVKvI7N8OdSt60mYJJtE1xldGbVq\nFADPNXyOQbcN0gbZItmAhl+ym4MHUzbgqFjRreGe15s9TROSEmg3qR3TfpgGwMAmA3m+0fOexBKR\n4FBRz66mT4dWbgle3ngj5UNVD8QnxHPnR3fy5Y4vARjefDgPX/uwZ/FEJHAq6tmZtdCsGXz+uWtv\n3QqRkZ6FOxp/lPrD67PxwEYAprabSuvqrT2LJyIXT0U9FPz0kxuKAbj9dpgzB8K8+yx7//H9VHun\n2pmnUhc+sJCoK6I8iyciF05FPZS88w48+aQ7njIF7r3X03A7j+ykwrAKZ9oru6+kTqk6nsYUkbR5\nVtSNMSOAZsA+a23t5NfaAi8B1YF61tqVaZyvoh6IkyfhqqvcMgMQtD1S07Jh/wZqvlvzTHtbj21U\nLFrR05gicm5ezlMfBdx51mtrgdbAlxcbUC5Qnjzw449uoTCAUqXcEr8e/oKsEVEDG235pss3AES+\nGUnhQYXZf3y/ZzFFJLjSLerW2sXAobNe22St3QxowrPX6tRxhfzll91QTFgYTJ3qacgG5Rpgoy3T\nO0zn9/jfKTG4BLXeq8Xxk8c9jSsiGacnSrOLF1+EP/5ws2LatHGLg/38s6chTy/3+99m/2Xdr+so\nOLAgd398NwlJCZ7GFZHABW//tTS89NJLZ46joqKIioryI2zoyZfPTXdcvRquuQZKl3Zz3KdO9WQF\nyNO6X9ed7td1J3phNP0W9SN3/9x0v7Y7/2n2Hz2dKhIksbGxxMbGZvg6FzT7xRhTAZh5+oPSVK8v\nBJ7RB6WZZNAgeOEFdzxxItx3n+chrbV0ndGVkatGAjDglgH0btzb87giOY2nUxqNMVfginqts15f\nCPzTWvtdGueqqHspPh5q13YfqgLs3p2yjruHTiWe4p5x9zB/23wAxrYey/217/c8rkhO4eWUxnFA\nFFAM2AdE4z44fQsoDhwGVllrm57nfBV1P6xd64o7uKdTp0/39MGl047GH6XOf+uw9dBWAP7X+X/c\nWvFWz+OKhDo9fCTO4MHw7LPueMIEaN/el7B7j+6lzJCUvxDWPraWq0pc5UtskVCkoi4pTp50H6Ru\ndGu7sGsXlC3rS+h1v66j1nspo3R7nt5D6UKlfYktEkpU1OWv1q93T6WCL2vJpLZg+wKajGkCQMUi\nFVn16CouzXupL7FFQoF2PpK/qlnTPbg0dCjMnw+5csGIEb6EvrXirdhoy9jWY9l+eDuFBxXmro/u\n0hx3EY+pp55TnDoFN97odloCNzRTrZpv4QcsGkCfhX0A6FW/F0PuHKI57iJp0PCLXJjt21PWaq9c\nGdascQ81+SDJJtF5Wmc+XvsxAO/d8x6P1n3Ul9gi2Y2KulyciRNTZsb07g0DBvgW+o9Tf9BoVCNW\n/uyeWZvz9zncWensNeNEcjYVdbl41kLnzvDRR64dGws33+xb+F+P/0qFYRWIS4gDYM2ja6h1ea10\nzhLJGVTUJXCHDsHll7txd4ADB1I2xfbBxv0bqfFujTPtvU/vpVShUr7FF8mKNPtFAle0qJvb/tVX\nrl28uBua8emXcfWI6thoyxedvgCg9JDS1Hy3JidOnfAlvkgoUVGXFI0auUL+0ktuzD0sDMaN8y18\nk8gm2GjLiBYj2LB/AwVeKUDHKR1Jskm+5SCS3Wn4Rc4tPh6uu849wASwZQtceaVv4a21PDf/OQYv\nGQzAwCYDeb7R877FF8lsGlMXb/z4I1St6o6vvhq+/dZtteeThKQE7hl3D/O2zgNgaruptK7e2rf4\nIplFRV28NWYMPPCAO375ZbcTk48Oxx2m2tvV2Hd8HwAru6+kTqk6vuYg4icVdfFeUpLb/HraNNde\nvNg9peqjrQe3UumtSmfamikjoUpFXfxz4ABERKS09+93M2Z8tGjHIm7+0M2prxlRk+XdlnNJ7kt8\nzUHES5rSKP4pXtzNklm0yLUjIqB5c9eT90njCo2x0ZaRLUayfv968r+Snwc/fRB1ICSnU1GXwN10\nkyvuAwfCrFluFch33vE1hYfqPERS3yR61e/F6NWjCesXxpvL3vQ1B5GsRMMvEhwJCdCkSUrv/dtv\noV49X1M4mXiSJmOasHjnYgDmd5rPbZG3+ZqDSLBoTF2yhr17Uza+zpfPtYsW9TWF/cf3U25oOeIT\n4wH48ckfqVyssq85iGSUxtQlayhd2g3JfPEFxMXBZZdBmza+jrdHFIgg7t9xrHl0DQBV3q5CycEl\nORJ3xLccRDKLirp4o0mTlCUHpk514+3vv+9rCrUur4WNtkxrP419x/dR5NUiNBvXjMSkRF/zEPGT\nhl/Ee6dOQePGsHSpa69cCXX8f3Ao9e5L/aL60efmPr7nIHKhNKYuWd/u3VCunDsuVAh27oQiRXxN\nITEpkXsn3suMTTMAmNFhBs2rNvc1B5ELoaIu2cfcuXDXXe64Qwe3EqTP+5UejjtM5bcqc+DEAQA2\nPrGRasX927NVJD0q6pL9/Pvfbh0ZgJEj4aGHfE/hhwM/UP2d6gAUz1+cLU9toXC+wr7nIXI2FXXJ\nnuLjoUED+P571161yq0G6bOZm2bSYkILAJpXac609tPIFZbL9zxETlNRl+xtxw644gp3nEnj7QD9\nv+xP39i+gNZwl8yloi6hIfV4e+vWMHmy24HJR4lJiTQf35zZW2YDMPvvs7mr0l2+5iDiWVE3xowA\nmgH7rLW1k18rCnwCVAB+AtpZa8/5ZIeKugTkpZcgJsYdDx0KvXr5nsKhPw5xxRtX8Hv87wBs7bGV\nyKKRvuchOZOXRb0RcAwYk6qovwr8Zq19zRjzL6Cotfacf6eqqEvAEhLgjjtg4ULXjo2Fm2/2PY21\n+9ZS+z+1Abiy6JWseWwN+XPn9z0PyVk8HX4xxlQAZqYq6j8AN1tr9xljSgKx1tpzzgdTUZcM+/VX\nKFnSPaEKsGsXlC3rexqfrPuEDlM6ANDlmi4MbzEc4/NUTMk5/F77pYS1dh+AtfYXoESA1xFJX4kS\nbu2Yb7917XLl3BOp8fG+ptH+qvYk9U2iZ/2ejFw1krB+YYz6fpSvOYikJ9Ce+kFr7WWpvv6btbbY\nec610dHRZ9pRUVFERUVlNG/JyUaOhIcfdsc9e8KwYb6nEJ8Qzw0jbmDVL6sA+Lbrt9Qr4+9SwxJa\nYmNjiY2NPdOOiYnxdfhlIxCVavhlobW2+nnO1fCLBJ+10K0bjBjh2uPHu6dTfbb7992UG+qWPjAY\n9v1zHxEFItI5SyR9Xg+/mOR/p80AHkw+fgCYfrGBRTLEGBg+HI4dgyuvhL/9zb22Zo2vaZS9tCw2\n2hL7QCwWS4nBJbhtzG0kJCX4mofIaRcy+2UcEAUUA/YB0cCnwCSgHLADN6Xx8HnOV09dvLdtmyvu\n4B5e2rHD9805AIYsGcIz854B9PCSZIwePhIBmDMHmjZ1x82awaefurXcfZSYlEjLCS35bPNnACzo\nvIBbKt7iaw6S/amoi6Q2YAD0SV4v/bXX4NlnfU/htxO/UXpIaU4mngRg9z92U+bSMr7nIdmTirrI\n2RIS4O67Yf58154/H27zfyPqFXtXUO8DNzOmXul6LO6ymDy58vieh2QvKuoi5/Pbb25u+x9/uPaW\nLSnj7z56/7v3eWTWIwC82OhFXm7ysu85SPahoi6SntWr4Zpr3HH58rBunftQ1UfWWv4+9e+MXzce\ngM86fsbdle/2NQfJHlTURS7UxInQvr07vvdemDTJ95Ugj8QdIfLNSA7+cRCAHb12UL5weV9zkKzN\n72UCRLKvdu3cw0u9e8PUqW52zMCBvqZQOF9hfnvuN1Y94p5IrTCsAnXfr3vmQ1WRQKmnLjlbQgI0\nb+6mQgLMmOHaPhuxcgRdZ3YFoE/jPvS7pZ/vOUjWouEXkYw4fBgqV4YDbiNq1q2DmjV9TcFaS8ep\nHZmwbgIA8+6fx+1X3u5rDpJ1qKiLBMOmTVAteRXp/PndtnrFzrlWnWcOxx2m3NByHDt5DND89pxK\nY+oiwVC1qhtvnzsXTpyA4sXdxhynTvmWQpF8RTj6wlFWdFsBQNmhZWk4oiGnEv3LQbIvFXWRc7nj\nDlfchwyBRYsgTx745z99TeG60tdhoy3v3fMeS3YvIc+APAxaPMjXHCT70fCLSHqshQcegLFjXXvM\nGOjUydcUkmwSrT9pzYxNMwD4usvXNCzX0NccxF8aUxfx2okTULcubNzo2kuWwA03+JrC/uP7KTHY\nbTQWZsL49Z+/Uiy/v2P+4g8VdRG/7N7tlh04LRP2TF20YxE3f+g24W5TvQ2T7puk/VJDjD4oFfFL\n2bJuSOabb1y7XDk3Y+b4cd9SaFyhMTbaEhMVw5SNUwjrF8aHqz70Lb5kXeqpi2TU2LHQubM7btfO\nba3n47IDJxNP0mBEA1b+vBKA9Y+vp0ZEDd/iizc0/CKS2V54AQYlz06JiYG+fX0Nv/3QdiLfjAQg\nsmgk6x5bxyW5L/E1BwkeFXWRrCAhAVq1gs/crkd8/DF07OhrClM2TKHtpLYAvNDoBV5p8oqv8SU4\nVNRFspKjR90yv9u2ufZXX0GjRr6Ft9by0PSHGL16NAALH1hI1BVRvsWXjFNRF8mKzp4ps3kzVKrk\nW/hDfxyi1OuliE+MB2D/s/spnr+4b/ElcJr9IpIVnZ4p8913rl25MhQtCgcP+hK+6CVFift3HN90\ncTN1Iv4vgrYT26KOVuhSURfxw7XXuuI+Y4ZbEbJYMWjYEOLjfQnfoFwDbLSl/y39z0yB/GjNR77E\nFn9p+EUkM7zxBvTq5Y67dIHhw8Gnh4dOJp6k3gf1WLNvDQBbe2wlsmikL7HlwmlMXSS7sRaefBLe\nfde1X3sNnn3Wt/Cbf9tMlberAHB9mev5usvXhIeF+xZf0qaiLpJdnToFd94JCxe69pQpbu9Un4z6\nfhRdZnQBYNidw+h5Q0/fYsv5qaiLZHeHD0P16vDLL669bBlcf70voZNsEk0/bsq8rfMAWPPoGmpd\nXsuX2HJuKuoioWL7dohMNca9bRtUrOhL6D2/76HsULc4WcUiFVn/+Ho9lZpJNKVRJFRUrOjG25cs\nce3ISLcDkw/TIMtcWgYbbZnabirbD28n/yv5eSn2Jc/jSvBkqKdujOkJdE1ufmCtffMc71FPXSQj\npk1LGWOvW9ftxHSJ971nay33T7ufcWvHAbC823Lqlq7reVxxfB9+McbUBMYD9YAEYDbwqLV221nv\nU1EXCYa334annnLHrVvDpEmQK5fnYVNvzFGyYEm29thK/tz5PY+b02XG8Et1YJm1Nt5amwgsAvz7\nyF4kp3nySUhKgn/9y/Xew8OhZ083VOOhiAIR2GjLrL/N4pdjv1DglQL0WdDH05gSuIz01KsBnwIN\ngHjgC2C5tbbnWe9TT10k2JKS3BruH3/s2j7NcbfW8uD0BxmzegwAy7ou4/oy/szQyWkyZfaLMeYh\n4AngGLAeiLfWPn3We2x0dPSZdlRUFFFRUQHHFJFU4uPhjjvcODu4DTvuv9/zsAdOHCDi/yIAiMgf\nwfae2ymQp4DncUNZbGwssbGxZ9oxMTGZO6XRGPMysMta+5+zXldPXcRrR45AnTpuOiTA3Lmu2Hts\n9ubZ3D3ubgB639SbAbcO8DxmTpFZPfUIa+1+Y0x5YA5wg7X297Peo6Iu4peff4by5d1mHeBWh7z2\nWk9Dnr12+4puK7iu9HWexswJMquoLwIuA04B/7DWxp7jPSrqIn7btMlthn3a1q1/fqDJA6lnyZQv\nXJ5NT24iX3g+T2OGMj1RKiJ/tXQpNGjgji+9FLZsgYgIT0NO/2E6rT5pBcCrt73Kczc+52m8UKWi\nLiLnN2sWNG/ujqtXh2+/hYIFPQtnreXeiffy6Q+fArDh8Q1Uj6juWbxQpKIuIukbPhy6dXPHUVEw\nZw7kzetZuNRrydQvU5+vu3xNrjDvH5gKBVr7RUTS17Wre1jplVcgNhby5YM2bVI+WA2y02vJfNjy\nQ5btWUZ4/3DtuOQx9dRFcipr4bnnYPBg137oIdeTD/Omr5eQlMANw2/gu5/dfq0/P/MzJQuW9CRW\nKNDwi4gExlp45BH44APX/sc/4PXXPdteb92v66j1nlur/f7a9zOm1RiMT1v5ZScq6iKSMYmJ8Le/\nuYXCAPr1gz7erfHS78t+RMe6p83/1/l/3FrxVs9iZUcq6iISHCdPQrNmMH++aw8b5hYO88Cxk8co\nM6QMv8f/TsE8BfnlmV+03EAyfVAqIsGRJw/MmwfHj7v123v1ckMxH34Y9FAF8xTkyPNHWNB5AcdO\nHqPgwIK8uvjVoMfJSdRTF5G0HTniivuWLa49ebKbMRNk1lo6TOnAxPUTAdjaYyuRRb19CjYr0/CL\niHhr/36oUQMOHHDtOXPgzjuDHmbnkZ1UGFYBgGZVmjGjw4wc+UGqirqI+GPPHrjiipS57V9+CY0b\nBz3Mm8vepOccN5a/8IGFRF0RFfQYWZmKuoj4a+tWqFQppb1oEdx0U1BDnDh1gtKvl+ZI/BFKFyrN\nth7byBvu3ROwWYk+KBURf115pZvj/sMPrt24sftAdfHioIXInzs/h58/zOy/z2bv0b3kezkfw1cO\nD9r1Q5F66iISHD/84BYLO23xYrjxxqBdPskmcfvY21mwfQEA+/65jxIFSgTt+lmNhl9EJGvwuLhv\n2L+Bmu/WBKBn/Z4Mu2tY0K6dlaioi0jWsnGjmy1z2tdfQ8OGQbv8s/OeZfASt27ND0/8QNXiVYN2\n7axARV1EsqYNG6BmzZT2N9+kbNyRQb8e/5XLB18OQMuqLZnWflrITH9UUReRrM3D4v7WsrfoMacH\nAN91/45rS3m7L6sfVNRFJHtYvx6uuiqlHaQx96PxR7l00KUA3FrxVr7o9EW27rVrSqOIZA81a7qp\nkGvXunajRm4q5Ny5GbpsobyFsNGWkS1GsmD7AsL6hbFk15IgJJy9qKcuIpnrxx+haqoPOSdNgrZt\nM3TJE6dOUOy1YsQlxFG3dF2WdV1GmMlefVj11EUke6pSxfXcd+6EQoXgvvtcz33kyIAvmT93fv7o\n/QcT2kxgxd4V5OqXi0U7FgUx6axLPXURyVr273erQu7c6dpDhrjdmAIUnxBPqddLcSjuEPXL1GfJ\nw0uyxVi7euoiEhoiImDHDjh8GGrXhqefdj33vn1dj/4i5Q3Py8F/HWRCmwks27OMsH5hrNm3xoPE\nswb11EUkaztxAu65B2JjXfuJJ+DNNwPaIPvYyWMUGlgIcPujjm09NoiJBpemNIpIaDt5Etq3h08/\nde2mTWHKFLjkkou+1JAlQ3hm3jMA7Oi1g/KFywcz06BQUReRnCEpyY2xv/mma0dGugeZLr/8oi6z\n79g+Sr5eEoC+jfsSc0tMsDPNEBV1Ecl53n4bnnoqpb1+/Z/Xm7kAPWb34K1v3wIgrndcllmvPVOK\nujHmH8DDQBKwFnjIWnvyrPeoqIuItz7/3I27nzZvHtx++wWfvvqX1Vzz32sA+KbLNzQoF5zlCzLC\n99kvxpjSwFPAtdba2kA40CHQ64mIBOzuu93MmFWrXPuOO9yMmdNDNOm4uuTVnOpzikJ5CtFwZEO6\nz+zuYbLeyuiUxlxAAWNMOJAf2JvxlEREAnT11a6479kDFStCz56uuLdq5WbRpCE8LJzfX/id1257\njQ9WfoCJMcQlxPmUePBkdPilB/AycAKYZ63tdI73aPhFRDLHyZPQrRuMGePaBQvCihV/XpbgHLYc\n3ELltyoDsOnJTVQpVsXrTP8iM4ZfigAtgQpAaaCgMaZjoNcTEQm6PHlg9GjXe//wQzh2DKpVc733\n0aPPe1qlyyoR19v10qu+XZXJGyb7lHDGBdxTN8a0Be601nZLbncC6ltrnzzrfTY6OvpMOyoqiqio\nqIATFhHJkPXr3TIEcclDK1FRbhGx4sXP+fYmY5qwYPsC+kX1o8/NfTxLKzY2ltjTD1gBMTEx/s5+\nMcZcD4wA6gHxwChgubX2nbPep+EXEcl64uLgscdcD/60kSPhwQddTz6VF//3IgMXD+TFRi/ycpOX\nfUkvs6bvefc+AAAFtklEQVQ0RuNmvJwCvge6WmtPnfUeFXURydq+/tqt635aRIR7cjXVnqoDFg2g\nz8I+fNL2E9rVbOd5Snr4SEQkoxIToX9/iDnr6dKBA6FLF2pPupW1B9Zjo72vaSrqIiLBdOgQdO8O\nk1M+JD2eG+ZUgjYbVNQ9jSEi4rmEBDhyBIoV8yWcirqISAjRJhkiIqKiLiISSlTURURCiIq6iEgI\nUVEXEQkhKuoiIiFERV1EJISoqIuIhBAVdRGREKKiLiISQlTURURCiIq6iEgIUVEXEQkhKuoiIiFE\nRV1EJISoqIuIhBAVdRGREKKiLiISQlTURURCiIq6iEgIUVEXEQkhKuoiIiFERV1EJISoqIuIhBAV\ndRGREBJwUTfGVDHGfG+MWZn83yPGmB7BTE5ERC5OwEXdWvujtbaOtfZa4DrgODAtaJllE7GxsZmd\ngqdC+f5C+d5A95dTBWv45TZgq7V2V5Cul22E+g9WKN9fKN8b6P5yqmAV9fbA+CBdS0REApThom6M\nyQ20ACZlPB0REckIY63N2AWMaQE8bq296zxfz1gAEZEcylprLvac8CDE/RtpDL0EkpSIiAQmQz11\nY0x+YAcQaa09GrSsREQkIBkefhERkazDkydKjTFtjTHrjDGJxphrz/pabWPMN8lfX22MyeNFDl5K\n6/6Sv17eGHPUGPN0ZuSXUee7P2PMbcaYFcnft+XGmFsyM89ApfPz+YIxZrMxZqMx5o7MyjFYjDFX\nG2OWJD8g+K0xpm5m5xRsxpinkr9fa40xgzI7n2AzxjxjjEkyxlx2Ie8Pxpj6uawFWgP/Tf2iMSYX\nMBb4u7V2nTGmKHDKoxy8dM77S+V14HP/0gm6893ffqCZtfYXY0xNYC5Q1u/kguB8P5/VgXZAddx9\nfWGMqWyz95+zrwHR1tp5xpimwP8B2fKX8bkYY6KA5kAta22CMaZ4JqcUVMaYssDtuGHuC+JJUbfW\nbkpO6OwPSe8AVltr1yW/75AX8b2Wxv1hjGkJbMM9YZstne/+rLWrUx2vN8bkM8bkttZmq1/MaXz/\nWgITrLUJwE/GmM3A9cAyn1MMpiSgcPJxEWBPJubihceAQcnfM6y1BzI5n2AbCjwLzLjQE/xe0KsK\ngDFmTvKf8c/6HN9TxpgCwHNADBDSs36MMW2BldmtoKejDJD6qeg9ya9lZ/8ABhtjduJ67S9kcj7B\nVgVobIxZaoxZGErDS8nTxXdZa9dezHkB99SNMfOBy1O/BFigt7V2ZhrxbgTqAnHA/4wxK6y1CwPN\nwysB3t9LwFBr7YnkTmCWLewB3t/pc2sCA3F/FmZJGbm/7Cate8Ut4dHTWvtp8i/ikWTh79u5pHF/\n/8bVlKLW2huMMfWAiUCk/1kGJp17e5E/f68uqJ4EXNSttYH8YOwGFp0edjHGfA5cC2S5oh7g/dUH\n2hhjXgOKAonGmD+ste8GN7uMC/D+To/xTQU6WWt/CmpSQRTg/e0ByqVqlyUbDFekda/GmLHW2p7J\n75tsjBnhX2bBkc79PYr7ecRauzz5A8Vi1trffEswA853b8aYq4ArgNXJw4Rlge+MMddba39N65p+\nDL+k/u0yF6iVPBYbDtwMbPAhBy+duT9rbWNrbaS1NhIYBrySFQv6RTpzf8aYwsAs4F/W2qWZl1JQ\npf75nAF0MMbkMcZUBCoB32ZOWkGzxxhzM4AxpgnwYybnE2yfAreCWw4cyJ1dCnparLXrrLUlk+tJ\nRVyHuE56BR28m9LYyhizC7gBmGWMmZ2c6GFgCLACWAmssNbO9iIHL53v/kJFGvf3JHAl0DfVWvrZ\nbrZBGj+fG3B/vm/AzV56PJvPfAHoBrxujPkeGAB0z+R8gm0UEGmMWQuMAzpncj5esVzg8IsePhIR\nCSHazk5EJISoqIuIhBAVdRGREKKiLiISQlTURURCiIq6iEgIUVEXEQkhKuoiIiHk/wGyBThUVaHW\nvgAAAABJRU5ErkJggg==\n", "text/plain": "<matplotlib.figure.Figure at 0x7fd3c5c45748>"}, "metadata": {}}], "metadata": {"collapsed": false}}], "nbformat": 4, "metadata": {"kernelspec": {"display_name": "Python 3", "name": "python3", "language": "python"}, "language_info": {"mimetype": "text/x-python", "nbconvert_exporter": "python", "version": "3.5.1", "name": "python", "file_extension": ".py", "pygments_lexer": "ipython3", "codemirror_mode": {"version": 3, "name": "ipython"}}}}