maop_steady-state.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from qutip import *"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [],
   "source": [
    "# Parameters\n",
    "gOP = 2 * np.pi # Optical pumping rate\n",
    "W_R = 0.01 * gOP # Coupling strength in terms of the Rabi frequency\n",
    "w_L = gOP * 2e2 # Etalon Larmor frequency\n",
    "gL = 2 * w_L  # Cavity linewidth\n",
    "gT = 1.0 * gOP # Thermal relaxation rate\n",
    "delta = np.linspace(-6*w_L,6*w_L,320)# Microwave detuning"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [],
   "source": [
    "# Initialize the quantum system\n",
    "\n",
    "states = []  # list with the basis states\n",
    "for i in np.arange(8):\n",
    "    states.append(basis(8, i))\n",
    "\n",
    "psi1_1 = states[0] # (1,-1) ground state\n",
    "psi10 = states[1] # (1,0) ground state\n",
    "psi11 = states[2] # (1,+1)ground state\n",
    "psi2_2 = states[3] # (2,-2)ground state\n",
    "psi2_1 = states[4] # (2,-1)ground state\n",
    "psi20 = states[5] # (2,0)ground state\n",
    "psi21 = states[6] # (2,+1)ground state\n",
    "psi22 = states[7] # (2,+2)ground state\n",
    "\n",
    "# Projector operators onto the basis states\n",
    "P = []\n",
    "for psi in states:\n",
    "    P.append(ket2dm(psi))"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Relaxation\n",
    "\n",
    "We model the optical pumping as relaxation from F=2 to F=1 at the same rate for all sublevels. The process is described by the following collapse operator:\n",
    "$$    \\hat{L}_{\\rm OP}(m_{\\rm F},\\widetilde{m}_{\\rm F})=\\sqrt{\\Gamma_{\\rm OP}}|F=1,m_{\\rm F}\\rangle\\langle F=2,\\widetilde{m}_{\\rm F}|. $$\n",
    "\n",
    "For the thermal relaxation, we assume that each Zeeman sublevel mixes with every other sublevel at the same rate:\n",
    "$$     \\hat{L}_{\\rm}=\\sqrt{\\Gamma_{\\rm th}}|i\\rangle \\langle j|,$$\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [],
   "source": [
    "# Collapse operators\n",
    "\n",
    "c_ops = [] # List with collapse operators\n",
    "\n",
    "# Optical pumping relaxation\n",
    "for psi1 in states[0:3]:\n",
    "    for psi2 in states[3:]:\n",
    "        c_ops.append(np.sqrt(gOP) * psi1 * psi2.dag())\n",
    "\n",
    "# Thermal relaxation between the Zeeman states\n",
    "for i in np.arange(8):\n",
    "    for j in np.arange(8):\n",
    "        if i != j:\n",
    "            c_ops.append(np.sqrt(gT) * states[i] * states[j].dag())"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Hmiltonian\n",
    "\n",
    "The Hamiltonian is given in the rotating-wave approximation"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Diagonal part"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "The diagonal part of the Hamiltonian is given by\n",
    "$$\\hat{H}_0 = \\hbar\\sum_{m_F}\\left(\\dfrac{\\delta}{2}-m_{\\rm F}\\omega_{\\rm L}\\right)| F=1,m_{\\rm F}\\rangle\\langle F=1,m_{\\rm F}|+\\hbar\\sum_{\\widetilde{m}_{\\rm F}}\\left(\\widetilde{m}_{\\rm F}\\omega_{\\rm L}-\\dfrac{\\delta}{2}\\right)|F=2,\\widetilde{m}_{\\rm F}\\rangle\\langle F=2,\\widetilde{m}_{\\rm F}|,$$\n",
    "\n",
    "where $\\delta$ is the microave detuning from the clock transition, and $\\omega_L$ is the Larmor frequency."
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [],
   "source": [
    "def H_d(d, wL):\n",
    "    \"\"\"Diagonal part of the Hamiltonian as a function of detuning d form the clock transition and the Larmor frequency\"\"\"\n",
    "    H = 0*qeye(len(states))\n",
    "    H += (wL+0.5*d)*P[0]\n",
    "    H += 0.5*d*P[1]\n",
    "    H += (-wL+0.5*d)*P[2]\n",
    "    H += -(0.5*d+2*wL)*P[3]\n",
    "    H += -(0.5*d+wL)*P[4]\n",
    "    H += -0.5*d*P[5]\n",
    "    H += -(0.5*d-wL)*P[6]\n",
    "    H += -(0.5*d-2*wL)*P[7]\n",
    "    return H"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Interaction Hamiltonian (off-diagonal)\n",
    "\n",
    "The interacion Hamiltonian is given by\n",
    "\n",
    "$$\\hat{V}=\\dfrac{\\hbar}{2}\\sum_{\\widetilde{m}_{\\rm F}}\\sum_{m_{\\rm F}}\\Omega_{m_{\\rm F},\\widetilde{m}_{\\rm F}}|F=1,m_{\\rm F}\\rangle\\langle F=2,\\widetilde{m}_{\\rm F}|+ {\\rm h.c.}.$$\n",
    "\n",
    "For the $\\pi$-transitions, the coupling strength is given by\n",
    "$$\\Omega_{m_{\\rm F}m_{\\rm F}}=\\Omega_{0}\\cos\\theta S_{m_{\\rm F}m_{\\rm F}}\\sqrt{L(\\delta_{m_{\\rm F}m_{\\rm F}})}$$,\n",
    "where $\\Omega_0$ is the full Rabi frequency,  $\\theta$ is the angle between the microwave field and the DC field, S is the transition strength ('Sl' in the code), and $L$ is the cavity line profile.\n",
    "\n",
    "For the $\\sigma$-transitions, the expression is the same, except $\\cos\\theta$ is replaced with $\\sin\\theta$.\n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [],
   "source": [
    "def L(f, g):\n",
    "    \"\"\"Cavity line-width\"\"\"\n",
    "    if g == 0:\n",
    "        return 1 # For simplicity, I set this condition fot the case when there is no cavity\n",
    "    return g**2/(g**2+f**2)\n",
    "\n",
    "def So(W_R, theta, wL):\n",
    "    \"\"\"Coupling strength for sigma transitions with cavity\"\"\"\n",
    "    transition_strength = np.array([0.5*np.sqrt(3),  0.25*np.sqrt(2), 0.25*np.sqrt(6), 0.25*np.sqrt(6), 0.25*np.sqrt(2),  0.5*np.sqrt(3)])\n",
    "    if wL == 0:\n",
    "        return 0 # To avoid double counting in the degenerate case\n",
    "    return transition_strength * W_R * np.sin(theta) * np.sqrt(L(np.array([-3,-1,-1,1,1,3])*wL,gL))\n",
    "\n",
    "def Sl(W_R, theta, wL):\n",
    "    \"\"\"Coupling strength for pi transitions with cavity\"\"\"\n",
    "    transition_strength = np.array([0.5*np.sqrt(3), 1, 0.5*np.sqrt(3)])\n",
    "    return transition_strength * W_R * np.cos(theta) * np.sqrt(L(np.array([-2, 0, 2])*wL,gL))\n",
    "\n",
    "def V(W_R, theta, wL):\n",
    "    \"\"\"Interaction part of the Hamiltonian as a function of the Larmor frequency\"\"\"\n",
    "    V_l = 0*qeye(len(states))\n",
    "    V_o = 0*qeye(len(states))\n",
    "\n",
    "    W_par = W_R * np.cos(theta)\n",
    "    W_ort = W_R * np.sin(theta)\n",
    "\n",
    "    for i in range(3):\n",
    "        V_l += 0.5*Sl(W_R, theta, wL)[i]*(states[i]*states[i+4].dag())+0.5*Sl(W_R, theta, wL)[i]*(states[i]*states[i+4].dag()).dag() # Pi transition\n",
    "        V_o += 0.5*So(W_R, theta, wL)[2*i]*(states[i]*states[i+3].dag())+0.5*So(W_R, theta, wL)[2*i]*(states[i]*states[i+3].dag()).dag() # Sigma- transitions\n",
    "        V_o += 0.5*So(W_R, theta, wL)[2*i+1]*(states[i]*states[i+5].dag())+0.5*So(W_R, theta, wL)[2*i+1]*(states[i]*states[i+5].dag()).dag() # Sigma+ transition\n",
    "    return V_o+V_l"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Total Hamiltonian"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [],
   "source": [
    "def H(d, W_R, theta, wL):\n",
    "    return H_d(d, wL) + V(W_R, theta, wL)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Finding steady state\n",
    "\n",
    "We find the steady-state population by solving the Lindblad equation\n",
    "$$    \\dfrac{d}{dt}\\hat{\\rho}=-\\dfrac{i}{\\hbar}\\left[\\hat{H},\\hat{\\rho}\\right]+\\dfrac{1}{2}\\sum_{n} \\left(2\\hat{L}_n\\hat{\\rho} \\hat{L}_n^{\\dagger}-\\hat{\\rho}\\hat{L}_n^{\\dagger}\\hat{L}_n-\\hat{L}_n^{\\dagger}\\hat{L}_n\\hat{\\rho}\\right).$$"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [],
   "source": [
    "theta = np.pi * 0.5\n",
    "\n",
    "ss_1 = [] # Steady-state population in F=1, mF=-1\n",
    "ss0 = [] # Steady-state population in F=1, mF=0\n",
    "ss1 = [] # Steady-state population in F=1, mF=1\n",
    "for d in delta:\n",
    "    ss = steadystate(H(d, W_R, theta, w_L), c_ops, drop_tol=1e-10) #steady state\n",
    "    ss_1.append(expect(P[0], ss)) # Steady-state populations of psi1\n",
    "    ss0.append(expect(P[1], ss)) # Steady-state populations of psi2\n",
    "    ss1.append(expect(P[2], ss)) # Steady-state populations of psi3\n",
    "\n",
    "# Convert to arrays\n",
    "ss_1, ss0, ss1 = map(np.array, (ss_1, ss0, ss1))\n",
    "rho = ss_1 + ss0 + ss1 # Total population of the F=1 ground state"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [
    {
     "data": {
      "text/plain": "Text(0, 0.5, 'Population')"
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEYCAYAAACtEtpmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5sklEQVR4nO2deZgc1XXof6eX2Uf7aJeQQGxiEWBZZonZjcVisPOcBMJmbH8KCTi2E8cxTkJiO3mJ4xcb7zLG2BAbY0wgYMJqOxhvLBIgkBCLEAKtaJdmNJqeXs77o6p6alo909MzPd19S+f3ffNNddW91ed2V99zz3LvFVXFMAzDMAYjVmsBDMMwjPrHlIVhGIZRElMWhmEYRklMWRiGYRglMWVhGIZhlMSUhWEYhlGSSCsLEblVRLaKyMoK3e/fRWSViKwWka+JiAyx3pkiskdEnvf/bixR/usi0jXU+iISF5HnROSB0Ll/EpGNoToXFNSZLSJdIvKp0LnHReSVUJ3J/vlGEfmJiKwRkadEZE6oztUi8pr/d3Xo/Fy/7Gt+3Qb/vPif3RoReUFEThrKZ1ji8zpdRJ4VkYyIfHCk9zMM40AirSyAHwCLK3EjETkVOA04HjgWeCdwRpFy6wa4xa9V9QT/7/ODvM9CYFyZ9T8OrC5S5yuhOg8WXgMeKlLn8lCdrf65jwC7VHWeX++LvqwTgH8E3gUsAv5RRMb7db7ov//hwC7/HgDnA4f7f0uAbxeRoVzeAj4E3FGBexmGUYRIKwtVfQLYGT4nIoeJyMMislxEfi0iRw31dkAT0AA0Akng7UrKKyJx4EvAp8uoMxO4ELiljDrvB9YCq4ZY5RLgNv/4buAc36p6L/CYqu5U1V3AY8Bi/9rZfln8uu8P3et29XgSGCci03y5rhCRp32r5jv+51ESVV2nqi8AuSG2xzCMMom0shiAm4GPqeo7gE8B3xpKJVX9PfC/wGb/7xFVLTaaH4hTRGSFiDwkIscMUOZ64H5V3VxG/ZvwlEuxjvJ639VzazDiF5FW4G+Bzw0gw/f9zvofQm62GcB6AFXNAHuAieHzPhv8cxOB3X7Z8HkGqiMiRwN/ApymqicAWeDyAWQ0DKPKJGotQDURkTbgVOCnoXBDo3/tD4Fi7qGNqvpeEZkHHA3M9M8/JiKnq+oTIvJNPBcVwHQRed4//qmq/gvwLHCIqnb5sYP/xnPDhGWbDvwRcGYRGYrWF5GLgK2qulxECut9G/gCnkX0BeA/gA/jKYmv+PcqfJ/LVXWjiLQD/wVcCdwOFIvN6DDOM8i1c4B3AM/4cjUDWwFE5HagWGzjW6o6JGVvGMbIOKiUBZ4ltdsfufZDVe8B7hmk7geAJ1W1C0BEHgJOBp5Q1euCQiKyrvD+qro3dPygiHxLRCap6vZQsROBecAav7NsEZE1qjpvoPp4CupiX4E0AWNE5IeqeoWq5l1kIvJdIAh+vwv4oIj8O15sJCciPar6DVXd6L9Hp4jcgReHuB1v9D8L2CAiCWAsnntvA/2V20zgcWA7nnsp4VsXM4FNfpngXuE6m/CUyG2qekPhB6+qVxWeMwyjuhxUbii/031DRP4I8pk5C4ZY/S3gDBFJiEgSL7g9JDeUiEwNXDoisgjvc99RINv/qOpUVZ2jqnOAbj+gPGB9Vb1BVWf65S8FfqmqV/jlpoVu/wFgpf8+7w69x03A/1XVb/jtmuTXTQIXBXWA+4Eg0+mD/vso8AhwnoiM991c5+G55xTPZRdkJl0N3Be611X+Z38ysMd3u/0CT4kFGVgTROSQoXy+hmGMPpG2LETkx3gj30kisgEvc+dy4Nsi8vd4Qeo7gRVDuN3deEHbF/HcJg+r6s+GKMoHgT8XkQywH7jU71ARkQeBj6rqpuHUH4R/F5ETfFnXAX9Wonwj8IivKOLAz4Hv+te+B/yniKzBsyguBVDVnSLyBeAZv9znVTVIKPhb4E4R+WfgOf8eAA8CFwBrgG7gGv9eL/nfyaMiEgPSwHXAmyXkRkTeCdwLjAfeJyKfU9WB4kKGYQwDsSXKDcMwjFIcVG4owzAMY3hE1g01adIknTNnTq3FMAzDcIbly5dvV9WOYtciqyzmzJnDsmXLai2GYRiGM4jIgDFCc0MZhmEYJTFlYRiGYZTElIVhGIZRElMWhmEYRklMWRiGYRglMWVhGIZhlMSUhWEYhlGSyM6zGC5f+8VrZLI5EEEAERB/VW3v2P8vfefAK9PvOgXXQ/fzyvvn/DrBhXD9/LXwPaWvbrh8TIR4TEjGhUQsRiIuJOMxEjEhEY/lzyfjQlMyTmtjgrbGBA0JGy8Y0SeVybIvlaWrJ0Mqk6U3myOTVdLZHGn/fyaXozejZHNKTr0/VVCUXA7vNaCq5JS+6/7r8HnAv4f3/uFVldRfrb9wpaXw0kvBoR7weuAywYmWxgTXnnHYCD6t4piyKGDpr16nuzdbazGqRkMiRntjgoltDUwd28y0MU0cPqWN+dPHcNLs8TQlh7RZXUl6Mzle2dLJqk172LR7P5v39LC1M0V3b4Z9qSz701m6ezOkMn17OBX+mBIxT9E1JWP+/zjNyTjjWpJMamuko937O3xyG0dMaae1cfQf70w2x8bd+3lrZzfrd+5ne1eKvfvT7NmfZm9Pmt5MjkxOyWS9Tiidy5FTSMY8ZZ5MxGiIx2hICM3JBBNak4xvbWBCSwPTxzVz+JQ2po5pyg8ORpOd+3pZu62LdTu6Wb+zm13dvezqTrO7u5e9PRlS6Sy9mRypTI7ebI5UOktO+wYrMX8QFfxvTsZpaYjT1pigrSnB5PZGpo1tZsa4ZuZPH8ORU9tJxiszWOlJZ1m2bherNu3h9W1dbN7Tw+Y9Pezc10tXT4be7MGxiaIITGprHBVlEdmFBBcuXKgjncGt/sgByI8ovP+ehu83Wgid03B9/xraf0RReD/vPdQv13e/A8oP8J459TujrJLJeaOlTNbrqNL+KCo435POsi+VoSuVoTOVobMnw46uFFv29LBx9362d/UC0NaY4D3zp7Dk9EM5etqYYX2Gy9/cxY+ffosHX9ycV8IxgcntTUwe00hrQ4LWxjjNDQlaknEak7F+uyMFnaSqks55sqfSOVKZLD3pHN29GXZ1p9nemaIzlen33od2tHL2kZO58PhpnDh7PJViR1eKe5/byC9f3srz63cfMLhoaYgztjlJe1OCpmSceEw8C8+3+ESEbC5HbyZHb1ZJ+51vd8pry/50//u1NSY46ZDxnHVkBxceP43J7U0Va8sLG3bz02Ub+O3r21m7bV/+vAiMbU4yvqWBsc1JxjQnaUrEaEjEaEzE/f8xYiL9nsFgtJ1TZX9vln293rPW2ZPh7b09bOtK5Z/h1oY4586fwjWnzeWEWeOGLf93f/0Gv1j9dv576GhvZPo4b+Azsa2B9ibvu2ht8Czq5oY4iZinoD1r27e8QxZ4PNbnDfAUoW/pS/HX4XKI94znyxR4GcKEPROFZfq8ClLwOly/soMIEVmuqguLXqumsvD3VF6Gt/vcRQXXzsTb8+AN/9Q9qvp5/9pi4Kt4S2ffoqr/Vuq9KqEsDla2d6VYsX43j6zawkMvbqGrN8MnzjmCvzxn3pAfzlQmy9/fu5KfLt9AW2OCC4+bxruPmMTxM8YxfVwTiQqNKMP0pLNs2dPDK2938sqWTpa/uYvfv76D3myORXMn8IVLjuXIqe3Dvv/u7l6++PAr/HTZejI55aip7SyaO4Fjp49l9sQWZk1oYXJ744hHy/t7s+zs7mX9zm5e29rFq1s68515Mi5c/q5D+OvzjqC9KTns91i5cQ833reSZ9/aTVMyxqmHTWLR3AkcOaWdOZNamTGueVRclOlsjg279vPixj38/vXtPPDCZjp7Mlx1yiHceNH8IT8X2ZzypUdeYemvXmdsc5ILjpvGecdMYcHMcUxobai43AcL9aQs/gpYCIwZQFl8qsj5OPAq8B68XdaeAS5T1ZcGey9TFpVhT3eaz/1sFfc8t5FPnnsEHz/38JJ1cjllyX8u5+er3+YvzjyM68+eR0tDbTyenT1p/mv5Br76i9fYn87yjctO4tz5U8q+zxvb93H1rU+zcfd+/nTRbK485RCOmDJ8xTMc1mzt4nu/Wcudz6znsI42bvvwImaMay77Pj9bsYm/vmsFY5qTfOzseXzgpBmMGYHiGQldqQxffvRVbv3tG1xywnRu+pMThjQg+Zf/eYnv/voNLls0m89ecNSIFKfRx2DKomrRTRGZCVwI3FJm1UXAGlVdq6q9eJsVXVJp+YzijG1J8h9/vID/c9JMvvLzV3l+/e6Sde54+i1+vvpt/uGi+Xx68VE1UxQA7U1JPnTaXB795BkcOaWdv7jjWV7csKese+za18uHvv80XakMd/3ZKXzh/cdWXVEAzJvcxr/+4fHc8dGT2bq3h2u+/zSdPemy7vHU2h188ifPc8KscTz2ydO5+tQ5NVMU4LnYbnzffP7qPUdw3/ObuO/5wfYA8/jtmu1899dvcNUph/Cvf3icKYoqUc1UmJuATwODRZpOEZEVIvKQiAQ7nc0A1ofKbPDPHYCILBGRZSKybNu2bZWQ2cDzi37+kmOY0NrAVx57ddCyqUyW/3j0FU49bCIfPm1OdQQcAh3tjXz/mkV0tDXyiZ8852W8DZF/e+hlNu3ezy1XL+Qdh1Qu9jFcTjlsIkuveAevbe3ipp+/NuR6qUyWv7prBbMntHDLhxYyvo7cNdedNY8TZo3jiw+/TDY3sLdDVfnyY68ybWwTf3fh0VWU0KiKshCRi4Ctqrp8kGLPAoeo6gLg68B/B9WLlC36NKnqzaq6UFUXdnQUXZLdGCatjQk++u65/OrVbazd1jVguZ+/tJVd3Wn+7IzDqpLBUw4TWhv4h4vm8/q2fdzz3MYh1XllSyd3LV/PVafM4aQKBslHyqnzJvHH75jF7b9fx1s7uodU546n3mLj7v18/pJja2pNFCMeE5acfiib9/TwmzXbByy3cuNelr+5i7848zAaE5XJ1DOGRrUsi9OAi0VkHZ4b6WwR+WG4gKruVdUu//hBICkik/AsiVmhojOB0raqUXEuOm46AL96dWCr7e7l65k2tok/mDepWmKVxXuPmcIx08fwg9+uG1L5O595i2QsxsfOnje6gg2DT77nCLI55e7l60uWVVV+8Lt1LJozgdPmTayCdOVzztGTGd+S5O7lGwYs8/grWwE4/7hp1RLL8KmKslDVG1R1pqrOAS4FfqmqV4TLiMhU8YeiIrLIl20HXkD7cBGZKyINfv37qyG30Z/ZE1uYO6mVJwZQFrmc8sy6XZxz9GTisfqyKgJEhPefMIOXNu9l/c7BR+TZnPLAC5s566gOxrXUj8smYOrYJk45bCL3r9hEqUSVl7d08uaObt5/4oy6s/gCGhNxzjxyMk+/sWPAMk+8to3jZoxlUltjFSUzoMbLfYjItSJyrf/yg8BKEVkBfA24VD0ywPXAI8Bq4C5VXVUbiY0zjujg92t30Js50Oe/bsc+ulIZjp8xrvqClcF7j5kKwCOrtgxabtm6nWzrTHHxgqIhsrrg4gXTWbejm1Wb9g5a7uGVWxCB844pPxOsmhw7Yyxv702xtbPngGtdqQzPvrWb04+oT6s16lRdWajq40F6rKouVdWl/vE3VPUYVV2gqier6u9CdR5U1SNU9TBV/Zdqy2z0cfzMsfSkc6zfdeCo/MWNXpbRsTPGVlusspg9sYUjprTxxGsD+8YBnn1rNwCnHlafbhuAUw/zOs7nSmSp/WbNdk6cNa7uR+TH+c/Oyo0HZqy9vrWLbE5ZMHNclaUywBYSNMrkkIktAEWDqis37qEhEePwKW3VFqtsjp0xlpc3Dz4aX7lxDzPHN9dV1lAhM8c3M64lycpB0oFzOeWVLZ35jrieOWb6GETgxQ0Hfjdv+W7DQya2VlssA1MWRpnMmuApizd37Dvg2kub93J0Bdf7GU2OnjqGrZ0pdnSlBizz4sY9HD+zvjtYEeG4GWPzVl0xNu7eT1cqw1HDXLKlmrQ2Jjh0UiurNh3YnkBZzJpQ/kREY+TU/6/aqCs62hppaYjzZpHg8OY9Pcwc31IDqcrnqGnepLpXtnQWvb6nO81bO7vr3qUGnuvm1bc76UkXXwBztW9BHTWCpU6qyYzxLby998CYxZs79tHR3ljTSZ4HM6YsjLIQEWZPaCnqhtremWJSW/26bMIcNdUbZa8eQFm8vMXrYI+ZXv/K4pjpY8nklDVbi89/Wb25ExFqMut8OExqa8gvZhnmzR3dHDLBjcFIFDFlYZTN7AkteZdAQCqTZW9Ppu4DqAEd7Y1MbG3g1QGUxRZ/ZDuctZeqzYzxnozFRuMAr27tZNb4lqos2V4JOtoa/dVp+6cDv7Wzm9kTTVnUClMWRtkEyiL8Y965zxsJTnREWYA3T6FYiib0dbxTxtR/ewIZ395bPP6ydW8P08ZWblnz0WZSWyO9mRxdoSXnU5ksW/b2MNssi5physIom472RlKZXL99HLZ3esrCFTcUeJ1SMXcHeB1vsHFPvTOprRGRgS2L7V29TGqvf6UXMKnde4bC383u7jSq3rNn1AZTFkbZtDV5HWh45LfdzypyqlNqa8zLXcjbe3uYUqUd6kZKMh5jYmvjgFbS9s4UHQ5ZfBNbPVnD301nj/esuaC8o4opC6Nsgh9s8AMG2Ob/sF3qlCa1N7Cjq7foUhlb96aY7JDimzKmsagbqiedpTOVcc7iA0/JBQQDk/YmUxa1wpSFUTbBDza8l0IwCpzoUKfU0dZIbzbH3v2ZA6693elZFq4wZUxTUTdU3uJzTIlDoWXhPWttjfW1Wu7BhCkLo2yCzWb6uaE6e2lpiDuVAx90oNsKXFGq6ruh3OlgB7IsAr+/S8piQksDIrAtFLPo6jHLotaYsjDKJnBDdYXcUDv2pZzqkCDk7ihQFnt7MvSkc05ZFpPbm9ixL0W6YFOnwJXjUiwpEY8xoaWh3+z6zpTFLGqNKQujbPIxi5BlsXNfb12voVSMYu4O8FJNwa3MmyljmlCFbZ3929LnhnLruxnf2pBPxwazLOoBUxZG2QQ/2LBl0ZPO0px063EqFkgF2LPf84+Pr8M9LAZiXIvnGtxbsCe3izELgOZkvN/yJYHL05WJhVHErV+3URcUy4ZKZXI0Jd3a5nJ8SwMx4YC5FsH8kZYGd9rT7MsanvsCXtvaGhPOfTdNyRip0J4pnT1pmpNxJxapjCr2yRtlk4jHaE7G6Ur1jWJ70lkaE249TvGYMLY5ya7u4sqi2SFl0eIrg/0FymJ3dy/jW93LIGpMHGhZtJkLqqZU9dctInEReU5EHihy7XIRecH/+52ILAhdWyciL4rI8yKyrJoyG8Vpa0oULMfgnmUB0NKQoCfdPyjc3eu7PBzK7ArcM/tS/dOA96eztCTdaUfAgZZFhnZzQdWUan/6H8fbGrXYwvpvAGeo6i4ROR+4GXhX6PpZqjr41mZG1WhvTPRzQ/WkszQl3FMWjcnYAUt7u+yG2l/Qlv3pHE2OxZIAGovELMyyqC1Ve4pEZCZwIXBLseuq+jtV3eW/fBKYWS3ZjPJpb0ocELNodLBTKgykQp8rxyk31AAxi5501kmLrzER62fxdfZkLBOqxlTz130T8GkgV6IcwEeAh0KvFXhURJaLyJKBKonIEhFZJiLLtm3bNiJhjcEpdEO52ik1J+MHjMb7LAt3OqfA1VRMWbik9AKakvF+bqiunozNsagxVVEWInIRsFVVlw+h7Fl4yuJvQ6dPU9WTgPOB60Tk9GJ1VfVmVV2oqgs7OjoqIboxAG2NiXzqrKrSk845F+AGr1M6UFlkaEzEiMfqfxHBgHw2VGHMotdR92AiRqrQDWVLfdSUav26TwMuFpF1wJ3A2SLyw8JCInI8npvqElXdEZxX1U3+/63AvcCiaghtDExbYzJvWfT6s4ZdtCyakvEiAe6sU/EKgIZEjGRc6C5QfD2ZaFgWnT1pc0PVmKooC1W9QVVnquoc4FLgl6p6RbiMiMwG7gGuVNVXQ+dbRaQ9OAbOA1ZWQ25jYNqbEvkJYEFn66Jl0dxwYMzCUxbudUzNyfgBqbP7e93MUmtKxOnN5sjmFFWlK2Uxi1pT009fRK4FUNWlwI3AROBb/h4CGVVdCEwB7vXPJYA7VPXh2khsBIQ72VTG+9/oZKdULBsq45xlAV76bGHqbCqddTQbypM5lckSEyGnblquUaLqykJVHwce94+Xhs5/FPhokfJrgQWF543akowJ6aw36kv5lkWTo5ZFsQC3i8qiuSF+gBtqfzpLs4OdbPAspdI5kv5xMu5ODCmKmF1nDIuEv+xCNqdOWxbFXTdu+vlbGvq3JZ3Nkcmpk8oieJZ6MlmCzQoTMfcGI1HCPn1jWARr9KSzmo9ZuGhZNPqB1Fyub7e8fb0Zp2ZvB7Q09HdDBe41F903geusJ53LJ1AkHXy+ooR9+sawCFwC6Vwu3ym5alkA/TJvnLYsQm6o4LjJwbY0JoLvJUsm6ynypEOpzFHElIUxLBL+DzeT1XxH66JlESyrHu5kXY1ZtDTE+03K6+n1vhcX3VBhyyJQFglbcbam2KdvDIvgh5vJ5hx3d/i+8ZCy2NebcTR1NtFvUl5PJvhe3PuZBxMJU+ks6ZzvhrIAd01x7yky6oIGX1n0ZnN5y8LJtaGKLMC331HLorWxfzZUfo0rB5V48Cz1ZHL5rWJtL4vaYp++MSwS8T43VN6ycHBZiaaCfSB6M14GkYvKornADRUoQCeVRaLP4su7oSxmUVNMWRjDIu+GyrltWTQl+wKp0LeXRbODbqiWZMJTdv5I3OXEg6b8pDyzLOoF+/SNYdHgWxa9Gbcti+a8ZeF1SMHIvNVBy6K10V9M0P8+eiJiWaSDbChTFjXFPn1jWAQTpDK5XN88Cwc7peaCAHfedeOgsih0qQXfi8ttSYUspYQFuGuKKQtjWAQ/3HQ2NIPbwdTZpoLU2YzDo9iGvGvQa0N+noWD7sH82lDpLOlc8J2Ysqgl7j1FRl3QkJ/B7VkWDfEYMQcDkPnReKAs/DRNl/ayCAhkzvoKz+VsqKawGypjMYt6wD59Y1j0zbPwLAsXrQroc9EEG+1kc+5m3uQz1HyF1zfPwj1lkYwLIr4bym+PrQ1VW+zTN4ZFot9yHzknM26gmGXhKQunLQu/DT293iJ8LipyEaEpES8IcLv3nUQJ954ioy7Iu6EyOVIZN/dMgL4lSoJsqL6cfvfaE1hDQee6P+1tqSriZifblIxZ6mwdYZ++MSz6XB7efhYujl7Bc6cl45J32bgds+hbNh7wLT43vxfw0mf7Tcozy6KmuPskGTUlGHmns97Iz+VRXyIWy3ewWYczbwpjFpmcOv29JBNCJquhtaHcbUsUqOqnLyJxEXlORB4ock1E5GsiskZEXhCRk0LXFovIK/61z1RTZqM4DaH9LDI5dXrUF49JfvTqcswiURCzyOZyxB11QYGnxNM5tWyoOqHan/7HgdUDXDsfONz/WwJ8GzwFA3zTvz4fuExE5o++qMZg9K0N5a2l5KKPPyAeE3Lqd7AOxywCBZfJKws3lV5AIiZkc7l8e1wekESBqv0iRGQmcCFwywBFLgFuV48ngXEiMg1YBKxR1bWq2gvc6Zc1akhfNpSSzeWcTDUNiMck5LpxN2aRn1WfDVkWDrYjILD48tlQDirwKFHNT/8m4NNAboDrM4D1odcb/HMDnT8AEVkiIstEZNm2bdtGLLAxMOFsqExWne+U/IQbp0exfZaF15isujlfJCARF7I5DWVDuduWKFAVZSEiFwFbVXX5YMWKnNNBzh94UvVmVV2oqgs7OjqGIakxVMKrzjofsxDP3QFuT8oLOtNwzMLFWfUBcT9mEawN5fKAJApUax3m04CLReQCoAkYIyI/VNUrQmU2ALNCr2cCm4CGAc4bNSSc0x+FmEVgUbg8z+LAmIU6qfQCgphFOqf+jG532xIFqvKLUNUbVHWmqs4BLgV+WaAoAO4HrvKzok4G9qjqZuAZ4HARmSsiDX79+6shtzEwydDaUK7HLBJxIVeQOht30FJKFMyzyOaUmMMdbCKIWWTcTs2OCjXd4UVErgVQ1aXAg8AFwBqgG7jGv5YRkeuBR4A4cKuqrqqNxEZAPCbExBuJOx+zkJBl4bAbqtCycN09mIgLqXSQbeduO6JC1ZWFqj4OPO4fLw2dV+C6Aeo8iKdMjDoiEY+RjkLMIpQ663Y2VF86M3iWhYvtCIjHYmRyWecnfUYF+waMYdMQj5HOqO8bd/dR6jcpz+H9novFLNyelOelNJuyqA/sGzCGTSIufjaU2zGLfpPy8qmz7v00gg41HLNw2bIIYhaZrNuWa1Rw7xdh1A2JWIx0Vsm6HrOIRTNm4byyCOZZOL7GVVSwb8AYNg1x8RYSjEDMIjw3ITjnGvm1oYKYhbqtLOL+Ao9eNpS77YgKpiyMYZOIx8hkc87HLBIhZREsLeGirz8ej5hlERM/gSLn9PMVFewbMIZNIi75GbYud0qxUOqsNzcBJ2c+J4q4oVx0pwXEY0LWXxvKLIvaY8rCGDZeNlTO+U4pPCnPSwN282cRtUl5ybinxC0bqj6wb8AYNongx5xTJ2c8B/S3LNzN7OqbZxGNSXlBLMmyoeoDUxbGsPGyoSJgWfSblOeunz8WE0T6gvS5nOa3WnWR4PlK58yyqAfKnsEtIpOBtvA5VV1bMYkMZ2iIh5WFuz/m8KS8KCi+cBqwywPysGVhyqL2DFlZiMhi4HvAVPovG654azYZBxmJuNDdm/WOHe5g4wXZUC6PxsNzRrKuWxahmIXLz1dUKOdJ+ibwBaBNVWOhP1MUBymJeIyetKcsXI5ZJGIxsup+zAK8toStJJcH5EFKswW464Ny3FDjge/4C/4ZBg1xySsLlzvYWMiycD0o7M16Dk/Kc7eT9RYStNTZeqGcJ+l7+MuGGwZ4o9j9eWXhbqcUnpQXpZhFFCwLgFQm62w6c5Qox7I4GfhLEfkMsCV8QVVPr6hUhhMk4kJPOpc/dpWY9LcsXM2GgsKlS9xOPAieqf29WbMs6oBylMUt/p9hAN4qp4Fl4XIH28+yyDrewfquG3B/Ul5gWfRkbLmPemDIykJVbxvum4hIE/AE0Oi/592q+o8FZf4GuDwk19FAh6ruFJF1QCeQBTKqunC4shiVIxkXejO+ZeGwsoj1Szd1e+kSLw3Y+04yuZzTFl8Qb+m1bVXrgrLmWYjINcCVwAxgI/Cfqvr9IVRNAWerapeIJIHfiMhDqvpkUEBVvwR8yX+f9wGfVNWdoXucparby5HXGF3CfmSXR36eZRF0sG4HU4N0U4Bczn2LL8Dl7yQqlDPP4u+Aq4D/AN4EDgE+LSLTVfVfBqvrZ1B1+S+T/t9gWVWXAT8eqmxGbWgIKwuHf8yFfn7XO9hs2Epy2A0V/h4aEu4ORqJCOZbFR4EzVfXN4ISIPILnXhpUWfhl48ByYB7wTVV9aoByLcBi4PrQaQUeFRHFS9+9eYC6S4AlALNnzx5Km4wR0JTsm2LjcgcbVhYZx2MWQbqpqpJTt7+XsDURftaM2lDOr6IV2FZwbgfQPJTKqppV1ROAmcAiETl2gKLvA35b4II6TVVPAs4HrhORotlXqnqzqi5U1YUdHR1DEcsYAY2JsBvK3U4pEZPQpLxoWBaB8nO5LeE5Io1mWdSccr6Bh4EficiRItIsIkcBtwGPlPOGqrobeBzPeijGpRS4oFR1k/9/K3AvsKic9zRGh8ZkNGIW4Ul5aeeDwt7uhYHyc1lZhAcgpixqTznfwPV4GUkr8OIPzwP7gI+VqigiHSIyzj9uBs4FXi5SbixwBnBf6FyriLQHx8B5wMoy5DZGicZEyA3lcAdbOCnP9Q42OpZFWFmYG6rWlJM6uxe4SkQ+BEwCtqtqbojVpwG3+XGLGHCXqj4gItf6917ql/sA8Kiq7gvVnQLcK16gLgHcoaoPD1VuY/RoSkbDDRUTIafekt6uxyyCbKggI8rl7yUcswhbsUZtGFRZiMgcVV3nHx9acLnN78BLLlGuqi8AJxY5v7Tg9Q+AHxScWwssGOz+Rm3oZ1k43CkFHWpWNQLLfXgTJXORsCzCMQuzLGpNKcviRaDdP16Dl5VU+PTZEuUHKWE/ssuTpgIXWjanXrqpwy61YInyTASURb+YhVkWNWdQZaGq7aFj+7aMfoSVhcudUjAXIfD1u21ZeBMMo2FZWIC7nhjyNyAiXxvg/E0Vk8ZwisZQ7rvLHWw85IbyNj9yuy2ZbMiycHhSXjgrzdxQtaccdf2hAc5fWQE5DAdpioplESiLrGdZJB0OcCfj3qS8KGRDhRMNmswNVXNKZkOJyIeDsqHjgEMBW6/pICVsWbgcswgHuDM5dT5mYamzxmgwlNTZwHJooL8VocDbwNWVFspwg6jELGKxcMzC9W1VhUzOJuUZlaekslDVswBE5J9V9e9HXyTDFaK03AeQzyJyuYONxyTvTgteu0rcsqHqinIm5eUVhXgTLCR0baiT84wI0RiZhQS9jqhvUp67bUnEhbTfDnBbiYddm+aGqj3lZENNF5F7RWQHkAHSoT/jIKQpKvMsfNGDwHDc4QB3ELPI5d1QbrclwALctaecb+A7QC9wDt7aUCcB9wPXjoJchgNEzbIIJuW5vNFOIhYjk82FJuXVWKARELaKGlxuSEQoZz+LU4HZqrpPRFRVV4jIR4DfAd8dHfGMeiYqMYtgLkIml3N+D4gDFxJ0t5MNfw/i8HyRqFDOk5TFcz8B7BaRDrxVZ2dUXCrDCcIKwuUONpA9CvuJx/2FBLMRmJTnsmszipTzbTwFXOAfPwL8BLgHWFZpoQw3CI/2XP5hB8oh5SsLl0fjUV2i3Kg95bihrqRPuXwC+BTQBtxUWZEMF3H5hx3InkpHwLKIRWkGt7uyR5FyUmd3h473A18YDYEMN3HZ3ZFXFpksgNM75SUDl1rWa4vLysLlmfRRpNR+Fp8fyk1U9cbKiGO4SszlTqnADeXyiDboYAMryWVl4fL3EEVKWRazqiKFYdSQwgC3y4qvMP7icofr8o6FUaTUfhbXVOJNRKQJeAJo9N/zblX9x4IyZ+Ltvf2Gf+oeVf28f20x8FW8TZZuUdV/q4RchgHRyoaKSf+2mGVhVIohxyyKbKuap9S2qkAKOFtVu0QkCfxGRB5S1ScLyv1aVS8qeN848E3gPcAG4BkRuV9VXxqq7IYxGHllkfU6WJdz+vPKIuu+snDZwosi5WRDFdtWVf3/gy7coqqKN+sbIOn/6cA1+rEIWBMoJBG5E7gEMGVhVIREgWURhWB9OgLKwqgvhuwUVNWYqsb9/zFgOnAzQ9z8SETiIvI8sBV4TFWfKlLsFBFZISIPicgx/rkZwPpQmQ0MMBFQRJaIyDIRWbZt27Yhtsw42AlG41HoYGOFysJhxWfUF8OOIKnqFrz5Fv86xPJZVT0BmAksEpFjC4o8CxyiqguArwP/7Z8v9rQXtUpU9WZVXaiqCzs6OoYilmHkU2XTeTdULaUZGYGeS2fdn2dh1BcjTTc4Emgpp4I/X+NxYHHB+b2q2uUfPwgkRWQSniURzsqaCWwavshGJRnbnKy1CCMmnvfzu9/BxiMU4A44dFJrrUUwKC/A/Wv6j+hbgGOAknMx/HWk0qq6W0SagXOBLxaUmQq8raoqIovwFNkOYDdwuIjMBTYClwJ/OlS5jdHlV39zJl2pTOmCdUyhnz/msGlR6FJzPaPoib85i3Gt7g9IokA5Ae5bCl7vA1ao6mtDqDsNuM3PbIoBd6nqAyJyLYCqLgU+CPy5iGSA/cClfmA8IyLX461HFQduVdVVZchtjCLjWhoY19JQazFGRF5ZZCKgLAqD9Y4ri9kTy3JcGKNIOct93DbcN1HVF4ATi5xfGjr+BvCNAeo/CDw43Pc3jMGIUgZRsJ5jFNpi1Bfl7JTXICKfF5HXRGSf//8L/oQ7w3CWYKZwb94NVUtpRkYsQvEXo74oxw31bbyA9l8CbwKHADfgpbF+uPKiGUZ1CFaV6M2o/9rdDjZKacBGfVGOsng/cFho9dmXROQpvMl6piwMZwksiygGuE1ZGJWinNTZLRyYJtsMbK6cOIZRfeKFHazDyuKAmIXDbTHqi3Isi/8EHhaRr9M39+E64HYROTsopKq/rKyIhjG6xAsm5bm82GmUFhI06otylMWf+f8/W3D+Wv8PvHkYAy44aBj1SOGkvCi4oXqzSkzcXhTRqC/KSZ2dO5qCGEatKJxn4fJoPNwW2w/CqCTlWBaISAI4FS8DagPwe1V1e/qucdBTuES5w7oiv65VOptz2p1m1B/lLPdxFPAzvKD2eryYRY+IvE9VV4+SfIYx6gTKIROBbKi8ZZFTC24bFaWcsce38JYkn6Wqp6jqTGCpf94wnEVEiEk0JrKFFxJ0eb6IUX+UoyxOAL7sr9cUcJN/3jCcJhGLRWKehYTSgF1fRNCoL8pRFpuAMwrOvRtbLtyIALFYNCblBVZRJptz2kIy6o9yAtw3APeLyAP0LfdxIXDFaAhmGNUkLkImSJ11ODAc3vzIZaVn1B8llYWItAB/DxwL3AesxttSdSVwo6q+OqoSGkYViMUknw3lcmA4vK2qWRZGJRmKZfEN4J3AQ8AFwC5V/YtRlcowqkw8Jn3ZUA53suGlS5LxsjLjDWNQhmJwnw+cp6qf9o8vHF2RDKP6xEXy+1a77L7pW0hQzbIwKspQlEWrqm4GUNX1wNhy30REmkTkaRFZISKrRORzRcpcLiIv+H+/E5EFoWvrRORFEXleRJaV+/6GUYpYTCKx+F4stJCgZUMZlWQodmpCRM4CZIDXQ1k8MAWcrapdIpIEfiMiD6nqk6EybwBnqOouETkfb07Hu0LXz1LV7UOQ1zDKJhGKWTisK0JrQ9k8C6OyDEVZbAVuDb3eUfC65OKB/tyMLv9l0v/TgjK/C718Epg5BNkMoyLERAhmELnsvglkV3XbQjLqj5LKQlXnVOKNRCQOLAfmAd9U1acGKf4RvIB6XgzgURFR4DuqenMlZDKMgLCCcFlZhOMtZlkYlaRq6RKqmgVOEJFxwL0icqyqriws57u4PgL8Qej0aaq6SUQmA4+JyMuq+kSRukuAJQCzZ88ejWYYESWsIFwekIf1Q9zh+SJG/VH1x8nflvVxYHHhNRE5HrgFuERVd4TqbPL/bwXuBRYNcO+bVXWhqi7s6OiovPBGZOnXyTqsLfpZSA63w6g/qqIsRKTDtygQkWbgXODlgjKzgXuAK8MT/USkVUTag2PgPLwJgYZRMcKdbBRSZ8HcUEZlqZYbahpwmx+3iAF3qeoDInItgKouBW4EJgLf8hdDy6jqQmAKntsqkPcOVX24SnIbBwlR6WTDslvqrFFJqqIsVPUF4MQi55eGjj8KfLRImbXAgsLzhlFJEvGwZVFDQUZIWHaXLSSj/rAQmGHQ37/vcjZUVNph1B+mLAyD/u4bl0fksYikABv1hykLw6D/iNxpZRGRdhj1hykLwyA6I3JzQxmjhSkLw6DQsqihICNEQr9osyyMSmLKwjDoG4WL9O1j7SJhpWeps0YlMWVhGPS5oVyf9RyVNa6M+sOUhWHQNwp33XUTFt/lyYVG/WHKwjDoUxIxx38R/QLcpiuMCuL4T8MwKkOwQqvrlkVUli0x6g9TFoZBn3/f9ZhFzFadNUYJUxaGQdgN5X4HGyi+hPmhjApiysIw6OtgI6Ar8m1w3aVm1BemLAyDPpdNFNJNYxFqi1E/mLIwDMKT8tzvYPMutQi0xagfTFkYBtEJcEOoLWZZGBXElIVhEJrBHYEONmhCFNpi1A/V2oO7SUSeFpEVIrJKRD5XpIyIyNdEZI2IvCAiJ4WuLRaRV/xrn6mGzMbBRWBRRMCwyCs+c0MZlaRalkUKOFtVFwAnAItF5OSCMucDh/t/S4BvA/j7dn/Tvz4fuExE5ldJbuMgIUqum0Dx2UKCRiWpirJQjy7/ZdL/04JilwC3+2WfBMaJyDRgEbBGVdeqai9wp1/WMCpGlILCEqE5I0b9ULWYhYjEReR5YCvwmKo+VVBkBrA+9HqDf26g88XeY4mILBORZdu2bauY7Eb06Vvuo7ZyVIKgLVEI1hv1Q9WUhapmVfUEYCawSESOLShS7MnWQc4Xe4+bVXWhqi7s6OgYkbzGwUXcX0EwSm6ouKWvGBWk6o+Tqu4GHgcWF1zaAMwKvZ4JbBrkvGFUjKgsJAjmhjJGh2plQ3WIyDj/uBk4F3i5oNj9wFV+VtTJwB5V3Qw8AxwuInNFpAG41C9rGBUjHqGYRZTmjBj1Q6JK7zMNuM3PbIoBd6nqAyJyLYCqLgUeBC4A1gDdwDX+tYyIXA88AsSBW1V1VZXkNg4S8ummEXDd2DwLYzSoirJQ1ReAE4ucXxo6VuC6Aeo/iKdMDGNUyPv5IzAaj9IEQ6N+iMA4yjBGTp9l4X4HG6VFEY36wZSFYRBeotz9DjZKc0aM+sGUhWFgbijDKIUpC8MgvER5jQWpAPkAdxQaY9QNpiwMg4itDRWh+ItRP5iyMAyitVKr2EKCxihgysIwCE3Ki0AHGw/24I5AW4z6wZSFYRBefK+2clQCm8FtjAamLAyDaKWbii0kaIwC9jgZBtEKCkdpnSujfjBlYRiEJ+XVWJAKEKxvFYXMLqN+MGVhGEQrdTYWoWC9UT+YsjAMouW6CRSepc4alcSUhWEQrXkWsQgtXWLUD6YsDINordRqbihjNDBlYRhEdG0oUxZGBanK5kciMgu4HZgK5ICbVfWrBWX+Brg8JNfRQIeq7hSRdUAnkAUyqrqwGnIbBw+xCE1ki9Jy60b9UK1tVTPAX6vqsyLSDiwXkcdU9aWggKp+CfgSgIi8D/ikqu4M3eMsVd1eJXmNg4xElLKhItQWo36oihtKVTer6rP+cSewGpgxSJXLgB9XQzbDgL5RuERgNG4BbmM0qHrMQkTm4O3H/dQA11uAxcB/hU4r8KiILBeRJaMupHHQ0TfPosaCVIC4xSyMUaBabigARKQNTwl8QlX3DlDsfcBvC1xQp6nqJhGZDDwmIi+r6hNF7r8EWAIwe/bsCktvRJlASUTBzx+LUGaXUT9UbRwlIkk8RfEjVb1nkKKXUuCCUtVN/v+twL3AomIVVfVmVV2oqgs7OjoqI7hxUBClhQRjEbKSjPqhKo+TeI7g7wGrVfXLg5QbC5wB3Bc61+oHxRGRVuA8YOXoSmwcbERpuY8ozUY36odquaFOA64EXhSR5/1znwVmA6jqUv/cB4BHVXVfqO4U4F4/8JgA7lDVh6shtHHw0GdZ1FiQCmALCRqjQVWUhar+Bij55KrqD4AfFJxbCywYFcEMwydKS5RHyaVm1A/m1TQM+uZZRKGDDdqQiMK2f0bdYMrCMIjmDO4otMWoH0xZGAahoHCU3FARaItRP5iyMAwitlNeMCnPLAujgpiyMAyi6YYyy8KoJKYsDINouaHEZnAbo4ApC8Ogb25CFLKhgpnbUbCSjPrBlIVhAAlfW0RhiYwo7fpn1A8R+GkYxsiJ0oZBtp+FMRpUddVZw6hXxjQl+NR5R3De/Km1FmXEXHT8NFoa4qYsjIpiysIw8ILC1599eK3FqAjzJrczb3J7rcUwIoa5oQzDMIySmLIwDMMwSmLKwjAMwyiJKQvDMAyjJKYsDMMwjJKYsjAMwzBKYsrCMAzDKIkpC8MwDKMkoqq1lmFUEJFtwJvDrD4J2F5BcWpJVNoSlXaAtaUeiUo7YGRtOURVO4pdiKyyGAkiskxVF9ZajkoQlbZEpR1gbalHotIOGL22mBvKMAzDKIkpC8MwDKMkpiyKc3OtBaggUWlLVNoB1pZ6JCrtgFFqi8UsDMMwjJKYZWEYhmGUxJSFYRiGURJTFoMgIh8TkVdEZJWI/Hut5RkJIvIpEVERmVRrWYaLiHxJRF4WkRdE5F4RGVdrmcpBRBb7z9MaEflMreUZLiIyS0T+V0RW+7+Nj9dappEiInEReU5EHqi1LCNBRMaJyN3+72S1iJxSqXubshgAETkLuAQ4XlWPAf5fjUUaNiIyC3gP8FatZRkhjwHHqurxwKvADTWWZ8iISBz4JnA+MB+4TETm11aqYZMB/lpVjwZOBq5zuC0BHwdW11qICvBV4GFVPQpYQAXbZMpiYP4c+DdVTQGo6tYayzMSvgJ8GnA6m0FVH1XVjP/ySWBmLeUpk0XAGlVdq6q9wJ14gxHnUNXNqvqsf9yJ1yHNqK1Uw0dEZgIXArfUWpaRICJjgNOB7wGoaq+q7q7U/U1ZDMwRwLtF5CkR+ZWIvLPWAg0HEbkY2KiqK2otS4X5MPBQrYUogxnA+tDrDTjcwQaIyBzgROCpGosyEm7CG0zlaizHSDkU2AZ833ep3SIirZW6eaJSN3IREfk5MLXIpb/D+2zG45nZ7wTuEpFDtQ5zjUu047PAedWVaPgM1hZVvc8v83d4rpAfVVO2ESJFztXds1QOItIG/BfwCVXdW2t5hoOIXARsVdXlInJmjcUZKQngJOBjqvqUiHwV+AzwD5W6+UGLqp470DUR+XPgHl85PC0iObwFurZVS76hMlA7ROQ4YC6wQkTAc9s8KyKLVHVLFUUcMoN9JwAicjVwEXBOPSruQdgAzAq9nglsqpEsI0ZEkniK4keqek+t5RkBpwEXi8gFQBMwRkR+qKpX1Fiu4bAB2KCqgZV3N56yqAjmhhqY/wbOBhCRI4AGHFuVUlVfVNXJqjpHVefgPUwn1auiKIWILAb+FrhYVbtrLU+ZPAMcLiJzRaQBuBS4v8YyDQvxRh7fA1ar6pdrLc9IUNUbVHWm//u4FPilo4oC/3e9XkSO9E+dA7xUqfsf1JZFCW4FbhWRlUAvcLVjI9ko8g2gEXjMt5SeVNVrayvS0FDVjIhcDzwCxIFbVXVVjcUaLqcBVwIvisjz/rnPquqDtRPJ8PkY8CN/QLIWuKZSN7blPgzDMIySmBvKMAzDKIkpC8MwDKMkpiwMwzCMkpiyMAzDMEpiysIwDMMoiSkLwzAMoySmLAzDMIySmLIwjBEgIo+LyC4RaSyjznQR2TCachlGpTFlYRjDxF9x9Xi8JbovLqPqBcDDoyGTYYwWpiwMY/hcBdwL/AC4uox6FwAPAohIu4hsFJEzwgX83ehURCZWSljDGAm23IdhDBMRWQNcCzyLt1fFoar6dok6SWAzMFdVO0Xkn4BTVPW9Rcp2ARep6uOVlt0wysUsC8MYBiLyB0Ar8L+quhP4JfCnQ6h6OrDCVxRxvB0Zb/Hv2SEih4bKZoBm/9o9/urHhlETTFkYxvC4GviJqmb91z9maK6ovAsKOBaYjLcSLcBfAUsARKQZaAeC7XwPA14fudiGMTxMWRhGmfgd+R/jKYiA+4B5IrJARGL+jn9B+cdD5S4A/sc/ngHsCu0yt5g+5XAGsAt4TkSagExIMRlG1TFlYRjl835gJ94OhE1+Z57FsxiuwtsLeS2AiIwDdvvHc4FGVX3Zv89OvJ3Z5orIZXgbbM336/wTcJOq5oCjgKCOYdQE2/zIMMrnamAOsL/ItbeBJ4EX/NcnAC/6xxfS54ICb/e8O4Hn/eOLgZ8BrwF3AF/0yx1LBXc8M4zhYMrCMMpEVRcPdl1EbgRe8V9eAvzWP74Ab7e/4D5ZoHALz/lFbnkM8FSR84ZRNUxZGEblOQ7PnfTHeFuoZoG7gceB/x3G/Y4BzvG3ZQX4E1XdUQlBDWOo2DwLw6gwIvILVT2n1nIYRiWxALdhVBA/2J2qtRyGUWnMsjAMwzBKYpaFYRiGURJTFoZhGEZJTFkYhmEYJTFlYRiGYZTElIVhGIZRElMWhmEYRklMWRiGYRgl+f/rpu1R3vk67QAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(delta/w_L, rho)\n",
    "ax.set_xlabel(r'$\\Delta_{\\mu}/\\omega_{L}$',fontsize=12)\n",
    "ax.set_ylabel('Population',fontsize=12)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Calculate absorption coefficient"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "outputs": [],
   "source": [
    "# Fractions of atoms in the velocity classes resonant with different transitions\n",
    "c = 2.99792458e8 # Speed of light in m/s\n",
    "f = 384.230484468562e12-302.073888e6 # Optical transition frequency in Hz\n",
    "kB = 1.380650324e-23 # Boltzmann constant in  J/K\n",
    "T = 300 # Vapor temperature in K\n",
    "m = 87*1.6605387313e-27 # 87Rb mass in kg\n",
    "v = np.sqrt(2*kB*T/m) # Most probable thermal speed in m/s\n",
    "\n",
    "N0 = 1. # The ensemble fraction in  F'=0\n",
    "N1 = N0 * np.exp(-((72.2e6*c)/(f*v))**2)#ensemble fraction in F'=1\n",
    "N2 = N0 * np.exp(-((72.2e6+156.9e6)*c/(f*v))**2)#ensemble fraction in F'=2\n",
    "\n",
    "# Total sigma+ absorption as a function of sublevel populations\n",
    "def Asp(P_1, P0, P1):\n",
    "    A = P_1*(N0*1./6+N1*5./24+N2*1./24)+P0*(N0*0+N1*5./24+N2/8)+P1*N2/4\n",
    "    return A\n",
    "\n",
    "# Same for pi\n",
    "def Api(P_1, P0, P1):\n",
    "    A = P_1*(N0*0+N1*5./24+N2*1./8)+P0*(N0/6+N1*0+N2/6)+P1*(N0*0+N1*5./24+N2*1./8)\n",
    "    return A\n",
    "\n",
    "# Same for sigma-\n",
    "def Asm(P_1, P0, P1):\n",
    "    A=P1*(N0*1./6+N1*5./24+N2*1./24)+P0*(N0*0+N1*5./24+N2/8)+P_1*N2/4\n",
    "    return A\n",
    "\n",
    "# Total absorption for light split between polarization components as:  π : 0.5,  σ+ :0.25,  σ− :0.25\n",
    "def A(P_1, P0, P1):\n",
    "    return 0.5*Api(P_1,P0,P1)+0.25*Asp(P_1,P0,P1)+0.25*Asm(P_1,P0,P1)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "outputs": [
    {
     "data": {
      "text/plain": "Text(0.5, 1.0, '$0.25\\\\sigma^{+}$')"
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEeCAYAAABsaamyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0T0lEQVR4nO2deZhkVXn/P28tvc10z9qzMM3MALKPMMJIUJSgRGRz1GASjRiiJkZjjFnckMSo0V9MTKIkrgRcIu5GiKAIKAxuYZkBZlhmWAaGWWCYnn2m16p7398f997qW9XVPdXd1VV1br2f56mnb92t3tvVfb7nXc45oqoYhmEYzUmq3gYYhmEY9cNEwDAMo4kxETAMw2hiTAQMwzCaGBMBwzCMJsZEwDAMo4kxETCMBkJE1tTbBqO5MBEwDMNoYkwEjKZGROaKyA0i0iciz4jIH5Y5p1VErguPHxKRB0TkopJz1ojIoIgcDl+P1e4pDGPymAgYzc7ngWFgIfBm4IsicmrJORlgG/DbwCzg74HvicjykvP+QlVnhq8TKzVARJaGIrIGWBlti8jSyT2SYVSO2LQRRrMiIjOAfcAKVX083PcNYIeqfugI124APqaq/xO+XwNcr6rXjnF+FvgogdAsIRAWgA2qenrsvDWqet4UHsswJoR5AkYzcwLgRQIQsh4o9QSKEJGF4bWPlBz6JxHZLSK/FpHzSo59AjgfeDkwG/g5cAPw+skabxjVwETAaGZmAgdK9h0AOse6IOzRfxP4uqpuih36IHAsQS//GuAmETkuvKYT+EvgLaq6TVX7gP8B5qrqU/H7mxdg1BonRUBEviIiu0Tk4Srd759F5OHw9QcTuG6WiNwkIutF5BEReesY531NRJ4WkQfD18qS4y8WEU9E3hDbN1tEfiAim0Rko4i8JNz/6XDfhjChOTvcv1xEBmKf8aXYvT4pIttE5HDJ57aKyHdF5EkRuSce4xaRK0TkifB1RWz/MeG5T4TXtoT7RUT+I7zXBhE5o9Lf41iIyLkicr+I5OO/mypyGOgq2dcFHBrDnhTwDYIcwl/Ej6nqPap6SFWHVPXrwK+Bi8PD5wJPqeoTsUvmADun/giGMTWcFAHga8CF1biRiFwCnAGsBH4LeL+IlDYMiMiWMpe/G3g0jOmeB/xb1CiW4f2qujJ8PRi7bxr4Z+DWkvOvBn6qqicBpwMbw/23E8SwTwMeB66MXbM59hnvjO2/CTirjE1vB/ap6guAz4R2ICJzgX8g+H2cBfyDiMwJr/ln4DOqejxBPP3t4f6LgOPD1zuAL47xe5gIW4E/Br5VhXuV43EgIyLHx/adzugwDyIiwHUECeTLVDV3hHsrIOF2N8HvKn6v1wM3T950w6gOToqAqv4C2BvfJyLHichPRWSdiPxSRE6q8HanAHepaj5009dTucAo0Bn+U88MbcpXeG3EewhCA7uiHaEInUvQ6KCqw6q6P9y+TVWjz7gb6Dmikap3q+pzZQ69Fvh6uP0D4PzwWV4N3K6qe1V1H4HwXBgee2V4LuG1r4vd67814G5gtogsDp/nchG5N/RQvhwK3xFR1S2qugHwKzl/ooTf9w+Bj4vIDBE5J3yOb5Q5/YvAycBrVHUgfiD02l4tIm0ikhGRNxN8f5GwPwycISIrRaQd+CeCv53vTsdzGcZEcFIExuAa4D2qeibwPuALFV63HrhIRDpEZD7wCuDoCq/9HEHD8CzwEPBeVR2rwfpkGCb5jIi0AojIEoIe4ZdKzj0W6AW+KkFN+rVhJUspbwNuib0/Jjz/LhF5eQX2LyEofSQUlgPAvPj+kO3hvnnA/pgIRfsZ6xoRORn4A+AcVV0JeAQVMo3CnwPtBCL8beBdqvoIgIjcIiIfFpFlwJ8ReIs7ZWQsQPQcWYLEby+wm0DYX6eqjwGo6lrgk8BPgKeARcDFFXgThjHtZI58SuMjIjOBlwLfDzqrAEQN7e8CHy9z2Q5VfbWq3iYiLwZ+Q/BP/H+EvXkR+TxwTnj+USLyYLj9fVX9JEGP+UGC3vFxwO0i8ktVPVjyWVcSxH9bCMTqg6FNnwU+qKpezG4IvpczCETtHhG5GvgQQX169MxXhXZ+M9z1HLBUVfeIyJnAjSJyahlbin51ZfbpJPaPd6/zgTOB+8JnjBpcROS/w+cs5QuqWqmITwlV3cuIN1N6LD4grNzzRef1Ai8+wud8kkAIDKOhSIQIEHg0+8OeZhGq+kMCl39M4v+gIvIt4Ilw/7ujc0RkS5n7vxX4lAaDLZ4UkaeBk4B7S+4fhWKGROSrBJ4KwCrgO2HjOB+4WETyBGGe7ap6T3jeDwhEILLlCuBS4Pzws1HVIWAo3F4nIpsJyhjXjvPo2wm8nu0ikiEYCLU33H9e7LweYA1BL3e2iGRCb6CHwAuK3yt+zbMEjefXVTWeu4h+L380jm2GYdSARISDwt7u0yLye1CoVDn9CJcRnpsWkXnh9mnAacBtFX70VoKeblQ7fiKBu1/6GVFsXAh6nQ+Hdh+jqstVdTlBQ//nqnqjqu4EtolINOr0fODR8B4XEngSq1W1P/YZ3VGsXUSOJUjQjrKlhB8BUeXPG4A7QlG5FbhAROaECeELgFvDY3eG5xJe+7+xe/1R+Ls/GzgQit/PgTeIyILQtrlheMUwjEZAVZ17EcRunwNyBD3QtwPHAD8liPE/Cnykwnu1hec/StADXznGeVvK7DuKQDAeImjYL48d+wlwVLh9R+yc64GZZe71NeANsfcrCXrxG4AbgTnh/icJYu8Phq8vhfsvI6hqWQ/cT5DAjO71L+HvyQ9/fjT27N8P73kvcGzsmreF+58E3hrbf2x47pPhta3hfiGYgmFz+KyrYtf8QWjrBmAdcHaF382LQ3v7gD3AI/X+27OXvZL2smkjDMMwmphEhIMMwzCMyWEiYBiG0cQ4Vx00f/58Xb58eb3NMAzDcIp169btVtXu0v3OicDy5ctZu3a8qkfDMAyjFBF5ptx+CwcZhmE0MSYChmEYTYyJgGEYRhNjImAYhtHEmAgYhmE0MSYChmEYTYyJgGEYRhNjImAknvu27OVnjz5fbzMMoyFxbrCYYUyU//rFU2zd28/vnLKw3qYYRsNhnoCRePK+4ttsuYZRFhMBI/F4vuL5JgKGUQ4TASPx+KqYBhhGeUwEjMST9ywcZBhjYSJgJB5PLRxkGGNhImAkHt9XzBEwjPKYCBiJxzwBwxgbEwEj8fhWImoYY2IiYCQeT00EDGMsTASMxOP5WDjIMMbARMBIPJ7v2zgBwxgDEwEj8Xi+4psKGEZZTASMxOMrlhMwjDEwETASj+crnomAYZTFRMBIPEE4qN5WGEZjYiJgJB7fSkQNY0xqKgIikhaRB0Tk5pL97xMRFZH5tbTHaA4sHGQYY1NrT+C9wMb4DhE5GngVsLXGthhNghfOHaQmBIYxipqJgIj0AJcA15Yc+gzwAcD+Q41pIfICrErUMEZTS0/gswSNfSFFJyKrgR2qun68C0XkHSKyVkTW9vb2Tq+VRuKIRgvbqGHDGE1NREBELgV2qeq62L4O4CrgI0e6XlWvUdVVqrqqu7t7Gi01kkg0UMySw4YxmkyNPuccYLWIXAy0AV3AN4BjgPUiAtAD3C8iZ6nqzhrZZTQBI+EgEwHDKKUmIqCqVwJXAojIecD7VPWy+DkisgVYpaq7a2GT0TxEYwQsGmQYo7FxAkbiiTwBywkYxmhqFQ4qoKprgDVl9i+vtS1G8tHYqmI2iZxhjMY8ASPRxNt9ywkYxmhMBIxEEw8B2ahhwxiNiYCRaOK9f9MAwxiNiYCRaIo8AcsJGMYoTASMRBMPAZkIGMZoTASMRBOvCLJwkGGMxkTASDSWGDaM8TERMBJNXASsRNQwRmMiYCSaeO/fBosZxmhMBIxEY+EgwxgfEwEj0cQXmLfF5g1jNCYCRqIpCgeZJ2AYozARMBKNJYYNY3xMBIxE49tgMcMYFxMBI9HkvbgnUEdDDKNBMREwEo1vOQHDGBcTASPR2ARyhjE+JgJGorHqIMMYHxMBoyy+r/zs0edRxxvO+ChhGydgGKMxETDK8u37tvIn/72W76/bXm9TpoSViBrG+JgIGGXZsW8AgF0HB+tsydQoWk/ARMAwRmEiYJQlajBTKamzJVOjyBOwxLBhjMJEwChL1GCmJUEiYBpgGKMwETDK4oVJ1LTjnoCNGDaM8ampCIhIWkQeEJGbw/f/KCIbRORBEblNRI6qpT3G2ESNZ8p5T2Bk2/VKJ8OYDmrtCbwX2Bh7/2lVPU1VVwI3Ax+psT3GGES9Ztc9AVtPwDDGp2YiICI9wCXAtdE+VT0YO2UGYP+lDUJSEsNJCwd5vppHY1SVWnoCnwU+ABQN2RGRT4rINuDNjOEJiMg7RGStiKzt7e2ddkONZCaGXW87Dw7mOO2jt/LrJ/fU2xQjQdREBETkUmCXqq4rPaaqV6nq0cA3gb8od72qXqOqq1R1VXd39zRba8BI4+m4I5AoT2B/X46+YY9t+/rrbYqRIGrlCZwDrBaRLcB3gFeKyPUl53wLuKxG9hhHIGovXQ8HFU8l7bYIRCG6vGfzXxjVoyYioKpXqmqPqi4H3gjcoaqXi8jxsdNWA5tqYY9xZFxvMCOSNIGcF05+lHfcozEai0ydP/9TInIiQZ7gGeCddbbHCIlCJ66PsvWLppKuoyFVILLf9bCW0VjUXARUdQ2wJty28E+DEvWgXS+rTJInkA89gZzn9nMYjYWNGDbK4ifQE3BdBCIPwLM5sY0qYiJglGWkwUlGwwnuC1r0LJYTMKpJxeEgEZkFnAjMjO9X1TuqbZRRf/xCOKjOhkyRuP3OP0skAq4/iNFQVCQCIvLHwOeBw0C8SFmBY6tvllFvkpIYjodO3H8W8wSM6lOpJ/BJ4A2qest0GmM0DlFn0/nEcCx8bjkBwxhNpTmBDHDbdBpiNBbR/DSu5wT8BK0sFtlv1UFGNalUBP4Z+DsRsURyk5CccFBy5g7KJyRZbzQWlYaD/hpYBHxARIpmr1LVpVW3yqg7hdCD4y1n0VTSjjeenmc5AaP6VCoCl0+rFUbDkRRPwFdFJPACnBcBmzvImAYqEgFVvWu6DTEai5yXjHlqPF/JpIR8AubhT8rYDaOxqCjGLyJZEfmYiDwlIoPhz4+JSMt0G2jUh2EvOeGglAgpkUQ8C7gvzEZjUWk46F+AswgmeHsGWAb8PdBFkC8wEkbkCbgeDvJ8JZ0SVEemx3aVERGwcJBRPSoVgd8DTlfVKCn8mIjcD6zHRCCRRCLgevjZUyUtgp/SRAga2Ihho7pUWvI51soibq84YoxJLh96Ao6HUHxfSaWCcJDrz2I5AWM6qFQEvg/cJCKvFpGTReRC4Ebge9NmmVFXCjkBxxscT4NwUFokEV4NWE7AqC6VhoM+APwdwfxBRwHPAt8GPjFNdhl1Joo7u59MhZQIIu57NXnLCRjTQKUlosPAR8KX0QQUwkGO9zp9X0mnIJ1KQDgoKtu1nIBRRcYUARE5V1V/EW6/cqzzbCrpZJJLSDgo7yuZVIq8qPPPUpjUz/HnMBqL8TyBLwArwu3rxjjHppJOIKrKsJeMcJCvSioFqZQkoEQ0XF7S9QcxGooxRUBVV8S2j6mNOUYjEE88uh4O8vygRDQtkoBniX5aTsCoHpOaFVREXiEiL6+2MUZjkIuV0bgefo6qg1LivlcTNf6WEzCqSaXTRtwlIueE2x8EvgN8R0Q+PJ3GGfUhl0+QJ+CFIpCExHCozVYialSTSj2BFcDd4fafAucBZxNMI2EkjOG4J+B4g5P3lXQqFQwWc/xZIk/A9e/EaCwqFYEUoCJyHCCqulFVtwFzJvJhIpIWkQdE5Obw/adFZJOIbBCRG0Rk9oSsN6aF4nCQ2w2Or8EsoukkJIbVxgkY1adSEfgV8DngX4EbAEJB2D3Bz3svsDH2/nZghaqeBjwOXDnB+xnTQDzm7HrvOV+YNsJ9Qcvb3EHGNFCpCPwxsB/YAPxDuO8k4OpKP0hEeoBLgGujfap6m6rmw7d3Az2V3s+YPoaT5AmE6wkkIhxkK4sZ00ClI4b3AB8u2ffjCX7WZwmmn+gc4/jbgO9O8J7GNJBLVE7AD0pEk5AY1mQM4DMai0qrg/5GRFaG22eLyNZwYZmXVHj9pcAuVV03xvGrgDzwzTGOv0NE1orI2t7e3ko+0pgCcRFwveH0/WDKCEnCBHJh459z/UGMhqLScNBfA0+H2/8E/DvwSYLefSWcA6wWkS0E5aWvFJHrAUTkCuBS4M06xvp/qnqNqq5S1VXd3d0VfqQxWeLhBtd7nXnfD2YRTWHLSxpGGSoVgVmqekBEOoHTgf9U1euAEyu5WFWvVNUeVV0OvBG4Q1UvD6ek/iCwWlX7J2G/MQ14RSOG62hIFYhWFrPlJQ2jPJVOJb1NRF4KnAr8QlU9EekCvCl+/ueAVuB2EQG4W1Vt7EGdiapP0qkENJw6khh2vQc9srKY48psNBSVisD7gR8Aw8Bl4b5LgXsn+oGqugZYE26/YKLXG9NP1Ni0pFPON5x5LygRjdYZdpnou/B1ZMU0w5gqlVYH/YRgMZk43w9fRsKIBiO1ZNwXAb/gCbgfS497ZZ4qKVvd1agC460nsFxVt4Tb400X/VS1jTLqS8ETSIAIjAwWc79ENJ4LyHtKNl1HY4zEMJ4n8BAjNf1PEqwdUNr1UMD+FBNGPBzkesMZHyzm+nQLXmykcPAs9q9nTJ3x1hPojG1Pasppw00iEWhNiCeQFiGTFgbzbj9LUTjI8e/FaByscTdGkY+Hgxz3BKIS0XQqAdNGxOzP2fxBRpWoKDEsIksJ5gx6ETAzfkxVT5gGu4w6Es8JDOamWgVcXzxfyaSFTEqcr6+Pi4B5Aka1qLRE9PvAJuAjwMD0mWM0AlFjmU0nwxNIhXMHud5wxu13Pb9hNA6VisBJwEtU1f7ymgA/nhh2/BuPBotlUinn59wxT8CYDirNCdwE/PZ0GmI0DvkElYh6scFizj+L5QSMaaBST+Avgd+IyGbg+fgBVX1b1a0y6ooXHyzmejgo8gTSCcgJWHWQMQ1UKgJfJZgnaCOWE0g8cU/A9YqaaI3hTAI8gbzlBIxpoFIReCVwlKoemk5jjMagME4gIYnhdArSqZT7nkCs4bclJo1qUWlOYAMwbzoNMRoHL14d5HDDqaqhCCTDE/D8wDsDm07aqB6VegJ3ALeJyFcZnRP4StWtMupKUsJBkenR8pLuVwf5tGZSDOd95wXNaBwqFYGXATuAC0r2K2AikDCKJpBzOBwUPUcmLWTTSfAElNZMmkPkbU0Bo2pUOpX0K6bbEKNxiA8Wczn/GDX6wWCxJOQElFYLBxlVptKF5rtFZGa4nRaRt4rIW0TE5h5KIF5sXV6nPYHQ9mCwWAI8AR0RAdefxWgcKm3EbwaOD7f/H/A+4G+Bf5sOo4z6ko8mXXN8SUYvtkxmNFjM5cXmPU8LiWHX8xtG41CpCJwAPBhuvxm4iKBs9I3TYJNRZwpz8IfLF7qaHI5q6dOhJxDsc/NZIPAE2sKVZFwWZ6OxqDQx7AEtInICcEBVt4ahoJlHuM5wkLgnAO4uZRiFg9IpIZ0On8V3d0WueE4gZyJgVIlKReAW4HsEYwW+E+47haBiyEgYXokn4GrDGfWW0ykhm3I/oZr3ldbwi7DqIKNaVCoCfwJcAeSAb4T75gMfnQabjDoTTbWQjsJBjsbR4yIQPYvn8Ehbqw4ypoNKS0SHgGtK9q2ZDoOM+uN54VQLMuIJuEhBBMLlJcHtOXeKRMBhMTMai0pXFptLUBG0ktEri51bfbOMehLMvJmKJYbrbNAkiQ8WS6fcFjQIbI8Swy6LmdFYVBoO+hbQSpAX6J/sh4lIGlgL7FDVS0Xk9whCSicDZ6nq2sne26gehXV5w1ywqw1OPBwUVQe5nFAtSgybJ2BUiUpF4KVAdxgWmgrvJZiOuit8/zDwu8CXp3hfo4rkw8RwOh0OTHI0J5CPh4PCxLDTOYFYiaglho1qMZFZRHum8kEi0gNcAlwb7VPVjar62FTua1SfwohhSUY4KJ1yPyfg+4oqtGUtMWxUl4nMIvrTcBbRnfEDE5hF9LPAB4DOiq0z6kLe08K0EeCuJ1C2OsjRxjP6DlozgSdgI4aNalGpCLwc2A68qmR/RbOIisilwC5VXSci503EwPD6dwDvAFi6dOlELzcmSJQTSInbI4bjg8VcHzEcT3KnxKqDjOpRq1lEzwFWi8jFQBvQJSLXq+rlFX7+NYQlqqtWrbK//mkmWpfX+d5zkSfgdmllJF7BZHjuz4hqNA4TngVUAlLRq5JrVPVKVe1R1eUE8w3dUakAGLWnUB2UGpk2wkXyXjlPwM0wyoigpcikxRLDRtWodCrpJSJyg4jsAfIEI4ej16QRkdeLyHbgJcCPReTWqdzPqA55Lxwn4Hg4yC9MJZ0qJIad92ok8AbMEzCqRaWewJeAYeB84DBwBvAj4J0T/UBVXaOql4bbN4QeQquqLlTVV0/0fkb1SYwnUOg9U3gWVxvPKBGczaTIplOWGDaqxkTGCSxV1T4RUVVdLyJvB34D/Nf0mWfUg7zv05LJFDwBV3vPfjyEElU6Ofosw/lQBNJROMjN5zAaj0o9AY8gDASwX0S6gT5gybRYZdQVT4vLKh0NoxcNFnPdE4jsbkmnyKRS5Fz9UoyGo1IRuAe4ONy+Ffgu8EOCKSCMhOH5flCFko6mWnCzwSk3bYSrCdVCOCidImuegFFFKg0HvYURwfgrgqUlOwkGgBkJIxos1hKOFsvl3Ww44yKguO0JROGgTFrIpFPOVjkZjccRRSCc9O1qwsFaqjoAfGKa7TLqiOcrmbTE1rN1s+GMLy8Z9WFczQlEnkAQDhJnvxOj8TiiCKiqJyIXANb1aBI8X0mJkI08AUdDKCMloiNLY7rqCUR2R4lhV8XMaDwqzQl8BviYiGSn0xijMYhmEc2GOYFhR0Wg3GAxz9EwSq5QHRSMGHZVmI3GY1wREJE3hZvvAd4PHBKRbSKyNXpNu4VGzfHC5SVbEuIJxCudXA2jREKcscSwUWWOFA76MvBtwKZ4aCK8gicQiMCwo4nhfCwxHA55cDaMEolXVCJqiWGjWhxJBARAVe+qgS1Gg5D3lXRayGYc9wRiIpBSt6uD8oURw0Hp7kDOzecwGo8jiUBaRF5BKAblUNU7qmuSUW883yctIyWiw46GHuKDxaK/YM9RQRsuGidgnoBRPY4kAq3AdYwtAgocW1WLjLqT9xM2TiDtfnVQFA7KpoISUcsJGNXiSCLQp6rWyDcZfpQTyETJVLdFIBEloiXhIFe/E6PxqHTEsNFEFHICjlcHRQ1+SpKQGI5NIJdKOfscRuNRUWLYaC6i6qCoB+1qTsCPeQLRjKiuhlGG4+GgtI0YNqrHuCKgqrYofJOhqmFOIIWEyeFklIgG3oCzg8Vi4aCslYgaVWTCy0saySaKMqTDnnPW4fizr0pKQMJncXlFrnw8HGSDxYwqYiJgFBH1MKNppFsy7k5REFU5RaRT7s65E4WDokF8rn4nRuNhImAUEUUZosbT5QbHKxGBYM4dN0Ug5/lk00FYy2WPxmg8TASMIgqeQEwEhvNuNjhBgnvkTzzwBNwUtFzeL1RrZdIpCwcZVcNEwCgivhALuB0OCqbEHnmfTbvbg877WhCBbFqcXe3NaDxMBIwi8iUikE2Ls9VBweI4pZ6AmyIw7I14AumUoDpSAmsYU8FEwCii1BNwOSeQ9/1ROQFXPYEgHDTynYC7az8bjYWJgFFE6VQL2XTK2UVlBnM+bdliT8DlheYLOYGU2wPfjMaipiIgImkReUBEbg7fzxWR20XkifDnnFraY4xmxBMI/jRczgkM5jzas+nCe5eranK+FjyBKMRlImBUg1p7Au8FNsbefwj4uaoeD/w8fG/UkXyJJ9CSdresciDn0RYXAYfX5o1XB0ViYOEgoxrUTAREpAe4BLg2tvu1wNfD7a8Dr6uVPUZ5ohLKVMr9EcMDw8UikHY5J1AUDjJPwKgetfQEPgt8AIi3KAtV9TmA8OeCGtpjlGEwF3w9bZmo1+nu3EGDeb/YE3C4OijnxcNBbk/xbTQWNREBEbkU2KWq6yZ5/TtEZK2IrO3t7a2ydUacgZwHQHtL0HhmMw4nhoc92ksSw642nHFPIBIDV70ao7GolSdwDrBaRLYA3wFeKSLXA8+LyGKA8Oeucher6jWqukpVV3V3d9fI5OZkYDgQgY5QBFocLhEdKEkMt2fThedzjVzROIEoHOTm92I0FjURAVW9UlV7VHU58EbgDlW9HPgRcEV42hXA/9bCHmNsIk8gCqNk00LO0WkjBksSwzNa0/Q7KwIj4aBsyjwBo3rUe5zAp4BXicgTwKvC90YdiXrKUQ/a5RLR0uqg9pZMQeRcoygxbCWiRhWp+fKSqroGWBNu7wHOr7UNxtiMygk4PVjMKzwHQEc2Td9Qvo4WTZ6c55PNRCJgJaJG9ai3J2A0GIWcQDboH7iaE8h7PjlPacvERKDV5ZyAFsJAWSsRNaqIiYBRRCEn0DJSIuriYLHBsKy1vWXkT7yjJU1/zkPVvecpDgdF00a4J85G42EiYBQxmPNISeABQCACnq/O1deX5jYAOloyeL46Gd7KeVoIB42MGHbrOzEaExMBo4j+4aCsMlqXN5txc2DSYOjRtBaJQLDdP+ReSCjn+UXCDDg7iM9oLEwEjCIGSpKpUcPjWu85EoH2ciLgYIVQzvML8znNaA3yNf3Dbia5jcbCRMAoYnC4RATCEETOsV5nYfqLkhJRgAEHG894ddDMUAQOO1rpBMGCOK55l0nFRMAoonSUbWEBE8eSwwNlPIEZobj1ORYOUtVwsFjwXUSewOFBd0Xgyh8+xPFX3VJvMwxMBIwSxhYBt3ptI+MdRv7EIw/HtVHD0cjgqES0IxuJmbsi8N212wC3nyEpmAgYRfSXTL8chYOG8m41nIXEcKa4OghgIOdWwxM1lJEHkEoJM1szHHbMoynH8wcH621C02MiYBQxmPMKCVSAzrag4TnoWOhhsGTkM7gbDjo4EPzuo+8CgnmQXO5FR8t+Pn9wqM6WGCYCRhEDJYnhrrYsAAcHcvUyaVKUGycQPZdro4YPDga/+672bGHfjNYMhx1McEdEye1dh8wTqDcmAkYRpZOuzWp32xNoy44OB7lWWhmJQNwTmNmacToxHIW2dh4wEag3JgJGEaWLs0e9zwOueQJhiWgSxgkcChv7yCuDQARcDgdFM3dYOKj+mAgYRUQjhiOcDQcVEsMjf+KtmRQpcW/EcPS7j4vAjNaM0+MEIgF73sJBdcdEwCigqgyUJIbbsmlaMqlCSMIVBnNe0OiHZZUAIkJHS8a5EtGCJ9BeHA7qcyysFedQKAK7rDqo7pgIGAWG8j6q0BYTAYBZ7VnnPIF9fcPM7siO2t/eknauRDQS4CiZClF1kFtiFjGc9wvzHu00Eag7JgJGgXLz7QB0tWUKZYqusPvwEPNnto7aP6PFvcbz0GCeGS3pwopiADNbs86Gg6JQUEpgX59bnYskYiJgFOgrU1YJQXLYtXDQ7sPDZUWgsy3rXJL74ECOzrZir2Zma7qoR+0SkXgtmdPO4aG8rYtQZ0wEjAJ7DgeVGvNKGs9Z7e41nGN5Agu7Wtl1yK2KlEOD+aJ8AIyUWLpYIRSJwFGz2gH3yo+ThomAUWBXWK63oLO48exqcysnoKrsOTzM/M6WUccWdLU5l4w8ODjaE5jh8EyicU8A3Cs/ThomAkaBqIe8oKtEBNozTvXWDg7kGfZ8ust5Ap1t7OkbdiqMcmgwT1dbsSfQGXkCDlYIFURgtolAI2AiYBTYdWgQEUaFUaJwkCtr8/aGYa2xwkHxc1xgXE/AIXGOiGyORGB//3A9zWl6TASMArsODTG3o6UwfXREV1sWz1dn6ut3jysCbYBbs1eWywl0hyE7F0fc9lk4qKEwETAK7Do4VGhc4kT19nv73OixFUSgbE4geD5X8gKerxwcyBWNFoaRBnT7vv56mDUlConh0BNwKd+URGoiAiLSJiL3ish6EXlERD4W7j9dRP5PRB4SkZtEpKsW9hjl6T00yIKwpxzn6DkdAGzd60aDs/tQJZ6AGz3o5w4MkPeVnvA7iOhqy9LVlmHH/oE6WTZ5Dg5aTqCRqJUnMAS8UlVPB1YCF4rI2cC1wIdU9YXADcD7a2SPUYZdh4ZGVQYBLJ8/A4Ate/pqbdKk2HlwiExKmNMx2hMIwl3iTDho655AeJfP6xh1bMmcDrbvc08Eeg8NMn9mC23ZNO3ZtIlAnamJCGjA4fBtNnwpcCLwi3D/7cBltbDHGI3vK71jiMCirjZaMym27HZDBB7beZAXLJhJOjZvUEQqJSzobHOmB70lFIFloRDH6ZnTzg4HRWDngcGCR+biGJSkUbOcgIikReRBYBdwu6reAzwMrA5P+T3g6FrZYxSz69AQeV9ZNGt0OCiVEpbN6+Dp3W6EgzbtPMRJizrHPH7Sok4effZgDS2aPM/s6aMlk2JxmTDdktntbN/X70zVVsRzBwZZPMtEoFGomQioqqeqK4Ee4CwRWQG8DXi3iKwDOoGymUcReYeIrBWRtb29vbUyual4eMcBAE5eXD4ts3zeDJ5xIBy0v3+Y5w4MctIYzwGwYsksNvcedmJxmWf29LN0bkfRbKgRPXPa6Rv2nGtEnz9onkAjUfPqIFXdD6wBLlTVTap6gaqeCXwb2DzGNdeo6ipVXdXd3V07Y5uIh3YcICVwylgiMH8Gz+ztx/cbu9e5aechgHE9gRcumYWvOOENbNnTx7K5o/MBQCFZ/LQjYToIJinc158reAJd7VkOODY5YdKoVXVQt4jMDrfbgd8BNonIgnBfCvg74Eu1sMcYzcM7DnBc98zCIKRSTljYyXDeLzSyjcqRPBqAF/bMAmDD9gM1sWmyDOY8nt7dx7Hdo/MBACuPng3Aumf21dCqqREl5BeF8wZ1d7ay88CAcyGtJFErT2AxcKeIbADuI8gJ3Ay8SUQeBzYBzwJfrZE9RgxVZcOOA4XGsRwve8F8AH75RGOH4+7YtIsXLJhZCDeUY2FXG4tntfF/T+2poWUTZ90z+xjK+7zkuHlljy+a1cayeR3c/dTeGls2eZ4L1xReFH4/Jy6cyb7+HL2OTeqXJGpVHbRBVV+kqqep6gpV/Xi4/2pVPSF8fUitO1AXHn3uIL2Hhjhz2Zwxz1k0q40TF3byyyd219CyibG3b5h7nt7LhacuOuK5F61YzF2P9TZ0PPoXT/SSTQtnH1teBADOPmYe923Z2/BhuogRTyAUgUWBx7axwT3MJGMjhg1uuH8H2bRw8YrF45537gnzuffpvQ07cvjHG57F85ULVxxZBFavPIphz+eWh56rgWUTR1VZs6mXVcvm0tFSPkQH8JLj5nFgIMd9W9zwBraEFWaRCES5m8d2Nn5+JqmYCDQ5A8MeNz74LK84cQFzZoweXBXn91cdzbDn87VfP10j6ypnMOfx+Ts3c+ayOZx61JEHnp/eM4uTFnXyhTWbGco33pxIv3pyN489f4jXnH7UuOe9+tRFzO7I8pUG/E7K8csnelmxpKuwVOacGS0s6mpj03PmCdQLE4Em53N3PsHuw0P86bnHHvHc4xd2csEpC/nqr7c03Jw1/3rrY+w8OMj7LjgRkdHllKWICFddcjJb9/bzHz9/ogYWVk7e8/n32x9nUVcbl525ZNxz21vS/OFZS7nt0edZv21/bQycJAf6c9y/dR/nnbCgaP/JiztZt3WfMyGtpGEi0MTc+shOvrhmM7/7oiW8ePnciq656pKTAXjn9evY1wBhIVXlS3dt5tpfPc0VL1k2ZhK1HC8/vpvfX9XD5+/czHW/erohKlRUlU/8eCMPbN3PBy86kdZM+ojX/Nm5x7G4q42//M4DDT0x3s82Po+vcN6JxWXeq1cexTN7+vlFgxcdJBUTgSbk8FCeT9+6iXddv47Tj57NP75uRcXXLps3g6vftJLHnz/Maz73K2556Lm6rBHr+8qvn9zNW667l0/dsolLT1vM3116yoTv84nXvZALTlnIP978KG/92n3c+/TeuonBlt19/Ol/r+Nrv9nCW89Zzutf1FPRdbM6svzHm15E76EhXv+F33DnY7saQtDiHBjI8S+3buKkRZ28aGlxAcIlLzyK7s5WPvuzJxjMNV5oLulIo/2xHIlVq1bp2rVr622GU/i+8szefjZs388dm3Zx+6PP0z/s8btnLOHjr11RiM9OhPu37uP931/P5t4+FnS2cuGKRZy5bA4nLuqkZ07HpO45FqrKnr5hntnTx1O9fdy/dT+/erKXbXsHmN2R5W9edQKX/9aysqNqK8H3la/8+mmu/tkTHBrKs3xeBy87fj5nLJ3Dsd0zWTq3gzkd2YrCTBP5zD19w2zuPcz6bfv5+cZd3LtlL23ZFO+74ETe/rJjJvx567ft56+++2AwtmD+DF67cglnHTOXFUu6Ri1KUys8X7nn6T189EePsLm3j/9510sL4xvi3LT+Wd7z7Qd46XHz+NjqUzl+4diD/YzJISLrVHXVqP0mAuVRVVSDWe4K7yHcFxyj5L2WnMc4xzQ4WPy+sB00EjnPJ+cFP/OF9z55b+RY3vcZzvscHspzaDDPwcEchwaD7X19w2zf18+z+wcZDnvrszuyXLRiMW988dGcXuafcSLkPJ87N+3ie2u385vNu4sWnZnVnmXejBY62zJ0tmXpbMuQSafIpISUSPAzJYgEMfCcpwx7Prl88Ix9wx4H+nPsHxhmf3+OodhykJ1tGVYtm8PqlUdx0YrFtGWPHDKphP7hPD95aCc/Wv8s67bspS/2PB0taeZ0tNDVnmV2e5YZrRlaMymyaSGbTpHNpGgJF+PxfMVTxfPCn77SN5TncPg6MJDjuQODRUtcvmDBTF7/oiVcdkZP2fmbKmUw5/HjDc/x7Xu3sjY2iGxOR5aj53awoLOVWe0tzGrPMqs9S0dLmpZMipZMitbCzzTZtCAiCCACgoQ/AYGUCJ6vDOd9hvI+w57PUM5jKO+zt2+Y5w8Osn3fAA9s3cfBwTzzZ7Zy9RtXck443qQc37tvGx+76RH6hj2OmT+D03pmsWxuB92drcyf2cqM1gxt2TRt2VTwM5MmlYJ0+DeVEiGdEtIiSArS4XuRYDsVE9Vos5rC3ug0vQh8/KZH+da9z4w07OM13g7Tkk7R1R40vLPas/TMaWfJnHaOmTeD03pmc8LCmWTS1Y8C5j2fx54/xFO9fezYP8D2ff3s6w8E6fBgjsNDefKekveDRtHzg21QMqkU2UzQmLakU2TTKdqzaWZ3ZJnT0cKsjiyLZ7WxfN4Mls7rYPm8GWVnCK0mnq9s7j3Mlt19bN0bCOmBgRwHQlE6PJQfEea8z7CnDIdVRumUjLwkELuZrZng1Zahqy3L4tltLJndztFzOzhtySzmlVn7YKrs7x/m/q37eOL5w2zd28/Wvf3sPjzMwYEcBwZy07pI/dwZLSzsauP0nlmcfew8LlyxqCKx3n14iBsf2MG9T+/l4R0HeO7gYE3/JwviULRPivaNnFN8cvE1xeeU3jcuPmP+JZc58MU3n8nLjh9bSMej6UXgpw/v5P6t+wo9mXjPZnRPZ3QPKLws2CcjX+xY9wmOx+8Tvh/nM+J/ZCkRWjIpMqkUmbQUepyZ1EjvMxP1QtMpZrZm6GzLVK1XbCSfvOczmPfD3rw30qsPf1LUSSrxZFVJpyTmRaQLnsScjhZaMtXpaOQ9n739w+w+NMxALs/AsM9gzmMw7zGY8/FDr8tXDbZ9xVfwdfR2wXuPed3B+5BwR7xFPNI1pcfjB+NRhLKfx9idzqL7xXjzby3lBQsmFyobSwSqF7htcC5csaiiQUSG0Sxk0ilmplNQfSekamTSKRZ0trGgc/IhMmN8rDrIMAyjiTERMAzDaGJMBAzDMJoYEwHDMIwmxkTAMAyjiTERMAzDaGJMBAzDMJoYEwHDMIwmxrkRwyLSCzwzycvnA427PuLEsGdpPJLyHGDP0qhM5VmWqWp36U7nRGAqiMjacsOmXcSepfFIynOAPUujMh3PYuEgwzCMJsZEwDAMo4lpNhG4pt4GVBF7lsYjKc8B9iyNStWfpalyAoZhGEYxzeYJGIZhGDGaUgRE5D0i8piIPCIi/1Jve6aKiLxPRFREJrfkUJ0RkU+LyCYR2SAiN4jI7HrbNFFE5MLwb+pJEflQve2ZLCJytIjcKSIbw/+P99bbpqkgImkReUBEbq63LVNBRGaLyA/C/5ONIvKSat276URARF4BvBY4TVVPBf61ziZNCRE5GngVsLXetkyB24EVqnoa8DhwZZ3tmRAikgY+D1wEnAK8SUROqa9VkyYP/K2qngycDbzb4WcBeC+wsd5GVIGrgZ+q6knA6VTxmZpOBIB3AZ9S1SEAVd1VZ3umymeAD8AY69E5gKrepqrRgrd3Az31tGcSnAU8qapPqeow8B2CjoZzqOpzqnp/uH2IoLFZUl+rJoeI9ACXANfW25apICJdwLnAdQCqOqyq+6t1/2YUgROAl4vIPSJyl4i8uN4GTRYRWQ3sUNX19balirwNuKXeRkyQJcC22PvtONpwxhGR5cCLgHvqbMpk+SxBB8mvsx1T5VigF/hqGNq6VkRmVOvmiVxjWER+BpRbUPgqgmeeQ+Dqvhj4nogcqw1aJnWEZ/kwcEFtLZoc4z2Hqv5veM5VBOGIb9bStiogZfY15N9TpYjITOB/gL9S1YP1tmeiiMilwC5VXSci59XZnKmSAc4A3qOq94jI1cCHgL+v1s0Th6r+zljHRORdwA/DRv9eEfEJ5uPorZV9E2GsZxGRFwLHAOtFBIIQyv0icpaq7qyhiRUx3ncCICJXAJcC5zeqII/DduDo2Pse4Nk62TJlRCRLIADfVNUf1tueSXIOsFpELgbagC4RuV5VL6+zXZNhO7BdVSOP7AcEIlAVmjEcdCPwSgAROQFowcHJpVT1IVVdoKrLVXU5wR/KGY0oAEdCRC4EPgisVtX+etszCe4DjheRY0SkBXgj8KM62zQpJOhRXAdsVNV/r7c9k0VVr1TVnvB/443AHY4KAOH/9DYROTHcdT7waLXun0hP4Ah8BfiKiDwMDANXONjzTBqfA1qB20Ov5m5VfWd9TaocVc2LyF8AtwJp4Cuq+kidzZos5wBvAR4SkQfDfR9W1Z/UzyQDeA/wzbCT8RTw1mrd2EYMG4ZhNDHNGA4yDMMwQkwEDMMwmhgTAcMwjCbGRMAwDKOJMREwDMNoYkwEDMMwmhgTAcM4AiJyuYg8KiIHRGSPiLwvduwoEdleT/sMYyo042Axw6gYEZlFMMDwTFV9KJy4Kz453MXAT+tinGFUAfMEDGN8+gmG6H9eRD4OLFHVx2PHLwZ+AiAih8P1HRCRd4YL/SwM379fRJye0thIJiYChjE+rwG+parnEjT2a8LJ+6KJ1s4lWBQHYB/QGc6/827gSWBO+P7PgP+stfGGcSRMBAxjDERkGfBx4N8AVPVu4C7gsvCUc4H14eIrAPuBmQTTez8JPAzMBi4ktu6DiPwwnLzQMOqOiYBhjM3lwP+pqhfb1wX0hduFUFDIPgIR+CuC5QAPEqxd8efAf8TOOw7YPD0mG8bEMBEwjLFZStC7B0BE5gK/zUgi+GLgx7Hz9wOrgMWquoZABE4HVhBMYY6ItAH5EmExjLphImAYY/M4cKmIdIvIAoIVz24Mq4SOAVpVdVPs/H3AXzPS6z9IsND5l2ON/klA/BrDqCsmAoYxNl8C7gceA+4Nt/8kPHYJxaEgCEQgA3wrfH+QICfwX7FzVlDFBUEMY6rYegKGMQlE5CfA5ya62IqI/BNwj6reOC2GGcYEMU/AMCbHGuDOSVx3KvBhEflZ+JpXXbMMY2KYJ2AYhtHEmCdgGIbRxJgIGIZhNDEmAoZhGE2MiYBhGEYTYyJgGIbRxJgIGIZhNDEmAoZhGE2MiYBhGEYTYyJgGIbRxPx/TqPgpvNOXlMAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot transmission for sigma+ polarization\n",
    "alpha_p = Asp(ss_1,ss0,ss1)\n",
    "Tsigma = 1-alpha_p\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(delta/w_L, Tsigma,label=r'$-1$')\n",
    "ax.set_xlabel(r'$\\delta/w_L$',fontsize=12)\n",
    "ax.set_ylabel(r'Transmission',fontsize=12)\n",
    "ax.set_title(r'$0.25\\sigma^{+}$')"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "outputs": [
    {
     "data": {
      "text/plain": "Text(0, 0.5, 'Transmission')"
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEWCAYAAACAOivfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0cElEQVR4nO2de7xcZXnvv8+amX1J9s49JCEXEjRcJEAqgWJRvOAFAbGt2uIRRW0Px9bjoe3xhlptPeUc255WaatVRK2tWlos4gVUUAx4PuViAgQCQSAYkpBA7jvZ15lZ6zl/rLVmr9l79s7svWfPmnft5/v5zCcza9Za87yZ2e/vfZ7nfZ9XVBXDMAxjZuKlbYBhGIaRHiYChmEYMxgTAcMwjBmMiYBhGMYMxkTAMAxjBmMiYBiGMYNxUgRE5Ksisk9Etjbofn8pIlujx+9O4Lq5IvJ9EdkiIo+JyHvGOO+fRORXIvJw9Fg/4v1zRcQXkbcmjs0TkW+LyBMisk1EXhYd/+vo2CMi8h0RmRcdXy0iA4nP+GLiXteJyC4R6R3xue0i8m8i8rSI3C8iqxPvXSUiT0WPqxLH10TnPhVd2xYdFxH5u+hej4jIS+v9fxwLEblQRB4UkXLy/8YwjMbhpAgA/wRc3IgbicilwEuB9cCvAx8SkTk1zttR4/L3A4+r6tnAq4C/iTvFGnxIVddHj4cT980Bfwn8eMT51wM/UtXTgLOBbdHxO4F1qnoW8CRwbeKa7YnPeF/i+PeB82rY9HvAYVV9MfDZyA5EZAHwKcL/j/OAT4nI/OiavwQ+q6prgcPRPQDeCKyNHlcD/zjG/8NE2Am8G/hWA+5lGEYNnBQBVb0HOJQ8JiIvEpEfichmEfm5iJxW5+1eAtytqmVV7QO2UL/AKNAtIgJ0RTaV67w25gPAfwD74gORCF0IfAVAVYuqeiR6foeqxp9xH7DiuEaq3qeqe2u89Wbg69HzbwMXRW15A3Cnqh5S1cOEwnNx9N5ronOJrv3NxL3+WUPuA+aJyLKoPVeKyAORh/KlSPiOi6ruUNVHgKCe8w3DmDhOisAY3AB8QFXPAT4IfKHO67YAbxSRWSKyCHg1sLLOa/8BOB3YAzwKXKOqY3VY10Vhks+KSDuAiCwHfgv44ohzTwb2A18TkYdE5EYRmV3jnu8Ffph4vSY6/24ReUUd9i8HdgFEwtIDLEwej9gdHVsIHEmIUHycsa4RkdOB3wUuUNX1gA+8ow7bDMNoAvm0DWgEItIF/AZwczhYBSDuaH8b+HSNy55T1Teo6h0ici7wn4Qd771Eo3kR+TxwQXT+iSLycPT8ZlW9jnDE/DDh6PhFwJ0i8nNVPTris64FngfaCMXqI5FNnwM+oqp+wm4Iv5eXEora/SJyPfBR4E8Tbf54ZOc3o0N7gVWqelBEzgFuFZEzathS9V9X45hO4vh497oIOAf4RdTGTiKvR0T+OWrnSL6gqvWKuGEYUyATIkDo0RyJRppVqOotwC3jXRx16NcBiMi3gKei4++PzxGRHTXu/x7gMxoWYHpaRH4FnAY8MOL+cShmSES+RuipAGwAboo6x0XAJSJSJgzz7FbV+6Pzvk0oArEtVwGXARdFn42qDgFD0fPNIrIdOAXYNE7TdxN6PbtFJA/MJQxp7SbMccSsADYCBwjDPPnIG1hB6AUl75W8Zg+hOHxdVZO5i/j/5V3j2GYYRhPIRDgoGu3+SkTeBpWZKmfXc62I5ERkYfT8LOAs4I46P3on4UgXEVkCnAo8U+Mz4ti4EMbQt0Z2r1HV1aq6mrCj/0NVvVVVnwd2icip0S0uAh6P7nExoSdxuar2Jz5jcRxrF5GTCRO0o2wZwfeAeObPW4G7IlH5MfB6EZkfJYRfD/w4eu9n0blE1343ca93Rf/35wM9kfj9FHiriJwQ2bZARE46jl2GYTQLVW3aA8gBDwE/GHH8g4Shg0V13udfCcMfJcIR6O8Ba4AfEcb4Hwc+Wee9OqLzHyccga8f47wdNY6dSCgYjxJ27Fcm3rsdODF6flfinG8AXTXu9U/AWxOv1xOO4h8BbgXmR8efJoy9Pxw9vhgdfwvwWNT+B4E3Je71V9H/UxD9+2eJtt8c3fMB4OTENe+Njj8NvCdx/OTo3Keja9uj4wJ8HtgetXVD4prfjWx9BNgMnF/nd3NuZG8fcBB4rJm/V3vYYyY8RLV5paRF5E8IQyBzVPWy6NhK4EbCMMo5qnqgaQYZhmHMcJoWDhKRFcClhB1+ks8CH2Y4wWgYhmE0iWYmhj9H2Nl3xwdE5HLCWTpbRsyOqUJEriZcgMTs2bPPOe20epcAGIZhGACbN28+oKqLRx5vigiIyGXAPg1nrbwqOjYL+Dhh0nFcVPUGwqmVbNiwQTdtGm/Ci2EYhjESEXm21vFmeQIXAJeLyCWEycg5wL8QJnNjL2AF8KCInKfh7BjDMAxjmmmKCGg4R/xagMgT+KCqviV5TlSbZ4Mlhg3DMJpHJtYJGIZhGJOj6SuGVXUj4erTkcdXN9sWwzCMmY55AoZhGDMYEwHDMIwZjImAYRjGDMZEwMg8v9hxiJ88/kLaZkwZVeWGe7ZzsHcobVOMDGEiYGSeL9/zDP/3jl+mbcaU2dMzyP++/Ql+um3f8U82jDoxETAyTzlQgiYWSpwufD9sQzlwvy1G62AiYGQeP1D8DHScfiRkfgYEzWgdTASMzBOokgENqAiZ74+1jbVhTBwTASPzlP1shIO04gmkbIiRKUwEjMzja7bCQUEG2mK0DiYCRuYJAiUDjgBBFAWynIDRSEwEjMyTFU8gDmlloS1G62AiYGSeICNTRAMLBxnTgImAkXl8zYYIVGYHZaAtRutgImBkHj/IRgjFwkHGdGAiYGQePwgysU4gboOJgNFITASMzOMHmok4uoWDjOnARMDIPIGSiZyAJYaN6cBEwMg8fqCZGD1X1glY1QijgZgIGJknDAelbcXUqawYzoCgGa2DiYCReYKMTBG12UHGdGAiYGSe7ISDbD8Bo/GYCBiZx49qB6njQhD3/ZYYNhqJiYCReYZj6SkbMkVsiqgxHZgIGJmn0nk6rgI2RdSYDkwEjMwTd5quJ4criWHH22G0FiYCRubJytTKrHg0RmthImBknniNgOt9Z6xhJgJGIzERMDKPn5H59eYJGNOBiYCRaTSxq5jrCdWshLWM1qKpIiAiORF5SER+EL3+XyLyiIg8LCJ3iMiJzbTHyD7Jft/1zlMz4tEYrUWzPYFrgG2J13+tqmep6nrgB8Anm2yPkXGSHabrs2r8ykbz6dphZIumiYCIrAAuBW6Mj6nq0cQpswH7eRsNJTn6d1wDbJ2AMS3km/hZnwM+DHQnD4rIdcC7gB7g1bUuFJGrgasBVq1aNa1GGtmiyhNwvPO0AnLGdNAUT0BELgP2qermke+p6sdVdSXwTeC/17peVW9Q1Q2qumHx4sXTbK2RJZIhINc7T5sdZEwHzQoHXQBcLiI7gJuA14jIN0ac8y3gLU2yx5ghJEMn7oeDwn9dz20YrUVTREBVr1XVFaq6GrgCuEtVrxSRtYnTLgeeaIY9xswhS4nhwDwBYxpoZk6gFp8RkVOBAHgWeF/K9hgZI9lhuj5FNLB1AsY00HQRUNWNwMbouYV/jGklOfp3fVZNVlY+G62FrRg2Mo2FgwxjfEwEjEyT3GDe9c3mKzuLOS5mRmthImBkmqpwkOOdZ+wB2B7DRiMxETAyTZYSw2orho1pwETAyDRBlhaL2c5ixjRgImBkmrKf9ARSNKQBVHICjuc2jNbCRMDINEGGcgI2O8iYDkwEjEyTpQJyldpBjouZ0VqYCBiZJkuzgyq1gxwXM6O1MBEwMk1yJo3rsXQrJW1MByYCRqbJ0hRR21TGmA5MBIxMU7WfgOMiYDkBYzowETAyTZUn4PgI2sJBxnRgImBkmupwUIqGNIA4p+F6WMtoLUwEjEyTyRXDjrfDaC1MBIxM4ydmBKnjI+jhTWXcb4vROpgIGJkmi/sJgHkDRuMwETAyTbbCQcnnbrfFaB1MBIxMk+z4Xe83q+ogOb7wzWgdTASMTJMlTyDIUGjLaB1MBIxMU11K2u2OM0uCZrQOJgJGpslSATm/ar9kt9titA4mAkamqZ5Rk6IhDSApYrbPsNEoTASMTJMlTyBLG+QYrYOJgJFpggxVEc3SBjlG62AiYGSaLBWQS2qYiYDRKEwEjExTvcAqPTsaQZb2RjBah3y9J4rIXOBUoCt5XFXvarRRhtEo/MSqKtc9Ad+miBrTQF0iICLvBj4P9AL9ibcUOLnxZhlGY6iaVun46FlNBIxpoF5P4Drgrar6w6l8mIjkgE3Ac6p6mYj8NfAmoAhsB96jqkem8hmGkSTI4M5i4H5bjNah3pxAHrijAZ93DbAt8fpOYJ2qngU8CVzbgM8wjArZqh00/Nw8AaNR1CsCfwl8QkQmnUgWkRXApcCN8TFVvUNVy9HL+4AVk72/YdQiS9MqrYCcMR3UGw76Y2Ap8GEROZh8Q1VX1XmPzwEfBrrHeP+9wL/VekNErgauBli1qt6PM4yw4xQJvQDXRcAPEm1x3a0xWoZ6ReDKqXyIiFwG7FPVzSLyqhrvfxwoA9+sdb2q3gDcALBhwwb79Rt14wdKTgQfdX43rkChkPMolgPnBc1oHeoSAVW9e4qfcwFwuYhcAnQAc0TkG6p6pYhcBVwGXKSu/5UaLYcfKDlP0MD90XMQKG2RCLg+08loHeqK8YtIQUT+XESeEZHB6N8/F5G2eq5X1WtVdYWqrgauAO6KBOBi4CPA5araP+5NDGMSxCKQE8H1wXOgSiEnQHWJbMOYCvWGg/4KOA94H/AscBLwp8AcwnzBZPkHoB24U0QA7lPV903hfoZRha9hOCjwNBOLxfK5cNxmnoDRKOoVgbcBZ6tqnBT+pYg8CGxhgiKgqhuBjdHzF0/kWsOYKEGgeJ7gqTjfcapCWyQClhMwGkW9Uz5lgscNoyXwdTgc5Pp+An4wHA5yPb9htA71isDNwPdF5A0icnoUy78V+Pdps8wwGoAfgCeCiPshFD9IhIPMEzAaRL3hoA8DnyCsH3QisAf4V+Avpskuw2gIQaDkPMh5WQgHKQULBxkNpt4pokXgk9HDMJzBVyXveZRFne84fVXaonCQ64JmtA5jioCIXKiq90TPXzPWeVZK2mhl/EDxPPC8LEwRpeIJ2B7DRqMYzxP4ArAuev6VMc6xUtJGSxOvGM6JOB9HDwILBxmNZ0wRUNV1iedrmmOOYTSWeHaQJ+7PqPFVKeRtnYDRWCZVFVREXi0ir2i0MYbRaHw/EoEMJIbDshG2YthoLPWWjbhbRC6Inn8EuAm4SUQ+Np3GGcZUKQdKzvPwshAOSuQEXBc0o3Wo1xNYR1jvH+C/Aq8CzicsI2EYLUugSt6TaIpo2tZMjUCVtrwlho3GUu86AQ9QEXkRIKq6DUBE5k+bZYbRAMpx2Ygs5AQsMWxMA/WKwP8jLPa2DPgOQCQIB6bJLsNoCEEQegLZCAeZCBiNp95w0LuBI8AjwKeiY6cB1zfeJMNoHOUgCKeIZiExrFQSwyYCRqOod8XwQeBjI47dNi0WGUYDCYKwZIRkpoCceQJGY6l3dtCfiMj66Pn5IrIz2ljmZdNqnZEaT+87xuqP3sb2/b1pmzIlykEQVhH1cHp7ydj2QkYSw1++5xku+IwVG2gF6g0H/THwq+j5/wH+FriOcPN4I4Pc8uBzAPzw0b0pWzI14p3FPBGnE8PxyL+QkSqi192+jeeODKRthkH9ieG5qtojIt3A2cBrVdUXkb+ZRtuMFIk7nZw3qfWELYOvw4lhl0MosYBVFos53JYksUgb6VGvCOwSkd8AzgDuiQRgDuBPn2lGmsSdTN7xP9CyH04RzXmCw45AxfbKdNeMiECxHNDZlkvbjBlNvSLwIeDbQBF4S3TsMuCB6TDKSJ9hT8BtEQgqnoDbHWfl+xAh73lOh7aSDJV9E4GUqXd20O2Em8kkuTl6GBmkHIRTaeLtDF1leLGY21NEY9vDOkhuC1qSobLjU7YywHj7CaxW1R3R8/HKRT/TaKOM9Ik7Gc91TyCxWCwWNheJTZfYE8iKCJTc/U6ywniewKNAd/T8acK9A0b2CAqYL5dB/KzkBKL9BPI5YbDsbscZh39ygvOhrSRDZUsrps14+wl0J567PUXEmDBxqWJP3BaBePZJznO7bETsxeRzHvmc+55A3hPKgTJonkDqWOdu1CSeHeRyHB1CEcjnpNLpuErSMwtDW+62BYYnHJgnkD51JYZFZBVhzaBfA7qS76nqKdNgl5EycadTcnzzEj9QvKh2kMuj59gzy3mhoLns1UAoZkNYYrgVqHeK6M3AE8AnAVvmNwOIww8ud5wwvFgs73mUHC4eVFm3kQsFzTwBo1HUKwKnAS9TVXf/iowJEfeXrnc2fmKxmMuC5keinPO8qC1u/ynmo/IXNjsoferNCXwfeOV0GmK0Fn7FE3D7j7TiCeTcHj3HtheicJDjUbqEJ+D27ysL1OsJ/A/gP0VkO/BC8g1VfW/DrTJSJ+50XO44YXixWN5xTyCZE/Cy4AlEIjBYsnBQ2tQrAl8jrBO0DcsJzAjiTqfs+JDTD+I9hj2nBc1P5ARcFzQwT6CVqFcEXgOcqKrHpvJhIpIDNgHPqeplIvI24M+A04HzVHXTVO5vNI44Mexyx6mq0ToBz/mOM/4ecp7nfEVUGPYELDGcPvXmBB4BFjbg864h9CZitgK/DdzTgHsbDaRYdj8nEPeT8faSTs8OimyP8xuui0DFE7DEcOrU6wncBdwhIl9jdE7gq/XcQERWAJcSbkbzJ9G126L36rXXaBKxm+6yJ5AMoRQc7ziztlgsxsJB6VOvCLwceA54/YjjCtQlAoS7kH2Y4XpEdSMiVwNXA6xatWqilxuToCICDucEKkXwxP2cQHlETiALK7nBwkGtQL2lpF89lQ8RkcuAfaq6WUReNdHrVfUG4AaADRs2uP3rd4R41obTo2cdHj27nhNI7vTmeeK0OMOwqFntoPSpd6P5xSLSFT3Pich7ROSdIlJvTuEC4HIR2QHcBLxGRL4xKYuNpjAcDnL3j9SPi+AlFou5utl8cqc31wUNhj1M8wTSp95O/AfA2uj5/wY+CPxPoK49hlX1WlVdoaqrgSuAu1T1ygnaajSRoQx4ApXKm1HHGR5zsz1xYjiuiOr6zmLlSjjI3UFGVqhXBE4BHo6evwN4I+G00Sum8uEi8lsisht4GXCbiPx4KvczGkcmcgKJ3bhy0Q5propa0hNwvQQGDM86s9lB6VNvYtgH2kTkFKBHVXdGoaCu41w3ClXdCGyMnn8H+M5E72FML2U/yMSK4eQ+yQUvHO+42p7hmU7ur3mApCdg4aC0qVcEfgj8O+FagZuiYy8hnDFkZIxiYj69q50mVItAPC/dd9SzKY+YIuq6CPgWDmoZ6hWB3weuAkrAv0THFhGu9jUyRnLGhsuLxSoiEG0vCe4muoeriGZjsVgcZrTaQelT7xTRIaIpmoljG6fDICN9ki660zmBETX4k8dcI2ueQCzG5gmkT707iy0gnBG0ntE7i13YeLOMNEkm67ISDopnB5Ucbc/IncVcnh0UBFop6WEikD71hoO+BbQT5gX6p88coxVI/mG6LALlZDgoSgw7nxPIeeHqZ0fbAVQJmCWG06deEfgNYHEUFjIyTvIPMxM5AS87OYFwiqi7YS2ott1WDKfPRKqIrphOQ4zWIVlt0+kRZ63ZQY52nuWqtnhOh4Oqf18mAmkzkSqiP4qqiD6ffKPeKqKGOxTLceE1dztNqF4s5vqK4TiMlSVPwBMoOTzIyAr1isArgN3A60Ycn0gVUcMR4pBJRyHnbCIVRnoC0WIxRzudUlWS23NaBGIh7ijknN7jISs0pYqo4RZxR9lZyGUnJ+C5nxPIeYKI+1NE/YQI9AyUUrbGqDcnUEFCvPgxHUYZ6RKPzjoKOWdHzpDciMWrJIZd7TzL0V7JgPOLxSqeQD70aAKH25IF6i0lvVxEviMiB4Ey4crh+GFkjDhO21HIRtgh5w1vZ+hyTiAWAdc9gTgZ3NGWA6DkqHeWFeodyX8RKAIXAb3AS4HvAe+bJruMFEnmBFzubIKKCHjD6wQcbU850IqQ5T1xNqwFw0LcWQhFwGVvMwtMZJ3AKlXtExFV1S0i8nvAfwJfnj7zjDQoJXICRwfddfaSi8Wc9wQCJZ8LhSznCYGCqjq5P3cyJwAmAmlTryfgE4aBAI6IyGKgD1g+LVYZqRLnBDrbcs6usIUxEsOOzkYpR4lhwP01D361J1B09DvJCvWKwP3AJdHzHwP/BtwCbJoOo4x0iTvK9nzO2ZEzjLVi2M32lH2lMFIEHF0wNuwJxHs8mAikSb3hoHcyLBh/RLi1ZDfwucabZKRNVhLDycVi8c/X1fb4gVZ2R3PdEyglck5g4aC0Oa4IiEgOuB64GkBVB4C/mGa7jBSpWizmsKuerLcT46wnEGgluZ13XARG5gRc/o1lgeOKgKr6IvJ6wL6pGUKparGYmx0NVJdfjnF18ZufmB3kidsiMDInYKUj0qXenMBngT8XkcJ0GmO0BsnEsKsjZ4AgudF8vJ+Aox1OyQ+qFouBuyIwMidgnkC6jCsCIvL26OkHgA8Bx0Rkl4jsjB/TbqHRdMq+4gkUMrIyNZkYdrU9mfIERuYEHG1HVjheOOhLwL8CVzbBFqNFKPkBhXjzkkCdnY8eJETAU8dnByXWCbheEbXsW06glTieCAiAqt7dBFuMFqHkK4WcV5WAjEfSLpFcLEZkvu9oh+Mnage5Pjto5IphE4F0OZ4I5ETk1VT+hEajqnc11iQjbcpBQD5XPbc+n0vZqElQWSeQy8LsoOwsFrMVw63F8USgHfgKY4uAAic31CIjdUp+OB0xK1MRszBF1A9C7wzcXyw2nBOwxWKtwPFEoE9VrZOfYZT8gLZcYiMWRzvO2G5PhDil4aqglXylo5AtT6BSNqLsZjuyQr0rho0ZRNkPyCdyAq7W2wkSnkA8o8bV0EPSE3DdQ6skhtvi2UFu/r6yQl2JYWNmUYoSwa6POJNTRCXyBlxdLJYsJe3+FNF4UxnLCbQC464TUNXuZhlitA6lckBbzqPgeNG1QMP1DvH01rAOv5tt8YMsLRYbXowIVkU0bZq6PaSI5ETkIRH5QfR6gYjcKSJPRf/Ob6Y9Rm3KFU/A7aJrydEzhB5BFtpSCW053BZIJIbNE0iVZu8RfA2wLfH6o8BPVXUt8NPotZEyYYkCz/lFSf4IEch7nrNlI8qJ7SVd3yXNHxkOcjRElxWaJgIisgK4FLgxcfjNwNej518HfrNZ9hhjU/aVtpw3vBuXo+66n6i8CbEn4G5bYs/M9VxNpUBhmxWQawWa6Ql8Dvgw1dVIl6jqXoDo3xNqXSgiV4vIJhHZtH///mk3dKZT8qPFYhnwBBKOAIWcuzmBchBUcjSui4A/onaQrRhOl6aIgIhcBuxT1c2TuV5Vb1DVDaq6YfHixQ22zhhJKapTE9eqcTVmm9yXF9zOCSRDW+4vFhuZEzARSJNmrRO4ALhcRC4BOoA5IvIN4AURWaaqe0VkGbCvSfYY41D2AwqeVEaers7eKAdBJYkKYSzdXU8gmRNwO0xXLIczndpycSlpN7+TrNAUT0BVr1XVFaq6GrgCuEtVrwS+B1wVnXYV8N1m2GOMT1xFNHbXh0p+yhZNjsFSQGdbtSfgasdZ9odzAvH3Mlhysy2DpYCOQg6RMORo4aB0afbsoJF8BnidiDwFvC56baRM2Q+niLbnw5/HUNnNP9LBkl8pTQBurxOIi/rBcLmFQUfFeajsV35beYfzNFmh6WUjVHUjsDF6fhC4qNk2GONTCqo9AVc7m4GSX2kDhB2OyzmBOAwUx9IHHP1eYk8AoOB55gmkTNqegNGClMpKIQOewECxWgRyGckJxDV3XBXnobJPeyRkhbyJQNqYCBijCEMP7nsCg+Wg2hNwdHZQECiqDOcE8o5/L6WA9qgNeU+cnX2WFUwEjFGUfKXgue8JDBZ9OgvViWEXR52laF59nBMoRMX9XA0HDZX9SkirkHN3FXdWMBEwRhGXknbdExgYkRjuLOQYKLrXltjm+PsQkagt7gkawFApGJEYdrMdWcFEwBhFvMew857AiMTw7PYc/Q6KQGzz7LbhtnQUPAbL7rUFYk8gSgznLCeQNiYCxihKUYmCeGMZlz2BpAh0tuWdDKHEItBZJQI5Bh0UNIhzAsMb5Fg4KF1MBIwq/CgJGRdea897znoCQ6WgquOcVcjRN1RO0aLJ0V8MbZ7VNjyju7OQy4wn4OoCvqxgImBUEbvmhXw8Jz3npCdQ9gOKflCZSQMwq93NnEDtcJCbbYFqT6CQM08gbUwEjCoqIpAoUeBieYLBcrx71fBPfFZbjv6SjzpWeG2gRjios5BzMrQFMJjwBPKWE0gdEwGjinjOdjwdMQwHudfZxN5LMicwqy2PH6hzBfFiTyAZDmoveE6KM1TPDnK5vHdWMBEwqhiekx7lBBz1BEZOq4TQEwDoH3JL1PoqOYFqT8DFMJ2qVnsCnuUE0sZEwKgijs+2ZcQT6KwlAo51ngMVTyA508nNcFDJDyceJBPDRcsJpIqJgFFFPCrLV3ICHkMOegKx9zJyiijAQNGtGUK1wkEdeTc9gXhGU1U4yDyBVDERMKoojcoJ5Jz0BAZqeALx7Jo+x8JB/cUyIsPVQyHyBBycHRQPKNptsVjLYCJgVBHPSZ8djTo7HE1AVkSgrbrjBJxbNdxf9JkVbcIS4+ysrVK1JzCrzc1V3FnCRMCooncwFIGujlAEXPUEhjub6tlBAAMl98JBnW3VW390FDyKfuBcVdR44WEcputqz9Pr4AK+LGEiYFRxLPqD7Gp32xOoJIbbshEOmt2eqzrm6u5iIz2Bro48/UXfOTHLEiYCRhWxJ9DtuCdQa4poLAiuxdL7i9XVUGG4Xa7NEKrlCQDmDaSIiYBRRW/WPIHC6HBQv2OzgwaKftX0UHDXExga4QnEgw0TgfQwETCqqIjACE/AuVILkXBlYZ1AGA4akRNwdIvJ0Z5AARj2QI3mYyJgVHFssExbzqskVDsKHoHiXJGvgREjzvi5J+6tGK4ZDora5drGMsPlPIZzAgC9Q6XUbJrpmAgYVRwbLFVcdBieXeNaXmCw5NOW9/C84WmVIsKstrxzUxL7a4WD2tzOCcS/q/i3dtQ8gdQwETCq6B0qV0ZnMDxicy0vcLivyPxZhVHHw3ILbnU4taaIzutsA+BwfzENkybNSE+gO04MmwikhomAUUXvYLmSFIbhEZtrsecDvUMs6mofdXx2W87NKaIjPIFF3aEIHOgdSsOkSTNy/UaXJYZTx0TAqOLYULUIuJqAPNBbrCkC3R0FegbciT8XywH9RZ/ujmqvZuHssG0HjrnlCYxM2HeZJ5A6JgJGFb2D5aqcwMLZ4YjzYJ9bnc1YnsCSOe3sO+bO6Hl/NNJfMqe6LW15j7mdBec8gUN9Q3QWcpWcRlye5Jh5AqlhImBU0TvCE1jY5V7YQVU52FushEySnDCng31HB1OwanK8ENm6ZE7HqPcWdbU59b1A6KHFvykAz5OwdIR5AqlhImBUcWywVJUYjkfTB3vd8QSODpQp+gGLa3kC3R0c7CtSLLuR6I4F64Q5o9uyqKvdQREY7aF1tec5NuhOiC5rmAgYFVSV3qFyVfx5/qw2PHHLE4hDKGOFg5LntDovHI3DQTU8ge52DjgkzlA7V9PdYUXk0sREwKgwVA4o+VoVDsp5woLZbU51NgfGFYGwM33BkZDQC0cHyXvCglmjQ1uLu9o54FB+A2JPoLotXSYCqdIUERCRDhF5QES2iMhjIvLn0fGzReReEXlURL4vInOaYY9Rm/gPMZkYBvfCDnHoqnZOIBQGV/ICLxwd4oTu9qpFbzGLuto4NlR2ZuZWECiH+kZ7AmE4yEQgLZrlCQwBr1HVs4H1wMUicj5wI/BRVT0T+A7woSbZY9RgfzSqnD9i1Lmoq52DDonAeJ7A0oon4EZ79h0bZMnc0aEgSORrHJm5dWSghB/oKE9g4ey2ym/PaD5NEQEN6Y1eFqKHAqcC90TH7wTe0gx7jNrsPNQPwEkLZ1UdX9jlVjhob08YQhkpZhAKXCEnToWDlnTXFoHYq3m+Z6CZJk2aWJwXjhDnVQtmsbdnwJlkfdZoWk5ARHIi8jCwD7hTVe8HtgKXR6e8DVg5xrVXi8gmEdm0f//+ptg7E9l5MBKBBbOrjrsWDnryhWO8+IQucjVCKJ4nnNDdwXNHWr/jVFX2HBkctUYgZu0J3QD88vnemu+3GmN5aKsWziZQnPhOskjTREBVfVVdD6wAzhORdcB7gfeLyGagG6g53FTVG1R1g6puWLx4cbNMnnE8e6iPOR155o6oubOoq53+ou9MHf4n9h7ltKXdY75/2tJuHt9ztIkWTY6dh/rpHSpz6tLaqbIV8zvpas/zxPOt3xag4k0uHpGriT3PZw/2Nd0mI4XZQap6BNgIXKyqT6jq61X1HOBfge3NtscYZuehAU5aOHvU8aVzw5Hb7sOtP1Lr6S+xp2eQ05aNPcdg3fK5bN/fS1+Lz0h59LkeAM5cPrfm+yLCaUu7eWLvsWaaNWl2Hw49zRNGTHddtSAUgV1RONJoLs2aHbRYROZFzzuB1wJPiMgJ0TEP+ATwxWbYY9Rm58G+yh9kkjNODDuhrVGn1MrEo+LxPIEzl88lUHh8b2uPoB99roe2nMcpS7vGPOe0Zd1se/6oE5v+PPbcUVYtmMWcEXWQTuhup6Pg8exBE4E0aJYnsAz4mYg8AvyCMCfwA+DtIvIk8ASwB/hak+wxRuAHyu7DA6xaOFoETl40m46Cx9bnWrvTBNi6JxaBsT2BM1eEovbo7tYWtUd393Dq0u5Kxc1anL5sDscGy+w61Ppe2qPP9bBu+ejvRURYtWBWZWKC0VyaNTvoEVX9NVU9S1XXqeqno+PXq+op0eOj6sJwJqNs23uUcqCcsmT0qDOf83jJsjlOeAJ3PfECL1o8m6VjTKuEcMHYiXM7uPeZg020bGIMFH0e3nWE9SvnjXvey05eCMBPn3ihCVZNnp7+EjsP9bNujNDW2iXdPLK7xwmPJmvYimEDgLufDGddvfzFtRPvZy6fy2N7evCD1v0jPdxX5L5nDnHxuqXHPfeNZy7j7l/up6e/NWvW/GTbC/QXfS45c9m45528uItTl3Tzw63PN8myybF1z/j5jVe8eBHPHx3kqX1uzHTKEiYCBgD3PLmflyybw+Lu2tMRz12zgL6iz73bW3f0fNuje/ED5Q1nHF8ELj/7RIp+wO1b9zbBsonz3YefY8mcds5bs+C4575h3VJ+seNQJfHaitz5+Au05TzOWjGv5vsXnhIOPu550qaANxsTAYPnewbZ/OxhXnnq2NNvX3v6EuZ05Ll5864mWlY/gyWfz//sadavnDfmaDPJWSvmcvqyOXxh49MtV3Zhy64j/GTbPn5nw8qaax1GcsW5Kyl4Hn//06ebYN3EGSr7fPfh53jdGUuY2zl6y0+AE+d1svaELn7wyF4LCTUZEwGDL2wMO4//ct6qMc/pKOR48/rl/HDr8y05n/uzP3mSvT2DfOgNpyJy/I5TRPj4Jaez69AA1//0qSZYWB+DJZ9P3LqVRV1tXH3hyXVdc+K8Tt5x/ipu3ryrJT21b2/ezeH+Em87Z8W4573rN1bz8K4j/PypA02yzAATgRnPvdsP8q37d/K2DStZWWN6aJI/fPWLaM95fOjmR1pm9Kyq3PjzZ/jS3c/wjl9fxQUvXlT3tS9fu4grzl3JP27czpfveSb1EWjvUJn3f/NBtu7p4brfOnPUlpLj8SevO4U1i2bzh9/czEM7D0+jlRPjVwf6+MztT3D+yQu4cO34Cz1/Z8MKls/r5JPf3cqRfnfKlLiOicAMRVW5/dG9XP3Pmzhp4SyuveS0416zbG4nn/7NM3hgxyHe/bUHUl3co6ps2nGI9/7TL/iL27Zx8RlL+dSbzpjwfT795nVcfMZSrrt9G+/66gP85/YDBE1Ofg+VfW596DneeP09bHxyP//rzevqymsk6e4o8NV3n0t3R4Hf+dK9XHfb46nnCH7+1H5+50v3ks8Jf/WWs2tWQk3Sns/xd29fz54jg7z9y/ezfb8liZuBpD36mSgbNmzQTZs2pW2Gk/iB8qsDfWz85T6+v2UPW3b3cObyuXzxneewfF5n3fe55cHdfOLWrQyVA155ymIuOXMZ61fOZc2i2vV6GkHfUJkdB/t4Zn8fD+48zM+fOsDT+3rpbs9zzWvX8t4L1hy3kxmLIFC+fu8OPnvnkxwdLLNificvf/Eifm3VPNYs6mLVglljlnOeKIMln31Hh3h6/zGeeqGXB3ce5v89dYC+os+pS7r59JvP4NejaZ+T4VBfketu28YtD+0G4Pw1Czn/5IWctWIuqxfNZsX8Tgq56Rn7DRR9Ht97lId2Hua2R/fy0M4jrFk0mxveeQ5rl4y9eG8kG3+5j2tuepjeoTKvO30Jl5y1jPPXLBi10tiYGCKyWVU3jDpuIjB5VBXVsBxq/P8YPgclfI/E6+Hn0bXRa+Lz0JrvB6qU/ICyr5SDcOOXsq+UguiYH1AKlFI5oK9Y5thgmd6hMr2DZQ73F3m+Z5C9PYPsONhHfzEM47xk2Rzeft5KfvfcVbTlJ94p7DkywL/c9yy3PLi7UpY55wlLuttZNq+TeZ0FZrfnmd2eY1ZbnrwneJ7gCeREEBFECNvhBxT9oNLGvqJPz0ApfPQXOTJQ4khiKmdHweOlq+bz5vUn8qazT2RWW34sMyfEQNHnx489z/e27GHTjkMcTdS4z3kS1lXqLDC3sxC2KSfkPSHneeQ9IZ8TAlWKZQ3bUw7bNFQOONxf5HBfkb5idRht+bxOXnnqYt5wxlJe8eJFDREaCIux/dsDO/nJtn3RiuLhdiyc3cb8WW3MnVVg/qwCnYUcbXmP9nz8r0c+50H0GwxUCaLfoSoUywF9Q2X6iz69Q2UO9A7xfM9gVUnrU5Z0ccW5q/gvv76KjsLYi93GYv+xIb5493a+t2VPpcz0oq52Vi7oZNncDpbM6WBOR4Hujjxd7Xm6OvLMastRyHnRQ0Y9z+c8chL+BuPfnxe/RhAPhPhY+H58TnxcomtdZMaLwHW3Pc5Nv9gFWulzqzraZAccvlm7UyY63wU8gbmdBZbO7eTEuR2sXDCLl5w4h3NXL2DNotE1giaDHyjb9/eyZdcRdhzsY2/PIHuPDHJ0sETfUJm+ok//UJlyEP5f+qqVzgTCTqmQEwqeRyEf/sF2FnJhZzurLep08yyb28maRbM5aeEs1p7QPSnhmghBoOw81F957O0ZoGegxNGBMkcHS/QP+ZSDgHKgFXEu+xq1J2xLe86jkBfach7zZ7Uxf3YbC2a3sairjRct7mLtCd2jivVNBz0DJZ564Rg7Dvaz40Af+48Ncbi/yJH+EkcGigyWAorlgKGyT7EcC3L4BXkyulMs5Dy62sNOd3Z7noWz2yq/sbVLujl75VyWza3fsxwPP1Ae3nWYLbt62Lb3KHt6BtjbM8i+o0Op7kbmjRCD+Fl8SJDqN8Y5R2qeIyMvB4F/fMc5vHxt/XmvqsvHEIHGDKEcYP3K+ZR8DX/ISOKLoKLuEh1IfjmV96NrKl9KdH6t95KjhTE/L3lMpPra+Fj03PPCTjKfE/I5j4IX/pvPDR+PRzvxqKi7vUBHwZv2UUvOE05Z0s0pE3D3YdiLatTIt9F4nrB60WxWN0gs02RuZ4ENqxewYfXx1xzEqGpLjHhznnDOSQs456TRtvuB0lcMPd7eoTJ90WCjFIlYqRxQDgKKsbfsB/hBOLgLdPg3GHs6w559teejlfejY1CVN0p6+eFrql4nz2HUOYn71HF9XMyxkcwYEbj0rGVcetb4qy+N5hGLnNGatIIAHI8wRFcYVZDOmBg2O8gwDGMGYyJgGIYxgzERMAzDmMGYCBiGYcxgTAQMwzBmMCYChmEYMxgTAcMwjBmMiYBhGMYMxrmyESKyH3h2kpcvArJSrNza0npkpR1gbWlVptKWk1R1VD1v50RgKojIplq1M1zE2tJ6ZKUdYG1pVaajLRYOMgzDmMGYCBiGYcxgZpoI3JC2AQ3E2tJ6ZKUdYG1pVRrelhmVEzAMwzCqmWmegGEYhpHARMAwDGMGMyNFQEQ+ICK/FJHHROSv0rZnqojIB0VERWRy+86ljIj8tYg8ISKPiMh3RGRe2jZNFBG5OPpNPS0iH03bnskiIitF5Gcisi36+7gmbZumgojkROQhEflB2rZMBRGZJyLfjv5OtonIyxp17xknAiLyauDNwFmqegbwf1M2aUqIyErgdcDOtG2ZAncC61T1LOBJ4NqU7ZkQIpIDPg+8EXgJ8HYReUm6Vk2aMvA/VfV04Hzg/Q63BeAaYFvaRjSA64EfqeppwNk0sE0zTgSAPwA+o6pDAKq6L2V7pspngQ8zvDWpc6jqHaoa7xp+H7AiTXsmwXnA06r6jKoWgZsIBxrOoap7VfXB6Pkxws5mebpWTQ4RWQFcCtyYti1TQUTmABcCXwFQ1aKqHmnU/WeiCJwCvEJE7heRu0Xk3LQNmiwicjnwnKpuSduWBvJe4IdpGzFBlgO7Eq9342jHmUREVgO/BtyfsimT5XOEA6QgZTumysnAfuBrUWjrRhGZ3aibZ3KjeRH5CbC0xlsfJ2zzfEJX91zg30XkZG3RubLHacvHgNc316LJMV47VPW70TkfJwxHfLOZtjWAWruyt+TvqV5EpAv4D+CPVPVo2vZMFBG5DNinqptF5FUpmzNV8sBLgQ+o6v0icj3wUeBPG3XzzKGqrx3rPRH5A+CWqNN/QEQCwqJM+5tl30QYqy0iciawBtgiIhCGUB4UkfNU9fkmmlgX430nACJyFXAZcFGrCvI47AZWJl6vAPakZMuUEZECoQB8U1VvSdueSXIBcLmIXAJ0AHNE5BuqemXKdk2G3cBuVY09sm8TikBDmInhoFuB1wCIyClAGw5WGFTVR1X1BFVdraqrCX8oL21FATgeInIx8BHgclXtT9ueSfALYK2IrBGRNuAK4Hsp2zQpJBxRfAXYpqp/m7Y9k0VVr1XVFdHfxhXAXY4KANHf9C4ROTU6dBHweKPun0lP4Dh8FfiqiGwFisBVDo48s8Y/AO3AnZFXc5+qvi9dk+pHVcsi8t+BHwM54Kuq+ljKZk2WC4B3Ao+KyMPRsY+p6u3pmWQAHwC+GQ0yngHe06gbW9kIwzCMGcxMDAcZhmEYESYChmEYMxgTAcMwjBmMiYBhGMYMxkTAMAxjBmMiYBiGMYMxETCM4yAiV4rI4yLSIyIHReSDifdOFJHdadpnGFNhJi4WM4y6EZG5hAsMz1HVR6PCXcnicJcAP0rFOMNoAOYJGMb49BMu0f+8iHwaWK6qTybevwS4HUBEeqP9HRCR90Ub/SyJXn9IRJwuaWxkExMBwxifNwHfUtULCTv7jVHxvrjQ2oWEm+IAHAa6o/o77weeBuZHr/8b8PfNNt4wjoeJgGGMgYicBHwa+BsAVb0PuBt4S3TKhcCWaPMVgCNAF2F576eBrcA84GIS+z6IyC1R8ULDSB0TAcMYmyuBe1XVTxybA/RFzyuhoIjDhCLwR4TbAR4l3LviD4G/S5z3ImD79JhsGBPDRMAwxmYV4egeABFZALyS4UTwJcBtifOPABuAZaq6kVAEzgbWEZYwR0Q6gPIIYTGM1DARMIyxeRK4TEQWi8gJhDue3RrNEloDtKvqE4nzDwN/zPCo/yjhRudfSnT6pwHJawwjVUwEDGNsvgg8CPwSeCB6/vvRe5dSHQqCUATywLei10cJcwJfTpyzjgZuCGIYU8X2EzCMSSAitwP/MNHNVkTk/wD3q+qt02KYYUwQ8wQMY3JsBH42ievOAD4mIj+JHgsba5ZhTAzzBAzDMGYw5gkYhmHMYEwEDMMwZjAmAoZhGDMYEwHDMIwZjImAYRjGDMZEwDAMYwZjImAYhjGDMREwDMOYwfx/06PIJXfuxMwAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot transmission for the mixed polarization\n",
    "alpha_mixed = A(ss_1,ss0,ss1)\n",
    "Tmixed = 1-alpha_mixed\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(delta/w_L,Tmixed,label=r'$-1$')\n",
    "ax.set_xlabel(r'$\\delta/w_L$',fontsize=12)\n",
    "ax.set_ylabel(r'Transmission',fontsize=12)"
   ],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}