Asymmetric spectrum in QuTiP

AsymmetricSpectrum.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# Asymmetric spectrum\n",
        "An example that results in an asymmetric spectrum. Note the Hamiltonian is pulled out of thin air and may not be physical. It's motivated by the work of Smart and Swain in the early 90's where asymmetric spectra arise when squeezed light interacts with a two-level atom. Additionally, we use the notion that the squeezing operator can be viewed as fictitious evolution under a Hamiltonian of the form $\\hat{a}^2 - (\\hat{a}^\\dagger)^2$ [Lvovsky](https://arxiv.org/abs/1401.4118). I simply added a weak term to the [Jaynes-Cummings QuTiP example](http://qutip.org/docs/latest/guide/guide-correlation.html?highlight=correlation), the rest is directly from said example.\n",
        "\n",
        "The ultimate goal is simply to show that QuTiP can compute asymmetric spectra using the `spectrum` function.\n",
        "\n*AMC Dawes ([email protected])*"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "from qutip import *\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline"
      ],
      "outputs": [],
      "execution_count": 1,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "N = 4                   # number of cavity fock states\n",
        "wc = wa = 0.1 * 2 * np.pi  # cavity and atom frequency\n",
        "g  = 0.1 * 2 * np.pi       # coupling strength\n",
        "kappa = 0.75            # cavity dissipation rate\n",
        "gamma = 0.25            # atom dissipation rate\n",
        "\n",
        "# Jaynes-Cummings Hamiltonian\n",
        "a  = tensor(destroy(N), qeye(2))\n",
        "sm = tensor(qeye(N), destroy(2))\n",
        "H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag()) + 0.1*(a*a - a.dag()*a.dag()) \n",
        "\n",
        "# collapse operators\n",
        "n_th = 0.25\n",
        "c_ops = [np.sqrt(kappa * (1 + n_th)) * a, np.sqrt(kappa * n_th) * a.dag(), np.sqrt(gamma) * sm]\n"
      ],
      "outputs": [],
      "execution_count": 43,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# calculate the power spectrum using spectrum, which internally uses essolve\n",
        "# to solve for the dynamics (by default)\n",
        "wlist = np.linspace(-1, 1, 200) * 2 * np.pi\n",
        "spec = spectrum(H, wlist, c_ops, a.dag(), a)\n"
      ],
      "outputs": [],
      "execution_count": 44,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# plot the spectra\n",
        "fig, ax = plt.subplots(1, 1)\n",
        "ax.plot(wlist / (2 * np.pi), spec, 'r-', lw=2)\n",
        "ax.legend()\n",
        "ax.set_xlabel('Frequency')\n",
        "ax.set_ylabel('Power spectrum')\n",
        "ax.set_xlim(wlist[0]/(2*np.pi), wlist[-1]/(2*np.pi))"
      ],
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 47,
          "data": {
            "text/plain": [
              "(-1.0, 1.0)"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<matplotlib.figure.Figure at 0x118f527f0>"
            ],
            "image/png": [
              "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XmcHHWd//HXJ5MTc5BjgCQEEkIgCRAJDIgilyJEVwK4KMHVBcFlRVl0XdnFZb1w/Xmt6yobVyMKrheXsgaMBpRDQMAkQIAQAwkhZCA35CTXTD6/P75VTmfSM1PdU9XV3fN+Ph79qOrq6qrP1Mz0p79Hfb/m7oiIiCTRK+8ARESkdihpiIhIYkoaIiKSmJKGiIgkpqQhIiKJKWmIiEhiShoiIpKYkoaIiCSmpCEiIon1zjuAUo0YMcLHjh2bdxgiIjVlwYIF6929sbvHqbmkMXbsWObPn593GCIiNcXMVqRxHFVPiYhIYkoaIiKSmJKGiIgkpqQhIiKJKWmIiEhiShoiIpKYkoaIiCSmpCFSL+bOhcsug82b845E6piShki9+MIX4Ic/hG99K+9IpI4paYjUA3d49tmw/p3vwM6d+cYjdUtJQ6QerFoFmzaF9dWr4ZZb8o1H6paShkg9WLw4LHtHw8mpikoyoqQhUg/iqqn3vx8aGuDxx2HXrnxjkrqkpCFSD+KkMXUqjBoV1l9+Ob94pG4paYjUg7h6avJkOPjgsN7cnF88UreUNETqQVzSmDQJxowJ6ytX5heP1C0lDZFat25deAwcGEoZKmlIhjJNGmY2zcyWmNlSM7umyOuXmNk6M3syenw4y3hE6lJh1ZSZShqSqcymezWzBmAm8A6gGZhnZrPd/dl2u97i7ldmFYdI3XvuubCcODEsVdKQDGVZ0jgRWOruL7j7LuBm4NwMzyfSM61dG5YjR4alkoZkKMukMRooLB83R9va+2sze8rMbjezMRnGI1Kf1q8PyxEjwlLVU5KhLJOGFdnm7Z7fCYx19ynA74AfFT2Q2eVmNt/M5q9bty7lMEVq3IYNYTl8eFgedFC4wW/NGt3gJ6nLMmk0A4Ulh4OBVwp3cPcN7h6PrPZ94PhiB3L3We7e5O5NjY2NmQQrUrPalzQaGnSDn2Qmy6QxD5hgZuPMrC8wA5hduIOZjSx4Oh1YnGE8IvWpfUkD1K4hmcksabh7C3AlMJeQDG5190Vmdp2ZTY92u8rMFpnZQuAq4JKs4hGpW+1LGqB2DclMZl1uAdx9DjCn3bbPFqx/Gvh0ljGI1D2VNKSCdEe4SC1raYGNG6FXL9h//7btKmlIRpQ0RGrZq6+G5dChoQE8FjeEr1pV+ZikrilpiNSyYu0ZAHEvw/h1kZQoaYjUsmLtGdCWRHRfk6RMSUOklqmkIRWmpCFSy+Kk0L6kET/fsAH27KlsTFLXlDREallcPdW+pNGnT+hN1doaeleJpERJQ6SWdVTSgLZEoioqSZGShkgt66ikUbhNjeGSIiUNkVrWUUM4qDFcMqGkIVLLOupyCyppSCaUNERqmUoaUmFKGiK1LElJQ0lDUqSkIVKrWlvhtdfALIw91V5c0lD1lKRISUOkVr32GriH+zF6F5nlQCUNyYCShkitim/aK1bKADWESyaUNERqVZw0hgwp/roawiUDShoitWrTprAsnHypkEoakgElDZFaFZc0OkoagweHMai2boUdOyoXl9Q1JQ2RWtVV9ZRZW2kj7por0k1KGiK1qqvqKVAVlaROSUOkVnVV0gA1hkvqlDREalVXbRqgkoakTklDpFYlqZ5SSUNSpqQhUquSVE+ppCEpU9IQqVWlNISrpCEpUdIQqVVqCJccKGmI1Co1hEsOMk0aZjbNzJaY2VIzu6aT/S4wMzezpizjEakragiXHGSWNMysAZgJvBOYDFxkZpOL7DcIuAp4LKtYROrOnj1tSWPw4I73U0lDUpZlSeNEYKm7v+Duu4CbgXOL7PdF4GuABscRSWrr1jCXxsCBxefSiBU2hLtXJjapa1kmjdHAyoLnzdG2vzCzqcAYd78rwzhE6k+S9gyAvn1DSaS1te09It2QZdKwItv+8lXHzHoB3wT+qcsDmV1uZvPNbP46FbNFkvWciqnbraQoy6TRDIwpeH4w8ErB80HA0cD9ZvYicBIwu1hjuLvPcvcmd29qjBv2RHqyJI3gMTWGS4qyTBrzgAlmNs7M+gIzgNnxi+6+yd1HuPtYdx8LPApMd/f5GcYkUh/KKWmolC4pyCxpuHsLcCUwF1gM3Orui8zsOjObntV5RXqEpG0aoJKGpKqTbhfd5+5zgDnttn22g31PzzIWkbpSSvWUShqSIt0RLlKL1BAuOVHSEKlF5TSEq6QhKVDSEKlFpbRpqKQhKVLSEKlFpVRPqSFcUqSkIVKL1BAuOVHSEKlFKmlITrrschvdoX0tcGi0vwHu7lMyjk1EOlJKSWPIEGhogM2bYedO6Ncv29ikriW5T+OnwNXA08CebMMRkURKaQg3C1VUa9bAhg0walS2sUldS5I01rn77K53E5GKKaV6CkIV1Zo1oV1DSUO6IUnS+JyZ3QD8HtgZb3T3X2YWlYh0bMcO2LUrDHvev3+y96jbraQkSdL4EDAR6ENb9ZQDShoieSgsZVixGQiKUGO4pCRJ0nijux+TeSQikkwp7RkxdbuVlCTpcvtosbm9RSQnpfSciql6SlKSpKTxVuBiM1tOaNNQl1uRPJXaCA4af0pSkyRpTMs8ChFJTiUNyVGSpOFd7yIiFVNOm4ZKGpKSJEnj14TEYUB/YBywBDgqw7hEasfu3fCTn8Bpp8Fhh2V/vnKqp1TSkJR02RDu7se4+5RoOQE4EXgo+9BEasS//ztceilMnQqzK3AfbDnVU+pyKykpecBCd38cOCGDWERqz+rV8I1vhPXNm+G882D+/GzPWU5JY/jwsFy/Hlw1zlK+JAMWfrLgaS/gOEAVoyIAX/wibNsG55wTvvn/+Mdw553Q1JTdOcspafTvDwMHwtatIbmVknBECiQpaQwqePQjtHGcm2VQIjVh40aYNQt69YIvfxnOPz9sf/DB7M8LpSUNUGO4pKLTkoaZNQAD3f3qCsUjUjsefhhaWuAtb4GjjoIDDgjbH3mkbWyoLJRTPQWhMXz58lBFdfjh6cclPUKnJQ13byVUR4lIew88EJannRaWjY0waVIYUHDBguzOW071FKikIalIUj31pJnNNrMPmtl74kfmkYlUuz/8ISxPPbVt2ymn7P1aFsqtnlK3W0lBkqQxDNgAvA04J3q8O8ugRKre1q2hl1SvXnDyyW3b46SRZbtGd6qnQElDuiXJzX03uPvDhRvM7OSOdhbpER55BFpbQy+pQYPatseljoceCl1bkw5dnlRLS0hYZnufNwlVT0kKkpQ0rk+4TaTnaN+eETvkkFBttGlTmCkvbZs3h+XgwaGUUwqVNCQFHZY0zOzNwFuAxnb3agwGGrIOTKSq/fGPYRlXRxUaNw6eeCL0VDrooHTPW24jOLSVNNauTS8e6XE6+6rSFxhISCyF92psBi7IPjSRKvbMM2E5deq+r40bF5Yvvpj+ecttBIe2BLZ6dXrxSI/TYUnD3R8AHjCzm9x9RTkHN7NpwLcIJZMb3P0r7V7/CPAxoBXYClzu7s+Wcy6Rilm7NrQLDBoEY8bs+3qcNJYvT//c5TaCg5KGpCJJpegNZvaXrzVmNtTM5nb1pujGwJnAO4HJwEVFZgD8WTQQ4rHA14D/TB66SE4WLQrLyZOLN3RnmTS6Uz0VJ401a2DPnvRikh4lSdIY4e4b4yfu/hpwQIL3nQgsdfcX3H0XcDPthh9x980FT9+A5u6QWhAnjaOPLv56tZY0+vWDoUNDD6wNG9KNS3qMJEljj5kdEj8xs0NJ9uE+GlhZ8Lw52rYXM/uYmS0jlDSuKnYgM7vczOab2fx16i4oeYuTxlEdTClTiZJGuQMOjhwZlqtWpROP9DhJksa1wENm9mMz+zHwB+DTCd5XrIP6PsnG3We6+3jgX4B/K3Ygd5/l7k3u3tQY9wARyUtXSWPs2LB86aVwL0eaXn01LIcNK+/9cdJQu4aUKckkTL8ljD91C3ArcLy7d9mmQShZFLYSHgy80sn+NwPnJTiuSH7c23pOdVQ9NWBAaD9oaYHm5nTPHyeNeH6MUsXtGippSJm6TBpmZsA04Dh3vxPYz8xOTHDsecAEMxtnZn2BGcBe05qZ2YSCp38FPJ84cpE8rF4Nr70WGqLjb+3FZFVFpZKG5CxJ9dR3gDcDF0XPtxB6RXXK3VuAK4G5wGLgVndfZGbXmdn0aLcrzWyRmT0JfBK4uNQfQKSiCqumOhsipFqThkoa0k1Jxp56k7sfZ2ZPQOg9FZUcuuTuc4A57bZ9tmD946UEK5K7Z6PbiDpqz4hVa9JQSUO6KUlJY3d0z4UDmFkjoE7e0jMtXRqWRxzR+X7VmjRU0pBuSpI0vg3cARxoZl8CHgL+X6ZRiVSr56Nmt65mvjsk6qX+8svpnl8lDclZl9VT7v5TM1sAvD3adJ67L842LJEqFZc0JkzofL8s7ofYsyc0wkN5d4SDShrSbUnHVt6PMH5UL2BAduGIVLHdu0N1kxkcdljn+2bx4bxxY+jyO2QI9E7SHFnE/vuHO8O3bIFt29KLTXqMJF1uPwv8iDCD3wjgRjMrehOeSF1bsSLcrDdmDPTv3/m+w4dDnz7hDu7t29M5f3erpiAkPA1cKN2QpKRxEXCCu3/e3T8HnAT8TbZhiVShpO0ZsPeHc1qljTSSBmgoEemWJEnjRaDwa1U/YFkm0YhUs6TtGbG0P5y7ezd4TO0a0g1JKkZ3AovM7B5Ct9t3EMai+jaAuxcdZFCk7pRS0oD0eyqlVdIYHY0b+kpno/qIFJckadwRPWL3ZxOKSJWrlpJGd5NGPHHUSy917zjSIyXpcvujeN3MhgJj3P2pTKMSqUblljSqNWmsXNn5fiJFJOk9db+ZDTazYcBCQu8pzbAnPcvu3WHObzMYPz7Ze6o1acQ3HippSBmSNIQPiWbYew9wo7sfD5yZbVgiVWbFijDU+cEHd93dNlatSUPVU9INSZJGbzMbCbwPuCvjeESqU6ntGVC9XW5HjQolplWrQglKpARJksZ1hOHNl7r7PDM7DM17IT1Nqe0ZUL0ljT59QuJwVw8qKVmSmftuc/cp7v7R6PkL7v7X2YcmUkXKKWkceGD4Rr9uXaja6q60kgaoikrKlnTsKZGerZySRu/e0NgYvtGvXdv9GLJIGmoMlxIpaYgkUU5JA9KronJX0pCq0GnSMLNeZva+SgUjUpVaWtomU+pqdNv20koamzeHwRIHDQptEt0Vd7tV9ZSUqNOk4e57CPN8i/RccXfbMWNgQIkzA6SVNDZsCMuhQ7t3nJhKGlKmJNVT95jZp8xsjJkNix+ZRyZSLeKqqVLaM2JpJY24TeSAA7p3nJiShpQpydhTl0bLjxVsc6DEcrpIjYobwUttz4D0k8aBB3bvODFVT0mZkow9Na4SgYhUre6UNNK6wW/NmrBMK2k0NkLfvqFxfds2eMMb0jmu1L0kY0/tZ2b/ZmazoucTzOzd2YcmUiWqoaQRJ420qqd69YJx0ffBF15I55jSIyRp07gR2AW8JXreDPx7ZhGJVJvnngvLeqqegraBF+OSlEgCSZLGeHf/GrAbwN23A5ZpVCLVYteu0N22V6/uNYSvXh3utShX2tVT0PbzLNNEnJJckqSxy8wGEBq/MbPxhNn8ROrfsmXh/oixY6Ffv9LfP2AADBkSBgaMb84rR9rVU6CShpQlSdL4PPBbYIyZ/RT4PfDPWQYlUjX+/OewPPLI8o+RRhVVFtVTKmlIGZIMWHg3YS6NS4CfA03ufn+2YYlUiSVLwjLvpJFFSSNOGippSAmS9J76MSFpLHP3u9x9fdKDm9k0M1tiZkvN7Joir3/SzJ41s6fM7Pdmdmhp4YtkrBqSRly11asXDB9efhztjR0bjvnSS6HtRiSBpL2nRgLXm9kyM/uFmX28qzeZWQMwE3gnMBm4yMwmt9vtCULJZQpwO/C1kqIXyVo1JI24aqqxERoayo+jvb59w01+e/aEqWxFEkhSPXUv8CXgM8ANQBNwRYJjn0iYuOkFd98F3Ayc2+7Y97n769HTR4GDS4hdJFvu6bRpdPcGv7SHECmkxnApUZLqqd8DDwMXAkuAE9x9YoJjjwYKB7ZpjrZ15DLgNwmOK1IZ69fDa6/BwIFtpYVydLekkUV325gaw6VESaqnniLc3Hc0MAU4OuqC25Vi93IU7ahuZh8glGC+3sHrl5vZfDObv27dugSnFklBYdWUdePWpMJ7NcqRZdJQSUNKlKR66h/d/VTgfGADoY1jY4JjNwNjCp4fDOwzIbGZnQlcC0x396L3f7j7LHdvcvemxsbGBKcWSUEa7RmQXptGFtVT6kElJepywEIzuxI4BTgeWAH8EHgwwbHnARPMbBzwMjADeH+7Y08FvgdMc/cU5sMUSdGzz4blpEndO041V0/FCXHx4vSPLXUpydDoA4D/BBa4e0vSA7t7S5Rw5gINwA/dfZGZXQfMd/fZhOqogcBtFor/L7n79FJ/CJFMPP10WB5zTPeOM2RIuDN861bYsiXMvleKLJPGhAlhJsDly0N8Awemfw6pK0mGRv+6mb0R+Ej0wf6guy9McnB3nwPMabftswXrZ5YWrkgFPfNMWB59dPeOYwajR4cqoObm0ksuWVZP9ekTShvPPBNKGyeckP45pK4k6T11FfBT4IDo8RMz+4esAxPJ1YYNoTrpDW9oG0K8O7ozU17cgJ5FSQPgqKPCctGibI4vdSVJ9dSHgTe5+zYAM/sq8AhwfZaBieQqrpo66qhw13R3dSdpxO8ZM6bz/cqlpCElSPLfYEBrwfNWNDS61Lu0qqZi8Qd+c3Np79uyJdwr0q9fuCM8C/HPqKQhCSQpadwIPGZmd0TPzwN+kF1IIlUgrUbwWLkljXj/Qw7p3r0inYlLGnGiFOlEkobw/zSz+4G3EkoYH3L3J7IOTCRX8QdoWknj4GiEnFKTxksvheUhh6QTRzHjx4eSzMqVsHkzDB6c3bmk5nWYNMysP/AR4HDgaeA7pXS5FalZ7uknjXJLGpVIGg0NMHEiLFwY7k056aTsziU1r7M2jR8RhvZ4mjBS7X9UJCKRvL30UvjG3diYXjfXcts0KpE0oK1dQ1VU0oXOksZkd/+Au38PuAA4tUIxieRr/vywPO649I65//6h++6WLbBpU/L3VSppHHtsWC5YkO15pOZ1ljR2xyuqlpIeJU4aTU3pHdOsvCqqSiWN+GeNf3aRDnSWNN5oZpujxxZgSrxuZpsrFaBIxc2bF5Zp3x1dTmN4pZLGcceFxLZwIewsOm6oCNBJ0nD3BncfHD0GuXvvgnV1r5D6tGdP27fttJNGqSWNPXva2kCyurEvNnhwGE5k9+627sYiRaRwq6tIHVm2LLQ5jBwJo0ale+xSG8PXrAkf4o2NYcDDrMVVVHFJS6QIJQ2RQllVTUHpJY1KVU3F4p9Z7RrSCSUNkUJZVU1B24f/ihXJ9q900lBJQxJQ0hApFH9gptlzKhbPkvfcc8n2j5NLpZLGsceGG/0WLYJt2ypzTqk5ShoisR07QtIwgxNPTP/4hx4ahut4+eUw4VFX4uRy2GHpx1LMfvvBlCmhAf7RRytzTqk5ShoiscceC91NjzkGhg1L//gNDaWVNuLpZuMBBSvh9NPD8oEHKndOqSlKGiKxP/whLE87LbtzHHFEWC5Z0vl+7vkkjfhnv//+yp1TaoqShkgs/nadZdI48siw7CpprF4d5tEYOjS7GfuKOeWUUD332GOwfXvlzis1Q0lDBGDXLvjjH8P6Kadkd544aXRVPRWXMiZPzm4ejWKGDQvtGrt2qV1DilLSEIEwUN/27TBpUnoj2xaTtKRRmDQqTe0a0gklDRFoq8PPsmoK9i5puHe8Xx7tGbE4adx7b+XPLVVPSUMEYM6csDzzzGzPM2wYDB8euty+8krH+8XzdedV0mhoCNV1GzdW/vxS1ZQ0RF59NXxA9ukD73hH9ufrqorKPd+ksf/+oV2ntRXmzq38+aWqKWmIzJ0bbmg75ZTKzI89cWJYxomhvXXrQiIbMiT9QROTeve7w/Kuu/I5v1QtJQ2R+IPxr/6qMueL7zZ/+OHir8e9lqZMqWzPqUJx0pgzB1o0B5u0UdKQnq21FX7727Aef1BmLe7S++CDxRvD77svLM84ozLxFHPEEeHu9VdfVddb2YuShvRsDz0UPhgPP7ztbu2sTZoUGsNfeQWWL9/39bjXUp5JwwzOOSes//KX+cUhVSfTpGFm08xsiZktNbNrirx+qpk9bmYtZnZBlrGIFPWzn4Xle99buXOawVvfGtYffHDv19atg6eegv794aSTKhdTMe97X1jefHMokYmQYdIwswZgJvBOYDJwkZm17wryEnAJ8LOs4hDp0K5dcNttYf3976/suTtKGvENdW95S0gceXrTm2D8eFi1SmNRyV9kWdI4EVjq7i+4+y7gZuDcwh3c/UV3fwrYk2EcIsXNnRvGdzrmGDj66MqeO27XeOihvbfHVVNve1tl4ynGrC2Z/vSn+cYiVSPLpDEaKJzXsjnaJlId4qqpSpcyAI47LsxfsWRJ2zhUu3fDr38d1vNszygUX5tf/EIDGAqQbdIo1lewk3ETOjmQ2eVmNt/M5q9bt66bYYkA69fDHXeEb9MzZlT+/H36tH0gf/GLYXnTTWGK1yOPDFVD1WDixDCL4ebNcOuteUcjVSDLpNEMjCl4fjDQybgJHXP3We7e5O5NjY2NqQQnPdwPfhAmXHrnO2Hs2HxiuPZa6N07lHiefLIteXzhC2EYj2pxxRVhOXNmvnFIVcgyacwDJpjZODPrC8wAZmd4PpFkWlvhf/4nrH/sY/nFMXYsfOhD4W70qVNh5crQvlLJnlxJXHRRGDNr3jz405/yjkZyllnScPcW4EpgLrAYuNXdF5nZdWY2HcDMTjCzZuC9wPfMrINxFURS9Otfw4oVYe7tadPyjeUznwlDhfTqFcZ8+uY3w3o1GTAALrssrF9/fb6xSO7MOxueuQo1NTX5/Pnz8w5DapU7nHwyPPIIfOMb8MlP5h1RiCmv4UKSWr483ADZq1douB83Lu+IpERmtsDdm7p7nCr7SiOSsXvuCQlj+HC4/PK8owmqPWFASBJ/8zdhHKqvfCXvaCRHShrSc7jD5z8f1q++GgYOzDWcmnPttaGkceONoZeX9EhKGtJz3HlnWykjzwbwWnXkkaF78u7d8LnP5R2N5ERJQ3qGHTvgE58I65/5jEoZ5frCF8I9JjfdFHpTSY+jpCE9w9e/Hhpzjz5apYzuOPxw+Md/DOtXXRW6C0uPoqQh9e+ZZ+BLXwrr118fbqiT8l17LRx4YJhnY9asvKORClPSkPq2c2fo9bNzZ7jX4PTT846o9g0e3Ha/xqc+VXxOEKlbShpS366+OsxPMX48/Nd/5R1N/Xjve8N8G9u2wcUXa0rYHkRJQ+rX978fvhH36ROG9lbjd7pmzoSDDgpzgvzrv+YdjVSIkobUp7lz4aMfDeuzZlXPqLH1ZMSIMPJtQ0PoaBBPaCV1TUlD6s/998N554Uqk099Ci65JO+I6tcpp4SEAfDBD+47E6HUHSUNqS9z5sC73hXuy7j8cvja1/KOqP594hOhVLdzJ0yfDgsW5B2RZEhJQ+rHjTeGD63t2+HSS8Pw57UwrlOtM4Nvfxve8x7YuBHe/nZ47LG8o5KMKGlI7du9O4xWe+mlYa6Ma66BG26oviHG61lDA/z85yFxbNoUEsddd+UdlWRA/1VS2158EU47LcxD0adPaPT+8pdVwshD375w882hbWPbNjj33NDeUWPTL0jnlDSkdt1yC7zxjWEQwoMPhnvvhb/7u7yj6tn69IEf/SiMJrxnD/zzP4cqw9Wr845MUqKkIbVn1Sp4//vDiKubN8P558PChfDWt+YdmUAo5X3uc/CrX4XZCO+6CyZNCvOyq9RR85Q0pHa0tIQG14kTQ/35gAHw3e/CL34R5rCW6jJ9Ojz5JJx9dmgg//CH4YwzwlhgUrOUNKT67dkTEsNxx8HHPx5KF+ecA88+C3//92q/qGaHHgq/+U24I7+xER54AKZMCaXExYvzjk7KoKQh1au1NdxlfOyxcMEF8PTT4UPoV7+C2bNh7Ni8I5QkzEJ14uLFcOWVod3jllvgqKPgwgvDDYGqtqoZShpSfZqbw2Q/Y8eGQfGefhpGj4b//m/4859DtYfUnuHDw1hgS5fCFVeEIepvvRVOPTWUPmbOhPXr845SumBeYxm+qanJ58+fn3cYkrYNG+D//g9uvx3uvrttcp/x4+Gf/incg9GvX74xSrpWroTvfS/cU7NmTdjW0ABnnhlKINOnh0QjqTCzBe7e1O3jKGlILtxDg+jdd8Nvfwv33ReqoyB8Az3//NBeccYZukmv3u3aBXfcEe7o/93v2v4OevWCE0+EadPgrLPg+OPDvSBSFiUNqS07d8ITT4R7Kh55BB56KHSdjTU0hLuIL7ggJIwRI/KLVfKzfj388pehLeuBB8Ld/rH+/eGEE+Dkk8PjpJP0d1ICJQ2pXps3h1LEU0+F9ognnwyD2O3cufd+I0eGb5Bnnx2WqoqQQlu3hhLob34TRi4u1ttq9OjQUeKNbwyPY44JVZoqkexDSUPytW0bLFu27+O552DFiuLvmTQJ3vzmtsfkyeouK8lt2BBKqQ8/HEqqjz8Or7++7369esEhh8CECXs/xo2DMWNg0KDKx14FlDQkG3v2hH/OV16Bl18Oy/brzc2wdm3Hx+jXLySEY44JjylTQrXC0KGV+zmk/rW2hi8qCxe2PRYtCl9a4o4UxQwZEpLHmDFh+JkxY0Kp94AD9n7U2UyPShrSMfcwPPjmzXs/Nm4MCWHDhlB3XGz9tdc6/4eL9e0bvrmNH9/2OPzwtmXv3tn/nCLF7NoFL7wQuvY+/3zbY8WK0GNrx45kx9lvv7YEMmJE+NKz//5hGT/aPx8yJCSbhoZsf8YypJU0Mv3PNrNpwLeABuAGd/8UBKmnAAALKElEQVRKu9f7Af8LHA9sAC509xezjKkqtLaGP+xdu8If8PbtoZj9+ut7r7d/3v6111/fOyls2dK2HvdAKcfQoTBqVHiMHl18feTIqvzHEKFv3zDUzMSJ+77mHr4crVwZSswrV4bH2rWh2+/atW3rr78eRlF+8cXSYxgwICSPpI8BA0JDf7xMst67dy7Vu5klDTNrAGYC7wCagXlmNtvdny3Y7TLgNXc/3MxmAF8FLuz0wLt3h6qSlpa9H62t+27raHsp+8bbd+9u+6DfubP4emevFa535wM9qf79YfDgvR9DhoTG5vgxYsS+68OGqZQg9css/K2PGAFTp3a8n3toiI+TyPr1oRS+cWNYxo/C5xs3hse2beEL3vbtsG5ddj9Lr157J5N+/ULCjB/tn6cks+opM3sz8Hl3Pzt6/mkAd/9ywT5zo30eMbPewGqg0TsJqsnM66Jyql+/tl/qfvu1PQYM6Ph5sfX2iWHw4NDQp94jIvmIq4e3bk322LIl1DjEtQ5J11taSgrLoOqrp0YDKwueNwNv6mgfd28xs03AcKDjsQR69w51jL17tz0aGvZ+3tn27uzbPnt3tN7Vfg0N6jUkUq/M2r7cHXBAdudpaWlLJDt2dF3Tcd55qZw2y6RR7FOxfQkiyT6Y2eXA5QCHHHJIx106RUR6it6929pEKijL8RmagTEFzw8GXulon6h6agjwavsDufssd29y96bGxsaMwhURka5kmTTmARPMbJyZ9QVmALPb7TMbuDhavwC4t7P2DBERyVdm1VNRG8WVwFxCl9sfuvsiM7sOmO/us4EfAD82s6WEEsaMrOIREZHuy7RvpbvPAea02/bZgvUdwHuzjEFERNKjMadFRCQxJQ0REUlMSUNERBJT0hARkcRqbpRbM9sCLMk7jgRG0Nmd7dVDcaanFmIExZm2WonzSHfv9mQitTgy3ZI0xk/JmpnNV5zpqYU4ayFGUJxpq6U40ziOqqdERCQxJQ0REUmsFpPGrLwDSEhxpqsW4qyFGEFxpq1HxVlzDeEiIpKfWixpiIhITqoyaZjZe81skZntMbMOeyWY2TQzW2JmS83smoLt48zsMTN73sxuiUbZzSLOYWZ2T3See8xsaJF9zjCzJwseO8zsvOi1m8xsecFrx+YVZ7Rfa0Esswu2Z349E17LY83skehv4ykzu7DgtUyvZUd/awWv94uuzdLoWo0teO3T0fYlZnZ2mnGVEecnzezZ6Pr93swOLXit6O8/pzgvMbN1BfF8uOC1i6O/k+fN7OL2761wnN8siPE5M9tY8FpFrqeZ/dDM1prZMx28bmb27ehneMrMjit4rfRr6e5V9wAmAUcC9wNNHezTACwDDgP6AguBydFrtwIzovXvAldkFOfXgGui9WuAr3ax/zDCaL77Rc9vAi6owPVMFCewtYPtmV/PJDECRwATovVRwCpg/6yvZWd/awX7fBT4brQ+A7glWp8c7d8PGBcdpyHHOM8o+Pu7Io6zs99/TnFeAvx3kfcOA16IlkOj9aF5xdlu/38gjOZd6et5KnAc8EwHr78L+A1h0ruTgMe6cy2rsqTh7ovdvasb+E4Elrr7C+6+C7gZONfMDHgbcHu034+AdOY53Ne50fGTnucC4Dfu/npG8XSk1Dj/ooLXs8sY3f05d38+Wn8FWAtUYlauon9r7fYpjP924O3RtTsXuNndd7r7cmBpdLxc4nT3+wr+/h4lTI5WaUmuZ0fOBu5x91fd/TXgHmBalcR5EfDzjGLpkLv/gSKT1xU4F/hfDx4F9jezkZR5LasyaSRUbA7y0YQ5xje6e0u77Vk40N1XAUTLriYEnsG+f1RfioqM3zSzflkESfI4+5vZfDN7NK5Co3LXs6RraWYnEr79LSvYnNW17Ohvreg+0bWK57tP8t5KxlnoMsI30Fix338Wksb519Hv83Yzi2cBrcrrGVXzjQPuLdhcqevZlY5+jrKuZW53hJvZ74CDirx0rbv/KskhimzzTraXpbM4SzzOSOAYwqRUsU8DqwkffrOAfwGuyzHOQ9z9FTM7DLjXzJ4GNhfZr6zrmfK1/DFwsbvviTandi2LnbLItqTz3af699iFxOcysw8ATcBpBZv3+f27+7Ji769AnHcCP3f3nWb2EUIp7m0J35uWUs41A7jd3VsLtlXqenYl1b/N3JKGu5/ZzUN0NAf5ekLxq3f0ja/Y3OSJdRanma0xs5Huvir6IFvbyaHeB9zh7rsLjr0qWt1pZjcCn8ozzqjKB3d/wczuB6YCvyCl65lGjGY2GPg18G9RUTs+dmrXsohS5rtvtr3nu0/y3krGiZmdSUjUp7n7znh7B7//LD7kuozT3TcUPP0+8NWC957e7r33px5h27mS/u5mAB8r3FDB69mVjn6Osq5lLVdPFZ2D3EMLz32E9gMIc5AnKbmUo3CO867Os099Z/ThGLcbnAcU7f2Qgi7jNLOhcZWOmY0ATgaereD1TBJjX+AOQv3sbe1ey/Jadme++9nADAu9q8YBE4A/pRhbSXGa2VTge8B0d19bsL3o7z/HOEcWPJ0OLI7W5wJnRfEOBc5i79J7ReOMYj2S0JD8SMG2Sl7PrswG/jbqRXUSsCn6klXetaxE636pD+B8QhbcCawB5kbbRwFzCvZ7F/AcIXtfW7D9MMI/5lLgNqBfRnEOB34PPB8th0Xbm4AbCvYbC7wM9Gr3/nuBpwkfcD8BBuYVJ/CWKJaF0fKySl7PhDF+ANgNPFnwOLYS17LY3xqh+mt6tN4/ujZLo2t1WMF7r43etwR4Z8b/O13F+bvofyq+frO7+v3nFOeXgUVRPPcBEwvee2l0nZcCH8ozzuj554GvtHtfxa4n4cvoquh/o5nQVvUR4CPR6wbMjH6GpynokVrOtdQd4SIiklgtV0+JiEiFKWmIiEhiShoiIpKYkoaIiCSmpCEiIonV4hzhIiUzs1ZCd8PYee7+Yk7hiNQsdbmVHsHMtrr7wE5ej+94F5FOqHpKeiwLczbcZmZ3AndH2642s3nRQHlfKNj3WgvzKvzOzH5uZp+Ktt9v0ZwvZjbCzF6M1hvM7OsFx/r7aPvp0XtuN7M/m9lPo7vYMbMTzOyPZrbQzP5kZoPM7EErmBvEzB42symVukYi7al6SnqKAWb2ZLS+3N3Pj9bfDExx91fN7CzCMB8nEu6inW1mpwLbCENITCX8zzwOLOjifJcRhms4IRpO4mEzuzt6bSpwFGH8n4eBk83sT8AtwIXuPi8aY2s7cANhbolPmNkRhLvxn+rWlRDpBiUN6Sm2u3ux2fzucfd4LoKzoscT0fOBhCQyiDDY5OsAlmwWtrOAKWYWj9k1JDrWLuBP7t4cHetJwjAzm4BV7j4PwN03R6/fBnzGzK4mDPlwU9IfWCQLShrS020rWDfgy+7+vcIdzOwTdDxkdAtt1bz92x3rH9x9rwHgzOx0wphqsVbC/6EVO4e7v25m9xAm0nkfYSwukdyoTUOkzVzgUjMbCGBmo83sAOAPwPlmNsDMBgHnFLznReD4aP2Cdse6wsz6RMc6wsze0Mm5/wyMMrMTov0HRUOsQ6ii+jYwr6BUJJILlTREIu5+t5lNAh6J2qa3Ah9w98fN7BbCqLArgAcL3vYfwK1m9kH2nrXtBkK10+NRQ/c6Opkm1913mdmFwPVmNoDQnnEmYZ7pBWa2GbgxpR9VpGzqcitSIjP7POHD/D8qdL5RhMlxJnrbTIUiuVD1lEgVM7O/BR4jzOWghCG5U0lDREQSU0lDREQSU9IQEZHElDRERCQxJQ0REUlMSUNERBJT0hARkcT+P8sNFjGqvwVxAAAAAElFTkSuQmCC\n"
            ]
          },
          "metadata": {}
        }
      ],
      "execution_count": 47,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "outputs": [],
      "execution_count": 32,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "outputs": [],
      "execution_count": null,
      "metadata": {
        "collapsed": false,
        "outputHidden": false,
        "inputHidden": false
      }
    }
  ],
  "metadata": {
    "kernel_info": {
      "name": "python3"
    },
    "kernelspec": {
      "name": "python3",
      "language": "python",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python",
      "version": "3.6.3",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "nteract": {
      "version": "0.8.4"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}