{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Incoherent Type 4 FFL with AND gate\n",
    "\n",
    "Implement the incoherent type 4 FFL with an AND gate, as described by Mangan & Alon (PNAS, 2003).\n",
    "\n",
    "The parameters are\n",
    "- B_y: basal rate of Y's transcription\n",
    "- Beta_y: maximum rate of Y's transcription\n",
    "- B_z: basal rate of Z's transcription\n",
    "- Beta_z: maximum rate of Z's transcription\n",
    "- K_xy: the activation threshold for X's activation of Y\n",
    "- K_xz: the activation threshold for X's activation of Z\n",
    "- K_yz: the activation threshold for Y's activation of Z\n",
    "- H: the Hill coefficient (exponent) for all regulation\n",
    "- alpha_y: the degradation rate constant for Y\n",
    "- alpha_z: the degradation rate constant for Z\n",
    "- SX: amount of active X in system (X*=SX)\n",
    "- SY: whether or not Y is activated (Y*=Y if SY=1, Y*=0 if SY=0)\n",
    "\n",
    "The model is\n",
    "\\begin{eqnarray*}\n",
    "\\frac{dY}{dt} &=& B_y + \\beta_Y \\frac{1}{1+(X^*/K_{xy})^H} - \\alpha_y Y\\\\\n",
    "\\frac{dZ}{dt} &=& B_z + \\beta_Z \\left(\\frac{(X^*/K_{xz})^H}{1+(X^*/K_{xz})^H}\\right) \\left(\\frac{(Y^*/K_{yz})^H}{1+(Y^*/K_{yz})^H}\\right) - \\alpha_z Z\\\\\n",
    "\\end{eqnarray*}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def incoherent_type_4_and_model( t, statevars, params ):\n",
    "    Y = statevars[0];\n",
    "    Z = statevars[1];\n",
    "    B_y = params[0]\n",
    "    Beta_y = params[1]\n",
    "    B_z = params[2]\n",
    "    Beta_z = params[3]\n",
    "    K_xy = params[4]\n",
    "    K_xz = params[5]\n",
    "    K_yz = params[6]\n",
    "    H = params[7]\n",
    "    alpha_y = params[8]\n",
    "    alpha_z = params[9]\n",
    "    SX = params[10]\n",
    "    SY = params[11]\n",
    "    \n",
    "    if SX == 0:\n",
    "        Xstar = 0\n",
    "    else:\n",
    "        Xstar = 1\n",
    "    if SY == 0:\n",
    "        Ystar = 0\n",
    "    else:\n",
    "        Ystar = Y\n",
    "    \n",
    "    dydt = np.zeros( (2,) )\n",
    "    yreg = 1/(1 + (Xstar/K_xy)**H)\n",
    "    zregx = ((Xstar/K_xz)**H)/(1 + (Xstar/K_xz)**H)\n",
    "    zregy = ((Ystar/K_yz)**H)/(1 + (Ystar/K_yz)**H)\n",
    "    zreg = zregx * zregy\n",
    "    dydt[0] = B_y + Beta_y*yreg - alpha_y*Y;\n",
    "    dydt[1] = B_z + Beta_z*zreg - alpha_z*Z;\n",
    "\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate the model\n",
    "\n",
    "Simulate the model. In the paper, they specify that for most simulations,\n",
    "\n",
    "- B_y = 0\n",
    "- Beta_y = 1\n",
    "- B_z = 0\n",
    "- Beta_z = 1\n",
    "- H = 2\n",
    "- alpha_y = 1\n",
    "- alpha_z = 1\n",
    "\n",
    "The activation thresholds are given in the figure captions, as are SX and SY."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "B_y = 0\n",
    "Beta_y = 1\n",
    "B_z = 0\n",
    "Beta_z = 1\n",
    "K_xy = 0.1\n",
    "K_xz = 0.1\n",
    "K_yz = 1 # 1, 0.3, 0.1\n",
    "H = 2\n",
    "alpha_y = 1\n",
    "alpha_z = 1\n",
    "\n",
    "time_x_on = 0.5\n",
    "time_x_off = 2.5\n",
    "params_off = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,0,1);\n",
    "params_on = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,1,1);\n",
    "\n",
    "mf_off = lambda t,y : incoherent_type_4_and_model(t,y, params_off)\n",
    "mf_on = lambda t,y : incoherent_type_4_and_model(t,y, params_on)\n",
    "\n",
    "sol1 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(0,time_x_on), y0=(0,0), \n",
    "                                t_eval=np.arange(0,time_x_on,.01) )\n",
    "sol2 = scipy.integrate.solve_ivp( fun=mf_on, t_span=(time_x_on,time_x_off), y0=sol1.y[:,-1], \n",
    "                                t_eval=np.arange(time_x_on,time_x_off,.01) )\n",
    "sol3 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(time_x_off,6), y0=sol2.y[:,-1], \n",
    "                                t_eval=np.arange(time_x_off,6,.01) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the simulation\n",
    "Plot X over time in an upper subplot and Z over time in a lower subplot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqc0lEQVR4nO3deXhV5bn+8e+zdwYggTAkDCZAmBQBFSQOaLVWK+KpFdvaOnSwPf5qB23tcE6rPa1a23NO7aTtUVup2trJoaiVttRZq7aCBAEREIgRSsIUpkAIkOn5/bEXNoQAQfZirb25P9e1r73Xu4b9LAXurLXevK+5OyIiInGTiLoAERGRziigREQklhRQIiISSwooERGJJQWUiIjEUk7UBaRLcXGxl5eXR12GiIgcpLlz525w95KO7VkTUOXl5VRWVkZdhoiIHCQzW9lZe2i3+MzsXjNbb2av72O9mdlPzazKzF4zsxPbrbvCzJYHryvCqlFEROIrzGdQvwKm7Gf9+cCo4HUV8DMAM+sL3AicApwM3GhmfUKsU0REYii0W3zu/oKZle9nk6nArz01lMUsM+ttZoOAs4Cn3H0TgJk9RSro7g+r1ij8fvY/eWx+7WH9zqnjS7n8lCGH9TtFRN6pKHvxlQKr2i3XBG37at+LmV1lZpVmVllXVxdaoWF4bH4ti9dsPWzft3jN1sMeiCIihyKjO0m4+zRgGkBFRUXGDSo4ZlAvHvzMpMPyXZfc9fJh+R4RkXSJ8gqqFhjcbrksaNtXu4iIHEGiDKgZwCeC3nynAvXuvgZ4AphsZn2CzhGTgzYRETmChHaLz8zuJ9XhodjMakj1zMsFcPefAzOBfwOqgEbgU8G6TWb2HWBOcKibd3eYEBGRI0eYvfguO8B6B67ex7p7gXvDqEtERDKDxuITEZFYUkCJiEgsKaBERCSWFFAiIhJLCigREYklBZSIiMSSAkpERGJJASUiIrGkgBIRkVhSQImISCwpoEREJJYUUCIiEksKKBERiSUFlIiIxJICSkREYinUgDKzKWa21MyqzOy6Ttbfambzg9cyM9vSbl1ru3UzwqxTRETiJ8wZdZPAHcC5QA0wx8xmuPvi3du4+5fbbf8FYEK7Q+xw9/Fh1SciIvEW5hXUyUCVu1e7exPwADB1P9tfBtwfYj0iIpJBwgyoUmBVu+WaoG0vZjYUGAY82665m5lVmtksM7toH/tdFWxTWVdXl6ayRUQkDuLSSeJSYLq7t7ZrG+ruFcDlwG1mNqLjTu4+zd0r3L2ipKTkcNUqIiKHQZgBVQsMbrdcFrR15lI63N5z99rgvRp4nj2fT4mISJYLM6DmAKPMbJiZ5ZEKob1645nZaKAP8HK7tj5mlh98LgZOBxZ33FdERLJXaL343L3FzK4BngCSwL3uvsjMbgYq3X13WF0KPODu3m73Y4G7zKyNVIh+r33vPxERyX6hBRSAu88EZnZou6HD8k2d7PcP4LgwaxMRkXiLSycJERGRPSigREQklhRQIiISSwooERGJJQWUiIjEkgJKRERiSQElIiKxpIASEZFYUkCJiEgsKaBERCSWFFAiIhJLCigREYklBZSIiMSSAkpERGIp1IAysylmttTMqszsuk7Wf9LM6sxsfvD6f+3WXWFmy4PXFWHWKSIi8RPafFBmlgTuAM4FaoA5Zjajk4kHH3T3azrs2xe4EagAHJgb7Ls5rHpFRCRe9nkFZWZD9rPujC4c+2Sgyt2r3b0JeACY2sW6zgOecvdNQSg9BUzp4r4iIpIF9neL73kz+1pwJQSAmQ0ws98Ct3bh2KXAqnbLNUFbRx8ys9fMbLqZDT6Yfc3sKjOrNLPKurq6LpQkIiKZYn8BNREYAcw3s7PN7FrgFeBlUldH6fAnoNzdjyd1lXTfwezs7tPcvcLdK0pKStJUkoiIxME+n0EFt9Y+EwTT08Bq4FR3r+nisWuBwe2Wy4K29t+xsd3i3cD32+17Vod9n+/i94qISBbY3zOo3mZ2F/ApUs9/pgN/NbOzu3jsOcAoMxtmZnnApcCMDt8xqN3ihcCS4PMTwGQz62NmfYDJQZuIiBwh9teL71XgTuBqd28BnjSz8cCdZrbS3S/b34HdvcXMriEVLEngXndfZGY3A5XuPgP4opldCLQAm4BPBvtuMrPvkAo5gJvdfdM7PksREck4+wuoMzveznP3+cBpZvbprhzc3WcCMzu03dDu8/XA9fvY917g3q58j4iIZJ993uLb37Mmd/9FOOWIiIikaKgjERGJJQWUiIjEkgJKRERiSQElIiKxpIASEZFYUkCJiEgsKaBERCSWFFAiIhJLCigREYklBZSIiMSSAkpERGJJASUiIrGkgBIRkVhSQImISCyFGlBmNsXMlppZlZld18n6r5jZYjN7zcyeMbOh7da1mtn84DWj474iIpLd9jdh4SExsyRwB3AuUAPMMbMZ7r643WbzgAp3bzSzzwHfBy4J1u1w9/Fh1SciIvEW5hXUyUCVu1e7exPwADC1/Qbu/py7NwaLs4CyEOsREZEMEmZAlQKr2i3XBG37ciXw13bL3cys0sxmmdlFne1gZlcF21TW1dUdcsEiIhIfod3iOxhm9jGgAnh3u+ah7l5rZsOBZ81sobu/2X4/d58GTAOoqKjww1awiIiELswrqFpgcLvlsqBtD2b2XuC/gAvdfdfudnevDd6rgeeBCSHWKiIiMRNmQM0BRpnZMDPLAy4F9uiNZ2YTgLtIhdP6du19zCw/+FwMnA6071whIiJZLrRbfO7eYmbXAE8ASeBed19kZjcDle4+A/gBUAj8wcwA/unuFwLHAneZWRupEP1eh95/IiKS5UJ9BuXuM4GZHdpuaPf5vfvY7x/AcWHWJiIi8aaRJEREJJYUUCIiEksKKBERiSUFlIiIxJICSkREYkkBJSIisaSAEhGRWFJAiYhILCmgREQklhRQIiISSwooERGJJQWUiIjEkgJKRERiSQElIiKxFGpAmdkUM1tqZlVmdl0n6/PN7MFg/WwzK2+37vqgfamZnRdmnSIiEj+hBZSZJYE7gPOBMcBlZjamw2ZXApvdfSRwK3BLsO8YUjPwjgWmAHcGxxMRkSNEmBMWngxUuXs1gJk9AExlz6nbpwI3BZ+nA7dbamrdqcAD7r4LeMvMqoLjvRxWsd/+0yIWr94a1uH3snjNVsYM6nXYvm/3d15yV2j/CUXkCDTmqF7c+P6xoRw7zIAqBVa1W64BTtnXNsEU8fVAv6B9Vod9Szt+gZldBVwFMGTIkLQVfjiMGdSLqeP3OqXQHM7vEhFJh1CnfA+bu08DpgFUVFT4oRwrrJ8A4uLyU4Zw+SmZFeIicmQLs5NELTC43XJZ0NbpNmaWAxQBG7u4r4iIZLEwA2oOMMrMhplZHqlODzM6bDMDuCL4fDHwrLt70H5p0MtvGDAKeCXEWkVEJGZCu8UXPFO6BngCSAL3uvsiM7sZqHT3GcA9wG+CThCbSIUYwXYPkepQ0QJc7e6tYdUqIiLxY6kLlsxnZnXAykM8TDGwIQ3lxJXOL3Nl87mBzi/THer5DXX3ko6NWRNQ6WBmle5eEXUdYdH5Za5sPjfQ+WW6sM5PQx2JiEgsKaBERCSWFFB7mhZ1ASHT+WWubD430PllulDOT8+gREQklnQFJSIisaSAEhGRWFJABQ40d1UmM7N7zWy9mb0edS3pZmaDzew5M1tsZovM7Nqoa0onM+tmZq+Y2YLg/L4ddU3pZmZJM5tnZn+OupZ0M7MVZrbQzOabWWXU9aSbmfU2s+lm9oaZLTGzSWk9vp5BvT131TLgXFIjp88BLnP3xfvdMUOY2ZlAA/Brdx8XdT3pZGaDgEHu/qqZ9QTmAhdl0f87AwrcvcHMcoGXgGvdfdYBds0YZvYVoALo5e4XRF1POpnZCqDC3bPyl3TN7D7gRXe/OxjSroe7b0nX8XUFlfL23FXu3gTsnrsqK7j7C6SGkso67r7G3V8NPm8DltDJ1CyZylMagsXc4JU1P1WaWRnwPuDuqGuRg2NmRcCZpIasw92b0hlOoIDarbO5q7LmH7kjhZmVAxOA2RGXklbBLbD5wHrgKXfPpvO7Dfga0BZxHWFx4EkzmxvMX5dNhgF1wC+DW7R3m1lBOr9AASVZwcwKgYeBL7n74Zsa+TBw91Z3H09q2pmTzSwrbtOa2QXAenefG3UtIXqXu58InA9cHdxuzxY5wInAz9x9ArAdSOvzewVUiuafymDBs5mHgd+5+yNR1xOW4PbJc8CUiEtJl9OBC4PnNA8AZ5vZb6MtKb3cvTZ4Xw88SupxQraoAWraXdFPJxVYaaOASunK3FUSQ0EngnuAJe7+46jrSTczKzGz3sHn7qQ68rwRaVFp4u7Xu3uZu5eT+jv3rLt/LOKy0sbMCoKOOwS3viYDWdOT1t3XAqvM7Jig6RxSUySlTUZP+Z4u+5q7KuKy0sbM7gfOAorNrAa40d3vibaqtDkd+DiwMHhOA/ANd58ZXUlpNQi4L+hpmgAecves646dpQYAj6Z+hiIH+L27Px5tSWn3BeB3wQ/21cCn0nlwdTMXEZFY0i0+ERGJJQWUiIjEkgJKRERiSQElIiKxpIASEZFYUkCJiEgsKaBERCSWFFAiIhJLCigREYklBZSIiMSSAkpERGIpawaLLS4u9vLy8qjLEBGRgzR37twN7l7SsT1rAqq8vJzKysqoyxARkYNkZis7a9ctPhERiSUF1BHA3Vm+fhv1jc1RlyIi0mUKqCPA/z1bxeQfv8B/Tl8QdSkiIl2mgMpyTy1ex4+fWsaoAYU8uXgds6o3Rl2SiEiXKKCy2KpNjXzlofkcV1rEg1dNoqxPd775x9dpammLujQRkQNSQGWpltY2vvTgfADu/OiJ9CnI49sXjqVqfQP3vPRWtMWJiHSBAipL3f3SW8xduZnvXjSOwX17AHDOsQOYPGYAP31mOTWbGyOuUERk/0INKDObYmZLzazKzK7rZH2+mT0YrJ9tZuVB+8lmNj94LTCzD4RZZ7ZZuXE7tz61jMljBnDhCUftse7GC8cCcNOMxVGUJiLSZaEFlJklgTuA84ExwGVmNqbDZlcCm919JHArcEvQ/jpQ4e7jgSnAXWaWNb9UHLabZiwiN5ng5qnjMLM91pX27s617x3F00vW8eSitRFVKCJyYGFeQZ0MVLl7tbs3AQ8AUztsMxW4L/g8HTjHzMzdG929JWjvBniIdWaV595Yz3NL67j2nFEMLOrW6TZXvmsYxwzoyQ2PLWLrTv1ulIjEU5gBVQqsardcE7R1uk0QSPVAPwAzO8XMFgELgc+2C6y3mdlVZlZpZpV1dXUhnEJmaWlt439mLmFYcQFXnFa+z+1ykwluufh41m/byf/OXHL4ChQROQix7STh7rPdfSxwEnC9me11OeDu09y9wt0rSkr2GmfwiPPovFqWr2/ga+cdQ17O/v/Xjh/cm0+fMZz7X1nFS8s3HKYKRUS6LsyAqgUGt1suC9o63SZ4xlQE7PGbpO6+BGgAxoVWaRZobm3jp88u57jSIqaMG9ilfb587tEMKy7gukdeY/uuvS5QRUQiFWZAzQFGmdkwM8sDLgVmdNhmBnBF8Pli4Fl392CfHAAzGwqMBlaEWGvGe3ReLas27eDL547aq2PEvnTLTfL9i4+ndssOfvDE0pArFBE5OKEFVPDM6BrgCWAJ8JC7LzKzm83swmCze4B+ZlYFfAXY3RX9XcACM5sPPAp83t11H2ofWtucnz//JuNKe/GeY/of1L4nlffliknl/OofK3jlrU0hVSgicvDMPTs6yFVUVPiROh/U46+v4bO/fZXbL5/ABccfdeAdOti+q4Xzf/Iibe789doz6NktN4QqRUQ6Z2Zz3b2iY3tsO0lI19394lsM7tud88cNekf7F+TncOslJ7B6yw79Aq+IxIYCKsMtrKmncuVmrphUTjLRtWdPnZk4tC/XvGckD79aw8yFa9JYoYjIO6OAynC/fnkFPfKSfOSkwQfe+AC+cM4oTigr4huPLmRt/c40VCci8s4poDJYfWMzMxasZur4Unql4blRbjLBrZeMZ1dzG/85fQFtbdnxfFJEMpMCKoP9cX4tu1ra+OgpQ9J2zOElhXzrgjG8uHwDd79UnbbjiogcLAVUBntgzirGlfZiXGlRWo972cmDOX/cQG55fClzV6rruYhEQwGVoRatrmfJmq18pOLQnz11ZGbccvHxlPbuzjW/n8em7U1p/w4RkQNRQGWoR16tJS+Z2Gu+p3Tp1S2XOz96IhsbmvjqQ/P1PEpEDjsFVAZqaW1jxoLVnD26P7175IX2PeNKi/jWBcfy3NI67npBz6NE5PBSQGWgWdWbqNu2i6njw7l6au9jpw7lguMH8cMnl/LymxsPvIOISJoooDLQnxaspjA/h/eMPrhx994JM+N/P3gcw4oLuPr3r1KzuTH07xQRAQVUxmlubePxRWs5d8wAuuUmD8t39uyWy7SPT6S5tY2rfj2XHU2th+V7ReTIpoDKMP94cyP1O5r5t+Pe2bh779TwkkL+77IJLFm7lf+cvoBsGWRYROJLAZVhHn99DQV5Sc4YVXzYv/usY/rztfNG8+fX1vDzv6nThIiEKyfqAqTrWtucpxav46zR/Q/b7b2OPvvu4Sxes5XvP/EGI0oKmDy2a7P3iogcLF1BZZD5qzazoaGJ8yIMBTPj+x86nuPLevPFB+axYNWWyGoRkeymgMogTy1eT07COOuYkkjr6J6X5O5PVFDSM58r75vDqk3q2Sci6aeAyiDPLFnHycP6pmXk8kNV0jOfX37yZJpa2vjUr+ZQ39gcdUkikmUUUBli1aZGlq9v4OzD8LtPXTWyfyHTPlHByo3b+cxvK9nZrO7nIpI+CqgM8fzS9QCxCiiAU4f34wcXn8Cs6k188f55tLS2RV2SiGSJUAPKzKaY2VIzqzKz6zpZn29mDwbrZ5tZedB+rpnNNbOFwfvZYdaZCZ5fWsfgvt0ZVlwQdSl7uWhCKTe+fwxPLl7HdY8s1MCyIpIWoXUzN7MkcAdwLlADzDGzGe6+uN1mVwKb3X2kmV0K3AJcAmwA3u/uq81sHPAEUBpWrXHX1NLGy9Ub+dCJZZhZ1OV06lOnD6N+RzO3Pb2cou65fPN9x8a2VhHJDGH+HtTJQJW7VwOY2QPAVKB9QE0Fbgo+TwduNzNz93nttlkEdDezfHffFWK9sTV35WYam1o58+hoe+8dyLXnjGJLYzP3vPQWvbrlcu17R0VdkohksDADqhRY1W65BjhlX9u4e4uZ1QP9SF1B7fYh4NXOwsnMrgKuAhgyJH3TnsfNi8vryEkYpw7vG3Up+2Vm3HDBGLbtbOHWp5eRkzSufs/IqMsSkQwV65EkzGwsqdt+kztb7+7TgGkAFRUVWfvg4+9VGxg/uDc9Y9C9/EASCeP7Fx9Pa1sbP3hiKe7ONWfrSkpEDl6YAVULtJ+PvCxo62ybGjPLAYqAjQBmVgY8CnzC3d8Msc5Yq9/RzMLa+oz6Rz6ZMH70kfEA/PDJZQAZVb+IxEOYATUHGGVmw0gF0aXA5R22mQFcAbwMXAw86+5uZr2BvwDXufvfQ6wx9mZXb6TN4bQR/aIu5aDsDikz44dPLsMdvnCOQkpEui60gAqeKV1DqgdeErjX3ReZ2c1ApbvPAO4BfmNmVcAmUiEGcA0wErjBzG4I2ia7+/qw6o2rl6s3kp+TYMKQ3lGXctCSCeOHHz4BA3701DJ2tbTx1clHq3efiHRJqM+g3H0mMLND2w3tPu8EPtzJft8FvhtmbZliVvUmJg7tQ35ONKOXH6pkwvjBh08gLyfB7c9Vsamxie9MHUcyoZASkf3TSBIxtqWxiTfWbuXU4Zl1e6+jZCI1bfznzhrB72f/ky/c/yq7WjQskojsX6x78R3p5qzYjDucMize3cu7wsz4+pTR9CvI47t/WUL9jjnc9fEKCvP1R1BEOqcrqBibs2ITeckEJwzuHXUpafP/zhjOjz6cGrvvsmmzWLd1Z9QliUhM7TegzOy8/azb69mRpNfstzZxwuCiyGbPDcuHJpbxi09M5M26Bi664+8sWl0fdUkiEkMHuoKaaWbPmVln4+BdH0ZBkrKjqZVFtfVUlGf+7b3OnD16AH/47CQAPvzzl3l68bqIKxKRuDlQQL0G/B6YZWYXd1inblghmr9qCy1tzknlfaIuJTRjjyrisatPZ2T/Qj79m0rufrEa96wdEEREDtKBAsrd/RfAOcDXzeyXZtZj97pwSzuyvfrPzQCcOCR7Awqgf69uPHjVJM4bM5Dv/mUJX/3DAk18KCJAFztJuPsyYBKwDphnZh0HfZU0m7tyMyP7F9K7R17UpYSue16SOz96IteeM4pHXq3lg3f+g39ubIy6LBGJ2IEC6u3beO7e4u7XAZ8B7gc0bk1I3J15/9zMiRk4esQ7lUgYXz73aO79ZAU1mxt5/+0v8dzSI27gEBFp50AB9e2ODe7+PDAR+O8wChJYsbGRzY3NTMjy23udOXv0AP70hXdxVO/u/Puv5nDrU8s0jbzIEWq/AeXuf9xH+2Z3/14oFQnzgudPmTj+XjoM7VfAI587jQ9MKOUnzyznsl/MonbLjqjLEpHDTL+oG0PzV22hIC/JqP49oy4lMt3zkvz4I+P58UdOYPHqrZx/2wvMXLgm6rJE5DBSQMXQglVbOK6sSAOqAh88sYyZ157BsJJCPv+7V7n+kdfYvqsl6rJE5DBQQMXMrpZWFq/ZmlXDGx2qof0KmP7ZSXzurBE8MGcV5932Av94c0PUZYlIyBRQMfPGmm00tzrjy3pHXUqs5CYTfH3KaB76zCRyEsblv5jNN/+4UFdTIllMARUzr9VsAeC4sqJoC4mpk8r78tdrz+TKdw3jd7P/mbqaqtLVlEg2UkDFzGs19fQtyKO0d/eoS4mt7nlJvnXBGP7wmUnkJhNcfvds/uMPC9jYsCvq0kQkjRRQMbOwtp7jSos0LXoXVJT3ZeYXz+BzZ43gj/NqOftHf+N3s1fS2qZRuESygQIqRnY2t7J8fQPHler2Xld1z0vy9Smj+eu1Z3DsoJ7816Ov88Gf/YOFNZrCQyTTKaBiZMmarbS2OeMUUAdt1ICe3P/pU7ntkvHUbt7BhXe8xNenv6YJEUUyWKgBZWZTzGypmVWZ2XWdrM83sweD9bPNrDxo7xfMQ9VgZreHWWOcvL56KwDjSntFXElmMjMumlDKM199N/9++jAemVfDWT94ntueXkZjk3r7iWSa0ALKzJLAHcD5wBjgMjMb02GzK4HN7j4SuBW4JWjfCXwL+I+w6oujRbX1FHXPVQeJQ1TUPZdvXTCGp7/ybt4zuoTbnl7OWT94nofmrNLzKZEMEuYV1MlAlbtXu3sT8AAwtcM2U4H7gs/TgXPMzNx9u7u/RCqojhiL12xl7FG91EEiTYb2K+DOj07k4c9NorRPd7728Guce+vfmLFgNW0KKpHYCzOgSoFV7ZZrgrZOt3H3FqAe6NfVLzCzq8ys0swq6+rqDrHcaLW0tvHG2m2MPUq399Jt4tC+PPK50/j5x04kN5Hgi/fPY8pPUmP7KahE4iujO0m4+zR3r3D3ipKSkqjLOSTVG7bT1NLG2KPUQSIMZsaUcYP467VncPvlE2hz+PzvXuXffvoif1qwWlN6iMRQmAFVCwxut1wWtHW6jZnlAEXAxhBriq3FQQeJYwfpCipMiYRxwfFH8cSXzuQnl46nqbWNL9w/j/f86Hl+/fIKdjRpunmRuAgzoOYAo8xsmJnlAZcCMzpsMwO4Ivh8MfCsux+R91yWrNlKXjLB8JKCqEs5IiQTxtTxpTz95Xdz18cnUlyYzw2PLeK07z3DbU8vY9P2pqhLFDni5YR1YHdvMbNrgCeAJHCvuy8ys5uBSnefAdwD/MbMqoBNpEIMADNbAfQC8szsImCyuy8Oq96oLV6zlVEDCslNZvRd14yTSBjnjR3I5DEDqFy5mbv+9ia3Pb2cnz3/Ju8/4Sg+MWkox2vgXpFIhBZQAO4+E5jZoe2Gdp93Ah/ex77lYdYWN2+s3caZozL7OVomMzNOKu/LSeV9WbZuG/f9YwWPzqtl+twaThjcm4+fOpQLjh9Et9xk1KWKHDH043oMbGzYRd22XRw76MidQTdOjh7Qk//+wHHM+sY53PT+MTTsbOY//rCASf/7DN/582KWrNkadYkiR4RQr6Cka5au3QbA6IHqIBEnvbrl8snTh3HFaeW8/OZGfjNrJb9+eQX3vPQWY4/qxcUTy5g6vpS+BXlRlyqSlRRQMfBGEFDHDNQVVByZGaeNLOa0kcVs3t7EjAWrmT63hm//aTH/M3MJZ4/uzwcmlHLWMf11C1AkjRRQMbBs3Tb6FeRR0jM/6lLkAPoU5HHFaeVccVo5b6zdysNza3h0Xi1PLFpHj7wk7z12AO87fhDvPrpEYSVyiBRQMfDG2m0cPUBXT5lm9MBe/Nf7xvD1KaOZ/dYm/vzaGh5/fQ0zFqymIC/Je8cM4PxxAzljVAkF+fqrJnKw9LcmYu7O8nXb+HDF4ANvLLGUk0xw+shiTh9ZzHemjuXl6o385bU1PL5oLY/NX01eMsGpI/rx3mP7c/bo/pT16RF1ySIZQQEVsdotO9je1MqoAYVRlyJpkJNMcMaoEs4YVcJ3LhpH5YrNPLNkHc+8sZ4bHlvEDY8tYvTAnpxzbH/efXR/xg/uTV6OOtOKdEYBFbHl6xoAdIsvC+UmE0wa0Y9JI/rxzQvGUF3XwDNL1vP0knX8/G/V3PHcm/TIS3LKsL6cPrKYd40q5pgBPTWavUhAARWxZetSPfiO7q+AynbDSwoZXlLIp88cTv2OZmZXb+TvVRt4qWoDz/1lCQDFhXmcNqKYSSP6cVJ5H0aUFCqw5IilgIrY8vUN9O+ZT1GP3KhLkcOoqHsuk8cOZPLYgQCsqd/B36v+FVgzFqwGoE+PXCYO7UNFeV9OKu/DuNIi8nPUO1CODAqoiC1f36DnT8Kgou5cPLGMiyeW4e6s2NjInBWbqFyxicoVm3l6yXoA8nISHF9axHFlRZxQ1pvjyooY1q+AREJXWZJ9FFARcneq1INPOjAzhhUXMKy4gI8EfzY2NOxi7srNzHlrE/NXbeH+V/7JL/++AoDC/BzGlfbi+LLeHFdaxNijejG0XwFJhZZkOAVUhNbU72R7Uysj+usKSvavuDCf88YO5LzglmBLaxtv1m1nQc0WFtbU81ptPb/6+wqagokXu+UmGNW/J8cM7Mnogan3Ywb2pKQwX8+0JGMooCJUtT7Vg29kiQJKDk5OMvF26Oy+ympqaWPZum0sXrOVpWu3sXTtNp5fWsf0uTVv79e3II+jB6Q6awwvLmB4SQHDigsZ3Kc7OZrqRWJGARWhN+uCgNIVlKRBXk6CcaVFjCst2qN9Y8Mulq7dxhtBaC1bv42ZC9ewpbH57W1yEsaQfj0YXlzI8JICyvsVMLhvd8r69OCo3t3UMUMioYCKUNX6Bnp1y6G4UKNhS3j6FeZz2sh8ThtZvEf75u1NVG9ooLpuO9UbtvNW3XaqNzTwwrK6t28VApjBgJ7d3g6swX1S72V9u1PWuwf9e+Vr3EEJhQIqQm/WNTCiv37PRaLRpyCPiQV9mTi07x7trW3Omvod1GzewapNjdRsDj5vbuSVtzbx2PwdtHmHY/XIZUCvbgws6sbAXt06/dynR67+rMtBUUBFqLpuO2cerVl0JV6SCUtdIfXpwanD++21vqmljbX1O1m1uZHazTtYt3Una7fufPv99dqtbNy+C+8QYnk5CYoL8uhXmE+/wjz6FeRTXJj39ud+hXkUB+v6FuTptqIooKKybWcz67ftYoQ6SEiGyctJMKRfD4b02/egt82tbazftou19UFw1e9k3badbGxoYmPDLjZub2L5ugbqGnbR1NLW6TEK83Mo6p5LUfdcevdIvVLLef9qC9YX9cild488enbLoSAvR13ss4QCKiLVddsBGF5SEHElIumXm0xQ2rs7pb2773c7d2d7UysbG3axoV147X6v39FMfWMz9TuaWbaugS2NzdTvaKK51fd73IK8JAX5ORR2y6Fn8F6Q96/l9usK8nMoDN675yXpnhu88pJvL+eqh2MkQg0oM5sC/ARIAne7+/c6rM8Hfg1MBDYCl7j7imDd9cCVQCvwRXd/IsxaD7fqDakefCMUUHIEMzMKg4AY2q9rfxfcnR3NrWxpbA4CKxVaWxqb2bazhW27Wti+q4WGnS007PrXa8O2xj2WWzs+SNuPnIT9K7w6vndoy0smyMtJkJ+TDN5Ty7s/5+fsY30yQX5u+/1TryP5uV1oAWVmSeAO4FygBphjZjPcfXG7za4ENrv7SDO7FLgFuMTMxgCXAmOBo4Cnzexod28Nq97D7a267SQMhvRVQIkcDDOjR14OPfJyOOoAV2j74u7samlj2+4Q29lCY1MLjc2t7GxqZUdz8GoKXs2tNDa1srN9e/Bev6P5X8vNrTS1tNHU0kbLQQTg/iQTRk7CyE0mSCaM3KSRk0iQk0y15SSMnGQiaG//ObHHth3X7d5/9/ETFrwnbM+25J7rcoL1u9v6FOR1+qwyHcK8gjoZqHL3agAzewCYCrQPqKnATcHn6cDtlvpxYSrwgLvvAt4ys6rgeC+HVey3/7SIxau3hnX4vazbupPBfXtoLiCRCJgZ3XKTdMtNUtIzP5TvaG1zmlra2NXSGrynXk0tbTS1trGruTV4Ty3va9uWtjZaWp3mVqelrY3mVqd1d1ub09La9va61HZt7Gxuo6W1Zc/2tjZa2+3zdlub09rme/XM7KoJQ3rz6OdPT+9/vECYAVUKrGq3XAOcsq9t3L3FzOqBfkH7rA77lnb8AjO7CrgKYMiQIWkr/HDITSb2+oVKEckeyd23BfMyozeieyqoWtqcNg/eg/BqbXNa3Wlp7bDOnbwQn89ldCcJd58GTAOoqKg4pOvpG98/Ni01iYhkIgtu58Wpd3+Y95dqgfbDdJcFbZ1uY2Y5QBGpzhJd2VdERLJYmAE1BxhlZsPMLI9Up4cZHbaZAVwRfL4YeNbdPWi/1MzyzWwYMAp4JcRaRUQkZkK7xRc8U7oGeIJUN/N73X2Rmd0MVLr7DOAe4DdBJ4hNpEKMYLuHSHWoaAGuzqYefCIicmDmHccjyVBmVgesPMTDFAMb0lBOXOn8Mlc2nxvo/DLdoZ7fUHffa9y3rAmodDCzSneviLqOsOj8Mlc2nxvo/DJdWOenX8IREZFYUkCJiEgsKaD2NC3qAkKm88tc2XxuoPPLdKGcn55BiYhILOkKSkREYkkBJSIisaSACpjZFDNbamZVZnZd1PWkk5nda2brzez1qGtJNzMbbGbPmdliM1tkZtdGXVM6mVk3M3vFzBYE5/ftqGtKNzNLmtk8M/tz1LWkm5mtMLOFZjbfzCqjrifdzKy3mU03szfMbImZTUrr8fUM6u25q5bRbu4q4LIOc1dlLDM7E2gAfu3u46KuJ53MbBAwyN1fNbOewFzgoiz6f2dAgbs3mFku8BJwrbvPOsCuGcPMvgJUAL3c/YKo60knM1sBVLh7Vv6SrpndB7zo7ncHQ9r1cPct6Tq+rqBS3p67yt2bgN1zV2UFd3+B1FBSWcfd17j7q8HnbcASOpmaJVN5SkOwmBu8suanSjMrA94H3B11LXJwzKwIOJPUkHW4e1M6wwkUULt1NndV1vwjd6Qws3JgAjA74lLSKrgFNh9YDzzl7tl0frcBXwPaIq4jLA48aWZzg/nrsskwoA74ZXCL9m4zS+sU4QooyQpmVgg8DHzJ3Q/f1MiHgbu3uvt4UtPOnGxmWXGb1swuANa7+9yoawnRu9z9ROB84Orgdnu2yAFOBH7m7hOA7UBan98roFI0/1QGC57NPAz8zt0fibqesAS3T54DpkRcSrqcDlwYPKd5ADjbzH4bbUnp5e61wft64FFSjxOyRQ1Q0+6KfjqpwEobBVRKV+aukhgKOhHcAyxx9x9HXU+6mVmJmfUOPncn1ZHnjUiLShN3v97dy9y9nNTfuWfd/WMRl5U2ZlYQdNwhuPU1GcianrTuvhZYZWbHBE3nkJoiKW0yesr3dNnX3FURl5U2ZnY/cBZQbGY1wI3ufk+0VaXN6cDHgYXBcxqAb7j7zOhKSqtBwH1BT9ME8JC7Z1137Cw1AHg09TMUOcDv3f3xaEtKuy8Avwt+sK8GPpXOg6ubuYiIxJJu8YmISCwpoEREJJYUUCIiEksKKBERiSUFlIiIxJICSiQCZtYvGOF6vpmtNbPa4HODmd0ZdX0icaBu5iIRM7ObgAZ3/2HUtYjEia6gRGLEzM7aPS+Smd1kZveZ2YtmttLMPmhm3w/mF3o8GOIJM5toZn8LBiR9IpiCRCTjKaBE4m0EcDZwIfBb4Dl3Pw7YAbwvCKn/Ay5294nAvcB/R1WsSDppqCORePuruzeb2UJSw3DtHipnIVAOHAOMA54KhtRJAmsiqFMk7RRQIvG2C8Dd28ys2f/10LiN1N9fAxa5e1qn2haJA93iE8lsS4ESM5sEqalHzGxsxDWJpIUCSiSDuXsTcDFwi5ktAOYDp0ValEiaqJu5iIjEkq6gREQklhRQIiISSwooERGJJQWUiIjEkgJKRERiSQElIiKxpIASEZFY+v9LkDgH9WG/lwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplot( 2, 1, 1)\n",
    "plt.plot( [sol1.t[0],time_x_on,time_x_on,time_x_off,time_x_off,sol3.t[-1]], [0,0,1,1,0,0] )\n",
    "plt.ylabel( \"X\")\n",
    "plt.subplot( 2, 1, 2 )\n",
    "plt.plot(sol1.t, sol1.y[1,:].T); #Simulation before signal on\n",
    "plt.gca().set_prop_cycle(None) # this resets the color cycle\n",
    "plt.plot(sol2.t, sol2.y[1,:].T); #Simulation after signal on\n",
    "plt.gca().set_prop_cycle(None) # this resets the color cycle\n",
    "plt.plot(sol3.t, sol3.y[1,:].T); #Simulation after signal on\n",
    "plt.xlabel( 'Time')\n",
    "plt.ylabel( \"Z\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the right-hand side of Figure 3 in Mangan & Alon, 2003."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAmt0lEQVR4nO3de5RV9X338feHmeE2XBWwyICDIibERDQTLzGi1RoxpuJTo6IrF2364Eq1NbpWE21rzOVxGZ8njTarVoNoQoy3VqqlWRo1VTGmNTIoXgAxI4nMIMgAwx1mmOH7/DHnnB7GuXk5Z+8z5/Naa699+80+3z0D85mz92//jiICMzOztBmUdAFmZmbdcUCZmVkqOaDMzCyVHFBmZpZKDigzM0ulyqQLeK/GjRsXtbW1SZdhZmbv07JlyzZFxPi+2pVcQNXW1lJfX590GWZm9j5Jeqs/7XyJz8zMUqlgASXpbkkbJb3Ww35J+pGkBkmvSDquULWYmVnpKeQ7qJ8Cs3vZfzZwZGaaB9xewFrMzKzEFOweVEQ8K6m2lyZzgJ9F51hLz0saI2liRKwvVE0fxPz587nvvvuSLsOK5JJLLmHevHlJl2FW1pK8BzUJaMxbb8psexdJ8yTVS6pvbm4uSnFd3XfffSxfvjyR17biWr58uf8YMUuBkujFFxHzgfkAdXV1iY1uO3PmTJ555pmkXt6K5LTTTku6BDMj2XdQ64DJees1mW1mZmaJBtRi4MuZ3nwnAtvSev/JzMyKr2CX+CTdD5wGjJPUBNwAVAFExB3Ao8DngAZgN3BZoWoxM7PSU8hefBf3sT+AKwr1+mZmVto8koSZmaWSA8rMzFLJAWVmZqnkgDIzs1RyQJmZWSo5oMzMLJUcUGZmlkoOKDMzSyUHlJmZpZIDyszMUskBZWZmqeSAMjOzVHJAmZlZKjmgzMwslRxQZmaWSg4oMzNLJQeUmZmlkgPKzMxSyQFlZmap1GNASZrSy75TClOOmZlZp97eQT0j6RuSKrIbJB0i6efALf05uKTZklZLapB0bTf7L5XULGl5ZvqL934KZmY2EPUWUJ8EjgCWSzpd0lXAC8B/A8f3deBMsN0GnA3MAC6WNKObpg9GxMzMtOA9n4GZmQ1IlT3tiIgW4PJMMP0KeBs4MSKa+nns44GGiFgDIOkBYA6w8oOVbGZm5aC3e1BjJP0YuAyYDTwEPCbp9H4eexLQmLfelNnW1fmSXpH0kKTJ/Ty2mZkNcL1d4nsR+B1QFxFPRMTXgS8B/0fS/R/S6/8HUBsRnwCeBBZ210jSPEn1kuqbm5s/pJc2M7M06y2gZkXEDyKiPbshIpZHxKeBp/px7HVA/juimsy2nIjYHBGtmdUFdN73epeImB8RdRFRN378+H68tJmZlboeA6q3e00RcWc/jr0UOFLSVEmDgbnA4vwGkibmrZ4LrOrHcc3MrAz02Enig4qIdklXAo8DFcDdEbFC0neB+ohYDPy1pHOBdmALcGmh6jEzs9JSsIACiIhHgUe7bPtW3vJ1wHWFrMHMzEqThzoyM7NUckCZmVkqOaDMzCyVHFBmZpZKDigzM0slB5SZmaWSA8rMzFLJAWVmZqnkgDIzs1RyQJmZWSo5oMzMLJUcUGZmlkoOKDMzSyUHlJmZpZIDyszMUskBZWZmqeSAMjOzVHJAmZlZKjmgzMwslRxQZmaWSg4oMzNLpYIGlKTZklZLapB0bTf7h0h6MLP/t5JqC1mPmZmVjoIFlKQK4DbgbGAGcLGkGV2afRVoiYhpwC3AzYWqx8zMSktlAY99PNAQEWsAJD0AzAFW5rWZA3w7s/wQ8E+SFBHR00FXr17NaaedVpCCe7N8+XJmzpxZ9Ne1ZCxfvjyRf2dm9j8KGVCTgMa89SbghJ7aRES7pG3AwcCm/EaS5gHzMqutS5Ysea0gFfdhyZIlSEripQHG0eX7UiYSO+8lS5Yk8bLgn3W5KcfzPqo/jQoZUB+aiJgPzAeQVB8RdQmXVHQ+7/JRjucMPu+k6ygmSfX9aVfIThLrgMl56zWZbd22kVQJjAY2F7AmMzMrEYUMqKXAkZKmShoMzAUWd2mzGPhKZvkLwFO93X8yM7PyUbBLfJl7SlcCjwMVwN0RsULSd4H6iFgM3AXcI6kB2EJniPVlfqFqTjmfd/kox3MGn3c56dc5y29YzMwsjTyShJmZpZIDyszMUqmkAqqvoZMGIkl3S9ooKZFnv5IgabKkpyWtlLRC0lVJ11QMkoZKekHSy5nz/k7SNRWLpApJL0n6RdK1FIukP0h6VdLy/na7HggkjZH0kKTXJa2SdFKPbUvlHlRm6KQ3gDPpfOh3KXBxRKzs9QtLnKRZwE7gZxFxdNL1FIOkicDEiHhR0khgGXBeGfysBVRHxE5JVcBzwFUR8XzCpRWcpGuAOmBURHw+6XqKQdIfgLqIKKuHdCUtBH4dEQsyPbyHR8TW7tqW0juo3NBJEdEGZIdOGtAi4lk6eziWjYhYHxEvZpZ3AKvoHHVkQItOOzOrVZmpNP6C/AAk1QDnAAuSrsUKS9JoYBadPbiJiLaewglKK6C6GzppwP/SKneZEe6PBX6bcClFkbnUtRzYCDwZEeVw3rcC3wD2J1xHsQXwhKRlmeHcysFUoBn4SeaS7gJJ1T01LqWAsjIjaQSwCPh6RGxPup5iiIiOiJhJ58grx0sa0Jd1JX0e2BgRy5KuJQGfiYjj6PzEhysyl/MHukrgOOD2iDgW2AX02J+glAKqP0Mn2QCRuQezCLg3Iv4t6XqKLXPZ42lgdsKlFNrJwLmZ+zEPAKdL+nmyJRVHRKzLzDcCD9N5G2OgawKa8q4MPERnYHWrlAKqP0Mn2QCQ6SxwF7AqIn6YdD3FImm8pDGZ5WF0dgh6PdGiCiwirouImoiopfP/9FMR8cWEyyo4SdWZDkBkLnF9FhjwPXUjYgPQKCk7mvkZHPgRTAcoidHMoeehkxIuq+Ak3Q+cBoyT1ATcEBF3JVtVwZ0MfAl4NXM/BuBvI+LR5EoqionAwkyP1UHAv0RE2XS7LjOHAA9nPr6nErgvIn6ZbElF81fAvZk3GmuAy3pqWDLdzM3MrLyU0iU+MzMrIw4oMzNLJQeUmZmlkgPKzMxSyQFlZmap5IAyM7NUckCZmVkqOaDMzCyVHFBmZpZKDigzM0slB5SZmaVSyQwWmzVu3Liora1NugwzM3ufli1btikixvfVruQCqra2lvr6+qTLMDOz90nSW/1pV7BLfJLulrRRUrefcaJOP5LUIOkVST1+aJWZmZWfQt6D+im9fxro2cCRmWkecHsBazEzsxJTsICKiGeBLb00mQP8LDo9D4yRNLFQ9ZiZWWlJshffJKAxb70ps83MzKw0uplLmiepXlJ9c3Nz0uWYmVkRJBlQ64DJees1mW3vEhHzI6IuIurGj++zZ6KZmQ0ASQbUYuDLmd58JwLbImJ9gvWYmVmKFOw5KEn3A6cB4yQ1ATcAVQARcQfwKPA5oAHYDVxWqFrMzKz0FCygIuLiPvYHcEWhXt/MzEpbSXSSMDOz8uOAMjOzVHJAmZlZKjmgzMwslRxQZmaWSg4oMzNLJQeUmZmlkgPKzMxSyQFlZmap5IAyM7NUckCZmVkqOaDMzCyVHFBmZpZKDigzM0slB5SZmaWSA8rMzFLJAWVmZqnkgDIzs1RyQJmZWSr1GFCSpvSy75TClGNmZtapspd9z0i6A/iHiOgAkHQI8A/AR4C6ItRnJWjv3r20tLSwe/duWltbaW1tpa2tjdbWVvbt20dFRcUBU2VlJUOHDmXEiBGMHDmSkSNHUlVVlfRpmFnCeguoTwLfB5ZLugr4OHAN8H+BL/fn4JJmA/8IVAALIuL7XfZfCvw/YF1m0z9FxIL3cgJWHBHBpk2baGhooKmpKTc1NjayYcMGtmzZwpYtW2hpaWHv3r0f+PUGDx7MyJEjGTNmDOPHj2fChAm5+YQJEzjkkEOYPHkyU6ZM4dBDD6Wysrd/ymZWinr8Xx0RLcDlmXD6FfA2cGJENPXnwJIqgNuAM4EmYKmkxRGxskvTByPiyvdVvRXE+vXrefHFF3n55ZdZvXp1btq6dWu/vr6qqoqxY8dSXV3NkCFDGDJkCIMHD2bIkCFUVVWxf/9+2tvb6ejooKOjg/b2dvbu3cuOHTtyU1tbG5s3b2bz5s28+eabvb7eoEGDmDRpElOmTOGwww7jsMMOY/r06Rx55JFMnz6dcePGIelD+M6YWTH1GFCSxgA3AycAs4HPAY9JuioinurHsY8HGiJiTeZ4DwBzgK4BZQnasWMH//Vf/8VvfvMbXnzxRZYtW8aGDRu6bTty5EimT5/OlClTqKmpyU2HHnooBx10EGPHjuWggw5i+PDhHygQIiIXWC0tLTQ3N9Pc3MzGjRvZuHEjzc3NrF+/nsbGRtauXZtbbmxs5De/+c27jjdmzJhcWE2fPp0ZM2bw8Y9/nCOOOMLvvMxSrLf/nS8C/wxcERHtwBOSZgL/LOmtiLi4j2NPAhrz1pvoDLuuzpc0C3gDuDoiGrs2kDQPmAcwZUqPfTesH7Zt28bTTz/Ns88+y7PPPstLL73E/v37D2gzatQojjvuOGbOnMlHP/pRjjrqKI466igOOeSQorwTkcSwYcMYNmwYEyZM4Kijjuq1fWtrK+vWrWPt2rWsXbuWNWvW8Lvf/Y433niDN954g61bt7J06VKWLl16wNcNGTKEj370oxx99NG56eMf/zjjxo2jqqrK98HMEqaI6H6HVNPT5TxJ/zsi7uz1wNIXgNkR8ReZ9S8BJ+RfzpN0MLAzIlolXQ5cFBGn93bcurq6qK+v7/Wk7H9EBCtWrOCxxx7j0Ucf5bnnnqO9vT23v7Kykk9+8pOccsopHH/88Rx33HFMnTqVQYMGxhMIEUFzc3MurFavXs1rr73Ga6+9xtq1a7v9miFDhiCJb37zm8ybN49DDz20yFWbDWySlkVEnx3tegyoD6GAk4BvR8RZmfXrACLiph7aVwBbImJ0b8d1QPUtInjhhRd48MEHWbRo0QG/iAcNGsRJJ53EGWecwSmnnMKJJ57IiBEjEqw2Odu2bWPlypW5wHr11Vd59dVX2bRpU65NZWUl559/PldccQWf+cxnfC/L7EOQhoCqpPOy3Rl09tJbClwSESvy2kyMiPWZ5f8FfDMiTuztuA6o7kUEy5Yt48EHH+Rf//Vfeeutt3L7JkyYwNlnn83nPvc5zjzzTMaOHZtgpen3zjvv8Ktf/YpHHnmEhx9+mI6ODgA+8YlPcMUVV3DJJZeUbaibfRgSD6hMEZ8DbqWzm/ndEXGjpO8C9RGxWNJNwLlAO7AF+FpEvN7bMR1QB9q0aRM///nPueuuu3jttddy2ydNmsQFF1zABRdcwIknnjhgLtkVW1NTEz/+8Y+ZP38+GzduBKC6upoLL7yQyy67zO+qzN6HVARUITigOt8t/epXv2LBggU88sgjtLW1ATBu3Djmzp3LRRddxKc//WmH0oeotbWVRYsWcfvtt/Pcc8/ltk+bNo3LLruML3/5y9TU1CRYoVnpcEANQHv27OGee+7h1ltvZdWqVUDnPaWzzjqLr371q/zpn/4pgwcPTrjKge+NN97gpz/9KQsXLuTtt98GOn8Op59+OnPnzuXP/uzPfBnVrBcOqAFk/fr13Hbbbdxxxx1s3rwZ6LyEd/nll3PppZcyefLkhCssTx0dHTzxxBP85Cc/4ZFHHmHfvn1A54PKs2fPZu7cuZx77rm+X2XWhQNqAGhsbOT73/8+CxYsyF3Gq6ur4+qrr+aCCy7wczopsmXLFh5++GEeeOABnnrqqdyzZcOGDePss8/mvPPO45xzzuGggw5KuFKz5DmgSlhjYyM33XQTd911F21tbUjivPPO45prruHkk0/2TfmUe+edd3jooYe4//77DxjZoqKiglmzZjFnzhzmzJlDbW1tckWaJcgBVYI2bNjA9773Pe6880727duHJC688EKuv/56PvaxjyVdnr0PjY2N/Pu//zuPPPIIS5YsOeAh6WOOOYZzzjmHz372s5x00km+f2hlwwFVQnbt2sUPf/hDbr75Znbt2oUkLrroIq6//npmzJiRdHn2IWlpaeGxxx7jkUce4bHHHmPnzp25fSNGjOD000/nrLPO4qyzzuKII45IsFKzwnJAlYCOjg4WLlzI9ddfn+sNdu6553LTTTc5mAa41tZWnnnmGR5//HEef/xxVq48cAzlww8/nDPPPJNTTz2VU0891cMt2YDigEq5+vp6vva1r5E9l7q6On7wgx9w6qmnJlyZJaGxsZEnn3ySxx9/nCeffJKWlpYD9k+bNo1Zs2blAuuwww5LqFKzD84BlVItLS383d/9HXfccQcRQU1NDTfffDNz5871g7UGdL6zrq+v5+mnn2bJkiU899xzB1wOhM5R/U866SROPPFETjjhBI499liGDh2aUMVm740DKmUignvvvZdrrrmG5uZmKisrufrqq/nWt75VkOdkOjo62L17N3v37qWtrY22tjb27duXm+/fv5+IyE3Z9UGDBuU+ij1/uaKigsGDB+c+eDD/Awjdq7Cw2tvbWb58OUuWLGHJkiX8+te/fteHR1ZVVTFz5sxcYJ1wwgkcccQR/tlYKjmgUuTtt9/m8ssv5xe/+AUAp556Krfddtv77pm3Z88eWlpa2Lp1K9u2bcvNd+3axe7du9m1axf79u1j+PDhDB06lMGDB1NVVZWbV1VVMWjQICQh6YDl/fv309HRkZtnl9vb23NB19raekDoZcMq+xlO+dPw4cPfNa+urmbYsGH+5fk+7d+/nxUrVvD888/z29/+lueff56VK1fS9f/yqFGjmDlzJsceeywzZ85k5syZzJgxw70FLXEOqBSICBYuXMjVV1/N1q1bGT16NLfccguXXnppv34579u3jw0bNvDOO+8c8Kmy7e3tjB07ltGjRzNmzBjGjBnD6NGjqa6uprq6OhdMxQiA/fv350Jrz549uWn37t3drmcDtK2tLVfviBEjepyPGDHCYdYP27dvZ+nSpbnAeuGFF3jnnXfe1a6qqoqPfexjHHvssRxzzDHMmDGDGTNmcOihh/p7bEXjgErYO++8w5//+Z/z6KOPAnDOOefw4x//mEmTJvX4NVu3buWtt96iqamJdevW0dzczPjx4znkkEOYMGEC48ePZ8KECYwcObLkf5m0t7eza9cudu7c2e08u7xz585cmI0YMYKRI0fmgiu7nL+toqIi6VNLjfXr17N8+XJeeuml3LyhoaHbtqNGjcqFVf40efJk3xu1D50DKkG//OUv+cpXvsLGjRsZM2YMP/rRj/jiF7/4rlDZs2cPv//971mzZg2///3v2bt3L7W1tdTU1FBTU8Mf/dEfeTgjOu+nZcNqx44dB8zzl3ft2sXQoUP7DLKRI0eW7fd1x44dvPLKK7z00ku8+uqrrFy5khUrVryr12DW8OHDmTZtGkcccQTTpk07YKqpqXF42fvigEpAa2sr1157LbfeeisAf/zHf8w999xzwLum7du3s3r1al5//XWampqYMmUKhx9+OIcffjgTJkwo+XdGSdq/fz+7d+/uM8h27NhBZWVlj0E2fPhwJDFhwgRGjBgx4H8mEcHGjRtZuXIlq1atYuXKlbmpu8uEWUOGDGHq1Km5ADvssMOYMmVKbj5+/PgB/72z98cBVWS/+93vuPDCC1m+fDmVlZV873vf42/+5m+oqKhgz549vPbaa7zyyits2rSJ6dOn85GPfIQjjjjCN6wTEBHs3bu3xyBraWmhtbWV3bt309rayujRo3NT9n5fdj5q1KgBfVmxpaWFN998k4aGhgOmN998kw0bNvT6tUOHDmXKlCm5KRtckyZN4tBDD2XixImMHTvWIVaGHFBFtHjxYr70pS+xfft2Dj/8cO6//37q6upoaGjg5Zdf5s0332TatGkcc8wxHH744QP6F9pAs2/fvgN6Smbn2eWdO3cyfPhwRo4c2euUfVc2kOzcuTMXXmvWrGHt2rWsXbuWt956i7Vr1/Z42TDfkCFDmDhxYi6wsvP85QkTJnDwwQdTWVlZhLOyYnBAFUFHRwc33HADN954IwDnn38+t912Gw0NDSxdupQRI0Ywc+ZMjj76aD9EOUB1dHSwa9cuduzY0eO0c+dO9u7dm7uMmO1pme1yn13Pnw+Ee2Q7duzIhVZ+cL399tu8/fbbrF+/nu3bt/f7eGPHjmX8+PHvmsaNG/eubWPHjqW6unrA/VEwUDigCmzz5s1ccsklPPHEEwwaNIi///u/p66ujpUrVzJ9+nROOOEEj59mOe3t7blLifnPq/U0HzRo0AGBNXToUIYOHZp7viy73HVeWVlZUr+Ud+3axfr163OB1XW+fv16mpub2bx587ue8+pLZWUlY8aMYezYsYwdO/aA5a7r2eVRo0YN6He9aeGAKqDXX3+dc845hzVr1nDwwQdz5ZVXUl1dzac+9Snq6ur8Car2gUQEbW1tBwTWnj172Lt3b26ev5w/Bw4IreyD2tmHqfNHAeltOW1B19HRwZYtW2hubmbTpk255wKzU/62TZs20dLSwu7duz/Qaw4aNCj3rre/U7aTTfah9O6mwYMHp+p7m4T+BlRBL+pKmg38I1ABLIiI73fZPwT4GfBJYDNwUUT8oZA1fVD/+Z//yfnnn8+2bduYOnUql1xyCX/yJ3/Cpz71KYYMGZJ0eTYASMoFxXv9BN59+/YdEFqtra25KftA9c6dO9myZUtuPX9fdrmjoyM3lFV2yh+JpKqqisrKygPWu2uT37aiooLKysp3Lfenq3pFRUXu8l1/tba2snXrVrZu3UpLS0tu6m09/9Lsnj172L59+3u6DNkfgwYN6jPEhg0blvs3kJ2GDh36rm09be+pbfbnU1FRURIhWbCAklQB3AacCTQBSyUtjoj8zxX4KtASEdMkzQVuBi4qVE0f1J133slf/uVf0t7eztFHH80tt9zCrFmzeuyJt3///gN+UWSHB+o6ZND+/fv7/Y8lIhL7h5X/2tl33u+1lux4f/2ZssMwfdD23W3PDu000GQDYeTIkR/oOB0dHQcMZ9WfadeuXWzdurXbfe3t7bS3t9PR0ZFbzq4DvQZYd8vZeXa8yL7mw4YNo7q6msMOO6zfX5N9bCH76ELXB8rzH1/ITtn2+SOndJ3a2tpyX5uk3v7Q+DD2dffzy079VbBLfJJOAr4dEWdl1q8DiIib8to8nmnz35IqgQ3A+OilKEmJX5M8//zzufHGG2lvb2f79u3s3bu321/cknKXWrJ/vXQdcDX710y5yB+cNjtAbXb5/U49HaOv7b0ZMmQIra2tvZ4HHBjQ+dv6uz8ru63r2IjdTdn2SU75dXc3721f1zbdjfuYv5y/Pxtw+e27m/L3dW3X0778MSi7+7qu/26zP8OufwTl/wy72xYRdHR05IK7ra2N9vb2d4V6fpjnh3x2uad5T8fL7s+eV8ISv8Q3CWjMW28CTuipTUS0S9oGHAxsym8kaR4wL7PaCrxWiIL7a9GiRSxatKjYLzuOLt+XMlGO512O5ww+73JyVH8alcSDBRExH5gPIKm+P8k70Pi8y0c5njP4vJOuo5gk9aunWyEH0loHTM5br8ls67ZN5hLfaDo7S5iZWZkrZEAtBY6UNFXSYGAusLhLm8XAVzLLXwCe6u3+k5mZlY+CXeLL3FO6Eniczm7md0fECknfBeojYjFwF3CPpAZgC50h1pf5hao55Xze5aMczxl83uWkX+dccg/qmplZefCHuZiZWSo5oMzMLJVKKqAkzZa0WlKDpGuTrqcYJN0taaOkRJ/9KiZJkyU9LWmlpBWSrkq6pmKQNFTSC5Jezpz3d5KuqVgkVUh6SdIvkq6lWCT9QdKrkpb3t9v1QCBpjKSHJL0uaVVmUIfu25bKPajM0ElvkDd0EnBxl6GTBhxJs4CdwM8i4uik6ykGSROBiRHxoqSRwDLgvDL4WQuojoidkqqA54CrIuL5hEsrOEnXAHXAqIj4fNL1FIOkPwB1EVFWD+lKWgj8OiIWZHp4D4+Ird21LaV3UMcDDRGxJiLagAeAOQnXVHAR8SydPRzLRkSsj4gXM8s7gFV0jjoyoEWn7ABtVZmpNP6C/AAk1QDnAAuSrsUKS9JoYBadPbiJiLaewglKK6C6GzppwP/SKneSaoFjgd8mXEpRZC51LQc2Ak9GRDmc963AN4DEB4grsgCekLQsM5xbOZgKNAM/yVzSXSCpuqfGpRRQVmYkjQAWAV+PiA/3Mw9SKiI6ImImnSOvHC9pQF/WlfR5YGNELEu6lgR8JiKOA84Grshczh/oKoHjgNsj4lhgF9Bjf4JSCqj+DJ1kA0TmHswi4N6I+Lek6ym2zGWPp4HZCZdSaCcD52buxzwAnC7p58mWVBwRsS4z3wg8TOdtjIGuCWjKuzLwEJ2B1a1SCqj+DJ1kA0Cms8BdwKqI+GHS9RSLpPGSxmSWh9HZIej1RIsqsIi4LiJqIqKWzv/TT0XEFxMuq+AkVWc6AJG5xPVZEv6UhmKIiA1Ao6TsaOZnAD12fiqJ0cyh56GTEi6r4CTdD5wGjJPUBNwQEXclW1XBnQx8CXg1cz8G4G8j4tHkSiqKicDCTI/VQcC/RETZdLsuM4cAD2c+G6sSuC8ifplsSUXzV8C9mTcaa4DLempYMt3MzcysvJTSJT4zMysjDigzM0slB5SZmaWSA8rMzFLJAWVmZqnkgDIzs1RyQJmZWSo5oMzMLJUcUGZmlkoOKDMzSyUHlJmZpVLJDBabNW7cuKitrU26DDMze5+WLVu2KSLG99Wu5AKqtraW+vr6pMswM7P3SdJb/WnnS3xmZpZKBQsoSXdL2iip2w/hUqcfSWqQ9IqkHj9V0czMyk8h30H9lN4/rvps4MjMNA+4vYC1mJlZiSlYQEXEs8CWXprMAX4WnZ4HxkiaWKh6zMystCR5D2oS0Ji33pTZ9i6S5kmql1Tf3NxclOLMzCxZJdFJIiLmR0RdRNSNH99nz0QzMxsAkgyodcDkvPWazDYzM7NEA2ox8OVMb74TgW0RsT7BeszMLEUK9qCupPuB04BxkpqAG4AqgIi4A3gU+BzQAOwGLitULWZmVnoKFlARcXEf+wO4olCvb2Zmpa0kOkmYmVn5cUCZmVkqOaDMzCyVHFBmZpZKDigzM0slB5SZmaWSA8rMzFLJAWVmZqnkgDIzs1RyQJmZWSo5oMzMLJUcUGZmlkoOKDMzSyUHlJmZpZIDyszMUskBZWZmqeSAMjOzVHJAmZlZKjmgzMwslXoNKEln9bLvgg+/HDMzs059vYN6VNLTkiZ1s++6vg4uabak1ZIaJF3bzf5LJTVLWp6Z/qK/hZuZ2cDWV0C9AtwHPC/pC132qbcvlFQB3AacDcwALpY0o5umD0bEzMy0oJ91m5nZANdXQEVE3AmcAXxT0k8kDc/u6+NrjwcaImJNRLQBDwBzPli5ZmZWLvrVSSIi3gBOAt4BXpJ0Qj++bBLQmLfelNnW1fmSXpH0kKTJ3R1I0jxJ9ZLqm5ub+1OymZmVuL4CKncZLyLaI+Ja4HLgfuDID+H1/wOojYhPAE8CC7trFBHzI6IuIurGjx//IbysmZmlXV8B9Z2uGyLiGeCTwI19fO06IP8dUU1mW/6xNkdEa2Z1Qea4ZmZmvQdURDzSw/aWiPh+H8deChwpaaqkwcBcYHF+A0kT81bPBVb1WbGZmZWFykIdOCLaJV0JPA5UAHdHxApJ3wXqI2Ix8NeSzgXagS3ApYWqx8zMSosi+uqMly51dXVRX1+fdBlmZvY+SVoWEXV9tfNQR2ZmlkoOKDMzSyUHlJmZpZIDyszMUskBZWZmqeSAMjOzVHJAmZlZKjmgzMwslRxQZmaWSg4oMzNLJQeUmZmlkgPKzMxSyQFlZmap5IAyM7NUckCZmVkqOaDMzCyVHFBmZpZKDigzM0slB5SZmaWSA8rMzFLJAWVmZqlU0ICSNFvSakkNkq7tZv8QSQ9m9v9WUm0h6zEzs9JRsICSVAHcBpwNzAAuljSjS7OvAi0RMQ24Bbi5UPWYmVlpqSzgsY8HGiJiDYCkB4A5wMq8NnOAb2eWHwL+SZIiIno66LJly5BUmIrNzCw1ChlQk4DGvPUm4ISe2kREu6RtwMHApvxGkuYB8zKrrcBrhSg45cbR5ftSJsrxvMvxnMHnXU6O6k+jQgbUhyYi5gPzASTVR0RdwiUVnc+7fJTjOYPPO+k6iklSfX/aFbKTxDpgct56TWZbt20kVQKjgc0FrMnMzEpEIQNqKXCkpKmSBgNzgcVd2iwGvpJZ/gLwVG/3n8zMrHwU7BJf5p7SlcDjQAVwd0SskPRdoD4iFgN3AfdIagC20BlifZlfqJpTzuddPsrxnMHnXU76dc7yGxYzM0sjjyRhZmap5IAyM7NUKqmA6mvopIFI0t2SNkoqm2e/JE2W9LSklZJWSLoq6ZqKQdJQSS9Iejlz3t9JuqZikVQh6SVJv0i6lmKR9AdJr0pa3t9u1wOBpDGSHpL0uqRVkk7qsW2p3IPKDJ30BnAmnQ/9LgUujoiVvX5hiZM0C9gJ/Cwijk66nmKQNBGYGBEvShoJLAPOK4OftYDqiNgpqQp4DrgqIp5PuLSCk3QNUAeMiojPJ11PMUj6A1AXEWX1kK6khcCvI2JBpof38IjY2l3bUnoHlRs6KSLagOzQSQNaRDxLZw/HshER6yPixczyDmAVnaOODGjRaWdmtSozlcZfkB+ApBrgHGBB0rVYYUkaDcyiswc3EdHWUzhBaQVUd0MnDfhfWuUuM8L9scBvEy6lKDKXupYDG4EnI6IczvtW4BvA/oTrKLYAnpC0LDOcWzmYCjQDP8lc0l0gqbqnxqUUUFZmJI0AFgFfj4jtSddTDBHREREz6Rx55XhJA/qyrqTPAxsjYlnStSTgMxFxHJ2f+HBF5nL+QFcJHAfcHhHHAruAHvsTlFJA9WfoJBsgMvdgFgH3RsS/JV1PsWUuezwNzE64lEI7GTg3cz/mAeB0ST9PtqTiiIh1mflG4GE6b2MMdE1AU96VgYfoDKxulVJA9WfoJBsAMp0F7gJWRcQPk66nWCSNlzQmszyMzg5BrydaVIFFxHURURMRtXT+n34qIr6YcFkFJ6k60wGIzCWuz1IGn9IQERuARknZ0czP4MCPYDpASYxmDj0PnZRwWQUn6X7gNGCcpCbghoi4K9mqCu5k4EvAq5n7MQB/GxGPJldSUUwEFmZ6rA4C/iUiyqbbdZk5BHg489l2lcB9EfHLZEsqmr8C7s280VgDXNZTw5LpZm5mZuWllC7xmZlZGXFAmZlZKjmgzMwslRxQZmaWSg4oMzNLJQeUmZmlkgPKzMxS6f8DmN7/MWkTpWsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "B_y = 0\n",
    "Beta_y = 1\n",
    "B_z = 0\n",
    "Beta_z = 1\n",
    "K_xy = 0.1\n",
    "K_xz = 0.1\n",
    "H = 2\n",
    "alpha_y = 1\n",
    "alpha_z = 1\n",
    "\n",
    "time_x_on = 0.5\n",
    "time_x_off = 2.5\n",
    "\n",
    "colors = ['gray', 'gray', 'black']\n",
    "thickness = [0.5,1,2]\n",
    "K_yzs = [1, 0.3, .1]\n",
    "fig, axs = plt.subplots( nrows=3, ncols=1 )\n",
    "\n",
    "# X on top\n",
    "axs[0].plot( [sol1.t[0],time_x_on,time_x_on,time_x_off,time_x_off,sol3.t[-1]], \n",
    "            [0,0,1,1,0,0], color='k' )\n",
    "axs[0].set_xlim( [0,6])\n",
    "axs[0].set_ylabel( \"X\")\n",
    "axs[0].set_ylim( [0,1.1])\n",
    "\n",
    "\n",
    "# Z when SY=1\n",
    "for i in range(3):\n",
    "    K_yz = K_yzs[i]\n",
    "    params_off = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,0,1);\n",
    "    params_on = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,1,1);\n",
    "\n",
    "    mf_off = lambda t,y : incoherent_type_4_and_model(t,y, params_off)\n",
    "    mf_on = lambda t,y : incoherent_type_4_and_model(t,y, params_on)\n",
    "\n",
    "    sol1 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(0,time_x_on), y0=(0,0), \n",
    "                                    t_eval=np.arange(0,time_x_on,.01) )\n",
    "    sol2 = scipy.integrate.solve_ivp( fun=mf_on, t_span=(time_x_on,time_x_off), y0=sol1.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_on,time_x_off,.01) )\n",
    "    sol3 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(time_x_off,6), y0=sol2.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_off,6,.01) )\n",
    "\n",
    "    ax = axs[1]\n",
    "    ax.plot(sol1.t, sol1.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation before signal on\n",
    "    ax.plot(sol2.t, sol2.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "    ax.plot(sol3.t, sol3.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "axs[1].set_ylabel( \"Z\")\n",
    "axs[1].set_xlim( [0,6])\n",
    "axs[1].set_ylim( [0,1.1])\n",
    "\n",
    "\n",
    "# Z when SY=0\n",
    "for i in range(3):\n",
    "    K_yz = K_yzs[i]\n",
    "    params_off = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,0,0);\n",
    "    params_on = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,1,0);\n",
    "\n",
    "    mf_off = lambda t,y : incoherent_type_4_and_model(t,y, params_off)\n",
    "    mf_on = lambda t,y : incoherent_type_4_and_model(t,y, params_on)\n",
    "\n",
    "    sol1 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(0,time_x_on), y0=(0,0), \n",
    "                                    t_eval=np.arange(0,time_x_on,.01) )\n",
    "    sol2 = scipy.integrate.solve_ivp( fun=mf_on, t_span=(time_x_on,time_x_off), y0=sol1.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_on,time_x_off,.01) )\n",
    "    sol3 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(time_x_off,6), y0=sol2.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_off,6,.01) )\n",
    "\n",
    "    ax = axs[2]\n",
    "    ax.plot(sol1.t, sol1.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation before signal on\n",
    "    ax.plot(sol2.t, sol2.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "    ax.plot(sol3.t, sol3.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "axs[2].set_ylabel( \"Z\")\n",
    "axs[2].set_xlim( [0,6])\n",
    "axs[2].set_ylim( [0,1.1])\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "                                 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Incoherent Type 1 FFL with AND gate\n",
    "\n",
    "Implement the incoherent type 1 FFL with an AND gate, as described by Mangan & Alon (PNAS, 2003).\n",
    "\n",
    "The parameters are\n",
    "- B_y: basal rate of Y's transcription\n",
    "- Beta_y: maximum rate of Y's transcription\n",
    "- B_z: basal rate of Z's transcription\n",
    "- Beta_z: maximum rate of Z's transcription\n",
    "- K_xy: the activation threshold for X's activation of Y\n",
    "- K_xz: the activation threshold for X's activation of Z\n",
    "- K_yz: the activation threshold for Y's activation of Z\n",
    "- H: the Hill coefficient (exponent) for all regulation\n",
    "- alpha_y: the degradation rate constant for Y\n",
    "- alpha_z: the degradation rate constant for Z\n",
    "- SX: amount of active X in system (X*=SX)\n",
    "- SY: whether or not Y is activated (Y*=Y if SY=1, Y*=0 if SY=0)\n",
    "\n",
    "The model is\n",
    "\\begin{eqnarray*}\n",
    "\\frac{dY}{dt} &=& B_y + \\beta_Y \\frac{(X^*/K_{xy})^H}{1+(X^*/K_{xy})^H} - \\alpha_y Y\\\\\n",
    "\\frac{dZ}{dt} &=& B_z + \\beta_Z \\left(\\frac{(X^*/K_{xz})^H}{1+(X^*/K_{xz})^H}\\right) \\left(\\frac{1}{1+(Y^*/K_{yz})^H}\\right) - \\alpha_z Z\\\\\n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def incoherent_type_1_and_model( t, statevars, params ):\n",
    "    Y = statevars[0];\n",
    "    Z = statevars[1];\n",
    "    B_y = params[0]\n",
    "    Beta_y = params[1]\n",
    "    B_z = params[2]\n",
    "    Beta_z = params[3]\n",
    "    K_xy = params[4]\n",
    "    K_xz = params[5]\n",
    "    K_yz = params[6]\n",
    "    H = params[7]\n",
    "    alpha_y = params[8]\n",
    "    alpha_z = params[9]\n",
    "    SX = params[10]\n",
    "    SY = params[11]\n",
    "    \n",
    "    if SX == 0:\n",
    "        Xstar = 0\n",
    "    else:\n",
    "        Xstar = 1\n",
    "    if SY == 0:\n",
    "        Ystar = 0\n",
    "    else:\n",
    "        Ystar = Y\n",
    "    \n",
    "    dydt = np.zeros( (2,) )\n",
    "    yreg = ((Xstar/K_xy)**H)/(1 + (Xstar/K_xy)**H)\n",
    "    zregx = ((Xstar/K_xz)**H)/(1 + (Xstar/K_xz)**H)\n",
    "    zregy = (1)/(1 + (Ystar/K_yz)**H)\n",
    "    zreg = zregx * zregy\n",
    "    dydt[0] = B_y + Beta_y*yreg - alpha_y*Y;\n",
    "    dydt[1] = B_z + Beta_z*zreg - alpha_z*Z;\n",
    "\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate the model\n",
    "\n",
    "Simulate the model. In the paper, they specify that for most simulations,\n",
    "\n",
    "- B_y = 0\n",
    "- Beta_y = 1\n",
    "- B_z = 0\n",
    "- Beta_z = 1\n",
    "- H = 2\n",
    "- alpha_y = 1\n",
    "- alpha_z = 1\n",
    "\n",
    "The activation thresholds are given in the figure captions, as are SX and SY."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "B_y = 0\n",
    "Beta_y = 1\n",
    "B_z = 0\n",
    "Beta_z = 1\n",
    "K_xy = 0.1\n",
    "K_xz = 0.1\n",
    "K_yz = 0.3 # 0.01, 0.3, 0.1\n",
    "H = 2\n",
    "alpha_y = 1\n",
    "alpha_z = 1\n",
    "\n",
    "time_x_on = 0.5\n",
    "time_x_off = 2.5\n",
    "params_off = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,0,1);\n",
    "params_on = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,1,1);\n",
    "\n",
    "mf_off = lambda t,y : incoherent_type_1_and_model(t,y, params_off)\n",
    "mf_on = lambda t,y : incoherent_type_1_and_model(t,y, params_on)\n",
    "\n",
    "sol1 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(0,time_x_on), y0=(0,0), \n",
    "                                t_eval=np.arange(0,time_x_on,.01) )\n",
    "sol2 = scipy.integrate.solve_ivp( fun=mf_on, t_span=(time_x_on,time_x_off), y0=sol1.y[:,-1], \n",
    "                                t_eval=np.arange(time_x_on,time_x_off,.01) )\n",
    "sol3 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(time_x_off,6), y0=sol2.y[:,-1], \n",
    "                                t_eval=np.arange(time_x_off,6,.01) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the simulation\n",
    "Plot X over time in an upper subplot and Z over time in a lower subplot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAEYCAYAAADrpHnMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjJ0lEQVR4nO3dfZRddX3v8fdnzjxl8jAzSQYIeSABAoKCBGYFlV6kKhirEnvVRXDZorUrtgUf6rr1QntXQWy70Hqrt1dqSSEtViVYRBu9KFCBqqsgmUAEMiEQApIZAplk8jSZZB7OfO8fZyecTGbyYM7OmX3m81prr9n799v7nO8mJJ/Z+/zObysiMDMzy5KqchdgZmZ2rBxeZmaWOQ4vMzPLHIeXmZlljsPLzMwyp7rcBZTK9OnTY+7cueUuw8zMSmT16tVbI6JlpL6KCa+5c+fS1tZW7jLMzKxEJP16tD7fNjQzs8xxeJmZWeY4vMzMLHMcXmZmljkOLzMzy5zUwkvScklbJD0zSr8k/b2kDZKeknRhUd81kp5PlmvSqtHMzLIpzSuvfwEWHab/PcD8ZFkKfANA0lTgRuBiYCFwo6TmFOs0M7OMSe17XhHxM0lzD7PLYuCbUXgmy2OSmiTNAC4DHoyIbgBJD1IIwbvSqrVcvvPLl/n3NZ3lLsOO0+ILZvKRi+eUuwyzcaWcn3nNBDYVbXckbaO1H0LSUkltktq6urpSKzQt/76mk/bNu8pdhh2H9s27/AuIWRlkeoaNiFgGLANobW3N5FM1z50xhbs/+dZyl2G/oatue7TcJZiNS+W88uoEZhdtz0raRms3MzMDyhteK4HfT0YdvgXYGRGbgfuBKyQ1JwM1rkjazMzMgBRvG0q6i8Lgi+mSOiiMIKwBiIh/BO4DfgfYAPQCH0/6uiV9EViVvNTN+wdvmJmZQbqjDa8+Qn8A147StxxYnkZdZmaWfZ5hw8zMMsfhZWZmmePwMjOzzHF4mZlZ5ji8zMwscxxeZmaWOQ4vMzPLHIeXmZlljsPLzMwyx+FlZmaZ4/AyM7PMcXiZmVnmOLzMzCxzHF5mZpY5qYaXpEWS1kvaIOn6Efq/KmlNsjwnaUdRX76ob2WadZqZWbak+TDKHHArcDnQAayStDIi2vfvExF/WrT/p4AFRS+xNyIuSKs+MzPLrjSvvBYCGyJiY0T0AyuAxYfZ/2rgrhTrMTOzCpFmeM0ENhVtdyRth5B0GjAPeKiouV5Sm6THJH1glOOWJvu0dXV1lahsMzMb68bKgI0lwD0RkS9qOy0iWoGPAF+TdMbwgyJiWUS0RkRrS0vLiarVzMzKLM3w6gRmF23PStpGsoRhtwwjojP5uRF4hIM/DzMzs3EszfBaBcyXNE9SLYWAOmTUoKQ3AM3Ao0VtzZLqkvXpwCVA+/BjzcxsfEpttGFEDEq6DrgfyAHLI2KtpJuBtojYH2RLgBUREUWHnwPcJmmIQsDeUjxK0czMxrfUwgsgIu4D7hvW9pfDtm8a4bj/As5LszYzM8uusTJgw8zM7Kg5vMzMLHMcXmZmljkOLzMzyxyHl5mZZY7Dy8zMMsfhZWZmmePwMjOzzHF4mZlZ5ji8zMwscxxeZmaWOQ4vMzPLnFHDS9Kcw/T9t3TKMTMzO7LDXXk9IunzknL7GySdLOlbwFfTL83MzGxkhwuvi4AzgDWS3iHpM8DjFB4aufBEFGdmZjaSUcMrIrZHxCeB24H/AP4MuCQibo2IoaN5cUmLJK2XtEHS9SP0f0xSl6Q1yfKHRX3XSHo+Wa459lMzM7NKNerDKCU1AV8CLgYWAb8D/FjSZyLioSO9cHK78VbgcqADWCVp5QhPRL47Iq4bduxU4EagFQhgdXLs9qM+MzMzq1iHu234BPA80BoRD0TEZ4HfA/5K0l1H8doLgQ0RsTEi+oEVwOKjrOvdwIMR0Z0E1oMUAtTMzOyw4XVpRHwlIgb3N0TEmoh4G3DEKy9gJrCpaLsjaRvug5KeknSPpNnHcqykpZLaJLV1dXUdRUlmZlYJDveZV8dh+v6pRO//Q2BuRJxP4erqzmM5OCKWRURrRLS2tLSUqCQzMxvr0vyScicwu2h7VtJ2QERsi4i+ZPN2CiMcj+pYMzMbv9IMr1XAfEnzJNUCS4CVxTtImlG0eSWwLlm/H7hCUrOkZuCKpM3MzGz00YbHKyIGJV1HIXRywPKIWCvpZqAtIlYCn5Z0JTAIdAMfS47tlvRFCgEIcHNEdKdVq5mZZUtq4QUQEfcB9w1r+8ui9RuAG0Y5djmwPM36zMwsmzwxr5mZZY7Dy8zMMsfhZWZmmePwMjOzzHF4mZlZ5ji8zMwscxxeZmaWOQ4vMzPLHIeXmZlljsPLzMwyx+FlZmaZ4/AyM7PMcXiZmVnmOLzMzCxzUg0vSYskrZe0QdL1I/R/TlK7pKck/VTSaUV9eUlrkmXl8GPNzGz8Su15XpJywK3A5UAHsErSyohoL9rtSaA1Inol/THwZeCqpG9vRFyQVn1mZpZdaV55LQQ2RMTGiOgHVgCLi3eIiIcjojfZfAyYlWI9ZmZWIdIMr5nApqLtjqRtNJ8Afly0XS+pTdJjkj6QQn1mZpZRqd02PBaSPgq0Am8vaj4tIjolnQ48JOnpiHhh2HFLgaUAc+bMOWH1mplZeaV55dUJzC7anpW0HUTSu4C/AK6MiL797RHRmfzcCDwCLBh+bEQsi4jWiGhtaWkpbfVmZjZmpRleq4D5kuZJqgWWAAeNGpS0ALiNQnBtKWpvllSXrE8HLgGKB3qYmdk4ltptw4gYlHQdcD+QA5ZHxFpJNwNtEbES+FtgEvBvkgBejogrgXOA2yQNUQjYW4aNUjQzs3Es1c+8IuI+4L5hbX9ZtP6uUY77L+C8NGszM7Ps8gwbZmaWOQ4vMzPLHIeXmZlljsPLzMwyx+FlZmaZ4/AyM7PMcXiZmVnmOLzMzCxzHF5mZpY5Di8zM8sch5eZmWWOw8vMzDLH4WVmZpnj8DIzs8xxeJmZWeakGl6SFklaL2mDpOtH6K+TdHfS/0tJc4v6bkja10t6d5p1mplZtqQWXpJywK3Ae4BzgaslnTtst08A2yPiTOCrwJeSY88FlgBvBBYB/5C8npmZWapPUl4IbIiIjQCSVgCLgfaifRYDNyXr9wBfl6SkfUVE9AEvStqQvN6jaRX7hR+upf2VXWm9/IjaN+/i3BlTTuh7Wum1b97FVbel9r+mWSade+oUbnz/G1N7/TTDayawqWi7A7h4tH0iYlDSTmBa0v7YsGNnDn8DSUuBpclmj6T1x1nzdGDrcb7GMXkG+O4fnch3PPHnWAZl+XM8gfxnWBkq/hxvOv5zPG20jjTDK3URsQxYVqrXk9QWEa2ler2xyOeYfZV+fuBzrBRpnmOaAzY6gdlF27OSthH3kVQNNALbjvJYMzMbp9IMr1XAfEnzJNVSGICxctg+K4FrkvUPAQ9FRCTtS5LRiPOA+cDjKdZqZmYZktptw+QzrOuA+4EcsDwi1kq6GWiLiJXAHcC/JgMyuikEHMl+36UwuGMQuDYi8mnVWqRktyDHMJ9j9lX6+YHPsVKkdo4qXOiYmZllh2fYMDOzzHF4mZlZ5ji8EkeayirrJC2XtEXSCf5K0okhabakhyW1S1or6TPlrqnUJNVLelzSr5Jz/EK5a0qDpJykJyX9qNy1pEHSS5KelrRGUlu560mDpCZJ90h6VtI6SW8t+Xv4M68DU1k9B1xO4QvRq4CrI6L9sAdmiKRLgR7gmxHxpnLXU2qSZgAzIuIJSZOB1cAHKuzPUMDEiOiRVAP8AvhMRDx2hEMzRdLngFZgSkS8r9z1lJqkl4DWiKjYLyhLuhP4eUTcnow2b4iIHaV8D195FRyYyioi+oH9U1lVjIj4GYURnRUpIjZHxBPJ+m5gHSPMypJlUdCTbNYkS0X99ilpFvBe4PZy12K/GUmNwKUURpMTEf2lDi5weO030lRWFfUP33iSPJ1gAfDLMpdScskttTXAFuDBiKi0c/wa8HlgqMx1pCmAByStTqa4qzTzgC7gn5Pbv7dLmljqN3F4WUWRNAn4HvDZiDixMy2fABGRj4gLKMw6s1BSxdwClvQ+YEtErC53LSn7rYi4kMITN65NbulXkmrgQuAbEbEA2AOUfByBw6vA01FVgORzoO8B346Ie8tdT5qS2zAPU3hkUKW4BLgy+UxoBfAOSd8qb0mlFxGdyc8twPcpfGxRSTqAjqK7AvdQCLOScngVHM1UVjaGJYMZ7gDWRcTflbueNEhqkdSUrE+gMMDo2bIWVUIRcUNEzIqIuRT+Dj4UER8tc1klJWliMqCI5FbaFZzwhxKkKyJeBTZJOjtpeicHPwqrJDI9q3ypjDaVVZnLKilJdwGXAdMldQA3RsQd5a2qpC4Bfg94OvlMCODPI+K+8pVUcjOAO5PRsVXAdyOiIoeTV7CTge8XfteiGvhORPykvCWl4lPAt5OLgY3Ax0v9Bh4qb2ZmmePbhmZmljkOLzMzyxyHl5mZZY7Dy8zMMsfhZWZmmePwMjOzzHF4mZlZ5ji8zMwscxxeZmaWOQ4vMzPLHIeXmZllTsVMzDt9+vSYO3duucswM7MSWb169daIaBmpr2LCa+7cubS1tZW7DDMzKxFJvx6tz7cNzcwscxxeZmaWOQ4vMzPLHIeXmZllTmrhJWm5pC2SnhmlX5L+XtIGSU9JurCo7xpJzyfLNWnVaGZm2ZTmlde/AIsO0/8eYH6yLAW+ASBpKnAjcDGwELhRUnOKdZqZWcakFl4R8TOg+zC7LAa+GQWPAU2SZgDvBh6MiO6I2A48yOFD0MzMxplyfuY1E9hUtN2RtI3WfghJSyW1SWrr6upKrVAzMxtbMj1gIyKWRURrRLS2tIz4JWwzM6tA5QyvTmB20faspG20djMzM6C84bUS+P1k1OFbgJ0RsRm4H7hCUnMyUOOKpM3MzAxIcW5DSXcBlwHTJXVQGEFYAxAR/wjcB/wOsAHoBT6e9HVL+iKwKnmpmyPicAM/zMxsnEktvCLi6iP0B3DtKH3LgeVp1GVmZtmX6QEbZmY2Pjm8zMwscxxeZmaWOQ4vMzPLHIeXmZlljsPLzMwyx+FlZmaZ4/AyM7PMcXiZmVnmOLzMzCxzHF5mZpY5Di8zM8sch5eZmWWOw8vMzDLH4WVmZpmTanhJWiRpvaQNkq4fof+rktYky3OSdhT15Yv6VqZZp5mZZUuaT1LOAbcClwMdwCpJKyOiff8+EfGnRft/ClhQ9BJ7I+KCtOozM7PsSvPKayGwISI2RkQ/sAJYfJj9rwbuSrEeMzOrEGmG10xgU9F2R9J2CEmnAfOAh4qa6yW1SXpM0gdGOW5psk9bV1dXico2M7OxbqwM2FgC3BMR+aK20yKiFfgI8DVJZww/KCKWRURrRLS2tLScqFrNzKzM0gyvTmB20faspG0kSxh2yzAiOpOfG4FHOPjzMDMzG8fSDK9VwHxJ8yTVUgioQ0YNSnoD0Aw8WtTWLKkuWZ8OXAK0Dz/WzMzGp9RGG0bEoKTrgPuBHLA8ItZKuhloi4j9QbYEWBERUXT4OcBtkoYoBOwtxaMUzcxsfNPBmZFdra2t0dbWVu4yzMysRCStTsY+HGKsDNgwMzM7ag4vMzPLHIeXmZlljsPLzMwyx+FlZmaZ4/AyM7PMcXiZmVnmOLzMzCxzHF5mZpY5Di8zM8sch5eZmWXOYcNL0rsP0/fh0pdjZmZ2ZEe68rpP0sOSRnoC8g1pFGRmZnYkRwqvp4DvAI9J+tCwPqVTkpmZ2eEdKbwiIv4JeCfwPyX9s6SG/X1HenFJiyStl7RB0vUj9H9MUpekNcnyh0V910h6PlmuOZaTMjOzynZUAzYi4jngrcBrwJOSLj7SMZJywK3Ae4BzgaslnTvCrndHxAXJcnty7FTgRuBiYCFwo6Tmo6nVzMwq35HC68CtwYgYjIjrgU8CdwHzj3DsQmBDRGyMiH5gBbD4KOt6N/BgRHRHxHbgQWDRUR5rZmYV7kjh9YXhDRHxCHAR8NdHOHYmsKlouyNpG+6Dkp6SdI+k2cd4rJmZjUOHDa+I+MEo7dsj4pYSvP8PgbkRcT6Fq6s7j+VgSUsltUlq6+rqKkE5ZmaWBWl+SbkTmF20PStpOyAitkVEX7J5O4UruqM6Njl+WUS0RkRrS0tLyQo3M7OxLc3wWgXMlzRPUi2wBFhZvIOkGUWbVwLrkvX7gSskNScDNa5I2szMzKhO64UjYlDSdRRCJwcsj4i1km4G2iJiJfBpSVcCg0A38LHk2G5JX6QQgAA3R0R3WrWamVm2KOKIX9fKhNbW1mhrayt3GWZmViKSVkdE60h9npjXzMwyx+FlZmaZ4/AyM7PMcXiZmVnmOLzMzCxzHF5mZpY5Di8zM8sch5eZmWWOw8vMzDLH4WVmZpnj8DIzs8xxeJmZWeY4vMzMLHMcXmZmljkOLzMzy5xUw0vSIknrJW2QdP0I/Z+T1C7pKUk/lXRaUV9e0ppkWTn8WDMzG79Se5KypBxwK3A50AGskrQyItqLdnsSaI2IXkl/DHwZuCrp2xsRF6RV33gWEezoHaBzx14279zHqzv3snPvAD19eXr6Bujty5OPQMn+kshViYm1OSbWVReW2hyT6muYNrGWaZNqmTapjmkTa6mvyZX13MxsfEgtvICFwIaI2AggaQWwGDgQXhHxcNH+jwEfTbGecatrdx+PbtzGU5t28Oyru1m3eRfb9vQfsl9tddWBgMpVFaIrAoJgMB/09ufZ0zfI4NDoT9+eVFfNtEm1TJ9Ux4zGek5tmsCMxnpmNE7g1KbCz+mTapE06muYmR1JmuE1E9hUtN0BXHyY/T8B/Lhou15SGzAI3BIRPxh+gKSlwFKAOXPmHG+9FeW513bzgyc7+em6Lax/bTcAddVVnHXyZN55zkmcdfJkZjVPYEZjIVyaGmqprT66u8h9g3l6+/Ls3jfItj19bOvpZ9uePrb29B9Y37Krj2c6d/JA+2v0Dw4ddHxtdRUzGuuZ1TyBOVMbmD21gTnJMru5gaaGGoebmR1WmuF11CR9FGgF3l7UfFpEdEo6HXhI0tMR8ULxcRGxDFgG0NraOvrlwDixtz/P957o4FuP/ZpnX91NleAtp0/jzy44m0vOnM6bTp1Cde74P+asq85RV52jeWItc6Y1HHbfiKB7Tz+bd+4r3KZMblV27thLx/a9PLD2tUOuAifXVR8ItNlTDw64mc0TqKv2rUmz8S7N8OoEZhdtz0raDiLpXcBfAG+PiL797RHRmfzcKOkRYAHwwvDjDfYN5LnjFy9yxy9epHtPP+fNbOSm95/Le88/lZbJdWWtTVLh87BJdbxpZuOI++zpG2TT9l5e3tbLy929dGzfy8vdvWzo6uHh9VvoK7pyk+CUKfXMTq7S5kxtYM60CQfWWybX+arNbBxIM7xWAfMlzaMQWkuAjxTvIGkBcBuwKCK2FLU3A70R0SdpOnAJhcEcViQiWPmrV/jSj5/llZ37+O2zW/ijt5/BwnlTM/UP+MS6at5wyhTecMqUQ/qGhoKtPX283N17YNnUvZdN3b381wtbuffJfUTRNXddddXrV23NE4qu4ArLpLoxcbPBzI5Tan+TI2JQ0nXA/UAOWB4RayXdDLRFxErgb4FJwL8l/9i+HBFXAucAt0kaojCc/5ZhoxTHvS279nHDvU/z02e3cN7MRv7uqgt4y+nTyl1WyVVViZOm1HPSlHpa5049pH/fQJ7OHYUw21QUbi9397LqxW529w0etP/UibUHhVvxLckZjfUlua1qZulTRGV8VNTa2hptbW3lLuOEePjZLXz27jXsG8jz+UVv4ONvm0tVVXautE6UiGDn3oGDAq1wW7Lws3P73oNGTuaqdGCE5MymwujIU5smHNie0VjP5PqaMp6R2fgiaXVEtI7U53soGRIR/MMjL/CVB9ZzzilT+PpHFnB6y6RylzVmSaKpoZamhlrOn9V0SP9gfohXd+1Lwi25Hbm9l8079vH4i928umsf+WFfC5hcX50EWyHc+geHmN3cwAcvmsWpTRNO0JmZmcMrIwbzQ3z+e09x7xOdXPnmU/nSB89nQq1H3R2P6lwVs5obmNXcAGcc2p8fCrp299G5Yy+vFC2dO/bxyo69PPHydnb0DgDwvx98jjNaJnLpWS1cOr+Fi0+fSkOt/3qZpcW3DTOgbzDPp+96kvvXvsafvussPv3OMzM1IKOS7ekb4IWuPTz+Yjc/e34rv9y4jb7BIWpzVbTObeaSM6ezcN5Uzp/V6CH+ZsfocLcNHV5j3EB+iKXfbOPh9V3c+P5z+fgl88pdkh3GvoE8q17q5ufPb+Vnz3Xx7KuFL4jXVlexYHYTF8+bysJ507jwtCZfmZkdgcMro4aGgs99dw0/WPMKf/O75/GRiz2LSNZs39PPqpe6+eWL3Tz+YjdrX9nJUEB1lXjjzEYWzG7igtlNvHl2E3OnNfiK2qyIwyuj/upH7dz+ixf5s3efzbW/fWa5y7ES2L1vgNW/3s7jL3bT9tJ2nu7cyd6BPACNE2p48+wmLpjVyJtnN3H+rKayf8ncrJw82jCD7l71Mrf/4kU+9ra5/MllI4wmsEyaXF/DZWefxGVnnwQUBuI8v6WHX23awZpk+frDXewf5NgyuY43nDKZc2dM4ZxkOb1lIjX+PpqNc77yGoN+tWkHH/7HR1k4byp3/sHCAzO82/jQ2z/IM527eKpjB+s2F54CsGFLD/35wjRZtbkq5p88iXNmTOGskydxRkthmdU8wV+ytoriK68M2dbTxx99azUtk+v4+6sXOLjGoYbaahbOm8rCea/PKDKQH2Jj1x7Wbd7Fus27aN+8i/98rot7Vncc2Kc2V8Xc6Q0HwuyMkyZyRsskTps2kcYJ/nK1VRaH1xgSEVx/79Ns6+nn3j95G1Mn1pa7JBsjanJVnH3KZM4+ZTIfWDDzQPvO3gFe2NrDC1t6eKFrDy909bD+td080P7aQV+wbpxQ8/pjZ4oeQTNnagMzmup9G9Iyx+E1hty9ahMPtr/G/3rvOaPOwG5WrLGhhgvnNHPhnOaD2vsHhwoz82/p4eXuPcnUWHtZt3kXD7S/ykB+hGmxGidwSmM9M5rqmTGlnlOSB4ie0ljP9Il1noLMxhSH1xjx0tY93Pyjdi45cxp/4O9y2XGqra7izJMmceZJh04flh8KXkumxXp9aqxeXtm5jzWbdvCTZ/Yd+Hxtv5qcOHlKPTMa6zl5Sj0tk+uYPqmOlsl1tEx6fX3apFpfxdkJ4fAaA4aGgs/f8xTVVeIrH36zf8O1VOWqdGDC4ZGeRBARbNvTz6s797F55z427yw8QPTVnYVpsZ7p3MnWnn56hs3Yv19zQw3ThwVa04RamifWFOaanFBDc0MtTQ01NDXUMKmu2t9vs2Pm8BoD7m7bxOMvdfPlD53PjEZP7mrlJelA+Bzu9vXe/jxbe/ro6ulj6+79P/sLbbv72NrTx1MdO9jW03/Io2mK1eRE44RCmDU31BxYn1xfzeS6aibXF9Yn1b++XtzeUJtz+I1DDq8y27JrH39z3zrecvpUPnzRrHKXY3bUJtTmDjzk80gG8kPs3DvAjt5+tvcOsH1PPzuKtnf09rOjd4Dtvf10bO9l7SsD7N43OOrVXbEqwaSiMJtYVwi0CTW5ws/awnZhPUdDTY6GuqK2mmH9tdVMqMlRW13l0b5jWKrhJWkR8H8oPIzy9oi4ZVh/HfBN4CJgG3BVRLyU9N0AfALIA5+OiPvTrLVcvvCjdvoGh/ib3z3Pvz1axarJVR24mjsWQ0NBT/8gPfsG2b1vkN37BtjdV1jvSbZ7ku1d+wqB19tfCL2u3X309ufp7c+zt3+Q3oE8x/q11pqcqKvOUVddVVhqitarc9TVFK1XVyXbudd/JvtWV4ma6ipqqqqoqRY1uSqqq6qorRbVVVXU5KqoySXtOVGbqzpovbqovybnUIUUw0tSDrgVuBzoAFZJWjnsicifALZHxJmSlgBfAq6SdC6wBHgjcCrwH5LOioh8WvWWw0/Xvcb/e2oz/+OKs/xcLrMRVFWJKfU1TCnBQ0Ajgn0DQ/T2DxYCbSCfhNsgew+EXJ49/YP0DQ7RNzBE32C+sD6YZ9/AUNL+eltP3yDbeor3K+4fOnJRvyGp8AtBTZXIHVgKIZkbvmiEtqS9OieqJKqrRFWVDj1+hH32v15VlagS5CS0v00k7WJm0wTe/+ZTU/tvkOaV10JgQ0RsBJC0AlgMFIfXYuCmZP0e4OsqXH4sBlZERB/woqQNyes9mlaxX/jhWtpf2ZXWyx8iPxSs27yLs06exNJLPf2TWdokMSG5NXjoMJXSiwj684UQG8wHA/mhZAkG80P054vbD+4fyA8xODTEwGAwMDTEwOAQg0Nx0DH71/NDhWVwKBja/zOKt4fID0F+aIh8JD+T9r7B5PiI118rXn/Ng5YI8vnX+yNgKArbI13RvvX0aZkNr5nApqLtDuDi0faJiEFJO4FpSftjw46dOexYJC0FliabPZLWH2fN04Gtx/kax6QdqPvciXzHE3+OZVDp51jp5wc+x8z7NbDik8d9jqeN1pHpARsRsQxYVqrXk9Q22jxalcLnmH2Vfn7gc6wUaZ5jmt8m7ARmF23PStpG3EdSNdBIYeDG0RxrZmbjVJrhtQqYL2mepFoKAzBWDttnJXBNsv4h4KEoTHO/ElgiqU7SPGA+8HiKtZqZWYakdtsw+QzrOuB+CkPll0fEWkk3A20RsRK4A/jXZEBGN4WAI9nvuxQ+EhoErj1BIw1LdgtyDPM5Zl+lnx/4HCtFaudYMc/zMjOz8cMzaJqZWeY4vMzMLHMcXglJiyStl7RB0vXlrqfUJC2XtEXSM+WuJQ2SZkt6WFK7pLWSPlPumkpNUr2kxyX9KjnHL5S7pjRIykl6UtKPyl1LGiS9JOlpSWsktZW7njRIapJ0j6RnJa2T9NaSv4c/8zowldVzFE1lBVw9bCqrTJN0KdADfDMi3lTuekpN0gxgRkQ8IWkysBr4QIX9GQqYGBE9kmqAXwCfiYjHjnBopkj6HNAKTImI95W7nlKT9BLQGhEV+wVlSXcCP4+I25PR5g0RsaOU7+Err4IDU1lFRD+wfyqrihERP6MworMiRcTmiHgiWd8NrGOEWVmyLAp6ks2aZKmo3z4lzQLeC9xe7lrsNyOpEbiUwmhyIqK/1MEFDq/9RprKqqL+4RtPJM0FFgC/LHMpJZfcUlsDbAEejIhKO8evAZ8H0pvVtvwCeEDS6mSKu0ozD+gC/jm5/Xu7pImlfhOHl1UUSZOA7wGfjYgTN9PyCRIR+Yi4gMKsMwslVcwtYEnvA7ZExOpy15Ky34qIC4H3ANcmt/QrSTVwIfCNiFgA7AFKPo7A4VXg6agqQPI50PeAb0fEveWuJ03JbZiHgUVlLqWULgGuTD4TWgG8Q9K3yltS6UVEZ/JzC/B9Ch9bVJIOoKPorsA9FMKspBxeBUczlZWNYclghjuAdRHxd+WuJw2SWiQ1JesTKAwwerasRZVQRNwQEbMiYi6Fv4MPRcRHy1xWSUmamAwoIrmVdgVQUSOAI+JVYJOks5Omd3Lwo7BKItOzypfKaFNZlbmskpJ0F3AZMF1SB3BjRNxR3qpK6hLg94Cnk8+EAP48Iu4rX0klNwO4MxkdWwV8NyIqcjh5BTsZ+H7y1PRq4DsR8ZPylpSKTwHfTi4GNgIfL/UbeKi8mZlljm8bmplZ5ji8zMwscxxeZmaWOQ4vMzPLHIeXmZlljsPLrMwkTUtmGF8j6VVJncl6j6R/KHd9ZmORh8qbjSGSbgJ6IuIr5a7FbCzzlZfZGCXpsv3PtJJ0k6Q7Jf1c0q8l/XdJX06eC/WTZGosJF0k6T+TSV/vTx4VY1ZxHF5m2XEG8A7gSuBbwMMRcR6wF3hvEmD/F/hQRFwELAf+ulzFmqXJ00OZZcePI2JA0tMUpjHbP63Q08Bc4GzgTcCDyfRDOWBzGeo0S53Dyyw7+gAiYkjSQLz+gfUQhb/LAtZGRMkfuW421vi2oVnlWA+0SHorFB4RI+mNZa7JLBUOL7MKERH9wIeAL0n6FbAGeFtZizJLiYfKm5lZ5vjKy8zMMsfhZWZmmePwMjOzzHF4mZlZ5ji8zMwscxxeZmaWOQ4vMzPLnP8PsdOZgE5qRrgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplot( 2, 1, 1)\n",
    "plt.plot( [sol1.t[0],time_x_on,time_x_on,time_x_off,time_x_off,sol3.t[-1]], [0,0,1,1,0,0] )\n",
    "plt.ylabel( \"X\")\n",
    "plt.ylim([0,1.1])\n",
    "\n",
    "plt.subplot( 2, 1, 2 )\n",
    "plt.plot(sol1.t, sol1.y[1,:].T); #Simulation before signal on\n",
    "plt.gca().set_prop_cycle(None) # this resets the color cycle\n",
    "plt.plot(sol2.t, sol2.y[1,:].T); #Simulation after signal on\n",
    "plt.gca().set_prop_cycle(None) # this resets the color cycle\n",
    "plt.plot(sol3.t, sol3.y[1,:].T); #Simulation after signal on\n",
    "plt.xlabel( 'Time')\n",
    "plt.ylabel( \"Z\")\n",
    "plt.ylim([0,1.1])\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuzklEQVR4nO3deXxV9Z3/8dcnCyQsYQ17IGgRRbRsgyjVoojiSh2x7tQuQxcptp1pf7XTGe0yVR+2Y8dqq4gitloX3KjFrSpYx42AqcgSjBElbAkBAwlk//z+yM2dS8jGcu+5N/f9fDzO42zfnPs5JOSdc873fq+5OyIiIvEmJegCREREWqKAEhGRuKSAEhGRuKSAEhGRuKSAEhGRuJQWdAGHqn///p6bmxt0GSIicphWrVq1092z22uXcAGVm5tLXl5e0GWIiMhhMrNPOtJOt/hERCQuRS2gzOwBMysxsw9a2W9mdqeZFZrZ+2Y2IVq1iIhI4onmFdSDwMw29p8HjApNc4E/RLEWERFJMFF7BuXur5tZbhtNZgEPeeNYS2+bWW8zG+zu26JV05FYsGABjzzySNBlSIxcddVVzJ07N+gyRJJakM+ghgKbI9aLQ9sOYmZzzSzPzPJKS0tjUlxzjzzyCPn5+YG8tsRWfn6+/hgRiQMJ0YvP3RcACwAmTZoU2Oi248aNY/ny5UG9vMTItGnTgi5BRAj2CmoLkBOxPiy0TUREJNCAWgrMCfXmmwKUx+vzJxERib2o3eIzsz8D04D+ZlYM3ASkA7j7PcAy4HygENgHfDVatYiISOKJZi++K9vZ78D10Xp9ERFJbBpJQkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4lKrAWVmw9vYd3p0yhEREWnU1hXUcjP7kZmlNm0ws4Fm9ifgjo4c3MxmmlmBmRWa2Y9b2H+dmZWaWX5o+sahn4KIiHRGbQXUROBYIN/MzjKzG4B3gbeAye0dOBRsdwPnAWOAK81sTAtNH3P3caFp4SGfgYiIdEppre1w993AN0PB9DdgKzDF3Ys7eOzJQKG7FwGY2aPALGDdkZUsIiLJoK1nUL3N7F7gq8BMYAnwvJmd1cFjDwU2R6wXh7Y1d6mZvW9mS8wsp4PHFhGRTq6tW3yrgQ+BSe7+krt/D7gW+KWZ/fkovf5fgFx3Pxl4GVjcUiMzm2tmeWaWV1paepReWkRE4llbAXWGu//a3euaNrh7vrufBrzagWNvASKviIaFtoW5e5m7V4dWF9L43Osg7r7A3Se5+6Ts7OwOvLSIiCS6VgOqrWdN7n5fB469EhhlZiPNrAtwBbA0soGZDY5YvRhY34HjiohIEmi1k8SRcvc6M5sHvAikAg+4+1oz+zmQ5+5LgflmdjFQB+wCrotWPSIikliiFlAA7r4MWNZs239GLN8I3BjNGkREJDFpqCMREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLCigREYlLUQ0oM5tpZgVmVmhmP25hf1czeyy0/x0zy41mPSIikjiiFlBmlgrcDZwHjAGuNLMxzZp9Hdjt7p8D7gBui1Y9IiKSWNKieOzJQKG7FwGY2aPALGBdRJtZwM2h5SXAXWZm7u6tHbSgoIBp06ZFpeC25OfnM27cuJi/rgQjPz8/kJ8zEfk/0QyoocDmiPVi4JTW2rh7nZmVA/2AnZGNzGwuMDe0Wr1ixYoPolJxO1asWIGZBfHSAP1p9u+SJAI77xUrVgTxsqDvdbJJxvMe3ZFG0Qyoo8bdFwALAMwsz90nBVxSzOm8k0cynjPovIOuI5bMLK8j7aLZSWILkBOxPiy0rcU2ZpYG9ALKoliTiIgkiGgG1EpglJmNNLMuwBXA0mZtlgJfCS3PBl5t6/mTiIgkj6jd4gs9U5oHvAikAg+4+1oz+zmQ5+5LgfuBP5pZIbCLxhBrz4Jo1RzndN7JIxnPGXTeyaRD52y6YBERkXikkSRERCQuKaBERCQuJVRAtTd0UmdkZg+YWYmZBfLeryCYWY6ZvWZm68xsrZndEHRNsWBmGWb2rpn9I3TePwu6plgxs1Qze8/Mngu6llgxs01mtsbM8jva7bozMLPeZrbEzDaY2XozO7XVtonyDCo0dNJGYAaNb/pdCVzp7uva/MIEZ2ZnABXAQ+4+Nuh6YsHMBgOD3X21mfUEVgFfSoLvtQHd3b3CzNKBN4Ab3P3tgEuLOjP7ATAJyHL3C4OuJxbMbBMwyd2T6k26ZrYY+Lu7Lwz18O7m7p+11DaRrqDCQye5ew3QNHRSp+bur9PYwzFpuPs2d18dWt4LrKdx1JFOzRtVhFbTQ1Ni/AV5BMxsGHABsDDoWiS6zKwXcAaNPbhx95rWwgkSK6BaGjqp0//SSnahEe7HA+8EXEpMhG515QMlwMvungzn/VvgR0BDwHXEmgMvmdmq0HBuyWAkUAosCt3SXWhm3VtrnEgBJUnGzHoATwLfc/c9QdcTC+5e7+7jaBx5ZbKZderbumZ2IVDi7quCriUAX3D3CTR+4sP1odv5nV0aMAH4g7uPByqBVvsTJFJAdWToJOkkQs9gngQedvengq4n1kK3PV4DZgZcSrRNBS4OPY95FDjLzP4UbEmx4e5bQvMS4GkaH2N0dsVAccSdgSU0BlaLEimgOjJ0knQCoc4C9wPr3f2/g64nVsws28x6h5YzaewQtCHQoqLM3W9092Hunkvj/+lX3f2agMuKOjPrHuoAROgW1zlAp++p6+7bgc1m1jSa+XQO/AimAyTEaObQ+tBJAZcVdWb2Z2Aa0N/MioGb3P3+YKuKuqnAtcCa0PMYgJ+4+7LgSoqJwcDiUI/VFOBxd0+abtdJZiDwdOjje9KAR9z9hWBLipnvAg+HLjSKgK+21jBhupmLiEhySaRbfCIikkQUUCIiEpcUUCIiEpcUUCIiEpcUUCIiEpcUUCIiEpcUUCIiEpcUUCIiEpcUUCIiEpcUUCIiEpcUUCIiEpcSZrDYJv379/fc3NygyxARkcO0atWqne6e3V67hAuo3Nxc8vLygi5DREQOk5l90pF2UbvFZ2YPmFmJmbX4GSfW6E4zKzSz982s1Q+tEhGR5BPNZ1AP0vangZ4HjApNc4E/RLEWERFJMFELKHd/HdjVRpNZwEPe6G2gt5kNjlY9IiKSWILsxTcU2ByxXhzaJiIikhjdzM1srpnlmVleaWlp0OWIiEgMBBlQW4CciPVhoW0HcfcF7j7J3SdlZ7fbM1FERDqBIANqKTAn1JtvClDu7tsCrEdEROJI1N4HZWZ/BqYB/c2sGLgJSAdw93uAZcD5QCGwD/hqtGoREZHEE7WAcvcr29nvwPXRen0REUlsCdFJQkREko8CSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4pICSkRE4lKrAWVmw9vYd3p0yhEREWnU1hXUcjP7kZmlNm0ws4Fm9ifgjuiXJiIiyaytgJoIHAvkm9lZZnYD8C7wFjC5Iwc3s5lmVmBmhWb24xb2X2dmpWaWH5q+cTgnISIinU9aazvcfTfwzVAw/Q3YCkxx9+KOHDh05XU3MAMoBlaa2VJ3X9es6WPuPu+wqhcRkU6rrWdQvc3sXuCrwExgCfC8mZ3VwWNPBgrdvcjda4BHgVlHWrCIiCSHtm7xrQY+BCa5+0vu/j3gWuCXZvbnDhx7KLA5Yr04tK25S83sfTNbYmY5LR3IzOaaWZ6Z5ZWWlnbgpUVEJNG1FVBnuPuv3b2uaYO757v7acCrR+n1/wLkuvvJwMvA4pYaufsCd5/k7pOys7OP0kuLiEg8azWg2nrW5O73deDYW4DIK6JhoW2Rxylz9+rQ6kIaO2aIiIhE9Y26K4FRZjbSzLoAVwBLIxuY2eCI1YuB9VGsR0REEkirvfiOlLvXmdk84EUgFXjA3dea2c+BPHdfCsw3s4uBOmAXcF206hERkcRi7h50DYdk0qRJnpeXF3QZIiJymMxslbtPaq+dxuITEZG4pIASEZG4pIASEZG4pIASEZG4pIASEZG4pIASEZG4pIASEZG4FLU36sqR2717Nx999BFFRUVs27aNkpKS8LRr1y4qKyvZt28flZWVVFZWUlNTA0DTe9ua5ikpKWRkZLQ69ezZk6ysLHr16nXAvPly37596devH5mZmYH9m4hI8lBAxYH6+nrWrl3L6tWrWbVqFatXr2bdunV89tlnR+01Kisrj9qxMjIy6NevXziwIpfbmnfp0uWo1SAinZ8CKiBbtmzhueee429/+xuvvvoqu3btOqhN9+7dOfbYYxk5ciTDhg1j4MCBDBgwgAEDBtC3b1+6d+9+wNSlSxfMDOCAeX19PdXV1VRVVYWn/fv3h+d79+5lz5497Nmzh/Ly8haXP/vsM3bv3k1ZWRlVVVVs2bKFLVu2HFRzW3r06BEOtI5OWVlZ4XMRkeSigIqhsrIyHn/8cR599FH+/ve/EznM1PDhwznllFOYMGECEydO5KSTTmLgwIFH7Zdzz549j8px3J19+/ZRVlbGrl272p03LZeVlVFRUUFFRQWffPJJh18vLS3tgCu1lqaW9utqTSTxKaBiYM2aNdx555386U9/oqqqCoCuXbty3nnnMXPmTM4++2yOOeaYhLhSMLPwFdvw4cM7/HXuzp49e8Jh1dGpoqIi/NztUHT0aq1Pnz706dOH3r1707t3b7p06YK7J8T3QqSzU0BF0bvvvst//Md/8NJLL4W3nXPOOVxzzTXMmjWLrKysAKuLLTOjV69e9OrVi2OOOabDX1dTU3PAVdihBNuhXq0BdOvWjZSUFNydM844g9GjR3P88ceH50fzqlZE2qbRzKOgsLCQH/7whzzzzDNA41/z1113HfPmzWP06NHBFpcEDuVqrenZWtO8vr6+zWNnZWVx/PHHh0Prc5/7HMcccwzHHHMMffr0UXiJdEBHRzNXQB1F+/fv55ZbbuG2226jpqaGzMxM5s+fzw9/+EP69esXdHnSDnensrKSHTt2sGnTJvbu3cuGDRsoKChgw4YNbNiwoc2elU1Xh82nkSNHkpOTQ0ZGRuxORiSOKaBibOXKlVxzzTVs3LgRgK985SvceuutDBo0KODK5Ghxd0pLS8OhVVBQQFFREUVFRXz00UdUVFS0+fXZ2dkMGzaMnJycg+Y5OTkMHTqUrl27xuhsRIKjgIqR+vp6fvWrX/Gzn/2M+vp6xowZwz333MPpp59OVVUVFRUV7Nu3j/3797N//37q6upoaGigvr6ehoYGGhoaSElJITU1ldTU1AOW09PT6dKlS3jeNKWnp5Oenq7bSXHE3SkrKwsHVvNpy5Yt1NXVtXucAQMGMHjwYAYOHMigQYPCU/N13U6URKaAioHdu3dz1VVX8cILLwCNV02zZ8+mpKSEsrIyamtrycrKIjMzk27dupGRkUFaWlo4hFJSUkhJSQkHVmRo1dXVUVdXR01NTXiqra0NL9fV1R0UXJFT165dSU9Pp2vXrgdsa94mcj0tLU2/9KKkvr6ekpISNm/ezObNmykuLj5oeevWre0+A2uSnp7OwIEDGThwIP3796dfv37heeRy5LZu3bpF+SxFOkYBFWVr165l1qxZfPTRR/Ts2ZPLLruMadOmMXz4cHJycsjOzqZHjx5R+4Xf0NAQDqza2lqqq6vD86YQi1yOnJpvb1pvaGg4KLiah1hHw64pIFNSNNxjR9XX17N9+3a2b9/Ojh072lwuLy8/5ONnZmaGw6p379706tWrQ/OmZb23TI6WjgaUupkfhueff57Zs2ezb98+hg4dyu9//3tmzJgR0zHqUlJS6Nq161F9ZlFfX99uiDVN5eXlHQq/yKu85iGWnp5OWlpa+JZl0/KhzlNSUjrFlV9qaipDhw5l6NCh7bbdv38/JSUlbN++PdwjcefOnQctR873799PcXExxcXFh1VfZmYmWVlZ9OjRIzz17NnzgPWOTN26dSMzMzM86Y8YaU1UA8rMZgL/A6QCC9391mb7uwIPAROBMuByd98UzZqO1OLFi/n6179OfX0955xzDkuWLDlqozQELTU1NfxL42hw94Ou6iJDrLa2ltraWurq6sLL+/btC9/ejNzX2ry2thbgoNBKTU0lLS2tzXlr+wDq6uro2bMnaWlpLU5NrxU5xfIXbWZmJiNGjGDEiBEdat/UQ7EpwMrLy/nss88Oad70HHXHjh1H9Vy6dOlyQGA1Tc2DrKUpIyPjgD9+Wrr678g+3d6OT1ELKDNLBe4GZgDFwEozW+ru6yKafR3Y7e6fM7MrgNuAy6NV05Fwd2677TZuvPFGAL7//e/zm9/8Rj/UbTCz8C+EaKqvrz8guJqmpu0dnVdXV1NfX09VVRW1tbWkp6cfcLyWXiNyAloNr8hQbN4ZJnK5tTbtbW/aZmbhZ5uRU9P2/v37M2DAgIP2NU2taRriqry8nMrKyvAboZtPe/fubXVf0/7ITkNVVVUHXJUHpelnNTLEIq/sW/q+Hu19kd/H5vMg9jX9zLQ2j8Wdi6g9gzKzU4Gb3f3c0PqNAO5+S0SbF0Nt3jKzNGA7kO1tFGVmgT80+853vsNll10W1SFxIn9pxGo61Ndtat/SvK19idDmcDQFXWtTU0eYyCmyg0xb25pvb+3rmiZ3P2C9ve3AQb98Wgq41rZ39GempUBsuhJumiKvriM7CEWuR25v/odD8yvspvWmzkXNX6+2trbDnVPkYM2/t5Hf49aCraysLNhOEmY2G5jp7t8IrV8LnOLu8yLafBBqUxxa/yjUZmezY80F5oZWxwIfRKXo+NYf2Nluq84nGc87Gc8ZdN7JZLS7t/tsJCE6Sbj7AmABgJnldSR5Oxudd/JIxnMGnXfQdcSSmXWoK3Y0n+puAXIi1oeFtrXYJnSLrxeNnSVERCTJRTOgVgKjzGykmXUBrgCWNmuzFPhKaHk28Gpbz59ERCR5RO0Wn7vXmdk84EUau5k/4O5rzeznQJ67LwXuB/5oZoXALhpDrD0LolVznNN5J49kPGfQeSeTDp1zwo0kISIiyUFv4RYRkbikgBIRkbiUUAFlZjPNrMDMCs3sx0HXEwtm9oCZlYTeM5YUzCzHzF4zs3VmttbMbgi6plgwswwze9fM/hE6758FXVOsmFmqmb1nZs8FXUusmNkmM1tjZvkd7XbdGZhZbzNbYmYbzGy9NQ7q0HLbRHkGZY1DJ20kYugk4MpmQyd1OmZ2BlABPOTuY4OuJxbMbDAw2N1Xm1lPYBXwpST4XhvQ3d0rzCwdeAO4wd3fDri0qDOzHwCTgCx3vzDoemLBzDYBk5oPTNDZmdli4O/uvjDUw7ubu3/WUttEuoKaDBS6e5G71wCPArMCrinq3P11Gns4Jg133+buq0PLe4H1QPtDfCc4b9T0sbzpoSkx/oI8AmY2DLgAWBh0LRJdZtYLOIPGHty4e01r4QSJFVBDgc0R68UkwS+tZGdmucB44J2AS4mJ0K2ufKAEeNndk+G8fwv8CGgIuI5Yc+AlM1sVGs4tGYwESoFFoVu6C82se2uNEymgJMmYWQ/gSeB77r4n6Hpiwd3r3X0cjSOvTDazTn1b18wuBErcfVXQtQTgC+4+ATgPuD50O7+zSwMmAH9w9/FAJdBqf4JECqiODJ0knUToGcyTwMPu/lTQ9cRa6LbHa8DMgEuJtqnAxaHnMY8CZ5nZn4ItKTbcfUtoXgI8TeNjjM6uGCiOuDOwhMbAalEiBVRHhk6STiDUWeB+YL27/3fQ9cSKmWWbWe/QciaNHYI2BFpUlLn7je4+zN1zafw//aq7XxNwWVFnZt1DHYAI3eI6hyT4lAZ33w5sNrPRoU3TgVY7PyXEaObQ+tBJAZcVdWb2Z2Aa0N/MioGb3P3+YKuKuqnAtcCa0PMYgJ+4+7LgSoqJwcDiUI/VFOBxd0+abtdJZiDwdOgzsdKAR9z9hWBLipnvAg+HLjSKgK+21jBhupmLiEhySaRbfCIikkQUUCIiEpcUUCIiEpcUUCIiEpcUUCIiEpcUUCIiEpcUUCIiEpcUUCIiEpcUUCIiEpcUUCIiEpcUUCIiEpcSZrDYJv379/fc3NygyxARkcO0atWqne6e3V67hAuo3Nxc8vLygi5DREQOk5l90pF2usUnIiJxKWoBZWYPmFmJmbX4IVzW6E4zKzSz982s1U9VFBGR5BPNK6gHafvjqs8DRoWmucAfoliLiIgkmKgFlLu/Duxqo8ks4CFv9DbQ28wGR6seERFJLEE+gxoKbI5YLw5tO4iZzTWzPDPLKy0tjUlxIiISrIToJOHuC9x9krtPys5ut2eiiIh0AkEG1BYgJ2J9WGibiIhIoAG1FJgT6s03BSh3920B1iMiInEkam/UNbM/A9OA/mZWDNwEpAO4+z3AMuB8oBDYB3w1WrWIiEjiiVpAufuV7ex34Ppovb6IiCS2hOgkISIiySfhxuITiYaGhga2bt1KUVEReXl5vP766zz++ON06dIl6NJEkpYCSpKGu7Nt2zYKCwv58MMP2bBhA+vXr2fjxo18+umnVFdXH9D+tNNO49lnn2Xo0BbfniciUaaAkk6nrq6OwsJC1q5dywcffEB+fj7r16/nk08+oaqqqtWv6927N7m5ufTp04eVK1eyatUqTjrpJJ566immTZsWuxMQEUABJQmsoaGBoqIi1qxZwwcffMB7773HBx98wMcff0xdXV2LX9OzZ0+OOeYYRo8ezZgxYzjuuOMYNWoUn/vc5+jdu3e4XWlpKf/8z//MG2+8wfTp07n11lv5t3/7N8wsRmcnItbYmS5xTJo0yfV5UMmnrq6ODRs2sHr1at59911WrlzJBx98wL59+1psP3DgQMaMGcO4ceM4+eSTOeGEExg1ahR9+/Y9pNe88cYb+fWvfw3Al770JRYvXkxWVtZROSeRZGVmq9x9UrvtFFASbxoaGigoKOCtt97izTffZOXKlWzYsIGampqD2vbr14+xY8eGg+jEE09kzJgx9OzZ86jV89RTT3Httdeyb98+hg4dypIlS5gyZcpRO75IslFAScIoLy/nnXfe4X//939ZsWIFq1atoqKi4qB2Q4YMYcKECUyZMoUJEyYwfvx4Bg0aFJMaN27cyCWXXMK6detISUnh5ptv5ic/+QmpqakxeX2RzkQBJXFr+/btLF++nJdffpkVK1ZQVFRE85/D7OxspkyZwumnn86kSZMYN24cffr0CajiRtXV1fzoRz/izjvvBODUU0/liSeeUC8/kUOkgJK4sWPHjnAgvfLKK2zatOmA/WlpaYwdO5YzzjiDqVOncuqpp5KTk9PyweLA888/z9VXX83u3bvp2bMn99xzD1deeaU6UIh0kAJKAlNVVcXrr7/OX//6V5577jmKiooO2J+RkcEpp5zCueeeyxe/+EUmTJhARkZGQNUenh07dnD11VfzyiuvAHDBBRdw//33M3DgwIArE4l/CiiJGXensLCQ559/nmeeeYY333zzgDe9du3alcmTJzNz5kzOPPNMJk2aRHp6eoAVHx3uzoIFC/j+97/P/v37ycrK4p577uGKK67Q1ZRIGxRQElX19fW8/fbbLFmyhCVLllBcXHzA/hNOOIGLLrqI888/n1NPPbVTDxn06aefMmfOHFasWAE0Xk3de++9ejYl0goFlBx11dXVvPrqqzz++OMsXbqUXbt2hfdlZWUxffp0Zs2axbnnnhuz3nXxovnVVGZmJv/1X//Fd7/7XdLS9H54kUgKKDkqamtreemll3jwwQdZtmzZAW+MHTJkCJdeeimXXXYZp512mrpcA5s3b+Zb3/oWy5YtA2DMmDEsWrSIyZMnB1yZSPxQQMlha2ho4I033uChhx7iiSeeYM+ePeF9xx13HJdffjmXXnopJ598sp61tOLZZ5/lW9/6Ftu3b8fM+NrXvsatt95K//79gy5NJHAKKDlkBQUFLFy4kIceeoiSkpLw9mOPPZbrrruOK6+8kmOPPTbAChNLZWUlN998M3fccQf19fV0796dm2++mfnz53fqZ3Ii7VFASYdUVlby+OOP8/vf/57If9fBgwdzzTXXcM0113DSSSfpSukIrF27luuvvz7ciSInJ4ff/va3XHLJJfp3laSkgJI2rVy5kj/84Q889thj4edKGRkZfPnLX+Zf/uVfOO2000hJ0QcuHy3uzvPPP8+8efP4+OOPgcaRKO68804mTWr3/6lIp9LRgGrzN5CZndvGvssOpzAJTnV1NX/84x/5/Oc/z+TJk1m0aBH79u1j/PjxLFy4kNLSUhYvXswXvvAFhdNRZmacf/75FBQU8Lvf/Y6srCzeeust/umf/okLLriANWvWBF2iSPxx91YnoB54DRjawr7VbX1tqM1MoAAoBH7cwv7rgFIgPzR9o71jTpw40eXQbN261f/93//d+/bt64AD3rNnT58/f76vXbs26PKS0q5du/xf//VfvUuXLg64mfmll17qBQUFQZcmEnVAnrfzu97d2w2o94B/ATYDs5vva+drU4GPgGOALsA/gDHN2lwH3NWRQpsmBVTHrV692mfPnu2pqanhYBo9erTfd999XllZGXR54u7btm3zb3/7256WlhYOqquvvlpBJZ1aRwOqvfs47u73AdOB/2dmi8ysW9O+dr52MlDo7kXuXgM8Csxq52vkCLk7y5cv58wzz2TChAksWbKEhoYGLrroIpYvX8769ev5xje+Qbdu3do/mETdoEGD+P3vf89HH33Eddddh5nx8MMPc/zxxzNr1ixWr14ddIkigenQgwZ33wicCuwA3jOzUzrwZUNpvPJqUhza1tylZva+mS0xsxaHsDazuWaWZ2Z5paWlHSk56TQ0NPDss88yceJEzjzzTJYvX05GRgbz589n06ZNLF26lC9+8YvqNRanhg8fzqJFi9i4cSPXXXcdqampLF26lIkTJzJt2jSWL19+0EeSiHR6bV1e0cJtPGAaUATsbedrZwMLI9avpdntPKAf0DW0/E3g1fYu+XSL70ANDQ3+5JNP+qhRo8K38bKysvw///M/vaysLOjy5DAVFxf7DTfc4BkZGeHv69ixY33x4sVeVVUVdHkiR4Sj9AzqS61s70MLnR6atTkVeDFi/UbgxjbapwLl7RWsgGrU0NDgf/3rX/3EE08M/wIbMGCA33HHHV5RURF0eXKUlJWV+U033eRZWVnh73Pfvn39pz/9qW/dujXo8kQOy1EJqCOZgLTQldZI/q+TxInN2gyOWL4EeLu94yqg3F955RWfMGFC+BdW//79/Xe/+53+su7E9u3b5wsXLvTRo0eHv+9paWk+e/ZsX7FihTc0NARdokiHBR5QjTVwPrCRxt58/x7a9nPg4tDyLcDaUHi9Bhzf3jGTOaDWrFnj06ZNC/+C6tWrl99+++2+b9++oEuTGGloaPAVK1b4hRde6GYW/lkYMWKE33LLLb5t27agSxRpV0cDSiNJJICSkhJuvPFGFi1ahLvTvXt3fvzjH/O9732PHj16BF2eBOSTTz7h3nvv5b777mPnzp0ApKSkMHPmTL7zne9w7rnn6qM+JC5pqKNOoKqqijvuuINf/OIX7N+/n5SUFObOncsvfvELjYotYXV1dbz44ovcfffdvPjiizQ0NADQp08frrjiCubMmcMpp5yiHpwSNxRQCe65557jW9/6Flu2bAHg7LPP5ne/+x3HH398wJVJPNu+fTuLFi1iwYIFbNq0Kbx96NChzJkzhzlz5uhnSAKngEpQn376Kd/+9rfDH3h33HHHcddddzFjxoyAK5NE4u689957LF68mIcffpiysrLwvlGjRnH55Zcze/ZsfaaXBEIBlWBqa2v5zW9+w80330x1dTWZmZn88pe/ZP78+XqOIEekvr6eFStWsHjxYp588kkqKyvD+4YMGcKXv/xlZs+ezamnnqpBgiUmFFAJJC8vj6uvvpqNGzcCMGvWLO6++26GDm1p4A2Rw1dTU8Nrr73GY489xjPPPMPu3bvD+/r27cv555/PxRdfzIwZM+jdu3dwhUqnpoBKANXV1dx0003cfvvtNDQ0MGzYMO677z5mzpwZdGmSBOrr63nrrbdYsmQJTzzxBFu3bg3vS0lJYcKECVx88cVccMEFjBs3TldXctQooOLcypUrueqqqygsLMTMmDdvHrfeeqsGcZVAuDtr1qxh2bJlPPPMM+Tl5VFfXx/e37dvX84++2zOPfdczjzzTEaOHBlgtZLoFFBxqra2lp/+9KfcfvvtuDsjRozg4YcfZurUqUGXJhK2Z88e/va3v/GXv/yFZcuWUVJScsD+QYMGMX36dGbMmMGZZ57J8OHDA6pUEpECKg4VFhZy6aWX8v7772NmzJ8/n1/96le6apK45u6sXbuWV155hRdffJE33niDvXv3HtBmyJAhTJ06lTPOOIOpU6dy0kknqXOPtEoBFUfcnQcffJDrr7+e/fv3M2jQIB5//HFOP/30oEsTOWQNDQ384x//4LXXXuOFF17gzTffPKBnIEBGRgbjxo3ji1/8IlOnTmXKlClkZ2cHVLHEGwVUnCgvL+frX/86Tz75JACXXHIJDzzwgHpISadRX1/P+++/z5tvvsnrr7/Om2++SXFx8UHtBg4cyIQJE5gyZQoTJ05kwoQJDB48OICKJWgKqDiQn5/PRRddRHFxMRkZGdx111187Wtf0xsjpdMrKSnh7bff5o033uD1118nPz+f6urqg9r17duXz3/+80yZMoXx48dz4oknMmrUKNLT0wOoWmJFARWwBx98kG9+85vU1NRw4okn8tRTT3HccccFXZZIIOrr69m4cSPvvfceK1eu5J133uH9998/6NYgQFpaGrm5uYwdO5bx48czduxYTjzxRI499lg91+okFFABqa6uZt68eSxcuBCAOXPmcO+995KRkRFwZSLxxd35+OOPWb16Ne+++y6rV69m3bp1bNu2rcX26enp5Obmctxxx3HCCSdw/PHHM2rUKEaNGsWgQYN0ZyKBKKACUFxczEUXXUR+fj7p6encddddzJ07N+iyRBJKRUUF69evZ+3atbz//vvk5+ezbt06duzY0erXZGRkMGLECI477jjGjBnDqFGjyM3NZcSIEeTk5NC1a9cYnoG0RwEVY++88w7nn38+u3btYvDgwfzlL39h4sSJQZcl0mns2bOHgoICCgsL2bhxI+vWraOgoICioqKDur1HMjP69etHTk4Oubm5jBo1ipEjRzJixAhGjBjBsGHDyMrKiuGZiAIqhh599FHmzJlDbW0tp59+Ok8//TT9+vULuiyRpLFr1y4KCwv58MMPKSgo4MMPP6SoqIhPP/2UkpKS8GdktSYzM5Ps7GyGDBnCsGHDGD58ODk5OQwZMoShQ4cyZMgQBg8erFv1R4kCKgbcnZtvvpmf//znAHzta1/jnnvuUQ8kkThSV1fHli1b+OSTT9i0aRObNm3iww8/5OOPP+aTTz6htLS0xR6GLenZsyf9+vUjOzubgQMHMnDgQIYMGcKAAQMYMGAA2dnZ4Xm/fv1ITU2N8tklJgVUlFVVVTFnzhyeeOIJzIzbb7+dH/zgB3pQK5Jg3J3y8nK2bt0anrZs2cKnn37K5s2bKS4uZvv27ezcufOA8QnbY2ZkZWXRq1cv+vTpQ9++fenXrx/9+/cnOzubvn37hrdHzvv06dPpr9Q6GlDqs3kYdu/ezcyZM3n33XfJzMzkscce46KLLgq6LBE5DGZG79696d27N2PGjGm1XUNDA2VlZZSWllJSUhKe79ixg23btrF9+3ZKSkooKSmhrKyMvXv3Ul5eTnl5OZ9++ukh1dS1a1e6d+9Ojx496NGjB7169SIrK4usrKxwrb169aJXr1707NmTrKysg+Y9evSgW7duCX1HRwF1iLZs2cJZZ53Fxo0bGTBgAC+//DInn3xy0GWJSJSlpKSQnZ1NdnZ2m0HWpLa2lrKyMnbt2sXu3bvZvXt3eLmsrIydO3eyc+dOysrKwvvLy8vZu3cv1dXVVFdXs2vXriOuOy0tja5du5KRkUFmZibdunWjW7dudO/ePTxvCsKmKXJfRkZGeGo6TtO8pW1H87ZmVAPKzGYC/wOkAgvd/dZm+7sCDwETgTLgcnffFM2ajkRBQQFnnnkm27Zt49hjj+W1114jJycn6LJEJA6lp6czaNAgBg0adEhf5+5UVlayZ88e9uzZw969ew+YNy1/9tlnfPbZZ+GrtD179lBeXk5FRQUVFRXs37+fqqoq6urqqKura/FN0dGQmppKeno6Xbp0oUuXLnTt2jU8NYVaR0UtoMwsFbgbmAEUAyvNbKm7r4to9nVgt7t/zsyuAG4DLo9WTUdi5cqVzJgxg/LycsaPH8/LL7+snnoictSZWfhKZsiQIUd0LHentraWffv2UVlZecC8tW1NARcZdPv376e6upqqqiqqqqrCyzU1NVRXV1NbW0tNTQ01NTXU19dTX19PVVXVkf9bRKuThJmdCtzs7ueG1m8EcPdbItq8GGrzlpmlAduBbG+jKDMLtFfH9OnTefbZZ+nevXuQZYiIxB13p66u7oAQi1xumk+fPj3wThJDgc0R68XAKa21cfc6MysH+gE7IxuZ2VygaUiGauCDaBTcEa+88go9evQI4qX70+zfJUkk43kn4zmDzjuZjO5Io4ToJOHuC4AFAGaW15Hk7Wx03skjGc8ZdN5B1xFLZtah9wqlRLGGLUBkD4JhoW0ttgnd4utFY2cJERFJctEMqJXAKDMbaWZdgCuApc3aLAW+ElqeDbza1vMnERFJHlG7xRd6pjQPeJHGbuYPuPtaM/s5kOfuS4H7gT+aWSGwi8YQa8+CaNUc53TeySMZzxl03smkQ+eccEMdiYhIcojmLT4REZHDpoASEZG4lFABZWYzzazAzArN7MdB1xMLZvaAmZWYWWDv/Yo1M8sxs9fMbJ2ZrTWzG4KuKRbMLMPM3jWzf4TO+2dB1xQrZpZqZu+Z2XNB1xIrZrbJzNaYWX5Hu113BmbW28yWmNkGM1sfGtSh5baJ8gwqNHTSRiKGTgKubDZ0UqdjZmcAFcBD7j426HpiwcwGA4PdfbWZ9QRWAV9Kgu+1Ad3dvcLM0oE3gBvc/e2AS4s6M/sBMAnIcvcLg64nFsxsEzDJ3ZPqTbpmthj4u7svDPXw7ubun7XUNpGuoCYDhe5e5O41wKPArIBrijp3f53GHo5Jw923ufvq0PJeYD2No450at6oIrSaHpoS4y/II2Bmw4ALgIVB1yLRZWa9gDNo7MGNu9e0Fk6QWAHV0tBJnf6XVrIzs1xgPPBOwKXEROhWVz5QArzs7slw3r8FfgS0/bnsnY8DL5nZqtBwbslgJFAKLArd0l1oZq0ObJpIASVJxsx6AE8C33P3PUHXEwvuXu/u42gceWWymXXq27pmdiFQ4u6rgq4lAF9w9wnAecD1odv5nV0aMAH4g7uPByqBVvsTJFJAdWToJOkkQs9gngQedvengq4n1kK3PV4DZgZcSrRNBS4OPY95FDjLzP4UbEmx4e5bQvMS4GkaH2N0dsVAccSdgSU0BlaLEimgOjJ0knQCoc4C9wPr3f2/g64nVsws28x6h5YzaewQtCHQoqLM3W9092Hunkvj/+lX3f2agMuKOjPrHuoAROgW1zkE+CkNseLu24HNZtY0mvl0oNXOTwkxmjm0PnRSwGVFnZn9GZgG9DezYuAmd78/2KqibipwLbAm9DwG4Cfuviy4kmJiMLA41GM1BXjc3ZOm23WSGQg83fi3GGnAI+7+QrAlxcx3gYdDFxpFwFdba5gw3cxFRCS5JNItPhERSSIKKBERiUsKKBERiUsKKBERiUsKKBERiUsKKBERiUsKKBERiUv/H8Lq8d9rluDnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "B_y = 0\n",
    "Beta_y = 1\n",
    "B_z = 0\n",
    "Beta_z = 1\n",
    "K_xy = 0.1\n",
    "K_xz = 0.1\n",
    "H = 2\n",
    "alpha_y = 1\n",
    "alpha_z = 1\n",
    "\n",
    "time_x_on = 0.5\n",
    "time_x_off = 2.5\n",
    "\n",
    "colors = ['gray', 'gray', 'black']\n",
    "thickness = [0.5,1,2]\n",
    "K_yzs = [0.01, 0.1, .3]\n",
    "fig, axs = plt.subplots( nrows=3, ncols=1 )\n",
    "\n",
    "# X on top\n",
    "axs[0].plot( [sol1.t[0],time_x_on,time_x_on,time_x_off,time_x_off,sol3.t[-1]], \n",
    "            [0,0,1,1,0,0], color='k' )\n",
    "axs[0].set_xlim( [0,6])\n",
    "axs[0].set_ylabel( \"X\")\n",
    "axs[0].set_ylim( [0,1.1])\n",
    "\n",
    "\n",
    "# Z when SY=1\n",
    "for i in range(3):\n",
    "    K_yz = K_yzs[i]\n",
    "    params_off = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,0,1);\n",
    "    params_on = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,1,1);\n",
    "\n",
    "    mf_off = lambda t,y : incoherent_type_1_and_model(t,y, params_off)\n",
    "    mf_on = lambda t,y : incoherent_type_1_and_model(t,y, params_on)\n",
    "\n",
    "    sol1 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(0,time_x_on), y0=(0,0), \n",
    "                                    t_eval=np.arange(0,time_x_on,.01) )\n",
    "    sol2 = scipy.integrate.solve_ivp( fun=mf_on, t_span=(time_x_on,time_x_off), y0=sol1.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_on,time_x_off,.01) )\n",
    "    sol3 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(time_x_off,6), y0=sol2.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_off,6,.01) )\n",
    "\n",
    "    ax = axs[1]\n",
    "    ax.plot(sol1.t, sol1.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation before signal on\n",
    "    ax.plot(sol2.t, sol2.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "    ax.plot(sol3.t, sol3.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "axs[1].set_ylabel( \"Z\")\n",
    "axs[1].set_xlim( [0,6])\n",
    "axs[1].set_ylim( [0,1.1])\n",
    "\n",
    "\n",
    "# Z when SY=0\n",
    "for i in range(3):\n",
    "    K_yz = K_yzs[i]\n",
    "    params_off = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,0,0);\n",
    "    params_on = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,1,0);\n",
    "\n",
    "    mf_off = lambda t,y : incoherent_type_1_and_model(t,y, params_off)\n",
    "    mf_on = lambda t,y : incoherent_type_1_and_model(t,y, params_on)\n",
    "\n",
    "    sol1 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(0,time_x_on), y0=(0,0), \n",
    "                                    t_eval=np.arange(0,time_x_on,.01) )\n",
    "    sol2 = scipy.integrate.solve_ivp( fun=mf_on, t_span=(time_x_on,time_x_off), y0=sol1.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_on,time_x_off,.01) )\n",
    "    sol3 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(time_x_off,6), y0=sol2.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_off,6,.01) )\n",
    "\n",
    "    ax = axs[2]\n",
    "    ax.plot(sol1.t, sol1.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation before signal on\n",
    "    ax.plot(sol2.t, sol2.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "    ax.plot(sol3.t, sol3.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "axs[2].set_ylabel( \"Z\")\n",
    "axs[2].set_xlim( [0,6])\n",
    "axs[2].set_ylim( [0,1.1])\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "                                 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzOUlEQVR4nO3deXxU9fX/8ddJIIGwJhAghGgQqAoWZFdRBKmKShGtVaiiYFuwbmCpLVhr1bYPF/zVpVoLsihu2ELVUHBpCwJa/EIQZAsBCkgCBELYEgLZ5vz+yCRNkJBtJvfemfN8PO5jlntn5lyYk/fMZ+58RlQVY4wxxm0inC7AGGOMORMLKGOMMa5kAWWMMcaVLKCMMca4kgWUMcYYV2rkdAG11bZtW01OTna6DBOm1q5de0hV452uIxCsl4xTatpHnguo5ORkUlNTnS7DhCkR+cbpGgLFesk4paZ9ZEN8xhhjXCloASUic0TkoIhsqmK9iMhLIrJDRDaISJ9g1WKMl1kvmXAVzHdQrwPDz7L+OqCbf5kAvBrEWozxstexXjJhKGifQanqChFJPssmNwLztHSupS9FpLWIJKjq/kDWMXPmTN55551A3qVr/OhHP2LChAlOl2GCzHopuKyP3MvJz6ASgYwKlzP9132LiEwQkVQRSc3Ozq7Vg7zzzjusX7++zkW61fr160Pyj4WpE+ulOrI+cjdPHMWnqjOBmQD9+vWr9ey2F198MZ999lmgy3LUkCFDnC7BeJD1UmXWR+7m5DuovUBShcud/NcZY2rHesmEJCcDKgW4038E0iXAsUCPmRsTJqyXTEgK2hCfiLwLDAHaikgm8FugMYCq/gVYAlwP7ADygfHBqsUYL7NeMuEqmEfxjalmvQL3BevxjQkV1ksmXNlMEsYYY1zJAsoYY4wrWUAZY4xxJQsoY4wxrmQBZYwxxpUsoIwxxriSBZQxxhhXsoAyxhjjShZQxhhjXMkCyhhjjCtZQBljjHElCyhjjDGuZAFljDHGlSygjDHGuJIFlDHGGFeygDLGGONKFlDGGGNcyQLKGGOMK1lAGWOMcaUqA0pEzjnLuiuCU44xocd6yZi6Ods7qM9E5JciEll2hYi0F5G3gOdrcuciMlxE0kVkh4hMPcP6cSKSLSLr/ctPar8LxrhevXrJ+siEq7MFVF+gC7BeRK4SkUnAamAVMKC6O/Y34yvAdUB3YIyIdD/Dpu+p6sX+ZVat98AY96tzL1kfmXDWqKoVqnoEmOhvpn8B+4BLVDWzhvc9ANihqjsBRGQ+cCOwpX4lG+Mt9ewl6yMTts72GVRrEZkBjAeGAwuAj0TkqhredyKQUeFypv+60/1ARDaIyAIRSarhfRvjGfXsJesjE7bONsT3FbAd6Keqn6rqZGAs8HsReTdAj78ISFbVnsA/gTfOtJGITBCRVBFJzc7ODtBDG9Nggt1LNeojsF4y3nK2gBqsqs+panHZFaq6XlUvA5bW4L73AhVfyXXyX1dOVXNUtcB/cRalY/XfoqozVbWfqvaLj4+vwUMb4yr16aWA9ZF/W+sl4xlVBtTZxsdV9bUa3PcaoJuIdBaRKGA0kFJxAxFJqHBxJJBWg/s1xlPq2UvWRyZsVXmQRH2parGI3A98AkQCc1R1s4g8CaSqagrwoIiMBIqBw8C4YNVjjBdZH5lwFrSAAlDVJcCS0657rML5acC0YNZgjNdZH5lwZVMdGWOMcSULKGOMMa5kAWWMMcaVLKCMMca4kgWUMcYYV7KAMsYY40oWUMYYY1zJAsoYY4wrWUAZY4xxJQsoY4wxrmQBZYwxxpUsoIwxxriSBZQxxhhXsoAyxhjjShZQxhhjXMkCyhhjjCtZQBljjHElCyhjjDGuZAFljDHGlSygjDHGuJIFlDHGGFcKakCJyHARSReRHSIy9Qzro0XkPf/6/xOR5GDWY4xXWS+ZcBS0gBKRSOAV4DqgOzBGRLqfttmPgSOq2hV4HngmWPUY41XWSyZcNQrifQ8AdqjqTgARmQ/cCGypsM2NwOP+8wuAl0VEVFWrutP09HSGDBlS4yLWr1/PxRdfXKvCvWL9+vW1+rcwnmW9FETWR+4VzIBKBDIqXM4EBla1jaoWi8gxoA1wqOJGIjIBmOC/WLB8+fJNtSlk+fLliEhtbtLQ2nLaPtfU8uXLA1xKQNV5v1zsfAce03qp5ur0nLM+anA16qNgBlTAqOpMYCaAiKSqaj+HSwqoUNwnCM39EpFUp2uoD+sl7wnVfarJdsE8SGIvkFThcif/dWfcRkQaAa2AnCDWZIwXWS+ZsBTMgFoDdBORziISBYwGUk7bJgW4y3/+FmDp2cbMjQlT1ksmLAVtiM8/Dn4/8AkQCcxR1c0i8iSQqqopwGzgTRHZARymtPGqMzNYNTsoFPcJQnO/GnyfrJdqxfbJG2q0T2IvsowxxriRzSRhjDHGlSygjDHGuJKnAqq66V68RkSSRGSZiGwRkc0iMsnpmgJFRCJFZJ2I/MPpWgJBRFqLyAIR2SoiaSJyqdM11ZX1kXeEWh9B7XrJM59B+ad72QZcTekXFdcAY1R1y1lv6GIikgAkqOpXItICWAuM8vI+lRGRnwP9gJaqOsLpeupLRN4AVqrqLP+RdDGqetThsmrN+shbQq2PoHa95KV3UOXTvahqIVA23Ytnqep+Vf3Kfz4XSKN0RgBPE5FOwA3ALKdrCQQRaQUMpvRIOVS10Ivh5Gd95BGh1kdQ+17yUkCdaboXzz8Jy/hnn+4N/J/DpQTCC8AvAZ/DdQRKZyAbmOsfbpklIs2cLqqOrI+84wVCq4+glr3kpYAKWSLSHFgITFbV407XUx8iMgI4qKprna4lgBoBfYBXVbU3cALw/Gc3ocb6yBNq1UteCqiaTPfiOSLSmNKmeltV/+50PQEwCBgpIrspHT66SkTecrakessEMlW17FX5AkqbzIusj7whFPsIatlLXgqomkz34ilSOi30bCBNVf/odD2BoKrTVLWTqiZT+n+0VFXvcLiselHVLCBDRMpmYB5G5Z+68BLrIw8IxT6C2veSJ2Yzh6qne3G4rPoaBIwFNorIev91j6jqEudKMlV4AHjb/0d9JzDe4XrqxPrIuECNe8kzh5kbY4wJL14a4jPGGBNGLKCMMca4kgWUMcYYV7KAMsYY40oWUMYYY1zJAsoYY4wrWUAZY4xxJQsoY4wxrmQBZYwxxpUsoIwxxriSBZQxxhhX8sxksWXatm2rycnJTpdhwtTatWsPqWq803UEgvWScUpN+8hzAZWcnExqaqrTZZgwJSLfOF1DoFgvGafUtI+CNsQnInNE5KCIbKpivYjISyKyQ0Q2iIhXfwDOmKCyXjLhKpifQb0ODD/L+uuAbv5lAvBqEGsxxstex3rJhKGgBZSqrgAOn2WTG4F5WupLoLWIJASrHmO8ynrJhCsnj+JLBDIqXM70X2eMqR3rJROSPHGYuYhMEJFUEUnNzs52uhxjPMt6yXiJkwG1F0iqcLmT/7pvUdWZqtpPVfvFx4fEEb7GBJL1kglJTgZUCnCn/wikS4BjqrrfwXqM8SrrJROSgvY9KBF5FxgCtBWRTOC3QGMAVf0LsAS4HtgB5APjg1WLMV5mvWTCVdACSlXHVLNegfuC9fjGhArrJROuPHGQhDHGmPBjAWWMMcaVLKCMMca4kgWUMcYYV7KAMsYY40oWUMYYY1zJAsoYY4wrWUAZY4xxJQsoY4wxrmQBZYwxxpUsoIwxxriSBZQxxhhXsoAyxhjjShZQxhhjXMkCyhhjjCtZQBljjHElCyhjjDGuZAFljDHGlSygjDHGuNJZA0pErj3Luh8GvhxjQo/1kTF106ia9UtEZAVwh6ruPW3dNOBvwSnLmJBifRQAqkp2djYHDhwgNzeXU6dO0bRpU2JiYmjfvj3t27dHRJwu0wRQdQG1AXgH+FJEHlLVBRXWVftMEJHhwItAJDBLVZ8+bf04YDpQ1rQvq+qsGtZujFdYH9VBbm4uS5cuZfny5XzxxRds3bqV48ePV7l9dHQ0ycnJ9O7dm/79+zNw4EAGDBhA48aNG7BqE0jVBZSq6msishx4W0RuAO5T1XxAz3ZDEYkEXgGuBjKBNSKSoqpbTtv0PVW9v471G+MF1kc15PP5+OSTT5g7dy6LFi3i1KlTlda3atWKxMREWrRoQZMmTSgoKODEiRPs27ePnJwc0tPTSU9PZ/78+QC0bNmSa665huuuu44bb7yRNm3aOLFbpo6qCygAVHWbiFwK/B5YJyJ31uBmA4AdqroTQETmAzcCpzeWMWHB+qhqJSUlvPXWWzz77LNs2fK/XbvssssYNmwYgwcPplevXrRt27bKYby8vDy2b99Oamoqq1evZuXKlaSnp7NgwQIWLFjAPffcw3XXXcftt9/O97//fZo2bdpQu2fqSlWrXIB1Z7huCLATyK3mtrdQOhxRdnkspUMPFbcZB+yndAhkAZBUxX1NAFKB1HPOOUeNcQqQqmd53p9pcUsfqUt7aenSpXrxxRcrpe8mNTExUZ966inds2dPve97586d+sorr+g111yjERER5Y/RunVrnTRpkm7evFkPHDgQgL0wtVHTPqqusUZVcX0sMLWa29aksdoA0f7zE4Gl1RXct2/fYP2bGVOtOgaU6/pIXdBLx44d0/Hjx5eHRlJSkr7++utaWFgYlMfbv3+/vvDCC9qvX7/yxwQ0Ojpav/jii6A8pjmzmvbRWQ8zV9UPqrj+iJ72Qe0Z7AWSKlzuxP8+xC27nxxVLfBfnAX0reY+jfEc66NvW7lyJT179mTu3Lk0adKE3/3ud6Snp3PXXXcF7aCGDh06MGnSJNasWcNXX33F3XffTUREBAUFBfzgBz+oNLRo3CGYX9RdA3QTkc4iEgWMBlIqbiAiCRUujgTSgliPMV4UUn2kqvz5z39m6NChfPPNN/Tt25d169bx6KOPNuhnQr1792b27NlkZmbSv39/srKyGDJkCBs2bGiwGkz1ghZQqloM3A98QmnD/FVVN4vIkyIy0r/ZgyKyWUS+Bh6kdCzdGOMXSn1UVFTEPffcw3333UdJSQkPP/wwq1at4oILLnCspoSEBJYvX87w4cPJzs5myJAhfP75547VYyqT0uFA7+jXr5+mpqY6XYYJUyKyVlX7OV1HIDRkLxUUFDB69Gg++OADmjRpwqxZs7j99tsb5LFroqCggNtuu40PP/yQ6Oho3njjDW677TanywpZNe0jm4vPGBNU+fn5jBw5kg8++IDWrVuzbNkyV4UTlH7Jd8GCBdx7773lYfr000/jtRfwocYCyhgTNIWFhYwaNYpPP/2U+Ph4PvvsMy655BKnyzqjRo0a8fLLL/Pcc88BMG3aNG699VZyc3Mdrix8WUAZY4LC5/Nx11138c9//pN27dqxYsUKevXq5XRZZyUiTJkyhffff58WLVqwYMEC+vfvz+bNm50uLSxZQBljguKhhx5i/vz5tGjRgo8//tjRgyFqa9SoUaSmpnLRRReRnp7OgAEDmDlzpg35NTALKGNMwM2YMYOXXnqJqKgoPvzwQ3r37u10SbX2ne98hy+//JI77riD/Px8Jk6cyA033MD+/fudLi1sWEAZYwLqiy++4IEHHgBg1qxZDB061OGK6q5Zs2a8+eabvPvuu8TGxvLRRx/Ro0cP5s6di8/nc7q8kGcBZYwJmH379nHLLbdQVFTE5MmTGTt2rNMlBcTo0aPZtGkTw4cP58iRI9x9990MHjzYvtgbZDWazTzcHDlyhGXLlrFy5Uq2bt3Kf//7X44ePUpeXh4+n4+YmBhiYmKIi4sjMTGxfOnSpQsXXHAB559/Pq1atXJ6N4xpUD6fjzvvvJOsrCyGDh3K9OnTnS4poDp27MiSJUt45513mDJlCl988QV9+vThvvvu49FHHyU+Pt7pEkNPTSbsc9MSzAkuP/vsM73ppps0MjKy0mSSdVnat2+vgwcP1okTJ+pLL72kS5cutVmTQwB1mCzWrUuge+m5555TQOPj4zUrKyug9+02R44c0QceeKB8hvTmzZvr448/rsePH3e6NE+oaR/ZTBLAnj17ePDBB/nwww8BiIyM5PLLL2fYsGH07NmTbt260aZNG5o3b46IcPLkSU6cOMGhQ4fYu3cve/fuJTMzk+3bt5Oens62bds4efLkGR8rPj6eHj16cNFFF9GjR4/yJS4uLqD7ZILDZpI4s/Xr1zNgwACKiopYtGgRI0aMOOv2Pp+PI0eOkJubS25uLgUFBfh8Pnw+H40bN6ZJkyY0adKE5s2b07p1a6KjowNSZ6Bt2LCBadOmsWTJEqC0vydPnsy9995L69atnS3OxWraR2EfUB9//DG33347hw8fplmzZvziF79g4sSJJCQkVH/jKvh8PjIzM9m6dStpaWls2rSJzZs3s2nTpiq/9JeQkFAptC666CK6d+9Oy5Yt61yHCTwLqG8rLCykT58+bN68mXvvvZdXXnnlW9uUlJTwzTffsH37djIyMjh48CDNmjWjZcuWtGjRgujoaCIjIxERioqKKCgo4NSpU+Tm5nL06FEaNWpE69atiY+Pp3379uVLs2bNqvwBw4a0YsUKpk6dyqpVqwBo3rw599xzD5MnTyYxMdHh6tzHAqoG3n33XcaOHUtJSQnXXXcds2fPrlcwVUdVyczMrBRYmzdvZsuWLeTn55/xNuecc06l0OrRowcXXHABkZGR9ougDrCA+rbf//73/OY3v6Fr1658/fXXxMTElK87dOgQqampbNy4kdjYWLp160ZycjIdOnSo8bsiVSU/P58jR46QnZ3NgQMHypeIiAg6depEYmIinTp1omPHjo6921JV/vWvf/HMM8/w73//GyidneLGG29k4sSJDBs2jIgIOy4NLKCqlZKSwqhRo1BVpk6dyh/+8AfHnjw+n4/du3dXCq1NmzaxdetWCgoKzngbEeGKK67g0Ucf5aqrriIyMrKBqw5PFlCVpaen07NnTwoLC1m6dGn5IeWHDh1i+fLl7Nq1i969e9OnTx9iY2MDUXY5VeXYsWPlQ+x79+4lKyuL2NhYEhMTSUpKIikpiTZt2jT4u6zU1FSeffZZ/v73v1NSUgJAly5d+PGPf8yYMWNITk5u0HrcxgLqLL7++msuu+wy8vPzeeyxx3jiiScCVF1gFRcXs3Pnzm+940pLS6v0HYyEhARGjx7NbbfdRv/+/e1VWhBZQP2Pz+dj6NChrFixgvHjxzNnzhyKi4tZsWIFa9eu5ZJLLmHAgAEN+o6mpKSEAwcOkJmZSWZmJhkZGRQUFNCpU6fywOrYsSNRUVENUs++ffuYPXs2r732GhkZGeXXX3rppYwePZof/vCHQR21cSsLqCrk5+fTp08f0tPTGTt2LG+88YYrxrBro7CwkH/961+sWLGChQsXsmPHjvJ17du3Z8SIEYwcOZLvfe97lYZbTP1ZQP3P22+/zR133EF8fDxpaWkUFxfzt7/9jfbt2zN8+HBatGgRwGrrLjc3l4yMDDIyMsjMzOTAgQO0bduWpKQkOnToQGxsLOeee25Q/w6UlJTw8ccf89Zbb5GSklI+pC8i9O/fnxtuuIEbbriB3r17h8ULTAuoKjzwwAO8/PLLdO/endTUVM9/jqOqrFmzhrfffpsPP/yQb775pnxddHQ0l112GUOHDmXo0KEMGDCgwV45hioLqFInTpzg/PPPZ+/evcyZM4c+ffrw6aefcu2119KzZ88AVxpYxcXF7Nu3j4yMDLZv3052djaNGjWiS5cudO3alfPOO48mTZoE7fFPnDjBP/7xD959910+/vjjSsP4HTp04Oqrr+bKK69k8ODBdO3a1XMvoGvCAuoM1q1bR9++fYmMjGTVqlV06NCBI0eOkJeXV3rMvQhRUVHExMTQrFmz8qOMGjduHOC9CA5VZePGjaSkpLBo0SJWr15daX1MTAyXXXYZAwcOpH///gwYMCAshxfqwwKq1BNPPMHjjz9Onz59mD59Ops2bWL06NG0a9cuwFUGn6py+PBhduzYwY4dO9izZw9t2rThnHPOKV+aN28elMc+ceIES5cuZfHixSxevJjMzMxK6xMSEhg8eDCDBg2ib9++XHzxxSExKmIBdRpV5corr2TlypWMGjWKgQMHEhsbS9u2bYmJiSEiIgJVpbCwkPz8fE6cOEFeXh65ubk0adKEVq1aVbnExMS48lVOdnY2y5cvZ9myZSxbtoy0tLRvbdOpUyf69+9Pz549y48U7Nq1q2dCuaFZQEFmZibf+c53OHnyJM8++ywtW7bk9ttvD9of8YZW9g5rz5497Nmzh4yMDJo2bUpiYiIdOnQgISGBDh06BDwoVJVNmzaxbNkyVqxYwYoVK8jOzq60TUREBD169KBv37706dOn/Oso7dq1c+XfoKpYQJ1mwYIF/PCHP6R58+a8/vrrXHvttTVqKFUlLy+PY8eOVbkUFRXRokULWrVqRcuWLWnZsmV5eJWdj46OdvwJlJWVxapVq1i9ejWrV68mNTWV48ePf2u7xo0bc/7559O9e3e6du1Kly5d6NKlC+eddx6JiYlhMUZeFQsoGD9+PK+//jqXX34548eP50c/+lFQh8ScpqpkZ2ezb98+9u/fT1ZWFgcOHCA6OpoOHTrQtm3b8qVNmzYBCy5VZevWrSxfvpzVq1ezdu1aNm/eXH5UYEVxcXFceOGFdO/enQsuuIDOnTuXL26cds0CqgKfz8dFF11EWloav/nNb3jyyScDWlNRURHHjh3j+PHj3zotOw+lX95r1qwZzZs3JyYmptLlsiHFZs2aNViY+Xw+tm3bRmpqaqUjBXfv3l3lbaKioujcuTPJyckkJibSsWPHb522a9cuZA97D/eASk9Pp3v37ogIv/3tb3n44YdDOpyqoqocOXKErKwscnJyOHToEIcOHSInJ4fIyEjatGlD69atadWq1bdO6zM6kZ+fz9dff83atWtZt24daWlpbNmypfxvzJnExsaSnJxM586dSUpKIiEhgYSEBDp27Fh+PjY2tkFfQLsioERkOPAiEAnMUtWnT1sfDcwD+gI5wG2quvts91mXpnr//fe5+eabiY+PJyMjo8G/yKeqFBQUkJeXx4kTJ8qHD8tO8/PzK60rKioqn+qladOmlU7Lzh87doxzzjmHmJgYmjZtWn7auHHjej/R8vLy2LJlC2lpaezcuZOdO3fy3//+l507d3LgwIFqbx8ZGUlcXFz5K8ozncbGxpa/2yybTaBly5auHS4t41RAuaWXxowZw/z587n00kv55JNPXHOknluUjbjk5ORw9OhRjh07Vn5a9qI1KiqK5s2bly8VX6SWnZb1eVRUVLX9oKpkZWWxZcsWtmzZwvbt29m1a1f5UtUkABVVfDcYFxdHmzZtzngaFxdHixYtypfmzZvX6cArxwNKRCKBbcDVQCawBhijqlsqbHMv0FNV7xGR0cBNqnrb2e63tk2lqvTs2ZNNmzbxwgsvMGnSpLrsToPy+XycOnWKkydPlp9WPH/q1CkyMjJo2bIl+fn5nDx5svxUVSsFVnR0NFFRUURFRZWfj46OrnQ+KiqKxo0b07hxYxo1alS+NG7cuHz6mTJ5eXns3LmTjIwM9u7dy759+yqd7t27l0OHDtV53yMiIiqFVosWLcr3paZLdHQ0jRs3rrRfZcvp153p8un7XJETAeWWXtq4cSO9evUiIiKCNWvWlP8IYVFREYcPHy4fLThx4kRAfnm27MClsyl7zlZ3evp5p4apVfVbL1DLloovWMv6vLi4mOjo6EovVMue4xWf56cvFZ/bR48eLT/MPisri6ysrPLhyrLlTEP9NRUVFVUptE7v29P7NyYmhl/96leOB9SlwOOqeq3/8jQAVX2qwjaf+LdZJSKNgCwgXs9SlIjUqeBWrVoxf/78SsMRZ2oAVSUiIoLIyEgiIiK+db66yxXPi0itl5rczv/vcMb9LCoqqhRYBQUFFBYWVnladr6oqIji4uLypexySUlJpdCq2Ohl+3r6qc/nq/SusKwByyYGzc3N5cSJE+UHo+Tn55efr2rmjIYWGRlZaSnb35ycHCcCylW9dNNNNzFlyhQKCwuB0ul84uLiyj93bdasWYMEgM/n+9bztarT089X/KJ7QUFB+ahKVX11+j9jTf4u1HSbs/W6z+er1KenTp0q792ioqLy5UyXy64rLCws7+WSkpJKz+eyni4uLiY/P59Tp06Vn5b15ul9XPHFcn5+fl1/uLFGfRTM34NKBDIqXM4EBla1jaoWi8gxoA1Q6SW4iEwAJvgvFgCbalvMsWPHuO6662p7s4bSltP2OUR4dr/KmvkMzm/oWnBZL73//vu8//77tb1ZQ/Hsc+4sQnGfatRHnvjBQlWdCcwEEJHUUPmQukwo7hOE5n6JSGB/66WBWS95T6juU022C+b78L1AUoXLnfzXnXEb/7BEK0o/4DXG/I/1kglLwQyoNUA3EeksIlHAaCDltG1SgLv8528Blp5tzNyYMGW9ZMJS0Ib4/OPg9wOfUHpo7BxV3SwiT1L6c78pwGzgTRHZARymtPGqMzNYNTsoFPcJQnO/GnyfrJdqxfbJG2q0T577oq4xxpjwEL5z1hhjjHE1CyhjjDGu5KmAEpHhIpIuIjtEZKrT9dSXiCSJyDIR2SIim0XE/dNc1JCIRIrIOhH5h9O1BIKItBaRBSKyVUTS/F+e9STrI+8ItT6C2vWSZz6Dqsl0L14jIglAgqp+JSItgLXAKC/vUxkR+TnQD2ipqiOcrqe+ROQNYKWqzvIfSRejqkcdLqvWrI+8JdT6CGrXS156BzUA2KGqO1W1EJgP3OhwTfWiqvtV9Sv/+VwgjdIZATxNRDoBNwCznK4lEESkFTCY0iPlUNVCL4aTn/WRR4RaH0Hte8lLAXWm6V48/yQsIyLJQG/g/xwuJRBeAH4J1GmSLhfqDGQDc/3DLbNEpJnTRdWR9ZF3vEBo9RHUspe8FFAhS0SaAwuByapa92mFXUBERgAHVXWt07UEUCOgD/CqqvYGTgCe/+wm1FgfeUKteslLAVWT6V48R0QaU9pUb6vq352uJwAGASNFZDelw0dXichbzpZUb5lApqqWvSpfQGmTeZH1kTeEYh9BLXvJSwFVk+lePEVK5/afDaSp6h+dricQVHWaqnZS1WRK/4+WquodDpdVL6qaBWSISNkMzMMAr34Ab33kAaHYR1D7XvLEbOZQ9XQvDpdVX4OAscBGEVnvv+4RVV3iXEmmCg8Ab/v/qO8ExjtcT51YHxkXqHEveeYwc2OMMeHFS0N8xhhjwogFlDHGGFeygDLGGONKFlDGGGNcyQLKGGOMK1lAGWOMcSULKGOMMa5kAWWMMcaVLKCMMca4kgWUMcYYV7KAMsYY40qemSy2TNu2bTU5OdnpMkyYWrt27SFVjXe6jkCwXjJOqWkfeS6gkpOTSU1NdboME6ZE5BunawgU6yXjlJr2kQ3xGWOMcaWgBZSIzBGRgyKyqYr1IiIvicgOEdkgIl79hVJjgsp6yYSrYL6Deh0Yfpb11wHd/MsE4NUg1mKMl72O9ZIJQ0ELKFVdARw+yyY3AvO01JdAaxFJCFY9xniV9ZIJV05+BpUIZFS4nOm/7ltEZIKIpIpIanZ2doMUZ4yHWC+ZkOSJgyRUdaaq9lPVfvHxIXGErzGOsF4yXuJkQO0Fkipc7uS/zhhTO9ZLJiQ5GVApwJ3+I5AuAY6p6n4H6zHGq6yXTEgK2hd1ReRdYAjQVkQygd8CjQFU9S/AEuB6YAeQD4wPVi3GeJn1kglXQQsoVR1TzXoF7gvW4xsTKqyXTLjyxEESxhhjwo8FlMeoKqmpqXzve99jxYoVTpdjjDFB47nJYsNRWSi99957/O1vf2PPnj0AZGZmsnXrVoerM8aY4LCAcrE9e/Ywb948Zs+eze7du8uvb9GiBSdPniQ9PZ1PPvmEa6+91rkijTEmSGyIz2WKi4t5//33ueKKKzj33HP5zW9+w+7du4mNjeVnP/sZy5Yt4/Dhw/zud78D4MEHH6SkpMThqo0xJvAsoFwiJyeHp59+mqSkJG6++WY+//xzoqKiuPnmm1m8eDEHDx7kz3/+M0OGDKFRo0ZMnjyZhIQEtm3bxuzZs50u3xhjAs4CymEHDhzg5z//OYmJiUybNo2srCySkpJ4/vnnOXDgAAsXLuT666+nUaPKo7FNmjTh//2//wfAI488Ql5enhPlG2NM0FhAOSQrK4sHHniAc845h+eff56CggKGDBnC4sWL2b17N5MnT6Z169ZnvY/Ro0fTq1cvcnJyeOyxxxqmcGOMaSAWUA0sPz+fxx9/nM6dO/Pyyy9TWFjI8OHDSU1NZdmyZVx//fVERNTsv0VEeO211xARXnzxRTZu3Bjk6o0xpuFYQDUQVWXevHkkJyfzxBNPcOrUKYYPH866dev46KOP6Nu3b53ut3///vz0pz/F5/Nx99134/P5Aly5McY4wwKqAWzfvp3LL7+cu+66i+zsbLp3786yZcv46KOPuPjii+t9/8888wxxcXGkpqbaARPGmJBhARVERUVF/OEPf6BHjx785z//oVWrVsyZM4eNGzcyZMiQgD1O69at+dOf/gTAlClT2LdvX8Du2xhjnGIBFSQ7d+5k4MCBPProoxQVFTFmzBh27NjB+PHja/wZU22MGTOGYcOGkZuby7hx4yidP9QYY7zLAioI3nrrLS666CLWrVtHu3bt+PTTT3nnnXdo27Zt0B5TRHjjjTdo3rw5//znP22ozxjjeRZQAVRQUMC4ceMYO3YsJ0+e5Pvf/z5paWlcffXVDfL4iYmJvPrqq0DpDBPffPNNgzyuMcYEgwVUgOzfv59BgwbxxhtvEBUVxYwZM/jwww+Ji4tr0Dpuv/12RowYwcmTJ7n11lspLi5u0Mc3xphAsYAKgDVr1tCrVy/Wrl1L+/btWbVqFRMmTEBEGrwWEWHOnDnEx8ezevVqHnnkkQavwRhjAsECqp6WLFnC5ZdfTnZ2NgMHDmTDhg306dPH0Zri4+NZsGABIsL06dP56KOPHK3HGGPqwgKqHubNm8fIkSMpLCzk9ttvZ8WKFbRr187psgAYPHgwTzzxBFA6JVLZb0gZY4xXnDWgRKTKHxoSkR8GvhzveOGFF7jrrrsoKSnhF7/4BW+++SZRUVFOl1XJr3/9a4YOHcrx48cZPnw4J06ccLqksGR9ZEzdVPcOaomILBORxDOsm1bdnYvIcBFJF5EdIjL1DOvHiUi2iKz3Lz+paeFOevHFF3nooYcAePbZZ5k+fbojnzdVJyIigoULF3LuueeSlpbGbbfdZlMhOcP6yJg6qC6gNgDvAF+KyC2nrTvrX2QRiQReAa4DugNjRKT7GTZ9T1Uv9i+zali3Y1555RUmT55cfv7hhx92tqBqxMbG8vHHH9OsWTMWL17MtGnV/j00gWd9ZEwdVBdQqqqvAcOAX4nIXBGJKVtXzW0HADtUdaeqFgLzgRvrV66zZs2axf333w+Uvou69957Ha6oZi644AIWLlyIiPDss8/y0ksvOV1SuLE+MqYOanSQhKpuAy4FDgDrRGRgDW6WCGRUuJzpv+50PxCRDSKyQESSznRHIjJBRFJFJDU7O7smJQfcP/7xDyZOnAjAc889x4MPPuhIHXV17bXXMmPGDAAmTZrEW2+95XBF4cfpPgJ39JIxNVVdQJUPP6hqsapOBSYC7wLdAvD4i4BkVe0J/BN440wbqepMVe2nqv3i4+MD8LC1k5qayi233ILP52PatGlMmTKlwWsIhJ/+9Kc89dRTAIwbN47Fixc7XFHYcEUf+R/f0V4ypjaqC6gnTr9CVT8D+gJ/qOa2e4GKr+Q6+a+reF85qlrgvzjLf7+usmvXLq655hoKCgoYM2YMf/hDdbvtblOnTmXKlCmUlJRw0003sWjRIqdLCgdh30fG1MVZA0pVP6ji+iOq+nQ1970G6CYinUUkChgNpFTcQEQSKlwcCaRVW3EDysvLY/jw4Rw5coTBgwfz+uuvu/JovdqaPn06EydOpKioiJtuuomFCxc6XVJIC/c+MqaugvZFXVUtBu4HPqG0Yf6qqptF5EkRGenf7EER2SwiXwMPAuOCVU9tqSp33nkn27Zto3PnzqSkpLjue051JSK8+uqrTJo0iZKSEm699VbefPNNp8syZ+D1PjKmXlTVU0vfvn21ITz99NMKaNOmTTUtLa1BHrOh+Xw+nTp1qlJ6JJk++eST6vP5nC7L1YBUdUEfBGJpqF4y5nQ17SOb6ugMli9fXv59obfffpsLLrjA4YqCQ0R46qmnyr9o/NhjjzFu3DiKioqcLs0YYyygTnf48GFuvfVWVJWHH36Ym266yemSgu4Xv/gFCxYsICoqinnz5nHllVeSlZXldFnGmDBnAVWBqjJu3DgOHjxI7969PX/EXm3cfPPNrFy5kri4OFatWsV3v/tdPv/8c6fLMsaEMQuoCmbNmsWiRYto2rQpCxYsoHHjxk6X1KAGDBjApk2bGDhwIIcOHeLKK6/kmWeesfn7jDGOsIDy2717d/nsEDNmzOC8885zuCJnJCQk8Pnnn/PQQw/h8/mYOnUqgwYNYvfu3U6XZowJMxZQlA7tjR8/nlOnTjFq1CjGjh3rdEmOatSoEX/84x9JSUkhNjaWL7/8ku7du/Paa69RegCOMcYEnwUUMHfuXD777DNatmzJX/7yF6fLcY3vf//7pKenM2LECE6ePMmECRO45JJL2LRpk9OlGWPCQNgHVFZWFpMmTQLgT3/6E+3bt3e4IneJj48nJSWFN998k9jYWFavXk2vXr2YPHkyx48fd7o8Y0wIC/uAuu+++8jLy2PYsGFhP7RXFRHhjjvuYOfOnUyYMAFV5cUXX+Scc87hj3/8IwUFBdXfiTHG1FJYB9Ty5cv5+9//TlRUFLNnzw6JefaCqXXr1syYMYM1a9bQr18/jh07xpQpU0hOTmbOnDkUFhY6XaIxJoSEbUAVFxeX/77TtGnTOPfccx2uyDv69u3L6tWrWbRoEd26dSMrK4sf//jHJCUl8fzzz5OXl+d0icaYEBC2AfWXv/yF9PR0OnbsyK9+9Suny/EcEWHEiBGkpaUxb948unTpwsGDB/n5z39Ox44d+eUvf8nOnTudLtMY42FhGVA5OTk88sgjQOmBEU2bNnW4Iu+KjIxk7NixbNu2jQ8//JC+ffuSm5vL9OnT6dKlC0OGDGHhwoU2v58xptbCMqCefPJJcnNzufzyy8Nirr2GEBERwciRI0lNTeU///kPo0ePJioqiuXLl3PLLbfQtm1b7r77bpYtW0ZJSYnT5RpjPCDsAmrPnj28+uqrALz00kt2YEQQXHrppbz77rvs37+f559/nq5du3L8+HHmzp3LVVddRfv27bnnnnv46KOPOHnypNPlGmNcKuwC6te//jVFRUXcfPPN9O7d2+lyQlpcXByTJ09m27ZtbNiwgalTp5KYmEhOTg4zZszg+uuvJzY2lmuuuYY///nPbNu2zWaqMMaUE6/9QejXr5+mpqbW6bZpaWn06NGDiIgItm7dSteuXQNcnamOqrJ27Vo++OAD3n//fbZs2VJpfVxcHFdccQVXX301Q4YM4cILLyQiwj2vo0Rkrar2c7qOQKhPLxlTHzXto7AKqBEjRrB48WJ+8pOf8NprrwW4MlMX+/btY8mSJaSkpLBy5UqOHj1aaX1MTAzf/e53ueSSSxgwYAD9+vWja9eujoWWBZQx9WcBdZoNGzbQq1cvoqKi2L17NwkJCUGoztSHqrJ161Y+++wz/v3vf7Ny5UoOHjz4re2io6Pp2rUrPXr0oFevXvTo0YMLL7yQzp07B/0nUiygjKm/mvZRo4Yoxg0ef/xxAH7yk59YOLmUiHDhhRdy4YUX8rOf/QwonStx7dq1rFmzhlWrVvHVV19x6NAhNm/ezObNm/nrX/9afvuIiAjatWtHcnIyXbp0oVu3bpx33nl07tyZjh07kpCQYF8pMMZDwuIdVHp6OhdeeCGRkZHs2rWLTp06Bak60xAOHz5MWloaaWlpbNy4ka+//pqtW7dy8ODBag+yaN68OfHx8XTs2JGOHTuSlJREYmIicXFxtGnTptJpXFwcjRpVfg1n76CMqT9XvIMSkeHAi0AkMEtVnz5tfTQwD+gL5AC3qeruQNfxu9/9DlXljjvusHAKAXFxcQwaNIhBgwZVur6goIA9e/awa9cudu3axY4dO9i+fTu7d+9m//79HD58mLy8PPLy8ti1a1eNHismJobWrVvTqlUrYmNjg7E7NeKWXjKmIQUtoEQkEngFuBrIBNaISIqqVjxs68fAEVXtKiKjgWeA2wJZxzfffMP8+fMRER599NFA3rVxmejoaLp160a3bt3OuN7n83H48GH279/Pvn37yk/37dvHoUOHyM7OJicnh8OHD3Ps2DHy8vLIz88nPz+fffv2NfDe/I9besmYhha0IT4RuRR4XFWv9V+eBqCqT1XY5hP/NqtEpBGQBcTrWYoSkToVfOutt/Lee+/V5aYmTPl8Po4dO0ZOTg5Hjx7l6NGjXH311Q0+xOe2XjImABwf4ksEMipczgQGVrWNqhaLyDGgDXCo4kYiMgGY4L9YANT6J13/+te/VvpA3WXacto+h4hQ3K/zHXhMV/WSy4Xicy4U96lGfeSJo/hUdSYwE0BEUkPlQ+oyobhPEJr7JSKePqrAesl7QnWfarJdML/tuBdIqnC5k/+6M27jH5ZoRekHvMaY/7FeMmEpmAG1BugmIp1FJAoYDaSctk0KcJf//C3A0rONmRsTpqyXTFgK2hCffxz8fuATSg+NnaOqm0XkSSBVVVOA2cCbIrIDOExp41VnZrBqdlAo7hOE5n41+D5ZL9WK7ZM31GifPPdFXWOMMeHBPdNEG2OMMRVYQBljjHElTwWUiAwXkXQR2SEiU52up75EJElElonIFhHZLCKTnK4pUEQkUkTWicg/nK4lEESktYgsEJGtIpLm//KsJ1kfeUeo9RHUrpc88xmUf7qXbVSY7gUYc9p0L54iIglAgqp+JSItgLXAKC/vUxkR+TnQD2ipqiOcrqe+ROQNYKWqzvIfSRejqkcdLqvWrI+8JdT6CGrXS156BzUA2KGqO1W1EJgP3OhwTfWiqvtV9Sv/+VwgjdIZATxNRDoBNwCznK4lEESkFTCY0iPlUNVCL4aTn/WRR4RaH0Hte8lLAXWm6V48/yQsIyLJQG/g/xwuJRBeAH4J+ByuI1A6A9nAXP9wyywRaeZ0UXVkfeQdLxBafQS17CUvBVTIEpHmwEJgsqoed7qe+hCREcBBVV3rdC0B1AjoA7yqqr2BE4DnP7sJNdZHnlCrXvJSQNVkuhfPEZHGlDbV26r6d6frCYBBwEgR2U3p8NFVIvKWsyXVWyaQqaplr8oXUNpkXmR95A2h2EdQy17yUkDVZLoXTxERoXQsNk1V/+h0PYGgqtNUtZOqJlP6f7RUVe9wuKx6UdUsIENEymZgHgZ49QN46yMPCMU+gtr3kidmM4eqp3txuKz6GgSMBTaKyHr/dY+o6hLnSjJVeAB42/9HfScw3uF66sT6yLhAjXvJM4eZG2OMCS9eGuIzxhgTRiygjDHGuJIFlDHGGFeygDLGGONKFlDGGGNcyQLKGGOMK1lAGWOMcaX/DxU+SeuPJb8AAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Put it all together. \n",
    "# Figure 3!\n",
    "B_y = 0\n",
    "Beta_y = 1\n",
    "B_z = 0\n",
    "Beta_z = 1\n",
    "K_xy = 0.1\n",
    "K_xz = 0.1\n",
    "H = 2\n",
    "alpha_y = 1\n",
    "alpha_z = 1\n",
    "\n",
    "time_x_on = 0.5\n",
    "time_x_off = 2.5\n",
    "\n",
    "colors = ['gray', 'gray', 'black']\n",
    "thickness = [0.5,1,2]\n",
    "fig, axs = plt.subplots( nrows=3, ncols=2 )\n",
    "\n",
    "# For type 1\n",
    "K_yzs = [0.01, 0.1, .3]\n",
    "\n",
    "# X on top\n",
    "axs[0,0].plot( [sol1.t[0],time_x_on,time_x_on,time_x_off,time_x_off,sol3.t[-1]], \n",
    "            [0,0,1,1,0,0], color='k' )\n",
    "axs[0,0].set_xlim( [0,6])\n",
    "axs[0,0].set_ylabel( \"X\")\n",
    "axs[0,0].set_ylim( [0,1.1])\n",
    "\n",
    "# Z when SY=1\n",
    "for i in range(3):\n",
    "    K_yz = K_yzs[i]\n",
    "    params_off = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,0,1);\n",
    "    params_on = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,1,1);\n",
    "\n",
    "    mf_off = lambda t,y : incoherent_type_1_and_model(t,y, params_off)\n",
    "    mf_on = lambda t,y : incoherent_type_1_and_model(t,y, params_on)\n",
    "\n",
    "    sol1 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(0,time_x_on), y0=(0,0), \n",
    "                                    t_eval=np.arange(0,time_x_on,.01) )\n",
    "    sol2 = scipy.integrate.solve_ivp( fun=mf_on, t_span=(time_x_on,time_x_off), y0=sol1.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_on,time_x_off,.01) )\n",
    "    sol3 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(time_x_off,6), y0=sol2.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_off,6,.01) )\n",
    "\n",
    "    ax = axs[1,0]\n",
    "    ax.plot(sol1.t, sol1.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation before signal on\n",
    "    ax.plot(sol2.t, sol2.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "    ax.plot(sol3.t, sol3.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "axs[1,0].set_ylabel( \"Z\")\n",
    "axs[1,0].set_xlim( [0,6])\n",
    "axs[1,0].set_ylim( [0,1.1])\n",
    "\n",
    "# Z when SY=0\n",
    "for i in range(3):\n",
    "    K_yz = K_yzs[i]\n",
    "    params_off = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,0,0);\n",
    "    params_on = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,1,0);\n",
    "\n",
    "    mf_off = lambda t,y : incoherent_type_1_and_model(t,y, params_off)\n",
    "    mf_on = lambda t,y : incoherent_type_1_and_model(t,y, params_on)\n",
    "\n",
    "    sol1 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(0,time_x_on), y0=(0,0), \n",
    "                                    t_eval=np.arange(0,time_x_on,.01) )\n",
    "    sol2 = scipy.integrate.solve_ivp( fun=mf_on, t_span=(time_x_on,time_x_off), y0=sol1.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_on,time_x_off,.01) )\n",
    "    sol3 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(time_x_off,6), y0=sol2.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_off,6,.01) )\n",
    "\n",
    "    ax = axs[2,0]\n",
    "    ax.plot(sol1.t, sol1.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation before signal on\n",
    "    ax.plot(sol2.t, sol2.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "    ax.plot(sol3.t, sol3.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "axs[2,0].set_ylabel( \"Z\")\n",
    "axs[2,0].set_xlim( [0,6])\n",
    "axs[2,0].set_ylim( [0,1.1])\n",
    "\n",
    "# For type 4\n",
    "K_yzs = [1, 0.3, .1]\n",
    "\n",
    "# X on top\n",
    "axs[0,1].plot( [sol1.t[0],time_x_on,time_x_on,time_x_off,time_x_off,sol3.t[-1]], \n",
    "            [0,0,1,1,0,0], color='k' )\n",
    "axs[0,1].set_xlim( [0,6])\n",
    "axs[0,1].set_ylabel( \"X\")\n",
    "axs[0,1].set_ylim( [0,1.1])\n",
    "\n",
    "\n",
    "# Z when SY=1\n",
    "for i in range(3):\n",
    "    K_yz = K_yzs[i]\n",
    "    params_off = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,0,1);\n",
    "    params_on = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,1,1);\n",
    "\n",
    "    mf_off = lambda t,y : incoherent_type_4_and_model(t,y, params_off)\n",
    "    mf_on = lambda t,y : incoherent_type_4_and_model(t,y, params_on)\n",
    "\n",
    "    sol1 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(0,time_x_on), y0=(0,0), \n",
    "                                    t_eval=np.arange(0,time_x_on,.01) )\n",
    "    sol2 = scipy.integrate.solve_ivp( fun=mf_on, t_span=(time_x_on,time_x_off), y0=sol1.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_on,time_x_off,.01) )\n",
    "    sol3 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(time_x_off,6), y0=sol2.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_off,6,.01) )\n",
    "\n",
    "    ax = axs[1,1]\n",
    "    ax.plot(sol1.t, sol1.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation before signal on\n",
    "    ax.plot(sol2.t, sol2.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "    ax.plot(sol3.t, sol3.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "axs[1,1].set_ylabel( \"Z\")\n",
    "axs[1,1].set_xlim( [0,6])\n",
    "axs[1,1].set_ylim( [0,1.1])\n",
    "\n",
    "\n",
    "# Z when SY=0\n",
    "for i in range(3):\n",
    "    K_yz = K_yzs[i]\n",
    "    params_off = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,0,0);\n",
    "    params_on = (B_y,Beta_y,B_z,Beta_z,K_xy,K_xz,K_yz,H,alpha_y,alpha_z,1,0);\n",
    "\n",
    "    mf_off = lambda t,y : incoherent_type_4_and_model(t,y, params_off)\n",
    "    mf_on = lambda t,y : incoherent_type_4_and_model(t,y, params_on)\n",
    "\n",
    "    sol1 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(0,time_x_on), y0=(0,0), \n",
    "                                    t_eval=np.arange(0,time_x_on,.01) )\n",
    "    sol2 = scipy.integrate.solve_ivp( fun=mf_on, t_span=(time_x_on,time_x_off), y0=sol1.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_on,time_x_off,.01) )\n",
    "    sol3 = scipy.integrate.solve_ivp( fun=mf_off, t_span=(time_x_off,6), y0=sol2.y[:,-1], \n",
    "                                    t_eval=np.arange(time_x_off,6,.01) )\n",
    "\n",
    "    ax = axs[2,1]\n",
    "    ax.plot(sol1.t, sol1.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation before signal on\n",
    "    ax.plot(sol2.t, sol2.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "    ax.plot(sol3.t, sol3.y[1,:].T, color=colors[i], linewidth=thickness[i]); #Simulation after signal on\n",
    "axs[2,1].set_ylabel( \"Z\")\n",
    "axs[2,1].set_xlim( [0,6])\n",
    "axs[2,1].set_ylim( [0,1.1])\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "                                 "
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}