SIMPSON_COLAB.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5492a77b",
   "metadata": {},
   "source": [
    "<a href=\"https://githubtocolab.com/alsinmr/ConferenceExamples/blob/master/SIMPSON_COLAB.ipynb\" target=\"_blank\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32db91f0",
   "metadata": {},
   "source": [
    "# SIMPSON Template\n",
    "Simple template for creating a SIMPSON script, running it, and plotting the results using Jupyter notebooks, including installation and running of SIMPSON in Google Colab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9ffb32f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Setup for Google Colab\n",
    "import sys\n",
    "if 'google.colab' in sys.modules:\n",
    "    from urllib.request import urlretrieve\n",
    "    import os\n",
    "    url=\"https://inano.au.dk/fileadmin/_migrated/content_uploads/Simpson_Setup_Linux_4.1.1b.tar.bz2\"\n",
    "    filename=os.path.split(url)[1]\n",
    "    out=urlretrieve(url, filename)\n",
    "    !tar xf Simpson_Setup_Linux_4.1.1b.tar.bz2\n",
    "    %cd Simpson\\ Setup\\ Linux\\ 4.1.1\n",
    "    !./install.sh\n",
    "    !cd ..\n",
    "    !simpson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2d94bcaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "from subprocess import Popen, PIPE\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "bca549a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Class to read and run a SIMPSON script\n",
    "class RunRead():\n",
    "    def __init__(self):\n",
    "        self.stdout=None\n",
    "        self.stderr=None\n",
    "    def __call__(self,script):\n",
    "        with open('SIMPSON.tcl','w') as f:f.write(script)\n",
    "        process = Popen(['simpson', 'SIMPSON.tcl', '-d'], stdout=PIPE, stderr=PIPE)\n",
    "        self.stdout, self.stderr = [x.decode() for x in process.communicate()]            \n",
    "\n",
    "        if len(self.stderr):print(self.stderr)\n",
    "        # Collect the results\n",
    "        with open('SIMPSON.fid','r') as f:\n",
    "            while True:\n",
    "                if 'DATA'==f.readline().strip():break\n",
    "            S=[]\n",
    "            while True:\n",
    "                line=f.readline()\n",
    "                if 'END'==line.strip():break\n",
    "                S.append(complex(line.strip().split()[0])+complex(line.strip().split()[1])*1j)\n",
    "        S=np.array(S)\n",
    "        return S\n",
    "run_read=RunRead()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9ab3656f",
   "metadata": {},
   "outputs": [],
   "source": [
    "script=\"\"\"\n",
    "spinsys {\n",
    "  channels 13C 15N\n",
    "  nuclei   13C 15N\n",
    "  dipole   1 2 895 10 20 30\n",
    "  shift    1 10p 1p 0.5 50 20 10\n",
    "}\n",
    "\n",
    "par {\n",
    "  variable index   1\n",
    "\n",
    "  np               128\n",
    "  spin_rate        10000\n",
    "  proton_frequency 400e6\n",
    "  start_operator   I1x\n",
    "  detect_operator  I1p\n",
    "  method           direct\n",
    "  crystal_file     rep144\n",
    "  gamma_angles     100\n",
    "  sw               spin_rate/2\n",
    "  variable tsw     1e6/sw\n",
    "  verbose          1101\n",
    "  variable rfF1    75000\n",
    "  variable rfF2    50000\n",
    "  variable t180F1  0.5e6/rfF1\n",
    "  variable t180F2  0.5e6/rfF2\n",
    "  variable Del_t180 0.5*(t180F2-t180F1)     \n",
    "  variable tr2     0.5e6/spin_rate-t180F2\n",
    "  variable tr1     0.5e6/spin_rate-0.5*t180F2\n",
    "}\n",
    "\n",
    "proc pulseq {} {\n",
    "  global par\n",
    "  store 4\n",
    "  store 5\n",
    "\n",
    "  reset\n",
    "  delay $par(tr2)\n",
    "  pulse $par(t180F2) 0 x $par(rfF2) x\n",
    "  delay $par(tr2)\n",
    "  pulse $par(t180F2) 0 x $par(rfF2) y\n",
    "  store 1\n",
    "\n",
    "  reset\n",
    "  pulse $par(t180F2) 0 x $par(rfF2) y\n",
    "  delay $par(tr2)\n",
    "  pulse $par(t180F2) 0 x $par(rfF2) x\n",
    "  delay $par(tr2)\n",
    "  store 2\n",
    "\n",
    "  reset\n",
    "  delay $par(tr1)\n",
    "  pulse $par(Del_t180) 0 x $par(rfF2) x\n",
    "  pulse $par(t180F1) $par(rfF1) x $par(rfF2) x\n",
    "  pulse $par(Del_t180) 0 x $par(rfF2) x\n",
    "  delay $par(tr1)\n",
    "  store 3\n",
    "\n",
    "  acq\n",
    "  for {set i 1} {$i < $par(np)} {incr i} {\n",
    "    reset\n",
    "    prop 5\n",
    "    prop 2\n",
    "    store 5\n",
    "\n",
    "    reset\n",
    "    prop 4\n",
    "    prop 1\n",
    "    store 4\n",
    "\n",
    "    prop 3\n",
    "    prop 5\n",
    "    acq\n",
    "  }\n",
    "}\n",
    "proc main {} {\n",
    "  global par\n",
    "\n",
    "  set f [fsimpson]\n",
    "  fsave $f $par(name).fid\n",
    "}\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1db28b6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "S=run_read(script)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4e659b4c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmUAAAGsCAYAAACcgNDGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/v0lEQVR4nO3dd3yV5cHG8evOBkISIAlksWfYEIbgwLpw4nhrUetqrVq1dr6tbfXtsm/tsNXW9Vq3rdpqFVFRrFYQ2RkghBmyE0L2Ijvnfv9IaBHZJOd5zjm/7+fDh5yR5NInJ+fiue/nvo21VgAAAHBWkNMBAAAAQCkDAABwBUoZAACAC1DKAAAAXIBSBgAA4AIhTgc4VbGxsXb48OFOxwAAADimjIyMSmtt3OEe8/lSNnz4cKWnpzsdAwAA4JiMMQVHeozhSwAAABeglAEAALgApQwAAMAFKGUAAAAuQCkDAABwAUoZAACAC1DKAAAAXIBSBgAA4AKUMgAAABeglAEAALgApQwAAMAFKGUAAAAuQCkDAABwAUoZAACAC4Q4HQBwUnNbp4prmhQSHKSwkCCFBhtFRYQqIjTY6WgAgABDKUNAsdYqp7xRK3dVaOWuCq3Pq1Zbh+czzwkLCdIZo2N1wcQhOmdCvAZFhjuUFgAQSChlCBh7Khr14ze2aF1utSRp7OBI3XjaME1Kipa1UluHR22dHu2paNT72fv04Y5yBRlp/uhY/WDheE1Kinb4vwAA4M+MtdbpDKckLS3NpqenOx0DLtbS3qnHV+zR4yv2KCI0SHefM0YXTk5QUkyfI36OtVbb9tZr+dYy/WV9oWqa2nTF9CT99wXjlBB95M8DAOBojDEZ1tq0wz5GKYM/21xUq2//bZNyK/fr8mmJ+vHFqYrrf2LDkfUt7Xr0oxw9+0m+goKkOxaM1p1nj1ZwkOml1AAAf0UpQ0Ban1ulrzy3UTF9w/TAVZN1xpi4U/p6RdVNeuDdHXpny16dMSZWDy+eroH9wnooLQAgEBytlLEkBvzS6pxK3fTsRg2JjtAbd8w75UImSSkD++rR62bogSsna31utS790yfaWlLXA2kBAKCUwQ+t2Fmurzy3UcMG9dXfbjtN8VERPfr1F88eqldvP03WWl31+Br9I6O4R78+ACAwUcrgV1buqtCtL2RodHykXv7aXMX20nIWU1Ni9NY3TtfMYQP03Vc36+lP8nrl+wAAAgelDH6jpLZZd7+cpVHxkXrplrka0MvzvQZFhuu5m2frwklD9Iu3t+mxFTm9+v0AAP6NUga/0N7p0TdeylSnx+rx62Youm+oV75vWEiQ/nTNdC2alqjfvLdTf/jnLvn6xTMAAGeweCz8wu/e36nMwlo9cu10DY/t59XvHRIcpN9fPU1hwUF6+MPdau/06PsLx3s1AwDA91HK4PM+2lGu/1uZq+vmDNUlUxIdyRAcZPTrq6YoNCRIj63Yo4H9wnTLGSMdyQL3staqvrlDVftbVbW/TVWNrer0SJERIYoMD1H/iBANiY5QVIR3zvQCcBevlTJjzDOSLpFUbq2ddJjHjaSHJV0kqUnSTdbaTG/lg2/aW9es7/x9kyYkROm+S1IdzRIUZHT/okmqa2rX/e9sV1z/cC2aluRoJjirqa1DmwprlVFQo8zCGmUV1aq2qf2Yn5c8oI9SE6I0ISFKs0cM1OwRAxUazGwTwN9580zZc5IekfTCER6/UNKY7j9zJD3e/TdwRPct2arWDo8evXa6IkKDnY6joCCjB6+eqsrGVn3v1c2KjQzX/NGxTseCF3V0erR6T5WWZJVoeXaZmto6JUlj4iO1cOIQjY6PVGxkuAb2C9PAfmEKCTba39qhhpauP4XVTdq2t17b99brn9v3yVopKiJE50wYrPNTB+uscXHqG8YgB+CPvPbKttZ+bIwZfpSnLJL0gu2aJb3OGBNjjEmw1u71TkL4mhU7y/XB9nL98MLxGhkX6XScf4sIDdaTN6Tp6ifW6rYXM/T3205TamKU07HQyyobW/XMJ3l6NaNYFQ2tiooI0aJpSbpg4mBNHzpA0X1OfEhyf2uHVudUann2Pn24Y5/eyCpRVESIFs8equvnDlPKwL698F8CwCle3Wapu5S9fYThy7clPWCt/aT79oeSfmCt/dweSsaYWyXdKklDhw6dWVBQ0Ku54T5tHR4tfPhjWSst/9aZCgtx39DO3rpmXfnYGnms1dK7TtfgHl7EFu5Q3tCiJ1fm6q/rC9XS0anzJgzWlTOSdfb4OIWH9NzZ245OjzbkVeuvGwr13tYyWWt1Xupg3XrmKM0cNqDHvg+A3nW0bZbcdA78cLs7H7YxWmuflPSk1LX3ZW+Ggjs9vyZfuRX79exNs1xZyCQpIbqPnr15lq58bI1ufTFDf7t1riuGWNEzGls79McPd+v5Nflq7/To8ulJuvPs0RrVS2dtQ4KDNG90rOaNjlVpbbP+sq5AL28o1PLsfTovdbC+f8E4jRncv1e+NwDvcNO7WbGklINuJ0sqdSgLXKyioVV//HC3zh4Xp7PHxzsd56jGD4nS76+eps1FtfrRG1tYw8wPWGu1bMtenfvgSj35ca4umZKof313gX5/9bReK2SHSozpo+8vHK/V93xB3zt/rNbtqdIFD32s77+2WXvrmr2SAUDPc1MpWyrpBtNlrqQ65pPhcH67fIdaOjodv9ryeC2cNETfOneMXs8sYTsmH1dY1aSbnt2oO/6aqQH9wvT6HfP04NVTvb423gF9w0J01xfGaOX3z9bN80doSVapzn1wpZ5dnadOD/8AAHyNN5fEeFnSAkmxxphiST+RFCpJ1tonJC1T13IYOepaEuNmb2WD79hcVKu/pxfrtjNHumpy/7Hc/YUx2lnWoP9dtl1jBvfXWWPjnI6EE2Ct1T8yS/Q/b25VkDH6n0tSdcNpwxTikmUqBvYL032XpOrG04br3je36mdvbdOSrBL96sopXGQC+BCvTvTvDWlpaTY9/XPXAsBPXf/0em0rrdeK/16g/j62wOb+1g5d9fgaldY26527z+DKOR/R0NKue5ds1ZubSjV7xEA99KVpSozp43SsI7LWaunmUv38rW2qbW7X7WeN1DfPGevauZdwH4/HqrSuWbkV+1Xe0KqqxlZV729T1f42tXd65LGSx1pZa9U3LERREaGK6tP1d0J0hFIG9lXKwL4ndcVxIPCVif7AUX1aXKtVuyv1g4Xjfa6QSVK/8BA9eX2aLv7TKt35UqZevf20Hr06Dz1vU1GtvvFypkprW/Td88bqjrNHKzjocNckuYcxRoumJemssXG6/53tevSjPfp4V6X+8KVpGh3vO2eX4R3WWuVW7teGvGql59do57565Vbs//f6egeEhQRpUL8whYcEKcgYGdP1s9bc1qn65nY1tHZ87mtHRYRo7OD+mpgYpYmJ0UpNjNK4If1ZCPkoOFMGn3Hbi+lau6dKq+/5gk+WsgPe21qm2/+SoRtOG6afL/rc6jBwidcyivWj17corn+4Hl48TWnDBzod6aS8t3Wv7nl9i1raO3Xvxam6bs5QdW2ggkBV19Suf+3cpw+2lWt9XpUqG9skSYP6hSk1MUqj4yM1Oj5So+IilRjdRwMjw9QvLPioPzedHqv65naV1DaruKZJRdXNKqjerx17G7R9b732d5e8vmHBmjF0wL93qpg+NCbg/nHKmTL4vF37GrQ8e5/u/sJony5kUtfE/1tOH6GnPsnTrOEDdelUZ/brxOF1dHr0wLs79NQneZo3apAevXaGBvQLczrWSVs4KUHThw7Q917drHuXbNXKXRX63RenMrQUYOqa2rX001It31qmdblV6vBYxfcP15lj4zR7+EDNGjFQI2P7nXRhDw4yGtAvTAP6hWlSUvRnHvN4rAqqm7SlpE4Z+dVan1etP3ywS9Z2lbR5o2K1YFycFoyLU/KAwJ7WwZky+IRv/22TlmeXafUPvuDTb5AHtHd6tPjJddqxt15Lv3G615ZSwNHVNbXrrpcztWp3pW6aN1w/vniC3wy1eDxWz6zO0wPv7lBiTB89dt2Mz715wr94PFbrcqv0t/Qivbu1TG0dHo2M66cLJg7R+amDNTU5RkEODcfXNrVpQ161Pt5doRU7K1Rc07WUS2pClC6ekqCLJidohENXNfe2o50po5TB9QqrmnT2gyt087zhutdHlsE4HnvrmnXRw6sU3z9CS+6crz5hgXUK321Kapt1w9PrVVjdpJ8vmqRrZg91OlKvyCio0V0vZapqf5t+ftlEfWlWCsOZfqaprUOvZRTrmU/ylF/VpKiIEF0+PUlXp6W4sohba7WnYr8+2lGud7fuVWZhrSRpQkKULp48RBdNTvCpq+2PhVIGn/bD17foHxnFWvWDs/1uq6KVuyp007MbdNWMZP3ui1OdjhOwdpTV68ZnNqiprVN/viFNc0cOcjpSr6pqbNW3/rZJq3ZX6oszk3X/FZMCbl6PP6poaNULa/P14roC1Ta1a/rQGN142nAtnDTEp3YTKa1t1rIte7Vsy38K2vgh/XXx5ARdOjXRsXUBewqlDD6rrK5FZ/7mI30xLVm/vGKy03F6xe/f36k//itHv7lqiq6elXLsT0CPWp9bpVteSFffsGA9/5XZGj8kMNb16vRYPfzBLv3xXzmaMTRGT1w/U/H9/esfPYGisrFVT6zYoxfXFait06PzUwfr1jNHauYw37w45WCltc16d2uZlm3Zq4yCGknS1JQYXT4tURdPSfDJn1lKGXzWr97drqdW5emj7y7Q0EH+OQG002N1wzPrlZ5fozfumM9in160PLtM33g5S8kD+uiFr8wOyEnGy7bs1Xf/vlkxfUP15PVpmpzsvuEtHF5dU7ueXLVHz67OV0t7p66Ynqw7zx7lV0N9B9tb16y3NpdqSVaptu2tV5CR5o+O1aJpSbpg4mCfuQiMUgaf1NrRqbn/+6HmjBikJ66f6XScXlXR0KqL/7hK/cJDtPSu+T7zy8WXLckq0Xdf3azJSdF69qZZfnEBycnKLq3TrS9kqLKxVb/94lRdxhXBrtbW4dELa/P1xw93q76lQ5dOTdQ3zxkTUOvQ7d7XoDc3lerNzSUqqm5WeEiQzp0wWIumJeqscXGuHo6nlMEnvbmpRN98ZZNe/OpsnTHG/7clWp9bpWufWq8LJg7Wo9fOYPJ1L3p5Q6F+9MYWzR0xSE/dmKZ+4awOVNnYqq//JUMb82t0x4JR+t754xy7Mg+HZ63V8ux9+tW721VQ1aSzxsbpngvHa0JC4J5dt9Yqs7BWSzeV6O1P96pqf5ui+4TqoslDdNnUJM0ZMdB1P8eUMvikq59Yq30NLfrouwtc96LqLY+v2KNfv7dDP7k0VTfPH+F0HL/0zCd5+vnb27RgXJye+PJMn5oA3dvaOjz6ydKtenlDkc6dEK8/fGkaZ21dYmdZg/7nza1an1etsYMj9aOLJmjBuHinY7lKe6dHq3Mq9eamUi3PLlNTW6eGREXosmmJWjQtUakJUa74xy6lDD5n174Gnf+Hj/XDC8frtrNGOR3Hazweq6+9kK6Pd1fo77edpulDBzgdya8cKL0XTBysP14z3dVDHE6x1urFdQX62VvbNDK2n/58Q5rPX+3myxpbO/TwB7v0zOp89Y8I0ffOH6fFs1IU4ifr5/WW5rZO/XP7Pi3dVKIVOyvU4bEaEx+pCycn6IKJgx0taJQy+JyfLs3WS+sLte5H52hggM31qW1q0yV/+kQej9U7d58R0HOdetITK/fogXd36LKpifr91VN5UzuGNTmVuuOlTEnSY9fN0LxRsQ4nCizWWr27tUw/f2ubyupbtHhWir6/cHzA/T7sCTX727Rs6169ualUG/OrZa2UPKCPLpg4ROeMj1fa8IEKC/He7wNKGXxKU1uH5vzvhzpnfLweWjzd6TiO+LS4Vv/1+FrNGz1Iz9w4K2CGb3vLkx/v0f8uo5CdqIKq/frq8+nKr9yvny2aqOvmDHM6UkDYV9+ie5ds1T+37VNqQpR+cfkkzRzGWfOeUNnYqg+27evaISanSm2dns9s9XRe6uBeXw+TvS/hU97evFcNLR26bm7gvgFMSY7RfZdM0H1vZuvxlXt059mjnY7ks55alav/XbZDl0xJoJCdoGGD+un1O+bp7pez9OM3tmr3vkbde/EE/h/2Emut/p5epPvf2a62Do9+eOF4ffX0Efz/7kGxkeFaPHuoFs8eqv2tHVqzp0ord5Vrxc4KfbB9n0KCjBY7uJsHpQyu85f1BRo7OFJpAf4vwy/PHaYN+TV68P2dmj40huGjk/D0J3m6/53tunhygh760jTe3E5CVESonr5xln61bLue+iRPeyoa9cg1MxTdlwsAelJhVZPuef1TrdlTpTkjBuqBq6b47d6PbtEvPETnpQ7WeamDZa1VbuV+xfYLdzQTv6HgKp8W1+rT4jpdN2eYK66ScZIxRr+6crKGx/bT3S9vUnl9i9ORfMqzq/P0i7e36cJJQ/TQYgrZqQgOMrr3klT95qopWpdbpSseW63cikanY/mFTo/VU6tydf5DK/VpcZ1+ecUkvfy1uRQyLzPGaFRcpOP/2OC3FFzlbxuLFBEapMunJzkdxRUiw0P0xJdnan9rh+56OUsdnR6nI/mE59fk62dvbfv3VZahFLIecfWsFP31lrmqbW7X5Y+u1ie7K52O5NN2ljXoqsfX6P53tmveqFj98ztn6ro5w5hDGsD4TQXXaOvw6J0te3V+6hBF92Fo5ICxg/vrl1dM0oa8aj34z11Ox3G9F9bm6ydLs3V+6mD96ZoZFLIeNnvEQL1553wlRPfRjc9u0HOr8+TrF4x5W1uHRw99sEuX/GmVCqr26+HF0/T0jWlKiO7jdDQ4jDllcI2Pd1Wotqldi6axxcuhrpyRrI351Xp8xR6lDRugcyYMdjqSK724rkD/82a2zp0wWI9cO8Orl7kHkpSBffWPO+bpW69s0k/f2tY97DZZfcJY9+1YNhfV6vuvfaqd+xp02dRE/eTSVA2KdHYeE9yD31hwjSWbSjSgb6jOHOv/WyqdjJ9cOlGpCVH69t82Ka9yv9NxXOfFtfm6b8lWnTM+Xo9dRyHrbZHhIXry+pn6znlj9camEl31+BoVVTc5Hcu1mto69Mt3tumKx1arrrldT92Qpj9eM51Chs/gtxZcobG1Qx9s36eLpyQw3HQEEaHB+r/rZyo4yOhrL6SroaXd6Uiu8cLafN33ZrbOnRCvx75MIfOWoCCju88Zo2dunKXimiZd+sgnWrmrwulYrvPxrgqd/4eP9edVefrSrBS9/50zdW4qZ7vxefzmgiss31qmlnaPLp/GBP+jSRnYV49eN0N5lfv1rVc2qdPDXJ4X1ub/e8jysetmsnWSA84eH6+ld52uIVERuvGZDfrNezu4KEVSVWOrvv23TbrhmQ0KCw7SK7fO1a+unKIo9hPFEVDK4ApLNpUoeUAfVq0+DvNGxeonl6bqwx3levD9nU7HcdTza7oK2XmpgxmydNjw2H564475WjwrRY+t2KNr/rxOe+uanY7lCGut/pFRrHN/v1Jvf1qqu78wWsu+eYbmjhzkdDS4HL/B4LiKhlatzqnUommJAb822fG6fu4wXTO7681v6eZSp+M44rnVef++yvJRJvW7Qp+wYD1w1RQ9vHiatpXW66KHV+mDbfucjuVVhVVNuuGZDfruq5s1Iraf3rn7DH3n/HGKCOUMLo6Nqy/huLc/LZXHiqHLE2CM0c8um6Sc8kb996ublRgdobThA52O5TXPrs779zpkf7qGQuY2i6YlaXJStO58KUu3vJCuq9OSdd8lqervx8N2bR0ePbM6Tw99sEshQUH6xeWTdN3soaw5hhPCbzI4bsmmUqUmRGnM4P5OR/EpYSFB+r/r05QY00e3vJCuPQGywvozn/ynkLHshXuNjIvUkjvn6esLRum1jGItfGiV1uzxz8VmV+ws18KHPtYD7+7QmWPi9MF3ztL1c1kEFieO32ZwVF7lfm0uqtXl01mb7GQM7Bem52+erZAgoxuf2aDyBv/eiumpVbn6+dvbtHDiED1yLQvDul14SLB+sHC8Xr39NIUGG1375/W6b8lW1TX7x5XDBVX7dcvz6brp2Y2SpGdvnqUnb0jTkOgIh5PBV/EbDY5auqlUxkiXTqWUnayhg/rq6RtnqaqxTV95bqP2t3Y4HanHWWv1+/d36v53tuuiyUP0p2vZOsmXzBw2UMu+eYZumjdcf11foHMeXKklWSU+uxNAZWOrfro0W+f+fqXW7qnUDy8cr/e+dabOHhfvdDT4OOOrL4oD0tLSbHp6utMxcJIufHiV+oUF67Wvz3M6is/7aEe5bnkhXfNGDdJTN6b5zdIQHo/Vz97K1vNrC/SltBT98opJbC7uw7YU1+neJVu0ubhOp40cpJ9eNlHjhvjG1IXG1g79+eNcPbUqVy0dHl2dlqJvnztG8VGcGcPxM8ZkWGvTDvcYv9ngmKLqJm3fW68LJg5xOopfOHt8vH515WSt2l2p21/MUGtHp9ORTll7p0ff+fsmPb+2QLeeOVIPXDWZQubjJidH6/U75uuXV0xSdmmdFj78sb79t00qqHLvLhW1TW16+IPdOuPX/9LDH+7WWePi9P63z9SvrpxMIUOP4upLOGZ5dpkkUcp60NVpKerotPrRG1v09b9k6vEvz/DZM2YNLe2666UsrdxVoe8vHKevnzWKJVP8RHCQ0XVzhumiSQl64uM9en5Nvt7aXKqrZ6XozrNHKynGHRtz761r1tOr8vTShkI1tXXq3Anx+sYXxmhqSozT0eCnGL6EY774xBo1tHTovW+d6XQUv/PX9QX68Rtbu7Ydum6mz12hWFTdpK8+v1G5Fft1/+WTtHj2UKcjoReV17fokY9y9PKGQnV6rM5PHaIb5w3X3JEDvV7EOzo9WrmrQi9vKNJHO8slSZdNTdRtZ43U+CFRXs0C/3S04UvOlMERFQ2tSi+o0d1fGON0FL903Zxh8ljpviVb9fW/ZOiRa2eoT5hvnDHLKKjWrS9kqL3To+e/MlvzR8c6HQm9LD4qQj9fNEm3nzVKL64r0CsbCvVedpnGDe6vL6Yla+GkIUoe0LfXvr/HY/VpSZ3ezy7T65klKqtvUWxkuL52xkhdN2eoUgb23vcGDsaZMjji5Q2F+uHrW7Ts7jOUmsi/PnvLX9YV6L43t2pKcoyevjFNsZHhTkc6qtczi3XP61uUGB2hp2+apVFxkU5HggNa2ju1dHOpXlxboC0ldZKkyUnRWjhpiM4YE6vxQ6JO+exvZWOrsgpr9dHOcn2wbZ/KG1oVZKQzx8Zp8ayhOmdCPFf4olcc7UwZpQyOuOnZDdpT0aiP//ts5gn1svezy3T3K1mK6x+u526e7cqi09LeqZ8uzdYrG4s0Z8RAPfHlmRrQL8zpWHCBvMr9Wp5dpve2lmlTUa0kKTwkSJOSojUtJUaj4iKVEB2hwVERGhIdoT6hwbKyslaykmr2t6m0tll761pUUtusbXvrtbmoVsU1Xfty9g0L1llj43Re6mCdPS6enzv0OkoZXKWhpV0zf/GBbjhtmO69JNXpOAFhU1Gtbnl+ozo8Vv/35Zma46KNkfdUNOrOv2ZqR1mD7lgwSt85byxXWOKwyupalFFQo6zCGmUV1WprSZ1aOzwn9DWSYvpoWkqMpqZEa1rKAE1JjmZfSngVc8rgKh/trFBbp0cXTOKqS2+ZlhKj178+Xzc9t0HXPrVe3zpnjO44e7SCHdwGxlqr1zKK9ZOl2QoPCdJzN8/SAhbfxFEMiY7QxVMSdPGUBEldS6aUN7SqrK5F++pbVFbXorbOrpJmJBkjxfQNU0J0hBKi+yghOkL9wnnbg3vx0wmvW55dptjIMM0YOsDpKAFl6KC+WnLnfN23ZKse/Ocufby7Qn/40rRenUB9JHmV+/XjN7ZozZ4qzR4+UA9fM00J0e5YBgG+IzQ4SEkxfVyzhAZwqhgjgFe1tHdqxY5ynZc62NGzNIEqKiJUDy+erj98aaq2723QhQ+t0msZxfJ4vDONoa3Doz99uFsXPPSxtpTU6f7LJ+mVW+dSyABAnCmDl63ZU6n9bZ06nwVjHXXF9GSlDRuob/1tk7736mY9tyZPP7xwQq8tP9HW4dHrmcV6dEWOiqqbdfGUBP3kklRWQweAg1DK4FX/3FaufmHBmjfKPRPNA1XKwL569bbTtGRTiR58f5eue2q9zhwbp++eN1ZTkqN75KrY1o5OvZZRrMc+2qOS2mZNSY7W/ZdP1llj43rgvwAA/AulDF5jrdXKneU6fUysz27942+CgoyunJGsiyYn6MW1BXrkoxwtenS1xg6O1BXTk3X59MQTHlps7/RodU6l3vl0r5Znl6m+pUPTUmJ0/xWTtGBsHEugAMARUMrgNbv2Naq0rkV3n8Mq/m4TERqsr505UlfPStFbm0v1emaxfv3eDv1m+Q5NSY7RxMQoTUiIUmpCfyXF9FWQkWQkI6Oapjbt2tegXWUN2rWvUWtzq1TX3K7+4SE6L3WwrpyRrPmjB1HGAOAYKGXwmgP7yJ01jqErt4ruE6ovzx2mL88dpvzK/Xojq0Rrc6v01uZSvbS+8KifG2Sk4YP66Zzx8bpocoLOGMsZUQA4EZQyeM2KneUaP6Q/V9r5iOGx/fTt88bq2+oaei6pbdb2vQ0qb2j592rpslaRESEaO7i/RsVFsggnAJwCShm8oqGlXen5NframSOdjoKTYIxR8oC+jqxpBgCBgnXK4BWrcyrV4bFawFV3AAAcFqUMXvHRjgr1jwjRjGGs4g8AwOFQytDrrLVasatcZ46JUygbTQMAcFi8Q6LXbd/boH31rVrAVZcAABwRpQy9jqUwAAA4NkoZet2KneWalBSl+P7scwgAwJFQytCr6pralVlYq7PHxTsdBQAAV6OUoVetyqlQp8cynwwAgGOglKFXfbyrQlERIZqWwlIYAAAcDaUMvcZaq9U5VZo3KlbBQWxGDQDA0Xi1lBljFhpjdhpjcowx9xzm8WhjzFvGmM3GmGxjzM3ezIeeVVDVpJLaZs0fE+t0FAAAXM9rpcwYEyzpUUkXSkqVdI0xJvWQp90paZu1dqqkBZIeNMaEeSsjetYnOZWSpNNHU8oAADgWb54pmy0px1qba61tk/SKpEWHPMdK6m+MMZIiJVVL6vBiRvSg1TmVSorpo+GD2MQaAIBj8WYpS5JUdNDt4u77DvaIpAmSSiVtkfRNa63n0C9kjLnVGJNujEmvqKjorbw4BZ0eqzV7qjR/9CB1dWwAAHA03ixlh3tntofcvkDSJkmJkqZJesQYE/W5T7L2SWttmrU2LS6OpRbcKLu0TnXN7ZrP0CUAAMfFm6WsWFLKQbeT1XVG7GA3S3rddsmRlCdpvJfyoQcdmE82bxSlDACA4+HNUrZR0hhjzIjuyfuLJS095DmFks6RJGPMYEnjJOV6MSN6yOqcSo0f0l9x/cOdjgIAgE/wWimz1nZIukvScknbJf3dWpttjLndGHN799N+IWmeMWaLpA8l/cBaW+mtjOgZLe2d2phfw1WXAACcgBBvfjNr7TJJyw6574mDPi6VdL43M6HnpefXqK3Dw/pkAACcAFb0R4/7JKdSocFGs4cPdDoKAAA+g1KGHrc6p1LTUwaoX7hXT8QCAODTKGXoUTX727S1tI6lMAAAOEGUMvSotblVslY6fcwgp6MAAOBTKGXoUWv2VKpfWLCmJMc4HQUAAJ9CKUOPWpdbrVkjBio0mB8tAABOBO+c6DGVja3KKW/UnBEMXQIAcKIoZegx63OrJUlzR7IUBgAAJ4pShh6zLrdK/cKCNSkp2ukoAAD4HEoZesy63CqlDWc+GQAAJ4N3T/SIysZW7S5v1ByGLgEAOCmUMvSIDXkH5pMxyR8AgJNBKUOPWJdbpb5hwZrMfDIAAE4KpQw9Yn1utWYOG8B8MgAAThLvoDhlVY2t2rmvgaFLAABOAaUMp4z5ZAAAnDpKGU7Zutwq9QkN1pRk5pMBAHCyKGU4ZevzqpU2nPlkAACcCt5FcUqq97dpRxnzyQAAOFWUMpySDXlVktjvEgCAU0UpwylZn1etiNAgTU6KcToKAAA+jVKGU7Ixv1rTUwYoLIQfJQAATgXvpDhpja0d2lZar1nDBzgdBQAAn0cpw0nLKqyRx0qzRjCfDACAU0Upw0nbmFetICNNH8qZMgAAThWlDCdtY36NUhOjFBke4nQUAAB8HqUMJ6Wtw6OsohrNGs7QJQAAPYFShpOSXVqnlnYPpQwAgB5CKcNJSc+vkSSlceUlAAA9glKGk7Ihv1rDB/VVfP8Ip6MAAOAXKGU4YdZapedXK42hSwAAegylDCdsT0WjapraNZtSBgBAj6GU4YRtZD4ZAAA9jlKGE7Yxv1qxkWEaEdvP6SgAAPgNShlO2Mb8aqUNGyhjjNNRAADwG5QynJCyuhYVVTczdAkAQA+jlOGEpBdUS5Jmswk5AAA9ilKGE5KeX6M+ocGakBDldBQAAPwKpQwnJKOgRtNSYhQazI8OAAA9iXdWHLemtg5t21uvmcOYTwYAQE+jlOG4bS6qU6fHUsoAAOgFlDIct8zCrkVjpw+NcTYIAAB+iFKG45aeX60x8ZGK6RvmdBQAAPwOpQzHxeOxyiysZegSAIBeQinDccmtbFRdc7tmUMoAAOgVlDIcl4yCrvlknCkDAKB3UMpwXDIKajSgb6hGsgk5AAC9glKG45JeUKOZwwawCTkAAL2EUoZjqt7fptyK/cwnAwCgF1HKcExZ3euTzRxKKQMAoLdQynBMGQU1CgkympIc43QUAAD8FqUMx5RRUKOJiVHqExbsdBQAAPwWpQxH1d7p0ebiWuaTAQDQyyhlOKptpfVqafcobdhAp6MAAODXKGU4qgOLxs4YFuNsEAAA/BylDEeVVVSrxOgIJUT3cToKAAB+jVKGo8oqrNF0lsIAAKDXebWUGWMWGmN2GmNyjDH3HOE5C4wxm4wx2caYld7Mh88qb2hRcU2zpg+NcToKAAB+L8Rb38gYEyzpUUnnSSqWtNEYs9Rau+2g58RIekzSQmttoTEm3lv58HmbCmsliVIGAIAXePNM2WxJOdbaXGttm6RXJC065DnXSnrdWlsoSdbaci/mwyGyimoVGmw0MTHa6SgAAPg9b5ayJElFB90u7r7vYGMlDTDGrDDGZBhjbjjcFzLG3GqMSTfGpFdUVPRSXGQV1ig1IUoRoSwaCwBAb/NmKTOHuc8ecjtE0kxJF0u6QNJ9xpixn/ska5+01qZZa9Pi4uJ6PinU0enRp8V1mpYS43QUAAACgtfmlKnrzFjKQbeTJZUe5jmV1tr9kvYbYz6WNFXSLu9ExAG79jWqqa2TKy8BAPASb54p2yhpjDFmhDEmTNJiSUsPec6bks4wxoQYY/pKmiNpuxczoltWUdeisUzyBwDAO7x2psxa22GMuUvScknBkp6x1mYbY27vfvwJa+12Y8x7kj6V5JH0lLV2q7cy4j82FdZqYL8wDR3Y1+koAAAEBG8OX8pau0zSskPue+KQ27+V9Ftv5sLnZRXVanpKjIw53FRAAADQ01jRH59T19yunPJGhi4BAPAiShk+Z3NRrSQxyR8AAC+ilOFzsgprZYw0JZlFYwEA8BZKGT4nq6hGY+P7q39EqNNRAAAIGJQyfIa1VlmFtcwnAwDAyyhl+Iy8yv2qa26nlAEA4GWUMnxGVmGtJCb5AwDgbSe8Tpkx5lpJl0nqVNd+lm9Za1/u6WBwRlZRjfqHh2h0XKTTUQAACCgns3jsWdbaxQduGGMelUQp8xNZhbWamhKjoCAWjQUAwJtOZvgy3BhzsTFmijHmIkl9ejoUnNHU1qEdZQ3MJwMAwAEnU8rukDRA0kXdf9/Vo4ngmC3Fder0WEoZAAAOOOHhS2ttk6S/9EIWOCyreyX/qckxjuYAACAQHfNMmTHmQ2PMxINuX2aMudcYM7t3o8HbsgprNGxQXw2KDHc6CgAAAed4hi+TrbXZkmSMmSfpRUlDJT1njLmiN8PBe6y1yiys1fSUGKejAAAQkI6nlNUf9PENkp6w1t4qaYGkH/RGKHhfaV2LKhpaWZ8MAACHHE8pyzHG/JcxJl7S5ZLelCRrbbkkxrn8xKZ/Lxob42gOAAAC1fGUsm9Luk1SiaQsa+0aSTLGhEpihVE/kVVYo/CQII0fEuV0FAAAAtLxXH05QtL5koy11nPQ/WdL+qhXUsHrsopqNTkpWmEh7LwFAIATjucd+EZJGZJeMsbcZIwZIknW2ve755bBx7V1eLSlpI6hSwAAHHTMM2XW2tslyRgzXtKF6rrqMlpdZ8nek7TaWtvZqynRq7bvrVdbh4dJ/gAAOOi4x6qstTustX+w1i6U9AVJn0j6oqT1vRUO3pFVWCOJSf4AADjpZDYkl7W2WdKy7j/wcVlFtRoSFaGEaLYxBQDAKaeyov+c3o0Gb8kqrOUsGQAADjuVFf2fZUV/31fZ2KrC6iZKGQAADmNF/wD3n0VjmeQPAICTWNE/wGUV1SgkyGhSYrTTUQAACGgnuqJ/Jiv6+5eswlpNSIhSn7Bgp6MAABDQjmedsjJJ5xljgg6s6N+9gOwUsaK/T+v0WG0uqtVVM5OdjgIAQMA77iUxDtliaZm1doak93s+Erxld3mD9rd1alpKjNNRAAAIeCe70aHp0RRwRBaT/AEAcI2TLWV/7tEUcERWYY1i+oZq+KC+TkcBACDgnVQps9Y+JknGmJt7Ng68aVNRraanxMgYTnwCAOC0kz1TdsDPeiQFvK6+pV27yxsZugQAwCWOOdHfGPPpkR6SNLhn48BbPi2qk7VsQg4AgFscz9WXgyVdIKnmkPuNpDU9nghekVVYI2OkqVx5CQCAKxxPKXtbUqS1dtOhDxhjVvR0IHhHVlGtRsdFKioi1OkoAABAx7d47FeP8ti1PRsH3mCtVVZhjc5LZfQZAAC3ONWJ/vBBBVVNqmlqZ5I/AAAuQikLQFlFXdMDmeQPAIB7UMoCUFZhrfqFBWtMfH+nowAAgG6UsgCUVVirqSkxCg5i0VgAANyCUhZgmts6tX1vPUOXAAC4DKUswGwtrVOHx2p6CpP8AQBwE0pZgMkq7JrkP40zZQAAuAqlLMBkFtRq6MC+io0MdzoKAAA4CKUsgFhrlVlYw3wyAABciFIWQErrWlTe0KoZLBoLAIDrUMoCSGZB13wyShkAAO5DKQsgWYW1iggN0vgEFo0FAMBtKGUBJLOwRlOSYhQazGEHAMBteHcOEC3tncourdP0YTFORwEAAIdBKQsQ2aX1au9k0VgAANyKUhYgDiwaO4MzZQAAuBKlLEBkFtYoeUAfxfePcDoKAAA4DEpZgMgqrGUpDAAAXIxSFgD21jVrb10LK/kDAOBilLIAkFlQK4lFYwEAcDNKWQDIKqxReEiQJiREOR0FAAAcgVdLmTFmoTFmpzEmxxhzz1GeN8sY02mM+S9v5vNXmYU1mpwUrbAQOjgAAG7ltXdpY0ywpEclXSgpVdI1xpjUIzzv15KWeyubP2vt6NTWknrNGMbQJQAAbubNUyezJeVYa3OttW2SXpG06DDP+4akf0gq92I2v7WttF5tnR7NYJI/AACu5s1SliSp6KDbxd33/ZsxJknSFZKeONoXMsbcaoxJN8akV1RU9HhQf5JZWCtJms4kfwAAXM2bpcwc5j57yO2HJP3AWtt5tC9krX3SWptmrU2Li4vrqXx+KbOwRkkxfTQ4ikVjAQBwsxAvfq9iSSkH3U6WVHrIc9IkvWKMkaRYSRcZYzqstUu8ktAPbSqsZX0yAAB8gDfPlG2UNMYYM8IYEyZpsaSlBz/BWjvCWjvcWjtc0muS7qCQnbx99S0qqW1m6BIAAB/gtTNl1toOY8xd6rqqMljSM9babGPM7d2PH3UeGU5cZkH3JuScKQMAwPW8OXwpa+0yScsOue+wZcxae5M3MvmzrKJahYUEaWJitNNRAADAMbCaqB/LLKjRpMQoFo0FAMAH8G7tp9o6PPq0pI79LgEA8BGUMj+1fW+92jo8rOQPAICPoJT5qczCA5P8KWUAAPgCSpmfyiysVUJ0hIZEs2gsAAC+gFLmp7IKazhLBgCAD6GU+aHyhhYV1zSzkj8AAD6EUuaHMgtqJbEJOQAAvoRS5oeyimoUFhykSUlRTkcBAADHiVLmh7IKajUxKUrhIcFORwEAAMeJUuZn2js9+rSkVtNTGLoEAMCXUMr8zI69DWpp92jGsBinowAAgBNAKfMzBxaNZZI/AAC+hVLmZ9ILapQQHaGkmD5ORwEAACeAUuZnMvKrNZP9LgEA8DmUMj9SUtus0roWpVHKAADwOZQyP5KeXy1JShs+0OEkAADgRFHK/EhGQY36hgVr/JD+TkcBAAAniFLmR9LzuzYhDwnmsAIA4Gt49/YTja0d2lFWzyR/AAB8FKXMT2QV1shjpbThlDIAAHwRpcxPpOfXKMiwaCwAAL6KUuYnMgpqNH5IlCLDQ5yOAgAATgKlzA90dHqUWVjD0CUAAD6MUuYHdpQ1qKmtk0n+AAD4MEqZHziwaOwsFo0FAMBnUcr8QHpBjRKjI5TIJuQAAPgsSpkfyCio0UzOkgEA4NMoZT6upLZZe9mEHAAAn0cp83EH5pMxyR8AAN9GKfNxG/KqFRkewibkAAD4OEqZj9uQV6204WxCDgCAr+Od3IdVNbZqd3mj5owY5HQUAABwiihlPmxj93yy2SO48hIAAF9HKfNh63Kr1Sc0WJOTop2OAgAAThGlzIetz6vWjGExCgvhMAIA4Ot4N/dRdU3t2lFWz3wyAAD8BKXMR23Mr5a1zCcDAMBfUMp81Ib8aoUFB2laSozTUQAAQA+glPmo9blVmpYSo4jQYKejAACAHkAp80GNrR3aWlqvOSMZugQAwF9QynxQRkGNOj2W+WQAAPgRSpkP2pBXpZAgwybkAAD4EUqZD1qfW61JSdHqGxbidBQAANBDKGU+pqW9U5uLa5lPBgCAn6GU+ZjMwhq1d1rNZdFYAAD8CqXMx6zLrVaQkWYOZz4ZAAD+hFLmY9buqdTkpGhFRYQ6HQUAAPQgSpkPaWrrUFZhrU4bFet0FAAA0MMoZT5kY36NOjxW80cznwwAAH9DKfMha3IqFRpslDaMKy8BAPA3lDIfsmZPlaYPHaA+Yex3CQCAv6GU+Yi6pnZtLa3TvFEMXQIA4I8oZT5iXV6VrJXmj2aSPwAA/ohS5iPW5FSqT2iwpibHOB0FAAD0AkqZj1izp0qzRgxUWAiHDAAAf8Q7vA8ob2jR7vJG5pMBAODHvFrKjDELjTE7jTE5xph7DvP4dcaYT7v/rDHGTPVmPrdau6dKkihlAAD4Ma+VMmNMsKRHJV0oKVXSNcaY1EOelifpLGvtFEm/kPSkt/K52ZqcKkVFhGhiYrTTUQAAQC/x5pmy2ZJyrLW51to2Sa9IWnTwE6y1a6y1Nd0310lK9mI+11qTW6m5IwcpOMg4HQUAAPQSb5ayJElFB90u7r7vSL4q6d3DPWCMudUYk26MSa+oqOjBiO5TVN2koupmhi4BAPBz3ixlhzvNYw/7RGPOVlcp+8HhHrfWPmmtTbPWpsXFxfVgRPdZs6dSkjSP9ckAAPBrIV78XsWSUg66nSyp9NAnGWOmSHpK0oXW2iovZXOtlbsqNDgqXGPiI52OAgAAepE3z5RtlDTGGDPCGBMmabGkpQc/wRgzVNLrkq631u7yYjZX6uj0aNXuSp01Nk7GMJ8MAAB/5rUzZdbaDmPMXZKWSwqW9Iy1NtsYc3v3409I+h9JgyQ91l1COqy1ad7K6DaZhbVqaOnQgnHxTkcBAAC9zJvDl7LWLpO07JD7njjo41sk3eLNTG62cle5goMM+10CABAAWNHfxVbsrNDMoQMU3SfU6SgAAKCXUcpcqryhRdml9TprnH9fXQoAALpQylxq5c6u9dcWUMoAAAgIlDKXWrGrQvH9w5WaEOV0FAAA4AWUMhfq6PRo1a4KlsIAACCAUMpcaFNRrepbOphPBgBAAKGUudDKXRUKMtIZoyllAAAECkqZC63YWaEZQwcoui9LYQAAECgoZS5T0dCqLSV1XHUJAECAoZS5zMe7DiyFwdZKAAAEEkqZy/xz2z4NjmIpDAAAAg2lzEVa2ju1cleFzksdrKAglsIAACCQUMpc5JPdlWpu79QFE4c4HQUAAHgZpcxF3t9Wpv4RIZozYpDTUQAAgJdRylyi02P1wfZyfWF8vMJCOCwAAAQa3v1dIqOgRtX723R+KkOXAAAEIkqZSyzPLlNYcBBbKwEAEKAoZS5grdX728o0f/QgRYaHOB0HAAA4gFLmAjvKGlRU3azzueoSAICARSlzgfez98kY6ZwJrOIPAECgopS5wPvbyjRj6ADF949wOgoAAHAIpcxhxTVNyi6t1/mpg52OAgAAHEQpc9h7W8skiflkAAAEOEqZw97IKtHkpGiNiO3ndBQAAOAgSpmDdpY1KLu0XlfOSHI6CgAAcBilzEGvZxUrJMjo0qmJTkcBAAAOo5Q5pNNjtSSrRGeNjVNsZLjTcQAAgMMoZQ5Zu6dK++pbdeWMZKejAAAAF6CUOeT1zGL1jwhhwVgAACCJUuaI/a0dei+7TJdMSVBEaLDTcQAAgAtQyhywPLtMTW2dDF0CAIB/o5Q54I2sEqUM7KO0YQOcjgIAAFyCUuZlZXUt+iSnUldMT5Yxxuk4AADAJShlXvaPzGJZK10xnQVjAQDAf1DKvKij06O/rivQ6aNj2VYJAAB8BqXMiz7Yvk+ldS264bRhTkcBAAAuQynzoufW5Csppo/OmTDY6SgAAMBlKGVesrOsQetyq3X9acMUHMQEfwAA8FmUMi95YW2+wkOC9KW0FKejAAAAF6KUeUFdc7tezyzRommJGtAvzOk4AADAhShlXvBaRrGa2zt1w2nDnY4CAABcilLWyzweqxfX5mvmsAGalBTtdBwAAOBSlLJetmJXufKrmnTjvOFORwEAAC5GKetFHo/V75bvUvKAPlo4cYjTcQAAgItRynrRm5tLtG1vvf77gnEKC+F/NQAAODKaQi9pae/U75bv0qSkKF06JdHpOAAAwOUoZb3kxbUFKqlt1g8vnKAgFosFAADHQCnrBXVN7XrkoxydOTZO80fHOh0HAAD4AEpZL3hsZY7qW9p1z8LxTkcBAAA+glLWw0pqm/Xs6nxdMT1JqYlRTscBAAA+glLWgzweq3v+8amCjPTd88c5HQcAAPgQSlkPem5NvlbtrtS9F6cqKaaP03EAAIAPoZT1kB1l9XrgvR06d0K8rpsz1Ok4AADAx1DKekBLe6e++fImRUWE6IGrpsgYlsAAAAAnJsTpAP7gN+/t1M59DXr25lmKjQx3Og4AAPBBnCk7Rf/IKNYzq/N042nDdPa4eKfjAAAAH8WZspPU0enRr97doac/ydOcEQP1w4smOB0JAAD4MErZSaje36a7XsrUmj1VumnecP344gkKDeakIwAAOHlebRLGmIXGmJ3GmBxjzD2HedwYY/7Y/finxpgZ3sx3LPUt7fpHRrEue+QTpRfU6Lf/NUU/vWwihQwAAJwyr50pM8YES3pU0nmSiiVtNMYstdZuO+hpF0oa0/1njqTHu/92TFNbh/61o1xvbS7VRzsr1Nbh0fBBffX3207TtJQYJ6MBAAA/4s3hy9mScqy1uZJkjHlF0iJJB5eyRZJesNZaSeuMMTHGmARr7V4v5vyMVzYU6edvb1N8/3BdN2eoLp2aqOkpMSx7AQAAepQ3S1mSpKKDbhfr82fBDvecJEmfKWXGmFsl3SpJQ4f27kKti6YlakJClGaPGKjgIIoYAADoHd6cDHW4RmNP4jmy1j5prU2z1qbFxcX1SLgjGRQZrtNGDaKQAQCAXuXNUlYsKeWg28mSSk/iOQAAAH7Hm6Vso6QxxpgRxpgwSYslLT3kOUsl3dB9FeZcSXVOzicDAADwFq/NKbPWdhhj7pK0XFKwpGestdnGmNu7H39C0jJJF0nKkdQk6WZv5QMAAHCSVxePtdYuU1fxOvi+Jw762Eq605uZAAAA3IBVTwEAAFyAUgYAAOAClDIAAAAXoJQBAAC4AKUMAADABShlAAAALkApAwAAcAFKGQAAgAtQygAAAFyAUgYAAOAClDIAAAAXMF3bTfouY0yFpIJe/jaxkip7+Xvg+HE83IXj4S4cD3fheLiLG47HMGtt3OEe8PlS5g3GmHRrbZrTOdCF4+EuHA934Xi4C8fDXdx+PBi+BAAAcAFKGQAAgAtQyo7Pk04HwGdwPNyF4+EuHA934Xi4i6uPB3PKAAAAXIAzZQAAAC5AKQMAAHABStkxGGMWGmN2GmNyjDH3OJ0n0Blj8o0xW4wxm4wx6U7nCTTGmGeMMeXGmK0H3TfQGPNPY8zu7r8HOJkxkBzhePzUGFPS/RrZZIy5yMmMgcIYk2KM+cgYs90Yk22M+Wb3/bw+HHCU4+Hq1wdzyo7CGBMsaZek8yQVS9oo6Rpr7TZHgwUwY0y+pDRrrdOL/wUkY8yZkholvWCtndR9328kVVtrH+j+h8sAa+0PnMwZKI5wPH4qqdFa+zsnswUaY0yCpARrbaYxpr+kDEmXS7pJvD687ijH42q5+PXBmbKjmy0px1qba61tk/SKpEUOZwIcY639WFL1IXcvkvR898fPq+sXH7zgCMcDDrDW7rXWZnZ/3CBpu6Qk8fpwxFGOh6tRyo4uSVLRQbeL5QMH1c9ZSe8bYzKMMbc6HQaSpMHW2r1S1y9CSfEO54F0lzHm0+7hTYbLvMwYM1zSdEnrxevDcYccD8nFrw9K2dGZw9zHeK+z5ltrZ0i6UNKd3cM3AP7jcUmjJE2TtFfSg46mCTDGmEhJ/5D0LWttvdN5At1hjoerXx+UsqMrlpRy0O1kSaUOZYEka21p99/lkt5Q1xAznLWve/7GgXkc5Q7nCWjW2n3W2k5rrUfSn8VrxGuMMaHqKgB/tda+3n03rw+HHO54uP31QSk7uo2SxhhjRhhjwiQtlrTU4UwByxjTr3vCpowx/SSdL2nr0T8LXrBU0o3dH98o6U0HswS8AwWg2xXiNeIVxhgj6WlJ2621vz/oIV4fDjjS8XD764OrL4+h+3LZhyQFS3rGWvtLZxMFLmPMSHWdHZOkEEkvcTy8yxjzsqQFkmIl7ZP0E0lLJP1d0lBJhZK+aK1l8rkXHOF4LFDX0IyVlC/ptgNzmtB7jDGnS1olaYskT/fdP1LXPCZeH152lONxjVz8+qCUAQAAuADDlwAAAC5AKQMAAHABShkAAIALUMoAAABcgFIGAADgApQyAAAAF6CUAQAAuAClDIDfM8YkG2O+dJTH/88YM9+bmQDgUJQyAIHgHEkzjvL4HEnrvJQFAA6LUgbAr3Vvt/J7Sf9ljNlkjBlxyOMTJO2y1nYedN9wY8wOY8xTxpitxpi/GmPONcasNsbsNsbM7n5eP2PMO8aYzd3PO+LZOAA4FkoZAL9mrf1E0kZJi6y106y1eYc85UJJ7x3mU0dLeljSFEnjJV0r6XRJ31PXHnqStFBSqbV2qrV20hG+DgAcF0oZgEAwTtLOIzx2gQ5fpvKstVustR5J2ZI+tF2bBW+RNLz7OVsknWuM+bUx5gxrbV0P5wYQQChlAPyaMWaQpDprbfthHusrKcZaW3qYT2096GPPQbc9kkIkyVq7S9JMdZWzXxlj/qcnswMILCFOBwCAXjZC0uFKlySdLemjk/3CxphESdXW2r8YYxol3XSyXwsAKGUA/N0OSbHGmK2SbrXWrjnosQslvXYKX3uypN8aYzyS2iV9/RS+FoAAZ7qmSABA4DHGZEqac7ihTQDwNkoZAACACzDRHwAAwAUoZQAAAC5AKQMAAHABShkAAIALUMoAAABcgFIGAADgApQyAAAAF/h/LaWGkc4SAT0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the results\n",
    "x=2/10000*1e3*np.arange(len(S))\n",
    "fig,ax=plt.subplots()\n",
    "ax.plot(x,1-S.real)\n",
    "ax.set_xlabel(r'$t$ / ms')\n",
    "ax.set_ylabel(r'1-$S/S_0$')\n",
    "fig.set_size_inches([10,7])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba534529",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}