3D_Energy_Spectra_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1e0c9509",
   "metadata": {},
   "source": [
    "# 3D Turbulence Energy Spectra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "id": "156044ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import numpy, and matplotlib\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "351478e7",
   "metadata": {},
   "source": [
    "This piece of code generates an artificial flow filed which represent the velocity components in a homogeneous and isotropic field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "id": "f0b34c42",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 166,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAGKCAYAAAB+eFQVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAACFVElEQVR4nO29e5Bd13Xe+R10E000gVbjQQDdJKiLIUCCMAlCfAg0KL6pkGaZiqWyLdGVyPGMrUx54lSlKpV4ylNJRhlNJimnPHGSSqJkLMV2jVnj2ExERaIsWoYIERRFQGxRFF4EiCZBogkIjyZBEg81cOePc/ZZ37n323efS3ZD5OH6qlh39T6Pvc85P/D2+XrtvbJ2uw2Xy+VyuZqmeT/tAbhcLpfLNRfyLziXy+VyNVL+BedyuVyuRsq/4Fwul8vVSPkXnMvlcrkaqcFeG38H/6hMsfx7+L2yfdmjb9pOx4rPa6zpxKYFZTyJVhkP4FwZj+NQfq4pOheN5sClY2W8DZsBADuxXp5rM7aV8f0/2GIn+TfF57etqXLFt1L8KQvP3JF/Hh1eWrYN4Yzs+zCW0zjzEz6KB8u2redvK+MTPx4t4/UrdgIAPoFHy7ZfwCO2/czOMp4eWmznw21dfTxZ3B8AOPLGChvnoI1z6OJ8/MdeW1a24fEhix+2EFso/knx+XFq+/sW3nn3Y2UcngPfH35m22icU3tX20lm8o+FraNl0/jwIbsOOt/kG60yPvWd4r78kMZ2qYULP2Pn+7vDv1/GgWXJMZBkWXEMaJYVx0CaZckxkGZZcAykWVYcA5plxTGQZllx3NmHYllxDNRgeYs1lRwDkmXFMZBmWXEMpFmWHANJlk8OL8swR/rfs2zWUur/cbs9Z+PsR/4G53K5XK5GqucbnMvlcrk+GGril0HWa6L3USwqNy77Z2S//DntdLz4XEttbP1xO73C42Dx+XKk82spvif/eG3dh8qmQbIOln2LxvZvLXyyGOfjkS7up3jTLfRDGP9yajtHMV/HOMV35x9/tfpny6YncHsZT2O0jNcjt20+SVbOsofpOibovEsovi//eOZ6u0FbqY8jNGi2V+bjLADgGMyqehqbyvh7O+0clRsWrvVea3pww5+W8d/Cl8o4WDuDEVuH+9uDq8v4HAYAVO2+VSUgVbHl/Sw+0tW2CCfL+AF8rYz/LsyiLFlWHANplhXHgGZZcAzUYFlwDKRZlhwDaZYFx4BmWXEM1GBZcAykWVYcAzVYVhwDkmXFMZBmWXEMpFlWHHe2M8v34RsAgC/g83Nm/f2LWbQo/4FblC6Xy+VyzZ2a+Fbqcrlcrj7VxC+DntdUyTIjm2TLdouDs3Pdi9a21t7Oq3YP6+WOz069QHHRycqNr1sbWw7fsvDVGnZO0GMUz3zX4luPFAHbOtzfGYrZ2imGd9dDT5VNC8ZOlbGydmL3uGJRLqX4rfzj5qHny6ZF68zKOALKoiR7JcRs6yyl1MGB9bbvjpU3lvHZ03l22kfHny7b2Pq778w3yviSvzifB3Svbt/4vTIeXT1dxi1MolNs63BcyaIkC2c5jnS1KVsH0PdZcQz0yfLLkThIcAykWe6HY8BYlhwDaZYFx4BmOWZRJlkWHANplhXHQJplxTGgWZYcA0mWFcdAmmXFcWf7KKbL2Fj+vOxvNnTRnJ35pye3KF0ul8t1QZVl2R9kWXYky7Lnqe2fZFn2apZlE8V/D7zbfpr4VupyuVyuPnWBvwy+jHyG5x92tP9eu93+3dnqpPc18eTX4zLE4eJzBbWtJYuD7Y6KNVKc5DidbJBGM8J9v97x2XkuaufD+hFf00+Kk1xEc0crtiT3fQnFxzo+AYyOnZD9jZ6b7toXb1DM1wrRTsctgllDpzBcxtXss/wCONOL7bxKPGrxzMzb+XjJLuF9LzlOdk4YG98fuiY+xyjsvpwrMOTz8r4stU/sOG5XLCuOgRosC44BzbLkuDMWLM8mx0ANlhXHHbFiueS487gUy9yWYFlxDKRZVhwDmhHJMZBkWXHcOR7Fcmx7jN8Ky3OkC2lRttvtJ7Isa811P25Rulwul2tWlWXZ57Is207/fa7moX8ny7LnCgtzcXr33vIvOJfL5XJhcBb/a7fbX2y32zfRf1+sMYR/B+BKABsBTAH4l7NxTXHRmnycQcZZZsHOueESvW80i/Li/GMJWyc8GnWO1KRxABsoU+uRPnye6yi+KFw3Z57FLErOPivG99Za+71hP9aU8UksKuPhgdwy2XDN3q7jAYCcmmr2WbFPm54NZ17xmoKDYoIsb+fjOD62+zI7eXGtkxts+wHa99UxG9xla49VjgGAM3RN+3El9Wdr+AWr6Szm23EU83UcEGOOZZ5xu7rPimOgT5YvtlCyHDs+wfKscgykWRYcA5plxTFQg2XBMZBmWXHM2zuPC7HiGNAsS447jlMsK46BNMuKYz4voFnegLnTTzuLst1ul38pyLLsPwL46rs9p7/BuVwul+unrizLxujHTwJ4PrZvXXkWpcvlcrku6JdBlmV/AuBOAMuyLHsFwD8GcGeWZRsBtAFMAvjb77afntfEpUIW32oTlnnya5llxpbEJopjFmWwRFZFtrPXUqyvd+J6G8/AObMtRt4y24L7vqOYj8wVRli0JB1aG+mHcA72rdhGOk0xWzsfyz92DN1UNj0LOzFPkH27yBBbtcn81cUv2D0G211s7RT3YmKJ+ToT1Mfh6ETv/AKOwUqM8Nj27iXz4zvUX3Gtewdt+8R6WztvB9X3GN60Ne+Lns3OAVu/j9fc24Or0Ckee+w62M7ZjhuLNrOJONuMJ9DeVowNMJYlx0B/LPPzVywLjoEaLAuOgTTLkmMgzbLgGNAsK46BGiwLjoE0y4pjoAbLgmNAs6w4BtIsK4557J2xmugdOM7be7P8Cdnb7OgCZ1E+JJr/n9nuxy1Kl8vlcjVSblG6XC6X66eeZDIX6vkFx6/Ri9fusg1s1SiLso8syooVx6LjQvbSPsowGhowW3LDWp29taG4ukmyZPghXkcZcHLMdSzKKyw8sSa3naoZUmY5sLUzjFPFdtu3co9prT62dmaKsXEfHPezFiWPDbuh43Ctl1vT5HrddygLMjCgLcXJyH0JOkswcBYaS52DqyMvGjkZ2dfi8j4rjjvb+8iilCwLjoEaLAuOgRosx8abYllwDGiWFcf5dttXsiw47uxDsVxnLUrJsuIYkCwrjoE0y4pjIM1y7N9CHZbnSk1823GL0uVyuVyNVBO/tF0ul8vVpz5wFiXbARVbgxVrV1KVhLltQGwnDVZ2fveaiY1ddVPj+kPG1eCAzvqqxucqn70Hmt4lpdB3rf5Ygx2fXZvf2TPhe3Eu8XtWzKIaEDcmdn1Jlt8tx53tA2I76YKwHOsicf2cOahYVhx3xrX77VOxvpMajMRl09xznO/f/e9ecdy5z4VQE9923KJ0uVwuVyPVxC9tl8vlcvWpD5xFyZVoeb3HSuXiUJ+Ds8l40ii3q/UjYxW9aT3AoWKfsXU2noqlcECcF8Cuor+XIl28QOPZpK6PM69q2FIjL+fZcOOruZrvVBlzxlm4t5V7zNfB8RILB4v2VdfbgLmPwYhlFNbwG6K1/LjvV9bYmpmVIsXhWlvWxMeNiZjHwBl3nJ3GE4SDVlDRmlUV4EyckRYm0J4dsYw1nhzL55AsK46BNMuRdVAly4JjoAbLgmMgzbLkGEizLDgGNMuK485Ysiw4BtIsK46BGixP0hj4mbUsDMcpjnkMgGZZcQykWVYcA/VYnis18W3HLUqXy+VyNVJN/NJ2uVwuV5/6wFmUy6aoZgvbHRSHKsaVUiG8Jh+3szXycsdnpz7Uve/KpTQbd6Z7e2f8QuTUQfso3sTnCOOPlciJxYWLMLZa2zbDeLurvc49rkz6LtrHrzcrh/vgjCxVZmSILoqPW9g6WsZvtmyNv3C6BS2rXDwGbcGu/PHrlWMAYHxM3wsuuRLKjETtLtIZZVFSG9s6fA55nwXHQA2WFcedcZDgGKjBch8cA8ay5BhIsxyxXRXLimOgBsuCYyDNcqxcToplxTGgWZYcdxynWFYcd45Hsaw4BuqxPFdq4tuOW5Qul8vlaqSa+KXtcrlcrj71gbMo63z9DSYmAlc0kNjO50jtGzuO9I4fWDhfYuJ5l4r9ecLnuciFlO2xCaixezHY+7wpzcTGMxM53096b+9nHLG+U/vGJuGGvqPXxO0JPgcTE4ErqnMZiqE6En2/a447x5FimfZNsRy9xyoWHMfOW0d9sfyT3tv7HcNsssx912J5jtTEt50mXpPL5XK5+lQT3+D8b3Aul8vlaqR6vsEduHSsjFdfa1lGnNY1cqwIuDQHVzHmdrZGwgRYzjLjt/BrKC7Ot+/Sy2lXe9VffR2N7VoLP7o9/7Qpl9XfUj5KcWXMoW+e5BubKEv7tIvj9lMpFC6LMo3FZRzGf+BSqyRcucd8r/geFdfHlYT34OoyPoLlXX0AwPwi44yrIPNxpyZsbJig/opxnB0bKZv23m7H8fWtuPRwV7/cx36soXYbf7DBOCMtNoGWy4bsK87H93sU02Xcopm+kmXFMZBmWXEMaJYFx0ANlgXHQA2WFcdAmmXBMaBZVhwDNVgWHANplhXHQA2WJ6g/emaKZcVxZ9+KZcUxkGZZcZyfN83yXKmJb3BuUbpcLperkV8GblG6XC6Xq5Hq+aW9DZvLePU9f2YbeN0+VdH7FgvPqCrIoHX5YhO9yWp5Ziz3M3ZifdlWWXPx+q1lfMUdR8p47fP554Nk8fBreIvGiTsovr34jFVBpkmzbZpMu23Jjfkn3bdtuLWMeS27EPOadcvveKSML/nQeTsxVUJ+blVuiTyNTWXb02S2TpHXxJNiVUXvvc9tsBM/Dh2H7DO6/m3X2vVdtWRPGau1A39ID5LHzHZPEE9mPViZpWxia2fi0EfyYLfNzD4yaqWpR2+YtjErlhXHQF8sD6UmeguOgTTLimOgBsuKYyDJsuIY0CwrjvMuEiwLjoE0y4pjoAbLimNAsqw47uxPsaw4BtIsS46BWizPlS5q4CtcAy/J5XK5XP1qsIHfBm5Rulwul6uR6vmdzTbKa+vsfX/lRvJzhEV54voFZczZSfzqH8qFVNbkI3GWWRjHLhoPWwdLYSlwV9xk1g42FkPjLrgUykaKzZXBmcIxODpsFkhsjUdeRy6Mk+/bzrOWkvb6UcpULJzEq2HWyFVDFq/fuLOMp4fsOHUvuL+p42RRXkzWzmA+/tdfo/X5yO6qxM9QHC511JrObrcstJ1/zfpeXGR98bPZGRvnztXo1HTLOpkeHu3aDgCTx1s05sLO2UU70C3euY76Hu5mWXIMJFlWHAOaZcVxPuTeLCuOgRosC46BNMuK4664YFlxDKRZVhwDaZYVx0ANlhXHgGRZcQykWVYcA2mWJcdAmmWdXDwrumju55JfcDXwpdTlcrlc/cotSpfL5XK53ifq+Z3Nr+eVNdQS5TYGztm+QwNnoVSeL7Iunup7IJLdxJaLLAvC9sRpsb3juIGZ7vNWSnacI8tkwPYZEOMcHLR4XmXfme59OT53vmvfaB98LwY5pvGHdtpe+ZWNSVAzPrmtclj3M4llvXGsyOPxVvat7CPOERlb7HySvUicYjn676JXv6jBcqy8TYplwXHlvNAsK467xl/ce8Vx176C5di+KZYlx/mJKRYsx2YuC15izybJcuT/oCmWJceRsXWeb67kWZQul8vlaqYa+Dc4tyhdLpfL1Uj1fIPbjG1lvOxbVK33W7RTyOqiKsAjb9lr/Ya1e/XJDxSfPDmWRsPrS4bJr5wtyXbH5jM2Tvwlna8Y5xbqgx2AW9kGojUFgxuwcjmlrEUso5Hxl+yHYh4vT0A9N89+LTq5wtanu6rIOON7/JGDlEI1YeHKJTaOzbdu6+qDdXjEsuH4Hg0VVsvhcdu+4+ct5e7sm5YZWanCHq77Xmu6/G6rhc7jvw1bu/odhVVP5vb9V9mk4FAKJFYpnI9bNWKg7bg/H/+x1mVl27ylVjZ683wbm2RZcQz0x/IBigXLimOgBsuCY6AGy4JjIM2y4hjQLCuOgTTLiuPOPliBZcUxUINlxTEgWVYcd/atWFYcA2mWFcdAHZY/hTlTA/28Bl6Sy+VyufpWA78N3KJ0uVwuVyOVtdvt+NYfZLbx89b86p9bHIyWDewybKI4shZlaQPF1qIkm6RcX+8mamPLge2cf2XhFwqXIJZ/tIDiX6V4RehvOTXGsteuoPie/OPoLy4sm57AbWXMZUbWI5/8essPJuz4P6JzUTOWUHxf/nHmIWvaNvyzZXwouhZlfhe4xAivAfh4GDyAl76/zk5eXOvlm82W/Cz+sIx/DV8q4zU/eCUP6Ia/sXG+jXPA1jXkMiNBbOXwWn4sXsPv2WJ28yRssu0inCzjB/C1Mv6lHzxqJylYVhwDfbJMdqZkWXEMpFkWHANpliXHQJplwTGgWVYcAzVYFhwDaZYVx0CaZcUxoFmWHANJlhXHQJplxXHe3pvlX8JXMtnhbOjKrMeXQZ/a3567cfahBr6Uulwul6tveRaly+VyuVzvD/W2KH/DXlmf/E/W/Ljal8TOyIbIO+Ku4tWfiipHq22vDXbORmpkr4ayzL7wst6lH91QfF5Gbaco5sobLY5DZRyyXPCAhW2qaJxNFMGfWNsbf2zxk2SpcKWTGzaKPshSqlRuFhPZz9D27cOWvfUN8oy204KGoUrxJjxdtrFFueZhsnOCvcb90nqIuNvCN64xu+fcQP6r4+KDdJfZ7mNLjWzC58Zye4jtHq6CfPvB79nOZLEHllMcA2mWd9G1KpYlx0CS5dnkGEiz3KK21q30g2BZcgwkWZYcd/ahWI4syJBiWXEMADfRwpXBmpQcd/YtWFYcAzVYFhwDNVheNYfW3/pZtCh3ukXpcrlcrveKGvht4Baly+VyuRqp3t/Z37awjp0jDsNkxF95STeXOkxxqGJcKRVCttWWWbBzWN/v+OylH1L8y0/mnyvYJqSJtxlZO2VmGd2sR8mW3AetBcVx13Af/BQTFuUQbd90644ynh4ZLWNeO+8scguGbZ01u8nO4YcdYrYU36CY/og9coTWdQzjZyuHsxP5OighbcOt+cTr5auMlsVvEyRqbJhdlt8Jx0Ca5dnkuDNWUhwDmmXJMdAXy4FjoAbLEYsyxbLiGKja7SXLEVZSLEuOgTTLgmOgBst/A3OnBiaZNPCl1OVyuVx9q4HfBm5Rulwul6uR6v2dPQvf6LGKFf0cV8Zcwfi03nc2xZdfZ7J4WW6Cx8knuVjEkYo1sXEM9mxEdf09JbIhzgzZ7ze8jp6KZ/hA7oNjdQEDke1qnLx9IBKL/jhD7syQWVFDgxH76B3qnXAmOQbecyxLjgHNsuK4oxN1u2PIzhXLMaYlyymOO/or94mNMcVypL9aLM+VGvgG18BLcrlcLlffauDf4NyidLlcLlcj1fsNjiZ83k/VLx5LnJSqUeA6sjBmeFJsEceyBXmCbOuWIthIjXQuLnvzJGUv8YTWfqQmevPkbrZ4WhQvCfeLJ8pSPENrHA6G7MoXre3nqPTKUrKteAm/tRtFHzwbOZZFWZzvDM223TFkM4+30+KIvK5fsHY4I239alt/cN3HKI8wPAfudyPFNOYT15opFvpYNkUlmei+VERL/z13afdE70UDtn7fXXc8JfsOLKc4BjTLimNAsyw5BpIszybHQJrlFrUtifAbYskxkGRZctzZh2JZcAykWVYcA5plyXFn32LMimOgBsuCY6Amy3OlBvp5Dbwkl8vlcvWtBn4buEXpcrlcrkaq93c2FY/d9LzFM9+1+HjxeR0d1tpIP0TK5WwqJj9uipXL4RMG2+JGaotU4/5VKjPy7yOnVrqfxxbG/A7L5fB9e/JSG/Q0Rst4/e25NbL6LSurMWIFfHHvBJ1XlRmhPn605soyPoSxMh6qlBnJL4CrJ7MtWVmL8rCN+fyZPN1r+gobO5f/GP+k1UUZuaToL2LrPLvqmjLeg6vRqfGxQzJmHSALZ6JYHLBi61CJkYExe2i3f8rWpQwsK46B/lgOHAMRlhXHQJLlWeUYeMflchTLimOgBsuRcjkplhXHQJplxTGgWZYcA0mWFcdAmmXFMZBm+XbZ2yypgW9wDbwkl8vlcvUtz6J0uVwul+v9oZ5vcGfIUhmiTKdbj1j8k6IU8kXmPsWrILM1EuycVZHR8PmK9/IzVK5igKyDQYpXTFh8Q7F+W2wdPs4yq9g5qqJ3LLOKsxaLe/TspTZ4LtnBlZDfxnB++B1WaXqILS5e74+tnWJsu9d8uGzibLEjNGi2c+YXFxCrgvy9Q/TQttAs1OJad9xrFtD6cbOl1g9YvLnI9OJns3PEMsS4v30wKypMbq1YnxTzdahKyPuxpmzjEiNLqU73pjvMogwsK46BGiwrjgHNsuAYSLOsOAbSLEuOgTTLgmNAs6w4BmqwLDgG0iwrjoEaLAuOAc2y4hhIs6w4BtIsxyp6p1h2i7I/+Rucy+VyufIvuNn6L6Esy/4gy7IjWZY9T21Lsiz7ZpZlLxSfi3udo478C87lcrlcF1pfRjUnCgB+G8BfttvttcjLzv72u+2k53ft0WHLUrpsOXk4ZHdcNNTdVindyzErZG+xXcJ/5GTLZEX3eAbIc1i5nMpK0DjC5Fa2dfiCefJrZfwh5rFz5uQ5sS+N+TAdeITik1hE7fmBde4xJYuV7WzPcDbZUYoHaaDBGuHtfBxeIzvnNeovXPdB23503I7jaw3Xws8mNrYjAowhgoGtKL4OdQ5u48m2lb7VfVYcd7QnWWYuFMuC487xSJYFx0ANlhXHHX1LlgXHgGZZcQzUYFlwDKRZVhzz9s7jSpYVx4BkWXGc99ebZcUxkGY59m+hDstzpguYZNJut5/IsqzV0fzXAdxZxP8ZwBYA//Dd9NNA19XlcrlcfWsWvw2yLPscgM9R0xfb7fYXE4etaLfbUwDQbrensixbntg/Kf+Cc7lcLtesqvgyS32hzbl6fsHxa3bFllMVdmNZhrH4XGL76e52tgvYtqj0Tc1qDT/ugtfkq1xfOF/MllTXQftX7Qk9STW0V+4xXwfHwgZj24PPG7NzrD89nmSJnwW8WWe1hecTGwP3zfcltS9rsI9r4j4ky5FK0UmWY89f7Ss4zseZYDnBcWd3JcuK485xqDFHFjJQLNe6x+p+RuzcFMvx/hIsx8r+CJYVx7y9s4/Qt+I4ti9Lcdy5b+y650w//dedw1mWjRVvb2MAjiSPSMiTTFwul8uV/w1utv57Z/oKgF8t4l8F8N/e8ZkK+Recy+VyuS6osiz7EwBPAbg6y7JXsiz7nwD8XwA+nmXZCwA+Xvz8rtTzpbTy2h+zc5RtczoSs9RxLGG1VF7rz0XWi6P+KhZkjyF0Hidtm4R9yvsPVGyNmUjcbb8k7VxqH6xYINriqdqjM2I7PciYnXO6ezufd1BYMQPnqI8BbbkMiAcfs3NZ6hwDERspFpdd17HVUyzHjguKWIZJlvvguNJ1gs2u9tO999X3u5vjzljel8gYUiwrjjv7kyxHqqYrlhXHQJplxXH3mLtZjv8/oga/c6ULaFG22+2HIpvuibS/I/30XVeXy+Vy/fTVwG8DtyhdLpfL1Uj1/M4+TLMxF49TxVuevHqJaOOYy3CksihZdI72qjAemxA5QIvEjfDYqL9W8flDOi0lUFWqGFfGGfqOlRiJrEV5JjHRm8uMLCvWloveY84f4jX8RB+xWNk9PGF0igf/CvUxSfFPurcfWm/HcX+HivPxs+HyPTwpWE2QHUxYmJ39hXMcOW/nfXveAtq+nI4T91lx3NmuWK5j/YnjA8f5eBIsC46BGiwrjoE0y4LjznGG+604zvdNsCw47uxDxTELPsnyJPXHPq9gWXEMpFmOTfROsRz7f0QdludMDXyDa+AluVwul6tvebkcl8vlcrneH+r5BreN6masu5ssB1r6sXQouEzHxyw8scZeszkjaeTlIsuIKiLzbxBtKjOybUlepmMn1tOu5LNca+G6e2ycrcLP+eUnbftFdMVLqCxIJXcntLPFE5koy3bOtuGfzT+x2dooZmsnxMvJi1x2t9k9y5a/aSem0jnfHdvYdd6tuK2MK9YXWSMhM4zXunvp6XV24q9aiC0Uiwy4p661vq8e34NO8bPZSfViuKTJXqqEfO58/uDH5llZkVUVMEyTWF3Gu17cmAfPZ2XbMVvWEE9storIN2JHGZcsK46BJMuSY0CyrDgG0iwrjoEaLCuOgSTLimNAs6w4BmqwLDju7EOxrDgGarC8xZpitnJgWXGc992bZcUxkGZZcgzUYnnO1EA/r4GX5HK5XK6+1cBvA7coXS6Xy9VI9fzOfhQPlvHq1QfK+K6HrOJtsHbeWmvflTuGrGLuAcoBGxyw1/3x1fkr/Nhqe5Xnirj7qVJusDBitg5nUy37xS9ZfDi3RlawPcOTP9nC+ZSFoYpxKiOxc58wTr5v+7f9jJ34qIX7b8or9w6M62yrq6/dW8YnsbCMtxY1fbmPiZ232Ik5G1JN3uYSIo9T/DDFJ16gH4oxvUJe2zIrN/Lob9s4js0L5XLY1rFntve5DXaOie4ujrWsMMzza6jkNWs3XdR3is9dtH3UwudP32zjvLubZcUxkGZZcQxolhXHQJplxTFQg2XBMZBmWXEMRFgWHANplpnjLbirjL+GB8pYshxbhCDFsuIYkCwrjvPrSLA8obtIsqw4BtIs3919qllTA5NMGvhS6nK5XK6+1cBvg56XtPW8/cF347yJMl4wZuubj46dAADsh/0m9yw2ljH/MZV/qxvHVPHJv/XarxD7Kr/55r+e7jxrv3kNDtpvVufm2XFXYl8Zf+qBr+cBz3HiK6Y3uCcvtQSA7chjnp8SWy7oiPjNt/LWxskbP6a4KMq49dftHo/Om5bn5QKTf1XUA5x4jn7T/Tqd91WK+3mDO8GlNPmExbWesN9w8dUbyvDYtfab6hP359fCz+bYM1Su8xt0Wp7QFXZvUVuLBs/XsZvi8Jsvn+tSiu22Yeud3SwrjoE0y4rjPO5mWXEMpFmWHANJlhXHQJplxTEQYVlwDKRZ5oSUJyiZJMlyv29wJcuCY0CyrDgGarCsOAbSLCuOO8+nWJ7LN7gGfsH53+BcLpfL1Ug18Dvb5XK5XH2rgd8GPS/pxI9Hy3h6BcX8189CbKPF5smwtTNclHEcxtuy72ks7jrH60etbR4toXNyxSJ5XLuYd5PR/Bu2OGbsb8kdY87PwdcUW/1bXiv9Eb5i5/Bcq2IfvsfV66B2GsebIZ7uPldXzE93KLEvjlOs1pw6bCFfE50jPB9+NpVrnoZuD7svpLZRMYSO/sqYhxsZm2JZcQykWVYc53E3y4pjIM2y4hhIs6w4BtIsx/7NSpYFx0CaZclxvoPuL8SK49i+AIzl2DqAgmXBMVCDZcUxkGY5NvYaLM+ZGphk4haly+VyuRqpBr6Uulwul6tvNfDboOclrV+x02LoePTcNABgeMDsmbcxXMZs4bA1EjLOYlmUvG9pmdAcILaJroIts8NjyyaKYAImsnUGKSNt/e12XBg/r+AdK2DIq6mHcfLcIM4yq9gMxfQqvsdXR66DrZ2wYvs+6uPNKVrTxxY8T2dRTlO89aP0g/JdKOX0Ztps08RwzXg+Zn42O2+0uUPnp+mGsx2isijpFlbE1s9p0Taqx6ZYVhwDaZYVx51xYFlynO9cSrEsOQaSLCuOgTTLimMgwrLgGEizrDgGarBcJ4tymuKS5Zh/KFgWHAM1WFYcA2mWFced7aMU0/jmTA38gnOL0uVyuVyNVAO/s10ul8vVtxqYZNLzC+4TeLSMP4lHynjZo7TSfeECbLjGlpZatclWz56kd3Vl7SybonPRaA5cajMeVxRZT2x78Lk2Y1sZ3/KDCTvJnxSf36aL4it+0cLVb9GE3Tvy6z46bGmW8aW6zGoJq6nzkkU8+ZWzzIKdw/f4F+gerz9jNsn0kGV1jRZezOCw9fHkp21i7pE3qJoATVIdujgf/7HXyAJqkX26cMTiLZ+2OBSK/Lg14dctvHPDY2UcngPfn6tXGBfbfsXGOXWTLQAQXKCFLfO+xofN7uPzTX6sVcan1hX3JTI5duFn7HyKZcUxkGY5ZlEqlhXHQJplyTGQZFlxDKRZVhwDmmXFMZBmWXEMpFlWHAM1WFYcA5JlxTGQZllxDKRZlhwDNVima55tNfB1xy1Kl8vlcjVSDfzOdrlcLlffauC3Qc9LYpth2cNkv/w57fRG8UlFIhe/YJmTi9fy8tikUJzgZT2a1dea1bL8jnwcVw3toV3NOvjIQerjj2hof5x/PkpZSnzBP2cFEjDylsVDxZguW06+lTkjFSti8bgVpgxFHjnzitfk48mvwaLie3zzjuftxBMWXrLExvHJ+/5r3gfZOmwpHRkxq4ntlZA5d2Tctm//rKVmTS0kq4WzvsKjvNeaNt793TLm8Qdrh58NW3Fs0e27yjoJGYexjMSKRTnSsvE/kI9/8hZrG10yXcb30eKXkmXFMdAfy8SQYllxDNRgWXAMpFlWHANplhXHgGZZcQykWVYcA2mWFcdADZYVx4BkWXEMpFlWHANplhXHQB2WfwNzpgb+Dc4tSpfL5XI1Ug18KXW5XC5X32rgt0Hvid6UyVeZYMrx68UnOZiV9eLIMqkoYVGyDXjJh87n49lo4xk8dz45ticLO2cftJaS3XMvnyOs97ec2iIWJTkqWLY8vwlcrJTLkPAE2nKycZ17vMTCoWJM62+lLEs671KahTtYWXMwTOi17TwReOqmVdbJm/Qgwj26ud01dgC4Gnat4Voqz8a66FiXlGfF5mIrZ4xivo4FtN5jGP/wEmvj7LzK4gTqPk9YU8kx0B/LCYtScQzUYJnG9iRxmmJZcgykWRYcA5plxTFQg2XBMZBmWXEM1GBZcQxIlhXHQJplxTGQZllxDNRjec7UwC84tyhdLpfL1Ug18Dvb5XK5XH2rgd8GPS+JJ2ZyBhSWip2X1ojZ2gt2Bds+nMXD9kooBULj4eyulUvIXyIbxMxBE1/wEkR+CDGPnS2Oc2JfoBzzSVpQjtffU3H0HvN5xf2MlTfheFBkn52IlkWhO8NrDYYJstNZ2XQiUsonXEusnFDsuJB9xuVm2MLh61BlhE7IeiQd90XdZ8VxZ7uKFceAZllw3DkeyXKCYyDCsuK4o2/JsuAY0CzHmE6yHLmvKZYVx0ANlhXHgGQ5dh0plhXHQJrlWDmsOizPmTyL0uVyuVyu94ca+FLqcrlcrr7VwG+Dnpe0FbaO4mdoYmbFignuIE2OxS0WznA7d5zKorzWwudWXQUA2AkrV1FZv+9Wm5h52W6zRm54If9cMKG7WLuRfriP4juKzzpZlFT25Ltj+Qm34vay7a9wZxm/KcqFcKYUT37lLDO2c3607sqiD3s222Dr9x0mQ4vvUVh/kLfv2vsRO/FXqb8tFAc7i5751nV2fWuG95exKhHDz4zHuR9XlvFMOdGb1lGsTI61Gz4Jm5C+/fCNAIDze6wMz0u0VN+C9Ta7dyOl9ZUsK46BvlgeTGVRCo6BNMuKY6AGy4pjIM2y4BjQLCuOgRosC47zPnqzrDjm7UCE5S3WVLFlBcuK486+FcuKYyDNsuIYqMfynKmBX3BuUbpcLperkWrgd7bL5XK5+lYDk0x6fsE9igfL+MrrbYrpzUO0zlzhpLSvsaaJJfbDgUppW9Oq6/MyJOPX26s8WwN7YHbO09gEANgVsXWOkffxaw99uYyHirFdQ/ZL5YqpsC8+ZeHuNR8uzmsewXyyRjiri22SYFvwfZt4jjyuaQtDFWMuFcJr8sUmvwY75yv4RNn21CGzdXCQZiZfZGFZCfkVattC8cMUHzhOPxTpZ6/Ydb650u7LV37TrjU8h5its+NFGufzlskWbLKXWuvKpvktWxySy6Wc2kWlRcL4bYnAig321EN3lfHqG8xLDCwrjoH+WA4cA5plxTGQZllxDNRgWXCcn7c3y4pjIMLytPXB1bhTLCuOgRosK46BNMuKY0CyrDgGarAsOAbSLEuOgTTLN2Du1MDXHbcoXS6Xy9VINfA72+VyuVx9q4HfBj0v6UmyKm7EjjJetO6kxcXCfVy5ewIbyzhmUYYsI84wYu3B1WX8ND4KoGoRDEbWgLt62N7x77rnqbCziWPKMvvRGsuGClZSzKqIZXIF22ViJ9mSX6f+aOLpm1O5NcIVjLlUSGzya7CPKlbOY2RLsm3D1xrc38PUtoXiA3zgYxQX2VsHH6TNrTJ86VqzYkLCHT+b/S/aM8NjZOc8i25RxZ6zLaowzn8b+BHF3yk+d1PbKMXkAj15QzfLimOgP5Zj2XJBimMgzbLkON9ZxwXLimMgzbLiGIiwLDgG0iwrjoEaLCuO80GbtlBcsiw4BiTLiuO86wTLimMgzbLiGEizPJcWZQP/BucWpcvlcrkaqQa+lLpcLperbzXw26DnJR15w2wLrhTNJWBOFaUeeMIn2x1HIivphVf/gcqsaeqbzjdVzEKdOm4pZJxZd3jE+jjEM1bHOz6BavkTaj+Esa6+j5KtEyvZwddaxuz2vUoxr41XdBe7x1z2hq2dso+DEVtykmLOOAtPeoraeF+8RDEPWuw82ZJ9h+cz/+KztG8m95VdcLZczC7hc4QERr6mN/W+6j4rjoH+WK5y0c2y4hhIsyw57owFy4pjIM2y5BjQz0xwDKRZlhwDaZYVx0ANlhVkHTsHlgXHQA2WY12kWFYcA7VYnjM18AvOLUqXy+VyNVIN/M52uVwuV99q4LdBz0ti62QgYtGFOGbhxSrequPUdsBKZLBdMDA4Q/vOdO0LwCZeahe00j5Ex6WuifvgjMpyHLFMNxHH7vFgJC73YQuEzxuzc0IcmzSbpJsOZGuM7JfwfPjZJMcTG0Kd6wv7kHNUGRsdp+5zjNN+WI4dp7YzNymWJcedMbrbFcedY1NjkhznB3bHkeeYYjk2huSzjnEzmywLjoEaLMe6SrGsOAZqsTxn8ixKl8vlcrneH2rgS6nL5XK5+lYDvw16XtLQxWZbzI9YH2Fdu+r2GbmvshIrVgxpfsUyKSwlsgsGyQ4ZivRdnoLLg7CohIa6ppi9FLu+chxsP3A81N1e7x4LS5TPy9ZCzM65uOOzc3vUJwrXukDvSs3h+fCzqYwtdi/C7jErkpW6pshlqPusnjlvz/fpftYxS1yxrDgG0ixLjjtj1unuPmJ9q+uTHAOaZcExkGY5Zu0nWY5xkWRZcQxIlgXHQA2WFce9xpzaXoPlOVMDv+DconS5XC5XI9XA72yXy+Vy9a0GJpn0/II79pqtM3ds3CaK8rp2oSwIT47l0hzxNfC6s7e4Iq46x+s0HvDk2PEV8rgzxXzNIZ4cSw/xDM015XGGc/Dk2Fh2Gh9X7v8a9cfx0e52vsdHxu0eLqOdT6gJsjzxk9fk44miKsvstci+lRmrKlWPtr9GC+LR+crnw7YOjy025tBdHfuJxxz6PkWlUE5f1L0dmmXFcT603izH7DzFcuzfQoplxTGQZllxDKRZlhwDmmXBMZBmWXIMpFmOZUsmWY6lnAqWFcdAmmXFcec4Uwsu8HXUYHnO1MDXHbcoXS6Xy9VINfA72+VyuVx96wJ+G2RZNgngJPLUnJl2u33TXPTT+5IetxShpz9rpTeWUvnjRchLjnCJkWepxMgk140ghdIidcrl7H1uQx5s55Hb0Hf8/I1lvH6JVcJeP5zHm261Uj9nhuyldceQ3dPtsDhVLoez03h9wZeeLkpuPE7j5Jitnenis2X3ePtnbQxv09qIvIbfrr0fyYMtdC6OJylOWSMzT9MP30dv0fYDVvYF3ySrKTyTiwmr70bGuUt0MUlxKzIMLidyIhzwQ2trL6H+qGS7YFlxnA/DOk+xzPymyuWUHANJlhXHQJplxTGQZllyDGiWFcdAkmXJMZBmuY5d/W5ZVhwDaZYVx0CaZckxUIvludKFf925q91uH03v9s7lFqXL5XK5Gin/gnO5XC4X2gOz91+WZZ/Lsmw7/fe5zu4A/EWWZTvEtllT1m634xsfgG38XWv/2fV/VcbKoty7l6wYfhVnrck/FrbsDfXcjKWFnZqgcszBGqnYOhTfb+GHf9M6/Cz+CACwCWZfcLYc2znfwH1l/L1DhbXzGs3ijNkknAH21eLzy9R2gu2S4xQX1sjPUbXfX6fNN1FK1lHqUPVRqcbNZW/UgnicLZmycmKiWbH4Bxb+nOiWn9nUC/TD98R511C8lmI+IVk4Ff9MicD4ObPrAsuKY6BPlmnIimXJMZBkWXEMpFmWHANplhXHQIRlwTGQZllx3NmHZDm2sONssiw47uxasqw4BtIs98MxEMBotzdliR3fsWbeyOJfBn1qcKTdc5xZlo232+1DWZYtB/BNAL/VbrefmK3+g/wNzuVyuVwXVO12+1DxeQTAI6j8pjR78ixKl8vlcuHcLH4b9DpVlmWXAJjXbrdPFvFfA/D52eu93jiqGU30Fr1jpWUtLhrNrZ1juy+zHb5Dx8Usysn8483WMr19QvT9DLWxU0Huy0u3WAbY9hvycXLWGNs6nGW2/bBdE7YUJ+TJlbG14yYp3lJ8nmAr7usU80zQIoNvy6etiV2NN6lDzjMKfRxgm+gximMlhsP5YpNf+9Epip+08JtFphc7mCcn6QceJ48/aF8kju2TEvW3hey6gifFMdAny5MWSpYnKGYnKsGy4hiosnwW822fwqKUHANpliepbQvFkmXBMZBmWXEM9MlybH3Jd6rAsuAYqMGy4hhIs9wPx9zfpp57vRvNDMyeoTfUe/MKAI9kWQbkD/T/bbfbj/U+5J3J3+BcLpfLdcHUbrdfBHD9hejLv+BcLpfLhXODzfs66H1FtCwauwFnT9sL6MzM213buQxNJVbuAh8X6a9s18W/qyVEqL9zRYds5bBFyfH5M/RSHfrmsbMtye1yzH2UYObj2fmL9XFaNfKBMc2GnaNE4wjP51RkeyXu47yzIfGcJMe0HYBmOeaSKWZjfKdYFhwDaZYlxx3nkyzH/u0lWaa2FMuS414HKl1AjoH3JstzpHMDzVtt2bMoXS6Xy9VINe+d1OVyuVx9ix2Bpqj3F9zHKb7Xwo+O22TT0WIxuskNrbJt7yBNjr080lux+4LWibKJJ3qfHaNJo8GVGKXjOfOMxnb5ZstOCpNib6IZmlyShzPSpq+wk++4t5gAfpDsHs6mik30Dqd75RprO/Eg7cC1NYpMrcg9xs0053Ka5ky+FfqgtfMOch+TFKvS27MxOZZvwJ0W3iO6fYYmuR7lcfIE2WDh8IRYuocVPUfxt3sNEpUbKu6z4hiowbLgGNAsS47zDk2CZeaY+eWJ3oplyTGQZllxDERYFhwDaZYVx0ANliMl5GeV5Tut6R7anGRZcQykWe6HY6B6Q+dGMx+4LziXy+VyfSB0roFfB/43OJfL5XI1Ur3XovyWrUX54N1/WrY/gK+VcVjD7wB5NROwUhiTkbonobTIWKRczl4qM7Lt+GYAwNntZPfQLxuX3212zmfxh2X8a/gSAGDNbvJfyKnZvfrDZfyXZAGECeBc2XiwUmLEKiIfojIjTx3Kx4kvUye85t6PKb65+KQ1+zbebfU41sNKpHAl5K1v3w4AePPLNKmYp0hOUsyzLYPVwhN+D7DV9J8pVpls/Nvd71j496g5LP3IWXpPUcwTnSdEF61IzNpDcbmcH9s9VGLkJvLH/7mFgWXFMdAfy1wiR7GsOAbSLCuOgTTLimMgzbLkGNAsK46BJMuKY6AGy4pjoAbLsYxMwbLiGEizPBHpohWJgyTHQIrl9jOYs7UoX8aKWVuL8gocnrNx9qPmvZO6XC6Xq281McnELUqXy+VyNVI9Lcq78Fi58bfw+2X7fWe+UcaXHD8PAHh1zCyQHbD18NjuYWtkrKzobWV5+TeIfbiyjLcht0x2Yr0812ZsK+OKnfNwYedwkhLbHR+z8I1P2gTaHQP5+A/Dsru4ijf3zfuEcT563jLBjn2V1jXktfiKBLc7N5gn8wt4pIyvxt4yPolFZbwVtwEAvgLr46UnqAIzZ8PxL2Qh+YxtnW9S/F8o5nX5ysywO62J7Jx5/9DS4TavyJ8D3x9+Zke+dYUdOCG6aFEbJ6GxeD3I74i2UYppacQ7P233ObCsOAbSLCuOAc2y4hhIsyw5BpIsK46BNMuKYyDCsuAYSLOsOAZqsKw4BmqwrDgGFMuKY6AGyxORLloUK5YVx53toxQXLLc/PXcW5T6smjWLcg0OukXpcrlcrveGmjhNwC1Kl8vlcjVSPd/g2Prj+JK/MDsHr+cfl621shnDm7aW8SoclOcO1s7KH78ut6+41LKigk2wmCbj8iTt22D9rfkB2Tl/WXyyrcNXTGtYjlximZGb78jTpY4Om1U1QLN0hyJZlEHH5tlxT9xvVszrR6268zXjeWZZ7B6vP2NZlNNDdlywvo5RVhwsIQ1Tx2088y+2cQ4M5uN//TXKWOPFVXniLZcLCbeZJ79SlhnbObcXz4GfzVIqp7LtbrO+plauRqe4Ivb4sM6unWy1yvjswiIT8cO0AxXQXvignU/dZ8Ux0B/LbFEqlhXH+TCny1ixLDkGkiwrjoE0y4pjQLOsOAbSLCuOgTTLimOgBsuKY0CyrDgG0iwrjoE0y5JjoAbLkfJis6AmzoNr3hW5XC6Xq295FqXL5XK5XO8T9XyD49dztldkBQ1qGzhn+w4M6Log5fki1S+47xAPRCZbs+UixxYrTRKpBDIw031evv7q9dk+apyDgxbPq+w7070vx+fOd+0b7YPvS8TOKcdB48HF9PiZBM5UC3NleYItTX4dTDwnNfau/sR4o+cYFOfgsdF5Y+eT7MVYSLAc/Xch+o3zK641NrYUy4LjynmhWVYcd8aBIcVx176C5di+KZYlx/mJLVYsK44BybLiuDsW44/8HzTFsuS4c2wRludKTXyDc4vS5XK5XJ5F6XK5XC7X+0U93+B4YiPHt2+kEhFv5B9naDLjzgHbN7YW5XQxi3F8TGfL7aE1/H6I67rGwK/9o7AyJZtobCMfOVsZY3GgaaOOd45cBaCa3VWxlwY4i3LMjivKYvA4jz1DE70tCQs7b8z3uXqFTei+mheoG7ZwmmZ8hnNzH/tftBiTNL+S19EL181L9n2X4u0Un5ykH4rZq1wqhNbk23m79R2yzNjK4XG+snONHTiBLr1+dKXFrcXdOwDAbprdHMYcmRz7+jo7387N3SwrjoH+WOZno1hmjvleBKYBzbLkuGOckuWN1hQ4BtIsK447x1yyLDgG0iwrjjtjybLiGEizrDgGJMuK47y7BMsTkEqyrDgG0izTMqGzLc+idLlcLlcj1cS/wblF6XK5XK5GqucbHK9Jx5WER1dPW1xMWN1Pa+49W6NcTpg0y+VG+I+c+2F2Vij7EbMo2Ua4esBskvvv3hJ2MPEV0zzQZ1eZLaPK5QxFst6OYHnXcXufoyrQttwhaG4vzk9fAgDY9it2j1eQ58J2DpcZCc9kx4vkVTxGtiSvRcnXerHYvoXiqRfoB65ZUlg7XI37cbN4jtxs60uGya98fyq2JJ/2eYrD42tRW4ssHL4OtnC2iDaaM8xzYrdt7mZZcQz0xzJP/lYsM8extSgVy5LjfAeTYFlxDKRZVhwDEZanrSlwDKRZVhwDNVhWHPN2IMKy4BiQLCuOgRosK46BNMuK4852xfKcWpTNe4Nzi9LlcrlcjfyCc4vS5XK5XI1Uzze4qb22ztqeqywbrEVlo0PW1yRoX1j2Frezlfh2kVrFJTRYfI6QiTa1k9Z9o5Hvv+ow7WvHbb4mX19u5AhloZFbcOJamwnK2W6hxMmRSokRPbGc9ymrN0/QhXC1Xso+C78sTd1k17TvqjW02awRvkelffY8WTnP0nlfpZifbrhuzjzbRTEooxDH0S3aPkFZaBMWluvycb+0vWLncHu4ndPU9qYYAlC1cML5TlJq4Qla14/HJlhWHANplhXHgGZZcQykWVYcA2mWFcdAmmXJMaBZFhwDaZYlx0CaZcUxUINlxTFvh7E8YU2V9SVTLHMbW5TTFCuWFcdALZbnSk2cB+cWpcvlcrkaOU3ALUqXy+VyNVK9v7LplTv2B8jwrV/nD5TqN4TUeesoeo6Boj1yqtnou3Lc+eJ8sbUDVXute2ztpY3Q79J0or/+9JO+mqVi6yum9o09jpmuoNb51H3mZ55iOcZH6rx1JM8xQG19sPyuOQaSzMb41uMRHHeeLzk43V9/EtD2wzH3XWcMKZZnoj9cUDUxyaR576Qul8vl6ltN/IJzi9LlcrlcjVTPN7hKVVqaxMrxIpwEAJzF/LLtMGVkna2kPZnCRFA+F4szrsI+063Rso3LR4xVxjZVxosPFjUyXqYT0xUvm7L0Jl5HMPQ3RCW/Y6VOuMzG2Lz8uGMtWn+yRX0vRFd7nXs8jLepPb++l1rr7FxcVDhSbqOMeaLsJMVTNIkV+9AtypxsQcbhWipVl2lNvspx06KLyHkr13Fa7LN7ibUtFduh77PiGOiPZZ7QrFhWHANpliXHQJJlxTGQZllxDERYFhwDaZYVx0ANlhXHQA2WFceAZLllTXwdSZanI120InEYv+IYqMXyXMmzKF0ul8vVSHkWpcvlcrlc7xP1/MoeHzabIbbmXljD70zU1rF2VmotSp5AexCrAADTw6NlG9ssbHdUbKJg59jQq2v5vWihsnbmk60TW4uSJ/2Ga3p+ja1fiBbdYhs+whKFfI957Gy7LqhYO8XYWjYh9GyLJoTy9akq3GzrtCieItumYu2E9DJb47ByHB0WroXvT6VUCK/Jx5Nf5VqUFDOlXJk5OFFs90RsHcWy4hjoj+XUWpSKYyDNsuQ4P4lJsByzKFMsK46BCMs2dNBSm0mWFcdADZYj1eTTLCuOAcmy4DgfQoJlxXHneDgOLCuOgVosz5WamGTSvHdSl8vlcvWtJn7BuUXpcrlcrkaq5xvcQGVWZe/2wYptV/+42Dn60UBscqSamFrjl5Qwznc6nooS78gDkfs2GInDtQ4M0thi16T6jo5n9l7mo88j1sWg2N7PcGpcv2IvxmOM5RS/sXP0I3nvYpObEyzHeHrHepcsK46BGiynuOl3oH2oL5Zj/M4yy3Mlz6J0uVwuVyPlWZQul8vlcr1P1PMre/KNlsUjFFNKT5ggewB6O8dsW4SMtDORLEt1jsnjdC6yNVaNHOzaFwBuXbsjD3j+Lc/VtYokcvxcBTlmW3GWXVlmZTfdVi6LcZTiYrLs5MeoX7rHnHHG1b1DH6d2UUbXj+i8XOWYs8zCkKaojcdWqeujJsg+Z+GeO+U5JlstAB2W0+4huW+170KcQXZKbAeAPRSHEik8XK5Yw2MTLCuOgTjLocSL4hjQLMf+LaRYlhwDSZZjY0+xLDkGNMuCYyDNsuIYqMGy4hiowXJsordgWXAM1GBZcQykWVYcA7VYnis1Mcmkee+kLpfL5epbTfyCc4vS5XK5XI1Uzze4U98x6+DZB2zC53IcKeMwQZbtkO24sYxjFb2DJcLWCP8GsY9mP04cKvreThYBjXzH/dbfdfPMfrhuLLcqNty6Vx733KXmUU7Aru9ZbATQf0XvXS/mx+E71kclZmunsDBOrbN7vP2Bm8o4Vil6++HiWrdE+uCJwGqC7GvUdmKSfngcvfVtC394p+z77MJiki5TtZ3iLRRzFWM10Zsnv7Iqdk7wfr5vbSeXU382e1exrDgG0iwrjjvjwLLkGEiyrDgG0iwrjoH+KnqXHAOaZcExkGZZcgykWY5N9J5NlhXHQJplxTGQZllyDNRhea7UxDc4tyhdLpfL1chpAm5Rulwul6uR6v0GR4l1kw+0LBbZZ9VsMbMlOXuNdXYkt2g4C41/gwgZawAse4lf68m24JIek+up72JMy1dZSROe6xHLcNtfeAqceRZbi/LIebIRns+6x8nJiWxhLOzePnmLjWF4iWVRnqDss/N7LskDzsLiDCvOLMsoDtbOKV6TjwfXjygLbfcGiz9cfLKlFMucPPkG/VDcGC4Vwllo/ItlJTEu2Dk0HrLasJtsHcFyLIuyH5YDx4BmWXIMJFlWHANplhXHQJplyXHnOMM9VBzzdmiWJcdAmmXFMTC7LCuOgTTLimMgzbLkmMYDIMryHKmJ8+Cad0Uul8vl6ltN/BucW5Qul8vlaqR6v8FdaiFbOByH7DPOQqvsO2IxW3thH96Xxec7MnpFHtB8UB75vKVv9Rzb4rdfL9vODJmNtGig9zXxbzSxdfbenregjI8tK09gonuIH1M82r19dMl01xg69VLog0tpcH9cvoMnAger5TR5Lm2yUfoSHcd9h+fDVPF2HvMJylRT25d2bwZQnfxaZpmRlYNIRWTBsnrmnbFiWXHcGatzlRwDSZZjY0uxHDsuxbLkOD+JKdxDxTFvR5rll7iPFMuKY2B2WeZ+I89Gsqw45u2dcZDkGKjF8hypiW9wblG6XC6Xy7MoXS6Xy+V6t8qy7P4sy/ZkWbYvy7Lfnqt+er7BLfyMzeh8AF8r4/vwjTJWWZQ8gXayMuPRpKoq8ytyC5NlPHrDNABg57r1ZdvAoGUsbZ6/TY7z9oPfywOa1zk0aBlkd93xlJ1vzGybpTgGADhGvsD8WBYlzF54YnN+L54/fbN1yFYET5At5sHyPeb7uh47y5jX8FuwPp8U+tRDd9m52FLhtShVyQ6eHLvlVouneOfH0K17aeyXW/xpCxc+mF8LP5vX1620HdiWmhBdtCIxD40z2cLkV84wYyvnIRqbuM/MSiyLMsVyrKJ3YFlxDKRZlhwDSZYVx0CaZcUxEGFZcAykWVYcAzVYjpWeSbKsOAYky4JjoAbLE5EuWpE4DE1xDNRiea50obIosywbAPBvAXwc+VN+Jsuyr7Tb7Z29j+xfblG6XC6X60L+De6jAPa12+0XASDLsocB/HUAs/4F5xaly+VyuS6kLkN1UcFXirZZV9Zut6Mbfwf/qNz49/B7ZfuyRylVL7gg11jTiU2WkRUrlxPsnGVTdC56nzxw6VgZb8NmAMBOkK1D59oMsyjv/8EWO8m/KT7J1qm8s5KrgU9ZeOaO/PPoME+OPSP7PkzWzrbihI/iwbJt6/nbyvjEj0fLeP2K/JeVT+DRsu0X8IhtP0MW5ZD5NltxW1cfTxb3BwCOvGFZWFzqY+jifPzHXiOf8HFKT3vYwsragGEu7cep7e9beOfdZgOF58D3h5/ZNhrn1F4qyVK4QAtbZg2ND5vdx+fjydbl+pI8x5cy+dgy+7vDv1/GgWXJMZBkWXEMaJYVx0CaZckxkGZZcAykWVYcA5plxTGQZllx3NmHYllxDNRgeYs1geeEC5YVx0CaZcUxkGZZcgwkWT45vIynvc+q+P/371b/Z/ZP/zaAz1HTF9vt9hcBIMuyXwJwX7vd/vXi578J4KPtdvu3Zqv/ILcoXS6XyzWrFmXxZfbFyOZXAKyiny9HtdLhrMktSpfL5XJdSD0DYG2WZauzLJsP4DMAvjIXHfW0KI9iUblx2T8j++XPaafjxScvlcbWH7fzGnbBgX050vm1FN+Tf7y27kNlE1clXvYtGtu/tfDJYpyx4hn3U7zpFvohjJ/nX1Jh38p1jFN8d/7xV6t/tmx6AreXMWeRhcyyT5KVs+xhuo4JOi/PYb0v/3jmertBW6kPzoZjeyVkznE23dPYVMbf22nnqNywcK2UePbghj8t47+FL5VxsHYGI7YO97cHV5dx+M2R7b5VFYvexJb3s0VpmNiakpzJV7HYA8uKYyDNsuIY0CwLjoEaLAuOgTTLkmMgzbLgGNAsK46BGiwLjoE0y4pjoAbLimNAsqw4BtIsK46BNMuK4852xfIX8Pk5syj/Pv6PWbMofxf/W89xZln2AID/G/nqnH/Qbre/MFt9s9yidLlcLtcFXWy53W5/DaC5MHOknldU+SM8/Ra5hQr/hV98r3vR2taymxpbBPvljs9OvUBx0cnKjbZMUeU3sm9Z+GqN33aDeJbMzHctvjVMfeLferm/MxTzb77F8O56yObXLRiz+T7qN9/YPa68wfF8mGJVspuHrNLionX2mx4XrlRLMvFvvTxPamC97btjpRWjPHs6/+P9R8efLtsqcyLP2FvSJX9xPg/oXt2+0eZwja6eLmOeHxbEv/VyXEkyEXPUYr/18jjVfVYcA32y/HIkDhIcA2mW++EYMJYlx0CaZcExoFmOvcElWRYcA2mWY8vkpVhWHAOaZckxkGRZcQykWa4z15KXODM34vOyv9lQE5fq8r/BuVwul6uRcovS5XK5XI18g+v9Bcdzg47LEKH8Iq+BvZYsDrY7KtZIcZLjdLJBGs0I9/16x2fnuaidD+tHfE0/KU5yEa9izrYk9011G8vOaRCjYydkf6Pnprv2BddO5GuFaKfjFlEJgVMYLuPqH+fzC2CIYyvhLxq1eGYmL1YZW2H/kuNk54Sx8f2ha6qu2G/3JXj/sZXwWWqf2HGV1f0Fy4pjoAbLgmNAsyw57owFy7PJMVCDZcVxR6xYLjnuPC7FMrclWFYcA2mWFceAZkRyDCRZVhx3jkex3Hcli0jVldmUL7bscrlcLtf7RG5Rulwul+uCZlFeKPW+IlqyiDPIOMss2Dk3XKL3jWZRXpx/LGHrhEejzpGaUwdgA2VqPdKHz3MdxReF6+bMs5hFydlnxfjeWmsvxvuxpoxPUmmB4YHcMtlwzd6u4wFUC5dy9lmxT5ueDWde8ZJLg2L+EG/n4zg+tpuWhSuudXKDbT9A+746ZoO7bO2xyjEAcIauaT+upP5siaNgNZ2FFfA8QzFfxwEx5ljmGber+6w4Bvpk+WILJcux4xMszyrHQJplwTGgWVYcAzVYFhwDaZYVx7y987gQK44BzbLkuOM4xbLiGEizrDjm8wKa5Q2YOzXxb3BuUbpcLperkWreO6nL5XK5+lYT3+B6fsHxSuqLb7UJyzz5tcwyY0tiE8UxizJYIqsi29lrKZYfOnG9jWfgnNkWI2+ZbcF931HM8+UF2Fm0Yg9aG+mHcA72rdhGOk0xWzsfyz92DFkVyGdhJ+YJsm8XGWKrNpm/uvgFu8dgu4utneJeTCwxX2eC+jgcneidX8AxqjrKY9u7l8yP71B/xbXuHbTtE+ttaaEdtPz58KateV/0bHYO2PJGvCTRHlyFTvHYY9fBds523Fi0mU3E2WY8gfa2YmyAsSw5BvpjmZ+/YllwDNRgWXAMpFmWHANplgXHgGZZcQzUYFlwDKRZVhwDNVgWHAOaZcUxkGZZccxj74zVRO/Acd7em+VPyN5mR038gnOL0uVyuVyNlFuULpfL5WrkPLieX3D8Gr147S7bwFaNsij7yKKsWHEsOi5kL+2jDKOhAbMlN6zV2VsbiqubJEvmIuriOsqAk2OuY1FeYeGJNbntVM2QMsuBrZ1hnCq2276Ve0xr9bG1M1OMjfvguJ+1KHls2A0dh2u93Jom1+u+w6rpAwPaUpyM3JegswQDZ6HxdahMTC4euWjkJG3XfZf3WXHc2d5HFqVkWXAM1GBZcAzUYDk23hTLgmNAs6w4zrfbvpJlwXFnH4rlOmtRSpYVx4BkWXEMpFlWHANplmMZxXVYnis1cZqAW5Qul8vlaqSa95Xtcrlcrr7VxCSTnl9wbAdUbA1WrF1JFVrktgGxnTRY2fndayY2dtVNjesPGVeDAzrrqxqfq3z2Hmh6l5RC37X6Yw12fHZtfmfPhO9FP9ZILKNObY+2q/v5bjnubB8Q20kXhOVYF4nr58xBxbLiuDOu3W+fivWd1GAkLpveWxx37nMh1MQvOLcoXS6Xy9VIuUXpcrlcrg9eFiVXouX1HiuVi0N9Ds4m40mj3K7Wj4xV9Kb1AIeKfcbW2XgqlsIBcV4Au4r+Xop08QKNZ5O6Ps68qmFLjbycZ8ONr+ZqvlNlzBln4d5W7jFfB8dLLBws2lddbwPmPgYjllFYw2+I1vLjvl9ZY2tmVooUh2ttWRMfNyZiHgNn3HF2Gk8QDlpBRWtWVYAzcUZamEB7dsQy1nhyLJ9Dsqw4BtIsR9ZBlSwLjoEaLAuOgTTLkmMgzbLgGNAsK447Y8my4BhIs6w4BmqwPElj4GfWsjAcpzjmMQCaZcUxkGZZcQzUY3mu5FmULpfL5XK9T9S8r2yXy+Vy9a0mJpn0/IJbNkU1W9juoDhUMa6UCuE1+bidrZGXOz479aHufVcupdm4M93bO+MXIqcO2kfxJj5HGH+sRE4sLlyEsdXathnG213tde5xZdJ30T5+vVk53AdnZKkyI0N0UXzcwtbRMn6zZWv8hdMtaFnl4jFoC3blj1+vHAMA42P6XnDJlfAPK2p3kc4oi5La2Nbhc8j7LDgGarCsOO6MgwTHQA2W++AYMJYlx0Ca5YjtqlhWHAM1WBYcA2mWY+VyUiwrjgHNsuS44zjFsuK4czyKZcUxUI/luVITv+DconS5XC5XI+UWpcvlcrk+eFmUdb7+BhMTgStK3T8+Rz/3OtL3Rbq5/vkSE8+7VOzP2Uix1/6yPTYBNXYvBnufN6UYxOdmIuf7Se/t/Yyjn39AvG9sEm7oO3pN3J7gczAxEbiiOpehGKoj0fe75rhzHCmWad8Uy9F7rGLBcey8ddQXyz/pvb3fMcwmy9x3LZbnSJ5F6XK5XC7X+0TN+8p2uVwuV99qYpJJzy+4A5eOlfHqay3LiNO6Ro4VAZfm4CrG3M7WSJgAy1lmfH+vobg4375LL6dd7VV/9XU0tmst/Oj2/NOmXFbtno9SXBlz6Jsn+cYmytI+7eI4LoXBZVGmsbiMw/gPXGqVhCv3mO8V36Pi+riS8B5cXcZHsLyrDwCYX2SccRVkPu7UhI0NE9RfMY6zYyNl097b7Ti+vhWXHu7ql/vYjzXUbuMP1ghnpMUm0HLZkH3F+fh+j2K6jFs001eyrDgG0iwrjgHNsuAYqMGy4BiowbLiGEizLDgGNMuKY6AGy4JjIM2y4hiowfIE9UfPTLGsOO7sW7GsOAbSLCuO8/OmWZ4rNfELzi1Kl8vlcjVSblG6XC6Xq5FvcD2/4LZhcxmvvufPbAOv26cqet9i4ZlIRe+h1ERvslqeGcv9jJ1YX7ZV1ly8fmsZX3HHkTJe+3z++SBZPGzrtGicuIPi24vPWBVkmjTbpsm025bcmH/SfduGW8uY17ILMa9Zt/yOR8r4kg+dtxNTJeTnVuWWyNPYVLY9TWbrFHlNPClWVfTe+9wGO/Hj0HHIPqPr33atXd9VS/aUsVo7kJ8Z3xe2e4J4MuvByixlE1s7E4c+kge7bWb2kVErTT16w7Tsu2RZcQz0xfJQaqK34BhIs6w4BmqwrDgGkiwrjgHNsuI47yLBsuAYSLOsOAZqsKw4BiTLiuPO/hTLimMgzbLkGKjF8lypidME3KJ0uVwuVyPlFqXL5XK5GjkPrucV8Sv5a+vsfX/lRvJzhEV54voFZczZSfzqH8qFVNbkI3GWWRjHLhoPWwdLYSlwV9xk1g42FkPjLrgUykaKzZXBmcIxODpsFkhsjUdeRy6Mk+/bzrOWkvb6UcpULJzEq2HWyFVDFq/fuLOMp4fsOHUvuL+p42RRXkzWzmA+/tdfo/X5yO6qxM9QHC511JrObrcstJ1/zfpeXGR9saX0Q/LoKuPcuRqdmm5ZJ9PDo13bAWDyeIvGXNg5u2gHusU719E9Gu5mWXIMJFlWHAOaZcVxPuTeLCuOgRosC46BNMuK4664YFlxDKRZVhwDaZYVx0ANlhXHgGRZcQzELcoQK46BNMuSYyDNsk4unhU18W9wblG6XC6Xq5Fq3jupy+VyufpWE9/gen7B8et5ZQ21RLmNgXO279DAWSiV54usi6f6HohkN7HlIsuCsD1xWmzvOG5gpvu8lZId58gyGbB9BsQ4BwctnlfZd6Z7X47Pne/aN9oH34tBjmn8oZ22VxZgZBLU4ofcVjms+5nEnxP33d0Fj7eyb2UfcY7I2GLnk+xF4hTL0X8XvfpFDZZj5W1SLAuOK+eFZllx3DX+4t4rjrv2FSzH9k2xLDnOT0yxYDm2iKfgJfZsYhmcZRz5P2iKZclxZGyd55sreRaly+VyuVzvE7lF6XK5XK4PXhblZmwr42Xfomq936KdQlYXVQEeecte6zes3atPfqD45MmxNBpeXzJMfuVsSbY7Np+xceIv6XzFOLdQH+wA3Mo2EK0pGNyAlcspZS1iGY2Mv2Q/FPN4eQLquXn22n9yha1Pd1WRccb3+CMHKYVqwsKVS2wcm2/d1tUH6/CIZcPxPRoqrJbD47Z9x89byt3ZNy0zslKFPVz3vdZ0+d1WC53Hfxu2dvU7CquezO37r7JJwcH7j1UK5+NWjRhoO+7Px3+sdVnZNm+plY3ePN/GJllWHAP9sXyAYsGy4hiowbLgGKjBsuAYSLOsOAY0y4pjIM2y4rizD1ZgWXEM1GBZcQxIlhXHnX0rlhXHQJplxTFQh+VPYa7UxL/BuUXpcrlcrkaqee+kLpfL5epbTXyDy9rtdnzrDzLb+HlrfvXPLQ5GywZ2GTZRHFmLsrSBYmtRkk1Srq93E7Wx5cB2zr+y8AuFSxDLP1pA8a9SvCL0t5waY9lrV1B8T/5x9BcXlk1P4LYy5jIj65FPfr3lBxN2/B/RuagZSyi+L/8485A1bRv+2TI+FF2LMr8LXGKE1wB8PAwewEvfX2cnL6718s1mS34Wf1jGv4YvlfGaH7ySB3TD39g438Y5wGtR2lqEQWzl8Fp+LF7D79lidvMkbLLtIpws4wfwtTL+pR88aicpWFYcA32yTHamZFlxDKRZFhwDaZYlx0CaZcExoFlWHAM1WBYcA2mWFcdAmmXFMaBZlhwDSZYVx0CaZcVx3t6b5V/CVzLZ4Sxo2flXe3wZ9Kej8y6bs3H2I7coXS6Xy9VIuUXpcrlcLszMzKJFOT+9y4VQb4vyN8yifPI/WfPjal8SOyMbIl+hu4pXfyqqHK22vTbYORupkb0ayjL7wst6l350Q/F5GbWdopgrb7Q4DpVxyHLBAxa2qaJxNlEEf2Jtb/yxxU+SpcKVTm7YKPogS6lSuVlMZD9D27cPW/bWN8gz2k4LGobU4Ztosb+KnfMw2TnBXuN+aT1E3G3hG9fYv4BzA/k/rMUH6S6z3ceWGtmEz43l9hDbPVwF+faD37OdyWIPLKc4BtIs76JrVSxLjoEky7PJMZBmuUVtrVvpB8Gy5BhIsiw57uxDsRxZkCHFsuIY0CxLjjv7FiwrjoEaLAuOgRosr2rPmfW36O2js2ZRnhxe5haly+VyuVxzJbcoXS6Xy4Vzs2lRvkfU+wvu2xbWsXPEYZiM+Csv6eZShykOVYwrpULIttoyC3YO6/sdn730Q4p/+cn8cwXbhDTxNiNrp8wso5v1KNmS+6C1oDjuGu6Dn2LCohyi7Ztu3VHG0yOjZcxr550tzPRNeLpsW7Ob7Bx+2CFmS/ENiunfz8gRWtcxjJ+tHM5O5OughLQNt+YTr5evMloWv02QqLFhdll+JxwDaZZnk+POWElxDGiWJcdAXywHjoEaLEcsyhTLimMgwnKElRTLkmMgzbLgGKjB8t/AnKmJX3BuUbpcLperkXKL0uVyuVyY+Unz3uB6f8HNwtdfrGJFP8eVMVcwPq33nU3x5deZLH5ROIDHySe5WMSRijWxcQz2bER1/T0lYvjMkL3A8yoGbOeE9kopDe6DY3UBA5Htapy8fSASi/44Q+7MkI19aDBiH71DvRPOJMfAe45lyTGgWVYcd3SibncM2blimZnmWLKc4rijv3Kf2BhTLEf6q8XyHOn8uea977hF6XK5XK73jLIs+ydZlr2aZdlE8d8D6aO0mveV7XK5XK7+9d5KMvm9drv9u+/2JL2/4GjC5/1U/eKxxEmpGgWuIwtjhifFFnEsW5AnyLZuKYKN1Ejn4rI3T1L2Ek9o7UdqojdP7maLp0XxknC/eKIsxTO0xuFgyK580dp+jkqvLCXbipfwW7tR9MGzkWNZlMX5ztBs2x1DNvN4O3QcrB3OSFu/2tYfXPcxyiMMz4H73UgxjfnEtWaKhT6WTVFJJrovFdHSf89d2j3Re9GArd931x1Pyb4DyymOAc2y4hjQLEuOgSTLs8kxkGa5RW1LIvyGWHIMJFmWHHf2oVgWHANplnl9SrYoFcuS486+xZgVx0ANlgXHQE2W50rvrS+4WZFblC6Xy+UCZrJZ+y/Lss9lWbad/vtcn6P5O1mWPZdl2R9kWbY4vbuWW5Qul8vlmlW12+0vAvhibHuWZY8DWCk2/Q6AfwfgnwJoF5//EsD/+E7G0fsLjorHbnre4pnvWny8+LyODmttpB8i5XI2FZMfN8XK5fAJg21xI7VFqnH/KpUZ+feRUyvdz2MLY36H5XL4vj15qQ16GqNlvP723BpZ/ZaV1RixAr64d4LOq8qMUB8/WnNlGR/CWBkPVcqM5BfA1ZPZiqysRXnYxnz+TJ7uNX2FjZ3Lf4x/0uqijFxS9BexdZ5ddU0Z78HV6NT42CEZsw6QhTNRLA5YsXWoxMjAmD202z9l61IGlhXHQH8sB46BCMuKYyDJ8qxyDLzjcjmKZcUxUIPlSLmcFMuKYyDNsuIY0CxLjoEky4pjIM2y4hhIs3y77G2WNBurC9RUu92+N70XkGXZfwTw1Xfaj7/BuVwul+uCfsH1UpZlY+12O/y29EkAz/fav5f8C87lcrlc7yX9iyzLNiK3KCcB/O13eqKeX3BnyFIZokynW49Y/JOiFPJF5j7FqyCzNRLsnFWR0fD5ivfyM1SuYoB+2xikeMWExTcU67fF1uHjLLOKnaMqescyqzhrsbhHz15qg+eSHVwJ+W0M54ffYZWmh9ji4vX+2NopxrZ7zYfLJs4WO0KDZjtnfnEBsSrI3ztED20LzUItrnXHvWYBrR83W2r9gMWbi0wvfjY7RyxDjPvbB7OiwuTWivVJMV+HqoS8H2vKNi4xspTqdG+6wyzKwLLiGKjBsuIY0CwLjoE0y4pjIM2y5BhIsyw4BjTLimOgBsuCYyDNsuIYqMGy4BjQLCuOgTTLimMgzXKsoneK5aZYlL3Ubrf/5mydy9/gXC6Xy1WdP9IQ+TQBl8vlcjVSPd/gjg5bltJly8nDIbvjoqHutkrpXo5ZIXuL7RKeZ8iWyYru8QzQ+/TK5VRWgsYRJreyrcMXzJNfK+MPMY+dMyfPiX1pzIfpwCMUn8Qias8PrHOPKVmsbGd7hrPJjlI8SAMN1ghv5+PwGtk5r1F/4boP2vaj43YcX2u4Fn42sbEdEWAMEQxsRfF1qHNwG0+2rfSt7rPiuKM9yTJzoVgWHHeOR7IsOAZqsKw47uhbsiw4BjTLimOgBsuCYyDNsuKYt3ceV7KsOAYky4rjvL/eLCuOgTTLsX8LdVieM51L7/J+k1uULpfL5XrP/A1uNuUWpcvlcrkaqZ5vcPyaXXl9VRV2Y1mGsfhcYvvp7na2C9i2qPRNzWoNP+6i8jdVvr5wvpgtqa6D9q/aE3qSamiv3GO+Do6FDca2B583ZudYf3o8yRI/C3izzmoLzyc2Bu6b70tqX9ZgH9fEfUiWI5WikyzHnr/aV3CcjzPBcoLjzu5KlhXHneNQY44sZKBYrnWP1f2M2LkpluP9JViOlf0RLCuOeXtnH6FvxXFsX5biuHPf2HXPmRr4BucWpcvlcrka+QXnFqXL5XK5Gqmeb3CV1/6YnaNsm9ORmKWOYwmrpfJafy6yXhz1l5rWUelaWUkx+yliO4X9Byq2xkwk7rZfknYutQ9WLBBt8VTt0RmxnR5kzM453b2dzzsorJiBc9THgLZcBsSDj9m5LHWOgYiNFIvLruvY6imWY8cFRSzDJMt9cFzpOsFmV/vp3vvq+93NcWcs70tkDCmWFced/UmWI1XTFcuKYyDNsuK4e8zdLMf/H1GD37lSA9/g3KJ0uVwuVyO/4NyidLlcLlcj1fMN7jDNxlw8ThVvefLqJaKNYy7DkcqiZNE52qvCeGxC5AAtEjfCY6P+WsXnD+m0lEBVqWJcGWfoO1ZiJLIW5ZnERG8uM7KsWFsueo9pncTKGn6ij1is7B6eMDrFg3+F+pik+Cfd2w+tt+O4v0PF+fjZcPkenhSsJsgOJizMzv7COY6ct/O+PW8BbV9Ox4n7rDjubFcs17H+xPGB43w8CZYFx0ANlhXHQJplwXHnOMP9Vhzn+yZYFhx39qHimAWfZHmS+mOfV7CsOAbSLMcmeqdYjv0/og7Lc6YGvsG5RelyuVwuX4vS5XK5XK73i3q+wW2juhnr7ibLgZZ+LB0KLtPxMQtPrLHXbM5IGnm5yDKiisi8FmWbyoxsW5KX6diJ9bQr+SzXWrjuHhtnq/BzfvlJ234RXfESKgtSqWIc2tniiUyUZTtn2/DP5p/YbG0Us7UT4uXkRS672+yeZcvftBNT6Zzvjm3sOu9W3FbGFeuLPIeQGcZr3b309Do7MdfM3UKxyIB76lrr++rxPegUP5udVC+GS5rspUrI587nD35snpUVWVUBwzSJ1WW868WNefB8VrYds2UN8cRmq4h8I3aUccmy4hhIsiw5BiTLimMgzbLiGKjBsuIYSLKsOAY0y4pjoAbLguPOPhTLimOgBstbrClmKweWFcd5371ZVhwDaZYlx0AtludMvhaly+VyuRqpBv4Nzi1Kl8vlcjVSPd/gHsWDZbx69YEyvushq3gbrJ231tp35Y4hq5h7gHLABgfsHXh8df4KP7baXuW5Iu5+qpQbLIyYrcPZVMt+8UsWH86tkRVsz/DkT7ZwPmVhqGKcykjs3CeMk+/b/m0/Yyc+auH+m/LKvQPjOtvq6mv3lvFJLCzjrUVNX+5jYuctdmLOhlSTt7mEyOMUP0zxiRfoh2JMr5DXtszKjTz62zaOY/NCuRy2deyZ7X1ug51joruLYy0rDPP8Gip5zdpNF/Wd4nMXbR+18PnTN9s47+5mWXEMpFlWHAOaZcUxkGZZcQzUYFlwDKRZVhwDEZYFx0CaZcVxZx+S5dgiBCmWFceAZFlxnF9HguUJ3UWSZcUxkGb57u5TzZoa+AbnFqXL5XK5GvkF5xaly+VyuRqpnm9wW89bRtPGeRNlvGDMCniMjp0AAOyHWRXPYmMZc7YQ2xbjmCo+2daxLKR9FWsn9192njVrYXDQrINz8+y4K7GvjD/1wNfzgCfx8hWTRfnkpZbhth15zBMwY+vhHRHWTsWW5OzEH1NcVB3e+ut2j0fnTcvzcgXlv8KdAICJ58jK+Tqd91WK+7EoT3CtaD5hca0nzMLBV28ow2PXmhXzxP35tfCzOfYM1aP+Bp2WZyyH3VvU1qLB83XspjhYO3yuSym224atd3azrDgG0iwrjvO4m2XFMZBmWXIMJFlWHANplhXHQIRlwTGQZllxDNRguV+LsmRZcAxIlhXHQA2WFcdAmmXFcef5FMtuUfYltyhdLpfL1cgvOLcoXS6Xy9VI9XyDO/Hj0TKeXkExp/cUYvshNhGUrZ3hok7xMN6WfU9jcdc5Xj9qbfNojbiTKxbJ49rFxNKMJpiyxTFjyVIdY87PwdcUK28hr5WyzCp2Dk8mLvbhe1y9DmqncbwZ4unuc3XF/HSHEvviOMXqV7nDFvI10TnC8+FnU7nmaej2sPtCahsVQ+jor4x5uJGxKZYVx0CaZcVxHnezrDgG0iwrjoE0y4pjIM1y7N+sZFlwDKRZlhznO+j+Qqw4ju0LwFiOvZIIlgXHQA2WFcdAmuXY2GuwPGdq4BucW5Qul8vl8rUoXS6Xy+V6v6jnG9z6FTstho5Hz00DAIYHzJ55G8NlzBYOWyMh4yyWRcn7lpYJTXJlm+gq2DpyPLZsoggmYCJbZ5Ay0tbfbseF8XOJiliFXi4XEsbJk185y6xiMxTzh/keXx25DrZ2QkmSfdTHm1O0aJ1V9EhnUU5TvPWj9IPyXSjl9GbabPOgcc14PmZ+NjtvtMmx56fphtO6ozKLkm5hRWz9nBZto3psimXFMZBmWXHcGQeWJcf5zqUUy5JjIMmy4hhIs6w4BiIsC46BNMuKY6AGy3WyKKcpLlmO+YeCZcExUINlxTGQZllx3Nk+SjGNb87ka1G6XC6Xq5Fq4N/g3KJ0uVwuVyPV8w3uE3i0jD+JR8p42aNUyqVwATZcY2snrtpk5SEm6V1dWTvLpuhcNJoDl9qMxxVF1hPbHnyuzdhWxrf8YMJO8ifF57fpoviKX7Rw9Vs0YfeO/LqPDluaZXwtSrNaQrkQXpOPJ79yllmwc/ge/wLd4/VnzCaZHrKsrtHCixkctj6e/LRNzD3yBpXLoUmqQxfn4z/2GllALbJPF45YvOXTFoc/PH/cmvDrFt654bEyDs+B78/VK4yLbb9i45y6yRYACL85LmyZ9zU+bHYfn2/yY60yPrWuuC+RybELP2PnUywrjoE0yzGLUrGsOAbSLEuOgSTLimMgzbLiGNAsK46BNMuKYyDNsuIYqMGy4hiQLCuOgTTLimMgzbLkGKjBMl3zbKuBb3BuUbpcLperkV9wblG6XC6Xq5Hq+QbHNsOyh8l++XPa6Y3ik6ogL37BMicXr+X6D6RQfedlPZrV15rVsvyOfBxXDe2hXc06+MhB6uOPaGh/nH8+SllKfME/ZxWAMPKWxUPFmC5bTr6VOSOV33QWj1vl5VDFmDOveE0+nvwaLCq+xzfveN5OPGHhJUtsHJ+877/mfZCtw5bSkRGzmtheCZlzR8Zt+/bPWmrW1EKyWjjrKzzKe61p493fLWMef7B2+NmwFccW3b6rrJOQcRjLSKxYlCMtG/8D+fgnb7G20SXTZXwfLX4pWVYcA/2xTAwplhXHQA2WBcdAmmXFMZBmWXEMaJYVx0CaZcUxkGZZcQzUYFlxDEiWFcdAmmXFMZBmWXEM1GH5NzBnauA8OLcoXS6Xy9XIaQJuUbpcLperkeo90Zsy+SoTTDl+vfgkB7OyXhxZJhUlLEq2AS/50Pl8PBttPIPnzifH9mRh5+yD1lKye+7lc4T1/pZTW8SiJEcFy5bnN4GrcXMZEp5AW042rnOPl1g4VIxp/a2UZUnnXUqzcNlesQm9tp0nAk/dtMo6eZMeRLhHN7e7xg4AV8OuNVxL5dlYFx3rknb/ushWzhjFfB0LaL3HMP7hJdbG2XmVxQnUfZ6wppJjoD+WExal4hiowTKN7UniNMWy5BhIsyw4BjTLimOgBsuCYyDNsuIYqMGy4hiQLCuOgTTLimMgzbLiGKjH8pypgUkmblG6XC6Xq5FfcG5Rulwul6uR6vkGxxMzOQMKS8XOS2vE/BtCsCvY9uF13dheCaVAaDyc3bVyCflLZIOYOWjiC16CyA8h5rGzxXFO7AuUYz5JC8rFyt6EOHqP+bzifsbKm3A8KLLPTkTLotCd4bUGQ2bVdFY2nYhcU7iWWDmh2HEh+4zLzbCFw9eh7ucJWY+k476o+6w47mxXseIY0CwLjjvHI1lOcAxEWFYcd/QtWRYcA1WWw/1WHAM1WI7c1xTLimOgBsuKY0CyHLuOFMuKYyDNcuz/C3VYnjN5FqXL5XK5GinPonS5XC6X6/2hnm9wW2HrKH6GJmZWrJjgDtLkWNxi4Qy3c8epLMprLXxu1VUAgJ2wchWV9ftutYmZl+02a+SGF/LPBRO6i7Ub6Yf7KL6j+KyTRUllT747lp9wC+4q256ge/imKBfCmVI8+ZWzzNjO+dG6KwFUn8022Pp9h8nQ4nsU1h/k7bv2fsRO/FXqbwvFwc6iZ7513e1lvGZ4fxmrEjH8zHic+3FlGc+UE71pHcXK5Fi74ZOwCenbD98IADi/x8rwvERL9S1Yb7N7N1JaX8my4hjoi+XBVBal4BhIs6w4BmqwrDgG0iwLjgFgK+xZ/xXuBKA5BmqwLDjO++jNsuKYtwMRlrdYU8WWFSwrjjv7ViwrjoE0y4pjoB7Lc6YGJpm4RelyuVyuRn7BuUXpcrlcrkaq5xvco3iwjK+83qaY3jxE68wVTkr7GmuaWGI/HKiUtjWtuj4vQzJ+vb3KszWwB2bnPI1NAIBdEVvnGHkfv/bQl8t4qBjbNWS/VK6YCvviUxbuXvPh4rzmEcwna4SzutgmCbbF1/BA2TbxHHlc0xaGKsZcKoTX5ItNfg12zlfwibLtqUNm6+AgzUy+yMKyEvIr1LaF4ocpPnCcfihSq16x63xzpd2Xr/ymMRKeQ8zW2fEijfN5y2QLvzm+1FpXNs1v2eKQXC7l1C4qLRLGb0sEVmywpx4yq3j1DeYlBpYVx0B/LAeOAc2y4hhIs6w4BmqwLDjOz9ubZcUxUP1/QMnytPXB1bhTLCuOgRosK46BNMuKY0CyrDgGarAsOAbSLEuOgTTLN2Du5FmULpfL5WqkPIvS5XK5XK73h3q+wT1JVsWN2FHGi9adtLhYuI8rd09gYxnHLMqQZcQZRqw9uLqMn8ZHAVQtgsHIrxtXD9s7/l33PBV2NnFMWWY/WmPZUMFKilkVsUyuYLtM7CRb8uvUH008fXMqt0a4gjGXColNfg32UcXKeYxsSbZt+FqD+3uY2rZQfIAPfIziInvr4IO0uVWGL11rVkxIuONns/9Fe2Z4jOycZ9EtqthztkUVxnkBgB9R/J3icze1jVJMLtCTN3SzrDgG+mM5li0XpDgG0ixLjvOddVywrDgG0iwrjoEIy4JjIM2y4hiowbLiOB+0aQvFJcuCY0CyrDjOu06wrDgG0iwrjoE0y3NpUTYwycQtSpfL5XJ98L7gjrxhv9VxIU1eIf9UsRI2z4fh3waPRBYaCr8ZDUTu6hE631QxSWfquP2FnRMPDo9YH4d4Qs94xydQXR2e2g9hrKvvo/Rbb2xFc77WMuaXoVcp5qWDiu5i95irAvBvvmUfByNvbZMU8x/kw5OeojbeFy9RzIMWO0+2ZN/h+cy/+Cztm8l9ZRecTDAgtneeI+R38DW9qfdV91lxDPTHcpWLbpYVx0CaZclxZyxYVhwDaZYlx4B+ZoJjIM2y5BhIs6w4BmqwrCDr2DmwLDgGarAc6yLFsuIYqMXynKmBSSb+NziXy+VyNVL+BedyuVyuPItytv57F8qy7JeyLPtRlmXnsyy7qWPb/5pl2b4sy/ZkWXZf7BxBPS1Ktk4GIhZdiGMWXqwgoDpObQdsBXG2CwYGZ2jfma59AZinHPOWqX2IjktdE/fBCSflOGKJACKO3ePBSFzuwxYInzdm54Q4Nqco+edYOpCtMbJfwvPhZ5McT2wIda4v7EPOUWVsdJy6zzFO+2E5dpzaztykWJYcd8boblccd45NjUlynB/YHUeeY4rl2BiSzzrGzWyyLDgGarAc6yrFsuIYqMXynOm98ze455HP6PwP3Jhl2XoAnwHwM8hN+cezLLuq3W5Hv1L9Dc7lcrlc7xm12+1d7XZ7j9j01wE83G63z7Tb7QPIC9x/VOxXyrMoXS6XyzWrb3BZln0OwOeo6YvtdvuL7/K0lwH4Lv38StEWVc8vuKGLzbaYH7E+wrI/1e0zcl9lJVasGNL8imVSWEpkFwySHTIU6bs8Ba+ezqIVxtU1xeyl2PWV42D7geOh7vZ691hYonxeztKK2TkXd3x2bo/6ROFaF+hdqTk8H342lbHF7kXYPWZFslLXFLkMdZ/VM+ft+T7dzzpmiSuWFcdAmmXJcWfMOt3dR6xvdX2SY0CzLDgG0izHrP0kyzEukiwrjgHJsuAYqMGy4rjXmFPba7A8Z5rFLMriyyz6hZZl2eMAVopNv9Nut/9b7DDVVa9x+Bucy+VyuS6o2u32ve/gsFcArKKfLwciK4UU8r/BuVwul+s9k0XZQ18B8Jksy4ayLFuNvHLj93od0PMN7thrtgzPsXGbKMrL/oRV03lyLK9cHl8iqDt7iwsGqnO8TuMBT44dXyGPO1PM1xziybFkM5yhuaY8znAOnhwby07j48r9X6P+OD7a3c73+Mi43cNltPMJNUGWJ37ykkU8UVRlmb0W2bcyY1WZ8bT9NVoviM5XPh+2dXhssTGH7urYTzzm0Pcp8lZOX9S9HZplxXE+tN4sx+w8xXLs30KKZcUxkGZZcQykWZYcA5plwTGQZllyDKRZjmVLJlmO/VFJsKw4BtIsK447x5lacIGvowbLc6b3SBZllmWfBPCvAVwK4L9nWTbRbrfva7fbP8qy7P8DsBP5aP+XXhmUgFuULpfL5XoPqd1uPwLgkci2LwD4Qt1z+Recy+Vyud4zb3Czqd5fcI9bitDTn7WVyZdSdchFyFdk5xXYn6UV2Cd5WW1SWHm9TjWBvc9tyIPtPHIb+o6fv7GM1y+xQqHrh/N4061WCeHMkP3ZcceQTZLfDotT1QQ4O43XF3zp6WJF8sdpnByztTNdfLbsHm//rI3hbVobkdfw27X3I3mwhc7F8STFKWtk5mn64fvoLdp+gKaefJOspvBMLiasOKmXx7lLdDFJcSsyDF5t/UQ44IfW1l5C/VFFW8Gy4jgfhnWeYpn5TVUTKDkGkiwrjoE0y4pjIM2y5BjQLCuOgSTLkmMgzXIdu/rdsqw4BtIsK46BNMuSY6AWy3MlX4vS5XK5XK73h9yidLlcLlcjK3r3/oJ72MLv3WRVAAfW251QFuXevWTF8Ks46ZU1awAAC1vmd5ybsbSwUxNUrTJYIxVbx8Kzb1pBwcd/854yDpbR9Mio9UGpZ2znfAO2buf3DhXWzms0izNmk3AG2FeLT7pvOMF2yXELtxbWyEIb+9RCs8CmbqLpHkepQ9VHpVgpl71RC+JxtmTKyonpP1v4X/6BxW+JbvmZTb1AP4js3qk1FK+lDXxCsnAq/pnQFB33sNl1gWXFMdAfy4FjQLMsOQaSLCuOgTTLkmMgzbLiGIiwLDgG0iwrjjv7kCzHFnacRZYVx51dS5YjWepJlvvgGCCWN/Xc7V2pgX+Dc4vS5XK5XI2UW5Qul8vlauQbXO8vuC0U01v0jpWWtbhoNLd2ju2mNS+/Q8dFLMqQZfRma5nePiH6foba2Kkg9+WlWywDbPsN+Tg5a4xtHc4y237YrglbihPy5MrY2nGTFG8pPk+wFfd1ipmgIoNvy6etiVwNvEkdctZa6OMA2UR4jOJYieFwvtmg+BTFT1r4zSLTi9b1w8lJ+oHHyeMP2heJY/ukRP1tIWun4ElxDPTJ8qSFkuUJitmJSrCsOAbSLEuOgTTLk9S2hWLJsuAYSLOsOAb6ZDm2vuQ7VWBZcAzUYFlxDKRZ7odj7m8OLUrPonS5XC6X6/0htyhdLpfL9QHMouRXVnIDzp4262Nm5u2u7VyGphIrd4GPi/RXtsceAJcQof7OFR2exXxqG5Dx+TNk54S+eexsS3K7HHMfJZj5eHb+Yn2cVo18YExzZbDTOMLzORXZ3pcHMst+iXhOkmPaDkCzHHPJFLMxvlMsC46BNMuS447zSZZj//aSLFNbimXJca8DlS4gx8B7k+W5UgP/BucWpcvlcrkaKbcoXS6Xy9XIN7jeX3Afp5jK03103NZ9Gy0Wo5vc0Crb9g7S5NjLI70Vuy9onSibeKL32TGbNFre+FE6njPPaGyXb7bspJuKmZmbYOPlkjyckTZ9hZ18x73FBPCDZPdwNlVsonc43SvXWNuJB2kHrq1RZGpF7jFupkK101TINkxCfYXWzjvIfUxSrEpvz8bkWL4Bd1oY5iZzt8/QJNejPE6eIBssHJ4QS/ewouco/navQaJyQ8V9VhwDNVgWHAOaZclx3qFJsKw4BtIsS46BNMuKYyDCsuAYSLOsOAZqsBwpIT+rLN9pTffQ5iTLimMgzXI/HAPVGzpHep84qf3ILUqXy+VyNVJuUbpcLperkVmUWbvdjm/8FsqND979p2X7A/haGYc1/A6QVzMBK4UxGal7EtbXG4uUy9lLZUa2Hd8MADi7newe+mq+/G6zcz6LPyzjX8OXAABrdpP/Qk7N7tUfLuO/JAsgTJrlysaDlRIjVhH5EJUZeepQPk58mTrhNfd+TPHNxeevW9PGu60ex3pYiRSuhLz17XwdxTe/TJOKeW7sJMU0jNJq4Qm/B9hqovUlZSYb/y70Oxb+PWq+v/jkLL2nKOaJzhOii1YkZu2huFzOj+0eKjFyE/nj/9zCwLLiGOiPZV4nUrGsOAbSLCuOgTTLimMgzbLkGNAsK46BJMuKY6AGy4pjoAbLsYxMwbLiGEizPBHpohWJgyTHQIrl9jPIMEfKMsS/DPpUuz134+xHblG6XC6Xq5HyLziXy+VyNVI9Lcq78Fi58bfw+2X7fWe+UcaXHD8PAHh1zCyQHbD18NjuYWtkrKzobWV5eeLqPlxZxtuQWyY7sV6eazO2lXHFznm4sHM4SYntjo9Z+MYnbQLtjoF8/Idh2V1cxZv75n3COB89b5lgx75K6xryWnxFgtudG8yT+QU8UsZXY28Zn8SiMt6K2wAAX4H18dITVIGZs+EGKA7JZ2zrfJPi/0Ixr8tXplbdaU1k58z7h5YOt3lF/hz4/vAzO/KtK+zACdFFi9o4CY3F60F+R7SNUkxLI975abvPgWXFMZBmWXEMaJYVx0CaZckxkGRZcQykWVYcAxGWBcdAmmXFMVCDZcUxUINlxTGgWFYcAzVYnoh00aJYsaw47mwfpbhguf1ptyj7kb/BuVwul6uR8i84l8vlcjVSPacJsPXH8SV/YXYOXs8/LltrZTOGN20t41U4KM8drJ2VP35dbl9xqWVFBZtgMU3G5Unat8H6W/MDsnP+svhkW4evmNawHLnEMiM335GnSx0dNqtqgGbpDkWyKIOOzbPjnrjfrJjXj1p152vG88yy2D1ef8ayKKeH7LhgfR2jrDhYQhqmjtt45l9s4xwYzMf/+muUsTZIN4Mn3nK5kHCbefIrZZmxnXN78Rz42Sylcirb7jbra2qlVXwO4orY48M6u3ay1Srjs6GC9IdpByqgvfBBO5+6z4pjoD+W2aJULCuO82FOl7FiWXIMJFlWHANplhXHgGZZcQykWVYcA2mWFcdADZYVx4BkWXEMpFlWHANpliXHQA2WI+XFZkWzOdP7ovQuF0D+BudyuVyuRsonertcLpcLs7sY5XvjDa7nFxy/nrO9IitoUNvAOdt3YEBPjy/PF7mn3HeIByKTrdlykWOLzdCPVAIZmOk+L19/9fpsHzXOwUGL51X2nenel+Nz57v2jfbB9yVi55TjoPHgYnr8TAJnqoW5sswrTX4dTDwnNfau/sR4o+cYFOfgsdF5Y+eT7MVYSLAc/Xch+o3FkuXY2FIsC44r54VmWXHcGQeGFMdd+wqWY/umWJYc5ye2WLGsOAYky4rj7liMP/J/0BTLkuPOsUVYnjvNpkW5IL3LBZBblC6Xy+VqpNyidLlcLheaWC+n5xccT2zk+PaNVCLijfzjDE1m3Dlg+8bWopwuZjGOj+lsuT20ht8PcV3XGPi1fxRWpmQTjW3kI2crYywONG3U8c6RqwBUs7sqltIAZ1GO2XFFWQwe57FnaKK3JWFh5435PlevsAndV/MCdcMWTtOMz3Bu7mP/ixZjkuZX8jp64bp5yb7vUryd4pOT9ENhW3CpEFqTb+ft1nfIMmMrh8f5ys41duAEuvT60ZUWtxZ37wAAu2l2cxhzZHLs6+vsfDs3d7OsOAb6Y5mfjWKZOY79e1IsS447xilZ3mhNgWMgzbLiuHOcJcuCYyDNsuK4M5YsK46BNMuKY0CyrDjOu0uwPAGpJMuKYyDNMi0TOvtqXr0ctyhdLpfL1Ui5RelyuVwuNPENrucXHK9Jx5WER1dPW1xMWN1Pa+49W6NcTpg0y+VGuELxfpidFcp+xGwdthGuHjCb5P67t4QdTHzFNA/02VVmy6hyOUOVTDfr+wiWdx239zmqAm3LHYLm9uL89CUAgG2/Yvd4BXkubOdwmZHwTHa8SF7FY2RL8lqUfK0Xi+1bKJ56gX7gmiUF9FyN+3GzeI7cbOtLhsmvfH8qtiSf9nmKw+NrUVuLLBy+DrZwtog2mjPMc2K3be5mWXEM9McyT/5WLDPHsbUoFcuS43wHk2BZcQykWVYcAxGWp60pcAykWVYcAzVYVhzzdiDCsuAYkCwrjoEaLCuOgTTLiuPOdsXynFqUzfsbnFuULpfL5Wqk3KJ0uVwuFz5wFuXUXltnbc9Vlg3WorLRIetrErQvLHuL29lKfLtIreISGiw+R8hEm9pJ677RyPdfdZj2teM2X5OvLzdyhLLQyC04ca1NRuRst1Di5EilxIieWM77lNWbJ+hCuFovZZ8Fq2nqJrumfVetoc1mjfA9Ku2z58nKeZbO+yrF/HTDdXPm2S6KQRmFOI5u0fYJykKbsLBcl4/7pe0VO4fbw+2cprY3xRCAqoUTzneSUgtP0Lp+PDbBsuIYSLOsOAY0y4pjIM2y4hhIs6w4BtIsS44BzbLgGEizLDkG0iwrjoEaLCuOeTuM5QlrqqwvmWKZ29jhm6ZYsaw4BmqxPHdyi9LlcrlcrveF3KJ0uVwuFz5wFiW/sZ6rpG9x+2DP7Wrfalvv89ZR9BwDRXvkVLPRd+W488X5YmsHqvZa99jay0zTft0E0V9/isDfz7+J2PqKqX1jj2OmK6h1PnWf+ZmnWI7xkTpvHclzDFBbHyy/a46BJLMxvvV4BMed50sOTvfXnwS0/f6/Xaz9mdwX0M9vJvrDBZZblC6Xy+VyvS/kFqXL5XK58IGzKCtVaWkSK8eLcBIAcBbzy7bDlJF1tpL2ZAoTQflcLM64CvtMt0bLNi4fMVYZ21QZLz5Y1Mh4mU5MV7xsytKbeB3B0N8QlfyOlTfhMhtj8/LjjrVo/ckW9b0QXe117vEw3qb2/Ppeaq2zc3FR4Ui5jTLmibKTFE/RJFbsQ7coc7IFGYdrqVRdpjX5KsdNiy4i561cx2mxz+4l1rZUbIe+z4pjoD+WeUKzYllxDKRZlhwDSZYVx0CaZcUxEGFZcAykWVYcAzVYVhwDNVhWHAOS5ZY18XUkWZ6OdNGKxGH8imOgFstzJ7coXS6Xy+V6X8gtSpfL5XLhA2dRjg+bzRBbcy+s4XcmautYOyu1FiVPoD2IVQCA6eHRso1tFrY7KjZRsHNs6NW1/F60UFk788nWia1FyZN+wzU9v8bWL0SLbrENH2GJQr7HPHa2XRdUrJ1ibC2bEHq2RRNC+fpUFW62dVoUT5FtU7F2AvS2xmHlODosXAvfn0qpEF6Tjye/yrUoKWZKuTJzcKLY7onYOoplxTHQH8uptSgVx0CaZclxfhKTYDlmUaZYVhwDEZZt6KClNpMsK46BGixHqsmnWVYcA5JlwXE+hATLiuPO8XAcWFYcA7VYnjs1z6L0NziXy+VyoYlvcP43OJfL5XI1Uj3f4AYqsyp7tw9WbLv6x8XO0Y8GYq/WamJqej56Oc53Op6KEu/IA5H7NhiJw7UODNLYYtek+o6OZ/Ze5qPPI9bFoNjez3BqXL9iL8ZjjOUUv7Fz9CN572KTmxMsx3h6x3qXLCuOgRosp7jpd6B9qC+WY/zOMstzJ7coXS6Xy9VIuUXpcrlcLtf7Qj3f4CbfaFk8QjGl9IQJsgegt3PMtkXISDsTybJU55g8TuciW2PVyMGufQHg1rU78oDn3/JcXatIIsfPVZBjthVn2ZVlVnbTbeWyGEcpLibLTn6M+qV7zBlnXN079HFqF2V0/YjOy1WOOcssDGmK2nhslbo+aoLscxbuuVOeY7LVAtBhOe0ekvtW+y7EGWSnxHYA2ENxKJHCw+WKNTw2wbLiGIizHEq8KI4BzXLs30KKZckxkGQ5NvYUy5JjQLMsOAbSLCuOgRosK46BGizHJnoLlgXHQA2WFcdAmmXFMVCL5blT897g3KJ0uVwuF5r4Nzi3KF0ul8v1nlGWZb+UZdmPsiw7n2XZTdTeyrLsVJZlE8V//z51rp5vcKe+Y9bBsw/YhM/lOFLGYYIs2yHbcWMZxyp6B0uErREup7GPZj9OHCr63k4WAY18x/3W33XzzH64biy3Kjbculce99yl5lFOwK7vWWwE0H9F710v5sfhO9ZHJWZrp7AwTq2ze7z9gfJZRitFbz9cXOuWSB88EVhNkH2N2k5M0g+Po7e+beEP75R9n11YTNJlqrZTvIVirmKsJnrz5FdWxc4J3s/3re3kcurPZu8qlhXHQJplxXFnHFiWHANJlhXHQJplxTHQX0XvkmNAsyw4BtIsS46BNMuxid6zybLiGEizrDgG0ixLjoE6LM+d3jMW5fMAPgXgP4ht+9vt9sa6J3KL0uVyuVx4r1iU7XZ7FwBkWfauz+UWpcvlcrlmVVmWfS7Lsu303+dm6dSrsyx7Nsuyb2dZdltq595vcJRYN/lAy2J6/w7WTsgwy7ebLcnZa6yzI7lFw1lovBYln6/MXuLXerItuKTH5Hrquxjn8lVW0oSrHMcy3PYXngJnnsXWojxynmyE57PucXJyIv+CtLB7++QtNobhJZZFeYKyz87vuSQPOAuLM6w4s4x/AQrWzim2IXhw/Yiy0HZvsPjDxSdbSrHMyZNv0A/FjeFSIZyFxpNfK4lxwc6h8ZDVht1k6wiWFcdAfywHjgHNsuQYSLKsOAbSLCuOgTTLkuPOcYZ7qDjm7dAsS46BNMuKY2B2WVYcA2mWFcdAmmXJMY0HQJTlOdPsWZTtdvuLAL4Y255l2eMAVopNv9Nut/9b5LApAFe02+1jWZbdCOC/Zln2M+12+43I/m5Rulwulwu4kBZlu92+9x0ccwbIVw1vt9s7sizbj3yCzPbYMW5Rulwul+s9ryzLLs2ybKCI/wfkNSBe7HVM7ze4Sy3kibBs54R21QYAi0YsZmsv7MP7svh8R0avyAOaD8ojn7f0rZ7jXPz262XbmSGzkRYNnOx5HGd1xtbZe3vegjI+tqw8gYnuIX5M8Wj39tEl011j6NRLoQ8upcH9cfkOnggcrJbT5Lm0yUbpS3Qc9x2eD1PF23nMJyhTTW1f2r0ZQHXya5llRlYOIhWRBcsxZlMsK447Y3WukmMgyXJsPCmWY8elWJYc5ycxhXuoOObtSLP8EveRYllxDMwuy9xv5NlIlhXHvL0zDpIcA7VYnjO9N7Iosyz7JIB/jZyo/55l2US73b4PwO0APp9l2QzyFVr/53a7fbzXudyidLlcLhfeQ1mUjwB4RLT/GYA/6+dcblG6XC6Xq5Hq+Qa38DM2o/M+fEPGwRKpZHrRBNrJyoxHk6qqzDZKC5NlPHrDNABg57r1ZdvAoP22sXn+tjJ+AF8r49sPfi8PaF7n0KBlkN11x1N2vjGzbZbiGADgGPkC82NZlDB74YnN+b14/vTN1iFbETxBtpgHG7vH67GzjHkNvwXr80mhTz10l52LLRVei1KV7ODJsVtutXiKd34M3aK/Cd90ucWftnDhg/m18LN5fR0lSrEtNSG6aEViHhpnsoXJr5xhxlbOQzQ2cZ8Vx0B/LMcqegeWFcdAmmXJMZBkWXEMpFlWHAMRlgXHQJplxTFQg+VY6Zkky4pjQLIsOAZqsDwR6aIVicPQFMdALZbnTu8Ni3I25Raly+VyufBesShnU25Rulwul6uRytrt9k97DC6Xy+Vyzbr8Dc7lcrlcjZR/wblcLperkfIvOJfL5XI1Uv4F53K5XK5Gyr/gXC6Xy9VI+Recy+VyuRqp/x9jn2VtLqlKOQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x504 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dim=128\n",
    "\n",
    "#Domian, i.e., x1, x2, x3 are in the range of [0,L]\n",
    "L=2*np.pi\n",
    "dx = 2*np.pi/L;\n",
    "\n",
    "x = np.linspace(0, L, dim+1)\n",
    "x=x[:-1]\n",
    "X,Y,Z=np.meshgrid(x,x,x)\n",
    "\n",
    "# Assume that there are three dominant wave-number with different amplitudes\n",
    "w1 = 2.0; Aw1 = 6.0;\n",
    "w2 = 8.0; Aw2 = 4.0;\n",
    "w3 = 16.0; Aw3 = 3.0;\n",
    "\n",
    "# Assuming isotropic and homogeneous conditions, the following artificial velocity field is considerd.\n",
    "u = (Aw1*(np.sin(w1*X)+np.sin(w1*Y)+np.sin(w1*Z))\n",
    "    + Aw2*(np.sin(w2*X)+np.sin(w2*Y)+np.sin(w2*Z)) \n",
    "    + Aw3*(np.sin(w3*X)+np.sin(w3*Y)+np.sin(w3*Z)));\n",
    "v = u;\n",
    "w = u;\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 7))\n",
    "sns.heatmap(u[:,:,0], cmap='jet', xticklabels=False, yticklabels=False, robust=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0eb71b2",
   "metadata": {},
   "source": [
    "Apply fftn to perform discrete Fourier transform in 3 dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "id": "e5091cb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ufft=np.fft.fftn(u)\n",
    "vfft=np.fft.fftn(v)\n",
    "wfft=np.fft.fftn(w)\n",
    "\n",
    "Euu=(np.abs(ufft)/dim**3)**2\n",
    "Evv=(np.abs(vfft)/dim**3)**2\n",
    "Eww=(np.abs(wfft)/dim**3)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8786fe33",
   "metadata": {},
   "source": [
    "Now it is time to define scalar $k$ as $(K_i K_i)^{1/2}$ and integrate $E_{ii}$ over a sphere of radius $k$ to obtain $E(k)$. \\\n",
    "$E(k) = \\frac{1}{2} \\oint E_{ii}(\\vec{K}) dS(k)$,\\\n",
    "where $S(k)=4\\pi k^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "id": "485d2ece",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define radius which start from zero, reaches to maximum number of dim/2 then again reduces to zero\n",
    "ri = np.concatenate((np.array(range(int(dim/2+1))), np.array(range(int(dim/2-1)))-dim/2+1))\n",
    "\n",
    "rx,ry,rz = np.meshgrid(ri, ri, ri)\n",
    "r = np.sqrt(rx**2+ry**2+rz**2)*dx\n",
    "\n",
    "# define 1D wave number\n",
    "kf = int(dim/2);\n",
    "k = (np.array(range(kf))+1)*dx\n",
    "\n",
    "# define bins to discretize k-domian which helps calculating the integral\n",
    "bins = np.zeros(kf+1)\n",
    "bins[1:-1] = (k[:-1]+k[1:])/2\n",
    "bins[-1] = k[-1]\n",
    "\n",
    "Ek = np.zeros(k.shape[0]);\n",
    "\n",
    "# find bin indices to which the r values belong\n",
    "bin_ind = np.digitize(r, bins)\n",
    "bin_samples = np.zeros(k.shape[0]);\n",
    "\n",
    "for i in range(kf):\n",
    "    Ek[i] = np.sum(Euu[bin_ind==i+1])+np.sum(Evv[bin_ind==i+1])+np.sum(Eww[bin_ind==i+1]);\n",
    "    bin_samples[i] = np.sum(bin_ind==i+1);\n",
    "Ek = Ek/bin_samples;\n",
    "\n",
    "# The for-loop calculated the avergae Ek over a sphere of radius k\n",
    "# Accordingly, it should be multiplied by the area\n",
    "Ek = Ek*2*np.pi*(k**2)/(dx**3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83383702",
   "metadata": {},
   "source": [
    "Plot the energy spectra, i.e., $E(k)$ as a function of $k$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "id": "8da77b77",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2c4035f50d0>]"
      ]
     },
     "execution_count": 202,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAE/CAYAAAB4uRluAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvYklEQVR4nO2dfbRU5X3vv1/eVMAXkAMqoIAKKZSal1OTNDZNWmqhMSVVr0L0JmlcWk0xqet2rdquGmNsbtrb3toYNQavxpgohoXeSBWjNdFLbPAF0QgEQUQCJwgcPSLvZ87M/O4fZwbm7LNfnj2z9+z9PPv3WYvFmT179vz2d579fX7P79kvFBEoiqIo7WdI1gEoiqIUFTVgRVGUjFADVhRFyQg1YEVRlIxQA1YURckINWBFUZSMGJZ1AO1m3LhxMmXKlKzDUBTFMV566aW3RaQjzmcKZ8BTpkzB6tWrsw5DURTHIPnruJ/REoSiKEpGqAEriqJkhBqwoihKRqgBK4qiZIQasKIoSkZYb8AkP0PyLpKPkDw/63gURVFMyaUBk7yH5G6S6zzL55LcSHIzyesBQER+LCJXAvgCgEszCFdRFKUpcmnAAO4FMLdxAcmhAG4HMA/ATAALSc5sWOUfau8riqJYQS4NWERWAujxLD4XwGYR2SIiJQAPApjPfv4ZwOMisqbdsSqKojRLLg04gIkAtje87qotuxbAHAAXk7za74MkryK5muTq7u7u9CNVFEUxwKZLkemzTETkVgC3hn1QRBYDWAwAnZ2dqTyDqVIVDB3iF6KiKIo/NmXAXQAmN7yeBGBHRrEM4mCpDFANWFEUc2wy4BcBnE1yKskRABYAWJ5xTEc4WKpkHUJu6KtUga99LeswFCX35NKASS4BsArADJJdJK8QkTKARQCeALABwFIRWZ9lnI0c6C1nHUJuONhbAW66KeswFCX35LIGLCILA5avALCizeEYcbBUAS64IOswcsH+UhknqhapUqkKhr68BvjQh7IORWmBXGbANnKgtwz8x39kHUYuOKhapM6BUlnN1wHUgBPiYKkCfPrTWYeRC/b3llWLlDlUqkAuUI1tRw04IQ6WKsCjj2YdRi5QLdKnVK6Cj6nGtqMGnBClip4FUWe/TkimTl+lmnUISgKoASdEXzmV6zus5GBJDThtSmrATqAGnBClShV7DvRmHUYuONxXRbmsI4I06SsL3tpzMOswlBZRA06IvkoVQ+66K+swckFfpYrqdxdnHYbTlCpVjLj77qzDUFpEDTghyhXBCX+9KOswckGpXMWIv7om6zCcplSu4uT/cW3WYSgtogacEFqTO0pfRevhaaOTcG6gBpwAIqIHRA0RQamsWqSNtjc3UANOgHJVIAL8+ntLsg4lcypVQVUEu3+4NOtQnKZUruL1u+7POgylRdSAE6Cejeyd9TsZR5I99fLDgdnnZByJ25QqVeyZqe3NdtSAE6B+DvDsc2dlHEn21GvhU8+ZkXEkbtNXEfzux2ZnHYbSImrACaATcEfR+m970BqwG6gBJ4AeDEdRLdqDdnRuoAacAHXT6br4sowjyZ66Fjsv/e8ZR+I2fZUqtv75Z7MOQ2kRNeAEqJvOL2/4l4wjyZ56ZvbaTf+acSRuU6pU8cLf/1PWYSgtogacAKXaJNzHLv2TjCPJnno9/IMX/nHGkbhL/VzrP7r8T7MORWkRNeAEqGfAJ722NuNIsqd+GtoJv3o140jcpX7e+ckb12UditIiasAJoBNPR1Et0kc1dgc14ASoD7sPnjw+40iyp14DPjROtUiL+nnnB7S9WY8acALUh933P7wq40iyp94ZPfzoCxlH4i69taevLF7yc4jojY9sRg04AfpqWd85d/97toHkgLoWs777bxlH4i71Dv8j931bLwKyHDXgBKjX5M6979uFz0jqhvD+792acSTuUu/kPvrD2/TWn5ajBpwAjVlIuVrsA0IniNKnsb316RVxVqMGnACNWUjRDajx4aRFHw2kReNlyFqCsBs14ASom+79tz1U+Kcjlxq10OFxKjS2N70nhN2oASdAY9Zb9Iyk0RCKPhpIC9XYHaw3YJK/RfJOkstIZvIkyPoBcdmiiwp/QNT3/7JFF2l2lhL1kcVliy4qfIdvO5kaMMl7SO4muc6zfC7JjSQ3k7w+bBsiskFErgZwCYDONOMNQmvAR2nc/6JrkRYDNC54yct2ss6A7wUwt3EByaEAbgcwD8BMAAtJziQ5m+Sjnn/ja5/5MwDPAvhpe8PvR02nn2pVBnRGmp2lQ0lLXs4wLMsvF5GVJKd4Fp8LYLOIbAEAkg8CmC8i3wRwQcB2lgNYTvIxAA+kGPIgKlVBpXbq2arLF+HEAmckjWaw6vJFOEUn4VKhXtpZdfkiQMs8VpN1BuzHRADbG1531Zb5QvITJG8l+V0AKwLWuYrkapKru7u7Ew22MeN97nPXFjoDVi3aQ11X1dh+8mjA9FkWmEqJyDMi8mUR+UsRuT1gncUi0ikinR0dHYkFCgzM+q5ccF6hD4jG8sOVC87TSbiUqLexorc3F8ijAXcBmNzwehKAHRnFEknjlUije7oLXZMrqRZtoa7z6J5uNWDLyaMBvwjgbJJTSY4AsADA8oxjCsR7sUGRLz7wmoFeJpsOpYY21qsaW03Wp6EtAbAKwAySXSSvEJEygEUAngCwAcBSEVmfZZxhNJrOrrNmFdp0Sl4tCtwZpUm9zanG9pP1WRALA5avQMCEWt5oNJ0H7ngY7yvwkLDPo8X7C6xFmtQ7+QfueBinqMZWk8cShFU0ms6cW24odN2zsQZcdC3SolKVI3fcm3PLDTrRaTlqwC1SbhgCzn58aaGHhI2d0ezHl6o5pIBXY52Esxs14BbxZnlFPiBKZe+EZHG1SAtve9NRht2oAbeId9KtyKYz6CyIAmuRFt5RRalc1fsuW4wacIs0lhwWL1lZ6GF3474vXrJywOlSSjI0dmqLl6yEiD6FxWbUgFuk8YCYsGm91oBrTNi0vtCn5KVF493PJmzqPztTRxr2ogbcIo01uPk3XlPog0G1SB+vxoDektJm1IBbxGsyjXdHKxp+9UklWfw07a1UMohESQI14Bbxy/KKmvl5yy/lqqBa0M4oLfzbm2psK2rALdI4/HvqK18HUNxTgxrNoehapIWfxlprtxc14BZpNJi1n7oUQHEPiD4/LdSAE6WxBFHXWDs5e1EDbpFGg7nu/Bm1ZcUcEjbemavoWqRFyae9aa3dXtSAW0RrwP2IiO9+qzkki1+HVsT25gpqwC3id0AUcUhYrgr8LshSc0gWvw5NOzl7UQNuAREZ0Pi3fPiTAIppOt59rmtRxM4oTRp1PtretMxjK2rALeBt+I/cfGf/8gKeGO/Nwo5ooQacKI0dmmpsP2rALeBt+PNvuBpAMbM+7z7XtShiZ5QmjW2urrE+lshe1IBbwGvA055/2nd5EfCOBupaFLEzSpPGkUaR25srqAG3QJC5FPGACJoI0gmiZNGzbtxCDbgFgiY/inhABO1zEbVIEz0NzS3UgFvAe8XbLU9uBDD4yRBFwJvp1rVQc0gO71k3R9ubamwrasAt4DWX2Y/9yHd5EVAt0sdb8qprrDe+txc14BbwHhBzvvVVAMU0HW8WVtdCzSE5vOWHIrc3V1ADbgGtAR8lSAsdHidH0E2eVGN7UQNugSCjLWLWp5Nw6RN01k2RHwJgO2rALeDNSB656Tu+y4uA1xyOaKEGnBiDrjasaQyozraiBtwCXtPZNX0WgGIeDF5zqGuhw+Pk8LarusaAXg1nK9YbMMnTSS4neQ/J69v53d6651ULP15bXryDwbvPR7XQoXFSeDv8usZAMducC2RqwDXT3E1ynWf5XJIbSW42MNXpAB4TkS8CmJlasD4E1z2L9yw0rQGnT9h9NVRnO8k6A74XwNzGBSSHArgdwDz0G+pCkjNJzib5qOffeAAvA1hA8mcAnm5n8GGNvq9arAMiqNRQqQrKag6JEHZfDb3pkZ0My/LLRWQlySmexecC2CwiWwCA5IMA5ovINwFc4N0Gyb8BcGNtW8sAfC/lsI/gNZ218y458ndfRXBMpuq2F++ZH14thg1td0TuEdbeSvpoeivJOgP2YyKA7Q2vu2rLgvgJgC+TvBPAVr8VSF5FcjXJ1d3d3YkF6q1vPnXdzUffK9ikiHc00KiF3hEtGUI11gzYSvJowPRZFti6RGSdiFwsIleLyN8ErLNYRDpFpLOjoyOxQMueMsNnv3Thkb+LVpPzdjhF1iItvDqqxvaTRwPuAjC54fUkADsyiiUU75Bwwub1R98r0AFRqQrKnknHRi3UHJLBq2NR25tL5NGAXwRwNsmpJEcAWABgecYx+RJ2ilWRTr+KMlg9FzgZws71LVrJyxWyPg1tCYBVAGaQ7CJ5hYiUASwC8ASADQCWisj6sO1khdd49o/tCHzPZfyyr6JqkSbeTr1RY82A7STrsyAWBixfAWBFm8OJhd/193c9+OyRv4uU9fllXwO1KM5oIE28HVlR25tL5LEEYQV+Wd1H7vt26Puu4pd9FVWLNPGa7ECNtZOzETXgJvEznY/+8LYjfxfpgPC7CGCgFmrASeDVUTW2HzXgJoma9CjSARFVf9T6ZDKE6aglCDtRA26SqAy3SKYTdfCrObSOiIReblyk9uYSasBN4pfh3n/bQ0ffL5DpRGpRoHJMWpSrgqoM1HGgxsVpby6hBtwkURlHkUwn6uBXc2gd1dhN1ICbxK/BX7bootD3XcWvMyqqFmnhV35o1FjLPHaiBtwkUbf/K1JNTmvA6dMbcbezIt6D2gXUgJskugRRHNPRCcn0MSlpFe0e1C6gBtwkfga76vJFR98vUNanWqSPn4aNGgM60rARNeAm8TOd5z53bcP7xRkOqhbpE6Vx/zqqs22oATeJ3wFx5YLzjvxdpGG33126iqpFWkRpDBSr7OUKasBN4neDmdE9R5+20VepQqQYGYnfgV9ULdIiSmNASxA2ogbcJFHZhggG3aTcVaJqvCI6PG4VE/10pGEfasBN4mfAu86aFbmOi/gd+EXVIi38slvV2H7UgJvEr7E/cMfDA9cpyH1w/bKzQVqoObSEtjc3UQNuEu9j2AFgzi03eNZx33RExDc7K6IWaeKn32CN9dH0tqEG3CR+dc/Zjy8duE4BTCeoNjlYC83OWsGvk/NqrE8esQ814CYxMdcizEqbZrZ6MUZrmLS3InT4rqEG3CR6QPRjaqxagmgN7fDdRA24CapV8R1SL16ycsDrIphOkDEM0kLNoSX85hy8Ghehw3cNNeAmCLrpyYRN6weuV4C6Z1AnM1gLNYdW8OvAvBoXocN3DTXgJggy1vk3XuNZz/0DIiizHayF+51Rmvi1Ja/GOsqwDzXgJjCtexZh4snUWIvQGaWJSVvSTs4+1ICbwNRMijAkNNaiAJ1RWlSqYnRZu3Zy9qEG3ARBxvrUV74+4HURMhK/u3QBg7UoQmeUFkHGOkhj7eSsQw24CcoBxrr2U5cOeF2EjCRoH4uoRVoEdV5ejbWTsw814CYIMpPrzp9htJ5LqBbpE1T/9dNYb/tpF1YZMMlpJO8muaxh2W+RvJPkMpLXhH0+KUwzjSIMCU2NVW8U0zym7a1It0B1hbYZMMl7SO4muc6zfC7JjSQ3k7w+bBsiskVErvAs2yAiVwO4BEBn8pEPxnzm3/2DwbST0eFx88TpvIrQ6btEOzPgewHMbVxAciiA2wHMAzATwEKSM0nOJvmo59/4oA2T/DMAzwL4aXrhHyUo69vy4U8arecSfldoAYO1UGNonqDOy6sxUIw25xLD2vVFIrKS5BTP4nMBbBaRLQBA8kEA80XkmwAuiLHt5QCWk3wMwAMJhRxIUE3ukZvvHLheAQ4G1SJ9gjovr8aAjjRsI+sa8EQA2xted9WW+ULyZJJ3AvgAyb+rLfsEyVtJfhfAioDPXUVyNcnV3d3dfqvEIqiRz7/haqP1XCLIWL1aqAE3j6nGgI40bKNtGXAA9FkWWPASkXcAXO1Z9gyAZ8K+REQWA1gMAJ2dnS0XZoNqu9Oef3rgemWBiID02003COpkBmlREVSrgiFD3NUiLYIM2Ktx/7ruzzu4RNYZcBeAyQ2vJwHYkVEsxphmc1URVByflY6TcQXdxEgJJ5bGOtKwitgGTHJUbfIsCV4EcDbJqSRHAFgAYHlC206NOI3c9Ywkzv7p8Lg54pSyVGO7iDRgkkNIfpbkYyR3A3gNwFsk15P8F5Jnm3wRySUAVgGYQbKL5BUiUgawCMATADYAWCoi68O2kweCGvktT24cvK7DGUnQ8+AAfy1c74zSIki3orU3FzHJgJ8GcCaAvwNwiohMFpHxAH4fwHMA/onk5VEbEZGFInKqiAwXkUkicndt+QoRmS4iZ4rIN1rYl7YR+By0x37ks667B0SlKqgGXHlVNC3SJEg3X401A7YKEwOeIyI3i8irInLk1xWRHhF5SEQuAjC4JThM0AEx51tfNV7XBcIyWj8tdHjcHEG6+bc3HWXYRKQBi0gfAJD8R+979VpwfZ2iEKsG7PAluHGHuy53RmkSqwasj6a3ijiTcBNJLqy/qF2Z9lTyIeWfeAeEu6YTN6PV7Kw54nRc+mh6u4hzHvBfAniC5BvoP1f3ewD+NpWocoyIBGa1j9z0nUHLXM76wvbNTwstQTRH4JVwBWtvLhJpwCTvA7AGwMsA/gr9l/qWAXxGRDanG17+CJt42jV91qBlLh8QYfvmp4XLo4E0CdLZV2Pt5KzCpATx/dp6X0S/+U4B8C6Ay0lenF5o+SRsGH3Vwo/7rO/uARG2b0XTIk2C2pxqbD+RGbCI/BQNdxkjOQz9dy47B8BHACwL+KiTxM3iXK7JBT2OKAg1h/iEnWvth2psFyYlCErDbfZrF0+8Wvv3A791XCZuA3f5gIg7qeayFmkRt8OP2ykq2WJ0IQbJa0me3riQ5AiSf0jy+wA+n054+SPMRNbOuyTW+rYTVwutT8YnrJPzb2+FyIOcweQsiLnor/8uITkVwB4AxwIYCuBJALeIyCtpBZg3ws7rfeq6mwev77ABhxmqnxZBN29Xggm7sq1o7c1FTDLgM0XkDhH5GIAzAPwRgA+KyBkicmWRzBcIHxJ+9ksX+qzvrunE1UIvk41PXI0rVffvwOcSJgb8g4a/Py8ib4nIHgAgOTKVqHJMWIYxYfPg+wi5bDph++arhWZnsQkbZfhpHPUZJV+YGHDjHbS/5Hnv5wnGYgU6CXcUnYRLn2Y00/Ot7cHEgBuPMu/jDLK+oXvbCTsg9o/tiLW+7YTdd8BPC5fLMWkR1sn5adz/GXfbnGuYTMKdQvILAH6JwQZcuCMq7Lzeux58dvD6DptO2ISkrxY6NI5NmGZ+Gkd9RskXJhns1wB0Avh3AJNqN2J/iOQ3AIxLMbZcEpZdfOS+b8da33bChrpF0yIt4moMqM42YXI7ysUiskhE/kBExgH4EwB3A9gLYGXaAeaNsMb90R/eNnh9h7ORuFroDH184moc9RklX8R+KrKIdKH/YZq+j4B3nbiNu1x192nAzQx1+ypVDB2S1CMF3acZjfVqOHso3CRaqzRT03V1Vlpn6NOnGY31ajh7UAOOSVhJ4f7bHvL/jIOmU61K6IEepIVOEMUjrO0Uqb25ihpwTMpVzUiA5jNZNYd4NFNO0E7OHtSAYxJWgrhs0UW+y100nah9CtTC4dtzpkFY5x2ksZZ57EENOCbNnNXgYkbSbFav5hCPpmrADrY3V1EDjklzkyLuHRDN7pOLWqSJTsK5jRpwTMIOiFWXLwr4jHsHRFRWH6yFGnAcwnQO0lgfTW8PasAxiJr5f+5z1/oud9F0okoJQVq4WI5JkzCdA9ub1tmtQQ04Bn0RZ0BcueA83+Uu1j2jOpUiaZEWIhJqpqqx/VhlwCSnkbyb5LKGZaeTXE7yHpLXp/n9UaWE0T3d/p9zMOuLymQDtXCwHJMW5aqgGvKoxSCNdZRhD20z4JpB7ia5zrN8LsmNJDdHGaiIbBGRKzyLpwN4TES+iP6nNadGs0bqouk0PQmn5mCMTnS6Tzsz4HvR/3y5I5AcCuB2APPQb54LSc4kOZvko55/4wO2+zKABSR/BuDpFOOPbNi7zprV1OdsJOy2nECxtEiLqFquamw/sW/G0ywispLkFM/icwFsFpEtAEDyQQDzReSbAC4w3PRfALixtv1lAL6XVMxeomprD9zxcFOfs5Gog7xIWqRFs+2tr+LuDaBcI+sa8EQA2xted9WW+ULyZJJ3AvgAyb+rLf4JgC/Xlm8N+NxVJFeTXN3d7V83MyGqlDDnlhsCPuee6UTVGYO00PqkOVEGHKSxyWeVfJC1Aft10YEuJyLviMjVInJmLUuGiKwTkYtry/8m4HOLRaRTRDo7Ovwf42JClJHOfnxpU5+zkea1cK8enhZR9fIgjQE325yLZG3AXQAmN7yeBGBHRrFE0mz25uJ5mXoznvRpRSvt6OwgawN+EcDZJKeSHAFgAYDlGccUSLMHhIvDwWYPcDVgc1q5sbqWeuygnaehLQGwCsAMkl0krxCRMoBFAJ4AsAHAUhFZ366Y4hJlOouX+D+hyUXTiTrAg7TQpzWYE9VugjQ2+aySD9pmwCKyUEROFZHhIjJJRO6uLV8hItNrdd1vtCueZohq1BM2+fcdLh4MrWghIRcXKEeJ6vCDNAbcHHW5SNYlCKuIatTzb7zGd7mL9bgoAw7SQqT/Ci8lmqhRRpDGJp9V8oEacAyavYqrVHYv62ullODiiCANWpuEU41tQA04Bq1ksi5lwSLSmjk4eFZIGrRSRlADtgM14BhENeqnvvL1pj9rE+WqICqhD9NC65NmtNLedLLTDtSAYxBlHGs/dWngey4ZsMm+hGmhBmxGVB03vL3pKMMG1IBjEGU8150/I/A9l0zHZIInTAu9I5oZrbQ31dgO1IBj0EqjdikjabUzcWk0kCZhT+CO/qxqbANqwDFoaRLOoYyk1c5EzcGMVk4l007ODtSAYxBlHFs+/MnA91w6IEw6k3At3BkNpElUmwnTWM8DtgM1YENMTr165OY7A99zKesz2ZdQLdQcjIjq6MI01k7ODtSADTE59Wr+DVcHvufSAWFioOFaqAFHUalK5BWDYRqXyvpoehtQAzbExDSmPR/8RCSXTKdVLVwaDaRF6+3NnQ7fZdSADWn16i2XJuFaLSG4pEVatNpJaSdnB2rAhugBcZRWsyvNzqJptZPSu87ZgRqwIeVq9AFxy5MbA99zyXRMhsfhWrjTGaWFSYcdprGIW23OVdSADTEpQcx+7EfBn3fIdEzMIUwLPQsimlbbG+BWm3MVNWBDTExnzre+GvieSweDiYGGaeFSOSYtWm1vgHZ0NqAGbEirBurSwdCqFi51RmmRhEaqc/5RAzakddNxpx6nBpw+SXTYOtLIP2rAhpiYxiM3faelz9uCyU1iwrVwpzNKi1bbW/82VOe8owZsSMlgUmTX9FmB7zllwAbZWZgWLj6iKWlMstcwjQG3yl6uogZsiImBXrXw44HvuTQcbFULwC090sDEPKM0dqnTdxU1YENarnuWxZmsL4kr2XR4HE4S+mgnl3/UgA1p1YCrIqg48Dh2k5vEmKCXI4eTRPaqJYj8owZsiMnE09p5l4S+70LWZ2oM0VqoOYRhYp6qsf2oARtikrE9dd3Noe+7MCQ03YciaJEmJvpEaawGnH/UgA0xacyf/dKFLW8j75iWDqK00OFxOEm0N9U4/6gBG2JyQEzYvL7lbeQd08w1Wgv7yzFpYtLRRWncykM9lfZglQGT/AzJu0g+QvL82rLfInknyWUkr0nru5NozK3eUzgPJLUPLnRGaZJEiUYnOvNP2wyY5D0kd5Nc51k+l+RGkptJXh+2DRH5sYhcCeALAC6tLdsgIlcDuARAZ0rhGxnG/rEdoe+XKvY/JsbUGKK1UHMIw2SEEKWxdnL5p50Z8L0A5jYuIDkUwO0A5gGYCWAhyZkkZ5N81PNvfMNH/6H2ufp2/gzAswB+mlbwJtnEXQ8+G/q+ydV0ecf0oI7SQrOzYETEqH4b2d7UgHNP2wxYRFYC6PEsPhfAZhHZIiIlAA8CmC8ia0XkAs+/3eznnwE8LiJrGra9XER+D8BlacReNTz39SP3fTv0fRcyEtOJnSgt1ByCMdUmUmPt5HJP1jXgiQC2N7zuqi0L4loAcwBcTPJqACD5CZK3kvwugBV+HyJ5FcnVJFd3d3fHDtL0gPjoD28Lfd8FAzbdhyJokRamE5TRGts/4nKdYRl/P32WBbYaEbkVwK2eZc8AeCbsS0RkMYDFANDZ2Rm7VSZlFi5kfUntgwvlmLRIqjxTv+kR6XeYKXkg6wy4C8DkhteTAOzIKJZAksokXMhIkhrWagYcTFLauHL5u8tkbcAvAjib5FSSIwAsALA845gGYXpA3H/bQ+HbcaAmZ9qJRGqhBhxIr2E7idIYcKPTd5l2noa2BMAqADNIdpG8QkTKABYBeALABgBLRST87PIM0KzvKEntgwtapEWS2uhEXL5p51kQC0XkVBEZLiKTROTu2vIVIjJdRM4UkW+0K544mB4Qly26KPR9F2rAiWmhxhCIadYapTHgRptzmaxLEFagNeCjmA6Po9DLZINJsnPSkUa+UQM2QIfdR1Et0ifJrFVHGvlGDdgA0wNi1eWLQt93wXRMJxIjtVBjCMS0nURpHGdbSjaoARtgahbPfe7a0PddyEZMyyhRWpSreopUEKbtJEpjQGvAeUcN2ABT07lywXmJbCfPmB7QUVoAmp0FYaqLmcb2tzmXUQM2wPSAGN0Tfpmz7YZjepMYIFoLQLOzIJJqb4Aboy6XUQM2ICmjqFg+7E46m9I6sD9JniFie6fvOmrABpg24l1nzUpsW3kkTkdkpoW9nVGamGatJhrrKCPfqAEbUDY0igfueDhyHZsPiDgZq4kWNndGaWKqi1F701FGrlEDNsDUNOfcckPkOjYPu+MYpokWSV3U4RqmOhu1N+3kco0asAGmjXj240sNtmXvsDtO9m6mhZqDH6ZZq2psP2rABiSZtdp8QCQ9nLVZizTRK+GKgxqwAUlmrVbXgJM+C8JiLdJCRBJ9erbecyPfqAEbYGqai5esjFzHZtOJE7uJFvpUjMFUqoKqmOli1N40A841asARiIix8UzYFH0r4ySzm3YTJ3s30sLizigtVONioQYcQbkqMExIMP/GayLXsbkEEaeeaKSFZmeDiNNBq8b2owYcQdIZhM0ZiWqRPkl30OWqoGrx1ZeuowYcQdIlA5tNJ+nYbR4NpEUamqjO+UUNOII4jfepr3w9ch2bDTjOcNZMC83MvMSZNDPRGFADzjNqwBHEMcy1n7o0ch2bZ/7jnNJkooXNnVFaJN3eAD0TIs+oAUcQ54C47vwZiW4vb8Q5kE200AmiwcTJVk00BnSkkWfUgCPQiaejqBbpk0anpB1dflEDjiDpkoHNppN0LVFrk4NJI1tVnfOLGnAEcQxzy4c/GbmOzZeGxsmkTLToKwvE9CTrgpB0e4u7TaW9qAFHEKfxPnLzndHbs3g4GCc7M9GiKnY/ISQN4nRyJhrH3abSXtSAI4gzfJt/w9WR69icjcSJ3USL/m2qATeSdHsD7G5zrqMGHEEcg5j2/NMG27PzYChXqrGyVRMtAM3OvMRpH8YaW9rmioBVBkzyMyTvIvkIyfNry04nuZzkPSSvT/o7ky4Z9FXsvDQ0rUxVzWEgaXRIOsrIL20z4JpB7ia5zrN8LsmNJDdHGaiI/FhErgTwBQD1s9CnA3hMRL4IYGbScaeRsfZV7TOdtIzS1hFBWqShh44y8ks7M+B7AcxtXEByKIDbAcxDv3kuJDmT5GySj3r+jW/46D/UPgcALwNYQPJnAMzGZDGIYzy3PLnRaD0bM5K4B7G5FmoOjcQ5S0Y1tp+2GbCIrATQ41l8LoDNIrJFREoAHgQwX0TWisgFnn+72c8/A3hcRNbUtvEXAG4UkT8E8Kmk445jlrMf+5HZNi3MSOIexMZaqDkMIE7bUI3tJ+sa8EQA2xted9WWBXEtgDkALiZZnwL+CYAvk7wTwFa/D5G8iuRqkqu7u7tjBRjrScDf+mri28wLcWM21UKfjDyQOCMu1dh+hmX8/fRZFphyisitAG71LFsH4OKwLxGRxQAWA0BnZ2es8X8qNbkCGLD5du0rx6RJGqMjGzv8opB1BtwFYHLD60kAdmQUiy86K91PWlmUmsNRqlVBOYUzZGwseRWFrA34RQBnk5xKcgSABQCWZxzTAGJd/XXTdwy3ad8BEbfTMNZCzeEIcUdGphrbOOIqCu08DW0JgFUAZpDsInmFiJQBLALwBIANAJaKSPSTBttIHLPcNX2W0Xo2nhYUt9Mw1kLN4QhxtTDV2MYRV1Fo51kQC0XkVBEZLiKTROTu2vIVIjJdRM4UkW+0Kx4TKtV49yq4auHHjdazMQM+VKrEWt9cCzWHOnFHA6Yal8pVvelRTsl6Ei7XFH3iqVoVvNG9H69s34Oudw+l8h02jgbSIs120Vuu4tjhQ1PbvtIcasAhFPXqr4OlMtZ2vYe1v3kP+w6XU/2uvGvRLg73VbD2N++ltv17/utNTD15FKZ2jMKUk0epGecENeAQyjEzkrXzLjFaL691z53vHcYr2/dg0659Ld8m0nYt2kVfpYpfbt+DF7b2oLcvnhamGgNAb18Vr+3ch9d27sMQEqeddCymdYzCtHGjMWbUiLhhKwnBotWGOjs7ZfXq1Ubr7tp7GA88vy3xGGaddgLOn3VK4ttthnKlitd395cZdr53uO3fP2LYEPzemSfjnEknYcgQv9PC3aRaFWzYuRer3ngn9VFGFGNGDsfUjtGYNm4UTjvpOAwt0O8QRrUqONhXwehjzPJUki+JSGec78j6NLRcE7c++dkvXWi0Xh5qwIdKFfxi89u4+9k38ZN1OxM3X1MtSuUqntnYjftf2IbtPQcTjSGPiAjefPsA7n9hG55cv6sl8zXVOIp3D/Zhza/fxbKXuvDdlW/gJ+vewptvHyj0zfJFBE9t2IUlz29Dz4FSat+jJYgQ4tYnJ2w2O4Muy7rn4b4K1mx7Fy9v25PqBJipFnXe3teLZS91YfqE4/H708fhhGOHpxRZdux87zB+/np3YhOacTU2obevig1v7cOGt/Zh5IihmD7heMw45XiceuKxIIuRGYsIfvbabqzfsRcAsOyl7bjwg5MwbvQxiX+XGnAIWd4Dt1SuYn9vGfsO92H0McNwcos/fm+5gle27cFL296NXWtsJ5t27cObb+/H704Ziw+dMQbDhto/SNtzsIT/2vwONu3al3UosThYquCV7XvwyvY9OOG44XjfKf1m3IwR1YfzIgKSIIAhJEj0/0P/30Nq7/Uvb7/hiwie2dSNV7uOToge6K1g2UtduPADEzH+hGMT/T414BDiZqr7x3YYrddorvsO9/+/93D5yN/7DpcHnXc7dtQInDV+NM4ePxodxx9j3Djrkzyrf/1u7HN5W8FUCz/6KoJfvPEO1u/Yi49P78CZHaOsy75EBL/Zcwgb3tqHX+3Yi2oKcy2taByXvYf68MKbPXjhzR50HH8M3nfK8Zh+yvFHRiqlcrWhHQ9s13sP9eFAbyWWBiTQcfwxmDxmJCaNOQ6nnXRc6mduiAhWvv42Xtm2Z9B7h0oVLFvThT//wESceuJxiX2nTsKFsGbbu/h/G+PdPa0dnHjc8H4znjAap5zgPzQsV6p49TfvYfXWHhzobZ/xpsGUcSPxB9PHY2zOZ+tFBN37e7Fx5z5s3Lkv88m1djB21AgcLFVwuC/dNkYCE044FpPGHIdJY0bitJOOxTHDkjNkEcF/bX4HL2713jF3ICOGDcH895+GSWNG+sQYfxJODTiE57e8g1+88Y7xtj9y37fx3OeubTa0pjj+2GE4c/xonNUxGhNPOg4CYP2O9/DCmz2ZGkDSWgwh8YHTT8LksSMhIhAA/U1XIIIjr6X+Wvq1mXjScamfXbHnYAmv7dyHTbv24Z396U3YeMmiveWFISQmnHAMJo0Zicljj8PkMSNb+p1XvfEOnttidqwPH0p8+pzTcMbJowYsVwM2II4BP7flHayKYcDXnT/D+CkFaTByxFAMHcJcZF5Za1Hn2OFDMXXcSEzrGI0zTh6ZWNZ0oLeMjbv6M90sTt8D8qNxHjjxuOH40BljMOu0E2LPG8RNtABg2BDiU79zKqZ1jD6yrBkD1hqwQxxsY43XFg73VY7M6g8dQkweexzO7BiNaR2jjc7vFBEc6qug50DpyL/ufb34zZ5DKFjukmveO9SHn722G8+/+Q4+ePoYzJ50olFnu3prT2zzBYByVfDoq29h3m+fgrMnHN9MyADUgJUCUakKtr59EFvfPoifbtiNCSccizM7RmFax2icPGoE9h7uQ8+BEt49WELPgT70HOhFz4G+1OubSnIc6K3g56+/jRe29uD9k07C+08/CSNH+Nvcmm3v4uevv930d1WqghVrd+JPRPC+U05oahtqwAly/20PZR1CbrBBi117D2PX3sP4xRvvYAiZypkKaWKDxlnR21fF82/2YM22dzFr4on40BljBpxb/svtexKZYK+K4Cfrdsa+bUEdNWBFAawzX8WMvorglW178Or29/C+U4/H704Zix17DuFnr+1O7DtEgP/81a6mPmv/We454rJFF2UdQm5QLdJHNTanKoJf7diL+1ZtxVMbmjPLNNAMWFGUwpC3gY5mwIqiKBmhBpwgqy5flHUIuUG1SB/V2H7UgBOkqFcl+aFapI9qbD9qwAly5YLzsg4hN6gW6aMa248acIKM7snfjXuyQrVIH9XYftSAFUVRMkINOEF2nTUr6xByg2qRPqqx/agBJ8gDdzycdQi5QbVIH9XYftSAE2TOLTdkHUJuUC3SRzW2HzXgBJn9+NKsQ8gNqkX6qMb2owasKIqSEWrAiqIoGaEGnCCLl6zMOoTcoFqkj2psP4V7JhzJfQBcfpDWOADN3+Y//+j+2Y3L+zdDRGI9n6iIt6PcGPfBeTZBcrXun73o/tkLSbOn/TagJQhFUZSMUANWFEXJiCIa8OKsA0gZ3T+70f2zl9j7VrhJOEVRlLxQxAxYURQlFxTKgEnOJbmR5GaS12cdT9KQ3EpyLclXmpmRzRsk7yG5m+S6hmVjSf4nyddr/4/JMsZmCdi3r5H8Te33e4Xkn2YZYyuQnEzyaZIbSK4n+ZXacld+v6D9i/UbFqYEQXIogE0A/hhAF4AXASwUkV9lGliCkNwKoFNEnDjPkuTHAewHcJ+I/HZt2f8C0CMi/1TrRMeIyN9mGWczBOzb1wDsF5F/zTK2JCB5KoBTRWQNyeMBvATgMwC+ADd+v6D9uwQxfsMiZcDnAtgsIltEpATgQQDzM45JCUFEVgLo8SyeD+D7tb+/j/5Gbx0B++YMIvKWiKyp/b0PwAYAE+HO7xe0f7EokgFPBLC94XUXmhAs5wiAJ0m+RPKqrINJiQki8hbQfxAAGJ9xPEmziOSrtRKFlcNzLySnAPgAgOfh4O/n2T8gxm9YJAOmzzLX6i8fE5EPApgH4K9qw1zFHr4D4EwA7wfwFoD/nWk0CUByNICHAPy1iOzNOp6k8dm/WL9hkQy4C8DkhteTAOzIKJZUEJEdtf93A/i/6C+7uMauWv2tXofbnXE8iSEiu0SkIiJVAHfB8t+P5HD0m9P9IlJ/fIczv5/f/sX9DYtkwC8COJvkVJIjACwAsDzjmBKD5KjaZABIjgJwPoB14Z+ykuUAPl/7+/MAHskwlkSpG1ONP4fFvx9JArgbwAYR+beGt5z4/YL2L+5vWJizIACgdkrIvwMYCuAeEflGthElB8lp6M96gf6bLD1g+/6RXALgE+i/g9YuADcC+DGApQBOB7ANwH8TEesmswL27RPoH7oKgK0A/rJeL7UNkucB+DmAtQCqtcV/j/46qQu/X9D+LUSM37BQBqwoipInilSCUBRFyRVqwIqiKBmhBqwoipIRasCKoigZoQasKIqSEWrAiqIoGaEGrCiKkhFqwIriA8k5JH+QdRyK26gBK4o/5wB4OesgFLdRA1YUf84B8DLJY0jeS/J/1q7/V5TEGJZ1AIqSU85B/526ngDwf0TkhxnHoziI3gtCUTzUbjP4NoBfo/9mKqsyDklxFC1BKMpgZqL/9qVlAJWMY1EcRg1YUQZzDoBfoP+e0d8jOSHjeBRHUQNWlMGcA2CdiGwC8LcAltbKEoqSKFoDVhRFyQjNgBVFUTJCDVhRFCUj1IAVRVEyQg1YURQlI9SAFUVRMkINWFEUJSPUgBVFUTJCDVhRFCUj/j/Ne0ubHJPv4wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "ax.fill_between(k,Ek, alpha=0.5)\n",
    "#ax.plot(k,Ek)\n",
    "\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlabel('$k$')\n",
    "ax.set_ylabel('$E(k)$')\n",
    "ax.set_xlim([0,25])\n",
    "\n",
    "ke = np.zeros((k.shape[0],2))\n",
    "ke[:,0] = k\n",
    "ke[:,1] = Ek\n",
    "ke = ke[ke[:, 1].argsort()]\n",
    "\n",
    "ax.plot([ke[-1,0],ke[-1,0]], [0,ke[-1,1]], 'r--', linewidth=1.0 )\n",
    "ax.plot([ke[-2,0],ke[-2,0]], [0,ke[-2,1]], 'r--', linewidth=1.0 )\n",
    "ax.plot([ke[-3,0],ke[-3,0]], [0,ke[-3,1]], 'r--', linewidth=1.0 )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ba2aaaa",
   "metadata": {},
   "source": [
    "It is observed that the only the wave numbers (w1, w2, and w3) that were used to create the artifical velocity field contain energy. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff37423f",
   "metadata": {},
   "source": [
    "# Calculation of the kinetic energy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33a5511a",
   "metadata": {},
   "source": [
    "The blue area under the curve, $\\int_{0}^{\\infty} E(k) dk$, is calculated as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "id": "da7463ca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "123.1812731964925\n"
     ]
    }
   ],
   "source": [
    "blue_area = np.trapz(Ek,k)\n",
    "print(blue_area)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "id": "5970a366",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "137.24999999999983\n"
     ]
    }
   ],
   "source": [
    "kinetic_energy = 1/2*np.average(u**2+v**2+w**2)\n",
    "print(kinetic_energy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9addcb0f",
   "metadata": {},
   "source": [
    "There is about 11% deviation between $\\langle u_{i}^2\\rangle$ and $\\int_{0}^{\\infty} E(k) dk$. The main reason is that the resolution of the wave number (1 $m^{-1}$) is poor. This can be enhanced by increasing the domain size, i.e., $L$ along with the sampling resolution, i.e., dim."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "id": "cf187422",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax.plot?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "id": "86c4ac90",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m~\\AppData\\Local\\Temp/ipykernel_36364/656958259.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mk\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mEk\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'x' is not defined"
     ]
    }
   ],
   "source": [
    "x[:,0] = k\n",
    "x[:,1] = Ek"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "id": "92510425",
   "metadata": {},
   "outputs": [],
   "source": [
    "del x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "id": "957fe23b",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m~\\AppData\\Local\\Temp/ipykernel_36364/32546335.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mx\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'x' is not defined"
     ]
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "id": "29bd863f",
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "'type' object cannot be interpreted as an integer",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m~\\AppData\\Local\\Temp/ipykernel_36364/2384889276.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0margsort\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mreversed\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m: 'type' object cannot be interpreted as an integer"
     ]
    }
   ],
   "source": [
    "x = x[x[:, 1].argsort()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "id": "7e13d02b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2.80000000e+01, 8.53095101e-32],\n",
       "       [6.40000000e+01, 1.07878904e-31],\n",
       "       [4.40000000e+01, 1.16781817e-31],\n",
       "       [2.60000000e+01, 1.31258257e-31],\n",
       "       [6.10000000e+01, 1.46476643e-31],\n",
       "       [6.30000000e+01, 1.55567721e-31],\n",
       "       [3.70000000e+01, 1.58118411e-31],\n",
       "       [2.30000000e+01, 1.74837882e-31],\n",
       "       [2.70000000e+01, 1.87136567e-31],\n",
       "       [6.20000000e+01, 1.95237896e-31],\n",
       "       [5.70000000e+01, 2.04433781e-31],\n",
       "       [5.00000000e+01, 2.08700632e-31],\n",
       "       [2.50000000e+01, 2.14084901e-31],\n",
       "       [4.30000000e+01, 2.26758392e-31],\n",
       "       [5.30000000e+01, 2.60947743e-31],\n",
       "       [5.60000000e+01, 2.77440345e-31],\n",
       "       [4.60000000e+01, 2.83754465e-31],\n",
       "       [3.50000000e+01, 2.96072625e-31],\n",
       "       [2.90000000e+01, 2.96394909e-31],\n",
       "       [5.20000000e+01, 3.20842591e-31],\n",
       "       [1.20000000e+01, 3.26652575e-31],\n",
       "       [5.90000000e+01, 3.48840571e-31],\n",
       "       [5.50000000e+01, 3.51332666e-31],\n",
       "       [4.80000000e+01, 3.52226167e-31],\n",
       "       [5.10000000e+01, 5.31530806e-31],\n",
       "       [5.80000000e+01, 5.41820991e-31],\n",
       "       [1.30000000e+01, 5.61221947e-31],\n",
       "       [6.00000000e+01, 5.84655148e-31],\n",
       "       [5.40000000e+01, 7.26882048e-31],\n",
       "       [5.00000000e+00, 7.65081770e-31],\n",
       "       [3.40000000e+01, 7.86358086e-31],\n",
       "       [3.00000000e+01, 8.07802620e-31],\n",
       "       [2.20000000e+01, 8.82700158e-31],\n",
       "       [1.00000000e+01, 9.09626997e-31],\n",
       "       [1.00000000e+00, 9.13028866e-31],\n",
       "       [1.10000000e+01, 9.36372338e-31],\n",
       "       [3.60000000e+01, 1.26893488e-30],\n",
       "       [4.50000000e+01, 1.29355667e-30],\n",
       "       [1.90000000e+01, 1.44661723e-30],\n",
       "       [3.00000000e+00, 1.55138155e-30],\n",
       "       [4.20000000e+01, 1.58558285e-30],\n",
       "       [1.40000000e+01, 1.86602071e-30],\n",
       "       [3.80000000e+01, 2.06881305e-30],\n",
       "       [2.10000000e+01, 2.09157929e-30],\n",
       "       [2.00000000e+01, 2.17077905e-30],\n",
       "       [4.00000000e+00, 2.31547034e-30],\n",
       "       [4.70000000e+01, 2.56807915e-30],\n",
       "       [4.90000000e+01, 2.74299316e-30],\n",
       "       [3.10000000e+01, 2.80099895e-30],\n",
       "       [3.30000000e+01, 3.28863864e-30],\n",
       "       [4.10000000e+01, 3.60850092e-30],\n",
       "       [6.00000000e+00, 4.41105244e-30],\n",
       "       [3.20000000e+01, 4.61911529e-30],\n",
       "       [1.80000000e+01, 4.83745800e-30],\n",
       "       [3.90000000e+01, 5.42864508e-30],\n",
       "       [2.40000000e+01, 5.80176943e-30],\n",
       "       [9.00000000e+00, 6.80023014e-30],\n",
       "       [1.50000000e+01, 1.19370384e-29],\n",
       "       [1.70000000e+01, 1.34775565e-29],\n",
       "       [4.00000000e+01, 1.75411267e-29],\n",
       "       [7.00000000e+00, 1.87720438e-29],\n",
       "       [1.60000000e+01, 1.95158973e+01],\n",
       "       [8.00000000e+00, 3.79959552e+01],\n",
       "       [2.00000000e+00, 6.56694206e+01]])"
      ]
     },
     "execution_count": 192,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adb4b10e",
   "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
}