{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Damped Rabi Flopping\n",
    "\n",
    "This example covers damped Rabi flopping as calculated with the optical Bloch equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pylcp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the problem\n",
    "\n",
    "As always, we must define the laser beams, magnetic field, and Hamiltonian.  Here, we will make a two-state system that is magnetic field insensitive and connected only by $\\pi$ light.  In this particular case, we show how we can define the rotating frame such that the excited state of the Hamiltonian can rotate and the laser beams can rotate, or some combination of the two.  The total detuning is the sum of `ham_det` and `laser_det`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[((<g|H_0|g> 1x1), (<g|mu_q|g> 1x1)) (<g|d_q|e> 1x1)]\n",
      " [(<e|d_q|g> 1x1) ((<e|H_0|e> 1x1), (<e|mu_q|e> 1x1))]]\n"
     ]
    }
   ],
   "source": [
    "ham_det = -2.\n",
    "laser_det = -2.\n",
    "\n",
    "laserBeams = pylcp.laserBeams(\n",
    "    [{'kvec':np.array([1., 0., 0.]), 'pol':np.array([0., 0., 1.]),\n",
    "      'pol_coord':'cartesian', 'delta':laser_det, 's':20.}]\n",
    ")\n",
    "\n",
    "magField = lambda R: np.zeros(R.shape)\n",
    "\n",
    "# Now define the extremely simple Hamiltonian:\n",
    "Hg = np.array([[0.]])\n",
    "mugq = np.array([[[0.]], [[0.]], [[0.]]])\n",
    "He = np.array([[-ham_det]])\n",
    "mueq = np.array([[[0.]], [[0.]], [[0.]]])\n",
    "dijq = np.array([[[0.]], [[1.]], [[0.]]])\n",
    "\n",
    "gamma = 1\n",
    "\n",
    "hamiltonian = pylcp.hamiltonian(Hg, He, mugq, mugq, dijq, gamma=gamma)\n",
    "hamiltonian.print_structure()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create the governing equation\n",
    "\n",
    "In this example, we create both the rate equations and optical Bloch equations to compare.  We also print out the decay rates given the chosen $\\Gamma$ to ensure the decay matrix evolution is being constructed properly for a two-level system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0. +0.j  0. +0.j  0. +0.j  1. +0.j]\n",
      " [ 0. +0.j -0.5+0.j  0. +0.j  0. +0.j]\n",
      " [ 0. +0.j  0. +0.j -0.5+0.j  0. +0.j]\n",
      " [ 0. +0.j  0. +0.j  0. +0.j -1. +0.j]]\n"
     ]
    }
   ],
   "source": [
    "# First the OBE:\n",
    "obe = pylcp.obe(laserBeams, magField, hamiltonian,\n",
    "                     transform_into_re_im=False)\n",
    "rateeq = pylcp.rateeq(laserBeams, magField, hamiltonian)\n",
    "print(obe.ev_mat['decay'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evolve the state\n",
    "\n",
    "We are going to evolve for $T=4\\pi/\\Gamma$, to see multiple oscillations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_eval = np.linspace(0, 4*np.pi/gamma, 501)\n",
    "rho0 = np.zeros((hamiltonian.n**2,), dtype='complex128')\n",
    "rho0[0] = 1.\n",
    "obe.set_initial_rho(rho0)\n",
    "obe.evolve_density([t_eval[0], t_eval[-1]], t_eval=t_eval)\n",
    "\n",
    "N0 = np.zeros((rateeq.hamiltonian.n,))\n",
    "N0[0] = 1\n",
    "rateeq.set_initial_pop(N0)\n",
    "rateeq.evolve_populations([t_eval[0], t_eval[-1]], t_eval=t_eval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot it all up:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAFECAYAAAAHuRaaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd5xU1f3/8ddnZntf2lIXkCoIAhYQDCp20agxRaMxmkQTU4xpxm9M8ftNfqZ8041J9KvRmKgx9tiDFaMiRZHeO8tStvcyc35/nFlYtsCWmWUW38/HY7jsvefee6Z+7jn3FHPOISIiIvErcKQzICIiIoemYC0iIhLnFKxFRETinIK1iIhInFOwFhERiXMK1iIiInFOwVpERCTOKViLiIjEOQVrERGROKdgLSIiEucUrEVEROKcgrWIiEicU7AWERGJcwrWIiIica7XBWszm2hmT5uZM7PPHibtDDN7w8xKzexNMzutp/IpIiISLb0qWJvZFGAJ0NCBtNOAV4H5wEnAYuBFM5sa00yKiIhEmTnnjnQeOszMUoFhQB2wBbjGOffXdtI+Bgx1zs2I/G3ASmChc+6aHsmwiIhIFPSqkrVzrsY5t+5w6cwsEbgEeLLZvg54CZgbuxyKiIhEX8KRzkCMDAGCwKYW67cB/cwsxTlXe6gDmNnKdjaNBWqA7d3OpYiIHA2GAdXOuYGxOsHRGqybXrDSFusrI8scoLCLxw4QTMwcMCR/Qr+M5C4eQkREjhYbN26krq4upuc4WoP1vsgyq8X6bCDcbHu7nHMT21pvZisTcwdNSL3idzz81VlMHprTvZyKiEivNnHiRFatWhXT2tZedc+6E3YCDhjRYv1wYJdzrrE7Bx+Q6UvUn7l3IYVlh6xNFxER6bajMlg752qAecDFTevMLACcDzze3eNnpyYyNi+DspoGLr/7HWrqQ909pIiISLuOmmBtZueYWZGZnRdZ9XNgppn9wMxGA78D+gF3RON8n505AoAzxg8gNSkYjUOKiIi06agJ1vh70XVEBkxxzr0KXAp8FFgEHAuc5pzbEI2TnTdxIAGDpm7qawsrWLe7IhqHFhEROUivbGDmnNsKWIt1LwODW6x7BngmFnnom5HM8cNymL9+L+Gw4+v/eJ+dJTX86aoTOHVMv1icUkREPqSOppJ1j/vI6H5s2ltFYXktd1wxlazURK65byH/XKQu2CIiEj0K1t1w6pj+ALy9sYgxeZk89ZVZTBycxc2PL+O/nlhObYManomISPcpWHfD5KHZJAaN97eVANA/M5lHvngKl580jBdW7KK4qv4I51BERI4GCtbdkJIYZMKgLN7fVnrQup9dNpmXbprN4JxUAN7esI/eNGGKiIjEFwXrbpqan8uawnKq6w8eZyUvKwWA+ev28ul73uXqvyxkZ2lNTPIQCjtWFZSzvbg6JscXEZEjS8G6m6YMyyHsYNmOsja3nzKqL984aywLNhVx7m/mc8+bm6hvDEft/Fv2VTH3929ywe/f5CO/eI2vPPgeVXXdGqBNRETijIJ1N03N92ODN68Kby4xGODrZ43hma+dypi8DH7y3GrO+918Kmobun3usuoGPnvfQrYVV/Pd88bz2VOG89zyXXz5wfcIh1XtLiJytOiV/azjSX6fNPqkJ+1vZNae8QOzePxLM3lq6U6W7SgjMyURgMKyWgZmp3Tp3D97cQ1bi6q575qTOGP8AAD6ZSTzq3nr+Mtbm/nCR47p0nFFRCS+qGTdTWbG1GE5LN3edsm6uUDA+Ni0odz2UT+hV1FlHWf+6nU++5eFLNxc3KlGaBv2VPKPRdv42LQh+wM1wJfPGM0Jw3P57cvrKVFrdBGRo4KCdRRMHJLNnoo69lV2bj7ThGCAq2eOYOHmYj551ztcfOdbPPX+ThpCh7+n/afXNxIw4xtnjT1ofTBg3HL+eCrrGrn7zU2dyo+IiMQnBesomDAoE4DVu8o7tV92aiLfPW88b98yh++cO47d5bXc9MhSHluy45D7ldc28OyyAi6YNIhhfdJabT9pRB8+MqYf97+1pdMXECIiEn8UrKPg2EFZAKwq6FywbpKbnsRXzhjNmzfP4XeXT+HiKX6I85dWFvKJP7/NPxZuO6hK+4Xlu6hrDHPZtCHtHvMbZ4+lpiHEgwu2dSlP0dIQClPXqJHcRES6Qw3MomBYbhoZyQmdLlm3lJQQ4OIpBwJwdX0j24qrueWJ5dz61ApOHtGHcyfm8cwHBfTPTObU0e1PGDItP5fjh2bz8MJtfOWMUSQEe/a6bP66vdzx6nqWbC0h7GDCoCyunTWCy6YNJRCwwx9ARET2U8k6CgIBY/zATFZ1M1i3dOnUobx9y5n8/fPT+fTJ+WzcW8ltz6xiybZSLj5+MGsKK/j3ykLK2+kGduWM4RSW1/Lqmj1RzdehOOf4xYtruPovC9m8r4orTs7n2lkjqKpv5DuPLeNLf18SlW5rIiIfJipZR8mEwVk8+O42ahtCpCQGo3bcYMA4dUw/Th3Tj//+6ERuf3419/xnM5dMHcLDC7fx4LvbCBgcNySbafm5TBmWw5xjB5CVkshFkwfzk2dX8fd3t3HOxIFRy9Oh/OHVDfzx9Y2cN3Egv/zk8WQk+4/YrRccy29fXs8fXtvAJ/78Do9cfwrZaYk9kicRkd5OJesoOXZQFqGwY/3uypidIxAwVheWMzArxc/udd547v7MCVw1YzgNIccD72zhpkeWUlbtS65Lt5cyrE8a89ft5fW1e2I+C9gzHxTwq3nrmDN+AHd8eur+QA2+5fu3zx3H7y6fwtrdFXzur4taDdHaU2rqQ2zZV8WqgnKKq+o1bruIxD2VrKOkqZHZ6l3lTBqaHZNzVNU1snBzMR+bOhQzIzs1kXMmDtxfaq6qa2TFzjKG5voJRJZuL2VlpNHbNfctImAwvG86v/j4ZE4a0YeSqnqWbC1heN80hvVJ61aNwJ7yWr7/1ApG9U/nzk9PI7Gde+QXTxlCVV2I7z25nO89sZzffGoKZrG/hx0OO15YUcg/Fm1jwaYiGkIHAvTQ3FTOnpDH1aeMYGS/9JjnRUSksxSso2RcXiZmsLowuvetm3trwz4aQo7Tx/Vvc3t6cgLTj+m7/+8bTh/FVTPyOfvXb2BmnDtxIBv2VJKd6qufl24v5QsPLN6fvk96EgMyk/ncrJF88qRhOOd4eOF2BmQmk5ueRG5aIn3Sk8hKSWzVSOx7T66goraBv37uZFKTDh30Pz09n3W7K7j/7S2cMKIPn5kxvKsvSYesLCjj5seWsbKgnH4ZSXxs6lDGD8okNTFIQWkNb28s4v63t3D/21u4+PjB3HL+sV0eVa6rnHMUVdVTXFVPY8iRkhhgYHYKaUn6ioqIgnXUpCYFGdE3nbWFFTE7x+vr9pIQMGaNab8VeEuZKYlcfnI+v315PVdOz2dMXub+bVPzc7j/2pPYVlzN9uJqCsvr2F1eS2KCD8TltY1878nlrY45c1RfHrpuBgA/fHoFy3eU8f72UiYOzuLfKwt5a8M+bjhtFIGAsWlvJQWltaQnB8lITiAlMUhyQoBvnzOOZTtK+fEzq5g0JJspw3K6+eq05pzjgXe28pPnVpGSGORHF03gipPzW9UgfBM/Icpd8zfyyKLtzFu1mxvPHMPnTh3Zbg1BNGwvrubZZbt4d3MRH2wvpaS6dcO7ITmpTM3P4cThucwZn0d+39b96mOpMRSmuLqeytpGAmYEA0ZiMEDfjKSYvjYicjAF6ygal5fJoi3FMTv+G2v3csLwXLJSOtcw65IpQ/jty+t5aulOvnPu+P3rc9KSOH3cgHb3S08K8u9vzGZ3eS0l1Q2UVNVTUl3P4OzU/Wn2VdSxfGcZhr8FsLKgnKRggK+cMRqAR5fs4E+vb2x17Ae/MJ07r5zGOb+ZzyV3vkVmcgLJkUCekhjgipPz949tfu19Cwk5SAwcCBbD+qRxy/n+uTz47lbWFVYQDARIDBoJQSNgRlFlHQ8t3M5xg7M4Y/wAwg4eXriNgBkBg1NG9WP0gAwqI7cXpuXnMjA7hafe28lPX1jDQwu3cccVU5k8NIc31u2lpr4RM9u/f2IwwOyxvpZjw55KCiJToDbV6hvG5GHZZKUkUlJVz+pd5YTCjne3FPPSil2s31MFQF5WMieNyCU7NYmctEQSAgFqG0IUV9ezvbiaN9bt5dllu7jtmVUMy03lxBF9OGlELmPzMgkGjElDskkIBthbUceuspr9527Ky7DcNLLTEqltCLFhT+s2FTlpiQzKTmXT3kpeXFnImsIK1u2uYFdZLZW17bcryE5JoE9GMgOzkjl2UBbH9M+gf0YyiUGjX0YyweCB2pcBmSn0z0wmHHZt1j6lJyUwInILYvO+KqrrG3FhKK1poLiqjuKqepISAhRV1rOvqo7S6gbqGkM0NDoCBsGgkZmcSFpSkGDACIUdGSkJZCQnkJmcSGZqAkNyUhmTl0lmcgIb9lbSGDq4rUJSgjF6gL+YLSitoaym9cVTfm4qYfxtn61F1dQ2hqlrCFHX6McTSElIIBg0GkJhiirrCASMAEYw4NucpCclkJedQkLA2FNRh3OOhECAYMBICBjBoDE+L4uEoFFW4y+SgsEAQTMSg4ZzkJedQlpSkIraBrYVVRNyEHaOcBhCLkxmciK56UmEwmE27q2iMeRwztHoHOGwIxgwBmSm4HAUV9VT2xAmaICx/7M9NDeN5IQgVfUNFFc1EIhsM/yyf2YyOelJuLDzU/OagXOE8XlJSQiSl5VCyPntdQ0hwvgL6FDYAcaQ3FTCzlESuRhs+txa5HObl5VCSmKQusYQu8trMYtsj+Q1JyWJPhlJGLC1qIqWLU+SgkGGRG4J7i6vpbq+ZZsdR36fNIIBo6ymgaLK1sMz98tIIjMlkYZQmB0lrac4Tk8OMiAzJaozKbZHwTqKxg3M5MWVheytqKN/ZnJUj11QWsPO0ho+ceLQTu87ol86xw/L4an3C/jW2eM63M85IRhgbF4mY5uVxluac2wez68o5MeXHMdV0/OpaQhRVXfgS/GJE4ZyQn4uVfWNVNWFqG0IUdsYYkS/dAZlp/K9C8bzX0+sIDM1gWn5uZEfvTCZKQc+moXlddQ1hGgMOxpDYRrC/kemyfx1e5m3ajdtTTR22bShTBqSxW3PrGq17ZefOJ7RAzIoqarn5seXtdq+s6SGS+58i8/NGsnLq3ezpejg+cKzUhJYdtu5APx9wVbuf3tLq2M8+eWZTM3PZdGWYq7/25JW26+cns9PLjkOgJH/9Xyr7VPzc/jgh+ewoqCMrz70vq8FKdnJk+/v3J/m9kuP46xj8/jXBwX85LnVrY5xxxVTuej4wewoqebCO/7Tanu/jCRq6kNUtfox89KTgnz3/PGEw45/fVDAe5EZ5spqGymrbWTzvire2XToi9RZo/oye2x/EoPG/zzbOo9DclI4dXR/9lXWsWBTUbt5OZRgAMIODtdeMODjSqsf98SgD9YNoTC7ymoO+hyLHEpBi9+GWDC1hO0cM1s5YcKECStXrmy17YXlu7jhwff4++enc2onqqo74l8fFHDjw+93+dj3v7WZ255ZxaNfOoWTRvSJSp5CYcdZv36DUNjx6rdO6/LAKz9/cQ1/en0jt100gWtmjexyfsJhR21jiJsf/YBnlxdyzczh/OiiiVTVhyiraSAcdrimUohz9MtMJislkfrGMBv3VhJ2fntTmsZQmHvf2szzywsZkJnMNTNHMDU/J7Ld/+jPjAxMs353BTtKanCREND0tRo7IJN/LSvg3jc3UVzdQF5WMudNHMjJI/uQEAwwqn8Gowdk4JzjxRWF+wNI0/45aYnMipxjwaYi9lb4UuaqXeWsLChjXWEF9ZESYr/0JPKyU+iTnkRuWhLJCb7ENjQ3DXBsL65h0ZZiCluUMgZnp3DyyD5MHppDdX0jfdOTSUw48F6mJAa4cLIfVW/J1hI276tq9drPHNWH8tpG3ttawtsb97Gv0t9/31fpa2M68jOTnBCgX0YySUEjLTmBrJQEslITyU5JJDM1kUunDiEvK4WdJTWs3V3eqmHirNH9GJKTSmF5Dc8t20V1fYiapkdDiNTEIJmpiZRV17OioIzqOl/aC4f95wH8rHWJwQBVdQ3UN/pSaDBgBAK+ZDsuL5Ps1ERqGkPsLqslMRggKSFAYtA/jhucxZi8TBpDYV5ds4eQc/s/L2EXJic1iQmDfc+RJVtLKKtpIBT2pc1wpNQ5eWgOjWFfkisordlfYg6FfYkzv28a2SmJ1DeG2FJU7Uu8BkEzzIy8rGRG9EsnYMYH20txDgKBSMnYjPTkICcO7+Pb2BSUs7fSl/DDAM6XficNzSEpIcCe8lo27a3y343I59I5x5DcNPqmJ1HXGGLFzvL9n/umvOSkJjImL5OAGat3lVPbENq/zfBdUk8Z1Q8z2FJURUFpbeSDDw7/mk0YlEVWaiKlNQ2sLig78N2I/DMwO4WhuWk4HIu3lOzPQ1OitKQEJgz2DX/X76mgrLoBX6Y/YEp+DolB/zy3Flfvr5FqMrJvOv0yk6lrCLFsZ1mrz2yf9CRG9c/g5ivOYufm9auccxMP/0nvGgXrTjpUsN60t5I5v3qD7889NurTU/7o6RX8bcFWlt127kFdojpqX2Ud029/hU+dNIzbL50UlTw980EBX3v4fX76sUlccXJ+l4/TGApz+d0LWLajjMdvmNnl1vSNoTDfevQDnl5awOdmjeQHFx4blZbm/15ZyA+fXklheS1ThuXw9TPHMHtsf4KHqKHYuLeSv769hceW7KC6PsTxw3K44bRRnDMhL6ojuNU1hli0uYRFW4pZtqOUVbvKI9WrrdMGA8bIfumMzctg4mDfTuC4Idn7GxzGSmMoTFVdiKr6xv0B1CK3ERKD/rZGbnoS6UnBHukZIBJtEydOZNWqVTEN1qoGj6LhfdNJSQzEpJHZ4q0lHDsoq0uBGnyJ4dTR/Xhu2S5+dNEEkhO6N3CLc447X9vAwKwUPnaIMco7IiEY4PdXTOWC37/JVx56j2dvPLXT9+VjFagBzpk4kFmj+/HAO1u5e/5Grr1/EXlZycwZn8eUYdnkZaUQDJgv8RaU8+b6fazaVY4ZnDk+j8+fOpIZx/SJSSBKTgjuHzSnSX1jmMKyWmoaQjSEwiQEbX9p+0g0CksIBshOC2gQHJFuULCOomDAGDMgk7W7oxusK+saWb2rvNtdnD42bQhvrNvLa2v2cN5xg7p1rFdW72FNYQU/vLD7gR9gcE4qv/7k8Xzu/sV8/eH3ufvqEzscWOoaQ9z48Pu8tHJ31AN1k/TkBG44fRSfOWU4Ty/dydNLC3hsyXYeXth6opSxeRncOGc0l50wlOF9e77fdlJCoMdbjYtIbClYR9m4gZk8u6yAUKTVZTR8sL2UsIMTunmv+dyJA8lMTuCxJTu6Faydc/zhtQ30TU/qVvV3S3PG5/Gts8fyq3nr+K8nlvOLyyYftsq4ur6RL/5tCW+u38dXzhjFt88ZF9Oq1IzkBK6cPpwrpw+ntiHEut0V7KusIxyG3PRERg/IjHm1soh8+ChYR9n4gZk8tiTM1qIqjumfEZVjLt5SAsAJw3O7dZyUxCAXHj+Ify7e0a0W629vLGLp9lK+c+64ww6A0llfnTOavZV1PPDOVqrrG/n1J6e0O7Lapr2VfPnB91hTWMEt54/nS6eNimpeDiclMcjkodHvHy4i0pJGNYiycQN9N6do3rdevLWYQdkpDMlJPXziw7hs2lBCYcfTS3cePnE7/vDqBjJTEvjMKdEfeczMuO2iidxw+iieX17IRXf8h/+s33fQ+N3ltQ387uX1XHjHf9hWXM3vLp/S44FaRKQnqWQdZeMH+q4CaworOH9S9+4Lg+8e9f620naHGO2sE4bnMqJvGo+/t7NLLdaXbC3mnU1FfG3O6E43AuuoQMD47nnjGZeXyU+eW8VV977LoOwUJg/NprS6gfe2ldAQcpxyTF9+fMlxjB4QnRoMEZF4pWAdZf0zk+mbnhS1kvXawgoq6xo5sZtV4E3MjI+fMJRf/nsdK3aWcdyQznWT+sOrG0hLCnJtN/pDd9QlU4cw59gBPLJwO6+u2cPS7aWkJyVw2bShzJ08iFNH91NXHxH5UFCwjoFxA6PXInzJVj8y1IlRGsgE4GPThvKbl9fz9wVb+dllkzu834qdZby2di/XfWQkfdKTopafQ8lKSeS62cdw3ezo9lsXEelNdM86BsYNzGRLURU1XRgysaXFW0tISwoyfmD7Q3521uCcVM6bOJAn3t9JUWVdh/e787UNJCUEuC7KA76IiMihKVjHwPiBmTjnh7jrrsVbSpian9PloTzb87lTR1DfGOahd1v3E27Lsh2lvLCikMtPGsaArJ6dPlJE5MNOwToGxjVrZNYdhWW17Cyt4YTh0asCbzItP5fjh+Xw13e2UFXX/sxK4PtV/7/nVpOeFORrc8ZEPS8iInJoCtYxMDYvA7Pud99asjU6/avbYmbcdOYY9lXWc9f8TYdM+9LK3by7uZgbTh8V9dnERETk8HpdsDazi8xsgZmVmNk8M2tz4HTzbjKzNWZWZWbLzeyzPZHHtKQE8vukdTtYL95ajJmfJjEWTh/Xn1mj+3L3/I0UltW2maa4qp7vP7WCITmpfP5U3asWETkSelWwNrO5wBPA34AZwD5gvpkNbiP5N4EfAt8CxgL/B9xvZuf1RF7H5WV2uxp8ydYSxuVlxqw/s5lx6wUTqGsM8/2nltNyBjbnHN97YjlFVXX8+pPHR320MhER6ZheFayB24DHnXN3OufWAtfin8P1baS9HPirc+4559xO59zvgU1AjwTr8QMz2VdZ16nW1s1V1zeysqCcE0dEvwq8uQmDs/jqGaN5efUefvbimv0BOxx2/OyFNby4spDrZx/D9GP6xjQfIiLSvl7Tz9rMRgInAr9sWuecqzWz14G5+EDe3FZgUrP904BMYGms8woHGpmtLaxg5ujO3+ddur2UUNhxYgwal7V001lj2bi3krve2MTyHWXMGt2PN9btZeHmYuZOGsTN546PeR5ERKR9vSZYA00DUbdsDbUNXyXe0m3AM2a2CHgAuA5YADzckZOZ2cp2NnVoEOqmMcLXFFYwc3S/w6RubUmUJu/oiGDA+P3lUxmbt4F73tzM2xuLyExJ4DvnjuNLp42K2uxhIiLSNb0pWA+MLEtbrK8E2opo64CngS8A3wEG4e93d3+kkg4Y0TeNpIRAlxuZLd5aQl5WMkNzuz95R0ckBAPcdNZYbjh9FPsq68nLTI56324REema3hSs90WWWS3WZwOFbaR/CDgWGAkUA58E/gwcA1x9uJM559prZb4SmHC4/ROCAcYMyGBNF4YdDYcd720rYfaY/j0+9nVyQjAqs3uJiEj09Kai047IckSL9cObbQPAzEYBlwF3OOf2OudCzrmHgV8DV5lZ9MbuPIRxAzNZv7uCcNgdPnEz6/ZUUFHb2CNV4CIiEv96U7BeG3lc3LTCzNKB04HHW6RtatHVsil2Q2Rd252Ko+zYgVlU14fYWlzdqf0W9+D9ahERiX+9Jlg736foduBKM7vOzMYB9wFlwH1mdo6ZFUX6Ua8G3gF+YmZzzWyQmV2M73t9t3OuoSfyPGmon35y2Y6Wt9kPbeHmYtKSgkwc3LLGX0REPox6TbAGcM49gO9TfSPwFpACzHTOlQJhfKm5IRLYL8SXuH8PbMQH+tvxg6T0iElDsgkYvL+t48HaOcfCzcWcMDxXDbxERAToXQ3MAHDO3Qvc28b6l4HBzf4uxgf1G3sudwdLT05gbF4mS7d3PFhvL66hsLyWq2bkxzBnIiLSm6joFmNT83NYVVBOXWPHeoy9u7kIQCOGiYjIfgrWMXb80BzqQ2FW7+pYF65FW4pJSggwOXK/W0RERME6xqZEZsx6f1tJh9K/u7mYqcNySE7QpBkiIuIpWMfYmAGZpCcFO9TIbGtRFVuLqpk5qvPDk4qIyNFLwTrGggHjhBF9eGdTUaspKFt6fe1ewM8zLSIi0kTBugfMHtOPvRV1rD3M0KOvr91D3/QkJg3R/WoRETlAwboHzB7rS8rz1+1tN01tQ4h3NhVx2tj+BDTLlYiINKNg3QPGDMggLyuZN9fvazfNgk1F1DaEOU1V4CIi0oKCdQ8wMz4ypj/vbi6mtqHt/tbzVu0mGDBmj1GwFhGRgylY95DZY/tT3xjmrQ2tS9cNoTDPL9/FrNH9yE1POgK5ExGReKZg3UPmjB9ASmKAp5YWtNr2+tq9lFQ38NHjB7exp4iIfNgpWPeQjOQEzpkwkH+vLKSkqv6gbQ+8s4WslAQumDTwyGRORETimoJ1D7r6lOHUNYa5/+0t+9ct3V7Km+v38amThpGW1OvmVRERkR6gYN2DThzRh+kj+/B/b25ia1EVdY0hfvSvlWQkJ/DF00Yd6eyJiEicUrDuYf/v0uMIhR0f//M7XHrn23ywvZTvXXAs/TKSj3TWREQkTilY97DRAzL56+dOJislgT0Vtfzgwgl8errmrhYRkfbpJukRMOOYvrzyrdOPdDZERKSXUMlaREQkzilYi4iIxDkFaxERkTinYC0iIhLnFKxFRETinIK1iIhInFOwFhERiXMK1iIiInFOwVpERCTOKViLiIjEOQVrERGROKdgLSIiEucUrEVEROKcgrWIiEicU7AWERGJcwrWIiIicU7BWkREJM4pWIuIiMS5XheszewiM1tgZiVmNs/MJh4mfY6Z/dDM3jKzqsOlFxERiTe9Klib2VzgCeBvwAxgHzDfzAa3k342sAK4ALgXOB5Y0zO5FRERiY5eFayB24DHnXN3OufWAtfin8P1LROa2SjgWeAFYJZz7i/OuQ3OuVBPZlhERKS7ek2wNrORwInAk03rnHO1wOvA3DZ2+W9gD/BVBWgREenNEo50BjpheGS5qcX6bfgq8f3MLAm4BPgj8DMzOxt/YfIM8N/OuerDnczMVrazaVRnMi0iItJdvaZkDQyMLEtbrK8EclusGw2kA18AtgJXAv8LfBG4L4Z5FBERibreVLLeF1lmtVifDRS2WNeU5kfOuTsi///AzPoDPzezrznn9hzqZMJa1ugAACAASURBVM65NluNR0rcEzqebRERke7pTSXrHZHliBbrhzfb1mRrZLm9xfoVkeWw6GVLREQktnpTsF4beVzctMLM0oHTgcebJ3TO7QI2Ame2OMZxQEPkOCIiIr1CrwnWzjkH3A5caWbXmdk4/P3nMuA+MzvHzIrM7LzILt8FrjOzr5pZvpl9NLLuf51zlUfkSYiIiHRBb7pnjXPuATNLBG4Cfgq8Dcx0zpWaWRiow5eccc49bma1wPcjabcAPwF+eyTyLiIi0lW9KlgDOOfuxY9G1nL9y8DgFuueA57roayJiIjERK+pBhcREfmwUrAWERGJcwrWIiIicU7BWkREJM4pWIuIiMQ5BWsREZE4p2AtIiIS5xSsRURE4pyCtYiISJyLarA2s7OjeTwRERGJ/nCjNwPzAMxsI/Ae8H7TwznXct5pEREROYyolqydc81L1rcA64CPAPcDO81sl5k9H81zioiIHO2iWrI2s7Odc/MAnHOPAo822zYYOAGYGs1zioiIHO16qhr8PWCpc+4Z4Jkon1NEROSo1lPV4H/lQDW4pqwUERHphJjNZ61qcBERkeiISrA2MwMuAaYAZcByYIFzrqIpjXOuAChA1eAiIiKd0u1gbWbZwCv4ErNFVjugyMxuBe5xzrnunkdEROTDKhr3rH+OD9S/Bk4Hzge+AxQDfwbuiMI5REREPrSiUQ1+LnCXc+47zda9ZGa/BX4CfNfM3nTOPRKFc4mIiHzoRKNknQksbLnSORd2zn0PeAm4OgrnERER+VCKRrBeBsw4xPanD7NdREREDiEawfpe4PNm9ql2tucB26JwHhERkQ+lbgdr59yDwCPAQ2b2rJldbma5AGY2G7gJeLy75xEREfmwiko/a+fcVWa2EPhv4ALAmVktkAJsBt4zs3xgu7pxiYiIdE7Uhht1zv0e6I/vunU3sAUIA8fgB0LZDFSb2Soz+1e0zisiInK0i+pwo865Rnzr75cAzCwFOA44PrLMB4bihx0VERGRDojZ2OAAzrlaYHHkISIiIl0Q1Vm3REREJPoUrEVEROKcgrWIiEicU7AWERGJcwrWIiIicU7BWkREJM71umBtZheZ2QIzKzGzeWY2sQP7TDazejN7rSfyKCIiEk29Klib2VzgCeBv+Jm89gHzzWzwIfZJAO4HKnoijyIiItHWq4I1cBvwuHPuTufcWuBa/HO4/hD7/BeQBDwX++yJiIhEX68J1mY2EjgReLJpXWSEtNeBue3sMwm4Ffg6EIp9LkVERKIvpsONRtnwyHJTi/Xb8FXiB2lW/X2Xc+4VM7uqMyczs5XtbBrVmeOIiIh0V68pWQMDI8vSFusrgdw20n8XSAduiWWmREREYq03laz3RZZZLdZnA4XNV5jZccD3gDOcczVdOZlzrs1W5pES94SuHFNERKQrelOw3hFZjgCWNFs/vNm2Jl8HUoFXzKxpXQpgZlYBnO+c+0/ssioiIhI9vakafG3kcXHTCjNLB04HHm+R9lZgDDCl2eNf+Kk6p6ApO0VEpBfpNcHaOeeA24Erzew6MxsH3AeUAfeZ2TlmVmRm5znn9jjnNjZ/AOVATeTv2iP4VERERDql1wRrAOfcA/g+1TcCb+Grtmc650qBMFAHNBy5HIqIiERfb7pnDYBz7l7g3jbWvwy0O5KZc+7aWOZLREQkVnpVyVpEROTDSMFaREQkzilYi4iIxDkFaxERkTinYC0iIhLnFKxFRETinIK1iIhInFOwFhERiXMK1iIiInFOwVpERCTOKViLiIjEOQVrERGROKdgLSIiEucUrEVEROKcgrWIiEicU7AWERGJcwrWIiIicU7B+kgo3gSPXAX3Xwjr5x3p3IiISJxTsO5pFbvhnrNh3UuwewU89ElY//KRzpWIiMQxBeue9uJ3oaYYPvsMfO09yB4GT38F6quOdM5ERCROKVj3pL1rYeWTcMK1kD8D0vrA+b+AykJ4964jnTsREYlTCtY96d27wAIw6+sH1o09FwZNgUX3QKjxyOVNRETiloJ1T2mshxWPw6g5kDv8wHozmP5FKN8J6186cvkTEZG4pWDdUza9BrWlMPlTrbdNuBgS02D5oz2fLxERiXsK1j1lw8u+CnzMOa23JaXD+Lmw9kU1NBMRkVYUrHvKxldhyImQmtP29vEXQmMNbJ7fs/kSEZG4p2DdE0q2QtEGGH1m+2lGnQEW1CApIiLSioJ1T9j4ql+OmtN+mpRs351rwzxwrmfyJSIivYKCdU/Y+CokZ8PgaYdON/osKN0G+9b1TL5ERKRXULCONedgy39gxKkQTDh02jFn+6WqwkVEpBkF61gr2eyHFx120uHT5h0HGXm+m5eIiEiEgnWs7Vjil0NOPHxaMxg+C7a9C+FQbPMlIiK9hoJ1rO1c7PtXD57asfTDZ0J9BRQuj22+RESk11CwjrUdi6H/sZCc0bH0w2f55da3Y5cnERHpVRSsY6mxDgqXwdAOVIE36T8eUnNh61uxy5eIiPQqvS5Ym9lFZrbAzErMbJ6ZTWwn3XQze8LMtjVLO6lHM1u4AkL1nQvWgQDkz/Qla/W3FhERelmwNrO5wBPA34AZwD5gvpkNbiP5FcBC4FxgNuCA182sTw9lFwre88shJ3Ruv+EzfQvyvWuinycREel1elWwBm4DHnfO3emcWwtci38O17dM6Jy7yTn3M+fcaufccuBLQB/gkh7L7e4VEEyCfmM7t1/+KX65bUH08yQiIr1OrwnWZjYSOBF4smmdc64WeB2Y24FD7IssM6OeufbsXgX9xkEwsXP7DZwECSmwfWFs8iUiIr3KYYbUiivDI8tNLdZvw1eJH06kmTXvd+RkZraynU2jOrI/zsGe1TD+gg4lP0hCkq86366StYiI9KKSNTAwsixtsb4SyD3UjmYWBH4CLHLO9cwclKXbfH/pARO6tv+w6VC8CSr3RDdfh7LySbhrNtw+BO4+A9Y813PnFhGRdvWmYN1UjZ3VYn02UHiYff8PGIm/x90hzrmJbT2AjR06wJ5VfpnXZmP1w8uPVBZsf7dr+3fWKz+GR6+BmhIYPxcqd8M/Pg2v/r+eOb+IiLSrN1WD74gsRwBLmq0f3mxbK2b2C3zL8LOdc+1VbUff7sipulqyHhoZS3zbAjj2oujkqT3vPQBv/hLGXwiX3QOJqVBfBU9cD/N/AZl5cNIXYpsHERFpV28qWa+NPC5uWmFm6cDpwONt7WBm3wa+CVzhnPtPD+TxgD2r/BzVWW31KuuAtD5+gJRYl6xLt8PzN8OgKfDxv/hADZCU7v8eehK8dCvs1bSdIiJHSq8J1s45B9wOXGlm15nZOOA+oAy4z8zOMbMiMzsPwMxuAH4OfB14zcyyI4/0Hsnw7lUwYKKfnKOrhk2HgqXQUBO9fLX0yv9AqA4u+RMkJB+8LSEZLr3Lj23+5PUQaoxdPkREpF29JlgDOOcewPepvhF4C0gBZjrnSoEwUAc0RAZJuRP//P6Ab5TW9Hg25hltrIOi9ZDXxSrwJvkzINwABR1qwN55O9+D5f+EaVe3n9e+o+Ds//F5WPpgbPIhIiKH1KuCNYBz7l7n3CTnXD/n3Eedc9si6192zg12zr3inCtwzgWcc9bG44yYZ3LfOgg3dv1+dZNh0/0yVoOjvHY7JGXAGbceOt0J1/oq+dduh7rK2ORFRETa1euCda+wu5stwZv0OQbS+sVmcJS9a2HDPJj2WcgYcOi0wQQ4+8dQWQjv3Bn9vIiIyCEpWMfCnqaW4Md27zhmvip8+7vRn9Tj3T/7e9HTW43U2rYxZ/vpOxfcCbXl0c2LiIgckoJ1LOxeBdnDfGvw7ho23U/qsW9994/VpLoYlj4M4y6A3BEd28cMZn8bastg8b3Ry4uIiBxWb+pn3XvsWQV5x0XnWPsHR1kA/Ts5IUh7VjwOjTUw/Yud2++YM2DwNHj7D3DyFyEpLTr5OZRwGDa+Auv/7S9YktL9rGTHXgQ5+bE/v4hIHFDJOtpqSqB8Z/dbgjcZdDwEk2FbFPtbL33QB7rhp3ZuPzOY/R2o3ucHUom1HYvh/86ABz8Oi+6Bsu2+sd1L34PfTYHnvu1rCUREjnIqWUfbntV+OaCbjcuaJCTDkGnRm9Rjz2rfDeu070KgC9dqY8/zz+2t38GJ17bumx0tSx+Cf33Nl6TP+zkc/ylIzfX37nd94M+/6P9g3Yvwib/C0E7OGS4i0ouoZB1tTcOMRqtkDf6+ddEGqNp3+LSHs/Qhvzz+8q7tHwjAR74JFQXwwT+6n5+2vHsXPHWDv5Xw5Xdhxpd8oAZfuh88BT5xH3zmSWiohvvOg9Wx7z4vInKkKFhH255VEEiAvmOid8xoTeoRaoRlj0D+TN8trKsmXgp9RsF/fhP9Uc2WPwYv3Oyr6K99HrIGtZ921Bz44nzIHQn/vBo+eCS6eemI+ioo2epnSKtpOSGciEh0qBo82navgn5j/ZzU0TL0ZL/c/q6fEaurNr3mZ9Oa84Pu5ScQhFO/Af/6qp9Wc/Inune8JgVLIyXqSXDFQ74K/HCyh8K1L8DfL/VDojZU++r5WHEOtr4NKx6DzfN9jUdzmYN8TcixF/nW9j3RCE9EjnoK1tHknL8nPObs6B43va8vqXe3kdkHD0NCKky4+PBpD2fyp+D1n8Gbv4LjLuva/e/masv9FJ2JaXDFw53r9pbeFz77DDz4CXj2JsDBiZ/rXn5acs7fH3/9p/6euQV9UB57np+sxQJQtRcKl8OGV2DVU77q/oRr4KTrIHtIdPPTEeEQhBp8TU9QX3WR3kzf4Ggq2wF1ZdG9X90kfzose9SPO96VRl11FbDmeRh/AaS0nBK8CxKSYNbX4YXvwNrnujeNp3PwzNehZDNc8Q/IGdb5Y6Rkw1WPw98vg2e/4Y950ue7nqfm9qz2+dv+LqT392OlH/9pyOjfdvqGWt/VbOHd/lbB23+AyZ/0r1f/cdHJU3PlBbDtHT/W++6VULrN90horD2QJiUHMvJ8978BE2DgZBhywqFvM8SCcz5fFoBAYvcv8kQ+JBSso2lPZJjRaLUEb27YDHj/776qOH965/df/azvWz35U9HL07TPwPz/hfmRubC7OsPYkvtg5RNwyldh3Pldz09y5oGA/dw3/bqOBGznfCnUhf0D55ehBljwR3jz1750Ovs7MOUqP41ouBHKdx0IdtXF/oKIyEhzAyfBR++A6iIftJc+7LvMjTzNj7U+MNIPP7WPrxlwzo8pv3+kusgyIQX6jPT/L9nix2YPh2D3Cl8dv+0dKN544LlkDobswTBgvK+lsARISvX7le2AzW/C6mcOpE/rC3mTYeQs34c+IfnANKlNAokwaHIkD1t9172W+h/rb48UrvAXNeUFkcdOv2yo9kG6vuXY8gEIJh6YTja9v89DcqYfajetL6RHlqPmQGKKf93LC1rnoe8xvjajvspfYIXD0FDpn3t9hT9XQpJ/r/as9nkJh8CF/BL8+YOJPr+hRv++N9VMBJP9e5GS7d//mjKf14Rk/z4lJPuJb3JH+OMVvOfThSOfKxfy70nuCL9+7xqorQDX6P8ONfr/9xnl/64ohKo9/ljhxkgezb9OyZkQqvdDABPwFz2W4C+C0vtGanuCfkIhZ5HtAb8uKc2PrmgBKN7su5u6sP/suTAQ9s8jkOg/v+U7DzyHpu9G+gB/0d9Y78+B83mzgP8dSM6EnOH+75It/rtkgQPbA0H/HbGgfy+r9kaO3fTZd/51Skr371/Jlsh3wx1YpvX1F6BEajQPGuXRRb47kbY5pdva+OzhayyDif4zUbm79W9Y5iBIzfHPs/n3rElylr8VF8uZESMUrKMpFi3BmzQfHKUrwXr5Pw/84EVLYirM/BrM+wGsebbt0nVtuf9xr6v0X5bGWl/yHDbd/6hseMX3l84a4n/IXr7N1x6MnH0gcD/2eb9fuNF/6cMN/ov80Tv89udv9vePww2RNPX+h+a5b0JdOaTkwkv/dSAYN/0oXXwnTLnC/xD8fkrrvCdl+h/50Wf5AVnm/69/NEnJhlu2+f+/8XM/hGtLX3gFPna3fz5PfwU2v+EfTU79Jpz1I///O09uvf/Qk+ALL/uR4/7x6QOfseZmfMXfetm5BF79sW+p39zH/+JvVexd2/oc1UWw/W3Y/Frr4+5/njnwqb/51+ztP/gx5VvKGOh/7HCtt4EvzQ+e5oPl4r802xD2U7S6kA92RRugdGskMLQhtY8PWnVlrbdlD/eBqKYUKne1/3wOJZDoP0cinVEa+wmOFKyjac+qyJVWF6pxD6fvaF/S2PKWr07tjIrdsOl1fx83mNi5fUON/qq3usg/aop9kGoK+pW7/RXsY5/3pY76Kh8wv7XGb3/zl75PdEtX/wsGT/V9qV3IX73P/8WB7cGkA8F6xyL/Ax2MlHICif5Hu0lypp+MJJDgn18g8rGuLfPBf9RZMPESf8z9V/8B/5qCfz6nfNWva6z1JdbdKwAHH7sHJkUGZaktPXj/hJQDeRh/oS9JQLOrczvwWRg5Gy74pb9o2fqWfzRUwzt/8CWwwSf4SVWSs3ypI1Tv819dDPec5fvGhxv9cx9wrJ8kJu84X5Kc9AkfBDMH+YuelnkYcqL/b+ZAuPTu1u9F7nC/X8F7sOyfvqRUss3XxIB/3n9t5zaHJfiS74AJ/jVOyfUXZen9/fqmz9vAST7foUbfG6Gl9H4wKjIh3sbX/AA4NWX+3LVlvvSXMcB/Fit2+wvAUH3knnwQCPjXIinNP8/c4b5UlpTmLwKTMnwpa+hJPt32hb60GUiIvKfm042/wF/MbVvgS1JNJdpwpPQ75MRIaXEz7FrmLzRCDT4vjXWQ1sd/LsKN/iJwf6k38kjJ9p+7QBCKNvr9LOj/DgR9fkbO9svynVC2069vStP0WqbmQHWJHziIsC/5Ni0z8iA33+dh+7uREnukZO/CEEjyt0Oc8xdGdeU+nxZ5HQj4wZiS0vzrXLLl4M895kveWUP8hc3OJX7d/lJv2H8O+o315yxY6j8T+0vGYZ9+0PF+e9lOX4PQ9D40fW77j/Pfh7oKX+t0UKk3UsOQM8z/f8ciWl0oJqb71wr8hWptG701Bk/1vwsVu6FkU+vtfUb5210NtbBraevtaX39+/n3G2Hv9tbbo8hctCeIOMqZ2coJEyZMWLmyjRLOH0/xgePz/47NyR+/zldhfneLrw7sqHf+6EuWn58Hw5qVrMJh/4NQstlXb1ZEqhfHz/Ultepi+MUxtPoSHHM6XP20///fP+4bXFXt8R/s/uMhOQMu+bP/kdr6jt+enOF/PBPTfAlq4GTfGGzV0zD3Nz4wJ6b4H7pgcnTuZYYafDewxX/x1V0X/K/Pe1vV9fVVflS213/mv9RTrvL3ptP7dj8f7eVt42u++n/rW76ark3mA2F+pDHbyNmtq6ljIRyOXKTt8/376yoiQSPgf9wy8nxQTM3t+u0PkaPExIkTWbVq1SrnXAzugXoqWUdLY72/+pv6mdidY8zZvjp761sw+syO7/fBw/6HddcyPy74jC/7ksfa5+GRKw9OG0j0JeQxZ/sf4ulf8iWa9H6+NJvW1x+ryVWP+R/2e8+CPWv8PeOme6wAw0/xj5be+aMP1DNvhJOi3HK7STARLvyND9DP3AR/u8SXrCZc4kumiWm+BLd5vr8Iqin2V/sXPHrwRU2s8jb2HP8AqNzjPz9Ve31JLiHF3w/vOyY6DQI7KxCAzDz/EJEjTsE6Wvat89VO3Z3D+lBGnQkYrJ/XfrCuKjpQGlzyV3j9dt9QBeD5b/mS0ZizfbAePAXOus0PKpI7wjeUSO1zoFRrBuf/7PD5CgTg4j/CXR+Bp7/qu1EdqmS88VV/n3v4LDjzRx177t0x4WIfsBf8yZee/33rwdsDiTDyI35ykjHnHJkWyhkDDj+vuIh8aClYR8vuFX7ZdI8kFtL7+u42G+YBP/P3gIo3+cFOtr7j79uUboVvr/c//Kk5kBAZlOPs//H3mfuNPdD1K3uoH9wkGgaMhznfh3k/hBdvgfN/3nb16PZF8I+r/P2uj9/Xc/1/U7Lh9Fv8mOh7Vvlq/4Zqf98rb2J0pjMVEYkRBetoKVzul7EsWQMcc5ofiKRoo+8S1NQCOSnTT/gx6eMHujCMmwsv3OIDfGcbpXXFKV/zI7gtvMs3Tpnzw4NLqcse9Q3KUrLg6qeOTBWrWaRxVozfJxGRKFKwjpbC5b46OTkz+scuL4A1z8Hqf/lWquDvPY8915egjzkdBk050Fq0yboXfTeeObe2PGJsBAJw8R98y8///AY2veEnDAkm+fxvmOdbMH/q7wff1xYRkUNSsI4G53w1+PBZ0T3uhld80Nvypv87YyBMvQq2vO1nvJr9nUP3m178F1+9O/Fj0c3XoQQT4ZN/8yX+N3/lW2ODL/nP/Bqc/j2Nly0i0kkK1tFQscv3QY7G/epdy3x/1GCiP+7eNX7Qi4mX+D6egYAfmOLft/o+lsNOavs4RRth4ysw/YaeD46BAJzyZTj5+gMN7/qN7Vx3MxER2U/BOhoKu9m4LByCtS/4oS23vnVgxKlJn/TDg7YcyGTSJ3xr6g8eaj9Yv/FzP5DCydd1LU/REEyIzWhuIiIfMhpFPxp2NzUuO65z+4Ua/ZjRfzjR93feswo+8q0DIzwlJLU94lhmHow518/fXF3cevvetX4kqilX+JGGRESkV1OwjoZdy/z4ydlDO7df+Q4/XrQL+8E7vrEKzvxhx2ZCOvUmaKiCt37beturP45MPHFz5/IjIiJxSdXg0bBzie8e1ZFhF3cu8VMZnnydH4jkmmf9qFqdHbM7f4YffvKdP8JxHz8wK9KyR/1oXKd81Q98IiIivZ5K1t1VvssPWTm0nXvHTaqK4KmvwP/Ngdf+n5+NCmD4zM4H6ibn/8KPt/3Qp/w970X3wNNf9mNJz/l+144pIiJxRyXr7tq52C8PFaxXPgXPfcu3GJ/6GV/VHY3xnnOHwxUPwyNXwcOX+3V5k+DKf/bMZA8iItIjFKy7a8civxwyre3t6+fBo5/1M1Jd/uCBeamjZfhM+OpiP952ciYcc4ZvmCYiIkcNBevu2rHYz4yU1ufg9Y31PmiOOhPO/18/mEms+jun9fHDjIqIyFFJ96y7I9ToG4s1rwJ3zg9a8udZUFvmBwiZfr1G7RIRkS5Tybo7di+HxhoYeoL/u6EWnroBVj4BIz7iS9ciIiLdpGDdHZte98sRs31r739cAdvf9d2mzvrvnpv+UUREjmqKJt2x8TU/L3PfUXD36X4yj7m/hpM+f6RzJiIiRxEF666qq/DTVU76hJ+a8qwfQTgMY8850jkTEZGjTK9rYGZmF5nZAjMrMbN5ZjbxEGlnmNkbZlZqZm+a2WlRy8i6lyBUd2CI0dFnKVCLiEhM9KpgbWZzgSeAvwEzgH3AfDMb3EbaacCrwHzgJGAx8KKZTY1KZhbfDxgs+QvUV0flkCIiIm3pVcEauA143Dl3p3NuLXAt/jlc30ba7wHLnHM/cM6tB74JbAa+3u1cNFTD1jchmASfflTdskREJKZ6TbA2s5HAicCTTeucc7XA68DcFmkTgUtapHXASy3TdknZDr/82N0weEq3DyciInIovamBWdMUUptarN+GrxJvbggQbCdtPzNLiQT6dpnZynY2jd9Y1MjEe4BHbwVuPWzGRUTk6LVx40aAYbE8R28K1gMjy9IW6yuB3E6kBcgBCruYj0BdiPCqneVr2Lmqi4eQbhoVWW48orkQvQ9Hnt6D+DAeiOnsSb0pWO+LLFtOV5VN68B7qLThZtvb5Zxrs5V5U4m7ve0Se3oP4oPehyNP70F8OERNbNT0mnvWQORGMSNarB/ebFuTnYBrJ+0u51xjtDMnIiISK70pWK+NPC5uWmFm6cDpwOPNEzrnaoB5LdIGgPNbphUREYl3vSZYR1pz3w5caWbXmdk44D6gDLjPzM4xsyIzOy+yy8+BmWb2AzMbDfwO6AfccSTyLyIi0lW9JlgDOOcewPepvhF4C0gBZjrnSvH3ouuAhkjaV4FLgY8Ci4BjgdOccxuOQNZFRES6rDc1MAPAOXcvcG8b618GBrdY9wzwTA9lTUREJCbM1y6LiIhIvOpV1eAiIiIfRgrWIiIicU7BWkREJM4pWIuIiMQ5BWsREZE4p2AtIiIS5xSsRURE4pyCdRvM7CIzW2BmJWY2z8zandHGzGaY2RtmVmpmb5rZaT2Z16NVJ9+DjWbmWjy+2JP5PVqZ2UQzezrymn72MGn1XYiBTr4H+i7EgJmdZ2YvmtkuM9trZk+a2fBDpO/w71dHKVi3YGZzgSeAvwEz8NNpzjezwW2knQa8CswHTgIWAy+a2dSey/HRp5PvQRDIx0/oktPscV9P5fdoZWZTgCVEhvA9TFp9F2Kgk++Bvgux81ngaeBUYC5+Bsd5ZpbYMmFnfr86QyOYtWBmi4CNzrnLI3+nALuA3znnbmuR9jFgqHNuRuRvA1YCC51z1/Rkvo8mnXwPRgCbgdzIGPESJWaWCgzDj7m/BbjGOffXdtLquxADnXwPRqDvQo8wsznAK8CZkXkomm/r8O9XZ6hk3YyZjQROBJ5sWuecqwVex19NNU+bCFzSIq0DXmqZVjquM+9BxEigSD9O0eecq3HOrTtcOn0XYqej70GEvgs9Z19kmdl8ZRd+vzpMwfpgTfcgNrVYvw0Y2mLdECDYTtp+kasp6bzOvAcAxwDZkftI683sz2Y2KKY5lJb0XYgP+i70nFmAA5a2WN/Z368OU7A+2MDIsuWVaSWQ24m04O8VSed15j0AeAE4HjgL+BlwPv7+0ICY5VBa0nchPui70APMLAO4FXjMObe1xebO/n51WK+bIjPGmqo2slqszwYKO5E23Gy7dE5n3gOccwVAQeTPD8xsAbAC+Dzw01hlUg6i70Ic0Hch9iK3fB4HGoGb2kjSxPZ2WwAABU9JREFUqd+vzlDJ+mA7IssRLdYPb7atyU58NUhbaXc55xqjnbkPic68B60451YCxYCq/3qOvgtxSN+F6DKzAPAAcAJwbuTiqKVu/X4dioL1wdZGHhc3rTCzdHxXiMebJ3TO1QDzWqQN4KueDkorndLh96AtZnYM0Ad4P0b5kxb0XYhP+i5E3e+BjwJznXNr20nTrd+vQ1GwbibSgvV24Eozu87MxuH7KJYB95nZOWZWZGbnRXb5OTDTzH5gZqOB3wH9gDuORP6PBp19D8zsZjM7xcwGRbpTPI3vl/rgkXoOHwb6Lhx5+i70HDP7KXAd8BlgjZllRx6pZnatme0zs8mH+/3qTh4UrFtwzj0AXA/cCLwFpAAzI90hwvj+jg2RtK8Cl+KvthYBxwKnOec2HIGsHzU6+h6YWRIwCXgY3/ryT8CLwBnOufojkfcPEX0Xjjx9F3qAmZ0C3AIk4UvHpc0ef8S/D7X4+9iH+/3qej40KIqIiEh8U8laREQkzilYi4iIxDkFaxERkTinYC0iIhLnFKxFRETinIK1/P/27iXUqiqO4/j3l4olQkFE0Qt7UmAaGUg0iKgGRjQIc1BmNbk1ylE4VIoIGjYobSAVYgoVNTDCCA2MiMoeJvSgICgrGwT5CE36N9j7yuFw7jXusdy4vx+4LM5aa6+79ujHXmvtcyRJHWdYS5LUcYa1JEkdZ1hLktRxhrUkSR1nWEuS1HGGtaTjkqxOciDJ/CRzkzyW5N0kPyf5Lck7SW4f83+sTVLT/J17su5HOl3MPtUTkPT/SrII2FOjf8VnAthcVQeTPAesovmloVeB89r27UlWVdWmMafyPPDeiPoDY44rnXb81S2pR5K8D1xUVQtGtN0M7AKWVNXuJNcAR6vq+4E+FwLfAn8Al1TVsRnMYS2wDnioql6a0Y1IPeMyuNQTSeYBS4FPpugyAXxcVbsBquqrwaBu6/YBO4ELgCv+u9lKGmRYSz2QZAtwCJgF3NPuDe8baD8HuBfY8C+Gm1yOO76NluT1JPuSHDnBfvSyk3hbUm8Y1lI/bAAml5zXAsuB+wbaHwCOAa9MN0iSM4AlwBHg64Gmj9px/wa+aMdfDjzStu9sP4/ao5Z0Ah4wk3qgqnYkub/9uL6q9g91mQA2VdWhEwz1MM0S+MbB/eqqejrJVcCZwFtV9RpAktvaLm8P1I15N1L/GNZSf9wA/DQc1EluAhbSPF1PqT1c9hTwJ/DEiC63tOXOgbrFbfnZDOYrqeUyuNQDSebSBPKnI5ongA+raspATXIW8CZwPjBRVT+M6HY3zWtXOwbqrm/Lz2cyb0kNw1rqh+uAOQyFdZKzgRVMc7AsyWxgC3Aj8Myo96uTzAfuALZV1dGBpsXA/qr6Zew7kHrMsJb6YVFb7h2qXwn8BWwddVGSWcBmmqfmF6pqzRTjL6PZr35j4No5wLW4BC6NzbCW+mFBW/44VD95sOzw8AXtye+XaV7pehF4dJrxHwR+B7YN1F1J8zT/3YxmLOk4D5hJ/TD5bvTqJBdX1dYkS2meuFcOd05zZHsjzetd3wDbgRVDJ7m/rKq9SRYCdwJPVtXBgfbL29IlcGlMhrXUD+uBW4G7gMtolr0ngA+qas+I/pfSPC0DXE2zFD5sHc2y+hrgMPDsUPu8tvx1nIlL8rvBpd5K8jhNWO861XORND3DWpKkjvOAmSRJHWdYS5LUcYa1JEkdZ1hLktRxhrUkSR1nWEuS1HGGtSRJHWdYS5LUcYa1JEkdZ1hLktRxhrUkSR1nWEuS1HGGtSRJHWdYS5LUcYa1JEkdZ1hLktRx/wDHyAemTzSv/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 487.5x300 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def final_value(s, det):\n",
    "    return s/2/(1+s+4*det**2)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, num='evolution', figsize=(3.25, 2.))\n",
    "ax.plot(obe.sol.t*gamma/2/np.pi, np.abs(obe.sol.rho[0, 0, :]), linewidth=0.75,\n",
    "         label='$\\\\rho_{00}$')\n",
    "ax.plot(obe.sol.t*gamma/2/np.pi, np.abs(obe.sol.rho[1, 1, :]), linewidth=0.75,\n",
    "         label='$\\\\rho_{11}$')\n",
    "ax.plot(rateeq.sol.t*gamma/2/np.pi, np.abs(rateeq.sol.y[0, :]), linewidth=0.75,\n",
    "         label='$\\\\rho_{00}$ (rate eq.)', color='C0', linestyle='--')\n",
    "ax.plot(rateeq.sol.t*gamma/2/np.pi, np.abs(rateeq.sol.y[-1, :]), linewidth=0.75,\n",
    "         label='$\\\\rho_{11}$ (rate eq.)', color='C1', linestyle='--')\n",
    "# ax.plot(obe.sol.t[-1]*gamma/2/np.pi,\n",
    "#          final_value(len(laserBeams)*laserBeams[0].beta(np.array([0., 0., 0.])),\n",
    "#                      ham_det+laser_det) ,'o')\n",
    "# ax.legend(fontsize=6)\n",
    "ax.set_xlabel('$t/2\\pi\\Gamma$')\n",
    "ax.set_ylabel('$\\\\rho_{ii}$')\n",
    "ax.set_xlim((0., 2))\n",
    "ax.set_ylim((0., 1))\n",
    "fig.subplots_adjust(bottom=0.2)"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}