Drag_hd.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T13:42:07.420528Z",
     "start_time": "2020-05-05T13:42:02.587285Z"
    }
   },
   "outputs": [],
   "source": [
    "import qutip as qt\n",
    "import numpy as np\n",
    "import matplotlib.pylab as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T13:47:02.552384Z",
     "start_time": "2020-05-05T13:47:02.547398Z"
    }
   },
   "outputs": [],
   "source": [
    "from qutip import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 283,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T15:40:24.384397Z",
     "start_time": "2020-05-05T15:40:24.380408Z"
    }
   },
   "outputs": [],
   "source": [
    "N=2\n",
    " # Define time vector\n",
    "wq=6\n",
    "alpha=0.2\n",
    "adag_q1=create(N)\n",
    "a_q1=destroy(N)\n",
    "H0=wq*adag_q1*a_q1+alpha/2*adag_q1*adag_q1*a_q1*a_q1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 284,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T15:40:25.043850Z",
     "start_time": "2020-05-05T15:40:25.039896Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\\begin{equation*}\\left(\\begin{array}{*{11}c}0.0 & 0.0\\\\0.0 & 6.0\\\\\\end{array}\\right)\\end{equation*}"
      ],
      "text/plain": [
       "Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\n",
       "Qobj data =\n",
       "[[0. 0.]\n",
       " [0. 6.]]"
      ]
     },
     "execution_count": 284,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 285,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T15:40:27.752840Z",
     "start_time": "2020-05-05T15:40:27.747918Z"
    }
   },
   "outputs": [],
   "source": [
    "def projector(N,n,m):\n",
    "    # |m><n| operator\n",
    "    a = np.zeros((N,N),dtype = np.float64)\n",
    "    if(n >= N or n < 0 or m >= N or m < 0):\n",
    "        return Qobj(a)\n",
    "    a[n,m] = 1\n",
    "    return Qobj(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 286,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T15:40:28.115648Z",
     "start_time": "2020-05-05T15:40:28.111659Z"
    }
   },
   "outputs": [],
   "source": [
    "psi0 = ket2dm(basis(N,0)) # Define initial state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 330,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T16:13:11.927773Z",
     "start_time": "2020-05-05T16:13:11.918798Z"
    }
   },
   "outputs": [],
   "source": [
    "def H1_coeff(t, Hargs):\n",
    "    w = Hargs['w']\n",
    "    A = Hargs['A']\n",
    "    B = Hargs['B']\n",
    "    T = Hargs['T']\n",
    "    D = Hargs['D']\n",
    "    q = Hargs['q']\n",
    "    return A * np.exp(-B * (t / T - 0.5) ** 2) +1j* A * np.exp(-B * (t / T - 0.5) ** 2)*(B*(T-2*t))/T/T/D*q\n",
    "\n",
    "H1 = (a_q1+adag_q1) \n",
    "H1=[H1, 'A * sin(w * t) * exp(-B*(t/T-0.5)**2)-0*A*cos(w*t)*( B*(T-2*t)* exp(-B*(t/T-0.5)**2)/T/T)/D*q']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 346,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T16:15:46.390140Z",
     "start_time": "2020-05-05T16:15:46.263507Z"
    }
   },
   "outputs": [],
   "source": [
    "# H = [H0,[H1,H1_coeff]]\n",
    "# for wq in np.linspace(5,6,100):\n",
    "H = [H0,H1]\n",
    "q=1\n",
    "tlist = np.linspace(0, 60, 100)\n",
    "Hargs = {'w': wq,'A': 0.14,'B':10,\"T\": 40, \"D\": alpha,\"q\":q}\n",
    "opts = Options(rhs_reuse=True)\n",
    "# eops=[projector(N,0,0),projector(N,1,1),projector(N,2,2),projector(N,3,3)]\n",
    "eops=[projector(N,0,0),projector(N,1,1)]\n",
    "output = mesolve(H, psi0, tlist, c_ops=[],e_ops= eops, args=Hargs,options=opts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 347,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T16:15:47.826355Z",
     "start_time": "2020-05-05T16:15:47.394006Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x26875187320>"
      ]
     },
     "execution_count": 347,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4VOX99/H3d2Yy2QMhJCxJSFjCvgQIyFZREUSL4FJF3LDaYm2tPq31qb3a2v6q/dnFR22ttVKl7riACyKKilpB2cIe9h0CBEJYE7LNzP38cUYNkJAhmeTMTL6v65przjlz5pzvgcknJ/c5c99ijEEppVRkcdhdgFJKqeDTcFdKqQik4a6UUhFIw10ppSKQhrtSSkUgDXellIpAGu5KKRWBNNyVUioCabgrpVQEctm147Zt25rs7Gy7dq+UUmFpxYoVh40xqfWtZ1u4Z2dnk5+fb9fulVIqLInI7kDW02YZpZSKQBruSikVgTTclVIqAmm4K6VUBNJwV0qpCFRvuIvIDBE5JCIFdbwuIvJ3EdkmImtFZFDwy1RKKXU+Ajlzfx4Yf47XLwdy/I9pwNONL0sppVRj1HufuzHmCxHJPscqk4AXjTVe3xIRaS0iHYwxB4JU42mW7zrCwi3Fpy8UQawnHP5ph0NwOgSnWM9RTsHldOByCG6XA7fTgdvlICbKSbT/Oc7tJCHaRXy0izi3ExFpikNQKriMAZ8HvFXgqQRvNfiqrXmvxz9dDT6vtZ7PA8YLxmctM8aa/vpBzXlTY94/JOdpy78eprPmtH/+zBpre63OYT4DGP4znIcI7TQc0no26S6C8SWmdGBvjflC/7Kzwl1EpmGd3dOpU6cG7Wzl7qM8+dm2b+ab6v/X7XSQHB9FcpybtKQY0lvH0LFVLJlt4uiSGk/ntvEkxkQ1zc5Vy+LzQVkxlB789nGqxP84AhXHoaoUKk9CVRlUl1sPTzl4qsBTQUBhqELHdx8Li3Cv7fS21k+aMWY6MB0gLy+vQZ/GO0d35c7RXWt9zeczGMBnjPXwgcfnw+szVHsNHp+Pao+h2uejymM9Kj0+Kqq9lFd7Ka/yUlrpobTSw/Hyao6UVlFSVsWhkxVs2H+cw6VVp+0vLTGaHu0T6d4ukZ7tE+mb3oqctARcTr1OrWrweeH4XijZDkd2WNPH9sLxQjixH0qLrLPpMzndENsGYltDdCK4EyChHUTFQVQMuGLAFW09O6PB5bbe44zyP7vB4bLmHVHgcFrzDieI0//s8D+c1p++It8uo8Z8bdPgf65l+pv5mup67VzvqU+Y/nUdk9TkuwhGuBcCmTXmM4D9QdjueXM4rP9o52n/4c6gbb+i2sveI6fYcbiMHcVlbDtUypaDJ3l5yW4qPT4Aol0O+nRMYliXFIZ3TSEvqw2x7uDVoMLAkR2waxHsX2U9Dm4Ab+W3rzvd0CoDktIhexQkdbQeCe0gsT3Ep1oPd3wDwk4pSzDCfQ5wt4i8BlwAHG+q9na7xUQ5yWmXSE67xNOWe32GnYfLKNh3nHX7jrNyz1Gmf7GDf36+HbfLwdhe7bhqYDqju6fidulZfcTxeWHXQtg4F7YvsMIdILoVdOgPQ38IbbtDSjdo08UKcYd+DlTTqjfcRWQmcBHQVkQKgd8BUQDGmH8B84ArgG3AKeD7TVVsqHI6hG5pCXRLS+CqgekAlFV6WL7rCJ9tOsTctQd4f90B2sS7uW1ENlNHZNMqVtvrw97+1bBmJqx/22onj4qDzhfCBXdB14uhTVcNcWUbMTZdcc7LyzMtpVfIaq+PRVsP8/KS3SzYdIjEGBffH5HNtNFdSYi2rWNO1RAVx2Hdm7DiBShaazWx5IyDftdB98sgKtbuClWEE5EVxpi8etfTcG9eBfuO849Pt/Hh+iLSEqP59Xd7MXFAR73tMtQdWAv5z8HaN6D6FLTrB4OnQr/vQWyy3dWpFkTDPcSt3nuMB98tYG3hcS7o3IY/Xt2PbmkJdpelajIGtn0Cix6H3V+CK9YK87zboeNAvdipbKHhHga8PsPry/fy5w83UV7t5d4xOUy7sAtReiulvYyBTXPh8z/BwQJIyoBhd8HAm/QsXdlOwz2MFJ+s5HdzCpi3rog+HZN49LoB9OrQ9PfBqlqUbId591t3vbTtDqN+Bn2/Z91DrlQICDTc9RQxBKQmRvPPmwbzr5sHcfBEJRP/sYinP9+O16ffOmw2Ph988Vf45zDYuwzG/xnuWgy5N2qwq7Ckt2qEkPF9OzC0cwq/fnsdf/5wEws2HuTxyblktomzu7TIVnEC3r4TNs+D3lfB5X+2vkykVBjTM/cQ0ybezT9vGsQTk3PZfPAkE55cxH/P7ChNBU/JdnhuLGyZb52tX/e8BruKCBruIUhEuGpgOnN/OooOrWK47T/L+MenW/FpM01wFRXAs5daX0C65W0Y9iO9A0ZFDA33EJaVEs/bPx7JxAEdefSjLfxi1hpthw+WogJ44Uqr060fLIAuo+2uSKmg0jb3EBfrdvLE5Fy6pibw2MdbcIjwl2v7f9NJmmqAg+vhxYlWsN82F1Jq72VUqXCm4R4GRIR7xuTgM4YnPtmKU4RHrumnAd8QJdvhhYlWF7ka7CqCabiHkXvH5OD1GZ78dBtJsS5+/d3edpcUXk4WwUtXAQamvqfBriKahnsYERF+PrY7x05V8++FOxmclcz4vh3sLis8VByHl6+FshK47T1o283uipRqUnpBNcyICL+Z0IsBGa24/8217C4ps7uk0FddATNvhOJNMPklSB9sd0VKNTkN9zAU7XLyjxsHIQI/fmUlFdVeu0sKXT4vvD0Ndi+Cq/4F3cbYXZFSzULDPUxltonjsetzWb//BI/M22h3OaHJGPjwV7DhXRj3R+h/nd0VKdVsNNzD2KW923HbiGxeWLyb5buO2F1O6PnyCVj2DAy/G0bcbXc1SjUrDfcwd/9lPUhvHcsvZ6/V5pma1s2CT35v9eg49iG7q1Gq2Wm4h7n4aBf/e00/dhSX8dRn2+wuJzTsXgzv3AVZI+Gqf+o4pqpF0k99BBjdPZVrBqbz9Ofb2XjghN3l2KtkO7x2I7TuBJNfBle03RUpZQsN9wjx2wm9aRUbxQNvrWu5HYyVH4NX/BdNb3wD4trYW49SNtJwjxDJ8W5+M6EXa/YeY+byPXaX0/yMgffugWO74YZX9dunqsXTcI8gV+WmM6xLG/78wSYOl1baXU7zWvmCdcvjJb+BrOF2V6OU7TTcI4iI8PBVfSmv9vLIvE12l9N8Dm2CDx6ALhfBiHvtrkapkKDhHmG6pSXyw+90YfbKQpbuKLG7nKZXXQGzbgd3PFz9jN4Zo5Sf/iREoJ9ekkNGciy/eaeAKo/P7nKa1ld/h0Pr4ep/6fB4StWg4R6BYt1O/jCpD1sPlfLvhTvsLqfpHC+EhY9B70mQM9buapQKKRruEeqSnu24ol97/r5gK7sOR2jPkZ/8HoxPv4GqVC003CPY767sQ5TTwW/fLcCYCLv3fc9SWPcmjLwHkrPsrkapkKPhHsHaJcVw/2U9WLj1MHPW7Le7nODx+eDDX0JiBxj5f+yuRqmQFFC4i8h4EdksIttE5IFaXu8kIp+JyCoRWSsiVwS/VNUQNw/LYkBGKx6au4Hj5dV2lxMcBbNh/yq49PcQnWB3NUqFpHrDXUScwFPA5UBvYIqInDl452+AN4wxA4EbgH8Gu1DVME6H8Mer+1FSVsXjH2+xu5zG8/lg4aOQ2gv6XW93NUqFrEDO3IcC24wxO4wxVcBrwKQz1jFAkn+6FRBBbQDhr296K26+IIsXF++iYN9xu8tpnM3vW8Plfec+vaddqXMI5KcjHdhbY77Qv6ym3wM3i0ghMA/4aVCqU0Hzi3E9SI5z8+C7BeHbsZgx8MWjkNwZ+lxtdzVKhbRAwl1qWXZmOkwBnjfGZABXAC+JyFnbFpFpIpIvIvnFxcXnX61qsFZxUTxweU9W7jnGrBWFdpfTMNsXwIHVMOpn4HTZXY1SIS2QcC8EMmvMZ3B2s8sdwBsAxpjFQAzQ9swNGWOmG2PyjDF5qampDatYNdi1gzIYnJXMnz7cRGmlx+5yzt8X/w+S0mHAFLsrUSrkBRLuy4EcEeksIm6sC6ZzzlhnDzAGQER6YYW7npqHGIdDeHBCb46UVTFj0U67yzk/uxfDnq9gxD3gcttdjVIhr95wN8Z4gLuB+cBGrLti1ovIH0Rkon+1+4AfisgaYCZwm4m4b81EhgGZrRnbux3//mIHx05V2V1O4L76O8SlwKBb7a5EqbAQ0O0Gxph5xpjuxpiuxpg/+pc9aIyZ45/eYIwZaYwZYIzJNcZ81JRFq8b5+djunKz0hE+/M4e3wuYPYMgPwB1ndzVKhQW9l6wF6tUhiQn9O/CfL3eFx6Aei58CpxuG/NDuSpQKGxruLdTPxnanotrL059vt7uUcys7DGtmwoDJkKAX4ZUKlIZ7C9U1NYFrBmXw0pLdoX32vvw58FTA8LvtrkSpsKLh3oL9aHRXqjw+Xl++t/6V7VBdAcv/DTmXQWoPu6tRKqxouLdg3dISGNWtLS8v2Y3HG4IjNq19HcqKYYSetSt1vjTcW7hbh2dx4HgFH284aHcpp/N5rdsfOwyA7O/YXY1SYUfDvYUb06sd6a1jeWHxLrtLOd2m96Fkm9Vfu9TWA4ZS6lw03Fs4p0O4ZXgWS3YcYVPRCbvLsRgDix63OgjrfWYHpEqpQGi4KybnZRLtcvDi4t12l2LZtRD2r4QRPwWH0+5qlApLGu6K5Hg3k3I78vbKfZyoCIHRmhY9AfFpkHuT3ZUoFbY03BUAN12QRXm1l/fXHrC3kANrra59h/0IomLsrUWpMKbhrgDon9GKbmkJ9vf1vuhxcCdA3h321qFUmNNwVwCICNcOymDF7qPsPFxmTxHFW2D92zD0hxDb2p4alIoQGu7qG1cPTMch8NZKm87eFz0GrhgY9hN79q9UBNFwV99o3yqGUTmpvLVyX/OPs3pkJ6x9A/Ju1w7ClAoCDXd1mmsHpbPvWDlLdpY0744XPQ4Ol3X7o1Kq0TTc1Wku69OexGgXs1fsa76dHtsLq1+FQbdAUofm269SEUzDXZ0mJsrJd/t34IOCA5Q11yDaXz0JGKurAaVUUGi4q7NcOziDU1VePiwoavqdlZXAyheh/2Rondn0+1OqhdBwV2fJy0qmU5s4ZjfHXTPLnwVPOYy4p+n3pVQLouGuziIiXDMoncU7Sth3rLzpdlRdDsumW4NxpPVsuv0o1QJpuKtaXTMwA2PgnVVNeGF19atw6jCMvLfp9qFUC6XhrmrVKSWOodltmL2yEGOa4J53n9e6kJo+GLJGBH/7SrVwGu6qTtcOTmdHcRmr9x4L/sY3zYWjO622dh2MQ6mg03BXdbq8XweiXQ7eWtkETTNLnrYG4+h1ZfC3rZTScFd1S4qJ4rI+7ZmzZj+VHm/wNly8BfYshrzv62AcSjURDXd1TtcMSud4eTWfbDgUvI2uetHqamDAlOBtUyl1Gg13dU7fyUklvXUsM5ftCc4GPVWweib0uBwS0oKzTaXUWTTc1Tk5HcLkIZks2naY3SVB6Od98zzr9sdBUxu/LaVUnTTcVb2uz8vE6RBmLtvb+I2tfBGSMqDrJY3fllKqThruql7tW8VwSc80Zq3YS5XH1/ANHdsD2z+FgTfrhVSlmlhA4S4i40Vks4hsE5EH6ljnehHZICLrReTV4Jap7Hbj0E4cLq3i4w0HG76RVa9YzwNvCk5RSqk61RvuIuIEngIuB3oDU0Sk9xnr5AC/AkYaY/oA2ndrhLmweyMvrBoDa16FLhdB607BLE0pVYtAztyHAtuMMTuMMVXAa8CkM9b5IfCUMeYogDEmiPfNqVBQ88LqroYMoL1vhdUs0++64BenlDpLIOGeDtS8klboX1ZTd6C7iHwpIktEZHxtGxKRaSKSLyL5xcXFDatY2WbykExcDuGVpbvP/80Fb4HTDT2/G/zClFJnCSTca+v448yepFxADnARMAV4VkRan/UmY6YbY/KMMXmpqToIcrhplxTDZX3b8/ryvZyqOo9Rmnw+WP82dLsUYs/6WCilmkAg4V4I1BwiJwPYX8s67xpjqo0xO4HNWGGvIsz3R2RzosLD2+fTFfDeJXByP/S5pukKU0qdJpBwXw7kiEhnEXEDNwBzzljnHeBiABFpi9VMsyOYharQMDgrmb7pSTz/5a7AuwIueAtcMdCj1tY6pVQTcNW3gjHGIyJ3A/MBJzDDGLNeRP4A5Btj5vhfGyciGwAvcL8xpqQpC1f2EBGmDs/m/llr+Wp7CSO7tT33G7we2PAO5IyD6MTmKVKpelRXV1NYWEhFRYXdpdQpJiaGjIwMoqKiGvT+esMdwBgzD5h3xrIHa0wb4Of+h4pwVw7oyCMfbOL5r3bVH+67F0FZMfTVJhkVOgoLC0lMTCQ7OxsJwfEEjDGUlJRQWFhI586dG7QN/YaqOm8xUU6mDM3kk40H2Xvk1LlXLngLouKtcVKVChEVFRWkpKSEZLCD9RdySkpKo/6y0HBXDXLzsCwcIry85By3RXqrYeN7Vlu7O675ilMqAKEa7F9rbH0a7qpBOrSK5dJeaby5orDugTx2LYTyI9Dn6uYtTiml4a4a7qYLsjhSVsWHBUW1r7DhXatJptulzVuYUkrDXTXcqG5tyUqJ45UltfQ34/VYTTLdL4Oo2OYvTqkQ99BDD9GzZ0/Gjh3LlClTePTRR4O6/YDullGqNg6HcOPQTjzywSY2F52kR/satzru/hJOlUCfq+wrUKkA/M9769mw/0RQt9m7YxK/u7JPna/n5+cze/ZsVq1ahcfjYdCgQQwePDioNeiZu2qU7w3OwO108OqZ/c1seAei4qDbWHsKUyqELVq0iEmTJhEbG0tiYiJXXnll0PehZ+6qUVISormiX3veWrmPX17ekzi3C3xeq0kmZ5zeJaNC3rnOsJtKwN/ubgQ9c1eNdtOwLE5Wenhvjb/Lod1fWV9c0iYZpWo1atQo3nvvPSoqKigtLeX9998P+j403FWj5WUl0zU1ntkr/Z2JbXgXXLHWmbtS6ixDhgxh4sSJDBgwgGuuuYa8vDxatWoV1H1ouKtGExGuHpjOsp1HKCw5CRvnQM6l4I63uzSlQtYvfvELNm/ezDvvvMPmzZv1gqoKTZNyrfFbVnz+LpQe1BGXlKrHtGnTyM3NZdCgQVx77bUMGjQoqNvXC6oqKDLbxDEkO5nYTf/GRCch2peMUuf06quvNun29cxdBc01fZMZUbWYo9lXQFSM3eUo1aJpuKuguTJmDQlSwVy+Y3cpSrV4Gu4qaBI2z6bEmcZT29Pw+pr+Pl6lVN003FVwlBbDtgUc6zqJg6XVLN6uA3EpZScNdxUc698C4yVj9FQSo128s/o8BtBWSgWdhrsKjrWvQ7t+RKf3Y2zvdny0vogqj8/uqpRqsTTcVeMd3AD7VsCAyQBc0a8DJyo8fLX9sM2FKRW6Xn75ZYYOHUpubi533nknXm8dg940kN7nrhovfwY43TDgRgBG5bQlIdrFB+uKuKhHms3FKVWPDx6AonXB3Wb7fnD5n+p8eePGjbz++ut8+eWXREVF8eMf/5hXXnmFW2+9NWglaLirxqkshTWvQe+rID4FsAbQvrRXGvM3FPGwty9RTv0DUamaFixYwIoVKxgyZAgA5eXlpKUF90RIw101TsFsqDoJQ+44bfHl/Trwzur9LN5ewoXdU20qTqkAnOMMu6kYY5g6dSqPPPJIk+1DT6lUwxkD+c9BWm/IvOC0l0Z3TyXe7eSDggM2FadU6BozZgyzZs3i0KFDABw5coTdu3fX867zo+GuGm7/SjiwBvJuB5HTXoqJcnJJr3bMX38Qj1fvmlGqpt69e/Pwww8zbtw4+vfvz9ixYzlwILgnQtosoxoufwZExUP/ybW+/N1+7XlvzX6W7DjCqJy2zVycUqFt8uTJTJ5c+89OMOiZu2qY8qOwbjb0+x7EJNW6yujuacRGOXl/nTbNKNXcNNxVw6x6GTzlMOQHda4S63Yyplca89cXUa1NM0o1Kw13df58Xlg2HTqNgA79z7nqhP4dOVJWxVfa14xSzUrDXZ2/LfPh2B644M56V72oRyoJ0S7mfj14tlIhwpjQ7rm0sfUFFO4iMl5ENovINhF54BzrfU9EjIjkNaoqFdqW/guS0qHnhHpXjYlyMq53O+ZrXzMqhMTExFBSUhKyAW+MoaSkhJiYhg96U+/dMiLiBJ4CxgKFwHIRmWOM2XDGeonAPcDSBlejQt+hTbDzvzDmQXAGdrPVhAEdeGvVPhZuLWZMr3ZNXKBS9cvIyKCwsJDi4mK7S6lTTEwMGRkZDX5/ID+dQ4FtxpgdACLyGjAJ2HDGeg8BfwF+0eBqVOhb9gw4o2HQbQG/ZVS3VFrFRjF37QENdxUSoqKi6Ny5s91lNKlAmmXSgb015gv9y74hIgOBTGPM3CDWpkJN+TGrH5l+133Tj0wg3C4H4/u056P1RVRUB7fnO6VU7QIJd6ll2TcNVSLiAB4H7qt3QyLTRCRfRPJD+c8hVYe1b0D1KRha9+2PdZkwoANlVV4+33yoCQpTSp0pkHAvBDJrzGcANW99SAT6Ap+LyC5gGDCntouqxpjpxpg8Y0xeaqp2JhVWjIEV/4GOA63HeRreJYWUeDdz9K4ZpZpFIOG+HMgRkc4i4gZuAOZ8/aIx5rgxpq0xJtsYkw0sASYaY/KbpGJlj8LlcGgDDL6tQW93OR1MzO3IJxsOcaSsKri1KaXOUm+4G2M8wN3AfGAj8IYxZr2I/EFEJjZ1gSpE5P8H3AnQ99oGb2LykEyqvD7eXqXjqyrV1AK6l80YMw+Yd8ayB+tY96LGl6VCSvlRawDsAVMgOrHBm+nZPonczNa8vnwPt4/MRqS2yzlKqWDQb6iq+q19AzwVDW6SqemGIZlsOVjKyj3HGl+XUqpOGu7q3IyBFc/7L6TmNnpzEwZ0JM7t5PXlexpfm1KqThru6twK860LqYOmBmVzCdEuJg7oyHtrDnCyojoo21RKnU3DXZ3bmlfBFduoC6lnmjwkk/JqL3PXaj/vSjUVDXdVt+oKawDsXhPqHJCjIXIzW9OjXSKvLN0dsh03KRXuNNxV3bZ8CBXHrbtkgkhEmDoim4J9J1i8Q/t5V6opaLiruq2ZCYkdoMtFQd/0NYPSaZvgZvoXO4K+baWUhruqS2kxbP0Y+l8PDmfQNx8T5eS2Edl8vrmYTUUngr59pVo6DXdVu4JZYLxBb5Kp6eZhWcS5nXr2rlQT0HBXtVv9KnTIhbReTbaL1nFuJg/JZM7q/Rw4Xt5k+1GqJdJwV2c7uB6K1kLujU2+qztGdcYAMxbtbPJ9KdWSaLirs61+FRyuoN7bXpeM5Dgm9O/Aq0v3cOyU9hapVLBouKvTeT1WXzI5l0F822bZ5V0XdaWsysuML3c1y/6Uagk03NXptn8KZYeapUnmaz3bJ3FZn3Y8/+VOTmiXBEoFhYa7Ot2aVyG2DeSMa9bd3n1xDicqPLy0eHez7lepSKXhrr5VfhQ2zbMGwHa5m3XX/TJacXGPVJ5duIOySk+z7lupSKThrr5V8BZ4KyG36e5tP5e7L8nh6KlqXl2q3QEr1Vga7upba2ZCai/r/nYbDM5KZmS3FJ75YgcV1V5balAqUmi4K8vhbdYg2LlTwMbh735yUTcOl1by1kodZ1WpxtBwV5YV/wFxQv/JtpYxvGsK/TNaMf2L7Xh92h2wUg2l4a6g6hSsegl6T4TE9raWIiL8aHRXdpWcYv76IltrUSqcabgrWPeG1W/70Gl2VwLAZX3a07ltPE9/vl0H81CqgTTcWzpjYNm/oV1f6DTc7moAcDqEaRd2Yd2+43y1XQfzUKohNNxbuj2L4WCBddZu44XUM109MJ3UxGj+9d/tdpeiVFjScG/plj4DMa2tLy6FkJgoJ7eP7MzCrYdZv/+43eUoFXY03FuyE/th43sw6BZwx9ldzVluvKAT8W4nzy3U7oCVOl8a7i3ZkqfB+CDvDrsrqVWr2CiuH5LJnDX7KTpeYXc5SoUVDfeW6mSRdSG1//XQprPd1dTp9pGd8RnD81/tsrsUpcKKhntLtfAx8FbB6F/aXck5ZbaJ4/K+HXhl6W5KtUMxpQKm4d4SHdtrfSN14M2Q0tXuaur1g+905mSFhzeW77W7FKXChoZ7S/TFX6znC++3t44ADeyUTF5WMjO+3InH67O7HKXCQkDhLiLjRWSziGwTkQdqef3nIrJBRNaKyAIRyQp+qSooSrbDqldg8Pehdabd1QRs2oVdKDxazty1B+wuRamwUG+4i4gTeAq4HOgNTBGR3mestgrIM8b0B2YBfwl2oSpIPvkdON3wnZ/bXcl5ubRXO7q3S+Cpz7bh0w7FlKpXIGfuQ4Ftxpgdxpgq4DVgUs0VjDGfGWNO+WeXABnBLVMFxdZPrPvaL7zP9g7CzpfDIfzk4m5sPVTKRxsO2l2OUiEvkHBPB2peySr0L6vLHcAHtb0gItNEJF9E8ouLiwOvUjVedQV8cD+kdIMR99hdTYN8t18HslLieOqzbdqhmFL1CCTca+twpNafLBG5GcgD/lrb68aY6caYPGNMXmpqauBVqsb76kk4sgMu/wu4ou2upkFcTgd3je7Kun3H+WLrYbvLUSqkBRLuhUDNK28ZwP4zVxKRS4FfAxONMZXBKU8FxdFdsPBR6D0Juo2xu5pGuXpQOu2TYnjq0212l6JUSAsk3JcDOSLSWUTcwA3AnJoriMhA4BmsYD8U/DJVgxkDH/zSGmXpskfsrqbRol1O7hzdhWW7jrBwqzbtKVWXesPdGOMB7gbmAxuBN4wx60XkDyIy0b/aX4EE4E0RWS0ic+rYnGpuG9+DLR/Cxb+CVue6VBI+pgztRKc2cfzhvQ1637tSdRC7Lkzl5eWZ/Px8W/bdYlScgKeGQlxbmPY5OF12VxQ089cXcedLK/ifiX2YOiLb7nKUajYissIYk1ffevoN1Uj22R+tDsKu/FtEBTvAuN7tGNkthcc+3sLRsiq7y1Eq5Gi4R6p9K6yBOIb8ADIG211N0IkID07ow8mP8gQRAAAM+UlEQVSKap74ZIvd5SgVcjTcI5HPC3N/BgntYMxv7a6myfRon8hNF2Tx8tI9FOzT0ZqUqknDPRItfxYOrIHxj0BMK7uraVL3jetOakI098xcRZl2CazUNzTcI83JIvj0Yeh6CfS52u5qmlzrODePT85lZ0kZv5uz3u5ylAoZGu6RZv6vwVMJVzwKUtuXiyPP8K4p3H1xN2atKOTd1fvsLkepkKDhHkm2fwYFs2DUz8JiEI5gundMDoOzkvn12wXsKC61uxylbKfhHimqy+H9n0NyZyvcWxiX08Hfbsglyinc8UK+3h6pWjwN90jx+Z+sjsGufAKiYuyuxhYZyXH8+9Y89h0t586XVlDp8dpdklK20XCPBAfWWL0+5t4MXS6yuxpb5WW34a/X9WfZriM8MHuddg2sWqzI+tpiS+T1wJyfQlwKjHvI7mpCwqTcdPYeOcWjH20hIzmW+8b1sLskpZqdhnu4W/KUdeZ+3fMQ18buakLGTy7uRuHRcp78dBtpSTHcMkyH9VUti4Z7ONu30rqnvecE6H2V3dWEFBHh4av6UnyykgffLSA1IZrxfcNraEGlGkPb3MNV+VF4cyrEp8HEJ1vMPe3nw+V08OSNAxmQ0Zp7XlvF8l1H7C5JqWaj4R6OjIF3fgwn9mtzTD3i3C5m3DaEjNax3PH8cjYXnbS7JKWahYZ7OPrqSdg8D8Y+BJlD7K4m5LWJd/PiHUOJdTu5dcZSCo+esrskpZqchnu4WT0TPn4Qek2EYXfZXU3YyEiO44Xbh3KqysutM5ZxRL/kpCKchns4WfsGvHMXdL4Qrpmu7eznqWf7JJ6bOoR9R8u55bmlHDulAa8il4Z7uFg3C96+E7JHwZTXICrW7orC0tDObXjmlsFsPVTKTc9qwKvIpeEe6nw++O9fYfYPIHMY3Pg6uOPsriqsXdQjjem3DGbrwVJu1jN4FaE03ENZxXF47Ub47GHofz3cPBvc8XZXFREu6pHGM7cOZktRKdc+/ZX2JKkijoZ7qCrMh+kXwbaP4fK/wtXP6Bl7kF3cI40Xbh/K0VPVTHrqSz7ffMjukpQKGg33UOOpsr51+txY8FbD1LlwwTS9eNpEhndN4d2fjCQjOY7vP7+cxz7eQkW19iapwp+GeygpWgfPjoEv/goDpsBdX0LWcLuriniZbeKYfddwrspN5+8LtjL+iS/475Ziu8tSqlE03EOBp9I6W59+EZw8AJNfhqv+GfGDW4eSOLeLxyfn8tIdQxERps5Yxg9eWM6qPUftLk2pBhG7+rvOy8sz+fn5tuw7pOxaBO/fB8WbrLP1y/5XuxOwWaXHy7MLdzL9ix0cL69mWJc23Dm6K6NzUnE4tHlM2UtEVhhj8updT8PdJicOwMe/hXVvQqtOMOExyBlrd1WqhrJKDzOX7eHZhTspOlFBVkocN13QiesGZ5Ic77a7PNVCabiHqrISqw/2pc9YF0xH3muNeap3woSsKo+P+euLeGnJbpbtPILb6eCyvu25YUgmw7uk6Nm8alYa7qHmyA7InwHLZ0D1Keg9ES79PbTpYndl6jxsLjrJzGV7eGtlIScqPGQkxzJxQEcm5abTo32i3eWpFkDDPRQc3Q1bP4K1r0PhchAH9L0WvvMLSOtpd3WqESqqvXxYUMTslYV8tb0Er8+Qk5bAxT3TuDAnlbzsZGKinHaXqSJQUMNdRMYDfwOcwLPGmD+d8Xo08CIwGCgBJhtjdp1rmxEX7p5KOLQB9q+2voC06ws4tsd6La2P9Q3Tft+DVhn21qmC7nBpJfPWHeDDgiLydx2lyuvD7XLQq0MS/dNb0Tc9iS6pCWSnxNM2wY3odxZUIwQt3EXECWwBxgKFwHJgijFmQ411fgz0N8b8SERuAK42xkw+13bDNtw9VXB0J5Rsg+LNcGij9SjeBL5qa53YZMgaafXe2Hm0nqW3IKeqPCzZUcLi7SWs23ecgn0nKK30fPN6QrSLjORYOraOpWPrGDq0iqV9UgztW8WQlhhNWmIMSbEu/QWg6hRouAcyhupQYJsxZod/w68Bk4ANNdaZBPzePz0L+IeIiLGrzac+xlgXM71V4KmA6nLruaoUKkuh8iRUHIOyw3DqMJwsguOFcHyv9Wx8324rKQPSekG3S6BDLnTMheTO+o3SFirO7eKSnu24pGc7AHw+w54jp9hZUsauw9Zj37EK9h0rZ8Xuoxwvrz5rG26ngzbxblrHRdEm3k1STBQJMS4Sol3EuZ3ERjmJiXISE+XA7XIQ7XLidjlwOx1Eff3sFKKcDlxfPzsEl8OB0yk4RXA4sJ5FcDgEh4BDBKnxLNSYx/pI6y+d8BFIuKcDe2vMFwIX1LWOMcYjIseBFOBwMIo8zcqXrJGIvmGssP5m2vftw+cDn8f/qAbv18/n0Qug0w0J7aBVptUrY/8sSMmBlG7Qtpt+0Uidk8MhZLeNJ7ttPPQ4+/XyKi9FJyooOl5BcWklh05Yz0fLqjhSVs3RU1XsOFxKWaWXExXVlFd58fjsPWeygv/raakxbf1CsGa+Xafm+76ZrvHq6cvr2mf9v1TC6dfObyf05vohmU26j0DCvbZ/szM/XYGsg4hMA6YBdOrUKYBd1yIuxTpTPn3D35bgcFrT4gCHCxwOEKcV0s4oa5kr2pp2Rlv9oruiwRVr9bgYnQDuBKtpJS4FohP1LFw1mVi3k85t4+ncNvDePqu9PsqrvVRW+6jy+qis9lLl9VHtMVR5vVR5DB6fj2qvjyqPwev7et7g8xk8PoPXGIyxXvP6f1n4jMHnP1fy+V83xvpB9tWY/voPcmu+5jRnLf+GqXWSmn/c1/V3fiC/ykK0jaBOXVKbvnfXQMK9EKj5KyYD2F/HOoUi4gJaAWcNNW+MmQ5MB6vNvSEF0/MK66FUCxXldBDldECM3ZWoUBZI3zLLgRwR6SwibuAGYM4Z68wBpvqnvwd8GrLt7Uop1QLUe+bub0O/G5iPdSvkDGPMehH5A5BvjJkDPAe8JCLbsM7Yb2jKopVSSp1bIM0yGGPmAfPOWPZgjekK4LrglqaUUqqhtMtfpZSKQBruSikVgTTclVIqAmm4K6VUBNJwV0qpCGRbl78iUgzsbuDb29IUXRvYQ48l9ETKcYAeS6hqzLFkGWNS61vJtnBvDBHJD6RXtHCgxxJ6IuU4QI8lVDXHsWizjFJKRSANd6WUikDhGu7T7S4giPRYQk+kHAfosYSqJj+WsGxzV0opdW7heuaulFLqHMIu3EVkvIhsFpFtIvKA3fWcDxGZISKHRKSgxrI2IvKxiGz1PyfbWWMgRCRTRD4TkY0isl5E7vUvD8djiRGRZSKyxn8s/+Nf3llElvqP5XV/d9chT0ScIrJKROb658P1OHaJyDoRWS0i+f5lYff5AhCR1iIyS0Q2+X9mhjfHsYRVuPsH634KuBzoDUwRkd72VnVengfGn7HsAWCBMSYHWOCfD3Ue4D5jTC9gGPAT//9DOB5LJXCJMWYAkAuMF5FhwJ+Bx/3HchS4w8Yaz8e9wMYa8+F6HAAXG2Nya9wyGI6fL4C/AR8aY3oCA7D+f5r+WMw3w2mF/gMYDsyvMf8r4Fd213Wex5ANFNSY3wx08E93ADbbXWMDjuldYGy4HwsQB6zEGiP4MODyLz/tcxeqD6xR0hYAlwBzscaeDLvj8Ne6C2h7xrKw+3wBScBO/Nc3m/NYwurMndoH6063qZZgaWeMOQDgf06zuZ7zIiLZwEBgKWF6LP6mjNXAIeBjYDtwzBjj8a8SLp+zJ4D/C/j88ymE53GANXTqRyKywj/2MoTn56sLUAz8x99c9qyIxNMMxxJu4R7QQNyqeYhIAjAb+D/GmBN219NQxhivMSYX68x3KNCrttWat6rzIyITgEPGmBU1F9eyakgfRw0jjTGDsJpgfyIiF9pdUAO5gEHA08aYgUAZzdScFG7hHshg3eHmoIh0APA/H7K5noCISBRWsL9ijHnLvzgsj+VrxphjwOdY1xFa+wd7h/D4nI0EJorILuA1rKaZJwi/4wDAGLPf/3wIeBvrl244fr4KgUJjzFL//CyssG/yYwm3cA9ksO5wU3Nw8alY7dchTUQEa9zcjcaYx2q8FI7Hkioirf3TscClWBe8PsMa7B3C4FiMMb8yxmQYY7Kxfi4+NcbcRJgdB4CIxItI4tfTwDiggDD8fBljioC9ItLDv2gMsIHmOBa7Lzg04ALFFcAWrHbRX9tdz3nWPhM4AFRj/Ua/A6tddAGw1f/cxu46AziOUVh/3q8FVvsfV4TpsfQHVvmPpQB40L+8C7AM2Aa8CUTbXet5HNNFwNxwPQ5/zWv8j/Vf/5yH4+fLX3cukO//jL0DJDfHseg3VJVSKgKFW7OMUkqpAGi4K6VUBNJwV0qpCKThrpRSEUjDXSmlIpCGu1JKRSANd6WUikAa7kopFYH+P7KKGvmD3zh/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(tlist,output.expect[0],label='g')\n",
    "plt.plot(tlist,output.expect[1],label='e')\n",
    "plt.plot(tlist,output.expect[2],label='f')\n",
    "# plt.plot(tlist,output.expect[3],label='d')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 341,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T16:15:18.633152Z",
     "start_time": "2020-05-05T16:15:18.627285Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9994618394934353"
      ]
     },
     "execution_count": 341,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "output.expect[1][-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 336,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T16:14:47.863331Z",
     "start_time": "2020-05-05T16:14:47.335939Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2687da585f8>]"
      ]
     },
     "execution_count": 336,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAD8CAYAAABkbJM/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8VFX6x/HPk0mnBFLooSZ0qaGLIB1UsICIDZUVG2t3Rf2thdVd24q9UMSKiFhARSMqoIACoROKJNRQA4EACenn98e9uDEGE8gkNzPzvF+veWXm3nNnvleHPLnlnCPGGJRSSqmS+DkdQCmllGfQgqGUUqpUtGAopZQqFS0YSimlSkULhlJKqVLRgqGUUqpUtGAopZQqFS0YSimlSkULhlJKqVLxdzqAO0VGRprGjRs7HUMppTzKqlWrDhtjokpq55aCISJDgJcAFzDNGPN0kfUXAC8C7YCrjDFzCq3LBzbYL3cbY4bby5sAs4BwYDVwnTEm569yNG7cmISEBHfsklJK+QwR2VWadmU+JSUiLuA1YCjQGhgjIq2LNNsN3ADMLOYtThljOtiP4YWWPwNMNsbEAkeBcWXNqpRS6ty54xpGVyDJGLPdPgKYBYwo3MAYs9MYsx4oKM0biogA/YDTRyLvApe6IatSSqlz5I6CUR/YU+h1ir2stIJFJEFEfhWR00UhAjhmjMk7x/dUSinlZu64hiHFLDubMdMbGmP2iUhT4EcR2QAcL+17ish4YDxAw4YNz+JjlVJKnQ13HGGkANGFXjcA9pV2Y2PMPvvndmAR0BE4DNQQkdMF7YzvaYyZYoyJM8bERUWVeJFfKaXUOXJHwVgJxIpIExEJBK4C5pVmQxGpKSJB9vNIoBewyVizOi0ERtpNxwJz3ZBVKaXUOSpzwbCvM0wA4oHNwGxjTKKITBKR07fIdhGRFGAU8JaIJNqbtwISRGQdVoF42hizyV73IHCviCRhXdOYXtasSimlzp140xStcXFxRvtheLY9aZnsOpLJ/vRTHDqRTWigi7phwdQJC6FlnWoEB7icjqiU1xGRVcaYuJLaeVVPb+WZkg6d4Kv1+/lmwwG2HjxxxnahgS76tazFsPPq0q9lLS0eSlUwLRjKMTsOZ/Bc/BbmbziACHRpHM6jF7emdb3q1A0Lpnb1YDKy89ifnkXK0VP8tC2V+I0H+Gr9fmpVC+Kegc0Z1bkB/i4dEk2piqCnpFSFO5GVy3PxW5m5fDeB/n78rXdTru3WkFrVg0vcNr/AsCz5MJMX/Mbq3ceIqVWVxy5pTe9YvUNOqXOlp6RUpbRxbzoTZq5mz9FTjOkazZ39Y6lVreRCcZrLT+gdG8X5MZHEJx7k2W+3cN30FdzWtxn3DWyuRxtKlSMtGKpCGGP4YPlu/vXlJsKrBPLx+O7ENQ4/5/cTEYa0rUPfFlFM+moTbyxKZuWONF65uiN1w0LcmFwpdZr+OabKXUGB4YkvN/HPLzbSMyaC+Xf1LlOxKCw4wMW/LzuPl67qwOb9xxn+6lJ++4sL50qpc6cFQ5Wr/ALDPz5dzzvLdvK385vw9tguhFcJdPvnjOhQn8/v6IUAV771C+tTjrn9M5TydVowVLnJySvgzo/WMGdVCncPiOWRi1rh51fc0GPu0bx2Nebc2pOqQf5cPXU5K3akldtnKeWLtGCoclFQYLj/k3V8vWE//3dRK+4e0Bxr1Pry1TAilE9u7UGt6kHcOGMFG/eml/tnKuUrtGCocvFs/FbmrdvHA4Nb8LfeTSv0s+uGhfDRzd0JCwngxndWsicts0I/XylvpQVDud17v+zkzcXJXNu9Ibf3beZIhtrVg3n3pq5k5+Zzw4wVHMv8y9l9lVKloAVDudUPmw/y2LxEBrauzRPD21bIaagzia1djanXx7En7RQ3v5dATl6pJnxUSp2BFgzlNruOZHD3x2tpWy+Ml6/qiKscL3CXVremETw3qh0rdx7lP99sdjqOUh5NC4Zyi6zcfG77YDV+Irx+TSdCAivPwIAjOtTnhp6NmbF0J1+tL/XcXkqpIrRgKLd4fF4im/YfZ/Lo9kSHhzod508eHtaKTg1r8OCc9SSnnnQ6jlIeSQuGKrNPV6Uwa+UeJlwYQ7+WtZ2OU6xAfz9eu6YTQQEubvtgFVm5+U5HUsrjaMFQZZJyNJPH5iXStUk49wxs7nScv1Q3LITJozvw28GTPBe/1ek4SnkctxQMERkiIltFJElEJhaz/gIRWS0ieSIystDyDiLyi4gkish6ERldaN07IrJDRNbajw7uyKrc53TnPID/jmpfKS5yl6RP8yiu79GI6Ut2sCz5sNNxlPIoZS4YIuICXgOGAq2BMSLSukiz3cANwMwiyzOB640xbYAhwIsiUqPQ+geMMR3sx9qyZlXu9fbSHfy6PY1HL2ldKa9bnMnEoS1pElmFBz5Zz/GsXKfjKOUx3HGE0RVIMsZsN8bkALOAEYUbGGN2GmPWAwVFlv9mjNlmP98HHAJ0JhwPsO3gCZ6N38qAVrUZ1bmB03HOSmigP/+9sj37008x6ctNTsdRymO4o2DUB/YUep1iLzsrItIVCASSCy1+yj5VNVlEgsoWU7lLfoHhgTnrqRrkz9NXnOdo57xz1alhTW7vG8OcVSks3HrI6ThKeQR3FIzifluc1byvIlIXeB+40Rhz+ijkIaAl0AUIBx48w7bjRSRBRBJSU1PP5mPVOfpw+S7W7jnGoxe3JrKq59bxv/ePoVlUFf75xUYyc/KcjqNUpeeOgpECRBd63QAode8oEakOfA38nzHm19PLjTH7jSUbmIF16utPjDFTjDFxxpi4qCg9m1XeDqRn8ey3W+kdG8mIDvWcjlMmQf4u/nN5O1KOnuLF77c5HUepSs8dBWMlECsiTUQkELgKmFeaDe32nwPvGWM+KbKurv1TgEuBjW7IqsrosXkbySso4KlLPfNUVFFdm4Qzpms005fsIHGfDoWu1F8pc8EwxuQBE4B4YDMw2xiTKCKTRGQ4gIh0EZEUYBTwlogk2ptfCVwA3FDM7bMfisgGYAMQCTxZ1qyqbL5LPEB84kHu6t+chhGec1dUSSYOaUXN0EAe+mwD+QVndTZVKZ8ixnjPP5C4uDiTkJDgdAyvdConnwEvLKZasD9f/v18Alze1edz3rp93PnRGp68tC3Xdm/kdBylKpSIrDLGxJXUzrv+1aty88biZPYeO8UTw9t4XbEAuKRdXbo3Def577ZyNEPnzlCqON73L1+53Z60TN5cnMwl7evRrWmE03HKhYjw+PA2nMjK44UFvzkdR6lKSQuGKtGTX2/CJcLDw1o6HaVctaxTneu6N+LD5bvYtO+403GUqnS0YKi/9PO2VOITDzKhXwx1w0KcjlPu7hnQnLCQAB6fl4g3Xd9Tyh20YKgzyssvYNKXm2gYHsq485s4HadChIUG8MDglqzYmcZX6/c7HUepSkULhjqjjxP2sO3QSR4e1pLggMozg155G90lmpZ1qvFs/Bay83TeDKVO04KhinUyO4/JC36jS+OaDG5Tx+k4FcrlJzw8rBV70k7x3rJdTsdRqtLQgqGK9eaiZA6fzOHhYa28okf32bqgeRR9mkfxyo/b9DZbpWxaMNSf7E8/xdSft3NJ+3p0bFjT6TiOeXhYK05m5/HKj0lOR1GqUtCCof7k+fjfMAb+MbiF01Ec1aJONUZ3ieb9X3ey83CG03GUcpwWDPUHm/cf57M1KYzt2cijZtErL/cMaI6/nx/Pf6dzgCulBUP9wfPxW6ka5M8dF8Y4HaVSqFU9mHHnN+Gr9fvZkKKj2SrfpgVD/W7lzjR+2HKIW/s0o0ZooNNxKo3xfZpSMzSAZ+O3OB1FKUdpwVAAGGN4+pst1KoWxE29fKOTXmlVDw7gjgtj+HnbYZYmHXY6jlKO0YKhAPh+8yFW7TrKXQNiCQn0nU56pXVt90bUCwvmmW+36JAhymdpwVDkFxiei99Ck8gqXBkXXfIGPig4wMU9A5uzPiWd+RsOOB1HKUe4pWCIyBAR2SoiSSIysZj1F4jIahHJE5GRRdaNFZFt9mNsoeWdRWSD/Z4viy/2Hqsg89bt5beDJ7lvUHOvnOvCXS7v1IDYWlV5YcFWnZlP+aQy/3YQERfwGjAUaA2MEZHWRZrtBm4AZhbZNhx4DOgGdAUeE5HTPcXeAMYDsfZjSFmzqj/LzS9g8oJttKpbnWFt6zodp1Jz+Qn3DmxOcmoGn6/Z63QcpSqcO/6c7AokGWO2G2NygFnAiMINjDE7jTHrgYIi2w4GFhhj0owxR4EFwBARqQtUN8b8YqwTxu8Bl7ohqyrik4QUdqdlcv+g5vj56UFcSYa0rUPb+tV58fvfyMkr+nVWyru5o2DUB/YUep1iLyvLtvXt5+fynqqUsnLzeeXHbXRsWIN+LWs5HccjiAj3DWpBytFTfJywp+QNlPIi7igYxf1ZWtoTvGfattTvKSLjRSRBRBJSU1NL+bEK4MPlu9mfnsUDg1r45ACD56pv8yjiGtXk1R+3kZWrw58r3+GOgpECFL61pgGwr4zbptjPS3xPY8wUY0ycMSYuKiqq1KF9XWZOHm8sSqJnswh6xkQ6HcejiAj3D27BwePZfPCrDn+ufIc7CsZKIFZEmohIIHAVMK+U28YDg0Skpn2xexAQb4zZD5wQke723VHXA3PdkFXZ3l22i8Mnc7hvUHOno3ik7k0jOD8mkjcWJZORned0HKUqRJkLhjEmD5iA9ct/MzDbGJMoIpNEZDiAiHQRkRRgFPCWiCTa26YB/8IqOiuBSfYygNuAaUASkAx8U9asynIiK5e3fkqmb4soOjcKdzqOx7pnYHOOZOTw7i87nY6iVIXwd8ebGGPmA/OLLHu00POV/PEUU+F2bwNvF7M8AWjrjnzqj95ZupNjmbncM0CPLsqic6Oa9G0RxZSftnNd90ZUCw5wOpJS5Up7afmY9FO5TP15OwNa1aZ9dA2n43i8ewc251hmLjOW7nQ6ilLlTguGj5n+83aOZ+Vx70A9unCHdg1qMLB1bab+vJ30zFyn4yhVrrRg+JCjGTm8vXQnw86rQ+t61Z2O4zXuGdCcE1l5TFuy3ekoSpUrLRg+ZOrP28nIyeOu/np04U6t61Vn2Hl1mLF0J0czcpyOo1S50YLhI46czOadZTu5uF09WtSp5nQcr3NX/+Zk5OQx9Wc9ylDeSwuGj5jy03aycvO5q3+s01G8Uos61bi4XT3eWbaTIyeznY6jVLnQguEDUk9k894vuxjRoT4xtao6Hcdr3dU/lqzcfKboUYbyUlowfMBbi5PJzsvn7/1inI7i1WJqVWV4+3q8t2wXqSf0KEN5Hy0YXu7Q8Sze/3UXl3VsQNMoPboob3f2jyU7L583Fyc7HUUpt9OC4eVeX5RMXoHhzv56dFERmkZV5bKODfjg110cOp7ldByl3EoLhhfbn36KmSt2M7JTAxpFVHE6js+4s38MeQWG1xfpUYbyLlowvNjrC5MpKDBM0GsXFapRRBVGdmrAzBW72Z9+yuk4SrmNFgwvtffYKWat3M2VXaKJDg91Oo7PmdAvhoICw+sL9ShDeQ8tGF7q1R+TEIQJF+rRhROiw0O5sks0s1buZu8xPcpQ3kELhhfak5bJJwl7GNM1mno1QpyO47MmXBiDILz64zanoyjlFlowvNBLP2zD5Sfc1lePLpxUr0YIY7pG80lCCruOZDgdR6ky04LhZZJTT/LZ6hSu7d6IOmHBTsfxeXdcGIPLT3jpBz3KUJ7PLQVDRIaIyFYRSRKRicWsDxKRj+31y0Wksb38GhFZW+hRICId7HWL7Pc8va6WO7J6u5e+30aQv4vb+jZzOooCalUP5voejfhizV6SDp10Oo5SZVLmgiEiLuA1YCjQGhgjIq2LNBsHHDXGxACTgWcAjDEfGmM6GGM6ANcBO40xawttd83p9caYQ2XN6u22HDjOl+v3cWOvxkRWDXI6jrLd2qcZwQEuXvz+N6ejKFUm7jjC6AokGWO2G2NygFnAiCJtRgDv2s/nAP1FRIq0GQN85IY8Pmvygt+oGujP+AuaOh1FFRJRNYgbezXmq/X72bz/uNNxlDpn7igY9YE9hV6n2MuKbWOMyQPSgYgibUbz54Ixwz4d9c9iCowqZENKOvGJBxnXuwk1QgOdjqOKGN+7GdWC/XlhgR5lKM/ljoJR3C9yczZtRKQbkGmM2Vho/TXGmPOA3vbjumI/XGS8iCSISEJqaurZJfciz3+3lRqhAdx0fhOno6hihIUGcHPvpizYdJC1e445HUepc+KOgpECRBd63QDYd6Y2IuIPhAFphdZfRZGjC2PMXvvnCWAm1qmvPzHGTDHGxBlj4qKiosqwG55r+fYjLP4tldv7NqN6cIDTcdQZ3HR+EyKqBPJ8/Fanoyh1TtxRMFYCsSLSREQCsX75zyvSZh4w1n4+EvjRGGMARMQPGIV17QN7mb+IRNrPA4CLgY2oPzHG8Px3W6ldPYjrezR2Oo76C1WD/Ln9whiWJB1mWdJhp+ModdbKXDDsaxITgHhgMzDbGJMoIpNEZLjdbDoQISJJwL1A4VtvLwBSjDGFpykLAuJFZD2wFtgLTC1rVm+06LdUVu48yt/7xRIc4HI6jirBNd0aUjcsmOe+24r9N5NSHkO86UsbFxdnEhISnI5RYQoKDJe8uoQTWXl8f28fAv21H6YnmLViNxM/28C06+MY0Lq203GUQkRWGWPiSmqnv2E82PyN+0ncd5x7BsZqsfAgV3RuQJPIKjwXv5X8Au/5g015P/0t46Fy8wt4Ln4rLetUY3j7oncxq8oswOXHfYOas/XgCb5Ys9fpOEqVmhYMDzVrxW52HcnkH0Na4PLTLiqeZljbupxXP4wXFvxGVm6+03GUKhUtGB4oIzuPl35IomuTcC5soUNseSI/P2Hi0JbsPXaKD37d5XQcpUpFC4YHmr5kB4dPZjNxaEu0A7zn6hUTSe/YSF5dmMTxrFyn4yhVIi0YHubIyWzeWpzM4Da16dSwptNxVBk9OKQlxzJzeWuxTuWqKj8tGB7m5R+2kZVXwAODWzodRblB2/phDG9fj+lLdrA/XadyVZWbFgwPkpx6kg+X7+aqLtHE1KrqdBzlJg8MbkFBAfz3Ox2YUFVuWjA8yNPfbCE4wMU9A5s7HUW5UXR4KDf2asynq1NI3JfudBylzkgLhof4dfsRFmw6yG19m+nkSF7o9gtjqBESwFNfb9YhQ1SlpQXDAxQUGJ76ejP1woIZp8OXe6WwkADu6h/LsuQjLNyqk0uqykkLhgeYu24vG/am88CQFjrAoBe7ulsjmkRW4d/zt5CbX+B0HKX+RAtGJZeRncfT32yhXYMwRugQIF4t0N+Ph4a2JOnQSe3MpyolLRiV3BuLkjl4PJvHLmmDnw4B4vUGtq7N+TGRTF7wG2kZOU7HUeoPtGBUYruPZDLl5+1c1rE+nRtpJz1fICI8dklrMnLy+e93OjOfqly0YFRiT83fhEuEB4doJz1fElu7Gtd1b8RHK3azad9xp+Mo9Tu3FAwRGSIiW0UkSUQmFrM+SEQ+ttcvF5HG9vLGInJKRNbajzcLbdNZRDbY27wsPjZo0tKkw8QnHuSOC5tRJyzY6Tiqgt0zoDlhIQE88WWi3marKo0yFwwRcQGvAUOB1sAYEWldpNk44KgxJgaYDDxTaF2yMaaD/bi10PI3gPFArP0YUtasniInr4DH5iUSHR7C33o3dTqOckBYaAD3DWrB8h1pzFu3z+k4SgHuOcLoCiQZY7YbY3KAWcCIIm1GAO/az+cA/f/qiEFE6gLVjTG/GOvPq/eAS92Q1SNMW7KdpEMneWJ4G72N1oeN6dqQdg3CePLrzTqaraoU3FEw6gN7Cr1OsZcV28YYkwekAxH2uiYiskZEFotI70LtU0p4T6+UcjSTl3/YxuA2tenXUud79mUuP+GpS8/j8MlsXtBxplQl4O+G9yjuSKHoSdcztdkPNDTGHBGRzsAXItKmlO9pvbHIeKxTVzRs2LDUoSurJ77chCA8ekmbc3uD9BTYtQxSVsKxPXB8L5w8BH4uCAiBgFCIaAZ121uP6O4QGOrenVBuc16DMK7r3oj3ftnJyM4NaFs/zOlIyoe5o2CkANGFXjcAip50Pd0mRUT8gTAgzT7dlA1gjFklIslAc7t9gxLeE3u7KcAUgLi4OI++OvjD5oMs2HSQiUNbUr9GSOk3PLQZ1n0EiV/AMbvDV2A1qNkYqte1CgMGck9B9gnYuxoSP7fbVYVWl0C7K6FJH6uwqErlvkEtmL9hP498sZHPb+up/XGUY9xRMFYCsSLSBNgLXAVcXaTNPGAs8AswEvjRGGNEJAqrcOSLSFOsi9vbjTFpInJCRLoDy4HrgVfckLXSOpmdxz+/2Ehsrarc1KsU40Xl58L62bDiLdi/DsQFMf2h+23QqCfUbvvXv/wz02Dfatg0FxLnWgUnsgX0ewRaDQffuimtUgsLCeD/LmrN3R+v5f1fdzG2Z2OnIykfVeaCYYzJE5EJQDzgAt42xiSKyCQgwRgzD5gOvC8iSUAaVlEBuACYJCJ5QD5wqzEmzV53G/AOEAJ8Yz+81nPfbmH/8Szm3NqTQP+/uLSUmwWr34Vlr0D6HqjVBoY8A22vgKpRpf/A0HCIGWA9hj4HW76Cxc/C7OuhXkcY/G+r8KhKYUSHeny2Zi/PfruFAa1rn90RqFJuIt50j3dcXJxJSEhwOsZZW7UrjZFv/sLYHo15fPgZrl0YAxs/he+fgPTd0LAHnH8vxA5039FAQT6smwWL/mNdC+k5Afr9E/x1OPXKYE9aJoNf/ImuTcKZcUMXnc9duY2IrDLGxJXUTnt6Oyw7L58HP91AvbAQHhjcovhG+9bC24Ph03EQHAbXz4WbvoXmg9x76sjPBR2vgTuWQ+cbrKOYKRfCwU3u+wx1zqLDQ7l/UAsWbU1l7lrtm6EqnhYMh732YxJJh07y1GVtqRJU5AxhTiYseBSm9oO0HTD8FbhlMTTtW76hAqvAJS/C1bMhIxWmDYDfvivfz1SlMrZnYzpE1+CJLxM5fDLb6TjKx2jBcNC6Pcd4bVEyl3eqT98Wtf64cudSeKMnLH3J+qt/wkrodH3F3sXUfDDc+jNExsBHo2HltIr7bFUsl5/w7Mh2ZGTn88jnG3TYEFWhtGA4JCs3n3tnr6VWtSAeK9znIi8Hvn8c3rkIMDD2S+vIIqSGM0Gr1YEb5kPMQPj6PljwmHU9RTmmee1q3DeoOfGJB/l8zV6n4ygfogXDIc9+u5Xk1AyeG9mesJAAa+HhbTB9ACyZDJ2ug1uXQpMLnA0KEFQVrpoJnW+EpS/CD09o0XDY33o3pUvjmjw2N5F9x045HUf5CC0YDliWfJi3l+5gbI9GnB8baS1cPxve6mP1zh79gXVUEVTV2aCFufzh4slW0VgyGRY/U/I2qty4/IT/jupAvjE8MGcdBQVawFX504JRwY5l5nD/7HU0iazCxKGtrN7X8+6Ez262emTfttTqeV0ZicBFL0CHa61bb5dMdjqRT2sYEcr/XdSapUlHeHvpDqfjKB/gjp7eqpSMMfxjznpST2bz6W09CTm5Gz6+Dg5ugPPvgQv/z/pLvjLz84PhL0N+tnWtJSwazhvpdCqfNaZrNIu2HuKZb7fQtUk47Ro4dK1L+QQ9wqhA7/+6i+82HeTBIS1pd2olTOlrdcK7ejYMeLzyF4vT/Fww4nVo1Au+uB32rHQ6kc8Sse6aiqoaxN8/WsMJHQZdlSMtGBUkcV86T361mf4tIhlXMAc+HGX9dT5+kXX7qqfxD4Qr34fq9WDWGDi6y+lEPqtGaCAvj+lIytFTPPL5Rr3VVpUbLRgV4HhWLn+fuYb6oXm8HvgisvApa3TYcd9BuAfPqFclwjo6ysuBj8ZATobTiXxWXONw7hkQy7x1+/hoxZ6SN1DqHGjBKGcFBYZ7P16HpG3n69AnCEqKhyFPw2Vvecc8FFHNYdQMOLQJ5v/D6TQ+7ba+MVzQPIrH5yWyZvdRp+MoL6QFo5y9tjCJnK3fMT/0UUJzjsB1n1tDkHvTwHEx/eGCB2DtB7B2ptNpfJbLT3j5qg7UDgvitg9Wk3pChw5R7qUFoxwt3HKQEwtf4J3A5wiMaAjjF0LTPk7HKh99J0Lj3vDVvdaETsoRNUIDefPazhw7lcMdM1eTm1/gdCTlRbRglJPkfYfImHUTD/vPpKDlJci4BdYMeN7KzwVXTLM6G35ygzVwonJEm3ph/Ofy81ixI43H5yXqRXDlNlowykHavmTypg5mGEs51uNB/Ee/a40A6+2q1YHLp0DqFvhhktNpfNplHRtwS5+mfLh8N9OXaKc+5R5aMNwsO/lnXFMvpF7BfnYOnE6NwQ971/WKkjTrB13Hw/I3YPtip9P4tAcHt2Ro2zo8NX8z8YkHnI6jvIBbCoaIDBGRrSKSJCITi1kfJCIf2+uXi0hje/lAEVklIhvsn/0KbbPIfs+19qNW0fetVIyh4Jc38H9/BIfzQ1k7+FOa9rrC6VTOGPAERMTA3DsgK93pND7Lz0+YPLoD7RrU4K5Za1i355jTkZSHK3PBEBEX8BowFGgNjBGR1kWajQOOGmNigMnA6ZHrDgOXGGPOA8YC7xfZ7hpjTAf7caisWctNTibm8/H4xU/kx/wO/NR3Fr179nI6lXMCQ+HSN+H4Xvj2IafT+LTgABfTro8jsmoQN8xYwbaDJ5yOpDyYO44wugJJxpjtxpgcYBYwokibEcC79vM5QH8REWPMGmPM6bkmE4FgEfGsCaTTdsDbg2D9JzyfO4pVPV7lxn4dnE7lvOgu1vhYaz+EbQucTuPToqoF8cG4bvi7/Lh2+nL2pOkNCercuKNg1AcKdy1NsZcV28YYkwekAxFF2lwBrDHGFL55fIZ9OuqfcoYZ70VkvIgkiEhCampqWfbj7G2ZD2/1IfvwTm7MuZ+0uLt4cGiris1QmfV5ECJbWLfaZp90Oo1PaxxZhffHdSUrt4Brpy/n0PEspyMpD+SOglHcL/Ki9/H9ZRsRaYN1muqWQuuvsU9V9bYf1xX34caYKcaYOGNMXFRU1FkFP2f5edZIrbPGkBraPoimAAAZ8UlEQVRYl/4Z/6LaeRfxrxFtOUNd803+QXDJS9YAiwv/7XQan9eyTnVm3NiF1BPZjJn6Kwe1aKiz5I6CkQJEF3rdANh3pjYi4g+EAWn26wbA58D1xpjk0xsYY/baP08AM7FOfTkvPQXevRiWTCax7uWcnzqR9ue154Ur2+Py02LxJ416WJMuLX8D9q52Oo3P69SwJu/c2JUD6VmMfusXna1PnRV3FIyVQKyINBGRQOAqYF6RNvOwLmoDjAR+NMYYEakBfA08ZIxZerqxiPiLSKT9PAC4GNjohqxls2U+vHk+5sAG5jd/kot2jGRYxya8NLoDAS69Q/mMBj4BVWrBl3dCvg6/7bSuTcJ5b1w3jpzMYfSUX/Sahiq1Mv+Ws69JTADigc3AbGNMoohMEpHhdrPpQISIJAH3AqdvvZ0AxAD/LHL7bBAQLyLrgbXAXmBqWbOes5xM+Po+mDUGU6MhLzabzu3rmzKqcwOeH9Uefy0Wfy04DIY9Cwc2wArn/jeq/+ncqCYf/K0b6Zm5XPHGMhL36e3PqmTiTcMGxMXFmYSEBPe+6b418OnNcGQbuV1v587US/hm81FuuaApDw5piZ+ehiodY+DDkbBnBUxYafUKV47beuAEN8xYwfFTubxxbWcuaF5B1wFVpSIiq4wxcSW10z+NzyQvBxY9DdMGQE4Gx0Z9ypU7L+bbLUd5YngbHhrWSovF2RCBoc9CXhYseNTpNMrWok41Pr+9F9Hhodz0zkpmr9S5NNSZacEozv51MLUfLPoPtLmM9cO/YchcYfP+47xxTWfG9mzsdELPFNEMet4J6z+GnUtLbq8qRJ2wYD65tQc9mkXwj0/X8+jcjeTk6Si36s+0YBSWk2ndLju1H2QcgqtmMiv6UUa+s5kAf+Gz23oxpK2eSimT3vdZU9POf8C6PVlVCtWCA5hxQxdu7t2E937ZxTXTfuXQCb3tVv2RFozTti2A17vDksnQbjTHxy3hnnX1mfjZBro1DefLCefTul51p1N6vsBQGPIfOJQICW87nUYV4u/y45GLWvPymI5s2JvOsJeWsGhr5R2RR1U8LRgAP//XuiDrHww3fE1ChycZNmUj89bt4+4BsbxzY1dqhAY6ndJ7tLwYmvSBRf+GzDSn06gihrevx9w7zieiSiA3zFjJE18mkpWb73QsVQlowQBofSn0+yenxi3mP5sjuPKtXxCB2bf04O4BzbVDnruJWEcZWenWjQWq0mlRpxpzJ/Tihp6NmbF0J5e8soRVu3SecF+nBQMgohk/1x3L4Fd+5a3F2xnVOZr5d/amc6OaTifzXrXbWD3AV06DQ1ucTqOKERzg4vHhbXjnxi5kZOcx8s1lPDp3IyeytPOlr9KCAby5OJnrpq/A30/46ObuPDOyHdWCA5yO5f0ufMSa0jX+YaufhqqU+raoxXf39uGGno15/9dd9P/vYj5J2ENBgf4/8zVaMICBrWtzZ78Y5t/Vmx7Nig6iq8pNlQjoMxGSf4Bt3zmdRv2FqkH+PHZJGz6/vRf1aoTwwJz1XPzKEpYlHXY6mqpA2tNbOSs/17o7TVxw2zJw+TudSJXAGMOX6/fzzDdb2HvsFN2ahHPXgFh6NI3Q0Zo9lPb0Vp7BFQADJ8HhrbD6HafTqFIQEYa3r8cP9/Xh0Ytbs+NwBldPXc6oN3/h240HyNdTVV5LjzCU84yBdy6C1K1w5xoI1v4uniQrN5+PV+5hyk/b2XvsFA1qhjC2R2Mu71SfiKqeNYGmryrtEYYWDFU57F0NUy+0pnUd8LjTadQ5yMsvYMGmg8xYupMVO9Pw9xP6tazFqLhoLmgeSZC/y+mI6gy0YCjP89l4SPwC/p4ANRo6nUaVwdYDJ5izag+fr9nL4ZM5VAv2Z2Cr2gxpW4fesVGEBGrxqEy0YCjPk54Cr3SGVsPhCp03wxvk5hfw87ZUvtlwgO82HST9VC6B/n50axJOn+ZRdG8aQau61bVzrMO0YCjP9P0TsOQFuHkh1O/kdBrlRrn5BSzfnsairYdY9FsqSYdOAlAt2J8ujcPpGF2DdtE1aFc/jJpVdCieilShBUNEhgAvAS5gmjHm6SLrg4D3gM7AEWC0MWanve4hYByQD9xpjIkvzXsWRwuGF8hKh5c7Qq3WMPZLaxgR5ZX2HTvFih1pLN+RxoodR0hOzfh9Xb2wYFrUqUaLOtWJrVWVJlFVaBpZRcd0KycVVjBExAX8BgwEUrDm+B5jjNlUqM3tQDtjzK0ichVwmTFmtIi0Bj4CugL1gO+B5vZmf/mexdGC4SWWT4FvHoAxH0OLIU6nURXkeFYuG1PSWZeSzpYDx9l64ARJh06SV+g23bCQAOrXCKFBzRDq1QihdvVg6oQFUbtaMBFVg4ioGkjN0EA9xXWWSlsw3NFLqiuQZIzZbn/wLGAEUPiX+wjgcfv5HOBVsXr4jABmGWOygR32nN9d7XYlvafyVnE3wvI3rZn5YgZoZz4fUT04gJ4xkfSMifx9WU5eAXuOZrIjNYMdhzPYnZZJytFMth/OYFnyEU5m/3lOFRHrvWqEBhAWEkC1YH+qBVk/qwT5UyXIRWigPyEBLkICXYQEuAgO8CPI30WQvx8B/n4EuvwI9PcjwOVHgEvwd/kR4Ce4/AR/Pz/8/Pj9p0us5b7QadEd/xLrA4XndUwBup2pjTEmT0TSgQh7+a9Ftq1vPy/pPZW3cgVYt9bOvg7WvG8VEOWTAv39aBZVlWZRVYtdfzI7j4PHszh0PJsjGdkcOZnDkYwc0jNzOHYql2OZuZzIyuXQ8WxOZOWRkZNHRnYef+xbaAgmh2qcoqqcIpQsqpBFqGQRQg7B5BAiOQSRQxC5BJFLoOQSQD6BWD8DyCNA8vEnH38pIIB8/CUfFwY/CnBRgEusn2Iv88P8/pBCjz++tvIJ/L4MwXpd5OzQvk730HX4reXy/+E0dxSM4spq0fNcZ2pzpuXF9UAv9tyZiIwHxgM0bKi3YnqNVpdAdHdY+G84b5Q1SKFSRVQN8qdqcQUlJxOO74Xj++Dkwf89MtMwGYcxmUcwmUch6yh+2ceRgrOb/TFf/CmQAPL97J/iby/zp0Bc5OOiQOwHfhSIC0MABeKHMYKR0yUDq0SI2KVECpWHP5YKhD+UEX5fb6kSXu+c/zuWljsKRgoQXeh1A2DfGdqkiIg/EAaklbBtSe8JgDFmCjAFrGsY57YLqtIRgUFPwvQBsOwVuPAhpxOpyiYzDQ7/Zj2OJMPRndbj2G44VczEXK4gqBKJhEYgoeFQIxpCakJwmPUIqgaB1aw/TgKrQEAVa4bIgBDwD4GAYGuSNVcQLj8/XICvjWntjoKxEogVkSbAXuAq4OoibeYBY4FfgJHAj8YYIyLzgJki8gLWRe9YYAXWkUdJ76m8XXQXa3KrZS9bp6Wq6XzqPskYSNsOe1fBgfVwYCMcTISMQtPH+gVYnT1rNrZux65eH8IaQPV6ULUOVKsNQdX1rrsyKnPBsK9JTADisW6BfdsYkygik4AEY8w8YDrwvn1ROw2rAGC3m411MTsPuMMYkw9Q3HuWNavyQAMegy1fW6emhr/sdBpVEfKyraFidi2FXcusQpF1zFrnCoJaLSF2INRqBZEtIDIGajQCP+09Xt60456q/L6ZCCvesoY/r9XK6TTK3YyxTisl/WDNjbJzKeSdstZFtYLorlC/s/WIaql3zZWDirytVqny1ecfsHamdZvtNZ84nUa5Q0EBpKyELV9ZR5BpydbyyObQeSw0uQAa9oDQcGdzqj/QgqEqv9BwuOB+WPBPSF4IzS50OpE6VwcTYcMnsOFTSN9tXXtocgH0uB1iB+mgk5WcFgzlGbrdAiunwXf/B7f8pOerPUn2CdgwB1a/C/vWWLMrNusH/R6BFkOtO5SUR9CCoTyDf5DVmW/OjdbpqU7XOZ1IleRwktVjf+1MyM2wxgcb8gycNxKqRJa8vap0tGAoz9HmMvj1dfjxSeu5duarnHYuhaUvwbZ4cAVC25HQZZx10Vpva/VoWjCU5xCBQU/B24OsvhkXPux0InWaMdYdTj/9F3Yvg9BI6PsQxN0EVWs5nU65iRYM5VkadrOOLpa+DJ2utzpnKWftXAI/TII9y60Oc0OfhY7XWb2klVcpbswmpSq3gZMAAwseczqJbzuwET64At65yBqO4+LJcOda6wYFLRZeSY8wlOep0RB6/h1+eg66jreOOlTFyThsXUda/a413MbASdb/h4AQp5OpcqZHGMoz9bobqtWFbx+0OoGp8leQD8vfgpc7wer3rCJx5xrodZcWCx+hBUN5pqCq1m22+9bA+llOp/F+e1fD1H7wzT+gQWe4/RcY+oz2xPYxWjCU5zrvSqgfZ13LyEp3Oo13ysmAbx+Caf3hxH4YOQOu/QyiWjidTDlAC4byXH5+MOw5yEiFRU87ncb77FwKb/Sy+r50vhEmrIS2l2tfCh+mBUN5tvqdoPMN1rn1gzoCvlvknoJvHoR3hgEGxn4FF7+gQ3goLRjKC/R/1Ppl9vX9Vgcyde72r4O3+lhDenS9xRpSvklvp1OpSkILhvJ8oeHWREu7l8H62U6n8UwFBdZwHlP7W9eDrv0Mhj1rTVWqlE0LhvIOHa+3xir67hE4ddTpNJ7lZCrMHGXNN9J8sHVUEdPf6VSqEipTwRCRcBFZICLb7J81z9BurN1mm4iMtZeFisjXIrJFRBJF5OlC7W8QkVQRWWs//laWnMoH+PnBxS9CZhp890+n03iOHT/Dm+dbPy/6L4z+AKpEOJ1KVVJlPcKYCPxgjIkFfrBf/4GIhAOPAd2ArsBjhQrL88aYlkBHoJeIDC206cfGmA72Y1oZcypfULcd9JwAa963fgGqMysogJ+eh/eGW31abv4BuvxN74BSf6msBWME8K79/F3g0mLaDAYWGGPSjDFHgQXAEGNMpjFmIYAxJgdYDehIcqps+kyEmo3hq7shN8vpNJVTZhp8dBX8+C9rIMfxi6HOeU6nUh6grAWjtjFmP4D9s7hxjOsDewq9TrGX/U5EagCXYB2lnHaFiKwXkTkiEn2mACIyXkQSRCQhNTX1XPdDeYvAUGsQvCNJ1lhT6o/2r4cpfSD5Rxj6HFwxXecVUaVWYsEQke9FZGMxjxGl/IzijnF/v/dRRPyBj4CXjTHb7cVfAo2NMe2A7/nfUcyf38iYKcaYOGNMXFRUVCkjKa/WrB+0HwNLX7SGtFCW9bNh+iDIz4Mbv4Fu4/UUlDorJRYMY8wAY0zbYh5zgYMiUhfA/nmomLdIAQofITQA9hV6PQXYZox5sdBnHjHGZNsvpwKdz263lM8b8h+oUgs+v9XqiObL8nPhm4nw2c1WR8dbFkN0F6dTKQ9U1lNS84Cx9vOxwNxi2sQDg0Skpn2xe5C9DBF5EggD7i68wekiZBsObC5jTuVrQmrCiFfh8FZrKG5flXEE3r8Mlr8B3W6D6+fqDHjqnJW1YDwNDBSRbcBA+zUiEici0wCMMWnAv4CV9mOSMSZNRBoAjwCtgdVFbp+9077Vdh1wJ3BDGXMqXxTTH+LGwS+vWbPC+Zr962FKX9izAi59E4Y+Da4Ap1MpDybGi4ZSiIuLMwkJCU7HUJVJToY1gF5BPtz6k3Xk4Qs2fgpf3GH1gh/9gXUqSqkzEJFVxpi4ktppT2/l3QKrwBXT4MQ+mDvB+8eaKsiH75+AOTdB3fYwfpEWC+U2WjCU92sQZ00juuUr+PUNp9OUn6x0+GgMLHnBGsF37Jd6vUK5lc7prXxD99ut+R0WPArR3axZ47xJ6laYdTUc3WkN8dFFR9NR7qdHGMo3iMClr1nzgH8yFjIOO53IfbbM/98os2O/1GKhyo0WDOU7QmrCle9aM/TNutrzhw4pKICF/4ZZYyCimXW9olFPp1MpL6YFQ/mW+p3g0jdgz3L48k7PvQiemQYzr4TFz0D7q+GmbyFMh2JT5UuvYSjf0/ZyOJIMC5+EiFjo84DTic7O/nUw+3pI3wsXvQBxN+kQH6pCaMFQvumC++HINqtoVImEuBudTlQyY2DlNIh/GKpEWeNB6RAfqgJpwVC+SQSGv2LNzvfVPeAfBB2udjrVmWUdhy/vgsTPIGYgXPaWTnSkKpwWDOW7/IPgyvfho9Ew9w5wBcJ5I51O9We7foHPx0N6CvR/FHrdY80wqFQF02+d8m0BwXDVR9CwJ3w2HhJmOJ3of/JzrYET3xkG4gc3xUPv+7RYKMfoN0+pwFC4+mNrHo2v7obvH7duWXXS3tUw5UJrEqj2V8OtSyC6q7OZlM/TU1JKgTXr3JhZMP9+WDIZju6yhkcPrFKxOXIyrL4Vv75uzecx+kNodXHFZlDqDLRgKHWay9+a3jW8CSx4DA6sh8umVMwwIgUFsGG2NXDgiX3WWFADnoCQGuX/2UqVkp6SUqowEeh1F4ydZ/UEnz4QFj0NeTnl83nGwPbFMK0/fH6LNVjgjd/CJS9psVCVjhYMpYrT5AK4bSm0vQIW/QdejYN1s6zhw92hoAA2f2UViveGw4n91iRHNy+ERj3c8xlKuVmZCoaIhIvIAhHZZv8sdnYaERlrt9kmImMLLV8kIlvt2fbWikgte3mQiHwsIkkislxEGpclp1LnJKQGXDEVrvkUgqtbRwBvnm/dSXXq6Lm955Fk+PEpeLkDfHwNZB6xToPduRY6jNE7oFSlVqYZ90TkWSDNGPO0iEwEahpjHizSJhxIAOIAA6wCOhtjjorIIuB+Y0xCkW1uB9oZY24VkauAy4wxo0vKozPuqXJTUACbPofFz0HqZqvPRuwgaNrXmqiodps/XyDPy7GuR+xfB7uWWcOrH9wACDTtAx2vg9aXWtdOlHJQaWfcK+s3dQTQ137+LrAIeLBIm8HAAntub0RkATAE+KiE933cfj4HeFVExHjTfLLKs/j5Waen2lxuFYD1s61e11u+staLHwRVh4BQq29H9glrVNzT/EOs22IHToLzRkH1es7sh1JlUNaCUdsYsx/AGLP/9CmlIuoDewq9TrGXnTZDRPKBT4En7aLw+zbGmDwRSQciAC+axEB5JBGo18F6DH4KjttHEAfWW6eXcjMh9xQEVoXq9aF6XYhqCXU7gH+g0+mVKpMSC4aIfA/UKWbVI6X8jOKG0Tx9pHCNMWaviFTDKhjXAe+VsE3RfOOB8QANGzYsZSSl3EAEwupbj5bDnE6jVLkrsWAYYwacaZ2IHBSRuvbRRV3gUDHNUvjfaSuABlinrjDG7LV/nhCRmUBXrIKRAkQDKSLiD4QBaWfINwWYAtY1jJL2Ryml1Lkp6y0Z84DTdz2NBeYW0yYeGCQiNe27qAYB8SLiLyKRACISAFwMbCzmfUcCP+r1C6WUclZZr2E8DcwWkXHAbmAUgIjEAbcaY/5mjEkTkX8BK+1tJtnLqmAVjgDABXwPTLXbTAfeF5EkrCOLq8qYUymlVBmV6bbaykZvq1VKqbNX2ttqtZeQUkqpUtGCoZRSqlS0YCillCoVLRhKKaVKxasueotIKrDrHDePxHt6kuu+VD7esh+g+1JZlWVfGhljokpq5FUFoyxEJKE0dwl4At2Xysdb9gN0XyqritgXPSWllFKqVLRgKKWUKhUtGP8zxekAbqT7Uvl4y36A7ktlVe77otcwlFJKlYoeYSillCoVLRiAiAyx5xZPsqea9Rgi8raIHBKRjYWWlWqu9cpERKJFZKGIbBaRRBG5y17uifsSLCIrRGSdvS9P2Mub2HPUb7PnrPeIGZVExCUia0TkK/u1p+7HThHZICJrRSTBXuZx3y8AEakhInNEZIv9b6ZHReyLzxcMEXEBrwFDgdbAGBFp7Wyqs/IO1pS3hU0EfjDGxAI/2K8ruzzgPmNMK6A7cIf9/8ET9yUb6GeMaQ90AIaISHfgGWCyvS9HgXEOZjwbdwGbC7321P0AuNAY06HQ7aee+P0CeAn41hjTEmiP9f+n/PfFGOPTD6AHEF/o9UPAQ07nOst9aAxsLPR6K1DXfl4X2Op0xnPYp7nAQE/fFyAUWA10w+pU5W8v/8P3rrI+sCY8+wHoB3yFNRumx+2HnXUnEFlkmcd9v4DqwA7sa9AVuS8+f4RByXOOe6I/zLUOFDfXeqUlIo2BjsByPHRf7NM4a7FmoVwAJAPHjDF5dhNP+Z69CPwDKLBfR+CZ+wHWNM/ficgqe2pn8MzvV1MgFZhhnyqcZs8vVO77ogXjLOYPV+VPRKpize9+tzHmuNN5zpUxJt8Y0wHrL/SuQKvimlVsqrMjIhcDh4wxqwovLqZppd6PQnoZYzphnX6+Q0QucDrQOfIHOgFvGGM6AhlU0Kk0LRj/mz/8tAbAPoeyuMtBe451/mKu9UrHnn3xU+BDY8xn9mKP3JfTjDHHsOaw7w7UsOeoB8/4nvUChovITmAW1mmpF/G8/QDAGLPP/nkI+ByrkHvi9ysFSDHGLLdfz8EqIOW+L1owrKljY+07PwKxpoOd53CmsirNXOuViogI1tS8m40xLxRa5Yn7EiUiNeznIcAArIuSC7HmqAcP2BdjzEPGmAbGmMZY/y5+NMZcg4ftB4CIVBGRaqefA4OAjXjg98sYcwDYIyIt7EX9gU1UwL5oxz1ARIZh/eXkAt42xjzlcKRSE5GPgL5YI1UeBB4DvgBmAw2x51o3xqQ5lbE0ROR84GdgA/87X/4w1nUMT9uXdsC7WN8nP2C2MWaSiDTF+ks9HFgDXGuMyXYuaemJSF/gfmPMxZ64H3bmz+2X/sBMY8xTIhKBh32/AESkAzANCAS2Azdif9cox33RgqGUUqpU9JSUUkqpUtGCoZRSqlS0YCillCoVLRhKKaVKRQuGUkqpUtGCoZRSqlS0YCillCoVLRhKKaVK5f8B+U2lbE8Luv4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(tlist,np.real(np.array([H1_coeff(t, Hargs) for t in tlist])))\n",
    "plt.plot(tlist,np.imag(np.array([H1_coeff(t, Hargs) for t in tlist])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 307,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T16:06:39.660902Z",
     "start_time": "2020-05-05T16:06:39.657232Z"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.integrate import quad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 345,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-05-05T16:15:34.463965Z",
     "start_time": "2020-05-05T16:15:34.454959Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0986049370977187"
      ]
     },
     "execution_count": 345,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "quad(H1_coeff,0,50, args=Hargs)[0]/np.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}