{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6b1727a0",
   "metadata": {},
   "source": [
    "## Importing libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "fe9aea2d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/cbrun/Documents/UFR/M1_TURB/TP/DATA_cylindre\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import netCDF4 as nc\n",
    "from scipy import stats, arange, misc, signal\n",
    "import netCDF4 as nc\n",
    "import datetime as datetime\n",
    "import sys\n",
    "Path2 = \"/Users/cbrun/Documents/UFR/\"\n",
    "sys.path.append(Path2)\n",
    "import os\n",
    "from pathlib import Path\n",
    "print(Path.cwd())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f867714f",
   "metadata": {},
   "outputs": [],
   "source": [
    "Path2='/Users/cbrun/Documents/UFR/M1_TURB/TP/DATA_cylindre/'\n",
    "sys.path.append(Path2)\n",
    "import Fonction_structure\n",
    "import Autocorrelation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ede1d056",
   "metadata": {},
   "source": [
    "## Reading the data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7f58c36",
   "metadata": {},
   "source": [
    "## NETCDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c19849d1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'netCDF4._netCDF4.Dataset'>\n",
      "root group (NETCDF4 data model, file format HDF5):\n",
      "    dimensions(sizes): time(37632)\n",
      "    variables(dimensions): float32 u(time), float32 v(time), float32 w(time)\n",
      "    groups: \n",
      "[10.733789 11.020116 10.271388 ... 10.222719  9.902984 10.832353]\n",
      "[ 0.09184415  0.29724777 -0.02563373 ... -0.08954611 -0.10164868\n",
      "  0.06464213]\n",
      "[ 0.19265087  0.07004414  0.03188846 ... -0.37784773 -0.40864572\n",
      "  0.02214655]\n"
     ]
    }
   ],
   "source": [
    "file0='u_10mps_y-5,7'\n",
    "path='/Users/cbrun/Documents/UFR/M1_TURB/TP/DATA_cylindre/tp2_oct2023/'\n",
    "filepath=path+file0+'.nc'\n",
    "r3 = nc.Dataset(filepath, 'r', format='NETCDF4')\n",
    "print(r3)\n",
    "print(r3['u'][:])\n",
    "print(r3['v'][:])\n",
    "print(r3['w'][:])\n",
    "u=r3['u'][:]\n",
    "v=r3['v'][:]\n",
    "w=r3['w'][:]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5f81f3f",
   "metadata": {},
   "source": [
    "# Spectres\n",
    "## statistiques utiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2480bece",
   "metadata": {},
   "outputs": [],
   "source": [
    "um=np.sqrt(np.mean(u)**2+np.mean(v)**2+np.mean(w)**2)\n",
    "up=u-np.mean(u)\n",
    "nu=15.e-6\n",
    "Fs = 1250\n",
    "nfft = 1024\n",
    "C=2\n",
    "t = np.linspace(0.00001,4096-1, 4096)/Fs\n",
    "Kx = 2*np.pi*Fs/um\n",
    "tau, autocorrelation = Autocorrelation.autocorr_lisse(u, 1024*4)\n",
    "rho=tau*um/Fs\n",
    "#Integral_cobra[i] = np.trapz(autocorrelation[0:np.where(autocorrelation<0)[0][0]], rho[0:np.where(autocorrelation<0)[0][0]])\n",
    "s2 = Fonction_structure.structure_f(up,2,1024*4)\n",
    "r = t*um\n",
    "epsilon_structure_courbe = (s2*r**(-2/3)/C)**1.5\n",
    "epsilon_structure_cobra = np.nanmax(epsilon_structure_courbe)\n",
    "Kolmogorov_cobra = (nu**3/epsilon_structure_cobra)**(1/4)\n",
    "k_Kolmogorov_cobra = 2*np.pi/Kolmogorov_cobra\n",
    "U_Kolmogorov_cobra = (epsilon_structure_cobra*nu)**(1/4)\n",
    "#Taylor_cobra = np.sqrt(15*nu*Varu_cobra/epsilon_structure_cobra)\n",
    "#k_Taylor_cobra = 2*np.pi/Taylor_cobra\n",
    "#k_integral_cobra = 2*np.pi/Integral_cobra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "39a5ee25",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f961f977dc0>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAApvUlEQVR4nO3de5hVdd3+8fcNBCGiFKByEEHl0cgDyogoRpKpoCZqmhDkU/lIqPzSTPOQdrAsMyzzlGEZeWXykKZimIfKR81DMZAHBDVCyREVREHAFMHP74/vJrbDDOyZ2XvW3nvu13Xta2avw56bdY1+5rvW96CIwMzMrNy0yzqAmZlZQ1ygzMysLLlAmZlZWXKBMjOzsuQCZWZmZalD1gGy0qNHj+jfv3/WMczM2rw5c+a8FhE9629vswWqf//+1NbWZh3DzKzNk7S4oe2+xWdmZmXJBcrMzMqSC5SZmZUlFygzMytLLlBmZlaWXKDMzKwslUWBkrSjpPslLZD0tKQzGjhGkq6UtFDSk5L2zds3StKzuX3ntW56MzMrhbIoUMA64KsR8RFgGHC6pEH1jhkNDMy9JgI/BZDUHrgmt38QMK6Bc83MrMKUxUDdiHgZeDn3/SpJC4A+wPy8w8YAN0ZawOoxSd0k9QL6AwsjYhGApOm5Y/PP3cSiZWs48WePNi/wU0+CBNtsC9tsA9t0hXbtm/dZZmbWoLIoUPkk9Qf2Af5ab1cf4MW893W5bQ1t37+Rz55Ian2xda9dmh/yg51hxQpY/nrug4EuW8O2uYK17TbQ6YPN/3wzMyuvAiVpa+BW4MyIeLP+7gZOic1s33RjxFRgKkBNTU3875cOaGbS3HlvvAGPPQaPPAKP3Ad//SusWZP29e4NBx4IBxyQvu6zD3Tq1MyfZ2ZWvWZManh72RQoSR8gFaebIuJ3DRxSB+yY974vsATo2Mj20vvQh2D06PQCWLcOnnoqFaxHH01fb7kl7evUCWpqUrHaULi2375VYpqZVSKlRzoZh5AE/Ap4PSLObOSYI4HJwBGkW3hXRsRQSR2A54BDgJeA2cBnI+Lpzf3MmpqaaJXJYl9+eWOxeuQRmDMH1q5N+3bZZWML68ADYY89oL2fZZlZ2yJpTkTUbLK9TArUQcBDwFPAe7nNFwD9ACLiulwRuxoYBbwFfCEianPnHwFcAbQHboiIS7b0M1utQNX3zjswd+7GgvXww/Dqq2nf1lvDsGEbW1jDhkG3bq2f0cysFZV1gcpCZgWqvgh44YWNBeuRR+DJJ+G991JPwUGD4MILYezYrJOamZVEYwWqbJ5BtVkSDBiQXuPHp22rVsHs2alY3XorTJiQnncdfni2Wc3MWpFbUOVu1Sr42Mfgn/+Ehx6CwYOzTmRmVlSNtaDKZSYJa0zXrjBrVnoWdeSR8OKLWzzFzKwauEBVgj59UpFatSoVqZUrs05kZlZyLlCVYq+90vOoBQvghBPg3XezTmRmVlIuUJXk0EPh+uvhvvvgS19KPQDNzKqUe/FVms9/Hp5/Hi6+OPX8u+iirBOZmZWEC1Ql+ta30tipb3wDdtoJTjop60RmZkXnAlWJpHSrr64OTj4Z+vaFT3wi61RmZkXlZ1CVqmPH1Glit93guOPg6c1OPWhmVnFcoCpZt25w112w1VZpRvUlrTOJu5lZa3CBqnT9+qUxUq+/DkcdBatXZ53IzKwoXKCqwT77wIwZaZLZE09M61KZmVU4F6hqccQRcO216Zbf5MkeI2VmFc+9+KrJxIlpjNSll6YxUueem3UiM7Nmc4GqNpdcksZInXdeej41blzWiczMmsUFqtq0awfTpsFLL6VZJ/r0gREjsk5lZtZkZfEMStINkpZKmtfI/nMkPZ57zZO0XtKHc/tekPRUbl8FLPDUCjp1gttvT7f5jjkG5s/POpGZWZOVRYECpgGjGtsZET+MiMERMRg4H3ggIl7PO2Rkbv8mC161WR/+MPzhD2lA7yc/mRY8NDOrIGVRoCLiQeD1LR6YjANuLmGc6jFgAPzxj7B2LRxyiBc7NLOKUhYFqlCStiK1tG7N2xzAvZLmSJqYTbIytscecM898MYbqUi98krWiczMClJRBQr4FPBwvdt7wyNiX2A0cLqkRnsESJooqVZS7bJly0qdtXwMGZLGR730UlpTavnyrBOZmW1RpRWosdS7vRcRS3JflwK3AUMbOzkipkZETUTU9OzZs6RBy87w4TBzJvzjH3D44V423szKXsUUKEnbAh8H7sjb1kVS1w3fA4cBDfYENNItvltugSeeSPP2rVmTdSIzs0aVRYGSdDPwKLCbpDpJJ0uaJGlS3mHHAvdGRP7/VbcH/iLpCeBvwKyIuLv1klego46C3/wGHnkkdUF/++2sE5mZNagsBupGxBanO4iIaaTu6PnbFgF7lyZVFTvhBHjrrTSQ9zOfSetKfeADWacyM3ufsmhBWQb++7/hmmvgzjvhc5+D9euzTmRm9j5l0YKyjJx2WmpJnXMOdO4Mv/hFmirJzKwMuEC1dWefnRY5/Pa3oUsXuOoqkLJOZWbmAmXAN7+ZevRNmZKK1KWXukiZWeZcoCwVo8suS0Xqssuga1e48MKsU5lZG+cCZYkEV1+ditRFF6WW1Fe+knUqM2vDXKBso3btUkeJf/8bzjorFamJnt7QzLLhAmXv16ED/PrXqXffpEmw1VYwYULWqcysDXKfYttUx47w29/CyJFpvNStt275HDOzInOBsoZ17gx33AHDhsG4cWnxQzOzVuQCZY3bemuYNQv23BOOOw7+7/+yTmRmbYgLlG1et25pwcNddkkTzT76aNaJzKyNcIGyLevRA+67D3r1gtGjYZ5XNDGz0nOBssL06gV/+lPq1Xf00fDaa1knMrMq5wJlhevXD26/HV5+GT79aVi7NutEZlbFXKCsaYYOhRtugAcfTLOhR2SdyMyqlAfqWtONGwfz58N3vwsf/ainRDKzknALyprn299OXc/PPttjpMysJMqiQEm6QdJSSQ12D5N0sKSVkh7Pvb6Rt2+UpGclLZR0XuulbuPatYMbb4S99oKxY1OLysysiMqiQAHTgFFbOOahiBice10MIKk9cA0wGhgEjJM0qKRJbaMuXWDmzDTrxNFHw/LlWScysypSFgUqIh4EXm/GqUOBhRGxKCLWAtOBMUUNZ5u3446pZ19dHRx/vHv2mVnRNKtASdpd0iGStq63fUutoJY4QNITkv4g6aO5bX2AF/OOqctta5CkiZJqJdUuW7ashFHbmGHD4Oc/T1Mhff/7WacxsyrR5AIl6cvAHcD/A+ZJym+xfK9YweqZC+wUEXsDVwG3b4jTwLGN9nuOiKkRURMRNT179ix+yrZswgQ48cRUoP7xj6zTmFkVaE4L6hRgSEQcAxwMXCTpjNy+hgpGi0XEmxGxOvf9XcAHJPUgtZh2zDu0L7CkFBmsAD/6EXTqBKef7vFRZtZizSlQ7fOKxQukIjVa0o8oUYGStIMk5b4fSsq9HJgNDJQ0QFJHYCwwsxQZrAC9e8Mll6R5+2bMyDqNmVW45hSoVyQN3vAmV6yOAnoAezYnhKSbgUeB3STVSTpZ0iRJk3KHHE+6nfgEcCUwNpJ1wGTgHmABMCMinm5OBiuSU0+FIUPgzDNh5cqs05hZBVM08VaMpL7Auoh4pYF9wyPi4WKFK6Wampqora3NOkZ1mjMnTYl02mlw1VVZpzGzMidpTkTU1N/e5BZURNQ1VJxy+yqiOFmJDRmSitM114D/CDCzZmrROChJn8j/avYf3/0ubL89fOlLsH591mnMrAK1dKDulHpfzZJtt4UrroC5c+Haa7NOY2YVqFgzSZSk955VuM98Bg47DL7+dVji3v9m1jRlMdWRVSkpPYd65x248MKs05hZhXGBstLadVc44wyYNg3+/ves05hZBXGBstK74ALo3h2++lXPMGFmBWtpgVqd+7qqpUGsinXrlhY4vP9+uPPOrNOYWYVo8kDdauGBuq1s3bq0uOG6dTBvHnTsmHUiMysTRRuoK2lXScMb2P4xSbs0N6BVuQ4dYMqUNNP5dddlncbMKkBzbvFdQcO39P6d22fWsNGj4dBD4Vvfgtebsz6lmbUlzSlQ/SPiyfobI6IW6N/iRFa9JLj88jSJ7He/m3UaMytzzSlQH9zMvs7NDWJtxJ57wsknw9VXe2FDM9us5hSo2ZJOqb9R0snAnJZHsqp38cVpYcNzz806iZmVsQ7NOOdM4DZJ49lYkGqAjsCxRcpl1WyHHeD889MUSA88AB//eNaJzKwMNbubuaSRwB65t09HxJ+LlqoVuJt5xv79b9htN+jZE2bPhnYeM27WVhWtm/kGEXF/RFyVe1VUcbIy0LkzXHppmu3817/OOo2ZlSH/2WrZGTs2rbx7wQWwZk3WacyszJRFgZJ0g6SlkuY1sn+8pCdzr0ck7Z237wVJT0l6XJLv2VWSdu3gRz+Cl15K3c/NzPI0p5MEAJI6AZ8mjX36z+dExMXN+LhpwNXAjY3sfx74eES8IWk0MBXYP2//yIh4rRk/17I2fDiccAL84AfwP/8DvXtnncjMykRLWlB3AGOAdcCavFeTRcSDQKNTC0TEIxHxRu7tY0Df5vwcK1OXXprm6POaUWaWp9ktKKBvRIwqWpLCnQz8Ie99APdKCuBnETG1sRMlTQQmAvTr16+kIa0Jdt45rRk1ZQpMnAjDhmWdyMzKQEtaUI9I2rNoSQqQ69p+MpA/wnN4ROwLjAZOlzSisfMjYmpE1ERETc+ePUuc1prkoougT59UoN59N+s0ZlYGmjOb+VOSngQOAuZKejbXeWHD9pKQtBfwc2BMRCzfsD0iluS+LgVuA4aWKoOVUNeuaXn4p55KLSkza/Oac4vvqKKn2AJJ/YDfAZ+LiOfytncB2kXEqtz3hwHN6aRh5eDoo+HTn06LGx5/PAwcmHUiM8tQk1tQEbE4IhYDp234Pn9bc0JIuhl4FNhNUp2kkyVNkjQpd8g3gO7AtfW6k28P/EXSE8DfgFkRcXdzMliZuPLKNE/fpEleHt6sjWvJVEdzc89+8rc9GRF7FSVZiXmqozJ23XVw6qkwbRr8939nncbMSqyYK+qeKukpUmvnybzX80DJnkFZGzJxYhofddZZ8JqHt5m1Vc3pxfcb4FPAzNzXDa8hETGhiNmsrWrXLrWiVqxIA3jNrE1qzjOolRHxQkSMq/cMymt4W/HssQdMmJAWNlyyJOs0ZpaB5tzi+0vu6ypJb+a9Vkl6s/gRrc365jfTDBNeHt6sTWpOC+qg3NeuEbFN3qtrRGxT/IjWZu28M5xyClx/PSxalHUaM2tlzZ5JQtKNkk6RtHsxA5m9z4UXQocOaWyUmbUpLZnqaBrQC7hK0j8l3SrpjOLEMsvp3RsmT06LGs6fn3UaM2tFLVlR98/AJcBFpCmIaoBTi5TLbKNzz4UuXeAb38g6iZm1opbc4vsT8DBwIvAssF9E+HafFV+PHmm281tvdSvKrA1pyS2+J4G1wB7AXsAekjoXJZVZfWecAVttBZddlnUSM2slLbnF95WIGAEcCywHfgmsKFIus/fr0SP16LvpJvjXv7JOY2atoCW3+CZL+l/gceAY4AbSmkxmpXHWWenr5Zdnm8PMWkVLVtTtDPwImBMR64qUx6xx/frB+PFpXNRFF6VWlZlVrZbc4vthRPzVxcla1bnnwr//DVddlXUSMyuxlnSSMGt9H/kIHHNMKlCrVmWdxsxKyAXKKs9558Ebb6RbfWZWtVygrPLsvz+MHJk6S7zzTtZpzKxECi5QkkZI2lnSryXNkDSiWCEk3SBpqaR5jeyXpCslLcwtjrhv3r5Rkp7N7TuvWJmszJ13XlqG49e/zjqJmZVIU1pQ44ALgbOACRR3WqNpwKjN7B8NDMy9JgI/BZDUHrgmt38QME7SoCLmsnJ16KGwzz5p4O769VmnMbMSaEqB+iiwfUQsjYi1wMpihYiIB4HNLXg4BrgxkseAbpJ6AUOBhRGxKJdpeu5Yq3YSnH8+PPcc3H571mnMrASaUqAuAvLX376nyFk2pw/wYt77uty2xrY3SNJESbWSapctW1aSoNaKjjsOdt0Vvv99iMg6jZkVWcEFKiIeyLV0Nti+BHkaowa2xWa2NygipkZETUTU9OzZs2jhLCPt28PXvgZz5sC992adxsyKrCW9+B4vVogC1AE75r3vCyzZzHZrK046Ka28e+aZsHZt1mnMrIhaMpPEY8UMsgUzgZNyvfmGASsj4mVgNjBQ0gBJHYGxuWOtrejUKQ3afeYZ+NGPsk5jZkXUkslizy1WCEk3A48Cu0mqk3SypEmSJuUOuQtYBCwErgdOA8hNszSZ9DxsATAjIp4uVi6rEEccAcceC9/5DixenHUaMysSRYEPlyXNyH8LDI6IgSVJ1QpqamqitrY26xhWLP/6V5oG6bDD4Lbbsk5jZk0gaU5E1NTf3pQW1JsR8Znc6wTgj8WLZ9ZC/fqlJeFvvx3uvjvrNGZWBE1pQQ2IiOfz3n84IjY3dqmsuQVVhdauhUGDoHNnePzx1MvPzMpei1tQG4qTpHaSegM7STpQ0kBJXYuY1ax5OnZMY6LmzYNf/SrrNGbWQgW1oCSJ1ENuDNCLNOvDSuBtoBvwodyhTwA3RcQTpQhbTG5BVakIOOAAePHFNMtEly5ZJzKzLWisBbXFFXUljQQOBx4CvhgRb23m2J2BoyV9Frg8Ipa2ILNZ00kwZQp87GPw4x/DhRdmncjMmmmzLShJnySNOZrdpA9NLa7PAX+OiLqWRSwNt6Cq3LHHwh//mMZH9Wl09iszKwPNegYVEX8spDhJet900rlJXW8s1+JkbcCUKfDuu3D22VknMbNmKtaChQ3NiWeWnV12SWtGTZ8Of/pT1mnMrBmKVaA8lbSVn3PPTfP0TZ7sefrMKtBmC5Skk3KvyZJOaa1QZkXRufPGefp+/OOs05hZE22pBbU9sHtEXA30boU8ZsV1xBFwzDFw8cVpOiQzqxhbKlACfi+pHTC3FfKYFd8VV6TxUV/5StZJzKwJttSL7zJgBfB5oLekfVshk1lx7bRTGg/1u995nj6zClLwXHyb/RBpfURU1MRnHgfVxrzzDuy1F6xZA3/+M/zXf2WdyMxymjUOStJ4SWdL2kbSqNLFMyuxTp3gt79NvflGjEjz9ZlZWdvSM6idgZ8CZwEuUFbZ9toLHnggzXJ+8MGpd5+Zla0tFag5EbEG+DbwWivkMSutj3xkY5EaPRpeeSXrRGbWiC11krgr9+2JEfHdVshjVnq77gqzZsGyZakb+urVWScyswYUOpPEYSVNAUgaJelZSQslndfA/nMkPZ57zZO0XtKHc/tekPRUbp97PtiW1dTAjBnwxBPwhS+kbuhmVlYKLVAlnWtPUnvgGmA0MAgYJ2lQ/jER8cOIGBwRg4HzgQfqreg7Mrd/k54gZg064gj4wQ/gllvSQodmVlYKLVC3lzIEMBRYGBGLImItMJ20OGJjxgE3lziTtQVf/SqMH5/GSc2YkXUaM8tTUIGKiDsk9ZS0i6QPlCBHH+DFvPd1uW2bkLQVqUfhrfkRgXslzZE0sbEfImmipFpJtcuWLStCbKt4Elx/PRx0EEyYAPfdl3UiM8vZYoHK/U/9PuAS4EvANEnXShpQxBwN3UJs7KHAp4CH693eGx4R+5JuEZ4uaURDJ0bE1IioiYianj17tiyxVY/OnWHmTBg0CI4+Gr73vTSw18wyVUgLqha4LyImRsTXImI8cA4FLBffBHXAjnnv+wJLGjl2LPVu70XEktzXpcBtpFuGZoXr1g3uvReOPBK+/nUYMgT+8Y+sU5m1aVssUBExFxiTW3Jjt9y2NRGxyX+9kr4vqUPu+3aStikwx2xgoKQBkjqSitDMBj5/W+DjwB1527pI6rrhe1KPQ08TYE233Xapw8SsWWl81H77pWmRzCwThXaSOA5YCBwr6eebOW7riFgHEBHvAQUtwpM7ZzJwD7AAmBERT0uaJGlS3qHHAvfmBg9vsD3wF0lPAH8DZkWEZwS15jviCJgzB/r0gbFj4dVXs05k1iYVdbJYST+JiDPytl8REWe2+AeUgCeLtS2aPz/d6jvkELjzztShwsyKrlmTxeZO7JL7unUBP+dhSVMk9ZXUC9iu6VHNysSgQWmc1KxZMHVq1mnM2pxCbvF9SNJk4KAtHRgRM4BbgAuB7wGeHskq2+TJcOihcNZZ8NxzWacxa1MKKVCHkBYs3FnSFltEEfFYREyKiC9ExPyWBjTLVLt28MtfpuU6xo2Dt97KOpFZm1FIgfob8EVgca4bt1nb0qcP3Hgj/P3vaTDv+vVZJzJrEwrpZr4gIp6MiFmtEcisLB11FFxxBdx2G5x7btZpzNqEYg62NatuX/4yLFwIl18Ou+wCp56adSKzqlZIL75DJV0vaXDufaNz3ZlVvR//OM02ccYZaakOMyuZQp5BnUaa2miCpE8Ag0uayKyctW8Pv/oV9OiRnke9/XbWicyqViEFallErIiIs0nTCO1X4kxm5a17d/jFL2DevHS7z8xKopAC9Z/OERFxHnBj6eKYVYjRo+G449LM53V1Wacxq0qFdjP/j4i4qkRZzCrLlCmpy/mnPw0vv5x1GrOqU0iBugRA0nhJD0s6ssSZzCrDgAFw883w9NOw//6weHHWicyqSiEFakXu62Gk6Y6OLVkas0pz7LHw0EPw5ptpSqQ33sg6kVnVKKRAdZB0IfBipKnP12zpBLM2ZZ994Pe/h0WL4Mwzs05jVjUKKVBfB/4CfCf3vj2ApA+WKpRZxTnoILjggjQl0l13ZZ3GrCpscT0oSVcA7wHzSbf7BuReN0XEw7lj1kdE+5ImLTKvB2VFt3Yt7L03rFuXuqB36pR1IrOK0Oz1oHILDn4H+CcQwN0RcdqG4rTh84sV1KxideyY5utbuDDN17duXdaJzCraZguUpH0BIuKNiLg/Im6NiKfqHxcRm3yOpN6FLM+Rd/woSc9KWijpvAb2HyxppaTHc69vFHquWas5/HCYOBF+8hMYMSIVKzNrli21oN6QdK6kjxT6gZK2lnQKsHehy3NIag9cA4wGBgHjJA1q4NCHImJw7nVxE881ax3XXQc33QQLFqQl4x9+eMvnmNkmNjubeUQ8L+ly4DRJFwAvAI8BL5OeR70DbAt0B/YGhgDLgSkR8UoTcgwFFkbEIgBJ04ExpOdepTzXrPgk+OxnYfhwOOyw1Kp69FHYc8+sk5lVlC0utxER64ArASTtCRwOfArYDtgKWEYqWH8BboiI5iw52gd4Me99HbB/A8cdIOkJYAlwdkQ83YRzzVrXTjvB/fdDTQ0cfTTMmAH7eSpLs0I1aT2o3POnTZ5BFUFDnSzqdy+cC+wUEaslHQHcDgws8Nz0Q9JSIRMB+vXr1+ywZgXr3RvuuAOOOQaGDoVPfSpNkfRf/5V1MrOyV8g4qNZQB+yY974vqZX0HxHxZkSszn1/F/ABST0KOTfvM6ZGRE1E1PTs2bOY+c0at99+MH8+fOtb8OCDqSv6X/+adSqzslcuBWo2MFDSAEkdgbHAzPwDJO0gSbnvh5KyLy/kXLPMbbstfPObqeNEr14wZkyax++997JOZla2yqJA5Z5zTQbuARYAMyLiaUmTJE3KHXY8MC/3DOpKYGwkDZ7b+v8KswL06gV33gnbbZc6Uuy9N9x2W9apzMrSFmeSqFaeScIy9d57MH16Wk9qwQK45x745CezTmWWiWbPJGFmJdCuXWpBPfoofOQjqRPFLbf4lp9ZHhcosyx17Qr33Qe77gonnAD9+6d5/MzMBcosc716QW1t6jTx7rtpjaklDXZENWtTXKDMykGHDjB2LNx6K7zyChxwAFx1FaxcmXUys8y4QJmVkwMPhD/9CdasgS9/GUaOhKUFTWlpVnVcoMzKzdCh8OKL8NvfwjPPpDn8xoyBuXOzTmbWqlygzMpR585w/PFp5omaGnjggXTb7/LL4fXXs05n1ipcoMzKWU0NzJqV1pU67DA4+2zYeedUuMyqnAuUWSXo0QNmzkxd0nv1SrOj+5afVTkXKLNKIaXZJu6+Gz74wbQY4mmnwfr1WSczKwkXKLNKs9NO8NhjcPrp8NOfwqmnwttvZ53KrOhcoMwqUf/+cPXVcM45cP318NGPwnPPZZ3KrKhcoMwq2Q9+AL//PaxalXr5nXUW/P3vWacyKwoXKLNKJsGRR6bBvd26pdknhg1L3dFXr846nVmLuECZVYM994R//hNefhkGDUrd0Q85BN58M+tkZs3mAmVWTXr0SN3PZ8xIE9B26wbjx8PatVknM2syFyizaiOlpTv+7/9g8mT4zW/SzBSHHAI33ZTm+TOrAC5QZtXqYx+DK69My3icdFJqWU2YAPvtl1pYb72VdUKzzSqbAiVplKRnJS2UdF4D+8dLejL3ekTS3nn7XpD0lKTHJXkdd7N8Y8fCL38Jy5alHn91dXDiiWkC2mXLsk5n1qiyKFCS2gPXAKOBQcA4SYPqHfY88PGI2Av4DjC13v6RETG4oXXtzYy05tSRR8LixXDJJfDnP8OAATBuXJpCyazMlEWBAoYCCyNiUUSsBaYDY/IPiIhHIuKN3NvHgL6tnNGsOnzoQ3DBBfD003DQQTB9epqI9uCDU9GKyDqhGVA+BaoP8GLe+7rctsacDPwh730A90qaI2liYydJmiipVlLtMt/asLZu993TvH4rVsApp8ATT6SOFBMm+PmUlYVyKVBqYFuDf8ZJGkkqUOfmbR4eEfuSbhGeLmlEQ+dGxNSIqImImp49e7Y0s1l12HZbmDo1PZv62tdSr7/ttoPhw9N7r+hrGSmXAlUH7Jj3vi+wpP5BkvYCfg6MiYjlG7ZHxJLc16XAbaRbhmbWFF26bJw6aeRIeOGFNCPFLrvA17+eCphZKyqXAjUbGChpgKSOwFhgZv4BkvoBvwM+FxHP5W3vIqnrhu+Bw4B5rZbcrNoceSTceSe89FJ6TjV8OHz/+7DjjvCJT8CXvpT2m5VYWRSoiFgHTAbuARYAMyLiaUmTJE3KHfYNoDtwbb3u5NsDf5H0BPA3YFZE3N3K/wSz6rThOdWTT6blPR5+ON0OPPpo+Oxn03azElG00R47NTU1UVvrIVNmTbJ+PSxZAj/5Sbr9B9A316F2zBi47DLYaqvs8llFkjSnoSFCLlBm1jwLF8KvfpXWoXr+eZg9G9q1g898BnbbDb74RejXL+uUVgFcoOpxgTIrsvvvh+99D555JnWo6NQpzaq+776w//7QZ3MjR6wta6xAdcgijJlVoZEj0wvSs6mzzkoFKyLNYnHOOTBkSBocvP322Wa1iuACZWbFt9de8Mc/wmuvpRV+zz479QTc4LOfhVGj0lL1e+4JH/hAdlmtbPkWn5mV3rvvwr33pg4WM2fCH/6QOlxssPvuaTaLE0+Erl3TSw2N37dq5GdQ9bhAmWVozZrUueLhh2HWrLRk/bvvvv+YIUPStEsTJ7pnYJVzgarHBcqsjLz7LtxzDzzyCMybl77PXwV4/Hg47ri0nP2uu6ZnWlY1XKDqcYEyK3MR8OCDMGVKmn4p3+GHpymYdtgBRoyAAw/0c6wK5gJVjwuUWQVZvjy1rh59FH72M3j7bXjvvfQV0vOqT34Shg2DoUPT+Kvddktd3a3suUDV4wJlVgWeeSatFvzMM6lr+wsvbNzXvj0MHAhHHZXGYO26a2ppffjDmcW1hnkclJlVn913TzOwb/Cvf8HvfpeK1fLlqcfgM8+8/5yOHdN5n/hEevbVu3fqiNGjR+tmty1yC6qMHHzwwZx//vkcfvjh/9l2xRVX8Nxzz3HttddmmMysQkXA66/D3/4G//hHWoixthZuvXXTY7t3T/MKDhkCe+yRWlp9+8LgwWmflYxbUBVg3LhxTJ8+/X0Favr06fzwhz/MMJVZBZNScRk9Or02iICnnkrPsWpr4dVX4fHHYfFiuOGGTT+nY8c0ZVNNDWy9dZpz8JBDUm/C4cPT7UQrOregysjy5cvZfffdqauro1OnTrzwwguMGDGCxYsXIw9aNGsd772XCtW//pUmwX3++dSbsK4Oli2DlSvff3yHDmlV4k6d0jOv7t03dtTo3Rv22w86d/bA481wC6qpzjwz/UVVTIMHwxVXNLq7e/fuDB06lLvvvpsxY8Ywffp0TjzxRBcns9bUrh0MGJBeH//4pvvffjs945o9OxUyKb1/6KE08HjduvQcrL4uXdJnvv566rQxenSak3CHHdL+vfeGXr08KDmPC1SZ2XCbb0OBuqGh2w1mlp0PfjC1kIYObXj/u+/Cyy+nAceLFsFjj6XitGABrFiRCtzs2enV2Od3755eO+2UnoNtt10qnN27p6LWrl3a161bKnBdupTqX5sp3+IrM6tXr2bnnXfm7rvvZty4cTz77LNZRzKzYouAd95JPQxXr07juzp3TlNAvfZa6oH46qvp9uLSpen95myYv7BTp3S7sX//VCj790/bVqxILbQuXVKB69kzFbcNr4zHi/kWX4XYeuutOfjgg/niF7/IuHHjso5jZqUgpZbS4MHp/UEHbf74detSkVm3Lj0be+65NBVU+/bw0kupiK1YkQrZa6/B/PmpFXfXXakYFpJnq61Sgdxll1Tkttkmvd9uu9RKa9cuHdOjRyqG3bqlotqtGxx6aEmesZVNgZI0CvgJ0B74eURcWm+/cvuPAN4CPh8Rcws5t9KMGzeO4447junTp2cdxczKQYcOG8dp7bBD47cXGxKRCtjSpamDR/v2qZi98Ub6umJFasW9+mr6funSVDxXrUrfz52bfnZE6iTSUMF76KEtF9lmKIsCJak9cA1wKFAHzJY0MyLm5x02GhiYe+0P/BTYv8BzK8qxxx5LW731amZFJqXbesVYJDIitdxWr07FbuVKeOWVkhQnKJMCBQwFFkbEIgBJ04ExQH6RGQPcGOn/3I9J6iapF9C/gHPNzKylpPS8qlOnVhm83K7kP6EwfYAX897X5bYVckwh5wIgaaKkWkm1y5Yta3FoMzMrnXIpUA09Xat/j6uxYwo5N22MmBoRNRFR07NnzyZGNDOz1lQut/jqgB3z3vcFlhR4TMcCzjUzswpTLi2o2cBASQMkdQTGAjPrHTMTOEnJMGBlRLxc4LlmZlZhyqIFFRHrJE0G7iF1Fb8hIp6WNCm3/zrgLlIX84WkbuZf2Ny5GfwzzMysiDyThJmZZaqxmSTK5RafmZnZ+7hAmZlZWWqzt/gkLQMW523aFljZwKGFbq//vgfwWgtjFqqxjKU4v5BjN3dMU65zIdva6nVu7v5CrqmvcWHH+He5eMcOjIhtN9kaEX6lIj21JdsbeF+bdfZSnF/IsZs7pinXuZBtbfU6N3d/gdfU17iAY/y7XPrr7Ft8G93Zwu2NHdcaWvqzm3J+Icdu7pimXOdCt7WWcrrOzd1fyDX1NS7sGP8uF+/YBo9ps7f4Sk1SbTTQK8WKy9e59HyNW4ev86bcgiqdqVkHaCN8nUvP17h1+DrX4xaUmZmVJbegzMysLLlAmZlZWXKBMjOzsuQCZWZmZckFKgOSjpF0vaQ7JB2WdZ5qJWlnSb+QdEvWWaqJpC6SfpX7HR6fdZ5q5d9fF6gmk3SDpKWS5tXbPkrSs5IWSjpvc58REbdHxCnA54ETSxi3YhXpOi+KiJNLm7Q6NPF6HwfckvsdPrrVw1awplxn//66QDXHNGBU/gZJ7YFrgNHAIGCcpEGS9pT0+3qv7fJOvTB3nm1qGsW7zrZl0yjwepNWrX4xd9j6VsxYDaZR+HVu88piwcJKEhEPSupfb/NQYGFELAKQNB0YExHfB46q/xmSBFwK/CEi5pY4ckUqxnW2wjXlegN1pCL1OP4jt0maeJ3nt3K8suNfruLow8a/KCH9B9xnM8f/P+CTwPEbVg22gjTpOkvqLuk6YB9J55c6XBVq7Hr/Dvi0pJ+S7Xxy1aLB6+zfX7egikUNbGt0io6IuBK4snRxqlZTr/NywH8ANF+D1zsi1gBfaO0wVayx69zmf3/dgiqOOmDHvPd9gSUZZalmvs6ty9e7dfg6N8IFqjhmAwMlDZDUERgLzMw4UzXydW5dvt6tw9e5ES5QTSTpZuBRYDdJdZJOjoh1wGTgHmABMCMins4yZ6XzdW5dvt6tw9e5aTybuZmZlSW3oMzMrCy5QJmZWVlygTIzs7LkAmVmZmXJBcrMzMqSC5SZmZUlFygzMytLLlBmFU6J/1u2quNfarMKJKm/pAWSrgXm8v653MyqgmeSMKtAuTWFFgEHRsRjGccxKwm3oMwq12IXJ6tmLlBmlWtN1gHMSskFyszMypILlJmZlSV3kjAzs7LkFpSZmZUlFygzMytLLlBmZlaWXKDMzKwsuUCZmVlZcoEyM7Oy5AJlZmZl6f8Ds//DV0IRnhsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#ex de jeu de donnéees statistiques:\n",
    "eta = 1.\n",
    "plt.semilogx(r[1:]/eta, (epsilon_structure_cobra*r[1:])**(-2/3)*s2[1:], c='red', label='V')\n",
    "plt.ylabel(r'$\\left(\\frac{(S2*r)^{-2/3}}{C}\\right)^{1.5}$ with C=2')\n",
    "plt.xlabel('r')\n",
    "plt.axhline(C)\n",
    "#plt.axhline(C*4/3)\n",
    "plt.tight_layout()\n",
    "plt.legend(edgecolor = 'white', framealpha =1, loc=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ac2b944b",
   "metadata": {},
   "outputs": [],
   "source": [
    "nfft = 1024\n",
    "um=np.mean(u)\n",
    "up=u-um\n",
    "Kx = 2*np.pi*Fs/um\n",
    "f_1 = np.zeros((1), dtype=np.ndarray)\n",
    "SU_1 = np.zeros((1), dtype=np.ndarray)\n",
    "f_1, SU_1 = signal.welch(up, Fs, nperseg =nfft, return_onesided = True, detrend=False, window = 'hanning')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "86438782",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f961ff230d0>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEBCAYAAAB7Wx7VAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAweklEQVR4nO3dd3gUVdvH8e+dCoQWiiEQOkhRxBSqCChVQFCagC+PiqKg2ECUR4THAoIdQaQoioA0UUQQFAHpAdKoIgSpoQZCCUlIPe8fwRBCAskmZGeT+3Nde5mdnZm9lzG/PTlz5owYY1BKKVXwOdm7AKWUUvlDA18ppQoJDXyllCokNPCVUqqQ0MBXSqlCQgNfKaUKCQ18pZQqJDTwlVKqkMi3wBeRGiIyQ0QW5dd7KqWUuiZbgS8i34jIGRHZnWF5RxHZJyIHRGTEzfZhjDlojHk6N8UqpZSynUs215sJfAHM+neBiDgDk4F2QAQQJCK/AM7AuAzbDzDGnMl1tUoppWyWrcA3xqwXkWoZFjcGDhhjDgKIyHygmzFmHNDF1oJE5FngWQAPDw//unXr2rorpZQqlEJCQs4aY8pnXJ7dFn5mKgHH0j2PAJpktbKIlAXGAr4i8t+rXww3MMZMB6YDBAQEmODg4FyUqJRShY+IHMlseW4CXzJZluXUm8aYc8CgXLyfUkqpXMjNKJ0IoHK65z7AidyVk0pEHhaR6RcvXsyL3SmllCJ3gR8E1BaR6iLiBvQBfsmLoowxS40xz5YqVSovdqeUUorsD8ucBwQCdUQkQkSeNsYkAUOA34G9wEJjzJ68KEpb+EoplffEyne80pO2SimVcyISYowJyLhcp1ZQSqlCwpKBr106jssYw6xZs0hMTLR3KUqpDCwZ+HrS1nH9+eefPPHEEzz44IOcOJEng7aUUnnEkoGvHNeDDz7I999/T2hoKL6+vvz555/2LkkpdZUlA1+7dBxbv3792LZtG56enrRt25Zx48aRkpJi77KUKvQsGfjapeP47rrrLoKCgujVqxdvvvkm3bp14/z58/YuS6lCzZKBrwqGEiVKMG/ePCZNmsTvv/+On58fISEh9i5LqUJLA1/dViLCkCFDWL9+PcnJyTRv3pzp06dj5es/lCqoLBn42odf8DRt2pTQ0FBat27Nc889xxNPPEFsbKy9y1KqULFk4GsffsFUrlw5li9fzv/+9z/mzJlDkyZN2L9/v73LUqrQsGTgq4LL2dmZt99+mxUrVnDy5EkCAgJYtEhvc6xUftDAV3bRoUMHwsLCqF+/Pr169eLVV1/Vq3OVus0sGfjah184VK5cmfXr1/Piiy8yYcIEWrduTUREhL3LUqrAsmTgax9+4eHm5sbEiROZN28eO3bswM/Pj1WrVtm7LKUKJEsGvip8+vTpQ1BQEOXKlaN9+/aMGTNGr85VKo9p4CvLqFevHtu2baNv376MGjWKLl26cO7cOXuXpVSBoYGvLKV48eLMmTOHL7/8ktWrV+Pn50dQUJC9y1KqQNDAV5YjIgwePJiNGzcC0KJFC6ZMmaJX5yqVS5YMfB2lowAaNWpEaGgobdq04fnnn6d///7ExMTYuyylHJYlA19H6ah/lS1blmXLlvHee+8xd+5cGjduzN9//23vspRySJYMfKXSc3Jy4q233mLlypVERkYSEBDAggUL7F2WUg5HA185jLZt2xIaGkrDhg3p06cPL730EgkJCfYuSymHoYGvHIqPjw9r167llVdeYdKkSbRs2ZJjx47ZuyylHIIGvnI4rq6ufPbZZyxcuJA9e/bg6+vLypUr7V2WUpanga8cVq9evQgODsbb25uOHTvyzjvvkJycbO+ylLIsDXzl0OrUqcOWLVv4v//7P95++206derE2bNn7V2WUpakga8cnoeHB9999x3Tpk1j7dq1+Pn5sXXrVnuXpZTlWDLw9cIrlVMiwrPPPsvmzZtxdnbm/vvvZ9KkSXp1rlLpWDLw9cIrZSt/f39CQkLo0KEDL730Ev369SM6OtreZSllCZYMfKVyo0yZMixZsoT333+fhQsX0rhxY/bs2WPvspSyOw18VSA5OTnx3//+l1WrVhEVFUXjxo2ZO3euvctSyq408FWB9sADDxAWFoafnx+PP/44L7zwAvHx8fYuSym70MBXBV7FihVZs2YNw4YN48svv+T+++/nyJEj9i5LqXynga8KBVdXVz7++GN+/PFH9u3bh5+fHytWrLB3WUrlKw18Vah0796d4OBgfHx86NSpE6NGjdKrc1WhoYGvCp3atWuzZcsWnnrqKcaMGUPHjh2JjIy0d1lK3XYa+KpQKlq0KN988w1ff/01GzZswNfXl82bN9u7LKVuKw18Vag9/fTTBAYG4u7uTqtWrZgwYYJenasKrHwNfBF5RES+EpElItI+P99bqaz4+voSEhJC586defXVV+nduzeXLl2yd1lK5blsB76IfCMiZ0Rkd4blHUVkn4gcEJERN9uHMeZnY8xA4EngMZsqVuo2KF26NIsXL+aDDz5g8eLFNGrUiF27dtm7LKXyVE5a+DOBjukXiIgzMBl4CKgP9BWR+iLSQESWZXjckW7Tt65up5RliAivv/46q1ev5tKlSzRp0oTZs2fbuyyl8ky2A98Ysx6IyrC4MXDAGHPQGJMAzAe6GWN2GWO6ZHickVQfACuMMaGZvY+IPCsiwSISrCMnlD20atWK0NBQGjVqxH/+8x8GDRrElStX7F2WUrmW2z78SkD6G4pGXF2WlReBtkBPERmU2QrGmOnGmABjTED58uVzWZ5StvH29mb16tW88cYbTJs2jRYtWnDo0CF7l6VUruQ28CWTZVkOcTDGTDTG+BtjBhljpubyvZW6rVxcXBg/fjxLlizhwIED+Pn5sWzZMnuXpZTNchv4EUDldM99gBO53KfeAEVZSteuXQkJCaFatWo8/PDDvPnmmyQlJdm7LKVyLLeBHwTUFpHqIuIG9AF+yW1RegMUZTU1a9Zk8+bNPPPMM4wbN4727dtz+vRpe5elVI7kZFjmPCAQqCMiESLytDEmCRgC/A7sBRYaY3J9pwlt4SsrKlq0KF999RXffvstgYGB+Pr6snHjRnuXpVS2iZWvKgwICDDBwcH2LkOpG+zYsYOePXty6NAhPvjgA4YOHYpIZqe0lMp/IhJijAnIuFynVlDKBg0bNiQ4OJhu3brx2muv0aNHD/QvUmV1lgx87dJRjqBUqVIsWrSITz75hF9++YWAgAB27Nhh77KUypIlA19P2hZsmw+c5cV5YRw4E23vUnJNRBg6dChr164lJiaGpk2bMnPmTHuXpVSmXOxdwM2cvHiF95b9ZdO2Ajg7Cy5OgouTEy5OgrOz4OrkhLOT4OKcbnm6585OV7dJ/zzdfpydBFfnq9s4OaW9lv65q7MTbi6W/C5NczEukWNRsXi4u1D86qOIq9Nt7YfeGXGBj37fx4bwswAci4rlx8HNcXZy/L7vFi1aEBYWRt++fXnqqafYtGkTEydOpGjRovYuTak0lgx8EXkYeNi9Qi0WBB275fqZSTGGpBRD8tVHfivj4UZlz6L4lClGZc9i+HgWpXKZYlT2LEolz6K4uzjne02XriSy6q/T/LrzJOvDI0lMvv7fxdlJ8HBzpkQRVzzcnfEs5oZfVU+a1ihLQFVPPNxt+98l/HQ0n63az/Jdpyjj4cZbnetRoogLb/y4i7lbj9C/WbU8+HT25+XlxR9//MHo0aN5//33CQ4OZtGiRdSsWdPepSkFFJJROsakhn7S1UdysiExJeXasuSUtC+HpGRDUsq154nJ6dczJF/32rXnqdulPo9PTOHExStEnI8l4nwcx8/HkZCccl1NXiXdqexZLO1LwMezGD5lilLZsxjepYrg4pw3fyFcjk9i1V+nWbbzJOv3R5KQnELFUkXo1MAb/6qexCUmExOfRHR8EjHxSVy+ksTl+GQuxydy+lI8u49fJCnF4Owk3ONTiibVy9K0RhkCqpWh+E2+AE5ejGPZjpMs3XmCnREX8XBz5pn7a/DM/dUpUcQVYwz9Z2xjx7ELrBrWCq+SRfLk81rFsmXL6N+/P8YYvvvuO7p162bvklQhktUonUIR+PaWkmI4HX2FY1FxRJyP5VhUHMfOx3IsKvUL4eTFONL/EeLsJFQsXQSf0sXwKulOqaKulCrqSsmr/017FLv2c1FX57TumJj4JFb/fYZfd57gz32RJCSlUKFkash3vscb38qlccpmN0psQhIhR86z5eA5thyMYsexCyRdLbashxvepYvgXaoo3qVS/+vqLKz86zRBh6MwBhpUKkXXhhXp7leJssXdr9v34bMxtJ+wnnb1vJj8uF+26klMTmH13tPU9ipBzfLFs7WNvRw6dIiePXsSGhrK66+/ztixY3FxseQf1aqA0cC3sMTkFE5euMKx87E3fCGcvZzAxbhELl1J5GaHytVZ0r4Ujp+PIz4phTtKuNOpgTdd7vHGr4pntkP+Zv79Ath+9AInLl7h5MU4Tl28wokLcVy6kjrdQM3yHnRtWImHG3pT4xah/MWacD5euZ9vngzgwbpeWa53OT6J+duO8s3GQ5y4eAWvku4sHdKCO27jXwY7Iy6w//Rlevr72LyPK1eu8MorrzBt2jRatWrF/PnzqVChQh5WqdSNHCrw/+3Dr1Wr1sDw8HB7l2MJKSmG6PgkLsUlcjEukQuxqf/N+LgUl0i54m50auBNQLUy+XpCNCY+iUtXEqlQski2T/4mJKXQeeIGYhOS+WNoS4q53dgC3hh+lue/D+HSlSSaVC9Dt3sr8d6yv6jnXYJ5zzbF3cWZyOh4Ri7eReTlePyreNKoehna1fPK9Esu5EgU+09fpl19L8pl+Kvj2jrn6T9jK7EJybzT9S6eaF4NgCuJySzZfpw1f5+hw10V6O6XvS+D2bNn89xzz1GqVCnmz59Pq1atsrWdUrZwqMD/V2Fp4Rd2QYej6DU1kE4NKjDu0XsoVcw17bW1+87w7OwQapTz4IMe99CwcmkAlu86yfPfh9KnUWX6NK7CoNkhnI9N4O5Kpdh1/CIJSSkMbXcnL7Wpfd17xSYk0fLDtZy9HI+TQPv6FZjY1/e6UVV7T16i99RAyhZ3o1o5D9bvj+Sjng1pW8+Lp78LIvjIeURSu6t+GdIi259z9+7d9OjRgwMHDvD+++/z+uuv69W56rbQK22VZTWqVobhHerw+57TtPl0LUu2H8cYw5q/T/PsrBBqlS/OvIFN08IeoFMDb154oCbzg47RY8pmXJyFn55vzo+Dm7Pr7fa0rlOemZsPcyUx+br3+nbTYc5ejufjXg15ukV1fttziilr/0l7PS4hmSFzQynm7sz3A5vy5eN+NPApzbAfdtDo/VVsP3aBiX19GdyqJn+duERcwvX7v5m7776boKAgunfvzogRI3jkkUe4cOFCbv/5lMo2DXxlCS88UItfhtxHpdJFeXn+dnpPC+S52SHUqVCCuQOb4OnhdsM2w9rV4VHfSrStdwdLh7TgroqpF+q5uzgzqFVNomIS+Cn0eNr6F2MTmbbuHx6sewc9/X0Y2bk+XRtW5Is/wwk/nXoR2LgVe/knMoZPet1LpdJFKebmwo+DmvFFP19a1CrHV/8JoGvDivhX9SQpxbAz4kKOPmfJkiVZuHAhEyZMYPny5fj7+xMWFmb7PxzYZdixckwa+Moy7qpYip+ev4+3H67PXycuUb9iKeY804TSxW4MewAnJ+Gzx+5lWv+AG74QmlQvQ4NKpfh640FSrgbitPX/cOlKEq+1r5O23v8erk9xdxeemxNCn+mBzAo8wtMtqtOidrm0dVycnehyT0W+ebIRD9RNvTWzbxVPAEKOns/x5xQRXn75ZdatW0d8fDzNmjVjxowZ2NK9umT7cWq+uZxTF/UWjOrWLBn4OpdO4eXsJDx5X3W2vNmGRYOaUaqo6603yoSI8Mz91TkYGcMfe08zff0/fL3hEF0bVqR+xZJp65Ut7s6YRxpwNjqeuMQUBt5fneEd6txkz6nKeLhRo5wHoUcuAKkjrcJyGP7NmzcnLCyM+++/n2eeeYYBAwYQGxubo33M3HwYIMd/aajCyZKBr3PpqBJFXHHN5cVnnRp4U7FUEV74PpT3l/9NqzrlebvrXTes1/keb3a+3YElL9zHyM71KeKavaug/ap6Enr0PMYYJqzaz6NfbmbrwXPZ2jY5xTB6yW5OJ7jx22+/MWrUKGbOnEmzZs242ci0A2eiuRCbkPY86erV0uFnLmfrfVXhZsnAVyovuDo78UrbO/EqWYRJfX2Z3t+fMpmcC7CVXxVPomISmLvtKDM2pt7gfNKaAzesFxWTei1FeuvDI5kVeIQfgo/h7OzMu+++y/Lly4mIiCAgIICffvrphv38E3mZDhM2EDBmFYfOxmCM4ci5GAD+PuX4E9Gp208DXxVovRtVZtOIB3m4YcU8HwLZ+R5v7qpYkpGLd5OcYujftCobD5wl5Mh5EpJSCPznHMkphl5TN9N3+haS0k2vsWBb6hxRO49f67Z86KGHCA0NpU6dOvTo0YPXXnuNxMTUL4qk5BRGL9mdNs3HxgNnORoVm3ax298nL+XpZ1MFk17nrZSNShV1ZeFzzRi1ZDd1K5Tg8SZVWbH7JMN/2EE975L8uuskbet58U9kaiv8tR92UL6EOwNaVGfV3tO4uzjx14lLJCanpHVfVa1alQ0bNjBs2DA++eQTtm7dyoIFC5i/O5pNB84xrnsD3v5lD0fPxbDt6rUDreuUZ8vBcxhjdFy/uilt4SuVCx7uLnza+16ebVkTD3cXJvfz42hULL/uOolXSXdW7T1NpdJFaV2nPD9vP8FXGw4xeE4oSSmGQa1qEp+Uwv7T0VxJTCYlxXD4bAzJ4swXX3zB1BkzCQkNxdfXl29/XEbrOuXp27gKVcsW4/C5WIIOR1G6mCut7yzPlcQUIi/H2/ufQ1mcJVv46aZWsHcpSuVIkxpl+aKfHwfORNPD34deUwN5qU1tutzjzamLV3hmVjDbj12gUTVPHvWtxOerw1n11xme3hZMMTdnjkTF0r6+F3UrlOSz/eUo3fdjWPUJu796nXoyjJQnAqhSxoOj52L5JzKFgKqeVClbDICI83HcUaJgzTqq8pYlW/g6Skc5so53V2DIg7XxLlWUDa8/QO+AyhRzc6FG+eI83aI6AI81Sm2p1yjnwWer9hMVm8ClK0mUK+7Git2n+GzVfjo1qIBbuSq4PDqOYnXu44dpH9GtWze83BPZdzqag5ExNKhUmsqeqYF/LCpnQzpV4WPJFr5SBUXGPvXHAipTppgb7e+qgIiw4LlmvLN0D23q3UHXhpW4HJ9Eh8/W07pOecZ1b8Bnq8KZuDqcoeMm4xWxjqFDh7IlZDsu7Ybh7l2bet4lqOSZeletiPNx9viIyoHo5GlKWUz6k7jxSclsDD/LA3XuwMlJ2Lp1K9179OTEqVOUafscYd9/QJWyHjQeu4oWtcvxae977Vu8sgSdPE0pB5H+gjN3F2fapJvmuUmTJuzYHkaRKvcQ9ftk3np1MDExMdxVsSR7jl8/NNMYw6ifd7Nuf2S+1q+sSwNfKQdTrlw5Qjas5okhw5kzZw5NmzbljpRzhJ+J5sCZaPpMD+TAmcv8fSqa2VuO8OS32+xdsrII7cNXygHVr1SamZM+pG+XNjz++ON88XJvirUdQttPU19ftvMEcVenhjYm9QY6eXHHM+XYtIWvlAPr0KEDYWFh3HX3XZxdMp6o1V9RzMXw2+5TfL/laNp6B8/G2LHK2yM+KZnEdFcvq1uzZODrbJlKZV/lypXZuH49TwwcTHTwEmJ+HMXu/Ye4HJ/EmEfuBqDtp+tYsv34LfbkWOq89RsPT9po7zIciiUDX8fhK5Uzbm5uzJz+JfPnz+fssQOcnPkScYe38/A9FdPWeXn+dpvm3LcynTQuZywZ+Eop2zz22GOEBAdRplx5IheOYuIn4zHmWrfH2csJN9laFXQa+EoVMHXr1uXw3h3069eP0aNHc2bROyTHpQ7ZDD8dzYKgo+w+rt2lhZEGvlIFkIeHB7Nnz2bKlCkkR+wi4YfhxJ/YR7+vt/LGj7t4cd61++jO3Xo07Z6+qmDTwFeqgBIRBg0axKaNGynq5sKp798gOvRXjDGcuBDH4DkhLNt5gjcX7+K52SH2LlflAw18pQq4Ro0aERISQotWDxD1xxSKbf6SuNgYVuw+xZC5qS39lAJ2MldlTgNfqUKgbNmyrFv1G++99x77Nv3GqVnDSDx7LO11NxeNgsJAj7JShYSTkxNvvfUWH369gOS4i5yc9Soxf60D4Pj5uOuGbIafjib06HmORcUyc9OhtNdSUowlhnampNi/Bkekga9UIdPqwTZ4P/k5bndU5+zSj4j6YyqX464QMGYVH/72NwDtPltP9y83M2lNOG8v/YsVu09hjKHTxA28MDfUzp8AkjTwbaJz6ShVyNzjU4reLRuyu8YXbJk/kejgJcSfDCf5kTf4cm0CF+MS09YNP3MZgOe/D8W3Smn+PhVtiYud9JyDbfKthS8i9URkqogsEpHB+fW+Sqnrubs48+lj91K+VDHKtBlImyHjSTx3lJMzXyHuYAjfb702B0/Y0QuZ/pz+S8EeNO9tk63AF5FvROSMiOzOsLyjiOwTkQMiMuJm+zDG7DXGDAJ6AzdMzK+Uyl/31SoHwOSRz+P9xASci5fhzA9vc2Hj92x8vRU+V++klZkDV1v+9qItfNtkt4U/E+iYfoGIOAOTgYeA+kBfEakvIg1EZFmGxx1Xt+kKbARW59knUErZZFDLmux5pwN3ehXHtUwlRk77kfseepSLm+bxdN8eVCySlOW2i8Mi8rHSG2ng2yZbffjGmPUiUi3D4sbAAWPMQQARmQ90M8aMA7pksZ9fgF9E5Fdgrs1VK6VyzclJ8HBPjYDD4zsDMPqRRXz99de8+OKLuAbvoHin4bhXqnvDtnO2HKVccXeeb12LVxaEUcerJC+3rZ1vtes5W9vkpg+/EnAs3fOIq8syJSKtRWSiiEwDlt9kvWdFJFhEgiMj9dZsSuUnEWHgwIFs2rQJcXbm1NwRBFzect1QzOn9/QGYsCqcIXNDWb7rFJ+t2s/l+CRW/XUagKTkFMJPRxN8OOq21GmFoaGOKDeBn9ntc7I8CsaYtcaYl4wxzxljJt9kvenGmABjTED58uVzUZ5Sylb+/v70eHcORav78uPkMfiETSclPhaA+2tf+71ceTXgAQbNDuGZWcGEHInio5X7aPfZenpODbwt8/RoC982uQn8CKByuuc+wInclZNKb4CilP3VreZN+R6jeH74KLas/pWTs4aSEHmEIq6Zx8bGA2cBCDp8no3hZ9OW/3XyUqbr54b24dsmN4EfBNQWkeoi4gb0AX7Ji6L0BihK2d+rbe/kmycbM/nDd5n941JSrlzm1OyhzJ1789NvkdHxVCx9bYTPuatz8C/dcYL9edTa18C3TXaHZc4DAoE6IhIhIk8bY5KAIcDvwF5goTFmT14UpS18pezPzcWJNvW8AOjeuUPq1bletfi///s/zq38EpOU+Vj8WYGHKeF+bTzI+dgEYuKTeHFeGI9O3pQntWne2yZbgW+M6WuM8TbGuBpjfIwxM64uX26MudMYU9MYMzavitIWvlLWUsTVmQf96zJz0VJee+01Loctp8jv75B08cwN6yYmm+u6cc7FJLD92AUAYhKS86QebeHbRqdWUEply6wBjQHo0+QjmjdvzpNPPsnlmS9Rrssw1nz+Kj8ERxAZHc+vu07y96loGlcrw/nYBKIuJ7Dt0LXROonJKbg65+4ifz1paxtLTp6mXTpKWdujjz5KSEgId9aszplF7zDvy48Y1bkub3Wpl7ZOnQolAPhtzyn+3HftL4Fdxy9y+tKVtOcpKYa1+87kaKilzpZpG0sGvnbpKGV9tWrVIjRoK0899RRjx46lQ4cOSNy1rpwe/j5pk6/tjLhIjXIeAHT/cjNN3k+92P5M9BW+3XyYJ78NYt62Yze+SRa0R8c2lgx8pZRjKFq0KN988w0zZsxg06ZNBAT4cyViLwA1y3vQtt4daev2Cqh83bbHomJpPHY17y37C4Bth85l+T6bD5wlPula/7/24dvGkoGvXTpKOZYBAwYQGBiIu7s7kfP/S60z6yju7sLkx/3S1ulyj/d129z/4Z/XPb+QxQych87G0O/rrby1+NrcjRr4trFk4GuXjlKO59577yUkJISHu3Rm9bcf0bt3b+JjY9jy3zb8+VpryhZ3u+n2a/dFUm3Er3zw299UG/ErMzYe4sCZyyQkpQAQlG6aBu3Ct40lA18p5ZhKly7N4sWL+fDDD1m8eDEBAQFEHg2nejkPirllb1DglLX/APDesr9o++k6EpNTA//wudi0dXQuHdto4Cul8pSIMHz4cNasWUN0dDRNmjRh1qxZNu9vy8Frfft/nUg9KawtfNtYMvC1D18px9eyZUvCwsJo3LgxTzzxBM899xwmKSHH+xnz6960nztN3ABAUkrKTbfZc+KiDt3MhCUDX/vwlSoYKlSowKpVq3jjjTeYPn06Jf94lwV9axA6qh2Vrs63M+aRu3Fxymzy3cwFH46i88SNWb6+M+ICnSduZMq6f3Jdf0FjycBXShUcLi4ujB8/niVLlnDsyCE6tm7OpjW/s2hwM8Y+ejf9Glfh8z6+t9zPv18QPacGXrc8KTmF6CuJpKQYwk9Hp13VG3LkfN5/GAcnVj75ERAQYIKDg+1dhlIqjxw8eJCePXsSFhbGiBEjeO+993Bxuf5kbrURv2a6rauzkJicdV693KY2n68OT3veolY5pvX3p8n7q5nY914erOuVNx/CAYhIiDHmhnuHawtfKZVvatSowebNmxk4cCDjx4+nffv2nD59+rp1JvbNvLV/q/l3Av+5/sKtU5euMH7F31yOT+LzVeFZbFW4WDLw9aStUgVXkSJFmD59OjNnziQwMBBfX182bNiQ9nrXhhXZ+MYDdG7gTeiodgxuXZNyxd1vebHVtgy3Uzxw5jKztxwBwMPdhV0RmifapaOUspudO3fSo0cPDh06xPjx4xk2bBgimZ/A/Sk0gqELd+Tq/T7u1ZCe/j4cPReLCFQuUyxX+7Mq7dJRSlnOPffcQ3BwMN26dWP48OF0796drP6y7+7nk6PRPJl57YcdJCWn0PKjP2+Y2qEw0MBXStlVqVKlWLRoEZ988glLly7F39+fHTsyb8n//MJ9aT8Pbl3TpverNXJF2s8BY/64blK2gk4DXylldyLC0KFDWbt2LXFxcTRt2pRvv/32hvXurlSK+2uXA6CyZ+67Y85eTmDyn/8QFZPA7MDDafP2FFQa+Eopy2jRogWhoaE0b96cAQMG8MwzzxAXF3fdOpHR8UDqPXf/tfj55mk/H3y/E9tHt6Osx80na/vXxNXh+L33B6OW7OGrDQfZfuwCIUei2Hcqb264biWWDHwdpaNU4eXl5cXKlSsZOXIkM2bMoHnz5vzzz7WrZv+JTL2pSqXSRXm9Yx2Gd6iDbxXPtNednITSxdwIGdXuupupZ8eZS1d4ZPImekwJpMOE9fSfsTXLdRcGHeOpb7fl8NPZlyUDX6dWUKpwc3Z2ZsyYMSxbtowjR47g7+/Pzz//DIBv5dRwb1y9DM+3rsULD9QCYNmLLZjcz++6/bi7pkZck+plsvW+YVdvtv6vDeFnmbftKACLQiLoNXUzcQnJnLscz+s/7uTPfZG2fkS70GGZSilLO3z4MD179iQkJIThw4fzxqh3iIxJpG6FkrfcNuzoeSauDmdqf3/CT1+my6Ss5+C5mbWvtab1x2sB8PEsSsT5a91Mh8d3tmmft5MOy1RKOaRq1aqxceNGBg0axEcffUT3Lh0pZWKyta1vFU++faox7i7O1LqjOOWKu/HSg7XYPrpdjmr4N+yB68IeUi/wchQa+EopyytSpAhTpkxh1qxZBAcH4+vry9q1a3O2D1dngt9qx9D2dShdLHsndLNjwMygPNvX7aaBr5RyGP3792fr1q2ULl2aNm3a8MEHH5Byi7nxb7dLVzK/F68VaeArpRzK3XffTVBQED169GDEiBE8+uijnD+f86mQ5w1sSpU8mFohLsFxLtzSwFdKOZwSJUqwYMECJkyYwPLly/H39ycsLCxH+2hWsyxLXriPPo0q3/Da3ZVufUL4XxYe93IDSwa+jsNXSt2KiPDyyy+zfv16EhISaNasGV9//XWObnDu6eHG+B73cHh8Z/yrXhvLv/v4pWzv41azeFqJJQNfx+ErpbKrWbNmhIWF0bJlSwYOHMiAAQOIjY3N8X5+HNycA2MfInTUtRE82ZmsTQNfKaXyUfny5VmxYgWjR4/mu+++o1mzZoSH5/ymJy7OTpTxcEu7nWKzmmVvWCd87EM0rnbtQq4WtcvbXng+08BXShUIzs7OvPPOOyxfvpyIiAj8/f356aefbNrX108E0LdxFXw8U4P/pTa1aV/fiz9ebYmrsxOfPtaQZjXK4ubilPbl4Ag08JVSBUrHjh0JDQ2lbt269OjRg2HDhpGYmLOhk/W8SzKue4O0aRyaVi/D9P8EUNurBAA+nsWY92xTPIu55uicgb1p4CulCpyqVauyYcMGXnjhBT799FMeeOABjh8/nuP99ArwYdXQljSvVS7T1wXRPnyllLI3d3d3vvjiC+bOncv27dvx8/NjzZo1OdqHiFDrjhJZvu4kOixTKaUso2/fvmzbto2yZcvSrl073n///Ty7OldESNHAV0op66hfvz7btm2jd+/ejBw5kq5duxIVFZXr/To5oX34SillNcWLF2fu3LlMmjSJlStX4u/vT26nX9c+fKWUsigRYciQIWzYsIHk5GTuu+8+pk2bZnMr3UnQLp2siIiHiISISJf8fF+llEqvSZMmhIWF8eCDDzJo0CD+85//EBOTvTn203MSwYHyPnuBLyLfiMgZEdmdYXlHEdknIgdEZEQ2dvUGsNCWQpVSKi+VLVuWX3/9lXfeeYfvv/+eJk2asG/fvhztQ6RgTq0wE+iYfoGIOAOTgYeA+kBfEakvIg1EZFmGxx0i0hb4Czidh/UrpZTNnJycGD16NL/99hunTp0iICCAH374IfvbixS8k7bGmPVAxlPajYEDxpiDxpgEYD7QzRizyxjTJcPjDPAA0BToBwwUET1/oJSyhPbt2xMWFsbdd99N7969eeWVV0hISLjldiJg5/uv5EhuQrcScCzd84iryzJljBlpjHkFmAt8ZYzJ9J9JRJ4VkWARCY6MdKw7wiulHFflypVZt24dL7/8Mp9//jmtW7cmIiLiptuk9uEXsBZ+FjKbN/SWn9wYM9MYs+wmr083xgQYYwLKl3ecWeiUUo7Pzc2NCRMmsGDBAnbt2oWvry9//PFHlusXpguvIoD0t4rxAU7krpxUegMUpZQ99e7dm6CgILy8vOjQoQPvvfdeplfnpk6t4DiJn5vADwJqi0h1EXED+gC/5EVRegMUpZS91a1bl61bt9KvXz9Gjx5N586dOXfu3HXrOBXEFr6IzAMCgToiEiEiTxtjkoAhwO/AXmChMWZPXhSlLXyllBV4eHgwe/ZspkyZwpo1a/Dz82Pbtm1prxfIYZnGmL7GGG9jjKsxxscYM+Pq8uXGmDuNMTWNMWPzqiht4SulrEJEGDRoEJs2bUJEaNGiBZMnT8YYg4jobJlKKVXQBAQEEBoaSrt27RgyZAiPP/44KfGxBa+Fn9+0S0cpZUVlypRh6dKljBkzhgULFrDuo2e5cOKQvcvKNksGvnbpKKWsysnJiZEjR7Jy5UoSLl/gt7EDmD9/vr3LyhZLBr5SSlldmzZtaDtyJqV9atO3b1+GDBlCfHy8vcu6KUsGvnbpKKUcQbHSd/DAsC8YOnQokydPpmXLlhw9etTeZWXJkoGvXTpKKYcgIE4ufPLJJyxatIi9e/fi5+fH77//bu/KMmXJwFdKKUeQfn6ZHj16EBwcTMWKFXnooYf43//+R3Jyst1qy4wGvlJK5UL6QZl33nknW7ZsoX///rz77rt06tQJK00CacnA1z58pZQjEOGGKSOLFSvGzJkzmT59OuvWrcPPz4/AwEC71JeRJQNf+/CVUo5AyHx6ZBFh4MCBbN68GVdXV1q2bMnEiRPtPtGaJQNfKaUcgQg3nVrBz8+PkJAQHnroIV5++WX69OlDdHR0/hWYgQa+UkrZSOTWNwHx9PTk559/Zvz48SxatIhGjRqxZ0+ezDOZYxr4SillIyF797R1cnLijTfeYPXq1Vy4cIHGjRszZ86cfKgwQx35/o7ZoCdtlVKOIDst/PRat25NWFgY/v7+9O/fn8GDB+fr1bmWDHw9aauUKqi8vb1Zs2YNw4cPZ+rUqbRo0YLDhw/ny3tbMvCVUspR2DLwxsXFhQ8//JDFixezf/9+/Pz8WL58ed4Xl4EGvlJK2UhEctSlk9EjjzxCSEgIVapUoXPnzrz11lu39epcDXyllLKRgG1N/HRq1apFYGAgAwYMYOzYsbRv354zZ87kSX0ZaeArpZSNcnrSNitFixZlxowZzJgxg82bN+Pr60tISEge7Pl6lgx8HaWjlHIEQq4b+NcZMGAAgYGBVK1aFS8vr7zb8VWWDHwdpaOUcgSpffh5O13Cvffey6ZNm/Dx8cnT/YJFA18ppRyB3HoV2/Yrt2fPGvhKKZULdp4PLUc08JVSyka3mjzNajTwlVLKZrkbh5/fNPCVUspGqS18x4l8DXyllLLR7Tppe7tYMvB1HL5SyhFoH34e0HH4SilHIA7Wxrdk4CullKPI6wuvbicNfKWUspF26SilVCGRV5On5RcNfKWUslF272lrFRr4SillK23hK6VU4ZB6AxR7V5F9GvhKKWWj2zWr5e2iga+UUrngQA18DXyllLJV6h2vHCfyNfCVUspGOiwzCyLSWkQ2iMhUEWmdX++rlFK3S17f0/Z2y1bgi8g3InJGRHZnWN5RRPaJyAERGXGL3RjgMlAEiLCtXKWUso7bcU/b28klm+vNBL4AZv27QEScgclAO1IDPEhEfgGcgXEZth8AbDDGrBMRL+BT4PHcla6UUvblaC38bAW+MWa9iFTLsLgxcMAYcxBAROYD3Ywx44AuN9ndecA9qxdF5FngWYAqVapkpzyllLKPQjSXTiXgWLrnEVeXZUpEuovINGA2qX8tZMoYM90YE2CMCShfvnwuylNKqdvL0aZHzm6XTmYy+6RZftcZY34CfsrWjkUeBh6uVauWjaUppZTKKDct/AigcrrnPsCJ3JWTSm+AopRyBIXpnrZBQG0RqS4ibkAf4Je8KUsppaxPKIDj8EVkHhAI1BGRCBF52hiTBAwBfgf2AguNMXvyoii9p61SyhE42g1QsjtKp28Wy5cDy/O0otT9LgWWBgQEDMzrfSulVF4RHGscviWnVtAWvlLKEThaC9+Sga8nbZVSjsDBZke2ZuArpZSjcKAGvga+UkrZTrRLJ7e0D18p5QjEwe5xaMnA1z58pZQjcLTJ0ywZ+Eop5Qj0Bih5QLt0lFKOoIyHO5VKF7V3GdkmVp4HIiAgwAQHB9u7DKWUcigiEmKMCci43JItfKWUUnlPA18ppQoJDXyllCokLBn4etJWKaXyniUDX8fhK6VU3rNk4CullMp7GvhKKVVIaOArpVQhYekLr0QkErgAZDx7Ww44m+8FXVOKG2vKr/1kd5tbrXez17N6LSfL7XmMHOH4ZGfdvDpGBfX42LovR/gdyu3xqWqMKX/DUmOMpR/A9EyWBVutpvzaT3a3udV6N3s9q9dystyex8gRjk9+HqOCenxu9zGy5+/Q7To+jtCls9TeBWQir2qyZT/Z3eZW693s9axey+lye3GE45OddfPqGBXU42Prvgrt75Clu3SyIiLBJpN5IpR16DGyNj0+1na7jo8jtPAzM93eBahb0mNkbXp8rO22HB+HbOErpZTKOUdt4SullMohDXyllCokNPCVUqqQ0MBXSqlCokAEvoh4iMh3IvKViDxu73rU9USkhojMEJFF9q5FZU5EHrn6+7NERNrbux51PRGpJyJTRWSRiAy2dT+WDXwR+UZEzojI7gzLO4rIPhE5ICIjri7uDiwyxgwEuuZ7sYVQTo6PMeagMeZp+1RaeOXwGP189ffnSeAxO5Rb6OTw+Ow1xgwCegM2j8+3bOADM4GO6ReIiDMwGXgIqA/0FZH6gA9w7OpqyflYY2E2k+wfH2UfM8n5MXrr6uvq9ptJDo6PiHQFNgKrbX1Dywa+MWY9EJVhcWPgwNUWYwIwH+gGRJAa+mDhz1SQ5PD4KDvIyTGSVB8AK4wxoflda2GU098hY8wvxpjmgM3d1o4WjpW41pKH1KCvBPwE9BCRKVhv3pDCJNPjIyJlRWQq4Csi/7VPaeqqrH6HXgTaAj1FZJA9ClNA1r9DrUVkoohMA5bbunOX3FaXzySTZcYYEwM8ld/FqBtkdXzOARoi1pDVMZoITMzvYtQNsjo+a4G1ud25o7XwI4DK6Z77ACfsVIu6kR4f69NjZG239fg4WuAHAbVFpLqIuAF9gF/sXJO6Ro+P9ekxsrbbenwsG/giMg8IBOqISISIPG2MSQKGAL8De4GFxpg99qyzsNLjY316jKzNHsdHZ8tUSqlCwrItfKWUUnlLA18ppQoJDXyllCokNPCVUqqQ0MBXSqlCQgNfKaUKCQ18pZQqJDTwlVKqkPh/pOHuhCN2aOIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(f_1, SU_1)\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "x = np.linspace(10**(1), 10**(3), 1000)\n",
    "a = 3\n",
    "plt.plot(x, a*x**(-5/3), c='black')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f474e56e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}