{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "98ba5968",
   "metadata": {},
   "source": [
    "# Introduction to LLG equations and the Heun algorithm\n",
    "\n",
    "### © Branislav K. Nikolić and Jalil Varela-Manjarres, University of Delaware\n",
    "[PHYS660: Computational Methods of Physics](https://wiki.physics.udel.edu/phys660) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eab2eaf6",
   "metadata": {},
   "source": [
    "## What is covered in this notebook\n",
    "\n",
    "* Introduction to LLG equations:\n",
    "    * Main contribution to the energy \n",
    "    * Dynamics of LLG equations (damping and driving terms)\n",
    "    * Analytical solution for a single magnetic moment under a magnetic field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a9efbf96-5d6c-4110-8b0a-9b225bcb28cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Libraries and global parameters \n",
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rcParams\n",
    "### Parameters for plotting \n",
    "rcParams['axes.linewidth'] = 2\n",
    "rcParams['text.usetex'] = True  ### \n",
    "fs = 25                         ### Font size"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "509f7742-a95a-4d06-adae-64226a0fddf6",
   "metadata": {},
   "source": [
    "# LLG Theory"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45e41edf-a8bb-435a-adfa-5ce9e4d6c63f",
   "metadata": {},
   "source": [
    "The LLG equations describes the classical evolution of a set of localized magnetic moments, inside an effective field that accounts for the interactions between spins, anisotropies or in general coupling with localized field, for an arbitrary number of spins, the LLG equations reads:\n",
    "\n",
    "- $\\frac{\\partial \\mathbf{S}_i}{\\partial t}=-\\frac{\\gamma}{1+\\lambda^2}\\left[\\mathbf{S}_i \\times \\mathbf{H}_{\\mathrm{eff}}^i+\\lambda \\mathbf{S}_i \\times\\left(\\mathbf{S}_i \\times \\mathbf{H}_{\\mathrm{eff}}^i\\right)\\right]$\n",
    "\n",
    "Being $\\mathbf{S}_i$ the normalized and localized magnetic moment at site i, and $\\mathbf{H}_{\\mathrm{eff}}^i$ the effective field produced by the coupling  to external fields or interactions. The parameter $\\lambda$ is the Gilbert damping and  leads to the magnetization aligning with \n",
    "the effective field and therefore energy minimization. of the spin and $\\gamma$ is the gyromagnetic ratio and is related to the amount of angular momentum of the moments around the effective field. \n",
    "\n",
    "Microscopically, the effetive fields comes from the energy functional, and is calculated with its derivate\n",
    "- -$\\mathbf{H}_{\\mathrm{eff}}^i=-\\frac{1}{\\mu_s} \\frac{\\partial H}{\\partial \\mathbf{S}_i}$\n",
    "\n",
    "\n",
    "The energy functional then contain the different types of interactions or couplings between the localized magnetic moments \n",
    "-  $H=H_{\\mathrm{exc}}+H_{\\mathrm{sd}}+H_{\\mathrm{B}}+H_{\\mathrm{ani}}+H_{\\mathrm{dem}}$.\n",
    "\n",
    "where, \n",
    "\n",
    "- $H_{\\mathrm{exc}}=-\\sum_{i j} J_{\\mathrm{exc}} \\mathbf{S}_i \\cdot \\mathbf{S}_j$\n",
    "\n",
    "- $H_{\\mathrm{sd}}=-J_{\\mathrm{sd}} \\sum_i\\left\\langle\\hat{\\mathbf{s}}_i\\right\\rangle \\cdot \\mathbf{S}_i$,\n",
    "- $H_{\\mathrm{B}}=-\\sum_i \\mu_{\\mathrm{S}} \\mathbf{S}_i \\cdot \\mathbf{H}_{\\mathrm{ext}}^i$.\n",
    "- $H_{\\mathrm{ani}}=-J_{\\mathrm{ani}} \\sum_i\\left(\\mathbf{S}_i \\cdot \\mathbf{e}_{\\mathrm{ani}}\\right)^2$\n",
    "- $H_{\\mathrm{dem}}=J_{\\mathrm{dem}} \\sum_i\\left(\\mathbf{S}_i \\cdot \\mathbf{e}_{\\mathrm{dem}}\\right)^2$.\n",
    ".\n",
    "\n",
    "\n",
    "Taking the derivate of this elements, the effective field reads: \n",
    "\n",
    "- $\\mathbf{H}_{\\mathrm{eff}}^i=\\frac{1}{\\mu_B}\\left[J_{\\mathrm{exc}} \\sum_j \\mathbf{S}_j+J_{\\mathrm{sd}}\\left\\langle\\hat{\\mathbf{s}}_i\\right\\rangle+\\mu_B \\mathbf{H}_{\\mathrm{ext}}^i+2 J_{\\mathrm{ani}}\\left(\\mathbf{S}_i \\cdot \\boldsymbol{e}_{\\mathrm{ani}}\\right) \\boldsymbol{e}_{\\mathrm{ani}}-2 J_{\\mathrm{dem}}\\left(\\mathbf{S}_i \\cdot \\boldsymbol{e}_{\\mathrm{dem}}\\right) \\boldsymbol{e}_{\\mathrm{dem}}\\right]$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ed286ab-6be3-447d-9029-8a3ef7253c91",
   "metadata": {},
   "source": [
    "## Heun's scheme "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9726158-f54f-46e7-9041-9961719c2a90",
   "metadata": {},
   "source": [
    "The Heus method is a two step method of second order to solve differential equations. it is commonly used to solve long dynmics with larger time stepd due to its predictor corrector nature.  the\n",
    "Heun scheme, despite its relative simplicity, is sufficientl \n",
    "computationally efficient that it is still the most widely us d\n",
    "integration scheme for stochastic magnetization dynam. In the first step, we calculate the effecrive field on each site, and then we compute the change of local spins: \n",
    "\n",
    "- $\\Delta \\mathbf{S}_i=-\\frac{\\gamma}{\\left(1+\\lambda^2\\right)}\\left[\\mathbf{S}_i \\times \\mathbf{H}_{\\mathrm{eff}}^i+\\lambda \\mathbf{S}_i \\times\\left(\\mathbf{S}_i \\times \\mathbf{H}_{\\mathrm{eff}}^i\\right)\\right] \\Delta t$\n",
    "\n",
    "which is applied to the initial spins \n",
    "\n",
    "- $\\mathbf{S}_i^{\\prime}=\\mathbf{S}_i+\\Delta \\mathbf{S}_i$\n",
    "\n",
    "And the new spins $\\mathbf{S}_i'$ are normalized. After this step, the effective field $\\mathbf{H}'^i_{eff}$  is remcomputed and using this new field and $\\mathbf{S}_i'$, the $\\Delta S'$ is computed again. Then, the next step is calculated through: \n",
    "\n",
    "- $\\mathbf{S}_i^{t+\\Delta t}=\\mathbf{S}_i+\\frac{1}{2}\\left[\\Delta \\mathbf{S}_i+\\Delta \\mathbf{S}_i^{\\prime}\\right]$\n",
    "\n",
    "and finally the resultant magnetic moment is normalized.ail."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b41a6f8f-3b41-4951-8acd-4a3970aaa2a8",
   "metadata": {},
   "source": [
    "## Analitycal solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d46af7c3-51c5-454f-9f45-facdf1972d7e",
   "metadata": {},
   "source": [
    "In order to test the implementation of the LLG equations we review the analitycal solution for a single spin interaction with a magnetic field in the z direction, with the magnetic momemnt initally at the x direction. The analytical solution of this is: \n",
    "\n",
    "$\\begin{aligned} & S_x(t)=\\operatorname{sech}\\left(\\frac{\\lambda \\gamma H}{1+\\lambda^2} t\\right) \\cos \\left(\\frac{\\gamma H}{1+\\lambda^2} t\\right) \\\\ & S_y(t)=\\operatorname{sech}\\left(\\frac{\\lambda \\gamma H}{1+\\lambda^2} t\\right) \\sin \\left(\\frac{\\gamma H}{1+\\lambda^2} t\\right) \\\\ & S_z(t)=\\tanh \\left(\\frac{\\lambda \\gamma H}{1+\\lambda^2} t\\right) .\\end{aligned}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "03a00f0e-1826-47e5-8cae-7ee5dcd148b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "GAMMA_R = 1.#1e-2#1.760859644e-4           ### Gyromagnetci ratio\n",
    "MBOHR = 5.788381e-5                ### Bohrs magneton\n",
    "g_lambda = 0.2                     ### gilbert damping factor\n",
    "n_sites = 1\n",
    "t0,dt,tf =0.0, 0.01, 50\n",
    "h0 = np.array([0., 0., 1.])        ### external magnetic field vector\n",
    "fs = 25\n",
    "\n",
    "times=np.arange(t0,tf+dt,dt)#[0:-1]\n",
    "H = np.linalg.norm(h0)\n",
    "\n",
    "omega_0 = GAMMA_R*H/(1+g_lambda**2)\n",
    "damp = omega_0*g_lambda\n",
    "mx_a = 1/np.cosh(damp*times)*np.cos(omega_0*times)\n",
    "my_a = 1/np.cosh(damp*times)*np.sin(omega_0*times)\n",
    "mz_a = np.tanh(damp*times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "347b4e04-c31f-4747-8b9b-53cfafc5b708",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fe52bfb84c0>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHSCAYAAAB8Y3pCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABmIUlEQVR4nO3deXhb1Z0//rfi7CGxbEMgEEosh7W0TWQHSr+lQCxBWgptsWTTTqdD21iiyzztMI0V9+nGdFpHJp12thLJoZ2WbrZF+LWdLkRygJZShlhKKFC2SDY0CxAiXycs2Wz9/jg5shxb+726V9b79Tx6rqzl3uNEid8+y+eYEolEAkRERERUMWbp3QAiIiIiKi0GQCIiIqIKwwBIREREVGEYAImIiIgqDAMgERERUYVhACQiIiKqMAyARERERBWGAZCIiIiowszWuwFUOosWLcLRo0dRVVWFpUuX6t0cIiIiysGrr76KsbExzJ8/H2+88YYq5zRxJ5DKUVVVhfHxcb2bQURERAWYNWsWxsbGVDkXewAriAyAs2bNwrJly/RuDhEREeXgwIEDGB8fR1VVlWrnZACsIEuXLsW+ffuwbNky7N27V+/mEBERUQ6WL1+Offv2qTp9i4tAiIiIiCoMAyARERFRhWEAJCIiIqownANYIEVR0NjYiHA4DLPZrOq5/X4/gsEgamtrAQDxeBxutxs2m03V6xAREVFlYgDMg6IoGBwcRDAYRHd3NwARztQKgIqioLm5GU1NTejv75/2cZ/Pp8q1iIiIqHJxCDhHNTU1aGxsRH9/P9ra2mC1WlW/Rnt7O2pra6eEPLPZjIGBAfT19cHv96t+XSIiIqos7AHM0cjIiKbnD4VCCAQCCAaD0z5vNpvhcrngdrvR2tqq+rAzERERVQ72ABqE1+uF2WzOOM+vra0NANgLSEREREVhADQARVEQCoVgsVgyvk4+39vbW4pmERER0QzFAGgAg4ODAJA1AMph30gkonWTiIiIaAZjADQAGehk2ZdMZEhkCCQiIqJCMQAawKFDhwAgr4Ud8Xhco9YQERHRTMdVwAagKErOr5UhMRaLZX3tsWPHcOzYseTXiUQi36bl5vnnAbcb2LAB+MAHtLkGVaSxMeDEiey3kyenf3xsDBgfn7id/nUxjyUSEzdg8jGX+8U+n+trS43XndnXpexsNuDv/k7vVmTHAGgAWvXmdXV14c4779Tk3JP84AfAQw8BDz8sju97n/bXJMMYHwdGR8XtyBHg8OHJt9THjhwB3noLOHpUHE+/nf742Jje3x0RUX4WL2YApBzJHsC6ujpVz9vZ2Yk77rgj+fWll16K/fv3q3oNAMDnPgc88wzwq18B3/wmkKaWIZWPEyeAAweAffuAvXuB/fuBV18FXntt6u3QodIGtdmzgTlzcrvNng3MmjX5VlVV+GPycZNJ3AfEfXlMd1/tx/J9j5Z4jcq7BmWmwT4RmmAANBA5F1At8+bNw7x585Jfm7T6n+H884HvfU8EwIEB4JVXgLPP1uZapIqxMRHs9uyZuEWjwEsvidD3yiv5DzHNnw8sWTJxW7x46tdnnAEsXAgsWCBu8+dPf19+PX/+1EDHH3BERMVjADSAQnb1yGXFcEnV1wONjUA4DPzud8Btt+ndIoIYnn3xReDJJ4G//EUcn3xShL3jxzO/d84c4NxzgeXLgfPOE5n+zDOnv9XVASm/axARkcExABpAPmFODhcbcis4m00EwEceYQDUyauvAo89Jm5//rP46zhyZPrXzp0LWCzAypXi1tAAXHCBCHvLl4tgJ4c5iYhoZmEANICGhgYAua3slQtGDNcDCAD/7/+J45/+pG87KsjBg2LUPRgU62+m+wjNnQtceinwzneK2zveAVxyiQh5VVUlbzIRERkAA6AByOLOuZSDka+xGnGW6VVXieOzz4oln0uW6NueGWh8HNi5E7j/fuCBB4Ddu6e+5rLLxF/Fu98NXHmlCHtz5pS8qUREZGAMgAZgs9kAZO8BlOEv25ZxujnzTDFpbP9+4OmnJwIhFWV8HPjDH4D77hPBb9++yc+/852A3Q40N4s/ciPODiAiImNhADQAs9kMq9WadXs3uWew2+0uRbMK8453iAD45JMMgEWKRoEf/UjcXnpp4vEzzgA++EHgpptE6OOCayIiyheneJdQKBRK+1xnZycAIBAIpH1Nf38/AMDlcqnbMDVdfrk4PvWUvu0oU2NjwLZtwLXXioUZ3/ymCH/V1cAnPwn87/+KeX8//znwsY8x/BERUWEYAEvE6XTCbrejpqZm2rl+DocDVqsVXV1d075fURT4/X74fD5jrgCWZAB8+ml921FmDh8GvvtdEfpaWsSmKiYTcMMNIuy9/LLYcOXGG0VtPCIiomIwABYgFoslh2sz9dilkr1/iqKknes3MDAAYOoQr6IoaG5uRkdHh7F7/wDgwgvFMRrVtx1l4sgR4FvfEuVX7rgDGB4GamuBL39Z1O/7/e+BW29l6CMiInWZEgluKZ0Lu92enIM3XQ+e7JULh8PTLtLw+/3weDxoampCMMtWad3d3di5c2ey1Es8Hofb7U4uFinU8uXLsW/fPpx33nnYu3dvUedK68ABsRBk1iyxmevcudpcp8y99RbwH/8B3HWX2EoNAC6+WITAj39c7JZBREQEaPPzmwGwgpQkACYSwKJFIuE8//xEjyABEH88gQDwpS9NLOy4+GLg618HWltZl4+IiKbS4uc3h4BJXSaT2F4CmL4qcQV78knguutE0HvpJbGF8o9/LKZLfvSjDH9ERFQ6DICkvlM7m3AeoHD8OHDnnWKr5IcfBhYsAL7xDVEv++//nsGPiIhKj3UASX2yB5ABELt3i22Rn3hCfP2hD4m5f297m56tIiKiSsceQFJffb04vviivu3QUSIB/Od/AldcIcJfXZ0o53L//Qx/RESkP/YAkvqWLxfH0/csqxCKAnz606KgMyB6/Xw+Fm0mIiLjYAAk9ckAqNVKYwN76ing5puBoSFgzhxg82bgH/9RrI0hIiIyCgZAUp8MgAcOiL3NKmSVw29/K4o2HzkiRsH7+oCmJr1bRURENBXnAJL6zj5bhL6xMeCVV/RuTUn8+78DN90kwt811wA7dzL8ERGRcTEAkvqqqsRuIMCMHwZOJICNG4EvfhEYHxdz/7ZvF4s+iIiIjIoBkLRx3nniOIMD4Pg48LnPAV6v+HrTJqCnh7vfERGR8XEOIGljhi8EGRsDPvlJ4N57xQKPLVsAl0vvVhEREeWGAZC0MYNLwcih3nvvBWbPFtu5ffSjereKiIgodwyApA05BDzDAmAiIcq6/OhHYqpjXx/wkY/o3SoiIqL8cA4gaUNWPZ5hq4A3bgS+/30x7PujHzH8ERFReWIAJG3MwAD4X/8FdHeL+z4f8Hd/p297iIiICsUASNqYYQHwN78BvvAFcf/b3wba2/VtDxERUTEYAEkbMgC+9ppYMlvGdu8G2tomFn9s3Kh3i4iIiIrDAEjaOPNMMVFufFyEwDJ18KDY4eONN4DmZuDuu7mvLxERlT8GQNLG7NkT22G8+qq+bSnQ2Jgo77J3L3DxxUAgAMyZo3eriIiIiscASNop83mAX/saMDAALFoEbNsGmM16t4iIiEgdDICknTIOgL/+tVjsAQD33ANcdpm+7SEiIlITAyBpp0wD4L59wD/8g7j/hS+IBSBEREQzCQMgaacMA+D4OHDbbcDICNDUNFH3j4iIaCZhACTtlGEA/I//AEIhYMEC4Cc/AebO1btFRERE6mMAJO0sXSqOZRIAn3pqosbfd74jVv4SERHNRAyApJ0y6gE8eVLM+zt2DPjAB4Dbb9e7RURERNphACTtnHmmOB46pG87cvC97wGRiCj1cs89LPZMREQzGwMgaUcWgjZ4AIzFRM0/ANi8GTjnHH3bQ0REpDUGQNKODIBvvAEcPapvW9JIJAC3G3jrLeC664BPfUrvFhEREWmPAZC0U10NVFWJ+wbtBfzJT8Sq33nzAJ+PQ79ERFQZGABJO7NmAbW14r4BA+Dhw8CGDeL+178OXHihvu0hIiIqFQZA0paB5wH+67+KBcoXXgj88z/r3RoiIqLSYQAkbcmVwK+9pm87TvPCC2LlLwD827+x4DMREVUWBkDSlkF7AO+4AzhxAli3DrjxRr1bQ0REVFoMgKQtA9YCfOAB4H//F5g9G/jud7nwg4iIKg8DIGlL9gAaZAh4fBzweMT9z38euOQSfdtDRESkBwZA0pbBhoB7e4EnngCWLAG+8hW9W0NERKQPBkDSloGGgE+cAL76VXF/w4aJbEpERFRpGABJWwYaAr7nHiAaBZYuBb74Rb1bQ0REpB8GQNKWQYaA33wT+Jd/Efe/+lXgjDN0bQ4REZGuGABJWwYZAr77buDAAWDFCsDl0rUpREREumMAJG3JHkBFAU6e1KUJb70F3HWXuP+1r7HoMxEREQMgaaumZuJ+PK5LE+65R2z5dsEFwMc/rksTiIiIDIUBkLQ1ezZgNov7OgwDHzsGeL3i/saNwJw5JW8CERGR4TAAkvZkL6CilPzSP/4xsHcvcO65wG23lfzyREREhsQASNqTAXBkpKSXPXEC6OoS9zs6gPnzS3p5IiIiw2IAJO3pFAB7e4GhIVH3r729pJcmIiIyNAZA0p4OATCRAP7t38T9L3wBWLiwZJcmIiIyPAZA0p4OAfDhh4Fdu0Twu/32kl2WiIioLDAAkvbkKuASBkDZ+3fbbUBtbckuS0REVBYYAEl7JV4F/NxzwK9/DZhMYviXiIiIJmMAJO2VeAj4e98Tx5tuAi66qCSXJCIiKisMgKS9EgbA114DfvQjcf+OOzS/HBERUVliACTtlTAA3nOP2PvXagXe9z7NL0dERFSWGABJeyUKgGNjgM8n7n/+82IOIBEREU3FAEjaK1EAfOABUfi5pgZoa9P0UkRERGWNAZC0J8vAvP662J9NI3ffLY633cbCz0RERJkYIgB2dnbq3QTSkgyAADA6qsklhoeB3/xG3GfhZyIioswMEQD9fr/eTSAtzZ4NLF4s7ms0DOz3i+3fmptZ+oWIiCib2VqdeHh4GLFYLOvrIpEIlBIVCCYd1dQAR45oEgCPHxerfwHgs59V/fREREQzjuoBcOvWrfB4PAx1NFlNDfDSS5oEwPvvB159FTj3XODmm1U/PRER0YyjagAcGBiAy+WCxWKB0+mEOXXuVxrhcBg7duxQsxlkRBquBP7BD8TxU58So81ERESUmao/Ljdu3Ij+/n60tLTk9b66ujo1m0FGpFEAfOklIBgU9z/5SVVPTURENGOpugikpqYm7/AHAPX19Wo2g4xI9garPDXgxz8Wiz+uvRawWFQ9NRER0YylagC0Wq0FvW9wcFDNZpARadADOD4O/PCH4j57/4iIiHKnagDkwg9KS4MA+Mc/ArGYqDBTQMczERFRxVI1ADqdTmzbti3v97EQdAXQIADK3r+2NmDRItVOS0RENOOpGgCbm5sRjUbzDoHlVgja7/fD6XTC7XbD7XbD6XQiFAqpcm673Z7xzyMQCMBut6tyrZKqrhZHlXYCOXIE6O8X9zn8S0RElB9VVwEPDw/D6XSiv78fa9asgc1mQ0NDAywZZueXUyFoRVHQ3NyMpqYm9Mv0cdrjPp+vqGvE43G43W54vV54PB40NTXBbDYjEomgq6sLkUhk0rXLhsoBsK8PePNN4OKLgauuUuWUREREFUPVAGi1WjF66gd8IpFAOByGyWRS8xK6am9vR21t7ZSQZzabMTAwgPr6ejQ2NsLlchV9rVgsBrfbPeVxn88Hh8NR9PlLTuUA+JOfiONttwEz6CNGRERUEqoGwNraWgBAa2trTkWggfIpBB0KhRAIBBCURedOYzab4XK54Ha78/r+p+NyuTA4OIhYLAZFUWC1WtHU1ASPx5OxN9XQVAyAe/cCDz8s7n/sY0WfLjeJBPCnPwGzZokuR6ZOIiIqY6oGQIvFgo0bN2L9+vV5va8cCkF7vV6YzWbYbLa0r2lra0N3dzf8fj86OjoKus50PYwzgooB8Be/EHns6quBt72t6NNll0gAH/848LOfia8//nHgRz8SYZCIiKgMqfoTzGKxFNRDZfRC0IqiIBQKZf3e5PO9vb2laFZ5kQHw6FHg+PGiTiVzWMl6/+65Z+KigBh//p//KdHFiYiI1KdqANyyZQvWrl2b9/uMXghati9bAJTDvpFIROsmlZ8lSybuF9EL+MwzwK5dYs/fkkyFPHEC+Na3xP277gI2bxb3v/Y14OTJEjSAiIhIfRzDyoEMdHKOYyYyJDIEnqaqCjjjDHG/iAD485+L47p1wJlnqtCubAYGgOFhcbHPfhb4/OeBs84C9u0DfvnLEjSAiIhIfarOATzd8PAwfD4fIpEI4vE4amtrYbFY4Ha7sWrVKi0vrapDhw4BQF4LO+LxeNHXjcViCIVCiEajaGhoKHpxie6qq4HXXy84ACYSOgz/9vWJY2srsHChuP/pTwObNgE//Sm3ICEiorKkWQ9gZ2cnGhoa4PV6EQwGEQ6HEQwG4fP50NjYiFtvvRWHDx/W6vKqyqdOoQxosVis4OvFYjF4PB7EYjG0trbC6/WitrYWjY2NCAQCBZ9Xd0UuBNm5E4hGRQ67+WYV25XO2Bjw//1/4r7TOfF4a6s4/v73ohghERFRmVG9B3B0dBSNjY2Ix+NoaWnBmjVrkqFIURQcOnQI4XAYfX19CIVCGBoawuLFi9VuhqrU6M3L51qBQABer3fS4w6HA1arFQ0NDejo6Jjy/HSOHTuGY8eOJb9OJBKqtzcvRQbAn/5UHD/84RJt/fbEE2LrusWLgfe+d+LxVauAFSvE0HAwCHzoQyVoDBERkXpUD4BOpxM2mw1btmzJ+LrR0VF0dHRg7dq12Llzp9rNUJXsASxFuZqenh5YrdZpn7NYLOjo6EB3dzfcbnfWRSldXV248847tWhmYYoIgGNjE6OxJRv+lcUGr75arDqRTCbgAx8Avv99YMcOBkAiIio7qg4B9/T0wGKxZA1/AFBdXZ0cDt66dauazdCMnAuopXThT5L7AE+3S8jpOjs7MTo6mryde+65qrSxYEUEwD/9CXj5ZcBsBkq2FfJDD4njtddOfU4+9uCDJWoMERGRelQNgKFQKKfwl2rLli2G39u2kIUXuawYLkRTUxMA8Wedzbx587BkyZLkTfdt+YoIgPfdJ44f+hAwd66KbUonkQAefVTcv/rqqc9fc404Pvkk8NprJWgQERGRelQNgDU1NQW9z+iFoPMJc3K4WKvVuqnnLWahiS4KDIDj4xMBsGTbIO/bJ4JdVZWY83e6pUuByy4T92VQJCIiKhOqBsBCe5h075nKoqGhAUBugUsuGMm3B1BRFHR3d+fUsydVSgD8v/8TeWzx4hIO/+7aJY6XXgrMnz/9a668UhwNPoeViIjodKoGwEJXyxo9yMjFFrmUg5GvyTaX73R+vx8ejwfO1HIjWWg1zKyZAgOgrHxz883AvHkqtykdGQBXr07/mjVrxJEBkIiIyoyqAbC+vh73339/Xu/ZvHkzVmf6IWsANpsNQPagKsNfIfshn36tbNco9jq6KCAAJhITAbBkw79A/gFQ7xI7REREeVC1DMymTZuwcuVKAMBHPvKRrK/v6elBV1dXSVbXFsNsNsNqtWbd3k3uGZzLCt3TWa1W+Hw+uFyujK+TQ8Q2m638dgUpIAAODgIvvSTq/t1wg0btms7u3eKYKQC+851iRUo8DgwNAeUWyImIqGKpvhNIX18fWlpasG7dOmzduhW7d+/G8PAwhoeHsXv3buzYsQOdnZ248MILcfvtt+c1501PnZ2dAJBxJw65mjlTiEv3/dpstpxWQ/t8vknHslJAAJR/3DfeCCxYoEGbpvPGG6LIMwC84x3pXzd3LvCud4n7p8I/ERFROVA9AFqtVgwODuLgwYNwuVxobGxEQ0MDGhoa0NjYCLvdDq/Xi0Qige3btxt++FeSO3F0dXVN+7yiKPD7/fD5fGl75pxOJ+x2O2pqaqadT5htDmAgEEAoFILX6y2/4V8g7wCo2/DvCy+IY12duGUiA+CTT2rbJiIiIhVpshew1WpFOBzGli1bsGrVKiQSieStvr4eXq8Xe/bsQXNzsxaX18zAwACAqUO8iqKgubkZHR0dOfX+KYoy7XxCm80Gu90Ou90+Zbi5u7sb7e3t8Pl86OjoKPZb0UeeAXD3biAWEz1/73+/ds2a4tlnxfHii7O/VvYQMgASEVEZUX0ruFQulyvrnLZyYjabEQ6H0d3dDafTmVyFG4/H4fV6sy7g8Hq98Hg8aGpqSrtK2OVywWazwePxTAqBNpsN4XC4PHv+JBkA33oLOHECmDMn48u3bRPHdeuAM87QuG2pnntOHHMJgJdfLo5PPaVde4iIiFSmaQCcqQrtgcs1EFssFsPvjlKQJUsm7o+OAmeemfHlv/qVOH74w9o1aVr5BEDZAxiLibmDixZp1y4iIiKVaDIEnK9tsquHZrbZsycCUpZh4OFh4C9/AWbNEgtASiqfAHjWWWJXkEQCePppbdtFRESkEkMEwPb2dr2bQKWS4zxA2fv33vdmX4ehqkQCeP55cT+XAAhM9AJyGJiIiMpE3kPAu3fvxqrp9kY99Vw+4vE4YrFYTjts0AxRXQ3s359zAPzQh0rQplQHDgCvvy66Hk9tAZjV5ZcDAwNcCEJERGUjrwC4ceNG3HXXXWhsbMTjjz8+5fm1a9diNM9tvqjC5NADqCjAww+L+zffrH2TJhkaEse3vU3U+cvFZZeJoxw6JiIiMri8AmAsFkMikcDIyMi0z8tVsXKhQ12WsbtDhw4hGo1yDmAlySEA/u53wMmTIled2limdGQB6BUrcn/PRReJIwMgERGVibwCYE9PD+x2e9pyJxaLBRs3bsT69evzakS2oEgzSA4B8Je/FMeS9/4BhQVAOVdweBg4dgyYN0/lRhEREakrr0Ug1dXVaG9vR319/bTP2+12NDU15d2IdOejGShLADx+XPQAAjrM/wMKC4DnnCMKFY6PA9GoFq0iIiJSlaqrgDds2JB2gUgmg9xHtXJkCYB/+ANw+DBw9tnAFVeUsF1SIQHQZJroBZQriImIiAzMEGVgqIJkCYBy+Pemm8RC3JIrJAACnAdIRERlpaQ7gQwMDCAYDMJkMqGhoQFNTU0F9RhSGcsQABOJifIvusz/Gx8HXnxR3M83ALIHkIiIyoiqAbCurg6HDh1K+3xzczOam5sxOjqKeDyOUCgEn88Hr9eLJanbhNHMJQPg4cNTnvrLX4CXXgIWLACybKusjQMHxB7FVVXAeefl914GQCIiKiOqBsBEIpHT66qrq5MLSkZHR+HxeHD33Xer2RQyKhn0jxyZ8tRvfyuOzc0iBJacHP49/3yxbV0+OARMRERlRNVZViaTKe/3VMseIaoMMgBO0wMoA2DJ9/6VCp3/B0wEwIMHgTR1MomIiIxC90Ugo6OjXAVcSRYvFsfTAuDICPDnP4v7739/idskyfl/F1yQ/3vPOANYtkzcZykYIiIyuLyHgIdlL8lp5PDviy++mHUoWFEUxONxRCIRdHV1obW1Nd9mULlK0wMYDAJjY2L3j0Lylyr27RPH5csLe7/FIuYRxmJAAfUwiYiISiXvABgOh7Fz507EYjGEQqEpe/9aLJacz5VIJGCz2Tj/r5KkBsBEQtTQw8Tw7wc+oFO7gIkAmO8CEMliAf70JxEAiYiIDCzvANjS0oKWlpbk15FIBB6PBwMDAzCZTJOey8RiscBut6O5uTnfJlA5kwFwfBx4801g0SKMjwO//714WLfhX0CdAAgwABIRkeEVvQrYarUiGAzC7XYjEAigr69PjXbRTLVwoajwPD4uegEXLcKuXcArr4hpdO99r45tYwAkIqIKodoiEJ/Pxz19KTuTaUopGDn8a7cDc+fq1K6TJ0UKBRgAiYhoxlN1FXBnZ6eap6OZ6rSFIL/7nfhS1/l/L78seiWrqoClSws7hwyAL70kCkoTEREZlKqFoHOd/3f77bcnawba7XbccsstajaDjC6lFMxrrwGPPSa+XLdOvyYlh3+XLSt8E+Jly4D584GjR0UIbGhQr31EREQqKulewNKWLVuS9wcGBrB161asX79ej6aQHlJ6ALdvF4uB3/nOwquvqKLY+X+AGN62WIC//lUMAzMAEhGRQeleCLq5uRk+n0/vZlAppQRAQ5R/AdQJgADnARIRUVnQrAdw9+7diGX4IRiPx6EoCnp7e1FbW6tVM8iITgXAMeUIHnhAPKRr+ReAAZCIiCqK6gFweHgYdrs9Y/hLVV9fj/7+frWbQUZ2KgAOPrMIr70GVFcDV12lc5v27xdHBkAiIqoAqgdAm82GeDyOlpYWrFmzBmazGf39/bDb7TCbzQDEVnDBYBBms5l1AyvRqQD4wFMibNntwJw5ejYI7AEkIqKKomoA7OnpgcViQTgcRnV19aTn2trasETO/QKwYcMGDA0NobOzE11dXWo2g4zu1Crg7XtEWLr+ej0bcwoDIBERVRBVF4EEAgEEAoEp4c9isWBwcHDK6+vr6+FyubB161Y1m0FGt2QJRrEEj72yAoDoAdSdWgFQFkNXFGBkpLhzERERaUTVAFhfXz+pl0+qra1FMBhM+55oNKpmM8jolizBQ7gWY4kqXHQRsGKFzu15/XVxA0Qtv2IsXAicdZa4/+KLxZ2LiIhII6oGwJqammkfX716NUKhkJqXonK2ZAm2Q4z7GqL3T24Bt3Ch2JC4WBdcII4MgEREZFCqBsBEIpH2ucbGRtx///3TPpfrimGaIVICoCHm/8kAePbZophzsWSX5vBw8eciIiLSgKoB0G63Y9u2bQCArVu3oq2tDYdP7ffqcDjgcDjw4IMPTnrPwMAAA2CFGXr9LOzBhZiNE7j2Wr1bg8kBUA3sASQiIoNTNQA2NzcjGo3irrvugsvlQiAQSJZ5sdlsWLt2LWw2G9atW4fOzk60tbXh+uuvh9vtVrMZZHDBp84BALx7ThjTTBktvZdfFke1AiB7AImIyOBUrwO4YcMG7Nq1Cxs2bMDo6ChaW1uTzwUCAaxduxbbt29HMBhEIpGA1WrlPsAVZvvjYq7o9YntAN6tb2MA9gASEVHF0WQruNWrV2P16tVTHq+urkY4HEYgEEAoFEJDQwNcLpcWTSCDOnkSGPjzAgDA9Sd/A5zo1L8KNAMgERFVGM32As5EzgekyjM4CCijs2DGCJowCBw5Aui9F7RWAfDQIVFeRo2VxURERCpSdQ4gUTayHGTzrIdQhXHg1CIhXakdAKurgVPbHrIXkIiIjEjVHsBs27oNDQ3B4/HAZDKhtrYWZrMZDQ0NqK2txS233KJmU8igtm8Xx+sX/Qk4AmMFwHPOUe+cK1YAu3eLAPj2t6t3XiIiIhWo2gPo9Xpx0UUX4Yknnpj2+fr6evT19aG3txculwsjIyNwuVyTForQzHX4MPDnP4v79rrIxIN6U7sHEJgYBuZKYCIiMiBVA6DFYsGePXtgtVrx5S9/OeNrV69ejS1btmDTpk0ZC0jTzPHgg8DYGLByJVBfdyr4HTmib6PeeGNiGzg1A6AsBcMhYCIiMiBVA6DD4UA0GsWKFSvg9XpxxRVX4MUsPwA7OjpQXV2tZjPIoOT8v+uvB7B4sfhC7x5A2fs3f/5Em9TAHkAiIjIwVQOgyWRCfX09otEovvSlL2FwcBAWiwX33HNPxvc1NTWp2QwyqOT8v+uBZAVoowRAtbaBk1gKhoiIDEyzVcBerxfbt2/HkiVL4HK5sG7duuS2cKczyxWTNGMNDQEvvABUVUFs/2bEAKgmDgETEZGBqRoAFUWZ9LXNZsPw8DBuueUWbN++HRaLZcpewABQq3cdONLcwIA4XnmlqJIy4wOg7AF8+WXg6FF1z01ERFQkzesAVldXo7+/H1u2bEE8HofNZsNnP/vZSa8xqTn0RoYkA6DNduoBowVANUvAAKK4tSwA/dJL6p6biIioSKoGwFgslvY5l8uFPXv2YNWqVdiyZcukcjFcBTyzJRLAjh3i/tq1px40WgBUuwfQZOJCECIiMixVA2AoFMq46tdisSAcDmPDhg3JcjGbN29mD+AM99e/Aq++CixYALz73acelAFQ7zIwL78sjmoHQIALQYiIyLBUDYCJRAI2m23aeX6pvF4vBgcHsWLFCnR0dKCnp0fNZpDByOHf974XmDfv1INGKwOjRQAsx4Ugo6PA3XcDn/sc8C//Ajz3nN4tIiIiDagaANvb27F69Wp0dXXhhhtuQFtbG3bIsb/TWK1WRKNRtLe3Y3x8XM1mkMFMGf4FjDMEfPCgOJ51lvrnLrch4IcfBi65BPjsZ4Hvfx/4+teBSy8FOjtFBW8iIpoxVN0L2OfzFfSejRs3qtkMMpCTJ4GHHhL3m5tTnjBKADx0SBzPPFP9c5dTD+CjjwLr1okVyxdeCNxyC/Dkk8Bvfwts2gTs3w/88IfALM3XjRERUQmoGgALVV9fr3cTSCO7dolRxepqwGpNecIIAXBsDBgZEffr6tQ/f7n0AMbjIvAdPQp84ANAfz+wcKF47mc/Az7xCeDHPwYaGoCvfU3fthIRkSoM8ev8tm3b9G4CaUQO/157rSgCnWSEAKgogJx+oEUAlD2A+/cDx4+rf361fPGLYi7kpZcCfX0T4Q8APvYxYOtWcf8b3xDDxEREVPYMEQDb29v1bgJpRC4AmTT/D5i8ClivMkBy+HfxYmDuXPXPv3Sp2GN4fBzYu1f986shHAbuvVeUrfnhD4FFi6a+5rbbgE99Svw9ud0sbE1ENAPkPQS8e/durFq1Ku1z+YjH44jFYlN2EKGZ4dgx4JFHxP20AXB8HHjzzemDh9ZkANSi9w+YqAX43HNiHqDFos11ivHVr4rjxz8utmlJ5zvfEfMBn3sOuOuuifcREVFZyisAbty4EXfddRcaGxvx+OOPT3l+7dq1GB0dVa1xVN4eewx46y3REfb2t5/25IIFYkHB+LgYBtYjAL72mjhqsQBEWrFChCYjzgPcvRv43e/E2Hy2uX1mM/Bv/yaGhO+6C7j9dm1WThMRUUnkFQBjsRgSiQRG5MT508g9fV0uFwCgLkvPyqFDhxCNRjkHcIZKLf8ypda3ySR6ARVFBMBly0rdPO17AAFjLwT57/8WR4cDWLky++vb2oDNm4FIBPj2t4Hvflfb9hERkWbyCoA9PT2w2+2wJTd0ncxisWDjxo1Yv359Xo3IFhSpPMkAOKn8S6rUAKiHUgRAuRDEaAFQUYCf/lTc/9zncnvPrFmiJMz114s6gV/6EnDeeZo1kYiItJPXIpDq6mq0t7enLdtit9vR1NSUdyNYBmbmef11MQQMTDP/T9J7JXCphoAB49UCvPdeMT7/jneILVpyZbMBV18tVjV/73uaNY+IiLSl6irgDRs2pF0gksng4KCazSADeOQRUQR6xYoMax/0DoCV3AP4s5+J4/r104zPZ2AyAR6PuO/ziZ5EIiIqO4YoA0MzT9ryL6lSS8HooZRzAPfuFYnYCIaGRPfsrFlAa2v+73//+8WqniNHRAgkIqKyU9KdQAYGBhAMBmEymdDQ0ICmpqaCegzJ+Kbd//d0evcAlmII+JxzRI3B48eBffsmAqGeenvF8brrRPvyNWsWsGGDqA/47/8O/NM/aVNHkYiINKNqD2C2xRzNzc3YtGkTNm7ciObmZuzcuROf+cxncFjv/WBJVfG42AIOyBIAFy8Wx5k8BDxrlvFWAv/iF+J4662Fn+OjHxXh8cAB4P771WkXERGVjKoBMJHjjg7V1dWor69He3s7Nm3aBI+cU0QzwkMPiU0jLr00S3UXvXsASxEAAWMFwJdeAp54QgTTj3yk8PPMnQvIHXy+/3112kZERCWjagA05TOZ/JTq6mo1m0AGkLX8i6RnAEwkSjMEDBhrIchvfiOO73lP8cHX7RZFpP/wB+DJJ4tvGxERlYzui0BGR0e5CniGyWkBCKBvADxyZGJRhtY9gEYqBSMD4Ac+UPy5zjsP+PCHxX32AhIRlZW8F4EMp+nFkMO/L774YtahYEVREI/HEYlE0NXVhdZCViKSIe3fDzz7rKgWcu21WV6sZwCUw7/z5wMLF2p7LaP0AL711kT37I03qnPOz38euO8+UVdw0yaAPfpERGUh7wAYDoexc+dOxGIxhEKhKXv/WvLY8D6RSMBms+Huu+/OtxlkUA8+KI5WK1BTk+XFepaBKdXwL2CcOYAPPihC4PLlogC0Gq65RpSEefpp4Mc/Bv7xH9U5LxERaSrvANjS0oKWlpbk15FIBB6PBwMDAzCZTJOey8RiscBut6M560QxKic5D/8CxugBLMU2hLIH8G9/A8bGxLw5Pfz2t+J44435FX/OxGQCPvtZsZ3c3XeLHkG1zk1ERJopug6g1WpFMBiE2+1GIBBAX1+fGu0yNL/fj2AwiNraWgBAPB6H2+1Ou0eyUa+htkRiIgDmlOv1LANTygC4bBkwZw5w4oQYIz//fO2vOR35l3PDDeqe9+MfF7uDPPMM8PDDOYz9G8T4uKjN+PrrwKJFwLnnArNLWhqViEg3qv1v5/P5EA6H1TqdISmKgubmZjQ1NaG/v3/ax31F7oxQimtoJRYTVUbmzMlxe1nZA3jaNIKSKOUQcFWVCH2xmBgG1iMA5jU5M09LlogQuGWLWAxi5AB49Kiog/iLXwB/+pMIf9Ls2cAVV4iA/A//YIyi3UREGlF1FXBnZ6eapzOc9vZ21NbWTglgZrMZAwMD6Ovrg9/vN/w1tCLXF7z73aJDJSs95wCWsgcQ0H8hyEMPieOqVTlMzizAZz4jjvffL8Km0SQSwD33ACtXAp/8JPDAAyL8zZ4t/jzmzBGrwh99FPj614H6euCWW4C//EXvlhMRaULVADjd/L/Dhw+nXTlcTkKhEAKBQNqi1WazGS6XC263G4qiGPYaWspp+7dUcsXo0aNiq7RSqrQAKFfnXHedNud/5ztFt+/Jk8DWrdpco1D79wPr1gHr14sh3+XLgW9+U4S7t94SW9ccPQpEo0BPj/gAJxIizL7rXcAnPgG8+qre3wURkao0qQM4PDyMz3zmM6iqqkJNTQ1Wrlw56fmhoSG0trZiq9F+UGTg9XphNpszzsFra2sDgIJ76EpxDa0kEnkUgJbkHECg9PMASzkEDOhfC1AGwJzTeQE++1lx9PnEfEcj2L1bDOtu3y5K/tx1F7BnD/CVr4iV0HLO36xZgMUiQuLAAPDXvwKn/q3h3nuBSy4RwTbH3Y6IiIxO9QC4Y8cOWCwW+Hw+1NfXo7m5GfX19ZNeU19fj76+PtTX12Pz5s1qN0F1iqIgFAplLXEjn+/t7TXkNbT09NOik2TBAuDKK3N8U1XVxFhxqQNgqXsA9SwF87e/id6tqirg6qu1u84ttwBnnSV63H79a+2uk6sdO0Sv5L59Yl/CXbuAL30JmDcv+3svvVTME/y//xPD5iMjYuu7G28EXn5Z86YTEWlN1QA4NDQEh8MBl8uFaDSKPXv2YPv27WlLw8hwuEN2HRmU3KkkWzgzm80ARGkcI15DS3KB6dVXi21ic6ZXKZhKGgKWvX+NjRN/3lqYN0/0oAGiJIye/vhH4KabgDfeEF3Sjz4qevHydcUVwM6doudw3jzgd78DLr9cDA+Xo5Mn2YtJRABUDoAbN26E1+vFli1bJvX6ZdojuKWlZdJqVyOSYUuWZMlEBrh8A1oprqGlvId/Jb0CoF5DwC+9JMqPlJLW8/9Sud1ipXEoBDz3nPbXm044LLa6e/NNMffvN78BTv3iVJDZs0XPYTgsegMPHRK9nZ/6lD4ljDJ57TWxwOXb3wY+9jHRA3rBBcAZZ4jvY84c0RO8ZInYyq+pCWhpAf7pn8QK7kceAQw4v5iI1Kdq0auRkRG0t7fn/b5sW8fp7dCp3iJzHj9E4vG44a6hlZMnJxaZ5j3FrFJ6AM89V/zgPXECOHBA/PAtlUceEcdrrtH+WhdcAHzwg2IIeMsW4Lvf1f6aqQ4eBD7yEbHC97rrgG3bchvyzcXb3w489phYJdzdDfzwh6Lr+wc/KOA3H5WcPCnK2fz+9+K2e3f29yQSYuX9kSNiuH668l3nny8W9qxaJRbCrFoFNDSIuZJGMT4uQv6RI+Lv+/XXxdfj4+J7lEeTSXwGFiwQt/nzJ+4vXKhfYXYinakaAPPZBi7VyMiIms1QXT4rbmWAi8Viul/j2LFjOHbsWPJrrYJ2JCLym9kMrF6d55vlSuBSBsA33xSrP4HSBcDZs8UP1eFhcStVADx4UCx6AER9nlL4zGdEAPzhD4F//dccawKp4ORJ4KMfFXMeL7pIDNMuWKDuNebNE3se33ijqBU4NATYbGIBjNcretq0lkiI+Yz33gv8/OfAK69Mfv6ii8Rw/6pVopzN+ecDS5eK4DNvnvgl5MgR8W9u3z7RK/3ii6KQ95NPiq//9jdx+81vJs67aNHUUHj55er+/SYSoi7ogQPi9vLLE8fU24EDYvW2Gv+nzZ8vvodFi8Tfn7yf6Tbd60ymycFTHhMJ8Wd+7JhYbX7sWG63kycn306cmPrYyZPiezCZRDiXt9SvU+9XVeV+zOe1cpQv9e+jku+/731ilMDgVO8BLES+YanUStHTpsU1urq6cOedd6p+3tPJ4d9rry3gl2k9ikHL3r/Zs7WdE3e6FSsmAuD/+3+lueZjj4njJZdoU/9vOjfcIILH0JBYSPHpT5fmul/5iuiRW7RI9PzJXy60cPXVooyMxyOGTr//fTE/8LvfBW6+WZvt8IaHgZ/9DPjJT0RYk+rqxJD3DTcAdrsIe9mcc444NjZOfU5RgKeeAp54YuL2l7+I+ZR//rO4SSaT6PU9//yJ21lniX9XixeLm8kktkAcHxchZnR04jYyIgLs/v0i1O3fP/HLWa5MJhHIFi8WPXqpwUeGsmPHxHnlLeUXYxw9Km7y/wWiYiUSlRcAa2pq8OCDD+K6POYadXZ2Gn4/YNk7V6dhb5EW1+js7MQdd9yR/PrSSy/Ffg2K9N5wg/jZsGpVAW/WYwg4dfi3lPvW6lEKRv6wvuqq0l1z1izRC9jRAfznf4q5clr/OW/bJnrgADEk+/a3a3s9QISO//5vMeT86U+LwPvhD4vh4O5uwGot/hojI0B/vwh9f/zjxOPz5wMf+pDYgeWGG8TcPrWYzWLuYOp2PidPAi+8IMLg7t0Tx5dfnvilRk01NSKkLlsmjtPdzjprIvTl+/kaHxeh7403pt5ef336x7M9n0hMDp6p9+fMEb2vshc23U0+P3eueM+cOeIX1dRb6mOy9218fHLP4+m3RGIihGc65vKa6d6TKvXvolLvl/L/2yKoGgA7Ojpw/fXXo6enJ6cQuHnzZgQCAbzwwgtqNkMzh0rwG6Ka15g3bx7mpcx/yrQYpxirVxcw9CvpHQBLSY9SMDIAvuc9pbsmIELfN74hgsIDD4jFGFp59lkxHAsAd9wBtLZqd63p2Gyix2zTJuA73xG9kI2N4vEvfhG4/vr8Alo8Lubz3Xcf8L//O1Ek3WQS8xr//u9F70Ipe69nzxalcS69FLj11onHX3lFlBhKHTaOxyeGmOVWe6lDkNXVImSazeL+2WeLObLnnisC37Jl6g/dn27WLBEcFy4UQZKoAqk+B3DLli3JPWvb2tqwevVqKIqC4eFhKIqCeDyOSCQCn8+HWCxWFvsH57MwQ8plNW+pr2FIegTAUq8AlmQP4NBQaa538iTw+OPifql/I62rEyuCv/tdoKtLuwB45MjEoo9rrpnoBSy1xYuBb31LlMH56lfF0HcoJG51dWJY+L3vFb8pnX++6OFKJCbm4T33nCg38+ij4pa6Uvwd7xCh76MfFbuYGMnZZ4tbqX/BIKKiqRoAAcBms6Gvrw8ulwsdHR3Jx1N3rkgkErBYLBgcHMSqgsYNSyufoCWHcvMNdKW4hiFVUg+gLI1UqgD45JNiwUt1tei5KbU77gD+67+AP/xBrFRVe95jIiH29X32WbGoprd3YmcPvdTXi+Hab30L+N73xJy9V18VC2J++MPcz3P55WI19Uc/KhZeEBGpTJM1/Q6HA/F4HF1dXVi1ahUSiUTyVl9fD6/Xiz179mB1weOGpdXQ0AAgt8UqcjFHvr1zpbiGIemxClivAHjq7xgvvliardLk8O+VV+pTvmP58omh2W9/W/3zb94shknnzAECAdETZRQXXCB6P/ftA4JBYONGsTJwugUadXViyNjtFtvNDQ+L8N7VxfBHRJrR9Nfljo6OSb2A5UqWt8mlVIt8jTXPCeCluIYh6bEKWK8h4HPPFZO8jx4Vc6ZkINTKo4+Ko54Tkjs6xKKM3/5W1CNMXVhQjB07RKgCgH//99KVuMnX7NliLmDq/t7Hjomha5NJzEHTer4bEdE0DFTV07hsp/7zztY7J4NZIfUQS3ENQ6qkIeBZswD59xaNan89PVYAn+7CCye2h9uwQZ26bS+9JBYijI8Dt90G3H578ecspXnzxC8fdXUMf0SkG0MEwM7OTr2bkJHZbIbVas0azuR+vm6325DXMCQ9A2CpewABYOVKcZTFmbXy6quA/CxdeaW218rmG98QPV2PPSaGbIvx5pui1MrBg2JBxfe/X9pSPkREM4QhAmDqAhGjkiE1EAikfY3c09jlcqV9TSgU0vwaZUXPVcCl7gEEJoZ9te4BlL1/l11W3D64ali2TOylC4iyKIX+XScSojdx1y4R3rXY6YOIqEJoNgdweHg4pwUNkUgkr23Q9OJwOGC1WtHV1QWHwzHleUVR4Pf74fP50q7OdTqdCAQCMJvNGBoamvI6Na5RdippCBiY6AEsVQA0SkFSjwf46U/F9+3xAHffnf85vv1tsfXZ7Nli0Yesq0hERPlLqKynpydRW1ubmDVrVl63cjAyMpKwWq0Jl8s17eMdHR0Z3282mxMAEgAS4XBYk2tkct555yUAJM4777yCz6G6eFzulplIHD9emmsuWSKu99xzpbleqt//Xlz77W/X9jrve5+4ztat2l4nHzt2TPxd/+pX+b3X55t473//tzbtIyIyKC1+fpsSCTVmZQsDAwOw2+2wWCyw2Ww59VKFw2Hs2LEDY6dvJ2Ng3d3d2LlzZ7IMSzweh9vtTi7kSMfv98Pj8aCpqQnBYFCTa2SyfPly7Nu3D+eddx727t1b8HlUdfLkxC4Jr72mfa/c8eNiEn6prne6aFT0Ai5YIIoXa1Ge5cQJUV7nrbeAp58Ww8BG8YUvAP/xH6Ln9//+T+xRnM3//I/Yam18HPjyl0WNPSKiCqLFz29VA+CaNWuwceNGtLS05PW+urq6kmyzVukMGQABYNEiMbk/FpsolqyVl18Wc9JMJhGUqqq0vd7pTpwQ4W9sDNi7VxQwVls4DDQ1ibl/hw7pUwMwnePHxV65jzwivvcHHxQrhaczPi5q4X3lK+Jrt1sMHXPRBxFVGC1+fqv6k6Gmpibv8AcA9Vr/0CdjK+U8QPmLRm1t6cMfIHo75dw1reYB6l0AOpO5c8VK4MsuE0WSr7xSfH3676HPPiu2j5Phb8MGhj8iIhWpugik0MLEsrQJVaglS0TPXCkCoJ4rgKWVK0Vv5549YncItRltAcjpli4VPX8f/KDY/9bhAN71LlEsed488djAgOgBXLBAFHpub9e71UREM4qqAbAcVvOSAenRA6hnANS6FIzRAyAgQuAjjwB33im2THviCXFLddNNQHd3bvMEiYgoL6qODzmdTmzbti3v9xm9EDRprJT7AetZBFrSshTMK68AQ0NiqFTvAtDZzJ0rFnS8+KLYLu6LXwQ+9zkRCJ9/HvjVrxj+iIg0omoAbG5uRjQazTsElkMhaNJQKfcDNsIQsOwB1GI3kNQC0DJYG91ZZwGf/KQIfv/1XyIIplsYQkREqlB1CHh4eBhOpxP9/f1Ys2YNbDYbGhoaMu5bWy6FoElDlTYEfNFF4vj882Lxg5oLGx59VBzf8x71zklERDOO6otARk/14iQSCYTDYZi4ao+y0SMA6jkE3NAgViAfOQLs369uKZhymP9HRES6UzUAyqLFra2tOW9VJgtBUwUrZQA0whDw3LkiBD7/vCh3olYAPHECkCvqGQCJiCgDVQOgxWLBxo0bsX79+rzeV6fnD2PSX6UNAQPApZeKAPjMM6Iwshp27waOHgVqaiaGmYmIiKah6iIQi8WScb5fOiwEXeEqbRUwMLG69Zln1DunnP931VXGKwBNRESGomoP4JYtWwp6HwtBV7hKWwUMiB5AQAwBq0XO/+MCECIiykLVboLdu3cXVAeQKlyphoDHxoCREXFf7wCodQ8gERFRBqoGQIfDAY/HU9B7N2/ejLq6OlRVVaGurg7333+/mk0jIytVAFSUiT1njRIADxxQp+dz717gb38TQ79XXFH8+YiIaEZTNQC63W688MILeb+vtbUVHo8HTqcTfX19WL9+PdavX48HH3xQzeaRUZUqAMrh3yVLgDlztL1WNtXVwLJl4r4aw8By+Pdd7wLOOKP48xER0Yym6hzAeDye93t27dqFQCAAt9uNu+++GwDQ0tKCtrY2tLW1FRQoqcyUKgAaZQWwdOmlogfw2WeL37aNw79ERJQHVXsAm5qasGPHDuzYsQM33HBDckh33bp1aXvzent7YTKZ4PV6Jz1utVrhcrk4p7ASyFXAb74JnDyp3XVkD6DeK4AlOQz8178Wfy4uACEiojyoGgBbWlqwadMm2O12BINBjIyMIJFIYPv27bDZbPjyl7885T2hUAhWqxVLZC9Qig0bNqC3t1fNJpIRLV48cV/LXkCj9QBefrk4Pvlkcec5ehSIRMR9BkAiIsqBqgFwdHQUg4OD2LBhA8LhMMbHxzE+Po5oNIq7774bW7ZsmdITGIvF0NTUlPacNTU1ajaRjGjOHGDBAnG/kgLgu94ljk88Udx5wmGxC8jZZwMrVhTdLCIimvlUDYCbNm3CwMAANm3ahNWrVycfr6+vh8vlQiwWm1IrUFGUjNvGcS/hClGKeYBGGwJ+xzvEcf9+4ODBws/zpz+J41VXAfz3QkREOVA1AI6MjEwKfqczm83s0aPplSIAGq0HcPFiYOVKcb+YXsCHHxbH972v+DYREVFFUDUA5hLu2KNH06rEAAgUPwx88iTwyCPi/jXXqNMmIiKa8VQNgIlEAk9k+EF2+PBhHJI/hHOU7+upTJViP2CjDQEDEwFw9+7C3r97t/gzq66eOBcREVEWqgZAl8uF6667Dt/5zncwPDycfPzw4cPYunUr6uvrcfvttycf37VrFwCxEGQ627Ztw5o1a9RsIhlVKfYDNmIP4KpV4lhoD2Dq8G9VlSpNIiKimU/VQtAWiwWbNm3C7bffjo6OjinPr169GpFIBLW1tdi5cye8Xi+sVisA4MUXX8QFF1yQfO3o6Cg8Hg/C4bCaTSSjqvQh4GeeEeVc5s/P7/0PPSSOHP4lIqI8qNoDCIhewMHBQaxduxaJRAKJRAL19fXYsmULBgcH0dLSgl/84hcIh8PweDwYHByEx+OB1WrFPffcg+HhYezYsQNNTU1oaWmZtj4gzUBaB8BEYiIAGmkI+PzzgaVLxVw+WcsvV2NjwB/+IO5fe63qTSMioplL1R5AyWq1IhgMTvtcfX09Nm3aNOX1Pp8Pra2tyUUi7e3tU15HM5jWAfDw4YldRozUA2gyifItv/yl2M0jn0LOTzwhvq8lSyaGkomIiHKgSQAshMPhwPj4OAYGBmCxWFBfX693k6iUtA6AsvdvwYKJotNG8e53iwD42GP5ve+BB8Txmms4/4+IiPKiSwC8/fbbkz19drsdt9xyS/K55uZmPZpEetM6ABpxBbB01VXiKPfzzdXvfieO73+/uu0hIqIZT5cAmLobyMDAALZu3Yr169fr0RQyClkGRqtVwEZcACI1NYkevH37gL17geXLs79ndBR49FFxf906bdtHREQzjuqLQPLV3NwMn8+ndzNIb6UaAjZiAFy0CHjnO8V9ua1bNqGQWARy8cUAp0sQEVGeNOsB3L17d9r6fgAQj8ehKAp6e3tRW1urVTOoXFTyEDAg5vHt2gUMDABtbdlfL4d/2ftHREQFUD0ADg8Pw263Zwx/qerr69Hf3692M6jcVHIPIADY7cD3vgds3y5K1mTaMvHkSeBXvxL3b7yxJM0jIqKZRfUAaLPZEI/H0dLSgjVr1sBsNqO/vx92ux1msxkAoCgKgsEgzGYz+vr61G4ClaNKD4DXXAPMmQO8+CKwZw9w4YXpX/uHPwAHD4rvhfX/iIioAKoGwJ6eHlgsFoTDYVTLSf2ntLW1TSrqvGHDBgwNDaGzsxNdXV1qNoPKkfxsvPGGmNumdlkTow8BL1oEvPe9wIMPivIumQKg7DH/yEdEaCQiIsqTqotAAoEAAoHAlPBnsVgwODg45fX19fVwuVzYunWrms2gcpS644sWvYBG7wEEJubzbduW/jUnTgD33SfuO53at4mIiGYkVQNgfX39tFu31dbWZtwZJBqNqtkMKkdz507sg1upAbC1VRwfegg4cGD61/zyl2L495xzgOuuK1nTiIhoZlE1ANbU1Ez7+OrVqxEKhdS8FM1EWtYCNPoQMACsWCGKQicSQLq5sbKG5qc/zeFfIiIqmKoBMJFIpH2usbER999//7TP5bpimGY4LQNgOfQAAsDHPiaOd98NjI9Pfm73blEmxmQC2ttL3jQiIpo5VA2Adrsd207NX9q6dSva2tpw+NRwnsPhgMPhwIMPPjjpPQMDAwyAJJxaJQ5FUfe8b74JHD0q7hu5BxAA/uEfRBB+7jng17+e/NzXviaOt94KXHBB6dtGREQzhqoBsLm5GdFoFHfddRdcLhcCgUCyzIvNZsPatWths9mwbt06dHZ2oq2tDddffz3cbreazaBypVUAlMO/c+YAZ5yh7rnVtngx8JnPiPv//M/AW2+J+/fdJwLhrFkTQZCIiKhAqtcB3LBhA3bt2oUNGzZgdHQUrXJiO8Qq4bVr12L79u0IBoNIJBKwWq3cB5gErQJg6vBvpgLLRtHZCfz4x0A0Kkq93Hwz0NEhntu4EbjkEn3bR0REZU+TreBWr16N1atXT3m8uroa4XAYgUAAoVAIDQ0NcLlcWjSBypHWAdDow7/SkiXAz38OXH+9qAn4wAPicbudvX9ERKQKzfYCzkTOBySaROshYKMvAEn1vvcBjz4KfPObwP79wE03ARs2APPm6d0yIiKaAXQJgETTKsUQcDmxWoE0K+eJiIiKkfMikG2ZdicokpbnpjLCIWAiIqKSyDkAtmtYd0zLc1MZ4RAwERFRSeQcAEdGRjRrRKYC0lRBOARMRERUEnnNAbzwwgthtVpRW1urysXj8TgikQhGtdj5gcqP3AmEAZCIiEhTeQXAaDSKaDSqeiNM5VCbjbQnewDV/oWgHPYBJiIiKqGcA2B9fT0aGhpgt9tVbcDjjz/ORSAkpA4BJxLqFW1mACQiIpok5wAYjUYxMDCAYDAIk8mEtrY2rFq1SpVG1HFojoCJADg2BrzxhnrbtnEVMBER0SR5DQE3NzejubkZAHDfffdhy5YtWLlyJVwuF5YsWVJwI7gIhAAACxaI/XpPnBC9gGoEwKNHgddfF/cZAImIiAAUUQi6paUFLS0tGB0dhc/nQywWQ2NjY0H7+g4MDBTaDJpJTCbRC3jwoAiAy5cXf07Z+1dVNbHIhIiIqMIVvRNIdXU1NmzYAAAYGhrCxo0bYTKZYLfbsXbt2pzOMd2+wVShUgOgGlJrAHKxEREREQCVt4Krr6/Hpk2bAIhePRkG1ZwvSDOc2rUAOf+PiIhoCs32Aj59vuAvfvELnHnmmXA4HFixYoVWl6Vyp3YA5ApgIiKiKTQLgKlS5wv29fUhGo1i5cqVaG1tLWrxCM1AaheDZgAkIiKaIuet4NRQXV2N9vZ2bNq0CXv27EFNTQ1uvfVW7Nixo5TNICPTqgeQpYaIiIiSStIDKO3YsQM+nw+BQACAKP/S19eHmpqanBeM0AzHIWAiIiLNaR4Ad+/eDZ/Ph76+PiiKkqz5Z7PZ4Ha70dLSonUTqJyovR0cF4EQERFNockQ8PDwMDZv3owLL7wQjY2N8Pl8GBkZwerVq5P3t2/fzvBHU7EHkIiIykAgEIDdbkdNTQ1MJhNqamrgdDoRCoWSr1EUBU6nU8dWpqdaD+Dhw4fR19cHn8+HSCQCQAzxWiwWuN1uOBwO1NfXq3U5mqk4B5CIiAwsFoslQ11nZyf6+/thNpuhKApCoVAy83i9XjidTlgsFp1bPL2iA+C2bdvQ29s7aV6f2WxGa2sr3G43izxTftgDSEREBuX3++F2u+Hz+eByuSY9Zzab4XA44HA40NjYiMbGRkQiEfT39+vU2swKCoDp5vU5HA643e5k/T+ivLEQNBERGZDs3fN6vVPC3+l6enrQ2NgIQKx5MKKcA+Dw8DB8Ph/8fv+UxRxOpxPt7e0FN6KqqgpjY2MFv59mEDUD4FtvAW+8Ie4zABIRURGcTidsNhs6OjqyvtZqtSaHfs3y55rB5BwALRYLTCbTpHl9LpcL1bJwb4F27dpV1PtphkkNgIlEcfv3yt6/2bMBFhwnIqICyc4vt9ud83ssFoth5/8BeQ4BW61WuN1uNDY2wmQyYWhoqOALx+NxKIqSU5KmCiJ/oTh5EnjzTWDRosLPlboApJggSUREFU3O48sn0JnNZtjtdq2aVLS8AqDFYsHIyAjC4XByCLhQiqLg0KFDiMViMJXRD2e/349gMIja2loAIsi63W5VxvjtdjucTmfauQWBQAA+nw/BYLDoaxnWokVAVRUwNgaMjKgTADn8S0RERYjFYsmj1WrN6T2Kohh2/h+QRwA0mUzo7e1VvQE2mw033HCD6udVm6IoaG5uRlNT06QVPamP+3y+oq4hw6TX64XH40FTUxPMZjMikQi6uroMvZpINSYTUFMjwpuiAMuXF34uLgAhIiIVWCwWxGIxdHV1weFw5PQeWR7GqHIOgMXO9UvHZrMV3ZtYCu3t7aitrZ0S8sxmMwYGBlBfX4/GxsasK4NyEYvFpp1n4PP5cv7glbXaWhEA4/HizsMeQCIqF4mEmPZC6S1cqNt0HrvdjlAohEgkgsbGRni93qy9e0YOf0AeATAcDmvWiGg0qtm51RAKhRAIBNIOvZrNZrhcLrjdbrS2thb1l+5yuTA4OIhYLAZFUWC1WtHU1ASPx2PoyaSqkkWbZQ9eoVgEmojKxZtvAmecoXcrjO3114ubFlSEjo4O+Hw+xGIxRCKR5Nw+q9WKtrY22Gy2nIeGjSLnreC03MXD6DuEeL1emM3mjGm/ra0NgJgjWCjZwxgOhzEyMoJEIoFwOAyfz1c54Q8QPYAAewCJiMgwwuHwlFG4SCQCj8eDxsZG1NTUoLu7W6fW5U+1reBmKrm1S7ZkLwNab28vVzYXS60AyDmARFQuFi4UPVyU3sKFul7ebDajv78fiqKgr68PwWAQoVAIyqm6tYqiwOPx4NChQ/B6vbq2NRcMgFkMDg4CyL70Ww77yn2QqQhqDwEzABKR0ZlMug1vUn7ktC855z8WiyEQCMDj8QAAuru74Xa7DT9yl/MQcKWSgU6WfclE/mUzBBZJ7SFgzgEkIiKNWCwWdHR0TFrPEAqFdGxRbhgAszh0qhcqn4Ud8WKDC8RvFH6/Hx6PJ1mBvGJwDiAREZUZi8WS7BXUcuGsWhgAs8gneMmQKAtGFiIWi8Hj8SAWi6G1tRVerxe1tbVobGxEIBAo+LxlRY0h4EQCOHhQ3GcAJCKiPAUCgbwXdjY0NAAwfgkYgHMAs1KjNy+fawUCgSmTRx0OB6xWKxoaGtDR0ZHz5NJjx47h2LFjya/Lod4iAHV6AF9/HTh6VNw/++zi20RERBWlt7c3r71/U61Zs0bl1qiPATAL2QNYV4J5ZD09PWlXG8s5BvlMLu3q6sKdd96pdjO1J/+siwmAr74qjgsXcmI1ERHlLRQKobOzM6/3yHrB6TZtCAQC2LlzZzJbyNqCPp8PdXV1iEajRe8qlisOAefoULErUnOQrdSMLDyZ628knZ2dGB0dTd7OPffcottYErIHsJg/cxkAly4tvj1ERFRR5GYMshJILkKhEEKhUNpROjmNy+v1JoNfd3c3vF4vvF4vrFYr/H5/yRaQlHUPYGr9HTVYLJYpIayQcfxcVgwXoqmpCUDuq4vmzZuHefPmJb826bSFTt7kn99bb4nbggX5n4MBkIiICiR/zvp8vpy2eFUUBU6nE1arddpawIqiYOfOnZPCocVigcfjSS4YiUQiyd2/SqGsA6DT6VQ1AFqt1ikrd/IJc7ItWk3+TD1vLBYzfI2hgi1ZAlRVAWNjwMgIAyAREZVU6tavTqcT/f39aV8bi8Vgt9tRW1uLgYGBaV/T19c3ZfRucHAQVqs12fHU0dFR0o0kynoIWG6XptZtumXbckVPLit75YKRfHsAFUVBd3d3Xt2+xaw0NjyTqfhhYAZAIiIqkBzKlbmgoaEB3d3dyaFhYKJqR0NDA2w2G6LRaNoOIJfLNaXTJhKJZNxiVmtl3QNYCvIvLJeeRvmafDeElvX+zGYzRkZGcnqPVsPMhlFbK8q4FLoQRJaAYQAkIqI8xGIx2Gy2ZG9cf38/IpEIfD4f7HY74vE4FEWBxWJJBr98R+Rkh4+c268HBsAsZDrP1uMmw18xw7LZfhNIDaEzdvhXKnYlMHsAiYioABaLZcqQr9VqVXV1rhxi1rMHsKyHgEvBbDbDarVmDYBypVAhNYPkByvTHANg4jcGm81WFkUmi8IhYCIimqFCodCU0UJFUfIuPF0MBsAcyDpAmXbikOEt02qhdHP8bDZb1vAHIPnbR6lqBOmq2GLQDIBERGQQdrs9Ge5isRgikciU1b5+vx+tra0laxMDYA7kThxdXV3TPi9Tu8/nS9sz53Q6YbfbUVNTM+18Qo/HA6fTmbYNgUAgOSl1xg//AhwCJiKiGSESiUzqAAoEAujo6JhUYzAUCsFisZR0dI8BMEdyaffpQ7yKoqC5uRkdHR059f4pijLtcLLNZoPdbofdbkckEpn0XHd3N9rb2+Hz+Uq6RFxXxfQAjo9zEQgRERmCrA2oKAo8Hg9sNhu8Xi+amprg8Xjg8XigKEra3UO0wkUgOTKbzQiHw+ju7obT6Uyuwo3H4/B6vVkncnq9Xng8HjQ1NaVdJexyuWCz2eDxeCaFQJvNhnA4XBk9f1IxcwDjcRECAeDMM9VrExERUQGm2x1E7+lcDIB5KrQHzuVy5VRNfLrVRxVJDgEXEgDl8G9tLTBnjnptIiIimiE4BEzGJHvuXnst//dy/h8REVFGDIBkTDK8yTCXDwZAIiKijBgAyZhkeHvtNbEncD4YAImIiDJiACRjknMAE4n8VwLLAHjWWeq2iYiIaIZgACRjmj17IgTmOwz88svieM456raJcuZ2u1FTUwOn05l1Fx0iIio9BkAyrkLnAe7fL47LlqnbHsqZz+dDT08PAoEAGhsbpy1+PlP5/X44nU643W643W44nc60uwAVIhaLwe12o6GhASaTCQ0NDXA6nRl3KiIiOh0DIBmXHMLNNwAeOCCODIC6cjgcsFgsUBQF7e3tejdHc4qioLGxEeFwGP39/fD5fMkg7PF4Cton/HShUCgZLKPRKBKJRLJslNPpRGNjI3tciSgnDIBkXLIHUO7qkSsGQMOQoScQCMz4XsD29nbU1tZOKe5qNpsxMDCAvr6+ojZ6j0Qi8Pl8CAaDk4rJW61W9Pf3w+v1IhKJwG63F3wNIqocDIBkXIUMAY+NAa+8Iu4zAOoutfh5X1+fji3RVigUQiAQgMfjmfZ5s9kMl8sFt9tdcBD2eDwZi8R3dHTAYrEkh4iJiDJhACTjKmQI+OBBsQ2cycQyMAZgNpuT2yTO5B1uvF7vpO91Om1tbQBQUC9gJBLJaStIud1UMT2NRFQZGADJuAoZApbDv0uXipXEpDvZGxUKhWbkMLCiKAiFQlkDmny+t7c372sMDg4iFAql7WGUUoeGU/cTJyI6HQMgGVchQ8Cc/2c4DocjeX8m9kwNDg4CQNYAaDabARQWzKLRKGKxGLq7uzMu8qitrU3ej+dbP5OIKgoDIBlXIUPADICGJENgIb1fRicDXWr4SkeGxHxDYFtbW3KIOVPQTA19uQwZE1HlYgAk4ypmCJgB0FDk/LdIJDLjypQcOnQIwEQPXy7y7Z2zWq0YGRlBMBjM+LrUYMkASESZMACScckAODICHD+e23sYAA0pdRj49DIp5S6feY0yJGoVgmUPa+rqayKi6XCWPBlXTQ1QVSVKuxw8CJx3Xvb3MAAalsvlgt/vRyAQSK5WLUYsFlN1UYnZbC6o18woc+0URUnuBqLGny8RzWwMgGRcs2YBZ58ttnY7cIABsMzZ7Xb4/X7EYjFEIpFJK1bzFQqFVC94bDabMTIykvf7ZAitk3tX60TuttLf35/XcDQRVSYOAZOxydAn9/fNhgHQkBRFgcfjSQaTYoeBbTYbEomEqrdCwl8qORdQD7IQtc/nmzTcTkTaCQQCsNvtqKmpgclkQk1NzZS9vxVFgdPp1LGV6TEAkrGde6445hIAx8aAffvE/eXLtWsT5UXukevxeJJz02bSriCF9LblsmI4V7FYDE6nEx0dHZz7R1QCsVgMjY2N6OrqgtvtxtDQEBKJBIaGhtDW1ga3252s2el0OlX9964mBkAytnwC4CuvACdPiqFj9gAaggx/DocDLpcruRpYFk+eCfL5z10OF6s5RGu329Ha2sp5f0Ql4Pf70dDQALfbjXA4DIfDkfz3bDab4XA4EI1GEQqF0NjYqMl0FbUwAJKxyQAoe/Yy+dvfJt7DXUB0J8Of1WpNhhOr1ZpcaDFTVgM3NDQAyG1lr1wwolaPgN1uh81mmzF/lkRGFgqF4Ha74fV6s/a29/T0JMsyZdoiUk/8KUnGls8cwJdeEse3vU279lBOZPgzm81T9gB2OBzo7u5OrlgtlFFWAcv35NIW+ZpiFsBIbrcbFouF4Y+oRJxOJ2w2Gzo6OrK+NvWXXaMuymIAJGPLZwhY9gCef7527aGsZPiLx+MIh8NTnne73eju7gYgJlEXsmghEomgsbGx6LamKnQVsPztPlsPoAx/ahRoln9+6cKfFkPNRJXM7/dDUZTk3ua5sFgshi7IziFgMrZCAiB7AHUjw18sFsPAwMC0AcRisSR7wArtvbJarYZZBWw2m2G1WrMGQLlncD4/QKYTCAQQjUYz/tn19fXNuB1XiPQkRzLyCXRms9mw8/8ABkAyOhkAX3sNOHYs82vZA6ir1PDX39+fcZhThqBQKKTqMK5eOjs7ASDjsLb8AZJp7lC2hTGRSATBYDBrcA6Hw6oMMxORIH+hyucXK0VRDDv/D2AAJKOrrQXmzRP3ZY2/dOQcQAbAkksNf16vN+uwbmtra/L+TCgJ43A4YLVa0dXVNe3ziqLA7/fD5/OlHZZ1Op3JmmLTheJYLIaurq6s4S8WiyV7G4lIHbLnL92/8ekYvSg7AyAZm8mU+zAwh4B1I8Ofw+HIaYK0LJcAzJxtywYGBgBMHeJVFAXNzc1Z6/TJ3j9FUab0MiiKArvdjkAgAJPJlPHW0NBg2LpjROVKDuXK+ce5lLEycvgDAFMikUjo3QgqjeXLl2Pfvn0477zzsHfvXr2bk7v3vQ/44x+Bn/8cuPXW6V9z/Dgwfz6QSACvvgqcdVZp21jhampq0NTUhGAwmPN7ZDFVRVEQDAYNPVSSj+7ubuzcuTMZwuLxONxud9bvz+/3w+PxTPvn6Ha74ff7c26Dy+Xi6uAyk0gAb76pdyuMbeFC0Segl4aGhim/nFmtVrS1tcFms2k67UKLn98MgBWkbAPgJz4B3Hsv8O1vA6fmWk3x/PPAxReL/yFef13f/yWIiPL0xhvAGWfo3Qpje/11YNEi/a6vKAra29vTzvU1m83o7OzMaRQkX1r8/OYQMBlffb04Dg2lf000Ko4NDQx/RESkOlnXdGRkJLnvduowr9zzXG4DZ3SsA0jGJwPg8HD616QGQCKiMiMHLyi9hQv1boFgNpvhcrmSc3pjsRgCgUAy+HV3dycLtRsZAyAZ34oV4phrDyARUZkxmfQd3qTCWSwWdHR0wOFwJLeGDIVCWbeL0xuHgMn4ZA/giy8CY2PTv0YGwJUrS9MmIiKiFBaLJRn6ptsFyWgYAMn4li8HZs8GTpxIXwpmzx5xZA8gERGpIBAI5LX6HkCyB9DoJWAABkAqB1VVE7X9ppsHOD4+MTzMAEhERCro7e0teB7fmjVrVG6N+hgAqTzIYWA51JvqpZeAo0eBuXNZBJqIiFQRCoXyLqoua3hm2g3J7/ejpqYmWQdVKvXqYQZAKg+XXCKOzzwz9bmnnpp4zWyuayIiouLEYjEoipLXtoqhUAihUCjj7kaBQABerxetra3JXYIURUF3d3fJew3505LKw2WXieNf/zr1uaefFse3v7107SEiohlLbvXm8/lyWs2rKAqcTiesVmvGQtDBYBDRlJGsSCQCj8cDs9msSQHpTNgDSOVBhjsZ9lLJHsDLLy9de4iIaMZK3Y7R6XRmfK3c1rK2tja5J3g6p2/RaDabYTabddkTnQGQyoPsARweFnsmpWIAJCIiFcmhXFnOpaGhAd3d3cmhYUAEP4/Hg4aGBthsNkSj0bxW/0YikeSQsB4YAKk8nHUWcOaZYsf0Z5+dePz48Yl5gQyARERUpFgsBpvNlhyS7e/vR39/P6LRKOx2O+rr62EymWC326EoCqLR6JSevWxCoRAGBwdLPuybinMAqXy8613AwAAwOAg0NorHnngCOHYMqKubWClMRERUIIvFgv7+/kmPWa3WvENeOoFAILmdnJ7YA0jl46qrxPHPf5547LHHxPHKK8VeSkRERAbV3d0Ni8UCm8026XG/349YLFbStjAAUvmYLgA++qg4vvvdpW8PERFRjvx+P3bu3Dkl6IVCIfT39xdcdLpQDIBUPmTIe/554OWXxb7A27eLx669VrdmERERZRKJRBAOh9HT04OdO3eipqYGbrcbTqcTTqdTteHlfHAOIJWP2lrgiiuAxx8HfvUrsegjHgfM5oneQSIiIoPp7e1Nhjy56tfv98NisSAcDpe89w9gDyCVm498RBx/9CPgnnvE/Ztu4g4gRERkWKeXevF6vRgZGdEt/AEMgFRuPvEJsefvo48CP/iBeOz22/VtExERUZlhAKTycu65wJe/PPH1xz4GvOc9+rWHiIioDHHcjMrP174m6gAeOQI4HHq3hoiIqOwwAFL5MZmAD35Q71YQERGVLQ4BExEREVUYBkAiIiKiCsMASERERFRhGACJiIiIKgwDIBEREVGFYQAkIiIiqjAMgEREREQVhgGQiIiIqMIwABIRERFVGAZAIiIiogrDAEhERERUYRgAiYiIiCoMAyARERFRhTElEomE3o2g0pg7dy5OnDiBWbNmYdmyZXo3h4iIiHJw4MABjI+PY86cOTh+/Lgq52QArCBVVVUYHx/XuxlERERUgFmzZmFsbEyVc81W5SxUFubPn4+jR4+iqqoKS5cuVfXciUQC+/fvx7nnnguTyaTquYlKjZ9nmkn4eS5/r776KsbGxjB//nzVzskeQFLF4cOHUV1djdHRUSxZskTv5hAVhZ9nmkn4eabpcBEIERERUYVhACQiIiKqMAyApIp58+bh61//OubNm6d3U4iKxs8zzST8PNN0OAeQiIiIqMKwB5CIiIiowjAAEhEREVUYBkAiIiKiCsMASERERFRhuBMIFczv9yMYDKK2thYAEI/H4Xa7YbPZdG4Z0WSxWAxerxehUAixWAwWiwVWqxVtbW1wOBw5nYOfdyoHfr8fiqKgo6Mjp9fyM13BEkR5GhkZSVit1oTL5crpcSI9BYPBhM1mS4TD4eRj4XA44XA4EgASVqs1EY1G076fn3cqFyMjIwkAiY6Ojqyv42eaGAApbw6HI2Gz2aZ9bmRkJGE2mxM+n6/ErSKaSga9dLxebwJAwmKxpH0NP+9ULuQvNdkCID/TlEgkEpwDSHkJhUIIBALweDzTPm82m+FyueB2u6EoSmkbR3Qaj8eD/v7+tM93dHTAYrEgFovB7XZPeZ6fdyoXgUAgp88gP9MkMQBSXrxeL8xmc8Y5Im1tbQDE/BIivUQiEVgslqyv83q9AKb/vPLzTuUiGAxO+0vM6fiZJokBkHKmKApCoVDWH6ry+d7e3lI0i2hag4ODCIVCaXs6JKvVmrwfiUSS9/l5p3Lh8Xiyfs4BfqZpMgZAytng4CAAZP3Pw2w2A5j8w5So1KLRKGKxGLq7uxGLxdK+Tq6ABMQqSImfdyoHkUgEdXV1OfV28zNNqRgAKWfyP4PUH5jpyP9g+B8I6aWtrS051JXpB15q6Et9HT/vVA66urpyKvkC8DNNk7EOIOXs0KFDACZ+O8xF6g9XolKyWq0YGRnJ+rrUH3CpAZCfdzK67u5udHZ25vx6fqYpFXsAKWf5rAiT/8FkGnojMgI5z8nlck16nJ93MrJYLIZDhw5NmsOaDT/TlIo9gJQz/iZIM42iKAgEAgAmVgNL/LyTkWUrcTQdfqYpFXsAKWfyt8e6ujp9G0Kkkvb2dgBAf3//lGExft7JqPx+f04lX07HzzSlYgCkvMl5JETlTBbE9fl8GfcD5uedjERRFITD4aL26+VnmgAGQMpDPhOHpVxWmxGVWiwWg9PpREdHx5S5fxI/72REHo9nynSFXPEzTakYACln+fxHIIcaCvkPh0hrdrsdra2tGX+Q8vNORhMKhdDY2Fjw54yfaUrFRSCUs4aGBgC5rQqTk4352yMZjd1uh81mg8/ny/g6ft7JaHw+X94LP1LxM02pGAApZ7JGWi6lBORr8ilRQKQ1t9sNi8WSNfwB/LyTsYRCIYRCIdTU1GR9bXd396R9fGU9TH6mKRWHgClnctJxtt8e5X8cuWxNRFQq3d3dAJA2/CmKMukHIz/vZCQ2mw0jIyMZb3JKQ0dHx6THU88B8DNNAgMg5cxsNsNqtWb9z0PuN1lImQIiLQQCAUSj0Yw9f319fZM+2/y800zDzzSlYgCkvMhth2Tx3OnIOSrpVlcSlVIkEkEwGMw67BsOh6cMd/HzTjMNP9MkMQBSXhwOB6xWK7q6uqZ9XlEU+P1++Hw+rh4j3cViMXR1dWUNf7FYLNnrkYqfdyonuczt42eaJAZAytvAwACAqcMDiqKgubk5Y201olJRFAV2ux2BQAAmkynjraGhIe1qR37eqVzs3LkTgFgwkgk/0wQApkQikdC7EVSeuru7sXPnzuQPzng8DrfbXVSFeiK1uN3uSSshs3G5XBl7Cvl5JyPyeDzJz3lqD6DsvWttbU37ueZnurIxABIRERFVGA4BExEREVUYBkAiIiKiCsMASERERFRhGACJiIiIKgwDIBEREVGFYQAkIiIiqjAMgEREREQVhgGQiIiIqMIwABIRERFVGAZAIiIiogrDAEhERERUYRgAiYiIiCrMbL0bQERUzvx+Pzwej2rn83q9cLlcqp2PiGg67AEkIipCNBqFoiiwWCzo6elBOBzGyMhI8jY0NITa2trka1KfGxkZQTgcRk9PDywWCxRFgaIok86vKAqcTqeqIZOIiD2ARERFiMVisNlsCAaD0z5vNpthtVoRi8VgsVimPG+xWGCxWOBwONDY2IhoNDrp+b6+PgQCAQBAW1sbrFar+t8EEVUc9gASERVBURR4vd6cXltbW5vx+c7OTsTj8UmP2Ww2WCwWWK1Whj8iUg0DIBFRkdQKZg6HY8oQsMViQTQaRTgcVuUaREQAAyARUVHMZrPeTSAiyhsDIBGRgTBQElEpMAASERVhuoUdRjofEdF0GACJiIqQ6wIQvc5HRDQdloEhIjKo7u5uBINBxGIxxONxDAwMTFlw4na7EYvFkq8Jh8OwWCwIhULo7++H2WxGJBJBPB5HZ2cnHA5H8r3yNQCS7/d6vbDZbFnbFovF4PV6EY/HUVtbi1gsBrPZDK/Xy15MojLAAEhEZFA2mw1msxkej2fK6mDJ6XQiFotNek13dzcsFgt8Pl/ydYFAAE6nEz6fDy6Xa9rX+P1+2O129Pf3TwqKpwuFQnC73fD5fJPCot/vR2NjI/r7+3MKkUSkHw4BExEZlNVqhcvlyrg1nM1mg8vlQmtrKwAR9GRh6VQOhwMWiwUejwehUGja17hcLlgsFrS3t6e9XiAQgN1unxL+5PtbW1tht9sRi8Xy/XaJqIQYAImIZgC5eri3tzdt753NZksWrk73GqvVCkVR0ga49vZ2WK3WtD18css6zmUkMjYGQCIig6urq8v5tZmGXmVIzFS4Wu5WMt2Qc3d3NxRFQVtbW9r3y63t/H5/bg0mIl0wABIRzSANDQ1ZX7NmzZqCzt3b2wsg+84nchFIunmLRKQ/LgIhIppBclmBW2ix6UgkAgDo7+9HMBhM+zo5fByPx1nYmsigGACJiCgvHo8nY9Dk/D8i4+MQMBER5YVDu0TljwGQiIhyInv94vG4zi0homIxABIRUU7cbjeAibmAmbAOIJGxMQASEVFOOjo6puweko6sB0hExsQASERUIoUOnR46dKio9xci3bX6+/sRi8XQ3d2d9r2hUAh2u12rphGRChgAiYg0pChKcsg0EokUtIBCDqdmGlZNvUYxr8l2LavVinA4DJ/PN20vXyAQQCQSybh9HRHpz5RIJBJ6N4KIaCbxeDzJnTCmC3yyNl5nZyc6OjrSnqexsRGxWGzSOcxmM3p6epJbuaV7jcvlSpZjUes1032fcl9huYOI0+nMuBsJERkDAyARERFRheEQMBEREVGFYQAkIiIiqjAMgEREREQVhgGQiIiIqMIwABIRERFVGAZAIiIiogrDAEhERERUYRgAiYiIiCoMAyARERFRhWEAJCIiIqowDIBEREREFYYBkIiIiKjCMAASERERVRgGQCIiIqIKwwBIREREVGH+f7Vwy/uu4bekAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "times=np.arange(t0,tf+dt,dt)#[0:-1]\n",
    "\n",
    "tep=100\n",
    "\n",
    "plt.plot(times, mx_a,color='red', label = r'$S_x$')\n",
    "plt.plot(times,mz_a,color='blue', label = r'$S_z$')\n",
    "plt.annotate(r'$\\mathrm{\\lambda = 0.2 }$', xy=(0.5, 0.5), xytext=(0.3,0.1), xycoords='axes fraction',fontsize=fs )\n",
    "plt.ylabel(r\"$\\rm{ Magnetization }  $\",fontsize=fs)\n",
    "plt.xlabel(r\"$\\rm{Time} $\",fontsize=fs)\n",
    "plt.tick_params(axis= \"both\", labelsize =fs)\n",
    "plt.legend(fontsize=fs,frameon = False,loc= (0.7,  0.01))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "221d8d5e-1309-4665-a25e-ef97af22cdc2",
   "metadata": {},
   "source": [
    "## Computational solution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab7bb482-dfef-400a-bd8e-ebd900cde970",
   "metadata": {},
   "source": [
    "**Excersice:** Implementing the heun's scheme, find the dynamics of a single spin under a localized field as in the analycal solution. \n",
    "Make a plot comparing your answer "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c0bf1a2",
   "metadata": {},
   "source": [
    "Send corrections to <code>bnikolic@udel.edu</code>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}