jgalazm/nwogus-spectral-radius.ipynb

## nwogus-spectral-radius.ipynb
{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "import numpy as np"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "dx = 1\n",
        "dt = 0.6\n",
        "mu = 1\n",
        "z = -0.531\n",
        "N = 100"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "mu1 =  mu * np.sqrt(-(z*z/2+z+1/3))\n",
        "mu2 = mu*np.sqrt(-z*(0.5*z+1))\n",
        "\n",
        "kappa = np.linspace(0,2*np.pi*(N-1)/N, N)\n",
        "m1 = np.exp(1j*2*kappa)+np.exp(-1j*kappa)-(np.exp(1j*kappa)-1)\n",
        "m1 = m1 *(1-np.exp(-1j*kappa))\n",
        "m2 = 4*np.sin(kappa/2)*np.sin(kappa/2)\n",
        "m3 = (1-np.exp(1j*kappa))*(1-np.exp(-1j*kappa))\n",
        "m4 = (np.exp(1j*kappa)+1)*(1-np.exp(-1j*kappa))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "g11 = 1 - dt/(2*dx) * m3\n",
        "g12 = -dt/(2*dx) * m4 + dt/(2*dx*dx)*mu1*mu1*m1\n",
        "g21 = (-dt/(2*dx) * m4) / (1 + mu2*mu2/(dx*dx) * m2)\n",
        "g22 = 1 - (dt/(2*dx) * m3 ) / (1 + mu2*mu2/(dx*dx) * m2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "G = np.array([[[g11[i], g12[i]],[g21[i], g22[i]]] for i in range(N)])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "from IPython.display import display, Latex"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "rhos = np.zeros((N,))\n",
        "for i in range(N):\n",
        "    GG = np.matmul(G[i], np.transpose(np.conjugate(G[i])))\n",
        "    eigenvalues = np.linalg.eigvals(GG)\n",
        "    spectral_radius = np.max(np.abs(eigenvalues))    \n",
        "    # rhos[i] = np.sqrt(spectral_radius)\n",
        "    rhos[i] = np.abs(eigenvalues).max()\n",
        "\n",
        "    # latex_code = r'$\\begin{bmatrix} %f & %f \\\\ %f & %f \\end{bmatrix}$' % (\n",
        "    #     GG[0, 0], GG[0, 1], GG[1, 0], GG[1, 1]\n",
        "    # )\n",
        "    # display(Latex(latex_code))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "Text(0.5, 0, '$\\\\kappa$')"
            ]
          },
          "execution_count": 9,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAHPCAYAAACBRNrVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABeQ0lEQVR4nO3deVzUdf4H8NccMMM5KMgpyOGBioInouZJWZplp2W7HpVpWdbadriVbrvbWvZbszY3y9KsNM3KDi1NySOVREEOFVG8QI4BVGa4j5nP74+BKRKU+zvH6/l4zB9+5/udefMVndd8TpkQQoCIiIhIInKpCyAiIiL7xjBCREREkmIYISIiIkkxjBAREZGkGEaIiIhIUgwjREREJCmGESIiIpIUwwgRERFJimGEiIiIJMUwQkRERJJiGCEiIiJJKaUugMhWCSEglzPvN4ZbYhHR78m4UR5Rx0hISMC+ffvw3HPPSV0KEZFF49c2og6ybds23H777VKXQURk8RhGiDpIeno6+vbtK3UZREQWj2GEqAPk5uYiICBA6jKIiKwCwwhRB2ivLprly5cjPDwcRqOxHarqHKtXr0ZQUBCqqqpadN2RI0cwcuRIuLi4QCaTITk5uWMKJCKLwzBC1AH27duHsWPHtuk19Ho93njjDbzwwguNzso5e/YsnnzySfTu3RtqtRouLi6IjIzEX//6V2RmZjb6mq25pt7evXsxe/bsG9Y9e/ZsVFdX4/3332/WzwkANTU1uO+++3DlyhW89dZb+PTTT9GjR49mX09E1o1Te4naWWVlJZRKJRwcHNr0OmvXrkVtbS0efPDBa57bsGED5s6dCxcXFzzwwAOIiIiA0WjEsWPH8NFHH2H16tXQ6/UNQkxrrtHpdEhPT8eIESMavH9xcTEyMjIQHR19TW1qtRqzZs3CihUr8NRTT0Emk93wZz179iwuXryINWvW4NFHH23JbSIiG8AwQtTO4uLiMGHChDa/zrp163DHHXdArVY3OP7dd99h1qxZGD9+PDZt2gRPT88Gz7/xxht4++23G4SK1lwDABcvXsSsWbNw66234uabbwYAfPXVV3juueewcOHCRsMIANx///1Yvnw59uzZ06x7UVBQAADw8PC44bnNVVZWBhcXl3Z7PSLqQIKI2tWCBQtEQUHBdc/5+eefxeDBg4Wzs7NwdnYW999/vyguLjY/f+7cOQFAfPzxxw2uKy4uFt7e3qJv376itLS0WfW05prfq6qqEm+++abw9/cXLi4u4v777xcXLly44XVdu3YVCxcuvOF5s2bNEgAaPMaOHWt+PikpSdx6663Czc1NuLi4iAkTJoj4+PgGr7F06VIBQJw4cUI8+OCDwsPDQ0RFRd3wvUNCQsRDDz10zfFx48aJMWPG3PD6toiNjRUjRowQhw4dEmPHjhXOzs4iLCxMbN++XQghxPbt20V0dLRwdnYWkZGR4ujRow2u9/LyEgsWLLjmdYcMGSImT57cobUTtTeOGSFqZ4WFhejWrVuTz3/44YeIjY1F37598eabb2LKlCn44osv8Mwzz5jPOXToEABg8ODBDa5duXIlCgoKsHz58mZ/62/NNb8nk8kgl8vN3S0ymaxZXS+DBw/GwYMHb3jevHnz8Le//Q0AsHDhQnz66ad46aWXAAAnTpzATTfdhJSUFDz//PN45ZVXcP78eYwbNw6HDx++5rXuu+8+lJeX49///jfmzp173fctLS3FhQsXEBkZec1zqampGDhw4DXHa2pqUFRU1KzHjQYdp6amQqfTYfr06Rg/fjz+/e9/o6SkBDNmzMD777+PhQsX4q677sLLL7+Ms2fP4uGHHzZfm5ubi6KiomtqNxgMOHHiRKO1E1k0qdMQkTVKTU0VK1euvOZ4cnKyeO2115q8Lj09XSiVSvHOO+80OD5mzBihVqtFTU2NEEKIl19+WQAQJSUlDc4LDw8X/v7+wmAwNDiu0+lEYWGh+fH7FpDWXPP7nzM8PFw89dRT4vvvvxezZs0SW7ZsESEhIY3+/L/32GOPCScnp+ueU2/Pnj0CgNiyZUuD49OmTROOjo7i7Nmz5mO5ubnCzc2tQctFfcvIgw8+2Kz3E0KI+Ph4AUDs3LmzwfHs7GwBQHzwwQdN1tmcx/nz55t8b61WKwAIHx8fkZubaz7+zjvvCAAiPDxc6HQ68/FFixYJmUwmKisrhRBC/PjjjwKAOHz4cIPXPX78uAAgNmzY0Oz7QGQJOGaEqIX27duHZ555BoWFhXj66acbPLdt2zZMnTq1yWv//ve/Y+DAgXjyyScbHB8zZgz279+PK1euwNvbG5cvX4ZSqYSrq6v5nIKCApw6dQoPPPDANWM7hg8fjoyMDPOf33zzTfz1r39t1TW/FxQUhHXr1mHEiBHYu3cvAODee+9FbGxsg2sb06VLF1RUVKC8vBzOzs7XPbcxBoMBP/30E6ZNm4bQ0FDzcT8/P8yYMQNr1qyBXq+Hu7u7+bn58+c3+/WPHz8OANe0LqSkpABAo60LkZGR2LVrV7Ne39fXt8nnUlNTAZh+H/z8/MzH6/++33zzzQY/l0ajgVwuN/8dpqamQi6XIyIiotHaBwwY0KwaiSwFwwhRC40dOxZ79uyBj48PkpKSGnSlpKWlmbsY/qi2thY//PADFi9efE03R1lZGWQyWYMPoD+6dOkSACAkJOSa59577z0YDAbExcXh9ddfN3+Qtuaa39NoNNfMpAFMA02bGrxaT9Rte9WcLp3GFBYWory8HH369Lnmub59+8JoNCI7Oxv9+/c3H2/s52xKWloafHx84OPj0+B4Ux/0gClgxcbGtuCnaPq9AeCOO+5ocDwjIwNOTk7mwcL1Tp8+jbCwMPMMrZSUFPTs2fOakJecnAwHBweEh4e3uUaizsQwQtQKHh4eGDduHL777jtzGLnRWJGkpCSUlJQgKirqmueSk5MRGRlpnjnj6emJ2tpalJSUwM3NDQBQXV0NwBRq/mj8+PEATK02wG/f9ltzTVPGjRuHcePGXfec37t69SqcnZ3h5OTU7GvaqiXvdfz48UZ/5uTkZISGhjY6vqa6uhpXrlxp1ut369YNCoWi0edSU1Ph5+cHf3//BsdTUlIQEREBlUp1zfHfh8W0tLRGaz9y5Aj69OnT5mnlRJ2NA1iJWmnq1Kn4/vvvzX/+4YcfMHny5CbPr19R9I8fcnl5eThw4ADuuusu87H6b7bnz583H6vvqjhx4kST7/HHb/utuaa9nD9/vk1783Tr1g3Ozs6NdgedOnUKcrkcgYGBrX79tLQ09O7du8Exo9GIn3/+uckBoIcOHYKfn1+zHtnZ2U2+d2pqaqNhIiUl5ZrjNTU1yMjIMNdkNBqRkZFxzb0tKCjAgQMHOHiVrBJbRohaaerUqVi4cCFycnIQEBCAPXv2YPXq1U2eXz9OYN++fRgzZgwAU4vF448/Do1Gg3nz5pnPjYmJAQAcPXrU/OHi7e2NMWPG4Mcff8SBAwcwevToBq9fU1NzzTfo1lzTXpKSkvDQQw+1+nqFQoFbbrkF3377LS5cuIDg4GAAgFarxcaNGzF69OjrdmtdT0FBAQoLC5GXl9fg+DvvvIOioqImx1y0x5gRg8GAkydPXtMVU1RUhLy8vGvCSHp6Ompqasx/RwaDATU1NSgvLzefU1tbi3nz5qG2tpbjRcgqMYwQtVKPHj0QERGBbdu24eGHH4bRaLxmgbLfS0tLQ79+/fDaa6+hrKwMfn5+2LRpE44cOYIvv/yyQctEaGgoIiIisHv37gZTOletWoWbbroJ48ePx8yZMzFo0CAIIXDmzBls2bIFWq0WCxYsaPC+rbmmrRITE3HlyhXceeedbXqdf/3rX9i1axdGjx6NJ554AkqlEu+//z6qqqqwfPnyVr9u/ZiNn376CU888QTCw8Px66+/YufOneb6Dx8+fM24mPYYM3LmzBlUVlY2OXD2j8frQ2x9GHFwcMDAgQPx3nvvwcnJCU5OTtiyZYu5i4phhKySxLN5iKza3/72NzFlyhSxa9cusWbNmuue26VLF/Hyyy+LDz74QAQGBgqVSiViYmJEXFxco+evWLFCuLq6ivLy8gbHz549K2bPni38/f2FUqkUGo1GDB06VDz//PMiLS2t0ddqzTVt8cILL4igoCBhNBqbdX5TU3uFMC16NmnSJOHq6iqcnZ3F+PHjxaFDhxqcUz+1t7CwsFnv99ZbbwmFQiG2b98uwsLChFqtFjfffLNIS0sTYWFhonv37iIxMbFZr9VSX3zxhQAgjh8/3uD4ihUrBIAGi98JIcTzzz8v3N3dG9zLpKQkMWTIEKFWq0X//v3FBx98ID766CMBQFy8eLFD6ibqSDIh6oa8E1GLHT58GOPGjcPMmTOvmab5e9nZ2QgKCsLGjRsb3WumMTqdDqGhoVi+fDkeeeSR9iy7Q1VVVSE4OBgvvvjiNVOfLcWjjz6K/fv34/Tp01KXQkTgAFaiNhk+fDg0Gg1SUlKaDCLAb90Cv5+GeiMajQbPP/883nzzzRuu5mlJ1q1bBwcHhxat+dHZ6rvMiMgyMIwQtYFMJsOUKVNw2223Xfe81NRUKBSKRtfMuJ4XXnjBPHPEWsyfPx9ZWVnXTE+1FEIInDx5kmGEyIJwACtRG82fP/+Gu82mpaUhLCzMYj+g7cn58+dRWlrKMEJkQThmhIiIiCRlPW2/REREZJMYRoiIiEhSDCNEREQkKasYwGo0GpGbmws3N7dW7wBKREREnUsIgZKSEvj7+193VqBVhJHc3Nw2bYhFRERE0snOzkb37t2bfN4qwkj9FurZ2dmt3hiLiIiIOpder0dgYKD5c7wpVhFG6rtm3N3dGUaIiIiszI2GWHAAKxEREUmKYYSIiIgkxTBCREREkmIYISIiIkkxjBAREZGkGEaIiIhIUgwjREREJCmGESIiIpIUwwgRERFJimGEiIiIJNXiMLJ//35MnToV/v7+kMlk+Oabb254zd69ezF48GCoVCr07NkTH3/8cStKJSIiIlvU4jBSVlaGyMhIrFq1qlnnnz9/HlOmTMH48eORnJyMZ555Bo8++ih27tzZ4mKJiIjI9rR4o7zbbrsNt912W7PPX716NUJCQvCf//wHANC3b18cOHAAb731FiZNmtTSt+9QQghcLa+Bg0IGN7WD1OUQERF1uKpaA7Iul8PPwwmuKmn2z+3wd42Pj0dsbGyDY5MmTcIzzzzT5DVVVVWoqqoy/1mv13dIbZ/9ehEJ568gX1eJfL3pUV1rBAB4uToi2NMFwV4uCOvmiikD/BDk6dwhdRAREXUGXUUNtqXm4lReCS5cLsO5wjLk6iogBLB29lBMCPeRpK4ODyP5+fnw8Wn4w/n4+ECv16OiogJOTk7XXLNs2TK8+uqrHV0aEs5fwXcpuY0+V1RajaLSahy9eBUAsHznKcT29cHDo0IwIrTrDbdDJiIishTnCkvx8aEL+DLxEsqrDdc876pSQl9RK0FlJtK0x9zA4sWLsWjRIvOf9Xo9AgMD2/19pkb6Y0CABr4atenhroa3uwpVtUZcLCrH+ctluFBUhiMXruCXM0XYdVKLXSe1CPd1w+PjwnBHpD9DCRERWaxjWVfxTtwZ7MkoNB/r4+OGCX29EeLpgpBuLgj2dIGXq6Okn2cdHkZ8fX2h1WobHNNqtXB3d2+0VQQAVCoVVCpVR5eGm/s13hylUiowoLsGA7przMcyC0qw7uAFfJ2Ug1P5JXh6UzJ2pxdg2d0DJOtjIyIiaozRKLB6/1n856fTMBgFZDJgYrg35owKwcgwT4v7It3hn6IxMTH44YcfGhzbtWsXYmJiOvqt21VPbze8dtcAPDepDz4+dAH//TkT36fk4niODqtmDEY/f3epSyQiIsKVsmos+iIZe+taQ6ZG+uPZm3sj2MtF4sqa1uKpvaWlpUhOTkZycjIA09Td5ORkZGVlATB1scycOdN8/vz583Hu3Dk8//zzOHXqFP73v//hiy++wF/+8pf2+Qk6mYezI56J7Y0v5o2Av0aN80VlmPa/g9h4OAtCCKnLIyIiO3b0whVMeecX7M0ohEopxxv3DMA7D0RZdBABWhFGjh49ikGDBmHQoEEAgEWLFmHQoEFYsmQJACAvL88cTAAgJCQE27dvx65duxAZGYn//Oc/+PDDDy1uWm9LDenRFdsX3oTxfbqhutaIv21Nw6vfn2QgISIiSXyfkovpH/yKPF0lQr1c8M2CUZg+LMjiumQaIxNW8Omp1+uh0Wig0+ng7m5Z3SFGo8AHv5zDGztOQQhgwfgwPDcpXOqyiIjIjvx8SovHPklErVHg9oF+eP2egRYxnrG5n9/cm6aN5HIZ5o8Nwz/vjAAArNpzFqv3nZW4KiIishfxZy/j8c+SUGsUuDPKH+88MMgigkhLMIy0kz+N6IEXbjW1iLz+4ylsOHxR4oqIiMjWpWQX49H1R1BVa0RsXx/8332RkMstv1vmjxhG2tHj48LwxLgwAMDL3xzHt8k5EldERES2KiO/BLPWJaCs2oCYUE+8O2MQHBTW+bFunVVbsOcm9cGfR/SAEMCzX6Rg3+nCG19ERETUArnFFZi59jCKy2sQFeiBNbOGQu2gkLqsVmMYaWcymQyv3tEfd0b5o9Yo8MRniTiRq5O6LCIishG6ihrMXpcArb4Kvbxd8fGcYVY3RuSPGEY6gFwuw5v3RiIm1BNl1QY8/PER5BZXSF0WERFZuepaI+Z/mojT2lJ4u6nw8cPD4eHsKHVZbcYw0kEclXKs/vMQ9PJ2hVZfhTnrjkBfWSN1WUREZKWEEHjhq1TEn7sMF0cF1s4ehgCPxrdVsTYMIx1I4+SAjx8eDm83FTK0JXj8s0RU1xqlLouIiKzQil2nsfVYDhRyGVY9NBgRAZobX2QlGEY6WICHE9bOHgYXRwUOZl7G05uOobLm2u2biYiIGiOEwPv7zuK/P2cCAP59VwTG9fGWuKr2xTDSCSICNFj10GA4KGT48Xg+/vzRYRSXV0tdFhERWTiDUeDv353Ash9PAQAWTuyF6cOCJK6q/TGMdJJxfbyxfs5wuKmVOHLhKu5+7xCyLpdLXRYREVmoimoD5n2aiPXxpkU0X57SF3+J7SVxVR2DYaQTjezphS/nj4S/Ro1zhWW4+72DSMkulrosIiKyMEWlVXhgza/Yna6Fo1KOVTMG49GbQq1i07vWYBjpZH183bB1wSj083NHUWk1pn8Qz5VaiYjI7HiODtNWmb6sejg7YOOj0Zgy0E/qsjoUw4gEfNzV+GJ+DMb27obKGiOe3pSMpd8e50wbIiI798WRbNz93iFculqBoK7O+PrxkRga3FXqsjocw4hEXFVKrJ09DE+O7wkAWB9/EdM/iOfiaEREdqiyxoAXvkzF81+lorrWiAnh3vj+ydEI7eYqdWmdgmFEQgq5DH+d1AcfzRoKd7USx7KKcft/D+CXM9zPhojIXly8XIZ73juEzUezIZeZ9jj7cOZQaJwdpC6t0zCMWICJfX2w7amb0N/fHVfKqjFzbQKW/ZjObhsiIhu39dglTHnnAE7k6tHVxRGfPByNBeN7Qi63zYGqTWEYsRBBns746vGRmBEdBCGA9/edwz3vHcL5ojKpSyMionZWUlmDZzYdw182p6C0qhbDgrtg21OjMbqXl9SlSUImhBBSF3Ejer0eGo0GOp0O7u7uUpfT4XYcz8MLX6VBV1EDZ0cFXr2jP+4d0t1mp3QREdmTY1lX8fSmZGRdKYdcBjw9sTcWjA+DUmF77QPN/fxmGLFQeboK/GVzMn49dwUAMKm/D/41bQC6uakkroyIiFqjqtaA/8Zl4r19Z2EwCgR4OOHtB6JserYMw4gNMBgFVu87i7d2nUatUaCLswP+OS0Ctw/0l7o0IiJqgeM5Ovx1SwpO5ZcAAKZG+uNf0yKgcbLtQaoMIzbkRK4Of92SivQ8PQBg8gBf/PPOCHi6spWEiMiSVdca8e6eTPxvTyZqjQJdXRzxr2kRmDzAthcxq8cwYmOqa41YtScTq+p+obs4O2Dx5L64j2NJiIgs0tELV/DS1uPI0JpaQyYP8MU/7oyAlx19kWQYsVF/bOobHtIV/74rAj293SSujIiIAKC4vBqv/3gKm45kAwC6ujjiH3f2t8sudoYRG1ZjMGLtgfNYufsMKmoMcFDIMG9MGBaM7wknR4XU5RER2SUhBLYey8Fr29NxuawaADB9aCBevC0cXVwcJa5OGgwjduDS1XIs/fYE4k4VAAD8NWq8OLkvpg70Y9cNEVEnSs4uxj++P4GkrGIAQC9vV7x21wAMD7HdmTLNwTBiJ4QQ2HlCi39uO4mcun1thvTogiW390NkoIe0xRER2bh8XSWW7ziFr4+Zdl93dlTgyQk98ejoUDgqbW/dkJZiGLEzlTUGrNl/Dv/bexYVNQYAwN2DAvCXm3sjsKuzxNUREdmW0qpafPjLOby/75z5/9x7BnfH87f2gY+7WuLqLAfDiJ3S6ivxxo5T+DrJlNIdFDI8FN0DT07oaVcjuImIOkJljQEbDmdh1Z5MXKkbFzK0RxcsmdoPA7t7SFucBWIYsXOpl4qxfEcGDmQWATA1HT46OgSPjgmFu9q2F9khImpvtQYjth7LwcrdZ8xd4qFeLnj2lj6YPMCX4/SawDBCAICDmUVYvuMUUi7pAADuaiUeHh2COaNCbH7lPyKitqoPIav2ZOLC5XIAgK+7Gs/E9sK9Q7rb5H4y7YlhhMxMg1zz8Z+fTuNMQSkAwE2lxJxRwXh4dAg8nO1zyhkRUVNqDEZsTcrBu3sykXXFFEK6ujhi/thQzIwJhtqByyg0B8MIXcNoFPjxeD7eiTtjXhHQVaXEQyOC8MjoEHi7cdAVEdm3yhoDvjiajff3nTN3x3i6OOKxMaH404gecFEpJa7QujCMUJOMRlNLydtxZ8wruToq5bh/aHfMGxPG2TdEZHf0lTX4NP4i1h08j6JS08BUL1dHzBsThodGBMHZkSGkNRhG6IaMRoGfTxXgf3szzQv1KOQy3D7QD/PGhKGfP+81Edm2An0l1h26gM/iL6KkqhYAEODhhPljQ3Hf0EB2x7QRwwg1mxACh89fwao9mfjlTJH5+Nje3TBvbChiQj05UpyIbMq5wlKs+eUcvkrMQbXBCADo6e2KJ8aFYWqkPxw4MLVdNPfzu1V3e9WqVQgODoZarUZ0dDQSEhJueH7fvn3h5OSEPn364JNPPmnN21IHkclkGBHqiU8fica2p0bj9oF+kMuAfacLMWPNYUxbdRA/puXBYLT43EpEdF0p2cWY/2kiJq7Yh88TslFtMGJwkAc++PMQ/PTMGNw9uDuDiARa3DKyefNmzJw5E6tXr0Z0dDRWrlyJLVu2ICMjA97e3tec/9577+GFF17AmjVrMGzYMCQkJGDu3LnYuHEjpk6d2qz3ZMtI58u6XI41v5zDF0ezUVVr+tYQ6uWCeWNDMW1QAFRKNl0SkXUQQuBAZhHe23sWh85eNh+P7euN+WPDMDTYvveP6Ugd1k0THR2NYcOG4d133wUAGI1GBAYG4qmnnsKLL754zfkjR47EqFGj8Oabb5qPPfvsszh8+DAOHDjQrj8Mtb/LpVVYf+gCPj50AfpKU3+qj7sKc28KxYxoDuoiIstVP1h/1d5MHM/RAwCUchnujArAvLGh6O3jJnGFtq+5n98t+iSprq5GYmIiFi9ebD4ml8sRGxuL+Pj4Rq+pqqqCWt1wyqiTkxMSEhJQU1MDBwcuvGXJPF1VWHRLHzw2NgybErKw5pdz0Oqr8K/t6Xhv71k8elMo/hzTA66c7kZEFsJgFNielod3fz6D01rT2kpODgo8MDwQj94UigAPJ4krpD9q0SdIUVERDAYDfHx8Ghz38fHBqVOnGr1m0qRJ+PDDDzFt2jQMHjwYiYmJ+PDDD1FTU4OioiL4+fldc01VVRWqqqrMf9br9S0pkzqAq0ppDh5bk3Lwv71nkXWlHG/sOIX395/Fo6NDMHtUCEMJEUnGYBT4LiUH//05E+cKywCYFnicPSoYc0aFoKsLF3i0VB3+yfHKK68gPz8fI0aMgBACPj4+mDVrFpYvXw65vPFBQsuWLcOrr77a0aVRK6iUCjwwPAj3DumOb5Nz8e6eTJwvKsP//XQaaw9ewILxPfFQdBCnwxFRpxFCYNdJLf7vpwxzS4jGyQGPjA7BrJHB3PrCCrRozEh1dTWcnZ3x5ZdfYtq0aebjs2bNQnFxMb799tsmr62pqYFWq4Wfnx8++OADvPDCCyguLm40kDTWMhIYGMgxIxbIYBTYlpqLt3efwbki0zcRf40az8T2xt2DA7hvAxF1qPizl7F85ykcq1sryV2txLyxYZgZ0wNu3BRUch06gHX48OH473//C8A0gDUoKAhPPvlkowNYGzN27FgEBARg48aNzTqfA1gtX63BiC8TL2Hl7jPI11cCAHp5u2Lp1P4Y3ctL4uqIyNacLyrDP7edxM+nCgCYxoTMGRWMeWPCoHFmCLEUHRZGNm/ejFmzZuH999/H8OHDsXLlSnzxxRc4deoUfHx8sHjxYuTk5JjXEjl9+jQSEhIQHR2Nq1evYsWKFdi1axcSExMRHBzcrj8MSa+yxoBP4y/if3szcbW8BgBwa39fvDSlL5eZJ6I2K62qxbs/Z+KjA+dQYxBQymWYER2EJyf05P5aFqhDZtMAwPTp01FYWIglS5YgPz8fUVFR2LFjh3lQa15eHrKyssznGwwG/Oc//0FGRgYcHBwwfvx4HDp0qNlBhKyL2kGBuWNCcf/QQLy1+zQ+/fUidpzIx56MAswbE4onxvfkeBIiajEhBL5JzsGyH06hoMTUjT+2dzcsmdoPYd1cJa6O2orLwVOHysgvwavfnzAvNBTi5YLl9w7EMC4yRETNlFNcgRe/SjVvVxHs6YxXbu+HCeHe3KrCwnFvGrIYQgj8eDwfr35/Alp9FWQyYFZMMJ6/tQ8XTSOiJhmNAhsTsrDsh3SUVRugUsqxcGIvPHpTCFeBthIMI2RxdBU1eG37SXxx9BIAILCrE5bfE4mYME+JKyMiS5N9pRzPf5mK+HOmVtWhPbpg+b0DEcouGavCMEIWa//pQiz+Og05xRWQyYCnxvfE07G9oZCzuZWIgB/T8vD8l6koqaqF2kGO5yeFY9bIYP4fYYUYRsiilVTW4LXt6dh0JBsAMCK0K955YBC83TkansheVdUasOyHU/j40AUAwJAeXbDi/kj08HSRtjBqNYYRsgrfJufgb1+noazaAC9XR7z9wCCM6sl1SYjsTdblcizYmIS0HB0AYP7YMDx7S284cOFEq9bcz2/+LZOk7owKwHdPjUa4rxuKSqvxp48OY9WeTFhBRiaidrI3owBT/vsL0nJ08HB2wLrZw/DibeEMInaEf9MkubBurvhmwSg8ODwQQgBv7szAC1+losZglLo0Iupgn/16EY+sP4qSyloMDvLADwtvwvhwb6nLok7GMEIWQe2gwLK7B+Ifd/aHXAZ8cfQSZq9LgK6iRurSiKgDGI0Cr20/iZe/OQ6DUeCewd2x6bEY+Hs4SV0aSYBhhCzKzJhgfDRrGFwcFTiYeRn3vHcI2VfKpS6LiNpRRbUBj29IxJpfzgMAnr25N/7vvoFwVPIjyV7xb54szvhwb3wxPwa+7mpkFpTirv8dxKl8vdRlEVE70FXUYMaHv2LnCS0cFXK8/UAUnprYiyup2jmGEbJI/f01+GbBKPT1c0dRaTUe/OBXnMxlICGyZrryGvz5o8M4llUMD2cHbJgbjTujAqQuiywAwwhZLF+NGpvmjkBkdw2ulpu+TR2vm/ZHRNblalk1Znz4K1Iv6dDVxREbHx3BParIjGGELJrG2QGfPhqNQUEeKC6vwYw1vyL1UrHUZRFRC1wpq8aMDw/jRK4eni6O+HzuCPTz55pR9BuGEbJ47moHfPLwcAzp0QX6ylo89OFhHMu6KnVZRNQMl0urMGPNr0jP08PLVYVNj41AH183qcsiC8MwQlbBTe2A9Q8Px/DgriiprMWcj48gs6BU6rKI6DrKqkz/Vk/ll8DbzRREevkwiNC1GEbIariqlPj44WGICjR12cxam4ACfaXUZRFRI2oNRjy5MQmpl3To4uyAzx8bgZ7e3HGXGscwQlbF2VGJj2YNRbCnM3KKKzB73RGUVHJhNCJLIoTAS1uPY09GIdQOcnw0exjCujGIUNMYRsjqeLqqsP7h4fBydcTJPD0e/ywJ1bVcOp7IUqzcfQabj2ZDLgPefXAwBgd1kboksnAMI2SVeni6YO3sYXB2VOBAZhFe+CqVm+sRWYDPE7LwdtwZAMC/pg1AbD8fiSsia8AwQlZrYHcPrHpoMBRyGbYey8H/9p6VuiQiu/bruct4+ZvjAICFE3piRnSQxBWRtWAYIas2vo83/nFnfwDAf37KwC9nCiWuiMg+5esq8eTGJBiMAndG+eMvN/eWuiSyIgwjZPVmDA/C/UO7wyiAhZ8fw6Wr3FiPqDNV1xrxxIZEFJVWI9zXDa/fPZB7zVCLMIyQ1ZPJZPjHnREYEGBaNv7xz5JQWWOQuiwiu/Gv7SeRlFUMN7US7/95CJwcFVKXRFaGYYRsgtpBgff+NBhdnB2QlqPD0m9PSF0SkV34OukSPom/CABYOT0KPTxdJK6IrBHDCNmM7l2c8c6DgyCXAZuPZmNTQpbUJRHZtJO5eiz+Og0AsHBiL0zsy5kz1DoMI2RTburVDc/e0gcA8PfvT3DJeKIOUlljwMJNx1BVa8S4Pt3wzMReUpdEVoxhhGzO42PDcFMvL1TWGLHoi2TUGLggGlF7e2PHKWQWlKKbmwor7o+CXM4Bq9R6DCNkc+RyGd68NxIaJwekXtLhvz9nSl0SkU05cKYI6w5eAAC8ee9AdHVxlLYgsnoMI2STfDVq/GtaBABg1Z5MJGVdlbgiItugK6/BX7ekAAD+PKIHxvXxlrgisgUMI2Szpkb6484ofxiMAos2J6O8ulbqkois3ivfHke+vhKhXi5YPDlc6nLIRjCMkE37xx0R8NOoceFyOV7bni51OURW7buUXHyXkguFXIYV06Pg7KiUuiSyEQwjZNM0zg74v/siAQAbDmdh32kuF0/UGgX6SrxSt+/MUxN6IirQQ9qCyKYwjJDNG9XTC7NHBgMAXv4mDRXVXJ2VqKX+se0kdBU1GBCgwYLxPaUuh2wMwwjZhecm9YGfRo3sKxX4789npC6HyKrsO12Ibal5kMuAZXcPgIOCHx3UvvgbRXbBRaXE3+8w7e77wf5zOK0tkbgiIutQWWMwd8/MGRWCiACNxBWRLWIYIbsxqb8vYvv6oNYo8NLWNBiNQuqSiCzef38+g6wr5fDTqLHo5t5Sl0M2qlVhZNWqVQgODoZarUZ0dDQSEhKue/6GDRsQGRkJZ2dn+Pn54eGHH8bly5dbVTBRW7x6Z384Oypw5MJVbEnMlrocIot2RluCD/afAwD8/Y7+cFFx9gx1jBaHkc2bN2PRokVYunQpkpKSEBkZiUmTJqGgoKDR8w8ePIiZM2fikUcewYkTJ7BlyxYkJCRg7ty5bS6eqKUCPJzwl1jTt7tlP57C5dIqiSsiskxGo8BLW4+jxiAQ29cHk/r7Sl0S2bAWh5EVK1Zg7ty5mDNnDvr164fVq1fD2dkZa9eubfT8+Ph4BAcHY+HChQgJCcHo0aMxb968G7amEHWUOaOC0dfPHcXlNXjtB649QtSYL5MuIeHCFTg7KvDqnf2lLodsXIvCSHV1NRITExEbG/vbC8jliI2NRXx8fKPXxMTEIDs7Gz/88AOEENBqtfjyyy8xefLktlVO1EpKhRz/visCMhnwdVIOkrOLpS6JyKKUVNZg+Y5TAIBnYnshwMNJ4orI1rUojBQVFcFgMMDHx6fBcR8fH+Tn5zd6zahRo7BhwwZMnz4djo6O8PX1hUajwapVq5p8n6qqKuj1+gYPovY0KKgL7h7UHQDwr20nIQQHsxLVW73vLIpKqxHq5YI5o0KkLofsQIfPpjl58iSefvppLFmyBImJidixYwcuXLiA+fPnN3nNsmXLoNFozI/AwMCOLpPs0HOT+kDtIMfRi1ex43jjYZrI3uQUV+DDX84DAF68LZxrilCnaNFvmZeXFxQKBbRabYPjWq0Wvr6ND25atmwZRo0aheeeew4DBw7EpEmT8L///Q9r165FXl5eo9csXrwYOp3O/MjO5qwHan++GjUeuykUAPD6jlOoquXKrERv7jiFqlojokO64uZ+Pje+gKgdtCiMODo6YsiQIYiLizMfMxqNiIuLQ0xMTKPXlJeXQy5v+DYKhQIAmmwaV6lUcHd3b/Ag6gjzxoahm5sKFy+X49P4i1KXQySplOxifJOcCwB4eUo/yGQyiSsie9Hi9rdFixZhzZo1WL9+PdLT0/H444+jrKwMc+bMAWBq1Zg5c6b5/KlTp+Lrr7/Ge++9h3PnzuHgwYNYuHAhhg8fDn9///b7SYhawUWlxF9vMU31fSfuDK6WVUtcEZE0hBDmna3vHhSAAd250ip1nhavYDN9+nQUFhZiyZIlyM/PR1RUFHbs2GEe1JqXl4esrCzz+bNnz0ZJSQneffddPPvss/Dw8MCECRPwxhtvtN9PQdQG9w4JxLqDF3AqvwTv/HwGS6dyGiPZn50n8pFw4QrUDnL8dVIfqcshOyMTVjCNQK/XQ6PRQKfTscuGOsSBM0X400eHoZTLsGvRWIR4uUhdElGnqa414pa39uHC5XI8NaEnnr2FYYTaR3M/vzlMmgjA6F5eGN+nG2qNAit2nZa6HKJOtfloNi5cLoeXqwrzxoZJXQ7ZIYYRojrP3xoOANiWmouMfO7qS/ahssaAVT9nAgAWTuwJV+4/QxJgGCGq09fPHVMG+EEIYOVuto6QfdiUkIV8fSX8NGpMH8Y1nUgaDCNEv/N0bC/IZMCPx/NxIlcndTlEHaqyxoBVe88CAJ6c0BMqpULiisheMYwQ/U5vHzdMHWiacr5y9xmJqyHqWJ/9ehGFJVUI8HDCfUPYKkLSYRgh+oOFE3tBLgN2ndQi7RJbR8g2lVfX4r26VpGFE3vCUcmPA5IOf/uI/qCntyumRQUAAFbsypC4GqKO8Un8RVwuq0ZQV2fcPbi71OWQnWMYIWrEUxN7QSGXYU9GIZKyrkpdDlG7Kq2qxfv7TK0iT0/sxc3wSHL8DSRqRIiXC+4eZGodeYvrjpCNWX/oAq6W1yDUywV3RnFbDpIewwhRE56a0AtKuQy/nClCcnax1OUQtYvy6lqs+eUcANP4KCVbRcgC8LeQqAlBns64o+5b4+q6gX5E1m7zkWwUl9egh6czpkayVYQsA8MI0XXMr1sae+fJfJwtLJW4GqK2qTEY8eEv5wEAj40JhUIuk7giIhOGEaLr6O3jhti+3hACWLP/nNTlELXJ9ym5yCmugJerCvdwBg1ZEIYRohuobx35OikHWn2lxNUQtY4QAu/vMwXqOaOCoXbgaqtkORhGiG5gaHBXDAvugmqDEWsPnJe6HKJW2ZNRgAxtCVxVSvxpRA+pyyFqgGGEqBnqW0c2HM6CrqJG4mqIWq5+tdUZ0UHQODlIXA1RQwwjRM0wvo83evu4orSqFhsOX5S6HKIWSbx4BUcuXIWjQo5HRodIXQ7RNRhGiJpBLpeZW0fWHriAyhqDxBURNd97e01jRe4aFAAfd7XE1RBdi2GEqJmmRvojwMMJRaVV+DopR+pyiJrljLYEu9O1kMmAx8aGSl0OUaMYRoiayeF3TdwfHTgHIYTEFRHd2NqDFwAAN/f1QVg3V2mLIWoCwwhRC9w3tDtcVUqcLSzDgcwiqcshuq7i8mpsPXYJADhWhCwawwhRC7ipHXDvENNiUevqvnESWapNR7JRWWNEXz93DA/pKnU5RE1iGCFqodkjgyGTAT+fKsD5ojKpyyFqVK3BiE/jTTO/5owKhkzGpd/JcjGMELVQsJcLxvfxBgB8En9B2mKImrA7XYuc4gp0dXHEHdwQjywcwwhRK8weGQwA2HL0EkoquQgaWZ76gasPDg/k0u9k8RhGiFrhpl5eCOvmgtKqWnyVeEnqcogaOJGrQ8L5K1DIZVz6nawCwwhRK8hkMsweZZqdsD7+IoxGTvMly7H+0AUAwG0RvvDTOElbDFEzMIwQtdLdgwLgplbifFEZ9p0ulLocIgDA5dIqfJOcC8A0cJXIGjCMELWSi0qJ6UMDAQDr6r6JEklt05FsVNcaMSBAg8FBXaQuh6hZGEaI2mBW3TTf/acLcbawVOpyyM5xOi9ZK4YRojYI7OqMCXXTfDclZElcDdm7PRmFyNdXoquLI6YM9JO6HKJmYxghaqMZ0UEAgC8TL3E3X5LUxsOmVpF7h3SHSsnpvGQ9GEaI2mhcH2/4a9S4Wl6DnSfypS6H7NSlq+XYWzeQ+sHhQRJXQ9QyDCNEbaSQyzB9mOk//w2H2VVD0th8JBtCACPDPBHi5SJ1OUQtwjBC1A6mDwuEQi5DwvkryCwokbocsjO1BiM2H8kG8Fu3IZE1YRghage+GjUmhJsGsm48nC1xNWRv4k4VoKCkCp4ujriln6/U5RC1WKvCyKpVqxAcHAy1Wo3o6GgkJCQ0ee7s2bMhk8muefTv37/VRRNZohl1/fRfJXEgK3WujXXdg/cO7Q5HJb9jkvVp8W/t5s2bsWjRIixduhRJSUmIjIzEpEmTUFBQ0Oj5b7/9NvLy8syP7OxsdO3aFffdd1+biyeyJGN6d0OAhxN0FTX4IS1P6nLITmRfKcf+M3UDV4exi4asU4vDyIoVKzB37lzMmTMH/fr1w+rVq+Hs7Iy1a9c2er5Go4Gvr6/5cfToUVy9ehVz5sxpc/FElkQhl+GBYaYVWT/nmiPUSTYdyYIQwOieXgjmwFWyUi0KI9XV1UhMTERsbOxvLyCXIzY2FvHx8c16jY8++gixsbHo0YM7SZLtub9uIOuRC1dxWsuBrNSxagxGfHHUtGs0B66SNWtRGCkqKoLBYICPj0+D4z4+PsjPv/H6Crm5ufjxxx/x6KOPXve8qqoq6PX6Bg8ia+DjrkZs3/qBrGwdoY4Vl65FYUkVvFxVuLmfz40vILJQnTrSaf369fDw8MC0adOue96yZcug0WjMj8DAwM4pkKgd1C849U1yDqpqOZCVOk59q8h9Q7vDQcGBq2S9WvTb6+XlBYVCAa1W2+C4VquFr+/1p5MJIbB27Vr8+c9/hqOj43XPXbx4MXQ6nfmRnc2pkmQ9burVDb7uahSX1yAuvfGB3URtVaCvxN4M0+/XfUO6S1wNUdu0KIw4OjpiyJAhiIuLMx8zGo2Ii4tDTEzMda/dt28fMjMz8cgjj9zwfVQqFdzd3Rs8iKyFQi7D3YMDAABbjjJIU8f4+lgOjAIY0qMLQru5Sl0OUZu0uF1v0aJFWLNmDdavX4/09HQ8/vjjKCsrM8+OWbx4MWbOnHnNdR999BGio6MRERHR9qqJLNy9dd9U950uRIG+UuJqyNYIIcxBl60iZAuULb1g+vTpKCwsxJIlS5Cfn4+oqCjs2LHDPKg1Ly8PWVkNB+7pdDp89dVXePvtt9unaiILF9rNFUN7dMHRi1fx9bEczB8bJnVJZEOOZRfjbGEZ1A5yTBnoJ3U5RG0mE0IIqYu4Eb1eD41GA51Oxy4bshqbErLw4tdpCOvmgt2LxkImk0ldEtmIv21Nw8bDWbhrUADemh4ldTlETWru5zeHXxN1kCkD/aB2kONsYRmSs4ulLodsRGWNAd+n5AJgFw3ZDoYRog7ipnbA5AhTE/qWxEsSV0O2YueJfJRU1iLAwwkjQj2lLoeoXTCMEHWg+oGs36fkcvM8ahdb6tYWuWdId8jl7Poj28AwQtSBRoR6onsXJ5RU1mLniRuvUkx0PTnFFTh4tggAu2jItjCMEHUguVyGewabPjS+ZFcNtdHXiZcgBDAitCsCuzpLXQ5Ru2EYIepg9V01BzKLkFNcIXE1ZK2EEPgyqW759yHcIoNsC8MIUQcL7OqMEaFdIQTwbXKO1OWQlUrKuoqLl8vh7KjAbQOuv/0GkbVhGCHqBHcPMrWObE3KgRUs7UMWaOsxU5C9NcIXzo4tXq+SyKIxjBB1glsH+MJRKceZglKczNNLXQ5ZmepaI7an5gEApkUFSFwNUftjGCHqBO5qB8T29QYAfHOMXTXUMvtPF+JqeQ26uakwMoxri5DtYRgh6iT132i/S8mFwciuGmq+rXVjje6I9IdSwf+2yfbwt5qok4zr4w0PZwdo9VX49dxlqcshK6GvrMHuk1oAwF2D2EVDtolhhKiTOCrlmDzAtDz8VnbVUDPtOJ6Pqlojenq7or8/Nwol28QwQtSJ6r/Z7jiej4pqLg9PN1Y/HfyuQQHc+ZlsFsMIUScaEtQF3bs4obSqFrvTtVKXQxYuX1eJQ2dNXXp3RPpLXA1Rx2EYIepEcrkMd0aZPlS4ABrdyHcpORACGBbchcu/k01jGCHqZPWzavZmFOJKWbXE1ZAl23osFwAwjQNXycYxjBB1sl4+bogIcEetUWB7aq7U5ZCFysgvQXqeHg4KGabUDXwmslUMI0QSqG8d4awaaso3dd144/t4w8PZUeJqiDoWwwiRBO6I9IdMBiRlFSP7SrnU5ZCFEULg+xRTq9mdXP6d7ADDCJEEvN3ViA7pCgDYnpYncTVkaZKzi3HpagWcHRWYEO4tdTlEHY5hhEgiU+umatZ/Ayaq932KKaDe3M8HTo4Kiash6ngMI0QSuS3CDwq5DCdy9ThbWCp1OWQhDEaBbXUDm6cO5NoiZB8YRogk0tXFEaN7egEAtqWwq4ZMjly4goKSKrirlbipt5fU5RB1CoYRIgmZu2pScyEEd/IlmFtFbo3whUrJLhqyDwwjRBK6pb8PHBVyZBaUIkNbInU5JLFagxE/pOUD+C2oEtkDhhEiCbmrHTCuTzcAHMhKwKGzl3GlrBqeLo6ICfWUuhyiTsMwQiSx282zavLYVWPn6gPpbQN8oVTwv2eyH/xtJ5JYbF9vODkokHWlHKmXdFKXQxKpqjVgx4m6LhrOoiE7wzBCJDFnRyUm9jUtbMWuGvu1/3QRSipr4eOuwrDgrlKXQ9SpGEaILED9YMVtqXkwGtlVY4/qg+jtA/0hl8skroaoczGMEFmAsb27wU2lRL6+EkcvXpW6HOpkFdUG7E7XAuAsGrJPDCNEFkDtoMDN/X0AANtT2VVjb/ZkFKC82oDuXZwQ2V0jdTlEnY5hhMhCTBngBwD48Xg+u2rsTP1miVMG+EEmYxcN2R+GESILMbqXF9xUShSUVLGrxo5UVBvwc3oBAGByXSAlsjcMI0QWQqVU4OZ+pq6aH9K4V4292JtRgIoaUxfNQHbRkJ1iGCGyIJPNXTWcVWMv6rtoJrOLhuxYq8LIqlWrEBwcDLVajejoaCQkJFz3/KqqKrz00kvo0aMHVCoVgoODsXbt2lYVTGTLbupt6qrR6quQlMWuGltXWWPAz6fYRUPU4jCyefNmLFq0CEuXLkVSUhIiIyMxadIkFBQUNHnN/fffj7i4OHz00UfIyMjA559/jj59+rSpcCJbpFIqEFvXVbOdXTU2b2/dLJoAD86iIfvW4jCyYsUKzJ07F3PmzEG/fv2wevVqODs7N9nSsWPHDuzbtw8//PADYmNjERwcjJiYGIwaNarNxRPZInNXTRpn1di6+h16Jw/wZRcN2bUWhZHq6mokJiYiNjb2txeQyxEbG4v4+PhGr/nuu+8wdOhQLF++HAEBAejduzf++te/oqKiosn3qaqqgl6vb/Agshc39fKCa90CaMey2VVjqyprDIirW+iMXTRk71oURoqKimAwGODj49PguI+PD/Lz8xu95ty5czhw4ACOHz+OrVu3YuXKlfjyyy/xxBNPNPk+y5Ytg0ajMT8CAwNbUiaRVVM7KBBbt1fN9tTG/12R9dt3uhBl1Qb4a9SICvSQuhwiSXX4bBqj0QiZTIYNGzZg+PDhmDx5MlasWIH169c32TqyePFi6HQ68yM7O7ujyySyKJxVY/vqp2/fxlk0RC0LI15eXlAoFNBqtQ2Oa7Va+Pr6NnqNn58fAgICoNH8Njirb9++EELg0qVLjV6jUqng7u7e4EFkT8b07gZXlRJ5ukokXyqWuhxqZ6YuGs6iIarXojDi6OiIIUOGIC4uznzMaDQiLi4OMTExjV4zatQo5ObmorS01Hzs9OnTkMvl6N69eyvLJrJtagcFJtZ11fyQylk1tmb/6UKUVtXCT6PGIHbRELW8m2bRokVYs2YN1q9fj/T0dDz++OMoKyvDnDlzAJi6WGbOnGk+f8aMGfD09MScOXNw8uRJ7N+/H8899xwefvhhODk5td9PQmRjJv9urxoh2FVjS8xdNBF+kMvZRUOkbOkF06dPR2FhIZYsWYL8/HxERUVhx44d5kGteXl5yMrKMp/v6uqKXbt24amnnsLQoUPh6emJ+++/H//617/a76cgskFje3eDs6MCOcUVSL2kQyS/QduEqtrfd9E03r1NZG9kwgq+cun1emg0Guh0Oo4fIbuyYGMStqfmYf7YMLx4W7jU5VA72JNRgDnrjsDbTYVfF09kywjZtOZ+fnNvGiILdluE6ZvzjuN57KqxETvqFjqb1N+XQYSoDsMIkQUb38cbjko5LlwuR4a2ROpyqI1qDUb8dNIURuqDJhExjBBZNBeVEmN6dQNgWh6erFvC+Su4Wl6DLs4OGB7SVepyiCwGwwiRhfutq4ZhxNrtOGH6O7y5nw+UCv73S1SP/xqILFxsXx8o5TJkaEtwrrD0xheQRTIahTlQ3hbBhc6Ifo9hhMjCaZwdEBPmCeC3b9ZkfY5lX0VBSRXcVEqM7OkpdTlEFoVhhMgK1H+TZleN9aof8zOhrzdUSoXE1RBZFoYRIitwS38fyGVA6iUdLl0tl7ocaiEhhLlVi7NoiK7FMEJkBbxcVRgWbJp9wdYR63MiV49LVyugdpBjbG9vqcshsjgMI0RW4lbOqrFaPx437UUzrrc3nBzZRUP0RwwjRFaiPowkZl1Fgb5S4mqouYQQ+LF+Fg33oiFqFMMIkZXw0zghKtADQgA7T2qlLoeaKbOgFOcKy+CokGNCOLtoiBrDMEJkRepbR3ayq8Zq1HerjerpCTe1g8TVEFkmhhEiKzKpvymM/HruMorLqyWuhppjZ91eNLdyFg1RkxhGiKxIiJcL+vi4odYoEJdeIHU5dAPZV8pxPEcPucy0ki4RNY5hhMjKTKrvquFqrBbvp7qxPcOCu8LTVSVxNUSWi2GEyMpM6m/6hr3vdCHKq2slroaup35sT333GhE1jmGEyMr083NHYFcnVNUasf90odTlUBMKS6pw5OIVAL+1ZhFR4xhGiKyMTCbDpH5cAM3S7U7XQghgQIAGAR5OUpdDZNEYRoisUP3MjLhTBaiuNUpcDTWmfkwPZ9EQ3RjDCJEVGhzUBV6uKpRU1iL+3GWpy6E/0FfW4GBmEYDfxvgQUdMYRoiskFwuwy11H3KcVWN59pwqQI1BIKybC3p6u0ldDpHFYxghslL1MzR+OqGFwSgkroZ+76cTpim9nEVD1DwMI0RWKibUE25qJYpKq3As66rU5VCdyhoD9mSYFqTjeBGi5mEYIbJSjko5JtZtvMauGstx4EwRyqsN8NeoMSBAI3U5RFaBYYTIitV/895xIh9CsKvGEuyoC4a39PeFTCaTuBoi68AwQmTFxvTuBrWDHNlXKpCeVyJ1OXav1mBEXDrHixC1FMMIkRVzdlTipl7dALCrxhIkXLiCq+U16OLsgGHBXaQuh8hqMIwQWbn6b+AMI9Krn0UT29cHSgX/eyVqLv5rIbJysX29oZDLcCq/BFmXy6Uux24JIfDTCW6MR9QaDCNEVs7D2RHRIV0BsHVESmk5OuTqKuHsqMDoXl5Sl0NkVRhGiGwAu2qkV3/vx/XpBrWDQuJqiKwLwwiRDahfGj4x6yoKSiolrsY+7eSqq0StxjBCZAP8NE6I7K6BEMCuk1qpy7E7mQWlyCwohYNChvF1C9ERUfMxjBDZiFvMXTUMI52tvosmJswL7moHiashsj4MI0Q2or57IP5sEfSVNRJXY19+m0XjI3ElRNapVWFk1apVCA4OhlqtRnR0NBISEpo8d+/evZDJZNc88vM50I6oPfX0dkVYNxfUGAT2nCqQuhy7kaerQMolHWQy4OZ+DCNErdHiMLJ582YsWrQIS5cuRVJSEiIjIzFp0iQUFFz/P7+MjAzk5eWZH97e7Fclam+cVdP56hc6GxzUBd5uaomrIbJOLQ4jK1aswNy5czFnzhz069cPq1evhrOzM9auXXvd67y9veHr62t+yOXsISJqb/VhZG9GISprDBJXYx92souGqM1alAiqq6uRmJiI2NjY315ALkdsbCzi4+Ove21UVBT8/Pxw88034+DBg9c9t6qqCnq9vsGDiG5sYHcN/DRqlFcbcOBMkdTl2LyrZdU4fP4KAE7pJWqLFoWRoqIiGAwG+Pg0/Abg4+PT5BgQPz8/rF69Gl999RW++uorBAYGYty4cUhKSmryfZYtWwaNRmN+BAYGtqRMIrslk8lwS924BXbVdLy4UwUwGAXCfd3Qw9NF6nKIrFaH95X06dMH8+bNw5AhQzBy5EisXbsWI0eOxFtvvdXkNYsXL4ZOpzM/srOzO7pMIpsxKcL0DX13uha1BqPE1di2ndyLhqhdtCiMeHl5QaFQQKttuI6BVquFr2/z/zEOHz4cmZmZTT6vUqng7u7e4EFEzTM8uCu6ODvgankNEi5ckbocm1VeXYv9pwsBMIwQtVWLwoijoyOGDBmCuLg48zGj0Yi4uDjExMQ0+3WSk5Ph5+fXkrcmomZSKuSI7WvqqvmJC6B1mH0ZhaiqNSKoqzP6+rlJXQ6RVWtxN82iRYuwZs0arF+/Hunp6Xj88cdRVlaGOXPmADB1scycOdN8/sqVK/Htt98iMzMTx48fxzPPPIOff/4ZCxYsaL+fgoga+P0UXyGExNXYpt/PopHJZBJXQ2TdlC29YPr06SgsLMSSJUuQn5+PqKgo7NixwzyoNS8vD1lZWebzq6ur8eyzzyInJwfOzs4YOHAgdu/ejfHjx7ffT0FEDYzu5QVnRwXydJVIvaRDZKCH1CXZlOpaI+LqFpa7NYJdNERtJRNW8LVJr9dDo9FAp9Nx/AhRMy3YkITtaXl4YlwYnr81XOpybMq+04WYtTYB3dxUOLx4IuRytowQNaa5n99ceYzIRtXPqtnBKb7tbsdx0z29pZ8PgwhRO2AYIbJR4/t0g6NCjnOFZcgsKJG6HJthMArsOmkaGMxZNETtg2GEyEa5qR0wqqcngN++yVPbJWVdRVFpFdzVSowI9ZS6HCKbwDBCZMN+m1XDKb7tZWddsJvY1weOSv4XStQe+C+JyIbF9vOBXAak5ehw6Wq51OVYPSGEeQwOu2iI2g/DCJEN83JVYWhwVwBcAK09nMzT49LVCqgd5Bjbu5vU5RDZDIYRIht3a3/Oqmkv9V00Y3t3g5OjQuJqiGwHwwiRjbulv2lBwqMXrqCotEriaqxb/dgbdtEQtS+GESIb172LMwYEaGAUME9JpZY7V1iKDG0JlHIZJob7SF0OkU1hGCGyA/VLlv/IKb6tVn/vRvb0gsbZQeJqiGwLwwiRHbitLowcyiyCrrxG4mqsU/1aLbdxLxqidscwQmQHQru5oo+PG2qNArvT2VXTUtlXypGWo4NcZloCnojaF8MIkZ1gV03r7aybiTQ8pCs8XVUSV0NkexhGiOzEbQNMYWT/mUKUVtVKXI11+dHcReMncSVEtolhhMhO9PFxQ4iXC6prjdhzqkDqcqyGVl+JxItXAXBKL1FHYRghshMymczcVcON85qvvotmcJAHfDVqiashsk0MI0R2pH4myJ6MAlTWGCSuxjr8mMYuGqKOxjBCZEcGBGgQ4OGE8moD9p0ulLoci3e5tAqHz18G8NsAYCJqfwwjRHaEXTUts+ukFkYBRAS4I7Crs9TlENkshhEiO1PfVbM7XYvqWqPE1Vg2zqIh6hwMI0R2ZnBQF3i7qVBSWYuDZ4ukLsdi6SpqcKju/rCLhqhjMYwQ2Rm5XGaeorojjV01TYlL16LGINDbxxVh3VylLofIpjGMENmh+gXQdp7MR42BXTWN+SEtDwBwK7toiDocwwiRHYoO8YSXqyOKy2tw6OxlqcuxOLqKGuw/beqiuX0gwwhRR2MYIbJDCvlvs2q2p+ZKXI3l2X1Si2qDEb28XdHbx03qcohsHsMIkZ2aMsAfALDzhJZdNX9Q30Uzha0iRJ2CYYTITg0P6QovVxV0FTU4mMlZNfV0FTXYf8a0INyUAQwjRJ2BYYTITinkMkweUN9VkydxNZZj10nTLJo+Pm7oxS4aok7BMEJkx+q/+e88kc8F0OrUj6FhFw1R52EYIbJjQ4O7opubCvrKWnbVANCV1+CXM6b7MJldNESdhmGEyI4p5DJMrptVs41dNdh5Mh+1RoFwXzf09OZCZ0SdhWGEyM5NGWiaVfPTSXbVmGfRsFWEqFMxjBDZuaE9ftur5kBmodTlSKa4vBoH6rtoOF6EqFMxjBDZOblcZh4fYc9dNT+d0KLWKNDXz5170RB1MoYRIjLPHNl1QouqWoPE1Uhjm7mLhjv0EnU2hhEiwpCgLvB1V6Okqhb7Muyvq+ZKWTUOZXIWDZFUWhVGVq1aheDgYKjVakRHRyMhIaFZ1x08eBBKpRJRUVGteVsi6iByuczcOvJdiv3tVbM9LQ+1RoGIAHeEsouGqNO1OIxs3rwZixYtwtKlS5GUlITIyEhMmjQJBQUF172uuLgYM2fOxMSJE1tdLBF1nGlRAQCA3elalFbVSlxN5/r2WA6A3+4BEXWuFoeRFStWYO7cuZgzZw769euH1atXw9nZGWvXrr3udfPnz8eMGTMQExPT6mKJqONEBLgj1MsFlTVG/HQiX+pyOk32lXIcvXgVMhlwe900ZyLqXC0KI9XV1UhMTERsbOxvLyCXIzY2FvHx8U1et27dOpw7dw5Lly5t1vtUVVVBr9c3eBBRx5LJZLizrmXg22T76aqp75aKCfWEr0YtcTVE9qlFYaSoqAgGgwE+Pj4Njvv4+CA/v/FvUmfOnMGLL76Izz77DEqlslnvs2zZMmg0GvMjMDCwJWUSUSvdEWVqGTiQWYSi0iqJq+kc39UFrzuj2CpCJJUOnU1jMBgwY8YMvPrqq+jdu3ezr1u8eDF0Op35kZ2d3YFVElG9EC8XRHbXwGAUdrGTb3qeHhnaEjgq5Lg1grNoiKTSvKaKOl5eXlAoFNBqtQ2Oa7Va+PpeOze/pKQER48exbFjx/Dkk08CAIxGI4QQUCqV+OmnnzBhwoRrrlOpVFCpVC0pjYjayZ1RAUi5pMO3yTmYNTJY6nI6VH131PjwbtA4OUhcDZH9alHLiKOjI4YMGYK4uDjzMaPRiLi4uEYHprq7uyMtLQ3Jycnmx/z589GnTx8kJycjOjq67T8BEbWr2yP9IJcBSVnFyLpcLnU5HcZoFPgu2TSL5k7OoiGSVItaRgBg0aJFmDVrFoYOHYrhw4dj5cqVKCsrw5w5cwCYulhycnLwySefQC6XIyIiosH13t7eUKvV1xwnIsvg7abGqJ5e+OVMEb5LycGTE3pJXVKHOHrxKnJ1lXBTKTEh3FvqcojsWovDyPTp01FYWIglS5YgPz8fUVFR2LFjh3lQa15eHrKystq9UCLqPHdE+uOXM0X4JjkXC8b3hEwmk7qkdvdtXavIpAhfqB0UEldDZN9kQgghdRE3otfrodFooNPp4O7uLnU5RDavpLIGQ/61G9W1RmxfOBr9/TVSl9SuqmuNGP7v3Sgur8Fnj0RjdC8vqUsisknN/fzm3jREdA03tQNi+5q6Lr6zwTVHfjlTiOLyGni5qhAT5il1OUR2j2GEiBp1R6RpUOc3yTkwGC2+AbVFvq5b/n1qpB8UctvrgiKyNgwjRNSoCeHe6OLsAK2+CvvP2M5OvsXl1dh1wrQ8wT2Du0tcDREBDCNE1ARHpRzTBplaR7YctZ2FB79NzkW1wYh+fu6ICLCtsTBE1ophhIiadN8Q01YMu05qcaWsWuJq2scXdcHq/qFsFSGyFAwjRNSkfv7uiAhwR41BmKfCWrMTuTqcyNXDUSHnQmdEFoRhhIiu6/6hptaRL45ekriStttS9zPc3M8HXVwcJa6GiOoxjBDRdd0R6Q9HpRzpeXocz9FJXU6rVdUa8E1d68597KIhsigMI0R0XR7Ojriln2mFZWseyBqXXoDi8hr4uqtxU69uUpdDRL/DMEJEN1TfVfNNci4qawwSV9M69QNX7xkSwLVFiCwMwwgR3dConl7w16ihq6jB7nSt1OW0WL6uEvtPm9ZKqZ8hRESWg2GEiG5IIZfhniGmcRbWOJD1q6RLMApgeHBXBHu5SF0OEf0BwwgRNcu9dWHklzOFyC2ukLia5hNCmMe6cOAqkWViGCGiZunh6YIRoV0hBLDpiPUMZI0/exkXLpfDxVGByQP8pC6HiBrBMEJEzfanET0AABsPZ6Gq1joGsn586AIA4O7B3eGiUkpbDBE1imGEiJptUn9f+LqrUVRahR/S8qQu54ayr5SbB9zOGtlD4mqIqCkMI0TUbA4KOf40IggA8PGhixJXc2Of/XoRRgHc1MsLPb3dpC6HiJrAMEJELfLA8CA4KuRIyS7GsayrUpfTpIpqg3lsy6yYYGmLIaLrYhghohbxclXh9kjTQND1deMxLNE3yTnQVdQgsKsTxod7S10OEV0HwwgRtdickSEAgO1peSgoqZS4mmsJIcxBaVZMMFdcJbJwDCNE1GIDumswOMgDNQaBjYezpC7nGr+eu4JT+SVwclDgvqFccZXI0jGMEFGrzB5lah3ZcDgL1bVGiatpaL15Om8ANE4O0hZDRDfEMEJErXJbhC+83VQoLKnCj8ctZ5rvpavl+OlkPgBg1shgaYshomZhGCGiVnFQyPFQtGntjnUHL0AIIXFFJp/WTecdGeaJ3j6czktkDRhGiKjVZkQHQaWUIzm7GAcyi6QuB1fKqvFZvGn9kzl13UhEZPkYRoio1bq5qTAj2rQI2opdpyVvHflg/zmUVRvQ398dsX05nZfIWjCMEFGbPD4uDGoHOY5lFWPf6ULJ6igqrcIn8RcAAH+J7Q2ZjNN5iawFwwgRtYm3mxp/qhs78paErSMf7D+H8moDBnbXYCJbRYisCsMIEbXZvLFhcHJQIOWSDnsyCjr9/QtL2CpCZM0YRoiozbq5qTAzpr515Eynt46s3ncWlTVGRAV6YFyfbp363kTUdgwjRNQuHhsTCmdHBdJydNid3nmtIwX6Snz2q2kGzV9uZqsIkTViGCGiduHpqsLsukXGOnPsyP/2nkVVrRFDenTBmF5enfKeRNS+GEaIqN3MvSkUriolTubp8UNafoe/X05xBTYmmPbG4VgRIuvFMEJE7aaLiyMeHm1abOzV709AV17TYe8lhMDfvk5Dda0R0SFdMaqnZ4e9FxF1LIYRImpXT4wLQ6iXCwpKqvDP7Sc77H2+TLyEfacL4aiU47W7BrBVhMiKMYwQUbtSOyiw/N6BkMlMgWFvB0z11eor8c9tpqDzl9je6Ont2u7vQUSdp1VhZNWqVQgODoZarUZ0dDQSEhKaPPfAgQMYNWoUPD094eTkhPDwcLz11lutLpiILN/Q4K7mwayLv05DSWX7ddcIIfDS1jToK2sxsLsGc2/iHjRE1q7FYWTz5s1YtGgRli5diqSkJERGRmLSpEkoKGj824+LiwuefPJJ7N+/H+np6Xj55Zfx8ssv44MPPmhz8URkuZ6b1AdBXZ2Rp6vEsh9PtdvrfpeSi93pBXBQyPDmvZFQKtjAS2TtZKKF8++io6MxbNgwvPvuuwAAo9GIwMBAPPXUU3jxxReb9Rp33303XFxc8OmnnzbrfL1eD41GA51OB3d395aUS0QSij97GQ+u+RUAsPHRaIzs2bapt4UlVbj5rX0oLq/Bopt7Y+HEXu1RJhF1kOZ+frfoK0V1dTUSExMRGxv72wvI5YiNjUV8fHyzXuPYsWM4dOgQxo4d2+Q5VVVV0Ov1DR5EZH1iwjzxpxGmXX2f+zIVOcUVrX6tyhoDnt2SguLyGvT1c8fj48Laq0wikliLwkhRUREMBgN8fHwaHPfx8UF+/vXXFOjevTtUKhWGDh2KBQsW4NFHH23y3GXLlkGj0ZgfgYGBLSmTiCzIi7f1RbCnM3KKK3D/6nhkXS5v8WtUVBsw95Oj2H+6ECqlHG/eOxAO7J4hshmd9q/5l19+wdGjR7F69WqsXLkSn3/+eZPnLl68GDqdzvzIzs7urDKJqJ25qpT4/LERCPFyQU5xBe57/xAyC0qbfX1pVS1mr0vAL2eK4OyowLo5wxARoOnAiomosylbcrKXlxcUCgW0Wm2D41qtFr6+vte9NiTENOJ9wIAB0Gq1+Pvf/44HH3yw0XNVKhVUKlVLSiMiC+anccLmeSPwpw8P47S2FA98EI/PHo1GuO/1x4DpKmowe10CjmUVw02lxLo5wzA0uGsnVU1EnaVFLSOOjo4YMmQI4uLizMeMRiPi4uIQExPT7NcxGo2oqqpqyVsTkZXzdlNj02Mx6OfnjqLSajzwwa/44kg2isurrzm3qtaAuHQtZqz5FceyiqFxcsCGudEMIkQ2qkUtIwCwaNEizJo1C0OHDsXw4cOxcuVKlJWVYc6cOQBMXSw5OTn45JNPAJjWJAkKCkJ4eDgAYP/+/fi///s/LFy4sB1/DCKyBl1dHPH53BGYtS4BydnFeP6rVPxtqwyjenph8gBfdHF2xI7j+dh1UouSqloAgKeLIz57NBp9/TiTjshWtTiMTJ8+HYWFhViyZAny8/MRFRWFHTt2mAe15uXlISsry3y+0WjE4sWLcf78eSiVSoSFheGNN97AvHnz2u+nICKroXF2wIZHo7Hu4HlsS83DqfwS7DtdiH2nCxuc5+Ouwm0RfnhkdAgCuzpLVC0RdYYWrzMiBa4zQmS7zhaW4se0PPx4PB8llbWY2NcbUwb4YXBQF8jl3G+GyJo19/ObYYSIiIg6RIcsekZERETU3hhGiIiISFIMI0RERCQphhEiIiKSFMMIERERSYphhIiIiCTFMEJERESSYhghIiIiSTGMEBERkaQYRoiIiEhSDCNEREQkKYYRIiIikhTDCBEREUmKYYSIiIgkpZS6gOYQQgAwbUVMRERE1qH+c7v+c7wpVhFGSkpKAACBgYESV0JEREQtVVJSAo1G0+TzMnGjuGIBjEYjcnNz4ebmBplM1m6vq9frERgYiOzsbLi7u7fb69K1eK87F+935+G97jy8152nve61EAIlJSXw9/eHXN70yBCraBmRy+Xo3r17h72+u7s7f7E7Ce915+L97jy8152H97rztMe9vl6LSD0OYCUiIiJJMYwQERGRpOw6jKhUKixduhQqlUrqUmwe73Xn4v3uPLzXnYf3uvN09r22igGsREREZLvsumWEiIiIpMcwQkRERJJiGCEiIiJJMYwQERGRpOw6jKxatQrBwcFQq9WIjo5GQkKC1CVZvWXLlmHYsGFwc3ODt7c3pk2bhoyMjAbnCCGwZMkS+Pn5wcnJCbGxsThz5oxEFduG119/HTKZDM8884z5GO9z+8rJycGf/vQneHp6wsnJCQMGDMDRo0fNz/N+tw+DwYBXXnkFISEhcHJyQlhYGP75z3822NuE97p19u/fj6lTp8Lf3x8ymQzffPNNg+ebc18rKyuxYMECeHp6wtXVFffccw+0Wm3bixN2atOmTcLR0VGsXbtWnDhxQsydO1d4eHgIrVYrdWlWbdKkSWLdunXi+PHjIjk5WUyePFkEBQWJ0tJS8zmvv/660Gg04ptvvhEpKSnijjvuECEhIaKiokLCyq1XQkKCCA4OFgMHDhRPP/20+Tjvc/u5cuWK6NGjh5g9e7Y4fPiwOHfunNi5c6fIzMw0n8P73T5ee+014enpKbZt2ybOnz8vtmzZIlxdXcXbb79tPof3unV++OEH8dJLL4mvv/5aABBbt25t8Hxz7uv8+fNFYGCgiIuLE0ePHhUjRowQI0eObHNtdhtGhg8fLhYsWGD+s8FgEP7+/mLZsmUSVmV7CgoKBACxb98+IYQQRqNR+Pr6ijfffNN8TnFxsVCpVOLzzz+XqkyrVVJSInr16iV27dolxo4daw4jvM/t64UXXhCjR49u8nne7/YzZcoU8fDDDzc4dvfdd4uHHnpICMF73V7+GEaac1+Li4uFg4OD2LJli/mc9PR0AUDEx8e3qR677Kaprq5GYmIiYmNjzcfkcjliY2MRHx8vYWW2R6fTAQC6du0KADh//jzy8/Mb3HuNRoPo6Gje+1ZYsGABpkyZ0uB+ArzP7e27777D0KFDcd9998Hb2xuDBg3CmjVrzM/zfrefkSNHIi4uDqdPnwYApKSk4MCBA7jtttsA8F53lObc18TERNTU1DQ4Jzw8HEFBQW2+91axUV57KyoqgsFggI+PT4PjPj4+OHXqlERV2R6j0YhnnnkGo0aNQkREBAAgPz8fABq99/XPUfNs2rQJSUlJOHLkyDXP8T63r3PnzuG9997DokWL8Le//Q1HjhzBwoUL4ejoiFmzZvF+t6MXX3wRer0e4eHhUCgUMBgMeO211/DQQw8B4O92R2nOfc3Pz4ejoyM8PDyaPKe17DKMUOdYsGABjh8/jgMHDkhdis3Jzs7G008/jV27dkGtVktdjs0zGo0YOnQo/v3vfwMABg0ahOPHj2P16tWYNWuWxNXZli+++AIbNmzAxo0b0b9/fyQnJ+OZZ56Bv78/77UNs8tuGi8vLygUimtGAGu1Wvj6+kpUlW158sknsW3bNuzZswfdu3c3H6+/v7z3bZOYmIiCggIMHjwYSqUSSqUS+/btwzvvvAOlUmn+dsP73D78/PzQr1+/Bsf69u2LrKwsAPy9bk/PPfccXnzxRTzwwAMYMGAA/vznP+Mvf/kLli1bBoD3uqM05776+vqiuroaxcXFTZ7TWnYZRhwdHTFkyBDExcWZjxmNRsTFxSEmJkbCyqyfEAJPPvkktm7dip9//hkhISENng8JCYGvr2+De6/X63H48GHe+xaYOHEi0tLSkJycbH4MHToUDz30EJKTkxEaGsr73I5GjRp1zRT106dPo0ePHgD4e92eysvLIZc3/GhSKBQwGo0AeK87SnPu65AhQ+Dg4NDgnIyMDGRlZbX93rdp+KsV27Rpk1CpVOLjjz8WJ0+eFI899pjw8PAQ+fn5Updm1R5//HGh0WjE3r17RV5envlRXl5uPuf1118XHh4e4ttvvxWpqanizjvv5LS8dvD72TRC8D63p4SEBKFUKsVrr70mzpw5IzZs2CCcnZ3FZ599Zj6H97t9zJo1SwQEBJin9n799dfCy8tLPP/88+ZzeK9bp6SkRBw7dkwcO3ZMABArVqwQx44dExcvXhRCNO++zp8/XwQFBYmff/5ZHD16VMTExIiYmJg212a3YUQIIf773/+KoKAg4ejoKIYPHy5+/fVXqUuyegAafaxbt858jtFoFK+88orw8fERKpVKTJw4UWRkZEhXtI34YxjhfW5f33//vYiIiBAqlUqEh4eLDz74oMHzvN/tQ6/Xi6effloEBQUJtVotQkNDxUsvvSSqqqrM5/Bet86ePXsa/f951qxZQojm3deKigrxxBNPiC5dughnZ2dx1113iby8vDbXJhPid8vaEREREXUyuxwzQkRERJaDYYSIiIgkxTBCREREkmIYISIiIkkxjBAREZGkGEaIiIhIUgwjREREJCmGESIiIpIUwwgRERFJimGEiIiIJMUwQkSSKC4uhkwmw8GDBwEAmZmZCA8Px8svvwzuUkFkX5RSF0BE9ik1NRUymQyRkZE4cOAApk+fjmXLlmHmzJlSl0ZEnYxhhIgkkZKSgrCwMHz77bd4/vnn8fnnn2PMmDFSl0VEEmA3DRFJIiUlBfn5+Zg9ezZ8fX0ZRIjsGMMIEUkiJSUFQ4cOxd69e3Hs2DF88803UpdERBKRCY4UI6JOZjAY4Orqis2bN+OOO+7A/fffj9OnT+PYsWOQyWRSl0dEnYwtI0TU6U6fPo3KykpERUUBAF555RWkpqbiq6++krYwIpIEwwgRdbqUlBR4eHggKCgIADBgwADcfffd+Pvf/w6j0ShxdUTU2dhNQ0RERJJiywgRERFJimGEiIiIJMUwQkRERJJiGCEiIiJJMYwQERGRpBhGiIiISFIMI0RERCQphhEiIiKSFMMIERERSYphhIiIiCTFMEJERESSYhghIiIiSf0/vMVEa+8lDM0AAAAASUVORK5CYII=",
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "plt.plot(rhos)\n",
        "plt.title(r\"$\\sqrt{\\rho(GG*)}$ for $\\mu={mu}$\")\n",
        "plt.xlabel(r\"$\\kappa$\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "(1.0010439001657412, 0.28761639607360895)"
            ]
          },
          "execution_count": 10,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "rhos.max(), rhos.min()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "(1.0005218139379777, 0.5362987936529495)"
            ]
          },
          "execution_count": 11,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "np.sqrt(rhos.max()), np.sqrt(rhos.min())"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "start with $dt=1$ and reduce it by `x10` each time\n",
        "\n",
        "1 (, 1.0016857181744783)\n",
        "\n",
        "so apparently it is $1+O(\\Delta t)$ "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "rho_max = np.array([1.177993438292792,\n",
        "\n",
        "1.0028383249771309,\n",
        "                    \n",
        "1.0002834709998514\n",
        "          ])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "dt = [1, 0.1, 0.01]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "Text(0.5, 1.0, '$max_\\\\kappa \\\\sqrt{\\\\rho(GG*)}$')"
            ]
          },
          "execution_count": 14,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhAAAAEiCAYAAACoWuuyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAujUlEQVR4nO3deXRUZZ4+8KeyLySVhOxJVbEJAQIJhCSVVhEU2wVRQBYhCaijto7rMLbbeEan+9j+Rs44HKdz2u4RdKyAgnSDikqriNBKKiSBLMi+3ewbIZV9q3p/f1RSSQxgilRya3k+53COeeveqm8lyH1y61Y9CiGEABEREZEV3OQegIiIiBwPAwQRERFZjQGCiIiIrMYAQURERFZjgCAiIiKrMUAQERGR1RggiIiIyGoMEERERGQ1BggiIiKyGgMEERERWY0BgoiIiKzmIfcARORchBBwc+PvJlfD+iFyFgqWaRGRLR0+fBgHDhzAb3/7W7lHIaJRxF8TiMim9uzZg3vuuUfuMYholDFAEJFNnThxAtOnT5d7DCIaZQwQRGQzlZWViImJkXsMIhoDDBBEZDO2evnirbfeQlxcHEwmkw2mGhvvvvsu1Go1Ojs75R6FaEwwQBCRzRw4cAC33HLLiO6jqakJ//mf/4kXX3zxqu/mOHfuHJ566ilMnToVPj4+8Pf3R0JCAp5//nmcPXvWZvv0+f777/Hggw9ec5sHH3wQXV1d+POf/zys50nk6Pg2TiKyiY6ODnh4eMDT03NE97Nlyxb09PRgzZo1V7x969atePTRR+Hv748HHngA8fHxMJlMOHr0KDZv3ox3330XTU1Ng8LH9exjMBhw4sQJaLXaQY/f2NiIU6dOITU1ddC6j48P1q9fj7fffhtPP/00FArFiL4PRHZPEBHZwJ49e8QHH3ww4vuZPXu2yMjIuOJtn376qXB3dxeLFi0S9fX1Q25vaGgQr7322oj3EUKIoqIiMXXqVPHMM8+Izz//XKxfv17s3LlTTJw4Ufz3f//3FefLz88XAMS+fft+8XkSOToGCCKZLVq0SGi1WnHo0CFxyy23CD8/PzF58mTxxRdfCCGE+OKLL0Rqaqrw8/MTCQkJIj8/f9D+x48fF7/5zW/E1KlTha+vr4iIiBBr164VlZWVlm3Ky8uFt7e3eOihhwbt+8033wgPDw/x3HPPjfh5PPnkk6K2tvYXt/vuu+/E3LlzhZ+fn/Dz8xOrVq0SjY2NQgghzp8/LwBcMYg0NjaK8PBwMX36dNHS0jKsma5nn4E6OzvFxo0bRXR0tPD39xerVq0SFy9evOY+ISEh4plnnrH6sYgcDa+BIJJZcXExDAYDVq9ejYULF+IPf/gDmpubsXbtWvz5z3/GM888g2XLluHVV1/FuXPn8PDDDw/af8uWLbhw4QIyMzPxzjvvYNWqVdi5cycyMzMt28TExOCRRx5BdnY2JEkCAJw8eRIrV67EXXfdhf/6r/8a8fOoq6tDWFjYNbd57733sGjRIkyfPh0bN27E4sWLsWPHDjz33HMAgEOHDgEA5s6dO2TfTZs2oba2Fm+99Rb8/f2HNdP17DOQQqGAm5ub5eUIhULxiy9NzJ07Fz/++KPVj0XkcOROMESurKamRgAQERERg84YvPPOOwKAiIuLEwaDwbK+YcMGoVAoREdHh2WttbV1yP2++uqrws3NTbS3t1vW+s5CPPHEE6K+vl5MnjxZJCYmWvWbeXFxsdi0adOQ9cLCQvHGG29cc98TJ04IDw8P8c477wxanz9/vvDx8RHd3d3i1VdfFQBEc3PzkP3j4uJEdHS0MBqNg9YNBoOoq6uz/Bn4fK5nn4HPNS4uTjz99NOWlzA++eQTMXHixCt+D/o89thjwtfX95rfCyJnwDMQRDIqLi4GALz++uuIioqyrI8bNw4AsHHjRgQGBlrWlUol3NzcBl3s5+fnZ/nvxsZG1NfXIzg4GCaTCT09PZbbYmJi8Oijj2LLli1YvHgx2tvbsWfPnmH/Zn7gwAGsW7cOGzduHHLbcN6++frrr2P27Nl46qmnBq3Pnz8fHR0daGhowKVLl+Dh4WF5/n1qa2tx8uRJzJ8/f8g7M1JSUhAWFmb586c//em69xlIrVbj/fffxzvvvGOZZ8WKFThy5MiQCysHCg4ORnt7O9ra2q75/SBydAwQRDIqKSkBANx7772D1k+dOgVfX1/cfvvtg9ZPnz6NyZMnW97pYDQasWXLFsydOxd+fn4IDg5GWFgY/vVf/xWhoaFDDsTPP/88Ojs7UVxcjM8++8yqD3265ZZbsH//ftTV1eHIkSNDnsfs2bOvum9PTw++/PJLrFixYshLAK2trVAoFIOC0s+Vl5cDACZOnDjktj/96U/45ptv8NJLLwGAZY7r2WcgpVJ5xaAQFBQ05B0YA4neeiG+C4OcHd/GSSSj4uJiREVFITo6etB6UVER4uPj4e3tPWR94MFu3bp12LFjB9avX4/nnnsOoaGh8PLywjPPPDPojEafN954A4D5gB4SEmL1vEFBQViwYAE+++wzy3UKw7n24ciRI2hubkZiYuKQ2woLC5GQkAAfHx+MHz8ePT09aG5uRkBAgGWbrq4uy9w/t3DhQgDmMyQAkJCQcN37XM2CBQuwYMGCa27T5/Lly/Dz84Ovr++wtidyVDwDQSSj4uLiKx68ioqKhqx3d3fj1KlTlgBx/PhxbNu2DRs3bsR7772HdevW4e6778bMmTNx+vTpIfv3bffHP/4RHh4eljAx0ObNm3HfffcBMF9kOWvWLOzdu3fQNkuWLMHnn39u+frLL7/E3Xfffc3nWVhYCABDXi6pqqrCDz/8gGXLlgEA4uLiAAAXLlwYtN2kSZMAAD/99NNVH6OkpAQRERGIiIi47n1s4cKFC+wCIZfAAEEkE6PRiOPHjw850NfX16OqqmrI+okTJ9Dd3T3kFP2UKVMs23R3d+ORRx6B0WgcdKZi9+7deOmll/D73/8eTz75JB577DF8+OGHQw7UfS9FHDp0CPfffz8++OAD3HnnnYO2WbJkCY4ePYqKigoAwP79+y2/0V9N37Uefb/xA+YzA0888QSUSiV+85vfAADS0tIAAPn5+YP2Dw8Px/z58/HVV1/hhx9+GHL/3d3dQ87OXM8+tnDkyBH86le/sul9EtkjvoRBJJMzZ86go6NjSFAoKioCMPS0et9BuO+Al5iYCD8/Pzz77LO4cOEC2tvbodPp4O7uPmj/goICpKenIz09Hf/2b/8GAHjhhRfw7rvv4o033sB7771neYySkhJERkbivvvuQ25uruW3+IE0Gg3i4+OxZ88ePPzwwzCZTPDx8bnmcy0pKcGMGTPwxhtvoLW1FVFRUfj444+Rl5eHnTt3DjprEB8fj2+//XbI21WzsrJw8803Y+HChVi3bh3mzJkDIQTOnDmDTz75BDU1NXjyySdHvM9IFBQUoKGhwXIWh8ipyfwuECKXtWPHDgFAHDt2bND622+/LQBYPlypzwsvvCACAwOFyWSyrO3Zs0fExcUJHx8fMXv2bPH++++L119/XXh4eIiOjg5RVlYmoqKixI033jjorZ9CCPHEE08IT09Pcf78ectaWFiYuPnmm4VSqRQHDhy46uyvvPKKWLx4sfjmm2/E//7v//7icw0ODhavvvqq+Mtf/iJUKpXw9vYWaWlpV/zExrfffluMGzdOtLW1Dbnt3Llz4sEHHxTR0dHCw8NDKJVKMW/ePPHCCy+IkpKSKz729exzvV588UWhVqsH/YyInJVCiN5LhonIpdXU1GDChAlobGzE7373Oxw+fBjffPPNFbfNzc3FggULsG7duiFvQf25srIyqNVqbNu27ar9FgMZDAZMmjQJb731Fv7pn/7pup/PWOvs7MSECRPw0ksv4dlnn5V7HKJRx2sgiAiA+WWG6dOnw9vbG8899xwOHToEvV5/xW1TUlKgVCpRVFR0zfDQd78AMHPmzGHNoVQq8cILL2Djxo0OVef9/vvvw9PTE48//rjcoxCNCQYIIgJgPtDHx8cDAMLCwvDII4/g97///RW3VSgUWLx4Me66665fvN/i4mK4u7tj2rRpw57lxRdfxMmTJ69a522PHn/8cZSWlg556y2Rs+JLGER0XfLy8hAUFIQbbrjhmtulp6cjPz8fp06dGqPJiGgsMEAQERGR1Rzn/CARERHZDQYIIiIishoDBBEREVnN6T6J0mQyobKyEgEBAWzDIyIisoIQAs3NzYiOjv7Fd0E5XYCorKyESqWSewwiIiKHVVZWhtjY2Gtu43QBoq8CuKysDIGBgTJPQ0RE5DiampqgUqksx9JrcboA0feyRWBgIAMEERHRdRjOJQC8iJKIiIis5jQBIisrCzNmzEBycrLcoxARETk9p/skyqamJiiVShgMBr6EQURETqvK0I4L9a2YGOqPKKWvTe7TmmOo010DQURE5Oy255Xi5b+VwCQANwXw5vJZWJ2sHtMZnOYlDCIiIldQZWjHS73hAQBMAnjlb8dQZWgf0zl4BoKIiMgBdPYYsfdYNbL2n8XPLz4wCoGL9W02eyljOBggiIiI7Fj55TZsyy3F9rwyXGrtuuI27goFJoT6jelcDBBERER2xmQSOHCmDlv1EvadrLWccYgM9MHaVDV8Pd3x/746CaMQcFco8Ifl8WN69gFggCAiIrIbl1u7sCO/DFtzS1Ha0GZZv2lKKDK0GiyaHg4Pd/Pli/ckROFifRsmhPqNeXgAGCCIiIhkJYRAYVkjdHoJe4qr0NVjAgAE+HhgZZIK6Vo1JoeNG7JflNJXluDQhwGCiIhIBu1dRnxWVAGdXsKxiibL+szoQKxL02BJQjT8vOz3MG2/kxERETmh83UtyNaXYmdBGZo6egAAXh5uuGd2FDK1GiSqgobVRSE3BggiIqJR1mM04dsTtcjWS/jhbL1lXRXii4xUDVbOUyHE30vGCa3HAEFERDRKaps68HFeGbbllqK6qQMAoFAAt04LR0aaBrfcEAY3N/s/23AlDBBEREQ2JIRA7oUG6PQS/n6sGj29HxkZ4u+F1ckqrE1RQxUytp/ZMBoYIIiIiGyguaMbu45WQJcj4Uxti2U9SROMTK0Gd82KhLeHu4wT2hYDBBER0QicrG6CLkfCrqMVaOsyAgB8Pd2xdE4MMrRqzIxWyjzh6GCAICIislJfL0W2XkLexcuW9clh/sjUarA8KRaBPp4yTjj6GCCIiIiGqaKxHdtyJWzPK0N9i7mXwt1NgTtmRiBDq0HapPEO8RZMW2CAICIiugaTSeAfZ+uhy5Hw3ckaS412RKA31qSosSZFjYhAH3mHlAEDBBER0RVcbu3CzoJyZOdKkC7191LcOGU8MrUa3DY9Ap69vRSuiAGCiIhogKLeXorPiyrR2ddL4e2B+5NikaHVYEr40F4KV8QAQURELq+9y4jPiyuRrZdQXG6wrM+IMvdS3Jto370UcuB3g4iIXNaF+lZs1Uv4pKAchvZuAICXu7mXIiNNgzkO0kshBwYIIiJyKT1GE747WQudXsI/zvT3UsQG+yJDq8HKpFiMH+ct44SOgQGCiIhcQm1zB3b09lJUGvp7KRZOC0emVoP5U8Pg7qC9FHJggCAiIqclhEDexcvQ6SXsPVaFbqP5PZjBfp5YnaxGeqpz9FLIgQGCiIicTnNHN3YfrUC2vhSnapot63PVQchM0+Cu+Cj4eDpPL4UcGCCIiMhpnKpuhk5/EbuOVKB1UC9FNNJTNYiPcc5eCjkwQBARkUPr6jFh70/VyM6RcPhig2V9Ul8vxdxYKH2du5dCDgwQRETkkCoa2/FRbik+zitDfUsnAHMvxa9nRCBTq0HaZNfppZADAwQRETkMk0ngh7P10Okl7DvR30sRHtDfSxGpdL1eCjk4TYDIyspCVlYWjEaj3KMQEZGNNbb19lLoJVwc0EuRNmk8MtM0uH2Ga/dSyEEhhBByD2FLTU1NUCqVMBgMCAwMlHscIiIageLyRuhyJHx2xV4KNaaEB8g8oXOx5hjqNGcgiIjIOXR0G/F5kbmXomhAL0VcZADWpU3AfYnR8Pfm4Utu/AkQEZFduFjfiq25EnbkD+6luHtWJDLTNJirDuZFkXaEAYKIiGRjNAlLL8XB03WW9ZggX6Rr1Vg1T4VQ9lLYJQYIIiIac3XNndiRb+6lqGhsB2DupbhlahgytRosmBbOXgo7xwBBRERjQgiBfOkydDkSvvpZL8WqZBXSUzRQj2cvhaNggCAiolHV0tnT20sh4WR1fy9FoioImVoNFs9mL4UjYoAgIqJRcbqmGdl6CX87UoGWzh4AgI+nG5YmxiBDy14KR8cAQURENtPVY8LXx6uhy5GQe2FAL0WoPzK0Gtw/NxZKP/ZSOAMGCCIiGrEqg7mX4qO8MtQ19/dS3D49AplpGvyKvRROhwGCiIiui8kkcOjcJej0F/HtiVoYe4spwiy9FCpEKX1lnpJGCwMEERFZxdDWjZ1HyrFVL+F8fatlXTspBJnaCfj1TPZSuAIGCCIiGpaScgOy9RI+LapAR7e5l2KctwfunxuDdK0GUyPYS+FKGCCIiOiqOrqN2FNcBZ1eQlFZo2U9LjIAmWkaLE2MYS+Fi+JPnYiIhpAutWJbbim255ehsc3cS+HprsDds6KQqdUgScNeClfHAEFERADMvRT7+3opztRBmK+JREyQL9amqrE6mb0U1I8BgojIxdW3dGJ73uBeCqC/l2JhHHspaCgGCCIiFySEQIF0GTq9hC9L+nspgvw8sWqeCmtT1JgQ6i/zlGTPGCCIiFxIa2cPdhdWQJczuJciobeX4h72UtAwMUAQEbmAM729FH8d0Evh7eGG+xKjkaHVYHZskLwDksNhgCAiclLdRhO+/qkGOv1F6M/391JMDPVHeqoaK5JiEeTnJeOE5MgYIIiInEyVoR0fHS7Dx4dLUdvbS+GmABb19lLcODkUbrwokkaIAYKIyAkI0dtLkSPhmxM1ll6K0HHeWJOiwpoUNaKD2EtBtsMAQUTkwAzt3fhrQTmycyWcr+vvpUiZGIJMrQZ3zIyElwd7Kcj2GCCIiBzQsQpzL8Xuwv5eCn8vdyyfG4sMrQbTItlLQaOLAYKIyEF0dBvxZYm5l+JoaaNlfVpEADLSNFg2Jwbj2EtBY4R/04iI7FzppTZsPSxhR14ZLg/opbgrPgqZaRrMYy8FyYABgojIDhlNAgdO10KXI+H70/29FNFKH6RrNVg1T4WwAPZSkHwYIIiI7Millk7syC/H1lwJ5Zf7eynm9/VSTAuDhzsviiT5MUAQEclMCIEjpY3I1kv4orgKXUbzRZFKX0+smheLtakaTGQvBdkZBggiIpm0dfXg08JK6HIkHK9qsqwnxCqRodVgSUI0eynIbjFAEBGNsbO1LeZeioJyNA/opbg3wdxLkaAKkndAomFggCAiGgPdRhO+PV4DnV7CoXOXLOsTxvshQ6thLwU5HAYIIqJRVG3owEeHS/FxXilqmvp7KW6bHoFMrQY3TWEvBTkmBggiIhsTQiDn3CXo9BK+Pj6wl8ILDySrsSZVjRj2UpCDY4AgIrIRQ3s3/nakHNl6CecG9lJMCEFGmgZ3speCnAgDBBHRCP1U2dtLcbQS7d1GAOZeimVzY5Ch1SAuMlDmCYlsjwGCiOg6dHQb8dWxKuhyJBwZ0EsxNWIcMrUaLJ0TgwAfT/kGJBplDBBERFYoa2jD1txS7MgvQ0NrFwDAw02BO+MjkanVIGViCHspyCUwQBAR/QKjSeDg6Tro9BL2n6q19FJEKX2wNkWN1SkqhAf4yDsk0RhjgCAiuoqG1i7syC/D1lwJZQ39vRQ33xCKDK0Gt8WFs5eCXBYDBBHRAEIIHC1rRHaOhD0lVejqMfdSBPp4YOU8FdJT1ZgUNk7mKYnkxwBBRARzL8VnhZXQ6SX8VNnfSxEfE4h12glYkhANXy/2UhD1YYAgIpd2rs7cS7GzoBzNHeZeCi8PNyyZHY3MNA0SYpW8KJLoChggiMjl9BhN+PaEuZfix7P9vRTqED9kaNVYmaRCsD97KYiuhQGCiFxGbVMHPjpcho8Ol6K6qQOAuZfi1rgIZGjVmH9DGHspiIbJ7gJEY2MjFi1ahJ6eHvT09ODZZ5/Fo48+KvdYROSghBDQn29Atl7C33+qRk9vL8V4fy88kKLCmhQ1YoP9ZJ6SyPHYXYAICAjAwYMH4efnh9bWVsTHx2P58uUYP3683KMRkQNp6ujGriMV0OklnK1tsawnTwhGhlaDO+Mj4e3BiyKJrpfdBQh3d3f4+Zl/G+js7IQQAqLvU1uIiH7B8comZOdK2H20Am1d5l4KPy93LJtj7qWYHsVeCiJbsPoTUA4ePIglS5YgOjoaCoUCu3fvHrJNVlYWJkyYAB8fH6SmpuLw4cNWPUZjYyMSEhIQGxuL3/72twgNDbV2TCJyIZ09RnxaWIEVfzqEu9/5B7bllqKty4gbwsfhd/fNRO4rt+GNZbMYHohsyOozEK2trUhISMDDDz+M5cuXD7l9+/bt2LBhA959912kpqZi06ZNuOOOO3Dq1CmEh4cDABITE9HT0zNk36+//hrR0dEICgpCUVERampqsHz5cqxYsQIRERHX8fSIyJmVX27DttxSbM8rw6UBvRR39PZSpLKXgmjUKMQIXh9QKBTYtWsXli5dallLTU1FcnIy/vjHPwIATCYTVCoVnn76abz00ktWP8Y///M/49Zbb8WKFSuueHtnZyc6OzstXzc1NUGlUsFgMCAwkL9tEDkbk0ngwJk6bNVL2Heyv5ciMtAHa1PVeCBZhfBA9lIQXY+mpiYolcphHUNteg1EV1cXCgoK8PLLL1vW3NzcsGjRIuTk5AzrPmpqauDn54eAgAAYDAYcPHgQTzzxxFW3f/PNN/Ef//EfI56diOzbZUsvRSlKG9os6zdNMfdSLJrOXgqisWTTAFFfXw+j0Tjk5YaIiAicPHlyWPchSRIee+wxy8WTTz/9NGbNmnXV7V9++WVs2LDB8nXfGQgicnxCCBSWNUKnl7CnuL+XIsDHAyuTVEjXqjGZvRREsrC7d2GkpKSgsLBw2Nt7e3vD29t79AYiojHX3mXEZ0Xmt2Aeq+jvpZgZHYh1aRosSYiGn5fd/fNF5FJs+n9gaGgo3N3dUVNTM2i9pqYGkZGRtnwoInJC5+takK0vxc6CMjQN6KW4Z3YUMrUaJKqCeFEkkZ2waYDw8vJCUlIS9u3bZ7mw0mQyYd++fXjqqads+VBE5CTMvRS1yNZL+OFsvWVdFeKLjFQNVs5TIYS9FER2x+oA0dLSgrNnz1q+vnDhAgoLCxESEgK1Wo0NGzZg/fr1mDdvHlJSUrBp0ya0trbioYcesungROTYaps68HFeGbbl9vdSKBTArdPCkZGmwS3spSCya1YHiPz8fCxcuNDydd8FjOvXr8cHH3yA1atXo66uDv/+7/+O6upqJCYmYu/evfwcByKCEAK5Fxqg00v4+7H+XooQfy+sTlZhbYoaqhD2UhA5ghF9DoQ9suY9rEQ0Npo7urHraAV0ORLODOilSNIEI1OrwV2z2EtBZA9k+xwIOWVlZSErKwtGo1HuUYio18nqJuhyJOwa0Evh6+mOpXNikKFVY2a0UuYJieh68QwEEdlUZ48Re49VI1svIe/iZcv65DB/ZGo1WJ4Ui0AfTxknJKKrcckzEEQkr4rGdmzLlbA9rwz1LeZeCnc3Be6YGYEMrQZpk8bzLZhEToQBgoium8kk8I+z9dDlSPjuZA16r4lERKA31qSosSZFjQj2UhA5JQYIIrLa5dYu7CwoR3auBOlSfy/FjVPGI1OrwW3TI+DJXgoip8YAQUTDVtTbS/F5USU6+3opvD1wf1IsMrQaTAlnLwWRq2CAIKJrau8y4vPiSmTrJRSXGyzrM6LMvRT3JrKXgsgV8f96IrqiC/Wt2KqX8ElBOQzt3QAAL3dzL0VGmgZz2EtB5NIYIIjIosdowncna6HTS/jHmf5eithgX2RoNViZFIvx49h+S0ROFCD4QVJE16+2uQM7enspKg39vRQLp4UjU6vB/KlhcGcvBRENwA+SInJRQgjkXbwMnV7C3mNV6Daa/ykI9vPE6mQ10lPZS0HkavhBUkR0Vc0d3dh9tALZ+lKcqmm2rM9VByEzTYO74qPg48leCiK6NgYIIhdxqroZOv1F7DpSgdZBvRTRSE/VID6GvRRENHwMEEROrKvHhL0/VSM7R8Lhiw2W9Ul9vRRzY6H0ZS8FEVmPAYLICVU0tuOj3FJ8nFeG+pZOAOZeil/PiECmVoO0yeylIKKRYYAgchImk8APZ+uh00vYd6K/lyI8oL+XIlLJXgoisg0GCCIH19jW20uhl3BxQC9F2qTxyEzT4PYZ7KUgIttjgCByUMXljdDlSPjsir0UakwJD5B5QiJyZgwQRA6ko9uIz4vMvRRFA3op4iIDsC5tAu5LjIa/N/+3JqLR5zT/0vCTKMmZXaxvxdZcCTvyB/dS3D0rEplpGsxVB/OiSCIaU/wkSiI7ZTQJSy/FwdN1lvWYIF+ka9VYNU+FUPZSEJEN8ZMoiRxYXXMnduSbeykqGtsBmHspbpkahkytBgumhbOXgohkxwBBZAeEEMiXLkOXI+Grn/VSrEpWIT1FA/V49lIQkf1ggCCSUUtnT28vhYST1f29FImqIGRqNVg8m70URGSfGCCIZHC6phnZegl/O1KBls4eAICPpxuWJsYgQ8teCiKyfwwQRGOkq8eEr49XQ5cjIffCgF6KUH9kaDW4f24slH7spSAix8AAQTTKqgzmXoqP8spQ19zfS3H79AhkpmnwK/ZSEJEDYoAgGgUmk8Chc5eg01/EtydqYewtpgiz9FKoEKX0lXlKIqLrxwBBZEOGtm7sPFKOrXoJ5+tbLevaSSHI1E7Ar2eyl4KInAMDBJENlJQbkK2X8GlRBTq6zb0U47w9cP/cGKRrNZgawV4KInIuThMg+FHWNNY6uo3YU1wFnV5CUVmjZT0uMgCZaRosTYxhLwUROS1+lDWRlaRLrdiWW4rt+WVobDP3Uni6K3D3rChkajVI0rCXgogcEz/KmsjGjCaB/X29FGfq0Be7Y4J8sTZVjdXJ7KUgItfCAEF0DfUtndieN7iXAujvpVgYx14KInJNDBBEPyOEQIF0GTq9hC9L+nspgvw8sWqeCmtT1JgQ6i/zlERE8mKAIOrV2tmD3YUV0OUM7qVI6O2luIe9FEREFgwQ5PLO9PZS/HVAL4W3hxvuS4xGhlaD2bFB8g5IRGSHGCDIJXUbTfj6pxro9BehP9/fSzEx1B/pqWqsSIpFkJ+XjBMSEdk3BghyKVWGdnx0uAwfHy5FbW8vhZsCWNTbS3Hj5FC48aJIIqJfxABBTk+I3l6KHAnfnKix9FKEjvPGmhQV1qSoER3EXgoiImswQJDTMrR3468F5cjOlXC+rr+XImViCDK1GtwxMxJeHuylICK6HgwQ5HSOVZh7KXYX9vdS+Hu5Y/ncWGRoNZgWyV4KIqKRYoAgp9DRbcSXJeZeiqOljZb1aREByEjTYNmcGIxjLwURkc3wX1RyaKWX2rD1sIQdeWW4PKCX4q74KGSmaTCPvRRERKPCaQIE2zhdh9EkcOB0LXQ5Er4/3d9LEa30QbpWg1XzVAgLYC8FEdFoYhsnOYxLLZ3YkV+OrbkSyi/391LM7+ulmBYGD3deFElEdL3YxklOQwiBI6WNyNZL+KK4Cl1G80WRSl9PrJoXi7WpGkxkLwUR0ZhjgCC71NbVg08LK6HLkXC8qsmynhCrRIZWgyUJ0eylICKSEQME2ZWztS3mXoqCcjQP6KW4N8HcS5GgCpJ3QCIiAsAAQXag22jCt8droNNLOHTukmV9wng/ZGg17KUgIrJDDBAkm2pDBz46XIqP80pR09TfS3Hb9AhkajW4aQp7KYiI7BUDBI0pIQRyzl2CTi/h6+MDeym88ECyGmtS1YhhLwURkd1jgKAxYWjvxt+OlCNbL+HcwF6KCSHISNPgTvZSEBE5FAYIGlU/Vfb2UhytRHu3+UO+/L3csWxuDDK0GsRF8rM6iIgcEQME2VxHtxFfHauCLkfCkQG9FFMjxiFTq8HSOTEI8PGUb0AiIhoxBgiymbKGNmzNLcWO/DI0tHYBADzcFLgzPhKZWg1SJoawl4KIyEkwQNCIGE0CB0/XQaeXsP9UraWXIkrpg7UpaqxOUSE8wEfeIYmIyOYYIOi6NLR2YUd+GbbmSihr6O+luPmGUGRoNbgtLpy9FERETowBgoZNCIGjZY3IzpGwp6QKXT3mXopAHw+snKdCeqoak8LGyTwlERGNBQYI+kVtXT34rLASOr2Enyr7eyniYwKxTjsBSxKi4evFXgoiIlfCAEFXda7O3Euxs6AczR3mXgovDzcsmR2NzDQNEmKVvCiSiMhFOU2AyMrKQlZWFoxGo9yjOLQeownfnjD3Uvx4tr+XQh3ihwytGiuTVAj2Zy8FEZGrUwjRd928c2hqaoJSqYTBYEBgID+kaLhqmjrw8eEyfHS4FNVNHQDMvRS3xkUgQ6vG/BvC2EtBROTkrDmGOs0ZCLKeEAL68w3I1kv4+0/V6OntpRjv74UHUlRYk6JGbLCfzFMSEZE9YoBwQU0d3dh1pAI6vYSztS2W9eQJwcjQanBnfCS8PXhRJBERXR0DhAs5XtmE7FwJu49WoK3LfK2In5c7ls0x91JMj+JLPkRENDwMEE6us8eIvceqocuRkC9dtqzfED4OmWkaLGMvBRERXQcGCCdVfrkN23JLsT2vDJcG9FLc0dtLkcpeCiIiGgEGCCdiMgkcOFOHrXoJ+07291JEBvpgbaoaDySrEB7IXgoiIho5BggncNnSS1GK0oY2y/pNU8y9FIums5eCiIhsiwHCQQkhUFjWCJ1ewp7i/l6KAB8PrExSIV2rxmT2UhAR0ShhgHAw7V1GfFZkfgvmsYr+XoqZ0YFYl6bBkoRo+Hnxx0pERKOLRxoHcb6uBdn6UuwsKEPTgF6Ke2ZHIVOrQaIqiBdFEhHRmGGAsGPmXopaZOsl/HC23rKuCvFFRqoGK+epEMJeCiIikgEDhB2qberAx3ll2Jbb30uhUAC3TgtHRpoGt7CXgoiIZMYAYSeEEMi90ACdXsLfj/X3UoT4e2F1sgprU9RQhbCXgoiI7AMDhMyaO7qx62gFdDkSzgzopUjSBCNTq8Fds9hLQURE9ocBQiYnq5ugy5Gwa0Avha+nO5bOiUGGVo2Z0UqZJyQiIro6Bogx1NdLka2XkHexv5dicpg/MrUaLE+KRSB7KYiIyAEwQIyBisZ2bMuVsD2vDPUt5l4KdzcF7pgZgQytBmmTxvMtmERE5FAYIEaJySTwj7P10OVI+O5kDXqviUREoDfWpKixJkWNCPZSEBGRg2KAsLHLrV3YWVCO7FwJ0qX+Xoobp4xHplaD26ZHwJO9FERE5OAYIGykqLeX4vOiSnT29VJ4e+D+pFhkaDWYEs5eCiIich5OEyCysrKQlZUFo9E4Zo/Z3mXE58WVyNZLKC43WNZnRJl7Ke5NZC8FERE5J4UQQsg9hC01NTVBqVTCYDAgMDDQJvdZZWjHhfpWTAz1R5TSFxfqW7FVL+GTgnIY2rsBAF7u5l6KjDQN5rCXgoiIHJA1x1D+evwLtueV4uW/lcAkAAWAKeHjBn3gU2ywLzK0GqxMisX4cd7yDUpERDSGGCCuocrQbgkPACAAS3i4NS4cmVoN5k8Ngzt7KYiIyMUwQFzDhfpWS3gY6J0HEnFvYszYD0RERGQn+H7Ca5gY6o+fn1xwVyiQPDFEnoGIiIjsBAPENUQpffHm8llw770g0l2hwB+WxyNK6SvzZERERPLiSxi/YHWyGvOnhuFifRsmhPoxPBAREYEBYliilL4MDkRERAPwJQwiIiKyGgMEERERWY0BgoiIiKzmdNdA9H0yd1NTk8yTEBEROZa+Y+dwWi6cLkA0NzcDAFQqlcyTEBEROabm5mYolcprbuN0ZVomkwmVlZUICAhASkoK8vLyRnyfTU1NUKlUKCsrs1lBFzmO5ORkm/w9cgXO+L1yhOdkDzOO9Qxj8Xij9Ri2ut/RODYJIdDc3Izo6Gi4uV37KgenOwPh5uaG2NhYAIC7u7tND/iBgYEMEC7I1n+PnJkzfq8c4TnZw4xjPcNYPN5oPYa9H5t+6cxDH6e+iPLJJ5+UewRyAvx7NHzO+L1yhOdkDzOO9Qxj8Xij9Rj28POyBad7CWM0WNOPTkRENBbkPjY59RkIW/H29sZrr70Gb29vuUchIiICIP+xiWcgiIiIyGo8A0FERERWY4AgIiIiqzFAEBERkdUYIIiIiMhqDBBERERkNQYIGysrK8OCBQswY8YMzJ49G5988oncIxERkYtbtmwZgoODsWLFCpvdJ9/GaWNVVVWoqalBYmIiqqurkZSUhNOnT8Pf31/u0YiIyEV9//33aG5uxv/93/9h586dNrlPnoGwsaioKCQmJgIAIiMjERoaioaGBnmHIiIil7ZgwQIEBATY9D5dLkAcPHgQS5YsQXR0NBQKBXbv3j1km6ysLEyYMAE+Pj5ITU3F4cOHr+uxCgoKYDQaWS1ORERXNZbHJVtyuQDR2tqKhIQEZGVlXfH27du3Y8OGDXjttddw5MgRJCQk4I477kBtba1lm8TERMTHxw/5U1lZadmmoaEB69atw1/+8pdRf05EROS4xuq4ZHPChQEQu3btGrSWkpIinnzyScvXRqNRREdHizfffHPY99vR0SFuvvlm8eGHH9pqVCIicgGjdVwSQoj9+/eL+++/3xZjCiGEcLkzENfS1dWFgoICLFq0yLLm5uaGRYsWIScnZ1j3IYTAgw8+iFtvvRWZmZmjNSoREbkAWxyXRgsDxAD19fUwGo2IiIgYtB4REYHq6uph3cePP/6I7du3Y/fu3UhMTERiYiJKSkpGY1wiInJytjguAcCiRYuwcuVKfPnll4iNjbVJ+PAY8T3QIDfddBNMJpPcYxAREVl8++23Nr9PnoEYIDQ0FO7u7qipqRm0XlNTg8jISJmmIiIiV2XPxyUGiAG8vLyQlJSEffv2WdZMJhP27duHtLQ0GScjIiJXZM/HJZd7CaOlpQVnz561fH3hwgUUFhYiJCQEarUaGzZswPr16zFv3jykpKRg06ZNaG1txUMPPSTj1ERE5Kwc9rhks/dzOIj9+/cLAEP+rF+/3rLN//zP/wi1Wi28vLxESkqK0Ov18g1MREROzVGPS+zCICIiIqvxGggiIiKyGgMEERERWY0BgoiIiKzGAEFERERWY4AgIiIiqzFAEBERkdUYIIiIiMhqDBBERERkNQYIIiIishoDBBEREVmNAYKIRlVOTg4UCgUWL158xdv/5V/+BcuXLx/jqYhopBggiGhUbd68GWvWrMG+fftQWVk55PbDhw9j3rx5MkxGRCPBMi0iGjUtLS2IiorCvn378Nprr+Hmm2/GK6+8AgDo6uqCv78/enp6LNunpqZCr9fLNS4RWYFnIIho1OzYsQORkZFISUlBeno6tmzZgr7fWTw8PPDjjz8CAAoLC1FVVYW9e/fKOS4RWYEBgohGzebNm5Geng4AWLp0KaqqqnDgwAEAgJubGyorKzF+/HgkJCQgMjISQUFBMk5LRNZggCCiUXHq1CkcOnTIEiDGjRuH++67D5s3b7Zsc/ToUSQkJMg1IhGNAAMEEY2KzZs3Izk5GTfccINlLT09HX/9619hMBgAmF+6YIAgckwMEERkcz09Pfjwww+xdu3aQeu//vWv4efnh48++ggAUFJSgsTERBkmJKKR8pB7ACJyPnv27EFNTQ3i4+Nx7NixQbfNnz8fmzdvxuOPPw6TyYRTp06hsrIS/v7+UCqVMk1MRNbi2ziJyOaWLFmCPXv2XHOboqIiFBcX48UXX0RlZSWef/55bNy4cYwmJKKRYoAgIiIiq/EaCCIiIrIaAwQRERFZjQGCiIiIrMYAQURERFZjgCAiIiKrMUAQERGR1RggiIiIyGoMEERERGQ1BggiIiKyGgMEERERWY0BgoiIiKz2/wGyz5BNdf6+LAAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 600x250 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "plt.figure(figsize=(6,2.5))\n",
        "plt.loglog(dt[1:],rho_max[1:]-1,\".-\")\n",
        "plt.xlabel(r\"$\\Delta t$\")\n",
        "plt.gca().set_xscale(\"log\")\n",
        "plt.gca().set_yscale(\"log\")\n",
        "plt.title(r\"$max_\\kappa \\sqrt{\\rho(GG*)}$\")\n",
        "# plt.savefig(\"asdf\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([6.4, 4.8])"
            ]
          },
          "execution_count": 15,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 0 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "plt.gcf().get_size_inches()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python (Pyodide)",
      "language": "python",
      "name": "python"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "python",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.8"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}