{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "37e17f59-4312-4476-86b1-090baf77b82a",
   "metadata": {},
   "source": [
    "# Example-17: Line (element)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ff8535ea-9a8d-460d-849f-c6e7bf7248a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this example line usage is illustrated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d69df16e-0fa8-440d-b2c5-c5514c1f7d90",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "from os import system\n",
    "\n",
    "from pprint import pprint\n",
    "\n",
    "import matplotlib\n",
    "from matplotlib import pyplot as plt\n",
    "matplotlib.rcParams['text.usetex'] = True\n",
    "\n",
    "import torch\n",
    "\n",
    "from model.library.drift import Drift\n",
    "from model.library.quadrupole import Quadrupole\n",
    "from model.library.sextupole import Sextupole\n",
    "from model.library.dipole import Dipole\n",
    "from model.library.line import Line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d434635d-f8d4-4c24-8ec5-96312a79629f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define unique elements\n",
    "\n",
    "QF = Quadrupole('QF', 0.5, +0.25)\n",
    "QD = Quadrupole('QD', 0.5, -0.20)\n",
    "SF = Sextupole('SF', 0.25)\n",
    "SD = Sextupole('SD', 0.25)\n",
    "DR = Drift('DR', 0.25)\n",
    "BM = Dipole('BM', 3.50, torch.pi/8.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f46db43e-2c35-405c-b141-d413f0cf18d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a line\n",
    "# With propagate flag, values for other flags will be propagated to all root elements (but not to lines!)\n",
    "\n",
    "FODO = Line('FODO', \n",
    "            [QF, DR, SF, DR, BM, DR, SD, DR, QD, QD, DR, SD, DR, BM, DR, SF, DR, QF], \n",
    "            propagate=True, \n",
    "            dp=0.0, \n",
    "            exact=False, \n",
    "            output=False,\n",
    "            matrix = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "62bc353f-a7c1-410d-ae7c-bb0f38ab8a09",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([0., 0., 0., 0.], dtype=torch.float64)\n",
      "tensor([0., 0., 0., 0.], dtype=torch.float64)\n"
     ]
    }
   ],
   "source": [
    "# Similar to elements, lines have to variants of call method\n",
    "\n",
    "# 1) without deviation variables\n",
    "\n",
    "state = torch.tensor([0.0, 0.0, 0.0, 0.0], dtype=torch.float64)\n",
    "print(FODO(state))\n",
    "\n",
    "# 2) without deviation variables (defaut dictionary)\n",
    "\n",
    "state = torch.tensor([0.0, 0.0, 0.0, 0.0], dtype=torch.float64)\n",
    "print(FODO(state, data=FODO.data()))\n",
    "\n",
    "# In both cases, it is possible to apply alignments errors using alignment flag\n",
    "# Alignment errors can be applied to elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f2ac0d8f-0a5c-4af8-8485-c0886bd9afaf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'BM': {'dl': tensor(0., dtype=torch.float64),\n",
      "        'dp': tensor(0., dtype=torch.float64),\n",
      "        'dw': tensor(0., dtype=torch.float64),\n",
      "        'e1': tensor(0., dtype=torch.float64),\n",
      "        'e2': tensor(0., dtype=torch.float64),\n",
      "        'kn': tensor(0., dtype=torch.float64),\n",
      "        'ks': tensor(0., dtype=torch.float64),\n",
      "        'mo': tensor(0., dtype=torch.float64),\n",
      "        'ms': tensor(0., dtype=torch.float64)},\n",
      " 'DR': {'dl': tensor(0., dtype=torch.float64),\n",
      "        'dp': tensor(0., dtype=torch.float64)},\n",
      " 'QD': {'dl': tensor(0., dtype=torch.float64),\n",
      "        'dp': tensor(0., dtype=torch.float64),\n",
      "        'kn': tensor(0., dtype=torch.float64),\n",
      "        'ks': tensor(0., dtype=torch.float64)},\n",
      " 'QF': {'dl': tensor(0., dtype=torch.float64),\n",
      "        'dp': tensor(0., dtype=torch.float64),\n",
      "        'kn': tensor(0., dtype=torch.float64),\n",
      "        'ks': tensor(0., dtype=torch.float64)},\n",
      " 'SD': {'dl': tensor(0., dtype=torch.float64),\n",
      "        'dp': tensor(0., dtype=torch.float64),\n",
      "        'ms': tensor(0., dtype=torch.float64)},\n",
      " 'SF': {'dl': tensor(0., dtype=torch.float64),\n",
      "        'dp': tensor(0., dtype=torch.float64),\n",
      "        'ms': tensor(0., dtype=torch.float64)}}\n"
     ]
    }
   ],
   "source": [
    "# Values of deviation variables can be passed only to unique elements\n",
    "\n",
    "pprint(FODO.data(alignment=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "652863b1-3d6e-48b0-ab0b-e0c41503e31c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([512, 4])\n"
     ]
    }
   ],
   "source": [
    "# Mapping over a set of initial conditions\n",
    "\n",
    "state = 1.0E-3*torch.randn((512, 4), dtype=FODO.dtype, device=FODO.device)\n",
    "print(torch.vmap(FODO)(state).shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2a39ba62-32bc-49d3-93a7-f86ba15ce44c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([512, 4])\n"
     ]
    }
   ],
   "source": [
    "# Mapping over a set of initial conditions and parameters\n",
    "\n",
    "state = 1.0E-3*torch.randn((512, 4), dtype=FODO.dtype, device=FODO.device)\n",
    "dknqf = 1.0E-3*torch.randn(512, dtype=FODO.dtype, device=FODO.device)\n",
    "\n",
    "def wrapper(state, dknqf):\n",
    "    data = FODO.data()\n",
    "    data['QF']['kn'] = dknqf\n",
    "    return FODO(state, data=data)\n",
    "\n",
    "print(torch.vmap(wrapper)(state, dknqf).shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "bcd03909-f051-46cc-9e93-47f4d8978b02",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([[-0.4395, 15.4433,  0.0000,  0.0000],\n",
      "        [-0.0522, -0.4395,  0.0000,  0.0000],\n",
      "        [ 0.0000,  0.0000,  0.4963,  5.3596],\n",
      "        [ 0.0000,  0.0000, -0.1406,  0.4963]], dtype=torch.float64)\n"
     ]
    }
   ],
   "source": [
    "# Differentiability (state)\n",
    "\n",
    "state = torch.tensor([0.0, 0.0, 0.0, 0.0], dtype=torch.float64)\n",
    "\n",
    "print(torch.func.jacrev(FODO)(state))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "3939a00c-0e10-4f3e-ac20-8c294186f70f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([0., 0., 0., 0.], dtype=torch.float64)\n"
     ]
    }
   ],
   "source": [
    "# Differentiability (parameter)\n",
    "\n",
    "state = torch.tensor([0.0, 0.0, 0.0, 0.0], dtype=torch.float64)\n",
    "\n",
    "dknqf = torch.tensor(0.0, dtype=torch.float64)\n",
    "\n",
    "def wrapper(dknqf):\n",
    "    data = FODO.data()\n",
    "    data['QF']['kn'] = dknqf\n",
    "    return FODO(state, data=data)\n",
    "\n",
    "print(torch.func.jacrev(wrapper)(dknqf))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "768e19cc-c1f0-4324-bf38-3a1d6ef20da5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([[-0.4395, 15.4433,  0.0000,  0.0000],\n",
      "        [-0.0522, -0.4395,  0.0000,  0.0000],\n",
      "        [ 0.0000,  0.0000,  0.4963,  5.3596],\n",
      "        [ 0.0000,  0.0000, -0.1406,  0.4963]], dtype=torch.float64)\n",
      "tensor([[-7.5649, -3.8172,  0.0000,  0.0000],\n",
      "        [ 0.4176, -7.5649,  0.0000,  0.0000],\n",
      "        [ 0.0000,  0.0000,  2.7423,  1.3262],\n",
      "        [ 0.0000,  0.0000,  0.5426,  2.7423]], dtype=torch.float64)\n"
     ]
    }
   ],
   "source": [
    "# Differentiability (composed)\n",
    "\n",
    "state = torch.tensor([0.0, 0.0, 0.0, 0.0], dtype=torch.float64)\n",
    "\n",
    "dknqf = torch.tensor(0.0, dtype=torch.float64)\n",
    "\n",
    "def wrapper(dknqf):\n",
    "    data = FODO.data()\n",
    "    data['QF']['kn'] = dknqf\n",
    "    return torch.func.jacrev(lambda state, data: FODO(state, data=data))(state, data)\n",
    "\n",
    "print(wrapper(dknqf))\n",
    "print(torch.func.jacrev(wrapper)(dknqf))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "228bb219-b5b7-4ad9-b441-47ee0484a6ce",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 8.17 s, sys: 2.67 ms, total: 8.18 s\n",
      "Wall time: 8.18 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# Tracking\n",
    "\n",
    "# Note, in general tracking is slow\n",
    "# This is especially the case for exact integration of elements with large number of slices\n",
    "\n",
    "qx = torch.linspace(0.0, 0.01, 16, dtype=torch.float64)\n",
    "px = torch.zeros_like(qx)\n",
    "qy = torch.zeros_like(qx)\n",
    "py = torch.zeros_like(qx)\n",
    "\n",
    "state = torch.stack([qx, px, qy, py]).T\n",
    "orbit = []\n",
    "\n",
    "for _ in range(2**10):\n",
    "    state = torch.vmap(FODO)(state)\n",
    "    orbit.append(state)\n",
    "\n",
    "qx, px, *_ = torch.stack(orbit).swapaxes(0, -1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "57a0935c-aa05-45ab-abf9-cb60d19ae51b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAFeCAYAAACB7binAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABV9ElEQVR4nO19zW4cR5J/dLcoeagZiWzdxuNug8Be5qhpvgH9BoQX1H05WsxtB+DA8MHYw8DbB98WaNIPIGFEvYH0BiS1T9DESJr1zWyStrTSWGT8D/5XOzuZVRWREVlV3YwfUNBHZX10VWX84jtbiIhgMBgMBsP/R7vuGzAYDAZDs2DEYDAYDIYZGDEYDAaDYQZGDAaDwWCYgRGDwWAwGGZgxGAwGAyGGRgxGAwGg2EGRgwGg8FgmMGNum+gabi8vITvvvsOfvOb30Cr1ar7dgwGg0EMRIQffvgBfvvb30K7XW4PGDF4+O677+CTTz6p+zYMBoNBHa9fv4bf/e53peOMGDz85je/AYCfH+CdO3dqvhuDwWCQ4/z8HD755JOpfCuDEYOHzH10584dIwaDwbBQoLrHLfhsMBgMhhkksxiOj4/h6dOnsLa2BsfHx7C9vQ0rKyvssbH7AACeP38Ox8fHsLa2BgAAGxsbqX6uwWAwLA4wEe7fvz/9+3g8xs3NzaixsfuePXuG29vb031ra2uk+z47O0MAwLOzM9J4g8FgaDq4ci2JxXB8fDzz77W1NXj+/Dl7bOw+AIA//vGPcHR0NN337NmziF9iMBgM1w9JYgzPnz+Hbrc783/dbhdevHjBGhu77/j4GE5OTmBlZQVevHgBp6enU3eSwWAwGIqRhBhOT0+D/39ycsIaG7vvxYsX0O12p/GHb7/9Fp4+fUq4c4PBYDBUmq6aJ8y5Y8v2nZycwPHxMWxsbMDKygpsb2/D6uoqYGAV0/fv38P79++n/z4/Pyffo8FgMCwiklgMKysrV6yDzLXDGRu7b21tbTomuwYABF1ZX3/9Ndy9e3e6WdWzwWC47khCDHlpoYPBgDU2dh8nnvDFF1/A2dnZdHv9+jX5WIPBYFhEJHEl+YL5+PgYBoPBjOaeafZFY30Lg7NvMBjA6ekprKysTGsZ7t+/f+Veb926Bbdu3RL/ZoPBYFgUJIsx7O/vw1/+8hdYX1+Hg4MD2N/fn+77+uuvYX19HXZ2dkrHSvf94Q9/gKOjI0tXNVSC9fV1ODw8hOXlZXj79i0MBgP4l3/5F3jy5AncuXMHJpMJrK6uwvn5+fTfvV4PfvzxRwAA+PWvfw2vXr0qHPPXv/4VHj58WOfPNCw4WhiKyF5jnJ+fw927d+Hs7Mx6JRlmkAn9IkF+cXFx5bhOpxP8/1gsLS3BTz/9dIU8XDI6ODhQu55h/sGVa9YryWDw8ODBA7hx4wZ0u11otVrQ7Xbhxo0bcHh4CAAAr169gpOTEzg5OYFXr14BAMBkMpkR/svLywDwc8zr888/h06nA6urqwAAsLq6OvPvXq8H3W4Xut0u9Hq90jE//fTTzDUnkwkAALx9+xYAAA4PD+HevXtw79496Pf70Gq1YH19PekzMywWzGLwYBbD9cKDBw/gyZMn8PHHH8M//vEP+Oijj+Ddu3dweXmZe0yRxfD555/Do0ePkt5zZrnkWQyZReGj0+lccWdVcb+G+sGVa0YMHowYFhsZEWQCMoTl5WV4//793ArR3d1d+PLLLwHgF+LKQ6vVglarBR999BF88803FrtYUBgxCGHEsFjY3d2F//qv/wJEzBWQnU5nxmJYRAHpE2JGdpeXl9PCT58QLVaxODBiEMKIYTGQCUJEvOIWcl0q1134PXjwAP72t7/lutDsWS0GLPhsuHbY3d2FTz/9FNbX16HT6cDt27fhb3/7G1xcXFwRdIPBAD58+AAnJyeAiNde0D169AguLi7gzZs38K//+q8zAW8AmAluHx4ewo0bN+DBgwd13a6hIhgxGOYW6+vr0Gq14E9/+hO8fPkSDg8P4fLyEt6+fQsfffQRdDod6PV60Ol0YGtry4igBI8ePZohza2trSBRPH78eCZby4hi8WDEYJg7ZOmkWfpoZhUsLS1Bu92G5eVl+Oabb+DDhw/w8uVL+PDhQ61B47z0V//f/X4fOp0O3Lx5E27fvg337t2D3d3d2u47jygyZOmyjx8/NoJYNKguE7QAsBXcmovBYIAAcGXr9XrY7/dxNBpVej9bW1vY6XRwdXUVAQBXV1dn/r28vBy8X87WarWw1Wphr9crvNZgMKjsd2fvodVqXbnXdruNy8vLlb8LQzG4cs2IwYMRQ3MwGo2w3+/jYDDATqczI4Q6nQ5ubW0lv4etrS1st9u4tLSEy8vLuLy8jDdv3pwKbKqALxLq2b97vd7MtdrtNptIqiQM//24z2N5ebmyd2QohxGDEEYM9SPTxPMEYwphl10zE87ZViT8W60WyWKIvd/RaITLy8ski4FKSimJIiNRn9RarZYRRM0wYhDCiKE+5LmKUlgIvhuoSPjnWQxNEnbZs6MQRhVEsbW1FSTRJj2z6wQjBiGMGKpHJqR94ZVpyRrChEIEIYth3gUZ9XenIArX4vEtrypjIgYjBjGMGKpDnoWgKTRCpFO1i6VJqIso8sh/3ol3XmDEIIQRQ3qECKEqyyDTXK8LEZSBShRaArwKZcBwFUYMQhgxpEOe9q4hFLLAZ55gGwwGtaS0zhvKiEIaUPevYxZENTBiEMKIQR+ZEND0M49GI+x2u7i0tBTMHrpuLqJUKHPFSQW5WRDVwIhBCCMGfbhCWypMMkLIyx6qq7jKzQrKMnCye8sC2WV1DNx/V6ldZ78vr2gvhQXRtMyveYYRgxBGDHrIJns24ZeXl0Xny9Muq8geyksHzf6d55+vYsuI9+bNm9jtdishxjxXkNRSC71jq6SWw4hBCCMGOfzJ3Wq1RP79kBBaXl5OKgR9IqAK6SothrL7qCLrKkXtSZbmau4lPRgxCGHEEA/tgGJV/uetra2p1p0JcarArcOtk8HvWXTz5s1CskhJFGVB65jrhRIKLEAdByMGIYwY+NAW4KkzVsr85SEimBdh1JQK6NA7jLXwtC3Q6wgjBiGMGOgIBYJjhGjWjK3X6yW1EIoybHyLYV6IgApKvUL2DDR9+iGlQaowaMWsrhOMGIQwYqCj3++rTPiQ5q6hpVeVkz+PoBCFljWhpUBk0Mxyuy4wYhDCiIGGTBPsdDpRLoJMWGgHGbN4Qcrc+0WEG2fJ26QEoeVyTFEXs+gwYhDCiKEYIXdMDEKEIBHYRdkxVuzGQ5E1oeFu0qqAt+ppOowYhDBiCENb29MISuadL7MY6iiQyiyhLJsmI8BQMLgoHdUdU7ewKyLdFGTO/ab883Q6neh7WlQYMQhhxDCLLDAs1cpCmSWxAs/v9V+1O8EX/pxFc6RbVifhau9VkUeRu0ny3LUC1O53YdbhLIwYhDBimEWm8cZ2JdU09/M0zJRCkZP+6W7aFgOHPNzV71IuxBOy1CQuppDmz32vVX0X8wYjBiGMGH6GX33abrdFx0s1wTzhmGLtBq7wr6LK2O/FFLIYiogi1b3lkXWsG0/qXrKmfGEYMQhhxBCeXJxJHnL11CEkypBHXqGtqhYTEriZYkXV2ymsCW03k9TS1EqSWBQYMQhhxDBrjnM0v5Ag19IctdwCruZdpl3HElqTUGYFaZNdRhA+SWi5D7mBZd8VOO/vMxZGDEJcZ2LwhQhVUGjGEULn0hJYZb76JloBmvCD5nmWkdaz0OyUKg0su+9+0d9zCEYMQlxHYvCblXG1Mg1BnsJtVKQtX3cNEpEeU9GoQNd4rxLFw7+H6wYjBiGuIzH4gp0z4VyBHitA/EkrEdhZ3CDP332diaAIFJKQknTIUuGcUyswfR1boRgxCHHdiMENGnLM/NAk5UKzwVp2rpDwqTtgnJGVuyaDGyB2M42a0sSvqKhN4z1J3EshtyDneV03FyKiEYMY14UYYoN6GpqfZkwib6lPiT+be3231iNUp1DWg6hsC9U3VKn1pqgfkWa++efguD81lJp5gxGDENeFGGImuIavWMNKKHJ79Hq9JL36Q9lMMZXOVIuBc06/niKVtZFXcR6beaapZMRYhe63uOhWQ2OIYTwe43A4xP39fRwOhziZTKLGxu5zsbOzU3h9F4tODH5mCmdCuJpv1eZ/hjztNcUyn25dAEVA51kMMQLb/Z1FFkMRAaVYXyFDnpIQQ8oa30asJaNlATUdjSGG+/fvT/8+Ho9xc3MzamzsvgxHR0cIAEYMeFWo9vt90nGSBVI0tMLQvbuCUjO1kiJw3b/XKUxCPZso915FgaAW6UsqninHXpcGfI0ghvF4PCO0ERFXVlbYY2P3udjf38e1tTUjBpzVjqgatjtxuEsqalRA58UQNMmgTKDGrjlRB0KtM4pIYmlpSfy7NFuWSK2H2LTURW/A1whi2Nvbw42NjZn/W1tbw6OjI9bY2H0Z9vf3p/9/nYlB4oeNmaBahU2aWUv+OfMEmasBL4pbgVL9rPFs86qepXUtXIUkNi3VfyaL8v4R+XKtDQlwenoa/P+TkxPW2Nh92XlXVlZK7hTg/fv3cH5+PrMtEnZ3d+Hx48dwcXEB5+fngIhwcHBQetyDBw/gxo0bsLq6CgAAg8EAHj16VHrc+vo6HB4eTv/d6XRga2sL3rx5Aw8fPiTdc3Zt9zzdbhdGoxHp3n3s7u7C7du3odPpTM95cXEx3b+6ujpzn/izwkT6vfOAR48ewYcPH+Dk5AS2trag0+lM3yvAL8/i8PAQOp0OdDodePDgQdR1Li8v4fLyEgaDwfT/Dw8PodVqwfr6OvlcBwcH03ttt9vw8uVL+POf/0w+FhHh7du30+tTru3e88XFBTx58oR8v4uGJMSQhzxhzh1L2ffkyRPY2Ngovc7XX38Nd+/enW6ffPIJ+R7nAe5k+vzzz0nHrK+vs8kkJMwHgwF8+PCBJGB3d3fh008/hX6/P722e57vv/+eTCz+Pf3pT3+Ct2/fwuXl5XRfp9OZks3JyQn5PucdLkmMRiPo9/vQ6/Wm+zPB/vjxYxFJHBwczAhagJ8F9O3bt2F3d5d1r9l7e/v2Ldy7d498vE9ON27cKPwtGaFsbW1Bq9WCi4sL1v0uFFKYLXt7e0H//7Nnz1hjY/c9e/ZsxnVU5Ep69+4dnp2dTbfXr18vhCvJ99VS2mZL/LsS90FeGqSk11Ke22SR3AOayJIE8lJlNddZ4L4HSZq0f21q+3h3HlCTNJqMRsQY8gLDIeFcNDZ2X0Yc2QYAuLOzE4xx+FiUGEP2YXNWSovteeROPm7GUqjQjUssoftwt16vV9kKZ3mL8Wj+PXVwNAVJhLLTOOcJJSFwEgK4sSP3GTShEl2KRhAD4tVUUjdQfHR0hOPxmDQ2dp8LAJi5XhEWgRhi+hfFHBMS7JKMpVarFd1507+PlG0wygigis0VsimFV166scSik5wnNiFhNBpFHed+V1VU0qdCY4hhPB7jzs4O7u/vXykw29zcxOFwSBobuw8RcTKZ4HA4RADA7e3thbcYQlpVGWIXNJFkDPmTNKYwKhNYoZTMVKuUlS2Ak9JKoBJPKssibyEeLeuBc68x33l2Xe637hPZvLqVGkMM84p5JgZ/snA1oljti6PxjUajK1oYF5oxiRBCVcdUAqjDzUMlKy13Wt7zj6lZcI+XVIfHtsK4Lm4lIwYh5pUY3I+93W6zi9eoH7sv2GMnZKvViioa8y2VTIvVqm8o2lyLoWlFUNRiPS0SLXIxSVtpx9bZcI51v2Fqr6fsmHmsjjZiEGJeicGdHGXCVsN91Ol0WELdFwCcY9221ZrCjdILqe62FxK4GnnKdRbqtB5iiSV0Pc4xTVMMymDEIMQ8EgPXpI5xH/kTSbJuA3VSZQSm0WvJv58iQTmPJEABxSqSCDz//Nz4g2avJAr8lumU9+5eI0UX31QwYhBiHomBMyFi/KudTmcm8Ji6sVnesbFCu6jZW6YtzksvJC0UkURsl9S883Lfma+EcPt6ca08V1EqO8a3MrvdLukadcOIQYh5IgZuD6SYBmMuIVAnWyhrRFLwJqlDyHN1zLOLSBt5zyiWiEPxB+46GRKlghsDc38/tQDO/b7nwa1kxCDEvBCD+zFTg2Fcrd+9BjV7qI7WyT78dtTuFruoDBd5QXJXALnvz/+3lnuHg7wuqbHZONL4Q0jBoBCLf13KfWcWZbvdJikMfmZe02HEIMS8EIMraKjmr7twSxligsWxLoDQsTHa/Gg0wqWlJTXNl3LP7spreYJdY/NJLnUBX17L7th1NCSJAzGV9X6iRBlGoxH2+/3pfVKOiU2XrQNGDELMAzFwNXl3QlLNee4kDhWtUeGTUMwkC2mn2qml3JRQLYuBsqUobMtrjRHj2ssjbeq9xqweGNNunpt55P6WJsOIQYh5IAZXcJRNEK5WEzOZRqPRjPDgtNSQpqDmBVE1hKPf15+q0WsGs/PcYtTV2mKeadGzkD7n2EaNkkZ63Gtxxseu/VA1jBiEaDIxcLMuuMFmrvkdugY1S0MaTwi5O7JFYiSCsMgdlVI7jwG1TkGLrPLShzUaHlKVlhhFgusWjZkHnPPXASMGIZpMDFyzlZOdEdPmwp9wMcuFZtoWJxYRcrVoFWmFAtbztKpbWdtxjeeVR57cuoWYIDHi1e+Ha3FQ3K9cyzkmDlIljBiEaCoxcHyfMfncrvCLKXjjuI9iBZRmbUPoXjTP3RSUVXdrr7MgeaecFvExKaauhUlVYjiWACchpGoYMQjRVGLg9GnhWBa+eR6zkAlFEGi4jvxjY4qwysiAm28/TyirfI5NDQ7FYDiCMWQFUp4/N8U0lPnGcSuVWQIxKeRVwYhBiCYSA6damWNZxPTGj+k0qdlFU+IGKav2vS4oSkeVEITvXuJaIy7BUNwxfoop5d792gjKdTj1P+6zbdL6DUYMQjSNGDiFNNygGTcF0BfwlN70sUI9m2A3b94UuXbyUi6rJgNKqitHwGmhzIKSNsPjvLPYlOcYtybH7cNN4mjisqBGDEI0jRg4pfdczZ+jNcUEjGPXbdDwXRcJvBRC179nyhoJ1C1EJtq/QbMtRih7KFbIx1TNczulAuimfcdY1alhxCBEU4gh87tmGnNZGijHheQL7FQCPqafTIgUuGmWoXNoTtIiKyRGyEvIJIWVIW1lkXeelFlH/vVStHxx74mCJq3fYMQgRFOIwTV1KcFQzkfITb/kZiwhzk7u2NqGmzdviglBS2CWBW4pQj5Wa46ptNYiQP+80lbaqftmcZ8BJ37AifUhNmv9BiMGIZpADBI3D8eFRMlAcs9NGa/V+IwzkULX1JiMZf537YWDOPcV6tHkbxorzhURYmyKMrelBeeaXJenJH5AKabknDsljBiEaAIxcFpecLs8xvSejx1PFUiSrKMUtQ1FGTtVkwAVlJoM6XORWEGh90StUud+H9wkCT8llUMkZRY618pIBSMGIeomBs6HxF1/mfuRSiYjJS4Q0vQlAerYlg9UN1HTyCAPblFbnkUR+1tCz4pjOYRINyblmTKeEwDmWCbcld+aEGswYhCibmKgfvx+FgY31bQMMZXWEiuBo32GtGONfj3+1qQ8dAk03EGUc8a6lqjHxmjfVMWJawm45y4b24RYgxGDEHUSAyclzv2IKcKLM0H8TqFlkGp/HE1fIpDc62sXeM0DtAlCQtC+C1T7O867x6JrcC0BjsDnKGUpYMQgRJ3EQP14uKmp1I89NFk5cYWq4wkxwkw7JjGPKCIIadYR51n631u73WZbvlxLg5LMwSEf6pytO9ZgxCBEXcTAEbAcn6Xrwy8b77cK4BT9UDQhrbz2GK2+SkKgBoK5m/a9avZOiiX80Wh0pQI/xbW4bS1SCPw6Yw1GDELURQwpPkR/0hWN59YdcCenpDVGrMDKa/AGwK+RyAO3vkFz01q/usi1Jqk2p1ofMQWUkpgZ5RoxAr9sbJ2xBiMGIeoghtjYQhH8eoUy36p7Xq45z9XAqrIU8khBKkzrJIOyTfLbiog09TuLqWPhume4NTkxweWy83LITBNGDELUQQzUjyXG3QRQbgFwW1dIFgCKtRQ4Ak+SM5/3GyjLe8YSGAXU+gotkqgjayzmnXOD0Zyqf46GTz1vXbEGIwYhqiYG6sfHKWRzWyAvLS0VWgASFxK3nUas1ikteNNMydQWwLHgkFXZN5AHSYuRWHKQFrNxFRVO5TI1caPMwqgj1mDEIETVxED9SDj9ityxRVWfnFS+0PgycE330DUkLowYoZ26gjglNIPJZeeNfS8x8QPtCn2O68evF6Kel1ozUVWswYhBiCqJgap9cwQsJ14QWzUdMz5F64O842ImXBkhaAV6q8BoNLqyjoXUgpAkAWgUQHK+N4oSksqlRFWcOAqWBowYhKiSGGI+orKP2BUIZaZvRgy9Xo91r1wNLpV/OXQtrjZfpmXPCxnkoYzwuL/Pj3OkjBlxj+FY1RKXUhFilL0qrAYjBiGqIgZObCFWgBeBU9/A0Zik2SWxRMKdYEWEsIiVz2UEIQksp4wdcSxPrrCPdSlp1StUaTUYMQhRFTFQP55+vz/9eDhdIlNkUXDbY6SyFLQzZlISgrTYTfueQmszZxuntkMrLZVrOXAEONfKKBvrWkvUDqxFz6RKq8GIQYgqiIGjffR6vemE0OorT81E4k5iSXC6CkshdfVzVfUNq6urSe+X+jyqyjxC5LmJOETiWhmtVqtwrL94VhGo84AzXyQwYhCiCmLgfAzUsdSS/1TBae54rmYX0nQlVblSLS1V24uYLUXBXtUFbZT3T/1uEXmWQIouqZlCV3b9qqwGIwYhUhNDrL+emv5W1muGqnlxYxscKyim6C12tbQ8AZ66tkGDgGKuRXlXIUifUxVZaO7YpaWlwrGcb4zaNgYxzgVcZom495kKjSGG8XiMw+EQ9/f3cTgc4mQyiRobu+/o6AiHwyEOh0Pc3NwsvL6L1MTAKW6hjqWu+MZJ6XMJhKNJUT5uP5WyDDEBbb9PlIQQONZBal8xhyy4v1VS0BY6PsZy4IzVcm2647SsBk7BWxXV0I0hhvv370//Ph6PcXNzM2ps7L7hcDjzd3dsEVITQ7awe5mfOCbtTauKE5GuRUkrVbmpr1SNOBRg5QptihCuO52Vco8av5sa1/CrsSnX5oz3763MpUSdb7G1CkXXT2GJxKIRxDAej68I4pWVFfbY2H1HR0cz1xuPxwgAOB6PS+89NTFQtSPtcRxfJpVspEJem0RCx3CEWgbXP9xEMsiDVkpqhhj3UMyaHlyLkOM6jV1bgaoMFV0/VTpsDBpBDHt7e7ixsTHzf2tra3h0dMQaG7sPEXF/f3/6/0dHRwgAJHdSSmLQLpmPGcchkLKxnEmXmkRCx3CF4SIVu2Wasr+1Wi1W1XPoPJTiMW6iAPd9U7Vx/7ujulrLNHdqsgc1CO7GJFJYDY0ghuFwGBTaz549Y42N3edjZ2fnytgM7969w7Ozs+n2+vXrZMTAjRlQP86y4Fbsurdlk5NqpscEmzUsBWka67yRQQgadRuxz5YjlBF55MBRNKjp2f55NZQtaqyBm+zBRaOJwdXiKWNj97mYTCa4traWay189dVXwUmjTQxbWz+3TS7ruRMbMygCVYCnWuXK1ZrKSCykbZYhJLioyCOF2I6kMdeqgozyLAgJeVLu032XZZlEoesUgfNdcVrLU69PTUmlkk3KOEMjiGFvby/o/w9p80VjY/e52N7eLowtVGUxULV76kdJ9bNy2nVT75Gj2XH81DGaaazASlnslieEpZu0sE3DeuAeF7M6mzueusAUReGgnpc6XzjuH8o5OQFwLhpBDHmB4ZDWXjQ2dl+G4XA4JYXJZFJbjMFNnSwLllG1e2qKKlVTiknt4wp6jgUUQwrUpSS1SaGuYrdY4ZF3v7FxnJRJB2XzhSNMXeulqHI5Ra0PtYMBlZS4aAQxIF5NJXXdPkdHRzNafNHY2H37+/tT62EymeDe3h7pvlMQQyacOUtsFmkgMamsZb5VahyC426SuKa4mSzUiRRa3CZ1bQP1eq6g0Ton9zekKmjj1rpwUj2p90EV5BwrW7vnWarspMYQw3g8xp2dHdzf38ednZ0ZbX1zc3OmzqBobMy+LD3V3fLSZX2kJIYy4UzNYKB+tO44ajYGtb6hzHTn3CdXo4zRWvOax3EnH9VNpO0KoMYnuNcNkQM1+Mm9NmdsbDYRZ05QXaAamYEc6yJFrKExxDCvSEEMmRlJ1VKKcp61P0SOdsSZgBx3AFfAcwV7KLe+7J58lBGCVmM7rfvhEl6MJRUiK854jsXJ+YZSkAhVUdNQ6KjuZA6MGIRIQQyUD8KNBRR9XDFLgRaZrq6JW5Q1wkll5VSnciYqVxDlHcPRqssEcOoWGGUosySo95fnxuIGizUrl2O/uTLLnGoNUOsaYsZpEAgHRgxCaBMD9SPM+gfdvHmz8HxUbYI6+ahBsRQVpKmD0yFLQaO2oW4yCKEsPkG951DVdxn891j2bXJSWDmxiRQpqdT5xk0aoRCIZsq0EYMQ2sRA/QAp46huH6obhzPpYno8cX4zV9OMCU7HFso1mRB8FBFEypoFjptIsqxmWXBZc76lGEdVxKhxSSqMGITQJIaYAFaRIKdqRNRU1tiK6DxIKlfLUCcpzGP1s7SaO8Zt57vdOM0Xi8DpOeSeUyPhQnucO+eKXLzUTEYqjBiE0CQG6odPWTLQ1Ya03E2UcZx2FtRxXO2SUwkdSwopCt6omURlm8RSkZBdTHyG86449xOzFC0nYK0xjuIm4iSPUGufKDBiEEKTGLKPmerOKRpHXfgjtlAtD66VkmKilQkEzgI9sT19Qhk5McJYiwiKtthMlVAQnXIuruXAcRNJ3E9FoNbPaBMD12rQ6oVGgRGDEFVbDNQFxql9WajuIeo4d1EdbdO8bKJxm/r5Aoxigsd0D6Wco4qNe58hIU/RRrnBZck7LgL1t8e4PjXcRNQWGVSrgZLmToURgxBaxKDdOIvyoVOvGTMhtHo8xQr6ssnBzZFHvCrQOSY7p/I5dasNrhUR87slJM3JoOO4IKnnLAJXkdLqkEq5P81CNyMGIbSIQVPgU4PY1MBWzOTSsAI47gCOJumPpWhYvnDkCFft4jKt63J+g+8+oxzrX59Tn8JxARYh5pwa1gVX0St7PpT4nmZTPSMGIbSIgfvii7Q26qTJrlm2GAv3fByNr2iyxsYruC6kMsSSQpFgrjKVVatmIYYc3PFlVqT/vKjnpVoNnAwlyn1yvvM8UOe0puJIgRGDEFrEQHmhVO2CEsSmauMxbqQy85n68VJ/L8dl4BdklWnrMaTQ5GK3PLKKJTtu/ECrBoFjNWQJCRxioGr52j3K8kC1BihznwIjBiE0iMFNsdQILlE+NGqr4JiPu2iyUNMIqwgylgnpGA05jxSq7o1UhlBmFZW4OPUH3PEpCh7dTDVtN5FGnzJqFwPKvVGfXRmMGITQIAZq1SIlHY36MWpqM6nHpQgwlpnk/nkpgj0kbJtGCC7ygtUx8QNuexLqWA0SSVGxH2OF5IFaiEqxBqjxxTIYMQihQQya9QuaH6x7zSJtRjtriTqROcV0vguJ45qimOUaLbpdcFNaJQQUuhblN/vHcAg8hdWg0dGUu84INW7CKQql/IYiUMcVwYhBCA1iyIRKUYMwqr+dSgyUD9Y1wYusmZSTmOoiKLo2JzVVo7VDjLamXdvAJaWQ9cCNH3AJtwix74vy+6jun6LzUV2iVALRtAY04gxGDEJoEANFmGvGFxBpWQ6u4NXIWuJOYI6wpxINJ6uj7LxuG4Jsoyxgn4FT2xC7cS0JKTlwnm/RuVNYDZpKE3X9aO1sIso46rmKYMQghJQYtra2sNVqYavVUtE8KOOorh/KB+YKR+pE17AW3HEcsuT08+f6zalCWLsQjWJtUM8VsgI0M49SWw1F3z01CE29P8pco7bboGr5HMtC0obbiEEIKTFQNQpK5oKmwEekffhuC4yi36CpEXHGUXPjQ8K6CDEZS6H7cbcyTZaCsv5LqdJSOc8uRuCXnZPbkK7ofFRFglt7RL1mEajjpG24jRiEkBIDNQuCoiloflwxJCMlI+0sE056qi/YOAKL4j4qshI0etuEEAqIU95VBskz4QjLqtyU1HHub6HG1qrKJqLKCyOGmiElBsqHOhqNpmZw0YumfDTaVoWmwI8hGaqw5/izObnuZc8HMb9moKpUVj8OQvmdiPw22rFuoiKLlPpNUF02VAudEo+I+a6LQBlHvX9pQz0jBiGqsBgywVLWuoLy0aQS+BqCnBpHSalFcsdyi7vKhE1K+FlcZVosIp8IY8lYOo5qRboEovHNVm0NaNcp5cGIQYgqYgxc7aXog9dclIf68Wn2ltEWKBxrwRfyZRp3bH1ASsSssyBxExWB+izdcVK3U2zWm+Sa1HGaGUzUOZcHIwYhJMRAzUiimIXUD5RCMnUU3KR0N1HHcQQKlVip95GHomByrOURQw6xDe401hCnni/mnReh6rgedX1nitCXtuA2YhBCQgxU1wnlJWuW8lNL9FNlZUhN5NjgtAaBIMpXeMsjAspGRUxaKnWs9joITdbyqYVuHGEOUNwGX3NJ0DwYMQghIQZq7rKWhkCt/tQ0aTVzx2MmPtUVxrkut3KaYinkBahjN4ow4C5tyumPFPMOqurtRfX5U87lWtfS2B5VmGu6ZvNgxCCEhBiyD5SzfkAeuBpJVRlJFAulDqsi9nycc1IFNEXQ++81FEgObWXgtgCJJdM8UNua1EEgmvE46rkoXgTqb5TEGYwYhJAQg6ZJWPWHrpn2Wsekj3VPcHo3UWtTQhuntQZiMblwg+RVtqqgxrO4WUcamW3Uhnqcc5VZ4ZRx1LkniTMYMQghIQbNIJJWFhHV3aQppFNoeZxrUl0AHAIpE+x5Gr+0tiHvvGXuSv84jfoQ6rNL9S1VFavinktD+eEQiMUYaoCEGDiCrkxgaAWoXS1Jmo1B1Qa5Foq2YKCmQGoFp/MqkjURQzqxJCgdpynMqdesw01EuS+qxUPJVJR0WTViECJ18Jk7GaSmJTVlTjO7iSs8qoqPUK0PThO5KhfzCZFDkSDh1Glovo86+glp3j9Vyag6u5BKMiEYMQiR2mLQzKXmrMGgsdCPZnZTHQI/RnhwBGBKUsgQqlkoAnVszLgmKQeaqa3U3kqUuUV9T5Q1XCRtMYwYhIglBs2UOc2AG7cATpo9oTlBm1AzwSGQqqqgQ3GHPMRaSXnwV87LQ6r3WwRN5UZL6FOtf4rCaMHnGhFLDNSPl/IBaBIDZQy1qKdqLVBT4McQFvWZlQmQvHPnbZQMJv8Y7eAy9dqUMZqWRRGqtqA144Xc5YC5AWgjBiFiiUEzYMXtpVTVxFuUc1HdeFq1EGVrKuRtZYuyUO/BD45Tzid1dWpaA6kUl6raVFDORVmfxb9/DowYhIglBs21YDWFpmZMQ6tYJ1XaYFUk47tTKBpq7FZ0bqrAp/6umFiN9JqUMZq9vmLWZshDHdlLsUVuRgxCxBIDNZVMW4uq6nrUcVzy0CRIqZvEdw9RzsURtGXHcMZy78Ud01R3YRE0z6UVP3DHSK0PamA5Ns7AlWttMKjg97///cyfIezu7kKn0wEAgF6vlzsuG5P9GcLbt29n/pSg1WrN/BnC7u4u6VyIOPNnGb755hvSOApOTk5Ex7vPcnV1lXRM0bjJZHLl//BnZSw4Pm9f0Xsp2peHx48fs4/Ju6b0XJqgzBsAgMvLy5k/Q8jea9H7HQwG079/+eWXuePOz89n/gzhf//3f2f+zMPHH38882cqtJA6g5k4Pj6Gp0+fwtraGhwfH8P29jasrKywx6bYV4Tz83O4e/cunJ2dwZ07d8i/98aNG3BxcQGdTgc+fPgQHPPpp5/Cy5cvAQCg3+/D3//+9+C4drsNiAitViv34+12uzCZTGB1dbVQIFLuy53oeZ/DvXv3ptdZWlqCf/7zn9Hnooxp6rnW19fh8PCQdS7K/ZUdv7y8DG/evCkdm3eNfr8Pr169Kh1HOdfu7i78+7//u8q53DGDwQAODg5UziX9FijzBuBnErq8vIR2uw0XFxfBMZS5Sp3PlOuFwJZrLHuEgfv370//Ph6PcXNzM2psin1FSBl81uzPTvVJasUFsuAYlJjNlHNRxmieK8ZFRBmjMa4IMdeaN9eOZoyIeu9VVz9r3jul3iGERsQYxuPxjGBGRFxZWWGPTbGvDCnTVauedIh6k0Dz3qtOx3XHSH9jzLliKlX9c2iTlWRMHeeqeplb6rkosT6tMYj07CUfjYgxPH/+HLrd7sz/dbtdePHiBWtsin2pQPFJUsZQx1H9qZmfO+Tv5ozRBDLjEJqQxiFizlXkyy4C9flQ4yHzjMxtUuQ+ocTKtEGJ9WmNAQD46aefZv5MhSTEcHp6Gvz/0EQqGptin4/379/D+fn5zBYDSoCJMoY6jhJAqwPtdnvmT0M+Wq3WdJNAk+zmGVSFo2olTvN6n3zyycyfqVDp7M0T2Nyxmvu+/vpruHv37nSLfeBZQKcosEMZQx2XBdLLAupVf7h1aG0GAwcUCzmTFWUyizJOU2mkZi9JkYQYVlZWrmgxJycnQSFWNDbFPh9ffPEFnJ2dTbfXr18zfukv0HTZNPFc1OtRTH6DoemgWh8UC7lqpVEDSYhhY2Mj+P9u3i9lbIp9Pm7dugV37tyZ2WKQ1SUU1SdUbZ4CmGunycD/X7cgjbWUxZmuC6hxN8rcoc4biiJUtaKngSTSYm1tbebfx8fHMBgMphr7ixcv4Pj4uHRsin2p8OOPP878GYJmjOHs7GzmzzxQXDtUklkE3L59W+1cN2/eVB3ng+qOa1qcqS5QtXxN4ao1vzSVRhWwcp4YGI/HuLOzg/v7+7izs4OTyWS6b3NzE4fDIWlsin1FiE1XzRaxKeq/Qu2OWHWXSMqYVGtCV3Hvmj2AslRbzrlippnfH4h6rXlrqKi5QpvmvWuuCEdJz9ZcQzqERtQxzDNiiYHa60RrQlEnQdUTqonCp+5zxZCDf6z0/Wmuk0Edxx1TlcJR9ZzQXOktdrGeRtQxXEdk7QbctgMhaJmV1OwfStaEpnnKPdf6+rr4ming18FwxyGz51HZuLyUVOo53TYeUlCfDRePHj0SHa9Z20PBgwcPpn8PxS8zaMUhAAD+8Y9/zPyZCkYMSqAKRK3UNWq6aiagQoIqg6bPlXKura2t6d81BVa/3xcdv7S0NP079VkUjQt9C0V1C3n7it6dC2qCgfv8Y+D+Zum5NKGZEUeZz0+ePAGAn4kor8cTwC8xpqJYE1V+fPTRRzN/poIRgxKogWVKd0RKShpVmFO1KAq0LAupZujC1dTKrLUy5DUG9OFnnuUFtE9OTmbIxoVb3FZEFkXP2r9unkD0O+PmPf8YS0DzXebB1cw1z1Wk5VPmFzV1lKKcUeUHpYuzCliOqmuA1Cu4aTW1owaDKderumFYzLk0G/dRfcvUcWW/NXYFtzJQx8eMk/rfUwX8m7RMKPVcmnPQgs81IfWaz5QmWJofCfcD15x40nO5AkG6oA9VUFFXRQstwlMG/9x5G2V1Lv8YauC2jiSEqoLKVTfaowpprYWBqMpgCEYMQsQSA/WlUbooUlNRtbQazYyVupftpAoh6jjOus/c5RZjEVouNA9UQqQ+6zqW/6TOB8o4Teu46vviKCE+jBiEiCUGRJoGoUkMmhO0iT3qqeNihJXGuNDY1OTAIYXQ/WmOk1oCdSsRlHNJrV6qwqiZ0hqCEYMQEmLIPoKij4my0IY7+aWTL1VhWtG5qnZfUcdRLSNfyFCKEasgB7+4jvoeuM9Z+j784rwqrqkZh6CuF6Kp5VPOFRtfQDRiEENCDBRrgPLRuYJJI2CVIuBdh4anTSAcN1ERQlp82TEc5AWvi9576Jg8aLub3BXzpN8cdZymAKaeS1PL57wnsxhqgIQYKC4gqpuI8tFRNQjKR5fKfSXNJsp+I+eaGtlEvmAtey55gWXuEow+8kiniNQQ44PT1CA/9RlrWrxF4yjKj3uuojmjGeOjzFFN11UejBiEkBCDZo+jqlPcsl5PZR8npzcMR9DkgeqWSOEm8oU9tW9OaOMs71mW3lpGCv59cCyLImi+r5SJBdJzUeYo9Xdquohi22EgGjGIISEG6oRotVrYarVU/fRFoJBMjF84DyndRNpZRxxhSBHKoWNCm2tJ+M8rbytTJvLORb3XovNTCcR1I2m/e+m9UZ8JZe5RFSnKNakCn9KoMw9GDEKkthgQ9fKtqYHlVNlEVWWluJo71WfNuW7RcwkJWgo5xBa15W2UvPXQNYuO88dT3X5VJhVwFSSpJUAlGU13E/VcRgw1InW6KiLtY8m0iE6no1r1KyUQd+JQhbSUQGLdSZxsoiKEXEScAqNYMuDEJ2Kyo2KJlHrOPMS6/fKQSvHRIDaKTOAQ1tLSUqE8yIMRgxASYtAMIvX7/emYfr8vOhc16KY5+amZVTETUSubiCNEQ4KXYjmEkBdQjg1Ux5CC7/LRqJyOqSfRLjQsgibJUDIQqRa9VmV0EYwYhJAQA6LexzcajaYfTJHvkfohV52mRx3HXaiFq2Vy8vzLrIAqaxaoCJFM2T1l1ihlfGxwWmolUsdRrQrtzB8uARbdl1bMsQhGDEJIiUEzC4HyYWmeizp5NIuJUo/jjo0hB2lKagzyYhhlVowkOK1BIFRBzk1TrfJbpM4TzS4HkuI2RCMGMaTEoJm3TCEZ6rk06xSakHVUZEX5QoUTgC0THIj5aalVWQ++G4hKaiFSKDrGtyyKEPOutb8J6v3lgasYlaUgU66pGcQughGDEFJi0FziU9Ma0Eyvo46jan2x4zgpqUUBO/+8FAugKOuoSEOUII+QKGmsiFefSdE36gf8q6yF4IzTFKyaAp/qCta8ZhGMGISQEoNmnxVNyyKm8Eh6TUR9zY9aCc2tXPaFLlW4FxW1lf1mCspqHCjFTqPR6EqfJW7lNHVs0XOjWpoxWWjUFFSNTgIUMnKtuqLkEUrvNPf+YxMejBiEkBIDxa+ImCbOkEKYU8ZQ/KxlZEk9H1VoIPLSVxH5lc5518nbWq1Waaoh9VxUK8F/XpTfxbHMfGFPTa0uehcx1q124WMRKMTgEnHRM6HIC+p9FcGIQQgpMYxGI1xeXlbLMuBaFhoEwtXWNLS/WAFTRsC+UCwTzCH/PVXr1y5qC22cQHfI2ij73vxjOM+3yHrh1EJotqRApCth3H5LGm01snFFRWvS+AKiEYMYUmJA1F1JStOHqtmz3h+nrdkV/V5OSqo/lqJpSwvaEPVJgutCCF2/TMhzM5Y4RYLU3xLTeoOaXFAkfDXb0yPy3E1llqQ0IwnRiEEMDWLgVhAXfYiue0Pj46HcmzvpiiZxTOohtd142bWpPXlC56VoXiHLIda/i3g1w6dsi2mUliGm3iLkcuIQrsZ3gpi2qE3aQQBRd25Trkn1QJTBiEEIDWKgfmTcrqdFGg83bbXoo+WY6ZpxC8T4FEluURcluBwSsJQ4QV2IrW1AvBpf4QSnOe9ew1rQbtFBHafdQYCy/js1ZlkGIwYhNIhBs7mdK6Q5/t48uJqwxsdNjVtQTX//2hxBwi1Oo1gOMYvk1IG87KiY5ntl1grHhZS6clqjJUWMhq8xrynxBWqWYxmMGITQIAZE3ThDplkAFJvD1IBX9rFppMpqpyByzonIE1KIV11EFJdNUcqoxL2kgTxCoGqYIRdXETiZXhx3U2wyg4Z7iPutaXRPpsYXJGswuDBiEEKLGDR9kW7qW1FOdMo0vyJQJz9Hk40dS9HkfWFKnXRFcYJ2u12pi6mofoJKVtwgOzdWwyEc6j1Qv/EU7iZu7y/Nuh1J4BnRiEEMLWLQfPFU7T2mMIhq0WgHobmaKWcsxX3iF8pxXENlRW2prIiy60p+Q1l7d0SeoOdkjnHeNTVLjxJL86+t7W4qAiW+QL0mBUYMQmgRA/WlUtddoGoO1GCVZnorYlwvJk2rgdspNZSFwxHqnHTUWKIoI4IYQgidNybWohUn8sdyLJYiaLeaoI7TjC9Qr0mBEYMQWsSASBPmruZa5CaKaRtM9b9qaEpUE5/jo5YII4pACrWKiBHk3FRUrY3rdw7dZwwplD2fVJYh9R7qylpylQ2N+IJGYVsGIwYhNImBEjiiCnxEfkodtaBHoy1HquAyR2iHtHiKCc7tOso5j/YWKyRi15KI6T5LfV/+N6NtLXC+Fw3liJp5R/0dWvEFRCMGMTSJgWrOaruJNIvdEOOC1VruA6lgooxHzBfqkniBBlFI89fzXF6x2Vgc4V327N3aibJWH1Rh796DRkM/RH7L7rL5XnV8AdGIQQxNYqAK8jpWl/LHUYUzNUjOFfhl4AjrkDCkCHhJgVgTIVlxToOQtZrv+ecuArV5nRucphJN0TjXJVz0u+uILyAaMYihSQyccnbKh6Bt+nLcSVTrIkVBEyKPdEL3wRHwkmKxJiDv/qm/P8a1xqlZQOQRPdUKoLZeQdRduAqRv1Z0lfEFRCMGMTSJAVHfTeRqOpQsJo7mlqIGgbNADtcS4AorjnBvclFbHjRqG0K/m/teUgantRMqtIiBMs6NqWi5l6kwYhBCmxi03URugIvaO0kjZY8jlN2x3I6eZeAKLWlgOfVCPFLkpd1q/dYySMlaKxaFGLfGCfU71nA3uZl7Zcpaq9USN85z0QhiGI/HOBwOcX9/H4fDIU4mk6ixsfuOjo5wOBzicDjEzc3Nwuv70CYGRF03kSsIioJXiLqLASHyyIZaYORfPyZ+QIF/DNdEz1tnuQ6SoAS1OVZNbOCdG2xGlL1r6liNZpP+vRaBOs8oq7Uh0t1cHDSCGO7fvz/9+3g8xs3NzaixsfuGw+HM392xZUhBDNrZRNRsJ+o4zmSh9m3iWALc+AEi33LQCixTi8403U2c7CYuOcXWNiDygs2h38EZqxXIjrm+RsorIm0+urU1mmuH104M4/H4iiBeWVlhj43dd3R0NHO98XiMAIDj8Zh0/3VZDIhxjes0xvn3SC1kK9N8OEKD40sOjadoq4iyTJ2y61O3vHchOSeXECTpuTHPXlLYqNWoDzFtZl8RqAH0bH5p996qnRj29vZwY2Nj5v/W1tbw6OiINTZ2HyLi/v7+9P+Pjo4QAMjupBTEwNEqUn64GhlKvk9ba4IjzqYaUtY94AqFvPviHk89X8pNkq0Sul9KryTEuOA01yLkjOW4LFO4kbi1SgD5Sg81lTUGtRPDcDgMCu1nz56xxsbu87Gzs3NlrIt3797h2dnZdHv9+rU6MSDqF52lSEml9m3iaPcSlwB3snMFfMh6kLY3RkxHFNIYhpQQuau8ZeBci1NBj0j3x3PmQYxFrrFONDWVNQaNJQZXi6eMjd3nYjKZ4NraWqG18NVXXwUnizYxxBTKlH0cMZOi6MN0r13m36QG+/yxZZM9JOhTuZXyrscRlhxQYgWpAth5/ZyoRBjzXkLHlYHacwsxPu2UqpyUEQg1JZz6+7UW5QkhGTHs7e3hzs5O7pZp63t7e0H/f0ibLxobu8/F9vZ2aWyhKosBkV9aT61pKJsU/X6fNIE4WUeI9EA0oqyYrSz7KnQMR8i6JnwVBFEVioiIQ0AxPahGoxF7uVCOshFriZbdQ8p23WXX1lqUJ4TaLYa8wHBIay8aG7svw3A4nJLCZDKpNcaQgfoxUYvYqB+xa56XBYw5k8iNCXAsDMq5YzKVJG6lvOOr0OZTIM9tRI0l5J2Hkluv4RLkuCeL3jGnqIyjlFGtAOo4RP2iNhe1EwPi1VRS1+1zdHQ0o8UXjY3dt7+/P7UeJpMJ7u3tke89JTEg8hfmoeZkl31MHBOdOum4wV9pW4vYjBmuMKxrIR4JRqPRjAWXwkqgCDfEWYWB8rw4CQ2IvICz+91rBIgRecqTdswwFo0ghvF4jDs7O7i/v487Ozsz2vrm5uZMnUHR2Jh9WXqqu+Wly4aQmhioH0DMEoZUzalsgnOqWTmaHvfcofNTJ02M+4NyDuk5tZEibpF3Tsp5uO8XcbZ4kKPVl1m//rk1XEOIvFiAdpFpLBpBDPOM1MTg+vs5aaFFiBH4HG297KPm+IYR5ZYDNVUzFHCN0fY5BWaprYnU9xKq8I4l4xRE4lojnKSHsrHUOcRxN6WIQ8TCiEGI1MTgBjrLgkwxH1bZxHLNZU4wkGOyUz5ujpaI2JyaBY5gjtXYY69D+abyILGwQsdSnq9UQeDEIapu6805J3WcBEYMQqQmBkRekIk6NqabatmH6Fo3ZcKbsxpXaDxFoHHXdXaRRxCxrqCiTKaqN4mWmfdcqO0YYgP+qZMRUilA2m1mOPEKCYwYhKiCGDgfg3axDSLddOUK+yoyj/xrLC0tkQPLqWsWqqyAlt5z0b1KrLGY7DFunElT0HMUpRjLnKrQpchEcmHEIEQVxIDIMx+pY6mBNv+cVAsjxuVTJrT9OECMC4IrKLXdSzHXomzcbCrp/XCsJ62kAMoz52QhcYopEXmWhTaBpM5EcmHEIERVxBAbCKbWK3ACxtRFhKgfsFvYFNPzKMZykAq21ARRN7SK3ULnomQHhY7jWhfctT1SWRZaBMIhMCmMGISoihgQ08QaUripuHnm/nhKplKMJhmqruUK9rIg77yTxGg0mtGKJYSAGCZkilWjUZNSdh2O9h/rQkrRZ6mKb8yIQYgqiYHzcXC0G87kp47luog4ml7eNaiCK6bS1keZy4cTy2gCNGIIZeejrjAWk8YqyViinD/WhUSdI02JLWQwYhCiSmJApJuT1PoHRN6i6CnrFbgprP41JMelIgiOQKwaZdZPbMwiNh01ZNFR+l5Jkxg459d0IWn3WdKEEYMQVRMDJ37gLg1YNsE5cYFU9QpcF1SGWMshdGyshpx3Li0NXAOcWgcttxH1XO5KZC4paMebuCQSG5zmuJC0+ixpwohBiKqJAZH+8fldUosmGScQ7Z6XQlCSiZsyeyXvWKkAz4tlaAvimN+U4h62trauCHUJsVOP5b5z7reIGL/mNKduQqvPkiaMGISogxg4QWBOSwDOeV0Lg+NSitHSYgXM8vKyqDMogNwN5D+rpm1SYSOt9fCPpy5RGWMlcl2VsYkZANXXN2jDiEGIOogBkW5e+uSgGYhONXFC4yXtFjTSUrU0+iqL2lIQQQaN5+Q/CyqRSzOWuFlvnCxAyv2kIBBtGDEIURcxxPr5NesVEHlEwv3QfQHQbrejMlRiBGLqqmcXqchC+16LrB+OdSYhb42MpbL7lLg/66hvSAEjBiHqIgZEupnpTmjuwilctw/HtxpbtUwVBtw+/z6KBOG81ytwoFm7ISFtbkEaYpx1wXE5cSxhThyvrthCBiMGIeokBu7Hw/FXcjQsThuMmMkdEtCxaakxMYMyjb6JqagaKPrdWs8xZXKBtEKe4nJyU07L7omT+VdXbCGDEYMQdRIDIs/cjLUEKBMkNt01O3+qgGPoeqkIoqn1ChykqOyWFhTGvnepdapZlMlRnuqMLWQwYhCibmJwPyLKx8/58FPWIMQWQUkyliQplf65ioSn5Nx1IVWrj9j3jKhnKaYiEmrKqT+WkwBSF4wYhKibGBB5HxLHEuCsCofIzzwajUYz5EOZYKHrcLSqvLiBpGaBko7KTZ1NCU6hm4TcNOMJnPcckwHXarWw0+mQY3Cc8amC0ylhxCBEE4gh+/BWV1ex0+moViNLMoliXD0Ut5V/TMwkkqa1Us9Zt0XhZ5lVQQYhy4zTOyo2Gy10LDd43O/3Vccjxi0DWldsIYMRgxBNIIYMnCwljkYlMbGpQiamT5LEcggdnwkhasyDe+68TdLN1L+Ov2pdVSRVVO1dRTwhREjcWBfFqnPfVdn47J6y+yp7xpx1JFLDiEGIJhFDbHCZY95SiCSmAjmmXYH/OwB+Duxx3TZFQryphW1U7T/vGE1rJe+3SSvPU7fU4L5nzpzxv2eOJULtLpwSRgxCNIkYEHnxBk7wLJsUVI3an+QUk1tTMHCFuvv7tLXpvHvONMkYIZ+3hSyGlEV5WtfTyjijkhHXl89VjlwLmGuJNCEWZcQgRNOIIdYSoGgpfvM8iv+WE6Tz7ykTmpTjNIPKIQ1Yy8VERVlQO5X2X4Yiy4e7BoU0U0yDTKi+fM7z5sTMYuZIFTBiEKJpxIDIS41zNeQyIRzKIuLeT0xKKldYcH3NIZR1SG3KBK4C2fsIERU14SHvnLHPNca6DH2/KRIq3PNzLAWKVV0VjBiEaCIxcCuRM22FqkFxJ4okoBgrjCXkEjpPyMXUbrdxeXm5UkuiSpQF0GPJMfRuOC4g3/VGtSj94ygWr38M5f5iLIsmpTMjGjGI0URiQJytRKYIrtj2GrFWQGxLC25fnpD/XrNYS+vcTYD7+0Kxg+ybktR8+Bo7x/0kqYvgzgdEjCKSbCy37UzTYMQgRFOJwdeQuPnW3JWzYgrTqIiNO2QICXSJll9kRcwbUVAL9LQL3WJIXsN6LLOgQ8dQLRLO/TWlkC0PRgxCNJUYEH/R0DKC0O6P5I6nLrDCLcYLXUtDsGhNyrwAqrtx1sNIiaJMIn/Lxkifj5SUQ8+Xk84ZQyjS9GzKNepuklcGIwYhmkwMGTgfOTfNLiaYhxhXBIcY153VRcidQf2tMfdXJHS1hXCGUApn2T1lBKZBXEUuNymZxy7mQ712zPfFTTWN6Q5QNYwYhJgHYkhZmBMStDETUOIainEL5QkvbsplEfzKV+rmE4cv2POIJfb8msJJyzKLrU8I3QMn5Zn7XLjp4dxgdl0wYhBiHogBkWe6atQfUCBxD2lkHRW5gVJkieS5cmLbWFA3l1hSFevlXZtLOqGYB+cc3C6/eb8hxTfPDWbXCSMGIeaFGFwhSBV6sZlHnGKw2PRDxPi89KJ719B2Y+AKGanFUIV7YmtrC9vtdjAAH1PXgKhD9px1QfKuS1UIYovemuo+cmHEIMS8EAPiVY1Fu0cSIr86OoPEdaBVs1BkQVAD+NcBKUhU0lEV8aqCQc0+irEUYhIouPOobhgxCDFPxOALPk4xW6vVIk1WiRYf8g1z/MqhmoVYl1CRe0Q7WNx0lKW0coW4i9B7k8QjMhdoTEwhRc8k/5h5+WaMGISYJ2LIwP1QfSuAG5TjpOT55EVdnyF0XQ33ip/yGxKKEsHYVPjvL/TbJYH6vJYjHMGpqYSkzKZrempqCEYMQswjMSDyTFtfY6R+4Jk2GKNhx6zP4F5XaxlP/7zu7wlt1N77TQSlsltihRVdQ6PXUqzbMlX9TWzNThNgxCDEvBJDzEcbaxLHCADfzI8Vtr5mmPU30liIJy/46gpQqcWSEu6zKfodvV5PRbBppbL650md6CB1H1GPaRIaQQzj8RiHwyHu7+/jcDjEyWQSNTZ2n4udnZ3C6/uYV2LIwOnEijgr5LmaE1cYhI7TynbR1OxdC4VSs+BaHCkJYzQaYb/fx16vFySrok3rvvJcRjFuKGmSQWwsQ+I+mlfrsRHEcP/+/enfx+Mxbm5uRo2N3Zfh6OgIAeBaEQOnEyvi1cnJdStxTfjQNbWsh1QTt6hNddHma+0uybh/d4mF8nfqtTW7xGbPwH/nWqmskn5ZnDgEt55nnt1HLmonhvF4PCO0ERFXVlbYY2P3udjf38e1tbVrRQyIs7nf1EBqrFuJ24/JvV5svUPoPCHBmXKhlLLOpSk3l3BSFrtlFkpeHCHmeho9smKqjf3rUtdK4F6nqaidGPb29nBjY2Pm/9bW1vDo6Ig1NnZfhv39/en/Xzdi8DuxUq2AWA3ODSxzM3q0rIe8c8XcUyz8oH4Ki0GzxUcR8p5lr9fDfr8fdQ/Sdy053j0m1XKhTQZXrrVBGaenp8H/Pzk5YY2N3Zedd2VlpeROf8b79+/h/Px8Zpt3PHz4EP77v/8blpeXAQDg8vISHjx4UHrcYDCY/v3i4gKePHlCut5f//pX6HQ602s9fvyYdD0AgIODA9ja2poeDwBweHgI/X4fbty4QT5Pdq7RaDT93Rmye+p2u9BqtWB9fZ18Tg4ePnwIb968AfxZ4YKLi4vp3xERLi8vg38/OTkh//2f//wnPHz4MMn9r6+vQ6vVglarBYeHhzP7Op0ObG1twcuXL+Hvf/876x4ePHgAN27cmDnnYDAARISDg4Oo47e2tkjH7+7uwr1796bf12AwgDdv3hTef/YcHj9+DAA//3bKfS4UVGkJEYfDYVCbzzR46tjYfYg/WyLu/xdZDF999VVQM5pni8FFTIAuRlPSWKM5Rc1CXlzguhW15YGS0iqxtlK0V49NY41xH81zXMFFMlfS3t4e7uzs5G7Pnj2bjgv5/7P9/jnzxsbue/bs2QwRlBHDu3fv8OzsbLq9fv16oYgh1icbOzmkla/Z8ZortZXVK1wnkvADsHlkIHkWeV1oYwSs646TfIOcFFiNZ9A01B5jyAsMh4Rz0djYfRlxZBsA4M7OTjDGEcIixBhC4GZxhDQ9jq9emnWSafsawsUFlSS0s3rqRF4XWE3LAPGXWhANosneU3Y+zmI+sVYyd630eULtxIB4NZXUdfscHR3heDwmjY3d5wIAZq5XhkUlBsRfJgyn/0yo3xEVWoFlzQC1C0rlcyaU5kWD9Ft+5P02TspmGfLcUTHnDwXwOcHu2K6nLqEtivvIRSOIYTwe487ODu7v718pMNvc3MThcEgaG7sPEXEymeBwOEQAwO3t7WtvMSDGd6x0J2vMZA+RS8w5fE1Qu6iNQhKuRdGEDq2uICyrs0iRi69ZiS4pgIx1mfpExLFM5gmNIIZ5xiITQ4bY1g6S2EFo0muu1KYtoN0KX0phm+uCaLVawVRTyv25tQM+UYVIKW9zLQZNMihqQqiZZkx1W0qSHvy0bgB6kea8wYhBiOtADG7fIq5bSVqU5muYrVaLPRnzgpsx98O5piuoKcK5aHOFWXZOaaFcapdXHilLiEfqJozpFJxBUn8zbzBiEOI6EANivFsJcXYyc33AiD8Tk+vTXVpaUhcuVbYwcHsrFVkM3I1iMaR0YxWl+0qvr7Fug9v4MCYxIvsGua3g5xFGDEJcF2LIkE38LEWU69ONyRpBvEoOUu0zJGjqIIk8hHzZIYshs4Tq1l5Tuew0WqjH9vdCvNrlN7aKe95gxCDEdSOGzKftTtbYbA7uBM+uffPmTTXhQ0nNnIfsoqrhJxikeG6hYrWqExkk64LMM4wYhLhuxJDBn3CSJThjUhRDzdqkWr7fGdO3Iq6bcPBRFYmGrA9u5pJWn6Xr5D5yYcQgxHUlBsT45RER09YspBRQmbBY1CU9ffgCMrRpuLLyEgS4ZB9SPKSV+NfFfeTCiEGI60wMiLKV1qSTuOxcVdUrZC6VeSpsC8FNuQ1ZTa71pBWD0ax+RpS5KhFlKdaLBCMGIa47MSDGFwtlkLbTKDuX5uR2M1vyiCIjOG0hqgWf7Ch1F6kW8tGyPvwgNTe5IWZRnkWGEYMQRgy/QLp6lUaw0b2PlJquex3fYigTsFWQhi/8V1dXSYSW6ln59+U/E8k63KFYE+dc/ndH7aq6yDBiEMKI4Sok691Ke9+4KEqhlAiisvv3ew9RC9uy9NNMcw6RTohY/PRVyrVC10i1qE/2HqRppz5CJMN5r34jvJTfxbzBiEEII4arkPTEd88hqXtwUZRNU5XLwA3ixpBG7BayGFqtViUCsIiYtddsiPnGFr0RngRGDEIYMYQRmrwx0J68eQVTdcYCXOLSshjq0HzLAvWtVkuUyaWVrKCpdCwqjBiEMGIohkbcIbRIjEbxVJG2XtVayYuAolYYWu9Kww3lKysSN+Wiw4hBCCMGGjT63GgHLRGLi9q0BNsiosg9p1kMqFX9LK2PuG4wYhDCiIEGLTdA6jbaRR1Br3P1M4VANQv+NKqfEXXiXdcRRgxCGDHwoFWzEJrwmhPfLfbK2zL3SRNrFaSgEEGK365Z9OgXX1rGER1GDEIYMcRBw0VQtAiMpmZIXa0NAGoN/sbCTSdttVqFhJiKCPKerdZiPotE2lXAiEEII4Z4aNYsIIatkRSavF/YVlQ74GYM1eWOyoRvr9djVzynsoaKXIKS64Usjnkh6CbBiEEIIwY5tNMHUzbWC8Ht9sopMPNTTUNC2yWfTJvP2o77i/yE0lc5dRLuNTTbX7jII28pYVb9zhcdRgxCGDHoIUXNQshfXYXWHmpJQXVHpdqyGgiffKoQnqkEt2ZMwvALjBiEMGLQQ6qahVRuCwnKiEPTYhgMBpXn62fxn9CiSinfqVkIOjBiEMKIQR+ptMCyQOciZhdVjSqK3cxCSA8jBiGMGNIhtVZYVv1s2mc5MiLIrBm/e2tmMWgVIto7qgZGDEIYMaRHngDXEgpFOfvXvbDNRxkRZG4vrWI3sxDqgRGDEEYM1aEqrbEoJtFut68VUWTtJG7evInLy8u5RKCdzWQWQr0wYhDCiKF6VCU0qIVtddYoaMMngtDvTZnWahZCM2DEIIQRQ31IlRMfQuZCcS2Goi2UHdQk4eYX6ZUt8pMRRZX1DfNOsvMMIwYhjBjqR13pqH5sglqjkJdWKiG0vHuJLXpLTQT+fZuF0CwYMQhhxNAc5Ll+MtdHVf2L8moUKKThE0jm1skCvX5bC+kqcL7FkPoZ+ZaKWQjNhBGDEEYMzUWoYK7dbifXgvNA1eolG9ViqFojt9Tg+YIRgxBGDM1HFlB1CSLrINpk/z/HYmjS/btN+zIrxI9fVNmOw8AHV661EBHBMMX5+TncvXsXzs7O4M6dO3XfjqEAu7u78OWXXwIAwMnJSXDMYDCAg4ODKm9rYbC+vg6Hh4fBfcvLy/D+/Xv4/PPP4dGjRxXfmYELrlxrV3BPBkMSPHz4EL7//nv4/vvvYTAYAADA6uoqdDqd6ZjDw0O4ceMGdLtdaLVasL6+XtftNh7r6+vQarWg2+3CjRs3Zkih0+lAr9eDVqsFy8vL8M0338CHDx+MFBYURgyGhcDBwQEgIpycnMCHDx+mRAEAcHFxAZPJBAB+JopWqwX9fh9u374NnU4HHjx4UNdt14bd3d3p7+/3+zNEMJlM4OLiYjp2MBjAhw8f4OXLl3B5eQlv3ryBhw8f1nXrhgpgxGBYSGREsbW1BZ1OB1ZXV2f2v3r1Ct6+fQuXl5fw+PHjGU15EYnCtwb+/Oc/T3//q1evZoggs7q2trYAEc0Vdw1hxGBYaDx69Ag+fPgAJycngIhTS6LX68Hy8jK0279MgUxT9omi3+9Du92G27dvw+7ubl0/hYSMAG7fvj3zp28NvHv3bvr7e73eDBFkVpe5ia4vLPjswYLP1w9ZkHV1dRXOz89ntGcfWdD1zp07MJlMpsdk/04Z7PbvM7tmr9eDH3/8EQDyg/AZsmMtaHy9wJVrSYjh+PgYnj59Cmtra3B8fAzb29uwsrLCHhu7DwDg+fPncHx8DGtrawAAsLGxQbp3IwaDL4A//vhjeP36NfzqV7+C//u//4OyKdPtdgEA4Ne//jW8evUKlpeX4e3bt9M/V1dX4fT0FH71q1/B73//e/if//mfXKLhCH0AgKWlJfjpp5+uXNOys6432HItQcos3r9/f/r38XiMm5ubUWNj9z179gy3t7en+9bW1sj3bnUMhiKUrdSW1VOk2PxajeyavV4Pu91ubYV+huaj9gK38Xg8I7QREVdWVthjY/chIq6treFkMpm5DhVGDAYJsiUwu90u9nq9meKvULuKwWBQSDQm9A0a4Mq1G0IL5QqeP38+NaUzdLtdePHiBdy/f5889vDwMGrfysoKnJycwMrKCrx48QLW1tam7iSDITUePnxoqZyGuYd6VtLp6Wnw/0P+0aKxsftevHgB3W53Gn/49ttv4enTp7n3+/79ezg/P5/ZDAaD4TpD3WLIQ54w544t23dycgLHx8ewsbEBKysrsL29Daurq7kBw6+//hr+8z//k3xvBoPBsOggE8O3334L4/E4d/9nn302Fca+dZC5dnwUjZXsy7bsGgAQdGUBAHzxxRfwH//xH9N/n5+fwyeffJL7Ow0Gg2HRQSaG7e1t0riNjQ3Y29u78v9uiwLK2LW1tah9lJQ+F7du3YJbt26xjjEYDIZFhroryQ/0Hh8fw2AwmNHcV1ZWgkFhd6xvYXD2DQYDOD09hZWVlWktQ8haMBgMBsNVJIkx7O/vw1/+8hdYX1+Hg4MD2N/fn+77+uuvYX19HXZ2dkrHSvf94Q9/gKOjI3j27FmKn2kwGAwLCWuJ4cEqnw0Gw6LB1mMwGAwGgwhGDAaDwWCYgRGDwWAwGGZQWYHbvCALuVgFtMFgWBRk8owaUjZi8PDDDz8AAFiRm8FgWDj88MMPcPfu3dJxlpXk4fLyEr777jv4zW9+A61Wi3xcVjH9+vVry2ZyYM8lH/ZswrDnko/YZ4OI8MMPP8Bvf/vbmVUL82AWg4d2uw2/+93voo+/c+eOfcwB2HPJhz2bMOy55CPm2VAshQwWfDYYDAbDDIwYDAaDwTADIwYl3Lp1C7766itryOfBnks+7NmEYc8lH1U9Gws+GwwGg2EGZjEYDAaDYQZGDAaDwWCYgaWrEnB8fDxdQ/r4+Bi2t7eDK9JRxr548QL+7d/+DY6OjqKv0RRoPZeifS9evAAAgPv378Px8TGcnp42bm2NKp7DPH4fAPaN5KHxMgUNpbh///707+PxGDc3N6PG7u/v49HREYYeO+caTYHWcynat729jQCAAIAbGxs4mUyU7l4PVTyHefw+EO0byUPTZYoRQwnG4/HMA0ZEXFlZEY31XyLnGk2B1nMpO8/e3h5OJpPGTvYqnsM8fh+I9o3kYR5kisUYSvD8+XPodrsz/9ftdqfma+xYjePqhNZzoZwntJxrU1DFc5jH7wPAvpE8zINMsRhDCU5PT4P/f3JyIhqrcVyd0HouZec5PT2Fp0+fAgDAwcEB/PGPf7yyVnidqOI5zOP3AWDfSB7mQaYYMUQi78FLx2ocVye0nku2zw2Yra2twWeffQbj8Tj+BiuC9nOQXqNJsG8kjCbJlGtLDN9++23hx/PZZ5/BxsYGrKysXGHZk5OToNnKGatxXApU/VzKznN8fDzNMMmyK46PjxujEVbxHJr0fXBg30gYcyFTyNGIa4q8IE4o0EUd6z92zjWaAq3nUrTv6OhoJmA2mUwQABr1XKp4DvP4fSDaN5KHeZApFnwuga91HB8fw2AwmMmhPj4+Jo114Zp0nOOaAq3nUrZvOBxO9z1//hw2Nzcb9Vyqeg5F12gq7BsJYy5kCok+rjnG4zHu7Ozg/v4+7uzszLDu5uYmDodD0thnz57hzs4OAsB0DOW4pkLruRTtOzo6wuFwiHt7e7izs1PBr+Kjiucwj98Hon0jeWi6TLEmegaDwWCYgbmSDAaDwTADIwaDwWAwzMCIwWAwGAwzMGIwGAwGwwyMGAwGg8EwAyMGg8FgMMzAiMFgMBgMMzBiMBgMBsMMjBgMBoPBMAMjBoPBYDDMwIjBYDAYDDMwYjAYDAbDDP4fVXtD1L0XGhQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot trajectories\n",
    "\n",
    "plt.figure(figsize=(4, 4))\n",
    "plt.scatter(qx.cpu().numpy(), px.cpu().numpy(), s=1, color='black')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "9ee4e13a-fc88-4fcf-8181-1b563b5cfe83",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArkAAADECAYAAACbfTmCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAslUlEQVR4nO2dz24UybbuV9mNORcLUy7pTFptLNV+A7B1H2C734AegM4Ub5juARajPbSaB7jA7vEF7W2/AfUCR8b1Bq4tNxKzYxtz4Rwa4bwDU+XMVSsyMzIjIuPP95NKwiYrMx0ZK+LLFSvW6mVZlhEAAAAAAAARsdD1DQAAAAAAAGAaiFwAAAAAABAdELkAAAAAACA6IHIBAAAAAEB0QOQCAAAAAIDogMgFAAAAAADRAZELAAAAAACiAyIXAAAAAABExw9d34AtLi4u6P3793Tz5k3q9Xpd3w4AAAAAAGBkWUYfP36kH3/8kRYWzPpeoxW579+/p7W1ta5vAwAAAAAAVPDu3Tv66aefjJ4zWpF78+ZNIrpstJWVlY7vBgAAAAAAcM7Pz2ltbW2m20wSrcidhiisrKxA5AIAAAAAeIyN0FJsPAMAAAAAANEBkQsAAAAAAKIDIhcAAAAAAERHtDG5AAAAAABW4XGkWdbNfQAReHIBACAler35DwBAH8l2Fhfd3wdQApELwoBPyg8edH1HAAAAQJGLi67vAOSAyAX+I70tv37t/j6ILt/S4QUDoaLqr+jHAOhR5rGFPXkDRC4Ily6WhaS39PV19/cBgGk2N91fEy+MIFTgsQ0CiFzgNy9eqP/P9SCjmoR//93tfQBggtXV4s9v37q9vmRPELogRG7fJrp2rfg7132ZvzAuL7u9vqdA5AK/efy4/P9dDSRlYhuAEOC2cnIyf8xg4OZeuvAaA2AKbkvHx0R//NHNvRDJq4mfP7u/Dw+ByAVh0VV6Fi62uRcMO2pBiHB7Oj11c90yr/HSkpt7AMA0fF5w5YRRrSbCmwuRCzyGe3umE7LrZSFp0uVeMMRnAZ/hKxFlL4sLlqeFKnv9+tXu9QEwyf37V/+WVkdsU2ZP8OZC5AKPUXl7pGUhmynF+KQ7FQi3b9u7JgAmKQv74YLX9WpJls17wADwFS4qX70q/ryxUX68bWy/pAYGWgOEAR8o8m/PRPZSipUNUMfHxZ8xuIBYsDUxq87LPWCwJRAqBwfzv7PlhOH21OsRfftW/F3itpT2Xw/8hW+A4eEA/O2ZyE0sX5mXC+UcgY+own6qfmcaaVOb6rqwJeAjPOzn+XP5ON5/bThhJPEqhc0lbksQucBP6myA4cZrOpaPvyXzWGCieY8yAL5RNzWY7Vh3btN8WZf/DIBv8LCfR4/qf9d0PvWyMCOE0s3oZVmcMv/8/Jxu3bpFHz58oJWVla5vB+iSn2CvXVOnZ5EmYhNdejCYn5RV5+X3EKdJgVDJ989er3yTpK2+XNdOYUvAZ3T7Z5f2FJAt2dRr8OQC/+BpT8ryD9oyXC5wVctSAPhMVdgPx9WmGY8nXABE6oT9VGEipE4n7AdA5AIP0U17wmOT2k7M0vfLlqUggIGv6Oa9lTbNtF1m5fZUthGG2xKKsABfaFIR0EZIHbdpVWgCbImIEK4AfCQ/Kd64QfTpk953iC4NXCdequxcdUwkoKUhkBB1w37KvkfUvD8vLalT8Nm+NgAm0Qn7UX1P97tV5yIqt49AbAnhCiAdeOWwOgKXaP5ttqocsAppQAIgRPjSaJuyo02XWbnAxaoHCBG+mqEjUk3loZbSkHkqWn0CIhf4RdM3XJ6zlki/pCGPudK5Hz7Y2CxOAUAd2iyNmlhm1Q37UV1bsksAXKIqm1sXE84TnoZMyvbDgS1B5AKP0U2Dwg1aN7aXx1y1SWlkqzgFAE0wUVFMZ2KW4v+aep2axEICYIs64pLTtuy7ZHtNVmYStCWIXOAP3JAl72zbc6qQNsNIm3AACAEe9sMritWhzVIoDxfSFQYIEwK+oJPtpwyTeah1bDNxW4LIBXHRdGI2ETfFv2M6+TcAdWnrOZrSJHMJF9hE+sKA37+UNgkAF+iuCKpoKo7bilRuS7phfIEDkQv8pE0lMd34J1tvum3juAAwQZvqR9++6X+HT6omNsfopkIDwAY3brT7Pg8bqpp7pBfGtvZkSrQHAkQu8ANu7K9eNT9XWy9Wm0GkLAcoAC4wHfajMzGbfGFsEvsIgEmaZvtRoRs2xOeyptlJ2orzgMGMDOKEe69Uky//fdvBgHu+ElsaAhEiTcxS9hBp53abF0a+vGuiWhQAOpgK+8nDVynrzk1EzXO/c3EueYgjBSIXdA/fiW1ieVPyXvG4PkmAtn1T5yS2NAQ8o03YTx7uQZKyh/Cd2yYyOuQxUS0KgKa0CfvJI61S8pdG26V7bYh3T7FW8WwymdD+/j4Nh0OaTCa0vb1N/X6/0bHj8ZgePnxIh4eHta+PimcBYbMqS9m5bV13ebkobpGwG7jCpS3lqxHqVmKqy2BQjMeFLQFXLCwU+5vpvqczN92+3T7saH29uE/EI1uyqdd+MHq2HL/88stMlE4mE3r48CHt7e1pHzsVv+Px2NatgpTo9S6N22ZalU+fiudfXGy2gQcAn+B2U7ZKYWoCPTkpXnNhISkvFOgQ1yJwMJjv71NMpNM8Pp4vTeyR0LWFFZE7mUwKPw+HQxqNRo2OvXfvnvkbBP7A4/hMG11dQWvT2DEpAxfwJU8XExj3dtkmgUkZeIipsJ88fG46PTVbRAUQkaWY3NFoRAMWUzIYDERvrM6xIEK6qMDCRa+NjAim4rcAqIuLKnt18kmbnpTbVB4EoAkms/3owIuomJ6bbIh1z7Eics/OzsTfnwi7dHWOLePLly90fn5e+LhkcfHSLqYf0ABbDVc16doIJeDLS+gUIBbK+rINQcorD8KWQCy4npu4WE/AlpxmV1AJ2rbHEhHt7u7SrVu3Zp+1tTW9m2sJX5FOoO+0h+8gtbmsv7pKfDjJiJJ8swURYjvsJ0+ZnaIUNggdHjLQNDdtWxCmYAQrIrff7895Yk9OTsTsCjrHlvH06VP68OHD7PPu3Tvd224FagA0wGEVowGdUEY09+m9trgMBQENXOE67Efq2zYnZX5uKU8vACbgIQNNc9PWRbIbm14yfj0pv3VEWJFmW1tb4u83hKUsnWPLuH79Oq2srBQ+LpFWFeDN1cBydaPTU6JFymgh91n87tuVYv2NkODSkAl46I+15xMrLvrZq1eXk2X+4xIX8ccAuILbksvNyl3si3GIFZE7HA4LP08mE9rY2Jh5Z8fj8SyrQtWxeXRDGFwjjfPr6+7vIwh4IQZe3cggVV52/uIOumN9fX58x/OpgA8yyOYBQDNchv10SUIOF2uL7Ht7e7Szs0P7+/v08uXLQo7c3d1d2t/fr3XsaDSinZ0d8Xs+wlfx8rmXQQ6HlcDqbAi3ZvNYZtVCZS8Jjcn6pDLIcFuCBwGYJnKv5gz+IixVWIsEaxXPuqbLime2iv9ERb6R8pWTLF5myvRZ2CwOVXoT6AwiS0vVlVvRdAK8WELMhUdgS8AmvFhCzKsiHtmSTb2G7VIWwLhbweJi8WdLAlei7NnAW9gtXOBKm5rxjBg87MehwM3HTTt7LtjhC2yRWtiP5X0wvoARwxGYnHM4Gjx4m9+4UfzZWdgCllkrkdr90SNZ00S8sqaPw7CfPNLz6vUcRONwEc9FPgBNSSXsZwrfBxOpLUHkWgLe3Jo4rAwmOYw7eU6pDaYNmD4XyTHpMPNcWPC3uA54/drxO1xHIh9EToorBpHaUoJP0h2rq8Wf4c2l+UbglcEsXUYnI52V55TioFkT3t51QsVgS9RZ2E9V2//+u2WnkAdiHkRGh2E/nZKALWHmtYhUmRib6+0j5VUtK8QkiSjj+bGxzFobKZoFQlfAk5hB6dl8/ny5kdAKXMxzsQ+ALpF6MStJwJYgci3DN88gh3kOSxXBeF7VOlUZedSE9UwyqQ6qDC5Uy/ZCSN54vDR+x1HYD39e08tKQvfrV0fvcp6IfRAJCXg3lURoS0gh5gBpYrC0Su83DlKWPHgw/yJR9zL89oxnY1peLorbOE1PC90ugfR831lYKP7hDhphc3P+5Y9fVno+Vsa79fVibHuSnaAa3k2m3L8/X5AxWRYXi+Iutb7kgS3Z1GsQuY7wKCVddzhohLaXsH6LKeU0rYC39eqqHOJT9T2iBO2pgwGFX/L588sMGFXHEVm6PQyqlVSF9KDJCP2IqPM2QJ7cCEkuJJOvK1vy4prGatxnhEtDbagjcIkcxVCHhM6uyoZIadskgUuE+GlfqNPmTnMch4DDbD/ADRC5juADf3IhmQ6CkZuGKbT9jhYYRIlIHdtZF56sIpVqnEQ033hluyoNwdO2VdmJE6FrKaY/BqQ0biiEI+Ao24/38I0r0u7tQIHI7ZDkvLkWMWmTfN+B0QmAD6LJzi5FdOcWKcoDTWmHplkSrAtdHlSKDjCDp+KeapgsU2cz7PUSXxFJGb4sw3dvBwxErkOS9ebyWdKCu5TbZJtLSOlGsYvfHG29uFOkZxyRA0KG21Kd1CEt4eWWdWxLerbQonZRVQ+c8u2b+hm+fZvQ82mTTD1CLoho2i0yiieUBSK3Y5Lw5vJZMgD4SqjRaAt+8sTdJ21WCLlXKiIHhAy3JVVgrCHaps08PpYnSmM5dLlaw9voHCpBW+XVTQ4HYT9lTEUl/7i65iIVO8p/UhzzErIrdEBSmzmXlooTs4WMArba02rqt6Q6wRU22lSqmBblnj7rOe6qL9m0m1rNuJCoLUk0bWeVmIq2KTuwpTJ0xGzbZ1J2rQvq0fS/MyJaoMxJH0B2hciJ2pvLPU8Bpczixs3j3Iwi7RRJABMvDfw5RTkxSzOTZVviXr42XiVkXOiGuragOi7KZ9SBLZWh28ZTz6uU8aTO98q4oN7Mn9ujOMZSiNwOSCY218EIKeXutIm1ZVarCtoPpLy4rq4VNNIfo1pnNgjvom2949aEbp2KFAnQ1r6iF7oPHniXZFuaT+qGBp+e1hOudY7Z2LhshsXsggqHRvDwIXI9QfetzHs6GkxMhyjyW7YaXhzBgKJD3by4dZC6VhTOcVWfcOzFNaWpnXl0oxtQy5H6ehP7UsXp9nqBhzuvr8sbKxxs3CyDzye3b1+GBmdZ8VOFJGTriNvp+TsOR7YKRG5H8I7L81AGjaq2p+XLuNoca8yJltAaLt/AZKM8PO9iwTvHOwyU5Jcwqakl+2ktoLhYiWpArYb39TZd5Ns3ebh+/drgSpZLBgN5MMgy6xs3y5DMWxW+VVfw1t2wVnquyLJMQOR6RLSep17PSZJtW2+jVmM+VW6TyOBL3VKatrZIXSzYplS9STkQuLa8uFMkwdw6e0nd+sIR0jYDhsTxsezk/PpVnQVA9+OEpSX5hafjYFNpoUEnfrrp7df67sHB5cNfX+/c020CZFfomKg2Bjv0PC0sFE+7ump2+bvqekQWd4cT2f+DHLK8XIw7t72ROXibWlyUA2Ad/SGu2s/KcOFZzKULbD8vm4LU6qPp2I7KkOKnmw73dZ6PB39yKciukBDBpkx1vLTKT2tbD1pNSSW1UUTLrXxjpeuNzEE59DqemG17cfOUbXTq9Rp6KKVSv0F1AD1c/Gk2u561+1flEfRA7Ul/c5v5q8w728brGwsQuR3DO+Dbt93cRytUgVqWrIunXLt2zcpl5jCZUmmOSONzuwjBCbYpl5Y6n5htxuJKlK2GXlxcCd7aMbuvXsnKPMhgUn1sdRW+EarpR8L4GOFx0l+pH5u6rbrtnRoQuR4SVGnS9XU55YBFC+OewT/+sHapAtKEb3THcbDqTI3JDTE6qJrS2x3iy8vO7YjDPacOspTRo0f1/sTXrzViOSVDDbDqYhVSPQPfkZ610Q2iHgtcovnYc1cOmpQJwCzih9tfMKVJX7xQ71q1RNfhHDwrgNFyv0RRCd2uX9akCWQqlrxifV1Olu14YuZOZJdhJW3TJYknlL4YMaHU2bH2aDwXuNLtuXLQpAxErqd463HKI6lxywMKD+dwPX5JWQGMz51S/p4AJ2jePVw/qz/+UHu3vGnOzU3nL4oSfDW/y/bRTZdUeiLpSxHA/wwbKflsYvzReC5wJTy+taiAyPUE3uGNewhN08Eu5q49g1Osz53Hx+pUbECLb988ruS0uSkH4Xcw+/HVfKsbLTXQEbvi+BCx0M1jIyWfbaR3+UYrdQEI3Ai7XDBA5HqMjfyHRugoTU/XnsE81tPbqlQG33XnKb6l8SoTup3sSXrwwBuBG0KXmordslK1jx8rbFByc4YQwKqgqyI4ppHyWmtvvA5A4Er25dHtRU+4lh4hvOP74k0pIA0qESSM1kUV/2a0mqg0En7+HEgsi3+oJpZpgntnvHghL9V0NPPxcGCfJ+CTk8v7Kxty5p6l5ObMsu4D/BsgeatDLsnaytEegMAlmrcvKcsdsAeKQXiG17nMVeV6HVQz880zOMXJ8/K6U8zDb/f5806rZ87x4IE6HMhJs3r0PKWICY+7lkhtreNRuzfFd9tqivajCUTgRtDlnGBTr1kTuZPJhPb392k4HNJkMqHt7W3q9/vax+qcJ0+oIpdo3jCuXfNgF6ZkrY5uTBIlPg0UELpFfH0h4XQyT3r2HEN5VlVIFQmJ4hK6MbyQqOBVEYlKBHwgApcoHvuyTZAi9+7du3R4eEhEl0J1Z2eH9vb2tI/VOU+ekEXu+np3+UVFOh5UQhgoIHQv4bfoyNHfGFXXtuIh8+z5+f7yqIvKQx+L0A1hHGxDrceytNR5Pum6eOms8pTgyvpOJpPCz8PhkEajkfaxOueJCUkUdLY7s2OBazTG1SJSjKDxZybtMPF8267PApdI3Y0fPza8Ic1DYRWTwCW6LHZWy0QCTNEXyjjYhlrxuYEIXGnTOARuN1gRuaPRiAbMKgeDAY3HY61jdc4TG5LdOt8FrdqB7HBQOT3t7NJaPHokb+I2OnceHDhI69CcUHN3Zpl8r8Y2pHm4WdOXdHymOTiooWFVb14eK8lQxsG2SH/X7GXTwxdFFXzTuKe3mQRWRO7Z2Zn4+5OTE61jdc7z5csXOj8/L3xCh6fLkYoiWWN9vWaQmz1Cy/KjylVp1COoSuvQSR6sckLK3fnpk6V8uqrNmh3vFvIpHZ9pjo/lVGOFRyH9wVxJekIIKd5Mwsf9r18pKIHric8BfMepjFCJVt1jpf/b3d2lW7duzT5ra2v6N+gZgpZ3Z0AeVGHil/N0TCsg3ePXr4azfqku0mFKJD4xhVqT3ajQld7SFhb8j+GIgJOTGisrgRSKCCnFmwn4e/w36tHcn+xpI0jjvKe3mgxWRG6/35/ztp6cnIhZEcqO1TnP06dP6cOHD7PPu3fvWv8dPiAZiPVVNQ/emj2ca2ojNZXxCnbSRbQzqZuD307I8WdlQrf2Mr9qJUTliXdI7BuYpnz6JL9shSR0vS0IZJnpY/lIy9QjosIT8fgNmo/zZcVLgBusiNytrS3x9xvCroCyY3XOc/36dVpZWSl8YoF7JKyuqnkgcCU8uAUtnMydnkzQfCIOLcxEomxDWq0m9mAlBFy+bFVWyJay83sidFOO7dzYIFqmzzOBmxHRBZG3b9DSuCetxgK3WJmOhsNh4efJZEIbGxszD+x4PJ5lTig7tuo8qSDFNloZgz0RuKFuYOKkInT5ROyBs9IIWaYW7KVN7IkdSaTixc2jqhw5a4tXr7zc0BnDy2IbDt5etX/2/bM4H7jgDSGG16XAD7ZOvLe3Rzs7O7S5uUkHBweF3La7u7u0ublJT548qTy27P9SIsvmx9zBwOCbojSgd1AUXYppCmkDE+f27XmnXq9neAB8/nx+J5Hxi8hwL64nDjBjTAV7mTew0MweC9yUkcZPopyZfPtWcYB7khZN359F3os7FbgdPhIlsY17MYGyvgEhVYUx8vQ6rGZWdSv37186WkKmdjWmNgwGchyLZfNOyTNYNpFlmeIAjzpwjLalS2Xab0+qaS0uFj3QvZ7aIx0dwjP4v3Sf/oOKndWXsSa2oipdEGTFs66JUeQSWXAUeTKoE1kU8R4gNfONG4a91FIDEllrxBQnYqkaIdHlDnBxg4xH8YMpvZCUEYLQTfZZSQbW69Hy/7rwdm4IrcqjjwRX8QzYw2i2BVVu1Y5GD18HMRNIf8vnz4aT8n/65DS2kAva2AUu0eXkxZ/lhNbnBS6R1wJXKpiQCpVp4qwkTK4PH8+TWgqX3iAvLtztS9FEytoIgesX8OQGiFS+W/spbm7KKac66g6phDI6+TsdeKJ4H0zBi8uZtsHFdy8uUXGDjE/9N1nPYAm+enSTfVY1Bkff5olkn5Vh4MkFBSQHkfZbrUcCV8KjWzGKKhmC0YJlDjxR/CUrNYFLdGmHWc5/y3eAa+XVtQjvW6FmKzGNjx5do0VjQqKmevVpXkj2WQUGPLkBI6XaqhXj6dnrML+dFLyClZuYbF6k5QV4LK6BU4YJa98LUqc4grfJX3zy6Cb5rFSZfQ4Oah/eRTsl+awsAU8uEOEeGWnP0Ry+jBAlxC5wieT881N6vatPK0+gJU8UBC6JbbhQ8lCnz9M1Scd31kTVf2fxlh3H6EaNyh2qELhEiY43oDHw5AaOlmaVDn7+nOjRI6P3pENjb3Qk6M6TjXbuGvREpZhRYY4aKfeceOprAG9TfSrHUsse3SSfVUOni7SlBHYVLvDkAiWScYmVcqTBZHW1U4G7vj7/u5QELtHl88uy+uXYf/+96OmdfkqrI6lG4AYllVLMqFBAlcqEBcqXTXquHICSfQE1lcUD4dE1S4tVxRJHLwAFIHIjgKcDmhsnVJV8Oi6szbPFPH/ezX34wB9/XAneJl4BVUWnwgHS7zR2vHFNnOTcrlFwo+xZumg7bl/wNlVT2UaWHij/evRjoYGwOd5GrsYjeHHDAiI3AqTl65khqixfcMFJHkL+MYV0rg6dyt6RF7w6g6i20P36tXaiZf715Ly4DSfmLJOd5jZ3Z2Pnd3OkjbAFHLy5RD0WSi/Wq6vap4m6jYAxIHIjQRp3/19vufTgJgLW1tsy3obL4aK3sZdQ+tLpqZzVvOScDSIdwqal5+nbt/nf8VKgJkGZ0ebUenkzKHSTK9TB8w8SNV5V5Hs9bXtzpdLYwG9Sm6qiJm9wf6EXtExX6RYyukxx1KOstVe27UCS5DK3JVRiV1vovn2r5f6TRFu0SIq+wXqy9nNqiA+5eUOHOxbF52TJoxt1xSzD2X1evWpxLwbo+vqgGojciMgb3HN6LFZhqksdb2GTsVxyGMLL1J6NjfnfaQtdhWuRn6fuJrlo4G1140bjtVKp2U2HFjx+XH1NUI7kWBQXO1oK3eQLdRjonHzss+VESc7jHglIIRYh/9Nbouv0lXpUX+BW9QJTKZGkgSJqz4VDpHLPRBXPp0EpzThHDAUWKl8sLMyfwmSbJv28DPLihcYLQ8P0Ykk9K25LBvNFumjHpJ6VY5BCDGjxb98FLhFRj4gWskwZ01l3Y1OWqVdo63p1pWMgcM0hlXsmauDRzbmspDzGSWGh8oUU82nK+4SJ2BySs16ZjKRBQy+zLRPRh3Hxjm8wX6Rtby4/X4N9cqAjIHJjo2IDURseParO/6mzIQkTsHkarZ7yGYJnWc+RVB5jHtxqME5Dek7Ia+sf/DlJKyXKg4lKDY9XqEwqW4nht2XXeXM7zr4JNIDIjQ0HZWCyTB2PNM3Xyj0e0XspPEK79kPNGvHJxaDxtWqVq7whvH15Xtu258NLpB20V0YEambtiwfeaBbelm1lWkh+T0LgQOSCRhwfl4/nX7+WpyfDBGwXVe0HpdBlX7gQHhpCS8xiM2wBmEN7rOJfEB4qrymC8bA9rjIdGH7XBZaByI0JHuTlYORsUqFLygQAzKMSurXip9nP0Vdg4vCwH0u2JJ12cVH/PPDi2oV779q8jCQXlsLTh1gcTExXQYMXN3yQXSEmPJjp6gwqcfY4f6m98Tu3nTwjosf0nF7SI/nY2HFoS8vL8/GZupfzwPSjR9qEWbrqnv9CLo1Mcs/K8R/ML7ex0TxmN7ln1RHIrgD06WgbfJVnF4OEe2pvRsttJ+8R0f+hx6XfTwbLMQSSUNK5JCZiN3AvIX8xKeV7wHVyXtwO4P2/ZB9tKfDixgFEbizwnV4db4NvmqoM2KGu0L0gmmVUTjY8lIf9ONj23rQamulCEkCNlFKs9BkJAfB8c2H0YyJX9R39wdykm4BY3DCByI2F0tw2ANQTurxoSLZhLyWdt2i56MwhhSpWCV1epC560dQxWhXrWO3r/+qlllKB2qcMaQh/TromDS9uPEDkxggyVQMFOmWae0TN1/piweHs9uiRLGpV6aaw9N0NvEsoqmHPMaBiSoXkXkh0kqhboE3UEby44QKRGwN8OzYyVYMSqibXC+pVFIGOGG5Ljmc3KTKCp5uaktzStydIXUIpoFJ2AfIYAebZtk1Te0i+ymNkQOTGQFKlcoAJyso0/0CsP6WUud4DW6oTn6ssLwucID0jXiCPiOYU8X/TkvL70dFR2E+eudWpBt7cpKo8RghEbmwkV5YKNKWsTHNhLlC5EmOnQ1uqEro8BD8J0eQ5vEDelPxGzuuU6N6Jjtyh0jur+DLyneSrPEYIRG7ocKtEWSqgiZgFI8VlVs9sSSV0URXND+pmxPgvWi2E/yTxQsLDfjp0h/ItKqqXEV7/hajzIQAYACIXADAPDzw0kYMHaFPHk5SEaPIUqe3zYqnXI/p3YnskOt6A5QQPwn6mSFtUpKqCfI9tclUeI+UHGyedTCa0v79Pw+GQJpMJbW9vU7/fb3TseDymhw8f0uHhoY1bjYv797u+AxArHsTXOcUTWzo+vpyQVZoBXl3/KEtI0iNK763EgzV/Xs6c25MU4y7lRQbhYaWs7927d2eidDKZ0M7ODu3t7WkfOxW/d+/eJd3bTKKsL0odAZvwerMx9y/Pbal2aWbQCVUvG/9Jm7RJb69i3WN+cJ7akvSMprfm6S0ng029ZtyTO5lMCj8Ph0MajUaNjr13757p2wMA1OXTJ1YpYtF5GiBwSZZdrnLnJ19MxP7APYWc/00HlOW3c/Z6eICOkZ7R5iZSgceO8eCg0WhEA5ZyaDAY0Hg8bnUsYPAyOxgwgW08irMzCt9u7aktXVygRLbPlCUQ2Nhwdx+dwm3Jk7CfKXwTmiRwYVtxYVzknp2dib8/EaK/dY6t4suXL3R+fl74OGVx8Wrrs7RN0zR1y+wA0IYu4unytuQi6FS13RoADcoSCBwcUHc7mfK2pKxBbAhuS69e2b2eJlXSIsWkMrHjbJunStC2PXbK7u4u3bp1a/ZZW1vTPkcr8l6ut2/dCF0AbMNz6NgWnUtL8x5j7K4CgSB5AWe/4zuZXPRrfg04R0o9tSjfGx+1Y3L//ve/09HRkfL/f/75Z9ra2qJ+vz/niT05ORGzK+gcW8XTp0/pr3/96+zn8/Nzt0J3YWFe6NqCF63H+gqIBV7lYMr6up2klfxlFLYEWuJNF1KJaFvxwAHZktQEnkVWAEMYz64wmUzol19+KaT8Wl1dpX/9619z4rXusb1eL4zsCnxQWViws1EHW0GBSx48KHqAbPW3Ks+WjevCloBr8n3u/n17S/pl9nTtmnm3JWwJNMSmXjMerjAcDgs/TyYT2tjYmInW8Xg8y6pQdWyeJiEMzuFG7WKjTgqJxUG38Em4q/AB29dFWARwja3wAcnhkke1YmLr+gB0hJViEHt7e7Szs0Obm5t0cHBQyJG7u7tLm5ub9OTJk8pjR6MRvXnzpvC94NKKmV4a4pWnkNIJxIBqFcTmZMkyu0SbPQIAyZZMzk2wJeApVopB+ECnxSBsLttgSQh0hc1lVlW/LsvgbuuaANgm3/du3zYbb8779cbG9/QOFf9n8pqwJaBBUOEKgOYTJtryRpUlZgTAJiaXWaWJd4o0WfKNlyZA7iDQFb//bu5cUoqwvIitkyi2LbAl4BEQuTaQEiaamJgXF6uvA0BISHbBPUs8k74JUcCL1SN3EHCJrb0U/OWTvyRKiWLbOmFgS8BjIHJtwRN/m5iYEecEuoRPmCZyQXO7kDy30nJq24nZ9sYbAMrgeyl4TGsT6p7DdCgBbAl4DESuLXjibyKzYQtdVKICIE/bpU6+MlGGzRg/voQLgGtOT82fQ8dmTM1NsCXgGRC5NrG54cxGYnwAqjD5osZXJnTtpem98KXiBmXEAWiNydhVnRdGItnWeNhBk+vCloBnQOS6BvkDQchwYdp0mbWJHZjahIad38AHeOxqE5E5pckLIy/x1STsACF0wHMgcm0jDTYvXrQ7J+oPAl8wscxKVF94mt6Exs8HQFc0jW1t6jiRUgC2ccIghA54CESuC3ic0uPHet/nA4+tMpAA1KHtMivvzzrna7sJjR9rIkcoAE2xEcOqs1LRxgmDEDoQABC5LpDilHjlMgBCoc0yqzSB6qYcQrgBiAU+N+jG1rZ5YVR9R9cJA4DHQOS6gk/Mnz/X+x5P7o0JHviGzjIrn0BNhQvU8eZygc3T/AHQNW1jXJvkqJW+U2VP3JYQQgc8BWV9XdKkRCnKJQIfGQyK8bh1+iX/Tt3vqeC2UVUeFbYEfGR9vRhbXrdf8v68utouu4GOfcCWgEFQ1jcWMBCAWOCTaZ0KTiYFLpGdSmgAuIa/mDXd/GU6fRcyAYEIgMh1DRcDZQMJT48EkQx8papv2ihjqrMJjVdngy2BkOH93ETYj2QTUopA2BIICIhc1/ByjkTq8qjwTAGf0ZlY+URoamKse5621dkAsIlOTKs0X5jKEsJtWkoRCFsCAYGY3C7Y3JwfKKTHkH9bX1iQBTIAXVInNq9JLHqbe5DOnz+m10MSe+AfdeNcbcfDVtkTbAkYBjG5sSG9dfPUMXyZCAIXxILpSVk6X97bxW0JkzIIFRepJ8vsE7YEAgMityv4QMIHC1OVpACwSVUaLu4VunHDzn3wakv5lRLYEggBbktSTmmeetLWQqxq7whsCQQGwhW6pGxZKP9/1641y38IgAtUy6cPHhC9fi3/n4v7mC6lwpZAKJSFIiws2Ittr3MvvV7xerAlYAiEK8SKaoDiFaQwkIAQ4QLXdsJ4SQBwjxRsCYSKS4Fb53qwJRAAELldw5dvez29ClIAdA2f/F68kGMHX71ycz954lyoArHC++s0tlzyqgIAKvmh6xtInk+fyges1VV39wKACXjpXiJ3YjPL1PYEWwKhoUrX5WrDl8qeTJXjBsAy8OT6QNnmHdNVbACIHb4JbQpsCYRAldf22jV390IkC1pTeXkBsAw2nvmCyvsU5+MBMeJTH+b3sroKkQvCoWx1D3MCiAybeg3hCr4gLQthMAOh01XsIGwHxIhqlQIAIIJwBZ/IsuIHgJCQ+iySxQOgjyqf9PGx2/sAIHDgyQUAmAMvZwC0R9qQDNsCQBuIXAAAAMA3IGoBaA3CFQAAAAAAQHRA5AIAAAAAgOiAyAUAAAAAANERbUzuNP3v+fl5x3cCAAAAAAAkpjrNRtmGaEXux48fiYhobW2t4zsBAAAAAABlfPz4kW7dumX0nNFWPLu4uKD379/TzZs3qecgIf35+Tmtra3Ru3fvwqiw1hFop2rQRtWgjeqBdqoGbVQPtFM1aKNqpDbKsow+fvxIP/74Iy0smI2ijdaTu7CwQD/99JPz666srKBz1wDtVA3aqBq0UT3QTtWgjeqBdqoGbVQNbyPTHtwp2HgGAAAAAACiAyIXAAAAAABEB0SuIa5fv05/+9vf6Pr1613fitegnapBG1WDNqoH2qkatFE90E7VoI2qcd1G0W48AwAAAAAA6QJPLgAAAAAAiA6IXAAAAAAAEB3RphBzyWQyof39fRoOhzSZTGh7e5v6/X7Xt+WE8XhMo9GIiIgODg7ot99+m/3t4/GYiIju3LlDk8mEzs7O6M6dO0RU3maxtaeNdoitjYiI9vf3aWtri4ho7m9JuS+Nx2N6+PAhHR4eFn5vo++E2l6qNsL4dEVZGxFhfJqiaieMT1eU2ZV341IGWnPnzp3Zv4+OjrJ79+51eDdu+fXXXwv/zrfF9vZ2RkQZEWVbW1vZ6enp7P/K2iy29rTRDrG1UZZlszbKf6b9K9W+tLe3lx0eHmbSUG2j74TYXmVthPHpkrI2wvh0RVk7YXy6osyufBuXIHJbcnR0VHgAWZZl/X6/o7txy+HhYeFvPTo6yogoOzo6yrIsy16+fJmdnp4WDH56nKrNYmxP0+0QYxudnp5me3t7hd/lB9LU+xKfdG30ndDbi7cRxqd5JPGG8Wke3k4Yn64osysfxyXE5LZkNBrRYDAo/G4wGMyWL2Lmzp079Ntvv81+Pjs7IyIqtEe/359bVihrs1jb02Q7xNpG9+7dm/17f3+/8DMR+lIeG30ntvbC+FQfjE/VYHy6pMyufByXEJPbkukD5pycnLi9kY7IG/o//vEP2tramhn62dkZ7e/vE9Fl3M5f/vIXGg6HpW0WY3uabocY2yg/OZydndHJyQkNh8PC79CXrrDRd2JsL4xP1WB8qgbjUxGVXfk4LkHkWkL1YGJlauT5YP18cPhwOKSff/6Zjo6OSs/R5P98x1U7hNxGeXZ2dujXX38t/A59qR42+k4M7YXxSQ3GJz0wPl0h2ZXqOJf/lwfhCi3p9/tzbxQnJydB7JA0yc7ODr1586bwd08mk9m/p7siJ5NJaZvF2J6m2yHGNppydnZGo9Fo7m9BXypio+/E3F4Yn9RgfKoPxqci3K58HJcgclsyTSnC2djYcHwn3fHs2TPa2dmZLc+cnZ3ReDymP//5z3PHDgaD0jaLrT1ttENsbZTn7du3c4MX+tI8NvpOrO2F8UkNxic9MD5dIdmVj+MSwhVako/LIbp8o9vY2AjqbawN+/v7dOfOnVlH/+c//0nb29s0HA4LSzqj0Yju3bs3ezPLk2+zsv8LERvtEFsb5RmPx3ObDNCXLjk7O5vde9m4k3K/yrcREcYnCd6PMD7J8L5EhPFpisqubPSZtu3Uy7Is0/vzAGcymdDLly9pc3OTDg4O6OnTp0F01LZMJhP605/+VPhdv9+n09NTIrpKGN3v9+no6KgwEJS1WWztaaMdYmujKc+ePaOjoyN6+fJl4fep9qXRaERv3ryhZ8+e0ZMnT2hzc3O26cNG3wmxvVRthPHpirJ+hPHpirJ2IsL4RFQ97/s2LkHkAgAAAACA6EBMLgAAAAAAiA6IXAAAAAAAEB0QuQAAAAAAIDogcgEAAAAAQHRA5AIAAAAAgOiAyAUAAAAAANEBkQsAAAAAAKIDIhcAAAAAAEQHRC4AAAAAAIgOiFwAAAAAABAdELkAAAAAACA6IHIBAAAAAEB0/H/t3Jj4Bb9DogAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 800x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 14 s, sys: 159 ms, total: 14.2 s\n",
      "Wall time: 14 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# Output can be collected at each integration step\n",
    "# Note, container is overwritten at each call\n",
    "\n",
    "FODO.ns = 0.01\n",
    "FODO.output = True\n",
    "\n",
    "state = torch.tensor([+0.01, 0.0, -0.01, 0.0], dtype=torch.float64)\n",
    "orbit = []\n",
    "for _ in range(16):\n",
    "    state = FODO(state)\n",
    "    orbit.append(FODO.container_output)\n",
    "qx, _, qy, _  = torch.vstack(orbit).T\n",
    "\n",
    "plt.figure(figsize=(8, 2))\n",
    "plt.scatter(range(len(qx)), qx.cpu().numpy(), s=1, color='blue')\n",
    "plt.scatter(range(len(qy)), qy.cpu().numpy(), s=1, color='red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "c1938592-11da-488a-8213-be052456ac2b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 17.8 s, sys: 5.1 s, total: 22.9 s\n",
      "Wall time: 24.8 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "tensor([0., 0., 0., 0.], dtype=torch.float64)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# Functions can be compiled, but note that dynamo unrolls loops completely (torch 2.4)\n",
    "# This leads to very long compilation times\n",
    "\n",
    "FODO.ns = 1\n",
    "\n",
    "state = torch.tensor([0.0, 0.0, 0.0, 0.0], dtype=torch.float64)\n",
    "table = 1.0E-3*torch.randn((512, 4), dtype=FODO.dtype, device=FODO.device)\n",
    "\n",
    "fodo = torch.compile(FODO)\n",
    "fodo(state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "66909d17-533d-4f58-a02b-ddd5f946277e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.71 ms ± 57.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "\n",
    "FODO(state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "74810df4-e8b9-4f9e-a518-7c4482c976a1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "458 µs ± 2.08 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "\n",
    "fodo(state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "8596d5b3-8591-41d7-947e-5c920ba4198d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.5 ms ± 46.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "\n",
    "_ = torch.vmap(FODO)(table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "d21f04d3-14c4-411f-b195-f407bc0646c1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.89 ms ± 514 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "\n",
    "_ = torch.vmap(fodo)(table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "34cfa8a9-6808-4b70-9e34-63ee7b55e0b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note, compositional operations (map or jacobian seems to break/ignore compiled version)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "b05bf635-b14b-4fb1-b81b-362afb47a92c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'FODO'"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# name (property)\n",
    "# This property can be used to get/set line name\n",
    "\n",
    "FODO.name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "fe8e9207-318c-4663-a455-313ed530f488",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# sequence (property)\n",
    "# Contains ordered sequence of elements\n",
    "# Elements can be added or removed from it\n",
    "\n",
    "len(FODO.sequence)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "300488b2-180f-4519-9366-4d528d1c86fd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'BM': ('Dipole',\n",
      "        tensor(3.5000, dtype=torch.float64),\n",
      "        tensor(0.3927, dtype=torch.float64)),\n",
      " 'DR': ('Drift',\n",
      "        tensor(0.2500, dtype=torch.float64),\n",
      "        tensor(0., dtype=torch.float64)),\n",
      " 'QD': ('Quadrupole',\n",
      "        tensor(0.5000, dtype=torch.float64),\n",
      "        tensor(0., dtype=torch.float64)),\n",
      " 'QF': ('Quadrupole',\n",
      "        tensor(0.5000, dtype=torch.float64),\n",
      "        tensor(0., dtype=torch.float64)),\n",
      " 'SD': ('Sextupole',\n",
      "        tensor(0.2500, dtype=torch.float64),\n",
      "        tensor(0., dtype=torch.float64)),\n",
      " 'SF': ('Sextupole',\n",
      "        tensor(0.2500, dtype=torch.float64),\n",
      "        tensor(0., dtype=torch.float64))}\n"
     ]
    }
   ],
   "source": [
    "# unique (property)\n",
    "# name: (type, length, angle) data for all unique elements\n",
    "\n",
    "pprint(FODO.unique)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "c5ae2b68-a332-4da5-b2cb-c039a2dd2764",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor(12., dtype=torch.float64)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# length (property)\n",
    "\n",
    "FODO.length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "691b3bf5-dbe1-4e90-ba7f-d9c930cf9092",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor(0.7854, dtype=torch.float64)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# angle (property)\n",
    "\n",
    "FODO.angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "157fbe77-ba7e-48f0-8b27-a2d7bbf39b7b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'QF': 10, 'DR': 10, 'SF': 10, 'BM': 10, 'SD': 10, 'QD': 10}\n",
      "{'QF': 5, 'DR': 3, 'SF': 3, 'BM': 35, 'SD': 3, 'QD': 5}\n",
      "{'QF': 5, 'DR': 1, 'SF': 25, 'BM': 35, 'SD': 25, 'QD': 5}\n"
     ]
    }
   ],
   "source": [
    "# ns (property)\n",
    "# This property can be used to get/set number of integration steps to unique elements\n",
    "\n",
    "# Set value integration steps to all elements\n",
    "\n",
    "FODO.ns = 10\n",
    "print(FODO.ns)\n",
    "\n",
    "# Set ceil(element.length/value) integration steps to each element\n",
    "\n",
    "FODO.ns = 0.1\n",
    "print(FODO.ns)\n",
    "\n",
    "# Set by name or type\n",
    "\n",
    "FODO.ns = (('DR', 1), ('Sextupole', 0.01))\n",
    "print(FODO.ns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "ef28fdcc-107d-49b1-9be0-c15428430d0a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'QF': 0, 'DR': 0, 'SF': 0, 'BM': 0, 'SD': 0, 'QD': 0}\n",
      "{'QF': 0, 'DR': 0, 'SF': 1, 'BM': 1, 'SD': 1, 'QD': 0}\n"
     ]
    }
   ],
   "source": [
    "# order (property)\n",
    "# This property can be used to get/set integration order to unique elements\n",
    "\n",
    "# Set value integration steps to all elements\n",
    "\n",
    "FODO.order = 0\n",
    "print(FODO.order)\n",
    "\n",
    "# Set by name or type\n",
    "\n",
    "FODO.order = (('BM', 1), ('Sextupole', 1))\n",
    "print(FODO.order)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "5fefd15f-c3b5-4df4-84c3-1ac926cfd7ce",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'QF': 1, 'DR': 1, 'SF': 1, 'BM': 1, 'SD': 1, 'QD': 1}\n",
      "{'QF': 0, 'DR': 0, 'SF': 0, 'BM': 0, 'SD': 0, 'QD': 0}\n",
      "8\n"
     ]
    }
   ],
   "source": [
    "# Nested lines\n",
    "\n",
    "FODO = Line('FODO', [QF, DR, SF, DR, BM, DR, SD, DR, QD, QD, DR, SD, DR, BM, DR, SF, DR, QF], propagate=True, dp=0.0, exact=False, output=False, matrix = False)\n",
    "RING = Line('RING', 8*[FODO], propagate=True, dp=0.0, exact=False, output=False, matrix = False)\n",
    "\n",
    "RING.ns = 1\n",
    "print(RING.ns)\n",
    "\n",
    "RING.order = 0\n",
    "print(RING.order)\n",
    "\n",
    "# Note, sublines are not flattened\n",
    "\n",
    "print(len(RING.sequence))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "55555b52-dfc1-4434-8c1b-a7fa5091cea0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([0., 0., 0., 0.], dtype=torch.float64)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Tracking, mapping and differentiation is similar to flat lines\n",
    "\n",
    "state = torch.tensor([0.0, 0.0, 0.0, 0.0], dtype=torch.float64)\n",
    "RING(state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "c402d469-4119-4b52-ab7c-ea57c02ad00c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'FODO': {'BM': {'dl': tensor(0., dtype=torch.float64),\n",
      "                 'dp': tensor(0., dtype=torch.float64),\n",
      "                 'dw': tensor(0., dtype=torch.float64),\n",
      "                 'e1': tensor(0., dtype=torch.float64),\n",
      "                 'e2': tensor(0., dtype=torch.float64),\n",
      "                 'kn': tensor(0., dtype=torch.float64),\n",
      "                 'ks': tensor(0., dtype=torch.float64),\n",
      "                 'mo': tensor(0., dtype=torch.float64),\n",
      "                 'ms': tensor(0., dtype=torch.float64)},\n",
      "          'DR': {'dl': tensor(0., dtype=torch.float64),\n",
      "                 'dp': tensor(0., dtype=torch.float64)},\n",
      "          'QD': {'dl': tensor(0., dtype=torch.float64),\n",
      "                 'dp': tensor(0., dtype=torch.float64),\n",
      "                 'kn': tensor(0., dtype=torch.float64),\n",
      "                 'ks': tensor(0., dtype=torch.float64)},\n",
      "          'QF': {'dl': tensor(0., dtype=torch.float64),\n",
      "                 'dp': tensor(0., dtype=torch.float64),\n",
      "                 'kn': tensor(0., dtype=torch.float64),\n",
      "                 'ks': tensor(0., dtype=torch.float64)},\n",
      "          'SD': {'dl': tensor(0., dtype=torch.float64),\n",
      "                 'dp': tensor(0., dtype=torch.float64),\n",
      "                 'ms': tensor(0., dtype=torch.float64)},\n",
      "          'SF': {'dl': tensor(0., dtype=torch.float64),\n",
      "                 'dp': tensor(0., dtype=torch.float64),\n",
      "                 'ms': tensor(0., dtype=torch.float64)}}}\n"
     ]
    }
   ],
   "source": [
    "# Deviation table contains unique lines with unique elements (unique by names)\n",
    "\n",
    "pprint(RING.data(alignment=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "202f4feb-5c02-4837-ac84-0b1f94b560e3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([-0.0047,  0.0005, -0.0045, -0.0014], dtype=torch.float64)\n",
      "tensor([-0.0055,  0.0005, -0.0048, -0.0014], dtype=torch.float64)\n"
     ]
    }
   ],
   "source": [
    "# If this can be used to pass different values of deviation variables to unique elements\n",
    "# But values should match for consistent differentiation\n",
    "\n",
    "LA = Line('LA', [QF, DR, SF, DR, BM, DR, SD, DR, QD, QD, DR, SD, DR, BM, DR, SF, DR, QF], propagate=True, dp=0.0, exact=False, output=False, matrix = False)\n",
    "LB = Line('LB', [QF, DR, SF, DR, BM, DR, SD, DR, QD, QD, DR, SD, DR, BM, DR, SF, DR, QF], propagate=True, dp=0.0, exact=False, output=False, matrix = False)\n",
    "\n",
    "RING = Line('RING', [LA, LB], propagate=True, dp=0.0, exact=False, output=False, matrix = False)\n",
    "\n",
    "state = torch.tensor([0.01, 0.0, 0.01, 0.0], dtype=torch.float64)\n",
    "\n",
    "kn = torch.tensor(0.01, dtype=torch.float64)\n",
    "\n",
    "data = RING.data()\n",
    "\n",
    "data['LA']['QF']['kn'] = kn\n",
    "data['LB']['QF']['kn'] = kn\n",
    "\n",
    "print(RING(state, data=data))\n",
    "\n",
    "kna = torch.tensor(+0.01, dtype=torch.float64)\n",
    "knb = torch.tensor(-0.01, dtype=torch.float64)\n",
    "\n",
    "data['LA']['QF']['kn'] = kna\n",
    "data['LB']['QF']['kn'] = knb\n",
    "\n",
    "print(RING(state, data=data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "c6537d61-ed50-4d9e-9029-a1bd08b36b73",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAFeCAYAAACB7binAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABmy0lEQVR4nO19TW4jSZJ1UPqoAgRUtcRlo0sa6Aga9g1UN9AmDzDqXs+Cg14lelXNE7BVV2DdQHmDHuYNKGCQQC1TmVWb6k37t6g25qPJzN38J8igZA8gUskIenhEuNuzPzcfhRBC53A4HA7Hv3G07w44HA6HY1hwYnA4HA7HFpwYHA6Hw7EFJwaHw+FwbMGJweFwOBxbcGJwOBwOxxacGBwOh8OxBScGh8PhcGzh/+27A0PDv/71r+6nn37qvv766240Gu27Ow6Hw1GNEEL3yy+/dL///e+7o6O0PeDEwPDTTz9133777b674XA4HM3x4cOH7g9/+EPyPCcGhq+//rrrut8e4DfffLPn3jgcDkc9fv755+7bb7/dyLcUnBgYyH30zTffODE4HI4XBat73IPPDofD4diCE4PD4XA4tuDE4HA4HI4t9BZjeHx87H788cfu6uqqe3x87O7u7rqzs7Psc0uPdV3XvXv3rnt8fOyurq66ruu6m5ubvm7X4XA4Xg5CT7i+vt78vV6vw+3tbdG5pcceHh7C3d3d5tjV1ZWp358/fw5d14XPnz+bznc4HI6hI1eu9WIxPD4+bv3/6uqqe/fuXfa5pce6ruv+9Kc/davVanPs4eGh4E4cDofj9aGXGMO7d++6yWSy9d1kMunev3+fdW7pscfHx+7jx4/d2dlZ9/79++7Tp08bd5LD4XA44uiFGD59+iR+//Hjx6xzS4+9f/++m0wmm/jDDz/80P3444/i+f/85z+7n3/+eevjcDgcrxk7XeCmCfPcc1PHPn782D0+PnY3Nzfd2dlZd3d3152fn3chhGfnf//9991f//pXc78cDofjpaMXi+Hs7OyZdUCunZxzS49dXV1tzqFrdF0nurL+8pe/dJ8/f958Pnz4kHm3DofD8bLQCzFoaaHT6TTr3NJjOfGEr776alP+wstgOErw97//vfuP//iP7u9///u+u+JwNEEvxMAF8+PjYzedTrc0d8oqip1bc2w6nW5cTrSW4fr6uuFdOhy/4W9/+1v3f//3f93f/va3zXdOFo6DRl95s+v1Osxms7BcLsNsNgtPT0+bY7e3t2E+n5vOLT329PQU7u7uwv39fbi7uwvr9drUb1/H4MjFYrEIl5eXYbFYbL67vLwMXdeF4+Pjre+tvy+5ZmlbjpePXLnWGzEcKpwYHIQaIbtYLMLx8XHoui5cXl5G2z89PQ1d14WjoyPztYh4eNs5hOR4PXBiqIQTg4OAwrcPrX4ymYSu67Y+GonktJ0iJMfrgxNDJZwYHAQUvpqGXtImCe7JZBJGo9GGFI6OjsKbN2+a9jv2neP1wImhEk4Mhw2rAMTz6O83b96EyWQSJpOJquW/efMmS8Dy/nBXD1kNR0dHvbqAJBeTk8XrgRNDJZwYDhskACeTSVTooQWAQpO09+Pj4zCdTsPx8fGWFh+zHGJBaDqfn4OEgy6gmNDOCTzH2neX0+uBE0MlnBgOG9x/nwr8ShYDae/o4uHnSgJbIo0crTzlukrdm/QbHichiwjjG6enp71kTjmGAyeGSjgxvAzkCjKJKC4uLkLXdWE8Hps06xbCM+aykqwhqd+x2AK1QZlQ1qB3Lek59gsnhko4MbwepDR0+u7o6Cgad2gpGGtdVSlwq2M0GpmC3tK1qY3JZGK6tmN/cGKohBPD6wF3s0jadswPnxLKuXECdPWUBretoGtRVlRJnMGJ4XDgxFAJJ4ZhQgva1vjGLW3EzqGArqZtE3HwTCCMY5BATpFQK+tEspJKM6HclXQ4cGKohBPDMMH966itxoRTrqslBLvAs1gMXNjTb7hARiEdiy9Y4hwxqyNlJeXcv+Nw4MRQCSeGYYL7xjlBlLhztPNimj4K3FTb2m8kwS2llZasOUDiuby8DNPpNHRdF6bTqbkdTsIlbir6nZPMMODEUAknht1DchOR8MRFZ9PpNBwdHW3SK0t88tL1Qtj2l6OmzwlI86mnAsMlgl0jvJgLiz8TbEs6Hns+KavMEpfxOMQw4MRQCSeG3YO7Sbi7hdcTOj4+Fn+Xez3UyLkAk4QjxQUkISf1JZX1JCGWrkpAi0L6Lf6GWwxIcNa+lKyb4Cu7R6NRNoE72sGJoRJODLtHjsWAmnKpy0Ly/Vv87TWB7JYuFc1isJAPCurazKdYZtXp6enmgwsGverrfuDEUAknhuEiJlxzrAfuTnkpfvCcLCtLzafS58LdYaenp0lXnKNfODFUwomhHWoFLv99TPjnCDzNdZXrkjpk0D3HAszac0lZEGTl0RoJIgKPN+wPTgyVcGJoB4vw0ZBy92jCSYsVSBlAuRk/lj4fiuWRih/QOVKgWnuvnEhKn/MhPcdDgRNDJZwY2sEifDTwIKa1TU4Mu7QGUoKxBVq3mdL+peesZS3lliTXgO++VZuvHU4MlXBiaI+SwKz2G9RWMUCtWRK71D4115cWcCWN/OTkJIxGI9MmPRLRlabtWq4Ts/ZqiD8GtBa9NHgbODFUwolh97Bo9TyGwBeC7UJ45JKM5A5D0H3zNNzcPmA7LYWz9V6RmKQ9LGqu7xZDGzgxVMKJoQw1mnnqt2/evHkm+PbhbrBo0RypVFyevVP6/FpbDNbr0r2VkJzUjqMfODFUwomhDCm3iQWagMBFbqUuopxUTnJJUT4+T2u11mmS2uXaPa4raJmxs4uUXLT0qP2YxaC5DN1d1D+cGCrhxJAPEkJSCeccgaS5lFJVTC2waPt4fa4B474FKd96LM1WE9i0cQ59L/XRGotBYdtnaYrYu5WsGOn9xhSKnNiTIw4nhko4MeSDJre0oY0m7CX06RKxpmdKFgOu3LVs88nv2Wqt8G03+e+lvseELdWVmkwmm13oUiueWwlebhnxti0xBOneLOnIjudwYqiEE0M+0FfONfLcIGJMg2wBqyDhQowEq7TNZ6tMKM2q4d/HsrD49bnlkyLpElKTkHJlaQqDRB54bSKG09PTJNE7vsCJoRJODHHEJitZDNJEtVoO+/Q5SwKV+lBqFeUIVqsLyfos6be8ECEV1LP0oZULSiKc2OI57d649XR6etokC+qlw4mhEk4Mz4GTmHzh6JLgJastwk37LvZ9rH+5Wq30GxRKkhWgEVbs+poGLgVpc2IIJc+HYkA5GUOtiMEac7DEk2JuN4cMJ4ZKODE8B7oj0N/OhV1O2mhLy6DE/SEJpZQQLiEyLhDxWaKQxuch7ewWe0bWfklCN/WsWvrwLe/JYjHgc8J1LftI2T0UODFUwonhOVC4xSag1Q2AArJFLIETkyWd1CoQpbY0YS/dNydAzWKQCCNGRgh8lm/evAmnp6dhNBptrLsYqZS8BwtZpn6Xc04sZkL3ht+79fAcTgyVcGKwQZvAGmloArJUGHGiImFIfmctfdZ6P1qQ8/j4+FkJ6ZigJIElxSa056MJOM3NwrVoTLHlazFi15WINOVS20XFVE58qATgeJDWnLj18BucGCrhxGCDljaIAlRyEZRaCKiVo3Yo7fBGAg6FY+qavH/S/aEA1XZz09wlFl94yiUn7dzGf4MWA5FmSpOO9ZHfD39OZJWMx+PeXDm0C512D5LFFrPiXiOcGCrhxGCDVsV0MplUbWavwWIxxFYppwQEL7vROt7QQoOVLIaUAMSV1Sm3mURIksVAz5tSd3nsqbX1gG1b3YJuMWzDiaESr5UYalw7/P+LxSK6P3Jffa4REJK74hCEijVWIrlgCBK5xNpFi40+b968af7OqQ/WxXkOHU4MlXitxGAxvWPCQvPRW4VELD6hnc+zmtBq4f1Ebdsa4HxJbggkPokELKmkeC5ZDFQmhI5ZSMoq3CUr9CW9k13CiaESr40YUn5tRExYxHZbswC1UMvkR0FH1yDhL2XjoHsrJfhfohtCe8+ayw3Pt2ZgpcAD+qnna3HdOWxwYqjEayOGnImuTU5JSOciZTGk3EY8OMyzVFIWw2tFzMrC4zRGagoaYryjhFysVqXl/b62MeDEUInXRgy1sYWSNkqgERi3Ekho8eqiu8QhCZ1UX9+8ebNZE8E3R8q9vxzrVELKqoxZjFpbfdXkGhoGQwzr9TrM5/OwXC7DfD4PT09PReeWHkPMZrPo9REvlRhaCSsU0KVt5sYT6DfStUhQjUajrT5Z1g70hRbPiGNfZIP3goJZcynGrLpapMYNjQWsgBtrq4bkDg2DIYbr6+vN3+v1Otze3hadW3qMsFqtQtd1r54YuMZd6hJoEaSVNL9S87+kzEMLWAPxMc10sVhsrVLOie+UkGsJJD8/pgjHAtbWBIQaBQN/Z62zhPfzWoLagyCG9Xq9JbRDCOHs7Cz73NJjiOVyGa6url49MfBJlKNdxdos0fwlYWghGTonVupiV7D2JRaYx8VnKQHK31/KrdI3YuROriLcfAjB76V0j4USxUQiWE58L9F6GAQx3N/fh5ubm63vrq6uwmq1yjq39BhhuVxuvn/txMCR8sf2NTl5OQlJG01p4prA6QOa4NCyebQ2kEAxO8dKDKk2pf/3KeRi15NIM2ZpauVHUmMqJ95lGV8lRHMoGAQxzOdzUWg/PDxknVt6LIQQnp6eNn/HiOHXX38Nnz9/3nw+fPjwYomBC7ZYHR3rJEkJIH4cy3ZbSjBIaFmfJ2X1oJDjdYgww0lLA5UsIuw/ZurQSu6UsJP+j32z5P23Jg68P+mZ8uPS/SBBlGY+aePHMq5abCE7VAyaGEiDt55beiyE3ywR/F4jhrdv325pbfR5icQgCSdNwNaY2CgYaNXq6elpCGFbO5QEskResXNqQc9EExqa2wefnaQF8/+nLA8kTC5AkUgkjZquRVt5YgkJbVOe1toxF6q8/ZjQ5UqK1WKQkLIYcuI4LwmDIIb7+3vR/y9ZDLFzS489PDxsEYFbDL8BJ0eKGBC5EwaFLX2Ojo62+mCZ/Cmh3QLWOAkJtouLi43wxr7VunKkxV9orUgEIVkQ/Nlrm/Lgb1poyjEilI5Lv6X+asHtFKykoJ1Hz2E6nTZTPIaCQRCDFhiWhHPs3NJjRBz06bouzGYzMcbB8ZJiDDHhZBWIqXa086ltaaey2PXx+yFsvqIRGf3/6OgoKfyk58eFsSW7yvIeuFvp4uIieY9S1dZc5LoUQ9gWxJPJdsVavg6lRuNPWXP8PLr2S9o2dBDEEMLzVFJ0+6xWq7Ber03nlh5DdF23db0YXhIx4EDPEaxU5ji2N3AONG2ahKu0Y1lLK8HqFpOOca09JqQ1zZu/h8VisaUhozDHe87tK8Y7uMUQa2sfvnV8Bpo1ietQLFZrK4vh4uIiHB8fF22FOlQMhhjW63WYzWZhuVw+W2B2e3sb5vO56dzSYyH8FoCez+eh67pwd3f3aiwG9NmSFpYzodAV0QJc2KOw5RU5cyyZ3OvjJjsSWUrCB4Ut/atlA2nuMV7SG58HtUv/p30NUsIQXU/Sjm2az9+S6lsTX7ICXZkUf5KuH/uuBlJ7vGSHWwyODQ6ZGKRAqUUQcAFktRisk5ULexRaOBlrJz0KcXRlYTCWKoJKsZWYZolaLGrjXMBJApWeL2nAUgaTtHsb/Za74xaLxUZ4YbwjtmNbzErj9y1lTaHLrMU7wlTdXQleicTx+lKc56XAiaESh0wMqBVyrTsm9Eq19JR5r2l+fIvPVqtP0TLgwW90W1gFALZHgoxvGypZPfzZaHEK/qw0oc5jANiv3B3bpOfN36OWgdWidASO0V2vOkbLDO8R3W2xedCX9bQLODFU4pCJITZwtclPAiVnknL3SkqwYpso1HKCqhZwiwH3DDg/P99YQRZBHcK2C4hvEZpyt7R0h3C3UKwdFHg8u0Zyi0nvEf3//FO7dqSFQhJrMwbaTIjcdVJihETo1G5KERoynBgqccjEEIPmLkD3iFUjtE4QaeLXkoBVmHDNWOrzYrHYxGGkctyoYR6KQMD7TlkG2jO5/Lf7CsdHbLV5zTuVFAWCNX0UXV/Wa8UsKiQFPM8thleMl0QMMSG6WHzZfjPlUom5oVJA10GLSYXCOiaoucDRhExsAx+rZTQ0aP3mz0AKqKMVlbKKCDW++ZgrUYq7SLDGqXLf5yETAYcTQyVeEjGkhKh14KPAzE1ttLptrEiRnSTIUPhwMkR3wksSBBI0jRnHiRSfSoFn8+S+Y+29pSwGFPSvpUpqKZwYKvGSiKEPHy5qcTn+89K+lFgnKBx4Tj9uLP/aBEnMksDNeDCwnmoLhXkLy0p6h9q1S6wb7R4s3x8ynBgq8ZKIoQ+goNWycDRYztF+g4vDtEkrkQ8nspecklgKfC9oPWjptyGk/fqlwtXitqQ+au5J6zjT7qFknA4dTgyVGDox5E64FtoPN+8lYrBobCV9x9RQFFzapOXHczJ6Xiv4u6OMLum5kzDmwehc4rBCet+pd2h1d2p9lCygQx83TgyVGDox5GozLbQfbAMXVGmTpJXGhZlTqTLXhMVisZXXv08csjBBS00KThNB0HuI1ZKqcTGlLETpWG4w2jKOU2N66O/aiaESQyeGUoshtV9A7Br4f1zo1aqPCFx1zTeywUmJE5Vfr0VRuBbgfcyNsexT2OC18T6spTakgH+u9SDdv9YvRO7Kf8v1c0hkiHBiqMTQiaEUfODGBrI22UPQJyuuBygBZpbgpOYWA++HlJGy66JwFtcZPVOeDktCU6qCyt/RvjaSQVLjloHWp8VisQn0Y7XSHCEtjdGYQkCgAPpoNNr0IUVKte5PtxheOIZKDLUDzxLUo+9JyGvbfvLfo9Ar1dI5KaRWKPNjLco1lCJGpAQUrpQiyy2ilFBCS6jEAunjHrnw5u+FSJ2XESkVtpZxzK9tIQZNUcqxBIa8A5wTQyWGSgw0QC3ZOaVA81/yI+N1eapgK4tBWk9ApQxOT083faB+tlo4x6EJXrRUJpMvO9R1XbelIUuCi/4mKwEL8kkaNf8txU6IWPB3sd3RWjwfq0Z9cnKyeQY4JtDv32JNC9+zvOueF8MjBcOqENVYDCEMx4UpwYmhEkMkBhJSuDhL8rvXDm4iH14XCLVFFMh9aawY4ORklSqd3QpoBUmuDC7U6cPLnKPGSX+jIF8sFuqOcChoUJhiH5CYJKvDYs20BPYFwd9pbV+kQolYDE9THKwWQE4f6VxtU6ohwImhEkMkBj65tcyglFlfOhkk9wBqgS2Civg9pkKipol18vsWcprFQJoqF8hc65cK1WlEyhfhcfccF3D4/KS4TOqZ9vW80BriY7O1hYu1nCSLSRv7uUqSxaqRzh1azMGJoRJDJAY+yHDj+Njkq9EWNf8uZpjkEkOKqHh/yYVUej0r0D2UsrikZ8qfCwpsbjlgqWzuZqE9D5CUYrn05F6KEQOhr2dHwOeyi/IU+Dxj4zpXQKfGgtSudI0cYtkFnBgqMRRiiA08SfCk2siFJkhqXAIpotKIrc+4Cl6Ha+38OJIDtwRoJzVcGIbtcKsA/0bhLmX+ENCF1nXdxp9P8Rd6J5K1w8eTRDw1wPZ5TCh2ruX71PVawirQtfOsFWF3DSeGSgyFGFIaWN+mampVaMl1U7/VhFere+TCi/b3lVxC0rNGdxAJcy3WQNYct5JIYNI+EXxzHSROXuCPxyLw+viekOikQDi/nmR91gCvj8+RB4w14svdLa4F0eT+RjsP40JDylJyYqjEUIhB8+vv+vo1wkISqLH2UAvrwxTHNmO7vPHgO0FKg+Q7unGXDreStNRbIhu8rlQSHMns4uJCTMfkLia6FvZFipHEFi3mQLNIcMEbH1uc+HLeu6bEtB5DXFmR5giSAWaf7RtODJUYAjHw1E3JYtC0kZSWsktzHYVZTIOV/l+rbaXcKTRxUeO/uLjYEt50ffoNBsKlfZjRxWS5xxRiKbz4fLkriT97KRAuEWIfAgyFs/ROJYsq18UVs25buszQEiKijxGPpUrAruDEUIl9E8NisdgSovQdH9RaznQqlzpHi+Labu7kQn8r96On+lKr7WnuDAISLmqyJERIe8X0U3wPfVg0uZCsSiQxIgsudLnlc3FxoSoZsVhBTh81YuRjjH9veb4xsuXKiTYWLOMarRp6JjELom93bw6cGCqxb2JAgRbTlrWqoalc6hKLgQc9rcIQJ7dEDDym0HJSaTEK6Z4w9bHrfgvqpna3G9Kk5/fKhSHGqeieLNYY/qYVAXKBrz3H0iCuNI4kN12sT7G2uesOf0t/U1LAdDpV+7VrODFUYt/EkCu4caFZanCXDk7NTZK6Bv6t5bcTWmjgGhlI18EArxZr2GeZjRpI74usIBJaFoJuYTHwfmkatmZJ5BKTNI5SY4tiMqmV+0i4kgKDpMzHeovxXQMnhkrsmxhi4JMHXQKWtQroi+1rnwJp8oQgWzg4GVtYCFI5D00TpX7gb8jaGlqqYQtwK4LcSHx9CM9yignv0j6QcJS0bbQkUJmwllGX+mtRaCzlLKztobKRsox2BSeGSgyZGPjk0Ra6SeBmMM+YkDSa3MFMwj7l279kfvEWtWW435z6oa1JIH/7kMsYtAQKNXz/6FrCd0cuNLREMfhekkasuQtxLQhaEzzTK3ecWLX0XCXF0u6+iYDDiaESQyYGzWKwrGZFjRGDq5QxIWlaFvcU9o1nH0mrdjFTo6bwHmn+5+fnG2GFQoSv/0Ah1Ifv/NCAbh16J7S2ggQ0t7ysq6xzgYSDpE3XoeyxHMvS4lYkxCwZre0hCX0LnBgqMWRi4LAMUCkojX5QSSjTcYt7iv8GtUlsR6okWoPYOgSswkr95yu2Wy+eO2Rg1hIJRCm1lBNDqeXAwdvlxI7j0GoF5Aj3mOWSgyEThhNDJQ6JGCxASwG15xbpe9pvuC9WWmlbq21yiwHXJGAJDa4Vt9JyXxIk3zkKVvKZa6u8cxQICdwNyIPTaLnmWAy57iBUmFL3JLVvJa19wImhEkMkhhpNhIQjlmCobTN1PU46GEhsXRkVr8c1T54iuy9tbsiapAbsMycBXhIkx+WoXStmwbV6fngdaTU+Hcd6V5oiIZHAkN+zE0MlhkIMOMhimkhsMEqbmVAbkqugRW0XtFCwTy2yjvjvkRSkHdEw62afkxVdajxdk0pTaDnv+663w4mBB45bZCvhtaR33ErYShaBNK/4eda+DhlODJUYCjHgoI2llmqDGyc0rnilgczT8/D8VPZHbFK09t/HAuFIQlibplUco6bPUr4+ukwkvzq/r8vLLzWdRqPRXgQRCkmrK65UaEpKhaZolAAtVyxgyEtmvMQYlBNDJYZCDDHhj//XJiFOaEnb5GSTOl9qW3Md8OM12hVq27wNvAdOEn0JUUlooLuOC3se7yCXnlQpldrE+yL/PgbbeSZXn9ortS2lRmsL4HLGEr+W5IbUdmPjfbTcP9aJury83Bo3Xfdl1fJ4PBato9gz6tMNVgsnhkrsihhyBgw/1/Jbi1mOEztHmFonCg8+a0QS+30spRVJY1caHgo9TtT8g/sy8LIi6Mfm5CBt6SkFfrlW3WccRRLQ2CcU5DVrVPjYQVLl26YSUoqKdC7f44Nn62mKmIRY2nhO3/qEE0MldkUM+xgw/JrSxI4hV/spFVhci8NJJ/m2Uei2FooSKfOURq1aKS0cQ/I6OTkJo9Fok03FCwxyAqBnF3OVca2+r8wr/rzJQiLtGsmD7r10FTkfq1yg835Z18Sk3ERkoeVYDDFicIvhhWCIFkNf10SBZOkHaq+0J8B0OlWtk1ItnqeiohCUyEZyQbQCXY+XtkbBhZomX1+BlUultF3pg+6w6XQajo6O1HUb+IypfepnaT5+ClrcRxLeOYswpWvwkuHavVhKWhBaK2XaXJasn33BiaESQ4kxSGhNJrntaYINJ1osWGwFN+kl7Zi33+rZUDvT6TSMRiNRICM5kLvLIuzRT0/aP49LkJYtVX3FD5b9QDKRFhPS82stmHicSqqImlO2RUJKiOP7slooKUHe6jlxyzd1D30Sx2CIYb1eh/l8HpbLZZjP5+Hp6ano3NJjq9UqzOfzMJ/Pw+3tbfT6iF0TgyUWQCjVvlqBtFL0jV9cXGz1t9R9FMKXmMd4PA6TyWRLMJIwJpdBX4XuJFKSiIEgkSU+n/F4/Mw3zxdsUbyBhLdERvw7cptwcohZJa3HTczdQ99z0swNRmtWZ8pqSSkk0hzry5KIWQx9WruIwRDD9fX15u/1eh1ub2+Lzi09Np/Pt/7Gc2PYBTGgpiUJe22A5hKDdWBaNSw067XBXKL98IAltwQwuIxCsAW48KH+a3tB0wprIiapfpDWptVnre1FsFhsF5bDQDV9sKy2Rgytymlr8ReygtAyKnlvMUEtKSBcsYjdm9R2a6vT0g4qIn2uVRkEMazX62eC+OzsLPvc0mOr1Wrreuv1OnRdF9brdbLvuyAGKesETe1Wpi7X4HgmBveVp9qRXCktNHbePloM6IKhSV9aeC91bXwm1ngABof7dgeEoBMlPR+pj5xg+iwimHpuWkxKu1d+HpKmtjbHYiH14b5JuVE1ItWyrVpiEMRwf38fbm5utr67uroKq9Uq69zSYyGEsFwuN9+vVqvQdZ3JnbQLYiAXyXQ6fZZX3RKo6aOGj2Rg0V553j0tCOKL7kqDzfhbXv2Uf1o9I7rOxcXFlvZMz0Zba6AJ232AW3KaUMZANr/n1oTGLSi6Jr4/izKigWeqYTs8OyzXsq59BpIVIx2n/qbOb4lBEMN8PheF9sPDQ9a5pcc4ZrPZs3M17IIY0CUkWQw5KNG+cn6D2igKfm6FtNBCebs8+NoqHRXvi7fN74VcR5xAhggkipjWju8Sg8MtMpgwc41rxpprLQc4JnBPa4q3jEajjcVp3XGwlVs05TLU/r8La3PQxIBavOXc0mOIp6encHV1pVoLv/76a/j8+fPm8+HDh70Qg6bdpAZnjfaVArpZpL2lybXD6+iXDnB0fVEbrfzhKWFJgoYTUqtYxq4hxW3oE7OEagPUOD7IouzL7UiaNr0nvI9cN6mkdNTOLatFkGPBl2IQxHB/fy/6/yVtPnZu6THE3d1dNLbw9u1bcYK0JgZ84eiGqRX8tVpN7rnYH54RU9uHGndUCjH3FFkDeC/7Ll7XCmiREnnTv1IJDxJktc8fLT8Uvlw7Lhm/0m5ruPlUjpCVzsE2a8ajNYaQG/MrwSCIQQsMS1p77NzSY4T5fL4hhaenJ/H6u7IY8IVrL19y95SamdqkqB142C5OxJz+xfoQEyI14BZDTGs+dDJA8HGE/8asp1rBxIPE/HlT+yXjMTZ/WigWuXNVA3eNWrID+3IrDYIYQnieSopun9VqtaXFx84tPbZcLjfWw9PTU7i/vzf1u68Yg8XHjytXuQ8/Ba0oHg7o1hp5iVad6gOfTK1IjP6vCcJWPvZDAXczcf98q+dA4xBTa3FMcwvAgpiW3yoGZZ2rMdcbKjkpa6AvS4EwGGJYr9dhNpuF5XIZZrPZlrZ+e3u7tc4gdm7JMUpPxY+WLsuxr5XPqIFr2hUHCmZ0lWgDuu/BZwEGO7vuefoiEkeN9iQFFXe5+Gvo4M+CXJw8nlPrUrO4ZXLKWUhtt3TBWC0BdMlp58bK5ZdctwaDIYZDRR/EYNHUaWBjTZ6UYEQiwYCpNjlKBl9qYue2KblwJK2yxkrgRemob+hvtz7jlwruCjw+Pn5WVBEtrKOjo6prxTJ/SlezIxnUClb6vbQGgbdN103Nt1T/cKyORqOtulqt4cRQiT6IAQdSbJDkLnahVa4nJyfVaa+pvvP88VxNje6PqotywY2TsiYtlU9aXJ+AllUflVgPDfRO6H1gtVRuhY7H42LBy10qdG1859Iug9b+t3B9oXJCY4WPGeoX9R0zr7T+xeYJH6tI0q3hxFCJXVgMmiYhTaAYaDKR9qu5RKzalDVDg9f6sbSNk4CsIt7fFua0tMgJn2usnEff6NtdUApJW8b3RRo9EURp++gqRFcM/p8UHesz0sZ97rNOWQWxdNZYOirdr1TpVrMY+lg86cRQidbEkAqUlZjDdB4KVxTgvA2rVp+j/ecMXKlKKO513ELzk1Zod133rOZRznNuAeldkzCksuJkNVk10b6APnG+Ih+fYYs9wTl5c0K3krZGDJbgcAzYH01Ri7mfCHi/1nuyzsMcODFUojUxpF5yiZCKaSpcO84RNNbAW25/uRBADZFXBC2ZDFrqJV+70DcZ8BRNaW8GqZ+kMWrlu0ejUfOFYhJwrMYC9UdHR8X9wLgFWiR47zkZYvyZSEpT7PzUs4ilmcaUMTyHKz2pPvTxnp0YKrELi6HPNlFIorluFbraROMB3BxNDGtDoVYs9bPkOcUWr9FnOp1mt9tHP2o+qH32YU3gu47tMUHWXqlGy+tuaSRU0n5MaYoFwTmk/cT5b3I0+1S8IXZ+CzgxVGLIG/WEkA64obsCNbLUIKOJQJMeg8I4iGtNdGyvlTZPq17H47Eo0Posa0HXJh9xbEVx7Ye7yfqs3cSFNVo7tUXfuJuK+9pLfeyxuYGkYCFWHKfShkT4Tsiis7ZnEfqt3UlODJXYBTHUaAMoUHGQYRYFd9VYrqOlkaYsiFzg7zELphb4XPBeLJM2FzGNOuaC6fPT2iIiIUtZZBR3aBX7QCFemzIc8/XzY1ZrJKbhaxspxZD77NxiGBh2QQwpbSA2KLhWRG1JZbVzUjKlDV5au7+4NlcqaLSAPq3epYqoLQkB+5/yv+cIdNxrIsdakD6tCg4ipL1D0AJtoeTUZomh+4iTDVq4l5mBfW3DJLQYzs/PTZlErS2AXDgxVGLXO7gRrD5ITYNHrStnstK53P1Rm/7HgRp9bdqo5s6S3FMtwIPblj2ekSj4YjpNQ+aC3VIRFq/DF6i1unfqE79vVEKImHJAZCulcua2w+cOjmN81kgiqfHMV2TH5kCNsrcLODFUYpcWg7ZQrGYAlvaFhAppg1yIaVkeEjSNnru4cq0aaofvOU3XytnaMed6XHNPBWZjwj8XqWJ0sc94PG7yDGiMSNV0OXnlXo+TfAtXJY0xKbuJj8PYPOIWQywtNRX72zecGCqxyxhDyUKx1pqHlnLHhR8JJAsxWCyeEsHJSYFrhi0XrqGASX2Oj483hNG3UCABZCWHFivhJVLXdo2rsf4wQFyzMp3X4+J9ovFpvYbVymittLWEE0MldpmV1ELI17QR03K4ALZkNml9ahXA5m4CFEYtymlYrsU/u0iDlUDvzmJBtHgmKEylVGO04HLvQ3MDkYDNHS80Fsbj8da+JzyOZn131ljYvt1FMTgxVGKf6aolA4trKTkaOU5EXiMGJ3+O1kmTCAuCcTdZiWaPWUz4OT8/3wSdW2js6D44hD0b0E8fc3HVppji2ORZOfw51dwHz1JCK/H09DSrLXQXSURvTWMesiVghRNDJVoQQ0rAa8dLBiDl0fPN7C3CFwOLqHFzLTxnzYIkTHGi55rxBIsGX7NegbKaUBPWXDZDIQUO7RnxNS21kEqclFiWBElZ4O5Wug8r6HdYboQsBspcs66X2IUl0Pc1nBgq0YIYNC2eZ0ZwoV1jMaBGnmsxaL8r6Y+U9orPoTRAR4IbC/BxrbVGYEskILnThugmIHBy68vawTFHCwu1d57THgWLMWNtMpls1dWygscuWilipUjNr5hMaBHUdmKoRCuLQaumSsdbaQe1wrallsIzeHCyxwZ/yXUuLy+rKn5iW5aVyqUrvfcBS4A6NxYgtU+WJrndSlct0/vkpMZ3Pstxh+G80LKJWi7YS0Gy5LliJs1jrviVwomhEq1iDJI2vgtts/RaLfqIg5hcCpghgqtoS4BWw3Q6zd4SUkJKgB7yRj6pe2v1vltlg3G3EQlJKWaQ057m3hqSxaD1p9WiRSeGSrQihtZWgdaW1U2VgpZPbhWM6MagQYyBw1I3QyxrpYUmFROe+8o4aoXYM2vx7OidU4yrRWkLag8DztwCzc1O0qy92jnah8WNVkNpooYEJ4ZK1BIDN7MtAjXlCkLNjE++Vm4qrl1JwegYUADR/fBUxpxyyrxduifctL5mxSwSnyQwWwRphwIt7tBCmNP4oIBurdsNiZpvN5pr6bQS3JLAxjmCc6M2JoDjHedUrbvLiaEStcTANbTxeBw1BS0+RNQctOyNFoMfBztPX021L2U4Sf7iXFAGDK5wLsmW4veKz5ELzlYrhocCKfZDwhytvBJgdVlUDEpdIOh6fPPmzbN3Q//fhfuHwOco/R9jGFxRs/QxVRrHLYYBoZXFYHWhWDWMXNdOKUqvwzV7KSc9F1JaJCevEnBftlQD6KVBsox4Blmp5UXCCxUDHP85z5PeDVmX/Pe7DBgTNItBSqigMWtRLngtJu3aLea7E0MlWscYptNpVfAIA66tVl5K59XGKrim02IwS1YH7j1d+jxjPvcWJSQsfeiT3GPXjd17SXYXjROetorEkDP2F4vFs7pdUpZOKXn3HRewCHuCtBi0LzgxVKLvBW6okVsEHHdNaSt8c8xOKVuDT7iaCcTbL62LJAltdFeUQEtNLYl/WEHulphQ5u6SvogjVbE195o0TrjVxVeOl1gN0jvOGZeS5drCCokpUbhboRU5ZFIKJ4ZK9LHADcELfNGA0AY6ZvvE3B10TYsvH1eFovunlTDiMYYSPynesyTMS4khJhRbI7ZC2Pop3e40hpjlUKqJpxbW5SoFLcYiKlW4YI4rZLmCWcvgkxSsGNCrEIs1tIATQyX6thh4eQKtAqTUZqwmUI4bidrRFiPVDkpuUpdYDPicpOym0r7hFpVcw22FXAvB+mnZR4lsa3LlU2sm9hHM12JdqEShZWO1GPBe+W9y5o6mzOH8OT09beJmcmKoRIvgc2xgSPvdpgaSdk5JHAMHNVafxDZzNHwtXlETfNaspBoNmmdd8TZbIWcTn5LP+fl5k34uFguz1WptL9bvFu3XKisYA5O258yJg+B9lfSR5oi0aA8tHXxHNc/PiaESrdJVY0I1d5BLbfKJjWQTg1YTiF8rt1a9Vvoix7TmbWJtpNpsLGyzlfWBSAlG/PA1GJTauWvLgY8fShEtfc48fkHt5LoSpTgcxgZKMpM09w+6+2rdPyjQU32jc6UyH9xioDZrKuQ6MVSihhh4hgI/VvJStTZxEJLWa7EYSBvHrSAphz3W/1h7qH1xK6jkvnlAE4urlT5DzbWDpFiCWNs1wlxzeSG51BJaLBBduq4B002pf6WKEI4rdKlgP3P7pVnJ1v5pZTa4lZuKV1g8C9RP+rtmDYcTQyVqiIFrx1JQKpf1NY0bYwW5/kfeF9J+yO+co0WTxoW1bah/tVUhWwSxpXa6rquu20SIkUKt8KZnGCOI2v5LLhX61LjsWriNJIuB99cKEq70vrBAX05fpcSNEJ4rarQVa8sV0G4x7BG1FgO+OHy53MdtNVulwVUz+aQ2JeFgFb5c++b7QuQKcrw3CuKenJxsNqIpmWhaimpN6QbqaysLIYaYZl9r8cTaLtFMES1IAjdOwvIbuVYNEiwlXpQEeJG08N74Sm+8nvYcYy4xyequcfM5MVSilBikl6cFZq0vVyKWGuuDt0ngwiFH+EouD3JVjcfjbEGOGhmSqGY5pRAT3jXEoAnUPhfIadesJSEtvlFijZBQxIw7cgOWZNeg0kJZdCUCkgtqGk/83q0VCPg5fHxKLjXt3naxk5wTQyVKiYFeWssN6XEy1VofIcjWQiyYZgH1i9fLKXFHcCGO9ZpKzXJu4tOnxs0TI5u+oV23Zn8FzV1VQpza86ZP7iIuDA6T+6ek3Dof+/T/3AQOnOcxhS/lSpVcwSmlscYCc2KoREuLgR/LFZJIMrXWRwhx11aJ9SEN/hp3BBcqJLxraiPRfZaWaJCguab6shQQMddPKTjRUZCz9HljGmatdYOkdXFx8axd6zOXLBYU8pjoEGszJ908puFLx2osghScGCrRqlYSQnrhKaFu1UxyIFkgRAi51kcIcgogTdzz8/Mioolpm7mCCp8V1wxLJ5/Wx10WddPIobZN8uXXJg2EsB2UR2Ge+9yRGKRYmHVdi+S24X78nBI1fHc56b5iczVmvfehYDgxVKIPYpBeeEo7sPoyrYMpZYHUZGhIeziUBsZjWT6lbiQpX7908mnWwq4hxXZqt+vEdEt6ZiVtxgg+l0DJlXRycrK13iO3bhbWMLLOLe3etEqruci1ENyVtEfsihhKXzL/nXVwSRZIzUBDHym6HmpiLBh45nWGSgKsWq2immBtC8JqBUmDLr037sarIb1YW7mkTAQwGo3CJbhqkSQsxMBdptIYtc4H65xrHTOocTU5MVRiV64kQu7gKNX0U1ZL6SClyVZTAoPAF7TVCjvNz11TxVISdPuCppm3bKukPe5S5BZbjkuQ7xmBcyhnzCIxaL+rmUsSagR5zXUlDIYY1ut1mM/nYblchvl8Hp6enorO7eNYDLuyGAi5g6dW20FoWU8WcJO6hmSoPV44L3fzdw5t8VZpPEASnruMLUho6dZqFVTnz4nvxJZDDNpKY7qOdZxZzh2SIG+NwRDD9fX15u/1eh1ub2+Lzu3jWAx9EAOBB7vQNC61GOi7EheOJMxLNX6JZHJiDfQbXqqjprKkJORq/PAthXArnJ+fN7NgNCLNDfpzVxKNyZK1DC2VqpprHToGQQzr9XpLMIcQwtnZWfa5fRxLoU9isPg5ragVxCF80fppRTG1WdIvvB+e0mftC01KXhumNJjdMugcwjCJgQvhEg2foLmScusmoaKC7+4lCHJJuSuN0eFv+76XQRDD/f19uLm52fru6uoqrFarrHP7OJbCroihdiDwiVfaFgoWaTJbQcFo3EO5VKBLwq4kp17aF6Fm45tagdkHJPdWab/QlURKQ8mz5wpHqfBDd2VtdV2pTyXzBpUwTPEujdFJ6eJ9kMQgiGE+n4uC+eHhIevcPo5x/Prrr+Hz58+bz4cPH6qIIfZSS33v0m9ameU1k4+3h5vClwxsdGe1IAbUXCXXRi5axxdaWR68X6X7S+Bzryl1zom9hYKAO7C16NPlZVmZFRxT2Jfctug30gLT1pZVCAMnhuVymXVuH8c43r59K5rPpUX0rK4YC0lY2ivNzKhph0CWAWmouEiqZPJie5SOSJ8SISxtmlNTCqMVwUjttWqn68o38+HCczweFxM8J+QSK1lSWjA1Ojf4bLEYrPPS+rtUe7l9KLUmBkEM9/f3oo9f0thj5/ZxjKOlxYATy7pRR0y4WNpD0zY1WLggL22HQEFirqGWuiDQ7UMrcGuEsET4LVNVa039voihtD3JLVX67LEce9dt19LKbQ+FIY7hUi095toq1dY1gV1qlWjCv7R/gyAGLfgrpYzGzu3jWAq1ZbetFoPFFWS1GKzXJEE7Ho+fbeyTGzBeLBYbnzRNVlpVmrPwCNuThBIJlxIhLLXVMitpKNkrvF+lmUlaqY0Saw0tBh7nyX1uqGjExnAKXPmRFBirRs7P0xSrEoshpqQdtMUQwvN0UXTtrFarsF6vTef2cSyGVns+a/swc9+ihUCOjo42m4xIO09ZfMFoSkvaYK7FwN0OkmZeEtSTPqXun74thvF4XNxWK0iEWhpj0MqR1FgMLUgLV8TTWMcV0VZwJQqJITdmJ9UJq0kZ5xUJWlZpDmFAxLBer8NsNgvL5TLMZrMtbf329jbM53PTuX0ci6FVdVWcZPhyaSCkgmicQLjAzSEXvC4XbFTSOncTnJiGT1ZEbpAdF7jxWkAlpr20YKsmYCytpN43pJLZLdd9lFpGWinv3PbQfUREgP/PzcBCApCSHayKQ8rasBKNRgLcMijdz4IwGGI4VLTajwFTNzHQZc30oPZoIlC7FNglDcoqzBcLucIoJ4ySDdaphAXXEkvWVaAbiqdP5kCzQGpKbrcSmi0hafmlkLT80owiXsqilORjBJM7/i9Z8Jq+R7eXdfxjYFy6fs4mPBb3UO6mPhxODJWosRj4YjGLDzLWHg4GnueMqaHWySGVF0BfK7Vn9WFK2g6fyFYhQP3AbRtLYhXYnlYrqaS9EGRi2PdaBqlPpUCi1yxeKzS3VO4iQ+wTbe86Ho+zg9lcUaCMKxyvOVt7SjEKVP5y92NPkYNbDHtGTYwBhT9385QIc/wNdz1xTccyOTQth8x1XJyWK9CpzdKUVZy4+PxyiE/qG8V7uFAogZT+WiOIayEFi0sDz5JrkCzB3GfPrdPpdLq1H3gO6N1dXFyI79KamMDnCyodnCxi7dE4ReuD7gvnt0UBiZEMHq9d3BeCE0M1WmUlSTGEEmGuaQopC0UDDmz+HZbYyDHPkQSl9i1Acnrz5k21W4og7VtQKsy5vzvX/dAa0n2VCg+uUZeSArbFA6qWjXA4NKuvxPrDPuB2sTmbOGEb6JbCe+Yptdr90m8mk8kWmWjxxRLLjeDEUIlWWUmxoHLOSmNN0PLvrQJZ0kJKsiAkEqzZm5nfQ6sVy5pAqUErLb0GkoZfkyWlpaqWPHOcA5KQy2mThCxp87Q3ea5FKgWdueAlq8SqXKHyJ7mRUtY3yQIiBb4r3GQy2XInu8WwR7SolRQjBzxmEeZaW1wDswSwEPzapb9HjZlPmJoSG1Tn6PT0dCMMSjRzzdddM8mk9mo2/ymBpEnX3JPmJit14dHY5GsQcjPfahQNBCoa2KdcCz6EuBXPr2NV/LQNtChdvaZESQhODNVoQQwxgY/HcoSxdG6O20YjktLBJlkZkvun1gQOIb5i2/rbGheEtc1duZSk65+cnFS12cqy4vGFmhhRbHznkgZaRFzTz2knZl2jBWINPKdkACpgNfPJiaESfVsMpdkF0iTJIRYcYNhGKVFJ10dtKTdeobWJE7pk4ZYmxGuyiSQ3TgsXlQVazKRWo251PzjOcLzljq0Q4uNbG8+l7VnBNXwpgGzdYMriWUD3rFsMe0TLsts1Wn7MVVQ6yTC4y/uIezFYTGAJ6P4p7ackWOhTUs5Cy4Ov3ZtB88m3ENIaNFKoWc0dglyevPR5c9Is3SGQ2sLxw+eAVdO3jkPLefy6kmJlFeCxmmJ4rAWhOTFUoiUx1Gj5MY2oZJJJv5N8wVKqbSlKhYEmcEuCvNzv28pqCEEv+9B1bWMOi4W+JqOWiGLPu+T9LxZfVpzjAkVrYBfbkfYcaTHuU9dLtV9jZSO0dRBofWDGU818dGKoRG26qqbhlLSlaUSl1+HnIflIGRa1mi8Jzlwh2VoI0v1IGnctYpZDrVVCVp7Wfov+axZVaRZMbKVyDrjlSGnZpesqUuOZu4hi7bRYV0DtpayPVvPRiaESNSufYxpHax8rHufaRY5W1iLrQ+pnSf0ZAhcuSBQ1WlNfbp+Y5UCac+51Um12XZtgt1RPqpQUFovnZVdKrT1OiqUlu2Pt5ypXuZp7bputFDIJTgyVaFErKRYsyxnYmitKSgVEgqidPLkDFCcwxTDQfZPrqyaft7Rncw2J9SlgLYKcC18MKNLqYOtvWwiP1oHsmLVQ4gbEMYQpy7njU9L0S92c1utaXVN9kgHCiaEStdVVU9q9FkiTfit9r7l/LH2wInfSSK4fJItSi4EWChHh8AJouYgJ2xYTM+ZWavVpGbvQrlEKXtcIia4knoOuSFS8rAoQt6hRUPctkGOK4q7IAOHEUImWMYbYeThIcwRxK/dPDLmxAS2rhSps5gozviiK/q5ZzxBCXKtv6aKQ3DO1n9bF+mJB59L2Wixqo7a4BYxuKqvlSPMKS1jUpn3G+ovtxWRBibWSajMFJ4ZKtCiix1+4FmSScqFrUBqE5t/TBLRq+kh0XfclA8VaM0ZqT6ppQ+3VBHU1d01NOW4NLSyImh3ncvp1fHxcVb0TyUAqP5HzbKktqqJK7watUKsSJbmOWgvm3PZK442+wG2P6MNiQPdPn2asdYDSeTTpiLD4Ah2rpo/3h9dH10LugJYmMhY8K3UnxbJ8+iyhLVlVmguubxeDdO+l26eGsG19kDZfI4RRAZD6mqsYEKlQ/Ksk8ByCvsnWLlxDMdeUBU4MlWi5joGwC/cPXcdqMfBJRwSBE9IqfOn+qF4+3WdpCW5+L9gO7tdQAmpXC7zu2ve7a8SsmFJ3GlcMqAxGaVFFLJwnBbRzyTOlmKXIiwevWyV5xK7VOi7hxFCJPohhV8gptyFNOO6yydHKUaOh37YK9nHB00K71xa91VgjQ0cqc6pU6FCht/F4vFU+utT1wWNJUr9z2uQWA1rMlpgDEgdX8krdQtpvSt1cKTgxVGIfxFCrDZTEBtB1xN09JfWcUKufTPRSwrkDXtNwa1wfsXZLXBWHgFiwOWenMQkoyGkc0Hclrg+0XHHsUGZaTsE7VEpQuOfsUKe1E0JdGrq07qgvt5QTQyX2QQyWwRUbMBisyxHoqImju6ekTgveg+QPLS0eyDV7XN9g2XErhlj2UB9ugn0hFQQvBQ8IkwuJSKL03fBAa6mwxL7wOmEliQzc8qDvsG+xvu7SJcXhxFCJoVoMKMTPz8+3hCwNsJOTE7PwJQKQKkGW1GnBe5AmS6lbgQu1o6OjZ/GRUtfPYhGvQbSLQHCf4Jak9KlZF8ED+XidksJ+sTFU0ha9W75yHo+VuEtxHOfEK/ixXQStCU4MldgHMVi0ackVQGWoS1xJ3JxFP2tpuqGmPWG2U0kAHl1etENWC2KgPuKOW7ztklLfQ0Gs5HiLRIjUSu1cxIrK5YJbwxQXwHhFrjtSSiLJEfYp4uuTKJwYKrFLYsgV6DyThmrwSAI9BX6u5ArK1e757zj54ArU0mAdCvLxeNxEyPFN3PkntzLovsFdPK1dZTyz6/z8PBwdHT17frng7p0aDVuzUtE1WVJiRVN+cncZlOZY6byzwImhEq2IwTKIc2MDNHDQDMbBlDNxtEknpRha2035W9EtVOJSwqAmtcG1zBJIbUufQwhKp9xjJcTMQWOHx39o/OTGkmic8Hdb65OX+lKS4RQLPtPxHBeaFKvg12kNJ4ZKtCIGGkAU4JS0ztxAHZ8o3CTNIQlp0mhaVqkmw8kGLaPcwY9uCySuFsSA/Y0J1KGTg4UUaq0r3IyJ0lRL3DIItCxxxTsf17nVg5HAuOvHWq4D54T2G+7qtPYrl5Rq4MRQidYWA5+YmH/PNaWUcJOIAYHCPiXMY6YsF9y125HyT0lqJM9tJ1KguEArv7lEDq12N+sD6FpLkVqLld14HVzNXRvn0Xz3pfWWqF2u6OQKWpwTPCmDFLPcQLZ0v1q/SpUyDieGSrSOMXDTVQpmWrVebsLzwZJjMViDX5oVkcJisV1MDrOJSiwGLfujD+G3WDwv2KbtW1C7BqC0f0SKKULAleItLB1eNZf+rS1Mx4mfu/esYy/lzsyFlIiBFkLNeg2ERgBuMQwErYmBB+ooPTCVoRBriyYhn4ylmn0MmhWRAl9ZTa4HFLg1k1VKW6VPC0h+75jwbTF5rZBWrUuf2r2gOUjJoeyy6XS6teisRrNFl1QI8UVgMbTSsHl7OP65VVlS9iVGYJJFUQsnhkq0JgZNy8/R7gma1iy1qf0m1Sb/LuW+0iAJL65t5U56QqycRUvtWLuP1Ke0qqgEIvvz8/NnlpeFsFpAKgCIliQ989z7Rcvn5ORky/oiRSrHZacJ1Jp3oVnMWsJCjJBSQWx+vZYE58RQib4sBi5YY4PE6m+0aBkpDUo6bv0udd84YSg4KfnEc/3TfIe4ruueuXpa1TtKWQ5aNhOuzj4+Pg4XFxdbQdsUWWPNqhJiamk58md9cXHxjDBLnre0ZoRIAr+3QrMyShQmPC9GNjgWU1Y19k+zCOic0jU/GpwYKtGSGKzuotiagth5iJgGUmMx4LWxNr4F2m+or7VCHDNwRqPRVjyndfYQdwmipmzV4CXyoH7z0hIl7XVd2x3e6L5jBFVTlgSfm3bP5+fn5n4SmfLV/NJ4tio6FrKxkAwSjBTEJvThDg7BiaEaLYkBB09sIJZaDISYT7LUjOakVmveStZN7Y5avFwG195LXVWW61rWPVCfSoV86nN0dCRao7WQtGFMN25B6iFs7xSIFVrRYrDGSWhucO1d8+VbnxkPOKfmqNaGVPtJarNmM54YnBgqsQuLIRZ4sraHyCEdK3Ci1ZS04JNAEtalfQxBDnK3CoimQASXKuFN8YFca8JCDH2AW3T03qXaSDUKR0xxytWcNYHP2y0dazVxCnye0oI2WjVNhFhbOViCE0MldlESw2ohWBFzB7XQxmtyyeleKe1T0pJqMjD4fWL76AbpuyAeCTJeSC5V2TRW4VWLM9C+B63dDZKlgCSFKbKlzxPHPr77FlZPyi1qHWu181FqL2YJcMWipjKtBieGSvQVfI758DUtRhuglu9rNKPJRN9PoaTOUayEco21wK/DiXCxeJ5W2Cc5UD9Q+HDtGz9IZNrxWoLPAb4LyU3XQmjzmkKl40oCXwuB182ZE63GJMISY6RifyVb4abgxFCJvtJVtZcc02K031q+j7muYuCCCuvOl65B0Mpv8DhDjeXAFwnyuAhq4buE5m7i61koqE0B6X2U/CYyoOujBltLTEgA3GKQSr+XgK+FwOtKSomEPuJ1qd9K5NVaiXFiqERrYkj5SmPEUWMxWK/BIaUQ8rpEuZNYuj5q0pfgYy4V3pwYYpr4dDpt7i44ZKDLjVsJ9Exrg6I0ftANRu+8VZKApoDktB0jEm0eWa4Ru9dc8irBIIhhvV6H+XwelstlmM/n4enpqejc0mOr1SrM5/Mwn8/D7e1t9Pocu7QYav3rVuRcBwUq3/6Qt1MTNOf+ZV5moOQe8RqphWnWMiQvFfj8tXUYLYofIrgywN9ZrlCUyIYfpxiOJY2XW7HWTCSL4kW/lRIidqGkDIIYrq+vN3+v1+twe3tbdG7psfl8vvU3npvCLmIMhNSAajVgUsEv7Vwt6EyaGU26GuGK18O6PrXgPvLxeLxpn/YPIIvoNVoNKKTxw0miVS0gXk5DUixyLNsQZPKndiWysz4THP8W90/OPN2Xpbp3Yliv188E8dnZWfa5pcdWq9XW9dbrdei6LqzXa1P/97FRj+YSKkm3lLQw68rMVL8I3F9eunIZF4tRDjsJo9oJpAk+6bOLoPTQgGUuKE1SErSt6kBxAY3j8vT0tKjekFZ2RXr3ltIasUB/LokNzVW5d2K4v78PNzc3W99dXV2F1WqVdW7psRBCWC6Xm+9Xq1Xous7sTtrH1p4cNb5XmiyU8kb/rxF+fJBjFkXJ4NeCwtRXCsDW+rQv/50FwwUdPQ9My8S1D0OZzH0ChSelEkub2NRYb6S54/sklw5q5zlb0iJ45VN0TfIV6aXjiMDjACnXbK710zf2Tgzz+VwU2g8PD1nnlh7jmM1mz85F/Prrr+Hz58+bz4cPH3ohBmtWQurcFHAxGpJLTVZJa+1IshguLi7ETKIWufp8LQMnuZbCcIggAX10dLQhc04Cl5eX1UXxOLjmjsoJCvWcxWxSBg+NS+3/OUqRJbHDIvTdYmDQhDZq8ZZzS48hnp6ewtXVVdRaePv2reheaE0MscEUO5YT4I2VmSjRYCTTWupPjXbEhZMUiK6dYLFniFYKt2B2sXagbywWC3EVNSkOuO8xlaWoTZXFccM1d7qedXMqDq654/iMFaq0IrUWwtJu7fE+0Bsx3N/fh9lspn5IW7+/vxf9/5I2Hzu39Bji7u4uGVsYosWAsArdlDC1mL9amzh5eX9K2kVIVVEvL79UM+WWT4vJhBpqLHspd2euIYALS+3+6PlSUL7luglpzGp1pXKfbWvFhCO1FsICibys/e2LNPZuMWiBYUlrj51beowwn883pPD09HRQMQaEVfBKaabchEbt2DLwJGKozSTR7o+nCHIXWMu6R3znMY0YWgbDdwVuAdFzRQKmVNTLy8tmfnh8PpJrSCKo3GQIy5a2tWixFgLnbG56akuSQ+ydGEJ4nkqKbp/VarWlxcfOLT22XC431sPT01O4v783932IxGBJNUU3APqP8Td8gVpqwFsmXIv4SMw90CLuwmG1GNDdQX0kAXt+fr71bPZBHHhN+ptr5uPxWIzd4BaaLSwGej6xvZF5DMN6PWxbmgc1wpQrJrkWvXYe9SnX0n2xFkMIvwnq2WwWlstlmM1mW9r67e3t1jqD2Lklxyg9FT9auqyEfRGDNCBwMlk3AcGJyX/DrY/UhNKOx4R36f7QFnO7NhtKQowYcFJLmwPRB1fxttb0pGfN15Gg24tSP/F76R5bV2bFd0994Ps1lwq9lLJQGguSCFMa67hXRExJkhScWBA79l1rDIIYDhn7IgZJw0Bhb9FWrAtziBx4QTPpXKkdTYOj77uufsMYSRvmm7m0EMJ4HV4fSCIqqRoqpr5SITR8vkhkuUIAnyndc6zMNwZP6dpS8LlVZVbJhUR9wMB2C1g3tLIC3W6np6dbVgPNEz7mYim1lncr9bkvpQLhxFCJfVoMvKBYiW/TOjBjprn1OjwbBP31pRvRx6wHynA5OTnppfT0YrHY6j8nJy2ASmRAQkYT2lzAoDClNk5PTzfkQ7u7IRHQbm+Uesrb5cIY33cfxflQSUAXUh8bzqD1HEK99s1jaCig+XOjf1PjLXV9yfW6i7I4TgyV2GeMQbMQWpqaOBBb1cBHIVAbCJQC3tgm16BLCUiDdC0S3JqQPzo6ehbsJQFPcR/pdyhAUwSCwomeNfaJu5CIAPp0U3CyxIBrKsU59zp0v6m4To72HYuP4TzJSa3N1f53YS2E4MRQjV0Tg0VQ1wyeksBZCtICNXJPcGLIFQqcGHicRVud25fg0wQ3rublQjyE5+42jRhwS0tp/2i6N/wes114e/j/voSNZj1psZaa8YtWKP89tcuD3NbxZrU4cgou0rPRivrVWDg1cGKoxK6Jgbt2UnGB3AGEkzI2CHMmrybsLsH/TfsR56aZ8j5yCwFTcdH/m7q/UqC2iu9JE+AEKQ0W+4vF5Pg9dt120T9uiaUyqfrcz0EiS8n9iftrlFimUlyJC9RcjV4bWykSyx1X3OUVa3tXcGKoxK6JgdIFU6WNc9YgSIFbFEDSoM+ZvGgx8OAqz5ipXbm8WCyeBQC1lbo5ml1pX/A+Je2czpGEN67VQMKMWSZEEpIGi+TUMlNLumdpwRz53Pl44C6/XEHIn4fk28exbH3vfPxLY7+FcqG5U3dlHUhwYqjEviwG9DdLWUV8DUIM2kSp9clKIE2aZyFRH1qUtOBCltrkExpTNncxASVLjp6nFlcgISf5szViaEHmpaD70TbtkSwjUm5qrVx8Zvw+eQE9y3PQztunJr8rODFUYh8WAy5Mk/K/Q5CFkDbQc32iNcIFJ7BkqWD6ZOl1eFsogOhfXLuhrePYBWIWA/YbQYKJW0aam2qXkKwBfM+SpVOaqswJj88HzLjKtQ5j47xmDlhihISWK7Rz4cRQiX3HGDSLQULMMthFClwIX2odjcdjUfOiPrbY13exWGwJIFx0pq2yzV152gpcyGFBOc3FUOuX7wPcYuDKirSuo/Qda5o7TzhAsrDGFVKxrhJy4MSopehyi7D14kILnBgqMcSsJO382OSImceaK6LERZESaq39uDyuoLXVOtbxGiA9Swy+88qr3OWT+4z59bR3KVkl1mtRHyUFIRV/SwFjPTTuUSHTLN19lHV3YqjE0GolceBkTK1w1o5pWmBMO7T26/LyMjnJSiZh7u/pHClw7dgGF2CYLSUFnem5k5JSusERvsfYeJXiGLkuJC60c66vQbLYpXnAs9HclXSAGDoxtHATpSwGWp+Q4yuW4gua5VEb17D8HvuAGie5vFploBw6tDgBbfeJLrnW7kkcyzxLS9LuW5AQJ7bW4zA2D0r3oGgBJ4ZKDIkYcgZuy0Gu5WHnQtPsW8dAUi6JxeJ5ymuOn/qlQBJasSD55eXz1catQWMErd9ai1JCzpizziXtPO37vtOpY3BiqMSQiCFngpROJik41yp7QpsgqL21CA5Te1jfCNsmNwJflHaIG/HUQHJz4DO6vPyyQPHk5OTZb1qU6OZjgqed1lgFtX0hxBQaPF87z/r7XcKJoRJDIgZrcC51DMEnopbVU6PRp/qJk79FaineB11XWkdwcnKyKT8xmTyvLvvSXUv47jXypwwgWpmN7kB0OZVq81xoSmTVylpIERleW7KmUoSRazHsE04Mldg3MViCxjWaiJT7H6vImZoEEmLWCz/GK7K20kSlNQFoMeD9cyFFtW72uX6gBhbBxt0a9JtYsL6VxSBl7rRQSDhiRMavFRuz+JtaxWxfcGKoxL6JIeaHzDV9JUg1d7AmD12HtHqu2edeg/efT3zu426hKYawvccwCTuKK9D9cjcW9RsJhKqWDqm8AQd/rrif8+Xl9loKblFR6iQqAyXZPzn95GtaWghcHuildSN0j6k5U+s+zZmD+4ATQyWGTAwaSoVUyqeLwiInFqBNEvoe28LSGVpVyhJIwXT6aHWFJIuBWxnSJjG7JAkuAFHbpf5IgXbtudM404r+1d4T9U9apJby1aeOSeeldmTj/bL0wfJ+h6QoSHBiqMS+iWHXy+a55ohATbRk0Rs/l1wRSDZ4vzipLXEVK2hRHAola78pbZP/Vtu3YbHY3nyHt1d7P1yzR0FP5I5CmFw/WkqotE6h761TuSIScyeVWAzoLrWmTEvXoe8s1rITwwvHvomhlUlqHah8P1vtN5K2nwsu1Gj1Kt2vVHu/ZYqfhXS5m4UHsXnQnlsll5fPt95Ed5SkqfMaRHQebvaDAXIUfLF3hW40bR1Ci/cqPUNumYzH42dEGrMQcuZATMgjaeP4toyrnGfTat72BSeGSuybGHLMVkmQ52g5ITwP/mq/IQEmFfizglsMJPBiGmKfud+SX1raZwFdXJpgRk0X93pGYuHbX4awnQyAQkxy62C/cXtPLhSxfhW3CPreHZCPPYmMpcKRSHy5AeiYUMbxzUkyNa5Sz4ZbOm4xvGDsmxgs4AIEJwQd04KmEiSi0XywPK01ZzJgGzhJY21LgdNW4M9Rq8mDQgQ12tgeASE899tjrj4JPnQ7oduOWyq4D8disXjWT+obd93gb+iZt4gdxJ5narGalEpcs5gyNg55ggX9nesajVlbQ7USEE4MlTgEYrBYDLmrelPBOGny5U4MycUgbV4zmTzfTD7HrLdCIh/SxCm2MJlMtjRw7jbi7aEAQaLBon8xQuf3ix/pGFkM+Dwldx0+876EmRRbkcYNt8roPiywxAYQUowj9/7xeePvhm4lIJwYKnEIxCCh1gVjIYLUNUv6KxVuk9we/FhLgkBIBMTjMLH7xmC1tL8GtY9be9L9YYowpgpbFuLx58PdJ0gOrdYJaM8uRpyXEIM5Pj7OLtioEaj2e+76K7E8S57b0EjDiaESh0YMWkwhd2DGzu9zkEtZUZa+tFgxLUG6No/DxMC1Ui5UYjEgbiWQO9AS10Hrp1YQloL6oNVX4veHykDOOI0tkJPa6mv8xtpNEdau4cRQiUMjBknDDaFMy0ldQxrktdfJsWy4ENjlZjyxjCZJKHddt9n1TIpJaFq/5PYgv/jJyYn6nPEd9UnkFlcOWg78/lIZVRZo47GlMLas8k7NC7cYXhCGTgy83g1ld3DNELUzHLglA9aiGZUKpZzz+UQkUsKVzPsACkJNsJNbRwNq27E29iWItNiVREhIkLyvLcaHZCVI+y3UwFIXKjfesU84MVRi6MQg+Wcx1ZELTb44raVWhX5xdAv0ZULHtFMeaNXO7wNcOOF2l3wXOQ343IgcptPps9hDX/GB1L1pC9W0sdXq2UvjCdvmpFwz7rDdmMWgkVWLPvQFJ4ZKDJ0YtAqZMf8qX4xlmbCW87RJ21eAU5uQVLqCLAdyU/S5BiKnr5ZnOTRtkxMCZopJ1ktf7z2mDHArpeT5aSSTO+6tv90XnBgqMXRikGBx9eQOWIvmL5nS6NppWS0V+0QxBamcAwqsfRLDoYE/b/6skeAo8wrXVsQs1haCMpdorcB+SyvvU32JfTckODFU4tCIITWYS/2gqfIRKbeOZYLFIBEa12KlbSdji/Ry7v+1QRLs2nOShD4/v8WaiVTMItXHXMsYx1zMAjnEMeTEUIlDIwae9ZIatJIZHrM0tEmtaZQooGs0Rslni/528tvTKlZNQ9U04VJh9VLBFYGU0KX1ETyuQyBiyF1lLZGBlskUe5elFi8nCfo/ji1JaRk6nBgqcWjEwC2G1ISQqoJaLA1+jCaKtECNC2EiLUvbsT5IFgm/b05WMU0Yg7yvATmKgOSG4xZbzDIs1aqpH1i+Q1MyYoK9NEaGbfMV+UgUQw40S3BiqMShEUMIeRPCajHEgNYCxRN4pVTqF2bnSH3IvU+6HlUdpQ/WGsLFYzHhgYLNsk7hULTDEOQ+Y9kM0sqlcYMWARIDvTfK1qF30ar2kqZw4LW1MVMirFPKD7VH94rrSA5tTDgxVOIQiYFgiSdoAzpnoEuao7aKVauuWbrnBL8GapaWe0EBgxYDCbnRaBT9zZAQUwioz+juIGGLyQESgUq/DUGuVFoS3Nf6HVNarAqP5N4pEeLYHrrFDjWZwYmhEodMDFyA5Qi0XOGH/n6eY58iH6kcRClI2NG6jhTZaH3D0hySv10rw5CD1lomvjPJWuMaNMVlaN1L7hoEOgc/JX52rd9anIjuJxa/ynGTWd4Dxjo0K+aQ4MRQiUMmhlwtK/bbFKTJzf/PS1VI39cKS/o9asGl7UjrIWL3Gwu8xwRrC0KkPseuyZ8tWkXWsSJp9y0SC7T3Lz0jadtRLbaRek6cBGP3zvtzaO4jhBNDJQ6ZGFJoObD55MYV0Hx9QWwCaiSSi1z3FAo5ui5p2KPRyGQx4P0hOUpptJw86D5j2UDcKpMEcQ7ZoFVkFXitySyEODFIfcHMOyzZHatjFLsPCxlJfT1kDIIY1ut1mM/nYblchvl8Hp6enorOLT2GmM1m0etzvGRiSE3ymClvAc9kkYSvdE2JRFogpgXzVbxEaFR7SiIv3jdJiKOrA4UY/Za3hamY/DiSDQpGbhVIpCw9a22l/K4Fo3SPsfeO5SlidZhiiN3HSxH+MQyCGK6vrzd/r9frcHt7W3Ru6THCarUKXdc5Mfwb1gBeqemsrTS2TH4LieRCskYWi4W4UUysKJxWRhr7HfPHY9YOP99iMeAeD9JzlKwVPIfawb0Z+D2gNdG3gIxZDBxE1vTOSsZlqg+vAXsnhvV6vSW0Qwjh7Ows+9zSY4jlchmurq6cGAzgmmdKk4y1kco3j4GuSUHNmgmccvlI+zlrC5m05yA9IxJm4/FYDaRqlgz/Hq0bsmZ4m1qVXX7/3AWD17NUE+0DKSsV+2UJdFufraQ0vGTsnRju7+/Dzc3N1ndXV1dhtVplnVt6jLBcLjffOzGkIQm4Ei0+l0zwOiSsUBCXtJdqHwkQSYP+JndNjitMilWQINMEtab54/ecKGPnaO1JeyzzZ2vZf6DkeWvvADX/GCFp7j0NqWdB3/fpwhwi9k4M8/lcFNoPDw9Z55YeCyGEp6enzd8pYvj111/D58+fN58PHz68KGKwatyWAGSNxhaDNmmn02mWUMiBpDFyzdoayMY9FMjPL60ux2taV/KmgtXS7/j3/L448bV4ttiW5FKUngOeXxrX0vrByTzHSnuJGCwxkAZvPbf0WAi/WSL4fYwY3r59+0ybeinEwLXvXODm7n2WANAmJw/MWn9Xek1ODtbUVylWIQkoSVPmz5P6gHs6WKyD2D3yukOthSC3ktBCkdxWfQvj2nH/EsmiN2K4v78Ps9lM/ZCGfn9/L/r/JYshdm7psYeHhy0ieM0Wg0XTj00CzSXSGlrbfMMaLmhzhGQOclNfLy4utgLNFiLV7hldK/TB4HeOhs01877eH+8zt1B4MkJffeEWUukq6L7G1T6xd4tBCwxLwjl2bukxIg76dF0XZrOZGOOQ8JJiDLmTgJ+fU3Gzpj/aRETXC6+0iZpwSWmNWkjaOP+Qzz7HTUKCH0tOYBZVjvbfSgBLLi7p/0jauCitJEstp1983UwsoGy5tlsMO0pXRbfParUK6/XadG7pMUTXdVvXS+ElEYMFKFy1dQg0QVAzLBHGWjqr5nbhRfF4nfx9anZojdG/FBOhbCReSyoFLSiLxMBJaBf3zp+z9blb/Pr0nku2LeXv4OjoyFyWJRbbeYkYBDGs1+swm83Ccrl8tsDs9vY2zOdz07mlx0L4LQA9n89D13Xh7u7uVVoMFuDkSlW0XCzkaql0LDW5UsQgpZXG3CCaW2UXE10LcvJ+5xRe47/lqbR4vZbrPVJIWQyx36XcXkiAOUSHbUuL3nLI69BKaJdgEMRwyHhtxIAT4/T0NOk60tw3mtDn15IECk1iykCijXi09QXSb3FS79OSKBWiuefuEiVaNo6rmHun1GKQlJaS546K0dCee0s4MVTiNRIDd9mUCFQLMWhAf7TkKkGtkgsYi9B6bRvytIYkhFNaNid7bdV1CtztJP1t/e2+rMshwImhEq+NGLj7gmtvqcmouVRK+oCxDiwch9U8S8gLiSaEsv2sXyqscQA8h94B36BH+412vgU4PlNWAkLb2fClu4w0ODFU4rURQyzAzCej1W2TK1i187F+T07bKYshtT5in66ovqAROK+RRMdizyBVEwtJnFZUl1qTMYKK9RED9ly5eA2Ez+HEUInXRgwc3GrQFmjh+fw7FBA1mre0O5sFMYFBQd3RaLS1HeihuxdSPna0yvDZ4DPm1qP23lPJAFIWVakPP2aRxt4T363utdRE0uDEUInXTgwEFCQxN4MEOl5bh4ZcPhcXF0V7LUj9o/uivaIl90JJsLWmTy3AyVBzu3ABiy6Xi4sL8ZwcC4rOHY1GyTLgOfeVmzmEq/ZbjMVDhxNDJZwYfkMswGgVFLnBQg0p108OeLYM5r4TpPuTrCAe0LQQUl+CiQv+HIHMYzAaqaQUgdr3bLmvVLt0Pk+9PkQrsCWcGCrhxPAFWpC2ZJJJgrW2HxZIgkXLUOG/kSwGrYQ3/l+6x1SgNLUnA2rfMSFZQkA8BpMr6PsmPYKlL1ghtzSu8RLhxFAJJ4YvaDnhJcHaQouLtcGtAy4oYr7yWJspi4FvoZnyzVt3cYvdC/aPk17Nc7aMgZL2a5SLWLyCxhd/9q8dTgyVcGL4gj786qjBt1hxqmnp2gbyWilo7GtKGKaIAwO6Mb8/Icdi4IvGYs8ktocDXTe1B0PKuiqF1KdUkDvm3tTacPwGJ4ZKODG0Q0zA8k1xSqEFulGAavs6SGSg7ZVguS++UE/aE6BUyNJvcX2HRkxYhC9mmfHsoZQWXrIOgfdNs7ZCsNVkcsFfBieGSjgxtENsEteU0LCcq22Wo7WrEUxOv5CMckkl1Ta6lGLCXrOgpPb5zm7adqpYDLDEuuPPlq7DlYOUxeAohxNDJZwYdgMU3FKAOIT6GEeq3a57XiW2FRnx41pfOLR75hp27Dxr+5gOLAXRkUhr3ElIVugKQsLBRYyO9nBiqIQTgx3WIK0EFD7oqkCBVLvfQkx4kkAajUZF18A2Tk5Okte3WC3a99bvrP1Gi8NiMeQipvnzvzFtONWOoxxODJVwYrAjptHHjnENVMskaWUxSIKFB3Jz10hwoSpdU/s7dW90LmnyJycnTXztsWBu7cpk/G3Oe4uRf4xMnTTy4MRQCScGO2IuhhyhmKM9t0apVbJYLDYCG6u2WoVi7N446VB7uMVnCWlKv6l9xn20qbVD1zo+Pn71K5lz4cRQCSeGPFg0X752IcdfvU/SiF0PYwaYTtoitZNbDOPxWNTuLVlKsXtogVZuJ8vKZnwGr732US6cGCrhxJAHi+aLG/BQiiQSSY47IfV9X8DrcSGN7qhSAiTw54mL5d68ebOVdssti31ozzXELT3HlMB3F1IZnBgq4cTQDtoGPNyfzTVhSTNPCR4uNFtDCtpq+0dg3KTW3cMrzHKCKtn9rCU0V5Jl8aL0HEs39HHE4cRQCSeGdpAmviTAYmsOrEDiod/2pV1aAq6lC/gkwkttt9oHrNeJPYtUZV7+d20g3KHDiaESTgztENP2ya1Uuo7AYjFYtxttIWytls6QEOujtoGP5b4srkHuMoqRiaMeTgyVcGLoH+gbLy2lnbIsFosvNYtSxLDreEXfsJKSJqRD0DfwsTyr2Dk8ISFlTTrawImhEk4M/QMtBtpIJVcgpIQf10BjWuwhaPc50Pz+/B6lrDECFtiL7eIXazcWQMZ4Qotiio44nBgq4cSwO1iDlKVtX0LuP2VDSRZK3ymx+0qtRYFu0eK5cOd7GlhiKxZwS6V2lzdHGk4MlXBi2B124Vcm4omVvtCEm+Zq4amoVutl1xoxXjeXnNDdR8RgtUQ4DjH+8tLgxFAJJ4bdoQ+NPJXVk9MPzdWCQtPid0/dp9RHSZjm7ORmuW4MnPy0Plna3xcxOr7AiaESTgyHBXRHHR0dZefBl2q8msUgCdTUtaQ9rbkwRTJCK6gk66oVIVsFvlsI+4cTQyWcGIYFq5sGP9w1FWujlTZL18A9DjTXFP++1GKgWk2pktXcZRe7Zyk+kbpnF/jDhxNDJZwYhgXNzx/C9tqFWMojF4Qo0KxF9KwEhWmefQeziYQsFgMv15EiSuuqbXp+pdlljt3AiaESTgzDQiylkrtgrNlFSBRWiyGlZZ+enoajo6OkgKz1++OzsBSSy7EA6Hzc8jN1P7ji3OMIw4UTQyWcGIYJSaBybV8S3qk8e6ugtmrZNQRjvY6UJVXjKuO/xzpPqTIjbjEcBpwYKuHEcHiIacWlNYss17Omr+YQkcX/n3KVWX6j/Z7uhSwgHrOxWliOYcGJoRJODIeHmLBC7bd2vQS5jHh+f6o/Vkum5PvUMet90e+pr1jU0BrIdwwXTgyVcGI4PKSE5mQy2Vr5XCrU0GXEiSFmFUj9s7qUaOU2EVKpVWAFD1Q7XgacGCrhxPAy0ULgkcVANYRyBL3UlsWlxIlIIp4a904qBuN4GXBiqIQTw8tFC4EnxTN4jKFlX6fT6bMAMC8nXnNfHjN4HXBiqIQTg4MQcwFhSW8uXFtr3IvFIozH4036qLTPBAaNcwjKrYPXASeGSjgxOAg8WwctBVxcpqV7phaeSdAENbqUJKuFu57cAnAgBkEM6/U6zOfzsFwuw3w+D09PT0Xnlh4LIYSHh4dwf38fHh4ewsPDg7nvTgwOgpStw60CKXMHdz/LhebaIZfSdDoVzy21GByvA4Mghuvr683f6/U63N7eFp1beuzh4SHc3d1tjl1dXZn77sTgkJCzhqAPi6H2XMfrxt6JYb1ebwntEEI4OzvLPrf0WAghXF1dPbMurHBicOTAs3och4Bcufb/usZ49+5dN5lMtr6bTCbd+/fvu+vra/O5//u//1t07OzsrPv48WN3dnbWvX//vru6uuqurq4a3qHD8QV//vOfuz//+c/J7xyOQ8JR6wY/ffokfv/x48esc0uPvX//vptMJt2PP/7YXV1ddT/88EP3448/qv395z//2f38889bH4fD4XjNaG4xaNCEee65qWMfP37sHh8fu5ubm+7s7Ky7u7vrzs/PuxCC+Jvvv/++++tf/2rum8PhcLx0mInhhx9+6NbrtXr8u+++2whjbh2Qa4cjdm7NMfrQNbquE11ZXdd1f/nLX7r//u//3vz/559/7r799lv1Ph0Oh+Olw0wMd3d3pvNubm66+/v7Z99Pp9Osc6+uroqOSS6rGL766qvuq6++yvqNw+FwvGQ0dyXxQO/j42M3nU63NPezszMxKIzncgsj59h0Ou0+ffrUnZ2ddY+Pj93V1ZVoLTgcDofjOXqJMSyXy+5//ud/uj/+8Y/dP/7xj265XG6Off/9990f//jHbjabJc+tPfaf//mf3Wq16h4eHvq4TYfD4XiRGAUtKvtK8fPPP3e/+93vus+fP3fffPPNvrvjcDgc1ciVa83TVR0Oh8Nx2HBicDgcDscWnBgcDofDsYWdLXA7FFDIxVdAOxyOlwKSZ9aQshMDwy+//NJ1XeeL3BwOx4vDL7/80v3ud79LnudZSQz/+te/up9++qn7+uuvu9FoZP4drZj+8OGDZzMB/Lno8Gcjw5+LjJrnEkLofvnll+73v/99d3SUjiC4xcBwdHTU/eEPfyj+/TfffOODWYA/Fx3+bGT4c5FR+lwslgLBg88Oh8Ph2IITg8PhcDi24MTQCF999VX39u1bL8jH4M9Fhz8bGf5cZOzyuXjw2eFwOBxbcIvB4XA4HFtwYnA4HA7HFjxd1YDHx8fNHtKPj4/d3d2duCOd5dz37993//Vf/9WtVqviawwFrZ5L7Nj79++7ruu66+vr7vHxsfv06dPg9tbYxXPw8XG440PC4GVKcCRxfX29+Xu9Xofb29uic5fLZVitVkF67DnXGApaPZfYsbu7u9B1Xei6Ltzc3ISnp6dGvW+HXTwHHx/ysUMYHxKGLlOcGBJYr9dbDziEEM7OzqrO5S8x5xpDQavnkmrn/v4+PD09DXbC7+I5+Pg43PEh4RBkiscYEnj37l03mUy2vptMJhsTtvTcFr/bJ1o9F0s70nauQ8EunoOPj8MdHxIOQaZ4jCGBT58+id9//Pix6twWv9snWj2XVDufPn3qfvzxx67ruu4f//hH96c//enZXuH7xC6eg48PvZ2hjw8JhyBTnBgKoT342nNb/G6faPVc6BgGzK6urrrvvvuuW6/X5R3cEVo/h9prDAU+PnQMSaa8WmL44YcfogPou+++625ubrqzs7NnLPvx40fRdM05t8Xv+sCun0uqncfHx02WCWVXPD4+DkYr3MVzGNL4sMLHh46DkCnmaMQrhRbEkYJd1nP5Y8+5xlDQ6rnEjq1Wq62A2dPTU+i6blDPZRfPwcfH4Y4PCYcgUzz4nADXPB4fH7vpdLqVR/34+Gg6F4EmXc7vhoJWzyV1bD6fb469e/euu729HdRz2dVziF1jiPDxoeMgZIqJPl451ut1mM1mYblchtlstsW6t7e3YT6fm859eHgIs9ksdF23Ocfyu6Gi1XOJHVutVmE+n4f7+/swm812cFf52MVz8PFxuONDwtBlihfRczgcDscW3JXkcDgcji04MTgcDodjC04MDofD4diCE4PD4XA4tuDE4HA4HI4tODE4HA6HYwtODA6Hw+HYghODw+FwOLbgxOBwOByOLTgxOBwOh2MLTgwOh8Ph2IITg8PhcDi28P8BB+rAGUYer2IAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Modulation\n",
    "\n",
    "FODO = Line('FODO', [QF, DR, SF, DR, BM, DR, SD, DR, QD, QD, DR, SD, DR, BM, DR, SF, DR, QF], propagate=True, dp=0.0, exact=False, output=False, matrix = False)\n",
    "\n",
    "dkf = 1.0E-3*torch.randn(2**10, dtype=torch.float64)\n",
    "dkd = 1.0E-3*torch.randn(2**10, dtype=torch.float64)\n",
    "\n",
    "qx = torch.linspace(0.0, 0.01, 8, dtype=torch.float64)\n",
    "px = torch.zeros_like(qx)\n",
    "qy = torch.zeros_like(qx)\n",
    "py = torch.zeros_like(qx)\n",
    "\n",
    "state = torch.stack([qx, px, qy, py]).T\n",
    "orbit = []\n",
    "\n",
    "data = FODO.data()\n",
    "\n",
    "for i in range(2**10):\n",
    "    data['QF']['kn'] = dkf[i]\n",
    "    data['QD']['kn'] = dkd[i]\n",
    "    state = torch.vmap(lambda state: FODO(state, data=data))(state)\n",
    "    orbit.append(state)\n",
    "\n",
    "qx, px, *_ = torch.stack(orbit).swapaxes(0, -1)\n",
    "\n",
    "plt.figure(figsize=(4, 4))\n",
    "plt.scatter(qx.cpu().numpy(), px.cpu().numpy(), s=1, color='black')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "6de9c74c-7559-403f-ae0f-ed985f82325c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([[-0.4395, 15.4433,  0.0000,  0.0000],\n",
      "        [-0.0522, -0.4395,  0.0000,  0.0000],\n",
      "        [ 0.0000,  0.0000,  0.4963,  5.3596],\n",
      "        [ 0.0000,  0.0000, -0.1406,  0.4963]], dtype=torch.float64)\n",
      "\n",
      "tensor([[-0.4395, 15.4433,  0.0000,  0.0000],\n",
      "        [-0.0522, -0.4395,  0.0000,  0.0000],\n",
      "        [ 0.0000,  0.0000,  0.4963,  5.3596],\n",
      "        [ 0.0000,  0.0000, -0.1406,  0.4963]], dtype=torch.float64)\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApcAAAFeCAYAAADQYzRwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAokUlEQVR4nO3dS29bR7ru8YfaSQTIJxStngXRhPkGukx7EucbKNkRsXEmB5HpNLAnOohg9KA7A8OwcTRMTNtzum3zG5ij3kPK/AZiDwxk1rrQsAA5gNYZrBQvS7UoXta9/j+AiETR8gotFR9W1ftWyfM8TwAAAEAEltK+AAAAABQH4RIAAACRIVwCAAAgMoRLAAAARIZwCQAAgMgQLgEAABAZwiUAAAAiQ7gEAABAZD5J+wIk6erqSr/99ps+//xzlUqltC8HAAAAAZ7n6f379/riiy+0tBQ+P5mJcPnbb79pfX097csAAADADd69e6cvv/wy9OuZCJeff/65JP9iy+VyylcDAACAoH6/r/X19UFuC5OJcGmWwsvlMuESAAAgw27awkhBDwAAACJDuAQAAEBkCJcAAACIDOESAAAAkSFcAgAAIDKESwAAAESGcAkAAIDIEC4BAAAQGcIlAAAAIkO4BAAAQGQIlwCwqEZDunVLKpUm327d8h8LAAVGuASAeW1v+6Hx3j3p4uLmx19c+I8laAIoMMIlAMxidJby6Gj+72OC5vIyIRNAoRAuAWBatVr4LOXKivTkieR59tuTJ/5jgj5+9L9nrRb/9QNAAgiXAHCTWs2fqXzxYvz+Ukna3fXD44cPUr0e/j3qdf8xJmh++un411+8kJaWCJkAco9wCQCT1GrXQ6WZpby6kprN2b9nve7PWAZnMz3P/7sImAByjHAJAGG2t8eDpZmpvGmWclpmNnN3d/z+Fy8o+AGQW4RLALCp1cYLdnZ355+pvEmz6c9ajobMiwvpL38hYALIHcIlAIwy1eCjM5ZbW/GEyqBmczxgXl1JP/5IwASQK4RLADAaDX+2cLQafHdX6nSSu4Zm09+LaXgeM5gAcoVwCQDSMFheXQ3v291NZsYyqF6/PoNJwASQE4RLAAgGy6Ulf/YwjWBpmBnMpT+GaQImgJwgXAJwmy1Y/vJLNNXgi6rX/WshYALIEcIlALf99a/ZDJaGLWD+9a/pXhMATEC4BOCuWk06OfE/zmKwNIIB8+SERusAMuuTtC8AAFLRaIy3G8pqsDTMtd275//XXHua+0IBwIKZSwDuMfssjZWVbAdLI1hFzlGRADKIcAnAPcF9loeH6V7PLIKN1v/xj/SuBQAsCJcA3JKXfZaTNJv+bKvkN1ln9hJAhhAuAbijVsvXPstJRmdbWR4HkCGESwBuCBbw7O7mN1hK9v2X9L8EkAGESwDFFyzgSetYx6gF91/SYB1ABhAuARRb8ASetbViBEuj2fT/nyRO8AGQCYRLAMUWrAx/8CDd64nDgwec4AMgMwiXAIqrCJXh0+AEHwAZwgk9AIopbyfwLIoTfABkRMnzPC/ti+j3+1pdXdX5+bnK5XLalwOgCP70p+Gs5cqK9OFDuteTlNF2S6XScEsAACxo2rzGsjiA4gkuh+fpBJ5F0WAdQMoIlwCKpUiN0udFg3UAKSJcAiiOojVKnxcN1gGkiHAJoBiK2ih9XjRYB5ASwiWA/Ct6o/R50WAdQAoIlwDyLRgsi9oofV7BBusETAAxI1wCyLfgCTwuFvBMEmywzgk+AGJGuASQX42GdHrqf1wqESzDcIIPgARxQg+AfDLL4eYciNu3CZaTcIIPgIQwcwkgf9hnOR9aFAFIAMc/Asif0aMd2Wc5u1u3pIsL/+O1Nenf/073egDkAsc/Aiim4NGOBMvZHR6y/xJAbNhzCSA/ONoxGuy/BBAjZi4B5ANHO0aL/ZcAYkK4BJB9HO0YD46IBBCDmQp6ut2u2u22JKnT6ej58+eqVCqSpF6vp1arpWq1ql6vp729vcHXbkJBD4BQtqMdKUCJFgVSAKYwbV6bac9lu93WTz/9JEl6/Pixvv76a719+1aS9O233w4+7vV6+uGHH/T69et5rx8AfPv7tByK24MHwwBvjoiUCJgA5jL1sni329XDhw8Hn+/s7Kjb7arX66nX6409tlqtDmY4AWButdqwZY7EjFpcbEdEskQOYE5Th8uNjQ09f/588PnZ2ZkkaW1tTe12W2tra2OPX1tbU7fbjeYqAbiHAp5kcQY5gIjMVNCzs7Mz+Pjly5e6c+eOKpXKIGgGnZg9PAGXl5fq9/tjNwAYoIAnHSZglkr+56enzF4CmNlc1eJnZ2dqtVo37qkMC50PHz7U6urq4La+vj7PZQAoIlsBD8EyOfW6f0675J/bzvI4gBnNFS4PDg705s2bQTV4pVK5Nkt5cnISWi1+//59nZ+fD27v3r2b5zIAFA1nhmfDgwfsvwQwt5nPFn/8+LF2dnZUrVYHM5MnJydj1eKSdPv2bf3rX/+aqh0RrYgASKIlTpbQAgpAQCxni7daLW1sbAyC5atXr1SpVFStVsce1+v1tLW1NXWfSwDgzPCMCRb4cAY5gClNPXPZ6/X01Vdfjd1XqVR0eno6+PrTp0+1vb2tTqej+/fv00QdwHSCZ4Y/eUKwzIpGY3gGuURxFeCwafPazMvicSBcAg4LBkvCS/bwbwRAMS2LA0CkbL0sCS3ZEzyD/MULlsgBhCJcAkjPaJPulRWCZZbZAiYV5AAsCJcA0hEs4Dk8TPd6cLNgwNzfT+9aAGQW4RJA8oJ7+KgMz49mc1hBfnHB8jiAawiXAJJlKw4hWObLf/7n8GP2XwIIIFwCSA5Vx8VAgQ+ACQiXAJJBZXixUOADIAThEkD8zFGCBsGyGIIBkzPIAYhwCSButjOqCZbF0Wz6/6aS/29MwAScR7gEEK/9/WGwXFqSHjxI93oQvQcPhhXkBEzAeYRLAPGp1fx2NQYth4qpXvf/bUcDJj0wAWcRLgHEg5ZDbjEB06AHJuAswiWA6NFyyE31Oi2KABAuAUSMlkNuo0UR4DzCJYDo0HIIEi2KAMcRLgFEg5ZDGEWLIsBZhEsAiwsGS1oOQaJFEeAowiWAxQV7WdJyCJK9RREBEyg8wiWAxdDLEpMQMAHnEC4BzI9elpgGTdYBpxAuAcyHXpaYBU3WAWcQLgHMjl6WmIetyTrL40DhEC4BzIZellgEPTCBwiNcApgevSwRBXpgAoVGuAQwHXpZIkr0wAQKi3AJ4Ga2YEnLISyCFkVAYREuAdyMJumIAwETKCTCJYDJaJKOOBEwgcIhXAIIR5N0JIEm60ChEC4B2NEkHUmiyTpQGE6Gy1pNKpX8261brL4A1xAskQZbk3UCJjCm0fCzS6nkT/Zn8VfEyXD56tXw44sL6d49wiYwQLBEmoJN1gmYcNxomCyV/MxitsF73nimyQonw+V334V/bTRsEjThHIIlsoCACceNBsrRMBlUKk3ONGlxMlw2m37af/JE+vTT8MeZoEnIhBMIlsgSAiYcVatNDpSSHyp3d/3atywO006GS6Nelz5+9IOmCZsrK9cfd3Eh/fgjARMFRrBEFhEw4RAzWzk6FBsrK35GMXklq6HScDpcBtXr0ocP9qDpecxioqAaDYIlsouAiYIzoTI4W2lmJz3PzyZ56gJHuAxhgubomCb5//D090VhmGMdjZUVgiWyh4CJgjJDcHAJPMtL3tMgXN6g2bw+i3l15b/DYGxDrtnOCz88TPeagDC2gMm7fOSY2VtphmBpuPyd11BpEC6nYGYxnzwZHiAh8eYZOWYLlhzriKxrNsff6bOMhJwKbnNfWvIzRt6Wv8MQLmdgDpAYHdtevGAfJnKGYIk8OzzkHHLkWjBYrqwUbwgueZ7npX0R/X5fq6urOj8/V7lcTvtypmL74fjwIb3rAaZCsEQR8HOMnMp7dpg2rzFzOafg9h+OwUXm8YKMojDLSMxgIkdswbKo29w/SfsC8sxsuDU/LOa/ed+IiwIiWKJozM+u+bk2AXP0a0BGuNZKmJnLBdEhA7mwv0+wRPEwg4kccC1YSoTLSBAwkWm12ngTNYIlisQWMPf3070m4A8uBkuJcBkZAiYyyTayESxRNCZgGmyCRwa4GiwlwmWk6PGLTHF5ZIN76nXe4SMzXD9Vl3AZsWDAZPsPUkGwhItYQkJGjO7McHH4JVzGoNmU1tb8j9lfjsQRLOEyAiZSNrrNvVRyc/glXMbkwQP2lyMFBEuAgInUBIfg779P71rSRLiMCfvLkTiCJTBEwETCGIKHCJcxYn85EuP67nHAhoCJhBAsxxEuY8bYhtiZ03eMlRW3RzVgFIMwYkawvI5wmQDGNsTGdqxjUQ+rBeZFnzjEhGBpR7hMCGMbIsd54cD0mk1/Vt+gjQcWxG6kcITLBNEDE5EhWAKzOzzkHHJEIrgbiWA5bqZw2e12tbm5ab2/2+1Kknq93uBjXEcPTCyMYAnMx3YOOYMw5rC/PxyC19YIlkFTh8tWqyVJ1uD49OlTbW5uqlQq6e7du6pWq9FdYQHRAxNzq9Wke/cIlsC8bAHz3j02wmNqo03SJf81HeM+mfaBOzs7oV/b3NzU6empJKlSqSx8UUVncsC9e/5/TQ9M3vlgouDOcYIlMB/zOzO6AmB+txiIMYGtgIch+LrI9lxWKhWC5QzogYmZBHeOr6wQLIFFmBnM0SIfBmJMQGX49KaeuZzk7OxssGze6XRuXBq/vLzU5eXl4PN+vx/FZeSO+aE0P6y8cYaVrY/lhw/pXQ9QFPW6fxtNDQzEsCBYziaScLm3tzeYtaxWq/rmm290fHwc+viHDx/q559/juKvzj0CJiaijyUQP9tA/Oc/szIASQTLeZQ8z/Nm+gOlkoJ/pNvtamNjQ5I/i3n79m0dHx+Hzl7aZi7X19d1fn6ucrk86/9DIfDDi2uoCgeSdevWsFKD3zeI1+agfr+v1dXVG/Pawnsuu92uvv7662v3r5l+OxbLy8sql8tjN9dxig/GECyB5NEHEyMIlvObK1yenZ0NPq5Wq3r06NHg83a7rZ2dHYp75kDAhCSCJZAW+mDiD5y+s5ip91y22229efNGkr9ncnt7exAit7a29PjxY1UqFR0fH+v169exXXDRsQfTcbQbAtIVbFNk+mD+858MxI6w1VDyTz+bmfdcxmHaNXyXjGaMUmk4iYUCI1gC2RFcQZCYvnLEn/4knZz4HzMMj0tszyXi0WwO2695HsvjhUcfSyBb6IPppFqNYBkFwmWGjXacYUwrsLA+loxoQLrqdf93kc3wTgguHhEs50e4zDBO8XEAfSyB7LNVW1LkUygc6xgtwmXGUUFeYLWaXyhAVTiQfaN7lST/d5fBuBBoORQ9wmUOEDALiOIdIH9G+2BKDMYFQLCMB+EyJwiYBRIczSjeAfKBIp9CIVjGh3CZIwTMArCNZhTvAPlBkU8h0CQ9XoTLnCFg5hhvk4HiYDDOLZqkx49wmUMULuYQwRIoHttgvL2d3vXgRjToSAbhMqeChYscf5thtj2WBEugGIIB8+iIGcyMsgVLtrvHg3CZY6OFi1dXBMxMsgVL3iYDxdJsSltbw89ZIs8cgmWyCJc5ZgoXRwMmrdcyotGQbt2ieAdwRadzfYn81i3e8WcALYWT90naF4DFmF+O0XdkJs+w8pqS4FtkiT2WgAvM77gZhC8uhpUjJJlU0FI4HcxcFgCt1zKEYAm4LbgHkz1LqaGlcHoIlwVB67WM2N8fX3t58oRgCbim2fR/99kUnxpbH0t2JSWHcFkwtF5LidljeXExvI+3yIC72BSfGvpYpo9wWUD0wUyYGclGg+XuLsEScF0wYEoMyDGjj2U2EC4Lij6YCWGPJYBJTMAslYb3MSDHgnZD2UG4LDD6YMbM1t+CPZYAgup16ddfWSKPEcEyWwiXBWbb8kPAjAj9LQDMImyJnIC5MIJl9hAuC46AGQP6WwCYB33jIkewzCbCpQMImBHa3qa/BYD5hfWN4zSfmREss4tw6QgCZgRqNenoaPg5hTsA5hVs62FO82FQnlqwrTDBMjsIlw6xBcz9/XSvKTeCS+FbWwRLAIuxnebDoDyVWo22wllGuHSMCZjGxQWrMROZ5ujBpfBOJ71rAlAc5jQfg0F5orAhmWCZLYRLB9XrrMZMJaw5OjOWAKLEoDwVhuT8IFw6yrYaw1g2ItjDUmIUAxAf26BML8wBzqvIF8Klw8xqDEU+AbYeljRHBxC34KAs0apI9qpwhuRsI1w6jiryAHpYAkgTvTDH0G4onwiXIGAawWBJD0sAaaAX5gDthvKJcAlJBEzrjCVrLgDSZOuF6VCrItoN5RfhEgNO9sG09bVYWZEOD9O7JgAwbAHTgRlM20ISwTI/CJcY41QfzLC+FiyFA8gSh07zCetjyUJSvhAucY0TLddoNQQgTxxoVUQfy+IgXMKq0H0waTUEII8K3KqIPpbFQrhEqEL2wdzeptUQgPwKa1WU4/1L9LEsHsIlJipUFXmtJh0dDT9nfyWAPLK1KsppJTl9LIuJcIkb5T5g2naIb23xthhAvhWgkpw+lsVEuMRUbAEzF3vJw3aIdzrpXRMARCWnleTmPT99LIuJcImpBQOm5E8GZnYMoyIcgAtyVkke9p6fYFkchEvMxATMUml4XybfJFMRDsAlOakkpyrcDYRLzKxel379NcN7MKkIB+CijFeSUxXuDsIl5pLJIh+ziYeKcACuCqskT3mZnKpwtxAuMbdMBUzbJh4qwgG4KrgPU0ptkzzB0j2ESywkE1XkYYU7VIQDcJnZhzm6ST7hATo4PBMs3UC4xMLCqsgTGb8o3AGAcMFN8lJiA3SjcX14Jli6gXCJSITtI49t/LI1RqdwBwCuS6HQxyyFG6USw7NLCJeIjG0feSxbfMw6S7BJGoU7AGCXYKGPbY/lr78yPLuEcInINZvjb5AjHbuC6ywSTdIAYFphhT4RDdLssYREuERMDg9j2OJjRi2jVGJ/JQDMyhT6BJfJt7cX+ra2LfAESzcRLhGLSPdg2vZXss4CAPOzLZMfHc29DzMYLNkC7zbCJWITtgdzpoBp21/JqAUA0Wg2/Z7Axhz7MIPBki3wIFwidsEtPlMXKQZHLIlRCwCi1unMvQ8zeNouW+AhES6RkGDAvLiQ9vdDHhzWZoj9lQAQj7B9mBNmAmq166ftMkRDIlwiQbaAee2NcaMh/fgjbYYAIGkztCsKLixx2i5GES6RKNsS+WDMMvsrPW/4AN4KA0CywtoV/VFNbttjyWm7GDVTuOx2u9rc3Lx2f6/X0+PHj9VqtfT48WOdnZ1FdX0oIFvAfPkfNXn0rwSAbLAsk3tHR3r5HzX2WOJGJc8bnSYK12q1VK1Wtbm5qeAf2dzc1Nu3byX5QfPg4ECvX7+e+iL6/b5WV1d1fn6ucrk8w+Ujz2o1qfyiof+nfd3ShUqSPEmllRW/USbL4ACQvu1teUdHgzH6g1b0f3Wo/m6dYOmYafPaJ9N+w52dHev9vV5v7PNqtap2uz3tt4XDmqrJ0wuV/vjck/QP7er8sEmuBICMaPyfjlaPavr+j/H6f+lCT3RPJf1TEukS1y2857LdbmttbW3svrW1NXW73UW/NYpqpBp8NFg2tauamnEcdQsAmIOpsaypqaZ2ZdYtS9IMfeXgmoXDZdj+ypOTk9A/c3l5qX6/P3aDI2zV4Csr+p/dJ/rfS8N3wBEedQsAmEOjIf3lL8May/9SU/+zG2hXNEfTdRRfbNXik4p6Hj58qNXV1cFtfX09rstAlgRHKmnQZujPzbr1uEjeEANA8sxwfXU1vG93V/pz09KuSGLAxpiFw2WlUrk2S3lycqJKpRL6Z+7fv6/z8/PB7d27d4teBrLOtBkKjlQju8FNi7XRgMkbYgBIVnC4XlqynGFhqslLpeF9DNj4w8Lh8s6dO9b7t0bPKg1YXl5WuVweu6GgbKftWEeqocND/yEGS+QAkIxgD8ulJemXX0Kad9Tr0q+/Xh+w2YfpvLnC5eiSd7VaHftar9fT1tbWxJlLOMK8/Q3srwwfqXz1uqxL5ARMAIhPMFhOMVzbB2z2YTpv6nDZbrd1cHAgyd8z2Wq1Bl97/fq1Dg4O1Gq19PTp05l6XKKgtrfHRylppmMcbaeQ8YYYAKJnW2Ca6dRd24AtMWg7bOom6nGiiXqBNBrS/v712coFmqLPtEwDAJjJrVvjQ/ZCp+7YXgMW/qbIimnzGmeLIzq2ZfCtrRne/toFj4u8uvKrGHkzDACLqdUiDJbS5FlMlsmdQbjE4mxrKpI/uHQ6kfwVpjDR7Bu/uvLbZRIwAWB2YUvhkU0u2qrJWSZ3BuESiwkr2plQDT4vs2/c8DxmMAFgVqaHZaQzljammpxiH+cQLjGfSbOVCy6DT1Kvs0QOAPMKa44e23ZIin2cRLjE7EKOcIxjttLGtkTOG2EAmGyq5uhxMQM3s5hOIFxiNmZ0shzhmGT5tlkip9k6ANwsE103Js1ibm8neCGIG+ES05m0DJ5Se4mwZuustADAULDt8FTN0eNkm8U8OmLwLhDCJW6WYNHOrGxvhC8u2IcJAJI/fB8dDT9PYaHJzgzeo0dFs0xeGIRLhKvV/DYSCRftzMPWC5NWRQBcZVts2tpKfT7guk7Hvky+tETIzDHCJexsxzdmZLYyjFlpMWhVBMBFYa2GImo7HD3bMrnnsc8pxwiXGGfe7o6uo5RKmZyttLG1KmKVBYArghXhUk5OXgwr9mGpPJcIlxgKO77x6ioHI9NQsFWRRCU5gOKzVYRneLHJrtn0Zy05PjLXCJdI5PjGpFFJDsAVtiE89YrwRXF8ZK4RLl1mRqSMVoIvikpyAC7Y37++vzIHu5huNun4SEJmphEuXWVbApcKNCoNUUkOoIjM/EDsZ4Snib2YuUS4dFFwY45UmNnKMFSSAyiSsIrwgg7h9opyidN9Mopw6RJb38pSyf+FLdhspQ2V5ACKILcV4Ysys5i2033oi5kphEsXhBXsrKz4+1kKHipHUUkOIM8KURG+KNvpPvTFzBTCZdHZ9lbmqG9lHKgkB5A3hawIX5TtdB/2YmYC4bKoGg3ps8/s7YVy1rcyDlSSA8iTwlaEL2pSX0yWylNDuCwiM1v5++/D+wpesDMvWyU5b3oBZIUTFeFRmHSEJAU/iSNcFomtYEfy96XwFjdU2D5MZjABpMm5ivBFhbUtouAncYTLIphUsPPkSW5P2UmS2Yc5ehgEM5gA0uJsRXgUzFI5BT+pIVzmWdgJO44X7MzLHAYRnMFkLAKQFNtcgZMV4VGYVPDDwB4rwmVeNRr+MTO2E3Yo2JmbrZKcQh8ASbAtgztfEb6osIIfqspjRbjMm9HZSs8b3k/BTmRs23Yo9AEQp7BlcBagIjLphJ/lZWYPIka4zIuwJXCJESgmYYU+rKYAiBKN0RMSdsLPx48slUes5Hmj01/p6Pf7Wl1d1fn5ucrlctqXkz22s8Al/5fj8JBQGbNG43qPuZUVf4wCgEUEh/dSybmD09Kzve1XkgdRORVq2rzGzGWWhbUWMkvgzFYmIqzhOm9yAcwr7MQdgmWCOp3wpXJaFy2EcJlFN7UWIlSmIthwnf3gAOZhO5WX3U0pCVsqN62LCJlzIVxmCfsqMy8YMCX2YQKYTti8AauwGRDWgJ3+mHMhXGaF7a2sRBV4BtmKDmlXBGCSsDZDDO8ZQ+uiSBAu08a+ylyiXRGAaZm2xLQZypFJrYtYKr8R4TItjYb02WeEypwLa1e0vZ3eNQHIDrMoNdqXhdnKnGA/5twIl0kb3Vf5++/D+zmyMbdsp/ocHbFFB3DZpP2VDPE5w37MmREukzKpWGdriyMbc86MPVtbw/vYhwm4if2VBXXTfkxC5gDhMm6TQqUZbTqddK4Nket02IcJuIxjHB0Qth+TkDlAuIzTTRXgjDaFxLGRgHtsy+Ac41hgYfsxJUKmCJfRMyMMFeBOs+3DpJMFUEymGjy4DP7LLwz1hUfItCJcRmXS8jfFOk4K2wP+4gUBEygKWzU4w72DCJljCJeLMn0qJ4VKinWcZpbJS6XhfSyTA/nGaTuwCptVkIYh04EWRoTLeUxa+paGy9+ESvyhXpd+/ZVlcqAIbNvpqQbHGFNZbpvJNC2MSiVpebmQswyEy2mNBkrbLKXEnkpMNGmZnFlMIPsmzVYy7MNq0nK5JH386GeKUqlQLwSEyzCjYXJSoJQIlZhJ2NnkzGIC2RVWtMNsJaZiQqatT6ZhXghM7sjx8rmb4dLsk5x0mxQmpeGo4nmESsxs0iwmR0cC2ULRDiI1acl81OjyedgtowHUzXD56tV8f45AiYjZZjE5OhLIBop2EKvR2cybgmYYz5s/08TIzXD53XfTPW40TBIoEZOwoyNZJgfSQ9EOEjUaNGcJm6XS9JkmQSXPG53oT0e/39fq6qrOz89VLpfTvhwgNbWavff+4SHva4AkNBrS/v71XVHMVgLT5zU3Zy6BjKLYB0gPs5VANAiXQMbQsghInm3VgKIdYD6ESyCjmMUE4meah4wGy1KJ2UpgEYRLIMOYxQTiEVYJvrLin6bFbCUwP8IlkAPN5vWAySwmMB/b3spSiWVwICqESyAnTN9dZjGB+UzqW3l1xTI4EBXCJZAzk/ZiEjIBOyrBgeQQLoEcCtuLeXHhn39MwAR8k2YrWQIH4hFZuOx2u+p2u5KkXq83+BhAfGyzmJ7HXkxAmry3ktlKID6RhcunT59qc3NTpVJJd+/eVbVajepbA5hgUkX50hIhE+6xtReS2FsJJCWycLm5uanT01Odnp7qzZs3qlQqUX1rAFMIm8Wk4AeumNReiL2VQHIi3XNZqVQIlUCKJu3FvHdP2t5O57qAuNFeCMiOyMLl2dmZWq2WWq2WDg4O1Ov1Qh97eXmpfr8/dgMQnbC2RUdHzGKiWGgvBGRPyfM8L4pvdHZ2Npi17Ha7+vbbb3V8fGx97N///nf9/PPP1+4/Pz9XuVyO4nIA/KHRkP77v6Xffx+/f2VFOjxkRgf51GhI+/vjM5WSP1v5/feESiAO/X5fq6urN+a1yMJlt9vVxsaGJD9o3r59W8fHx9bCnsvLS11eXo5d7Pr6OuESiFGtdn12R6JyFvnDzzKQjmnDZSTL4t1uV19//fW1+9fW1qyPX15eVrlcHrsBiJet4Eeiqhz50WhIn31GwQ6QdZGEy2q1qkePHg0+b7fb2tnZobgHyBhT8BNWVb68zH5MZI/ZV3nv3vj2DhMqKdgBsiXSZfF2u61KpaLj4+OxsHmTaadZAUSL5UVkXdjP6NaW1Okkfz2AyxLfc7kIwiWQru1tv5J8FIURSFNYqKQQDUhPonsuAeRbpxO+VM5+TCQp7HQdlsCB/CBcApAU3oCdkIkkhPWrpBE6kD+ESwBjwhqwc5Qk4jBarBPsWUkjdCCfCJcArMJCpjlKkpCJRUwKlbQWAvKNcAlgorD+mIRMzGOaUMkSOJBvhEsANwrrjykNQyZ7MjEJoRJwB+ESwNQmhUz2ZMKGUAm4h3AJYGbTzGQSMt02KVRSAQ4UG+ESwNxYLkeQ6VM5aaaSCnCg2AiXABY2zXI555YXW1ifSonlb8A1hEsAkZkUMj9+ZCaziKaZqSRUAm4hXAKInAmZnidtbY1/zcxklkrsy8wrM0tpO6ZRIlQCriNcAoiV7dxyg32Z+TKpSEciVALwES4BxG50JjN44o/E+eVZN2np21R+ex6hEoCPcAkgUeZYyUnFPyyZp++mpW8TKqn8BhBEuASQiknFPxJL5mlpNKTPPrt56ZtQCSAM4RJAqqZdMmc2Mz6js5T37km//z7+dZa+AcyCcAkgMyYtmUvD2UyC5uKCgdI2S/nZZ8xSApgd4RJA5kxqZWQQNOczqTjHMEvfl5fMUgKYHeESQKZ1OpNnMyX2Z97EBMqw4hxpGChZ+gawKMIlgFwYnc0MC5qj+zNdntEcXfKeFCjNsjeBEkCUCJcAcuemIiDDpaXzafZQSuPFOSx7A4gD4RJArpkioFmCZlGWz0eXu6cNlBTnAIgb4RJAYYwGzSdPpE8/tT8uuHxeKuUjbI6GyUnL3RKBEkB6CJcACqlelz5+vLkYyAiGzbRnN4NB8qYwKY0X5RAoAaSFcAmg8Eb3aN60fG7YZjejDp62ADltkDS2tob/XxTlAMgCwiUA54wun08bNkdNCp6z3KYNkMboUre5dTqzfQ8AiBvhEoDzgmFznsAZh9Flbpa6AeQF4RIALGyBM47gObqsHbyxzA0gjwiXADCjScFz1hvL2gCKhnAJAACAyBAuAQAAEBnCJQAAACJDuAQAAEBkCJcAAACIDOESAAAAkSFcAgAAIDKESwAAAESGcAkAAIDIEC4BAAAQmU/SvgBJ8jxPktTv91O+EgAAANiYnGZyW5hMhMv3799LktbX11O+EgAAAEzy/v17ra6uhn695N0UPxNwdXWl3377TZ9//rlKpVLsf1+/39f6+rrevXuncrkc+98HH897Onje08Hznjye83TwvKcjjefd8zy9f/9eX3zxhZaWwndWZmLmcmlpSV9++WXif2+5XOYXIQU87+ngeU8Hz3vyeM7TwfOejqSf90kzlgYFPQAAAIgM4RIAAACRcTJcLi8v629/+5uWl5fTvhSn8Lyng+c9HTzvyeM5TwfPezqy/LxnoqAHAAAAxeDkzCUAAADiQbgEAABAZDLRiigpvV5PrVZL1WpVvV5Pe3t7qlQqaV9W4XW7XbXbbUlSp9PR8+fPed4TdnBwoPv37/O8J6TdbqvX66larUqS7ty5k/IVFV+v11O73dba2pp6vZ52dnYGzz+i1e129cMPP+jt27dj9/MaG5+w5zyzr6+eQzY2NgYfHx8fezs7OylejTsePXo09vHovwPi9/btW0+Sd3p6mvalOOHNmzfe3t6e53n+OFOtVlO+IjeMjjOe5w3+DRCt169fD8aUIF5j4zHpOc/q66szy+K9Xm/s82q1Okj7iE+329XDhw8Hn+/s7Kjb7V7790B8RmfQEL+7d+/q0aNHkvxx5s2bNylfkRtevnyZ9iU4YWdnRxsbG9fu5zU2PmHPeZZfX50Jl2a5ZNTa2pq63W5KV+SGjY0NPX/+fPD52dmZJF37t0A8Wq2WdnZ20r4MZ/R6PZ2cnKhSqajb7ers7Ixgn5C1tTVtbm4Olse/+eabtC/JKbzGJi/Lr6/OhEvzpAednJwkeyEOGg03L1++1J07d7KxJ6Tgzs7OeJ4T1u12tba2Nth39uzZM7VarbQvywmvX7+WJH311Vd6/fo1b6oSxmtsOrL6+upUQY9N2C8Eond2dqZWq3VtQzLi8erVK+3t7aV9GU45OTlRr9cbDPB7e3u6ffu2PNoJx67dbuvRo0fq9Xq6e/euJOnp06cpXxV4jU1G1l5fnZm5rFQq195BmeUrJOPg4EBv3rzhOU9Au93Wd999l/ZlOKdarapSqQx+xs1/WRqMV6/XU6fT0Z07d7S3t6fj42O9evUqE3vPXMFrbLqy9vrqTLgMawWytbWV8JW46fHjxzo4OFC1WtXZ2RnvZhPw6tUrPXv2TM+ePVOv19PDhw8JOTFjf2U6ut2utre3B59Xq1Xdv3+fcSZBvMamJ4uvr86Ey+Cg3+v1tLW1lZmUX2StVksbGxuDH/xXr17xvMfMzOCYm+RXMdsqDhGdarWqra2tweBuKvV53uO1sbGhTqczdt+///1vnveYjYYYXmOTEQyOWX19deps8V6vp6dPn2p7e1udToem0gno9Xr66quvxu6rVCo6PT1N6YrccnZ2pmfPnung4EB7e3sEzAScnZ3p4OBAm5ubevv27WBGAfFqt9vqdruDMf3OnTs87zFot9t68+aNHj9+rJ9++knb29uDohJeY+MR9pxn+fXVqXAJAACAeDmzLA4AAID4ES4BAAAQGcIlAAAAIkO4BAAAQGQIlwAAAIgM4RIAAACRIVwCAAAgMoRLAAAARIZwCQAAgMgQLgEAABAZwiUAAAAiQ7gEAABAZP4/xjqpL/pxoFwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Matrix output and twiss parameters\n",
    "\n",
    "from twiss import twiss\n",
    "from twiss import propagate\n",
    "from twiss import wolski_to_cs\n",
    "\n",
    "QF = Quadrupole('QF', 0.5, +0.25)\n",
    "QD = Quadrupole('QD', 0.5, -0.20)\n",
    "SF = Sextupole('SF', 0.25)\n",
    "SD = Sextupole('SD', 0.25)\n",
    "DR = Drift('DR', 0.25)\n",
    "BM = Dipole('BM', 3.50, torch.pi/8.0)\n",
    "\n",
    "FODO = Line('FODO', [QF, DR, SF, DR, BM, DR, SD, DR, QD, QD, DR, SD, DR, BM, DR, SF, DR, QF], propagate=True, matrix=True)\n",
    "FODO.ns = 0.01\n",
    "\n",
    "state = torch.tensor([0.0, 0.0, 0.0, 0.0], dtype=torch.float64)\n",
    "\n",
    "print(torch.func.jacrev(FODO)(state))\n",
    "print()\n",
    "\n",
    "out, *ms = FODO.container_matrix\n",
    "for m in ms:\n",
    "   out = m @ out\n",
    "print(out)\n",
    "print()\n",
    "\n",
    "*_, w = twiss(out)\n",
    "ws = [w]\n",
    "\n",
    "for m in FODO.container_matrix:\n",
    "    w = propagate(w, m)\n",
    "    ws.append(w)\n",
    "\n",
    "ws = torch.stack(ws)\n",
    "\n",
    "_, bx, _, by = torch.vmap(wolski_to_cs)(ws).T\n",
    "\n",
    "s = torch.linspace(0.0, 12.0, len(bx), dtype=torch.float64)\n",
    "\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.scatter(s.cpu().numpy(), bx.cpu().numpy(), s=1, color='blue')\n",
    "plt.scatter(s.cpu().numpy(), by.cpu().numpy(), s=1, color='red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "49f36929-f210-4150-9030-678f82080a70",
   "metadata": {},
   "outputs": [],
   "source": [
    "# More line properties\n",
    "\n",
    "# Define simple FODO based lattice using nested lines\n",
    "\n",
    "DR = Drift('DR', 0.25)\n",
    "BM = Dipole('BM', 3.50, torch.pi/4.0)\n",
    "\n",
    "QF_A = Quadrupole('QF_A', 0.5, +0.20)\n",
    "QD_A = Quadrupole('QD_A', 0.5, -0.19)\n",
    "QF_B = Quadrupole('QF_B', 0.5, +0.20)\n",
    "QD_B = Quadrupole('QD_B', 0.5, -0.19)\n",
    "QF_C = Quadrupole('QF_C', 0.5, +0.20)\n",
    "QD_C = Quadrupole('QD_C', 0.5, -0.19)\n",
    "QF_D = Quadrupole('QF_D', 0.5, +0.20)\n",
    "QD_D = Quadrupole('QD_D', 0.5, -0.19)\n",
    "\n",
    "SF_A = Sextupole('SF_A', 0.25, 0.00)\n",
    "SD_A = Sextupole('SD_A', 0.25, 0.00)\n",
    "SF_B = Sextupole('SF_B', 0.25, 0.00)\n",
    "SD_B = Sextupole('SD_B', 0.25, 0.00)\n",
    "SF_C = Sextupole('SF_C', 0.25, 0.00)\n",
    "SD_C = Sextupole('SD_C', 0.25, 0.00)\n",
    "SF_D = Sextupole('SF_D', 0.25, 0.00)\n",
    "SD_D = Sextupole('SD_D', 0.25, 0.00)\n",
    "\n",
    "FODO_A = Line('FODO_A', [QF_A, DR, SF_A, DR, BM, DR, SD_A, DR, QD_A, QD_A, DR, SD_A, DR, BM, DR, SF_A, DR, QF_A], propagate=True, dp=0.0, exact=False, output=False, matrix=False)\n",
    "FODO_B = Line('FODO_B', [QF_B, DR, SF_B, DR, BM, DR, SD_B, DR, QD_B, QD_B, DR, SD_B, DR, BM, DR, SF_B, DR, QF_B], propagate=True, dp=0.0, exact=False, output=False, matrix=False)\n",
    "FODO_C = Line('FODO_C', [QF_C, DR, SF_C, DR, BM, DR, SD_C, DR, QD_C, QD_C, DR, SD_C, DR, BM, DR, SF_C, DR, QF_C], propagate=True, dp=0.0, exact=False, output=False, matrix=False)\n",
    "FODO_D = Line('FODO_D', [QF_D, DR, SF_D, DR, BM, DR, SD_D, DR, QD_D, QD_D, DR, SD_D, DR, BM, DR, SF_D, DR, QF_D], propagate=True, dp=0.0, exact=False, output=False, matrix=False)\n",
    "\n",
    "RING = Line('RING', [FODO_A, FODO_B, FODO_C, FODO_D], propagate=True, dp=0.0, exact=False, output=False, matrix=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "cd9be56f-a79a-4bce-927b-b0b5e48b8db3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FODO_A\n",
      "['QF_A', 'DR', 'SF_A', 'DR', 'BM', 'DR', 'SD_A', 'DR', 'QD_A', 'QD_A', 'DR', 'SD_A', 'DR', 'BM', 'DR', 'SF_A', 'DR', 'QF_A']\n",
      "\n",
      "RING\n",
      "['FODO_A', 'FODO_B', 'FODO_C', 'FODO_D']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Name and names\n",
    "\n",
    "print(FODO_A.name)\n",
    "print(FODO_A.names)\n",
    "print()\n",
    "\n",
    "print(RING.name)\n",
    "print(RING.names)\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "c338c46d-e07f-4916-891b-6f2bbbefb79e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "18\n",
      "72\n",
      "2\n",
      "8\n"
     ]
    }
   ],
   "source": [
    "# Scan (recurcivly traverse lattice and yeild elements with matching attribute)\n",
    "\n",
    "# All elements, since all elements have a name\n",
    "\n",
    "print(len([*FODO_A.scan('name')]))\n",
    "print(len([*RING.scan('name')]))\n",
    "\n",
    "# All dipoles\n",
    "\n",
    "print(len([*FODO_A.scan('angle')]))\n",
    "print(len([*RING.scan('angle')]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "84806bea-4732-47aa-92aa-8d22c1dcec48",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4\n",
      "72\n"
     ]
    }
   ],
   "source": [
    "# Flatten\n",
    "\n",
    "print(len(RING))\n",
    "\n",
    "RING.flatten()\n",
    "print(len(RING))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "06e0e482-d538-41aa-a226-deb2a78e1d02",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4\n",
      "['QF_A', 'DR', 'SF_A', 'DR', 'BM', 'DR', 'SD_A', 'DR', 'QD_A', 'QD_A', 'DR', 'SD_A', 'DR', 'BM', 'DR', 'SF_A', 'DR', 'QF_A']\n"
     ]
    }
   ],
   "source": [
    "# First positon of an element\n",
    "\n",
    "print(FODO_A.position('BM'))\n",
    "print(FODO_A.names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "38e602f1-09b8-42dc-864c-28c9281a0407",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "QF_A\n",
      "['QF_A', 'DR', 'SF_A', 'DR', 'BM', 'DR', 'SD_A', 'DR', 'QD_A', 'QD_A', 'DR', 'SD_A', 'DR', 'BM', 'DR', 'SF_A', 'DR', 'QF_A']\n"
     ]
    }
   ],
   "source": [
    "# Get first element\n",
    "\n",
    "print(FODO_A.start)\n",
    "print(FODO_A.names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "72c98e77-88a3-4985-8b9d-a0657cca8808",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['BM', 'DR', 'SD_A', 'DR', 'QD_A', 'QD_A', 'DR', 'SD_A', 'DR', 'BM', 'DR', 'SF_A', 'DR', 'QF_A', 'QF_A', 'DR', 'SF_A', 'DR']\n"
     ]
    }
   ],
   "source": [
    "# Set first element (rotate sequence)\n",
    "# Note, first mathced occuranve is used\n",
    "\n",
    "FODO_A.start = 'BM'\n",
    "print(FODO_A.names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "9f284c25-3c10-4fd6-b9ef-8d0914c8c3ce",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['SF_A', 'SD_A', 'SF_B', 'SD_B', 'SF_C', 'SD_C', 'SF_D', 'SD_D']\n"
     ]
    }
   ],
   "source": [
    "# Itemize (list of all elements with matching kind)\n",
    "\n",
    "print(RING.itemize('Sextupole'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "801aa567-b8d4-4eb3-9055-acfb4a567e05",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "72\n"
     ]
    }
   ],
   "source": [
    "# Number of first level elements/lines\n",
    "\n",
    "print(len(RING))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "0363c0e0-fc6c-4adb-b4de-9b292653ac3f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Drift(name=\"DR\", length=0.25, dp=0.0, exact=False, ns=1, order=0)\n",
      "Drift(name=\"DR\", length=0.25, dp=0.0, exact=False, ns=1, order=0)\n"
     ]
    }
   ],
   "source": [
    "# Get by index or name (first matched)\n",
    "\n",
    "print(FODO_B[1])\n",
    "print(FODO_B['DR'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "88c6c89e-d259-478e-af00-9656620e40ec",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['QF_B', 'DR', 'SF_B', 'DR', 'BM', 'DR', 'SD_B', 'DR', 'QD_B', 'QD_B', 'DR', 'SD_B', 'DR', 'BM', 'DR', 'SF_B', 'DR', 'QF_B']\n",
      "['QF_B', 'DR', 'DR', 'DR', 'BM', 'DR', 'SD_B', 'DR', 'QD_B', 'QD_B', 'DR', 'SD_B', 'DR', 'BM', 'DR', 'SF_B', 'DR', 'QF_B']\n",
      "['QF_B', 'DR', 'DR', 'DR', 'BM', 'DR', 'SD_B', 'DR', 'QD_B', 'QD_B', 'DR', 'SD_B', 'DR', 'BM', 'DR', 'DR', 'DR', 'QF_B']\n"
     ]
    }
   ],
   "source": [
    "# Set by index or name (first matched)\n",
    "\n",
    "print(FODO_B.names)\n",
    "\n",
    "FODO_B['SF_B'] = FODO_B['DR']\n",
    "print(FODO_B.names)\n",
    "\n",
    "FODO_B['SF_B'] = FODO_B['DR']\n",
    "print(FODO_B.names)"
   ]
  }
 ],
 "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
}