{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "262a5ec8-2553-4237-ab62-319b6ca22089",
   "metadata": {},
   "source": [
    "# Example-42: Tune (Tune spread from quafrupole errors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "af09c916-090f-4f9a-80ae-e51f29125d4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import\n",
    "\n",
    "from random import random\n",
    "from pprint import pprint\n",
    "\n",
    "import torch\n",
    "from torch import Tensor\n",
    "\n",
    "from pathlib import Path\n",
    "\n",
    "import matplotlib\n",
    "from matplotlib import pyplot as plt\n",
    "matplotlib.rcParams['text.usetex'] = True\n",
    "\n",
    "from model.library.line import Line\n",
    "\n",
    "from model.command.external import load_sdds\n",
    "from model.command.external import load_lattice\n",
    "from model.command.build import build\n",
    "from model.command.tune import tune"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "633298af-f0ce-4ffd-a7e9-9c4d5cfe963f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load ELEGANT twiss\n",
    "\n",
    "path = Path('ic.twiss')\n",
    "parameters, columns = load_sdds(path)\n",
    "\n",
    "nu_qx:Tensor = torch.tensor(parameters['nux'] % 1, dtype=torch.float64)\n",
    "nu_qy:Tensor = torch.tensor(parameters['nuy'] % 1, dtype=torch.float64)\n",
    "\n",
    "psi_qx:Tensor = torch.tensor(parameters['dnux/dp'], dtype=torch.float64)\n",
    "psi_qy:Tensor = torch.tensor(parameters['dnuy/dp'], dtype=torch.float64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e9072e91-882a-4b54-8994-9c163292ad13",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build and setup lattice\n",
    "\n",
    "# Note, sextupoles are turned off and dipoles are linear \n",
    "\n",
    "# Load ELEGANT table\n",
    "\n",
    "path = Path('ic.lte')\n",
    "data = load_lattice(path)\n",
    "\n",
    "# Build ELEGANT table\n",
    "\n",
    "ring:Line = build('RING', 'ELEGANT', data)\n",
    "ring.flatten()\n",
    "\n",
    "# Merge drifts\n",
    "\n",
    "ring.merge()\n",
    "\n",
    "# Turn off sextupoles and set linear dipoles\n",
    "\n",
    "for element in ring: \n",
    "    if element.__class__.__name__ == 'Sextupole':\n",
    "        element.ms = 0.0\n",
    "    if element.__class__.__name__ == 'Dipole':\n",
    "        element.linear = True\n",
    "\n",
    "# Set number of elements of different kinds\n",
    "\n",
    "nb = ring.describe['BPM']\n",
    "nq = ring.describe['Quadrupole']\n",
    "ns = ring.describe['Sextupole']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "efbd9cba-c422-486f-be49-e4b4807889b7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(0.0012, dtype=torch.float64)\n",
      "tensor(0.0008, dtype=torch.float64)\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAHpCAYAAABTH4/7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAlcklEQVR4nO3dQW4aW9o/4DetSB61XU3PWp8ntQQuXkGI9E08g5sV2JlnEORR646seAdJVpDAzEPTK/A1O3BNIvWQYO7I+qTmP8jfdLjGDsRwKOB5JCRTp4ATVyp18uOtc56NRqNRAAAAAEBCf1l1BwAAAADYPkIpAAAAAJITSgEAAACQnFAKAAAAgOSEUgAAAAAkJ5QCAAAAIDmhFAAAAADJPV91Bx7zn//8J/7973/HX//613j27NmquwMAbKHRaBR//PFH/OMf/4i//GU9vs8zhgIAVmnW8VOpQ6l///vfsb+/v+puAADEly9f4n/+539W3Y2ZGEMBAGXwo/FTqUOpv/71rxHx7Q+xu7u74t4AANtoOBzG/v7+eFyyDoyhAIBVmnX8VOpQ6q7cfHd314AKAFipdboNzhgKACiDH42f1mNiBAAAAAA2ilAKAAAAgOSEUgAAAAAkJ5QCAAAAIDmhFAAAAADJzbX6XrfbjYiIwWAQl5eX8erVq6hWqxERURRFdDqdyPM8iqKI4+PjyLLsh20AAAAAbJ+5Qqlmsxn/+te/ol6vR7/fj2azGdfX1+O2q6uriPgWQh0dHUW73f5hGwAAAADbZ67b99rt9rgyKiImKqG+l+f5uKrqsTYAAAAAttNcoVS9Xh//3G634/Xr1xHx7ba+SqUysW+lUoler/doGwAAAADbaa7b9yIier1efPr0KV6+fBnHx8cR8W2OqWn6/f6jbX92e3sbt7e34+fD4XDe7gEAAACwBuZefa9arcbJyUlcX19Hp9N5dN+HAqmH2k5PT2Nvb2/82N/fn7d7AAAAAKyBuUOpiG9zSTWbzWg2mzEYDCLLsnuVT/1+P7Ise7Ttz05OTuLm5mb8+PLly890DwAAAICSmzmU6na78be//W38PM/ziPg2kfn3c019r1arPdr2Zzs7O7G7uzvxAAAAAGDzzDynVKVSmQiYer1eZFk2sRrfnaIoolarjSulHmoDAAAAYDvNHEpVq9V49epVfPjwISIiLi4u4urqatzebrej1WrFwcFBXF5eRrvdnqkNAAAAgO3zbDQajVbdiYcMh8PY29uLm5sbt/IBACuxjuORdewzALA5Zh2L/NRE5wAAAADwFEIpAAAAAJKbeU4pgFU6PLy/7fw8fT8AALaBsReQgkopAAAAAJITSgEAAACQnFAKAAAAgOSEUgAAAAAkJ5QCAAAAIDmhFAAAAADJCaUAAAAASE4oBQAAAEByQikAAAAAkhNKAQAAAJCcUAoAAACA5IRSAAAAACQnlAIAAAAgOaEUAAAAAMkJpQAAAABITigFAAAAQHJCKQAAAACSe77qDgAs2+Hh/W3n5+n7AQAAwH+plAIAAAAgOaEUAAAAAMkJpQAAAABITigFAAAAQHJCKQAAAACSE0oBAAAAkJxQCgAAAIDknq+6AwBlcXh4f9v5efp+AAAAbAOVUgAAAAAkp1IKACCRVqsVJycnkWVZREQURRGdTifyPI+iKOL4+HimNgCATSCUAgBIoNfrxdnZWZycnIy3NZvNuLq6iohvIdTR0VG02+0ftgEAbAK37wEAJFAUReR5PvH8e3meR7fb/WEbAMCmEEoBACxZp9OJRqMxsa3b7UalUpnYVqlUotfrPdo2ze3tbQyHw4kHAEDZCaUAAJZoMBhMnQtqMBhM3b/f7z/aNs3p6Wns7e2NH/v7+z/ZWwCAdIRSAABL9Pnz56jX6zPv/1Ag9VjbyclJ3NzcjB9fvnyZs5cAAOmZ6BwAYEm63W78+uuvU9uyLLtX+dTv9yPLskfbptnZ2YmdnZ2F9BkAIBWhFLBRDg9X3QOASZ8/fx7/XBRFnJ6exqtXr6Jer8f79+/v7V+r1SLP8wfbAFZl2jjr/Dx9P4DNIZQCAFiSP9+29/r163j9+vXEKnx3iqKIWq02rpR6qA0AYFMIpQAAlmwwGMSHDx8iIuLdu3fx+vXrqFar0W63o9VqxcHBQVxeXka73R6/5rE2AIBN8Gw0Go1W3YmHDIfD2Nvbi5ubm9jd3V11d4AVmrVcfNbb92Z9rZJ0YB3HI+vYZ6BcnjKmAph1LGL1PQAAAACSE0oBAAAAkJw5pQAAADaU6QmAMlMpBQAAAEByQikAAAAAkhNKAQAAAJCcUAoAAACA5IRSAAAAACQnlAIAAAAgueer7gCwPSxJDAAAwB2VUgAAAAAkp1IKWKlp1VMAAKRlTAasglAKWFsGTwAAAOvL7XsAAAAAJCeUAgAAACA5oRQAAAAAyQmlAAAAAEhOKAUAAABAckIpAAAAAJJ7vuoOAKybw8P7287P0/cDAABgnQmlAAAANsC0L84AymyuUKrX60W3242IiMvLy/j48WNkWTZui4ioVqtRFEUMBoOoVqsREVEURXQ6ncjzPIqiiOPj4/HrAAAAWE8qyIGnmCuU6na78fbt24iIODs7ixcvXsTV1VVERLx//z4+fPgQERH1ej3a7fb4dc1mc7xfURRxdHQ00Q4AAMBmEFQBs5p5ovNerxenp6fj541GI3q9XhRFERERv/zyS3z9+jW+fv0aFxcX40qou/Y7eZ6Pq60AAAAA2E4zV0pVq9X4+PHj+PlgMIiIiEqlMt427Za8brc7sc/da3q93vj2vju3t7dxe3s7fj4cDmftHgAAAABrZOZKqYhv1VF3Pn36FPV6fRxEDQaD6HQ60el0otVqjSuk7sKrP+v3+/e2nZ6ext7e3vixv78/T/cAAAAAWBM/tfreXQB1N09URExMXp7nebx8+TKur68ffY8/Ozk5iTdv3oyfD4dDwRQAAMACWaUPKIu5KqXutFqtiXmjIibnjrpbZa8oisiy7F5VVL/fn3qr387OTuzu7k48AAAAANg8c1dKnZ2dRavVijzPx9VORVHEixcv4uvXrxP7ViqVqNfr8f79+3vvU6vVfq7HAAn5JhEAAGA55qqU6nQ6Ua1Wx4HU58+fI8uyyPM83r17N96v2+1Go9EYt32vKIqo1WpTK6UAAAAA2A4zV0oVRRHNZnNiW5Zl47mkarVanJ2dRZZlcX19He12e7xfu92OVqsVBwcHcXl5OdEGAAAAwPaZOZTK8zxGo9GD7dVqNarV6oOvvauk+n4FPwAAAAC2009NdA4AAAAATyGUAgAAACA5oRQAAAAAyc08pxTAJjk8XHUPAAAAtptKKQAAAACSE0oBAAAAkJxQCgAAAIDkhFIAAAAAJCeUAgAAACA5oRQAAAAAyQmlAAAAAEhOKAUAAABAckIpAAAAAJITSgEAAACQnFAKAAAAgOSEUgAAAAAkJ5QCAAAAIDmhFAAAAADJCaUAAAAASO75qjsAbKbDw1X3YPWm/Q7Oz9P3AwAAoIxUSgEAAACQnEopgIRUTwEAAHyjUgoAAACA5IRSAAAAACQnlAIAAAAgOaEUAAAAAMkJpQAAAABITigFAAAAQHJCKQAAAACSE0oBAAAAkJxQCgAAAIDknq+6AwCb4PBw1T0AALaJsQewCYRSwFymDYDOz9P3AwAAgPXm9j0AAAAAkhNKAQAAAJCcUAoAAACA5IRSAAAAACQnlAIAAAAgOavvAQAAsFRWcAamEUoBAACQnKAKcPseAAAAAMkJpQAAAABITigFAAAAQHJCKQAAAACSE0oBAAAAkJxQCgAAAIDkhFIAAAAAJCeUAgAAACA5oRQAAAAAyQmlAAAAAEju+ao7AJTX4eGqewAAAMCmUikFAAAAQHJCKQAAAACSE0oBAAAAkJxQCgAAAIDkhFIAAAAAJGf1PQCAJep2uxERMRgM4vLyMl69ehXVajUiIoqiiE6nE3meR1EUcXx8HFmW/bANAGATCKUAAJao2WzGv/71r6jX69Hv96PZbMb19fW47erqKiK+hVBHR0fRbrd/2AZsl8PDVfcAYDncvgcAsETtdntcGRURE5VQ38vzfFxV9VgbAMCmUCkFALBE9Xp9/HO73Y7Xr19HxLfb+iqVysS+lUoler1e/P777w+2fR9w3bm9vY3b29vx8+FwuMg/AgDAUqiUAgBYsl6vF61WK16+fBnHx8cR8W2OqWn6/f6jbdOcnp7G3t7e+LG/v7+IbgMALJVQCgBgyarVapycnMT19XV0Op1H930okHqs7eTkJG5ubsaPL1++PKG3AABpzHX7Xq/XG89ncHl5GR8/fpxphRirxwAA2y7Lsmg2m/Hy5cv4+vVrZFl2r/Kp3+9HlmWPtk2zs7MTOzs7y+o6AMBSzFUp1e124+3bt/H27ds4ODiIFy9ejNuazWa8ffs2Go1GNBqNODo6mqkNAGBTdbvd+Nvf/jZ+nud5RHz7wu77uaa+V6vVHm0DANgUM1dK9Xq9OD09jbdv30ZERKPRiFardW91mAirxwAARHybnPz7gKnX60WWZVMnKy+KImq12rhS6qE2YLMdHq66BwDpzBxKVavV+Pjx4/j53ZwGlUolPn/+vJDVY6wcAwBskmq1Gq9evYoPHz5ERMTFxUVcXV2N29vtdrRarTg4OIjLy8tot9sztQEAbIK55pRqNBrjnz99+hT1ej2yLFvY6jGnp6fx22+/zdMlAIBS+378dLfy3p08z+Pdu3f39vtRGwDAJvip1fcGg0F0Op0ffmM37+oxVo4BAAAA2A5zVUrdabVacXFxMZ7XYFGrx1g5BgAAAGA7zF0pdXZ2Fq1WK/I8j8FgEIPBwOoxAAAAAMxlrkqpTqcT1Wp1HEh9/vw5jo+PH10hxuoxsPmsEgMAwCJMG1een6fvB5DGzKFUURTRbDYntmVZNp6w0+oxAAAAAMzq2Wg0Gq26Ew8ZDoext7cXNzc3sbu7u+ruwNZRAZWGb/+g3NZxPLKOfQa+Mf66z1gJ1s+sY5GfWn0PAAAAAJ5CKAUAAABAckIpAAAAAJKba/U9AFbHajQAAMAmUSkFAAAAQHJCKQAAAACSE0oBAAAAkJxQCgAAAIDkhFIAAAAAJGf1PSAipq/sBgAAAMsilAJYsWmB4Pl5+n4AAACk5PY9AAAAAJJTKQVQQm6nBAAANp1KKQAAAACSE0oBAAAAkJxQCgAAAIDkhFIAAAAAJCeUAgAAACA5oRQAAAAAyQmlAAAAAEhOKAUAAABAckIpAAAAAJITSgEAAACQnFAKAAAAgOSEUgAAAAAk93zVHQDg5x0e3t92fp6+HwAAAPMSSgEAAKzAtC+XALaJUAoAAIDSUhkOm8ucUgAAAAAkJ5QCAAAAIDmhFAAAAADJCaUAAAAASE4oBQAAAEByVt8DAABgrViRDzaDSikAAAAAklMpBbBhfHMIAACsA5VSAAAAACQnlAIAAAAgOaEUAAAAAMkJpQAAAABITigFAAAAQHJW34MNZyU2AAAAykilFAAAAADJCaUAAAAASE4oBQAAAEByQikAAAAAkjPROWyhaZOfAwAAQEoqpQAAAABITigFAAAAQHJCKQAAAACSE0oBAAAAkJxQCgAAAIDkrL4HsKWmrcJ4fp6+HwAAwHZSKQUAAABAckIpAAAAAJJz+x5skGm3YwEAAEAZqZQCAAAAIDmhFAAAAADJCaUAAAAASE4oBQAAAEByQikAAAAAkpsrlOr1evHLL79M3d7r9SIioiiK8c93z8/OzqLT6cTZ2VkMBoOn9RgAAACAtfd81h07nU7keT4RON15//59fPjwISIi6vV6tNvtcVuz2Yyrq6uI+BZQHR0dTbQDAAAAsH1mDqUajcaDbb/88kt8/fo1IiKyLBtvL4piYr88z6Pb7c7ZRQAAAAA2zcyh1I98H0bd6Xa7UalUJrZVKpXo9XpRrVbv7X97exu3t7fj58PhcFHdAwAAAKBEFjLR+WAwiE6nE51OJ1qt1rhC6qH5o/r9/tTtp6ensbe3N37s7+8vonsAAAAAlMxCKqWOj4/HlVJ5nsfLly/j+vr6wf0fCqtOTk7izZs34+fD4VAwBQAAALCBFhJKFUUxvh0vz/MoiiKKoogsy+5VRfX7/am3+kVE7OzsxM7OziK6BMB3Dg9X3QMAAIBJT759r9frxYsXL+5tr1QqUa/Xp76mVqs99WMBAAAAWGM/VSk1GAwmbtd79+7duK3b7Uaj0Ygsy+5VRBVFEbVa7cFKKQAAgE2kahngvplDqW63GxcXFxHxbULyg4ODcfhUq9Xi7OwssiyL6+vraLfb49e12+1otVpxcHAQl5eXE20AAAAAbKdno9FotOpOPGQ4HMbe3l7c3NzE7u7uqrsDpecbOJ7q/HzVPYDyWcfxyDr2GTadcdryGcdAecw6FnnynFIAAAAAMC+hFAAAAADJCaUAAAAASE4oBQAAAEByQikAAAAAkhNKAQAAAJDc81V3AAAAYJMcHq66BwDrQaUUAAAAAMkJpQAAAABIzu17AABL1Ov1otvtRkTE5eVlfPz4MbIsi4iIoiii0+lEnudRFEUcHx/P1AYAsAmEUgAAS9TtduPt27cREXF2dhYvXryIq6uriIhoNpvjn4uiiKOjo2i32z9sAwDYBG7fAwBYkl6vF6enp+PnjUYjer1eFEURRVFM7Jvn+bii6rE2AIBNIZQCAFiSarUaHz9+HD8fDAYREVGpVKLb7UalUpnYv1KpjG/3e6htmtvb2xgOhxMPAICyE0oBACxRo9EY//zp06eo1+uRZdk4oPqzfr//aNs0p6ensbe3N37s7+8/tdsAAEsnlAIASGAwGESn0/nhvFAPBVKPtZ2cnMTNzc348eXLlyf0FAAgDROdAwAk0Gq14uLiYryCXpZl9yqf+v1+ZFn2aNs0Ozs7sbOzs4xuAz9weLjqHgCsL5VSAABLdnZ2Fq1WK/I8j8FgEIPBIOr1+tR9a7Xao20AAJtCKAUAsESdTieq1eo4kPr8+XNkWRZ5nk/sVxRF1Gq1H7YBAGwKt+8BACxJURTRbDYntmVZFsfHxxER0W63o9VqxcHBQVxeXk7MN/VYGwDAJng2Go1Gq+7EQ4bDYezt7cXNzU3s7u6uujtQeuY04KnOz1fdAyifdRyPrGOfYV0Zf5WHcQyUx6xjEbfvAQAAAJCc2/dgTflWDgAAgHWmUgoAAACA5IRSAAAAACQnlAIAAAAgOaEUAAAAAMkJpQAAAABIzup7AAAArL1pq1Ofn6fvBzA7lVIAAAAAJCeUAgAAACA5oRQAAAAAyQmlAAAAAEjOROcAjJkgFAAASEWlFAAAAADJCaUAAAAASM7te7AGpt1SBQAAAOtMpRQAAAAAyQmlAAAAAEhOKAUAAABAcuaUAgAAmIF5PgEWS6UUAAAAAMkJpQAAAABITigFAAAAQHJCKQAAAACSE0oBAAAAkJxQCgAAAIDkhFIAAAAAJCeUAgAAACA5oRQAAAAAyQmlAAAAAEhOKAUAAABAckIpAAAAAJITSgEAAACQ3PNVdwAAAACW4fDw/rbz8/T9AKYTSgEAALA1BFVQHm7fAwAAACA5oRQAAAAAyQmlAAAAAEhOKAUAAABAckIpAAAAAJITSgEAAACQnFAKAAAAgOSez7Nzr9eLo6OjuLq6mtheFEV0Op3I8zyKoojj4+PIsuyHbQAAAABsp5lDqbtgqdfr3WtrNpvjoKooijg6Oop2u/3DNgAAAAC208yhVKPRmLq9KIqJ53meR7fb/WEbAAAAANtrrtv3pul2u1GpVCa2VSqV6PV68fvvvz/YVq1W773X7e1t3N7ejp8Ph8Ondg8AAACAEnryROeDwWDq9n6//2jbNKenp7G3tzd+7O/vP7V7AAAAAJTQkyulHvJQIPVY28nJSbx582b8fDgcCqYAVuzw8P628/P0/QAAADbLk0OpLMvuVT71+/3IsuzRtml2dnZiZ2fnqV0CAAAAoOSefPtevV6fur1Wqz3aBgAAAMD2+qlKqcFgMK52yvN8oq0oiqjVauNKqYfaAAAAANheM4dS3W43Li4uIuLbhOQHBwfRaDQiIqLdbker1YqDg4O4vLyMdrs9ft1jbQAAAABsp2ej0Wi06k48ZDgcxt7eXtzc3MTu7u6quwMrM22iaVglE52zTdZxPLKOfYZ1YEy2uYxtYLFmHYs8eU4pAAAAAJiXUAoAAACA5IRSAAAAACQnlAIAAAAguZlX3wPSMIEmAAAA20ClFAAAAADJCaUAAAAASE4oBQAAAEByQikAAAAAkhNKAQAAAJCcUAoAAACA5IRSAAAAACQnlAIAAAAgOaEUAAAAAMkJpQAAAABITigFAAAAQHJCKQAAAACSE0oBAAAAkNzzVXcAttnh4ap7AAAAAKshlAJgbtMC1fPz9P0AAADWl1AKElIZBQAA5eMLN1gNoRQAAAD8iaAKlk8oBcBC/HngZtAGAAA8xup7AAAAACQnlAIAAAAgOaEUAAAAAMmZUwoAANhqJrQGWA2VUgAAAAAkJ5QCAAAAIDmhFAAAAADJCaUAAAAASM5E5wAshUljAVhn065jACyWSikAAAAAkhNKAQAAAJCcUAoAAACA5MwpBQCwRL1eL46OjuLq6mpie1EU0el0Is/zKIoijo+PI8uyH7YBAGwKoRQAwJLcBUu9Xu9eW7PZHAdVRVHE0dFRtNvtH7YBAGwKoRQAwJI0Go2p24uimHie53l0u90ftgEAbBKhFABAYt1uNyqVysS2SqUSvV4vfv/99wfbqtXq1Pe7vb2N29vb8fPhcLj4TgMALJhQChbg8HDVPYD1MO1cOT9P3w9YtcFgMHV7v99/tO0hp6en8dtvvy2gZwAA6Vh9DwCgJB4KpH7UdnJyEjc3N+PHly9fFt85AIAFUykFAJBYlmX3Kp/6/X5kWfZo20N2dnZiZ2dnGV0FAFgalVIAAInV6/Wp22u12qNtAACbRKUUAEACg8FgXO2U5/lEW1EUUavVxpVSD7UBAGwSoRQAwJJ0u924uLiIiG+TkR8cHESj0YiIiHa7Ha1WKw4ODuLy8jLa7fb4dY+1AQBsimej0Wi06k48ZDgcxt7eXtzc3MTu7u6quwMPsvoe/Dyr71F26zgeWcc+wyoZyzEr4xaYzaxjEXNKAQAAAJCc2/dgTr5JAwAAgKcTSgEAAMAMpn1B7ZY++HlCKQBWyuAOAAC2kzmlAAAAAEhOKAUAAABAckIpAAAAAJITSgEAAACQnFAKAAAAgOSEUgAAAAAkJ5QCAAAAIDmhFAAAAADJCaUAAAAASE4oBQAAAEByz1fdAQAAgGU4PLy/7fw8fT8AmG5hoVSv14uIiGq1GkVRxGAwiGq1GhERRVFEp9OJPM+jKIo4Pj6OLMsW9dEAAAAArJmFhVLv37+PDx8+REREvV6Pdrs9bms2m3F1dRUR3wKqo6OjiXYAAAAAtsvCQqlffvklvn79GhExUQVVFMXEfnmeR7fbXdTHAgAAALCGFjqn1LRb8rrdblQqlYltlUoler3e+Pa+O7e3t3F7ezt+PhwOF9k9mNu0eQgAAACAp1tYKDUYDKLT6URExOXlZbx+/TryPI/BYDB1/36/f2/b6elp/Pbbb4vqEgAAwARfOgKUx8JCqe8nL8/zPF6+fBnX19cP7j8trDo5OYk3b96Mnw+Hw9jf319UFwEAAAAoib8s6o2+nzvqbpW9oigiy7J7VVH9fn/qrX47Ozuxu7s78QAAAABg8ywklOr1evHixYt72yuVStTr9amvqdVqi/hoAAAAANbQQm7fy/M83r17N37e7Xaj0WhElmX3KqKKooharTa1UgoAAACA7bCQUCrLsqjVanF2dhZZlsX19XW02+1xe7vdjlarFQcHB3F5eTnRBgAAAOtq2uT55+fp+wHraGETnVer1ahWq1Pbvq+kajQai/pIAAAAANbUwiY6BwAAAIBZCaUAAAAASE4oBQAAAEByQikAAAAAkhNKAQAAAJCcUAoAAACA5IRSAAAAACQnlAIAAAAgOaEUAAAAAMk9X3UHoCwOD1fdAwAAANgeKqUAAAAASE4oBQAAAEBybt8DoHRmvZ32/Hy5/QAAAJZHpRQAAAAAyamUAgAA1p5FawDWj0opAAAAAJITSgEAAACQnNv3ANgo027fMCE6AACUj0opAAAAAJITSgEAAACQnNv3AAAAYIFMJwCzEUoBsLYs/w0AAOvL7XsAAAAAJCeUAgAAACA5oRQAAAAAyZlTCoCNZ7JRAAAoH6EUAAAALJkvyeA+oRQbzz/+AAAAUD5CKTaK5eEBAABgPZjoHAAAAIDkhFIAAAAAJCeUAgAAACA5oRQAAAAAyZnonK1kQnQAAABYLaEUAACwVnzBCLAZ3L4HAAAAQHIqpVhbviEDANh8xnwAm0soxVowGAEAADbNtP/nnJ+n7wesilAKgK1kEAgAAKtlTikAAAAAkhNKAQAAAJCcUAoAAACA5IRSAAAAACRnonMAAKAUrLgMFmNhuwilAOD/MwgEAIB03L4HAAAAQHIqpSgdZdvAOlJlBTAfYz4AVEoBAAAAkJxQCgAAAIDk3L4HAAAAa8bUAWwCoRQAPMKADwAAlkMoxUqZ4BLYZAItAGAR/L+JTSWUAgAAlsp/qAGYRijFk6kEAAAAAOYllAKAOfnGHwAAnk4oBQAAABvAXSysG6EUyagsAAAASEtQRZkJpVgKARQAAADwGKEUAACwML6cBGBWQikepMwTYPGe8m+rf5cBANgkQinm4psvAIDtZBwIm8MXXZTFX1bdAQAAAAC2T5JKqaIootPpRJ7nURRFHB8fR5ZlKT56qyz6lhAA0vBtJdMYPwEAmy5JKNVsNuPq6ioivg2wjo6Oot1up/jorec/OgCb7SlfKrgelJvxE/Mw5gNgHS09lCqKYuJ5nufR7XaX/bEAAGvL+IlVU0UP22fR4bawnFksPZTqdrtRqVQmtlUqlej1elGtVie2397exu3t7fj5zc1NREQMh8Nld/On/frrbPt9/vzz7zfra//v/2bbb9qvc9bXApDG//7v8j/jKZfXRV//Zv2Mp7zfz7obh4xGo2SfOc/4KWI9x1As1rSx3KyHf9bzGdg+08YjT/n/qcvS9ph1/LT0UGowGEzd3u/37207PT2N33777d72/f39RXcrub291bw2xfsBsJ5SXA826Rr2xx9/xF6iDswzforY7DEUP8+YD1iGMv3flvL70fgpyZxS00wbbJ2cnMSbN2/Gz//zn/9Ev9+Pv//97/Hs2bOEvVs/w+Ew9vf348uXL7G7u7vq7mw9x6M8HIvycCzKxfGY3Wg0ij/++CP+8Y9/rLorD4ZVxlA/z7lQHo5FeTgW5eJ4lIdjMbtZx09LD6WyLLv3rV6/35+6eszOzk7s7Ozcez2z293ddXKUiONRHo5FeTgW5eJ4zCZVhdSdecZPEcZQi+BcKA/Hojwci3JxPMrDsZjNLOOnvyy7E/V6fer2Wq227I8GAFhLxk8AwDZYeiiV5/nE86Ioolar+fYOAOABxk8AwDZIMqdUu92OVqsVBwcHcXl5Ge12O8XHbpWdnZ345z//ea90n9VwPMrDsSgPx6JcHI/yM35Kw7lQHo5FeTgW5eJ4lIdjsXjPRinXNwYAAACASHD7HgAAAAD8mVAKAAAAgOSEUgAAAAAkl2Sic56mKIrodDqR53kURRHHx8cPrr7z2L6zvk+r1YqTkxMr/EyR4lj0er3odrsREXF5eRkfP350LCLN736ez9h2zoXycI2A6Zwb5eK6sTrGUOXhPCgX14mSGFF61Wp1/PP19fWo0Wj81L6zvM/V1dUoIkZfv359Yq83U4pj8e7du4mfv993m6X43c/zGdvOuVAerhEwnXOjXFw3VscYqjycB+XiOlEOQqmSu76+vvcPSZZlc+876/u02+1RnudOlilSHIurq6uJ97y+vh5FxOj6+vrJ/V9nKX7383zGtnMulIdrBEzn3CgX143VMYYqD+dBubhOlIc5pUqu2+1GpVKZ2FapVKLX68217yzv0+l0otFoLLD3myXFsahWq/Hx48fx9sFgMG7fZil+9/N8xrZzLpSHawRM59woF9eN1TGGKg/nQbm4TpSHUKrk7v4h+bN+vz/Xvj96n8Fg4N7WH0h1LL7/B+vTp09Rr9e3/tik+N3P8xnbzrlQHq4RMJ1zo1xcN1bHGKo8nAfl4jpRHkKpNfXQX/55971r+/z5c9Tr9ad1akst+lh8/7zT6US73f65jm2BZf3uf/Yztp1zoTxcI2A650a5uG6sjjFUeTgPysV1Ij2r763Ihw8f4vr6+sH2ly9fjlPtP6e1/X5/atr62L6PtXW73fj1119//g+z5sp0LL7XarXi4uJCsh6P/z7n2XcRx4U0x+N7zoWHuUawbcp0zXZulOt4fM9147+MocrDeVAurhMlsupJrXjcQxOnTZsg7bF9H2u7uLgYvX//fvyIiNHbt29HV1dXC//zrLMUx+LOu3fvxhMSfv36desnxEvxu5/nM7adc6E8XCNgOudGubhurI4xVHk4D8rFdaI8VEqVXJ7nE8+LooharTZOcHu9XmRZFnmeP7rvnxPf79v+XE74+vXreP369b3323YpjkXEt4nwqtVq5Hkeg8EgPn/+HMfHx8v6Y62FFL/7Hx0X/su5UB6uETCdc6NcXDdWxxiqPJwH5eI6USKrTsX4sevr69Hbt29H7XZ79Pbt24n0ttFojN69ezfTvo+1jUbfUvR3796NImJ0fHwswZ1i2cfibtnW7x+ZJXVHo1Ga8+BH5wj/5VwoD9cImM65US6uG6tjDFUezoNycZ0oh2ej0WiUNgYDAAAAYNtZfQ8AAACA5IRSAAAAACQnlAIAAAAgOaEUAAAAAMkJpQAAAABITigFAAAAQHJCKQAAAACSE0oBAAAAkJxQCgAAAIDkhFIAAAAAJCeUAgAAACC5/weSkKeLZbpF8wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Compute tune spread using MC\n",
    "\n",
    "kns = 5.0E-3*torch.randn((8192, nq), dtype=torch.float64)\n",
    "\n",
    "nuqx, nuqy = torch.vmap(lambda kn: tune(ring, [kn], ('kn', ['Quadrupole'], None, None), matched=True, limit=1, epsilon=None), chunk_size=1024)(kns).T\n",
    "\n",
    "# Spread\n",
    "\n",
    "print(nuqx.std(-1))\n",
    "print(nuqy.std(-1))\n",
    "print()\n",
    "\n",
    "# Plot samples\n",
    "\n",
    "fig, (ax, ay) = plt.subplots(1, 2, figsize=(12, 5))\n",
    "ax.hist((nuqx - nuqx.mean()).cpu().numpy(), bins=100, range=(-5.0E-3, +5.0E-3), color='blue', alpha=0.7)\n",
    "ay.hist((nuqy - nuqy.mean()).cpu().numpy(), bins=100, range=(-5.0E-3, +5.0E-3), color='blue', alpha=0.7)\n",
    "plt.tight_layout() \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "20b74300-2cc6-44c1-ab1f-48138f07c3e2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(0.0012, dtype=torch.float64)\n",
      "tensor(0.0008, dtype=torch.float64)\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAHqCAYAAADVi/1VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAArv0lEQVR4nO3dP24bW5o34NeNCyhqqS4na4yTmh3Q1ApMA1/ijGyvQFJ+AxOKBo4EaQe2V2CTmUKxV8AWd6BKBHRIU7wRMUDzCzyqES1Klizx8N/zAARUdYrkkeuWeO6P7zn1YjKZTAIAAAAAEvrLojsAAAAAwOYRSgEAAACQnFAKAAAAgOSEUgAAAAAkJ5QCAAAAIDmhFAAAAADJCaUAAAAASE4oBQAAAEByvy26A/f597//Hf/617/ir3/9a7x48WLR3QEANtBkMok///wz/va3v8Vf/rIa3+cZQwEAi/TQ8dNSh1L/+te/4uXLl4vuBgBAXF5exn/+538uuhsPYgwFACyDn42fljqU+utf/xoR33+J7e3tBfcGANhEo9EoXr58WY5LVoExFACwSA8dPy11KHVdbr69vW1ABQAs1CpNgzOGAgCWwc/GT6uxMAIAAAAAa0UoBQAAAEByQikAAAAAkhNKAQAAAJCcUAoAAACA5IRSAAAAACQnlAIAAAAgud8ec3C3242IiOFwGL1eL969exfVajUiIoqiiE6nE3meR1EUsb+/H1mW/bQNAAAAgM3zqFCq2WzGP/7xj6jX6zEYDKLZbMbFxUXZdn5+HhHfQ6i9vb1ot9s/bQMAAABg8zxq+l673S4royJiqhLqpjzPy6qq+9oAAAAA2EyPCqXq9Xr5c7vdjoODg4j4Pq2vUqlMHVupVKLf79/bBgAAAMBmetT0vYiIfr8fX758iTdv3sT+/n5EfF9japbBYHBvGwAAAACb6dF336tWq3F4eBgXFxfR6XTuPfauQOqutvF4HKPRaOoBAAAAwPp5dCgV8X0tqWazGc1mM4bDYWRZdqvyaTAYRJZl97b96OjoKHZ2dsrHy5cvf6V7AAAAACy5B4dS3W43fv/993I7z/OI+L6Q+c21pm6q1Wr3tv3o8PAwrq6uysfl5eVDuwcAAADACnnwmlKVSmUqYOr3+5Fl2dTd+K4VRRG1Wq2slLqr7UdbW1uxtbX18N4DAAAAsJIeHEpVq9V49+5dfPr0KSIizs7O4vz8vGxvt9vRarVid3c3er1etNvtB7UBAAAAsHleTCaTyaI7cZfRaBQ7OztxdXUV29vbi+4OALCBVnE8sop9BgDWx0PHIg+ulAJYVW/f3t53epq+HwAAq8L4CUjhl+6+BwAAAABPIZQCAAAAIDmhFAAAAADJWVMKACCRVqsVh4eHkWVZREQURRGdTifyPI+iKGJ/f/9BbQAA60AoBQCQQL/fj5OTkzg8PCz3NZvNOD8/j4jvIdTe3l602+2ftgEArAOhFABAAkVRRJ7nU9s35Xke3W73p20A8zDrbnsA82ZNKQCAOet0OtFoNKb2dbvdqFQqU/sqlUr0+/172wAA1oVKKQCAORoOhzPXghoOhzOPHwwG97bNMh6PYzwel9uj0eix3QQASE6lFADAHH39+jXq9fqDj78rkLqv7ejoKHZ2dsrHy5cvH9lLAID0hFIAAHPS7Xbj73//+8y2LMtuVT4NBoPIsuzetlkODw/j6uqqfFxeXj5L/wEA5sn0PQCAOfr69Wv5c1EUcXR0FO/evYt6vR4fP368dXytVos8z+9sm2Vrayu2traer9MAAAkIpQAA5uTHaXsHBwdxcHAwdRe+a0VRRK1WKyul7moDAFgXQikAgDkbDofx6dOniIg4Pj6Og4ODqFar0W63o9Vqxe7ubvR6vWi32+Vz7msDAFgHLyaTyWTRnbjLaDSKnZ2duLq6iu3t7UV3B1hRb9/e3nd6mr4fwGpaxfHIKvYZWKxZ46UfGT8BD/XQsYiFzgEAAABITigFAAAAQHJCKQAAAACSE0oBAAAAkJxQCgAAAIDkhFIAAAAAJCeUAgAAACA5oRQAAAAAyQmlAAAAAEhOKAUAAABAckIpAAAAAJITSgEAAACQnFAKAAAAgOR+W3QHAB7i7dvb+05P0/cDAGBTGY8Bz02lFAAAAADJCaUAAAAASM70PQAAgDVlyh2wzFRKAQAAAJCcUAoAAACA5IRSAAAAACRnTSkAAIA1MGv9KIBlJpQClo4BFQAAwPozfQ8AAACA5IRSAAAAACRn+h6wUKbqAQAAbCaVUgAAAAAkp1IK4B6zKrlOT9P3AwAAYN0IpQAAAPglvsADnsL0PQAAAACSE0oBAAAAkJxQCgAAAIDkrCkFrKxlWsNgmfoCALBIxkXAQ6mUAgAAACA5lVIAAADMleopYBahFDAXBh4AAADcx/Q9AAAAAJITSgEAAACQnOl7wFqZNW0wxXMBAAB4HJVSAAAAACSnUgoAAGCDqA4HloVKKQAAAACSUykFADBH3W43IiKGw2H0er149+5dVKvViIjo9/sREVGtVqMoihgOh2VbURTR6XQiz/MoiiL29/cjy7KF/A4AAPOgUgoAYI6azWZUKpVoNBrxX//1X9FsNsu2jx8/xqtXr+LFixdxcHAQeZ5PPe/9+/fRaDSi0WjE3t7eIroPADA3KqUAAOao3W6X1U8RMVXt9OrVq/j27dut/UVRTL1GnudlxRUAwLoQSgE8ksVBgceo1+vlz+12Ow4ODqbaZ03J63a7UalUpvZVKpXo9/tTARcAwCoTSgEAzFm/348vX77EmzdvYn9/v9w/HA6j0+lERESv1yun8A2Hw5mvMxgMZu4fj8cxHo/L7dFo9HydBwCYk0eFUv1+vywd7/V68fnz5/LbPQt1AgDMVq1WI8/zaLVa0el0otFoRERMjYnyPI83b97ExcXFna9zV1h1dHQUHz58eO5uAwDM1aMWOu92u/H+/ft4//597O7uxuvXr8s2C3UCANwty7JoNpvRbDbLcOnm2lHXX94VRRFZlt2qihoMBnd+qXd4eBhXV1fl4/Lycl6/BgDAs3lwKNXv9+Po6KjcbjQa0e/3y8HU9UKd3759i7Ozs3LQZKFOAGBTdbvd+P3338vt6y/tiqKIfr8/9QXftUqlMrUO1U21Wm3m/q2trdje3p56AAAsuwdP36tWq/H58+dy+/obvpuLcFqoEwDg//wYMPX7/ciyLKrVagyHwzg+Pi7but1uNBqNyLLs1piqKIqo1WqWPwDWyqybx5yepu8HsDiPWlPqev2DiIgvX75EvV4vB0fPtVAnAMC6qFar8e7du/j06VNERJydncX5+XlEfP8yr1arxcnJSWRZFhcXF9Fut8vnttvtaLVasbu7G71eb6oNAGAd/NLd964DqOtBVcTzLNTpzjEAwLq5+aXezTvvRXwPre6qHM/zvKykuvkaAADr4lELnV9rtVpT60ZFPM9CnUdHR7Gzs1M+Xr58+SvdAwAAAGDJPbpS6uTkJFqt1tTUvKIo4vXr1/Ht27epY6/XUfj48eOt15m1UOfh4WH88ccf5fZoNBJMAXMxaw0DAAAA0nlUpVSn04lqtVoGUl+/fo0sy6bKyyOmF+q8vsvMtfsW6nTnGAAAAIDN8OBKqaIootlsTu3LsqxcS8pCnQAAAAA81INDqTzPYzKZ3NluoU4AAAAAHuqXFjoHAAAAgKcQSgEAAACQnFAKAAAAgOSEUgAAAAAkJ5QCAAAAIDmhFAAAAADJCaUAAAAASE4oBQAAAEByQikAAAAAkhNKAQAAAJDcb4vuALBa3r69ve/0NH0/AAAAWG0qpQAAAABITigFAAAAQHJCKQAAAACSs6YUAADAipm1zifAqlEpBQAAAEByQikAAAAAkhNKAQAAAJCcUAoAAACA5IRSAAAAACQnlAIAAAAgOaEUAAAAAMkJpQAAAABI7rdFdwBgXb19+7DjTk/n2w8AAIBlpFIKAAAAgORUSgEAALAUZlWaqyqH9aVSCgAAAIDkhFIAAAAAJCeUAgAAACA5oRQAAAAAyQmlAAAAAEhOKAUAAABAcr8tugMAAADc7e3bRfcAYD5USgEAAACQnEop4Ml8ewcAAMBjqZQCAAAAIDmVUgAAc9TtdiMiYjgcRq/Xi3fv3kW1Wo2IiKIootPpRJ7nURRF7O/vR5ZlP20DAFgHQikAgDlqNpvxj3/8I+r1egwGg2g2m3FxcVG2nZ+fR8T3EGpvby/a7fZP2wAA1oHpewAAc9Rut8vKqIiYqoS6Kc/zsqrqvjYAgHWhUgq4kwXMAZ6uXq+XP7fb7Tg4OIiI79P6KpXK1LGVSiX6/X7885//vLPtZsAFALDKhFIAAHPW7/fjy5cv8ebNm9jf34+I72tMzTIYDO5tm2U8Hsd4PC63R6PRk/oLLI4vBYFNIpQCWLBZg8/T0/T9AOanWq1GnufRarWi0+lEo9G489i7Aqn72o6OjuLDhw9P7CUAQFrWlAIASCDLsmg2m9FsNmM4HEaWZbcqnwaDQWRZdm/bLIeHh3F1dVU+Li8v5/VrAAA8G6EUAMCcdLvd+P3338vtPM8j4vtC5jfXmrqpVqvd2zbL1tZWbG9vTz0AAJad6XsAAHNSqVSmAqZ+vx9Zls1crLwoiqjVamWl1F1tAADrQigFJGPhTmDTVKvVePfuXXz69CkiIs7OzuL8/Lxsb7fb0Wq1Ynd3N3q9XrTb7Qe1AQCsgxeTyWSy6E7cZTQaxc7OTlxdXSlDhwUQIi2Ohc5heazieGQV+wx8Z/x1m3ERrJ6HjkWsKQUAAABAckIpAAAAAJITSgEAAACQnFAKAAAAgOSEUgAAAAAkJ5QCAAAAIDmhFAAAAADJCaUAAAAASE4oBQAAAEByQikAAAAAkhNKAQAAAJCcUAoAAACA5IRSAAAAACQnlAIAAAAgOaEUAAAAAMn99piD+/1+dLvdiIjo9Xrx+fPnyLIsIiKKoohOpxN5nkdRFLG/v/+gNgAAAAA2z6NCqW63G+/fv4+IiJOTk3j9+nWcn59HRESz2Sx/Looi9vb2ot1u/7QNAAAAgM3z4Ol7/X4/jo6Oyu1GoxH9fj+KooiiKKaOzfO8rKi6rw0AAACAzfTgUKparcbnz5/L7eFwGBERlUolut1uVCqVqeMrlUo53e+uNgAAAAA206Om7zUajfLnL1++RL1ejyzLyoDqR4PB4N62H43H4xiPx+X2aDR6TPcAAAAAWBG/dPe94XAYnU7np+tC3RVI3dV2dHQUOzs75ePly5e/0j0AAAAAltyjKqWutVqtODs7K++gl2XZrcqnwWAQWZbd2/ajw8PD+OOPP8rt0WgkmIJE3r5ddA8AAADYJI+ulDo5OYlWqxV5nsdwOIzhcBj1en3msbVa7d62H21tbcX29vbUAwAAAID186hQqtPpRLVaLQOpr1+/RpZlkef51HFFUUStVvtpGwAAAACb6cHT94qiiGazObUvy7LY39+PiIh2ux2tVit2d3ej1+tNrTd1XxsAAAAAm+fFZDKZLLoTdxmNRrGzsxNXV1em8sGcWVNquZyeLroHwLVVHI+sYp+B74zJbjMugtXz0LHIL919DwAAAACe4pfuvgcAAAApzKoeUz0F60GlFAAAAADJCaUAAAAASE4oBQAAAEByQikAAAAAkrPQOQAAwALMWsAbYJOolAIAAAAgOaEUAAAAAMkJpQAAAABITigFAAAAQHJCKQAAAACSE0oBAAAAkJxQCgAAAIDkflt0BwC47e3b2/tOT9P3AwAAYF6EUgAAc9Tv96Pb7UZERK/Xi8+fP0eWZWVbRES1Wo2iKGI4HEa1Wo2IiKIootPpRJ7nURRF7O/vl88DAFgHQikAgDnqdrvx/v37iIg4OTmJ169fx/n5eUREfPz4MT59+hQREfV6Pdrtdvm8ZrNZHlcURezt7U21AwCsOmtKAQDMSb/fj6Ojo3K70WhEv9+PoigiIuLVq1fx7du3+PbtW5ydnZWVUNft1/I8L6utAADWhVAKAGBOqtVqfP78udweDocREVGpVMp9WZbdmpbX7Xanjrl+zvV0PwCAdWD6HgDAHDUajfLnL1++RL1eL0Oo4XAYnU4nIr6vN3VwcBB5npfh1Y8Gg8G8uwsAkIxQCmCFuUsfrI7rAOp6naiImFq8PM/zePPmTVxcXNz7GrOMx+MYj8fl9mg0epY+AwDMk+l7AAAJtFqtqXWjIqbXjrq+y15RFJFl2a2qqMFgcOfd946OjmJnZ6d8vHz5ch6/AgDAsxJKAQDM2cnJSbRarXJq3nA4jH6/H69fv751bKVSiXq9PvN1arXazP2Hh4dxdXVVPi4vL5+1/wAA8yCUAgCYo06nE9VqtQykvn79GlmWRZ7ncXx8XB7X7Xaj0WiUbTcVRRG1Wu3OSqmtra3Y3t6eegAALDtrSgEAzElRFNFsNqf2ZVlWriVVq9Xi5OQksiyLi4uLaLfb5XHtdjtarVbs7u5Gr9ebagMAWAdCKQCAOcnzPCaTyZ3t1Wo1qtXqnc+9rqS6eQc/AIB1IZQCAABgpbgDMawHa0oBAAAAkJxQCgAAAIDkhFIAAAAAJCeUAgAAACA5oRQAAAAAyQmlAAAAAEhOKAUAAABAckIpAAAAAJITSgEAAACQnFAKAAAAgOSEUgAAAAAkJ5QCAAAAIDmhFAAAAADJCaUAAAAASE4oBQAAAEByQikAAAAAkvtt0R0A4GHevl10DwAAAJ6PUAo2kHADACAt4y+A20zfAwAAACA5oRQAAAAAyQmlAAAAAEhOKAUAAABAckIpAAAAAJITSgEAAACQnFAKAAAAgOSEUgAAAAAkJ5QCAAAAIDmhFAAAAADJCaUAAAAASE4oBQAAAEByQikAAAAAkvtt0R0AAACAp3r79va+09P0/QAeTqUUAAAAAMmplAIAAHhGsyp2ALjtUZVS/X4/Xr16NXN/v9+PiIiiKMqfr7dPTk6i0+nEyclJDIfDp/UYAAAAgJX34EqpTqcTeZ5PBU7XPn78GJ8+fYqIiHq9Hu12u2xrNptxfn4eEd8Dqr29val2AAAAADbPg0OpRqNxZ9urV6/i27dvERGRZVm5vyiKqePyPI9ut/vILgJPoXwcAACAZfRsC51nWTYVSEVEdLvdqFQqU/sqlcrMaisAAAAANsezLHQ+HA6j0+lERESv14uDg4PI8/zO9aMGg8HM/ePxOMbjcbk9Go2eo3sAAAAALJlnCaX29/fLKqk8z+PNmzdxcXFx5/F3hVVHR0fx4cOH5+gSADfMmsZ5epq+HwAAANeeZfrezbWj8jyPoiiiKIrIsuxWVdRgMLg1ze/a4eFhXF1dlY/Ly8vn6B4AAAAAS+bJlVL9fj9ev35dLnR+rVKpRL1ej48fP956Tq1Wm/laW1tbsbW19dQuAQAsjX6/X97opdfrxefPn8sv6IqiKO9wXBTFVPX5fW0AAOvgl0Kp4XA4NV3v+Pi4bOt2u9FoNGYufF4URdRqNQMqAGBjdLvdeP/+fUREnJycxOvXr+P8/DwiIprNZvlzURSxt7cX7Xb7p20AAOvgwaFUt9uNs7OziPi+9tPu7m4ZPtVqtTg5OYksy+Li4mJqwNRut6PVasXu7m70ej2DKQBgY/T7/Tg6OipDqUajEa1Wa2rpg2t5npcVVT+232wDAFgXDw6l6vV61Ov1qaqoa9VqNarV6szn3aykajQav9hNAIDVU61W4/Pnz+X29c1eKpVKfP36NSqVytTxlUol+v1+/POf/7yz7a4xFwDAqnmWu+8BADDbzS/lvnz5EvV6PbIsu/NuxIPB4N62WcbjcYzH43J7NBr9cn8BAFJ5lrvvAQBwv+FwGJ1O56dLGdwVSN3XdnR0FDs7O+Xj5cuXT+gpAEAaQikAgARarVacnZ2VN3zJsuxW5dNgMChvFnNX2yyHh4dxdXVVPi4vL+fxKwAAPCuhFADAnJ2cnESr1Yo8z2M4HMZwOIx6vT7z2Fqtdm/bLFtbW7G9vT31AABYdtaUAqD09u3tfaen6fsB66TT6US1Wi0Dqa9fv8b+/v6tqqeiKKJWq5WVUne1AQCsC6EUAMCcFEURzWZzal+WZbG/vx8REe12O1qtVuzu7kav15tab+q+NgCAdSCUAgCYkzzPYzKZ3Nt+fHwcEdN36ftZGwDAOhBKwRqZNfUKAAAAlpFQCmDNCCcBAIBV4O57AAAAACQnlAIAAAAgOdP3ADaUaX4AAMAiqZQCAAAAIDmhFAAAAADJCaUAAAAASE4oBQAAAEByQikAAAAAkhNKAQAAAJCcUAoAAACA5IRSAAAAACQnlAIAAAAgud8W3QHg17x9u+geAAAAwK9TKQUAAABAciqlAAAAWEuzZhecnqbvBzCbSikAAAAAkhNKAQAAAJCcUAoAAACA5IRSAAAAACQnlAIAAAAgOXffAwAA+EWz7u4GwMOolAIAAAAgOaEUAAAAAMkJpQAAAABITigFAAAAQHJCKQAAAACSE0oBAAAAkJxQCgAAAIDkhFIAAAAAJCeUAgAAACA5oRQAAAAAyQmlAAAAAEhOKAUAAABAckIpAAAAAJITSgEAAACQ3G+L7gDwc2/fLroHAAAA8LxUSgEAAACQnEopAACAB1C9DvC8VEoBAAAAkJxQCgBgjvr9frx69Wrm/n6/HxERRVGUP19vn5ycRKfTiZOTkxgOh6m6CwCQjOl7AABz0ul0Is/zqcDp2sePH+PTp08REVGv16PdbpdtzWYzzs/PI+J7QLW3tzfVDgCwDoRSAABz0mg07mx79epVfPv2LSIisiwr9xdFMXVcnufR7Xbn0j8AgEUyfQ8AYEGyLJsKpCIiut1uVCqVqX2VSmVmtRUAwCpTKQUAsADD4TA6nU5ERPR6vTg4OIg8z+9cP2owGNz5WuPxOMbjcbk9Go2eta8AAPMglIIl41bDAJthf3+/rJLK8zzevHkTFxcXdx5/32LnR0dH8eHDh2fuIcB6mjXePj1N3w/A9D0AgIW4uXZUnudRFEUURRFZlt2qihoMBrem+d10eHgYV1dX5ePy8nJe3QYAeDZCKQCAxPr9frx+/frW/kqlEvV6feZzarXana+3tbUV29vbUw8AgGVn+h4AQALD4XBqut7x8XHZ1u12o9FozFz4vCiKqNVq91ZKAQCsIqEUAMCcdLvdODs7i4jv6z7t7u6W4VOtVouTk5PIsiwuLi6i3W6Xz2u329FqtWJ3dzd6vd5UGwDAungxmUwmi+7EXUajUezs7MTV1ZUydDaGhc5ZNhb+ZNOt4nhkFfsMq8A4bX0Z78DzeuhYxJpSAAAAACQnlAIAAAAguUetKdXv92Nvby/Oz8+n9hdFEZ1Op7yd8f7+frkY531tAAAAAGymB4dS18FSv9+/1dZsNsugqiiK2NvbKxfkvK8NAAAAgM304FCq0WjM3F8UxdR2nufR7XZ/2gYAAADA5nrymlLdbjcqlcrUvkqlEv1+/942AAAAADbXo9aUmmU4HM7cPxgM7m2bZTwex3g8LrdHo9FTuwcAAADAEprb3ffuCqTuazs6OoqdnZ3y8fLly/l0DgAAAICFenIolWXZrcqnwWAQWZbd2zbL4eFhXF1dlY/Ly8undg8AAACAJfTkUKper8/cX6vV7m2bZWtrK7a3t6ceAAAAAKyfX1pTajgcltVOeZ5PtRVFEbVarayUuqsNAAAAgM314FCq2+3G2dlZRHxf+2l3dzcajUZERLTb7Wi1WrG7uxu9Xi/a7Xb5vPvaAAAAANhMLyaTyWTRnbjLaDSKnZ2duLq6MpWPjfH27aJ7ANNOTxfdA1isVRyPrGKfYRUYp60v4x14Xg8di/zS9D0ANsesAbiBGwAA8FRPXugcAAAAAB5LKAUAAABAckIpAAAAAJITSgEAAACQnFAKAAAAgOSEUgAAAAAkJ5QCAAAAIDmhFAAAAADJCaUAAAAASO63RXcAAABg2bx9u+geAKw/lVIAAAAAJCeUAgAAACA50/dggZSFs05+/O/59HQx/QAAAFaDSikAAAAAkhNKAQAAAJCcUAoAAACA5IRSAAAAACQnlAIAAAAgOaEUAAAAAMkJpQAAAABI7rdFdwA2ydu3i+4BAAAALAeVUgAAAAAkJ5QCAAAAIDmhFAAAAADJCaUAAAAASE4oBQAAAEBy7r4HwFzMutvk6Wn6fgAAAMtJKAUAMEf9fj/29vbi/Px8an9RFNHpdCLP8yiKIvb39yPLsp+2AQCsC6EUAMCcXAdL/X7/Vluz2SyDqqIoYm9vL9rt9k/bAADWhVAKgEebNTUPuK3RaMzcXxTF1Hae59Htdn/aBgCwToRSAACJdbvdqFQqU/sqlUr0+/345z//eWdbtVpN2U2AjWEtTFgMoRQAQGLD4XDm/sFgcG/bXcbjcYzH43J7NBo9pXsAAEkIpWBOTG8C4LHuCqR+1nZ0dBQfPnx4/g4BAMzRXxbdAQCATZNl2a3Kp8FgEFmW3dt2l8PDw7i6uiofl5eX8+g2AMCzUikFAJBYvV6Pjx8/3tpfq9Uiz/M72+6ytbUVW1tbz9pH2CQq3JnFOlMwf0IpAIAEhsNhWe2U5/lUW1EUUavVykqpu9oAANaJUAqAZHzjyKbpdrtxdnYWEd/Xfdrd3Y1GoxEREe12O1qtVuzu7kav14t2u10+7742AIB18WIymUwW3Ym7jEaj2NnZiaurq9je3l50d+BRlIHDwwilWHarOB5ZxT7DIhm38VDGLfAwDx2LWOgcAAAAgOSEUgAAAAAkJ5QCAAAAIDmhFAAAAADJCaUAAAAASO63RXcAAAAAVsGsOzW6Ix/8OqEUAACwMWaFCgAshul7AAAAACQnlAIAAAAgOaEUAAAAAMkJpQAAAABITigFAAAAQHJCKQAAAACSE0oBAAAAkJxQCgAAAIDkhFIAAAAAJCeUAgAAACC53xbdAVg1b9/e3nd6mr4fAAAAsMpUSgEAAACQnFAKAAAAgOSEUgAAAAAkJ5QCAAAAILlnW+i83+9HRES1Wo2iKGI4HEa1Wo2IiKIootPpRJ7nURRF7O/vR5Zlz/XWAAAAAKyYZwulPn78GJ8+fYqIiHq9Hu12u2xrNptxfn4eEd8Dqr29val2WHWz7sgHAAAA3O3ZQqlXr17Ft2/fIiKmqqCKopg6Ls/z6Ha7z/W2AAAAAKygZ11TKsuyW9Pyut1uVCqVqX2VSqWc7gcAAADA5nm2SqnhcBidTiciInq9XhwcHESe5zEcDmcePxgMbu0bj8cxHo/L7dFo9FzdAwAAAGCJPFsodXPx8jzP482bN3FxcXHn8bPCqqOjo/jw4cNzdQkAAACAJfVsoVRRFOXd9q7vslcURWRZdqsqajAYzLz73uHhYfzxxx/l9mg0ipcvXz5XFwFYQrNuFHB6mr4fAABAWs+yplS/34/Xr1/f2l+pVKJer898Tq1Wu7Vva2srtre3px4AAAAArJ9nqZTK8zyOj4/L7W63G41GY+bC50VRRK1Wm1kpBQAAAMBmeJZQKsuyqNVqcXJyElmWxcXFRbTb7bK93W5Hq9WK3d3d6PV6U20AAAAAbJ5nW1OqWq2Wa0r96GYlVaPReK63BAAAAGBFPVsoBQAAsExm3UwDgOXxLAudAwAAAMBjCKUAAAAASM70PQBWwqwpGKen6fsBAAA8D5VSAAAAACSnUgqApWNhWgAAWH8qpQAAAABITigFAAAAQHKm78E9TCECYJ76/X5ERFSr1SiKIobDYVSr1YiIKIoiOp1O5HkeRVHE/v5+ZFm2wN4CADwvoRQAwIJ8/PgxPn36FBER9Xo92u122dZsNuP8/DwivgdUe3t7U+0AAKtOKAUAsCCvXr2Kb9++RURMVUEVRTF1XJ7n0e12U3YNAGDurCkFALBAWZbdmpbX7XajUqlM7atUKuV0PwCAdaBSCgBgQYbDYXQ6nYiI6PV6cXBwEHmex3A4nHn8YDBI2DsAHmLWOrSnp+n7AatIKAX/y6LmAKR2c/HyPM/jzZs3cXFxcefxd4VV4/E4xuNxuT0ajZ6zmwAAcyGUYu355gKAZVUURXm3veu77BVFEVmW3aqKGgwGd9597+joKD58+DDv7gIAPCtrSgEALEC/34/Xr1/f2l+pVKJer898Tq1Wm7n/8PAwrq6uysfl5eWz9hUAYB5USgEALECe53F8fFxud7vdaDQaMxc+L4oiarXanZVSW1tbsbW1NcfeAvAYZmvAwwilAAAWIMuyqNVqcXJyElmWxcXFRbTb7bK93W5Hq9WK3d3d6PV6U23AbdYHBVg9QikAgAWpVqvlmlI/ullJ1Wg0UnYLACAJa0oBAAAAkJxQCgAAAIDkhFIAAAAAJCeUAgAAACA5oRQAAAAAybn7HhvJLYMBAABgsVRKAQAAAJCcUAoAAACA5IRSAAAAACRnTSkAAACYs1nr2p6epu8HLBOhFAAry+AOAABWl1AKAABYKe6kDLAehFIArJWH/o+KiiqA5aMCFmCzWOgcAAAAgOSEUgAAAAAkZ/oea8X6AgAAwKowZZVNp1IKAAAAgOSEUgAAAAAkZ/oeK8tUPQAAAFhdKqUAAAAASE4oBQAAAEByQikAAAAAkrOmFAAAsLSsIwqwvoRSAAAAsCRmBbGnp+n7ASkIpQDgfxkEAgBAOkIpADbSU6aDCK8AAODphFIAAEBy1ooCQCgFAPfwP00AADAff1l0BwAAAADYPEIpAAAAAJITSgEAAACQnDWlWAnWdAEAAID1IpQCAACAJTbrS/rT0/T9gOcmlCIZf0gBAACAa0IpAABgrizFAMAsQimWjkELAAAArD933wMAAAAgOZVSLJSqKGCdPfRvnPX1AIDnYB1fVo1QirkQNgGbxt89YBP52wfAUwilAAAAYMOpsmIRhFIAsGAGgQDAvBhnsMyEUivKHxYAAABglSUJpYqiiE6nE3meR1EUsb+/H1mWpXhrAICVZPzEsrF+FCwX1yTrIEko1Ww24/z8PCK+D7D29vai3W6neGtmeMrdoPzhA4A0jJ8AgHU391CqKIqp7TzPo9vtzvtt+V9CJID19pQvGlhexk/My0OXgDCGBCKe9jfD2IOH+Mu836Db7UalUpnaV6lUot/vz/utAQBWkvETALAJ5l4pNRwOZ+4fDAa39o3H4xiPx+X21dVVRESMRqO59G1V/P3vDzvu//2/533f5349AB5u1t/gr19v7/uf/3nY6836KH3o58us95313FnHPbdFvO/1OGQymcz3jW54zPgpwhhq0zzlOpj1N8OYDzbPU677WR8ts/62+AjabA8dPy3s7nuzBltHR0fx4cOHW/tfvnyZoEcAsNx2dpb7uU95j6dI9b5//vln7Czql/xfd4VVxlAs+D9NYIMs+7iA5fKz8dPcQ6ksy259qzcYDGbePebw8DD++OOPcvvf//53DAaD+I//+I948eLFvLu60kajUbx8+TIuLy9je3t70d3ZeM7H8nAulodzsVycj4ebTCbx559/xt/+9rdk7/mY8VOEMdRTuBaWh3OxPJyL5eJ8LA/n4uEeOn6aeyhVr9fj48ePt/bXarVb+7a2tmJra2tqn1sfP8729raLY4k4H8vDuVgezsVycT4eJnWF1GPGTxHGUM/BtbA8nIvl4VwsF+djeTgXD/OQ8dPcFzrP83xquyiKqNVqBkoAAHcwfgIANkGSNaXa7Xa0Wq3Y3d2NXq8X7XY7xdsCAKws4ycAYN0lCaXyPI/j4+OIiGg0GinecuNsbW3Ff//3f98q3WcxnI/l4VwsD+diuTgfy8/4KQ3XwvJwLpaHc7FcnI/l4Vw8vxeTlPc3BgAAAIBIsKYUAAAAAPxIKAUAAABAcknWlOJpiqKITqcTeZ5HURSxv79/59137jv2oa/TarXi8PDQHX5mSHEu+v1+dLvdiIjo9Xrx+fNn5yLS/Ns/5j02nWthefiMgNlcG8vF58biGEMtD9fBcvE5sSQmLL1qtVr+fHFxMWk0Gr907ENe5/z8fBIRk2/fvj2x1+spxbk4Pj6e+vnmsZssxb/9Y95j07kWlofPCJjNtbFcfG4sjjHU8nAdLBefE8tBKLXkLi4ubv0hybLs0cc+9HXa7fYkz3MXywwpzsX5+fnUa15cXEwiYnJxcfHk/q+yFP/2j3mPTedaWB4+I2A218Zy8bmxOMZQy8N1sFx8TiwPa0otuW63G5VKZWpfpVKJfr//qGMf8jqdTsctp++R4lxUq9X4/PlzuX84HJbtmyzFv/1j3mPTuRaWh88ImM21sVx8biyOMdTycB0sF58Ty0MoteSu/5D8aDAYPOrYn73OcDg0t/UnUp2Lm3+wvnz5EvV6fePPTYp/+8e8x6ZzLSwPnxEwm2tjufjcWBxjqOXhOlguPieWh1BqRd31H/9jj71u+/r1a9Tr9ad1akM997m4ud3pdKLdbv9axzbAvP7tf/U9Np1rYXn4jIDZXBvLxefG4hhDLQ/XwXLxOZGeu+8tyKdPn+Li4uLO9jdv3pSp9o9p7WAwmJm23nfsfW3dbjf+/ve///ovs+KW6Vzc1Gq14uzsTLIe9/97PubY5zgvpDkfN7kW7uYzgk2zTJ/Zro3lOh83+dz4P8ZQy8N1sFx8TiyRRS9qxf3uWjht1gJp9x17X9vZ2dnk48eP5SMiJu/fv5+cn58/+++zylKci2vHx8flgoTfvn3b+AXxUvzbP+Y9Np1rYXn4jIDZXBvLxefG4hhDLQ/XwXLxObE8VEotuTzPp7aLooharVYmuP1+P7IsizzP7z32x8T3ZtuP5YQHBwdxcHBw6/U2XYpzEfF9IbxqtRp5nsdwOIyvX7/G/v7+vH6tlZDi3/5n54X/41pYHj4jYDbXxnLxubE4xlDLw3WwXHxOLJFFp2L83MXFxeT9+/eTdrs9ef/+/VR622g0JsfHxw869r62yeR7in58fDyJiMn+/r4Ed4Z5n4vr27befGRuqTuZTNJcBz+7Rvg/roXl4TMCZnNtLBefG4tjDLU8XAfLxefEcngxmUwmaWMwAAAAADadu+8BAAAAkJxQCgAAAIDkhFIAAAAAJCeUAgAAACA5oRQAAAAAyQmlAAAAAEhOKAUAAABAckIpAAAAAJITSgEAAACQnFAKAAAAgOSEUgAAAAAkJ5QCAAAAILn/D72Be+R38LDBAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1200x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Compute tune derivatives and estimate spread from linear surrogate model using MC\n",
    "\n",
    "# Compute derivative\n",
    "\n",
    "kn = torch.zeros(nq, dtype=torch.float64)\n",
    "dnu_dk = torch.func.jacrev(lambda kn: tune(ring, [kn], ('kn', ['Quadrupole'], None, None), matched=True, limit=1, epsilon=None))(kn)\n",
    "\n",
    "# Sample\n",
    "\n",
    "kns = 5.0E-3*torch.randn((8192, nq), dtype=torch.float64)\n",
    "nuqx, nuqy = dnu_dk @ kns.T\n",
    "\n",
    "# Spread\n",
    "\n",
    "print(nuqx.std(-1))\n",
    "print(nuqy.std(-1))\n",
    "print()\n",
    "\n",
    "# Plot samples\n",
    "\n",
    "fig, (ax, ay) = plt.subplots(1, 2, figsize=(12, 5))\n",
    "ax.hist((nuqx - nuqx.mean()).cpu().numpy(), bins=100, range=(-5.0E-3, +5.0E-3), color='blue', alpha=0.7)\n",
    "ay.hist((nuqy - nuqy.mean()).cpu().numpy(), bins=100, range=(-5.0E-3, +5.0E-3), color='blue', alpha=0.7)\n",
    "plt.tight_layout() \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0249c95c-97bf-4c38-a034-9bddd5eb2e6e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([0.0012, 0.0008], dtype=torch.float64)\n"
     ]
    }
   ],
   "source": [
    "# Compute spread using error propagation\n",
    "\n",
    "print((dnu_dk @ (5.0E-3*torch.eye(nq,  dtype=torch.float64))**2 @ dnu_dk.T).diag().sqrt())"
   ]
  }
 ],
 "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
}