{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "262a5ec8-2553-4237-ab62-319b6ca22089",
   "metadata": {},
   "source": [
    "# Example-48: Advance (Computation of phase advances)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "19c9ce37-2486-4548-aeb0-259ec857b96a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this example phase advances between BPMs are computed and compared with ELEGANT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "70bf6e2e-b566-495b-a0ab-edefefce88d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import\n",
    "\n",
    "from pprint import pprint\n",
    "\n",
    "import torch\n",
    "from torch import Tensor\n",
    "\n",
    "from pathlib import Path\n",
    "\n",
    "import matplotlib\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib.patches import Rectangle\n",
    "matplotlib.rcParams['text.usetex'] = True\n",
    "\n",
    "from model.library.corrector import Corrector\n",
    "from model.library.line import Line\n",
    "\n",
    "from model.command.util import chop\n",
    "from model.command.util import select\n",
    "\n",
    "from model.command.external import load_sdds\n",
    "from model.command.external import load_lattice\n",
    "\n",
    "from model.command.build import build\n",
    "\n",
    "from model.command.tune import tune\n",
    "from model.command.advance import advance\n",
    "\n",
    "from model.command.layout import Layout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "77bdc7c3-ad18-44fe-b8d7-81c022068e83",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load ELEGANT twiss\n",
    "\n",
    "path = Path('ic.twiss')\n",
    "parameters, columns = load_sdds(path)\n",
    "\n",
    "nu_qx:Tensor = torch.tensor(parameters['nux'] % 1, dtype=torch.float64)\n",
    "nu_qy:Tensor = torch.tensor(parameters['nuy'] % 1, dtype=torch.float64)\n",
    "\n",
    "# Set phase advances from lattice start to each BPM\n",
    "\n",
    "kinds = select(columns, 'ElementType', keep=False)\n",
    "\n",
    "mu_qx = select(columns, 'psix', keep=False)\n",
    "mu_qy = select(columns, 'psiy' , keep=False)\n",
    "\n",
    "mu_qx:Tensor = torch.tensor([value for (key, value), kind in zip(mu_qx.items(), kinds.values()) if kind == 'MONI'], dtype=torch.float64)\n",
    "mu_qy:Tensor = torch.tensor([value for (key, value), kind in zip(mu_qy.items(), kinds.values()) if kind == 'MONI'], dtype=torch.float64)\n",
    "\n",
    "positions = select(columns, 's', keep=False).items()\n",
    "positions = [value for (key, value), kind in zip(positions, kinds.values()) if kind == 'MONI']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "743ea4ef-dd8c-4e47-a7ed-5b42e94a3857",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build and setup lattice\n",
    "\n",
    "# Load ELEGANT table\n",
    "\n",
    "path = Path('ic.lte')\n",
    "data = load_lattice(path)\n",
    "\n",
    "# Build ELEGANT table\n",
    "\n",
    "ring:Line = build('RING', 'ELEGANT', data)\n",
    "ring.flatten()\n",
    "\n",
    "# Merge drifts\n",
    "\n",
    "ring.merge()\n",
    "\n",
    "# Split BPMs\n",
    "\n",
    "ring.split((None, ['BPM'], None, None))\n",
    "\n",
    "# Roll lattice start\n",
    "\n",
    "ring.roll(1)\n",
    "\n",
    "# Set linear dipoles\n",
    "\n",
    "for element in ring:\n",
    "    if element.__class__.__name__ == 'Dipole':\n",
    "        element.linear = True\n",
    "\n",
    "# Split lattice into lines by BPMs\n",
    "\n",
    "ring.splice()\n",
    "\n",
    "# Set number of elements of different kinds\n",
    "\n",
    "nb = ring.describe['BPM']\n",
    "nq = ring.describe['Quadrupole']\n",
    "ns = ring.describe['Sextupole']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "333d487b-3cb8-4e43-806e-adb71779d979",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "# Compute tunes (fractional part)\n",
    "\n",
    "guess = torch.tensor(4*[0.0], dtype=torch.float64)\n",
    "nuqx, nuqy = tune(ring, [], alignment=False, matched=True, guess=guess, limit=8, epsilon=1.0E-9)\n",
    "\n",
    "# Compare with elegant\n",
    "\n",
    "print(torch.allclose(nu_qx, nuqx))\n",
    "print(torch.allclose(nu_qy, nuqy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "23a2d707-1943-4bfa-bc39-3fcedddff0a0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([16])\n",
      "torch.Size([16])\n"
     ]
    }
   ],
   "source": [
    "# Compute phase advances\n",
    "\n",
    "muqx, muqy =  advance(ring, [], alignment=False, matched=True, guess=guess, limit=8, epsilon=1.0E-9).T\n",
    "\n",
    "print(muqx.shape)\n",
    "print(muqy.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "04602bd2-7281-4751-bbf0-ec479da4941c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "# Check total tunes\n",
    "\n",
    "print(torch.allclose(muqx.sum()/(2*torch.pi), torch.tensor(parameters['nux'], dtype=torch.float64)))\n",
    "print(torch.allclose(muqy.sum()/(2*torch.pi), torch.tensor(parameters['nuy'], dtype=torch.float64)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2d4156d5-a926-4504-aee6-70b20443b6fc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABjUAAAGGCAYAAAAzegNcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABf40lEQVR4nO3df3Rb933f/xf4C/whkiAlUaRsSRYo25KyNA5FJbaanbVH0nbWnXXNmRR3XetKTSTtbN1Z53bWvHYnR1tP9bW6dd3O1kZyYqlc08Ym++MkaZPGdNukmexaEZM6qSTLIizJtihREgmAEkXwB+73j08vgQtcgAAI4ufzcQ4O8ePi8gLkJS7xwvv99liWZQkAAAAAAAAAAKDE1RR7AwAAAAAAAAAAADJBqAEAAAAAAAAAAMoCoQYAAAAAAAAAACgLhBoAAAAAAAAAAKAsEGoAAAAAAAAAAICyQKgBAAAAAAAAAADKAqEGAAAAAAAAAAAoC4QaAAAAAAAAAACgLBBqAAAAAAAAAACAspBzqHH06FEFg8HFy4FAQCdOnNDQ0JBOnDjhuA0AAAAAAAAAAGC5PJZlWdneaWRkRDt27NDk5KR8Pp8kaceOHTp//rwkE3AcPXpUg4ODed1YAAAAAAAAAABQvXKq1AgEAvL7/Y7L8fx+v4aHh5e3ZQAAAAAAAAAAAHGyDjWGhoa0b98+x3XDw8Pq7Ox0XNfZ2amRkZHlbR0AAAAAAAAAAMDfqctm4WAwuNhuKvF6NxMTE67XRyIRRSKRxcvRaFQTExNavXq1PB5PNpsEAAAAAAAAAADKnGVZmpqa0vr161VTk7oeI6tQ45VXXtHhw4czXj5V2HH8+HEdO3Ysm28NAAAAAAAAAAAq3HvvvaeHH3445e0ZhxrDw8P61Kc+5Xqbz+dLqsqYmJhwreqQpOeff17PPvvs4uVQKKSNGzfqvffeU1tbW6abBAAAAAAAAAAAKkA4HNaGDRvU2tqadrmsKzVsgUBAx48f19NPP609e/bo5MmTScv39/e7rsfr9crr9SZd39bWRqgBAAAAAAAAAECVWmpERcahxp49exyXjxw5oiNHjsjv9yctGwgE1N/fn7JSAwAAAAAAAAAAIFupp22kEAwGdeLECUnSCy+8oJGREUnS4OCgjh49qqGhIZ08eVKDg4P53VIAAAAAAAAAAFDVPJZlWcXeiHA4rPb2doVCIdpPAQAAAAAAAABQZTLNCbKu1AAAAAAAAAAAACgGQg0AAAAAAAAAAFAWCDUAAAAAAAAAAEBZINQAAAAAAAAAAABlgVADAAAAAAAAAACUBUINAAAAAAAAAABQFgg1AAAAAAAAAABAWSDUAAAAAAAAAAAAZYFQAwAAAAAAAAAAlAVCDQAAAAAAAAAAUBYINQAAAAAAAAAAQFkg1AAAAAAAAAAAAGWBUAMAAAAAAAAAAJQFQg0AAAAAAAAAAFAW6oq9AUAlGB+XQqFib0Xm2tulrq5ib0XhlNvPR6q+n1E54veq9PEzAoCVUW5/X6vtb2u5/Xyk6vsZlSN+r0pfuf2Mqu3nAyC/CDWAZRoflz75SWliothbkrnOTumP/qg6DiDK8ecjVdfPqBzxe1X6+BkBwMoox7+v1fS3tRx/PlJ1/YzKEb9Xpa8cf0bV9PMBkH+EGsAyhULmwKGhQWpsLPbWLG1mxmxvKFQdBw/l9vORqu9nVI74vSp9/IwAYGWU29/XavvbWm4/H6n6fkbliN+r0lduP6Nq+/kAyD9CDSBPGhul5uZib0VmZmeLvQWFV04/H6k6f0bliN+r0sfPCABWRjn9fa3Gv63l9PORqvNnVI74vSp95fQzqsafD4D8YVA4AAAAAAAAAAAoC4QaAAAAAAAAAACgLBBqAAAAAAAAAACAskCoAQAAAAAAAAAAygKhBgAAAAAAAAAAKAuEGgAAAAAAAAAAoCzUFXsDkKPxcSkUWno5n09au3bFN2fFrcTjzWSdlfL8JejslOrrk68fvz+u0EzsOfE1+rS2JfXjX3L5HJ7jvG9DmUr1M0qynN/jct0HVvoxp7t/mSvIvs9+vywZ7fsrse9W076Rz8daTc9bLniuc5PpY5U49s0Qx76lL2+vfxLHvtneN939y5zb71XiPiSl34+WXJ79PmcF+Z93OevMZb2liOOxwuH5KShCjXI0Pi598pPSxET65To7pc99TurokOrK+Ee9Eo83k3VWyvOXoKND+rVfk3p6nNeP3x/XJ7/0SU3MmOeks6lTn/snn1NHU4fqapIf/5LL5/Ac530bylSqn1GS5fwel+s+sNKPOd39y1xB9v07E+z3y5DRvr8S+2417Rv5fKzV9Lzlguc6N5k+Volj3wxx7Fv68vb6J3Hsy7HvIrffq8R9SEq/Hy25PMe+OSvI/7zLWWcu6y1FHI8VDs9PwfHslaNQyOwkDQ1SY2Pq5To7peZmaX6+vHeUlXi8mayzUp6/BKtWSe3tkmU5rw/NhDQxM6GG2gY11jWqs6lTzQ3Nikajro3qUi0/H503B1c5PMdLrjPbbShTqX5GSZbze1yu+8BKP+Z09y9z2e77Oe137PfLktG+vxL7bjXtG/l8rNX0vOWC5zo3mT5WiWPfDBXk9Y9j32XJ2+ufxLEvx76L3H6vEvchSWn3oyWXZ7/PWUH+513OOnNZbynieKxweH4KjmevnDU2mh0hFa/XfF3yVaJMrMTjTbfOSnv+EqR6WI11jWqub5a31jx+S+kff9LyiSvO4Tlecp3ZbkOZyvhhLOf3uFz3gZV6zJncv8xlvO8vZ79jv1+WjB7KSuy71bRv5POxVtPzlgue69ws9Vgljn2zVJDXP4nXwGVY9uufxLFvtvfN5P5lzu1h2fuQpIz2oyWXZ7/PWUH+581lnctZbynieKxweH4KhkHhAAAAAAAAAACgLBBqAAAAAAAAAACAspBV+6nh4WFJUjAY1Llz5/T000+rr69PkjQyMiJJ6uvrUyAQUDAYXLwNAAAAAAAAAABgubKq1Ni/f786Ozu1b98+9fb2av/+/Yu3nTx5Ujt27JDH49GRI0fk9/vzvrEAAAAAAAAAAKB6ZVWpMTg46Ki+8Pl8i+d37NihycnJpOsBAAAAAAAAAADyIatQY8+ePYvnBwcHdeTIEcfthBkAAAAAAAAAAGClZBVqSGZ2xssvv6y9e/fq8OHDi9cHg0ENDQ1Jks6dO0cLKgAAAAAAAAAAkFdZhxp9fX3y+/06evSohoaGtG/fPknS4cOHFys1/H6/9u7dq9HRUdd1RCIRRSKRxcvhcDiHTYdmZtLfvmqV+erxrPy2FMJKPN5068xyfUttXqlY6mHNzJsHsqrBLOhR+seftHziinN4jpdcZ7bbsMRmlJqsf5WX83ucx32goFbqMWdy/yxWVUqy3veXs9+VyH6/1KaUmqx2uZXYd6vpuCKfj7Wanrdc8FznJpM/Xhz7ZqSgr39SybwGlsvPR8rj618mK+PYN/v7Z7GqUpLuYdn7kJTZfrfk8iWy3y+1KaWkoP/z5rLO5ay3FHE8Vjg8PwWTdaghmTZT+/fv1969ezU5OSmfz6dAILA4b8Pv9ysQCCgQCLhWaxw/flzHjh1b3pZXs/Z2qbNTmpiQZmdTL1dXJ01Pm6/lbCUebybrzHB9mW5eqUj1sNob29XZ2KmJmQnNLsyqrqZO07PTqqmpcV1PquXrav5uxTk8x0uuM9ttyHAzSk3Gv8rL+T3O4z5QUCv9mNPdP4dVlZJs9/2c9rsS2e+lyvoZOazEvltNxxX5fKzV9Lzlguc6N9n88eLYNyMFef2TSuY1sNx+PlIeX//SrYxj39Tfo4qOfRP3IUlp97slly+R/V4qv59RQf7nXc46c1lvKeJ4rHB4fgrOY1mWlcmCw8PD2r9//+Iw8EAgoN7eXp0/f16StHv37sXbgsGgOjo6FgOPRG6VGhs2bFAoFFJbW9tyH1N1GB+XQqGll/P5pLVrV3xzVtxKPN5M1pnh+jLdvFKR6mGN3x9XaCb2QHyNPq1tSf34l1w+h+c479uQ4WaUmox/lZfze5zHfaCgVvoxp7t/DqsqJQXZ90tkv890U0pNRrvcSuy71XRckc/HWk3PWy54rnOTzR8vjn0zwrFv6cvb61+6lXHsm/39c1hVKXF7WIn7kJR+v1ty+RLZ7zPdlFJSkP95l7POXNZbijgeKxyen7wIh8Nqb29fMifIONQYGRnR8ePHNTg4KEkaGhrSoUOHNDk5qWAwqFdeeWVxxsbQ0JBefvnlxWXztbEAAAAAAAAAAKDyZJoTZFzr0tfXp6efflqnTp2SJL366quLVRo+n0/9/f06ceKEfD6fRkdHMw40AAAAAAAAAAAAMpFxpcZKolIDAAAAAAAAAIDqlWlO4D4BGAAAAAAAAAAAoMQQagAAAAAAAAAAgLJAqAEAAAAAAAAAAMoCoQYAAAAAAAAAACgLhBoAAAAAAAAAAKAsEGoAAAAAAAAAAICyQKgBAAAAAAAAAADKAqEGAAAAAAAAAAAoC4QaAAAAAAAAAACgLBBqAAAAAAAAAACAskCoAQAAAAAAAAAAygKhBgAAAAAAAAAAKAuEGgAAAAAAAAAAoCwQagAAAAAAAAAAgLJAqAEAAAAAAAAAAMoCoQYAAAAAAAAAACgLhBoAAAAAAAAAAKAsEGoAAAAAAAAAAICyQKgBAAAAAAAAAADKAqEGAAAAAAAAAAAoC4QaAAAAAAAAAACgLBBqAAAAAAAAAACAskCoAQAAAAAAAAAAygKhBgAAAAAAAAAAKAuEGgAAAAAAAAAAoCwQagAAAAAAAAAAgLJAqAEAAAAAAAAAAMoCoQYAAAAAAAAAACgLhBoAAAAAAAAAAKAsEGoAAAAAAAAAAICyQKgBAAAAAAAAAADKQl02Cw8PD0uSgsGgzp07p6efflp9fX2SpEAgoKGhIfn9fgUCAR0+fFg+ny/vGwwAAAAAAAAAAKpTVqHG/v379dprr2nPnj2amJjQ/v37NTo6unjb+fPnJZmA49ChQxocHMz/FgMAAAAAAAAAgKqUVfupwcHBxcoMSYuVGIFAwLGc3+9frOoAAAAAAAAAAADIh6wqNfbs2bN4fnBwUEeOHJFk2lJ1dnY6lu3s7NTIyIgjBAEAAAAAAAAAoCzNz0uzs9LMjBSJOE/ZXvfII9KhQ8V+RGUpq1BDkkZGRvTyyy9r7969Onz4sCQzY8PNxMTEsjYOAAAAAAAAALDCzpyRamqkZ55Jvm1gQIpGpQMHCr1VmbGsWFiQGDjk+7qFhfxt9717+VtXlck61Ojr65Pf79fRo0c1NDSkffv2pVw2VdgRiUQUiUQWL4fD4Ww3AwAAAAAAAACQDzU10unT5nx8sDEwYK4/eDC79VmWNDe3dMVCPgKHubn8PQ+Z8nikhgapsVHyep2nTK9bu7bw210hsg41JDNLY//+/dq7d68mJyfl8/mSqjImJiYWZ24kOn78uI4dO5bLtwYAAAAAAAAA5JMdZHz+89L770sf/aj0p38qvfaatGuXND0t/eZvZt5eaXbWBBuFVl9vAoTEwCHf19XXm2ADReGxrMx+u4aHh7V//35NTk5KMsPBe3t7df78+cWQ4/z584vLd3R06N1333UNNtwqNTZs2KBQKKS2trZlPiQAAAAAAAAAQEqWJY2PS6Oj0pUr5jQ6Ko2MSDdvmjfsLUvq7jan5aitza6CIZtl469raDAVJyhb4XBY7e3tS+YEGVdqdHZ2OgaFj4yMyOfzuQ4CDwQC6u/vT1mp4fV65fV6M/3WAAAAAAAAAIBczM5KV686A4xAwH2mQ3e3NDERq0Y4eDC7CobEsKGhQarLqVkQkFLGv1F9fX16+umnderUKUnSq6++6qjMGBwc1NGjR7Vz506dO3dOg4OD+d9aAAAAAAAAAIC7ycnk6ovr182g70S1tdKmTdKWLebU2yu9+ab08ssmiJifNyGH2/BwoIgybj+1kjItKwEAAAAAAACAqrewIL33XnKA8XejA5K0tZnQwg4vtmwxgUZ8FUX8UPBnnkm+DKywvLefAgAAAAAAAAAU2L17JrCwA4zRUendd6W5ueRlPR7poYec1Re9vdKaNekHW7sFGPbX06edl4EiI9QAAAAAAAAAgGKzLGlsLLn64tYt9+WbmiS/3xlgbN5s5llkKxp1r8iwL7u1rwKKhPZTAAAAAAAAAFBIkYgZ1h1ffREISNPT7st3dSVXX6xfn776AigztJ8CAAAAAAAAgGKyLOnu3eTqi/ffN7clqq+XHnnEGWD4/VJra8E3HShVhBoAAAAAAAAAsFzz89K1a87qiytXpHDYfXmfL7n6YsMG5/BuAEnYQwAAAAAAAAAgG+FwcvXFtWsm2EhUU2PCisQAo7Oz8NsNVABCDQAAAAAAAABwE41KH3yQXH1x54778i0tJrCIDy8eeUTyegu62UAlI9QAAAAAAAAAgOlp6d13ndUXgYAZ6u2mpye5+mLdOoZ3AyuMUAMAAAAAAABA9bAsaXw8uX3UjRvuy3u90ubNsfBiyxYzvLu5ubDbDUASoQYAAAAAAACASjU7a2ZdxIcXo6PSvXvuy69enVx98fDDZi4GgJJAqAEAAAAAAACg/E1OJldfXL9u5mIkqq2VNm1yVl/09krt7YXfbgBZIdQAAAAAAAAAUBxnzpgqiGeeSb5tYMAEEgcOOK9fWJDef98ZXly5YkINN62tydUXmzZJ9fX5fjQACoBQAwAAAAAAAEBx1NRIp0+b8/HBxsCAuf5f/AvprbecAcbVq6atVCKPR3rooeTqizVrGN4NVBBCDQAAAAAAAADFYQcZp09LwaD00Y9Kv/M70quvSj090u//vjklamx0Bhe9vWZ4d2NjQTcfQOERagAAAAAAAAAovPl56fvfN0O7Z2elz37WVFRYltTdLTU1meW6upKrL9avp/oCqFKEGgAAAAAAAAAKIxyW/vqvpddfl958U7p/31zf0GBaUXm9UkuL9Ku/GqvAaG0t7jYDKCmEGgAAAAAAAABWhmVJ169LZ8+aIOMHPzDX2drbpSefNGHHwoIJNebnTdjxxBNF22wApYtQAwAAAAAAAED+zM1Jf/M3JsR4/XVpbMx5u98vPfWUtGuXtHWr9Lu/K/3Zn0mf+YyZsWEPCZecw8MBQIQaAAAAAAAAAJYrGJTeeMOEGOfOSQ8exG6rr5f6+kyQ8eST0rp1sdvsAOPgwViAET88PP4yAIhQAwAAAAAAAEC2LEsKBGLVGBcvOttKdXaaEOOpp0ygYQ/9ThSNOgMNm305Gl2Z7QdQtjyWFf/XpjjC4bDa29sVCoXU1tZW7M0BAAAAAAAAkGh2Vvrud2NBxvi48/ZHHzUtpZ56SnrsMcnjKc52AihLmeYEVGoAAAAAAAAAcHf3rmkrdfasdP68FInEbvN6pR07Ym2l1qwp3nYCqBqEGgAAAAAAAAAMy5LeecdUYpw9K12+7Lx97dpYW6mPftQEGwBQQIQaAAAAAAAAQDWbmZFGRmJtpe7ejd3m8Uhbt8aCjN5e2koBKCpCDQAAAAAAAKDajI/HQozvftfMy7A1NUn9/bG2Uh0dxdtOAEhAqAEAAAAAAABUOsuSLl0yLaVef10aHXXe3t0dG/L9kY9I9fXF2U4AWAKhBgAAAAAAAFCJpqfNcO+zZ82w72AwdpvHI33oQybE2LVL2rSJtlIAygKhBgAAAAAAAFApxsZibaW+9z1pfj52W3Oz9LGPmSDj4x+X2tuLtpkAkCtCDQAAAAAAAKBcRaPS3/6tCTHOnpWuXXPe/tBDsbZSH/6wVMfbgQDKG3/FAAAAAAAAgHJy/7705psmyHjjDWlqKnZbTY0JL+wgY8OG4m0nAKwAQg0AAAAAAACgWM6cMUHEM88k3zYwYCoxDhyQ3n8/1lbqrbekhYXYcq2tpq3Url3Szp3mMgBUqKxCjZGREQ0PD0uSzp07pxdffFE+n2/xNknq6+tTIBBQMBhUX19ffrcWAAAAAAAAqCQ1NdLp0+Z8fLBx+rT0f/6PGeb92msm1Ii3aZOpxHjqKbNMbW3hthkAiiirUGN4eFjPPfecJOnEiRPavXu3zp8/L0k6efKkTp06JUnas2ePBgcH87ypAAAAAAAAQIWxg4zTp6V796THHpN++7elb35TWrtWun7d3F5XJ33kI7EgY/364m0zABSRx7IsK5MFR0ZGtHv3bk1OTkqSAoGAent7NTo6Kr/fr1OnTulTn/qUJC1Wb2QqHA6rvb1doVBIbW1t2T0CAAAAAAAAoNzMzEiXL0sXLkgXL0pf+5r0zjuSxyNZltTdLT3+uPTkkybE2LlTam4u9lYDwIrJNCfIuFKjr69PL7744uLlYDAoSers7Fy8LtswAwAAAAAAAKh40ah07ZoJLy5eNEHGu++a8MLW0mICDa9X6uiQBgelbdtMeyoAwKKs2k/t27dv8fzLL7+sPXv2LAYZwWBQQ0NDksy8jSNHjsjv97uuJxKJKBKJLF4Oh8PZbjcAAAAAAABQmu7ejVVgXLwovf229OBB8nJr1pjgYvt26coVM/zb65Xm56Xz582sDACAQ1ahhs0OMOx5GpJ0+PDhxYDD7/dr7969Gh0ddb3/8ePHdezYsVy+NQAAAAAAAFA6HjwwbaTsCoyLF6U7d5KXa2oy7aS2bYud1qwxtw0MmGHgn/mMmbExMOA+PBwAkPlMjXhHjhzR0aNHHZUYIyMj6uvrk2RCj46OjsV5G4ncKjU2bNjATA0AAAAAAACUrmhUuno1VoFx8WJyGynJtJHavNlUYNgBxqZN7q2k7ADj4EFngJHqegCoUHmfqWE7ceLEYqBhz9UIBAKOIeK2+Hkb8bxer7xeb7bfGgAAAAAAACicO3ecFRhvv20GfCdauzYWXmzfLj36qKnMyEQ06h5c2Jej0eU9BgCoMFmFGkNDQ+rr61sMNF555RUdPnxYfr9fL7zwwuJyw8PD2rdvH4PDAQAAAAAAUB4ePDChRXwVRqo2Ulu3OttIrV6d+/c9cCD1bVRoAECSjNtPBQIB9fb2Oq7z+XyL1RkjIyMaHh6Wz+fT6OioI+RYSqZlJQAAAAAAAMCy2W2kLlyQLl1K30bK73cGGKnaSAEAliXTnCCnmRr5RqgBAAAAAACAFXP7trMCI1Ubqa4uE1xs3WraSD32mNTYWPjtBYAqtGIzNQAAAAAAAICSNT1tQgu7AuPCBenu3eTlmpulxx/PXxspAEBBEGoAAAAAAACgPC0smDZS8VUYV68mt5GqqTFtpOwKjG3bpI0baSMFAGWIUAMAAAAAAAClz7JibaQuXTIVGG+/LUUiyct2dZnwwh7oTRspAKgYhBoAAAAAAAAoPXYbqfgqjFRtpOzwwp6HQRspAKhYhBoAAAAAAAAorvg2UhcumK/Xrrm3kertjYUY27ebNlIeT1E2GwBQeIQaAAAAAAAAKJz4NlL2aak2UnYVxmOPSV5v4bcZAFAyCDUAAAAAAACwcuw2UnYFxsWL0sRE8nLxbaTseRidnYXfXgBASSPUAAAAAAAAQH4sLEjvvhsLLy5ckK5fT91Gyq7A2LaNNlIAgIwQagAAAAAAACB7dhup+AqMy5fd20itWxerwNi2TXr0UdpIAQByQqgBAAAAAABQjc6cMRUTzzyTfNvAgBSNSgcOxK6bnpYuXYpVYFy6lLqNVHwFxrZtUkfHSj0KAECVIdQAAAAAAACoRjU10unT5nx8sDEwIL30kvRjPyZ9+csmvEjVRqq2VvL7YxUYW7fSRgoAsKIINQAAAAAAAKrRT/2UdP++9Nu/Lb33nrRrl6neeO01M6D7a18zp3jd3c4KDNpIAQAKjFADAAAAAACgnFmWdO+eFA5nd5qeNvcPBqX/9t9MdYVlmeBi9WqppSUWXmzdShspAEBJINQAAAAAAAAoBZYlzcxkH05MTSW3hcqUx2OqLYJB045q1SrpN3/TBBgbNtBGCgBQcgg1AAAAAAAA8m1uzgQOoZAJHdKdjz/Nz+f+PZuapLa27E6rVkm/+7tmW+rqzPe/dUv6R/8of88FAAB5RKgBAAAAAACQysJCLHyYmsospAiFpEgk9+9ZX599ONHaau6XrYEBMyz84EEzLNy+LDmHhwMAUCIINQAAAAAAQGk4c8a0QHJ7M31gQIpGpQMHclu3ZZmh2HYwER9SpDt//37uj8fjyT6caGszg7cL0fYpMdCQYl8JNgAAJYpQAwAAAAAAlIaaGvc30+PffLcsUwWR7cyJcNiEIrlqack+nGhpKe2ZFNGoM9Cw2ZeX83wBALBCPJaV6ySp/AmHw2pvb1coFFJbW1uxNwcAAAAAABTLb/2W9OKL0o4d0tat0re/LX3nO5LfL3V1mXBibi739Xu9ubV2qq3N32MEAABJMs0JqNQAAAAAAADFMTMjvf22dPFi7HT7trn+D/7AVDlYltTdbUKFu3dj962tjYUO7e0meEh33j41NBTv8QIAgGUj1AAAAAAAACsvGpWuXYuFFxcuSO++a0KLeB6PtGuX9K1vSXV1JoT4H/8jOaRoairt1k4AAGBFEGoAAAAAAID8u3vXBBd2iPH229KDB8nLrVkjbdsmbd9u2k09/rg0OCgFAibUmJ+X3n+fgdUAAEASoQYAAAAAAFiuBw+ky5djFRgXL0p37iQv19hogoutW02IsW2bCTXixQ8Ff+aZ2GWJYAMAABBqAAAAAACALESj0tWrzjkYqdpIbd5sggv79MgjUk1N6nUnBhpS7CvBBgAAEKEGAAAAAABI584dZ4CRSRupbdukxx4zcy+yEY06Aw2bfTkaze0xAACAiuGxrMSPUhReOBxWe3u7QqGQ2trair05AAAAAABUp0zbSDU1mdkX8VUYiW2kAAAAspBpTkClBgAAAAAA1SjbNlJ2Bca2bdKmTenbSAEAAKwQQg0AAAAAAKqB3UbKrsB4+21pZiZ5ubVrnRUYubSRAgAAWCGEGgAAAAAAVJoHD0xoEV+FkaqN1Nat5mTPw1i9uvDbCwAAkCFCDQAAAAAAyhltpAAAQBUh1AAAAAAAoJzk0kZq+3bTRqqxsfDbCwAAkEdZhRojIyMaHh6WJJ07d04vvviifD6fJCkQCGhoaEh+v1+BQECHDx9evA0AAAAAAOQgsY3UhQvS3bvJy9ltpOJnYdBGCgAAVKCsQo3h4WE999xzkqQTJ05o9+7dOn/+vCRp//79i+cDgYAOHTqkwcHBPG8uAAAAAAAVym4jZVdgXLqUuo2U3x+rwNi6lTZSAACgangsK/HoyN3IyIh2796tyclJSSa46O3t1ejoqCRnqCFJHR0di8suJRwOq729XaFQSG1tbdk+BgAAAAAAys/t2845GKnaSHV1OSswaCMFAAAqUKY5QcaVGn19fXrxxRcXLweDQUlSZ2enXnnlFXV2djqW7+zs1MjIiPr6+rLcdAAAAAAAKsz0tHT58tJtpJqbpccfp40UAABAClm1n9q3b9/i+Zdffll79uyRz+dbDDgSTUxMuF4fiUQUiUQWL4fD4Ww2AwAAAACA0pXYRuriRXM5sVFCTY20eXOsjdS2bdLGjbSRAgAASCOrUMMWDAY1NDTkaDeVajk3x48f17Fjx3L51gAAAAAA5N+ZMyZMeOaZ5NsGBkxQceCA+33j20hduGAqMmgjBQAAsCJyCjWOHj2qV199VT6fT5Lk8/mSqjImJiYWb0/0/PPP69lnn128HA6HtWHDhlw2BQAAAACA5aupkU6fNufjg42BAXP9wYPm8vS0mX1hD/Jeqo2UXYGxdSttpAAAAPIg61DjxIkTOnr0qPx+/2Ilxp49e3Ty5MmkZfv7+13X4fV65fV6s/3WAAAAAACsDDvIiA82zpyRPvc5adcu6dYt6ed+LnUbKb/fWYVBGykAAIAVkVWoMTQ0pL6+vsVA45VXXtHhw4eTKjICgYD6+/tTVmoAAAAAAFByfvqnpQ8+kH7t16Rf+RXp3j1p3Trpu981J1tXl7MCgzZSAAAABeOxrMSPmLgLBALq7e11XOfz+TQ5Obl4+8mTJ7Vz506dO3dOzz//fMahRjgcVnt7u0KhkNra2rJ7BAAAAAAA5Or+fencOen116W//mspFJL+5m9MNYbHIz31lAku4qswOjuLvdUAAAAVJ9OcIONQYyURagAAAAAACub9902I8frr0ltvSQsLsdsmJ02w4fOZ6ot/9a+kn/3Zom0qAABAtcg0J8hpUDgAAAAAAGVjft6EF3aQ8cEHzts3bjQVGbdvS8PD0rPPmpka9pBwj8c5PBwAAABFQ6gBAAAAAKg8k5OmndTrr0vf+Y40PR27ra5OeuIJE2Q8+aS0fr0JMP78z6VPfzoWYLgNDwcAAEBREWoAAAAAAMqfZUlXrsSqMS5dct7e0WECjKeeknbskJqbnbdHo9LBg8nBhX05Gl25bQcAAEDGmKkBAAAAAChPMzPS+fMmxHjjDenuXeftjz1mQoynnjLnPZ7ibCcAAACWxEwNAAAAAEDlGRszAcbrr0vf+540Nxe7rbFR6u83IcbHPy6tXl20zQQAAMDKINQAAAAAAJSuhQXpb/821lbq2jXn7T09sWqMj3xEqq8vznYCAACgIAg1AAAAAAClZWpKevNNE2K8+aa5bKupkT784diQ740baSsFAABQRQg1AAAAAADFZVnS1auxtlI/+IG5ztbaGhvy3d9vLgMAAKAqEWoAAAAAAApvdtbMxLDbSt265bzd74+1ldq2zVRoAAAAoOoRagAAAAAACuPOnVg1xvnzUiQSu62hQfroR2NBRldX8bYTAAAAJYtQAwAAAACwMixLunQpFmS8847z9jVrYiFGX5/k9RZnOwEAAFaIZUkTE9KNG7HTBx+YsWDPPFPsrStPhBoAAAAAgPyZnpbOnTNBxhtvSMFg7DaPx7SSsoMMv58h3wAAoOzNz0s3b5qw4sYNaWzMeT6+ONX2oQ8RauSKUAMAAAAAsDzvv28qMd54Q3rrLfOfva25WfrYx0yI8bGPST5f0TYTAAAgV/fuOcOK+NP4uKnISMXjkbq7pfXrpZ4e6aGHpM2bC7ftlYZQAwAAAAAgnTljhnG7fWRwYECKRqUDB8zl+Xnp+9+PDfl+/33n8hs2SE8+aYKMD39YquNfTwAAUNosy4z/im8RNTYWOz81lf7+Xq8JK9avTz6tW8fhUD7xVAIAAAAATKBx+rQ5Hx9sDAyY6z/1KenP/syEGOfOmTZTtro66SMfMUHGk09KDz9c2G0HAADIwOysCSriwwq7RdSNG9LcXPr7d3QkBxZ2kOHz0VWzUAg1AAAAAACxIMMONn7mZ6Rf/3Vzef16aXDQ2VfB5zOVGE8+KfX3mzZTAAAARTY15d4i6sYNU4mRrk1Uba2pqrDDip4e5/mmpsI9DqRGqAEAAAAAMJUXfX3SyIh07Jj0H/6D+Thjd7dpPSVJjz4aG/L9+ON8HBEAABRcNCrdvp0cWNjtou7dS3//pqbkKgv71NVlgg2UNkINAAAAAKgmliXdvClduSKNjsZOY2OxZSYnzXK1tdI//+exiozVq4u33QAAoGpEIu4toj74wBzGzM+nv//q1e6zLdavl9rb+VxGuSPUAAAAAIBKNTMjBQLmZIcYgYBzHka8NWukBw+ktWvNf/xNTdJjj0n/5J8UdrsBAEBFsywpFHKGFfFVF3fvpr9/XZ0pJk2cbWG3i/J6C/M4UByEGgAAAABQ7izL9GFIrL744AP3xtF1ddLmzVJvr/P0x39sZmj88i+bGRv2kHDJOTwcAABUlDNnpJoa95f7gQHT8unAgezWubAgjY+7z7a4cSP1ZyxsLS2ph3KvXWu2F9WJUAMAAAAAykkkIl296gwvAgEzFdNNR4cJLLZsiYUXGzaYYCOeHWAcPBh7RyNxeDjBBgAAFammxv3lPv7wwM3MjPtsixs3pFu3TLCRzpo1zrCipyd2vrWVNlFwR6gBAAAAAKXIsqSJieTqi+vX3asvamulTZuSqy86OjL7ftGoM9Cw2ZftYeEAAKDiuH2OYWBAeukl6emnpR07pG98I7ld1ORk+vXW1SWHFfapu5s2UciNx7LcjoYLKxwOq729XaFQSG1tbcXeHAAAAAAorPn55OqL0VHTbNpNW1ty9cWmTVJ9fUE3GwAAlKdIRAoGzaFGKBQ7/7WvSd/8pqmwmJkxbZ5Wr06/rtbW1EO5166l2gKZyzQnoFIDAAAAAAopGEyuvrh2zb0/g8cjbdzorLzYskXq7OQdAgAAIMkcQkxNOcOJdOdDIRNqpDIxYYpCPR4TaHg8Jpxwm23R02NCDaCQCDUAAAAAYCXMz0vvvRcLLq5cMbMvJibcl1+1Krl11COP0JcBAIAqYlnSgwfOAMItnIi/7t49986US6mrk9rbzcnnM18vX5Zu3pQaGsztP/7j0s//PMWgKC2EGgAAAACwXFNTydUX775rgo1EHo/5eGNi9QX9GQAAqDjz886AIl04YZ93O3zIRGtrckiR7nxTk/PQY2BA+ou/kJ5/PjZT4/RpU62ROHILKCZCDQAAAADIVDQqvf9+cvXF7dvuyzc3S35/LLjw+82psbGw2w0AAJbNsqT795NbOaVr9XT/fm7fy+uNhQ+ZhBStrVJtbe6PzQ4wDh6MBRhuw8OBUkCoAQAAAABu7t9PHtz97rupm1D39CRXX3R3U30BAEAaZ85INTXub5gPDJjPExw4sDLfe3bWfd5EqqAiHHYfgbUUj8cZUKQLJ+yvhe4+GY06Aw2bfTkaLez2AOkQagAAAACobpYl3bjhrL4YHZVu3XJf3utNrr7o7TVVGQAAICs1Ne6VAPGVA5mwLBM6pBqM7db2aWYmt21ubs48pLCrKEr9Mw7pgiMqNFBqCDUAAAAAVI/paVNtEV99EQiYiZxuurqSqy/WrzfvwAAAgGVza3H0O78jfeEL0v790sc/Lr355tIhRTic27Ds2lpnALFUSNHWZoZoAygej2XlsrvnVzgcVnt7u0KhkNra2oq9OQAAAADKnWWZSovE6osbN9yXr6+XNm+OBRe9vaYCo7W1sNsNAEAFsts82af4igr78tmz0ve+Z9oczc5K69aZLo7ZWrUqOZBIF1I0N5d+FQVQLTLNCbKq1BgZGdGhQ4d0/vz5pOslqa+vT4FAQMFgUH19fTlsNgAAAABkKRIx1Rf20G47yEg1mXP16uTqi4cfXt50TQAAqsTsrDOUSAwoEk/hcOZtnmZmzOcSPB4TaNTXJ8+aSBdUtLVJdfSlASpexrv50NCQ/H7/YoAR7+TJkzp16pQkac+ePRocHMzfFgIAAACoXNlMB7Us6c6d5OqL99937zdRVydt2pRcfeHzreADAgAgO8UclD03lz6UcKuoSNWxcSm1tcmBRFtb7Pybb5r12wOy/+W/lD79aaooACTLONTYt29fytt27NihyclJSZKPfxAAAAAAZCrVdNCXXpJOnpQ+8Qnpf/9vU4Fx5Yo0NeW+Hp8vufpi40Y+rgkAKHn5GpQ9P790IJF4eXo6921OFU6kOjU1pQ4oBgakc+ekf/fvzHNgP/aGBoZUA0iWtyN8wgwAAAAADpZlelTMzJgWUZGIOR9/3YYN0o4d0q//ummm/fDD0p/9mfSDH5i+EyMj5mSrqTFhRXz1RW+v1NlZvMcJAMAyuA3KPn3aDMr+5CelH/oh6VvfSh9QhELLCyjShRKJt7W1SS0t+augiA9v7OfC7TkBAFteQo1gMKihoSFJ0rlz53TkyBH5/f6Uy0ciEUUikcXL4XA4H5sBAAAA5Fcx+0GsFDtoSAwb4k+JwUOq6+Jvc1vP7Gzm2zU/L738snmHxLJMoPHoo862UVu2mHZSDQ0r9/wAAFAA9+9LY2PmdPOmFAyaoODoUekXf9G8LHZ3S1/9qjllyuNJXTmR6rpVq4rb4ikadQYaNvtyNFr4bQJQ2vISahw+fHixUsPv92vv3r0aHR1Nufzx48d17NixfHxrAAAAYOXkqx/EUizLvHuxVFCQS7CQeJ/ZWff5Eyutrk5qbDSNsr1ec76hwXndF79orlu1SvqTP5HWrKGRNgCgLEUiJqy4eTMWXtgBxs2bqbspRiKxQdk9PVJra2bhhH0qdkCRi3SfD6FCA4Abj2Vl9x+Nx+NR4l1GRkbU19cnyVRtdHR0aHR0NGW1hlulxoYNGxQKhdTW1pbtYwAAAABWjh1g/PRPS3v2SL//+9If/IH0Yz9mLucSLLhdV4ygobbWGTLY59MFD27XpVuHfaqpSb8t9vNcV2cCHrePbAIAUCLm56Xbt90Di7ExaWJi6XW0t5vgorvbfL140bSZam42t3/60+VXEAoAyxEOh9Xe3r5kTrDsSo2RkRHt3r17cVC4rTNNT1uv1yuv17vcbw0AAADkn2VJd+5Io6PmdPWqNDcn/dIvOVsjfetb5pRvHk92gUGuYURDQ+kM0U5spm1flgg2AABFYVnS3bvOwCL+6/j40p9HaG6OBRb2yb7c3W0GZ9sGBqTvfS95UHaqLpgAUM1y+i8mGAw62k298MILi7cNDw9r3759DA4HAABA6ZudNaGFHWDYp8SeEPX1JmzweEwY0N+//OqFVMFDbW359Y1YDqaDAgCKwLLMsG23wGJsTLp1y3ymIZ36emdokXi+tTWzl3ReCgEgOxmHGsPDw3r11VclmZkYO3fuXAwv+vv7deLECfl8Po2OjmpwcHDFNhgAAADImv1xy8Tw4r333KdP1tRIGzaYwdR+v/T22yYAaWqSFhakvXt5dyFfmA4KAFgh09POoCIxvHjwIP39a2qkri73wKKnR+rszM/nEHgpBIDsZD1TYyVk2isLAAAAWNLsrHTtmnTlihQIxAKMcNh9+dZWE1709poAY8sWadMmU0UhpW6NxMwHAACKanbWVFSkqrZI9dIfb/Vq9xZR3d0m0KitXfnHAQAwCjZTAwAAACgKyzJTOOPDiytXUldfeDzSxo3O8KK317ybkepjlvSDAACgaBYWYsO444dw26e7d5deR2ure2DR0yOtW2e6PwIAyguhBgAAAErf3Fxs9kUgEAsyQiH35VtbTWARf9q0Kft3LugHAQDAirE/n+AWWNy8aYZxLyykX0djo3tgYZ+amwvzWAAAhUP7KQAAAJQO+92NxNkX16+nrr7YsCE5wFizprqGbQMAUIIsS7p3zz2wsL/OzqZfR12dM6xInG3R3s5LPgBUCtpPAQAAoLTNzZnZF4kBRqrqi1WrYqGFPcD7kUfoGwEAQI7OnDHDsN06KQ4MmM8THDiQfh0zM+7DuO3T9HT6+3s80tq1yUO47a98TgEAkIhQAwAAACvLsqTJSffqC7eeEh6P9PDDzvCit9e848G7GgAA5E1NjfuIqPiRUvPzycO44wOMYHDp79PR4R5Y9PSYYdx1vDsFAMgCLxsAAADIn/l59+qLVO94rFrlHNrt90ubN1N9AQBAATzzjPnswcmT5rMGO3dKQ0PSa6+Zl+avftVUcyzVuHzVqtQzLbq7eVkHAOQXoQYAAAByY1df2EO7R0dNoLFU9UX8ieoLAABW3P37seqKGzec1RZjY2ac1X//7+Yl2bJMEDE3J92+be7v9aYOLHp6TKgBAEChEGoAAAAgvfl58/HNK1dMcBEImPOpqi9aWpLDC6ovAABYMfPz0vh4cmBhX56aSn//9evNZxXq66WmJulXfsUZYPh8fAYBAFA6CDUAAAAQEz/7wg4vrl8375Yk8nikhx5KDjC6unjnAwCAPLIsKRRyhhbx58fHl24R5fM5qyx6ekyY0dMjff3r0u/8jpltMT9vii537y7IQwMAIGuEGgAAANXIrr5InH0xOem+fHOzc2j3li3SI49IjY0F3WwAACpVJBILK27eTA4wIpH0929ocIYV3d2x0KKnx1RguBkYMIHGwYNmxoY9JFxyDg8HAKBUEGoAAABUumAwOby4dm3p6ov4Ad5UXwAAsCzRqHTnTnJYYQcYqT5XYPN4zCiqxLDCDjA6O7N/qbYDDDvQkGJfCTYAAKWKUAMAAKDUnDkj1dS4v4swMGDeFTlwIPk2u/rCHtptnyYm3L9Pc7P77AuqLwAAyMnUVOpqi1u33D9PEK+lJRZYJFZbrFtnZl7kUzTqDDRs9uVoNL/fDwCAfCDUAAAAKDU1Ne4fj4z/OGUwGAsvrlwx569eTV19sX59coCxbh3VFwAAZGFuzoQTiQO57SDj/v3096+rMy+/8fMs4k+trYV5HDa3z0jYqNAAAJQqQg0AAIBSE9/3YWFB+pEfkT73OemP/9i0g/ryl2OhR6Lm5tjci/jqi1SNtAEAwCLLMgWOiWGFff7OnaUHcnd0JAcW9uU1a8xnFwAAQO4INQAAAErN3bumMXZrq/SLv2iusyzTg2J62pykWPVF/ADv7m6qLwAASGN62rSGSgwsxsbM9UsN5G5sdA8senrMyzBdHAEAWFmEGgAAAMVmWaaF1Nmz5nT5svN2j8c02T5yJDbA2++n+gIAUHFyHSsVb2FBGh9Pnmlhnw+F0t/f45G6upJnWtgBRns7nx8AAKCYCDUAAACKIRKRRkZMiPH666Y6w+bxSFu3mndlpqdNxcbCgqne+Kf/tHjbDADACstkrJRlSeGwM6yIDzBu3Vp6wHVra3JYYQcYXV1m9gUAAChNvEwDAAAUyt27JsB4/XXp/Hlnf4vGRqm/X9q1S3rySekrXzHv3vybf2Pe1bHfzZGY3AkAqFj2S9wXvmDmV/zwD0tf/KL09a9L27ZJ3/ym9PLLsU6MqdTVJYcV8S2jWlpW/rEAAICVQagBAACwUixLeucdE2K4tZXq6pKeesoEGU88ITU0mOvjP45qv7sTPzw8/jIAAGVofl66fTs2xyJ+psXYmGkf9V//qyletMdKhcPmZFu9OnW1xerVtIgCAKBSEWoAAADkUyZtpXbtMqfNm93fcYlGnYGGzb68VE8NAACKzLKkiYnkIdz2+du307+cdXebZerrzQipf/tvnQHGunWS11u4xwMAAEoHoQYAAMBy2W2lzp41gUZiW6mdO01FxpNPSh0dS68v3QRUKjQAACXAsqSpKffAwp5rMTeXfh319Sa8sIOK7u7Y5b/8S/M96utNVceqVdJP/EQhHhkAACh1hBoAAADZsttKnT1rTu+847w9VVspAADKyIMHqdtD3by59FyLmhrzkmgHFfEBRk+P1NnpXrA4MCB96UvSz/0cY6UAAEAyQg0AAIBMRCJmuLc96DuXtlIAAJSQuTlTUeEWWIyNSaHQ0uvo7HQPLLq7pbVrzcDubDBWCgAALIVQAwAAIJU7d5xtpWZnY7fl0lYKAIACikbNXAq3SouxMZPPW1b6dbS2ugcWKzXXgrFSAABgKR7LWuoQZuWFw2G1t7crFAqpra2t2JsDAACqlWVJly/Hggy3tlJ2NcYTT5hG3wAAFIllScFg6rkW4+PSwkL6dXi9yYFFfIjR0lKQhwIAAJBxTkClBgAAqG5LtZXats2EGE89RVspAEDB3bvnHljcvGlOkUj6+9fVmYqKVHMt2tt5aQMAAOWFUAMAAFSfdG2lmpqk/n4TZDz5pOTzFW0zAQCVLxJJHViMjZlQIx2PR1qzJvVcizVrzMBuAACASkGoAQAAKp/dVursWRNmJLaVWrcuVo1BWykAqEpnzpg3/92GUA8MmFkOBw5kv975edMGym0Q99iYNDm59Dp8vtRzLbq6eNkCAADVhVADAABUJrut1Nmz0htvJLeV2r7dhBi0lQIAyAQap0+b8/HBxsCAuf7gQff7WZZ5iUk11+L27aWHcTc3p55r0d1tiggBAABgEGoAAIDKcft2bDYGbaUAAFmwg4z4YGNgQHrpJeknf1L62Mekv/xL9zZR8/Pp193Q4AwpEqsuVq0iWwcAAMiUx7KW+szIyst0qjkAAKgimfQB+dmfld5+OxZkpGortWuX9JGP0J8DACDJVE7MzEihkDmFw7HzX/uaCS8WFswyXV1SZ2f69dXUmOVSzbXo7CS0AAAAWEqmOUFWlRojIyM6dOiQzp8/77g+EAhoaGhIfr9fgUBAhw8flo9PPwIAgOVI1Qfk85+Xfuu3pMcek77yFWliInZbfFupXbukRx7hXSQAqHDpAor48/GncFiam0u9zokJs16PJxZorF6deq7F2rVSbW1hHi8AAEC1yzjUsEOLkZGRpNv279+/GHQEAgEdOnRIg4OD+dtKAABQfeL7gITD0saN0qlTZkbGunXSrVvm9qYmaedOE2J8/OO0lQKAHK3UoOxsrERAkU59vdTebk5tbebr5cumtZTXK9XVST/+49K//tfmMgAAAIov41Bj3759rtcHAgHHZb/fr+Hh4eVtFQAAqF7T06al1MWLpp3UvXvSf/7P5uOylmU+FvvEEybEeOop2koBQJ7kOig7lfiAIjGQWOmAIj6kSDzFX+/1Ogv6BgZM66n/9J9iMzVOnzZVGm5hDwAAAApv2YPCh4eH1ZnQYLSzs1MjIyPq6+tb7uoBAEAlW1iQAgHp0iUTYly8KF27Zt4Js9nTU5uapI4O03KKtlIAkHepBmWfPm0qNPbvN0OxixFQZBJSJAYU2YoPb+znwu05AQAAQHEtO9QIBoOu10/E97dOEIlEFIlEFi+Hw+HlbgYAACh1liWNj5vg4sIFE2RcvizFHRMsWrdO2rbNnN55x4QfXq80Py/91V9JmzcXfvsBoILMzZnQIRg0AYT9dWFBWr9e+uxnTZHc7Kz08MPSF79o2lPlor7edAaMDyVWOqDIRTTqDDRs9uVotLDbAwAAAHfLDjVSSRV2SNLx48d17NixlfrWAACgFNy7Z4KL+CqMycnk5VpapK1bYyHG1q2xqawDA9LwsPSZzzg/MizxcVkA+DuWZf7kxldJ2CFF/Pn466an068zFIoNyl61KlZx0dDgHkykCymKEVDkIt28EF5yAAAASseyQw2fz5dUlTExMSFfmiGdzz//vJ599tnFy+FwWBs2bFjupgAAgGKZn5dGR02AceGCCTDeey95udpaqbc3FmBs2yZt2OD+bhd9QABUqbk591AiVWARDpsKi2zV1JjgweeLhRA+n/SDH0i3b8cGY//jf2ze8C+ngAIAAACVa9mhxp49e3Ty5Mmk6/v7+1Pex+v1ymsfIQMAgPJiWdLYWKz6wh7o7dY4vadH2r49VomxZUvsXbKl0AcEQAWIr6Jwq5hwCywePMjtezU3O8MJt6/xJ3tkUbyBAZNRP/ecs0Cuu5scGQAAAKUhp1AjGAwuVmL4/X7HbYFAQP39/WkrNQAAQBmZmnIGGJcumXfeErW2xtpHbd8uPf64eRctV/QBAVCCZmeTg4h0gUU4nFsGW1ubOoxwCyza2szsiuWgQA4AAADlIONQY3h4WK+++qokMxNj586d2rdvnyRpcHBQR48e1c6dO3Xu3DkNDg6uzNYCAICVNTcnXbniDDE++CB5ubo6U3WxfXssyHjoIXqSAMjamTOmDZLbm+UDAyYQSJdxLodlmdw2VSjhFljMzOT2vZqbM6uesK9raWFQNgAAAODGY1mWVeyNCIfDam9vVygUUltbW7E3BwCA6mBZJrCIDzCuXDHzMRI9/LCzCqO3d/kfCQYAuVcHpLs+nUgkffWEW2CRy39DdXWZVU/E31637Ma/AAAAQGXLNCfg0BoAgGoRDJrWUfFtpKamkpdrb3cO8t661bSWAoAV4NbeaGBAeukl6Sd/UvqRH5Heemvpdk/BoAk1ctHSkln1hP21uZnCNAAAAKBYqNQAAKASRSKm6uLChViQMTaWvFxDg/Too84Qo7ubd+sArLipKVMsduOG+frVr0r/7/+ZFkezs1JXl7RuXfbrja+iyKTdE1UUAAAAQGmgUgMAgGphWdL167HqiwsXpEBAWlhIXnbTJlN5YQcYfj/v5gFYEZZlqic++MAZXtgnt0Kx6WlzP48nFmisWpX5sOz2dqmpiVwWAAAAqGS8iwEAQLmZmIiFF3aQMT2dvFxHhwkutm83QcbWrabHCgDkiWVJt287w4r48GKpodqrV0sPPSStX2+y2MlJyes1ocT+/dLhw+SuAAAAAJz4FwEAgFI2MyNdvuwc5j0+nryc1ys9/rhzmPfatXxcGcCyzc9Lt265hxZjY9LcXOr72hUX69fHwouHHoqdb2w0yw0MSF//uvTss7GZGqdPmyqNTIeEAwAAAKgOhBoAAJSKaFS6ds0ZYLz7rrk+nscjPfJILLzYts1crq0txlYDqACRiAko3NpE3bqV/GcoXl2dGcVjhxXxoUV3t1Rfn/572wHGwYOxAMNteDgAAAAASIQaAAAUz507zgDj7belBw+Sl1uzxjnI+7HHpObmwm8vgLI2PZ26TdSdO+nv6/UmV1nY57u6pJqa3LcrGnUGGjb7crpABQAAAED18ViWZRV7IzKdag4AQNmanjahxaVLsRDD7V3EpqZYGyn7tGZN4bcXQNmxLCkcTg4t7K/BYPr7t7QkBxb2qbOTbnYAAAAAVlamOQGVGgAA5NvCgnT1qgkuLlwwQcbVq+Ydx3gej+T3OwOMTZuW95FnABXNsqS7d5NbRNmn6en09/f5kltE2edbWwkuAAAAAJQ+Qg0AAOKdOWNCBbcG7gMDpg/KgQOx6yxLun07Fl7YbaQikeT7d3WZGRhbt8baSNlTcgHg7ywsSOPjqQdzu/15ibd2rftQ7oceonMdAAAAgPJHqAEAQLyaGvfJtPYk25/6Ken8eRNgXLhgQozJyeT1NDfHwgs7yOjsLMxjAFDy5uZMQBEfWNjnx8ZMsJFKTY20bl1yi6iHHjKDub3ewj0OAAAAACg0Qg0AAOLZQcZLL5mh3T/6o9Jv/Zb05S+bjzr/3u+ZU7zaWtNGavt2E2Js3Spt3EgfF6ACZFu8FW9mJnWbqNu3kzvSxaurS11tsW6duR0AAAAAqhH/DgEAKptlmSbz4bA5TU1JoZD5al/ndhofl/7jfzTBhGWZjz/X15t19vQ452Bs2cJHo4EKtVTx1k/+pCncchvMPTGRft1NTe6hxUMPmRZS5KIAAAAAkIxQAwBQHizLNJJ3CyASg4rE0CIazf77rVsn3bxp3tFctUp67rlYFYbPl/eHB6C0zM+bPyWf+IR07Zr0G79hOs995CPSn/yJ9MYbZkzOl75kTqm0tia3iLLDC5+P4AIAAAAAskWoAQAovLk59/AhXVARDpv75crrldrazDuM7e3O862t5nL86U//1IQhDQ3m3c3aWunJJ/P3HAAoGMuS7t83f1JCISkYjJ1PvGyfv3/fuY5IxFRnxBdv2fnm6tXuFRfr15s/LwAAAACA/CHUAADkbn5euncvfZWEW1ARieT+PevqkgOIxHAiMahobc2uPdTAgPTKK9KnP236zdh9ZiT3xvoACsrORe0AIvGrW2CRbvB2Kh6P+XPS3i790A9JX/mKKd5qajKjdh56yHSja2rK7+MDAAAAAKRGqAEApWw5E2qzEY2acCLTdk72aXo69+/p8SSHE+mCCvvU2Liy/VrsAOPgwdjzbn8l2ADyzrJi2Wi6yon463L909PYaP6s+HzOr/Gn+NtWrTJ/giXzp+Gtt0yuOj8vvfuu9Pf/fh6eAAAAAABAVgg1AKCULTWh9uBB5/KJQ7EznTtx7565b65aW91bOKULKlpaSrOZfDTqDDRs9uVc5nMAVWR2NnWbp1QhRS67VXwVhVtI4RZYZFOwFS8x66R4CwAAAACKh1ADAEqZ/W7Z5z9vJtXaE2r/4i+kJ56Q3n5b+vmfX/5QbFtzc+oWTm7hRFub86PMlSBd5QvvXqLK2FUUS1VOxN/24EFu36upKTmIcAsp7K+rVhUmF6V4CwAAAABKC6EGAJSahQUpEJAuXZIuXjSnO3ek3/gN54Ta27fNyU3iUOxMqyjqeFkASsFKdZ6LRDKrnLDP55qT1tRk1+aprU1qaMj++xQCxVsAAAAAUFp49woAismypFu3YuHFxYvSO+8kD9Jet06amDBhRUOD9Au/kL6KItceKwBKQiad5ywr1kku01ZPMzO5bU9z89KVE/GnQlVRFALFWwAAAABQWgg1AKCQpqZMy6gLF2KVGMFg8nLNzdK2bdLWrdL27dL3vicNDsYm1DY1ST/1U4XeegArJH5YdjgsPfaY1N8vvfCC9Jd/KX3oQ9Jf/ZX0ne9ImzdLf/iHppojl1E4tbVLV07Ef21rk+rr8/loAQAAAADIHaEGAKyUuTlpdNRZhfH++8nL1dZKvb0mxLBPGzbEPuY8MGACDSbUAmVhft7kl3ZAkfg11XVuAUU0Kv3RH0l//MexznN1deZ+tpaWzNo82eebmyunigIAAAAAUH0INQAgHyxL+uADE1xcumQqMa5cMe9uJlq/PhZebN0qbdmSul0UE2pRxVZqrkQ24mdQ2OFDfMsnt+ump3P/fk1NsU5ydhDxf/+veR4aGqT/+T+TQwtG4QAAAAAAqgn/BgNALoJB5yDvS5fMR7MTtbaa9lFbt8ZCjPb2zL8PE2pRxTKZK5Epy5Lu33evkkisoIi/PDub27Z7PGauhN2+yf4af97tusQ2TwMDJge1O89dv06OCQAAAACoboQaALCUSMQM744PMMbGkperr5cefdTZRqqnZ3l9XphQiyrmVpQ0MCC99JL09NPSP/gH0ltvpa+giD+fawZYVxcLH+IrKNJd19pqQpnlSCzUovMcAAAAAACEGgDgZFnStWvOAGN01P3d0I0bY4O8t241czHoAwMsy/y8dPu2yQ3HxkylRGen9Pzz0nPPmVE1XV3SK6+YU7YaG5eulki8rqmp8DMo6DwHAAAAAIA73n0DUN3u3jXzL+xWUm+/7d4Qv6PDWYHx+OOmtwyArFiW2e3GxqSbN51fx8ak8XH3gdkPHpjrPR4zLFsy1RDZVlA0NBT28eaKznMAAAAAALgj1ABQPaanpcuXY1UYFy9Kd+4kL+f1mtDCnoGxfbu0dm3hP6oNlCHLMuNl7JAiPrAYG5Nu3TLVFunU15vgoqfHnC5fNmNsGhvN7Z/6lHT4sFRbu+IPp2joPAcAAAAAgDtCDQCVaWFBevddU4FhV2JcvZr8EXCPR9q8OTbIe9s26ZFHKvvdUmCZHjxwr7KwL7sVO8WrqTEtpHp6THgRH2D09Jh2U3aGODBg8sd//++dcyVaWnhzHwAAAACAakSoAaD8WZbpWRNfgXH5shnwnairyxlgPPaYaZgPYNH8vAkn3AKLmzdN1cRSOjtjoUV8YNHdbQqfMhk/w1wJAAAAAACQKG+hxsjIiCSpr69PgUBAwWBQfX19+Vo9AMTcuxebgWFXYri9y9rc7Awwtm6VVq8u+OYCpSYaNZ3X3EKLsTFzm9tci3itre5VFvZ1Xm9+tpO5EgAAAAAAIF7eQo2TJ0/q1KlTkqQ9e/ZocHAwX6sGUM3m56UrV2IhxsWL0nvvJS9XWyv5/Wb+hR1kbNzIHAxUJcuSQiH3Kgt7rsX8fPp1eL3ugYX9ddWqlX8czJUAAAAAAACJ8hZq7NixQ5OTk5Ikn8+Xr9UCqCaWJd244Wwj9c477u++9vTEKjC2bZO2bMnPR8OBMjE97V5lYZ+fmUl//9ra2FyLxNCip0fy+cgEAQAAAABA6cnrTA3CDABZCQadFRiXLklTU8nLtbbG2kdt3y49/rh5xxWoYLOzpqLCLbQYG3PfVeJ5PKbbWmJbKPvymjUm2AAAAAAAACgneQs1gsGghoaGJEnnzp3TkSNH5Pf7XZeNRCKKxA3wDYfD+doMAKUqEjFtpC5ciAUZY2PJy9XVSY8+6qzCWL+ej4yj4iwsSLdvuwcWN29Kd+8uvY62Nvcqi54ead06qb5+5R8HAAAAAABAIeUt1Dh8+PBipYbf79fevXs1Ojrquuzx48d17NixfH1rAKXGsqTr151tpAIB8y5uog0bnIO8t2wxwQZQJGfOSDU17jMbBgbMcOp0sx5sliVNTqYOLcbH3XeJeE1N7lUW9nXNzbk8QgAAAAAAgPKVt3cOA4GA+vr6JJlQIxAIKBAIuFZrPP/883r22WcXL4fDYW3YsCFfmwJgubJ9V/fu3Vj7qAsXpLffNg3/E/l8zkHejz9uWksBJaSmRjp92pyP3wUGBsz1Bw/Grpuacg8s7FNcUaKrurrUgUVPj6nEoEgJAAAAAAAgJi+hxsjIiHbv3r04KNzW2dnpurzX65WXgb5A6Ur3ru7nPy/9w38ofelLsSqM27eT1+H1mtDCDjC2bTNTiXmHFiXO/pU/fVqam5P27JG+8AXpD/9Q2rHDdFE7dMgEGPfvp1+XxyOtXZs6tFizhl0CAAAAAAAgG3kJNfx+v1544YXFy8PDw9q3bx+Dw4FyFf+ubjAo+f3S7/6u9M1vSh0d0muvmZPN45EeeSQ2yHvbNnOZKcQocdGodOeOqaq4cSNWbXHjhnTvnvRLv2R+vS3LBBHXr5tTPJ/PPbDo6TE5Ht3UAAAAAAAA8icvb7X4fD719/frxIkT8vl8Gh0d1eDgYD5WDaDQolFTfRGJSPPz0mc/63xXd90689FzO8DYulV67DGa+6Nk3b/vDCzs0MJuFTU/736/VavMr77HY2Zb7NuXHFp0d0uNjYV9PAAAAAAAANUsb58f7evrW5ypAaDMPHggfec70tmz0htvmOoMyXzE3OMxgUVbm3TqlKnCWL26qJsLxJufl27dcgYW8ZUXU1Pp719ba7K69etNSLF+vQks/vqvzbq9XjPQ+9FH3cfMAAAAAAAAoHBoigFUq/FxE2KcPSt997vOj6s3N0sf/7j5iPv8vPko+vy8FAhIn/hE8bYZVcmypFAouT2UfX583CyTjt0iyg4s7NP69WauRWKntIEB6etfN7MznnkmNiRcItgAAAAAAAAoJkINoFpYlnTpUizICASct69fL+3aZU4f/rD0e79n3sXlXV0UQCTiXmVhXxeJpL+/1+ussrADC7tVVFNT5tti/6ofPBj7VY8fMxN/GQAAAAAAAIVFqAFUspkZZ1upycnYbR6P9KEPxYKMjRvNdRLv6iLv7IHciVUW9uX4X003Ho8Z5eIWXPT0mPn19q9vPrY1/lffZl+ORvPzfQAAAAAAAJA9Qg2g0ty+Lb3+ugkyRkakubnYbc3N0s6dJsT4+Mel9nb3dfCuLnJw7557aLHUQG5bS4t7e6ieHjPzor6+MI/jwIHUt5HlAQAAAAAAFJfHspbqRL7ywuGw2tvbFQqF1NbWVuzNAcqLZUmXL5sQ4/XXpXfecd7e3W1CjKeekp54wgz/BnIwP2/CCbf2UGNjJtRIp7Y21g7Kbb5Fa2thHgcAAAAAAABKT6Y5Ae9uAuUoEpHOnzchxuuvS3fvxm7zeKTt202IsWuX9Mgj+evLg4pmWVIw6N4eamzMFAFlMpDbDisS20StWSPV1BTikQAAAAAAAKBSEWoA5eLu3VhbqfPnpdnZ2G2NjbG2Uk8+ad5ZRkU4c8YEAW5tjwYGTCewdO2SEs3MJAcW8ZUXmQzkjh/CnVht0diYzaMDAAAAAAAAskOoAZQqy5KuXDEhxtmzpsVUvK6u2JDvJ54o3MABFFRNjfts9vhZ7vHsgdxu7aHGxjIfyO3WHmr9epOXUfgDAAAAAACAYiHUAErJ7Kz03e/G5mPcvu28fevWWJDh9/PuchWwg4z4YOPkSemll6Tdu6WGBuk3fiMWWty6tfRA7lWr3NtDdXcXdiA3AAAAAAAAkC0GhQPFNjEhvfGGCTK+8x1n/x+vV+rvj7WV6uws3naiYCIR0xIqvk3UN74hvfmmNDdnQovubnNyU1cXu90tuGAgNwAAAAAAAEoNg8KBUmVZUiAQayt16ZLz9jVrYtUYH/2o+Sg+Ksr8vDQ+Hgss4gOMmzdTt4iamTG/Ph6PtG2be3uonh4GcgMAAAAAAKByEWoAhTA352wrNT7uvP3xx02I8dRT0pYttJUqc/Zci5s3kysubt40XcWWqpFrbnYGFpcuSaGQub6mRvqJn3AfHg4AAAAAAABUMkINYKUEg7G2UufOmY/Z27xeaccOE2I89ZS0enXRNhPZsywTMCRWWthfM5lr4fXGWkR1dzsDjMQWUQMD0ve/L/3CL5ggwx4SLhFsAAAAAAAAoLoQagD5YlnS1auxaowLF5wfx1+92gQYu3ZJfX3mXW2UrPv3U4cWN286Myo3tbVSV1dyWGEHGB0dmRXk2AHGwYOxAMNteDgAAAAAAABQDQg1gHTOnDG9ftzeNR4YMG2lPvIRE2KcPWve7Y736KMmyPjhHzbnaStVMuxh3G7toW7elKam0t/f4zE5VXxgEX9+zRoTbCxXNOoMNGz25Wh0+d8DAAAAAAAAKBeEGkA6NTXJH4cPhaRf/VVpcFDy+cxH7m319aYKw56PsXZtwTcZhj2MO3EIt30+1TDueD5f6vZQ69aZH/dKO3Ag9W1UaAAAAAAAAKDaEGoA6TzzjDQ9Lf2v/2WqMerqpNdeM++Kd3ebQKOjI9ZWascOqbGx2FtdFSzLDOOOr66Ir7jIZRh3fIDR3S01NRXmsQAAAAAAAADIDKEGYJufl65flwIBaXQ09vXuXTNg4UtfMj2HLEv6oR+SfvZnTZDx+OO0lVoB9jBut0oLO8RYahh3Q4MzpEhsE7VqFT86AAAAAAAAoJwQaqA6TUwkhxfXrqV+l7yvz1RqtLRInZ3SN79Z2O2tUNPTyfMs4r8uNYy7psa0gXJrD9XdbX5UhBYAAAAAAABA5SDUQGWbm5OuXk0OMIJB9+Wbm6XeXsnvl7ZsMV83bzbzM27cMO2n5ufNkPAKHmiw1Hz0aDT9rAdbJCLdupU8z8K+nM0wbrdKi3wN4wYAAAAAAABQHgg1UBksy7SJGh11hhfXr5t34BN5PNLDD8cCjN5ec+rqSv5o/8CAGRZ+8KB5l9++LFVssOE2H11yPhWSyXdu33YPLMbGTEHMUtrbU4cWhRrGDQAAAAAAAKA8EGqg/EQiseqLK1diAUaqj/23tiaHF488Inm9S3+vxEBDin2t4GDjmWdMTvTiiyaY2L3bPBVf/aoZJ/Ld70p/+qcm0HDLjOI1N6duD9XdbW4HAAAAAAAAgEwQaqB0WZY0Pu6svBgdld5/39yWqKZG2rgxOcBYvTr3wQrRqDPQsNmXl3pHv0RYlnTvnhm8HQpJ4bD7+cTLd+5Ix45J/+W/mHV0d5sfyfh4bN32MG63SguGcQMAAAAAAADIJ49lub07XFjhcFjt7e0KhUJqa2sr9uagGB48kN591xlejI6aSdJufD5neOH3S5s2mXfYK1x8QJEukEg8n2v+8v3vm7kVjY3Sz/+8M7RgGDcAAAAAAACAfMg0J6BSA4VlWWbYQmL1xY0b7svX1ZmwIjHA6OwsyObma2B2KqkCiqUqKXINKJqbzQyL9naprc39fPzlL3/ZPE57Pvq6ddLP/EzujxcAAAAAAAAAloNQAytnetoZXgQC5vTggfvyq1cnhxcbN5p31Isk04HZUiygyLS1Uz4CilRhRKrgIpuh2wMD5lRF89EBAAAAAAAAlDhCDSxfNGoqLeKHdgcC0s2b7svX15tB3Vu2xAKMzZtNS6kSYVkme9m71zy03/xN04Zp507p61+XvvUt6e/9Pen8eenP/3z5AUVTU2bBRK4BRbaqdD46AAAAAAAAgBJHqIHsTE0520YFAmYWRiTivnxXl3Not98vPfywGdJQQPPzzgqJUEgKBpPPx183Px+7/4MH0uc/L33hC7GB2aGQ9NZbyd/LDigyCScKEVDkokLmowMAAAAAAACoMAwKh7uFBen9951DuwMB6fZt9+W9XlNtYQcXW7aYy62ted80y5Lu33cGEOmCilDILJ8LrzcWPnzjG2YgdkOD9Mu/7B5OtLVVxaxyAAAAAAAAAMgrBoUjc8Ggs23U6Kh09ao0N+e+fE+Pc+5Fb6+0fr0ZQJGD2Vn3IMItrAgGTZunhYXsv4/H414p4fM5v8af93rNfQcGpHfeiQ3MlqR/9s9yergAAAAAAAAAgBzlLdQIBAIaGhqS3+9XIBDQ4cOH5SuhGQmQeTf++nVneDE6Kk1MuC/f1JQcXvj9ZkJ1CpZlOlSlq5pIvC7V3PCl2G2e4sOIxFAi/nxrqwk2spU4X4KB2QAAAAAAAABQHHkLNfbv36/z589LMgHHoUOHNDg4mK/VIxuWZYKKxOqL69edgyJsHo/00EPJAUZ3tyKzHmcocdZ9/kT8KZeGZjU1qasmUoUVhWjzxMBsAAAAAAAAACgdeQk1AoGA47Lf79fw8HA+Vl11znzm26qpkZ459Ymk2wYOf1vRqHTg83G3zc5K1645w4tAwCQOCaKWR2Fvl0LrtynU9aiCHZsVatugUOM6BacbTCjxAyn47VhAkWr+91JaWtJXTSQGFS0tuVVRrDQGZgMAAAAAAABA6chLqDE8PKzOzk7HdZ2dnRoZGVFfX18+vkXVqKmRTg+2SPq2I9gYOPxXOv1yiw7+6FXpi1+UNRrQzOXrCl4NKjTbpNDCKoXmWxSab1Vw/ofN5eYeBZt6FGpYo6B8umetkuWpl656pKuZb1NdXfqAwi2oqKuQaS0HDqS+jQoNAAAAAAAAACisvLz1HHSpCpCkiVSzGpCSCTK+rVNfatXId76pvg/N6tXXPPr2rUe1tema/urNBn3lbKdC849pzvq7H19trdTULDU1mkETq5qkxsaUg7tbW5eePxF/vqmpNKsoAAAAAAAAAADVZUU/T58q7IhEIorE9TUKh8MruRll55lTn9DY9df0/31jhwa+Z8mSR931dzVT36Yr1loTWjQ1SY1Nqm9rlG9tvXw+T0bDstvaTAYCAAAAAAAAAEC5yUuo4fP5kqoyJiYm5PP5XJc/fvy4jh07lo9vXbE+8+JTOrn5jmo8lhpq53Xs11vU3rtJvjV1jqDC66WKAgAAAAAAAABQHdz7E2Vpz549rtf39/e7Xv/8888rFAotnt577718bEZF+ZP/OqJNrXe1ue2O1jWF9eD7V/QPf6xOH/uY9PjjUne36TBFoAEAAAAAAAAAqBZ5CTX8fr/jciAQUH9/f8pKDa/Xq7a2NscJMQOHv63Tgy06uP++Xr27Qwf339fpwRYNHP52sTcNAAAAAAAAAICiydtMjcHBQR09elQ7d+7UuXPnNDg4mK9VV5X4QMMMDY8NDz892CLp24vXAwAAAAAAAABQTTyWZVnF3ohwOKz29naFQqGqr9o485lvq6ZGrsHFwOFvKxqVDnyeUAMAAAAAAAAAUDkyzQkINQAAAAAAAAAAQFFlmhPkZaYGAAAAAAAAAADASsvbTI3lsItFwuFwkbcEAAAAAAAAAAAUmp0PLNVcqiRCjampKUnShg0birwlAAAAAAAAAACgWKamptTe3p7y9pKYqRGNRnXjxg21trbK4/EUe3NKRjgc1oYNG/Tee+8xawSoMuz/QPVi/weqG38DgOrF/g9UL/Z/wLAsS1NTU1q/fr1qalJPziiJSo2amho9/PDDxd6MktXW1sYfNKBKsf8D1Yv9H6hu/A0Aqhf7P1C92P8Bpa3QsDEoHAAAAAAAAAAAlAVCDQAAAAAAAAAAUBYINUqY1+vVZz/7WXm93mJvCoACY/8Hqhf7P1Dd+BsAVC/2f6B6sf8D2SmJQeEAAAAAAAAAAABLoVIDAAAAAAAAAACUBUINAAAAAAAAAABQFuqKvQFIFggENDQ0JL/fr0AgoMOHD8vn8xV7swAUyMjIiCSpr69PgUBAwWBQfX19Rd4qACtlZGREhw4d0vnz5x3XczwAVL5U+z/HAkDlGxkZ0fDwsCTp3LlzevHFFxdf5zkGACpbuv2fYwAgM4QaJWj//v2L/9gEAgEdOnRIg4ODRd4qAIVy8uRJnTp1SpK0Z88e9n+ggtlvWNj/vMTjeACobOn2f44FgMo3PDys5557TpJ04sQJ7d69e/F1n2MAoLKl2/85BgAyw6DwEhMIBBwHMJLU0dGhycnJIm4VgEI6deqUPvWpT0kSn8gCqoTH41H8IRnHA0D1SNz/JY4FgEo3MjKi3bt3L76uBwIB9fb2anR0VJI4BgAqWLr93+/3cwwAZIiZGiVmeHhYnZ2djus6OztdP8EFoHL5fD4OYIAqxvEAAI4FgMrV19enF198cfFyMBiUZF7rOQYAKlu6/d/GMQCwNNpPlRj7j1miiYmJwm4IgKIJBoMaGhqSZPprHjlyRH6/v8hbBaCQOB4AqhvHAkDl27dv3+L5l19+WXv27JHP5+MYAKgCqfZ/iWMAIFOEGmUi1YENgMoTPwjQ7/dr7969i6XoAKobxwNAdeBYAKge9huY8e2mUi0HoLK47f8cAwCZof1UifH5fEmfwJiYmKDsDKgigUBg8bzf71cgEHBcB6DycTwAVDeOBYDqcfToUb366quLr/EcAwDVI3H/lzgGADJFqFFi9uzZ43p9f39/gbcEQDHYQ8MSJfbVBVDZOB4AqhfHAkD1OHHihI4ePSq/369gMKhgMMgxAFAl3PZ/jgGAzBFqlJjEPnmBQED9/f18KgOoEn6/Xy+88MLi5eHhYe3bt4+/AUAViG8rwfEAUF0S93+OBYDKNzQ0pL6+vsU3NF955RX5fD6OAYAqkG7/5xgAyIzHsiyr2BsBp0AgoJMnT2rnzp06d+6cnn/+ef6AAVVkZGREw8PD8vl8Gh0ddRzUAKgsw8PDevXVV3XixAk999xz2rlz5+LgQI4HgMqWbv/nWACobIFAQL29vY7rfD6fJicnF2/nGACoTEvt/xwDAJkh1AAAAAAAAAAAAGWB9lMAAAAAAAAAAKAsEGoAAAAAAAAAAICyQKgBAAAAAAAAAADKAqEGAAAAAAAAAAAoC4QaAAAAAAAAAACgLBBqAAAAAAAAAACAskCoAQAAAAAAAAAAygKhBgAAAAAAAAAAKAuEGgAAAAAAAAAAoCwQagAAAAAAAAAAgLJAqAEAAAAAAAAAAMoCoQYAAAAAAAAAACgL/z9pnVlpa/fbrgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot accumulated phase advances\n",
    "\n",
    "muqx_sum = torch.cat([torch.tensor([0.0], dtype=torch.float64), muqx.cumsum(-1)])\n",
    "muqy_sum = torch.cat([torch.tensor([0.0], dtype=torch.float64), muqy.cumsum(-1)])\n",
    "\n",
    "# Compare\n",
    "\n",
    "print(torch.allclose(muqx_sum[:-1], mu_qx))\n",
    "print(torch.allclose(muqy_sum[:-1], mu_qy))\n",
    "\n",
    "# Plot\n",
    "\n",
    "layout = Layout(ring)\n",
    "_, _, lengths, *_ = layout.slicing_table()\n",
    "rectangles = layout.profile_1d(scale=5.0, shift=35.0, text=False,  exclude=['BPM', 'Corrector'])\n",
    "\n",
    "plt.figure(figsize=(16, 4))\n",
    "plt.errorbar(positions, mu_qx.cpu().numpy(), fmt=' ', color='red', alpha=0.75, marker='x')\n",
    "plt.errorbar(torch.cat([ring.locations(), ring.length.unsqueeze(-1)]).cpu().numpy(), muqx_sum.cpu().numpy(), fmt='-', color='red', alpha=0.75)\n",
    "plt.errorbar(positions, mu_qy.cpu().numpy(), fmt=' ', color='blue', alpha=0.75, marker='x')\n",
    "plt.errorbar(torch.cat([ring.locations(), ring.length.unsqueeze(-1)]).cpu().numpy(), muqy_sum.cpu().numpy(), fmt='-', color='blue', alpha=0.75)\n",
    "for rectangle in rectangles:\n",
    "    plt.gca().add_patch(Rectangle(**rectangle))\n",
    "plt.ylim(-1.0, 40.0)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [
    "myt0_gMIOq7b",
    "5d97819c"
   ],
   "name": "03_frequency.ipynb",
   "provenance": []
  },
  "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.12.1"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}