examples/spherical_spheroidal_mixing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ab4aab22",
   "metadata": {},
   "source": [
    "# Example usage of `postmerger.spherical_spheroidal_mixing`\n",
    "\n",
    "Spherical spheroidal mixing coefficients $\\mu_{ml'ln}$ are defined by\n",
    "$$\n",
    "h_{lmn}=\\sum_{l'}\\mu_{mll'n}^*~^{^S}h_{l'mn}~_{-2}{Y}{_{lm}}\n",
    "$$\n",
    "where $^{S}h_{l'mn}$ are the strain components defined in the _spheroidal-harmonics_ basis, while $h_{lmn}$ are the strain components defined in the _spherical-harmonics_ basis. The symbol $^*$ denotes complex conjugation.\n",
    "\n",
    "The above definition follows from Eq. (5) of https://arxiv.org/abs/1408.1860 . However, we apply a correction factor $(-1)^{l+l'}$ to compensate for a different definition of spin-weighted spherical harmonics in the SXS catalog: this is obtained by setting `enforce_sxs_convention=True` (default)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8d1039bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('..')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "79e0e757",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a2749ba2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import postmerger as pm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5efc2fbe",
   "metadata": {},
   "source": [
    "## Read documentation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e549fe79",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function spherical_spheroidal_mixing in module postmerger._utils:\n",
      "\n",
      "spherical_spheroidal_mixing(lm, mode, spin, method='BK14', prograde=1, qnm_method='interp', enforce_sxs_convention=True, s=-2, _fitting_coeffs={})\n",
      "    Returns spherical-spheroidal mixing coefficients mu_{m,l,lt,n} as defined by Eq. (5) in https://arxiv.org/abs/1408.1860 .\n",
      "    Note that they are the complex conjugates of the mixing coefficients C_{lt,l,m,n} defined in https://arxiv.org/abs/1908.10377 .\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    lm : tuple\n",
      "        Indices (lt,m) of the spherical harmonic.\n",
      "    \n",
      "    mode : tuple\n",
      "        Indices (l,m,n) of the spheroidal harmonic.\n",
      "    \n",
      "    spin : float or array_like\n",
      "        Dimensionless spin of the Kerr black hole.\n",
      "    \n",
      "    method : str. Default='BK14'.\n",
      "        The method used to compute mixing coefficients. Allowed options: ['BK14','PT73'].\n",
      "        \n",
      "        If 'BK14', it uses the fitting coefficients presented in https://arxiv.org/abs/1408.1860 and provided at https://pages.jh.edu/eberti2/ringdown/ .\n",
      "        \n",
      "        If 'PT73', it uses the leading order expressions in perturbation theory, see Press & Teukolsky 1973, ` Perturbations of a rotating black hole. II. Dynamical stability of the Kerr metric`.\n",
      "    \n",
      "    prograde : int. Default=1.\n",
      "        Allowed options: [-1,1]. If 1, return mixing coefficients for prograde modes. If -1, return mixing coefficients for retrograde modes.\n",
      "    \n",
      "    qnm_method : str. Default='interp'.\n",
      "        The method used to approximate the Kerr spectrum. Allowed options: ['interp','L18'].\n",
      "        \n",
      "        If 'interp', it interpolates linearly from the numerical tables provided at https://pages.jh.edu/eberti2/ringdown/ . They are only defined for spin in [-0.998,0.998] and any use outside this range is not guaranteed to produce sensible results. Note that we only support 2<=l<=5, but original tables are also available for l=6 and 7.\n",
      "        \n",
      "        If 'L18', it uses the fits in https://arxiv.org/abs/1810.03550 . They are defined for spin in the whole physical range [-1,1].\n",
      "    \n",
      "    enforce_sxs_convention : bool. Default=True.\n",
      "        If True, applies a sign correction to match the conventions of SXS. See footnote 4 of https://arxiv.org/abs/1902.02731 .\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "        mu_re : float or array_like\n",
      "            Real part of the mixing coefficients.\n",
      "        mu_im : float or array_like\n",
      "            Imaginary part of the mixing coefficients.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(pm.spherical_spheroidal_mixing)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5d74bed",
   "metadata": {},
   "source": [
    "## Evaluate mixing coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "66ab4647",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu_{2320} = (0.0665939069543019+0.011046238081249502j)\n"
     ]
    }
   ],
   "source": [
    "## evaluate mu_{2320}\n",
    "\n",
    "## spherical-harmonic indices\n",
    "lm = (3,2)\n",
    "\n",
    "## spheroidal-harmonic indices\n",
    "mode = (2,2,0)\n",
    "\n",
    "## final spin\n",
    "spin = 0.68\n",
    "\n",
    "mu_re, mu_im = pm.spherical_spheroidal_mixing(lm,mode,spin)\n",
    "\n",
    "## results\n",
    "print('mu_{2320} =',mu_re+1j*mu_im)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05031e8e",
   "metadata": {},
   "source": [
    "Note that spherical-spheroidal mixing with $l=l'$ is not null:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "dc506439",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu_{2330} = (0.9899656358674489-0.022394437936684233j)\n"
     ]
    }
   ],
   "source": [
    "## evaluate mu_{2330}\n",
    "\n",
    "## spherical-harmonic indices\n",
    "lm = (3,2)\n",
    "\n",
    "## spheroidal-harmonic indices\n",
    "mode = (3,2,0)\n",
    "\n",
    "## final spin\n",
    "spin = 0.68\n",
    "\n",
    "mu_re, mu_im = pm.spherical_spheroidal_mixing(lm,mode,spin)\n",
    "\n",
    "## results\n",
    "print('mu_{2330} =',mu_re+1j*mu_im)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "182a0d44",
   "metadata": {},
   "source": [
    "Nonetheless, the amount of mixing is small, i.e., the real part is close to 1. \n",
    "\n",
    "The default `method` to compute mixing coefficients is `'BK14'`, which provdes fitting coefficients to numerical tables. Instead, `method='PT73'` resorts to a perturbation theory expansions and can be written in closed form. In this limit, the mixing for $l=l'$ is exactly null:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "63be0cf3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu_{2330} = (1+0j)\n"
     ]
    }
   ],
   "source": [
    "## evaluate mu_{2330} with method='PT73'\n",
    "\n",
    "## spherical-harmonic indices\n",
    "lm = (3,2)\n",
    "\n",
    "## spheroidal-harmonic indices\n",
    "mode = (3,2,0)\n",
    "\n",
    "## final spin\n",
    "spin = 0.68\n",
    "\n",
    "mu_re, mu_im = pm.spherical_spheroidal_mixing(lm,mode,spin,method='PT73')\n",
    "\n",
    "## results\n",
    "print('mu_{2330} =',mu_re+1j*mu_im)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "063b6ec0",
   "metadata": {},
   "source": [
    "while for $l\\neq l'$ it gives results similar to `'BK14'`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2c429898",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu_{2320} = (0.06691801252539002+0.010410156367386333j)\n"
     ]
    }
   ],
   "source": [
    "## evaluate mu_{2320} with method='PT73'\n",
    "\n",
    "## spherical-harmonic indices\n",
    "lm = (3,2)\n",
    "\n",
    "## spheroidal-harmonic indices\n",
    "mode = (2,2,0)\n",
    "\n",
    "## final spin\n",
    "spin = 0.68\n",
    "\n",
    "mu_re, mu_im = pm.spherical_spheroidal_mixing(lm,mode,spin,method='PT73')\n",
    "\n",
    "## results\n",
    "print('mu_{2320} =',mu_re+1j*mu_im)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3ea24a3",
   "metadata": {},
   "source": [
    "## Batch evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d6c2ecec",
   "metadata": {},
   "outputs": [],
   "source": [
    "## evaluate mu_{2320}\n",
    "\n",
    "## spherical-harmonic indices\n",
    "lm = (3,2)\n",
    "\n",
    "## spheroidal-harmonic indices\n",
    "mode = (2,2,0)\n",
    "\n",
    "## final spin\n",
    "spin = np.linspace(-1,1,100)\n",
    "\n",
    "mu_re, mu_im = pm.spherical_spheroidal_mixing(lm,mode,spin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c3750593",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAFICAYAAAAI6BoFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABGMUlEQVR4nO3deZwcZbXw8d/pZZZMksm+JxBIAoQtQNgEZMcECBHcQAVENHAF9UX0CqJX71UQ79WroghGREAFLi6RRMMaNgGBJBBCFkIWCNnInskkma27zvtHVc/UTLp7umd6uno5Xz7FdNfScyrTM6efp546j6gqxhhjTLEIBR2AMcYYkw1LXMYYY4qKJS5jjDFFxRKXMcaYomKJyxhjTFGxxGWMMaaoWOIyxpgyICL3isgWEVmSYruIyB0iskpEFovIsb5tU0RkhbftpvxFnZwlLmOMKQ/3AVPSbJ8KjPeWGcBdACISBu70tk8ELhORiT0aaScscRljTBlQ1ReAHWl2mQ48oK5XgH4iMhw4AVilqmtUtRl42Ns3MJa4jDHGAIwE1vmer/fWpVofmEiQ37yQicg0YFqfPn2+OGHChKDDMcYUgYULF25T1cHdeY2PnFmj23fEs/u+i5uWAo2+VTNVdWaW31qSrNM06wNjiSsFVZ0DzJk8efIXFyxYEHQ4xpgiICJru/sa23fEee2JMVkdEx6+slFVJ3fzW68HRvuejwI2AhUp1gfGugqNMaaAKOBk+V+OzAau8EYXngTUqeomYD4wXkTGikgFcKm3b2CsxWWMMQVFiWvOklErEXkIOAMYJCLrge8CUQBVvRuYC5wPrAL2AVd522Iicj3wBBAG7lXVpTkPMAuWuIwxpoC4La7cX0JS1cs62a7AdSm2zcVNbAXBEpcxxhSYHHb/lSRLXMYYU0AUJW4T/KZlicsYYwpMT3QVlhJLXMaUGUebaI5vJhoaRDjUq3W9agux+DbC4QGEpDLACMubAnFLXGlZ4jKmRKnG2d20CEcb6Ft5LCGpZtPuX7Fp953udhwG11zG6H63sG3PvWyu+ylKDBAG9b6K4f2+iYjdMRMEa3GlZ4krhUTljHHjxgUdijFZ29O8nKWbrybu7AUEJc7QXudT1zgbRxta99u292GaW96jofkl1L9+z+8ISSVDa7/Kzvp72LnnXhytp6bydAb3+xbRyOgk39XkgoJd4+qEJa4U/JUzgo7FmHRUlc1757Bu9wPEnT0MrD6TbXv/RMzZ1W6/HfseQaT9H0RHG9jbNA/p8AlftYGt9b9B4+uob5jdmtTqG+awt+l5xg57nki4W5WNTBo2pjA9S1zGFLlVO3/IhvpHWltSjfXvE6YlyZ6pPsUnX+/oPur3zUJp9q9FnX3s3HMvg2u/6R6tceLOdsKhWsSujXWbonaNqxOWuIwpMvXNq1i58052Nr1JVXgQzS1LvWtTLqUFRferjOrglj3oSKiAdsnJFQ31RWjGncmijdJMQ5Nbv7N+zx/ZuftW1NkHEqJPzeX0r/0OIvanpcsU4pa30rIrr8YUkfrmlby88dN8sG8eTfEt1DcvwfElLSBJynI59PKSVNv2kFQxpO8NiFS321ekmsF9v4ImbbmFqYgczL6Gx9lR9x0cZydKE6oN1O/9PTvrvt/d0yxrbuWM7JZyYx+LjClwLfF63t/zGHta1rK78TXi2kCie08T/2uXqwSHEGFC4CW1kFRTUzGRgwfczMa6O9jXspSqyIGMqP0qfas+RN/qU/hg149oaHmbysgYhtXeSJ/q02lofIyGptfxt8hEKhjQ54ts3XFtuwEdQGvy6l/7Les27DIhnuLDh3FZ4jKmgNU3v8tzGz6Ho83EtZFKibUbYKG+YRXtc1c14wbcws6GF4g7exhccwFDek8nJBVMGHLvft+npvIYDh768H7rRw66n807vsGehsdRlGhkFMP6/5iK6MHE4ylmtlAl7tQRCQ/p+omXMQUc6ypMyxKXMQVK1WHhlu/S4tTT1sLq2BEoNBOhAocwEUQihKSSwwbdxqBeZzG8zye7FUM41IcRg+7GcRpQbSQU6oeIG0FF9Agam17Y7xgJVREODWx3HvHYuyCVRCKjuhVPubAWV3qWuIwpMO/u/jtvbf8VDfEtRKX9TLgxDRMljvj+roWkmqG9pzG+3wzizl56RccikmwYRteFQtVA++tg/Wu/xQdb57frLhSppn/fW1q/f1PTv9i54zpUd6MaJxIZx4CB9xCJHJDT+Ex5scRlTAF5b/djLNz6I+LqzsKuSrsk5RAipkpElIhUocQZUXM+hw/8FiGJ5jXWyoqjGTZ4Fjvrfkhzy2LC4RH063sjNdVTAYjHN7Fj++Wo7ms9JhZbzratFzN02Gs28jAFt+STtbjSsXdOClY5w+RTi7OXLQ1v8eb2n7cmLRLXrzokL5VqRvf+BOP6fYqK8AAioZoAIna5yWv/a2MAe/c+hGq8w1oH1T00Nb1AVdVZPR9gkXLUElc6lrhSsMoZJl9W1P2N+Vt/TkgihHRXu8/accJEiKMKYalAJExtxTgmDryOSKg65WsWAie+HmhKtoV4/IN8h1M0rMXVOUtcxgRoe+MK5m/9OXFtIq5N7l1WHYa2xwhTFe7HEQOuprZiAgOrjm0dIFHIKipPpaFhdruuQnAHa1RUHB9QVIVPEeJ2i21alriMCUDcaWbt3pdYtvOPxH2VKWKEiKrTLnmFpZpJA2/kgL5TA4i066qrL2BP/S+Ixd6jreVVTVX1R4hGxwcYWeGzrsL0LHEZk2e7mzcyZ92/EXMacXQPId99WY5XZTCKe19WTWQ4Rw38EmP6nBdUuF0mUsmgwXPYs+duGhoeRaiipuYKetV8OujQCpp1FXbOEpcxefbcB/9FQ3wnieJMst/IwTAxKvj42L9RFe4XUJS5EQr1pm/fr9O379fT7qfOPpobn0C1jkjFKUTKukUmxNW6CtOxxGVMntQ1b2DprrlsbVxO2w3FgqLtRg6GpYqjB1xV9EkrUy3Nr1O//TMoDqg7kWVlr49RU3t7UVzLyzW3VqElrnQscRmTB8t2Pcbzm3+Oaoyo+K9huRfiBSVKhAN6f5hDai9hWK9jggw3b1Tj1O/4PKr17dY3NcwiWnkGldXFdV0vV6yrMD1LXMb0sIZYHc9v/nnrIIz9780SRCqYUHsJJw25PqAogxFreQNtvW/NR/fRtO+hskxcqtZV2BlLXMb0oA8a3uHVrfej2jb5RItGqJBYa1WMiFTTJzqcYwdeFWCkAdEWhORTWXacB6ycONbiSssSVwpWOcN01z8338OinY8S10bCvvqCitCkEcIoQyrHc9ygyxlT8yFCZVgCKVJxLCT9I92Lql4fz3c4BcEdVZj7FpeITAF+jjuf6D2qenuH7d8APuM9jQCHAYNVdYeIvAfUA3EgpqqTcx5gFqw9moKqzlHVGbW1tUGHYorQ1sbVLNr5KDFtSjEJuyBSzalDv8aBvT9clkkL3CHzvfv/AqgCKryVNUQrjqei+qMBRhYkt6swm6XTV3SrHt8JTAUmApeJyET/Pqr6P6o6SVUnATcDz6vqDt8uZ3rbA01aYC0uY3JuS+O7PPPB3cQ0cdOtOwAjrA4I7hSPEuaYAZ9iaPVhgcZaCCqqzqHf0Bdo2vcX1NlBtPIMopWnIVKen6t7aFThCcAqVV0DICIPA9OBZSn2vwx4KNdB5IolLmNy6Pkt9/Pa9r+BNhAGX/dgiBhChDAH9f4QHxryBWorRgQZakEJh0fQq8+Xgw6jYMRzXzljJLDO93w9cGKyHUWkFzAF8I8UUuBJcWcx/bWqzsx1gNmwxGVMjmxuXMNr22cR02YEYf8ZsQQlwoeHXU9NZEAAEZpi0MVahYNEZIHv+cwOySVZJkw1z/I04KUO3YSnqOpGERkCPCUib6vq/rOI5oklLmNyYHfLdp7ZfB8xbQHcPz4xhIgqCIQlghDi3OE3WNIyPWFbJ9ee1gOjfc9HARtT7HspHboJVXWj93WLiMzC7XoMLHEVbCeyiEwRkRUiskpEbkqyXUTkDm/7YhE51rftXhHZIiJLOhwzQESeEpGV3tf++TgXU9oW7XyGO965ljV73kS17UNsnDBNhHGoYFyf0/nCuD9wWO3ZAUZqioWjoayWDMwHxovIWBGpwE1OszvuJCK1wOnAo751NSLSJ/EYOA9Y0vHYfCrIxJXJCBhv23hvmQHc5dt2H24fbUc3AfNUdTwwz3tuTJftie3i7xt/RUybiWmynhfBIczpQ75Ir0i/fIdXcuLNi9izdTq7N46lftMkGuvvSDJZZXFLDIfPZun0NVVjuNesngCWA4+o6lIRuVZErvXtejHwpKru9a0bCrwoIm8CrwH/UNXHc3W+XVGoXYWZjICZDjyg7kfcV0Skn4gMV9VNqvqCiByY5HWnA2d4j+8HngO+2TOnYEpdzGnhuS3/59YaxLs2oSHCuDcbhyRMWMJ8ZPh19I5a92B3xVtWsHf7p0AbAFDdQXP9nWh8C9X9fhBwdLnjvo9yfwOyqs4F5nZYd3eH5/fhfvD3r1sDHJ3zgLqhUBNXJiNgku0zEtiU5nWHquomAFXd5F1oNCZrHzSs5bdrvktM96Eaa730HSfsJS/lgJpJTB91A73tmlZONO25E7TjjMoNtOz7P6r6fh0J9QsirB5hRXbTK9TElckImGxGyWT3zUVm4HY/MmbMmFy8pCkhqsoD793KvvhuQIlIxzejEJIqzhp6lSWtHHJalgLO/hukAie2lnBFv3yH1CNUsVqFnSjUf51MRsBkM0omYbOIDAfwvm5JtpOqzlTVyao6efDgwVkFbkrfst3z2Rvb7T1zu3VU3T84ABGp4MSBFzK8+qDAYixFocihJP2Tpc1IpJQ+YApOlku5KdQWV+sIGGAD7giYjtOmzgau965/nQjUJboB05gNXAnc7n19NP3uxrRpcZr5/dqf8U79G4RpKwCrhIipEkKpjQzi8rHfYWjVAQFGWpoq+1xPrPFpoMG3tppo9UcJhUpngLBiLa7OFOS/ToYjYOYCa4BVwG+ALyWOF5GHgH8Bh4jIehG52tt0O3CuiKwEzvWeG5ORuZse4p36xcQ0zv4DCIWwVHPq4I9Z0uoh4ehh9Br4AKGIN8BY+lBRczVV/W4NNrAekOtRhaWmUFtcnY6A8UYTXpfi2MtSrN8O2I00JmtxjfPK9qe9G4yFGGEi3jBsEagIVTG0agyTB5wTbKAlLlJ5Ir2HPI6qU7K1DBXB6YFRhaWkYBOXMYVi8a75PLT217T4ugcdQjQj7nQlKnxs1PUcVnsiYdm/0JPJvVJNWgnl2IrKhiUuY9LY2LCOB977JS3aTNgbPeifuThOhINqDuOIfh8KMEpTShQyrYZRtixxGZPGc1vmEtMY4Fbsjoi2zlwcIkw0VMHFIz8fcJSmtAjxMhwpmA1LXMYk0Rhv4I/v/45FO18iJIl7h4SYQkiUMCEO63ssF4/8HP0qBgUaqykt1uLqnCUuY5L41ar/Zc3eVajqft2Djgohonxs1BfoG+0XXJCmZFmLKz1LXCmIyDRg2rhx44IOxeTZpoYNvLt3jddFGCJEHLQteVWEKjll4NmWtEyPUBVrcXXC/nVSUNU5qjqjtrY26FBMHr1T/w6/WHUHzU5iBKHQomHiCI5CZagXl42ZwfSRnw00TlPa4hrKaik31uIyxrN6z2p+8s5PaHaaqGjXUyPENYxIhFMGTeXY/jaC0JggWeIyxvOX9X/xWlpu/beQamv3oCBUhCo5c8i5gcZosqPOLuL7HkabF0JkApGazyLh4UGHlZZCWdYfzIYlLlP2HHV4dssLrKhf0bou5s2rFUZB4Ii+R/OpMZfTN2pdx8VC4xto3jYdnH1AIzT9k+Z99xMd8CChiqOCDi8NKcvuv2xY4jJl785VM1m0azExdQi1u7k4TFyhUiq55uCvEA1FgwzTZCm2+3ZwdtE2FUozaDOxupupGPyPACNLzx0Oby2udCxxmbK2bt96Fu16k2anBSGE4viGvkNFqIJzhpxjSasIOU0vkGz+Lo2tRJ29SKgm/0FlyEo+pWeJy5St+pY9/Gndo7Q4bmUMbb1Hyy39HgmFmTpsKheNuCjIME1XSRVofbINIIX7QcSK7HbOEpcpS9uadnDT4ttojO8D2gZhKCHiqkQkwoXDpvHRkRcGGqfpunCvTxPf82ug0bc2SqjqPEQqggorI461uNKyxGXK0oNrZ7E3thcHh6jQWn/QJYQlwhlDTg0yRNNN4d5fQluWu12GEgEcJDKBSG1hz9+l6tbFNKlZ4krBKmeUrm1NO5i/cxEOSuIG44i41TEA+lfUct24a+hX0S/IME03iUSJDrgbJ7YGbXkbiYwhFD0i6LAyYl2F6VniSkFV5wBzJk+e/MWgYzG58+LWBdy5+vc42uJbK8Q0AigRCfHjo26jMlIZVIgmx0KRgyByUNBhZMy9xmVdhelY4jJlY2+sgTtX/94bQejWO/OPIIxIlOP6H2lJywTOiuymZ4nLlI3nt76Get2BiqBoa/dgWEKMrRnNNQdfHlyAxmD3cWXC2qOm5Kkqv1z5EPes/gstrV2EgkPILZ6LcEz/o/jBkd+kJtIr0FiNwesqzGbJ6FVFpojIChFZJSI3Jdl+hojUicgib/mPTI/NN2txmZL37Jb5PLdlAS3qENX23YMgVIQquGD4WUGFZ8x+cl2rUETCwJ3AucB6YL6IzFbVZR12/aeqXtjFY/PGWlym5D264VmavOK5MSeEqjvkGIWoRDhv6GkcXjsh6DCNAdqGw2ezZOAEYJWqrlHVZuBhYHqGIXXn2B5hLS5TsvbGGvjh8vtYs3d96zpFaHFChESJSoQbJlzNiQMLueCqKUc9MKpwJLDO93w9cGKS/U4WkTeBjcDXVXVpFsfmjSUuU7JuX34/b+5aRVw7jiB0S+pURXoxeUBx3NdjykcXSz4NEpEFvuczVXWm73myF9QOz18HDlDVPSJyPvA3YHyGx+aVJa4U7Abk4rajeTdv7lpJTGO4v3farjpGZSjKDRMuJyzWW25KwjZVnZxm+3pgtO/5KNxWVStV3e17PFdEfiUigzI5Nt/stzYFVZ2jqjNqa23+pWKzsWE7//HWb2n2iue6U5S4owcdhV7hXvzi2Js5bsDEQOM0JhUn8X7NcMnAfGC8iIwVt1DjpcBs/w4iMkzE/WgnIifg5oftmRybb9biMiVlX6yRLy/8OXUte4i0+1gmKBAizIcHT2Z49eCAIjQmvZ64j0tVYyJyPfAEEAbuVdWlInKtt/1u4OPAv4lIDGgALlVVBZIem9MAs2SJy5SUZ7e8QWO8GcUdmRX2VX4PEaJXpIpPjTk30BiN6UxPlHxS1bnA3A7r7vY9/iXwy0yPDZIlLlMyVu/ZxB/fe4ZGpxlwf/lVlXDIQRAOrR3LzYddwaDKfsEGakw6avNxdcYSlykJK+s3cN2CX9HsNBESfPNrCTEnTHW4ks8eMNWSlil4Su5vQC41lrhMSZi5+jGanGbvOlb7+bUiEmZk9SAm9bMRoqY4WIsrPUtcpuit2fMBr+9Y7d1YIsQ0RFgc8JLXmUMm8eUJH0PE/hiYzqmzA919GzQ+6a6omoL0vRkJ9c/P98cSV2cscZmi9s8ty/juWw/SorF2NxjHNQxApUS54ZBPUhmOBhajKR6qzej2T0B8E+DdTtH4d7RlEQyai0h+/mRa4krP7uMyRSvmxLlt2Z9oclpwVFqnLEmoDEW5YMTxlrRM5pqeAWcHrUkL3MfOVmh6Li8hJCpnZLOUG0tcKYjINBGZWVdXF3QoJoXXd66hIdbkPZPW4rmJBDZl+HFcP/6iwOIzRSi2EnTv/uu1wd2WJz1wA3JJscSVglXOKGx/W/caNyx8gBYn3rpOCbV+Ap3Qeww3HvoxIqFwgFGaohMeC5JkTjaphshB+YlBsRZXJyxxmaKzubGOHy//Oy1OHAc6dBEKVaFKLj3gtICiM0Wt6jyQvrgFIhLCILVQmZ852xKDMyxxpWaJyxSdpzYtRr1sFXfct3CiizBEiOmjTuDsoUcHGaIpUiIVyMBHoPJ03OQVhsozkYGPIJK/a6WWuNKzUYWmqDyy9hV+ueJJHBJdhO7kkAKEED5x4Cl85ZALggzRFDkJD0P6342q4z7P8wwCXZzWpKxY4jJFY9HOtfz87cdp0fh+82spEAlFuGBkupkdjMlcvhOWn1riSqtguwpFZIqIrBCRVSJyU5LtIiJ3eNsXi8ixnR0rIt8TkQ0isshbzs/X+Zjue/i9f9HktJBIVP4RhNFQmCsPOp2D+wwNMkRjcsJGFaZXkC0uEQkDdwLn4k5iNl9EZqvqMt9uU3Fn5xyPO430XcCJGRz7U1X9cZ5OxeTI7HVv8MwHy1qnXVUERREgKhG+cdhFfHS0tbZM8VO1G5A7U6gtrhOAVaq6RlWbgYeB6R32mQ48oK5XgH4iMjzDY00ReeaDZdz61hyanfh+Iwjd1CWcPeyIgKIzxuRboSaukcA63/P13rpM9uns2Ou9rsV7RSQ/xcdMt/xqxTM0Oi2t/f7+5FUVinLdhHPpE60KKDpjck9VslrKTaEmrmQ/Cc1wn3TH3gUcDEwCNgE/SfrNRWaIyAIRWbB169aMAjY945Wta1hdv8V7JsQdwVFwFFDhh8dcymfGnhpkiMbkmJV86kyhJq71wGjf81HAxgz3SXmsqm5W1bi641x/g9utuB9Vnamqk1V18uDBNsV7UOauf4vrXn2IWLv+QUE1hOOE6B2p4dTBEwKLz5ieYi2u9Ao1cc0HxovIWBGpAC4FZnfYZzZwhTe68CSgTlU3pTvWuwaWcDGwpKdPxHSNow4/WvIEjfEWHCe0XwHdqnCUrxxyDqEAhywb0xOsckbnCnJUoarGROR64Anc29fvVdWlInKtt/1uYC5wPrAK2Adcle5Y76X/W0Qm4b433gOuydtJmaws2/UBO5oTxU6FuBMiFHIH/oZE+K+jL+YjI44MMkRjeoZ2LGNmOirIxAWgqnNxk5N/3d2+xwpcl+mx3vrLcxym6QHPf7CSr776CI6q74ql4Dhu/bgDeg+ypGVKWjnem5WNgk1cpjzFHIdvLvgrTY43MaT6K2S4XYTXHPLhwOIzpqe5N9db4krHEpcpKG/t3EBj3J3ET72kleg2CYtw/SFncOGoowKM0JieVp7XrbJhicsUjCU7N/GFlx6iyUnMPpuY1djNXEf3H8NV408JKjxj8saucaVnicsUhJjjcPWLD1Lf0kjI+7DpL6JbHY5y2cHHBxWeMXllXYXp2VhiUxBe27aWZieGv5ukbY4t4aLRR3H+yMODDdKYPHDf97m/jyuDwuWf8aoKLRaRl0XkaN+290TkLa84+YIcnm6XWIsrBRGZBkwbN25c0KGUvK2Ne7j9zafY29LsjSJ0q2O4j+CEwQfxvWMuDDBCY/Ir19e4Mixc/i5wuqruFJGpwEzcAuYJZ6rqtpwG1kXW4kpBVeeo6oza2tqgQylpjiqfefYBVtRt6VDTyy2eWx2u4GMHTgokNmOCkuhtyHTJQKfFx1X1ZVXd6T19BbfqUEGyxGUC9drWtWxp3OPVHuwwx5aEOaL/CKaMnBhojMbkWw90FWZSuNzvauAxf0jAkyKyUERmZH1COWZdhSYwMcfhkdWL2Bdr9tZ4Y9+9L5MGjuK+0z5LJGSfr0z5ULpUf3BQh2tPM1V1pu95JoXL3R1FzsRNXP7q1aeo6kYRGQI8JSJvq+oL2QaZK5a4TGC++vIsnt20EkfbjyBEoVc4yrTRR1jSMmWpC6Pht6lquplUMylcjogcBdwDTFXV7a3xqCYKlW8RkVm4XY+BJS77q2ACsWLXFp7btJqmuAO076cPIfSrrOaiA2xySGNypNPC5SIyBvgrcLmqvuNbXyMifRKPgfMIuEC5tbhMIJ7ZsIqY4yUtR0DaMteBfQbw4FlX0CtSEVR4xgRHc38fV4aFy/8DGAj8StwukJjXihsKzPLWRYAHVfXxnAaYJUtcJu/+svotfvbWP4mp0zr8HXWrZFSGwkwfczQDK2uCDtOY4PRA5YwMCpd/AfhCkuPWAEd3XB8kS1wmr+qaGrhl/uM0x+NIiP2K6IZDIT5+kNUiNOXNKmekZ4nL5NU/N71LyLtHSx1FQm3Xt2oiUX53+qUMqe4daIzGBM1qFaZniSsFq5yRexv37ub2N56jId7irRHUAXcAMFw05ggmDx6d+gWMKQM2rUnnbFRhClY5I7dUlSvm/R8b9+5O8mlSqApHueQgmxzSGDdzSXZLmcm4xSUiQ3GHQR4N9AN2AW8CT6nqBz0RnCkdS3ZsZuPe3a01CBO5S4BQSLjq0OM5bnDBVpgxJq+sqzC9TltcInKYiPwZWAZcDkSBD7yvlwNLReTPImJ1eUxSqsqc95bREGubZwtfiaej+o/gG5POCDBCYwqMZrmUmUxaXPcB/wN8RlWbOm70bmabDvwWODmn0ZmS8Oulr/HA8jdQr5yTy31QFY5wwQGHBhabMYWnSyWfykqniUtVT+xkezPwJ28xpp3GeIw73nyJxnisbaX3OxlCGFLdm0vHF9QtIsYErwxbUdmwUYWmR62p20Hcq5DRWj3XUxOt4O8XfI7e0cpggjOmEPVA5YxS06XEJSJVtFUP7g/sAF7CLSPSkLvwTDHbsm8Pn3/qzzQ5cd9aaf00ecSAYfStqAokNmMKmrW40urqcPgHgBbgm8CngJuAZm+9MQD84LVn2dawL+kF5OpwhC8f9aFA4jKm8EmWS3npalfhsA5zvdQBvxGRy3MQkykRT61b5dYjBPd3y5e8bv/QVD40/IBA4jKm4FmLK62uJq4XRORB3Bky64Ba4CPAi7kKLGhWOaN75r2/moZYW4UM9xfR/W0MS4gLD7SRhMakZIkrrS51Farqt4FbgRrgcO/r7ar6rRzGFiirnNF17+zcxpfmzUadjjdSChEJc87ocYRtgkhjkrPKGZ3q6uCMMHAI7o3I96hqzFv/CVW1YfFl7t4lC2lx4u4voG8goQgM7lXDrR86L9D4jDHFrTuDM44BJgEvikiiP+3fchGUKV6qyqKtm4irl7Uc3EWhSqL814nnMqja5toyJh3V7JZy09VrXCNU9TMAInI/cJ+IfDd3YZli9ZMFL7Jyx/bW1lbriCcFR5WjBw8LMDpjikQZJqNsdDVxVYpIpao2qeq7InIh8BBwRA5jM0Vm6769/GbxfGKO0zaK0Mtb0VCYzxw6iSG9bK4tYzpVhtetstHVxPU13BuPPwBQ1XoRmQ5clqvATPFZuHkDkVAIicdJ9BQmPjlO6DeQ/zjpzCDDM6ZoiLW40upS4lLVV5KsiwN/6HZEpig5qsx6Zxl7W9wh8AKtSSsswjGDhyNinyKN6VSZVnzPRldHFQ7FLfm0FHgd+BawF/gfVd2cu/BMsXho2WKeX/tu+4m2PNFQmM9OnBREWMYUofIc4p6Nro4q/COwARgNPAnM877em6O4TJG5d/FCGuNxdwQhtJsr6HunnM1hA4cEGJ0xRcbm40qrq4kroqr3q+ovgTpV/bOqPgmUTMVUEZkmIjPr6uqCDqUobNu313vkGwLvQIWEOXP0QQFGZkwRssSVVlcTV9T3+Brf45KZJsUqZ2Tu7oWvsae52fcL1Fb8c2B1DUN62X1bxmTFEldaXU00F4uIqOtNaJ0J+Zu5C80Ug91NjfzstZeJx9X9GOQbAh8JhfjRGR+xQRnGZCNR8smk1NUW13bc5PVREYlA60zIo3MWmSkKi7dspiIcRhJdhIlPgA4c2n8wHx59YLABGlOERLNbMnpNkSkiskJEVonITUm2i4jc4W1fLCLHZnpsvnWn5NOxWMmnstcci7G3OTEEXhAVxBFCKozt2z/g6IwpUjnuKvTqy94JTAUmApeJyMQOu00FxnvLDOCuLI7Nq1yVfPqdiHwvZ1GZotAUi/Hv857AQdt1EQJURiJcdfSxKY81xuTVCcAqVV0DICIPA9OBZb59pgMPqKoCr4hIPxEZDhyYwbF51dUWV6WIVAKo6rvANODr5LDkU080a0VkgIg8JSIrva/WJOiGJ9esojEWa/vE5/v0d+G4Qzh2+IigQjOmqPVAV+FIYJ3v+XpvXSb7ZHJsXnU1cSVKPgFuySfcDPy1XATVg83am4B5qjoe996zwPtqi9n63XVu4upIYXjvPvkPyJhSkf18XINEZIFvmdHhFZON9uiY8lLtk8mxeZVVV6GIVKhqcx5KPvVUs3Y6cIZ3/P3Ac9hIyC5RVV5Zt96bvqS9mmiUI4cMDSAqY0pA14a4b1PVyWm2r6f94LlRwMYM96nI4Ni8yrbFtcfrlvu9iNwoIueIyDgR+V2O4+qpZu1QVd0E4H1NWs5BRGYkPrls3bq1yydRyp55910WbNiQ9JdsZJ++nHmg3XRsTAGZD4wXkbHerUuXArM77DMbuMK7DHMSbnGJTRkemzXvcs5bIjLMe36EiNSJyGGdHZtt4hoCfAVYiHs96zfACuDULF+nM4E2a1V1pqpOVtXJgwcPzubQsvHX5UtpiMXaD4FXiBLiayecQiTU1V5oY0yuRxV6s9RfDzwBLAceUdWlInKtiFzr7TYXWAOswv3b/qV0x+bgLCcBe4Bp3m1V9wO3qeryzg7MqqtQVXfhdq89l1gnIj8ANmXzOhnoqWbtZhEZrqqbvG7FLTmNukxs37ePV9a5jVpB2v3yVFVEqYqWTAEVYwLRE9OaqOpc3OTkX3e377EC12V6bA7iaRCRJ4AjgW8DMeDHmRybi4/FPyD3gxx6qlk7G7jSe3wl8GiO4y4LMx59lN2NTUk/6akqJ42y+9CN6ZbyKfm0BLgQuBH4nDdWolPZDs64E1jkLW+paiMwHNiXzet0RlVjIpJomoaBexPNWm/73bjZ/3zcZu0+4Kp0x3ovfTvwiIhcDbwPfCKXcZeDd3fu5O2tW4k72m6iSIBQSLjj/AuojFiLy5huKe5klI1lwFjgm5l0ESZk+xdmA3AW8FXgIBHZgHvd6wkRuRi3/3NlplkznZ5o1qrqduDs7sZWzl59fx0t8XjrRJH+369DBwzkzLE2KMOY7simjFMJuML7mtUAv2yvcd2WeCwiUdz7pI70li96XwdTQtObmDYrtm7lB/OedYfAe0NgEiNhKsJhzj7o4MBiM6aklEGRXRE5AXcAyHrgMCDjIdydXuMSka8kqmT4qWqLqr6pqn8Avgs8oaqjgWEZR26Kyh0v/4umeDxpv3ptZSVXHmslnozJiRK/xuXllN8Bt+AO9jsym+MzGZwxDFglIr8WkU+LyHEiMsH7epmI3A28g9vSSow8NCVGVVmwfqN7v4FX/d0/BP4XF1zIgOrqYIM0pkT0RHX4AvOfwGbgl8Bi4AJvMF1GOk1cqvot4BhgJXA18BjuSJC5wOdx7+M6RlW/nXXopmj8+LkX2bmvofXTnQDiuEsY4dAhdr+bMTlTwi0uETkedwLiq7yxCn8GxgHPZ/oaGV3jUtVtuOPrMxpjXwpEZBowbdy4cZ3uW+o27a7nvgVv4CQmi/SJhkJccsTh9KncrzfZGNMVxduKyoiqzqd9rdt3gQnZvIaVN0hBVeeo6oza2tqgQwncK++vc1tYHSeLVBhT24/vnnNWsAEaU2pKuMWVC9nex/VR3DucI8DbuPdzvQksStQANKVld2Mj//vsS+6gDNqSl6KERfjwgQdaeSdjcq0Mk1E2sr2PayZupYwluEPhjwIuBg4Heuc2NFMI7njhFbbv3dv2i9Q6DF6oCEf41KSsBgMZYzJQyl2FuZBt4moA7vRuMH4msVJESv+mgzI1d/kKYo5vRIYvgX3nnDMYN2hgUKEZY8pUtn08P8GrGOznjQwxJealNWvZttet5pWolNE2BF6YemhW11ONMZmya1xpZdviegp4TETOBf6Oe41rsVez0JSQvc3NXP/IHHcwhjsyo7VKRkiESSNH0KfKRhIak3MlPqowF7Jtcc0CFuAmrPOAB4F6EVmW7iBTfF5Y9Z77IPEL5Pt0VxWJ8D8XTQkmMGNM2cu2xTUE+IS/a1BE+uAO0jAlwlHlwfmL2NfS4l7WcrwN3jWuCw6dwMjavgFGaEyJsxZXWtm2uB4GTvevUNV6VX0pdyGZoD204E0Wb/igfZUM3O6LXtEoFxx+SJDhGVP67BpXWtkmroOAP4nIDSJS0lfmRWSaiMysq6sLOpS8e+DVN2iMeTPTdPjlOH3cgZw8dkxQoRlT8hIfEku8VmG3ZJu4/oR7XWs68IqI1IvIyyLyq9yHFqxyrZyxdvsu1u1MnqwjoRDfOu907O4HY3qYtbjSyugal4iMVtV1qvrbDusPAI7GrnGVBMdRrrr/z+7sxtA2jNAztHdvhvSx+8yN6VFl2orKRqaDM9aKyA7c0YSLgDe85W1VXQvM7pHoTF4tfH8Duxoa3SHwIdxPct6AjHBYuHXaudbaMiYfLHGllWniGgBM8pZjgH/HnbEyJiJLgDdU9ZqeCNDkz6zXl7Kv2TeSMJGjFE4cPYqTD7JrW8bkhSWutDKd1mQX7iyVzyXWeZN+HYGbyCblPDKTV/9c+R6PvbWitZXlr5RRHY1wwRGHBhugMWXEugrTy/Y+rlaq2gy87i2myP3uxQXuSMJEPUJfj+Cwvn248ChLXMbkjSWutLqcuEzp2LFnH6+/781K0+EXpiIc4tbp51EVtbeKMXlRpiMFs2ETKRm+/Mc5NLfE3NFMtL+PJBoKc/iIIUGHaExZyfd9XCIyQESeEpGV3tf+SfYZLSLPishyEVkqIl/1bfueiGwQkUXecn73o0rNEleZ27CzjmUbNreVdfLVJgyJ8M0pH6bSWlvG5Ff+7+O6CZinquOBed7zjmLAjap6GHAScJ2ITPRt/6mqTvKWuTmJKgVLXCmUS+WMvy1cRlMs7l7Scmj3yzCiTx8+ebzdomdMvgVQOWM6cL/3+H7gox13UNVNqvq697geWA6MzMl3z5IlrhTKoXLGO5u28dtn57evSaggDlSEQkw5sqSrehlTuPLf4hqqqpvATVC4BdVTEpEDcUeUv+pbfb2ILBaRe5N1NeaSJa4y9qfXFtMSjyf9BehVEeWqU48LKjRjyle2Scv9vR0kIgt8y4yOLysiT4vIkiTL9GzCE5HewF+A/6equ73VdwEH494atQl30uEeYxcvytTuhkb+/vrbON6ADPUNgY+EhK+ddxoDevcKMkRjylJigFSWtqnq5HQ7qOo5Kb+nyGYRGa6qm0RkOLAlxX5R3KT1R1X9q++1N/v2+Q3uRMM9xlpcZerWWc+yt7E5aTdhmBBnHDY20PiMKWv57yqcDVzpPb4SeLTjDuLWe/stsFxV/7fDtuG+pxcDS3ISVQqWuMqQ4yhPvPlOWzFd30hCgKtPP57Bfa2YrjFl5HbgXBFZCZzrPUdERohIYoTgKcDlwFlJhr3/t4i8JSKLgTOBG3oyWOsqLEOzXltCzHHaugh9ohLi+vNODiQuY4wr3yWfVHU7cHaS9RuB873HL5KiF1NVL+/RADuwxFVmduzZxw//+mxbF6GvkG5IhNOti9CY4FnljLQscZWZeYtX4ShJpy6JRsJ886IzggzPGAOWuDphiauM7G5o5JePvUxLPN5+6hLvl+SSyRMZ0b9vgBEaY2wiyc7Z4IwUSrFyxh+ee509DU37jyRUqI5E+OgJhwcanzHGk/9RhUXFElcKpVg549HXltES94oSdnjTX3DcoRw+elhQoRljfAIo+VRUrKuwTDz8z0Vs3lUPtOsdBKAiHOYLZx0fSFzGmCTKMBllwxJXGWiJx/nFnJdw4rhtbGk/pnXcsIGMGtQvmOCMMfspx1ZUNixxlYE3Vm1kb1OL+6TDL0Q4JPz88xflPyhjTHJlet0qG5a4SlxTS4yb7v0H6qh7z5a2v+n48NFDGda/T2DxGWOSsMSVVsENzshkJk5vvykiskJEVonITZ0dLyIHikiDr1TJ3fk6pyA98+YqGppb2n2KSxTxrAiHuHbqSQFGZ4zpyD/a1wZnJFdwiYsMZuIUkTBwJzAVmAhc5puJM93xq30zdF7bkydRCBxHeeDJBTQ2x9w3d4fllEMP4NSJVinDmIJjw+HTKsTENZ1OZuIETgBWqeoaVW0GHvaOy/T4svDoy0tYs2m7e0MjEPKqv4sDNZEIl3zoyKBDNMYkIapZLeWmEBNXJjNxjgTW+Z6vp20K6XTHjxWRN0TkeRE5LfehF5Y/zHud5pjvvi3augmHDejLqYdba8uYgtO1iSTLSiCDM0TkaSDZ3a63ZPoSSdZ19uPbBIxR1e0ichzwNxE53DeDpz++GcAMgDFjxmQYUmF5Zdla3vtgR/tJIr1/oUhI+MEVUwiHCvFzizGmHK9bZSOQxJWDmTjXA6N9z0cBG73HSY9X1SagyXu8UERWAxOABUnimwnMBJg8eXLRvYVUlf+6/6nWWoTS4Y7j/jW9OHR0soasMaYgFN1fnfwqxI/cnc7ECcwHxovIWBGpAC71jkt5vIgM9gZ1ICIHAeOBNT1yBgFbuX4rW3ft2W+CyMTXH1w1hVCoC5ODG2NMASjExNXpTJyqGgOuB54AlgOPqOrSdMcDHwYWi8ibwJ+Ba1V1R57OKW8cR/n3u+bgqLr9qQ5uwvK+jhrQlxMOLc7uT2PKhQ2HT6/gbkDOZCZO7/lcYG6S/VId/xfgLzkNtgC9tnwt23btbb3LWERaW1rRcIgrzp0cYHTGmIyUYTLKRsElLtN1qspv57xKY3PMXdE6MsM1Zkg/LvnwUcEEZ4zJTJm2orJhiauEPP/Gapa9+4H7xk+Ud/K2VVVE+Oy5x9m1LWOKgSWutCxxlZAH5s6nqTnuXrn0hsAn0lSvyijnHX9IgNEZYzKRKPlkUrPEVSIWLHufpas3ub2D3jD4hHBIuPmz51BVGQ0qPGNMNsqwGkY2LHGViP/9w7P737flvferK6OcepRVyTCmWFiLK71CHA5fEERkmojMrKurCzqUTm3cWsfq9duT37el8J9fnEpF1D6jGFMUrORTpyxxpaCqc1R1Rm1tbdChpKWqfO1/ZrUW0k0U0U3ct1VbXcnpxxwccJTGmGwkfo8zXcqNJa4i9/a7m9mwZVfrJ69EEd2QQkSEKy44IdgAjTHZy3OLK4t5EN8Tkbe8OQ0XZHt8rljiKmKOo/zonqdoifk+cvnezAP71vDZ8+2GY2OKTQCVMzqdB9HnTG9OQ/8fl2yO7zZLXEVswZK1vLdhO3jlnfyTRUZDIS77yLF235YxxUZxRxVms3TfdLo3j2F3j8+KJa4i9ru/vuret+V7oye6Cnv3quSjZ1mVDGOKUQAtrkzmQQT3r82TIrLQm/4p2+NzwoaaFamXFq5m6cqNrde1WufcQomEQ1z78VOoqa4INkhjTNdkn4wG+a85ATO96Zla5WAeRIBTVHWjiAwBnhKRt1X1hayj7SZLXEXqt3/6Fy0tTmub2ctZoFBZEea8Uw4NMDpjTFd1sXLGtg7XnPaTg3kQE8XOUdUtIjILOAF4gRTzIPYU6yosQgvfep8Vaza3n7bEN6rwG1edQ68qa20ZU5Syvb6Vm2tcnc6DKCI1ItIn8Rg4D1iS6fG5ZIkrhUK+Afmn98xDnbZrWjhtS3U0ylknTwg2QGNMsel0HkRgKPCiN6fha8A/VPXxdMf3FOsqTEFV5wBzJk+e/MWgY/Fb/d5W3lu/3b2u5V3ccufcckvCf+3zZ1mVDGOKXL5LPmUyD6KqrgGOzub4nmJ/4YqIqnLTbbPapi1xvOTlzV/St1clF5xxRMBRGmO6rQzLOGXDElcRWbRkHVu317e+qUXa7t0SEa685MRA4zPG5IYV2U3PEleRiMcdfvDTfxB3fO9o38MB/XrxqWlWJcOYoqeAY5krHUtcReLl+avZtWuf75pWW96KhIXLLz7RqmQYUyosb6VliasItLTE+cXMeW01CUPuuzoxKKNfnxqmnWdVMowpFdZVmJ4lriIw74XlbN+xt7W11TrDsSoV0TDXfPY0KivsR2lMybAZkNOyv3YFLh53uOe+F2hpjkOYdl2FAL2qKjjrNKuSYUwpsRZXepa4Ctyjf3+D7Tv2IKpoXLxbxt13dTgU4sZ/O9daW8aUkjKd1Tgb9hcvBRGZBkwbN25cYDGoKvf/8SWcuFvLqbVKhqe2TxWnnTQ+qPCMMT3ArVVomSsdK/mUgqrOUdUZtbW1gcXw8COvUre7wVdAt31tsltuvIBw2H6ExpQcJ8ulzFiLq0A5jvL7P7zk3s8h7ghC/4ew4UP6MPnYAwOLzxjTc6zFlZ4lrgL14EMv09jYkrIm4U03XuA+N8aUFrvG1SlLXAWooaGZ39//YluyctqSF8DY0QOZdNSYYIM0xvSQnE1VUrIscRWg++59gVjMaZ0dUkRah8eKwJe/lHI+OGNMCbDh8OlZ4iowO3fs4dG/LvDm2vKua4Xa3sWHHTKCY+3aljGlzVpcadmQtALzm7uecVtb0DpRpDjuEkGY8cUzA43PGGOCZi2uArJu7XaeeWqp+2nL8W42bv3kJUyadABHHW3Xtowpaep+UDWpWeIqIDN/+RSxljiEpfVm47YK8DDjmrMCjM4YkzfWVZiWdRWmICLTRGRmXV1dXr7f8iXrmf+vVe79G07bjcbidRee/KEJjBs/NC+xGGMCplkuZcYSVwr5rpxx10+fIB73pitRIO4lMEeJAP92/dl5icMYEzxJfGjNcCk31lVYAOa/tJIVyza2rfBaWSggwjnnHcHQYf2CCc4Yk39lmIyyYYmrANz548fQuIOEBC9lkfgSDYe49ivnBRecMSa/lLKsP5gNS1wBe/ofi9i0fmfrc/FKOiX6rS/+5An07lMVUHTGmHwTyrP7LxuWuAKkqtz948dbuwY1jjengbu9sjLC566x+7aMKTuWuNIquMEZIjJARJ4SkZXe1/4p9psiIitEZJWI3ORb/wkRWSoijohM7nDMzd7+K0TkIz19Lp35830vUr+7ofVN6l5oBXEUiTt84bqziUbDAUdpjMk7/xRGmSxlpuASF3ATME9VxwPzvOftiEgYuBOYCkwELhORid7mJcAlwAsdjpkIXAocDkwBfuW9TiCam1p44Ffz2t50/q+q9O1XzfRPnRhUeMaYoCSuceVxPq5MGgwicoiILPItu0Xk/3nbviciG3zbzu9+VKkVYuKaDtzvPb4f+GiSfU4AVqnqGlVtBh72jkNVl6vqihSv+7CqNqnqu8Aq73UCcedtc2hpjrfdswXtPkF96etTgwrNGBOwAIbDd9pgUNUVqjpJVScBxwH7gFm+XX6a2K6qc3MRVCqFmLiGquomAO/rkCT7jATW+Z6v99al05VjesTObfU8+bc3vBuMaXfDMaqMHNWfM6ccFURoxphCkP+uwkwaDH5nA6tVdW0uvnm2AhmcISJPA8OSbLol05dIsq6zn17Gx4jIDGAGwJgxua8NeOsND6JOW/teEPeGY9wBhdffPC3n39MYUywCuW7VrsEgIskaDH6XAg91WHe9iFwBLABuVNWd+x+WG4EkLlVNOaGUiGwWkeHeP95wYEuS3dYDo33PRwEbk+zXpWNUdSYwE2Dy5Mk5fQe9s2Q9y15f640elP3eoBMnjeHYkw7O5bc0xhQTpSuJa5CILPA9n+n9HWuVgwZD4nUqgIuAm32r7wK+jxv994GfAJ/P5nWzUYjD4WcDVwK3e18fTbLPfGC8iIwFNuBm/09n8LoPisj/AiOA8cBruQo6E6rKbV/+A+oohAS3DHRb8gqFhOu+Za0tY8pe9gMutqnq5HQ75KDBkDAVeF1VN/teu/WxiPwG+HvmoWevEK9x3Q6cKyIrgXO954jICBGZC6CqMeB64AlgOfCIqi719rtYRNYDJwP/EJEnvGOWAo8Ay4DHgetUNZ7PE3tuzhts2bjTm7ZE2z5Zectp5x7BQYcMz2dIxpgCFMDgjESDAVI3GBIuo0M3oZfsEi7GHd3dYwquxaWq23Ev/HVcvxE43/d8LrDfyBVVnUX7kS7+bbcCt+Ys2Cw0NjTzi2//BY07EA4jaNu0JapEomG+cOOUIEIzxpjbgUdE5GrgfeAT4DYYgHtU9XzveS/cBsU1HY7/bxGZhPsn7b0k23Oq4BJXqXr4l0/TsKfJ7SJ0HK+rMHGpS/jEVacxeHi/QGM0xhSIPA/OyKLBsA8YmGS/y3s0wA4sceXBtg928edfP9O2QtXrw3bfnJW9KviklXYyxoB3A3L5VcPIhiWuPPjFzX8i3hKHUMhNWOINzAAkJHzpP6ZTXVMZaIzGmEJRnmWcsmGJq4eteGMtC55Z2nbHWGuVDPfLgCF9Ofvi4wKJzRhToCxxpWWJqwc5jsPtX/odTtxxW1lxr9XlERG+fecVhEKFOLjTGBMYS1xpWeJKQUSmAdPGjRvX5dd4+v9eYfO67bQW7VB1B2Z4Jh4/lkOPOaB7gRpjSotd4+qUfdRPQVXnqOqM2traLh2/d/c+7rzpYXf4u7+emPdYBL7z6x67sdwYU7QU1MluKTPW4uohv/3Pv9Lc2OINxKB94hJh2uWnUjuoT3ABGmMKl3UVpmWJqwesfXsjjz3wAopAKITE424C85JYVVWEGd+7JOAojTEFyboKO2WJqwfc/oVfu/UIBVBF8R46DiLCV/77MsIRm9nYGJOCtbjSssSVY8//9VXeXbqhbRIVxx1RqF4x3ZHjhnLmxwKbv9IYUwwscaVliSuHGvc18ZNr7nFbVq1D3H0DM0T4j/uuDSw+Y0wxsBuQO2OJK4fu/sYfaGpoAgmhIm6jq7WfUDnj4ycy5pARwQZpjClsSrvbZsz+LHHlyLtL1vHYvc+1zbMVj6Nh7zqWQkVFmBvuuDLdSxhjjMtaXGlZ4soBVeX7n/wJqoo4uHfHafvm/pd+fgWV1VaP0BiTAUtcadkNyCmIyDQRmVlXV9fpvnPveZr1Kz9wRxAmJomE1gkjR40bytQrT+/hiI0xpcH7G5LNUmYscaWQTeWM3//nn7yEFU8cC462Vs34zoNf6elwjTGlQkHVyWopN9ZVmAM7Nu0CvHEY8Zg3OANQZcrnz+LAiaOCDM8YY0qKtbhyYNhBQ1qvaakq6sTReJzK6ijX/cwGZBhjsmRdhWlZ4sqBL/7ociqrK9qqvzsOFRVhvnnfdVRUVQQdnjGm2CQGd2W6lBnrKsyB0y45kUj0Bu791oNsWrOZ4QcP5erbPsNJF9oEkcaYLHWY/sjszxJXjpw8bTInT5scdBjGmFJQhq2obFjiMsaYAqPW4krLEpcxxhSU8rxulQ1LXMYYU0hsPq5O2ajCFLKpnGGMMTmlTnZLmbHElUI2lTOMMSZXFFBHs1q6S0Q+ISJLRcQRkZSjzERkioisEJFVInKTb/0AEXlKRFZ6X/t3O6g0LHEZY0whUQ2ixbUEuAR4IdUOIhIG7gSmAhOBy0Rkorf5JmCeqo4H5nnPe4wlLmOMKTD5bnGp6nJVXdHJbicAq1R1jao2Aw8D071t04H7vcf3Ax/tdlBp2OAMY4wpNIV53WoksM73fD1wovd4qKpuAlDVTSIypCcDscTViYULF24TkbUpNtcCmY7eGARsy+JbZ/PahbAvZHeOhRJzT/0MCyVme4+219M/wwOyiCWpenY+8bT+eVCWh1WJyALf85mqOtO/g4g8DQxLcuwtqvpoBt9DkqwLZvijJgrD2pL1gvvmyHTfBT342oHvm+05FlDMPfIzLKCY7T1agD/DQl2A54DJKbadDDzhe34zcLP3eAUw3Hs8HFjRk3HaNa7umVMgr10I+2arUGLuqXMslJhL/WdYjOdXrOYD40VkrIhUAJcCs71ts4HEVBhXApm04LpMvAxpepiILFDVki5mWOrnaOdX/MrhHLtCRC4GfgEMBnYBi1T1IyIyArhHVc/39jsf+BkQBu5V1Vu99QOBR4AxwPvAJ1R1R0/Fa9e48mdm57sUvVI/Rzu/4lcO55g1VZ0FzEqyfiNwvu/5XGBukv22A2f3ZIx+1uIyxhhTVOwalzHGmKJiiSuHMi17kqZsyvdEZIOILPKW85Mdn2+p4vVtFxG5w9u+WESOzfTYQtDN83tPRN7yfl4LOh5bKDI4x0NF5F8i0iQiX8/m2ELQzfMrip+h8Ql6+GUpLcB/Azd5j28CfpRknzCwGjgIqADeBCZ6274HfD3o88g0Xt8+5wOP4d7ncRLwaqbHBr105/y8be8Bg4I+jxyc4xDgeOBW/3uwhH6GSc+vWH6GtrRfrMWVW5mUPUlXNqUQZRLvdOABdb0C9BOR4RkeG7TunF+x6PQcVXWLqs4HWrI9tgB05/xMEbLElVvtyp7gfsrrKFnZlJG+59d73VH39nSF5Qx1Fm+6fTI5NmjdOT9wKwc8KSILRWRGj0XZPd35OZTKzzCdYvgZGh8bDp+ldGVTMn2JJOsSQzvvAr7vPf8+8BPg89nGmGOZlHlJtU/hlIhJrTvnB3CKqm70arM9JSJvq2rKCtsB6c7PoVR+hukUw8/Q+FjiypKqnpNqm4hsFpHh6haZHA5sSbLbemC07/koYKP32pt9r/Ub4O+5ibpbUsabwT4VGRwbtO6cH+re54KqbhGRWbjdVoX2Ry+Tc+yJY/OlWzEWyc/Q+FhXYW5lUvYkZdmUDtdNLsadIydo6cq8JMwGrvBG350E1HldpZkcG7Qun5+I1IhIHwARqQHOozB+Zh115+dQKj/DpIroZ2j8gh4dUkoLMBB3ErWV3tcB3voRwFzffucD7+COhLrFt/73wFvAYtxfvOFBn1OqeIFrgWu9x4I7wdxqL/7J6Y4ttKWr54c7iu1Nb1laqOeX4TkOw2257MYt+bMe6FtCP8Ok51dMP0Nb2harnGGMMaaoWFehMcaYomKJyxhjTFGxxGWMMaaoWOIyxhhTVCxxGWOMKSqWuIwxxhQVS1zGGGOKiiUuYzoQkR+KyP8LOIbXROTwIGMwplBZ4jIlR0T6i4iKyB4R2Scia0Xk6gyPHQxcAfy6Z6Ps1I+B/wo4BmMKkiUuU4omAdtUtbeq9gJuBn4tIoMyOPZzuOW5GnowvkzMBs4ssnm/jMkLS1ymFE0CXvc9fx53ltzW+c1E5IsiskxE6kTkMW9KC4Cp3v749r1FRO7yPe8vIi0iUtVTJ6CqjcBC3KKvxhgfS1ymFB2D+0cfEekH/NB7vspb9y3cAqwXAYOBDcAPvGOPBFZ0eL0jgUW+55OAFV5y6UnLgaN7+HsYU3QscZlSNAn4qojsBnbizkQ9RVXVa1l9G7hMVVepO9X7b4HjvWP7AfUdXi9Z4noTQEROFpF/icjzIvKQiERFZKiIvOyte0ZEhifbL/FiInKriPxTRP4sIr1836fei8cY42OJy5QUEakEDgOOUtW+wMeBk4AWb5ezcSe4fE1EdonILuBxoM7bvhPo43u9CuBg3OlMEo6mLZGtBc5S1dOBNcB0YBtwqrfuAeDqFPshIkcAB6vqacDTtJ/xug/uFBzGGB9LXKbUHAE04SYHVPUvwPvAx7ztA4BZqtrPt9Sq6hne9sXABN/rTQQ2qOo+ABER4Ay8FpeqbvQN5IgBjqrGVdXx1vUBlibbz3t8GvCY9/gx4FTf9z4s8X2MMW0scZlScwywRNtPNDcX93oWuIM2zhSRYwFEpK+ITPcSUmLf033HHgkMEZGDRaQa+D5wAPCe/5uKyFjcgR1/955PEpFXgevxDRTpuB/ugJFEa68ON7EmWo7HAU914d/AmJJmicuUmkm4rSa/x4FzRaRKVf+Fe3/UX0RkD7AM7/qXt+8DwPlekgI3cT2B2xpaBWzGbc3dknhxEekL3A9c7l0zQ1UXqeqJwHdwh+Mn3Q+3a7LWe1wL7PAeXwQ8p6obu/FvYUxJshmQjelARG4Dtqjqz0TkMeAer8sx2b4R4FHgJ6r6jLeuUlWbvMcfAT4C/HvH/bztRwI3q+qnRWQGUKmqv/Baa1er6pIePFVjipIlLmPSEJH1wHmquizF9suBnwKJBHMX7jW1HwFxoBF3wMU5HfdT1f/zXuOHuNe2tgBXqOrenjkbY0qDJS5jUhCR/rhdgzWq2tLZ/saY/LDEZYwxpqjY4AxjjDFFxRKXMcaYomKJyxhjTFGxxGWMMaaoWOIyxhhTVCxxGWOMKSqWuIwxxhQVS1zGGGOKyv8HdUht/Rn96fIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1,1,figsize=(6,5))\n",
    "t = axs.scatter(mu_re,mu_im,c=spin)\n",
    "cbar = plt.colorbar(mappable=t,ax=axs)\n",
    "cbar.set_label(r'$\\chi$',rotation=0,fontsize=12,labelpad=10)\n",
    "axs.set_xlabel(r'$Re(\\mu_{2320})$',fontsize=12)\n",
    "axs.set_ylabel(r'$Im(\\mu_{2320})$',fontsize=12)\n",
    "axs.minorticks_on()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}