Testing_make_nonlinear_interaction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# In this ipython notebook we will consider our method of constructing nonlinear interactions. We examine the resulting terms and confirm they have the expected behavior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import Potapov_Code.Roots as Roots\n",
    "import Potapov_Code.Potapov as Potapov\n",
    "import Potapov_Code.Time_Delay_Network as Time_Delay_Network\n",
    "import Potapov_Code.Time_Sims as Time_Sims\n",
    "import Potapov_Code.functions as functions\n",
    "import Potapov_Code.tests as tests\n",
    "\n",
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "from scipy.integrate import ode\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def contour_plot(Mat,func = abs):\n",
    "    '''\n",
    "    Make a simple plot to view a matrix\n",
    "\n",
    "    Args:\n",
    "        Mat (compelx-valued matrix): a matrix to view\n",
    "        func (optional[function]): a function to apply to each component\n",
    "\n",
    "    Generates a plot of the matrix.\n",
    "    '''\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "    cax = ax.matshow(func(Mat), interpolation='nearest')\n",
    "    fig.colorbar(cax)\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 400,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Ex = Time_Delay_Network.Example1(r = 0.7, max_linewidth=35.)\n",
    "Ex.run_Potapov()\n",
    "E = Ex.E\n",
    "roots = Ex.roots\n",
    "M1 = Ex.M1\n",
    "delays = Ex.delays\n",
    "modes = functions.spatial_modes(roots,M1,E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 401,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "vecs = Ex.vecs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 402,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "roots.sort(key = lambda z: z.imag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 403,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(-1.1889164797957745-20.943951023931955j),\n",
       " (-1.1889164797957745+0j),\n",
       " (-1.1889164797957745+20.943951023931955j)]"
      ]
     },
     "execution_count": 403,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " roots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 463,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "root_z = lambda z: 1j*z  ## a fake root we will vary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 471,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "roots_to_use = lambda z: [root_z(z).imag,roots[0].imag,roots[-1].imag]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 472,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "modes_to_use = modes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 473,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "delay_indices = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 474,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plus_or_minus_arr =  [1,1,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 488,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.linspace(-10000, 10000,4000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 489,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<module 'Potapov_Code.functions' from 'Potapov_Code/functions.pyc'>"
      ]
     },
     "execution_count": 489,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reload(Potapov_Code.functions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 490,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f = lambda z: functions.make_nonlinear_interaction(roots_to_use(z),modes_to_use,Ex.delays,0,\n",
    "                                                   0,0.01,plus_or_minus_arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 491,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12929b850>]"
      ]
     },
     "execution_count": 491,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEACAYAAACK+7BGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH9JJREFUeJzt3XmUXOV95vHvo6W1IWSJRYolsRhhNtsstmU5jo/bDtaC\nY4R9PLLIidlMQoKVhZkcIzk+QcyQwTrjbYAQkC0fC8eJ0ExiI8/IIGPoHMeOgRgwONoaPMhIIAkQ\nEqANdes3f7y31dXVS1V339slVT2fc+5R1a33vfetq+p6+n3ve28rIjAzMyvKsFo3wMzM6puDxszM\nCuWgMTOzQjlozMysUA4aMzMrlIPGzMwKVVXQSJoraaOkzZJu7KXMbZJaJT0p6YJKdSVNlLRO0iZJ\nD0iakK2fJOkhSa9Luq1sHxdJeirb1tcH9pbNzGwoVQwaScOAO4A5wHnA5ZLOLiszDzgjIs4ErgPu\nqqLuYuDBiDgLeAhYkq0/AHwR+C89NOfvgM9GxNuBt0ua04/3amZmNVBNj2Ym0BoRWyLiELAKmF9W\nZj5wD0BEPAJMkDS5Qt35wMrs8Urgsqz+voj4GXCwdAeSpgDjI+KxbNU9HXXMzOzoVU3QTAWeL3m+\nNVtXTZm+6k6OiB0AEbEdOLmKdmyt0A4zMzvKFDUZQAOo43vhmJnVoRFVlNkGnFLyfFq2rrzM9B7K\nNPVRd7ukyRGxIxsW21lFO3raRzeSHFpmZgMQEQPpKPSpmh7NY8AMSadKagIWAmvKyqwBrgCQNAvY\nnQ2L9VV3DXBV9vhK4L4e9n3kDWfDa3skzZSkbH891eko7yWH5aabbqp5G+plmTIlgJt49tnat6Ve\nFn8+812KUrFHExHtkhYB60jBtCIiNki6Lr0cyyNiraRLJD0D7AWu7qtutullwGpJ1wBbgAUd+5T0\n/4DxQJOk+cDsiNgIfA74NjAaWBsR9+dwDMwK99praZk+HZ5/Ht72tlq3yGzoVDN0RvaFflbZurvL\nni+qtm62fhdwcS91Tu9l/S+Ad1bTZrOjybZtMG0ajBwJW7dWLm9WT6oKGmtczc3NtW5CXXjpJTjp\nJJg6tZltPZ5ZtIHw5/PY4FvQWJ/8g5yPl1+GE0+ECy9s5pVXat2a+uHP57HBQWM2BDp6NBMnwq5d\ntW6N2dBy0JgNgY4ezaRJ8OqrtW6N2dBy0JgNgY4ezaRJ7tFY43HQmA2BXbtSyEyc6B6NNR4HjdkQ\neP11OP5492isMTlozIbAa6/B+PGeDGCNyUFjNgRefz0FzfjxsG8ftLfXukVmQ8dBYzYEOobOhg2D\nsWNh795at8hs6DhozIZAR48G4Ljj4I03atses6HkoDEbAh3naMBBY43HQWNWsIg0VHbccen5+PGp\nh2PWKBw0ZgXbuxdGj4YR2S1s3aOxRuOgMStY6fkZcNBY43HQmBWs9PwMeOjMGo+DxqxgHVObO7hH\nY43GQWNWMA+dWaNz0JgVbO9eGDeu87mHzqzROGjMCnbgAIwZ0/ncPRprNA4as4IdOACjRnU+HzfO\nQWONxUFjVrADB9J1NB3GjoX9+2vXHrOh5qAxK1h50IwZ46CxxuKgMSvYwYPdezT79tWuPWZDzUFj\nVjD3aKzROWjMClY+GcA9Gms0DhqzgrlHY43OQWNWMAeNNToHjVnBPBnAGp2Dxqxg7tFYo3PQmBWs\npws23aOxRuKgMStY+ayzjh5NRO3aZDaUHDRmBSvv0YwcCRIcOlS7NpkNJQeNWcHKgwZ8nsYai4PG\nrGDls87AQWONxUFjVrCeejSeEGCNpKqgkTRX0kZJmyXd2EuZ2yS1SnpS0gWV6kqaKGmdpE2SHpA0\noeS1Jdm2NkiaXbL+cklPZftYK2nSwN622dApnwwA7tFYY6kYNJKGAXcAc4DzgMslnV1WZh5wRkSc\nCVwH3FVF3cXAgxFxFvAQsCSrcy6wADgHmAfcqWQ48HXgQxFxAfA0sGgQ791sSLhHY42umh7NTKA1\nIrZExCFgFTC/rMx84B6AiHgEmCBpcoW684GV2eOVwGXZ40uBVRHRFhHPAa3ZdpS9Pl6SgOOBF/rz\nZs1qwZMBrNFVEzRTgedLnm/N1lVTpq+6kyNiB0BEbAdO7mVb24CpEdEGXE/qyWwl9XhWVNF+s5rq\naTKA/8qmNZIRBW1XlYt00+fla5JGAH8CnB8Rz0m6HfgC8Dc9lV+6dOmRx83NzTQ3Nw+gSWaD11uP\nxkNnVmstLS20tLQUvp9qgmYbcErJ82nZuvIy03so09RH3e2SJkfEDklTgJ0VtnUBENlwGsBqoMeJ\nCdA1aMxqJQLefLP7ZIDRo1MAmdVS+S/hN998cyH7qWbo7DFghqRTJTUBC4E1ZWXWAFcASJoF7M6G\nxfqquwa4Knt8JXBfyfqFkpoknQ7MAB4lhc25kk7Iyn0U2NCfN2s21A4e7LwTQKnRo9NrZo2gYo8m\nItolLQLWkYJpRURskHRdejmWR8RaSZdIegbYC1zdV91s08uA1ZKuAbaQZpoREeslrQbWA4eA6yMi\ngBcl3Qz8RNKbWZ2rcjoOZoXoadgMUg/HPRprFIo6vLOfpKjH92XHnu3b4fzzYceOruv//M/h9NPh\nL/6iNu0y64kkImIg59j75DsDmBWopxln4KEzaywOGrMCeejMzEFjVqjegsazzqyROGjMCtTTfc7A\nQ2fWWBw0ZgXy0JmZg8asUH1NBnDQWKNw0JgVqK9zNB46s0bhoDErkIfOzBw0ZoXqazKAg8YahYPG\nrEAeOjNz0JgVykNnZg4as0J51pmZg8asUB46M3PQmBXKQ2dmDhqzQnnWmZmDxqxQHjozc9CYFaq3\nyQAeOrNG4qAxK5D/TICZg8asUL0FTVMTtLXB4cND3yazoeagMStQb5MBpLTe52msEThozArUW48G\nfJ7GGoeDxqxAfQWNZ55Zo3DQmBWot1ln4AkB1jgcNGYF8tCZmYPGrFAeOjNz0JgVqrdZZ+ChM2sc\nDhqzAnnozMxBY1aoSpMBPHRmjcBBY1agSudo3KOxRuCgMStIe3u6zczIkT2/7qEzaxQOGrOCHDyY\nwkTq+XUPnVmjcNCYFaSvYTPw0Jk1DgeNWUEqBY2HzqxROGjMCtLXjDPw0Jk1DgeNWUGqGTrbv3/o\n2mNWK1UFjaS5kjZK2izpxl7K3CapVdKTki6oVFfSREnrJG2S9ICkCSWvLcm2tUHS7JL1IyXdndVZ\nL+kTA3vbZsWrZujMPRprBBWDRtIw4A5gDnAecLmks8vKzAPOiIgzgeuAu6qouxh4MCLOAh4ClmR1\nzgUWAOcA84A7pSPzdv4K2BERZ0XEucC/DPSNmxWtr9vPgIfOrHFU06OZCbRGxJaIOASsAuaXlZkP\n3AMQEY8AEyRNrlB3PrAye7wSuCx7fCmwKiLaIuI5oDXbDsA1wK0dO42IXdW+UbOh5skAZkk1QTMV\neL7k+dZsXTVl+qo7OSJ2AETEduDkXra1DZhaMrR2i6RfSLpX0klVtN+sJjwZwCwpajJAL5eo9Skq\nvD4CmAb8a0S8G/g58JUB7MdsSLhHY5aMqKLMNuCUkufTsnXlZab3UKapj7rbJU2OiB2SpgA7+9pW\nRLwiaW9EfC9b/79IQ2k9Wrp06ZHHzc3NNDc391bUrBDVzDpzj8ZqqaWlhZaWlsL3U03QPAbMkHQq\n8CKwELi8rMwa4HPAvZJmAbuzAHm5j7prgKuAZcCVwH0l678r6WukYbQZwKPZaz+Q9OGIeBi4GFjf\nW6NLg8asFipNBnCPxmqt/Jfwm2++uZD9VAyaiGiXtAhYRxpqWxERGyRdl16O5RGxVtIlkp4B9gJX\n91U32/QyYLWka4AtpJlmRMR6SatJIXIIuD4iOobVFgPfyULopY79mB2NfAsas6SaHg0RcT9wVtm6\nu8ueL6q2brZ+F6lX0lOdWymZXVay/jfAh6pps1mt+Toas8R3BjArSDWzztyjsUbgoDEriHs0ZomD\nxqwgPkdjljhozApSzawz92isEThozAriHo1Z4qAxK4hvQWOWOGjMCuJb0JglDhqzgvgWNGaJg8as\nIJUmA4wYARHQ1jZ0bTKrBQeNWUEq9WjAvRprDA4as4JUEzQ+T2ONwEFjVpBKs87APRprDA4as4K4\nR2OWOGjMClLtORoHjdU7B41ZQSrNOgPfhsYag4PGrCDu0ZglDhqzglQzGcA9GmsEDhqzgrhHY5Y4\naMwK0HG1/4gKfyzd05utEThozApQzUQA8PRmawwOGrMCVDNsBu7RWGNw0JgVoNqgcY/GGoGDxqwA\n1cw4A/dorDE4aMwK4B6NWScHjVkBfI7GrJODxqwAnnVm1slBY1YA92jMOjlozApQ7WQA92isETho\nzArQnx6Ng8bqnYPGrAD9mXXmoTOrdw4aswJUOxnAPRprBA4aswK4R2PWyUFjVgCfozHr5KAxK4Bv\nQWPWyUFjVgDfgsask4PGrAC+YNOsU1VBI2mupI2SNku6sZcyt0lqlfSkpAsq1ZU0UdI6SZskPSBp\nQslrS7JtbZA0u4d9rZH0VP/eqtnQ8S1ozDpVDBpJw4A7gDnAecDlks4uKzMPOCMizgSuA+6qou5i\n4MGIOAt4CFiS1TkXWACcA8wD7pSkkn19AnhtoG/YbCi4R2PWqZoezUygNSK2RMQhYBUwv6zMfOAe\ngIh4BJggaXKFuvOBldnjlcBl2eNLgVUR0RYRzwGt2XaQNA64Abilv2/UbCj5FjRmnaoJmqnA8yXP\nt2brqinTV93JEbEDICK2Ayf3sq1tJXX+G/BlYH8V7TarGfdozDoVNRlAlYt0E31uUDqfNDy3Jtv+\nQPZhNiQ868ys04gqymwDTil5Pi1bV15meg9lmvqou13S5IjYIWkKsLPCtt4PvFvSr4GRwMmSHoqI\nj/TU6KVLlx553NzcTHNzc9/v0ixH/bkFjXs0VistLS20tLQUvh9F9NmRQNJwYBPwu8CLwKPA5RGx\noaTMJcDnIuJjkmYBX4+IWX3VlbQM2BURy7LZaBMjYnE2GeC7wPtIQ2Y/As6MkoZKOhX4QUS8q5c2\nR6X3ZVak5ma46Sb48Icrlx0+PIXNiGp+7TMrkCQiIvfRooof7Yhol7QIWEcaaluRBcV16eVYHhFr\nJV0i6RlgL3B1X3WzTS8DVku6BthCmmlGRKyXtBpYDxwCrndq2LGm2qEz6LwNzXHHFdsms1qp2KM5\nFrlHY7V24YXwrW+lfyuZNAlaW+GEE4pvl1lfiurR+M4AZgUYSI/GrF45aMwKUO1kAPCfCrD656Ax\nK8CBAzBmTHVl3aOxeuegMStAf4fO3KOxeuagMSvA/v3VB40v2rR656Axy1lE6qFUe47GPRqrdw4a\ns5x1hMywKn+63KOxeuegMctZf87PgHs0Vv8cNGY568/5GXCPxuqfg8YsZ/2Z2gzu0Vj9c9CY5ay/\nQ2fu0Vi9c9CY5ay/Q2e+YNPqnYPGLGf9HTrzLWis3jlozHI2kFln7tFYPXPQmOVsIOdo3KOxeuag\nMcvZ/v39n3XmHo3VMweNWc7cozHrykFjljOfozHrykFjlrOBTG92j8bqmYPGLGcDmd7sHo3VMweN\nWc58U02zrhw0ZjnzLWjMunLQmOVsINOb3aOxeuagMcuZezRmXTlozHLmczRmXTlozHLmWWdmXTlo\nzHLmPxNg1pWDxixnvgWNWVcOGrOc9XfobMyY1Asyq1cOGrOc9XfobOxYB43VNweNWc4GMnT25pvQ\n3l5cm8xqyUFjlrP+Bo3k4TOrbw4as5z19xwNpOGzffuKaY9ZrTlozHK2b1//ejTgoLH65qAxy9m+\nfTBuXP/qOGisnjlozHIUkQJj7Nj+1fPMM6tnVQWNpLmSNkraLOnGXsrcJqlV0pOSLqhUV9JESesk\nbZL0gKQJJa8tyba1QdLsbN0YSf8nW/e0pP8+8LdtVoyDB2HkSBg+vH/13KOxelYxaCQNA+4A5gDn\nAZdLOruszDzgjIg4E7gOuKuKuouBByPiLOAhYElW51xgAXAOMA+4U5KyOv8jIs4BLgR+R9Kcgb5x\nsyIMpDcDDhqrb9X0aGYCrRGxJSIOAauA+WVl5gP3AETEI8AESZMr1J0PrMwerwQuyx5fCqyKiLaI\neA5oBWZGxP6I+JdsH23A48C0/r5hsyLt3dv/8zPgoLH6Vk3QTAWeL3m+NVtXTZm+6k6OiB0AEbEd\nOLmXbW0r35+ktwAfB35cRfvNhox7NGbdjShou6pcpJuoasPScOAfgK9nPZ4eLV269Mjj5uZmmpub\nB9Aks/5x0NixpKWlhZaWlsL3U03QbANOKXk+LVtXXmZ6D2Wa+qi7XdLkiNghaQqws8K2OiwHNkXE\n7X01ujRozIbK3r0OGjt2lP8SfvPNNxeyn2qGzh4DZkg6VVITsBBYU1ZmDXAFgKRZwO5sWKyvumuA\nq7LHVwL3laxfKKlJ0unADODRbNu3AMdHxA39fqdmQ2Ag19CAg8bqW8UeTUS0S1oErCMF04qI2CDp\nuvRyLI+ItZIukfQMsBe4uq+62aaXAaslXQNsIc00IyLWS1oNrAcOAddHREiaCnwB2CDpCdJQ2x0R\n8a28DobZYHnozKy7qs7RRMT9wFll6+4ue76o2rrZ+l3Axb3UuRW4tWzdNnyBqR3lBjN0tnt3/u0x\nOxr4i9ssRx46M+vOQWOWIw+dmXXnoDHLkWedmXXnoDHL0WCGznxTTatXDhqzHHnozKw7B41Zjjx0\nZtadg8YsR4Pp0ezdm397zI4GDhqzHA30HM24cQ4aq18OGrMc7d0LY8b0v9748fDaa/m3x+xo4KAx\ny9Hrr6fQ6K/jj091zeqRg8YsR6+/nkKjv0aNgoj0p6DN6o2DxixHr702sB4NpHru1Vg9ctCY5Wig\nPRpw0Fj9ctCY5SRicD0an6exeuWgMcvJwYMwbFg63zIQ7tFYvXLQmOVkoDPOOniKs9UrB41ZTl57\nbeDnZ8BDZ1a/HDRmOcmjR+OgsXrkoDHLyWB7NB46s3rloDHLiXs0Zj1z0JjlxOdozHrmoDHLyWCu\noYEUNHv25Nces6OFg8YsJ7t2wQknDLz+pElpG2b1xkFjlpNdu1JYDJSDxuqVg8YsJ6+8MrigOeGE\ntA2zeuOgMcvJYIfOHDRWrxw0Zjlxj8asZw4as5wMtkczZky6A/S+ffm1yexo4KAx64edO2Hbtp5f\nG+xkACkFVU8TAnbvhueeG/i2zWrJQWNWpV//Gt75TnjHO+Dxx7u+FjH4oIGeh89eeQXOPz8tP/zh\n4LZvVgsOGrMq3XYbfPazsGwZLFnS9bVdu2DcOBg9enD7OOEEePnlruu++lWYPRu+9z24/npoaxvc\nPsyG2ohaN8DsWLBvH/z938MvfgFTpqSgef55mD49vf7CC/DWtw5+P1Onwtatnc8j4Dvfgfvvh3PP\nTftYswY++cnB78tsqLhHY1aF1ath1iw49dT0FzQ/8Qm4997O1/MKmunTuwbNL38JTU1wzjnp+aJF\n8Ld/O/j9mA0lB41ZFe66C/74jzuff+pTaSirQ15BM21a6il1+MEP4OMfTxMFIPVknnoKnn128Psy\nGyoOGrMKnngiBcm8eZ3rPvIRWL8etm9Pz7duza9H01PQdBg1Cj7zGVixYvD7MhsqVQWNpLmSNkra\nLOnGXsrcJqlV0pOSLqhUV9JESeskbZL0gKQJJa8tyba1QdLskvUXSXoq29bXB/aWzbp67DH4p3/q\n/Y+O3X03/NEfwfDhneuammDuXLjvvvR882Y488zBt2X6dPjNb9LjF16AZ56BD36wa5lrr4Vvf7vn\nSQFtbbB2LaxbB4cPD749ZrmIiD4XUhg9A5wKjASeBM4uKzMP+L/Z4/cBP69UF1gGfD57fCPwpezx\nucATpIkKp2X1lb32CPDe7PFaYE4vbQ7Lx8MPP1zrJhTm0KGIa66JOO20iDlzIqZPj3j66a5lXn01\nYuLEiBde6F7/3nsj5s5Nj9/73oif/rTyPisdzzfeiBgzJuLNNyOWL49YuLDncr/92xHf/37Xdbt2\nRbzvfREzZ0ZcdFHE+98fsWNH5TYdy+r581kL2XdnxVzo71JNj2Ym0BoRWyLiELAKmF9WZj5wT/YN\n/wgwQdLkCnXnAyuzxyuBy7LHlwKrIqItIp4DWoGZkqYA4yPisazcPSV1rCAtLS21bkK/7d8P99wD\nixenE+c7d3Yvc+BAOs/ywgvwH/+RZnV96UtpGvHmzZ3lli+Hj30Mfuu3um9j7lz46U/T1OaNG+Hs\nsyu3rdLxHDcuTTjYsAH++Z/TpIOe/OEfwje+0fl83z74vd+D978ffv5z+Pd/T8N7H/hAzxd6trWl\nyQyLF8Ptt6fjcCw6Fj+fjaia6c1TgZJRY7aSAqRSmakV6k6OiB0AEbFd0skl2/q3kjrbsnVtWf3y\nfdgQikhf5IcPw3HH9V7m1VfTFfQHD8LkyWlK8MiRXcsdPpzKPPNMOrk9ejS8/e1pGm/ptvfsgYcf\nhh/9KF0oGZEumrz00hQMHdeuvPJKOml/++3wnvekWWKPPAJf/CJ8+tPwl38JM2ak61QWLICTT06z\nyZqaUv3f//30hT17NvzkJzBiBHz5y/DjH/f8Po8/Pg1r3XJLOok/2Is1O7z3vWko72c/S+3ryYIF\n8IUvpHbOnJlCc8YM+MpXOicO3HJLOvYf/GC60PMd74BDh+Af/zG9NmUKzJmTjulNN6XHN9yQtgcp\njH7yk1T3kUfSsTnjDLj44u7hG5HC6umnobU1/Z+ccQZceCFMnNi17R0Xt27dmgJ//Hg46aR0DdGw\nHn71bWtLn6cRI9IxLx3C7ElE5zGwo0NR19EM5L858mzAJZdkGy3bakRaDh/ufYlIH/iOZfjw9O/h\nw9Denj747e1pkdLrHWWGD0/r2tq6LocPpx+UjqWjTltb+uHvWNra0nZGjuy6HD6cvrTffDP9e/Bg\naueoUemLsqkpPR4+PP3wHjiQAmH//lS2qSn98JcubW3wxhuwd2/69403UtvHju1c9uxJv+3v2ZN+\n2F99NR3H4cPT+3jrW9MXbMffu9+5M33hjBqVrglpakrrdu5MX+zTp6d7er34YvpNe+LE9AV5xhmp\nzV/9auodTJkCb3tbCoVnn02/qc+eDQsXpjY+/jh87WspHN71rnTMfvWr1AP48Y/hvPM6/8937kzh\nM2tWmia8fn268PLWW7t/aV17bTpX8573pON+ww3pbgC9+YM/SG245ZZBf2SP+Mxn0nv9sz/r/S92\njh0L3/xmmoU2YUJq74oV3b+o//RP0xf4hz6Ugqa1NfW8li+H5ubOcnv2pPoLFqTtTZqUQuO001Ko\n/PVfp97Wxo3p8/D5z6cyp52WPjebN6fPwzvfmX5ZePPN1Kv85S/TLw0zZqTPwiuvpF8sRoxIn48x\nY1L9l15Kx/3EE9P//YQJsGNH+iy98UZ63t6eHo8fn8qdeGL6v123Lt2iZ/fu9D727UvbfctbUr3j\nj0/L4cPpM/r666nMvn3pZ2D8+FRu/Pj0uR01Kh2Tjp/38gU6f+aamtLPYXt7Z/nSf0t/lkeM6Fxf\nWqbjXFrpd4PUvVx7e9pXx3dHx9LR1vJF6vo9NmxY53df6Xeg1PleClNpbA2YBdxf8nwxcGNZmbuA\nT5c83whM7qsusIHUqwGYAmzoafvA/aTzPkfKZOsXAn/XS5vDixcvXrz0fyniHE01PZrHgBmSTgVe\nJH3BX15WZg3wOeBeSbOA3RGxQ9LLfdRdA1xFmhRwJXBfyfrvSvoaaWhsBvBoRISkPZJmZm26Arit\npwZHhDvOZmZHiYpBExHtkhYB60izyFZExAZJ16WXY3lErJV0iaRngL3A1X3VzTa9DFgt6RpgC7Ag\nq7Ne0mpgPXAIuD6bDQEpzL4NjAbWRsT9ORwDMzMrkDq/w83MzPJ3TNwZQNKnJP1KUruki8pe69fF\nnZKaJK3K6vybpFNKXrsyK79J0hVD8+5qS9JNkrZKejxb5pa8ltuxteoufDaQ9JykX0p6QtKj2Tpf\n4F0lSSsk7ZD0VMm63I7fgH7Oizjxk/cCnAWcCTwEXFSy/hz6eXEn8CfAndnjT5Ou2QGYCDwLTADe\n0vG41u99CI7tTcB/7mF9bsfWS0AVFz57OXKsfg1MLFtX2AXe9bYAvwNcADxVxPEbyM/5MdGjiYhN\nEdFK92nT8+n/xZ2lF4r+b+Aj2eM5wLqI2BMRu0nnlY78dl/nepo8kcex/d3imnzMqebCZ0tE99EW\nX+BdpYj4V+DVstV5Hr9+/5wfE0HTh/ILQjsu7pxK7xd3HqkTEe3AHkmT+thWI1iU3aPumyVd6jyO\n7e7s2FrvFzVbdwH8SNJjkq7N1nW5wBsovcC7v5/TRnRyjsev3z/nR80fPpP0I9K1N0dWkT5wfxUR\nPyhy1wVu+6jQ17EF7gT+a0SEpFuArwDXdt/KwHad03assXwgIl6UdBKwTtIm0ue1lGcxDU6ex6/i\nz/lREzQR8dEBVNsGTC95Pi1b19v60jovSBoOHB8RuyRtA5rL6jw8gDYddfpxbL8BdIR6bsd2oO2u\nM9uA0pOmpcfNSkTEi9m/L0n6PmnYcYekyZGuz5sCdNzBbiCf00aU5/Hr98/5sTh0Vpqea4CF2SyI\n0+m8uHM7aUhspiSRLu4svSD0yuzxfyJNMAB4APiopAmSJgIfzdbVtexD1+GTwK+yx3keWyu58FlS\nE+ni5TU1btNRR9JYScdlj8cBs4Gn6bzAG7pf4N3fz2kjEN2/K6/KHg/2+PX/57zWMySqnEVxGWlM\ncD/pDgM/LHltCWmmxAZgdsn6d5M+oK3A/yxZPwpYna3/OXBayWtXZes3A1fU+n0P0bG9B3iKNAvq\n+2S3Bcr72HoJSJNLNmXHZ3Gt23M0LsDp2WfxiewztjhbPwl4MDt+64C3lNTp1+e03hfgH4AXgIPA\nb0gX0E/M6/gN5OfcF2yamVmhjsWhMzMzO4Y4aMzMrFAOGjMzK5SDxszMCuWgMTOzQjlozMysUA4a\nMzMrlIPGzMwK9f8BgIlb89H0tr8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x122275150>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, [abs(f(z))**2 for z in x])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 492,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Ex = Time_Delay_Network.Example3(r1 = 0.7, r3 = 0.7, max_linewidth=35.)\n",
    "Ex.run_Potapov()\n",
    "E = Ex.E\n",
    "roots = Ex.roots\n",
    "M1 = Ex.M1\n",
    "delays = Ex.delays\n",
    "modes = functions.spatial_modes(roots,M1,E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 493,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vecs = Ex.vecs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 494,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "roots.sort(key = lambda z: z.imag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 495,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(-1.1867857239576904-53.731290378625815j),\n",
       " (-1.1714495298385854-40.60438955513678j),\n",
       " (-1.2178438327806835-32.542503045139362j),\n",
       " (-1.1500606861696614-20.178105255320496j),\n",
       " (-1.236142384757039-10.901754569501422j),\n",
       " (-1.1431816829945924+1.2586240840022705e-167j),\n",
       " (-1.2361423847570399+10.901754569501421j),\n",
       " (-1.1500606861696616+20.178105255320492j),\n",
       " (-1.2178438327806838+32.542503045139362j),\n",
       " (-1.1714495298385852+40.60438955513678j),\n",
       " (-1.1867857239576902+53.731290378625815j)]"
      ]
     },
     "execution_count": 495,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " roots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 496,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-2.358235253796275+94.33567993376259j)"
      ]
     },
     "execution_count": 496,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "roots[-1]+roots[-2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 497,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "root_z = lambda z: 1j*z  ## a fake root we will vary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 563,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "roots_to_use = lambda z: [root_z(z).imag,roots[0].imag,roots[1].imag]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 564,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "modes_to_use = [modes[0],modes[0], modes[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 565,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "delay_indices = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 566,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plus_or_minus_arr =  [1,1,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 582,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.linspace(-800,1000,4000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 583,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<module 'Potapov_Code.functions' from 'Potapov_Code/functions.pyc'>"
      ]
     },
     "execution_count": 583,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reload(Potapov_Code.functions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 584,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = lambda z: functions.make_nonlinear_interaction(roots_to_use(z),modes_to_use,Ex.delays,0,\n",
    "                                                   0,0.01,plus_or_minus_arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 585,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### by default, indices_of_refraction=[1,1,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will vary the fake root we introduced to obtain the phase-mismatch diagram. That is, the phase mismatch $\\delta k$ is going to be some linear function of $z$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 586,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12a91b210>]"
      ]
     },
     "execution_count": 586,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEGCAYAAACevtWaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYVMW9xvHvj01BFMEoIhBxQ8UNMSJuYRRXRLgiIKBI\niLhEjcYtgiaKxhuUuJC4oyyCAooKKKIiyniFALIKyCCbKBgkKosLoDBT94/q0XEclp453dV9+v08\nTz/T9Jzp89LT8zunq+pUmXMOERGJl0qhA4iISPRU3EVEYkjFXUQkhlTcRURiSMVdRCSGVNxFRGIo\nZcXdzAaa2RozmxfBc+WZ2Rwzm534usnM2kaRU0QkjixV49zN7BTgW2Coc+7oCJ+3NrAEaOCc2xzV\n84qIxEnKztydc5OBdSUfM7MDzex1M5thZu+aWeNyPHUH4HUVdhGRbUt3m/sA4Frn3PHALcDj5XiO\nzsCISFOJiMRMlXTtyMx2A04CRpmZJR6umvjeBcDdQMk2IgNWOefOLfEc+wJHAm+mJbSISJZKW3HH\nf0pY55xrVvobzrnRwOideI5OwGjnXGHU4URE4iTpZhkzq2Vmo8yswMw+NLMTtrd54oZz7hvgYzPr\nUOK5ku1o7YKaZEREdqg8be7/BMY75w4HjgEKytrIzIYD/wYam9mnZtYDuBi4zMzmmtkCYKeHM5rZ\n/vgRMu+WI7OISE5Jaiikme0BzHHOHZS6SCIiUlHJnrkfAHxpZoMTFxQNMLPqqQgmIiLll2xxrwI0\nAx5NdIxuBHpFnkpERCok2dEyq4CVzrmZiX+/CNxaeiMz0/JOIiLl4JyzHW+1Y0mduTvn1gArS1xZ\n2gpYuI1tM+p25513Bs+QDZkyNZcyKVMu5IpSeca5Xwc8Z2ZVgeVAj0gTiYhIhSVd3J1zHwDHpyCL\niIhEJGfmc8/Lywsd4RcyMRNkZi5l2jnKtPMyNVdUUjLlr5m5VDyviEicmRkuRIeqiIhkBxV3EZEY\nUnEXEYkhFXcRkRhScRcRiSEVdxGRGFJxFxGJIRV3EZEYUnEXEYkhFXcRkRhScRcRiSEVdxGRGFJx\nFxGJIRV3EZEYUnEXEYkhFXcRkRgqzxqqIlnBOVi6FN5/H1auhE2boE4dOOwwaNECatUKnVAkdXTm\nLrHz3XfwwANwxBFw2mkwejSsXQtmvtj37QsNG8IFF8A774ROK5IaOnOX2HAOhgyB226DU0+FAQPg\n5JN9US/t66/hhRfgyivh17+Gxx+Hxo3THlkkZbSGqsTCF19Ajx7w2WcwcCA0a7ZzP7dlCzz2GPzt\nb3D33fCHP5R9MBBJhyjXUFVxl6xXUABt2kCHDr5IV6uW/HMsWQIXXujb4h99FKpWjT6nyI6ouIsk\nTJ8ObdtCv37QvXvFnuubb6BrVygqgpdegl13jSajyM5ScRcBZs6E1q1h8GA477xonnPLFrj4Yt8p\nqwIv6abiLjlv8WL47W/hySehXbton3vLFujc2TfNDB8OlTSmTNIkyuKut61knbVrfRv7PfdEX9jB\nF/XnnoNPP4U77oj++UXSQcVdssrWrdCxI5x/PvTsmbr97LorjB3rz9xHjkzdfkRSJelmGTNbAWwA\nioAtzrnmZWyjZhlJiTvugKlT4Y03oHLl1O9vzhw46yz497/hkENSvz/JbUHb3M1sOXCcc27ddrZR\ncZfI5edDly6+4O67b/r2+/jj/oKoqVPVwSqpFbrN3cr5cyLltnYtdOsGgwalt7ADXHUVHHSQ2t8l\nu5T3zH09UAgMcM49VcY2OnOXSPXoATVrwsMPh9n/F1/AUUf5dvgTTgiTQeIvyjP38swtc7JzbrWZ\n7Q28ZWYFzrnJUYQRKcvEiX6CrwULwmXYe2/45z/9QWb2bDXPSOZLurg751Ynvn5hZqOB5sAvinuf\nPn1+vJ+Xl0deXl65Q0ru2rjRT+71xBOw++5hs3TqBM8/D//7v36aA5GKys/PJz8/PyXPnVSzjJnV\nACo55741s92ACcBdzrkJpbZTs4xEolcvP958+PDQSbxVq+CYY2DGDDjwwNBpJG6CjZYxswOA0YDD\nn/U/55y7t4ztVNylwpYt8+3b8+dDvXqh0/ykb1+YNs23v4tESdMPSE5o3x6OPx569w6d5Oe+/x6O\nPBL+9S8499zQaSROVNwl9iZNgt//3k/nm4mdl6+9BjfeCB9+CFW05I1EJPQ4d5GUKiyEG27w0/hm\nYmEHPxtl/fp+RkqRTKTiLhnn+eehRg2/+EamMoN774W77vILb4tkGhV3yShbt8Kdd/oZHzN9ubvm\nzX2Hb6gLq0S2R23uklEGDYJnn/UXLWWDRYv8YtyLF0Pt2qHTSLZTh6rE0vffw6GH+rnUTz45dJqd\nd9llvv397rtDJ5Fsp+IusfTYYzBuHIwfHzpJcorH4y9bBrVqhU4j2UzFXWLnhx/8zIsvv+zHtmeb\nbt3g8MPhtttCJ5FspuIusTN4sJ9i4K23Qicpn4ICyMvzZ+81a4ZOI9lK49wlVoqK4L77/Dwy2erw\nw6FlS79gt0gmUHGX4MaOhT32gNNPD52kYv7yF3jgAdi8OXQSERV3Ccw5PxFXr16ZP659R44+Gpo2\n9aN9REJTcZegJk2Cr7+G//mf0EmicdNN8OCD/qAlEpKKuwTVrx/8+c9QKSbvxNNP9xOJvflm6CSS\n6zRaRoJZtMh3Qn7ySeZOEFYeQ4fCsGHZO/JHwtFoGYmFhx+GK66IV2EH6NwZFi6EDz4InURymc7c\nJYj16/0ydQsWwH77hU4Tvb594aOPYMiQ0Ekkm+giJsl6Dz4Is2bFd2TJ2rX+ittFi6Bu3dBpJFuo\nWUayWmEhPPIIXHdd6CSpU6cOXHghDBwYOonkKhV3Sbtx42CfffxkW3F29dXwxBP+YCaSbiruknb/\n+hdcf33oFKnXrJnvT3jttdBJJBepuEtaLVrkF5W+8MLQSdLjmmvg0UdDp5BcpOIuafXUU9CjB1Sr\nFjpJenTsCHPmwJIloZNIrtFoGUmbzZuhYUOYPt0Pg8wVt97q14Z94IHQSSTTabSMZKWXX4Zjj82t\nwg5w1VXwzDOwcWPoJJJLVNwlbZ58Eq68MnSK9DvgAGjeHEaNCp1EcomKu6TFokWweDG0bRs6SRg9\ne2rMu6SXirukxYABviO1atXQScJo08ZPR/DRR6GTSK5Qh6qkXHFH6vvv+yaKXHXLLX5q4/vuC51E\nMlXwDlUzq2Rms83slShCSLy99JK/oCeXCzvAZZf5jtUtW0InkVxQ3maZ64GFUQaR+Bo4EC6/PHSK\n8A47DA4+WFesSnokXdzNrAHQGng6+jgSNytWwPz5cP75oZNkhp494Wn95UgalOfM/SHgFkCN6rJD\nQ4f6xSt22SV0kszQsSNMmQKffRY6icRdlWQ2NrPzgDXOublmlgdss+G/T58+P97Py8sjLy+vfAkl\naxUV+cUqNL77J7vtBp06+dfl9ttDp5HQ8vPzyc/PT8lzJzVaxsz+DlwCbAWqA7sDLzvnLi21nUbL\nCO++C9deC/PmgUXS/x8P778PXbrA0qV6XeTngo2Wcc7d5pz7tXPuQKAz8E7pwi5SbMgQ+N3vVMBK\nO/54v27slCmhk0ic6SImSYlvv4XRo+Hii0MnyTxmcOmlvj9CJFV0EZOkxDPPwIsvwquvhk6SmVat\ngqOP9h2r1auHTiOZIvhFTCI7UtwkI2Vr0ACOO04HP0kdFXeJ3Mcfw4IFGtu+I2qakVRScZfIDR3q\nR4PkympL5XXBBTB5MqxZEzqJxJGKu0TKOd/e3r176CSZr2ZNaNcORowInUTiSMVdIjV1qh/m16xZ\n6CTZQU0zkioq7hKpESOga1eNbd9ZeXnwxRd+/h2RKKm4S2S2boUXXvBzycjOqVwZLrkEhg0LnUTi\nRsVdIvPOO9CokZ/WVnZet27w7LNQWBg6icSJirtEZvhwP0pGktOkCdSvD2+/HTqJxImKu0Ri0yYY\nOxYuuih0kux08cUaNSPRUnGXSIwf70fI1KsXOkl2uugiGDPGHyRFoqDiLpEoHiUj5VOvnp+OYPz4\n0EkkLlTcpcI2bIC33oL27UMnyW5du/p+C5EoqLhLhY0ZA6edBrVrh06S3dq3h4kTYf360EkkDlTc\npcI0SiYae+4JrVr5efBFKkrFXSpkzRqYPl0zQEalSxc1zUg0VNylQkaNgjZtoEaN0EnioU0bmDkT\nVq8OnUSynYq7VIhGyUSrenVo29ZP4yBSESruUm4rVsDixXDmmaGTxItGzUgUVNyl3EaOhA4doGrV\n0EnipVUrf+BcujR0EslmKu5SbholkxpVqkCnTpqOQCpGxV3KZcECWLcOTjkldJJ46toVnnvOr2wl\nUh4q7lIuI0b4edsr6R2UEi1awPffw9y5oZNIttKfpiTNOY2SSTUzdaxKxai4S9KmT4dq1aBp09BJ\n4q1rV38QLSoKnUSykYq7JE3rpKbHEUfAXnvBe++FTiLZSMVdkrJ1Kzz/vEbJpEuXLho1I+Wj4i5J\nyc+Hhg3hkENCJ8kNnTvDSy/BDz+ETiLZJqnibma7mNl0M5tjZvPN7M5UBZPMpLHt6dWoETRu7OfL\nF0mGuSQH0ppZDefcRjOrDEwBrnPOvV9qG5fs80rm27wZ9tsP5s/3CzpLejzyCEybBs8+GzqJpJqZ\n4ZyLpDcr6WYZ59zGxN1dgCqAqniOeP11OOYYFfZ069gRxo2DjRt3vK1IsaSLu5lVMrM5wOfAW865\nGdHHkkykse1h1K0LJ5wAr74aOolkkyrJ/oBzrgg41sz2AMaYWRPn3MLS2/Xp0+fH+3l5eeTl5VUg\npoT29dfw5pvwxBOhk+Sm4jHvF10UOolEKT8/n/z8/JQ8d9Jt7j/7YbO/At855x4s9bja3GNm2DC/\nMMcrr4ROkpu+/tqPUlqxQmvVxlmwNncz+5WZ1Urcrw6cCSyKIohkNo2SCWuPPeCMM+Dll0MnkWyR\nbJt7PWCSmc0FpgNvOufGRx9LMskXX8DUqX6FIAlHc81IMirULLPNJ1WzTKw89hhMnqzCEtqmTX4o\n6sKFUK9e6DSSCkGHQkru0SiZzFC9OrRr56d/ENkRFXfZrk8+gYICOOus0EkENNeM7DwVd9mukSPh\nwgv9FL8SXvH6qsuWhU4imU7FXbZrxAiNkskkVar4K1Z19i47ouIu27RwIXz5JZx6augkUlKXLr5z\nW2MWZHtU3GWbiq+IrFw5dBIp6cQT/Twz8+aFTiKZTMVdyuScPzvUKJnMU6mSn+ddTTOyPSruUqYZ\nM/wZe7NmoZNIWbS+quyIiruUqXi6Aa2TmpmOOgpq1vRXDouURcVdfqGw0F8ooyaZzGX209m7SFlU\n3OUXJk3yC3IcemjoJLI9nTvDCy/4RctFSlNxl19QR2p2OOggOOAAePvt0EkkE6m4y89s3gxjxviz\nQsl8milStkXFXX7mtdfg2GP97IOS+Tp18guobNoUOolkGhV3+Rk1yWSXevXguONgvFZVkFJU3OVH\n69fDxIl+ojDJHsXTEYiUpOIuPxo92s86uOeeoZNIMtq39wflDRtCJ5FMouIuP3ruOTXJZKPateG0\n0/zBWaSYirsAsHo1zJoF550XOomUhxbxkNJU3AXwV6S2a+eXcpPsc/75MH06rFkTOolkChV3ATRK\nJtvVqAFt2sCoUaGTSKZQcReWLIFPP4XTTw+dRCpCc81ISSruwvDhflGOKlVCJ5GKOPNMWLzYr7Eq\nouKe47QoR3xUrQodOvhFzUVU3HPcrFl+VsHmzUMnkSjogiYppuKe44YNg27dtChHXJxyCqxbBwsW\nhE4ioam457AtW3wHXLduoZNIVLS+qhRTcc9hb7wBjRv7ecElPoqnAdb6qrktqeJuZg3M7B0z+9DM\n5pvZdakKJqk3dChcemnoFBK1pk39+qrvvRc6iYRkzrmd39hsX2Bf59xcM6sJzALaOecWldrOJfO8\nkn7r1vlVfFas0ERhcXT//bBwIQwaFDqJJMPMcM5F0gOW1Jm7c+5z59zcxP1vgQKgfhRBJL2efx7O\nPluFPa4uvthPJPbdd6GTSCjlbnM3s0ZAU2B6VGEkfdQkE2/16sFJJ2mmyFxWrmsSE00yLwLXJ87g\nf6FPnz4/3s/LyyMvL688u5IUWLIEli2Ds84KnURSqXt3GDAALrkkdBLZlvz8fPLz81Py3Em1uQOY\nWRVgHPC6c+6f29hGbe4Z7I474Jtv4KGHQieRVNq8GerXh7lzoWHD0GlkZwRrc08YBCzcVmGXzFZU\n5C9cUpNM/O26K3Ts6H/fknuSHQp5MnAxcLqZzTGz2WZ2TmqiSSpMnuyHyTVtGjqJpMPvfgdDhvg5\nhCS3JNXm7pybAlROURZJg+KOVE03kBtOOMH/rqdNgxNPDJ1G0inpNvedelK1uWekb7/1ba8LF/rR\nFJIb/v53P1//E0+ETiI7ErrNXbLUCy/AqaeqsOeabt38Ck2bN4dOIumk4p5DBg6Enj1Dp5B0a9gQ\nmjWDsWNDJ5F0UnHPEQUF8PHH0Lp16CQSQo8e/uAuuUPFPUcMHOg7UrWUXm5q3x7mzPEHeMkN6lDN\nAT/84D+aT54MhxwSOo2E8qc/+WGw99wTOolsizpUJSnjxsFhh6mw57rLL4fBg/2yihJ/Ku454Omn\n4bLLQqeQ0I44AvbfH157LXQSSQcV95hbtcpfwNKhQ+gkkgmuuAKeeip0CkkHFfeYGzIELroIatQI\nnUQyQceO8O9/+4O+xJuKe4wVFvqztMsvD51EMsVuu/kFtLVCU/ypuMfY+PGw777+AhaRYpdf7vth\nCgtDJ5FUUnGPsccfh6uvDp1CMs2xx0LdujBhQugkkkoq7jG1bBnMmAGdOoVOIpnoiiv8wV/iSxcx\nxdSf/+wX5rj//tBJJBN9950fFjlzJjRqFDqNFIvyIiYV9xjavBl+/WuYMkUXLsm23XSTn47ivvtC\nJ5FiukJVtmvUKN+uqsIu2/OHP/hRM5s2hU4iqaDiHkOPP+7/cEW25+CDoXlzGDkydBJJBRX3mJkz\nB1auhDZtQieRbPDHP8LDD2uN1ThScY+Z/v3hmms0ta/snLPOgm++galTQyeRqKlDNUZWr4YmTfww\nyDp1QqeRbNG/P7z/PgwfHjqJaLSMlOmvf4WvvoLHHgudRLLJ+vVwwAHw4Yew336h0+Q2FXf5hU2b\n/Ljl996DQw8NnUayzR//CLvvDn//e+gkuU3FXX7hqaf8AsjjxoVOItlo+XI/cubjj32RlzA0zl1+\nxjnfbnrDDaGTSLY68EA4/XTNFhknKu4xMGECVK7s/zhFyuvmm+Ghh7QMX1youMdAv37+UnKL5MOc\n5KrmzX2/zahRoZNIFFTcs9y0abB0KXTtGjqJxMEtt/jJ5tRllv2SLu5mNtDM1pjZvFQEkuT07ev/\nIKtWDZ1E4qB1a9i4ESZNCp1EKqo8Z+6DgbOjDiLJW7AApk+Hyy4LnUTiolIlP120hkRmv6SLu3Nu\nMrAuBVkkSffeC3/6E1SvHjqJxMkll/irnKdMCZ1EKkJt7llq+XJ4/XXN/ijRq1oVbrsN/va30Emk\nIlTcs9S998JVV0GtWqGTSBx17w4FBb7ZT7JTyuYO7NOnz4/38/LyyMvLS9Wucs6yZfDSS7B4cegk\nElfVqkGvXnD33fDaa6HTxFd+fj75+fkpee5yTT9gZo2AV51zR23j+5p+IIW6d/cTPZU4fopE7vvv\n/YIeo0fDb34TOk1uCDq3jJkNB/KAvYA1wJ3OucGltlFxT5FFi+DUU/3YdjXJSKo98ojv29HZe3po\n4rAc1rkzHHMM9O4dOonkgu+/h8MPh8GDoWXL0GniT8U9R82b51fOWboUatYMnUZyxXPP+aX4pk7V\nFBepplkhc1SvXv6mwi7p1KULbN4MY8aETiLJUHHPEhMmwJIlcPXVoZNIrqlUyQ+97d1bM0ZmExX3\nLFBY6Kdj7dfPD1ETSbezz4Z69Xzbu2QHtblngaefhqFD4d131eYp4cye7ScWW7QI9twzdJp4Uodq\nDvnmG78m6iuvaKyxhHfllX4uo/79QyeJJxX3HHLzzfDllzBkSOgkIv692KQJvPMOHHlk6DTxo+Ke\nI+bNgzPO8FP77rNP6DQi3qOPwssvw8SJaiaMmoZC5oCiIj8x2D33qLBLZrnySn8G/8ILoZPI9qi4\nZ6hBg/xSZz17hk4i8nNVqsCTT/q1BL76KnQa2RY1y2Sg1auhaVN4803/VSQT3XijP4MfOjR0kvhQ\nm3uMOQdt2/qirsUSJJN99x0cfbSfmqB169Bp4iHK4p6y+dylfIYMgZUr/XztIplst91gwADo0cN3\n/mvse2bRmXsG+fRTOO44ePttf0Ykkg2uvdY3z4wYodEzFaXRMjG0datfmPjGG1XYJbv84x/w4Ye6\nFiPT6Mw9Q/TuDbNm+YURKlcOnUYkOQsWwGmnwZQp0Lhx6DTZS2fuMTN+PAwbBs8+q8Iu2enII/16\nqx06+I5WCU9n7oEtXw4nnQSjRvnl80SylXPw+9/74v7882p/Lw+ducfE+vXQpg389a8q7JL9zODx\nx+GTT6Bv39BpRGfugWzd6scGH3qoHycsEheffQYnnAAPPQQdO4ZOk100zj3LFRXBFVf4y7gfeih0\nGpFo1a8P48b59X7r1IFWrUInyk1qlkkz5+D662HxYt/OXkWHV4mhpk39xGKdO8PMmaHT5KaUNct8\n9ZWjRg3YZRd1rBRzDm691c+F/fbbUKtW6EQiqTV2rP+U+sorvqlGvKIi3/H8zTd+GU3n/K1Royxo\nljnoINi4EbZsgdq1/fqL++7rv+6/vx8L27ixb3OuXTtVKTJHYSFcc40/i3njDRV2yQ3t2vlPp+ef\nDy++CL/9behEqVdU5K82X7jQL0m4ciWsWuVvn30GGzbAt9/6Fa1q1vSvj1n0J8Ep71DduhXWrfMz\nHX7+uf+6YoVvlli8GD76yM9Rcdxxfhm54q/16kUeK5hNm+DSS2HtWhgzBnbfPXQikfR6+23fRPPo\no9CpU+g00Skq8nVs2jR/mzULCgr8PDtNmsBhh/mT2QYNoGFD2G8/fzJbs2bZ17TEalZI5/xRbtYs\nf5s5E2bMgF/9Clq2hLw8/7VBg8hjpsUnn8CFF8Ihh/jLs3fZJXQikTDmzvVn8j16wB13QKUs7PEr\nLPT/j4kTYdIkmD7ddxq3aOGbnX7zG1/UyzuJWqyKe1mKivzlzPn5/vZ//+ePdmefDeee6y9zrlEj\nsrgp8/rr/o18yy1+zhj1PUiuW7MG2rf3zZKDBvmm2ky3fLkv5m+95fvL6tb1y1+2agUnnhjtSmmx\nL+6lFRX5KUXfeMMXzNmz/VWd557rb40bZ1bh/PpruOkm/2YYMsR/+hARb8sWuOsuePppeOwxuOCC\nzPr73bwZ3n3XTwsyfrxvHz/jjJ9u9eunbt9Bi7uZnQP0xw+jHOicu6+MbVJ6EdOGDf5IWlzsq1WD\nc87xt9NOC9emXVTkV6X5y1/8QeeBB2CPPcJkEcl0kyfD5Zf7Nun+/X37dCgrVvhaMn68L+xHH+0v\nMmzdGo45Jn0HnyiLO865nb7hC/pSYH+gKjAXOKyM7Vy6FBU5N2+ec//4h3OtWjlXs6ZzLVs617ev\nc7NnO1dY6LebNGlSyjJs2eLcyJHONW3q3IknOjd16s79XCozVUQm5lKmnZNtmX74wbkHH3Rur72c\n697duYKC9GT67jvn7r9/krv5ZueaNHFu772du/RS/3f81VfpyVCWRO1Mqi5v65Zsl0ZzYIlz7hPn\n3BZgJNAukqNMOZnBUUfBzTf7s/nPP/dt3P/5D1x0ke+dvvRSePDBfBYu9GfXUVmyxM+Ed8ghfhTA\nXXf5KU9btNi5n8/Pz48uTIQyMZcy7Zxsy1S1KtxwAyxd6v+OWrb0n8BHjPBDqaOyaZPvu7vrLr+P\nffaB/v3zqVnTN51+/jk884yvGXXqRLffkJId514fWFni36vwBT9j7LYbnHeev4HvDJkwwbftnX++\nH47YogU0awaHH+57tg891P/c9mzd6p9rwQL/sW3SJPjvf/2bYdQo30suIuWz555w++2+0I8d6wvt\nFVf4vrUzz/R/X0cdBXvttf3nKSryY8mXLPHDrGfP9qPvFi+GI47wzba9e8Mpp8D998Odd6bn/xdC\n7C9+P/BAuOoqf2Tu08f31k+dCh98AK++Cv36+V/8Lrv4XvB99vH3i4dprVvnDwirV/ue/SOO8G+M\np57yY/I1fYBIdGrUgC5d/G39ej865Z13/PUh8+f7seF168Lee/u+tkqV/PDEtWv97b//9SPrDj7Y\nfxJo1gx69vTt5rvuGvp/l15JdaiaWQugj3PunMS/e+HbiO4rtZ2mhBQRKQcXYrSMmVUGPgJaAauB\n94EuzrmCKMKIiEg0kmpUcM4Vmtm1wAR+Ggqpwi4ikmFSchGTiIiEVaHZHczsGDObamZzzOx9M/tN\nie/1NrMlZlZgZmeVeLyZmc0zs8Vm1r8i+99Btj8m9j3fzO7NoFw3mVmRmdUp8ViQTGbWL7HPuWb2\nkpntUeJ7QV+nEvs7x8wWJfZ3a6r3V2K/DczsHTP7MPEeui7xeG0zm2BmH5nZm2ZWq8TPlPmapSBb\nJTObbWavZFCmWmY2KrGfD83shNC5zOwGM1uQeL8+Z2bV0p3JzAaa2Rozm1fisaQzlOvvriKD5IE3\ngbMS988FJiXuNwHm4Jt9GuEvfCr+lDAdOD5xfzxwdlSD9kvkysM3HVVJ/PtXia+HB87VAHgD+Bio\nEzoTcAZQKXH/XqBvJvz+SuTbqYvmUrTvfYGmifs18X1NhwH3AX9OPH4rcO+OXrMUZLsBeBZ4JfHv\nTMg0BOiRuF8FqBUyF7AfsByolvj380D3dGcCTgGaAvNKPJZ0hvL83VV0XraixC8RYE/gs8T9tsBI\n59xW59wKYAnQ3Mz2BXZ3zs1IbDcU+J8KZijLH/Av2FYA59yXicfbBc71EHBLqceCZXLOTXTOFV/W\nNQ1/8IHwv79iwS6ac8597pybm7j/LVCAf33aAc8kNnuGn/7/Zb5mUecyswZAa+DpEg+HzrQHcKpz\nbjBAYn/JH84ZAAADVElEQVQbQucCKgO7mVkVoDq+PqU1k3NuMrCu1MNJZSjv311Fi/sNwP1m9inQ\nD+ideLz0xU6fJR6rj7/wqdiqxGNRawz81symmdkkMzsudC4zawusdM7NL/Wt0K9Vsd/jzwgyKVNZ\nF82lcn9lMrNG+LOvaUBd59wa8AcAoHhOwG29ZlErPkEo2VkWOtMBwJdmNjjRXDTAzGqEzOWc+w/w\nAPBp4vk3OOcmhsxUwj5JZijX390OR8uY2VtA3ZIP4d9Yt+M/1l/vnBtjZh2AQcCZO3rOKGwn11/w\n/6/azrkWZnY8MAo4MHCm20jTa7OTmW53zr2a2OZ2YItzbkS682U6M6sJvIh/n39rv7yGI20jEszs\nPGCNc26umeVtZ9N0j5KoAjQDrnHOzTSzh4BeZeRI52u1J/4MeX9gAzDKzC4OmWk7UpJhh8XdObfN\ngmRmw5xz1ye2e9HMij8qfgY0LLFpg8Rj23o8aTvIdRXwcmK7GWZWaGZ7Jfb161Tl2lYmMzsS34b2\ngZlZ4vlnm1nzUJlKZPsd/mP+6SUeTvnvbydt67VJi8TH+ReBYc65sYmH15hZXefcmsTH5f+WyJrq\n1+ZkoK2ZtcY3M+xuZsOAzwNmAn8mudI5V7wU9kv44h7ytToDWO6cWwtgZqOBkwJnKpZshvJlq2Bn\nwYdAy8T9VsCMUh0D1fAf2Up2DEzDt2UZvhngnIp2WpSR6wrgrsT9xsAnmZCrRL6P8Z8sgmYCzkn8\nDvcq9XimvE6V+alDtRq+Q/XwVO2vjP0PBR4s9dh9wK2J+2V1hv3iNUtRtpb81KHaL3Qm4F2gceL+\nnYnXKdhrlXiPzgd2TbxXhwDXhMiEP7GbX5H3UHn+7ioa+iRgZiLQVODYEt/rnQhXQGJETeLx4xIv\n+hLgnyl6o1UFhiX2M5PEASh0rhL7Wk5itEzITInn/QSYnbg9FjpTGRnPwY9UWQL0SvX+Suz3ZKAQ\nf0CZk3h9zgHqABMTmSYAe+7oNUtRvpLFPXgm4BhgRuL1ehk/0CJoLvxBpgCYh++4rJruTMBw4D/A\n9/j2/x5A7WQzlOfvThcxiYjEUBYuUSsiIjui4i4iEkMq7iIiMaTiLiISQyruIiIxpOIuIhJDKu4i\nIjGk4i4iEkP/D75srivS7VPFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x12a28f250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, [abs(f(z))**2 for z in x])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens when we change the indices of refraction for the different modes? The phase-mismatch will shift depending on where the new $\\delta k = 0$ occurs. The width of the peak may also change if the indices of refraction are large."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 605,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "indices_of_refraction=[3.,5.,10.]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 606,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f = lambda z: functions.make_nonlinear_interaction(roots_to_use(z),\n",
    "                modes_to_use,Ex.delays,0,0,0.1,plus_or_minus_arr,indices_of_refraction=indices_of_refraction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 607,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12b4a1910>]"
      ]
     },
     "execution_count": 607,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEACAYAAAB/BTv2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHsNJREFUeJzt3X+UXGWd5/H3p9NpCMhE/JFmJIBAYkhYR5A14Khrjz9C\nEkeambOyyegosOdsVsjujOsRkoE5hF2PGH+sLsthIg4zExzdEHRdMntYCBEbxxkhKEQiCaEZJJBA\ngoCABNDuznf/eG6RotLdVd1Wdz23+LzO6VO3bj3Pre+tru5PP/e5t1oRgZmZWTN1tLoAMzNrPw4X\nMzNrOoeLmZk1ncPFzMyazuFiZmZN53AxM7OmayhcJC2UdL+kByRdPEKbKyX1S9oi6ZR6fSUdKWmj\npB2SbpE0veqxlcW2tktaULV+qqSvFX22Sfqj8e22mZlNpLrhIqkDuAo4EzgZWCrppJo2i4ATI2I2\nsAxY00DfFcCmiJgD3AasLPrMA84B5gKLgKslqehzCbA3IuZExDzg9vHuuJmZTZxGRi7zgf6I2BkR\nA8A6oLemTS9wHUBE3AlMl9Rdp28vsLZYXgucXSyfBayLiMGIeBjoL7YDcD5wReVJI+LpRnfUzMwm\nTyPhcjTwaNX9XcW6RtqM1rc7IvYCRMQeYMYI29oNHF112Oyzkn4i6XpJb2ygfjMzm2QTNaGv+k0O\nUu9zaDqBmcAPI+I04A7gy+N4HjMzm2CdDbTZDRxbdX9msa62zTHDtOkape8eSd0RsVfSUcATo20r\nIp6StC8ivlusv4F0mOwgkvyBaWZm4xAR4xkcHKSRkctdwCxJx0nqApYAG2rabAA+DiDpDOCZ4pDX\naH03AOcWy58Abqxav0RSl6TjgVnA5uKxf5D0B8XyB4BtIxUdEVl9XXbZZS2vwTW1V12uyTU1+6uZ\n6o5cImJI0nJgIymMro2I7ZKWpYfjmoi4SdJiSQ8C+4DzRutbbHo1sF7S+cBO0hliRMQ2SetJwTEA\nXBAH9noF8A1JXwF+UXkeMzPLSyOHxYiIm4E5Neu+VnN/eaN9i/VPk0Yfw/W5gqqzwqrWPwK8t5Ga\nzcysdXyF/iTp6elpdQkHcU2Ny7Eu19QY19QaavZxthxIinbcLzOziSSJmMQJfTMzszFxuJiZWdM5\nXMzMrOkcLmZm1nQOFzMzazqHi5mZNZ3DxczMms7hYmZmTedwMTOzpnO4mJlZ0zlczMys6RwuZmbW\ndA4XMzNrOoeLmZk1ncPFzMyazuFiZmZN53AxawP79sH117e6CrMDHC5mbeDmm2HJklZXYXaAw8Ws\njezf3+oKzBKHi1kbePHFdPvss62tw6zC4WLWBp5/Pt2+8EJr6zCrcLiYtYFKuLz0UmvrMKtwuJi1\nAYeL5cbhYtYGHC6WG4eLWRvYty/dVib2zVrN4WLWBgYH061HLpaLhsJF0kJJ90t6QNLFI7S5UlK/\npC2STqnXV9KRkjZK2iHpFknTqx5bWWxru6QFVeu/X2zrHkl3S3rD+HbbrL0MDaVbh4vlom64SOoA\nrgLOBE4Glko6qabNIuDEiJgNLAPWNNB3BbApIuYAtwEriz7zgHOAucAi4GpJqnq6pRFxakS8PSKe\nHN9um7WXysjFh8UsF42MXOYD/RGxMyIGgHVAb02bXuA6gIi4E5guqbtO315gbbG8Fji7WD4LWBcR\ngxHxMNBfbGcsNZu9qlRGLpWQMWu1Rn5RHw08WnV/V7GukTaj9e2OiL0AEbEHmDHCtnbXPN/fFYfE\nLm2gdrNXhaEh6OyEgYFWV2KWTNQoQPWbHCQaaPMnEfFW4D3AeyR9bBzPY9Z2Bgfh0EMdLpaPzgba\n7AaOrbo/s1hX2+aYYdp0jdJ3j6TuiNgr6SjgiTrbIiIeL273SfoW6XDZ3w9X9KpVq15e7unpoaen\nZ7R9NCu1oaEULj4sZmPR19dHX1/fhGxbEaMPGCRNAXYA7wceBzaTJtW3V7VZDFwYER+SdAbw1Yg4\nY7S+klYDT0fE6uIssiMjYkUxof9N4HTS4bBbgdmkUdZrI+IpSVOBbwG3RsQ1w9Qc9fbLrJ18+MPw\n05/CRRfB8uWtrsbKShIRMZ4jTwepO3KJiCFJy4GNpF/w1xbhsCw9HNdExE2SFkt6ENgHnDda32LT\nq4H1ks4HdpLOECMitklaD2wDBoALIiIkHQLcIqkTmAJsAr7ejBfBrOw8crHc1B25lJFHLvZqs2AB\n7NkDf/qn8JnPtLoaK6tmjlx8Wq9ZG/DIxXLjcDFrAz5bzHLjcDFrAx65WG4cLmZtwCMXy43DxawN\nVEYuDhfLhcPFrA1URi4+LGa5cLiYtQGPXCw3DhezNuCRi+XG4WLWBjxysdw4XMzagEculhuHi1kb\n8MjFcuNwMWsDvojScuNwMWsDvojScuNwMWsDPixmuXG4mLUBT+hbbhwuZm3AIxfLjcPFrA145GK5\ncbiYtYGhIejqSrdmOXC4mLWBwUE45BCHi+XD4WJWchEeuVh+HC5mJbd/P0gwdarnXCwfDhezkhsa\ngs5OmDLFIxfLh8PFrOQGB1OwOFwsJw4Xs5LzyMVy5HAxKzmPXCxHDhezkhsacrhYfhwuZiU3OJgO\ni3V2+mwxy4fDxazkPHKxHDlczErOE/qWo4bCRdJCSfdLekDSxSO0uVJSv6Qtkk6p11fSkZI2Stoh\n6RZJ06seW1lsa7ukBcM81wZJ945tV83akyf0LUd1w0VSB3AVcCZwMrBU0kk1bRYBJ0bEbGAZsKaB\nviuATRExB7gNWFn0mQecA8wFFgFXS1LVc/0R8Nx4d9is3XjkYjlqZOQyH+iPiJ0RMQCsA3pr2vQC\n1wFExJ3AdEnddfr2AmuL5bXA2cXyWcC6iBiMiIeB/mI7SDoc+BTw2bHuqFm78sjFctRIuBwNPFp1\nf1exrpE2o/Xtjoi9ABGxB5gxwrZ2V/X5b8CXgBcbqNvsVaEycunsdLhYPjonaLuq3+QgMeoGpbeR\nDr39F0lvrvccq1atenm5p6eHnp6ecZRklr/qkYtPRbax6Ovro6+vb0K23Ui47AaOrbo/s1hX2+aY\nYdp0jdJ3j6TuiNgr6SjgiTrbeidwmqSHgKnADEm3RcT7hiu6OlzM2pnnXGy8av/wvvzyy5u27UYO\ni90FzJJ0nKQuYAmwoabNBuDjAJLOAJ4pDnmN1ncDcG6x/Angxqr1SyR1SToemAVsjog1ETEzIk4A\n3g3sGClYzF5NPOdiOao7comIIUnLgY2kMLo2IrZLWpYejmsi4iZJiyU9COwDzhutb7Hp1cB6SecD\nO0lniBER2yStB7YBA8AFETHqITOzV7PKyKWjI/3jsIj0/13MWknt+HtbkvPIXjV+8AO45BL4x39M\no5df/zqFjdlYSSIimvKnia/QNyu5ysgFfMaY5cPhYlZylTkX8Bljlg+Hi1nJVY9cPKlvuXC4mJVc\n7cjF4WI5cLiYlVzlI/fB4WL5cLiYlVzln4WBw8Xy4XAxK7nqkYvPFrNcOFzMSq52Qt9ni1kOHC5m\nJecJfcuRw8Ws5HwqsuXI4WJWch65WI4cLmYl55GL5cjhYlZy1SMXny1muXC4mJWczxazHDlczErO\ncy6WI4eLWcl5zsVy5HAxKzmPXCxHDhezkvPIxXLkcDErOZ8tZjlyuJiVnM8Wsxw5XMxKznMuliOH\ni1nJ+Z+FWY4cLmYl538WZjlyuJiVnEculiOHi1nJ+VRky5HDxazkak9F9tlilgOHi1nJeeRiOXK4\nmJWcT0W2HDUULpIWSrpf0gOSLh6hzZWS+iVtkXRKvb6SjpS0UdIOSbdIml712MpiW9slLaha//8k\n3SNpq6SrJWl8u23WPjxysRzVDRdJHcBVwJnAycBSSSfVtFkEnBgRs4FlwJoG+q4ANkXEHOA2YGXR\nZx5wDjAXWARUh8hHIuLUiHgrMAP4yHh33KxdeORiOWpk5DIf6I+InRExAKwDemva9ALXAUTEncB0\nSd11+vYCa4vltcDZxfJZwLqIGIyIh4H+YjtExPMAkqYCXUCMbXfN2o9HLpajRsLlaODRqvu7inWN\ntBmtb3dE7AWIiD2kkchw29pd/XySbgb2AM8B326gfrO25rPFLEedE7Td8cyFNDQKiYiFkrqAbwLv\nA743XLtVq1a9vNzT00NPT884SjLLn0cuNl59fX309fVNyLYbCZfdwLFV92cW62rbHDNMm65R+u6R\n1B0ReyUdBTxRZ1svi4jfSNpAOrRWN1zM2pnnXGy8av/wvvzyy5u27UYOi90FzJJ0XDFiWAJsqGmz\nAfg4gKQzgGeKQ16j9d0AnFssfwK4sWr9Ekldko4HZgGbJR1ehBCSOoEPAfePdYfN2o1HLpajuiOX\niBiStBzYSAqjayNiu6Rl6eG4JiJukrRY0oPAPuC80foWm14NrJd0PrCTdIYYEbFN0npgGzAAXBAR\nIelwYEMRUh3A9ynOSjN7NfM/C7McNTTnEhE3A3Nq1n2t5v7yRvsW658GPjBCnyuAK2rWPUFx1piZ\nHeB/FmY58hX6ZiXnORfLkcPFrOT8kfuWI4eLWclV/7Mwz7lYLhwuZiVXO3LxnIvlwOFiVnI+Fdly\n5HAxKzmfimw5criYlZxPRbYcOVzMSs6nIluOHC5mJec5F8uRw8Ws5DznYjlyuJiVnOdcLEcOF7OS\n85yL5cjhYlZynnOxHDlczErOcy6WI4eLWcl5zsVy5HAxKznPuViOHC5mJec5F8uRw8Ws5DznYjly\nuJiVWATs3++P3Lf8OFzMSmxoCDo6QEr3fVjMcuFwMSux6n8UBg4Xy4fDxazEqv/FMaRlHxazHDhc\nzErMIxfLlcPFrMSqT0MGh4vlw+FiVmLVpyGDT0W2fDhczEpsuJGL51wsBw4XsxKrHbn4sJjlwuFi\nVmKec7FcNRQukhZKul/SA5IuHqHNlZL6JW2RdEq9vpKOlLRR0g5Jt0iaXvXYymJb2yUtKNZNk/R/\ni3VbJX1u/Ltt1h4852K5qhsukjqAq4AzgZOBpZJOqmmzCDgxImYDy4A1DfRdAWyKiDnAbcDKos88\n4BxgLrAIuFqqXH/MFyNiLnAq8G5JZ453x83agedcLFeNjFzmA/0RsTMiBoB1QG9Nm17gOoCIuBOY\nLqm7Tt9eYG2xvBY4u1g+C1gXEYMR8TDQD8yPiBcj4vbiOQaBu4GZY91hs3biORfLVSPhcjTwaNX9\nXcW6RtqM1rc7IvYCRMQeYMYI29pd+3ySXgt8GPheA/WbtS3PuViuOus3GRfVb3KQaGjD0hTgW8BX\ni5HNsFatWvXyck9PDz09PeMoySxvnnOx30ZfXx99fX0Tsu1GwmU3cGzV/ZnFuto2xwzTpmuUvnsk\ndUfEXklHAU/U2VbFNcCOiPifoxVdHS5m7cpzLvbbqP3D+/LLL2/aths5LHYXMEvScZK6gCXAhpo2\nG4CPA0g6A3imOOQ1Wt8NwLnF8ieAG6vWL5HUJel4YBawudj2Z4HfiYhPjXlPzdqQ51wsV3VHLhEx\nJGk5sJEURtdGxHZJy9LDcU1E3CRpsaQHgX3AeaP1LTa9Glgv6XxgJ+kMMSJim6T1wDZgALggIkLS\n0cBfANsl3UM6jHZVRPxNs14Ms7LxnIvlShENTXWUiqRox/0yq3X77fCXfwk/+EG6/9xzMHNmujUb\nK0lExHjmzA/iK/TNSmy4j9z3nIvlwOFiVmK1/yzMh8UsFw4XsxKrHbn4VGTLhcPFrMRGmtD3lKO1\nmsPFrMRqT0WW0tf+/a2ryQwcLmalVjtyAc+7WB4cLmYlVjtyAc+7WB4cLmYlVjuhDz4d2fLgcDEr\nscFBmDr1let8WMxy4HAxK7Ha61zA4WJ5cLiYldhw4eI5F8uBw8WsxEYauXjOxVrN4WJWYj4sZrly\nuJiVmMPFcuVwMSsxX+diuXK4mJXYSFfoe87FWs3hYlZiPixmuXK4mJWYw8Vy5XAxK7GRrnPxYTFr\nNYeLWYl55GK5criYlZjDxXLlcDErMX/8i+XK4WJWYv74F8uVw8WsxIa7iNKHxSwHDhezEvOci+XK\n4WJWYsNdoe85F8uBw8WsxDznYrlyuJiVmA+LWa4aChdJCyXdL+kBSReP0OZKSf2Stkg6pV5fSUdK\n2ihph6RbJE2vemxlsa3tkhZUrf+spEckPTe+3TVrLw4Xy1XdcJHUAVwFnAmcDCyVdFJNm0XAiREx\nG1gGrGmg7wpgU0TMAW4DVhZ95gHnAHOBRcDVklT02QC8Y9x7a9ZmfJ2L5aqRkct8oD8idkbEALAO\n6K1p0wtcBxARdwLTJXXX6dsLrC2W1wJnF8tnAesiYjAiHgb6i+0QEZsjYu/Yd9OsPXnOxXLVSLgc\nDTxadX9Xsa6RNqP17a4ERUTsAWaMsK3dwzyfmTF8uEydCgMDranHrGKiJvRVv8lBoulVmLU5h4vl\nqrN+E3YDx1bdn1msq21zzDBtukbpu0dSd0TslXQU8ESdbY3JqlWrXl7u6emhp6dnrJswy95wV+g7\nXKxRfX199PX1Tci2FTH6gEHSFGAH8H7gcWAzsDQitle1WQxcGBEfknQG8NWIOGO0vpJWA09HxOri\nLLIjI2JFMaH/TeB00uGwW4HZUVWopF9FxBGj1Bz19susHbznPfC5z6XbigsvhLlzYfny1tVl5SSJ\niBjPkaeD1B25RMSQpOXARtJhtGuLcFiWHo5rIuImSYslPQjsA84brW+x6dXAeknnAztJZ4gREdsk\nrQe2AQPABZWkKALpT4Bpkh4B/joi/mszXgizMhrusFhXF/zmN62px6yi7siljDxysVeLd7wDrr46\n3VZcdBG8/vVw8bBXpJmNrJkjF1+hb1ZintC3XDlczErM4WK5criYlZjDxXLlcDErsYEBT+hbnhwu\nZiU2MJDCpJpHLpYDh4tZif3mNw4Xy5PDxazEHC6WK4eLWYk5XCxXDhezEnO4WK4cLmYlFZHCZerU\nV6732WKWA4eLWUlVrnHpqPkp9sjFcuBwMSup4Q6JgcPF8uBwMSsph4vlzOFiVlIOF8uZw8WspIab\nzAdP6FseHC5mJeWRi+XM4WJWUg4Xy5nDxaykhvvQSnC4WB4cLmYlNdLIxXMulgOHi1lJjRQuhx0G\nL744+fWYVXO4mJXUSOEybRq88MLk12NWzeFiVlIvvTTyyMXhYq3mcDErmRdeSB9a+cILcPjhBz9+\n6KHw61/D/v3pdnBw8ms0c7iYlUhECpQ1a2DfvuHDpaMjBcxLL8Gpp8IFF0x+nWYOF7MS2bUr3d56\n68gjF0iHxp55BrZvh1tumbz6zCocLmYlsnUrdHfD/fenkcthhw3fbtq0FCzHHANPPOGzx2zyOVzM\nSmTrVvjjP4aHHoJf/Wr0kcvPfgZz5sCsWbBjx+TWaeZwMSuRrVth/nx44xtTYNQLl+OPh3nz4L77\nJrdOs4bCRdJCSfdLekDSxSO0uVJSv6Qtkk6p11fSkZI2Stoh6RZJ06seW1lsa7ukBVXr3y7p3mJb\nXx3fLpuV19at8Na3wgknwL33jnxY7LDDUqCccAKcfDJs2za5dZrVDRdJHcBVwJnAycBSSSfVtFkE\nnBgRs4FlwJoG+q4ANkXEHOA2YGXRZx5wDjAXWARcLUlFn78C/n1EvAV4i6Qzx7vjk62vr6/VJRzE\nNTWuVXU9/3wafXzxi+nzwvr700jkxBPhvvv6OOKI4fu95jXw05+mviefnEYxAL298NGPTly9OX7/\nXFNrNDJymQ/0R8TOiBgA1gG9NW16gesAIuJOYLqk7jp9e4G1xfJa4Oxi+SxgXUQMRsTDQD8wX9JR\nwBERcVfR7rqqPtnL8c3kmho3mXXdfjvs3JmW16yBGTPgiivgn/85hcW0aSlcoI8ZM4bfxowZ6Wyy\nE09M4XLfffCTn8A//RNs2gR3353a3Xsv/Mu/NK/2HL9/rqk1GgmXo4FHq+7vKtY10ma0vt0RsRcg\nIvYAlR+T2j67q7a1q04dZlmIOLD81FPw6KNpXQTccQd873tpJPLUU/DpT8Nll6XrUr70JVi6FN75\nzjSn8uUvw9e/DgsXwic/CaedlraZwiWdOTacSujMnp0m9B9/HC69FFauhIsuSiOhTZvgXe+C009P\nwbN3Lyxfnh4bHEzh9N3vpjPTKn7xi1RnhT8g00bSOUHbVf0mB4n6TRr3h3/4228jmljRAw/A5s3j\n768GX9Gx1Nzf/8qaKn0b2YZ0oKba2+Fqqv2qtO/oSLf796f1/f3pr+vqx4eG0i+7KVPSV+39/fvT\nV3VtEWmdlNoMDh7o09mZfikODqaPp+/sTKfqDg6miw87OtJpvkND6fBSRPr+fec7MH166vvkk+mj\nV97wBnj22fTLe/r09Et9z550Pcqb3pTW7dyZtjttWnquKVPS+kceSc/xkY+k5VmzUr2bN8P118Pc\nubBkCfze78Gf/3kKgUsvTft4UnFw+YQThn/Nf/d30+30Yibzgx+EG2+EdevSc6xeDTffDDfdlAJu\nwYK0/mMfS8G3Zk3ar7e9LY16jjkGfvnLtN8Raf1jj8HDD6eTC+bMSa/Bz38ON9yQQu+xx9I2Dj0U\njjgiva7PPZdeg0MOSbcDA2n58MPT61/5SJtDD02vTeW/bR5ySLo/NPTK7/vQ0IHXtPIeqrw3K++D\nhx6CH/7wlev270/1VPoNDr7yfTnae3ykddXv/9qfpdrbBx+EH/3olX1rtzGcRn8PZCEiRv0CzgBu\nrrq/Ari4ps0a4N9V3b8f6B6tL7CdNHoBOArYPtz2gZuB06vbFOuXAH81Qs3hL3/5y1/+GvtXvUxo\n9KuRkctdwCxJxwGPk36pL61pswG4ELhe0hnAMxGxV9KTo/TdAJwLrAY+AdxYtf6bkr5COuw1C9gc\nESHpWUnzi5o+Dlw5XMERUaZ8NzNrO3XDJSKGJC0HNpLmaK6NiO2SlqWH45qIuEnSYkkPAvuA80br\nW2x6NbBe0vnATtIZYkTENknrgW3AAHBBxMsHai4E/g44FLgpIm5uwmtgZmZNpgO/t83MzJqj1Ffo\nS3qbpB9JukfSZkn/uuqxll2IKek/Fc+7VdLnc6ipeJ5PS9ov6XU51CTpC8XzbpH0HUm/k0NdNTXW\nvYB4gp53pqTbJN1XvI/+c7F+zBcfT0BtHZLulrQhh5okTZd0Q/Ec90k6PYOaPiXpZ8V79ZuSulpR\nk6RrJe2VdG/Vusm5gL1Zkzet+AJuARYUy4uA7xfL84B7SIf93gw8yIFR2p3AO4rlm4Azm1xTD+kw\nYGdx/w3F7dxW1VRsdybp5IifA6/LpKYPAB3F8ueBK1r9/aupr6N47uOAqcAW4KRJem8fBZxSLL8G\n2AGcRDqcfFGx/mLg8/Veswmo7VPA3wMbivstrYl0qPy8YrkTmN7KmoA3AQ8BXcX960nzypNeE/Bu\n4BTg3qp1Y65jPD93pR65APtJbySA15KuiYHWXoj5SdI3axAgIp4s1ve2sCaArwCfqVnX0poiYlNE\nVE4ivoMUgJDPhbSNXEA8ISJiT0RsKZafJ51dOZMxXnzc7LokzQQWA39dtbplNRWj3fdExN8CFM/1\nbCtrKkwBDpfUCUwj/W6a9Joi4ofAL2tWT8oF7GUPl08BX5L0CPAFio+QobUXYr4F+DeS7pD0fUmn\ntbomSWcBj0bE1pqHcrpg9XzSX0Q51dXIBcQTTtKbSX993sHYLz5utsofKdWTta2s6XjgSUl/Wxyq\nu0bSYa2sKSIeA74MPFJs/9mI2NTKmmrMGGMd4/q5m6iLKJtG0q2ka2ZeXkV6Y19COqzyZxHxfyT9\nW+BvgA+2sKZLSa/pkRFxhqR3ADcAI1zqNmk1/QWT8LqMsa5LIuIfijaXAAMR8b9aUGLWJL0G+Dbp\nff68pNozcCbtjBxJHwL2RsQWST2jNJ3Ms4Q6gbcDF0bEj5UuYVgxTA2T+Tq9ljQ6OA54FrhB0kdb\nWVMdE1JH9uESESP+UpT0jYj4s6LdtyVVhuq7gWOqms4s1o20vpk1/Ufgfxft7pI0JOn1xfMcO9k1\nSfpXpOOnP5WkYvt3K10vNKE1jVZXVX3nkg6zvK9q9YR+/8ZgpNdnUhSHVL4NfCMiKteB7ZXUHek6\nsqOAJ6pqnejX5l3AWZIWkw71HCHpG8CeFta0izQq/3Fx/zukcGnl6/QB4KGIeBpA0neB329xTdXG\nWsf46mvmRNZkfwH3Ae8tlt8P3BWvnJjqIg2bqyem7iAdzxTpMMzCJtf0H4DLi+W3ADtbXVNNfT8n\njaxaXhOwsPgevr5mfS6v1RQOTOh3kSb0507i+/s64L/XrFvNgU+5GG4y9qDXbIJqey8HJvS/0Mqa\ngNuBtxTLlxWvUctep+L9uZV0PZ5IJxxc2KqaSH9cbv1t3kPj+bmblB+SCXyD/z7w4+IF+RFwatVj\nK4sXZzvFGWXF+tOKb3w/8D8moKapwDeK5/gxRfi1sqaa+h6iOFus1TUV294J3F18XZ1DXTU1LiSd\nqdUPrJjo56t63ncBQ6RAu6d4fRYCrwM2FTVtBF5b7zWboPqqw6WlNQFvI31qxxbSUYPpGdR0WbH9\ne0mT5lNbURPwLeAx4NekOaDzgCPHWsd4fu58EaWZmTVd2c8WMzOzDDlczMys6RwuZmbWdA4XMzNr\nOoeLmZk1ncPFzMyazuFiZmZN53AxM7Om+//ffJOfubf36QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x12b042250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x, [abs(f(z))**2 for z in x])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating a Hamiltonian from  a model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " In this section we will use example 3 to generate a Hamiltonian with nonlinaer coefficients resulting when inserting a nonlinearity in a circuit. We will assume that the nonlinearity is inserting at the delay line of index 0 corresponding to $\\tau_1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 611,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "import itertools\n",
    "from qnet.algebra.circuit_algebra import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 612,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Ex = Time_Delay_Network.Example3(r1 = 0.9, r3 = 0.9, max_linewidth=35.,max_freq=25.)\n",
    "Ex.run_Potapov()\n",
    "E = Ex.E\n",
    "roots = Ex.roots\n",
    "M1 = Ex.M1\n",
    "delays = Ex.delays\n",
    "modes = functions.spatial_modes(roots,M1,E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 613,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(-0.36514939263783197-10.901513851604911j),\n",
       " (-0.33769377381866061+2.0232070709165708e-21j),\n",
       " (-0.36514939263783225+10.901513851604911j),\n",
       " (-0.33972600185569901+20.178467995423109j),\n",
       " (-0.3397260018556984-20.17846799542311j)]"
      ]
     },
     "execution_count": 613,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "roots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 614,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## nonlinearity information\n",
    "\n",
    "delay_index = 0 \n",
    "start_nonlin = 0.\n",
    "duration_nonlin = .1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 806,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "NONLIN_WEIGHT = 10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 615,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m = len(roots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 616,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "indices = range(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 617,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "chi_order = 3  ## i.e. chi-3 nonlinearity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 618,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plus_minus_combinations = list(itertools.combinations(range(chi_order + 1), 2))  ## pick which fields are annihilated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 619,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "list_of_pm_arr = []\n",
    "for tup in plus_minus_combinations:\n",
    "    ls = [1]*(chi_order+1)\n",
    "    for i in tup:\n",
    "        ls[i]=-1\n",
    "    list_of_pm_arr.append(ls)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 620,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = [sp.symbols('a_'+str(i)) for i in range(m)]\n",
    "a_H = [sp.symbols('a^H_'+str(i)) for i in range(m)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 621,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A,B,C,D = Potapov.get_Potapov_ABCD(Ex.roots,Ex.vecs,Ex.T,z=0.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 622,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Omega = (A-A.H)/(2j) #### closed dynamics only. i.e. not damping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 623,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Omega = -1j*A  ## full dynamics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 624,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "H_lin_sp = 0\n",
    "## with sympy only\n",
    "for i in range(m):\n",
    "    for j in range(m):\n",
    "        H_lin_sp += a_H[i]*a[j]*Omega[i,j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 625,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_nonlin_term_sp(combination,pm_arr):\n",
    "    '''\n",
    "    Make symbolic term\n",
    "    With sympy only\n",
    "    '''\n",
    "    r = 1\n",
    "    for index,sign in zip(combination,pm_arr):\n",
    "        if sign == 1:\n",
    "            r*= a_H[index]\n",
    "        else:\n",
    "            r *= a[index]\n",
    "    return r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's impose a large 'index of refraction'. In the future we will replaces this by better conditions for phase-mismatch, including realistic values. For now, this will narrow the gain versus $\\Delta k$ function so that few interaction terms remain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 626,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def weight(combination,pm_arr):\n",
    "    roots_to_use = np.array([roots[i].imag for i in combination])\n",
    "    modes_to_use = [modes[i] for i in combination]\n",
    "    return functions.make_nonlinear_interaction(roots_to_use, modes_to_use, delays, delay_indices,\n",
    "                                start_nonlin,duration_nonlin,pm_arr,\n",
    "                                indices_of_refraction = [1000.]*len(combination),\n",
    "                                eps=1e-12,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 627,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## TODO: add a priori check to restrict exponential growth\n",
    "weights = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 628,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "delta_k is less than epsilon, approximating linearly.\n",
      "420\n"
     ]
    }
   ],
   "source": [
    "count = 0\n",
    "\n",
    "for pm_arr in list_of_pm_arr:\n",
    "    field_combinations = itertools.combinations_with_replacement(range(m), chi_order+1)\n",
    "    for combination in field_combinations:\n",
    "        count += 1\n",
    "        weights[tuple(combination),tuple(pm_arr)] = weight(combination,pm_arr)        \n",
    "print count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 629,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEACAYAAAC+gnFaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEkJJREFUeJzt3WGMXNd53vH/Q9MrKBVWke1KDKhEUqvYpg24lNGybl20\nDNLIUhY1BX8QlASoVEGoENWJ0RSNSRfFqkEAW0HtIkWgL7UdMEUFmXWbWk5jSRakAdZFRCmxFNMm\no7KoKYtsSLOGkVYxSnDLtx/2EhpvV9wZzlzu7pz/DxjgzrnnzD0vl7zP3HNnuKkqJEnt2rbRE5Ak\nbSyDQJIaZxBIUuMMAklqnEEgSY0zCCSpcesGQZKrkhxO8lKSI0kWu/bFJCeTfL173DE05kCS40mO\nJbm9zwIkSZPJKN8jSPIjVfWDJG8B/gvwy8CdwP+uqs+s6rsLeAz4a8CNwDPAT5ZfWJCkTWmkpaGq\n+kG3eRWwHbh4Us8a3fcBj1fVclWdAI4DeyacpySpJyMFQZJtSV4CTgNfraoXu10fTfJyks8mubZr\n2wm8NjT8VNcmSdqERr0iuFBVt7Gy1LMnyXuAR4G/VFW7WQmIT/c3TUlSX7aP07mq/leSAXDHqnsD\n/wb4crd9CvjxoX03dm0/JIn3DCTpMlTVWsvyl22UTw294+KyT5KrgZ8B/iTJjqFuHwG+2W0/AdyT\nZC7JLcCtwAtrvXZVzexjcXFxw+dgfdbXYn2zXFtVP++fR7ki+DHgYJJtrATHF6rq95P8TpLdwAXg\nBPBgd3I/muQQcBQ4DzxUbzL7733vewBcddVVXHPNNZPWIkm6DOsGQVUdAd6/Rvvfv8SYTwKfXO+1\nd+58JwBvfWvx3e+e4uqrr15viCRpysa6RzBt586tXBFUzXP+/PmZCoK9e/du9BR6ZX1b2yzXN8u1\n9WWkL5T1cuCkLn4dYW5unrNnTzI/P78hc5GkrSIJdaVvFkuSZptBIEmNMwgkqXEGgSQ1ziCQpMYZ\nBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEg\nSY0zCCSpcQaBJDXOIJCkxq0bBEmuSnI4yUtJjiRZ7NqvS/J0kleSPJXk2qExB5IcT3Isye19FiBJ\nmsy6QVBV54CfqqrbgN3AnUn2APuBZ6rqXcCzwAGAJO8B7gZ2AXcCjyZJT/OXJE1opKWhqvpBt3kV\nsB0oYB9wsGs/CNzVbX8YeLyqlqvqBHAc2DOtCUuSpmukIEiyLclLwGngq1X1InBDVZ0BqKrTwPVd\n953Aa0PDT3VtkqRNaPsonarqAnBbknngd5O8l5Wrgh/qNv7hHwZgefkcS0tLLCwsjP8SkjTDBoMB\ng8Gg12Okarzzd5J/DvwAeADYW1VnkuwAnquqXUn2A1VVj3T9nwQWq+rwqtepi9kxNzfP2bMnmZ+f\nn7wiSZphSaiqqd53HeVTQ++4+ImgJFcDPwMcA54A7uu63Qt8qdt+ArgnyVySW4BbgRemOWlJ0vSM\nsjT0Y8DBJNtYCY4vVNXvJ3keOJTkfuBVVj4pRFUdTXIIOAqcBx6qcS87JElXzNhLQ1M7sEtDkjS2\nDVkakiTNNoNAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklq\nnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1bt0gSHJj\nkmeTfCvJkSS/1LUvJjmZ5Ovd446hMQeSHE9yLMntfRYgSZrM9hH6LAO/UlUvJ7kG+KMkX+32faaq\nPjPcOcku4G5gF3Aj8EySn6yqmubEJUnTse4VQVWdrqqXu+3XgWPAzm531hiyD3i8qpar6gRwHNgz\nnelKkqZtrHsESW4GdgOHu6aPJnk5yWeTXNu17QReGxp2ijeCQ5K0yYyyNARAtyz0ReBjVfV6kkeB\nX6uqSvLrwKeBB8Y7/MMALC+fY2lpiYWFhfGGS9KMGwwGDAaDXo+RUZbuk2wHfg/4SlX95hr7bwK+\nXFXvS7IfqKp6pNv3JLBYVYdXjSlYOfbc3Dxnz55kfn5+4oIkaZYloarWWpa/bKMuDX0eODocAkl2\nDO3/CPDNbvsJ4J4kc0luAW4FXpjGZCVJ07fu0lCSDwK/ABxJ8hIrb+M/Afx8kt3ABeAE8CBAVR1N\ncgg4CpwHHvITQ5K0eY20NNTLgV0akqSxbeTSkCRpRhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgk\nqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIa\nZxBIUuMMAklqnEEgSY0zCCSpcesGQZIbkzyb5FtJjiT55a79uiRPJ3klyVNJrh0acyDJ8STHktze\nZwGSpMmMckWwDPxKVb0X+BvAP0rybmA/8ExVvQt4FjgAkOQ9wN3ALuBO4NEk6WPykqTJrRsEVXW6\nql7utl8HjgE3AvuAg123g8Bd3faHgcerarmqTgDHgT1TnrckaUrGukeQ5GZgN/A8cENVnYGVsACu\n77rtBF4bGnaqa5MkbULbR+2Y5Brgi8DHqur1JLWqy+rnI3gYgOXlcywtLbGwsDD+S0jSDBsMBgwG\ng16Pkar1z99JtgO/B3ylqn6zazsG7K2qM0l2AM9V1a4k+4Gqqke6fk8Ci1V1eNVr1sXsmJub5+zZ\nk8zPz0+zNkmaOUmoqqnedx11aejzwNGLIdB5Ariv274X+NJQ+z1J5pLcAtwKvDCFuUqSerDu0lCS\nDwK/ABxJ8hIrb+M/ATwCHEpyP/AqK58UoqqOJjkEHAXOAw/VKJcdkqQNMdLSUC8HdmlIksa2kUtD\nkqQZZRBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgk\nqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGrRsEST6X5EyS\nbwy1LSY5meTr3eOOoX0HkhxPcizJ7X1NXJI0HaNcEfw28KE12j9TVe/vHk8CJNkF3A3sAu4EHk2S\nqc1WkjR16wZBVX0N+P4au9Y6we8DHq+q5ao6ARwH9kw0Q0lSrya5R/DRJC8n+WySa7u2ncBrQ31O\ndW2SpE1q+2WOexT4taqqJL8OfBp4YPyXeRiA5eVzLC0tsbCwcJnTkaTZNBgMGAwGvR4jVbV+p+Qm\n4MtV9b5L7UuyH6iqeqTb9ySwWFWH1xhXsHLsubl5zp49yfz8/GTVSNKMS0JVTfXe66hLQ2HonkCS\nHUP7PgJ8s9t+ArgnyVySW4BbgRemMVFJUj/WXRpK8hiwF3h7ku8Ai8BPJdkNXABOAA8CVNXRJIeA\no8B54KEa5ZJDkrRhRloa6uXALg1J0tg2cmlIkjSjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLU\nOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0z\nCCSpcQaBJDXOIJCkxhkEktS4dYMgyeeSnEnyjaG265I8neSVJE8luXZo34Ekx5McS3J7XxOXJE3H\nKFcEvw18aFXbfuCZqnoX8CxwACDJe4C7gV3AncCjSTK96UqSpm3dIKiqrwHfX9W8DzjYbR8E7uq2\nPww8XlXLVXUCOA7smc5UJUl9uNx7BNdX1RmAqjoNXN+17wReG+p3qmuTJG1S26f0OnV5wx4GYHn5\nHEtLSywsLExpOpI0GwaDAYPBoNdjpGr9c3iSm4AvV9X7uufHgL1VdSbJDuC5qtqVZD9QVfVI1+9J\nYLGqDq/xmnUxP+bm5jl79iTz8/NTK0ySZlESqmqq915HXRpK97joCeC+bvte4EtD7fckmUtyC3Ar\n8MIU5ilJ6sm6S0NJHgP2Am9P8h1gEfgU8O+T3A+8ysonhaiqo0kOAUeB88BDNcolhyRpw4y0NNTL\ngV0akqSxbeTSkCRpRhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINA\nkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSp\ncdsnGZzkBPBnwAXgfFXtSXId8AXgJuAEcHdV/dmE85Qk9WTSK4ILwN6quq2q9nRt+4FnqupdwLPA\ngQmPIUnq0aRBkDVeYx9wsNs+CNw14TEkST2aNAgK+GqSF5M80LXdUFVnAKrqNHD9hMeQJPVoonsE\nwAer6k+T/EXg6SSvsBIOw1Y/H/IwAMvL51haWmJhYWHC6UjSbBkMBgwGg16PkapLnKfHeaFkEXgd\neICV+wZnkuwAnquqXWv0r4sZMTc3z9mzJ5mfn5/KXCRpViWhqjLN17zspaEkP5Lkmm77LwC3A0eA\nJ4D7um73Al+acI6SpB5NsjR0A/C7K+/s2Q78u6p6OskfAoeS3A+8Ctw9hXlKknoytaWhsQ/s0pAk\njW1TLQ1JkmaDQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaB\nJDXOIJCkxhkEktQ4g0CSGmcQaFPaseNmkpCEHTtu3ujpSDNt0l9eL/XizJlXufiLi86cmerv4JC0\nilcEV5jvdKXNq9V/n/6qyissCRfrhrBRf/6bnX9O2ghb4e+dv6pSkjR1BoG2lOFL99Yu37U5zcJy\nkktDV9hWuPTcDN7sz+mH2394nzSpy/n3eaX/TW+ppaEkdyT5kyT/NcnH+zqOJGkyvQRBkm3AbwEf\nAt4L/FySd/dxrM1qMBhs9BR6ZX1b2yzXtxVq22zLSX1dEewBjlfVq1V1Hngc2NfTsTalrfCXcRLW\nt7XNcn1bobY3vidT3fbG6isIdgKvDT0/2bVJkjaZDf1m8fz83wPgz//8/7Btmx9gkqSN0MunhpJ8\nAHi4qu7onu8HqqoeGerjRz0k6TJM+1NDfQXBW4BXgJ8G/hR4Afi5qjo29YNJkibSy9JQVf3fJB8F\nnmblPsTnDAFJ2pw27AtlkqTNYWp3aEf5AlmSf53keJKXk+xeb2yS65I8neSVJE8luXZa8x1HT7X9\nRpJjXf//kGTDvlbdR31D+/9JkgtJ3tZnDZfSV31Jfqn7GR5J8qm+63gzPf39/CtJ/iDJS0leSPJX\nr0Qta7mM+m4bav9ckjNJvrGq/1Y9t4xS2/jnlqqa+MFKoPw34CbgrcDLwLtX9bkT+M/d9l8Hnl9v\nLPAI8Kvd9seBT01jvpuktr8LbOu2PwV88krX1md93f4bgSeBbwNvm6X6gL2sLH1u756/Y8bqewq4\nfWj8c1utvu753wJ2A99YNWZLn1vWqW3sc8u0rghG+QLZPuB3AKrqMHBtkhvWGbsPONhtHwTumtJ8\nx9FLbVX1TFVd6MY/z8pJcyP09bMD+FfAP+27gHX0Vd8vsnLyWO7G/c/+S1lTX/VdAC6+S/5R4FS/\nZbypSeqjqr4GfH+N193q55Y3re1yzi3TCoJRvkD2Zn0uNfaGqjoDUFWngeunNN9x9FXbsPuBr0w8\n08vTS31JPgy8VlVHpj3hMfX183sn8LeTPJ/kuQ1cOumrvn8M/Msk3wF+AzgwxTmP43LqO7VGn9Wu\n36LnllFqGzbSuWUjv8V1OZ+D3Sp3tkeuLck/A85X1WM9zmfaLllfkquBTwCLo47ZZEaZ63bguqr6\nAPCrwKF+pzRVo9T3i8DHquonWAmFz/c7pQ23Vc4tIxvn3DKtIDgF/MTQ8xv5/y8lTwE/vkafS409\nffEyKMkO4LtTmu84+qqNJPcBPwv8/PSmO7Y+6vvLwM3AHyf5dtf+R0k24l1XXz+/k8B/BKiqF4EL\nSd4+vWmPrK/67q2q/wRQVV9kZRljI0xS36Wc2eLnlksa+9wypZseb+GNmx5zrNz02LWqz8/yxk2P\nD/DGDas3HcvKDZ2Pb/ANnb5quwP4FvD2K13Tlahv1fhvs/LueWbqAx4E/kW3/U7g1Rmp7+LN4m8B\nf6fb/mngxa1W39D+m4Ejq9q29LllndrGPrdMs6g7WPk28XFgf9f2IPAPh/r8Vlf4HwPvv9TYrv1t\nwDPdvqeBH92gv4x91HYceBX4evd4dCNq66u+Va//39mgTw31+PN7K/BvgSPAH9KdNGeovr/Z1fUS\n8AfAbVu0vseA/wGcA74D/IOufRbOLW9W29jnFr9QJkmN87/8lKTGGQSS1DiDQJIaZxBIUuMMAklq\nnEEgSY0zCCSpcQaBJDXu/wFIbVDBmwJqGgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x12b464350>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist([abs(x) for x in [weights[key] for key in weights] ],bins=100);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we see above, most of the interactions are negligible. Let's drop them out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 630,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "significant_weight_keys = [key for key in weights if abs(weights[key]) > 1e-4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 631,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "significant_weights = dict((key,weights[key]) for key in significant_weight_keys)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 632,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "significant_weights =  {k:v for k,v in weights.iteritems() if abs(v) > 1e-4} ## more elegant "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 633,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "76"
      ]
     },
     "execution_count": 633,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(significant_weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 782,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "H_nonlin_sp = 0  ## with sympy only\n",
    "\n",
    "for combination,pm_arr in significant_weights:\n",
    "    H_nonlin_sp += make_nonlin_term_sp(combination,pm_arr)*significant_weights[combination,pm_arr]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 807,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "H_sp = H_lin_sp + H_nonlin_sp*NONLIN_WEIGHT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 808,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_sp_conj(A):\n",
    "    '''\n",
    "    Returns the symbolic conjugate of A.\n",
    "    Args:\n",
    "        A (symbolic expression in symbols a[i] and a_H[i])\n",
    "    Returns:\n",
    "        The complex conjugate of A\n",
    "    '''\n",
    "    A_H = sp.conjugate(H_sp)\n",
    "    for i in range(len(a)):\n",
    "        A_H = A_H.subs(sp.conjugate(a[i]),a_H[i])\n",
    "        A_H = A_H.subs(sp.conjugate(a_H[i]),a[i])\n",
    "    return A_H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 809,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def make_eq_motion(H_sp):\n",
    "    '''\n",
    "    Input is a tuple or list, output is a matrix vector\n",
    "    '''\n",
    "    A_H = make_sp_conj(H_sp)\n",
    "    diff_ls = [1j*sp.diff(H_sp,var) for var in a_H] + [-1j*sp.diff(A_H,var) for var in a]\n",
    "    fs = [sp.lambdify( tuple(a+a_H),expression) for expression in diff_ls ]\n",
    "    return lambda arr: (np.asmatrix([ f(* arr ) for f in fs])).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 810,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "    A_H = sp.conjugate(H_sp)\n",
    "    for i in range(len(a)):\n",
    "        A_H = A_H.subs(sp.conjugate(a[i]),a_H[i])\n",
    "        A_H = A_H.subs(sp.conjugate(a_H[i]),a[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 811,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "eq_mot = make_eq_motion(H_sp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 812,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def double_up(M1,M2=None):\n",
    "    if M2 == None:\n",
    "        M2 = np.zeros_like(M1)\n",
    "    top = np.hstack([M1,M2])\n",
    "    bottom = np.hstack([np.conj(M2),np.conj(M1)])\n",
    "    return np.vstack([top,bottom])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 813,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A_d,C_d,D_d = map(double_up,(A,C,D))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 814,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "B_d = -double_up(C.H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 815,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def make_f(eq_mot,B,a_in):\n",
    "    '''\n",
    "    Nonlinear equations of motion\n",
    "    '''\n",
    "    return lambda t,a: np.asarray(eq_mot(a)+B*a_in(t)).T[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 816,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_f_lin(A,B,a_in):\n",
    "    '''\n",
    "    Linear equations of motion\n",
    "    '''\n",
    "    return lambda t,a: np.asarray(A*np.asmatrix(a).T+B*a_in(t)).T[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 817,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a_in = lambda t: np.asmatrix([1.]*4).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 818,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = make_f(eq_mot,B_d,a_in)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 819,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f_lin = make_f_lin(A_d,B_d,a_in)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 820,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.33972600-19.11352881j]\n",
      " [-0.87314689+23.3997759j ]\n",
      " [-0.42266661+11.76116304j]\n",
      " [-1.60513022 +2.79281965j]\n",
      " [-1.09514712 -7.76191194j]\n",
      " [-0.33972600+19.11352881j]\n",
      " [-0.87314689-23.3997759j ]\n",
      " [-0.42266661-11.76116304j]\n",
      " [-1.60513022 -2.79281965j]\n",
      " [-1.09514712 +7.76191194j]]\n"
     ]
    }
   ],
   "source": [
    "eq_res = eq_mot([1.]*10)\n",
    "print eq_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 821,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.33972600-20.178468j  ]\n",
      " [-0.87314689+20.38977829j]\n",
      " [-0.42266661+10.69622386j]\n",
      " [-1.60513022 +0.1043035j ]\n",
      " [-1.09514712-10.4504281j ]\n",
      " [-0.33972600+20.178468j  ]\n",
      " [-0.87314689-20.38977829j]\n",
      " [-0.42266661-10.69622386j]\n",
      " [-1.60513022 -0.1043035j ]\n",
      " [-1.09514712+10.4504281j ]]\n"
     ]
    }
   ],
   "source": [
    "mat_res = A_d*np.asmatrix([1.]*10).T\n",
    "print mat_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 822,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 7.18115794])"
      ]
     },
     "execution_count": 822,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## compute L2 error between \n",
    "np.sqrt(sum(np.asarray(abs(eq_res - mat_res))**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 845,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "r = ode(f).set_integrator('zvode', method='bdf')\n",
    "r_lin = ode(f_lin).set_integrator('zvode', method='bdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 846,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "y0 = np.asmatrix([0.]*10).T\n",
    "t0=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 847,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<scipy.integrate._ode.ode at 0x12d9a89d0>"
      ]
     },
     "execution_count": 847,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r.set_initial_value(y0, t0)\n",
    "r_lin.set_initial_value(y0, t0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 848,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t1 = 100\n",
    "dt = 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 849,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Y = []\n",
    "\n",
    "while r.successful() and r.t < t1:\n",
    "    r.integrate(r.t+dt)\n",
    "    u = a_in(r.t)\n",
    "    Y.append(C_d*r.y+D_d*u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 850,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Y_lin = []\n",
    "\n",
    "while r_lin.successful() and r_lin.t < t1:\n",
    "    r_lin.integrate(r_lin.t+dt)\n",
    "    u = a_in(r_lin.t)\n",
    "    Y_lin.append(C_d*r_lin.y+D_d*u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 851,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcXGWd7/HPr/aqru5OOitZCIEIAQyEKAGJSAOyKzCD\nC4LrOOgdh+vIjIh4x0v0unIZBZxhFOGKjjIwhi2OIIrQIEJYhCQEEkIgZOmQTkg6vVZ3bc/945xe\n0pB0J6mu0931ffvqV6pOnTrnV0de9a3nPM95jjnnEBGRyhQKugAREQmOQkBEpIIpBEREKphCQESk\ngikEREQqmEJARKSCDRoCZnarmTWZ2cpB1jvezHJm9telK09ERIbTUFoCPwPO2tsKZhYCvgc8WIqi\nRESkPAYNAefc40DzIKv9T2AJsK0URYmISHkccJ+AmU0DLnTO/TtgB16SiIiUSyk6hq8Hrur3XEEg\nIjJKREqwjXcDd5iZAROBc8ws55xbOnBFM9NERSIi+8E5Nyw/sIfaEjD28AvfOXeo/zcbr1/gC28X\nAP3W159zXHPNNYHXMFL+dCx0LHQs9v43nAZtCZjZ7UA9MMHMNgLXADHv+9zdPGB1/dIXERlFBg0B\n59wlQ92Yc+5vDqwcEREpJ10xHJD6+vqgSxgxdCz66Fj00bEoDxvu80277czMlXN/IiJjgZnhAu4Y\nFhGRMUghICJSwRQCIiIVTCEgIlLBFAIiIhVMISAiUsEUAiIiFUwhICJSwRQCIiIVTCEgIlLBFAIi\nIhVMISAiUsEUAiIiFawUt5eUAYpFRybTTTIZp7s7S6FQpL09Qy5XoFh0ZLN5zIxw2MjlChQKRcC7\n61qh4N1JqFgsUixCsVikZ+LVgXcZ6j8ha//lxeLbL9cMriIykELA19zcRnNzO9lsgVdffZNMJs+r\nr3ayc2eR5mbo7jY6OkI0NSUwc+zalSQez9PSkqKrK0Zraw2RSB6zIh0daXK5OIWCIxw2zIx4PEQk\nUsDMiES8GWGLRS8IQiHDu/2yEQo5zLypY73H3ro9r3uv9dXdd9tmo/+N3frWsd3W0W2eRaS/irif\nQLHoCIWMe+99jvr6w/nyl5/l0ksnsXhxB8ce28nttx9JNhujrW080WgXuVyCdLqFmTMbmTixncmT\nu5k8uUAkAnPnRonHQ4wbFwOgri5JbW2Curqq3l/vtbUpJk2qJZ8vEImEy/55RWRsGc77CYzZECgW\nHb/61dOceuqhzJw5ibvueo6LLlrAKac8xaOPngBAKFSgWAxz9dUN3H9/HTffnGT16p1cdNExZLM5\n6upqylKriMjeKAT2QbHo+OpXH+WccyZx2mlH8/GPP8ovf3kKtbU7aGmZAMCXvtTAHXfMYdmyGM8+\nu4mLLnrXsNYkInIgFAKD2LRpO9XVSb785Wd517uifOELizj22JWsWHEMAEccsZqXXz6Se+55jnXr\nOvjyl08ueQ0iIsNFIbAHt9zyJOPHx7j88ikcckgTy5a9iylTGmlqmg7AOec8wQMPnMTttz/DCSfM\n4NBDDyrZvkVEykUhMMDXvtZANGp885unMH36BhobZ+32+qmnPsUjj5zAt7/9KGecMZ0FCw4lHNYl\nESIyOg1nCIyqIaLr1r1BNBrmu9+tx8wbDN/enu59vecU0LvfnWHmzEe58MLZHHXUwUGVKyIy4o2q\nEDj//GYAzCbinPfLvqezF2DevGZWrIBx44xrrz0lkBpFREaTURUCa9e+g0IhSjTaTS4X711+8MHr\n2bhxNtOmOX75y6c57zyN9hERGYpREQKf/WwDbW0hwuHjKRSiRCI5crk4yWQbmUw1c+duZePG2UyY\nEOLSSxcGXa6IyKgxKkLg7rvfya5dE0km2wAIh/MATJ3axPr11aTTeR588AVOOun4IMsUERl1RkUI\nRKM5ACIR78s/n48CUFXVBcCUKQXOPHNeMMWJiIxiIzoEfvKTJ1iyJEQ4PAOAWCwL9IVAOFzkD39Y\nxYIFCwKrUURkNBvRg+d/+tMUDz10IuFwAaD331Sq3f83y/vf/07N8SMisp9GZEugqyvL44+vparK\n++UfDnvXBGQyScCbPvmGGx7nhBMmB1ajiMhYMCJbArfd9ixnnPFOYrHibsvb2sYDXgh88Yvv5YQT\nDg+iPBGRMWNEhcCGDds49NB1NDd7HcCRSE8LILHberoxiohIaYyo00H3PL6c9evPZE3TKgB2tHkh\n0Nzcd1Xwvfc+TzIZASYFUaKIyJgyIloC6zZuYdrMV9mR7wCgucX7pb+zxWsB5PMx4gnvtQsuOE7D\nQUVESmREhMD9K17gjc2H0ZXxRv9k2r0GSqa1b9RPuqY5kNpERMayEREC/lxw7Grx+gI62rz793a2\n94VAIu0NC53/0C3lLU5EZAwLNARyuTxfvHYJXV3el39Lq/dvpt2bHK69tY5UVQsAdUduYNJBG1gR\nmRNMsSIiY1CgIfDfy5/jR1d9iEy39+Xf1uotz7SlAK8vIFnlzRf04pcitN2yJpA6RUTGqkBDIB7x\npn9o6egGoLsjDECmo6p3nUTK6xAuRsJ0pbwWws9fXcb9javKWaqIyJg0aAiY2a1m1mRmK/fw+iVm\ntsL/e9zMhjx0pyvrTQzX1ur9293hhUKmo+9uYbGqzre87zOvbeNDLywf6m5ERGQPhtIS+Blw1l5e\nfw14n3PuWOBbwE+HuvOMHwLtrd71ANlOLwQ622sIh73XoqmMt7Ir9L7PRWvIRMYPdTciIrIHg4aA\nc+5xYI/jM51zy5xzLf7TZcD0oe68pyWQ6fDun5zt8K4LyHTWUFW9yyvQv3dAJLtzt/cau08pISIi\n+67UfQJ/Czww2EqZTDf1F93XGwJdfgjk/OkhisVwb19A6p2bqK7ZQcp5zy3vDRXFOZbv3MSG9jdL\n/BFERCpHyaaNMLNTgc8A793beosXL2bj9u08evck5syLA2eQzXgdwrnOVO960ZjXWbzpoi7aPraC\nqdkCrUCk0E4u4vUZLPjLc6QKLbSf/clSfQwRkcA1NDTQ0NBQln2VJATM7BjgZuBs59xeL+1dvHgx\ndz7+Z3520yLqZv0agGzGK6O7XwiE/RvIxMxBKETUnzQu4vLkABeKQThBBk0mJyJjS319PfX19b3P\nv/GNbwzbvoZ6Osj8v7e+YHYwcBfwCefcq0PZWFfe6+Rta/NOB+W7vA7h/iEQ8VsCCX+v0d5Xer70\nvRfM5Yf0AURE5K0GbQmY2e1APTDBzDYC1wAxwDnnbga+DtQBN5mZATnn3MK9bTPrh0B7h9e5m8v4\nIZCpIhQqUCyGicS9lkB0wHvDfghYsQsX9q4bWL5zE1OTNUxN1g76gUVEpM+gIeCcu2SQ1y8DLtuX\nnWZzPaOCvC/0fJc3V1BXZ5pkqpWO9vGEY7nd3vP9uQu4q/Fljpl0EHdvfZVnQ94gpEK4muOWv8yE\n7EbePPtv9qUMEZGKF8j9BLpzXgugq9MLgZ5RQd3dKarHvUlH+3hC/v2Ep8YirAE+NGsBH5rl3VD+\n6qOh+ne/IFNIUQglIRyna2TdGkFEZFQI5Jsz54dArjvk/xvvfS2e8q4Qjpy/lmOmdPKrhReytnXb\nW7bxyLtPoeCKnPjCKwC4t++yEBGRvQgkBPIF71d+rssPga6+20dGol5fwM5jImw6Ps201Dimpca9\nZRvvnjjLexDaAEAmXEvid7dz8bgIt534keEsX0RkzCj7BHKnXLyErqw3oifX7V8f0JXsfT0c9aeL\nGOL2/irSyJeqW3DRGroT03ikpbWk9YqIjGVlbwk8dueHmHvC7QDku7zdZ/uHQMQLiI9OOYgfb12P\nNzBpz+5+76UAXO9fWNHkUhz++1v4wVEn8IEZ+38byi2du7hjw/NkiwVWt+2kJZeltZAn6xxF53BA\nxIxkKEQqFCYViZAOR0lHIqQjMWojcaqjcaoiMRKhCMlIlEQoSiwUwsyIWJiwhQiZdxqr6Lz+EeeP\nfup5Xux3HUTReafRnBu4DnvdRv/t7G0dEak8gZwO2tXinQ7Kd/sXib1NCHxn/nl8Z182Wsgwt9DI\nmsQcXgG+tvoplm5Zy3Xzz6Umltzj2/LFAjeve4JsocCNG9exy4Vpjs0AV6A610SEIhNDBapCRnU4\nRNyMUCiEAXnnaCsU2J7L05XpIuug20EWI0eIHGGKhLw/C+EI46zvkguHgVnfpQ8DLnyzt70QbsA6\nb/kSH/o2bMBzEak8gYRApt370ilke0KgikgkSz4fIxTdv4u/mk5axMR4mvBjjxHK7eKFyGxe6Ayz\n/em7+VN7N795Vz1Xv/gYPz3uHP766d9w9oQpXL+9i5DLk43UgIX4cLqKSfEEXz96IYlwlHGx1OA7\nFhEZZsZeR+ofkGBGB2W936D5bu/6gFwuTlW6mXx7jFCksLe37tHkpHc/4mld6/jIpElc3xaGYpZ7\n89MhAYuWr8JFD2HuU09QjM1hVSuclWymvVDg/kUXkorEiITCpfmAIiKjRDAh4I8K6gkBgFg8Q0f7\neBIH7TigbTee/bcA/Pb3t/DRKdP4VnOM020Df7SJfDKxjTWFdm49diFhC3HkuPoD2peIyGgXzBDR\nrPeLu9AvBOIJ7+YxuUt2MP7DvwE+eED7WHumFwafbtvGrKqTCeGdyxcRkT7BhIB/kVg+2xcC0UQX\nAL84+QO80r69ZPs6rHpyybYlIjLWBBICPR3C+X5XCveMCjp83BSOmTAjiLJERCpOIOdHek4H5bJ9\nIdAzV1AiPNTLxERE5EAFEgI9LYH+00VYaP9GBYmIyP4LKAS8X/vZbIJo1Lt5jGn+NxGRsgs0BHLd\nCWJxb9bQYlEjd0REyi2YPgF/aGh3d5JY3GsJ1E7fGkQpIiIVLZiWQM4LgWIxQjTuXR9wwWU5bn/2\nySDKERGpWAGNDuobART1WwJ1VSk+9q73BFGOiEjFCqgl0D8EvIvEUonYnlYXEZFhElBLoO/6gEjc\nu5PYxJp0EKWIiFS0YEIgFyUU8q4QDse8EJgyriaIUkREKlpAIRAjFvNOA0X8PoEZEyYEUYqISEUL\nKATiRP0WQCjs3Rxxcro2iFJERCpaICGQy8aIRL2WQM/dEWORQOayExGpaMGEQL+WQDLdFUQJIiJC\nUENEC1EiMa8vYM57drF5y5tBlCEiUvECm7AnEs15/0aM6QdNDKoMEZGKFlwI+KeDIlFNHyoiEpTA\nQiDc0xIIKwRERIISXAjEvBCIxRQCIiJBCTAEvNNBibiGhoqIBCXA00HetBGJeDioEkREKl6AHcPe\n6aBkQjeWFxEJStlDIBz2vvxDUe/G8sm4QkBEJChlD4Fi0btILBLzTgdVJRUCIiJBKXsIOOdfHxD3\nWgLpRKLcJYiIiC+APoGeloAXAhNrdTMZEZGgBBYC8bR3OuidM2aWvwQREQECDIFkdZE/rV3D9PF1\n5S9BRESAIYSAmd1qZk1mtnIv69xoZq+Y2XIzm7/3LXohEAob733H3H0sV0RESmkoLYGfAWft6UUz\nOwc4zDn3DuDzwI/3vjk/BAK7QkFERHoM+lXsnHscaN7LKhcAv/DXfQqoNbMpe16957aS+1CliIgM\ni1L8Hp8ObOr3vNFftgfencRCmjdORCRwgXUMqyUgIhK8Ukzh2Qj0H+c5w1+2B78AnmbVIytpeO9M\n6uvrS1CCiMjY0dDQQENDQ1n2Zc65wVcyOwT4jXNu3tu8di7w986588zsROB659yJe9iOgzuAj/Kp\n7/+S277y8QMqXkSkEpgZzrlhOYk+aEvAzG4H6oEJZrYRuAaIAc45d7Nz7n4zO9fM1gEdwGf2vkXv\ndFBYnQIiIoEbNAScc5cMYZ3Lh77LniGiCgERkaAF1jEcVsewiEjgAgwBtQRERIIWWAhEY7pkWEQk\naMHNIhrTDeZFRIIWQAh400bE1RIQEQlcAN/E3s1kQqY+ARGRoAUQAkXAu/hBRESCFVhLQBEgIhK8\nAELAm6Zi+pTa8u9aRER2E0hL4JSaD/LZU+vLv2sREdlNICFQ1bmz/LsVEZG3KHsIGI6C3zksIiLB\nKnsIxMNG0Tm2rngDVxx8GmsRERk+AYRAiCJFps+fwQ2LFpd79yIi0k9ALYEiRYps2LSh3LsXEZF+\nyh8CkVBvn0A4pKkjRESCFEwIOC8EYvFEuXcvIiL9lD0EYtEQGefNJOo0dYSISKDK3xKIhmkjA0C+\nkC/37kVEpJ+yh0A0Cq10AJDL5cq9exER6afsIVA3s5Z2OgHIFxQCIiJBKn9LIBmms+d0UF6ng0RE\nglT2EAiFwPkziapPQEQkWGUPgXCkb0RQoagQEBEJUvlbApG+XRYKhXLvXkRE+il7CER2awkoBERE\nglT+qaTDCgERkZGi/KeDon27zDuFgIhIkMo/g5s/aVyKpFoCIiIBK/91AnHvdFCSRO9EciIiEozy\ndwxHe1oCCQo6HSQiEqjyh0AsDECKuEJARCRgZQ+B0484BICUJXTDeRGRgJU9BOZMnQpAghib8tt4\nd2x+uUsQERFf+TuGUzEAUqEE69jIX3IrWPPAS+UuQ0RECCAEEmnvlpLJcJwCXp/A8rufLncZIiJC\nIC2BKADJcN/9hXds3FruMkREhCBaAtXel38qmuxd1trcVu4yRESEQG4q47UEUvFU77KO9vZylyEi\nIgRynUAEgKpUundZLttd7jJERIQAQiDmjw5KV1X1Lst1KwRERIIwpBAws7PNbI2ZrTWzq97m9Roz\nW2pmy83sBTP79J62FU/FmWdHcNARs3uX5bLZ/aldREQO0KAhYGYh4F+Bs4CjgY+Z2dwBq/098KJz\nbj5wKvAvZhbZ0zZXFtdQO2U8AFGiutewiEhAhtISWAi84pzb4JzLAXcAFwxYxwHV/uNqYIdzbq/f\n7NGkd1ooSYJcIbdPRYuISGkMJQSmA5v6Pd/sL+vvX4GjzGwLsAL4h8E2Gk3GAUgQJ1/I09XaNaSC\nRUSkdPZ4ymYfnQU875w7zcwOA/5gZsc4594y9nPx4sUArHtyLQBJ4jS0/YVkbRLnXInKEREZvRoa\nGmhoaCjLvmywL14zOxFY7Jw723/+VcA5577fb53/Br7rnPuz//yPwFXOuWcHbMv17O/R6x+i/ooz\nmMts1rAevI2W7IOJiIwVZoZzzgZfc98N5XTQM8AcM5tlZjHgYmDpgHU2AO8HMLMpwOHAa3vbaCzl\nnQ5KWnwfSxYRkVIZ9HSQc65gZpcDv8cLjVudc6vN7PPey+5m4FvAbWa20n/bV5xzO/e23ViV3ydg\nMa9bGSjmi4Qi5b/tsYhIpRpSn4Bz7nfAEQOW/aTf4zfw+gWGLFHjzR2U6NcS2PHaDiYdPmlfNiMi\nIgcgsJ/dPSHQfzbR7a80lXQf2fYsxXyRdfe9SOOTG0u6bRGRsaBUo4P2WaLG+/JPhPu1BF7dXtJ9\nrJtyEi3jZ/OexiUANNt4Vs84k+xBszj61n9k0junlHR/IiKjTWAtgbg/pXQy1jeldPOmHSXZdsP7\nv8Wjf3U9R3X+hfc0LuHR477E+ug7eOnv/pXjN93F+DVPMmneVF6PvYNVtywj16mL1USkMgXWEkiO\n8/sE+oVA69a99iUPWf0fvw7AY3MvI3XpX3PSP55ONPVDZgPuRx/j2JDx5JV3kWvczpGfO5/CZa38\n+f1fZ9IlZ3D0ZxaWpAYRkdEgsBCIp73TQLF4X59A67ZmNjz8KrNOO2y/tvnyncuJVMWZQYxl8z7H\nsfd9k3Gzx++2joW8obbv+b8XAdB106dpeq6RxMevZPJnb+S5f5pH6+xjmPh3H+HwixcQS8f2q5Ye\nuc4cnds76G7pwhUdrugNhXJFB67v8W7L/TpD0TChSIhQJISFQ72PRURKJbAQCMfCAFRV9d1XoOvx\nP+N+/z1wzfu1zRkXL6KKTlZWL+KUlT8a0nsS4xLMOu0wZm25m7YtbeT+7RF4cjmJyz9L4bLXWFV1\nLJlUHZkJMyim0hCLQyQCmQyhznbCmXYiXe1Eu9uJZduJ59tJ5ttJFtupcu1EyANpzOJgBhgOwwCH\n99wArG85/lq4IlDs+5ci5go4G5ZrRkSkAgUWAj2/yCfPOAhe9pYt6XyNf2IX+3vdcIIutoYOIvSD\n6/br/dXTqjnh2+cD5wP/m45tHRT+6zncthZ4ZRO0dUAu5/3V1kJNGmrTMC4NddVQl4a6KphcDROr\nYEqacE2c2pC+tEXkAAzjD7/AQgDg9ede56mfPw5/hBAhGtm/FgBA25Y2wsSZ3L2ZqSU6ZVI1uYpj\nLz+5JNsSERmJAg2BWcfNYsXdzwFQTZos3igdV3S9LYWhaLjwh9DWzoz4bObonLmIyJAFGgIACX+o\naDVVtNMJQNeuLpJ1yb29bTfzl36TcW4Xq6pOGJYaRUTGqsBDIDnOu9dwDSnexBsi2tHUvk8hkCPG\nc9c+RNWsicNSo4jIWBV4CKTqUgCkLUWX824437q5mYlHDm0Ooa5dXVS7FuZfcaqGT4qI7KPAvzWr\n6rwholXhvl/+W55cx2Mf+/dB37trfTN/+V9382Z4qgJARGQ/BP7NmZ7oh0Ak1bts50/v4H13fGHQ\n975w9pUsuulSon4LQkRE9s3ICYFoXwhku7NDe+827741E4qlnXhORKRSBN4n0DOHUDpeBW3esq5c\nlibSDDbH57iORp684r+omTeLo4e3TBGRMSnwlkA0EQWgtqqmd9l/t21jKu17nN3zjWc289iRn2NS\nrpGj//EsTfomIrKfAg+BUMgrYcbcvknjthRaANj67Gb+9Klb3vKetd9ZwvvW/JQ0HdTMqHnL6yIi\nMjSBhwCAc45pRxzc+7wbr0/g1cX/wcm/uGy3dTM7M8SferSs9YmIjFUjIgQAqidWA5AiSbc/fQSd\n3hXE7VvbefKflrDypsdJTkgxYee6oMoUERlTRk4ITKoFoIZq2v2WQKHNOy30wrUP8J4ffJhdT7wE\nwNTu11lRczLro4cHU6yIyBgR+OigHuNn1gFQa9Vsc28CUNzhDf3M7/DCoLhxMwDVtPOOV39HJDFi\nyhcRGZVGTEtg3EzvDmA14TTNeF/6odZtABRe3wRA/LWXetdPTUwd8F2/REQq3YgJgbqDvZZATbRv\ntI91efcXiGxeD8DUbS+UvzARkTFsxIRAosabUnry+Km9y6KulSJQ/aYXArNza4MoTURkzBoxIRAK\nhbj63Ks5au5RvcuWMJ4wMLn9td5l6+JH0UU8gApFRMaeERMCAN/57XeYedQhvc+fpwuAg4pbKPil\n2pIlbL73L0GUJyIy5oy44TWz5h8CQJQorfRNG9EUnsa0wmYO+8CRAVUmIjL2jKiWAMBBcw8CYLyN\no9VvCQB0RMYFVZKIyJg14kJgymHe3KFTUlNp9e85XAC2HryQzeFZAVYmIjL2jLjTQTWTvSGi0+um\n80KHNyR0RWw+xz/9b1jIgixNRGTMGXEtgVAoROOaRmbPnN27bNIzS0mMSxCv0aggEZFSGnEhADDt\niGmkq9K9zzvbOgOsRkRk7BqRIQCQru4LgTt/cCdmOhUkIlJqIzYE5i+aD0CSJA/96aGAqxERGZvM\nOVe+nZm5oe6vu7Obr1/4dR566iHasm2s61pHOWsVERkpzAzn3LCcDhmxIdDjlCmn8MS2J8iTp7Ol\nk2RNcpiqExEZmYYzBEbs6aAeqXiKPHkAmt9oDrgaEZGxZeSHQCIFwDgbx643dgVcjYjI2DLiQ6C1\noxWASfFJ7HxjJ02vNgVckYjI2DHiQyBf9E4FpWNprrz8SqbOmTrIO0REZKiGFAJmdraZrTGztWZ2\n1R7WqTez581slZk9UqoCr7/tepb+cCnpRJplO5eVarMiIsIQRgeZWQhYC5wObAGeAS52zq3pt04t\n8ARwpnOu0cwmOuffLX73be3z6KAe5x1yHvdvuB9Ao4REpKIEPTpoIfCKc26Dcy4H3AFcMGCdS4C7\nnHONAG8XAAcqHvXmDaq1Wra9vq3UmxcRqUhDCYHpwKZ+zzf7y/o7HKgzs0fM7Bkz+0SpCuzR0tEC\nwPjoeC479zJ+e8NvS70LEZGKU6qppCPAAuA0oAp40syedM6tG7ji4sWLex/X19dTX18/pB3s7NgJ\nQG28lj80/oENV2/gvH8474ALFxEZaRoaGmhoaCjLvobSJ3AisNg5d7b//KuAc859v986VwEJ59w3\n/Oe3AA845+4asK397hO48bIbeeSRR9jZtpPHtj1GFVW0u/b92paIyGgSdJ/AM8AcM5tlZjHgYmDp\ngHXuA95rZmEzSwEnAKtLWegXf/pF7ll3D6GQV3I33Vx19lU8dddTpdyNiEhFGfR0kHOuYGaXA7/H\nC41bnXOrzezz3svuZufcGjN7EFiJdzfIm51zLw1HwZlsBoDaUC3XPngtq9as4rcXqX9ARGR/DKlP\nwDn3O+CIAct+MuD5dcB1pSvt7XV0dwDeFcQ7Mjvo6OoY7l2KiIxZI/6K4YH++Sv/zBWnXUE84g0Z\n3da2jQsPu5Al314ScGUiIqPPiJ9Kek/mpeexqmMVSZJkyHD+oedz36v3lWTbIiIjSdAdwyPShPQE\nAJIh78rhTHeGxtWNZLuyQZYlIjKqjNoQ+NGvfsTDtz3cGwIbdm5gxlEzuPLsKwOuTERk9CjVxWJl\nN+/0eQDMuWoOjU2NrM+sB2Dj5o1BliUiMqqM2pZAjx8v+TF/uv1PRPw8a9zRyEeO/Ah3fOOOgCsT\nERn5Rm3H8EDXfeo67lx6Jy/teolOOrn4qIv5zxf/c1j2JSJSThV9o/l9USwWSYfTZMhwypRTmDxu\nMp/4H5/gg1/64LDtU0RkuA1nCIzaPoG3EwqF2LZjG+ueWsdx5x4HTVC4qUD9x+upnlgddHkiIiPO\nqO8TGChdl2beGfM4acJJ3Pj5G7n7lbupmVTD9y75HtvW6z4EIiL9janTQQNlu7LEk3E+PPfD/HrN\nrwG497p7SVWnOONzZ5StDhGRA6E+gQOw9ZWt1E6pJVWb4vxDz2fpa94EqHf+nztZ9cwqFt+zuHdm\nUhGRkUghUALb1m+jdkoti6Yt4n0L3scPH/khAKcfdDobmjfwwGMP8MbaNzj50pMDqU9EZE8UAsPg\ntzf8loldrFPZAAAFpUlEQVQzJnLt1deydedWlu1YRpEiiyYuIp1Mc9FFF1E3pY6FH1jI+GnjSdel\ngy5ZRCqUQqAM2ne2s/zB5Tx4+4MUCgUefOJBWrtbWde1jrSlqQpVMTU1lZkTZpJKpJg8cTI1tTUc\n8o5DqJ1QSzwZJxqPUjOhhmgsSjwVJ14V7/03UZUgno6TSCeIxMbUoCwRGWYKgYAU8gW62rvY/vp2\n1q9YT9PrTWx5fQstu1po3dXKmzveZHvzdpo7mim6Itl8lkw+Q9EVyRVz5FyOXDFH3uW9x/7/DCNC\nhChRIhbBMMIWfuu/ZoQIYWYYhlnffwPGHh7b2/93sqd19rSdkWxPn1FkrFrduVohMJbks3n+8Ls/\nsHD+Qro6uigWihTyBfK5PMV8kULOe+yKznteKFAsFHvf74ru7R+7t19edPv23nJ7fvXzHHfkcUNa\nt1gsDr7SKLZ89XLmHzk/6DJGBB2LPqd+6lRdLDaWRGIRnnruKc45/5ygSxkRHl78MO/7+PuCLmNE\naFjcQP0n64MuY0TQsejnU8O3aY2NFBGpYAoBEZEKVvY+gbLtTERkDBkTHcMiIjKy6HSQiEgFUwiI\niFSwsoWAmZ1tZmvMbK2ZXVWu/ZaLmc0ws4fN7EUze8HMvugvH29mvzezl83sQTOr7feeq83sFTNb\nbWZn9lu+wMxW+sfq+iA+TymYWcjMnjOzpf7zijwWZlZrZr/2P9uLZnZCBR+LK8xslf85fmVmsUo5\nFmZ2q5k1mdnKfstK9tn9Y3mH/54nzezgIRXmnBv2P7ywWQfMAqLAcmBuOfZdrj9gKjDff5wGXgbm\nAt8HvuIvvwr4nv/4KOB5vGs1DvGPT08fzVPA8f7j+4Gzgv58+3lMrgB+CSz1n1fksQBuAz7jP44A\ntZV4LIBpwGtAzH9+J94I+Io4FsB7gfnAyn7LSvbZgb8DbvIffxS4Yyh1laslsBB4xTm3wTmXA+4A\nLijTvsvCObfVObfcf9wOrAZm4H3On/ur/Ry40H98Pt7/SXnn3OvAK8BCM5sKVDvnnvHX+0W/94wa\nZjYDOBe4pd/iijsWZlYDnOyc+xmA/xlbqMBj4QsDVWYWAZJAIxVyLJxzjwPNAxaX8rP339YS4PSh\n1FWuEJgObOr3fLO/bEwys0PwEn8ZMMU51wReUACT/dUGHpNGf9l0vOPTY7Qeqx8CVwL9h59V4rGY\nDbxpZj/zT43dbGYpKvBYOOe2AP8CbMT7XC3OuYeowGPRz+QSfvbe9zjnCsAuM6sbrAB1DJeYmaXx\nUvgf/BbBwDG4Y35MrpmdBzT5LaO9jW0e88cCrzm/APg359wCoAP4KpX538U4vF+rs/BODVWZ2aVU\n4LHYi1J+9iFdV1CuEGgE+ndSzPCXjSl+E3cJ8B/Oufv8xU1mNsV/fSrQc6PjRmBmv7f3HJM9LR9N\nFgHnm9lrwH8Cp5nZfwBbK/BYbAY2Oeee9Z/fhRcKlfjfxfuB15xzO/1fqvcAJ1GZx6JHKT9772tm\nFgZqnHM7ByugXCHwDDDHzGaZWQy4GFhapn2X0/8DXnLO3dBv2VLg0/7jTwH39Vt+sd+jPxuYAzzt\nNwlbzGyhmRnwyX7vGRWcc19zzh3snDsU7//rh51znwB+Q+UdiyZgk5kd7i86HXiRCvzvAu800Ilm\nlvA/w+nAS1TWsTB2/4Veys++lL6p5j4MPDykisrYM3423oiZV4CvBtE7P8yfbxFQwBv59DzwnP+Z\n64CH/M/+e2Bcv/dcjdfrvxo4s9/ydwEv+MfqhqA/2wEel1PoGx1UkccCOBbvh9By4G680UGVeiyu\n8T/XSrxOzGilHAvgdmAL0I0XiJ8BxpfqswNx4L/85cuAQ4ZSl6aNEBGpYOoYFhGpYAoBEZEKphAQ\nEalgCgERkQqmEBARqWAKARGRCqYQEBGpYAoBEZEK9v8Bix64/cuFqroAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1223e4b50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(4):\n",
    "    plt.plot([(y).real[i][0,0] for y in Y ])\n",
    "    \n",
    "for i in range(4):\n",
    "    plt.plot([(y).real[i][0,0] for y in Y_lin ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### When Making the nonlinear terms above zero, we find agreement with the linear equtions of motion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the symbolic packages is kind of slow. For classical simulations maybe we can avoid that. We just need to extract the equations of motion, which should end up being sparse in the interaction terms.\n",
    "\n",
    "TODO: implement without sympy, e.g. with Julia"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ##############################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Testing different cases with make_nonlinear_interaction\n",
    "\n",
    "making sure different exceptions get caught"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "roots_to_use = np.array([roots[i] for i in combination])\n",
    "modes_to_use = [modes[i] for i in combination]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def call_make_non_lin():\n",
    "    return functions.make_nonlinear_interaction(roots_to_use, modes_to_use, delays, delay_indices,\n",
    "                            start_nonlin,duration_nonlin,pm_arr,\n",
    "                            indices_of_refraction,\n",
    "                            eps=1e-12,func=lambda z : z.imag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0021779902481650101+0.0036737045407247908j)"
      ]
     },
     "execution_count": 208,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "call_make_non_lin()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0051419161593478663+0j)"
      ]
     },
     "execution_count": 209,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "indices_of_refraction = 1000.\n",
    "call_make_non_lin()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "Exception",
     "evalue": "each element of start_nonlin must be greater than 0.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mException\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-210-699e1b01d691>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mstart_nonlin\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m  \u001b[0;31m## this shouldn't happen\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mcall_make_non_lin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-207-65a81927434a>\u001b[0m in \u001b[0;36mcall_make_non_lin\u001b[0;34m()\u001b[0m\n\u001b[1;32m      3\u001b[0m                             \u001b[0mstart_nonlin\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mduration_nonlin\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpm_arr\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m                             \u001b[0mindices_of_refraction\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m                             eps=1e-12,func=lambda z : z.imag)\n\u001b[0m",
      "\u001b[0;32m/Users/gil/Google Drive/Copy/Potapov/Potapov_4/Potapov_Code/functions.py\u001b[0m in \u001b[0;36mmake_nonlinear_interaction\u001b[0;34m(roots, modes, delays, delay_indices, start_nonlin, duration_nonlin, plus_or_minus_arr, indices_of_refraction, eps, func)\u001b[0m\n\u001b[1;32m    281\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mdelay_index\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mstart_loc\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mzip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdelay_indices\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mstart_nonlin\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    282\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mstart_loc\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 283\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'each element of start_nonlin must be greater than 0.'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    284\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mduration_nonlin\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mstart_loc\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0mdelays\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mdelay_index\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    285\u001b[0m             raise Exception('duration_nonlin + start_loc must be less than the '\n",
      "\u001b[0;31mException\u001b[0m: each element of start_nonlin must be greater than 0."
     ]
    }
   ],
   "source": [
    "start_nonlin = -1  ## this shouldn't happen\n",
    "call_make_non_lin()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "Exception",
     "evalue": "duration_nonlin + start_loc must be less than the delay of index delay_index for start_loc in start_nonlin and delay_index in delay_indices.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mException\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-211-4b091efebb23>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mstart_nonlin\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mroots_to_use\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mcall_make_non_lin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-207-65a81927434a>\u001b[0m in \u001b[0;36mcall_make_non_lin\u001b[0;34m()\u001b[0m\n\u001b[1;32m      3\u001b[0m                             \u001b[0mstart_nonlin\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mduration_nonlin\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpm_arr\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m                             \u001b[0mindices_of_refraction\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m                             eps=1e-12,func=lambda z : z.imag)\n\u001b[0m",
      "\u001b[0;32m/Users/gil/Google Drive/Copy/Potapov/Potapov_4/Potapov_Code/functions.py\u001b[0m in \u001b[0;36mmake_nonlinear_interaction\u001b[0;34m(roots, modes, delays, delay_indices, start_nonlin, duration_nonlin, plus_or_minus_arr, indices_of_refraction, eps, func)\u001b[0m\n\u001b[1;32m    285\u001b[0m             raise Exception('duration_nonlin + start_loc must be less than the '\n\u001b[1;32m    286\u001b[0m                            \u001b[0;34m+\u001b[0m\u001b[0;34m'delay of index delay_index for start_loc in '\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 287\u001b[0;31m                            +'start_nonlin and delay_index in delay_indices.')\n\u001b[0m\u001b[1;32m    288\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    289\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mindices_of_refraction\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mException\u001b[0m: duration_nonlin + start_loc must be less than the delay of index delay_index for start_loc in start_nonlin and delay_index in delay_indices."
     ]
    }
   ],
   "source": [
    "start_nonlin = [1]*len(roots_to_use)\n",
    "call_make_non_lin()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0050904969977543877+0j)"
      ]
     },
     "execution_count": 212,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "start_nonlin = 0.00001\n",
    "duration_nonlin = .099\n",
    "call_make_non_lin()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0050904969977543877+0j)"
      ]
     },
     "execution_count": 213,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "start_nonlin = [0.00001]*len(roots_to_use)\n",
    "duration_nonlin = .099\n",
    "call_make_non_lin()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Unused methods below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\nconsolidated_weightes = {}\\nfor key in significant_weights:\\n    if not key[0] in consolidated_weightes:\\n        consolidated_weightes[key[0]] = significant_weights[key]\\n    else:\\n        consolidated_weightes[key[0]] += significant_weights[key]\\n'"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## consolidated weights do not take into account which modes are createad or annihilated.\n",
    "\n",
    "consolidated_weightes = {}\n",
    "for key in significant_weights:\n",
    "    if not key[0] in consolidated_weightes:\n",
    "        consolidated_weightes[key[0]] = significant_weights[key]\n",
    "    else:\n",
    "        consolidated_weightes[key[0]] += significant_weights[key]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## QNET annihilation and creation operators\n",
    "\n",
    "a_ = [Destroy(local_space('fock', namespace = str(i))) for i in range(m)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Make linear Hamiltonian with QNET\n",
    "\n",
    "H_lin = sum([a_[i].dag()*a_[i]*Omega[i,i] for i in range(m)])  ## with QNET\n",
    "def make_nonlin_term(combination,pm_arr):\n",
    "    '''\n",
    "    Make symbolic term\n",
    "    With QNET\n",
    "    '''\n",
    "    r = 1\n",
    "    for index,sign in zip(combination,pm_arr):\n",
    "        if sign == 1:\n",
    "            r*= a_[index].dag()\n",
    "        else:\n",
    "            r *= a_[index]\n",
    "    return r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Make nonlinear Hamiltonian in QNET \n",
    "\n",
    "H_nonlin = 0  ## with QNET\n",
    "\n",
    "for combination,pm_arr in significant_weights:\n",
    "    H_nonlin += make_nonlin_term(combination,pm_arr)*significant_weights[combination,pm_arr]\n",
    "    \n",
    "H_qnet = H_lin+H_nonlin"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}