TLEFit - COE - JAX.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "likely-hearing",
   "metadata": {},
   "outputs": [],
   "source": [
    "%config Completer.use_jedi = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "indirect-brighton",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "indian-worry",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "successful-tobacco",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:jax._src.lib.xla_bridge:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
     ]
    }
   ],
   "source": [
    "from tlefit_coe_jax import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "final-resolution",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "national-necessity",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8d4758f-4617-41dc-ba1b-f3a86fda319e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "disabled-juvenile",
   "metadata": {},
   "outputs": [],
   "source": [
    "line1 = '1 25544U 98067A   14020.93268519  .00009878  00000-0  18200-3 0  5082'\n",
    "line2 = '2 25544  51.6498 109.4756 0003572  55.9686 274.8005 15.49815350868473'\n",
    "satellite = EarthSatellite(line1, line2, 'ISS (ZARYA)', ts)\n",
    "\n",
    "line1 = '1 40019U 14033K   21064.48089419  .00000027  00000-0  13123-4 0  9994'\n",
    "line2 = '2 40019  97.7274 245.3630 0083155 314.3836  45.0579 14.67086574359033'\n",
    "satellite = EarthSatellite(line1, line2, 'APRIZESAT 10', ts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "persistent-command",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial semi-major axis (a) = 7054.852 km\n",
      "COE elements (original) = [7048.523352797561, 0.0083155, 1.7056648994135064, 4.282392212570847, 5.487028934272845, 0.786408709034353, 1.3123000000000001e-05]\n",
      "Residuals at epoch time [ 5.60046126e-01 -3.55208784e+00 -1.46641965e+01 -2.73377804e-03\n",
      " -7.52656544e-03  1.66212876e-04]\n",
      "Residual magnitudes at epoch time 15.0987, 0.00800939\n",
      "\n",
      "\n",
      "#################### ITERATION 1 ####################\n",
      "\n",
      "Condition number (A): 50095.081\n",
      "Condition number (ATWA_acc): 13543803.64408722\n",
      "Lambda:  0.001\n",
      "Residuals after/before 196 < 8.2e+06\n",
      "Covariance a: 12.323 m\n",
      "dx  [-6.32427946e+00  2.91129898e-04  8.86588005e-05  1.50245701e-06\n",
      " -1.39220109e-01  1.39216623e-01 -8.27898944e-04]\n",
      "COE elements = [ 7.04852725e+03  8.25455910e-03  1.70566489e+00  4.28239372e+00\n",
      "  5.48210981e+00  7.91329584e-01 -8.26898944e-04]\n",
      "\n",
      "Current Solution 7048.52725359321 0.0082546 97.7274 245.3631 314.1018 45.3398 -8.2690e-04\n",
      "        Original 7048.523352797561 0.0083155 97.7274 245.3630 314.3836 45.0579 +1.3123e-05\n",
      "Residual (b) = [-7.14673731e-01  3.11148753e+00  1.46442084e+01  2.40917232e-03\n",
      "  7.13957301e-03  8.32241505e-04]\n",
      "Residuals (b) r = 15, v = 0.00758\n",
      "Convergence test: 0.998346, sigma_new(82.7) < sigma_old(5e+04)\n",
      "\n",
      "#################### ITERATION 2 ####################\n",
      "\n",
      "Condition number (A): 52688.495\n",
      "Condition number (ATWA_acc): 13460677.647120733\n",
      "Lambda:  0.001\n",
      "Residuals after/before 0.00341 < 196\n",
      "Covariance a: 12.306 m\n",
      "dx  [-3.89747238e-03  6.07801641e-05  1.03595921e-08 -1.50098852e-06\n",
      "  4.95057397e-03 -4.95231764e-03  8.39971132e-04]\n",
      "COE elements = [7.04852336e+03 8.31533927e-03 1.70566490e+00 4.28239221e+00\n",
      " 5.48706038e+00 7.86377266e-01 1.30721882e-05]\n",
      "\n",
      "Current Solution 7048.523356120826 0.0083153 97.7274 245.3630 314.3854 45.0561 +1.3072e-05\n",
      "        Original 7048.523352797561 0.0083155 97.7274 245.3630 314.3836 45.0579 +1.3123e-05\n",
      "Residual (b) = [ 1.80849849e-01  4.73284718e-01  1.65062797e-01 -5.61274612e-05\n",
      "  5.84284210e-05  5.49490803e-04]\n",
      "Residuals (b) r = 0.533, v = 0.000555\n",
      "Convergence test: 0.995116, sigma_new(0.404) < sigma_old(82.7)\n",
      "\n",
      "#################### ITERATION 3 ####################\n",
      "\n",
      "Condition number (A): 52604.441\n",
      "Condition number (ATWA_acc): 13456143.993964372\n",
      "Lambda:  0.001\n",
      "Residuals after/before 3.36e-09 < 0.00341\n",
      "Covariance a: 12.304 m\n",
      "dx  [-3.31994956e-06  1.60570309e-07  6.54888081e-12 -1.47485126e-09\n",
      " -3.14145371e-05  3.14123509e-05  5.07615134e-08]\n",
      "COE elements = [7.04852335e+03 8.31549984e-03 1.70566490e+00 4.28239221e+00\n",
      " 5.48702897e+00 7.86408678e-01 1.31229497e-05]\n",
      "\n",
      "Current Solution 7048.523352800877 0.0083155 97.7274 245.3630 314.3836 45.0579 +1.3123e-05\n",
      "        Original 7048.523352797561 0.0083155 97.7274 245.3630 314.3836 45.0579 +1.3123e-05\n",
      "Residual (b) = [-8.00064704e-04 -3.58253555e-04  4.17990680e-03  9.96707508e-07\n",
      "  2.09324536e-06 -2.51753382e-07]\n",
      "Residuals (b) r = 0.00427, v = 2.33e-06\n",
      "Convergence test: 0.995822, sigma_new(0.00169) < sigma_old(0.404)\n",
      "\n",
      "#################### ITERATION 4 ####################\n",
      "\n",
      "Condition number (A): 52605.026\n",
      "Condition number (ATWA_acc): 13456142.267325355\n",
      "Lambda:  0.001\n",
      "Residuals after/before 3.35e-15 < 3.36e-09\n",
      "Covariance a: 12.304 m\n",
      "dx  [-3.31465667e-09  1.64514835e-10  6.43527917e-15 -1.47288603e-12\n",
      " -3.07760699e-08  3.07738762e-08  5.05440578e-11]\n",
      "COE elements = [7.04852335e+03 8.31550000e-03 1.70566490e+00 4.28239221e+00\n",
      " 5.48702893e+00 7.86408709e-01 1.31230003e-05]\n",
      "\n",
      "Current Solution 7048.5233527975615 0.0083155 97.7274 245.3630 314.3836 45.0579 +1.3123e-05\n",
      "        Original 7048.523352797561 0.0083155 97.7274 245.3630 314.3836 45.0579 +1.3123e-05\n",
      "Residual (b) = [-7.78604317e-07 -3.16637852e-07  4.16414787e-06  9.85177495e-10\n",
      "  2.08363116e-09 -2.04010142e-10]\n",
      "Residuals (b) r = 4.25e-06, v = 2.31e-09\n",
      "Convergence test: 0.999008, sigma_new(1.67e-06) < sigma_old(0.00169)\n",
      "\n",
      "Stopped due to dx convergence\n",
      "\n",
      "#################### SOLUTION IN 4 ITERATIONS ####################\n",
      "\n",
      "Solution 7048.5233527975615 0.0083155 97.7274 245.3630 314.3836 45.0579 +1.3123e-05\n",
      "Original 7048.523352797561 0.0083155 97.7274 245.3630 314.3836 45.0579 +1.3123e-05\n",
      "Residuals (b) r = 4.25e-06, v = 2.31e-09\n",
      "Stopped in 4 iterations\n",
      "Residual at epoch     1.510e+01 km 8.009e-03 km/s\n",
      "Residual at new epoch 4.209e-06 km 0.000e+00 km/s\n",
      "Residual at the end   4.248e-06 km 2.314e-09 km/s\n"
     ]
    }
   ],
   "source": [
    "fit_method = test_tle_fit_normalized\n",
    "# fit_method = test_tle_fit\n",
    "\n",
    "iterations, sigma, sigmas, dxs, bs, lamdas, b_epoch, b_new_epoch, b, P, A = fit_method(satellite, lamda=1e-3, rms_epsilon=0.0001, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca8079b6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "separate-sleeping",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f81087a6080>]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8bUlEQVR4nO3deViVdf7/8dd9Dpso4IKiCO77BoqImooLZlaO1rRZuWdlZRpWX5tf5cw0M85MqWXSZrm1acuoLbYoJi65grjvkksKuCQgKgjn/v2BWZYaKHCf5fm4rvOHhxvO+9zX8Zwn9/3hHMM0TVMAAAAuwmb1AAAAACVBvAAAAJdCvAAAAJdCvAAAAJdCvAAAAJdCvAAAAJdCvAAAAJdCvAAAAJfiZfUApc3hcOjIkSMKCAiQYRhWjwMAAIrBNE3l5OQoNDRUNtvVj624XbwcOXJE4eHhVo8BAACuwaFDhxQWFnbVbdwuXgICAiQV3fnAwECLpwEAAMWRnZ2t8PDwi6/jV+N28fLzqaLAwEDiBQAAF1OcJR8s2AUAAC6FeAEAAC6FeAEAAC6FeAEAAC6FeAEAAC6FeAEAAC6FeAEAAC6FeAEAAC6FeAEAAC6FeAEAAC6FeAEAAC6FeAEAAC6FeCmmk7n5uu/tNdp48CerRwEAwKMRL8U0ZfFurdp7Qre//r2eX7hV2efOWz0SAAAeiXgpprFxjXV7u9oyTWnO6gOKm5SkRVuOyjRNq0cDAMCjEC/FVK2SrybfFakPHohR/eCKyszJ0yPvp2j4rPU6dPKM1eMBAOAxiJcS6twoWF+N6aoxvRrLx27Td7uOqfeUJL2ZtE/nCx1WjwcAgNsjXq6Bn7ddT/RuokVjuiqmflWdO+/QxK92qt+rK5XCgl4AAMoU8XIdGtWopLkPdtRLd0aoir+3dqbn6M+vf69nF2xR1lkW9AIAUBaIl+tkGIbuiApT4rjuuiMqTKYpvbfmoOImJ+mLzUdY0AsAQCkjXkpJ1Yo+eunOCH04sqMaVK+oYzl5euyDjRrGgl4AAEoV8VLKOjWspq/GdNUTcU3kY7dp2YUFva8vY0EvAAClgXgpA75edo2Ja6yvxnZVpwbVdO68Q//5eqdunbpSyQdY0AsAwPUgXspQw+qV9MHIGE26sKB3V0bRgt6/zN+irDMs6AUA4FoQL2XMMAz9OSpMS8d1113twyRJH6w9qF6Tk/TZJhb0AgBQUsRLOalS0Uf/vSNCcx/sqIbVK+r46Tw9/uFGDZm5XgdPsKAXAIDiIl7KWccG1bRoTFeN691EPl42Ld9dtKD3tWV7WdALAEAxEC8W8PWya3SvxvpmbDfd0Kia8goc+u/Xu3TL1BXa8MNJq8cDAMCpES8Wqh9cUe+NiNGUuyNUtaKPdmec1h1vrNYz/9vMgl4AAK6AeLGYYRi6rW2YEuNjdXf7cEnSh+sOqdfkZVqY+iMLegEA+A3ixUlUqeij/9zRRh891EmNalTS8dP5GjM3VYNnrNOBE7lWjwcAgNMgXpxMh/pVtejxrnryxqIFvSv2HNeNU5Yr4bu9yi9gQS8AAMSLE/Lxsumxno317dhu6tIoWHkFDr34TdGC3vUs6AUAeDjixYnVC66od0d00Cv3RCq4ko/2ZJ7WnW+s1vhPN+vUmXyrxwMAwBLEi5MzDEP9I2trSXysBnYoWtA7d/0h9ZqUpPkbD7OgFwDgcYgXF1HZ30cTb2+jjx/upMY1KulEbr6emLdJ97+zVmnHWdALAPAcxIuLia5XVV8+3lVP9WkqXy+bVu09oT4vL9eriXuUV1Bo9XgAAJQ54sUF+XjZ9GiPRvr2iW7q2jhY+QUOTVq8W7dMXam1+09YPR4AAGWKeHFhdatV1Jzhvyzo3Zt5Wne/tUZPf7JJP+WyoBcA4J6IFxf384LexPjuujemjiTpow2H1Wtykj5NZkEvAMD9EC9uIsjfW/+6rbU+HdVJTUMCdDI3X+M+3qT73l6r/cdOWz0eAAClhnhxM1F1q+qLx7vo/25qJj9vm77fd0I3vbxCryxhQS8AwD0QL27I227TqO4N9e3YWHVrUl35hQ5NWbJbfV9ZoTUs6AUAuDjixY3Vqeav2cOi9erAtgqu5Kv9x3J1z1tr9OTHm3SSBb0AABdFvLg5wzDULyJUieNidd+FBb2fJB9Wr0nL9AkLegEALoh48RBBFbz1z9ta69NRndWsZoB+OnNeT368SQOnr9E+FvQCAFyIU8bLbbfdpipVquiOO+6wehS3E1W3ij4f3UXj+xYt6F2z/6T6vrxCUxbv1rnzLOgFADg/p4yXMWPGaM6cOVaP4ba87TY9HNtQi5+IVfemRQt6X0nco5tfWaHv9x23ejwAAK7KKeOle/fuCggIsHoMtxde1V8zh0Yr4d52qh7gq/3Hc3Xv9LUa9xELegEAzqvU42X58uXq16+fQkNDZRiGFixY8LttEhISVK9ePfn5+SkmJkbr1q0r7TFQTIZh6JY2tZQ4LlaDOtaVYUifphxWz0nL9NGGQyzoBQA4nVKPl9zcXEVERCghIeGyX583b57i4+M1YcIEpaSkKCIiQn369FFmZmZpj4ISCPTz1gsDWul/Fxb0njpzXk9/sln3vLVGezNZ0AsAcB6lHi99+/bVP/7xD912222X/frkyZM1cuRIDRs2TC1atNAbb7whf39/zZgx45puLy8vT9nZ2ZdccO3a1ila0PuXm5upgrdda9NOqu8ryzX5210s6AUAOIVyXfOSn5+v5ORkxcXF/TKAzaa4uDitXr36mn7mxIkTFRQUdPESHh5eWuN6LG+7TQ92a6hvn+imHk2r63yhqalL96rvKyu0ai8LegEA1irXeDl+/LgKCwsVEhJyyfUhISFKT0+/+O+4uDjdeeedWrRokcLCwq4aNs8884yysrIuXg4dOlRm83ua8Kr+mjE0Wq/d1041AnyVdjxX9729VvHzUnXidJ7V4wEAPJSX1QNczpIlS4q9ra+vr3x9fctwGs9mGIZubl1LXRoHa9I3uzRnzQH9b+OPWrorU8/0baY7o8JlsxlWjwkA8CDleuQlODhYdrtdGRkZl1yfkZGhmjVrlucoKKFAP2/9rX8rzX/kBjWvFahTZ87r/z7donveWqM9GTlWjwcA8CDlGi8+Pj6KiopSYmLixescDocSExPVqVOn8hwF1ygyvLI+f+wGPXtLc1XwtmvdDyd189QVmsSCXgBAOSn1eDl9+rRSU1OVmpoqSUpLS1NqaqoOHjwoSYqPj9f06dM1e/Zs7dixQ6NGjVJubq6GDRtW2qOgjHjZbXqgawMtju+mXs1q6HyhqVeX7tVNLy/Xyj0s6AUAlC3DLOV3IVu2bJl69Ojxu+uHDBmiWbNmSZKmTZumF198Uenp6YqMjNTUqVMVExNTKrefnZ2toKAgZWVlKTAwsFR+Jq7MNE19sy1dEz7bpozsokW8AyJD9eytLRRcibVIAIDiKcnrd6nHi9WIF2vknDuvSd/u1uzVP8g0iz7F+pm+zXRXexb0AgD+GPFCvFhm06FT+sv8Ldp2pOjNAtvXraJ/3d5aTUL4rCoAwJWV5PXbKT+YEa4rIryyFj5atKDX38euDQd+0s2vrNCL3+xkQS8AoFQQLyh1vyzojVVc8xAVOEwlfLdPN05ZruW7j1k9HgDAxREvKDO1K1fQ20Pa681BUaoZ6KeDJ89o8Ix1evzDjcrMOWf1eAAAF0W8oMz1aVlTS8bFatgN9WQzpM82HVHcpCR9sPagHA63WnIFACgHxAvKRSVfL03o11ILHr1BrWoHKvtcgf4yf4vufHO1dqXzDr0AgOIjXlCu2oRV1oJHbtDzt7ZQRR+7kg/8pFumrtB/vt6ps/ks6AUA/DHiBeXOy27T8C71tTg+Vje2KFrQ+/qyfbrx5SQt25Vp9XgAACdHvMAyoZUr6K3BRQt6awX56dDJsxo6c71Gs6AXAHAVxAss16dlTS2Oj9WILvVlM6TPNx1Rr0lJem/NARb0AgB+h3iBU6jk66Xnbm2hzx7rota1g5RzrkDPLtiqO974XjvTs60eDwDgRIgXOJVWtYO04NEb9Nd+RQt6Uw6e0i1TV2riVzt0Jr/A6vEAAE6AeIHTsdsMDb2hvpaMi9VNLWuq0GHqzaT9unHKcn3Hgl4A8HjEC5xWraAKemNQlKYPbq/QID8d/umshs1cr0c/SFFmNgt6AcBTES9wer1bhGhxfKweuLCg98vNR9VrUpLeZUEvAHgk4gUuoaKvl569sKA3IixIOXkFem7BVt3++vfafoQFvQDgSYgXuJRWtYP0v0du0N/+1FKVfL2UeuiU+k1bqYmLWNALAJ6CeIHLsdsMDelcT0viY3Vz6wsLepfvV+/Jy7V0Z4bV4wEAyhjxApdVM8hPr90XpXeGtFftyhX046mzGj5rgx55P1kZLOgFALdFvMDl9WoeosXx3fRgtway2wwt2pKuXpOSNGf1DypkQS8AuB3iBW7B38dLf7m5uT577AZFhFfW6bwCPb9wm25/bZW2HcmyejwAQCkiXuBWWoYG6X+jOuuF/i0V4OulTYez9Kdpq/TPL7crN48FvQDgDogXuB27zdCgTvW0ZFysbmldS4UOU9NXpOnGKcu1ZDsLegHA1REvcFshgX5KuK+dZg6Nvrig94E5G/Twu8lKz2JBLwC4KuIFbq9HsxpaHN9ND8UWLej9elu64iYnadaqNBb0AoALIl7gEfx9vPRM3+b6YnQXRV5Y0PvXz7frttdWaeuPLOgFAFdCvMCjNK8VqP+N6qx/DGilAD8vbT6cpT9NW6kXvmBBLwC4CuIFHsdmM3R/x7pKjI/VrW1qyWFK76xMU+/JSVrMgl4AcHrECzxWjUA/Tbu3nWYOi1ZYlQo6knVOI+ds0INzNujIqbNWjwcAuALiBR6vR9MaWvxErB6ObSgvm6Fvt2eo9+QkzVjJgl4AcEbECyCpgo9d4/s20xePd1G7OpWVm1+ov3+xXQMSVmnLYRb0AoAzIV6AX2lWM1CfPNxZ/7ytaEHvlh+z1D9hpf72+TadZkEvADgF4gX4DZvN0H0xdZU4LlZ/igiVw5RmrvpBvScn6Ztt6VaPBwAej3gBrqBGgJ+mDmyr2cM7KLxqBR3NOqeH3k3WSBb0AoCliBfgD8Q2qa5vx8bqke5FC3oXb89Q3OQkvb1ivwoKHVaPBwAeh3gBiqGCj11P39RMXz7eVVF1q+hMfqH+8eUO9U9Ypc2HT1k9HgB4FOIFKIGmNQP08UOdNPH21gr089K2I9kakLBKf/1sm3LOnbd6PADwCMQLUEI2m6GBHeoocVx39Y8sWtA76/sfFDc5SV9vPSrT5L1hAKAsES/ANaoe4KtX7mmrOcM7qG41f2Vk5+nh91I0cs4G/ciCXgAoM8QLcJ26Namub8Z202M9GsnbbmjJjkz1ZkEvAJQZ4gUoBX7edj3Zp6m+fLyrouv9sqD3T9NWKfXQKavHAwC3QrwApahJSIDmPdhJ//lzawVV8Nb2o9m67bVVmrBwq7JZ0AsApYJ4AUqZzWbo7ug6ShwXq9va1pZpSrNXH1DvyUlatIUFvQBwvYgXoIwEV/LVlLsj9d6IGNW7sKD3kfdTNGL2Bh06ecbq8QDAZREvQBnr0jhYX4/tptE9ixb0Lt2ZqRunLNdby/fpPAt6AaDEiBegHPh52zXuxqb6akxXdahXVWfPF+pfi3bqT9NWaePBn6weDwBcCvEClKNGNQI098GO+u+f26iyv7d2HM3W7a9/r2cXbFHWGRb0AkBxEC9AObPZDN0VHa7E+Fjd3q5oQe97aw6q56Rl+iT5MAt6AeAPEC+ARapV8tXkuyL1wQMxali9ok7k5uvJjzfprjdXa8fRbKvHAwCnZZhu9mtedna2goKClJWVpcDAQKvHAYolv8Chd1amaWriHp09Xyi7zdCQTvX0RO/GCvDztno8AChzJXn95sgL4AR8vGwa1b2hloyLVd9WNVXoMDVjVZp6TUrSwtQfOZUEAL9CvABOpHblCnr9/ijNGhatetX8lZmTpzFzU3Xf22u1NzPH6vEAwCkQL4AT6t60hr4e203xvZvI18um7/edUN9XVujfX+3UmfwCq8cDAEsRL4CT8vO26/FejbX4iVj1alZD5wtNvZG0T3GTkvQVHzMAwIMRL4CTq1PNX+8Mjdb0we1Vu3IFHck6p1Hvp2jIzPVKO55r9XgAUO6IF8BF9G4RoiXxsXqsRyP52G1avvuY+kxZrsnf7tK584VWjwcA5YZ4AVxIBR+7nuzTVF+P7aqujYOVX+jQ1KV71XtKkhJ3ZFg9HgCUC+IFcEENqlfSnOEdlHBvO9UM9NOhk2c1YvYGPcAnVgPwAMQL4KIMw9AtbWppybhYPditgbxshpbsyFDvKUmatnSP8go4lQTAPfEOu4Cb2J2Ro+cWbNXatJOSpAbBFfW3/i3VtXF1iycDgD/GO+wCHqhJSNEnVr98d6SCK/lq//FcDXpnnR59P0VHs85aPR4AlBriBXAjhmFoQNvaWvpkrIZ2riebIX255ah6TUrSm0n7dL7QYfWIAHDdOG0EuLFtR7L03IKtSjl4SpLUuEYlvTCglTo2qGbtYADwG5w2AiBJahkapE8e7qz//rmNqlb00Z7M07rnrTUaO3ejMnPOWT0eAFwT4gVwczabobuiw7V0XKzujakjw5AWpB5Rr5eSNHNVmgo4lQTAxXDaCPAwmw6d0nMLt2rz4SxJUotagXphQCtF1a1i8WQAPFlJXr+JF8ADFTpMfbjuoF78Zpeyzp6XJN3VPkz/d1MzVavka/F0ADwRa14AXJXdZuj+jnW1dFys7owKkyR9tOGwek5K0vtrD6jQ4Va/0wBwMxx5AaANP5zUswu2amd6jiQpIixILwxopTZhla0dDIDH4LQR8QKUWEGhQ3NWH9Dkxbt1Oq9AhiHdF1NHT93YTEH+3laPB8DNcdoIQIl52W0a3qW+lo6LVf/IUJmm9N6ag+oxaZk+2nBIDk4lAXASHHkBcFnf7zuu5xdu097M05Kk9nWr6O/9W6lFKP+vAJQ+ThsRL0CpyC9waOaqNL2SuEdn8gtltxka3Kmu4ns3UYAfp5IAlB5OGwEoFT5eNj0U21BL4mN1c+uaKnSYmrnqB/WclKSFqT/KzX73AeAiiBcAfyi0cgW9dl+UZg/voPrBFXUsJ09j5qbq3ulrtScjx+rxAHgYp4yX2267TVWqVNEdd9xh9SgAfiW2SXV9PbarxvVuIl8vm1bvP6G+r6zQxK92KDevwOrxAHgIp4yXMWPGaM6cOVaPAeAyfL3sGt2rsZbExyqueQ0VOEy9mbRfcZOT9NWWo5xKAlDmnDJeunfvroCAAKvHAHAV4VX99faQaL09uL3CqlTQ0axzGvV+iobMXK+047lWjwfAjZU4XpYvX65+/fopNDRUhmFowYIFv9smISFB9erVk5+fn2JiYrRu3brSmBWAE4prEaLFT8RqdM9G8rHbtHz3MfWZslyTv92lc+cLrR4PgBsqcbzk5uYqIiJCCQkJl/36vHnzFB8frwkTJiglJUURERHq06ePMjMzL24TGRmpVq1a/e5y5MiRa78nACxTwceucTc21TdPdFPXxsHKL3Ro6tK9ipucpCXbM6weD4Cbua73eTEMQ/Pnz9eAAQMuXhcTE6Po6GhNmzZNkuRwOBQeHq7Ro0dr/Pjxxf7Zy5Yt07Rp0/TJJ59cdbu8vDzl5eVd/Hd2drbCw8N5nxfAIqZp6uut6fr7F9t1NOucJCmueQ1N6NdS4VX9LZ4OgLOy7H1e8vPzlZycrLi4uF9uwGZTXFycVq9eXZo3ddHEiRMVFBR08RIeHl4mtwOgeAzDUN/WtbQkPlYPxTaQl83Qkh2ZipucpFcT9yivgFNJAK5PqcbL8ePHVVhYqJCQkEuuDwkJUXp6erF/TlxcnO68804tWrRIYWFhVw2fZ555RllZWRcvhw4duub5AZSeir5eeqZvc301pqs6NqiqvAKHJi3erZteXqHlu49ZPR4AF+Zl9QCXs2TJkmJv6+vrK19f3zKcBsD1aBwSoA9HdtRnm47oH1/uUNrxXA2esU43t66p525toVpBFaweEYCLKdUjL8HBwbLb7crIuHSBXkZGhmrWrFmaNwXAhRiGof6RtbV0XKyG31BfdpuhRVvS1WtSkt5M2qf8AofVIwJwIaUaLz4+PoqKilJiYuLF6xwOhxITE9WpU6fSvCkALijAz1vP92uhzx/rovZ1q+hMfqEmfrVTN09dodX7Tlg9HgAXUeJ4OX36tFJTU5WamipJSktLU2pqqg4ePChJio+P1/Tp0zV79mzt2LFDo0aNUm5uroYNG1aqgwNwXS1CA/XRQ5304h1tVK2ij/ZmntbA6Ws0Zu5GZWafs3o8AE6uxH8qvWzZMvXo0eN31w8ZMkSzZs2SJE2bNk0vvvii0tPTFRkZqalTpyomJqZUBv4jJflTKwDWyzpzXi9+u1Pvrz0o05QCfL30RO8mGtyprrzsTvkm4ADKQElev6/rfV6cEfECuKbNh0/puQVbtelwliSpea1A/WNAS0XVrWrxZADKA/FCvAAuqdBhau76g/rv17uUdfa8JOnOqDCN79tM1SrxV4WAO7PsTeoA4HrYbYbui6mrpeNidVf7MEnSx8mH1XNSkt5bc0CFDrf6XQvANeLICwCnlXzgpJ5dsE07jmZLktqEBemF/q0UEV7Z2sEAlDpOGxEvgNsoKHTo3TUHNPnb3crJK5BhSAM71NHTfZqqsr+P1eMBKCWcNgLgNrzsNg27ob4Sn4zVbW1ryzSlD9YeVM9JSfpo/SE5OJUEeByOvABwKWv2n9BzC7ZqT+ZpSVJU3Sr6e/+WahkaZPFkAK4Hp42IF8CtnS90aOaqNL28ZI/O5BfKZkiDO9VT/I1NFOjnbfV4AK4Bp40AuDVvu00PdmuoxHGxuqV1LTlMadb3P6jXpCQt2Pij3Ox3MgC/QbwAcFm1gioo4b52mjO8g+oHV9SxnDyNnZeqgdPXaE9GjtXjASgjxAsAl9etSXV9Pbarnryxify8bVqz/6T6vrJCExftUG5egdXjAShlxAsAt+DrZddjPRtr8ROximseogKHqTeX71fc5CQt2nKUU0mAGyFeALiV8Kr+entIe70zpL3Cq1bQ0axzeuT9FA2esU77j522ejwApYB4AeCWejUP0eInYvV4r8by8bJpxZ7juunlFXrpm106m19o9XgArgPxAsBt+XnbFd+7ib4d202xTaorv9Chad/tVe8pSVq8PcPq8QBcI+IFgNurF1xRs4ZF64372yk0yE+HfzqrkXM2aMSs9Tp08ozV4wEoIeIFgEcwDEM3taqlJeNi9XBsQ3nZDCXuzFTc5CRNTdyjc+c5lQS4CuIFgEfx9/HS+L7N9PXYrurcsJryChyavHi3bnp5uZJ2H7N6PADFQLwA8EiNagTo/QdiNHVgW9UI8NUPJ85oyIx1GvVeso6cOmv1eACugngB4LEMw9CfIkKVOC5WI7rUl91m6Kut6YqbnKQ3kvYpv8Bh9YgALoMPZgSAC3YczdbzC7dq/Q8/SZIa1aikv/dvqc4Ngy2eDHB/fDAjAFyD5rUC9dFDnfTSnRGqVtFHezNP697pa/X4hxuVmX3O6vEAXEC8AMCvGIahO6LCtHRcdw3qWFeGIX226Yh6TkrSOyvTVFDIqSTAapw2AoCr2HI4S88u3KpNh05JkprVDNA/BrRS+3pVrR0McDOcNgKAUtI6LEjzR3XWv25rrcr+3tqZnqM73litJz/epOOn86weD/BIxAsA/AGbzdC9MXW0dFx33d0+XJL0SfJh9Xxpmd5dc0CFDrc6gA04PU4bAUAJpRz8Sc8t2KptR7IlSW3CgvRC/1aKCK9s7WCACyvJ6zfxAgDXoNBh6r01B/TSt7uUc65AhiEN7FBHT/dpqsr+PlaPB7gc1rwAQBmz2wwN6VxPS8d11+1ta8s0pQ/WHlSPl5Zp3vqDcnAqCSgzHHkBgFKwdv8JPbdwq3ZnnJYktatTWS8MaKWWoUEWTwa4Bk4bES8ALHC+0KHZ3/+gKYt3Kze/UDZDGtypnuJvbKJAP2+rxwOcGqeNAMAC3nabHujaQInjuuvWNrXkMKVZ3/+gni8laf7Gw3Kz3xUBy3DkBQDKyMo9x/X8Z1u1/1iuJCmmflW9MKCVmoQEWDwZ4Hw48gIATqBL42B9NaarnurTVH7eNq1NO6mbX1mhf365XafzCqweD3BZxAsAlCFfL7se7dFIS+JjdWOLEBU4TE1fkaa4SUn6YvMRTiUB14B4AYByEFbFX28Nbq+ZQ6NVp6q/0rPP6bEPNmrwjHXad+y01eMBLoV4AYBy1KNZDX37RDeN6dVYPl42rdhzXDe9vFwvfrNTZ/MLrR4PcAnECwCUMz9vu57o3USLn+im7k2r63yhqYTv9ilucpK+3ZbOqSTgDxAvAGCRutUqaubQaL05KEq1K1fQj6fO6sF3kzVi9gYdPHHG6vEAp0W8AICFDMNQn5Y1tTi+mx7p3lDedkNLd2aq95QkvbJkj86d51QS8FvECwA4AX8fLz19UzN9NaabbmhUTXkFDk1Zslt9Xl6uZbsyrR4PcCrECwA4kUY1Kum9ETF6dWBbhQT66sCJMxo6c70efjdZP546a/V4gFMgXgDAyRiGoX4RoUoc110PdKkvu83Q19vSFTcpSa8t26v8AofVIwKW4uMBAMDJ7UzP1vMLtmndDyclSQ2rV9QL/Vupc6NgiycDSg8fDwAAbqRZzUDNe6ijJt0ZoeBKPtp3LFf3vr1Woz/cqIzsc1aPB5Q74gUAXIBhGPpzVJgSx3XX4E51ZTOkzzcdUa9JSXp7xX4VFHIqCZ6D00YA4IK2/pilZxdsVeqhU5KkZjUD9MKAVoquV9XawYBrxGkjAHBzrWoH6X+jOuvft7dWFX9v7UzP0Z1vrNa4jzbp+Ok8q8cDyhTxAgAuymYzdE+HOlo6rrsGdgiXYUifphxWz5eW6d3VP6jQ4VYH1oGLOG0EAG5i48Gf9NzCrdr6Y7YkqU1YkN4a1F41g/wsngz4Y5w2AgAP1LZOFS18tIte6N9SAX5e2nw4SwOnr+EvkuB2iBcAcCN2m6FBnerpqzFdFValgtKO52rgWwQM3AvxAgBuKKyKvz4c2VG1K1fQ/uO5Gjh9jTIJGLgJ4gUA3FR4VX/NffBCwBy7EDA5BAxcH/ECAG4svGrREZjQIL+id+advlbHcvhTarg24gUA3Fydav768MGOqhXkp72Zp3Xv9DUEDFwa8QIAHqButYqa+2BH1Qz0054LAcOb2cFVES8A4CEuFzAnCBi4IOIFADxIveCK+vDBjgoJ9NXujNO6d/paAgYuh3gBAA9TP7iiPhzZUTUCfLUrI0f3vb1WJ3PzrR4LKDbiBQA8UIPqlfThg0UBszM9R/dOX6OfCBi4COIFADxUw+qV9MHIjqr+c8C8vZaAgUsgXgDAgzWqUUkfjuyo4Eq+2nE0W/e9vVanzhAwcG7ECwB4uEY1KmnugzEKruSj7QQMXADxAgBQoxoBF47A+GjbkWzd/85aZZ05b/VYwGURLwAASVLjkAB9MLKjqlX00dYfCRg4L+IFAHBRkwsBU7Wij7b8mKVBM9Yq6ywBA+dCvAAALtG0ZoA+GBmjqhV9tPlwlga/s1bZ5wgYOA/iBQDwO81qBur9B2JUxd9bmw5nadA76wgYOA3iBQBwWc1rBer9BzoWBcyhUxr8zjrlEDBwAsQLAOCKWoQG6r0HYlTZ31uph05p8AwCBtYjXgAAV9UyNEjvjYhRUAVvbTx4SkNmrNPpvAKrx4IHI14AAH+oVe0gvf9AUcCkEDCwGPECACiWVrWLjsAE+nkp+cBPGkrAwCLECwCg2FqHBem9B4oCZsOBnzR85nrlEjAoZ8QLAKBE2oRV1rsjYhTg56V1P5zUsFkEDMoX8QIAKLGI8AsB4+uldWknNXzWep3JJ2BQPpwuXk6dOqX27dsrMjJSrVq10vTp060eCQBwGZHhlTVnRAcF+HppLQGDcmSYpmlaPcSvFRYWKi8vT/7+/srNzVWrVq20YcMGVatWrVjfn52draCgIGVlZSkwMLCMpwUApBz8SYPfKVq826lBNc0YGq0KPnarx4KLKcnrt9MdebHb7fL395ck5eXlyTRNOVlfAQB+pV2dKpo9vIMq+Xpp9f4TGjF7vc7mF1o9FtxYieNl+fLl6tevn0JDQ2UYhhYsWPC7bRISElSvXj35+fkpJiZG69atK9FtnDp1ShEREQoLC9NTTz2l4ODgko4JAChHUXWraPbwaFX0sev7fSf0wJz1OneegEHZKHG85ObmKiIiQgkJCZf9+rx58xQfH68JEyYoJSVFERER6tOnjzIzMy9u8/N6lt9ejhw5IkmqXLmyNm3apLS0NH3wwQfKyMi4xrsHACgvUXWravbwDqroY9eqvSc0cs4GAgZl4rrWvBiGofnz52vAgAEXr4uJiVF0dLSmTZsmSXI4HAoPD9fo0aM1fvz4Et/GI488op49e+qOO+647Nfz8vKUl5d38d/Z2dkKDw9nzQsAWGT9Dyc1ZMY6nckvVNfGwZo+uL38vFkDg6uzbM1Lfn6+kpOTFRcX98sN2GyKi4vT6tWri/UzMjIylJOTI0nKysrS8uXL1bRp0ytuP3HiRAUFBV28hIeHX9+dAABcl+h6VTVrWAf5+9i1Ys9xPfhuMkdgUKpKNV6OHz+uwsJChYSEXHJ9SEiI0tPTi/UzDhw4oK5duyoiIkJdu3bV6NGj1bp16ytu/8wzzygrK+vi5dChQ9d1HwAA169D/aqaOTRaFbztWr77mB4iYFCKvKwe4Lc6dOig1NTUYm/v6+srX1/fshsIAHBNYhpU08xh0Ro2c72Sdh/Tw+8l681BUfL14hQSrk+pHnkJDg6W3W7/3QLbjIwM1axZszRvCgDgAjpeeN8XP2+blu06poffTVZeAUdgcH1KNV58fHwUFRWlxMTEi9c5HA4lJiaqU6dOpXlTAAAX0alhNc0YUhQw3+06plHvpRAwuC4ljpfTp08rNTX14qmdtLQ0paam6uDBg5Kk+Ph4TZ8+XbNnz9aOHTs0atQo5ebmatiwYaU6OADAdXRuFKx3hkTL18umpTsz9ej7KcovcFg9FlxUif9UetmyZerRo8fvrh8yZIhmzZolSZo2bZpefPFFpaenKzIyUlOnTlVMTEypDPxH+HgAAHBeK/cc14jZ65VX4FBc8xC9dl87+Xg53Zu9wwIlef12us82ul7ECwA4txV7jumB2RuUV+BQ7xYhSriXgIGLf7YRAMC9dW1cXdMHt5ePl02Lt2do9IcpOl/IKSQUH/ECACh33Zr8EjDfbMvQ6A82EjAoNuIFAGCJ2CbV9dagKPnYbfp6W7oe/5CAQfEQLwAAy3RvWkNvXgiYr7ama8xcAgZ/jHgBAFiqR7MaemNQO/nYbVq0JV1j56WqgIDBVRAvAADL9WwWotfvbydvu6EvNx8lYHBVxAsAwCn0ah6i1++Lkrfd0Bebj+qJjzYRMLgs4gUA4DTiWoTotQsB8/mmIxr3MQGD3yNeAABO5ec3rvOyGVqYekRPfrxJhQ63ej9VXCfiBQDgdG5sWVPTLgTMAgIGv0G8AACc0k2tamravW1ltxmav/FHPUXA4ALiBQDgtG5qVUvTBhYFzP82/qinP9lMwIB4AQA4t76ta2nqPUUB82nKYY3/dLMcBIxH87J6AAAA/sgtbWrJlKkxc1P1cfJhGYb079vbyGYzrB4NFuDICwDAJdzaJlQv3x0pmyF9tOGwnvnfFo7AeCjiBQDgMvpFhGrKhYCZt+GQ/jKfgPFExAsAwKX0j6x9MWDmrj+k/7dgKwHjYYgXAIDL6R9ZW5PvKgqYD9cd1LMLCRhPQrwAAFzSgLa1NemuCBmG9MHag3r+s60yTQLGExAvAACXdVvbML10R1HAvLfmoJ5fuI2A8QDECwDApf05KkwvXgiYd9cc0F8/I2DcHfECAHB5d0SF6b9/biPDkGavPqC/fb6dgHFjxAsAwC3c2T5c/7m9jSRp1vc/6O9fEDDuingBALiNu6LD9Z8/t5YkzVz1g174YgcB44aIFwCAW7k7uo4m3l4UMDNWpekfXxIw7oZ4AQC4nYEd6uhftxUFzDsr0/SvRQSMOyFeAABu6d6YOvrnba0kSdNXpOnfX+0kYNwE8QIAcFv3xdTVCwOKAubN5fv1768JGHdAvAAA3NqgjnX1Qv+WkqQ3k/brv9/sImBcHPECAHB7gzrV09/+VBQwry/bpxcJGJdGvAAAPMKQzvX0134tJEmvLdunSd/uJmBcFPECAPAYQ2+or+dvLQqYad/t1ZTFBIwrIl4AAB5leJf6eu5CwExdulcvL9lj8UQoKeIFAOBxRnSpr2dvaS5JeiVxj15estviiVASxAsAwCM90LWB/t/NRQHz8pI9eoUjMC6DeAEAeKyR3RroLzc3kyRNWbJbryYSMK6AeAEAeLQHuzXU+L5FATNp8W5NW0rAODviBQDg8R6Obainb2oqSXrp291K+G6vxRPhaogXAAAkPdK9kZ7qUxQwL36zS68v22fxRLgS4gUAgAse7dFIT97YRJL0n6936o0kAsYZES8AAPzKYz0ba1zvooD591c79dZyAsbZEC8AAPzG6F6N9URcUcD8a9FOTV++3+KJ8GvECwAAlzEmrrHGxjWWJP1z0Q69vYKAcRbECwAAVzA2roke71UUMP/4cofeWZlm8USQiBcAAK7qibjGGt2zkSTphS+2a+YqAsZqxAsAAFdhGIbiezfRYz2KAuZvn2/XLALGUsQLAAB/wDAMjbuxiR7p3lCS9NfPt2vO6h+sHcqDES8AABSDYRh6qk9TjboQMM8v3KZ3CRhLEC8AABSTYRh6uk9TPRTbQJL03MJtenfNAYun8jzECwAAJWAYhsbf1EwPdbsQMAu26v21BEx5Il4AACghwzA0vm8zjexaX5L0/+Zv1QdrD1o8lecgXgAAuAaGYegvNzfXiC5FAfOX+Vs0dx0BUx6IFwAArpFhGHr2luYafkNRwIz/3xbNW0/AlDXiBQCA62AYhp67tbmG3VBPUlHAfLT+kLVDuTniBQCA62QYhp6/tYWGdq4n05T+73+b9fEGAqasEC8AAJQCwzA0oV8LDelUV6YpPf3pZn2SfNjqsdwS8QIAQCkxDEN//VNLDepYFDBPfbJJnxIwpY54AQCgFBmGob/3b6n7O9aRaUpPfrJJ8zcSMKWJeAEAoJQZhqG//6mV7o0pCphxH23Sgo0/Wj2W2yBeAAAoAzaboX/0b6WBHerIYUrxH6VqYSoBUxqIFwAAyojNZuifA1rpnuhwOUzpiXmp+mzTEavHcnnECwAAZchmM/Sv21rr7vZFATN27kZ9TsBcF+IFAIAyZrMZmnh7a93VPqwoYOal6ovNBMy18rJ6AAAAPIHNZujft7eRw5Q+ST6sMXNTZcjQLW1qWT2ayyFeAAAoJzabof/8uY1MU/o05bAen7tRNkPq25qAKQlOGwEAUI7sNkP/vaONbm9bW4UOU6M/3Kivtx61eiyXQrwAAFDO7DZDL94Zodva1laBw9RjH2zU11vTrR7LZRAvAABYwG4z9NKdERoQGXohYFL07TYCpjiIFwAALGK3GZp0V6T6XwiYRz9I0eLtGVaP5fSIFwAALGS3GZp0Z4T6RYTqfKGpR95P1hIC5qqIFwAALOZlt2nKXRG6tU0tnS80Ner9ZCXuIGCuhHgBAMAJeNltevnuSN3S+kLAvJeipTsJmMshXgAAcBJedptevidSN7euqfxChx5+N0Xf7cq0eiynQ7wAAOBEvO02vXJPW/VtVRQwD72brGUEzCWIFwAAnIy33aapA9vqppY1lV/g0IPvJitp9zGrx3IaxAsAAE7I227Tq/e2VZ+WIcovcGjknA1aTsBIctJ4qVevntq0aaPIyEj16NHD6nEAALCEt92mVwe2U+8WvwTMij0EjFPGiyR9//33Sk1N1XfffWf1KAAAWMbHy6aEe9sprnmI8gocemD2Bq3cc9zqsSzltPECAACK+HjZ9Np97RTXvIbyChwaMXu9Vu313IApcbwsX75c/fr1U2hoqAzD0IIFC363TUJCgurVqyc/Pz/FxMRo3bp1JboNwzAUGxur6Ohovf/++yUdEQAAt+PjZVPCfe3Us9kvAfP9Ps8MmBLHS25uriIiIpSQkHDZr8+bN0/x8fGaMGGCUlJSFBERoT59+igz85c/84qMjFSrVq1+dzly5IgkaeXKlUpOTtZnn32mf/3rX9q8efM13j0AANyHr5ddr9/fTj2aVte58w4Nn7Veq/edsHqscmeYpmle8zcbhubPn68BAwZcvC4mJkbR0dGaNm2aJMnhcCg8PFyjR4/W+PHjS3wbTz31lFq2bKmhQ4de9ut5eXnKy8u7+O/s7GyFh4crKytLgYGBJb49AACc3bnzhXr4vWQt23VMFbztmjksWh0bVLN6rOuSnZ2toKCgYr1+l+qal/z8fCUnJysuLu6XG7DZFBcXp9WrVxfrZ+Tm5ionJ0eSdPr0aS1dulQtW7a84vYTJ05UUFDQxUt4ePj13QkAAJycn7ddb9wfpdgm1XX2fKGGzVyvtfs95whMqcbL8ePHVVhYqJCQkEuuDwkJUXp6erF+RkZGhrp06aKIiAh17NhRgwcPVnR09BW3f+aZZ5SVlXXxcujQoeu6DwAAuAI/b7veHBSlro2DiwJm1nqtSztp9VjlwsvqAX6rQYMG2rRpU7G39/X1la+vbxlOBACAc/Lztmv64PYX3v/luIbOXKfZwzsoul5Vq0crU6V65CU4OFh2u10ZGZd+CmZGRoZq1qxZmjcFAAD0S8B0aRSsM/mFGjpjnTb84N5HYEo1Xnx8fBQVFaXExMSL1zkcDiUmJqpTp06leVMAAOCCnwPmhkbVlJtfqCEz1in5gPsGTInj5fTp00pNTVVqaqokKS0tTampqTp48KAkKT4+XtOnT9fs2bO1Y8cOjRo1Srm5uRo2bFipDg4AAH5RwceutwdHq3PDnwNmvZIP/GT1WGWixH8qvWzZsst+3tCQIUM0a9YsSdK0adP04osvKj09XZGRkZo6dapiYmJKZeA/UpI/tQIAwN2czS8sev+X/SdUyddLc0Z0ULs6Vawe6w+V5PX7ut7nxRkRLwAAT3cmv0DDZ63Xmv0nFXAhYNo6ecBY9j4vAADAev4+XpoxNFod6ldVTl6BBr+zTqmHTlk9VqkhXgAAcEP+Pl6aOTRaHeoVBcygd9Zqk5sEDPECAICbqujrpZnDohVdr4pyzhXo/nfWavPhU1aPdd2IFwAA3FhRwHRQ+7oXAubttdr6Y5bVY10X4gUAADdXyddLs4Z3UFTdKso+V6D7XDxgiBcAADxAJV8vzRoWrXZ1Kivr7HmXDhjiBQAADxHg563Zwzuo7YWAuf+dtdp2xPUChngBAMCD/BwwkeGVdepM0RGY7UeyrR6rRIgXAAA8TKCft+aM6KCIiwGzRjuOuk7AEC8AAHigQD9vzRneQRFhQfrpwhGYnemuETDECwAAHiqogrfmjIhRm7AgnczN173T12pXeo7VY/0h4gUAAA8WVMFb7w6PUevaPwfMGu3OcO6AIV4AAPBwQf7eem9EjFrVDtSJCwGzx4kDhngBAAAXA6ZlaKCOn87XwOlrtTfTOQOGeAEAAJKkyv4+ev+BGLWoFajjp/N0z1trtTfztNVj/Q7xAgAALvo5YJpfCJiB09c4XcAQLwAA4BJVKhYFTLOaATqWUxQw+445T8AQLwAA4HeqVvTRByM7/hIwb63RficJGOIFAABcVtULR2CahgQo88IRmLTjuVaPRbwAAIArq1bJV++PjFGTkErKyC46AvODxQFDvAAAgKsKruSrD0Z2VOMalZSefU5Pf7pZpmlaNg/xAgAA/tDPAdOrWQ1NuTtShmFYNouXZbcMAABcSvUAX70zNNrqMTjyAgAAXAvxAgAAXArxAgAAXArxAgAAXArxAgAAXArxAgAAXArxAgAAXArxAgAAXArxAgAAXArxAgAAXArxAgAAXArxAgAAXArxAgAAXIrbfaq0aZqSpOzsbIsnAQAAxfXz6/bPr+NX43bxkpOTI0kKDw+3eBIAAFBSOTk5CgoKuuo2hlmcxHEhDodDR44cUUBAgAzDKNWfnZ2drfDwcB06dEiBgYGl+rPdDfuq+NhXxce+Kj72VfGxr0qmrPaXaZrKyclRaGiobLarr2pxuyMvNptNYWFhZXobgYGBPMCLiX1VfOyr4mNfFR/7qvjYVyVTFvvrj464/IwFuwAAwKUQLwAAwKUQLyXg6+urCRMmyNfX1+pRnB77qvjYV8XHvio+9lXxsa9Kxhn2l9st2AUAAO6NIy8AAMClEC8AAMClEC8AAMClEC8AAMClEC+/kZCQoHr16snPz08xMTFat27dVbf/+OOP1axZM/n5+al169ZatGhROU1qvZLsq1mzZskwjEsufn5+5TitNZYvX65+/fopNDRUhmFowYIFf/g9y5YtU7t27eTr66tGjRpp1qxZZT6nsyjp/lq2bNnvHleGYSg9Pb18BrbIxIkTFR0drYCAANWoUUMDBgzQrl27/vD7PPH56lr2lac+X0nS66+/rjZt2lx8A7pOnTrpq6++uur3WPG4Il5+Zd68eYqPj9eECROUkpKiiIgI9enTR5mZmZfd/vvvv9fAgQM1YsQIbdy4UQMGDNCAAQO0devWcp68/JV0X0lF78Z49OjRi5cDBw6U48TWyM3NVUREhBISEoq1fVpamm655Rb16NFDqampGjt2rB544AF98803ZTypcyjp/vrZrl27Lnls1ahRo4wmdA5JSUl69NFHtWbNGi1evFjnz5/XjTfeqNzc3Ct+j6c+X13LvpI88/lKksLCwvTvf/9bycnJ2rBhg3r27Kn+/ftr27Ztl93esseViYs6dOhgPvrooxf/XVhYaIaGhpoTJ0687PZ33XWXecstt1xyXUxMjPnQQw+V6ZzOoKT7aubMmWZQUFA5TeecJJnz58+/6jZPP/202bJly0uuu/vuu80+ffqU4WTOqTj767vvvjMlmT/99FO5zOSsMjMzTUlmUlLSFbfx5OerXyvOvuL56lJVqlQx33777ct+zarHFUdeLsjPz1dycrLi4uIuXmez2RQXF6fVq1df9ntWr159yfaS1KdPnytu7y6uZV9J0unTp1W3bl2Fh4dfteQ9mac+pq5XZGSkatWqpd69e2vVqlVWj1PusrKyJElVq1a94jY8tooUZ19JPF9JUmFhoebOnavc3Fx16tTpsttY9bgiXi44fvy4CgsLFRIScsn1ISEhVzx/np6eXqLt3cW17KumTZtqxowZWrhwod577z05HA517txZhw8fLo+RXcaVHlPZ2dk6e/asRVM5r1q1aumNN97Qp59+qk8//VTh4eHq3r27UlJSrB6t3DgcDo0dO1Y33HCDWrVqdcXtPPX56teKu688/flqy5YtqlSpknx9ffXwww9r/vz5atGixWW3tepx5XafKg3n1KlTp0vKvXPnzmrevLnefPNNvfDCCxZOBlfWtGlTNW3a9OK/O3furH379mnKlCl69913LZys/Dz66KPaunWrVq5cafUoTq+4+8rTn6+aNm2q1NRUZWVl6ZNPPtGQIUOUlJR0xYCxAkdeLggODpbdbldGRsYl12dkZKhmzZqX/Z6aNWuWaHt3cS376re8vb3Vtm1b7d27tyxGdFlXekwFBgaqQoUKFk3lWjp06OAxj6vHHntMX3zxhb777juFhYVddVtPfb76WUn21W952vOVj4+PGjVqpKioKE2cOFERERF65ZVXLrutVY8r4uUCHx8fRUVFKTEx8eJ1DodDiYmJVzzX16lTp0u2l6TFixdfcXt3cS376rcKCwu1ZcsW1apVq6zGdEme+pgqTampqW7/uDJNU4899pjmz5+vpUuXqn79+n/4PZ762LqWffVbnv585XA4lJeXd9mvWfa4KtPlwC5m7ty5pq+vrzlr1ixz+/bt5oMPPmhWrlzZTE9PN03TNAcNGmSOHz/+4varVq0yvby8zJdeesncsWOHOWHCBNPb29vcsmWLVXeh3JR0X/3tb38zv/nmG3Pfvn1mcnKyec8995h+fn7mtm3brLoL5SInJ8fcuHGjuXHjRlOSOXnyZHPjxo3mgQMHTNM0zfHjx5uDBg26uP3+/ftNf39/86mnnjJ37NhhJiQkmHa73fz666+tugvlqqT7a8qUKeaCBQvMPXv2mFu2bDHHjBlj2mw2c8mSJVbdhXIxatQoMygoyFy2bJl59OjRi5czZ85c3IbnqyLXsq889fnKNIv+jyUlJZlpaWnm5s2bzfHjx5uGYZjffvutaZrO87giXn7j1VdfNevUqWP6+PiYHTp0MNesWXPxa7GxseaQIUMu2f6jjz4ymzRpYvr4+JgtW7Y0v/zyy3Ke2Dol2Vdjx469uG1ISIh58803mykpKRZMXb5+/lPe315+3jdDhgwxY2Njf/c9kZGRpo+Pj9mgQQNz5syZ5T63VUq6v/7zn/+YDRs2NP38/MyqVaua3bt3N5cuXWrN8OXocvtI0iWPFZ6vilzLvvLU5yvTNM3hw4ebdevWNX18fMzq1aubvXr1uhgupuk8jyvDNE2zbI/tAAAAlB7WvAAAAJdCvAAAAJdCvAAAAJdCvAAAAJdCvAAAAJdCvAAAAJdCvAAAAJdCvAAAAJdCvAAAAJdCvAAAAJdCvAAAAJdCvAAAAJfy/wHwYf0JgTF+vgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(range(len(sigmas)), sigmas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "soviet-defensive",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f81085aa980>]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOT0lEQVR4nO3deVzU1f4/8NcsMMM6bLKDuG8gIAjhknqj1ExN07SszO713hbLsqy8fau7VFrdX9e61/LWTa3boi2mpWYqmkuiKAruKIqA7Ivs+8z5/THM6MgiyDDDB17Px4PHI2Y+M59zuFzm5fu8P+cjE0IIEBEREUmE3NoDICIiImoPhhciIiKSFIYXIiIikhSGFyIiIpIUhhciIiKSFIYXIiIikhSGFyIiIpIUhhciIiKSFKW1B2BuOp0O2dnZcHJygkwms/ZwiIiIqA2EECgvL4evry/k8tZrK90uvGRnZyMgIMDawyAiIqJbkJmZCX9//1aP6XbhxcnJCYB+8s7OzlYeDREREbVFWVkZAgICjJ/jrel24cWwVOTs7MzwQkREJDFtaflgwy4RERFJCsMLERERSQrDCxEREUkKwwsRERFJCsMLERERSQrDCxEREUkKwwsRERFJCsMLERERSQrDCxEREUkKwwsRERFJCsMLERERSQrDCxEREUlKt7sxY2c5crkYW5Kz4ahWwtXeFgO9nBDq7wKNvY21h0ZERNSjMLy00emsUnwWn27ymEIuQ1SQG2aM8MO0UF+obRRWGh0REVHPIRNCCGsPwpzKysqg0WhQWloKZ2dns71vYvpV7DmXj4raBuSV1eBsThkuF1UZn3d3sMUfxvbFgtFBDDFERETt1J7Pb4aXDsgsrsJPJ7LxRXw6sktrAAA+GjWW3DkQM0f4QyGXder5iYiIuguGFwuFF4MGrQ6bk7Lx3s7zyCqpBgAM9XHGK1OGYHR/D4uMgYiISMoYXiwcXgxq6rX47OBl/HtPKsprGgAAvxvsiWWTB2OAl9NNX9+g1SElrxzHMkpwPP0qMoqrIJfJ4Opgg0FeTrh9YC+EB7qyokNERN0Ow4uVwotBcWUdPoi7gC8OpaNBp//xjhvYCzPC/RDVxw1ezmoIIZBTWoMzOWVIzizBsYyrOHGlFFV12lbfu7e7Pf4wti/uj/SHSsneGiIi6h4YXqwcXgwuFVRgxc/nsONMnsnjMhnQ0k/dSaVEWKALwgNcMMhbP/68shokZZZgT0q+saLTt5cD3rg3GKP6cVmKiIikj+Gli4QXg/SiSmw4kom95wtwNqcMjcUY2Crk6OfpiOF+GoQHumBEb1f07+UIeQvLQlV1DfjmSCb+veciCitqAQB/GtcXL9w1CDYK7jdIRETSxfDSxcLL9bQ6geLKOijlMjiplVDeQugora7H29vP4avDGQCAiN6u+HDeCHg5q809XCIiIotoz+c3/7luYQq5DL2cVHB1sL2l4AIAGjsbvDUjBB/NGwEnlRKJ6Vcx5YMDSEgrNvNoiYiIuh6GFwmbHOKDLc+MwWBvJxRW1OLBTw5h3W9puJViWnlNPU5lleLo5WJkFFWhQavrhBETERF1HJeNuoGquga89P1J/JScDQCYGe6HN2eEwM629auRckqr8d3RK9h6MgcpeeUmTcTOaiUmB/tg4e190d/TsTOHT0RExJ6XnhZeAEAIgU8PpGH5z+eg1QkM9nbC3+8NxsggN5Pjahu02H02H98mXsGvKfnG5mEA8HBUwd5WgdyyGtQ16CsvchnwxPh+WHzHQNgqWagjIqLOwfDSA8OLQfzFIiz66hiKKusAAKEBLoju4wa1Uo7UggocvFiEkqp64/FRfdxwf2QAbh/oAU8nfcOvVieQkFaMTw9cwq6z+QCAkUGu+OSRSLjY21p+UkRE1O0xvPTg8AIA+eU1WLnrAjYcyYRW1/R/Xi9nFWaE++P+SH/07dX6ktDPJ3Pw4vcnUF7TgAGejtjwpxi4OTDAEBGReTG89PDwYlBQXotdZ/NwPq8cDVoBb40akb1dERnk1q5bDKTklmP+mgTkltUg1F+DrxbeBgeVshNHTkREPQ3DC8OL2aXmV2D26oO4WlWPScO88dFDIyCT8R5LRERkHtznhcyuv6cjPn10JGwUMmw/nYuP912y9pCIiKiHYnihNhsR6IrXpg4DALy9/RwS07kpHhERWR7DC7XLQ9GBmB7mC50AlnyTjMraBmsPiYiIehiGF2oXmUyGv00Pho9GjfSiKry57ay1h0RERD0MLxmhdtPY2eAfs0Mx77+H8dXhDNw5xAsTBnu2+fVVdQ3YnJSNfecLkF1aA2e1EqP7e+CBqEBo7Gw6ceRERNQd8GojumV//ek01v52Gb2cVNi+eCzcHVU3fU38xSIsXn8c+eW1TZ5zc7DFP+eEYdzAXp0xXCIi6sJ4tRFZxEuTBqO/pyMKymux5Jtk6JrZEM9ApxP4V9wFzPvvIeSX18Lf1Q4v3DUQHz8cgb9NH4b+no4orqzDo2sTsOVEtgVnQUREUsPKC3XI2Zwy3LvqN9Q26LB04iA8NaF/k2OKKmrx7IYk7L9QCACYFeGPv00fBnvba6uWNfVavPLDKXx/7ApsFDKs/+NtiOjt1uS9iIioe2LlhSxmiI8z/jZdf/n0P3ak4NujmSbPH7xYiLs/2I/9FwqhtpHj3VnD8Y/ZoSbBBQDUNgq8M2s4Jgd7o14r8NyGZFTwSiYiImoGG3apw+6PDMDZnHKsO3gZS787gSOXixHdxx2/pRZi4/EsAPpN7j6cNwIDvZxafB+FXIa3Zw3HiSulyCiuwvu7zuOVKUMtNQ0iIpIILhuRWQgh8Na2s/hkf1qT5x6MDsQrdw9p8/2Q9pzLx4J1R2CrkCPu+XEIcLM393CJiKiLac/nNysvZBYymQyvTBmKO4d6Y31CBrJLq9HHwxFzRwYgNMClXe81flAvjO7vjt9Si/BB3AW8Ozu0cwZNRESSxMoLdUnHMq5i5ocHYauQ48DLE+DppLb2kIiIqBOxYZckb0SgK0YEuqBOq8MX8enWHg4REXUhDC/UZT02pg8A4JujV6BtZQ8ZIiLqWRheqMu6c6gXnNVK5JbVIP5ikbWHQ0REXQTDC3VZKqUCU0N9AQDfH7ti5dEQEVFXwfBCXdrMEX4AgJ1n8lDXoLPyaIiIqCtgeKEuLTzAFR6OKlTUNuDQJS4dERERwwt1cXK5DLFDPAEAu87mWXk0RETUFTC8UJd351AvAMCuM3noZtsSERHRLWB4oS5vdH8P2CrkyC6twaXCSmsPh4iIrIzhhbo8tY0CI3q7AAAO8pJpIqIej+GFJGFUPw8AwCGGFyKiHo/hhSRhVD93AED8pSLouNsuEVGPxvBCkhAa4AI7GwWKK+uQWlBh7eEQEZEVMbyQJNgo5BjurwEAJGWUWHcwRERkVQwvJBlhgS4AgOOZJVYdBxERWRfDC0lGeIALAOB4xlXrDoSIiKyK4YUkIzzQFQBwPq8clbUNVh4NERFZC8MLSYaXsxo+GjV0AjiZVWrt4RARkZUwvJCkhBv6Xti0S0TUYzG8kKQM93cBAJzMKmn3a385nYv/7r/E+yMREUlclwwvM2bMgKurK2bNmmXtoVAXM9xPf7l0e5eN9pzLx5/+l4g3tp7F5aKqzhgaERFZSJcML4sXL8bnn39u7WFQFzSsMbxkFlejpKquTa+pqmvA0u9OGL+vrtN2ytiIiMgyumR4GT9+PJycnKw9DOqCNHY26O1uDwA4lVXWptesT8hEYUWt8Xsdl42IiCTN7OFl3759mDp1Knx9fSGTybBp06Ymx6xatQpBQUFQq9WIjo5GQkKCuYdB3VhwY/XlRBv6XoQQWHswrZNHRERElmT28FJZWYnQ0FCsWrWq2ec3bNiAJUuW4PXXX8exY8cQGhqKiRMnIj8/39xDoW7K0Pdyqg19L8czS5BZXA0HWwVc7W0AsPJCRCR1Zg8vkydPxhtvvIEZM2Y0+/x7772HhQsXYsGCBRg6dChWr14Ne3t7rFmz5pbOV1tbi7KyMpMv6t5C2tG0+1NyNgAgdqgX7G2VAABmFyIiabNoz0tdXR0SExMRGxt7bQByOWJjYxEfH39L77l8+XJoNBrjV0BAgLmGS11Ue5p295zTV/QmB/tAJtM/xuxCRCRtFg0vhYWF0Gq18PLyMnncy8sLubm5xu9jY2Mxe/ZsbNu2Df7+/q0Gm2XLlqG0tNT4lZmZ2Wnjp66hrU27V65W4XJRFRRyGUb1dzeGFy4bERFJm9LaA2jOrl272nysSqWCSqXqxNFQVxTip0F6URVOZJVgzACPZo85mFoEAAj118BZbQMZ9OmF2YWISNosWnnx8PCAQqFAXl6eyeN5eXnw9va25FBI4kLa0LR7OK0YADCqnz7cyA3LRkwvRESSZtHwYmtri4iICMTFxRkf0+l0iIuLQ0xMjCWHQhLXlqbdE1dKAABhAS4AAFnjuhGjCxGRtJl92aiiogKpqanG79PS0pCUlAQ3NzcEBgZiyZIlmD9/PiIjIxEVFYWVK1eisrISCxYsMPdQqBu7sWnXxd7W5PmK2gakFlQAAIYH6I81NuwyvRARSZrZw8vRo0cxYcIE4/dLliwBAMyfPx/r1q3DnDlzUFBQgNdeew25ubkICwvD9u3bmzTxErVGY2eDIHd7XC6qwsmsUowd0Mvk+VNZpRAC8NWo4emkBoDGjhc27BIRSZ3Zw8v48eNv2lOwaNEiLFq0yNynph4m2E/TYngxLBkZ7kINXLdsxOxCRCRpXfLeRkRt0VrT7okr+scMS0YAG3aJiLoLhheSrNaads9k6/d/Cfa9Fl6Ml0pbYGxERNR5GF5Isq5v2r1aeW2n3aq6BqQVVQIAhvg4Gx9nwy4RUffA8EKSpbGzQX9PRwDA/tRC4+MpueUQAvBwVKGX07UNDA09L2zYJSKSNoYXkrS7huqvUvv5ZI7xsbM55QCAIT5OJscarjZidCEikjaGF5K0u0N8AAB7UvJRXlMPADibo+93GXrdkhEAyBt/21l5ISKSNoYXkrRhvs7o18sBNfU6/JSsr74YwsuQG8KLoWGXpRciImljeCFJk8lkmDsyEACw/kgGdDqBc7mGZaMbKi/G7ML0QkQkZQwvJHn3RfjDViHHiSul2HIyBxW1DbBVytG3l4PpgYaGXZ0VBklERGbD8EKS5+Zgi7uG6Rt3n/n6OAAgItAVNgrTX2827BIRdQ8ML9QtPHRbb5PvR/Vzb3IMd9glIuoeGF6oW4ju42YMLEq5DFNDfZscc22fF4sOjYiIzMzsN2YksgaZTIZ/PzgCnx28jBG9XRHk4dDkGEPlhQtHRETSxvBC3Yabgy2eu3Ngi88bLpVm5YWISNq4bEQ9B+9tRETULTC8UI/BfV6IiLoHhhfqMbhsRETUPTC8UI9huLcRL5UmIpI2hhfqMQyVF2YXIiJpY3ihHkPGnhciom6B4YV6DMMmday8EBFJG8ML9RiGPerYsEtEJG0ML9Rj8N5GRETdA8ML9RhcNiIi6h4YXqjHMCwbsWG362N1jIhaw3sbUY/RGZUXIQROZpWitkGHsAAX2Cja9u+BE1dK8EFcKgrKa/C7wV54fHxfqJQK8w3sBkmZJXj753NIK6zEbX3d8Oe7h8DTWd1p5wOAugYdckqr4aOxg62ybT+Xy4WV+PMPJ3HoUhH6eDjgtanDMG5gr04dJxFJD8ML9RiGS6Wba9iNO5uHj369iIraBkwP88PCsX2gvEkQKa+pxx8/T0T8pSIAQH9PR6yZPxKB7vatvu7gxUIsWHsEtQ06AEDylVIkZlzF2kdHQnHt1tfNKqqoxWcHLyO3rAZ3DvXGnUO9Wj0eAI5cLsYjnyagul4LANiUlI1jGSXY/NRouDrYtvratMJKvLn1DM7mlCMyyBV/vnsIvNoQerafysGyjSdxtaoe7g62eGtmCCYO8271NfllNXjwk0PILq0BAFwsqMRj647gi99HI6af+03PSUQ9B5eNqMdo6d5G3ydewR8+P4qj6VdxLrccb28/h5e+P9nq0oUQAk9+eQzxl4qgUsrhpFIiNb8Cj65LQE1jSGhOaVU9nvn6OGobdBg3sBfenBEMe1sF9p0vwEe/prY6/qySakz792/4YHcqvjl6BQs/P4p3fznX6msqaxvw7PokVNdrMXaABz57LAoBbnbIKK7C/2061eprM4qqcP9/4rHrbD6ySqqxOSkbs1fH42plXauvO3ChEIu+Oo6rVfUAgKLKOjz55TEkpBW3+rr/23QK2aU16NvLAVufGYMpIT7Q6gQWrz+OytqGVl9LRD0Lwwv1GM3d2yijSP8hLgQwJzIAf502DAq5DN8fu4KtJ3NafK9NSVnYf6EQahs5vnt8FOKeHwdPJxUuFVTig7gLLb7un7vOo7CiDv09HfGfhyMwL7o3/jY9GACwas9FFJTXNvs6nU5gyYYkZJVUo4+HAx6ICjS+Zu/5ghbP98n+S8gqqYa/qx1WPxSBcQN74aN5EVDKZdh6MgeJ6c0HCiEEXvr+BArKazHY2wmfPBJpDD1/33KmxfNV12nx4nfJaNAJ3Bvmi3N/n4R7hutDyHMbkloMdglpxdhxJg8KuQyrH4rAMF8N3p09HIFu9sgvr8Xa39JaPCcR9TwML9RjyIwdu9fSy9u/nEN1vRa39XXD8pkhmD8qCIsm9AcAvLHlLOoal3aup9UJrNylDyhP/24AQvw18HRWG0PI2t8uN1udKKqoxfojGQCA16cOhdpG3+Ny3wg/hAW4oLpei4/3XWx27FtP5uBwWjHsbRX4bEEUls8MwYLRQQCA1zafQoO26Tgraxuw9rfLAICXJg2Gg0q/Shzsp8GsCH8AwAdxzVd7dp7JM1aVPnkkEncO9cL7c8MhkwEbj2fhVFZps69bezAN2aU18HOxw4r7hkNto8Db9w2Hj0aNrJJqfHk4o9nXfbzvEgDg/sgADPRyAgDY2yrx/F0DAQD/PZDWakWLiHoWhhfqMeSGht3G7y8WVGBbY3XlL9OGQd64rvTkhH7wclYht6wGm5KymrzPrrN5SC+qgsbOxhggAGDiMC8M83VGdb0W/zuU3uR1649koqZeh+H+Gozp72F8XCaT4enf6QPTt4lXmnxICyGwao8+ZPzp9n7GnpoX7hoENwdbpBdVNVsl2nIiG6XV9Qhyt8fdIT4mzz05Xn++vecLcOVqVZPXfrJfHyYeG9MHAW76840IdMW0UF8AwOq9TUNWXYMOnx28DABYcudAYzhzUCnxzB0DAAD/3X+pSdC6crUKcefyAAC/H9PH5Ll7hvvCV6NGSVU9fj7VciWMiHoWhhfqOQwNu43rRl8eyoAQQOwQTwz2djYeplIq8Nho/Yeo4cP4ehuOZAIAHogKhL3ttZ53mUyGP4zVv+67xCsmPTNCCPxwXB+EHorubbzyyWD8IE/jh/Qvp3NNnjP04jjYKvDoqCDj4w4qJebH6L83VFiuZzjf/SMDmjQCB7rbY1RjE+z3iaYB7Ux2GY5cvgobhczkfIA+PAHAz6dymyxx7TiTi7yyWvRyUmFqY8gxmDnCD24OtsgprcGeFNNlrh+TsyEEMKqfO/p7Opo8p5DLMLdxiezbo1eazJGIeiaGF+oxrq+81DXojFWVB6MDmxw7OzIASrkMp7PLcD6v3Ph4cWUd9jX2mBiWXq43cZg3HGwVyCiuwtH0q8bHz+aUIzW/ArZKOSaFNL3qRiGXYVZkAABg4zHTMLGpMYRMDvGBxt7G5LkHowOhkMuQlFmCiwUVxsezSqpx6JK+n+XeML9mfx6zI/Xj/+5YpknQ2tz4c7ljsFeTK4uG+jojNMAFWp3AT8nZN4xT//3sCP8ml0arlArMbvx5fZ1gunS0JVlfUZl2Q+AxmBGuH/+hS0Uoqmi+J4iIehaGF+oxrr+30Z6UfBRX1qGXkwq3D2i6j4ibgy3GD/IEcK2CAQDbTuagQScwzNe5SZUA0PdpGJZovk+8Vim4Fgg84ay2afI6ALg3TP/hfSC1EIWNH9J1DTpsOaH/cDd8iF9PP379EtQP14Uew/lu6+sGXxe7Zs83aZgPHGwVyCyuxvHMEgD6qtSPjaHk3vDmw8TMxnFc/3MprarH3vP5AIDpLYYlfTjbe74A+eX6y6EvFVTgTE4ZlHJZi5dSB7jZI8RPA50AdpzJa/YYIupZGF6ox5Bdd28jQ7CYGe7X4n4uM0foP4Q3H88yLjUZPthbqhLoX6evMGw9mYOaeq1JIJge1vLr+vZyxHB/DbQ6YezF+TUlH6XV9fB0UuG2vs3vdWI43w/XjzOpMYC0ECQAwM5WYdwnxlBFOXK5GDmlNXBSK43h7UZTQ32hlMtwMqsUqfn6qtT20zmo1woM9nbCIG+nZl/X39MR4YH6qo2hmmQIZqP7e7S658zkxmrVtlauACOinoPhhXoMw7LR1ao67EnRVwkMH/zN+d1gTziplMgurcHR9KvILa3Bkcv6pZgbezquF93HDT4aNcprGvBrSj4OpzUGAlXLgcDAEIo2N4YPQ3VjephvixvY3TnUC44qJbJKqnE0/SpS88txLrccNgoZJgf7NPsaA8M8tp7I0YeKxorN5GBvY8PtjfRVKX21yrBUZBhvaz8XAJgdoa++GHqCDBWie4a3Pk7DPA5eLLrpPjNE1P0xvFCPYfjo/zE5G/VagWA/5xarBACgtlFgYrD+X/ybk7Kw9WQOhAAie7u2uBQDAHK5DNMaKyw/HM/CxmP6Ks+U4T4tBgKDaaG+kMmAxPSrOHmlFHFn9SFrRnjLIUtto8DkxnH+cDwLPzX2kIwd0KtJj8yNxg7oBWe1Evnltdh3ocDYf9JaxQYA7m1cOtqUlIWskmrjLsOtVaQA4J5QH6iUcpzPq8D6I5m4WFAJW6Xc+HNuSR8PBwzxcYZWJ/DzqdxWjyWi7o/hhXoMwxU+mcXVAICZrQQCA8Myz9aTOfjqsP7y55tVCYBrH/6/nM7Dt4YlqlaqPAaezmrjVUAL1h1BnVaHYb7OGOrr3OrrDP0wm5Oy8Hn85TaP01YpN1Y1Fqw9gvLaBvi52LW4RGUQO0Rf7blytRr/98NJCAFEBbkZL6tuibPaBpMag8qyjScBtN4HdD1DMLqxUZiIeh6GF+oxbuvrZvxvD0cV5kYF3PQ1MX3dEeBmh5KqelwsqITaRn7TpREAGOLjjKiga+cb5uuMyN6ubRqnoeHV0LQ7d+TNxxnd1x19ezmgqk5rvJ/QzZaMDG78OTwYHWjc86YlahsF7m7sQzFc+my4eulmDEtHBo/dsLdLSwz9OYkZV29yJBF1dwwv1GPMivDH4+P64ba+bvjooREme7S0RKmQY+nEwcbvF47tC3dHVZvO9+aMYPi52MHLWYUVM4ffNBAYTA/zRViACwAg2M8Zc0Y2vZT7Rgq5DK/dM9TYF/PipEGws23bXarDA13x0G36cwz1cW6yUVxLnv7dANg3nmOgl2OzV0M1Z3R/dzwQFQCVUo4/3t4XI68Lea1xbNwhWNfcnTWJqEeRidbuPidBZWVl0Gg0KC0thbNz66V2oraKv1iE8pp6xA7xanMIAWDcP+XGTelupqZei+MZJRjurzFu698W6UWVaNAJ9OvV9DLum8kprYank/qmd7a+3pnsMhxOK8I9w33Ry6ltoc5ACNGun0teWQ2i34qDQi7Dxbfubte5iKjra8/nd9v/KhL1YDH9Wu8BaUl7Q4uB2kZxS+fs7e5wS+cDAB9Ny03ILRnahn6clrT3Z3Pt1lTd6t9bRHQLuGxERNJg2KfHuqMgoi6A4YWIiIgkheGFiCRB1lh64aoRETG8EJEk3GL7EBF1QwwvRCQJ12cXNu0S9WwML0QkCddfncTsQtSzMbwQkSSYVF6sNgoi6goYXohIEq7veeGyEVHPxvBCRJLD6ELUszG8EJEkyMDLjYhIj+GFiKTBZNnIesMgIutjeCEiSTDpeeHCEVGPxvBCRJJgus+L1YZBRF0AwwsRScKt3qGbiLofhhcikgRWXojIgOGFiCSBPS9EZMDwQkSSw8oLUc/G8EJEksB9XojIgOGFiCTBdNmIiHoyhhcikhze24ioZ2N4ISJJYOWFiAwYXohIEq7veWHhhahnY3ghIkkw2aOO4YWoR2N4ISLJ4T4vRD0bwwsRSQJ32CUiA4YXIpIE3tuIiAwYXohIEtjyQkQGXTK8zJgxA66urpg1a5a1h0JEXYTJpdJcNyLq0bpkeFm8eDE+//xzaw+DiLqQ65eNGF2IerYuGV7Gjx8PJycnaw+DiLooFl6IerZ2h5d9+/Zh6tSp8PX1hUwmw6ZNm5ocs2rVKgQFBUGtViM6OhoJCQnmGCsREQBeKk3U07U7vFRWViI0NBSrVq1q9vkNGzZgyZIleP3113Hs2DGEhoZi4sSJyM/PNx4TFhaG4ODgJl/Z2dm3PhMi6vaMK0fMLkQ9mrK9L5g8eTImT57c4vPvvfceFi5ciAULFgAAVq9eja1bt2LNmjV4+eWXAQBJSUm3Ntpm1NbWora21vh9WVmZ2d6biLoWGfS5hdmFqGcza89LXV0dEhMTERsbe+0EcjliY2MRHx9vzlMZLV++HBqNxvgVEBDQKechIuvjXi9EBJg5vBQWFkKr1cLLy8vkcS8vL+Tm5rb5fWJjYzF79mxs27YN/v7+rQafZcuWobS01PiVmZl5y+Mnoq7NuGrE0gtRj9buZSNL2LVrV5uPValUUKlUnTgaIuoqDIUXNuwS9Wxmrbx4eHhAoVAgLy/P5PG8vDx4e3ub81RE1APJGmsvrLwQ9WxmDS+2traIiIhAXFyc8TGdToe4uDjExMSY81RE1BMZKy9E1JO1e9mooqICqampxu/T0tKQlJQENzc3BAYGYsmSJZg/fz4iIyMRFRWFlStXorKy0nj1ERFRR/H2AEQ9W7vDy9GjRzFhwgTj90uWLAEAzJ8/H+vWrcOcOXNQUFCA1157Dbm5uQgLC8P27dubNPESEbUXG3aJCLiF8DJ+/Pib/qtn0aJFWLRo0S0PioioObxSmoiALnpvIyKi5rBhl4gAhhcikhBWXogIYHghIgm5dmsjll6IejKGFyKSDMPtAbhsRNSzMbwQkeQwuxD1bAwvRCQZ1y6VZnwh6skYXohIOrjDLhGB4YWIJISb1BERwPBCRBIik1273oiIei6GFyKSDO7zQkQAwwsRSQiXjYgIYHghIgkx7vNi5XEQkXUxvBCR5LDyQtSzMbwQkWTw9gBEBDC8EJGEGBp2WXkh6tkYXohIQnhvIyJieCEiCTFWXrhsRNSjMbwQkWRwmxciAhheiEhC2PNCRADDCxEREUkMwwsRSYaMDbtEBIYXIpIQNuwSEcDwQkQS0pF7G9U2aHGpoALVdVqzjomILE9p7QEQEbXVrd7baOOxK/j7ljO4WlUPB1sFnrtzIH4/po/x/YhIWhheiEhyRDtKL98czcSL350AACjlMlTWafHG1rPQCYE/3t6vU8ZXWl2Pf+48j73nC+DppMKi3/XH2AG9OuVcRD0Rl42ISDLaWyhJL6rEq5tOAQB+P6YPzv59El6aNBgAsOLnczidXWruIaK0uh73r47HuoOXkVZYicNpxXhkTQJ+Ss42+7kMhBBYn5CBuR/HY8HaBOxJye+0cxF1BQwvRCQZ1xp22+aNrWdR26DDqH7u+L8pQ2CjkOPxcX0xJcQHOgH87aczZh/jnzeeREpeOXo5qfDhvBGYFeEPIYCl3yUjrbDS7OcTQuDVzafw8saTOHSpGHtSCrBg7RH871C62c91/TnXHEjD2Hd2I+rNXVj+81nUNrCXiCyH4YWIJKctq0YpueXYeSYPMhnwt+nDjP0tMpkM/3fPENgq5DicVoxjGVfNNq4DFwqx9WQOFHIZPp0fibtDfPD2fcMxur87aup1eGvbWbOdy+C7xCv44lAGZDLghbsG4sHoQADA3346jTPZZWY/HwD8Y0cK/rblDDKLq5FfXov/7L2Ep748Dp2OV4GRZTC8EJFkyK5db3TTY9cdTAMATBrmjf6eTibP+WjscG+4LwBgzYE0s43vg90XAAAPRQdiuL8LAEAhl+Gv04ZBIZdh55k8swaKspp6rPj5HADghbsGYdHvBuDNe4MxcZgX6rUCf/nxtNnOZfBbaiFW7bkIAPjz3YPxrwfCoVLKsetsHr4+kmH28xE1h+GFiCSjrbcHqKnXYsuJHADAwzG9mz3mkZggAMCOM3koq6nv8NhOXilFQloxlHIZnhjf3+S5/p5OmBTsDQD47ODlDp/LYH1CBooq69C3lwP+eHtfAPrK0l+nBcNWIUfC5WIculRktvNpdQJ/36JfanvotkD88fZ+mBrqixcb+4je2Z6CytoGs52PqCUML0QkGW2tu8SdzUd5TQN8NWrc1se92WOG+TpjoJcj6hp0+PlkTofH9tMJfUPupGBveGvUTZ5/dFQQAODH5GxU1XX8A16rE8a+lj/d3hc2imt/zr01atw/0h8AsO63y62+T71WhwMXChF/sQjamyz77DmXj3O55XBSK/HCXYOMjz86Kgh9PBxQWl2PrxNarr7U1Gvx0a8XceJKyU1mR9Q6hhcikgzjPi83SS+GK3vuDfeDXN78JUoymQzTw/wAAL+czuvQuIQQ2HlG/x6Tg32aPSaytysC3OxQXa/FnnMFHTofAOw9n4/M4mpo7GwwLdSvyfMP3aavOMWdy0NRRW2z75FfVoN7PjiAhz49jAc+OYSZH/7W4rEA8Fn8ZQDAg1GBcLG3NT6ukMuMlZ8vDqW3eCn73I8P4e3t5/DGFvP3/lDPwvBCRJJxbYfdltNLvVaH31ILAQATh3m3+n6xQ7wA6Ps4OrLz7sWCCqQVVsJWIce4Qc3v5yKTyTAlRN9ns/Vk85dN12t1eG9HCsa+sxvTV/2GAxcKWzynYVlsRrgf7GwVTZ4f7O2MED8N6rUCW5upLDVodVj4+VGk5OkrKY4qJZKvlOLpr483+/PNK6vB/sbxGILR9aaF+sLORoHLRVVIyixp8vyBC4XGxxMuF7c4L6K2YHghIulowz4vxzNKUF7bAFd7GwT7aVo9dqCXI/xc7FDboMPBiy0HhZvZ0Vh1iennDkdVy3t/TgnRV2V2n8tvdunozxtP4oPdqcgsrkZyZgkeXZtgDGLXa9DqEHdWv5fL5OCWA9rUUP35djRTWfo6IQPJV0qhsbPBT4vG4IcnR0FtI8fBi0X4LvFKk+N/OZ0LABgR6IIAN/smzzuolJg4TB8GNx3PMnlOCIF3d6QYvw8LcGlxzERtwfBCRJLT2qrRvvP6JZkxA3pB0cKSkYFMJsMdQzwBALvONt3YrbpOi7/8eBpj39mNR9cm4GJBRbPvYwgHdzV+eLck2M8Z/q52qKnXYd9501Cy51w+vk28AoVchjdnBGNKiA8adAIvfneiSRNswuVilFbXw83BFhG9XVs8351D9cHm0KUilFZfa0qubdDi/bhUAPrLq4M8HDDAywmL7xgIAFi1J7VJ/8u2xurN3SHNL4sB+mU6QF8VatDqjI8nZZYguZlqDNGtYnghIsloy40Z913Qh5fbB3i06T0nDNaHl70p+SbLJUIIPPFlItYdvIzM4mr8mlKA2avjceVqlcnr88tqjMshhmWoFscvk+GuxkCx40yu8XGdTuDNxj1gHhsdhHnRvfHOrOHwc7FDVkk1Po833XDOEJbuGOwJpaLlP+N9PBzQ39MRDTqBX6/bdXf7qVwUVtTCy1mFuVGBxscfiekNF3sbXC6qMlZaAKCoohYJafqlntaW4kb394CLvQ2KKutMloYMTbxOan1VirvBUEcxvBCRZBgadqvqGprtyyiurMPJLP2W/+MGtu1eQjF93aFSypFdWoPzedcqK+uPZOLXlAKobeR4d9ZwDPVxRnFlHV7bbLp3iqFiExrgAi/nplcZ3chQndl9Lt9YnTiQWojU/Ao4qpR4+o4BAPTLMEvu1FdCPj2Qhpp6fU/O9c3Bd92kpwcA7hqqP9/1S0f/awxDD0b1NrlKyUGlxMON/SxfXLdD744zedAJYLi/ptklIwMbhRwTG8OZoVJTUduAn5JzGs8X2OJridqD4YWIJMNQ9fj9Z0dx9wcHcCGv3OT5/RcKIAQw2NsJnm0IEgCgtlFgVD/95dSGewLV1Gvx/3acB6Df/G12ZAD+/WA4lHIZdp/LN+lDMVRQDCHhZiJ7u8LV3gYlVfU4clm/u+/a3/Qb5c2K8Iez2sZ47LQwX/i52KGwohY/JumbfE9nlyGrpBp2NgqMbUN1yRBwfk3JR029Fim55TiafhVKuQwPRAU0OX5uVCDkMuDgxSLjMpkhiExqpb/GYHKI/pjtp/Kg1QlsO5GD6not+vZywMggN/1B7bixJlFzGF6ISDL69XI0/vfZnDLM++9hlFTVGR8z9JG0tepiMH6Qfulozzl9ePk28QoKK2rh52KH+Y37s/Tt5Wi8yuafO89DCIGK2gYcTNVvAtfW8KJUyHFH4/LSjjO5SCusxJ6UAshk1/aCMbBRyPFI4yZ7a35LgxACP5/SB4nbB3pAbdP0KqMbDffTwNtZjco6LQ5eLMSGI5kA9EtczQU8Pxc7/K5xKe3LQxkoqqjFwYv6ObZ0Gfj1Rvf3gMbOBoUVtThyudjY/Dsrwr/dN9YkagnDCxFJxgt3DcKC0UHY/uxY9O3lgPzyWry5Vd8rotOJa/0u7QwvExrDy9H0q7haWYeP9+m3v//jDZu/PTG+H2yVchxNv4r4S0XYm1KAOq0OQe726O/p2Ox7N+f6pRxD1WXCIE8EeTg0OXbuyEDY2ShwLrcc8ZeKsO2kvtLTWuPs9eRymXGp6vvELGw8rg8Tc5qpuhjMi9YHpu8SM7HxWBa0OoEQPw36NDO+G9ko5LizcX4f/noRCZeLIZPpL+k2YN2FOorhhYgkY8JgT7w+dRgGezvjnfuGA9BXSU5llSL5SgkKymvhqFIiMqjlK3CaE+huj769HKDVCSzekITM4mq4O9ji/kjTD3gvZzXmjtQ/9kHcBWxp3FV34jBvYz9OW4wd0AtqG7lJM+7vx/Rp9liNvQ3ui9B/8D+3IUm/n4zyWvWmLQzBYevJHJRU1cNXo8btA1oOeLcP7AV/VzuU1TQYG4mnhfq2+XyGS8INV35NGOQJH40dKy9kNgwvRCRJkUFuxg/Ut7efMzaxjhvUCyrlzZdTbmR4L8MH7mNj+jS7+dvj4/rBRiHDoUvF+PmUvgoytR0f7ABgZ6vAPcOvvSY80MXYd9OcR0fpg01emX732ykhPq3uJ3Oj8EBXxPS99v7P3zWo1cvIFXKZyUZ0Tmol7ovwb/P5bh/YCyMCXQAASrnM2HhswJYX6qi2//YTEXUxL9w1CD+fysH+C4XG3V/b2ntyo0dHBeGLQxkorKjFIC8nPDa6+UqIr4sdZkcG4KvD+st/Y/q6Y5ivc7vP9+qUoSipqkNlrRZv3ze81cpNf09HzIkMwIajmXBSK/Fc7MAWj23Jvx4Mx6cH0tDbzR4zRzS9ncCNHhvdBxfzK3DwYhFemTIEbg62N32NgUIuw3/nj8T3iVcQHuhi3CxQ1pZdBonaQCZa22dbgsrKyqDRaFBaWgpn5/b/QSEiafnLj6exrvFOzb4aNfYsHX9LlRcAyC2tweG0Iowf6AmNvU2LxxVV1OL5b5NRVafFu7OGo7f7zXtBOkqnE0i+UoIAN3t4OKo6/XydYc+5fCxYdwTBfs7Y8vRYaw+Hupj2fH6z8kJEkvbipEGorG3AqewyvHbP0FsOLoD+bsyGmzW2xt1RhXULom75PLdCLpchPLB9vTxdDgsvZCYML0Qkafa2Srw7O9Taw6B26F71frIGNuwSEZFFsPBC5sLwQkREFsXKC3UUwwsREVlEe/bCIWoNwwsREVkUCy/UUQwvRERkEay7kLkwvBARkUV1s+3FyAoYXoiIyCLY8kLmwvBCREREksLwQkREFsF7G5G5MLwQEZFFseWFOorhhYiILII9L2QuDC9ERGRRgju9UAcxvBARkUWw8ELmwvBCREQWxZ4X6iiGFyIisgyWXshMGF6IiMiiWHihjmJ4ISIii+A+L2QuDC9ERGRRvLcRdRTDCxERWQT3eSFzYXghIiKLYt2FOorhhYiILIKFFzIXhhciIrIsll6ogxheiIjIImRseiEz6XLhpaSkBJGRkQgLC0NwcDA++eQTaw+JiIjMiIUX6iiltQdwIycnJ+zbtw/29vaorKxEcHAwZs6cCXd3d2sPjYiIOoCFFzKXLld5USgUsLe3BwDU1tZCCME9AYiIuhH+TaeOand42bdvH6ZOnQpfX1/IZDJs2rSpyTGrVq1CUFAQ1Go1oqOjkZCQ0K5zlJSUIDQ0FP7+/li6dCk8PDzaO0wiIupiWHghc2l3eKmsrERoaChWrVrV7PMbNmzAkiVL8Prrr+PYsWMIDQ3FxIkTkZ+fbzzG0M9y41d2djYAwMXFBcnJyUhLS8NXX32FvLy8W5weERF1Nay7UEe1u+dl8uTJmDx5covPv/fee1i4cCEWLFgAAFi9ejW2bt2KNWvW4OWXXwYAJCUltelcXl5eCA0Nxf79+zFr1qxmj6mtrUVtba3x+7KysjbOhIiILIk9L2QuZu15qaurQ2JiImJjY6+dQC5HbGws4uPj2/QeeXl5KC8vBwCUlpZi3759GDRoUIvHL1++HBqNxvgVEBDQsUkQEVGnYssLdZRZw0thYSG0Wi28vLxMHvfy8kJubm6b3iM9PR1jx45FaGgoxo4di6effhohISEtHr9s2TKUlpYavzIzMzs0ByIi6iwsvZB5dLlLpaOiotq8rAQAKpUKKpWq8wZERERmJdj1Qh1k1sqLh4cHFApFkwbbvLw8eHt7m/NUREQkMex5IXMxa3ixtbVFREQE4uLijI/pdDrExcUhJibGnKciIiKJYs8LdVS7l40qKiqQmppq/D4tLQ1JSUlwc3NDYGAglixZgvnz5yMyMhJRUVFYuXIlKisrjVcfERFRz8TCC5lLu8PL0aNHMWHCBOP3S5YsAQDMnz8f69atw5w5c1BQUIDXXnsNubm5CAsLw/bt25s08RIRUc/Eygt1VLvDy/jx42+6tfOiRYuwaNGiWx4UERF1P7yrNJlLl7u3EREREVFrGF6IiMgiWHchc2F4ISIii+JdpamjGF6IiMgi2PJC5sLwQkREFsW6C3UUwwsREVmEjF0vZCYML0REZFFseaGOYnghIiKLYM8LmQvDCxERWRTvKk0dxfBCREREksLwQkREFsWeF+oohhciIrII9ryQuTC8EBGRRbHwQh3F8EJERBbBfV7IXBheiIjIotjzQh3F8EJERBbBnhcyF4YXIiKyMJZeqGMYXoiIyCJYeSFzYXghIiKLYs8LdRTDCxERWQSvNiJzYXghIiKLYuGFOorhhYiILII9L2QuDC9ERGRRgk0v1EEML0REZBEsvJC5MLwQEZFFse5CHcXwQkREFsGeFzIXhhciIrIotrxQRzG8EBGRhbD0QubB8EJERBbFq42ooxheiIjIItjzQubC8EJERBbFugt1FMMLERFZBAsvZC4ML0REZFksvVAHMbwQEZFFyNj0QmbC8EJERBbFwgt1FMMLERFZBOsuZC4ML0REZFHc54U6iuGFiIgsgi0vZC4ML0REZFGsu1BHMbwQEZFFyNj1QmbC8EJERBbFlhfqKIYXIiKyCPa8kLkwvBARkUUJdr1QBzG8EBERkaQwvBARkUWx54U6iuGFiIgsgj0vZC4ML0REZFEsvFBHMbwQEZFF8K7SZC4ML0REZFksvVAHMbwQEZFFsO5C5sLwQkREFsV9XqijGF6IiMgi2PJC5sLwQkREFsV9XqijGF6IiMgieFdpMheGFyIisigWXqijGF6IiMgi2PNC5sLwQkREFiXY9EIdxPBCREQWwcILmQvDCxERWRTrLtRRDC9ERGQZLL2QmTC8EBGRRbHlhTqK4YWIiCyC+7yQuTC8EBERkaQwvBARkUVwnxcyF4YXIiKyOO71Qh3B8EJERBbBwguZC8MLERFZHAsv1BEML0REZBEyNr2QmTC8EBGRxbHwQh2htPYAmhMUFARnZ2fI5XK4urpiz5491h4SERF1EOsuZC5dMrwAwMGDB+Ho6GjtYRARUSfQX23EOEO3hstGRERkEWx5IXNpd3jZt28fpk6dCl9fX8hkMmzatKnJMatWrUJQUBDUajWio6ORkJDQrnPIZDKMGzcOI0eOxJdfftneIRIRURfHnhfqiHYvG1VWViI0NBSPPfYYZs6c2eT5DRs2YMmSJVi9ejWio6OxcuVKTJw4ESkpKfD09AQAhIWFoaGhoclrd+zYAV9fXxw4cAB+fn7IyclBbGwsQkJCMHz48FuYHhERdRW8txGZS7vDy+TJkzF58uQWn3/vvfewcOFCLFiwAACwevVqbN26FWvWrMHLL78MAEhKSmr1HH5+fgAAHx8f3H333Th27FiL4aW2tha1tbXG78vKytozHSIisgLu80IdYdael7q6OiQmJiI2NvbaCeRyxMbGIj4+vk3vUVlZifLycgBARUUFdu/ejWHDhrV4/PLly6HRaIxfAQEBHZsEERF1DhZeyEzMGl4KCwuh1Wrh5eVl8riXlxdyc3Pb9B55eXkYM2YMQkNDcdttt+GRRx7ByJEjWzx+2bJlKC0tNX5lZmZ2aA5ERNT5BLteqAO63KXSffv2RXJycpuPV6lUUKlUnTgiIiIyB15tROZi1sqLh4cHFAoF8vLyTB7Py8uDt7e3OU9FREQSxp4X6gizhhdbW1tEREQgLi7O+JhOp0NcXBxiYmLMeSoiIpIYFl7IXNq9bFRRUYHU1FTj92lpaUhKSoKbmxsCAwOxZMkSzJ8/H5GRkYiKisLKlStRWVlpvPqIiIiIqCPaHV6OHj2KCRMmGL9fsmQJAGD+/PlYt24d5syZg4KCArz22mvIzc1FWFgYtm/f3qSJl4iIehbeVZrMpd3hZfz48Y33pGjZokWLsGjRolseFBERdW/seaGO4L2NiIjIIlh3IXNheCEiIovjPi/UEQwvRERkEWx5MT8hBKrrtDdt5+huutwmdURE1P11t8/aBq0Ou8/lY09KAfLKauBib4Pb+rjj7uE+cFSZ/6M2Mf0qPj1wCfsvFKK8pgH2tgqM6e+BP4zti6g+bmY/n4EQAmU1DVDKZXDohHm1FcMLERFZRHe9q/SxjKtY9v1JpOSVmzy+8VgWVmw/h1fvGYIZ4f5mOVddgw5vbD2Dz+PTTR6vqtNix5k87DiThzmRAfjr9GFQ2yjMck4ASM4swX/2XTSGpUdHBeEv01q+72BnY3ghIiKL6+zCS2lVPf574BK2n8pFVkk17G2ViOrjivkxQYju62628/xw/Ape/O4E6rUCLvY2mBnujwFejsgprcHmpCykF1XhuQ3JOJ5RgtenDoNCfusBrrpOi4WfH8WB1EIAwKwIfzwS0xtBHg7IKKrCl4czsOFIBjYczcTFggp89lhUh6sjWp3A/9uRgo/2XjSpljXodB16346SiW62UFZWVgaNRoPS0lI4OztbezhERNSopl6Lwa9uBwCc+uvETllOAYC95wvw7PrjuFpV3+zzM0f44W/Tgzt8/q0ncrDo62MQApg0zBsr7guBi72t8fl6rQ4f7rmIlXHnIQQwd2QAls8MuaX9bhq0Ovzpf4mIO5cPB1sF/v3gCEwY7NnkuAMXCvHkl4koq2nA6P7uWPPoSKiUt1aBadDqsOSbZPyYnA0AmB7mi8dG98EgbyeolHKz79vTns9vNuwSEZHFdda/mzcdz8KCtQm4WlWPAZ6OeH9uGPYuHY/vnxiFB6ICIZfpl3Me/OQQiivrbvk8yZkleHbDcQgBPBAViA/njTAJLgBgo5BjcewAfDA3HHIZsP5IJv7fjvO3dL43t51F3Ll8qJRyfPZYVLPBBQDGDPDQV1xsFfgttQh/33Lmls4HAH/fcgY/JmdDKZfh/blheH9uOEIDXKC2UVh9w0GGFyIisridZ/JQVFFr1vfcf6EAL3ybDJ3QV1e2PDMG08P80NvdARG9XbF8Zgg2/CkGbg62OHGlFI+uTUBVXUO7z1NTr8Xz3yajXitw11AvvHFvMOStLAdNDfXFipnDAQD/3pOKX07ntut8v6bkY+1vlwEA788NR2RQ6w254YGu+Pe8EZDJgC8OZWDT8ax2nQ8Avjycjs/i0yGTAf9+MBzTw/za/R6dieGFiIgsQn7dv9aXfJOMqLfisGzjCVztQAXEIL+sBovXJ6FBJ3BvmC/+MSu02eWSkUFu+OZPMXC1t8GJK6V45usk6HTtqwK9t/M8UvMr0MtJhbfvG96mPpb7RwbgsdF9AAAvfJOMtMLKNp2rqKIWL3x7AgDw6KggTAr2btPrJgzyxNO/GwAAeHXzKeSW1rTpdQCQVlhprNgsnTgIk4J92vxaS2F4ISIii7BVyrH4jgGI6uOGwd5O0OoEvk7IxNR/H0BqfvnN36AFQgg8/20yiivrMNTHGW/PGt5qJaS/pyP+Oz8Stko5dp3Nw38PXGrzuY5cLsYn+/XHr5gZAlcH25u84ppldw9GVJAbymsb8NSXx1BTr231eCEEXvr+BAorajHQyxEvTx7c5nMBwOI7BiA0wAXlNQ348w8n27RUp9UJLP02GTX1Oozp74EnxvVr1zktheGFiIgs5rk7B+KbP8Vg+7O349vHY9Db3R5XrlZj7seHkVFUdUvv+cvpPOy/UAi1jRwfPBDepgbViN5ueH3qUADAO9tTkJRZctPXVNU14IVvkyEEMDvCH3cMad8Nh20U+vG5OdjiTE4Z3tx6ttXjv0rIwK6z+bBVyPH+3PB2X/qskMvwj1nDYauQY/e5fGw8dvPlo7W/peFo+lU4qpRYcd+tNRdbAsMLERFZxcggN/zw5GgM8XFGYUUtHl2XgMra9vWg1Gt1eGf7OQDAwrF90d/Tsc2vfTAqEFNCfNCgE3j662Moq2n+6iSDFT+fQ3pRFXw1arzaGHzay1ujxnv3hwIA/ncoHdtO5jR7XGp+uXHp5sVJgzDE59aunh3g5YTFsfrlo7/+dBr5ZS0vH10sqMC7v6QAAF6ZMgT+rva3dE5LYHghIiKrcXOwxWcLRsJHo8algkq8/uPpdr1+fUIGLhVWwt3BFn+8vW+7XiuTyfDWzBD4u9ohs7gar2461eLSym+phcaN4d6ZFQpntU27znW98YM88cR4/XLMS9+dQHqRaf9LdZ0WT3153Lh0Y+iVuVV/ur0vQvw0KGtl+ahBq8ML3yajtkGHsQM8MHdkQIfO2dkYXoiIyKo8ndVYOScMchnwXeIV/NxCNeJGFbUNWLnrAgBgcewAON1CoNDY2eD9ueFQyGXYnJTd7NJKaXU9XvxO3zT78G29MWaAR7vPc6Pn7xyIyN6uKK9twKKvjhv7X3Q6gZc3nkBKXjk8HFX455ywVvt32kKpkOPd2cNho5Bh19l8bE7KbnLMh79exPGMEjiplHj7vuFddrnIgOGFiIisLrqvO54c3x8A8JefTqP8Jks4APDx3osoqqxDHw8HPBAVeMvnjujtimfv0C+tvLb5FC4VVBifE0Jg2cYTyCqpRm93+3Y3zbZE2dj/4mJvg5NZpZjz8SFsO5mDP32RiM1J2VDIZfhgbhh6OanMcr7B3s54pvHqo9d/PI3M4mv9RYnpxXg/Th8C/3bvMPi62JnlnJ2J4YWIiLqERb/rj97u9sgrq8V7O1vfzC2vrAaf7E8DALw0aRBsFB37OHtyQn9E93FDZZ0WD3+agPSiSuh0Aiu2n8O2k7lQymX4YG64WW9G6Otih08eiYSzWonkzBI8+eUx7DyTB1uFHP+cE4ZR/Tte4bne4+P7Ybi/BqXV9XhkTQLSCiuRmF6MP3x2FFqdwD3DfXBvF9vPpSW8PQAREXUZ+84X4JE1CZDLgB8XjUGwn6bZ45ZtPIGvEzIR0dsV3z0eY5ZljvzyGsz9zyFcKqyE2kYOL2c10huvgHprRggejL716k5rrlytwqo9F3HiSgmC3B3w5IR+GObb/Lw7Kq+sBjM/PIiskmqTx0P9Nfhq4W1WvVN0ez6/GV6IiKhLWfTVMWw5kYPQABdsfGJUk03gLuSVY+LKfdAJ4PsnYhDRu/UdZ9sjr6wGT399HAlpxQAAOxsF/jptGO7v4g2s7ZFVUo1lG09i3/kCyGXAvWF++Ov0YbfUM2RODC8ML0REkpVXVoM7/t9eVNQ24I17g/HQbb2NzwkhMH/tEew7X4CJw7zwn4cjzX5+IQSSMktQWFGHyN6u7dqITkquVtZBoZB16Mopc+KNGYmISLK8nNV44a6BAIB3tp8z2Zsk7mw+9p0vgK1CjmWTh3TK+WUyGcIDXXHnUK9uG1wAwNXBtssEl/ZieCEioi7n4Zgg494kT355DHUNOuSX12DZDycBAI+N6YMgDwcrj5KshctGRETUJV0qqMD0f/+G8toGDPFxRkVtPTKLqzHQyxGbnxoDO9v2bZdPXRuXjYiISPL69nLERw9FwEmtxNmcMmQWV8NXo8bHD0cyuPRw1rsmioiI6CbGDPDAjudux5bkHKhtFbg3zNfqV8WQ9TG8EBFRl+ajscPCdt63iLo3LhsRERGRpDC8EBERkaQwvBAREZGkMLwQERGRpDC8EBERkaQwvBAREZGkMLwQERGRpDC8EBERkaQwvBAREZGkMLwQERGRpDC8EBERkaQwvBAREZGkMLwQERGRpHS7u0oLIQAAZWVlVh4JERERtZXhc9vwOd6abhdeysvLAQABAQFWHgkRERG1V3l5OTQaTavHyERbIo6E6HQ6ZGdnw8nJCTKZzKzvXVZWhoCAAGRmZsLZ2dms790VcH7S193n2N3nB3T/OXJ+0tdZcxRCoLy8HL6+vpDLW+9q6XaVF7lcDn9//049h7Ozc7f9pQQ4v+6gu8+xu88P6P5z5PykrzPmeLOKiwEbdomIiEhSGF6IiIhIUhhe2kGlUuH111+HSqWy9lA6Becnfd19jt19fkD3nyPnJ31dYY7drmGXiIiIujdWXoiIiEhSGF6IiIhIUhheiIiISFIYXoiIiEhSGF7aaNWqVQgKCoJarUZ0dDQSEhKsPaQ2Wb58OUaOHAknJyd4enri3nvvRUpKiskxNTU1eOqpp+Du7g5HR0fcd999yMvLMzkmIyMDU6ZMgb29PTw9PbF06VI0NDRYciptsmLFCshkMjz77LPGx7rD/LKysvDQQw/B3d0ddnZ2CAkJwdGjR43PCyHw2muvwcfHB3Z2doiNjcWFCxdM3qO4uBjz5s2Ds7MzXFxc8Pvf/x4VFRWWnkoTWq0Wr776Kvr06QM7Ozv069cPf//7303ubyK1+e3btw9Tp06Fr68vZDIZNm3aZPK8ueZz4sQJjB07Fmq1GgEBAXjnnXc6e2oAWp9ffX09XnrpJYSEhMDBwQG+vr545JFHkJ2dbfIeUp3fjR5//HHIZDKsXLnS5PGuPD+gbXM8e/Yspk2bBo1GAwcHB4wcORIZGRnG5636t1XQTa1fv17Y2tqKNWvWiNOnT4uFCxcKFxcXkZeXZ+2h3dTEiRPF2rVrxalTp0RSUpK4++67RWBgoKioqDAe8/jjj4uAgAARFxcnjh49Km677TYxatQo4/MNDQ0iODhYxMbGiuPHj4tt27YJDw8PsWzZMmtMqUUJCQkiKChIDB8+XCxevNj4uNTnV1xcLHr37i0effRRcfjwYXHp0iXxyy+/iNTUVOMxK1asEBqNRmzatEkkJyeLadOmiT59+ojq6mrjMZMmTRKhoaHi0KFDYv/+/aJ///7igQcesMaUTLz55pvC3d1dbNmyRaSlpYlvv/1WODo6ivfff994jNTmt23bNvHKK6+IjRs3CgDihx9+MHneHPMpLS0VXl5eYt68eeLUqVPi66+/FnZ2duI///mPVedXUlIiYmNjxYYNG8S5c+dEfHy8iIqKEhERESbvIdX5XW/jxo0iNDRU+Pr6in/+858mz3Xl+Qlx8zmmpqYKNzc3sXTpUnHs2DGRmpoqNm/ebPK5Z82/rQwvbRAVFSWeeuop4/darVb4+vqK5cuXW3FUtyY/P18AEHv37hVC6P/Q2NjYiG+//dZ4zNmzZwUAER8fL4TQ/5LL5XKRm5trPOajjz4Szs7Oora21rITaEF5ebkYMGCA2Llzpxg3bpwxvHSH+b300ktizJgxLT6v0+mEt7e3ePfdd42PlZSUCJVKJb7++mshhBBnzpwRAMSRI0eMx/z8889CJpOJrKyszht8G0yZMkU89thjJo/NnDlTzJs3Twgh/fnd+MFgrvl8+OGHwtXV1eR39KWXXhKDBg3q5BmZau3D3SAhIUEAEOnp6UKI7jG/K1euCD8/P3Hq1CnRu3dvk/AipfkJ0fwc58yZIx566KEWX2Ptv61cNrqJuro6JCYmIjY21viYXC5HbGws4uPjrTiyW1NaWgoAcHNzAwAkJiaivr7eZH6DBw9GYGCgcX7x8fEICQmBl5eX8ZiJEyeirKwMp0+ftuDoW/bUU09hypQpJvMAusf8fvzxR0RGRmL27Nnw9PREeHg4PvnkE+PzaWlpyM3NNZmjRqNBdHS0yRxdXFwQGRlpPCY2NhZyuRyHDx+23GSaMWrUKMTFxeH8+fMAgOTkZBw4cACTJ08GIP353chc84mPj8ftt98OW1tb4zETJ05ESkoKrl69aqHZtE1paSlkMhlcXFwASH9+Op0ODz/8MJYuXYphw4Y1eb47zG/r1q0YOHAgJk6cCE9PT0RHR5ssLVn7byvDy00UFhZCq9Wa/PABwMvLC7m5uVYa1a3R6XR49tlnMXr0aAQHBwMAcnNzYWtra/yjYnD9/HJzc5udv+E5a1u/fj2OHTuG5cuXN3muO8zv0qVL+OijjzBgwAD88ssveOKJJ/DMM8/gs88+A3BtjK39jubm5sLT09PkeaVSCTc3N6vP8eWXX8bcuXMxePBg2NjYIDw8HM8++yzmzZsHQPrzu5G55tPVf28Nampq8NJLL+GBBx4w3sRP6vN7++23oVQq8cwzzzT7vNTnl5+fj4qKCqxYsQKTJk3Cjh07MGPGDMycORN79+41jtGaf1u73V2lqWVPPfUUTp06hQMHDlh7KGaTmZmJxYsXY+fOnVCr1dYeTqfQ6XSIjIzEW2+9BQAIDw/HqVOnsHr1asyfP9/Ko+u4b775Bl9++SW++uorDBs2DElJSXj22Wfh6+vbLebXk9XX1+P++++HEAIfffSRtYdjFomJiXj//fdx7NgxyGQyaw+nU+h0OgDA9OnT8dxzzwEAwsLCcPDgQaxevRrjxo2z5vAAsPJyUx4eHlAoFE06qPPy8uDt7W2lUbXfokWLsGXLFuzZswf+/v7Gx729vVFXV4eSkhKT46+fn7e3d7PzNzxnTYmJicjPz8eIESOgVCqhVCqxd+9efPDBB1AqlfDy8pL0/ADAx8cHQ4cONXlsyJAhxq5/wxhb+x319vZGfn6+yfMNDQ0oLi62+hyXLl1qrL6EhITg4YcfxnPPPWespEl9fjcy13y6+u+tIbikp6dj586dxqoLIO357d+/H/n5+QgMDDT+zUlPT8fzzz+PoKAg4/ikOj9A/7mnVCpv+nfHmn9bGV5uwtbWFhEREYiLizM+ptPpEBcXh5iYGCuOrG2EEFi0aBF++OEH7N69G3369DF5PiIiAjY2NibzS0lJQUZGhnF+MTExOHnypMn/GQ1/jG785ba0O+64AydPnkRSUpLxKzIyEvPmzTP+t5TnBwCjR49ucnn7+fPn0bt3bwBAnz594O3tbTLHsrIyHD582GSOJSUlSExMNB6ze/du6HQ6REdHW2AWLauqqoJcbvqnSKFQGP/1J/X53chc84mJicG+fftQX19vPGbnzp0YNGgQXF1dLTSb5hmCy4ULF7Br1y64u7ubPC/l+T388MM4ceKEyd8cX19fLF26FL/88gsAac8P0H/ujRw5stW/O1b/7OhQu28PsX79eqFSqcS6devEmTNnxB//+Efh4uJi0kHdVT3xxBNCo9GIX3/9VeTk5Bi/qqqqjMc8/vjjIjAwUOzevVscPXpUxMTEiJiYGOPzhsvd7rrrLpGUlCS2b98uevXq1WUuJb7R9VcbCSH9+SUkJAilUinefPNNceHCBfHll18Ke3t78cUXXxiPWbFihXBxcRGbN28WJ06cENOnT2/20tvw8HBx+PBhceDAATFgwIAucan0/PnzhZ+fn/FS6Y0bNwoPDw/x4osvGo+R2vzKy8vF8ePHxfHjxwUA8d5774njx48br7Yxx3xKSkqEl5eXePjhh8WpU6fE+vXrhb29vUUutW1tfnV1dWLatGnC399fJCUlmfzduf4KE6nOrzk3Xm0kRNeenxA3n+PGjRuFjY2N+Pjjj8WFCxfEv/71L6FQKMT+/fuN72HNv60ML230r3/9SwQGBgpbW1sRFRUlDh06ZO0htQmAZr/Wrl1rPKa6ulo8+eSTwtXVVdjb24sZM2aInJwck/e5fPmymDx5srCzsxMeHh7i+eefF/X19RaeTdvcGF66w/x++uknERwcLFQqlRg8eLD4+OOPTZ7X6XTi1VdfFV5eXkKlUok77rhDpKSkmBxTVFQkHnjgAeHo6CicnZ3FggULRHl5uSWn0ayysjKxePFiERgYKNRqtejbt6945ZVXTD7opDa/PXv2NPv/u/nz5wshzDef5ORkMWbMGKFSqYSfn59YsWKF1eeXlpbW4t+dPXv2SH5+zWkuvHTl+QnRtjl++umnon///kKtVovQ0FCxadMmk/ew5t9WmRDXbWNJRERE1MWx54WIiIgkheGFiIiIJIXhhYiIiCSF4YWIiIgkheGFiIiIJIXhhYiIiCSF4YWIiIgkheGFiIiIJIXhhYiIiCSF4YWIiIgkheGFiIiIJIXhhYiIiCTl/wM97zVMn7+xCQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(bs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "indoor-above",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "rising-adrian",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tles/22335.txt    tles/boom.txt         tles/gps-ops.txt  tles/spire.txt\n",
      "tles/active.txt   tles/cubesat.txt      tles/planet.txt   tles/stations.txt\n",
      "tles/analyst.txt  tles/engineering.txt  tles/ses.txt\n"
     ]
    }
   ],
   "source": [
    "ls tles/*.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "swedish-projection",
   "metadata": {},
   "outputs": [],
   "source": [
    "tle_filename = 'tles/TWOLINE.TXT'\n",
    "tle_filename = 'tles/22335.txt'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "studied-government",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 VANGUARD 1\n",
      "1 00005U 58002B   22334.74432056  .00000101  00000-0  15418-3 0  9999\n",
      "2 00005  34.2490 346.9265 1845722  66.0949 312.1720 10.85038156302601\n",
      "0 VANGUARD 1\n",
      "1 00005U 58002B   22335.11262076  .00000100  00000-0  15292-3 0  9993\n",
      "2 00005  34.2490 345.7915 1845718  67.7532 310.8033 10.85038245302646\n",
      "0 VANGUARD 2\n",
      "1    11U 59001A   22335.21011001  .00000606  00000-0  30430-3 0  9993\n",
      "2    11  32.8791 287.3922 1465359 275.6560  67.9281 11.86232147723203\n",
      "0 VANGUARD 2\n"
     ]
    }
   ],
   "source": [
    "!head $tle_filename"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "silver-allah",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(tle_filename, 'r') as f:\n",
    "    \n",
    "    tle_lines = f.readlines()\n",
    "    \n",
    "    tles = []\n",
    "    \n",
    "    for cnt, ix in enumerate(range(0, len(tle_lines), 3)):\n",
    "        #print(cnt, ix, tle_lines[ix])\n",
    "        tles.append((tle_lines[ix].strip(), tle_lines[ix + 1].strip(), tle_lines[ix + 2].strip()))\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "manufactured-sword",
   "metadata": {},
   "outputs": [],
   "source": [
    "if tle_filename == 'tles/TWOLINE.TXT':\n",
    "    with open(tle_filename, 'r') as f:\n",
    "\n",
    "        tle_lines = f.readlines()\n",
    "\n",
    "        tles = []\n",
    "        tle = []\n",
    "\n",
    "        for line in tle_lines:\n",
    "            if line.startswith('#'):\n",
    "                continue\n",
    "                \n",
    "            tle.append(line.strip())\n",
    "            \n",
    "            if line.startswith('2 '):\n",
    "                tles.append((tle[0][2:9], tle[0], tle[1]))\n",
    "                tle = []\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "willing-sight",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('0 VANGUARD 1',\n",
       "  '1 00005U 58002B   22334.74432056  .00000101  00000-0  15418-3 0  9999',\n",
       "  '2 00005  34.2490 346.9265 1845722  66.0949 312.1720 10.85038156302601'),\n",
       " ('0 VANGUARD 1',\n",
       "  '1 00005U 58002B   22335.11262076  .00000100  00000-0  15292-3 0  9993',\n",
       "  '2 00005  34.2490 345.7915 1845718  67.7532 310.8033 10.85038245302646'),\n",
       " ('0 VANGUARD 2',\n",
       "  '1    11U 59001A   22335.21011001  .00000606  00000-0  30430-3 0  9993',\n",
       "  '2    11  32.8791 287.3922 1465359 275.6560  67.9281 11.86232147723203')]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tles[:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "inside-timothy",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "fifty-freeware",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TLE                       Iter    Cov (m)   StdDev (m)   Res @ Epoch (m)   Res @ End (m)\n",
      "0 VANGUARD 1                 4 1172945.940     5.34e-08          1.29e-08        2.19e-08\n",
      "0 VANGUARD 1                 4 1160378.512     5.35e-08          6.32e-09        1.58e-08\n",
      "0 VANGUARD 2                 3 1459333.263     1.54e-05          2.32e-05        2.20e-05\n",
      "0 VANGUARD 2                 3 1473976.124     1.53e-05          2.21e-05        2.10e-05\n",
      "0 VANGUARD R/B               4 686523.319     7.40e-08          1.14e-08        3.52e-08\n",
      "0 VANGUARD R/B               4 674240.858     7.14e-08          2.44e-08        2.07e-08\n",
      "0 VANGUARD R/B               4 662323.305     7.40e-08          2.50e-08        3.45e-08\n",
      "0 VANGUARD R/B               4 1248879.038     4.87e-08          1.09e-08        2.95e-08\n",
      "0 VANGUARD R/B               4 1241693.806     4.51e-08          2.55e-08        2.67e-08\n",
      "0 VANGUARD 3                 3 1373927.036     1.80e-05          3.49e-05        3.65e-05\n",
      "0 VANGUARD 3                 3 1366386.578     1.78e-05          3.59e-05        3.74e-05\n",
      "0 VANGUARD 3                 3 1366386.578     1.78e-05          3.59e-05        3.74e-05\n",
      "0 EXPLORER 7                 4 893840.370     7.52e-07          1.72e-06        1.76e-06\n",
      "0 EXPLORER 7                 4 889234.928     7.67e-07          1.76e-06        1.79e-06\n",
      "0 EXPLORER 7                 4 883494.273     7.85e-07          1.80e-06        1.84e-06\n",
      "0 TIROS 1                    4 970207.386     3.39e-04          6.14e-04        5.86e-04\n",
      "0 TIROS 1                    4 965170.648     3.93e-04          6.96e-04        6.64e-04\n",
      "0 TIROS 1                    4 964822.023     3.97e-04          7.02e-04        6.70e-04\n",
      "0 TRANSIT 2A                 4 1221992.786     1.12e-07          2.70e-07        2.59e-07\n",
      "0 SOLRAD 1 (GREB)            3 682723.569     2.94e-04          3.83e-04        3.99e-04\n",
      "0 SOLRAD 1 (GREB)            3 681629.695     2.94e-04          3.83e-04        3.99e-04\n",
      "0 SOLRAD 1 (GREB)            3 681192.305     2.94e-04          3.82e-04        3.98e-04\n",
      "0 THOR ABLESTAR R/B          4 1260638.347     1.32e-07          2.79e-07        2.75e-07\n",
      "0 THOR ABLESTAR R/B          4 1260937.490     1.35e-07          2.79e-07        2.64e-07\n",
      "0 DELTA 1 R/B                4 1008159.356     5.70e-07          1.46e-06        1.49e-06\n",
      "0 DELTA 1 R/B                4 1002163.889     5.43e-07          1.40e-06        1.42e-06\n",
      "0 ECHO 1 DEB (METAL OBJ)     4 1198266.518     2.96e-07          8.54e-07        8.37e-07\n",
      "0 ECHO 1 DEB (METAL OBJ)     4 1197012.063     2.67e-07          7.66e-07        7.55e-07\n",
      "0 ECHO 1 DEB (METAL OBJ)     4 1000064.187     3.15e-07          4.28e-07        4.43e-07\n",
      "0 COURIER 1B                 4 1223532.800     1.78e-07          5.27e-07        5.19e-07\n",
      "0 COURIER 1B                 4 1231194.040     1.68e-07          5.05e-07        5.01e-07\n",
      "0 THOR ABLESTAR R/B          4 1076916.853     2.44e-07          6.82e-07        6.53e-07\n",
      "0 THOR ABLESTAR R/B          4 1080338.719     2.42e-07          6.77e-07        6.55e-07\n",
      "0 SCOUT X-1 R/B              4 842937.398     9.00e-08          3.89e-08        2.27e-08\n",
      "0 SCOUT X-1 R/B              4 811236.587     9.07e-08          2.55e-08        2.11e-08\n",
      "0 SCOUT X-1 DEB              3 1214131.655     6.36e-05          1.00e-04        9.88e-05\n",
      "0 EXPLORER 11                4 1160849.247     5.36e-08          1.04e-08        1.71e-08\n",
      "0 EXPLORER 11                4 1151600.656     7.07e-08          8.81e-09        8.45e-09\n",
      "0 EXPLORER 11                4 1139818.558     7.81e-08          2.19e-08        1.66e-08\n",
      "0 THOR ABLE DEB (YO)         4 984884.258     5.87e-05          1.50e-04        1.50e-04\n",
      "0 THOR ABLE DEB (YO)         4 984884.258     5.87e-05          1.50e-04        1.50e-04\n",
      "0 TRANSIT 4A                 4 991121.374     2.89e-06          7.06e-06        7.19e-06\n",
      "0 TRANSIT 4A                 4 991892.224     2.86e-06          6.98e-06        7.11e-06\n",
      "0 SOLRAD 3/INJUN 1           4 938506.257     3.43e-06          8.54e-06        8.66e-06\n",
      "0 THOR ABLESTAR R/B          4 972304.877     2.51e-06          6.12e-06        6.24e-06\n",
      "0 THOR ABLESTAR R/B          4 973634.924     2.47e-06          6.04e-06        6.15e-06\n",
      "0 THOR ABLESTAR DEB          4 895449.083     1.51e-06          3.73e-06        3.79e-06\n",
      "0 THOR ABLESTAR DEB          4 896912.488     1.51e-06          3.73e-06        3.78e-06\n",
      "0 THOR ABLESTAR DEB          3 912879.936     4.49e-04          7.26e-04        7.60e-04\n",
      "0 THOR ABLESTAR DEB          3 912879.936     4.49e-04          7.26e-04        7.60e-04\n",
      "0 THOR ABLESTAR DEB          4 871910.471     2.31e-06          4.93e-06        5.10e-06\n",
      "0 THOR ABLESTAR DEB          4 871440.714     2.32e-06          4.95e-06        5.11e-06\n",
      "0 THOR ABLESTAR DEB          4 845095.875     1.77e-06          3.65e-06        3.77e-06\n",
      "0 THOR ABLESTAR DEB          4 845399.966     1.66e-06          3.37e-06        3.50e-06\n",
      "0 THOR ABLESTAR DEB          4 845399.966     1.66e-06          3.37e-06        3.50e-06\n",
      "0 THOR ABLESTAR DEB          4 822040.866     7.13e-07          1.77e-06        1.79e-06\n",
      "0 THOR ABLESTAR DEB          4 1130139.536     2.18e-06          5.40e-06        5.47e-06\n",
      "0 THOR ABLESTAR DEB          4 1130051.515     2.21e-06          5.47e-06        5.55e-06\n",
      "0 THOR ABLESTAR DEB          4 1120296.582     3.09e-06          7.75e-06        7.87e-06\n",
      "0 THOR ABLESTAR DEB          4 1120153.946     3.10e-06          7.76e-06        7.85e-06\n",
      "0 THOR ABLESTAR DEB          4 1119603.654     3.10e-06          7.77e-06        7.88e-06\n",
      "0 THOR ABLESTAR DEB          4 903903.729     4.73e-06          1.18e-05        1.19e-05\n",
      "0 THOR ABLESTAR DEB          4 904587.703     4.79e-06          1.19e-05        1.21e-05\n",
      "0 THOR ABLESTAR DEB          4 1096442.730     3.70e-06          9.22e-06        9.33e-06\n",
      "0 THOR ABLESTAR DEB          4 1095948.604     3.68e-06          9.18e-06        9.31e-06\n",
      "0 THOR ABLESTAR DEB          4 1095268.506     3.66e-06          9.13e-06        9.24e-06\n",
      "0 THOR ABLESTAR DEB          4 902595.992     1.84e-06          4.60e-06        4.66e-06\n",
      "0 THOR ABLESTAR DEB          4 902841.083     1.84e-06          4.61e-06        4.67e-06\n",
      "0 THOR ABLESTAR DEB          4 903783.812     1.86e-06          4.62e-06        4.70e-06\n",
      "0 THOR ABLESTAR DEB          4 933750.581     2.65e-06          6.56e-06        6.68e-06\n",
      "0 THOR ABLESTAR DEB          4 934507.808     2.65e-06          6.57e-06        6.67e-06\n",
      "0 THOR ABLESTAR DEB          4 935497.066     2.66e-06          6.58e-06        6.66e-06\n",
      "0 THOR ABLESTAR DEB          4 1050806.382     4.41e-07          1.16e-06        1.14e-06\n",
      "0 THOR ABLESTAR DEB          4 1048487.947     3.99e-07          1.06e-06        1.04e-06\n",
      "0 THOR ABLESTAR DEB          4 1048487.947     3.99e-07          1.06e-06        1.04e-06\n",
      "0 THOR ABLESTAR DEB          4 1112602.683     1.26e-06          2.51e-06        2.61e-06\n",
      "0 THOR ABLESTAR DEB          4 1113025.881     1.22e-06          2.40e-06        2.51e-06\n",
      "0 THOR ABLESTAR DEB          4 1113025.881     1.22e-06          2.40e-06        2.51e-06\n",
      "0 THOR ABLESTAR DEB          4 1006321.478     4.63e-06          1.12e-05        1.14e-05\n",
      "0 THOR ABLESTAR DEB          4 1006802.484     4.53e-06          1.10e-05        1.12e-05\n",
      "0 THOR ABLESTAR DEB          4 1006802.484     4.53e-06          1.10e-05        1.12e-05\n",
      "0 THOR ABLESTAR DEB          4 1063653.139     2.90e-07          5.92e-07        6.20e-07\n",
      "0 THOR ABLESTAR DEB          4 1064018.802     3.07e-07          6.25e-07        6.64e-07\n",
      "0 THOR ABLESTAR DEB          4 920557.505     7.91e-06          2.00e-05        2.02e-05\n",
      "0 THOR ABLESTAR DEB          4 921530.589     8.06e-06          2.04e-05        2.06e-05\n",
      "0 THOR ABLESTAR DEB          4 921530.589     8.06e-06          2.04e-05        2.06e-05\n",
      "0 THOR ABLESTAR DEB          4 1028625.005     1.23e-06          2.93e-06        3.00e-06\n",
      "0 THOR ABLESTAR DEB          4 1028404.854     1.37e-06          3.27e-06        3.34e-06\n",
      "0 THOR ABLESTAR DEB          4 1028283.450     1.45e-06          3.45e-06        3.53e-06\n",
      "0 THOR ABLESTAR DEB          4 915814.062     9.47e-06          2.07e-05        2.13e-05\n",
      "0 THOR ABLESTAR DEB          4 915767.296     9.44e-06          2.07e-05        2.13e-05\n",
      "0 THOR ABLESTAR DEB          4 1010707.446     1.90e-05          4.95e-05        4.92e-05\n",
      "0 THOR ABLESTAR DEB          4 895425.279     5.74e-06          1.29e-05        1.32e-05\n",
      "0 THOR ABLESTAR DEB          4 895346.989     5.61e-06          1.26e-05        1.29e-05\n",
      "0 THOR ABLESTAR DEB          4 895346.989     5.61e-06          1.26e-05        1.29e-05\n",
      "0 THOR ABLESTAR DEB          4 1055305.981     3.77e-06          9.50e-06        9.60e-06\n",
      "0 THOR ABLESTAR DEB          4 1054449.375     3.65e-06          9.19e-06        9.29e-06\n",
      "0 THOR ABLESTAR DEB          4 1054104.820     3.60e-06          9.08e-06        9.17e-06\n",
      "0 THOR ABLESTAR DEB          4 1095531.182     1.89e-07          4.69e-07        4.57e-07\n",
      "0 THOR ABLESTAR DEB          4 1108716.680     4.95e-06          1.24e-05        1.26e-05\n",
      "0 THOR ABLESTAR DEB          4 1108568.571     4.95e-06          1.24e-05        1.26e-05\n"
     ]
    }
   ],
   "source": [
    "# LM Test\n",
    "\n",
    "limit = 100#9*1000000\n",
    "failed_tles = []\n",
    "\n",
    "print(f'{\"TLE\":24s} {\"Iter\":>5s} {\"Cov (m)\":>10s}   {\"StdDev (m)\":>10s}   {\"Res @ Epoch (m)\":>13s}   {\"Res @ End (m)\":>13s}')\n",
    "\n",
    "for idx, tle in enumerate(tles):\n",
    "\n",
    "    line1 = tle[1]\n",
    "    line2 = tle[2]\n",
    "    satellite = EarthSatellite(line1, line2, tle[0], ts)\n",
    "    \n",
    "    try:\n",
    "        iterations, sigma, sigmas, dxs, bs, lamdas, b_epoch, b_new_epoch, b, P, A = fit_method(satellite, fit_span=1, max_iter=25, lamda=1e-3, rms_epsilon=0.0001, debug=False)\n",
    "        \n",
    "        print(f'{tle[0]:24s} {iterations:5d} {np.sqrt(np.diag(P)[0]) * 1000:10.3f} {sigma * 1000:12.2e} {np.linalg.norm(b_new_epoch[0:3]) * 1000:17.2e} {np.linalg.norm(b[0:3]) * 1000:15.2e}')\n",
    "    except:\n",
    "        print(f'{tle[0]:24s} Failed')\n",
    "    \n",
    "    if idx == limit:\n",
    "        break\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e53cf3a3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.10.6 (main, Nov 14 2022, 16:10:14) [GCC 11.3.0]"
  },
  "vscode": {
   "interpreter": {
    "hash": "422d96f55d2b279cc3694c0d7101cd7135f83d32ba2b51a1c4e4fd4874965f43"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}