{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Shallow donor in Si\n",
    "Example of more complicated simulations, in which we compare the coherence predicted with point-dipole hyperfine couplings and one obtained using the hyperfines from model wavefunction of the shallow donor in Si (P:Si)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "import ase\n",
    "\n",
    "import pycce as pc\n",
    "\n",
    "seed = 8800\n",
    "np.set_printoptions(suppress=True, precision=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, as always, generate spin bath with `BathCell` instance. To get parameters we use ase interface. It allows to conveniently read structure files of any type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate unitcell from ase\n",
    "from ase import io\n",
    "s = io.read('si.cif')\n",
    "s = pc.bath.BathCell.from_ase(s)\n",
    "# Add types of isotopes\n",
    "s.add_isotopes(('29Si', 0.047))\n",
    "# set z direction of the defect\n",
    "s.zdir = [1, 1, 1]\n",
    "# Generate supercell\n",
    "atoms = s.gen_supercell(200, remove=[('Si', [0., 0., 0.])], seed=seed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculations with point dipole hyperfine couplings\n",
    "Here we compute Hahn-echo decay with point dipole hyperfine couplings. All of the parameters are converged, however it never hurts to check!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters of CCE calculations engine\n",
    "\n",
    "# Order of CCE aproximation\n",
    "CCE_order = 2\n",
    "# Bath cutoff radius\n",
    "r_bath = 80  # in A\n",
    "# Cluster cutoff radius\n",
    "r_dipole = 10  # in A\n",
    "\n",
    "# position of central spin\n",
    "position = [0, 0, 0]\n",
    "# Qubit levels (in Sz basis)\n",
    "alpha = [0, 1]; beta = [1, 0]\n",
    "# Mag. Field (Bx By Bz)\n",
    "B = np.array([0, 0, 1000])  # in G\n",
    "# Number of pulses in CPMG seq (0 = FID, 1 = HE etc)\n",
    "pulses = 1\n",
    "\n",
    "# Setting the runner engine\n",
    "calc = pc.Simulator(spin=0.5, position=position, alpha=alpha, beta=beta,\n",
    "                    bath=atoms, r_bath=r_bath, magnetic_field=B, pulses=pulses,\n",
    "                    r_dipole=r_dipole, order=CCE_order)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Time points\n",
    "time_space = np.linspace(0, 2, 201)  # in ms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For comparison, we compute both with generalized CCE and usual CCE coherence. Note a relatively large bath (`r_bath = 80`), so the calculations will take some time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "l_cce = calc.compute(time_space, method='CCE')\n",
    "l_gen = calc.compute(time_space, method='gCCE')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Hyperfine couplings of the shallow donor\n",
    "\n",
    "We compute the hyperfine couplings of the shallow donnor, following the formulae by Rogerio de Sousa and S. Das Sarma ([Phys Rev B 68, 115322 (2003)](https://doi.org/10.1103/PhysRevB.68.115322))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PHYSICAL REVIEW B 68, 115322 (2003)\n",
    "n = 0.81\n",
    "a = 25.09\n",
    "\n",
    "def factor(x, y, z, n=0.81, a=25.09, b=14.43):\n",
    "    top = np.exp(-np.sqrt(x**2/(n*b)**2 + (y**2 + z**2)/(n*a)**2))\n",
    "    bottom = np.sqrt(np.pi * (n * a)**2 * (n * b) )\n",
    "\n",
    "    return top / bottom\n",
    "\n",
    "def contact_si(r, gamma_n, gamma_e=pc.ELECTRON_GYRO, a_lattice=5.43, nu=186, n=0.81, a=25.09, b=14.43):\n",
    "    k0 = 0.85 * 2 * np.pi / a_lattice\n",
    "    pre = 8 / 9 * gamma_n * gamma_e * pc.HBAR_MU0_O4PI * nu\n",
    "    xpart = factor(r[0], r[1], r[2], n=n, a=a, b=b) * np.cos(k0 * r[0])\n",
    "    ypart = factor(r[1], r[2], r[0], n=n, a=a, b=b) * np.cos(k0 * r[1])\n",
    "    zpart = factor(r[2], r[0], r[1], n=n, a=a, b=b) * np.cos(k0 * r[2])\n",
    "    return pre * (xpart + ypart + zpart) ** 2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We make a copy of the `BathArray` object, and set up their hyperfines according to the reference above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "newatoms = atoms.copy()\n",
    "\n",
    "# Generate hyperfine from point dipole\n",
    "newatoms.from_point_dipole(position)\n",
    "\n",
    "# Following PRB paper\n",
    "newatoms['A'][newatoms.dist() < n*a] = 0\n",
    "newatoms['A'] += np.eye(3)[np.newaxis,:,:] * contact_si(newatoms['xyz'].T, newatoms.types['29Si'].gyro)[:,np.newaxis, np.newaxis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we set up a `Simulator` object. Because hyperfines in `newatoms` are nonzero, they are **not** approximated as the ones of point dipole."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "calc = pc.Simulator(spin=0.5, position=position, alpha=alpha, beta=beta,\n",
    "                    bath=newatoms, r_bath=r_bath, magnetic_field=B, pulses=pulses,\n",
    "                    r_dipole=r_dipole, order=CCE_order)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "shallow_l_cce = calc.compute(time_space, method='CCE')\n",
    "shallow_l_gen = calc.compute(time_space, method='gCCE')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare the results\n",
    "\n",
    "We find that the point dipole gives a poor agreement with the experimental data. Model wavefunction, on the countrary, produces great agreement with the experimental coherence time from work of Eisuke Abe et al. ([Phys Rev B 82, 121201(R) (2010)](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.82.121201))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABRVElEQVR4nO3dd3zN1//A8de5N3tPIRJJxCZmbLHVKqql9uigqLbor6pT+XZoq1R1oPYebZW2VKm9JWoTM4idiOx9z++Pz0VCEOTmZpzn43Ef7uecz/183veT677v53M+5xwhpURRFEUpvnTmDkBRFEUxL5UIFEVRijmVCBRFUYo5lQgURVGKOZUIFEVRijkLcwfwuDw8PKS/v7+5w1AURSlUwsLCoqSUnjnVFbpE4O/vT2hoqLnDUBRFKVSEEOcfVKcuDSmKohRzKhEoiqIUcyoRKIqiFHOFro1AUZTCIz09ncjISFJSUswdSrFhY2ODj48PlpaWuX6NSgSKophMZGQkjo6O+Pv7I4QwdzhFnpSS6OhoIiMjCQgIyPXrTHZpSAgxWwhxXQhx5AH1QgjxnRDitBDikBCitqliURTFPFJSUnB3d1dJIJ8IIXB3d3/sMzBTthHMBdo9pL49UN74GAz8ZMJYFEUxE5UE8teTHG+TJQIp5Vbg5kNW6QLMl5rdgIsQopSp4jm+5x/+HNmUV154ma1zV5FwPtJUu1IURSlUzHnXUGngYpblSGPZfYQQg4UQoUKI0Bs3bjzRztb9vpRPlu1gXdgm+p+woPbUUFo27c4f0+c90fYURSkc9Ho9NWvWvPOYMGGCSfe3evVqk+9j8+bN7Ny5M8+2Vygai6WUM4AZAMHBwU80k06DevVZs3w+U8cNI8PJkXU7jvNl2AbCbkYzv35rTvzwMp7evrw8bmaexq4oinnZ2tpy4MCBfNlXRkYGnTt3pnPnzibdz+bNm3FwcKBRo0Z5sj1znhFcAnyzLPsYy0yiQbsX+aefPQ39LQjp3JRPv3yN84tHM7xfZ95aeoBLR3fyyXfzyEhVt7kpSlEXGxtLxYoVCQ8PB6BXr178/PPPADg4ODBy5EiqVq1Kq1atuH0V4syZM7Rr1446deoQEhLCiRMnABg4cCBDhgyhfv36jB49mrlz5zJ8+PA7dUOHDqVBgwaULVuWzZs38/LLL1O5cmUGDhx4J55//vmHhg0bUrt2bbp3705CQgKgDakzduxYateuTVBQECdOnCAiIoJp06YxefJkatasybZt2576eJgzEawG+hvvHmoAxEopr5hsb7au4OQD147eKRr50adUuLCcxp7J9K0qmd3JCt3pDdlft3saLO4JW7+GdJUkFOWpNG9+/+PHH7W6pKSc6+fO1eqjou6vy4Xk5ORsl4aWLVuGs7Mz33//PQMHDmTp0qXExMQwaNAgABITEwkODubo0aM0a9aMcePGATB48GCmTp1KWFgYEydOZNiwYXf2ERkZyc6dO5k0adJ9+4+JiWHXrl1MnjyZzp07M3LkSI4ePcrhw4c5cOAAUVFRfPrpp2zYsIH9+/cTHBycbTseHh7s37+foUOHMnHiRPz9/RkyZAgjR47kwIEDhISE5Oo4PIzJLg0JIZYAzQEPIUQkMBawBJBSTgPWAB2A00AS8JKpYrltzKZMdh5fztYXZmJIT+NsdAq3rl/mtTYJLFqezhv1rWDxh/DZs3deM+z/PqRp6Uwqu/7JqU2H6DZ2ganDVBQlDz3o0lCbNm1YsWIFr7/+OgcPHrxTrtPp6NGjBwB9+/bl+eefJyEhgZ07d9K9e/c766Wmpt553r17d/R6fY7779SpE0IIgoKC8PLyIigoCICqVasSERFBZGQkx44do3HjxgCkpaXRsGHDO69//vnnAahTpw6//fbbEx6FhzNZIpBS9npEvQReN9X+cxJYvgLJNyKQ6SnoUuMY3ciaGf9dplKbs8w+kEZAYFkSj5xg6NwZuA0cTGZKIltOxlCyfAuWhYVz4MwSnhs1CQvHHEdyVRTlUTZvfnCdnd3D6z08Hl7/mAwGA8ePH8fOzo6YmBh8fHxyXE8IgcFgwMXF5YFtDfb29g/cj7W1NaAlmNvPby9nZGSg1+tp06YNS5Yseejr9Xo9GRkZuXlrj61YjTU06OUBTGlnjYg6CUnRpGZKktIM1He5RdwYR0441OPDTans+uJDiIlBH3ueo8Mc+HjUEKZOn8P+QfboF71A7wbeRBzcfv8OLv+nLh8pSiExefJkKleuzOLFi3nppZdIT08HtATxyy+/ALB48WKaNGmCk5MTAQEBrFixAtB68GY9i3gaDRo0YMeOHZw+fRrQLk2dPHnyoa9xdHQkPj4+T/YPxSwR4KWdkmVcOsT+vTtZeCidGc/a4HJ1O8LOjZ4dWnL9/xwI9fOGuDiI0hqS8KiAT+3WuDZ7jVMRl1iy5wq71/2q1aVqjTpc3AczmiNXvITBRFlbUZTHd28bwZgxYwgPD2fmzJl88803hISE0LRpUz799FNA+3W/d+9eqlWrxsaNG/n4448BWLRoEbNmzaJGjRpUrVqVVatW5Ul8np6ezJ07l169elG9enUaNmx4pyH6QTp16sTKlSvzrLFYaFdoCo/g4GD5pBPTyIx0AkvY8EKrejzXpQujxnzAvOdsWXYkA+cSPoyYth6+q8nHhtd4/6MvWPh+DzatW82CvTfR2TgAkBFzidgvKuHSdSL6al3YMbICF8v24tShMMqmH2f9mXTKVG/C+Pkb7+74+nFwKwsW1g+ITFGKpuPHj1O5cmVzh/FYHBwc7ty1U1jldNyFEGFSyuCc1i8U/QjyirCwpHfdEtQqY03jip7sedWBLkuTWB2ewYBmkGntwve7U7lS4hQr/j1C9OFQIpMs7yQBAAtHD9ztdJAeB7cu8NmWZM6tno2VTtK0bnU8y6QSd3oPGAyg00F6CnJaU/ZYNaLBmLz5BaEoipKXilUiAPi0d11IT4IkbfQLBwcnJrZJ5O13+iBtnfm/9am0rXuId0eP4FLDKN7tWD/7Biys+SFMEmgZRgvXivzQwQZPe4G9JSS/Mhe7G/vhj7cg9gK4+kPideaGJfDy6tXsqvATDZ4fmv9vWlGUXCvsZwNPotglApxKk3ZqM1MXrmX7zjRWftwZjq8GV3+ETsf1sYEcSvemzQ8nMHhaoD929b5NfLEtmbbJB7mW6crA7xI4vfD/CCzphF2Z6kAyAPL6CYSrP8f27+GX4xmMaGBN3UtzQA4BNQiXoigFSPFqLAY+//M09u+dIDM5Dr2lNWHRtjy/LIlTsVpOdPUoQcOK3gT3+T8azoln878n4Uj2kbSPf1ybn1+uRe2yrnzRypqSbUdAq48ASLL3xWdSPBOnaJ1kbl27wLkYA/2fb4s+OhySY/L1/SqKojxKsUsETerXZmwza96oI/jl9Tocj7Vl5YkMhJs2icOKo6nM/ucwbb0ScLYRiHRbmDgx2zYcXd3Rpd4iqJQtY0JssXe7O2iqnbs3zwU5U8VVu3OoUQUPjr3ugFelhkzelcqFo7vz780qiqLkQrFLBE2bt+LDptbYxp0GO3d6v/MN5w5so1ytJgCs2B/NDxvO8XwFHev72bO5RmuwsoIsd1etPBzPt38d5+b1yyTonLVG4Sy+f60pHf20S0QkXAfgqt6bUf+kErZrm9aQnBKbP29YURTlEYpdIsDJm5QMScjseCb8fQGdhQX+NZrcqV4wujOH3vKiqnMqKZmC0MqtYcaMbNf1/zgYxZRNlxjy3VqCv7+/DQHPSiRdOoHMzGTywjUMWJ1BUOO2RI925Lk63lxYM4mF/cuSEvew6RoURckrn332GVWrVqV69erUrFmTPXv24O/vT1RUVK63sXnzZp59Vht+JuvAckVBsWssjsURl8+0Hnk9nZzvq7d29oKUWwyduIwb1+Fkp7IkpmZgf+oEVK0KQjB9VBcs/pvLuhgvoqs53reNhfti6D/+KhHd9hIfG0NMuiWWTiVwc3aEuEi2/xtKvxU3CR6xh0qN2pv8PStKcbZr1y7+/PNP9u/fj7W1NVFRUaSlpZk7rAKl2J0ROJXw5fM2jux51Z7Xe90/k2bYhQRG/5OEs+EWPiXdSU5MYNnU+RAUBFu2AGDp6I7ISKZdWejToup926jTuAWfNLfGKu48H7f3ZvWoRiAEi0/ZMWvlZqo53OLYMHsCXYvd4VeUfHflyhU8PDzujNnj4eGBt7c3AFOnTs02xDPA3r17adiwIbVq1aJRo0Z3hqp+kIiICFq2bEn16tVp1aoVFy5cIDMzk4CAAKSU3Lp1C71ez9atWwFo2rQpp06dMuE7fnzF7oxA6HS816kCRJ8CO/f76sOvJvHdnjSODEuhbOveTBvwBp+Xr8bLjo6wcCE0b07YuViW/JNCv1qnCGwchMM926jcqB0fN7MGeVVrI/DSksXiA4nciLnC8EuJDK9rxddJOVxWUpQiatwfRzl2OS5Pt1nF24mxne7/MZbVM888w/jx46lQoQKtW7emR48eNGvWDLg7xPOPP/7IxIkTmTlzJpUqVWLbtm1YWFiwYcMG3n//fX799dcHbv+NN95gwIABDBgwgNmzZ/Pmm2/y+++/U7FiRY4dO8a5c+eoXbs227Zto379+ly8eJHy5cvn6XF4WsXyJ+m+aDucJ8Sx/fj9X8Q9n3+W5A8cKeemQ+fmT/tX3kFW60j6C91gxQpITubk1QR+3JdGzR9j+OS3I/fvwNaVdAdfLh7ZSaspR/h5u7afFR91Y2tfPT91tCFTwhbjLwSZfAsWvQhhc034rhWleHJwcCAsLIwZM2bg6elJjx49mGuc4yDrEM8RERGANmlN9+7dqVat2p25Ax5m165d9O7dG4B+/fqxfbs2IGVISAhbt25l69atvPfee2zfvp19+/ZRt25d07zRp1DszggA3lgSTlwqOHt631enc/DkyA0DbRcmMa/0LYYMHMqQhWGE1alGg7lz4M8/6fl8R3omzWb+wXQqt6mX4z56rYjjwJlVBDhnImy0dgRbr0CwgIE1rfD7Np7o1dup2fsiLzSvwWeNU6kvM6HOQFO+dUUxm0f9cjclvV5P8+bNad68OUFBQcybp81VntMQzx999BEtWrRg5cqVRERE0DyXE+Dcq2nTpvz0009cvnyZ8ePH8/XXX7N58+Y8mUgmrxXLM4Ivh3Ri/nM2BNVucF9dfKY17/2bio0FlC4XRJPyHhiizzPr9CXw9oYlSxC2bgghGFDTinp1aua4j9debMunzfWs72fPqz06AXDqlp4uS5PYdj6TbaPr8FPPshz9fRL69DimHLDm1e83mfJtK0qxFB4enu2a/IEDB/Dz83vg+rGxsZQuXRrgzpnDwzRq1IilS5cC2gilt7/o69Wrx86dO9HpdNjY2FCzZk2mT59O06ZNn+LdmEaxTATNOvWkX/MK4HT/GYHe0ZNN5zJ4u6E1les1x8HagrTts1n63Xj480+YO5eYVMHgP5KZ818aKbp7Wwg0bTp1p2c1S23BoQQANw0OrA7PoOncRMpUqI5dymUauUSxbngVatSpz76IWDJTit84J4piSgkJCQwYMIAqVapQvXp1jh07xieffPLA9UePHs17771HrVq1cjURzNSpU5kzZw7Vq1dnwYIFTJkyBdDONnx9fWnQQPvBGRISQnx8/J0ZygoUKWWhetSpU0eammF8CSk/LSWlwSCllPKTuWtk6WHzZGRMkpRSyqgLpyQgAfnvgkk5byTuqpzZyUa62Qp5cutvUkopM2OvyH2D7OWad5vL/T8Nlp+3spNpX1aUcsVLUh5aIeVYJymvHNZenxgt5fpPpExPNfn7VRRTOXbsmLlDKJZyOu5AqHzA92qxPCN4lG4rkmm9IPFOJ7Le7Zpg4ejO9lM3YOZM3KfP5srbDizsaktQzdo5b8TRi1f/SOFmssTB0xcAnaMXwZX9af9cD3afS+T9f5OwejecFeE68KigvS5aO4W9tns5R3/7mowzW03+fhVFKd5UIshBut4WJ2fXO8vlSjhgcek/fp45Cw4cgEnfUtLdhT7VLfH0LffA7fz1f004NMSeUoHVtAIh4K2DUP81+vfqTvhwe96sZ0Wlei0xuAbQfUUSP89dBMDCleuo9lMiB7evMeVbVRRFUYkgJ/P+2sO8v3bdWRZCoDuznR0r5yCffx6Skxn8WyxvrElG2ro+cDsderxMUHBjsLS5W6jTgxDYlypHBXc9Uzo6EtT8eXQ2jkSlWpFwIxKALnW0xqovZ/5mmjepKIpiVCxvH30UV2//+8reGfcln64/z5nKtSnn4sLPobcAmGpl9+ANBb+kPXKQ6VCK0etSaN/An9bGRLHp4zaQom23nHMma3rbUsrbGm6ehT9GQNfp4FQqx+0piqI8KXVGkEtta5dH6PTsirgF7dtz5HlHVg4s/cTb09u7MWl3GpN3Jt0t9KgAUadASg4cO0U5Nx017a+y5NPB7N72L+yf//RvRFEU5R7qjCCXfN1s0R1fx5STq+jXqytV14dStUHFp9rm1TNHcfK4+wt/3ek0xky9ytruh3hlxh5cLdJ5t5Gk98J/KeMsaH/sK37c+jY6C/VnUxQl76hvlFwSQmAff5HwCxcx/DIdXffuT71Nr7JVsi3befhR0kGQdPUUPz7nQYy1L898vY2Pm1ph51ebSb/t4d1dvxEQ8uJT71tRigu9Xk9QUBAZGRlUrlyZefPmYWdnd6c8PT0dCwsL+vfvz8iRI9Hpit+FkuL3jp/C2K+m4N5tHMevGgfOisvbAbRCnnmWtX3sKeuUQX2PRNq1ac6moX6MaFuOkT/8ReS7pQiI3pKn+1SUos7W1pYDBw5w5MgRrKysmDZtWrbyo0ePsn79etauXcu4cePMHK15qETwGBqX9wRg15lomDUL3Nzg+vW824Gz1t8g5cJ//HMykWup1jQfORPXPj9j5eiGpV9duHE87/anKMVMSEgIp0+fvq+8RIkSzJgxg++//x6ZZTbC4kJdGnoMpZxtMexZyFc7E3l14oeQmQkbN0LPnnmzA0sbnl2eifxtNmsOJ7Gg8mX6Ptv6TvV3m6+y7/BJFgzKm90pSr5aOwauHs7bbZYMgvYTcrVqRkYGa9eupV27++chAShbtiyZmZlcv34dLy+vvIyywFNnBI/J18OJqKRMMmrWAhcXWL8+T7dfwcedxl4pbH/JjjYtW2SrS8i0JDouOU/3pyhFXXJyMjVr1iQ4OJgyZcrwyiuvmDukAkedETym9z4ay/DF/3H4agK1WrbUEoGU2eY0fhqTXnsGDi8HdOBXIVvd+wPbQ8BhMGRqHdMUpTDJ5S/3vHa7LeBRzp49i16vp0SJEqYPqoBRZwSPqV6AGwB7z0ZDmzZw8SLk5bRzrn6ER2Wy9lQ6Gdb39Fq2dQMkJN/Ku/0pisKNGzcYMmQIw4cPR+TRj7rCxKSJQAjRTggRLoQ4LYQYk0N9GSHEJiHEf0KIQ0KIDqaMJy+UcLQhfvVnfD1mGHTsCN99pzUa55HVB6Op9EMiHRYnI22zb3fr8as0m5tIRPihu4XxarpLRXkSty8ZVa1aldatW/PMM88wduxYc4dlFia7NCSE0AM/AG2ASGCfEGK1lPJYltU+BJZLKX8SQlQB1gD+poopr1SuWZcLselIHx/EG2/k6bbLlK9Kl4oW9KnjgqWtfbY6nY0LAKmxNwBIizyExc8h6HouhMqd8jQORSkqEhJynuMjMzMznyMpuEx5RlAPOC2lPCulTAOWAl3uWUcCTsbnzsBlE8aTZwYNH4kuqCOXbiVDVBQsWgS5mMAiN2o2bs3vPe3o3tD/vromTUPYMtCeit7a1JcLZkzBZUIckeum5sm+FUUpnkyZCEoDF7MsRxrLsvoE6CuEiEQ7G8jx57UQYrAQIlQIEXrjxg1TxPpY6vi5Ig2ZbDt6Qbt9tG9fCA3Nm407+fDXKQM7Ludwsnb7UlFyDACVnZJoUsaCYT9t4ux+1dFMUZQnY+7G4l7AXCmlD9ABWCCEuC8mKeUMKWWwlDLY09Mz34O8V4US9lz6cSDfTvgftGyp3TG0YUPebFxvwbOLE2g3Oey+qiRpTePZicz7ZS0AjVxvMr1vBcKuZHLqn1l5s39FyWPFsYOWOT3J8TZlIrgE+GZZ9jGWZfUKsBxASrkLsAE8TBhTnrCytKB6x/4Iv9rg4QG1auVpf4KTu//m1KG995XbupbEzlJgJVMwZGRw5dRBfOt24OKP3Wir2wGpar5jpWCxsbEhOjpaJYN8IqUkOjoaGxubR6+chSn7EewDygshAtASQE+g9z3rXABaAXOFEJXREoH5r/3kQs+XhjBj61mS0zKxbdlSu3soJQUe8w+Qk/L12+ZYLnQ61r9WBiqX5mTYZip+eY359nH0e2k08ueWnFv5KWV7mudebUXJiY+PD5GRkRSES7rFhY2NDT4+Po/1GpMlAillhhBiOLAO0AOzpZRHhRDj0SZRXg28DfwshBiJ1nA8UBaSnw51/FxJuXWdzQdO0T4kBCZOhP37oVEj0+7Yzg2Sb+KSepFv21rTpE1n8Anm/QM+/PTVV5xr+hqu3gGmjUFRcsnS0pKAAPV5LOhM2rNYSrkGrRE4a9nHWZ4fAxqbMgZT8XeQXPrpJaZn/B/tvx0HFy6Ar++jX/iUXltxleSMK8wfX463GtlDnVYA9B31OWUNvbA7vAC8P37EVhRFUe4yd2NxoRXoW5JKL76LZWADsLPLlyQAUMrDmdIOsH/nZpKcyt2ZD7lqs+cY1Lsr1v/NgtT4fIlFUZSiQSWCp9DhhV6cTXPSGsK2bNFuI82j/gQP8km/ZoxvYUPI/zbz7r9p2epkk5EsC73Bsi/ztpOboihFm0oET6Galw2RR/eyP/w8XLmidSzLxeBWT8XOFZF0nd9etGXwywOyVQnfukw9aMNXP68wbQyKohQpKhE8BbfMm1xf9iHzVqyCkBCtcNs2k+5zzqbTlJ+aQFM/PUHP9L2vfvO3gwkb4mzSGBRFKVpUIngK7ZoE49v7M+zL14fSpaFsWZMnglKlvIlOkhzN8Af7+7tcWDh5QnoipKt5CxRFyR2VCJ6CtZUl9Zo051SMQSsICdESgQnvgG3ZtBGZEhacsMqx/sDFREb8ncKN8+GQkaYajhVFeSSVCJ6Sj2UC2/9YTEpqGjRvDj4+2kB0JmLl7sfRYQ6MeivnBuHzMenMOZDG1Yhw2PoVfF0euX+B6tmpKMoDqRnKntbVE1xb+wPrd/Wj08CBMHCgaffnE4z/qPXgWz/H6s6dniX2xnfg58b/jVxIRlQM1fYOxrr+bvp98INpY1MUpVBSZwRPaUDPFyg9ZDbJ9qXuFpry17cQUKbBA6fGFLfbDRKjSE+KJ93Bm4l7JNs2rMlxfUVRFHVG8JRqlPXG2bMURy7H8SLAZ5/BkiVw5IhZ4kmzdOLtNcm0c9jClA72UP4ZbmXa4nJkDqTEgo26o0hRlOzUGcFT0ukEnnEnWTl/hlZgbw9Hj8KlewdazR+WDu4sP5bB8fBTkHAdnLxxCe4GhnQ4lXcjpCqKUnSoRJAH0s/v59iaOaSkpd8ddG7XLrPEInQ6ro2rQK9qltSZEc+6I9HgE8zYHXreeleNQaQoyv1UIsgDw98eg88bizgblQw1a2pDUZspEQBg50HalWN42euwc/cGnZ6r+tJ4p52F6U3h0n7zxaYoSoGjEkEeaFDJF6HTc/jSLbCygjp1zJoIvt8Zx8wd11jTx46QFm0A+HrB37wy8mO4EY7hwFKzxaYoSsGjEkEe8HOzI3X/78z5ebpWMHAgdOhgtniOXk8n7EqmtuDkrf3jWRqPju/TcxW0HT3HbLEpilLwqESQB3Q6AZGHOBy2Wyt49VX48EOzxfPTiC6ElLGg0ewksHPPVtc8uArtAgymvcVVUZRCRSWCPDL8q1nYtXub1AzjL/GEBIiMNE8w9u6UdBBULOVwX3+DIQN68HZwhnZHkaIoCioR5JmaZdxJz5ScvGqcQL5mTRg50iyxbAu/yV+nMpjQo8b9lV5VyDBIks6HacsGA6TE5W+AiqIUKCoR5JFyrhbcWP01M+Yt1Arq1TNbg3GitOZktIE4C7f76mIsS2P/eTzTZswEIHTROP59PRCZfCufo1QUpaBQiSCPlC/tjuHmBcIjLmsFDRtqncouXsz3WNq2bUtiumT14dj76lx9AhnTwo0GWhsyU2YuZsCKaETkvnyOUlGUgkINMZFHdDod3T5bws1E4/SRDRtq/+7alW/zGd+WbunMM2UtCAwsm2P9uH4hkBwNwOSOLoQH2nJ295+ULd8mP8NUFKWAUGcEeai6jzPhV+NJSc+EGjXA1tYsl4esSlZk+v9G8Nzg93Newasq4ceOkJ4Qg0fKWd5Zn8qA/83P3yAVRSkw1BlBHrKLj+TinBH82mg6fZ5tAfPmQeXK+R+ITg8dvnpg9drz1nT47iZzHN4g9WASQ5uWpKRFPKSngKVNPgaqKEpBoM4I8lD9yn7orG0JvxyjFXTvDtWqmTeoHIR0e43vOzoQdXAdQ/5KoV7H/rQJAC6FmTs0RVHMQCWCPFSnSjmqvjqRZNdArSAhAZYvh7NnzRvYPRzcvHi95zO8XSeViLfdKffsWxy5nsn+DcvNHZqiKGagEkEeEkIQVNqZgxdvaQVxcdCjB6xebda4cpLi14o/Tmag8yiH3smL51akM2HWKq0y/iqkJpg3QEVR8o1KBHks+cgG/v2wM9G34sHbW7tjaM8ec4d1n+suNemyNJkFR7Wex/NeD+Hzti5Ig4Hlw+sRuXSUmSNUFCW/qESQx+pWr4xdxcb8d/aqVlC/Puzebd6gclAmqBEzxvSm66D3AGjcpAnldJFcCQ+jx/yLrPx7s3kDVBQl36hEkMf6PNce97avcynFUiuoXx8iIuB6wRvbZ9AXi6jctDMANy1Ls+JwEuLk3/zV2xZD/HWkwWDmCBVFyQ8qEeQxLydrPBysCTud5YwAYF/B7rl7NtGWF39JZvfaxZy/JRnxVyyRx9VdRIpSHOQqEQghvIQQs4QQa43LVYQQr+Tide2EEOFCiNNCiDEPWOdFIcQxIcRRIcTixwu/4BFCcGvNJOa821srqFdPu2vIjPMT5EZQk/bsH2xP7NUI3O0Eh4faU0ofY+6wFEXJB7k9I5gLrAOMI9RwEhjxsBcIIfTAD0B7oArQSwhR5Z51ygPvAY2llFUftc3CokmrdlhVbU1yWiZYW0NAwH3DQRc01k7u1KoUwLgtqfwe6Ua1EnosYk5DZgYYMs0dnqIoJpTbROAhpVwOGACklBnAo74d6gGnpZRnpZRpwFKgyz3rDAJ+kFLGGLdb8C6kP4E+vXrgGNyFY1eMg75t2QKvvKIN+VyA7Y71YFQDa74Z8xprL9jy15q/4Ofm8E0lZox6TrUZKEoRldtEkCiEcAckgBCiAXD/0JbZlQayDr0ZaSzLqgJQQQixQwixWwjRLqcNCSEGCyFChRChN27cyGXI5hNU2hlDWjLbD57WCiIiYPZsOHHCrHE9yorDSby7IQWvqo35YnsaL0xYw9njB9l9w5bXJq9i92o1xaWiFEW5TQSjgNVAoBBiBzAfeCMP9m8BlAeaA72An4UQLveuJKWcIaUMllIGe3p65sFuTcvLyZqrM4cwY9KnWsHtBuMCeBtpVmPeGcWVjwPR+dRh0ZiulHfT0WVFBuVfns7ffeyo4lGwL28pivJkcjXonJRyvxCiGVAREEC4lDL9ES+7BGQdf9nHWJZVJLDHuK1zQoiTaImhYN9i8whCCBr2eoskC2etoEIFcHHROpa9/LJZY3sYzwY9oEEPAHyr1OXvPsu5VHs07uVq07acJaRfNXOEiqKYQm7vGnodcJBSHpVSHgEchBDDHvGyfUB5IUSAEMIK6Il2VpHV72hnAwghPNAuFRWsgXmeUJduPYh2CiQpLQN0Ou3uoQJ+RpBNUHdKd/uCer3eAwtr9sd7sH3HTnNHpSiKCeT20tAgKeWt2wvGxt1BD3uBsUF5ONrdRseB5VLKo0KI8UKIzsbV1gHRQohjwCbgHSll9GO+hwKpSkkHUq6dY/shYzvB7YlqMjLMF9TjsHeHRsNBr500vrM2lnfmbDNzUIqimIKQUj56JSEOA9WlcWXjraGHjLd85qvg4GAZGhqa37t9bAdOnKFW5XL0GTGWhZM/ASkL/C2kD3Pk+77YnN9Mua8jzR2KoihPQAgRJqUMzqkut2cEfwPLhBCthBCtgCXGMuUBalQsS9keH2JbwXgmUIiTAEC1WnUpZx8PapJ7RSlycpsI3kW7dDPU+PgXGG2qoIoCIQTN2nXmfEqWGb8GD9b6ExRCN3WeLD6czuXje80diqIoeSxXiUBKaZBS/iSl7GZ8TJdSqu6mj+DvkMnBrX8TdStOK0hKgrVrtctEhUxkig19fktmx6Z15g5FUZQ8ltu7hhoLIdYLIU4KIc4KIc4JIYrE3T2mJK6d4vrvE1i1YYdWUL8+XLkCkYXvOnului04OsyezjU8APhm+HPM+6C3maNSFCUv5Hby+lnASCCMRw8toRj16tKWH/dOItM9QCto0ED7d88ebcKaQsTKwZUqZX0gQessPnH+n/SqZgkJN8Ch4HfyUxTlwXLbRhArpVwrpbwupYy+/TBpZEVAOZ8S+FQI4sSNFK2gRg1tELrC1J8gi71xXny7aC0y9jKLu1rTu5oe9s81d1iKojyl3CaCTUKIr4UQDYUQtW8/TBpZEVE64yprVyzQFqys4NVXtZ7GhdDfF60ZszKCxCNraBFgwYU0Zxr1/5i05ERzh6YoylPI7aUh42A5ZL0HVQIt8zacoif59B5O/PojUbfG4uHiBN9/b+6QntgbI0bxdsndHFo5BV2sxL7RIDLW/4+M5S9j1Ws2WNmbO0RFUZ5AbscaamHqQIqq14YO45h7EyJuZeDhYixMTtaGpLYvXF+crlVagJ0tY5cd53qaNQdmjKVJOSfstv0PVgyEPivMHaKiKE/ApDOUKdC4WgB6GwcORRpH7b5yBZycYP588wb2JCysWXnNFzdbwdxRz4IQ2LccxfnKQ+n75SpO7Fxr7ggVRXkCJpuhTNF4OdmgO72VJQuMY/mXLAluboW2wfhwgguhlzOpXP/uVUGb4D78czaDY2tmmDEyRVGelClnKFOMMk7vJHT979qCENptpHv2mDWmJ/X++C85/W4ZrCu2vlPmVbYql2b253mHMEhLMmN0iqI8CVPOUKYYDRs/FZcXPycx1TjyaP36EB4OMYVvcngL3zrwbgS4B2Yrt6w3EFJi4exmc4SlKMpTMPcMZcVC3fKlkAiOXTEONXF7xrK9RWfcnjibMlT7MYGZ8xebOxRFUR7TI+8aMg453cz4eJwZyhSjSiXsuLlhOnPsI6j70evaJDUTJ0LFiuYOLc84lvSnoqcFHlap5g5FUZTH9MhEIKXMFEL0klJOBo7mQ0xFTml3R9IiwjhwuLRW4OgIb79t3qDymNDp+HVwZfB2MHcoiqI8ptx2KNshhPgeWAbc6UYqpdxvkqiKoD4TV3L+ZvLdgqgo2LkTOnUq9HMV3OHsA3H3TkutKEpBl9tEUNP47/gsZapn8WMI8nFhY/gNElMzsLe2gJUrtfkJTp2CcuXMHV6eeP+PS/y16zgHVQ8TRSlUcjsfQYscHioJPAbX9CiuLh/LirWbtILbDcaF9DbSnFStGEgrP8Cg3Vn81dBOdK3ri6GwzNOsKMWU6lmcT+qU8yYzPppDZ69oBVWrakNMFNKOZTnp07U9k9paQ/xVSI7hTNhGwiMuYzjym7lDUxTlIVTP4nxSrbw/tUb8jCxdXSvQ66Fu3SJ1RoCzDwAyNhJCZzO9gwXH3quExZbPISPNzMEpivIgqmdxPgoq7czhS1n64dWvDwcOQEqK2WLKS+HX03CZEMfCBfP57eeJZAa0QD47mTV7T7Lmpw/MHZ6iKA+gehbno8yzu9nyaW+uRBl7FA8bBocPa5PVFAEly1WnX3VL9m9axQsLrrLDshmUa82b/2Qye4HqaKYoBVVu7xq6t2exJ9DNZFEVUUGBpbHy9Cf0ZCSdPFyhTBlzh5SnnL3KMLWLBxmpcbQtF0BIjzcROh3/TBpCmVPzIekm2LmZO0xFUe6R27uG9qP1LG4EvAZUlVIeMmVgRVHfrh3w7Po+1w1ZOl0tXw7Tp5svqLwkBDiXxiChXd83EHo9AGVbv4yFzgDH/zBzgIqi5CS3l4YA6gE1gNpALyFEf9OEVHSVcLTBy8mawxezDDb3yy8wYYL5gspjLu8ewPrTeKjZ525hyerMOenMm++NM19giqI8UG5vH10ATASaAHWNj+CHvkjJUeLmWcx6u/vdgvr1ISICrl0zW0x5afaEd5k0ojs4lLhbKAQRluUIO3GBtMvHzBecoig5ym0bQTBQRUopTRlMcVCrTjCRCQbiU9JxtLHM3rGsc2fzBpcHnn/j0xzLP5w8n0++r4UImw7eU/I5KkVRHia3l4aOACVNGUhx8dKAvrg07c/xK/FaQe3aWp+CotSfIAeWbj6IWr1J2ruY+CtnzR2OoihZPDQRCCH+EEKsBjyAY0KIdUKI1bcf+RNi0VKttDNSGtgTflErsLODGjXgwgXzBpYPblXpT5lvopk7cYy5Q1EUJYtHXRqamC9RFCMlHG24Pmc4P+yoyhvt1miFO3cWmb4ED+MSWIfBDZyp6aH6IipKQfLQRCCl3HL7uRDCC62RGGCvlPL6ozYuhGgHTAH0wEwpZY63xwghXgB+AepKKUNzGXuhVbt9b+J19ncLikESuO3zvvUBNa+xohQkub1r6EVgL9AdeBHYI4R4aIcy48xmPwDtgSpot5xWyWE9R+AtoGhfJM/ixX4vE1+yFgm35zCOjYUuXWDZMvMGlg+kezmunNXuHEqOj8WQqc4OFMXccttY/AHar/UBUsr+aH0KPnrEa+oBp6WUZ6WUacBSoEsO6/0P+BIoGgPu5EJ1H2fS42+y88g5rcDREbZtgw0bzBtYPvj238t4f3aBqLOHaV7Vi5p+LqSeKzpzNytKYZTbRKC751JQdC5eWxq4mGU50lh2hxCiNuArpfzrYRsSQgwWQoQKIUJv3LiRy5ALLi+rVC790J9Zc+dpBTodNGwI27ebN7B80OaZtvzQwQbLI8v57hk9z/gbsF7yPKQlPvrFiqKYRG4Twd/GO4YGCiEGAn8Ba55mx0IIHTAJeOTkvVLKGVLKYCllsKen59PstkCoUtYX/y4jsPKvfbcwJAROnIAikOgeplqjNgyra4Vz+BLq+9kz8edlkBYPF4vNlUFFKXAedftoOSFEYynlO8B0oLrxsQuY8YhtXwJ8syz7GMtucwSqAZuFEBFAA2C1EKJY9Fhu/XwfLpFlALYmTbR/d+40T0D5xcWPqBQLak86x4+nSiIDmvH2P2n88O035o5MUYqtR50RfAvEAUgpf5NSjpJSjgJWGuseZh9QXggRIISwAnqijWCKcXuxUkoPKaW/lNIf2A10Lg53DQFUcLPk+P49XI+J0wrq1tXOCowDtRVZOj2N5yTx31UDGW6BCBtHDt6y4/SJw+aOTFGKrUclAi8p5X3/Q41l/g97oXHymuFoM5sdB5ZLKY8KIcYLIQr/WApPKePSUa4uHsOvf2/WCqytYetWePZZs8aVH74Z0ICtA+14c/RYANZ/9yaTmyRAaryZI1OU4ulRicDlIXW2j9q4lHKNlLKClDJQSvmZsexjKeV9vZKllM2Ly9kAwIsdW1Oi21gyXP2yV6SmQhGf7P3Z7n0JqVcLSmrTdoqApiSlZbBrzodwQbUVKEp+e1QiCBVCDLq3UAjxKhBmmpCKh0CfEgTWDuFkjOFu4c6d4Oxc9O8eqvsKDN2u3S0FSJ969FmZRvPh38HsthBz3swBKkrx8qhEMAJ4SQixWQjxjfGxBXgFrROY8hQCrBP5Z/UvdwsqV4a0tKKfCO4hrO0ZO3UR6+d8gZQGOLrS3CEpSrHy0EQgpbwmpWwEjAMijI9xUsqGUsqrpg+vaEs5tZOTSz8n/FykVuDqCtWqaZ3LipmarV6gab8xLLpYmn7/94W5w1GUYiVX8xFIKTcBm0wcS7HzyksD2G0IJDLFgoq3C0NCYP58rZ3AIrfTRRQd1+wrcebSHyRdPIKdbzVzh6MoxcLjTFWp5LFmNcpj7ebNwYuxdwubNIGEBDhcPG+nHPW/79j5ij12p37XCi7uhbjLZo1JUYo6lQjMyN7aAvfYkyxfNPduYfPm8MUXUAR6UD8J4eILlTsTtX4K5//4iiUjWzCgXbHoY6goZiMK2+yTwcHBMjS06NxlWrNNN47u3kRq7A10OpWXATJvXaJ6RX86BOrwdrFi6aEkNvy7GceKIeYOTVEKLSFEmJQyx19V6pvHzN5492O8h8wiIjrLGP2xsfDHH2AwPPiFRZjepTTj330LnZ0rb07byJ5hpXA8Mt/cYSlKkaUSgZk1ruaP0Fty4OKtu4WrVmkT2RfTdgKAF0ZN5MsNV9GXqQu1+pJ5+FcyYtWNaopiCioRmFn5Eo6k/Lea2bNm3i1s0UL7d+NG8wRVwByxCMJtQgx/Lfje3KEoSpGkEoGZ6XUCeT6U/3ZtvVvo6wvlysEmdccuQGD9DvSpbo2PuKZdLrtyyNwhKUqRohJBAfDm13Oxa/9/pKRnmbaxZUvYsqXIjzuUG7ZOrvw4qBF17C7BwcUwPYTUU8Wv052imIpKBAVAnQAP0jMlx67E3S1s2RLi4uC//8wXWEHiU4/Lx/cRt2M2dWYk8M7bI8wdkaIUGSoRFACVS9gR9ddkps2ce7ewfXs4cgSC1T30ACtPGij9dRTOb2wiwyCo73AZDGrie0XJCyoRFAB+JZyRN89z7GyWKZ6dnKBqVRDCfIEVIA3avsjIBlZ0rWTB71PepU+FJAzntpN58zxkpps7PEUp1FQiKCBe+noZGVU6kq2DX2goDB6szVFQzJWqWIdJ3QP57Z3WBDz3HrEZ1tRp0orNI6rAZ6Xgn4/MHaKiFFoqERQQ9QPcuBqXwsWbyXcLIyPh559h717zBVZQCAG9l0HXaWDtwKKUltSqUh63NqM4ZAjki8k/IotpBzxFeVoqERQQgY4Grix4m2+mzbpb2KyZ9gWo+hNoSgaBcUa3YROXMvvf49Tq9z92plbkw3U3iTy608wBKkrhpBJBAVGvkh9Wtg6cj8lyvdvVFWrXhvXrzRdYIdDn1eHEvOuIr+GSuUNRlEJJJYICwsJCz4sf/kRcydrZK9q3h127ICbGPIEVAo4BdXBydIKL2nzH0mBg/dwvSY67aebIFKVwUImgAKkb4MaZG/Fcjkm4W9i+Pfj5wblz5gusoNPpWXW9NGOmLgfgz+njeealMRz/vpe6xVRRckElggLES94kcmpffpiz9G5hw4Zw5ox2iUh5oLCb9izdc5n0+Gjq6o/yRWs7gpJ3wdrRUMiGWleU/KYSQQHSpl4QjhUbci3T7m6hENojM1N9oT3ER++PIWKEI5b7Z1Py2hbGjHwdiyZv8OGX3/PNG13NHZ6iFGgqERQgDnY2tB/6CVesfbNXbN8OJUvC/v3mCawQsKzQCsq15tcfxvL11jhk9R7IVuM4mV6Kk7vWIDd/CZlq3CZFyYlKBAVM3QA3jpy+wPVb8XcLK1SA6GhYs8Z8gRV0Oj28MJMvdsPcozrwroXOwoJF/x5m2pg+iM2fw4oB5o5SUQoklQgKGJvoU1yY2pc5v/x1t7BECahbVyWCR7F1Zduh8+wMO4YwTvtpaeeI6D6HNbpnqDNqGXHXI80cpKIUPCoRFDDdnmmCW/OBxFh4ZK9o3x727IGoKPMEVkjYunjiXMr/vnK7wEa42ULUETXHg6LcSyWCAqaEmzNNu79KeJJt9ooOHbTG4n/+MU9ghVzzbq+yvr8TZbkASTdh60SY1ZbkvQvNHZqimJ1KBAVQQ39n9uzcQeT1LB2igoNh9GhtRFLl8Vk7QslqyPM7+V//EA4t+YQr509TrV1/Zn040NzRKYpZqURQADnGnePKonf5eemqu4U6HXz5JdSoYb7ACrkfDlqje3kNX64+xj+iBa7vhNGgsi81rv8CO74zd3iKYjYmTQRCiHZCiHAhxGkhxJgc6kcJIY4JIQ4JIf4VQviZMp7Col/nNpTu9iFJ7pWyVxgMsHMnHD1qnsAKuQq1GjEs2JIjb5dh2OezsXF0YdHW0wS36Q7rP+Kfz3qSnpJk7jAVJd+ZLBEIIfTAD0B7oArQSwhR5Z7V/gOCpZTVgV+Ar0wVT2HiYGdDmw6dCLtyzzwE6elaW8E335gnsEKuTa/X+aGLG/5dP8TOxVMr1FvCCzM5Xqo77T5cxsS+tdWwFEqxY8ozgnrAaSnlWSllGrAU6JJ1BSnlJinl7Z9guwEfE8ZTqNQsYcHBf38j9NiZu4XW1tC5M/z+u5YUlMfjUAL+7yTUG5S9XKen8mszWfnVEN6qeBlOb2DmR69gb6Xj2KchcGGPeeJVlHxiykRQGsgy9yKRxrIHeQVYm1OFEGKwECJUCBF648aNPAyx4KrobODm31OZsfi37BXdu2sjkao5Cp6MjdMDq7qM+g4711LI7VPoIP+lZQVHyltegz/eUsN7KEVagWgsFkL0BYKBr3Oql1LOkFIGSymDPT098zc4M3mmUW1qj5pFRmCz7BVt2oCjI6xYYZ7AijK9JXvsWlF59DoSY6P44+9/0bcZyx/bD7H7l6nmjk5RTMaUieASkHXQHB9jWTZCiNbAB0BnKaWanNdIp9PxbNN67DwbTUp6lmvWNjba5aHNm9WvVBPwe2YIGRKueTQCn2AyKnXhtTVpTP1yHHwbxP7v1DAVStFjykSwDygvhAgQQlgBPYHVWVcQQtQCpqMlgesmjKVQqlPSkivrfmbGinuGlpg0SbtzSAjzBFaElSxXnZUrllN9qDZlqJWdAxunjWFu+3R+ORBLnbfms2H2ZwAkx8WQGPGfOcNVlDxhskQgpcwAhgPrgOPAcinlUSHEeCFEZ+NqXwMOwAohxAEhxOoHbK5YalrFh4TD61mz6Z65eEuU0BqOFZMIavkCTiXunsxW6vk/LN8+Sucfj/HdiwE0jfyRqxOb4OHhzuyhjWD3NDNGqyhPT8hCdnkhODhYhoaGmjuMfPPyrJ0cv57MjndbotNlOQNYuRK++AK2bVNJIT9Fn4GNn0JyDOPWXaNjRTuCDWHw4jyo0uXRr1cUMxFChEkpg3OqKxCNxcqDdaxZhiuxKRyMvJW9ws4O9u2DP/80S1zFlnsgdJ8D/X9n7KJdBL/3J3/f9KVT975kpKbA3p9hxxTVfqMUKioRFHCtKpUgatUXvD3mg+wVrVuDtzfMm2eewBSNpS3J5TtzPiqJE0s+RK4ZTczqjxjXtzEZZ7ZCrBr2Win4LMwdgPJwLvZWlPRwI+JWBlJKxO0GYr0e+vWDiRPh2jXw8jJvoMVYl6Gf8FzGKp6bMIV0qadds/p8sngrb/s9i4OVYF18BaoNmkbpSnXMHaqi5EidERQC4yd+j6F6F45dicteMXCgNpfx/PlmiUvR6Cwsoe4gOlWwoHOnjrwxbROhaxdiO/A3Mpq+z8fLDvD1mNfMHaaiPJBKBIVAu2olsdAJFm85kr2iUiV4911tiGrFrET9wbz60fcMmbAYodNRp10f9BVaY9HyXV7v2ohhZS9Cqjb9aNT5kxzbptp2lIJDJYJCwM3eCqs9c/l6yHOkpd8zAfuECdCihXkCU+6ytIW6r4DF/Xdw9R89kQpOKcgDSyA1gZPfdSWkbWcWjWwDybfyP1ZFuYdKBIVE7x4v4Fi/GztPXbu/8vx5WLQo/4NScscnmBO6ivR9bQTzh9Sljt0lhnWsRXvrUFjUHdISzR2hUsypRFBIvD2wGyUbdGHNsRzmLP7pJ+jfH86dy//AlEcTAtH2c5KtPPDOvIB120/434ow3PrPI+38PrrV9eb3V8vCqfXmjlQpplQiKCRsrfS0q+LB4mUruHQ9Onvl8OHaDGZTppgnOOWRKtZvzW97L9B6xmVoMkIrrNKZXWVHkYgtVoZkbR5lRTEDlQgKkZp2t7j0y2eM+2529gofH+jZE2bNglu3zBKbkks2ztkWm730MWsPXqXDq+9z/vBOIkPXcWrxaPo3r8ji97reaWBWFFNSiaAQ6fNsS+oOm0ykR13uGxpk1ChISIAffzRPcMpTSSz/HFV+TOCL4S+QuPVH9h6PoHnaRpjRAnnlkLnDU4o4lQgKEZ1Oxxv9unLsaiIHLt7KXlmrFjz7rNa5TCl07D19mTfiGUY30FHzlSmcuJaK9/A/ITWOnm3q8vlLLSAxh/YhRckDatC5QiYhNYOK3d+hRGYU//21IHtlZqbW41gpnNISIeY8eN2d2js99hpDXmhOoDzH+y2coe3n2m2qivKYHjbonBpiopBxsLYg0DKWfQcOEHEjDn/PLFMv3k4CYWFQvjw4PXhaRqUAsrLPlgQALJ29mLXhOPL6CVj3Hlt+HEFG7YO0sjrAmehM/rNtxAtvfYnQqZN75cmpT08htHD6t5Tu+yXzd+cwoNm5c1CvHnz5Zf4HppiMKFEJQ7e5jNxgYN3S6ZB0k3lbz9L3nYmc/FNNo6k8HZUICqEyHk50ruHNwm0nOH3xnjaBgADo0UObxSxSjXxZlOhsnVm1diMVQ7rC0J189PtZDn9SnwrHJrNswht80KsxhpvnzR2mUgipNoJC6vC5q9QKqkz9Nl3YsXJu9sqICG0coueeg6VLzRCdkm+un4DpIdSddpOO5S34pH0pLpbtzdoDV+jcvTclg1qAXl0BVtTENEVSUEBJQroNItK1FtfiUrJX+vvDe+/BsmWwXvVWLdJKVIK+v7F6xWI+Xn4IXPw4tPJbXvt8LpE/doXpTeHWBXNHqRRwKhEUYnO/GYtVqfJ880/4/ZXvvgs1asDFi/kfmJK/AkIo1agHuhIV4bUtPDPzKuf+20qNV3+A2Ei+7lObk7vXmTtKpQBTiaAQ83O3p1dNT6ZPHM/yv7dkr7Sx0e4eevll8wSnmI2lnRP+NUOwrDeA2Bd/46stMcz4YAAHFo2jjKsVi997HlITzB2mUoCoRFDIDWoaQNLRTXw2ayUGwz3tPbdvJ122DNauzf/gFLNzDgxm75rFfNUklWrh31DWw4pGyevhqwA2vV2LKSO6k3TzirnDVMxMJYJCztfLgxm/byE2sA0Ldudwx0h6OnzxBfTuDWfP5n+AitkFNO2BrvssLNpPYHN4LP7/twnqD+G/C/F8PO0XDNNbQOwlc4epmJG6a6gIkFIycM4+toceYt6rDWlSu2r2Fc6c0WYxK1MGtm1THc2UO67s+4NS64dicPLlSMURVK5RG8uYM0ivqggXX3OHp+Shh901pBJBEXH2agwVy5fDtUwFLhzYjo3lPUNNrFsHHTtqs5n99RdYWZknUKXgObOJT17tyLjNyZx7ywF/Fx2Td6Wx/qozS1f+hVPZOuaOUMkD6vbRYqBsSVcm/DALqxbD+GT10ftHJ23bFmbOhA0btESgKLcFtmDQtB38Onk0fn2+hYF/4Vi9A44yHscFbWB2O6a/1ZHUxDhzR6qYiDojKGK+XneC7zee5oXSCUx6s+f9K+zbB3Xras+lBCHyN0Cl8EiMgm2TyLgYhtdb6xncMpAvJn3P1StXSXEth1+1hmqMo0JEnREUI2+3qUjFuDAmv9WLD6f9cv8Kt5NAWBg0bar6GSgPZu8B7T7HYtA6Zn02gveCbsDiF/n5k9cIrNmEqxOCYef3yPjr5o5UeUrqjKAISkhKod3r47lYohGj2lTkzVblEPf+8v/nH+jWDSwttUtGXbuaJ1il8LhxElJiORlxkcPb/+YFjzNwKZSXV6eQauXGok8GQs0+ULKauSNVcqAai4uhjEwDo389xPKth3E48ivrFs3Ar3SJ7CudPAm9esH+/VpS+OYb7c4iRcmt6ycY/84Q5K1IxtaNh4wUqvyUTLPyTvz0Tk8I6g4BTdUlyAJAXRoqhiz0Or7pXoNOXnGc3PUPz05cw+qDl7N3OqtQAXbtgk8/hT/+UAPUKY+vRCU+nreZsatOw6jj0PZz2jcKokuTIDi2mtgZnWhRyZXd41rB6jeI3z6T0DULSb6lLicVJOqMoBjYdvgsX22+zOFLsVge+IUX2zXl4+H9sbbIcovphQvg5gYODrBkCSxfrg1n3akT2NubL3il8MpIJXTR/xg+biqfdy5Dy5KxbDgWTZsFSWzsb0eL6mWIsA1ia7Q7z700AidnF9IzMrBwLqUaoU3AbJeGhBDtgCmAHpgppZxwT701MB+oA0QDPaSUEQ/bpkoETybTIFm+6zSvdG2FZWB9AtoPoml5D2wi9zG0d1cCsl42mjoVJkyAy5fBzk7re9CiBbz9tvnegFL4GQxEndzL1n9W0bKcAy5JZ5m6cDVvroriwggHfJ11TN6Vykeb04j8LBiX0uXYdlnPoauZvPbyACw8A0m1cMDC3g29peoH87jMkgiEEHrgJNAGiAT2Ab2klMeyrDMMqC6lHCKE6Al0lVL2eNh2VSJ4OhkZmfx7NJJ14TGs3bybEz8Nxf3Zt6nQqD2+Nilc2fEbHbv3pUblitgd+4/41b9R6ehhPOytcdi+Fb1OaO0JCQnaJDglS2qPKlUgJETbyYULWoc1OzvtbELNo6w8gCEjg2PbVlHN9gboLdgSepy/Nu7gy+f8ELciGLX4MDP2JRP/niNCCN5cm8KCQ2nEfOQNdm5M2plM2KV0Fo1+FuzcWb33PNfi0xj0QmuwtGPviYskpUPzxvXAwpZzl6OQemvKlisHFrYkpWWis7TBxtYWEKRnZCB0eiwsLEEI7Q5rna5ItHGYKxE0BD6RUrY1Lr8HIKX8Iss664zr7BJCWABXAU/5kKBUIsg76RmZLPlrI5czHbmSZs3OLZsInfEuXr2+wManMsln9nH9l3GU7PcN1t4VSTq1m6hVE/DrNxEHNx8SwncSuXEmVbp+gENJP2JunODcmhlU7zASG1tnrp/Zx/l9q6jZeTRWJX24FnmI8xvmU7vTO1ha23P1xHbO/7eW4O4fY+FRgsund3FxywrqPv8BOr0llw5vJPLwRur3+h86V1cuHN3IlT1rqN/1AwQQ8d/fXDuzj/rdPgJ3d86F/UXU0R3U7aiduZzZt4qYS+EEPzcaPD05vWMFcRFHqf3M6wCc2rWChOhL1Hp2BHiVIPzfeSTfuETNFq8CcGLrItKS46jediiUKsnxNdPISIwnqEk/AI5unA1SUrXVK+DtzZGVk9Cho0qDFwE4tG4aVjYOVGrWF0r7cHDJp9jYu1KxThcA/vtzCg5u3pRv1B18y7B/zvs4efpTrkY7AMJ+/xoX7woE1usCfv7s+2kE7n7VKVu1JQB7f/mUEmXr4F+7PQQEsHvyYEpVCcGvQhMAdi35GJ9qzfENaonBz489kwbhW+sZfMrWIzM9lb2/fIZfzWfwrhxChrcX+6YOx7/us5QqU5O05HjCfv+asnU74VWhAalujuyf/n8ENnieEqWrkBwfzYE/p1Cu4Qt4lgsmyU7PwTkfUqFJT9y9ypFw8zKH102jYtPeuAXUIF6XwpFFn1G5eT9c3P2IvXaOYxvnUKXlQJzLVOFW+i2Or5hI1Vav4ORSipuXThC+dTHV2gyipIcdbuIS+kMr0LuUJS4+ng4V7dDFnGTbFStOX09kYTdHnElgwO9JHLpm4PjrDgB0XZbEmZsGDg3VljsuTuJagoHQwdpymwWJJKfD9pe1y58hcxKx1MHGAdpy/ZkJuNoI/u5rj0FCvZ8T8XXWsay7HVJKGsxMpJKnntld7Mk0GGg0K4HapSyZ2sGWTIOBJrMTaeJvzeetbTFkZtB0dgKty9nwYVNrDAYDzefE82wVe0Y2tMaQmUnruXG8EGTPa7UtMBgkbRfE0au2C/2CdCB0xDYfT52Orz7R/3dzTV5fGsh6k3okUP9B60gpM4QQsYA7EJV1JSHEYGAwQBl1V0uesbTQ079Lm7sFPWuRMfUtbialcT0+jeNnfNlX14sKtZuAtT0nj2eyz+IlGrSqha2zBxdKpbM//SKNG1TF2tGV8xEJGCrVoEYZd+ys7bmQ7In0DSTI3RprD2vO670xVKhGNQcdlhaZODjbk+nlSxWrdCyc9Nj4epMZWInKulR0ugys7G3I9ChFJV0KwkGHrrQ30q8cFfSpABgc7NC5laC8PhUc9WT4lMIiLkBbBtIc7LFx9dCWnfSk+ngTnRFPeb02kU+yowOxmW7aspMlib6+JFrpKK/T6hOcnEix0mvLzlbE+fiQERd3p/6WkzNSoi27WnHTtwz69AzKC60+yskFa2tbbdnNmhu+ftjqbe7UX3Nxx9HBwVhvw5UyAbg6uN2pv+zigYe9nbbsbktkmQBKetytv+BagtL2xu172BNRJpAybi536s+6eVHG1ppyIgVDCQci/AIp4+JEoEghQ5fGeTcvythaEaBLIc3LkQt+gfg5OeAnUkjVpRPp5oWfjQW++lSSS/pwyS8QfydbSosUEvUZXHHzwt9aj7dFGvElS3HVLxB/eytKihTiLAxcd/MiwFJQwjKdGE9XbvgFEmBrgYdI4aYlRLt5EWAJ7lYZRHm4cdMvkLI2OlxFCjesdNxy8yLQCpxdvLlm5cWB1KM08mqDSxkHfr18hkOHUmnc5AXKNy/L4Mw4jmxYQUj9lrS2hAGHz3DxzF6q1KpPgKsTo1Pg+uFtVKxeikoyk0ln04m5cYaqlaqArTWzk12JvhBO7Sp26HU65l3Uc/PmFepVLYu1gz0LEh2JuxJBvSqWOFhb8MslHQnx0dStVBpXd2dWJduQFHWN4AoZeDnZ8vc1SEmOp3agO+6l3dmSbEFqfAw1ApJxcXdgd7SBtLRkqpWxx9GrBAfSBBnJiVT0tsTOwZXjcRkYDBkEltBh7Vqac5kZ6HR6PB3cTfJdYMozgm5AOynlq8blfkB9KeXwLOscMa4TaVw+Y1wnKqdtgjojUBRFeRLmun30EpB1+EIfY1mO6xgvDTmjNRoriqIo+cSUiWAfUF4IESCEsAJ6AqvvWWc1MMD4vBuw8WHtA4qiKEreM1kbgfGa/3BgHdrto7OllEeFEOOBUCnlamAWsEAIcRq4iZYsFEVRlHxkysZipJRrgDX3lH2c5XkK0N2UMSiKoigPp7rvKYqiFHMqESiKohRzKhEoiqIUcyoRKIqiFHOFbvRRIcQN4PwTvtyDe3otFxAqrsej4np8BTU2FdfjeZq4/KSUnjlVFLpE8DSEEKEP6llnTiqux6PienwFNTYV1+MxVVzq0pCiKEoxpxKBoihKMVfcEsEMcwfwACqux6PienwFNTYV1+MxSVzFqo1AURRFuV9xOyNQFEVR7qESgaIoSjFXZBKBEKKdECJcCHFaCDEmh3prIcQyY/0eIYR/lrr3jOXhQoi2+RzXKCHEMSHEISHEv0IIvyx1mUKIA8bHvUN4mzqugUKIG1n2/2qWugFCiFPGx4B7X2viuCZniemkEOJWljpTHq/ZQojrxsmUcqoXQojvjHEfEkLUzlJnkuOVi5j6GGM5LITYKYSokaUuwlh+QAiR5zM95SK25kKI2Cx/r4+z1D30M2DiuN7JEtMR42fKzVhnkmMmhPAVQmwyfg8cFUK8lcM6pv18SSkL/QNtmOszQFnACjgIVLlnnWHANOPznsAy4/MqxvWtgQDjdvT5GFcLwM74fOjtuIzLCWY8XgOB73N4rRtw1vivq/G5a37Fdc/6b6ANb27S42XcdlOgNnDkAfUdgLWAABoAe/LheD0qpka39wW0vx2TcTkC8DDj8WoO/Pm0n4G8juuedTuhzZFi0mMGlAJqG587Aidz+P9o0s9XUTkjqAecllKelVKmAUuBLves0wWYZ3z+C9BKCCGM5UullKlSynPAaeP28iUuKeUmKWWScXE32kxuppab4/UgbYH1UsqbUsoYYD3Qzkxx9QKW5NG+H0pKuRVtzowH6QLMl5rdgIsQohQmPF6PiklKudO4T8i/z9btfT/qeD3I03w28zqufPl8SSmvSCn3G5/HA8fR5nPPyqSfr6KSCEoDF7MsR3L/gbyzjpQyA4gF3HP5WlPGldUraFn/NhshRKgQYrcQ4rk8iulx4nrBeBr6ixDi9rSjBeJ4GS+hBQAbsxSb6njlxoNiN+Xxehz3frYk8I8QIkwIMdgM8QA0FEIcFEKsFUJUNZYViOMlhLBD+0L9NUuxyY+Z0C5Z1wL23FNl0s+XSSemUXJPCNEXCAaaZSn2k1JeEkKUBTYKIQ5LKc/kU0h/AEuklKlCiNfQzqZa5tO+c6Mn8IuUMjNLmTmPV4ElhGiBlgiaZCluYjxWJYD1QogTxl/L+WU/2t8rQQjRAfgdKJ+P+3+UTsAOKWXWsweTHjMhhANa4hkhpYzLq+3mRlE5I7gE+GZZ9jGW5biOEMICcAaic/laU8aFEKI18AHQWUqZertcSnnJ+O9ZYDPaL4V8iUtKGZ0llplAndy+1pRxZdGTe07bTXi8cuNBsZvyeD2SEKI62t+vi5Qy+nZ5lmN1HVhJ3l0OzRUpZZyUMsH4fA1gKYTwwMzHK4uHfb7y/JgJISzRksAiKeVvOaxi2s9XXjd8mOOBdmZzFu1Swe0Gpqr3rPM62RuLlxufVyV7Y/FZ8q6xODdx1UJrHCt/T7krYG187gGcIo8azXIZV6ksz7sCu+Xdxqlzxvhcjc/d8isu43qV0BruRH4cryz78OfBjZ8dyd6Yt9fUxysXMZVBa/NqdE+5PeCY5flOoF1eHqtcxFby9t8P7Qv1gvHY5eozYKq4jPXOaO0I9vlxzIzvez7w7UPWMennK0//8OZ8oLWqn0T7Uv3AWDYe7Vc2gA2wwvgfYy9QNstrPzC+Lhxon89xbQCuAQeMj9XG8kbAYeN/hMPAK/kc1xfAUeP+NwGVsrz2ZeNxPA28lJ9xGZc/ASbc8zpTH68lwBUgHe067CvAEGCIsV4APxjjPgwEm/p45SKmmUBMls9WqLG8rPE4HTT+jT/Iy2OVy9iGZ/l87SZLssrpM5BfcRnXGYh2A0nW15nsmKFdspPAoSx/qw75+flSQ0woiqIUc0WljUBRFEV5QioRKIqiFHMqESiKohRzKhEoiqIUcyoRKIqiFHMqESjFghDCPcuokleFEJeMzxOEED+aaJ8jhBD982A7S4UQBanXrVLEqNtHlWJHCPEJ2kilE024Dwu0YRRqS21sq6fZVjOgr5RyUJ4Epyj3UGcESrFmHBf/T+PzT4QQ84QQ24QQ54UQzwshvjKOQf+3cRgAhBB1hBBbjIOPrTOOAnmvlsD+20lACLFZaHMphAohjgsh6gohfjOOIf+pcR17IcRfxoHYjgghehi3tQ1obUwuipLnVCJQlOwC0b7EOwMLgU1SyiAgGehoTAZTgW5SyjrAbOCzHLbTGAi7pyxNShkMTANWoQ17Ug0YKIRwRxvt8rKUsoaUshrwN4CU0oDWa7RGnr5TRTFSvzAUJbu1Usp0IcRhtElS/jaWH0Ybo6Yi2pf3em06C/RoQxbcqxTauPJZ3Z417TBwVEp5BUAIcRZt4LDDwDdCiC/RJm3ZluW11wFv7k8uivLUVCJQlOxSQfsVLoRIl3cb0Qxo/18E2pd4w0dsJxltfKv7tm3cVmqWcgNgIaU8aZyCsAPwqRDiXynleOM6NsZtKkqeU5eGFOXxhAOeQoiGoA0fnGVSlayOA+UeZ8NCCG8gSUq5EPgabUrF2yoAOc6zqyhPS50RKMpjkFKmCSG6Ad8JIZzR/g99izYiZVZrgQWPufkg4GshhAFtdMyhAEIILyBZSnn1aWJXlAdRt48qiokIIVYCo6WUp55yOyOBOCnlrLyJTFGyU5eGFMV0xqA1Gj+tW2hThSqKSagzAkVRlGJOnREoiqIUcyoRKIqiFHMqESiKohRzKhEoiqIUcyoRKIqiFHP/Dwj1bc9ci71FAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "t2exp = 0.27 # Experimental T2 from PhysRevB.82.121201\n",
    "decay = lambda t: np.exp(-(t/t2exp)**2.4)\n",
    "plt.plot(time_space, decay(time_space), color='red', label='Experiment', ls='--')\n",
    "\n",
    "plt.plot(time_space, shallow_l_cce.real, label='Shallow')\n",
    "plt.plot(time_space, shallow_l_gen.real, ls=':', c='black')\n",
    "\n",
    "plt.plot(time_space, l_cce.real, label='PD')\n",
    "plt.plot(time_space, l_gen.real, ls=':', c='black')\n",
    "plt.legend();\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('Coherence');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interesting to note - the decay depends significantly on the orientation of the magnetic field. You can check it yourself!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.3"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}