{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#AUTHOR: Lisa Blum Moyse\n",
    "#        lisa.blum-moyse@inria.fr\n",
    "#\n",
    "# REFERENCE: Blum Moyse & Berry. Modelling the modulation of cortical Up-Down state switching by astrocytes\n",
    "#\n",
    "# LICENSE: CC0 1.0 Universal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import math as math\n",
    "import cmath\n",
    "import scipy.integrate as integrate\n",
    "from scipy.misc import derivative\n",
    "import scipy.special as sc\n",
    "import scipy as sci\n",
    "import mpmath as mp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def M(a,b,z):\n",
    "    return mp.hyp1f1(a, b, z)\n",
    "\n",
    "def U(y,latau):\n",
    "    return np.exp(y**2)/sc.gamma((1+latau)/2)*M((1-latau)/2,0.5,-y**2) + 2*y*np.exp(y**2)/sc.gamma(latau/2)*M(1-latau/2,1.5,-y**2)\n",
    "\n",
    "def partial_derivative(func, var=0, point=[]):\n",
    "    args = point[:]\n",
    "    def wraps(x):\n",
    "        args[var] = x\n",
    "        return func(*args)\n",
    "    return derivative(wraps, point[var], dx = 1e-6)\n",
    "\n",
    "def der(U,ytr,latau):\n",
    "    return partial_derivative(U, 0, [ytr,latau]) \n",
    "\n",
    "def R(latau,yr,yt,r0,sig,tau):\n",
    "    return r0/(sig*(1+latau))*(der(U,yt,latau)-der(U,yr,latau))/(U(yt,latau)-U(yr,latau))\n",
    "\n",
    "def S(la,D,taur,taud):\n",
    "    return np.exp(-la*D)/((1+la*taur)*(1+la*taud))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X= [0.05965212-6.6258458e-19j] lambda= [0.01+1.e-20j]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAEJCAYAAACwph1QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbHElEQVR4nO3df5BdZZ3n8fenu9OJ4WdI+BECGXEGQdwxK0ZAYRRGQZIZJ+rqLIgaXagUu0Jp7bgrFrXO1DpbheM4pdYKqV5kxBlLZlaykrWiLDCiM4tBAouBEIHwY0IngZAQfoaQ7tvf/eOcm7l07o/TOef2Pefm86o61fec89znPH3uzTdPP+f5oYjAzMzKb6DXBTAzs2wcsM3MKsIB28ysIhywzcwqwgHbzKwiHLDNzCoiV8CWdJSk2yQ9mv6c0yLdkZJ+KOk3kjZKelee65qZlZmkGyRtl/Rgi/OS9C1JmyStl3R6lnzz1rCvAu6IiJOBO9L9Zr4J/DQiTgUWARtzXtfMrMy+C1zY5vwS4OR0WwFclyXTvAF7GXBj+vpG4EOTE0g6HHgP8B2AiNgbEc/nvK6ZWWlFxC+A59okWQZ8LxJrgSMlze+U71DOch0bEdvSAm6TdEyTNG8CngX+WtIi4F7gcxHxSrMMJa0g+R+HQ2brHaf+znDOIppZv3vyqTF2PFdTnjw+cN4hsfO5Wqa0965/bQOwp+HQSESMTOFyC4CnGvZH02Pb2r2pY8CWdDtwXJNTV2cs2BBwOnBlRNwt6ZskTSf/pVni9JceAVi8aFb86tYTM17GzA5WZ3zgqc6JOtj5XI1f3bowU9rB+Y/uiYjFOS7X7D+XjvOEdAzYEfH+lleUnpE0P61dzwe2N0k2CoxGxN3p/g9p3dZtZtYTAUwwMV2XGwUaa6MnAFs7vSlvG/ZqYHn6ejlwy+QEEfE08JSkU9JD7wMeynldM7NCBcFY1DJtBVgNfCrtLXIW8EK9ebmdvG3Y1wB/L+lSYDPwMQBJxwPXR8TSNN2VwPclDQOPA5/JeV0zs8IVVcOW9APgXGCepFHgT4EZABGxElgDLAU2AbvJGBNzBeyI2ElSY558fGtamPr+/UCe9h4zs64KglpB001HxMUdzgfw2anmm7eGbWbWNyY6P/frqVIH7M1jh3DFljN7XQwzK7nNY+26PGcTQM0B+8C9PDbML5/+rV4Xw8xK7uWxYsZruIadQ+3VIV58aG6vi2FmJVd7NX8oC2Cs5EsmljpgD9Rg+IVcg5fM7CAwUEBPuyDcJJLHwDjMfrrcN9DMem9gvIBMAmolDzflDthjwezthXRSN7M+NjCWP9ImIx3LrdwBe/deDr13c6+LYWYlN7B7bwG5iFrTKT7Ko9QBm6EhJuYe2etS9BWV/KGKtRYqdzDpqV1FPXQs9z0udcCOoQHG5s3udTHMrOTiyfyrHSb9sB2wD9j4TPH878zsdTHMrOTGHygm0E64hn3gYghePabcN9DMei8KiGSuYed1SI2hM3b1uhRmVnY35e9NFoha7hmnu6vUAfuI4T1csPA3vS6GmZXcs8N7OifKwE0iOZwwYzdfO+7/9boYZlZyP5+xO3cegdgbgwWUpntKHbDNzKZLMnDGTSJmZpXgh45mZhUQIWrhGraZWSVMuIZtZlZ+yUPHcofEcpfOzGya+KGjmVmF1NwP28ys/DzS0cysQibcS8TMrPySyZ8csA/Y1vFZ/Omzb+11Mcys5LaOP5s7j0CMeWj6gdtdG+aBF47vdTHMrOR214Zz5xGBB87ksWdsBg9uccA2s/b2jM0oIBf198AZSUcBfwe8EXgS+OOIaDqBtaRBYB2wJSL+MEv+MS7GnveKM2bWXoznD7RB/9ewrwLuiIhrJF2V7n+xRdrPARuBw7NmrpqYsavcbUpm1nuqFVMz7veHjsuAc9PXNwJ30iRgSzoB+APgvwH/MWvmqsHwi+X+E8XMek/5F5whUN8vYHBsRGwDiIhtko5pke4bwH8GDuuUoaQVwAqAGYfNYWBvzhKaWf+LYrIYq/pcIpJuB45rcurqLBeQ9IfA9oi4V9K5ndJHxAgwAjDzxBNj9/wCPgkz62sTRTxzRNWfDzsi3t/qnKRnJM1Pa9fzge1Nkp0N/JGkpcAs4HBJfxsRn+h0bQUMFPCnjpn1NxVUw+73kY6rgeXANenPWyYniIgvAV8CSGvYX8gSrAFmvBwc/4vxnEU0s373zMvF/CVe+Rp2B9cAfy/pUmAz8DEASccD10fE0ly5B2giZwnNrP8VUcMOFVrDlnQh8E1gkCQeXjPp/BHA3wILSWLxX0bEX7fLM1fAjoidwPuaHN8K7BesI+JOkp4k2fIfFK8d4W59ZtZeDBbTD7uooenpuJNvA+cDo8A9klZHxEMNyT4LPBQRH5R0NPCwpO9HRMuuFqV+JFqbCS+eVO42pQPi56hWNeVuKaBWyPi6Qtd0PAPYFBGPA0i6iaQbdGPADuAwSQIOBZ4D2rYBlzpgT8wKdp/6Wq+LYWYlNzErfy0oeeiY+X+meZLWNeyPpD3c6hYATzXsjwJnTsrjv5M8B9xK0uX530ZE20bgUgfs3z18B7+64Du9LoaZldwZX99RSD5TGOm4IyIWtznfLPJP/l/lA8D9wO8Dvw3cJukfI+LFVpn2YXuDmdnU1Uc6ZtkyGAVObNg/gaQm3egzwKpIbAKeAE5tl6kDtplZaoKBTFsG9wAnSzpJ0jBwEUnzR6PNpJ02JB0LnAI83i7TUjeJmJlNlwgYmyimDhsR45KuAG4l6dZ3Q0RskHR5en4l8BXgu5IeIGlC+WJEtG3bccA2M6PeJFJco0NErAHWTDq2suH1VuCCqeTpgG1mlur3kY5mZn1hit36esIB28wMoOAmkW4odcDe8PTRLPqL/9DrYphZyW16+q8Kyaev13TsthiEvUf0uhRmVnZFTAGS9BIp99xF5Q7YQ7D3SE/XZ2btFbFQzMGwRFhXqQYzXir3DTSz3itiTUdwk0gumoDBPeW+gWbWe0XMm+9eIjlpAma83OtSmFnZFbXQiXuJ5BEwsNeTR5tZBwWtODPugH3gNAFDe3pdCjMru+Jq2G4SMTMrPbdh5xSCiVKX0MzKoKg464CdQ20WvPDmXpfCzMquNit/Hu6HndPw7DFOePvkRRrMzF5v++yxQvJxP+wcTpn1PD976y29LoaZldwZs57PnUcEjBe0gEG3lDpgm5lNJzeJmJlVgNuwzcwqJBywzcyqoewPHXO1sEs6StJtkh5Nf85pkuZEST+TtFHSBkmfy3NNM7NuiEjasLNsvZL3kehVwB0RcTJwR7o/2TjwJxHxFuAs4LOSTst5XTOzgonaxECmrVfyXnkZcGP6+kbgQ5MTRMS2iLgvff0SsBFYkPO6ZmaFi1CmrVfytmEfGxHbIAnMko5pl1jSG4G3A3e3SbMCWAGwcIGb2M1sevTFXCKSbgeOa3Lq6qlcSNKhwM3A5yPixVbpImIEGAFYvGiW51Y1s+kRSTt2mXUM2BHx/lbnJD0jaX5au54PbG+RbgZJsP5+RKw64NKamXVRX/cSAVYDy9PXy4H9xpFLEvAdYGNEFLMWvZlZweIgeOh4DXC+pEeB89N9JB0vaU2a5mzgk8DvS7o/3ZbmvK6ZWeEism29kuupXkTsBN7X5PhWYGn6+p+g5H9nmJnhkY5mZpWQ1J4dsM3MKqHy3frMzA4Wle/WZ2Z2MAjEhBcwMDOrhpJXsHN36zMz6w9R7Fwiki6U9LCkTZKaTYyHpHPTrs4bJP28U56uYZuZ1RVUxZY0CHybZHzKKHCPpNUR8VBDmiOBa4ELI2Jzp7mYwDVsM7N9CqxhnwFsiojHI2IvcBPJ7KaNPg6siojNybWj6dQejRywzcxIZ+ubUKYNmCdpXcO2YlJ2C4CnGvZH2X9a6TcDcyTdKeleSZ/qVEY3iZiZQRKxs/fD3hERi9ucb5bR5AaXIeAdJKPF3wD8UtLaiHikVaYO2GZmqQL7YY8CJzbsnwBsbZJmR0S8Arwi6RfAIqBlwHaTiJlZXWTcOrsHOFnSSZKGgYtIZjdtdAvwe5KGJM0GziRZkasl17DNzAAobvmviBiXdAVwKzAI3BARGyRdnp5fGREbJf0UWA9MANdHxIPt8nXANjOrK3DkTESsAdZMOrZy0v7XgK9lzdMB28wMkoEzE578ycysIhywzcyqoeSTiThgm5nVOWCbmVXA1AbO9IQDtplZygsYmJlVhXuJmJlVg1zDNjOrgOzDznvGAdvMDAD5oaOZWWW4hm1mVhETvS5Aew7YZmbgfthmZlVS9l4ihSxg0Gk5dyW+lZ5fL+n0Iq5rZlao4hYw6IrcAbthOfclwGnAxZJOm5RsCXByuq0Arst7XTOzg00RNewsy7kvA74XibXAkZLmF3Bt6yerXkTvfBIdvwm980lY9WKvS2QHGUW2rVeKaMNutpz7mRnSLAC2Tc4sXS5+BcDs4w7hU//8ngKKaGX3rp8+xmXXPMHMPbXkwOg4e/9kJyPb38rPL3gze2tD7J0Y5LXaEGO1QfbWBhmrDTBeG2QikqWdIkACKRhQMDgwwdDgBMODNWYM1pg5OM7wYI3hgXGGB2oMKBhQybsFWCZP7L0lfybBQTE0Pcty7lnSJAcjRoARgCNmHhc7Lz4yV+GsGj66+T5m1mqvOzZzT40//vN1rPne24HkSzQrglmMA+NTv4gGqTHIqwzzav4iW4mMbymo/0TJHzoW8VtmXc69U5r9TMwa4tU3H5O7gFZ+Rz/xUvPjtZf8HbCOJnYWE7DL3kukiN9y33LuwBaS5dw/PinNauAKSTeRNJe8EBH7NYdMNnao2PJ7M3IXsOwfgsG2tXNY8MKu/Y8fMYet5+T/Dlh/G3ugoKaMkseK3AE7y3LuJCsHLwU2AbuBz2TKfNYEOuXlvEW0Cvj68vP485X/m9l7x/Yd2z08g68vPw9O9XfAOphV0LOIfg/Y0Hk594gI4LNTzffQ4dc4e+ET+QtopbfrksP5m7ln8W+uu4+5z7zCzmMP4eZ/fzq7Ljycs/F3wNp7bvi13Hn0ugdIFqUe6Thn6BU+Mm9dr4th0+UT8MtPnLRv93h28RH8+Vtn/zT0SjEZHQS9RLrmDRrnd4d39LoYZlZyb9AB9BpqwjXsHIY1yMKhQ3tdDDMruWHt/8D6gDhgm5lVgNuwzcwqxAHbzKwayj5TQSHTq5qZWfe5hm1mVucmETOzCvBDRzOzCnHANjOrCAdsM7PyE+4lYmZWDRmXB8vazt1pcfKGdO+UVJP00U55OmCbmdUVtGp6xsXJ6+m+SjI9dUcO2GZmdQUFbLItTg5wJXAzsD1Lpg7YZmapKTSJzJO0rmFbMSmrVguP/8u1pAXAh4GVZOSHjmZmddl7ieyIiMVtzmdZePwbwBcjoiZlm4fbAdvMDJKHjsX1Esmy8Phi4KY0WM8Dlkoaj4gftcrUAdvMrK64ftgdFyePiH3LK0n6LvDjdsEaHLDNzPYpamh6xsXJp8wB28ysrsCRjp0WJ590/NNZ8nTANjODqXTZ6xkHbDMz0qHpDthmZtXggG1mVhUO2GZmFeGAbWZWARVYcaaQuUQ6TSMo6RJJ69PtLkmLiriumVmhipv8qSty17AbphE8n2Q45j2SVkfEQw3JngDeGxG7JC0BRoAz817bzKxIZV/AoIgmkX3TCAJIqk8juC9gR8RdDenXkoyrNzMrlYOhSaTjNIKTXAr8pIDrmpkVJ2tzSJWbRMg2jWCSUDqPJGCf0zKzZF7ZFQALF/iZqJlNo4Oghp1lGkEkvQ24HlgWETtbZRYRIxGxOCIWHz13sIDimZl1Vh/pWNSajt1QRBW24zSCkhYCq4BPRsQjBVzTzKxwmih3FTt3wM44jeCXgbnAtelk3eMdVmswM5teB8vkT52mEYyIy4DLiriWmVm3lL2XiJ/qmZnVOWCbmVWDa9hmZlXhgG1mVgHFrpreFQ7YZmZ4xRkzs2qJckdsB2wzs5Rr2GZmVXCwDJwxM+sHfuhoZlYRDthmZlUQ+KGjmVlV+KGjmVlVOGCbmZWfB86YmVVFRP8vYGBm1jfKHa8dsM3M6twkYmZWBQG4ScTMrCLKHa8Z6HUBzMzKQpFty5SXdKGkhyVtknRVk/OXSFqfbndJWtQpT9ewzcxSRfUSkTQIfBs4HxgF7pG0OiIeakj2BPDeiNglaQkwApzZLl/XsM3M4F9m68uydXYGsCkiHo+IvcBNwLLXXS7irojYle6uBU7olKlr2GZm1AfOZK5hz5O0rmF/JCJGGvYXAE817I/SvvZ8KfCTThd1wDYzq8s+W9+OiFjc5ryaHGv6v4Gk80gC9jmdLuqAbWaWmkINu5NR4MSG/ROArftdT3obcD2wJCJ2dsrUbdhmZlB0G/Y9wMmSTpI0DFwErG5MIGkhsAr4ZEQ8kiVT17DNzAAobi6RiBiXdAVwKzAI3BARGyRdnp5fCXwZmAtcKwlgvEMziwO2mdk+BS5gEBFrgDWTjq1seH0ZcNlU8iykSaRTB/GGdO+UVJP00SKua2ZWmEiWCMuy9UrugN3QQXwJcBpwsaTTWqT7KsmfCGZm5RORbeuRImrYHTuIp64Ebga2F3BNM7PiFffQsSuKCNjNOogvaEwgaQHwYWAlHUhaIWmdpHXP7qwVUDwzs2w0MZFp65UiAnaWDuLfAL4YER0jcESMRMTiiFh89NzBAopnZpZBkAycybL1SBG9RLJ0EF8M3JR2XZkHLJU0HhE/KuD6Zma5iShy4ExXFBGw93UQB7aQdBD/eGOCiDip/lrSd4EfO1ibWen0e8DO2EHczKz8+j1gQ+cO4pOOf7qIa5qZFarehl1iHuloZpbqZQ+QLBywzcwA6O2gmCwcsM3MIB0U44BtZlYN5W4RccA2M6s7GPphm5n1BwdsM7MKiIBaudtEHLDNzOpcwzYzqwgHbDOzCgigoDUdu8UB28wMSAbOuA3bzKz8Aj90NDOrDLdhm5lVhAO2mVkVePInM7NqCMDTq5qZVYRr2GZmVeCh6WZm1RAQ7odtZlYRHuloZlYRbsM2M6uACPcSMTOrDNewzcyqIIhardeFaMsB28wMPL2qmVmllLxb30CvC2BmVgYBxERk2rKQdKGkhyVtknRVk/OS9K30/HpJp3fK0wHbzAySB44xkW3rQNIg8G1gCXAacLGk0yYlWwKcnG4rgOs65euAbWaWilot05bBGcCmiHg8IvYCNwHLJqVZBnwvEmuBIyXNb5dpqduw713/2o7B+Zv+ueBs5wE7Cs7zQLkszZWpLFCu8rgszZ2SN4OX2HXr7fHDeRmTz5K0rmF/JCJGGvYXAE817I8CZ07Ko1maBcC2VhctdcCOiKOLzlPSuohYXHS+B8Jlaa5MZYFylcdlaW5S8DwgEXFhEWVJqdklDiDN67hJxMyseKPAiQ37JwBbDyDN6zhgm5kV7x7gZEknSRoGLgJWT0qzGvhU2lvkLOCFiGjZHAIlbxLpkpHOSaaNy9JcmcoC5SqPy9JcmcpCRIxLugK4FRgEboiIDZIuT8+vBNYAS4FNwG7gM53yVZR87LyZmSXcJGJmVhEO2GZmFdE3ATvDMND/JOn+dHtQUk3SUem5JyU9kJ7L3T0oQ1nOlfRCQ3m+nPW9XSrPJenQ2PWS7pK0qOHcdN+blsN1u3Fv0nyPknSbpEfTn3OapDlR0s8kbZS0QdLnGs79maQtDZ/n0m6XJ03X9LPJ+v6iyiLplIbf/X5JL0r6fHqusHsj6WPpvZ+Q1LI7YavvSZH3pWciovIbSaP+Y8CbgGHg18BpbdJ/EPiHhv0ngXnTVRbgXODHeX+PAsvzbmBO+noJcHcP781S4CckfVTPqpelG/em4Zp/AVyVvr4K+GqTNPOB09PXhwGP1K8P/BnwhQK/zx3L0+6zyfr+Issy6TN+Gvitou8N8BaSATJ3Aoun+h0r8r70auuXGnaWYaCNLgZ+UJKyFPXeA84zIu6KiF3p7lqS/qDdkGe4bjfuTeM1b0xf3wh8aHKCiNgWEfelr18CNpKMSuuGjuXp8vvz5PU+4LGIKHqEMhGxMSIe7pCs3fekyPvSE/0SsFsN8dyPpNnAhcDNDYcD+D+S7pW0YprK8i5Jv5b0E0lvneJ7u1GeuktJarh1031vWqXpxr2pOzbS/q/pz2PaJZb0RuDtwN0Nh69Im3BuKOBP7azlafXZTOn3KagsdRexf2WoyHvTSbvvSZH3pSf6pR/2VIZ4fhD4vxHxXMOxsyNiq6RjgNsk/SYiftHFstxH8ifjy2mb3o9IZuya8lDVgsqTJJTOIwnY5zQcnu570ypNrnsj6XbguCanrs6aR5rPoST/2X8+Il5MD18HfCUtz1eArwP/bhrKU8hnU+C9GQb+CPhSw+Ep3Zt2ZYmIW7IUo8mxvum73C8BeypDPPerAUTE1vTndkn/i+TPqgMNSh3L0vAPnYhYI+laSfOm+HsUVh4ASW8DrgeWRMTOhvJN671pk2Y4y+/RSkS8v9U5Sc9Imh8R29Lml+0t0s0gCdbfj4hVDXk/05DmfwA/no7ytPlsMr2/yLKklgD3Nd6Pqd6bdmXJqN13bEr3pYz6pUkkyzBQJB0BvBe4peHYIZIOq78GLgAe7GZZJB0nSenrM0g+h51Zf48ulGchsAr4ZEQ80nB82u8NrYfrduPeNF5zefp6OQ3fj7r08/oOsDEi/mrSucYpMT9MvnuUtTztPpuO7y+yLA32ezbUhXvTSbvvSZH3pTd6/dSzqI2kd8EjJE+Ir06PXQ5c3pDm08BNk973JpInyb8GNtTf282yAFek1/o1yUO+d7d77zSU53pgF3B/uq3r4b0RycTvjwEP0NAboBv3Js13LnAH8Gj686j0+PHAmvT1OSR/Wq9vuE9L03N/k5Z1PUlQmD8N5Wn52bR6f7fKku7PJql0HDHp/YXdG5KAPwq8BjwD3NqiLE2/J0Xel15tHppuZlYR/dIkYmbW9xywzcwqwgHbzKwiHLDNzCrCAdvMuiYd3bhdUiHd+ZRM2lafSKqobp2V4V4iZtY1kt4DvEwyP8y/KiC/lyPi0PwlqybXsM2sayIZKt84DQSSflvST9M5UP5R0qk9Kl7lOGCb2XQbAa6MiHcAXwCuncJ7Z0laJ2mtpA91pXQl1i9ziZhZBaSTZ70b+J/p7AwAM9NzHwH+a5O3bYmID6SvF0Yy4dWbgH+Q9EBEPNbtcpeFA7aZTacB4PmI+NeTT0Qyodaq/d7x+jT1Ca8el3QnyTS3B03AdpOImU2bSGaqfELSx2DfknCLOryNNO0cSfXa+DzgbOChrhW2hBywzaxrJP0A+CVwiqRRSZcClwCXSqpPXJV15aC3AOvS9/0MuCYiDqqA7W59ZmYV4Rq2mVlFOGCbmVWEA7aZWUU4YJuZVYQDtplZRThgm5lVhAO2mVlF/H/+Q87/lQjSMwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# without astrocytes\n",
    "Ne = 4000; Ni = 1000\n",
    "Vth = 20 ; taub = 1; taue=20; taui=10; Jee = 280/Ne*taue; Jei=70/Ni*taue; Jie=500/Ne*taui; Jii=100/Ni*taui;  Vli=6.5; Vle=7.6; Vr=14;\n",
    "delaye = 0.5; delayi = 0.25; taude = 23; taure = 8  ;taudi = 1; tauri = 1  \n",
    "\n",
    "beta = 1\n",
    "Ka = 600\n",
    "\n",
    "sig = 4.427 ; sige = sig; sigi = sig; \n",
    "\n",
    "#ri0 = 0.00029499999999978415; re0 = 0.0002879999999997962 #stable down\n",
    "ri0 = 0.0012129999999998197; re0 =  0.0008939999999997962 #unstable\n",
    "#ri0 = 0.0017624999999998384; re0 =  0.0011219999999997961 #stable up\n",
    "\n",
    "Lre0=[re0]; Lri0=[ri0]\n",
    "Lla = []\n",
    "for re0 in Lre0:\n",
    "    for ri0 in Lri0:\n",
    "    \n",
    "        Ie0 =  Vle + taub*(Ne*Jee*re0-Ni*Jei*ri0)-Ka*beta*re0\n",
    "        yte = (Vth-Ie0)/sige; yre = (Vr-Ie0)/sige\n",
    "        Ii0 = Vli + taub*(Ne*Jie*re0-Ni*Jii*ri0)\n",
    "        yti = (Vth-Ii0)/sigi; yri = (Vr-Ii0)/sigi\n",
    "        \n",
    "        LRe = np.concatenate([np.arange(-0.6,-1e-13,1e-2),np.arange(1e-13,0.6,1e-2)])\n",
    "        LIm = np.concatenate([np.arange(-1e-5,-1e-20,1e-6),np.arange(1e-20,1e-5,1e-6)])\n",
    "        \n",
    "        X = np.zeros((len(LRe),len(LIm)),dtype='complex')\n",
    "        Mla = np.zeros((len(LRe),len(LIm)),dtype='complex')\n",
    "        for i in range(len(LRe)):\n",
    "            for j in range(len(LIm)):\n",
    "                la = complex(LRe[i],LIm[j])\n",
    "                Mla[i,j] = la\n",
    "                X[i,j] = taub*Ne*Jee*R(la*taue,yre,yte,re0,sige,taue)*S(la,delaye,taure,taude)*(1+taub*Ni*Jii*R(la*taui,yri,yti,ri0,sigi,taui)*S(la,delayi,tauri,taudi)) - taub*Ni*Jii*R(la*taui,yri,yti,ri0,sigi,taui)*S(la,delayi,tauri,taudi) - taub*Ni*Jei*R(la*taue,yre,yte,re0,sige,taue)*S(la,delayi,tauri,taudi)*Ne*Jie*R(la*taui,yri,yti,ri0,sigi,taui)*S(la,delaye,taure,taude) -1 \n",
    "        eps = 1\n",
    "        id = np.where(np.abs(X.real)+np.abs(X.imag)==np.min(np.abs(X.real)+np.abs(X.imag)))#np.where(np.abs(X.real)+np.abs(X.imag)<eps)\n",
    "        print('X=',X[id],'lambda=',Mla[id])\n",
    "        Lla.append(Mla[id].real)\n",
    "        plt.figure()\n",
    "        plt.imshow(np.abs(X.real)+np.abs(X.imag),aspect='auto',extent=[(Mla.imag)[0,-1],(Mla.imag)[0,0],(Mla.real)[-1,0],(Mla.real)[0,0]])\n",
    "        plt.clim(0,1)\n",
    "        plt.colorbar()\n",
    "        plt.scatter(LIm[id[1]],LRe[id[0]],color='r')\n",
    "\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X= [-0.03078483-1.25337383e-19j] lambda= [0.1+1.e-20j]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f5afd731a90>"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAEJCAYAAAC5Tb0qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYsUlEQVR4nO3dfbBdVXnH8e/vXhJilFd5CwFqqClIrYhmEF9aRUgNqRra0RlQKdbYDDOi4ks7WGbU1ukMba1T7aDMFSJYGWhVlAwT5a0w1iKYgLwkRAmihUsiIYAgQsy95zz9Y++rh5tz79nn7nXOXTn5fZg9OfvlrP3cfQ/PXWfttdZWRGBmZvkZmu0AzMysPSdoM7NMOUGbmWXKCdrMLFNO0GZmmXKCNjPLVK0ELelASTdI2lz+e0CbY+ZJ+qGkuyVtlPT3dc5pZpY7SaslbZO0oWVbx3w5Wd0a9PnATRGxGLipXJ/sN8CbI+J44JXAMkkn1TyvmVnOLgOWTdpWJV8+T90EvQK4vHx9OXD65AOi8Ey5OqdcPDrGzAZWRHwPeGLS5o75crK9asZxaERsLQPaKumQdgdJGgbuAF4KXBQRt09VoKRVwCqAF87Xq4996dyaIZrZoPv5w2Nsf6KhOmW85eQXxuNPNCode8c9v9kI7GjZNBIRIx3eVilftuqYoCXdCBzWZtcFnd47ISIawCsl7Q98S9LLI2LDFMeOACMAS46fFz+87siqpzGzPdSJb3m4dhmPP9Hgh9cdVenY4QWbd0TEkton7aBjgo6IU6faJ+lRSQvKvwYLgG0dyvqlpFso2mbaJmgzs9kQQJNmL0/RVb6E+m3Qa4Czy9dnA9dMPkDSwWXNGUkvAE4FflzzvGZmSQXBWDQqLTPUMV9OVjdBXwgslbQZWFquI+lwSWvLYxYAN0u6B1gH3BAR19Y8r5lZcs2K/3Ui6UrgB8AxkkYlrWSKfDmdWjcJI+Jx4JQ227cAy8vX9wAn1DmPmVmvBUEj0fTLEXHmFLt2yZfTqduLw8xsYDQz6wHsBG1mRnGTsOEEbWaWJ9egzcwyFMBYZo8AdII2M6O8SegatJlZhgIaeeVnJ2gzM5gYSZgXJ2gzMwBEg1rzLSXnBG1mxsRNQidoM7PsFP2gnaDNzLLUdA3azCw/rkGbmWUqEI3aE3ym5QRtZlZyE4eZWYYCsTOGZzuM53GCNjNjYqCKmzjMzLLkm4RmZhmKEI1wDdrMLEtN16DNzPJT3CTMKyXmFY2Z2SzxTUIzs4w13A/azCw/HkloZpaxpntxmJnlp5gsyQnazCw7gRjzUG8zs/xE4IEqZmZ5UnYDVWr9uZB0oKQbJG0u/z2gzTFHSrpZ0iZJGyV9uM45zcx6IShq0FWWfql7pvOBmyJiMXBTuT7ZOPCxiHgZcBLwAUnH1TyvmVlyDYYqLf1S90wrgMvL15cDp08+ICK2RsSd5etfAZuAhTXPa2aWVCCaUW3pl7pt0IdGxFYoErGkQ6Y7WNJLgBOA26c5ZhWwCuCohW4iN7P+CGBsd5uLQ9KNwGFtdl3QzYkkvQj4JnBeRDw91XERMQKMACw5fl50cw4zs5nT7jcfdEScOtU+SY9KWlDWnhcA26Y4bg5Fcr4iIq6ecbRmZj0S5DeSsG40a4Czy9dnA9dMPkCSgEuBTRHxuZrnMzPrmUZZi+609EvdBH0hsFTSZmBpuY6kwyWtLY95PXAW8GZJd5XL8prnNTNLKkI0Y6jSUoWkj5RdizdIulLSvG5jqtUiHhGPA6e02b4FWF6+/j5k1rBjZjZJcZMwzVBvSQuBDwHHRcRzkv4LOAO4rJty8rplaWY2a5I/k3Av4AWSxoD5wJaZFGBmtscrbhJW/rJ/kKT1LesjZQ+0oqyIRyR9FngIeA64PiKu7zYmJ2gzs1IXowS3R8SSqXaW016sABYBvwS+Luk9EfG1buLJq0+JmdksSTyS8FTgZxHxWESMAVcDr+s2JtegzcxKCR8a+xBwkqT5FE0cpwDrp3/Lrpygzcwo5oMea6ZJ0BFxu6RvAHdSTBj3I8oR0t1wgjYzY6KJI12rb0R8CvhUnTKcoM3MSrvdXBxmZnuCLrvZ9YUTtJkZQOImjhScoM3MSrk9k9AJ2syMiV4caebiSMUJ2syM3w1UyYkTtJlZyU0cZmYZci8OM7OMuReHmVmGIsS4E7SZWZ7cxGFmliG3QZuZZcwJ2swsQ+4HbWaWMfeDNjPLUASMJ5qwPxUnaDOzkps4zMwy5DZoM7OMhRO0mVmecrtJWKtFXNKBkm6QtLn894ApjlstaZukDXXOZ2bWKxFFG3SVpV/q3rI8H7gpIhYDN5Xr7VwGLKt5LjOzHhKN5lClpV/qnmkFcHn5+nLg9HYHRcT3gCdqnsvMrKciVGnpl7pt0IdGxFaAiNgq6ZC6AUlaBawCOGqhm8jNrD92y7k4JN0IHNZm1wXpw4GIGAFGAJYcPy96cQ4zs11E0Q6dk44JOiJOnWqfpEclLShrzwuAbUmjMzPro4HqxQGsAc4uX58NXFOzPDOzWREDeJPwQmCppM3A0nIdSYdLWjtxkKQrgR8Ax0galbSy5nnNzJKLqLb0S627cBHxOHBKm+1bgOUt62fWOY+ZWT94JKGZWYaK2rETtJlZlna7bnZmZnuK3a6bnZnZniAQTU/Yb2aWp8wq0LW72ZmZDYZIOxeHpP0lfUPSjyVtkvTabkNyDdrMbELaKvTnge9GxDskzQXmd1uAE7SZWSlVNztJ+wJ/Ary3KDd2Aju7LcdNHGZmlLPZNVVpAQ6StL5lWTWpuKOBx4CvSPqRpEskvbDbmJygzcygyNChagtsj4glLcvIpNL2Al4FfCkiTgB+zdQPNJmSE7SZWSnhXByjwGhE3F6uf4MiYXfFCdrMbEJUXDoVE/EL4GFJx5SbTgHu6zYc3yQ0MwMg+eOsPghcUfbgeBD4q24LcII2M5uQsJtdRNwFLKlThhO0mRkUA1WanizJzCxTTtBmZnnKbDIOJ2gzswlO0GZmGZoYqJIRJ2gzs5In7Dczy5V7cZiZ5UmuQZuZZajiMO5+coI2MwNAvkloZpYt16DNzDLVnO0Ans8J2swM3A+6W/c+dRCLrv3r2Q7DzDL3i6e+kKScgezFIWkZxRNsh4FLIuLCSftV7l8OPAu8NyLu7FTuH+23nR++9cspQjSzAXbiv29PU1BmCbr2E1UkDQMXAacBxwFnSjpu0mGnAYvLZRXwpbrnNTMbdClq0CcCD0TEgwCSrgJW8PzHu6wAvhoRAdwmaX9JCyJi63QFPxdNNu58LkGItruaoybDBHsL5krsrSHmaS/21pzKZYxFg9/EGDuiwY4IxgLGEI0Qzcyml7SZeS7S3N0bxCaOhcDDLeujwGsqHLMQ2CVBl48vXwVw4OF7c9tzixKEaLujYTUZosmwgjkaZ57GmKsG8zRWro//NoEPEQwraJQ3ecYYYiyGGIthdjLMjuY8xhhmR3NusS2GaTJEI/xYzkHw6+ZT9QsJBnKod7ufaPLfoSrHFBuLx5ePAOx37KHxrUdPqBedDaShLqs6zczuzltaT47fn6agAaxBjwJHtqwfAWyZwTG7GGsMseXpfWsHaGaDbayR5pvQIDZxrAMWS1oEPAKcAbxr0jFrgHPL9unXAE91an8GaIwP8+T2fRKEaGaDrDE+nKagQUvQETEu6VzgOopudqsjYqOkc8r9FwNrKbrYPUDRza7a48cbYujprLtqm1kOGomasAYtQQNExFqKJNy67eKW1wF8oNty1YC9nnHboZlNT40EZcRgNnH0jJow92knaDObnlLNoTGAvTh6Rk0Y3jHbUZhZ7lIlaNegu6AmzHkmsytmZtlJVoPOLN1knaBpwvDOzK6YmeUnRYJ2G3R31IQ5z2Z2xcwsO65BzwJFMDSW2RUzs+wo0uSJZIk+kawTNE0Yfi5B/xkzG2yZJdZUsk7QAobGXYM2s+kl6xyXWbrJOkHTDIZ2ugZtZh00E2RW3yTsUjhBm1kFidqgXYPuggI0PqCNS2aWTLKarxN0FyLQmGvQZtZBghq0cC+O7o07QZtZHyRugy6f17oeeCQi3jqTMvJO0BFobHy2ozCz3OXZBv1hYBMw46eOZJ6ggWZm3znMLD+ZtUFLOgL4M+AfgY/OtJzME3S4icPMOks1krB6MQdJWt+yPlI+T3XCvwF/C9R6JFTeCZqAhhO0mXXS9yaO7RGxpN0OSW8FtkXEHZLeVCecvBN0QLiJw8w6SZGfI1kvjtcDb5e0HJgH7CvpaxHxnm4LyjtBuwZtZpXkc5MwIj4BfAKgrEF/fCbJGXJP0AE0XIM2sw5S3dzzQJVuuAZtZlXkU4N+XnERtwC3zPT9eSfogEjVv9HMBleiNmgP9e6Wa9Bm1gfCTRzdifBAFTPrrP/9oPsi6wQduInDzDpLliUySzdZJ2ggzUTcZmZVZJZu8k7QEYTboM2skxTftAf1iSqSlgGfB4aBSyLiwkn7jwW+ArwKuCAiPlu58HAbtJn1yaAl6HLO04uApcAosE7Smoi4r+WwJ4APAad3fQK3Qe8xTo6HWMkGDuZZHmM+l/JybtZRsx2W7UFym7B/KEEZJwIPRMSDEbETuApY0XpARGyLiHXAWILz2QA6OR7io9zBoTzLEHAoz/JR7uDkeGi2Q7M9iKLa0i8pEvRC4OGW9dFym1llK9nAPJ5/v2EeDVayYZYisj1OdLH0SYo2aLXZNuMfQdIqYBXAPOaD2hVvg+bgeLb9dp71Z8A6y2zC/lRSJOhR4MiW9SOALTMtrJz0egRgXx0YKEUl33L3GPM5lF2T9GPMx58B64dBHUm4DlgsaRHwCHAG8K4E5YJAQ6497QlWxyv4SHPd85o5djDM6qFX+DNgnSW6uafMxl3UTtARMS7pXOA6im52qyNio6Rzyv0XSzqM4um2+wJNSecBx0XE09OXLtee9hA373U0NMT7Gnf/thfH6uHjuXl40WyHZruFBH/EB3WypIhYC6ydtO3ilte/oGj6MJvSzcOLnJBtVg1iE0dv+eutmfWLE3R1AuQ7+GbWQaos4Rq0mVmunKC7IMGQbxKaWQcpvmmne6p3MnknaDOzPhnUftC95TZoM+uXzCZnyz9Bm5n1iWvQ3XIbtJn1w6AOVDEzGwS+SdgNuR+0mVWQKE04QZuZ5SjwTUIzs1z5JmFX5Lk4zKyCRHnCCdrMLD8eqGJmlquIwZuwv+c8Yb+Z9Ute+Xk3SNBmZn3iJg4zsxwF4CaOLgj34jCzzlKlibzyM27gNTMrKaotHcuRjpR0s6RNkjZK+vBM4sm7Bm1m1kcJe3GMAx+LiDsl7QPcIemGiLivm0JcgzYzg9/NZldl6VRUxNaIuLN8/StgE7Cw25Dyr0F7siQz64NioErlGvRBkta3rI9ExEjbcqWXACcAt3cbU/4J2sysX6rPZrc9IpZ0OkjSi4BvAudFxNPdhuMEbWZW6qIG3bksaQ5Fcr4iIq6eSRlO0GZmkPSJKiomsr8U2BQRn5tpOZknaHmot5lVkOJeVdK5OF4PnAXcK+muctvfRcTabgrJPEGbmfVRoiaOiPg+Cf5qJKmeSlom6SeSHpB0fpv975Z0T7ncKun4FOc1M0smikdeVVn6pXYNWtIwcBGwFBgF1klaM6lD9s+AN0bEk5JOA0aA19Q9t5lZUpk98ipFDfpE4IGIeDAidgJXAStaD4iIWyPiyXL1NuCIBOc1M0sr0UCVVFK0QS8EHm5ZH2X62vFK4DtT7ZS0ClgFMG/oRZ4sycz6Rs28HuudIkG3y6Bt/8ZIOpkiQb9hqsLK0TgjAPvNOSSv7xtmNriCbgaq9EWKBD0KHNmyfgSwZfJBkl4BXAKcFhGPJzivmVkyIpIOVEkhRRv0OmCxpEWS5gJnAGtaD5B0FHA1cFZE3J/gnGZm6UVUW/qkdg06IsYlnQtcBwwDqyNio6Rzyv0XA58EXgx8sRhgw3iVcexmZn2VWQ06yUCVcnTM2knbLm55/X7g/SnOZWbWEwPaBt07Anm6UTPrJFGaGMReHGZmA6C/7ctVOEGbmUE5CMUJ2swsT3m1cDhBm5lNyK0ftBO0mdkEJ2gzswxFQCOvNo78E7S72ZlZv7gGbWaWKSdoM7MMBZDumYRJOEGbmQHFQBW3QZuZ5SfwTUIzs2y5DdrMLFNO0GZmOfJkSWZmeQrA0412QzCU4qlcZjbYEg1ocw3azCxHHuptZpangHA/aDOzTHkkoZlZptwGbWaWoQj34jAzy5Zr0GZmOQqi0ZjtIJ7HCdrMDDzdqJlZ1jLrZudhemZmFBXoaEalpQpJyyT9RNIDks6fSUz516D9TEIz64dIN2G/pGHgImApMAqsk7QmIu7rppz8E7SZWZ8kvEl4IvBARDwIIOkqYAXQVYJWZNatpJWkx4D/S1zsQcD2xGXOlGNpL6dYIK94HEt7x0TEPnUKkPRdip+pinnAjpb1kYgYaSnrHcCyiHh/uX4W8JqIOLebmLKuQUfEwanLlLQ+IpakLncmHEt7OcUCecXjWNqTtL5uGRGxLEUspXZts13Xhn2T0MwsvVHgyJb1I4At3RbiBG1mlt46YLGkRZLmAmcAa7otJOsmjh4Z6XxI3ziW9nKKBfKKx7G0l1MsRMS4pHOB64BhYHVEbOy2nKxvEpqZ7cncxGFmliknaDOzTA1Mgu40rFLS30i6q1w2SGpIOrDc93NJ95b7anfXqRDLmyQ91RLPJ6u+t0fxvFvSPeVyq6TjW/b1+9pI0hfK/fdIelXV99aI6UBJN0jaXP57QJtjjpR0s6RNkjZK+nDLvk9LeqTl97m81/GUx7X93VR9f6pYJB3T8rPfJelpSeeV+5JdG0nvLK99U9KU3fum+pykvC59ExG7/ULRCP9T4GhgLnA3cNw0x78N+O+W9Z8DB/UrFuBNwLV1f46E8bwOOKB8fRpw+yxem+XAdyj6kZ40EUsvrk3LOf8ZOL98fT7wT22OWQC8qny9D3D/xPmBTwMfT/h57hjPdL+bqu9PGcuk3/EvgN9LfW2AlwHHALcAS7r9jKW8Lv1aBqUG/dthlRGxE5gYVjmVM4ErM4kl1XtnXGZE3BoRT5art1H02eyFKj/fCuCrUbgN2F/SgorvnakVwOXl68uB0ycfEBFbI+LO8vWvgE3AwkTn7zqeHr+/TlmnAD+NiNQjgImITRHxkw6HTfc5SXld+mJQEvRC4OGW9VGm+J9H0nxgGfDNls0BXC/pDkmr+hTLayXdLek7kv6wy/f2Ip4JKylqsBP6fW2mOqYX12bCoRGxFYpEDBwy3cGSXgKcANzesvncsklmdYKvzlXjmep309XPkyiWCWewa+Un5bXpZLrPScrr0heD0g+6m2GVbwP+NyKeaNn2+ojYIukQ4AZJP46I7/UwljspvgI+U7bJfRtYXPG9vYinOFA6mSJBv6Flc7+vzVTH1Lo2km4EDmuz64KqZZTlvIjij/t5EfF0uflLwGfKeD4D/Cvwvj7Ek+R3k/DazAXeDnyiZXNX12a6WCLimiphtNm22/YlHpQE3c2wyl3+wkfElvLfbZK+RfE1aaZJqGMsLf9jExFrJX1R0kFd/hzJ4gGQ9ArgEuC0iHi8Jb6+Xptpjplb5eeYSkScOtU+SY9KWhARW8vmlG1THDeHIjlfERFXt5T9aMsxXwau7Uc80/xuKr0/ZSyl04A7W69Ht9dmulgqmu4z1tV1ycGgNHFUGlYpaT/gjcA1LdteKGmfidfAnwIbehmLpMOkYqJrSSdS/B4er/pz9CCeo4CrgbMi4v6W7X2/NuX6X6pwEvBU+XW0F9em9Zxnl6/PpuXzMaH8fV0KbIqIz03at6Bl9c+pd42qxjPd76bj+1PG0mKXezs9uDadTPc5SXld+mO271KmWiju/t9PcQf3gnLbOcA5Lce8F7hq0vuOprjTezewceK9vYwFOLc8190UN+VeN917+xDPJcCTwF3lsn4Wr40oJjr/KXAvLXfre3FtynJfDNwEbC7/PbDcfjiwtnz9Boqvyve0XKfl5b7/KGO9hyIJLOhDPFP+bqZ6f69iKdfnU1Qy9pv0/mTXhiLBjwK/AR4Frpsilrafk5TXpV+Lh3qbmWVqUJo4zMwGjhO0mVmmnKDNzDLlBG1mliknaDPrmXL04DZJSbrXSTpK0vUqJq26rxzVObCcoM2sly6jmFohla8C/xIRL6MYmJP9YJM6nKDNrGeiGHreOq0Ckn5f0nfLOUT+R9KxVcqSdBywV0TcUJb9TEQ8mz7qfDhBm1m/jQAfjIhXAx8HvljxfX8A/FLS1ZJ+JOlfJA33LMoMDMpcHGa2Gygnm3od8PVytgOAvct9fwH8Q5u3PRIRb6HIV39MMZPgQ8B/UowOvrS3Uc8eJ2gz66ch4JcR8crJO6KYgOrqXd7xO6PAjyLiQQBJ36Z4qMPAJmg3cZhZ30Qxk+PPJL0TfvuIs+M7vG3COuAASQeX628G7utBmNlwgjaznpF0JfAD4BhJo5JWAu8GVkqamOip0pNxIqJB0WZ9k6R7KSbW+nJvIs+DJ0syM8uUa9BmZplygjYzy5QTtJlZppygzcwy5QRtZpYpJ2gzs0w5QZuZZer/AdsIdcojURD5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# with astrocytes\n",
    "Ne = 4000; Ni = 1000; Na = 2000; Ca = 0.1\n",
    "Vth = 20 ; taub = 1; tauba = 1; taue=20; taui=10; tau_astro = 160; Jee = 280/Ne*taue; Jei=70/Ni*taue; Jie=500/Ne*taui; Jii=100/Ni*taui; \n",
    "Cae = 0.5; Cai = 0.5\n",
    "Jia = 0.4*11/Na*10*20*taui; Jaa = 2/Na*tau_astro; Jea = 1*11/Na*10*20*taue; \n",
    "Jai = 3*0.2/22*100/Ni/1.5*(1/Cai)/10*tau_astro; Jae = 0.2/8*400/Ne/1.5*(1/Cae)/10*tau_astro\n",
    "   \n",
    "Vli=6.5; Vle=7.6; Vr=14; Vla = 7; Vra = 9; Vtha = 13\n",
    "\n",
    "siga = 3; delaye = 0.5; delayi = 0.25; delaya = 1e3\n",
    "taude = 23; taure = 8; taudi = 1; tauri = 1; taura = 8; tauda = 2\n",
    "beta = 1\n",
    "Ka = 600\n",
    "\n",
    "sig = 3; sige = sig; sigi = sig;\n",
    "re0 = 0.00121; ri0 = 0.000546225; ra0 = 0.00027 #unstable\n",
    "LRe = np.concatenate([np.arange(-3e-1,-1e-13,5e-3),np.arange(1e-13,3e-1,5e-3)])\n",
    "\n",
    "# re0 = 0.00294; ri0 = 0.0064484; ra0 = 0.00041 #stable up\n",
    "# LRe = np.concatenate([np.arange(-2e-2,-1e-13,5e-4),np.arange(1e-13,2e-2,5e-4)])\n",
    "\n",
    "\n",
    "Ie0 =  Vle + taub*(Ne*Jee*re0-Ni*Jei*ri0) -Ka*beta*re0 + tauba*Ca*Na*Jea*ra0\n",
    "yte = (Vth-Ie0)/sige; yre = (Vr-Ie0)/sige\n",
    "Ii0 =  Vli + taub*(Ne*Jie*re0-Ni*Jii*ri0) + tauba*Ca*Na*Jia*ra0\n",
    "yti = (Vth-Ii0)/sigi; yri = (Vr-Ii0)/sigi\n",
    "Ia0 =  Vla + taub*(Cae*Ne*Jae*re0-Cai*Ni*Jai*ri0) + tauba*Na*Jaa*ra0\n",
    "yta = (Vtha-Ia0)/siga; yra = (Vra-Ia0)/siga\n",
    "\n",
    "LIm = np.concatenate([np.arange(-1e-6,-1e-20,1e-7),np.arange(1e-20,1e-6,1e-7)])\n",
    "\n",
    "X = np.zeros((len(LRe),len(LIm)),dtype='complex')\n",
    "Mla = np.zeros((len(LRe),len(LIm)),dtype='complex')\n",
    "for i in range(len(LRe)):\n",
    "    for j in range(len(LIm)):\n",
    "        la = complex(LRe[i],LIm[j])\n",
    "        Mla[i,j] = la\n",
    "        Fee = taub*Ne*Jee*R(la*taue,yre,yte,re0,sige,taue)*S(la,delaye,taure,taude)\n",
    "        Fii = taub*Ni*Jii*R(la*taui,yri,yti,ri0,sigi,taui)*S(la,delayi,tauri,taudi)\n",
    "        Faa = tauba*Ca*Jaa*R(la*tau_astro,yra,yta,ra0,siga,tau_astro)*S(la,delaya,taura,tauda)\n",
    "        Fie = taub*Ne*Jie*R(la*taui,yri,yti,ri0,sigi,taui)*S(la,delaye,taure,taude)\n",
    "        Fei = taub*Ni*Jei*R(la*taue,yre,yte,re0,sige,taue)*S(la,delayi,tauri,taudi)\n",
    "        Fea = tauba*Ca*Jea*R(la*taue,yre,yte,re0,sige,taue)*S(la,delaya,taura,tauda)\n",
    "        Fia = tauba*Ca*Jia*R(la*taui,yri,yti,ri0,sigi,taui)*S(la,delaya,taura,tauda)\n",
    "        Fae = taub*Ne*Jae*R(la*tau_astro,yra,yta,ra0,siga,tau_astro)*S(la,delaye,taure,taude)\n",
    "        Fai = taub*Ni*Jai*R(la*tau_astro,yra,yta,ra0,siga,tau_astro)*S(la,delayi,tauri,taudi)\n",
    "        \n",
    "        X[i,j] = (Fee-1)*(Fii-1)*(Faa-1)+ Fei*Fia*Fae + Fea*Fie*Fai - (Fee-1)*Fia*Fai - Fei*Fie*(Faa-1) - Fea*(Fii-1)*Fae\n",
    "        \n",
    "eps = 1\n",
    "id = np.where(np.abs(X.real)+np.abs(X.imag)==np.min(np.abs(X.real)+np.abs(X.imag)))\n",
    "print('X=',X[id],'lambda=',Mla[id])\n",
    "plt.figure()\n",
    "plt.imshow(np.abs(X.real)+np.abs(X.imag),aspect='auto',extent=[(Mla.imag)[0,-1],(Mla.imag)[0,0],(Mla.real)[-1,0],(Mla.real)[0,0]])\n",
    "plt.clim(0,10)\n",
    "plt.colorbar()\n",
    "plt.scatter(LIm[id[1]],LRe[id[0]],color='r')\n"
   ]
  },
  {
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}