doc/source/examples/t2_non_markovian_noise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## $T_2$ Dephasing with Non-Markovian Noise\n",
    "In this notebook, we show the effect of non-Markovian noise sources on qubit\n",
    "coherence and their implementation in qopt. In terms of noise spectral densities\n",
    "(NSD), these noise sources are characterized by non-trivial frequency\n",
    "dependencies. While Markovian (or white) noise has a flat NSD, a typical\n",
    "frequency behaviour of a non-Markovian SND is described for example by a power\n",
    "law (referred to as colored noise).\n",
    "\n",
    "We will demonstrate the effect of such a noise source based on Hamiltonian\n",
    "motivated by the computational states of a singlet-triplet spin qubit with\n",
    "basis $\\{|S \\rangle, |T_0 \\rangle \\}$. The basis states correspond to the spin\n",
    "singlet and triplet zero states\n",
    "\n",
    "\\begin{align}\n",
    "|S \\rangle = \\frac{1}{\\sqrt{2}} \\left( |\\uparrow \\downarrow \\rangle\n",
    "- |\\downarrow \\uparrow \\rangle \\right) \\\\\n",
    "|T_0 \\rangle = \\frac{1}{\\sqrt{2}} \\left( |\\uparrow \\downarrow \\rangle\n",
    "+ |\\downarrow \\uparrow \\rangle \\right)\n",
    "\\end{align}\n",
    "\n",
    "where $|S \\rangle$ is defined as out $|0 \\rangle$ state.\n",
    "\n",
    "The Hamiltonian we employ is\n",
    "\\begin{align}\n",
    "H = \\frac{J(\\varepsilon + \\delta \\varepsilon)}{2} \\sigma_z +\n",
    "\\frac{\\Delta B_z}{2} \\sigma_x.\n",
    "\\end{align}\n",
    "\n",
    "Here, $J$ is the exchange coupling between the involved spins controled by an\n",
    "control parameter $\\varepsilon$. $\\Delta B_z$ describes a magnetic field\n",
    "difference between the two spins. $\\delta \\varepsilon$ is noise on the control\n",
    "parameter whose one-sided NSD we assume to follow $S_{\\varepsilon}(f) = S_0/f$.\n",
    "\n",
    "A common approximation for the exchange coupling that we will use is\n",
    "$J(\\varepsilon) = J_0 e^{\\varepsilon / \\varepsilon_0}$.\n",
    "\n",
    "As we are dealing with correlated noise sources here, we cannot use the\n",
    "Lindblad master equation but perform Monte Carlo simulations.\n",
    "If we wanted to optimize control pulses, we would use generalized filter\n",
    "functions for the calculation of infidelities.\n",
    "Furthermore, as\n",
    "our noise acts on a control parameter instead of a control amplitude, we\n",
    "setup our solver with the according qopt solver class that incorporates this\n",
    "usecase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "from qopt.matrix import DenseOperator\n",
    "from qopt.solver_algorithms import SchroedingerSMCControlNoise\n",
    "from qopt.amplitude_functions import CustomAmpFunc\n",
    "from qopt.noise import NTGColoredNoise\n",
    "\n",
    "pauli_x, pauli_z = DenseOperator.pauli_x(), DenseOperator.pauli_z()\n",
    "\n",
    "np.random.seed(0)\n",
    "\n",
    "S_0 = 1e-3\n",
    "n_traces = 1  # The stored data was simulated with 500 traces. We plot\n",
    "# the stored data and calculate only dummy values to comply with computational\n",
    "# time budget\n",
    "\n",
    "time_step = 0.1\n",
    "\n",
    "delta_Bz = np.pi / 2\n",
    "J_0 = 1\n",
    "eps_0 = 1\n",
    "\n",
    "def exchange_coupling(ctrl_param):\n",
    "    return J_0 * np.exp(ctrl_param / eps_0)\n",
    "    \n",
    "exchange_coupling_ampl_fctn = CustomAmpFunc(\n",
    "    value_function=exchange_coupling,\n",
    "    derivative_function=None\n",
    ")\n",
    "\n",
    "def pink_spectrum(f):\n",
    "    \n",
    "    return S_0 / f\n",
    "    \n",
    "def average_monte_carlo_traces(solver, initial_state, measurement_state):\n",
    "    \n",
    "    propagators = solver.forward_propagators_noise\n",
    "    probabilities = np.zeros(n_traces)\n",
    "    \n",
    "    for i in range(n_traces):\n",
    "        final_state = propagators[i][-1] * initial_state\n",
    "        probabilities[i] = np.abs(\n",
    "            (measurement_state.dag() * final_state)[0,0]) ** 2\n",
    "\n",
    "    return np.mean(probabilities)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Free Induction Decay (FID)\n",
    "The effect of $1/f$ noise becomes apparent when performing an FID. Starting in\n",
    "the $|0 \\rangle$ state, the procedure consists of two $X_{\\pi/2}$ pulses that\n",
    "are intercepted by a waiting period of lenth $\\tau$. During the latter, the\n",
    "qubit precesses around the z-axis with frequency $J$ and is sensitive to noise\n",
    "on $\\varepsilon$. At the end, the probability of the $| T_0 \\rangle$ (or the\n",
    "$|S \\rangle$) state is measured which reveals an approximately Gaussian\n",
    "$T_2^*$ decay following\n",
    "\n",
    "\\begin{align}\n",
    "P(T_0) = \\langle |\\langle T_0 | \\psi_f \\rangle |^2 \\rangle_{\\delta\n",
    "\\varepsilon} \\propto e^{-(\\tau / T_2^*)^2} %*\n",
    "\\end{align}\n",
    "\n",
    "where $|\\psi_f(\\tau) \\rangle$ is the state after a FID sequence with waiting\n",
    "time $\\tau$ and $\\langle \\cdot \\rangle_{\\delta \\varepsilon}$ is the average\n",
    "over an enseble of noise traces.\n",
    "\n",
    "In our simulation we use an idealistic pulse with vanishing exchange coupling\n",
    "$J$ during the $X_{\\pi/2}$ rotation. This is achieved by setting $\\varepsilon$\n",
    "to a large negative value. The $X_{\\pi/2}$ rotation is mediated by $\\Delta\n",
    "B_z$.\n",
    "\n",
    "To solely focus on noise related effects, we let $\\Delta B_z$ vanish during the\n",
    "precession periods. This is done by setting it to zero in the drift\n",
    "Hamiltonian. Alternatively, we could use a $\\Delta B_z$ that is negligible\n",
    "with respect to the pressesion frequency. This would result in much longer\n",
    "pulse sequences increasing the computation time drastically.\n",
    "\n",
    "Another important consideration to be made when constructing the noise trace\n",
    "generator is to specify the $\\texttt{low_frequency_extention_ratio}$ parameter.\n",
    "This will artificially increase the sequence's length during the generation of\n",
    "noise traces by multiplying with the specifed value. As the low-frequency\n",
    "cut-off of the noise spectrum is determined by the sequence duration, the\n",
    "cut-off frequency is decreased. In this way, we can include noise that is\n",
    "already quasi-static with respect to $\\tau$ but still fluctuates on comparable\n",
    "time scales. We chose a factor of 500 as the $T_2^*$-time roughly converges at\n",
    "this point.\n",
    "\n",
    "The computational time budget of readthedocs does not allow us to compute the\n",
    "proper simulation remotely. Therefore we only use a single noise trace for the\n",
    "calculation of dummy values and plot data calculated with 500 noise traces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAiYElEQVR4nO3de5hcVZnv8e+vqglogEEhoiTBxBjFHATEFrnpiIiTeCH4eAui4G3yZBQdj+PR6HjGcZwZZXQYxzlgTg6HQUeRg2g4USOX8XgbAScdhwGCBNoApg2Y5g6CJF31nj/2rqZSVKd3Vdeld+/f53n66dqXteutSqXeXmvttZYiAjMzs8mU+h2AmZnlgxOGmZll4oRhZmaZOGGYmVkmThhmZpaJE4aZmWXihGHWQ5L+UtLX+h0HgKSLJP11+vhlkrZ08NqvkDTSqevZ9OCEYT0n6Q5Jj0l6pO7nEEkLJIWkgfS8iyTtlPRw+nOTpM9K+oNJrv88Sd+UdI+kByXdIOnDkspTjLsnX4KSFkqqSjq/289VExE/jYjn18Vwh6RX9er5LR+cMKxfXh8R+9b9bJ/gvL+LiP2AOcC7gGOBn0ma3exkSYuAnwPbgBdGxB8AbwYGgf06/iqe/PwDHbjMmcD9wApJe3fgemYd4YRhuRARv4+IjcCpwIEkyaOZTwPXRMSHI+KutOyWiHhbRDwAIOlUSZslPSDpR5JeUCuc/mX9kbRW8qCk/yNpnzRBfR84pKFW9JeSLpP0NUkPAe9M96+XdJ+kYUl/3OLLPRP4JLALeH39gbQG9j5Jt6W1rs9IWiTpWkkPSbpU0qz03FdIGpH0ibS2dYekM5o9YX3tSdK/AIcC30lf50eb1a7qayGSnpLWCO+XdDPwkoZzD5H0LUmjkm6X9MEW3xObBpwwLFci4mHgauBlE5zyKuCyicpLeh7wDeBDJLWWDSRfjLPqTnsLsBRYCBwBvDMifgcsA7Y3qRUtT5/zAODr6fVHgEOANwF/K+nkLK9P0suAecAlwKUkyaPRUuDFJLWtjwJrgTOA+cDhwOl15z4TOAiYC5wFrJX0fPYgIt4B/JonaoF/lyH0TwGL0p8/Sp+r9ppKwHeA/0zjOBn4kKQ/ynBdm0acMKxfLk//wn9A0uUtlt0OPH2CYwcCd+2h7FuB70XE1RGxC/gC8BTg+LpzvhQR2yPiPpIvuqMmiefaiLg8IqokX84nAh9La0XXAxcA75jkGjVnAd+PiPuBi4Flkp7RcM45EfFQRGwGbgKuioitEfEgSS3oRQ3n//eIeDwifgx8jyQhdtpbgL+JiPsiYhvwpbpjLwHmRMRfRcTOiNgK/C9gRRfisC5ywrB+OS0iDkh/Tmux7FzgvgmO3Qs8aw9lDwHurG2kX/Lb0mvW3F33+FFg30ni2dZw/fvSmlDNnQ3Xb0rSU0j6W76exnYtyV/6b2s49bd1jx9rsl0f7/1p7ag+lkMmi6UNh7D7+3Bn3eNnkzTl1f5AeAD4BHBwF+KwLnLCsFyRtC9Js9NPJzjlX4E37uES20m+wGrXE0lTzm8yPP1EUzvX798OPF1SfQf7oRmv/wZgf+B8SXdLupsk0TRrlsrqaQ03CByaxjiZxtf6O+CptY30jrM5dcfvInkf65+nZhtwe90fCAdExH4R8ZpMr8CmDScMywVJe0t6MXA5yR1E/zzBqZ8Cjpf0eUnPTMs+N+2UPoCkX+C1kk6WtBfwZ8DjwDUZwvgtcOCebutNm2OuAT6bdpYfAbyHtNYwibOAC4EXkjSDHQWcABwl6YUZyk/k05Jmpf0jrwO+maHMb4Hn1G3fCuwj6bXp+/ZJoP4OrkuBj0t6mqR5wAfqjv078JCkj6Wd42VJh0varWPcpj8nDJvuPirpYZImqK8Cm4DjG5pZxkXEr4DjgAXAZkkPAt8ChoCHI2IL8Hbgn4B7SO5Cen1E7JwskIi4haRDe2vatDJR087p6fNvB9YBn4qIq/d0bUm1zuAvRsTddT+bgCuo60Ru0d0kCXY7SdJalb6OyXwW+GT6Oj+S9o+8j6Q/5jckNY76u6Y+TdIMdTtwFfAvtQMRUSF5n49Kj9+TXmeP42ls+pEXUDKbmSS9AvhaRMzrcyg2Q7iGYWZmmThhmJlZJm6SMjOzTFzDMDOzTDoxUdq0ddBBB8WCBQv6HYaZWW5s2rTpnoiY0+zYjE4YCxYsYGhoqN9hmJnlhqQ7JzrmJikzM8ukpwlD0lJJW9Ipn1fv4byXSKpIelOrZc3MrDt6ljDSuWfOI5kieglwuqQlE5x3DnBlq2XNzKx7elnDOAYYTqdh3kky3//yJud9gGQqhx1tlDUzsy7pZcKYy+7TH4/QMOVzOp/OG4A1rZY1M7Pu6mXCUJN9jaMGv0iy8EyljbLJidJKSUOShkZHR1uP0szMmurlbbUj7D5f/jyePC//IHBJskQBBwGvkTSWsSwAEbGWZMlKBgcHPYzdzKxDepkwNgKLJS0kmR55BQ0riUXEwtpjSRcB342IyyUNTFY2bx54dCdfu+5Odo5V+x2KmfXI3nuVOev4Bey7dz6HwPUs6ogYk3Q2yd1PZeDCiNgsaVV6vLHfYtKyvYi7W66++bd84apbAVCzBjczm1Fq0/YtmjObpYfvaRXh6aunaS4iNgAbGvY1TRQR8c7JyubZrkry6fn5J07m4P336XM0ZtZtwzse4VXn/pidlfy2lHukd59U0j83Sq5emBVCuZT8X69WnTCsRZVK0ndR+xCZ2cxWTv84HHPCsFbVaqVOGGbFUC67hmFtqn1onDDMiqFWw6jkeNE6J4w+qVVLy+7DMCuEUvpt6yYpa1k1XMMwK5KBNGO4ScpaVnGTlFmhjDdJOWFYq2rVUucLs2KoNUk5YVjLqtWgXBJyH4ZZIdSapNzpbS2rRLjD26xAXMOwtlWqMf4BMrOZz30Y1rZKNcarqGY289VucHHCsJZVquEOb7MCkURJT9xSn0dOGH1SjfAttWYFUy7JNQxr3VjVCcOsaEpywrA2VJ0wzApnwDUMa0el6ttqzYqmVJLHYWQlaamkLZKGJa1ucny5pBskXS9pSNKJdcfukHRj7Vgv4+6G5LZaJwyzIsl7H0bPlmiVVAbOA04BRoCNktZHxM11p/0AWB8RIekI4FLgsLrjJ0XEPb2KuZsqEQw4YZgVipuksjsGGI6IrRGxE7gEWF5/QkQ8EjFeX5sN5PednYRrGGbFU5J8W21Gc4Ftddsj6b7dSHqDpFuA7wHvrjsUwFWSNklaOdGTSFqZNmcNjY6Odij0znMfhlnxlEtirOKEkUWzb8cnvXMRsS4iDgNOAz5Td+iEiDgaWAa8X9LLmz1JRKyNiMGIGJwzZ04Hwu6Oiu+SMiucsju9MxsB5tdtzwO2T3RyRPwEWCTpoHR7e/p7B7COpIkrtzxwz6x4yiV5AaWMNgKLJS2UNAtYAayvP0HSc5XO9y3paGAWcK+k2ZL2S/fPBl4N3NTD2DvOA/fMiqcs5XqJ1p7dJRURY5LOBq4EysCFEbFZ0qr0+BrgjcCZknYBjwFvTe+YOhhYl+aSAeDiiLiiV7F3g5ukzIqnXMp3p3fPEgZARGwANjTsW1P3+BzgnCbltgJHdj3AHqp6PQyzwsn7OAyP9O6TsYpvqzUrGs8lZW2peuCeWeEMlJ0wrA3uwzArnpJEjodhOGH0S7KAkhOGWZEkfRjVfofRNieMPql4HIZZ4ZTdh2HtqFRxwjArmGTgXr+jaJ8TRp9UPZeUWeF4ahBry1i16hqGWcGUSvke6e2E0SfVcJOUWdEMeC4pa4dvqzUrHg/cs7b4tlqz4imXcMKw1lWqHultVjQDpZI7va11lfBcUmZFU3IfhrUj6cPodxRm1ktl4bukrHVJk5TffrMiKZdK7sOw1lXd6W1WOOUSuV5AqacJQ9JSSVskDUta3eT4ckk3SLpe0pCkE7OWzZsxN0mZFU7ZA/eykVQGzgOWAUuA0yUtaTjtB8CREXEU8G7gghbK5koy+aAzhlmRlN3pndkxwHBEbI2IncAlwPL6EyLikYjx+tpsILKWzZuqaxhmhVOW55LKai6wrW57JN23G0lvkHQL8D2SWkbmsmn5lWlz1tDo6GhHAu+GMU8+aFY4pZKo5HgFpV4mjGbfjk965yJiXUQcBpwGfKaVsmn5tRExGBGDc+bMaTfWrqpVST0Ow6xYXMPIbgSYX7c9D9g+0ckR8RNgkaSDWi073dU+MB7pbVYsZa/pndlGYLGkhZJmASuA9fUnSHqulLTTSDoamAXcm6VsnlRcwzArpLyvuDfQqyeKiDFJZwNXAmXgwojYLGlVenwN8EbgTEm7gMeAt6ad4E3L9ir2Tqt9YNyHYVYseV9AqWcJAyAiNgAbGvatqXt8DnBO1rJ5VfvAeHpzs2Ipl0QERATK4R+MvrGzD2qd3k4YZsVSa1XIa7OUE0YfjDlhmBVSrd8yr6O9nTD6wDUMs2Kq3RmZ1/mknDD6YLwPI4dtmGbWvtofiW6SsszGKr6t1qyISu7DsFZVPXDPrJAGyk4Y1qKK+zDMCmm8huE+DMtqfKS3+zDMCsV9GNYyD9wzKyaPw7CWuUnKrJhq/+er1T4H0iYnjD7wXFJmxVQeH7iXz4zhhNEHrmGYFVPJA/esVVX3YZgV0sB4p3efA2mTE0Yf1AbuOWGYFUvtzkg3SVlmtbukfFutWbG409taVvuw1EZ9mlkxjDdJuQ9jcpKWStoiaVjS6ibHz5B0Q/pzjaQj647dIelGSddLGupl3J1Wq466hmFWLKXxPox8VjF6tuKepDJwHnAKMAJslLQ+Im6uO+124A8j4n5Jy4C1wEvrjp8UEff0KuZucae3WTE9MXCvz4G0qZc1jGOA4YjYGhE7gUuA5fUnRMQ1EXF/unkdMK+H8fVM7cPiyQfNisVTg2Q3F9hWtz2S7pvIe4Dv120HcJWkTZJWTlRI0kpJQ5KGRkdHpxRwt3guKbNiKud8HEbPmqSAZt+OTd81SSeRJIwT63afEBHbJT0DuFrSLRHxkyddMGItSVMWg4OD0/JfxQP3zIqpnP6J7iVaJzcCzK/bngdsbzxJ0hHABcDyiLi3tj8itqe/dwDrSJq4csmTD5oVU7mUfOVWnTAmtRFYLGmhpFnACmB9/QmSDgW+DbwjIm6t2z9b0n61x8CrgZt6FnmHeU1vs2LK+2y1PWuSiogxSWcDVwJl4MKI2CxpVXp8DfAXwIHA+Ure2LGIGAQOBtal+waAiyPiil7F3mljnnzQrJBKOW+S6mUfBhGxAdjQsG9N3eP3Au9tUm4rcGTj/ryq1TBKHjZpVih57/T2V1YfVMbX9Pbbb1YkA76t1lo15hqGWSGVct6H4a+sPqi6D8OskDxwz1pW+7C4ScqsWMqefNBaVXGTlFkhuYZhLfPAPbNiyvs4DCeMPvDUIGbF5NtqrWUVd3qbFVItYdSWac6blhNGOk1HuRvBFIVrGGbFVJrpNQxJJUlvk/Q9STuAW4C7JG2W9HlJi7sf5sxSjaAkkGsYZoVShIF7PwQWAR8HnhkR8yPiGcDLSBY5+pykt3cxxhlnrBquXZgVUG3g3kyeS+pVEbGrcWdE3Ad8C/iWpL06HtkMVq2GF08yK6DxTu+cJoxJaxiNyaJZH0azhGITq7iGYVZI47fVug/DsqqEE4ZZEZVKQprBNQzch9FxrmGYFVdZmvl9GJKeHRHV2k73YbSvUg2PwTArqFJJM7dJqq5/Yl3jMUnHNpyzR5KWStoiaVjS6ibHz5B0Q/pzjaQjs5bNk6qbpMwKa6CkmdskJektkj4H7CfpBQ0d3muzPlFa7jxgGbAEOF3SkobTbgf+MCKOAD5Tu37GsrkxVnHCMCuqPDdJZenD+BlwM/A04FzgNkm/kPRd4LEWnusYYDgitkbETuASYHn9CRFxTUTcn25eB8zLWjZPKuHbas2KqpTjGsakfRgR8Rvgq5J+FRE/A5D0dGAhyR1TWc0FttVtjwAv3cP57wG+32pZSSuBlQCHHnpoC+H1TrUaDJSdMMyKaCDHfRiTJgxJisTPavvSDu/7Gs+Z7FJN9jUtI+kkkoRxYqtlI2ItaVPW4ODgtPxXGXOnt1lhlUqa2VODSPqApN3+XJc0S9IrJX0FOCvDdUaA+XXb84DtjSdJOgK4AFgeEfe2UjYvqhHjk5CZWbGUNbMTxlKgAnxD0nZJN0vaCtwGnA78Q0RclOE6G4HFkhZKmgWsANbXn5AmpW8D74iIW1spmyeVaoxPQmZmxVIuiUp18vOmoyx9GL8Hzpf05fT8g4DHIuKBVp4oIsYknQ1cCZSBCyNis6RV6fE1wF8AB6bPBzAWEYMTlW3l+aeTiueSMiusJGHkM2NkGbhXc01EHAfcBSBpX2BxRPxH1gtExAZgQ8O+NXWP3wu8N2vZvPJIb7PiKpdETtdPamkBpb0BJJ0LEBGPAOd3I6iZrhJePMmsqMo5vq22lYQhSc8A3q4nVv55ShdimvEq1aoThllBJQP3Zn6T1MeBfwMuBv5B0q14TfC2eC4ps+IqzeRO75qIuAJ4HoCk44A3k4yVsBZVq1ByqjUrpHIpv2t6Zx64V78vIq4Frt3TOTaxsWqVpw60Urkzs5miXCrN6LmkOjVwz1KVwAP3zAqqnOMFlLL8mbsUeDfJwL2FwAPAPiTjIa4iGbh3fbcCnImq1cBTSZkVUznHU4NkHrhHMphuL9ocuGdPSMZhuBPDrIjKM3zywXOBG9KfzRFxV9ejmuGShNHvKMysH8olsWtXPm+TytIkNQwcC/wx8AJJd/NEAtkI/CQiHu9eiDNPxSvumRVWKceTD2ZpktptNHfaj/FC4AjgT4D/KelPIuLK7oQ481TdJGVWWAMlzdzbahtFxO0kS6muB5D0LOC7JBMDWgZj7vQ2K6xySYzldDKpKf+Zm/ZpXNyBWAqjUvV6GGZFVVJ+axgdaReJiL/vxHWKohpeD8OsqAbK+e3DcEN6H4x5enOzwspzp7cTRh9UvYCSWWHleRxGTxOGpKWStkgalrS6yfHDJF0r6XFJH2k4doekGyVdL2mod1F3XsVNUmaFNaNHeneKpDJwHnAKMAJslLQ+Im6uO+0+4IPAaRNc5qSIuKergfZApeJOb7OiKrtJKpNjgOGI2BoRO4FLgOX1J0TEjojYCOzqYVw9Vwmvh2FWVHmuYfQyYcwFttVtj6T7sgrgKkmbJK2c6CRJKyUNSRoaHR1tM9Tu8preZsVVyvHAvV4mjGbfkK28aydExNHAMuD9kl7e7KSIWBsRgxExOGfOnHbi7DonDLPiGihpRq+H0SkjwPy67XnA9qyFI2J7+nsHsI6kiSuXPJeUWXH5ttpsNgKLJS2UNAtYQTq9yGQkzZa0X+0x8Grgpq5F2kXVahCBb6s1K6hySTN6AaWOiIgxSWeTzDlVBi6MiM2SVqXH10h6JjAE7A9UJX0IWEKyBsc6JV+yA8DF6RrjuVO7/9q31ZoVU56bpHq6sHREbAA2NOxbU/f4bpKmqkYPAUd2N7reqFVFfVutWTG509syq31Q3IdhVkweh2GZ1T4obpIyK6ZySVQDIoe1DCeMHhtvknKnt1kh1VoX8ljLcMLosdqHxE1SZsU0njBcw7DJVNyHYVZotf/71WqfA2mDE0aPuYZhVmy1eeTGcpgxnDB6bDxhuA/DrJBKrmFYVrUPicdhmBVTOf2v7z4Mm1StGurbas2KqVxOvnbdJGWTqg3ccw3DrJhqzdE5zBdOGL1WST8k7sMwK6a0guEmKZtcrRrqu6TMiqlcSr52KxUnDJtErRrqhGFWTK5hWGZPDNzrcyBm1he1aYE8NYhNqjLeJOW33qyIBmpNUk4YNhl3epsV23iTlBPGnklaKmmLpGFJq5scP0zStZIel/SRVsrmxRMLKPU5EDPri1qTVB4XUerZ15akMnAesIxk2dXTJS1pOO0+4IPAF9oomwtPrIfhjGFWRAPl2lxSThh7cgwwHBFbI2IncAmwvP6EiNgRERuBXa2WzQt3epsVmzu9s5kLbKvbHkn3dbvstFL1AkpmhTY+vbmbpPao2Tdk1ncsc1lJKyUNSRoaHR3NHFyvuEnKrNi84l42I8D8uu15wPZOl42ItRExGBGDc+bMaSvQbhpzp7dZoZXdJJXJRmCxpIWSZgErgPU9KDutVL3inlmh5bmGMdCrJ4qIMUlnA1cCZeDCiNgsaVV6fI2kZwJDwP5AVdKHgCUR8VCzsr2KvZO8gJJZsZVyvKZ3zxIGQERsADY07FtT9/hukuamTGXzyEu0mhVbbS0cTz5ok3LCMCu28dtqc1jDcMLosdqHxLfVmhXT+G21OezDcMLosfHbastOGGZFVGuS8khvm5Q7vc2KreSBe5aV1/Q2KzaPw7DMxiq1kd5OGGZFVHaTlGXlGoZZsbnT2zJzH4ZZsZVzPHDPCaPHxjwOw6zQ8jw1iBNGj1WdMMwKzZ3eltn4AkpukjIrpJJrGJbVE2t6O2GYFZGbpCyzSjV8S61ZgQ2409uyqkS4dmFWYLV55HxbrU2qWg33X5gV2BNNUn0OpA1OGD025iYps0Kr/fevVPOXMZwweqxadZOUWZFJolyS+zAmI2mppC2ShiWtbnJckr6UHr9B0tF1x+6QdKOk6yUN9TLuTqpEeAyGWcGVpVw2SfVsiVZJZeA84BRgBNgoaX1E3Fx32jJgcfrzUuDL6e+akyLinh6F3BWVqhOGWdGVS3KT1CSOAYYjYmtE7AQuAZY3nLMc+GokrgMOkPSsHsbYdRV3epsVXpIw+h1F63qZMOYC2+q2R9J9Wc8J4CpJmyStnOhJJK2UNCRpaHR0tANhd1al6mlBzIquJC+gNJlm35KN79iezjkhIo4mabZ6v6SXN3uSiFgbEYMRMThnzpz2o+2SSrXqhGFWcAPlEmNuktqjEWB+3fY8YHvWcyKi9nsHsI6kiSt3KuEahlnRlXLa6d3LhLERWCxpoaRZwApgfcM564Ez07uljgUejIi7JM2WtB+ApNnAq4Gbehh7x1SrgfOFWbGVS/kc6d2zu6QiYkzS2cCVQBm4MCI2S1qVHl8DbABeAwwDjwLvSosfDKxT0lk8AFwcEVf0KvZOGqtWGSh5+ItZkQ2USrlcorVnCQMgIjaQJIX6fWvqHgfw/ibltgJHdj3AHqhUPVOtWdGVSu70tgyqEZT9rpsVWjJwzwnDJjHmcRhmhVcqOWFYBlWP9DYrvAEnDMvCU4OYWUmefNAyqFRjfAEVMyumsmsYlkUlgoGyE4ZZkblJyjJxDcPMSiX5tlqbXNXrYZgVnm+rtUzGKl6i1azoyiXlcqS3E0aPVcNNUmZFVy4pl3NJOWH0mG+rNTOv6W2ZOGGYmW+rtUwq7vQ2Kzx3elsmXtPbzDyXlGXiJikz88A9y8QJw8xK7vSenKSlkrZIGpa0uslxSfpSevwGSUdnLZsX1QgvoGRWcGX5tto9klQGzgOWAUuA0yUtaThtGbA4/VkJfLmFsrng9TDMLK8D93q5ROsxwHC63CqSLgGWAzfXnbMc+Gq6VOt1kg6Q9CxgQYayHfP6f/o3fr+r0o1L89Bju9wkZVZw5ZK4+8Hfc8q5P+7K9Z/21Flcuuq4jl+3lwljLrCtbnsEeGmGc+ZmLAuApJUktRMOPfTQtgJdNGc2OyvVtspO5nkH78epRx3SlWubWT688eh5PLazQtCdWsb+++zVlev2MmE0+7O68d2a6JwsZZOdEWuBtQCDg4Nt/Wt8ccWL2ilmZpbJcYsO5LhFB/Y7jJb1MmGMAPPrtucB2zOeMytDWTMz66Je3iW1EVgsaaGkWcAKYH3DOeuBM9O7pY4FHoyIuzKWNTOzLupZDSMixiSdDVwJlIELI2KzpFXp8TXABuA1wDDwKPCuPZXtVexmZgaKHA4eyWpwcDCGhob6HYaZWW5I2hQRg82OeaS3mZll4oRhZmaZOGGYmVkmThhmZpbJjO70ljQK3Nlm8YOAezoYTj/5tUxPfi3T10x6Pa2+lmdHxJxmB2Z0wpgKSUMT3SmQN34t05Nfy/Q1k15PJ1+Lm6TMzCwTJwwzM8vECWNia/sdQAf5tUxPfi3T10x6PR17Le7DMDOzTFzDMDOzTJwwzMwsEyeMBpKWStoiaVjS6n7HMxWS5kv6oaRfStos6U/7HdNUSCpL+g9J3+13LFOVLj98maRb0n+fzq+n2SOS/mv6+bpJ0jck7dPvmLKSdKGkHZJuqtv3dElXS7ot/f20fsaY1QSv5fPpZ+wGSeskHTCV53DCqCOpDJwHLAOWAKdLWtLfqKZkDPiziHgBcCzw/py/nj8FftnvIDrkH4ErIuIw4Ehy+rokzQU+CAxGxOEkyw+s6G9ULbkIWNqwbzXwg4hYDPwg3c6Di3jya7kaODwijgBuBT4+lSdwwtjdMcBwRGyNiJ3AJcDyPsfUtoi4KyJ+kT5+mORLaW5/o2qPpHnAa4EL+h3LVEnaH3g58L8BImJnRDzQ16CmZgB4iqQB4KnkaDXMiPgJcF/D7uXAV9LHXwFO62VM7Wr2WiLiqogYSzevI1mttG1OGLubC2yr2x4hp1+wjSQtAF4E/LzPobTri8BHgWqf4+iE5wCjwD+nTWwXSJrd76DaERG/Ab4A/Bq4i2SVzKv6G9WUHZyu9En6+xl9jqdT3g18fyoXcMLYnZrsy/19x5L2Bb4FfCgiHup3PK2S9DpgR0Rs6ncsHTIAHA18OSJeBPyO/DR77CZt318OLAQOAWZLent/o7JGkv6cpIn661O5jhPG7kaA+XXb88hR9boZSXuRJIuvR8S3+x1Pm04ATpV0B0kz4Sslfa2/IU3JCDASEbXa3mUkCSSPXgXcHhGjEbEL+DZwfJ9jmqrfSnoWQPp7R5/jmRJJZwGvA86IKQ68c8LY3UZgsaSFkmaRdN6t73NMbZMkknbyX0bEuf2Op10R8fGImBcRC0j+Tf5fROT2r9iIuBvYJun56a6TgZv7GNJU/Bo4VtJT08/byeS0A7/OeuCs9PFZwP/tYyxTImkp8DHg1Ih4dKrXc8Kok3YOnQ1cSfKhvzQiNvc3qik5AXgHyV/k16c/r+l3UAbAB4CvS7oBOAr42/6G0560lnQZ8AvgRpLvlNxMqyHpG8C1wPMljUh6D/A54BRJtwGnpNvT3gSv5X8A+wFXp///10zpOTw1iJmZZeEahpmZZeKEYWZmmThhmJlZJk4YZmaWiROGmZll4oRh1kQ6m+z76rYPkXRZl55rL0kzZRS7zWBOGGbNHQCMJ4yI2B4Rb+rSc50IXNOla5t1jBOGWXOfAxalg50+L2lBbZ0BSe+UdLmk70i6XdLZkj6cTiR4naSnp+ctknSFpE2SfirpsAmeaylNJoWT9Iikc9Ly/yrpGEk/krRV0qnpOf9F0r+ncd4gaXGX3g8zJwyzCawGfhURR0XEf2ty/HDgbSRT4v8N8Gg6keC1wJnpOWuBD0TEi4GPAOdP8FwnAT9qsn828KO0/MPAX5OMPH4D8FfpOauAf4yIo4BBknmqzLpioN8BmOXUD9M1Rh6W9CDwnXT/jcAR6QzBxwPfTKZYAmDvxotIOgS4b4J5fnYCV9Rd9/GI2CXpRmBBuv9a4M/T9UK+HRG3Tf2lmTXnGoZZex6ve1yt266S/CFWAh5Iayi1nxc0uc4ykrnLmtlVN7vo+HNERO05iIiLgVOBx4ArJb1yCq/JbI+cMMyae5hk0ra2pOuO3C7pzZDMHCzpyCanNu2/yErSc4CtEfElkllWj2j3WmaTccIwayIi7gV+JukmSZ9v8zJnAO+R9J/AZhqW+03XkF8cEbdMIdS3AjdJuh44DPjqFK5ltkeerdasTySdCLw9Ilb1OxazLJwwzMwsEzdJmZlZJk4YZmaWiROGmZll4oRhZmaZOGGYmVkmThhmZpbJ/wddlKrStQhNrgAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X_pi_half_time = np.pi / 2 / delta_Bz\n",
    "n_steps_X_pi_half = int(X_pi_half_time / time_step)\n",
    "\n",
    "# Control parameter during X_pi/2 pulse\n",
    "eps_idle = -20  \n",
    "\n",
    "# Control param. resulting in precession freq. J(eps_work) = pi / 8\n",
    "eps_work = eps_0 * np.log(np.pi / 8 / J_0)  \n",
    "\n",
    "def fid_pulse(tau):\n",
    "    \"\"\"\n",
    "    Generates the control pulse (epsilon) for the FID sequence.\n",
    "    \"\"\"\n",
    "    n_steps_tau = int(tau / time_step)\n",
    "    n_steps_pulse = 2 * n_steps_X_pi_half + n_steps_tau\n",
    "    \n",
    "    ctrl_parameter = np.zeros((n_steps_pulse,1)) + eps_idle\n",
    "    \n",
    "    ind_start, ind_end = n_steps_X_pi_half, n_steps_X_pi_half + n_steps_tau\n",
    "    ctrl_parameter[ind_start:ind_end] = eps_work\n",
    "\n",
    "    return ctrl_parameter\n",
    "\n",
    "\n",
    "def create_monte_carlo_solver(tau, pulse_function):\n",
    "    \n",
    "    control_pulse = pulse_function(tau)\n",
    "    n_time_steps = len(control_pulse)\n",
    "    \n",
    "    active_dBz = np.where(control_pulse == eps_idle, 1, 0).flatten()\n",
    "\n",
    "    h_drift = [0.5 * delta_Bz * active * pauli_x  for active in active_dBz]\n",
    "\n",
    "    ntg_pink = NTGColoredNoise(\n",
    "        n_samples_per_trace=n_time_steps,\n",
    "        noise_spectral_density=pink_spectrum,\n",
    "        dt=time_step,\n",
    "        n_traces=n_traces,\n",
    "        low_frequency_extension_ratio=500\n",
    "    )\n",
    "    \n",
    "    solver = SchroedingerSMCControlNoise(\n",
    "        h_drift=h_drift,\n",
    "        h_ctrl=[.5 * pauli_z,],\n",
    "        tau=time_step * np.ones(n_time_steps),\n",
    "        noise_trace_generator=ntg_pink,\n",
    "        amplitude_function=exchange_coupling_ampl_fctn\n",
    "    )\n",
    "    \n",
    "    solver.set_optimization_parameters(control_pulse)\n",
    "    \n",
    "    return solver\n",
    "\n",
    "sample_solver = create_monte_carlo_solver(10, fid_pulse)\n",
    "sample_amplitude = sample_solver._ctrl_amps\n",
    "\n",
    "plt.plot(np.arange(len(sample_amplitude)) * time_step, sample_amplitude)\n",
    "plt.xlabel(\"time / ms\")\n",
    "plt.ylabel(\"$J(\\\\varepsilon(t))$\")\n",
    "plt.title(\"FID Control Amplitude\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "precession_times = np.arange(0,100,3)\n",
    "\n",
    "initial_state = DenseOperator(np.array([[1], [0]]))\n",
    "measurement_state = DenseOperator(np.array([[0], [1]]))\n",
    "\n",
    "probabilities_fid = np.zeros(len(precession_times))\n",
    "\n",
    "for i, tau in enumerate(precession_times):\n",
    "    \n",
    "    solver = create_monte_carlo_solver(tau, fid_pulse)\n",
    "    probabilities_fid[i] = average_monte_carlo_traces(\n",
    "        solver, initial_state, measurement_state)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [],
   "source": [
    "# np.save(r\"..\\stored\\precession_times_1\", precession_times,)\n",
    "# np.save(r\"..\\stored\\probabilities_fid_1\", probabilities_fid,)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [],
   "source": [
    "precession_times_loaded = np.load(r\"..\\stored\\precession_times_1.npy\")\n",
    "probabilities_fid_loaded = np.load(r\"..\\stored\\probabilities_fid_1.npy\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit parameters:  [46.70488819  1.96333833]\n"
     ]
    }
   ],
   "source": [
    "def T2_star_decay(t, T2, alpha):\n",
    "    omega = exchange_coupling(eps_work)\n",
    "    return 0.5 * (1 + np.exp(-(t / T2)**alpha) * np.cos(omega * t))\n",
    "\n",
    "popt, cov = curve_fit(T2_star_decay, precession_times_loaded,\n",
    "                      probabilities_fid_loaded, p0=[50, 2])\n",
    "print(\"Fit parameters: \", popt)\n",
    "\n",
    "precession_times_fine = np.linspace(min(precession_times_loaded),\n",
    "                                    max(precession_times_loaded), 500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEJCAYAAAB7UTvrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABJiUlEQVR4nO3deXxU5b348c8zWSArIQuBJEBYwxYgLG5oBbEC6m3R2op1qbbW2mpd2vIT2972tt4rXLXXrYtVr7W3WrFVRBSVior7whJ2CCBrEpYQyEISyDLP749nzjAJWSYzZ85Mku/79eJ1yJlzTp6TwPmeZ/s+SmuNEEIIAeAKdwGEEEJEDgkKQgghvCQoCCGE8JKgIIQQwkuCghBCCK/ocBcgWOnp6To3NzfcxRBCiC5l7dq1R7XWGS33d/mgkJuby5o1a8JdDCGE6FKUUvta2y/NR0IIIbwkKAghhPCSoCCEEMKry/cpCCHs1dDQQHFxMSdPngx3UYQNevfuTU5ODjExMX4dL0FBCNFMcXExSUlJ5ObmopQKd3FEELTWlJeXU1xczJAhQ/w6x7HmI6XUM0qpI0qpzW18rpRSjymldimlNiqlJoWqLEsLS5i26F2GLFjOtEXvsrSwJFTfSogu5+TJk6SlpUlA6AaUUqSlpXWq1udkn8KzwOx2Pp8DjPD8uQX4UygKsbSwhHuXbKKkog4NlFTUce+STRIYhPAhAaH76Ozv0rGgoLX+ADjWziFfB/5PG58BKUqpAXaX48EVRdQ1NHGF60O+7voIhZu6hiYeXFFk97cSQoguJ5JGH2UDB3y+LvbsO4NS6hal1Bql1JqysrJOfZPSijoAroz6kEdj/8iPo5Y22y+ECL//+q//YuzYsYwfP56JEyfy+eefA3DzzTezdetWW75Hbm4uR48ebfeY+++/v9nX5513ni3f+8Ybb2TIkCFMmDCBkSNHcsMNN1BS0nFrxSOPPEJtba0tZWhLJAWF1uo4ra4ApLV+Ums9RWs9JSPjjFna7cpKiQPghoYFvN50NrdFLyWNSu9+IUR4ffrpp7z++uusW7eOjRs3snLlSgYOHAjA008/zZgxYxwrS8ug8Mknn9h27QcffJANGzZQVFREQUEBM2bMoL6+vt1zelpQKAYG+nydA5Ta/U3mz8ojLiYKjYtHGr9BL9XIvNiPmD8rz+5vJUSPYPfAjYMHD5Kenk6vXr0ASE9PJysrC4Dp06d709okJiZyzz33MHnyZC6++GK++OILpk+fztChQ1m2bBkAzz77LLfffrv32pdffjmrVq0643vOnTuXyZMnM3bsWJ588kkAFixYQF1dHRMnTuTaa6/1fk8wo3rmz5/PuHHjyM/P58UXXwRg1apVTJ8+nauuuopRo0Zx7bXX0tHqlkop7r77bvr378+bb74JwA9/+EOmTJnC2LFj+fWvfw3AY489RmlpKTNmzGDGjBltHhesSAoKy4AbPKOQzgEqtdYH7f4mcwuyWXhlPtkpcXypc9imhvPdtC3MLWi1pUoI0Y5QDNy45JJLOHDgACNHjuRHP/oR77//fqvH1dTUMH36dNauXUtSUhK//OUvefvtt3nllVf41a9+1anv+cwzz7B27VrWrFnDY489Rnl5OYsWLSIuLo7169fz/PPPNzt+yZIlrF+/ng0bNrBy5Urmz5/PwYPmcVVYWMgjjzzC1q1b2b17Nx9//LFfZZg0aRLbt28HTPPZmjVr2LhxI++//z4bN27kjjvuICsri/fee4/33nuvzeOC5eSQ1BeAT4E8pVSxUup7SqlblVK3eg55A9gN7AKeAn4UqrLMLcjm4wUXsWfRZYz+ylWkHd8Ate31gQshWmMN3PAV7MCNxMRE1q5dy5NPPklGRgZXX301zz777BnHxcbGMnu2GdCYn5/PhRdeSExMDPn5+ezdu7dT3/Oxxx5jwoQJnHPOORw4cICdO3e2e/xHH33ENddcQ1RUFJmZmVx44YWsXr0agLPOOoucnBxcLhcTJ070uyy+NYp//OMfTJo0iYKCArZs2dJmP4q/x3WGY5PXtNbXdPC5Bm5zqDinDZsB7y+CfZ/A6Msd//ZCdGVtDdAIduBGVFQU06dPZ/r06eTn5/PXv/6VG2+8sdkxMTEx3uGWLpfL29zkcrlobGwEIDo6Grfb7T2ntfH6q1atYuXKlXz66afEx8czffr0Dsf1t9ckZJXDug+rLB0pLCxk5syZ7Nmzh4ceeojVq1fTt29fbrzxxlbL4+9xnRVJzUfhMWAiRMXCgc/DXRIhupy2BmgEM3CjqKio2Zv6+vXrGTx4cEDXys3NZf369bjdbg4cOMAXX3xxxjGVlZX07duX+Ph4tm/fzmeffeb9LCYmhoaGhjPO+cpXvsKLL75IU1MTZWVlfPDBB5x11lkBlVFrzWOPPcbBgweZPXs2VVVVJCQk0KdPHw4fPuztZwBISkqiuroaoN3jgiFpLmJ6m8Bw4Mx/LEKI9s2flce9SzY1a0KKi4kKauDGiRMn+PGPf0xFRQXR0dEMHz7c2/nbWdOmTWPIkCHk5+czbtw4Jk06M1HC7NmzeeKJJxg/fjx5eXmcc8453s9uueUWxo8fz6RJk5r1K1xxxRV8+umnTJgwAaUUDzzwAP379/f2Cfhj/vz53HfffdTW1nLOOefw3nvvERsby4QJEygoKGDs2LEMHTqUadOmNSvPnDlzGDBgAO+9916bxwVDddQzHummTJmig15kZ/nPYMNiuPcAyExO0cNt27aN0aNH+3380sISHlxRRGlFHVkpccyflScDNyJMa79TpdRarfWUlsdKTQEgcwzUV0PlAUgZFO7SCNGlzC3IliDQjUifAkA/z2SYI9vCWw4hhAgzCQoAGaPM9og90+cdsft9+N0o+P1ZUP5luEsjhOgmJCgAxKVAcnbXqSnUVcDL3wPlgupD5u9dvG9ICBEZJChY+o3uOjWF9c9DTRlc/TeYfT+UFsKud8JdKiFENyBBwdJvNJTtgCb/JpqEg5VjZvubf2KzGsnSI/0h/1sQlwobXgh38YQQ3YAEBUvaCGg6BVWRudiOlWPGVbmXUa4DvFJ/lskxs6kMxnwNit6EBllTV3QPUVFRTJw40ftn79693rTVe/fu5e9//3uYS9h9SVCw9M012+N7w1mKNlk5Zi50mYRXK92TTueYGTkbGmqgWCbgie7BSkRn/cnNzfWmrZagEFoSFCwRHhSsXDKTXDs5olPYpzNP7x88DVQU7F4VxhIKEVpW2uoFCxbw4YcfMnHiRB5++OEwl6r7kaBgSc4GV3TEBgUrl0yB2kmhezjWmkRZKXHQOxkGjIfi1WEsoRD2sdYxmDhxIldccUWzzxYtWsQFF1zA+vXrufvuuyktLeXSSy8NU0m7H5nRbImKNrOZj+8Jd0laNX9WHg8s+ZghrsMsbrgIaJFjJqsANr0Mbje4JNaLrs1qPvJHVlYWb7zxRmgL1IPI08NX39yIrSnMLcjmd+eZbI2F7hFkp8Sx8Mr80+kFsgrgVGXEBjUhRNcgNQVffXPNmP8IdW7sbnBF84/f/BBi45t/mFVgtqWFkDbM+cIJ4RDf9NHCflJT8NU3F+qOmxnDkejwZkgfeWZAAJOqI7p3RAc1Iewwfvx4oqOjmTBhAg8//LD0KdhMagq+fEcgxU0MY0HaUFYEWRNb/ywqBvrnS1AQ3cKJEyfa3BcTE8M77zSfwS99CvaRmoKvPgPNtrI4vOVoTUOdCVbp7SxekjkODm+RPEhCiIBJUPDVJ8dsq0rDW47WlO8CNGS0ExTSR8LJCqgtd6pUQohuRoKCr/h0cMVEZqqLsiKz7SgoABzdEfryCCG6JQkKvlwuSM6K3KCgXJA2vO1j0keYrQQFIUSAJCi0lJwdmc1HZduh7xCI7tX2MX0GQnQcHN3pXLmEEN2KBIWWIrWmcHTn6eahtrhckD5cagqi29u6dSvPPvssBw4ckDkLNpOg0FIfT03B7Q53SU7TGir2QerQjo9NHylBQXR5Tz31lDf3kcvl8v79Jz/5CQANDQ08/vjjvPLKK95EeXZoamqioKCAyy+/3LuvoqKCq666ilGjRjF69Gg+/fTTZucUFRU1S/OdnJzMI488AsBbb71FXl4ew4cPZ9GiRbaVM5RknkJLydnQVG9G8CRmhLs0Rk0ZNNRC38EdH5s+EjYvMWsrxPQOfdmECIHvf//7fP/736ekpITzzjvvjDxIBw4c4KabbmLo0KFUV1eTnJxsy/d99NFHGT16NFVVVd59d955J7Nnz+all16ivr6e2traZufk5eV5y9fU1ER2djZXXHEFTU1N3Hbbbbz99tvk5OQwdepUvva1rzFmzBhbyhoqUlNoKTnLbCOpCcnKx5TiR1DoOwTQULE/lCUSwhGbN28mPz//jP2XX345V111FZdeeqltAaG4uJjly5dz8803e/dVVVXxwQcf8L3vfQ+A2NhYUlJS2rzGO++8w7Bhwxg8eDBffPEFw4cPZ+jQocTGxjJv3jxeffXVM87Zu3cvo0aN4uabb2bcuHFce+21rFy5kmnTpjFixAi++OILampquOyyy5gwYQLjxo3jxRdftOWeWyM1hZaSPQnmqkranj3stOP7zNaacd0e31nZGR30QQgR4TZt2sS4ceNa/ax///7tnnvBBRe02t/w0EMPcfHFF5+x/6677uKBBx5ods7u3bvJyMjgpptuYsOGDUyePJlHH32UhISEVr/n4sWLueaaawAoKSlh4MCB3s9ycnL4/PPPWz1v165d/POf/+TJJ59k6tSp/P3vf+ejjz5i2bJl3H///Vx//fVkZWWxfPlyACorK9u992BIUGjJCgqVEVRTqNhrtimDOj7WamKq2Bey4oge5M0FcGiTvdfsnw9z/Gtf37x5M1/96lcD+jYffvih38e+/vrr9OvXj8mTJ7Nq1Srv/sbGRtatW8fjjz/O2WefzZ133smiRYu47777zrhGfX09y5YtY+HChQDoVjILKKVa/f5Dhgzx1ojGjh3LzJkzUUqRn5/P3r17yc/P52c/+xn33HMPl19+ORdccIHf99ZZ0nzUUkKGWcXsxKFwl+S043shoV/rifBaSsw0ifEiNAW4EJ3RXk2hIxdccEGzDmDrz8qVK8849uOPP2bZsmXk5uYyb9483n33Xa677jpycnLIycnh7LPPBuCqq65i3bp1rX6/N998k0mTJpGZaVZFzMnJ4cCBA97Pi4uLycrKavXcXr1ODzV3uVzer10uF42NjYwcOZK1a9eSn5/Pvffey29/+9uAfib+kJpCSy4XJPaDE4fDXZLTju/zr+kIQCnPYkF7Q1ki0VP4+UYfCm63m507dzJq1CjvvqVLl7J8+XKOHDnCbbfdxiWXXNLm+Z2pKSxcuND7hr9q1SoeeughnnvuOQAGDhxIUVEReXl5vPPOO212FL/wwgvepiOAqVOnsnPnTvbs2UN2djaLFy8OeG3p0tJSUlNTue6660hMTOTZZ58N6Dr+cCwoKKVmA48CUcDTWutFLT7vAzwHDPKU6yGt9V+cKl8zif3gxJGwfOtWHd8Hg872//i+uV2j+ejoLvj8CdO5f+5t7U/MEz3Orl27yMnJafYWPXfuXObOncvx48f52c9+1m5QsMvjjz/OtddeS319PUOHDuUvfzGPpUsvvZSnn36arKwsamtrefvtt/nzn//sPS86Oprf//73zJo1i6amJr773e8yduzYgMqwadMm5s+fj8vlIiYmhj/96U+23FtrVGvtXrZ/E6WigB3AV4FiYDVwjdZ6q88xPwf6aK3vUUplAEVAf611fXvXnjJlil6zZo29BX7+W1B9EG71/00jZJoa4D/7wfk/gZn/7t85y38GG1+EBftNzSESVRyAJ6eb9St0E4y6HK5+LnLL24Ns27aN0aNHh7sY7frpT3/Ktddey6RJk8JdlC6htd+pUmqt1npKy2Od6lM4C9iltd7tecgvBr7e4hgNJCnTE5MIHAMaHSpfc5FUU6g+BNoNKQM7PhZYWljC44UNcKqK2YuWsbQwgjrMff3rl2buxW2fw1fvg+2vw9Yzh+sJ4UtrzT333MOcOXMkIISIU0EhGzjg83WxZ5+v3wOjgVJgE3Cn1rrVacVKqVuUUmuUUmvKysrsL21ippkw5m6y/9qdZc2XSG754zrT0sIS7l2yic21KQBEV+3n3iWbIi8wlH8JW5eaJqP0EZ7tSPjgIVkLQrTr8ccfZ+XKlbz00ks88cQT4S5Ot+RUUGitTaDl//5ZwHogC5gI/F4p1eqsFK31k1rrKVrrKRkZIZh1nJhpmjRqj9l/7c7qRFB4cEURdQ1NFGvzM8lRR6lraOLBFUWhLGHnrfs/cEXDVM8kIVcUnPMjOLwJStaGt2wiot1xxx2sXbuWJ554gltvvTXcxemWnAoKxYBv+0cOpkbg6yZgiTZ2AXuAUYRDYj+zjYQRSNZ8ieTWh7L5Kq2oA+CgTgNggCpvtj9ibF8OQ74CSacnH73unsYpYnj6iYeYtujdyKvdCNFDOBUUVgMjlFJDlFKxwDxgWYtj9gMzAZRSmUAesNuh8jWXaMYZR0RQqCqF2ETo3afDQ7NS4gA4RhKndAz91bFm+8NlaWEJ0xa9y5AFy7n6/uegfCeMnN3s8/mv7ebTpjHMdK2jpKI2Mpu9hOgBHAkKWutG4HZgBbAN+IfWeotS6lallFUHvA84Tym1CXgHuEdrfdSJ8p0hKZKCQrGpJfgxKmf+rDziYqIAxUGdygB1jLiYKObPame1thCz+jlKKurQwNiazwD4V8ME7zFWs9dK9ySGuA4zTJVGZrOXED2AY/MUtNZvAG+02PeEz99LgdAPOvZHQgQ1H1WV+tWfADC3wBz34IoiDtWmMjj6OAu/lu/dHw7WA98yw1XIDnc2v/molkvON/us5q13mwog5i/MdK3jy6bsyGv2EqIHkDQXremVaJpsImFYaieCApjA8PGCizhn4ngmJNeGNSBA8/6MKJqY5NrJJ+6xzfZbzVulpLPdPZDzXZub7RdCOEeCQlsiIdVFU4OZp9AngAd7chZUl4Z9WK3vg32kKiZBnWKde0Sz/aebvWCNeyQFrl0kxKiwNnsJ0VNJUGhLYmb4awrVhwDt18ijMyRng7vRzLcII98HfoFrFwBbo0Y1e+DPLchm4ZX5ZKfEsdadR5Kq4/GZvcJeyxGiJ5KEeG1J7AeHt3Z8XCh1Yo7CGXxTgCe1n3c+lHz7OSbV7OQYfbj9iovOeODPLcg2+44PhUf/yEXxu/EMRhNCOEhqCm2JmJoCkDSg8+daTU4RsIKc1c9xVb+DpOZNY+6knLYPThls7vdA64uRCCFCS4JCWxIz4VQlNIRxBIzVpxHIm35y5AQFAE5WmfkJ2ZPbP04pGHgWHPjCmXIJIZqRoNAW7wS2MNYWqg+ZdBBxqZ0/Nz4NonpFTlA44mmK63/mertnGDDRpP6uOx7SIgkhziRBoS2REBROHDZzJlwB/JqUMh3UVS2ziYTJ4S1m26/1BUqaGeCZ2Gb3MpBCiA5JUGiLN/9RGJflrD50enZ1IPrkRM5a00e2Qq8+pkwdsYLCwY2hLZMQ4gwSFNqSkG62NeHJtAGYWkpiECOHIq2mkDnGv0V0EtIhKQsObgh9uYQQzUhQaEu8JyjUhjMoBFlTSM6OiAlsaG2G9/rTdGQZMB4OSU1BCKdJUGhLTG+ITYKa8vB8/6ZGU0sJtqbgbgz/0NrKYjOSK7MTQaH/eDi6A+prQ1cu0SYnlukVzujs71KCQnsS0sJXU6g5AujTfRuBsOY3hLNfBKBsu9l2pqaQOdYsQ1q+MzRlEm3q3bs35eXlEhi6Aa015eXl9O7d2+9zZEZze+LTw5cmwjtxLYiagnVudZhzOJWb9BakjfD/nAzP+kplRac7noUjcnJyKC4uJiRL3QrH9e7dm5wcPwZ4eEhQaMPSwhLSD0Jq4x6+v+hd5s/KczYXj9XkE0zzkXdYbZhrCkd3mkWCrM57f6QONXM0jmwLXblEq2JiYhgyZEi4iyHCRJqPWmEtDFNan0CqqqKkos75lcCsB3kwHc1WUAh7TWGnqSX4M/LIEh0LqcNMTUEI4RgJCq2wFoYpJ5lUqgDt/Epg1oM8IYg+hehYMxs63DWF8i8hbXjnz8vIO90fIYRwhASFVlgLwJTrZGJVE0nUNdvviBOHTKqK6NjgrpM0ILw1hfoak2ojPZCgMAqO74GGk/aXSwjRKulTaEVWShwlFXUc00kApKlKqnW8syuBVR8+3fwTjKTM8NYUyr802850Mlv6jfKMQNoF/cfZWy67NNTx4ap/8fjqE6yuSiErJc75/ichbCQ1hVZYC8McIxmAVKqJi4lydiWwEzYFhcT+p0cyhYM1pDSg5iNrBFKENiHt/4y6h/K54OMb+Ef9j7g/+ikOV1Q73/8khI2kptAK6y1v2ZsHoB5GJJ7khkvzHR59dBjSA3i7bikp01zL7Q4ssV6wynebberQzp+bNhyUKzI7m8u/hOe/SVl9Av9ZfzeTXTv4QfRyNPDzhu/z4IoiqS2ILkmCQhvmFmQzd8hseORn/PfsLHDyP7jW9tYU3I1Qd6xzQ0LtcnyvKUNsfOfPje5lgklZhA1L1RpeuxOU4pq6eykhnX+5p9JIFLdFL+Nf7qm8XzEx3KUUIiDSfNSehDDlP6o7Dk319vUpQPiakCr2Qd/cwM/PGAVHIqz5aM8HsPdDmPELSBno3f1I41XscWfy8+jnGdgnyAECQoSJBIX2xMRBTILzmVKtWdTBpLiwhDvVxfG90Hdw4OenjzAjkJoabStS0FYtMj/XSd/x9j8BNBDNosZrGOkq4ZHRERbIhPCTBIWOJKSFLyjY0dwTzglsjfVmOGowNYXUYab5q3K/bcUKyqHNsP8TOO/HENObuQXZLLwyn+yUOBSwOekrVCbnMengi6aZSYguRvoUOpKQ4XzzkRWEEjKCv5Y3/9HB4K/VWZUHzJDSlCBqCmnDzLb8y8A6q+224QVwxcD4ed5dcwuym3cqr/2R6XPY9wnkTgtDIYUInNQUOhKfHsaagg1BISbOrHh2Igw1hYp9ZhtMTcEaymrNdwinpkbY+A8YOcvUINuS/y3onQJr/texoglhFwkKHUkIR1A4CiiTosIOSZnh6Wg+vtdsg+lTSMgw61oci4CgsPcDk9J8wrz2j4uNh3HfgO1vwKlqZ8omhE0kKHQk3rOmgpPtwzVlEJ8KUTa17iVmhqemcHwfRMWe7uwOhFKmCclKvx1ORW9BdG8YNrPjY8d/CxrrYPvy0JdLCBtJUOhIQoYZHurkG19NmT1NR5akAeGrKfQZCK6o4K6TNiz8zUdaw443Yeh0/+Zc5JwFfQaZ5iYhuhAJCh0Jx1yFmqM2BwVPTcHp0TDBzlGwpA4zndaNp4K/VqDKtkPFfhg527/jXS4YdwXseR9OVoa2bELYyLGgoJSarZQqUkrtUkotaOOY6Uqp9UqpLUqp950qW7viPUHBybWaa4+aZiu7JPaHxpPOP5yCnaNgSRtuRjFZfRThsGOF2Y6c5f85eZea4bS7VoamTEKEgCNBQSkVBfwBmAOMAa5RSo1pcUwK8Efga1rrscA3nShbh6xRJk4uy2l785E1LNXBJqSTlWZmth01Bd9hqeGy5wNIz4PkLP/PyZlqgnvRm6ErlxA2c6qmcBawS2u9W2tdDywGvt7imG8DS7TW+wG01kccKlv74h1uPmpqMA9TO4NCOJblPO4ZjhrMHAWLNT8hXCOQmhpg/2cw5ILOneeKMs1NO/9lriFEF+BUUMgGDvh8XezZ52sk0FcptUoptVYpdYNDZWuf1afg1LDU2vLm39cO3pqCgyOQKjwzkFMGBX+t+FQzPDdcI5BK1kFDDQz5SufPzZtjak37P7O/XEKEgFNBobXFeVv2ekYDk4HLgFnAvyulRrZ6MaVuUUqtUUqtKSsLcbNObALExJ9+WIeanRPXLFZQcLKmUFlstnYEBQjvCKS9H5jt4PM7f+7QGWZYrjQhiS7CqaBQDAz0+ToHKG3lmLe01jVa66PAB8CE1i6mtX5Saz1Faz0lI8PGh2dbnJzVHIqg0CvJJPZzsqZQVWzG9NvVYZ42PHxBYd8n0G9M+7OY29IrEXIvkM5m0WU4FRRWAyOUUkOUUrHAPGBZi2NeBS5QSkUrpeKBs4HISKSfkOZcR7M375HNax8k9nO4plBiOmVVa5XEAKQOg+pSqK+153r+cruheC0MPCvwawybAUeLzM9EiAjnSEI8rXWjUup2YAUQBTyjtd6ilLrV8/kTWuttSqm3gI2AG3haa73ZifJ1KD696weFpP5wwsG++8piSLZxYSJrBNKx3Y6t17y0sIR/vrmS5+sr+a8NSYzNKQlsNbWhM8x293tQcJ29hQzCq2v38eWbj3NR/bsMdx3EFZdCfP7lcN4dzdaJED2LY/MUtNZvaK1Haq2Haa3/y7PvCa31Ez7HPKi1HqO1Hqe1fsSpsnUoId3ZPgVXtEmoZqdEh/MfVZWY2cx28Q5LdaazeWlhCfcu2URWjXkvebdmUOBrL2eOhYR+8OV7NpcycG98tpEhy67gJ41PodC81Hg+n9Rm07TmWfjjObBhcbiLKMLEr5qCUqov8G/AFZhRQnswzT3LtNZhSKrjsPg0Z4NCQoZ9zS6WxEznHkpNjSZVdx8bawoOD0t9cEURdQ1NFETvpEInsFsPQDc0Bbb2slKmCWnXO+FbK9tXzVHGrJhHJmXcXv9jXnefAyhohMnJ1byc9Ry88gMzC37aneEtq3Bch/86lVJLgCVAP+Aez8SyHwG9gb8ppVaFtISRID4NGmqdac+uORqatZSTMuFUJTTU2X/tlk4cMjOQ7Ww+6pVkAtux3fZdsx2lFebnNMm1i/Xu4WjPfxVrf6cNnWHmuhzeZFcRA+Nugpe+ywD3Eb5Tfw+vu8/Fd3DguqokuG4JjL0S3v4VrPtb+MoqwsKfmsJ3tdYVvjs8E8weBx73zETu3qwRNLXlgS1A3xl2z2a2eCewHbZnlnF7rOGofXLsvW7qMCh3JihkpcRRWVHOSFXMG01nN9sfiLdqRzEbWPT7P/Fa0reYPysvsP6JYH30MOx5n9/F3MYXp0af8XFWShxEx8KVT0Ldcdyv3cWtb1XzdtVgslLiwldu4ZgOawotA4JSKsGTtqLVz7slJ5Pi1ZSdnkVtp0QHJ7CFLCgMdaz5aP6sPKbG7MWlNIXaLPQTFxPF/Fl5nb7W0sIS7n7zMEXuHKa5NlFSURd4/0Qwyr+E9x+A0V9jzGW3edeWtjS7v6gYlo9aSKm7L7869TuSqAlfuYWj/Gk+cimlvq2UWq6UOgJsBw56ktY9qJQaEfpihpmTSfFqy0NUU+hntk4MS63yPDTsbD4CSBtqajoOpDGfW5DN/LFVAGxwDyM7JY6FV+YH9JZs9U985M7nLFcRvainztM/4ai3FkB0L5jzwBlrS7d2f/e/e5Af199OJse5P+ZpgPCUO0SWFpYwbdG7DFmwnGmL3pVg5+FP89F7wErgXmCz1toNoJRKBWYAi5RSr2itnwtdMcPMt/kolOprof5EiPoUrFnNDgxLrSyGXsnQO9ne66b6DEsd0Oq8RluNUXsgdSgb7/hWUNex+iE+dOfzveg3meoq4iN3fuD9E4HY+5HJwXTxbyDZLHp0xtrSLZRW1FHCCB5uvIr/F/MirzWtZoV7qrPlDhFrdFldQxOAtxYE9PjmMX+GQVystb5Pa73RCggAWutjWuuXtdbfAF4MXREjgDWTNdTNR9b1Q1FTiE8DFeXMsNTKEvtrCeB8ttSDG6H/+KAvY/VDfOEeRb2O4nzX5mb7Q05reOe3kJQFZ//A79Os8j3ZdBlb3YP5bcxfSKLWuXK3wY43fKv25qs71YKC4U+fQofpHf05pkvr1cc8UENdUwhFiguLK8pc14llOauK7e9PAGeHpdZVmEWCBgQfFObPyiMuJopaelOoRzDNtSng/omA7PsEDnwOF/wEYvx/oFvlbiSaBQ03k04lv4h90blyt8J6wy+pqENDwP0cbdV2ukMtKFidntGslPo28DWgCTOW7TWt9Qt2FyyiuFzmTTvU+Y9qQlhTgNMrsIVaZQlkFdh/3dgEs7SoEyOQDnmGjvYPvpnKao54cEURH1eP466Yl/nd7BwudaqZ4uNHTb9YJ2dT+5Z7U8UwXoq+nHlNr0Hqfs5McuwM3zf8VKqYHbWa89jCiGVH4ONocMVA6hBTwxvxVRgwsdV5IVkpcZS0EgDCXQvyx9LCEh5cUURpRV1IRoQFkubiQq31POsLpdQfgO4dFCDkE9iWFpawZfkH/AL4xv/t5Po5/e1v23RiVnNDnWkGSw5BTQE8I5CcCAobzdaGmgL4tN8fSIT/fYlLk3YCDqTrOLwVdq6AGb/oVC3B0qzf4dSF8Mf18Prd8IMPzNBVh5VW1JGjjnBX9BLmuj4iWrkp1unsaMwhr1+uWbL16E709uWoVfez153Ja7FzGHLJD7n8rFHe68yfldesTwECH13mJCf6QgIJCr2UUpdh1kfIASI/tNohhKkurF/0De5yiIGtVTGh6fRKzDTt5KFU5Ul+a+dsZl+pQ2HHW6G5tq+DG80wXmvUll2yJkFsEuxeBWOvsPfaHr5vkr9PeJpZUXFET705+Av3SoQ5D8Dia+DT35vmKCc1NbIgcTk3NSzGjYu/Ns3ipaavsE0PIjslnou+dRFg7v+BJR9znnsN34paxY8bn6Vi+T/ZsucWxl4xH2ITmtWCQvXG3ZI/b/gdHdNeX4hjQUEp9ZDW+mc+u34EXAnkYwLD7baUJNLFp8KR0CRttX7RadFV1Ope1NEbbP5FA2YEUs0RM6vVFdXx8YEI1RwFS9ow0/dyssr+0U2+Dm20rZbQTFS0WcFt9yr7r03zN8lkTjCz8UNe0hfSu6iOuXa06I26FEZdbuY7jLsy9BMhLVUH4Z/f4QeNn/OmPpdf11/HEfoCZ77hP7iiiNKGBF7iQl5qupBxajd3R7/MzG0Pw/7nYfoCmPSdDkdf2cmfN3x/jnGiL8SfmsJFvl9orWuB7jv8tC0hXFPB+oWmqUrKdfIZ+22TmGnST9SW2/8GbAnVHAWLd1jql6HptwDTBFZWBHmXhub6Q6dD0RtwbI9p/7aR75vk3KiP6a0a+Fv9RVTY+YIx57/hD2fD8p/Btf+0P09XS6WF8MI15kXgyqc51XQeMSuKUG28Tbf8f7NZD+V7DfOZ1LiDJYNWwPKfwOdPmOG5eXNsKX9n3vCTqCXftZth7lJOvv4c7HTByUpG7CnhBdVEQ2w0jTqaShIo030oW54B6kJIG05enya2V575QmdnX4gjqbO7hfg0s3ZyCN6yrU6vdKooJ7nZfltZqS6qD4UuKFg1hVAFBd9hqaEKCke2gm4KTU0BTqfS3vO+7UHh9ANRc03Uu2x0D2GLzkXZ+YLRJwdm/BxW/By2Lg1ZMxgAez6Ev19taurf+xf0H8dc2m9WbasT+XCfCXDTXSYgv/1r0ww2eBp89T7ImRxwETt8w6+vYWzVh9wWXchU1w6GqVJcyiw8WdcYC4cHQ1wKRxvj0CiiaSRWNTKUUs52baNv4wlYarpt3wKO9UriSz2A3e4s/qfxKqpiMmztC/EnKExQSu0BNgGbfbbbu/1QVF8J6YA2gcHmyWVWp1eaquKQbr1KbAvf/EehUllsRk/F9A7N9ft6HqKh7Gy2+l1smKPQqvQRZs7A7lUw+UZbL209EAvULka7DnBvw/e8+2111g9gwwvw5gJerxnNwndL7W+b//I9U0NIGQTfWXZ6AmYH2u1EVgpGXQYjLoF1/werFsLTF5kEgDN/FVCQbq2dP7XhELuWvwNb9sCeD3gy9hRVOo7V7lEsazqX9Xo4Re6BxPQZwMc/ngnAzxe922owy+0TzarvDzVp48t3UVm0gaj9W5mp1vF/CTdy7+zAZtq3xZ+gsBGzbvI4TD/CJcBPMSupHdBaO7PiSbhZs5pDkMXU+oX2e7WaLY25ZIeq0yvJgaBQFaKJax5LtxznPNL4aOX7/O7TKaH5OR3aaOamhKq9XCnThLTjLdtTaVsPxHn6PWp0L5Y1nReaF4yoaLj8UfTTMzn++q8pqb8BsHE0zM6VsPjbZhnWG16FRP+HafvViRwVA1O/B+O/BZ88bv5sWwajvwZTvgu55/vdrFRaUUcUTUxSO7koqpDprvWMdh2ARuDYUJh6Mx+5JvPDD3tR3XD6mnExUSyc3fGIqLtmj4P0bPMyAQw5D6zQtdzvn4r//Go+0lqXYtZU/pe1TymlgOEhKFNkCnGqi7kTs+D1aq6eNpGrv3pRxycEwrf5KFQqi81/5BCwqunPkMlgdTh0qQkOboT++aFtKx86HTb83QSgrIm2XXZuQTbRDdXMfONTXm08j5SU1NCNqsmZzMtRc7hWv8lSdTZrtQk8QY+GKXoT/nEDZIwyASE+tdOX8LsTuVeSaQqbfBN88hisfx62LIE+g2DkLBh8rqkxJmefzpDceMr8O6/YB6WF/DX+LSY0baaPqqVBR7HGncd/NlzL5oRzWXzH9QCcD9yX0X6/QzhGRLXGn6Dwh9Z2aq01sNPe4kSwUGdKPVUFTfWhm7gGZpx6rz6hzX9UWQJDLgzJpa1q+p7oTGa51gD2D8fD3QSHt8CUm+y5XluGen5Gu1fZGhQALlcfAaeY98N/Z1524G3l/vhNzZVMjV3Do7F/4NJTC6kiAQhikMTWZfDSTeZBfP0SiOtrY2nbkTwAZi+Ei/7d9JNsew0Kn4PVT50+JqoXoM3/Ux8FCbn8q/ps3m6YwEfucZwg3tQC5uQ3O86fQOXkiKi2+BMU/rejA5RSyhMkuq9QJ8UL9WxmS1Jm6DKlnqyE+uqQDUe1HjR7dX/SVDXJ1FBFgr2jtI7uhMa60PUnWJL6Q8Zos27z+XfZd12tYe2zpqaTNcm+67YhOSWNOytv55+xv+H+mKe5veEOQAXWh7HxH7hfuZUtDOfa3T8i6dFC59+UY+Nh4rfNn8Z6KNtmXhKqD5n+RFcUxMSbf+N9ciBzHEnxqUQVlrBpRRE1FXWha/51iF9ZUpVSLwOvehbXAUApFYupFX0Hk0n12ZCUMFJ4+xRCFRSsvEchyJDqKzEzdGsqVHqGo4Zo4prVibpXmw7HXHWIjXqYvZ2oNs9kbtfwmfDFkyYVeK8ke65ZWmhSdFz6UOiHimK1g9fzu8ZvsiBmMTvcr/CU65ud78NY+yz6tbv4Qo/hu6d+Si29qQp35tLoWJON14+MvJHwhm8Xf3q4ZmPyHL2glDqolNqqlNqNaTq6BnhYa/1sCMsYGaJ7mZmoIasphDAZnq/EEOY/8g5HDU1NwUrQtkeb1M+56pD9nagHN5hmgvSR9l2zLSNnm6YIOyeyrX0WouNMB6oDrHUZXkv8Fi81fYWfxLzEc1O+9P8B6XbDqkXw2p185irgO6fmU8vpkWuSudR5HdYUtNYngT8Cf1RKxQDpQF2PWHGtpYS00PUpOB0UtLb/TbLKms0cmjcm60Hz6FvgPqkYH1/ORZfZOxyPQxshc4wZnRJqg84xfTw73oLR/xb89U5Vw+aXYdw3oHef4K/nJ+9bcuMF8PxVTF7/S8iO77hfpvYYvH4XbH0VJl7LjZ9dwinO/LmHM3NpqJPPRSJ/Vl77jlLqqFLqGPA0cKJHBgQIbVI8q0/BaqYKlaRMaPAs5mO3yhKTYjzRv/HkgZhbkM17987G1Sebm0e77f0PqrVtayj4JSrGNCHt+Jd5Yw7W5pfN73Xyd4K/ViCiY+GaxTD8YvOwX3qbefC35Habsv7hbNj2Olzyn/D1P5Ce0nraknBlLrUrTXdX40/z0b8DXwVGAfuB+0NaokgWwlQX1Bw1b3ehzjwZyrWaq0pMausoBybKp4UgW2rlAThZ4Ux/gmXkbJOP6mBhcNfRGtY8A/3GQM5Ue8oWiNh4mPcCfOX/mSG3j+TDq7fBmr9A4fNmJvHvp8BL3zUvKLe8B+f9GJTyNg/6Cmfm0p66EI8//3urtNbWv9h/V0p9HsoCRbT4NDMSIRRqykLfdATN12pOt3k+QWVx6LKjtpQ6zAwdtJN3JnPol/r0Gn4xKBcUvQXBDB8tXWf6QxzqYG5XVDRc9AvTjPXRw7DVM7wTTE1y8HlmbsCYuc1eICJlnL6lpy7E409QGKCUugXYBmyHVhr9egqrTyEU7fFOBQXvWs0hqClUFocuH1FLacPMEMHaYwFNbmrVoY3mAZ051p7r+SMhDQadZwLcjJ8H/u9q9TMQkwDjr7a1eEHpNwqu/DM0NfLW5xt4ZtV2NlbFk3aoD/Mb85jbSo0ykkbxdOWFeILhT/PRfwDjgf8EioBxSqk3lFILlVLXhLJwESc+DRpPmjZ5u4UgfUarvLOabQ4KbrdZSyFUKbNb8mZLtbEJ6eBGSBtxeuaqU8ZdCUd3nF7trbPqjps2+vHfDG068QAt3XiYu98s44uqFE4S22Xa5iOtOcsp/gSFp7TWt2utL9RapwJDgd8DFcClnnQXPUO856Edin4Fp2oKcX0hKtb+mkLtUWg65VxQ8M2WapeDG5ztT7CMmQuuaHa882xgC9JvWGwm3E35bkiLGaiu2jZvDbfNTolDAdkpcSy80ubRbhHIn+ajt5VSR4GlwHKtdbFnJFIcJkneOsChNoMw885qPgp9B9t3XXeTGdXkRFBQKjRzFUKdMrullMGAMusq2KHmKFSX+jVRyXYJaRxKP5fEna9SemoGGpf/eZ2sDubsKeEpux+6ctt8JDVnOaXDmoLWeibwGyAXWK6U+gx4B5Mx9WGtdc8ICOCT/6iVYXbBqDsO6NM1kVBL7Gd/UKgK7WzmM8T0hpSBJi2FHQ5uMNswPVj/fHwSWeook9Tp+/HrbXrHCtP0dNYtIS5h4Npqg+/ubfNdlV85e7XWW7XWC7XWFwDTtdbnaq3/Q2u9LsTliyy+6bPt5FSKC0tif/v7FLwpLgbae932ZIwyD0Q7WEGhf377x4XIP6vHU6t7cXXUqmb7232b1ho++h+T0XPclSEtXzB6att8V9XpRO6eGc49U6iS4jk1m9mS2M/+pHiVB0x6iFBPvvOVMcrUFJoag7/WoY2mScqprJwt9ElJ5ZWm8/la1CekUO3d3+7b9L5P4MDnZpy/EzOwA9RT2+a7KseW41RKzQYeBaKAp7XWi9o4birwGXC11volp8rnl959wBVtf6oLp2sKSf1NYGtqsO9hUlViOpmdHHeQMcp0bh/fG/yci3B1MnvMn5XHM0tmc616h6ujVvHnpn9r9W3aN+3CC/EPUhCbSq+C68JT6E7oiW3zXZV9Sz61QykVhVmXYQ4wBrhGKTWmjeP+G1jhRLk6TanQpLrwps0O0brJLXmX5Qx+XYWlhSVMW/Qu6zZtYnVFgrPDDDM8q1aVbQ/uOierzNDWMHbUzi3I5rtXXsZq13hujn6DYX3UGW/TvmkXprk2cY67kMdOzmHpluNhK7fofhwJCsBZwC6t9W6tdT2wGPh6K8f9GHgZCOEqMEGKT7c/ffaJI2amp1NNFzat1ez7kMpS5eytT3F2/HmGJ5NpsEHBmh/g5EzmVswtyGbqjQ+RoSp55/ztZ7xZW0M7o2nkl9HPsc/dj6fqL4n4oZ2ia3EqKGQDB3y+Lvbs81JKZQNXAE90dDGl1C1KqTVKqTVlZWW2FrRD8amh6VNISLd1rd522bRWs+9Dqh8VlJLu7PjzXkmmY7ssyO/nXUMhAoZ0Djrb5EP68H/g+L5mH1mdzndFv8wo1wHua7yeemK6xNBO0XU4FRRaa2huuVLbI8A9WuumVo5tfqLWT2qtp2itp2RkONQ5a0lID02fglOdzGDbWs3Ww6i/Oo5LaUp0WrP9jsjIM6tjBePgBvMzsYJluM15wGxf+YFZD9gjKyWOma61/ChqGYsbp7PSPdm7Xwi7OBUUigHfsYo5QGmLY6YAi5VSe4GrMOs3zHWkdJ0Rkj4Fh4OC1XcRZJ+C9TDKwgTJUp3ebL8jrBFI7g7fJdp2cGNk1BIsfQfD5Q/D/k/h798y82Lcbh4bvZ0/xTzCJj2E3zbeAMjQTmE/p0YfrQZGKKWGACXAPODbvgdorYdYf1dKPQu8rrVe6lD5/BefbiabNTXalyL6xBEYOKTj4+wSHQtxqUEPSzVLMW4iq8kEyVKd5vxDKiPP5KOq2AepQzt/fkOd6ZPIm2N/2YIx/ptmZNWyO+ChkRAbz+STlZSnTeKemrupq4/q8msBi8jkSFDQWjcqpW7HjCqKAp7RWm9RSt3q+bzDfoSIYY3Drzt2Og11sGqOOltTADMsNciagvUwOrx8GTQCyTksnO3w+POM0WZbVhRYUDi8FXRTWIejtqngOpNOe+OLcLISBk8jbewVvOWK6vhcIQLk2DwFrfUbwBst9rUaDLTWNzpRpoAk+ExgsyMo1NdAQw0kOhwUEjOD7lMAT2Ao6QVbUnn3njC8bWd4aiWHtwT2tl+yxmyzp9hXJjv1Gw0X/0e4SyF6EKf6FLoPu1NdOD2b2ZI0AKoP2nOtymLnsqO21DsZ+uZ2Ou20Nb9i6etLOUIaS21exE2IrsqxmkK3YSWts6uz2emJa5bkLFNTsKNvpLIY+jrYJ9JS//zTw0r9YM2vqGtoYmLsl6xpGuZfRlIhegCpKXSWN1OqTTUFq13fqRQXlj7Zpi29xoZ5guGsKYCZdHZst5mZ7AdrfkUqVeS6DrPePaxL5PcXwgkSFDorzrP0o13ps63mI7s6rf1lrXtQ1XJkcCedrIRTVeENClYn8eHNfh1uzaOY6NoFQKF7RLP9QvRkEhQ6KzoWevWxsU/B86bu1FoKluQss60KMiVFpcPrKLSmvycoHPSvCcmaR1Hg2kWjdrHJMxpaJoEJIUEhMHamuqg5Cr2SzaIxTrKrpmCtuObkOgotJfU3QdVaE6EDVn7/iWoX2/QgTtJLJoEJ4SFBIRB2pro4ccT5kUfA0u21nCKWJ1//sHPrAbdU6UlpFc7mI6XMeH5reGkH5hZks/CKsRRE7Wa9e7jk9xfCh4w+CkR8WvDNLhanU1zgGX3zymYmqL70V8f8Xw+4NZXFZo2JxDDnDRo4FXauMH098akdHj53wDGgluu/dTXXj78o9OUToouQmkIg7EyfXVPm+MQ1a/TNIZ1Gf2U6zAMefVNZbPonwj3LduDZZlvsX22BPR+Ybe4FoSmPEF2UBIVAWH0KumWi1wCEoaZgjbI5SCoDOHbG/k6pKglvf4IlaxIoFxR/4d/xez6EtOGQPCC05RKii5GgEIiEDJOs7JR/4+Lb1NRomjscDgrWKJtDOpVMdQyFu9n+Tqk8cLrTOpx6JULmWDjgR1BoajTrG0stQYgzSFAIhDWn4ESQC/zUlgPa8aBgjb45qFOJVU2kUR3Y6Bt3kxm9FM7hqL4GnQvFq5utQdCqg+uhvhqGSFAQoiUJCoGwHuLBzgYOU96juQXZLLwyn/r4/gCMT64JbPRNVSm4G03uoUgw7CJoqDXrELRnx1umqWnoDGfKJUQXIkEhEIn2LFLjDSpOz2bGBIb/vmk2AM9cmRXYcMzje802UoJC7gXgioFd77R/XNGbplbhxyglIXoaCQqBsJLX1QTZfORNhuf8PAUg+AlsVlBIGWxLcYLWKxEGnQNfvtv2MRX7TTqMSFtUR4gIIUEhEAnppvkh2JpCuJLhWeLTzZt1oHMuju8FFRXeiWstDZ9pHvoV+1v/fPtys8271LkyCdGFSFAIhCvKTGCzo0/BFQO9U2wpVqe5XGZIZmWAQaFinwkIUTH2lisYY68w240vnvmZ1rDubzBgIqQNc7RYQnQVEhQCldDPhj4FzxwFpewpUyD6DDqdqqKzju81i8xHkr65pm9h/d/PnEdSvAaObIEpN4WlaEJ0BRIUApWYYU9QcHoZzpb6Dobj+wI79/jeyOlk9rG27xw4tpt5P/9d87xOq5+G2EQY943wFlCICCZBIVCJmcE3H4UpGV4zKYPMspwdje1vqb7GBLUICwpLC0u4eU0OZTqZu2JeoqSilnuXbOKd91fBpn/ApO9Ar6RwF1OIiCVBIVAJGWbyWjCpLk4chsT+9pUpECmDAH06Bba/rNpFhAWFB1cUcbwhmkcar+Ic1zZujXqNuIbjDF51B8T1hQt+Gu4iChHRJCgEKrEfNNZB/YnAznc3maCQFAlBgbZH67TFOxw1187SBM3K3/T3potY3nQWC2IWs673rQx0l8CVT0FCWphLKERkk9TZgUrwmcAWSHNEzVHQ7ggICp6O4opO9itURGZNISsljpKKOjQu7mj4Me+5P2KwOszn8dN5bvjMcBdPiIgnNYVAWR3EgXY2Vx8023AHhaQBZj2EQGoKsUkRNyvYyusE0EQULzVdyB/VNVw155Iwl0yIrkFqCoGyFpUJtLP5xGHPdcIcFKKizczmQIJC39zwDqdthZWu48EVRZRW1JGVEsf8WXmyqpoQfpKgEKiEIPMfRUpNAcyw1M4GhWN7IH1EaMoTpLkF2RIEhAiQNB8FKj4NUIHnP6q2agphXsYSTL/CsT3+H9/UCMd2m0VqhBDdigSFQEVFm8BgNQN11olD5vzoWHvLFYi0YaYZ7KSfiwZV7AN3Q8TWFIQQgZOgEIzEzMAX2qk+FP7+BIv1xn/sS/+OL9/lOU+CghDdjQSFYCRmBN7RXH0oMvoT4HRQKPczKBzdYbZSUxCi25GgEIzEzNN9A50VSUGh7xBAna4BdOToTtP0FWHDUYUQwXMsKCilZiulipRSu5RSC1r5/Fql1EbPn0+UUhOcKlvAkvqbvoHOprpwuyNjNrMlpjekDOxcUEgfGdoyCSHCwpGgoJSKAv4AzAHGANcopca0OGwPcKHWejxwH/CkE2ULSlIWNNVD7bHOnVd7FHRT5PQpAKQO8z8olO+UkUdCdFNO1RTOAnZprXdrreuBxcDXfQ/QWn+itT7u+fIzIIKW82qD9aZf3cnlLKsPNT8/EqQNN30KHdV6ao+ZYbhSUxCiW3IqKGQDviu5FHv2teV7wJttfaiUukUptUYptaasLMh1koORnGW21kPeX5EaFE5VdTzE9vBms80cG/oyCSEc51RQaC0XQquvpEqpGZigcE9bF9NaP6m1nqK1npKREcb1CKyHemcXvj8RgUEh09OaZz3023Jok9n2zw9teYQQYeFUUCgGBvp8nQOc8SRVSo0Hnga+rrUud6hsgbP6BKyUFf6yagqRMJvZkjnObA91FBQ2m3In9gt9mYQQjnMqKKwGRiilhiilYoF5wDLfA5RSg4AlwPVa6x0OlSs40bFmsZ1AgkJcKkT3Ck25AhGfahLjdVRTOLzpdAARQnQ7jiTE01o3KqVuB1YAUcAzWustSqlbPZ8/AfwKSAP+qEzmzUat9RQnyheUpP5Q1cmgUFV6uj8ikmSOa7+m0FgPR7bDuRc5VyYhhKMcy5KqtX4DeKPFvid8/n4zcLNT5bFNUlbnRx9VFkOfCBxc1X8c7Fpp1mturRZzdIfJedR/vPNlE0I4QmY0Byupf+dHH1UeiMygkDnOzJ8o297658WrzXbARMeKJIRwlgSFYCVnmXH7jfX+HX+qGk5WRGZQyCowW+vh39L+T00fStow58okhHCUBIUgLC0sYeFHJt30VQ+8zNLCko5PqvQc02dg+8eFQ99c0xy29+PWP9/3KQw+L+JWWxNC2EeCQoCWFpZw75JNbK5NBiD6RDH3LtnUcWCoLDbbSKwpKAW558Pej86c2VxxACr3w6DzwlM2IYQjJCgE6MEVRdQ1NFGi0wHIUWXUNTTx4Iqi9k+s9EzsjsSgACYo1BwxSe987f/UbAdLUBCiO5OgEKDSijqz9QSFbI4229+mymJQUZGVDM9X7vlmu/fD5vt3rYS4vpLeQohuzrEhqd1NVkocJRV11BPDYZ1Ctjrq3d+uymLTOR0VoT/61KGQMgiK3mBp9GweXFFEWUU1a3u/xrGBlzDYFRXuEgohQkhqCgGaPyuPuBjzgCzR6WSro8TFRDF/Vl77J1aVRG7TEZh+hXHfwP3le/xuyYeUVNTxVdcakqjlP/eO9q8zXQjRZUlQCNDcgmwWXplPdkocJTqdwdHlLLwyn7kF7SV/xSx6H8lBAaDgetBurtWvoXBzS/Tr7HdnsLJhbMd9JkKILi1C2zC6hrkF2SYIvP0+fLqWnAkD2j+hsd40H42f50wBA5U2jKVN07g56g3GqH1McO3m7vofonF13GcihOjSJCjYIWWgSf9w4lD7OY0q9oN2m3b7CPdE/K0MqD3Gua6tPNV4Ka+4TQd0h30mQoguTYKCHfoMMtuKA+0HhWO7zbYLBIUfzZ7Ed5f8mlOnGnB7Whn96jMRQnRp0qdgh765Zms99NvShYKC1WcyICUBBWSnxPnXZyKE6NKkpmCHvrlm7kFHC98f2w2xSZCQ7kixguXtMxFC9BhSU7BDdKwJDOU72z/u+B5IHSK5g4QQEUuCgl3ShsNRP2oKXaDpSAjRc0lQsEv6CDj2JbjdrX/uboLj+0xNQQghIpQEBbukDYfGk1BV3Prnx/aYYatpI5wtlxBCdIIEBbukex72LbOLWg5vMltJKCeEiGASFOxi1QBaGYG0tLCEZ5csp1G7mPHXQ5I/SAgRsSQo2CWxn0ktfXhzs93WYjzZ9V+yWw9gT2WTf4vxCCFEGEhQsItSZo3jksJmu63FeEa79rNNDwbwbzEeIYQIAwkKdsqaBEe2QsPppHGlFXUkU0OOOsp296Bm+4UQItJIULBTVgHoJji06fSulDhGqf0AbNMDm+0XQohII0HBTtmTzLZknXfX/Fl5nB1jOp83uocBklhOCBG5JPeRnZKzzNrLpaeDwtyCbA59VsKeIzkcJ5nslDjmz8qTnEJCiIgkQcFug86BPR+A1qbzubGe/scLYdJV7Pm3y8JdOiGEaJc0H9ltxCVQfRAObTRf7/sI6qth5KzwlksIIfwgQcFuIy4BVzRs/If5etNLEBMPQy4Mb7mEEMIPEhTslpgBoy6DwuegZK0JChOugdj4cJdMCCE6JEEhFKbfCw218NRFEN0Lzr8r3CUSQgi/OBYUlFKzlVJFSqldSqkFrXyulFKPeT7fqJSa5FTZbNdvNFz/Ckz9Ply/FFIGdXiKEEJEAkdGHymlooA/AF8FioHVSqllWuutPofNAUZ4/pwN/Mmz7Zpyzzd/hBCiC3GqpnAWsEtrvVtrXQ8sBr7e4pivA/+njc+AFKXUAIfKJ4QQAueCQjZwwOfrYs++zh4DgFLqFqXUGqXUmrKyMlsLKoQQPZlTQaG1lep1AMeYnVo/qbWeorWekpGREXThhBBCGE4FhWJgoM/XOUBpAMcIIYQIIaeCwmpghFJqiFIqFpgHLGtxzDLgBs8opHOASq31QYfKJ4QQAodGH2mtG5VStwMrgCjgGa31FqXUrZ7PnwDeAC4FdgG1wE1OlE0IIcRpjiXE01q/gXnw++57wufvGrjNqfIIIYQ4k8xoFkII4aXMC3rXpZQqA/YFeHo6cNTG4nQVPfG+5Z57jp5434Hc82Ct9RnDN7t8UAiGUmqN1npKuMvhtJ5433LPPUdPvG8771maj4QQQnhJUBBCCOHV04PCk+EuQJj0xPuWe+45euJ923bPPbpPQQghRHM9vaYghBDChwQFIYQQXj0yKHS0Clx3oZQaqJR6Tym1TSm1RSl1p2d/qlLqbaXUTs+2b7jLajelVJRSqlAp9brn655wzylKqZeUUts9v/Nzu/t9K6Xu9vzb3qyUekEp1bs73rNS6hml1BGl1GaffW3ep1LqXs/zrUgpNasz36vHBQWfVeDmAGOAa5RSY8JbqpBpBH6qtR4NnAPc5rnXBcA7WusRwDuer7ubO4FtPl/3hHt+FHhLaz0KmIC5/25730qpbOAOYIrWehwmr9o8uuc9PwvMbrGv1fv0/B+fB4z1nPNHz3PPLz0uKODfKnDdgtb6oNZ6nefv1ZiHRDbmfv/qOeyvwNywFDBElFI5wGXA0z67u/s9JwNfAf4XQGtdr7WuoJvfNyZ/W5xSKhqIx6Tb73b3rLX+ADjWYndb9/l1YLHW+pTWeg8myehZ/n6vnhgU/F7hrTtRSuUCBcDnQKaVltyz7RfGooXCI8D/A9w++7r7PQ8FyoC/eJrNnlZKJdCN71trXQI8BOwHDmLS7f+LbnzPLbR1n0E943piUPB7hbfuQimVCLwM3KW1rgp3eUJJKXU5cERrvTbcZXFYNDAJ+JPWugCooXs0m7TJ04b+dWAIkAUkKKWuC2+pIkJQz7ieGBR61ApvSqkYTEB4Xmu9xLP7sFJqgOfzAcCRcJUvBKYBX1NK7cU0DV6klHqO7n3PYP5dF2utP/d8/RImSHTn+74Y2KO1LtNaNwBLgPPo3vfsq637DOoZ1xODgj+rwHULSimFaWPeprX+H5+PlgHf8fz9O8CrTpctVLTW92qtc7TWuZjf7bta6+voxvcMoLU+BBxQSuV5ds0EttK973s/cI5SKt7zb30mpt+sO9+zr7bucxkwTynVSyk1BBgBfOH3VbXWPe4PZoW3HcCXwC/CXZ4Q3uf5mGrjRmC958+lQBpmtMJOzzY13GUN0f1PB173/L3b3zMwEVjj+X0vBfp29/sGfgNsBzYDfwN6dcd7Bl7A9Js0YGoC32vvPoFfeJ5vRcCcznwvSXMhhBDCqyc2HwkhhGiDBAUhhBBeEhSEEEJ4SVAQQgjhJUFBCCGElwQFIYQQXhIUhBBCeElQEMImSqk/K6WmhbscQgRDgoIQ9jkb+CzchRAiGBIUhGiHUirZk4p6i1KqVim1Xin1mVLK1eK40cAOrXWTz75czypoT3tWBnteKXWxUupjz2pZZ3mOS1BKLVdKbfAcd7XDtymEV3S4CyBEJNMm1XiB5wH+C611WwsyzQHeamX/cOCbwC2YZIzfxuSk+hrwc8zCKLOBUq31ZQBKqT523oMQnSE1BSH8Mw7Y0s7ns2g9KOzRWm/SWrs957+jTcKxTUCu55hNwMVKqf9WSl2gta60sdxCdIoEBSH8MwaTifMMSql4IEVr3VrO+lM+f3f7fO3GU1PXWu8AJmOCw0Kl1K/sKrQQnSVBQQj/ZAGH2vhsBvBeoBdWSmUBtVrr5zDLS04K9FpCBEv6FITwzwrgf5VSN2qt32/x2RzMSmeBygceVEq5MfnyfxjEtYQIiqynIESQlFLrgLO1WRJSiC5NgoIQQggv6VMQQgjhJUFBCCGElwQFIYQQXhIUhBBCeElQEEII4SVBQQghhJcEBSGEEF7/H35XyukWFLV/AAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(precession_times_loaded, probabilities_fid_loaded,\n",
    "            label=\"Simulation Data\")\n",
    "plt.plot(precession_times_fine, T2_star_decay(precession_times_fine, *popt),\n",
    "         color=\"tab:orange\",\n",
    "         label=\"Fit: \\n$T_2^*$ = {:.2f} ms\\n \".format(popt[0]))\n",
    "plt.xlabel(\"$\\\\tau$ / ms\")\n",
    "plt.ylabel(\"$\\langle P(T_0) \\\\rangle$\")\n",
    "plt.legend()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The simulated $T_2^*$ time lies in the expected order of magnitude when\n",
    "approximated by\n",
    "$(T_2^*)^2 = \\frac{2}{S_0 (dJ/d\\varepsilon)^2\\ln\\frac{T_{meas}}{2\\pi \\tau}}$\n",
    "(c.f., supplement of\n",
    "[Dial et al., PRL 110 (2013)](https://doi.org/10.1103/PhysRevLett.110.146804))\n",
    "leading to $T_2^* \\approx 54~\\text{ms}$ for our parameters. We approximated the\n",
    "measurement time by the number of noise traces\n",
    "$T_{meas} \\approx n_{traces} \\tau$.\n",
    "A slight deviation from the Gaussian decay is caused by the finite\n",
    "frequency cutoff, the non-linearity of $J$ and stochastic fluctuations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Spin-Echo\n",
    "The coherence time can be increased by intercepting periods, during which the\n",
    "qubit is susceptible to noise, with appropriate echo pulses. For example, a\n",
    "$X_{\\pi}$ pulse can be placed in the middle of the free precession period of\n",
    "the FID sequence. This allows to refocus over- or underrotation of the qubit\n",
    "that occurred due to noise that is quasi-static with respect to the precession\n",
    "time $\\tau$.\n",
    "\n",
    "In the noise-less case, the qubit would return to its initial position after\n",
    "the pulse sequence. We therfore expect the average probability of the $| S\n",
    "\\rangle$ state to decrease as $\\langle P(S) \\rangle = \\frac{1}{2} \\left(1 +\n",
    "e^{-(t/T_{2,\\text{echo}})^2} \\right)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAApfElEQVR4nO3dfZxcZX338c93JkQqYhUJKkkgKU2L1ALiFnygWlRsYluDLVUoClRsSm/Rtne9W6y2tU9Wa29qrWiaKlWriFaNd7QRUKs3VcAmUASCoDGAWQMSnpTnJDO//nHObGaXmd2zm50z+wvf9+u1r+zMnDPzy+yZuc51ne85lyICMzOzqTSGXYCZmeXgBsPMzCpxg2FmZpW4wTAzs0rcYJiZWSVuMMzMrBI3GDZQkr4g6YyaX/MWSS+p8zWHQdISSSFp3hyo5RckjXbd3iTpF2bx+R8Tf9O5zg2GTUnS8ZIul/RDSXdL+rqkn6uybkSsiIgPz/B1b5H0kKT7u37eO5PnmilJ8yW9TdJ3JD1Q1nSBpCWz8Ny1fAlK+qqkeyQ9btCv1RERPxMRXy1f/22SPlrXa9vguMGwSUl6IvB54B+BA4CFwJ8Dj9RUwq9ExBO6fs6p6XU7PgW8HPgN4MeBo4CrgBcP+oVno+dQNmw/DwTF/8Nsxtxg2FR+CiAiPh4RrYh4KCIujYhrASSdWfY4/rHsgdwoaezLtNy7fV3Xsl+T9HflHu/NklbMtDBJvyXpW5Luk3SDpGO6Hj5a0rVlTZ+QtO+E9TaXvaV1kg7u8/wvAU4EVkbEhojYFRE/jIjzI+KD5TIHl89xd/mcv9W1/tskfVLSR8oaN0kaKR/7V+AQ4HNlz+kPu4aYzpL0PeA/JDUkvVXSrZLuKJ/rx6fxNp0OXAl8CBg3NCjpQ5LeVw4b3l/+HZ8m6d3l3+dGSc/qWv4WSW8u3+t7JP1L9/s64blvkfQSScuBPwZeVb7GN7sfn/BefbTr9mvK//Ndkt4y4bkbks6V9N3y8U9KOmAa74nNkBsMm8q3gZakD0taIenJPZY5DtgCHAj8GfCZST7AxwE3lcv+LfBBSZpuUZJ+HXgbxRfiEyn2nu/qWuSVwHJgKXAkcGa53ouAvykffzpwK3BRn5d5CfBfEbF1klI+DowCBwMnA2/vbjDLui4CngSsA94LEBGvAb7H7h7U33at80LgGcAvlnWfCZwA/ATwhM5zVHQ68LHy5xclPXXC468E3krx93gEuAK4urz9KeC8CcufVtZ1GMXOxFsne/GIuBh4O/CJ8v951FQFSzoCeD/wGor39SnAoq5F3gicRPE+HQzcA5w/1fPannODYZOKiB8Bx1MMafwzsL3co+7+4rkDeHdE7IyIT1A0CL/U5ylvjYh/jogW8GGKL+2JX2LdPivp3q6fzh7864C/Lff8IyI2R8StXeu9JyK2RcTdwOeAo8v7TwMuiIirI+IR4M3Ac/sck3gKcFu/wiQtpnhv/igiHo6Ia4APUHzRdXwtItaX/99/pRjSmsrbIuKBiHiorPe8iNgSEfeX9Z5SZbhK0vHAocAnI+Iq4LsUQ2vd1kbEVRHxMLAWeDgiPlLW+wngWROWf29EbC3f178GTq3w/5muk4HPR8Rl5d/oT4B21+O/DbwlIkbLx98GnDwbQ3g2OTcYNqWI+FZEnBkRi4BnUuzVvbtrke/H+KtY3lou08vtXc/7YPnrEyZ5+ZMi4kldP/9c3r+Y4guwn9u7fn+w6zUOLuvr1HA/Rc9kYY/nuIuiQevnYODuiLiv675bJzzXxDr2rfDF1t2jGVdv+fs8Jm9kO84ALo2IO8vbFzJhWAr4QdfvD/W4PfFv013bZH/nPXFw9+tExAOM7z0eCqzt7EQA3wJaVHtPbA+4wbBpiYgbKcbDn9l198IJw0qHANsGXMpWimGR6dpG8YUDgKT9KHoS3++x7JeAYyUt6vFY57kOkLR/132H9HmuXvpdKrr7/nH1ls+/i/Ff7I8i6ccohpteKOl2SbcDvw8cJalKL6efxRNqqfJ37vX/fAB4fNftp3X9flv360h6PMXfqGMrsGLCjsS+EVH1fbcZcoNhk5J0uKQ/6HxplsMwp1IcSO04CHijpH3KYwvPANYPuLQPAG+S9GwVflLSoVOuVexl/6ako1XETN8OfCMibpm4YER8Cfgixd7ssyXNk7S/pLMlvbY8tnE58DeS9pV0JHAWxfGCKn5AcVxiMh8Hfl/SUklPYPfxgF1TrHcSxV73ERTDcUdT/F3+k+K4xky9XtKi8hjVH1MMW03lB8ASSd3fN9dQDK3tUwYBTu567FPAL6uIc88H/oLx31Wrgb/u/L0lLZC0cub/JavKDYZN5T6KA9XfkPQARUNxPfAHXct8A1gG3Ekxrn1yRNw18YlmqJMi6vysBYiIfytf68Kyxs9SxH4nFRFfphgT/zTFnuxhwCmTrHIyReP3CeCHFP/3EYreBxSN5xKKPe21wJ9FxBcr/t/+BnhrObTypj7LXEBx7OMy4GbgYeANFZ77DOBfIuJ7EXF754figPlpezDefyFwKUXIYQvwVxXW+bfy37skXV3+/icU7/09FDHtCzsLR8Qm4PXlfbeVy4ydFAj8A0WA4FJJ91Fsk8fN8P9j0yBPoGR7QtKZwOsi4vhh12KDJekWir/1l6Za1vZO7mGYmVklbjDMzKwSD0mZmVkl7mGYmVkle/WZkQceeGAsWbJk2GWYmaVx1VVX3RkRC3o9tlc3GEuWLGHjxo3DLsPMLA1Jt/Z7zENSZmZWSa0NhqTlkm5ScRnocydZ7ucktSSdPN11zcxsMGprMCQ1KS5BvILicgWnlpcx7rXcO4FLpruumZkNTp09jGOBzeVlmndQzBHQ6/ovb6C4bMMdM1jXzMwGpM4GYyHjL408yoRLSktaCLyC4uJi01rXzMwGq84Go9esahPPGnw3xWQ0rRmsWyworZK0UdLG7du3T79KMzPrqc5Y7Sjjr6W/iEdfS38EuKicWuFA4GWSdlVcF4CIWAOsARgZGfFp7GZms6TOBmMDsEzSUooJZk5hwnSREbG087ukD1FM0/jZ8lLMk66bwSc2fI/v3/PQsMuwil740wfx7EN7TWE+d1255S4u33zn1AvanLDkwP341WP6zc8199TWYETELknnUKSfmhTzKm+SdHb5+MTjFlOuW0fds+W+h3fyR5++DgD1GmCzOSUC/nvrvfzrWbmmWXjHF27kmq33ehtLIKL4Lnj5UQczr5njlLhaz/SOiPVMmImtX0MREWdOtW4mO1vF6Nifv/xnOON5S4ZbjE3plf90BTt2tYddxrTtbLV58eEH8cEzf27YpdgUzv/KZt51yU20ItJcciNHs7YXaJdXBW54zy+Fhoo9wGzaAXL3IoXOnynTduYGoybtdtlguMVIodkQrUyf5FK7HSQZ3XjMa5YtRqudZzvzplWTzpdP03t/KTSkVB/kjlYETe+UpND5O2XaMXGDUZOWexipNBsaG0bMpN0OGt4pSaHzd2on2jFxg1GTdnn81D2MHJruYdiAjfUwEm1nbjBqMjYk5Q9zCo1G0gajHd4pSaLhISnrpzO84c9yDllTUuGUVBoNp6Ssn844pXsYOWRNSbWckkrDKSnryympXBpSqoORHT6GkUfDxzCsH6ekcsnaw3BKKo/OzmOmNJ4bjJo4JZWLU1I2aE5JWV9OSeXSaCQdknIPI43OaIN7GPYoTknl0lBxXaZsInCDkURn3zHTduYGoyZOSeWS9RiGU1J5OCVlfXU2Ch/DyCFzSsrBihyckrK+Onur/jDnkLWH0faZ3mk4JWV9jaWk3GCk4KvV2qA5JTUFScsl3SRps6Rzezy+UtK1kq6RtFHS8V2P3SLpus5jddY9G8Z6GN77S6GZMCUVET7onUjGlFRtMwNKagLnAycCo8AGSesi4oauxb4MrIuIkHQk8Eng8K7HT4iIlDPce8a9XDKmpDr1usHIwSmpyR0LbI6ILRGxA7gIWNm9QETcHzHW3O4HJHorJ+eUVC6NhMcwxoIVHmhOwSmpyS0EtnbdHi3vG0fSKyTdCPw78NquhwK4VNJVklb1exFJq8rhrI3bt2+fpdL33NilQbz3l0IzYUqq7WBFKmNDUom2szobjF5b8aPeqYhYGxGHAycBf9n10PMj4hhgBfB6SS/o9SIRsSYiRiJiZMGCBbNQ9uxo+0zvVDKmpBzdzsVTtE5uFFjcdXsRsK3fwhFxGXCYpAPL29vKf+8A1lIMcaXRckoqlYZERHEgOQtffiaXhoekJrUBWCZpqaT5wCnAuu4FJP2kytlfJB0DzAfukrSfpP3L+/cDXgpcX2Pte8wpqVwyRh7bHvZMpemUVH8RsUvSOcAlQBO4ICI2STq7fHw18GvA6ZJ2Ag8BryoTU08F1pZtyTzgwoi4uK7aZ0M4JZVKxgTL7pTUcOuwasa2sfZw65iO2hoMgIhYD6yfcN/qrt/fCbyzx3pbgKMGXuAAtZySSiVjRt7bWC5jQ1KJtjEH8GrilFQuGSOPTknl0nRKyvpxSiqXjAkWp6RyybiNucGoiVNSuXR6gpn2/jwNcC5OSVlfTknlkjIlFe5hZJIxJeUGoya7I49DLsQqSZ2S8qc6hYwpKW9aNfExjFwyp6Tci83BKSnry+PLuWROSXmnJAenpKwvjy/nknH6TKekcnFKyvpySiqXjNNnuhebS8YknhuMmrSdkkrFKSkbtIzbmBuMmrSckkpFTknZgGVM4nnTqokPSOaSMSPvlFQuGZN4bjBq0m4HEsgf5hSckrJBy7iNucGoSSvCY8uJOCVlg+aUlPXVaju9kknGlFTbKalUnJKyvtruYaSSMcHiKVpz2b2NDbmQaai1wZC0XNJNkjZLOrfH4yslXSvpGkkbJR1fdd25rtUOJ6QSSZ2S8naWwu6UVJ6NrLYGQ1ITOB9YARwBnCrpiAmLfRk4KiKOBl4LfGAa685p7QgPFSSSMSXlOb1zkYSUaxurs4dxLLA5IrZExA7gImBl9wIRcX/E2Lu3HxBV153r2u3wUEEiGRMsnqI1n6aUahurs8FYCGztuj1a3jeOpFdIuhH4d4peRuV1y/VXlcNZG7dv3z4rhc8Gp6RyaSS8MJznXMmn0ZBTUn302oof9U5FxNqIOBw4CfjL6axbrr8mIkYiYmTBggUzrXXWOSWVS8bIY9s9jHSaUqqdkjobjFFgcdftRcC2fgtHxGXAYZIOnO66c1G77R5GJhmnz3RKKp9mQ05J9bEBWCZpqaT5wCnAuu4FJP2kylOhJR0DzAfuqrLuXNcKp6Qy6fytEnUwnJJKKNtB73l1vVBE7JJ0DnAJ0AQuiIhNks4uH18N/BpwuqSdwEPAq8qD4D3Xrav22eCUVC4Zz8NwSiqfZkNuMPqJiPXA+gn3re76/Z3AO6uum4lTUrlknD7TKal8nJKynlrha/xkknH6TKek8mkk62G4wahJu+0hqUyckrI6uIdhPbWckkrFKSmrg1NS1lMrArcXeWROSXk7y0OCSLSRucGoSYQPemeSOSXlnmweTZ/pbb20nJJKxSkpq4OPYVhPrXB6JZOMKalO2sbhijyckrKefB5GLhlTUp6iNR/3MKwnp6RyyTh9plNS+TSckrJenJLKZfdsaMOtYzrCKal0Gk5JWS8eksolY0rKQ1L5OCVlPbUdq01lbAKlRB9mp6TyafgYhvXilFQuGadobZfDnvJ2lka2q9W6waiJh6RyyZqS8nBULk5JWU+tdriHkUjWlJTPwcil0YC2U1I2Udsz7qWSNSXlbSyXhjwk1Zek5ZJukrRZ0rk9Hj9N0rXlz+WSjup67BZJ10m6RtLGOuueDb40SC5ZU1IeksolW0qqthn3JDWB84ETgVFgg6R1EXFD12I3Ay+MiHskrQDWAMd1PX5CRNxZV82zyVO05iIp3XzLLc+5kk5DSjXsWWcP41hgc0RsiYgdwEXAyu4FIuLyiLinvHklsKjG+gaq7Rn30sl2QNLR7Xyy9TDqbDAWAlu7bo+W9/VzFvCFrtsBXCrpKkmr+q0kaZWkjZI2bt++fY8Knk0eksqnkezD7CGpfIrzMIZdRXW1DUkBvbbknp9GSSdQNBjHd939/IjYJukg4IuSboyIyx71hBFrKIayGBkZmTOfdqek8mkmGy7wsGc+zUauJF6dPYxRYHHX7UXAtokLSToS+ACwMiLu6twfEdvKf+8A1lIMcaXhlFQ+DeVKSbXbTkll45RUfxuAZZKWSpoPnAKs615A0iHAZ4DXRMS3u+7fT9L+nd+BlwLX11b5LPCQVD7FlUTzfJhb4SGpbLINe9Y2JBURuySdA1wCNIELImKTpLPLx1cDfwo8BXhfeXmDXRExAjwVWFveNw+4MCIurqv22dAOT2yTTbbLNrSdkkon27BnnccwiIj1wPoJ963u+v11wOt6rLcFOGri/Zm0vfeXTraUVMspqXSckrKePCSVT7bpM52Syqc4D2PYVVTnBqMmbaek0snWw3BKKp9mI9fVBNxg1KTllFQ6TknZoDklZT15SCqfRiPXAclip8TbWCbZhj3dYNQknJJKJ9sBSc+5kk+2YU83GDVxRj6fbB9mp6TyaSY718cNRk18JdF8sg0X+PIz+RTHMIZdRXVuMGrQGQd3DyOXbD0MX602n70+JVVepqM5iGL2Vp1xcH+Wc5FTUjZge11KSlJD0m9I+ndJdwA3ArdJ2iTpXZKWDb7M3Dp7EB6SyqXplJQNWLZhzyo9jK8AhwFvBp4WEYsj4iDg5ykmOXqHpFcPsMb0OtuDhwtycUrKBi3bsGeVa0m9JCJ2TrwzIu4GPg18WtI+s17ZXqTzpeNjGLk0kn2YnZLKp+hhQESgBN8PU/YwJjYWvY5h9GpQbDcPSeWU8mq1Cb50bLfOTmSW/RIfw6jB7pTUkAuxack2XOAeRj7N8hs4y3bmYxg1GEtJ+cOcilNSNmga62Hk2NAqH8OQdGhEjF2I18cwquv0MDxckEuzIXbsynPt6bZTUul0eoRZGozpHMNYO/ExSc+ZsMykJC2XdJOkzZLO7fH4aZKuLX8ul3RU1XXnsrZTUillS0n5Apf5dI5h7DVDUpJeKekdwP6SnjHhgPeaqi9Urnc+sAI4AjhV0hETFrsZeGFEHAn8Zef5K647ZzkllVMj2fSZLc+HkU7n75VlEqUqxzC+DtwAPBk4D/iOpKslfR54aBqvdSywOSK2RMQO4CJgZfcCEXF5RNxT3rwSWFR13bms7ZRUStl6GG3PuJdOJwiTZTub8hhGRHwf+Iik70bE1wEkHQAspUhMVbUQ2Np1exQ4bpLlzwK+MN11Ja0CVgEccsgh0yhvcDrdzaav3JVKcR7GsKuozimpfDp/ryxDUlM2GJIUha937isPeN89cZmpnqrHfT3XkXQCRYNx/HTXjYg1lENZIyMjc+KvsPtaUv4wZ9JQcUJVFu12keyyPDopqSzbWaVYraQ3SBq3uy5pvqQXSfowcEaF5xkFFnfdXgRsm7iQpCOBDwArI+Ku6aw7VzkllVO2uQrannMlnbEexl7UYCwHWsDHJW2TdIOkLcB3gFOBv4+ID1V4ng3AMklLJc0HTgHWdS9QNkqfAV4TEd+ezrpz2dhBbw8XpNJIdgzDKal8sqWkqhzDeBh4n6T3l8sfCDwUEfdO54UiYpekc4BLgCZwQURsknR2+fhq4E+Bp5SvB7ArIkb6rTud1x+mTgLCPYxcmslSUm2npNLJlpKqcuJex+UR8VzgNgBJTwCWRcR/V32CiFgPrJ9w3+qu318HvK7qulm03cNIKVtKquWUVDpjlwZJsp1NJ7fzOABJ5wFExP3A+wZR1N7GKamcivMwhl1FdR6SyqeRbEhqOl9hknQQ8Grtvg7vjw2gpr1OZ+8hw+WLbbeG8lyyAYp5V7yJ5dJIlpKazpDUm4GvARcCfy/p23hO8Eo8p3dO2VJSLaek0smWkqrcYETExcBPAUh6LvDrFOdK2BR2D0n5w5xJtukzPSSVT7Yhqcon7nXfFxFXAFdMtozt1tkWnJLKJdt8GE5J5dNMlpKq88S9xyynpHJKNyTllFQ62VJSVYaklgOvpThxbylwL7AvxfkQl1KcuHfNoArcGzgllVNDSjOBUkTQDl/gMpu9bkiqc+Iexcl0+zDDE/cey5ySyilTSirGhj2HW4dNz16XkirPu7i2/NkUEbcNvKq9jFNSOWUakvKcKzntdVerBTYDzwF+C3iGpNvZ3YBsAC6LiEcGV2J+TknllCkl1fKcKymNDUkl2c6qDEmNO5u7PI7xs8CRwO8A/yTpdyLiksGUmJ9TUjllSkk5WJFTtpTUdE7cAyAibqaYSnUdgKSnA5+nuDCg9eAPc05FD6MYX57rx59aHvZMKVtKao9zO+UxjQtnoZa9llNSOXW+fDN0MsauiOydklQ6ow5Zroo8K19hEfF/Z+N59lZtp6RS6nz3ZjiO0anR7UUuYw1Ggm0MfC2oWni4IKdGogSLJ+nKKVtKyg1GDZySymnsgGSCvT9PA5yTexiTkLRc0k2SNks6t8fjh0u6QtIjkt404bFbJF0n6RpJG+ures+NnVTlBiOVTNNnuoeR0+4expALqWjaKamZktQEzgdOBEaBDZLWRcQNXYvdDbwROKnP05wQEXcOtNAB8ElVOWWaPtPDnjk95lJS03AssDkitkTEDuAiYGX3AhFxR0RsAHbWWNfA7T6pasiF2LQ0y+/eDB9mp6RyekympCpaCGztuj1a3ldVAJdKukrSqn4LSVolaaOkjdu3b59hqbNrd4LFH+ZMGpmOYTgllZKPYfTXa1Oezrv0/Ig4BlgBvF7SC3otFBFrImIkIkYWLFgwkzpnnYcLcsq09+djGDk5JdXfKLC46/YiYFvVlSNiW/nvHcBaiiGuFHydn5wyTZ/plFROmXqxUG+DsQFYJmmppPnAKZSXF5mKpP0k7d/5HXgpcP3AKp1lnW3Be3+5OCVlg7Z7GxtyIRXVlpKKiF2SzqG45lQTuCAiNkk6u3x8taSnARuBJwJtSb8HHEExB8fa8kzpecCF5RzjKTgllVPGlJR7GLk0kqWkamswACJiPbB+wn2ru36/nWKoaqIfAUcNtrrBcUoqp0yRx06j5h5GLs1Ex8nAZ3rXwuPLOWVKsDgllVOmbQzcYNTCQ1I5ZUxJOViRS6brlYEbjFq0nZJKKWNKyjsluWS6Xhm4wahFK8Jjywk1MqWkfIHLlLKlpNxg1KAd3vPLKNP0mS1fTSClThDGPQwb026HE1IJOSVlg5bpXB9wg1GLVju855eQEiVYnJLKySkpe5RWhIekEsqUkXdKKqfdJ4fO/W0M3GDUohiS8gc5m0wXhnNKKq9mQymGPcENRi2cksppLCWV4MPslFReTckpKdutHU6vZJQpJeU5V/JqNHwMw7q02zGWuLE8MqWkWk5JpVX0MOb+NgZuMGrhlFROTklZHRpSim0M3GDUohVuMDLKlJJqOyWVVqOhFNsYuMGoRTEk5Q9yNplSUp4GOC+npGycVnhsOaNMJ1U5JZVXwymp3iQtl3STpM2Szu3x+OGSrpD0iKQ3TWfduawd4bHlhHb3MIZcSAUeksqr2cgx7Ak1NhiSmsD5wAqKaVdPlXTEhMXuBt4I/N0M1p2zPCSVU8qUlIek0mnKQ1K9HAtsjogtEbEDuAhY2b1ARNwRERuAndNddy5zSiqnTkoqEnyYnZLKS05J9bQQ2Np1e7S8b9DrDl3bKamUMl1J1ENSeTWdkuqp15Zc9V2qvK6kVZI2Stq4ffv2ysUNUstDUik5JWV1KFJSw66imjobjFFgcdftRcC22V43ItZExEhEjCxYsGBGhc62VnjPL6OxK4kmGC5oeRrgtBryQe9eNgDLJC2VNB84BVhXw7pDFxE0/TlOJ9P0mZ1GzT3ZfJqNPJcGmVfXC0XELknnAJcATeCCiNgk6ezy8dWSngZsBJ4ItCX9HnBERPyo17p11b6nPCSVU8MpKatBI1FKqrYGAyAi1gPrJ9y3uuv32ymGmyqtm0WrHWOJG8ujkTAl5c0sn4aUYhsDn+ldi7Zn3EspVUrKZ3qnlWlIyg1GDTwklVMjU0oqnJLKquGUlHVzSiqnZqKUVNspqbSaTklZN6ekcsqUkvI0wHl5SMrG8ZBUTp2UVIYeRqvt4aisMqWk3GDUwCmpnBqJJlCKCCekknJKysZxSiqnsSGpBB9m92Lz8pCUjeMPc05jlwZJ8GFueackLaekbJy2U1JpZZk+s90Ob2NJOSVl47TaTkll1UwyfaZTUnl5SMrGaYf3/rJqNPKkpDznSk4NT6Bk3dqecS+thnJMbhOeNz4tNxg2jg9I5pVlvmUHK/LykJSN02r7oHdWjSTTZ7Y8DXBajYZIsIkBbjBq0Y6g6Xc6pUwpKfcwcmoqxwUuwQ1GLYqUlD/MGTXSpKR8afOsGh6Ssm5OSeXVbOTIyBfBimFXYTPR9EHv3iQtl3STpM2Szu3xuCS9p3z8WknHdD12i6TrJF0jaWOdde8pp6TyypJgafsYRlpZtjGocYpWSU3gfOBEYBTYIGldRNzQtdgKYFn5cxzw/vLfjhMi4s6aSp41PqkqryxXEnVKKq9iSGrYVVRTZw/jWGBzRGyJiB3ARcDKCcusBD4ShSuBJ0l6eo01DkTbJ1Wl1UySknIPI69mkpNDod4GYyGwtev2aHlf1WUCuFTSVZJW9XsRSaskbZS0cfv27bNQ9p5rOSWVVjPJheHcw8iruPxMgo2MehuMXlvzxHdpsmWeHxHHUAxbvV7SC3q9SESsiYiRiBhZsGDBzKudRU5J5dVIcmE4TwOcV5ZzfaDeBmMUWNx1exGwreoyEdH59w5gLcUQ15zXmRjFH+acspyF2/YFLtPKcjUBqLfB2AAsk7RU0nzgFGDdhGXWAaeXaannAD+MiNsk7SdpfwBJ+wEvBa6vsfYZ63zZeHw5pywJFh/DyKs403vub2NQY0oqInZJOge4BGgCF0TEJklnl4+vBtYDLwM2Aw8Cv1mu/lRgbTnN6Tzgwoi4uK7a90Rnz8HjyzllaTBang8jreICl8OuopraGgyAiFhP0Sh037e66/cAXt9jvS3AUQMvcAA6G4L3/nJKMyQVwbyGkxUZNRs5pgEGn+k9cLt7GEMuxGYky/SZTknl5ZSUjfExjNyyTJ/plFRemeaOd4MxYOFjGKmlGZJySiqtTuQ+w7CUG4wBcw8jNyU56O2UVF5jPYwE25kbjAFr+TyM1LJcSdQpqbw6DX2GpJQbjAHrbAQ+0zunNENSngY4rU4gxkNS5pRUck5J2aB1ehgZdkz8NTZgbR/DSC1LSqrtlFRaTaekrKPtlFRqWYakWk5JpdX5bvCQlDkllZxTUjZoklNSVmo7JZValpRU2ymptJpOSVlHyymp1NIMSTkllZZTUjam82XjlFROxaWnh13F1Fpt92Kz2n0extzf0Pw1NmBjQ1Le+0upqRxxx7anAU5r7KB3gu3Mm9iA7e5huMHIqJFlSMrTAKfllJSNcQ8jt4Y0dgHJuawdMZa2sVw6f7cM21mtDYak5ZJukrRZ0rk9Hpek95SPXyvpmKrrzlVOSeWWZb7lts/0TmvsarVOSe0mqQmcD6wAjgBOlXTEhMVWAMvKn1XA+6ex7pzklFRuxZDUsKuYWivcYGQ1lpJKMPRZ5xStxwKby+lWkXQRsBK4oWuZlcBHyqlar5T0JElPB5ZUWHfW/Mo/fo2Hd7Zm5bke3FE8j2fPzKnZgLsfeIQTz/v/wy5lUg/vbHvYM6nO3+23P7qRfec1Z+U5n/z4+Xzy7OfOynN1q7PBWAhs7bo9ChxXYZmFFdcFQNIqit4JhxxyyIwKPWzBfuyYxd3K5x32FH524Y/P2vNZfU46eiH3PLCTYG7v/f3U0/bnl3726cMuw2bg2Yc+mV89ZuGs7aQCPHHffWbtubrV2WD02v2Z+Cnst0yVdYs7I9YAawBGRkZm9Cl/9ynPmslqthcaWXIAI0sOGHYZthd7yhMex3mvPHrYZVRSZ4MxCizuur0I2FZxmfkV1jUzswGqc2R9A7BM0lJJ84FTgHUTllkHnF6mpZ4D/DAibqu4rpmZDVBtPYyI2CXpHOASoAlcEBGbJJ1dPr4aWA+8DNgMPAj85mTr1lW7mZmBMpwsMlMjIyOxcePGYZdhZpaGpKsiYqTXYw57mplZJW4wzMysEjcYZmZWiRsMMzOrZK8+6C1pO3DrDFc/ELhzFssZtEz1ZqoVXO+gud7BmUmth0bEgl4P7NUNxp6QtLFfUmAuylRvplrB9Q6a6x2c2a7VQ1JmZlaJGwwzM6vEDUZ/a4ZdwDRlqjdTreB6B831Ds6s1upjGGZmVol7GGZmVokbDDMzq8QNxgSSlku6SdJmSecOu57JSFos6SuSviVpk6TfHXZNVUhqSvpvSZ8fdi1TKacJ/pSkG8v3efbnvZwlkn6/3A6ul/RxSfsOu6aJJF0g6Q5J13fdd4CkL0r6Tvnvk4dZY0efWt9VbgvXSlor6UlDLHGcXvV2PfYmSSHpwD15DTcYXSQ1gfOBFcARwKmSjhhuVZPaBfxBRDwDeA7w+jleb8fvAt8adhEV/QNwcUQcDhzFHK1b0kLgjcBIRDyTYhqAU4ZbVU8fApZPuO9c4MsRsQz4cnl7LvgQj671i8AzI+JI4NvAm+suahIf4tH1ImkxcCLwvT19ATcY4x0LbI6ILRGxA7gIWDnkmvqKiNsi4ury9/sovswWDreqyUlaBPwS8IFh1zIVSU8EXgB8ECAidkTEvUMtanLzgB+TNA94PHNwVsqIuAy4e8LdK4EPl79/GDipzpr66VVrRFwaEbvKm1dSzP45J/R5bwH+HvhD+kxrPR1uMMZbCGztuj3KHP8C7pC0BHgW8I0hlzKVd1NsvO0h11HFTwDbgX8ph9A+IGm/YRfVS0R8H/g7ir3I2yhmq7x0uFVV9tRyZk3Kfw8acj1VvRb4wrCLmIyklwPfj4hvzsbzucEYTz3um/O5Y0lPAD4N/F5E/GjY9fQj6ZeBOyLiqmHXUtE84Bjg/RHxLOAB5s5wyTjluP9KYClwMLCfpFcPt6q9l6S3UAwJf2zYtfQj6fHAW4A/na3ndIMx3iiwuOv2IuZgt76bpH0oGouPRcRnhl3PFJ4PvFzSLRTDfS+S9NHhljSpUWA0Ijq9tk9RNCBz0UuAmyNie0TsBD4DPG/INVX1A0lPByj/vWPI9UxK0hnALwOnxdw+ke0wih2Ib5afuUXA1ZKeNtMndIMx3gZgmaSlkuZTHDRcN+Sa+pIkivH1b0XEecOuZyoR8eaIWBQRSyje2/+IiDm7FxwRtwNbJf10edeLgRuGWNJkvgc8R9Ljy+3ixczRA/Q9rAPOKH8/A/h/Q6xlUpKWA38EvDwiHhx2PZOJiOsi4qCIWFJ+5kaBY8rtekbcYHQpD2adA1xC8WH7ZERsGm5Vk3o+8BqKPfVryp+XDbuovcwbgI9JuhY4Gnj7cMvprewFfQq4GriO4rM95y5hIenjwBXAT0salXQW8A7gREnfoUjzvGOYNXb0qfW9wP7AF8vP2+qhFtmlT72z+xpzu0dlZmZzhXsYZmZWiRsMMzOrxA2GmZlV4gbDzMwqcYNhZmaVuMEw66G8Su3/6rp9sKRPDei19pGU5ex3ewxzg2HW25OAsQYjIrZFxMkDeq3jgcsH9Nxms8YNhllv7wAOK0/OepekJZ15BiSdKemzkj4n6WZJ50j63+UFCq+UdEC53GGSLpZ0laT/lHR4n9daTo+L2Em6X9I7y/W/JOlYSV+VtKW8qBySfkbSf5V1Xitp2YDeDzM3GGZ9nAt8NyKOjoj/0+PxZwK/QXFJ/L8GHiwvUHgFcHq5zBrgDRHxbOBNwPv6vNYJwFd73L8f8NVy/fuAv6I4E/oVwF+Uy5wN/ENEHA2MUFz+wWwg5g27ALOkvlLOQXKfpB8Cnyvvvw44sryC8POAfysu7QTA4yY+iaSDgbv7XJdoB3Bx1/M+EhE7JV0HLCnvvwJ4SznPyGci4jt7/l8z6809DLOZeaTr93bX7TbFjlgDuLfsoXR+ntHjeVZQXLusl51dV0Mde42I6LwGEXEh8HLgIeASSS/ag/+T2aTcYJj1dh/FReZmpJyX5GZJvw7FlYUlHdVj0Z7HL6qS9BPAloh4D8VVX4+c6XOZTcUNhlkPEXEX8HVJ10t61wyf5jTgLEnfBDYxYbrfcg75ZRFx4x6U+irgeknXAIcDH9mD5zKblK9WazYkko4HXh0RZw+7FrMq3GCYmVklHpIyM7NK3GCYmVklbjDMzKwSNxhmZlaJGwwzM6vEDYaZmVXyP73vy7wx/o/YAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def echo_pulse(tau):\n",
    "    \n",
    "    n_steps_half_tau = int(tau / 2 / time_step)\n",
    "    n_steps_pulse = 4 * n_steps_X_pi_half + 2 * n_steps_half_tau\n",
    "\n",
    "    ctrl_parameter = np.zeros((n_steps_pulse,1)) + eps_idle\n",
    "\n",
    "    ind_start = n_steps_X_pi_half\n",
    "    ind_end = ind_start + n_steps_half_tau\n",
    "    ctrl_parameter[ind_start:ind_end] = eps_work\n",
    "    \n",
    "    ind_start = 3 * n_steps_X_pi_half + n_steps_half_tau\n",
    "    ind_end = ind_start + n_steps_half_tau\n",
    "    ctrl_parameter[ind_start:ind_end] = eps_work\n",
    "\n",
    "    return ctrl_parameter\n",
    "\n",
    "echo_sample_solver = create_monte_carlo_solver(10, echo_pulse)\n",
    "echo_sample_amplitudes = echo_sample_solver._ctrl_amps\n",
    "plt.plot(np.arange(len(echo_sample_amplitudes)) * time_step,\n",
    "         echo_sample_amplitudes)\n",
    "plt.xlabel(\"time / ms\")\n",
    "plt.ylabel(\"$J(\\\\varepsilon(t))$\")\n",
    "plt.title(\"Spin Echo Control Amplitude\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [],
   "source": [
    "precession_times_2 = np.arange(0,100,3)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [],
   "source": [
    "probabilities_echo = np.zeros(len(precession_times_2))\n",
    "\n",
    "measurement_state_echo = initial_state\n",
    "\n",
    "for i, tau in enumerate(precession_times_2):\n",
    "    \n",
    "    solver = create_monte_carlo_solver(tau, echo_pulse)\n",
    "    probabilities_echo[i] = average_monte_carlo_traces(solver, initial_state,\n",
    "                                                       measurement_state_echo)\n",
    "    "
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "outputs": [],
   "source": [
    "# np.save(r\"..\\stored\\precession_times_2\", precession_times_2,)\n",
    "# np.save(r\"..\\stored\\probabilities_echo_2\", probabilities_echo,)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "outputs": [],
   "source": [
    "precession_times_2_loaded = np.load(r\"..\\stored\\precession_times_2.npy\")\n",
    "probabilities_echo_loaded = np.load(r\"..\\stored\\probabilities_echo_2.npy\")\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit parameters:  [135.30109118   1.90271031]\n"
     ]
    }
   ],
   "source": [
    "def T2_echo_decay(t, T2_echo, alpha):\n",
    "    return 0.5 * (1 + np.exp(-(t / T2_echo) ** alpha))\n",
    "    \n",
    "popt_echo, cov_echo_al = curve_fit(T2_echo_decay, precession_times_2_loaded,\n",
    "                                   probabilities_echo_loaded, p0=[136, 2])\n",
    "print(\"Fit parameters: \", popt_echo)\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEJCAYAAABlmAtYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0m0lEQVR4nO3deVxU1fvA8c8zw6AoKu4JWGiuKbkvSRrumppLlmnmll/LtM3y59K+U7Zpu5XaZmbmVpZbaqVWmmmpqYlLCZpLhivKAOf3xwyICM4gs8DwvF8vXzL3nnvnOYjzcJZ7jhhjUEoppS7G4u8AlFJKFXyaLJRSSrmkyUIppZRLmiyUUkq5pMlCKaWUS0H+DsBbKlSoYKKiovwdhlJKFRobNmw4YoypmNO5gE0WUVFR/PLLL/4OQymlCg0R+Su3c9oNpZRSyiVNFkoppVzSZKGUUsqlgB2zUKowstvtJCQkcObMGX+HogJY8eLFiYyMxGazuX2NJgulCpCEhARKlSpFVFQUIuLvcFQAMsbw77//kpCQQLVq1dy+zifdUCIyTUQOiciWXM6LiEwRkXgR+V1EGmc510VEdjjPjfdmnPM3JhITt4Jq4xcRE7eC+RsTvfl2Sl3gzJkzlC9fXhOF8hoRoXz58nluvfqqZTEDeB34MJfzXYGazj8tgLeAFiJiBd4AOgIJwHoRWWiM+cPTAc7fmMiEuZtJtqcBkJiUzIS5mwHo1SjivHKTluxgf1Iy4WEhjO1c+7zzSuWXJgrlbZfyM+aTZGGM+V5Eoi5SpCfwoXGsl/6TiISJSBUgCog3xuwGEJFZzrIeTxaTluwg2Z7G3da5CJCKlZT0IPYsWgZp9SCoOD/vtzP3xyOUSS2OSAmOJ5VgwtyzgCYUpVRgKyhjFhHAviyvE5zHcjreIrebiMgIYATA5ZdfnqcA9iclA3Bn0JeUlLPnTqQCXzm+bAG0sALWc6fPmiD+W1AWNkRBaGV2ny3Nnl3pNEstx99Sib1Jl2lCUYXKM888w8yZM7FarVgsFt555x1atGjB8OHDGTNmDFdddVW+3yPjodkKFSrkWubZZ59l4sSJma9btWrF2rVr8/3eQ4YM4bvvvqN06dIkJyfTsmVLnnvuOSIiLv7/79VXX2XEiBGUKFEi3zEURgUlWeTUJjIXOZ4jY8xUYCpA06ZN87SrU3hYCIlJydQ7Ox0L6dhIxUYqVcvY+GZ0S7An023S15SWU5TmdObfFeQYlSSJG4Ot8G88FQ4ncL/lJASfu/cJE8KBhZfBn9FQoRbrk6vw/s+Gg/ZKGIJy7fJSytd+/PFHvvrqK3799VeKFSvGkSNHSElJAeC9997zaSzZk4UnEkWGSZMm0bdvX4wxvPrqq7Rt25YtW7YQHByc6zWvvvoqAwcOLLLJoqA8Z5EAVM3yOhLYf5HjHje2c21CbI4mQzoWzhJMmq0Ud3RpBqUug3LVSCpTlx/T67EkvRmfp8Xyftr1PJ/an5dLjoFBC2DUzzQ4M5XaZ2bQ/uwkhqU8yBP225iT1obE1NJwcCusfoVmvzzAl9YH2VpsGN8Ej+cV2xsMTF/A0m/mwtmT58Wlg+7qYjz983HgwAEqVKhAsWLFAKhQoQLh4eEAxMbGZi6hExoayrhx42jSpAkdOnRg3bp1xMbGUr16dRYuXAjAjBkzGD16dOa9u3fvzqpVqy54z169etGkSRPq1avH1KlTARg/fjzJyck0bNiQW2+9NfM9wTGbZ+zYsdSvX5/o6Gg+++wzAFatWkVsbCx9+/alTp063HrrrbjaCVREuP/++7nsssv45ptvABg5ciRNmzalXr16PPbYYwBMmTKF/fv307ZtW9q2bZtruUBWUFoWC4HRzjGJFsAxY8wBETkM1BSRakAicAswwBsBZPxGf7GuobGda583CA4QYrMytnPtzNeOFgrsMhHsMueujQgLYc097cB+husffZ9akkAdyz5qy9+0sGyjt3UNpABxj0LFuhDRmE3p1Znxawn22yMwWLQFos7j7qSMvOjUqRNPPvkktWrVokOHDvTr14/rrrvugnKnTp0iNjaW559/nt69e/Pwww+zbNky/vjjDwYPHswNN9zg9ntOmzaNcuXKkZycTLNmzbjxxhuJi4vj9ddfZ9OmTReUnzt3Lps2beK3337jyJEjNGvWjDZt2gCwceNGtm7dSnh4ODExMaxZs4Zrr73WZQyNGzdm+/bt9OzZk2eeeYZy5cqRlpZG+/bt+f3337nnnnt4+eWXWblyZWbXWU7lrr76arfrXdj4JFmIyKdALFBBRBKAxwAbgDHmbeBr4HogHjgNDHWeSxWR0cASHCMF04wxW70VZ69GERf9T+aRhGIrzrEydZmfFAXp5+5djuO0K5XAi61SIXEDbF9Ew+SjzLfCMUsJ1qXX4ef0uvycWpeXF1t1/ENlTsrIKtmexqQlOy753z80NJQNGzbwww8/sHLlSvr160dcXBxDhgw5r1xwcDBdunQBIDo6mmLFimGz2YiOjmbv3r15es8pU6Ywb948APbt28fOnTspX758ruVXr15N//79sVqtVK5cmeuuu47169dTunRpmjdvTmRkJAANGzZk7969biWLrC2Q2bNnM3XqVFJTUzlw4AB//PFHjknA3XKBwlezofq7OG+AUbmc+xpHMikQvJVQkm1lufb6NpBRzhjaTJxOY9lJC8s2Wli20dH2KwAnzoTAJ63hyvYss0czYelxku2OzKOtj6IjY1KGu8fdZbVaiY2NJTY2lujoaD744IMLkoXNZsucfmmxWDK7rSwWC6mpqQAEBQWRnn7uN6Kc5vWvWrWK5cuX8+OPP1KiRAliY2Ndzv+/WNdSRhwZ9ciIxZWNGzfSvn179uzZw4svvsj69espW7YsQ4YMyTEed8sFkoLSDRVQPJFQECGtTBTzkyozP93xm1El/qOFZRvtQnbS+9942LmUjsASqch3QQ34Lr0Ba9PrcdpePF+/XarCIWNSRk7HL9WOHTuwWCzUrFkTgE2bNnHFFVdc0r2ioqJ48803SU9PJzExkXXr1l1Q5tixY5QtW5YSJUqwfft2fvrpp8xzNpsNu91+wZIUbdq04Z133mHw4MEcPXqU77//nkmTJrF9+/Y8x2iM4bXXXuPAgQN06dKFbdu2UbJkScqUKcPBgwf55ptviI2NBaBUqVKcOHGCChUqcPz48VzLBSpNFn7iKqHAhS2QQ5RlubU17bvf5WiBHN3Dwy9N5jrL7/Sx/sBtQctJMVZ+Tq/L0hPN4HhdKF0F0K6qQOTOGFpenTx5krvvvpukpCSCgoKoUaNG5qBzXsXExFCtWjWio6OpX78+jRs3vqBMly5dePvtt7n66qupXbs2LVu2zDw3YsQIrr76aho3bswnn3ySebx37978+OOPNGjQABHhhRde4LLLLstTshg7dixPPfUUp0+fpmXLlqxcuZLg4GAaNGhAo0aNqFevHtWrVycmJua8eLp27UqVKlVYuXJlruUClbiaLVBYNW3a1ATC5keuPuRj4laQmJRMMHaaWP7kOstvdLRs4ErLAUeByGZsKd2GBzdXZbu9UuZ1ITYrz/WJ1oRRwGzbto26deu6XV5/CVCXKqefNRHZYIxpmlN5TRaFXPYZMQAhNguvdyhBe1kH276EA78BsD29KovSWrAgPYa/TWXHDK3x7c67l37w+Fdek4VSlyqvyUK7oQq53MY/2jeKAK6DNmO5dvwMOlp/oYt1PQ/Y5vAAc9iYXoOFJ1rByfoQWskr0zCVUoFDk0UAcDX+YcIuZ3pSRaandaUK/9LDupae1rU8ZvsQXvoYqsey+e/6YG8IFM+8Lr/TMJVSgaOgPMGtvCjr0+kHKM/UtB70NS/wbduFcO0Y+HcXj9insK7YKJ4NepeGEk/Gqir5nYaplAoM2rIoAlx2VbV7mJHPvUG700voZV3LgKCV7EiPZHZaLD+Gdrzgfjq2oVTRo8miiLhoV5UInbv2ZsLcGjxxdhA9rD/Sz7qKR2wfk54yCz6fD83+B1e0Yv6m/Tq2oVQRpMlCAee3PmYltef7Ut156hqh3ekl8NtM2DoPKtVjx9E2YG+Gjm0oVbTomIXK1KtRBGvGt2NPXDfWjG9Hu+vaQtc4GLMdekwBsTAu9S1+Ljaah4M+4gr5J/Pa7GMbulpu4Wa1WmnYsGHmn71799KqVSsAkpKSePPNN13e491338283mKxZH49ZswYb4efb8OGDaNSpUrUr1//vONnzpyhefPmNGjQ4ILVZqOiooiOjqZhw4Y0bZrj7NOLXr948WJq165NjRo1iIuL807F8kGfs1DuM4Y7n3uT65O/pKtlHTZJ49u0Rryb1o19pRqzZkJ7ILdnP/QhQHcUlOcsQkNDOXnyZI7n9u7dS/fu3dmyZYtb90pMTKRVq1b89ddfngzRq77//ntCQ0MZNGjQefU0xnDq1ClCQ0Ox2+1ce+21TJ48mZYtW7q1oVNu1zdr1oxatWqxbNkyIiMjadasGZ9++qlHNprKTV6fs9CWhXKfCF269mIc99Hq7BReTe1DA8suZgU/zVchj8KWLyAt9aKroarCK2M/ifHjx7Nr1y4aNmzI2LFjXV63ZcsWoqOjvR2eR7Vp04Zy5cpdcFxEMr8Pdrsdu92ep/2sc7t+3bp11KhRg+rVqxMcHMwtt9zCggULLrh+79691KlTh+HDh1O/fn1uvfVWli9fTkxMDDVr1mTdunWcOnWKbt260aBBA+rXr5+530d+6ZiFypOsYxuTk/qyoOTNvFJnGw33fQxzhkHY5XQ6EctnxHI6y7gG5NxVpbOqCqaMjYcAqlWrlrmEOEBcXBxbtmzJca+JnGzevPmC7hx/aN26NSdOnLjg+IsvvkiHDh3cvk9aWhpNmjQhPj6eUaNG0aKFY6dnEaFTp06ICHfccQcjRoxw+/o5c+ZQteq5fd4iIyP5+eefc7w+Pj6ezz//nKlTp9KsWTNmzpzJ6tWrWbhwIc8++yy33XYb4eHhLFq0CHAs1ugJmixUnl04s6orpN8HO76Gta/xWNKH3Bc0hw/TOjEttQv/URo4fzVUfWLcDd+Mh382e/ael0U7xqFcCAkJcTsZuLJlyxY6djw3BXv+/PksWrSIQ4cOMWrUKDp16nRJ9x0+fHietnr94YcfLul9srNarWzatImkpCR69+7Nli1bqF+/PmvWrCE8PJxDhw7RsWNH6tSpk7kpk6vrcxoOyK3FkrE4I0C9evVo3749IpK5l0h0dDQPPvgg48aNo3v37rRu3doj9dZuKOUZFgvU7Q63L+G71jNZRz1GWRewuti9TAj6hKq2E+ethqpdVUVH9pZFr169ePfdd5kxY0aOXSR2u52JEydy3333cddddwGQmprKmDFjeOCBB5gyZQqnT58mPj6ehx56iN69e5OcnMyYMWMYNWoUDz/8cI5xtG7d+rxB+4w/y5cvv6R6hYWFERsby+LFiwEyt5+tVKkSvXv3znFJ9tyuj4yMZN++fZnnEhISMu+XXdY9O3LaS6RWrVps2LCB6OhoJkyYwJNPPnlJ9ctOWxbK465r34355Rpy2zfLuSl5FsODvmG4ZTnWf4bC8XuhdLjXNu4JKG60APwhY1+HrNq3b8+HH35IRMT5rcL09HR27txJnTp1LrjP008/zahRF+55NnXqVJKTkwkLC2P37t0AvPXWW/Ts2TNzi9fVq1fTtWtXxo0bx8CBA3nttdcYMGAATZs2pW/fvjnG7YmWxeHDh7HZbISFhZGcnMzy5csZN24cp06dIj09nVKlSnHq1CmWLl3Ko48+6vb1zZo1Y+fOnezZs4eIiAhmzZrFzJkzLynG/fv3U65cOQYOHEhoaCgzZszIZ60dNFkor3B0VQ0GBsO/u+CHl2Ddu/DLNGh0G43KNOfXY6UuuC4/G/co3yhfvjwxMTHUr1+frl278vzzzxMfH5/jgHB8fDyRkZHn/TZsjGH8+PF07do1xz0uNm7cyBtvvHHeNb/++isjR47MfL1+/Xp69OgBOLp1tm7dyr333ktKSgolSpTIdx379+/PqlWrOHLkCJGRkTzxxBPcfvvtHDhwgMGDB5OWlkZ6ejo333wz3bt3Z/fu3fTu3RtwtIIGDBiQue3s9ddfz3vvvUd4eHiu1wO8/vrrdO7cmbS0NIYNG0a9evUuKfbNmzczduxYLBYLNpuNt956K9/fD9Cps8qX/tsLq1+BjZ+QbgyfpbXllZReHKIsoNNroeBMnc2LLVu2MG3aNF5++WW3yk+ZMoUPPviAZs2a0bBhQ+68804OHjzIV199xe23386XX37JzJkzqVq1Ku3ataNLly4sWLCAhQsXUq5cOSZMmMC4ceN4++23+e+//4iLi6Nt27bMnj2bEiVKcM899xS676E/6H4WTposCrBjCfDDS6Rv+BC7sfBBakfml7yZEV2aFelEAYUzWXjC4sWLCQ4Opl27dq4LK4/Q/SxUwVcmErq/gqXVPRRbFceI3z9jhPkekkbBmVFQvLS/I1Q+ltFlowounQ2l/KdcNejzDtz1I1wZC9/FweQGsGYy2B0D3bpsiFIFg7YslP9Vqgv9PobEX2HF07DsUfj5HTbUuJuJ66/gtN3RVarPYijlP9qyUAVHRGO4bS4M/gpKVqTJr+OZLRO4xrI1s4g+i6GUf2iyUAVPtdbwv5Xcl3IXYXKST4Of4T3bJK4URxeUrnCrlO9pslAFk8XC+tIdaX/2ReLst9Dcsp0lweN4Ouh96pU5m1ksY9mQxKRkDOe6qgpzwgjUGYqq4LiUnzFNFqrAGtu5NhZbCG+n3UDs2Vf4OK0D/ayrmJc6Cn54GVLPBtyyIcWLF+fff//VhKG8xhjDv//+S/HixV0XzkIHuFWBdf7e4fBu6EgiYu6jY8Lr8O0TsPFjah3vQyKNLri2sC4bEhkZSUJCAocPH/Z3KCqAFS9enMjIyDxdo8lCFWg57x0eAzuXw+LxTA+exLdpjXgy9Tb+Mpdlliisy4bYbDaqVavm7zCUuoB2Q6nCqWYHGLmWLfUepKVlG0uD/4//C5pFCc4QYrOet8KtUir/NFmowisomPo3PcL3XZawwnotdwUt5LviD/JR87/o1TDn5Z2VUpdG14ZSgWPfOvh6LBzYBFGtodvLULGWv6NSqtDQPbhV0VC1OfxvBXR/Bf75Hd5qBSueyVw6RCl16TRZqMBisULTYTD6F6jfB75/Ad68Bnat8HdkShVqmixUYAqtBH2mwqAFIBb4qDfMuR1OHPR3ZEoVSjp1VgW26rEwci2sfoW0H17i9Javed7ej1Wh3XmwS11dkFApN2nLQgU+W3Hmhw2iu/0FfkurxtO26bySPIGpcxcX6mVBlPIlTRaqSJi0ZAfb7JUZaJ/ImJQ7qSmJzLP8HwcXPQNpdn+Hp1SB57NkISJdRGSHiMSLyPgczpcVkXki8ruIrBOR+lnO7RWRzSKySUR0PqzKs3PLfwhz09vQ4eyLLE9vzB2pn8C7beHAb4CuYKtUbnwyZiEiVuANoCOQAKwXkYXGmD+yFJsIbDLG9BaROs7y7bOcb2uMOeKLeFXgCQ8LITHLelFHKMMo+32sKb6JZ09+AFPbsqPGMB7b3oZjdiugmy0plZWvWhbNgXhjzG5jTAowC+iZrcxVwLcAxpjtQJSIVPZRfCrAje1cmxCb9bxjITYrza8fAqN+hgb9qb3zXebK/9FEzq1YW5hXsFXKk3yVLCKAfVleJziPZfUb0AdARJoDVwAZyyIaYKmIbBCREbm9iYiMEJFfROQXXbVTZdWrUQTP9YkmIiwEASLCQniuT7SjxRBSFnq9waCU8RQTO58HP8ljQR9QHMe+GYV1BVulPMlXU2clh2PZ1xmJAyaLyCZgM7ARSHWeizHG7BeRSsAyEdlujPn+ghsaMxWYCo7lPjwVvAoMOa9ge86u0i3olFSLsUGfMTRoCddZfuNB+50cLNPAh1EqVTD5qmWRAFTN8joS2J+1gDHmuDFmqDGmITAIqAjscZ7b7/z7EDAPR7eWUh41tnNtjK0kT6QOpn/KQwRLKp8HP8H0qosg9azrGygVwHyVLNYDNUWkmogEA7cAC7MWEJEw5zmA4cD3xpjjIlJSREo5y5QEOgFbfBS3KkKydlX9lF6PocUns++KPtTa+R5MbQsHfvd3iEr5jU+6oYwxqSIyGlgCWIFpxpitInKn8/zbQF3gQxFJA/4AbndeXhmYJyIZ8c40xiz2Rdyq6Lmwq6o7/NkPFt7tmGJ73XgWhN7MC8t2sT8pmfCwEMZ2rq2zpVTA0yXKlXLH6aPw9YOw5Qt+N1dyf8qd7DKOBBFis54bLFeqENMlypXKrxLloO80HrE9SCQHWRQ8kcHWJYDR6bWqSNBkoVQefHyiMZ3PvsDa9Ho8YfuA6bYXqMAxnV6rAp4mC6XyIDwshMOEMcw+loftQ7nG8geLi43jxlJb/R2aUl6lyUKpPDj3JLjwcVpHuqc8w2HK8qL9GVj0oO7KpwKWJgul8iD7k+DJZWqys8cCuGY0rH8X3rlOp9iqgKSzoZTylF0rYN5ISD4K7R+DlneB5eK/j83fmMikJTt0Gq4qEHQ2lFK+cGU7x658NTvB0ofg495w/ECuxedvTGTC3M0kJiVjOLfKrS6LrgoiTRZKeVLJ8tDvY+gxGfatg7djYOfyHItOWrKDZHvaecd0Gq4qqDRZKOVpItBkCIz4DkIvg09uhGWPXrAjX27TbXUariqINFko5S0Va8H/voUmQ2HNZJjeFZL+zjwdHhaS42W5HVfKnzRZKOVNthDo8Sr0nQ6Hd8Db18K2r4DcN2Qa27m2HwJV6uI0WSjlC/X7wB3fQdlq8Nmt8PX/0Su6Qu4bMilVwOjUWaV8KfUsLH8cfnoTqjRwtDjKX+nvqJQCdOqsUj43f2MiMXErqDZ+ETFxK85Nhw0qBl2eg1s+hf/+gnfawOY5/g1WKTdoslDKw9x6fqLO9XDnaqhcH764Hb68D+xn/BWyUi5pslDKw9x+fiKsKgxZBDH3wobpMK2zo7WhVAGkyUIpD8vT8xPWIOj4JNwyE47ucXRL/bnEyxEqlXeaLJTysEt6fqJON7hjlaO1MfNm+PYpSE/LvbxSPqbJQikPu+TnJ8pVh9uXQaPb4IcX4aPecPKwFyNVyn2aLJTysOzLmOfp+QlbCPR8HXq+Aft+dnRL/f2z12NWyhV9zkKpgurA7zB7EOlJ+3jNOohXT7YnPKyELmOuvEafs1CqMKpyNYuu+ZQV6Y24N3Uar9mm8F/Sf7qMufILTRZKFWDPrjjA8LP38ay9P10t65gX/CiVUxN0GXPlc5oslCrAHNNthalpPRhon0hFSWJh8CPUPr7W36GpIkaThVIFWNbptj+m16PH2Wf421TiveAX4bsXID0duMjyIkp5iCYLpQqw7NNwE6nIQPMUCZE9YOUz8NlAvlq3Q7dnVV6nyUKpAiynabiP92nC5bd/CF2ehz8XU/+b3oSn/n3edbo9q/I0nTqrVGG2dzVHpvenGHbG2EeyLP3crEcB9sR1819sqtDRqbNKBaqoa/lf8ZfYZarwbvDL3B/0OYJjHEO3Z1WepMlCqUJucJcYBpsnmJ16HfcGzeM920tUsp3R7VmVRwX5OwClVP5kPM09aXFJfj9ZncdsH/Ft6FOUqqqbKinP0TELpQLNX2vhs4GQlgo3TYca7f0dkSokdMxCqaLkilbwv5WO5c4/6Qs/vgkB+kuh8h3thlIqEJW9AoYtgXl3wJIJcHArCyMf5Pnle9iflEx4WIguSKjyxK2WhYiUFZFBIjJPRLaKyFci8j8RqeztAJVSl6hYKNz8EVw3DjZ9TOSXN3M26R99cE9dEpfJQkTmAnOBSsA4Y0w94C6gOPCRiKzyaoRKqUtnsUDbiTxse5C67GVBsYepJ3sBfXBP5Y07LYthxpi2xpgXjTF/Ahhj/jbGvGaM6QT08mqESql8++REY/qmPA7AnODHud7yE5D7fuFKZecyWRhjkrK+FpGSImLN7XxuRKSLiOwQkXgRGZ/D+bLObq7fRWSdiNR391ql1MWFh4Ww1UTR8+zTbDVRvBk8hfuD5hBRppi/Q1OFhDvdUBYRGSAii0TkELAdOOAcu5gkIjXduIcVeAPoClwF9BeRq7IVmwhsMsZcDQwCJufhWqXURWQsSHiEMgxIecj5AN9cZpd7C86e9Hd4qhBwpxtqJXAlMAG4zBhT1RhTCWgN/ATEichAF/doDsQbY3YbY1KAWUDPbGWuAr4FMMZsB6KcA+juXKuUuoisCxLasTG55L1srj+e8H9WwPQucEwHutXFuTN1toMxxp79oDHmKPAF8IWI2FzcIwLYl+V1AtAiW5nfgD7AahFpDlwBRLp5LQAiMgIYAXD55Ze7CEmpoqVXo4hsU2XbQ4Om8PlQeLcdDJjF/IOVmLRkh06vVRdwZ8zigkRxCWUkp8uyvY4DyorIJuBuYCOQ6ua1GXFMNcY0NcY0rVixoouQlFLU7Ai3LwVrMKnvd2HF3Pd0XwyVI3fGLL4VkXpZXt8gIg+LSI6/3eciAaia5XUksD9rAWPMcWPMUGNMQxxjFhWBPe5cq5TKh8pXwf++ZXv65UyxvsxI60Iyfh/T6bUqgztjFpHGmK0AItIK+Ai4HJguIr3dfJ/1QE0RqSYiwcAtwMKsBUQkzHkOYDjwvTHmuDvXKqXyKbQSfZMnsjDtGsbZZvFC0FRspAI6vVY5uDNmcTzL14OAt40x40SkEo4P7XmubmCMSRWR0cASwApMM8ZsFZE7neffBuoCH4pIGvAHcPvFrnW7hkopt5QPK8M9SaPZY6pwb9BcLrcc4s6U+ygZVsnfoakCwOWqsyLyKY6B7O+B34E+xpi1znMbjTGNvB7lJdBVZ5XKm/kbE5kwdzPJ9jR6Wlbzgm0qB6jIzg7v07F1jL/DUz6Q31Vn7wfuABKBX7MkChsQ6rEolVJ+lXV67cL0a7mn2FNUKXaWjmsGwJ7v/R2e8jO397MQEYsxJj3L605AX2PMCG8Flx/aslDKA/7bCzP7wb/x0P1VaHybvyNSXpSvloWICEDWROF8vTQjUWSUUUoFmLJRjqm11drAwtGw7FFIT3d5mQo8bj3BLSJ3i8h5T7mJSLCItBORD4DB3glPKeV3xcvAgM+h2XBYMxk+HwQpp/0dlfIxd2ZDdQGGAZ+KSHXgPxzLk1uBpcArxphNXotQKeV/1iC4/kUoXwMWT4DjPaD/LAjVh1+Lijztwe0c1K4AJLu72qy/6JiFUl6y7Sv4YjiUqgy3zoEKLtcSVYVEfscsBovIERE5CrwHnCzoiUIp5UV1u8OQrxyr1b7fEf5a6++IlA+4M2bxCNARqAP8DTzr1YiUUgVfZFMYvhxKVIAPe8KWL/wdkfIyd5LFcWPMRmPMIWPMIziWDFdKFXXlqjlmSkU0hTnDYPUrkIdubVW4uJMsqojICBFpLSIVAVfLkSuliooS5eC2eVD/Rlj+OHx1P6Sl+jsq5QXuzIZ6HLgauBWIBkJF5Gsc+0/8boz51HvhKaUKovkbE8/f96LTE/QKu9zRuji2D26aAcVK+TtM5UHutCzeNcaMNsZcZ4wpB1QHXgeSgOv1gTylipaMNaTO2/di3lbml/+f4ynvXSthelc4fsDfoSoPcidZLBORz0Skv4iUNsYkAKuAeOf5X70WnVKqwJm0ZAfJ9rTzjmXue9F0KAz4DI7ugffaw0FdIDpQuLNTXnvgCSAKWCQiP+HYKzsaxwN5BXLVWaWUd+S2v0Xm8ZodYeg3YNLh/c6wa4UPo1Pe4k7LAmPMH8aY54wxrYFYY8w1xpjHjTHaqlCqiAkPC3F9vMrVjqm1YVXhk5tg4yc+ik55i1vJIitjzBlvBKKUKhzGdq5NiM163rEQm5WxnWufd2z+bqHTsYmstteGBXfxx6xHLmlq7fyNicTEraDa+EXExK3QPcH9JM/JQilVtGXd90KAiLAQnusTTa9GEZllMgbB/zxmYah9HPPSYrhq+xR2f3AnpKflfvNschxMn7tZE4YfuDN1VimlztOrUcR5ySG7rIPgdoIYYx/JQVOOO/fOgtkn4cb3wJZzd1Zu98mQMZh+sfdXnqctC6WUx2UfBDdYiEvtzxP2QbB9EXxwA5w+muf7uDquvEeThVLK43IbBF9aqrfjgb0Dv8H7neC/vy7pPrkdV96jyUIp5XEXHQSv18uxRMipQ5x5ux1Dnn0/18FrdwfTlfdpslBKeZzLQfCoGL695kP+O2N47exDtLJsznHw2p3BdOUbedr8qDDRzY+UKthi4laQlpTAjOAXuFL286D9DhakX0tEWAhrxrfzd3hFUr42P1JKKW/Yn5TMP5Tn5pRH+SW9NpOD3+QO65fsT9L9vQsiTRZKKb/IGKQ+TkkG28fxZVpLJtg+5YWSM/P0LIbyDU0WSim/yDp4nYKNe+yjmZHejZvSFsGcoWDXxSIKEn0oTynlFxmD1Of2xShJWOdJkNwKlj4EJw9D/5kQUtbPkSrQAW6lVEG0eQ7MHwnlqsPAL6BMpL8jKhJ0gFspVbhE93UkiWOJjmXOD+/wd0RFniYLpVTBVK0NDF0EaSkwrTPsW+/viIo0TRZKqYKrSgO4fQkUD4MPesDOZf6OqMjSZKGUKtjKVYfbl0KFmvDpLfDbLH9HVCRpslBKFXyhlWDIIriiFcy7A9a+5u+IihxNFkqpwqF4abh1DlzVE5Y+DEsvbec9dWn0OQulVOERVAz6Toevx8LaKXDqCNwwBaw25m9MzPLMRghjO9fWBQc9SJOFUqpwsVih20sQWhlWPQun/+XLWs8yYWF85q56GSvYApowPES7oZRShY8IxI6Dbi/DzqVEfT2AYPux84pkbL+qPEOThVKq8Gp2O9z8AbXSd/F58BNcxr/nndbtVz3HZ8lCRLqIyA4RiReR8TmcLyMiX4rIbyKyVUSGZjm3V0Q2i8gmEdE1PJRS51zVkweDH6WKHOWLYo9zpZzbPEm3X/UcnyQLEbECbwBdgauA/iJyVbZio4A/jDENgFjgJREJznK+rTGmYW7rliiliq721/dlcPpjBJPKnOAnaCQ7dftVD/NVy6I5EG+M2W2MSQFmAT2zlTFAKRERIBQ4CqT6KD6lVCHWq1EEg/rcwKhiz3HMlGRmsWd5P+Y/Hdz2IF8liwhgX5bXCc5jWb0O1AX2A5uBe40x6c5zBlgqIhtEZERubyIiI0TkFxH55fDhw56LXilV4PVqFMHsibcS9X+rCbmsNq3WjYbfPvN3WAHDV8lCcjiW/WmazsAmIBxoCLwuIqWd52KMMY1xdGONEpE2Ob2JMWaqMaapMaZpxYoVPRK4UqqQOe9p7xGw9nV/RxQQfJUsEoCqWV5H4mhBZDUUmGsc4oE9QB0AY8x+59+HgHk4urWUUipn5z3t/RAsf1yf9s4nXyWL9UBNEanmHLS+BViYrczfQHsAEakM1AZ2i0hJESnlPF4S6ARs8VHcSqnCKuNp7yZDYfUr8OW9urd3PvjkCW5jTKqIjAaWAFZgmjFmq4jc6Tz/NvAUMENENuPothpnjDkiItWBeY5xb4KAmcaYxb6IWylVyFms0P0VKFEOfngJziRBn3cdiUTlic+W+zDGfA18ne3Y21m+3o+j1ZD9ut1AA68HqJQKTCLQ/lEIKefokjpzDPp9AsVC/R1ZoaJPcCulioZWo6HXW7DnB/jwBjh91N8RFSqaLJRSRUfDAdDvY/hnC0zr4tjjW7lFk4VSqmipcz3cNheO73fs7X0k3t8RFQqaLJRSRU/UtTDkK7AnOxLG/k3+jqjA02ShlCqawhvCsCVgC4EZ3WHvan9HVKBpslBKFV0VajgSRulw+KgPbP/a9TVFlO6Up5Qq2spEwNBv4JO+8NlANjR8mnu21dHtWbPRloVSSpUsD4MXcqh8M5psnECXE19gOLc96/yNOmtKk4VSSgEUK8XNJ8bwdVpzHrF9zANBswGj27M6abJQSimnv46lMtp+D5+mtuXuoPk8HTQNC+m6PSs6ZqGUUpnCw0IcXU+pw/mPUtwVtJAycooXS4zxd2h+py0LpZRyGtu5NiE2KyC8kHoLz9gH0MP6E7NLT4aUU/4Oz680WSillFOvRhE81yeaiLAQBPi61E382vBpKh/5ET7sWaTXkxIToBuCNG3a1Pzyyy/+DkMpVYDM35jIpCU78j4tdtuXMGcYlK8BA+dC6SreD9YPRGSDMaZpTud0zEIpVSTM35jIhLmbSbY7NkDKmBYLuE4YdXvAwC/g0/4wrRPLmrzD42vOFKlnMbQbSilVJExasiMzUWTI07TYam1g8JecPX2ChstvofSx7UXqWQxNFkqpIiG36a95mhYb0ZghPEkKQXwW/BRNZTuQx6RTSGmyUEoVCeFhIXk6npufjpen79nHOWzK8FFwHLGWjUAek04hpMlCKVUknJsWe06IzcrYzrXzdJ/wsBAOUJ6bUh4j3oTzru1leljW5jnpFDaaLJRSRUL2abERYSE81yc6zwPTGUnnKKXpn/IwG0wtJtve4I1aG70TeAGhU2eVUiqPsk7BjSpjYWbZd6jyz0po9wi0fgBE/B3iJdGps0op5UG9GkWc3yJJ6wALRsGKp+BMEnR8qtAmjNxoslBKqfyy2qDX21C8DKx9DZL/g+6TwRo4H7GBUxOllPIniwW6vgAhZeG75+HMMbjxfQgq5u/IPEIHuJVSylNEoO1E6BLnWCJk5s1w9qS/o/IIbVkopZSntRzp6JJaMNqxAOGtn0OJcjkWveT1qnxMWxZKKeUNDQdAv4/gn80woxscP3BBkYz1qhKTkgv80iGaLJRSylvqdIOBcyDpb5jWGY7uPu90vter8iFNFkop5U3V2sDghXD2BEzrAge3Zp7yyHpVPqLJQimlvC2iCQz9BsQK07vCvnWA59ar8gVNFkop5QuV6sCwxVCivGPQO/5bj61X5QuaLJRSylfKXgHDlkC5K2FmP3oFr/fIelW+oGtDKaWUryUnwcx+kLAOekyGxoP8HRFw8bWhtGWhlFK+FhIGt82DK9vDwrthzWR/R+SSJgullPKH4BJwy0yo1weWPQrLH4cC3NOjT3ArpZS/BAXDje85WhqrX3F0T3V7CSxWV1f6nCYLpZTyJ4sVur0MxcNg9cuOBQh7v+NIJE4FYUkQTRZKKeUFefqAF4EOjzlWrF32CJw9Djd/BMElMpcEyXjSO2NJEMCnCcNnYxYi0kVEdohIvIiMz+F8GRH5UkR+E5GtIjLU3WuVUqogueQ1n2LugRteg10r4KPekJxUYJYE8UmyEBEr8AbQFbgK6C8iV2UrNgr4wxjTAIgFXhKRYDevVUqpAiNfH/CNB0Hf6ZC4AWZ0JyXpnxyL+XpJEF+1LJoD8caY3caYFGAW0DNbGQOUEhEBQoGjQKqb1yqlVIGR7zWf6vWCW2fD0d3MDXmSSDl8QRFfLwniq2QRAezL8jrBeSyr14G6wH5gM3CvMSbdzWsBEJERIvKLiPxy+PCF31yllPIFj6z5dGU7GLSAykGn+CL4cWpIQuYpfywJ4qtkkdPO5dknFHcGNgHhQEPgdREp7ea1joPGTDXGNDXGNK1YseKlR6uUUvngsTWfqjYjePhiyhS38kWxJ2kgu3JdEmT+xkRi4lZQbfwiYuJWeHxPDF8liwSgapbXkThaEFkNBeYah3hgD1DHzWuVUqrA6NUownNrPlWuR/E7llMmrDwLSj3PmputOSYKb2+i5JO1oUQkCPgTaA8kAuuBAcaYrVnKvAUcNMY8LiKVgV+BBkCSq2tzomtDKaUCyvED8HEf+DfeMQBet3vmqZi4FSTmMB4SERbCmvHt3H4Lv68NZYxJBUYDS4BtwGxjzFYRuVNE7nQWewpoJSKbgW+BccaYI7ld64u4lVKqwChdBYYsgioNYPZtsGlm5ilfbKLks4fyjDFfA19nO/Z2lq/3A53cvVYppYqcEuXgtvnw2UCYP9LxtHfLkYSHheTYsvDkjCldSFAppQqTYqEw4DOoewMsHg8rnmFsp1pe30RJk4VSShU2QcXgphnQ6Db4/gV6/TOZ53rX8+omSro2lFJKFUYWq2NpkJAwWPsavaKT6DX2TbDavPJ2miyUUqqwEoGOTzkWIPz2SccChDfNAJvnn+7WbiillCrMRKD1A45lzv9cAh/fCCmnPP422rJQSqlA0Ox2R5fUrhUQ5PmWhSYLpZQKFPVvdPzxAu2GUkop5ZImC6WUUi5pslBKKeWSJgullFIuabJQSinlkiYLpZRSLmmyUEop5ZImC6WUUi75ZKc8fxCRw8Bfl3h5BeCIB8MpDIpinaFo1rso1hmKZr3zWucrjDEVczoRsMkiP0Tkl9y2FgxURbHOUDTrXRTrDEWz3p6ss3ZDKaWUckmThVJKKZc0WeRsqr8D8IOiWGcomvUuinWGollvj9VZxyyUUkq5pC0LpZRSLmmyUEop5ZImiyxEpIuI7BCReBEZ7+94vEVEqorIShHZJiJbReRe5/FyIrJMRHY6/y7r71g9TUSsIrJRRL5yvi4KdQ4TkTkist35b35NoNdbRO53/mxvEZFPRaR4INZZRKaJyCER2ZLlWK71FJEJzs+3HSLSOS/vpcnCSUSswBtAV+AqoL+IXOXfqLwmFXjAGFMXaAmMctZ1PPCtMaYm8K3zdaC5F9iW5XVRqPNkYLExpg7QAEf9A7beIhIB3AM0NcbUB6zALQRmnWcAXbIdy7Gezv/jtwD1nNe86fzcc4smi3OaA/HGmN3GmBRgFtDTzzF5hTHmgDHmV+fXJ3B8eETgqO8HzmIfAL38EqCXiEgk0A14L8vhQK9zaaAN8D6AMSbFGJNEgNcbx5bRISISBJQA9hOAdTbGfA8czXY4t3r2BGYZY84aY/YA8Tg+99yiyeKcCGBfltcJzmMBTUSigEbAz0BlY8wBcCQUoJIfQ/OGV4H/A9KzHAv0OlcHDgPTnd1v74lISQK43saYROBF4G/gAHDMGLOUAK5zNrnVM1+fcZoszpEcjgX0vGIRCQW+AO4zxhz3dzzeJCLdgUPGmA3+jsXHgoDGwFvGmEbAKQKj+yVXzj76nkA1IBwoKSID/RtVgZCvzzhNFuckAFWzvI7E0XQNSCJiw5EoPjHGzHUePigiVZznqwCH/BWfF8QAN4jIXhxdjO1E5GMCu87g+LlOMMb87Hw9B0fyCOR6dwD2GGMOG2PswFygFYFd56xyq2e+PuM0WZyzHqgpItVEJBjHQNBCP8fkFSIiOPqwtxljXs5yaiEw2Pn1YGCBr2PzFmPMBGNMpDEmCse/7QpjzEACuM4Axph/gH0iUtt5qD3wB4Fd77+BliJSwvmz3h7HuFwg1zmr3Oq5ELhFRIqJSDWgJrDO3ZvqE9xZiMj1OPq1rcA0Y8wz/o3IO0TkWuAHYDPn+u8n4hi3mA1cjuM/3E3GmOyDZ4WeiMQCDxpjuotIeQK8ziLSEMegfjCwGxiK4xfFgK23iDwB9MMx828jMBwIJcDqLCKfArE4liI/CDwGzCeXeorIQ8AwHN+X+4wx37j9XposlFJKuaLdUEoppVzSZKGUUsolTRZKKaVc0mShlFLKJU0WSimlXNJkoZRSyiVNFkoppVzSZKGUl4nIOyIS4+84lMoPTRZKeV8L4Cd/B6FUfmiyUOoSiEhp55LfW0XktIhsEpGfRMSSrVxd4E9jTFqWY1HOXevec+7k9omIdBCRNc7dzZo7y5UUkUUi8puzXD8fV1OpTEH+DkCpwsi5pHsj5wf7Q8aY3DbK6goszuF4DeAmYASORSwHANcCN+BYp6sXjt3M9htjugGISBlP1kGpvNCWhVL5Ux/YepHznck5Wewxxmw2xqQ7r//WOBZq2wxEOctsBjqIyPMi0toYc8yDcSuVJ5oslMqfq4AtOZ0QkRJAmDEmpz0Dzmb5Oj3L63ScLX5jzJ9AExxJ4zkRedRTQSuVV5oslMqfcOCfXM61BVZe6o1FJBw4bYz5GMc2oY0v9V5K5ZeOWSiVP0uA90VkiDHmu2znuuLYme5SRQOTRCQdsAMj83EvpfJF97NQyktE5FeghXNrT6UKNU0WSimlXNIxC6WUUi5pslBKKeWSJgullFIuabJQSinlkiYLpZRSLmmyUEop5ZImC6WUUi79P0D9Q3HFysOOAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(precession_times_2_loaded, probabilities_echo_loaded,\n",
    "            label=\"Simulation Data\")\n",
    "plt.plot(precession_times_fine, T2_echo_decay(\n",
    "         precession_times_fine, *popt_echo),\n",
    "         color=\"tab:orange\",\n",
    "         label=\"Fit, $T_{2,echo}$ = \" + \"{:.2f} ms\".format(popt_echo[0]))\n",
    "plt.xlabel(\"$\\\\tau$ / ms\")\n",
    "plt.ylabel(\"$\\langle P(S) \\\\rangle$\")\n",
    "plt.legend()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "We see that the spin echo increases the coherence time in comparison to free\n",
    "decay. Also, the coherence time matches the order of magnitude by the\n",
    "theoretically predicted (c.f., above mentioned reference)\n",
    "\n",
    "\\begin{align}\n",
    "T_{2,echo}^{1+\\alpha} = \\frac{1}{S_0 (dJ/d\\varepsilon)^2 (2-2^{-\\alpha+1})\n",
    "\\Gamma(-1-\\alpha) \\sin(\\pi \\alpha/2)}\n",
    "\\end{align}\n",
    "\n",
    "for $\\alpha = 1$ leading to $T_{2,echo} \\approx 137~\\text{ms}$."
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We see that the spin echo increases the coherence time in comparison to free\n",
    "decay. Also, the coherence time matches the order of magnitude by the\n",
    "theoretically predicted (c.f., above mentioned reference)\n",
    "\n",
    "\\begin{align}\n",
    "T_{2,echo}^{1+\\alpha} = \\frac{1}{S_0 (dJ/d\\varepsilon)^2 (2-2^{-\\alpha+1})\n",
    "\\Gamma(-1-\\alpha) \\sin(\\pi \\alpha/2)}\n",
    "\\end{align}\n",
    "\n",
    "for $\\alpha = 1$ leading to $T_{2,echo} \\approx 137~\\text{ms}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "outputs": [],
   "source": [
    "probabilities_echo = np.zeros(len(precession_times_2))\n",
    "\n",
    "measurement_state_echo = initial_state\n",
    "\n",
    "for i, tau in enumerate(precession_times_2):\n",
    "    \n",
    "    solver = create_monte_carlo_solver(tau, echo_pulse)\n",
    "    probabilities_echo[i] = average_monte_carlo_traces(solver, initial_state,\n",
    "                                                       measurement_state_echo)\n",
    "    "
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# np.save(r\"..\\stored\\precession_times_2\", precession_times_2,)\n",
    "# np.save(r\"..\\stored\\probabilities_echo_2\", probabilities_echo,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "outputs": [],
   "source": [
    "precession_times_2_loaded = np.load(r\"..\\stored\\precession_times_2.npy\")\n",
    "probabilities_echo_loaded = np.load(r\"..\\stored\\probabilities_echo_2.npy\")\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit parameters:  [135.30109118   1.90271031]\n"
     ]
    }
   ],
   "source": [
    "def T2_echo_decay(t, T2_echo, alpha):\n",
    "    return 0.5 * (1 + np.exp(-(t / T2_echo) ** alpha))\n",
    "    \n",
    "popt_echo, cov_echo_al = curve_fit(T2_echo_decay, precession_times_2_loaded,\n",
    "                                   probabilities_echo_loaded, p0=[136, 2])\n",
    "print(\"Fit parameters: \", popt_echo)\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEJCAYAAABlmAtYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0m0lEQVR4nO3deVxU1fvA8c8zw6AoKu4JWGiuKbkvSRrumppLlmnmll/LtM3y59K+U7Zpu5XaZmbmVpZbaqVWmmmpqYlLCZpLhivKAOf3xwyICM4gs8DwvF8vXzL3nnvnOYjzcJZ7jhhjUEoppS7G4u8AlFJKFXyaLJRSSrmkyUIppZRLmiyUUkq5pMlCKaWUS0H+DsBbKlSoYKKiovwdhlJKFRobNmw4YoypmNO5gE0WUVFR/PLLL/4OQymlCg0R+Su3c9oNpZRSyiVNFkoppVzSZKGUUsqlgB2zUKowstvtJCQkcObMGX+HogJY8eLFiYyMxGazuX2NJgulCpCEhARKlSpFVFQUIuLvcFQAMsbw77//kpCQQLVq1dy+zifdUCIyTUQOiciWXM6LiEwRkXgR+V1EGmc510VEdjjPjfdmnPM3JhITt4Jq4xcRE7eC+RsTvfl2Sl3gzJkzlC9fXhOF8hoRoXz58nluvfqqZTEDeB34MJfzXYGazj8tgLeAFiJiBd4AOgIJwHoRWWiM+cPTAc7fmMiEuZtJtqcBkJiUzIS5mwHo1SjivHKTluxgf1Iy4WEhjO1c+7zzSuWXJgrlbZfyM+aTZGGM+V5Eoi5SpCfwoXGsl/6TiISJSBUgCog3xuwGEJFZzrIeTxaTluwg2Z7G3da5CJCKlZT0IPYsWgZp9SCoOD/vtzP3xyOUSS2OSAmOJ5VgwtyzgCYUpVRgKyhjFhHAviyvE5zHcjreIrebiMgIYATA5ZdfnqcA9iclA3Bn0JeUlLPnTqQCXzm+bAG0sALWc6fPmiD+W1AWNkRBaGV2ny3Nnl3pNEstx99Sib1Jl2lCUYXKM888w8yZM7FarVgsFt555x1atGjB8OHDGTNmDFdddVW+3yPjodkKFSrkWubZZ59l4sSJma9btWrF2rVr8/3eQ4YM4bvvvqN06dIkJyfTsmVLnnvuOSIiLv7/79VXX2XEiBGUKFEi3zEURgUlWeTUJjIXOZ4jY8xUYCpA06ZN87SrU3hYCIlJydQ7Ox0L6dhIxUYqVcvY+GZ0S7An023S15SWU5TmdObfFeQYlSSJG4Ot8G88FQ4ncL/lJASfu/cJE8KBhZfBn9FQoRbrk6vw/s+Gg/ZKGIJy7fJSytd+/PFHvvrqK3799VeKFSvGkSNHSElJAeC9997zaSzZk4UnEkWGSZMm0bdvX4wxvPrqq7Rt25YtW7YQHByc6zWvvvoqAwcOLLLJoqA8Z5EAVM3yOhLYf5HjHje2c21CbI4mQzoWzhJMmq0Ud3RpBqUug3LVSCpTlx/T67EkvRmfp8Xyftr1PJ/an5dLjoFBC2DUzzQ4M5XaZ2bQ/uwkhqU8yBP225iT1obE1NJwcCusfoVmvzzAl9YH2VpsGN8Ej+cV2xsMTF/A0m/mwtmT58Wlg+7qYjz983HgwAEqVKhAsWLFAKhQoQLh4eEAxMbGZi6hExoayrhx42jSpAkdOnRg3bp1xMbGUr16dRYuXAjAjBkzGD16dOa9u3fvzqpVqy54z169etGkSRPq1avH1KlTARg/fjzJyck0bNiQW2+9NfM9wTGbZ+zYsdSvX5/o6Gg+++wzAFatWkVsbCx9+/alTp063HrrrbjaCVREuP/++7nsssv45ptvABg5ciRNmzalXr16PPbYYwBMmTKF/fv307ZtW9q2bZtruUBWUFoWC4HRzjGJFsAxY8wBETkM1BSRakAicAswwBsBZPxGf7GuobGda583CA4QYrMytnPtzNeOFgrsMhHsMueujQgLYc097cB+husffZ9akkAdyz5qy9+0sGyjt3UNpABxj0LFuhDRmE3p1Znxawn22yMwWLQFos7j7qSMvOjUqRNPPvkktWrVokOHDvTr14/rrrvugnKnTp0iNjaW559/nt69e/Pwww+zbNky/vjjDwYPHswNN9zg9ntOmzaNcuXKkZycTLNmzbjxxhuJi4vj9ddfZ9OmTReUnzt3Lps2beK3337jyJEjNGvWjDZt2gCwceNGtm7dSnh4ODExMaxZs4Zrr73WZQyNGzdm+/bt9OzZk2eeeYZy5cqRlpZG+/bt+f3337nnnnt4+eWXWblyZWbXWU7lrr76arfrXdj4JFmIyKdALFBBRBKAxwAbgDHmbeBr4HogHjgNDHWeSxWR0cASHCMF04wxW70VZ69GERf9T+aRhGIrzrEydZmfFAXp5+5djuO0K5XAi61SIXEDbF9Ew+SjzLfCMUsJ1qXX4ef0uvycWpeXF1t1/ENlTsrIKtmexqQlOy753z80NJQNGzbwww8/sHLlSvr160dcXBxDhgw5r1xwcDBdunQBIDo6mmLFimGz2YiOjmbv3r15es8pU6Ywb948APbt28fOnTspX758ruVXr15N//79sVqtVK5cmeuuu47169dTunRpmjdvTmRkJAANGzZk7969biWLrC2Q2bNnM3XqVFJTUzlw4AB//PFHjknA3XKBwlezofq7OG+AUbmc+xpHMikQvJVQkm1lufb6NpBRzhjaTJxOY9lJC8s2Wli20dH2KwAnzoTAJ63hyvYss0czYelxku2OzKOtj6IjY1KGu8fdZbVaiY2NJTY2lujoaD744IMLkoXNZsucfmmxWDK7rSwWC6mpqQAEBQWRnn7uN6Kc5vWvWrWK5cuX8+OPP1KiRAliY2Ndzv+/WNdSRhwZ9ciIxZWNGzfSvn179uzZw4svvsj69espW7YsQ4YMyTEed8sFkoLSDRVQPJFQECGtTBTzkyozP93xm1El/qOFZRvtQnbS+9942LmUjsASqch3QQ34Lr0Ba9PrcdpePF+/XarCIWNSRk7HL9WOHTuwWCzUrFkTgE2bNnHFFVdc0r2ioqJ48803SU9PJzExkXXr1l1Q5tixY5QtW5YSJUqwfft2fvrpp8xzNpsNu91+wZIUbdq04Z133mHw4MEcPXqU77//nkmTJrF9+/Y8x2iM4bXXXuPAgQN06dKFbdu2UbJkScqUKcPBgwf55ptviI2NBaBUqVKcOHGCChUqcPz48VzLBSpNFn7iKqHAhS2QQ5RlubU17bvf5WiBHN3Dwy9N5jrL7/Sx/sBtQctJMVZ+Tq/L0hPN4HhdKF0F0K6qQOTOGFpenTx5krvvvpukpCSCgoKoUaNG5qBzXsXExFCtWjWio6OpX78+jRs3vqBMly5dePvtt7n66qupXbs2LVu2zDw3YsQIrr76aho3bswnn3ySebx37978+OOPNGjQABHhhRde4LLLLstTshg7dixPPfUUp0+fpmXLlqxcuZLg4GAaNGhAo0aNqFevHtWrVycmJua8eLp27UqVKlVYuXJlruUClbiaLVBYNW3a1ATC5keuPuRj4laQmJRMMHaaWP7kOstvdLRs4ErLAUeByGZsKd2GBzdXZbu9UuZ1ITYrz/WJ1oRRwGzbto26deu6XV5/CVCXKqefNRHZYIxpmlN5TRaFXPYZMQAhNguvdyhBe1kH276EA78BsD29KovSWrAgPYa/TWXHDK3x7c67l37w+Fdek4VSlyqvyUK7oQq53MY/2jeKAK6DNmO5dvwMOlp/oYt1PQ/Y5vAAc9iYXoOFJ1rByfoQWskr0zCVUoFDk0UAcDX+YcIuZ3pSRaandaUK/9LDupae1rU8ZvsQXvoYqsey+e/6YG8IFM+8Lr/TMJVSgaOgPMGtvCjr0+kHKM/UtB70NS/wbduFcO0Y+HcXj9insK7YKJ4NepeGEk/Gqir5nYaplAoM2rIoAlx2VbV7mJHPvUG700voZV3LgKCV7EiPZHZaLD+Gdrzgfjq2oVTRo8miiLhoV5UInbv2ZsLcGjxxdhA9rD/Sz7qKR2wfk54yCz6fD83+B1e0Yv6m/Tq2oVQRpMlCAee3PmYltef7Ut156hqh3ekl8NtM2DoPKtVjx9E2YG+Gjm0oVbTomIXK1KtRBGvGt2NPXDfWjG9Hu+vaQtc4GLMdekwBsTAu9S1+Ljaah4M+4gr5J/Pa7GMbulpu4Wa1WmnYsGHmn71799KqVSsAkpKSePPNN13e491338283mKxZH49ZswYb4efb8OGDaNSpUrUr1//vONnzpyhefPmNGjQ4ILVZqOiooiOjqZhw4Y0bZrj7NOLXr948WJq165NjRo1iIuL807F8kGfs1DuM4Y7n3uT65O/pKtlHTZJ49u0Rryb1o19pRqzZkJ7ILdnP/QhQHcUlOcsQkNDOXnyZI7n9u7dS/fu3dmyZYtb90pMTKRVq1b89ddfngzRq77//ntCQ0MZNGjQefU0xnDq1ClCQ0Ox2+1ce+21TJ48mZYtW7q1oVNu1zdr1oxatWqxbNkyIiMjadasGZ9++qlHNprKTV6fs9CWhXKfCF269mIc99Hq7BReTe1DA8suZgU/zVchj8KWLyAt9aKroarCK2M/ifHjx7Nr1y4aNmzI2LFjXV63ZcsWoqOjvR2eR7Vp04Zy5cpdcFxEMr8Pdrsdu92ep/2sc7t+3bp11KhRg+rVqxMcHMwtt9zCggULLrh+79691KlTh+HDh1O/fn1uvfVWli9fTkxMDDVr1mTdunWcOnWKbt260aBBA+rXr5+530d+6ZiFypOsYxuTk/qyoOTNvFJnGw33fQxzhkHY5XQ6EctnxHI6y7gG5NxVpbOqCqaMjYcAqlWrlrmEOEBcXBxbtmzJca+JnGzevPmC7hx/aN26NSdOnLjg+IsvvkiHDh3cvk9aWhpNmjQhPj6eUaNG0aKFY6dnEaFTp06ICHfccQcjRoxw+/o5c+ZQteq5fd4iIyP5+eefc7w+Pj6ezz//nKlTp9KsWTNmzpzJ6tWrWbhwIc8++yy33XYb4eHhLFq0CHAs1ugJmixUnl04s6orpN8HO76Gta/xWNKH3Bc0hw/TOjEttQv/URo4fzVUfWLcDd+Mh382e/ael0U7xqFcCAkJcTsZuLJlyxY6djw3BXv+/PksWrSIQ4cOMWrUKDp16nRJ9x0+fHietnr94YcfLul9srNarWzatImkpCR69+7Nli1bqF+/PmvWrCE8PJxDhw7RsWNH6tSpk7kpk6vrcxoOyK3FkrE4I0C9evVo3749IpK5l0h0dDQPPvgg48aNo3v37rRu3doj9dZuKOUZFgvU7Q63L+G71jNZRz1GWRewuti9TAj6hKq2E+ethqpdVUVH9pZFr169ePfdd5kxY0aOXSR2u52JEydy3333cddddwGQmprKmDFjeOCBB5gyZQqnT58mPj6ehx56iN69e5OcnMyYMWMYNWoUDz/8cI5xtG7d+rxB+4w/y5cvv6R6hYWFERsby+LFiwEyt5+tVKkSvXv3znFJ9tyuj4yMZN++fZnnEhISMu+XXdY9O3LaS6RWrVps2LCB6OhoJkyYwJNPPnlJ9ctOWxbK465r34355Rpy2zfLuSl5FsODvmG4ZTnWf4bC8XuhdLjXNu4JKG60APwhY1+HrNq3b8+HH35IRMT5rcL09HR27txJnTp1LrjP008/zahRF+55NnXqVJKTkwkLC2P37t0AvPXWW/Ts2TNzi9fVq1fTtWtXxo0bx8CBA3nttdcYMGAATZs2pW/fvjnG7YmWxeHDh7HZbISFhZGcnMzy5csZN24cp06dIj09nVKlSnHq1CmWLl3Ko48+6vb1zZo1Y+fOnezZs4eIiAhmzZrFzJkzLynG/fv3U65cOQYOHEhoaCgzZszIZ60dNFkor3B0VQ0GBsO/u+CHl2Ddu/DLNGh0G43KNOfXY6UuuC4/G/co3yhfvjwxMTHUr1+frl278vzzzxMfH5/jgHB8fDyRkZHn/TZsjGH8+PF07do1xz0uNm7cyBtvvHHeNb/++isjR47MfL1+/Xp69OgBOLp1tm7dyr333ktKSgolSpTIdx379+/PqlWrOHLkCJGRkTzxxBPcfvvtHDhwgMGDB5OWlkZ6ejo333wz3bt3Z/fu3fTu3RtwtIIGDBiQue3s9ddfz3vvvUd4eHiu1wO8/vrrdO7cmbS0NIYNG0a9evUuKfbNmzczduxYLBYLNpuNt956K9/fD9Cps8qX/tsLq1+BjZ+QbgyfpbXllZReHKIsoNNroeBMnc2LLVu2MG3aNF5++WW3yk+ZMoUPPviAZs2a0bBhQ+68804OHjzIV199xe23386XX37JzJkzqVq1Ku3ataNLly4sWLCAhQsXUq5cOSZMmMC4ceN4++23+e+//4iLi6Nt27bMnj2bEiVKcM899xS676E/6H4WTposCrBjCfDDS6Rv+BC7sfBBakfml7yZEV2aFelEAYUzWXjC4sWLCQ4Opl27dq4LK4/Q/SxUwVcmErq/gqXVPRRbFceI3z9jhPkekkbBmVFQvLS/I1Q+ltFlowounQ2l/KdcNejzDtz1I1wZC9/FweQGsGYy2B0D3bpsiFIFg7YslP9Vqgv9PobEX2HF07DsUfj5HTbUuJuJ66/gtN3RVarPYijlP9qyUAVHRGO4bS4M/gpKVqTJr+OZLRO4xrI1s4g+i6GUf2iyUAVPtdbwv5Xcl3IXYXKST4Of4T3bJK4URxeUrnCrlO9pslAFk8XC+tIdaX/2ReLst9Dcsp0lweN4Ouh96pU5m1ksY9mQxKRkDOe6qgpzwgjUGYqq4LiUnzFNFqrAGtu5NhZbCG+n3UDs2Vf4OK0D/ayrmJc6Cn54GVLPBtyyIcWLF+fff//VhKG8xhjDv//+S/HixV0XzkIHuFWBdf7e4fBu6EgiYu6jY8Lr8O0TsPFjah3vQyKNLri2sC4bEhkZSUJCAocPH/Z3KCqAFS9enMjIyDxdo8lCFWg57x0eAzuXw+LxTA+exLdpjXgy9Tb+Mpdlliisy4bYbDaqVavm7zCUuoB2Q6nCqWYHGLmWLfUepKVlG0uD/4//C5pFCc4QYrOet8KtUir/NFmowisomPo3PcL3XZawwnotdwUt5LviD/JR87/o1TDn5Z2VUpdG14ZSgWPfOvh6LBzYBFGtodvLULGWv6NSqtDQPbhV0VC1OfxvBXR/Bf75Hd5qBSueyVw6RCl16TRZqMBisULTYTD6F6jfB75/Ad68Bnat8HdkShVqmixUYAqtBH2mwqAFIBb4qDfMuR1OHPR3ZEoVSjp1VgW26rEwci2sfoW0H17i9Javed7ej1Wh3XmwS11dkFApN2nLQgU+W3Hmhw2iu/0FfkurxtO26bySPIGpcxcX6mVBlPIlTRaqSJi0ZAfb7JUZaJ/ImJQ7qSmJzLP8HwcXPQNpdn+Hp1SB57NkISJdRGSHiMSLyPgczpcVkXki8ruIrBOR+lnO7RWRzSKySUR0PqzKs3PLfwhz09vQ4eyLLE9vzB2pn8C7beHAb4CuYKtUbnwyZiEiVuANoCOQAKwXkYXGmD+yFJsIbDLG9BaROs7y7bOcb2uMOeKLeFXgCQ8LITHLelFHKMMo+32sKb6JZ09+AFPbsqPGMB7b3oZjdiugmy0plZWvWhbNgXhjzG5jTAowC+iZrcxVwLcAxpjtQJSIVPZRfCrAje1cmxCb9bxjITYrza8fAqN+hgb9qb3zXebK/9FEzq1YW5hXsFXKk3yVLCKAfVleJziPZfUb0AdARJoDVwAZyyIaYKmIbBCREbm9iYiMEJFfROQXXbVTZdWrUQTP9YkmIiwEASLCQniuT7SjxRBSFnq9waCU8RQTO58HP8ljQR9QHMe+GYV1BVulPMlXU2clh2PZ1xmJAyaLyCZgM7ARSHWeizHG7BeRSsAyEdlujPn+ghsaMxWYCo7lPjwVvAoMOa9ge86u0i3olFSLsUGfMTRoCddZfuNB+50cLNPAh1EqVTD5qmWRAFTN8joS2J+1gDHmuDFmqDGmITAIqAjscZ7b7/z7EDAPR7eWUh41tnNtjK0kT6QOpn/KQwRLKp8HP8H0qosg9azrGygVwHyVLNYDNUWkmogEA7cAC7MWEJEw5zmA4cD3xpjjIlJSREo5y5QEOgFbfBS3KkKydlX9lF6PocUns++KPtTa+R5MbQsHfvd3iEr5jU+6oYwxqSIyGlgCWIFpxpitInKn8/zbQF3gQxFJA/4AbndeXhmYJyIZ8c40xiz2Rdyq6Lmwq6o7/NkPFt7tmGJ73XgWhN7MC8t2sT8pmfCwEMZ2rq2zpVTA0yXKlXLH6aPw9YOw5Qt+N1dyf8qd7DKOBBFis54bLFeqENMlypXKrxLloO80HrE9SCQHWRQ8kcHWJYDR6bWqSNBkoVQefHyiMZ3PvsDa9Ho8YfuA6bYXqMAxnV6rAp4mC6XyIDwshMOEMcw+loftQ7nG8geLi43jxlJb/R2aUl6lyUKpPDj3JLjwcVpHuqc8w2HK8qL9GVj0oO7KpwKWJgul8iD7k+DJZWqys8cCuGY0rH8X3rlOp9iqgKSzoZTylF0rYN5ISD4K7R+DlneB5eK/j83fmMikJTt0Gq4qEHQ2lFK+cGU7x658NTvB0ofg495w/ECuxedvTGTC3M0kJiVjOLfKrS6LrgoiTRZKeVLJ8tDvY+gxGfatg7djYOfyHItOWrKDZHvaecd0Gq4qqDRZKOVpItBkCIz4DkIvg09uhGWPXrAjX27TbXUariqINFko5S0Va8H/voUmQ2HNZJjeFZL+zjwdHhaS42W5HVfKnzRZKOVNthDo8Sr0nQ6Hd8Db18K2r4DcN2Qa27m2HwJV6uI0WSjlC/X7wB3fQdlq8Nmt8PX/0Su6Qu4bMilVwOjUWaV8KfUsLH8cfnoTqjRwtDjKX+nvqJQCdOqsUj43f2MiMXErqDZ+ETFxK85Nhw0qBl2eg1s+hf/+gnfawOY5/g1WKTdoslDKw9x6fqLO9XDnaqhcH764Hb68D+xn/BWyUi5pslDKw9x+fiKsKgxZBDH3wobpMK2zo7WhVAGkyUIpD8vT8xPWIOj4JNwyE47ucXRL/bnEyxEqlXeaLJTysEt6fqJON7hjlaO1MfNm+PYpSE/LvbxSPqbJQikPu+TnJ8pVh9uXQaPb4IcX4aPecPKwFyNVyn2aLJTysOzLmOfp+QlbCPR8HXq+Aft+dnRL/f2z12NWyhV9zkKpgurA7zB7EOlJ+3jNOohXT7YnPKyELmOuvEafs1CqMKpyNYuu+ZQV6Y24N3Uar9mm8F/Sf7qMufILTRZKFWDPrjjA8LP38ay9P10t65gX/CiVUxN0GXPlc5oslCrAHNNthalpPRhon0hFSWJh8CPUPr7W36GpIkaThVIFWNbptj+m16PH2Wf421TiveAX4bsXID0duMjyIkp5iCYLpQqw7NNwE6nIQPMUCZE9YOUz8NlAvlq3Q7dnVV6nyUKpAiynabiP92nC5bd/CF2ehz8XU/+b3oSn/n3edbo9q/I0nTqrVGG2dzVHpvenGHbG2EeyLP3crEcB9sR1819sqtDRqbNKBaqoa/lf8ZfYZarwbvDL3B/0OYJjHEO3Z1WepMlCqUJucJcYBpsnmJ16HfcGzeM920tUsp3R7VmVRwX5OwClVP5kPM09aXFJfj9ZncdsH/Ft6FOUqqqbKinP0TELpQLNX2vhs4GQlgo3TYca7f0dkSokdMxCqaLkilbwv5WO5c4/6Qs/vgkB+kuh8h3thlIqEJW9AoYtgXl3wJIJcHArCyMf5Pnle9iflEx4WIguSKjyxK2WhYiUFZFBIjJPRLaKyFci8j8RqeztAJVSl6hYKNz8EVw3DjZ9TOSXN3M26R99cE9dEpfJQkTmAnOBSsA4Y0w94C6gOPCRiKzyaoRKqUtnsUDbiTxse5C67GVBsYepJ3sBfXBP5Y07LYthxpi2xpgXjTF/Ahhj/jbGvGaM6QT08mqESql8++REY/qmPA7AnODHud7yE5D7fuFKZecyWRhjkrK+FpGSImLN7XxuRKSLiOwQkXgRGZ/D+bLObq7fRWSdiNR391ql1MWFh4Ww1UTR8+zTbDVRvBk8hfuD5hBRppi/Q1OFhDvdUBYRGSAii0TkELAdOOAcu5gkIjXduIcVeAPoClwF9BeRq7IVmwhsMsZcDQwCJufhWqXURWQsSHiEMgxIecj5AN9cZpd7C86e9Hd4qhBwpxtqJXAlMAG4zBhT1RhTCWgN/ATEichAF/doDsQbY3YbY1KAWUDPbGWuAr4FMMZsB6KcA+juXKuUuoisCxLasTG55L1srj+e8H9WwPQucEwHutXFuTN1toMxxp79oDHmKPAF8IWI2FzcIwLYl+V1AtAiW5nfgD7AahFpDlwBRLp5LQAiMgIYAXD55Ze7CEmpoqVXo4hsU2XbQ4Om8PlQeLcdDJjF/IOVmLRkh06vVRdwZ8zigkRxCWUkp8uyvY4DyorIJuBuYCOQ6ua1GXFMNcY0NcY0rVixoouQlFLU7Ai3LwVrMKnvd2HF3Pd0XwyVI3fGLL4VkXpZXt8gIg+LSI6/3eciAaia5XUksD9rAWPMcWPMUGNMQxxjFhWBPe5cq5TKh8pXwf++ZXv65UyxvsxI60Iyfh/T6bUqgztjFpHGmK0AItIK+Ai4HJguIr3dfJ/1QE0RqSYiwcAtwMKsBUQkzHkOYDjwvTHmuDvXKqXyKbQSfZMnsjDtGsbZZvFC0FRspAI6vVY5uDNmcTzL14OAt40x40SkEo4P7XmubmCMSRWR0cASwApMM8ZsFZE7neffBuoCH4pIGvAHcPvFrnW7hkopt5QPK8M9SaPZY6pwb9BcLrcc4s6U+ygZVsnfoakCwOWqsyLyKY6B7O+B34E+xpi1znMbjTGNvB7lJdBVZ5XKm/kbE5kwdzPJ9jR6Wlbzgm0qB6jIzg7v07F1jL/DUz6Q31Vn7wfuABKBX7MkChsQ6rEolVJ+lXV67cL0a7mn2FNUKXaWjmsGwJ7v/R2e8jO397MQEYsxJj3L605AX2PMCG8Flx/aslDKA/7bCzP7wb/x0P1VaHybvyNSXpSvloWICEDWROF8vTQjUWSUUUoFmLJRjqm11drAwtGw7FFIT3d5mQo8bj3BLSJ3i8h5T7mJSLCItBORD4DB3glPKeV3xcvAgM+h2XBYMxk+HwQpp/0dlfIxd2ZDdQGGAZ+KSHXgPxzLk1uBpcArxphNXotQKeV/1iC4/kUoXwMWT4DjPaD/LAjVh1+Lijztwe0c1K4AJLu72qy/6JiFUl6y7Sv4YjiUqgy3zoEKLtcSVYVEfscsBovIERE5CrwHnCzoiUIp5UV1u8OQrxyr1b7fEf5a6++IlA+4M2bxCNARqAP8DTzr1YiUUgVfZFMYvhxKVIAPe8KWL/wdkfIyd5LFcWPMRmPMIWPMIziWDFdKFXXlqjlmSkU0hTnDYPUrkIdubVW4uJMsqojICBFpLSIVAVfLkSuliooS5eC2eVD/Rlj+OHx1P6Sl+jsq5QXuzIZ6HLgauBWIBkJF5Gsc+0/8boz51HvhKaUKovkbE8/f96LTE/QKu9zRuji2D26aAcVK+TtM5UHutCzeNcaMNsZcZ4wpB1QHXgeSgOv1gTylipaMNaTO2/di3lbml/+f4ynvXSthelc4fsDfoSoPcidZLBORz0Skv4iUNsYkAKuAeOf5X70WnVKqwJm0ZAfJ9rTzjmXue9F0KAz4DI7ugffaw0FdIDpQuLNTXnvgCSAKWCQiP+HYKzsaxwN5BXLVWaWUd+S2v0Xm8ZodYeg3YNLh/c6wa4UPo1Pe4k7LAmPMH8aY54wxrYFYY8w1xpjHjTHaqlCqiAkPC3F9vMrVjqm1YVXhk5tg4yc+ik55i1vJIitjzBlvBKKUKhzGdq5NiM163rEQm5WxnWufd2z+bqHTsYmstteGBXfxx6xHLmlq7fyNicTEraDa+EXExK3QPcH9JM/JQilVtGXd90KAiLAQnusTTa9GEZllMgbB/zxmYah9HPPSYrhq+xR2f3AnpKflfvNschxMn7tZE4YfuDN1VimlztOrUcR5ySG7rIPgdoIYYx/JQVOOO/fOgtkn4cb3wJZzd1Zu98mQMZh+sfdXnqctC6WUx2UfBDdYiEvtzxP2QbB9EXxwA5w+muf7uDquvEeThVLK43IbBF9aqrfjgb0Dv8H7neC/vy7pPrkdV96jyUIp5XEXHQSv18uxRMipQ5x5ux1Dnn0/18FrdwfTlfdpslBKeZzLQfCoGL695kP+O2N47exDtLJsznHw2p3BdOUbedr8qDDRzY+UKthi4laQlpTAjOAXuFL286D9DhakX0tEWAhrxrfzd3hFUr42P1JKKW/Yn5TMP5Tn5pRH+SW9NpOD3+QO65fsT9L9vQsiTRZKKb/IGKQ+TkkG28fxZVpLJtg+5YWSM/P0LIbyDU0WSim/yDp4nYKNe+yjmZHejZvSFsGcoWDXxSIKEn0oTynlFxmD1Of2xShJWOdJkNwKlj4EJw9D/5kQUtbPkSrQAW6lVEG0eQ7MHwnlqsPAL6BMpL8jKhJ0gFspVbhE93UkiWOJjmXOD+/wd0RFniYLpVTBVK0NDF0EaSkwrTPsW+/viIo0TRZKqYKrSgO4fQkUD4MPesDOZf6OqMjSZKGUKtjKVYfbl0KFmvDpLfDbLH9HVCRpslBKFXyhlWDIIriiFcy7A9a+5u+IihxNFkqpwqF4abh1DlzVE5Y+DEsvbec9dWn0OQulVOERVAz6Toevx8LaKXDqCNwwBaw25m9MzPLMRghjO9fWBQc9SJOFUqpwsVih20sQWhlWPQun/+XLWs8yYWF85q56GSvYApowPES7oZRShY8IxI6Dbi/DzqVEfT2AYPux84pkbL+qPEOThVKq8Gp2O9z8AbXSd/F58BNcxr/nndbtVz3HZ8lCRLqIyA4RiReR8TmcLyMiX4rIbyKyVUSGZjm3V0Q2i8gmEdE1PJRS51zVkweDH6WKHOWLYo9zpZzbPEm3X/UcnyQLEbECbwBdgauA/iJyVbZio4A/jDENgFjgJREJznK+rTGmYW7rliiliq721/dlcPpjBJPKnOAnaCQ7dftVD/NVy6I5EG+M2W2MSQFmAT2zlTFAKRERIBQ4CqT6KD6lVCHWq1EEg/rcwKhiz3HMlGRmsWd5P+Y/Hdz2IF8liwhgX5bXCc5jWb0O1AX2A5uBe40x6c5zBlgqIhtEZERubyIiI0TkFxH55fDhw56LXilV4PVqFMHsibcS9X+rCbmsNq3WjYbfPvN3WAHDV8lCcjiW/WmazsAmIBxoCLwuIqWd52KMMY1xdGONEpE2Ob2JMWaqMaapMaZpxYoVPRK4UqqQOe9p7xGw9nV/RxQQfJUsEoCqWV5H4mhBZDUUmGsc4oE9QB0AY8x+59+HgHk4urWUUipn5z3t/RAsf1yf9s4nXyWL9UBNEanmHLS+BViYrczfQHsAEakM1AZ2i0hJESnlPF4S6ARs8VHcSqnCKuNp7yZDYfUr8OW9urd3PvjkCW5jTKqIjAaWAFZgmjFmq4jc6Tz/NvAUMENENuPothpnjDkiItWBeY5xb4KAmcaYxb6IWylVyFms0P0VKFEOfngJziRBn3cdiUTlic+W+zDGfA18ne3Y21m+3o+j1ZD9ut1AA68HqJQKTCLQ/lEIKefokjpzDPp9AsVC/R1ZoaJPcCulioZWo6HXW7DnB/jwBjh91N8RFSqaLJRSRUfDAdDvY/hnC0zr4tjjW7lFk4VSqmipcz3cNheO73fs7X0k3t8RFQqaLJRSRU/UtTDkK7AnOxLG/k3+jqjA02ShlCqawhvCsCVgC4EZ3WHvan9HVKBpslBKFV0VajgSRulw+KgPbP/a9TVFlO6Up5Qq2spEwNBv4JO+8NlANjR8mnu21dHtWbPRloVSSpUsD4MXcqh8M5psnECXE19gOLc96/yNOmtKk4VSSgEUK8XNJ8bwdVpzHrF9zANBswGj27M6abJQSimnv46lMtp+D5+mtuXuoPk8HTQNC+m6PSs6ZqGUUpnCw0IcXU+pw/mPUtwVtJAycooXS4zxd2h+py0LpZRyGtu5NiE2KyC8kHoLz9gH0MP6E7NLT4aUU/4Oz680WSillFOvRhE81yeaiLAQBPi61E382vBpKh/5ET7sWaTXkxIToBuCNG3a1Pzyyy/+DkMpVYDM35jIpCU78j4tdtuXMGcYlK8BA+dC6SreD9YPRGSDMaZpTud0zEIpVSTM35jIhLmbSbY7NkDKmBYLuE4YdXvAwC/g0/4wrRPLmrzD42vOFKlnMbQbSilVJExasiMzUWTI07TYam1g8JecPX2ChstvofSx7UXqWQxNFkqpIiG36a95mhYb0ZghPEkKQXwW/BRNZTuQx6RTSGmyUEoVCeFhIXk6npufjpen79nHOWzK8FFwHLGWjUAek04hpMlCKVUknJsWe06IzcrYzrXzdJ/wsBAOUJ6bUh4j3oTzru1leljW5jnpFDaaLJRSRUL2abERYSE81yc6zwPTGUnnKKXpn/IwG0wtJtve4I1aG70TeAGhU2eVUiqPsk7BjSpjYWbZd6jyz0po9wi0fgBE/B3iJdGps0op5UG9GkWc3yJJ6wALRsGKp+BMEnR8qtAmjNxoslBKqfyy2qDX21C8DKx9DZL/g+6TwRo4H7GBUxOllPIniwW6vgAhZeG75+HMMbjxfQgq5u/IPEIHuJVSylNEoO1E6BLnWCJk5s1w9qS/o/IIbVkopZSntRzp6JJaMNqxAOGtn0OJcjkWveT1qnxMWxZKKeUNDQdAv4/gn80woxscP3BBkYz1qhKTkgv80iGaLJRSylvqdIOBcyDpb5jWGY7uPu90vter8iFNFkop5U3V2sDghXD2BEzrAge3Zp7yyHpVPqLJQimlvC2iCQz9BsQK07vCvnWA59ar8gVNFkop5QuV6sCwxVCivGPQO/5bj61X5QuaLJRSylfKXgHDlkC5K2FmP3oFr/fIelW+oGtDKaWUryUnwcx+kLAOekyGxoP8HRFw8bWhtGWhlFK+FhIGt82DK9vDwrthzWR/R+SSJgullPKH4BJwy0yo1weWPQrLH4cC3NOjT3ArpZS/BAXDje85WhqrX3F0T3V7CSxWV1f6nCYLpZTyJ4sVur0MxcNg9cuOBQh7v+NIJE4FYUkQTRZKKeUFefqAF4EOjzlWrF32CJw9Djd/BMElMpcEyXjSO2NJEMCnCcNnYxYi0kVEdohIvIiMz+F8GRH5UkR+E5GtIjLU3WuVUqogueQ1n2LugRteg10r4KPekJxUYJYE8UmyEBEr8AbQFbgK6C8iV2UrNgr4wxjTAIgFXhKRYDevVUqpAiNfH/CNB0Hf6ZC4AWZ0JyXpnxyL+XpJEF+1LJoD8caY3caYFGAW0DNbGQOUEhEBQoGjQKqb1yqlVIGR7zWf6vWCW2fD0d3MDXmSSDl8QRFfLwniq2QRAezL8jrBeSyr14G6wH5gM3CvMSbdzWsBEJERIvKLiPxy+PCF31yllPIFj6z5dGU7GLSAykGn+CL4cWpIQuYpfywJ4qtkkdPO5dknFHcGNgHhQEPgdREp7ea1joPGTDXGNDXGNK1YseKlR6uUUvngsTWfqjYjePhiyhS38kWxJ2kgu3JdEmT+xkRi4lZQbfwiYuJWeHxPDF8liwSgapbXkThaEFkNBeYah3hgD1DHzWuVUqrA6NUownNrPlWuR/E7llMmrDwLSj3PmputOSYKb2+i5JO1oUQkCPgTaA8kAuuBAcaYrVnKvAUcNMY8LiKVgV+BBkCSq2tzomtDKaUCyvED8HEf+DfeMQBet3vmqZi4FSTmMB4SERbCmvHt3H4Lv68NZYxJBUYDS4BtwGxjzFYRuVNE7nQWewpoJSKbgW+BccaYI7ld64u4lVKqwChdBYYsgioNYPZtsGlm5ilfbKLks4fyjDFfA19nO/Z2lq/3A53cvVYppYqcEuXgtvnw2UCYP9LxtHfLkYSHheTYsvDkjCldSFAppQqTYqEw4DOoewMsHg8rnmFsp1pe30RJk4VSShU2QcXgphnQ6Db4/gV6/TOZ53rX8+omSro2lFJKFUYWq2NpkJAwWPsavaKT6DX2TbDavPJ2miyUUqqwEoGOTzkWIPz2SccChDfNAJvnn+7WbiillCrMRKD1A45lzv9cAh/fCCmnPP422rJQSqlA0Ox2R5fUrhUQ5PmWhSYLpZQKFPVvdPzxAu2GUkop5ZImC6WUUi5pslBKKeWSJgullFIuabJQSinlkiYLpZRSLmmyUEop5ZImC6WUUi75ZKc8fxCRw8Bfl3h5BeCIB8MpDIpinaFo1rso1hmKZr3zWucrjDEVczoRsMkiP0Tkl9y2FgxURbHOUDTrXRTrDEWz3p6ss3ZDKaWUckmThVJKKZc0WeRsqr8D8IOiWGcomvUuinWGollvj9VZxyyUUkq5pC0LpZRSLmmyUEop5ZImiyxEpIuI7BCReBEZ7+94vEVEqorIShHZJiJbReRe5/FyIrJMRHY6/y7r71g9TUSsIrJRRL5yvi4KdQ4TkTkist35b35NoNdbRO53/mxvEZFPRaR4INZZRKaJyCER2ZLlWK71FJEJzs+3HSLSOS/vpcnCSUSswBtAV+AqoL+IXOXfqLwmFXjAGFMXaAmMctZ1PPCtMaYm8K3zdaC5F9iW5XVRqPNkYLExpg7QAEf9A7beIhIB3AM0NcbUB6zALQRmnWcAXbIdy7Gezv/jtwD1nNe86fzcc4smi3OaA/HGmN3GmBRgFtDTzzF5hTHmgDHmV+fXJ3B8eETgqO8HzmIfAL38EqCXiEgk0A14L8vhQK9zaaAN8D6AMSbFGJNEgNcbx5bRISISBJQA9hOAdTbGfA8czXY4t3r2BGYZY84aY/YA8Tg+99yiyeKcCGBfltcJzmMBTUSigEbAz0BlY8wBcCQUoJIfQ/OGV4H/A9KzHAv0OlcHDgPTnd1v74lISQK43saYROBF4G/gAHDMGLOUAK5zNrnVM1+fcZoszpEcjgX0vGIRCQW+AO4zxhz3dzzeJCLdgUPGmA3+jsXHgoDGwFvGmEbAKQKj+yVXzj76nkA1IBwoKSID/RtVgZCvzzhNFuckAFWzvI7E0XQNSCJiw5EoPjHGzHUePigiVZznqwCH/BWfF8QAN4jIXhxdjO1E5GMCu87g+LlOMMb87Hw9B0fyCOR6dwD2GGMOG2PswFygFYFd56xyq2e+PuM0WZyzHqgpItVEJBjHQNBCP8fkFSIiOPqwtxljXs5yaiEw2Pn1YGCBr2PzFmPMBGNMpDEmCse/7QpjzEACuM4Axph/gH0iUtt5qD3wB4Fd77+BliJSwvmz3h7HuFwg1zmr3Oq5ELhFRIqJSDWgJrDO3ZvqE9xZiMj1OPq1rcA0Y8wz/o3IO0TkWuAHYDPn+u8n4hi3mA1cjuM/3E3GmOyDZ4WeiMQCDxpjuotIeQK8ziLSEMegfjCwGxiK4xfFgK23iDwB9MMx828jMBwIJcDqLCKfArE4liI/CDwGzCeXeorIQ8AwHN+X+4wx37j9XposlFJKuaLdUEoppVzSZKGUUsolTRZKKaVc0mShlFLKJU0WSimlXNJkoZRSyiVNFkoppVzSZKGUl4nIOyIS4+84lMoPTRZKeV8L4Cd/B6FUfmiyUOoSiEhp55LfW0XktIhsEpGfRMSSrVxd4E9jTFqWY1HOXevec+7k9omIdBCRNc7dzZo7y5UUkUUi8puzXD8fV1OpTEH+DkCpwsi5pHsj5wf7Q8aY3DbK6goszuF4DeAmYASORSwHANcCN+BYp6sXjt3M9htjugGISBlP1kGpvNCWhVL5Ux/YepHznck5Wewxxmw2xqQ7r//WOBZq2wxEOctsBjqIyPMi0toYc8yDcSuVJ5oslMqfq4AtOZ0QkRJAmDEmpz0Dzmb5Oj3L63ScLX5jzJ9AExxJ4zkRedRTQSuVV5oslMqfcOCfXM61BVZe6o1FJBw4bYz5GMc2oY0v9V5K5ZeOWSiVP0uA90VkiDHmu2znuuLYme5SRQOTRCQdsAMj83EvpfJF97NQyktE5FeghXNrT6UKNU0WSimlXNIxC6WUUi5pslBKKeWSJgullFIuabJQSinlkiYLpZRSLmmyUEop5ZImC6WUUi79P0D9Q3HFysOOAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(precession_times_2_loaded, probabilities_echo_loaded,\n",
    "            label=\"Simulation Data\")\n",
    "plt.plot(precession_times_fine, T2_echo_decay(\n",
    "         precession_times_fine, *popt_echo),\n",
    "         color=\"tab:orange\",\n",
    "         label=\"Fit, $T_{2,echo}$ = \" + \"{:.2f} ms\".format(popt_echo[0]))\n",
    "plt.xlabel(\"$\\\\tau$ / ms\")\n",
    "plt.ylabel(\"$\\langle P(S) \\\\rangle$\")\n",
    "plt.legend()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We see that the spin echo increases the coherence time in comparison to free\n",
    "decay. Also, the coherence time matches the order of magnitude by the\n",
    "theoretically predicted (c.f., above mentioned reference)\n",
    "\n",
    "\\begin{align}\n",
    "T_{2,echo}^{1+\\alpha} = \\frac{1}{S_0 (dJ/d\\varepsilon)^2 (2-2^{-\\alpha+1})\n",
    "\\Gamma(-1-\\alpha) \\sin(\\pi \\alpha/2)}\n",
    "\\end{align}\n",
    "\n",
    "for $\\alpha = 1$ leading to $T_{2,echo} \\approx 137~\\text{ms}$."
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}