{
"cells": [
{
"cell_type": "markdown",
"id": "2a88d805",
"metadata": {},
"source": [
"# REDOR"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"id": "bfadf6fc",
"metadata": {},
"source": [
"Rotational-Echo DOuble Resonance NMR (REDOR$^1$) is an important technique for dynamics measurements,$^{2,3}$ where applying REDOR to one-bond dipole couplings allows one to compare the known size of the rigid dipole coupling, to its measured size, which is reduced due to dynamics. This yields the order parameter, $S$, according to\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"|S|=\\left|\\frac{\\delta_{exp}}{\\delta_{rigid}}\\right|\n",
"\\end{equation}\n",
"$$\n",
"\n",
"We do not obtain the sign of the dipole coupling from REDOR, so we only get the absolute value of $S$. \n",
"\n",
"Here, we investigate how REDOR behaves as a function of correlation time.\n",
"\n",
"Note, this was the topic of a recent paper by Aebischer et al.,$^4$ where simulations were performed using Gamma$^5$. We hope SLEEPY makes this type of simulation more accessible.\n",
"\n",
"[1] T. Gullion, J. Schaefer. [*J. Magn. Reson.*](https://doi.org/10.1016/0022-2364(89)90280-1), **1989**, 196-200.\n",
"\n",
"[2] V. Chevelkov, U. Fink. B. Reif. [*J. Am. Chem. Soc.*](https://doi.org/10.1021/ja902649u), **2009**, 131, 14018-14022.\n",
"\n",
"[3] P. Schanda, B.H. Meier, M. Ernst. [*J. Magn. Reson.*](https://doi.org/10.1016/j.jmr.2011.03.015), **2011**, 210, 246-259.\n",
"\n",
"[4] K. Aebischer, L.M. Becker, P. Schanda, M. Ernst. [*Magn. Reson.*](https://doi.org/10.5194/mr-5-69-2024), **2024**, 5, 69-86\n",
"\n",
"[5] S.A. Smith, T.O. Levante, B.H. Meier, R.R. Ernst. [*J. Magn. Reson. A*](https://doi.org/10.1006/jmra.1994.1008), **1994**, 106, 75-105."
]
},
{
"cell_type": "markdown",
"id": "a00dac14",
"metadata": {},
"source": [
"## Setup"
]
}
,
{
"cell_type": "code",
"execution_count": 0,
"metadata": {"tags": [
"remove-cell"
]},
"outputs": [],
"source": [
"# SETUP SLEEPY\n",
"import sys\n",
"if 'google.colab' in sys.modules:\n",
" !pip install sleepy-nmr"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "0c25513a",
"metadata": {},
"outputs": [],
"source": [
"import SLEEPY as sl\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from time import time"
]
},
{
"cell_type": "markdown",
"id": "4c859836",
"metadata": {},
"source": [
"## Build the system\n",
"We construct a system undergoing 3-site symmetric exchange, in order to obtain a residual dipole coupling without asymmetry ($\\eta$=0). The `sl.Tools.Setup3siteSym(...)` function is used to create the three-site hopping motion out of the initial system, `ex0`."
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "62722ee4",
"metadata": {},
"outputs": [],
"source": [
"ex0=sl.ExpSys(v0H=600,Nucs=['15N','1H'],vr=60000,pwdavg=sl.PowderAvg('bcr20'),n_gamma=30)\n",
"# After varying the powder average and n_gamma\n",
"# a beta-average and 30 gamma angles were determined to be sufficient\n",
"delta=sl.Tools.dipole_coupling(.102,'15N','1H')\n",
"phi=35*np.pi/180\n",
"\n",
"ex0.set_inter('dipole',i0=0,i1=1,delta=delta)\n",
"L=sl.Tools.Setup3siteSym(ex0,tc=1e-9,phi=phi)"
]
},
{
"cell_type": "markdown",
"id": "ac44e08a",
"metadata": {},
"source": [
"## Generate and plot pulse sequences\n",
"We break the REDOR sequence up into three sequences (`L.Sequence()`), a first half (`first`) with $^1$H $\\pi$-pulses coming in the middle and end of the rotor period, a middle sequence (`center`) with a $^{15}$N $\\pi$-pulse to refocus the $^{15}$N chemical shift, and a second half (`second`) with $^1$H $\\pi$-pulses coming at the beginning and middle of the rotor period."
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "097096d8",
"metadata": {},
"outputs": [],
"source": [
"v1=120e3 #100 kHz pulse\n",
"tp=1/v1/2 #pi/2 pulse length\n",
"\n",
"t=[0,L.taur/2-tp,L.taur/2,L.taur-tp,L.taur]\n",
"first=L.Sequence().add_channel('1H',t=t,v1=[0,v1,0,v1],phase=[0,0,0,np.pi/2,0])\n",
"t=[0,tp,L.taur/2,L.taur/2+tp,L.taur]\n",
"second=L.Sequence().add_channel('1H',t=t,v1=[v1,0,v1,0],phase=[np.pi/2,0,0,0,0])\n",
"center=L.Sequence().add_channel('15N',t=[0,L.taur/2-tp/2,L.taur/2+tp/2,L.taur],\n",
" v1=[0,v1,0])\n",
"\n",
"rho=sl.Rho('15Nx','15Nx')"
]
},
{
"cell_type": "markdown",
"id": "5be0f664",
"metadata": {},
"source": [
"We plot the three sequences below."
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "5ac0da36",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEYCAYAAABRMYxdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAu40lEQVR4nO3dfZRV9X3v8fcXRBExPMiTMCweMoo8BKcKPoUYvAYhhEpNwOJiqXPRWlPtXe1tG+lSb0LUXu5N2vQBV6qminIRMuWWYigSCc2EaAgE7IiEh8JFCCCKICIGCTL87h9zhgzDDMycvff5/ea3P6+1WHNmnzPnfA7z2T9+7L3P3uacQ0RERCQm7XwHEBEREUmbJjgiIiISHU1wREREJDqa4IiIiEh0NMERERGR6JznO0CaevTo4QYOHOg7hrRx69evP+Cc6+k7Rz31WtISUrfVa0lLc72OaoIzcOBA1q1bd8byRYsWATBlypRSR5I2yMx2+c7QkHotaQmp2+q1pKW5Xkc1wWlOnz59fEcQSZ16LTFSryUtuZjgjBkzxncEkdSp1xIj9VrSooOMRUREJDq5mOBUVVVRVVXlO4ZIqtRriZF6LWnJxS6qsrIy3xFEUqdeS4zUa0lLLrbg3HDDDdxwww3MmDGDXr16MWLEiFP3feMb36Bfv35UVFRQUVHBsmXLAKiursbM+MEPfnDqsZMmTaK6urrU8UWapF5LjNRrSUsuJjj1KisrWb58+RnL//RP/5SamhpqamqYOHHiqeVlZWU88cQTpYwo0mrqtcRIvZakcrGLasGCBQDccccd7Ny5s8U/d+WVV/LJJ5+wYsUKxo0bl1E6keLE0utD36/iw6VLfcfgU5Mm0e33b/cdI/di6bX4l4stOIMGDWLQoEHN3j9nzhxGjhzJjBkzOHTo0Gn3PfLIIzz++ONZRxRptVh6/eHSpRzbssVrhmNbtgQxyZJ4ei3+5WILznXXXdfsfV/96ld59NFHMTMeffRR/uzP/oxnn3321P2f+9znAPjpT3+aeU6R1oip1x2vuIIB817w9vq77rzL22vL6WLqtfiViy04Z9O7d2/at29Pu3bt+IM/+APWrl17xmMefvhh7duVNkW9lhip19IauZjgzJ8/n/nz5zd53759+07dXrx48WlH7Ne75ZZbOHToEG+88UZmGUVaS72WGKnXkpZc7KK6/PLLgbqD1qqrqzlw4ABlZWXMmjWL6upqampqMDMGDhzIU0891eRzPPzww0yePLmUsUXOSr2WGKnXkhZzzvnOkJpRo0a5pq5OK9IaZrbeOTfKd456sfe6/viXEI7B8ZmhFELqduy9ltJprte52EUlIiIi+ZKLCc4LL7zACy/E/T8zyR/1WmKkXktacnEMzvDhw31HEEmdei0xUq8lLbmY4Fx99dW+I4ikTr2WGKnXkpZc7KISERGRfMnFBGfu3LnMnTvXdwyRVKnXEiP1WtJSsl1UZvYsMAnY75wbUVjWHfg+MBDYCdzunDtUuO8vgXuAWuC/Oed+WOxrV1RUJIkuEiT1WmKkXktaSnkMzlxgDtDw8PiZwErn3Gwzm1n4/iEzGwZMA4YDfYEfmdnlzrnaYl5YK4zESL2WGKnXkpaS7aJyzq0C3m+0eDLwfOH288DvNVi+0Dn3G+fcW8B24JpiX7u2tpba2qLmRiLBUq8lRuq1pMX3MTi9nXP7AApfexWW9wN2N3jcnsKyM5jZfWa2zszWvffee02+yLx585g3b156qUUypl5LjNRrKaVQPyZuTSxr8poSzrmngaeh7tTfTT3mqquuSi+ZSAmo1xIj9VpKyfcE510zu9Q5t8/MLgX2F5bvAfo3eFwZ8HaxLzJy5MgEEUXCpF5LjNRrSYvvXVQvAXcXbt8NLGmwfJqZXWBmg4DLgLXFvsgnn3zCJ598kiioSGjUa4mRei1pKeXHxBcAY4EeZrYH+DowG6gys3uAXwFTAZxzvzSzKmATcAJ4oNhPUAHMnz8fgMrKygTvQCQs6rXESL2WtJRsguOcu6OZu25u5vFPAE+k8dqjRp1xFXWRNk+9lhip15IW38fglMSIESN8RxBJnXotMVKvJS2+j8EpiWPHjnHs2DHfMURSpV5LjNRrSUsuJjgLFy5k4cKFvmOIpEq9lhip15KWFk9wzGylmU1stOzp9COl79prr+Xaa6/1HUNyLIv1R70W39RrCVlrjsEZRN11okY752YVlrWJo8GGDh3qO4JI6uuPei0BUK8lWK3ZRfUBdZ946m1mPzCzLtlESt/Ro0c5evSo7xiSbx+Q8vqjXksAPkC9lkC1ZoJjzrkTzrk/Av4v8Cq/vXZU0KqqqqiqqvIdQ/It9fVHvZYAqNcSrNbsovrH+hvOublmtgF4MP1I6bv++ut9RxBJff1RryUA6rUE65wTHDP7BwoXujSzv29090dZhErbkCFDfEeQnMpy/VGvxRf1WtqClmzBWdfg9izqLrHQpnz0Ud361rlzZ89JJIcyW3/Ua/FIvZbgnXOC45x7vv62mf1Jw+/bikWLFgG6tomUXpbrj3otvqjX0ha09lINLpMUGRszZozvCCKQ8vqjXksg1GsJUi6uRVVeXu47gkjq1GuJkXotaWnJQcZHqJuhG3ChmX1YfxfgnHOfyjBfKg4fPgxAly5t5tQ9Eoks1x/1WnxRr6UtaMkxOBeXIkiWFi9eDGifrpReluuPei2+qNfSFrR4F5WZXe2cW99o2e86536Qfqx03Xjjjb4jSM5lsf6o1+Kbei0ha82ZjJ8xs8/Uf2NmdwCPpB8pfYMHD2bw4MG+Y0i+pb7+qNcSAPVagtWag4ynAIvMbDowBrgLuCWTVCk7dOgQAN26dfOcRHIs9fVHvZYAqNcSrBZPcJxzO8xsGvCvwG7gFufcx1kFS9OSJUsA7dMVf7JYf9Rr8U29lpC15FNUb3L6eQ66A+2BNWaGc25kVuHSMnbsWN8RJKeyXH/Ua/FFvZa2oCVbcCZlniJjAwcO9B1B8iuz9Ue9Fo/UawleSz4mvqsUQbJ04MABAHr06OE5ieRNluuPei2+qNfSFrTmU1Rt1tKlS1m6dKnvGCKpUq8lRuq1pCUXl2q4+eabfUcQSZ16LTFSryUtLTnI+Hrg5865NnmhTYD+/fv7jiA5leX6o16LL+q1tAUt2UV1N7DezBaaWaWZ9ck6VNr279/P/v37fceQfMps/VGvxSP1WoLXkoOM7wcwsyuALwJzzawL8GNgOfCac64205QJLVu2DNB5FaT0slx/1GvxRb2WtqA1J/rbAmwBvmNmFwI3AVOBvwFGZRMvHePGjfMdQXIui/VHvRbf1GsJWVEHGRfOVLms8Cd4/fr18x1B5JS01h/1WkKiXktogvgUlZntBI4AtcAJ59woM+sOfB8YCOwEbnfOHSrm+d955x0A+vRpc4cPiTRLvZYYqdeSlpDOg3OTc67COVe/WXMmsNI5dxmwsvB9UZYvX87y5cvTyCgSDPVaYqReS1qC2ILTjMnA2MLt54Fq4KFinmjChAnpJBIJiHotMVKvJS2Jt+CYWVGTjkYc8IqZrTez+wrLejvn9gEUvvYq9sn79OnTos2dL7/8Mg8//DAnT54s9qVEWiXJ+qNeS6jUawlBqyc4ZlbV4M8/A/emkOOzzrmrqPu44QNmdmMr8txnZuvMbN17773X5GP27t3L3r17gbrNn0OGDKG8vJzZs2ef9rhVq1YxevRoVq9eXfQbETmblq4/6rW0Jeq1hKiYLTgfOuduL/yZCvwoaQjn3NuFr/uBxcA1wLtmdilA4WuTZ35yzj3tnBvlnBvVs2fPJp9/xYoVrFixgtraWh544AFefvllNm3axIIFC9i0adOpx7Vv35758+czdOjQpG9JpDktWn/Ua2lj1GsJTouPwTGznwGPAI83uuvhJAHM7CKgnXPuSOH2LcA3gZeoO1vm7MLXJcW+xsSJEwFYu3Yt5eXlDB48GIBp06axZMkShg0bBsDjjzd+ayLpyGL9Ua/FN/VaQtaaLTj3AQ8A/1S4DgkAzrn3E2boDbxqZm8Aa4F/c84tp25iM87MtgHjCt8XpVevXvTq1Yu9e/eedp2TsrKyU5tCRTKW+vqjXksA1GsJVmvOZLwR+IqZXQV808wAHnHO1SQJ4JzbAVzZxPKDQCqXld29e3f9c55xX+F9iGQqi/VHvRbf1GsJWTEfE98OPAb8V2Bdkc9RUitXrgRgyJAhp1YegD179tC3b19fsSSfUlt/1GsJiHotwWnNMTj/DlwGHAM2Ff5UZhMrXZMmTQKga9eubNu2jbfeeot+/fqxcOFCXnzxRc/pJA+yWH/Ua/FNvZaQtWaW/efA5sL1RtqUHj16nLo9Z84cxo8fT21tLTNmzGD48OEek0mOpL7+qNcSAPVagtWaY3BezzJIlnbu3AnAwIEDmThx4qmj9EVKJYv1R70W39RrCVlI16LKTHV1NdXV1b5jiKRKvZYYqdeSluAPEE7D5MmTfUcQSZ16LTFSryUtuZjgdOvWzXcEkdSp1xIj9VrSkotdVDt27GDHjh2+Y4ikSr2WGKnXkpZcbMFZtWoVwKlTfovEQL2WGKnXkpZcTHBuu+023xFEUqdeS4zUa0lLLiY4Xbp08R1BJHXqtcRIvZa05OIYnO3bt7N9+3bfMURSpV5LjNRrSUsutuC8+uqrAJSXl3tOIpIe9VpipF5LWnIxwZkyZYrvCCKpU68lRuq1pCUXE5zOnTv7jiCSOvVaYqReS1pycQzO1q1b2bp1q+8YIqlSryVG6rWkJRcTnNWrV7N69WqWL1/OkCFDKC8vZ/bs2U0+trKykkWLFp2xfOzYsaxbt+7U9zt37mTEiBGZZRY5F/VaYqReS1pysYvq9ttvp7a2loqKClasWEFZWRmjR4/m1ltvZdiwYb7jiRRFvZYYqdeSllxswenUqRMbN26kvLycwYMHc/755zNt2jSWLFly1p979NFHqays5OTJk2d93L333ktFRQUVFRX07NmTWbNmpRlfpEnqtcRIvZa05GILzubNm1mzZg39+/c/taysrIw1a9Y0+zNf+9rXOHz4MM899xxmBsD06dO58MILATh+/Djt2tXND7/3ve8BsGvXLsaPH09lZWVG70Tkt9RriZF6LWnJxRacNWvWsG3btjOW168IjT322GN88MEHPPXUU6c9Zv78+dTU1FBTU8OyZctO+5ljx44xdepU5syZw4ABA9J9AyJNUK8lRuq1pCUXE5xp06YxdepUdu/efWrZnj176Nu3b5OPHz16NOvXr+f9999v8Wvcf//9fPnLX+YLX/hC4rwiLaFeS4zUa0lLLiY4HTt2ZMyYMWzbto233nqL48ePs3DhQm699dYmHz9hwgRmzpzJl770JY4cOXLO53/yySc5cuQIM2fOTDu6SLPUa4mRei1pycUxOBs3bgRgzpw5jB8/ntraWmbMmMHw4cOb/ZmpU6dy5MgRbr311jM2bzb27W9/mw4dOlBRUQHU/e/g/vvvTy2/SFPUa4mRei1pMeec7wypGTVqlGt47oN6c+fOBdDBZNIiZrbeOTfKd456sfd61513ATBg3gu5zlAKIXU79l5L6TTX61xswZk+fbrvCCKpU68lRuq1pCUXE5wOHTr4jiCSOvVaYqReS1pycZDxhg0b2LBhg+8YIqlSryVG6rWkJRdbcF5//XUARo4c6TmJSHrUa4mRei1pycUE58477/QdQSR16rXESL2WtAS9i8rMJpjZVjPbbmZFn7Sgffv2tG/fPs1oIt6p1xIj9VrSEuwEx8zaA08CXwSGAXeYWVGXkq0/XbdITNRriZF6LWkJeRfVNcB259wOADNbCEwGNrX2iaqee46PPv6YtUteSjmitGWfKf80z7z2mu8YRYul1yePHqVdp050HDvWW4Zjm7fU5ejdx1uGtKjXEqNieh3sif7MbAowwTl3b+H7O4FrnXMPNnrcfcB9hW+HAFubecoewIGM4qZB+YqXdrYBzrmeKT5fq6nXJZO3fF67rV6XTN7yNdnrkCc4U4HxjSY41zjn/rjI51sXyhk8m6J8xQs5W9ZCf+/Kl0zo+bIS+vtWvmRKlS/YY3CAPUD/Bt+XAW97yiIiIiJtSMgTnF8Al5nZIDM7H5gGaKesiIiInFOwBxk7506Y2YPAD4H2wLPOuV8meMqn00mWGeUrXsjZshb6e1e+ZELPl5XQ37fyJVOSfMEegyMiIiJSrJB3UYmIiIgURRMcERERiU70E5y0LveQFTPbaWZvmlmNma0LIM+zZrbfzDY2WNbdzFaY2bbC126B5fuGme0t/B3WmNlEX/lKRb1udR71ug1Qr1udR70+i6gnOGle7iFjNznnKgI5b8FcYEKjZTOBlc65y4CVhe99mcuZ+QC+U/g7rHDOLStxppJSr4syF/U6aOp1UeaiXjcr6gkODS734Jw7DtRf7kGa4ZxbBbzfaPFk4PnC7eeB3ytlpoaayZc36nUrqddtgnrdSur12cU+wekH7G7w/Z7CspA44BUzW184jXmIejvn9gEUvvbynKcpD5rZhsImUW+bZEtEvU6Heh0W9Tod6nVB7BMca2JZaJ+L/6xz7irqNss+YGY3+g7UBn0X+DRQAewD/tprmuyp1/mgXqvXMSpZr2Of4AR/uQfn3NuFr/uBxdRtpg3Nu2Z2KUDh637PeU7jnHvXOVfrnDsJPEOYf4dpUq/ToV6HRb1Oh3pdEPsEJ+jLPZjZRWZ2cf1t4BZg49l/youXgLsLt+8GlnjMcob6lbngNsL8O0yTep0O9Tos6nU61OuCYC/VkIYMLveQtt7AYjODut/Fi8655T4DmdkCYCzQw8z2AF8HZgNVZnYP8CtgamD5xppZBXWbs3cCf+grXymo162nXodPvW499focr69LNYiIiEhsYt9FJSIiIjmkCY6IiIhERxMcERERiY4mOCIiIhIdTXBEREQkOprgiIiISHQ0wREREZHoaIIjIiIi0dEER0RERKKjCY6IiIhERxMcERERiU5UF9vs0aOHGzhwoO8Y0satX7/+gHOup+8c9dRrSUtI3VavJS3N9TqqCc7AgQNZt27dGcsXLVoEwJQpU0odSZoQwu/jbBnMbFep85yNei1pCanbIfc6hAzScs31OqoJTnP69OnjO4I0EMLvI4QMScXwHkQaC6HXIWSQ5HIxwRkzZozvCNJACL+PEDIkFcN7EGkshF6HkEGS00HGIiIiEp1cTHCqqqqoqqryHUMKQvh9hJAhqRjeg0hjIfQ6hAySXC52UZWVlfmOIA2E8PsIIUNSMbwHkcZC6HUIGSS5XGzBueGGG7jhhhuYMWMGvXr1YsSIEafuq6ysPHXEfL3OnTuXOmKu1P8+8p4hKfVaYhRCr2MYHyQnW3DqVVZW8uCDD3LXXXeV/LUPfb+KD5cuLfnrNvapSZPo9vu3e82gv4t0+ey1SNre+au/4jebtzDh3Xf5yqjR/PefvcauO+u6/dHPXuO9HW+xa8lLpx7vfvObU/enyff4EMo4Cf7/LoqViy04CxYsYMGCBdx44410797dS4YPly7l2JYtXl673rEtW4JYYRavX8crnrcmvNK5M4vXn3kOjrYkhF6LpG35iROsvLQP1/buTZcLLvCSIYTxIYR/MyCcfzeKkYstOIMGDTrr/X/xF3/B448/nnmOjldcwYB5L2T+Os3J4n85xbj06MfQvr3Xv4uyB/8Yjn7s7fXTEEqvRdJ0xS23ADDguutwO3fSYdKkU2NF58pK/vdPfsJTHx4+9Xi74ILUx5JQxgff/2ZAOP9uFCMXE5zrrrvurPd/61vfOu2MlTpWIVvDDh8+94NykCEp9VpiFEKvYxgfJCe7qERERCRfcrEFZ/78+QBMnz7dcxIBWHHppQDcm/MMSanXEqMQeh3D+CA5meBcfvnlANxxxx1UV1dz4MABysrKmDVrludk+dT/17/2HSGIDEmp1xKjEHodw/ggOZngjB49Gqj71Elj99xzzxnLPvroo8wz5dkVH37oO0IQGZJSryVGIfQ6hvFBdAyOiIiIRCgXW3BeeKHuY3Y6EVoYfti3bv/2fTnPkJR6LTEKodcxjA+SkwnO8OHDfUeQBgYFsKskhAxJqdcSoxB6HcP4IDmZ4Fx99dW+I0gDl394xHeEIDIkpV5LjELodQzjg+gYHBEREYlQySY4Zvasme03s40NlnU3sxVmtq3wtVuD+/7SzLab2VYzG5/ktefOncvcuXOTPIWk6OV+fXm5X9/cZ0hKvZYYhdDrGMYHKe0uqrnAHKDhhTVmAiudc7PNbGbh+4fMbBgwDRgO9AV+ZGaXO+dqi3nhioqKJLklZeUBbP4NIUNS6rXEKIRexzA+SAknOM65VWY2sNHiycDYwu3ngWrgocLyhc653wBvmdl24BpgdTGvHcIKI7912RH/g0cIGZJSryVGIfQ6hvFB/B+D09s5tw+g8LVXYXk/YHeDx+0pLCtKbW0ttbVFbfyRDJws/Ml7hqTUa4lRCL2OYXwQ/xOc5lgTy1yTDzS7z8zWmdm69957r8knmzdvHvPmzUsznyTww359+aHn/dshZDgb9Vpi1FZ6Hfr4IC3j+2Pi75rZpc65fWZ2KbC/sHwP0L/B48qAt5t6Aufc08DTAKNGjWpyEnTVVVell1gSu/yw/9Ogh5DhbNRriVFb6XXo44O0jO8JzkvA3cDswtclDZa/aGZ/Q91BxpcBa4t9kZEjRyaMKWn6dAAn0QohQ1LqtcQohF7HMD5ICSc4ZraAugOKe5jZHuDr1E1sqszsHuBXwFQA59wvzawK2AScAB4o9hNUAJ988gkAHTp0SPIWJCUnrKk9kPnLkJR6LTEKodcxjA9S2k9R3dHMXTc38/gngCfSeO358+cDUFlZmcbTSUIrCtd5+XTOMySlXkuMQuh1DOOD+N9FVRKjRo3yHUEauOLwYd8RgsiQlHotMQqh1zGMD5KTCc6IESN8R5AGBn30a98RgsiQlHotMQqh1zGMD5KTCc6xY8cA6Nixo+ckAnC8nf+zE4SQISn1WmIUQq9jGB8kJxOchQsXAjpWIRQrL+0D1H00Ls8ZklKvJUYh9DqG8UFaMcExs5XAXzvnljVY9rRz7r5MkqXo2muv9R1BGhj2gf/926XOkMX6o16Lb7H2OoQxSpJrzRacQdRdCHO0c25WYZn/o8FaYOjQob4jSAMDfu1//7aHDKmvP+q1BCDKXocwRklyrdnR+AF1H+nubWY/MLMu2URK39GjRzl69KjvGFJwrF07jnnex+0hwwekvP6o1xKAD4iw1yGMUZJca7bgmHPuBPBHZlYJvAp0yyRVyqqqqgAdqxCKHxf2bw/JV4bU1x/1WgIQZa9DGKMkudZMcP6x/oZzbq6ZbQAeTD9S+q6//nrfEaSBEYc+8B3BR4bU1x/1WgIQZa9DGKMkuXNOcMzsHyhcydvM/r7R3W3igh1DhmgeHpL+AexWKVWGLNcf9Vp8ib3XIYxRklxLtuCsa3B7FnXXkGpTPipcOK1z586ekwjA0fbtfUcoZYbM1h/1WjyKutchjFGS3DknOM655+tvm9mfNPy+rVi0aBGgYxVC8ZM+vQHw+VmJUmXIcv1Rr8WX2HsdwhglybX2RH8ukxQZGzNmjO8I0sBnDh3yHcFXhlTXH/VaAhFdr0MYoyS5XJzJuLy83HcEaaDs6Me+IwSRISn1WmIUQq9jGB+kZQcZH6Fuhm7AhWb2Yf1dgHPOfSrDfKk4XLgybJcubebUPVH79Xn+92+XKkOW6496Lb7E3usQxihJriXH4FxciiBZWrx4MaBjFUKxqnfd/u1hOciQ5fqjXosvsfc6hDFKkmvNtaiuds6tb7Tsd51zP0g/VrpuvPFG3xGkgSvf979/u9QZslh/1GvxLdZehzBGSXKtORf1M2b2mfpvzOwO4JH0I6Vv8ODBDB482HcMKej78cf0/djvPm4PGVJff9RrCUCUvQ5hjJLkWnOQ8RRgkZlNB8YAdwG3ZJIqZYcKR8R369YmriwRvSPn+T+23UOG1Ncf9VoCEGWvQxijJLkW/xadczvMbBrwr8Bu4BbnXJuY4i5ZsgTQsQqheLV3LwBG5ChDFuuPei2+xdrrEMYoSa4ln6J6k9PPc9AdaA+sMTOccyOzCpeWsWPH+o4gDfzOwfd9RyhZhizXH/VafIm91yGMUZJcS7bgTMo8RcYGDhzoO4I00OfYMd8RSpkhs/VHvRaPou51CGOUJNeSj4nvKkWQLB04cACAHj16eE4iAIc7dPAdoWQZslx/1GvxJfZehzBGSXKt+RRVm7V06VKWLl3qO4YU/KxXT37Wq2fuMySlXkuMQuh1DOOD5ORSDTfffLPvCNLA1QcP+o4QRIak1GuJUQi9jmF8kJYdZHw98HPnXJu80CZA//79fUeQBnod+43vCCXLkOX6o16LL7H3OoQxSpJryS6qu4H1ZrbQzCrNrE/WodK2f/9+9u/f7zuGFBw6/3wOnX9+XjJktv6o1+JR1L0OYYyS5FpykPH9AGZ2BfBFYK6ZdQF+DCwHXnPO1WaaMqFly5YBOl9IKH7es+7gwYocZMhy/VGvxZfYex3CGCXJteZEf1uALcB3zOxC4CZgKvA3wKhs4qVj3LhxviNIA6MP+N+/XeoMWaw/6rX4FmuvQxijJLmiDjIunKlyWeFP8Pr16+c7gjTQ4zf+92/7zJDW+qNeS0hi6nUIY5QkF8SnqMxsJ3AEqAVOOOdGmVl34PvAQGAncLtzrqhLvL7zzjsA9OnT5g4fitLBwr7tATnPkJR6LTEKodcxjA8S1nlwbnLOVTjn6jdrzgRWOucuA1YWvi/K8uXLWb58eRoZJQVre/ZgbU+/J6cLIUNS6rXEKIRexzA+SCBbcJoxGRhbuP08UA08VMwTTZgwIZ1Ekopr3jvgO0IQGZJSryVGIfQ6hvFBUtiCY2ZFTToaccArZrbezO4rLOvtnNsHUPjaq5nXv8/M1pnZuvfee6/JJ+/Tp0+LNne+/PLLPPzww5w8ebKoNyEtc8nx41xy/Hii50j6u2ppBjP7opk9YWaZbO1sbv1Rr6UtK0Wvz6ZU44PWrbC1etA2s6oGf/4ZuDeFHJ91zl1F3ccNHzCzG1v6g865p51zo5xzo3r2bPrU2nv37mXv3r1A3ebPIUOGUF5ezuzZs0973KpVqxg9ejSrV68u+o3IuR244AIOXHABAH/3d3/HiBEjGD58OH/7t3976jFn+z1B8t9VwwzneK0bgV8A1xf1Qo20dP1Rr6UtyarXoY8PWrfCVsz/Sj90zt1e+DMV+FHSEM65twtf9wOLgWuAd83sUoDC16LP/LRixQpWrFhBbW0tDzzwAC+//DKbNm1iwYIFbNq06dTj2rdvz/z58xk6dGii9yNn94sel/CLHpewceNGnnnmGdauXcsbb7zB0qVL2bZt2zl/T5D8d1WfoanXAjo2eGgtMB3YXOTbbSy19Ue9loCk3utQxwetW21Hi4/BMbOfAY8Ajze66+EkAczsIqCdc+5I4fYtwDeBl6g7W+bswtclxb7GxIkTAVi7di3l5eUMHjwYgGnTprFkyRKGDRsGwOOPN35rkoXrCvu3N2/ezHXXXUenTp0A+PznP8/ixYv53Oc+d9bfEyT/XdVnaKoTGzZs6Fr/OOfcI4leqCCL9Ue9Ft+y7PVPfvKTIMcHrVttR2u24NwHPAD8U+E6JAA4595PmKE38KqZvQGsBf7NObecuonNODPbBowrfF+UXr160atXL/bu3XvadU7KyspObQqV0ul2/Djdjh9nxIgRrFq1ioMHD3L06FGWLVvG7t27S/J7qs/Q1GsBWZyjPfX1R72WAGTW61DHB61bbUdrzmS8EfiKmV0FfNPMAB5xztUkCeCc2wFc2cTyg0Aql5XdvXt3/XOecV/hfUgJ7e9Yt2979NChPPTQQ4wbN47OnTtz5ZVXct5555Xk91SfoZlrBaZ+AcEs1h/1WnzLstdDAx0ftG61HcUcg7MdeAzYA6xLN042Vq5cycqVKykrKzu18gDs2bOHvn37ekyWT+svuYT1l1wCwD333MPrr7/OqlWr6N69O5dddllJfk/1GZp6LeCTVF/sdKmtP+q1BCT1XkOY44PWrbajNcfg/DtwGXAM2FT4U5lNrHRNmjQJgK5du7Jt2zbeeust+vXrx8KFC3nxxRc9p8ufG/b/9uOh+/fvp1evXvzqV7/iX/7lX1i9ejUXX3xx5r+n+gzDRo8+47WAD1J9MbJZf9Rr8S3LXkOY44PWrbajNSf6+3Ngc+F6I21Kjx6/PSPlnDlzGD9+PLW1tcyYMYPhw4d7TJZPXT757QaSr3zlKxw8eJAOHTrw5JNP0q1bNyD731N9hvPOO++M19qwYcOxVF+sTurrj3otAci01yGOD1q32o7WHIPzepZBsrRz504ABg4cyMSJE08dpS9+vNOx7lPYA4Cf/vSnTT4m699TwwyNX+uRR1L54NRpslh/1GvxLetehzg+SNsR0rWoMlNdXU11dbXvGFLwH5d05z8u6Z77DEmp1xKjEHodw/ggYV+LKjWTJ0/2HUEaGPNu0edsjCpDUuq1xCiEXscwPkhOJjj1+20lDBefOOE7QhAZklKvJUYh9DqG8UFyMsHZsWMHwKmzUYpfb194IVC3fzvPGZJSryVGIfQ6hvFBcjLBWbVqFaB/CELxRve6/6GlcvXKNpwhKfVaYhRCr2MYHyQnE5zbbrvNdwRp4MZ33/UdIYgMSanXEqMQeh3D+CA5meB06dLFdwRp4KITtb4jBJEhKfVaYhRCr2MYHyQnE5zt27cDUF5e7jmJAOzp5H//dggZklKvJUYh9DqG8UFyMsF59dVXAf1DEIo3C5+S+GzOMySlXkuMQuh1DOOD5GSCM2XKFN8RpIHPv+N//3YIGZJSryVGIfQ6hvFBcjLB6dy5s+8I0kCnWv/7t0PIkJR6LTEKodcxjA+Sk0s1bN26la1bt7J8+XKGDBlCeXk5s2fPbvKxlZWVLFq06IzlY8eOZd26dae+37lzJyNGjMgsc8x2d+rE7k6dcp8hKfVaYhRCr2MYHyQnW3BWr17NyZMneeKJJ1ixYgVlZWWMHj2aW2+9lWHDhvmOlzsbu3UFYEzOMySlXkuMQuh1DOOD5GQLzu23386gQYMoLy9n8ODBnH/++UybNo0lS5ac9eceffRRKisrOXny5Fkfd++991JRUUFFRQU9e/Zk1qxZacaPzk373uGmfe/kPkNS6rXEKIRexzA+SE624HTq1ImDBw/Sv3//U8vKyspYs2ZNsz/zta99jcOHD/Pcc89hZgBMnz6dCwun8D5+/Djt2tXND7/3ve8BsGvXLsaPH09lZWVG7yQOHc8xAOUlQ1LqtcQohF7HMD5ITrbgbN68md27d5+xvH5FaOyxxx7jgw8+4KmnnjrtMfPnz6empoaamhqWLVt22s8cO3aMqVOnMmfOHAYM0NkTzmbXRRex66KLcp8hKfVaYhRCr2MYHyQnE5w1a9Zw8ODB01aaPXv20Ldv3yYfP3r0aNavX8/777/f4te4//77+fKXv8wXvvCFxHljt6lrFzZ19Xu20hAyJKVeS4xC6HUM44PkZBfVtGnTOHHiBCNHjuStt96iX79+LFy4kBdffLHJx0+YMIHx48fzpS99iVdeeYWLL774rM//5JNPcuTIEWbOnJlF/OjcHMC+7RAyJKVeS4xC6HUM44PkZILTsWNHAObMmcP48eOpra1lxowZDB8+vNmfmTp1KkeOHOHWW289Y/NmY9/+9rfp0KEDFRUVQN3/Du6///7U8sfm/AD2b4eQISn1WmIUQq9jGB8EzDnnO0NqRo0a5Rqe+6Dexo0bAbye32PXnXcBMGDeC7nOAFD91a8CMPa73w0yg5mtd86NKnWm5oTca2lbQup2yL0OYYwKZbwOJcfZNNfrXGzBqV+J9A9BGLYUrhY8NucZklKvJUYh9DqG8UFyMsGZPn267wjSwLi39/mOEESGpNRriVEIvY5hfJCcTHA6dOjgO4I0cF4Au0VDyJCUei0xCqHXMYwPkpMJzoYNGwAYOXKk5yQC8P8KF9PzeVaVEDIkpV5LjELodQzjg+RkgvP6668D+ocgFP/Z5VMA/JecZ0hKvZYYhdDrGMYHyckE58477/QdQRoYv/dt3xGCyJCUei0xCqHXMYwPEviZjM1sgpltNbPtZlb02cbat29P+/bt04wmCbTDf/FCyJCUei0xCqHXMYwPEvAWHDNrDzwJjAP2AL8ws5ecc5ta+1w1NTUAp07sJH5tK5xp1Of+7RAyJKVeS4xC6HUM44MEPMEBrgG2O+d2AJjZQmAy0OoJTtVzz/HRxx+zdslLKUdsuZNHj9KuUyc6jh3rLcOxzVvqcvTu4y0DwDW/OwnateNxj38XIy+/nIs6dKAtX2EphF5LeD5T/mmeee013zGKtmjRIj766KNTEx0fRl5+OZw8ySMex8oQ/s2AcP7dKKbXwZ7J2MymABOcc/cWvr8TuNY592Cjx90H3Ff4dgiwtZmn7AEcyChuGpSveGlnG+Cc65ni87Wael0yecvntdvqdcnkLV+TvQ55gjMVGN9ognONc+6Pi3y+daGcorwpyle8kLNlLfT3rnzJhJ4vK6G/b+VLplT5Qj6Oag/Qv8H3ZYAObRcREZFzCnmC8wvgMjMbZGbnA9MAHWwgIiIi5xTsQcbOuRNm9iDwQ6A98Kxz7pcJnvLpdJJlRvmKF3K2rIX+3pUvmdDzZSX09618yZQkX7DH4IiIiIgUK+RdVCIiIiJF0QRHREREohP9BCetyz1kxcx2mtmbZlZjZusCyPOsme03s40NlnU3sxVmtq3wtVtg+b5hZnsLf4c1ZjbRV75SUa9bnUe9bgPU61bnUa/PIuoJToPLPXwRGAbcYWbD/KZq0k3OuYpAzlswF5jQaNlMYKVz7jJgZeF7X+ZyZj6A7xT+Diucc8tKnKmk1OuizEW9Dpp6XZS5qNfNinqCQ4PLPTjnjgP1l3uQZjjnVgHvN1o8GXi+cPt54PdKmamhZvLljXrdSup1m6Bet5J6fXaxT3D6AbsbfL+nsCwkDnjFzNYXTmMeot7OuX0Aha+9POdpyoNmtqGwSdTbJtkSUa/ToV6HRb1Oh3pdEPsEx5pYFtrn4j/rnLuKus2yD5jZjb4DtUHfBT4NVAD7gL/2miZ76nU+qNfqdYxK1uvYJzjBX+7BOfd24et+YDF1m2lD866ZXQpQ+Lrfc57TOOfedc7VOudOAs8Q5t9hmtTrdKjXYVGv06FeF8Q+wQn6cg9mdpGZXVx/G7gF2Hj2n/LiJeDuwu27gSUes5yhfmUuuI0w/w7TpF6nQ70Oi3qdDvW6INhLNaQhg8s9pK03sNjMoO538aJzbrnPQGa2ABgL9DCzPcDXgdlAlZndA/wKmBpYvrFmVkHd5uydwB/6ylcK6nXrqdfhU69bT70+x+vrUg0iIiISm9h3UYmIiEgOaYIjIiIi0dEER0RERKKjCY6IiIhERxMcERERiY4mOCIiIhIdTXBEREQkOprgtHFm1tXM/ugcj3nKzD5bqkwiSanXEiP1urQ0wWn7ugJnXWGAa4GfZx9FJDVdUa8lPl1Rr0tGE5y2bzbwaTOrMbNvNb7TzIYC/+mcq220vNrMhhRuX2JmGwu3LzKzfzOzN8xso5n9finehEgj6rXESL0uoaivRZUTM4ERzrmKZu7/ItDU9VLKgW2F2yOBNwu3JwBvO+e+BGBmXdKLKtJi6rXESL0uIW3Bid94Gq0wZjYA2Fu4XD3UrTAbCrffBL5gZv/LzD7nnDtcuqgiLaZeS4zU6xRpghMxM+sEdHXOvd3orgp+u4IAXF3/vXPuPwvfvwn8TzP7HyWIKtJi6rXESL1OnyY4bd8R4OJm7rsJ+HETy68EOgKY2WXAZAqbPM2sL3DUOfd/gG8DV6UdWKQF1GuJkXpdQprgtHHOuYPAa4UDzBoftNbc/twKoJ2ZvQH8D2AzcHfhvs8Aa82sBngYeDyL3CJno15LjNTr0jLnnO8MkhEzex241jn3SaPl24Hfcc4d8ZNMpHjqtcRIvU6fPkUVMefcGZsrzexi4KRWFmmr1GuJkXqdPm3BERERkejoGBwRERGJjiY4IiIiEh1NcERERCQ6muCIiIhIdDTBERERkehogiMiIiLR0QRHREREovP/AUiT66JdC/AsAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig,ax=plt.subplots(2,3,figsize=[8,4],sharey=True)\n",
"first.plot(ax=ax.T[0])\n",
"center.plot(ax=ax.T[1])\n",
"second.plot(ax=ax.T[2])\n",
"fig.tight_layout()"
]
},
{
"cell_type": "markdown",
"id": "5d145dee",
"metadata": {},
"source": [
"## Propagation"
]
},
{
"cell_type": "markdown",
"id": "d71bca3a",
"metadata": {},
"source": [
"We first generate a propagator for each of the three sequences (`Ucenter`,`Ufirst`,`Usecond`). Then below, we create two propagators for the first (`U1`) and second half (`U2`) of the sequence which are increased in length by one rotor period at every step of the rotor period. "
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "4fc71d91",
"metadata": {},
"outputs": [],
"source": [
"Ucenter=center.U()\n",
"Ufirst=first.U()\n",
"Usecond=second.U()"
]
},
{
"cell_type": "markdown",
"id": "0afcc7bd",
"metadata": {},
"source": [
"At each step below, the density matrix (`rho`), is set back to the initial time via `rho.reset()`, and then multiplied with propagators for the first half (`U1`), the center (`Ucenter`), and for the second half (`U2`), where `U1` and `U2` are each increased by one rotor period at every step. `(U2*Ucenter*U1*rho)()` propagates `rho` with the full REDOR sequence, followed by detection (don't forget the parenthesis!). "
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "8cb98c14",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"25.295255184173584\n"
]
}
],
"source": [
"rho=sl.Rho('15Nx','15Nx')\n",
"U1=L.Ueye()\n",
"U2=L.Ueye()\n",
"\n",
"t0=time()\n",
"for k in range(48):\n",
" rho.reset()\n",
" (U2*Ucenter*U1*rho)()\n",
" U1=Ufirst*U1\n",
" U2=Usecond*U2\n",
"print(time()-t0)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "f4bb1b96",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAo5ElEQVR4nO3deXxU9b3/8dcnO4EkQDZICATCJoiI7CAI9iqgrda2uLZWq7W4tb3to9X21/a2t/fe1vbWtlatW7W1tSp1r3WtgrLIEkD2LQRJQlgSIBASsn9/f2T0pgnLZJk5M5P38/HIw8zMceYtcvLO+X7P+R5zziEiItJSlNcBREQk9KgcRESkDZWDiIi0oXIQEZE2VA4iItJGjNcBukpaWprLzc31OoaISFhZs2ZNuXMuvfXzEVMOubm55Ofnex1DRCSsmNmekz2vYSUREWlD5SAiIm2oHEREpI2gl4OZPW5mB81s0yleNzO7z8wKzGyDmZ0X7IwiIt2dF0cOfwTmnub1ecAw39ctwO+DkElERFoIejk4594HDp9mk8uBJ12zFUBvM+sfnHQiIgKhOeeQDRS3eFzie64NM7vFzPLNLL+srCwo4UREuoNQLAc7yXMnXVfcOfeIc26Cc25Cenqbazj88ubm/Ty3pqRD/66ISKQKxYvgSoCcFo8HAKWB+CDnHE+tLGJF4SHO6p/E6KyUQHyMiEjYCcUjh1eA631nLU0Bjjrn9gXig8yMX185lr6Jcdz+1Foqa+oD8TEiImHHi1NZnwY+AEaYWYmZ3WRmC8xsgW+T14BCoAB4FLgtkHlSe8Xzu2vHUXzkBHc/vxHdGU9ExINhJefcNWd43QG3BykOABNz+/KdOSP4+evbmPRBX748LTeYHy8iEnJCcVjJE7fMGMKFIzP4r39sYX1xhddxREQ8pXLwiYoyfjV/LBlJCdz+17Ucrdb8g4h0XyqHFvr0jOP+a8dx4FgN3/7bes0/iEi3pXJoZdzAPtw97yz+ufUAjy3Z7XUcERFPqBxO4ivTc5kzOpN73tjGtv3HvI4jIhJ0KoeTMDPu+fw5JMRG8+CiXV7HEREJOpXDKfROjOO6yQN5dUMpew5VeR1HRCSoVA6n8ZXzBxMTFcUj7xd6HUVEJKhUDqeRmZzA58cP4G9rSjhYWeN1HBGRoFE5nMHXZg6hobGJPyzVmUsi0n2oHM4gN60nl4zpz1Mrijh6QhfGiUj3oHLww62z8jhe28BfVuzxOoqISFCoHPwwOiuFC4an8/jS3dTUN3odR0Qk4FQOfrptVh6HqupYmF985o1FRMKcysFPkwb35byBvXn4vULqG5u8jiMiElAqBz+ZGbfNGsreihO8uiEgdy0VEQkZKod2uHBkBsMze/H7xbtoatKKrSISuVQO7RAVZdw6K48dB47z7raDXscREQkYlUM7feacLAb06cHjy3RRnIhELpVDO8VER/G5cdmsKDxE+fFar+OIiASEyqEDLjmnP00O3ti03+soIiIBoXLogBGZSeSl9+QfG/Z5HUVEJCBUDh1gZlw6pj8rdx+irFJDSyISeVQOHXTpOVnNQ0ubNbQkIpFH5dBBwzN7kZfek9c0tCQiEUjl0EFmxqXnZGloSUQiksqhEy4d019DSyISkVQOnaChJRGJVCqHTtDQkohEqqCXg5nNNbPtZlZgZnef5PUUM/u7ma03s81mdmOwM7aHhpZEJBIFtRzMLBp4AJgHjAKuMbNRrTa7HdjinBsLzAJ+ZWZxwczZHsMzezE0oxf/0DLeIhJBgn3kMAkocM4VOufqgGeAy1tt44AkMzOgF3AYaAhuTP+ZGZeM6c+q3Yc1tCQiESPY5ZANtLzPZonvuZbuB84CSoGNwDeccyF967VPn6OhJRGJLMEuBzvJc63vmjMH+BDIAs4F7jez5JO+mdktZpZvZvllZWVdmbNdhmcmaWhJRCJKsMuhBMhp8XgAzUcILd0IvOCaFQC7gZEnezPn3CPOuQnOuQnp6ekBCeyvS31DSwcrazzNISLSFYJdDquBYWY22DfJfDXwSqttioBPAZhZJjACKAxqyg641De09KaW8RaRCBDUcnDONQB3AG8CW4GFzrnNZrbAzBb4NvspMM3MNgLvAHc558qDmbMjPhla2qgL4kQk/MUE+wOdc68Br7V67qEW35cCFwc7V1e4dEx/fvfuTsoqa0lPivc6johIh+kK6S409+x+NDlYtP2g11FERDpF5dCFRvZLol9yAou2qRxEJLypHLqQmTF7ZDpLdpZT3xjSl2aIiJyWyqGLzR6RwfHaBvI/OuJ1FBGRDlM5dLHpQ9OIjTYWa95BRMKYyqGL9YyPYfLgVN7VvIOIhDGVQwDMHpnBzoPHKT5c7XUUEZEOUTkEwOwRzUt5LN7h3XpPIiKdoXIIgMFpPRmUmqhTWkUkbKkcAsDMmD0ig+W7yqmpb/Q6johIu6kcAmT2yAxq6ptYUXjI6ygiIu2mcgiQyYP7khAbpaElEQlLKocASYiNZnpeGou2l+Fc6/sZiYiENpVDAM0emUHR4WoKy6u8jiIi0i4qhwCa5TulVUNLIhJuVA4BNKBPIsMze2kJbxEJOyqHAJs9IoNVuw9zvLbB6ygiIn5TOQTY7JEZ1Dc6lhWE/J1ORUQ+oXIIsPGD+pAUH6N5BxEJKyqHAIuNjmLG8DQWbT+oU1pFJGyoHIJg9ogMDhyrZeu+Sq+jiIj4ReUQBBd8fEqrzloSkTChcgiCjKQEzs5O5j0t4S0iYULlECTTh6axrugI1XU6pVVEQp/KIUim5aVR3+hY/dERr6OIiJyRyiFIJub2ITbaWK7rHUQkDKgcgiQxLoZxA/uwbJfKQURCn8ohiKbnpbG59BgV1XVeRxEROS2VQxBNH5qKc/DBLt0dTkRCm8ohiMbm9KZnXLSGlkQk5AW9HMxsrpltN7MCM7v7FNvMMrMPzWyzmb0X7IyBEhsdxaTBfVleoCMHEQltQS0HM4sGHgDmAaOAa8xsVKttegMPApc550YD84OZMdCmD02jsLyKfUdPeB1FROSUgn3kMAkocM4VOufqgGeAy1ttcy3wgnOuCMA5F1FrTkzLSwPQ0YOIhLRgl0M2UNzicYnvuZaGA33MbLGZrTGz60/1ZmZ2i5nlm1l+WVl4LE0xsl8SfXvGad5BREJasMvBTvJc63WsY4DxwKXAHOCHZjb8ZG/mnHvEOTfBOTchPT29a5MGSFSUMTUvleUFh7SEt4iErGCXQwmQ0+LxAKD0JNu84Zyrcs6VA+8DY4OULyim56Wx/1gNheVVXkcRETmpYJfDamCYmQ02szjgauCVVtu8DMwwsxgzSwQmA1uDnDOgpg9NBdBSGiISsoJaDs65BuAO4E2af+AvdM5tNrMFZrbAt81W4A1gA7AKeMw5tymYOQNtYN9Esnv3YJkmpUUkRMUE+wOdc68Br7V67qFWj38J/DKYuYLJzJg+NJU3Nx+gsckRHXWyqRgREe/oCmmPTB+axtET9WwpPeZ1FBGRNjpUDmY2wTdnIB00Na953kGntIpIKGp3OZhZf2A5cGXXx+k+MpISGJbRi2WalBaRENSRI4cvA38Cbu7iLN3O9KFprP7oMLUNjV5HERH5Fx0phy8B3wPizCyvi/N0K9PyUqmpb2JdUYXXUURE/kW7ysHMZgPbfBenPQHcFJBU3cTkIalEma53EJHQ094jh5uAP/i+fxaYb2Y646mDUnrEMmZAb5bp5j8iEmL8/sHuW0p7CvA6gHPuGLACuCQgybqJ6XmprC+u4Hhtg9dRREQ+4Xc5OOcqnHND3b+uFneDc+7VAOTqNqYPTaOhybF692Gvo4iIfKKzQ0IP+9Y/wsxmdkGebmf8oD7Ex0SxVPMOIhJCOrt8xn8AfzCzBuBDmldQlXZIiI1mQm4fXe8gIiGls0cOPwW203xPhoWdj9M9TctLY9v+SsqP13odRUQE6Hw5fNc592PgVpqPIqQDpg/13TpUZy2JSIjwqxzMbODJnvdd74Bzrgr4Whfm6lbGZKeQlBCj6x1EJGT4O+fwhpll0HwPho0032thI7DRd0orzjmtAdFB0VHG1CGpWoRPREKGX0cOzrlRQBbwdZqvbRgK/BDYZma7Axev+5g+NI3iwycoOlTtdRQRkXZd51DnnFsHvAisBPYDJ4D1AcrWrXw876CjBxEJBf7OOYwws2+Z2bs0L9c9FXgKOMs599kA5us28tJ7kpkcr1NaRSQk+DvnsBVYB/wceMU5p3Muu5iZMT0vjcU7ymhqckTp1qEi4iF/h5VuBT4A7gCKzWyrmS00sx+a2WcDlq6bmTY0jcNVdWzbX+l1FBHp5vw6cnDOPdzysZkNAM4BxgCfB17q8mTd0PShzbcOXb6rnFFZyR6nEZHurEMXwTnnSpxzrznn7nHOfamrQ3VX/VN6MCS9p9ZZEhHPnbEczGymmQ0xs7/4hpK0wF4ATc9LY9Xuw9Q1NHkdRUS6MX+OHK4BfgB8C/gizfMPEiDTh6ZRXdfI+pIKr6OISDfmTzmMBjKdcwedc3XA0QBn6tam+m4dunSnhpZExDv+lMMPgXtaPH4zQFkESEmM5ezsFJbrYjgR8dAZy8E5955zruV9GoYHMI/QvIT3uqIKqnTrUBHxiD8T0gtbfP0NuDkIubq18323Dl2lW4eKiEf8GVY65py70vc1H/hnoEN1dxNy+xAXE6WlNETEM/6Uw3+3evz/OvOBZjbXzLabWYGZ3X2a7SaaWaOZfaEznxeOEmKjGT+wD8t08x8R8Yg/cw67Acws0czGOuc+Gesws4Fmlu3vh5lZNPAAMA8YBVxjZqNOsd09dOPJ7/OHpbF13zHdOlREPNGeK6TrgRfMrGeL5x4D+rfjPSYBBc65Qt9psc8Al59kuzuB54GD7XjviDItr3kpjQ909CAiHmjP/Rzqab6Xw1Xwya1D051z+e34vGyguMXjEt9zn/AdiVwBPHSmNzOzW8ws38zyy8rK2hEj9H1861DNO4iIF9q7ttJjwI2+768Hnmjnv3+ydahdq8e/Ae7y57ajzrlHnHMTnHMT0tPT2xkltMVERzEtL5UlO8txrvUfkYhIYLWrHJxz2wDMbDjNy2r8uZ2fVwLktHg8AChttc0E4Bkz+wj4AvBgd10WfMawdPZWnKCwvMrrKCLSzXRkVdY/0HwEscE5d6Sd/+5qYJiZDTazOOBq4JWWGzjnBjvncp1zucBzwG3OuZc6kDPsXTC8+Wjo/R2RNWQmIqGvI+WwEBhLc0m0i3OugeYbBr1J893lFjrnNpvZAjNb0IEsES2nbyK5qYks0TpLIhJk/t4m9BPOuWogpaMf6Jx7DXit1XMnnXx2zt3Q0c+JFDOGpfPcmhJqGxqJj4n2Oo6IdBMdutmPBM/M4emcqG9kzZ72juCJiHScyiHETRnSl5go4/0dGloSkeBROYS4pIRYzhvUhyU7NSktIsGjcggDM4elsbn0GGWVWkpDRIJD5RAGZvpOadXV0iISLCqHMDA6K4U+ibG8r6ElEQkSlUMYiI4yzh+WrqU0RCRoVA5hYsawNMoqa9m2v9LrKCLSDagcwsTMYVpKQ0SCR+UQJvqlJDA8s5eW0hCRoFA5hJGZw9JZ9dFhTtSdcTVzEZFOUTmEkRnD06lraGLlbt0dTkQCS+UQRiYP7ktcTJSW0hCRgFM5hJGE2GgmD+6rpTREJOBUDmFmxrA0dh48TmnFCa+jiEgEUzmEmY+X0liqs5ZEJIBUDmFmRGYSGUnxvKehJREJIJVDmDEzZgxLZ+nOchqbtJSGiASGyiEMXTAinaMn6llXpLvDiUhgqBzC0KwR6cRGG29tOeB1FBGJUCqHMJScEMvUvDTe3Lxfq7SKSECoHMLUxaMy2XOomp0Hj3sdRUQikMohTF00KhOAtzbv9ziJiEQilUOYykxO4Nyc3pp3EJGAUDmEsYtHZ7Kh5KiulhaRLqdyCGMXj+oHwNs6ehCRLqZyCGNDM3qRl96Tt7Zo3kFEupbKIcxdPLofKwoPc7S63usoIhJBVA5h7uJRmTQ2Od7drqElEek6QS8HM5trZtvNrMDM7j7J69eZ2Qbf13IzGxvsjOFk7IDeZCTF89ZmlYOIdJ2gloOZRQMPAPOAUcA1Zjaq1Wa7gQucc+cAPwUeCWbGcBMVZVw0KpP3dpRRU697S4tI1wj2kcMkoMA5V+icqwOeAS5vuYFzbrlz7uMV5VYAA4KcMezMGd2P6rpGlhXoHg8i0jWCXQ7ZQHGLxyW+507lJuD1U71oZreYWb6Z5ZeVdd/7G0wZkkpSfIyGlkSkywS7HOwkz5105Tgzm01zOdx1qjdzzj3inJvgnJuQnp7eRRHDT1xMFLNHZvDPrQd0jwcR6RLBLocSIKfF4wFAaeuNzOwc4DHgcufcoSBlC2sXj87kUFUda/boHg8i0nnBLofVwDAzG2xmccDVwCstNzCzgcALwJecczuCnC9sXTA8nbjoKC3EJyJdIqjl4JxrAO4A3gS2Agudc5vNbIGZLfBt9iMgFXjQzD40s/xgZgxXSQmxTBuayltbDugeDyLSaTHB/kDn3GvAa62ee6jF9zcDNwc7VySYM7of33thI9sPVDKyX7LXcUQkjOkK6QjyqbMyMIM3NmloSUQ6R+UQQTKSEpgyOJUX1u6lSWctiUgnqBwizJUTB1B0uJqVuw97HSVsNDY5dhyoZOeBSooPV1NWWUtlTT0NjU1eRxPxTNDnHCSw5p3dnx+9vJmF+cVMzUv1Ok7IKq04wZKdZby/o5ylBeUcPXHyVW1jooyhGb24ckIOnzsvm96JcUFOKuINlUOESYiN5rKxWTy3poSfXD6a5IRYryOFjE17j/LC2r28v7OMgoPHAchMjufiUZlMzUslNjqKmvpG31cTNfWNVNc3srygnP98dQs/f2Mbc0f34+pJOUwZnEpU1Mmu6RSJDCqHCHTVxByeWlnEKx+W8sUpg7yO47nahkbue2cnD71XSEyUMXlIKldPzGHm8HSGZfTC7Mw/5LeUHuPZ1UW8uG4vr6wvZVBqIldNzOGGabkkxmk3kshjkXJO/IQJE1x+vi6JAHDOMe+3S4iPieLlO873Oo6nNu09yrcXrmf7gUrmjx/ADz49ipQeHT+aqqlv5PVN+3hmVTErdx8mp28P/ueKMcwY1n2Xb5HwZmZrnHMTWj+vCekIZGZcOSGH9SVH2bb/mNdxPFHX0MS9b+/g8geWUXGijidumMgv54/tVDFA87DdFeMG8OzXpvLsLVOIjYriS39YxbcXrqeiuq6L0ot4T+UQoa4Yl01cdBTPri4+88YRZkvpMS5/YBn3vbOTy8/N4q1vXsDskRld/jmTh6Ty2jdmcPvsPF7+cC//du97vLqhVFeoS0RQOUSoPj3juGhUJi+t20ttQ/e5CdALa0u4/IGllB+v5dHrJ3DvleeSkhi4SfmE2Gi+M2ckr9xxPv1TenDHX9fx1Sfz2X+0JmCfKRIMKocIduXEHI5U1/PPLQe9jhJwzjnue2cn31q4nkmD+/LWN2dy0ajMoH3+qKxkXrxtGj+49CyWFpRzyX1LWFmoBYUlfKkcItj5Q9PISklgYX5kDy3VNzZx9/MbufftHXzuvGyeuGESfXoG/3qEmOgobp4xhH98fQa9E2O57rGV/HnFHg0zSVhSOUSw6CjjC+MH8P7OMkorTngdJyCO1zZw05/yeTa/mK9/ahi/mj+WuBhv/1rnpffipdunM3N4Oj98aRPff3EjdQ262lrCi8ohws2fkINz8NyaEq+jdLkDx2q48qEPWFZQzj2fH8O3Lhru1zULwZCcEMuj10/gtll5PL2qmGsfXUFZZa3XsUT8pnKIcDl9E5mWl8rf1hRH1GJ8Ow5UcsUDy9hzqIrHb5jIVRMHeh2pjego47tzR/K7a8axqfQol92/lA0lFV7HEvGLyqEbuGpiDsWHT7AiQiZIiw9Xc91jK2locjz7talcMDy0L0D7zNgsnr91GlFmzH/oA97YtM/rSCJnpOv+u4E5o/uRnBDDwvxipg1N8zpOpxyuquPLj6+itr6R52+dxrDMJK8j+WV0Vgqv3DGdm5/M59an1vKTy0Zz/dRcr2OFJeccheVVfLDrEGv2HOFIdR3VtY0cr22guq6BqrpGqmobiI4ysnv3IKt3D7J6J5DVuwfZvXswoE8iY7JTPJ+bCnUqh24gITaay8/NZmF+Mf+vspb0pHivI3XIibpGbv7TakoqTvCXmyaHTTF8LLVXPH+9eQp3Pr2OH728mf1Ha/jOnBEhM08Sqpxz7DlUzQeFh/hg1yFWFB7ioG/+JiMpnn4pCSTGRZPVO4HEuBh6xsfQMy6a+sYm9lbUUFpxgrVFR6io/r+Vd5PiY7hgRDoXjcpk9sgMLVB5EiqHbuLG6bn8dVURDy4u4D8+M9rrOO3W0NjEnU+vY11xBb+/7jwmDe7rdaQO6REXzUNfPI8fvryZBxfvYv+xGu75/DnERuu32NYqa+p56cNSnlqxh237KwFIT4pnypBUpg5JZWpeKrmpiX6X6/HaBvZVnGBXWRWLth3knW0HeHXDPmKijKl5qVw0KpM5o/uRmZwQyP+ssKGF97qR7z63npfWlbL4O7PI6t3D6zh+c87x/Rc38fSqIn5y2Wi+PC3X60id5pzjd+8WcO/bO5gxLI3ff3E8veL1uxrA5tKjPLWyiJfX7aWqrpHRWcnMHz+A84elk5fes8uOtBqbHB8WH+GtzQd4e8sBCsuriDK4cGQG10wayKwRGUR3g2XZT7XwnsqhGyk5Us3s/13MF8bn8LPPjfE6jt9+985OfvX2Dm6dlcddc0d6HadLLVxdzPde3MhZ/ZN4/IaJZCR1z99aGxqbeGV9KX9esYd1RRXEx0Rx2dgsrpsyiLEDUoIy9FZwsJLn1+7lb/kllB+vpX9KAldOyOHKiTlkh9EvU+2lchAAfvTyJv66soh3vn0Bg1J7eh3njBauLua7z2/gc+Oy+dWVYyNyfH7RtoPc9tRa+vaM45HrxzM6K8XrSEHT1OR4deM+fv32DnaXVzEkvSdfnDyIz583IKBrYp1OfWMT72w9wF9XFbNkZxkAs4anc8vMPKYM6RtxfwdVDgLAwWM1zPzlIuad3Z9fX3Wu13FO663N+7n1qbVMy0vl8RsmRvS4/Ka9R/nqk/lUVNdz75VjmTemv9eRAso5xztbD/K/b21n2/5KRmQm8e2Lh3PRqMyQ+uFbfLiahfnFPL2qmPLjtUzM7cOdFw5jxrC0kMrZGSoH+cTPXtvKI0sKeeubM0P2jJ8lO8u46Y/5jMpK5i83T+4W4/EHK2v42p/XsK6ogm/+2zC+fuGwiLwV6fKCcn7x5nY+LK4gNzWRf79oOJ85Jyuk/1tr6ht5ZlURD71XyP5jNYzN6c3XLxzKhSMzwr4kVA7yicNVdcz8xaJPJkJDzeqPDvOlP6wkN7Unz94y1bPhBS/U1Dfy/Rc38sLavVwyph//O39sxNyGtOBgJf/56lbe31FG/5QEvvGpYXx+/ICwOiKsbWjkuTUl/H7xLkqOnGB0VjK3zx7KnNH9wnbyWuUg/+Let3dw3zs7efXO8zk7O3TGuDeWHOXaR1eQnhTPs1+bGrbXZHSGc47HluzmZ69vZWS/ZB798oSwnhA9Wl3Pr/+5gz+v2ENiXDTf+NQwvjhlEAmx0V5H67D6xiZeWreXBxfvYnd5FTl9e3DT9MHMn5BDzzA7ylU5yL84VlPPjHsWcd7A3jxx4ySv4wDN6yVd9fAHJMbF8LcFU8PqdNtAWLT9IF//6zriY6O498pzmRniy4S01tDYxNOri7n3re0cPVHPNZMG8q2LhpPaK3IKv7HJ8faW/Ty6ZDdr9hwhOSGG66YM4oZpuWFzvYTKQdp4cHEBv3hjO8/fOpXxg7y9qOyj8irmP/wBBvxtwdSwOJMqGAoOHue2p9aw48BxvjpjMN+ZMzIsln1YXlDOf766hW37K5kypC8/+vRoRmUlex0roNbsOcJjSwp5c/N+oqOMy8Zmc/3UQZwTpFNxO0rlIG1U1zUw8xeLGZrRk6e/OsWzv8ClFSeY/9AHVNc1sPBrU0N2ktwrNfWN/Nc/tvCXFUWcnZ3MfVePY0h6L69jndSGkgp++eZ2luwsZ0CfHvzg0rOYM7pfSP9w7GpFh6p5fNluFuYXU13XyNnZyVw3eRCXjc0KySGnkCkHM5sL/BaIBh5zzv281evme/0SoBq4wTm39kzvq3LomCeW7eYnf9/Cn74yyZPVTfcdPcHVj6zg8PE6nr5lSkjNf4SaNzfv567nN1DX0MSPLxvN/PEDQuaH7s4DlfzqrR28sXk/fRJjuX320LCfV+isYzX1vLxuL0+tLGLb/kp6xcdwxbhsrp08kLP6h85RVEiUg5lFAzuAi4ASYDVwjXNuS4ttLgHupLkcJgO/dc5NPtN7qxw6pqa+kbm/eZ/KmgZevG06A1MTg/bZLYvhTzdN4ryBfYL22eFq39ET/PuzH7Ki8DCfPqc///3ZMZ6ezVV8uJpf/3MHL67bS8+4GL46YwhfOT+XJC1k9wnnHGuLjvDUiiJe3biPuoYmRvZLYlpeGtPyUpk0pK+nC/+FSjlMBX7snJvje/w9AOfcz1ps8zCw2Dn3tO/xdmCWc+60i+CrHDpuV9lxPvfgclJ7xfHCrdPonRj4+y+rGDquscnx0Hu7uPftHSTGRnP9tEHcOH0waUGc6N267xiPL93NSx/uJcqMG6blsuCCPE/u3R1OjlTV8fzaEhZtP0j+R0eobWgiymBMdgpT89KYPLgvA1MT6Z+S4NcpzDX1jRyuqiMzOaHDp9KGSjl8AZjrnLvZ9/hLwGTn3B0ttnkV+Llzbqnv8TvAXc65Nj/5zewW4BaAgQMHjt+zZ08Q/isi06rdh/niYys5N6c3T940KaDDAaUVJ7jmURVDZ20pPcb9i3by+qb9xMdEcfXEgXx15pCAnfba1ORYtP0gf1i6m+W7DtEjNpr5EwZw26yh9EsJjzNzQklNfSPriip8S5GXs66ogoYWd2tMToihf0oP+vdOoH9KAs5B+fE6DlXVcriqjkPH6zhe2wDA0rtmM6BPx476Q6Uc5gNzWpXDJOfcnS22+Qfws1bl8F3n3JrTvbeOHDrvlfWlfP3pdXxmbBa/vercgFyx2rIYnrxpEuNUDJ1WcPA4D7+3ixfX7QXgs+OyWXDBEIZmdM3EflVtA8+vLeGJZR+xu7yK/ikJXD81l2sm5QTlKLO7qKptYNPeo+w7WuP7OsG+ozXs930fZUZqr3hSe8aR2iuO1J7xvn/GMW9Mf1J6dGxo6lTlEOyp8xIgp8XjAUBpB7aRALhsbBYlR6r5xRvbyenTg+928QqoKobAGJrRi1/OH8s3LxrOo+8X8szqIp5bU0JuaiIzhqVz/rA0pual+j2ufbCyhrV7KlhXdIS1RUfYUHKU2oYmxub05r5rxjHv7H5hdVVzuOgZH8PkIalex/hEsI8cYmiekP4UsJfmCelrnXObW2xzKXAH/zchfZ9z7oxXaenIoWu0vHfC/1wxhmsnD+yS912ys4zvvbCRo9X1KoYAO3S8lpc/LGVpQTkrCg9RXddIdJQxdkAK5w9LJz0pntr6Ruoam6itb6K2oYm6hibKjteyrugIJUdOABAXHcXo7GTOG9iHS8b0Z/wg/T+LRCExrOQLcgnwG5pPZX3cOfffZrYAwDn3kO9U1vuBuTSfynrjyeYbWlM5dJ2GxiZufjKfJTvLeezLE5g9IqPD71V0qJr/+scW3tpygEGpifz26nGcm9O768LKadU1NLGu6AhLC8pZsrOcDSUVNLXa5eNiooiPjiK5Ryxjc1I4b2Afxg3sw9nZycTHdN9TUbuLkCmHQFE5dK3jtQ1c9fAH7DhQydUTB3L77PZNOlbXNfDgol08sqSQmCjjjguHctP5g/XDxmOVNfXU1DcRHxtFXHTzVyivhiqBp3KQdjt0vJZ7397BwvxizIxrJw3ktll5ZJxmzRjnHH/fsI+fvbaVfUdr+Oy5Wdw97yydzSISolQO0mHFh6u5/90CnltbQkyU8aUpg1gwK4+UHrEUHDzOltJjbNl37JN/Hj1Rz+isZH5y2Wgm5Hq7ZpOInJ7KQTptz6Eq7nungBfXlRAbHYVzUNfYBEB8TBQj+yczqn8ykwf35TNjs8J2fXuR7kTlIF2msOw4f1r+ET3iYhiVlcyo/knkpvYkRqc3ioSdULnOQSLAkPRe/OTys72OISIBpF/1RESkDZWDiIi0oXIQEZE2VA4iItKGykFERNpQOYiISBsqBxERaUPlICIibUTMFdJmVga05z6haUB5gOJ0VqhmC9VcoGwdpWwdE0nZBjnn0ls/GTHl0F5mln+yS8ZDQahmC9VcoGwdpWwd0x2yaVhJRETaUDmIiEgb3bkcHvE6wGmEarZQzQXK1lHK1jERn63bzjmIiMipdecjBxEROQWVg4iItBHR5WBmc81su5kVmNndJ3ndzOw+3+sbzOy8EMp2nS/TBjNbbmZjQyVbi+0mmlmjmX0hlLKZ2Swz+9DMNpvZe6GSzcxSzOzvZrbel+3GIOV63MwOmtmmU7zu5X5wpmxe7genzdZiOy/2gzNm6/R+4JyLyC8gGtgFDAHigPXAqFbbXAK8DhgwBVgZQtmmAX18388LpWwttnsXeA34QqhkA3oDW4CBvscZIZTt+8A9vu/TgcNAXBCyzQTOAzad4nVP9gM/s3myH/iTrcX/96DuB37+uXV6P4jkI4dJQIFzrtA5Vwc8A1zeapvLgSddsxVAbzPrHwrZnHPLnXNHfA9XAAOCkMuvbD53As8DB4OUy99s1wIvOOeKAJxzwcrnTzYHJJmZAb1oLoeGQAdzzr3v+6xT8Wo/OGM2D/cDf/7cwJv9wJ9snd4PIrkcsoHiFo9LfM+1d5tAaO/n3kTzb3bBcMZsZpYNXAE8FKRMH/Pnz2040MfMFpvZGjO7PoSy3Q+cBZQCG4FvOOeaghPvtLzaD9ormPvBGXm4H/ij0/tBTABChQo7yXOtz9v1Z5tA8PtzzWw2zTvF+QFN1OIjT/Jc62y/Ae5yzjU2/xIcNP5kiwHGA58CegAfmNkK59yOEMg2B/gQuBDIA942syXOuWMBznYmXu0HfvNgP/DHb/BmP/BHp/eDSC6HEiCnxeMBNP/G1t5tAsGvzzWzc4DHgHnOuUNByOVvtgnAM74dIg24xMwanHMvhUC2EqDcOVcFVJnZ+8BYINDl4E+2G4Gfu+ZB4AIz2w2MBFYFONuZeLUf+MWj/cAfXu0H/uj8fhCsCZRgf9FcfIXAYP5vgnB0q20u5V8n4laFULaBQAEwLdT+3Fpt/0eCNyHtz5/bWcA7vm0TgU3A2SGS7ffAj33fZwJ7gbQg/dnlcurJS0/2Az+zebIf+JOt1XZB2w/8/HPr9H4QsUcOzrkGM7sDeJPmMwoed85tNrMFvtcfovkMg0to/stXTfNvdqGS7UdAKvCg7zeTBheEVSD9zOYJf7I557aa2RvABqAJeMw5d9pTEYOVDfgp8Ecz20jzD+K7nHMBX/bZzJ4GZgFpZlYC/AcQ2yKXJ/uBn9k82Q/8zOaZM2Xriv1Ay2eIiEgbkXy2koiIdJDKQURE2lA5iIhIGyoHERFpQ+UgIiJtqBxE/GRmvc3stjNs87CZTQ9WJpFAUTmI+K83cNpyACbTvECcSFhTOYj47+dAnm+N/F+2ftHMzgJ2OOcaWz3/RzP7vZktMrNCM7vAtx7/VjP7o2+baN92m8xso5n9e1D+i0ROIWKvkBYJgLtpXoLg3FO8Pg944xSv9aF5wb3LgL8D04GbgdVmdi7NV1VnO+fOhuYhrC5LLdIBOnIQ6TpzOHU5/N01L0ewETjgnNvompfr3kzzGjmFwBAz+52ZzQW8XqlVujmVg0gXMLNEoLdz7lSrmdb6/tnU4vuPH8e45hvajAUWA7fTvAqpiGc0rCTiv0og6RSvzQYWdfSNzSwNqHPOPW9mu2he5VPEMyoHET855w6Z2TLfTd1fd859p8XL84DnOvH22cATZvbx0fz3OvFeIp2mVVlFuoCZrQUmO+fqvc4i0hVUDiIi0oYmpEVEpA2Vg4iItKFyEBGRNlQOIiLShspBRETaUDmIiEgb/x+4HHtIUobfuAAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"_=rho.plot()"
]
},
{
"cell_type": "markdown",
"id": "1f119a79",
"metadata": {},
"source": [
"### Reducible sequence\n",
"The above sequence is computationally expensive, requiring a basis set of 48 elements. If we consider that we're mainly interested in the $S^+$ operator of the $^{15}$N spin, then its worth noting that the $^{15}$N $\\pi$-pulse in the middle of the sequence converts $S^+$ into $S^-$, $S^\\alpha$, and $S^\\beta$, so that if we could get rid of it, we could use a basis set 1/4 as big (12 elements). Since there is no $^{15}$N chemical shift included in our simulation, we can switch the channel of the middle pulse to $^1$H. An alternative approach would be to use a $\\delta$-pulse on $^{15}$N, although the basis set would then only be reduced to 24 elements, since $S^+$ would still be converted to $S^-$"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "21706ef8",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"State-space reduction: 48->12\n"
]
}
],
"source": [
"centerH=L.Sequence().add_channel('1H',t=[0,L.taur/2-tp/2,L.taur/2+tp/2,L.taur],v1=[0,v1,0])\n",
"\n",
"rho=sl.Rho('15Np','15Nx')\n",
"\n",
"rho,f,s,c,Ueye=rho.ReducedSetup(first,second,centerH,L.Ueye())"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "65deeecf",
"metadata": {},
"outputs": [],
"source": [
"Ufirst=f.U()\n",
"Usecond=s.U()\n",
"Ucenter=c.U()"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "60b97a9e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"13.82811188697815\n"
]
}
],
"source": [
"U1=Ueye\n",
"U2=Ueye\n",
"\n",
"t0=time()\n",
"for k in range(48):\n",
" rho.reset()\n",
" (U2*Ucenter*U1*rho)()\n",
" U1=Ufirst*U1\n",
" U2=Usecond*U2\n",
"print(time()-t0)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "ee96b5e6",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAApeUlEQVR4nO3deXxU9f3v8dcnG5AQIDtIwpawbwqRLYqgdQNb26ptcavbz2K16739advb7dfbW7vYWq8LLm1d6q+KWneUugKVRYLsIEvCkrAlISSBQPbv/WOiN01AJsvMmZm8n49HHmZmDnPeBs68c77nnO8x5xwiIiItRXkdQEREQo/KQURE2lA5iIhIGyoHERFpQ+UgIiJtxHgdoKukpqa6IUOGeB1DRCSsrFmzpsw5l9b6+YgphyFDhpCfn+91DBGRsGJme072vIaVRESkDZWDiIi0oXIQEZE2gl4OZvYXMysxs02neN3M7D4z22lmG8xsUrAzioh0d17sOTwOXPIZr18KDG/+uhV4KAiZRESkhaCXg3NuKVD+GYtcDjzpfFYC/cxsQHDSiYgIhOYxh4FAUYvHxc3PtWFmt5pZvpnll5aWBiWciEh3EIrlYCd57qTzijvnHnHO5TrnctPS2lzD4ZfFmw/y/JriDv1ZEZFIFYoXwRUDWS0eZwL7A7Ei5xxPr9rLysLDjOqfyLiBfQOxGhGRsBOKew6vANc3n7U0Dah0zh0IxIrMjHu/eiYpCXHM/9saKo/XB2I1IiJhx4tTWf8OrABGmlmxmd1sZvPNbH7zIouAQmAn8CjwzUDmSU6I44FrJnGwsobvL1xHU5PujCciEvRhJefcvNO87oDbgxQHgEmDkvjx3NH84tUtPLSkgNtn5wRz9SIiIScUh5U8ccOMIcydMIB7/rmN5QVlXscREfGUyqGZmfGbKyYwNDWBb/99LQcra7yOJCLiGZVDC717xPDQtZOprm3kjv/+iPrGJq8jiYh4QuXQyoiMRO6+Yjz5e47w2zc/9jqOiIgnVA4ncfmZA7lu2mAeXbaLVYWHvY4jIhJ0KodT+PHc0aQkxPHg+wVeRxERCTqVwyn0jI3mxrwhLNleypb9VV7HEREJKpXDZ7hu2hAS4qJ5eKn2HkSke1E5fIa+8bHMmzKI1zYcoKj8uNdxRESCRuVwGjefO5Qog8eWFXodRUQkaFQOpzGgby++eOZAns0v4vCxWq/jiIgEhcrBD984bxg19U08sXy311FERIJC5eCHnPRELhqTwRMr9lBd2+B1HBGRgFM5+Gn+rGwqT9TzzOqi0y8sIhLmVA5+mjQoiSlDk/nzskLqGjTnkohENpVDO9w2K5v9lTW8sj4gdy0VEQkZKod2mDUijVH9E3l4SYHuGCciEU3l0A5mxvzzstlRcox3Py7xOo6ISMCoHNrpsgkDGNivF39dvsvrKCIiAaNyaKeY6CiumDSQFQWHKdNFcSISoVQOHTB3whk0OXhz00Gvo4iIBITKoQNGZPQmJ703r2844HUUEZGAUDl0gJkxd/wAVu06TMnRGq/jiIh0OZVDB82dMIAmB4s1tCQiEUjl0EEjMhIZnt6b1zS0JCIRSOXQCXMnDODD3eWUVGloSUQii8qhE+aOH4Bz8IaGlkQkwqgcOmF4RiIjMxJ5faOGlkQksqgcOmnuhAGs3l3OIQ0tiUgECXo5mNklZrbNzHaa2V0neb2vmb1qZuvNbLOZ3RjsjO0x55OhJe09iEgECWo5mFk08ABwKTAGmGdmY1otdjuwxTk3EZgF3GNmccHM2R456b0Z1V9DSyISWYK95zAF2OmcK3TO1QHPAJe3WsYBiWZmQG+gHAjpe3POHT+A1buPcLBSQ0siEhmCXQ4DgZb32Sxufq6l+4HRwH5gI/Ad51xI33ptzoQBACzS3oOIRIhgl4Od5LnWd825GFgHnAGcCdxvZn1O+mZmt5pZvpnll5aWdmXOdslO683oAX00tCQiESPY5VAMZLV4nIlvD6GlG4F/OJ+dwC5g1MnezDn3iHMu1zmXm5aWFpDA/rpswgDW7DnC/ooTnuYQEekKwS6H1cBwMxvafJD5a8ArrZbZC1wAYGYZwEigMKgpO2DOeN/Qki6IE5FIENRycM41AHcAi4GtwELn3GYzm29m85sX+yUww8w2Au8AdzrnyoKZsyOGpiYw9ow+vL6h9Y6QiEj4iQn2Cp1zi4BFrZ5b0OL7/cBFwc7VFeaMH8DvFm9jf8UJzujXy+s4IiIdpiuku9DFY/sD8O7HJR4nERHpHJVDF8pOSyAruRfvb1M5iEh4Uzl0ITNj9sh0Pth5mJr6Rq/jiIh0mMqhi80emc6J+kY+3FXudRQRkQ5TOXSxacNS6BETxXsaWhKRMKZy6GK94qKZnp3Ckm3eXbEtItJZKocAmD0yncKyanaXVXsdRUSkQ1QOATB7ZDqAzloSkbClcgiAQSnxDEtL4D0NLYlImFI5BMjskemsKDzMiTqd0ioi4UflECCzR6ZT19DEisKQnxZKRKQNlUOAnD00ifi4aN77WENLIhJ+VA4B0iMmmrycVN7bVoJzre9nJCIS2lQOATRrZBrFR05QUHrM6ygiIu2icgigWc2ntGpoSUTCjcohgAb268XIjERNpSEiYUflEGCzRqWxenc5R2vqvY4iIuI3lUOAzR6ZTn2j44Odh72OIiLiN5VDgE0enERijxhNpSEiYUXlEGCx0VGcOyKV97eV6pRWEQkbKocgmDUynYNVNXx88KjXUURE/KJyCIJZI9IAdNaSiIQNlUMQpPfpydgz+ugGQCISNlQOQZKXk8ravRWapVVEwoLKIUhmZKdQ19jE6t3lXkcRETktlUOQTBmaTGy08UGBpvAWkdCncgiS+LgYzspKYrkuhhORMKByCKIZOSls2l9JxfE6r6OIiHwmlUMQ5eWk4hysLNTeg4iENpVDEE3M7Ed8XLTmWRKRkBf0cjCzS8xsm5ntNLO7TrHMLDNbZ2abzWxJsDMGSlxMFFOGJuugtIiEvKCWg5lFAw8AlwJjgHlmNqbVMv2AB4EvOOfGAlcFM2Og5WWnUlhazcHKGq+jiIicUrD3HKYAO51zhc65OuAZ4PJWy1wN/MM5txfAORdRc07MyEkB4IOd2nsQkdAV7HIYCBS1eFzc/FxLI4AkM3vfzNaY2fWnejMzu9XM8s0sv7Q0PKamGN2/D8kJcRpaEpGQFuxysJM813oe6xhgMjAXuBj4iZmNONmbOececc7lOudy09LSujZpgERFGdOHpbB852FN4S0iISvY5VAMZLV4nAnsP8kybzrnqp1zZcBSYGKQ8gXFjJwUDlbVUFhW7XUUEZGTCnY5rAaGm9lQM4sDvga80mqZl4FzzSzGzOKBqcDWIOcMqLzsVACW67iDiISooJaDc64BuANYjO8Df6FzbrOZzTez+c3LbAXeBDYAHwKPOec2BTNnoA1OiWdgv1663kFEQlZMsFfonFsELGr13IJWj38H/C6YuYLJzJiRncI/txyisckRHXWyQzEiIt7RFdIeyctJpfJEPVv2V3kdRUSkjQ6Vg5nlNh8zkA769HoHndIqIiGo3eVgZgOA5cBXuj5O95Ge2JMRGb11MZyIhKSO7Dl8HXgCuKWLs3Q7M7JTWb27nNoG3TpUREJLR8rhOuCHQJyZZXdxnm4lLyeVmvom1u6t8DqKiMi/aVc5mNls4OPmi9P+CtwckFTdxNRhyUSZrncQkdDT3j2Hm4E/N3//LHCVmemMpw7q0zOWCZn9+KBA1zuISGjx+4O9eSrtacAbAM65KmAlMCcgybqJvJwU1hVVcKy2wesoIiKf8rscnHMVzrkc9++zxd3gnHstALm6jbzsVBqbHB/u0t6DiISOzg4JPdw8/xFmNrML8nQ7kwYn0SMmSlNpiEhI6ez0GT8D/mxmDcA6fDOoSjv0jI0md0iSrncQkZDS2T2HXwLb8N2TYWHn43RPM7JT+fjgUcqO1XodRUQE6Hw5/Kdz7ufAbfj2IqQD8nKap/DWWUsiEiL8KgczG3Sy55uvd8A5Vw18owtzdSvjB/YlsWeMrncQkZDh7zGHN80sHd89GDbiu9fCRmBj8ymtOOc0B0QHRTffOlST8IlIqPBrz8E5NwY4A/g2vmsbcoCfAB+b2a7Axes+8nJSKSo/wd7Dx72OIiLSrusc6pxza4EXgVXAQeAEsD5A2bqVPE3hLSIhxN9jDiPN7Ptm9i6+6bqnA08Do51zXwxgvm4jO6036Yk9dEqriIQEf485bAXWAncDrzjndM5lFzMz8nJSWbq9lKYmR5RuHSoiHvJ3WOk2YAVwB1BkZlvNbKGZ/cTMvhiwdN1MXk4qh6vr2HboqNdRRKSb82vPwTn3cMvHZpYJTADGA1cAL3V5sm7o0+MOO8sYPaCPx2lEpDvr0EVwzrli59wi59xvnHPXdXWo7mpA314MS03QcQcR8dxpy8HMZprZMDP7W/NQkibYC6AZOSl8uKuc+sYmr6OISDfmz57DPOB/Ad8HrsV3/EECJC87leq6RtYXVXgdRUS6MX/KYSyQ4Zwrcc7VAZUBztStTc9OwQxN4S0invKnHH4C/KbF48UByiJAv/g4xp3RVxfDiYinTlsOzrklzrmW92kYEcA8gu+4w9q9Rzhep1uHiog3/DkgvbDF13PALUHI1a3lZadS3+j4cFe511FEpJvyZ1ipyjn3leavq4C3Ax2quzt7SDJx0VG6v4OIeMafcvhVq8c/7swKzewSM9tmZjvN7K7PWO5sM2s0sys7s75w1CsumrMG9dP1DiLiGX+OOewCMLN4M5vonPt0rMPMBpnZQH9XZmbRwAPApcAYYJ6ZjTnFcr+hGx/8zstJZcuBKsqr67yOIiLdUHuukK4H/mFmCS2eewwY0I73mALsdM4VNp8W+wxw+UmW+xbwAlDSjveOKHk5qTgHKzS0JCIeaM/9HOrx3cvhq/DprUPTnHP57VjfQKCoxePi5uc+1bwn8iVgwenezMxuNbN8M8svLS1tR4zQNzGzL717xPAvDS2JiAfaO7fSY8CNzd9fD/y1nX/+ZPNQu1aP7wXu9Oe2o865R5xzuc653LS0tHZGCW0x0VFMG5bCsh2lONf6RyQiElj+3s8BAOfcx2aGmY3AN63GOe1cXzGQ1eJxJrC/1TK5wDNmBpAKzDGzBufcS+1cV9ibOSKVt7ceYs/h4wxJTTj9HxAR6SIdmZX1z/j2IDY4546088+uBoab2VAziwO+BrzScgHn3FDn3BDn3BDgeeCb3bEYAGYO9+0NLd0RWUNmIhL6OlIOC4GJ+EqiXZxzDfhuGLQY393lFjrnNpvZfDOb34EsEW1wSjxZyb1Yul3HHUQkuNo1rATgnDsO9O3oCp1zi4BFrZ476cFn59wNHV1PJDAzzh2exstr91Hf2ERsdIduvyEi0m76tAlxM4enUV3XyEd72juCJyLScSqHEDc9O4XoKGPZDg0tiUjwqBxCXN9esZyZ1Y9lOigtIkGkcggDM4ensWFfpabSEJGgUTmEgXNH+KbS0ER8IhIsKocwMGFgX/r0jNHQkogEjcohDMRER3HO8FSWbi/TVBoiEhQqhzBx7vA0DlbVsLPkmNdRRKQbUDmEiXOHpwKwVKe0ikgQqBzCRGZSPMPSEli6XccdRCTwVA5hZObwNFbtOkxN/WlnMxcR6RSVQxg5d3gqNfVNrNFUGiISYCqHMDJtWAqx0aahJREJOJVDGEnoEcPkwUk6KC0iAadyCDPnDk9j64EqSo7WeB1FRCKYyiHMnDfCd3e4f2nvQUQCSOUQZsYM6ENyQpym8BaRgFI5hJmoKOOcnFSW7SijqUlTaYhIYKgcwtB5I9IoO1bLxn2VXkcRkQilcghD549KJ8rgrS2HvI4iIhFK5RCGkhLimDI0mcWbD3odRUQilMohTF08tj87So5RWKpZWkWk66kcwtSFYzIADS2JSGCoHMJUZlI8Y8/oo6ElEQkIlUMYu3hsf9YWVVBSpaulRaRrqRzC2EVjM3AO3t5a4nUUEYkwKocwNjIjkUHJ8RpaEpEup3IIY2bGxWMzWF5QxtGaeq/jiEgEUTmEuYvG9qe+0fH+Nt3jQUS6TtDLwcwuMbNtZrbTzO46yevXmNmG5q/lZjYx2BnDyaRBSaQkxGloSUS6VFDLwcyigQeAS4ExwDwzG9NqsV3Aec65CcAvgUeCmTHcREcZF47J4P1tpdQ26N7SItI1gr3nMAXY6ZwrdM7VAc8Al7dcwDm33Dn3yU2SVwKZQc4Ydi4am8Gx2gZWFBz2OoqIRIhgl8NAoKjF4+Lm507lZuCNU71oZreaWb6Z5ZeWdt8x9xnZqSTERbN4s66WFpGuEexysJM8d9KbEpjZbHzlcOep3sw594hzLtc5l5uWltZFEcNPz9hoZo1M560th3SPBxHpEsEuh2Igq8XjTGB/64XMbALwGHC5c05jJX64aGwGZcdqWVtU4XUUEYkAwS6H1cBwMxtqZnHA14BXWi5gZoOAfwDXOee2Bzlf2Jo1Mp2YKOOfOmtJRLpAUMvBOdcA3AEsBrYCC51zm81svpnNb17sp0AK8KCZrTOz/GBmDFd9e8UyPTuFxZsP4pyGlkSkc2KCvULn3CJgUavnFrT4/hbglmDnigQXje3PT17axM6SYwzPSPQ6joiEMV0hHUEuHO27x8ObmzS0JCKdo3KIIP379mTKkGReXLtPQ0si0ikqhwhzZW4mhWXVrNlz5PQLCwB1DU2sL6pg8/5K9h4+zpHqOuobm7yOJeKpoB9zkMCaO34AP39lMwvzi8gdkux1nJBVcbyO97aV8PbWEpZuK+VobUObZXrGRtG7Rywj+/fmikmZXDpuAL3ioj1IKxJ8KocIk9AjhssmDOC1DQf42efHktBDf8WfOFhZwyvr9/H21hLW7DlCY5MjtXcP5owfwMwRaURHGUdr6jlW28DRmgaO1TZQdaKe5QWH+f7C9fz05c3MHT+Aq3IzmTw4CbOTXdMpEhn0yRGBvpKbxcL8Yl7feICv5Gad/g9EuMYmx5MrdvO7xds4XtfIqP6JfHNWNheMzmDCwL5ERX32h3xTk+PD3eU8v6aYVzfs59n8IoalJnBlbiY3zBhCfJw2I4k8FikHLnNzc11+vi6JAHDOccE9S0jt3YOF86d7HcdT2w8d5c4XNrB2bwUzR6Txiy+MZWhqQoff71htA4s2HuD5/GI+3F1OVnIv7v7yBPJyUrswtUjwmNka51xu6+d1QDoCmRlX5Wbx4e5yCkuPeR3HE7UNjfzxre3MvW8Zu8uq+eNXJ/LEjWd3qhgAeveI8e2ZzZ/Os7dOI9qMax5bxQ//sYEq3Y1PIojKIUJdMWkg0VHG82uKvY4SdGv2HOGy+/7Fn97ZwaXjBvDW98/jS2dldvkxgqnDUnjzuzP5xsxhPLu6iIv+sJR3tmpmXIkMKocIld6nJ7NGpPHCR8U0dKPTMp9etYerFiynuraBv9yQy33zziK1d4+Ara9nbDQ/nDOaF7+ZR59eMdz8RD7ffWYt5dV1AVunSDCoHCLYVblZHKqqZdmOMq+jBJxzjj+9vYMfv7iJ80aksfh7Mzl/VEbQ1j8xqx+vfuscvn3BcF7bcIC59y1jY3Fl0NYv0tVUDhHs/FHppCTEsTC/6PQLh7HGJsfPXtnMH9/ezpcnDeSR63NJ7Bkb9Bw9YqL5/oUjeOn2PKLMuHLBcl5auy/oOUS6gsohgsXFRPHFswby9tZDETvMUdvQyLefWcuTK/Zw68xh/P7KicRGe/vPetzAvrx8Rx4Ts/rx3WfX8avXt3SroT2JDCqHCPeV3CzqG11E/gZ7rLaBmx5fzesbDvCjOaP40ZzRp71mIVhSe/fg6Vumcv30wTy6bBc3Pr6aiuORWdASmVQOEW5k/0QmZvZlYX5RRE3Gd/hYLfMeWcnKwnLuuWoit87M9jpSG7HRUfzX5eP4zRXjWVVYzhfu/4BtB496HUvELyqHbuCq3Cw+PniUTfuqvI7SJSqP13PNY6vYUXKUR6+fzBWTM72O9Jm+evYg/n7rNGrqG/nSgx/w3rYSryOJnJau++8GPj/xDH752hYW5hcxPrOv13E65URdIzc9sZrC0mr+csPZnDM8PK5Mnjw4iVe/dQ43Pb6aW57I5+4vj+cqTW3SYaVHa1m16zAf7angcHUtR2saOFpT3/xf35xYMdFGZlI8g5LjyUzuRVZSPFnJ8QxOjmdwSrzmxjoNlUM30LdXLJeO68/L6/Zx16WjwnYyvvrGJm57eg1r9x7hgasnhU0xfCKjT0+e/cZ0bvvbGn7w/AYOVdVw++wcfUj54WBlDat2HWZlYTmrdh2msLQagF6x0aQl9iCxZwyJPWPISo4nsWcMfXrGUt/YRNGRE2w5UMVbWw5R1+KkgAF9e3L+qHQuGJ3OjOxUesZqtt3WwvNTQtrthryhvLRuP3/9YBd3nD/c6zjt1tTk+J/Pref9baX8+svjuXT8AK8jdUjvHjH8+etnc+cLG/j9P7dzsKqGX3xhHNEhciA9lNTUN/LahgM8tXIP64sqAEjsEcPZQ5P5am4WU4elMO6MPsT4cXZaU5Pj0NEaispPsLPkGEu2l/Di2n08vWovPWOjmJGdyvmj0rloTAbpfXoG+P8sPGjivW7kP57MZ2XhYZb952z6xcd5Hcdvzjl+8eoWHl++mx9cPJLbZ+d4HanTmpocv128jQVLCrhoTAb3zTtLv70223O4mqdX7WVhfhEVx+vJSe/NVZMzyctJZfSAPl1WpLUNjawqLOfdj0t45+NDFJWfIDrKuHB0BtdNH8yM7JRusVd3qon3VA7dyMcHq7j0T8uYf142d14yyus4frvvnR384a3t3HLOUH48d3REbbCPf7CLX7y2hcmDknjs67lhVdpdyTnHe9tKeGL5HpZsLyUmyrh4bH+unTaYacOSA/537pxjZ8kxnl9TzML8Io4cr2dYWgLXTB3MlZMy6Rsf/Isqg0XlIAB895m1vLn5IEt/MDssdp+fWrGbn7y8mS9PGsjvr5wYMtcxdKXXNxzge8+uIzOpF499PZdhab29jhQ0zjmW7ijjnn9uY0NxJRl9enD1lMF8bUoWGR79+6ypb2TRxgP8beUePtpbQc/YKC6fOJBbzh3K8IxETzIFkspBANhdVs3n/rCEq6cO4r8uH+d1nM/06vr9fPuZtVwwKp2Hrp3s+ZXPgbR6dznzn1pDfWMT9189iZkj0ryOFHAf7irn94u38eHucjKTevGdC4bzxbMGhtTf8+b9lfxt5V5eWruPE/WNXDQmg2/OzuHMrH5eR+syKgf51I9e3Mhz+UW8+z9mkZUc73Wck1q2o5SbHl/NmVn9eOrmqd1iPL6o/Dj/8WQ+2w8d5SeXjeGGGUMiagjtExuKK/j9P7ezdHsp6Yk9+Nb5OXz17EHExYROKbRWXl3H48t38/gHu6iqaSAvJ4VvzsqJiOMSKgf51MHKGs773XtcNuEM7vnKRK/jtLG+qIJ5j65kUHI8z35jOn17Re54b2vVtQ1899l1vLXlEPOmZPGLL4wL6Q/N9thXcYJfL9rKaxsOkBQfy22zsrlu2hB6xYVP8R+rbeC/V+3h0WW7KD1ay8Ssftx2XjYXjckI2yFPlYP8m/+zaCuPLStk8XdnhtQ4akHpMa5asIL4uGj+cduMsDgu0tWamhz3vLWNB94rYMrQZBZcO5nkhPA9UH28roEFSwp5eEkBZnDrzGz+49yhnsyc21Vq6ht54aNiFiwpoKj8BNlpCXxjZjaXn3UGPWLCp+xA5SCtlFfXMfO373FOTioLrpvsdRzAt0dzxUPLqW1o5Ln5Mzp9S89w9/K6ffzg+Q2kJ/bg/qsnhd04t3OOV9bv5+43PuZAZQ2fn3gGd106ioH9enkdrcs0NDaxaNNBFrxfwJYDVWT06cHN5wxl3pRBYVN+Kgdp4963t3Pv2zt4+Xbf9NJeqjhex1ceXsH+ihqeuXUa4waG9zQfXWVdUQW3P/0Rh6pq+N6FI5h/XnZYXDC3obiC/3p1C/l7jjBuYB9+9vmxnD0k2etYAeOcY9mOMhYsKWB5wWESe8Zw7bTBXDN1EJlJoXlc7xMqB2njaE09M3/7HuMG9uWpm6d6luNEXSPX/nkVG4srefyms5mRHV7TYgRa5Yl6fvTiRl7fcIBpw5L541fPZEDf0Pztu7D0GPe8tZ3XNxwgtXcc/3nxKK6cnBm24/Edsb6ogoeXFvDGpoMAnDs8jXlnZ/G5MRkhdSbWJ0KmHMzsEuBPQDTwmHPu7lavW/Prc4DjwA3OuY9O974qh455dGkhv1q0lf++ZSozcoL/oXy8roFbn1zDBwVlPHj1pLCdFiPQnHM8v6aYn72ymdjoKH5zxXguGRc6P6sDlSf409s7eG5NMT1iorj5nKHcOnNY2AytBELxkeMszC/mufwiDlTWkNo7jismZ/K1sweF1JBpSJSDmUUD24ELgWJgNTDPObelxTJzgG/hK4epwJ+cc6f9tVbl0DE19Y187g9LqG9s4qXb84L6G+nRmnpufjyf/D3l/PbKiVwZ4lNvh4JdZdV855m1bCiuZN6ULH5y2Rji47ybIq28uo6H3t/JEyv24JzjmqmDuX12DmmJPTzLFGoamxxLtpfw9w+LePfjEhqbHBOz+nFuTirnDE9l0qAkT89IC5VymA783Dl3cfPjHwI4537dYpmHgfedc39vfrwNmOWcO/BZ761y6LitB6q48qHlDE5J4Ln504Mya2vl8Xqu/+uHbN5Xyb1fO5PLJpwR8HVGirqGJv749nYWLCkgJSGOW84dxrXTBtM7iLPt7q84wVMr9/DUij0cr2vgS2dl8t3PDQ/Z62ZCxaGqGp5fU8w7Ww+xvriSxiZHfFw0U4cmc87wNKYPSyEruRe9e8T4ff1EXUMTMVHW4aG7UCmHK4FLnHO3ND++DpjqnLujxTKvAXc75/7V/Pgd4E7nXJtPfjO7FbgVYNCgQZP37NkThP+LyPTexyXc/MRqzh+VwcPXTQ7oQc/Dx2q57s8fsrPkGPdffRYXje0fsHVFsjV7yrn37R0s21FG316x3JQ3lBtmDAnYPEDOOT7aW8FfPtjFm5sO4pzj4rH9+d6FIxgRQqdDh4uqmnpWFBzmXzvK+NfOMnaVVX/6WnxcNBl9epKe2OPT/wKUH6/jSHUd5cfrOVLt+/5obQPL7zqfMzp4FtipyiHY+6Mn+8Rp3U7+LON70rlHgEfAt+fQuWjd2+xR6fz0sjH8/NUt3P3GVn48d0xA1lNSVcPVj62iqPw4j349l/O6wTQRgTJ5cDJP3TyVdUUV3P/uTv749nYeXVbIddMHc/M5Q0nt3TVDO3UNTSzaeIC/frCL9cWVJPaM4aa8IVw/fYj2FDqhT89YLh7bn4ubfzkqPnKcj/ZWcKiyhoNVNRyqqqGkqpb1xRUcqqohyoyk+DiSE+JISohjaEo8SQlxJMfHBWQGgWCXQzHQ8vZXmcD+DiwjAXBD3lB2lVXz6LJdDE3tzdVTB3Xp+++rOME1j66k5GgtT9w0hWnDUrr0/burM7P68djXc9l6oIoH3tvJgiUFPLq0kLMG9WNGtm9ce2JmP7/HtWvqG9m4r5L83UdYs6ec/D1HqGiepfSXl4/ly5Myw/aGUaEsMyk+pE57DfawUgy+A9IXAPvwHZC+2jm3ucUyc4E7+P8HpO9zzk053XvrmEPXaGhs4pYn81m2o4wnbpzSZXdbW7q9lLte2MDR2gYev3EKkwcndcn7SlsFpcd4Lr+Y5QVlbNxXiXO+YYopQ5PJy04lLbEH9Y1N1Dc6GpqaqGtooqHJUXa0ljV7j7BpXyX1jb7PhWGpCUwenMScCQM4b3hatzoltbsIiWMOzUHmAPfiO5X1L865X5nZfADn3ILmU1nvBy7BdyrrjSc73tCayqHrHK2p58qHVrC/8gQvfnMGOekdH08uPVrL/359Cy+v28+w1ATum3eWLnALosrj9awoPMzygjI+2FlGQWn1KZeNi45ifGZfcgcnMbn5K6WLhqYkdIVMOQSKyqFrFR85zhcfWE50FHznghFcMXlgu+aMaWpyPJtfxK8XbaWmvonbZmVz26zsbjG7aigrOVrDsZoGYqOjiI2OIibaiI2KIjbGiIuO8uuWmxJZVA7Sbpv2VfLjFzeyvvkmLLecM4yrpw467XjzjkNH+dGLG1m9+whThybzqy+NJye9+9zARiScqBykQ5xzfLDzMA++v5PlBYfpFx/L16cP8Z0y2SuWg1U1FJZWU1h2jMLSagpKj7Gy8DAJPWL48ZzRXDk5M+znuxeJZCoH6bS1e4/w4PsFvLXlED1jozCME/WNn76eEBfNsLTeTBrUj29fMFzj1SJhIFSuc5AwdtagJB69PpdtB4/y9Ko9xEZHMSwtgWGpvRmWlkB6Yg/tJYhECJWDtNvI/okhf/9pEekcnZogIiJtqBxERKQNlYOIiLShchARkTZUDiIi0obKQURE2lA5iIhIGyoHERFpI2KmzzCzUqA99wlNBcoCFKezQjVbqOYCZesoZeuYSMo22DnX5paMEVMO7WVm+SebTyQUhGq2UM0FytZRytYx3SGbhpVERKQNlYOIiLTRncvhEa8DfIZQzRaquUDZOkrZOibis3XbYw4iInJq3XnPQURETkHlICIibUR0OZjZJWa2zcx2mtldJ3ndzOy+5tc3mNmkEMp2TXOmDWa23Mwmhkq2FsudbWaNZnZlKGUzs1lmts7MNpvZklDJZmZ9zexVM1vfnO3GIOX6i5mVmNmmU7zu5XZwumxebgefma3Fcl5sB6fN1untwDkXkV9ANFAADAPigPXAmFbLzAHeAAyYBqwKoWwzgKTm7y8NpWwtlnsXWARcGSrZgH7AFmBQ8+P0EMr2I+A3zd+nAeVAXBCyzQQmAZtO8bon24Gf2TzZDvzJ1uLvPajbgZ8/t05vB5G85zAF2OmcK3TO1QHPAJe3WuZy4EnnsxLoZ2YDQiGbc265c+5I88OVQGYQcvmVrdm3gBeAkiDl8jfb1cA/nHN7AZxzwcrnTzYHJJrvRtu98ZVDQ6CDOeeWNq/rVLzaDk6bzcPtwJ+fG3izHfiTrdPbQSSXw0CgqMXj4ubn2rtMILR3vTfj+80uGE6bzcwGAl8CFgQp0yf8+bmNAJLM7H0zW2Nm14dQtvuB0cB+YCPwHedcU3DifSavtoP2CuZ2cFoebgf+6PR2EBOAUKHCTvJc6/N2/VkmEPxer5nNxrdRnBPQRC1WeZLnWme7F7jTOdfo+yU4aPzJFgNMBi4AegErzGylc257CGS7GFgHnA9kA2+Z2TLnXFWAs52OV9uB3zzYDvxxL95sB/7o9HYQyeVQDGS1eJyJ7ze29i4TCH6t18wmAI8BlzrnDgchl7/ZcoFnmjeIVGCOmTU4514KgWzFQJlzrhqoNrOlwEQg0OXgT7YbgbudbxB4p5ntAkYBHwY42+l4tR34xaPtwB9ebQf+6Px2EKwDKMH+wld8hcBQ/v8BwrGtlpnLvx+I+zCEsg0CdgIzQu3n1mr5xwneAWl/fm6jgXeal40HNgHjQiTbQ8DPm7/PAPYBqUH62Q3h1AcvPdkO/MzmyXbgT7ZWywVtO/Dz59bp7SBi9xyccw1mdgewGN8ZBX9xzm02s/nNry/Ad4bBHHz/+I7j+80uVLL9FEgBHmz+zaTBBWEWSD+zecKfbM65rWb2JrABaAIec8595qmIwcoG/BJ43Mw24vsgvtM5F/Bpn83s78AsINXMioGfAbEtcnmyHfiZzZPtwM9snjldtq7YDjR9hoiItBHJZyuJiEgHqRxERKQNlYOIiLShchARkTZUDiIi0obKQcRPZtbPzL55mmUeNrO8YGUSCRSVg4j/+gGfWQ7AVHwTxImENZWDiP/uBrKb58j/XesXzWw0sN0519jq+cfN7CEze8/MCs3svOb5+Lea2ePNy0Q3L7fJzDaa2feC8n8kcgoRe4W0SADchW8KgjNP8fqlwJuneC0J34R7XwBeBfKAW4DVZnYmvquqBzrnxoFvCKvLUot0gPYcRLrOxZy6HF51vukINgKHnHMbnW+67s345sgpBIaZ2f81s0sAr2dqlW5O5SDSBcwsHujnnDvVbKa1zf9tavH9J49jnO+GNhOB94Hb8c1CKuIZDSuJ+O8okHiK12YD73X0jc0sFahzzr1gZgX4ZvkU8YzKQcRPzrnDZvZB803d33DO/aDFy5cCz3fi7QcCfzWzT/bmf9iJ9xLpNM3KKtIFzOwjYKpzrt7rLCJdQeUgIiJt6IC0iIi0oXIQEZE2VA4iItKGykFERNpQOYiISBsqBxERaeP/AX35aOvghBI0AAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"_=rho.plot()"
]
},
{
"cell_type": "markdown",
"id": "d7812428",
"metadata": {},
"source": [
"We see that the reduced sequence yields the same results as before, but with matrices only 1/16 as big as the previous calculation."
]
},
{
"cell_type": "markdown",
"id": "e793e4dd",
"metadata": {},
"source": [
"## Sweep the correlation time\n",
"We now re-run the above setup, but while varying the correlation time to observe the transition between having a fully averaged dipole coupling and having the rigid-limit coupling."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "10f1ed81",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"log10(tc /s) = -6.0, 12 seconds elapsed\n",
"log10(tc /s) = -5.6, 23 seconds elapsed\n",
"log10(tc /s) = -5.1, 34 seconds elapsed\n",
"log10(tc /s) = -4.7, 45 seconds elapsed\n",
"log10(tc /s) = -4.3, 56 seconds elapsed\n",
"log10(tc /s) = -3.9, 67 seconds elapsed\n",
"log10(tc /s) = -3.4, 78 seconds elapsed\n",
"log10(tc /s) = -3.0, 89 seconds elapsed\n"
]
}
],
"source": [
"rho_list=[]\n",
"legend=[]\n",
"t0=time()\n",
"for tc in np.logspace(-6,-3,8):\n",
" L.kex=sl.Tools.nSite_sym(n=3,tc=tc)\n",
"\n",
" rho_list.append(rho.copy_reduced())\n",
"\n",
" Ufirst=f.U()\n",
" Usecond=s.U()\n",
" Ucenter=c.U()\n",
"\n",
" U1=Ueye\n",
" U2=Ueye\n",
"\n",
" for k in range(28):\n",
" rho_list[-1].reset()\n",
" (U2*Ucenter*U1*rho_list[-1])()\n",
" \n",
" U1=Ufirst*U1\n",
" U2=Usecond*U2\n",
"\n",
" legend.append(fr'$\\log_{{10}}(\\tau_c)$ = {np.log10(tc):.1f}')\n",
" print(f'log10(tc /s) = {np.log10(tc):.1f}, {time()-t0:.0f} seconds elapsed')"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "43b12e5c",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAE8CAYAAABQLQCwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABvJElEQVR4nO3dd3hc5ZX48e+ZUe+S1at7lXvDFBtTbSCAaaEEQkJCSIDfZpNsQnaTbHbTN8lmUwglIZAQAiH0ENOrMbg33LtsybaK1XuZ8/tjxo4Qki3NjDQzmvN5nnks3bl37pF8NHPue98iqooxxhhjjAkfjkAHYIwxxhhjhpYVgMYYY4wxYcYKQGOMMcaYMGMFoDHGGGNMmLEC0BhjjDEmzFgBaIwxxhgTZqwANMYYY4wJM1YAGmOMMcaEGSsAvSQiB0XkggCc90ci8mUfjl8jIlP8GJLxA8sn40+WT8bfLKeGHysAQ4iIZAC3AA/48DI/A/7bPxGdnohcLyI7RKRJRPaJyDl97JcmIs969isRkRuHKsZwFWr5JCJvi0iriDR6HrtOs3+/cs/4x3DOJxG5S0TWiUibiDwyFPGZ0MupE0RknCe3/tzH8409Hl0i8uuhjBGsAAw1twLLVbXFh9d4AVgsIjn+CalvInIh8BPgM0AisBDY38fu9wLtQBZwE3CfXbUNulsJoXzyuEtVEzyPCX3tNMDcM/5xK8M0n4AjwPeBPwxRXMbtVkIvp8D9eba2rye75VwC7s+8FuBvQxXcCVYA+oGITPJcTdaKyDYRubzbc7NEZKOINIjI30TkryLyfS9PtRR4p8e5d3iuINo9jxNXFJN6ewFVbQXWAxd5GcNA/Bfw36q6SlVdqlqmqmU9dxKReOBq4Nuq2qiq7+H+o715CGIMOpZPftGv3AsHlk++U9VnVPU54HigYwkGllN9E5HrgVrgjX4ecg1QAawYrJj6YgWgj0QkEvg78CqQCdwNPCYiE0QkCngWeARIAx4HlvlwuqnAR25TqOokz1XEo8B/dbuy2HGK19kBTO/j53nR80fd2+PF/gYqIk5gDpAhIntFpFREfiMisb3sPh7oUtXd3bZtBsKuBdDy6bR+JCJVIrJSRM7t45wDyb1hzfLptE6bT+ajLKf6JiJJuG83f3UAh30a+JOq6kDO5Q9WAPruDCAB+LGqtqvqm8CLwA2e5yKAX6lqh6o+A6w5caCIJIu7g2qjiBR32/4DEVkhIk+JSFy3c6UADX3EMQ3Y2s+YGzyv9TGqepmqpvTxuKyfrw/uZu1I3Fc35wAzgJnAt3rZNwGo67GtDvetu3Bj+dS3bwCjgTzgQeDvIjKml/0GknvDneVT3/qbT+ajLKf69j3gIVU93J+dRaQQWAT8cYDn8QsrAH2XCxxWVVe3bSW431RygbIelX33xGgGLgWeOrHB80cxRlXPAV4HPttt/xp6KYpExAFMpv9/DIm4m6j9QkRu6tYM/5Jn84k+G79W1aOqWgX8L3BJLy/RCCT12JZE33/4w5nlU+/5hKquVtUGVW1T1T8CK+k9nwaSe8Od5ZPv+WQ+ynKql5wSkRnABcAvBvBStwDvqeoBf8U2EFYA+u4IUOBJyBMKgTLgKJAnItLtuYITX3iukCp7vN45wIk3qZeAs7s9twX37dKeCnH/X36kk7uIfFJE3vXc3riu21OTcN9i/RgReUk+PkLpY2+e3anqY92a4Zd6ttUApUB/mrV3AxEiMq7btunAtn4cO9xYPvWST31QQD62cWC5N9xZPvmYT+ZjLKd6z6lzgZHAIRE5BnwNuFpENvT2Gh63EKDWP7AC0B9WA03A10UkUtz9SD4BPAF8AHQBd4lIhIhcAcw7zeul8s/boXW4+1GcsBx3c3FPSZ4Yok5s8FyN3Aqcr6pn4e6XgYhEA7OB13o7uaou7ZbYPR+nevPszcPA3SKSKSKpwJdx3yroec4m4Bngv0UkXkTOAq7A3ccj3Fg+9UJEUkTkYhGJ8fzsN+Ee2ftKH4f0K/fCgOVTLwaaT559YgAn4DxxXH/PN8xYTvXuQWAM7i4nM4D7gX8AF/e2s4icibvVdMhH/55gBaCPVLUduBz3aKUq4LfALaq60/PcVcBtuJufP4X7Q6jtFC9ZAyR7vk4Gqrs99yfgEvl4Z/YduK9uakRkomfbZcD/qmqHJ84Oz/bLgbdV9cgAf1RvfA/3UPjdnhg3Aj+Ak1dd/95t3y8BsbhHQz0OfFFVw64F0PKpT5G4p+GoxP17uRu4UlV3Qa/51GfuhRPLpz4NNJ++hbtrwT24f08thGefUsupPqhqs6oeO/HA3bWp9USLZy859WngGVUNXFcnVbXHED5wXz19pse2R4Biz9dTgb94vr4duLvHvj8EvtyP8/wcuNjzdUSP8xcH+vdgD8snewTfw/LJHv5+WE4F70M8vyAzSERkEe5h7FW4Jzi+Hxitqkc9zy/H3VxcAjygqo+IyI9w94OowH1l1eTFeSfhvg3WBqxR1X/zw49jAszyyfiT5ZPxN8up0GEF4CATkdtx345KAPYB31TVfwQ2KhOqLJ+MP1k+GX+znAodVgAaY4wxxoQZGwRijDHGGBNmhnwYu4j8AfdonQpVLe7leQF+iXtCzmbgVlXd4Hluiec5J/B7Vf3x6c6Xnp6uI0eO9N8PYIaF9evXV6lqxkCPs3wyvbF8Mv5k+WT8qa98CsQ8Ro8Av8E9vLs3S4Fxnsd84D5gvrjX97wXuBD3JK9rReQFVd1+qpONHDmSdevW+Sl0M1yISIk3x1k+md5YPhl/snwy/tRXPg35LWBVfZePzvPT0xV4FkZW1VVAiojk4J5Mcq+q7lf3XENPePY1xpiQISK3i8g6EVlXWdlzUQRjBsbyyXgrGPsA5vHRtQNLPdv62v4x9gdh/MnyyfiTqj6oqnNUdU5GxoDv8hnzEZZPxlvBWAD2thZjX2s09jqE+XR/EKU1zfzx/YM+BWnCx+nyaV9lI0+sORSAyMxwtKW0lmc2lAY6DDNMvL+vihe3DMXCTybUBONahqV0WzwayMe9+HRUH9sH7M+rDnH/O/vITIxm6dQcrwM1BuCh9w7w17WHKc5Lpjgv+fQHGHMKv19xgNe2l7NwfAbpCdGBDseEuIdWHODDsjqWTMkmwhmMbT4mUIIxG14AbhG3M4A6zwzia4FxIjJKRKKA6z37DthXLhzPjIIU/u2pLRyoGvCE48Z8xDeWTGREfBRff2oLHV2uQIdjQty/XDCOts4u7n97X6BDMcPAJ+cWUNHQxlu7rPuK+aghLwBF5HHgA2CCiJSKyG0icoeI3OHZZTmwH9gL/A74EoCqdgJ3Aa/gXgj6SVXd5k0MUREO7r1pFpFO4Yt/Xk9Le5ePP5UJZ8mxkfz3FcVsP1rPg+/uD3Q4JsSNyUhg2cx8Hl1VQnl9a6DDMSFu8cRMMhKj+eta66ZiPioQo4BvUNUcVY1U1XxVfUhV71fV+z3Pq6reqapjVHWqqq7rduxyVR3vee4HvsSRlxLLLz45g13lDXz7+a0nFo02xitLirO5ZGo2v3xjD/sqGwMdjglx/3L+OLpcyr1v7Q10KCbERTodXDM7nzd3VnCszi4ozD8F4y3gIXPuhEzuPm8cT60v5cl1h09/gDGn8N3LpxAb6eSep7fgctkFhfFe4Yg4rp2TzxNrDlNW2xLocEyIu25OAS6Fp21wkekmrAtAcF9pnz02nW8/v41tR+oCHY4JYZmJMXzr0kmsPVjDY6u9msfVmJPuOm8cAL95c0+AIzGhblR6PGeMTuOvaw/bxak5KewLQKdD+OX1M0iKieS7L2yzW8HGJ9fMzueccen8+KWd1nJjfJKXEsv18wr427pSDh1vDnQ4JsTdMK+QQ9XNrNp/PNChmCAR9gUgwIiEaL58wTjWHqzhzZ0VgQ7HhDAR4YfLptLpUv7vtd2BDseEuDsXj3VfpL5hrYDGNxdPySY5NpIn1lp3J+NmBaDHJ+cWMCo9np+8vJMuayI3PihIi+OGeYU8u7GMw9XWcmO8l5UUw81nFPHsxlIbXGR8EhPpZNnMPF7eeoyapvZAh2OCgBWAHpFOB/928QR2lzfaLPzGZ19YNBoRuP8dm8vN+OaOc8cQ4XTwyMqDgQ7FhLhPzi2gvcvFsxvLAh2KCQJWAHaztDib6QUp/O9ru2ntsLkBjfdykmO5Zra7/5ZNvWB8kZ4QzWXTcnhmQykNrR2BDseEsEk5SUzPT+avaw9bf3djBWB3IsI3lkzgaF0rf/rgYKDDMSHuS+eOoUvVWgGNz25ZMJKm9i5ruTE+u35eIbvKG9hcarNehDsrAHs4c0w6i8ZncO9b+6hrtqtt472CtDiWzczj8TWHqGxoC3Q4JoRNz09mal4yf/qgxFpujE8unZZDpFN46cOjgQ7FBJgVgL34xpKJ1Ld2cJ+13BgffencMXR0ufj9ClsiznhPRLh5QRF7KxpZtb860OGYEJYUE8mZY9J5edsxu5gIc1YA9mJybhJXzsjj4ZUHbLSU8cnojAQum5bLo6tKqLZcMj64fHouKXGRPLrqYKBDMSFuSXE2Jceb2XmsIdChmACyArAPdywaQ1uny5aIMz6767yxNLd38fDKA4EOxYSwmEgn180p4JVt5TawyPjkgklZiMDLW48FOhQTQFYA9mFCdiJnjE7j0VUlNi+g8cn4rESWTMnmkfcP0tJuo8uN926aX4hLlb+sORToUEwIy0iMZm5RGq9sswIwnFkBeAqfXjCS0poW3rLVQYyPbj1rJA2tnfzDOl4bHxSNiGfR+AweX3OIji5XoMMxIezi4mx2HmvgYFVToEMxAWIF4ClcODmLnOQY/mhTwhgfzR+Vxuj0eB63lhvjo1sWFFHZ0GatN8YnF0/JArA8CmNDXgCKyBIR2SUie0Xknl6e/zcR2eR5bBWRLhFJ8zx3UEQ+9Dy3brBjjXA6uGl+ISv2VNkyTMYnIsIN8wpZX1LDLut4bXywaHwmBWmxPPpBSaBDMSEsPzWO4rwkXrYCMGwNaQEoIk7gXmApMBm4QUQmd99HVX+qqjNUdQbwTeAdVe0+78Fiz/NzhiLm6+cVEuV02Jut8dnVs/OJcjqsFdD4xOkQrp9byOoD1bbWtPHJkinZbDxUa4OKwtRQtwDOA/aq6n5VbQeeAK44xf43AI8PSWR9SE+I5tJpOTy9vpTGts5AhmJCXFp8FBcXZ/PMhlJbajCMicjtIrJORNZVVlZ69RpXzswDsJVBjE/5tKQ4G4DXtlsrYDga6gIwD+g+r0qpZ9vHiEgcsAR4uttmBV4VkfUicntfJ/HHG2x3tywooqGt095sw5Q/8+mGeQXUt3ay3AaDhC1VfVBV56jqnIyMDK9eIy8lljNGp/HMhlKbzDfM+ZJPYzMTGZMRb7eBw9RQF4DSy7a+3r0+Aazscfv3LFWdhfsW8p0isrC3A/3xBtvdjIIUpuUn86f3D9qbbRjyZz4tGD2CkSPi7Daw8dlVs/I5eLyZjYdrAx2KCWEXT8lm1f5qW/QgDA11AVgKFHT7Ph840se+19Pj9q+qHvH8WwE8i/uW8qATEW5ZMJI9FY18sO/4UJzSDFMnBoOsPVjDnnIbDGK8t7Q4m+gIB89sKA10KCaELSnOpsulvL6jPNChmCE21AXgWmCciIwSkSjcRd4LPXcSkWRgEfB8t23xIpJ44mvgImDrkEQNXDYth5S4SJuA1fjs6tn5RDqFx9fYKjPGe4kxkVw8JZsXtxylrdP6lBrvTM1LJjc5xqaDCUNDWgCqaidwF/AKsAN4UlW3icgdInJHt12XAa+qavcZKrOA90RkM7AG+IeqvjxUscdEOrlyRh6vbi+nttmayo330hOiuWhKNk/bYBDjo2Wz8qht7uCtnb73dTbhSUS4aEo2K/ZU2UpFYWbI5wFU1eWqOl5Vx6jqDzzb7lfV+7vt84iqXt/juP2qOt3zmHLi2KF07Zx82jtdvLC5r7vWxvTPjfMKqWvpsLU4jU/OGZtOekI0z26028DGe+dPyqSt08XKvVWBDsUMIVsJZACm5CYzJTeJJ9fZrTvjmwWjR5CfGstT6+2D23gvwungihm5vLmzwjrxG6/NHzWC+Cgnb+y0foDhxArAAbp2dj5by+rZfqQ+0KGYEOZwCFfNymflviqO1LYEOhwTwq6alUdHl/LiFrszYbwTFeFg4fgM3thRgctlM12ECysAB+iKGXlEOR38bb21AhrfXD0rD1WbzNf4ZnJOEhOyEnnG8sj44PxJWVQ0tLH1SF2gQzFDxArAAUqNj+LCKVk8t7GM9k5XoMMxIaxoRDxzR6bytE3ma3wgIlw1K4+Nh2rZb2uWGy8tnpCBCLyxoyLQoZghYgWgF66dnU9Ncwdv2LxJxkdXz8pnf2UTm2wyX+ODK2bkIQLPbbLbwMY7IxKimVWYav0Aw4hXBaCIzPHM4xeWzhmXQXZSjA0GMT67ZFoO0REOnrbJfI0PspNjWDB6BM9vKrPWZOO18ydlsrWsnmN1rYEOxQyBAReAIpIDvA9c5/9wQoPTIVwzO593dlfaH4rxSZJnMt+/b7bJfI1vrpyZR4ktDWd8cMGkLABrBQwT3rQAfhr4I/A5P8cSUq6ZnY9LsZYb47OrZ+dT19JhfW+MT5YUZxMV4eB5GwxivDQuM4GCtFh7LwoT3hSANwPfBKJEZIyf4wkZI9PjmTcqjafWWwd+45uzx6aTlRTN0zYnoPFBUkwkF0zK5MUtR+nosgFqZuBEhPMnZrFyr60KEg4GVACKyGJgp6pWAQ8Dtw1KVCHimtn5HKhqslsuxidOh3DlzDze3l1JZUNboMMxIezKGXkcb2rnvT22ooPxzgWTsmjrdPGerQoy7A20BfA24CHP138FrhWRsB1JvLQ4m+gIB8/ZLRfjo2tm5dPlUp7fZLlkvHfuhEySYyN5zvLIeGneqDQSoyNslosw0O/iTURSgDOAlwBUtR5YBVwyKJGFgMSYSC6cnMXfNx+xWy7GJ+OyEpmWn8zTG+yD23gvKsLBJVNzeHVbOU1tnYEOx4Sgk6uC7LRVQYa7fheAqlqrqmP1ox3eblXVFwchrpCxbGYeNc0dvLOrMtChmBB31cw8dhytZ+cxW2bQeG/ZzDxaOrp4dfuxQIdiQtT5kzKpbGjjwzJbFWQ48/X27QMiEgcgIgv9EE/IWTg+g7T4KJ61Wy7GR5dNz8XpEJ7baJP5Gu/NKUolLyXW8sh4bfGETJwO4eVtdhExnPlaAP4n8JCIPArM9UM8ISfS6eAT03J4fXs59a0dgQ7HhLD0hGgWjkvn+U1lduvFeM3hEC6fkct7e6tsUJHxSmp8FAtGj+DlrcdslothzNcC8HvALkCBJ/tzgIgsEZFdIrJXRO7p5flzRaRORDZ5Ht/p77GBcuXMPNo6Xby81a6WjG+unJnH0bpWVh+oDnQoJoQtm5lHl0t5cYu1AhrvLCnO5kBVE7vKGwIdihkkvhaAX1fV7wJfxN0aeEoi4gTuBZYCk4EbRGRyL7uuUNUZnsd/D/DYITejIIVR6fE8ax34jY8umpxNfJTTRpYbn4zPSmRSTpLlkfHaxVOyEYGXPrSGjeGqXwWgiBT2tt0zHyCq2gR8oR8vNQ/Yq6r7VbUdeAK4op+x+nLsoBIRrpyRx6oDxzlS2xLocEwIi41ysqQ4h+UfHqW1wyZiNd5bNjOXzaV17K9sDHQoZhCJyO0isk5E1lVW+m8wYkZiNHNHpvHS1qN+e00TXPrbAviyiFSJyAoR+a2I3CEiZ4lI0okdVLU/n1Z5wOFu35d6tvW0QEQ2i8hLIjJlgMcO2h/EqVw5MxdVeGGz3XIZboY6n5bNzKOhrZM3d9pyTMPRUOXT5dPzEIHnNtl70nCmqg+q6hxVnZORkeHX115anM3u8kb22UXEsNSvAlBVJwO5wP/DPfffWODbwE4ROTCA80lvL9/j+w1AkapOB34NPDeAY0/EO2h/EH0pGhHP7KJUnt1QZp1mh5mhzqcFY0aQmRjNs3b7blgaqnzKTo7hrDHpPLfR3pOMd5YUZwNY//ZhaiDzALar6kbgWWA1cAxoATYP4HylQEG37/OBj1yeqmq9qjZ6vl4ORIpIen+ODbQrZ+axq7yBHUet06zxntMhXDEjl7d3VVDT1B7ocEwIu3JmHoeqm9lwqCbQoZgQlJMcy4yCFLsNPEz1tw/gBBH5ioi8CbwPLAAeAyap6pUDON9aYJyIjBKRKOB64IUe58oWEfF8Pc8T4/H+HBtol03NIdIpPLOhNNChmBB35cw8OrqUf3xob7zGe0uKs4mJdFhrsvHaJVOz2VpWz+Hq5kCHYvysvy2AO4CbgPuAOar6FVV9zTMYo99UtRO4C3jF85pPquo2T5/COzy7XQNsFZHNwK+A69Wt12MHcv7BlhofxXkTM3lu0xE6bWk444PJOUmMz0qwUZzGJwnREVw0OZsXtxylvdPek8zALS3OAew28HDU3wLwi8AHuAuwwyKyQ0SeFJFvi8iVAzmhqi5X1fGqOkZVf+DZdr+q3u/5+jeqOkVVp6vqGar6/qmODTbXzC6gqrGNd/fY0nDGeyLClTPzWFdSY1fexifLZuZR29zB27tsUJEZuIK0OKbkJrHcbgMPO/0dBPKAqt6lqotUNRO4EHgEaAeuHsT4Qs65EzIYER/FU+vtNrDxzRUz3IPcn7H5JY0Pzh6Xzoj4KJ6z5SqNl5YWZ7PxUC1H62yas+HEq4mgVbXU0xr3E1W92d9BhbJIp4MrZuTx+vYKaputA7/xXl5KLGeOGcFTGw7b0nDGa5FOB5+YnsvrOyqoa7HlKs3ALZ3qvg38it0GHlZOWwCKyEIRGS0if/bc9l04FIGFsqtn59He5eLvNieg8dF1cwo4XN1iS8MZnyybmUd7p4uX7Tae8cKYjATGZyWw3FYFGVb60wJ4A/At4CvAp3D3BzSnMCU3mUk5SXYb2Pjs4inZJEZH8Ld1h0+/szF9mJafzOj0eOtOYLz2iWm5rDlYTWmN9UkeLvpTAE4BslS1wjPqt26QYxoWrpmdz+bSOvbYQtrGB7FRTj4xI5flW4/S0Gq374x3RIRlM/NYfaCaMluu0njhypnuPsk2M8Hw0Z8C8NvAT7p9/8ogxTKsXDEjlwiH8JTNCWh8dO3sfFo7XPxji92+M9478QH+tN2ZMF4oSItj3sg0nrGVZYaN0xaAqvqOqr7bbdP4QYxn2EhPiObcCZk8u6HM5gQ0PplRkMLYzASetNvAxgcFaXGcNXYEf11rg4qMd66alcf+yia2lNqNwOGgP4NAnuz2+BvwuSGIa1i4ZnY+FQ1trNhbFehQTAgTEa6bk8+GQ7XsrbBF2Y33bphXSFlti70nGa8snZpDVIStLDNc9OcWcL2qXud5XAu8PthBDRfnTcwkNS7SbrkYn105Mw+nQ/jbemsFNN67cHIWafFRPLHmUKBDMSEoOTaSCydl8ffNR+iwO1shrz8FYM8VN/5jMAIZjqIi3HMCvrq9nJommxPQeC8zMYbFEzJ4xroUGB9ERzi5elYer20vp7KhLdDhmBC0bGYex5vaeXe3rXYV6vrTB/AAgIjEich0VT05IZmIFIpI3mAGGOo+ObeA9k4XT9tgEOOja+cUUNlgywwa31w/r5BOl9p7kvHKogkZpMVH8YzdBg55A1kJpAN4RkTiu237PZDj35CGl0k5ScwqTOEvaw7ZyCnjk/MmZjIiPoon19oHt/HemIwE5o1K4wl7TzJeiHQ6+MS0HF7bXm4ry4S4fheAqtoBPAt8Etytf0CGqq4bpNiGjRvnF7G/solV+201B+O9SKeDq2bl8fqOcirqWwMdjglhN8wr4ODxZj7YfzzQoZgQdNWsfNo7Xbz0oU1NFcoGuhbw74HPeL6+BXjYv+EMT5dNyyEpJoLHVpcEOhQT4m6cX0SnS3l8jQ0GMd5bWux+T3rC8sh4YVp+MqMz4u02cIgbUAGoqjsBRGQ87iXiHh2MoIabmEgnV8/O55Vtx6hqtI7Xxnuj0uNZND6Dx1aX2Cg847WYSCdXzcrn5a3HbICaGTAR4aqZeaw5UM3halsaLlQNtAUQ4CHcLYFbVLVmoAeLyBIR2SUie0Xknl6ev0lEtnge74vI9G7PHRSRD0Vkk4iE1K3nm+YX0tGltj6w8dktC4qoaGjj1W3lgQ7FhLAb5hXS3mUD1Ix3ls3KxyHw2GqbUihUeVMAPglMx10IDoiIOIF7gaXAZOAGEZncY7cDwCJVnQZ8D3iwx/OLVXWGqs4ZcOQBNDYzkXmj0vjL6kM2C7/xybkTMslPjeVPHxwMdCgmhE3ITmSmZ4CavSeZgcpLiWVJcTaPrzlEc3tnoMMxXhhwAaiqzaqarKreTAg9D9irqvtVtR14Ariix+u/361lcRWQ78V5gtJN8ws5VN3Myn02C7/xntMh3HxGEasPVLPzWH2gwzEh7NYzR7K/som3d1cEOhQTgj571ijqWjp4ZoP1BQxF3rQA+iIP6N7ruNSzrS+3AS91+16BV0VkvYjc3tdBInK7iKwTkXWVlcEzZ9qS4mzS4qN4bJU1mYeSYMyn6+YUEB3h4NEPbGBRqAmmfLpkag45yTH87t0DAY3DeC+Q+TS7KJVp+ck8vPKAtSKHoKEuAKWXbb1mjYgsxl0AfqPb5rNUdRbuW8h3isjC3o5V1QdVdY6qzsnIyPA1Zr+JjnByzex8XrNpPEJKMOZTanwUn5iey7Mby6hvtbm4Qkkw5VOk08FnzxrFB/uP82FpXUBjMd4JZD6JCLedPYp9lU28YxPUh5yhLgBLgYJu3+cDR3ruJCLTcA80uUJVT05UpapHPP9W4J6TcN6gRjsIbphXSJdL+YutxWl89OkFI2lu77K1po1Prp9XQGJ0BL9bsT/QoZgQtLQ4h6ykaP7wnrUih5qhLgDXAuNEZJSIRAHXAy9038EzwfQzwM2qurvb9ngRSTzxNXARsHXIIveTUenxXDApk0feP0hTm3WcNd6bmp/MjIIUHv2gxG6/GK8lxkRy/bwC/vHhUcpqWwIdjgkxUREOblkwkhV7qthd3hDocMwADGkBqKqdwF3AK8AO4ElV3SYid4jIHZ7dvgOMAH7bY7qXLOA9EdkMrAH+oaovD2X8/vLFc8dS29zB49YKaHz06TOL2F/VxHt7bWCR8d5nzhqFAA9bK47xwg3zComOcPDwyoOBDsUMwFC3AKKqy1V1vKqOUdUfeLbdr6r3e77+nKqmeqZ6OTndi2fk8HTPY8qJY0PR7KJUzhidxu9W7KetsyvQ4ZgQdsnUHDITo/nNW3ttXVfjtdyUWC6dlsMTaw9bn1IzYGnxUVw1K49nNpRSbROLh4whLwCN252Lx1Je38azNnze+CA6wsmXzh3DmgPVfLDP1nU13vv8OaNpbOvkCbszYbzw2bNG0dbpsjtbIcQKwAA5e2w6U/OSue+dfXTakl7GB9fPKyQ7KYZfvL7bWgGN14rzklkwegQPrzxoywyaARuXlcjC8Rk89N4BGqwVOSRYARggIsKdi8dQcryZ5VuPBTocE8JiIp18afEY1h6sYeVeawU03rt94WiO1rXanQnjlX+7aALVTe088I6NKA8FVgAG0EWTsxmTEc9vrf+W8dEn5xZYK6Dx2bkTMphZmMLPXt1ly3uZAZuan8wVM3L5/Xv7OVZnc90GOysAA8jhEL547lh2HmvgrV22FJPxXnSEkzsXj2F9SQ0r9tiIYOMdEeFbl06ioqHNWnGMV7520QRcLvj5q7sCHYo5DSsAA+yKGbnkpcTymzetFdD45rq5BeQmWyug8c3sojQunZbDg+9aK44ZuIK0OD59ZhFPbSi1tcqDnBWAARbpdPClxWPYcKiW5R8Or76A7Z0uapra6bJJiodEdISTLy0ey8ZDtbxrrYDGB/csmUiXS/mZteIYL9y5eCyJ0RH8aPnOQIdiTiEi0AEYuH5uIY+tOsT3XtzOuRMyiI8Orf8WVWX70Xpe317BvspGSmuaKattoaKhDVVwCKTGRZEWH8WIhCim5LpHG84bnUZSTGSgwx9WrptTwH1v7+N/X93FOWPTcTh6W37bmFMrSIvjM2eN5MEV+7n1zJEU5yUHOiQTQlLiorj7vHH8YPkO3ttTxdnj0gMdkulFaFUaw5TTIXzvyilcfd8H/PrNvdyzdGKgQzotVWVLaR3Ltx7l5a3HKDnejEMgPzWOvJRYzhmXQV5KLEmxkdQ2t3O8qZ3jjW1UNrTx6KoSHnrvAA6BqXnJLByfwTWz8ykaER/oHyvkRUU4+OpF4/nKk5v58+oSblkwMtAhea2ioZVjda1UN7WffLR3uchIiCYzKYbMxGgyE6NJi49CxApdf/vS4rE8ue4wP/jHDv7y+fn2OzYDcvOCIh55/yA/emkHfx9ztl2MBiErAIPE7KI0rpmdz+9X7Oea2fmMzUwIdEi9UlVe31HBL9/YzdayeiIcwplj07lj0RgumpzFiITo075Ga0cXGw/V8sG+Kj7Yf5x739rLr9/cy1ljR3D93EIumpJFdIRzCH6a4WnZzDye3VjGT17ayXkTM8lPjQt0SP1S29zOB/uOs3JfFe/vPc7+qqZ+HZedFMO8UWnMHZXGvJFpjMtMsA8bP0iOjeRfLxzPd57fxus7KrhwclagQzIhJCbSydeXTOBfntjEI+8f5LNnjwp0SKYHKwCDyD1LJ/LKtmN894VtPHrbvKC64lZVXt1ezq/e2MO2I/UUjYjjh8umcsnUbFLiogb0WjGRThaMGcGCMSMAOFbXyt/WHeaJtYe5+/GNpMVHceO8Qm49ayTp/SgozUeJCD9cNpWL/+9d/uPZrTzymblBlUvddbmUN3dW8KcPDvLe3ipUIT7KyfzRI7hxfiFFI+JJi486+Yh0CpUN7pbkioY2jtS2sOlwLasPHOeFzUcASE+I5rJpOVw5M4/p+clB+7OHghvmFfLI+wf53ovbOWN0Gokh3GWjtaOL+tYOmtu6aGzrpLm9C4e48yUjMTrkut6Egsun5/L3zUf40Us7mDMylWn5KYEOyXQjw3204Jw5c3TdunWBDqPf/vj+Qf7zhW389qZZXDI1J9DhALDxUA3ffn4rW8vchd/d543jyhm5RDj9O4bI5VJW7K3isVUlvLajnCing2tm53P7wtF+vz0sIutPrDM9EKGUT4+sPMB3/76d/71uOlfNyg90OB9R09TOX9cd5tEPSiirbSE7KYZr5+SzaHwG0wtSiBxgbqkqh6tbWHOwmte3l/Pmzgrau1yMHBHH5TPyuHZ2PgVpg9cSOpzzae3Baq5/cBVLi7P59Q0zQ6KgbmjtYMOhWrYdqWPbkXq2H6nn4PEmTvVxFxvpJDMpmnGZiRTnJVGcm0xxXjJZSdFD/jMPp3yqaWrn0l+tIMLp4MX/d7b1+w6AvvLJCsAg09nl4vLfrKSmuZ3Xv7IooFel9a0d/OyVXTy6qoSsxBi+dvGEQSn8erOvspHfvbufZzaU0elysXRqDnefN5aJ2Ul+ef3h9AbbF5dLufaBD9hX2chr/7qIjMTAt6Y2tnXy4Dv7+N2KA7R0dHHG6DQ+vWAkF07O8mte1bV08MrWYzy/uYz3PWsknz02nZvmF3L+pKwBF5inM9zz6d639vLTV3bx/SuL+dQZRYEOp1e1ze28tr2cl7ceY8WeKto9y9nlp8ZSnJvMpJwk0hKiSIh2EhcVQXxUBF2qVDW0UdnYRlVDG0frW9l5tJ79Vf8sFnOSYzhnXDoLx2dw9tj0Ad/x8MZwy6d1B6v55IOrWFKczW9C5CJiOLECMISsL6nm6vs+4JNzCvjx1VOH/I9FVXl56zG++/dtVDS08ekFI/nqReMDcvunor6VP6w8yJ9XldDY1smSKdncff5YpuT6NipxuL3B9mVvRSOX/HIFF07J4t4bZwUsjo4uF0+sOcT/vb6H403tXDYth7vPG8eE7MRBP/fRuhaeXFvKX9ce4khdKxmJ0Vw3J58b5xeRlxLrl3MM93xyuZRbH1nLqv3HeeaLZwbNqGCXS3lrVwWPrirhvT1VdLqUvJRYlhRnc97ETIpzk0mOG/j7VlNbJzuO1rO1rI41B6t5b08V9a2dOASmF6Rw/sRMLpicxYSsxEF5fx6O+fTbt/fyPy/v4gfLirlpfnBeRAyUy6W0d7ncj04XThESYiL8foHpKysAQ8xPX9nJvW/t498vmcjtC8cM2XnL61v51nNbeW17OVNyk/jhsqlML0gZsvP3pba5nT+sPMjDKw/Q0NrJBZOy+OK5Y5hdlOrV6w3HN9i+nGi9+cUnp7Ns5tDeCj4xaOiHy3dwoKqJ+aPS+PdLJgUkp7pcyju7K/jL6kO8udO98s75k7K4ZUERZ49N9+mDPBzy6XhjG5f8agWxkU7+fvfZAe0P2NTWydMbSnl45UEOVDWRkxzD5TNyuaQ4h2mD0O+zs8vF5tJa3tldxTu7KthcWgdAXkosF07O4oJJWcwfnea3D/7hmE/dLyKev/MsJuX4527OYOpyKfsqG9lSWse+ykaO1rZwpK6Vo3UtVNS30dbp6vW4mEgHiTGRJMVEkJUUQ05yLDnJMWQnx1CYFsfYzARykmOGrHEnaApAEVkC/BJwAr9X1R/3eF48z18CNAO3quqG/hzbm2D+gzgVl0u5+/GNLN96lPtums2S4uxBPZ+q8tT6Ur734nbaOl189aLxfPasUUNyu3cg6lo6+OP7B3novQPUtXQwoyCFz549iqXF2QN68x2Ob7B96ehy8anfr2Z9SQ33f2o2FwzRaM59lY3819+38+7uSsZmJvDNpRM5b2JmUNz+Ka1p5i+rD/HE2sNUN7UzOiOea2cXcOXMXHKSB94qGC75tOZANdc/+AGXTM0JSH/A6qZ2fr9iP39eVUJ9ayczClK47exRLBng37+vKupbeWNnBW/sKGfFniraOl0kxURw/qQsLpqcxcLxvs3nOlzzqaqxjUt+uYIIh/Cn2+YH3WwXLe1drNp/nPf2VrGltJatZfW0dHQBEOEQspJiyE1xF3TZyTHERjqJinAQHeEgKsJBl0tpaO2kobWDhtZO6lo6KK93T2dV3tD2kUUR4qOcjM1MYExmAuOzEhnn+TcvJdbvsxgERQEoIk5gN3AhUAqsBW5Q1e3d9rkEuBt3ATgf+KWqzu/Psb0J9j+IU2nt6OL6B1ex81g9T35hwaCNoCqtaebfn93Ku7srmTcqjZ9cPY1R6cE9J19vLQCfOqOIq2bl9esDfLi+wfalobWDT/1+NTuONfCHT88d1IlZG9s6+fUbe/jDygPERDj58oXjuWVBUdDdFgH339jyD4/y2OpDrC+pQQQWjB7Bspl5LCnO7ncrVzjl04kW5S8sHM03lkwckil3Khpa+f2KAzz6QQmtnV0sLc7mtrNHe30HwJ9a2rtYsaeSV7eX88aOcmqaO4iKcDB3ZCpnjknnrLHpTM1LxjmA39NwzqetZXXc+vAaulzKw5+Zx4wA32E6UNXEGzvKeWd3JasPVNPe6SI6wsHUvGSm5iczLT+ZqXkpjEqPH9D/YU9dLqWyoY2S403sqWhkr+exp6KB8vq2k/vFRTlPFoMTshOZmJ3EhOxE0hO8n+80WArABcB3VfViz/ffBFDVH3Xb5wHgbVV93PP9LuBcYOTpju1NKPxBnEplQxvLfruStk4Xz915lt/6LIG7ZeiRlQf5xeu7Afc0NJ+aXxRSc6i5XMrbuyt46L0DrNx7HBE4a0w6V83K4+Ip2X1ehQ/nN9i+1Da3c/2Dqyg53syjt81jzsg0v75+l0t5ekMpP3tlFxUNbVw7O5+vL5kYFINP+uNgVRPPbizjuU1llBxvJsrpYP7oNC6YlMV5EzNPOYo4nPLJ5VK+/fxWHlt9iE9Mz+Vn104btHk7j9W18sC7+/jL6kN0dLm4YkYedy4ew9jMwe876o3OLhfrSmp4fXs57+2tYuexBgASYyKYOzKNGQUpzChIYXpBCsmxfV9cDPd8OljVxM1/WM3xxnbu/9RsFo7PGNLz76tsZPmWo/zjw6Mn/4/GZiZw7vgMFk3IYO7INGIih24u2rrmDvZWNrC7vJHd5Q3sLm9g17EGqhrbT+6TGB1BUXocRWnxFI2IozAtjtyUWM8jhriovlucg6UAvAZYoqqf83x/MzBfVe/qts+LwI9V9T3P928A38BdAJ7y2G6vcTtwO0BhYeHskpKSQf25Btue8gau+u37ZCXHcN9NsxiX5fub3+r9x/n281vZXd7I+RMz+e7lUwZ1moyhcOID/JmNpRyubiEuysk9Syf2uhrGQN5gh1M+VTa0cd0DH1DV0Mbjt5/ht8787++t4vv/2MH2o/XMLEzhO5dNZmZh4FtnvKGqbDhUy8tbj/LGzgr2V7onpB6flcBXLpzQa3eMcMsnVeWBd/fz45d2Mm9kGg/eMtuvo2MPVzdz3zv7eGpdKV2qXDUzjy8tHhv0dyZ6qmps4/19x3l/bxXrSmrYW9F48rnRGfF8+9LJLJ6Y+bHjwiGfKupb+fTDa9lb0cDPrp3OFTPyBvV8h6ubeWHzEf6++cjJom92USqXTM3h4ilZQTlhflVjG7uOuYvBkuNNlFQ3U3K8mcPVzXS6Plq7pcZF8usbZvV6dydYCsBrgYt7FHHzVPXubvv8A/hRjwLw68Do0x3bm1C5IjqdD/Yd586/bKCxrZOvXTSe284e7VVz9JHaFn726i6e2VBGXkos3718yrCb4V9VWVdSwzMbSllSnMOiXq4uh/sV9qmU1bZw3f0fUNvczl3njeOzZ4/0ugVn25E6fvHabl7fUUFeSiz3LJ3IZdNygqKfn7+cuEX05s4KPr9wNIsn+PaB3V2o59PfNx/hq09uJj8tlkdunUfhCN8+RHcda+B3K/bz7MYynCJcOyefOxaNCfmL0xPqWzv4sLSOTYdr2XS4li+dO6bXC6Vwyaf61g4+98d1rDlQzVUz8/jXC8f79f+6oqGVlz48xvObythwqBZwF32XTcthaXEO2ckxfjvXUOrscnGsvpWjda0cqW2hrLaFI7UtfPasUYzO+Hi/ymApAO0WsA+qGtv492c+5NXt5cwpSuVn105nZD+uiFWVNQeqeeT9g7y6vRyHwO0LR3PX4nHERoXnkmvh8gbbl7LaFv7z+W28vqOckSPi+M4nJnPexP5dCLR2dPHS1qP8eZW731xidARfWjyWz5w1ckhvmwSTcM6nNQeq+fyf1tHW2cXVs/L57NmjGNPLh1Bf6lo6eGHzEZ5ad5jNpXXERDq4cV4Rty8cHbIf0L4Kp3xq7ejiF6/t5pH3D+JS5cZ5hdx13jivuo6oKrvLG3l9RzmvbS9n0+FaACZmJ3LFjDw+MT0nKFv6BluwFIARuAdynA+U4R7IcaOqbuu2z6XAXfxzEMivVHVef47tTSj+QZyKqvLcpjK+8/w2OruUS6flMKcolTkj0xiTEY+IoKpUNraxv7KJHUfr+evaw+w81kBKXCSfnFvAzWcUheUfQXfh9AZ7Ku/sruS//76NfZVNnDshg09My2V8ViJjMxNOXhx0drk4UNXEjmMNbD5cy7Mby9wjZ9PjuemMIq6Zle/VXGvDSbjn06Hjzdz71l6e3VRGe6eL8ydmctvZoyjOTyYxOuIjLcLN7Z1sP1LPh2V1rDtYw2s7ymnvdDExO5FrZuezbGZev9YUH87CMZ+O1bXyyzf28OS6w0Q5HVw5M4+ZhSlMy09mXGZir3e8Wju62HG0ns2eFtV1JTWU1rQAMD0/mQsmZXFxcTbj/dBtKpQFRQHoCeQS4P9wT+XyB1X9gYjcAaCq93umgfkNsAT3NDCfUdV1fR17uvOF8h/EqRyta+FHy3eyYk8lNc0dgLsPQH5qHAePN9HQ2nly30k5Sdx6ZhFXzMgL2xaansLxDbYvHV0u/vj+QX75xp6TeSMCBalxxEdHsK+i8eSqChEO4YJJWdy8oIgzx4wYVrd6fWH55FbZ0MafV5Xw6KoSqpvcHdijIhykx0eRnhhNS3sX+yobOdF9KSMxmiVTsrluTgHFeUmWTx7hnE/7Kxv5v9f38ObOChrb3O9HsZFOJmS7i8DWji7aOl20dnRRXt9KR5c7mTITo5lRkMK5EzI5f1ImWUnh2Xrcm6ApAIfacPiDOBVVZX9VE+sP1rCupJqjda2MSo9ndHo8YzITGJ2RQO4QTjgZKsL5DbYvHV0uSo43nRyJtqeikaa2TiZkJTIxxz0dwZiMBKIigm86l0CzfPqo1o4u3txZQVlNC1VNbVQ1tFPV2EakU5iSm3xyig37kO6d5ZN7tPn+qia2lNaypbSOXccacDggJsJJdKSDmAgnGUnRzPSMqvZm/s5w0Vc+BW6hWeMXIsKYjATGZCRw3dyCQIdjQlik08HYzETGZiZyydScQIdjQlhMpNNyyPjE4RDGZiYwNjOBq2YN7QpG4cIu5Y0xxhhjwowVgMYYY4wxYcYKQGOMMcaYMDPsB4GISCXQfWr0dKAqQOF0FwxxBEMMEJg4ilR1wOsPWT4FfQxg+eQPwRBHMMQAlk/+EAxxBEMMEET5NOwLwJ5EZJ03o6uGYxzBEEMwxeGNYIk9GOIIhhiCKQ5vBEvswRBHMMQQTHF4I1hiD4Y4giGGYIoD7BawMcYYY0zYsQLQGGOMMSbMhGMB+GCgA/AIhjiCIQYInji8ESyxB0McwRADBE8c3giW2IMhjmCIAYInDm8ES+zBEEcwxADBE0f49QE0xhhjjAl34dgCaIwxxhgT1qwANMYYY4wJM1YAGmOMMcaEGSsAjTHGGGPCjBWAxhhjjDFhxgpAY4wxxpgwYwWgMcYYY0yYsQLQGGOMMSbMWAFojDHGGBNmrAA0xhhjjAkzVgAaY4wxxoQZKwC9JCIHReSCAJz3RyLyZR+OXyMiU/wYkvEDyyfjT5ZPxp8sn4YnKwBDiIhkALcAD/jwMj8D/ts/EfWPiIwTkVYR+fMp9vmziBwVkXoR2S0inxvKGMNRqOXTQHJERCaJyJsiUicie0Vk2VDEGM6Gcz51O+a072XGP4ZzPolImog8KyJNIlIiIjcORYw9WQEYWm4Flqtqiw+v8QKwWERy/BNSv9wLrD3NPj8CRqpqEnA58H0RmT3okYW3WwmtfOpXjohIBPA88CKQBtwO/FlExg9BjOHsVoZhPvXQn/cy4x+3Mnzz6V6gHcgCbgLuC0RLpRWAfuBpbXhbRGpFZJuIXN7tuVkislFEGkTkbyLyVxH5vpenWgq80+PcO0SkUUTaPY9Gz2NSby+gqq3AeuAiL2MYEBG5HqgF3jjVfqq6TVXbTnzreYwZ3OiCk+VT7waQIxOBXOAXqtqlqm8CK4GbBzvGYGT51LuBvuf0971suLN86l1/80lE4oGrgW+raqOqvoe7UB3y9ycrAH0kIpHA34FXgUzgbuAxEZkgIlHAs8AjuFsiHgd8uRU1FdjVfYOqTlLVBOBR4L9UNcHz2HGK19kBTO/j53nR84fd2+PFgQQrIkm4m9+/2s/9fysizcBO4CiwfCDnGw4sn06tnzkifWwrHuj5Qp3l06n19z1noO9lw5Xl06n1M5/GA12qurvbts2AtQCGoDOABODHqtruaW14EbjB81wE8CtV7VDVZ4A1Jw4UkWRxd1JtFJHibtt/ICIrROQpEYnrdq4UoKGPOKYBW/sZc4PntT5GVS9T1ZQ+Hpf18/VP+B7wkKoe7s/OqvolIBE4B3gGaDv1EcOS5dMp9DNHdgIVwL+JSKSIXAQsAuJ62Xe4s3w6hQG85wzovWwYs3w6hX7mUwJQ12Nbnee4IWUFoO9ygcOq6uq2rQTI8zxXpqra7bnubyDNwKXAUyc2eP4wxqjqOcDrwGe77V9DL0kiIg5gMv3/g0jEfSvDL0Tkpm5N8S95ts0ALgB+MZDX8tyyew/IB77orxhDiOVTL/nU3elyRFU7gCtx/y6O4W61eRIo9VeMIcTyycd88va9bJiyfPIxn4BGIKnHtiT6LnYHjRWAvjsCFHiS8oRCoAx3E3CeiHS/JVVw4gvPVVJlj9c7BziRVC8BZ3d7bgvu5uOeCnH/X+7vvlFEPiki74rIShG5rttTk3A3OX+MiLzULbl7Pj6W7J6f47FuTfFLPZvPBUYCh0TkGPA14GoR2dDba/QigvDsA2j51Hs+9abPHFHVLaq6SFVHqOrFwGi6tUaEEcsn3/PpXHx7LxtOLJ98z6fdQISIjOu2bTqw7RSvNSisAPTdaqAJ+LrndtO5wCeAJ4APgC7gLhGJEJErgHmneb1U/tk8XIe7L8UJy3HfyuopyRND1IkNnqvWW4HzVfUs3H0zEJFoYDbwWm8nV9Wl3ZK75+NUyd7Tg7iTf4bncT/wD+DinjuKSKaIXC8iCSLiFJGLcd9SeHMA5xsuLJ96MdAcEZFpIhIjInEi8jUgB3ffpHBj+dSLAeZTv9/LwoDlUy8Gkk+q2oT79vB/i0i8iJwFXIG7X+OQsgLQR6rajnvI91KgCvgtcIuq7vQ8dxVwG+4m6E/h7i9xqr5tNUCy5+tkoLrbc38CLhGR2B7H7MB9hVMjIhM92y4D/tdzO+zEbTE8sb6tqkcG+KMOiKo2q+qxEw/czd6tJ64APVde/35id9xN5aW4f/6fAV9W1ecHM8ZgZPnUp1PmSI98AveIuqO4+wKeD1yo/xyhFzYsn/rU73w63XtZOLF86tNA35++BMTifn96HPiiqg55CyCqao8hfOC+gvpMj22PAMWer6cCf/F8fTtwd499f4g7sU53np8DF3u+juhx/uJA/x7sYflkj+B7WD7Zw58Py6fgfojnl2QGiYgswj2UvQr3hI/3A6NV9ajn+eW4byuUAA+o6iMi8iPcfSEqcF9dNXlx3knAw7ivvtao6r/54ccxAWb5ZPzJ8sn4k+VTaLECcJCJyO24pxBIAPYB31TVfwQ2KhOqLJ+MP1k+GX+yfAotQ14AisgfcN+vr1DVj03M6hlB9EvgEtzDxm9V1Q2e55Z4nnMCv1fVHw9Z4MYYY4wxw0QgBoE8Aiw5xfNLgXGex+3AfQAi4sS9ft5S3HMA3SAikwc1UmOMMcaYYShiqE+oqu+KyMhT7HIF8Cd1N02uEpEUcS/kPBLYq6r7AUTkCc++2091vvT0dB058lSnM+Fo/fr1VaqaMdDjLJ9MbyyfjD9ZPhl/6iufhrwA7Ic8Pjp7eKlnW2/b5/f2Ap5+CLcDFBYWsm7dusGJ1IQsESkZwL6WT+aULJ+MP1k+GX/qK5+CcR7A3hZy11Ns//hG1QdVdY6qzsnIGPBFlDEfYflk/MnyyfiT5ZPxVjC2AJbSbfkY3OvpHcE963dv240xxhhjzAAEYwvgC8At4nYGUOeZQ2gtME5ERolIFHC9Z98Be2d3Jd99Yegn3TbD0yvbjvGTl3cGOgwzTDy/qYxfvr4n0GGYYeKJNYd48N19gQ7DBKEhLwBF5HHcawZOEJFSEblNRO4QkTs8uyzHvcjzXuB3uJdMQVU7gbuAV3AvBfOkerl0yp7yBh55/yAHqwY836QxH7OtrI4H3tlHXXPH6Xc25jTWl9TwwLv7aOvsCnQoZhh4f99xHnz3AC6XzflrPmrIC0BVvUFVc1Q1UlXzVfUhVb1fVe/3PK+qeqeqjlHVqaq6rtuxy1V1vOe5H3gbw5LibABe3nbM55/HmEUTMnEprNgbdkuDmkGwaHwGze1drDtYE+hQzDCwaHwGVY1tbD9aH+hQTJAJxlvAgy4/NY6pecm8tNUKQOO7GQUppMRF8tZOKwCN7xaMGUGU08HbuyoCHYoZBs4Znw64uz4Z011YFoDgbgXcfLiWI7UtgQ7FhDinQ1g0PoN3dlfYbRbjs7ioCOaNSrMPbOMXmYkxTMlN4p1dlk/mo8K6AAR41W4DGz9YPCGTqsZ2th6pC3QoZhhYND6D3eWNdoFq/OLcCRmsP1RDXYv1Uzb/FLYF4JiMBMZnJdhtYOMXC8dnIILdBjZ+ce4E93xu1gpo/GHR+Ey6XMr7e6sCHYoJImFbAAIsmZLN2oPVVDW2BToUE+LS4qOYnp/CW9Zvy/jB2MwEcpNj7Lad8YuZhSkkRkfYBYX5iPAuAItzcCm8tr080KGYYWDxhEw2l9ZS3dQe6FBMiBMRFk3IZOXeKjq6XIEOx4S4SKeDs8el887uSlStn7JxC+sCcFJOIoVpcbxst4GNHyyemIEqvGtX2cYPFo3PoKGtkw0lNh2M8d2i8RkcrWtld3ljoEMxQSKsC0ARYWlxNu/vq7LOscZnxbnJpCdE2W1g4xdnjR1BhEPstp3xi0Un+5Xa+5NxC+sCEODi4mw6upQ3d9ptYOMbh0NYND6Td3ZX0mXTwRgfJcZEMrsolbetH6Dxg5zkWCZkJdoFhTkp7AvAGfkpZCfF2G1g4xfnTsigtrmDTYdrAx2KGQbOnZDJ9qP1VNS3BjoUMwwsmpDB2gM1NLV1BjoUEwTCvgB0OISLp2Txzu5Kmtvtj8L4ZuG4DByCreJg/GLReJsOxvjPovEZtHe5+GDf8UCHYoJA2BeA4B4N3NrhslstxmfJce7bdtYP0PjDpJxEMhOjrQA0fjFnZCpxUU7LJwNYAQjA3JGpJMVE2OhN4xfnTshka1k9FQ122874RsS9zOCKPVV02nQwxkfREU7OHDOCt3dX2HQwxgpAgAing3mj0lh9oDrQoZhh4MQqDu/utln3je8WTcigrqWDzaW1gQ7FDAOLxmdwuLqFA1VNgQ7FBJgVgB7zR43gQFWTdbY2PpuUnURiTAQbDtn8bcZ3Z45JB2DdQcsn47szx1o+GTcrAD3mj04DsFZA4zOHQ5iWn8wWa7ExfpAWH0VeSiwfltUFOhQzDIwaEU9CdARbj1g+hbshLwBFZImI7BKRvSJyTy/P/5uIbPI8topIl4ikeZ47KCIfep5b58+4JuckkRAdweoDNjrK+G56fgo7jzbQ2tEV6FDMMDA1L5mtVgAaP3A4hMm5SZZPZmgLQBFxAvcCS4HJwA0iMrn7Pqr6U1WdoaozgG8C76hq92a5xZ7n5/gztging9lFqazeby2AxnfT8lPodCnbj9YHOhQzDEzNT+bg8WbqW23FIuO74txkth+tt4FFYW6oWwDnAXtVdb+qtgNPAFecYv8bgMeHJDLct4H3VDRyvLFtqE5phqnpBckAbLEJoY0fFOe588labYw/FOcl0drhYr8NBAlrQ10A5gGHu31f6tn2MSISBywBnu62WYFXRWS9iNze10lE5HYRWSci6yor+z+1y/xRIwBYe9BaAY1vspNiyEiMZkupfWAb3021AtD0wZvPO8snA0NfAEov2/qajOgTwMoet3/PUtVZuG8h3ykiC3s7UFUfVNU5qjonIyOj38FNzUsmNtLJKrsNbLrx5g1WRJien2xTd5iP8Saf/jkQxLoUmI/y5vNudEYCMZEOtlo+hbWhLgBLgYJu3+cDR/rY93p63P5V1SOefyuAZ3HfUvabqAgHs4pSbCSw+QhvLyim56ewv6rJ+m2Zj/A2n4rzrOO+8Q+nQ5icY/kU7oa6AFwLjBORUSIShbvIe6HnTiKSDCwCnu+2LV5EEk98DVwEbPV3gPNHjWDnsXrqmu1D2/hmWkEKqrDVbgMbP5ial8wBu6AwflKcl8y2I3W4XLYiSLga0gJQVTuBu4BXgB3Ak6q6TUTuEJE7uu26DHhVVbv3UM0C3hORzcAa4B+q+rK/Y5w/Kg1V6wdofDfN089msxWAxg9ODATZZrftjB8U5yXT1N7FweM2ECRcRQz1CVV1ObC8x7b7e3z/CPBIj237gemDHB7TC1KIinCw+sBxLpicNdinM8NYanwUhWlxNiG08YvuHfcXjBkR4GhMqCvOdefTh2V1jM5ICHA0JhBsJZAeYiKdzCiwfoDGP6YXpNhIYOMXIxKiyU2OsRVBjF+My0ogKsLBtiPWohyurADsxRmj0thaVkeD9bUxPpqen0xZbQuVDTa3pPFdsa0IYvwk0ulgUnai5VMYswKwF/NHj8ClsK7EFss2vpmWnwJgt4GNX0zNS2Z/VZNdnBq/mOK5oFC1gSDhyArAXswqTCXSKayx28DGR8V5STjEBoIY/zg5EMRu2xk/KM5Npr61k8PVLYEOxQSAFYC9iI1yMi0/hdX7jwc6FBPi4qIiGJeZaC2Axi9sSTjjTycHFh2xfApHVgD2Yd6oNLaU1tHS3hXoUEyIm16QzJZSu81ifJeRGE12kg0EMf4xPjuBCIfYBUWY8qoAFJE5nomch625I1PpdCmbDtcGOhQT4qblp1Dd1E5pjd1mMb4rzku2AtD4RXSEk/FZiZZPYWrABaCI5ADvA9f5P5zgMbswDYD1JdYP0PhmumcgiK0LbPzhxIogjW2dgQ7FDANT85LZdqTe7lCEIW9aAD8N/BH4nJ9jCSrJcZGMz0pg7UEbCWx8MyE7kSinw+YDNH4xNT8JVdhuA0GMHxTnJVHd1M7RutZAh2KGmDcF4M3AN4EoERnj53iCypyRaWw4VEOXrZVofBAV4WBybhKbrTuB8YMTA0Hstp3xhymWT2FrQAWgiCwGdqpqFfAwcNugRBUk5hSl0tDaye7yhkCHYkLc9Hz3fFu28LrxVWZiDFlJ0dZx3/jF5JwknA5hm+VT2BloC+BtwEOer/8KXCsiw3Yk8dyR7n6ANiG08dUUW3jd+NFUGwhi/CQm0snYjATLpzDU7+JNRFKAM4CXAFS1HlgFXDIokQWB/NRYMhOjWXfQBoIY35xYeH2r9dsyfjA5J4n9lY20dtg0VcZ3U3KT2HHU7nSFm34XgKpaq6pj9aNDhW5V1RcHIa6gICLMHZnGOhsIYnw0LiuBKKeDbTbhqvGDSTlJuBR2HbMPbeO7iTmJHKtvpaapPdChmCHk6+3bB0QkDkBEFvohnqAzuyiVstoWjtbZHG7Ge5FOBxOyE9lWZi2AxneTcpIA2HnM8sn4bmL2iXyyC4pw4msB+J/AQyLyKDDXD/EEnZP9AK0V0PhoSm4SW4/YiiDGd4VpccRHOe22nfGLiTmJgF1QhBtfC8DvAbsABZ7szwEiskREdonIXhG5p5fnzxWROhHZ5Hl8p7/HDoZJOYnERTmtH6Dx2ZS8ZGqbOzhi820ZHzkcwoTsRLYftQ9s47uMhGhGxEexw/IprET4ePzXVbVKROKBX3KayaFFxAncC1wIlAJrReQFVd3eY9cVqnqZl8f6VYTTwczCFBsJbHxWnOu+zbK1rI68lNgAR2NC3aScJF7YfARVRUQCHY4JYSLCpJwkuwUcZvrVAigihb1t98wHiKo2AV/ox0vNA/aq6n5VbQeeAK7oZ6y+HOuT2UVp7Dhab0svGZ9MzE7CIbDNRgIbP5iUk0RDaydltdY/2fhuYnYiu4412MIHYaS/t4BfFpEqEVkhIr8VkTtE5CwRSTqxg6r2Zz6CPOBwt+9LPdt6WiAim0XkJRGZMsBjEZHbRWSdiKyrrKzsR1inNndkKi6FjYesFdB4LzbKydjMBJtw1fjFiYEg1g/Q+MPEnCTaOl02V2kY6VcBqKqTgVzg/+Ge+28s8G1gp4gcGMD5ertP0fNyYwNQpKrTgV8Dzw3g2BPxPqiqc1R1TkZGxgDC693MwlQcgq0LbHxWnJvMVpsKJqz56wJ1YnYiIli/rTDnz3wC2GkXFGFjIPMAtqvqRuBZYDVwDGgBNg/gfKVAQbfv84EjPc5Tr6qNnq+XA5Eikt6fYwdLQnQEk3KSWF9iA0HCkT9blCfnJlFe30ZlQ5ufojOhxl8XqPHRERSlxVkBGOb8lU9jMxNwOsRGAoeR/vYBnCAiXxGRN4H3gQXAY8AkVb1yAOdbC4wTkVEiEgVcD7zQ41zZ4unRLCLzPDEe78+xg2lOUSobD9XS2eUaqlOaIOHPFuViz8LrNiG08YdJOUlWABq/iIl0Mjo93vIpjPS3BXAHcBNwHzBHVb+iqq95BmP0m6p2AncBr3he80lV3ebpU3iHZ7drgK0ishn4FXC9uvV67EDO74vZI9Nobu+y/jbGJ5M9I4FtIIjxh0k5SZRUN9NkA9SMH0zMsSXhwkl/p4H5IjAVdwF2r4gcBz488VDV5/p7Qs9t3eU9tt3f7evfAL/p77FDZe7IVADWHKxman5yIEIww0BSTCRFI+KsBdD4xaScJFTdKzjMLkoNdDgmxE3KSeTvm49Q39pBUkxkoMMxg6y/g0AeUNW7VHWRqmbinovvEaAduHoQ4wsaOcmxFI2IY+XeqkCHYkJccW4yW21JOOMHkzwrONhtO+MPkzxLwtka0+HBq5VAVLVUVZer6k9U9WZ/BxWsFo7L4IN9x2nr7M+MN8b0bnJuEoeqm6lr6Qh0KCbE5aXEkhQTYQWg8YuTS8JZPoWF0xaAIrJQREaLyJ9F5EkRWTgUgQWjheMzaOnoYr2tCmJ8cGIgyHbrB2h8JCKefluWS8Z32UkxJMdGst36AYaF/rQA3gB8C/gK8Cnc/QHD0oIxI4hwCO/uttvAxntTTg4EsX6AxneTPUt4uWwFB+MjEWFidqJNBRMm+lMATgGyVLXCM+o3bD+1EqIjmF2Uyru7fV9dxISv9IRospNibCSw8YtJOYk0t3dxqLo50KGYYWBSThK77IIiLPSnAPw28JNu378ySLGEhIXjM9h+tN4m8jU+Kc5LYqstCWf84J9LwtkFhfHdxGz3BcXhGrugGO5OWwCq6juq+m63TeMHMZ6gt3CceyLgFXusFdB4b0puMvsqG2lptwFFxjfjsxJx2JJwxk9sjenw0Z9BIE92e/wN+NwQxBW0puQmMSI+ihV7rB+g8d6U3CRcCjusr43xUUykk9EZCdZx3/jF+Cz3GtPWD3D4689E0PWqerLoE5H7BjGeoOdwCGePS2fFnkpcLsXhkECHZELQiZHAH5bWMavQJvA1vpmUk8QGm53A+EFslJNRI2xJuHDQnz6AP+jx/X8MRiChZOG4DKoa29lufyDGSznJMeQkx7D2YHWgQzHDwKScRMpqW2xuSeMXE3MS2WmTQQ97/ekDeABAROJEZLqqnvzEEpFCEckbzACD0Tnj0gF41/oBGi+JCHNHprH2YDWqNtrO+ObECg42ga/xh4nZSZQctzWmh7uBrATSATwjIvHdtv0eyPFvSMEvMymGidmJNh2M8cncUWmU17fZ9B3GZ1Py3AXghzay3PjBxGzPiiDWCjis9bsAVNUO4Fngk+Bu/QMyVHXdIMUW1BaNz2B9SY1dIRmvzR+VBsCaA3Yb2PgmMzGGgrRYW6XI+MXUfHcf5S2ltYENxAyqga4F/HvgM56vbwEe9m84oWPh+Aw6upRV+48HOhQTosZmJJASF2n9AI1fzC5MZV1JjXUpMD7LSY4lJznGLiiGuQEVgKq6E0BExuNeIu7RwQgqFMwZmUpspNNuAxuvORzCnKI0awE0fjF7ZBqVDW2U1rQEOhQzDMwqSmXjodpAh2EG0UBbAAEewt0SuEVVw/byIDrCyRmj03jX5gM0Ppg/Ko2Dx5upaGgNdCgmxM32TCdkrTbGH2YVplJW28KxOntvGq68KQCfBKbjLgQHTESWiMguEdkrIvf08vxNIrLF83hfRKZ3e+6giHwoIptEJOB9D8+blMWBqiabL8l4ba6nH+DaA/ahbXwzITuRhOgI1pVYi7Lx3ewi9wXFhkP23jRcDbgAVNVmVU1W1dcHeqyIOIF7gaXAZOAGEZncY7cDwCJVnQZ8D3iwx/OLVXWGqs4Z6Pn97dKpOUQ4hGc3lgU6FBOipuQmERvpZM0B60tqfON0CDMLU1hfUhvoUMwwMDkniegIh7UoD2PetAD6Yh6wV1X3q2o78ARwRfcdVPX9breWVwH5Qxxjv6XFR3HuhEye31RGl8s6XpuBi3Q6mF2UypqD9iZrfDerMJVdx+ppaLUJoY1voiIcTMtPthbAYWyoC8A84HC370s92/pyG/BSt+8VeFVE1ovI7X0dJCK3i8g6EVlXWTm4gzSWzcyjvL6N9/dZX0Djnbkj09h5rN5WcTA+mzMyFZfCpsO1gQ7FDAOzClPZWlZHa0dXoEMxg2CoC8DeFs7ttelMRBbjLgC/0W3zWao6C/ct5DtFZGFvx6rqg6o6R1XnZGRk+BrzKZ0/KZPEmAie3WC3gY135o5KRRXWW98t46MZBSk4xAaChJPBbPCYVZRKR5ey7YhNMD4cDXUBWAoUdPs+HzjScycRmYZ7pPEVqnqyc5SqHvH8W4F7Uup5gxptP8REOrl0ag4vbztGc7tNCm0GbmZBKpFOYY0NBAkLg/mBnRgTyYTsJCsAw8hgNnjM8ows32D9SoeloS4A1wLjRGSUiEQB1wMvdN/Bs8LIM8DNqrq72/Z4EUk88TVwEbB1yCI/hWUz82hu7+KVbccCHYoZBIPdpSA2ysnUvGSbEDpMDPYditlFKWw8VGv9ko3PMhKjKUyLswuKYWpIC0BV7QTuAl4BdgBPquo2EblDRO7w7PYdYATw2x7TvWQB74nIZmAN8A9VfXko4+/L3JFp5KXE8ozdBh6WhqJLwdxRaWwprbW+NsZnc4rSaGzrZJet42r8YFZhCusP2Qozw9FQtwCiqstVdbyqjlHVH3i23a+q93u+/pyqpnqmejk53Ytn5PB0z2PKiWODgcMhLJuZx8q9VVTU26SZZuDmjUyjo0tt5n3jsxPzt6230ZvGD2YVpdoKM8PUkBeAw9WyWXm4FJ7f9LEujcac1pyiNESw28DGZ/mpsWQmRrPecsn4wcl+gHZBMexYAegnYzISmJ6fzDM2KbTxQnJcJBOyEm1dYOMzEWF2Uaq1ABq/mJidSFyUkw3WD3DYsQLQj66cmceOo/XsPGZLw5mBO2P0CNaVVNPYZqPJjW9mF6VyuLrFuqQYn0U4HUzPT2GDdU8ZdqwA9KNPTM8lwiH86YOSQIdiQtBl03Jo7XDxylYbTW58c7IfoLXaGD+YVZTC9qP1NtXZMGMFoB+lJ0Rz4/xC/rr2MPsqGwMdjgkxs4tSKUiLtbWljc+m5CbbOq7Gb2YXpdLlUraUDs2E0KU1zZRb6/WgswLQz/7f+eOIiXDw05d3BToUE2JEhGUz8li5r4pjdfbmZ7wXFeG+bWeDiow/zCwY/BblA1VN3PvWXi779QrO/slbnP2TN/nuC9s43tg2aOcMd1YA+ll6QjS3LxzDy9uO2dW3GbBls/JRhRc2Wyug8c2iCRlsLq3jcHVzoEMxIS41PorRGfGDMhDkYFUTl/16BYt/9jY/fWUXEQ4H31w6kWtmF/DoqhIW/fRtfv3GHrv9PAisABwEnztnFOkJ0fz4pR02eaYZkFHp8cwoSLFJxY3Pls3MQwSe3lAa6FDMMHDO2HRW7K2irrnDb69Z2dDGLX9YQ1lNC9++bDIr7zmP5+48iy8sGsOPrprKK19eyJljRvDz13Zz7k/fZtPhWr+d21gBOCjioyP48gXjWHuwhtd3VAQ6HBNils3MY+exBnYctdHkxnu5KbGcNSadpzeU4rJl4YyPrp1TQHunixe2+Geu26a2Tm7741oqGlr5w61zue3sUeSlxH5kn7GZCTx4yxyeumMBMZFObnloNVvLhqYfYjiwAnCQfHJuAaPT4/nJyzvp7HIFOhwTQk6MJrfBIMZXV8/O43B1y6D3BSyvb+VIra0UMZxNyU1iYnYiT6077PNrdXS5uPMvG9haVsdvbpjFTM9k032ZMzKNv3x+Pokxkdz80Gp2l9syh/5gBeAgiXQ6+PqSieytaOSp9XYLxvRfWnwU507I4PlNZXRZy43xwcVTskmIjvD7e1BFfSvPbyrjm89sYfHP3mb+D9/gzB+/ydX3vc+fV5VQ29zu1/OZwBMRrp1TwObSOp8KMFXl35/5kLd3VfL9K6dyweSsfh2XnxrHY5+bT6TTwY2/W81+m2nDZ1YADqKLp2QxuyiVn7+224a0mwFZNjOf8vo2Pth3PNChmBAWFxXBpVNzWP7hUb91on/0g4Oc+eM3+ZcnNvHilqOMyYjnPy6ZxNeXTKC+pYNvPbeVeT94gy88uo491lIzrFw5w3134m8+tAL++s29/G19Kf/v/HHcOL9wQMeOTI/nL5+fj6py0+9X2wAnH1kBOIhEhO9dUUxzWyeffWStrfBg+u38SZkkRkfwzEZrPTa+uXp2Pk3tXbzs4wTjHV0uvvXch3z7+W0sGp/Bi3efzabvXMTvPz2Xzy8czZfOHcur/7qQF+8+m0+dUcTqA9Vc9uv3eHjlAeuDOEyMSIjm/EmZPLuxjA4vuja9v6+KX7y+m2Uz8/jXC8Z5FcPYzEQevW0+ze1d3Pj7VVQ2DN00MXUtHaw7WM1fVh/if17eybMbSyk53hSygz0jAh3AcDc5N4nf3DSLz/1xHXc+toGHPj2HCKfV3ebUYiKdXDI1h79vOcL3r+wkLsr+VI135o5MpTAtjqfWl3LVrHyvXqOmqZ0vPbaBD/Yf545FY/i3iyfgdMjH9hMRivOSKc5L5ovnjuEbT2/hv/6+nTd3VvCza6eTlRTj649jAuza2QW8sq2ct3ZWcNGU7H4fV93Uzr/+dROj0uP5wbJiRD6eP/01OTeJRz4zlxt/t5rPPLKGJ25fQEL04LxHHq5u5ldv7OHdPZWU1/+z2BSBE3VfekIUMwtTWTwhk6tm5RET6RyUWPqyt6KRpzeU8oWFo0mJi+r3cVaJDIHFEzL5/pXFvLO7km89tzVkrxbM0Fo2K4/m9i6etj6kxgciwtWz8vlg/3FKawZ+y2xvRQNX/nYl60tq+N/rpnPP0om9Fn89ZSRG89Cn5/CDZcWsO1jDxf/3Li99eNSbH8EEkXMnZJCeEM3fBvC+pKr82982U9PUwa9vmOmXC9qZhance9NMdhxt4It/Xk97p38HW1bUt/Kd57dy3s/f5vnNRzhj9AjuWTqRP9w6hxVfX8ye7y/lpX85hx8sK2bh+Ax2lzfw789+yNk/eYvfvr2X+lb/TZfTm7rmDh5dVcKV967kgv99hwff3c+6gwObp9GaFYbIDfMKKatp4Tdv7aUgLY47F48NdEgmyM0bmcYZo9P4wfIdzBmZxqScpECHZELUVbPy+MXru3l2Qxl3n9//W2/bj9TzqYdW4xB4/PYzTq4x3F8iwk3zi1gwegT/+tdNfPGxDXzqjEK+denkIW8lMf4R4XRw9aw8HnrvAFWNbaQnRJ/2mIdXHuSNnRV89xOTmZKb7LdYzpuYxY+WTeXrT2/hnqe38PPrpvvUsgjQ2NbJvW/t5eGVB+joUj45t4C7zxtLTnLsx/adlJPEpJwkbppfhKrywf7j3P/Ofv7n5V389q193HRGIbedPYrMRP+1fG87Usfv3t3P8q3HaO90MSErkf+4ZBJXzMwd8HmGvAAUkSXALwEn8HtV/XGP58Xz/CVAM3Crqm7oz7HB7qsXjaestoWfvrILVeULi8YQGYK3g1s7uqhuaqe2uYOWji7aOrpo7eyircOFS91N4yf+BCOcDuKjnMRFR5z8NzEmgoSoCBz9aEUIZw6H8KsbZnLZr97jjj+v54W7ziY5NjLQYZkQVJAWxxmj03h6Qyl3nTe2Xx+SW8vq+NRDq4mNdPL4589gZHq81+cfnZHAU188k5+9sosHPC0V9940izEZCV6/5qlUNbaxu7yB6qZ2jje2c7ypnbrmdpLjoshJjiE7OYac5BgKUuOIH6Rbh8PZtXPyeeDd/Ty3sYzPnTP6lPtuLavjRy/t4IJJWXz6zJF+j+W6uQWU17fy89d2k5kUwz1LJ3r9Wq9sO8Z/Pr+N8oZWLp+ey79eML7feS8inDkmnTPHpLO1rI7739nH797dz8MrD/LJOQV8YdFo8lPjvIpLVXl3TxW/e3c/7+2tIj7KyQ1zC7h2TgFTcpO8LnqHNPNFxAncC1wIlAJrReQFVd3ebbelwDjPYz5wHzC/n8cGNRHhJ1dPo6PLxc9e3c3yD4/xP9dMozjPf1dEvlJVKhra2F/ZxOGaZkprWiitdv9b2dhGVUMbDX4YzOIQSIyJJDk2kpS4SEbER5GeEM2IhGjSE6IYkRBFalwUI+KjSY2PZER8NLFR4ddikJkYw29vmsX1D67iq09u5sGbZ1vhbLxyzewCvva3zbyyrZwlxafuu7X5cC03P7SaxJhInrj9DArSvPvg6i7S6eCbl0zijNEj+MqTm/jEr9/j+1cWe90vsbvWji7WHqzmvT1VvLun6mOTqItAYnQEDW2ddO+B4xCYmJ3E7KLUk4/81FifW5F6crmUI3UtHKhq4kBVE/srm9hf1cTd541l7sg0v55rKIzNTGRGQQp/W1fKbWeP6vP3Vd3Uzt2Pb2REfDQ/vWaa33+vJ9x13liO1bdy/zv7KK9v5TuXTSY1vv994cpqW/jP57fx+o5yJmYn8ttPzWLWaeYmPJXivGR+c+MsDlQ18cA7+3hi7SEeX3OIK2bkceuZI5mSm9Sv9/G9FY28su0YL2w6wq7yBjITo/nGkoncOL/QL40BMpT90URkAfBdVb3Y8/03AVT1R932eQB4W1Uf93y/CzgXGHm6Y3szZ84cXbdund9/Fl+9vPUY33l+K8eb2vnc2aP48gXjh7TAUVUqG9rYeayBncfq2VPeyN7KRvZWNNLQ+s8CTwRykmLIT40jKznGU6i5i7WUuEhiIp0nH9ERDhwiKOo5B3R2Kc3tnTS1d9LU1kVTWycNrZ3UtXRQ39pBXUsHNc0dVDe1UdXQzvGmNjq6es/JmEgHaXFRpHmKw+TYSBJjIkiMiSQhOoKE6AjiPK2McZFO4qKdjBwRT27Kx5vuRWS9qs4Z6O8tUPn0h/cO8N8vbufrSybwpXP7331AVdlb0cg7uyspOd7MkdoWjtS1crSuhZb2LjISo8lMjCYzMYaspGgm5iQxPT+F8VkJNlhpAEIhn1rau1j225XsrWjkh8umct3cgl73W19Sw61/WENqfBR/+fx8r1stTuVoXQv/8vgm1hysZlZhCp8/ZzQXTcnuV9/CE+qaO3hjZzkvbz3Gu3sqae1wEekUZhelcs64DGYUpJCRGE1avPv9wukQOrpcVDS0cayuhaN1rewub2RDSQ0bD9XQ1N4FuPsuzi50F4OzilKYkps8oNvV9a0dlFQ1s+NoPduO1LHtSD07jtaffH2A+CgnozMS+NrFE1g0PuNjrxEK+fTY6hL+49mt/NflU7hlQdHHirudx+r53B/XUdHQxp8+O48zRo8Y1Hi6XMqv3tjDvW/tJSUuku9dUczSqTmnPKamqZ0/ryrhvnf2oQpfvmAcnz17lN/vzB2ta+F37x7gL2tKaO1wkRYfxZljRnDOuHTmj3L/XhrbOqlv7aChtZPNh2t5Zdsx9lU2ATA9P5lPnVHE5TNyiY4YeJ3QVz4NdQF4DbBEVT/n+f5mYL6q3tVtnxeBH6vqe57v3wC+gbsAPOWx3V7jduB2gMLCwtklJSWD+nN5q66lgx+/tIPH1xwmOymGS6flsLQ4m1mFqX5t5alv7WBPeQO7jjWyu7yB3eUN7Dzmvj1yQnpCNGMz4xmXmcjYzARGZ8RTmBZHTnIsURFDVwioKvUtnRxvaqOmuZ3qJndxeLypnZqm9pP/Vje109DaSX1rJ41tHbR29N4B+KsXju+1z1MovMF2p6rc/fhGln94lEdvm89ZY9P73Lejy8W6gzW8vqOc13eUU3Lc3fE/OTaSnOQYclNiyUmOIS7KSWVDGxWeR3ld68nW3dhIJ1PzkplekMyMglRmFKaQmxwzaFfwoS5U8qmhtYM7/7KRd3dXctfisXz1ovEn/0/L61v5w3sH+NMHJWQnx/CXz8/vtd+Tv3R2uXhs9SEeeu8Ah6qbKUiL5bNnjeKqmfkkx328daO1o4utZXVsOlzLu3uqeH9vFZ0uJTsphoumZLF4QibzRqV5dUu3y6XsOtbA+pJqNhyqZX1JDYe6zTGXmRhNfmosBWlx5KbE4hD3xW17l4uOLhe1zR0cqm7mUHUztd3Wyo2PcjIpJ4kpuUlMyE5idEY8o9PjyUiMPuXfUijkU3N7J7f/aT3v7a3inHHp/PjqaSeXcntl2zH+9a+bSIiO4MFb5jCjIGVIYgJ3v9WvP72ZrWX1XDI1m69dNIGCtLiPFHV7Kxp46L2DPLuxlNYOFxdMyuI/PzHZLy3dp1Ld1M5bOytYubeKFXur+pzCxukQzhidxsVTsrlgUlavjRgDESwF4LXAxT2KuHmqene3ff4B/KhHAfh1YPTpju1NsLYAdrdq/3EeeGcfK/cep73LRWZiNBdPyWZafjIFaXEUpMWRnRTT69WxqtLU3kV1o7v1rLy+jYPHmzjoudVw8HjTR4aux0U5GZeZwMTsJCZkJzIxJ5GJ2UmkDaC5PBi1d7poauukuaOLFk9rY3N718k37Z5C4Q22p6a2Tq64dyX7KhuZkJXIGaNHMH9UGrOKUik53szq/cdZfaCa9SU1tHR0EeV0sGDMCC6YnMX5EzNP+yaiqhyqbmbT4dqTj21H6k+OrstIjGZGQQqj0+PJS40lLyWWvNRYMhKiiYpwEBXhINLhwOEQ2jq7aGjt9Dw6qG/ppKa5nVpPUV/b0k59y4krXvdVb0tHF6rgUsWliipEOR1Ee1qXoyMcxEU5SYp1dx1IiokkKdbdApwUc6I12P19fLSTuCh3v1NfWjJVlS6X0uly/xvpdPR6QRRK+dTR5eLbz23libWHuXx6LncsGsMj7x/g2Y3ulWcunZbLty+dROYQTdnS5VJe236M3604wPoS9yjGxJgIcpJjyEmOJS0+ij0VDew82kCnZz7BUenxXDQliyVTspmenzIo3SIqGlrZUFLL7vIGDnu6wRyuaeZonXtS/0inEOlwEBnhIDEmgsK0OArT4iga4f53QnYSRWlxXsU2kHwKZIOHqvLY6kP8cPkOHCL8+yWTON7Yxs9f2830/GQevGVOQKb+6exy8eCK/fzfa3to73IhAhkJ0eSkxBLlFNYerCEqwsFVM/P4zFmjmJCdOOQxqiq7yxvZeKiGSKeDBM/7V1JMJAWpcb1eBHkrWApAuwV8CvWtHby1s4KXtx7j7V2VtHT885ZBpFMYEe8ebaWo54PSfUxvw9/TE6IoGhFP0Yg4xmYmMCErkfFZieSlxFofMkLrA7u78vpW/rr2MKsPHGd9Sc3HWj4nZicyf1QaC8akc/a4dJ/nxmrvdLHjaP3JgnBzaS2l1S20n2ISWKdDTruEXWJ0BEmx/yzakmIiiYly4hTBIe7+sgK0d7lo6/Q8OtxFfX1rB/UtHdS3dvZrqbyoCAdRTgdOhxDhECKcgsPT+nLi7U9xF3gdXUpnl4sOT8HX8/W/d2UxN59R9LFzhMoH9gmqyn3v7ON/Xt4FQHSEg+vmFPD5c0ZTOGJwW0FOZeOhGlbtrz55i/ZoXStVjW2MSo9nRkGK+1GY4tdRlQOlqoPeEh5q70+Hq5v5xtNbeN+zctGymXn86KqpAR/pXXK8iVX7j3Ok1t3t5WhdK7XNHVw0OYsb5xcyoh8jmIeDYCkAI4DdwPlAGbAWuFFVt3Xb51LgLtyjgOcDv1LVef05tjeB/sD2VnuniyO17ivOw9Xuf483ulvyBHGPtBVIiokkLT6KNM8giozEaApHxJEUY6NFTyXUPrB7097p4kPPLbGC1FjmjUob0CSg3nK5lKrGNkprWyiraaG6qZ0OT6HW4bklFhvpJDHmn300k2IiSIuPIiUuipS4SL/0sTnR+n2ihbGhteNkH5rmdnd/0+b2LpraO+nsOlHguU4Wdic+w8UzZt3pFCIdQoTTQYTTUyw6HEQ4BKfn+7PGpvc6jUWofWCf8Oq2Y+w61sAN8wv7NZ2HGRqhmE8ul/LU+lI6XcoN8wqsu0gQ6SufhnQUsKp2ishdwCu4p3L5g6puE5E7PM/fDyzHXfztxT0NzGdOdexQxj+UoiIcjEyP92nqBeMfqvog8CC432ADHM5JURGOkyMXh5LDIWQmxZCZFOPTSDlficjJwT85wTOQPqRcNCV7QKs5GNMXh0P6HFhkgtOQT4CkqstxF3ndt93f7WsF7uzvscYYY4wxZmBsngdjjDHGmDBjBaAxxhhjTJgZ0kEggSAilUD3XvvpQFWAwukuGOIIhhggMHEUqerHZ2A9DcunoI8BLJ/8IRjiCIYYwPLJH4IhjmCIAYIon4Z9AdiTiKzzZnTVcIwjGGIIpji8ESyxB0McwRBDMMXhjWCJPRjiCIYYgikObwRL7MEQRzDEEExxgN0CNsYYY4wJO1YAGmOMMcaEmXAsAB8MdAAewRBHMMQAwROHN4Il9mCIIxhigOCJwxvBEnswxBEMMUDwxOGNYIk9GOIIhhggeOIIvz6AxhhjjDHhLhxbAI0xxhhjwpoVgMYYY4wxYSZsCkARWSIiu0Rkr4jcM8jnKhCRt0Rkh4hsE5F/8Wz/roiUicgmz+OSbsd80xPbLhG52I+xHBSRDz3nW+fZliYir4nIHs+/qd3292scIjKh28+7SUTqReTLgfhd+JPlk+WTP1k+WT75k+WT5VO/qOqwfwBOYB8wGogCNgOTB/F8OcAsz9eJwG5gMvBd4Gu97D/ZE1M0MMoTq9NPsRwE0nts+x/gHs/X9wA/Gew4uv0/HAOKAvG7sHyyfLJ8snyyfLJ8snxyP8KlBXAesFdV96tqO/AEcMVgnUxVj6rqBs/XDcAOIO8Uh1wBPKGqbap6ANjriXmwXAH80fP1H4ErhyiO84F9qlpyin2G+nfhDcunj5/P8sl7lk8fP5/lk/csnz5+PsunXoRLAZgHHO72fSmnTlC/EZGRwExgtWfTXSKyRUT+0K0pejDjU+BVEVkvIrd7tmWp6lFw//ECmUMQB8D1wOPdvh/q34W/WD5ZPvmT5ZPlkz9ZPlk+9Uu4FIDSy7ZBn/9GRBKAp4Evq2o9cB8wBpgBHAV+PgTxnaWqs4ClwJ0isvBUIQ9WHCISBVwO/M2zKRC/C3+xfLJ88ifLJ8snf7J8snzql3ApAEuBgm7f5wNHBvOEIhKJ+4/hMVV9BkBVy1W1S1VdwO/4Z1PvoMWnqkc8/1YAz3rOWS4iOZ44c4CKwY4D9x/kBlUt98Qz5L8LP7J8snzyJ8snyyd/snyyfOqXcCkA1wLjRGSUpzK/HnhhsE4mIgI8BOxQ1f/ttj2n227LgK2er18ArheRaBEZBYwD1vghjngRSTzxNXCR55wvAJ/27PZp4PnBjMPjBro1hw/178LPLJ8sn/zJ8snyyZ8snyyf+sffo0qC9QFcgnt00j7gPwb5XGfjbsbdAmzyPC4BHgU+9Gx/Acjpdsx/eGLbBSz1UxyjcY8w2gxsO/FzAyOAN4A9nn/TBjmOOOA4kNxt25D+LiyfLJ8snyyfLJ8snyyf/vmwpeCMMcYYY8JMuNwCNsYYY4wxHlYAGmOMMcaEGSsAjTHGGGPCjBWAxhhjjDFhxgpAY4wxxpgwYwWgMcYYY0yYsQLQGGOMMSbMWAEYwkQkRUS+dJp9HhCRs4YqJhO6LJ+Mv1guGX+yfBocVgCGthTglH8UwHxg1eCHYoaBFCyfjH+kYLlk/CcFyye/swIwtP0YGCMim0Tkpz2fFJFJwG5V7eqx/W0RmeD5eoSIbPV8HS8i/xCRzSKyVUQ+ORQ/hAkalk/GXyyXjD9ZPg2CiEAHYHxyD1CsqjP6eH4p8HIv28fiXhcRYBruNQoBlgBHVPVSABFJ9l+oJgRYPhl/sVwy/mT5NAisBXB4u5gefxQiUgSUqarLs2ka7gWqwf3HcYGI/EREzlHVuqEL1YQAyyfjL5ZLxp8sn7xgBeAwJSJxQIqqHunx1Az++UcAMPvE96q62/P9h8CPROQ7QxCqCQGWT8ZfLJeMP1k+ec8KwNDWACT28dxi4K1etk8HYgBEZBxwBZ5mcRHJBZpV9c/Az4BZ/g7YBDXLJ+MvlkvGnyyfBoEVgCFMVY8DKz2dWHt2jO2rT8QMwCEim4HvADuAT3uemwqsEZFNwH8A3x+MuE1wsnwy/mK5ZPzJ8mlwiKoGOgYzCERkAzBfVTt6bN8LzFTVhsBEZkKR5ZPxF8sl40+WT96zUcDDlKp+rElbRBIBl/1BmIGyfDL+Yrlk/MnyyXvWAmiMMcYYE2asD6AxxhhjTJixAtAYY4wxJsxYAWiMMcYYE2asADTGGGOMCTNWABpjjDHGhBkrAI0xxhhjwowVgMYYY4wxYeb/AwsRmjtvoareAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig,ax=plt.subplots(2,4,figsize=[9,4.5])\n",
"ax=ax.flatten()\n",
"for a,l,r in zip(ax,legend,rho_list):\n",
" r.plot(ax=a)\n",
" a.set_title(l)\n",
" if not(a.is_first_col()):\n",
" a.set_ylabel('')\n",
" a.set_yticklabels([])\n",
" if not(a.is_last_row()):\n",
" a.set_xlabel('')\n",
" a.set_xticklabels([])\n",
"fig.tight_layout()"
]
},
{
"cell_type": "markdown",
"id": "712b2d16",
"metadata": {},
"source": [
"We observe above that when the correlation time is shorter than about 100 $\\mu$s, we see the fully averaged coupling, and when the correlation time is longer than about 1 ms, we obtain the rigid-limit value. In between, oscillations are damped due to dynamics being on the timescale of the coupling."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}