simple_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple example of diag_python\n",
    "**In the below, lines from `import sys` to `geom_set` are initialization for data analysis.  \n",
    "Users import a function of which he would like to analyze (`phiinxy` in this example).**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAGFCAYAAAASI+9IAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAilklEQVR4nO3dfdzt9Zzv8dfbTimiUtId1RFnCslsYXCm5CZ3kzBOxpjQPDqIycGhNOcxDM2EGWYYPGYPkZlmEiM14zaNhJHs7miXmz2IOiUhSpS9+5w/1u9a1t6u69rXda3733o9H4/f41q/7+9mfdaS9dnfm9/3m6pCkiSAO407AEnS5DApSJK6TAqSpC6TgiSpy6QgSeoyKUiSukwKkqQuk4L6luS7SR43hPs+JckXktyU5Pok70myfc/x9ye5PcktPduqRe73v5v7/CzJqUm26Tm2d5LPJrk1ydc3/zzjulYaNZOCJtk9gDcCuwO/BewBvGWzc95cVXfr2TbOd6MkTwROAA4D7gvsC7y+55R/AS4F7gmcBHw4yS7jvFYai6pyc1vxBvwjcAfwC+AW4NVDfK9nAF/r2X8/8MYlXvvPwF/07B8GXN+8vj9wG7B9z/HPAy8a57VubuPYrCmoL1X1POB7wNOq8y/1N29+TpL7NE1AC21/sMS3+x/Aus3KXpLkx0kuTvLMRa49ALi8Z/9yYNck92yOfbuqbt7s+AFjvlYaua3GHYDar6q+B+zQzz2SPB44Gnh4T/HbgVcCPwWeAHwwyfVV9cV5bnG35rw5c6+3n+fY3PE9xnytNHLWFDTxkjyCTjPMs6rqm3PlVXVJVf2oqjZU1ceB0+k0Mc3nFuDuPftzr2+e59jc8bl/wY/rWmnkTAoahEWn2m2aj25ZZHvuItceBJwDvLCqzltCHFng2DrgwJ79A4EfVNWPmmP79o5sao6vG/O10uiNu1PDbfo34ELg2CHc94HAD4D/ucDxZ9FpgrkTneajm4FDFjj3cOB6YH86TVn/AZyy2Wf4K+AuwJHATcAu47zWzW0c29gDcJv+DTiCTmfzTcCrBnjf99EZ2XRLz7au5/jn6bTB/4xOB+1RPcfu05x/n56yVzRJ5mfNvbfpObY3cD6dUVTfAB63WSxjudbNbdRbqlxkR5LUYZ+CJKnLpCBJ6jIpSJK6TAqSpK6JSgpJdkjy4Wa2yKuSPDLJTknOTfKt5u+O445TktpqokYfJTkN+HxVvSfJ1sB2wGuBH1fVKUlOAHasqtcsdp/tdt6udtj7HiOIWNK0u+7i62+sqr5mpr1fUrf2EwN8qqoO7yeGQZmYpJDkHsBlwL7VE1SSb9B5IOm6JLsB51fVAxa71+6rd6s/XvvCocYrqR3ekL+4uKpW93OPPZJ6SR/X/yn0HcOgTFLz0T7AD4H3Jbm0WVDlrsCuVXVdc871wK5ji1CSWm6SksJWwEOBd1fVQcDP6SxO0tXUIOat2iQ5NsnaJGtv/WE/FTlJWp4Ad+5jmySTlBSuAa6pqi83+x+mkyR+0DQb0fy9Yb6Lq2pNVa2uqtXb7bLdSAKWJOgkha362CbJxCSFqroe+H6Suf6Cw4Ar6cyQeXRTdjRw9hjCk6QFtammMGlJ6mXA6c3Io28DL6CTuM5McgxwNfDsMcYnSa02UUmhqi4D5uuBP2zEoUjSks01H7VBWz6HJI3NXPNRG5gUJKlP1hQkSV1tqilMzOgjSdL4WVOQpD7ZfCRJ6mpT85FJQZL6ZFKQJG2iLT+mdjRLkrraktwkaWxsPpIkdTn6SJLU1aaagn0KkqQuawqS1CebjyRJXW1qPjIpSFKfrClIkrraVFOwo1mS1GVNQZL6ZPORJKmrTc1HJgVJ6pNJQZK0ibb8mNrRLEnqaktyk6SxCXDnfn5NNwwqkv6ZFCSpTwlsZVKQJEEnKdx51bijGAz7FCRJXSYFSerTXPPRSrct3z+nJrkhyRU9ZTslOTfJt5q/Ow7is5gUJKlPcx3NK92W4P3A4ZuVnQCcV1X7Aec1+30zKUhSvwKs6mPbgqq6APjxZsVHAKc1r08Dnt7XZ2jY0SxJ/ep/8qOdk6zt2V9TVWu2cM2uVXVd8/p6YNe+ImiYFCRp/G6sqtUrvbiqKkkNIhCTgiT1azzTpP4gyW5VdV2S3YAbBnFT+xQkaRC26mNbmXOAo5vXRwNnr/hOPawpSFK/5jqah3X75F+AQ+j0PVwD/BlwCnBmkmOAq4FnD+K9TAqS1K8hNx9V1XMWOHTYoN/L5iNJUpc1BUnqV4vW42zJx5CkMWvJhHgmBUnqV4tqCvYpSJK6WpLbNlWEjW2py0mafC2qKbTkY0jSmLXk36ET13yUZFWSS5P8e7O/T5IvJ1mf5INJth53jJK0ibmawmifaB6KiUsKwPHAVT37bwLeVlX3A34CHDOWqCRpISaF4UiyJ/AU4D3NfoDHAh9uThnYnOGSpN80YTmKvwFeDWzf7N8TuKmqNjT71wB7jCEuSVpYizqaJ6amkOSpwA1VdfEKrz82ydoka2/94a0Djk6StmCIK6+N0iTltkcBv5fkycBdgLsDfwvskGSrprawJ3DtfBc3qxStAdht9e4DWWxCkpakRTWFifkYVXUicCJAkkOAV1XVc5N8CHgWcAZLnDM8FKvYOLxgJalXi5LCxDQfLeI1wCuSrKfTx/DeMccjSa01kbmtqs4Hzm9efxs4eJzxSNKihrzIzihNZFKQpKnSouajlnwMSRqzlvyaTkOfgiRpRFqS2yRpjOxTkCR12acw+XxOQdLImBQkSZtoya+pHc2SpK6W5DZJGiM7miVJXfYpSJK6TAqSpE20pPnIjmZJUpc1BUnql81HkqQuk4IkqatFScE+BUlSV0tymySNWUtGH5kUJKlfLWo+asnHkKQxMilIkjbRkuYjO5olSV2trCncwZ24mbuNOwxJs8LmI0lSl0lBktTVovUU7FOQJHVZU5Ckftl8JEnaREt+TVvyMSRpjKwpSJK6WtTR3Mqk8CvuzHXsPu4wJGnqtDIpSNJI2XwkSdpES35NW/IxJGmM7FOQJHW1qPnIJ5olSV0tyW2SNEYtqim05GNI0pgNuU8hyXeBm4GNwIaqWj2M92llUgjF1tw+7jAkzYrR1RQOraobh/kG9ilIkrpaWVOQpJHqv6awc5K1PftrqmrNZucU8OkkBfz9PMcHwqQgSf3qPyncuIQ+gkdX1bVJ7gWcm+TrVXVBX+86D5uPJGkQVvWxLUFVXdv8vQE4Czh4kOHPmZikkGSvJJ9NcmWSdUmOb8p3SnJukm81f3ccd6yStIm5msJKty3dPrlrku3nXgNPAK4Y8KcAJigpABuAV1bV/sAjgOOS7A+cAJxXVfsB5zX7kjRLdgW+kORy4CLgY1X1yWG80cT0KVTVdcB1zeubk1wF7AEcARzSnHYacD7wmjGEKEnzG/KQ1Kr6NnDg8N7h1yYmKfRKsjdwEPBlYNcmYQBcTydjLmoHbuJpnDO8ACW1xvsHdSMnxBuOJHcD/hV4eVX9LEn3WFVVMxxrvuuOBY4F2Pk+dxlFqJLU0aJpLiapT4Ekd6aTEE6vqo80xT9IsltzfDfghvmurao1VbW6qlbfY5etRxOwJMHQO5pHaWKSQjpVgvcCV1XVW3sOnQMc3bw+Gjh71LFJ0qyYpBz1KOB5wNeSXNaUvRY4BTgzyTHA1cCzxxOeJC2gRc1HE/MxquoLdL7a+Rw2ylgkabnKjmZJEkAFNrbk13TJfQrN08QjGScrSRqP5eS21wB/0yz08NqeZwcmzg4bfsaRP/7EuMOQNCtaVFNY8seoqkuAQ5M8E/hkko8Ab66qXwwtOkmaAhXYsKqfwZx3DCyWfi3rUzTDRr8BvBt4GfCtJM8bRmCSNC0qYeNWW614myRLjibJF4F9gHXAhcDzga8Dxyd5TFUdO5QIJWkKbFzVjuFHy0lRxwJXVtXm00y8rJm8TpI05ZbTp7BukcNPGUAskjSVirCxJTPiDaQxq5nWVZJmUhE2mBQkSXM2tuTndDkdzQGeC+xbVX+e5D7AvavqoqFFt0I/32pbvrzTA8YdhqSpcFnfd2hT89FyhqS+C3gk8Jxm/2bgnQOPSJI0Nsup7zy8qh6a5FKAqvpJEhcukDTz2lRTWE5S+FWSVUABJNmFSXoMT5LGaBaTwtuBs4B7JTkZeBbwf4cSlSRNkZkcfVRVpye5mM7aBgGeXlU+tCZJLbKc0UdvqqrX0JnaYvMySZpZnT6FdgxJXc7oo8fPU/akQQUiSdNsI6tWvE2SLaa2JC8GXgLsm+SrPYe2B/5zWIFJ0rSYtdFH/wx8AvhL4ISe8pur6sdDiapPX//l/jzym2vHHYakqbDQ0vBLVzA7Hc1V9VPgp8BzkuwI7AfcBSAJVXXBcEOUJI3Kcjqa/xg4HtiTznPhjwC+BDx2KJFJ0tSYzY7m44GHAVdX1aHAQcBNwwhKkqbJXJ/CTHQ09/hlVf0yCUm2qaqvJ3HWOUliNp9ovibJDsBHgXOT/AS4ehhBSdI0mbXRRwBU1ZHNy9cl+SxwDzqjkiRJLbGcjuZtgGcCe/dc9xDgzwcelSRNkZmc+wg4m87Q1IuB24YTzoB8Ezhk3EFImiVtGX20nE+xZ1UdPrRIJGlKtalPYTlDUv8zyYOGFokkaeyWU1N4NPD8JN+h03wUoKrqwUOJTJKmRJtqCstJCs6IKkkLmLmO5qrymQRJmkeb1lNYytTZX6iqRye5mc5kgOn9W1V3H3KMkjTRZqr5qKoe3fzdfvjhSJLGaSk1hVcsdryq3jq4cAbkvsDfjjsISVPhKYO5zczUFOissAbwADqzpJ7T7D8NuGgYQUnSNJmpJ5qr6vUASS4AHlpVNzf7rwM+NtToJGkKzFRHc49dgdt79m9vyiRp5s1S89GcDwAXJTmr2X868P5BByRJGp/lPKdwcpJPAI9pil5QVZcOJyxJmh4zNSQVIEnoTIh3CXDJcEOSpOkyc0mhqirJxwEnxJOkeczM6KMelyR5WFV9ZWjRLCLJ4XSePlgFvKeqTlno3N+ui1l7W0YWm6Tp5S/FppaTFB4OPDfJ1cDPGeEsqUlWAe8EHg9cA3wlyTlVdeWw31uStmRWh6Q+cWhRbNnBwPqq+jZAkjOAIwCTgqSxm7k+BejMkprkQH49+ujzVXX5cML6DXsA3+/Zv4ZOzUWSJkJbksKSV15LcjxwOnCvZvunJC8bVmDLleTYJGuTrP3hz8YdjaRZMjfNxUq3pUhyeJJvJFmf5IRhfZblNB8dAzy8qn4OkORNwJeAdwwjsM1cC+zVs79nU9ZVVWuANQCr75caQUySNBKj7FddTlIIsLFnfyOj67j/CrBfkn3oJIOjgD8Y0XtL0qJG0NE8sn7V5XyK9wFfbqa5SBPQewcd0HyqakOSlwKfojMk9dSqWjeK95akpRhyn8LI+lWX09H81iTnA4+ms/LaSKe5qKqPAx9fyrn/77/g9c8YckCS1BjA6KOdk6zt2V/TNImP3HLrOxvpJIQC7hh8OJI0nfpMCjdW1epFjm+xX3VQVjL6aGcmcPSRJLVYt181ydZ0+lXP2cI1KzIto48kaWINe+W1UfarTsvoI0maWKOY5mI5/ar9WOnoI+gssnPqwCOSpCnUlieaVzr6CFxkR5JaZ8lJIclpwPHNQjsk2THJqVX1wqFFJ0lTYCYnxAMeXFU3ze1U1U+SHDT4kPq3+z3gz3533FFImgavG8AYnmF3NI/ScpLCnZLsWFU/AUiy0zKvl6TWmsX1FP4a+FKSDzX7vw+cPPiQJGm6zGTzUVV9oHkM+7FN0TNc+UyS2mVZ9Z0mCZgIJKnHTNYUJEkLMylIkoB2jT5a8oR4kqT222JNIcnNdKbK/o1DQFXV3QcelSRNkVHMfTQqW/wUVbX9KAIZpEv3PZDt/+Uz4w5D0jS46y4Duc1M9ikk2RHYD7jLXFlVXTDooCRpmszk6KMkfwwcT2fFn8uAR9BZT+Gxi1wmSa03qx3NxwMPA66uqkOBg4CbhhGUJGk8ltN89Muq+mUSkmxTVV9P8oChRSZJU2RmOpp7XJNkB+CjwLlJfgJcPYygJGmazGSfQlUd2bx8XZLPAvcAPjmUqCRpisxkUuhVVZ8bdCCSNM3a0tG8lIfXXrHY8ap66+DCGYw7Lr2VW+562bjDkKSps5SawtzDaw+gM/pobp2ipwEXDSMoSZoms/ZE8+sBklwAPLSqbm72Xwd8bKjRSdIUmNU+hV2B23v2b2/KJGnmzWJS+ABwUZKz6EyG93TgtGEEJUkaj+UMST05ySeARzdFR1fVZUOJSpKmyEw1HyX5QlU9umcK7fQcc+psSTOvmKEhqVX16Obv1E2hLUmjMUOjj6bT3WGrx407CEnTYEP/t5ip5qNeSQ4EHtPsfr6qLh98SJKkcVny1NlJjgdOB+7VbP+U5GXDCkySpslGVq14myTLqSkcAzy8qn4OkORNdBbZeccwApOkadGmRXaWkxQCbOzZ30jPSCRJmlUzNc1Fj/cBX97s4bVThxGUJE2bSWsGWqnlPLz21iTnA49qinx4TZJaZikPr809tNYt6jnmw2uSZt5MDUmdxofWdv7tGzhi7TvHHYakKfDeAfSMFmHjHTOSFHpt9pzCBVX11cGHJElTpmDDhnYkhX6eUzjd5xQkqV0m4jmFJG+hs5Lb7cB/AS+oqpuaYyc2770R+JOq+lS/7ydJg1QVNm5ox5DUJdcUGO5zCucCD6yqBwPfBE4ESLI/cBRwAHA48K4k7aijSWqNTlJYteJtkqz0OQXoPKfw3kEEUVWf7tm9EHhW8/oI4Iyqug34TpL1wMF0aiiSNBmKiftxX6mVPKcwt8jOC6rq0iHE9ELgg83rPegkiTnXNGWSNDGqwoZfzUhSSHIX4EXA/YCvAe+qqmVPNpvkM8C95zl0UlWd3ZxzEp2JbE9fwf2PBY4FuOt9dlzu5ZIkllZTOA34FfB54EnAbwEvX+4bVdWiCxwkeT7wVOCwqpp7WO5aYK+e0/Zsyua7/xpgDcDqfVLvef9LlxuipBk0kDZwwh0b29HRvJRPsX9VPQggyXuBiwYdRJLDgVcDv1tVt/YcOgf45yRvBXYH9hvG+0tSXwqYoT6FX829qKoNyVAmRv07YBvg3Ob+F1bVi6pqXZIzgSvpNCsdV1UbF7mPJI1eZaaSwoFJfta8DrBtsx9gIHMfVdX9Fjl2MnByv+8hSUNTwIZ2rCSwlLmP2pH+JElb1I6eEUkat2WPyZxMJgVJ6ldhUpAkNUwKE+5XwA3jDkKSpk87k4IkjVLRM3h/ui1nllRJ0nyKzrzRK936kOR1Sa5NclmzPbmf+1lTkKRBGG+fwtuq6q8GcSOTgiT1q0UdzTYfSdL47Zxkbc927DKvf2mSryY5NUlf00RbU5CkfvVfU7ixqlYvdHCxpQeAdwNvaKJ4A/DXdNalWRGTgiT1a8jNR1taemBOkn8A/r2f92plUvjFvbfhq6++77jDkDQNXvPN/u8xxj6FJLtV1XXN7pHAFf3cr5VJQZJGbnwdzW9O8hA6qem7wP/q52YmBUmaYlX1vEHez6QgSf1q0RPNJgVJ6tfcE80tYFKQpH758JokqY2sKUhSv1pUUzApSFK/TAqT7Rr25AT+ctxhSJoKzxzMbUwKkiSgVTUFO5olSV3WFCSpXy2qKZgUJKlfPtEsSeryiWZJ0iZa0nxkR7MkqauVNYVfsC2XcdC4w5A0K+xoliR1mRQkSV0tGn1kn4IkqcuagiT1yyGpkqRN2KcgSQLsaJYk9WhRR3Mrk8LW3M7efHfcYUiaAteNO4AJ08qkIEkjZUezJKnLPgVJ0iZMCpIkoFUdzT7RLEnqmqikkOSVSSrJzs1+krw9yfokX03y0HHHKEm/Ya6jeaXbBJmY5qMkewFPAL7XU/wkYL9mezjw7uavJE0OO5qH4m3Aq4Gze8qOAD5QVQVcmGSHJLtV1aJDi+/GLfwOXxxiqJLa4kuDuEmLksJENB8lOQK4tqou3+zQHsD3e/avacrmu8exSdYmWXvrD38+pEglqd1GVlNI8hng3vMcOgl4LZ2moxWrqjXAGoB7r96j+rmXJC1Li0YfjSwpVNXj5itP8iBgH+DyJAB7ApckORi4Ftir5/Q9mzJJmiwT1mG8UmPvU6iqrwH3mttP8l1gdVXdmOQc4KVJzqDTwfzTLfUnSNLItahPYexJYQs+DjwZWA/cCrxgvOFI0jxMCsNTVXv3vC7guPFFI0mzZeKSgiRNHTuaJ9vW3M4+rqcgaVScOluStImW9ClMxMNrkqTJYE1Bkvrl6CNJUpcdzZKkLjuaJUldLWo+sqNZktTVyppCKFa1pS4naTq0pKbQyqQgSSNlR7MkqcuOZklSlx3NkqQ2sqYgSf1qUU3BpCBJ/WpRR7PNR5I0CBv72PqQ5PeTrEtyR5LVmx07Mcn6JN9I8sSl3M+agiRNtyuAZwB/31uYZH/gKOAAYHfgM0nuX1WLpqFWJoWNrOImdhh3GJJmSY3pbauuAkiy+aEjgDOq6jbgO0nWAwcDX1rsfjYfSdL47Zxkbc927ADuuQfw/Z79a5qyRbWypiBJU+bGqlq90MEknwHuPc+hk6rq7EEGYlKQpAlXVY9bwWXXAnv17O/ZlC3K5iNJaqdzgKOSbJNkH2A/4KItXWRNQZL6Nr4HFZIcCbwD2AX4WJLLquqJVbUuyZnAlXQerTtuSyOPwKQgSQMwvkeaq+os4KwFjp0MnLyc+5kUJKlv7XmkuZVJYSOruJntxx2GpJnRnsmP7GiWJHW1sqYgSaNl85EkqcukIEnahH0KkqSWsaYgSX2z+UiS1NWeIamtTQobWTXuECTNDGsKkqSu9tQU7GiWJHVZU5Ckvtl8JEnqak/zkUlBkvpmTUGS1NWemoIdzZKkLmsKktS39jQfTUxNIcnLknw9ybokb+4pPzHJ+iTfSPLEccYoSQvb0Mc2OSaippDkUOAI4MCqui3JvZry/YGjgAOA3YHPJLn/UhaflqTRsaYwaC8GTqmq2wCq6oam/AjgjKq6raq+A6wHDh5TjJLUepOSFO4PPCbJl5N8LsnDmvI9gO/3nHdNU/YbkhybZG2Stbf+8NYhhytJveZqCivdJsfImo+SfAa49zyHTmri2Al4BPAw4Mwk+y7n/lW1BlgDsNvq3au/aCVpOdozJHVkSaGqHrfQsSQvBj5SVQVclOQOYGfgWmCvnlP3bMokaYLYpzBoHwUOBUhyf2Br4EbgHOCoJNsk2QfYD7hoXEFK0vzmagqOPhqUU4FTk1wB3A4c3dQa1iU5E7iSzjd33FJHHq3CAUqStFwTkRSq6nbgDxc4djJw8mgjkqTlaE/z0UQkBUmabnY0S5K6rClIkrraU1OYlNFHkqQJYE1Bkvpm85Ekqas9zUfpPA7QLkl+CFw97jgWsTOdh/Mm3bTECcY6DNMSJ/QX632rapd+3jzJJ5sYVurGqjq8nxgGpZVJYdIlWVtVq8cdx5ZMS5xgrMMwLXHCdMU66exoliR1mRQkSV0mhfFYM+4Almha4gRjHYZpiROmK9aJZp+CJKnLmoIkqcukMEJJXpbk60nWJXlzT/mJSdYn+UaSJ44zxl5JXpmkkuzc7CfJ25tYv5rkoRMQ41ua7/SrSc5KskPPsYn6XpMc3sSyPskJ446nV5K9knw2yZXNf5/HN+U7JTk3ybeavzuOO1aAJKuSXJrk35v9fZrlfNcn+WCSrccd47QyKYxIkkOBI4ADq+oA4K+a8v2Bo4ADgMOBdyVZNbZAG0n2Ap4AfK+n+El0FjraDzgWePcYQtvcucADq+rBwDeBE2Hyvtfmvd9J5zvcH3hOE+Ok2AC8sqr2p7Ms7nFNfCcA51XVfsB5zf4kOB64qmf/TcDbqup+wE+AY8YSVQuYFEbnxcApVXUbQFXd0JQfAZxRVbdV1XeA9cDBY4qx19uAV9N5VHPOEcAHquNCYIcku40lukZVfbqq5h4lvZDOkq0wed/rwcD6qvp2s37IGU2ME6GqrquqS5rXN9P5wd2DToynNaedBjx9LAH2SLIn8BTgPc1+gMcCH25OmYg4p5VJYXTuDzymqeJ+LsnDmvI9gO/3nHdNUzY2SY4Arq2qyzc7NHGxbuaFwCea15MW66TFs6AkewMHAV8Gdq2q65pD1wO7jiuuHn9D5x8sdzT79wRu6vnHwcR+t9PAuY8GKMlngHvPc+gkOt/1TnSq5g8Dzkyy7wjD28QWYn0tnaajibBYrFV1dnPOSXSaQE4fZWxtk+RuwL8CL6+qn3X+Ed5RVZVkrMMVkzwVuKGqLk5yyDhjaSuTwgBV1eMWOpbkxcBHmrWnL0pyB525Uq4F9uo5dc+mbKgWijXJg4B9gMubH4Q9gUuSHMyExTonyfOBpwKH1a/HWI8l1kVMWjy/Icmd6SSE06vqI03xD5LsVlXXNU2FNyx8h5F4FPB7SZ4M3AW4O/C3dJoyt2pqCxP33U4Tm49G56PAoQBJ7g9sTWcCr3OAo5Jsk2QfOp24F40ryKr6WlXdq6r2rqq96VTFH1pV1zex/lEzCukRwE97mhbGIsnhdJoSfq+qbu05NFHfK/AVYL9mlMzWdDrBzxljPJto2uXfC1xVVW/tOXQOcHTz+mjg7FHH1quqTqyqPZv/No8C/qOqngt8FnhWc9rY45xm1hRG51Tg1CRXALcDRzf/ql2X5EzgSjrNH8dV1cYxxrmYjwNPptNpeyvwgvGGA8DfAdsA5zY1mwur6kVVNVHfa1VtSPJS4FPAKuDUqlo3rnjm8SjgecDXklzWlL0WOIVOU+cxdGYefvZ4wtui1wBnJHkjcCmdBKcV8IlmSVKXzUeSpC6TgiSpy6QgSeoyKUiSukwKkqQuk4IkqcukIEnqMiloWZJsTHJZkiuS/Fvv+gVLvP51SV7VvP7PJZx/ywpD7dtyY13gHjskeclgI5OGx6Sg5fpFVT2kqh4I/Bg4bqU3qqrfGVxY/Wmm7ljw/w99xLoDYFLQ1DApqB9fopmiOMlHk1zcrNp1bO9JSU5K8s0kXwAe0FN+S8/rBa+fT5I/alZbuzzJPzZlr2hqMFckeXnPuQuV753OSmgfAK4A9tpSrM01VyX5hybWTyfZdpHPcArw35ra1Vua8/4wyUVN2d/Pt/hPOqugPb55/cYk79jSdyINRFW5uS15A25p/q4CPgQc3uzv1Pzdls4P7D2b/d8GvgZsR2dGy/XAq3rvtYXrb5knhgPorLK289y1Pe9zV+BuwDo6awLMW95ctzedOfkfsdRYm2s2AA9p9s8E/nChz9Ccf0VP7L8F/Btw52b/XcAfzfMZ/wdwPvBc4GPAqnH/b+82G5sT4mm5tm0mTNuDzupc5zblf5LkyOb1XnRmJf0R8BjgrGpmME2y0MygC10/n8cCH6qqGwGq6sdJnte8z8+b9/lI895ZoPzS5l5XV2cVOZYR63eq6rLm9cV0fvgX+gzXb3btYXSSz1eaCfy2ZZ7pqKvqgmbm0lcAh9TkTpKoljEpaLl+UVUPSbIdnRk/j0vyVeBxwCOr6tYk59OZ635JmsVSVnx9n36+gmtu63m9kU6iPISlfYYAp1XViYu9QbOuxW7Aj6qzPKY0EvYpaEWaf03/CfBK4B7AT5ofw/9OZ3W5ORcAT0+ybZLtgafNc7vFrp/PfwC/n+SeAEl2Aj7fvM92Se4KHNmULVQ+n6XEupCFPsPNwPY9550HPCvJveZiT3Lf3huls5jN6XTWR76lWTNCGglrClqxqrq0qSXsAGyV5CrgG8CFPedckuSDwOV0mkm+Ms+tPgm8aL7rF3jfdUlOBj6XZCNwaVU9P8n7+fVCOu+pqksBFiqf575LiXUh836GqvpRki+ms47GJ6rq/yT5U+DTzWinX9EZwXV1E+t2wEeAV1bVVUneALypub80dK6nIEnqsvlIktRlUpAkdZkUJEldJgVJUpdJQZLUZVKQJHWZFCRJXSYFSVLX/weexCGoUAwbTwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "### Append diag_python packages to path ###\n",
    "import sys\n",
    "sys.path.append(\"./src/\")\n",
    "\n",
    "### Read NetCDF data *.nc by xarray ### \n",
    "from diag_rb import rb_open, rb_get_tri_filelist\n",
    "xr_phi = rb_open('../phi/gkvp.phi.*.nc')\n",
    "\n",
    "### Set geometric constants from GKV namelist etc. ###\n",
    "from diag_geom import geom_set\n",
    "geom_set(headpath='../src/gkvp_header.f90', \n",
    "         nmlpath=\"../gkvp_namelist.001\", \n",
    "         mtrpath='../hst/gkvp.mtr.001')\n",
    "\n",
    "# Plot phi[y,x] at t[it], zz[iz]\n",
    "from out_mominxy import phiinxy\n",
    "it = 250\n",
    "iz = 48\n",
    "phiinxy(it, iz, xr_phi, flag=\"display\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tips: To analyze interactively\n",
    "**1. xr_phi is just a xarray data set. `print(xr_phi)` shows the coordinate information.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.Dataset>\n",
      "Dimensions:  (kx: 1, ky: 13, t: 255, zz: 96)\n",
      "Coordinates:\n",
      "  * kx       (kx) float64 0.0\n",
      "  * ky       (ky) float64 0.0 0.05 0.1 0.15 0.2 0.25 ... 0.4 0.45 0.5 0.55 0.6\n",
      "  * zz       (zz) float64 -9.425 -9.228 -9.032 -8.836 ... 8.836 9.032 9.228\n",
      "  * t        (t) float64 0.0 0.1 0.2 0.3 0.4 0.5 ... 25.0 25.1 25.2 25.3 25.4\n",
      "Data variables:\n",
      "    rephi    (t, zz, ky, kx) float64 dask.array<chunksize=(255, 96, 13, 1), meta=np.ndarray>\n",
      "    imphi    (t, zz, ky, kx) float64 dask.array<chunksize=(255, 96, 13, 1), meta=np.ndarray>\n"
     ]
    }
   ],
   "source": [
    "print(xr_phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2. Please see the help of functions, which describes its feature, required inputs and outputs.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function phiinxy in module out_mominxy:\n",
      "\n",
      "phiinxy(it, iz, xr_phi, flag=None, nxw=None, nyw=None, outdir='./data/')\n",
      "    Output 2D electrostatic potential phi[y,x] at t[it], zz[iz].\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "        it : int\n",
      "            index of t-axis\n",
      "        iz : int\n",
      "            index of zz-axis\n",
      "        xr_phi : xarray Dataset\n",
      "            xarray Dataset of phi.*.nc, read by diag_rb\n",
      "        flag : str\n",
      "            # flag==\"display\" - show figure on display\n",
      "            # flag==\"savefig\" - save figure as png\n",
      "            # flag==\"savetxt\" - save data as txt\n",
      "            # otherwise       - return data array\n",
      "        nxw : int, optional\n",
      "            (grid number in xx) = 2*nxw\n",
      "            # Default: nxw = nxw in gkvp_header.f90 \n",
      "        nyw : int, optional\n",
      "            (grid number in yy) = 2*nyw\n",
      "            # Default: nyw = nyw in gkvp_header.f90\n",
      "        outdir : str, optional\n",
      "            Output directory path\n",
      "            # Default: ./data/\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "        data[2*nyw,2*nxw,3]: Numpy array, dtype=np.float64\n",
      "            # xx = data[:,:,0]\n",
      "            # yy = data[:,:,1]\n",
      "            # phixy = data[:,:,2]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(phiinxy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3. `flag=\"display\"` is convenient to see the result. On the other hand, choose `flag=None` to get data as a Numpy array and to draw figure as you like.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPSElEQVR4nO3df+xddX3H8eeL/gJFByg/attFDFVXjRPXEMyWLRMYSMyKi5qabdZJ0rhgpnGJAZu4mI1EZuKSOTdtollNmIg/SKvDKRAWsx+AnSK21NqKEkr4MYMgoBbbvvfHPcwbbPn223Nu7/f7/TwfyU3P/Zx7znm/09tXz/fccz/fVBWSpLacMO0CJEnHn+EvSQ0y/CWpQYa/JDXI8JekBhn+ktSg3uGf5MQkdyT5dpKdST7YjZ+d5PYke5N8NsnS/uVKkoYwxJn/fuB1VfWbwKuBS5KcD1wD/F1VnQP8GLh8gGNJkgbQO/xr5Inu6ZLuUcDrgM9341uAy/oeS5I0jMVD7CTJIuB/gHOAjwHfBx6tqgPdS/YBK46w7UZgI0CWLv2tJWecMURJktSMp/bt+1FVnT6bbQYJ/6o6CLw6ySnADcDLZ7HtZmAzwIkrV9WqK947REmS1Iy9V7333tluM+jdPlX1KHAr8FrglCRP/+eyErh/yGNJko7dEHf7nN6d8ZPkJOAiYBej/wTe1L1sA7C177EkScMY4rLPcmBLd93/BOD6qvpykruB65L8DfAt4JMDHEuSNIDe4V9VdwHnHmb8HuC8vvuXJA3Pb/hKUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktSgQeb2GcwhWPTzaRchSQufZ/6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGjS3pncIlP8dSdLEGbWS1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBvUO/ySrktya5O4kO5O8uxs/LclNSfZ0f57av1xJ0hCGOPM/APxlVa0BzgeuSLIGuBK4papWA7d0zyVJc0Dv8K+qB6rqm93y48AuYAWwDtjSvWwLcFnfY0mShjHo3D5JXgycC9wOnFlVD3SrHgTOnGn7WlzsP+vAkCVJkg5jsA98k5wMfAF4T1X9ZHxdVRVQR9huY5LtSbYffOLJocqRJD2LQcI/yRJGwX9tVX2xG34oyfJu/XLg4cNtW1Wbq2ptVa1ddPJzhyhHkjSDIe72CfBJYFdVfWRs1TZgQ7e8Adja91iSpGEMcc3/t4E/Bb6T5M5u7P3Ah4Drk1wO3Au8ZYBjSZIG0Dv8q+o/gBxh9QV99y9JGp7f8JWkBhn+ktQgw1+SGmT4S1KDDH9JatCg0zv0dcKSQzzvrCemXYYkLXie+UtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lq0Jya3uHQwRN4/NHnTLsMSVrwPPOXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaNMjcPkk+BbwBeLiqXtmNnQZ8Fngx8EPgLVX142ct5vFw+q1LhihJkppx7zFsM9SZ/z8Dlzxj7ErglqpaDdzSPZckzQGDhH9VfR145BnD64At3fIW4LIhjiVJ6m+S1/zPrKoHuuUHgTMP96IkG5NsT7L9wM+fnGA5kqSnHZcPfKuqgDrCus1Vtbaq1i4+8bnHoxxJat4kw/+hJMsBuj8fnuCxJEmzMMnw3wZs6JY3AFsneCxJ0iwMEv5JPgP8N/CyJPuSXA58CLgoyR7gwu65JGkOGOQ+/6p66xFWXTDE/iVJw/IbvpLUIMNfkho0yGWfoSx+4he88D8fmnYZkrTgeeYvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1aE7N7XPwOUt47Nwzpl2GJM0ve2a/iWf+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkho0p6Z3WHTGU5x6xb3TLkOS5pfrZ7+JZ/6S1KCJh3+SS5LsTrI3yZWTPp4kaWYTDf8ki4CPAa8H1gBvTbJmkseUJM1s0mf+5wF7q+qeqnoKuA5YN+FjSpJmMOnwXwHcN/Z8Xzf2/5JsTLI9yfanHv3ZhMuRJMEc+MC3qjZX1dqqWrv0lJOmXY4kNWHS4X8/sGrs+cpuTJI0RZMO/28Aq5OcnWQpsB7YNuFjSpJmMNEveVXVgSTvAr4KLAI+VVU7J3lMSdLMJv4N36q6Ebhx0seRJB29OTW9w+plj/Gll/7btMuQpHll0TFsM/W7fSRJx5/hL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBc2punwcPLuOaR86ZdhmSNM/smfUWnvlLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JatCcmt7hkYefz+c+etG0y5CkeeYrs97CM39JapDhL0kNMvwlqUG9wj/Jm5PsTHIoydpnrLsqyd4ku5Nc3K9MSdKQ+n7guwP4I+AT44NJ1gDrgVcALwJuTvLSqjrY83iSpAH0OvOvql1Vtfswq9YB11XV/qr6AbAXOK/PsSRJw5nUNf8VwH1jz/d1Y78iycYk25NsP/CzJydUjiRp3IyXfZLcDJx1mFWbqmpr3wKqajOwGeA5Z6yqvvuTJM1sxvCvqguPYb/3A6vGnq/sxiRJc8CkLvtsA9YnWZbkbGA1cMeEjiVJmqVed/skeSPwUeB04F+T3FlVF1fVziTXA3cDB4ArjuZOn8X/+yQv+MR/9SlJknQUeoV/Vd0A3HCEdVcDV/fZvyRpMvyGryQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QG9f1NXoPK0iUsXr5y2mVI0vxy7+w38cxfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoPm1PQO+1+4lHsu//VplyFJ88sHZr+JZ/6S1CDDX5IaZPhLUoMMf0lqkOEvSQ3qFf5JPpzku0nuSnJDklPG1l2VZG+S3Uku7l2pJGkwfc/8bwJeWVWvAr4HXAWQZA2wHngFcAnwj0kW9TyWJGkgvcK/qr5WVQe6p7cBT/8OxnXAdVW1v6p+AOwFzutzLEnScIa85v8O4Cvd8grgvrF1+7oxSdIcMOM3fJPcDJx1mFWbqmpr95pNwAHg2tkWkGQjsBFg8a+dOtvNJUnHYMbwr6oLn219krcDbwAuqKrqhu8HVo29bGU3drj9bwY2A5y4YlUd7jWSpGH1mtsnySXA+4Dfq6qfjq3aBvxLko8ALwJWA3fMtL+TnrefV16wp09JktSc7x3D3D59J3b7B2AZcFMSgNuq6p1VtTPJ9cDdjC4HXVFVB3seS5I0kF7hX1XnPMu6q4Gr++xfkjQZfsNXkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1KC+3/Ad1KIc4uTF+6ddhiQteJ75S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWrQnJre4fH9J/L1PUf8tcCSpIF45i9JDTL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoMMf0lqUK/wT/LXSe5KcmeSryV5UTeeJH+fZG+3/jXDlCtJGkLfM/8PV9WrqurVwJeBD3TjrwdWd4+NwD/1PI4kaUC9wr+qfjL29LlAdcvrgE/XyG3AKUmW9zmWJGk4vef2SXI18DbgMeD3u+EVwH1jL9vXjT1wmO03MvrpAGD/D//k/Tv61jSHvRD40bSLmCD7m78Wcm+w8Pt72Ww3mDH8k9wMnHWYVZuqamtVbQI2JbkKeBfwV7MpoKo2A5u7Y22vqrWz2X4+sb/5bSH3t5B7gzb6m+02M4Z/VV14lPu6FriRUfjfD6waW7eyG5MkzQF97/ZZPfZ0HfDdbnkb8Lburp/zgceq6lcu+UiSpqPvNf8PJXkZcAi4F3hnN34jcCmwF/gp8GdHub/NPeuZ6+xvflvI/S3k3sD+fkWqauZXSZIWFL/hK0kNMvwlqUFzIvwX+jQRST6c5LtdDzckOWVs3VVdf7uTXDzFMo9Jkjcn2ZnkUJK1z1g3r3t7WpJLuh72Jrly2vX0leRTSR5OsmNs7LQkNyXZ0/156jRr7CPJqiS3Jrm7e2++uxuf9z0mOTHJHUm+3fX2wW787CS3d+/RzyZZOuPOqmrqD+D5Y8t/AXy8W74U+AoQ4Hzg9mnXeoz9/QGwuFu+BrimW14DfBtYBpwNfB9YNO16Z9nbbzD6gsm/A2vHxud9b10fi7raXwIs7XpaM+26evb0u8BrgB1jY38LXNktX/n0e3Q+PoDlwGu65ecB3+vej/O+xy4LT+6WlwC3d9l4PbC+G/848Ocz7WtOnPnXAp8moqq+VlUHuqe3MfreA4z6u66q9lfVDxjdHXXeNGo8VlW1q6p2H2bVvO+tcx6wt6ruqaqngOsY9TZvVdXXgUeeMbwO2NItbwEuO541DamqHqiqb3bLjwO7GM0wMO977LLwie7pku5RwOuAz3fjR9XbnAh/GE0TkeQ+4I/55QRxR5omYj57B6OfZmBh9ve0hdLbQuljJmfWL7+L8yBw5jSLGUqSFwPnMjpDXhA9JlmU5E7gYeAmRj+ZPjp2gnlU79HjFv5Jbk6y4zCPdQBVtamqVjH6pvC7jlddQ5mpv+41m4ADjHqcN46mNy0cNbp2MO/vAU9yMvAF4D3PuLowr3usqoM1mkl5JaOfTF9+LPvpPbHb0aoFPk3ETP0leTvwBuCC7o0H86S/WfzdjZsXvR2FhdLHTB5KsryqHugurT487YL6SLKEUfBfW1Vf7IYXVI9V9WiSW4HXMrokvrg7+z+q9+icuOyz0KeJSHIJ8D7gD6vqp2OrtgHrkyxLcjaj339wxzRqnICF0ts3gNXd3RRLgfWMeltotgEbuuUNwNYp1tJLkgCfBHZV1UfGVs37HpOc/vTdgklOAi5i9JnGrcCbupcdXW/T/vS6Own+ArADuAv4ErBi7JPtjzG6pvUdxu4mmU8PRh923gfc2T0+PrZuU9ffbuD10671GHp7I6NrjPuBh4CvLpTexvq4lNEdI99nNJvt1Gvq2c9nGE2v/ovu7+5y4AXALcAe4GbgtGnX2aO/32F0SeeusX9zly6EHoFXAd/qetsBfKAbfwmjk6u9wOeAZTPty+kdJKlBc+KyjyTp+DL8JalBhr8kNcjwl6QGGf6S1CDDX5IaZPhLUoP+D5dqUH72HWAaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "it = 250\n",
    "iz = 48\n",
    "data = phiinxy(it, iz, xr_phi, flag=None)\n",
    "xx = data[:,:,0]\n",
    "yy = data[:,:,1]\n",
    "phixy = data[:,:,2]\n",
    "plt.pcolormesh(xx,yy,phixy,shading=\"auto\")\n",
    "plt.xlim(-30,30)\n",
    "plt.ylim(-30,30)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**4. If you are familiar with Python, you can analyze everything by yourself. However, be careful that you also need to be familiar with the formulation of local \n",
    "flux-tube gyrokinetic simulation.**  \n",
    "For example, flux-surface-averaged electric field spectrum is not $\\int \\frac{k_\\perp^2}{2}|\\phi_{k_x,ky}(z)|^2 dz/(\\int dz)$, but $\\int \\sqrt{g} \\frac{k_\\perp^2}{2}|\\phi_{k_x,ky}(z)|^2 dz/(\\int \\sqrt{g} dz)$ including the Jacobian of the employed flux-tube coordinates.  \n",
    "The squared wavenumber is not $k_x^2 + k_y^2$, but $k_\\perp^2 = g^{xx} k_x^2 + 2g^{xy}k_xk_y + g^{yy}k_y^2$, including the metrics.\n",
    "\n",
    "See GKV manual for theoretical discription, and use `diag_geom` module to pick up geometrical quantities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.DataArray (zz: 96, ky: 13, kx: 1)>\n",
      "dask.array<add, shape=(96, 13, 1), dtype=complex128, chunksize=(96, 13, 1), chunktype=numpy.ndarray>\n",
      "Coordinates:\n",
      "  * kx       (kx) float64 0.0\n",
      "  * ky       (ky) float64 0.0 0.05 0.1 0.15 0.2 0.25 ... 0.4 0.45 0.5 0.55 0.6\n",
      "  * zz       (zz) float64 -9.425 -9.228 -9.032 -8.836 ... 8.836 9.032 9.228\n",
      "    t        float64 25.4\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEXCAYAAACDChKsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABLJElEQVR4nO3dd3hUZfbA8e9Jp4YSihACSG+hdxEFC1ZUEMWGroplLavb3P3ZVrevXbC7YpciKrYVC6IYFBNKQu+QhN4CgYS08/vjTmIICczcmWQmyfk8zzzJ3Llz77kk4cx73/c9r6gqxhhjjK/Cgh2AMcaY6skSiDHGGFcsgRhjjHHFEogxxhhXLIEYY4xxxRKIMcYYVyyBGGOMccUSiDHGGFf8SiAiUk9EwgMVjDHGmOrDpwQiImEicpWIfCoiu4DVwHYRWSki/xGRjpUTpjHGmFAjvpQyEZH5wFfAR8ByVS3ybG8CnAlcBXygqm9VQqzGGGNCiK8JJFJV80VkPPC+lvPm4n0CGaQxxpjQ49MtrFKJ4U3gndL9HyJyQ5l9jDHG1GBuO9FXA/OB90Uk0rPtzsCEZIwxpjpwm0BUVV8AZgNzRKQOIIELyxhjTKiLcPm+/QCq+oaI5ACfAnUDFpUxxpiQ52sn+vnAUlXdVmb7hcA0VY0LcHzGGGNClK8J5GWgL9ACpx9kGbDU83WlqhZWQozGGGNCkE8JpORNIn8GWgMrgP7ARGCfqsYHNrzAiYuL03bt2gU7DGOMqVZSUlL2qGqz8l5z2wdyhar2Ln4iIs8Bv3d5rCrRrl07kpOTgx2GMcZUKyKypaLX3I7COigi/YufqGoK0NnlsYwxxlRDJ22BiEihqpYtmHgjMFtEfgZSgF6ATSA0xphaxJsWyHHzO1R1LdAP+BynQ30VcH5gQzPGGBPKvOkDKbeXXVXzgBmehzHGVCg/P5+MjAxyc3ODHYqpQExMDPHx8URGRp58Zw+3nejGGOO1jIwMGjRoQLt27RCxohWhRlXZu3cvGRkZtG/f3uv3uUogIrIOSKPUPBBV3ezmWMbUaqkz4OtHICsDYuNh9IOQOCHYUQVcbm6uJY8QJiI0bdqU3bt3+/Q+t6OwXgR2AHuB84DlIpImIo+UKq5ojDmR1Bnw8V2QlQ6o8/Xju5ztNZAlj9Dm5ufjNoFco6q3q+oUVb0VOA2YBxwEnnB5TGNql68fgfycY7fl5zjbjakG3CaQLBFJLH6iqkuBkar6GDA8EIEZU+NlZfi2vZZ4Yf4GkjbsOWZb0oY9vDB/g+tj3nPPPTz11FMlz88991xuuummkue//e1veeIJ++zrK7cJ5BbgNRF5VUTuFJEpwBHPa1GBCc2YGi62gso/FW2vJRLjY7njnSUlSSRpwx7ueGcJifGxro85fPhwkpKSACgqKmLPnj2sWLGi5PWkpCSGDRt2zHsKCgpcny8YghGvq050VV0tIoOAy4BEYD3wkIjUA94LYHzG1FyjH4QPbgEt+mVbZB1new32l49XsHLbwRPu07xBNNe9uogWDaPZefAoHZvX5+mv1vH0V+vK3b97q4Y8dFGPCo83bNgw7rnnHgBWrFhBz5492b59O/v376du3bqsWrWKfv36ccYZZ9CnTx8WLFjAxIkT6dOnD7/73e8oKChg4MCBPP/880RHR9OuXTsmTZrExx9/TH5+PjNnzqRr167s3r2bq666im3btjF06FC+/PJLUlJSiIs7tlD53Llzeeihhzh69CgdOnTgtddeo379+hUe9/Dhw9x5550sX76c/Px8Hn74YcaOHcu0adOYPXs22dnZFBYW8vnnn3P99dezfPlyunTpwrZt25g6dSqpqamkpqaWtMJefvllVq5cyZNPPunDT+54Xk0kFJFTRKRL6Y2qWqiqM1X1AVV9SlX3quphVf2rXxEZU1v0uhwi60FkqaV0ek+skaOwfBVbJ5IWDaPJPJBLi4bRxNbxb2xOq1atiIiIYOvWrSQlJTF06FAGDx7MwoULSU5OplevXkRFOTdP8vLySE5O5te//jXXX38906dPJy0tjYKCAp5//vmSY8bFxbF48WJuu+02HnvsMQD+8pe/MGrUKFasWMH48ePZunXrcbHs2bOHv/71r3z11VcsXryYAQMGHHP7rLzj/u1vf2PUqFEsWrSIefPm8fvf/57Dhw8DsHjxYmbNmsX8+fN57rnnaNy4MStXruTRRx8lJSUFgAkTJpQkJYDXXnuNX/3qV379m4IXLRBVDQMQkTgRudyzeRewwMq3G+OHrHTIOwTnPwYDboQXT4cN86AwH8Jr7mDGE7UUihXftrprVEfe+mkrd5/ViWEd/FtuaNiwYSQlJZGUlMS9995LZmYmSUlJxMbGMnz4L123V1xxBQBr1qyhffv2dO7slPmbNGkSU6dO5Te/+Q0Al112GQD9+/dn9uzZACxYsIAPPvgAgDFjxtC4cePj4vjxxx9ZuXJlyTnz8vIYOnRoyevlHXfu3LnMmTOnJKHk5uaWJKezzz6bJk2alJz/7rvvBqBnz54kJjpd1fXr12fUqFF88skndOvWjfz8fHr16uX+H9PD61tYqpqGM/cDEYkDLvQM2T0K/KCq+/yOxpjaJMNTHbp1fwgLgzP/DO9NhGXvQr/rghtbEBUnjylX9WVYhziGdGh6zHO3ivtB0tLS6NmzJ23atOHxxx+nYcOG3HDDDSX71atXz6vjRUdHAxAeHu5T/4OqcvbZZ/Puu+96fVxV5f3336dLl2NuBPHTTz95He9NN93E3//+d7p27XrM9frDVSe6qu5R1Y9UdRZOPayeIjJORK4VEavKa4w3MlMgPBpa9HSedzkPWvWD+f+BgrzgxhZEqRlZxySLYR3imHJVX1Izsvw67rBhw/jkk09o0qQJ4eHhNGnShAMHDrBw4cLjOtABunTpwubNm1m/fj0Ab775JiNHjjzhOYYPH86MGc48nrlz57J///7j9hkyZAg//PBDyXEPHz7M2rVrT3jcc889l2effZbi9ZuWLFly0vOvXLmStLS0ktcGDx5Meno677zzDhMnTjzh+bzldhRWaV1whu6eD3QDGgTgmMbUfJkpcEpviPAMXBSBM/8PsrbCkjeCG1sQ3Tqyw3EtjWEd4rh1ZAe/jturVy/27NnDkCFDjtkWGxt7XCc3OLWhXnvtNS6//HJ69epFWFgYt9566wnP8dBDDzF37lx69uzJzJkzadmyJQ0aHPtfYrNmzZg2bRoTJ04kMTGRoUOHsnr16hMe94EHHiA/P5/ExER69OjBAw88UO5+t99+O7t376Z79+7cf//99OjRg9jYX0avTZgwgeHDh5d7a80VVXX9AJ4B7gTa+3Ocqnj0799fjQkZBXmqj7ZQ/eyPx24vKlJ95RzVx7qq5uUEJ7ZKsHLlymCHUCVyc3M1Pz9fVVWTkpK0d+/eVXr+goICzclxfm/Wr1+v7dq106NHj5a8fsEFF+hXX31V4fvL+zkByVrB/6v+FlN8QFX9a1caUxvtWgkFORA/4NjtIjDq/+D1iyDlNRhyW3DiM65s3bqVCRMmUFRURFRUFC+//HKVnv/IkSOceeaZ5Ofno6o899xzREVFceDAAQYNGkTv3r0ZPXp0wM7nVwKx5GGMS5nO8Epa9z/+tfanQ7sR8P3jTmd6lHedpCb4OnXqVGH/RFVo0KBBuUt3N2rU6KT9LG6ctA9ERJqJyE0n288Y44OMFKjbFBq3K//1UffD4d2wqGo/wRrjC2860b8F+lZyHMbULpnJTuujogqoCUOgw2j44Wk4eqhqYzPGS97cwioCGorIH4HsMo9DZbep6pGKDmSMAXIPwu410OOyE+935v/BK6Pgxxdg5O+rJjZjfOBNAjkfuAsYgDNEt345j7p4lr4VkT3A71T1zcoI2Jhqb9sSQCG+nP6P0uL7Q5fzYeGzMOhmqNOoKqIzxmsnvYWlqumq+ntVvVxVx6jqaaraR1U7qmpLVa0PhAONgHjgNuCpSo3amOos09PJ2arfyfc988+QmwULp1ZuTKEmdQY82RMebuR8DcAiW+Hh4fTp04eePXty+eWXc+SI3SzxVyAmEuIZLnxIVbcD/wKmB+K4xtRIGSnQpAPUbXLyfVv2gu5j4cfn4fDeyo8tFFTSSo116tRh6dKlLF++nKioKF544YVjXq/McuiFhZVfNlBVKSoqOvmOARSQBFLGYFW9vRKOa0z1p+q0QMrO/ziRM/4EedmQ9HTlxVWVPr8PXrug4sdHd5S/UuNHd1T8ns/v8ymEESNGsH79er799ltGjBjBxRdfTPfu3cnNzeWGG26gV69e9O3bl3nz5gEwbdo0xo4dyxlnnEGnTp34y1/+UnKst956i0GDBtGnTx9uueWWkmRRv359fvvb39K7d28WLlx4zPk3bNjAmDFj6N+/PyNGjCiZiX799ddz1113MWzYME499VRmzZpV8p7//Oc/DBw4kMTERB566CEANm/eTJcuXbjuuuvo2bMn6enpPProo3Tp0oXTTjuNiRMn8thjj7Fhwwb69fulxbtu3bpjnrvlUwIR7xbNPb74izHGcTATsneWP/+jIs27Qa/xzpDe7F2VF1uoKDzq23YfFRQU8Pnnn5dUo128eDFPP/00a9euZerUqYgIaWlpvPvuu0yaNInc3FwAFi1axPvvv09qaiozZ84kOTmZVatWMX36dH744QeWLl1KeHg4b7/9NuDUuBo8eDDLli3jtNNOOyaGyZMn8+yzz5KSksJjjz3G7bf/8pl7+/btLFiwgE8++YT77nMS49y5c1m3bh2LFi1i6dKlpKSk8N133wFOMrj99ttZsWIFu3bt4v3332fZsmV8/vnnJXNCOnToQGxsLEuXLgWccu6BKKjo60TCeSLyPvCRqpYUuheRKJx10SfhrI0+ze/IjKmJSirw+tACARh5Hyx/HxY8CWP+Efi4qtJ5/zzx60/29Ny+KiO2DdzwqevT5uTk0KdPH8Bpgdx4440kJSUxaNAg2rdvDzjl0O+8804AunbtStu2bUsm4J199tk0bdoUcEquL1iwgIiICFJSUhg4cGDJOZo3bw44fS7jxo07Lo7s7GySkpK4/PLLS7YdPfpLcrzkkksICwuje/fu7Ny5E3ASyNy5c+nbt2/JMdatW0dCQgJt27Ytqe/1ww8/MHbsWGJiYoiJieGiiy4qOe5NN93Ea6+9xhNPPMH06dNZtGiR63/LYr4mkDHAr4B3RaQ9cACIwelEnws8parBm4ZpTKjLTIbwKGjZ07f3xXWE3lfBz6/CsDuhYavKiS8UjH7Q6fMofRsrACs1FveBlOVtOfSyN2BEBFVl0qRJ/OMfxyf1mJgYwsPDj9teVFREo0aNyo0FfinnDpRU31VV/vSnP3HLLbccs+/mzZu9jn/cuHElC17179+/JBn6w9c+kL7A86o6HGgLjAb6qWpbVb3ZkocxJ5G5GFomQkT0yfcta+TvQQudEic1WeIEuOgZp8WBOF8veqZKVmocMWJEyS2otWvXsnXr1pI1OL788kv27dtHTk4OH374IcOHD2f06NHMmjWLXbucW4v79u1jy5YtJzxHw4YNad++PTNnzgSc5LBs2bITvufcc8/lv//9L9nZ2QBkZmaWnLO04cOH8/HHH5Obm0t2djaffPJJyWsxMTGce+653HbbbQFbD8TXFsh1wFQRWQv8D/ifqh4ISCTG1HSFBc4ckL7Xunt/43bOe1Neh+F3Q6OEgIYXUhInBGVp39tvv53bbruNXr16ERERwbRp00paBIMGDWLcuHFkZGRwzTXXMGCAcxvyr3/9K+eccw5FRUVERkYydepU2rZte8LzvP3229x222389a9/JT8/nyuvvJLevXtXuP8555zDqlWrSlYurF+/Pm+99dZxLZyBAwdy8cUXk5iYSIsWLUrK1Re7+uqr+eCDDzjnnHNc/fuUJcVNJJ/eJNIVOA84F4jF6ff4H87KhCG5zO2AAQO0vCJjxlSZHWnwwmlw2cvu/3PMyoBn+kLiFTB2SmDjq0SrVq2iW7duwQ7DtWnTppGcnMyUKaH/b56dnU39+vU5cuQIp59+Oi+99FLJiKvHHnuMrKwsHn300XLfW97PSURSVLXcTju31Xh74vR3PCkidYAzgcuBJ3BmrBtjyjpRBV5vxcZD/xvg51fgtHugqX+LLJmaZ/LkyaxcuZLc3FwmTZpUkjwuvfRSNmzYwDfffBOwc7ltgeQAHwLXFLc4ROQGVX3Nx+OMAZ7G6YR/RVX/Web1aOANoD+wF7hCVTd71mJ/BeiHkwTfUNUTDk2xFogJuo/ugNWfwB82VVxE0RuHdsDTfZwJhpe9GLDwKlN1b4HUFr62QNxOJFwNzAfe9/xnDs7KhF4TkXBgKs6tsO7ARBHpXma3G4H9qtoReBJnljs4rZ1oVe2Fk1xuEZF2bi7EmCqTmXLiCrzeatASBt0EaTNgd+DXeKgsbj6smqrj5ufjNoGoqr4AzAbmeG5j+fpXMQhYr6obVTUPeA8YW2afscDrnu9nAaM9kxkVqCciEUAdIA846O5SjKkCRw/BrlW+z/+oyPDfQEQd+LZ6zAmJiYlh7969lkRClKqyd+9eYmJifHqf2z6Q/Z6TviEiR4BPcSry+qI1UHq2UAYwuKJ9VLVARLKApjjJZCyw3XPee1R1X9kTiMhkYDJAQkINHrFiQt+2pTgVeAOUQOrFwZBbnSG9p/8OWvQIzHErSXx8PBkZGezevTvYoZgKxMTEEB8f79N7XCUQVR1d6vtZIpJL1c4+HwQUAq2AxsD3IvKVqm4sE+dLwEvg9IFUYXzGHMuXCrzeGnqHU95k3t/hyrcDd9xKEBkZWTLb29QcvtbCOl9EjpsCq6qfqGqcj+fOBNqUeh7v2VbuPp7bVbE4nelX4cxByVfVXcAP2OgvE8oykqFxe6jn/+zfEnWbwNBfOx3z22wOr6l6vvaBXIrT55EuIl+KyGMico2I9PJ0ivviZ6CTiLT31NK6EphTZp85OPW1AMYD36hzE3UrMApAROoBQ3A69o0JTZmLA3f7qrQht0Gdxk4rxJgq5lMC8ZQrGQA8D6wFNuLMAfkJOPH8/eOPVQDcAXwBrAJmqOoKEXlERC727PYq0FRE1gP3AsU1m6cC9UVkBU4iek1VU305vzFV5uA2OLTNv/kfFYmJhWF3wbq5kO5/cTxjfOG2E/0KVS2Zdy8izwE+L9qsqp8Bn5XZ9mCp73NxhuyWfV92eduNCUluK/B6a9BkZ8XCeX+D6z6qnHMYUw63w3gPikjJxylVTQE6ByYkY2qYzGQIi3RWF6wM0fWdWekbv4XNP1TOOYwph9sEciPwuoi8JiJ3iMiLQH4A4zKm5shc7CSPSN/G2Ptk4I1QvyV881dn1UNjqoCrBKKqa3HKiHwOtMDpwzg/gHEZUzMUFTojpCqj/6O0yDow4rewNQk2zqvccxnj4XMCEZGuIjIaiFLVGar6gKo+BQwMeHTGVHe7VzvrmVfGCKyy+k+ChvHwzd+sFWKqhK/zQO4CPsKpe7VcREqXHrFxhMaUVVKBtwoSSES0s+hUZrIzKsuYSuZrC+RmoL+qXgKcATwgInd7XvOzQpwxNVBGsjPUtsmpVXO+Plc7C0/Ns1aIqXy+JpAwzxBaVHUzThI5T0SewBKIMccrrsAb5na8io/CI2HkH2H7MmeGujGVyNff6p0i0qf4iSeZXAjEAZU0RtGYaupoNuxaWTW3r0rrNQGadnRmpxcVVe25Ta3iawJ5GNhZeoOqFqjqdcDpgQrKmBph+zLQoqrpQC8tPALO+JOTvFZ+ULXnNrWKrwlkIvCZiLwnIteLSMviF1TVZjAZU1pxBd7KHsJbnh6XQfPuMO8fUFhQ9ec3tYKvtbBuU9W+OC2RxsA0EVkoIn8XkdNdFFQ0pubKSIZGbZ21O6paWJjTCtm7DtJmVv35Ta3gdiLhalV9UlXH4FTFXYBTm+qnQAZnTLVWWRV4vdXtImiZCPP/CYVWKMIEnt9DQ1Q1R1U/U9U7K1p43Zha59AOOJgRnNtXxUTgzP+D/ZvhsU7wcCN4siekzgheTKZG8SqBiMjZIvJy8Qgsz1KxxpiKVHYFXm8dPQgSBjn7AYWsdPj4LksiJiC8bYH8Cqdc+zUiMgroU2kRGVMTZKZAWASckhjcOL5+xBkJVlp+jrPdGD95m0AOqeoBVf0dcA5W98qYE8tMhhY9nSKHwZSV4dt2Y3zgbQL5tPgbVb0PeKNywjGmBigqhMwqqMDrjdh437Yb4wOvViRU1Y8ARCTBs+mjUt+X54CqHvQ3OGOqpT1rIe9QcEdgFRv9oNPnkZ/zy7bIOs52Y/zk65K2r3uxjwLTsFaKqa2qsgLvySROcL5+9TAczITohnDB479sN8YPPiUQVT2zsgIxpsbISIboWKceVShInOA8XhwJ0Q0seZiA8WseiIjUs9nnxpSRmQyt+1ZdBV5vJQxxWkc2qdAEiK8LSoWJyFUi8qmI7ALWADtEZKWI/EdEQuQjlzFBkncEdgahAq83EoZA/hHYkRrsSEwN4etHpHlAB+BPQEtVjVfVZsBpwI/Av0TkmgDHaEz1sX0ZaGFodKCX1WaI83WrVRwygeFrAjlLVR8FOuN0lgOgqvtU9X1VHQdMD2SAxlQrwazAezINT3GKO25dGOxITA3hazXe4punbwLvlO7/EJEbyuxjTO2TkQyxCVC/ebAjKV/CEEj/yZa7NQHhtpdvNTAfeF9EIj3b7gxMSMZUY5mLIT4EWx/FEoZA9k7YvynYkZgawG0CUVV9AZgNzBGROtia6Ka2y94FWVtD8/ZVMesHMQHkNoHsB1DVN4BXcUqd1A1UUMZUS6FSgfdEmnWFmFjrBzEB4etMdABUdXSp72eJSC7O7HNjaq/MFJBwOKV3sCOpWFgYtBns9IMY4ydf54GUe5tKVT9R1bgT7WNMjZeZDC16QFSIN8YThsDu1XBkX7AjMdWcz/NAROTOsoUURSRKREaJyOvApMCFZ0w1UVTkdKCHcv9HseJ+kPRFwY3DVHu+JpAxQCHwrohs88xA3wisAyYCT6nqtADHaEzo27vOWf0vFCcQltW6H4RFWj+I8ZuvxRRzgeeA5zzDd+OAHFU9UAmxGVN9hFIF3pOJrAOt+sDWH4MdianmXFd7U9V8Vd0O5Jx05wqIyBgRWSMi60XkvnJejxaR6Z7XfxKRdqVeSxSRhSKyQkTSRCTGbRzG+C0jGaIaQFynYEfinYQhsG0x5OcGOxJTjQWiXOhzInK2r2/yzGKfCpwHdAcmikj3MrvdCOxX1Y7Ak8C/PO+NAN4CblXVHsAZgM2AN8FTUoG3mhSnbjMECvNg+9JgR2KqMb8TiKreCCSIyDMiEufDWwcB61V1o6rmAe8BY8vsM5ZfFrGaBYz2jPI6B0hV1WWeGPaqaqFfF2KMW/k5sHNF9bh9VSyheEKh3cYy7vmdQETkXKA90BF4WUQu8fKtrYH0Us8zPNvK3UdVC4AsoCmeYo4i8oWILBaRP7i/AmP8tD0VigqqRwd6sXpxzoJXlkCMHwJxC6sd8F9VPV9VLwWqYtXCCJwS8ld7vl4qIqPL7iQik0UkWUSSd+/eXQVhmVoplCvwnkhxYcWiomBHYqqpQCSQ9qq6EUBEhgMPePm+TKBNqefxnm3l7uPp94gF9uK0Vr5T1T2qegT4DOhX9gSq+pKqDlDVAc2aNfPhkozxQUYyNIyHBi2DHYlv2gyBnH3OEGRjXAhEAvmf55P+1cB5qnrQy/f9DHQSkfYiEgVcCcwps88cfpmYOB74RlUV+ALoJSJ1PYllJLDS7ysxxo3MlNCuwFuRhKHOV7uNZVzyd030V4GLcSYR9lTV+719r6dP4w6cZLAKmKGqK0TkERG52LPbq0BTEVkP3Avc53nvfuAJnCS0FFisqp/6cy3GuHJ4DxzYUv1uXwE07QB14yyBGNdcFVMspqo3ekq59wMGisiLqnqLD+//DOf2U+ltD5b6Phe4vIL3voUzlNeY4KkOFXgrIuLpB7EEYtzxK4EAqGoO8IPnYUztUlyBt1WfYEfiTpvBsPoTZy2TUF1F0YSsQPSBlBCRG0TkjEAe05iQlpkMzbtDVL1gR+KO9YMYPwQ0geAsLBUd4GMaE5qKipwWSOvjBgBWH6f0hogYSyDGFb9vYZWmqrtwOsWNqfn2bYDcrOo1gbCsiChnAID1gxgXAn0L6xbP12o4JMUYH1WnCrwn0mYwbF8GeUeCHYmpZgJ9C2ufiCzAmRRoTM2WkQxR9aFZl2BH4p+EoU4pluKEaIyXAlELK6rU0304i0718ve4xoS8zGRoVY0q8FakzUDnq/WDGB8FogXyl1Lf56pqNvC3ABzXmNCVnws7llfPCYRl1WnsjCSzfhDjo0B0on8hIpOBw0BX4AdPuRFjaq4daVCUX7070EtLGAJps6CosPq3qEyVCXQpE28LKRpTvVXXCrwVaTPEWdN916pgR2KqEdcJxNP3cQfwZ+B+YKeIvBiowIwJaRnJ0KAVNGwV7EgCo2SBqYXBjcNUK64SiIjcDWwH1gFLgL6q+pQvdbCMqdaqawXeijRKgAanOOuDGOMlnxKIiDwtIpOAu4FuqhoPnA50F5FHKiNAY0LO4b2wf1PNuX0FvxRWtJFYxge+tkDmAR2AOCBJRBYD/wE2ABNFpHGA4zMm9NSUCYRltRkCWemQlRHsSEw14VMCUdUPPeXWfwTGAmcB04ACoAnwjYhsCHSQxoSUzBSQMGcOSE1S0g9irRDjHbed6L/GWYvjcZy1QHoCaaraF+gWoNiMCU2ZydCsG0TXD3YkgdWiJ0TWs34Q4zVXCURV1wGDgVlADJAKXOp5LS9g0RkTalSrfwXeioRHOLPSbSSW8ZLriYSeRPGp52FM7bBvI+TsrzkTCMtqMwS++zfkHoSYhsGOxoQ4nxKIiCR4uesBVT3oIh5jQltN7UAvljAEtAgyfoaOo4MdjQlxvrZAXgcUkBPsozgd62+4jMmY0JWRDJF1oVnXYEdSOeIHOAME0n+yBGJOyqcEoqpnVlYgxlQLxRV4wwO6FlvoiG4ALXtZP4jxSqDXAzGm5io46hRRrEkTCMvTZghkpEBhfrAjMSGuUhKIiFxWGcc1Jqh2LIfCvJqfQBKGQP5hJ1kacwKV1QIZX0nHNSZ4iivw1tQRWMWKJxTafBBzEl4lEBE5W0ReFpE+nueTy7z+KxH5g+fxR2wyoamJMpKhfkto2DrYkVSuhq2c4orWD2JOwtuewF8BtwH3i0gToE/pF1X1v6WfFycaY2qUzBTPKKUTDUKsIdoMgU3fORMna8P1Gle8vYV1SFUPqOrvgHOAgSfZf75/YRkTYo7sg30bauYM9PIkDIHsHbB/c7AjMSHM2wRSMttcVe8D3jzRzqr6gj9BGRNyMhc7X2vqBMKyrB/EeMGrBKKqH5XZ1LzsPiJiCymbmiszBZCaV4G3Is26QXSs9YOYE3I7Cqu1iEwsfiIizYGvAhOSMSEoMxmadQlofagX5m8gacOeY7YlbdjDC/NDYEWEsDBoMwi2WgvEVMxtArkFmCwig0RkIPAN8FjgwjImhKg6I7ACfPsqMT6WO95ZUpJEkjbs4Y53lpAYHxvQ87iWMAR2r3L6f4wph6/FFN8AFuOsg/5r4B2cxaQuUdX1gQ/PmBCwfzPk7Av4GujDOsTxr3GJXPfqIoZ3jCMtM4spV/VlWIe4gJ7HteJ+kIyfofO5wY3FhKSTtkBEpLDU02k4hRRvwFlQqh2wH7hGRGzyoKmZKrEC77pdhygoUuav3U3T+lEkxjcK+Dlca9UPwiKtH8RUyJtbWCWDwFX1G1V9UlWvV9V+OGuj3wusx1lgypiaJyMZIupA8+4BPezRgkJenL+RiDBhaIemrNuZzdlPzGftzkMBPY9rUXXhlN62xK2pkDcJRCt8QbVAVdNU9S1V/b2vJxeRMSKyRkTWi8h95bweLSLTPa//JCLtyryeICLZIvI7X89tjNcyk6FVn4BX4H187lqycvL5/bldePfmIdx/QTd2ZOVy4TMLmL04I6Dnci1hiDOEueBosCMxISho1Xg9w36nAucB3YGJIlL2I96NwH5V7Qg8CfyrzOtPAJ9XdqymFivIg+2pAS+gWFSkvJ+SQUKTukw+/VQAbhpxKlOv6kuLhtHcO2MZf5qdRm5+4UmOVMkShkDhUdi2NLhxmJDkKoGIyDoRmS0iD4nI2LItAy8NAtar6kbP8rjvAWPL7DMWZxErcNZfHy3i1FUQkUuATcAKN9dgjFd2Lnf+Aw1wAvl27S72Hs7jnrM7IaVKhZyf2Ip5vzuD287owLuLtjLu+SS27j0S0HP7pE3xhEK7jWWO57YF8iKwA9iL04JYLiJpIvKIiER6eYzWQHqp5xmebeXuo6oFQBbQVETqA38E/nKiE4jIZBFJFpHk3bt3exmWMaUUd6AHuALvS99t5JTYGC5MbHXcaxHhYfxxTFdeuW4A6fuOcMGz3zN3xY6Ant9r9ZtBkw7WD2LK5TaBXKOqt6vqFFW9FTgNmAccxLmtVNkeBp5U1ewT7aSqL6nqAFUd0KxZsyoIy9Q4mSlQrznEtgnYIVMzDvDjxn38anh7IsMr/hM8q3sLPr1rBO2a1mPymyn847NVFBQWBSwOryUMdRKIVtgdamoptwkkS0QSi5+o6lJgpKo+Bgz38hiZQOm/ynjPtnL3EZEIIBan1TMY+LeIbAZ+A/xZRO7w+SqMOZmM5IBX4H3xu400iI7gykEnT0ptmtRl5q1DuWZIAi9+t5GrXv6JnQdzAxaLVxIGO/Ng9qyr2vOakOfPTPTXRORVEblTRKYAxTdqo7w8xs9AJxFpLyJRwJXAnDL7zAEmeb4fD3yjjhGq2k5V2wFPAX9X1Skur8WY46XOgCe6w951sPkH53kApO87wudp27lqcAINYry72xsTGc5fL+nF01f2IS0ziwue+Z6k9XtO/sZASRjqfLV+EFOGqwSiqqtxOsH/h1NYcT1woYjUw+kM9+YYBcAdwBfAKmCGqq7w9KNc7NntVZw+j/U4802OG+prTMClzoCP74KDngbx0SzneQCSyKsLNhEeJtwwvL3P7x3bpzVz7hhOo7pRXPPqT0z5Zh1FRVVwW6lpR6jb1PpBzHFET3JfU0SKcDqzG6rqmiqJqhIMGDBAk5OTgx2GqQ6e7AlZ6cdvj20D9yx3fdj9h/MY9s9vOL/XKTw+obfr4xw+WsCfP0jjo6XbOKNLM56c0IfG9bxt+Lv07lWwezXctbhyz2NCjoikqGq5o0hO2gJR1TBV3Q5EicjlnsdIK99uaqysCibxVbTdS2//tIWc/EJuPt331kdp9aIjeOqKPjx6SU+S1u/lwmcXsGTrfr+OeVIJg50FtbJ3Ve55TLXi9S0sz4zzmao6E2fuxYUiMl5ELvIsc2tMzRAb79t2L+TmFzItaQsjOzeja0v/S8KLCNcOacus25z+iQkvLuT1pM2c7I6CayX9IFbe3fzCbR/IHlX9SFVn4cwE7yki40TkWhHpHNgQjaliox+Esg3syDrOdpc+XJLJnuyj3OKZdR4oifGN+PSu0zi9UzMemrOCO99dQvbRgoCeA3BqYoVHWz+IOUYgSpl0wRm6ez7QDWgQgGMaEzzdx0J4FETWBcTp+7joGUic4OpwRUXKS99vpEerhgzt0DSwsQKN6kbx8nUD+MOYLnyWtp2LpyxgzY4AF2SMiHZm41sCMaX4VR1ORJ4B1gHvqeqmwIRkTJCt+xIKcuDqWdDpbL8P9/XqXWzcfZinr+xzTNmSQAoLE24/oyN92zTmzneXMHbqAv5+aS8u6+f+tttxEgZD0rOQd8Sp1GtqPX9bIA+o6rOWPEyNkjod6sbBqWcG5HAvf7eR1o3qcEGvUwJyvBMZ2qEpn911Gr3jG3HvjGWMnbqAb9cc2/HtetnchKFQVADbbCSWcfiVQFQ1K1CBGBMScg7A2v9Br/EBKd++ZOt+Fm3ex69Oa0/ECcqWBFLzhjG8fdNgbjujA8vSs7jx9WQ+WOKMIPNr2dz4gc5XW2DKeJz0L0REmgFjVfWVKojHmOBa+REU5kHiFQE53EvfbaRhTARXDgxcLS1vFBdkHNC2MXe+s5h7py9j3urdLFi/x/2yuXWbQLNusNVGYhmHNx+JvgX6VnIcxoSG1BnQtBO08v9Xfsvew/xvxQ6uGdKWetGBXYzKW6O7teCLe0bSpF4Uc5ZtY8KAeP/WXE8YDOmLoCgIRR1NyPHmt7oIaCgifwSyyzwOld2mqkFcvMAYPxxIhy0L4Mz7A1I88ZXvNxEZFsb1w9r5H5sf0vcfId9Txff1pC2c3rmZ+ySSMBRSpsHuVdCiR+CCNNWSNy2Q83HW/hgAXATcDDyAs5rgbCAJpxbWduCQiOwUkWsrJ1xjKlGap9ZVr/F+H2rf4TxmpqRzSd9WNG8Y4/fx3Cru83jh2v6M6BRHZLjw67cXk7TBZTHGNoOdr9YPYvCiBaKq6cAJ1zv3rBJY3/MYCrwMvBmIAI2pEqqwbLqzAl8T/0qNALy5cAu5+UXcPCKwEwd9lZqRVdLnER0RxrjnF3LV4DakZmS5a4U0bgf1Wzr9IANvCni8pnoJyI1ZdeonHMJpgfwLmB6I4xpTZXakwp41cIH/66Hl5hfyxsLNjOranE4tgjuv9taRHUq+79+2CSM6xTF3xU6++4PLIcoiTj+ITSg0BGYmelmDVfX2SjiuMZUndQaERUKPS/0+1PuLM9h7OI/JAS5bEgh3j+7Enuw83v5xq/uDJAyFrK2QVXb9N1PbBDyBqOq+QB/TmEpVVAhpM6HTOc5QVT8UFimvfL+JxPhYBrcPvRqjA9o1YXjHprz43QZy8grdHaS4H8QWmKr1qmZmkzGhbNN8yN4Jvf2f+/Hlyp1s2nOYyaefWmllS/x19+jOTivkpy3uDtAyESLr2XwQYwnEGFJnQHQsdDrX70O9/P1G2jSpw5geLQMQWOUY1L4Jwzo05cXvNpKb76IVEh4B8f1tJJaxBGJqubzDsOpj6DEWIv0bbpuyZR8pW/Zz4/CqK1vi1t2jO7H70FHe/sllX0jCUNi5HI4GuOqvqVZc/ZaLyAAR+UBEFotIqoikiUhqoIMzptKt/gzysgNSuuTF+RtpVDeSCVVctsSNwac2ZcipTXhh/gZ3rZA2g0GLIOPnwAdnqg23H5PeBl4DxuFMLrzQ89WY6iV1OjSMh4Rhfh1m4+5svly1k2uHtKVuVHDKlvjq7tGd2X3oKO8uctEKiR8IEmb9ILWc2wSyW1XnqOomVd1S/AhoZMZUtuxdsOEbSLwcwvy75fTKgk1Ehodx3dB2gYmtCgzt0JTB7Zvw/LcuWiExDZ1SJtYPUqu5/at5SEReEZGJInJZ8SOgkRlT2ZbPBi30+/bVnuyjzErJYFy/1jRrEB2g4KrG3Wd1Yteho7znphWSMBQykqGwEpbQNdWC2wRyA9AHGINz66r4NpYx1UfqdGjZC5p38+swbyzcQl5BETeeFnoTB09m6KlNGdSuCc+76QtpMxjyD8POtMoJzoQ8twlkoKoOUNVJqnqD5/GrgEZmTGXas85ZWS/xSr8Ok5NXyJsLN3NWtxZ0bF4/QMFVHRHh7rM6sfPgUWYkp/v25oShzlfrB6m13CaQJBHpHtBIjKlKqTOcTuCe4/w6zKyUdPYfyeeWkdWv9VFsWIemDGzXmOfmbeBogQ+tkNjWENvG+kFqMbcJZAiwTETW2DBeU+2oOrev2o+Ehu7XKS8sUl5ZsIk+bRoxoG3jAAZYtUSEu0d3ZsfBXGb87GsrZAik/+T8m5pax20CGQN0BM7BhvGa6ib9Jziwxe/O8y9W7GDL3iPcEsJlS7w1vGNT+rdtzHPf+tgKaTMYDm13/j1NreM2gUyq4GFM6EudDhF1oJv7cR+qyovfbaRt07qcE8JlS7zltEI6sT0rl5nJGd6/0fpBajW3CeRwqUchcB7QLkAxGVN5CvKc4btdL4Bo92t1/Lx5P8vSD3DTae0JD6verY9iIzrF0TehEc/NW+99K6R5N4huaP0gtZSrBKKqj5d6/A04A6i+vYim9lj/JeQe8Pv21UvfbaBJvSjG9w/9siXeEhF+c1ZntmXlMivFy1ZIWDi0GeTcFjS1TqAqvtUF4gN0LGMqT+p0qBsHHVyuyAes35XNV6t2ce2QttSJCg9gcMF3eqc4+rRpxHPzNpBXUOTdm9oMgV0rIWd/5QZnQo7bYoppntFXqSKyAlgDPBXQyIwJtJwDsOZ/0Gs8hEe6Pswr328kOiKM64a2DVxsIaJ4XkjmgRzeX+xlKyRhiPM1fVHlBWZCks9V38QZbvJroHjYRQGwU1WtnoEJbavmQOFRSJzg+hC7DuUye3Emlw+Ip2n96lW2xFtndG5G7/hYpnyznnH94omKOMnnzNb9ISzCWSe9s/9rqpjqw+cWiKoqMLVUEcVMt8lDRMZ45pKsF5H7ynk9WkSme17/SUTaebafLSIpnpZQioiMcnN+U8ukzoCmHaFVP9eHeCNpC/lFRdw0ouZ2+RX3hWQeyGG2N62QqLpwSm8ngZhaxW0fyGIRGejPiUUkHJiKM4KrOzCxnNntNwL7VbUj8CTwL8/2PcBFqtoLZ/jwm/7EYmqBA+mw+Xun89zlnI3DRwt488ctnNO9Be3j6gU4wNByRpdmJMbHMmXeevILvegLSRjqlIYpOFr5wZmQ4TaBDAYWisgGP2aiDwLWq+pGVc0D3gPGltlnLPC65/tZwGgREVVdoqrbPNtXAHVEpGbeTzCBkTbT+drrcteHmJmcTlZOPpNP7xCgoEJX8byQjP05fLA48+RvaDMYCnJh+7LKD86EDLcr3wTiRmdroHTdhAycxFTuPqpaICJZQFOcFkixccBiVT3uo4+ITAYmAyQkJAQgZFMtFZcuaTMYmrR3dYiCwiJeWbCJ/m0b078aly3xxaiuzenVOpZn563j0n6tiTzRMr3Zu5yvr57t1Mca/aBffU2menDbArm99EJSnsWkbg9kYN4QkR44t7VuKe91VX3JUzV4QLNmzao2OBM6dqTB7tV+/Yf2+fIdZOzPYfLpNbfvo6ziVkj6vhw+WHKCVkjqDPjy/l+eZ6XDx3c5202N5jaBnF3OtvN8PEYmUHoWVrxnW7n7iEgEEAvs9TyPBz4ArlPVDT6e29QmqdMhLBJ6uFvzTFV56buNtI+rx9ndWgQ4uNA2ultzerZuyNR56ymoqC/k60cgP+fYbfk5znZTo/mUQETkNhFJA7qWmgeSKiKbAF9XlfkZ6CQi7UUkCrgSmFNmnzn8UmNrPPCNqqqINAI+Be5T1R98PK+pTYoKIW0WdDoH6jZxdYgfN+4jLTOLm0a0J6yGlC3xlohw16hObNl7hA+Xbit/p6wKRmpVtN3UGL62QN7Bqbr7Eb+sRHgR0F9Vr/blQJ6hv3cAXwCrgBmqukJEHhGRiz27vQo0FZH1wL1A8VDfO3CqAT8oIks9j+Y+XoupDTZ9B9k7/Lp99dJ3G2haL4px/WpnsYWzu7eg+ykNefabdeW3QmIr+HepaLupMXxKIKqapaqbgTwgq1T/h4rIf309uap+pqqdVbWDp6YWqvqgqs7xfJ+rqperakdVHaSqGz3b/6qq9VS1T6nHLl/Pb2qB1BlOsb/OY1y9fe3OQ8xbs5tJw9oRE1mzypZ4q3h2+pa9R/iovFbI6Achsk6ZN4XBmf9XNQGaoHHbB5KoqgeKn6jqfqBvQCIyJlDyjjizz7uPhcgYr9/2wvwNJG1wBvq9/N1GYiLD6NKyAS/Mr71dbed0b0G3Uxoypby+kMQJcNEzzugrBOo2BS2CPWuDEqupOm4TSJiIlIxlFJEmuB8SbEzlWPMZ5GX7XHk3MT6WO95Zwqep2/hwaSand2rGn2ankRgfW0mBhj5nRFZHNu05zMep5bRCEifAPcvh4QPwh43Q7zpY8CRs/LaqQzVVyG0CeRxnIuGjIvIokAT8O3BhGRMAqdOhYWtoO9yntw3rEMeUiX25Z/oy8guVnzbtY8pVfRnWIa6SAq0ezunekq4tG/Ds1+spLDrJErZj/gVxnWH2ZMjeXTUBmirndj2QN4DLgJ2ex2WqauVETOjI3g3rv3Zmnof5/mu+aPM+8jy3aiYNbVvrkwdAWJhw1+hObNxzmI+XVTAiq1hUXbj8NacC8oe3QZGXpeFNteJ6PRBVXamqUzyPlYEMyhi/rZgNWuhq4ahPU7fz1FfriIoI485RHXnrp60lfSK13ZgeLenSogHPfLPu5K2QFj3g3L85i3j9+FzVBGiqlNv1QERErhGRBz3PE0RkUGBDM8YPqdOhZS9oUbY+54mlZWTxm+lLiAgTXr6uP789pwtTrurLHe8ssSRCqVbI7sN8Ul5fSFkDb4KuF8JXD8O2JZUen6lablsgzwFDgYme54dwKusaE3x71kNmis+tj10Hc7n5jWTqRIbz7MS+jOzsTC0a1iGOKVf1JTUjqzKirXbO69mSzi3q88zXXrRCRODiZ6F+c5j1Kzh6qGqCNFXCdTVeVf01kAslw3ijAhaVMf5ImwEI9Bzv9Vty8wu5+c0UDubm897koZzX65RjXh/WIY5bR9b8KrzeCAsT7hzViQ27D/Np2vaTv6FuExj3CuzfDJ/+ttLjM1XHbQLJ96znoQAi0gywXjITfMWVd08dCQ1POfn+OLWu/vh+KsvSD/DkFX3o3qphJQdZ/Z3f6xQ6Na/Ps1+vo+hkrRCAtsNg5B+dn83Sdys/QFMl3CaQZ3AKGbYQkb8BC4C/BywqY9xKX+R80vXh9tVz327go6Xb+P25XTi3R8vKi60GCQ8T7hzdiXW7svlsuRetEIDTf+8Mqf70t85tRlPtuR3G+zbwB5yksQ24RFVnBjIwY1xJnQ4RdZyOWy/8b/kO/vPFGsb2acXtZ9gtKl9c0OsUOjSrxzPetkLCwuGylyEiCmbdYKsX1gBuR2HFAOcDZwGjgDGebcYET0GeM3y36/kQc/LbUCu3HeTeGUvp3aYR/xqXiLhc6ra2CveMyFq7M5vPl+/w7k2xrWHsVNiR6ozMMtWa21tYbwA9cG5lTcFZ09wmEprgWv8V5Oz36vbV7kNHuen1n2kYE8nL1/avtYUS/XVhYitO9aUVAtD1Ahg02ZkbsvaLyg3QVCq3CaSnqt6oqvM8j5txEooxwZM6HerGQYdRJ9ztaEEht76Vwr4jebwyaQDNG1rj2a3wMGe9kDU7D/HFCi9bIQBnPwotejmz1A962YdiQo7bBLJYRIYUPxGRwUByYEIyxoXcLFjzOfQcB+GRFe6mqvxpdhopW/bz+OV96Nm69hZIDJRtWTmc0jCGp0u1QpI27Dlx9eLIGBj/X2flwtk3Owt/mWrHbQLpDySJyGYR2QwsBAaKSJqIpAYsOmO8tXIOFB496e2rl77byOzFmfzmrE5ckOjdMF9zYn3aNOLg0XxW7zjE3JU7SNqwhzveWXLy6sXNOsN5/4bN38OCJ6omWBNQbkuwu1udx5jKkjodmnSA1v0q3OXrVTv55/9Wc0HiKdw9ulMVBlezDesQxwvX9GfSfxfxfx8sp0iVqVf3864AZd9rYOM8mPcPaHc6JAyu/IBNwLhtgQwC9nlWI7wWeBJoWmqFQmOqzoF051Ns4hVO6YxyrNlxiLveXULPVrE8Nr63jbgKsBGdmnFO95bsPZyHCNSP9vKzqQhc+KSz/O37NzrVe0214TaBPKCqh0TkNJyhvK8CzwcuLGN8sHyW8zXx8nJf3pt9lBtf/5l60RG8fN0A6kTZiKtAS9qwh0Wb93Fx71bsP5LPJVN/4PG5a8gr8KJARUys0x9yaDt8fJdTTcBUC24TSHGP1wXAS6r6KVYLywSDKiybDvGDoMmpx72cV1DEbW8vZveho7x03QBaxtqIq0Ar7vOYclVfnpnYl5euGUBEeBjPfrOei6csYHmmF0Uo4wfAqPth5UeQMq3SYzaB4TaBZIrIi8AVwGciEu3HsYxxb+dy2L3KWVK1DFXlgQ+Xs2jTPv49PpE+bRpVfXy1QGpG1jErNp7dowXTbhjIuH6t2Xs4j0um/sBTX60lv+xa6mUNuxtOPRP+dx/sWlUFkRt/uf1PfwLwBXCuqh4AmgC/D1RQxngtdTqERTjDd8v47w+bmZ6czh1ndmRsn9ZBCK52uHVkh+M6zId1iOPxCX348p7TuTDxFJ76ah1jp/zAqu0HKz5QWBhc+iJEN4CZNzhDfE1Ic1sL64iqzlbVdZ7n21V1bmBDM+YkigohbRZ0OscpGV7KvDW7+NunKzm3RwvuPbtzkAI0jepG8dSVfXnhmv7sOpTLxVMW8OzX6ypujTRoAZe+4LQqv/hz1QZrfGa3nUz1tfl7p+O1zO2r9bsOcdc7S+jSsiFPTOhDWJiNuAq2MT1bMveekZzboyWPf7mWy55LYs2OChaX6ngWDLsTkv/r9ImYkGUJxFRfy6ZDdEPo/Mu0pP2H87jx9WSiI8N4ZdIA6nk7nNRUuib1ophyVT+eu7ofmQdyuOjZBTz37XoKymuNjHoQWvWDOXfCga1VH6zxittqvHeKSONAB2OMV1JnwBM9YNk7zm2sVR8DkF9YxO1vL2b7gVxevHYArRvVCXKgpjzn9zqFufeczuhuzfn3/9Yw7oWFrN9VpjUSEQXjX4WiInj/ZigsCE6w5oTctkBaAD+LyAwRGSM2K8tUldQZzlyBgxnO8/zD8PFdaOoMHp6zgoUb9/KPy3rRv619vgllcfWjee7qfjwzsS9b9h7m/GcW8OL8Dceusd7kVGeSYfqPMP+fwQvWVMhtJ/r9QCecCYTXA+tE5O8iYivymMr11cPHj87Jz+HwZw/y9k9buWXkqYzrHx+U0IxvRISLe7di7j2nM7JzM/7x+WoufyGJDbuzf9kp8XLoczV89xhs+i54wZpyue4DUVUFdngeBUBjYJaI/DtAsRnjyM1y+jveuxoOZpa7S92cHZzVrTl/OLdrFQdn/NW8QQwvXdufp67ow4bdhzn/6e955fuNv7RGzvs3NO0AsyfD4b3BDdYcQ9RF2QARuRu4DtgDvAJ8qKr5IhIGrFPVkGuJDBgwQJOTreJ8tXF4L6z51Kmyu/FbKMqHBq3g6EHIyz5u953SjHr3rfa+BpMJSTsP5vLn2Wl8vXoXA9s15j/je9Murh5sXwavnAXNujqLhmVlOPWzRj9Y7iRSEzgikqKqA8p7zW0LpAlwmaqeq6ozVTUfQFWLAO8WozamrEM74edX4PWL4bFOzgicPWtg8C1w41dwzwq+6nAfheHHdo7nEMWPp/7akkcN0KJhDK9MGsDjl/dm9Y5DjHn6O6b9sImiFonQ41JnKdysdECdrx/f5fSLmaBw1QKpjqwFEqIOpDujqFbNga0/AgpNO0H3i6HbxXBK72Mq7CZt2MMnbz3Ng3VmEX1kO9u0Kc/IRMZe+xvvyoebamNHVi73zU7l2zW7Gdy+CS/tvZ7YvONXPTwU3ZIGf1oThAhrhxO1QFx9ZBORe8vZnAWkqOpSN8c0tci+jc6tqVVzIDPF2da8B5xxH3Qf69ymKJU0VJU92Xmk7z/C7kNHKegxnl5LhpJfqNSLCufl6wZY8qiBWsbG8Nr1A5mZnMGjn6ykATuhnPGe9Y/ugBdOc35v4ro4C1XFdXFGcUVYjdfK5LbNP8Dz+Njz/EIgFbhVRGaqqlcd6SIyBngaCAdeUdV/lnk9GngDZwXEvcAVqrrZ89qfgBtxKgPfpapfuLwWEyAvzN/AqLxv6bz8yZJ71Gt73sM3UWdwa7d8J2GsnAM705w3nNIHRj8E3S7mSMN2pO/LYeueI6Sv3czWfUfI2H+ErfuOkL4vh5z8Y5c8rRsVTn5hITee1t6SRw0mIkwY2IbhneLY92wz4gp3HbePRtansG5zwrcuRNJm/vJCWISTROI6O8mlWRfn+7jOEFW33POd8Hd4ZMh17Qad207074DzVTXb87w+8CnOSoUpqtrdi2OEA2uBs4EM4GdgoqquLLXP7UCiqt4qIlcCl6rqFSLSHXgXZ2GrVsBXQGdVrXBhZV9vYVXVL1JVnKeqrmXtl6/S5of7qENeybZ8wimKaUp0rvOHf7BZPzbEjeLnmNNYfqRRSaLYk513zLHqRYXTpkld59G4LglN6tCmSV0SmtQlc38O985cxjWDE3jrp63HVII1NZemzqDwwzuJKMot2XZEo7gv/ybmFJ1GdEQY8XUL6VlnN90ittOBDNoUptPi6BZiczMIK/XfQ1FsAtKsC1KcVJp1hWadWbtg9nG/wzlEkT78n3Q++8aAXEd1+7/lRLew3CaQ1UCv4s5zT0thmap2FZElqtrXi2MMBR5W1XM9z/8EoKr/KLXPF559FopIBM6Q4WbAfaX3Lb1fRefzNYGU95/hUSJZ2f4GYrqe5fVxTiZ39Vd03/Qa0eSXOc/1ATuPc45px51jVcJVRHUcgRQVgRZBURGqBc7sXy0ELUKKCj3fF0KRghYintfQIuf7Iuf5qWv/S2TB8fWNjhLJPwqu4bOCAezCmeAXHia0ahTjSQ51SyWLOiQ0qUuTelHlrhpYeu2JYR3ijntuara1X75K/R/+zinsZTtNSetyN1kdL2Xv4Tz2HT7KvsP5nq95nm15HMkrJJIC2skOOkqm8wjbRifJ5NSw7cSU+hsvIowwji+tkhfZkLVdb4ewcCAMxPMIC3dut0p4qW1hnufOdgkr3j8cwsLIWzef7unvHP83324SMV28+5tXL+ZuH139Fd03H/t37yYZVkYCeQC4FCiudHYRMAd4HGeBqau9OMZ4YIyq3uR5fi0wWFXvKLXPcs8+GZ7nG4DBwMPAj6r6lmf7q8DnqjqrovP53In+ZE/PaA/jL0X4z5AfSxJFQpO6tIyNITLc90GAL8zfQGJ87DHJImnDHlIzsuwWQw3n9sNDbn6hk0yy89h3xEk0e7Od5LI/Owey0ql/aCNNjmzi1rxp5XWz1CyxbeCe5V7vHtBOdE/ZkmnA58Bwz+ZbVbX4f+eTJo+qIiKTgckACQkJvr05K6PczQosH/WGn5H9ouc315X7C6vA8tGBOU/Prys6h7B8zPsln5gk3PM1LBwt+eQU7vmUFY6ECUiE5/Uwz+sRztewcDq8O4LI7OMn+klsPH8YE5gJfuUliWEd4qz1UQuUXbhqWIc4plzVl9SMrBP+/GMiw2ndqM4JaqOVumHy5NflfnDMr9eK9CvmUqQKRYWoFjot9qIi1NMCp6gIpRAtKm7N/9JS16JCT8u+iO6fj0c4/oO7AivOevPk/xBefujvUcHffUX/t7nhcwJRVRWRz1S1F+DPuNhMoE2p5/GebeXtk+G5hRWL05nuzXtR1ZeAl8BpgfgUXWx8ub9IEtuGXqeP9elQJ5TSpuLzjAjQeZIrOkc8vYaMDsw5gLW9f1v+/eOe92Archh/VcWHh7U97yn/d7jP7+ic0OYE7/RBUsX/t/Q87eLAnAMq/LsnNnClftxOJFwsIgP9PPfPQCcRaS8iUcCVOLfBSpsDTPJ8Px74xlNCZQ5wpYhEi0h7nLpci/yM5xhre95DTpll3nOIYm3PewJ5mio5T1VdyzdRZ5A+/J9OExmB2DakD/8n30SdEdDzGFNZquJ3uCb93+J2GO9g4BoR2Qwcxhmdraqa6O0BVLVARO7AWRo3HPivqq4QkUeAZFWdg1Os8U0RWQ/sw0kyePabAazEqcP16xONwHLjm6gzGDX8n8eMYEj3jGAI5KfpqjhPVV2L8wmxA5TqoOvseRhTHVTF73BN+r/FbSd62/K2q+oWvyOqJDYT3RhjfFcZtbC2AiOASZ6koThrhBhjjKkl3CaQ54ChwETP80PA1IBEZIwxplpw3Qeiqv1EZAmAqu73dIQbY4ypJdy2QPI9pUgUQESaQTnTN40xxtRYbjvRrwauwClyOA1niO39qjrzRO8LJhHZDbjt5I/DWTyrJrBrCU015VpqynWAXUuxtqrarLwXXK8HIiJdgeJZaN+o6iqXwYU8EUmuaBRCdWPXEppqyrXUlOsAuxZvuF0PJBrohzMzPAK4XERQ1UcCGZwxxpjQ5bYT/SM8C0gBRwMXjjHGmOrCbQKJV9UxAY0ktL0U7AACyK4lNNWUa6kp1wF2LSflthP9JeBZVU0LfEjGGGOqA7cJZCXQEdiEcwvL51pYxhhjqrdaUwvLGGNMYPk0kVBE/gAliWKQqm4pfgC3VEaAVUlExojIGhFZLyL3lfN6tIhM97z+k4i0C0KYXvHiWk4XkcUiUuBZHTJkeXEt94rIShFJFZGvK/qAE2xeXMetIpImIktFZIGIdA9GnN442bWU2m+ciKiIhOxwWC9+LteLyG7Pz2WpiNwUjDi94c3PRUQmeP5eVojIO36dUFW9fgCLy/u+vOfV7YFTUn4DcCoQBSwDupfZ53bgBc/3VwLTgx23H9fSDkgE3gDGBztmP6/lTKCu5/vbQvHn4uV1NCz1/cXA/4Idt9tr8ezXAPgO+BEYEOy4/fi5XA9MCXasAbqWTsASoLHneXN/zulrKROp4Pvynlc3g4D1qrpRVfOA94CySwKOBV73fD8LGO1Z4jfUnPRaVHWzqqYS+iVovLmWeap6xPP0R5wVKkONN9dxsNTTelDOuqehwZu/FYBHgX8BuVUZnI+8vZbqwJtruRmYqqr7AVR1lz8n9DWBaAXfl/e8umkNlF7/McOzrdx9VLUAZy5M0yqJzjfeXEt14eu13Ah8XqkRuePVdYjIr0VkA/Bv4K4qis1XJ70WEekHtFHVT6syMBe8/f0a57lFOktEArS2bcB5cy2dgc4i8oOI/Cgifk3H8DWB9BaRgyJyCEj0fF/8vJc/gRjjLxG5BhgA/CfYsbilqlNVtQPwR+D+YMfjhoiEAU8Avw12LAHyMdBOnVGmX/LLXYjqKALnNtYZOMtxvCwijdwezKcEoqrhqtpQVRuoaoTn++LnkW6DCBGZQOlPFvGebeXuIyIROKVc9lZJdL7x5lqqC6+uRUTOAv4PuFhVQ7E6gq8/k/eASyozID+c7FoaAD2Bbz3LXg8B5oRoR/pJfy6qurfU79QrOEVkQ5E3v2MZwBxVzVfVTcBanITiitty7jXRz0AnEWnvWdvkSmBOmX3mAJM834/HKSIZirfuvLmW6uKk1yIifYEXcZKHX/d0K5E311H6D/kCYF0VxueLE16LqmapapyqtlPVdjj9UheraiiuKe3Nz+WUUk8vBkK1cKw3f/cf4rQ+EJE4nFtaG12fMdgjB0LpAZyPk5E3AP/n2fYIzi8/QAwwE1gPLAJODXbMflzLQJxPI4dxWlErgh2zH9fyFbATWOp5zAl2zC6v42lgheca5gE9gh2z22sps++3hOgoLC9/Lv/w/FyWeX4uXYMdsx/XIji3F1cCacCV/pzPdTl3Y4wxtZvdwjLGGOOKJRBjjDGuWAIxxhjjiiUQY4wxrlgCMcYY44olEGOMMa5YAjHGGOOKJRBjgkBERovIm8GOwxh/WAIxJjh646zLYEy1ZQnEmODoDSzxrHI5TUT+HqJryxhToYhgB2BMLZUI7AK+AF5R1beCHI8xPrNaWMZUMRGJBPYAW4BbVHVhkEMyxhW7hWVM1euGU3q7ACgEEJEbROQ8cfxXROoENUJjvGAtEGOqmIhci7OIz1vAB8AooD7wK2ATkK6qXwQvQmO8Yy0QY6peb2C5qq7FWbp2BrAV6Af0seRhqgtrgRgTIkTkI+BmDd1VFY05hrVAjAkyEYkVkWeA1y15mOrEWiDGGGNcsRaIMcYYVyyBGGOMccUSiDHGGFcsgRhjjHHFEogxxhhXLIEYY4xxxRKIMcYYVyyBGGOMceX/AW1H9XxTkuMtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from diag_geom import ksq\n",
    "from diag_intgrl import intgrl_thet\n",
    "\n",
    "zz=xr_phi.zz\n",
    "ky=xr_phi.ky\n",
    "kx=xr_phi.kx\n",
    "phi=xr_phi.rephi[-1,:,:,:] + 1j*xr_phi.imphi[-1,:,:,:]\n",
    "print(phi)\n",
    "\n",
    "wrong_ksq = kx**2+ky**2\n",
    "wrong_energy = np.sum(0.5*np.abs(phi)**2*wrong_ksq,axis=0)/len(zz)\n",
    "proper_energy = intgrl_thet(0.5*ksq*np.abs(phi)**2)\n",
    "plt.plot(ky,wrong_energy,\"x-\",label=\"Wrong energy\")\n",
    "plt.plot(ky,proper_energy,\"o-\",label=\"Proper energy\")\n",
    "plt.xlabel(\"$k_y$\")\n",
    "plt.ylabel(\"Energy spectrum $\\int \\sqrt{g} k_\\perp^2|\\phi_{k_x,ky}(z)|^2/2 dz/(\\int \\sqrt{g} dz)$\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}