{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fb6264f7-57eb-4dd9-8799-0e9cded684af",
   "metadata": {},
   "source": [
    "## ELETTRA-40: ID linear optics distortion (transparent ID)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f9f12b5c-0667-43e7-8c8e-f88b4d677b82",
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this example the matched ID is inserted into the lattice\n",
    "# Given an ID liner transport matrix parameterization exp(S A) with A diagonal and diag(A) = (a, b, c, d)\n",
    "# The transparency condition (at a location where the alpha is zero) is bx = sqrt(b/a) and by = sqrt(d/c)\n",
    "# However, for realistic ID parameters this leads to very small beta-function values (< 1 m) at the ID, which are not feasible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e71ae06d-6ad2-43b6-9203-c888bdd455be",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import\n",
    "\n",
    "import torch\n",
    "from torch import Tensor\n",
    "\n",
    "from pathlib import Path\n",
    "\n",
    "import matplotlib\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib.patches import Rectangle\n",
    "matplotlib.rcParams['text.usetex'] = True\n",
    "\n",
    "from model.library.element import Element\n",
    "from model.library.line import Line\n",
    "from model.library.quadrupole import Quadrupole\n",
    "from model.library.matrix import Matrix\n",
    "from model.library.drift import Drift\n",
    "\n",
    "from model.command.external import load_lattice\n",
    "from model.command.build import build\n",
    "from model.command.tune import tune\n",
    "from model.command.orbit import dispersion\n",
    "from model.command.twiss import twiss\n",
    "from model.command.advance import advance\n",
    "from model.command.coupling import coupling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c641f373-4c51-458f-b4cc-71541b50aca3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set data type and device\n",
    "\n",
    "Element.dtype = dtype = torch.float64\n",
    "Element.device = device = torch.device('cpu')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f8e2c68a-fdb6-4676-b0a2-ce70826c8890",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load lattice (ELEGANT table)\n",
    "# Note, lattice is allowed to have repeated elements\n",
    "\n",
    "path = Path('elettra.lte')\n",
    "data = load_lattice(path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "02625ecc-4a7d-405f-aac5-77cfc5f71aee",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Marker': 12, 'Drift': 708, 'BPM': 168, 'Quadrupole': 360, 'Dipole': 156}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Build and setup lattice\n",
    "\n",
    "ring:Line = build('RING', 'ELEGANT', data)\n",
    "\n",
    "# Flatten sublines\n",
    "\n",
    "ring.flatten()\n",
    "\n",
    "# Remove all marker elements but the ones starting with MLL (long straight section centers)\n",
    "\n",
    "ring.remove_group(pattern=r'^(?!MLL_).*', kinds=['Marker'])\n",
    "\n",
    "# Replace all sextupoles with quadrupoles\n",
    "\n",
    "def factory(element:Element) -> None:\n",
    "    table = element.serialize\n",
    "    table.pop('ms', None)\n",
    "    return Quadrupole(**table)\n",
    "\n",
    "ring.replace_group(pattern=r'', factory=factory, kinds=['Sextupole'])\n",
    "\n",
    "# Set linear dipoles\n",
    "\n",
    "def apply(element:Element) -> None:\n",
    "    element.linear = True\n",
    "\n",
    "ring.apply(apply, kinds=['Dipole'])\n",
    "\n",
    "# Merge drifts\n",
    "\n",
    "ring.merge()\n",
    "\n",
    "# Change lattice start\n",
    "\n",
    "ring.start = \"MLL_S01\"\n",
    "\n",
    "# Describe\n",
    "\n",
    "ring.describe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "99db8165-52ba-4db9-afb8-dc052102c0b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute tunes (fractional part)\n",
    "\n",
    "nux, nuy = tune(ring, [], matched=True, limit=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c328bf74-8419-4886-8260-0df69ce51b21",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute twiss parameters\n",
    "\n",
    "ax, bx, ay, by = twiss(ring, [], matched=True, advance=True, full=False).T\n",
    "\n",
    "axi, *_ = ax\n",
    "bxi, *_ = bx\n",
    "ayi, *_ = ay\n",
    "byi, *_ = by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d1082a39-ebb6-49a8-8184-740f414ebe01",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute phase advances\n",
    "\n",
    "mux, muy = advance(ring, [], alignment=False, matched=True).T\n",
    "\n",
    "mux = mux.cumsum(-1)\n",
    "muy = muy.cumsum(-1)\n",
    "\n",
    "mux = (mux % mux.max()).roll(1)\n",
    "muy = (muy % muy.max()).roll(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4bafba70-94ea-4806-98ac-24d71d43628e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define rescaled ID\n",
    "\n",
    "ca, cb, cc, cd = -0.034441907232402175, -0.04458009513208418, 0.056279356423643276, 0.08037110220505986\n",
    "\n",
    "fx = bxi**2*ca/cb\n",
    "fy = byi**2*cc/cd\n",
    "ca, cb = ca/fx**0.5, cb*fx**0.5\n",
    "cc, cd = cc/fy**0.5, cd*fy**0.5\n",
    "\n",
    "A = torch.tensor([[ca, 0.0, 0.0, 0.0], [0.0, cb, 0.0, 0.0], [0.0, 0.0, cc, 0.0], [0.0, 0.0, 0.0, cd]], dtype=dtype)\n",
    "\n",
    "mask = torch.triu(torch.ones_like(A, dtype=torch.bool))\n",
    "\n",
    "Ax, Bx, _, _, Cx, _, _, Ay, By, Cy = A[mask].to(torch.complex128)\n",
    "\n",
    "ID = Matrix('ID', length=0.0, A=(A[mask]).tolist())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2e15bd15-0f0a-48e5-a323-494ca7fff06f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Matrix': 1,\n",
       " 'Drift': 708,\n",
       " 'BPM': 168,\n",
       " 'Quadrupole': 360,\n",
       " 'Dipole': 156,\n",
       " 'Marker': 11}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Replace marker by ID\n",
    "\n",
    "ring.replace('MLL_S01', ID)\n",
    "\n",
    "# Describe\n",
    "\n",
    "ring.describe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "85cb7bab-30c5-4a93-952f-55a4a8de5ea3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute tunes (fractional part)\n",
    "\n",
    "nux_id, nuy_id = tune(ring, [], matched=True, limit=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3d573930-474e-484f-9828-020944d99154",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute twiss parameters\n",
    "\n",
    "ax_id, bx_id, ay_id, by_id = twiss(ring, [], matched=True, advance=True, full=False).T\n",
    "\n",
    "axf, *_ = ax_id\n",
    "bxf, *_ = bx_id\n",
    "ayf, *_ = ay_id\n",
    "byf, *_ = by_id"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7234f798-f858-4803-a415-8718871372e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute phase advances\n",
    "\n",
    "mux_id, muy_id = advance(ring, [], alignment=False, matched=True).T\n",
    "\n",
    "mux_id = mux_id.cumsum(-1)\n",
    "muy_id = muy_id.cumsum(-1)\n",
    "\n",
    "mux_id = (mux_id % mux_id.max()).roll(1)\n",
    "muy_id = (muy_id % muy_id.max()).roll(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "e0af4b32-fd59-4e44-83e2-baff162df977",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(0.0062, dtype=torch.float64)\n",
      "tensor(-0.0107, dtype=torch.float64)\n"
     ]
    }
   ],
   "source": [
    "# Tune shifts\n",
    "\n",
    "print(nux - nux_id)\n",
    "print(nuy - nuy_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "c805d953-2fd9-498f-8a7e-18713de447ba",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(0.0062, dtype=torch.float64)\n",
      "tensor(0.0062, dtype=torch.float64)\n",
      "\n",
      "tensor(-0.0107, dtype=torch.float64)\n",
      "tensor(-0.0107, dtype=torch.float64)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Tune shifts (exact)\n",
    "\n",
    "def dnux(Ax, Bx, Cx, ax, bx, nux):\n",
    "    return nux - torch.arccos(torch.cos(2*nux*torch.pi)*torch.cosh(torch.sqrt(Bx**2 - Ax*Cx)) - ((Ax*bx**2 - 2*ax*bx*Bx + Cx + ax**2*Cx)*torch.sin(2*nux*torch.pi)*torch.sinh(torch.sqrt(Bx**2 - Ax*Cx)))/(2*bx*torch.sqrt(Bx**2 - Ax*Cx)))/(2*torch.pi)\n",
    "\n",
    "def dnuy(Ay, By, Cy, ay, by, nuy):\n",
    "    return nuy - torch.arccos(torch.cos(2*nuy*torch.pi)*torch.cosh(torch.sqrt(By**2 - Ay*Cy)) - ((Ay*by**2 - 2*ay*by*By + Cy + ay**2*Cy)*torch.sin(2*nuy*torch.pi)*torch.sinh(torch.sqrt(By**2 - Ay*Cy)))/(2*by*torch.sqrt(By**2 - Ay*Cy)))/(2*torch.pi)\n",
    "\n",
    "print((nux - nux_id))\n",
    "print(dnux(Ax, Bx, Cx, axi, bxi, nux).real)\n",
    "print()\n",
    "\n",
    "print((nuy - nuy_id))\n",
    "print(dnux(Ay, By, Cy, ayi, byi, nuy).real)\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "9081dd26-c671-4333-943b-b28b5d570c3e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(0.0062, dtype=torch.float64)\n",
      "tensor(0.0062, dtype=torch.float64)\n",
      "\n",
      "tensor(-0.0107, dtype=torch.float64)\n",
      "tensor(-0.0107, dtype=torch.float64)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Tune shifts (approximate)\n",
    "\n",
    "def dnux(Ax, Bx, Cx, ax, bx, nux):\n",
    "    cx = (1 + ax**2)/bx\n",
    "    return - ((Ax*bx - 2*ax*Bx + cx*Cx)/(4*torch.pi) - ((Ax**2*bx**2 - 4*ax*Ax*bx*Bx + 4*bx*Bx**2*cx + 2*Ax*(-2 + bx*cx)*Cx + cx*Cx*(-4*ax*Bx + cx*Cx))*torch.cos(2*nux*torch.pi)/torch.sin(2*nux*torch.pi))/(16*torch.pi))\n",
    "\n",
    "def dnuy(Ay, By, Cy, ay, by, nuy):\n",
    "    cy = (1 + ay**2)/by\n",
    "    return - ((Ay*by - 2*ay*By + cy*Cy)/(4*torch.pi) - ((Ay**2*by**2 - 4*ay*Ay*by*By + 4*by*By**2*cy + 2*Ay*(-2 + by*cy)*Cy + cy*Cy*(-4*ay*By + cy*Cy))*torch.cos(2*nuy*torch.pi)/torch.sin(2*nuy*torch.pi))/(16*torch.pi))\n",
    "\n",
    "print((nux - nux_id))\n",
    "print(dnux(Ax, Bx, Cx, axi, bxi, nux).real)\n",
    "print()\n",
    "\n",
    "print((nuy - nuy_id))\n",
    "print(dnux(Ay, By, Cy, ayi, byi, nuy).real)\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "1bfb7730-0af7-4ed4-9819-15c7a3656257",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([-7.2024e-15,  9.4152e+00], dtype=torch.float64)\n",
      "tensor([-7.1968e-15,  9.4152e+00], dtype=torch.float64)\n",
      "\n",
      "tensor([-2.8043e-15,  1.6398e+00], dtype=torch.float64)\n",
      "tensor([-2.8375e-15,  1.6398e+00], dtype=torch.float64)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Twiss at the observation point (exact)\n",
    "\n",
    "def csx(Ax, Bx, Cx, ax, bx, nux):\n",
    "    cx = (1 + ax**2)/bx\n",
    "    hax = (2*ax*torch.cosh(torch.sqrt(Bx**2 - Ax*Cx))*torch.sin(2*nux*torch.pi) + ((2*bx*Bx*torch.cos(2*nux*torch.pi) + (-(Ax*bx**2) + Cx + ax**2*Cx)*torch.sin(2*nux*torch.pi))*torch.sinh(torch.sqrt(Bx**2 - Ax*Cx)))/(bx*torch.sqrt(Bx**2 - Ax*Cx)))/(2*torch.sqrt(1 - (torch.cos(2*nux*torch.pi)*torch.cosh(torch.sqrt(Bx**2 - Ax*Cx)) - ((Ax*bx**2 - 2*ax*bx*Bx + Cx + ax**2*Cx)*torch.sin(2*nux*torch.pi)*torch.sinh(torch.sqrt(Bx**2 - Ax*Cx)))/(2*bx*torch.sqrt(Bx**2 - Ax*Cx)))**2))    \n",
    "    hbx = (bx*torch.cosh(torch.sqrt(Bx**2 - Ax*Cx))*torch.sin(2*nux*torch.pi) + ((Cx*torch.cos(2*nux*torch.pi) + (-(bx*Bx) + ax*Cx)*torch.sin(2*nux*torch.pi))*torch.sinh(torch.sqrt(Bx**2 - Ax*Cx)))/torch.sqrt(Bx**2 - Ax*Cx))/torch.sqrt(1 - (torch.cos(2*nux*torch.pi)*torch.cosh(torch.sqrt(Bx**2 - Ax*Cx)) - ((Ax*bx**2 - 2*ax*bx*Bx + Cx + ax**2*Cx)*torch.sin(2*nux*torch.pi)*torch.sinh(torch.sqrt(Bx**2 - Ax*Cx)))/(2*bx*torch.sqrt(Bx**2 - Ax*Cx)))**2)\n",
    "    return torch.stack([hax, hbx])\n",
    "\n",
    "def csy(Ay, By, Cy, ay, by, nuy):\n",
    "    cy = (1 + ay**2)/by\n",
    "    hay = (2*ay*torch.cosh(torch.sqrt(By**2 - Ay*Cy))*torch.sin(2*nuy*torch.pi) + ((2*by*By*torch.cos(2*nuy*torch.pi) + (-(Ay*by**2) + Cy + ay**2*Cy)*torch.sin(2*nuy*torch.pi))*torch.sinh(torch.sqrt(By**2 - Ay*Cy)))/(by*torch.sqrt(By**2 - Ay*Cy)))/(2*torch.sqrt(1 - (torch.cos(2*nuy*torch.pi)*torch.cosh(torch.sqrt(By**2 - Ay*Cy)) - ((Ay*by**2 - 2*ay*by*By + Cy + ay**2*Cy)*torch.sin(2*nuy*torch.pi)*torch.sinh(torch.sqrt(By**2 - Ay*Cy)))/(2*by*torch.sqrt(By**2 - Ay*Cy)))**2))    \n",
    "    hby = (by*torch.cosh(torch.sqrt(By**2 - Ay*Cy))*torch.sin(2*nuy*torch.pi) + ((Cy*torch.cos(2*nuy*torch.pi) + (-(by*By) + ay*Cy)*torch.sin(2*nuy*torch.pi))*torch.sinh(torch.sqrt(By**2 - Ay*Cy)))/torch.sqrt(By**2 - Ay*Cy))/torch.sqrt(1 - (torch.cos(2*nuy*torch.pi)*torch.cosh(torch.sqrt(By**2 - Ay*Cy)) - ((Ay*by**2 - 2*ay*by*By + Cy + ay**2*Cy)*torch.sin(2*nuy*torch.pi)*torch.sinh(torch.sqrt(By**2 - Ay*Cy)))/(2*by*torch.sqrt(By**2 - Ay*Cy)))**2)\n",
    "    return torch.stack([hay, hby])\n",
    "\n",
    "print(torch.stack([axf, bxf]))\n",
    "print(csx(Ax, Bx, Cx, axi, bxi, nux).real)\n",
    "print()\n",
    "\n",
    "print(torch.stack([ayf, byf]))\n",
    "print(csy(Ay, By, Cy, ayi, byi, nuy).real)\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "5aac5c72-b235-4b2f-a129-b7c50fa02f3a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor([-7.2024e-15,  9.4152e+00], dtype=torch.float64)\n",
      "tensor([-7.2003e-15,  9.4152e+00], dtype=torch.float64)\n",
      "\n",
      "tensor([-2.8043e-15,  1.6398e+00], dtype=torch.float64)\n",
      "tensor([-2.8364e-15,  1.6398e+00], dtype=torch.float64)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Twiss at the observation point (approximate)\n",
    "\n",
    "def csx(Ax, Bx, Cx, ax, bx, nux):\n",
    "    cx = (1 + ax**2)/bx\n",
    "    hax = ax + (-(Ax*bx) + cx*Cx - (ax*Ax*bx - 2*bx*Bx*cx + ax*cx*Cx)*(torch.cos(2*nux*torch.pi)/torch.sin(2*nux*torch.pi)))/2.0 + ((1.0/torch.sin(2*nux*torch.pi)**3)*(16*ax*(Bx**2 - Ax*Cx)*torch.sin(2*nux*torch.pi)**3 + 4*ax*torch.sin(2*nux*torch.pi)*(4*(Bx**2 - Ax*Cx)*torch.cos(2*nux*torch.pi)**2 + (Ax*bx - 2*ax*Bx + cx*Cx)**2 * torch.sin(2*nux*torch.pi)**2) + 4*(Ax*bx - 2*ax*Bx + cx*Cx)*(-2*Bx*torch.cos(2*nux*torch.pi) + (Ax*bx - cx*Cx)*torch.sin(2*nux*torch.pi))*torch.sin(4*nux*torch.pi) + 3*ax*(Ax*bx - 2*ax*Bx + cx*Cx)**2*(1.0/torch.sin(2*nux*torch.pi)) * torch.sin(4*nux*torch.pi)**2))/32.0\n",
    "    hbx = bx - bx*Bx + ax*Cx + (-0.5*(Ax*bx**2) + ax*bx*Bx + Cx - 0.5*(bx*cx*Cx)) * (torch.cos(2*nux*torch.pi)/torch.sin(2*nux*torch.pi)) + (4*bx*(Bx**2 - Ax*Cx) + bx*(Ax*bx - 2*ax*Bx + cx*Cx)**2 + 4*bx*(Bx**2 - Ax*Cx)*(torch.cos(2*nux*torch.pi)/torch.sin(2*nux*torch.pi))**2 + 3*bx*(Ax*bx - 2*ax*Bx + cx*Cx)**2*(torch.cos(2*nux*torch.pi)/torch.sin(2*nux*torch.pi))**2 + 2*(Ax*bx - 2*ax*Bx + cx*Cx)*(bx*Bx - ax*Cx - Cx*(torch.cos(2*nux*torch.pi)/torch.sin(2*nux*torch.pi)))*(1.0/torch.sin(2*nux*torch.pi)**2) * torch.sin(4*nux*torch.pi))/8.0\n",
    "    return torch.stack([hax, hbx])\n",
    "\n",
    "def csy(Ay, By, Cy, ay, by, nuy):\n",
    "    cy = (1 + ay**2)/by\n",
    "    hay = ay + (-(Ay*by) + cy*Cy - (ay*Ay*by - 2*by*By*cy + ay*cy*Cy)*(torch.cos(2*nuy*torch.pi)/torch.sin(2*nuy*torch.pi)))/2.0 + ((1.0/torch.sin(2*nuy*torch.pi)**3)*(16*ay*(By**2 - Ay*Cy)*torch.sin(2*nuy*torch.pi)**3 + 4*ay*torch.sin(2*nuy*torch.pi)*(4*(By**2 - Ay*Cy)*torch.cos(2*nuy*torch.pi)**2 + (Ay*by - 2*ay*By + cy*Cy)**2 * torch.sin(2*nuy*torch.pi)**2) + 4*(Ay*by - 2*ay*By + cy*Cy)*(-2*By*torch.cos(2*nuy*torch.pi) + (Ay*by - cy*Cy)*torch.sin(2*nuy*torch.pi))*torch.sin(4*nuy*torch.pi) + 3*ay*(Ay*by - 2*ay*By + cy*Cy)**2*(1.0/torch.sin(2*nuy*torch.pi)) * torch.sin(4*nuy*torch.pi)**2))/32.0\n",
    "    hby = by - by*By + ay*Cy + (-0.5*(Ay*by**2) + ay*by*By + Cy - 0.5*(by*cy*Cy)) * (torch.cos(2*nuy*torch.pi)/torch.sin(2*nuy*torch.pi)) + (4*by*(By**2 - Ay*Cy) + by*(Ay*by - 2*ay*By + cy*Cy)**2 + 4*by*(By**2 - Ay*Cy)*(torch.cos(2*nuy*torch.pi)/torch.sin(2*nuy*torch.pi))**2 + 3*by*(Ay*by - 2*ay*By + cy*Cy)**2*(torch.cos(2*nuy*torch.pi)/torch.sin(2*nuy*torch.pi))**2 + 2*(Ay*by - 2*ay*By + cy*Cy)*(by*By - ay*Cy - Cy*(torch.cos(2*nuy*torch.pi)/torch.sin(2*nuy*torch.pi)))*(1.0/torch.sin(2*nuy*torch.pi)**2) * torch.sin(4*nuy*torch.pi))/8.0\n",
    "    return torch.stack([hay, hby])\n",
    "\n",
    "print(torch.stack([axf, bxf]))\n",
    "print(csx(Ax, Bx, Cx, axi, bxi, nux).real)\n",
    "print()\n",
    "\n",
    "print(torch.stack([ayf, byf]))\n",
    "print(csy(Ay, By, Cy, ayi, byi, nuy).real)\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "803b534a-1078-4754-bf43-7a30aeba4ac8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABKYAAAGGCAYAAABBiol3AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAOcVJREFUeJzt3V9sXFl+H/hfieru6X9kkWpj/uyMMyp2Z/y0sSgxgNdjwJgmbWRjYIERSzLgzSRet1g9BvbNQzaDfVxYzep52wc3S404O44fRNYMkIckSLN6HMRjBxhJVM8uEMT28KphYz02YlKX7D/zxy3VPtRUiSVWkSyK5KWKnw8gSHX/nnvrHN7Lr849N1ev1+sBAAAAAEfsVNYFAAAAAOBkEkwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkInTWRfgsKRpGouLi7G0tBTLy8t7WqdcLkc+n2+tPzMz09N8AAAAAPauL3tMraysxOLiYqRpGuvr63tap1wuR0TE9PR0TE9Px9jYWJRKpT3PBwAAAKA3uXq9Xs+6EIelWq3G1atX49atW7suOzw8HHfu3Gn1iIqIyOVy0Tw9u80HAAAAoDd92WOqV0mSRJqmbaFTU61W23U+AAAAAL3r2zGmepEkScfp+Xw+0jTddX4n9+/fj7/+67+O559/PnK53EEVFQAAAOBYq9fr8f7778dnPvOZOHVq5z5RgqkdjIyMxPr6eseeUlvnd/LXf/3X8bnPfe4QSwcAAABwfP3VX/1VfPazn91xGcHUDnYbOH2n+c8//3xENL6EwcHBAy0XAAAAwHG1ubkZn/vc51rZyE4EUxFRKBQ6Tk/TNAqFwq7zO2k+vjc4OCiYAgAAAE6cvQxtZPDzaART+Xy+41hSExMTu84HAAAAoHd9HUx1e9QuSZIol8tt0+bm5tresFetVmN6enrP8wEAAADoTa5er9ezLsRBS5IkqtVqXL9+PVZWVmJmZibGx8djamoqIiIqlUrMz8/H6upq23rlcrn1aN6NGzdifn6+p/lbbW5uxtDQUGxsbHiUDwAAADgxeslE+jKYOg4EUwAAAMBJ1Esm0teP8gEAAABwfAmmAAAAAMiEYAoAAACATAimAAAAAMiEYAoAAACATAimAAAAAMiEYAoAAACATAimAAAAAMiEYAoAAACATAimAAAAAMiEYArgCFQqlSgWi1kXAwAA4FgRTAEcgVu3bkW1Ws26GAAAAMeKYAp+qlqtRpqmWReDPrWwsBD1ev1I96lOAwAAx51gCn6qWCzGzZs3sy4GHBh1GgAAOO4EU3CS/et/HfGNb7RP+8Y3GtMBAADgkAmmjrF6PeJHPzr+fx716aRqtRrnz5+PXC4X58+fj1qtFhERs7OzkcvlIkmStuUnJydjdHS09TlJkigWizE8PNzaxsrKyrb9rKysxOTkZAwPD8fw8HBMTk7GyspKFIvFyOVyrW3ncrnW57536lTE7//+g3DqG99ofD6V/Y+GvXyvlUql9T02v9vR0dFWHepluXK5HJOTkxERUSqVYnh4uPUY3MPrlcvl1np7raflcjmGh4fbPjeP6eH6n6Zp69hHR0ejUqn0dG72Uqe7tTsAAICjdDrrAtDdj38c8U/+Sdal2N1/+A8Rn/jE/tatVCpRKpViYWEhLly4ELVaLSYnJ2N1dTXm5+dbgcDq6mpr+Vqt1voc0fgFe2RkJJaWlqJQKMTs7Gy8/PLLcefOncjn8xERrV/+p6amYmlpKSIilpeX4/r163Ht2rWYm5uL8+fPx9LSUkxMTDzS+XisfOUrjb9///cj/uAPIj7+OOI3f/PB9Azt5XtN0zRqtVpcuXIl5ubmolQqxezsbKvOFAqFnpZLkiTOnz8faZrG3Nxc5PP5Vp2cmZmJ+fn5SJIkZmdn48aNG7G0tLTnevqwtbW1WFlZiStXrsT8/HxENAKxYrEYhUKhVc75+fkolUoxMTHRKudu52a3Or1Tu2vuAwAA4CgIpshMmqatX46np6cjImJsbCzW1tZifn4+FhYWYnl5OUZHR6NYLLZ+UW/+Mt40MzPTtt1r167F8PBwLC4utrZbLBbbQqmIaPtlvbm9fD7fCj1OjK985UEodfr0sQilIvb2vTbdunWr9e+JiYkYHh5u1aFelkuSJCYmJtqWK5VKrVAqolFHx8bGYnR0NKrVakxNTe2pnnYzPz/fqouzs7OtEGpqaioiIkZGRlo9oprb2+3c5PP5rnV6L+0OAADgqAimjrGnnmr0Rjrunnpqf+s1B2UulUpRKpXa5o2NjbX+3fylv1qtxszMTOsX9m6av4Q3e6skSRJJkviFu5tvfONBKPXxx43PxySc2urh73Wn5aampnZ9NK3bclvrycrKSiRJsq1+FgqFGBsbi+vXr7fqY6/1tOnChQtt242I1iOFEQ/awvr6+o7HErH7uYnYe7sDAAA4CtkPJENXuVzjEbnj/udRh2O6e/du1Ov1tj9be6wUCoVWz465ubmO26jValEsFmN0dLRtHJ+IaI394xGlDppjSv3mb0YsLzf+3jrmVMZ2+l53UigUdgxydlpuaz1phjgjIyMd1906rtRe6mknnXroddrfw/Z7bpp2a3cAAABHQTBFZraO67OTWq0Wi4uLMTU1FS+//PK2+cViMYrFYkxOTsby8nLcvXt3X/s5ke7fbx9T6itfaXy+fz/bcsXu3+tOkiTZUxD58HIPh0Q71Z2H192tnnbS7bHR3R4nfZRzoz0AAADHiWCKzDQfh7p69eq2ec23oTXfPnbt2rXW+FBbHz9K0zSq1Wpcu3YtpqenO4YRhUIhCoVCx0f5mvvZOpj2ifIv/sX2x/a+8pXG9Azt5XvtJkmSqFaruw5iv5flJiYmIp/Pb6s7KysrsbKyEpcvX25ta6d6epD2em661em9tDsAAICjYowpMnXt2rU4f/58a7yb9fX1WFpaav1dLBbbBoJeWlqK0dHROH/+fGuQ53w+H1evXo18Ph8jIyMdf+FeWFiIycnJ1lvP0jSN5eXluHnzZuvxpXw+H9evX2+FWPPz8ydvIPRjYq/fa9Pk5GTMzs5GmqZx5cqVyOfzHR+n2+tyW127di2KxWJENHoqNd/Kt7Ve7lZPD1Iv56Zbnd6t3QEAABwVPabI1NjYWKyurkaSJPHyyy+3AoBr165FqVSKJEni2rVrreWbv2CXSqVYWVlpLdvssXLlypWYnJyMiYmJGB0dba3XfNPa1uUe/iV8bm4uqtVqqwxkay/fa1OpVIrZ2dkoFotx4cKFuHXrVsdQca/LbTU1NdWqO5OTkzE/Px9zc3OxvLzc2uZe6ulB2uu56Vand2p3AAAARylXr9frWReiH21ubsbQ0FBsbGzE4OBg1sWBvlQul2N2djZ2+zG21+UAAAB4dL1kInpMAQAAAJAJwRQAAAAAmejrwc/L5XLbm6lmZmZ2XL5YLMbly5ejUCh0fG18rVZrDaJdKBRieXk5xsfHWwMeAwAAALB3fTvGVLlcjohohVG1Wi2Wlpa2vfZ9q9HR0UiSZNv0qampWFpaimq1GleuXIk0TaNQKMTs7GzXN24ZYwoAAAA4iXrJRPo2mBoeHo47d+609XzK5XI7Dn5cLpe39aqqVCqt8KlarcbExMSub/GKEEwBAAAAJ9OJH/w8SZJI07RjgFSr1bqu9/AjebVaLS5cuHDQxQMAAAAg+nSMqU6P40VE5PP5SNO063qFQqFtG0mSxMTERNsyi4uLMTIyEuvr67G6uhrz8/MHUmYAAACAk6Yvg6lumoHSXszPz28bj2psbCwiHgRYlUolisViLC0tHWxBAQAAAE6AvnyUr5u9hlIrKysdpxcKhbZeVZcuXYpqtbpjLywAAAAAOuvLYGpreLRV8216u1lYWIjR0dFt06vVatvn5hhW3R4dBAAAAKC7vg2m8vl8x8Do4TGjOqnVatsGTk/TNIrFYts2mz2l9hJ2AQAAANCuL4OpiIi5ubm2N/BVq9WYnp5ufU6SJMrlcsd1kyTZFjbl8/mYmZlpm16pVGJqaqrj2/8AAAAA2FmuXq/Xsy7EYSmXy60g6caNG21v0KtUKjE/Px+rq6vb1hsdHY2lpaXWYOdNaZpGpVJpfV5bW+v6Vr7Nzc0YGhqKjY2NGBwcPIjDAQAAADj2eslE+jqYypJgCgAAADiJeslE+vZRPgAAAACON8EUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJk4nXUBDlO5XI58Ph8REWmaxszMzI7L12q1WFhYiMnJySgUCrG8vBzj4+MxNTW1720CAAAA0Fnf9pgql8sRETE9PR3T09MxNjYWpVJpx3XSNI1arRalUilKpVKMjo5uC6V63SYAAAAAneXq9Xo960IchuHh4bhz506rd1NERC6Xi50Ot1qtxsTERNs6+93m5uZmDA0NxcbGRgwODu73MAAAAAAeK71kIn3ZYypJkkjTtGPAVKvVjs02AQAAAE6yvhxjKkmSjtPz+XykabrjuouLizEyMhLr6+uxuroa8/Pzj7xNAAAAALbry2Cqm2bg1M3Y2FhERBQKhYiIqFQqUSwWY2lpad/bBAAAAKCzExVM7RYgNQOppkuXLkWpVNqxR5RQCgAAAGB/+nKMqYcDpqY0TbvOi2gMfr5VczypJEn2vU0AAAAAOuvbYCqfz3ccF2piYqLjOmmaRrFYbFun2VOqUCjsa5sAAAAAdNeXwVRExNzcXNvb8qrVakxPT7c+J0kS5XK59Tmfz8fMzExb76dKpRJTU1OtnlO7bRMAAACAvcvV6/V61oU4LOVyuRU03bhxo/WGvYhG6DQ/Px+rq6utaWmaRqVSaX1eW1trW2e3bW61ubkZQ0NDsbGxEYODgwd2TAAAAADHWS+ZSF8HU1kSTAEAAAAnUS+ZSN8+ygcAAADA8SaYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMnE66wIcpnK5HPl8PiIi0jSNmZmZPa0TEbG6uhoREQsLC615tVotFhYWYnJyMgqFQiwvL8f4+HhMTU0dfOEBAAAA+lzfBlPNgGl6ejoiGqFSqVRqC5oeNjs7G/Pz863PpVIpJicnY3l5OSIa4VatVotqtRqFQiFmZ2eFUgAAAAD7lKvX6/WsC3EYhoeH486dO60eUxERuVwuuh1umqZRLBZjaWmptc7KykqcP38+VldXo1AoRLVajYmJibZtdrO5uRlDQ0OxsbERg4ODB3BEAAAAAMdfL5lIX44xlSRJpGnaMUCq1Wpd17t582YkSdL6XCgUIqIRWgEAAABwsPryUb6t4dJW+Xy+a8iUz+fj7t27bdOaIVYzoIqIWFxcjJGRkVhfX4/V1dW2R/8AAAAA2Lu+DKa6aQZKe3X16tVYWFho9bwaGxuLiAdBVaVSaT3+BwAAAEBv+vJRvm56CaVmZ2fj8uXLrcHTIxqB1NbeU5cuXYpqtepRPwAAAIB9OJJg6o033oiXXnop3n333da0rf8+aFvDo63SNO06b6tqtRqjo6MxMzOzbfpWzZ5U3R4dBAAAAKC7Iwmm8vl8vP766zE0NNSaVq/XY25uLjY3Nw98f4VCIfL5fMfAaGJiYsd1m+NKNXtKpWnaGky9WCy2bbPZU2ovYRcAAAAA7Y4kmBoZGYnz58/H2bNnW9POnTsXV69ejcXFxUPZ59zcXNsb+KrVattjeUmSRLlcbltnZWUlVlZWYmxsLJIkiSRJolKpxMjISOTz+ZiZmWkLoSqVSkxNTXV8+x8AAAAAO8vV6/X6UezoV37lV+Lu3btRKpXi0qVLMTg4GBERX//61+N3fud3DmWf5XK5FSTduHGj7Q16lUol5ufnY3V1NSIavZ/Onj3bcbyo5ilK0zQqlUpr+traWte38m1ubsbQ0FBsbGy0jhUAAACg3/WSiRxJMPXVr3619Ua75eXl1hhOERGlUunQgqksCaYAAACAk6iXTOT0URRoYmIiLl68GBERV65ciYjGWE63b99ue7wOAAAAgJPjSMaYitj+Fr6JiYn42te+1jYOFAAAAAAnx5EEUxcvXoy1tbX4+te/3pp2+/btGB8f7/jmPAAAAAD635ENft7JtWvX4sKFC3Hu3LmsinBojDEFAAAAnES9ZCL76jF16dKlGBgYiIGBgfjt3/7tfRUyojHeVD+GUgAAAADsrudgam5uLsbHx+PmzZtx/fr1+O53vxvj4+Nt81966aUYGBiIM2fOxK//+q/H9773vQMtNAAAAACPv54f5Xvttdfi9ddfb5tWLpfj1q1bkSRJDA8Px+TkZKytrUWSJFGr1WJjYyMqlUr81m/91oEW/jjzKB8AAABwEvWSiZzudeMvvPDCtmkzMzNx4cKFePXVV+OVV17ZNr9arcb09HRExIkKpwAAAADorudH+dbW1rrOGxkZ6Th9amoqkiSJ3/u934vNzc1edwkAAABAH+o5mCoUCvHee+91nJ4kSdf18vl8LC0txdWrV3vdJQAAAAB9qOdg6sqVK/Hmm29u6/mUpmnk8/kd1z179mzXXlUAAAAAnCw9jzEVEfH666/Hq6++GpcuXYovfelLERHx9ttv72ndXC63n10CAAAA0Gd67jHV9Oabb8bdu3fjq1/9anzrW9/a83o7jVEFAAAAwMmxrx5TTRcvXoyLFy/G7du347XXXmv1hpqcnGz1pNrq9u3bHuUDAAAAICIicvV6vX7QG33nnXdieXm5FVSdOXMmVldXI5/Pn5jBzzc3N2NoaCg2NjZicHAw6+IAAAAAHIleMpFDCaYe9s4778TKykqsra3FCy+8EGNjY3HhwoW+DmwEUwAAAMBJdOyCqYfdvn07arVarK2tRS6Xi/Hx8fjyl7981MU4VIIpAAAA4CTqJRN5pDGm9uvcuXNx7ty51uc7d+5kUQwAAAAAMrTvt/IdpLNnz2ZdBAAAAACO2J6DqXfeeScuX74cf/RHf3SY5QEAAADghNhzMPXyyy/HpUuX4urVq3HmzJmYm5uLd9999xCLBgAAAEA/6+lRvosXL8bbb78dSZJEoVCIV155JV566aX4+te/Hpubm4dVRgAAAAD60L7GmBoaGoorV67EzZs34+23346/+7u/i7GxsfjVX/3VeOuttw66jAAAAAD0oUce/Pzs2bPx+uuvx/e///14/fXX4+bNmzEyMhKXL1+Ob3/72wdRRgAAAAD60IG+le/cuXPx5ptvxvr6ely6dCnefPPNOHPmTHz1q181HhUAAAAAbQ40mNrq4sWLsbi4GEmSxNjYWNt4VO+9995h7RYAAACAx8ShBVNND49HVa/XY2JiIsbHx+Ott94yaDoAAADACXXowdRWZ8+eja997Wvx/e9/PyqVSty8eTM+//nPG48KAAAA4ATK1ev1etaFeOedd6JWq8XVq1cPdLvlcjny+XxERKRpGjMzM4+8zl63ubm5GUNDQ7GxsRGDg4P7PgYAAACAx0kvmcixCKYOQ7lcjohoBUe1Wi2WlpZiYWFh3+v0sk3BFAAAAHASCaYiYnh4OO7cudPq3RQRkcvlYqfD3W2dXrb52AdTr7wScepURKXyYNr0dMT9+xFvvfVgme9/P+IHP4jI5SI++cmI73434u//PuLJJyM+97mIzc2I//7fI+7dixgYiPiZn4lI04gf/7ix/VOnIk6fjvj448Z6uVzE889H/OQnjc+nfvq06dNPR3z4YWM7TzwRUShE/NIvRfzxH0f8zd9E1OuNdZ58srH8M89EfOELjfL88IeNac39/93fNfb3xBON7T7zTKOM9+839v3BB41yPPvsg22trbWX7+mnG/9eW2usNzra+FyvR3z60xEvvth+nnY7l83lvvOdiI2NiKGhxvn80z9tbPO55yI+8YnGvB//uDHtiScizpxplL15nu7fb0z/0Y8elPWjjxr/bmqey4jG9CeeiPjH/7jx+W/+pvFnc7Oxv8HBxjK/9msRf/iHjX3fv99YdnCwse16/cH+n3mmMe3jj9v3/fTTje/hmWcayzWXiWgs15y+ttb4u1BoHOtLL0X83M9tr3d7rZv1esTf/m1j2g9+0Phun3oq4md/tvH5hz9sL0dExPvvP6ibzfPVHAuvecydzuf9+w/OebEY8Rd/0TiXf/mXjfP23HPb6+ZPftKo01u/q+b5PH268d00v6tnn31Qtmbd/PSnG8fZbF/PPvvgXH/4YWPZX/iFxjlo1qsvfnF7vdvL+fzOd6LND37Q2MeTTza2m6aN42i2ty98IWJ1tTHt/v0H7WvruWwe//37jTI0y37/fuM8PPXU9rr54x83Pn/8ccQLLzTq5h//8YPzsJdtP/FEY/1crnHe7t3rXjc//enGuWu2s0KhMe+Tn2ys32tbby7z3/5bo44MDTXO5UcfNco1NNRYrvlzqlk3nnmm8T12qpsPt/VOx3z/fqNN/8ZvNNb9i79o/HxpttWH62bzZ/JO9X7r+fzJTxrn8umnG+sNDTXa8717D+r0Cy801ml+N2fONJZ71Lq59TpUrzfaXPNn/OBgY58//OGDn+H/8B8+qL9br01/+7eN9T/xiYh8/sF3cPp0o+wffdRog/V64zgfrpvNe4HmtejTn4741Kfar0Of+MSD77X5M+OFFxrnqrmvgYEH/753r/F38zw1t/OpTzWOrXmNjYj4xV98cAyPch36d/+u8e8vfKFxnpLkwXXok59sfG4ew1NPNc5VswzNujkwsP061O1aMTDQ+K7Gxxv7/M//ufEd/uhHD7bfbBfN77b5XT71VPv53K2t53IPzuXW7bz0UuM7bNbNp5+OOHcu4s/+rLHff/pPH71uHuQ9UvN8PvlkY97TTz+om5/4RMQ/+Ae93yM17ykervcR7pHcI/XHPdITTzR+nkRsv6Y3HdQ90upqo6w/8zONz2naaCPNsrz44uHcI73/fqOdP/lk41z+2q9FLC4+uAZu/Z6aPx/v32+cs63nsvk99vs9UvOYm+fzN36jsd96PeLGjcbPiuN6j/SY6CkTqR+Bcrlcf/HFF+u3b99uTdv674O2urpa73RoEVFfXl7e1zq9bnNjY6MeEfWNjY19HMExcOVKvZ7PN/7u9Lk57Zln6vUnnqjXc7nGn0bzOpo/AwMP/p3LtX/O4k8u1zgXzzyz/Tztdi6b0598sn17D2//sMuf5fnrVJ5Pfapzvdtr3Xz66WzqRZb1cad68+STnevdXs7nk09m084fPoaHz2sW5RkYaNSt/bT15rRPfSq7czk0lE3d3OnPo9TNrdehLOvmE080/hxlPXy4LAdxHWput1N7P8xzfByu452OdWDgYOqmeyT3SAddHvdIB1eeXK4/7pGOw5/H+R4p4sF3mfV5fLhMD9fNx0QvmciR9Ji6du1ajIyMxNjYWJw9ezYiIm7fvh2Li4sxNzd34D2KarVaTE5OxsOHNjw8HNeuXYupqame18nn8z1t87HvMRUR91+Zjh//m6W4Vz8VubgfNz5fjD/4YqVtmX/2nen4hdU/jCfufRS5LtuJiGietZ2W2Yt6h23czw3EH7/0v8X3P/nF+Gd/+mo8ee+HB7btre7lBuJU/d4Oy+TiJwNPx38Z/Y2O52n8vaWoR/dz2VzuF//i/46B+k+67me3cu5Vt+38ZODp+IP/6c24eHM28j/8m31t+36cilNxf5dlcpGLetdjuZc7HT9+4rmu9W638/mgbv4wctH9x9xe6+Zu573b/B8MfSH+/f/4Lw+1bjbOYvdzGRHx8akn409e/Odd691ezududXMvZd2L3ermi3/7nfilP38rTu3wvT5K+Xarm/XIxd8/QltvLvPUT96Pgbi3Yzl2O8a9HE+3Op4+/an45oX5HevmXup9p233ssxB1M29XYdyO/4seBTNuhkR+27rezmX93IDkavf36FeHMx16It//q92rJsHdR3qdq1oHMP/Gr+w+m8e6edmxKP97GzeXxx+3dy9rHvhHsk9UoR7pO3zj8c9UuNcPvp1aLe6+T//P78bn974s57W3ev8iP68R+q0zHG5R6qffjJO/eY/b+9p9hg5do/yffOb34zz58/H5z//+W3z3nrrrXjllVcOdH/dQqbR0dGYnZ2N6enpntcpFAo9bbMfgqn/+l8jPvuPzkQu7kc9TsX/8sW1jsv92++ciWfubcZA/eOO8+/lTsdHA4Px3Md3d/xhfC93+qc3NZ2XqUcu7ucGtu3n/dMjrbK98b3JOH+31nHb3crXnB8RXfdfj1x8cDofEbl4/uP1rtv4aGBwx/O027lsLrfTuXr/9Mgjn8udjvfW8ER87R8tR0TEt//T9h+Te/mePjw9FM9+vBERscP5HI6I6Ho+Pzid37XePWrdbJZjL+czInbczv3WTXm9bfrLv9y4+dypbu6l3p+qNy7ODy/XnP/RwGDXc9k8zt3q3V7O507nqtt52KqX492pbv7b75zZdryPsu2ty+xWNw+irTeXee7jtGs5VoZfjrG77+xYNz84nY9nP97Y9Zgjtrf3L/1y49871c2I3et9p21vXeaD08Nd2+FB1s3d23o+nvs43bWtb93Gw2FWtzq2tW7u5zq09Vzu9vPqz58f67j95j4O6jq028+Tg7imd7tWNM/no/7cjNi9bkbU4/mP73bcztb7i04Oom4e1D1St/a6n3ukTvU+wj1ShHukTtwjPVhmt7r5zL3NI7lHeuc/ndo2f7fvabd23tx/P94jdao/x+Ue6ePnh+PJze5187jrJRM5fRQFunjxYvzKr/xK3L17N0qlUly6dKlVsDRNj6IIERGxvt65AT3KOvvZ5uPixTem49Qn7se9+qkYyN2Pxfx0/OX/0Z7W/uz/OR3PDvwoch93D30G6h/Hcx+v75pY7xQcRTR+WHRa5rn6Rizmp+ODn/9i/Ox/+ZN9bXsv+36uvvlgDJuO27gXz57+Udfz9MwnfvpM8/3O57K53HPxwY4/ZLtdCNrLsv/jPffRn8S/+uVvxP/wf83ua9u5qMez9zZ3PIZc1OO5e3cbHVS7eLb+Qdx/+rl9n88HdbP7d5aL+oGcz251Mxf1+Pd3fi7+5jf/5b7rZrdtt82Pe/HcDseRi3o8l/tgx3q3l/O5W93crawRj3a8zbr53LvfadSfA9z21mV2q5sD8WhtvbnMqR++v2M5xtKdb7gioutN29btdDvm//j/fjr+v/99/tDqZnOZnX7+H1Td3O061CjHzqFUxPZjfnj53epmROzrfO75XNY3Yix9p+syB3Ydur+xczkO6Jre7Vpx7qM/icX8dIx8dMh1c5e23ry/OMy6eVD3SN3m7+ceaa/1fuvy7pH2tm33SO3TT9I90vMfr/90HKId7i8O4B7pU7//ux3L8aj3Z839u0dqOKp7pCd+9EFjLK7HtMdUTw7tgcItXn311XqlUqlXKpV6sVis53K5+osvvlh/8cUX62+88caB788YUwfA+Am9/zF+wsGXx/gJj/49GmPq4P88zuMnGGPq8OqmMaYe/VxmfR3vdKzGmDq48+ke6WDL4x7p4MpjjKmD+/M43yNFGGPqgPWSiRxJj6mJiYm4ePFiRERcuXIlIhqPzt2+fbvjY3WPqlAoRD6fjyRJovDQ2xYmJib2vU6v23ys3b/feGtGM51t/t1820jz3+Pjx+utfFvfEpPlW/kePk+7ncvm57NnvXEmovMbZ/ZTN+v1k/XGmb28la9TvdvL+fzp+IAtWb6Vr/kmmePwVr791M1isfHGmQhv5dv6xpn91s2t16F6Pfu38g0MHI+38u33OtR8i5S38j14K99B1M3j8Fa+ne6RjvKtfO6R3CMd17fyHdQ9krfyPR73SI/TW/kerpt96MjGmBodHY2f//mf3zbvW9/6Vnz5y18+8H2Wy+XI5/Ot4Ktarcby8nIsLCxERESSJFGtVmNmZmbP6+w2f6t+GGMKAAAAoFe9ZCKnjqJAFy9ejLW1tfj617/emnb79u0YHx+PJEkOZZ8zMzORpmlUq9WoVqtx48aNtgCpVqttC5R2W2e3+QAAAADs3ZH0mOrm2rVrceHChTh37lxWRTg0ekwBAAAAJ9Gh95i6dOlSDAwMxMDAQPz2b//2vgoZ0Rhvqh9DKQAAAAB213MwNTc3F+Pj43Hz5s24fv16fPe7343x8fG2+S+99FIMDAzEmTNn4td//dfje9/73oEWGgAAAIDHX8+P8r322mvx+uuvt00rl8tx69atSJIkhoeHY3JyMtbW1iJJkqjVarGxsRGVSiV+67d+60ALf5x5lA8AAAA4iXrJRE73uvEXmq+O3WJmZiYuXLgQr776arzyyivb5ler1dab7E5SOAUAAABAdz0/yre2ttZ13sjISMfpU1NTkSRJ/N7v/V5sbm72uksAAAAA+lDPwVShUIj33nuv4/QkSbqul8/nY2lpKa5evdrrLgEAAADoQz0HU1euXIk333xzW8+nNE0jn8/vuO7Zs2e79qoCAAAA4GTpeYypiIjXX389Xn311bh06VJ86UtfioiIt99+e0/r5nK5/ewSAAAAgD7Tc4+ppjfffDPu3r0bX/3qV+Nb3/rWntfbaYwqAAAAAE6OffWYarp48WJcvHgxbt++Ha+99lqrN9Tk5GSrJ9VWt2/f9igfAAAAABERkavX6/WD3ug777wTy8vLraDqzJkzsbq6Gvl8/sQMfr65uRlDQ0OxsbERg4ODWRcHAAAA4Ej0kokcSjD1sHfeeSdWVlZibW0tXnjhhRgbG4sLFy70dWAjmAIAAABOomMXTD3s9u3bUavVYm1tLXK5XIyPj8eXv/zloy7GoRJMAQAAACdRL5nII40xtV/nzp2Lc+fOtT7fuXMni2IAAAAAkKF9v5VvvzY3N7dNO3v27FEXAwAAAICMHVkwtbm5Ga+++qq38gEAAAAQEUfwKN97770X8/PzUalUol6vt97UBwAAAMDJdmg9pt577724fPlyjI6OxsLCQrz88ssxPT19WLsDAAAA4DFz4MHUu+++2wqklpaW4uLFi7G6uhpvv/12nD9//qB3BwAAAMBj6sAe5Xv33XdjdnY2arVa1Ov1mJqaivn5eQObAwAAANDRIwdT3/72t2N2djZWVlaiXq/H9PR0zM7OCqQAAAAA2NG+g6lOgdT8/HwMDQ0dZPkAAAAA6FM9B1Pf+ta3YnZ2NpIkiXq9HjMzMzE3NyeQAgAAAKAnPQVT3/zmN6NYLMbw8HDMzMzEa6+9JpACAAAAYF96eivfxYsX4+bNm1EsFmN0dFQoBQAAAMC+9RRMRUSMjY3Fm2++GS+//HK89tpr8dZbbx1GuQAAAADoc/se/Pzs2bPx+uuvx507d+K1116LF154IX7nd37nIMsGAAAAQB/bdzDV1AyoNjY2WgHV9PR0DA4OHkT5AAAAAOhTjxxMNQ0NDbUCqt/93d/NPKAql8uRz+cjIiJN05iZmdnTOhERq6urERGxsLDQmler1WJhYSEmJyejUCjE8vJyjI+Px9TU1MEXHgAAAOAEOLBgqmlrQLWwsBC5XO7IA6pmwDQ9PR0RjVCpVCq1BU0Pm52djfn5+dbnUqkUk5OTsby8HBGNcKtWq0W1Wo1CoRCzs7NCKQAAAIBHkKvX6/XD3skbb7wR6+vrUa/X44033oh79+4d6v6Gh4fjzp07rR5TERG5XC66HWqaplEsFmNpaam1zsrKSpw/fz5WV1ejUChEtVqNiYmJtm3uZHNzM4aGhmJjY8NjjQAAAMCJ0Usm0vNb+fbja1/7Wly9ejVGRkZiaGjoUPeVJEmkadoxQKrVal3Xu3nzZiRJ0vpcKBQiohFaAQAAAHDwDvxRvp3MzMzsaaynR7E1XNoqn893DZny+XzcvXu3bVozxGoGVBERi4uLMTIyEuvr67G6utr26B8AAAAAvTnSYCpLzUBpr65evRoLCwutnldjY2MR8SCoqlQqrcf/AAAAAOjdsQ+mqtVqXL9+fdfl5ubmWuFRJ72EUrOzs3H58uXW4OkR7T2nIiIuXboUpVKp62ODAAAAAOzs2AdTU1NTPb397uEAqSlN067ztqpWqzE6OtoWSjWnby1HM4xKkmTHQAwAAACAzo5k8POjVCgUIp/PdxxramJiYsd1m+NKNUOpNE1bg6kXi8W2bTbHq9pL2AUAAADAdn0XTEU0Huvb+ga+arXa1gMqSZIol8tt66ysrMTKykqMjY1FkiSRJElUKpUYGRmJfD4fMzMzbSFUpVKJqakpj/EBAAAA7FOuXq/Xsy7EYSiXy60g6caNG21v0KtUKjE/Px+rq6sR0ej9dPbs2Y5v7WuenjRNo1KptKavra3t+Fa+zc3NGBoaio2NjRgcHDyIQwIAAAA49nrJRPo2mMqaYAoAAAA4iXrJRPryUT4AAAAAjj/BFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkInTWRfgsJTL5cjn8xERkaZpzMzM7Lh8rVaLhYWFmJycjEKhEMvLyzE+Ph5TU1P73iYAAAAA3fVlj6lyuRwREdPT0zE9PR1jY2NRKpV2XCdN06jValEqlaJUKsXo6Oi2UKrXbQIAAADQXa5er9ezLsRBGx4ejjt37rR6N0VE5HK52OlQq9VqTExMtK3zKNvc3NyMoaGh2NjYiMHBwf0cBgAAAMBjp5dMpO96TCVJEmmadgyYarXasdkmAAAAwEnXd2NMJUnScXo+n480TXdcd3FxMUZGRmJ9fT1WV1djfn7+kbcJAAAAQGd9F0x10wycuhkbG4uIiEKhEBERlUolisViLC0t7XubAAAAAHR37IOparUa169f33W5ubm5VrjUyW4BUjOQarp06VKUSqUde0QJpQAAAAD279gHU1NTU21vx9vNwwFTU5qmXedFNAKwrftpjieVJMm+twkAAABAd303+HmhUIh8Pt9xXKiJiYmO66RpGsVisW2dZk+pQqGwr20CAAAAsLO+C6YiGo/1bX1bXrVajenp6dbnJEmiXC63Pufz+ZiZmWnr/VSpVGJqaqrVc2q3bQIAAADQm1y9Xq9nXYjDUC6XW0HTjRs3Wm/Yi2iETvPz87G6utqalqZpVCqV1ue1tbW2dXbb5sM2NzdjaGgoNjY2YnBw8ECOCQAAAOC46yUT6dtgKmuCKQAAAOAk6iUT6ctH+QAAAAA4/gRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGTidNYFOCzlcjny+XxERKRpGjMzMzsuXywW4/Lly1EoFFrrNRUKhajVarGwsBCTk5NRKBRieXk5xsfHY2pq6pCOAAAAAKC/9WUwVS6XIyJieno6IiJqtVqUSqVYWFjous7KykpUq9Vt06empmJpaSnSNI1arRbVajUKhULMzs4KpQAAAAAeQa5er9ezLsRBGx4ejjt37rT1fMrlcrHToZbL5W29qiqVSivcqlarMTExsa03VTebm5sxNDQUGxsbMTg42PMxAAAAADyOeslE+m6MqSRJIk3TjgFSrVbrut7DvZ9qtVpcuHDhoIsHAAAAwE/13aN8SZJ0nJ7P5yNN067rFQqFtm0kSRITExNtyywuLsbIyEisr6/H6upqzM/PH0iZAQAAAE6ivgumumkGSnsxPz+/bTyqsbGxiHgQYFUqlSgWi7G0tHSwBQUAAAA4IY59MFWtVuP69eu7Ljc3N9cKjzrZayi1srLScfrWHlUREZcuXYpSqdT1sUEAAAAAdnbsg6mpqame3n73cIDUlKZp13lbLSwsxOjo6Lbp1Wq1rRzNMCpJkh0DMQAAAAA667vBzwuFQuTz+Y5jTT08ZlQntVptWw+oNE2jWCy2bbM5XtVewi4AAAAAtuu7YCqi8Vjf1jfwVavVmJ6ebn1OkiTK5XLHdZMk2RY25fP5mJmZaZteqVRiamrKY3wAAAAA+5Sr1+v1rAtxGMrlcitIunHjRtsb9CqVSszPz8fq6uq29UZHR2NpaWnb43lpmkalUml9Xltb2/GtfJubmzE0NBQbGxsxODj4qIcDAAAA8FjoJRPp22Aqa4IpAAAA4CTqJRPpy0f5AAAAADj+BFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFN09eGHH0Yul4tcLhcffvhh1sWBI6cNgHYA2gBoB6ANHC7BFAAAAACZOJ11AQ5LmqaxuLgYS0tLsby8vKd1yuVy5PP51vozMzM9zQcAAABg7/qyx9TKykosLi5Gmqaxvr6+p3XK5XJERExPT8f09HSMjY1FqVTa83wAAAAAepOr1+v1rAtxWKrValy9ejVu3bq167LDw8Nx586dVo+oiIhcLhfN07Pb/Idtbm7G0NBQbGxsxODg4CMdR1Y+/PDDeO655yIi4oMPPohnn3024xLB0dIGQDsAbQC0A9AGetdLJtKXPaZ6lSRJpGnaFjo11Wq1XecDAAAA0DvBVDSCqU7y+XykabrrfAAAAAB617eDnx+EkZGRWF9f79hTauv8TpqP+G1ubh5W8Q7d1tdgbm5uxr179zIsDRw9bQC0A9AGQDsAbaB3zSxkL6NHHftgqlqtxvXr13ddbm5uLsbGxg5037sNnL7T/Pfffz8iIj73uc8daJmy8pnPfCbrIkCmtAHQDkAbAO0AtIHevP/++zE0NLTjMsc+mJqamoqpqalD3UehUOg4PU3TKBQKu87v5DOf+Uz81V/9VTz//PORy+UOrKwAAAAAx1m9Xo/3339/T0HesQ+mjkKhUIh8Ph9JkmwLmiYmJiIidp3/sFOnTsVnP/vZwykwAAAAwDG2W0+ppr4e/Lzbo3ZJkkS5XG6bNjc31/aGvWq1GtPT03ueDwAAAEBvcvW9jET1mEmSpDU21crKSszMzMT4+HjrkcBKpRLz8/Oxurratl65XG71iLpx40bMz8/3NB8AAACAvevLYIqDUS6XW28kTNM0ZmZmsi0QHKJarRYLCwsxOTkZhUIhlpeX2wLtCG2C/pOmaSwuLsbS0lIsLy9vm79bndcmeNzt1AZcFzgpmk+SNP/TfmFhYdt81wL63U7twPXg8Bljio6aDbP5uGKtVotSqbTtQgX9Ik3TqNVqUa1Wo1AoxOzs7LaLTYQ2Qf9YWVmJmzdvRpqmHR99363OaxM87nZrA64LnASzs7NtT4GUSqWYnJxsBbWuBZwEu7UD14PDp8cUHQ0PD8edO3daqW9ERC6XC9WFflWtVmNiYqKtzm+lTdCvqtVqXL16NW7dutU2fbc6r03QL7q1AdcF+l2aplEsFmNpaalVj1dWVuL8+fOxuroahULBtYC+t5d24Hpw+Pp68HP2J0mSSNO0Y8PbOgA8nBTaBCfNbnVem+Ck0wboFzdv3owkSVqfm+PppmnqWsCJsVM72I12cDA8ysc2WxvlVvl8fk+NEx5Xi4uLMTIyEuvr67G6utrq0qtNcNLsVue1CU4K1wX6WT6fj7t377ZNa/4iXSgU4ubNm13Xcy2gX+zWDppcDw6XYIo9azZE6EdjY2MR8eACVKlUWt16u9EmOGmadb5bV3Ztgn7iusBJdPXq1VhYWOj6cz7CtYD+93A7cD04fIIp9kzDop9t/R+RiIhLly5FqVTa8X86tAlOmt3qvDZBP3Fd4KSZnZ2Ny5cvtwZw7sa1gH7WqR24Hhw+Y0yxzcMNrylN067z4HFXrVbbPjf/hyRJEm2CE2e3Oq9NcBK4LnCSVKvVGB0dbXvFvWsBJ02ndtCcvpXrwcETTLFNoVCIfD7f8XnZiYmJDEoEh6v5No6tdb75PyDNGy9tgpNktzqvTdDvXBc4SZrj6TR7iDTHj3It4CTp1g5cD46GYIqO5ubm2t4iUK1Wd+3WC4+rfD4fMzMzbf+rUalUYmpqqvU/ItoE/apbV/Pd6rw2Qb/o1AZcFzgpVlZWYmVlJcbGxiJJkkiSJCqVSoyMjESEawEnw07twPXgaOTq9Xo960JwPJXL5VYDvHHjRuvNA9CP0jSNSqXS+ry2tratzmsT9JMkSaJarcb169djZWUlZmZmYnx8PKamplrL7FbntQkeZ7u1AdcF+l2apnH27NmO4+Rs/RXRtYB+tpd24Hpw+ARTAAAAAGTCo3wAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmTmddAAAADtbs7GycOXMm1tbWIk3TKBaLMTExkXWxAAC2EUwBAPSR8+fPx7Vr12JsbKw1rVgsxoULFyKfz2dXMACADjzKBwDQJyqVShQKhbZQamVlJarVaiRJkmHJAAA6E0wBAPSJ5eXlSNO0bVo+n4+ZmZm2sAoA4LgQTAEA9IlCoRC1Wi1GR0ejVCpFtVqNQqEQ8/PzWRcNAKCjXL1er2ddCAAADsbs7Gzbo3v5fD7u3LljfCkA4FjSYwoAoI/Mz8/H6upq1Ov1WFhYiDRNY3Z2NutiAQB0JJgCAHjM1Wq1GB4ejpWVlbbp09PTMTExobcUAHBsCaYAAB5zCwsLEdEYY+phSZJEqVQ66iIBAOyJYAoA4DE3OTkZS0tL23pGzc7ORqlU6hhYAQAcBwY/BwDoA5VKJVZXV+PMmTMREbG2thaTk5MxMTGRcckAALoTTAEAAACQCY/yAQAAAJAJwRQAAAAAmRBMAQAAAJAJwRQAAAAAmRBMAQAAAJAJwRQAAAAAmRBMAQAAAJAJwRQAAAAAmRBMAQAAAJAJwRQAAAAAmRBMAQAAAJAJwRQAAAAAmRBMAQAAAJAJwRQAAAAAmfj/ASVAjDiRz05wAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1200x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(4.6202e-12, dtype=torch.float64)\n",
      "tensor(4.6143e-14, dtype=torch.float64)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABKYAAAGGCAYAAABBiol3AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAM7NJREFUeJzt3V9sXNmdH/hfSe3x3yYvqTYyY4wDq+hOHj1NiVnsZh6CbhJBNg8LWFWSgUG8mdhidQ+wb+OiCOzjwlKV/TwWS0Ym8Ow+iCw3sA9BsM3qnodMJoAlVbexT0nMy8YY62SQIXXJtj2emW7VPtBVYonFP6Umdani5wM05Lr3nnvPvTxHZX51zrmFTqfTCQAAAAB4xs7lXQEAAAAAzibBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5eCHvCpykLMtieXk5VlZWYnV19Uhl6vV6JEnSK1+tVofaDwAAAMDRjOyIqXa7HcvLy5FlWWxubh6pTL1ej4iI+fn5mJ+fj+np6ahUKkfeDwAAAMDRFTqdTifvSpykZrMZN2/ejAcPHhx67MTERKyvr/dGREVEFAqF6D6iw/YDAAAAcHQjO2JqWGmaRpZlfaFTV6vVOnQ/AAAAAMMZ6TWmhpGm6cDtSZJElmWH7h/k0aNH8bOf/SxefPHFKBQKx1VVAAAAgFOr0+nEBx98EF/4whfi3LmDx0QJpg4xOTkZm5ubA0dK7d4/yM9+9rP44he/eIK1AwAAADidfvrTn8Zv//ZvH3iMYOoQhy2cftD+F198MSJ2fhBjY2PHWi8AAACA02h7ezu++MUv9nKRgwimfq1YLA7cnmVZFIvFQ/cP0p2+NzY2JpgCAAAAzpSjLGtk8fNfKxaLkSTJwLWkZmdnD90PAAAAwHBGPpjab6pdmqZRr9f7ti0uLva9Ya/ZbMb8/PyR9wMAAABwdIVOp9PJuxInIU3TaDabcffu3Wi321GtVmNmZiZKpVJERDQajajVarG2ttZXrl6v96bm3bt3L2q12lD7d9ve3o7x8fHY2toylQ8AAAA4E4bJQ0Y2mDoNBFMAAADAWTNMHjLyU/kAAAAAOJ0EUwAAAADkQjAFAAAAQC4EUwAAAADkQjAFAAAAQC4EUwAAAADkQjAFAAAAQC4EUwAAAADkQjAFAAAAQC4EUwAAAADkQjAF8Iw0Go0ol8t5VwMAAODUEEwBPCMPHjyIZrOZdzUAAABODcEU7NJsNiPLsryrwYhaWlqKTqfzTK+pTQMAAKeZYAp2KZfLcf/+/byrAcdGmwYAAE4zwRScdf/m30T84Af9237wg53tAAAAcIIEU6dcpxPxq1+d7v+OY2ZSs9mMS5cuRaFQiEuXLkWr1YqIiIWFhSgUCpGmad/xc3NzMTU11fucpmmUy+WYmJjonaPdbu+5Trvdjrm5uZiYmIiJiYmYm5uLdrsd5XI5CoVC79yFQqH3eeSdOxfxx3/8OJz6wQ92Pp/L/6+Ho/xcG41G7+fY/dlOTU312tAwx9Xr9Zibm4uIiEqlEhMTE71pcE+Wq9frvXJHbaf1ej0mJib6Pnfv6cn2n2VZ796npqai0WgM9WyO0qb363cAAADPygt5V4CD/c3fRPyzf5Z3LQ727/5dxKc+9fTlG41GVCqVWFpaisuXL0er1Yq5ublYW1uLWq3WCwTW1tZ6x7dard7niJ1fsCcnJ2NlZSWKxWIsLCzEa6+9Fuvr65EkSURE75f/UqkUKysrERGxuroad+/ejTt37sTi4mJcunQpVlZWYnZ29ulv6Hnz9a/v/PnHfxzxJ38S8eGHEb//+4+35+goP9csy6LVasX169djcXExKpVKLCws9NpMsVgc6rg0TePSpUuRZVksLi5GkiS9NlmtVqNWq0WaprGwsBD37t2LlZWVI7fTJ21sbES73Y7r169HrVaLiJ1ArFwuR7FY7NWzVqtFpVKJ2dnZXj0PezaHtemD+l33GgAAACdNMEWusizr/XI8Pz8fERHT09OxsbERtVotlpaWYnV1NaampqJcLvd+Ue/+Mt5VrVb7znvnzp2YmJiI5eXl3nnL5XJfKBURfb+sd8+XJEkv9Dgzvv71x6HUCy+cilAq4mg/164HDx70/vfs7GxMTEz02tAwx6VpGrOzs33HVSqVXigVsdNGp6enY2pqKprNZpRKpSO10/3UarVeW1xYWOiFUKVSKSIiJicneyOiuuc77NkkSbJvmz5KvwMAAHgWBFOn3Cc/uTMi6TT75Cefvmx3UeZKpRKVSqVv3/T0dO9/d3/pbzabUa1We7+w76f7S3h3tEqappGmqV+49/ODHzwOpT78cOfzKQmndnvy53rQcaVS6dCpafsdt7udtNvtSNN0T/ssFosxPT0dd+/e7bXHYdtp1+XLl/vOGxG9KYURj/vC5ubmgfcScfiziTh6vwMAADhp+S8iw4EKhZ1pcqf5v+NYiunhw4fR6XT6/ts9YqVYLPZGdiwuLg48R6vVinK5HFNTU33r+EREb+0fU5QG6K4p9fu/H7G6uvPn7jWncnbQz/UgxWLxwCDnoON2t5NuiDM5OTmw7O51pY7STgcZNEJv0PWe9LTPpuuwfgcAAHDSBFPkave6PgdptVqxvLwcpVIpXnvttT37y+VylMvlmJubi9XV1Xj48OFTXedMevSof02pr3995/OjR/nWKw7/uR4kTdMjBZFPHvdkSHRQ23my7GHtdJD9po0eNp304zwb/QEAADgtBFPkqjsd6ubNm3v2dd+G1n372J07d3rrQ+2efpRlWTSbzbhz507Mz88PDCOKxWIUi8WBU/m619m9mPaZ8i//5d5pe1//+s72HB3l57qfNE2j2Wweuoj9UY6bnZ2NJEn2tJ12ux3tdjuuXbvWO9dB7fQ4HfXZ7Nemj9LvAAAAngVrTJG7O3fuxKVLl3rr3WxubsbKykrvz3K53LcQ9MrKSkxNTcWlS5d6izwnSRI3b96MJElicnJy4C/cS0tLMTc313vrWZZlsbq6Gvfv3+9NX0qSJO7evdsLsWq12tlbCP2UOOrPtWtubi4WFhYiy7K4fv16JEkycDrdUY/b7c6dO1EulyNiZ6RS9618u9vlYe30OA3zbPZr04f1OwAAgGfBiClyNz09HWtra5Gmabz22mu9AODOnTtRqVQiTdO4c+dO7/juL9iVSiXa7Xbv2O6IlevXr8fc3FzMzs7G1NRUr1z3TWu7j3vyl/DFxcVoNpu9OpCvo/xcuyqVSiwsLES5XI7Lly/HgwcPBoaKRz1ut1Kp1Gs7c3NzUavVYnFxMVZXV3vnPEo7PU5HfTb7temD+h0AAMCzUuh0Op28KzGqtre3Y3x8PLa2tmJsbCzv6sBIqtfrsbCwEIf9VXbU4wAAAPh4hslDjJgCAAAAIBeCKQAAAAByMfKLn9fr9b43U1Wr1QOPL5fLce3atSgWiwNfG99qtXqLaBeLxVhdXY2ZmZnegscAAAAAHM1IrzFVr9cjInphVKvVipWVlT2vfd9tamoq0jTds71UKsXKyko0m824fv16ZFkWxWIxFhYW9n3jljWmAAAAgLNmmDxkpIOpiYmJWF9f7xv5VCgUDlz8uF6v7xlV1Wg0euFTs9mM2dnZQ9/iFSGYAgAAAM4ei59HRJqmkWXZwACp1WrtW+7JKXmtVisuX7583NUDAAAAOPNGdo2pQdPxIiKSJIksy/YtVywW+86RpmnMzs72HbO8vByTk5OxubkZa2trUavVjqXOAAAAAGfJyAZT++kGSkdRq9X2rEc1PT0dEY8DrEajEeVyOVZWVo63ogAAAAAjbmSn8u3nqKFUu90euL1YLPaNqrp69Wo0m80DR2EBAAAAsNfIBlO7w6Pdum/TO8zS0lJMTU3t2d5sNvs+d9ew2m/qIAAAAACDjXQwlSTJwMDoyTWjBmm1WnsWTs+yLMrlct85uyOljhJ2AQAAAPDYyAZTERGLi4t9b+BrNpsxPz/f+5ymadTr9YFl0zTdEzYlSRLVarVve6PRiFKpNPDtfwAAAADsr9DpdDp5V+Ik1ev1XpB07969vjfoNRqNqNVqsba2tqfc1NRUrKys9BY778qyLBqNRu/zxsbGvm/l297ejvHx8dja2oqxsbHjuB0AAACAU22YPGTkg6k8CaYAAACAs2aYPGSkp/IBAAAAcHoJpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIxQt5V+Ck1ev1SJIkIiKyLItqtXrg8a1WK5aWlmJubi6KxWKsrq7GzMxMlEqlpz4nAAAAAHuN9Iiper0eERHz8/MxPz8f09PTUalUDiyTZVm0Wq2oVCpRqVRiampqTyg17DkBAAAA2KvQ6XQ6eVfipExMTMT6+npvdFNERKFQiINuudlsxuzsbF+Zpz3n9vZ2jI+Px9bWVoyNjT3tbQAAAAA8N4bJQ0Z2xFSappFl2cCAqdVqnZpzAgAAAJxVI7vGVJqmA7cnSRJZlh1Ydnl5OSYnJ2NzczPW1taiVqt97HMCAAAA0G9kg6n9dAOn/UxPT0dERLFYjIiIRqMR5XI5VlZWnvqcAAAAAOx15oKpwwKkbiDVdfXq1ahUKgeOiBJKAQAAAAxvZNeYejJg6sqybN99ETuLn+/WXU8qTdOnPicAAAAAe410MJUkycB1oWZnZweWybIsyuVyX5nuSKlisfhU5wQAAABgsJENpiIiFhcX+96W12w2Y35+vvc5TdOo1+u9z0mSRLVa7Rv91Gg0olQq9UZOHXZOAAAAAI6m0Ol0OnlX4iTV6/Ve0HTv3r3eG/YidkKnWq0Wa2trvW1ZlkWj0eh93tjY6Ctz2Dl3297ejvHx8dja2oqxsbFjuycAAACA02qYPGTkg6k8CaYAAACAs2aYPGSkp/IBAAAAcHoJpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFy8kHcFTlq9Xo8kSSIiIsuyqFarRyoTEbG2thYREUtLS719rVYrlpaWYm5uLorFYqyursbMzEyUSqXjrzwAAADACBvpYKobMM3Pz0fETqhUqVT6gqYnLSwsRK1W632uVCoxNzcXq6urEbETbrVarWg2m1EsFmNhYUEoBQAAAPAUCp1Op5N3JU7KxMRErK+v90ZMRUQUCoXY75azLItyuRwrKyu9Mu12Oy5duhRra2tRLBaj2WzG7Oxs3zn3s729HePj47G1tRVjY2PHcEcAAAAAp9swecjIrjGVpmlkWTYwQGq1WvuWu3//fqRp2vtcLBYjYie0AgAAAOD4jOxUvt3h0m5JkuwbMiVJEg8fPuzb1g2xugFVRMTy8nJMTk7G5uZmrK2t9U39AwAAAOBoRjaY2k83UDqqmzdvxtLSUm/k1fT0dEQ8DqoajUZv+h8AAAAARzeyU/n2M0wotbCwENeuXestnh6xE0jtHj119erVaDabpvoBAAAADOnEg6n33nsvtre3T/oye+wOj3bLsmzffbs1m82YmpqKarW6Z/tu3ZFU+00dBAAAAGCwEw+mvv3tb0eSJPHyyy/HG2+8EW+++Wa8//77J33ZKBaLkSTJwMBodnb2wLLddaW6I6WyLOstpl4ul/vO2R0pdZSwCwAAAIDHTjyYWl5ejkePHsX3vve9GB8fj9u3b0exWIwLFy7EtWvX4p133jmxay8uLva9ga/ZbPZNy0vTNOr1el+Zdrsd7XY7pqenI03TSNM0Go1GTE5ORpIkUa1W+0KoRqMRpVJp4Nv/AAAAANhfodPpdPK4cLPZjOXl5ciyLNbX16PZbMZXvvKVY79OvV7vBUn37t3re4Neo9GIWq0Wa2trEbEz+unixYsD14vqPqYsy6LRaPS2b2xs7PtWvu3t7RgfH4+tra0YGxs7rlsCAAAAOLWGyUNyC6YiIu7cuRPXr1+Pdrsd8/Pz0Ww240tf+lJe1Tl2gikAAADgrBkmDznxqXw3btyIc+fOxczMTHz3u9/tW19qa2srIiKmp6fj/v37exYWBwAAAGB0nXgwNTMzE48ePYqFhYX40Y9+FMViMc6fPx8XLlzoHfPee+9FRMTFixdPujoAAAAAnBIvnPQFkiSJd955J0qlUpRKpYiIWF9fj8nJyRgfH4+IiFKpFOVyOWZmZk66OgAAAACcEs9kjamtra148OBBvPrqqyd9qVPFGlMAAADAWXOia0xdvXo1zp8/H+fPn48/+IM/OFKZ8fHxMxdKAQAAAHCwoYKpxcXFmJmZifv378fdu3fjRz/6Ud/0u8XFxXj55Zd7a0h97Wtfix//+MfHXmkAAAAAnn9DTeW7ceNG3Lp1q29bvV6PBw8eRJqmMTExEXNzc7GxsRFpmkar1Yqtra1oNBrxjW9849grf9qZygcAAACcNcPkIUMtfv7SSy/t2VatVuPy5cvx+uuvxze/+c09+5vNZszPz0dEnMlwCgAAAIDBhprKt7Gxse++ycnJgdtLpVKkaRrf+973Ynt7e7jaAQAAADCyhgqmisVivP/++wO3p2m6b7kkSWJlZSVu3rw5dAUBAAAAGE1DBVPXr1+P27dv7xn5lGVZJElyYNmLFy/uO6oKAAAAgLNnqDWmIiJu3boVr7/+ely9ejVeffXViIh46623jlS2UCgMezkAAAAARtRQI6a6bt++HQ8fPow33ngj3nzzzSOXO2iNKgAAAADOlqFHTHVduXIlrly5Eu+++27cuHGjNxpqbm6uN5Jqt3fffddUPgAAAAB6Cp1Op3OcJ3z77bdjdXW1F1RduHAh1tbWIkmSM7f4+fb2doyPj8fW1laMjY3lXR0AAACAEzdMHnLswdST3n777Wi327GxsREvvfRSTE9Px+XLl89EUCOYAgAAAM6aUxVMPendd9+NVqsVGxsbUSgUYmZmJr761a8+yyo8M4IpAAAA4KwZJg956jWmntYrr7wSr7zySu/z+vr6s64CAAAAAKfAU72V7zhdvHgx7yoAAAAAkIMjBVNvv/12XLt2Lf70T//0pOsDAAAAwBlxpGDqtddei6tXr8bNmzfjwoULsbi4GO+9994JVw0AAACAUXbkqXxXrlyJt956K9I0jWKxGN/85jfj5Zdfju9+97uxvb19knUEAAAAYAQNvcbU+Ph4XL9+Pe7fvx9vvfVW/NVf/VVMT0/HP/2n/zS+//3vn0QdAQAAABhBH2vx84sXL8atW7fiJz/5Sdy6dSvu378fk5OTce3atXjnnXeOq44AAAAAjKBjeyvfK6+8Erdv347Nzc24evVq3L59Oy5cuBBvvPGG9agAAAAA2OPYgqndrly5EsvLy5GmaUxPT/etR/X++++fxCUBAAAAeM6cSDDV9eR6VJ1OJ2ZnZ2NmZia+//3vWzQdAAAA4Aw70WBqt4sXL8a3vvWt+MlPfhKNRiPu378fX/rSl6xHBQAAAHBGFTqdTifPCrz99tvRarXi5s2bJ3L+er0eSZJERESWZVGtVj92maOec3t7O8bHx2NrayvGxsae+h4AAAAAnhfD5CG5B1MnqV6vR0T0gqNWqxUrKyuxtLT01GWGOadgCgAAADhrTlUw9d5770WxWMwlmJmYmIj19fXe6KaIiEKhEAfd8mFlhjnncx9MffObEefORTQaj7fNz0c8ehTx/e8/PuYnP4n4r/81olCI6HQi/uIvIj78MOITn4gYG4v467/e+e/v/m7nmH/wD3aO/8UvIj76KOL8+YjPfz7iL/9yp/ynPhWRJBF/9Vc753nhhYiXXor45S8jPvhg55hPfzriH/2jnTr8t/+289/f/E3Eb/zG47r+1m/t/PkXf7Fz/Yidc3X/7NbnM5/Z2f/o0c79nj//+LoffbTz5/h4xNbW4/P85m/u3Nv2dsR//+872/7xP358D7/1WxFf/nL/czrsWXaP+7f/dud//8N/uPOcfvKTnXp+/vM710zTnfoVChEvvrhT/24dzv16du758xG/+tXjY375y516dY/59Kcf3/OjRzvP/Pd+b2f/f/kvEX/+5zvXePHFx8/y7/29iHv3drZ3n1333IXC4+e3+9yf+ETE3/7tTn0+/emdcuPjERsbO8+209kpMzW18zPsPt8LF3bu/z/9p53P//yf731O2qa2qW1qm9qmtqltapvaprapbWqb3Trsbnvb2zvHfPKTj48f1H6ex7b5nBgqD+mcsHK53CkUCp0vf/nLnddff73zwx/+sLO+vn7Sl+2sra11Bt1eRHRWV1efqsyw59za2upERGdra+sp7uAUuH6900mSnT8Hfe5u+8xnOp1PfKLTKRQ6nZ3u9ez+233N8+efbR3On99bl0Jh51l85jN7n9Nhz7K7vXve7vkOuuZxP8vx8Wf/MzzsZ3v+/ODnpG1qm9qmtqltapvaprapbWqb2qa2edba5nNimDzkmU3la7Va0Wq1ot1uR6vViomJiZidnY1KpRKvvvrqiVxvbm4unry9iYmJuHPnTpRKpaHLJEky1Dmf+xFTEfE3/+t8fHR3JR51zkUhHsW9L5XjT3630XfMv/iz+fgf1/6v+MRHv4zCAefqRCEKcTLN7W/Pfzr+5H+6HV/+yz+L3/3P/zrOx0dDn6Nbs4Pu4aPC+Sh0HsW5fe+jEH97/tPxH6d+b+Bzmnl/JTqx/7PsHnfYPTyKwgF1OLr97jn79G/GDy/X4l/8+evxGx/99cByBz2ng849zDGPCufj37/8r/Z9Toc9T21zN23zsHMPc4y2uUPbfEzbPD7a5tPRNrv10Da1zeGO0TZ3aJuPnZa22Tl3Ps5941/1jzR7jpyqqXz7aTabsby8HFmWxfr6ejSbzfjKV75ybOffL2SampqKhYWFmJ+fH7pMsVgc6pyjEEytrkb8D//zhSjEo+jEufhffndj4HH/959diM98tB3nOx8O3N+JQvz8hSQ+92F24F/IHxVe6DvHk3+Bf1R4Ic51PtpzjgcTs/Gtr6z26vLih5uHnvvJ+j0qnI+IOOQeJuI/vzgdlx629q3/L8+PHficDnuWB91Dtx7tiddi+uHbhz7LQc9q93m69/zkca/+k53//Z0fzw2815+/kMRnP9x6qnPvPubnL0wc2G4+eGHy0OekbWqbu2mb2uZ+9dA2d2ibu+9B24zQNrvltM299dA2d2ibu+9B24x4dm3z7z43Gb/xwf7P6bQbJg954RnVaY9SqRQPHz6M69evR7vdjm984xvRbDbjS1/60oled3NzcAP/OGWe5pzPi9/9P+fjhU89ikedc3G+8CiWk/n4i/+9P7H9+//HfHz2/K+i8OHgDhURUYjOoX8RR+z9i/DJ4/frtK/88j/Ev/4nP4jPvfdn8blHW0c695PXOWh/95jPdbZiOnt732POdz6Kz77wq32f02c+9eu5x48GP8vucfvdQ7ce09nBfxHv1OXw+9nvmP/n//2t+P/+t1r8/f/4Hwbu/9yH2VOfe/cxn/tw88B/RfhcZ+vA53TY89Q2d19D2zzs3LuP0TYPLtu9jrb5mLa5l7apbfZdR9vUNp+gbe6lbWqbu33ir7d21uJ6TkdMDeXEJhT+2sLCQqdQKHQuX77c+c53vtO3vtR3vvOdvmOf/PxxWGPqGJhXffB/5lU/m5+tOf/aprapbWqb2qa2qW1qm9qmtqltapvPlWHykBMfMTUzMxOPHj3qTd2rVqtRKBQiSZJYXFyMiJ039/3O7/xOXLx48diuWywWI0mSSNM0isVi377Z2dmnLjPsOZ9rjx5FlMuPE9run48e9R8zM+NNFBF730Tx5HM67Fl2P3/+8zv/21tSHr+JYtBz0ja1TW1T29Q2tU1tU9vUNrVNbVPbHPW38j35nEbQia8x9fbbb0ehUOhb4Hx9fT0mJydjfHw8IiK+/OUvR7lcjpmZmfjqV796bNeu1+uRJElv7admsxmrq6uxtLQUERFpmkaz2YxqtXrkMoft320U1pgCAAAAGMapW/x8a2srHjx4cCJv3ztMvV7vjW66d+9e1Gq13r5GoxG1Wi3W1taOXOYo+7sEUwAAAMBZc6LB1NWrV+OHP/xhRERUKpX4oz/6o6ev6YgTTAEAAABnzTB5yLlhTry4uBgzMzNx//79uHv3bvzoRz+KmZmZvv0vv/xynD9/Pi5cuBBf+9rX4sc//vHT3QUAAAAAI22oEVM3btyIW7du9W2r1+vx4MGDSNM0JiYmYm5uLjY2NiJN02i1WrG1tRWNRiO+8Y1vHHvlTzsjpgAAAICzZpg8ZKi38r300kt7tlWr1bh8+XK8/vrr8c1vfnPP/maz2Vso/CyGUwAAAAAMNtRUvo2NjX33TU5ODtxeKpUiTdP43ve+F9vb28PVDgAAAICRNVQwVSwW4/333x+4PU3TfcslSRIrKytx8+bNoSsIAAAAwGgaKpi6fv163L59e8/IpyzLIkmSA8tevHhx31FVAAAAAJw9Q60xFRFx69ateP311+Pq1avx6quvRkTEW2+9daSyhUJh2MsBAAAAMKKGGjHVdfv27Xj48GG88cYb8eabbx653EFrVAEAAABwtgw9YqrrypUrceXKlXj33Xfjxo0bvdFQc3NzvZFUu7377rum8gEAAADQU+h0Op3jPOHbb78dq6urvaDqwoULsba2FkmSnLnFz7e3t2N8fDy2trZibGws7+oAAAAAnLhh8pBjD6ae9Pbbb0e73Y6NjY146aWXYnp6Oi5fvnwmghrBFAAAAHDWnKpg6knvvvtutFqt2NjYiEKhEDMzM/HVr371WVbhmRFMAQAAAGfNMHnIU68x9bReeeWVeOWVV3qf19fXn3UVAAAAADgFnuqtfE9re3t7z7aLFy8+yyoAAAAAcEo8k2Bqe3s7Xn/9dW/lAwAAAKDnRKfyvf/++1Gr1aLRaESn0+m9qQ8AAAAATmTE1Pvvvx/Xrl2LqampWFpaitdeey3m5+dP4lIAAAAAPKeONZh67733eoHUyspKXLlyJdbW1uKtt96KS5cuHeelAAAAAHjOHctUvvfeey8WFhai1WpFp9OJUqkUtVrNwuYAAAAA7OtjBVPvvPNOLCwsRLvdjk6nE/Pz87GwsCCQAgAAAOBQTxVMDQqkarVajI+PH3f9AAAAABhRQwVTb775ZiwsLESaptHpdKJarcbi4qJACgAAAIChHTmY+uEPfxjlcjkmJiaiWq3GjRs3BFIAAAAAPLUjv5XvypUrcf/+/SiXyzE1NSWUAgAAAOBjOXIwFRExPT0dt2/fjtdeey1u3LgR3//+90+qXgAAAACMuKda/PzixYtx69atWF9fjxs3bsRLL70Uf/iHf3jcdQMAAABghD1VMNXVDai2trZ6AdX8/HyMjY0dV/0AAAAAGFEfK5jqGh8f7wVU3/72twVUAAAAABzqWIKprt0B1dLSUhQKhVwDqnq9HkmSRERElmVRrVaPVCYiYm1tLSIilpaWevtarVYsLS3F3NxcFIvFWF1djZmZmSiVSsdfeQAAAIARd6zBVNf4+Hh861vfioiI73znO7G5uRmdTuckLrWvbsA0Pz8fETuhUqVS6QuanrSwsBC1Wq33uVKpxNzcXKyurkbETrjVarWi2WxGsViMhYUFoRQAAADAUyp0nlFiVK/X49atW7G5ufksLhcTExOxvr7eGzEVEVEoFPYNyLIsi3K5HCsrK70y7XY7Ll26FGtra1EsFqPZbMbs7GzfOQ+yvb0d4+PjsbW1ZVojAAAAcCYMk4ece0Z1imq1+sxCqTRNI8uygQFSq9Xat9z9+/cjTdPe52KxGBE7oRUAAAAAx+tEpvLlbXe4tFuSJPuGTEmSxMOHD/u2dUOsbkAVEbG8vByTk5OxubkZa2trfVP/AAAAADi6kQym9tMNlI7q5s2bsbS01Bt5NT09HRGPg6pGo9Gb/gcAAADAcJ6LYKrZbMbdu3cPPW5xcbEXHg0yTCi1sLAQ165d6y2eHtE/cioi4urVq1GpVPadNggAAADA/p6LYKpUKg319rsnA6SuLMv23bdbs9mMqampvlCqu313PbphVJqmBwZiAAAAAOz1zBY/f5aKxWIkSTJwranZ2dkDy3bXleqGUlmW9RZTL5fLfefsrld1lLALAAAAgH4jGUxF7Ezr2/0Gvmaz2TcCKk3TqNfrfWXa7Xa02+2Ynp6ONE0jTdNoNBoxOTkZSZJEtVrtC6EajUaUSiXT+AAAAACeQqHT6XTyrsRJqdfrvSDp3r17fW/QazQaUavVYm1tLSJ2Rj9dvHhx4Fv7uo8oy7JoNBq97RsbGwe+lW97ezvGx8dja2srxsbGjuOWAAAAAE61YfKQkQ6m8iaYAgAAAM6aYfKQkZ3KBwAAAMDpJpgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAABy8ULeFThJ9Xo9kiSJiIgsy6JarR54fKvViqWlpZibm4tisRirq6sxMzMTpVLpqc8JAAAAwGAjO2KqXq9HRMT8/HzMz8/H9PR0VCqVA8tkWRatVisqlUpUKpWYmpraE0oNe04AAAAABit0Op1O3pU4CRMTE7G+vt4b3RQRUSgU4qDbbTabMTs721fm45xze3s7xsfHY2trK8bGxp7mNgAAAACeK8PkISM5YipN08iybGDA1Gq1Ts05AQAAAM6ykVxjKk3TgduTJIksyw4su7y8HJOTk7G5uRlra2tRq9U+9jkBAAAA2Gskg6n9dAOn/UxPT0dERLFYjIiIRqMR5XI5VlZWnvqcAAAAAAz2XARTzWYz7t69e+hxi4uLvXBpkMMCpG4g1XX16tWoVCoHjogSSgEAAAA8necimCqVSn1vxzvMkwFTV5Zl++6L2AnAdl+nu55UmqZPfU4AAAAABhvJxc+LxWIkSTJwXajZ2dmBZbIsi3K53FemO1KqWCw+1TkBAAAA2N9IBlMRO9P6dr8tr9lsxvz8fO9zmqZRr9d7n5MkiWq12jf6qdFoRKlU6o2cOuycAAAAABxdodPpdPKuxEmp1+u9oOnevXu9N+xF7IROtVot1tbWetuyLItGo9H7vLGx0VfmsHM+aXt7O8bHx2NrayvGxsaO5Z4AAAAATrNh8pCRDqbyJpgCAAAAzpph8pCRncoHAAAAwOkmmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgFy/kXYGTVK/XI0mSiIjIsiyq1eqBx5fL5bh27VoUi8Veua5isRitViuWlpZibm4uisVirK6uxszMTJRKpRO6AwAAAIDRNbLBVL1ej4iI+fn5iIhotVpRqVRiaWlp3zLtdjuazeae7aVSKVZWViLLsmi1WtFsNqNYLMbCwoJQCgAAAOApFTqdTifvSpyEiYmJWF9f7xv5VCgU4qDbrdfre0ZVNRqNXrjVbDZjdnZ2z2iq/Wxvb8f4+HhsbW3F2NjY0PcAAAAA8LwZJg8ZyTWm0jSNLMsGBkitVmvfck+Ofmq1WnH58uXjrh4AAAAAMaJT+dI0Hbg9SZLIsmzfcsVise8caZrG7Oxs3zHLy8sxOTkZm5ubsba2FrVa7VjqDAAAAHDWjGQwtZ9uoHQUtVptz3pU09PTEfE4wGo0GlEul2NlZeV4KwoAAABwBjwXwVSz2Yy7d+8eetzi4mIvPBrkqKFUu90euH33iKqIiKtXr0alUtl32iAAAAAA+3sugqlSqTTU2++eDJC6sizbd99uS0tLMTU1tWd7s9nsq0c3jErT9MBADAAAAIC9RnLx82KxGEmSDFxr6sk1owZptVp7RkBlWRblcrnvnN31qo4SdgEAAADQbySDqYidaX2738DXbDZjfn6+9zlN06jX6wPLpmm6J2xKkiSq1Wrf9kajEaVSyTQ+AAAAgKdQ6HQ6nbwrcVLq9XovSLp3717fG/QajUbUarVYW1vbU25qaipWVlb2TM/LsiwajUbv88bGxoFv5dve3o7x8fHY2tqKsbGxj3s7AAAAAKfeMHnISAdTeRNMAQAAAGfNMHnIyE7lAwAAAOB0E0wBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EExxoF/84hdRKBSiUCjEL37xi7yrA8+cPsBZpw+AfgD6AGedPnCyBFMAAAAA5OKFvCtwkrIsi+Xl5VhZWYnV1dUjlanX65EkSa98tVodaj8AAAAARzOyI6ba7XYsLy9HlmWxubl5pDL1ej0iIubn52N+fj6mp6ejUqkceT8AAAAAR1fodDqdvCtxkprNZty8eTMePHhw6LETExOxvr7eGxEVEVEoFKL7iA7b/6Tt7e0YHx+Pra2tGBsb+1j3kZdf/OIX8bnPfS4iIn7+85/HZz/72ZxrBM+WPsBZpw+AfgD6AGedPjC8YfKQkR0xNaw0TSPLsr7QqavVah26HwAAAIDhCKZ+LU3TgduTJIksyw7dDwAAAMBwRnrx8+MwOTkZm5ubA0dK7d4/SHeK3/b29klV78TtfhXm9vZ2fPTRRznWBp49fYCzTh8A/QD0Ac46fWB43RzkKKtHPRfBVLPZjLt37x563OLiYkxPTx/rtQ9bOP2g/R988EFERHzxi1881jrl5Qtf+ELeVYBc6QOcdfoA6AegD3DW6QPD+eCDD2J8fPzAY56LYKpUKkWpVDrRaxSLxYHbsyyLYrF46P5BvvCFL8RPf/rTePHFF6NQKBxbXQEAAABOq06nEx988MGRgrznIph6ForFYiRJEmma7gmaZmdnIyIO3f+kc+fOxW//9m+fTIUBAAAATqnDRkp1jfzi5/tNtUvTNOr1et+2xcXFvjfsNZvNmJ+fP/J+AAAAAI6u0DnKSlTPoTRNe2tTtdvtqFarMTMz05sS2Gg0olarxdraWl+5er3eGxF17969qNVqQ+0HAAAA4GhGNpjieNTr9d4bCbMsi2q1mm+F4AS1Wq1YWlqKubm5KBaLsbq62hdoR+gTjJYsy2J5eTlWVlZidXV1z/7D2rv+wPPuoD7gO4GzpDuTpPuP9ktLS3v2+z5glB3UB3wfnDxrTLGvbufsTldstVpRqVT2fFHBqMiyLFqtVjSbzSgWi7GwsLDnCydCn2A0tNvtuH//fmRZNnDa+2HtXX/geXdYH/CdwFmxsLDQNwukUqnE3NxcL6z1fcCoO6wP+D44eUZMsa+JiYlYX1/vJb8REYVCITQZRlWz2YzZ2dm+Nr+bPsEoajabcfPmzXjw4EHf9sPau/7AqNivD/hO4CzIsizK5XKsrKz02nK73Y5Lly7F2tpaFItF3weMtKP0Ad8HJ2/kFz/n6aRpGlmWDex8uxeAh7NCn+AsOay96w+cdfoAo+T+/fuRpmnvc3c93SzLfB9wJhzUBw6jDxwPU/kYaHfH3C1JkiN1UHheLS8vx+TkZGxubsba2lpvWK8+wVlyWHvXHzgrfCcw6pIkiYcPH/Zt6/4yXSwW4/79+/uW833AKDisD3T5PjhZgimG0u2MMIqmp6cj4vGXUKPR6A3t3Y8+wVnSbe/7DWXXHxglvhM4q27evBlLS0v7/l0f4fuA0fZkH/B9cPIEUwxF52KU7f5XkYiIq1evRqVSOfBfO/QJzpLD2rv+wCjxncBZtLCwENeuXest4rwf3weMqkF9wPfBybPGFAM92fm6sizbdx8875rNZt/n7r+SpGmqT3CmHNbe9QfOAt8JnDXNZjOmpqb6XnPv+4CzZFAf6G7fzffB8RNMMVCxWIwkSQbOmZ2dnc2hRnCyum/k2N3mu/8K0v0/XvoEZ8Vh7V1/YNT5TuCs6a6p0x0l0l0/yvcBZ8V+fcD3wbMhmGJfi4uLfW8SaDabhw7rhedVkiRRrVb7/mWj0WhEqVTq/auIPsEo2m+o+WHtXX9gVAzqA74TOEva7Xa02+2Ynp6ONE0jTdNoNBoxOTkZEb4PGH0H9QHfB89GodPpdPKuBKdXvV7vdcJ79+713j4AoyjLsmg0Gr3PGxsbe9q8PsGoSNM0ms1m3L17N9rtdlSr1ZiZmYlSqdQ75rD2rj/wPDusD/hO4CzIsiwuXrw4cK2c3b8m+j5gVB2lD/g+OHmCKQAAAAByYSofAAAAALkQTAEAAACQC8EUAAAAALkQTAEAAACQC8EUAAAAALkQTAEAAACQC8EUAAAAALkQTAEAAACQC8EUAAAAALkQTAEAAACQC8EUAAAAALkQTAEAAACQC8EUAAAAALkQTAEAAACQixfyrgAAAMdrYWEhLly4EBsbG5FlWZTL5Zidnc27WgAAewimAABGyKVLl+LOnTsxPT3d21Yul+Py5cuRJEl+FQMAGMBUPgCAEdFoNKJYLPaFUu12O5rNZqRpmmPNAAAGE0wBAIyI1dXVyLKsb1uSJFGtVvvCKgCA00IwBQAwIorFYrRarZiamopKpRLNZjOKxWLUarW8qwYAMFCh0+l08q4EAADHY2FhoW/qXpIksb6+bn0pAOBUMmIKAGCE1Gq1WFtbi06nE0tLS5FlWSwsLORdLQCAgQRTAADPuVarFRMTE9Fut/u2z8/Px+zsrNFSAMCpJZgCAHjOLS0tRcTOGlNPStM0KpXKs64SAMCRCKYAAJ5zc3NzsbKysmdk1MLCQlQqlYGBFQDAaWDxcwCAEdBoNGJtbS0uXLgQEREbGxsxNzcXs7OzOdcMAGB/gikAAAAAcmEqHwAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkIv/HxB6/RxI9mf9AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1200x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(4.3448e-11, dtype=torch.float64)\n",
      "tensor(2.6894e-14, dtype=torch.float64)\n"
     ]
    }
   ],
   "source": [
    "# Beta beating (exact)\n",
    "\n",
    "def bbx(Ax, Bx, Cx, nux, ax, bx, axs, bxs, mux):\n",
    "    bxf = (bxs*(8*bx**2*(Bx**2 - Ax*Cx)*(torch.cos(2*mux) - torch.cos(4*nux*torch.pi))*torch.cosh(torch.sqrt(Bx**2 - Ax*Cx))**4 + 16*bx**2*(Bx**2 - Ax*Cx)*torch.cosh(torch.sqrt(Bx**2 - Ax*Cx))**2*torch.sin(mux)**2 + 8*bx*torch.sqrt(Bx**2 - Ax*Cx)*torch.cosh(torch.sqrt(Bx**2 - Ax*Cx))**3*(-4*Cx*torch.cos(mux)*torch.cos(2*nux*torch.pi) + 2*torch.sin(mux)*(2*(bx*Bx - ax*Cx)*torch.cos(2*nux*torch.pi) + (Ax*bx**2 - 2*ax*bx*Bx + Cx + ax**2*Cx)*torch.sin(2*nux*torch.pi)))*torch.sin(mux - 2*nux*torch.pi)*torch.sinh(torch.sqrt(Bx**2 - Ax*Cx)) + 16*(Cx*torch.cos(mux) + (-(bx*Bx) + ax*Cx)*torch.sin(mux))**2*torch.sinh(torch.sqrt(Bx**2 - Ax*Cx))**2 - 8*bx*torch.sqrt(Bx**2 - Ax*Cx)*torch.cosh(torch.sqrt(Bx**2 - Ax*Cx))*(-4*Cx*torch.cos(mux)*torch.cos(2*nux*torch.pi) + 2*torch.sin(mux)*(2*(bx*Bx - ax*Cx)*torch.cos(2*nux*torch.pi) + (Ax*bx**2 - 2*ax*bx*Bx + Cx + ax**2*Cx)*torch.sin(2*nux*torch.pi)))*torch.sin(mux - 2*nux*torch.pi)*torch.sinh(torch.sqrt(Bx**2 - Ax*Cx))**3 - 16*(-(bx**2*(Bx**2 - Ax*Cx)*torch.cos(nux*torch.pi)**2*torch.sin(mux)) + torch.sin(mux)*torch.sin(nux*torch.pi)*(-2*(bx*Bx - ax*Cx)*(Ax*bx**2 - 2*ax*bx*Bx + Cx + ax**2*Cx)*torch.cos(nux*torch.pi) + bx**2*(Bx**2 - Ax*Cx)*torch.sin(nux*torch.pi)) + (bx**2*Bx**2 - 2*ax*bx*Bx*Cx + (1 + ax**2)*Cx**2)*torch.cos(mux)*torch.sin(2*nux*torch.pi))*torch.sin(mux - 2*nux*torch.pi)*torch.sinh(torch.sqrt(Bx**2 - Ax*Cx))**4 + torch.sinh(2*torch.sqrt(Bx**2 - Ax*Cx))*(-16*bx*torch.sqrt(Bx**2 - Ax*Cx)*torch.sin(mux)*(-(Cx*torch.cos(mux)) + (bx*Bx - ax*Cx)*torch.sin(mux)) - 4*(-((-2*ax*bx*Bx*Cx + (1 + ax**2)*Cx**2 + bx**2*(2*Bx**2 - Ax*Cx))*torch.cos(mux)) + (bx*Bx - ax*Cx)*(Ax*bx**2 - 2*ax*bx*Bx + Cx + ax**2*Cx)*torch.sin(mux))*torch.sin(2*nux*torch.pi)*torch.sin(mux - 2*nux*torch.pi)*torch.sinh(2*torch.sqrt(Bx**2 - Ax*Cx)))))/(8*bx*torch.sqrt(Bx**2 - Ax*Cx)*(bx*torch.sqrt(Bx**2 - Ax*Cx)*torch.cosh(torch.sqrt(Bx**2 - Ax*Cx))*torch.sin(2*nux*torch.pi) + (Cx*torch.cos(2*nux*torch.pi) + (-(bx*Bx) + ax*Cx)*torch.sin(2*nux*torch.pi))*torch.sinh(torch.sqrt(Bx**2 - Ax*Cx)))*torch.sqrt((-4*bx**2*(Bx**2 - Ax*Cx)*torch.cos(2*nux*torch.pi)**2*torch.cosh(torch.sqrt(Bx**2 - Ax*Cx))**2 - (Ax*bx**2 - 2*ax*bx*Bx + Cx + ax**2*Cx)**2*torch.sin(2*nux*torch.pi)**2*torch.sinh(torch.sqrt(Bx**2 - Ax*Cx))**2 + bx*(4*bx*(Bx**2 - Ax*Cx) - (-(Ax*bx**2) + 2*ax*bx*Bx - Cx - ax**2*Cx)*torch.sqrt(Bx**2 - Ax*Cx)*torch.sin(4*nux*torch.pi)*torch.sinh(2*torch.sqrt(Bx**2 - Ax*Cx))))/(bx**2*(Bx**2 - Ax*Cx))))\n",
    "    return (100*(bxf - bxs)/bxs).real\n",
    "\n",
    "def bby(Ay, By, Cy, nuy, ay, by, ays, bys, muy):\n",
    "    byf = (bys*(8*by**2*(By**2 - Ay*Cy)*(torch.cos(2*muy) - torch.cos(4*nuy*torch.pi))*torch.cosh(torch.sqrt(By**2 - Ay*Cy))**4 + 16*by**2*(By**2 - Ay*Cy)*torch.cosh(torch.sqrt(By**2 - Ay*Cy))**2*torch.sin(muy)**2 + 8*by*torch.sqrt(By**2 - Ay*Cy)*torch.cosh(torch.sqrt(By**2 - Ay*Cy))**3*(-4*Cy*torch.cos(muy)*torch.cos(2*nuy*torch.pi) + 2*torch.sin(muy)*(2*(by*By - ay*Cy)*torch.cos(2*nuy*torch.pi) + (Ay*by**2 - 2*ay*by*By + Cy + ay**2*Cy)*torch.sin(2*nuy*torch.pi)))*torch.sin(muy - 2*nuy*torch.pi)*torch.sinh(torch.sqrt(By**2 - Ay*Cy)) + 16*(Cy*torch.cos(muy) + (-(by*By) + ay*Cy)*torch.sin(muy))**2*torch.sinh(torch.sqrt(By**2 - Ay*Cy))**2 - 8*by*torch.sqrt(By**2 - Ay*Cy)*torch.cosh(torch.sqrt(By**2 - Ay*Cy))*(-4*Cy*torch.cos(muy)*torch.cos(2*nuy*torch.pi) + 2*torch.sin(muy)*(2*(by*By - ay*Cy)*torch.cos(2*nuy*torch.pi) + (Ay*by**2 - 2*ay*by*By + Cy + ay**2*Cy)*torch.sin(2*nuy*torch.pi)))*torch.sin(muy - 2*nuy*torch.pi)*torch.sinh(torch.sqrt(By**2 - Ay*Cy))**3 - 16*(-(by**2*(By**2 - Ay*Cy)*torch.cos(nuy*torch.pi)**2*torch.sin(muy)) + torch.sin(muy)*torch.sin(nuy*torch.pi)*(-2*(by*By - ay*Cy)*(Ay*by**2 - 2*ay*by*By + Cy + ay**2*Cy)*torch.cos(nuy*torch.pi) + by**2*(By**2 - Ay*Cy)*torch.sin(nuy*torch.pi)) + (by**2*By**2 - 2*ay*by*By*Cy + (1 + ay**2)*Cy**2)*torch.cos(muy)*torch.sin(2*nuy*torch.pi))*torch.sin(muy - 2*nuy*torch.pi)*torch.sinh(torch.sqrt(By**2 - Ay*Cy))**4 + torch.sinh(2*torch.sqrt(By**2 - Ay*Cy))*(-16*by*torch.sqrt(By**2 - Ay*Cy)*torch.sin(muy)*(-(Cy*torch.cos(muy)) + (by*By - ay*Cy)*torch.sin(muy)) - 4*(-((-2*ay*by*By*Cy + (1 + ay**2)*Cy**2 + by**2*(2*By**2 - Ay*Cy))*torch.cos(muy)) + (by*By - ay*Cy)*(Ay*by**2 - 2*ay*by*By + Cy + ay**2*Cy)*torch.sin(muy))*torch.sin(2*nuy*torch.pi)*torch.sin(muy - 2*nuy*torch.pi)*torch.sinh(2*torch.sqrt(By**2 - Ay*Cy)))))/(8*by*torch.sqrt(By**2 - Ay*Cy)*(by*torch.sqrt(By**2 - Ay*Cy)*torch.cosh(torch.sqrt(By**2 - Ay*Cy))*torch.sin(2*nuy*torch.pi) + (Cy*torch.cos(2*nuy*torch.pi) + (-(by*By) + ay*Cy)*torch.sin(2*nuy*torch.pi))*torch.sinh(torch.sqrt(By**2 - Ay*Cy)))*torch.sqrt((-4*by**2*(By**2 - Ay*Cy)*torch.cos(2*nuy*torch.pi)**2*torch.cosh(torch.sqrt(By**2 - Ay*Cy))**2 - (Ay*by**2 - 2*ay*by*By + Cy + ay**2*Cy)**2*torch.sin(2*nuy*torch.pi)**2*torch.sinh(torch.sqrt(By**2 - Ay*Cy))**2 + by*(4*by*(By**2 - Ay*Cy) - (-(Ay*by**2) + 2*ay*by*By - Cy - ay**2*Cy)*torch.sqrt(By**2 - Ay*Cy)*torch.sin(4*nuy*torch.pi)*torch.sinh(2*torch.sqrt(By**2 - Ay*Cy))))/(by**2*(By**2 - Ay*Cy))))\n",
    "    return (100*(byf - bys)/bys).real\n",
    "\n",
    "plt.figure(figsize=(12, 4))\n",
    "plt.errorbar(ring.locations().cpu().numpy(), 100*((bx - bx_id)/bx).cpu().numpy(), fmt='-', marker='x', ms=0, color='blue', alpha=0.75, label='exact')\n",
    "plt.errorbar(ring.locations().cpu().numpy(), -bbx(Ax, Bx, Cx, nux, axi, bxi, ax, bx, mux), fmt=' ', marker='x', ms=5, color='red', alpha=0.75, label='approximate')\n",
    "\n",
    "plt.gca().tick_params(axis='x', length=6, width=1.5, direction='in', labelsize=12, bottom=True, top=False, labelbottom=True, labeltop=False)\n",
    "plt.gca().tick_params(axis='y', length=0, width=0, labelsize=12)\n",
    "plt.gca().set_xlabel(r'$s$', fontsize=18)\n",
    "plt.gca().set_ylabel(r'$\\Delta \\beta_x / \\beta_x$', fontsize=18)\n",
    "plt.legend(loc='upper left', frameon=False, fontsize=14, ncol=2)\n",
    "plt.ylim(-1, 1)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(100*(((bx - bx_id)/bx)**2).mean().sqrt())\n",
    "print((bbx(Ax, Bx, Cx, nux, axi, bxi, ax, bx, mux)**2).mean().sqrt())\n",
    "\n",
    "plt.figure(figsize=(12, 4))\n",
    "plt.errorbar(ring.locations().cpu().numpy(), 100*((by - by_id)/by).cpu().numpy(), fmt='-', marker='x', ms=0, color='blue', alpha=0.75, label='exact')\n",
    "plt.errorbar(ring.locations().cpu().numpy(), -bby(Ay, By, Cy, nuy, ayi, byi, ay, by, muy), fmt=' ', marker='x', ms=5, color='red', alpha=0.75, label='approximate')\n",
    "plt.gca().tick_params(axis='x', length=6, width=1.5, direction='in', labelsize=12, bottom=True, top=False, labelbottom=True, labeltop=False)\n",
    "plt.gca().tick_params(axis='y', length=0, width=0, labelsize=12)\n",
    "plt.gca().set_xlabel(r'$s$', fontsize=18)\n",
    "plt.gca().set_ylabel(r'$\\Delta \\beta_y / \\beta_y$', fontsize=18)\n",
    "plt.legend(loc='upper left', frameon=False, fontsize=14, ncol=2)\n",
    "plt.ylim(-1, 1)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(100*(((by - by_id)/by)**2).mean().sqrt())\n",
    "print((bby(Ay, By, Cy, nuy, ayi, byi, ay, by, muy)**2).mean().sqrt())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "03696dde-8e7c-42ea-aba3-d89cc79c59dd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABKYAAAGGCAYAAABBiol3AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAOcVJREFUeJzt3V9sXFl+H/hfieru6X9kkWpj/uyMMyp2Z/y0sSgxgNdjwJgmbWRjYIERSzLgzSRet1g9BvbNQzaDfVxYzep52wc3S404O44fRNYMkIckSLN6HMRjBxhJVM8uEMT28KphYz02YlKX7D/zxy3VPtRUiSVWkSyK5KWKnw8gSHX/nnvrHN7Lr849N1ev1+sBAAAAAEfsVNYFAAAAAOBkEkwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkInTWRfgsKRpGouLi7G0tBTLy8t7WqdcLkc+n2+tPzMz09N8AAAAAPauL3tMraysxOLiYqRpGuvr63tap1wuR0TE9PR0TE9Px9jYWJRKpT3PBwAAAKA3uXq9Xs+6EIelWq3G1atX49atW7suOzw8HHfu3Gn1iIqIyOVy0Tw9u80HAAAAoDd92WOqV0mSRJqmbaFTU61W23U+AAAAAL3r2zGmepEkScfp+Xw+0jTddX4n9+/fj7/+67+O559/PnK53EEVFQAAAOBYq9fr8f7778dnPvOZOHVq5z5RgqkdjIyMxPr6eseeUlvnd/LXf/3X8bnPfe4QSwcAAABwfP3VX/1VfPazn91xGcHUDnYbOH2n+c8//3xENL6EwcHBAy0XAAAAwHG1ubkZn/vc51rZyE4EUxFRKBQ6Tk/TNAqFwq7zO2k+vjc4OCiYAgAAAE6cvQxtZPDzaART+Xy+41hSExMTu84HAAAAoHd9HUx1e9QuSZIol8tt0+bm5tresFetVmN6enrP8wEAAADoTa5er9ezLsRBS5IkqtVqXL9+PVZWVmJmZibGx8djamoqIiIqlUrMz8/H6upq23rlcrn1aN6NGzdifn6+p/lbbW5uxtDQUGxsbHiUDwAAADgxeslE+jKYOg4EUwAAAMBJ1Esm0teP8gEAAABwfAmmAAAAAMiEYAoAAACATAimAAAAAMiEYAoAAACATAimAAAAAMiEYAoAAACATAimAAAAAMiEYAoAAACATAimAAAAAMiEYArgCFQqlSgWi1kXAwAA4FgRTAEcgVu3bkW1Ws26GAAAAMeKYAp+qlqtRpqmWReDPrWwsBD1ev1I96lOAwAAx51gCn6qWCzGzZs3sy4GHBh1GgAAOO4EU3CS/et/HfGNb7RP+8Y3GtMBAADgkAmmjrF6PeJHPzr+fx716aRqtRrnz5+PXC4X58+fj1qtFhERs7OzkcvlIkmStuUnJydjdHS09TlJkigWizE8PNzaxsrKyrb9rKysxOTkZAwPD8fw8HBMTk7GyspKFIvFyOVyrW3ncrnW57536lTE7//+g3DqG99ofD6V/Y+GvXyvlUql9T02v9vR0dFWHepluXK5HJOTkxERUSqVYnh4uPUY3MPrlcvl1np7raflcjmGh4fbPjeP6eH6n6Zp69hHR0ejUqn0dG72Uqe7tTsAAICjdDrrAtDdj38c8U/+Sdal2N1/+A8Rn/jE/tatVCpRKpViYWEhLly4ELVaLSYnJ2N1dTXm5+dbgcDq6mpr+Vqt1voc0fgFe2RkJJaWlqJQKMTs7Gy8/PLLcefOncjn8xERrV/+p6amYmlpKSIilpeX4/r163Ht2rWYm5uL8+fPx9LSUkxMTDzS+XisfOUrjb9///cj/uAPIj7+OOI3f/PB9Azt5XtN0zRqtVpcuXIl5ubmolQqxezsbKvOFAqFnpZLkiTOnz8faZrG3Nxc5PP5Vp2cmZmJ+fn5SJIkZmdn48aNG7G0tLTnevqwtbW1WFlZiStXrsT8/HxENAKxYrEYhUKhVc75+fkolUoxMTHRKudu52a3Or1Tu2vuAwAA4CgIpshMmqatX46np6cjImJsbCzW1tZifn4+FhYWYnl5OUZHR6NYLLZ+UW/+Mt40MzPTtt1r167F8PBwLC4utrZbLBbbQqmIaPtlvbm9fD7fCj1OjK985UEodfr0sQilIvb2vTbdunWr9e+JiYkYHh5u1aFelkuSJCYmJtqWK5VKrVAqolFHx8bGYnR0NKrVakxNTe2pnnYzPz/fqouzs7OtEGpqaioiIkZGRlo9oprb2+3c5PP5rnV6L+0OAADgqAimjrGnnmr0Rjrunnpqf+s1B2UulUpRKpXa5o2NjbX+3fylv1qtxszMTOsX9m6av4Q3e6skSRJJkviFu5tvfONBKPXxx43PxySc2urh73Wn5aampnZ9NK3bclvrycrKSiRJsq1+FgqFGBsbi+vXr7fqY6/1tOnChQtt242I1iOFEQ/awvr6+o7HErH7uYnYe7sDAAA4CtkPJENXuVzjEbnj/udRh2O6e/du1Ov1tj9be6wUCoVWz465ubmO26jValEsFmN0dLRtHJ+IaI394xGlDppjSv3mb0YsLzf+3jrmVMZ2+l53UigUdgxydlpuaz1phjgjIyMd1906rtRe6mknnXroddrfw/Z7bpp2a3cAAABHQTBFZraO67OTWq0Wi4uLMTU1FS+//PK2+cViMYrFYkxOTsby8nLcvXt3X/s5ke7fbx9T6itfaXy+fz/bcsXu3+tOkiTZUxD58HIPh0Q71Z2H192tnnbS7bHR3R4nfZRzoz0AAADHiWCKzDQfh7p69eq2ec23oTXfPnbt2rXW+FBbHz9K0zSq1Wpcu3YtpqenO4YRhUIhCoVCx0f5mvvZOpj2ifIv/sX2x/a+8pXG9Azt5XvtJkmSqFaruw5iv5flJiYmIp/Pb6s7KysrsbKyEpcvX25ta6d6epD2em661em9tDsAAICjYowpMnXt2rU4f/58a7yb9fX1WFpaav1dLBbbBoJeWlqK0dHROH/+fGuQ53w+H1evXo18Ph8jIyMdf+FeWFiIycnJ1lvP0jSN5eXluHnzZuvxpXw+H9evX2+FWPPz8ydvIPRjYq/fa9Pk5GTMzs5GmqZx5cqVyOfzHR+n2+tyW127di2KxWJENHoqNd/Kt7Ve7lZPD1Iv56Zbnd6t3QEAABwVPabI1NjYWKyurkaSJPHyyy+3AoBr165FqVSKJEni2rVrreWbv2CXSqVYWVlpLdvssXLlypWYnJyMiYmJGB0dba3XfNPa1uUe/iV8bm4uqtVqqwxkay/fa1OpVIrZ2dkoFotx4cKFuHXrVsdQca/LbTU1NdWqO5OTkzE/Px9zc3OxvLzc2uZe6ulB2uu56Vand2p3AAAARylXr9frWReiH21ubsbQ0FBsbGzE4OBg1sWBvlQul2N2djZ2+zG21+UAAAB4dL1kInpMAQAAAJAJwRQAAAAAmejrwc/L5XLbm6lmZmZ2XL5YLMbly5ejUCh0fG18rVZrDaJdKBRieXk5xsfHWwMeAwAAALB3fTvGVLlcjohohVG1Wi2Wlpa2vfZ9q9HR0UiSZNv0qampWFpaimq1GleuXIk0TaNQKMTs7GzXN24ZYwoAAAA4iXrJRPo2mBoeHo47d+609XzK5XI7Dn5cLpe39aqqVCqt8KlarcbExMSub/GKEEwBAAAAJ9OJH/w8SZJI07RjgFSr1bqu9/AjebVaLS5cuHDQxQMAAAAg+nSMqU6P40VE5PP5SNO063qFQqFtG0mSxMTERNsyi4uLMTIyEuvr67G6uhrz8/MHUmYAAACAk6Yvg6lumoHSXszPz28bj2psbCwiHgRYlUolisViLC0tHWxBAQAAAE6AvnyUr5u9hlIrKysdpxcKhbZeVZcuXYpqtbpjLywAAAAAOuvLYGpreLRV8216u1lYWIjR0dFt06vVatvn5hhW3R4dBAAAAKC7vg2m8vl8x8Do4TGjOqnVatsGTk/TNIrFYts2mz2l9hJ2AQAAANCuL4OpiIi5ubm2N/BVq9WYnp5ufU6SJMrlcsd1kyTZFjbl8/mYmZlpm16pVGJqaqrj2/8AAAAA2FmuXq/Xsy7EYSmXy60g6caNG21v0KtUKjE/Px+rq6vb1hsdHY2lpaXWYOdNaZpGpVJpfV5bW+v6Vr7Nzc0YGhqKjY2NGBwcPIjDAQAAADj2eslE+jqYypJgCgAAADiJeslE+vZRPgAAAACON8EUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJkQTAEAAACQCcEUAAAAAJk4nXUBDlO5XI58Ph8REWmaxszMzI7L12q1WFhYiMnJySgUCrG8vBzj4+MxNTW1720CAAAA0Fnf9pgql8sRETE9PR3T09MxNjYWpVJpx3XSNI1arRalUilKpVKMjo5uC6V63SYAAAAAneXq9Xo960IchuHh4bhz506rd1NERC6Xi50Ot1qtxsTERNs6+93m5uZmDA0NxcbGRgwODu73MAAAAAAeK71kIn3ZYypJkkjTtGPAVKvVjs02AQAAAE6yvhxjKkmSjtPz+XykabrjuouLizEyMhLr6+uxuroa8/Pzj7xNAAAAALbry2Cqm2bg1M3Y2FhERBQKhYiIqFQqUSwWY2lpad/bBAAAAKCzExVM7RYgNQOppkuXLkWpVNqxR5RQCgAAAGB/+nKMqYcDpqY0TbvOi2gMfr5VczypJEn2vU0AAAAAOuvbYCqfz3ccF2piYqLjOmmaRrFYbFun2VOqUCjsa5sAAAAAdNeXwVRExNzcXNvb8qrVakxPT7c+J0kS5XK59Tmfz8fMzExb76dKpRJTU1OtnlO7bRMAAACAvcvV6/V61oU4LOVyuRU03bhxo/WGvYhG6DQ/Px+rq6utaWmaRqVSaX1eW1trW2e3bW61ubkZQ0NDsbGxEYODgwd2TAAAAADHWS+ZSF8HU1kSTAEAAAAnUS+ZSN8+ygcAAADA8SaYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMiGYAgAAACATgikAAAAAMnE66wIcpnK5HPl8PiIi0jSNmZmZPa0TEbG6uhoREQsLC615tVotFhYWYnJyMgqFQiwvL8f4+HhMTU0dfOEBAAAA+lzfBlPNgGl6ejoiGqFSqVRqC5oeNjs7G/Pz863PpVIpJicnY3l5OSIa4VatVotqtRqFQiFmZ2eFUgAAAAD7lKvX6/WsC3EYhoeH486dO60eUxERuVwuuh1umqZRLBZjaWmptc7KykqcP38+VldXo1AoRLVajYmJibZtdrO5uRlDQ0OxsbERg4ODB3BEAAAAAMdfL5lIX44xlSRJpGnaMUCq1Wpd17t582YkSdL6XCgUIqIRWgEAAABwsPryUb6t4dJW+Xy+a8iUz+fj7t27bdOaIVYzoIqIWFxcjJGRkVhfX4/V1dW2R/8AAAAA2Lu+DKa6aQZKe3X16tVYWFho9bwaGxuLiAdBVaVSaT3+BwAAAEBv+vJRvm56CaVmZ2fj8uXLrcHTIxqB1NbeU5cuXYpqtepRPwAAAIB9OJJg6o033oiXXnop3n333da0rf8+aFvDo63SNO06b6tqtRqjo6MxMzOzbfpWzZ5U3R4dBAAAAKC7Iwmm8vl8vP766zE0NNSaVq/XY25uLjY3Nw98f4VCIfL5fMfAaGJiYsd1m+NKNXtKpWnaGky9WCy2bbPZU2ovYRcAAAAA7Y4kmBoZGYnz58/H2bNnW9POnTsXV69ejcXFxUPZ59zcXNsb+KrVattjeUmSRLlcbltnZWUlVlZWYmxsLJIkiSRJolKpxMjISOTz+ZiZmWkLoSqVSkxNTXV8+x8AAAAAO8vV6/X6UezoV37lV+Lu3btRKpXi0qVLMTg4GBERX//61+N3fud3DmWf5XK5FSTduHGj7Q16lUol5ufnY3V1NSIavZ/Onj3bcbyo5ilK0zQqlUpr+traWte38m1ubsbQ0FBsbGy0jhUAAACg3/WSiRxJMPXVr3619Ua75eXl1hhOERGlUunQgqksCaYAAACAk6iXTOT0URRoYmIiLl68GBERV65ciYjGWE63b99ue7wOAAAAgJPjSMaYitj+Fr6JiYn42te+1jYOFAAAAAAnx5EEUxcvXoy1tbX4+te/3pp2+/btGB8f7/jmPAAAAAD635ENft7JtWvX4sKFC3Hu3LmsinBojDEFAAAAnES9ZCL76jF16dKlGBgYiIGBgfjt3/7tfRUyojHeVD+GUgAAAADsrudgam5uLsbHx+PmzZtx/fr1+O53vxvj4+Nt81966aUYGBiIM2fOxK//+q/H9773vQMtNAAAAACPv54f5Xvttdfi9ddfb5tWLpfj1q1bkSRJDA8Px+TkZKytrUWSJFGr1WJjYyMqlUr81m/91oEW/jjzKB8AAABwEvWSiZzudeMvvPDCtmkzMzNx4cKFePXVV+OVV17ZNr9arcb09HRExIkKpwAAAADorudH+dbW1rrOGxkZ6Th9amoqkiSJ3/u934vNzc1edwkAAABAH+o5mCoUCvHee+91nJ4kSdf18vl8LC0txdWrV3vdJQAAAAB9qOdg6sqVK/Hmm29u6/mUpmnk8/kd1z179mzXXlUAAAAAnCw9jzEVEfH666/Hq6++GpcuXYovfelLERHx9ttv72ndXC63n10CAAAA0Gd67jHV9Oabb8bdu3fjq1/9anzrW9/a83o7jVEFAAAAwMmxrx5TTRcvXoyLFy/G7du347XXXmv1hpqcnGz1pNrq9u3bHuUDAAAAICIicvV6vX7QG33nnXdieXm5FVSdOXMmVldXI5/Pn5jBzzc3N2NoaCg2NjZicHAw6+IAAAAAHIleMpFDCaYe9s4778TKykqsra3FCy+8EGNjY3HhwoW+DmwEUwAAAMBJdOyCqYfdvn07arVarK2tRS6Xi/Hx8fjyl7981MU4VIIpAAAA4CTqJRN5pDGm9uvcuXNx7ty51uc7d+5kUQwAAAAAMrTvt/IdpLNnz2ZdBAAAAACO2J6DqXfeeScuX74cf/RHf3SY5QEAAADghNhzMPXyyy/HpUuX4urVq3HmzJmYm5uLd9999xCLBgAAAEA/6+lRvosXL8bbb78dSZJEoVCIV155JV566aX4+te/Hpubm4dVRgAAAAD60L7GmBoaGoorV67EzZs34+23346/+7u/i7GxsfjVX/3VeOuttw66jAAAAAD0oUce/Pzs2bPx+uuvx/e///14/fXX4+bNmzEyMhKXL1+Ob3/72wdRRgAAAAD60IG+le/cuXPx5ptvxvr6ely6dCnefPPNOHPmTHz1q181HhUAAAAAbQ40mNrq4sWLsbi4GEmSxNjYWNt4VO+9995h7RYAAACAx8ShBVNND49HVa/XY2JiIsbHx+Ott94yaDoAAADACXXowdRWZ8+eja997Wvx/e9/PyqVSty8eTM+//nPG48KAAAA4ATK1ev1etaFeOedd6JWq8XVq1cPdLvlcjny+XxERKRpGjMzM4+8zl63ubm5GUNDQ7GxsRGDg4P7PgYAAACAx0kvmcixCKYOQ7lcjohoBUe1Wi2WlpZiYWFh3+v0sk3BFAAAAHASCaYiYnh4OO7cudPq3RQRkcvlYqfD3W2dXrb52AdTr7wScepURKXyYNr0dMT9+xFvvfVgme9/P+IHP4jI5SI++cmI73434u//PuLJJyM+97mIzc2I//7fI+7dixgYiPiZn4lI04gf/7ix/VOnIk6fjvj448Z6uVzE889H/OQnjc+nfvq06dNPR3z4YWM7TzwRUShE/NIvRfzxH0f8zd9E1OuNdZ58srH8M89EfOELjfL88IeNac39/93fNfb3xBON7T7zTKOM9+839v3BB41yPPvsg22trbWX7+mnG/9eW2usNzra+FyvR3z60xEvvth+nnY7l83lvvOdiI2NiKGhxvn80z9tbPO55yI+8YnGvB//uDHtiScizpxplL15nu7fb0z/0Y8elPWjjxr/bmqey4jG9CeeiPjH/7jx+W/+pvFnc7Oxv8HBxjK/9msRf/iHjX3fv99YdnCwse16/cH+n3mmMe3jj9v3/fTTje/hmWcayzWXiWgs15y+ttb4u1BoHOtLL0X83M9tr3d7rZv1esTf/m1j2g9+0Phun3oq4md/tvH5hz9sL0dExPvvP6ibzfPVHAuvecydzuf9+w/OebEY8Rd/0TiXf/mXjfP23HPb6+ZPftKo01u/q+b5PH268d00v6tnn31Qtmbd/PSnG8fZbF/PPvvgXH/4YWPZX/iFxjlo1qsvfnF7vdvL+fzOd6LND37Q2MeTTza2m6aN42i2ty98IWJ1tTHt/v0H7WvruWwe//37jTI0y37/fuM8PPXU9rr54x83Pn/8ccQLLzTq5h//8YPzsJdtP/FEY/1crnHe7t3rXjc//enGuWu2s0KhMe+Tn2ys32tbby7z3/5bo44MDTXO5UcfNco1NNRYrvlzqlk3nnmm8T12qpsPt/VOx3z/fqNN/8ZvNNb9i79o/HxpttWH62bzZ/JO9X7r+fzJTxrn8umnG+sNDTXa8717D+r0Cy801ml+N2fONJZ71Lq59TpUrzfaXPNn/OBgY58//OGDn+H/8B8+qL9br01/+7eN9T/xiYh8/sF3cPp0o+wffdRog/V64zgfrpvNe4HmtejTn4741Kfar0Of+MSD77X5M+OFFxrnqrmvgYEH/753r/F38zw1t/OpTzWOrXmNjYj4xV98cAyPch36d/+u8e8vfKFxnpLkwXXok59sfG4ew1NPNc5VswzNujkwsP061O1aMTDQ+K7Gxxv7/M//ufEd/uhHD7bfbBfN77b5XT71VPv53K2t53IPzuXW7bz0UuM7bNbNp5+OOHcu4s/+rLHff/pPH71uHuQ9UvN8PvlkY97TTz+om5/4RMQ/+Ae93yM17ykervcR7pHcI/XHPdITTzR+nkRsv6Y3HdQ90upqo6w/8zONz2naaCPNsrz44uHcI73/fqOdP/lk41z+2q9FLC4+uAZu/Z6aPx/v32+cs63nsvk99vs9UvOYm+fzN36jsd96PeLGjcbPiuN6j/SY6CkTqR+Bcrlcf/HFF+u3b99uTdv674O2urpa73RoEVFfXl7e1zq9bnNjY6MeEfWNjY19HMExcOVKvZ7PN/7u9Lk57Zln6vUnnqjXc7nGn0bzOpo/AwMP/p3LtX/O4k8u1zgXzzyz/Tztdi6b0598sn17D2//sMuf5fnrVJ5Pfapzvdtr3Xz66WzqRZb1cad68+STnevdXs7nk09m084fPoaHz2sW5RkYaNSt/bT15rRPfSq7czk0lE3d3OnPo9TNrdehLOvmE080/hxlPXy4LAdxHWput1N7P8xzfByu452OdWDgYOqmeyT3SAddHvdIB1eeXK4/7pGOw5/H+R4p4sF3mfV5fLhMD9fNx0QvmciR9Ji6du1ajIyMxNjYWJw9ezYiIm7fvh2Li4sxNzd34D2KarVaTE5OxsOHNjw8HNeuXYupqame18nn8z1t87HvMRUR91+Zjh//m6W4Vz8VubgfNz5fjD/4YqVtmX/2nen4hdU/jCfufRS5LtuJiGietZ2W2Yt6h23czw3EH7/0v8X3P/nF+Gd/+mo8ee+HB7btre7lBuJU/d4Oy+TiJwNPx38Z/Y2O52n8vaWoR/dz2VzuF//i/46B+k+67me3cu5Vt+38ZODp+IP/6c24eHM28j/8m31t+36cilNxf5dlcpGLetdjuZc7HT9+4rmu9W638/mgbv4wctH9x9xe6+Zu573b/B8MfSH+/f/4Lw+1bjbOYvdzGRHx8akn409e/Odd691ezududXMvZd2L3ermi3/7nfilP38rTu3wvT5K+Xarm/XIxd8/QltvLvPUT96Pgbi3Yzl2O8a9HE+3Op4+/an45oX5HevmXup9p233ssxB1M29XYdyO/4seBTNuhkR+27rezmX93IDkavf36FeHMx16It//q92rJsHdR3qdq1oHMP/Gr+w+m8e6edmxKP97GzeXxx+3dy9rHvhHsk9UoR7pO3zj8c9UuNcPvp1aLe6+T//P78bn974s57W3ev8iP68R+q0zHG5R6qffjJO/eY/b+9p9hg5do/yffOb34zz58/H5z//+W3z3nrrrXjllVcOdH/dQqbR0dGYnZ2N6enpntcpFAo9bbMfgqn/+l8jPvuPzkQu7kc9TsX/8sW1jsv92++ciWfubcZA/eOO8+/lTsdHA4Px3Md3d/xhfC93+qc3NZ2XqUcu7ucGtu3n/dMjrbK98b3JOH+31nHb3crXnB8RXfdfj1x8cDofEbl4/uP1rtv4aGBwx/O027lsLrfTuXr/9Mgjn8udjvfW8ER87R8tR0TEt//T9h+Te/mePjw9FM9+vBERscP5HI6I6Ho+Pzid37XePWrdbJZjL+czInbczv3WTXm9bfrLv9y4+dypbu6l3p+qNy7ODy/XnP/RwGDXc9k8zt3q3V7O507nqtt52KqX492pbv7b75zZdryPsu2ty+xWNw+irTeXee7jtGs5VoZfjrG77+xYNz84nY9nP97Y9Zgjtrf3L/1y49871c2I3et9p21vXeaD08Nd2+FB1s3d23o+nvs43bWtb93Gw2FWtzq2tW7u5zq09Vzu9vPqz58f67j95j4O6jq028+Tg7imd7tWNM/no/7cjNi9bkbU4/mP73bcztb7i04Oom4e1D1St/a6n3ukTvU+wj1ShHukTtwjPVhmt7r5zL3NI7lHeuc/ndo2f7fvabd23tx/P94jdao/x+Ue6ePnh+PJze5187jrJRM5fRQFunjxYvzKr/xK3L17N0qlUly6dKlVsDRNj6IIERGxvt65AT3KOvvZ5uPixTem49Qn7se9+qkYyN2Pxfx0/OX/0Z7W/uz/OR3PDvwoch93D30G6h/Hcx+v75pY7xQcRTR+WHRa5rn6Rizmp+ODn/9i/Ox/+ZN9bXsv+36uvvlgDJuO27gXz57+Udfz9MwnfvpM8/3O57K53HPxwY4/ZLtdCNrLsv/jPffRn8S/+uVvxP/wf83ua9u5qMez9zZ3PIZc1OO5e3cbHVS7eLb+Qdx/+rl9n88HdbP7d5aL+oGcz251Mxf1+Pd3fi7+5jf/5b7rZrdtt82Pe/HcDseRi3o8l/tgx3q3l/O5W93crawRj3a8zbr53LvfadSfA9z21mV2q5sD8WhtvbnMqR++v2M5xtKdb7gioutN29btdDvm//j/fjr+v/99/tDqZnOZnX7+H1Td3O061CjHzqFUxPZjfnj53epmROzrfO75XNY3Yix9p+syB3Ydur+xczkO6Jre7Vpx7qM/icX8dIx8dMh1c5e23ry/OMy6eVD3SN3m7+ceaa/1fuvy7pH2tm33SO3TT9I90vMfr/90HKId7i8O4B7pU7//ux3L8aj3Z839u0dqOKp7pCd+9EFjLK7HtMdUTw7tgcItXn311XqlUqlXKpV6sVis53K5+osvvlh/8cUX62+88caB788YUwfA+Am9/zF+wsGXx/gJj/49GmPq4P88zuMnGGPq8OqmMaYe/VxmfR3vdKzGmDq48+ke6WDL4x7p4MpjjKmD+/M43yNFGGPqgPWSiRxJj6mJiYm4ePFiRERcuXIlIhqPzt2+fbvjY3WPqlAoRD6fjyRJovDQ2xYmJib2vU6v23ys3b/feGtGM51t/t1820jz3+Pjx+utfFvfEpPlW/kePk+7ncvm57NnvXEmovMbZ/ZTN+v1k/XGmb28la9TvdvL+fzp+IAtWb6Vr/kmmePwVr791M1isfHGmQhv5dv6xpn91s2t16F6Pfu38g0MHI+38u33OtR8i5S38j14K99B1M3j8Fa+ne6RjvKtfO6R3CMd17fyHdQ9krfyPR73SI/TW/kerpt96MjGmBodHY2f//mf3zbvW9/6Vnz5y18+8H2Wy+XI5/Ot4Ktarcby8nIsLCxERESSJFGtVmNmZmbP6+w2f6t+GGMKAAAAoFe9ZCKnjqJAFy9ejLW1tfj617/emnb79u0YHx+PJEkOZZ8zMzORpmlUq9WoVqtx48aNtgCpVqttC5R2W2e3+QAAAADs3ZH0mOrm2rVrceHChTh37lxWRTg0ekwBAAAAJ9Gh95i6dOlSDAwMxMDAQPz2b//2vgoZ0Rhvqh9DKQAAAAB213MwNTc3F+Pj43Hz5s24fv16fPe7343x8fG2+S+99FIMDAzEmTNn4td//dfje9/73oEWGgAAAIDHX8+P8r322mvx+uuvt00rl8tx69atSJIkhoeHY3JyMtbW1iJJkqjVarGxsRGVSiV+67d+60ALf5x5lA8AAAA4iXrJRE73uvEXmq+O3WJmZiYuXLgQr776arzyyivb5ler1dab7E5SOAUAAABAdz0/yre2ttZ13sjISMfpU1NTkSRJ/N7v/V5sbm72uksAAAAA+lDPwVShUIj33nuv4/QkSbqul8/nY2lpKa5evdrrLgEAAADoQz0HU1euXIk333xzW8+nNE0jn8/vuO7Zs2e79qoCAAAA4GTpeYypiIjXX389Xn311bh06VJ86UtfioiIt99+e0/r5nK5/ewSAAAAgD7Tc4+ppjfffDPu3r0bX/3qV+Nb3/rWntfbaYwqAAAAAE6OffWYarp48WJcvHgxbt++Ha+99lqrN9Tk5GSrJ9VWt2/f9igfAAAAABERkavX6/WD3ug777wTy8vLraDqzJkzsbq6Gvl8/sQMfr65uRlDQ0OxsbERg4ODWRcHAAAA4Ej0kokcSjD1sHfeeSdWVlZibW0tXnjhhRgbG4sLFy70dWAjmAIAAABOomMXTD3s9u3bUavVYm1tLXK5XIyPj8eXv/zloy7GoRJMAQAAACdRL5nII40xtV/nzp2Lc+fOtT7fuXMni2IAAAAAkKF9v5VvvzY3N7dNO3v27FEXAwAAAICMHVkwtbm5Ga+++qq38gEAAAAQEUfwKN97770X8/PzUalUol6vt97UBwAAAMDJdmg9pt577724fPlyjI6OxsLCQrz88ssxPT19WLsDAAAA4DFz4MHUu+++2wqklpaW4uLFi7G6uhpvv/12nD9//qB3BwAAAMBj6sAe5Xv33XdjdnY2arVa1Ov1mJqaivn5eQObAwAAANDRIwdT3/72t2N2djZWVlaiXq/H9PR0zM7OCqQAAAAA2NG+g6lOgdT8/HwMDQ0dZPkAAAAA6FM9B1Pf+ta3YnZ2NpIkiXq9HjMzMzE3NyeQAgAAAKAnPQVT3/zmN6NYLMbw8HDMzMzEa6+9JpACAAAAYF96eivfxYsX4+bNm1EsFmN0dFQoBQAAAMC+9RRMRUSMjY3Fm2++GS+//HK89tpr8dZbbx1GuQAAAADoc/se/Pzs2bPx+uuvx507d+K1116LF154IX7nd37nIMsGAAAAQB/bdzDV1AyoNjY2WgHV9PR0DA4OHkT5AAAAAOhTjxxMNQ0NDbUCqt/93d/NPKAql8uRz+cjIiJN05iZmdnTOhERq6urERGxsLDQmler1WJhYSEmJyejUCjE8vJyjI+Px9TU1MEXHgAAAOAEOLBgqmlrQLWwsBC5XO7IA6pmwDQ9PR0RjVCpVCq1BU0Pm52djfn5+dbnUqkUk5OTsby8HBGNcKtWq0W1Wo1CoRCzs7NCKQAAAIBHkKvX6/XD3skbb7wR6+vrUa/X44033oh79+4d6v6Gh4fjzp07rR5TERG5XC66HWqaplEsFmNpaam1zsrKSpw/fz5WV1ejUChEtVqNiYmJtm3uZHNzM4aGhmJjY8NjjQAAAMCJ0Usm0vNb+fbja1/7Wly9ejVGRkZiaGjoUPeVJEmkadoxQKrVal3Xu3nzZiRJ0vpcKBQiohFaAQAAAHDwDvxRvp3MzMzsaaynR7E1XNoqn893DZny+XzcvXu3bVozxGoGVBERi4uLMTIyEuvr67G6utr26B8AAAAAvTnSYCpLzUBpr65evRoLCwutnldjY2MR8SCoqlQqrcf/AAAAAOjdsQ+mqtVqXL9+fdfl5ubmWuFRJ72EUrOzs3H58uXW4OkR7T2nIiIuXboUpVKp62ODAAAAAOzs2AdTU1NTPb397uEAqSlN067ztqpWqzE6OtoWSjWnby1HM4xKkmTHQAwAAACAzo5k8POjVCgUIp/PdxxramJiYsd1m+NKNUOpNE1bg6kXi8W2bTbHq9pL2AUAAADAdn0XTEU0Huvb+ga+arXa1gMqSZIol8tt66ysrMTKykqMjY1FkiSRJElUKpUYGRmJfD4fMzMzbSFUpVKJqakpj/EBAAAA7FOuXq/Xsy7EYSiXy60g6caNG21v0KtUKjE/Px+rq6sR0ej9dPbs2Y5v7WuenjRNo1KptKavra3t+Fa+zc3NGBoaio2NjRgcHDyIQwIAAAA49nrJRPo2mMqaYAoAAAA4iXrJRPryUT4AAAAAjj/BFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkAnBFAAAAACZEEwBAAAAkInTWRfgsJTL5cjn8xERkaZpzMzM7Lh8rVaLhYWFmJycjEKhEMvLyzE+Ph5TU1P73iYAAAAA3fVlj6lyuRwREdPT0zE9PR1jY2NRKpV2XCdN06jValEqlaJUKsXo6Oi2UKrXbQIAAADQXa5er9ezLsRBGx4ejjt37rR6N0VE5HK52OlQq9VqTExMtK3zKNvc3NyMoaGh2NjYiMHBwf0cBgAAAMBjp5dMpO96TCVJEmmadgyYarXasdkmAAAAwEnXd2NMJUnScXo+n480TXdcd3FxMUZGRmJ9fT1WV1djfn7+kbcJAAAAQGd9F0x10wycuhkbG4uIiEKhEBERlUolisViLC0t7XubAAAAAHR37IOparUa169f33W5ubm5VrjUyW4BUjOQarp06VKUSqUde0QJpQAAAAD279gHU1NTU21vx9vNwwFTU5qmXedFNAKwrftpjieVJMm+twkAAABAd303+HmhUIh8Pt9xXKiJiYmO66RpGsVisW2dZk+pQqGwr20CAAAAsLO+C6YiGo/1bX1bXrVajenp6dbnJEmiXC63Pufz+ZiZmWnr/VSpVGJqaqrVc2q3bQIAAADQm1y9Xq9nXYjDUC6XW0HTjRs3Wm/Yi2iETvPz87G6utqalqZpVCqV1ue1tbW2dXbb5sM2NzdjaGgoNjY2YnBw8ECOCQAAAOC46yUT6dtgKmuCKQAAAOAk6iUT6ctH+QAAAAA4/gRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGRCMAUAAABAJgRTAAAAAGTidNYFOCzlcjny+XxERKRpGjMzMzsuXywW4/Lly1EoFFrrNRUKhajVarGwsBCTk5NRKBRieXk5xsfHY2pq6pCOAAAAAKC/9WUwVS6XIyJieno6IiJqtVqUSqVYWFjous7KykpUq9Vt06empmJpaSnSNI1arRbVajUKhULMzs4KpQAAAAAeQa5er9ezLsRBGx4ejjt37rT1fMrlcrHToZbL5W29qiqVSivcqlarMTExsa03VTebm5sxNDQUGxsbMTg42PMxAAAAADyOeslE+m6MqSRJIk3TjgFSrVbrut7DvZ9qtVpcuHDhoIsHAAAAwE/13aN8SZJ0nJ7P5yNN067rFQqFtm0kSRITExNtyywuLsbIyEisr6/H6upqzM/PH0iZAQAAAE6ivgumumkGSnsxPz+/bTyqsbGxiHgQYFUqlSgWi7G0tHSwBQUAAAA4IY59MFWtVuP69eu7Ljc3N9cKjzrZayi1srLScfrWHlUREZcuXYpSqdT1sUEAAAAAdnbsg6mpqame3n73cIDUlKZp13lbLSwsxOjo6Lbp1Wq1rRzNMCpJkh0DMQAAAAA667vBzwuFQuTz+Y5jTT08ZlQntVptWw+oNE2jWCy2bbM5XtVewi4AAAAAtuu7YCqi8Vjf1jfwVavVmJ6ebn1OkiTK5XLHdZMk2RY25fP5mJmZaZteqVRiamrKY3wAAAAA+5Sr1+v1rAtxGMrlcitIunHjRtsb9CqVSszPz8fq6uq29UZHR2NpaWnb43lpmkalUml9Xltb2/GtfJubmzE0NBQbGxsxODj4qIcDAAAA8FjoJRPp22Aqa4IpAAAA4CTqJRPpy0f5AAAAADj+BFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFN09eGHH0Yul4tcLhcffvhh1sWBI6cNgHYA2gBoB6ANHC7BFAAAAACZOJ11AQ5LmqaxuLgYS0tLsby8vKd1yuVy5PP51vozMzM9zQcAAABg7/qyx9TKykosLi5Gmqaxvr6+p3XK5XJERExPT8f09HSMjY1FqVTa83wAAAAAepOr1+v1rAtxWKrValy9ejVu3bq167LDw8Nx586dVo+oiIhcLhfN07Pb/Idtbm7G0NBQbGxsxODg4CMdR1Y+/PDDeO655yIi4oMPPohnn3024xLB0dIGQDsAbQC0A9AGetdLJtKXPaZ6lSRJpGnaFjo11Wq1XecDAAAA0DvBVDSCqU7y+XykabrrfAAAAAB617eDnx+EkZGRWF9f79hTauv8TpqP+G1ubh5W8Q7d1tdgbm5uxr179zIsDRw9bQC0A9AGQDsAbaB3zSxkL6NHHftgqlqtxvXr13ddbm5uLsbGxg5037sNnL7T/Pfffz8iIj73uc8daJmy8pnPfCbrIkCmtAHQDkAbAO0AtIHevP/++zE0NLTjMsc+mJqamoqpqalD3UehUOg4PU3TKBQKu87v5DOf+Uz81V/9VTz//PORy+UOrKwAAAAAx1m9Xo/3339/T0HesQ+mjkKhUIh8Ph9JkmwLmiYmJiIidp3/sFOnTsVnP/vZwykwAAAAwDG2W0+ppr4e/Lzbo3ZJkkS5XG6bNjc31/aGvWq1GtPT03ueDwAAAEBvcvW9jET1mEmSpDU21crKSszMzMT4+HjrkcBKpRLz8/Oxurratl65XG71iLpx40bMz8/3NB8AAACAvevLYIqDUS6XW28kTNM0ZmZmsi0QHKJarRYLCwsxOTkZhUIhlpeX2wLtCG2C/pOmaSwuLsbS0lIsLy9vm79bndcmeNzt1AZcFzgpmk+SNP/TfmFhYdt81wL63U7twPXg8Bljio6aDbP5uGKtVotSqbTtQgX9Ik3TqNVqUa1Wo1AoxOzs7LaLTYQ2Qf9YWVmJmzdvRpqmHR99363OaxM87nZrA64LnASzs7NtT4GUSqWYnJxsBbWuBZwEu7UD14PDp8cUHQ0PD8edO3daqW9ERC6XC9WFflWtVmNiYqKtzm+lTdCvqtVqXL16NW7dutU2fbc6r03QL7q1AdcF+l2aplEsFmNpaalVj1dWVuL8+fOxuroahULBtYC+t5d24Hpw+Pp68HP2J0mSSNO0Y8PbOgA8nBTaBCfNbnVem+Ck0wboFzdv3owkSVqfm+PppmnqWsCJsVM72I12cDA8ysc2WxvlVvl8fk+NEx5Xi4uLMTIyEuvr67G6utrq0qtNcNLsVue1CU4K1wX6WT6fj7t377ZNa/4iXSgU4ubNm13Xcy2gX+zWDppcDw6XYIo9azZE6EdjY2MR8eACVKlUWt16u9EmOGmadb5bV3Ztgn7iusBJdPXq1VhYWOj6cz7CtYD+93A7cD04fIIp9kzDop9t/R+RiIhLly5FqVTa8X86tAlOmt3qvDZBP3Fd4KSZnZ2Ny5cvtwZw7sa1gH7WqR24Hhw+Y0yxzcMNrylN067z4HFXrVbbPjf/hyRJEm2CE2e3Oq9NcBK4LnCSVKvVGB0dbXvFvWsBJ02ndtCcvpXrwcETTLFNoVCIfD7f8XnZiYmJDEoEh6v5No6tdb75PyDNGy9tgpNktzqvTdDvXBc4SZrj6TR7iDTHj3It4CTp1g5cD46GYIqO5ubm2t4iUK1Wd+3WC4+rfD4fMzMzbf+rUalUYmpqqvU/ItoE/apbV/Pd6rw2Qb/o1AZcFzgpVlZWYmVlJcbGxiJJkkiSJCqVSoyMjESEawEnw07twPXgaOTq9Xo960JwPJXL5VYDvHHjRuvNA9CP0jSNSqXS+ry2tratzmsT9JMkSaJarcb169djZWUlZmZmYnx8PKamplrL7FbntQkeZ7u1AdcF+l2apnH27NmO4+Rs/RXRtYB+tpd24Hpw+ARTAAAAAGTCo3wAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmBFMAAAAAZEIwBQAAAEAmTmddAAAADtbs7GycOXMm1tbWIk3TKBaLMTExkXWxAAC2EUwBAPSR8+fPx7Vr12JsbKw1rVgsxoULFyKfz2dXMACADjzKBwDQJyqVShQKhbZQamVlJarVaiRJkmHJAAA6E0wBAPSJ5eXlSNO0bVo+n4+ZmZm2sAoA4LgQTAEA9IlCoRC1Wi1GR0ejVCpFtVqNQqEQ8/PzWRcNAKCjXL1er2ddCAAADsbs7Gzbo3v5fD7u3LljfCkA4FjSYwoAoI/Mz8/H6upq1Ov1WFhYiDRNY3Z2NutiAQB0JJgCAHjM1Wq1GB4ejpWVlbbp09PTMTExobcUAHBsCaYAAB5zCwsLEdEYY+phSZJEqVQ66iIBAOyJYAoA4DE3OTkZS0tL23pGzc7ORqlU6hhYAQAcBwY/BwDoA5VKJVZXV+PMmTMREbG2thaTk5MxMTGRcckAALoTTAEAAACQCY/yAQAAAJAJwRQAAAAAmRBMAQAAAJAJwRQAAAAAmRBMAQAAAJAJwRQAAAAAmRBMAQAAAJAJwRQAAAAAmRBMAQAAAJAJwRQAAAAAmRBMAQAAAJAJwRQAAAAAmRBMAQAAAJAJwRQAAAAAmfj/ASVAjDiRz05wAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1200x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(4.6202e-12, dtype=torch.float64)\n",
      "tensor(2.1767e-14, dtype=torch.float64)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABKYAAAGGCAYAAABBiol3AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAM7NJREFUeJzt3V9sXNmdH/hfSe3x3yYvqTYyY4wDq+hOHj1NiVnsZh6CbhJBNg8LWFWSgUG8mdhidQ+wb+OiCOzjwlKV/TwWS0Ym8Ow+iCw3sA9BsM3qnodMJoAlVbexT0nMy8YY62SQIXXJtj2emW7VPtBVYonFP6Umdani5wM05Lr3nnvPvTxHZX51zrmFTqfTCQAAAAB4xs7lXQEAAAAAzibBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5eCHvCpykLMtieXk5VlZWYnV19Uhl6vV6JEnSK1+tVofaDwAAAMDRjOyIqXa7HcvLy5FlWWxubh6pTL1ej4iI+fn5mJ+fj+np6ahUKkfeDwAAAMDRFTqdTifvSpykZrMZN2/ejAcPHhx67MTERKyvr/dGREVEFAqF6D6iw/YDAAAAcHQjO2JqWGmaRpZlfaFTV6vVOnQ/AAAAAMMZ6TWmhpGm6cDtSZJElmWH7h/k0aNH8bOf/SxefPHFKBQKx1VVAAAAgFOr0+nEBx98EF/4whfi3LmDx0QJpg4xOTkZm5ubA0dK7d4/yM9+9rP44he/eIK1AwAAADidfvrTn8Zv//ZvH3iMYOoQhy2cftD+F198MSJ2fhBjY2PHWi8AAACA02h7ezu++MUv9nKRgwimfq1YLA7cnmVZFIvFQ/cP0p2+NzY2JpgCAAAAzpSjLGtk8fNfKxaLkSTJwLWkZmdnD90PAAAAwHBGPpjab6pdmqZRr9f7ti0uLva9Ya/ZbMb8/PyR9wMAAABwdIVOp9PJuxInIU3TaDabcffu3Wi321GtVmNmZiZKpVJERDQajajVarG2ttZXrl6v96bm3bt3L2q12lD7d9ve3o7x8fHY2toylQ8AAAA4E4bJQ0Y2mDoNBFMAAADAWTNMHjLyU/kAAAAAOJ0EUwAAAADkQjAFAAAAQC4EUwAAAADkQjAFAAAAQC4EUwAAAADkQjAFAAAAQC4EUwAAAADkQjAFAAAAQC4EUwAAAADkQjAF8Iw0Go0ol8t5VwMAAODUEEwBPCMPHjyIZrOZdzUAAABODcEU7NJsNiPLsryrwYhaWlqKTqfzTK+pTQMAAKeZYAp2KZfLcf/+/byrAcdGmwYAAE4zwRScdf/m30T84Af9237wg53tAAAAcIIEU6dcpxPxq1+d7v+OY2ZSs9mMS5cuRaFQiEuXLkWr1YqIiIWFhSgUCpGmad/xc3NzMTU11fucpmmUy+WYmJjonaPdbu+5Trvdjrm5uZiYmIiJiYmYm5uLdrsd5XI5CoVC79yFQqH3eeSdOxfxx3/8OJz6wQ92Pp/L/6+Ho/xcG41G7+fY/dlOTU312tAwx9Xr9Zibm4uIiEqlEhMTE71pcE+Wq9frvXJHbaf1ej0mJib6Pnfv6cn2n2VZ796npqai0WgM9WyO0qb363cAAADPygt5V4CD/c3fRPyzf5Z3LQ727/5dxKc+9fTlG41GVCqVWFpaisuXL0er1Yq5ublYW1uLWq3WCwTW1tZ6x7dard7niJ1fsCcnJ2NlZSWKxWIsLCzEa6+9Fuvr65EkSURE75f/UqkUKysrERGxuroad+/ejTt37sTi4mJcunQpVlZWYnZ29ulv6Hnz9a/v/PnHfxzxJ38S8eGHEb//+4+35+goP9csy6LVasX169djcXExKpVKLCws9NpMsVgc6rg0TePSpUuRZVksLi5GkiS9NlmtVqNWq0WaprGwsBD37t2LlZWVI7fTJ21sbES73Y7r169HrVaLiJ1ArFwuR7FY7NWzVqtFpVKJ2dnZXj0PezaHtemD+l33GgAAACdNMEWusizr/XI8Pz8fERHT09OxsbERtVotlpaWYnV1NaampqJcLvd+Ue/+Mt5VrVb7znvnzp2YmJiI5eXl3nnL5XJfKBURfb+sd8+XJEkv9Dgzvv71x6HUCy+cilAq4mg/164HDx70/vfs7GxMTEz02tAwx6VpGrOzs33HVSqVXigVsdNGp6enY2pqKprNZpRKpSO10/3UarVeW1xYWOiFUKVSKSIiJicneyOiuuc77NkkSbJvmz5KvwMAAHgWBFOn3Cc/uTMi6TT75Cefvmx3UeZKpRKVSqVv3/T0dO9/d3/pbzabUa1We7+w76f7S3h3tEqappGmqV+49/ODHzwOpT78cOfzKQmndnvy53rQcaVS6dCpafsdt7udtNvtSNN0T/ssFosxPT0dd+/e7bXHYdtp1+XLl/vOGxG9KYURj/vC5ubmgfcScfiziTh6vwMAADhp+S8iw4EKhZ1pcqf5v+NYiunhw4fR6XT6/ts9YqVYLPZGdiwuLg48R6vVinK5HFNTU33r+EREb+0fU5QG6K4p9fu/H7G6uvPn7jWncnbQz/UgxWLxwCDnoON2t5NuiDM5OTmw7O51pY7STgcZNEJv0PWe9LTPpuuwfgcAAHDSBFPkave6PgdptVqxvLwcpVIpXnvttT37y+VylMvlmJubi9XV1Xj48OFTXedMevSof02pr3995/OjR/nWKw7/uR4kTdMjBZFPHvdkSHRQ23my7GHtdJD9po0eNp304zwb/QEAADgtBFPkqjsd6ubNm3v2dd+G1n372J07d3rrQ+2efpRlWTSbzbhz507Mz88PDCOKxWIUi8WBU/m619m9mPaZ8i//5d5pe1//+s72HB3l57qfNE2j2Wweuoj9UY6bnZ2NJEn2tJ12ux3tdjuuXbvWO9dB7fQ4HfXZ7Nemj9LvAAAAngVrTJG7O3fuxKVLl3rr3WxubsbKykrvz3K53LcQ9MrKSkxNTcWlS5d6izwnSRI3b96MJElicnJy4C/cS0tLMTc313vrWZZlsbq6Gvfv3+9NX0qSJO7evdsLsWq12tlbCP2UOOrPtWtubi4WFhYiy7K4fv16JEkycDrdUY/b7c6dO1EulyNiZ6RS9618u9vlYe30OA3zbPZr04f1OwAAgGfBiClyNz09HWtra5Gmabz22mu9AODOnTtRqVQiTdO4c+dO7/juL9iVSiXa7Xbv2O6IlevXr8fc3FzMzs7G1NRUr1z3TWu7j3vyl/DFxcVoNpu9OpCvo/xcuyqVSiwsLES5XI7Lly/HgwcPBoaKRz1ut1Kp1Gs7c3NzUavVYnFxMVZXV3vnPEo7PU5HfTb7temD+h0AAMCzUuh0Op28KzGqtre3Y3x8PLa2tmJsbCzv6sBIqtfrsbCwEIf9VXbU4wAAAPh4hslDjJgCAAAAIBeCKQAAAAByMfKLn9fr9b43U1Wr1QOPL5fLce3atSgWiwNfG99qtXqLaBeLxVhdXY2ZmZnegscAAAAAHM1IrzFVr9cjInphVKvVipWVlT2vfd9tamoq0jTds71UKsXKyko0m824fv16ZFkWxWIxFhYW9n3jljWmAAAAgLNmmDxkpIOpiYmJWF9f7xv5VCgUDlz8uF6v7xlV1Wg0euFTs9mM2dnZQ9/iFSGYAgAAAM4ei59HRJqmkWXZwACp1WrtW+7JKXmtVisuX7583NUDAAAAOPNGdo2pQdPxIiKSJIksy/YtVywW+86RpmnMzs72HbO8vByTk5OxubkZa2trUavVjqXOAAAAAGfJyAZT++kGSkdRq9X2rEc1PT0dEY8DrEajEeVyOVZWVo63ogAAAAAjbmSn8u3nqKFUu90euL1YLPaNqrp69Wo0m80DR2EBAAAAsNfIBlO7w6Pdum/TO8zS0lJMTU3t2d5sNvs+d9ew2m/qIAAAAACDjXQwlSTJwMDoyTWjBmm1WnsWTs+yLMrlct85uyOljhJ2AQAAAPDYyAZTERGLi4t9b+BrNpsxPz/f+5ymadTr9YFl0zTdEzYlSRLVarVve6PRiFKpNPDtfwAAAADsr9DpdDp5V+Ik1ev1XpB07969vjfoNRqNqNVqsba2tqfc1NRUrKys9BY778qyLBqNRu/zxsbGvm/l297ejvHx8dja2oqxsbHjuB0AAACAU22YPGTkg6k8CaYAAACAs2aYPGSkp/IBAAAAcHoJpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIxQt5V+Ck1ev1SJIkIiKyLItqtXrg8a1WK5aWlmJubi6KxWKsrq7GzMxMlEqlpz4nAAAAAHuN9Iiper0eERHz8/MxPz8f09PTUalUDiyTZVm0Wq2oVCpRqVRiampqTyg17DkBAAAA2KvQ6XQ6eVfipExMTMT6+npvdFNERKFQiINuudlsxuzsbF+Zpz3n9vZ2jI+Px9bWVoyNjT3tbQAAAAA8N4bJQ0Z2xFSappFl2cCAqdVqnZpzAgAAAJxVI7vGVJqmA7cnSRJZlh1Ydnl5OSYnJ2NzczPW1taiVqt97HMCAAAA0G9kg6n9dAOn/UxPT0dERLFYjIiIRqMR5XI5VlZWnvqcAAAAAOx15oKpwwKkbiDVdfXq1ahUKgeOiBJKAQAAAAxvZNeYejJg6sqybN99ETuLn+/WXU8qTdOnPicAAAAAe410MJUkycB1oWZnZweWybIsyuVyX5nuSKlisfhU5wQAAABgsJENpiIiFhcX+96W12w2Y35+vvc5TdOo1+u9z0mSRLVa7Rv91Gg0olQq9UZOHXZOAAAAAI6m0Ol0OnlX4iTV6/Ve0HTv3r3eG/YidkKnWq0Wa2trvW1ZlkWj0eh93tjY6Ctz2Dl3297ejvHx8dja2oqxsbFjuycAAACA02qYPGTkg6k8CaYAAACAs2aYPGSkp/IBAAAAcHoJpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFwIpgAAAADIhWAKAAAAgFy8kHcFTlq9Xo8kSSIiIsuyqFarRyoTEbG2thYREUtLS719rVYrlpaWYm5uLorFYqyursbMzEyUSqXjrzwAAADACBvpYKobMM3Pz0fETqhUqVT6gqYnLSwsRK1W632uVCoxNzcXq6urEbETbrVarWg2m1EsFmNhYUEoBQAAAPAUCp1Op5N3JU7KxMRErK+v90ZMRUQUCoXY75azLItyuRwrKyu9Mu12Oy5duhRra2tRLBaj2WzG7Oxs3zn3s729HePj47G1tRVjY2PHcEcAAAAAp9swecjIrjGVpmlkWTYwQGq1WvuWu3//fqRp2vtcLBYjYie0AgAAAOD4jOxUvt3h0m5JkuwbMiVJEg8fPuzb1g2xugFVRMTy8nJMTk7G5uZmrK2t9U39AwAAAOBoRjaY2k83UDqqmzdvxtLSUm/k1fT0dEQ8DqoajUZv+h8AAAAARzeyU/n2M0wotbCwENeuXestnh6xE0jtHj119erVaDabpvoBAAAADOnEg6n33nsvtre3T/oye+wOj3bLsmzffbs1m82YmpqKarW6Z/tu3ZFU+00dBAAAAGCwEw+mvv3tb0eSJPHyyy/HG2+8EW+++Wa8//77J33ZKBaLkSTJwMBodnb2wLLddaW6I6WyLOstpl4ul/vO2R0pdZSwCwAAAIDHTjyYWl5ejkePHsX3vve9GB8fj9u3b0exWIwLFy7EtWvX4p133jmxay8uLva9ga/ZbPZNy0vTNOr1el+Zdrsd7XY7pqenI03TSNM0Go1GTE5ORpIkUa1W+0KoRqMRpVJp4Nv/AAAAANhfodPpdPK4cLPZjOXl5ciyLNbX16PZbMZXvvKVY79OvV7vBUn37t3re4Neo9GIWq0Wa2trEbEz+unixYsD14vqPqYsy6LRaPS2b2xs7PtWvu3t7RgfH4+tra0YGxs7rlsCAAAAOLWGyUNyC6YiIu7cuRPXr1+Pdrsd8/Pz0Ww240tf+lJe1Tl2gikAAADgrBkmDznxqXw3btyIc+fOxczMTHz3u9/tW19qa2srIiKmp6fj/v37exYWBwAAAGB0nXgwNTMzE48ePYqFhYX40Y9+FMViMc6fPx8XLlzoHfPee+9FRMTFixdPujoAAAAAnBIvnPQFkiSJd955J0qlUpRKpYiIWF9fj8nJyRgfH4+IiFKpFOVyOWZmZk66OgAAAACcEs9kjamtra148OBBvPrqqyd9qVPFGlMAAADAWXOia0xdvXo1zp8/H+fPn48/+IM/OFKZ8fHxMxdKAQAAAHCwoYKpxcXFmJmZifv378fdu3fjRz/6Ud/0u8XFxXj55Zd7a0h97Wtfix//+MfHXmkAAAAAnn9DTeW7ceNG3Lp1q29bvV6PBw8eRJqmMTExEXNzc7GxsRFpmkar1Yqtra1oNBrxjW9849grf9qZygcAAACcNcPkIUMtfv7SSy/t2VatVuPy5cvx+uuvxze/+c09+5vNZszPz0dEnMlwCgAAAIDBhprKt7Gxse++ycnJgdtLpVKkaRrf+973Ynt7e7jaAQAAADCyhgqmisVivP/++wO3p2m6b7kkSWJlZSVu3rw5dAUBAAAAGE1DBVPXr1+P27dv7xn5lGVZJElyYNmLFy/uO6oKAAAAgLNnqDWmIiJu3boVr7/+ely9ejVeffXViIh46623jlS2UCgMezkAAAAARtRQI6a6bt++HQ8fPow33ngj3nzzzSOXO2iNKgAAAADOlqFHTHVduXIlrly5Eu+++27cuHGjNxpqbm6uN5Jqt3fffddUPgAAAAB6Cp1Op3OcJ3z77bdjdXW1F1RduHAh1tbWIkmSM7f4+fb2doyPj8fW1laMjY3lXR0AAACAEzdMHnLswdST3n777Wi327GxsREvvfRSTE9Px+XLl89EUCOYAgAAAM6aUxVMPendd9+NVqsVGxsbUSgUYmZmJr761a8+yyo8M4IpAAAA4KwZJg956jWmntYrr7wSr7zySu/z+vr6s64CAAAAAKfAU72V7zhdvHgx7yoAAAAAkIMjBVNvv/12XLt2Lf70T//0pOsDAAAAwBlxpGDqtddei6tXr8bNmzfjwoULsbi4GO+9994JVw0AAACAUXbkqXxXrlyJt956K9I0jWKxGN/85jfj5Zdfju9+97uxvb19knUEAAAAYAQNvcbU+Ph4XL9+Pe7fvx9vvfVW/NVf/VVMT0/HP/2n/zS+//3vn0QdAQAAABhBH2vx84sXL8atW7fiJz/5Sdy6dSvu378fk5OTce3atXjnnXeOq44AAAAAjKBjeyvfK6+8Erdv347Nzc24evVq3L59Oy5cuBBvvPGG9agAAAAA2OPYgqndrly5EsvLy5GmaUxPT/etR/X++++fxCUBAAAAeM6cSDDV9eR6VJ1OJ2ZnZ2NmZia+//3vWzQdAAAA4Aw70WBqt4sXL8a3vvWt+MlPfhKNRiPu378fX/rSl6xHBQAAAHBGFTqdTifPCrz99tvRarXi5s2bJ3L+er0eSZJERESWZVGtVj92maOec3t7O8bHx2NrayvGxsae+h4AAAAAnhfD5CG5B1MnqV6vR0T0gqNWqxUrKyuxtLT01GWGOadgCgAAADhrTlUw9d5770WxWMwlmJmYmIj19fXe6KaIiEKhEAfd8mFlhjnncx9MffObEefORTQaj7fNz0c8ehTx/e8/PuYnP4n4r/81olCI6HQi/uIvIj78MOITn4gYG4v467/e+e/v/m7nmH/wD3aO/8UvIj76KOL8+YjPfz7iL/9yp/ynPhWRJBF/9Vc753nhhYiXXor45S8jPvhg55hPfzriH/2jnTr8t/+289/f/E3Eb/zG47r+1m/t/PkXf7Fz/Yidc3X/7NbnM5/Z2f/o0c79nj//+LoffbTz5/h4xNbW4/P85m/u3Nv2dsR//+872/7xP358D7/1WxFf/nL/czrsWXaP+7f/dud//8N/uPOcfvKTnXp+/vM710zTnfoVChEvvrhT/24dzv16du758xG/+tXjY375y516dY/59Kcf3/OjRzvP/Pd+b2f/f/kvEX/+5zvXePHFx8/y7/29iHv3drZ3n1333IXC4+e3+9yf+ETE3/7tTn0+/emdcuPjERsbO8+209kpMzW18zPsPt8LF3bu/z/9p53P//yf731O2qa2qW1qm9qmtqltapvaprapbWqb3Trsbnvb2zvHfPKTj48f1H6ex7b5nBgqD+mcsHK53CkUCp0vf/nLnddff73zwx/+sLO+vn7Sl+2sra11Bt1eRHRWV1efqsyw59za2upERGdra+sp7uAUuH6900mSnT8Hfe5u+8xnOp1PfKLTKRQ6nZ3u9ez+233N8+efbR3On99bl0Jh51l85jN7n9Nhz7K7vXve7vkOuuZxP8vx8Wf/MzzsZ3v+/ODnpG1qm9qmtqltapvaprapbWqb2qa2edba5nNimDzkmU3la7Va0Wq1ot1uR6vViomJiZidnY1KpRKvvvrqiVxvbm4unry9iYmJuHPnTpRKpaHLJEky1Dmf+xFTEfE3/+t8fHR3JR51zkUhHsW9L5XjT3630XfMv/iz+fgf1/6v+MRHv4zCAefqRCEKcTLN7W/Pfzr+5H+6HV/+yz+L3/3P/zrOx0dDn6Nbs4Pu4aPC+Sh0HsW5fe+jEH97/tPxH6d+b+Bzmnl/JTqx/7PsHnfYPTyKwgF1OLr97jn79G/GDy/X4l/8+evxGx/99cByBz2ng849zDGPCufj37/8r/Z9Toc9T21zN23zsHMPc4y2uUPbfEzbPD7a5tPRNrv10Da1zeGO0TZ3aJuPnZa22Tl3Ps5941/1jzR7jpyqqXz7aTabsby8HFmWxfr6ejSbzfjKV75ybOffL2SampqKhYWFmJ+fH7pMsVgc6pyjEEytrkb8D//zhSjEo+jEufhffndj4HH/959diM98tB3nOx8O3N+JQvz8hSQ+92F24F/IHxVe6DvHk3+Bf1R4Ic51PtpzjgcTs/Gtr6z26vLih5uHnvvJ+j0qnI+IOOQeJuI/vzgdlx629q3/L8+PHficDnuWB91Dtx7tiddi+uHbhz7LQc9q93m69/zkca/+k53//Z0fzw2815+/kMRnP9x6qnPvPubnL0wc2G4+eGHy0OekbWqbu2mb2uZ+9dA2d2ibu+9B24zQNrvltM299dA2d2ibu+9B24x4dm3z7z43Gb/xwf7P6bQbJg954RnVaY9SqRQPHz6M69evR7vdjm984xvRbDbjS1/60oled3NzcAP/OGWe5pzPi9/9P+fjhU89ikedc3G+8CiWk/n4i/+9P7H9+//HfHz2/K+i8OHgDhURUYjOoX8RR+z9i/DJ4/frtK/88j/Ev/4nP4jPvfdn8blHW0c695PXOWh/95jPdbZiOnt732POdz6Kz77wq32f02c+9eu5x48GP8vucfvdQ7ce09nBfxHv1OXw+9nvmP/n//2t+P/+t1r8/f/4Hwbu/9yH2VOfe/cxn/tw88B/RfhcZ+vA53TY89Q2d19D2zzs3LuP0TYPLtu9jrb5mLa5l7apbfZdR9vUNp+gbe6lbWqbu33ir7d21uJ6TkdMDeXEJhT+2sLCQqdQKHQuX77c+c53vtO3vtR3vvOdvmOf/PxxWGPqGJhXffB/5lU/m5+tOf/aprapbWqb2qa2qW1qm9qmtqltapvPlWHykBMfMTUzMxOPHj3qTd2rVqtRKBQiSZJYXFyMiJ039/3O7/xOXLx48diuWywWI0mSSNM0isVi377Z2dmnLjPsOZ9rjx5FlMuPE9run48e9R8zM+NNFBF730Tx5HM67Fl2P3/+8zv/21tSHr+JYtBz0ja1TW1T29Q2tU1tU9vUNrVNbVPbHPW38j35nEbQia8x9fbbb0ehUOhb4Hx9fT0mJydjfHw8IiK+/OUvR7lcjpmZmfjqV796bNeu1+uRJElv7admsxmrq6uxtLQUERFpmkaz2YxqtXrkMoft320U1pgCAAAAGMapW/x8a2srHjx4cCJv3ztMvV7vjW66d+9e1Gq13r5GoxG1Wi3W1taOXOYo+7sEUwAAAMBZc6LB1NWrV+OHP/xhRERUKpX4oz/6o6ev6YgTTAEAAABnzTB5yLlhTry4uBgzMzNx//79uHv3bvzoRz+KmZmZvv0vv/xynD9/Pi5cuBBf+9rX4sc//vHT3QUAAAAAI22oEVM3btyIW7du9W2r1+vx4MGDSNM0JiYmYm5uLjY2NiJN02i1WrG1tRWNRiO+8Y1vHHvlTzsjpgAAAICzZpg8ZKi38r300kt7tlWr1bh8+XK8/vrr8c1vfnPP/maz2Vso/CyGUwAAAAAMNtRUvo2NjX33TU5ODtxeKpUiTdP43ve+F9vb28PVDgAAAICRNVQwVSwW4/333x+4PU3TfcslSRIrKytx8+bNoSsIAAAAwGgaKpi6fv163L59e8/IpyzLIkmSA8tevHhx31FVAAAAAJw9Q60xFRFx69ateP311+Pq1avx6quvRkTEW2+9daSyhUJh2MsBAAAAMKKGGjHVdfv27Xj48GG88cYb8eabbx653EFrVAEAAABwtgw9YqrrypUrceXKlXj33Xfjxo0bvdFQc3NzvZFUu7377rum8gEAAADQU+h0Op3jPOHbb78dq6urvaDqwoULsba2FkmSnLnFz7e3t2N8fDy2trZibGws7+oAAAAAnLhh8pBjD6ae9Pbbb0e73Y6NjY146aWXYnp6Oi5fvnwmghrBFAAAAHDWnKpg6knvvvtutFqt2NjYiEKhEDMzM/HVr371WVbhmRFMAQAAAGfNMHnIU68x9bReeeWVeOWVV3qf19fXn3UVAAAAADgFnuqtfE9re3t7z7aLFy8+yyoAAAAAcEo8k2Bqe3s7Xn/9dW/lAwAAAKDnRKfyvf/++1Gr1aLRaESn0+m9qQ8AAAAATmTE1Pvvvx/Xrl2LqampWFpaitdeey3m5+dP4lIAAAAAPKeONZh67733eoHUyspKXLlyJdbW1uKtt96KS5cuHeelAAAAAHjOHctUvvfeey8WFhai1WpFp9OJUqkUtVrNwuYAAAAA7OtjBVPvvPNOLCwsRLvdjk6nE/Pz87GwsCCQAgAAAOBQTxVMDQqkarVajI+PH3f9AAAAABhRQwVTb775ZiwsLESaptHpdKJarcbi4qJACgAAAIChHTmY+uEPfxjlcjkmJiaiWq3GjRs3BFIAAAAAPLUjv5XvypUrcf/+/SiXyzE1NSWUAgAAAOBjOXIwFRExPT0dt2/fjtdeey1u3LgR3//+90+qXgAAAACMuKda/PzixYtx69atWF9fjxs3bsRLL70Uf/iHf3jcdQMAAABghD1VMNXVDai2trZ6AdX8/HyMjY0dV/0AAAAAGFEfK5jqGh8f7wVU3/72twVUAAAAABzqWIKprt0B1dLSUhQKhVwDqnq9HkmSRERElmVRrVaPVCYiYm1tLSIilpaWevtarVYsLS3F3NxcFIvFWF1djZmZmSiVSsdfeQAAAIARd6zBVNf4+Hh861vfioiI73znO7G5uRmdTuckLrWvbsA0Pz8fETuhUqVS6QuanrSwsBC1Wq33uVKpxNzcXKyurkbETrjVarWi2WxGsViMhYUFoRQAAADAUyp0nlFiVK/X49atW7G5ufksLhcTExOxvr7eGzEVEVEoFPYNyLIsi3K5HCsrK70y7XY7Ll26FGtra1EsFqPZbMbs7GzfOQ+yvb0d4+PjsbW1ZVojAAAAcCYMk4ece0Z1imq1+sxCqTRNI8uygQFSq9Xat9z9+/cjTdPe52KxGBE7oRUAAAAAx+tEpvLlbXe4tFuSJPuGTEmSxMOHD/u2dUOsbkAVEbG8vByTk5OxubkZa2trfVP/AAAAADi6kQym9tMNlI7q5s2bsbS01Bt5NT09HRGPg6pGo9Gb/gcAAADAcJ6LYKrZbMbdu3cPPW5xcbEXHg0yTCi1sLAQ165d6y2eHtE/cioi4urVq1GpVPadNggAAADA/p6LYKpUKg319rsnA6SuLMv23bdbs9mMqampvlCqu313PbphVJqmBwZiAAAAAOz1zBY/f5aKxWIkSTJwranZ2dkDy3bXleqGUlmW9RZTL5fLfefsrld1lLALAAAAgH4jGUxF7Ezr2/0Gvmaz2TcCKk3TqNfrfWXa7Xa02+2Ynp6ONE0jTdNoNBoxOTkZSZJEtVrtC6EajUaUSiXT+AAAAACeQqHT6XTyrsRJqdfrvSDp3r17fW/QazQaUavVYm1tLSJ2Rj9dvHhx4Fv7uo8oy7JoNBq97RsbGwe+lW97ezvGx8dja2srxsbGjuOWAAAAAE61YfKQkQ6m8iaYAgAAAM6aYfKQkZ3KBwAAAMDpJpgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAAByIZgCAAAAIBeCKQAAAABy8ULeFThJ9Xo9kiSJiIgsy6JarR54fKvViqWlpZibm4tisRirq6sxMzMTpVLpqc8JAAAAwGAjO2KqXq9HRMT8/HzMz8/H9PR0VCqVA8tkWRatVisqlUpUKpWYmpraE0oNe04AAAAABit0Op1O3pU4CRMTE7G+vt4b3RQRUSgU4qDbbTabMTs721fm45xze3s7xsfHY2trK8bGxp7mNgAAAACeK8PkISM5YipN08iybGDA1Gq1Ts05AQAAAM6ykVxjKk3TgduTJIksyw4su7y8HJOTk7G5uRlra2tRq9U+9jkBAAAA2Gskg6n9dAOn/UxPT0dERLFYjIiIRqMR5XI5VlZWnvqcAAAAAAz2XARTzWYz7t69e+hxi4uLvXBpkMMCpG4g1XX16tWoVCoHjogSSgEAAAA8necimCqVSn1vxzvMkwFTV5Zl++6L2AnAdl+nu55UmqZPfU4AAAAABhvJxc+LxWIkSTJwXajZ2dmBZbIsi3K53FemO1KqWCw+1TkBAAAA2N9IBlMRO9P6dr8tr9lsxvz8fO9zmqZRr9d7n5MkiWq12jf6qdFoRKlU6o2cOuycAAAAABxdodPpdPKuxEmp1+u9oOnevXu9N+xF7IROtVot1tbWetuyLItGo9H7vLGx0VfmsHM+aXt7O8bHx2NrayvGxsaO5Z4AAAAATrNh8pCRDqbyJpgCAAAAzpph8pCRncoHAAAAwOkmmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgF4IpAAAAAHIhmAIAAAAgFy/kXYGTVK/XI0mSiIjIsiyq1eqBx5fL5bh27VoUi8Veua5isRitViuWlpZibm4uisVirK6uxszMTJRKpRO6AwAAAIDRNbLBVL1ej4iI+fn5iIhotVpRqVRiaWlp3zLtdjuazeae7aVSKVZWViLLsmi1WtFsNqNYLMbCwoJQCgAAAOApFTqdTifvSpyEiYmJWF9f7xv5VCgU4qDbrdfre0ZVNRqNXrjVbDZjdnZ2z2iq/Wxvb8f4+HhsbW3F2NjY0PcAAAAA8LwZJg8ZyTWm0jSNLMsGBkitVmvfck+Ofmq1WnH58uXjrh4AAAAAMaJT+dI0Hbg9SZLIsmzfcsVise8caZrG7Oxs3zHLy8sxOTkZm5ubsba2FrVa7VjqDAAAAHDWjGQwtZ9uoHQUtVptz3pU09PTEfE4wGo0GlEul2NlZeV4KwoAAABwBjwXwVSz2Yy7d+8eetzi4mIvPBrkqKFUu90euH33iKqIiKtXr0alUtl32iAAAAAA+3sugqlSqTTU2++eDJC6sizbd99uS0tLMTU1tWd7s9nsq0c3jErT9MBADAAAAIC9RnLx82KxGEmSDFxr6sk1owZptVp7RkBlWRblcrnvnN31qo4SdgEAAADQbySDqYidaX2738DXbDZjfn6+9zlN06jX6wPLpmm6J2xKkiSq1Wrf9kajEaVSyTQ+AAAAgKdQ6HQ6nbwrcVLq9XovSLp3717fG/QajUbUarVYW1vbU25qaipWVlb2TM/LsiwajUbv88bGxoFv5dve3o7x8fHY2tqKsbGxj3s7AAAAAKfeMHnISAdTeRNMAQAAAGfNMHnIyE7lAwAAAOB0E0wBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EExxoF/84hdRKBSiUCjEL37xi7yrA8+cPsBZpw+AfgD6AGedPnCyBFMAAAAA5OKFvCtwkrIsi+Xl5VhZWYnV1dUjlanX65EkSa98tVodaj8AAAAARzOyI6ba7XYsLy9HlmWxubl5pDL1ej0iIubn52N+fj6mp6ejUqkceT8AAAAAR1fodDqdvCtxkprNZty8eTMePHhw6LETExOxvr7eGxEVEVEoFKL7iA7b/6Tt7e0YHx+Pra2tGBsb+1j3kZdf/OIX8bnPfS4iIn7+85/HZz/72ZxrBM+WPsBZpw+AfgD6AGedPjC8YfKQkR0xNaw0TSPLsr7QqavVah26HwAAAIDhCKZ+LU3TgduTJIksyw7dDwAAAMBwRnrx8+MwOTkZm5ubA0dK7d4/SHeK3/b29klV78TtfhXm9vZ2fPTRRznWBp49fYCzTh8A/QD0Ac46fWB43RzkKKtHPRfBVLPZjLt37x563OLiYkxPTx/rtQ9bOP2g/R988EFERHzxi1881jrl5Qtf+ELeVYBc6QOcdfoA6AegD3DW6QPD+eCDD2J8fPzAY56LYKpUKkWpVDrRaxSLxYHbsyyLYrF46P5BvvCFL8RPf/rTePHFF6NQKBxbXQEAAABOq06nEx988MGRgrznIph6ForFYiRJEmma7gmaZmdnIyIO3f+kc+fOxW//9m+fTIUBAAAATqnDRkp1jfzi5/tNtUvTNOr1et+2xcXFvjfsNZvNmJ+fP/J+AAAAAI6u0DnKSlTPoTRNe2tTtdvtqFarMTMz05sS2Gg0olarxdraWl+5er3eGxF17969qNVqQ+0HAAAA4GhGNpjieNTr9d4bCbMsi2q1mm+F4AS1Wq1YWlqKubm5KBaLsbq62hdoR+gTjJYsy2J5eTlWVlZidXV1z/7D2rv+wPPuoD7gO4GzpDuTpPuP9ktLS3v2+z5glB3UB3wfnDxrTLGvbufsTldstVpRqVT2fFHBqMiyLFqtVjSbzSgWi7GwsLDnCydCn2A0tNvtuH//fmRZNnDa+2HtXX/geXdYH/CdwFmxsLDQNwukUqnE3NxcL6z1fcCoO6wP+D44eUZMsa+JiYlYX1/vJb8REYVCITQZRlWz2YzZ2dm+Nr+bPsEoajabcfPmzXjw4EHf9sPau/7AqNivD/hO4CzIsizK5XKsrKz02nK73Y5Lly7F2tpaFItF3weMtKP0Ad8HJ2/kFz/n6aRpGlmWDex8uxeAh7NCn+AsOay96w+cdfoAo+T+/fuRpmnvc3c93SzLfB9wJhzUBw6jDxwPU/kYaHfH3C1JkiN1UHheLS8vx+TkZGxubsba2lpvWK8+wVlyWHvXHzgrfCcw6pIkiYcPH/Zt6/4yXSwW4/79+/uW833AKDisD3T5PjhZgimG0u2MMIqmp6cj4vGXUKPR6A3t3Y8+wVnSbe/7DWXXHxglvhM4q27evBlLS0v7/l0f4fuA0fZkH/B9cPIEUwxF52KU7f5XkYiIq1evRqVSOfBfO/QJzpLD2rv+wCjxncBZtLCwENeuXest4rwf3weMqkF9wPfBybPGFAM92fm6sizbdx8875rNZt/n7r+SpGmqT3CmHNbe9QfOAt8JnDXNZjOmpqb6XnPv+4CzZFAf6G7fzffB8RNMMVCxWIwkSQbOmZ2dnc2hRnCyum/k2N3mu/8K0v0/XvoEZ8Vh7V1/YNT5TuCs6a6p0x0l0l0/yvcBZ8V+fcD3wbMhmGJfi4uLfW8SaDabhw7rhedVkiRRrVb7/mWj0WhEqVTq/auIPsEo2m+o+WHtXX9gVAzqA74TOEva7Xa02+2Ynp6ONE0jTdNoNBoxOTkZEb4PGH0H9QHfB89GodPpdPKuBKdXvV7vdcJ79+713j4AoyjLsmg0Gr3PGxsbe9q8PsGoSNM0ms1m3L17N9rtdlSr1ZiZmYlSqdQ75rD2rj/wPDusD/hO4CzIsiwuXrw4cK2c3b8m+j5gVB2lD/g+OHmCKQAAAAByYSofAAAAALkQTAEAAACQC8EUAAAAALkQTAEAAACQC8EUAAAAALkQTAEAAACQC8EUAAAAALkQTAEAAACQC8EUAAAAALkQTAEAAACQC8EUAAAAALkQTAEAAACQC8EUAAAAALkQTAEAAACQixfyrgAAAMdrYWEhLly4EBsbG5FlWZTL5Zidnc27WgAAewimAABGyKVLl+LOnTsxPT3d21Yul+Py5cuRJEl+FQMAGMBUPgCAEdFoNKJYLPaFUu12O5rNZqRpmmPNAAAGE0wBAIyI1dXVyLKsb1uSJFGtVvvCKgCA00IwBQAwIorFYrRarZiamopKpRLNZjOKxWLUarW8qwYAMFCh0+l08q4EAADHY2FhoW/qXpIksb6+bn0pAOBUMmIKAGCE1Gq1WFtbi06nE0tLS5FlWSwsLORdLQCAgQRTAADPuVarFRMTE9Fut/u2z8/Px+zsrNFSAMCpJZgCAHjOLS0tRcTOGlNPStM0KpXKs64SAMCRCKYAAJ5zc3NzsbKysmdk1MLCQlQqlYGBFQDAaWDxcwCAEdBoNGJtbS0uXLgQEREbGxsxNzcXs7OzOdcMAGB/gikAAAAAcmEqHwAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkAvBFAAAAAC5EEwBAAAAkIv/HxB6/RxI9mf9AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1200x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(4.3448e-11, dtype=torch.float64)\n",
      "tensor(1.5944e-14, dtype=torch.float64)\n"
     ]
    }
   ],
   "source": [
    "# Beta beating (approximate)\n",
    "\n",
    "def bbx(Ax, Bx, Cx, nux, ax, bx, axs, bxs, mux):\n",
    "    t = 1.0\n",
    "    cx = (1 + ax**2)/bx\n",
    "    bxf = bxs + (bxs*t**2*(1.0/torch.sin(2*nux*torch.pi)**2)*(4*ax*Bx*Cx*(torch.cos(2*mux) + torch.cos(2*mux - 4*nux*torch.pi)) - 2*cx*Cx**2*(torch.cos(2*mux) + torch.cos(2*mux - 4*nux*torch.pi)) + Ax**2*bx**3*(1 + torch.cos(2*mux) + torch.cos(2*mux - 4*nux*torch.pi)) - 4*ax*Ax*bx**2*Bx*(1 + torch.cos(2*mux) + torch.cos(2*mux - 4*nux*torch.pi)) - 4*ax*bx*Bx*cx*Cx*(1 + torch.cos(2*mux) + torch.cos(2*mux - 4*nux*torch.pi)) + bx*cx**2*Cx**2*(1 + torch.cos(2*mux) + torch.cos(2*mux - 4*nux*torch.pi)) + 2*bx**2*cx*(2*Bx**2 + Ax*Cx)*(1 + torch.cos(2*mux) + torch.cos(2*mux - 4*nux*torch.pi)) - 2*bx*(2*Ax*Cx + (2*Bx**2 + Ax*Cx)*torch.cos(2*mux) + (2*Bx**2 + Ax*Cx)*torch.cos(2*mux - 4*nux*torch.pi)) - 4*Bx*Cx*torch.sin(2*mux) - 4*Bx*Cx*torch.sin(2*mux - 4*nux*torch.pi) + 2*Ax*bx**2*Bx*(torch.sin(2*mux) + torch.sin(2*mux - 4*nux*torch.pi)) + 6*bx*Bx*cx*Cx*(torch.sin(2*mux) + torch.sin(2*mux - 4*nux*torch.pi)) - 2*ax*cx*Cx**2*(torch.sin(2*mux) + torch.sin(2*mux - 4*nux*torch.pi)) - 2*ax*bx*(2*Bx**2 + Ax*Cx)*(torch.sin(2*mux) + torch.sin(2*mux - 4*nux*torch.pi))))/(8*bx) - (bxs*t*(1.0/torch.sin(2*nux*torch.pi))*((Ax*bx**2 - 2*ax*bx*Bx - 2*Cx + bx*cx*Cx)*torch.cos(2*mux - 2*nux*torch.pi) + 2*(bx*Bx - ax*Cx)*torch.sin(2*mux - 2*nux*torch.pi)))/(2*bx)\n",
    "    return (100*(bxf - bxs)/bxs).real\n",
    "\n",
    "def bby(Ay, By, Cy, nuy, ay, by, ays, bys, muy):\n",
    "    t = 1.0\n",
    "    cy = (1 + ay**2)/by\n",
    "    byf = bys + (bys*t**2*(1.0/torch.sin(2*nuy*torch.pi)**2)*(4*ay*By*Cy*(torch.cos(2*muy) + torch.cos(2*muy - 4*nuy*torch.pi)) - 2*cy*Cy**2*(torch.cos(2*muy) + torch.cos(2*muy - 4*nuy*torch.pi)) + Ay**2*by**3*(1 + torch.cos(2*muy) + torch.cos(2*muy - 4*nuy*torch.pi)) - 4*ay*Ay*by**2*By*(1 + torch.cos(2*muy) + torch.cos(2*muy - 4*nuy*torch.pi)) - 4*ay*by*By*cy*Cy*(1 + torch.cos(2*muy) + torch.cos(2*muy - 4*nuy*torch.pi)) + by*cy**2*Cy**2*(1 + torch.cos(2*muy) + torch.cos(2*muy - 4*nuy*torch.pi)) + 2*by**2*cy*(2*By**2 + Ay*Cy)*(1 + torch.cos(2*muy) + torch.cos(2*muy - 4*nuy*torch.pi)) - 2*by*(2*Ay*Cy + (2*By**2 + Ay*Cy)*torch.cos(2*muy) + (2*By**2 + Ay*Cy)*torch.cos(2*muy - 4*nuy*torch.pi)) - 4*By*Cy*torch.sin(2*muy) - 4*By*Cy*torch.sin(2*muy - 4*nuy*torch.pi) + 2*Ay*by**2*By*(torch.sin(2*muy) + torch.sin(2*muy - 4*nuy*torch.pi)) + 6*by*By*cy*Cy*(torch.sin(2*muy) + torch.sin(2*muy - 4*nuy*torch.pi)) - 2*ay*cy*Cy**2*(torch.sin(2*muy) + torch.sin(2*muy - 4*nuy*torch.pi)) - 2*ay*by*(2*By**2 + Ay*Cy)*(torch.sin(2*muy) + torch.sin(2*muy - 4*nuy*torch.pi))))/(8*by) - (bys*t*(1.0/torch.sin(2*nuy*torch.pi))*((Ay*by**2 - 2*ay*by*By - 2*Cy + by*cy*Cy)*torch.cos(2*muy - 2*nuy*torch.pi) + 2*(by*By - ay*Cy)*torch.sin(2*muy - 2*nuy*torch.pi)))/(2*by)\n",
    "    return (100*(byf - bys)/bys).real\n",
    "\n",
    "plt.figure(figsize=(12, 4))\n",
    "plt.errorbar(ring.locations().cpu().numpy(), 100*((bx - bx_id)/bx).cpu().numpy(), fmt='-', marker='x', ms=0, color='blue', alpha=0.75, label='exact')\n",
    "plt.errorbar(ring.locations().cpu().numpy(), -bbx(Ax, Bx, Cx, nux, axi, bxi, ax, bx, mux), fmt=' ', marker='x', ms=5, color='red', alpha=0.75, label='approximate')\n",
    "\n",
    "plt.gca().tick_params(axis='x', length=6, width=1.5, direction='in', labelsize=12, bottom=True, top=False, labelbottom=True, labeltop=False)\n",
    "plt.gca().tick_params(axis='y', length=0, width=0, labelsize=12)\n",
    "plt.gca().set_xlabel(r'$s$', fontsize=18)\n",
    "plt.gca().set_ylabel(r'$\\Delta \\beta_x / \\beta_x$', fontsize=18)\n",
    "plt.legend(loc='upper left', frameon=False, fontsize=14, ncol=2)\n",
    "plt.ylim(-1, 1)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(100*(((bx - bx_id)/bx)**2).mean().sqrt())\n",
    "print((bbx(Ax, Bx, Cx, nux, axi, bxi, ax, bx, mux)**2).mean().sqrt())\n",
    "\n",
    "plt.figure(figsize=(12, 4))\n",
    "plt.errorbar(ring.locations().cpu().numpy(), 100*((by - by_id)/by).cpu().numpy(), fmt='-', marker='x', ms=0, color='blue', alpha=0.75, label='exact')\n",
    "plt.errorbar(ring.locations().cpu().numpy(), -bby(Ay, By, Cy, nuy, ayi, byi, ay, by, muy), fmt=' ', marker='x', ms=5, color='red', alpha=0.75, label='approximate')\n",
    "plt.gca().tick_params(axis='x', length=6, width=1.5, direction='in', labelsize=12, bottom=True, top=False, labelbottom=True, labeltop=False)\n",
    "plt.gca().tick_params(axis='y', length=0, width=0, labelsize=12)\n",
    "plt.gca().set_xlabel(r'$s$', fontsize=18)\n",
    "plt.gca().set_ylabel(r'$\\Delta \\beta_y / \\beta_y$', fontsize=18)\n",
    "plt.legend(loc='upper left', frameon=False, fontsize=14, ncol=2)\n",
    "plt.ylim(-1, 1)\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "print(100*(((by - by_id)/by)**2).mean().sqrt())\n",
    "print((bby(Ay, By, Cy, nuy, ayi, byi, ay, by, muy)**2).mean().sqrt())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "7ff84266-00a1-410f-9518-0f3374e8af9e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tensor(4.6202e-12, dtype=torch.float64)\n",
      "tensor(2.9309e-08, dtype=torch.float64)\n",
      "\n",
      "tensor(4.3448e-11, dtype=torch.float64)\n",
      "tensor(0., dtype=torch.float64)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Analytical RMS beta-beting\n",
    "\n",
    "def rmsx(Ax, Bx, Cx, ax, bx, nux):\n",
    "    cx = (1 + ax**2)/bx\n",
    "    return 100*(((1.0/torch.sin(2*nux*torch.pi)**4)*((Ax**4*bx**6 + 16*bx**4*Bx**4*cx**2 - 4*Ax**3*bx**4*(2*ax*bx*Bx + Cx - bx*cx*Cx) - 8*bx**3*Bx**2*cx*(4*Bx**2 + 4*ax*Bx*cx*Cx - cx**2*Cx**2) + 4*(1 + ax**2)*Cx**2*(-4 + 4*Bx**2 + cx**2*Cx**2) - 4*bx*cx*Cx**2*(-4 + 4*(2 + ax**2)*Bx**2 + cx**2*Cx**2) + 2*Ax**2*bx**2*(4*bx**3*Bx**2*cx + 2*(1 + ax**2)*Cx**2 + 2*bx*Cx*(4*ax*Bx - 3*cx*Cx) + bx**2*(-2 + (-2 + 8*ax**2)*Bx**2 - 12*ax*Bx*cx*Cx + 3*cx**2*Cx**2)) + bx**2*(16*(1 + ax**2)*Bx**4 + 16*ax*Bx**3*cx*Cx + cx**2*Cx**2*(-4 + cx**2*Cx**2) - 8*ax*Bx*cx*Cx*(-2 + cx**2*Cx**2) + 4*Bx**2*(-4 + 3*cx**2*Cx**2 + 4*ax**2*(-1 + cx**2*Cx**2))) + 4*Ax*bx*(2*ax**2*Cx*(-2*bx*Bx**2 + 4*bx**2*Bx**2*cx + cx*Cx**2) + 2*ax*bx*Bx*(-4*bx**2*Bx**2*cx + 2*cx*Cx**2 + bx*(2 + 2*Bx**2 - 3*cx**2*Cx**2)) + Cx*(4*bx**3*Bx**2*cx**2 + 2*cx*Cx**2 + bx*(4 - 3*cx**2*Cx**2) + bx**2*cx*(-2 - 2*Bx**2 + cx**2*Cx**2))))*torch.cos(4*nux*torch.pi) + 2*(Ax**4*bx**6 + 8*bx**2*Bx**2 + 8*ax**2*bx**2*Bx**2 + 8*bx**2*Bx**4 + 8*ax**2*bx**2*Bx**4 - 16*bx**3*Bx**4*cx + 16*bx**4*Bx**4*cx**2 - 8*ax*bx**2*Bx*cx*Cx + 8*ax*bx**2*Bx**3*cx*Cx - 32*ax*bx**3*Bx**3*cx**2*Cx + 8*Cx**2 + 8*ax**2*Cx**2 + 8*Bx**2*Cx**2 + 8*ax**2*Bx**2*Cx**2 - 8*bx*cx*Cx**2 - 16*bx*Bx**2*cx*Cx**2 - 8*ax**2*bx*Bx**2*cx*Cx**2 + 2*bx**2*cx**2*Cx**2 + 6*bx**2*Bx**2*cx**2*Cx**2 + 16*ax**2*bx**2*Bx**2*cx**2*Cx**2 + 8*bx**3*Bx**2*cx**3*Cx**2 - 8*ax*bx**2*Bx*cx**3*Cx**3 + 2*cx**2*Cx**4 + 2*ax**2*cx**2*Cx**4 - 2*bx*cx**3*Cx**4 + bx**2*cx**4*Cx**4 - 2*Ax**3*bx**4*(4*ax*bx*Bx + (3 - 2*bx*cx)*Cx) + 2*Ax**2*bx**2*(4*bx**3*Bx**2*cx + (5 + ax**2)*Cx**2 + bx*Cx*(12*ax*Bx - 7*cx*Cx) + bx**2*(1 + (-1 + 8*ax**2)*Bx**2 - 12*ax*Bx*cx*Cx + 3*cx**2*Cx**2)) + 2*Ax*bx*(2*ax**2*Cx*(-2*bx*Bx**2 + 8*bx**2*Bx**2*cx + cx*Cx**2) - 4*ax*bx*Bx*(bx - bx*Bx**2 + 4*bx**2*Bx**2*cx - 3*cx*Cx**2 + 3*bx*cx**2*Cx**2) + Cx*(8*bx**3*Bx**2*cx**2 + 2*cx*Cx**2 - bx*(4 + 5*cx**2*Cx**2) + 2*bx**2*(cx - 5*Bx**2*cx + cx**3*Cx**2))) - 2*(Ax**3*bx**5 + Ax**2*bx**3*(-6*ax*bx*Bx + (-4 + 3*bx*cx)*Cx) + Ax*bx*(4*bx**3*Bx**2*cx + 4*ax*bx*Bx*(2 - 3*bx*cx)*Cx + 4*Cx**2 - 8*bx*cx*Cx**2 + 3*bx**2*cx**2*Cx**2 + 4*ax**2*(2*bx**2*Bx**2 + Cx**2)) + cx*(4*bx**3*Bx**2*cx*Cx + 4*Cx**3 - 4*bx*cx*Cx**3 + bx**2*cx**2*Cx**3 - 2*ax*bx**2*Bx*(4*bx*Bx**2 + 3*cx*Cx**2) + 4*ax**2*(2*bx**2*Bx**2*Cx + Cx**3)))*torch.sin(4*nux*torch.pi))))/(64*bx**2)).sqrt()\n",
    "\n",
    "def rmsy(Ay, By, Cy, ay, by, nuy):\n",
    "    cy = (1 + ay**2)/by\n",
    "    return 100*(((1.0/torch.sin(2*nuy*torch.pi)**4)*((Ay**4*by**6 + 16*by**4*By**4*cy**2 - 4*Ay**3*by**4*(2*ay*by*By + Cy - by*cy*Cy) - 8*by**3*By**2*cy*(4*By**2 + 4*ay*By*cy*Cy - cy**2*Cy**2) + 4*(1 + ay**2)*Cy**2*(-4 + 4*By**2 + cy**2*Cy**2) - 4*by*cy*Cy**2*(-4 + 4*(2 + ay**2)*By**2 + cy**2*Cy**2) + 2*Ay**2*by**2*(4*by**3*By**2*cy + 2*(1 + ay**2)*Cy**2 + 2*by*Cy*(4*ay*By - 3*cy*Cy) + by**2*(-2 + (-2 + 8*ay**2)*By**2 - 12*ay*By*cy*Cy + 3*cy**2*Cy**2)) + by**2*(16*(1 + ay**2)*By**4 + 16*ay*By**3*cy*Cy + cy**2*Cy**2*(-4 + cy**2*Cy**2) - 8*ay*By*cy*Cy*(-2 + cy**2*Cy**2) + 4*By**2*(-4 + 3*cy**2*Cy**2 + 4*ay**2*(-1 + cy**2*Cy**2))) + 4*Ay*by*(2*ay**2*Cy*(-2*by*By**2 + 4*by**2*By**2*cy + cy*Cy**2) + 2*ay*by*By*(-4*by**2*By**2*cy + 2*cy*Cy**2 + by*(2 + 2*By**2 - 3*cy**2*Cy**2)) + Cy*(4*by**3*By**2*cy**2 + 2*cy*Cy**2 + by*(4 - 3*cy**2*Cy**2) + by**2*cy*(-2 - 2*By**2 + cy**2*Cy**2))))*torch.cos(4*nuy*torch.pi) + 2*(Ay**4*by**6 + 8*by**2*By**2 + 8*ay**2*by**2*By**2 + 8*by**2*By**4 + 8*ay**2*by**2*By**4 - 16*by**3*By**4*cy + 16*by**4*By**4*cy**2 - 8*ay*by**2*By*cy*Cy + 8*ay*by**2*By**3*cy*Cy - 32*ay*by**3*By**3*cy**2*Cy + 8*Cy**2 + 8*ay**2*Cy**2 + 8*By**2*Cy**2 + 8*ay**2*By**2*Cy**2 - 8*by*cy*Cy**2 - 16*by*By**2*cy*Cy**2 - 8*ay**2*by*By**2*cy*Cy**2 + 2*by**2*cy**2*Cy**2 + 6*by**2*By**2*cy**2*Cy**2 + 16*ay**2*by**2*By**2*cy**2*Cy**2 + 8*by**3*By**2*cy**3*Cy**2 - 8*ay*by**2*By*cy**3*Cy**3 + 2*cy**2*Cy**4 + 2*ay**2*cy**2*Cy**4 - 2*by*cy**3*Cy**4 + by**2*cy**4*Cy**4 - 2*Ay**3*by**4*(4*ay*by*By + (3 - 2*by*cy)*Cy) + 2*Ay**2*by**2*(4*by**3*By**2*cy + (5 + ay**2)*Cy**2 + by*Cy*(12*ay*By - 7*cy*Cy) + by**2*(1 + (-1 + 8*ay**2)*By**2 - 12*ay*By*cy*Cy + 3*cy**2*Cy**2)) + 2*Ay*by*(2*ay**2*Cy*(-2*by*By**2 + 8*by**2*By**2*cy + cy*Cy**2) - 4*ay*by*By*(by - by*By**2 + 4*by**2*By**2*cy - 3*cy*Cy**2 + 3*by*cy**2*Cy**2) + Cy*(8*by**3*By**2*cy**2 + 2*cy*Cy**2 - by*(4 + 5*cy**2*Cy**2) + 2*by**2*(cy - 5*By**2*cy + cy**3*Cy**2))) - 2*(Ay**3*by**5 + Ay**2*by**3*(-6*ay*by*By + (-4 + 3*by*cy)*Cy) + Ay*by*(4*by**3*By**2*cy + 4*ay*by*By*(2 - 3*by*cy)*Cy + 4*Cy**2 - 8*by*cy*Cy**2 + 3*by**2*cy**2*Cy**2 + 4*ay**2*(2*by**2*By**2 + Cy**2)) + cy*(4*by**3*By**2*cy*Cy + 4*Cy**3 - 4*by*cy*Cy**3 + by**2*cy**2*Cy**3 - 2*ay*by**2*By*(4*by*By**2 + 3*cy*Cy**2) + 4*ay**2*(2*by**2*By**2*Cy + Cy**3)))*torch.sin(4*nuy*torch.pi))))/(64*by**2)).sqrt()\n",
    "\n",
    "print(100*(((bx - bx_id)/bx)**2).mean().sqrt())\n",
    "print(rmsx(Ax, Bx, Cx, axi, bxi, nux).real)\n",
    "print()\n",
    "\n",
    "print(100*(((by - by_id)/by)**2).mean().sqrt())\n",
    "print(rmsy(Ay, By, Cy, ayi, byi, nuy).real)\n",
    "print()"
   ]
  }
 ],
 "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.14.0"
  },
  "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
}