SkyCalibration-LeastSquare.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "described-carpet",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "import statistics as stat\n",
    "import math\n",
    "import scipy as sci\n",
    "from matplotlib import pyplot as plt\n",
    "import PrizmCalibration as cl\n",
    "import prizmatoid as pzt\n",
    "try:\n",
    "    reload(pzt)\n",
    "    reload(cl)\n",
    "except:\n",
    "    from importlib import reload\n",
    "    reload(pzt)\n",
    "    reload(cl)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "absolute-guest",
   "metadata": {},
   "outputs": [],
   "source": [
    "flow = 60\n",
    "fhigh = 120\n",
    "minperbin = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "meaning-problem",
   "metadata": {},
   "source": [
    "# Prepping Data (don't need to do if using saved arrays)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "younger-exclusive",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/kellyforan/Documents/Prizm/Calibration/PrizmCalibration.py:218: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  binned[i,:] = np.sum(tmp, axis=0)/len(tmp)\n"
     ]
    }
   ],
   "source": [
    "cl.GSM_calibration_data(minperbin, flow, fhigh, 100, 'EW')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "official-fourth",
   "metadata": {},
   "source": [
    "# Opening saved data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "patent-memorabilia",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('100MHz_EW_antenna_average.npy', 'rb') as f:\n",
    "    alldataavg = np.load(f)\n",
    "with open('100MHz_EW_GSM_average.npy', 'rb') as g:\n",
    "    allgsmavg = np.load(g)\n",
    "with open('100MHz_EW_antenna_allarrays.npy', 'rb') as h:\n",
    "    alldata = np.load(h)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "revolutionary-corrections",
   "metadata": {},
   "source": [
    "# LEAST SQUARES\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "dimensional-airport",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Tcd simulation values\n",
    "flow = 60\n",
    "fhigh = 120\n",
    "\n",
    "Tcd = np.array([-1.01879219e-02, -1.38655564e-02, -1.85246433e-02, -2.42953560e-02, -3.12793370e-02, -3.95323534e-02, -4.90465721e-02, -5.97345379e-02,\n",
    " -7.14172654e-02, -8.38188667e-02, -9.65697732e-02, -1.09219830e-01, -1.21261418e-01, -1.32161405e-01, -1.41399382e-01, -1.48508474e-01,\n",
    " -1.53114322e-01, -1.54967725e-01, -1.53966968e-01, -1.50167074e-01, -1.43774789e-01, -1.35129950e-01, -1.24675568e-01, -1.12920281e-01,\n",
    " -1.00397618e-01, -8.76265518e-02, -7.50773434e-02, -6.31455762e-02, -5.21360139e-02, -4.22565052e-02, -3.36209660e-02, -2.62595729e-02,\n",
    " -2.01338111e-02, -1.51539266e-02, -1.11965758e-02, -8.12093687e-03, -5.78213066e-03, -4.04138781e-03, -2.77289876e-03, -1.86766232e-03,\n",
    " -1.23487660e-03, -8.01511254e-04, -5.10689051e-04, -3.19421635e-04, -1.96125006e-04, -1.18212246e-04, -6.99443807e-05, -4.06259985e-05,\n",
    " -2.31641370e-05, -1.29654938e-05, -7.4827469e-06])\n",
    "\n",
    "start = int((flow-50)/2)\n",
    "end = 51-int((150-fhigh)/2)\n",
    "Tcd = Tcd[start:end]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "educated-virtue",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create data to fit model to (d = Tgsm + Tcd)\n",
    "d = allgsmavg + np.repeat(Tcd[np.newaxis,:], allgsmavg.shape[0], axis = 0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "satisfactory-hayes",
   "metadata": {},
   "source": [
    "# Dealing with N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "worthy-politics",
   "metadata": {},
   "outputs": [],
   "source": [
    "#find variance in power at one frequency over the 25 days of data\n",
    "tmp = []\n",
    "alldatavar = np.zeros((alldata.shape[1], alldata.shape[2]))\n",
    "\n",
    "for j in range(0, alldata.shape[1]):\n",
    "    for k in range(0, alldata.shape[2]):\n",
    "        for i in range(0, alldata.shape[0]):\n",
    "            if alldata[i,j,k] != 0:\n",
    "                tmp.append(alldata[i,j,k])\n",
    "        alldatavar[j, k] = np.nanvar(tmp)\n",
    "        tmp = []\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "nasty-boundary",
   "metadata": {},
   "outputs": [],
   "source": [
    "#find variance in power\n",
    "timevar = np.nanvar(alldataavg, axis = 0)\n",
    "freqvar = np.nanvar(alldataavg, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "strong-cliff",
   "metadata": {},
   "source": [
    "# Fit versus t at each f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "little-hurricane",
   "metadata": {},
   "outputs": [],
   "source": [
    "polydim = 8\n",
    "a = np.zeros((polydim, alldataavg.shape[1]))\n",
    "\n",
    "for j in range(0, alldataavg.shape[1]):\n",
    "        Ninv = np.identity(144)*(1/alldatavar[:,j])\n",
    "        P = alldataavg[:,j]\n",
    "        columns = np.zeros((len(P), polydim))\n",
    "        columns[:,0] = P\n",
    "        for k in range(1, polydim):\n",
    "            temparray = np.zeros(len(P))\n",
    "            for i in range(0, alldataavg.shape[0]):\n",
    "                temparray[i]= P[i]*(i*120)**k\n",
    "            columns[:,k] = temparray  \n",
    "        data = np.stack(columns, axis=0)\n",
    "        data[np.isnan(data)] = 0\n",
    "        \n",
    "        d[np.isnan(d)] = 0\n",
    "        AT = np.transpose(data)\n",
    "        first = np.matmul(np.matmul(AT, Ninv), data)\n",
    "        firstinv = np.linalg.inv(first)\n",
    "        second = np.matmul(np.matmul(AT, Ninv), d[:,j])\n",
    "        final = np.matmul(firstinv, second)  \n",
    "        a[:,j] = final       "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "parliamentary-estate",
   "metadata": {},
   "outputs": [],
   "source": [
    "ktmp = np.zeros((alldataavg.shape[0], alldataavg.shape[1], polydim))\n",
    "for i in range(0,alldataavg.shape[0]):\n",
    "    for j in range(0,alldataavg.shape[1]):\n",
    "        for l in range(0, polydim):\n",
    "            ktmp[i,j,l] = a[l,j]*((i*120)**l)\n",
    "k = np.nansum(ktmp, axis=2)\n",
    "kavg = np.nanmean(k, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "stable-enterprise",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([<matplotlib.axis.XTick at 0x128bfe410>,\n",
       "  <matplotlib.axis.XTick at 0x128c0d1d0>,\n",
       "  <matplotlib.axis.XTick at 0x128c0d2d0>,\n",
       "  <matplotlib.axis.XTick at 0x129111650>],\n",
       " [Text(0, 0, '60'), Text(10, 0, '80'), Text(20, 0, '100'), Text(30, 0, '120')])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAERCAYAAABrWly6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAlu0lEQVR4nO3deZhU9Z3v8fe3970buhtooOkFEaIgoIggajSrGCdOZpKreVzikmFMnEySmWRmcufJ7JPZcjPR5EZi3GI2M/GajJOoiUvijgqyiLIIzdLN2jT0vnd/7x91wLbplgKq+9TyeT1PPXXqnFOnvnDw469+9Tu/Y+6OiIgkvrSwCxARkdhQoIuIJAkFuohIklCgi4gkCQW6iEiSUKCLiCSJUAPdzO41s4NmtjFGx5thZr8xs01m9qaZVcfiuCIiiSDsFvr9wOUxPN4DwH+4+3uAxcDBGB5bRCSuhRro7v4scHjoOjObaWaPm9kaM3vOzOZEcywzOwvIcPcngmO3u3tn7KsWEYlPYbfQR3IX8Dl3Pw/4EvCdKN93JtBsZg+b2Voz+w8zSx+zKkVE4kxG2AUMZWYFwIXAz8zs6OrsYNsfAP8wwtv2uPuHifxZLgYWAruBnwI3AveMbdUiIvEhrgKdyDeGZndfMHyDuz8MPPwu720A1rp7HYCZ/QJYggJdRFJEXHW5uHsrsMPMPgFgEfOjfPurwAQzKw9evw94cwzKFBGJS2EPW/wJ8BIw28wazOwW4FrgFjNbD7wBXBXNsdx9gEif+1Nm9jpgwPfGpnIRkfhjmj5XRCQ5xFWXi4iInLrQfhQtKyvz6urqsD5eRCQhrVmz5pC7l4+0LbRAr66uZvXq1WF9vIhIQjKzXaNtU5eLiEiSUKCLiCQJBbqISJJQoIuIJAkFuohIklCgi4gkCQW6iEiSSLhA37y/lX99bDOt3X1hlyIiElcSLtB3N3Wy8pntbD/YHnYpIiJxJeECvba8AIC6xo6QKxERiS8JF+gzJuaRnmbsOKRAFxEZKuECPSsjjRkT86g7pC4XEZGhEi7QAWrK8tXlIiIyTEIGem1ZPjsOdTA4qJtziIgclZiBXl5AT/8ge1u6wi5FRCRuJGig5wMa6SIiMlRiBnrZ0UDXD6MiIkclZKCXF2ZTkJ2hoYsiIkMkZKCbGbXl+dQp0EVEjknIQIdIt4v60EVE3nbCQDez2Wa2bsij1cy+MGwfM7M7zGybmW0ws3PHrOJATVkBe5q76OodGOuPEhFJCCcMdHff4u4L3H0BcB7QCfx82G7LgVnBYwVwZ4zrPM7RkS47m9RKFxGBk+9yeT+w3d13DVt/FfCAR6wCSsysIiYVjkJDF0VE3ulkA/0a4CcjrJ8G1A953RCsewczW2Fmq81sdWNj40l+9DvVaOiiiMg7RB3oZpYFfBT42UibR1h33HX57n6Xuy9y90Xl5eXRVzmCvKwMKopzNNJFRCRwMi305cBr7n5ghG0NQOWQ19OBvadTWDQ0dFFE5G0nE+ifZOTuFoBHgBuC0S5LgBZ333fa1Z1AbVkBdY3tuGuSLhGRqALdzPKADwIPD1l3q5ndGrx8FKgDtgHfAz4b4zpHVFOWT1t3P4fae8fj40RE4lpGNDu5eydQOmzdyiHLDtwW29JO7OhIlx2HOigvzB7vjxcRiSsJe6UowMxj9xfVSBcRkYQO9KkluWRlpOmHUREREjzQ09OM6tI8tdBFREjwQIdgpIta6CIiSRDo5fnsbuqkb2Aw7FJEREKVBIFeQP+gU3+4M+xSRERClfCB/vacLup2EZHUlvCBPnPIWHQRkVSW8IFekpfFxPws6g5ppIuIpLaED3SIdLtsV5eLiKS4pAj02rJ8dbmISMpLjkAvL6CxrYe27r6wSxERCU2SBLpGuoiIJEegHx26qB9GRSSFJUWgzyjNI81gh1roIpLCkiLQszPSqZyYx3b9MCoiKSwpAh0i3S7qQxeRVJY0gV5TVsCOQ+0MDur+oiKSmpIm0GvL8+nuG2R/a3fYpYiIhCKpAh00dFFEUlfSBPqx+4tq6KKIpKikCfRJhdnkZ6WrhS4iKSuqQDezEjN7yMw2m9kmM1s6bPulZtZiZuuCx9+MTbnvWiM15fm6HZ2IpKyMKPe7HXjc3T9uZllA3gj7POfuV8autJNXW1bAa7uPhFmCiEhoTthCN7Mi4BLgHgB373X35jGu65TUlOWzp7mL7r6BsEsRERl30XS51AKNwH1mttbM7jaz/BH2W2pm683sMTM7e6QDmdkKM1ttZqsbGxtPp+6RCy3Pxx12Nen+oiKSeqIJ9AzgXOBOd18IdAB/NWyf14Aqd58PfAv4xUgHcve73H2Ruy8qLy8/9apHcWykS6NGuohI6okm0BuABnd/OXj9EJGAP8bdW929PVh+FMg0s7KYVhqFYzeM1g+jIpKCThjo7r4fqDez2cGq9wNvDt3HzKaYmQXLi4PjNsW41hPKz85gclE229VCF5EUFO0ol88BPwpGuNQBN5nZrQDuvhL4OPAZM+sHuoBr3D2USVVqywp0OzoRSUlRBbq7rwMWDVu9csj2bwPfjl1Zp662PJ9fbtiHuxN8aRARSQlJc6XoUbXlBbR09XG4ozfsUkRExlXyBbp+GBWRFJV8gR7Muqjb0YlIqkm6QJ8+IY/MdGO7Zl0UkRSTdIGenmZUlep2dCKSepIu0CHSj66hiyKSapIz0MsL2NXUQf/AYNiliIiMmyQN9Hz6BpyGI11hlyIiMm6SM9CPDV3UD6MikjqSM9CPzbqofnQRSR1JGegT87MoycvUxUUiklKSMtAh0u2iedFFJJUkbaDXlBWoy0VEUkrSBnpteT4H23po7+kPuxQRkXGRtIE+U3O6iEiKSdpAP2NSZKTL5v2tIVciIjI+kjbQa8oKKMjOYH1Dc9iliIiMi6QN9PQ045zpxayrbw67FBGRcZG0gQ6wcEYJm/e10d03EHYpIiJjLqkDfUHlBPoHnY17WsIuRURkzCV5oJcAqNtFRFJCVIFuZiVm9pCZbTazTWa2dNh2M7M7zGybmW0ws3PHptyTU16YzbSSXNYq0EUkBWREud/twOPu/nEzywLyhm1fDswKHhcAdwbPoVswo4R1u5vDLkNEZMydsIVuZkXAJcA9AO7e6+7Nw3a7CnjAI1YBJWZWEetiT8XCyhL2NHfR2NYTdikiImMqmi6XWqARuM/M1prZ3WaWP2yfaUD9kNcNwbp3MLMVZrbazFY3NjaectEnQ/3oIpIqogn0DOBc4E53Xwh0AH81bB8b4X1+3Ar3u9x9kbsvKi8vP+liT8XcacVkpBnr6o+My+eJiIQlmkBvABrc/eXg9UNEAn74PpVDXk8H9p5+eacvJzOdORWFaqGLSNI7YaC7+36g3sxmB6veD7w5bLdHgBuC0S5LgBZ33xfbUk/dgsoSNtS3MDh43JcGEZGkEe049M8BPzKzDcAC4GtmdquZ3RpsfxSoA7YB3wM+G+tCT8eCygm09fSzXTe8EJEkFtWwRXdfBywatnrlkO0O3Ba7smLr6A+ja+ubmTW5MNxiRETGSFJfKXpUbVk+hTkZ6kcXkaSWEoGelmYsqNQFRiKS3FIi0CHS7bLlQBtdvZp5UUSSU0oF+sCg87pmXhSRJJVSgQ6wdrcuMBKR5JQygV5akE3lxFz9MCoiSStlAh0i49EV6CKSrFIs0EvY19LNgdbusEsREYm5lAt0gLUavigiSSilAv3sqUVkppu6XUQkKaVUoOdkpvOeiiJNpSsiSSmlAh0i3S6vN7QwoJkXRSTJpGSgd/QO8NbBtrBLERGJqZQMdEDzuohI0km5QK8py6c4N1M/jIpI0km5QDcz5leWKNBFJOmkXKBDpNtl64E2Onr6wy5FRCRmUjLQF1aWMOiwoUEzL4pI8kjJQJ9/9IdRdbuISBJJyUCfmJ9FVWmeLjASkaSSkoEOkX50tdBFJJmkdKAfaO1hX0tX2KWIiMREVIFuZjvN7HUzW2dmq0fYfqmZtQTb15nZ38S+1NjSBUYikmwyTmLfy9z90Ltsf87drzzdgsbLWVOLyEpPY119M8vnVYRdjojIaUvZLpfsjHTeM7WItepHF5EkEW2gO/AbM1tjZitG2Wepma03s8fM7OyRdjCzFWa22sxWNzY2nlLBsbQwmHmxf2Aw7FJERE5btIG+zN3PBZYDt5nZJcO2vwZUuft84FvAL0Y6iLvf5e6L3H1ReXn5qdYcMwsqS+jqG2DrgfawSxEROW1RBbq77w2eDwI/BxYP297q7u3B8qNAppmVxbjWmFugC4xEJImcMNDNLN/MCo8uAx8CNg7bZ4qZWbC8ODhuU+zLja2q0jwm5GXqAiMRSQrRjHKZDPw8yOsM4Mfu/riZ3Qrg7iuBjwOfMbN+oAu4xt3j/pZAmnlRRJLJCQPd3euA+SOsXzlk+dvAt2Nb2vhYWDmBZ7Y20tbdR2FOZtjliIicspQdtnjUghkluMPrmnlRRBKcAn16CWbw8o7DYZciInJaUj7Qi/MyWVQ1gV+/sT/sUkRETkvKBzrA8rkVbN7fRl2jxqOLSOJSoAOXz50CwGMb1UoXkcSlQAemluSycEYJj23cF3YpIiKnTIEeuGJuBRv3tLK7qTPsUkRETokCPfB2t4ta6SKSmBTogcqJeZwzvZhH1Y8uIglKgT7E8rkVrK9vpuGIul1EJPEo0IdYHnS7PK5WuogkIAX6ENVl+ZxVUaThiyKSkBTow1wxbwprdh1hf0t32KWIiJwUBfowl8+N3DD6cY12EZEEo0Af5oxJBZw5uUCjXUQk4SjQR7B8bgWv7jzMwTZ1u4hI4lCgj+CKeRW4w6/fOBB2KSIiUVOgj+DMyQXUluerH11EEooCfQRmxhVzK1hVd5im9p6wyxERiYoCfRTL501hYNB54k11u4hIYlCgj+KsiiKqSvM02kVEEkZUgW5mO83sdTNbZ2arR9huZnaHmW0zsw1mdm7sSx1fZsbyuRW8uO0QzZ29YZcjInJCJ9NCv8zdF7j7ohG2LQdmBY8VwJ2xKC5sV8ybQr+6XUQkQcSqy+Uq4AGPWAWUmFlFjI4dmnnTiplWkqu5XUQkIUQb6A78xszWmNmKEbZPA+qHvG4I1r2Dma0ws9VmtrqxsfHkqx1nkW6XKTz3ViOt3X1hlyMi8q6iDfRl7n4uka6V28zskmHbbYT3+HEr3O9y90Xuvqi8vPwkSw3H8nkV9A04T21St4uIxLeoAt3d9wbPB4GfA4uH7dIAVA55PR3YG4sCw7awsoQpRTk8+rq6XUQkvp0w0M0s38wKjy4DHwI2DtvtEeCGYLTLEqDF3ZPiMsu0NOPyuVN4Zmsj7T39YZcjIjKqaFrok4HnzWw98ArwK3d/3MxuNbNbg30eBeqAbcD3gM+OSbUhuWJeBb39gzy9+WDYpYiIjCrjRDu4ex0wf4T1K4csO3BbbEuLH+dVTaC8MJvHXt/HR+dPDbscEZER6UrRKKSnGZefPYXfbjlIZ6+6XUQkPinQo7R83hS6+wZ5Zkv8D7cUkdSkQI/S4uqJlBVk8eCr9SfeWUQkBAr0KGWkp3HTshqe2drIuvrmsMsRETmOAv0kfOrCakryMrn9ya1hlyIichwF+kkoyM7gjy6u5bdbGlmvVrqIxBkF+km6YWlVpJX+1FthlyIi8g4K9JNUmJPJpy+q4enNB9nQ0Bx2OSIixyjQT8GnLqymODeT259UK11E4ocC/RQcbaU/tfkgrze0hF2OiAigQD9ln1oWtNKf0ogXEYkPCvRTVJSTyS0X1fDkpoNs3KNWuoiET4F+Gm5cVk1RTgbfVF+6SEobHHTae/qJzFMYnhPOtiijK8rJ5NMX1/KNJ7aycU8Lc6cVh12SiMRQW3cfT28+yKH2Xlq6+mjp7KW5q4/mzj6ag9ctXX20dPUx6JCbmc6MiXlUTsyjqjSPGRPzmBE8T5+QS3ZG+pjWq0A/TTcuq+bu5+q4/am3+N4Ni8IuR0RiYHdTJ/e/uJP/Wl1/7MY2ZpFGXEleJsW5kceMiXmU5EbW5Wdn0NjWw66mTuoPd/LCtkN09Q0cO6YZVBTlMKM0jz88dzqfWFQ52sefMgX6aYr0pdfyn0+qlS6SyNydV3Yc5t4XdvDEmwdIM+PKcyq4fmkVM8sLKMzJJD1tpNsnj368xvYe6g93squpk92Hg0dTJ529Ayc+wCmwsPp8Fi1a5KtXrw7ls2OtpauPi/7taZbWlnKXWukiCaW3f5BfbtjLvS/sYOOeVkryMrn2ghlcv6SaKcU5YZd3HDNb4+4jBo1a6DFQnBsZ8fLNJ9/ijb0tnD1VrXSReNfU3sOPX97NA6t20djWwxmTCvjax+bxsYXTyM0a277usaJAj5GbltVwz/M7uOOpt/ju9Wqli8Srzftbue/5nfxi3R56+gd575nl3PKJGi6eVYZZ9F0q8UiBHiPFuZncvKyG2596izf3tnLW1KKwSxKRwOCg8/Tmg9z7wg5e3N5ETmYaf3jedG5eVs0ZkwrDLi9mFOgxdPOyGu4NWukrrz8v7HJEUl5bdx8PrWng/hd3squpk4riHP7y8jl8cnElJXlZYZcXcwr0GCrOy+Smi2q446m32LSvlfdUqJUuEoZdTR3c/+JOfra6gfaefs6dUcKXPzybD589hcz05L2eMupAN7N0YDWwx92vHLbtUuC/gR3Bqofd/R9iVGNCuWVZDfc9v4N/f3wz9954fsL3yYkkAndnT3MXa3Yd4Zcb9vHkpgOkm/GRcyq4aVkNCypLwi5xXJxMC/3zwCZgtGbnc8ODPhUV52Xy+Q/M4p9+tYkfv7Kbay+oCrskkaTTNzDIpn2trN55hDW7j7Bm5xH2t3YDMCEvk9suPYPrl1YxuSj+hh2OpagC3cymAx8B/hn4szGtKAncHNxM+h9/+SaLqycya3Ly/OgiMl7cne6+Qdp7+uns7afuUAdrdh5h9a7DrK9vOXYV5tTiHM6vmciiqgmcVzWBOVMKyUjibpV3E20L/ZvAXwDvlkxLzWw9sBf4kru/MXwHM1sBrACYMWPGyVWaQNLSjP/ziflcfvtz/OmD6/j5Zy8kJzMxx7WKxEpv/yD7W7ppONJJQ3MXDUe62HOki+bOXjp6++noGQie++kMlgeHXfeYnmacVVHE1edXcl4Q4FNLcsP5A8WhEwa6mV0JHHT3NUFf+UheA6rcvd3MrgB+AcwavpO73wXcBZErRU+x5oQwqSiHr3/iHG6+fzX//vgW/ub3zgq7JJFxcbCtm+ffOsT2xnb2HAmCu7mL/a3dDL0w3QwmF+ZQWpBFflYGZQVZVGXnkZ+VQX52BvnZ6ZHnrMjzlOIc5k8vIT9bYzlGE83fzDLgo0FQ5wBFZvZDd7/u6A7u3jpk+VEz+46Zlbn7odiXnDjeN2cyN15Yzb0v7ODiM8u4bPaksEsSibmBQWddfTO/23KQ3245yMY9kThITzMqinOYPiGXC2eWMW1CLtMn5DK9JJfpE/KYUpxDVkZqdo2MlZOayyVooX9phFEuU4AD7u5mthh4iEiLfdSDJ9NcLu+mu2+Aq779Ak0dPTz2+UsoL8wOuySR03ako5dn32rk6c0HeXZrI0c6+0gzOK9qApfOnsSls8uZPTl1+7LH0pjM5WJmtwK4+0rg48BnzKwf6AKuebcwTyU5menc8cmFfPTbz/Oln63nvhvPJ+0kZmwTiRcHWrt5aE0DT206wNr6ZtyhND+Ly+ZM4rLZk7h4VllSXqyTSDTb4jj5wUs7+ep/v8FXrzyLWy6qCbsckai4Oy/vOMwPXtrFr9/YT/+gM396MZfOnsT75kxi3rRiNVDGmWZbjAPXLanima2H+LfHNrOkdqJmZJS41tHTz8Nr9/DDl3ax5UBbZK6ii2q49oIZVJXmh12ejEIt9HF0uKOXy7/5LIU5GfzP5y4iL0v/P5X4su1gOz9ctYuH1kQumZ87rYgbllTze/OnJuyUsslGLfQ4MTE/i/+8egHX3fMy//jLTfzLH8wLuyQROnv7eXZrIz9YtYsXtjWRlZ7GFfOmcMOF1SysLNH0FQlEgT7Olp1RxopLavnuM3W898wyLp9bEXZJkmK6egd4bfcRXtrexKq6JtY3NNM34EwtzuHLH57N1edXUlag0ViJSIEegj//4Gxe3NbEX/6/15lfWUJFsa50k7HT3TfAa7uOsKquiZfqmlhXHwnw9DRj3rRibrmolqUzS1k2s1TDDBOc+tBDUtfYzpXfep45Uwq576bFFOdmhl2SJIHuvgG2HWxn075WtuxvY8OeFtbtbqZ3YJA0g3nTilkys5QltaUsqppAYY7+3SWad+tDV6CH6Fcb9vH5B9dSW57PvTeez/QJeWGXJAnC3Wk40sXm/W1s2d/Kpv1tbN7Xys6mTgaCCVCyM9KYM6WQC2pLWVI7kfOrJyrAk4B+FI1THzmnggl5mfzxD9fwse+8yD2fWsQ500vCLkviREdPP3uau2g40kn94chzQzA3ys5DHbT19B/bd8bEPGZPKeQj8yqYPaWIORWFVJfmk64x4ilFLfQ48NaBNm6871UOd/RyxycX8sGzJoddksRI38Agv918kPojXfQPDNI/6AwM+rHl/kGnf8AZGBykb9Bp6eyLBPiRLg539L7jWDmZaUyfkMf0CblUTcw7FtxnTi6kQBNWpQx1uSSAg23d/NH3V7NhTwt/c+VZ3LRMV5MmsoNt3Tz4Sj0/fnn3sRsvDJeRZmSkGxlpaaSnGZnpRlFOJtMnRkI78sijMnguK8jSEEJRl0simFSYw4MrlvL5B9fy9//zJruaOvnqlWfpK3MCcXde232EB17axaOv76NvwLl4Vhn/+PtzWVw9kYx0Iz3NyEiLPCucJdYU6HEkNyudO687j689uol7nt9Bw5Eu7vjkAl1RGue6egd4ZP0eHnhpF2/sbaUwO4PrllRx/ZIqassLwi5PUoiSIs6kpxlfvfIsqkrz+LtH3uDq767inhsXMakwte6NmAjqD3fyg1W7+Omr9bR09TF7ciH/9Ptz+djCaboJg4RC/+ri1A1Lq5lWksuf/HgtH/u/L3Lvjecze4ruTRoP9jZ3ccdTb/GzNQ0AXH72FK5fWsUFNRPVjSKh0o+icW7jnhZuvv9V2rr7uXFZNSsurmVCvuacDsOh9h6+89vt/HDVLgCuXTKDP75kJlOK9e1Jxo9GuSS4fS1dfO3Rzfxyw17yszK4+aIabrmoRleXjpOWrj7ufq6Oe57fQXffAJ84r5I//cAspunmxBICBXqS2LK/jW8+uZXHNu6nKCeDP7q4lpsuqtEY5DHS2dvP91/cxcpnttPS1ceV51TwxQ+eyUz90CkhUqAnmY17Wvjmk2/x5KYDkStN3zuTG5ZWaTRMjPT2D/Lgq7v51tPbaGzr4X1zJvHnHzpTNyWRuKBAT1Lr65v5xhNbeWZrI2UFWdz63plct6SKnEzdiOBk7W/p5uUdkelkf7elkX0t3SyumchffHg2i6onhl2eyDEK9CS3ZtdhvvHEVl7Y1kRZQRYX1JQyv7KY+dNLmDutWEPoRrC/pZtVdU1BiB9mx6EOAApzMrigZiLXL63mklllGrUicUeBniJW1TXxo5d3s67+CPWHuwBIMzhzciHzp5cwv7KE+ZXFzJ5cmFLzXnf1DrDjUAdbDrTyct1hVtU1sbOpE3g7wJfURqaUfU9Fka7OlbgWk0A3s3RgNbDH3a8cts2A24ErgE7gRnd/7d2Op0AfW03tPWxoaGFdfTPrG5pZX9/Mkc4+IDLJ01kVRUwtyaU0P4vSgmxKC7LeXs7PojQ/m6LcjIRpobo7+1u7qWvsoK6xne2NHWxvbKeusYM9zV3H9ivKyWBxTWQ6WQW4JKJYzeXyeWATUDTCtuXArOBxAXBn8CwhKS3I5rI5k7hsziQgEni7D3dGAr6+hY17W3hzbyuH2nto7e4f8RiZ6caEvCyyM9NINyPNjLQ0I80gzSLzkQxdl52RRm5mOnlZGeRmpZOXlU5uZvqw5Qzys9LJz84IHunkZ2VQELzOynjnNwd3p72nn0PtvRxq76GxrYdD7T0cauuhsb2HxrZe9rd2saOxg47egWPvy8tKZ2Z5AYuqJ3B1eSW15fmcMamAWZMKFeCStKIKdDObDnwE+Gfgz0bY5SrgAY8091eZWYmZVbj7vtiVKqfDzKgqzaeqNJ+rFkx7x7be/kGOdEYCs6m9l8MdwXJHL4fbe+kdGGTQI9O+DrozOAgD7vixdTDoTk/fII3tPXT1dtLVO0Bn3wCdvQP09g9GXWdWelok5IN+/0PtPXT3Hf9+MyjNz6KsIJtJRTksqprIzPJ8ZpYXUFtewOSi7IT5diESK9G20L8J/AUw2rXn04D6Ia8bgnXvCHQzWwGsAJgxY8bJ1CljKCsjjclFOUwuGpsrHgcGna6+ATp7++nqHaCjJ7Lc3tNPR88AHT2R5ci6yOuOnn4cKCuIhHZZQTblhZHnssIsJuZlpdTvACLROGGgm9mVwEF3X2Nml4622wjrjuucd/e7gLsg0ocefZmSyNLTjILsDF0AJTLGomniLAM+amY7gQeB95nZD4ft0wBUDnk9HdgbkwpFRCQqJwx0d/+Ku09392rgGuBpd79u2G6PADdYxBKgRf3nIiLj65S/A5vZrQDuvhJ4lMiQxW1Ehi3eFJPqREQkaicV6O7+O+B3wfLKIesduC2WhYmIyMnRMAERkSShQBcRSRIKdBGRJKFAFxFJEqHNtmhmjcCuU3x7GXAohuVIbOi8xB+dk/h0Ouelyt3LR9oQWqCfDjNbPdpsYxIenZf4o3MSn8bqvKjLRUQkSSjQRUSSRKIG+l1hFyAj0nmJPzon8WlMzktC9qGLiMjxErWFLiIiwyjQRUSSREIEenBLu4fMbLOZbTKzpWY20cyeMLO3gucJYdeZSszsi2b2hpltNLOfmFmOzsn4M7N7zeygmW0csm7U82BmXzGzbWa2xcw+HE7VyW2Uc/IfQX5tMLOfm1nJkG0xOycJEejA7cDj7j4HmE/kZtV/BTzl7rOAp4LXMg7MbBrwp8Aid58LpBOZK1/nZPzdD1w+bN2I58HMziJyns4O3vMdM0sfv1JTxv0cf06eAOa6+znAVuArEPtzEveBbmZFwCXAPQDu3uvuzURuTP39YLfvA78fRn0pLAPINbMMII/IHap0TsaZuz8LHB62erTzcBXwoLv3uPsOIvcvWDwedaaSkc6Ju//G3fuDl6uI3NUNYnxO4j7QgVqgEbjPzNaa2d1mlg9MPnpXpOB5UphFphJ33wN8HdhN5EbgLe7+G3RO4sVo52G0m7nL+LoZeCxYjuk5SYRAzwDOBe5094VAB/oqH6qgT/YqoAaYCuSb2fDbEkr8iepm7jJ2zOyvgX7gR0dXjbDbKZ+TRAj0BqDB3V8OXj9EJOAPmFkFQPB8MKT6UtEHgB3u3ujufcDDwIXonMSL0c6DbuYeIjP7FHAlcK2/fQFQTM9J3Ae6u+8H6s1sdrDq/cCbRG5M/alg3aeA/w6hvFS1G1hiZnlmZkTOySZ0TuLFaOfhEeAaM8s2sxpgFvBKCPWlHDO7HPhL4KPu3jlkU0zPSUJcKWpmC4C7gSygjshNqNOA/wJmEAmYT7j78B+HZIyY2d8DVxP5+rgW+DRQgM7JuDKznwCXEpmO9QDwt8AvGOU8BF/5byZy3r7g7o8df1Q5HaOck68A2UBTsNsqd7812D9m5yQhAl1ERE4s7rtcREQkOgp0EZEkoUAXEUkSCnQRkSShQBcRSRIKdBl3ZjZgZuuGPKrDrilWzGyhmd0dLN9oZm5m7x+y/WPBuo8Hr39nZouGbK8eOkvfCMcvN7PHx/LPIIlLgS5h6HL3BUMeO49usIhE/nf5v4FvDXn9OvDJIa+vAdaf6sHdvRHYZ2bLTvUYkrwS+T8cSRJBq3STmX0HeA2oNLMvm9mrwfzRfz9k378O5o1+MpiH/UvB+mMtXTMrM7OdwXJ6MBf10WP9cbD+0uA9R+fZ/1Fw1Stmdr6ZvWhm683sFTMrNLPnggvcjtbxgpmdM+zPUQic4+5DA/s5YLGZZZpZAXAGsC7Kv5e7h3yLaTSzvw02/QK4Nsq/XkkhGWEXICkp18zWBcs7gC8Cs4Gb3P2zZvYhIpdALyYyedEjZnYJkYnZrgEWEvm3+xqw5gSfdQuR2SDPN7Ns4AUz+02wbSGReaj3Ai8Ay8zsFeCnwNXu/mowfXMXkSuVbwS+YGZnAtnuvmHYZy0ChneXOPAk8GGgmMil3jXD9vmRmXUFy1nAIIC7fxrAzKqAXxOZZxtgNfBPJ/hzSwpSoEsYutx9wdEXQR/6LndfFaz6UPBYG7wuIBLwhcDPj86FYWaPRPFZHwLOOdpnTSRUZwG9wCvu3hAcax1QDbQA+9z9VQB3bw22/wz4qpl9mchl2veP8FkVRKZ6Hu5BIjcEKQb+nEi3zFDXuvvq4HOqgV8e3WBmOcDPgD9x913B6oNEZrkUeQcFusSLjiHLBvyLu3936A5m9gVGn1q0n7e7EHOGHetz7v7rYce6FOgZsmqAyH8PNtJnuHunmT1BZNrg/0WkNT5c17DPPvreV8xsLpH/kW0NenaitRJ42N2fHLIuJ/gskXdQH7rEo18DNwd9zpjZNDObBDwLfMzMcoP+6t8b8p6dwHnB8seHHeszZpYZHOtMi9wgZTSbgalmdn6wf6FF7soEkW6XO4BXR5l0bBORPvKRfIXjW+bvysxuAwrd/V+HbTqT47t2RNRCl/jj7r8xs/cALwWt2XbgOnd/zcx+SuRHxV1EfnA86uvAf5nZ9cDTQ9bfTaQr5bXgR89G3uXWeO7ea2ZXA98ys1wiLeEPAO3uvsbMWoH7RnnvZjMrNrNCd28btu1UZtD7EtA35PeGle6+ErgM+NUpHE+SnGZblIRlZn9HJGi/Pk6fNxX4HTDH3QdH2eeLQJu73z2GdTwLXOXuR8bqMyQxqctFJApmdgPwMvDXo4V54E7e2Tcf6zrKgW8ozGUkaqGLiCQJtdBFRJKEAl1EJEko0EVEkoQCXUQkSSjQRUSSxP8HKK01zmsm7JAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(abs(kavg))\n",
    "plt.xlabel('Frequency (MHz)')\n",
    "labels=([\"60\", \"80\", \"100\", \"120\"])\n",
    "plt.xticks([0, 10, 20, 30], labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "destroyed-hepatitis",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = k*alldataavg\n",
    "obs=np.nanmean(obs, axis=0)\n",
    "v = []\n",
    "for j in range(flow, fhigh+2, 2):\n",
    "    v.append(j)\n",
    "v=np.array(v)\n",
    "fit = np.polyfit(np.log10(v/100), np.log10(np.abs(obs)), polydim)    \n",
    "power = np.zeros(int((fhigh-flow)/2)+1)\n",
    "i = 0\n",
    "while i <= polydim:\n",
    "    power = power + fit[polydim-i]*(np.log10(v/100))**i\n",
    "    i += 1\n",
    "y = 10**(power)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "confirmed-dependence",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Frequency (MHz)')"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABJr0lEQVR4nO2deXhjd3nvP68k78t4nX2xZ8ksCZnMZLJMFkhIAgkhSaG0CVugQGnKdksvvQV6721p6eW2hdt7KQ2BhrIUylqWQANJIAkkmRmSyTKTZTbP7llsy7LHlmRr/d0/pCPLtiQfLUce2e/nefSMdHR09PMc+7zn3b6vGGNQFEVR5i+u2V6AoiiKMruoIVAURZnnqCFQFEWZ56ghUBRFmeeoIVAURZnneGZ7AfnS0dFhurq6ZnsZiqIoFcWzzz7rNcZ0Znqv4gxBV1cXu3fvnu1lKIqiVBQicjzbexoaUhRFmeeoIVAURZnnqCFQFEWZ56ghUBRFmeeoIVAURZnnqCFQFEWZ56ghUBRFmeeoIVDmFMe8AX5zcGC2l6EoFYUaAmVOcd+vD/NfvvP8bC9DUSoKNQTKnMLrD3NuLIIOXFIU+6ghUOYUQ8EwcQOBcGy2l6IoFYMaAmVOMRQIAzA6HpnllShK5aCGQJlT+IKWIYjO8koUpXJQQ6DMGaKxOOfGEp6AegSKYh81BMqcIZEkTjwfUY9AUWyjhkCZMwwlw0KgoSFFyQc1BMqcwReYCAdpaEhR7KOGQJkz+ALqEShKIThqCETkZhE5ICI9IvLxHPtdJiIxEXmLk+tR5jaTQ0PqESiKXRwzBCLiBv4ZuAXYBLxVRDZl2e/vgIecWosyP7A8gtoqV9k9gmgszhcePZTqY1CUSsJJj+ByoMcYc8QYEwa+A9yRYb8PA/8B9Du4FmUeMBQIU1/tpr2hpuyGYOeRQT778EEeeaWvrN+rKKXASUOwDDiZ9ro3uS2FiCwD3gTcl+tAIvJ+EdktIrsHBlRZUsmMLximtb6aplpP2UNDOw8PAjCoHoFSgThpCCTDtqlKYP8X+HNjTE5hGGPMl40x24wx2zo7O0u1PmWOMRQI09ZQTXNtVdn7CHYkDYEvECrr9ypKKfA4eOxeYEXa6+XA6Sn7bAO+IyIAHcAbRCRqjPmxg+tS5ii+QJjWhmqqXMLpc+Nl+15/KMqLp84BMOhXj0CpPJw0BM8A60SkGzgF3AW8LX0HY0y39VxEvgb8TI2AUii+YJjujgYARvvKFxp65qiPWNzgcYmGhpSKxLHQkDEmCnyIRDXQPuB7xpiXReQeEbnHqe9V5i9DgQitDdU01VaVNVm847CXareLy7raJvUyFMvf/2I/P37+VMmOpyjZcNIjwBjzIPDglG0ZE8PGmHc7uRZlbhOKxvCHorTVVzMWSTw3xpAMOzrKziODbFnZwpKWWn57xFey4373mZNc1tXG72xZNvPOilIE2lmszAmGg4lQkOURxOKGYBmG0wwHw7x8eoSr1nTQ0ViD1x8qyXS0WNzgC4YZ0cY4pQyoIVDmBFZIpq0hUT4K5ZGZ+O1RH8bA9jXttDVUE4rGS2KAfIEwxqhUhlIe1BAocwKro9fqI4DyyEzsPDxIbZWLS1a00NZQDVCSPMFgsgxVpTKUcqCGQJkTWJPJrD4CKM9Mgp2HB7msq41qj4v2pCEoReWQVYaqHoFSDtQQKHOCoYyhIWfvpr3+EAf6Rtm+ph2A9sYaoDRNZV5/4hgj45GS5BwUJRdqCJQ5gTWLoKW+iqakR+D03fSuI4lu4u2rk4Yg6RF4S9BUZh0jEjOEovGij6couVBDoMwJhoJhmms9VLldZUsW7zg8SGONh1ctWwBQ2hyBf8Kr0MohxWnUEChzAl9SZwgoW2ho1+FBLu9uw+NO/BnVV7up8bhKZAh0yI5SPtQQKHOCoWBCZwigodqDiLMX0LPnxjniDXBVMj8AICJ0NNaURG/Im+YRqCFQnEYNgTIn8AXCtNUnDIHLJTTWOCtFvfOIF4ArV7dP2t7WUJ0q/SwGb3K2AsDImIaGFGdRQ6DMCYYCEx4BQLPDekM7egZZUFfFpiXNk7a3NVSXLEfQ1Z4U0FOPQHEYNQTKnMAXnMgRQCJP4GQfwc4jg1y5ug2Xa7KWUXtDddGhIWMMXn9oQklVk8WKw6ghUCqesXCM8Uic1vrJhsCpC+hJX5DeoTGuWtMx7b32xuI9gmDy55kwBOoRKM6ihkCpeCa6iqtS25wMDVljKbevaZ/2XltDDWORGMFw4d9teRQr2+sR0fJRxXnUECgVT7rOkEVTrYfRkDMX0B2HvXQ0VrNuYeO091IyE0WEhwaSFUOdTTXJpLd6BIqzqCFQKp505VELp4bTGGOS+YH2jLMOStFUZjWTdTTUJOcvq0egOIsaAqXiGUqGhlobpuYIoiXX6TnqDdA3EsqYHwBoayyBIUh+tqOpOvVzKIqTqCFQKp6UR1A/2SOIxQ1jkdIOp9mRIz8Aibt4mNwQli+WR2ApqWofgeI0agiUimcoEMYl0Fw3kSx2Sm9o55FBliyopau9PuP7pfAIvP4wTbUeajxu9QiUsqCGQKl4fMEwLfXVuNNq+p3QGzLGsOvwINuz5AcAGqrdVBepN+T1h+hISlo7mfRWFAs1BErF4wuEaa2vmrTNieE0B/v8DAbCWcNCkNAbam+oLmo4zaA/TEfSs2iuc7ZDWlFADYEyB0hXHrVwIjS043BCXyiXIYBEU9lgETkCrz9Ee0OaR+BA0ltR0lFDoFQ8Q4HIpB4CIG04TenCKjsPD7KyrZ7lrZnzAxZtDTVFVw21N1qS2omkdzBc2qS3oqSjhkCpeKbqDEHpPYJY3LDryGBqGlkuigkNRWNxhoLhSTkCUJkJxVnUECgVjTFmmvIolD5ZvO/MCCPj0RnDQlCcAqkvGMYYUjkCJzwbRZmKGgKlohkNRYnGzaQeAij9cBq7+QFI5AiC4RhjBYRzLGmK9qRH0Jw0aE4qqSqKGgKloknpDE3xCCaG05TmArrz8CCrOxtY1Fw7474pvaECBtRYhmAiNGRVP6lHoDiHGgKlopnQGaqa9l6pdHoisThPH/VNGkuZi7ZkxU8h4SGrI9lKFjdrjkApA2oIlIpmKCVBXTPtvVJ15b546hyBcIztqzPrC02lLeURFG4IOhomewSaI1CcRA2BUtH4AokL5NQcAZRuOI01f+DK1W229rcSvb4CpKgHA2Gq3EJzXcITsP5Vj0BxEjUESkUzkSOYHhoqlRT1zsODbFjclErgzkRbUTmCRDOZJWFRV+XG7RL1CBRHUUOgVDS+YOIOurHGM+29UoSGQtEYzxzz2aoWsmis8VDtdhUYGppoJoOEZEVTrYeRMfUIFOdQQ6BUNEOBMK311RlF4EoRGtrbe45QNG6rkcxCRBK9BIWEhvyhaZ6Hk/OXFQXUECgVTiadIQsrNFSMTk/vUBCAtRnGUuai0KYyb5rgnIWT85cVBdQQKBXOUDA8TWfIoqnWQzRuGI/ECz7+1AYvu7Q3VuPN0xAYYyZJUFvoTALFadQQKBXNTB4BFFd6mariqZ2eg8hFe0M1vjyTxYFwjFA0nmpIs2jSucWKw6ghUCqaoWAkY8UQlEaewefPnoPIRVtDTd45gtTQevUIlDKjhkCpWGJxw3AwnLGHAEojPJeQhM4vLASJ0FAgHGM8j5nJU7uKLUrVIa0o2XDUEIjIzSJyQER6ROTjGd6/Q0T2isgLIrJbRK5xcj3K3GJkLELcTNcZspgIDRV+Nz0YCE0L1dihvYDuYu8UnSGL5loP/lCUeFyH0yjO4JghEBE38M/ALcAm4K0ismnKbr8CNhtjLgHeA9zv1HqUuYcvJS8xk0dQRGgoEJ52h24Ha035hIcmEtPTcwTGgD+s4SHFGZz0CC4HeowxR4wxYeA7wB3pOxhj/Gaitq8B0FsexTapruKsoaESJIv92ZPRubAu5vl0F1s5gvaG6TkCUJkJxTmcNATLgJNpr3uT2yYhIm8Skf3Af5LwCqYhIu9Pho52DwwMOLJYpfKYUB51xiMIRWP4Q9GCQkOFKJB6/SGaaz1Ueyb/WarwnOI0MxoCEVmYvFh/UETeIyKXi4gdA5KpzGLaHb8x5kfGmA3A7wB/k+lAxpgvG2O2GWO2dXZ22vhqZT5gKY9myxE0pobTFHYBtS7ihSaL049hB28gPC0/ACo8pzhP1uJoEbke+DjQBjwP9AO1JC7Ya0TkB8DnjDEjWQ7RC6xIe70cOJ3t+4wxvxGRNSLSYYzx5vVTKPMSKxGbrWrI5RIaqz0Fl49aMftCQkNNNR6q3JJKANv7vunNZKAegeI8ubpk3gD8oTHmxNQ3RMQDvBG4CfiPLJ9/BlgnIt3AKeAu4G1TjrMWOGyMMSKyFagGBvP+KZR5yVAgTG2Vi7pqd9Z9iqnBtwzNVMkHO6T0hvLIEXj9YdZlkLKwQlwqPKc4RS5D8GljzLks720xxvw414GNMVER+RDwEOAG/tUY87KI3JN8/z7gd4G7RSQCjAF3mmKEYZR5hS8QyeoNWCT0hgoNDSUu4pmG3tihraEmr9DQoD+UceZBKfohFCUXuQzBr0TkJmPMUPpGEXkd8BUmh30yYox5EHhwyrb70p7/HfB3ea1YUZIMBcNZ8wMWTbWegpuxigkNQcKTsNtHEI3FGQpGplUMQaKhDHSAveIcuZK+XwIeE5FUdlZE3pbcfqvTC1OUmcilM2RRbGioEJ0hi7aG6pQxmQnLc+homm4IaqvcVLtdmixWHCPrb7gx5l9EZBx4NOkF3AncA1xvjDlWpvUpSlaGgmFWttXn3KeptorDA4GCjj/oD9HWkL/OkEU+UtSpruIcpbAqM6E4Rc5bHWPMvyWNwfPACeBqY4wmc5XzAvseQeHlo4XmByAhM+EPRQlFY9R4sie0YaLxLFupqgrPKU6Sq3z0RRJ1/wLUA+0kQkUCGGPMxeVZoqJMJxKLMzoezdpVbJE+nCbfO/vBwPQhMfmQ3lS2ZEFd7u/y565QKibprSgzkcsjeGPZVqEoeTKU0hnKLEFtkT6cJleZaSYG/TOHnnKRkpnwz2wIJpRHM3sEzXXqESjOkcsQnJiplFNERMs9ldlgKJC4O54pdNOcVnqZryHwBcIZq3jsko8CqdcfptrtypqYbqqpon/EX/BaFCUXuaqGHhORD4vIyvSNIlItIq8Vka8D73J2eYqSGSsJm20ojUVTgaWX45GkzlBRoSFLZmLmprLE0PrsiWnNEShOkssjuJmECNy3k93BwyQkJtzAw8A/GmNecHqBipKJoRkkqC0KbcaaSdDODpY3YaeE1Js0BNnQHIHiJLnKR8eBe4F7RaQK6ADGjDHDZVqbomTFN4POkEWhw2lSgnNFGILmOg8el9gKDQ3OEIZqrvMQCMeIxQ1uV2HlrIqSDVsy1MaYiDHmjBoB5XzBmkXQMqMhKEy5M9vYyHxI6Q3Z8AgG/bkH4FgGza/hIcUBdGaxUpH4gmGaaqZr90+l2NBQMcliSHYXz+ARGGMY8IfozCF3nRKe0/CQ4gBqCJSKZCgws84QFB4aSukMFeERQMKjmClZ7A9FCUfjOT2CZjUEioPYMgQiskpEbkw+rxORJmeXpSi58QUjtgxBY01hHoGlM9RUU5jOkEV7Q82MHkFqVnGuHEGBBk1R7GBnQtkfAj8gITYHiQEzP3ZwTYoyI0OBMG31uUtHAdwuobEm/+E0vkCI9oaagnWGLOzkCKx8RCbBOYtCPRtFsYMdj+CDwNXACIAx5hCw0MlFKcpM+GyGhqCwGvxCh9ZPpb2hmtGk3lA2vP6ZK5QmhtNoaEgpPXYMQcgYk7qlSU4n027iOcShvlHe9i+78hqiMhOPHegv6fGmMhQMz1g6alGI8NxgIHcVj12sHIPVCZ35u5IegY1ksfYSKE5gxxD8WkQ+CdSJyE3A94GfOrsspZx88deH2XF4kJ+/dKYkxzvpC/IHX32GLz7eU5LjTWU8EiMYjuXhEVTl7xEEQkX1EFikmspyJIztDMDR0JDiJHYMwZ8DA8CLwB+RmDj2351clFI+Bv0hfrYnYQAeeaWvJMe0jrPjsDOK5Xa7ii2aaj2MhvIsH/UXJ0FtkS48lw2vP8SCuqqcpbDVHhe1VS5GQ2oIlNKTsyRCRFzAXmPMRcC/lGdJSjn57u6ThGNxrl/fyVM9g/hD0VSlTaFYhuCVMyMMB8MzNn3lS0pnyHZoqIqjXvvDacYjMQLhWGlCQym9oeyGYKZmMguVmVCcIqdHYIyJA3umCs8pc4NY3PCtXSe4ak0797xmDeFYnCcODhR1zKFAmKeP+bi8uw1j4LdHfSVabfp3WMqjziSLB0sgL2FhR4HU6w/RYcP7aKr1MDKmHoFSeuyEhpYAL4vIr0TkAevh9MIU5/nVvj5ODY9x9/YuLl3VSkt9VdHhoUf39xOLG/7b69dT43Gx04HwkBVvn2kWgYWVLLarmG6Ve2abDZAPzbVVuF2Ss6nM6w/R0WTPI9CGMsUJ7MQAPuX4KpRZ4Rs7j7N0QS03blyIx+3itRsW8uj+fqKxOB53YU3nD79ylsXNtVy6qpVtXa3sOlJ6QzCUZ2ioubaKSMwQisaprZp5JoE3ZWiK9whcLplxiP1MgnMWzSpFrTjEjH/txphfZ3qUY3GKc/T0+3myx8vbr1yVuui/btMihoMRnjk2VNAxxyMxfnPQy02bFiEibF/dzv6zoyUvI/UFI4jAgjr7HgHYl2fw2ajrz4f2HHpDkVic4WDEVo6gWXMEikPY6SweFZGR5GNcRGIiMlKOxSnO8c1dx6l2u7jzshWpbdeu66Ta4+KX+woLDz15yMtYJMbrLlwEwPY17QD8tsRewVAgzIK6KtteS74KpCnBuRIkiyHZXZzFEFjbc/UQWDTV5t8hrSh2sOMRNBljmpOPWuB3gS84vzTFKfyhKD94tpdbL14y6QLUUOPh6jXtPPJKn+14ejoPv3KWploPV3QnDMDFy1uor3azs8SGwJdHMxkkxjyCfUPgDYSodruKrp6yyGUIUvIStqqG8m+MUxQ75B0INsb8GHht6ZeilIsfPX8KfyjKO7evmvbeTZsWc8IX5GBffvNxY3HDr/b1c/36hal6+Cq3i21dbSVPGNtVHrXItyvXlyznLFZnyKKjsSZ1wZ/KYB6J6ebaKsYjcSKxeEnWpSgWdkJDb057vEVE/jcqMVGxGGP4xo5jvGrZArasaJn2/o0bEzJSj7xyNq/jPndiiMFAOBUWsti+up1D/X4GRmee22sXXyBsO1EM+Xfl+gKl0RmyaGuoZnQ8ITU9FasCyk4+otAhO4oyE3Y8gtvSHq8HRoE7nFyU4hy7jvg41O/nndtXZbzjXdhcy+YVLTyyrz+v4z788lmq3S5ec0HnpO1WnqCU1UNDwbDt0lHI3yPwOmAIYKIjetJ3jSZzBDmURy0sg6bCc0qpsWMI7jfG/EHy8YfGmL8F1jm9MMUZvrHzGC31Vdy+eWnWfV63aRF7Tg7TNzJu65jGGB5+pY+r1ranLlYWFy1tprHGU7I8gTGGoYC9WQQW+Wr5+wIhW8lbu6SayjKUkFr5CDtzD9QjUJzCjiH4J5vblPOcM+fGePiVPu7ctiJnPf1NmxLhHbvVQ4f6/RwfDKY+l47H7eLy7jZ2lShPEAjHCMfieSWLG1Plo/YuoKWSoLaw4v+ZhOcG88hHTIS41CNQSkvW2xAR2Q5cBXSKyJ+mvdUMzNyVo5x3/PtvTxA3hndcOT1JnM66hY2sbKvnkVf6ePsVufeFRFgI4KaN0w0BJPIEj+7vp29knEXNtfkvPA2rmSyfC7U1nMbOBdRSNi1V6Sjk1hvy+u17H811+Rk0RbFLLo+gGmgkYSya0h4jwFucX5pSSkLRGN9++gSvXb+QFW31OfcVEW7atIgdPYMEbKhdPvxKH1tWtrAwy0XeyhOUonrIV4AhAPt6Q6XUGbLIFRqyKzgH6SEu9QiU0pLVECQ7iD8FXGmM+VTa4/8kp5QpFcQvXjqL1x/m7qu6bO1/06ZFhGNxfjODCN2Zc2Ps7T2XMSxksXFJM821ntIYgmTCNZ8cAdivwR/0W/ISpcsRLKiz9IYyGYKQLXkJSO+QVo9AKS12OmaCIvIPwIVA6pbPGKO9BBXE13cco6u9nmvXdtjaf1uaCN0tr1qSdb9fJkXqXrdpcdZ93C7hitXtJUkYp0JDeUpb2x1OM1jirmJI6A211ldPyxEYY/D6w7YE54BUg5t6BEqpsZMs/hawH+gmIUB3DHjGwTUpJealU+d47sQw79zehctlr0nK43bx2vULefRAQoQuGw+/0sfqzgbWLmzMebztq9s54Qtyangsr7VPJTWLwKnQUIl1hizaMwjPjYaihGNxWxLUkDgnDdVurRpSSo4dQ9BujPkKEEmGi94DXOnwupQS8o2dx6ircvOWS5fn9bmbkiJ0u49nFqE7NxZh5+HBnGEhi1LlCYaCYdwuobk2P/kHu0NdLLnoUkhQp5NJZmKiqzi/5jj1CJRSY8cQWL91Z0TkVhHZAuR3RVFmjeFgmJ+8cJrf2bLMtlqnxasv6KTa7co6o+DxA/1E4yZnWMhi/aImWuurijYEvkCE1vr85R/ySRZXexJ33qWkrTGTIcjf6OhwGsUJ7BiCT4vIAuC/Ah8D7gc+aufgInKziBwQkR4R+XiG998uInuTjx0isjmv1Ssz8r3dJwlF49ydQVdoJhpqPFy1NrsI3cOv9NHRWJNRqmIqLpdw5ep2dh0ZLEjQzmIokF9XsUU+oaH2htLpDFlkkqLOR3DOopD5y4oyEzkNgYi4gXXGmHPGmJeMMdcbYy41xsw4oSz52X8GbgE2AW8VkU1TdjsKvMYYczHwN8CXC/oplIzE4oZ/23Wcy7va2LikuaBj3LRpUUYRulA0xuP7+7lp00LbeYfta9o5NTzGSV/heQJfMD+dIYvm2irCsTjjkVju4wfsl3PmQ3tDDefGIpME47x++xLUFs119pLeipIPM80sjgG3F3jsy4EeY8wRY0wY+A5TNIqMMTuMMVYAehcaciopvzk4wEnfGHdflb83YHHjxsxdxjsPDxIIx2yFhSy2r07mCY54C17PUIE6QHblGQb9oZKWjlq0JY3LUJpXYOUI8vl57FY/KUo+2AkN7RCRL4jItSKy1XrY+Nwy4GTa697ktmy8F/h5pjdE5P0isltEdg8MFDdcfT7xzDEfHpfkdbGeyqLmWjYvX8DDU/IED7/SR0O1O5UEtsPahY10NFYXlScYCuYnQW1hV3huMBCmo8QVQ5B5iL3XH6KlvoqqPMaCJnIEGhpSSoud0ourkv/+ddo2w8wzCTLFCzIGh0XkehKG4JpM7xtjvkwybLRt2zaVwLZJ79AYS1vqUvMBCuWmTYv47MMHUxIR8bjhkVf6uG79QlszgC1EJvoJjDF5x+HjccNQMJJ3DwHYH05Taglqi0wyE4OBUN5lqnZzHYqSD3YmlF2f4WGnmawXWJH2ejlweupOInIxiQT0HcaY0k86n8ecHAqyvLWu6OPclPQorPDQC73DDIyGbJWNTmX76nb6RkIc9Qby/uzoeJRY3BTpEWS/iI6FEzpDbQ7kCKyEcPqAGq8/nHeZqt1ch6Lkg53BNItE5Csi8vPk600i8l4bx34GWCci3SJSDdwFTEoyi8hK4IfAO40xB/NfvpKL3qGxkhiCCxYlROisLuJHXunD4xKuX78w72Ol+gkK6DK2OnMLqxqaWafHOr7dBq98sPIOvimhoc68DYFKUSulx07M4GvAQ4AlYH8Q+JOZPmSMiQIfSn52H/A9Y8zLInKPiNyT3O1/Au3AvSLygojszm/5SjbGIzEGRkOsaM0tMGcHEeHGjYt46nBChO7hl89y5ep2FtTnf0Fe3dHAwqaagvIE1mCXQqqG7HgEhSRv7dJSV4VLpoSG8hCcs1ApasUJ7BiCDmPM94A4pC7wtvxSY8yDxpgLjDFrkgNtMMbcZ4y5L/n8fcaYVmPMJcnHtgJ/DmUKvUOJEs3lbcV7BJAUoYvG+dqOYxweCBQUFoKEUdm+pp1dR3x59xP4AomLXyEXaku5cyTHBdTngM6QxYTeUOI7wtE458YitgXnLFR4TnECO4YgICLtJBO9InIlcM7RVSlF0zsUBGB5CTwCgMu6WllQV8Xnf5UQni3UEEAiT+D1h+jp98+8cxpW6WUhHkGjHY8gJUFd+tAQJAyM1U1seTfqESjnA3YMwZ+SiO2vEZGngG8AH3Z0VUrRWB5BKUJDkBSh27CQUDTOq5YtYGlL4Z5GoXkCS4K6EI/A7ZIZBdtSEtQOeAQwWW9oYNTqKs4zR1CnOQKl9NipGnoOeA2JMtI/Ai40xux1emFKcZwcClLlFhbaGIpuF8sLKMYbAFjZVs/SBbV55wmGkjpA9QXqAM0k2OYLhKlxQGfIor2hJuV1WP/mIy8B6hEozjBjH4GI1AIfIFHjb4AnROQ+Y4y9yebKrNA7NMayljrb8g92uGHjQv74ujW87YqVRR1HRLhyTTuP7e8nHje21+gLhGkrQHDOYqYafK9DOkMW6R5BIYJzkJYjUOE5pYTYCQ19g8RQmn8CvkBCN+jfnFyUUjy9Q2MzjqTMlxqPmz+/eUPe4YxMbF/dzlAwwoG+UdufKbSr2GImwTZfIFRy+el02hurGQ4m9IYKEZwDaKz2IKIegVJa7BiC9caY9xpjHks+3g9c4PTClOLo9ZWmmcwpCplP4CtQedRiJp0ep7qKLawu4qFgmEF/IsxlTR2zi8slNNZ4zouqoZO+IG+69yn6RzU4UOnYMQTPJyuFABCRK4CnnFuSUizBcJTBQLhkFUNOsLy1nhVtdXkljIeCkYIqhizshoacIr2pzOtPaBoVEoZqPk+E557s8fL8iWGePZZ5cJFSOdi5HbkCuFtETiRfrwT2iciLgElKSCvnEaesHoLz2COARHjoFy+dJRY3uG3kCYq9Y7eTLHaih8AipTfkD+P1h+goMJGfMGizHxo6lJQmP1KAXIhyfmHHENzs+CqUknKyxD0ETrF9TTvf293LvjMjXLRsQc59o7FEA1YxhqC5NntIJRiOMhaJOSJBbWEZGW8gzGAgf3kJi6ZaT87GuHJxqD+R3zmmhqDisVM+ehwYARaQkINoJzHH+HjyPeU8Y6KH4Hz3CDqARIhhJobHCu8qtmiq9RCOxglFpzfGFzI/OF/aUx5BKCkvUZghOF9CQ1ZD4LFBNQSVjp3y0b8B3g0cZkJG2o4MtTJLnPQFqfG46CxhD4ETLF6QmHXwuYcP4HEJ772mO2vMvJiuYouJGvwoNY2TewVS8hIO5gha6qsRSfQQFKIzZNFU6+FQ/+wagpHxCGfOjSNCQUqyyvmFnWTx7wNrjDHX5SlDrcwSvUNjLGutc6wevpR8/T2Xc936hXz6P/fxh9/YPWmCVzrWhbpYjwAyd+VOKJs6ZwjcSb2hY4NBwrF4EaGh3LmOcmB5A9tWteL1h8+LUJVSOHYMwUtAi8PrUEpIQn76/M4PWLTUV/Pld17K/3zjJn59cIBbP/8Eu4/5pu1XjPKoRa6u3MEC5gcXQltDNYeSvRPFeAQj49G8RftKSU8yUWx1mWueoLKxYwg+Q6KE9CERecB6OL0wpXBODgXP+/xAOiLCe67p5j/++Co8bhd3fnkX9z7eQzw+caErRnnUIpdHUAqPww7tDdUcHvAnnxfuEcTihrFZHE5zqH+UGo+La9d1AhoeqnTsVA19Hfg74EWSUtTK+cvoeIThYKRiPIJ0Ll7ews8+cg2f+OGL/P0vDrDz8CD/eOcldDTWpDyClgJmIFjkmls8mNQZKlTHyC7tjdVEYgkDV6j3kS48V1+dX0NaqTjU72dNZyPdHQ2aJ5gD2Pkt8hpjPu/4SpSScGq4MnoIstFcW8UX3rqFq9a086mfvsIt/+8J/t9dl+ALhGmoduc1IznTsSGzlv+gP0xHY43jeZV0jyNfeQmL9BDXoubakqwrXw71+dnW1UptlZulC+o0NFTh2DEEz4rIZ0hIUacGriZVSZUS8uShRBnlNes6Cj7GSV+ydLTEOkPlRER4+xWr2LKilQ99+znecf9v6WyqKUpnCGZOFjsdFgIm9SkU+vNYP8e5WRKeC4SinBoe422LEuKD3R0N6hFUOHZyBFuAK4H/BXwu+fisk4uar/yvB/fx6f98pahjTAykqUyPIJ1NS5v56Yeu4XcuWUbfSPEXakvXJ1NoyOmuYgvLC2ipr6LKbefPbzrNOUJc5cCqGFq7sBGYMASzmbyuVM4FI7z1y7vY2zs8q+uY0SMwxlxfjoXMd4wxHPUGiMbjhKNxqj2FXSR6h8aoq3I7Wg9fThpqPPyfOy/hdRcuZkFd4fkBSAzXqa92Z5RwHvSHUxc2J7GMWTHVSc1p/RCzwaGkIViX/P/q6mhgZDyaNKbnd+/K+cZXnjrKziOD/OfeM1y8vGXW1jHj1UZEFonIV0Tk58nXm0Tkvc4vbX7RNxJiLBIjEjMc8eY3wjGdk0nV0UroIciHmy9anFIsLYZsOj2DgVBZjKdlCIr5rqZZNwSjVLtdrEyGH1d3NADaYZwvI+MRvvrUUQCeOzG7wn12bju/BjwELE2+Pgj8iUPrmbekX/wPnLWv0T+VRA9B5YeFnCKTFHUwHGU8Ei/L3axVMlqMR5Cr+qkc9PT5Wd3ZgCcZ2upKGoIjA2oI8uEbO44xOh5l++p29vaeIxydvaLMrIZARKywUYcx5nskS0eNMVFg9gqY5yjpybb9RRmCYEWWjpaLTMNprGayciSLrTxEoRVDAPXVbtwumbVu3kP9/klhtOWtdXhcoh5BHgRCUb7y5FFeu2Ehb79yJaFonH1nRmZtPbk8gqeT/wZEpJ2kzlByNsE5pxc23zjmDVDjcXHBokb2F/gLcW4swsh4lBVt6hFkI5NHUOj84EJora+mvaGatYuaCj6GiMw4W8EpxsIxTg4FuSBt/VVuFyva6rVyKA++ues4Q8EIH37tWraubAXg+VkMD+VKFltB5j8lUTq6RkSeAjqBtzi9sPnGUW+A7o4G1i9u4pmj0yUW7NBbIfLTs0lTrYeTvuCkbb6UzpDzoSG3S3jiz6+nxlNc49psGYLDA36MmUgUWyQqh4JZPqWkMxaO8S9PHOHadR1sSRqBxc21PHdimHdfPTtrymUIOkXkT5PPfwQ8SMI4hIAbgb0Or21eccQbYP2iJjYsbuYnL5zmXDDCgjy7aHsrZCDNbNKcIVns9TuvPJpOKbqBm2pmR3jOmkGwbtFkQ9DV3sDOw4MYY+ZcoUKp+fenT+D1h/nIDetS27auapnVhHGu0JAbaASagAYSRsMN1Ce3KSUiGotzYjBIV0cDGxYn/mvzGepuYd3prlCPICtNtVXTOotTEtRlCA2ViqZaT8YyWKc51OfH4xJWtTdM2t7d2cBYJEbfSCjLJxWA8UiML/36MFeubuOyrrbU9q0rW+kdGpu1+c+5bk3OGGP+umwrmcecGh4jGjd0dzSwYUnSEJwd4fLuthk+OZneoTEaqt1F6fHMdZpqJobTWOGZQX+I2irXrOn2FEJTbVUqFFhODvX76e5omNYM1500DEe8fhYvmB3Zi0rg+7tP0j8a4v/eecmk7VtWtgDw3PFhbr5ocdnXlcsjUP+uTFgzX1d3NLC4uZbmWg/7CqgcsuSn1TXPTiaZicFAuGAl0NmiuW52cgQ9/f5JiWKL7s5kL4HmCbISjsb54uOH2baqdVpPzIVLF1DlFp4/OTvhoVyG4IayrWKeczRZf51QchQ2LG4uqJegdyioFUMzkKkZq1zyEqWkeRaG04xHYhwfDGTswF7SXEuNx8XRIpoh5zr/8Vwvp8+N8+Eb1k27WautcnPh0gU8f3x4VtaW1RAYYworXVHy5qg3QHOtJ1XHvmFJEwfOjual3WKMqaiBNLNFpmasQX+4LD0EpaSp1oM/FJ00s8FpjgwEiJvpiWIAl0tY1V6vlUNZiMTi3Pt4D5uXL+DVWUQlt65sZe+pYSKx8jeWFSZoo5QUq3TUuktYv7gJfyiaqgKyw7mxCP5QVCuGZiCrR1BhoaGmWg9xA4Fw+cJDqYqhhZlrRRIlpOoRZOInL5zmpG+MD792ujdgsXVVC+OR2WksU0OQAV+WublOYRkCiw2Lm4H8pCYs+Wn1CHIz1SMwxuD1hyoyNATl1Rvq6ffjdgldHZl/x7o6GjjhCxIro5dSCcTihnsf62HTkmZu2Lgw635WT8Fzx8ufJ1BDMIUDZ0fZ9ulH+NW+vrJ833gkxulzY3R3TLjb65MlpPvP2r8zmEvy004ydThNMBwjFI1XnFrrbAjPHerzs6q9Pmsz3OqOBiIxw6k8PNn5wM/2nuaIN8BHblibs5Bj6YJaFjXX8PzJ4fItLokagik8/PJZ4ga+v7u3LN93fDCIMRNVF5DQzV/eWpeX5pAVRtIegtxMrRoq16ziUjMbwnMH+0e5IEtYCBJNZQBHVXMoRTxu+MKjPaxf1MTrNuUuCxURtq5snZXGMjUEU3jsQD8Ajx7oL4uolxVT7Z7SoLNhcXNehuDkUJCmWk/e3cjzjcYpF1CvP9EAVWmhIcsQlEt4LhSNcXwwmDFRbGHdzBwd0DyBxS9ePsuhfj8ffO1aXK6Zy7q3rmzlpG+MgdHyNuapIUjDFwjz/MlhXn1BJ+FonEdedj48ZFVZTI27bljcxFFvgFDUntCrVgzZo8rtoq7KPc0jqLRkcXNdeUNDx7yJ2H+u4T2djTU0VLs5NqiVQ5DIP/3Toz2s7mzg1lctsfWZratagPLPJ1BDkMYThwYwBj564zqWtdTxwJ7Tjn/nUa+fzqaaVMzXYsOSJmJxkxoLOBMJ+WnND9ghfTjNYIWHhqbKZTjFTBVDkAhtdHc2pBokC2H3MR9nzs2NHMMv9/Wz78wIH7xuLW4b3gBMNJbNKUMgIjeLyAER6RGRj2d4f4OI7BSRkIh8zMm12OGx/f20N1SzeXkLt21eypM9XscriKZWDFlYmkP7z8wcHjLGcNI3pvkBm6Qrd1qzCCotNDRRNVSe0NChPj8ugdWd039X0+lqb+BYgYZgOBjmzi/v4vX/+BseevlsQccolPFIjLPnSqfzk/AGDrGyrZ47Llk68weS1Fa52TQLjWWOGQIRcQP/DNwCbALeKiKbpuzmAz4CfNapddglFjf8+uAAr7mgE5dLuH3zUmJxw4MvnnH0e496A9PyA5D4g6r2uGyJz/kCYcYiMfUIbJI+k8AXCFFX5a4onSGAGo+LKreUTXiup9/PyrZ6aqtyy2ev7migdyhY0LStHYcHicUNzXVV/NG/Pcunfvqy7dBosXzmwX3c8LnHOTVcGm/koZf72Nt7jg9evyY1yc0uW1e2lL2xzEmP4HKgxxhzxBgTBr4D3JG+gzGm3xjzDDA7o5bS2NM7zFAwwnUbEnW+G5c0sXZhIz91MDw0Mh7B6w9Pqhiy8LhdrFvYaCthrPLT+TEpNFSBXcVgDacpn8zEwb5R1tkYptPV0UDcwAlf/nmCJw55aazx8PBHX80fXN3FV586xlu+uJMTDuccIrE4D+w5TSAc4y9/8nLRx/OHovzVAy+zYXETb966PO/Pb13ZyngkbisaUCqcNATLgJNpr3uT285LHj8wgEtItX+LCLddvJSnj/lK6jKmY7nQmUJDkOgnsDOt7GSyh2BFm4aG7NCc5hEMBsJlmUzmBM1lGk4TicU56g1MG0aTCet3uZBpZU/2DHDl6nbqqz385W0X8qV3XsrxwQC3fv4JRz3zJw95GQpGuHZdB7/c11d0WOpzDx+gb3Sc//XmV01TabXD1lXJxrIy5gmcNASZsiMFtRyKyPtFZLeI7B4YGChyWZl5/EA/W1e20lI/cVG4bfMSjEk0hDjB0TTV0UxsXNxM/2hoxjyF5REsU4/AFk21nlSS1ReoTI8AKJtHcHwwQDRucpaOWliGIN88wYnBICd9Y1ybpsPz+gsX8+B/uZY1Cxv5wLee43/8+CXGI6UPFT2w5zQL6qr40jsvZcPiJv7qgZfxhwozsHt7h/n6jmO844pVqRGU+WI1ls0VQ9ALrEh7vRwo6IpqjPmyMWabMWZbZ2dnSRaXzsBoiL2957h+w+T279WdjVy0rNmx8NCRgQAisLI985283Q7j3qEgC+qqUglEJTeTQ0Mh2hsrq3TUolzjKg/1JSrXclUMWbTUV9NaX5V35dATPYkbvKvXThZkW95az/fv2c77X72af9t1nDffu6Oks5HHwjEefvkst1y0mPpqD3/7pldxdmScf3zkYN7HisbifPJHL9LeWMOf3by+4DWJCFtWlLexzElD8AywTkS6RaQauIvE7OPzjl8fTPwSvuaC6Ubm9s1L2dN7ruBKiFwc9QZY1lKXtWU/Na1shjzBSd+Yyk/nQVNtFaFonHA0npxFUKkegacsDWWH+v2IwJrOmT0CSOQJ8v17efKQlyULalmTIV9W5XbxyTds5F/fvY3T58Z44+ef4CcvnMrr+Nl4dH8/gXCM2zcnKnsuXdXKWy9fyVefOspLp87ldayv7zzOS6dG+KvbLiz6pmzrqpayNpY5ZgiMMVHgQ8BDwD7ge8aYl0XkHhG5B0BEFotIL/CnwH8XkV4RaXZqTdl47EA/C5tquHDp9K++9eLEL4gT4aFjg5lLRy06m2poa6ie0RD0DgVZ3qL5AbtYNfh9I+OEovGKDQ2l5zqc5FC/nxWt9dRV564YskiokNo3BLG4YcfhQa5Z25FTi+e1Gxbx4EeuZeOSZv7Ld17g3sd7bH9HNh7Yc4qFTTVcsXpiUMyfv34DbQ3V/MWPX7ItoHd6eIzPPXyA69d38oZXFT9hzAorPV8mr8DRPgJjzIPGmAuMMWuMMX+b3HafMea+5POzxpjlxphmY0xL8nlZNVijsTi/OTjAdes7M/4SLmup47Ku1pI3lxljODoQyJofgISLuH5RU85pZRNzCNQjsIt1t3YsqYlTuaGhMhmCvlFbiWKL7vYGzo6MMxa2F89/6dQ5zo1FuCaLTn86S1vq+Pb7r+TGjYu497HDReVIRsYjPHZggFsvXjKp4WtBfRX//dZN7Dk5zL//9ritY/3lAy8TN4a/vuOikkwIvGiZ1Vg2XPSx7DDvO4ufOzHM6HiU69dnl4e9ffNSDvb581IDnQmvP8xoKJrTI4BEh/HBs6NZB5AM+EOEonGtGMoDyyOwwheVHBryh6KOyj5HY3GODARYayNRbJEaW2lTfO7JHi8wPT+QjSq3i4/csBZ/KMp3nzk58wey8NBLZwlH46mwUDp3XLKUa9Z28Pe/OED/SO6qwYdePssjr/Tx0RsvKNnfodVYVq48wbw3BI8d6MfjEq7OcTdyy6sSdwylTBpbrnPXTIZgcRNjkVjWumztIcifppRHkPg/rbSuYgvLoPkd9ApO+IKEY3FbiWKLlAqpzfDQk4e8bFzSTEcentnFy1u4vLuNrz51jGiBjVcP7DnNyrZ6LlnRMu09EeFvfuciQrE4f/2zV7Iewx+K8pc/SfQMvOea7oLWkY0tK1rY21uexjI1BPv72dbVmjO509FYw1Vr2vnpnjN5jY/MxbFU6WjuO631ySE12RrLJgyBegR2meoRVGyOoM6areBcwvhQv1UxZN8j6Mqjl2AsHOPZ40OTykbt8r5rujk1PMZDBYhDev0hdhwe5LbNS7KGcro7GvjgdWv52d4zqYKSqRTbM5CLravK11g2rw3BmXNj7D87mjMsZHHb5qWc8AXZ05tfJUE2jngDVLllxtr/CxY1IpK9cuikTwfS5Mu0HEGFKY9aNE+ZreAEluhhLtXRqTTWeFjYVGPLEPz26CDhWNx2WCidGzYuoqu9nvufPJL3Zx988QyxuOH2zbl7XO+5bjWrOxoy9jCUomcgF1tXtgDw/Ennw0Pz2hA8fiBh5af2D2Ti9Rcuptrt4oEXShMeOur1s6q9YUZVwvpqD6va6rPmJ3qHxmhrqKahprK0cmYTyyM46RujvtptuxrmfKOpDMJzh/pGWdZSl/fvl90S0qd6vFS7XVze1Zb32twu4T3XdPP8iWGezXO84wMvnGb9oqZUr042ajxuPv2mizjhC/KFRyeqlErVM5CLZS11LGyqKcvoynltCB7b38+yljpbbu+Cuipes76Tn+09XZLk3FFvIBVLnYkNi5uzegQqP50/1nCacKxyS0ehPFLUB/v8tjqKp7LaZgnpE4e8bOtqLdgYv+XS5Syoq+IreXgFvUNBdh8f4nabqqBXrengzVuW8aXfHKYnKcddyp6BbExMLBt25PjpzFtDEI7GearHm7VsNBO3b15K/2iIp4/6ivrueNxwbDA4o6SvxfrFTRwdDGQsx+sdUvnpfLGG00Dllo6C8x5BLG44PODPKz9g0dXRwGAgzLmx7GsbGA2x/+xoQWEhi/pqD2+7YiW/eOlsKkw6Ez/dk9Atuu1i+/LQn7x1I/XVHj75o5c4VeKegVxsXdXCCV8wNUnPKeatIdh9zEcgHLOVH7C4YeNC6qvdRfcUnD43Rjgan7F01GLD4iaMmRgOYhGPJwaFq0eQP9bddKWWjoLzOYLeoSChaH4VQxZ2NId2HE6UjRaSKE7nXdu7cInw1aeO2dr/gT2n2bKyJau0SyY6Gmv4xC0bePqoj9+/b2dJewZysSWZe3A6PDRvDcFjB/qpdru4am37zDsnqa/2cOPGRfz8pTNFlXQdnUF1dCoblmSuHBrwhwjH4moICsAyBJUdGnLWI7A0hvLpIbBIGYIcvQRPHPLSUl/FhUsXFLbAJIsX1HLb5qV895kTM1ZQ9fSPsu/MSMbegZn4/W0r2LaqlVPDYyXtGcjFq5YtwOMSnj857Oj3zGNDMMAVq9vyHkhy++alDAcjPHnIW/B352sIEgNBXNPKyFIVQ9pMljfWRbRSewgAqj0uajwuxzyCQwVUDFmsbKtHJCGsmAljDE8e8nLVmnbbYxxz8d5rugmEY3z36dwNZg+8cBqXwK0X25shnI7LJfzjnZfw325eX/KegWzUVrm5cGmzegROcNIXpKffz3V5hIUsrr2gg+ZaT1HNZUcGAtRXu1nYZC8+7XYlpCYO9E2uHLJ6CFaoR5A3cyE0BAmD5lQfwaH+UZYsqC0oGVpb5WbpgrqsHsHhgQBnR8a5Zm1p1IQvWraA7avb+epTR7M2mBljeGDPabavaWdhU21B37OirZ4PXLe25D0DudiyspW9vecKbpyzw7w0BI8f6Afg+vX5/xLWeNzcctESHnr5bMHa6JbYXD7xxfWLm6ZVDvUmB9IsU8G5vLEubpXaQ2DRXOdxrGroUJ+/IG/AYnVn9sqhJw8lSreLzQ+k875ruzl9bpyfv5R5sMyLp85xbDBYUFhoNtm6qpWxSMzWtMJCmZeG4LEDA6xqr7cdmpnKbZuXEgjHeHR/f0GfzzawPhfrFzfj9YcnydKe9I3R0VhTsXXws0kqR1DBoSFwTnguHjf09PsLShRbdLUnDEGmbvwne7ysbKsvaZz9+vULWd3RwP1PHMn4nQ+8cJoqt3DzhfmHhWaTLUkJDCd1h+adIRiPxNhx2Mv16xcWnPHfvqadjsaagsJD4Wick75gTtXRTGzMMJugd1h7CAplroSGmtOG7JSSU8NjjEViBfUQWHR3NDA6HmVwyoS9SCzOriM+W2qj+eBKNpjt6T03rcEsHjf8bO8ZXnPBQhbUV9YAp+WtdXQ21fC8g/0E884Q7DoyyHgkznUFhIUs3C7h1lct5lf7+/P+IzzhCxI3M4vNTSXTtDKVny6ciWRxZYeGmmo9jOSo1S+UngI0hqaSrYR0z8lh/KEo1xbRP5CN3926nJb6Ku5/4uik7U8f83F2ZNx2E9n5RKKxrEU9glLy+IEBaqtcXLnaftloJm6/ZCnhaDxvwat8K4Ys2htr6GisScUJY3HD6eExlZ8ukK0rW9mysoXOCjcETg2nsXpWigkNWb/jU8dWPnHIi0jCsy41ddVu3nHFKh565SzH0xLVD+w5TV2Vmxs35l8gcj6wdWUrxwedayybV4bAGMOj+/u5ak0HtVXFxdW3rmylu6OBb9kcXGFxrEBDALBxyUTCuG9knEjMqEdQINes6+BHH7iaak9l/wnkM7d4LByzPULyYJ+fhU01RYVRlrfW4XHJtO98qsfLxcsW0FLvTFju7u2r8LgmGswisTg/f/EMN21alHe5+PnC1lXWxLJhR45f2X8FeXLUG+CEL1hQtdBURIR3XrmK508Ms7d32PbnjngDtDVUF/RHsH5REwf7RonFjcpPK0AixDUWic3Y4BiJxXnXV5/mus8+zh1feJJ//+2JnGHNQ/2FaQyl43G7WNlWP6lyaHQ8wvMnh0ueH0hnYXMtt29exvd2n+TcWKLnZygYqbhqoXSsxjKnwkPzyhBYaqOF9A9k4i3bllNf7eYbO+17BUe9frryaG1PZ8OSZkLROMcGA6lmMu0hmN/YHU7zdz/fz9NHfbzjypWMRxLKmZf/7a/42Pf38Mwx36QqG2MMPX2jRYWFLLqmiM/tOuIjFjdF6QvZ4b3XdBMMx/jO0yd4YM9pFtRV8eoLStOzMBvUVrn53O9v5s1bcstmF0pl+kkF8tiBftYubCxZXL25too3b13G93b38sk3bLQlV3DUGyi4iWaDlTA+M5ryCJa2qCGYz1hJ75HxCK1Zfv9+tvc09z95lHdf1cVf3X4hxhj29J7ju8+c5Kd7TvODZ3tZ3dnAndtW8Oaty4nE4gTCsaJ6CCy6OxrYcdhLPG5wuYSnerzUVbm5dFXp9fvT2bS0mavXtvPVp44xOh7hts1LKz4MeMclzhgBmEceQTAc5bdHfCUJC6Vz9/YuwtG4rdmpgVCUvpGQbdXRqaxd2IhL4MDZEXqHgixsqik616FUNjMJzx3sG+W//WAvl65q5ZNv2AgkwpqXrGjhM29+FU//xQ38w1supq2+ms/8fD/bP/Mr/vibzwLFVQxZdHU0MB6J0zeamPv7xKEBLu9uo8bj/O/t+65ZzdmRcQLhWEWHhcrBvPEIdvQkJiHlozZqhwsWNbF9dTvf3HWc9796dU7dFKvdvtBGttoqN90dDew/O8rIeEQrhpRJHsFURscj3PNvz1Jf7eHet2/NeEdcX+3h97at4Pe2raCn38/3d5/kP57rpbbKxYbkmNRisPpljiY1hw4PBLjrspVFH9cOr7mgkzWdiV6GK4qsEpzrzBtD0N3ZwIeuX8u2AiYhzcS7rlrFPd98jl/u6+P1F2bXJ08NrLc5kCYTGxY38+Kpc8SNcdy9Vs5/mrJ4BMYYPvb9PRz3Bfn3913BouaZtXXWLmzkE2/YyMdev57hYKQkjVep+cWDAU4NJ8KZTiaK03G5hC+9cxvjkVhJhO3mMvMmNLSms5GPvX69I3HCGzcuYumCWr6x81jO/ay7oq6Owu/kNyxu4oQvyOlhbSZTJjSTpjaV3ffrIzz0ch+fuGVD3nfDVW4XnTYFEWdiSXMtNR4XRwcCPNnjpaOxOpXrKgdrFzZy0bLiZK7nA/PGEDiJx+3i7Veu4qmewdQou0wc9QZYsqC2qFpmq8M4btDJZArNddM9gqd6vPzDQ/t548VLeG+Z5JKz4XIJXe0NHPEGeKrHy9VrOxwf5qLkjxqCEnHXZSuodrtylpIeHcxfbG4q6XFb7SFQGmsmG4JTw2N8+NvPs3ZhI3/3uxefFxfdro56dhz24vWHucbhslGlMNQQlIj2xhreuHkJ//Fsb9ZGnUJUR6eyvLWOhqTaqIaGFI/bRX21m9HxCOORGB/45rOEo3Hue8elNNScHynA7o5GxiOJhrdr11VuLf9cRg1BCXnX9i4C4Rg/fO7UtPeGAmGGg5GiDYHLJVywuAkR7SFQEjTVehgZj/Cpn77Cnt5zfO73N7O6s/jSz1LRncyJrV3YyOIFhQ2EUZzl/LhlmCNsXtHC5hUtfH3nMe7evmqSW36kCI2hqVzR3c54JF7xDTJKaWiureKRV/oYCkb4wHVrclauzQbdHQmjpGGh8xe9kpSYd21fxZFkhUQ6xYjNTeXPXr+en3zw6qKPo8wNmmo9DAUjXLO2g//6uvWzvZxpbFrazNaVLbx5q3OdsUpxqCEoMbdevIT2hmq+vmNy0vioN4DbJSVpAnO7RL0BJcWSljqWtdTx+bduOS/r5RtrPPzwA1dz8fKW2V6KkgUNDZWYGo+bt16+kn9+vIeTvmDqwn/UG2BFa11Zh14r84O//92LicYNC+oqa/KWcv6gVyUHeNsVK3GJ8M20WQVHSlAxpCiZaKjxqBFQikINgQMsbanjdZsW8d1nTjIeiWGM4Zg3kEqaKYqinE+oIXCIu7d3MRyM8MCe0/SNhBiLxOguUHVUURTFSTRH4BBXrm7jgkWNfH3HMZYn6/27ixCbUxRFcQr1CBxCRLh7excvnx7hB8/1AqhHoCjKeYkaAgd505ZlNNV6+OFzp6jxuFhiQwpYURSl3DhqCETkZhE5ICI9IvLxDO+LiHw++f5eEdnq5HrKTUONh9+7dAWQaCRznYc13oqiKI4ZAhFxA/8M3AJsAt4qIpum7HYLsC75eD/wRafWM1u8c/sqoLhhNIqiKE7iZLL4cqDHGHMEQES+A9wBvJK2zx3AN4wxBtglIi0issQYc8bBdZWV7o4G/ucbN7FpafFj/xRFUZzASUOwDEif6N4LXGFjn2XAJEMgIu8n4TGwcmV55p2WkvfM8nAQRVGUXDiZI8gUEDcF7IMx5svGmG3GmG2dnapnriiKUkqcNAS9wIq018uB0wXsoyiKojiIk4bgGWCdiHSLSDVwF/DAlH0eAO5OVg9dCZybS/kBRVGUSsCxHIExJioiHwIeAtzAvxpjXhaRe5Lv3wc8CLwB6AGCwB84tR5FURQlM45KTBhjHiRxsU/fdl/acwN80Mk1KIqiKLnRzmJFUZR5jhoCRVGUeY4aAkVRlHmOJML0lYOIDADHZ9wxMx2Ad8a9lHKj5+X8Q8/J+Ukx52WVMSZjI1bFGYJiEJHdxphts70OZTJ6Xs4/9Jycnzh1XjQ0pCiKMs9RQ6AoijLPmW+G4MuzvQAlI3pezj/0nJyfOHJe5lWOQFEURZnOfPMIFEVRlCmoIVAURZnnzGlDkJx49gMR2S8i+0Rku4i0icgjInIo+W/rbK9zPiEiHxWRl0XkJRH5tojU6jkpPyLyryLSLyIvpW3Leh5E5BPJ2eIHROT1s7PquU2Wc/IPyevXXhH5kYi0pL1XsnMypw0B8P+AXxhjNgCbgX3Ax4FfGWPWAb9KvlbKgIgsAz4CbDPGXERClfYu9JzMBl8Dbp6yLeN5SM4avwu4MPmZe5MzyZXS8jWmn5NHgIuMMRcDB4FPQOnPyZw1BCLSDLwa+AqAMSZsjBkmMSf568ndvg78zmysbx7jAepExAPUkxhEpOekzBhjfgP4pmzOdh7uAL5jjAkZY46SkI2/vBzrnE9kOifGmIeNMdHky10khndBic/JnDUEwGpgAPiqiDwvIveLSAOwyBp+k/x34Wwucj5hjDkFfBY4QWIu9TljzMPoOTlfyHYess0WV8rLe4CfJ5+X9JzMZUPgAbYCXzTGbAECaMhhVknGnO8AuoGlQIOIvGN2V6XYwNZsccU5ROQvgCjwLWtTht0KPidz2RD0Ar3GmN8mX/+AhGHoE5ElAMl/+2dpffORG4GjxpgBY0wE+CFwFXpOzheynQedLT6LiMi7gDcCbzcTjV8lPSdz1hAYY84CJ0VkfXLTDcArJOYkvyu57V3AT2ZhefOVE8CVIlIvIkLinOxDz8n5Qrbz8ABwl4jUiEg3sA54ehbWN+8QkZuBPwduN8YE094q6TmZ053FInIJcD9QDRwhMRPZBXwPWEniwvR7xpipSTPFIUTkU8CdJNzc54H3AY3oOSkrIvJt4DoSssZ9wF8CPybLeUiGJt5D4rz9iTHm59OPqhRDlnPyCaAGGEzutssYc09y/5KdkzltCBRFUZSZmbOhIUVRFMUeaggURVHmOWoIFEVR5jlqCBRFUeY5aggURVHmOWoIlIpBRGIi8kLao2u211QqRGSLiNyffP5uETEickPa+29KbntL8vXjIrIt7f2udNXKDMfvFJFfOPkzKJWLGgKlkhgzxlyS9jhmvSEJKvn3+ZPAP6W9fhF4a9rru4A9hR7cGDMAnBGRqws9hjJ3qeQ/HGWek7wL3ici9wLPAStE5M9E5Jmkfvun0vb9i6Ru+y+TcxA+ltyeurMWkQ4ROZZ87k5qwVvH+qPk9uuSn7HmXHwr2SWNiFwmIjtEZI+IPC0iTSLyRLKx0VrHUyJy8ZSfowm42BiTfqF/ArhcRKpEpBFYC7xg8//l/jSvaUBE/jL51o+Bt9v871XmEZ7ZXoCi5EGdiLyQfH4U+CiwHvgDY8wHROR1JFrtLychyvWAiLyahODgXcAWEr/zzwHPzvBd7yWhjnqZiNQAT4nIw8n3tpDQgT8NPAVcLSJPA98F7jTGPJOUQR8j0dn+buBPROQCoMYYs3fKd20DpoZ1DPBL4PXAAhKSAt1T9vmWiIwln1cDcQBjzPsARGQV8BAJnXuA3cCnZ/i5lXmIGgKlkhgzxlxivUjmCI4bY3YlN70u+Xg++bqRhGFoAn5kabWIyAM2vut1wMVWTJ7ExXgdEAaeNsb0Jo/1AtAFnAPOGGOeATDGjCTf/z7wP0Tkz0jIAXwtw3ctISGZPpXvkBjkswD4ryTCR+m83RizO/k9XcDPrDdEpBb4PvAhY8zx5OZ+EqqvijIJNQRKpRNIey7AZ4wxX0rfQUT+hOwSvVEmQqS1U471YWPMQ1OOdR0QStsUI/F3JJm+wxgTFJFHSMhv/z6Ju/+pjE35buuzT4vIRSQM4MFkBMou9wE/NMb8Mm1bbfK7FGUSmiNQ5hIPAe9JxtQRkWUishD4DfAmEalLxuNvS/vMMeDS5PO3TDnWH4tIVfJYF0hisFE29gNLReSy5P5NkpjCBonw0OeBZ7KI6e0jkQPIxCeY7gnkREQ+CDQZY/73lLcuYHoISlHUI1DmDsaYh0VkI7AzeffsB95hjHlORL5LItl6nEQi1uKzwPdE5J3Ao2nb7ycR8nkumQweIMcITWNMWETuBP5JROpI3HnfCPiNMc+KyAjw1Syf3S8iC0SkyRgzOuW9QhQlPwZE0vIp9xlj7gOuB/6zgOMpcxxVH1XmHSLyVyQu0J8t0/ctBR4HNhhj4ln2+Sgwaoy538F1/Aa4wxgz5NR3KJWJhoYUxUFE5G7gt8BfZDMCSb7I5NxDqdfRCfwfNQJKJtQjUBRFmeeoR6AoijLPUUOgKIoyz1FDoCiKMs9RQ6AoijLPUUOgKIoyz/n/i7ZXoCMvuEMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plt.plot(y, label = 'Sky Model')\n",
    "#plt.plot(obs,lw = 3, ls = \":\", label = 'Observed Data')\n",
    "plt.plot(abs(np.abs(y) - np.abs(obs)))\n",
    "labels=([\"60\", \"80\", \"100\", \"120\"])\n",
    "plt.xticks([0, 10, 20, 30], labels)\n",
    "plt.ylabel('Temperature (K)')\n",
    "plt.xlabel('Frequency (MHz)')\n",
    "#plt.legend(loc = 'upper right', prop ={'size': 16})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "expensive-bathroom",
   "metadata": {},
   "source": [
    "# Fit versus f at each t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "sharp-deposit",
   "metadata": {},
   "outputs": [],
   "source": [
    "polydim = 8\n",
    "a = np.zeros((polydim,alldataavg.shape[0]))\n",
    "\n",
    "for j in range(0, alldataavg.shape[0]):\n",
    "    Ninv = np.identity(alldataavg.shape[1])*(1/alldatavar[j,:])\n",
    "    P = alldataavg[j,:]\n",
    "    columns = np.zeros((len(P), polydim))\n",
    "    columns[:,0] = P\n",
    "    for k in range(1, polydim):\n",
    "        temparray = np.zeros(len(P))\n",
    "        for i in range(0, alldataavg.shape[1]):\n",
    "            temparray[i]= P[i]*(i*2)**k\n",
    "        columns[:,k] = temparray  \n",
    "    data = np.stack(columns, axis=0)\n",
    "    data[np.isnan(data)] = 0   \n",
    "        \n",
    "    d[np.isnan(d)] = 0\n",
    "    AT = np.transpose(data)\n",
    "    first = np.matmul(np.matmul(AT, Ninv), data)\n",
    "    firstinv = np.linalg.inv(first)\n",
    "    second = np.matmul(np.matmul(AT, Ninv), d[j,:])\n",
    "    final = np.matmul(firstinv, second)  \n",
    "    a[:,j] = final "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "nervous-concert",
   "metadata": {},
   "outputs": [],
   "source": [
    "ktmp = np.zeros((alldataavg.shape[0],alldataavg.shape[1], polydim))\n",
    "for i in range(0,alldataavg.shape[1]):\n",
    "    for j in range(0,alldataavg.shape[0]):\n",
    "        for l in range(0, polydim):\n",
    "            ktmp[j,i,l] = a[l,j]*((i*2)**l)\n",
    "k = np.nansum(ktmp, axis=2)\n",
    "kaverage = np.average(k, axis=0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "plain-series",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([<matplotlib.axis.XTick at 0x129334850>,\n",
       "  <matplotlib.axis.XTick at 0x129335510>,\n",
       "  <matplotlib.axis.XTick at 0x129335090>,\n",
       "  <matplotlib.axis.XTick at 0x129210210>],\n",
       " [Text(0, 0, '60'), Text(10, 0, '80'), Text(20, 0, '100'), Text(30, 0, '120')])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAERCAYAAABrWly6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAl6ElEQVR4nO3deXRdZ3nv8e+j2ZpH27I1WZ6HxJbjOLZDQyCQkBAItGkJZUoKmASa20Bpb7kdKLQsem+7bgOhxIQw5UIphYaQkjkxmUd5SuLZlgcpHjRZgyXbmp77x9lyZEWKZflIZ/p91jrr7LP3Pvs88rZ/fvWed7/b3B0REYl9SZEuQEREwkOBLiISJxToIiJxQoEuIhInFOgiInFCgS4iEiciGuhm9kMzazSz18N0vAoze9TMtpvZNjOrCsdxRURiQaRb6D8G3hfG490D/LO7LwRWAo1hPLaISFSLaKC7+9NA69B1ZjbbzB42sw1m9oyZLRjLscxsEZDi7o8Fxz7u7t3hr1pEJDpFuoU+kruAW939IuDLwHfH+L55QJuZ3Wtmm8zsn80secKqFBGJMimRLmAoM8sG1gC/NLPB1enBtt8Hvj7C295w96sI/Sy/B9QAB4FfADcCP5jYqkVEokNUBTqh3xja3H3Z8A3ufi9w79u8twHY5O51AGZ2H7AKBbqIJIio6nJx9w5gn5n9IYCFLB3j218BCsysJHj9bmDbBJQpIhKVIj1s8efAC8B8M2sws08DHwM+bWZbgK3AdWM5lrv3E+pzf8LMXgMM+P7EVC4iEn1M0+eKiMSHqOpyERGR8YvYl6LFxcVeVVUVqY8XEYlJGzZsaHb3kpG2RSzQq6qqqK2tjdTHi4jEJDM7MNo2dbmIiMQJBbqISJxQoIuIxAkFuohInFCgi4jECQW6iEicUKCLiMSJmAv0nUc6+aeHdtBxsjfSpYiIRJWYC/SDrd2se2ovexuPR7oUEZGoEnOBPqs4C4B9zV0RrkREJLrEXKBXFGaSZAp0EZHhYi7Q01KSKC/MpE6BLiJyhpgLdAh1u+xrUqCLiAwVu4He3IVuziEi8qaYDPTq4ixO9PZztONUpEsREYkaMRnos4qzAahr1tBFEZFBsRnoJRq6KCIyXEwGemluBhmpSfpiVERkiJgM9KQko6ooSy10EZEhYjLQAapLFOgiIkPFbKDPKs7iYGs3vf0DkS5FRCQqnDXQzWy+mW0e8ugws9uG7WNm9m0z22Nmr5rZ8gmrODCrOJu+Aafh2ImJ/igRkZiQcrYd3H0nsAzAzJKBN4BfD9vtamBu8LgEuDN4njBvTtJ1/PSyiEgiO9culyuAve5+YNj664B7PORFIN/MSsNS4SiqgxCv00gXERHg3AP9BuDnI6yfCdQPed0QrDuDma01s1ozq21qajrHjz5TQVYa+Zmp+mJURCQw5kA3szTgg8AvR9o8wrq3TLTi7ne5+wp3X1FSUjL2KkcxOKeLiIicWwv9amCjux8dYVsDUD7kdRlw6HwKGwsFuojIm84l0D/KyN0tAPcDnwxGu6wC2t398HlXdxbVxVkcbj9Jd0/fRH+UiEjUG1Ogm1km8F7g3iHrbjazm4OXDwJ1wB7g+8Dnw1zniAYn6drf3D0ZHyciEtXOOmwRwN27gaJh69YNWXbgC+Et7eyG3l900Yzcyf54EZGoErNXigJUFWcCobHoIiKJLqYDPTMthdK8DN1fVESEGA900EgXEZFBCnQRkTgRF4He1t1La1dPpEsREYmomA/06pI3J+kSEUlkMR/op28YrUm6RCTBxXyglxVMISXJ1I8uIgkv5gM9NTmJisJMBbqIJLyYD3TQSBcREYizQB8YeMuMvSIiCSM+Ar0ki1N9AxzuOBnpUkREIiY+An1wki6NdBGRBBYXgV4dDF3UWHQRSWRxEejTctOZkpqsSbpEJKHFRaCbmUa6iEjCi4tAh9AUAAp0EUlk8RPoxVnUt3bT0zcQ6VJERCIibgJ9VkkWAw4HW3V/URFJTPET6KdHuqjbRUQSU/wEepGm0RWRxDamQDezfDP7lZntMLPtZrZ62PbLzazdzDYHj7+bmHJHl5eZSlFWmlroIpKwUsa437eAh939ejNLAzJH2OcZd782fKWdu1nFWZoXXUQS1llb6GaWC1wG/ADA3XvcvW2C6xoXjUUXkUQ2li6XaqAJ+JGZbTKzu80sa4T9VpvZFjN7yMwWj3QgM1trZrVmVtvU1HQ+dY9oVkkWjZ2nOH6qL+zHFhGJdmMJ9BRgOXCnu9cAXcBfDdtnI1Dp7kuBO4D7RjqQu9/l7ivcfUVJScn4qx5FdTBJ13610kUkAY0l0BuABnd/KXj9K0IBf5q7d7j78WD5QSDVzIrDWukYnL6/qAJdRBLQWQPd3Y8A9WY2P1h1BbBt6D5mNt3MLFheGRy3Jcy1nlVlUSZmmkZXRBLTWEe53Ar8LBjhUgfcZGY3A7j7OuB64BYz6wNOADe4+6TfPigjNZkZeVM0Fl1EEtKYAt3dNwMrhq1eN2T7d4DvhK+s8dMkXSKSqOLmStFBs4qzqGvuIgK/IIiIRFRcBnrnyT5aunoiXYqIyKSKy0AHTdIlIokn7gL99P1FNdJFRBJM3AX6zIIppCabxqKLSMKJu0BPTjIqi7I0dFFEEk7cBTpo1kURSUxxGejVxVkcaOmmf0BDF0UkccRloM8qzqKnf4BDbSciXYqIyKSJ20AHTdIlIoklPgO9JBiL3qQvRkUkccRloJdkp5OdnqKLi0QkocRloJvZ6TldREQSRVwGOuj+oiKSeOI60N9oO8HJ3v5IlyIiMiniNtCrS7Jwh4Ot3ZEuRURkUsRtoM8uCU3StfuoRrqISGKI20CfNy2HtJQktjS0RboUEZFJEbeBnpaSxJIZuWw6eCzSpYiITIq4DXSAmooCXm1op6dvINKliIhMuDgP9HxO9Q2w40hHpEsREZlwYwp0M8s3s1+Z2Q4z225mq4dtNzP7tpntMbNXzWz5xJR7bpZXFACw6WBbZAsREZkEY22hfwt42N0XAEuB7cO2Xw3MDR5rgTvDVuF5KM3LYFpuuvrRRSQhnDXQzSwXuAz4AYC797h727DdrgPu8ZAXgXwzKw13sefKzKgpL2BTfVukSxERmXBjaaFXA03Aj8xsk5ndbWZZw/aZCdQPed0QrDuDma01s1ozq21qahp30eeipiKfAy3dtBw/NSmfJyISKWMJ9BRgOXCnu9cAXcBfDdvHRnjfW24X5O53ufsKd19RUlJyzsWOR4360UUkQYwl0BuABnd/KXj9K0IBP3yf8iGvy4BD51/e+btgZh4pScamevWji0h8O2ugu/sRoN7M5gerrgC2DdvtfuCTwWiXVUC7ux8Ob6njMyUtmYWluWqhi0jcSxnjfrcCPzOzNKAOuMnMbgZw93XAg8A1wB6gG7hpAmodt5qKfP5rQwP9A05y0ki9QyIisW9Mge7um4EVw1avG7LdgS+Er6zwqqnI554XDrC7sZMF03MjXY6IyISI6ytFB9WUh74Y3XigLbKFiIhMoIQI9MqiTAqz0nSBkYjEtYQI9NAFRvm6wEhE4lpCBDqE+tH3NB6n/URvpEsREZkQCRTooX70LWqli0icSphAv7AsDzPYqH50EYlTCRPoORmpzJuaowuMRCRuJUygAyyvzGdzfRsDA2+ZZkZEJOYlVKDXlBfQfqKXfS1dkS5FRCTsEivQK/IBzbwoIvEpoQJ9dkk2Oekp+mJUROJSQgV6UpKxrCJfLXQRiUsJFegQGo++80gHXaf6Il2KiEhYJWCg5zPg8GpDe6RLEREJq4QL9GVl+QC6g5GIxJ2EC/SCrDSqi7PUjy4icSfhAh0Ivhg9Rui+HCIi8SEhA315RQHNx3toOHYi0qWIiIRNQgb64AVGGo8uIvEkIQN9/rQcpqQmqx9dROJKQgZ6SnISF5bl6Q5GIhJXEjLQIXSB0bZD7Zzs7Y90KSIiYTGmQDez/Wb2mpltNrPaEbZfbmbtwfbNZvZ34S81vGoq8untd7Ye0gVGIhIfUs5h33e5e/PbbH/G3a8934Imy9CZFy+qLIxsMSIiYZCwXS5TczIoK5iiL0ZFJG6MNdAdeNTMNpjZ2lH2WW1mW8zsITNbPNIOZrbWzGrNrLapqWlcBYdTTUUBmzR0UUTixFgD/VJ3Xw5cDXzBzC4btn0jUOnuS4E7gPtGOoi73+XuK9x9RUlJyXhrDpua8nwOtZ/kSPvJSJciInLexhTo7n4oeG4Efg2sHLa9w92PB8sPAqlmVhzmWsPuzX50tdJFJPadNdDNLMvMcgaXgSuB14ftM93MLFheGRy3JfzlhtfiGXmkpSRpPLqIxIWxjHKZBvw6yOsU4N/d/WEzuxnA3dcB1wO3mFkfcAK4wWNg5qu0lCSWzMhVC11E4sJZA93d64ClI6xfN2T5O8B3wlva5KipKOCnLx6gt3+A1OSEHfQjInEg4ROspiKfU30DbD/cEelSRETOiwK9ogCA2v3qdhGR2JbwgT4zfwpzpmbz6LYjkS5FROS8JHygA1yzZDov72ulqfNUpEsRERk3BTpwzYWlDDg8slWtdBGJXQp0Qje8qC7O4qHXD0e6FBGRcVOgA2bGNReU8sLeFlqOq9tFRGKTAj1w9QXTGXB4dNvRSJciIjIuCvTAotJcKosyefA1dbuISGxSoAcGu12e39vCsa6eSJcjInLOFOhDXLOklP4B5zF1u4hIDFKgD7FkZi7lhVN4QN0uIhKDFOhDmBnXLCnluT3NtHf3RrocEZFzokAf5uoLSukbcE0FICIxR4E+zNKyPGbmT+Gh1xXoIhJbFOjDmBlXL5nOM7ubaD+hbhcRiR0K9BFcc2Epvf3OE9s12kVEYocCfQTLyvIpzcvQRUYiElMU6CNISjKuXlLK07ua6TypbhcRiQ0K9FFcc8F0evoHWL+jMdKliIiMiQJ9FMsrCpiWm84Dr6rbRURiw5gC3cz2m9lrZrbZzGpH2G5m9m0z22Nmr5rZ8vCXOrkGu12e3NXE8VN9kS5HROSszqWF/i53X+buK0bYdjUwN3isBe4MR3GRdvWS6fT0qdtFRGJDuLpcrgPu8ZAXgXwzKw3TsSNmRVUhJTnpPKTRLiISA8Ya6A48amYbzGztCNtnAvVDXjcE685gZmvNrNbMapuams692kmWnGS8b/F0frezke4edbuISHQba6Bf6u7LCXWtfMHMLhu23UZ4j79lhftd7r7C3VeUlJScY6mRcc0FpZzsHeB3O6L/PyARSWxjCnR3PxQ8NwK/BlYO26UBKB/yugw4FI4CI23lrEKKstJ4UDeQFpEod9ZAN7MsM8sZXAauBF4fttv9wCeD0S6rgHZ3j4sETE4yrloynfXbGznR0x/pckRERjWWFvo04Fkz2wK8DDzg7g+b2c1mdnOwz4NAHbAH+D7w+QmpNkLef0EpJ3r7eWqXRruISPRKOdsO7l4HLB1h/bohyw58IbylRY9LZhVSmJXGA68d4X1LYn7wjojEKV0pOgYpyUlctXga67cf5WSvul1EJDop0Mfo6iWldPX089QujXYRkeikQB+j1bOLKMhM5Ze19WffWUQkAhToY5SanMSNa2bx+PZGXmtoj3Q5IiJvoUA/Bze9o4rcjBS+9cSuSJciIvIWCvRzkJuRymd/r5rHtzfyakNbpMsRETmDAv0c3XhpFfmZqdz++O5IlyIicgYF+jnKCVrp63c0srm+LdLliIicpkAfh0+tqaIgM5XbH1dfuohEDwX6OGSnp/DZy6p5cmcTGw8ei3Q5IiKAAn3cPrU61Er/lvrSRSRKKNDHKSs9hbWXzeapXU1sOKBWuohE3lkn55LRfXJ1Jd9/po7bH9/F//v0JZEuR0QioPNkL9sPd7K7sZP0lGSKstIoDB5F2Wlkpk1ezCrQz0NWegqfu6yabz60gw0HWrmosjDSJYnIBHF3jnacYtvhdrYd6mDroQ62He7gQEv3274vIzWJoqz0N0M+K40rF0+bkJlbFejn6RNBK/1fH9vNTz+jVrpIvOjtH+CVfa08vbuZrYdCId7S1XN6e1VRJotn5PKHF5WxeEYe86bn0Nc/QEtXD63He2jt6gktd50KnkOPPY3HmT01e0JqVqCfp8y0FD532Wy+8eB2XtnfysVVaqWLxKpjXT08uauRx7c38vTOJjpP9ZGWnMT86Tm8Z+E0Fs3IZdGMXBaW5pKdPnJ8VhZlTXLVb1Kgh8HHV1Xyvafr+NfHdvHvn10V6XJEZIzcnT2Nx3liRyNPbD/KhgPHGHAoyUnn/ReWcsXCaVw6p2hS+8HPR2xUGeWmpCVz8zur+ccHtvNSXQuXVBdFuiQRGUVf/wAv72/l0a1HWb+jkYOtoT7wxTNy+dN3z+WKBVO5YGYeSUkW4UrPnQI9TAZb6bc/vpufr1Wgi0STk739PLu7mUe2HuHx7Uc51t1LekoSl84p5nPvrObdC6ZSmjcl0mWeNwV6mGSkJnPLO2fz9d9u48W6FlaplS4SUZ0ne/ndziYe2XqEJ3c00tXTT056ClcsnMpVi6fzzvklMdOVMlbx9dNE2B9fUsG6p/byr4/t4hefWx3pckQSTlPnKdbvOMrDrx/huT0t9PQPUJydxgeXzeR9S6azurqItJT4vZ5yzIFuZslALfCGu187bNvlwG+AfcGqe93962GqMWZkpCbz+ctn8/f/vY3n9zazZnZxpEsSiWsDA87WQx2s39HI+h1H2RLcTaysYAqfWF3J+5ZMZ3lFAckx2B8+HufSQv8zYDuQO8r2Z4YHfSK6YWUF656q4xsPbOe/bllDRmpypEsSiSvHT/Xx7O5m1u84yu92NtHUeQozWFqWz5feO493L5jK4hm5mCVGiA81pkA3szLg/cA3gC9NaEUxLiM1mX/40BI+e08tX/vvrXzz9y+MdEkiMc3dqWvu4qmdTazf0chL+1ro7Xdy0lO4bF4J71owlcvnl1CcnR7pUiNurC3024G/BHLeZp/VZrYFOAR82d23Dt/BzNYCawEqKirOrdIY8t5F0/j85bP57pN7qSkv4I8uLo90SSIxpb61m+f3NvPC3hae39tCY+cpAGaXZHHjmiretWAqF1cVkpocv/3h43HWQDeza4FGd98Q9JWPZCNQ6e7Hzewa4D5g7vCd3P0u4C6AFStW+Dhrjgl/fuV8Xm1o529+8zqLZuSyZGZepEsSiVpH2k/yQl0zz+9p4YW6FhqOnQCgODuNVdVFrJldzKVziiJ6FWYsMPe3z1Uz+ybwCaAPyCDUh36vu3/8bd6zH1jh7s2j7bNixQqvra0dT80xo+X4KT5wx7MkJRm/vfUd5GemRbokkYgb7ELZeOAYGw+28VJdC3XNXQDkTUllVXUha2YXs3p2EXOnZidkX/jbMbMN7r5ixG1nC/RhB7qcUHfK8FEu04Gj7u5mthL4FaEW+6gHT4RAB9hc38YfrXuBNXOK+OGnLo7Jq89EzkfXqT621Lex8eAxNhw4xqb6Ntq6ewHIyUhhRWXB6QBfWJqbMCNSxuvtAn3c49DN7GYAd18HXA/cYmZ9wAnghrcL80SyrDyfv/vAIv7mvtf59vrd3PaeeZEuSWTCDAw4dc3H2VLfzqb6Y2w80MaOIx0MBGkwZ2o2Vy6axvKKApZXFjCnJFuNnDA6pxZ6OCVKCx1Cv2L++S+38OtNb/DDGy/mXfOnRrokkfPm7tS3nuDVN9p4taGdLfVtbD3UwfFTfQBkpSVTU1HA8op8aioLqCnPV7djGExIC13Gzsz4xocuYPvhTm77j8389tZ3UF6YGemyRMbM3TnScZLXGtpD4d3QxmtvtJ/uOklLTmJhaQ4frpnJBWV5LC3LZ87UbHWfTDK10CfRgZYurr3jWSoKM3XRkUS1xs43w/u1N0LPzcdDQweTk4y5U7NZWpZ/OrznT8+J60vqo4la6FGisiiL2z+yjE//pJa/ve91/s/1F+obfIm41q6eUIs7CPDX32jnSMdJAMxgTkk2l80r5sKZeVxQlsei0jympKkxEo0U6JPsioXTuPXdc7hj/R6WVxbw0ZXxe4GVRJ+Tvf1sPdTO5vpQn/fm+rbT84GbQXVxFquqC7mgLJ8Ly/JYVJpL1ih35pHoozMVAbe9Zx6b69v46m+2sqg0l6Xl+ZEuSeLQwICzt+k4m4Pg3tLQxo7DnfQFQ05m5GWwtDyfP76kgqVl+SyZmUtORmqEq5bzoT70CGnt6uEDdzxLT/8A3//kCpYp1OU8dff0sbm+jQ37j7Hh4DE2HjhGx8nQiJOc9BQuLA/1dy8rDz2m5mZEuGIZj7BdWBROiR7oALuPdnLTj1+hsfMU//sPLuDDNWWRLkliyJH2k9QeaKV2f+iCnW2HO+gPWt/zpmVzUWUByysKqKnIp7pY473jhb4UjVJzp+Vw/5++g1t+uoEv/mILO4508pdXLdBQL3kLd2d/Szcv1rXwUl0Lr+w/xhttoflOMlKTWFaezy3vnH06xPMy1XWSiBToEVaYlcZPP3MJX/vvrXzvqTp2HenkWx+tIVd9mQnNPdT//WJdKy/ta+WlujdnHCzOTueSWYV8+h2zWFFVwMLSXM06KIACPSqkJifxjx+6gAXTc/n7+7fy4X97jrs/dTGzijWzXKJwd/Y0HueFuhZeqmvlpX0tNB/vAWBabjqrqotYVV3EJdWFVBdnabirjEiBHkU+vqqSOVOzueWnG7juO8/ybx9bzu/NLYl0WTJBBuf8fm5PaM7vwQt3ZuRlcNncEi6pLuSSWUVUFmUqwGVM9KVoFKpv7eaz99Sy62gnf/P+Rdx0aZX+QceBxs6ToRs27Gnh+bpm6ltDfeAlOemsmV3EpcGMg2UFU3S+ZVT6UjTGlAdTA3zpPzfz9d9uY8eRDv7hQ0tIT9HVebGk/UQvL9WFWt/P7Wlmd+NxAHIzUlhVXcRn3lHNmtlFzNGc3xImCvQolZWewp0fu4jbn9jNt5/Yzc6jx/nqBxaxvKIg0qXJKE729rPhwDGe29PMc3tbeK2hjQEPjUK5uKqQP7iojDWzi1g8I08jmWRCqMslBjz42mH+9r7Xaenq4d0LpvKl987TLe2iQP+A89ob7Ty3p5nn9zbzyv5j9PQNkJxkLCvP59LZRayZU0xNRb5+u5Kw0YVFcaDrVB8/fn4/33tqLx0n+7h6yXS++N55zJv2dvftlnDqH3C2H+7gxbqW0Hjwfa10BldiLpiew5rZxbxjbhErZxWRrflPZIIo0ONI+4lefvDsPn747D66evq4bukMbnvPPKo0xDHsBgacHUc6eWEwwOtaTl9KPyuYxGr17GLWzC6iODs9wtVKolCgx6HWrh6+9/RefvL8fnr7neuXl3HrFXMoK9CNM8arr3+A7Yc7qT3QeroFPngDh8qiTFYHY8FXVRcxPU/zoEhkKNDjWGPnSb77u738+0sHcZyPXFzO7y8vY1lZvubuOIv27t7TNy7ecOAYm+vbONHbD0B54ZQzAnxG/pQIVysSokBPAIfaTvCd3+3hl7X19PY7JTnpvHfRNK5aPJ3V1UUJfzeZgQFnX0sXGw6EZiHccODY6WGEyUnG4hm5LK8o4KLK0EMBLtFKgZ5A2k/08uTORh7ZeoQndzbR3dNPTnoKly+YylWLp/HOeSVxP+f1qb5+dh89zrZDHWw91M62wx1sP9x5+ubFuRkpp4P7ospClpbnkZmmLzElNoQl0M0sGagF3nD3a4dtM+BbwDVAN3Cju298u+Mp0Cfeyd5+ntvTzKNbj/L49qO0dPWQlpzEmjlFXLloOkvL85gzNTtmh9S5Oy1dPextPM7WQx1sO9zB1kMd7GnspLc/9Pc6Ky2ZhaW5LJqRe7oVPrtEU8lK7ArXlaJ/BmwHckfYdjUwN3hcAtwZPEsEZaQmc8XCaVyxcBr9A86GA8d4ZOuRoPX+GhDqbphVnMX86TnMn5bD/Ok5LJieQ3lBZlSEnrvT2HmK/c1dHGjpZn/Lmc+DrW4IzUK4eEYul88vYfGMXBaV5lJVlBUVP4fIZBhToJtZGfB+4BvAl0bY5TrgHg819180s3wzK3X3w+ErVc5HcpKxclYhK2cV8jfvX8jepi62He5g15FOdhzp5NWGNh549c3TNSU1mXnTspk/PYfpuRnkZaZRkJlKfmYqeVMGl9PIzUgh5Rynbh0YcDpO9tLS1cOxrp4znluHLB/tOMmBlu7TX1QCpCQZ5YWZVBZlcnFVIZVFmVQVZ7G4NFd34JGEN9YW+u3AXwKjXcUyE6gf8rohWKdAj0Jmxpyp2cyZmg1L31zfdaqPXUc72Xmkk53B8/odTbR0neLteuZyM1LIz0wjJdkYGHAGPHQRzoD7Gc/9wbYTvf2n76wzXGZaMgWZaRRlpzEzfwqXzimmqiiTyqIsqoqymJGfcc7/gYgkirMGupldCzS6+wYzu3y03UZY95Z/sWa2FlgLUFGhu91Hm6z0FGoqCqgZNl/MwIDTebKPY909tJ3opa27h7bu0POx7l7aT/RyrLuHvgEn2YzkJMOM08tJSfbmshlT0pIozEqnKCuNgqy0M54zUmOzP18kGoylhX4p8EEzuwbIAHLN7Kfu/vEh+zQA5UNelwGHhh/I3e8C7oLQl6LjrlomVVKSkZeZqtuaiUS5s/7u6u5fcfcyd68CbgDWDwtzgPuBT1rIKqBd/eciIpNr3INvzexmAHdfBzxIaMjiHkLDFm8KS3UiIjJm5xTo7v4k8GSwvG7Iege+EM7CRETk3Gi4gIhInFCgi4jECQW6iEicUKCLiMQJBbqISJyI2PS5ZtYEHBjn24uB5jCWI+Gh8xJ9dE6i0/mcl0p3LxlpQ8QC/XyYWe1o00dK5Oi8RB+dk+g0UedFXS4iInFCgS4iEidiNdDvinQBMiKdl+ijcxKdJuS8xGQfuoiIvFWsttBFRGQYBbqISJyIiUAP7lH6KzPbYWbbzWy1mRWa2WNmtjt4Ljj7kSRczOyLZrbVzF43s5+bWYbOyeQzsx+aWaOZvT5k3ajnwcy+YmZ7zGynmV0Vmarj2yjn5J+D/HrVzH5tZvlDtoXtnMREoAPfAh529wWE7oK5Hfgr4Al3nws8EbyWSWBmM4H/Aaxw9yVAMqGbn+icTL4fA+8btm7E82Bmiwidp8XBe75rZrrnX/j9mLeek8eAJe5+IbAL+AqE/5xEfaCbWS5wGfADAHfvcfc24DrgJ8FuPwE+FIn6ElgKMMXMUoBMQrcc1DmZZO7+NNA6bPVo5+E64D/c/ZS77yN0Q5qVk1FnIhnpnLj7o+7eF7x8kdBtOiHM5yTqAx2oBpqAH5nZJjO728yygGmDt7kLnqdGsshE4u5vAP8CHAQOE7rl4KPonESL0c7DTKB+yH4NwTqZXH8CPBQsh/WcxEKgpwDLgTvdvQboQr/KR1TQJ3sdMAuYAWSZ2fD7zEr0sRHWadzyJDKzvwb6gJ8Nrhpht3Gfk1gI9Aagwd1fCl7/ilDAHzWzUoDguTFC9SWi9wD73L3J3XuBe4E16JxEi9HOQwNQPmS/MkJdZTIJzOxTwLXAx/zNC4DCek6iPtDd/QhQb2bzg1VXANuA+4FPBes+BfwmAuUlqoPAKjPLNDMjdE62o3MSLUY7D/cDN5hZupnNAuYCL0egvoRjZu8D/ifwQXfvHrIprOckJq4UNbNlwN1AGlAH3EToP6P/BCoIBcwfuvvwL4dkgpjZ14CPEPr1cRPwGSAbnZNJZWY/By4nNB3rUeCrwH2Mch6CX/n/hNB5u83dH3rrUeV8jHJOvgKkAy3Bbi+6+83B/mE7JzER6CIicnZR3+UiIiJjo0AXEYkTCnQRkTihQBcRiRMKdBGROKFAl0lnZv1mtnnIoyrSNYWLmdWY2d3B8o1m5mZ2xZDtHw7WXR+8ftLMVgzZXjV0lr4Rjl9iZg9P5M8gsUuBLpFwwt2XDXnsH9xgIbH89/J/AXcMef0a8NEhr28Atoz34O7eBBw2s0vHewyJX7H8D0fiRNAq3W5m3wU2AuVm9hdm9kowf/TXhuz718G80Y8H87B/OVh/uqVrZsVmtj9YTg7moh481ueC9ZcH7xmcZ/9nwVWvmNnFZva8mW0xs5fNLMfMngkucBus4zkzu3DYz5EDXOjuQwP7GWClmaWaWTYwB9g8xj+Xu4f8FtNkZl8NNt0HfGyMf7ySQFIiXYAkpClmtjlY3gd8EZgP3OTunzezKwldAr2S0ORF95vZZYQmZrsBqCH0d3cjsOEsn/VpQrNBXmxm6cBzZvZosK2G0DzUh4DngEvN7GXgF8BH3P2VYPrmE4SuVL4RuM3M5gHp7v7qsM9aAQzvLnHgceAqII/Qpd6zhu3zMzM7ESynAQMA7v4ZADOrBB4hNM82QC3wj2f5uSUBKdAlEk64+7LBF0Ef+gF3fzFYdWXw2BS8ziYU8DnArwfnwjCz+8fwWVcCFw72WRMK1blAD/CyuzcEx9oMVAHtwGF3fwXA3TuC7b8E/tbM/oLQZdo/HuGzSglN9TzcfxC6IUge8OeEumWG+pi71wafUwX8dnCDmWUAvwT+1N0PBKsbCc1yKXIGBbpEi64hywZ8092/N3QHM7uN0acW7ePNLsSMYce61d0fGXasy4FTQ1b1E/r3YCN9hrt3m9ljhKYN/iNCrfHhTgz77MH3vmxmSwj9R7Yr6NkZq3XAve7++JB1GcFniZxBfegSjR4B/iToc8bMZprZVOBp4MNmNiXor/7AkPfsBy4Klq8fdqxbzCw1ONY8C90gZTQ7gBlmdnGwf46F7soEoW6XbwOvjDLp2HZCfeQj+QpvbZm/LTP7ApDj7v80bNM83tq1I6IWukQfd3/UzBYCLwSt2ePAx919o5n9gtCXigcIfeE46F+A/zSzTwDrh6y/m1BXysbgS88m3ubWeO7eY2YfAe4wsymEWsLvAY67+wYz6wB+NMp7d5hZnpnluHvnsG3jmUHvy0DvkO8b1rn7OuBdwAPjOJ7EOc22KDHLzP6eUND+yyR93gzgSWCBuw+Mss8XgU53v3sC63gauM7dj03UZ0hsUpeLyBiY2SeBl4C/Hi3MA3dyZt98uOsoAf6vwlxGoha6iEicUAtdRCROKNBFROKEAl1EJE4o0EVE4oQCXUQkTvx/lORHqEbcMsgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(kaverage)\n",
    "plt.xlabel('Frequency (MHz)')\n",
    "labels=([\"60\", \"80\", \"100\", \"120\"])\n",
    "plt.xticks([0, 10, 20, 30], labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "hawaiian-fancy",
   "metadata": {},
   "outputs": [],
   "source": [
    "obsall = k*alldataavg\n",
    "obs=np.nanmean(obsall, axis=0)\n",
    "v = []\n",
    "for j in range(flow, fhigh+2, 2):\n",
    "    v.append(j)\n",
    "v=np.array(v)\n",
    "fit = np.polyfit(np.log10(v/100), np.log10(np.abs(obs)), polydim)    \n",
    "power = np.zeros(int((fhigh-flow)/2)+1)\n",
    "i = 0\n",
    "while i <= polydim:\n",
    "    power = power + fit[polydim-i]*(np.log10(v/100))**i\n",
    "    i += 1\n",
    "y = 10**(power)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "opposite-yukon",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([<matplotlib.axis.XTick at 0x12923c450>,\n",
       "  <matplotlib.axis.XTick at 0x129242110>,\n",
       "  <matplotlib.axis.XTick at 0x12923cc50>,\n",
       "  <matplotlib.axis.XTick at 0x12962a250>],\n",
       " [Text(0, 0, '60'), Text(10, 0, '80'), Text(20, 0, '100'), Text(30, 0, '120')])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEGCAYAAAB8Ys7jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABET0lEQVR4nO3deXRb93Uv+u/GTAIgKBLgqIHipNGy5dC2ZMW2PMaJ7Thp08aJkzo3yXNuptvkdkhy89q09/bd9r1mNW0z1nXcZN24TtLMdprY8WzLsiRqtGbOFEmRIAGQAAli/r0/Dg4IkgAxnQMC5P6s5SURBIEjg9z8Yf/2b28SQoAxxlj50az2BTDGGMsPB3DGGCtTHMAZY6xMcQBnjLEyxQGcMcbKlK6YT2a320VLS0sxn5Ixxsre8ePHp4QQjqW3ZwzgRPQEgPsBOIUQu5Nu/yyAzwCIAPi1EOLPMz1WS0sLuru7c7pwxhhb74hoKNXt2aRQvgfg3iUPdjuABwHsEULsAvDVQi+QMcZYbjIGcCHEqwDcS27+JIC/E0IE4/dxqnBtjDHGVpDvJmYngFuI6AgRvUJEN6S7IxE9SkTdRNQ9OTmZ59MxxhhbKt8ArgOwAcA+AH8G4MdERKnuKIR4TAjRJYTocjiW5eAZY4zlKd8APgLgZ0JyFEAMgF25y2KMMZZJvgH8FwDuAAAi6gRgADCl0DUxxhjLQjZlhE8BOAjATkQjAL4C4AkATxDRWQAhAI8IbmvIGGNFlTGACyE+kOZTH1L4WhhbM35xchR37KhDlUm/2pfC1jA+Ss+YwoZdfnzuR6fw6zNXV/tS2BrHAZwxhY3NzAMAvPPhVb4SttZxAGdMYRPeAABgNhhZ5Sthax0HcMYU5vQGAQC+AAdwpi4O4IwpTF6Bz/EKnKmMAzhjChvnFAorEg7gjClMTqFwAGdq4wDOmMImfLwCZ8XBAZwxBQkhMD4TD+C8iclUxgGcMQV55yMIRmIAeAXO1McBnDEFyekTu8XIAZypjgM4YwqSSwjbHGbMBiPgHm9MTRzAGVOQnP9uq7NACMAfiq7yFbG1jAM4Ywpy+qQSwla7GQAf5mHq4gDOmIImvAFUmXRwWI0AAB8HcKYiDuCMKWjCG0CDzQSLUWq1z6WETE0cwBlT0Lg3iPqqhQDOKRSmJg7gjCnI6Q2gzmqCOR7AOYXC1JQxgBPRE0TkjM+/XPq5PyUiQUQ8kZ6te7GYgNMXRIPNCKuJUyhMfdmswL8H4N6lNxLRJgB3AxhW+JoYK0uuuRCiMbE4hRLiAM7UkzGACyFeBeBO8amvAfhzAHxSgTEsHOJZlELhFThTUV45cCJ6N4BRIcTpLO77KBF1E1H35ORkPk/HWFmQA3h9lRFGnQZ6LfFxeqaqnAM4EVUC+DKAv8zm/kKIx4QQXUKILofDkevTMVY2JuJ9wBtsJhARLEYdV6EwVeWzAm8DsBXAaSIaBLARwAkialDywhgrN+PeAIikRlYAYDbqeBOTqUqX6xcIId4CUCd/HA/iXUKIKQWvi7Gy4/QGUGs2Qq+V1kUWo47LCJmqsikjfArAYQDbiGiEiD6m/mUxVn6kU5jGxMdWE6dQmLoyrsCFEB/I8PkWxa6GsTI24Q2i0WZKfGw26uCeC63iFbG1jk9iMqaQCW8AdVULAdzCOXCmMg7gjCkgFInBNRdCfdXiFAqXETI1cQBnTAGTs/ESwqQVuNnAAZypiwM4YwpYOMSTlEIx6eAPRRGN8WFlpg4O4IwpYCI+Sq0uKYXC/VCY2jiAM6aAlCtwHurAVMYBnDEFTPiC0GsJNZWGxG0WuaUs58GZSjiAM6aAiRlpkINGQ4nbEitwDuBMJRzAGVPAhC+wKP8NcAqFqY8DOGMKmPAGF5UQApxCYerjAM6YAia8gUUbmACnUJj6OIAzViB/KAJfIMIpFFZ0HMAZK5A8yKHeungFbuYVOFMZB3DGCiTXgDfYFgdwvVYDk17DLWWZajiAM1ag5FmYS/FQB6YmDuCMFSgxjX7JJibALWWZujiAM1agCW8QFXotrMbl81EsPJWHqSibkWpPEJGTiM4m3fb3RHSRiM4Q0c+JqFrVq2SshEmj1KRJ9EuZDZxCYerJZgX+PQD3LrntdwB2CyH2ALgM4EsKXxdjZcPpDaLOujz/DcSHOnAKhakkYwAXQrwKwL3ktueEEPJ35ZsANqpwbYyVhfEUh3hkFqOO28ky1SiRA/8ogN+k+yQRPUpE3UTUPTk5qcDTMVY6hBCJFEoqZt7EZCoqKIAT0ZcBRAA8me4+QojHhBBdQoguh8NRyNMxVnK88xEEI7G0KRSLiXPgTD3Lt82zRESPALgfwJ1CCJ4Zxdal8RSDHJJZjTqEIjGEIjEYdFz0xZSV13cUEd0L4AsA3i2E8Ct7SYyVj1STeJLJx+m5lJCpIZsywqcAHAawjYhGiOhjAL4BwArgd0R0ioi+o/J1MlaSEsfoV9jEBLgfClNHxhSKEOIDKW7+rgrXwljZcfqkRlZLOxHKrNwTnKmIk3KMFWB8JgBbhR4mvTbl57kjIVMTB3DGCiANcki9+ga4JzhTFwdwxgow4Qum3cAEOIXC1MUBnLECOFc4hQlwCoWpiwM4y9urlycRCEdX+zJWTTQm4PQFOYXCVg0HcJaXvslZ/NETR/Hzk6OrfSmrxjUXRDQm0pYQAlI3QoBX4EwdHMBZXi6N+wAAlyd8q3wlq8fplUsI0wdwjYZgNmg5gDNVcABneemZmAUA9E3OrfKVrJ7xmZVPYcos3FKWqYQDOMtLj1Naefc5Z1f5SlbPhC/9LMxkFqMOs9xSlqmAAzjLS288cI9Oz8O/ToPThDcIIsBhySKA8wqcqYADOMtZJBpD/+QcNtVUAAD612kaxekNwG4xQqdd+cfIYtJxDpypggM4y9kVzzxC0Rju3dUAQKpIWY/GM5zClFmMPNiYqYMDOMtZT7zy5O6dDdDQ+t3InPAGVywhlJmNOvg4hcJUwAGc5awnnv/e2VSFzTWV63YF7vQGViwhlFmNnEJh6uAAznLW65xFo80Ei1GHNodlXVaihCIxuOZCqLdmDuAWk5RC4cFVTGkcwFnOepw+tNdZAABtdRb0T80hGltfwcmZZQkhIKVQIjGBYCSm9mWxdYYDOMtJLCbQ65xFR50VANDusCAUiWHUM7/KV1ZcE/FTmPVpptEns8b7oXAenCmNAzjLyej0PALhGDrq5RW4GcD6q0RxyrMws0yhADwXkykvm5mYTxCRk4jOJt1WQ0S/I6Ke+J8b1L1MVirkE5gd8RRKq136s3ed5cEXhhlnkULhhlZMJdmswL8H4N4lt30RwAtCiA4AL8Q/ZuuA3ANFzoFvMBtQazasuxX4uDcIvZawodKQ8b7yCpxTKExpGQO4EOJVAO4lNz8I4Pvxv38fwHuUvSxWqnqds3BYjahOClxtDsu6C+BObwB1VhM0Gsp4X6tRD4BTKEx5+ebA64UQVwEg/mddujsS0aNE1E1E3ZOTk3k+HSsVPc7ZRPpE1lZnWXeHeSZ82Z3CBACzURp4zCkUpjTVNzGFEI8JIbqEEF0Oh0Ptp2MqEkKuQFkSwB1muOdCcM+FVunKim98ZuVRasksPBeTqSTfAD5BRI0AEP/TqdwlsVI17g1gNhhJ5L9lbfGP11MaxeldeZhxMjmFwgGcKS3fAP4rAI/E//4IgF8qczmslC1sYFoX3d7uiAfwdVKJMheMwBeMZB3ATXoNNMRzMZnysikjfArAYQDbiGiEiD4G4O8A3E1EPQDujn/M1ji5B4pcAy5rrq6AUadZNytwpy9+iCfLHDgRST3BeQXOFKbLdAchxAfSfOpOha+Flbhepw8bKvWoNS8undNoCK2O9bORme0otWRWk54DOFMcn8Qsc//51lX8r2fOF+W5eiakI/REy0vn2hzmdXOYZ6EPSvYB3GzUcgqFKY4DeBnrn5zFn/z4NL7/xiBiKjeTEkKgxzmL9iXpE1mbw4IrHj8C4aiq11EKcjmFKbMYdZhbp6PnmHo4gJepcDSGz//4NObDUURiAjPzYVWfb2o2hJn58LISQllbnQVCAIOutZ9GmfAGUWnQwmLMmIFMsJj0fBKTKY4DeJn6xou9OH1lOjHWbGo2qOrzLfRAsab8/EIlSukE8O++PoC7/+EVxU9ASqPUTClTSelYjFrOgTPFcQBXyBV38dIHJ4Y9+MZLvfi9vc34o/1bAACTKgfw3jQVKLKtdjOISqcWXAiBH7w5hB7nLL7xUq+ij+3MchZmMp6LydTAAVwBvkAYd3/tFTxxaED155oLRvD5H51CQ5UJf/XgLtitUiCZmlX3FGTPxCysRh3qrKkDV4VBi+bqipLZyDw35sXA1Bwaqkx4/LV+9Cv4i2Uih0M8MotRz5uYTHEcwBVwpN+NQDiGgSKU0f3Nr89j2O3HP/zhtagy6WG3xAO4T/0USnu9ZcW0QSk1tXr69Bh0GsIPPn4jTDot/vrp84qMNBNCJFIoubAYtZgN8Vg1piwO4Ao41DcFQMqNqul35yfw1NEr+MStbbiptRYAUF2hh1ZDcM2pn0JJt4Epa3NY0D85p3pFTCaxmMAzZ67i1k4H2uus+OO7OvDK5Uk8f6Hwjg8z82GEIrHcA7hJByEAf2jtV+mw4uEAroBDvVIAvzqjXgB3+gL4wk/PYGdjFf773Z2J2zUaQq3ZgCmfeikU91wIU7OhtBuYsvY6C+bDUVxV+RdZJieveDA6PY8Hrm0EADxycws66iz4n8+cK3ifIjFKLeccOPdDYcrjAF4gpy+AyxOzMOo0uDo9r8pbZCEEvvCTM5gLRvBPD10Hg27xy2a3GFWtQpHz2ulqwGVtDvOi+6+Wp09fhVGnwV076gEAeq0Gf/3uXbjinsdjr/YX9NgLNeC5rcDllrJcSsiUxAG8QG/0ugAA9+xqwFwoCp8KK6wnjwzjpUuT+OI7t6OjfvkquNZiKEoAz5hCqVv9plbRePrkju11sJr0idtvbrfjvmsa8c2XenHF7c/78cdzmIWZzMpzMZkKOIAX6FDvFKor9bhzuzTTYlzhNErf5Cz+5tfncUuHHY/sb0l5H4fFqGoVSo/Th0qDFk22ihXvV2s2wFahX9WNzCP9LkzNBvHAtU3LPvc/7tsBDRH+n19fyPvx5WHGdZxCYSWAA3gBhBA41DuF/a21aN4gBTcl8+DhaAyf/9EpmPRafPUPrk07vstuNWJyNqhahUOvcxbtdZaM48OICO11q1uJ8vSZMZgNWty+bfmQqObqCnzmjnb89tw4XuvJbzrUhDeI6ko9THptTl/HKRSmBg7gBRh0+TE2E8CBdjsabdJb6qvT84o9/tdf6MGZkRn87/des2LO1W4xIBSJqZK+AaQa8KVDHNJpc5hXrSthKBLDb86O4+6d9agwpA6wH79lK7bUVuKvfnUOoUgs5+eY8AZyTp8APBeTqYMDeAFej1efHGi3o85qApFyK/DjQ9Jpy9+/fiPedU3jiveVa8FdKqRRvIEwxr2BHAK4BZO+oOq9WVI51DuFaX84ZfpEZtRp8ZUHdqJvcg7ff2Mw5+eY8AZQb8s9gPNYNaYGDuAFeKN3Cs3VFWiprYRBp4HdYlQsB/6NF3vgsBrxV+/emfG+icM8KmxkLmxgrlxCKGuTe6KsQhrl6dNjqDLpcEvHyrNX79hejzu31+Efn7+cyGlna8IbRH2a06gr4cHGTA0cwPMUjQm80efCzW21idOJjTaTYjXQQy4/urbULKqkSEfN05i9E9lVoMhWqxIlEI7iufMTeOfuxmVllqn8xf07EY4K/O1vLmb9HNGYwORs7sfoAWnlb9BqOIAzRRUUwIno80R0jojOEtFTRJT7d3aZOj/mxcx8GG/vsCdua7SZFMmBCyEwOj2f2BjNxG6RJuSosQLvcfpg0GmwqaYyq/tv2lABg1ZT9Dz4y5ecmA1GVkyfJGuxm/Hora34+clRHBt0Z/U1rrkgojGR8yEemcWk434oTFF5B3Aiagbw3wB0CSF2A9ACeEipCyt1cv57f1tt4rZGW4UiKZTJ2SCCkRg2ZhnAa8wGEAGTKuTAe52zaHNYoM1QgSLTaTVosVcW/TDP06evwm4xYF9rTdZf86nb29BkM+Evf3kO0TTH/6MxgYGpOfz27FV866U+ALkf4pGZuaUsU1j2HenTf30FEYUBVAIYK/ySysMbfVPYVm9FXVJFQoPNBF8wAl8gnFXqI50Rj7SKzzaA67QabKhU5zBPj3MW12/ekNPXtDksuDTuU/xa0pkNRvDCxQn8Ydcm6LTZr0kqDTp8+b6d+PS/n8C/HxnCPbsacHHch8vjPunPCR96nD4EwlK1CpGUSrpmoy2v67QYeS4mU1beAVwIMUpEXwUwDGAewHNCiOeW3o+IHgXwKABs3rw536crKYFwFEcH3PjgTYv/PXIp4YQ3UFAAH40H8Obq7NIWgJRGUToH7g9FMOKZx/u7NuX0dW0OC547P4FQJJZVPrpQL1yYQCAcyzp9kuxd1zTg5rZa/MUvz+EvfnkucbvDasT2Bis+dNMWdDZYsb3BivY6CyoN+a95rEZOoTBl5f3dSEQbADwIYCuAaQD/QUQfEkL8IPl+QojHADwGAF1dXWuil+aJYQ+CkRje3m5fdHtj/KTi2HQA7VlWbaQir8CzzYED0kama07ZFIo8XSfdEId02ussiMYEht1zBf1/yNbTp8fQaDPhbTm+UwCkw0f/7+/vwROHBtBSa0ZnvRXbGqyoMRsUv06zUat633a2vhSSQrkLwIAQYhIAiOhnAG4G8IMVv2oNONQ7Ba2GcOPWxflWeQVeaB58xOPHhkp9TjMX7RYjTo9MF/S8S8lj1LKtAZfJpYS9TvUD+Iw/jFcuT+IjN7dkPCmazqaaSnzlgV0KX9lyFpMeQ678+7AwtlQh72+HAewjokqS6ujuBJB/k4kycqjXhes2VS9Lk8j9MQo9zDPimcfGDdmnT4B4R0KFUyg9zlnoNIQtteacvq413pWwGLXgz54bRzgq8kqfFJvFqFPttCxbn/IO4EKIIwB+AuAEgLfij/WYQtdVsmbmwzgzMo0DS9IngFTra7cYMO4trJRwdHoezdXZp08AqSPhXCiKeQUHBvRMzGKr3Qx9DhuDAGA26tBoMxWlFvzpM2PYUluJa5rz21gsJotRyzlwpqiCdpiEEF8RQmwXQuwWQnxYCKHuWJgScKTfhZgADiSVDyZrtFVgbDr/FbgQAiMef9YVKDKHCqcxe52+nPPfsmKMV5uaDeJQ7xQe2NOU04T41WIx6jEfjqYtWVRCr3MWJ4Y9qj0+Ky18EjNHh3qnUKHXYm+aDbMGm6mgHLhrLoRAOPsacJndKm26KTWdPhCOYtjtzzuHLXUlnFN1BuRv3rqKmEBZpE+A4vRD+fOfnMZHnjiq6DsxVro4gOfoUJ8LN26tSVse12gz4epM/imURAlhHjlwQLmGVgNTc4iJ7I/QL9XmMGM2GEmMIFPD06evorPegm0N6le6KMGicj+UCW8AJ4an4Q1E8MyZdXMkY13jAJ6D8ZkAep2zy8oHkzXYTPAGInm3Dc31EI9M6YZWPXITqwJSKIB6G5lXZ+ZxdNCNB/aUx+obWBjqoFZL2efOTwCQBms8eWRYledgpYUDeA7eiE+fv7k9df4bQGJqTb6VKCMeqcwslxpwQNrEBJRraNU74YOGgK323CpQZImmVioF8F+fuQoAuL9M0ifAQgpFraEOz50bx1a7GZ+6vR2nrkzj7OiMKs+TyrQ/hM/98CRcKo72Y8txAM/B671TqDEbsKOhKu19GgqsBR+dnkeVSYeqHE9yGnVaWE06RVfgLbVmGHW5TZ6R1VmNsBp1qlWiPH16DNc02/L+BbMa1EyhzPjDONznwjt2NeB912+EUafBvx8t3ir8ufMT+MWpMbx40Vm052QcwLMmhMAbvS7sb6td8cBIYjJPnnnwfGrAZUrOxuxxziZW0fkgIrTGNzKVNuSaw+mRGTxw7cqDLkqNmimUFy9NIBITeMeuetgq9Xjg2ib88uRo0XqvHBuQOjqeG/MW5fmYhAN4lvom5zDuDayY/wYWOtUVkkLJNf8ts1uMilShhCIxDE7N5b2BKWtzmFXpSvhMPH1yXxnlv4GkKhQVUijPnp1AfZUR126sBgA8fNNmzIWi+MXJUcWfKxW5Je+5seKlbRgH8KzJ+e8DbSsHcJNei1qzIa8ALtWAZ98HfCm7VZmOhEOuOURiIu8NTFmbw4Jxb0DxVeALFyZw7abqnA87rTZLvBGW0qcxA+EoXrk8iXt2NiTeHV63qRo7G6vw5JFhVUs5AcDpC2DQ5UeFXotzY17EVKxzZ4txAM/S6z1T2FRTgc21mdMbUi147imUaX8Y/lA07xSK3WJUpIywJ8cxaunIlSj9Cm5kBsJRvDU6k1Pf71Ihj1VTOoXy6uVJzIejeMeuhsRtRISH923GhatenLwyrejzLdU9KB0c+r3rm+EPRTHgWp2h1usRB/AsRGMCh/tdGVffMqkWPPcVeL4lhDK7xYiZ+XBe09aT9TpnQbQQgPPVrkIlyqkr0whHBW5sKb8ArtNqUKFXfqjDs+cmYKvQ46Ylv9QevK4ZZoMWT76p7mbm0QE3KvRavP8Gqe0w58GLhwN4Ft4anYEvEEnZ/ySVRltFngFcKiEsJIAD0uivQvQ4Z7FxQwUqDPlVoMi21FZCpyFF8+DyZtnbtuTeOrYUmI06RcsII9EYXrg4gTu31y3rWWMx6vCevc145swYpv3qtbE9NujG3s3V2NFYBYNWg3NFLF9c7ziAZ+FQfHzazWn6nyzVYDNhZj4Mfyi3H9TR+DzNjTkMcki2UAte2A9rz4Sv4PQJAOi1GmyurUz0FVfCsSEPttVbUV2pfL/uYrCadIqmUI4OuDHtD+OepPRJsodv2oJgJIafnlBnM9MXCOPCVS+6Wmqg12qwvdGKs7yRWTQcwLNwqHcKOxqrUGvJbphtvn3BRzzzsBp1qKrIr027EqcxI9EY+hWoQJEp2dQqGhM4MeRBV0t5rr4BaVWsZArl2XPjMOk1uK3TkfLzO5uqsHdzNZ48MqTKZubxIQ9iAomU1q6mKpwd9aq+ccokHMAzCISj6B7ypO0+mEq+h3lGPH40b6jIu7Oe3JGwkFLCK555hCKxnIc4pNNeZ8Ggaw6RaGF5eQC4cNWL2WAEN5Rh/ltmVrClbCwm8Oy5Cdza4Vgx3fXwTVvQPzmHN/vdijxvsmODbmg1hL2bqwEAu5psmJkPJ95NMnVxAM+ge9CDUCSGAx3Z5b+BheP0Y3mswPPNfwMLHQkLWYHLJycLOcSTrM1hQTgqcMVT+A90d7zW+Iat5RvAlRxsfGZ0BuPewKLqk1Tu39OIKpMOTx4ZUuR5kx0b8GB3UxXM8elRu+N92c+O8kZmMXAAz+D13inoNJRT1cPCCjz7oCWEwGgBpzABacp6pUFbUCnhYLwErCXHKTzptMWn8yixkXls0IMmm6ns6r+TWU3KpVCePTcOrYZw5466Fe9n0mvxvrdtwrPnxjGp4NSmYCSKUyPTi94RbW+wQqshPtBTJBzAM3ijbwrXb96QWGFkw6TXYkOlPqdKFO98BL5gpKAVOBAfrVbACnzY7YfVqMOGytx6saTTGi9FlOdr5ksIgWOD7rJefQPxFIqCAXxfa01WG7ofvGkzwlGB/zh+RZHnBoC3RmYQisQWvSYmvRbtDktRG2mtZxzAVzDjD+Ot0ZkVuw+m02CryCkHfkXuQljg6tJuKew05rDbj821lYpNuLFV6NHmMOONXldBjzPs9sPpC6KrjPPfgHIplF6nD/2TcxnTJ7L2Ogv2t9bi348MK3ZS8mg8pdW1pKRzV3MV14IXSUEBnIiqiegnRHSRiC4Q0X6lLqwUnBj2QAjgpq25B/AmmymnHHiihLCAFAoA1FqMBZURDrv82JLFadNc3LWzHm/2u+ANhPN+jGPx037leIAnmdWkQygSQzBS2MScZ89Jvb/v2ZldAAeAh/dtxohnHq/2TBb03LJjA260OczLqrN2N9ng9AXh9BY23JtlVugK/J8A/FYIsR3AtVhjU+lPDHug1RCu3ZT7wNxcj9MXegpTVkgKJRoTuOLxY3ONsi1a795Rj0hM4JVL+QeOYwNu2Cr0ipU3rhazQT5OX2gAH8d1m6oT+y3ZuGdnA+wWZYY9RGMC3UMe3JgipSVvZPIqXH15B3AiqgJwK4DvAoAQIiSEmFboukrCiWEPtjdYUWnIvS670WaCxx9GIJzdD+qIxw+zQYvqAnPPDosBbn8or7K9qzPzCEeF4ivwvZs3oMZswAsXJvJ+jGNDbnRt2bBiK99yYDEV3lJ2bHoeZ0Zmsk6fyAw6Df6waxNeuDBR0Ng/ALg07oMvkLqkc0ejdAiM8+DqK2QF3gpgEsC/EdFJInqciJYt3YjoUSLqJqLuyUll3roVQzQmcGp4GtenGV6cSUO8lDDbPPhovAthoblnu9UIIQCPP/d0xbBLysNvqVE2gGs1hDu21+HFi06E8/jF4poNon9yruzz34B0kAcobCrPc+fGAQDv2FWf89d+4MbNEAB+eLSwzUy5fWyqAG416bHVbuYVeBEUEsB1AK4H8G0hxF4AcwC+uPROQojHhBBdQoguhyP1abFS1OP0YS4UxfVbqvP6+qb4W9uxLFc6hQxySFbIacwhtxTAs+m4mKu7dtTDG4gkOtflQs5/31DGJzBlcgAvZCPz2XMT6KizJCp8crGpphK3dTrww2PDef0ylR0ddKPRZkqb8tvVVMVH6ougkAA+AmBECHEk/vFPIAX0NeHE0DQAYO+mfFfguZ3GLGSQQ7KCArjLD72W0GhTvs76lg47DDoNns8jjdI96IZBp8E1G3Pfiyg18lCHfFMonrkQjg66c06fJHv4pi2Y8AbxwoX8xp8JIdA96MYNLTVp3zHubrZhxDOvahMtVkAAF0KMA7hCRNviN90J4LwiV1UCTgx7UGM25J0PbrBlP5nHGwjDG4gockDFbsn/NOawew4bN1RCq0Ke2WzU4UBbLZ6/MJFzn4xjg25ct7E67/mcpSSRQskzgD9/YQLRmCgogN++zYFGmynvk5lX3POY8AZXrMnf1STNjeU0iroKrUL5LIAniegMgOsA/O+Cr6hEnBj24PrN1XnnpCsNOtgq9FmtwEc9ypQQAkiUdOVTSjjk8mOzwvnvZHftrMeQy5/TqUx/KIKzY17csLX80ydAUgolzxz4s+cm0Fxdgd3N6QdrZ6LTavDBGzfjtZ4pXBrP/YCVXP+9Uknnria5EoXTKGoqKIALIU7F89t7hBDvEULknuAsQdP+EPon57A3zw1MmTTYIXMOXKkSQgCoMulg0GpyXoELIVSpAU9253Zp0+13OaRRTg5PIxoTa2IDEygshTIXjOC1nkncvbO+4M3uD+/fgkqDFt9+uTfnr82mpLPGbEBzdQX3RFHZmjqJOe0P4TdvXS34ceQRVHKHtXxlO5lHHuSQ7yzMZEQEu8WQc0dCjz8MXzCi6gq8wWbCno02PH8++wB+bNANovId4LBUpV4LovxSKK9cnkQwEisofSKrrjTggzduxtNnruJKfPM6W8cGsyvp3MkbmapbUwH88dcG8MknTxTcOOnkkAcaQmLCd76yPU4/6pmHSa9BrVmZIQV2a+6zMYfiTay2KNTEKp27dtTj5JXprJsqHRt0Y3tDFapMyvRmWW0aDcFs0OWVQnn23Dg2VOoVq8b5+C2t0BDwL6/2Zf01k74g+qfmsupJs7vJhoGpOcVngLIFayqAvx6fnPPypfx212UnhqexvaEqpwZWqTTZTHDNhTIe5pFLCJXqP5LPaczh+CpMzRQKIAVwIYCXLmZ+jcLRGE4OT6+J8sFkFmPuU3lCkRhevOjEXTvqodMq82PbYDPhfW/biB93j8Dpy65aqnuF+u+ldjdXQQipjztTx5oJ4N5AGGdGpgEALxdwZDsaEzh1ZTrv+u9kciXKRIaeECPTypQQyvJpaDUUP8SjZgoFkE7pNVdXZJUHv3DVC38oWtYDHFKx5NFS9nC/C75ARJH0SbJP3NqGSDSG774+kNX9jw16YNJrcE1z5pLOhd7g6qZR/tcz5/H4a/2qPkepWjMB/Ei/GzEBXLvRhqMD7rzftvU6ZzEbjORd/51MrqfOlAcf8cwr2uPabpFSKLl0nRty+VFfZYRJr26pHhHhrh11eK1nMuM7k6MD2a/2yonZqMs5B/7suXFUGrR4ew6DRbLRYjfjXdc04sk3hzEzn/n07rFBN67bVA2DLnPoqLMaYbcYcFbFUsKzozP47usDqs38LHVrJoAf6p2CSa/B5+7qRCgaw+G+/NqXnhiWCmmuV2DTLJvDPLPBCKb9YUVKCGW1FiMiMZHVD6Rs2D2HLQo3sUrnrp31CIRjiWHR6XQPerCppiKnhk3lwJpHCuW1nknc0mFX5Rfspw62YzYYwf85PLji/WaDEZwbm8m6IyQRYVeTTdUV+Ndf7AEA9E3OIqpQm9xysmYC+OE+F25oqcHN7bUwG7R4Kc88+Ikh6QBPiwK54MYsDvOMKlhCKMvnMM+Qy6/KEfpUbtpaC4tRt+KpzMQAhy1ra/UNxAcb57CJOeMP44p7Htcp8K4wlZ1NVbh9mwNPHBrEfCj9u6IT8QHGuZR07m6uQq9zNuumbrm4OO7Fs+cm0Go3IxSJJfZx1pM1EcAnfUFcmvDh5jY7jDotbm634+VLk3lNxj55ZRp7N+V/gCeZ2ahDlUm3Yi24kiWEMkfiOH12lSjzoSicvqDiTazSMeg0uG2bA89fcKZN8wxMzcE1Fyr7CTypmHOcTC8fhpFPN6rhU7e3wz0Xwg+PpW81e2zQDQ3l9u50d5MNkZjA5YnCJjKl8vUXe2Ex6vCXD+wEAFWeo9StiQD+Rp/0VvxAfHLOwW0OjE7Po28yt3LCGX8Yvc7Zguu/kzXaKlZegU+rsAK35tYPZVjFJlbp3L2jHpO+IM6keXu90O1ubVWgALnPxTxbhAB+Q0sNbmypwb++2o9QJHWTq6MDbuxqsiVOk2ZDPpGp9IGeXqcP//nWVTxy85bEOwIl5q6Wm7URwHtdqDLpEt8sB7dJQ15fuphbNcrJK/H8d4EnMJNJgx3SB/ARzzyMOk1i1ayEXBtaFasGPNnBbQ5oNZT2UM+xQQ82VOrRlkfHvVJnia/As32HeHbUiyabadnkG6V98vY2jM0E8MtTyzcEg5EoTl2ZznlDeVNNBawmneIHer7xYi8q9Fp87O2tsBh1aK6u4BV4uTrUN4V9rbWJJkzN1RXorLfg5cu55cFPDE9LB3g2VSt2bU3VK5/GHPH40VxdeB/wZNUVemg1lPMKvFgpFEA6CXhDy4a0efBjg250rdDtrpyZjTpEYwKBcHbtXM+OzWBXFmV7hTrY6cDOxip8+5W+ZRuCZ0e9CEZiuDHHnjREhN1NNkWbWg1MzeFXp8fw4X1bUBM//NZeZ0HPBK/Ay84Vtx8jnnkcaF9cXnVwW13O5YQnhz3YpsABnmQNVRWYmg2mnYE4Eh/koCSNhlBjNmTd0GrI5YfVpCt4GlCu7tpRj4vjvmVHuZ3eAIZc/jWZPgEW+qFkk0aZDUYwMDWH3U3qB3Aiwqdub0P/5FxiaIRMTmnl05NmV1MVLlz1FtR/PNk3X+qFQafBx29pTdzWWW9Zl5UoZR/A5VK0A0smxx/sdCAcFRlL1WSx+AEeJfPfwEIlitObejU8qtAgh6VyOY055JaaWBV7tXv3Tqm51dJVePeQPMBh7W1gAlIZIZBdAL9w1QshUFD3wVy8c3cjttrN+NbLfYtSPMcG3Gh1mBPpuVzsbrYhFInlvCeVyrDLj5+fHMUHb9wCh3XhWjrqrAhGYjn3dSl35R/A+1yosxqX5Uq7WmpgNmjx8uXs8uC9k7PwBSKK5r+BlfuC+0MRuOZCim5gynI5jTnsKl4NeLIttWa011mWBfCjA26Y9JrEnsZaY86hpaxcQ727CCkUQBp/94lbW/HW6Axe65EWP7H4AON8SzrlXz5KbGR+6+Ve6Rpva110e0e99PO/3vLgZR3AhRA43DeFm9tql60eDToNDrTb8UqW5YQnhuQNzGpFr7GpWg7gy0sJ1agBlzksxqzKCKMxgRHPfFErUJLdtaMeR/rdiw4ddQ+5sXfThqxO+5WjXMaqnR31wm4xos6q7gZmsvde34z6KiO+FW81e9npw8x8OO+Szq12Cyr02oJ7g494/PjpiRE8dMMm1FctPtzVHm9t27POKlHK+ifk0oQPU7Mh3Nye+njxwW11GJ2ez+pFPTk8jepKaRirkhpWOE4/okIJocxulVIomX55jU3PIxITRd3ATHb3zjpEYgKvxN8p+QJhnB/zrtn8NyCVEQLZBfBzYzPY3VxV1PSWUafF/3VLK97sd+P4kAfHBjIPcFiJVkPY0WjFuQJX4N95Reqa+F9va1v2OatJjyabCT28Ai8fb/RKx+WXbmDKDm6Thihn053wxLBHsQM8ySxGHaxGXcpSwhEFJ/EsZbcYEIzEMgaJ1agBT3bdpg2oNRsS5YQnh6cRE1iTB3hkiRRKcOVWB4FwFD3O2aJsYC71gRs3o7pSj2+/3Iujgx7UVxmxqSb/hcbuZhvOjc3k1J8n2fhMAD8+NoI/6NqEpjR9gzrqrbi8zipRCg7gRKQlopNE9IwSF5SLN/qm0FJbmbYRVFN1BbbVWzN2J5yZD6PHOat4/lvWkGYyz4jHD4NW2RpwmT3L05hyF8Ji1oAn02oId2yvw0uXnAhHY4nTfoVOQyplCymUlY+XXxz3IRoTRdvATGY26vCRm1vw/AUnXr7oXHGAcTZ2N9kwF4piMH7mIFffeaUPMSHwyRSrb1lH3fqrRFFiBf7HAC4o8Dg5iURjONLvTps+kR3c5sCxQfeKK9FT8Qk8SjSwSqWxOvVpzFHPPJqqTRknm+SjNsvDPEPuORi0GjRUrV7DqLt21sMXiODYoBvHBt3Y2VSV02m/cpNIoWTYxJQ3MFdrM/cjN7eg0qCFLxjBjQW+I9rVnP+QY6cvgKeODuP3rm/GphVSfZ31668SpaAATkQbAdwH4HFlLid7Z0Zn4AtGcKAtUwCvy1hOeHLYA1L4AE+yxqrUh3lGVCohBJIaWmWYfDPs8mNjTYUqk+izdUuHHQadBr95azyv037lxqjTQKuhjCmUc2MzsFXoVdkjyUZ1pQEP37QZAAoO4B11Vui1lNeJzH99tR+RmMCnb29f8X7t9etvI7PQFfg/AvhzAGkr9InoUSLqJqLuycn8By0s9UY8IO9rXfkbq6tlAyxG3YpplBPD09hWb1Vt1ddgM2FqNrisx4TSfcCTObJdgbv8q7aBKas06PD2djt+1H0FgXBszQdwIopP5Vk5hXJ21Fv0DcylPn93J/71j7qwvaGwNI5Bp8G2htw3Mqdmg/jBm8N48NqmjGk+ecjyeiolzDuAE9H9AJxCiOMr3U8I8Vh8cn2Xw+HI9+mWeaPPhR2NVRn7Q+i1Ghxor8XLl5wpKzJiMYGTwx5Vc66NNhOEwKKxVYFwFFOzQdVWVzVmA4hWzoELITDs9q9a/jvZXTvqE7/gutZwBYrMYtTBt0IKJRSJ4dK4b1U2MJNVGnSJA1eF2t1kw9mxmZy6hD7+2gACkSg+fcfKq29AqkRptJnWVVOrQlbgBwC8m4gGAfwQwB1E9ANFriqDQDiK7iEPDrTVZr4zgNu31eHqTCDlDnVf4gBPtcJXuaCxenkpYaILYQE7+yvRaTXYULnyYR73XAizKk+iz9adO6QGZC21laizrq0BDqlIDa3Sp1B6nD6EorGi9EApll1NVZj2hzGWxaBvAPDMhfB/Dg/i/j1NWTc1kypReAWekRDiS0KIjUKIFgAPAXhRCPEhxa5sBceHPAhFYmnLB5e6bYVywpPD0wDU28AEUg92ULOEUJbpNOZQkQYZZ6O+yoT7rmnEe/Y2r/alFIXFtHIKRU417FaxhWyx7cpxRuYThwYwF4ris1msvmUddRb0OtdPJUpZ1oEf6p2CTkNZ1wo32iqwvSF1OeGJYQ9sFXpsVTGNsDBabaGUMDHIQaUcOCD3Q0mfQhl2lU4AB4BvPnw9PndX52pfRlFYMszFPDs2A7NBi5YSSG8pZUdDFTQEnMsQwE9dmcannjyOb77Ui3dd04DOemvWz9FZb0EwEkv8fOXjtZ5J/Pcfn8prIEyxKRLAhRAvCyHuV+KxsnGoz4VrN1XntOl4cFsdjg264Qssftt6YtiDvZurVSnlk1WZ9LAYdYtTKJ556DS07EiwkmozNLSSa8DVfBfAUpPGqqVPoZwdncGuJpuq35fFVmHQor3OknLIcSwm8OLFCbz/Xw7jPd88hNd6pvCJ29rwt+/dk9NzdMSDfSEHen7w5hB+dmJUkeZbaiu7Fbg3EMZbI9NZ579lB7c5EIkJHOpdGHbsDah7gCdZg82Eq9OLUyhN1eqW79kthhXLCIfcc2ioMqk+iZ4tt1IVSjQmcP6qN1E7vZbsarIt6okSjETx4+4reMc/voqPfq8bw24//u/7duDwl+7EF+7dDluOLY4XeqLklwePxQSOxFsH5DsYvZjK7rTEkX43YgIZD/As9bYtG2A16vDKZSfu3d0AADh9ZRpCKDuBJ51GmwlXvckB3K9q+gSQUihzoSjmQ1FUGJYH6eEiDjJmi1lWGKvWPzmLQDi26hUoatjVVIWfnxxF/+Qsnjs/gX87NIAJbxDbG6z42vuvxf17mqDX5r+urIpXouQ73OHCuBfTfumd0eF+Fz68vyXvaymGsgvgh3qnYNJrcu7brddq8PYOO166KHUnJCKcGJqOH+BR/welocqEyxMLOfgRzzxu61SurDKV5FrwVCfYhtx+HFT5Glhq8mDjWEwsS5PIh12K1UK2mOR/0z1fexWRmMCB9lr8f++7Frd22BWrd2+vs+S9ApdX3ftba/Fmvzvl61NKyi6F8kbfFG5oqYFRl/vb/oPbHBj3BnApXmZ0YtiDzjorrCb1J9E0VlfA6QsiHI0hGJGmwKude7Zb46cxU+TB/aEIJn3BktnAXG/koQ7+8PI0ytlRL4w6Ddoca2cDU3ZNsw3bG6x45zWNeOazb8eTH9+H2zodih5W6qy3otc5m1fjrDf7XdhSW4n3Xt8M91wIl/P8RVAsZbUCd/qkWu737t2Y19ff1inVGr98aRKddVacHPbgvj2NSl5iWguHeRZOZCo9Sm2plRpaLXQhXHtBohxYkvqhLN2MPzs6gx2NVdAVkEooVWajDr/93K2qPkdHnQWBcAxXPLkdUovG89/3XdOI/a3SHtvhPlfBp1DVVFbfIfLbm6Xj07LVYDNhR2MVXrroRP/ULLyBSNG63iWXEsolTmr3uFhpOn2iC2EJHOJZj9K1lI3FBM6PeVelA+FaIVei5JoHPz/mhS8Qwf62WmyqqcTGDRUlv5FZVgH8jV4Xqky6grqzHdzmwPEhD169LPVSUfMEZrLkwzxqTuJJJk/sTlWJUmo14OuNNU1L2WG3H75gZE1uYBaLXImSa/rjcL8UE+TV9/7WWhwZcOfdw7wYyiqAH+qbwv622oJK7w52SuWEj73ajyqTDq327I7oFqoxPplnfCaAEc88tBpSvYWrSa+F1aRLvQJ3z6HKpEN1pUHVa2CpWdK0lF3LG5jFYqvQo6HKhN4cV+CH+1xodZhRF/+53N9Wi5n5MC6MFz7LUy1lE8CHXX6MeOazPj6fzvXxcsJxbwB7N28o2g5zlUmHSoMWY9MBjHj8aKgyFSXH6bAYMTW3PAc+5CqNJlbrldmQOoVydtQLvZYSQ3pZfjrqLTmtwCPRGI4NehKrb0AK4EBp14OXTQA/1Ce9vbk5xwM8S+m1GtzSKf0SKEb9t4yI0GAzYdw7j9Hp+aL1eLZbjKlTKG6uAV9NC3MxF6dQzo3NoLPemleVFVvQUZdbJcpbozOYDUYSQRuQ3jW31FbizX4O4AV7o8+FOqsx665kKzm4TapGuX5LdcGPlYtGmzTYQc1BDkvZrcsbWkWiMYx65nkDcxUlxqolHacXQuDs6AznvxXQWS9VoshN4zI5HA/S+1oXLxD3t0l58FJtjlUWAVwIgcN9UzjQrkyx/3v3NuPrH9ibcZqP0hptFbji9mPcG1C9hFCWqqHV2HRAmkTPK/BVs1CFspADH5sJwOMPcwWKAjrqcztSf7jPhc56S6JyS7avtRa+QGTR8f9SUhYB/NKED1OzoYLTJzK9VoMHrm0q+gmrRpsJU7MhCKF+BYrMbjFiZj68aBrQkFsaLLu5hnPgq8Wg08Cg0yxKoSRmYPIGZsHa67JvahWKxNC9JP8tS64HL0VlEcDlBlS59j8pNXItOFC8AF4bn43pmltIowyXUB/w9cy6ZKjDudEZaEhqu8oKY6vQo77KiJ4shjucGZnGfDi6KP8tq6syoc1hTqRYSk1ZBHCnL4D2OovqzZ/U1pgcwKuLlAOXD/P4FtIowy4/DLrVnUTP4v1QksoIz4550V5nSdl4jOWus96a1YDjw30uEAE3bU39Dn9/Wy2ODbgRjqYd/btqyiKAf+mdO/DbP75ltS+jYHItuIYWr8bVlAjgSSvwIZcfmzZUlHSTnvVAGqu2OIXCG5jKaY9P58lUiXK4Xzouv8Gc+kzE/lY75kJRvJXlJKFiKosADmBN9IWQV+ANVSYYdMX59yQ6EiaVEg6VyCDj9U5qKSulUJzeAJy+IOe/FdRZb8V8OJqYP5tKMBLF8aHU+W/ZvlZp8lcp5sELmUq/iYheIqILRHSOiP5YyQtbi2wVepj0mqJOwFnoSCilUIQQGHbNlcQg4/XOYlzoCX5ubO3NwFxtnfFKlJWGHJ8cnkYwEkuZ/5bVWozYVm8tyXrwQpaBEQB/IoTYAWAfgE8T0U5lLmttIiLsbrIV9Zh0pUE6ASrXgrvmQpgLRXkDswQkT+WRK1B2cgBXjFyJslIe/HCfCxoCbswwX3d/Wy26Bz2LqrlKQSFT6a8KIU7E/+4DcAHA+hgpXoCnHt2HL9+3o6jPaU+ajTnETaxKhsWkgy++iXl2bAZb7eai9KZfL+RKlJVW4If7XdjVZIOtYuX/7/taazEfjuL0yLTCV1kYRRKxRNQCYC+AI0o83lqm12pUnYOZSq1l4TTmMNeAlwxLUhnh2VEvdvHqW3HykfpUAuEoTg1Pr5g+ke1rrQFR6eXBCw7gRGQB8FMAnxNCLGvbRUSPElE3EXVPTk4ufwCmOrvFCFc8Bz7k8oMI2FRT3iWZa4HFqEMgHMOkL4jR6XnuQKiCjnoLeiZSV6IcH/IgFI2tuIEpq640YEdD1doK4ESkhxS8nxRC/CzVfYQQjwkhuoQQXQ4Hz19cDckplGGXH41VJm6WVALkfihH41PQuYRQeR116StRDve5oNUQbsiQ/5btb6vF8WEPAinG4K2WQqpQCMB3AVwQQvyDcpfElOawGOCeCyEaExjiLoQlQw7gcnUDp1CU17lCT5TD/S5c02xbNtIunf2ttQhFYjg5PK3kJRakkBX4AQAfBnAHEZ2K//cuha6LKchuNSImAPdcSOoDzvnvkiAPdTgy4EJzdUXagyQsfx1peqLMBSM4fSW7/LfsxtYaaAgldaw+76HGQojXAfBRvjIgn8YcdvsxNRvkFXiJkFd+lydm8Y5d9at8NWuTrVKPOqtx2XzM7iEPIjGRVf5bVmXSY3ezDW/2uYC7lb7S/JT/8UaWUW18ZXdy2AOASwhLhTnprTvnv9Uj9URZnEI53OeCTkN425bchrrsb63FySsezIdKIw/OAXwdsFulFfjxoXgA5xRKSZCn8gA8A1NNqXqiHO534dpN1Yt+iWZjX1stwlGR+FlabRzA1wE5hSJ/03EKpTQkb57t4iEOqumst8IfWqhE8QXCODs6k1P6RHZDSw20GkpMsF9tHMDXgSqTDgatBk5fELYKfcZTZ6w45NVfndWIOiu39lWLPJ1HPtDTPehBNCZy2sCUWYw67NloK5l6cA7g6wARwR4f7MD579Ihr8A5faKujrrFTa0O97tg0Gpyzn/L9rfW4szIDOaSxuGtFg7g64ScB+cuhKVDqyHsbq7C7dvrVvtS1rTqSgMcVmOilPBwnwvXba6GSZ/fYbb9bbWIxASODbqVvMy8cABfJ+Q8OK/AS8szn70FH963ZbUvY83rrLeg1+nDzHwY58byy3/LurbUQK+lnOrBhVBnqj0H8HVCLiXkChS2HnXUSePVjvS7EBPIK/8tqzBocd2maqkePININIafHh/BPV97FX2Tmce75SrvgzysvCRSKLwCZ+tQR70F/lAUPz0xAqNOg72bqwt6vP2ttfjGS73wBsKoStECOBSJ4acnRvCtl3txxT2P7Q1WeOfDKR6pMBzA14kmmwlEQKudV+Bs/emsl47U/+78BPa11hbczG1fWy3++cVeHBtw484dC6doA+EofnTsCr7zSh+uzgSwZ6MNf3n/Lty5vU6VGbQcwNeJ33/bRnTWW1HHk+jZOiRXosQECsp/y67fvAEGnQaH+1y4c0c9/KEInnxzGI+91o9JXxBdWzbg735/D27tsEPq+6cODuDrRKVBh5sU+MZlrBzJlSiTvmBB+W+ZSa/F9Zur8VrPFL75Ui+++/oA3HMhHGivxT8/tDc+AEL9VlEcwBlj60JHnQWzgQj2bKxW5PH2t9rxtecv4++fvYTbtznwmTs68q4tzxcHcMbYuvDJg20Ym56HQadM8d1DN27CzHwY793bjGs2rs5hLA7gjLF14ZYOZSeC1VeZ8JcP7FT0MXPFdeCMMVamOIAzxliZ4gDOGGNlqtCp9PcS0SUi6iWiLyp1UYwxxjIrZCq9FsA3AbwTwE4AHyCi1c3oM8bYOlLICvxGAL1CiH4hRAjADwE8qMxlMcYYy6SQAN4M4ErSxyPx2xYhokeJqJuIuicnJwt4OsYYY8kKCeCpzokua3orhHhMCNElhOhyOJStw2SMsfWskIM8IwA2JX28EcDYSl9w/PjxKSIayvP57ABKY5IoS8avS+nh16Q0FfK6pJz6QflOiiAiHYDLAO4EMArgGIAPCiHO5XmBmZ6vWwjRpcZjs/zx61J6+DUpTWq8LnmvwIUQESL6DIBnAWgBPKFW8GaMMbZcQb1QhBD/CeA/FboWxhhjOSink5iPrfYFsJT4dSk9/JqUJsVfl7xz4IwxxlZXOa3AGWOMJeEAzhhjZapkAzgRVRPRT4joIhFdIKL9RFRDRL8jop74n8WdX7TOEdHniegcEZ0loqeIyMSvSfER0RNE5CSis0m3pX0diOhL8YZzl4joHatz1Wtbmtfk7+Px6wwR/ZyIqpM+p8hrUrIBHMA/AfitEGI7gGsBXADwRQAvCCE6ALwQ/5gVARE1A/hvALqEELshlY4+BH5NVsP3ANy75LaUr0O8wdxDAHbFv+Zb8UZ0TFnfw/LX5HcAdgsh9kA6M/MlQNnXpCQDOBFVAbgVwHcBQAgREkJMQ2qW9f343b4P4D2rcX3rmA5ARfwQVyWkk7f8mhSZEOJVAO4lN6d7HR4E8EMhRFAIMQCgF1IjOqagVK+JEOI5IUQk/uGbkE6rAwq+JiUZwAG0ApgE8G9EdJKIHiciM4B6IcRVAIj/WbeaF7meCCFGAXwVwDCAqwBmhBDPgV+TUpHudciq6RxT3UcB/Cb+d8Vek1IN4DoA1wP4thBiL4A58FvzVRXPqT4IYCuAJgBmIvrQ6l4Vy0JWTeeYeojoywAiAJ6Ub0pxt7xek1IN4CMARoQQR+If/wRSQJ8gokYAiP/pXKXrW4/uAjAghJgUQoQB/AzAzeDXpFSkex1ybjrHlENEjwC4H8DDYuHQjWKvSUkGcCHEOIArRLQtftOdAM4D+BWAR+K3PQLgl6tweevVMIB9RFRJRATpNbkAfk1KRbrX4VcAHiIiIxFtBdAB4OgqXN+6Q0T3AvgCgHcLIfxJn1LsNSnZk5hEdB2AxwEYAPQD+C+QfuH8GMBmSAHlD4QQSzdzmEqI6K8BvB/S28GTAD4OwAJ+TYqKiJ4CcBBSe9IJAF8B8AukeR3ib+E/Cul1+5wQ4jfLH5UVIs1r8iUARgCu+N3eFEL81/j9FXlNSjaAM8YYW1lJplAYY4xlxgGcMcbKFAdwxhgrUxzAGWOsTHEAZ4yxMsUBnKmOiKJEdCrpv5bVvialENFeIno8/vePEJEgojuTPv/e+G3vi3/8MhF1JX2+JbmDXYrHdxDRb9X8N7DyxQGcFcO8EOK6pP8G5U+QpJy/D/8HgK8nffwWgA8kffwQgNP5PrgQYhLAVSI6kO9jsLWrnH9wWJmKrzovENG3AJwAsImI/oyIjsV7J/910n2/HO+Z/Hy8B/mfxm9PrGSJyE5Eg/G/a+N9mOXH+kT89oPxr5F7zD8ZP1EKIrqBiN4gotNEdJSIrET0WvwwmXwdh4hoz5J/hxXAHiFEcoB+DcCNRKQnIguAdgCnsvz/8njSu5RJIvpK/FO/APBwlv972TpS0FR6xrJUQUSn4n8fAPB5ANsA/BchxKeI6B5Ix4lvhNTo51dEdCukJmYPAdgL6Xv1BIDjGZ7rY5A6Jd5AREYAh4joufjn9kLqwTwG4BCAA0R0FMCPALxfCHEs3sp4HtIp4I8A+BwRdQIwCiHOLHmuLgBL0x8CwPMA3gHABunY9NYl93mSiObjfzcAiAGAEOLjAEBEWwA8C6nHNAB0A/ibDP9utg5xAGfFMC+EuE7+IJ4DHxJCvBm/6Z74fyfjH1sgBXQrgJ/LfSSI6FdZPNc9APbIOWdIQbQDQAjAUSHESPyxTgFoATAD4KoQ4hgACCG88c//B4C/IKI/g3Tk+XspnqsRUtvjpX4IafiFDcCfQEqzJHtYCNEdf54WAM/InyAiE4D/APAZIcRQ/GYnpA6QjC3CAZytlrmkvxOAvxVC/EvyHYjoc0jfZjOChRSgacljfVYI8eySxzoIIJh0UxTS9z+leg4hhJ+Ifgephe4fQlptLzW/5Lnlrz1KRLsh/eK6HM/UZOs7AH4mhHg+6TZT/LkYW4Rz4KwUPAvgo/GcMYiomYjqALwK4L1EVBHPNz+Q9DWDAN4W//v7ljzWJ4lIH3+sTpKGgaRzEUATEd0Qv7+VpIlDgJRG+WcAx9I06LoAKcedypewfOW9IiL6NACrEOLvlnyqE8tTNYzxCpytPiHEc0S0A8Dh+Gp1FsCHhBAniOhHkDYBhyBtEMq+CuDHRPRhAC8m3f44pNTIifgm5SRWGPMmhAgR0fsBfJ2IKiCtdO8CMCuEOE5EXgD/luZrLxKRjYisQgjfks/l013uTwGEk/YLviOE+A6A2wH8Oo/HY2scdyNkZYOI/gpSYP1qkZ6vCcDLALYLIWJp7vN5AD4hxOMqXserAB4UQnjUeg5WnjiFwlgKRPRHAI4A+HK64B33bSzOrSt9HQ4A/8DBm6XCK3DGGCtTvAJnjLEyxQGcMcbKFAdwxhgrUxzAGWOsTHEAZ4yxMvX/A6v292+218HXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plt.plot(y)\n",
    "#plt.plot(obs)\n",
    "plt.plot(abs(np.abs(y) - np.abs(obs)))\n",
    "plt.xlabel('Frequency (MHz)')\n",
    "labels=([\"60\", \"80\", \"100\", \"120\"])\n",
    "plt.xticks([0, 10, 20, 30], labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "welcome-border",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}