TLEFit - COE - FD.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": [],
   "source": [
    "from tlefit_coe_fd 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": 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): 771.301\n",
      "Condition number (ATWA_acc): 95944535.69933206\n",
      "Lambda:  0.001\n",
      "Residuals after/before 196 < 8.2e+06\n",
      "Covariance a: 0.002 m\n",
      "dx  [-6.32411703e+00  2.91139243e-04  8.86586783e-05  1.50200078e-06\n",
      " -1.39219881e-01  1.39216166e-01 -9.28583744e-04]\n",
      "COE elements = [ 7.04852742e+03  8.25456845e-03  1.70566489e+00  4.28239371e+00\n",
      "  5.48211003e+00  7.91329127e-01 -9.27583744e-04]\n",
      "\n",
      "Current Solution 7048.527416026275 0.0082546 97.7274 245.3631 314.1018 45.3398 -9.2758e-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): 838.979\n",
      "Condition number (ATWA_acc): 87559236.39803033\n",
      "Lambda:  0.001\n",
      "Residuals after/before 0.00342 < 196\n",
      "Covariance a: 0.002 m\n",
      "dx  [-4.05968457e-03  6.07709790e-05  1.04805281e-08 -1.50053418e-06\n",
      "  4.95035141e-03 -4.95186098e-03  9.40579095e-04]\n",
      "COE elements = [7.04852336e+03 8.31533942e-03 1.70566490e+00 4.28239221e+00\n",
      " 5.48706039e+00 7.86377266e-01 1.29953509e-05]\n",
      "\n",
      "Current Solution 7048.523356341709 0.0083153 97.7274 245.3630 314.3854 45.0561 +1.2995e-05\n",
      "        Original 7048.523352797561 0.0083155 97.7274 245.3630 314.3836 45.0579 +1.3123e-05\n",
      "Residual (b) = [ 1.80663489e-01  4.73367185e-01  1.66587408e-01 -5.54445451e-05\n",
      "  5.99629261e-05  5.49625695e-04]\n",
      "Residuals (b) r = 0.533, v = 0.000556\n",
      "Convergence test: 0.995112, sigma_new(0.404) < sigma_old(82.7)\n",
      "\n",
      "#################### ITERATION 3 ####################\n",
      "\n",
      "Condition number (A): 829.520\n",
      "Condition number (ATWA_acc): 87393335.96840487\n",
      "Lambda:  0.001\n",
      "Residuals after/before 3.37e-09 < 0.00342\n",
      "Covariance a: 0.002 m\n",
      "dx  [-3.54049394e-06  1.60411088e-07  7.79815677e-12 -1.47296889e-09\n",
      " -3.14203686e-05  3.14122822e-05  1.27462590e-07]\n",
      "COE elements = [7.04852335e+03 8.31549984e-03 1.70566490e+00 4.28239221e+00\n",
      " 5.48702897e+00 7.86408678e-01 1.31228135e-05]\n",
      "\n",
      "Current Solution 7048.523352801215 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.93963225e-04 -3.58733145e-04  4.13769088e-03  9.78375391e-07\n",
      "  2.05201780e-06 -2.55677631e-07]\n",
      "Residuals (b) r = 0.00423, v = 2.29e-06\n",
      "Convergence test: 0.995825, sigma_new(0.00169) < sigma_old(0.404)\n",
      "\n",
      "#################### ITERATION 4 ####################\n",
      "\n",
      "Condition number (A): 829.526\n",
      "Condition number (ATWA_acc): 87392633.04161324\n",
      "Lambda:  0.001\n",
      "Residuals after/before 3.42e-15 < 3.37e-09\n",
      "Covariance a: 0.002 m\n",
      "dx  [-3.65036923e-09  1.64412275e-10  7.74967164e-15 -1.46935244e-12\n",
      " -3.07335704e-08  3.07640952e-08  1.86529948e-10]\n",
      "COE elements = [7.04852335e+03 8.31550000e-03 1.70566490e+00 4.28239221e+00\n",
      " 5.48702893e+00 7.86408709e-01 1.31230000e-05]\n",
      "\n",
      "Current Solution 7048.523352797564 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.12386588e-07 -3.16924343e-07  4.38957989e-06  1.08588760e-09\n",
      "  2.30885586e-09 -1.86742177e-10]\n",
      "Residuals (b) r = 4.48e-06, v = 2.56e-09\n",
      "Convergence test: 0.999007, sigma_new(1.68e-06) < sigma_old(0.00169)\n",
      "\n",
      "Stopped due to dx convergence\n",
      "\n",
      "#################### SOLUTION IN 4 ITERATIONS ####################\n",
      "\n",
      "Solution 7048.523352797564 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.48e-06, v = 2.56e-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.436e-06 km 0.000e+00 km/s\n",
      "Residual at the end   4.475e-06 km 2.558e-09 km/s\n"
     ]
    }
   ],
   "source": [
    "fit_method = test_tle_fit_normalized\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": 7,
   "id": "separate-sleeping",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fdbba743b50>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/G0lEQVR4nO3dZ3hUZcLG8fvMpBFIQgkkhATpIC2hhgQIoAiIooCIoiJiw4BI8dXVbbivvrKurhQJoFgAKyiCBQQUhQAJIiX03lsSAqQCaXPeD+6yi4ImkORM+f+uaz54cmbmnuc6ztzXwzzPGKZpmgIAAHARNqsDAAAAlAblBQAAuBTKCwAAcCmUFwAA4FIoLwAAwKVQXgAAgEuhvAAAAJdCeQEAAC7Fy+oAZc3hcOjkyZMKCAiQYRhWxwEAACVgmqZycnIUFhYmm+2351bcrrycPHlSERERVscAAADX4NixYwoPD//Nc9yuvAQEBEj6+cUHBgZanAYAAJREdna2IiIiLn2O/xa3Ky///qeiwMBAygsAAC6mJF/54Au7AADApVBeAACAS6G8AAAAl0J5AQAALoXyAgAAXArlBQAAuBTKCwAAcCmUFwAA4FIoLwAAwKVQXgAAgEuhvAAAAJdCeQEAAC6F8lJC5/IKdN+sddpyLNPqKAAAeDTKSwm9tnyPkg6c0cAZSXp12W7lFxVbHQkAAI9EeSmhp3s1Vb/IMBU7TCX8cEB3vLFW209kWR0LAACPQ3kpoeqVffTGkDaacX9b1ajsoz1pObozYa1e/3avCoocVscDAMBjUF5K6dZWtbV8XJz6tgpVscPU1BX7dGfCWu08mW11NAAAPALl5RrUqOKr6fe307T72qiav7d2ncrWHdPWaOqKfSosZhYGAIDyRHm5Dre3DtPycd3Uu0WIihymXv92rwZMX6s9qTlWRwMAwG1RXq5TzQBfzXygnabcG6WgSt7afiJbt7+xWgk/7FcRszAAAJQ5yksZMAxDd0bV0bfj4tTzxloqLDb16rI9umtGkvalMQsDAEBZoryUoVqBfpr1YHv98+5IBfh5acvxLN02dY1mrjqgYodpdTwAANwC5aWMGYahu9qF69tx3dSjaU0VFDv09292a9DMJB04nWt1PAAAXB7lpZyEBvnp3Yc66B+DWivA10ubj2aq75TVenv1QWZhAAC4DpSXcmQYhga3j9CycXHq2jhY+UUOvbR4l+55M1mHMvKsjgcAgEuivFSAsKqVNPfhjpo4sJUq+9i14cg53TolUe+uOSQHszAAAJQK5aWCGIahIR3ratm4OHVuVEMXCx3636936t5Z63TkDLMwAACUFOWlgoVX89cHj0Trpf4t5e9j1/pDZ9Vn8mrNTT7MLAwAACVAebGAYRh6oNMNWjY2Tp0aVNeFwmL99Ysduv/tH3Xs7Hmr4wEA4NQoLxaKqO6vjx7tpL/d0UKVvO1KPnhGfSYn6sMfj8g0mYUBAOBKKC8Ws9kMDYutp2/GdFWHetWUV1CsPy3crqHvrNeJzAtWxwMAwOlQXpxEveDKmvd4jP5ye3P5etm0Zn+Gek9K1LyfjjILAwDAf6G8OBGbzdAjXerrmzFd1bZuVeXmF+kPC7bpofd+0qksZmEAAJAoL06pQc0q+vSJWP2xbzP5eNm0au9p9ZqUqE83HGMWBgDg8SgvTspuM/R4XEMteaqLIiOqKudikZ75bKsembNBadkXrY4HAIBlKC9OrlGtAC14IkZ/6NNMPnabvt+drlteX6WFm48zCwMA8EiUFxfgZbcpvntDff1UF7WqE6Tsi0UaN2+LHn9/o9JzmIUBAHgWyosLaRISoM9Hxup/ejWRt93QtzvT1GtSor5IOcEsDADAY1BeXIy33aYnb2qsL5/soua1A5V5vlBjPklR/AeblJGbb3U8AADKHeXFRd1YO1BfPNlZ43o2kZfN0NIdqeo1KVGLt56yOhoAAOWK8uLCvO02jenZWItGdVaz0ACdzSvQqI826cmPNulsXoHV8QAAKBeUFzfQsk6Qvnyyi0bf1Eh2m6Gvt55Sr0mrtHR7qtXRAAAoc5QXN+HjZdPTvZpq4chYNa5VRRm5BXrig40a88lmnWMWBgDgRigvbqZ1eFV9/VQXjezeUDZD+iLlpHpNTtS3O9OsjgYAQJmgvLghXy+7nu3TTJ+P7KyGNSvrdE6+Hpu7QePnpyjrfKHV8QAAuC6UFzcWFVFVi5/qqhFxDWQY0uebTqjX5FX6YXe61dEAALhmlBc35+dt1/N9b9RnT8SofnBlpWXna/jsn/TMp1uUfZFZGACA66G8eIh2N1TXkqe66pEu9WUY0qcbj6v3pESt2nva6mgAAJSKU5aXAQMGqFq1aho0aJDVUdxKJR+7/nJ7c80fEaMbavjrVNZFDXt3vZ7/fKtymIUBALgIpywvY8aM0dy5c62O4bY61Kuub8Z01UOx9SRJH68/pj6TV2vNvgxrgwEAUAJOWV66d++ugIAAq2O4NX8fL71wRwt9/FgnRVSvpBOZF/TAOz/qz4u2KS+/yOp4AABcVZmXl8TERPXr109hYWEyDEOLFi361TkJCQmqV6+e/Pz8FB0drfXr15d1DJRQTMMaWjomTkM73SBJ+mDdUfWenKjkA2csTgYAwJWVeXnJy8tTZGSkEhISrvj3efPmafz48ZowYYI2bdqkyMhI9e7dW+npLN+1SmVfL73Yv6U+fDRadapW0vFzFzRk1jq98OUOnS9gFgYA4FzKvLzceuuteumllzRgwIAr/v3111/XY489puHDh6t58+aaOXOm/P399e67717T8+Xn5ys7O/uyG65N50bBWjq2q4Z0rCtJmp10WLdOWa31h85anAwAgP+o0O+8FBQUaOPGjerZs+d/Aths6tmzp5KTk6/pMSdOnKigoKBLt4iIiLKK65EC/Lw1cWArzX24o2oH+enImfO6561k/e9XO3WhoNjqeAAAVGx5ycjIUHFxsUJCQi47HhISotTU//wCcs+ePXX33XdryZIlCg8P/81i8/zzzysrK+vS7dixY+WW35PENampZePiNLh9uExTenftIfWdulobjzALAwCwlpfVAa7ku+++K/G5vr6+8vX1Lcc0nivQz1v/GBSpW1vV1nMLtupQRp4GzUzWY10baPwtTeTnbbc6IgDAA1XozEtwcLDsdrvS0i7/heO0tDSFhoZWZBSUQo+mtbR8bDfd1fbnWZi3Eg/qtqmrtfnoOaujAQA8UIWWFx8fH7Vr104rVqy4dMzhcGjFihWKiYmpyCgopSB/b/1zcKTefrC9agb46sDpPN01I0mvLN2t/CK+CwMAqDhlXl5yc3OVkpKilJQUSdKhQ4eUkpKio0ePSpLGjx+vWbNmac6cOdq1a5fi4+OVl5en4cOHl3UUlIOezUP07bg49Y8Kk8OUZqw8oNunrtHW45lWRwMAeAjDNE2zLB9w5cqV6tGjx6+ODxs2TLNnz5YkTZs2Ta+++qpSU1MVFRWlqVOnKjo6ukyePzs7W0FBQcrKylJgYGCZPCaubNmOVP1p4TZl5BbIbjMU362hRt/cSL5efBcGAFA6pfn8LvPyYjXKS8U6m1egCV/u0FdbTkqSmoUG6LW7I9WyTpDFyQAArqQ0n99O+dtGcB3VK/vojSFtNP3+tqpe2Ue7U3PUP2GtJn27VwVFDqvjAQDcEOUFZaJvq9paPi5Ot7YMVZHD1JQV+9Q/Ya12nmTHYwBA2aK8oMwEV/HV9Pvb6o0hbVTV31s7T2XrzoQ1emPFPhUWMwsDACgblBeUKcMw1C8yTMvHxemW5iEqLDb1z2/3auD0JO1JzbE6HgDADVBeUC5qBfjpraHtNPmeKAVV8ta2E1nq98YaJfywX0XMwgAArgPlBeXGMAz1b1NHy8fF6eZmtVRQ7NCry/borhlJ2p/OLAwA4NpQXlDuQgL99Paw9vrn3ZEK8PPSluNZ6jt1jd5cdUDFDrdaqQ8AqACUF1QIwzB0V7twfTuum7o3ramCIocmfrNbd89M0oHTuVbHAwC4EMoLKlRokJ/ee6iD/nFXa1Xx9dKmo5nqO2W13l59kFkYAECJUF5Q4QzD0OAOEVo2Lk5dGwcrv8ihlxbv0r1vJetwRp7V8QAATo7yAsvUqVpJcx/uqJcHtFJlH7t+OnxOfaYkavbaQ3IwCwMAuArKCyxlGIbui66rpWPjFNuwhi4WOvTCVzs1ZNY6HT1z3up4AAAnRHmBU4io7q8PHonWi3e2UCVvu348dFZ9piTq/eTDzMIAAC5DeYHTsNkMDY2pp2Vj4xRdv7rOFxTrL1/s0APv/KhjZ5mFAQD8jPICp1O3hr8+fqyTXujXXH7eNiUdOKM+kxP10Y9HZZrMwgCAp6O8wCnZbIYe6lxfS8fEqf0N1ZRXUKw/LtymB99dr5OZF6yOBwCwEOUFTq1ecGXNGxGjP992o3y9bFq9L0O9JyVq/k/HmIUBAA9FeYHTs9sMPdq1gZaM6ao2dasqJ79Izy7YquGzf1Jq1kWr4wEAKhjlBS6jYc0q+uyJWP2xbzP5eNm0cs9p3TJplT7beJxZGADwIJQXuBS7zdDjcQ215KkuigwPUs7FIv3Pp1v06JwNSs9mFgYAPAHlBS6pUa0ALYiP1bN9msrHbtOK3em6ZVKiFm0+wSwMALg5ygtclpfdppHdG+mr0V3Uqk6Qsi4Uauy8FI14f6NO5+RbHQ8AUE4oL3B5TUMD9PnIWD19SxN52w0t35mmXpNW6astJ5mFAQA3RHmBW/C22zT65sb6YlQXNa8dqHPnCzX6480a9dEmncllFgYA3AnlBW6leVigFo3qrDE3N5aXzdCSbanqNSlR32w7ZXU0AEAZobzA7fh42TTuliZaNKqzmoUG6ExegeI/3KQnP9qks3kFVscDAFwnygvcVss6QfryyS4afVMj2W2Gvt56Sr0mrdKyHalWRwMAXAfKC9yaj5dNT/dqqoUjY9W4VhVl5BZoxPsbNfaTzco8zywMALgiygs8QuvwqvpqdBfFd28omyEtSjmpWyYlasWuNKujAQBKifICj+Hnbdcf+jTTgvhYNaxZWadz8vXInA16ev4WZV0otDoeAKCEKC/wOG3qVtPip7rq8bgGMgxpwabj6jVplX7Yk251NABACVBe4JH8vO36Y98b9dkTMaofXFlp2fka/t5P+sNnW5V9kVkYAHBmlBd4tHY3VNeSp7rq4c71ZRjSvA3H1GdSolbvO211NADAVVBe4PEq+dj1137NNe/xGN1Qw18nsy5q6Dvr9ceF25SbX2R1PADAL1BegH/pWL+6vhnTVQ/F1pMkffTjUfWelKi1+zOsDQYAuAzlBfgv/j5eeuGOFvr4sU4Kr1ZJJzIv6P63f9RfFm1XHrMwAOAUKC/AFcQ0rKFlY+P0QKe6kqT31x1RnymJWnfwjMXJAACUF+AqKvt66aX+rfTBI9GqU7WSjp29oHvfWqcXvtyh8wXMwgCAVSgvwO/o0jhYS8d21ZCOP8/CzE46rL5TVuunw2ctTgYAnonyApRAgJ+3Jg5spTkPd1TtID8dPnNeg99M1otf79TFwmKr4wGAR6G8AKXQrUlNLRsXp8Htw2Wa0jtrDqnvlNXaeOSc1dEAwGNQXoBSCvTz1j8GReq9hzooJNBXBzPydPfMJE1csotZGACoAJQX4Br1aFZLy8d208C2deQwpTcTD+r2N9Yo5Vim1dEAwK1RXoDrEOTvrdcHR2nWg+1VM8BX+9NzNXD6Wv1j6W7lFzELAwDlgfIClIFbmodo+dg43RkVJocpTV95QP3eWKNtx7OsjgYAbofyApSRapV9NOXeNpr5QDsFV/HR3rRc9Z++Vq8v36OCIofV8QDAbVBegDLWp2Wolo/rpttb11axw9TU7/frjmlrtOMkszAAUBYoL0A5qF7ZR9Pua6uE+9qqemUf7U7N0Z3T1mryd3tVWMwsDABcD8oLUI5ua11by8fFqU+LUBU5TE3+bp/6J6zVrlPZVkcDAJdFeQHKWXAVX814oK2mDmmjqv7e2nEyW3dMW6Np3+9TEbMwAFBqlBegAhiGoTsiw7R8XJxuaR6iwmJTry3fq4EzkrQ3LcfqeADgUigvQAWqFeCnt4a206R7IhXo56Wtx7N0+9Q1mrHygIodptXxAMAlUF6ACmYYhga0Cde347vppma1VFDs0CtLd+ueN5N19Mx5q+MBgNOjvAAWCQn00zvD2usfg1qriq+XNhw5pz5TEvXx+qMyTWZhAOBqKC+AhQzD0OD2EfpmTFd1rF9d5wuK9fzn2/TonA1Kz7lodTwAcEqUF8AJRFT31yePddKf+t4oH7tNK3anq/ekRC3dfsrqaADgdCgvgJOw2Qw9FtdAX47urBtrB+rc+UI98cEmjZ+fouyLhVbHAwCnQXkBnEyz0EAtGhWrkd0bymZIn286oVsnr1bygTNWRwMAp0B5AZyQr5ddz/ZppvkjYlS3ur9OZF7QkFnr9OLXO3WxsNjqeABgKcoL4MTa16uub8Z01ZCOEZKkd9YcUr831mj7CX7kEYDnorwATq6yr5cmDmytd4a1V3AVX+1Lz1X/hLX8vAAAj0V5AVzEzTeGaNnYrpd+5PG15Xs1+M1kHc7IszoaAFQoygvgQmr860ceXx8cqQBfL206mqlbp6zWB+uOsLEdAI9BeQFcjGEYGtg2XEvHxSmmQQ1dKCzWnxdt1/DZPyk9m43tALg/ygvgoupUraQPH43WX25vLh8vm1buOa1ekxO1eCsb2wFwb5QXwIXZbIYe6VJfi0d3Ucs6gco8X6hRH23S2E82K+sCG9sBcE+UF8ANNA4J0OfxnTX6pkayGdKilJPqMzlRa/dnWB0NAMoc5QVwEz5eNj3dq6k+i49VvRr+OpV1Ufe//aNe+HIHG9sBcCuUF8DNtK1bTUvGdNUDnepKkmYnHdZtU1dr6/FMa4MBQBmhvABuyN/HSy/1b6XZwzuoVoCvDpzO08DpSZryHRvbAXB9lBfAjXVvWkvLxsbptla1VeQwNem7vbprZrIOns61OhoAXDPKC+DmqlX20bT72mjKvVEK9PPSlmOZ6jt1teYmH2ZjOwAuifICeADDMHRnVB0tGxenLo2CdbHQob9+sUMPvrteqVlsbAfAtThleRkwYICqVaumQYMGWR0FcCu1gypp7sMd9UK/5vL1smn1vgz1npyoL7ectDoaAJSYU5aXMWPGaO7cuVbHANySzWbooc71tfiprmodHqSsC4V66uPNGv3xZmWeL7A6HgD8LqcsL927d1dAQIDVMQC31qhWFS2Ij9WYmxvLbjP01ZaT6j05UYl7T1sdDQB+U6nLS2Jiovr166ewsDAZhqFFixb96pyEhATVq1dPfn5+io6O1vr168siK4Ay5m23adwtTbQgPlYNgisrLTtfD767Xn/9YrsuFLCxHQDnVOrykpeXp8jISCUkJFzx7/PmzdP48eM1YcIEbdq0SZGRkerdu7fS09MvnRMVFaWWLVv+6nbyJP/uDlghKqKqFj/VVcNibpAkzU0+otumrlbKsUxrgwHAFRjmdayVNAxDCxcuVP/+/S8di46OVocOHTRt2jRJksPhUEREhEaPHq3nnnuuxI+9cuVKTZs2TZ999tlvnpefn6/8/PxL/52dna2IiAhlZWUpMDCwdC8IgBL3ntazn21VavZF2W2GRvVopNE3NZK33Sn/lRmAm8jOzlZQUFCJPr/L9N2ooKBAGzduVM+ePf/zBDabevbsqeTk5LJ8qksmTpyooKCgS7eIiIhyeR7AU8Q1qallY+N0R2SYih2mpq7Yp7tmJGl/OhvbAXAOZVpeMjIyVFxcrJCQkMuOh4SEKDU1tcSP07NnT919991asmSJwsPDf7P4PP/888rKyrp0O3bs2DXnB/CzIH9vTR3SRm8MaaOgSt7aejxLt01drffWHpLDwcZ2AKzlZXWAK/nuu+9KfK6vr698fX3LMQ3gufpFhqlDvep6dsFWJe49rb99tVPf7UrTq4MiFVa1ktXxAHioMp15CQ4Olt1uV1pa2mXH09LSFBoaWpZPBaCChAb5ac7wDnrxzhby87Zp7f4z6j05UYs2n+DnBQBYokzLi4+Pj9q1a6cVK1ZcOuZwOLRixQrFxMSU5VMBqECGYWhoTD0teaqrIiOqKudikcbOS9GTH23WuTw2tgNQsUpdXnJzc5WSkqKUlBRJ0qFDh5SSkqKjR49KksaPH69Zs2Zpzpw52rVrl+Lj45WXl6fhw4eXaXAAFa9BzSpa8ESMxt/SRF42Q4u3nVLvyYn6YU/6798ZAMpIqZdKr1y5Uj169PjV8WHDhmn27NmSpGnTpunVV19VamqqoqKiNHXqVEVHR5dJ4N9TmqVWAK7d1uOZGjcvRQdO50mS7o+uqz/ddqP8fZzyq3QAnFxpPr+va58XZ0R5ASrOxcJivbJ0t95be1iSVK+Gv/45OErtbqhmbTAALseyfV4AeBY/b7sm9GuhDx+NVu0gPx0+c153z0zSa8v2qKDIYXU8AG6K8gLgunVuFKylY+M0oE0dOUxp2g/7NWD6Wu1Ly7E6GgA3RHkBUCaCKnlr0j1Rmn5/W1X199aOk9m67Y01env1QTa2A1CmKC8AylTfVrW1fGycujetqYIih15avEv3v/2jTmResDoaADdBeQFQ5moF+um9hzro/wa0VCVvu5IPnlGfSYlasPE4G9sBuG6UFwDlwjAM3R99g74Z01Vt61ZVTn6Rnv50i+I/2KSzbGwH4DpQXgCUq3rBlTV/RIye6d1UXjZDS3ekqtekRK3Ylfb7dwaAK6C8ACh3XnabRvVopEWjOqtxrSrKyM3XI3M26LkFW5WbX2R1PAAuhvICoMK0rBOkr0Z30aNd6sswpE9+Oqa+U1Zrw+GzVkcD4EIoLwAqlJ+3XX++vbk+erST6lStpKNnz2vwm8l6ZeluNrYDUCKUFwCWiGlYQ9+M7aq72obLYUozVh7QnQlrtTs12+poAJwc5QWAZQL9vPXPwZGa+UA7Va/so12nsnXHG2v1VuIBFbOxHYCroLwAsFyflqFaNjZOPW+spYJih15esltDZq3TsbPnrY4GwAlRXgA4hZoBvpr1YHv9fWArVfaxa/2hs7p1ymrN33CMje0AXIbyAsBpGIahezvW1Tdj4tT+hmrKzS/Ss59t1ePvb1RGbr7V8QA4CcoLAKdTt4a/5o2I0R/6NJO33dC3O9PUe1Kilu9ItToaACdAeQHglOw2Q/HdG+qLUV3ULDRAZ/IK9Pj7G/XsZ1uUc7HQ6ngALER5AeDUmocF6osnO2tEXAMZhjR/w3HdOmW1fjx4xupoACxCeQHg9Hy97Hq+742a93iMwqtV0vFzF3TvrHV6ecku5RcVWx0PQAWjvABwGR3rV9fSsXG6p32ETFN6K/Gg7py2VjtPsrEd4EkoLwBcShVfL70yqLVmPdheNSr7aHdqju5MWKPpK/ezsR3gISgvAFzSLc1DtGxcnG5pHqLCYlP/WLpH97yZrKNn2NgOcHeUFwAuK7iKr94a2k6vDmqtKr5e2nDknPpMSdTH64+ysR3gxigvAFyaYRi6u32EvhnTVR3rV9f5gmI9//k2PTpng07nsLEd4I4oLwDcQkR1f33yWCf9qe+N8rHbtGJ3unpPTtTS7aesjgagjFFeALgNm83QY3EN9NXoLrqxdqDO5hXoiQ82afz8FGWzsR3gNigvANxO09AAfTGqs0Z2byibIX2+6YRunbxayQfY2A5wB5QXAG7Jx8umZ/s00/wRMapb3V8nMi9oyKx1evHrnbpYyMZ2gCujvABwa+3rVdc3Y7pqSMcISdI7aw6p3xtrtP1ElsXJAFwrygsAt1fZ10sTB7bWuw+1V3AVX+1Lz1X/hLWa9v0+FRU7rI4HoJQoLwA8xk3NQrR8XJxubRmqIoep15bv1eA3k3U4I8/qaABKgfICwKNUr+yj6fe31euDIxXg66VNRzN165TV+mDdETa2A1wE5QWAxzEMQwPbhmvpuDjFNKihC4XF+vOi7Ro++yelZ1+0Oh6A30F5AeCx6lStpA8fjdZfbm8uHy+bVu45rV6TE7V4KxvbAc6M8gLAo9lshh7pUl+LR3dRyzqByjxfqFEfbdLYTzYr6wIb2wHOiPICAJIahwTo8/jOGn1TI9kMaVHKSfWZnKi1+zOsjgbgFygvAPAvPl42Pd2rqT6Lj1W9Gv46lXVR97/9o174cgcb2wFOhPICAL/Qtm41LRnTVQ90qitJmp10WLdNXa2txzOtDQZAEuUFAK7I38dLL/VvpdnDO6hWgK8OnM7TwOlJmvIdG9sBVqO8AMBv6N60lpaNjdNtrWuryGFq0nd7ddfMZB08nWt1NMBjUV4A4HdUq+yjaUPaaMq9UQr089KWY5nqO3W15iYfZmM7wAKUFwAoAcMwdGdUHS0bF6cujYJ1sdChv36xQw++u16pWWxsB1QkygsAlELtoEqa+3BHvdCvuXy9bFq9L0O9Jyfqyy0nrY4GeAzKCwCUks1m6KHO9bX4qa5qHR6krAuFeurjzRr98WZlni+wOh7g9igvAHCNGtWqogXxsRpzc2PZbYa+2nJSvScnKnHvaaujAW6N8gIA18HbbtO4W5poQXysGgRXVlp2vh58d73++sV2XShgYzugPFBeAKAMREVU1eKnumpYzA2SpLnJR3Tb1NVKOZZpbTDADVFeAKCMVPKx6293ttT7j3RUaKCfDmbk6a4ZSXr9270qZGM7oMxQXgCgjHVtXFPLxsbpjsgwFTtMTV2xT3fNSNL+dDa2A8oC5QUAykGQv7emDmmjN4a0UVAlb209nqXbpq7We2sPyeFgYzvgelBeAKAc9YsM07KxcYprUlP5RQ797audGvrujzqZecHqaIDLorwAQDkLDfLTnOEd9OKdLeTnbdPa/WfUe3KiFm0+wc8LANeA8gIAFcAwDA2NqaclT3VVZERV5Vws0th5KXryo806l8fGdkBpUF4AoAI1qFlFC56I0fhbmsjLZmjxtlPqPTlRP+xJtzoa4DIoLwBQwbzsNj11c2MtHNlZDWtWVnpOvoa/95P+tHCbzhcUWR0PcHqUFwCwSKvwIC1+qquGd64nSfrwx6PqO2W1Nh09Z20wwMlRXgDAQn7edk3o10IfPhqt2kF+OnzmvAbNSNJry/aooIiN7YArobwAgBPo3ChYS8fGaUCbOnKY0rQf9mvA9LXal5ZjdTTA6VBeAMBJBFXy1qR7ojT9/raq6u+tHSezddsba/T26oNsbAf8F8oLADiZvq1qa/nYOPVoWlMFRQ69tHiX7n/7R51gYztAEuUFAJxSrUA/vftQB708oJX8fexKPnhGfSYlasHG42xsB49HeQEAJ2UYhu6LrqslT3VV27pVlZNfpKc/3aL4DzbpLBvbwYNRXgDAydULrqz5I2L0TO+m8rIZWrojVb0mJer73WlWRwMsQXkBABfgZbdpVI9GWjSqs5qEVFFGbr4enr1Bz3++VXn5bGwHz0J5AQAX0rJOkL58sose61pfhiF9vP6Ybp2yWhsOn7U6GlBhKC8A4GL8vO36023N9dGjnVSnaiUdPXteg99M1itLd7OxHTwC5QUAXFRMwxr6ZmxXDWoXLocpzVh5QHcmrNXu1GyrowHlivICAC4s0M9br90dqZkPtFP1yj7adSpbd7yxVm8lHlAxG9vBTVFeAMAN9GkZqmVj49TzxloqKHbo5SW7NWTWOp3Oybc6GlDmKC8A4CZqBvhq1oPt9feBrVTZx671h85qwPS12p+ea3U0oExRXgDAjRiGoXs71tVXo7uoXg1/HT93QXfNSNKPB89YHQ0oM5QXAHBDDWpW0YL4WLWpW1VZFwo19J31+iLlhNWxgDJBeQEAN1Wjiq8+fqyT+rQIVUGxQ2M+SdGMlQf4bSS4PKcrL5mZmWrfvr2ioqLUsmVLzZo1y+pIAOCy/LztSri/rR7pUl+S9MrS3frTou0qKmY/GLguw3SyCl5cXKz8/Hz5+/srLy9PLVu21IYNG1SjRo0S3T87O1tBQUHKyspSYGBgOacFANfx3tpD+t+vd8o0pR5Na2rafW1V2dfL6liApNJ9fjvdzIvdbpe/v78kKT8/X6ZpMsUJAGVgeOf6mvlAO/l52/TDntO6561kpWdftDoWUGqlLi+JiYnq16+fwsLCZBiGFi1a9KtzEhISVK9ePfn5+Sk6Olrr168v1XNkZmYqMjJS4eHheuaZZxQcHFzamACAK+jdIlQfP9ZJNSr7aPuJbA2YnqS9aTlWxwJKpdTlJS8vT5GRkUpISLji3+fNm6fx48drwoQJ2rRpkyIjI9W7d2+lp6dfOuff32f55e3kyZOSpKpVq2rLli06dOiQPvroI6Wl8bPvAFBW2tStpoUjO6tBcGWdyPx5KXXS/gyrYwEldl3feTEMQwsXLlT//v0vHYuOjlaHDh00bdo0SZLD4VBERIRGjx6t5557rtTPMXLkSN10000aNGjQFf+en5+v/Pz/7CCZnZ2tiIgIvvMCAL/jXF6BHn9/g346fE7edkOv3NVaA9uGWx0LHsqy77wUFBRo48aN6tmz53+ewGZTz549lZycXKLHSEtLU07Oz1OYWVlZSkxMVNOmTa96/sSJExUUFHTpFhERcX0vAgA8RLXKPnr/kWjd1rq2CotNjZ+/RVNX7ON7hnB6ZVpeMjIyVFxcrJCQkMuOh4SEKDU1tUSPceTIEXXt2lWRkZHq2rWrRo8erVatWl31/Oeff15ZWVmXbseOHbuu1wAAnsTP26437m2jEd0aSJJe/3av/rBgqwpZSg0n5nRr5Dp27KiUlJQSn+/r6ytfX9/yCwQAbs5mM/T8rTcqvJq/JnyxXfM3HNeprIuafn9bBfh5Wx0P+JUynXkJDg6W3W7/1Rds09LSFBoaWpZPBQAoY0M73aBZD7ZXJW+7Vu/L0N0zk3Uq64LVsYBfKdPy4uPjo3bt2mnFihWXjjkcDq1YsUIxMTFl+VQAgHJw840hmjeik4Kr+Gp3ao4GJCRp16lsq2MBlyl1ecnNzVVKSsqlf9o5dOiQUlJSdPToUUnS+PHjNWvWLM2ZM0e7du1SfHy88vLyNHz48DINDgAoH63Dq2rhyFg1qlVFqdkXdffMZCXuPW11LOCSUi+VXrlypXr06PGr48OGDdPs2bMlSdOmTdOrr76q1NRURUVFaerUqYqOji6TwL+HnwcAgLKRdb5QIz7YoHUHz8rLZujlga00uD0rOlE+SvP57XS/bXS9KC8AUHbyi4r17Gdb9UXKz5uIPnVzY43r2ViGYVicDO7GpX/bCADgPHy97Jp8T5Se7NFIkjR1xT49/ekWFRSxlBrWobwAAH6TYRj6n95N9feBrWS3Gfp80wk99N56ZV0otDoaPBTlBQBQIvd2rKt3hrVXZR+7kg6c0d0zk3Qik6XUqHiUFwBAiXVvWkvzn4hRSKCv9qblakDCWm0/kWV1LHgYygsAoFRahAVp4cjOahoSoPScfA1+M1k/7Em3OhY8COUFAFBqYVUr6dP4GHVuVEPnC4r16JwN+ujHo1bHgoegvAAArkmgn7fee6ij7mobrmKHqT8u3KZ/LN0th8OtduCAE6K8AACumY+XTa/d3VpjezaWJE1feUBj56Uov6jY4mRwZ5QXAMB1MQxDY3s20auDWsvLZujLLSc19J31yjxfYHU0uCnKCwCgTNzdPkKzh3dUgK+X1h86q7tmJOnY2fNWx4IborwAAMpMl8bB+jQ+RrWD/HTgdJ4GTF+rrcczrY4FN0N5AQCUqWahgVo4srNurB2ojNwC3fPmOn23M83qWHAjlBcAQJkLDfLTp0/EKK5JTV0oLNbj72/Q3OTDVseCm6C8AADKRRVfL70zrL3u7RAhhyn99YsdennJLpZS47pRXgAA5cbbbtPEga30P72aSJLeSjyo0R9v1sVCllLj2lFeAADlyjAMPXlTY02+J0redkOLt53SA2//qHN5LKXGtaG8AAAqRP82dTTn4Y4K8PPShiPnNHBGko6cybM6FlwQ5QUAUGFiGwbr8/hY1alaSYcy8jRwepI2HT1ndSy4GMoLAKBCNQ4J0MJRsWpZJ1Bn8go05K11Wro91epYcCGUFwBAhasV4Kd5j8fopma1lF/kUPyHG/XumkNWx4KLoLwAACxR2ddLbw1tp/uj68o0pf/9eqf+9tUOFbOUGr+D8gIAsIyX3aaX+rfUc7c2kyS9t/awRn64URcKWEqNq6O8AAAsZRiGnujWUG8MaSMfu03LdqRpyKx1ysjNtzoanBTlBQDgFPpFhumDR6MVVMlbKccyNXB6kg6ezrU6FpwQ5QUA4DQ61q+uz0fGKqJ6JR09e14DZyRpw+GzVseCk6G8AACcSsOaVfR5fGdFRlRV5vlC3ff2j1q89ZTVseBEKC8AAKdTM8BXnzzWSbc0D1FBkUOjPtqktxIPyDRZiQTKCwDASVXysWvmA+30UGw9SdLLS3ZrwpcspQblBQDgxOw2QxP6Ndefb7tRhiHNTT6iEe9v0PmCIqujwUKUFwCAUzMMQ492baDp97WVr5dN3+1K171vrVN6zkWro8EilBcAgEu4tVVtffRYJ1Xz99bW41kaOD1J+9NzrI4FC1BeAAAuo90N1fT5yM6qV8Nfx89d0MDpSVp38IzVsVDBKC8AAJdSP7iyPh/ZWW3rVlX2xSI9+M56fZFywupYqECUFwCAy6le2UcfPdZJt7YMVUGxQ2M+SdH0lftZSu0hKC8AAJfk521Xwn1t9WiX+pKkfyzdoz8u3K6iYofFyVDeKC8AAJdlsxn68+3N9UK/5jIM6eP1R/Xo3A3KzWcptTujvAAAXN5DnevrzQfayc/bppV7TuueN5OVls1SandFeQEAuIVeLUL1yeMxCq7iox0nszUgYa32pLKU2h1RXgAAbiMqoqo+j++sBjUr62TWRQ2akaSk/RlWx0IZo7wAANxK3Rr++jw+Vh3rVVdOfpGGvbdeCzYetzoWyhDlBQDgdqr6+2juIx3VLzJMhcWmnv50i6Z8t4+l1G6C8gIAcEt+3nZNuSdKT3RrKEma9N1ePfvZVhWylNrlUV4AAG7LZjP03K3N9FL/lrIZ0qcbj+vh2T8p52Kh1dFwHSgvAAC390CnG/T2sPby97Fr9b4M3T0zWaeyLlgdC9eI8gIA8Ag3NQvRvMdjVDPAV7tTczQgIUk7T2ZbHQvXgPICAPAYrcKDtHBkrBrXqqLU7Isa/GayEveetjoWSonyAgDwKOHV/PVZfKw6Naiu3PwiDZ/9k+b/dMzqWCgFygsAwOMEVfLWnIc7qn9UmIodpp5dsFWvL9/DUmoXQXkBAHgkXy+7Jt0TpdE3NZIkTf1+v56ev0UFRSyldnaUFwCAxzIMQ0/3aqpX7molu83Q55tPaNi765V1gaXUzozyAgDwePd0qKt3H+qgyj52JR88o7tnJulEJkupnRXlBQAASd2a1NT8J2IUEuirvWm56p+wVttPZFkdC1dAeQEA4F9ahAVp4cjOahYaoNM5+Rr8ZrJ+2J1udSz8AuUFAID/Ela1kuY/EaMujYJ1vqBYj87doI9+PGp1LPwXygsAAL8Q6Oet94Z30KB24Sp2mPrjwm16ZeluORwspXYGlBcAAK7A227Tq4Naa1zPJpKkGSsPaMy8FOUXFVucDJQXAACuwjAMjenZWK/dHSkvm6GvtpzU0HfWK/N8gdXRPBrlBQCA3zGoXbjmPNxRAb5eWn/orAbOSNKxs+etjuWxKC8AAJRA50bB+jQ+RrWD/HTwdJ4GTF+rLccyrY7lkSgvAACUULPQQC0a1VnNawcqI7dA9761Tt/uTLM6lsehvAAAUAohgX6a/0SMujWpqQuFxRrx/gbNTT5sdSyPQnkBAKCUqvh66e1h7XVvhwg5TOmvX+zQ/y3eyVLqCkJ5AQDgGnjbbZo4sJWe6d1UkjRr9SE9+fEmXSxkKXV5o7wAAHCNDMPQqB6NNOXeKHnbDS3Zlqr73/5RZ/NYSl2eKC8AAFynO6PqaO7D0Qr089LGI+d014wkHTmTZ3Ust0V5AQCgDMQ0rKHPR8aqTtVKOpSRpwHTk7Tp6DmrY7klygsAAGWkUa0ALRwVq1Z1gnQ2r0BD3lqnpdtPWR3L7VBeAAAoQ7UC/PTJ4510c7Nayi9yKP7DTXpnzSGrY7kVygsAAGWssq+X3hzaTg90qivTlF78eqde+HKHillKXSYoLwAAlAMvu00v3tlSz9/aTJI0O+mw4j/YqAsFLKW+XpQXAADKiWEYGtGtoabd10Y+XjYt35mmIbPWKSM33+poLo3yAgBAObu9dZg+fDRaVf29lXIsUwOnJ+ng6VyrY7ksygsAABWgQ73qWhAfq4jqlXT07HkNnJGknw6ftTqWS6K8AABQQRrWrKKFIzsrMqKqMs8X6v63f9TXW09aHcvlOGV5qVevnlq3bq2oqCj16NHD6jgAAJSZ4Cq++uSxTurVPEQFRQ49+dFmvbnqgEyTlUgl5WV1gKtJSkpSlSpVrI4BAECZq+Rj14wH2unFr3dqdtJhTfxmt46fu6AJ/ZrLy+6U8wpOhRECAMACdpuhF+5oob/c3lyGIb2/7ohGvL9R5wuKrI7m9EpdXhITE9WvXz+FhYXJMAwtWrToV+ckJCSoXr168vPzU3R0tNavX1+q5zAMQ926dVOHDh304YcfljYiAAAu45Eu9TXj/rby9bJpxe503fPmOqXnXLQ6llMrdXnJy8tTZGSkEhISrvj3efPmafz48ZowYYI2bdqkyMhI9e7dW+np6ZfOiYqKUsuWLX91O3ny5y8trVmzRhs3btSXX36pl19+WVu3br3GlwcAgPPr07K2Pnqsk6pX9tG2E1kakJCk/ek5VsdyWoZ5Hd8QMgxDCxcuVP/+/S8di46OVocOHTRt2jRJksPhUEREhEaPHq3nnnuu1M/xzDPPqEWLFnrooYeu+Pf8/Hzl5/9ns5/s7GxFREQoKytLgYGBpX4+AACscjgjT8Nn/6RDGXkK9PPSWw+2V6cGNayOVSGys7MVFBRUos/vMv3OS0FBgTZu3KiePXv+5wlsNvXs2VPJyckleoy8vDzl5PzcNnNzc/X999+rRYsWVz1/4sSJCgoKunSLiIi4vhcBAIBF6gVX1oL4WLW7oZqyLxZp6Ds/6ouUE1bHcjplWl4yMjJUXFyskJCQy46HhIQoNTW1RI+RlpamLl26KDIyUp06ddKDDz6oDh06XPX8559/XllZWZdux44du67XAACAlapX9tGHj0arb6tQFRabGvNJihJ+2M9S6v/idEulGzRooC1btpT4fF9fX/n6+pZjIgAAKpaft13ThrTV36vt1luJB/Xqsj06fu68XryzJUupVcYzL8HBwbLb7UpLS7vseFpamkJDQ8vyqQAAcGs2m6E/9r1Rf7ujhWyG9PH6Y3pkzgbl5rOUukzLi4+Pj9q1a6cVK1ZcOuZwOLRixQrFxMSU5VMBAOARhsXW05tD28vP26ZVe09r8MxkpWV79lLqUpeX3NxcpaSkKCUlRZJ06NAhpaSk6OjRo5Kk8ePHa9asWZozZ4527dql+Ph45eXlafjw4WUaHAAAT3FL8xDNezxGwVV8tPNUtgYkrNWeVM9dSl3qpdIrV6684u8NDRs2TLNnz5YkTZs2Ta+++qpSU1MVFRWlqVOnKjo6ukwC/57SLLUCAMCVHDt7XsPeW6+Dp/MU4OulmUPbqXOjYKtjlYnSfH5f1z4vzojyAgBwZ5nnC/T4+xu1/tBZedkM/f2u1hrULtzqWNfNsn1eAABA+arq76P3H+mofpFhKnKY+p9Pt2jKd/s8aik15QUAABfj62XXlHuiFN+9oSRp0nd79cxnW1VY7LA4WcWgvAAA4IJsNkN/6NNM/zegpWyG9NnG4xr+3k/KvlhodbRyR3kBAMCF3R99g94Z1kH+Pnat2Z+hwTOTdSrrgtWxyhXlBQAAF9ejWS3NHxGjmgG+2p2ao/4Ja7XzZLbVscoN5QUAADfQsk6QFo6MVeNaVZSWna/BbyZr1d7TVscqF5QXAADcRHg1f30WH6uYBjWUm1+kh2f/pHk/HbU6VpmjvAAA4EaCKnlrzsMdNbBNHRU7TP1hwTb9c/ket1pKTXkBAMDN+HjZ9M/BkXrqpkaSpDe+36/x87eooMg9llJTXgAAcEOGYWh8r6b6x12tZbcZWrj5hIa9u15ZF1x/KTXlBQAANza4Q4Tee6iDqvh6KfngGQ2akaTj585bHeu6UF4AAHBzcU1qav6IGIUG+mlfeq4GTE/S9hNZVse6ZpQXAAA8QPOwQC0cFatmoQE6nfPzUuofdqdbHeuaUF4AAPAQtYMqaf4TMeraOFjnC4r1yJyf9OGPR6yOVWqUFwAAPEign7fefaiD7m4XLocp/Wnhdv39m91yOFxnKTXlBQAAD+Ntt+kfg1pr/C1NJEkzVx3QmHkpyi8qtjhZyVBeAADwQIZh6KmbG+ufd0fKy2boqy0nNfTt9co8X2B1tN9FeQEAwIPd1S5ccx7uqABfL60/fFYDZyTp6BnnXkpNeQEAwMN1bhSsz+JjFRbkp4On8zRwxlptOZZpdayrorwAAAA1DQ3QwlGd1SIsUBm5BbrnrWQt35FqdawrorwAAABJUkign+aNiFG3JjV1sdChER9s1Jykw1bH+hXKCwAAuKSKr5feGdZeQzrWlWlKE77cof9bvNOpllJTXgAAwGW87Da9PKClnu3TVJI0a/UhPfnxJl0sdI6l1JQXAADwK4ZhaGT3Rppyb5R87DYt2Zaq+2at09k865dSU14AAMBV3RlVR3Mf6ahAPy9tOpqpgdPX6nBGnqWZKC8AAOA3dWpQQ5+PjFV4tUo6fOa8/rBgq6V5KC8AAOB3NaoVoM9HxurmZrX0z8GRlmbxsvTZAQCAy6gV4Kd3HupgdQxmXgAAgGuhvAAAAJdCeQEAAC6F8gIAAFwK5QUAALgUygsAAHAplBcAAOBSKC8AAMClUF4AAIBLobwAAACXQnkBAAAuhfICAABcCuUFAAC4FLf7VWnTNCVJ2dnZFicBAAAl9e/P7X9/jv8WtysvOTk5kqSIiAiLkwAAgNLKyclRUFDQb55jmCWpOC7E4XDo5MmTCggIkGEYZfrY2dnZioiI0LFjxxQYGFimj+1uGKuSY6xKjrEqOcaq5Bir0imv8TJNUzk5OQoLC5PN9tvfanG7mRebzabw8PByfY7AwEAu8BJirEqOsSo5xqrkGKuSY6xKpzzG6/dmXP6NL+wCAACXQnkBAAAuhfJSCr6+vpowYYJ8fX2tjuL0GKuSY6xKjrEqOcaq5Bir0nGG8XK7L+wCAAD3xswLAABwKZQXAADgUigvAADApVBeAACAS6G8/EJCQoLq1asnPz8/RUdHa/369b95/qeffqpmzZrJz89PrVq10pIlSyooqfVKM1azZ8+WYRiX3fz8/CowrTUSExPVr18/hYWFyTAMLVq06Hfvs3LlSrVt21a+vr5q1KiRZs+eXe45nUVpx2vlypW/uq4Mw1BqamrFBLbIxIkT1aFDBwUEBKhWrVrq37+/9uzZ87v388T3q2sZK099v5KkGTNmqHXr1pc2oIuJidE333zzm/ex4rqivPyXefPmafz48ZowYYI2bdqkyMhI9e7dW+np6Vc8PykpSUOGDNEjjzyizZs3q3///urfv7+2b99ewckrXmnHSvp5N8ZTp05duh05cqQCE1sjLy9PkZGRSkhIKNH5hw4d0m233aYePXooJSVFY8eO1aOPPqply5aVc1LnUNrx+rc9e/Zcdm3VqlWrnBI6h1WrVmnUqFFat26dvv32WxUWFqpXr17Ky8u76n089f3qWsZK8sz3K0kKDw/X3//+d23cuFEbNmzQTTfdpDvvvFM7duy44vmWXVcmLunYsaM5atSoS/9dXFxshoWFmRMnTrzi+YMHDzZvu+22y45FR0ebI0aMKNeczqC0Y/Xee++ZQUFBFZTOOUkyFy5c+JvnPPvss2aLFi0uO3bPPfeYvXv3Lsdkzqkk4/XDDz+Yksxz585VSCZnlZ6ebkoyV61addVzPPn96r+VZKx4v7pctWrVzLfffvuKf7PqumLm5V8KCgq0ceNG9ezZ89Ixm82mnj17Kjk5+Yr3SU5Ovux8Serdu/dVz3cX1zJWkpSbm6sbbrhBERERv9nkPZmnXlPXKyoqSrVr19Ytt9yitWvXWh2nwmVlZUmSqlevftVzuLZ+VpKxkni/kqTi4mJ98sknysvLU0xMzBXPseq6orz8S0ZGhoqLixUSEnLZ8ZCQkKv++3lqamqpzncX1zJWTZs21bvvvqsvvvhCH3zwgRwOh2JjY3X8+PGKiOwyrnZNZWdn68KFCxalcl61a9fWzJkztWDBAi1YsEARERHq3r27Nm3aZHW0CuNwODR27Fh17txZLVu2vOp5nvp+9d9KOlae/n61bds2ValSRb6+vnriiSe0cOFCNW/e/IrnWnVdud2vSsM5xcTEXNbcY2NjdeONN+rNN9/Uiy++aGEyuLKmTZuqadOml/47NjZWBw4c0KRJk/T+++9bmKzijBo1Stu3b9eaNWusjuL0SjpWnv5+1bRpU6WkpCgrK0ufffaZhg0bplWrVl21wFiBmZd/CQ4Olt1uV1pa2mXH09LSFBoaesX7hIaGlup8d3EtY/VL3t7eatOmjfbv318eEV3W1a6pwMBAVapUyaJUrqVjx44ec109+eST+vrrr/XDDz8oPDz8N8/11PerfyvNWP2Sp71f+fj4qFGjRmrXrp0mTpyoyMhITZky5YrnWnVdUV7+xcfHR+3atdOKFSsuHXM4HFqxYsVV/60vJibmsvMl6dtvv73q+e7iWsbql4qLi7Vt2zbVrl27vGK6JE+9pspSSkqK219XpmnqySef1MKFC/X999+rfv36v3sfT722rmWsfsnT368cDofy8/Ov+DfLrqty/Tqwi/nkk09MX19fc/bs2ebOnTvNxx9/3KxataqZmppqmqZpDh061Hzuuecunb927VrTy8vLfO2118xdu3aZEyZMML29vc1t27ZZ9RIqTGnH6m9/+5u5bNky88CBA+bGjRvNe++91/Tz8zN37Nhh1UuoEDk5OebmzZvNzZs3m5LM119/3dy8ebN55MgR0zRN87nnnjOHDh166fyDBw+a/v7+5jPPPGPu2rXLTEhIMO12u7l06VKrXkKFKu14TZo0yVy0aJG5b98+c9u2beaYMWNMm81mfvfdd1a9hAoRHx9vBgUFmStXrjRPnTp16Xb+/PlL5/B+9bNrGStPfb8yzZ//H1u1apV56NAhc+vWreZzzz1nGoZhLl++3DRN57muKC+/8MYbb5h169Y1fXx8zI4dO5rr1q279Ldu3bqZw4YNu+z8+fPnm02aNDF9fHzMFi1amIsXL67gxNYpzViNHTv20rkhISFm3759zU2bNlmQumL9eynvL2//Hpthw4aZ3bp1+9V9oqKiTB8fH7NBgwbme++9V+G5rVLa8XrllVfMhg0bmn5+fmb16tXN7t27m99//7014SvQlcZI0mXXCu9XP7uWsfLU9yvTNM2HH37YvOGGG0wfHx+zZs2a5s0333ypuJim81xXhmmaZvnO7QAAAJQdvvMCAABcCuUFAAC4FMoLAABwKZQXAADgUigvAADApVBeAACAS6G8AAAAl0J5AQAALoXyAgAAXArlBQAAuBTKCwAAcCmUFwAA4FL+Hx+vlDayOU21AAAAAElFTkSuQmCC",
      "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": 8,
   "id": "soviet-defensive",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fdbb84d6770>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOKElEQVR4nO3deXhU1cE/8O+dmcxk38kyWdj3QBIICSgI1CggsopSpRax1VqlLvhata3a5VVs+/58sS2K+lZxX7CACwhKQAIaQAgJOyQkQPaVZLJnMnN+f0zmhiGTzYSZ3Mz38zx5NPfeM+eciLlfznKvJIQQICIiIlIIlbMbQERERNQTDC9ERESkKAwvREREpCgML0RERKQoDC9ERESkKAwvREREpCgML0RERKQoDC9ERESkKBpnN6Cvmc1mFBYWwsfHB5IkObs5RERE1A1CCNTU1ECv10Ol6nxsZcCFl8LCQkRFRTm7GURERPQj5OXlITIystNrBlx48fHxAWDpvK+vr5NbQ0RERN1hMBgQFRUl38c7M+DCi3WqyNfXl+GFiIhIYbqz5IMLdomIiEhRGF6IiIhIURheiIiISFEYXoiIiEhRGF6IiIhIURheiIiISFEYXoiIiEhRGF6IiIhIURheiIiISFEYXoiIiEhRGF6IiIhIURheiIiISFEG3IsZr5WK2iaUGJrgppbgplahscWE6nojqhuMMDS2wNBgRE1jC2oajdCoVfBwU2OQjw7DBnlhnN4Xvu5uzu4CERHRgMDw0k3bjxfhmc9O/qiyKgmIjfLHzePCcNukCIT4uvdx64iIiFwHw0s36TRqhPjoYDSZYTQJaDUq+Hu4wdf65a6Br4cbfHQatJgF6ptNKKpuQFZJLQqqGnD0UhWOXqrC/3x9FjeOCcGvZw1HfHSAs7tFRESkOJIQQji7EX3JYDDAz88P1dXV8PX1dXZzAAAFVQ349mwpNqcX4MjFy/Lx6SOC8dDsEZg6LBCSJDmxhURERM7Vk/s3w4uDZZXU4LXUHGw5WgCT2fKjnxTtj4dmj8BPxoQwxBARkUtieOnH4cUqr7IeG/aex6bD+Wg2mQEAY8J8sCIpGrdMCEeQt65bn9PUYsLlOiNUKiDYSweViuGHiIiUh+FFAeHFqtTQiH/vz8V7By6irtkEAFCrJIwJ88G4cF8E++igUUloaDahusGIyrpmVNQ1o7L1q7apRf4sN7WEcXo/TB0WiEWxERin7//9JyIiAhheFBVerKrrjdh0JA+fZxbiWH51j8pqVBJMQuDq/5JxUf5Yc9MozBgZzOkoIiLq1xheFBherlRU3YDMvCpkldSioq4ZQgi4u6nh6+GGQC8tAr20CJL/qYOvh2WHU1FVI45cqsQ3p0rwzakSGE2W/7Q3jBqE5xfHICrQ08k9IyIiso/hReHhpS+U1TRhw97zePfARTS3mOGpVePJuWNw99TBXBdDRET9DsMLw4ssp6wWT20+jkO5lQAs27P/tmwi9P4eTm4ZERFRm57cv/luowFu2CBvfHTfVPxxwTi4u6mwP7scc9alYnN6PgZYbiUiIhfB8OICVCoJ91w/FNsfnoG4KH/UNLZgzSeZ+PV76aiobXJ284iIiHqE00YupsVkxoa957FuVxZazALB3lo8v2QCbh4X2uGOJCEESgxNyCqtQVFVI5paTPDSaTA4yAsTI/3gpmYGJiKi3uGaF4aXLp0oqMaaTzJwrqQWADAq1Bs3jwvDmHAfeGrVqKwzIqesFscLqnG8oBpV9Ua7n+Oj02BBnB6/umEYBgd5ObILREQ0gDC8MLx0S6PRhJdTsrDxuwtoMJo6vVatkjA4yBPRgZ5w16hhaDTiTHENKuuaAVgekPfwT0bigVnDORJDREQ9xvDC8NIjhkYjdhwvxoHcCuRXNqDe2IIATy0iAzwQE+GHiRH+GBXmDZ1GbVPOZBY4mFuBV789j31Z5QCA60cE4ZUVk+Hn4eaMrhARkUIxvDC8OJQQAp9lFOJ3W46jvtmEMWE++PC+qQjw0jq7aUREpBDcKk0OJUkSFsdHYNMD0xDio8OZ4hrcs/EHm/cuERER9RWGF+oz4/V+eP+XSQjwdENmXhV+80E6zOYBNbBHRET9AMML9amRoT7YuCoROo0Ke86W4R+7s5zdJCIiGmAYXqjPxUb54/klEwAAL6dkYc+ZUie3iIiIBhKGF7omlk2OxIqkaAgBPPZJBoqqG5zdJCIiGiAYXuiaeXbBOMRE+KKq3oiHPzyKFpPZ2U0iIqIBgOGFrhmdRo1/3TkJ3joNfrhwGet2cf0LERH1HsMLXVNDgr3wwlLL+pf132ZjX1aZk1tERERKx/BC19zCWD3uTGxd//JxBkoNjc5uEhERKZjG2Q0g1/DcgnE4eukyzhTX4NGPM/DuL5KgVtl/izUANDSb8O3ZUuzLLkdxdSN83DWIj/LH4vgI+Hvyyb1ERK6Mrwcgh8kurcXCf+1HfbMJv7phGJ6aNwaSZBtgahqNePfARfx7Xy4qWl/6eCU/Dzf8ZXEMFsbqHdVsIiJygJ7cvznyQg4zIsQbzy+JwWMfZ+K11BxU1DXj97eMRYCXFqWGRrx/8BLe+i4XhkbLawUi/D0wZ3wYRoZ6o7KuGZ9lFOBcSS0e/vAoymuacO/0oU7uEREROQNHXsjh3k27gGc/PwkhAI1KQrC3DsVXrIMZPsgLD80egYWxemjUbcuyWkxm/H3nWbyWmgNJAl5dMRlzY8Kc0QUiIupjHHmhfu3uaUMwIsQH/73tFE4WGuTgEh/tj19OH4a5MWF218No1Co8NW8MGowmvJN2EU9tPobJgwMwyEfn6C4QEZETceSFnCr/cj0q65qh9/dAsHf3QojRZMbCf32H00UGLIzV4x93xl/jVhIR0bXWk/s3t0qTU0UGeGJipH+3gwsAuKlV+NttEwEAn2cW4nSR4Vo1j4iI+iGGF1KkCZF+mD8xHACwbtc5J7eGiIgcieGFFOux5JEAgK9PlSCvst7JrSEiIkdheCHFGhHigxkjgyEE8P7BS85uDhEROQjDCynaiqTBAIBPDufByLdWExG5BIYXUrTksSEI9taisq4Z32WXO7s5RETkAAwvpGgatQrzYiwLd7cfL3Jya4iIyBEYXkjxrLuOdp4sQXMLp46IiAY6hhdSvClDAhHsrUV1gxGHL1Y6uzlERHSNMbyQ4qlVEm4YNQgAsPdcmZNbQ0RE1xrDCw0IM63h5SzDCxHRQMfwQgPCjJGDIEnAmeIalFzxhmoiIhp4GF5oQAj00mJihB8AIJVTR0REAxrDCw0Y1nUvqVl83gsR0UDG8EIDxvQRwQCA77LLYTYLJ7eGiIiuFYYXGjDiowPgqVWjsq4Zp4sNzm4OERFdIwwvNGBoNSpMHRYEANjPqSMiogGL4YUGlOtbp472X/Weo+p6I3759mEsWv8dGo0mZzSNiIj6iMbZDSDqSzNGWsLLodxKNBpNcHdTo6nFhLvfPIhj+dUAgJyyOozT+zqzmURE1AsceaEBZWSIN0J9dWhqMePIxcsAgHW7suTgAgACXMxLRKRk/TK8LFmyBAEBAVi2bJmzm0IKI0mSPHWUcroU6Zcu47W9522uEcwuRESK1i/DyyOPPIJ33nnH2c0ghZo/wfKW6fcOXMT97xyGWQCL4/QI83V3csuIiKgv9MvwMmvWLPj4+Di7GaRQs0aHYGiwF5pNZpTXNiPUV4c/LYyBJFnOc+SFiEjZ+jy8pKamYsGCBdDr9ZAkCVu3bm13zfr16zFkyBC4u7sjKSkJhw4d6utmkAtTqyT88854RAd6YvggL/x75RT4ebpB1ZpeuOaFiEjZ+jy81NXVITY2FuvXr7d7/uOPP8aaNWvw3HPPIT09HbGxsZgzZw5KS0t/VH1NTU0wGAw2X0QxEX5I/e1spDw+CzGt7zyy4sN3iYiUrc/Dy7x58/Df//3fWLJkid3zL730Eu677z6sWrUK48aNw4YNG+Dp6Yk333zzR9W3du1a+Pn5yV9RUVG9aT4NYG3TRkwvRERK5tA1L83NzThy5AiSk5PbGqBSITk5GWlpaT/qM59++mlUV1fLX3l5eX3VXBpg5PDi3GYQEVEvOfQhdeXl5TCZTAgNDbU5HhoaijNnzsjfJycnIzMzE3V1dYiMjMSmTZswbdo0u5+p0+mg0+muabtpYJDQuuaF6YWISNH65RN2d+3a5ewm0ACkkqz/xvRCRKRkDp02Cg4OhlqtRklJic3xkpIShIWFObIp5IKk1nkjLtglIlI2h4YXrVaLyZMnIyUlRT5mNpuRkpLS4bQQUV+xDrxw2oiISNn6fNqotrYW2dnZ8ve5ubnIyMhAYGAgoqOjsWbNGqxcuRIJCQlITEzEunXrUFdXh1WrVvV1U4hscbcREdGA0Ofh5fDhw5g9e7b8/Zo1awAAK1euxMaNG7F8+XKUlZXh2WefRXFxMeLi4rBjx452i3iJ+po88uLUVhARUW/1eXiZNWtWl3+zXb16NVavXt3XVRN1yrrmhQMvRETK1i/fbUR0Lag4bURENCAwvJDLkJ/z4uR2EBFR7zC8kMvgW6WJiAYGhhdyOXyrNBGRsjG8kMvggl0iooGB4YVchoovZiQiGhAYXshlWNe8mDn0QkSkaAwv5DIkcOiFiGggYHghlyHvNmJ6ISJSNIYXchl8MSMR0cDA8EIug7uNiIgGBoYXchlcsEtENDAwvJDL4FuliYgGBoYXchmcNiIiGhgYXshlSPK/Mb0QESkZwwu5DBVHXoiIBgSGF3Id8oJd5zaDiIh6h+GFXEbbgl2mFyIiJWN4IZchP2GX2YWISNEYXshlWN9txOxCRKRsDC/kMlStf9oFh16IiBSN4YVchjzywuxCRKRoDC/kMvhWaSKigYHhhVwOR16IiJSN4YVcBl8PQEQ0MDC8kMtQydNGRESkZAwv5DKsD6kzc+iFiEjRGF7IZUgSh16IiAYChhdyGXw9ABHRwMDwQi6DrwcgIhoYGF7IZci7jZzcjt4wmwXMDnwttsks0NRiclh9RETdoXF2A4gcpaMFu/uzyvH9+XIMCfbC4rgIaDXdy/RHL13GF5lF8HHXYEVSNEJ83fu4xW2MJjP+vvMs3j9wEQLAiqRo/HbuGLipr83fP4QQeHXveby65zzqmlswNyYMa5dMhJ+n2zWpDwBKDY14Y18OCqoaMGtUCJZNjoTKukWMiOgKDC/kMq6eNhJC4MWvzuC11Bz5mvcPXsK7v0iEr3vnN+lPfsjDk5uPyZ/1TtoFfHj/VIwJ8+20nNFkxj93Z2PbsUL4e2rx8I0jMXPUoE7LCCHw9Obj+PRIvnzsjX25uFxvxP/cHttpWQDILq3BJ4ctZe9IiMKIEO8uy6zblYWXU7Lk77cfL0ZhVSM2PTCty8B0srAar+3NQVWDEbdODMftkyPbFkt3IK+yHne8loai6ka5vrScCrx0R2yXZYnI9XDaiFzG1W+V/vRIvhxc5k8Mh7+nGzLzqvCHLSc6/Zzj+dV4estxCAEkjw3FmDAfXK434sH30tFo7HiKxWwWWPNJJv6RkoXzZXU4cvEyVr11CDtOFHVa386TJfj0SD7UKgkv/zQO/7orHirJ0v6uyu49V4Zb/rEfr6fm4PXUHNzyj33Ye66sy/79Y7cluDxz6zhsefA6+LprkJFXhVf2nO+07MGcCix95Xt8nlmI1HNl+O2nx/DnL091WsZkFlj9QTqKqhsxbJAXHpw1HBqVhC1HC/DxD3mdliUi18TwQi5D/gu8EKhuMGLtV2cAAI/fNArr75qEN++ZArVKwueZhdiXZf8GL4TA77ceh8ksMH9CON74+WR8dP9UhPjokFNeh//bl2O3HAB8mp6PLzIL4aaW8MKSCbhtUiTMAnjyP8dRYmi0W6bFZMbz2y03/1/PHI5FcRG4daIev541HADwt51nYepgDUypoRGPfnQUzS1mTBsWhKnDAtHcYsYjHx1FaY39+gDgv7edghDAojg9fjF9KOKjA/CXxTEAgNdTz+NyXbPdctUNRjz6cQaaWsyYMTIYv/nJCEgS8NZ3F7DrVEmH9X1yOA+Z+dXwcdfgg19OxW/njsFv544GAKz96gxqm1o6LEtEronhhVzGlY95ef/gRVTWNWNEiDceaA0Ck6IDcPfUwQCA//3mHISdbUn7sspxLL8aHm5q/GnReEiSBH9PLX53y1gAwL/356Khuf3oS0OzCX/bYQlL/3XzaNyVFI0Xb5uAiZF+qG4w4p+7s9qVAYAdJ4uRV9mAIC8tHpo9Qj7+wMzh8Pd0Q05ZHXacKLZb9u87z+JyvRHj9b54+95EvHNvEsbrfVFVb8Rfvzprt8yx/CoczK2ERiXhyblj5OMLY/UYF+6LumYTNn5/wW7ZN1JzUFTdiCFBnnjt7sl4/ObRuH/GMADAH784CaPJ3K6MySywYa9lNOeRG0cizM+ybugX04dh+CAvVDcY8W7aRbv1EZHrYnghl2FdO9HcYsY731tuiA/MHG6zhuPBWcOh1aiQfqkKJwsN7T7j3/tzAQB3JkYj2FsnH791YjiiAj1wud6IrRkF7cr9Jz0f5bXNiAr0wL3ThwIA3NQqOfR88kO+3dEXaztXTB0MD61aPu7j7iYHrY9+uNSuXP7lemw5amnHXxbHQKtRQatR4fklEwAAWzMKkFdZ366cNSgsiNVD7+8hH5ckCb+aaQkimw7ntRvtqW4wyqHmqXlj4am1LKd7NHkUgry0yL/cgM8zCtvVl3K6BBcr6uHv6Ya7kqLl42qVhAdmWkLlewcuOnSHFRH1fwwv5DKss0bfn69AsaERQV5aLIgNt7kmxNcdc8aHAbBMZ1yprKYJ+7PLAQB3Txtsc06jVuFnSZZjW9Jtw4sQAm+33thXXTfUJixNHRaEyYMD0Gwy2yzIBYDCqgYculAJSQLuTIxq1587EizH9meXI/+ybRD58NAltJgFrhsehEnRAfLxuCh/zBgZDJNZ4N0DtiMaTS0m7DhpGcX56ZT29c0ZHwY/DzcUVjfKPwerbceKUNvUghEh3rh5XKh83EOrlsPa1fUBkIPeHQlRcuCxWhCrh49Og4KqBhzMrWxXlohcF8MLuQzryMvuM6UAgHkTwqDTqNtdd0dCJABg69ECmwW4244VwmQWiI3yx9Bgr3blFsbpIUnAoQuVNmHiVJEBWaW10GpUWNb62Vda3hoUNh3Os5mq2nbMshh3ypBAhPt5tCsXFeiJqcMCIQTw5bG2hbtms8DWo5ZRjitHM6ysIzab0wvQcsVUTuq5ctQ0tiDUV4cpQwLblXN3U2NhrB4A8NlR24C25agleN1uZ3vz8ilR0KgkZORVIbu0Rj5e29SClNOW/xbWz726vvkTLeHyP+n57c4TketieCGXcfWG21sntr9hAsD1w4MR4e8BQ2OLHHQAYGvrtMciOzdaAAj388DUoUEAgM8z26ZIrNMlN44JsbsFe/6EcHhq1bhQUY/DFy/Lx784Zim3oIP6AGBhbITl2ivqO3zxMgqqGuCt0yB5bGi7MrPHhCDIS4vy2iakXrEw2foZt07Ud/h8lUVxlrbsPFksB7u8ynr8cOEyJMkS4K4W7K3DrNEhAIBPj7SFnl2nStDUYsawYC+M19vfYn7bZEvY++p4Uac7uYjItTC8kMu48nEhIT72RxcAQKWS5Juwdd3IhfI6ZORVQSUBt1411XQl6839i0zLSIjZLORQYG90AQC8dBrcMsHymZvT2+o7ll8NtUrCvJiwDuubGxMGjUrCyUIDzpfV2rR5XkwY3N3ajyy5qVVYFGcJPdapqvrmFnzTuiOos7A0KToAEf4eqGs2ycFua2t91w0PsjtCBADLJkfI11rXy1gD3oJYfYfPcplspz4iIoYXchmqK26Q8yeGQ93J01sXt97cvz1biqr6ZnltxvUjghHi0/GTdK1h4nSRJUwcvngZhdWN8NZpMHtMSIfllsZb6tt2rBCNRpMceK4bHmSzMPhqgV5aTB8ZDMAyctLUYsK21hGbJa2fac9trWFi1ylL/1JOl6LBaEJUoAdiI/06LKdSSXJ4+zyjEEIIOSwtiW8/JWY1e0wI/DzcUGxoxIGcClyua0Zq6/NmOgtLKpUkn//MzkJoInJNDC/kMuqb254X0tkNEwBGh/lgbLgvjCaB7ceL5dGFzgIBAPh7toWJLzOL5LUgc8bbHwWxmjosCOF+7jA0tmDPmVI5EHQ0WnMl6zVfZBbi65MlMDS2INzPHUnDgjosM17vhzFhPmg2mfHFsSJ8ltE2OtTVE22t9e0+W4rUrHLklNfBw02NuZ2MEOk0betXNqcX4MvjRWgxC4zX+3b5xF/raNaeM2WobjB2ei0RuQaGF3IZO0+2PSgtPsq/y+sXt940f7flOC5U1MPDTS3vROqMdS3Nxz9ckqePrCMdHblyquqJT48hp7wO7m4qzJvQ8RSV1U3jQqHVqHC+rA5//PwkAGDZ5MhOR5as1wDAhm/PY+85y5SMdcSpM+PCfTF8kBeaW8y4d+MPACxTVN66zt82Yh1d2nGiCO+mXQDQdRgEgDFhPhgV6o1mkxn/2p2FwqqGLssQ0cDG8EIuY9X1QwBYdhN15305S+Ij4O7W9r/IHQmR8OriBg0AN48PhY9Og8LqRtQ2tSA60FNeyNuZ5QlRUKsk+Ymyi2IjugwEgOWZLze2TklV1DVDJbUFk84sjNNDq1ahoKoBRpPApGh/jAz16bKcJEnymhnr+pWfJrbf1XS1yYMDMDTYC3XNJpwrsey+srfA1159SydZ+vPGvlzc8Vpal2WIaGBjeCGX8dhNo/DmPQl4cenEbl0f4uuO5xaMh1olYVSoNx5JHtWtcr7ubngkeaT8/R/mj+3W25GHDfKWn6Ib6qvD4zd3rz4AePzm0XLQuf+G4Rgc1H4r99VCfNzx6E2WdrqpJfx+/rhu13fP9UMwKtQy3XPbpEgkDrW/+PlKkiThmVvHyt8/cMOwTtcPXWnV9UOgb336bv5ljrwQuTpJ2HsGuoIZDAb4+fmhuroavr6dv+GXqDsajSa4qVVdTsNcSQiBtPMV8NJpENuNKaor5VXWY5CPrtM1MvYUVjWgqLoBk6IDuv0mZiEEjhdUw99Di+ggzx7V12g0Ibe8DqNDfboVzqwy86pQWdeMmaMG9ahcQVUDrn9xN9QqCedfuKVHbSWi/q8n9++ux6SJXFxPQwRgGWW4bkTwj6ovKrBnIcJK7+9h80j/7pAkCRMj/X9Ufe5uaowN7/lfEHoa5qw0PQg6RDSwcdqIiBSh7aXgA2qwmIh+BIYXIlKGK94KTkSujeGFiBRBak0vHHghIoYXIlKEbq5BJiIXwPBCREREisLwQkSKcOXACxftErk2hhciUoQrn13D7ELk2hheiEgRbEZenNYKIuoPGF6ISBGuXLDLaSMi18bwQkSKIIHbjYjIguGFiBSH4y5Ero3hhYiUwWbayHnNICLnY3ghIkWwWfPCsRcil8bwQkSKYPucF6c1g4j6AYYXIlIEie8HIKJWDC9ERESkKAwvRKQInDYiIiuGFyJSBC7YJSIrhhciUoQrH1LHkRci18bwQkSKYDvyQkSujOGFiBSH7zYicm39MrwsWbIEAQEBWLZsmbObQkRERP1MvwwvjzzyCN555x1nN4OI+hFOGxGRVb8ML7NmzYKPj4+zm0FE/QgX7BKRVY/DS2pqKhYsWAC9Xg9JkrB169Z216xfvx5DhgyBu7s7kpKScOjQob5oKxG5MJsH7DK8ELm0HoeXuro6xMbGYv369XbPf/zxx1izZg2ee+45pKenIzY2FnPmzEFpaal8TVxcHGJiYtp9FRYW/vieENGAZptdmF6IXJmmpwXmzZuHefPmdXj+pZdewn333YdVq1YBADZs2IBt27bhzTffxFNPPQUAyMjI+HGttaOpqQlNTU3y9waDoc8+m4j6J04bEbm2Pl3z0tzcjCNHjiA5ObmtApUKycnJSEtL68uqZGvXroWfn5/8FRUVdU3qISLn4osZiciqT8NLeXk5TCYTQkNDbY6HhoaiuLi425+TnJyM22+/Hdu3b0dkZGSnwefpp59GdXW1/JWXl/ej209E/ReXvBCRVY+njRxh165d3b5Wp9NBp9Ndw9YQUX9gs1Wa80ZELq1PR16Cg4OhVqtRUlJic7ykpARhYWF9WRURuZgrp40YXYhcW5+GF61Wi8mTJyMlJUU+ZjabkZKSgmnTpvVlVUTkwjjwQuTaejxtVFtbi+zsbPn73NxcZGRkIDAwENHR0VizZg1WrlyJhIQEJCYmYt26dairq5N3HxERERH1Ro/Dy+HDhzF79mz5+zVr1gAAVq5ciY0bN2L58uUoKyvDs88+i+LiYsTFxWHHjh3tFvESEfWUJFlGXficFyLX1uPwMmvWrC4Xy61evRqrV6/+0Y0iIrJHQut6F2YXIpfWL99tRERkj3XRLrMLkWtjeCEixeGCXSLXxvBCRIph3SzNNS9Ero3hhYgUg28IICKA4YWIFERqHXvhtBGRa2N4ISLlaB15YXYhcm0ML0SkGPKaFw69ELk0hhciUhxmFyLXxvBCRIrBBbtEBDC8EJGCSGB6ISKGFyJSEOvIC6eNiFwbwwsRKQYfUkdEAMMLESkQR16IXBvDCxEpBl/MSEQAwwsRKQiX6xIRwPBCREoiL9jl2AuRK2N4ISLFaFuwS0SujOGFiBRDXvPC9ELk0hheiEiBmF6IXBnDCxEpBh9SR0QAwwsRKQh3GxERwPBCRArC57wQEcDwQkQKIu82YnohcmkML0SkOHy3EZFrY3ghIsXggl0iAhheiEhR+JwXImJ4ISIFkbjdiIjA8EJECtL2egAOvRC5MoYXIlIcThsRuTaGFyJSDE4bERHA8EJECiJxwS4RgeGFiBRE3irNNS9ELo3hhYgUg7NGRAQwvBCRgsjvNuLAC5FLY3ghIsVhdiFybQwvRKQ4gkMvRC6N4YWIFKNtwS4RuTKGFyJSDD7nhYgAhhciUhA+54WIAIYXIlIkphciV8bwQkSKIa95YXYhcmkML0SkGG1vlSYiV8bwQkSKwYfUERHA8EJECsLNRkQEMLwQkQLxIXVEro3hhYiUgw+pIyIwvBCRgsgLdpleiFwawwsRKYa8YLebYy9CCOzLKsO/9+ciM6/qGrasTV5lPb7ILMTZ4hqH1NfUYkL6pcvIv1zvkPoAoLKuGTWNRofVR3Q1jbMbQETUXfKC3W5kF6PJjEc/ysC240XysV/NHIan5o6RQ1BnymqakF1ai7HhPvD31HarfW9/fwF//vIUTGZLA1dOG4znFoyHStV1ffmX63H4wmUMG+SFiZH+3arvyMVKPPh+OkoMTQCAu5Ki8eeF46FRd/330uzSWnx9qhjB3josjNXD3U3dZZnKumas+SQD354tg1ol4a7EaDxz6zhoNV3Xd6KgGpsO50GlkrDquqGIDvLsuoNEHWB4ISLFsGaOi5X1uK6La9duP4Ntx4vgppaQODQQ32VX4LW9ORgxyBu3J0R1WvbdtAv4y5en0Wwyw8NNjbVLJ2BxfESnZb46XoTnPj8JABgR4o3zZbV4O+0iQv3c8eCsEZ2W3Xq0AE98mgmjyRJ67kyMxvOLYzoNPbnldVj11g8wNLbAx12D2qYWfHDwEnx0Gjx9y9hO69t+vAgPf3gULa0h69/7cvHBfUkI8tZ1WKbRaMKqjT/II1gms8C7By6ixWzG2qUTO63vq+NFWP3hUTnUZZfW4t1fJHVahqgznDYiIsUYGeoDAPh/X59FqaGxw+uO51fjre9zAQD/umsS3v/lVDx+0ygAwJ+/OIXLdc0dlt19pgTPfHYSzSYz/D3d0GA0Yc0nGTiYU9FhmeoGI5757AQA4JfTh+Kbx27A84snAADWfZOF3PK6DsseuViJNZ9kwGgSGBHiDZUEfHjoEjaknu+wjBACz2w9AUNjCyZF++PQ75Lxj5/GAwBeS81BRidTZFklNXjkI0twmTw4AMHeOpwtqcGjH2d0uovr9dQcZOZVwd/TDV89MgOv3z0ZkgR8eCgPKadLOixXUNWA/9qUKQcXAKhtaunweqLuYHghIsV4YfEEjA71QXltM17ccabD615OyYIQwKI4PeaMDwMAPDh7BMaG+6KmqQVvfZdrt1xDswlP/ec4AODn0wYj/Q83YUl8BMwCeHrLcTS3mO2WezftAsprmzF8kBeemDsakiThzsQo3DBqEJpNZryyJ9tuOZNZ4HebT8AsgIWxenz96A1Yu9QSel7elYW8SvvrWPZnl2N/djm0GhXWLY+Hh1aNBbF6LJ1kGR36n51nO/zZ/GHrCRhNArNGD8Inv5qGj+5Pglajwr6scmw/Xmy3THltE1751tKHPy+KwdhwX9w8Pgy/nD4UAPD3nWdhNtsPPi9sO426ZhMSBgdgw88mA+CCa+o9hhciUgw/Tzf8bZllimJzegHOFBvaXXO+rBa7TpdAkoCHbxwpH1erJDz8E8v0zdtpF9HUYmpX9v2DF1Fa04SoQA/87paxUKkk/HHBeAR7a5FTVocvMgvblWlqMeHN7y4AAH7zk5HQaSxrRyRJwmPJlvq3HC2wO1K082QxzpbUwN/TDX9aaFkbc0dCFKYNC0JTixmvdTD68n/7LOHrrsRom7UjjyWPgkYlYX92OU4Xtf/ZHL5QiYO5ldCqVXh+yQSoVRJGhPjggRuGAQD+tSfb7ujL+wcuodFoxsRIPyyYGC4ff2j2CPjoNDhTXIO958ralcspq8X2E5Y1R39ZHAN1N9b+EHUHwwsRKUpslD9umWAZTXn7+wvtzm9OzwcAzB4dguGDvG3O3Tw+DGG+7qhuMCLldKnNObNZ4K3WEPLgrBHyAlY/Tzfc2zrC8Ma+nHY392/PlqGyrhmhvjrcesWNHQDiowMwKdofLWaBzzLaB5839uUAAH4+bQgCvCyLgiVJkkPXpsP57aa4LlXUY++5MkgSsOr6ITbnogI9cfP4UADABwcvtavv3/stoWfppAhE+HvIx++dPhQebmqcLjIg7arpsRaTGe8dvAgA+OWMYTaLnf09tfL6oQ8Pta/v3QMXIQRw45gQjA33lY9z4IV6i+GFiBTnnussYWLL0QJUN7Rt2RWiLSQssbPAVq2S5KmVT4/k25w7dKESBVUN8NFp2pVdkTQY7m4qnCmuwbH8aptzW48WAAAWxUXY3eVz2+RIAMB/0vNtgk92aS2OXqqCRiXh7qmDbcpMHRaIseG+aGox48tjtqHn80xLfdcPD8bgIK929d2VaPmsrRkFaDS2jS5V1TdjV+valHuuCj3+nlos6eDnkpZTgbKaJgR4umFeTFi7+u5MtISXlDOlKKtpko+3mMz4ItMy6vKz1v617RZjfKHeYXghIsWZMiQAo0N90Gg020zlHLl4GfmXG+Ct0yB5bKjdstYwsfdcGSpq22621hByy4TwdtuG/TzccPM4y417S+t1AGxGcBbH2d+NdOsEPbQaS/A5dcVUzuet7b5h1CAM8rHd5SNJEm5rDRObr6hPCIGtreFsUZzebn3XDQ+C3s8dNY0tNqNL244XwWgSGBvuizFhvu3KWevbeaIY9c1tC2qtYXD+xHC42QlnI0N9EBvlD5NZ2ASttJwKlNdaQs/0kcGt/bLbZKIeY3ghIsWRJAm3J1hCyKbDefLxrRmWG/2c8WHw0Np/bsnwQd6YEOEHk1ngqxOWBaqNRpP8PBjrCMTVrKMxX2QWwmiyLNzdcaIIzSYzRof6YGy4j91yfp5uuHFMCIC2wGIZIbKO2NgPIQtj9VBJwNFLVbjQulvpdFENsktrodWoMMfOKAgAqFSS3AfrFBrQFs6WxNuvb1J0AKICPVDXbMI3pywjNI1GE3a0/owWdRDOAGBxax+2XjE11lno4bgL9RbDCxEp0pL4CLipJWTmV+NMsQEtJjO+Om690dq/QVtZ16ZYRwr2nClFTWML9H7uSBwSaLfM9JHBCPLSoqKuGfuzywFAnhZZGKfv9MF3C2It7fkyswhCCGTmV+NiRT083NQdjhCF+Lrj+hGWEQtrEPiitb2zRw+Cr7tbh/UtibcEu2/PlaG8tgkFVQ344cJlSFJbW64mSZI8evR5a33fni1FbVMLIvw9MDk6oMP6bp2oh1olITOvCrnldTahZ2FsW+jhyAv1FYYXIlKkIG+dfOP/5Id8pOVUoKKuGYFeWlw3PKjTsvNbw8vB3EqUGBrxn/TWUZD4iA4fDOemVsk3/q1HC1Be24Tvz1tCzIKJnYeln4wJgZdWjYKqBqRfuiyPgtw0LhReuo6fFWoNE59lFECItmmZjgKI1YgQb8RGWkaXvsgslKfWEocEItzPo8Ny1tC391wZLtc1y+Hs1onhnT4wb5CPTg5aW48WYNfpEtQ2WcJgwuD2oYdLXqi3GF6ISLHuaN3psuVoPra0BpC5MWFdPh4/MsATk6L9IQTwbtpFfHvWsjbktg6mjKysN/evT5bg0yP5MAtgYqRfl4+6d3dT4+bW581sTi+Qw0RHU1RWc2LCoNOokFNeh/cOXEReZQM8tWr8pHUaqjPWaa7N6QXyyE1nUz8AMCLEB+P1vmgxC2w6koeUM5bpo67CEtA2dfRZRoG86HfJJNswKIFDL9Q3GF6ISLFmjAxGqK8Ol+uN8sLWRd240QKWqQ7A8myTFrNAbKQfRoTYX7diFRflj6HBXmgwmvDiV5aH5HUVCKwWxFpGe94/eAkVdc0I9tZhRutoRUe8dRrcNM4yuvTMZ5ZXD8wZHwZPbddvdlkQq4dGJeF4QTVOFxmg1ajs7ha6mnW054XtZ9BoNGP4IC+M17df4Hu1m8eHwd1NhQsV9fj2rOWZL7dNirS9qDW7dPfFmkQdYXghIsXSqFVYNrntBhkX5Y/EofbXrFxtWUIkAr3aXrj46y7ePwRY1oX8euZw+Xs/Dzcsn9L5e5Ksbhg5CKND28LRAzOHdesFig9cUR8A3DdjWLfqC/LW2YzQ3JUYLT9LpjPLE6MQ4Nm2nuZXM4d360WW3jqNTViZMz4Uw656zo4Vp42ot/hiRiJStF/NHI4L5fUwNBrxp4Xju3WjBQBfdzesv2sS/rk7C9cND8Kc8fYXzl5t2eRInCupwdG8KjyWPArenaxZuZJGrcI/74rHi1+dwZAgL6y8bki3ysVE+OH5JTH45HA+7kiIxLhujIJY/WnReLipVfDSqfFfc0Z3q4yvuxv+ddckvPTNOUweHIBlV4+edOIP88fBZBZoMJrwzK3j2p3npBH1FUl09iYuBTIYDPDz80N1dTV8fbv/PzkREV1b354txT1v/YBx4b7Y/sgMZzeH+pme3L85bURERA7R3VExoq4wvBARkUMNqOF+cgqGFyIicgjruMsAW61ATsDwQkREDsFZI+orDC9ERESkKAwvRETkEHzCLvWVfhdeqqqqkJCQgLi4OMTExOCNN95wdpOIiKgPcckL9Va/e0idj48PUlNT4enpibq6OsTExGDp0qUICur8RWtERNS/cc0L9ZV+N/KiVqvh6Wl5yVlTUxOEEFyZTkQ0AMi7jbhZmnqpx+ElNTUVCxYsgF6vhyRJ2Lp1a7tr1q9fjyFDhsDd3R1JSUk4dOhQj+qoqqpCbGwsIiMj8cQTTyA4uPOXlxERkXLw76PUWz0OL3V1dYiNjcX69evtnv/444+xZs0aPPfcc0hPT0dsbCzmzJmD0tJS+RrreparvwoLLa9t9/f3R2ZmJnJzc/HBBx+gpKTkR3aPiIj6DU4bUR/p8ZqXefPmYd68eR2ef+mll3Dfffdh1apVAIANGzZg27ZtePPNN/HUU08BADIyMrpVV2hoKGJjY7Fv3z4sW7bM7jVNTU1oamqSvzcYDN3sCREROQMHXqi3+nTNS3NzM44cOYLk5OS2ClQqJCcnIy0trVufUVJSgpqaGgBAdXU1UlNTMXp0x29DXbt2Lfz8/OSvqKjuvZ6eiIgci1ulqa/0aXgpLy+HyWRCaKjtq+VDQ0NRXFzcrc+4ePEiZsyYgdjYWMyYMQO/+c1vMGHChA6vf/rpp1FdXS1/5eXl9aoPRER0bVh3G3ETBvVWv9sqnZiY2O1pJQDQ6XTQ6XTXrkFERETUr/TpyEtwcDDUanW7BbYlJSUICwvry6qIiEhh2rZKE/VOn4YXrVaLyZMnIyUlRT5mNpuRkpKCadOm9WVVRESkVEwv1Es9njaqra1Fdna2/H1ubi4yMjIQGBiI6OhorFmzBitXrkRCQgISExOxbt061NXVybuPiIjINUl8xC71kR6Hl8OHD2P27Nny92vWrAEArFy5Ehs3bsTy5ctRVlaGZ599FsXFxYiLi8OOHTvaLeIlIiLXxIEX6q0eh5dZs2Z1uVJ89erVWL169Y9uFBERDTwceKG+0u/ebURERAOTvGCXW6WplxheiIiISFEYXoiIyCHkh9Q5txk0ADC8EBGRQ3HWiHqL4YWIiByEK3apbzC8EBGRQ7RNG3HohXqH4YWIiIgUheGFiIgcom2rtFObQQMAwwsRETkUwwv1FsMLERE5BN9tRH2F4YWIiIgUheGFiIgcguMu1FcYXoiIyCHkrdJc9EK9xPBCREREisLwQkREDiG1Thxx3IV6i+GFiIgcirNG1FsML0RE5BDcKU19heGFiIgciu82ot5ieCEiIiJFYXghIiKHaNsq7dx2kPIxvBAREZGiMLwQEZFDcKs09RWGFyIicihOG1FvMbwQEZFDcKs09RWGFyIicoi28MKhF+odhhciIiJSFIYXIiJyCHnBLgdeqJcYXoiIiEhRGF6IiMgh5IfUObcZNAAwvBARkUNY1+sKzhtRLzG8EBERkaIwvBARkUNw2oj6CsMLERERKQrDCxEROQi3SlPfYHghIiIiRWF4ISIih5DXvHDohXqJ4YWIiBxC3irt1FbQQMDwQkRERIrC8EJERA4hca809RGGFyIiIlIUhhciInIIrnmhvsLwQkREDsHdRtRXGF6IiIhIURheiIjIISTrE3ad3A5SPoYXIiIiUhSGFyIicoi2NS/ObQcpH8MLERERKQrDCxEROZTgqhfqJYYXIiJyCE4bUV9heCEiIiJFYXghIiKHsL7biAMv1FsML0RERKQoDC9EROQQ1ncbceiFeovhhYiIHEKSur6GqDsYXoiIyKG4VZp6i+GFiIgcQn63EbML9RLDCxERESkKwwsRETmE/JA65zaDBgCGFyIiIlIUhhciInII62YjwUUv1EsML0RE5BjcKk19hOGFiIgciuMu1FsML0RE5BDcKk19heGFiIiIFIXhhYiIHIKvB6C+wvBCREQOwexCfYXhhYiIHI7bpak3GF6IiMghJM4bUR9heCEiIofjwAv1BsMLERE5xJXjLswu1BsaZzfAniFDhsDX1xcqlQoBAQHYs2ePs5tERERE/US/DC8A8P3338Pb29vZzSAioj5y5ZIXy4JdroGhH4fTRkRE5BASwwr1kR6Hl9TUVCxYsAB6vR6SJGHr1q3trlm/fj2GDBkCd3d3JCUl4dChQz2qQ5IkzJw5E1OmTMH777/f0yYSEVE/xzUv1Bs9njaqq6tDbGws7r33XixdurTd+Y8//hhr1qzBhg0bkJSUhHXr1mHOnDk4e/YsQkJCAABxcXFoaWlpV/brr7+GXq/H/v37ERERgaKiIiQnJ2PChAmYOHGi3fY0NTWhqalJ/t5gMPS0S0RE5Ag200bOawYpX4/Dy7x58zBv3rwOz7/00ku47777sGrVKgDAhg0bsG3bNrz55pt46qmnAAAZGRmd1hEREQEACA8Pxy233IL09PQOw8vatWvxpz/9qafdICIiIoXq0zUvzc3NOHLkCJKTk9sqUKmQnJyMtLS0bn1GXV0dampqAAC1tbXYvXs3xo8f3+H1Tz/9NKqrq+WvvLy83nWCiIiuCZsFu5w4ol7o091G5eXlMJlMCA0NtTkeGhqKM2fOdOszSkpKsGTJEgCAyWTCfffdhylTpnR4vU6ng06n+/GNJiIih+ByXeor/W6r9LBhw5CZmensZhAR0TXENS/UG306bRQcHAy1Wo2SkhKb4yUlJQgLC+vLqoiISGH4biPqK30aXrRaLSZPnoyUlBT5mNlsRkpKCqZNm9aXVREREZGL6vG0UW1tLbKzs+Xvc3NzkZGRgcDAQERHR2PNmjVYuXIlEhISkJiYiHXr1qGurk7efURERK7J5t1GnDaiXuhxeDl8+DBmz54tf79mzRoAwMqVK7Fx40YsX74cZWVlePbZZ1FcXIy4uDjs2LGj3SJeIiIioh+jx+Fl1qxZre+k6Njq1auxevXqH90oIiIaeLhVmvoK321EREQOwXcbUV9heCEiIofjmhfqDYYXIiJyCO6Upr7C8EJERA7HgRfqDYYXIiJyuK42flD/VN1gRKPR5Oxm9L/XAxAR0cDkyGmjgqoG/Gt3Fi6U1+P6EUH45YxhcHdTX7P6SgyN+OuOM8jMq0JMhB9+O3cMIvw9rll95bVNePazE0g9V47IAA88fctYzBw16JrVV1bThP/alIm958qgVavw82mD8Ydbx12z+rrC8EJERA6hliToNCo0tZjx8q6sbt38ahqNeO/AJeRfrses0SFIHhvS5WsG8irrsWzD9ygxNAEA0nIqkJZTgY2rEuGm7nzCwWwW+PpUCU4VViMu2h+zR3ddX2lNI+54LQ0XK+oBAOfL6nAgpwKfPTQdYX7uXfbxfFktDuVWYmiwF5KGBnZZn6HRiDtfP4Cs0loAwJniGty78Qe8ec+UbgWYRqMJRy9VIdhbi5GhPl1e39Bswr0bf8DxgmoAQLPJDF8Pty7LXUsML0RE5BAatQpPzRuDP31xCv+3PxczRg3q9GZbamjE8tcPILe8DgDw/sFLuPf6oXh2Qcehx2wWePijoygxNGFkiDd+mhiNl74+i++yK/DKnvN4JHlkh2WNJjNWf5COnSfb3s9326RI/H3ZRKhUHQeKZ7eexMWKekQFeuCpuWOxbtc5ZJXW4rf/OYa3V03pNIxs/C4Xf/7yFMyts2gLYvV46Y7YTkPWX786g6zSWoT66rBueTw+PHQJn2cW4olNmfjmsZnw8+w4WBzLr8L97xxBsaFRru9/bp8InabjUalXvs3G8YJqBHpp8dH9U9FiEgj20XZ4vSNwzQsRETnMquuHYtX1QwAAz312Ai0ms93rhBB4fFMmcsvroPdzx52JUZAk4M3vcvHlscIOP3/L0QIcvVQFb50G7/wiEb+YPhQvLJ0AAFj/bTZKWm/a9vxrdzZ2niyBVqPCLRPCoFZJ+E96Pt5Ou9Bhme/Pl2PHyWKoVRLe+HkC5k8Mx6s/mwStRoXUc2X49lxZh2X3Z5Xjj19YgktslD/c1BK+yCzEy7uyOiyTXVqD9w9eAgCsWx6PacOD8LdlEzEs2AulNU349/6cDsuW1TRh1Vs/oNjQiABPN6hVlvr+8uWpDsuUGBrxeqrlM19YEoNRoT4Yp/dFiE/XI0rXEsMLERE51OM3j0aQlxYXKuqx/USx3Wu2HS/Cvqxy6DQqvPOLJKxdOhGrZ48AAPzpi1N2F40KIfBa6nkAwIOzhyPcz7LmZGGsHgmDA9DcYpZvxFe7UF6Hf+2xvLfv78sm4pUVk/Fc6wjPX3ecQVlNk91yr+21fN5didEYE+YLABgR4oOfTx0MAFi/O9tuueYWM578zzEAwJ2J0dj64HV46Y44AMCre88jq6TGbjlr+28aF4ppw4MAAO5uajwxZzQAYOP3F1Db1GK37J+/PIWKumaMCfPBvid/gv/7eQIA4L0Dl/DDhUq7Zd4/cBFNLWZMivbHnPFhdq9xBoYXIiJyKG+dBndPs9zc39yf2+68EAL/ar3p/2rmcIwI8QYA/OYnIxHh74GymiZ8dOhSu3Lfn6/AuZJaeGnV+FlreAAASZLwUGvw2XQ4z27weeXbbJjMAjeMGoRFcREAgJ8lDUZslD8ajWb8n50RjdzyOuw9VwaVBNw3Y5jNuftuGAaNSsLhi5dxzk4Q2Zyej4KqBgzy0eGZW8dCkiQsiNUjeWwoTGaBV749366ModGIrRmWUacHZtrWN2d8GIYGe8HQ2ILtx4valT1bXCOPWP2/O2LhrdNg9pgQLE+IAgC7oz0tJjM+OJQHALh3+tAu1+I4EsMLERE53IqkwVCrJGTkVSGnrNbm3MHcSpwproGXVo17W6eYAECrUeFXrTftdw5cbLfdeuvRAgDAovgI+Lrbrvu4YdQghPu5w9DYgl2nS2zOXa5rxpbWso/c2LYmRqWS8JvW0PPBwUvtQs+XmZYwMH3kIEQHedqcC/V1x+wxIQCA/6Tn25wTQuDN7yyh7Vc3DIOntm35qbX+zzMLUVxtO8X19ckSNLeYMXyQFyZFB9icU6kkLJscCQDYkl6Aq725PxdCAPNiwjBe7ycf/82NI6BWSdifXY4TrQtyrdJyKlBe24RAL22/GnUBGF6IiMgJBvnoMH1EMADgswzbNSzWm++tE/Xw97RdGLokPgLubirklNXhaF6VfLy5xYydJy1TUAsm6tvVp1ZJWBJvGVH5zxHbMPHl8SIYTQJjw30xebBtKPjJmBBE+HugprEFX5+yDT1ftI5kLJgYbrePt02y1Lf1aAFM5ragdaLAgHMltdBpVLhjSpRNmQmRfpgyJAAms5ADldXnrWFpYWyE3VGQxa39S8upQP7levl4Q7MJ21pHY+65bohNmcgAT8yLsQSTT6/6uVhHcObGhHW5S8vR+ldriIjIZSyOt4SMzzIK5FGURqNJvmkuab35X8nH3Q23xFjCwqbDbTfbfVllMDS2YJCPDolDA+3Wd1vryERqVjkq65rl41taR0Zus1OfSiXJx68MPedKanCupBZatQo3dzAqMXtMCPw83FBiaMLBnAr5uHUk5ubxYe1GiCztsLRzc3q+/HOpqG3Cd9nlAICFce3DGQBE+Htg6jBL37881jZ19PWpYtQ2tSAq0ANThrT/2VhHbD7LKEBzi2UBtdFkxo7W9UjzJ9gPZ87E8EJERE5x87gwuLupcKGiHicKDACAXadLUNPUggh/DyTaudECbTfbL48VylM51pv1/AnhUHewrXn4IG+M1/vCZBbyjTm3vA7pl6qgkiwLe+2xhp59WWXybiXrlNENo4Lh18EzT3QaNea2BpsvWttnNJnlEZSl8e3DEgDcMjEcOo0KWaW18rNVtp8ohsksMCHCD0ODveyWAyxbnwHY7MiyjqgsiY+0u+V7+ohgDPLR4XK9EXvOlgIADuRU4HK9EYFeWiR1EAadieGFiIicwkunwY1jQwG0TcFYp4wWx+s7fLbK1GFBCPdzR01jC749W4pGownftE7pLIjtfJTg1om2N3fr1MyMkYMQ4mt/++/gIC8kDA6AWbSNElnDyK12pqhs6mttz44TRTCazNh7tgyVdc0I9tZixshgu2V83d3k0RzraM826xRVF/2bO96yxftEgQEXyutQXN0oj9jYG1kCLM/fWdw6mmNdN2Qd/ZozPgyafjZlBDC8EBGRE1nXp3yZWYhSQ6P8XJQl8ZEdllGpJHnqZMvRAuw6XYLa1tGa+KiADssBwK2t61MO5FSgxNAoh4OlHdzYraxTWJvTC3CiwIDc8jq4u6lw07jQTstNGxaEIC8tLtcb8f35Cmw6Ytm9syguotNQYA0an2cWorCqAYdyLVuZb+liCifIW4frWrdQbztehM1H82EWQMLgAAwO6njExvrzTjlditKaRmxrDWcdredxNoYXIiJymlmjB8Fbp0FhdSN+t+UETGaB+Gh/eXt0Rxa3bmfec6YMG7+7AMCyFqSzJ+ECQFSgJ2Kj/GEWwB+2nkBBVQN83DVd7qa5dYIeWrUKZ4pr8ML20wCA5LGh8NJ1/qB6jVqFua0LYt/+/gJSTlumZW5P6DicAa0jQa1TOY99nAGzAOKi/BEZ4NlpOaAtoH2eUYhPW9cF3ZEQ1VkRjNP7YkyYD5pNZqz5OBOGxhaE+bojaVhQl/U5A8MLERE5jbubWt7tYt3C3NWNFgDGhvtiXLgvmk1mHL54GUDHa0iuZp0isU41LYzVd/nSRj9PNySPs2x9TmtdfGtdC9Nlfa3t2n2mFC1mgdhIP/mBdh1RqyR5tOdg66hLd+ubMz4MWo0KZ0tqkFNeBy+tGrd0YwTl9taf+/7WaabF8REdrh9yNoYXIiJyqodvHAl3N8vtaGSIt7zbpitPzRsj//vtkyO79ZJBwPJE26hAy9N3vXUa/OYnHb/v6EprbhoFrcbSzsQhgZjVzbc4JwwOwOzRbdc+OXdMJ1e3+eX0YQjysmwVjw70xO3dDC/+nlrcf8VD835z40h4dzFCBAArktp+LoFeWtx/w7AuSjiPJK5+yo/CGQwG+Pn5obq6Gr6+nSdbIiLqH04VGvD9+XIsiNUjtIOFs/YcyKnApcr6bo2eXKmgqgHbjhVi5qgQjA7rXugBgBMF1cjMr8KCWL3dbc4dMTQaselwPkaH+mB6Bwt17blQXod92eW4eVxoj34uLSYzNh3Jh5tahaXxEV1Op1mVGBqx82QxZowc1OmupmuhJ/dvhhciIiJyup7cvzltRERERIrC8EJERESKwvBCREREisLwQkRERIrC8EJERESKwvBCREREisLwQkRERIrC8EJERESKwvBCREREisLwQkRERIrC8EJERESKwvBCREREisLwQkRERIqicXYD+pr1JdkGg8HJLSEiIqLust63rffxzgy48FJTUwMAiIqKcnJLiIiIqKdqamrg5+fX6TWS6E7EURCz2YzCwkL4+PhAkqQ+/WyDwYCoqCjk5eXB19e3Tz+7P2D/lG+g93Gg9w8Y+H1k/5TvWvVRCIGamhro9XqoVJ2vahlwIy8qlQqRkZHXtA5fX98B+4cSYP8GgoHex4HeP2Dg95H9U75r0ceuRlysuGCXiIiIFIXhhYiIiBSF4aUHdDodnnvuOeh0Omc35Zpg/5RvoPdxoPcPGPh9ZP+Urz/0ccAt2CUiIqKBjSMvREREpCgML0RERKQoDC9ERESkKAwvREREpCgML920fv16DBkyBO7u7khKSsKhQ4ec3aRuWbt2LaZMmQIfHx+EhIRg8eLFOHv2rM01jY2NeOihhxAUFARvb2/cdtttKCkpsbnm0qVLmD9/Pjw9PRESEoInnngCLS0tjuxKt7z44ouQJAmPPvqofGwg9K+goAA/+9nPEBQUBA8PD0yYMAGHDx+Wzwsh8OyzzyI8PBweHh5ITk5GVlaWzWdUVlZixYoV8PX1hb+/P37xi1+gtrbW0V1px2Qy4ZlnnsHQoUPh4eGB4cOH4y9/+YvN+02U1r/U1FQsWLAAer0ekiRh69atNuf7qj/Hjh3DjBkz4O7ujqioKPztb3+71l0D0Hn/jEYjnnzySUyYMAFeXl7Q6/X4+c9/jsLCQpvPUGr/rvbAAw9AkiSsW7fO5nh/7h/QvT6ePn0aCxcuhJ+fH7y8vDBlyhRcunRJPu/U362CuvTRRx8JrVYr3nzzTXHy5Elx3333CX9/f1FSUuLspnVpzpw54q233hInTpwQGRkZ4pZbbhHR0dGitrZWvuaBBx4QUVFRIiUlRRw+fFhMnTpVXHfddfL5lpYWERMTI5KTk8XRo0fF9u3bRXBwsHj66aed0aUOHTp0SAwZMkRMnDhRPPLII/JxpfevsrJSDB48WNxzzz3i4MGDIicnR+zcuVNkZ2fL17z44ovCz89PbN26VWRmZoqFCxeKoUOHioaGBvmauXPnitjYWHHgwAGxb98+MWLECHHnnXc6o0s2nn/+eREUFCS+/PJLkZubKzZt2iS8vb3Fyy+/LF+jtP5t375d/P73vxebN28WAMSWLVtszvdFf6qrq0VoaKhYsWKFOHHihPjwww+Fh4eHeO2115zav6qqKpGcnCw+/vhjcebMGZGWliYSExPF5MmTbT5Dqf270ubNm0VsbKzQ6/Xif//3f23O9ef+CdF1H7Ozs0VgYKB44oknRHp6usjOzhafffaZzX3Pmb9bGV66ITExUTz00EPy9yaTSej1erF27VonturHKS0tFQDE3r17hRCWXzRubm5i06ZN8jWnT58WAERaWpoQwvKHXKVSieLiYvmaV199Vfj6+oqmpibHdqADNTU1YuTIkeKbb74RM2fOlMPLQOjfk08+KaZPn97hebPZLMLCwsTf//53+VhVVZXQ6XTiww8/FEIIcerUKQFA/PDDD/I1X331lZAkSRQUFFy7xnfD/Pnzxb333mtzbOnSpWLFihVCCOX37+obQ1/155VXXhEBAQE2f0affPJJMXr06GvcI1ud3dytDh06JACIixcvCiEGRv/y8/NFRESEOHHihBg8eLBNeFFS/4Sw38fly5eLn/3sZx2WcfbvVk4bdaG5uRlHjhxBcnKyfEylUiE5ORlpaWlObNmPU11dDQAIDAwEABw5cgRGo9Gmf2PGjEF0dLTcv7S0NEyYMAGhoaHyNXPmzIHBYMDJkycd2PqOPfTQQ5g/f75NP4CB0b/PP/8cCQkJuP322xESEoL4+Hi88cYb8vnc3FwUFxfb9NHPzw9JSUk2ffT390dCQoJ8TXJyMlQqFQ4ePOi4zthx3XXXISUlBefOnQMAZGZmYv/+/Zg3bx4A5ffvan3Vn7S0NNxwww3QarXyNXPmzMHZs2dx+fJlB/Wme6qrqyFJEvz9/QEov39msxl33303nnjiCYwfP77d+YHQv23btmHUqFGYM2cOQkJCkJSUZDO15OzfrQwvXSgvL4fJZLL54QNAaGgoiouLndSqH8dsNuPRRx/F9ddfj5iYGABAcXExtFqt/EvF6sr+FRcX2+2/9ZyzffTRR0hPT8fatWvbnRsI/cvJycGrr76KkSNHYufOnfj1r3+Nhx9+GG+//TaAtjZ29me0uLgYISEhNuc1Gg0CAwOd3sennnoKP/3pTzFmzBi4ubkhPj4ejz76KFasWAFA+f27Wl/1p7//ubVqbGzEk08+iTvvvFN+iZ/S+/fXv/4VGo0GDz/8sN3zSu9faWkpamtr8eKLL2Lu3Ln4+uuvsWTJEixduhR79+6V2+jM360D7q3S1LGHHnoIJ06cwP79+53dlD6Tl5eHRx55BN988w3c3d2d3Zxrwmw2IyEhAS+88AIAID4+HidOnMCGDRuwcuVKJ7eu9z755BO8//77+OCDDzB+/HhkZGTg0UcfhV6vHxD9c2VGoxF33HEHhBB49dVXnd2cPnHkyBG8/PLLSE9PhyRJzm7ONWE2mwEAixYtwmOPPQYAiIuLw/fff48NGzZg5syZzmweAI68dCk4OBhqtbrdCuqSkhKEhYU5qVU9t3r1anz55ZfYs2cPIiMj5eNhYWFobm5GVVWVzfVX9i8sLMxu/63nnOnIkSMoLS3FpEmToNFooNFosHfvXvzjH/+ARqNBaGioovsHAOHh4Rg3bpzNsbFjx8qr/q1t7OzPaFhYGEpLS23Ot7S0oLKy0ul9fOKJJ+TRlwkTJuDuu+/GY489Jo+kKb1/V+ur/vT3P7fW4HLx4kV888038qgLoOz+7du3D6WlpYiOjpZ/51y8eBGPP/44hgwZIrdPqf0DLPc9jUbT5e8dZ/5uZXjpglarxeTJk5GSkiIfM5vNSElJwbRp05zYsu4RQmD16tXYsmULdu/ejaFDh9qcnzx5Mtzc3Gz6d/bsWVy6dEnu37Rp03D8+HGb/xmtv4yu/sPtaDfeeCOOHz+OjIwM+SshIQErVqyQ/13J/QOA66+/vt329nPnzmHw4MEAgKFDhyIsLMymjwaDAQcPHrTpY1VVFY4cOSJfs3v3bpjNZiQlJTmgFx2rr6+HSmX7q0itVst/+1N6/67WV/2ZNm0aUlNTYTQa5Wu++eYbjB49GgEBAQ7qjX3W4JKVlYVdu3YhKCjI5ryS+3f33Xfj2LFjNr9z9Ho9nnjiCezcuROAsvsHWO57U6ZM6fT3jtPvHb1a7usiPvroI6HT6cTGjRvFqVOnxP333y/8/f1tVlD3V7/+9a+Fn5+f+Pbbb0VRUZH8VV9fL1/zwAMPiOjoaLF7925x+PBhMW3aNDFt2jT5vHW728033ywyMjLEjh07xKBBg/rNVuKrXbnbSAjl9+/QoUNCo9GI559/XmRlZYn3339feHp6ivfee0++5sUXXxT+/v7is88+E8eOHROLFi2yu/U2Pj5eHDx4UOzfv1+MHDmyX2yVXrlypYiIiJC3Sm/evFkEBweL3/72t/I1SutfTU2NOHr0qDh69KgAIF566SVx9OhRebdNX/SnqqpKhIaGirvvvlucOHFCfPTRR8LT09MhW207619zc7NYuHChiIyMFBkZGTa/d67cYaLU/tlz9W4jIfp3/4Touo+bN28Wbm5u4vXXXxdZWVnin//8p1Cr1WLfvn3yZzjzdyvDSzf985//FNHR0UKr1YrExERx4MABZzepWwDY/XrrrbfkaxoaGsSDDz4oAgIChKenp1iyZIkoKiqy+ZwLFy6IefPmCQ8PDxEcHCwef/xxYTQaHdyb7rk6vAyE/n3xxRciJiZG6HQ6MWbMGPH666/bnDebzeKZZ54RoaGhQqfTiRtvvFGcPXvW5pqKigpx5513Cm9vb+Hr6ytWrVolampqHNkNuwwGg3jkkUdEdHS0cHd3F8OGDRO///3vbW50Suvfnj177P5/t3LlSiFE3/UnMzNTTJ8+Xeh0OhERESFefPFFp/cvNze3w987e/bsUXz/7LEXXvpz/4ToXh///e9/ixEjRgh3d3cRGxsrtm7davMZzvzdKglxxWMsiYiIiPo5rnkhIiIiRWF4ISIiIkVheCEiIiJFYXghIiIiRWF4ISIiIkVheCEiIiJFYXghIiIiRWF4ISIiIkVheCEiIiJFYXghIiIiRWF4ISIiIkVheCEiIiJF+f/83MzepyKk1wAAAABJRU5ErkJggg==",
      "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": 10,
   "id": "swedish-projection",
   "metadata": {},
   "outputs": [],
   "source": [
    "tle_filename = 'tles/TWOLINE.TXT'\n",
    "tle_filename = 'tles/22335.txt'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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": 12,
   "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": 13,
   "id": "manufactured-sword",
   "metadata": {},
   "outputs": [],
   "source": [
    "if tle_filename == '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": 14,
   "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": 14,
     "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": 16,
   "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    183.904     4.40e-08          1.63e-08        1.79e-08\n",
      "0 VANGUARD 1                 4    181.921     5.58e-08          5.74e-09        4.88e-09\n",
      "0 VANGUARD 2                 3    228.803     1.54e-05          2.30e-05        2.18e-05\n",
      "0 VANGUARD 2                 3    231.099     1.53e-05          2.19e-05        2.07e-05\n",
      "0 VANGUARD R/B               4    107.637     7.22e-08          1.18e-08        2.15e-08\n",
      "0 VANGUARD R/B               4    105.711     7.06e-08          2.46e-08        2.11e-08\n",
      "0 VANGUARD R/B               4    103.844     6.95e-08          1.35e-08        3.40e-08\n",
      "0 VANGUARD R/B               4    195.807     4.64e-08          1.27e-08        3.14e-08\n",
      "0 VANGUARD R/B               4    194.680     4.62e-08          2.91e-08        1.30e-08\n",
      "0 VANGUARD 3                 3    215.412     1.82e-05          3.42e-05        3.59e-05\n",
      "0 VANGUARD 3                 3    214.231     1.80e-05          3.57e-05        3.73e-05\n",
      "0 VANGUARD 3                 3    214.231     1.80e-05          3.57e-05        3.73e-05\n",
      "0 EXPLORER 7                 4    140.142     7.55e-07          1.76e-06        1.80e-06\n",
      "0 EXPLORER 7                 4    139.420     7.71e-07          1.80e-06        1.83e-06\n",
      "0 EXPLORER 7                 4    138.520     7.89e-07          1.85e-06        1.88e-06\n",
      "0 TIROS 1                    4    152.113     3.40e-04          5.74e-04        5.49e-04\n",
      "0 TIROS 1                    4    151.325     3.94e-04          6.56e-04        6.27e-04\n",
      "0 TIROS 1                    4    151.270     3.98e-04          6.62e-04        6.32e-04\n",
      "0 TRANSIT 2A                 4    191.590     1.17e-07          2.43e-07        2.32e-07\n",
      "0 SOLRAD 1 (GREB)            3    107.042     2.94e-04          3.84e-04        4.00e-04\n",
      "0 SOLRAD 1 (GREB)            3    106.870     2.94e-04          3.83e-04        3.99e-04\n",
      "0 SOLRAD 1 (GREB)            3    106.801     2.94e-04          3.83e-04        3.99e-04\n",
      "0 THOR ABLESTAR R/B          4    197.650     1.40e-07          2.62e-07        2.58e-07\n",
      "0 THOR ABLESTAR R/B          4    197.701     1.36e-07          2.62e-07        2.48e-07\n",
      "0 DELTA 1 R/B                4    158.073     5.66e-07          1.51e-06        1.53e-06\n",
      "0 DELTA 1 R/B                4    157.132     5.41e-07          1.44e-06        1.45e-06\n",
      "0 ECHO 1 DEB (METAL OBJ)     4    187.870     2.91e-07          8.62e-07        8.51e-07\n",
      "0 ECHO 1 DEB (METAL OBJ)     4    187.681     2.65e-07          7.87e-07        7.83e-07\n",
      "0 ECHO 1 DEB (METAL OBJ)     4    156.790     3.06e-07          4.30e-07        4.42e-07\n",
      "0 COURIER 1B                 4    191.924     1.71e-07          5.05e-07        5.10e-07\n",
      "0 COURIER 1B                 4    193.137     1.92e-07          6.43e-07        6.19e-07\n",
      "0 THOR ABLESTAR R/B          4    168.825     2.53e-07          7.12e-07        6.84e-07\n",
      "0 THOR ABLESTAR R/B          4    169.432     2.45e-07          6.46e-07        6.25e-07\n",
      "0 SCOUT X-1 R/B              4    132.156     8.46e-08          3.04e-08        2.67e-08\n",
      "0 SCOUT X-1 R/B              4    127.187     8.37e-08          3.47e-08        3.04e-08\n",
      "0 SCOUT X-1 DEB              3    190.360     6.15e-05          1.02e-04        1.01e-04\n",
      "0 EXPLORER 11                4    182.006     5.30e-08          1.31e-08        2.06e-08\n",
      "0 EXPLORER 11                4    180.555     6.78e-08          1.47e-08        1.01e-08\n",
      "0 EXPLORER 11                4    178.709     7.70e-08          1.82e-08        1.52e-08\n",
      "0 THOR ABLE DEB (YO)         4    154.416     5.87e-05          1.48e-04        1.48e-04\n",
      "0 THOR ABLE DEB (YO)         4    154.416     5.87e-05          1.48e-04        1.48e-04\n",
      "0 TRANSIT 4A                 4    155.387     2.90e-06          7.46e-06        7.59e-06\n",
      "0 TRANSIT 4A                 4    155.509     2.86e-06          7.37e-06        7.49e-06\n",
      "0 SOLRAD 3/INJUN 1           4    147.157     3.43e-06          8.99e-06        9.10e-06\n",
      "0 THOR ABLESTAR R/B          4    152.443     2.51e-06          6.43e-06        6.55e-06\n",
      "0 THOR ABLESTAR R/B          4    152.653     2.48e-06          6.34e-06        6.44e-06\n",
      "0 THOR ABLESTAR DEB          4    140.395     1.51e-06          3.86e-06        3.91e-06\n",
      "0 THOR ABLESTAR DEB          4    140.623     1.51e-06          3.85e-06        3.90e-06\n",
      "0 THOR ABLESTAR DEB          3    143.127     4.49e-04          7.36e-04        7.71e-04\n",
      "0 THOR ABLESTAR DEB          3    143.127     4.49e-04          7.36e-04        7.71e-04\n",
      "0 THOR ABLESTAR DEB          4    136.709     2.31e-06          4.89e-06        5.05e-06\n",
      "0 THOR ABLESTAR DEB          4    136.624     2.31e-06          4.91e-06        5.07e-06\n",
      "0 THOR ABLESTAR DEB          4    132.499     1.76e-06          3.74e-06        3.88e-06\n",
      "0 THOR ABLESTAR DEB          4    132.547     1.65e-06          3.46e-06        3.59e-06\n",
      "0 THOR ABLESTAR DEB          4    132.547     1.65e-06          3.46e-06        3.59e-06\n",
      "0 THOR ABLESTAR DEB          4    128.886     7.14e-07          1.82e-06        1.81e-06\n",
      "0 THOR ABLESTAR DEB          4    177.190     2.19e-06          6.00e-06        6.09e-06\n",
      "0 THOR ABLESTAR DEB          4    177.176     2.22e-06          6.09e-06        6.18e-06\n",
      "0 THOR ABLESTAR DEB          4    175.648     3.10e-06          8.31e-06        8.41e-06\n",
      "0 THOR ABLESTAR DEB          4    175.625     3.10e-06          8.31e-06        8.40e-06\n",
      "0 THOR ABLESTAR DEB          4    175.537     3.10e-06          8.30e-06        8.41e-06\n",
      "0 THOR ABLESTAR DEB          4    141.721     4.74e-06          1.24e-05        1.26e-05\n",
      "0 THOR ABLESTAR DEB          4    141.827     4.80e-06          1.26e-05        1.27e-05\n",
      "0 THOR ABLESTAR DEB          4    171.907     3.70e-06          9.69e-06        9.79e-06\n",
      "0 THOR ABLESTAR DEB          4    171.830     3.68e-06          9.64e-06        9.74e-06\n",
      "0 THOR ABLESTAR DEB          4    171.724     3.66e-06          9.56e-06        9.67e-06\n",
      "0 THOR ABLESTAR DEB          4    141.519     1.84e-06          4.77e-06        4.83e-06\n",
      "0 THOR ABLESTAR DEB          4    141.553     1.85e-06          4.79e-06        4.85e-06\n",
      "0 THOR ABLESTAR DEB          4    141.702     1.86e-06          4.81e-06        4.87e-06\n",
      "0 THOR ABLESTAR DEB          4    146.402     2.66e-06          6.88e-06        6.99e-06\n",
      "0 THOR ABLESTAR DEB          4    146.523     2.66e-06          6.88e-06        6.97e-06\n",
      "0 THOR ABLESTAR DEB          4    146.672     2.66e-06          6.89e-06        6.97e-06\n",
      "0 THOR ABLESTAR DEB          4    164.751     4.43e-07          1.11e-06        1.08e-06\n",
      "0 THOR ABLESTAR DEB          4    164.390     4.03e-07          1.01e-06        9.88e-07\n",
      "0 THOR ABLESTAR DEB          4    164.390     4.03e-07          1.01e-06        9.88e-07\n",
      "0 THOR ABLESTAR DEB          4    174.438     1.32e-06          2.98e-06        3.09e-06\n",
      "0 THOR ABLESTAR DEB          4    174.519     1.28e-06          2.85e-06        2.96e-06\n",
      "0 THOR ABLESTAR DEB          4    174.519     1.28e-06          2.85e-06        2.96e-06\n",
      "0 THOR ABLESTAR DEB          4    157.777     4.65e-06          1.20e-05        1.22e-05\n",
      "0 THOR ABLESTAR DEB          4    157.854     4.56e-06          1.17e-05        1.20e-05\n",
      "0 THOR ABLESTAR DEB          4    157.854     4.56e-06          1.17e-05        1.20e-05\n",
      "0 THOR ABLESTAR DEB          4    166.767     2.93e-07          6.26e-07        6.33e-07\n",
      "0 THOR ABLESTAR DEB          4    166.823     3.09e-07          6.46e-07        6.75e-07\n",
      "0 THOR ABLESTAR DEB          4    144.330     7.93e-06          2.13e-05        2.14e-05\n",
      "0 THOR ABLESTAR DEB          4    144.483     8.08e-06          2.17e-05        2.19e-05\n",
      "0 THOR ABLESTAR DEB          4    144.483     8.08e-06          2.17e-05        2.19e-05\n",
      "0 THOR ABLESTAR DEB          4    161.274     1.23e-06          2.74e-06        2.80e-06\n",
      "0 THOR ABLESTAR DEB          4    161.239     1.38e-06          3.04e-06        3.11e-06\n",
      "0 THOR ABLESTAR DEB          4    161.220     1.46e-06          3.21e-06        3.28e-06\n",
      "0 THOR ABLESTAR DEB          4    143.587     9.45e-06          2.07e-05        2.13e-05\n",
      "0 THOR ABLESTAR DEB          4    143.580     9.42e-06          2.06e-05        2.12e-05\n",
      "0 THOR ABLESTAR DEB          4    158.466     1.92e-05          5.37e-05        5.34e-05\n",
      "0 THOR ABLESTAR DEB          4    140.393     5.73e-06          1.29e-05        1.33e-05\n",
      "0 THOR ABLESTAR DEB          4    140.377     5.60e-06          1.26e-05        1.30e-05\n",
      "0 THOR ABLESTAR DEB          4    140.377     5.60e-06          1.26e-05        1.30e-05\n",
      "0 THOR ABLESTAR DEB          4    165.458     3.77e-06          9.57e-06        9.69e-06\n",
      "0 THOR ABLESTAR DEB          4    165.323     3.65e-06          9.25e-06        9.36e-06\n",
      "0 THOR ABLESTAR DEB          4    165.270     3.60e-06          9.13e-06        9.22e-06\n",
      "0 THOR ABLESTAR DEB          4    171.772     1.90e-07          4.52e-07        4.40e-07\n",
      "0 THOR ABLESTAR DEB          4    173.831     4.98e-06          1.34e-05        1.36e-05\n",
      "0 THOR ABLESTAR DEB          4    173.812     4.97e-06          1.33e-05        1.35e-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": "color-accommodation",
   "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
}