2017-03-17_mg_bs_feedback_setup.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from collections import OrderedDict\n",
    "from copy import copy\n",
    "import single_sided_network\n",
    "import qdyn_model\n",
    "from notebook_plots import display_hamiltonian\n",
    "from sympy import Symbol\n",
    "import QDYN\n",
    "from QDYN.pulse import blackman\n",
    "from QDYN.shutil import mkdir\n",
    "import matplotlib.pylab as plt\n",
    "import numpy as np\n",
    "import sympy\n",
    "from sympy import pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "circuit = single_sided_network.network_circuit(n_nodes=2, topology='bs_fb')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAABZCAIAAACR5Js8AAAJNmlDQ1BkZWZhdWx0X3JnYi5pY2MA\nAHiclZFnUJSHFobP933bCwvssnRYepMqZQHpvUmvogJL7yxLEbEhYgQiiog0RZCggAGjUiRWRLEQ\nFBSxoFkkCCgxGEVUUPLDOxPn3vHHfX49884755yZA0ARBQBARQFSUgV8Pxd7TkhoGAe+IZKXmW7n\n4+MJ3+X9KCAAAPdWfb/zXSjRMZk8AFgGgHxeOl8AgOQCgGaOIF0AgBwFAFZUUroAADkLACx+SGgY\nAHIDAFhxX30cAFhRX30eAFj8AD8HABQHQKLFfeNR3/h/9gIAKNvxBQmxMbkc/7RYQU4kP4aT6edi\nz3FzcOD48NNiE5Jjvjn4/yp/B0FMrgAAwCEtfRM/IS5ewPmfoUYGhobw7y/e+gICAAh78L//AwDf\n9NIaAbgLANi+f7OoaoDuXQBSj//NVI8CMAoBuu7wsvjZXzMcAAAeKMAAFkiDAqiAJuiCEZiBJdiC\nE7iDNwRAKGwAHsRDCvAhB/JhBxRBCeyDg1AD9dAELdAOp6EbzsMVuA634S6MwhMQwhS8gnl4D0sI\nghAROsJEpBFFRA3RQYwQLmKNOCGeiB8SikQgcUgqkoXkIzuREqQcqUEakBbkF+QccgW5iQwjj5AJ\nZBb5G/mEYigNZaHyqDqqj3JRO9QDDUDXo3FoBpqHFqJ70Sq0ET2JdqFX0NvoKCpEX6ELGGBUjI0p\nYboYF3PAvLEwLBbjY1uxYqwSa8TasV5sALuHCbE57COOgGPiODhdnCXOFReI4+EycFtxpbga3Alc\nF64fdw83gZvHfcHT8XJ4HbwF3g0fgo/D5+CL8JX4Znwn/hp+FD+Ff08gENgEDYIZwZUQSkgkbCaU\nEg4TOgiXCcOEScICkUiUJuoQrYjexEiigFhErCaeJF4ijhCniB9IVJIiyYjkTAojpZIKSJWkVtJF\n0ghpmrREFiWrkS3I3uRo8iZyGbmJ3Eu+Q54iL1HEKBoUK0oAJZGyg1JFaadco4xT3lKpVGWqOdWX\nmkDdTq2inqLeoE5QP9LEado0B1o4LYu2l3acdpn2iPaWTqer023pYXQBfS+9hX6V/oz+QYQpoifi\nJhItsk2kVqRLZETkNYPMUGPYMTYw8hiVjDOMO4w5UbKouqiDaKToVtFa0XOiY6ILYkwxQzFvsRSx\nUrFWsZtiM+JEcXVxJ/Fo8ULxY+JXxSeZGFOF6cDkMXcym5jXmFMsAkuD5cZKZJWwfmYNseYlxCWM\nJYIkciVqJS5ICNkYW53txk5ml7FPsx+wP0nKS9pJxkjukWyXHJFclJKVspWKkSqW6pAalfokzZF2\nkk6S3i/dLf1UBiejLeMrkyNzROaazJwsS9ZSlidbLHta9rEcKqct5ye3We6Y3KDcgryCvIt8uny1\n/FX5OQW2gq1CokKFwkWFWUWmorVigmKF4iXFlxwJjh0nmVPF6efMK8kpuSplKTUoDSktKWsoByoX\nKHcoP1WhqHBVYlUqVPpU5lUVVb1U81XbVB+rkdW4avFqh9QG1BbVNdSD1Xerd6vPaEhpuGnkabRp\njGvSNW00MzQbNe9rEbS4Wklah7XuaqPaJtrx2rXad3RQHVOdBJ3DOsOr8KvMV6Wualw1pkvTtdPN\n1m3TndBj63nqFeh1673WV9UP09+vP6D/xcDEINmgyeCJobihu2GBYa/h30baRjyjWqP7q+mrnVdv\nW92z+o2xjnGM8RHjhyZMEy+T3SZ9Jp9NzUz5pu2ms2aqZhFmdWZjXBbXh1vKvWGON7c332Z+3vyj\nhamFwOK0xV+WupZJlq2WM2s01sSsaVozaaVsFWnVYCW05lhHWB+1Ftoo2UTaNNo8t1WxjbZttp22\n07JLtDtp99rewJ5v32m/6GDhsMXhsiPm6OJY7DjkJO4U6FTj9MxZ2TnOuc153sXEZbPLZVe8q4fr\nftcxN3k3nluL27y7mfsW934Pmoe/R43Hc09tT75nrxfq5e51wGt8rdra1LXd3uDt5n3A+6mPhk+G\nz6++BF8f31rfF36Gfvl+A/5M/43+rf7vA+wDygKeBGoGZgX2BTGCwoNaghaDHYPLg4Uh+iFbQm6H\nyoQmhPaEEcOCwprDFtY5rTu4bircJLwo/MF6jfW5629ukNmQvOHCRsbGyI1nIvARwRGtEcuR3pGN\nkQtRblF1UfM8B94h3qto2+iK6NkYq5jymOlYq9jy2Jk4q7gDcbPxNvGV8XMJDgk1CW8SXRPrExeT\nvJOOJ60kByd3pJBSIlLOpYqnJqX2pymk5aYNp+ukF6ULMywyDmbM8z34zZlI5vrMHgFLkC4YzNLM\n2pU1kW2dXZv9ISco50yuWG5q7uAm7U17Nk3nOef9tBm3mbe5L18pf0f+xBa7LQ1bka1RW/u2qWwr\n3Da13WX7iR2UHUk7fiswKCgveLczeGdvoXzh9sLJXS672opEivhFY7std9f/gPsh4YehPav3VO/5\nUhxdfKvEoKSyZLmUV3rrR8Mfq35c2Ru7d6jMtOzIPsK+1H0P9tvsP1EuVp5XPnnA60BXBaeiuOLd\nwY0Hb1YaV9YfohzKOiSs8qzqqVat3le9XBNfM1prX9tRJ1e3p27xcPThkSO2R9rr5etL6j8dTTj6\nsMGloatRvbHyGOFY9rEXTUFNAz9xf2pplmkuaf58PPW48ITfif4Ws5aWVrnWsja0Latt9mT4ybs/\nO/7c067b3tDB7ig5BaeyTr38JeKXB6c9Tved4Z5pP6t2tq6T2VnchXRt6prvju8W9oT2DJ9zP9fX\na9nb+aver8fPK52vvSBxoewi5WLhxZVLeZcWLqdfnrsSd2Wyb2Pfk6shV+/3+/YPXfO4duO68/Wr\nA3YDl25Y3Th/0+LmuVvcW923TW93DZoMdv5m8lvnkOlQ1x2zOz13ze/2Dq8ZvjhiM3LlnuO96/fd\n7t8eXTs6/CDwwcOx8DHhw+iHM4+SH715nP146cn2cfx48VPRp5XP5J41/q71e4fQVHhhwnFi8Ln/\n8yeTvMlXf2T+sTxV+IL+onJacbplxmjm/Kzz7N2X615OvUp/tTRX9KfYn3WvNV+f/cv2r8H5kPmp\nN/w3K3+XvpV+e/yd8bu+BZ+FZ+9T3i8tFn+Q/nDiI/fjwKfgT9NLOcvE5arPWp97v3h8GV9JWVn5\nBy6ikLxSF1/9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAAHXRFWHRTb2Z0d2FyZQBHUEwgR2hvc3Rz\nY3JpcHQgOS4xOJQFEHMAAA8xSURBVHic7Z1fiB1XHcd/qwGpkJBzMQ81pVnOpgFrIKmzlIKoQeYS\nayiIeq9KSyBE5uqjFHuHVmlRLDNIHyr45x4ofajSegdFKrHKHqR5EEPt7CZQG4jdwyYk7kPMzppU\nQtNAfDjJZHbv3Ls7Z+/Ozrn3+3kI587OzP6yd+Y7Z77nd35n4tatW1SQ2Z9KIuINZ+ck0+3PfN9N\nt3/m++7yQqKiuF/7nk+yTz3uZNsXTqr/vKWy7U88zO//As+2z/4mvv7vpF87G0y/9ppBriewokHq\nALYqyPs/VMQ5cU5KkbrTjmNKkhVtxyHGVrSlJKK+bdclorvtJKE4XtFmjBxnRVsH0K+9SUH2BlY0\nyPUEVjTINf96Qw3y6pvxe/9Irt/LHzzCmYpn/5pcfcA59BVG8d32tT/If/2Lsu1bX3Qdh+62ebJ4\nIl5czG9fuMZ2HXY4S869FifEdh12OKlzf1HZ9uXt/MEjnCWZtopn/7oisOv38s8+wQcHVjTInV9z\nOcsE2RNY0SDX89crFGT611s8Ed/7hFtUjT9S9AAiev/lyOAoUDKvvbbVEYAKcPYsnTu31UGAdXDh\nBRNdnTDog5PvUxAY/DJQJkKQ5211EACA9WF2w5oo+OKvpUFvHwBQPklCRMTYVscB1kSbYAXZZvCL\nXnzH7e2BK6WEEAZnAxvH8zze891Ledtx3SC+7w/hLKA4nHNvw69Rqa8OKs7CM2Ly1cLehomCP3xG\nEK2+sJRSSZI0Gg2DE4KNEEWRUmrzFFxKGcA02wrCMNy4ggNbWFigyeJHmSj4Vw+o3O2MMRfP+tKR\nOquhh+IvZPnga90qwjDc+Enw1dnCuWPBoeJHmeSiLO7HRWEB6L0BIkqS21Y4qDiem98zHoyJgr/4\nDhTcAvp0zcF4Ece3rXBQcRaeMRlHNFHwh89gxNICoOAAWMTCgslRw/TBQaWAAQoIl4E93PyxSb6A\nSR/8wqNwWC0Aty4g+OD24PKyfPBfvDGkLIdNQCnl+/709HS9Xk83+r5fq9WazaZSa/+NfN+fmJjY\nzBhLAi4KIPjg9lCeD77/7z9IqvpY55wHQeC6LmMsnYoSBIHned1utzdpupeRyX0uquBKqX6JicBe\nGMOETDuYXTbpGRsp+MfjWq3WarWqfMN3Op0oiqoc4WZT1EXhnLdarampqTAMK/uEBkVxHHKcrQ4C\nrIMdT5q40yYK/r/Hf845F0LU6/Xp6WkhRAVveMZYp9NpNpurYkuSJAzDKIqiKMrOmPB9PwxDKWUU\n3a0QFsex7/tSylarFdv2Lmrgg7fbbW1DTU1N2fhfBr3oorOg+pj54INyUY4dO3bx4sXe7e+9d/ye\nez6m23Ect1ot3/cPHz68Z88egwg2D9d1Pc9rNpszMzPpxnq9ntopQgjf94Mg0H5Lu91edYZms9nt\ndh3H0Z5M9jyVot1u12q1VRsvXfry7t1/6t35vvvuy/1aNbt27bp8+XKSJEIIIYTruktLS0MOF6yP\nS5cuZYdz1iT3m11aeoiIarW5YUYGNoFvnL7124M5I3D9bljHcYIgGKTgFy9ezNWs04/6v7z/c2fP\nntUfOeeNRuORRx45deqUUeSbSBAE09PTaV87juM4jlM3nHOuFVwI0e12Vx0rpdTFXrQVo//VLxyO\n48Rx7LquU40X1DAMeye++z4FwfcKnUcpNTU1lX70PK/RaAxlbjcwYPfu3RvvNOj3qGpcp2AQv39U\nzLxR2EgxyQefnPyv+JUgosYdqH91ji2n2+1OT09rget1ewb4P4wxIkqVUZfhDcNwfn5eb9c99E0K\ne+MYuCharPUjOa13CAW3Gmi3LZTngz97aU+73Z6fn+92u9UvRqizU7S7rXvN6cNGSqmdE8/zUs83\n1XTHcVbtTERavvXH9WS2bCFFFVwppZTqdDrz8/NBEFT8fwfWCXxwW3DJpBNs0ge/50G/shl3eiBO\nq22aF+h5XirEMzMz6Q7pPtoKz75GaHcl3Zkxln1WxXFc/ZqrRZf84JxX1ugHxmj5xuO4+px+QR4s\n/uJsouDfPOMTVVS8OOe5zka6Ueeo9O6QleN03aLcnYUQjDHdr6/yKwh6XoCwOo89KOIHix9louAH\nDxgcNCLoLHidTFnxqtmoLgsIPrg9HPy5yR1rouDSDSotXZtJbv+9muDFGRBcFHvgShIvrKwmI5lV\nzToBK8CqpYAwkmkPp18wEVYTBf/OeSx9awG4bwERcY4OuB28vdPE2jDKB6/W1EuQD3xwQPBP7OGL\nPzFRcCMXxa1oIgrIglsXEFwUe+CqrHzw2d8p182RhyiKUAupfJRSuVkxvk/VTlgHZYCRTFt484fy\nUDn54M1l0ZsP7rpuOl8RAFARoN228N6ke6j4UfDBRxb44ICg4PZQng8uOF7OLQC3LiCssmYPXJrk\n/5r0wZfnFBHkoeqMmw/u+34URZzztORvHMdBEGQrAKfLD8VxnCRJu92ucl2EoVC9xVdAPm++rA4V\nf28emg8OwNYSBIEuNpmO6+p1l9KSOL7v64I2+mO9Xh95+Sa8itkDfHCwgnHzwZMkWZWWE8dxtgMu\nhMiuN+SNxx8ICm4L8MHBCsbt1o2iKJVvpVSr1XIcJ7tyHmMsuzz3OHTACT64PZj54CYKvjyHGQIW\n4I9Z7QO9fp6UUkoZxzFjrN1us0xxVW2ntFqtiYkJf2z+OkkCK9wO3nx52Csd98Mj+OCgckgp9bLU\n6ZZ6vf7222/rtl7dVJeWjKKo2Wymy8iNNqguawtmPrhJH3znQ6N/3Y8AY5WIonNLsvKtbXHdVkpl\nZwtr/4SNx9oHjGGRBzv49qul+eA0FkNAwCKklNkxzCRJwjBM005W1XsIw3CVwTLCwAe3BqN60MgH\nH1nGJB9cKSWE0GuZpu62zgRPxyp1Drj+aZIk2ZzCkQcmuC2Ulw8OHxxUB875mnI8PnrdC3xwW3jr\noHeo+FEmCg4f3ArGWLXAXcbDKxoFnvqlia6a+OC+gg9uAagKDYhISiyLaAeXny/LB997XhKN7VrH\n1iAEuuGgL77v68GDrQ5k7FBKdTqd3pr+//yjOvR04bOZKPjXd0LBAbCDAT54EAS5a4OATaXfbDL4\n4GAF6IADgg9uD/DBwQrggwOCD24Py981qfRgouB7z+OKsACj+QEAgK3h9GmTo+CDjyxjUPMDrA3y\nwW1h8TETb8PIB/8a5NsCxqP8NVgD+OC28K2nS/PBJRTcAuCDA6q8Dy6l9H1/amoqDEO9RSnVbDZr\ntdo6KwD7vj8xMbGZMZYEfHCwgqI+eIIKGiPK++9f2+oQ+uK6rq5gE8exrj7GOe92u57nrbMWwsiU\nTDDzwU0U/Nt7oOAWYOCDT09P+74PKR8lXJdOnfpJdnGiahIEwfgsu5HLuWMmjyITBSfMArCBoj44\nY8x13TAMa7Va9W94sH48zxNC1Ot1bVZU8wnNOXddt1fEdaHgKIqiKEqdFrqzjLWUMoqidGMcx3qu\naavVim0rqms2cLWN7lTH7/3Z8vJy7m380kt7jxNu76ozN1d76KGl3u3btm27efNm7iH79+/XDSGE\nEMJxHM/z+u0MNpvr168Xeo7mfrNzczUiOnLkyIkTJ5RSWvgajcaNGzeGGeswaLfb9Xp9Van3er3e\n7Xb1akpCCN/30956dhFUTbPZ1Os06VLDMzMzZca/fubm5no3bn/+N9eefrx3O2PM6Z9RdFvBVd6w\nVz8FP/DWn+UekyQWUCYnT37pypWcr2/Hjh1Xr17NPeSDDz7Yvn37tWu3bVNtTX744YebGCXoT1EF\nz/1mz5/fS0RXrlxJtyRJIqXct2/fUIIcLp1OJ7s2nr4C08XwOOdawYUQetXTLFJKpZT+3+mPRKSU\nUkqFYVgpNZ+dnc1+I5pPn7vxz7yv+/XXX3/33Xf7nWobEXl9uu+6TH7eT8Zj7QDLEYK8goUWfN/X\n8s059zzP8zzGWL1e35T4wFrUarWhDNMlSTI19ZRue57XaDRy/YoqoC+8VqulVbvXGxhgAekqXWn/\n/datW/qEnPOs/VIFjh8/3luRRnAK8pR4sB1k4oMv7ocPbgEGtloURY1Go9vtzs/Pj886ZKNNktCP\nfvQzxli73V5aWsqtilcp2u12KtOu6zqOk76ISCm1c+J5Xqpr6c6O46zaudS4N4znlrVW/YvvuOiB\nVx8piw05K6VmZmbGYfn2sSKOad++b83PP7vVgeQjpdSL5Onet97Y6XTS8cmZmRk9OKk/6pcSbYVn\nNVq7K+nOjLF0jT1bWHhGTL5aWFlNFPzhM4Kw2HHlKarg0O5R5YEHqmh5a1zX7X0n4Jyno5SMsU6n\n03tg1l/ShsmAna1gYYEmix9l4qJ89QBm+1kABBkQkeuOdfavUkpKmQ5vVhmzfHCTPvjifvdeg8NA\nuaAuCqA7a9WP7YiGHslMk1uqjOcqosLdLpM++IvvjPEz3R4q3+cAZRDHZNvUljFl4RmTetAmCv7l\nS6g8bQFQcEBEjI1vB9wuFhZMjjJxUT6/Gz64BYyz+wlSUB/cFm7+uKy6KBcehcNqAVBwQERJQpWs\ngwJW4/Ky8sF/8QYPvmBwHCiVotmE/dDFNIZwIrAVaBM890qIoqj6GRqjRxzHuZOqyssH/+YZnwhz\neqrOsBTc3gRb2xnKnJR+JrjnebnVkMBmoyea9m6fXeaTxc9mouAHDxgcBMpmWC5KxSdhg8H088F1\nml25sYBB7HjSxJ028cH/th8+uAVAeAERKYX19uzANSrZ/dHnnnuu389eeeWVo0eP9m4//iQ7epSU\nojC8/aqu20qR46xoS0lC9G0TEecr2kKQLoeQtjknxla0fZ+k7NvWspW2xyrIVYEJQY89ZnBJgJEi\njm9fCWaXGa3j+qeC98Ka17+lQW5QSS4/L/Z+J6fb1U+HNRNpSYFe6vV6flFdH9VlAQBgqAiRO5G6\nrw4TkeEqawAAAIaLUR0MIwVHxQ0AABguRpmdg3JRkAgMAAAlcfJkrogPTvocpOB9E4FfeomOHy8S\nGgAAgIHcuJGbQDY4nXfQSGZfMJIJAAAVwEjBAQAADBejWdTIRQEAgApgNJJppOAY3gQAgOFiNIsa\nfXAAAKgARgoOHxwAACpAeT44KuUAAMBwKc8HF1gnEwAAhgp8cAAAsBX44AAAYCt9ahMOZlB98EG/\nCRWCSwsSBdERZDZI3Cy2B9nvKzaq6I8+OAAA2Mr/AXQgzXV9XM6HAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "circuit.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "slh = single_sided_network.network_slh(n_cavity=2, n_nodes=2, topology='bs_fb')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\begin{pmatrix}1\\end{pmatrix}$"
      ],
      "text/plain": [
       "[[1]]"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "slh.S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\begin{pmatrix}- \\frac{\\sqrt{2} \\sqrt{\\kappa} \\cos{\\left (\\theta \\right )}}{\\sin{\\left (\\theta \\right )} - 1} \\hat{a}^{(c_{1})} - \\frac{\\sqrt{2} \\sqrt{\\kappa} \\cos{\\left (\\theta \\right )}}{\\sin{\\left (\\theta \\right )} - 1} \\hat{a}^{(c_{2})}\\end{pmatrix}$"
      ],
      "text/plain": [
       "[[-sqrt(2)*sqrt(κ)*cos(θ)/(sin(θ) - 1) * â^(c₁) - sqrt(2)*sqrt(κ)*cos(θ)/(sin(θ) - 1) * â^(c₂)]]"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "slh.L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_bs_decay():\n",
    "    from qnet.algebra.pattern_matching import pattern, wc\n",
    "    from qnet.algebra.operator_algebra import ScalarTimesOperator\n",
    "    θ = Symbol('theta', real=True)\n",
    "    κ = Symbol('kappa', positive=True)\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "    theta = np.linspace(0, 2*np.pi, 100)\n",
    "    L = slh.L[0, 0]\n",
    "    k1_form, k2_form = [op.coeff for op in pattern(ScalarTimesOperator).findall(L)]\n",
    "    k1_form = k1_form / sympy.sqrt(2*κ)\n",
    "    k2_form = k2_form / sympy.sqrt(2*κ)\n",
    "    k1 = np.array([k1_form.subs({θ: val}) for val in theta])\n",
    "    k2 = np.array([k2_form.subs({θ: val}) for val in theta])\n",
    "    ax.plot(theta/np.pi, k1)\n",
    "    ax.plot(theta/np.pi, k2)\n",
    "    ax.set_ylim(-3, 3)\n",
    "    ax.set_xlabel(r'BS mixing angle $\\theta$ ($\\pi$ rad)')\n",
    "    ax.set_ylabel(r'decay rate relative to open network')\n",
    "    plt.show(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAF9CAYAAACkiQodAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XmcjXX/x/HXZ8aMYTCjkDa39lt1Z0uljdJKSZtMi1KJ\nFpUWEi20EgmltNjCSKKUtSiVFkUlkbpVaJFdWWfM+f7+OMf9Gxpm5prrnGvOmffz8TgP53zPdc71\n1umYz3yv72LOOURERERiJSnoACIiIlK2qPgQERGRmFLxISIiIjGl4kNERERiSsWHiIiIxJSKDxER\nEYkpFR8iIiISUyo+REREJKZUfIiIiEhMqfgQERGRmIq74sPMOprZN2a2MXL7xMzOCzqXiIiIFI3F\n294uZtYCyAN+BAy4DrgXqOecWxxgNBERESmCuCs+CmJma4F7nHPDgs4iIiIie1cu6AAlYWZJQGug\nIvBpwHFERESkCOKy+DCzYwkXG2nA38DFzrnvg00lIiIiRRGXl13MrBxQC8gALgPaA6cXVICY2b7A\nucAvwLYYxhQREYl3aUBtYLpzbq1fbxqXxcfuzOxd4L/OuZsLeO5KYHTsU4mIiCSMq5xzY/x6s7i8\n7FKAJKD8Hp77BWDUqFHUqVMnZoFKq1P6XUMGBzLl7ieDjlIsjfq34LgKLXil4y107tyZ/v37Bx1J\nfKTPNLHo80wcixcv5uqrr4bIz1K/xF3xYWaPA1OB5UBl4CqgCXDOHl6yDaBOnTo0aNAgJhlLq4U/\n/8m2IxdxTbWucfffwvYrT6X0/WnQoAEZGRlxl1/2Tp9pYtHnmZB8HbYQd8UHUAMYAewPbAQWAOc4\n52YFmioODJgyDZzR+cLzg44iIiJlWNwVH865G4POEK+mLZ1MujuBOrWqBx1FRETKsLhbXl282bIt\nl1/TpnNCZvOgo4iISBmn4qOMeGn6J1D+L64/tUXQUUosKysr6AjiM32miUWfpxQm7i67iDej504m\nKbcmbZrWDzqKZzunhesftsSjzzSx6POUwqjno4xYsHUyh7nmlEuOz4/csKAjiIiIT+LzJ5EUy1uf\nfMf2jEVcftxFQUcRERHRZZey4Klp2VhuJl0vOTfoKCIiIio+El0o5Ph8SzZH2mVUSd/TIrAiIiKx\no8suCW7Yu3PZUfkn2jfWADARESkd1POR4J79IJukvP3pdGGToKOIiIgA6vlIaDm5eSzIe426yVeQ\nmpIcdJwSc8T/DswiIqKej4Q2YNIHhNJX0ikhLrloqq2ISKJQ8ZHAXvpsDOVCh3HtWY2CjiIiIvI/\nuuySoP7avJ0fy73BSZWySEpSr4GIiJQe6vlIUE+MnwppG7n33ES45CIiIolExUeCGv1NNmmuLi1P\nOjroKCIiIrvQZZcE9Pvav1lRcRJNq6nXQ0RESh8VHwnokXFvQco2erRqE3QUX+3c1VZEROKbLrsk\noHFLRlCZUzjlmH8FHcVHGjQrIpIoVHwkmHfn/ci6qu/RvtqIoKOIiIgUSMVHgun2xgtYaB/6XNs6\n6CgiIiIF0piPBLLur63Md8NokNyOzEppQccREREpkHo+EkjXV8fh0tbTu2XHoKOIiIjskYqPBJL9\n3+fZ151Ds/qHBx1FRERkj1R8JIjRs+azOfNzOtWeGHSUqNGutiIiiUHFR4J4ZOrzJHMQD2VdEHSU\n6HCaaisikig04DQBLF+1kSWpYzi9UnvSUlVPiohI6Vbs4sPM9rhmt5k9VbI44sVdI0ZCcg59r7wx\n6CgiIiKF8tLz8byZnb97o5n1B64ueSQpjlDI8fYfz3PgplY0OOKAoOOIiIgUykvxcRWQbWan7mww\ns0FAa+AMv4JJ0QycNJucjMV0PvXmoKOIiIgUSbEHCDjnJpvZLcAkMzsbuAG4CDjDOfeD3wFl7x6d\n/SRp9h86t1LdJyIi8cHT6ETn3BgzywTmAKuBJs65//qaTAo1auY81mZOp9P+2SQlaTaIiIjEhyIV\nH2b29B6eWg3MB24xC//wc87d5U80Kcx9k58gxQ6n7/2XBx0lJpzTOh8iIomgqD0f9ffQ/l+gSr7n\n9dMhRt75fDG/VZlA26ovkpqSHHScqDPUsyMikiiKVHw45zSgoJS58/XeJNkBDLr3mqCjiIiIFEux\nZruYWYqZ7TCzY6MVSAr38cJfWJo+ipbV7qFKevmg44iIiBRLsYoP51wusBxI/H7+UuyWUX2x7ZkM\nual90FFERESKzcs6H48Bj5vZPn6HkcIt+Gkl36a8zJkV76RG1fSg44iIiBSbl6m2twGHA7+b2TJg\nc/4nnXMN/AgmBWs/tD+EUnnplluDjiIiIuKJl+LjTd9TSJEsXr6auW4wJyXfxiH7Vw06Tsw5TaYS\nEUkIXlY47RmNIFK4rOcfBUtiaPuyuJSKptqKiCQKz/uvm1lDoE7k4XfOua/8iSQFmfX1Ur5JeZ5z\nUntSp1b1oOOIiIh4VuwBp2ZWw8xmAV8AAyO3eWY208yi/lPRzLqZ2Vwz+8vM/jSziWZ2ZLTPG7Tr\nRnYnaWsNRne6I+goIiIiJeJltssgoDJwjHNuH+fcPsCxhFc6HehnuD04LZLhROAsIAWYYWYVYnDu\nQAybMZcVGa9xba1HqJZRMeg4IiIiJeLlsst5wFnOucU7G5xzi8zsVmCGb8n2wDnXPP9jM7sOWAU0\nBD6O9vljLRRydJ7ShfJJx/LC/W2DjiMiIlJiXoqPJCC3gPZcvPWklFQm4T1l1gVw7qjrlT2FjVVn\n8/ARk8vEHi4iIpL4vBQLs4ABZnbAzgYzOxDoD8z0K1hRWHgr3WeAj51zi2J57ljYlrODJ+d1IXP9\nGTzQ5vyg4wROu9qKiCQGr4uMTQJ+MbMVkbaDgYXA1X4FK6LBwNHAKYUd2LlzZzIyMnZpy8rKIisr\nK0rRSq7jCyPYnrGIF08ZQVJSWZ9qWtb//iIi0ZWdnU12dvYubRs3bozKuczLb5ORHoezgH9HmhY7\n597zM1gRMjwLXAic5pxbvpfjGgDz5s2bR4MG8bP46s9/rOewZ/5NrR3N+KXfmKDjBC717qOol3YR\ncx/rE3QUEZEyY/78+TRs2BCgoXNuvl/vW+yeDzNLc85tA96N3GIuUnhcBDTZW+ERz1r074Yrt40J\nN/YNOoqIiIivvFx22WBmc4HZwPvAp865rf7G2jMzGwxkAS2BzWa2X+SpjZGiKO69OPVTFqcP4bKK\ng2hwxAGFv0BERCSOeBlwehYwjfA6G5OA9Wb2sZk9ZmZn+5quYB0JrynyAfB7vlvrGJw76rZsy+WO\ndztQccPxjL7z5qDjiIiI+M7L3i4fE15P43EzKwc0AjoAXYD7gKjOB3XOBTGdN2Za9x/AtsrfMarJ\nF5paKyIiCcnT3i6R5cyb5ruVB94h3BshHs35bhmTNz1EvVAnrjozfgbHxop2tRURSQxeBpz+BlQg\nXGh8APQGFjgtwlBil718O0nlqjL5vkeCjlLqmKbaiogkDC+XMFYDFYGakdt+hIsRKYEuwyawMnMS\nd9UZwAH7Vg46joiISNQUu/hwztUjXHQ8Sfhyy+PAGjP7xMwe8zlfmbDgp5X0XdKB/Te0ovd1lwQd\nR0REJKo8jflwzm0AJpnZHOATwmtuZBGeAdPdv3iJLxRynD3oRkhNYmbnF7WSqYiIJDwvYz4u4f8H\nmh5NeEO3j4G7Ca/9IcXQduBLrMqczAOHvU2dWtWDjiMiIhJ1Xno+XgA+BF4EZjvnvvU3Utnx7rwf\nGb26M//OvYleV18QdBwREZGY8LLOR41oBClrtuXs4OJR11AuaX/e794v6DhxQVNtRUQSg5fLLnnA\n/s65Vbu17wuscs5pZawiaPHkE2yu8gVDTvqYmvtUCjpO6ec0FkZEJFF4mWq7p58C5YGcEmQpM56f\nPIdZeT05lfu56fzGQccRERGJqSL3fJjZ7ZG7DrjRzDblezoZOB343sdsCWnhz39y2+zWVMlpzPQn\nHww6joiISMwV57JL58ifRnhzt7x8z+UAv0TaZQ+25ezgtAFtoHweszuNo2JaStCRREREYq7IxYdz\n7hAAM3sfuMQ5tz5qqRJUk17d2ZDxEQMbzqLeYfsHHUdERCQQXlY4PcM5t97MUs3sqMjOtlKIrsMn\nMjelDy3SnqRTy9ODjiMiIhKYYhcfZlbBzF4BtgDfAbUi7YPM7D6f8yWEd+f9SJ8fruOAjZcwqevd\nQccREREJlJfZLk8CdQmvcLotX/t7wBU+ZEooy1dt5MLRF5OyvSafdxum5dNLQhsni4gkBC+XTFoB\nVzjnPjOz/D8NvgMO8ydWYtiyLZcGj7dme4VfebvVpxxUvUrQkeKYijYRkUThpfioDqwqoD0dtATl\nTqGQ4/gHO7G2yiz61pvOBSfWCTqSiIhIqeDlssuXQIt8j3cWHDcCn5Y4UYJo2bsfi9OHcF21Idx9\nyZlBxxERESk1vPR83A9MNbOjI6+/I3L/ZKCJn+HiVZdhE5i8vQsnh+5n2O3XBx1HRESkVPEy1fZj\noB7hwuNb4BzCl2EaO+fm+Rsv/gybMZenll7NwX+1ZvZDjwQdR0REpNTxtEaHc24p0N7nLHFv4pyF\n3DDrfCpvq8/XvYZRLtnLVS0REZHE5qn4MLMk4HCgBrv1njjnPvQhV9yZ+dV/ueytsymfezDfdpvM\nPlUqBB0p4TiNZxYRSQjFLj7M7CRgDPAv/jn/0RHeZK5M+WLJr5w7+iySXRW+vGMG/9ovM+hICcc0\n1VZEJGF46fl4gf+f8fIHZXx67eLlqzl1yNmQ7Pjohvc4pnaNoCOJiIiUal6KjyOAy5xz//U7TLz5\n+Y/1NOx/LjtS1zOt9UecWOfgoCOJiIiUel5GRH5OeLxHmbZkxRqO6dOMbeWX8doFMzi74RFBRxIR\nEYkLXno+BgH9zKwm4am2ufmfdM4t8CNYafbdL6s4fuBZbE9dybgW73PZaccFHUlERCRueCk+3oj8\nOTRfmyM8+DThB5x+vfQPTnyuGTvKrefNiz+g5UlHBx1JREQkrngpPg7xPUWc+GLJr5zy4pnkJW9h\nSuvZnHv8kUFHKlOcdrUVEUkIxS4+nHPLohGktJv+5Q9cMPZcnDlmXv0hTeseGnSkMkZTbUVEEoWW\n4CyCV6Z/zvnjTyEpVIGPrlfhISIiUhIqPgrx8OjJ3PjhmVTadiSL7vmYxkfXCjqSiIhIXPO0vHpZ\n0W7gUIavvYn9t1zIwl5jtGS6iIiID1R8FCAUcjR7pBcf8DBHb+vIV32eJTUloSfxiIiIxIzn4sPM\nUil4Y7nlJQ0VpDUbt9CgZztWZIzj7KTHmPZEN5KSNNhRRETEL142ljuC8BofJ+/+FHG+zseXP/xG\nk+cvYkvFxdxb6w36tLsk6EiSj3a1FRFJDF56PoYDO4ALSKCN5YbNmMuN77XCypVj7DlzuKJJvaAj\nyS7U+yQikii8FB/1gIbOue/9DhOU9s+N4OU/OlIptx5zOk3kuENrBh1JREQkYXkpPhYB1fwOEoQN\nm7bRuNcdfJ/+Ikdsb8fcRwaTWSkt6FgiIiIJzcs6H12BPmbW1Mz2NbMq+W9+B9ydmZ1mZpPM7Dcz\nC5lZSy/v8+GCnznggVP4vvwIrq36Mj88NVSFh4iISAx46fl4L/LnzN3aYzXgNB34GngFmODlDXqO\nmULPb68mOTmTUWd8wlVnNvA1oIiIiOyZl+LjDN9TFINzbhowDcDMijUK8a/N2znzsfuZV/5pamy7\ngM+6jOSQ/atGJaeIiIgUzMvGcrOjESTa5nz3C6e+eiNbKy+kZfl+vNH3Tsola3X5eKKptiIiicHT\nT9/IuItRZvaJmR0YabvGzE71N55/bv/gKvKStjDmzM956767VHjEGdNUWxGRmAuFovNLX7F/ApvZ\npcB0YCvQACgfeSoDuN+/aP6qlXsuK3rMI6tp/aCjiIiIlHqLl6+mef+uUXlvL2M+egAdnXMjzaxN\nvvY5kedKpdrL1nPjtVm7tGVlZZGVlbWHV4iIiJQd2dnZZGdnA7BoxR8s/XsBbA9F5Vxeio+jgA8L\naN8IZJYsTvT079+fBg00q0VERKQgWVlZnNT0PM7s24lfGr7N/hta8cyZt3BFy3N8P5eXgQ8rgcML\naD8V+KlkcQpnZulmVtfMdq5/fmjk8cHRPreIiEiieuDVtzm8/7H8kvoOHaqP5Nd+Ezj8wH2jci4v\nPR8vAQPM7HrC63ocYGaNgb7AI36G24Pjgfcj53ZAv0j7COD6GJxfREQkYSz4aSUtnr2dXzNep3ru\n+Uy+4UUaHXVQVM/ppfh4knCPyUygIuFLMNuBvs65QT5mK1Bkqq+mqoiIiJRAKORoN+gVXl15L6Sk\ncFvNMQx4oA1JSdGfXehlnQ8HPGZmTxG+/FIJWOSc2+R3OJH8wv/riYhISU36bBFtx97CxqqzOTz3\nOqbc2ZcjDorOJZaCeOn5AMA5l2NmfwN/q/CQqHNa50NEpKRWrtvERU8/wtykp0lJrk3vY96ly2Vn\nxTxHsYsPMysHPATcTrjXAzPbBAwCejrncn1NKCIiIiUSCjnuGTaeAd93JlR+Lc3KPcT4nvcEtqGq\nl56PQcAlQBfg00hbY+BhYF/gZl+SiYiISIlNnLOQG1/vzLqq71Fzx0W8fu0znHps7UAzeSk+rgTa\nOOem5mtbYGYrgGxUfIiIiARuyYo1XPrsg3yXNoSU5MN46PB3ePiqFkHHArwVH9uBXwpo/xnIKVEa\nERERKZFNW3O4ZuBg3trQE5fkuKhCX0bddyuVKqQGHe1/vBQfzwIPmFk759x2ADMrD3SPPCciIiIx\nFgo57h76Os8uvp8dlX7mmFAHXr+1J3VqVQ862j94KT7qA82AX83sm0hbXSAVmGlmE3Ye6Jy7pOQR\nRcIcmmorIlKQAW/NpvsHXdicOZcaoRa8cM6bXHzKsUHH2iMvxccG4I3d2lb4kEVkLzTVVkRkd+M/\nWsCtE7qzKvMd0mlE/7rvc2erpkHHKpSXRcbaRSOIiIiIFM30L3/gpuyHWF75NcolHcqdB75Gvwcu\nj8nqpH7wvMiYmVUnvMMtwBLn3Gp/IomIiEhBPl20nOuG9eKHisNJTtqfazKH8MJ911ExLSXoaMXi\nZZGxdMJrfbTl//dYyTOzkUAn59wWH/OJiIiUeZ8uWs6NI55gUflXsKRMWlXsy7B7Owa2SFhJeen5\neBpoAlwIzIm0nQoMJLzDrNb5EBER8cGc75bRfsQTLE4bilkG55d/lKF33kLNfSoFHa1EvBQflwKX\nOec+yNc2xcy2AuNQ8SEiIlIiM7/6L7dl9+H7tOFYUgbN0x7jlTtvjvuiYycvxUdF4M8C2ldFnhOJ\nEk21FZHENnHOQjpPeIJllcaSZNUTrujYyUvx8SnQ08zaOue2AZhZBcKbzX2611eKeGSaaisiCezl\naZ/xwPQnWZn5FslJtbiiyiAG392OfapUCDpaVHgpPu4ApvPPRca2Aef6FUxERCSR7cgL0St7CgO+\n7MNfVT8iNekobth3GAO7XhV3s1eKy8s6HwvN7AjgKuDfkeZsYLRzbquf4URERBLNX5u3c9ewbEb9\n9BTbMxZRyU7ivtoTeeTqlpRLTir8DRKAp3U+ItNpX/I5i4iISMJasmINtw0fwqy/nyWUvpIa7gIe\najiEjs1PiZvFwfzieZExERERKdyUud9z7xsDWJQyAnDU4Vr6nn8nzU/4d6GvTVQqPkRERHy2Iy/E\n4+OmM3DuANZmTifJ1aRZ6v08164jRx1cLeh4gVPxIXHDOU21FZHS7fe1f3PH0BG8tXIQuVV+oCIN\n6VB9JH3uaU2V9PJBxys1VHxInChb10NFJL5M+mwR3Sc+x8JyI6HcVg52l9K90TDan9e4zI3nKIqS\nbCzXEKgTebjIOTffn0giIiKl35ZtuTw0ZhIvL3iODVXfJ4n9OK1cZ/pdeRONjjoo6HilmpeN5WoA\nY4GmwIZIc6aZvQ+00e62IiKSyD5fvIJ7x77EnK0vE0r/g8qcwm01x/BE20upVCE16HhxwUvPxyCg\nMnCMc24xgJkdDYwgvLlcln/xREREgpeTm8eT42fw/NwXWFnlHdhRkWPsGh48owOtT68bdLy446X4\nOA84a2fhAeCcW2RmtwIzfEsmIiISsC+W/ErXsUP5cNMr5FVaThp1uTJzME+1vZID9q0cdLy45aX4\nSAJyC2jPjTwnIiISt7bl7OCx16by0vyX+LPKZNhRgaPIomvj9lx7ViMNIPWBl+JjFjDAzLKcc78D\nmNmBQH9gpp/hRPJz2tVWRKJo5lf/5YEJQ/l8+3BC6X9QgfpcmTmY3ldncVD1KkHHSyheio/bgEnA\nL2a2ItJ2MLAQuNqvYCK70m8aIuK/NRu30GP0G7z2w1A2VP0AdmTwn6Sr6dbkBrKa1g86XsLysrHc\nCjNrAJzF/28st9g5956vyURERKIgFHK8PP0z+r8/jO/LjYXyf5NJU27ebxSPXnlJwm5jX5p4mWrb\nFnjNOfcu8G6+9lTCU21H+phPRETEF1/+8BsPvD6KWWuHk5PxPcmuFqendqZnq2tpWvfQoOOVKV4u\nuwwDpgGrdmuvHHlOxYeIiJQK6/7aysNj32LM4uGsrfIu5KVS2y6lw9EDueeSZmVmC/vSxkvxYVDg\nyL+DgI0liyMiIlIyoZDjhSlzGDh7BEtSxkH5v6jsTqHtPkN4JOtyatXICDpimVfk4sPMviJcdDhg\nppntyPd0MnAI4R4RERGRmJv19VJ6vfkqcza9yo7KP5Hs/sWpKbfTo8U1nHv8kUHHk3yK0/PxZuTP\nesB0YFO+53KAX4A3/IklIiJSuKW/r6NH9mu8s+JVNlX9FHIqc4Rdzq31hnLrBafpskopVeTiwznX\nE8DMfiE84HRbtEKJFETrfIgIwIZN23j89SmMWvAqf1SeDBaiup3DdTXH8FCbi6iWUTHoiFIIL1Nt\nR0QjiMjemNb5ECnTduSFeH7yxwz+eBRLkl/HpW2gIg1pld6HR1pncewh+wUdUYrBy4BTERGRmHjj\n4295atpovtyeTV6l5ZQL1eaUCrfR7fyraH7Cvwt/AymVVHyIiEip8umi5fSamM0Ha0ezLeNbbMe+\nHJ3UmluPv4oO55+svVUSgIoPEREJ3JIVa+j5+ngmLx/DX1U/gtwK1LKLuPawx+ly6TlUqpAadETx\nkefiI7Ki6SHAUufcjsKO95OZ3QrcA9QEvgE6Oee+iGUGEREpmZXrNvHIuEmMX5LNqsrTwBzV7Gw6\nVB/Jg1e00pb1CczL8uoVgUHAtZGmI4GfzGwQ8Jtz7kkf8xV0/iuAfsBNwFygMzDdzI50zq2J5rlF\nRKRkNm3Noc8bMxj51RiWpb0FqVuoxElcVrk/D112uQaOlhFeej6eAOoCTdl1UbH3gIeBqBYfhIuN\nITv3kDGzjkAL4HqgT5TPLQFyTlNtReLRjrwQz779IUM+yWZJ0nhchXWUd0dzdlp37m/ZRvuqlEFe\nio9WwBXOuc/MLP9Pg++Aw/yJVTAzSwEaAo/vbHPOOTN7D2gczXNL0DTATCSehEKOMe/P55mZY/gq\n9zVClX4jOfQvTix/E3efmcUlp/xHA0fLMC/FR3X+uakcQDoF7/nip2qEl3L/c7f2P4GjonxuEREp\nxNQvltB7cjaf/D2G3Co/YnnVOTb5Cm5tlEX78xqr4BDAW/HxJeHLHIMij3cWHDcCn/oRSkRE4sfn\ni1fwyJuvMWtVNlsz50NOFQ63S7i+zrN0bnUmaamaWCm78vJ/xP3AVDM7OvL6OyL3Twaa+BmuAGuA\nPGD3EUn7ASv39sLOnTuTkbHrToZZWVlkZWX5GlBEpCzYOTV2yvJsNlb9EHaU50C7gCtrdef+y5uT\nWSkt6IhSTNnZ2WRnZ+/StnFjdDar97K8+sdmVg+4D/gWOAeYDzR2zn3rc77dz51rZvOAZsAkADOz\nyOOBe3tt//79adCgQTTjiYgktJXrNtFr3FuMXzKG1ZVngDn2tbO4sdpwHri8lbaqj3MF/UI+f/58\nGjZs6Pu5PPWFOeeWAu19zlJUTwPDI0XIzqm2FYHhAeUREUlYm7bm8MT4abz61RhWVJwEKVup7E7m\n8irP8NBll3NM7RpBR5Q45GWdj/eAUcAE59xf/kfaO+fcODOrBvQifLnla+Bc59zqWGeR2NKutiKx\n8f9TY8ewJHk8Lm09afyHc9MepEerNpx6bO2gI0qc89Lz8R3htT4Gm9lkwoXIFOdcrq/J9sI5NxgY\nHKvzSfC0q61IdIVCjtc+/Jqn3x3D/JxsQpV+o1yoNo3Tbuaes7K4+JRjg44oCcTLmI87zKwzcBZw\nJTASyDOz8cBo59xsnzOKiEiUfPDNTzz61hg+2jCGnIzFWF41jklqzW2NrtLUWIkar2M+QsAMYEZk\nhdELge7ADYTX4RARkVJqyYo1PDRuHFNWjObvqp9ATjqH2MVcd2Q/7rn4LCqmpQQdURJciSZfm1lN\noA1wNXAc4QGgIiJSyqz7ays9X5vE2EWj/reJW3U7l2trjuGB1i2pUTU96IhShngZcFoFuJTwJZem\nwE/AaMJLri/1NZ2IiHiWk5vHoLdnM+TTUfyYMh7K/00lTuLyKs/Qs3Vr6tSqHnREKaO89Hz8CawH\nXgO6Oee+9DeSiIiUxFuffMcTU17ly+2jyav0K+VCh9G0/N30aHkVzeofHnQ8EU/FR0tgZmTch0gM\naaqtyJ5898sqHhiXzbQ/RrI1cz62Yx+OTrqC20+4hhvPPUkDR6VU8TLb5d1oBBHZO/3DKbK7vzZv\np+fYtxm1cASrqkwFl8T+XMAdtR+ge+vmVKqQGnREkQIVqfgws/lAM+fcejP7ir38Cuqc0xrmIiJR\nEgo5Rrz3Bf1mDmdR0lhc2nrSOYHLKw/gsTZtOOKgfYOOKFKoovZ8vAVsz3df/d8iIjG04KeVdBv7\nKjPXDGd7xiKSQgdyYmoHup/XlgtOrBN0PJFiKVLx4Zzrme/+w1FLIyIi/7Npaw69xr7DiAVDWVVl\nGoTKUcsupmOd/tx9cTNSU7SsksQnL1NtfwIaOefW7taeCcx3zh3qVzgRkbLorU++45F3hjI/71Vc\nxdWkcwK+bhFlAAAgAElEQVRZGc/yxFVt+Nd+mUHHEykxL7NdalPwKqblgYNKlEZEpIxauW4TXV8d\nyxu/vMzmzM+xUDXqJ19Dj7Ou174qknCKXHyYWct8D881s435HicDzYCf/QomIpLodg4e7f3uSyxJ\nHQspm6nOedx88HgeyrpQs1UkYRWn5+PNyJ8OGLHbc7nAL8DdPmQSKZBzGucsiWH5qo10eXU0b/06\nhG2ZC0imFk1T7+HJ1tdzYp2Dg44nEnVFLj6cc0kAZvYz4TEfa6KWSuQftM6HxLdQyPHqzC95fMYQ\nfkjNhnLb2Z+W3HJUb7pcerYGj0qZ4mWRsUOiEUREJBGt2biFe0Zk8/rPz7Mlcx7JHEyz8t3ok3U9\nDY44IOh4IoHwtKutmaUDTYBawC4XJZ1zA33IJSIS16Z+sYTuE17ga4bjym+kBs25+7C36XHF+erl\nkDLPy1Tb+sAUoCKQDqwDqgFbgFWAig8RKZNycvN4ZOwUBs8bxLqq72KuGieU60CfSzpw+nHqNBbZ\nyUvPR3/gbaAjsBE4ifCA01HAAP+iiYjEh5//WM/tw4cybe1z7Kj8M+nWiA7VR/Jk28vJrJQWdDyR\nUsdL8VEP6OCcC5lZHlDeOfeTmXUhPAtmgq8JRURKqalfLOHe8QP4rtwISM7lEHcFDzQeS7tzTgg6\nmkip5qX4yAVCkfurCI/7WEy4F0RzxCRqnLYUklIgFHL0eeM9nv7kGVZnTiGJ/TgjtSsD23bg2EP2\nCzqeSFzwUnx8BTQCfgRmA73MrBpwDbDQx2wi/2OaaisB27Q1h85Ds3l1aV+2ZyykAvW4sdpw+l3X\nhirp5YOOJxJXvBQf9wOVI/e7AyOB5wkXI9f7lEtEpFRY9ucGOr48hBkbBxJK/50a7gK61R3E7S2b\nkJSkoljECy/rfHyZ7/4q4DxfE4mIlAJf/vAbHYc/zTx7EZJyOMq1pe95d2n7ehEfeFrnQ0QkUU39\nYgl3jnuKHyqMBFeRk8vdzvPXdeK4Q2sGHU0kYRSp+DCzr6Boo/2ccw1KlEhEJACvzf6au996jN+q\nvEGS1aRFhcd5odNNHFS9StDRRBJOUXs+3iz8EBGR+DNsxly6TX2UPzPfplzSoVyTOYRn722rQaQi\nUVSk4sM51zPaQUQKo6m24qchUz6h+4yerK06g9Sko+hQfSTPdMsiLVVXo0WizeveLpnAZcBhwFPO\nuXVm1gD40zn3m58BRQBwmlUg/hg2Yy73Tn2QtZnTKZ90LLcfMJan7r9M+62IxJCXvV2OA94jvKhY\nbeAlwvu7XEJ4wbG2PuYTEfHF6Fnzuevth1iV+Q6pVofOB42jT49LKZecFHQ0kTLHS8/H08Bw51wX\nM/s7X/sUYIw/sURE/DH1iyW0z+7BbxnjSUk6klv2G03/+69QT4dIgLwUH42ADgW0/wZoLpqIlApf\nLPmVa17uyZKKw0hOPpAb9h3Gs92u1pgOkVLAy7dwO1DQ3LMjgdUliyMiUjI//7GeK559gi9sEJZc\niYvT+zG8S0fNXhEpRbwUH5OAB82sdeSxM7NaQG/gDd+SiYgUw6atOVw76AUmruuJS95Ok3JdGXXH\nXVqnQ6QU8lJ83A2MJ7yjbQXCm8vVBD4lvNeLSJRoqq38Uyjk6DbyTfp/25XcSks5Ku96xrbvRb3D\n9g86mojsgZe9XTYCZ5vZKUBdoBIw3zn3nt/hRHbSrrZSkPEfLeDGN25nY9XZ7Jt3DkPOfoNLT/1P\n0LFEpBDFKj7MLAWYBnR0zs0B5kQllYjIXiz9fR0XD3yQb8s/T2ryEfQ8cgoPZp0fdCwRKaJiFR/O\nudzIOh8iIjGXk5tHu2dfJvvP7rikHC5I60P2fZ2oVCE16GgiUgxexnyMAm4A7vM5i4jIHo2aOY8O\n73RgS+Y8Dt9xHW/c8oR2mhWJU16Kj3LA9WZ2FjAP2Jz/SefcXX4EExEB+HX1X1zw9AN8k/osaXYM\nLzSaQ4fmJwcdS0RKwEvxcSwwP3L/yN2e03QEEfFFKOToOnwCT39/O6GUDbRI6824bndQMS0l6Ggi\nUkJeZrucEY0gRWFm9wMtgHrAdufcPkFlkdhzTrVtWbHgp5WcN+gW/sicSM3clky4bhCNj64VdCwR\n8Um87aiUAowDng86iMSaptqWBaGQ46bnRlLvpaNZmTqHzgeN44/+b6nwEEkwcbXJgXOuJ4CZXRt0\nFhHx1xdLfqXFCx1YnTmF2jlXMe3OZzjq4GpBxxKRKIir4kNEEk8o5LjtxTE8v+xWklLS6XHoJB65\n5sKgY4lIFKn4EJHALP19HWf0vZkVGeOonXMVs+4ZxCH7Vw06lohEWeDFh5k9AXTdyyEOqOOc+6Ek\n5+ncuTMZGRm7tGVlZZGVlVWStxURj/qMf49un1+HK7+Z2w8Yy4CHrgg6kkiZlp2dTXZ29i5tGzdu\njMq5PBUfZnYN0BE4BGjsnFtmZncCPzvn3irm2/UFhhVyzE8eYu6if//+NGjQoKRvIyIltGlrDk0f\nvZ95qf2omtuM6TcOp9FRBwUdS6TMK+gX8vnz59OwYUPfz1Xs4sPMbgZ6Ac8Q3sU2OfLUBuBOoFjF\nh3NuLbC2uDlEJP58uOBnmg9rw+bKX3FR+acZ3+8OyiXH26Q7ESkpL9/6TkB759xjQF6+9i+BqG4n\naWYHm1ld4F9AspnVjdzSo3leKR2c1rCLa12HT6RJdn22J69i+GlzePO+zio8RMooL5ddDgG+KqB9\nOxDtIqAX0Dbf450rrZ4BfBjlc0uATOt8xK1NW3M4tde9fJM2kAO2X8InXV/hX/tlBh1LRALk5deO\nnwmvMLq784DFJYuzd865ds655AJuKjxESqGvl/7BAd3O4JtyL3BZxUGs6DtehYeIeOr5eBp4zszS\nCC87eYKZZQHdgBv9DCci8evFqZ9y8/uXQqrx0imzufG8k4KOJCKlhJe9XV42s63Ao0BFYAzwO3CH\nc26sz/lEJA5d/cyLjF53G5VzTuDj28dz3KE1g44kIqWIp6m2zrnRwGgzqwhUcs6t8jeWiMSjLdty\nOf7BTixOH8Ix22/msyeeoVKF1KBjiUgpU+wxH2Y2y8wyAZxzW3YWHmZWxcxm+R1QROLD8lUbOfi+\n5ixOG0rbzJdY2HuwCg8RKZCXno+mQEH/oqQBp5UojcheaKpt6TXnu2U0e6UFORV+o1/96dx18RlB\nRxKRUqzIxYeZHZfv4dFmlv8ibjLh2S6/+RVMZFeaaltajXzvS9rNuJCkpApMavUJF5xYJ+hIIlLK\nFafn42vC+6w4oKDLK1sJL0AmImVE91cn8fiSNqTnHsfnd0zimNo1go4kInGgOMXHIYR//fwJOAFY\nne+5HGCVcy6voBeKSOK56bmRvLTqeg7YchHf9hzFPlUqBB1JROJEkYsP59yyyF2thyxSxl3aZyAT\ntt7BkVtu4NveQ0hNSS78RSIiEZ6m2gKY2dFALXYbfOqcm1TSUCJSOoVCjmaP9OIDHub4nHv4vHcf\nkpI0HkdEisfLrraHAhMJbyLn+P+RgDunIuhXIJEEtCMvRKMed/F12gDOSX6cqY/cp8JDRDzxcgll\nAOH9XWoAW4BjgNMJ72rb1LdkIv+gqbZB2ZEX4rhut/B1+YG0qTSY6T26qfAQEc+8XHZpDJzpnFtj\nZiEg5Jz72My6AQOB+r4mFEG72gYpFHLUv78Tiyu+SLt9X2Fop3ZBRxKROOel5yMZ+Dtyfw1wQOT+\nMuAoP0KJSOkQCjkadL+ThRUH07bqiyo8RMQXXno+FgJ1CV96+RzoYmY5wE2Ep+GKSAIIhRyNetzN\nN2kDubLKC4y4Q5tWi4g/vBQfjwLpkfsPAu8AHwFrgSt8yiUiAQqFHI0fvI/55ftzefqzjO7cIehI\nIpJAil18OOem57v/X+DfZrYPsN45pxGBIgnggieeYm5KH1ql9WfcPbcGHUdEEkyxxnyYWYqZ7TCz\nY/O3O+fWqfAQSQw3PTeSqTu6cmqoBxO73hl0HBFJQMXq+XDO5ZrZcrSWhwRA9W30PTp2Gi/9eQNH\nbr2B2b17BR1HRBKUl9kujwGPRy61iMSIptpG24h3v+CBby9jv7/P55vHXtA6HiISNV4GnN4GHA78\nbmbLgM35n3TONfAjmIjEzrvzfqTdzOZU2nYcC3uOJS3V884LIiKF8vIvzJu+pxCRwPz8x3pajGlO\niqvG/HveplpGxaAjiUiC8zLbpWc0gohI7G3L2cHxva9gR9o6ZmV9wREH7Rt0JBEpA9S3KlKGnfJw\nF9ZlzKJv3Rk0rXto0HFEpIxQ8SFSRrV/bgTzy/fn0ooDufuSM4OOIyJliJfZLiKBcNrV1jcvT/uM\nl1fexJGbbmDc3bcFHUdEyhj1fEhc0K62/vnyh9/oMOtiKucczxePPacptSISc8Xu+TCzM6IRRESi\nb1vODpoNzgKXzCd3TKBKevmgI4lIGeTlsss0M1tqZj3M7GDfE4lI1Jz3+KP8lTGHAadnc+wh+wUd\nR0TKKC/Fx4HAs8BlwE9mNt3MWptZqr/RRMRPz7z5AbPdI5yZ/DC3XXha0HFEpAwrdvHhnFvjnOvv\nnKsHnAj8AAwmvOLpQDOr63dIESmZxctXc/cnV5G5sQlTu90fdBwRKeNKNNvFOTcfeIJwT0gl4Hpg\nnpl9ZGbH+JBPREooFHI06X8dLimHmbeOIjVF+0KKSLA8FR9mlmJml5nZFGAZcC7hPV/2I7zvyzLg\ndd9SiohnF/fpz+rMKTx03AgaHHFA0HFERIo/1dbMBgFZhLcZfRXo4pxbmO+QzWZ2D/C7PxFFwrTO\nR/GN+/AbJm25j4Z5d/PQlc2DjiMiAnhb5+NooBMwwTm3fQ/HrAE0JVd8pLUoimvLtlyum9iO8nYU\nHzz2eNBxRET+x8vGcs2KcMwOYLanRCLii4v79mVrlW8YcdrnVKqgyWgiUnp4XuHUzI4GagG7/Kvm\nnJtU0lAiUjKTPlvEjO0Pc6K7l7ZnHR90HBGRXXgZ83EoMBH4D+D4//7wnRfkNZReJEA5uXlc+dr1\npCYdwpSHHgo6jojIP3iZ7TIA+BmoAWwBjgFOB74EmvqWTEQ8ubzfADZnzGVgs6HsU6VC0HFERP7B\ny2WXxsCZzrk1ZhYCQs65j82sGzAQqO9rQhEpsnfn/cikTd2pn3cHHZqfHHQcEZECeen5SAb+jtxf\nA+xcOGAZcJQfoQpiZv8ys5fN7Ccz22JmP5rZw2aWEq1zSmmjqbZ7Ewo5Wr96E+W2HcC0Lo8GHUdE\nZI+89HwsBOoSvvTyOdDFzHKAm4CffMy2u38THl/SHlgKHAu8DFQEukTxvFIKmKbaFqrzK+PYUPUD\nHq8znRpV04OOIyKyR16Kj0eBnf+yPQi8A3wErAWu8CnXPzjnpgPT8zX9YmZ9gY6o+JAybs3GLTz3\nw73U3NGSbq3PCTqOiMheeVnnY3q++/8F/m1m+wDrnXOx7hfPBNbF+Jwipc6l/XuTl/Ynoy99Ougo\nIiKFKvaYDzPLiBQb/+OcWwdUNbMqviUrPMfhhPeTeSFW5xQpjeZ8t4wPd/ThJO7izHqHBR1HRKRQ\nXgacjgVaF9DeOvJcsZjZE2YW2sstz8yO3O01BwJTgdecc0M9/B1EEsYVr9xDUk5V3rjz/qCjiIgU\niZcxHycCdxbQ/gHwmIf36wsMK+SY/w1kNbMDgFnAx865DkU9SefOncnIyNilLSsri6ysrGJEFSld\nnp74Pr9ljKdD9ZEcsG/loOOISBzLzs4mOzt7l7aNGzdG5VxW3GEaZrYZOMk59+1u7f8BPnfOVfQx\n3+7nPpBw4fEFcE1RxpiYWQNg3rx582jQoEG0okmU7XNnMyon12BZv+zCDy4jtuXsIPO++qS4Sqzv\nO4dyyV46MkVE9mz+/Pk0bNgQoKFzbr5f7+vlX6u5hKfV7q4jMK9kcfYs0uPxAeH1RLoANcxsPzPb\nL1rnlNJEU213127QS2zPWMizzQeq8BCRuOLlsksP4D0zqwvMjLQ1AxoB0ZzjdzZwaOS2ItJmhFee\n0n4yUqas+2sr4/58hNp5V3Ht2Y2CjiMiUizF/nXJOTeH8BLrvxIeZHoh8F/gOOfcR/7G2+W8I5xz\nybvdkpxzKjykzLl+8AuEKqzi5asfDjqKiEixeen5wDn3NXClz1lEpAhWrd/MpPVPckToWprVPzzo\nOCIixebpQrGZHWZmj5rZGDOrEWk738yO8TeeiOzuusHP4cqv4+VregQdRUTEEy+LjDUBviU85fZS\noFLkqbpAT/+iicjufl/7N9P+7kOdbTdw+nGHBB1HRMQTLz0fTwI9nHNnAzn52mcBJ/mSSqQAsV+9\nv/Rp+9wAXMrfDGvXPegoIiKeeSk+/gNMLKB9FVCtZHFECqZdbWHZnxuYtbUfx+V24MQ6BwcdR0TE\nMy/FxwZg/wLa6wO/lSyOiOxJ2+f745K3MaJ9t6CjiIiUiNe9XXqbWU3Ca2wkmdkphJdJH+lnOBEJ\nW/r7Oj7M6U9Ddyv1Diuo9hcRiR9eio/7ge8JL/RVCVgEfAh8AjzqXzQR2emWoUMgOYcR7bsEHUVE\npMSKvc6Hcy4HaG9mjwDHEi5AvnLO/eh3OBGBTVtzeO+vQfybthxTu0bQcURESszTImMAzrnlwHIf\ns4hIAe4Z9hqh9D/ofW5Bm0mLiMSfIhUfZvZ0Ud/QOXeX9zgie1P2ptqGQo6RP/anGufR8qSjg44j\nIuKLovZ81N/tcYPIa5dEHh8J5BHFXW2lrCubU20HTprN1syveKDO9KCjiIj4pkjFh3PujJ33zewu\n4G/gWufc+khbVWAYELWN5UTKot6z+1PejqHrZWcHHUVExDdeZrvcDXTbWXgARO73iDwnIj54d96P\nrMx4mza1O5OUVDZ7fkQkMXkZcFoFqF5Ae3WgcsniiMhOd78+AHPVeOb6q4KOIiLiKy89HxOBYWZ2\niZkdFLldCrwCTPA3nkjZtPT3dXybPIwmFW8ls1Ja0HFERHzlpeejI+HVTMcAKZG2HYSLj3t9yiVS\npt027CVIyuPZa28OOoqIiO+8LDK2BbjFzO4FDos0L3XObfY1mUgZlZObx7sbnuUIrtKiYiKSkEqy\nyNhmYIGPWUT2ypWRdT76vPEueZV+pfvJHYOOIiISFZ6LD5FYsjK0zseQuUMpzzFc0+z4oKOIiESF\nlwGnIhIlS39fx6/pb3Fujes1vVZEEpaKD5FS5L7RY8BCPN7m6qCjiIhEjS67iJQik38fRk1aaKCp\niCQ09XyIlBLjPvyGrZnzaVe/XdBRRESiSj0fIqXEE1OHYaEa9GjdPOgoIiJRpZ4PiRuJPNV209Yc\nvnGjaZB8NRXTUgp/gYhIHFPPh8QHS+yZH73GvoOrsIbuZ+mSi4gkPhUfIqXAiAXDSKcRF59ybNBR\nRESiTpddRAK24KeVrKoylQsPUq+HiJQNKj5EAtZt7KsQKsfjV7YJOoqISEzosotIwD5YPY4D7UIO\n2b9q0FFERGJCPR8iAfp00XK2ZH7JxUddEnQUEZGYUfEhcSTxpto+9c5E2JFK10taBB1FRCRmdNlF\n4kKi7mo76/cJVLezOKh6laCjiIjEjHo+RAKy8Oc/2Zj5Ec0P0SUXESlbVHyIBKTPW5PAGV1btQw6\niohITOmyi0hApi6bQKY7nTq1qgcdRUQkptTzIRKAZX9uYE3lmTQ7UJdcRKTsUfEhEoDeEydDci73\nXtAq6CgiIjGnyy4SNxJpV9tJP0wk3TXixDoHBx1FRCTm1PMhcSGRptqu2biF3ypM5dR9dclFRMqm\nuCo+zOwtM1tmZlvN7HczG2lm+wedS2IrOzs76Agl0u/NGZC6hc7nqfjYKd4/U9mVPk8pTFwVH8As\n4HLgSOAS4DDg9UATSczF+z9s476dQPmNx3Du8UcGHaXUiPfPVHalz1MKE1djPpxzA/I9XGFmTwIT\nzSzZOZcXVC6Rotq0NYefU97m1NROQUcREQlMXBUf+ZnZPsBVwBwVHhIvBk/+EJe2gZubXBx0FBGR\nwMTbZRfM7Ekz2wSsAQ4GNFdR4sb4r2aQtLkmV5xeL+goIiKBCbznw8yeALru5RAH1HHO/RB53Ad4\nGfgX8BDwKnDBXl6fBrB48eKSh5XA5K7aiGMb8+fPZ+PGjcyfPz/oSJ4s+GUSNa0+X3/9VdBRSpV4\n/kzln/R5Jo58PzvT/Hxfcy7YtRPMbF9g30IO+8k5t6OA1x4IrAAaO+c+38P7XwmMLnFQERGRsusq\n59wYv94s8J4P59xaYK3HlydH/iy/l2OmEx4b8guwzeN5REREyqI0oDbhn6W+Cbzno6jM7ASgEfAx\nsB44HOgFVAeOdc7lBhhPREREiiieBpxuIby2x3vA98BLwNdAUxUeIiIi8SNuej5EREQkMcRTz4eI\niIgkABUfIiIiElMJUXyY2a1m9nNkw7nPzKxRIcc3NbN5ZrbNzH4ws2tjlVUKV5zP08yamFlot1ue\nmdWIZWYpmJmdZmaTzOy3yGfTsgiv0fezlCru56nvZ+lmZt3MbK6Z/WVmf5rZRDMrdNMpP76jcV98\nmNkVQD/CC47VB74BpptZtT0cXxt4B5gJ1AUGAC+b2dmxyCt7V9zPM8IBRwA1I7f9nXOrop1ViiSd\n8MDwWwh/Tnul72epV6zPM0Lfz9LrNGAQcCJwFpACzDCzCnt6gV/f0bgfcGpmnwGfO+fuiDw2wguP\nDXTO9Sng+N7A+c654/K1ZQMZzrnmMYote+Dh82xCeLfjqs65v2IaVorFzEJAK+fcpL0co+9nnCji\n56nvZxyJ/JK3CjjdOffxHo7x5Tsa1z0fZpYCNCRcgQHgwtXUe0DjPbzspMjz+U3fy/ESIx4/TwAD\nvjaz381shpmdHN2kEkX6fiYefT/jRybhnqp1eznGl+9oXBcfQDXCq5z+uVv7n4S79wpScw/HVzGz\nva2UKtHn5fP8A+gAXEp4HZgVwAdmpp3b4pO+n4lF3884Eellfgb42Dm3aC+H+vIdDXx5dZGSiGw4\n+EO+ps/M7DCgM6CBiiIB0vczrgwGjgZOicXJ4r3nYw2QB+y3W/t+wMo9vGblHo7/yzm33d94Ukxe\nPs+CzCW8/L7EH30/E5++n6WMmT0LNCe8YvgfhRzuy3c0rouPyLLq84BmO9siXUfNgE/28LJP8x8f\ncU6kXQLk8fMsSD3C3b0Sf/T9THz6fpYikcLjIuAM59zyIrzEl+9oIlx2eRoYbmbzCFfUnYGKwHAA\nM3sCOMA5t7OL7wXg1siI3aGE/yNeRrjqk+AV6/M0szuAn4HvCO++2B44A9DUzFLAzNIJ/5ZrkaZD\nzawusM45t0Lfz/hS3M9T38/SzcwGA1lAS2Czme3s0djonNsWOeZx4EC/v6NxX3w458ZFpgf1Itz1\n8zVwrnNudeSQmsDB+Y7/xcxaAP2B24FfgRucc7uP3pUAFPfzBFIJrwtyAOHNBxcAzZxzH8YutezF\n8cD7hEfQO8KfFcAI4Hr0/Yw3xfo80feztOtI+HP8YLf2dsDIyP39icJ3NO7X+RAREZH4EtdjPkRE\nRCT+qPgQERGRmFLxISIiIjGl4kNERERiSsWHiIiIxJSKDxEREYkpFR8iIiISUyo+REREJKZUfIiI\niEhMqfgQERGRmFLxISIiIjGl4kMkjpnZ+2b2dLSOL+2C+PuY2b5m9qeZ1YrlefeQZZe/v5llm9ld\nQWYSKQoVHyLFYGbDzCyU77bGzKaa2X/yHVPNzJ43s2Vmts3M/ogc0zgKkS4GHoji8fJP3YE3nXPL\ngw5SgEeB7mZWOeggInuj4kOk+KYC+xHePvxMYAfwdr7nJwB1gWuAI4ALCW9Zva/fQZxzG5xzm6N1\nfFlhZhXN7D4zu3ZvPSlmVoHw1vEv+3juFL/eyzn3HbAUuNqv9xSJBhUfIsW33Tm32jm3yjm3AHgS\nODjSHZ8BnAp0dc596Jxb4Zz70jnX2zn3zp7eMNJ9PtDM+pvZOjNbaWY3RH4oDjWzv8zsRzM7r4DX\nPR25Xy3Sy3JfvudPNrPtZnbG7sfnezzAzHqb2drI6x/K93wlMxttZpvMbIWZdSrKpQ4zO9fMPjKz\n9ZHeobfN7NCintfruS2sm5n9ZGZbzOwrM7t0b1kjXgMmOedGAP82syP3cFwLYJtz7osivOeeMr5v\nZoMin/VqYFqkvbD/ZhXNbKSZ/W1mv+3l8srbQBuv+URiQcWHSAmYWSXCPRw/OufWApsit1ZmllrM\nt2sLrAYaAQOBF4DXgTlAfWAGMNLM0gp6sXNuDeHfynuaWYNItpHAQOfc+4WcdxNwAtAFeNDMmkWe\n6w80Bi4AzgWaRrIUJh3oBzQg3DuUB0wsxnm9nvt+wr/13wQcHXmPV83stD29wMxuBNKcc4siTRWB\nQ/dw+KnAvEIyFEVbYDv8Xzv3FyJVGcZx/Pv4J/sj2ookFCLlBmFaEqlQiV1EiBBoSFIX1kWaN0pJ\nUjdJkhqYGdZFiPinlAq6sCQqiZJIpCQ2S7OE1Is0rcxt3SxT8Oniecc9jXPOzJid2fL3gcPOnDnv\n+7zzzi7nmffPchswJ52r12fLgYnESNrdRH/cUqPuHcD4CzmiInLBubsOHToaPIB1wGmgOx1ngIPA\n2Mw104CjwO/ANmAJMKZOvVuBjzPP+6T612fODUvxxleVW1FV10vAt8BGYCfQP+/66rjp3GfAUmAg\ncYOclnltEJEwrCh6PzXe39DU9lH14qbHDcXOvh/gkvT6hKp6VwMbC9r2HXB/pt87q+vIXLsJWF11\nbiWwCzgC7AEOp5/rc+rYCnzeTJ8RiclJ4N7M623AiRqf/xgicRne6r8XHTryDo18iDTvI+AmYl3H\nOGAL8L6ZDQdw903A1cQ31PeASUCHmc2sU+9XlQfufgb4hbipVc79mB5eVaeeBUA/YDrwgLufbjRu\ncjjFuC7Vc3aKwd2PA3vr1IeZtZvZa2a2z8y6gAOAA9kdInlxOc/Y7cSoxQdpaqLbzLqJkamROe28\nFRCikYwAAAOLSURBVBgOtJvZE8RISX/gy5wYlxFJQKX8TGL9x83As8CNwEp3H+XuDxW09ZzRkzp9\nNjK1a0flenfvpHZ//AEY0RcivVK/VjdA5D/ohLsfqDwxs1lAFzALWAjg7qeAD9OxxMxWA4uIaZA8\n1UmC1zgH9adL24nkpw9wLfEtvEituJUYVqdsnneIm+fDwA+pvq+J0YlG4p6PgennlBQz68+cMuOA\n7e7+DICZLQRed/eTOdcfJUYcAHD3V1O5ycQo02jgeANtrbXot7rP+gK76emzRj+LIURf/tzg9SKl\n08iHyIXhxLfiPN8QQ+f/qjTPvwF4g9hSu8bMhp5ndfuJBGFcpv7BQN5izMo1Q9I1i919q7vvpfmd\nPvuJXUTNxN5DJBkj3H1/1XEop8wgImnAzPoQoyQrC2J8QUyDVHsQ2E4ke00nbDl9NiRzyT7is5iQ\nKdNG7f4YDRx092PNtkOkLBr5EGneADMblh63AXOJIe7N6SbyJrCWmFboJm6gC4C3SmjbUuKGOpdY\nczKFWKdyT7MVuftvZvYKsNzMOolv0k8T6wm8oGgnMWU028yOACOIKYmiMv84diqzHHjBzPoS620G\nA7cDXe6+oUaxfUR/QYxcbXT33QVN2wIsNbPB7t4FYGbXA1e4++m0GHhMQfk8hX3m7ifMbA3wnJkd\nI/pjMdEf1SYSi5NFei0lHyLNm0zPsH43sbhzurt/kna4fAo8Ss88/ffAKuJmkqfWDbWRc2efm9kk\nYB5wp6f/5ZHWJOw0s0fcfVWDMbLmAy8T2zePA8uINRJ50xK4u5vZDGLHzi5iXcI84n+dVOI1kog8\nRuz4KYr9t3rc/Skz+wl4klg38ivQQSRltWwC7kpTZ5e6+6KiBrn7bjPrAO4jFrICzCC26kKMcM2u\n877Oee91+qxiATF6tpn4vXuensQJADMbAEwldsOI9Frm3vCXERG5yJnZ5cAhYL67r7tYYle1Ywqw\nzN1Ht6oNecxsDjDV3SfXvVikhTTyISK5zGwscAOxy+JKYkGtA2//n2MXcfd3086UawrWkrTKKWLK\nTaRXU/IhIvU8TixsPEVsEb2jxMWMrYydy91fbHUbanH3ta1ug0gjNO0iIiIipdJWWxERESmVkg8R\nEREplZIPERERKZWSDxERESmVkg8REREplZIPERERKZWSDxERESmVkg8REREplZIPERERKZWSDxER\nESnVX0UycOWw3iIBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f9202626cc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_bs_decay()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\begin{pmatrix}\\mathbb{0}\\end{pmatrix}$"
      ],
      "text/plain": [
       "[[0]]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "θ = Symbol('theta', real=True)\n",
    "slh.L.substitute({θ: 3*pi/2}).simplify_scalar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "\\begin{align}\n",
       "  \\hat{H} &= \\hat{H}_0 + \\hat{H}_{\\text{int}} + \\hat{H}_{d_1} + \\hat{H}_{d_2}\\\\\n",
       "  \\hat{H}_0 &= - \\frac{g_{1}^{2}}{\\Delta_{1}} \\hat{a}^{(c_{1})\\dagger} \\hat{a}^{(c_{1})} - \\frac{g_{2}^{2}}{\\Delta_{2}} \\hat{a}^{(c_{2})\\dagger} \\hat{a}^{(c_{2})}\\\\\n",
       "  \\hat{H}_{\\text{int}} &= \\frac{g_{1}^{2}}{\\Delta_{1}} \\hat{\\Pi}_{g}^{(q_{1})} \\otimes \\hat{a}^{(c_{1})\\dagger} \\hat{a}^{(c_{1})} + i \\kappa \\hat{a}^{(c_{1})\\dagger} \\otimes \\hat{a}^{(c_{2})} + \\text{c.c.} + \\frac{g_{2}^{2}}{\\Delta_{2}} \\hat{\\Pi}_{g}^{(q_{2})} \\otimes \\hat{a}^{(c_{2})\\dagger} \\hat{a}^{(c_{2})}\\\\\n",
       "  \\hat{H}_{d_1} &= - \\frac{i \\Omega_{1} g_{1}}{2 \\Delta_{1}} \\left(\\hat{\\sigma}_{e,g}^{(q_{1})} \\otimes \\hat{a}^{(c_{1})} - \\text{c.c.}\\right)\\\\\n",
       "  \\hat{H}_{d_2} &= - \\frac{i \\Omega_{2} g_{2}}{2 \\Delta_{2}} \\left(\\hat{\\sigma}_{e,g}^{(q_{2})} \\otimes \\hat{a}^{(c_{2})} - \\text{c.c.}\\right)\\\\\n",
       "\\end{align}"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display_hamiltonian(slh.H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical Values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "controls = sorted([sym for sym in slh.H.all_symbols() if sym.name.startswith('Omega')], key=str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# from algebra import generate_num_vals_code\n",
    "# generate_num_vals_code(slh.all_symbols())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Delta =  5000.0  # MHz\n",
    "g     =    50.0  # MHz\n",
    "kappa =     0.5  # MHz\n",
    "E0    =    70.0  # MHz\n",
    "t0    = -4    # microsec\n",
    "T     =  4    # microsec\n",
    "nt    = 2001\n",
    "theta = 1.5 * np.pi\n",
    "num_vals = {\n",
    "    Symbol('Delta_1', real=True):   Delta,\n",
    "    Symbol('Delta_2', real=True):   Delta,\n",
    "    Symbol('g_2', positive=True):   g,\n",
    "    Symbol('g_1', positive=True):   g,\n",
    "    Symbol('kappa', positive=True): kappa,\n",
    "    Symbol('theta', real=True):     theta,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Magnitude of coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from qnet.algebra.pattern_matching import pattern\n",
    "from qnet.algebra.operator_algebra import ScalarTimesOperator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'-I*Omega_1*g_1/(2*Delta_1)': -0.005*I*Omega_1,\n",
       " '-I*Omega_2*g_2/(2*Delta_2)': -0.005*I*Omega_2,\n",
       " '-I*kappa': -0.5*I,\n",
       " '-g_1**2/Delta_1': -0.500000000000000,\n",
       " '-g_2**2/Delta_2': -0.500000000000000,\n",
       " 'I*Omega_1*g_1/(2*Delta_1)': 0.005*I*Omega_1,\n",
       " 'I*Omega_2*g_2/(2*Delta_2)': 0.005*I*Omega_2,\n",
       " 'I*kappa': 0.5*I,\n",
       " 'g_1**2/Delta_1': 0.500000000000000,\n",
       " 'g_2**2/Delta_2': 0.500000000000000}"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "{str(t.coeff): t.coeff.subs(num_vals) for t in pattern(head=ScalarTimesOperator).findall(slh.H)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time Grid and Guess Pulse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tgrid = QDYN.pulse.pulse_tgrid(t0=t0, T=T, nt=nt) # microsec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p = QDYN.pulse.Pulse(\n",
    "    tgrid, amplitude=(E0 * blackman(tgrid, t0, T)),\n",
    "    time_unit='microsec', ampl_unit='MHz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "controls = OrderedDict([\n",
    "    (Symbol('Omega_1'), p.copy()),\n",
    "    (Symbol('Omega_2'), p.copy()),\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OCT for forward transfer (hs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = qdyn_model.make_qdyn_oct_model(\n",
    "    slh, num_vals, controls, energy_unit='MHz',\n",
    "    mcwf=False, non_herm=True, oct_target='excitation_transfer_fw',\n",
    "    lambda_a=1e-4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model.write_to_runfolder(\"./BS_FB/fw_oct_hs001/\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  OCT for $\\sqrt{SWAP}$ (hs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sqrt_SWAP = QDYN.gate2q.Gate2Q([\n",
    "    [1,             0,             0,   0],\n",
    "    [0, -1/np.sqrt(2),  1/np.sqrt(2),   0],\n",
    "    [0, -1/np.sqrt(2), -1/np.sqrt(2),   0],\n",
    "    [0,             0,             0, -1j]]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = qdyn_model.make_qdyn_oct_model(\n",
    "    slh, num_vals, controls, energy_unit='MHz',\n",
    "    mcwf=False, non_herm=True, oct_target='gate',\n",
    "    lambda_a=1e-4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model.write_to_runfolder(\"./BS_FB/sqrtSWAP_oct_hs001/\")\n",
    "sqrt_SWAP.write(\"./BS_FB/sqrtSWAP_oct_hs001/target_gate.dat\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  OCT for SWAP (hs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = qdyn_model.make_qdyn_oct_model(\n",
    "    slh, num_vals, controls, energy_unit='MHz',\n",
    "    mcwf=False, non_herm=True, oct_target='gate',\n",
    "    lambda_a=1e-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model.write_to_runfolder(\"./BS_FB/SWAP_oct_hs001/\")\n",
    "QDYN.gate2q.SWAP.write(\"./BS_FB/SWAP_oct_hs001/target_gate.dat\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  OCT for PE (hs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = qdyn_model.make_qdyn_oct_model(\n",
    "    slh, num_vals, controls, energy_unit='MHz',\n",
    "    mcwf=False, non_herm=True, oct_target='gate',\n",
    "    lambda_a=1e-2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model.write_to_runfolder(\"./BS_FB/PE_oct_hs001/\")\n",
    "QDYN.gate2q.CPHASE.write(\"./BS_FB/PE_oct_hs001/target_gate.dat\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  OCT for sqrtSWAP (LI, hs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = qdyn_model.make_qdyn_oct_model(\n",
    "    slh, num_vals, controls, energy_unit='MHz',\n",
    "    mcwf=False, non_herm=True, oct_target='gate',\n",
    "    lambda_a=1e-4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model.write_to_runfolder(\"./BS_FB/sqrtSWAP_oct_hsLI001/\")\n",
    "QDYN.gate2q.sqrt_SWAP.write(\"./BS_FB/sqrtSWAP_oct_hsLI001/target_gate.dat\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  OCT for sqrtSWAP (hs, dissipative)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "slh = single_sided_network.network_slh(n_cavity=2, n_nodes=2, topology='bs_fb')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sqrt_SWAP = QDYN.gate2q.Gate2Q([\n",
    "    [1,             0,             0,   0],\n",
    "    [0, -1/np.sqrt(2),  1/np.sqrt(2),   0],\n",
    "    [0, -1/np.sqrt(2), -1/np.sqrt(2),   0],\n",
    "    [0,             0,             0, -1j]]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "num_vals_diss = num_vals.copy()\n",
    "θ = Symbol('theta', real=True)\n",
    "num_vals_diss[θ] = 0.25 * np.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model = qdyn_model.make_qdyn_oct_model(\n",
    "    slh, num_vals_diss, controls, energy_unit='MHz',\n",
    "    mcwf=True, non_herm=True, oct_target='gate',\n",
    "    lambda_a=1e-4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:QDYN.model:Removing 'add_to_H_jump' attribute\n"
     ]
    }
   ],
   "source": [
    "model.write_to_runfolder(\"./BS_FB/sqrtSWAP_diss_oct_hs001/\")\n",
    "sqrt_SWAP.write(\"./BS_FB/sqrtSWAP_diss_oct_hs001/target_gate.dat\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  OCT for forward transfer (hs, dissipative)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = qdyn_model.make_qdyn_oct_model(\n",
    "    slh, num_vals_diss, controls, energy_unit='MHz',\n",
    "    mcwf=True, non_herm=True, oct_target='excitation_transfer_fw',\n",
    "    lambda_a=1e-4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "L = model.lindblad_ops()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:QDYN.model:Removing 'add_to_H_jump' attribute\n"
     ]
    }
   ],
   "source": [
    "model.write_to_runfolder(\"./BS_FB/fw_diss_oct_hs001/\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  OCT for dicke state initialization (hs, dissipative)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model = qdyn_model.make_qdyn_oct_model(\n",
    "    slh, num_vals_diss, controls, energy_unit='MHz',\n",
    "    mcwf=True, non_herm=True, oct_target='dicke',\n",
    "    lambda_a=1e-4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:QDYN.model:Removing 'add_to_H_jump' attribute\n"
     ]
    }
   ],
   "source": [
    "model.write_to_runfolder(\"./BS_FB/dicke_diss_oct_hs001/\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example optimization commands\n",
    "\n",
    "*   hilbert space optimization\n",
    "\n",
    "    ```OMP_NUM_THREADS=1 qdyn_optimize --internal-units=GHz_units.txt --debug --J_T=J_T_re .```\n",
    "\n",
    "*   trajectory optimization (model with `mcwf=True`)\n",
    "\n",
    "    ```OMP_NUM_THREADS=1 qdyn_optimize --internal-units=GHz_units.txt --n-trajs=20 --J_T=J_T_re .```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  },
  "toc": {
   "nav_menu": {
    "height": "123px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}