{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.019588012009272404 -0.12974197660060582\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa7873995f8>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD4CAYAAAAZ1BptAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5Sd9V3v8fdn78kVci0hDZc0gQZqQQl2ehcsBWnNUlO0dtXTQ1MKK6Dl2GM9HujpWYpFXYhFvFRJA6Icl1U5TVtoilaIFT1qixNJQ1IKJNxCyI1cSTKZ2/6eP/Zvz+xM9s7MZO9n9jwzn9dae83z/J7nt+f3zG9mf+d3eX6PIgIzM7NaCq0ugJmZjV0OEmZmVpeDhJmZ1eUgYWZmdTlImJlZXW2tLkAznXHGGbFo0aJWF8PMLFfWr1//WkTMq3VsXAWJRYsW0dHR0epimJnliqSX6h1zd5OZmdXlIGFmZnU5SJiZWV0OEmZmVpeDhJmZ1dVwkJB0m6Ttkjak17JBxxdKOizpf9TJ/35J/ylpk6QHJLWl9FmSviHpe5I2S7qu0bKamdnINKslcXdELE2vRwYd+33g72plklQAHgA+GhEXAy8BK9LhTwHfj4hLgPcBd0ma3KTympnZMGTa3STpQ8ALwOY6p7wB6I6IZ9P+o8DPpe0AZkgScDqwD+jNsLgTzuPP7uGF1460uhhmNoY1K0jcLGmjpPslzQGQdDpwC/CbJ8n3GtAmqT3tfxg4N21/Efgh4FXgKeDTEVEa/AaSVkrqkNSxZ8+eJl3OxLDi/ie44gv/1OpimNkYNqwgIemxNGYw+LUcuAc4H1gK7ADuStluo9wNdbje+0b5iUcfBe6W9ATwOtCXDn8A2ACcld77i5Jm1niP1RHRHhHt8+bVvKvczMxO0bCW5YiIq4ZznqR7gbVp953AhyXdCcwGSpKORcQXB733vwOXpfxXAxekQ9cBd6RAskXSC8BbgCeGUxYbvmM9fUydVGx1McxsDGrG7KYFVbvXAJsAIuKyiFgUEYuAPwB+Z3CASPnPTF+nUO6eWpUOvQxcmY7NBy4Enm+0vHail/YebXURzGyMasaYxJ2SnpK0EbgC+JWhMkh6RNJZaffXJD0NbAS+ERH/mNJvB94j6SlgHXBLRLzWhPLaIB68NrN6Gl4FNiKuHcY5tw3aX1a1/WvAr9XI8ypwdaPls6G9uNdBwsxq8x3XE1ixIABedEvCzOpwkJjA+koBuLvJzOpzkJigSilAADy763XKk8jMzI7nIDFB9aWgMH/mFPYf7eHVg8daXCIzG4scJCaoSlfTj5wzG4BN2w+2sjhmNkY5SExQpdSSuOismRQEmx0kzKwGB4kJqtKSOG1yG28+83Q2OkiYWQ0OEhNUKS2VWCiI9kVz6XhxP719J6yfaGYTnIPEBFUZuC4K3nv+GRzu6uV7r7g1YWbHc5CYoCrdTcWCePf5bwDg37Z41RMzO56DxARVGbguFMTc0yZzyTmzePTpXS0ulZmNNQ4SE1R/S0LlpTl++pKz2PjKQZ7fU/fxH2Y2ATlITFCVIFEoDAQJCdb85yutLJaZjTEOEhNUKY5vScyfOZWf+KH5/OW/v8ThLj9K3MzKHCQmqOqB64pPXfFmDh3r5b5/8bOdzKzMQWKCqh64rrjk3Nn81I8s4E+/vZUf7DzUqqKZ2RjiIDFBVe6bq3Q3Vdz2Mxcxc9okrv+LDl490NmCkpnZWNKMZ1zfJmm7pA3ptSylL5LUWZW+qk7+uZIelfRc+jonpUvSH0naImmjpB9ttKw2YKC76fj0M06fwp9/4u0c6uxh+Z/8K/+weaeXETebwJrVkrg7Ipam1yNV6Vur0m+qk/dWYF1ELKH8LOtbU/pPAkvSayVwT5PKalR1Nw1qSQD88Dmz+Movvoc50yex8i/X81N//P9Y/c9b+d62AxzxoLbZhNLwM66bYDnwvrT9APBPwC0p/f9E+d/Y70iaLWlBROxodgFe2X+U1f88sQZr9x7uBo4fuK524Rtn8M1fvowHO7bxV995md955Af9x+bNmMLc6ZOZOa2NGVMnMblYoFgUbQVRLIhJhQKFgqgRfwCok3yS82sfqHe+jR1XXHgmV7zlzKa937GePl7ce4SdB49xsLOHA0d7ONrdR29fiZ5S0NtXorcU9PYFwUALeKjG8ODWch7bzhefNYuPvP3cpr9vs4LEzZI+DnQAvxoR+1P6YklPAoeA/x0R/1Ij7/yqD/6dwPy0fTawreq8V1LacUFC0krKLQ0WLlx4SoU/cLSHb3zv1VPKm2dnz57G4jNOq3t8UrHAx975Jj72zjex69Ax1r+0nxdeO8LLe49yoLObQ5297Dp0jN6+oLc08MfZVwp6S/X+zGqn1/sjrvsu7gIb814/1sum7QcbDhLP7Hydr6zfxuPP7mHL7sPU/dUC2gqirSjaCoUT/okY/D+FBp0w1PljXXdvqXVBQtJjwBtrHPoc5W6g2yn/Pd8O3AV8kvKH+cKI2CvpbcDXJV0UEXWnzURESBrRX39ErAZWA7S3t5/SJ8fFZ8/iyV+/+lSyThjzZ05l2Q8vaHUxLEc+8edPsP9I9ynn33u4i994eDNrN+5gcrHAO8+bywcvXsCbzzyds2dPZfb0ycyaNonTJrelwKATPvitccMKEhFx1XDOk3QvsDbl6QK60vZ6SVuBCyi3NqrtqnQjSVoA7E7p24HqsHhOSjOzHChK/asNj9SW3Yf5r/d9l71Huvjl97+Z6967mDmnTW5yCW04mjG7qfrfy2uATSl9nqRi2j6P8gB0rY7/h4EVaXsF8FBV+sfTLKd3AQezGI8ws2wUCuJUHlGy93AX1/7Zd+ktBV/7pffymasvdIBooWaMSdwpaSnl7qYXgRtT+uXA5yX1ACXgpojYByDpPmBVRHQAdwAPSroeeAn4SMr/CLAM2AIcBa5rQlnNbJQUJUonG0Co45Y1G9l3pJs1v/geLj57VgYls5FoOEhExLV10tcAa+ocu6Fqey9wZY1zAvhUo+Uzs9YoFkRvaWRNiW8/s5vHnt7N/1r2FgeIMcJ3XJtZJgoFnXQmUi1/tO45Fs6dzifeszibQtmIOUiYWSaKGrizfzi+t+0AT758gE++dxGT2/zRNFa4JswsE+WB6+EHib/t2Mb0yUV+7m3nZFgqGykHCTPLRFtB/cu/DKWvFPzD5l1cceGZzJg6KeOS2Ug4SJhZJoojaEk8+fJ+XjvcxQcurnXPrrWSg4SZZaKg4bckHn92D8WCeN+F8zIulY2Ug4SZZWIkLYn/eHEfF501k5nuahpzHCTMLBMFDS9IdPeW2LDtAO1vmjsKpbKRcpAws0wMtyWx6dWDHOsp8fZFc0ahVDZSDhJmloliYXgL/G3cdgCASxc6SIxFDhJmlomCxHBW5Xhm12FmT5/E/JlTsi+UjZiDhJllolhgWC2JZ3Ye4sL5M/wsiDHKQcLMMlEcxsB1RPDsrsO85Y0zRqlUNlIOEmaWiUJ6fvrJlgvffqCTw129XOAgMWY5SJhZJoqp++hkXU7P7T4MwAXzHSTGKgcJM8tEsZiCxElaEtv2HQXgTXOnj0qZbOQcJMwsE5WWxMmW5ti27yhT2grMm+GZTWNVQ0FC0m2StkvakF7LUvoiSZ1V6avq5J8r6VFJz6Wvc1L6xyRtlPSUpH+TdEkj5TSz0VdMYxK9J21JdHLOnGme2TSGNaMlcXdELE2vR6rSt1al31Qn763AuohYAqxL+wAvAD8eET8M3A6sbkI5zWwUFTT0wPW2/Uc5111NY1qru5uWAw+k7QeADwFExL9FxP6U/h3ATyExy5lKS2KoMYmFDhJjWjOCxM2pa+j+SndRsljSk5Iel3RZnbzzI2JH2t4JzK9xzvXA3zWhnGY2ilKMoF6IONjZw6FjvZw7x0FiLGsb6gRJjwG1ngTyOeAeyt1Bkb7eBXwS2AEsjIi9kt4GfF3SRRFxqN73iYiQdNzvk6QrKAeJHztJ+VYCKwEWLlw41OWY2WgZYuB6+/5OAM6ZM23UimQjN2SQiIirhvNGku4F1qY8XUBX2l4vaStwAdAxKNsuSQsiYoekBcDuqvf7EeA+4CcjYu9JyreaNGbR3t4+/Afqmlmm+lsSdf4qd71+DIAzZ04dpRLZqWh0dtOCqt1rgE0pfZ6kYto+D1gCPF/jLR4GVqTtFcBDKc9C4KvAtRHxbCNlNLPWqAxc1wsSe17vAuBMT38d04ZsSQzhTklLKXc3vQjcmNIvBz4vqQcoATdFxD4ASfcBqyKiA7gDeFDS9cBLwEdS/l8H3gD8aZoa1xsR7Q2W1cxGUaUlUa+7qRIkfI/E2NZQkIiIa+ukrwHW1Dl2Q9X2XuDKOufcMDjdzPJDQ4xJ7Hm9ixlT25g6qTiaxbIRavUUWDMbpyq3x9Xrbtr9+jF3NeWAg4SZZWKoMYndh7rc1ZQDDhJmlolC+nSp2910uIszZ3hm01jnIGFmmSicZEwiItySyAkHCTPLxMDA9YnHjnT30dnT5zGJHHCQMLNMDAxcnxgl9h4uT3+de9rkUSyRnQoHCTPLRP/AdY1jB472ADBnuoPEWOcgYWaZONnNdAc6y0Fi9vRJo1kkOwUOEmaWif4xidKJxw4c7QZgtlsSY56DhJll4qQtiaNuSeSFg4SZZUInuZmuEiRmTXOQGOscJMwsEwMPHao1JtHN6VPamFT0R9BY5xoys0wUTnKfxMGjPe5qygkHCTPLhE4yJrH/aLeDRE44SJhZJgYW+Ks9BXb2NM9sygMHCTPLxEBL4sRjB4/2MMstiVxwkDCzTJxsqfByS8JBIg8cJMwsE/XGJCKCg50euM6LhoOEpNskbZe0Ib2WpfRFkjqr0lfVyT9X0qOSnktf5ww6/nZJvZI+3GhZzWz01Fsq/Eh3H32lYOZUB4k8aFZL4u6IWJpej1Slb61Kv6lO3luBdRGxBFiX9gGQVAR+F/iHJpXTzEZJve6mw8d6ATh9attoF8lOwVjobloOPJC2HwA+VHXsvwFrgN2jXSgza0y97qbDXSlITHGQyINmBYmbJW2UdP+g7qLFkp6U9Liky+rknR8RO9L2TmA+gKSzgWuAe072jSWtlNQhqWPPnj2NXoeZNUn/HdeDWxIOErkyrCAh6TFJm2q8llP+ED8fWArsAO5K2XYACyPiUuAzwJclzTzZ94nyhOrKr9QfALdERI01JI/Lszoi2iOifd68ecO5HDMbBao3JuEgkSvDqqWIuGo450m6F1ib8nQBXWl7vaStwAVAx6BsuyQtiIgdkhYw0LXUDvxN+kU7A1gmqTcivj6csphZa9Ubk3g9jUmc5iCRC82Y3bSgavcaYFNKn5cGnpF0HrAEeL7GWzwMrEjbK4CHACJicUQsiohFwFeAX3KAMMuPekuFV1oSMzxwnQvNqKU7JS2l3E30InBjSr8c+LykHqAE3BQR+wAk3QesiogO4A7gQUnXAy8BH2lCmcysxUTtBf48JpEvDddSRFxbJ30N5ZlJtY7dULW9F7hyiO/xiQaKaGYtoP6B69qzm9zdlA9jYQqsmY1D9ZYKP9zVy6SimNLmj588cC2ZWSYK6dPlhJbEsV5Om9LWP/vJxjYHCTPLRL2WxJGuXo9H5IiDhJllotJOGDy76XUHiVxxkDCzTFS6kwavFH74mINEnjhImFkmCnVmNx3p7vXifjniIGFmmai3VHhl4NrywUHCzDLRHyQGrb52uKuXGQ4SueEgYWaZqLdU+JGuXqZPdpDICwcJM8tE/x3XVWkRQWdPH9MnF1tSJhs5Bwkzy8TAKrADYaK7r0QpYJqDRG44SJhZJmrdTHesuzxAMW2Sg0ReOEiYWSZqLRV+tKe8uJ9bEvnhIGFm2egPEgNJnd19gFsSeeIgYWaZKOjEh1x39qQg4ZZEbjhImFkmao1JuCWRPw4SZpaJWmMSbknkT0NBQtJtkrZL2pBey1L6IkmdVemr6uSfK+lRSc+lr3Oqjr0v5d0s6fFGymlmo09uSYwLzbjt8e6I+EKN9K0RsXSIvLcC6yLiDkm3pv1bJM0G/hT4YES8LOnMJpTTzEZRrceXuiWRP63ubloOPJC2HwA+lLb/C/DViHgZICJ2t6BsZtaAgZvpBtLcksifZgSJmyVtlHR/dXcRsFjSk5Iel3RZnbzzI2JH2t4JzE/bFwBzJP2TpPWSPt6EcprZKDrZmISX5ciPIbubJD0GvLHGoc8B9wC3U16e5XbgLuCTwA5gYUTslfQ24OuSLoqIQ/W+T0SEpMpvUxvwNuBKYBrw75K+ExHP1ijfSmAlwMKFC4e6HDMbJTVnN6UgMdUtidwYMkhExFXDeSNJ9wJrU54uoCttr5e0lXLroGNQtl2SFkTEDkkLgEq30ivA3og4AhyR9M/AJcAJQSIiVgOrAdrb2wc/BMvMWqTWKrCd3X1IMKWt1T3dNlyNzm5aULV7DbAppc+TVEzb5wFLgOdrvMXDwIq0vQJ4KG0/BPyYpDZJ04F3Ak83UlYzG13ixAX+Orv7mD6p2D/zyca+Rmc33SlpKeXupheBG1P65cDnJfUAJeCmiNgHIOk+YFVEdAB3AA9Kuh54CfgIQEQ8LenvgY0p/30RsanBsprZKCqceMM1nT19ntmUMw0FiYi4tk76GmBNnWM3VG3vpTzuUOu83wN+r5HymVnr1Lvj2uMR+eKOQTPLRM0xCT9wKHccJMwsE5KQTryZzvdI5IuDhJllRhzf3XS022MSeeMgYWaZKUhE1VOuj7klkTsOEmaWmYJ0wsC1WxL54iBhZpmRjh+4Ptbbx9Q2B4k8cZAws8yUB64H9rt7S0yZ5I+dPHFtmVlmChKlqv6mrt4Sk4v+2MkT15aZZaY8cD2gq6fEFA9c54qDhJllZvCYRFdvnxf3yxnXlpllpiD1j0n09pUohVeAzRvXlpllprol0dVbAmCyg0SuuLbMLDPVLYlKkJjiKbC54iBhZpkpHNeSKD+Vzt1N+eLaMrPMqOqO6+5KS8L3SeSKa8vMMlOoWgW2f0yi6O6mPHGQMLPMCA10N/VUxiT8sZMnri0zy0yhalmO/jEJdzflSsO1Jek2SdslbUivZSl9kaTOqvRVdfLPlfSopOfS1zkpfZakb0j6nqTNkq5rtKxmNrpqjUl4WY58aVZt3R0RS9Prkar0rVXpN9XJeyuwLiKWAOvSPsCngO9HxCXA+4C7JE1uUnnNbBQUCieOSXhZjnwZCyF9OfBA2n4A+FDaDmCGJAGnA/uA3tEvnpmdqvLzJDwFNs+aVVs3S9oo6f5Kd1GyWNKTkh6XdFmdvPMjYkfa3gnMT9tfBH4IeBV4Cvh0RJQGZ5a0UlKHpI49e/Y06XLMrBmqH186cDOdg0SeDKu2JD0maVON13LgHuB8YCmwA7grZdsBLIyIS4HPAF+WNPNk3yfK7dLKamAfADYAZ6X3/mKt/BGxOiLaI6J93rx5w7kcMxsl1avAelmOfGobzkkRcdVwzpN0L7A25ekCutL2eklbgQuAjkHZdklaEBE7JC0Adqf064A7UuDYIukF4C3AE8Mpi5m1Xq21m7wsR740Y3bTgqrda4BNKX2epGLaPg9YAjxf4y0eBlak7RXAQ2n7ZeDKlH8+cGGd/GY2RpXXbqrcJ+EpsHk0rJbEEO6UtJRyN9GLwI0p/XLg85J6gBJwU0TsA5B0H7AqIjqAO4AHJV0PvAR8JOW/HfgLSU9R7tq8JSJea0J5zWyUlJ9MV97u7vMU2DxqOEhExLV10tcAa+ocu6Fqey+pxTDonFeBqxstn5m1znHdTb7jOpdcW2aWmeqb6bp6S0xuK1Ce1W554SBhZpkpHncznR9dmkeuMTPLTEGiLwWJ7t6Sg0QOucbMLDOFQd1Nnv6aPw4SZpaZgqBUGrhPwi2J/HGNmVlmioXq50n0+W7rHHKNmVlmJNGXWhLdfW5J5JFrzMwyU5QGHjrU4zGJPHKQMLPMFAr0z27q6u3zkhw55Bozs8wc/zyJkpfkyCHXmJllprx2U9V9Em5J5I5rzMwyU57dVN72fRL55CBhZpkpiP7ZTV6WI59cY2aWmRPGJBwkcsc1ZmaZqQ4SXrspn1xjZpaZYkFV3U0ek8gjBwkzy4wEEdDbV6KvFO5uyqGGakzSbZK2S9qQXstS+iJJnVXpq+rk/3lJmyWVJLUPOvZZSVskPSPpA42U08xao1goLxVeeXSpu5vypxnPuL47Ir5QI31rRCwdIu8m4GeBL1UnSnor8FHgIuAs4DFJF0REXxPKa2ajpJjGJPzo0vxqaY1FxNMR8UyNQ8uBv4mIroh4AdgCvGN0S2dmjZJEqVQejwCYMsljEnnTjCBxs6SNku6XNKcqfbGkJyU9LumyEb7n2cC2qv1XUpqZ5UixQLkl0VvuBPCyHPkzZI1JekzSphqv5cA9wPnAUmAHcFfKtgNYGBGXAp8BvixpZhYXIGmlpA5JHXv27MniW5jZKSqkpcK7+1sSDhJ5M+SYRERcNZw3knQvsDbl6QK60vZ6SVuBC4COYZZrO3Bu1f45Ka1W+VYDqwHa29tjmO9vZqOgkJbl6O9u8hTY3Gl0dtOCqt1rKA9EI2mepGLaPg9YAjw/grd+GPiopCmSFqf8TzRSVjMbfQUN6m7ywHXuNDq76U5JS4EAXgRuTOmXA5+X1AOUgJsiYh+ApPuAVRHRIeka4I+BecA3JW2IiA9ExGZJDwLfB3qBT3lmk1n+eHZT/jUUJCLi2jrpa4A1dY7dULX9NeBrdc77beC3GymfmbVW5fGlXb5PIrdcY2aWmWKh/PjSgZaExyTyxkHCzDJTWSrcYxL55Rozs8yUZzdF1ewmf+TkjWvMzDJTWSrc90nkl2vMzDJTrAxc+z6J3HKQMLPMVG6mO9ZTHpNwd1P+uMbMLDNtBQEDQcJrN+WPa8zMMlNMQeJodx+TiwUKad/yw0HCzDLT1h8kej39Nadca2aWmUpL4khXn8cjcsq1ZmaZmZTGII529zpI5JRrzcwyc1xLwk+lyyUHCTPLzHFjEp7ZlEuuNTPLTKUlcbirl6m+2zqXXGtmlpm2orub8s5BwswyUyyUP2LKLQkHiTxykDCzzLRVdzd5dlMuudbMLDNtVXdYuyWRTw0FCUm3SdouaUN6LUvpiyR1VqWvqpP/5yVtllSS1F6V/hOS1kt6Kn19fyPlNLPWqIxJAB64zqmGnnGd3B0RX6iRvjUilg6RdxPws8CXBqW/Bvx0RLwq6WLgW8DZjRfVzEZTZUwC3JLIq2YEiVMWEU9D+WHpg9KfrNrdDEyTNCUiukaxeGbWIHc35V8z2n83S9oo6X5Jc6rSF0t6UtLjki5r4P1/DvjPegFC0kpJHZI69uzZ08C3MbNmKzpI5N6QQULSY5I21XgtB+4BzgeWAjuAu1K2HcDCiLgU+AzwZUkzR1o4SRcBvwvcWO+ciFgdEe0R0T5v3ryRfgszy9DxLQmPSeTRkN1NEXHVcN5I0r3A2pSnC+hK2+slbQUuADqGWzBJ5wBfAz4eEVuHm8/Mxo7jWhJ+dGkuNTq7aUHV7jWUB6KRNE9SMW2fBywBnh/B+84GvgncGhH/2kgZzax1JhU9cJ13jbb/7kzTVDcCVwC/ktIvBzZK2gB8BbgpIvYBSLqvMt1V0jWSXgHeDXxT0rdS/puBNwO/XjWN9swGy2pmo6y6i8ndTfnU0OymiLi2TvoaYE2dYzdUbX+NcpfS4HN+C/itRspmZq03paqLyS2JfHJoN7PMVAcGtyTyybVmZpk5rrvJA9e55CBhZpmpbkl4qfB8cpAws8xUz26aNW1SC0tip8pBwsxGxezpDhJ55CBhZqNitlsSueQgYWajoq3oj5s8aukqsGY2/j1447v5wc5DrS6GnSIHCTPL1DsWz+Udi+e2uhh2itz+MzOzuhwkzMysLgcJMzOry0HCzMzqcpAwM7O6HCTMzKwuBwkzM6vLQcLMzOpSRLS6DE0jaQ/w0ilmPwN4rYnFyQNf88Tga54YGrnmN0XEvFoHxlWQaISkjohob3U5RpOveWLwNU8MWV2zu5vMzKwuBwkzM6vLQWLA6lYXoAV8zRODr3liyOSaPSZhZmZ1uSVhZmZ1OUiYmVldDhKApA9KekbSFkm3tro8zSLpXEnflvR9SZslfTqlz5X0qKTn0tc5KV2S/ij9HDZK+tHWXsGpkVSU9KSktWl/saTvpuv6W0mTU/qUtL8lHV/UynI3QtJsSV+R9ANJT0t693iuZ0m/kn6nN0n6a0lTx2M9S7pf0m5Jm6rSRlyvklak85+TtGIkZZjwQUJSEfgT4CeBtwK/IOmtrS1V0/QCvxoRbwXeBXwqXdutwLqIWAKsS/tQ/hksSa+VwD2jX+Sm+DTwdNX+7wJ3R8Sbgf3A9Sn9emB/Sr87nZdXfwj8fUS8BbiE8vWPy3qWdDbwy0B7RFwMFIGPMj7r+S+ADw5KG1G9SpoL/AbwTuAdwG9UAsuwRMSEfgHvBr5Vtf9Z4LOtLldG1/oQ8BPAM8CClLYAeCZtfwn4harz+8/Lyws4J/3hvB9YC4jyXahtg+sb+Bbw7rTdls5Tq6/hFK55FvDC4LKP13oGzga2AXNTva0FPjBe6xlYBGw61XoFfgH4UlX6cecN9ZrwLQkGfuEqXklp40pqYl8KfBeYHxE70qGdwPy0PR5+Fn8A/E+glPbfAByIiN60X31N/debjh9M5+fNYmAP8Oepm+0+SacxTus5IrYDXwBeBnZQrrf1jP96rhhpvTZU3w4SE4Ck04E1wH+PiEPVx6L8r8W4mAct6aeA3RGxvtVlGWVtwI8C90TEpcARBroggHFXz3OA5ZSD41nAaZzYJTMhjEa9OkjAduDcqv1zUtq4IGkS5QDxVxHx1ZS8S9KCdHwBsDul5/1n8V7gZyS9CPwN5S6nPwRmS2pL51RfU//1puOzgL2jWeAmeQV4JSK+m/a/QjlojNd6vgp4ISL2REQP8FXKdT/e67lipPXaUH07SMB/AEvSzIjJlAfAHm5xmZpCkoA/A56OiN+vOvQwUJnhsILyWEUl/eNplsS7gINVzdoxLyI+GxHnRMQiyvX4jxHxMeDbwIfTaYOvt/Jz+HA6P3f/bUfETmCbpNhXuToAAADkSURBVAtT0pXA9xmn9Uy5m+ldkqan3/HK9Y7req4y0nr9FnC1pDmpFXZ1ShueVg/KjIUXsAx4FtgKfK7V5Wnidf0Y5aboRmBDei2j3B+7DngOeAyYm84X5ZleW4GnKM8eafl1nOK1vw9Ym7bPA54AtgD/F5iS0qem/S3p+HmtLncD17sU6Eh1/XVgzniuZ+A3gR8Am4C/BKaMx3oG/pryuEsP5Rbj9adSr8An0/VvAa4bSRm8LIeZmdXl7iYzM6vLQcLMzOpykDAzs7ocJMzMrC4HCTMzq8tBwszM6nKQMDOzuv4/hmMsDYz8nckAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import nest\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "params =  {\n",
    "    \"C_m\": 250.0,\n",
    "    \"E_L\": -70.0,\n",
    "    \"t_ref\": 2.0,\n",
    "    \"tau_m\": 10.0,\n",
    "    \"tau_minus\": 20.0,\n",
    "    \"tau_syn_ex\": 0.2,\n",
    "    \"tau_syn_in\": 2.0,\n",
    "    \"V_reset\": -70.0,\n",
    "    \"V_th\": -55.0\n",
    "}\n",
    "\n",
    "g_L = params['C_m'] / params['tau_m']\n",
    "I_th = (params['V_th'] - params['E_L']) * g_L\n",
    "Vth = params['V_th']\n",
    "\n",
    "nest.ResetKernel()\n",
    "\n",
    "n = nest.Create('iaf_psc_alpha', params=params)\n",
    "\n",
    "sg = nest.Create('spike_generator')\n",
    "vm = nest.Create('voltmeter')\n",
    "\n",
    "nest.Connect(sg, n, syn_spec={'weight': 10.0})\n",
    "nest.Connect(vm, n)\n",
    "\n",
    "nest.SetStatus(n, {'V_th': 0., 'I_e': I_th, 'V_m': Vth})\n",
    "nest.SetStatus(sg, {'spike_times': [200., 700.], 'spike_weights': [1., -1.]})\n",
    "nest.Simulate(1000.)\n",
    "\n",
    "events = nest.GetStatus(vm, 'events')[0]\n",
    "times, V_m = events['times'], events['V_m']\n",
    "\n",
    "print(np.max(V_m) - Vth, np.min(V_m) - Vth)\n",
    "plt.plot(times,V_m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}