vct_b.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_levels = '/scratch/e1000/meteoswiss/scratch/swester/input_icon/cosmo_levels.txt'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "vct = pd.read_csv(file_levels, sep=' ', skipinitialspace=True, skiprows=1, names=['level', 'vct_a', 'vct_b'], nrows=81)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x149cab71b0a0>]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA7lUlEQVR4nO3deXxU1f3/8XfWSYAshJANQgh7ZCcRBMQNxSKldanixiJqS9UKYq0i/dXK1xr7bbXWKrjhwldUqiB1oUpcWBQVCGGRNUAgEBJCgKxknTm/PwKjIUFJSHJneT0fj3nYnLl38jlNMvPm3HPP8THGGAEAAFjE1+oCAACAdyOMAAAASxFGAACApQgjAADAUoQRAABgKcIIAACwFGEEAABYijACAAAs5W91AWfD4XDo0KFDCgkJkY+Pj9XlAACAs2CMUUlJieLi4uTre+bxD7cII4cOHVJ8fLzVZQAAgCY4cOCAOnfufMbn3SKMhISESKrtTGhoqMXVAACAs1FcXKz4+Hjn5/iZuEUYOXVpJjQ0lDACAICb+akpFkxgBQAAliKMAAAASxFGAACApQgjAADAUoQRAABgKcIIAACwFGEEAABYijACAAAsRRgBAACWanQYWbVqlcaPH6+4uDj5+Pho6dKlP3nOypUrlZycrKCgIHXr1k3PP/98U2oFAAAeqNFhpKysTAMHDtSzzz57VsdnZWXpqquu0qhRo5SRkaGHH35Y9957rxYvXtzoYgEAgOdp9N40Y8eO1dixY8/6+Oeff15dunTR008/LUlKSkrS+vXr9fe//13XXXddY789AADwMC2+Ud7XX3+tMWPG1Gm78sorNX/+fFVXVysgIKDeOZWVlaqsrHR+XVxc3NJlAgDgFWrsDuUVV+hQYYVyCk/oUGGFDh4v15QRXdU75sd3120pLR5G8vLyFB0dXactOjpaNTU1KigoUGxsbL1zUlNT9eijj7Z0aQAAeISKartyiyp0qLBcOcfLlVNYrqNllSour1FJRbVKKmpUUlGjovJq5ZdUyGHqv8YF3SI8N4xI9bcONsY02H7KrFmzNHPmTOfXxcXFio+Pb7kCAQBwUcYYFZVX6+Dx8tqwUfj9f3MKK5RzvFwFpZU//UI/EOjnq9jwIMWFBatT+2DFhQerZ5Q1QURqhTASExOjvLy8Om35+fny9/dXhw4dGjzHZrPJZrO1dGkAAFjC4TAqKK1UXnGFispPjVxUO0cyCsqqlPOD8HGiyv6Trxkc4KdO7YPVKbw2XHQMsSk0yF+hQQEKDfZXSFCAQoL8FRMapMh2Nvn6NjwgYIUWDyPDhw/XBx98UKdt+fLlSklJaXC+CAAA7q6i2q68ooqToxfldYJFTmG5cgsrVGV3NOo1I9sFqlP4yZGMH4xodDr5CG8TcMYrDq6u0WGktLRUu3fvdn6dlZWljRs3KiIiQl26dNGsWbOUk5OjBQsWSJKmTZumZ599VjNnztSdd96pr7/+WvPnz9dbb73VfL0AAMACJ6pqtDOvRDvySrQzr0Tbc4u1t6BMR0p++rKJr48UFRKk8Da1IxahJ0cuQoIC1L5toDqfHOHo1D5YsWFBCgrwa4UeWaPRYWT9+vW69NJLnV+fmtsxefJkvfbaa8rNzVV2drbz+cTERC1btkz33XefnnvuOcXFxemZZ57htl4AgMurrDk5wnH8+1GNUyMcB46V68DxEzINTAaVvr9sUjt6EVRvVCM6NEgBfiyELkk+xpzp/0bXUVxcrLCwMBUVFSk0NNTqcgAAHqS8yq7M/BLtyC3R7iOldS6rHCmtPGPYOCUqxKbeMSFKig1Vn5gQ9YwKUef27n3ZpLmc7ed3q9xNAwCA1RwOo4PHy7U9r/jkpZVi7cgtUdbRsh8NHEEBvs65GZ1Pm6/RM6qdOrTjhotzRRgBAHicovJqZ+DYnluinScDSNkZ7kqJaBuopNjvRzU6/2ByaETbQK8f4WhphBEAgNuqsTuUVVCm7Xkl2pFbrB0n/3uoqKLB4wP9fNUzul3tZZWYUPWJDVHvmBB1bGcjcFiIMAIAcAs1doe2HirW+v3Hte1QsXbkFSszv1RVNQ3fItspPFh9YkLqzOdIjGwrfyaNuhzCCADAJZVX2ZVx4LjWZR3Xun3HtCH7eIOLf7UN9FOvk4EjKSZEfWJD1Ss6RGHBrGXlLggjAACXUHiiSuv31QaPtfuO6bucIlXb684sDQsOUEpCew2MD1efkwGkU3iwS60misYjjAAALJFbVK61Wce0bt8xrcs6rp2HS+odExMapPMTIzS0a3udnxihXlEhBA8PRBgBALQ4Y4z2HCk7GTxqRz4OHi+vd1y3jm01tGuEzu8aoaGJEercPpiJpV6AMAIAaHY1doe25RY7Rz7W7zuuo2VVdY7x9ZH6xoWdDB7tldI1QpGs2eGVCCMAgGZRVF6tVbuO6LPth7Vi1xEVnqiu87zN31eD4sM1NLF25GNIQnu1s/ExBMIIAOAcZBWU6bPth/XZ9nyt23dMNY7vJ5yGBPkrJaG9hiZ20NDE9urXKUw2f8/d7A1NRxgBAJy1artD6/Yd0+fb8/X5jnztLSir83yPqHYa3SdKo5OiNaRLOGt64KwQRgAAP6ra7tDqzCP6z8ZD+nxHvkoqapzPBfj5aGhihEb3idbopCgldGhrYaVwV4QRAEA9xhhlHCjU0owcfbg5V8d+MPk0om2gLu0dpcuTonRhz0iFBLG4GM4NYQQA4LT3SKmWbjyk/2zM0f6jJ5ztke0CNX5gnH4+IE6D4sPlx1ofaEaEEQDwckdKKvXh5kNampGjTQeLnO3BAX76Wb8Y/XJQnC7sEcn8D7QYwggAeKGyyhot35anpRmH9OXuAtlP3gXj5+ujUT0jdc3gTrrivGi1CeRjAi2P3zIA8BI1dodW7y7Q0owcLd96WOXV3286NzA+XNcMitPPB8ax8BhaHWEEADzc3iOleuObbL2/KUcFpd9PRO3aoY2uHtxJvxzUSYmR3AUD6xBGAMAD1dgd+nxHvv7vm/1anVngbO/QtnYi6tWDO2lg5zD2fYFLIIwAgAc5XlalN9dm681vs5VTWLsRnY+PNLpPlG4e1kWjenZUABNR4WIIIwDgAfKKKvTS6r1689ts51yQ9m0CNOH8LrplWBfFR7SxuELgzAgjAODG9hWU6YVVe7Q4PUdVdockqW9cqKaOTNS4AbEKCmAvGLg+wggAuKEdecWat2KPPth0SKf2phuaGKG7L+2hi3pGMhcEboUwAgBuJCP7uJ77Yo8+3X7Y2XZp746669IeOr9rhIWVAU1HGAEAF2eM0Zo9R/XcF7u1Zs9RSbWTUq/qF6vfXtJd/TqFWVwhcG4IIwDgohwOo0+3H9ZzK/Zo04FCSZK/r4+uGdxJ0y7pru4d21lbINBMCCMA4GJq7A59uDlXc1fs1q7DpZIkm7+vbhraRXde1E2dwoMtrhBoXoQRAHARlTV2LU7P0fMr9yj7WO2OuSE2f906PEFTRyaqYwjLtMMzEUYAwGIOh9F/NuXoyeW7dPB47UJlEW0DNXVkV00c3lVhwQEWVwi0LMIIAFjEGKMVu47ofz/eqe25xZKkqBCbpl3cXTcOjWfHXHgNftMBwAIbDxTqif9u1zd7j0mSQoL8Ne3i7po6MlHBgSxUBu9CGAGAVrTnSKn+9vFOfbw1T5IU6O+rycMTdNclPdS+baDF1QHWIIwAQCvIK6rQPz/bpX+vPyi7w8jXR7p2SGfdd0Uv7o6B1yOMAEALKiqv1vMr9+jVr7JUUV27d8zlSdH6w896q1d0iMXVAa6BMAIALcDhMHo3/aD++vEOHS2rkiSlJLTXQ2P7KIVl24E6CCMA0Mw2HyzUn/6zVRtPrpraI6qdHvpZH41OimIDO6ABhBEAaCbHy6r0t+U79dbabBkjtbP5a8blPTV5RFcF+PlaXR7gsggjAHCOjDFavCFHj320TYUnqiVJVw+K08NXJSkqNMji6gDXRxgBgHOQU1iuh5ds0cpdRyRJfWJC9Ogv+mpYtw4WVwa4D8IIADSBw2H01rpspS7bodLKGgX6+2rG5T1156huXJIBGokwAgCNlH30hB5cvFlf7z0qSRrSJVz/+6uB6hHVzuLKAPdEGAGAs2R3GL2+Zp/+9slOlVfbFRzgpweu7K3JI7rKz5e7ZICmIowAwFnYnV+qBxdvVvr+45KkC7pF6K/XDVBCh7YWVwa4P8IIAPyIGrtDL67eq6c/zVRVjUPtbP56+Kok3Xh+vHwZDQGaBWEEAM5gZ16Jfv/OJm3JKZIkXdK7ox6/pr/i2EsGaFaEEQA4jTFGr63Zp9T/7lBVjUOhQf56ZHxfXTukEyuoAi2AMAIAP3CkpFIPvLtJK3bWrhtyWZ8oPXFtfxYvA1oQYQQATvpiR74eeHeTCkqrZPP31exxSZp4QQKjIUALI4wA8HoV1XY98d8dem3NPkm1q6g+c9Ng9YoOsbYwwEsQRgB4tYPHT+iuhRu0+WDtJNUpI7rqobF9FBTgZ3FlgPcgjADwWit3HdH0tzNUeKJa4W0C9I8bBunSPlFWlwV4HcIIAK/jcBj96/PdevqzXTJGGtA5THNvGaLO7dtYXRrglQgjALzK8bIqzVi00bnL7i3DuuhP48+TzZ/LMoBVCCMAvMaWg0Wa9ka6cgrLZfP31ePX9Nd1yZ2tLgvweoQRAB7PGKO31x3QI//Zqiq7Qwkd2mjeLck6Ly7U6tIAiDACwMNVVNv1x6Xf6d30g5Kky5Oi9eQNAxUWHGBxZQBO8W3KSXPnzlViYqKCgoKUnJys1atX/+jxCxcu1MCBA9WmTRvFxsbqtttu09GjR5tUMACcrf1Hy3TN3DV6N/2gfH2kP/yst16cmEwQAVxMo8PIokWLNGPGDM2ePVsZGRkaNWqUxo4dq+zs7AaP//LLLzVp0iTdfvvt2rp1q9555x2tW7dOd9xxxzkXDwBn8sWOfP38X19qe26xOrQN1Bu3D9Ndl/Rgp13ABTU6jDz11FO6/fbbdccddygpKUlPP/204uPjNW/evAaP/+abb9S1a1fde++9SkxM1IUXXqjf/OY3Wr9+/TkXDwCnM8bo5dV7NfX1dSqpqNGQLuH66N5RGtEj0urSAJxBo8JIVVWV0tPTNWbMmDrtY8aM0Zo1axo8Z8SIETp48KCWLVsmY4wOHz6sd999V+PGjTvj96msrFRxcXGdBwD8lKoahx5cvFmPfbRdxkg3DY3X278erpgwNrkDXFmjwkhBQYHsdruio6PrtEdHRysvL6/Bc0aMGKGFCxdqwoQJCgwMVExMjMLDw/Wvf/3rjN8nNTVVYWFhzkd8fHxjygTghY6VVenW+d/q3+tr54f86efn6fFr+ivQv0lT4wC0oib9lZ6+g6Ux5oy7Wm7btk333nuv/vSnPyk9PV0ff/yxsrKyNG3atDO+/qxZs1RUVOR8HDhwoCllAvASuw6X6JfPfam1WccUYvPXK1PO19QLE9ltF3ATjbq1NzIyUn5+fvVGQfLz8+uNlpySmpqqkSNH6oEHHpAkDRgwQG3bttWoUaP02GOPKTY2tt45NptNNputMaUB8FJfZhZo2hvpKq2sUZeINnplSop6RLHbLuBOGjUyEhgYqOTkZKWlpdVpT0tL04gRIxo858SJE/L1rftt/Pxql102xjTm2wNAHf/ZmKPbXlur0soaDU2M0H/uHkkQAdxQoxc9mzlzpiZOnKiUlBQNHz5cL774orKzs52XXWbNmqWcnBwtWLBAkjR+/Hjdeeedmjdvnq688krl5uZqxowZGjp0qOLi4pq3NwC8gjFGL63eq8eX7ZAk/XxArJ68YSD7ywBuqtFhZMKECTp69KjmzJmj3Nxc9evXT8uWLVNCQoIkKTc3t86aI1OmTFFJSYmeffZZ3X///QoPD9dll12mv/71r83XCwBew+Eweuyj7XrlqyxJ0u0XJmr2VUmsHwK4MR/jBtdKiouLFRYWpqKiIoWGspcE4K0qa+y6/9+b9OHmXEnS7KuSdOdF3SyuCsCZnO3nN3vTAHALZZU1unPBeq3Zc1QBfj76268G6urBnawuC0AzIIwAcHlF5dW67dW12pBdqLaBfnphYoou7MmKqoCnIIwAcGlHSys16ZW12nqoWGHBAXp96lANig+3uiwAzYgwAsBlHS6u0K0vf6vM/FJFtgvU/90+TEmxzBsDPA1hBIBLOnj8hG55+VvtP3pCMaFBWnjnMHXv2M7qsgC0AMIIAJez90ipbn35Wx0qqlB8RLDevOMCxUe0sbosAC2EMALApezIK9atL69VQWmlundsq4V3XMCuu4CHI4wAcBmbDxZq0itrVXiiWkmxofq/24cqsh37VAGejjACwCWs23dMt726TqWVNRoUH67XbxuqsDYBVpcFoBUQRgBYbnXmEd25YL0qqh0alhih+VPOVzsbb0+At+CvHYClPt12WHct3KAqu0MX9+qo529NVnAgG94B3oQwAsAyH3+Xp3ve3KAah9GVfaP1zE2D2XkX8EKEEQCWSNt22BlEfjEwTk/dMFD+fr5WlwXAAvzlA2h1n20/rLsWpqvGYTSeIAJ4Pf76AbSqL3bk67dvbFC13WjcgFj9gyACeD3eAQC0mpW7jug3b6Sryu7Q2H4xenrCIIIIAMIIgNbxZWaBfr1gvapqHM7JqgEEEQAijABoBen7j+nOBetVWePQ5UnR+tdNQwgiAJx4NwDQorbnFuu2V9epvNqui3p11HO3DFagP289AL7HOwKAFrOvoEwT569VcUWNUhLa64Vbk1lHBEA9hBEALSKvqEK3zv9WBaWVSooN1fwp57OyKoAGEUYANLvjZVWaOP9bHTxerq4d2mjB1KEKC2bTOwANI4wAaFallTWa8to6ZeaXKiY0SP93+zB1DLFZXRYAF0YYAdBsqu0O3bVwgzYdKFT7NgF6446hio9oY3VZAFwcYQRAszDGaNaSLVq164iCA/z02m1D1SMqxOqyALgBwgiAZvGPTzP1bvpB+fpIz90yWAPjw60uCYCbIIwAOGdvr83WM59lSpL+ck1/XdYn2uKKALgTwgiAc/LFjnzNXvqdJOney3ropqFdLK4IgLshjABoss0HC3XXwg2yO4yuG9JZ913Ry+qSALghwgiAJjl4/ISmvrZe5dV2jeoZqSeu6y8fHx+rywLghggjABqttLJGd7y+3rm66txb2PgOQNPx7gGgUewOo+lvZWhHXok6htg0f3KKQoJYXRVA0xFGADTKXz/eoc925Mvm76uXJqUoLjzY6pIAuDnCCICztmhdtl5ctVeS9LfrB2oQa4kAaAaEEQBn5Zu9RzX7vdpbeKeP7qlfDIyzuCIAnoIwAuAn7T9apmlvpKvGYTRuQKymj+5pdUkAPAhhBMCPKiqv1tTX1qnwRLUGdg7Tk9cPlK8vt/ACaD6EEQBnVGN36J43N2jPkTLFhAbppUkpCgrws7osAB6GMALgjOZ8uE2rMwsUHOCnlyenKCo0yOqSAHggwgiABi34ep8WfL1fkvSPCYPUr1OYxRUB8FSEEQD1rM48okc/2CZJ+sPPeutn/WIsrgiAJyOMAKhjd36pc/O7a4d00m8v7m51SQA8HGEEgNPxsird/vo6lVTUKCWhvVKvZfM7AC2PMAJAkuRwGN33743af/SEOrcP1gsTk2Xz584ZAC2PMAJAkjRv5R6t2HlENn9fvTAxWR3a2awuCYCXIIwA0Jo9BXpy+U5J0pxf9lXfOO6cAdB6CCOAl8svrtC9b22Uw0jXDemsG1LirS4JgJchjABerMbu0O/eylBBaaX6xITosav7MWEVQKsjjABe7Km0Xfo265jaBvrpuVuGKDiQCasAWh9hBPBSn20/rLkr9kiS/vqrAeresZ3FFQHwVoQRwAsdKizX/e9skiRNGdFVPx8QZ3FFALwZYQTwMjV2h2a8vVGFJ6o1oHOYHr4qyeqSAHg5wgjgZZ79YrfW7jumdjZ//eumwQr0520AgLV4FwK8yLd7j+qZzzIlSX+5pp8SOrS1uCIAIIwAXuN4WZVmLPp+PZFfDupkdUkAIIkwAngFY4z+sHizcosq1C2yreb8sq/VJQGAU5PCyNy5c5WYmKigoCAlJydr9erVP3p8ZWWlZs+erYSEBNlsNnXv3l2vvPJKkwoG0HhvfLNfadsOK9DPV8/cNFhtbf5WlwQATo1+R1q0aJFmzJihuXPnauTIkXrhhRc0duxYbdu2TV26dGnwnBtuuEGHDx/W/Pnz1aNHD+Xn56umpuaciwfw07bnFut/PtouSXpwbB/168S+MwBci48xxjTmhGHDhmnIkCGaN2+esy0pKUlXX321UlNT6x3/8ccf68Ybb9TevXsVERHRpCKLi4sVFhamoqIihYaGNuk1AG90oqpGv3j2K+3OL9VlfaI0f3IKy70DaDVn+/ndqMs0VVVVSk9P15gxY+q0jxkzRmvWrGnwnPfff18pKSn63//9X3Xq1Em9evXS73//e5WXl5/x+1RWVqq4uLjOA0Djzflgm3bnlyoqxKa//WoAQQSAS2rUZZqCggLZ7XZFR0fXaY+OjlZeXl6D5+zdu1dffvmlgoKC9N5776mgoEB33XWXjh07dsZ5I6mpqXr00UcbUxqA03y4+ZDeXndAPj7S0xMGqUM7m9UlAUCDmjSB9fR/XRljzvgvLofDIR8fHy1cuFBDhw7VVVddpaeeekqvvfbaGUdHZs2apaKiIufjwIEDTSkT8FoHjp3QrCVbJEl3XdJdI3pEWlwRAJxZo0ZGIiMj5efnV28UJD8/v95oySmxsbHq1KmTwsK+nzSXlJQkY4wOHjyonj171jvHZrPJZuNfcUBTVNsduvftDJVU1GhIl3DNuLyX1SUBwI9q1MhIYGCgkpOTlZaWVqc9LS1NI0aMaPCckSNH6tChQyotLXW27dq1S76+vurcuXMTSgbwY/6RtksZ2YUKCfLXP28crAA/lhMC4Noa/S41c+ZMvfzyy3rllVe0fft23XfffcrOzta0adMk1V5imTRpkvP4m2++WR06dNBtt92mbdu2adWqVXrggQc0depUBQcHN19PAOir3QWat3KPJOmJawcoPqKNxRUBwE9r9DojEyZM0NGjRzVnzhzl5uaqX79+WrZsmRISEiRJubm5ys7Odh7frl07paWl6Xe/+51SUlLUoUMH3XDDDXrsscearxcAdKysSvct2ihjpJuGxmvcgFirSwKAs9LodUaswDojwE+7+80N+mhzrnpEtdMH91yo4EA/q0sC4OVaZJ0RAK7pw82H9NHmXPn5+ugfNwwiiABwK4QRwM0dKanU/1v6nSTp7ku6q39nlnsH4F4II4AbM8bo/y39TsdPVKtPTIjuuaz+rfIA4OoII4Ab+2Bzrj7emid/Xx89ecNABfrzJw3A/fDOBbip/JIK/ek/tZdnfndZT/WN4/IMAPdEGAHckDFGs9/7ToUnqtU3LlR3Xdrd6pIAoMkII4AbWroxR2nbDivAz0d/v34gq6wCcGu8gwFu5nBxhf78/jZJ0vTRPZUUy9o7ANwbYQRwI8YYPbxki4rKq9W/U5imXczlGQDujzACuJHFG3L02Y58Bfr56skbBsqfyzMAPADvZICbyCuq0KMfbJUkzbiip3pFh1hcEQA0D8II4AaMMXpoyWaVVNRoYHy4fj2qm9UlAUCzIYwAbuCd9Qe1YucRBfr76snrB3B5BoBH4R0NcHE5heX6nw9r7565/4pe6hHF5RkAnoUwArgwY4weWrxZJZU1GtwlXHdweQaAByKMAC5s0boDWp1ZIJu/r/5+/UD5+fpYXRIANDvCCOCijpRU6vFl2yVJvx/TW907trO4IgBoGYQRwEX95aNtKq6oUd+4UN02sqvV5QBAiyGMAC7oq90FWrrxkHx8pMev6c/dMwA8Gu9wgIupqLbrj0u/kyRNuiBBA+PDrS0IAFoYYQRwMXNX7FFWQZmiQmy6/8reVpcDAC2OMAK4kD1HSvX8ij2SpEfG91VoUIDFFQFAyyOMAC7CGKPZ721Rld2hS3p31FX9Y6wuCQBaBWEEcBFLNuTom73HFBTgq//5ZT/5+LCmCADvQBgBXMDxsir95eSaIveO7qn4iDYWVwQArYcwAriAJ/67Q8fKqtQrup3uZMl3AF6GMAJYbG3WMS1af0BS7ZoiAawpAsDL8K4HWKiqxqHZ722RJN14frxSukZYXBEAtD7CCGChl1bvVWZ+qTq0DdRDY/tYXQ4AWIIwAlgk++gJPfNZpiRp9rgkhbcJtLgiALAGYQSwgDFGf/zPd6qscWhE9w66ZnAnq0sCAMsQRgALfLg5V6t2HVGgn6/+52rWFAHg3QgjQCsrKq/WnA+3SZJ+e0l3de/YzuKKAMBahBGglT25fKeOlFQqMbKtfntJd6vLAQDLEUaAVrQ9t1hvfLNfkvTY1f0UFOBncUUAYD3CCNBKjDGa88E2OYw0rn+sRvaItLokAHAJhBGglXyyNU9f7z0qm78va4oAwA8QRoBWUFFt12Mf1W6E95uLurERHgD8AGEEaAUvr96rg8fLFRMapGlMWgWAOggjQAvLK6rQc1/skSTNuqqP2gT6W1wRALgWwgjQwv768Q6VV9uVnNBevxgYZ3U5AOByCCNAC0rff1zvZeRIkh4Zfx4rrQJAAwgjQAtxOIzmfLBVknR9cmcN6BxubUEA4KIII0ALWZKRo00Hi9TO5q8Hftbb6nIAwGURRoAWcKKqRn/7ZIck6e5LeygqJMjiigDAdRFGgBbw0qosHS6uVKfwYN02sqvV5QCASyOMAM0sv7hCL6yqvZX3wbF92H8GAH4CYQRoZk8u36UTVXYN7hKu8QNirS4HAFweYQRoRttzi/Xv9AOSpD+OS+JWXgA4C4QRoJkYY/T4su0yJ3flTU6IsLokAHALhBGgmazYdUSrMwsU6OerB3/GrrwAcLYII0AzqLE79PjJXXknj0hQlw7sygsAZ4swAjSDResPKDO/VOFtAnTPpT2tLgcA3AphBDhHpZU1+kfaLknS9NE9FdYmwOKKAMC9EEaAc/Tiqr0qKK1SYmRb3TIswepyAMDtEEaAc5BfUqGXV++VJD1wZW8F+vMnBQCN1aR3zrlz5yoxMVFBQUFKTk7W6tWrz+q8r776Sv7+/ho0aFBTvi3gcp75LFMnquwaGB+usf1irC4HANxSo8PIokWLNGPGDM2ePVsZGRkaNWqUxo4dq+zs7B89r6ioSJMmTdLo0aObXCzgSvYeKdVba2sXOJs1tg8LnAFAEzU6jDz11FO6/fbbdccddygpKUlPP/204uPjNW/evB897ze/+Y1uvvlmDR8+vMnFAq7k78t3yu4wuqxPlC7o1sHqcgDAbTUqjFRVVSk9PV1jxoyp0z5mzBitWbPmjOe9+uqr2rNnjx555JGz+j6VlZUqLi6u8wBcSUb2cS3bkicfH7HAGQCco0aFkYKCAtntdkVHR9dpj46OVl5eXoPnZGZm6qGHHtLChQvl7+9/Vt8nNTVVYWFhzkd8fHxjygRalDFGT/x3hyTpuiGd1TsmxOKKAMC9NWkC6+nXxo0xDV4vt9vtuvnmm/Xoo4+qV69eZ/36s2bNUlFRkfNx4MCBppQJtIgVO4/o26xjCvT31cwrzv73GgDQsLMbqjgpMjJSfn5+9UZB8vPz642WSFJJSYnWr1+vjIwM3XPPPZIkh8MhY4z8/f21fPlyXXbZZfXOs9lsstlsjSkNaBV2x/ejIreN6Kq48GCLKwIA99eokZHAwEAlJycrLS2tTntaWppGjBhR7/jQ0FBt2bJFGzdudD6mTZum3r17a+PGjRo2bNi5VQ+0svcycrTzcIlCg/x11yU9rC4HADxCo0ZGJGnmzJmaOHGiUlJSNHz4cL344ovKzs7WtGnTJNVeYsnJydGCBQvk6+urfv361Tk/KipKQUFB9doBV1dZY3cu+373pT1Y9h0Amkmjw8iECRN09OhRzZkzR7m5uerXr5+WLVumhITaZbBzc3N/cs0RwB29+W22cgrLFRMapMkjulpdDgB4DB9jjLG6iJ9SXFyssLAwFRUVKTQ01Opy4IXKKmt08d++UEFplR6/pr9uHtbF6pIAwOWd7ec3G2kAZ+G1NftUUFqlhA5tdH1KZ6vLAQCPQhgBfkLRiWo9v3KPJGnmFb0U4MefDQA0J95VgZ/wwqo9KqmoUZ+YEI0fEGd1OQDgcQgjwI/IL6nQq1/tkyTdP6a3fH3ZDA8AmhthBPgRz32+W+XVdg3uEq7Lk6KsLgcAPBJhBDiDA8dO6M21tbepP3Bl7wa3PAAAnDvCCHAG//wsU9V2owt7RGpE90irywEAj0UYARqwO79USzYclCT9/sreFlcDAJ6NMAI04F+fZ8phpCvOi9ag+HCrywEAj0YYAU6z50ipPth0SJI04/KeFlcDAJ6PMAKc5rnPd8thpMuTotU3LszqcgDA4xFGgB/IKijT0o05kqTpoxkVAYDWQBgBfuC5L2pHRS7rE6X+nRkVAYDWQBgBTso+ekLvZdSOivzush4WVwMA3oMwApz03Be7ZXcYXdSrowZ3aW91OQDgNQgjgGpXW118cl0R5ooAQOsijACS5q3coxpH7WqryQmMigBAayKMwOvlFJbrnfUHJEn3MioCAK2OMAKv9/yKPaq2Gw3v1kFDEyOsLgcAvA5hBF4tt6hci9YxKgIAViKMwKu9sHKvquwOnd+1vS7oxqgIAFiBMAKvlV9coTfXZkuSpo/uJR8fH4srAgDvRBiB13ph1V5V1TiUnNBeI3t0sLocAPBahBF4pSMllVr47X5JtXNFGBUBAOsQRuCVXlq9VxXVDg2MD9dFPSOtLgcAvBphBF7naGml/u/r2lGRGYyKAIDlCCPwOi+tzlJ5tV39O4Xpkt4drS4HALweYQRe5VhZlRZ8vU8Sc0UAwFUQRuBVXvkySyeq7DovNlSXJ0VZXQ4AQIQReJHCE1V6bc0+SYyKAIArIYzAa7z61T6VVtaoT0yIxpwXbXU5AICTCCPwCmWVNc5RkXsu6yFfX0ZFAMBVEEbgFd5ed0BF5dVKjGyrsf1irS4HAPADhBF4vKoah15evVeSdOeobvJjVAQAXAphBB7v/U2HlFtUoch2Nl07pJPV5QAATkMYgUdzOIxeWLlHknT7hYkKCvCzuCIAwOkII/Bon+/IV2Z+qUJs/rrlgi5WlwMAaABhBB5t3slRkZsv6KLQoACLqwEANIQwAo+1bt8xpe8/rkA/X90+MtHqcgAAZ0AYgcd6fkXtqMh1yZ0UFRpkcTUAgDMhjMAj7cwr0Wc78uXjU3s7LwDAdRFG4JFO3UHzs74x6taxncXVAAB+DGEEHiensFzvbzokSZp2cXeLqwEA/BTCCDzOK19mqcZhNLxbBw2MD7e6HADATyCMwKMUlVfr7bXZkqRfX8xcEQBwB4QReJS31marrMquXtHtdEmvjlaXAwA4C4QReIyqGode/SpLUu0dND4+bIgHAO6AMAKP8cGmQzpcXKmoEJt+MSjO6nIAAGeJMAKPYIzRS6v3SpKmjOwqmz8b4gGAuyCMwCOszizQjrwStQn00y1DE6wuBwDQCIQReIRToyITzo9XWBs2xAMAd0IYgdvbdqhYqzML5OsjTWVDPABwO4QRuL2XT46KXNU/VvERbSyuBgDQWIQRuLXcou+Xfv/1RSxyBgDuiDACt/baV/tU4zAalhihAZ3DrS4HANAEhBG4rbLKGr15cun3O0cxKgIA7qpJYWTu3LlKTExUUFCQkpOTtXr16jMeu2TJEl1xxRXq2LGjQkNDNXz4cH3yySdNLhg4ZcmGgyqpqFFiZFtd1ifK6nIAAE3U6DCyaNEizZgxQ7Nnz1ZGRoZGjRqlsWPHKjs7u8HjV61apSuuuELLli1Tenq6Lr30Uo0fP14ZGRnnXDy8l8Nh9NqafZKkycMT5OvL0u8A4K58jDGmMScMGzZMQ4YM0bx585xtSUlJuvrqq5WamnpWr9G3b19NmDBBf/rTn87q+OLiYoWFhamoqEihoaGNKRceatWuI5r0ylq1s/nr61mXKSSItUUAwNWc7ed3o0ZGqqqqlJ6erjFjxtRpHzNmjNasWXNWr+FwOFRSUqKIiIgzHlNZWani4uI6D+CHTo2K/Cq5M0EEANxco8JIQUGB7Ha7oqOj67RHR0crLy/vrF7jySefVFlZmW644YYzHpOamqqwsDDnIz4+vjFlwsNlFZTp8x358vGRJo/oanU5AIBz1KQJrKdvzW6MOavt2t966y39+c9/1qJFixQVdeYJh7NmzVJRUZHzceDAgaaUCQ/1+slRkUt7Rykxsq21xQAAzpl/Yw6OjIyUn59fvVGQ/Pz8eqMlp1u0aJFuv/12vfPOO7r88st/9FibzSabzdaY0uAlSiqq9W76QUnSFEZFAMAjNGpkJDAwUMnJyUpLS6vTnpaWphEjRpzxvLfeektTpkzRm2++qXHjxjWtUkDS4vSDKq2sUfeObTWqZ6TV5QAAmkGjRkYkaebMmZo4caJSUlI0fPhwvfjii8rOzta0adMk1V5iycnJ0YIFCyTVBpFJkybpn//8py644ALnqEpwcLDCwsKasSvwdA6H0etf75dUOypyNpcGAQCur9FhZMKECTp69KjmzJmj3Nxc9evXT8uWLVNCQoIkKTc3t86aIy+88IJqamp099136+6773a2T548Wa+99tq59wBeY2XmEWUVlCnE5q9rh3S2uhwAQDNp9DojVmCdEUjS5FfWauWuI7r9wkT9v5+fZ3U5AICf0CLrjABW2XOkVCt3HZGPjzRpeILV5QAAmhFhBG5hwcnbeUf3iVJCB27nBQBPQhiByyuucztvosXVAACaG2EELu+d9QdVVmVXz6h2Gtmjg9XlAACaGWEELs3uMM4VV6eM5HZeAPBEhBG4tBU785V97IRCg/x1zeBOVpcDAGgBhBG4tFO78944tIvaBDZ6WRwAgBsgjMBlZR4u0erMAvn6SBMv4HZeAPBUhBG4rNe/3idJuuK8aMVHtLG2GABAiyGMwCUVlVdrcXqOJG7nBQBPRxiBS3pn/QGVV9vVJyZEF3SLsLocAEALIozA5dgdxnmJht15AcDzEUbgcj7fka8Dx8oV3iZAvxzE7bwA4OkII3A5r63JkiTdeH4XBQf6WVwNAKClEUbgUjIPl+ir3Ufl5+ujiezOCwBegTACl/Lm2mxJtbvzdgoPtrgaAEBrIIzAZVRU27VkQ+3tvDcP62JxNQCA1kIYgctYtiVXReXV6hQerIt6drS6HABAKyGMwGW8dfISzU1D4+Xry+28AOAtCCNwCZmHS7Ru33H5+fro+pR4q8sBALQiwghcwltrD0iqnbgaHRpkcTUAgNZEGIHlKqrtWrzhoCTpJiauAoDXIYzAcv/9jomrAODNCCOw3Fvf1l6iufH8ePkxcRUAvA5hBJbKPFyitfuOMXEVALwYYQSWOjVx9bI+UYoJY+IqAHgjwggs88OJqzcPZeIqAHgrwggsU2fiai8mrgKAtyKMwDKnJq5OYOIqAHg1wggssTu/duKqr490fUpnq8sBAFiIMAJLfD9xNVqxYcEWVwMAsBJhBK2uzsTVYdzOCwDejjCCVvfxd3kqPFGtuLAgXdwryupyAAAWI4yg1b25NluSNOH8LkxcBQAQRtC6dueXam1W7cTVG85n4ioAgDCCVvb2yVERJq4CAE4hjKDVMHEVANAQwghazSdb83SciasAgNMQRtBq3vyWiasAgPoII2gVe46U6lsmrgIAGkAYQav4fuJqFBNXAQB1EEbQ4ipr7Ho3/dTE1S4WVwMAcDWEEbS4T7Ye1vET1Ypl4ioAoAGEEbS4f6+r3RTvhpR4Jq4CAOohjKBFHTx+Ql/tKZAk/SqZiasAgPoII2hRi9NzZIw0skcHxUe0sbocAIALIoygxTgcRu+k116iuT6ZFVcBAA0jjKDFfJN1VAePlyskyF8/6xdjdTkAABdFGEGLeWd97e28vxgYp6AAP4urAQC4KsIIWkRxRbWWbcmVJF2fwiUaAMCZEUbQIj7clKvKGod6RbfTwM5hVpcDAHBhhBG0iH+v/37iqo8Pa4sAAM6MMIJml3m4RBsPFMrf10dXD+5kdTkAABdHGEGze+fkPjSX9YlSxxCbxdUAAFwdYQTNqtru0JINtWGEiasAgLNBGEGz+mJHvgpKqxTZzqZLene0uhwAgBtoUhiZO3euEhMTFRQUpOTkZK1evfpHj1+5cqWSk5MVFBSkbt266fnnn29SsXB9py7RXDukkwL8yLoAgJ/W6E+LRYsWacaMGZo9e7YyMjI0atQojR07VtnZ2Q0en5WVpauuukqjRo1SRkaGHn74Yd17771avHjxORcP17HpQKF+vWC90rYdliRdz6Z4AICz5GOMMY05YdiwYRoyZIjmzZvnbEtKStLVV1+t1NTUesc/+OCDev/997V9+3Zn27Rp07Rp0yZ9/fXXZ/U9i4uLFRYWpqKiIoWGhjamXLQgY4y+2XtMc1fs1urMAmf71JGJ+tP48yysDADgCs7289u/MS9aVVWl9PR0PfTQQ3Xax4wZozVr1jR4ztdff60xY8bUabvyyis1f/58VVdXKyAgoN45lZWVqqysrNOZlrA4/aC+O1TUIq/tDTYdKNSG7EJJkp+vj345KE6/vbi7ekaHWFsYAMCtNCqMFBQUyG63Kzo6uk57dHS08vLyGjwnLy+vweNrampUUFCg2NjYeuekpqbq0UcfbUxpTbJy1xG9v+lQi38fTxbo76sJKfH69UXdFB/RxupyAABuqFFh5JTTV9Q0xvzoKpsNHd9Q+ymzZs3SzJkznV8XFxcrPr75bxO94rxoxUcEN/vreovQoABdM6STokKCrC4FAODGGhVGIiMj5efnV28UJD8/v97oxykxMTENHu/v768OHTo0eI7NZpPN1vKLZY0fGKfxA+Na/PsAAIAza9TdNIGBgUpOTlZaWlqd9rS0NI0YMaLBc4YPH17v+OXLlyslJaXB+SIAAMC7NPrW3pkzZ+rll1/WK6+8ou3bt+u+++5Tdna2pk2bJqn2EsukSZOcx0+bNk379+/XzJkztX37dr3yyiuaP3++fv/73zdfLwAAgNtq9JyRCRMm6OjRo5ozZ45yc3PVr18/LVu2TAkJCZKk3NzcOmuOJCYmatmyZbrvvvv03HPPKS4uTs8884yuu+665usFAABwW41eZ8QKrDMCAID7OdvPb9brBgAAliKMAAAASxFGAACApQgjAADAUoQRAABgKcIIAACwFGEEAABYijACAAAsRRgBAACWavRy8FY4tUhscXGxxZUAAICzdepz+6cWe3eLMFJSUiJJio+Pt7gSAADQWCUlJQoLCzvj826xN43D4dChQ4cUEhIiHx+fJr1GcXGx4uPjdeDAAY/c34b+uT9P76On90/y/D7SP/fX2n00xqikpERxcXHy9T3zzBC3GBnx9fVV586dm+W1QkNDPfaXTKJ/nsDT++jp/ZM8v4/0z/21Zh9/bETkFCawAgAASxFGAACApbwmjNhsNj3yyCOy2WxWl9Ii6J/78/Q+enr/JM/vI/1zf67aR7eYwAoAADyX14yMAAAA10QYAQAAliKMAAAASxFGAACApbwijMydO1eJiYkKCgpScnKyVq9ebXVJTbZq1SqNHz9ecXFx8vHx0dKlS+s8b4zRn//8Z8XFxSk4OFiXXHKJtm7dak2xTZCamqrzzz9fISEhioqK0tVXX62dO3fWOcad+zhv3jwNGDDAueDQ8OHD9d///tf5vDv3rSGpqany8fHRjBkznG3u3sc///nP8vHxqfOIiYlxPu/u/ZOknJwc3XrrrerQoYPatGmjQYMGKT093fm8u/exa9eu9X6GPj4+uvvuuyW5f/9qamr0xz/+UYmJiQoODla3bt00Z84cORwO5zEu10fj4d5++20TEBBgXnrpJbNt2zYzffp007ZtW7N//36rS2uSZcuWmdmzZ5vFixcbSea9996r8/wTTzxhQkJCzOLFi82WLVvMhAkTTGxsrCkuLram4Ea68sorzauvvmq+++47s3HjRjNu3DjTpUsXU1pa6jzGnfv4/vvvm48++sjs3LnT7Ny50zz88MMmICDAfPfdd8YY9+7b6dauXWu6du1qBgwYYKZPn+5sd/c+PvLII6Zv374mNzfX+cjPz3c+7+79O3bsmElISDBTpkwx3377rcnKyjKffvqp2b17t/MYd+9jfn5+nZ9fWlqakWS++OILY4z79++xxx4zHTp0MB9++KHJysoy77zzjmnXrp15+umnnce4Wh89PowMHTrUTJs2rU5bnz59zEMPPWRRRc3n9DDicDhMTEyMeeKJJ5xtFRUVJiwszDz//PMWVHju8vPzjSSzcuVKY4xn9rF9+/bm5Zdf9qi+lZSUmJ49e5q0tDRz8cUXO8OIJ/TxkUceMQMHDmzwOU/o34MPPmguvPDCMz7vCX083fTp00337t2Nw+HwiP6NGzfOTJ06tU7btddea2699VZjjGv+DD36Mk1VVZXS09M1ZsyYOu1jxozRmjVrLKqq5WRlZSkvL69Of202my6++GK37W9RUZEkKSIiQpJn9dFut+vtt99WWVmZhg8f7lF9u/vuuzVu3Dhdfvnlddo9pY+ZmZmKi4tTYmKibrzxRu3du1eSZ/Tv/fffV0pKiq6//npFRUVp8ODBeumll5zPe0Iff6iqqkpvvPGGpk6dKh8fH4/o34UXXqjPPvtMu3btkiRt2rRJX375pa666ipJrvkzdIuN8pqqoKBAdrtd0dHRddqjo6OVl5dnUVUt51SfGurv/v37rSjpnBhjNHPmTF144YXq16+fJM/o45YtWzR8+HBVVFSoXbt2eu+993Teeec53wTcuW+S9Pbbb2vDhg1at25dvec84ec3bNgwLViwQL169dLhw4f12GOPacSIEdq6datH9G/v3r2aN2+eZs6cqYcfflhr167VvffeK5vNpkmTJnlEH39o6dKlKiws1JQpUyR5xu/ogw8+qKKiIvXp00d+fn6y2+36y1/+optuukmSa/bRo8PIKT4+PnW+NsbUa/MkntLfe+65R5s3b9aXX35Z7zl37mPv3r21ceNGFRYWavHixZo8ebJWrlzpfN6d+3bgwAFNnz5dy5cvV1BQ0BmPc+c+jh071vm/+/fvr+HDh6t79+56/fXXdcEFF0hy7/45HA6lpKTo8ccflyQNHjxYW7du1bx58zRp0iTnce7cxx+aP3++xo4dq7i4uDrt7ty/RYsW6Y033tCbb76pvn37auPGjZoxY4bi4uI0efJk53Gu1EePvkwTGRkpPz+/eqMg+fn59RKhJzg1o98T+vu73/1O77//vr744gt17tzZ2e4JfQwMDFSPHj2UkpKi1NRUDRw4UP/85z89om/p6enKz89XcnKy/P395e/vr5UrV+qZZ56Rv7+/sx/u3MfTtW3bVv3791dmZqZH/AxjY2N13nnn1WlLSkpSdna2JM/4Gzxl//79+vTTT3XHHXc42zyhfw888IAeeugh3Xjjjerfv78mTpyo++67T6mpqZJcs48eHUYCAwOVnJystLS0Ou1paWkaMWKERVW1nMTERMXExNTpb1VVlVauXOk2/TXG6J577tGSJUv0+eefKzExsc7zntDH0xljVFlZ6RF9Gz16tLZs2aKNGzc6HykpKbrlllu0ceNGdevWze37eLrKykpt375dsbGxHvEzHDlyZL3b6Xft2qWEhARJnvU3+OqrryoqKkrjxo1ztnlC/06cOCFf37of735+fs5be12yj5ZMm21Fp27tnT9/vtm2bZuZMWOGadu2rdm3b5/VpTVJSUmJycjIMBkZGUaSeeqpp0xGRobzVuUnnnjChIWFmSVLlpgtW7aYm266ya1uSfvtb39rwsLCzIoVK+rcenfixAnnMe7cx1mzZplVq1aZrKwss3nzZvPwww8bX19fs3z5cmOMe/ftTH54N40x7t/H+++/36xYscLs3bvXfPPNN+bnP/+5CQkJcb6nuHv/1q5da/z9/c1f/vIXk5mZaRYuXGjatGlj3njjDecx7t5HY4yx2+2mS5cu5sEHH6z3nLv3b/LkyaZTp07OW3uXLFliIiMjzR/+8AfnMa7WR48PI8YY89xzz5mEhAQTGBhohgwZ4rxN1B198cUXRlK9x+TJk40xtbdsPfLIIyYmJsbYbDZz0UUXmS1btlhbdCM01DdJ5tVXX3Ue4859nDp1qvN3sWPHjmb06NHOIGKMe/ftTE4PI+7ex1PrMQQEBJi4uDhz7bXXmq1btzqfd/f+GWPMBx98YPr162dsNpvp06ePefHFF+s87wl9/OSTT4wks3PnznrPuXv/iouLzfTp002XLl1MUFCQ6datm5k9e7aprKx0HuNqffQxxhhLhmQAAADk4XNGAACA6yOMAAAASxFGAACApQgjAADAUoQRAABgKcIIAACwFGEEAABYijACAAAsRRgBAACWIowAAABLEUYAAIClCCMAAMBS/x9RBWIX99sy6gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(vct['level'],vct['vct_b'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_vct = '/scratch/e1000/meteoswiss/scratch/swester/experiments/losvec_case_studies/exp3/out_galchen_laydist3/vgrid_DOM01.nc'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: b'/scratch/e1000/meteoswiss/scratch/swester/experiments/losvec_case_studies/exp3/out_galchen_laydist3/vgrid_DOM01.nc'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "File \u001b[0;32m/scratch/e1000/meteoswiss/scratch/swester/miniconda3/envs/vert_coord/lib/python3.10/site-packages/xarray/backends/file_manager.py:209\u001b[0m, in \u001b[0;36mCachingFileManager._acquire_with_cache_info\u001b[0;34m(self, needs_lock)\u001b[0m\n\u001b[1;32m    208\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m--> 209\u001b[0m     file \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_cache[\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_key]\n\u001b[1;32m    210\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mKeyError\u001b[39;00m:\n",
      "File \u001b[0;32m/scratch/e1000/meteoswiss/scratch/swester/miniconda3/envs/vert_coord/lib/python3.10/site-packages/xarray/backends/lru_cache.py:55\u001b[0m, in \u001b[0;36mLRUCache.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m     54\u001b[0m \u001b[39mwith\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_lock:\n\u001b[0;32m---> 55\u001b[0m     value \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_cache[key]\n\u001b[1;32m     56\u001b[0m     \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_cache\u001b[39m.\u001b[39mmove_to_end(key)\n",
      "\u001b[0;31mKeyError\u001b[0m: [<class 'netCDF4._netCDF4.Dataset'>, ('/scratch/e1000/meteoswiss/scratch/swester/experiments/losvec_case_studies/exp3/out_galchen_laydist3/vgrid_DOM01.nc',), 'r', (('clobber', True), ('diskless', False), ('format', 'NETCDF4'), ('persist', False)), '6fa0a022-28de-473e-8d09-040b1f04b6d8']",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[9], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m vct \u001b[39m=\u001b[39m xr\u001b[39m.\u001b[39;49mopen_dataset(f_vct)\u001b[39m.\u001b[39msqueeze()\n",
      "File \u001b[0;32m/scratch/e1000/meteoswiss/scratch/swester/miniconda3/envs/vert_coord/lib/python3.10/site-packages/xarray/backends/api.py:540\u001b[0m, in \u001b[0;36mopen_dataset\u001b[0;34m(filename_or_obj, engine, chunks, cache, decode_cf, mask_and_scale, decode_times, decode_timedelta, use_cftime, concat_characters, decode_coords, drop_variables, inline_array, backend_kwargs, **kwargs)\u001b[0m\n\u001b[1;32m    528\u001b[0m decoders \u001b[39m=\u001b[39m _resolve_decoders_kwargs(\n\u001b[1;32m    529\u001b[0m     decode_cf,\n\u001b[1;32m    530\u001b[0m     open_backend_dataset_parameters\u001b[39m=\u001b[39mbackend\u001b[39m.\u001b[39mopen_dataset_parameters,\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    536\u001b[0m     decode_coords\u001b[39m=\u001b[39mdecode_coords,\n\u001b[1;32m    537\u001b[0m )\n\u001b[1;32m    539\u001b[0m overwrite_encoded_chunks \u001b[39m=\u001b[39m kwargs\u001b[39m.\u001b[39mpop(\u001b[39m\"\u001b[39m\u001b[39moverwrite_encoded_chunks\u001b[39m\u001b[39m\"\u001b[39m, \u001b[39mNone\u001b[39;00m)\n\u001b[0;32m--> 540\u001b[0m backend_ds \u001b[39m=\u001b[39m backend\u001b[39m.\u001b[39;49mopen_dataset(\n\u001b[1;32m    541\u001b[0m     filename_or_obj,\n\u001b[1;32m    542\u001b[0m     drop_variables\u001b[39m=\u001b[39;49mdrop_variables,\n\u001b[1;32m    543\u001b[0m     \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mdecoders,\n\u001b[1;32m    544\u001b[0m     \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs,\n\u001b[1;32m    545\u001b[0m )\n\u001b[1;32m    546\u001b[0m ds \u001b[39m=\u001b[39m _dataset_from_backend_dataset(\n\u001b[1;32m    547\u001b[0m     backend_ds,\n\u001b[1;32m    548\u001b[0m     filename_or_obj,\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    556\u001b[0m     \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs,\n\u001b[1;32m    557\u001b[0m )\n\u001b[1;32m    558\u001b[0m \u001b[39mreturn\u001b[39;00m ds\n",
      "File \u001b[0;32m/scratch/e1000/meteoswiss/scratch/swester/miniconda3/envs/vert_coord/lib/python3.10/site-packages/xarray/backends/netCDF4_.py:572\u001b[0m, in \u001b[0;36mNetCDF4BackendEntrypoint.open_dataset\u001b[0;34m(self, filename_or_obj, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta, group, mode, format, clobber, diskless, persist, lock, autoclose)\u001b[0m\n\u001b[1;32m    551\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mopen_dataset\u001b[39m(\n\u001b[1;32m    552\u001b[0m     \u001b[39mself\u001b[39m,\n\u001b[1;32m    553\u001b[0m     filename_or_obj,\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    568\u001b[0m     autoclose\u001b[39m=\u001b[39m\u001b[39mFalse\u001b[39;00m,\n\u001b[1;32m    569\u001b[0m ):\n\u001b[1;32m    571\u001b[0m     filename_or_obj \u001b[39m=\u001b[39m _normalize_path(filename_or_obj)\n\u001b[0;32m--> 572\u001b[0m     store \u001b[39m=\u001b[39m NetCDF4DataStore\u001b[39m.\u001b[39;49mopen(\n\u001b[1;32m    573\u001b[0m         filename_or_obj,\n\u001b[1;32m    574\u001b[0m         mode\u001b[39m=\u001b[39;49mmode,\n\u001b[1;32m    575\u001b[0m         \u001b[39mformat\u001b[39;49m\u001b[39m=\u001b[39;49m\u001b[39mformat\u001b[39;49m,\n\u001b[1;32m    576\u001b[0m         group\u001b[39m=\u001b[39;49mgroup,\n\u001b[1;32m    577\u001b[0m         clobber\u001b[39m=\u001b[39;49mclobber,\n\u001b[1;32m    578\u001b[0m         diskless\u001b[39m=\u001b[39;49mdiskless,\n\u001b[1;32m    579\u001b[0m         persist\u001b[39m=\u001b[39;49mpersist,\n\u001b[1;32m    580\u001b[0m         lock\u001b[39m=\u001b[39;49mlock,\n\u001b[1;32m    581\u001b[0m         autoclose\u001b[39m=\u001b[39;49mautoclose,\n\u001b[1;32m    582\u001b[0m     )\n\u001b[1;32m    584\u001b[0m     store_entrypoint \u001b[39m=\u001b[39m StoreBackendEntrypoint()\n\u001b[1;32m    585\u001b[0m     \u001b[39mwith\u001b[39;00m close_on_error(store):\n",
      "File \u001b[0;32m/scratch/e1000/meteoswiss/scratch/swester/miniconda3/envs/vert_coord/lib/python3.10/site-packages/xarray/backends/netCDF4_.py:376\u001b[0m, in \u001b[0;36mNetCDF4DataStore.open\u001b[0;34m(cls, filename, mode, format, group, clobber, diskless, persist, lock, lock_maker, autoclose)\u001b[0m\n\u001b[1;32m    370\u001b[0m kwargs \u001b[39m=\u001b[39m \u001b[39mdict\u001b[39m(\n\u001b[1;32m    371\u001b[0m     clobber\u001b[39m=\u001b[39mclobber, diskless\u001b[39m=\u001b[39mdiskless, persist\u001b[39m=\u001b[39mpersist, \u001b[39mformat\u001b[39m\u001b[39m=\u001b[39m\u001b[39mformat\u001b[39m\n\u001b[1;32m    372\u001b[0m )\n\u001b[1;32m    373\u001b[0m manager \u001b[39m=\u001b[39m CachingFileManager(\n\u001b[1;32m    374\u001b[0m     netCDF4\u001b[39m.\u001b[39mDataset, filename, mode\u001b[39m=\u001b[39mmode, kwargs\u001b[39m=\u001b[39mkwargs\n\u001b[1;32m    375\u001b[0m )\n\u001b[0;32m--> 376\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mcls\u001b[39;49m(manager, group\u001b[39m=\u001b[39;49mgroup, mode\u001b[39m=\u001b[39;49mmode, lock\u001b[39m=\u001b[39;49mlock, autoclose\u001b[39m=\u001b[39;49mautoclose)\n",
      "File \u001b[0;32m/scratch/e1000/meteoswiss/scratch/swester/miniconda3/envs/vert_coord/lib/python3.10/site-packages/xarray/backends/netCDF4_.py:323\u001b[0m, in \u001b[0;36mNetCDF4DataStore.__init__\u001b[0;34m(self, manager, group, mode, lock, autoclose)\u001b[0m\n\u001b[1;32m    321\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_group \u001b[39m=\u001b[39m group\n\u001b[1;32m    322\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_mode \u001b[39m=\u001b[39m mode\n\u001b[0;32m--> 323\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mformat \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mds\u001b[39m.\u001b[39mdata_model\n\u001b[1;32m    324\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_filename \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mds\u001b[39m.\u001b[39mfilepath()\n\u001b[1;32m    325\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mis_remote \u001b[39m=\u001b[39m is_remote_uri(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_filename)\n",
      "File \u001b[0;32m/scratch/e1000/meteoswiss/scratch/swester/miniconda3/envs/vert_coord/lib/python3.10/site-packages/xarray/backends/netCDF4_.py:385\u001b[0m, in \u001b[0;36mNetCDF4DataStore.ds\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    383\u001b[0m \u001b[39m@property\u001b[39m\n\u001b[1;32m    384\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mds\u001b[39m(\u001b[39mself\u001b[39m):\n\u001b[0;32m--> 385\u001b[0m     \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_acquire()\n",
      "File \u001b[0;32m/scratch/e1000/meteoswiss/scratch/swester/miniconda3/envs/vert_coord/lib/python3.10/site-packages/xarray/backends/netCDF4_.py:379\u001b[0m, in \u001b[0;36mNetCDF4DataStore._acquire\u001b[0;34m(self, needs_lock)\u001b[0m\n\u001b[1;32m    378\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39m_acquire\u001b[39m(\u001b[39mself\u001b[39m, needs_lock\u001b[39m=\u001b[39m\u001b[39mTrue\u001b[39;00m):\n\u001b[0;32m--> 379\u001b[0m     \u001b[39mwith\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_manager\u001b[39m.\u001b[39macquire_context(needs_lock) \u001b[39mas\u001b[39;00m root:\n\u001b[1;32m    380\u001b[0m         ds \u001b[39m=\u001b[39m _nc4_require_group(root, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_group, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_mode)\n\u001b[1;32m    381\u001b[0m     \u001b[39mreturn\u001b[39;00m ds\n",
      "File \u001b[0;32m/scratch/e1000/meteoswiss/scratch/swester/miniconda3/envs/vert_coord/lib/python3.10/contextlib.py:135\u001b[0m, in \u001b[0;36m_GeneratorContextManager.__enter__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    133\u001b[0m \u001b[39mdel\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39margs, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mkwds, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mfunc\n\u001b[1;32m    134\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m--> 135\u001b[0m     \u001b[39mreturn\u001b[39;00m \u001b[39mnext\u001b[39;49m(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mgen)\n\u001b[1;32m    136\u001b[0m \u001b[39mexcept\u001b[39;00m \u001b[39mStopIteration\u001b[39;00m:\n\u001b[1;32m    137\u001b[0m     \u001b[39mraise\u001b[39;00m \u001b[39mRuntimeError\u001b[39;00m(\u001b[39m\"\u001b[39m\u001b[39mgenerator didn\u001b[39m\u001b[39m'\u001b[39m\u001b[39mt yield\u001b[39m\u001b[39m\"\u001b[39m) \u001b[39mfrom\u001b[39;00m \u001b[39mNone\u001b[39m\n",
      "File \u001b[0;32m/scratch/e1000/meteoswiss/scratch/swester/miniconda3/envs/vert_coord/lib/python3.10/site-packages/xarray/backends/file_manager.py:197\u001b[0m, in \u001b[0;36mCachingFileManager.acquire_context\u001b[0;34m(self, needs_lock)\u001b[0m\n\u001b[1;32m    194\u001b[0m \u001b[39m@contextlib\u001b[39m\u001b[39m.\u001b[39mcontextmanager\n\u001b[1;32m    195\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39macquire_context\u001b[39m(\u001b[39mself\u001b[39m, needs_lock\u001b[39m=\u001b[39m\u001b[39mTrue\u001b[39;00m):\n\u001b[1;32m    196\u001b[0m     \u001b[39m\"\"\"Context manager for acquiring a file.\"\"\"\u001b[39;00m\n\u001b[0;32m--> 197\u001b[0m     file, cached \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_acquire_with_cache_info(needs_lock)\n\u001b[1;32m    198\u001b[0m     \u001b[39mtry\u001b[39;00m:\n\u001b[1;32m    199\u001b[0m         \u001b[39myield\u001b[39;00m file\n",
      "File \u001b[0;32m/scratch/e1000/meteoswiss/scratch/swester/miniconda3/envs/vert_coord/lib/python3.10/site-packages/xarray/backends/file_manager.py:215\u001b[0m, in \u001b[0;36mCachingFileManager._acquire_with_cache_info\u001b[0;34m(self, needs_lock)\u001b[0m\n\u001b[1;32m    213\u001b[0m     kwargs \u001b[39m=\u001b[39m kwargs\u001b[39m.\u001b[39mcopy()\n\u001b[1;32m    214\u001b[0m     kwargs[\u001b[39m\"\u001b[39m\u001b[39mmode\u001b[39m\u001b[39m\"\u001b[39m] \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_mode\n\u001b[0;32m--> 215\u001b[0m file \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_opener(\u001b[39m*\u001b[39;49m\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_args, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[1;32m    216\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_mode \u001b[39m==\u001b[39m \u001b[39m\"\u001b[39m\u001b[39mw\u001b[39m\u001b[39m\"\u001b[39m:\n\u001b[1;32m    217\u001b[0m     \u001b[39m# ensure file doesn't get overridden when opened again\u001b[39;00m\n\u001b[1;32m    218\u001b[0m     \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_mode \u001b[39m=\u001b[39m \u001b[39m\"\u001b[39m\u001b[39ma\u001b[39m\u001b[39m\"\u001b[39m\n",
      "File \u001b[0;32msrc/netCDF4/_netCDF4.pyx:2463\u001b[0m, in \u001b[0;36mnetCDF4._netCDF4.Dataset.__init__\u001b[0;34m()\u001b[0m\n",
      "File \u001b[0;32msrc/netCDF4/_netCDF4.pyx:2026\u001b[0m, in \u001b[0;36mnetCDF4._netCDF4._ensure_nc_success\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: b'/scratch/e1000/meteoswiss/scratch/swester/experiments/losvec_case_studies/exp3/out_galchen_laydist3/vgrid_DOM01.nc'"
     ]
    }
   ],
   "source": [
    "vct = xr.open_dataset(f_vct).squeeze()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 22000.00 0.000\n",
      "2 21548.48 0.000\n",
      "3 21101.44 0.000\n",
      "4 20658.86 0.000\n",
      "5 20220.75 0.000\n",
      "6 19787.11 0.000\n",
      "7 19357.94 0.000\n",
      "8 18933.23 0.000\n",
      "9 18513.00 0.000\n",
      "10 18097.23 0.000\n",
      "11 17685.94 0.000\n",
      "12 17279.11 0.000\n",
      "13 16876.75 0.000\n",
      "14 16478.86 0.000\n",
      "15 16085.44 0.000\n",
      "16 15696.48 0.019\n",
      "17 15312.00 0.043\n",
      "18 14931.98 0.067\n",
      "19 14556.44 0.090\n",
      "20 14185.36 0.113\n",
      "21 13818.75 0.136\n",
      "22 13456.61 0.159\n",
      "23 13098.94 0.181\n",
      "24 12745.73 0.203\n",
      "25 12397.00 0.225\n",
      "26 12052.73 0.247\n",
      "27 11712.94 0.268\n",
      "28 11377.61 0.289\n",
      "29 11046.75 0.310\n",
      "30 10720.36 0.330\n",
      "31 10398.44 0.350\n",
      "32 10080.98 0.370\n",
      "33 9768.00 0.390\n",
      "34 9459.48 0.409\n",
      "35 9155.44 0.428\n",
      "36 8855.86 0.447\n",
      "37 8560.75 0.465\n",
      "38 8270.11 0.483\n",
      "39 7983.94 0.501\n",
      "40 7702.23 0.519\n",
      "41 7425.00 0.536\n",
      "42 7152.23 0.553\n",
      "43 6883.94 0.570\n",
      "44 6620.11 0.586\n",
      "45 6360.75 0.602\n",
      "46 6105.86 0.618\n",
      "47 5855.44 0.634\n",
      "48 5609.48 0.649\n",
      "49 5368.00 0.664\n",
      "50 5130.98 0.679\n",
      "51 4898.44 0.694\n",
      "52 4670.36 0.708\n",
      "53 4446.75 0.722\n",
      "54 4227.61 0.736\n",
      "55 4012.94 0.749\n",
      "56 3802.73 0.762\n",
      "57 3597.00 0.775\n",
      "58 3395.73 0.788\n",
      "59 3198.94 0.800\n",
      "60 3006.61 0.812\n",
      "61 2818.75 0.824\n",
      "62 2635.36 0.835\n",
      "63 2456.44 0.846\n",
      "64 2281.98 0.857\n",
      "65 2112.00 0.868\n",
      "66 1946.48 0.878\n",
      "67 1785.44 0.888\n",
      "68 1628.86 0.898\n",
      "69 1476.75 0.908\n",
      "70 1329.11 0.917\n",
      "71 1185.94 0.926\n",
      "72 1047.23 0.935\n",
      "73 913.00 0.943\n",
      "74 783.23 0.951\n",
      "75 657.94 0.959\n",
      "76 537.11 0.966\n",
      "77 420.75 0.974\n",
      "78 308.86 0.981\n",
      "79 201.44 0.987\n",
      "80 98.48 0.994\n",
      "81 0.00 1.000\n"
     ]
    }
   ],
   "source": [
    "for i, var in enumerate(vct['vct_b'].values): print(f'{i+1} {vct[\"vct_a\"].values[i]:.2f} {var:.3f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x149717f0a860>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9uklEQVR4nO3deVhVdeLH8c9ldwNTFEQRwX1fQE3T9nDMnHGs1JpcsvqNZbm1qtOiU1EzU9NiauaWo6WZVjpZI21qqakIamguoYIKIi6AKNu9398fFDO4FCBw7r28X89znyfPPYf7+T4Q98O55/s9NmOMEQAAgEU8rA4AAACqN8oIAACwFGUEAABYijICAAAsRRkBAACWoowAAABLUUYAAIClKCMAAMBSXlYHKA2Hw6Fjx46pTp06stlsVscBAAClYIxRdna2QkJC5OFx+fMfLlFGjh07ptDQUKtjAACAckhJSVGTJk0u+7xLlJE6depIKhqMv7+/xWkAAEBpZGVlKTQ0tPh9/HJcooz88tGMv78/ZQQAABfzW5dYcAErAACwFGUEAABYijICAAAsRRkBAACWoowAAABLUUYAAIClKCMAAMBSlBEAAGApyggAALBUmcvI+vXrNXDgQIWEhMhms+njjz/+zWPWrVunyMhI+fn5KSIiQrNnzy5PVgAA4IbKXEZycnLUuXNnzZgxo1T7Hzx4ULfeeqv69u2r+Ph4TZkyRePGjdOKFSvKHBYAALifMt+bpn///urfv3+p9589e7aaNm2q1157TZLUtm1bbdu2Tf/4xz90++23l/XlAQCAm6n0a0Y2bdqk6OjoEtv69eunbdu2qaCg4JLH5OXlKSsrq8QDAABUrNwCu1bvOKb7Fm7VqZx8y3JU+l1709LSFBQUVGJbUFCQCgsLlZGRoUaNGl10TExMjKZNm1bZ0QAAqHYcDqNth09r5fYj+nRnqrLzCiVJ/955TCN6NbMkU6WXEeniWwcbYy65/ReTJ0/WpEmTiv+dlZWl0NDQygsIAICbSz55Th9uP6KP4o8o5dT54u2N69bQ4G6NdV2rBpZlq/QyEhwcrLS0tBLb0tPT5eXlpfr161/yGF9fX/n6+lZ2NAAA3Nr5fLs+T0zVB1uPaFPSyeLttX29dGvHYA3u1kQ9mtWTh8elTw5UlUovI7169dLq1atLbFu7dq2ioqLk7e1d2S8PAEC1YozRrqOZWrY1RasSjhV/DGOzSX1aBOqOyCaKbhesGj6eFif9rzKXkbNnz+rAgQPF/z548KASEhJUr149NW3aVJMnT9bRo0e1aNEiSdKYMWM0Y8YMTZo0SQ888IA2bdqkefPm6f3336+4UQAAUM1l5xbok4Rjen9LshKP/XfiR5OramhIVKhuj2yixnVrWJjw8spcRrZt26Ybbrih+N+/XNsxcuRILVy4UKmpqUpOTi5+Pjw8XGvWrNHEiRP11ltvKSQkRG+88QbTegEAqAA7j5zRe98na9WOYzqXb5ck+Xh66HcdgjWse6iujqhv+ccwv8Vmfrma1IllZWUpICBAmZmZ8vf3tzoOAACWyi2wa9WOY1q8+bB2Hsks3h7RoJbu7tFUt3droqtq+ViYsEhp37+rZDYNAAC4codP5mjx5sP6YNsRZZ4vWqvLx9ND/TsG6+4eTdUjvN5lZ6o6M8oIAABOzOEwWr//hBZuPKR1+07ol88zGtetoXuuDtPQ7qGq5wRnQa4EZQQAACeUk1eolduPaMHGQ0o6kVO8/bpWDTSiV5iub91Qnk5+LUhpUUYAAHAiKafOadGmQ1q6NUXZuUXTcmv7emlIVKhG9ApTs8BaFieseJQRAACcQELKGb2zIUmf7UqV4+ePYprVr6lRvZvpjqhQ1fZ137ds9x0ZAABOzuEw+mLPcc3dcFBbDp0q3t6nRaBG92mm61s1dPppuRWBMgIAQBXLLbBr5fajemdDkg5mFF0P4u1p08DOIbq/T4TahVSvZSwoIwAAVJGs3AIt2Zys+d8d1InsPEmSv5+X/nR1mEb2aqbgAD+LE1qDMgIAQCVLz8rV/O8Oacnmw8X3igkJ8NN9fSM0rHuoarnx9SClUb1HDwBAJTpy+pxmr/tJH2w7ovxChySpZcPaGnNdc/2+S4i8PT0sTugcKCMAAFSwpBNnNeubn/RR/FEV/jw1plvTunro+ha6sU31uCi1LCgjAABUkH3HszXjqwP6985jxdNz+7QI1MM3tlBPF12qvSpQRgAAuEJ707L1xlf7tWZXavFy7Te1aaixN7ZQt6ZXWRvOBVBGAAAop33Hs/X6lyVLyO/aB+vhG1uoQ+MAa8O5EMoIAABldCA9W//8omQJ6d8hWONuaqm2jarXGiEVgTICAEApHT6Zo9e/2K+PE44WXxNCCblylBEAAH7DsTPn9eZX+/XBtiOy/9xC+rUP0oSbW1FCKgBlBACAy8g4m6cZXx3Qe98nK99etE7I9a0b6NFbWqtjE64JqSiUEQAALpCdW6C5Gw5q7oYk5eTbJUlXR9TTY9GtFdWsnsXp3A9lBACAn+UV2rVkc7JmfH1Ap3LyJUmdmgToiX5tdE2L+qwTUkkoIwCAas/hMPo44aheWbtPR8+clyRFBNbSY/1aq3+HYEpIJaOMAACqtW/3Z+jFNXu0OzVLkhTk76sJN7fSnZFN5MW9Y6oEZQQAUC39mJalmDU/at2+E5KkOr5eevCG5rq3d7hq+HhanK56oYwAAKqV41m5emXtXn0Yd0QOI3l52HTP1WEad1NL1avlY3W8aokyAgCoFs7n2/XOhiTNXveTzv08Q+bWjsF6ol8bNQusZXG66o0yAgBwa8YYrdpxTC9/9qOOZeZKkro2rau/DGinyDBuYucMKCMAALe1Pfm0pq/erYSUM5KkxnVr6Mn+bTSwUyNmyDgRyggAwO0cz8rVy5/9qJXxRyVJtXw89dANLXRfn3D5eXNxqrOhjAAA3EZeoV3zvj2oGV8d0Ll8u2w26Y5uTfR4v9Zq6O9ndTxcBmUEAODyjDH6ck+6/vrpbh0+eU6S1K1pXT33+/bq1KSuteHwmygjAACXdigjR8+tTtQ3e4vWC2lYx1dP9W+jQV0ay8OD60JcAWUEAOCSzufbNfObA3p7XZLy7Q55e9p0f98Ijb2hhWr78vbmSvhuAQBcijFGsbuPa9rq3cX3kbm2VQM9N7CdIhrUtjgdyoMyAgBwGYdP5ui5VYn6+uePZBrXraGnb2unfu2DmKrrwigjAACnl1do1zvrk/TmVweUV+iQj6eHHrg2XGNvaKGaPryVuTq+gwAAp7bxpwz95eMflHQiR5LUp0Wgpv+hPR/JuBHKCADAKWWczdOLn+4pXrgssLavnr6trX7fOYSPZNwMZQQA4FSMMVq+7YheWLNHmecLZLNJ9/QM02P9WiughrfV8VAJKCMAAKeRdOKspny0S5uTTkmS2jXy1wt/7KCuTbmhnTujjAAALJdf6NCc9T/pja8OKL/QIT9vD026pZVGXxMuL08Pq+OhklFGAACWik8+radW7NLe49mSpL4tA/XiHzsqtF5Ni5OhqlBGAACWOJdfqFfW7tP87w7KGKleLR89c1s7/aELF6hWN5QRAECV++5Ahp5auVMpp4pWUB3crbH+MqCd6tXysTgZrEAZAQBUmczzBXrx0z1ati1FkhQS4KcXBnfUDa0bWpwMVqKMAACqxJd7jmvyyl1Kz86TJI3oFaYnfteGm9qBMgIAqFyZ5wo07d+JWrm9aPGy8MBaevn2TuoRXs/iZHAWlBEAQKX537MhNpt0f59wPRrdWn7enlZHgxOhjAAAKtyFZ0MiAmvp73d2UmQYZ0NwMcoIAKBCrdt3Qk98uEPHs4rOhjzQN0KTbmnF2RBcFmUEAFAhcvIK9cKaPXrv+2RJv5wN6azIMJZyx6+jjAAArtiWg6f02PIdSj51TpI0qnczPfm7Nqrhw9kQ/DbKCACg3HIL7Ho1dp/e2ZAkY6TGdWvo73d0Uu8WgVZHgwuhjAAAymVPapYmLE0ovqfMnZFN9PTAdvL387Y4GVwNZQQAUCYOh9Hcb5P0j//sU77docDaPooZ3Em3tAuyOhpcFGUEAFBqR8+c16MfJGhz0ilJ0s1tG+ql2zspsLavxcngyigjAIBS+SThqP7y8Q/Kzi1UTR9PPX1bOw3rHsoddnHFKCMAgF+VnVugZz5J1EfxRQuYdQmtq9eGdlGzwFoWJ4O78CjPQTNnzlR4eLj8/PwUGRmpDRs2/Or+S5YsUefOnVWzZk01atRI9957r06ePFmuwACAqhN3+LRufWODPoo/Kk8Pmybc3FIfjulFEUGFKnMZWbZsmSZMmKCpU6cqPj5effv2Vf/+/ZWcnHzJ/b/99luNGDFC9913nxITE7V8+XJt3bpV999//xWHBwBUDrvD6I0v92vI25uUcuq8mlxVQx/8+WpNuLmVvDzL9XcscFk2Y4wpywE9e/ZUt27dNGvWrOJtbdu21aBBgxQTE3PR/v/4xz80a9Ys/fTTT8Xb3nzzTf3tb39TSkpKqV4zKytLAQEByszMlL+/f1niAgDK6OiZ85q4NEFbDhVdpPqHLiH666AOTNlFmZX2/btM9TY/P19xcXGKjo4usT06OlobN2685DG9e/fWkSNHtGbNGhljdPz4cX344YcaMGDAZV8nLy9PWVlZJR4AgMr3+Q+p6v/aem05dEq1fb30z6Gd9fqwrhQRVKoylZGMjAzZ7XYFBZWcSx4UFKS0tLRLHtO7d28tWbJEQ4cOlY+Pj4KDg1W3bl29+eabl32dmJgYBQQEFD9CQ0PLEhMAUEa5BXb95eNdGrN4u7JyC9U5tK4+HddHf+zaxOpoqAbK9cHfhdO4jDGXndq1e/dujRs3Ts8884zi4uL0+eef6+DBgxozZsxlv/7kyZOVmZlZ/CjtxzkAgLI7kJ6tQW99p8Wbi679+/N1EfpwTC+F1eciVVSNMk3tDQwMlKen50VnQdLT0y86W/KLmJgYXXPNNXr88cclSZ06dVKtWrXUt29fPf/882rUqNFFx/j6+srXlwV0AKAyGWO0fNsRPbsqUecL7Aqs7aNXhnTRda0aWB0N1UyZzoz4+PgoMjJSsbGxJbbHxsaqd+/elzzm3Llz8vAo+TKenkV3cSzjtbMAgAqSk1eoicsS9MSKnTpfYFffloFaM74vRQSWKPOiZ5MmTdLw4cMVFRWlXr16ac6cOUpOTi7+2GXy5Mk6evSoFi1aJEkaOHCgHnjgAc2aNUv9+vVTamqqJkyYoB49eigkJKRiRwMA+E17UrM09r3tSjqRI08Pmx6NbqUx1zaXhwcrqcIaZS4jQ4cO1cmTJzV9+nSlpqaqQ4cOWrNmjcLCwiRJqampJdYcGTVqlLKzszVjxgw9+uijqlu3rm688Ua9/PLLFTcKAMBvMsZo6dYUPbcqUXmFDjUK8NObd3VVVLN6VkdDNVfmdUaswDojAHBlzuYVaupHu/RJwjFJ0g2tG+iVIV1Ur5aPxcngzkr7/s29aQDAze1Ny9aDi+OUlFH0sczj/Vrr//pG8LEMnAZlBADc2MrtRzTlo13KLXAo2N9PM+7mYxk4H8oIALih3AK7pq3erfe3FF3D17dloF4f1pWPZeCUKCMA4GaST57Tg0vilHgsSzabNO7Glhp3U0t58rEMnBRlBADcyJd7jmvisgRl5Rbqqpreem1YV9YOgdOjjACAG7A7jF77Yp/e/OqAJKlr07p66+5uCqlbw+JkwG+jjACAizudk6/xyxK0ft8JSdLIXmGaOqCdfLzKdfsxoMpRRgDAhf1wNFNjFsfpyOnz8vP2UMzgjtxpFy6HMgIALmr5thT95eMflFfoUNN6NTX7nki1C2FhSLgeyggAuJj8Qoem/ztRizcXTdu9qU1DvTqkiwJqelucDCgfyggAuJD07FyNXbJdWw+dls0mTbiplR65sQWrqcKlUUYAwEXEJ5/Wg4u3Ky0rV3V8vfTasC66qW2Q1bGAK0YZAQAX8MHWoutD8u0OtWhYW3OGRyqiQW2rYwEVgjICAE6swO7Q9NW79a/NhyVJ/doH6ZUhXVTbl1/fcB/8NAOAkzp5Nk8PLdmu7w+eks0mTbq5lcbewPUhcD+UEQBwQruPZemBRdt09Mx51fb10mtDu+jmdlwfAvdEGQEAJ/PZrlRN+mCHzhfY1ax+Tb0zIkotg+pYHQuoNJQRAHASDofRa1/u1xtf7pck9W0ZqBl3dWP9ELg9yggAOIFz+YWatGyHPk9MkyTd1ydck/u3kZcn95eB+6OMAIDFjp05r/vf3abdqVny8fTQC3/soDujQq2OBVQZyggAWCg++bQeWBSnjLN5ql/LR3NGRCoyrJ7VsYAqRRkBAIt8knBUj3+4U/mFDrUJrqO5I6PU5KqaVscCqhxlBACqmMNh9GrsPs34+oAk6ea2QXptGAuZofriJx8AqtD5fLseXZ6gNbuKLlT983UReqJfG3mykBmqMcoIAFSR9OxcPfDuNu04kilvT5tiBnfSHZFNrI4FWI4yAgBV4Me0LN23sGhF1bo1vfX2PZHqGVHf6liAU6CMAEAl+3pvuh55L15n8woVEVhL80Z1V3hgLatjAU6DMgIAlWjRpkN6blWiHEa6OqKeZt8Tqbo1fayOBTgVyggAVAK7w+iFT/do/ncHJUl3RjbRC3/sKB8vVlQFLkQZAYAKdj7frvFL47V293FJ0uP9Wuuh65vLZmPGDHAplBEAqEAnsvN0/6Jt2pFyRj6eHvrHkM76fecQq2MBTo0yAgAV5EB6tkYt2Kojp4tmzMwZHqUe4SztDvwWyggAVIDNSSf1f4u2KSu3UGH1a2rBqO6KaFDb6liAS6CMAMAVWr3jmB79YIfy7Q51a1pX74yIUv3avlbHAlwGZQQAyskYo7kbDuqFNXskSb9rH6zXhnWRn7enxckA10IZAYBysDuM/vrv3Vq48ZAkaVTvZnr6tnbcYwYoB8oIAJRRboFdE5cl6LMfim52N/XWtrq/bzhTd4FyoowAQBmcOZev+9/dpm2HTzN1F6gglBEAKKWjZ85r5PwtOpB+VnX8vPTOiChdzc3ugCtGGQGAUtiblq2R87coLStXwf5+end0D7UOrmN1LMAtUEYA4Dd8n3RSD/y8hkiLhrW1aHQPhdStYXUswG1QRgDgV3z+Q6rGLU1QfqFDUWFXae7IKO66C1QwyggAXMbizYf19Cc/yBjplnZBevOurqwhAlQCyggAXMAYoze+PKB/frFPknRXj6b66x/ay8vTw+JkgHuijADA/3A4jKatTtS7mw5Lksbd1FITb27JGiJAJaKMAMDP8gsdenT5Dq3ecUyS9NzAdhp1TbjFqQD3RxkBAEk5eYUaszhOG/ZnyMvDpleGdNYfujS2OhZQLVBGAFR7p3Pyde/CrUpIOaMa3p6aPTxS17VqYHUsoNqgjACo1o5n5Wr4vO+17/hZ1a3prQWjuqtr06usjgVUK5QRANXW4ZM5umfe90o5dV5B/r5afF9PtQxiVVWgqlFGAFRLe9Oydc+873UiO09h9Wtq8X09FVqvptWxgGqJMgKg2olPPq1RC7Yq83yB2gTX0aLRPdTQ38/qWEC1RRkBUK18dyBDDyzapnP5dnVrWlcLRvVQQE1vq2MB1RplBEC1sTYxTQ+/F698u0N9Wwbq7eGRqunDr0HAavxfCKBa+CThqCZ9sEN2h1G/9kF6466u8vXiPjOAM6CMAHB7729J1pSPdskYaXDXxvrbHZ24zwzgRMr1f+PMmTMVHh4uPz8/RUZGasOGDb+6f15enqZOnaqwsDD5+vqqefPmmj9/frkCA0BZzN2QpMkri4rIPVc31T/u7EwRAZxMmc+MLFu2TBMmTNDMmTN1zTXX6O2331b//v21e/duNW3a9JLHDBkyRMePH9e8efPUokULpaenq7Cw8IrDA8DlGGP05lcH9Gps0Z13/3xthJ7q34Yb3gFOyGaMMWU5oGfPnurWrZtmzZpVvK1t27YaNGiQYmJiLtr/888/17Bhw5SUlKR69eqVK2RWVpYCAgKUmZkpf3//cn0NANWHMUYvff6j3l6XJEl69JZWevjGFhQRoIqV9v27TOcq8/PzFRcXp+jo6BLbo6OjtXHjxkses2rVKkVFRelvf/ubGjdurFatWumxxx7T+fPnL/s6eXl5ysrKKvEAgNJwOIymrd5dXESevq2dHrmpJUUEcGJl+pgmIyNDdrtdQUFBJbYHBQUpLS3tksckJSXp22+/lZ+fnz766CNlZGTooYce0qlTpy573UhMTIymTZtWlmgAIIfDaOrHu/T+lhRJ0gt/7KA/9QyzOBWA31Kuq7gu/AvDGHPZvzocDodsNpuWLFmiHj166NZbb9Wrr76qhQsXXvbsyOTJk5WZmVn8SElJKU9MANVIod2hx5bv0PtbUuRhk/5xZ2eKCOAiynRmJDAwUJ6enhedBUlPT7/obMkvGjVqpMaNGysgIKB4W9u2bWWM0ZEjR9SyZcuLjvH19ZWvr29ZogGoxgrsDk1YmqBPd6XK08Om14Z20cDOIVbHAlBKZToz4uPjo8jISMXGxpbYHhsbq969e1/ymGuuuUbHjh3T2bNni7ft27dPHh4eatKkSTkiA8B/5RXa9dCS7fp0V6q8PW166+5uFBHAxZT5Y5pJkyZp7ty5mj9/vvbs2aOJEycqOTlZY8aMkVT0EcuIESOK97/77rtVv3593Xvvvdq9e7fWr1+vxx9/XKNHj1aNGjUqbiQAqp3cArvG/CtOsbuPy8fLQ3OGR+l3HYKtjgWgjMq8zsjQoUN18uRJTZ8+XampqerQoYPWrFmjsLCiz2ZTU1OVnJxcvH/t2rUVGxurRx55RFFRUapfv76GDBmi559/vuJGAaDayS2w6//+Faf1+07Iz9tDc0d0V5+WgVbHAlAOZV5nxAqsMwLgf53Pt+uBRdv07YEM1fD21LxRUerdnCICOJvSvn9zbxoALuVcfqFGL9yqzUmnVNPHUwtGdVfPiPpWxwJwBSgjAFzG2bxCjV6wVVsOnVJtXy8tvLe7opqVb2VnAM6DMgLAJZzNK9So+Vu07fBp1fH10rv39VC3pldZHQtABaCMAHB6Z/MKNXL+FsUdPi1/Py/9676e6hxa1+pYACoIZQSAU7uwiCy+v6c6NalrdSwAFahcy8EDQFWgiADVA2UEgFOiiADVB2UEgNOhiADVC2UEgFPJ+XnWDEUEqD4oIwCcxrn8Qt27cGvR9F2KCFBtUEYAOIXz+Xbdt3Cbthw8pTq+RdN3KSJA9UAZAWC53IKie81sSjqp2j8vaNaFdUSAaoMyAsBSv9x999sDGarp46mF93ZnZVWgmqGMALBMXqFdDy3ZrvX7TqiGd9FN77jXDFD9UEYAWKLA7tAj78Xrqx/T5eftofncfReotigjAKqc3WE0cVmC1u4+Lh8vD80d0V29mlNEgOqKMgKgSjkcRo9/uEP/3pkqb0+bZt/TTX1aBlodC4CFKCMAqowxRlM/3qWV24/K08OmN+/qphvbBFkdC4DFKCMAqoQxRtNW79b7W1LkYZNeG9pFv+sQbHUsAE6AMgKg0hlj9NJnP2rhxkOy2aS/39FZAzuHWB0LgJOgjACodG98eUBvr0+SJL0wqKNuj2xicSIAzoQyAqBSvbM+Sf/8Yp8k6Znb2ununk0tTgTA2VBGAFSaxZsP64U1eyRJj0W30ug+4RYnAuCMKCMAKsXK7Uf09Cc/SJIevL65xt7QwuJEAJwVZQRAhftsV6oeW75Dxkgje4XpiX6tZbPZrI4FwElRRgBUqG/2pmvc0ng5jHRnZBM9O7A9RQTAr6KMAKgwWw+d0pjFcSqwGw3o1Egv3d5JHh4UEQC/jjICoEL8cDRToxdsVW6BQze0bqB/DukiT4oIgFKgjAC4YgfSz2rE/C3KzitUj/B6mvmnSPl48esFQOnw2wLAFTly+pyGz/tep3Ly1bFxgOaNjFINH0+rYwFwIZQRAOWWnp2re+Z+r9TMXLVoWFvvju6hOn7eVscC4GIoIwDKJfNcgUbM26JDJ8+pyVU1tPi+nqpXy8fqWABcEGUEQJmdz7dr9Ltb9WNathrU8dWS+3sqOMDP6lgAXBRlBECZFNgdenBJnOIOn5a/n5f+dV8PhdWvZXUsAC6MMgKg1BwOo8eW79A3e0/Iz9tDC+7trjbB/lbHAuDiKCMASsUYo2mrE/VJwjF5edg0+55IRYbVszoWADdAGQFQKq9/uV/vbjosm016ZUhnXd+6odWRALgJygiA3/TuxkN67Yv9kqTpv2+vP3RpbHEiAO6EMgLgV63ecUzPrU6UJE28uZWG92pmbSAAbocyAuCyNuw/oUkfJMgYaWSvMI27qYXVkQC4IcoIgEvaeeSM/vyv/96B95mB7WWzceM7ABWPMgLgIkknzmrUgq06l2/XNS3q69UhnbkDL4BKQxkBUMLxrFwNn7el+MZ3bw+Pkq8XN74DUHkoIwCKZZ4v0Mj5W3T0zHmFB9bSgnu7q7avl9WxALg5yggASVJugV0PLNpWfL+ZRaN7KLC2r9WxAFQDlBEAsjuMJn2QoC0HT6mOr5fevbeHQuvVtDoWgGqCMgJUc8YYTV+dqDW70uTj6aG3R0SqXQj3mwFQdSgjQDU3a91PenfTYUlFy7z3bh5ocSIA1Q1lBKjGPow7or99vleS9Mxt7TSwc4jFiQBUR5QRoJr6Zm+6nlyxU5L052sjNLpPuMWJAFRXlBGgGtp1JFMPLdkuu8NoUJcQPfm7NlZHAlCNUUaAaibl1Dndu7BoddU+LQL1tzs6y4PVVQFYiDICVCOnc/I1csEWZZzNU9tG/pp1Tzf5ePFrAIC1+C0EVBO5BXbdv2ibkk7kKCTATwvv7a46ft5WxwIAyghQHdgdRhOWJiju8Gn5+3lp4egeCvL3szoWAEiijABuzxijv/57tz5PLFrUbM6IKLUKqmN1LAAoVq4yMnPmTIWHh8vPz0+RkZHasGFDqY777rvv5OXlpS5dupTnZQGUw7xvD2rhxkOSihY1uzqivrWBAOACZS4jy5Yt04QJEzR16lTFx8erb9++6t+/v5KTk3/1uMzMTI0YMUI33XRTucMCKJvPdqXqhTV7JElTbm3DomYAnJLNGGPKckDPnj3VrVs3zZo1q3hb27ZtNWjQIMXExFz2uGHDhqlly5by9PTUxx9/rISEhFK/ZlZWlgICApSZmSl/f+6ZAZRG3OHTuvudzcordGhErzBN+3172WxM4QVQdUr7/l2mMyP5+fmKi4tTdHR0ie3R0dHauHHjZY9bsGCBfvrpJz377LOlep28vDxlZWWVeAAovUMZOXpg0TblFTp0c9uGenYgRQSA8ypTGcnIyJDdbldQUFCJ7UFBQUpLS7vkMfv379dTTz2lJUuWyMvLq1SvExMTo4CAgOJHaGhoWWIC1dqpnHzdu3CrTuXkq1OTAL1xV1d5sqgZACdWrgtYL/wLyxhzyb+67Ha77r77bk2bNk2tWrUq9defPHmyMjMzix8pKSnliQlUO7kFdv3fom06mJGjxnVraO7IKNX0Kd0fAQBglTL9lgoMDJSnp+dFZ0HS09MvOlsiSdnZ2dq2bZvi4+P18MMPS5IcDoeMMfLy8tLatWt14403XnScr6+vfH19yxINqPYcDqNHl+/QtsOnVcfPSwvv7a6GdVhLBIDzK9OZER8fH0VGRio2NrbE9tjYWPXu3fui/f39/bVr1y4lJCQUP8aMGaPWrVsrISFBPXv2vLL0AIr9fe1efbozVd6eNr09PFItWUsEgIso8/nbSZMmafjw4YqKilKvXr00Z84cJScna8yYMZKKPmI5evSoFi1aJA8PD3Xo0KHE8Q0bNpSfn99F2wGU37KtyZr1zU+SpJdv76TezQMtTgQApVfmMjJ06FCdPHlS06dPV2pqqjp06KA1a9YoLCxMkpSamvqba44AqDjfHcjQ1I9+kCSNv6mlBndrYnEiACibMq8zYgXWGQEubf/xbA2etVHZuYUa1CVE/xzahSm8AJxGpawzAsB5nMjO070Ltyo7t1Ddm12ll+/oRBEB4JIoI4ALyi2w64FF23Tk9Hk1q19Tbw+Pkq+Xp9WxAKBcKCOAi3E4jB79YIcSUs6obk1vzR/VXfVq+VgdCwDKjTICuJhXY/fp011FU3hn3xOpiAa1rY4EAFeEMgK4kBVxRzTj6wOSpJjBnXR1RH2LEwHAlaOMAC5iy8FTemrlTknSQ9c31x2RTOEF4B4oI4ALOHwyR3/+1zYV2I36dwjWY9GtrY4EABWGMgI4uczzBRq9cKtOnytQpyYBenVIF3lwF14AboQyAjixArtDD7+3XT+dyFGjAD/NHRGlGj5M4QXgXigjgJMyxmja6kRt2J+hmj6emjsySg39uQsvAPdDGQGc1KJNh7V4c7JsNun1YV3VPiTA6kgAUCkoI4ATWr/vhKb/e7ck6cnftdEt7YIsTgQAlYcyAjiZA+lnNfa97bI7jG7v1kR/vjbC6kgAUKkoI4ATOXMuX/e/W3Tzu6iwq/Ti4A7c/A6A26OMAE6iwO7QQ0u269DJc2pct4ZmD4/k5ncAqgXKCOAEjDF6blWiNv50UrV8PDVvVJQCa/taHQsAqgRlBHAC/9p8WEu+/+/MmTbB/lZHAoAqQxkBLPbdgQxNW/3fmTM3M3MGQDVDGQEsdPhkjh5aUjRzZnDXxsycAVAtUUYAi2TnFui+d7cp83yBOofW1YuDOzJzBkC1RBkBLGB3GI1fmqAD6WcV5O+rd4ZHys+bmTMAqifKCGCBv/9nr776MV2+Xh6aM5x7zgCo3igjQBX7OP6oZq/7SZL0tzs6qXNoXWsDAYDFKCNAFdp55IyeWLFTkvTg9c31hy6NLU4EANajjABVJD07V/+3KE75hQ7d2KahHotubXUkAHAKlBGgCuQV2vXg4u1Ky8pVRINaem1YF3l6MHMGACTKCFDpjDF69pNExR0+rTp+XnpnRJT8/bytjgUAToMyAlSyxZsPa+nWFNls0ht3dVXzBrWtjgQAToUyAlSizUkni5d6f6JfG93QuqHFiQDA+VBGgEpy5PQ5PbRkuwodRgM7h2jMdSz1DgCXQhkBKsH5fLv+/K84ncrJV/sQf/3t9k4s9Q4Al0EZASqYMUaTV+5U4rEs1avlo7eHR6qGD0u9A8DlUEaACjbv24P6OOGYPD1seuvubmpyVU2rIwGAU6OMABVo44EMxXz2oyRp6q1t1at5fYsTAYDzo4wAFeTI6XMa+9522R1Gg7s11r3XNLM6EgC4BMoIUAF+uWD19LkCdWwcoBf/2JELVgGglCgjwBX63wtW69fy0ezhkfLz5oJVACgtyghwhRZ8d6j4gtUZd3dT47o1rI4EAC6FMgJcgc1JJ/XCmj2SpClcsAoA5UIZAcopLTNXD/98weofuoRoNBesAkC5UEaAcsgrtGvM4jhlnM1Xm+A6emkwK6wCQHlRRoBymLZ6txJSzsjfz0tzhkexwioAXAHKCFBGy7Ym673vk2WzSa/f1VVN67PCKgBcCcoIUAY7Us7o6U8SJUmTbm6lG1o3tDgRALg+yghQSqdy8vXg4jjlFzp0c9sgjb2hhdWRAMAtUEaAUrA7jMa9H69jmbkKD6ylV4d2locHF6wCQEWgjACl8M/Yffr2QIZqeHtq9j2R8vfztjoSALgNygjwG2J3H9eMrw9Ikl66vaNaB9exOBEAuBfKCPArDmXkaNIHCZKkUb2b6Q9dGlsbCADcEGUEuIzz+UULm2XnFioy7CpNubWt1ZEAwC1RRoBLMMZo6ke79GNatgJr+2rmn7rJx4v/XQCgMvDbFbiExd8na2X80Z/vxNtVQf5+VkcCALdFGQEusCPljP66erck6cnftdbVEdyJFwAqE2UE+B+nc/L10JLtyrc71K99kB7oG2F1JABwe5QR4GcOh9HEDxJ09Mx5NatfU3+/szN34gWAKkAZAX424+sD+mbvCfl5e2gWC5sBQJUpVxmZOXOmwsPD5efnp8jISG3YsOGy+65cuVK33HKLGjRoIH9/f/Xq1Uv/+c9/yh0YqAwb9p/QP7/YJ0l6flBHtW3kb3EiAKg+ylxGli1bpgkTJmjq1KmKj49X37591b9/fyUnJ19y//Xr1+uWW27RmjVrFBcXpxtuuEEDBw5UfHz8FYcHKsKxM+c1fmmCjJHu6hGqOyKbWB0JAKoVmzHGlOWAnj17qlu3bpo1a1bxtrZt22rQoEGKiYkp1ddo3769hg4dqmeeeaZU+2dlZSkgIECZmZny9+cvVlSc/EKHhs3ZpO3JZ9Shsb8+HNNbft6eVscCALdQ2vfvMp0Zyc/PV1xcnKKjo0tsj46O1saNG0v1NRwOh7Kzs1WvXr3L7pOXl6esrKwSD6AyvPTZj9qefEb+fl6a9adIiggAWKBMZSQjI0N2u11BQUEltgcFBSktLa1UX+OVV15RTk6OhgwZctl9YmJiFBAQUPwIDQ0tS0ygVD7/IVXzvzsoSXplSBeF1qtpcSIAqJ7KdQHrhdMdjTGlmgL5/vvv67nnntOyZcvUsGHDy+43efJkZWZmFj9SUlLKExO4rEMZOXp8+U5J0p+vjdAt7YJ+4wgAQGXxKsvOgYGB8vT0vOgsSHp6+kVnSy60bNky3XfffVq+fLluvvnmX93X19dXvr6+ZYkGlFpugV0PLdmu7LxCdW92lR7r19rqSABQrZXpzIiPj48iIyMVGxtbYntsbKx69+592ePef/99jRo1Su+9954GDBhQvqRABZm2erd2p2apXi0fvXlXN3l7stwOAFipTGdGJGnSpEkaPny4oqKi1KtXL82ZM0fJyckaM2aMpKKPWI4ePapFixZJKioiI0aM0Ouvv66rr766+KxKjRo1FBAQUIFDAX7bR/FH9P6WZNls0uvDuig4gBvgAYDVylxGhg4dqpMnT2r69OlKTU1Vhw4dtGbNGoWFhUmSUlNTS6w58vbbb6uwsFBjx47V2LFji7ePHDlSCxcuvPIRAKW0/3i2pqz8QZI07saW6tuygcWJAABSOdYZsQLrjOBKncsv1B9mfKf96WfVp0Wg3h3dQ54e3HcGACpTpawzAriqZz5J1P70s2pYx1f/HNqFIgIAToQyAre3fFuKPow7Ig+b9MZdXdWgDjO1AMCZUEbg1vYdz9bTnxRdJzLplla6OqK+xYkAABeijMBtncsv1Ngl25Vb4FDfloF66PoWVkcCAFwCZQRu68LrRDy4TgQAnBJlBG7pwutEAmtznQgAOCvKCNzO/uPZeuaTRElcJwIAroAyArdyPt+use9t1/kCO9eJAICLoIzArUxbnah9x8+qAdeJAIDLoIzAbXyScFRLt6YU3XdmaBeuEwEAF0EZgVs4mJGjKSt3SZIeubGlercItDgRAKC0KCNweXmFdj383nbl5NvVM7yext/U0upIAIAyoIzA5cWs+VGJx7JUr5aPXh/WlfvOAICLoYzApf0nMU0LNx6SJL1yZ2cFB/hZGwgAUGaUEbisI6fP6fHlOyRJf742Qje0aWhxIgBAeVBG4JIK7Q5NWJqgrNxCdQ6tq8f6tbY6EgCgnCgjcEmvf7lf2w6fVh1fL824q6u8PflRBgBXxW9wuJyNBzI04+sDkqQXB3dUaL2aFicCAFwJyghcysmzeZqwLEHGSMO6h2pg5xCrIwEArhBlBC7D4TB6bPkOpWfnqUXD2np2YHurIwEAKgBlBC5j/ncH9fXeE/Lx8tCMu7uqho+n1ZEAABWAMgKXsOtIpl7+/EdJ0tO3tVObYH+LEwEAKgplBE7vbF6hHnl/uwrsRv3aB+menk2tjgQAqECUETi9Zz9J1KGT5xQS4KeXb+8km43l3gHAnVBG4NQ+STiqFduPyMMmvTasq+rW9LE6EgCgglFG4LRSTp3TXz76QZL08I0t1SO8nsWJAACVgTICp1Rgd2jc0nhl5xUqKuwqjbuxhdWRAACVhDICp/T6F/sVn3xGdfy89NqwLvJiuXcAcFv8hofT2fhTht76pmi595cGd1KTq1juHQDcGWUETuV0Tr4mLdshY6ShUaEa0KmR1ZEAAJWMMgKnYYzRUyt3Ki0rVxENaunZ37ezOhIAoApQRuA0lm5N0X8Sj8vb06Y3hnVVTR8vqyMBAKoAZQRO4acTZzV99W5J0uP9WqtD4wCLEwEAqgplBJbLL3Ro/NJ4nS+w65oW9XV/nwirIwEAqhBlBJZ7JXavfjiapbo1vfXKnV3k4cFy7wBQnVBGYKmNBzI0Z32SJOnl2zspOMDP4kQAgKpGGYFlTufka+IHCTJGuqtHU/VrH2x1JACABSgjsIQxRpNX7tLxrDxFNKilp29ra3UkAIBFKCOwxAfbUvR5YhrTeAEAlBFUvYMZOZr28zTeR6OZxgsA1R1lBFWqwO7QhGUJOpdv19UR9fRAX6bxAkB1RxlBlXrzqwPakXJG/n5eenVIF3kyjRcAqj3KCKpM3OFTmvHVfknSC3/sqJC6NSxOBABwBpQRVIns3AJNWJYgh5H+2LWxBnYOsToSAMBJUEZQJZ5btVspp86rcd0amvaH9lbHAQA4EcoIKt2aXalasf2IPGzSa8O6yN/P2+pIAAAnQhlBpTqelaspH+2SJD14fXN1b1bP4kQAAGdDGUGlcTiMHlu+Q2fOFahj4wBNuLmV1ZEAAE6IMoJKs2jTIW3YnyFfLw/9c2hneXvy4wYAuBjvDqgU+49nK+azHyVJU25tqxYN61icCADgrCgjqHD5hQ5N/CBBeYUOXduqgUb0CrM6EgDAiVFGUOFe/3Kffjiapbo1vfX3OzrJZmOVVQDA5VFGUKG2HTqlWd/8JEl68Y8dFeTvZ3EiAICzo4ygwpzKydfED4pWWR3crbFu7djI6kgAABdAGUGFOJ9v1/3vblXKqfNqclUNPfd7VlkFAJQOZQRXzO4wGr80XtuTi+7Gu2BUd1ZZBQCUWrnKyMyZMxUeHi4/Pz9FRkZqw4YNv7r/unXrFBkZKT8/P0VERGj27NnlCgvnY4zRtNWJWrv7uHw8PTR3ZHe1DGIaLwCg9MpcRpYtW6YJEyZo6tSpio+PV9++fdW/f38lJydfcv+DBw/q1ltvVd++fRUfH68pU6Zo3LhxWrFixRWHh/XmrE/Sok2HJUn/HNpFPcJZ7h0AUDY2Y4wpywE9e/ZUt27dNGvWrOJtbdu21aBBgxQTE3PR/k8++aRWrVqlPXv2FG8bM2aMduzYoU2bNpXqNbOyshQQEKDMzEz5+/uXJS4q0ScJRzV+aYIk6S8D2ur+vhHWBgIAOJXSvn97leWL5ufnKy4uTk899VSJ7dHR0dq4ceMlj9m0aZOio6NLbOvXr5/mzZungoICeXtffG1BXl6e8vLySgymMqyIO6IfjmVWytd2d4V2o6Vbi86G3dcnnCICACi3MpWRjIwM2e12BQUFldgeFBSktLS0Sx6TlpZ2yf0LCwuVkZGhRo0unv4ZExOjadOmlSVauazbd0Krdhyr9NdxZwM6NtLUW9taHQMA4MLKVEZ+ceGKmsaYX11l81L7X2r7LyZPnqxJkyYV/zsrK0uhoaHlifqrbmkXpNB6NSr861YXwQE1dGdkE3l4sMIqAKD8ylRGAgMD5enpedFZkPT09IvOfvwiODj4kvt7eXmpfv36lzzG19dXvr6+ZYlWLgM7h2hg55BKfx0AAHB5ZZpN4+Pjo8jISMXGxpbYHhsbq969e1/ymF69el20/9q1axUVFXXJ60UAAED1UuapvZMmTdLcuXM1f/587dmzRxMnTlRycrLGjBkjqegjlhEjRhTvP2bMGB0+fFiTJk3Snj17NH/+fM2bN0+PPfZYxY0CAAC4rDJfMzJ06FCdPHlS06dPV2pqqjp06KA1a9YoLKzoNvGpqakl1hwJDw/XmjVrNHHiRL311lsKCQnRG2+8odtvv73iRgEAAFxWmdcZsQLrjAAA4HpK+/7NvWkAAIClKCMAAMBSlBEAAGApyggAALAUZQQAAFiKMgIAACxFGQEAAJaijAAAAEtRRgAAgKXKvBy8FX5ZJDYrK8viJAAAoLR+ed/+rcXeXaKMZGdnS5JCQ0MtTgIAAMoqOztbAQEBl33eJe5N43A4dOzYMdWpU0c2m61cXyMrK0uhoaFKSUlxy/vbMD7X5+5jdPfxSe4/Rsbn+qp6jMYYZWdnKyQkRB4el78yxCXOjHh4eKhJkyYV8rX8/f3d9odMYnzuwN3H6O7jk9x/jIzP9VXlGH/tjMgvuIAVAABYijICAAAsVW3KiK+vr5599ln5+vpaHaVSMD7X5+5jdPfxSe4/Rsbn+px1jC5xASsAAHBf1ebMCAAAcE6UEQAAYCnKCAAAsBRlBAAAWKpalJGZM2cqPDxcfn5+ioyM1IYNG6yOVG7r16/XwIEDFRISIpvNpo8//rjE88YYPffccwoJCVGNGjV0/fXXKzEx0Zqw5RATE6Pu3burTp06atiwoQYNGqS9e/eW2MeVxzhr1ix16tSpeMGhXr166bPPPit+3pXHdikxMTGy2WyaMGFC8TZXH+Nzzz0nm81W4hEcHFz8vKuPT5KOHj2qe+65R/Xr11fNmjXVpUsXxcXFFT/v6mNs1qzZRd9Dm82msWPHSnL98RUWFuovf/mLwsPDVaNGDUVERGj69OlyOBzF+zjdGI2bW7p0qfH29jbvvPOO2b17txk/frypVauWOXz4sNXRymXNmjVm6tSpZsWKFUaS+eijj0o8/9JLL5k6deqYFStWmF27dpmhQ4eaRo0amaysLGsCl1G/fv3MggULzA8//GASEhLMgAEDTNOmTc3Zs2eL93HlMa5atcp8+umnZu/evWbv3r1mypQpxtvb2/zwww/GGNce24W2bNlimjVrZjp16mTGjx9fvN3Vx/jss8+a9u3bm9TU1OJHenp68fOuPr5Tp06ZsLAwM2rUKPP999+bgwcPmi+++MIcOHCgeB9XH2N6enqJ719sbKyRZL7++mtjjOuP7/nnnzf169c3//73v83BgwfN8uXLTe3atc1rr71WvI+zjdHty0iPHj3MmDFjSmxr06aNeeqppyxKVHEuLCMOh8MEBwebl156qXhbbm6uCQgIMLNnz7Yg4ZVLT083ksy6deuMMe45xquuusrMnTvXrcaWnZ1tWrZsaWJjY811111XXEbcYYzPPvus6dy58yWfc4fxPfnkk6ZPnz6Xfd4dxnih8ePHm+bNmxuHw+EW4xswYIAZPXp0iW2DBw8299xzjzHGOb+Hbv0xTX5+vuLi4hQdHV1ie3R0tDZu3GhRqspz8OBBpaWllRivr6+vrrvuOpcdb2ZmpiSpXr16ktxrjHa7XUuXLlVOTo569erlVmMbO3asBgwYoJtvvrnEdncZ4/79+xUSEqLw8HANGzZMSUlJktxjfKtWrVJUVJTuvPNONWzYUF27dtU777xT/Lw7jPF/5efna/HixRo9erRsNptbjK9Pnz768ssvtW/fPknSjh079O233+rWW2+V5JzfQ5e4UV55ZWRkyG63KygoqMT2oKAgpaWlWZSq8vwypkuN9/Dhw1ZEuiLGGE2aNEl9+vRRhw4dJLnHGHft2qVevXopNzdXtWvX1kcffaR27doV/xJw5bFJ0tKlS7V9+3Zt3br1oufc4fvXs2dPLVq0SK1atdLx48f1/PPPq3fv3kpMTHSL8SUlJWnWrFmaNGmSpkyZoi1btmjcuHHy9fXViBEj3GKM/+vjjz/WmTNnNGrUKEnu8TP65JNPKjMzU23atJGnp6fsdrteeOEF3XXXXZKcc4xuXUZ+YbPZSvzbGHPRNnfiLuN9+OGHtXPnTn377bcXPefKY2zdurUSEhJ05swZrVixQiNHjtS6deuKn3flsaWkpGj8+PFau3at/Pz8LrufK4+xf//+xf/dsWNH9erVS82bN9e7776rq6++WpJrj8/hcCgqKkovvviiJKlr165KTEzUrFmzNGLEiOL9XHmM/2vevHnq37+/QkJCSmx35fEtW7ZMixcv1nvvvaf27dsrISFBEyZMUEhIiEaOHFm8nzON0a0/pgkMDJSnp+dFZ0HS09MvaoTu4Jcr+t1hvI888ohWrVqlr7/+Wk2aNCne7g5j9PHxUYsWLRQVFaWYmBh17txZr7/+uluMLS4uTunp6YqMjJSXl5e8vLy0bt06vfHGG/Ly8ioehyuP8UK1atVSx44dtX//frf4HjZq1Ejt2rUrsa1t27ZKTk6W5B7/D/7i8OHD+uKLL3T//fcXb3OH8T3++ON66qmnNGzYMHXs2FHDhw/XxIkTFRMTI8k5x+jWZcTHx0eRkZGKjY0tsT02Nla9e/e2KFXlCQ8PV3BwcInx5ufna926dS4zXmOMHn74Ya1cuVJfffWVwsPDSzzvDmO8kDFGeXl5bjG2m266Sbt27VJCQkLxIyoqSn/605+UkJCgiIgIlx/jhfLy8rRnzx41atTILb6H11xzzUXT6fft26ewsDBJ7vX/4IIFC9SwYUMNGDCgeJs7jO/cuXPy8Cj59u7p6Vk8tdcpx2jJZbNV6JepvfPmzTO7d+82EyZMMLVq1TKHDh2yOlq5ZGdnm/j4eBMfH28kmVdffdXEx8cXT1V+6aWXTEBAgFm5cqXZtWuXueuuu1xqStqDDz5oAgICzDfffFNi6t25c+eK93HlMU6ePNmsX7/eHDx40OzcudNMmTLFeHh4mLVr1xpjXHtsl/O/s2mMcf0xPvroo+abb74xSUlJZvPmzea2224zderUKf6d4urj27Jli/Hy8jIvvPCC2b9/v1myZImpWbOmWbx4cfE+rj5GY4yx2+2madOm5sknn7zoOVcf38iRI03jxo2Lp/auXLnSBAYGmieeeKJ4H2cbo9uXEWOMeeutt0xYWJjx8fEx3bp1K54m6oq+/vprI+mix8iRI40xRVO2nn32WRMcHGx8fX3Ntddea3bt2mVt6DK41NgkmQULFhTv48pjHD16dPHPYoMGDcxNN91UXESMce2xXc6FZcTVx/jLegze3t4mJCTEDB482CQmJhY/7+rjM8aY1atXmw4dOhhfX1/Tpk0bM2fOnBLPu8MY//Of/xhJZu/evRc95+rjy8rKMuPHjzdNmzY1fn5+JiIiwkydOtXk5eUV7+NsY7QZY4wlp2QAAADk5teMAAAA50cZAQAAlqKMAAAAS1FGAACApSgjAADAUpQRAABgKcoIAACwFGUEAABYijICAAAsRRkBAACWoowAAABLUUYAAICl/h/Je1b/LO1t2gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(vct['height'].values,vct['vct_b'].values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "1\n",
      "1.0\n",
      "22000.0\n",
      "\n",
      "2\n",
      "0.9875\n",
      "21453.4375\n",
      "\n",
      "3\n",
      "0.975\n",
      "20913.750000000004\n",
      "\n",
      "4\n",
      "0.9625\n",
      "20380.937499999996\n",
      "\n",
      "5\n",
      "0.95\n",
      "19855.0\n",
      "\n",
      "6\n",
      "0.9375\n",
      "19335.9375\n",
      "\n",
      "7\n",
      "0.925\n",
      "18823.75\n",
      "\n",
      "8\n",
      "0.9125\n",
      "18318.4375\n",
      "\n",
      "9\n",
      "0.9\n",
      "17820.0\n",
      "\n",
      "10\n",
      "0.8875\n",
      "17328.4375\n",
      "\n",
      "11\n",
      "0.875\n",
      "16843.75\n",
      "\n",
      "12\n",
      "0.8625\n",
      "16365.9375\n",
      "\n",
      "13\n",
      "0.85\n",
      "15895.000000000002\n",
      "0.006562499999999886\n",
      "\n",
      "14\n",
      "0.8375\n",
      "15430.937499999998\n",
      "0.03556640625000011\n",
      "\n",
      "15\n",
      "0.825\n",
      "14973.75\n",
      "0.064140625\n",
      "\n",
      "16\n",
      "0.8125\n",
      "14523.4375\n",
      "0.09228515625\n",
      "\n",
      "17\n",
      "0.8\n",
      "14080.0\n",
      "0.12\n",
      "\n",
      "18\n",
      "0.7875\n",
      "13643.437500000002\n",
      "0.14728515624999988\n",
      "\n",
      "19\n",
      "0.775\n",
      "13213.749999999998\n",
      "0.1741406250000001\n",
      "\n",
      "20\n",
      "0.7625\n",
      "12790.9375\n",
      "0.20056640625\n",
      "\n",
      "21\n",
      "0.75\n",
      "12375.0\n",
      "0.2265625\n",
      "\n",
      "22\n",
      "0.7375\n",
      "11965.937500000002\n",
      "0.2521289062499999\n",
      "\n",
      "23\n",
      "0.725\n",
      "11563.750000000002\n",
      "0.2772656249999999\n",
      "\n",
      "24\n",
      "0.7125\n",
      "11168.437499999998\n",
      "0.30197265625000014\n",
      "\n",
      "25\n",
      "0.7\n",
      "10779.999999999998\n",
      "0.3262500000000001\n",
      "\n",
      "26\n",
      "0.6875\n",
      "10398.4375\n",
      "0.35009765625\n",
      "\n",
      "27\n",
      "0.675\n",
      "10023.750000000002\n",
      "0.3735156249999999\n",
      "\n",
      "28\n",
      "0.6625\n",
      "9655.937500000002\n",
      "0.39650390624999987\n",
      "\n",
      "29\n",
      "0.65\n",
      "9294.999999999998\n",
      "0.4190625000000001\n",
      "\n",
      "30\n",
      "0.6375\n",
      "8940.937499999998\n",
      "0.44119140625000014\n",
      "\n",
      "31\n",
      "0.625\n",
      "8593.75\n",
      "0.462890625\n",
      "\n",
      "32\n",
      "0.6125\n",
      "8253.437500000002\n",
      "0.48416015624999986\n",
      "\n",
      "33\n",
      "0.6\n",
      "7920.000000000001\n",
      "0.5049999999999999\n",
      "\n",
      "34\n",
      "0.5875\n",
      "7593.437499999999\n",
      "0.52541015625\n",
      "\n",
      "35\n",
      "0.575\n",
      "7273.749999999998\n",
      "0.5453906250000001\n",
      "\n",
      "36\n",
      "0.5625\n",
      "6960.9375\n",
      "0.56494140625\n",
      "\n",
      "37\n",
      "0.55\n",
      "6655.000000000002\n",
      "0.5840624999999999\n",
      "\n",
      "38\n",
      "0.5375\n",
      "6355.937500000001\n",
      "0.60275390625\n",
      "\n",
      "39\n",
      "0.525\n",
      "6063.749999999999\n",
      "0.621015625\n",
      "\n",
      "40\n",
      "0.5125\n",
      "5778.437499999998\n",
      "0.6388476562500001\n",
      "\n",
      "41\n",
      "0.5\n",
      "5500.0\n",
      "0.65625\n",
      "\n",
      "42\n",
      "0.4875\n",
      "5228.437500000001\n",
      "0.67322265625\n",
      "\n",
      "43\n",
      "0.475\n",
      "4963.750000000001\n",
      "0.689765625\n",
      "\n",
      "44\n",
      "0.4625\n",
      "4705.937499999999\n",
      "0.70587890625\n",
      "\n",
      "45\n",
      "0.45\n",
      "4455.0\n",
      "0.7215625\n",
      "\n",
      "46\n",
      "0.4375\n",
      "4210.9375\n",
      "0.73681640625\n",
      "\n",
      "47\n",
      "0.425\n",
      "3973.7500000000005\n",
      "0.751640625\n",
      "\n",
      "48\n",
      "0.4125\n",
      "3743.437500000001\n",
      "0.76603515625\n",
      "\n",
      "49\n",
      "0.4\n",
      "3519.999999999999\n",
      "0.78\n",
      "\n",
      "50\n",
      "0.3875\n",
      "3303.4374999999995\n",
      "0.79353515625\n",
      "\n",
      "51\n",
      "0.375\n",
      "3093.75\n",
      "0.806640625\n",
      "\n",
      "52\n",
      "0.3625\n",
      "2890.9375000000005\n",
      "0.81931640625\n",
      "\n",
      "53\n",
      "0.35\n",
      "2695.0000000000005\n",
      "0.8315625\n",
      "\n",
      "54\n",
      "0.3375\n",
      "2505.9374999999995\n",
      "0.84337890625\n",
      "\n",
      "55\n",
      "0.325\n",
      "2323.7499999999995\n",
      "0.854765625\n",
      "\n",
      "56\n",
      "0.3125\n",
      "2148.4375\n",
      "0.86572265625\n",
      "\n",
      "57\n",
      "0.3\n",
      "1980.0000000000002\n",
      "0.87625\n",
      "\n",
      "58\n",
      "0.2875\n",
      "1818.4375000000002\n",
      "0.88634765625\n",
      "\n",
      "59\n",
      "0.275\n",
      "1663.7499999999998\n",
      "0.896015625\n",
      "\n",
      "60\n",
      "0.2625\n",
      "1515.9374999999998\n",
      "0.90525390625\n",
      "\n",
      "61\n",
      "0.25\n",
      "1375.0\n",
      "0.9140625\n",
      "\n",
      "62\n",
      "0.2375\n",
      "1240.9375000000002\n",
      "0.92244140625\n",
      "\n",
      "63\n",
      "0.225\n",
      "1113.7500000000005\n",
      "0.930390625\n",
      "\n",
      "64\n",
      "0.2125\n",
      "993.4374999999995\n",
      "0.93791015625\n",
      "\n",
      "65\n",
      "0.2\n",
      "879.9999999999998\n",
      "0.945\n",
      "\n",
      "66\n",
      "0.1875\n",
      "773.4375\n",
      "0.95166015625\n",
      "\n",
      "67\n",
      "0.175\n",
      "673.7500000000001\n",
      "0.957890625\n",
      "\n",
      "68\n",
      "0.1625\n",
      "580.9375000000003\n",
      "0.96369140625\n",
      "\n",
      "69\n",
      "0.15\n",
      "494.9999999999997\n",
      "0.9690625\n",
      "\n",
      "70\n",
      "0.1375\n",
      "415.93749999999994\n",
      "0.97400390625\n",
      "\n",
      "71\n",
      "0.125\n",
      "343.75\n",
      "0.978515625\n",
      "\n",
      "72\n",
      "0.1125\n",
      "278.4375000000001\n",
      "0.98259765625\n",
      "\n",
      "73\n",
      "0.1\n",
      "220.0000000000002\n",
      "0.98625\n",
      "\n",
      "74\n",
      "0.0875\n",
      "168.4374999999998\n",
      "0.98947265625\n",
      "\n",
      "75\n",
      "0.075\n",
      "123.74999999999993\n",
      "0.992265625\n",
      "\n",
      "76\n",
      "0.0625\n",
      "85.9375\n",
      "0.99462890625\n",
      "\n",
      "77\n",
      "0.05\n",
      "55.00000000000005\n",
      "0.9965625\n",
      "\n",
      "78\n",
      "0.0375\n",
      "30.937500000000075\n",
      "0.99806640625\n",
      "\n",
      "79\n",
      "0.025\n",
      "13.74999999999995\n",
      "0.999140625\n",
      "\n",
      "80\n",
      "0.0125\n",
      "3.4374999999999876\n",
      "0.99978515625\n",
      "\n",
      "81\n",
      "0.0\n",
      "0.0\n",
      "1.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x149ca1e82b00>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+BElEQVR4nO3deVhVdeLH8c/lAheQRQVFUURwRXEFNSzbo1Fz2jMrl6wmp02zVW2yHBuapmmaFm0sl5zM/JUtVlTS5pLmQrjviuICIqjs673n9wfJRGoJAod77/v1PPcpzz0HPt9H5X4853u+x2IYhiEAAACTeJgdAAAAuDfKCAAAMBVlBAAAmIoyAgAATEUZAQAApqKMAAAAU1FGAACAqSgjAADAVJ5mBzgXDodDR44cUUBAgCwWi9lxAADAOTAMQ/n5+QoLC5OHx9nPfzhFGTly5IjCw8PNjgEAAGrh4MGDatu27Vnfd4oyEhAQIKlyMIGBgSanAQAA5yIvL0/h4eFVn+Nn4xRl5NSlmcDAQMoIAABO5vemWDCBFQAAmIoyAgAATEUZAQAApqKMAAAAU1FGAACAqSgjAADAVJQRAABgKsoIAAAwFWUEAACYqsZlZPny5Ro2bJjCwsJksVj08ccf/+4xy5YtU2xsrHx8fBQVFaU33nijNlkBAIALqnEZKSwsVK9evfTaa6+d0/5paWkaMmSIBg0apNTUVE2ePFkPPfSQFi9eXOOwAADA9dT42TSDBw/W4MGDz3n/N954Q+3atdPLL78sSYqOjtb69ev14osv6sYbb6zptwcAAC6m3h+Ut3r1aiUkJFTbdvXVV2v27NkqLy+Xl5fXaceUlpaqtLS06td5eXn1HRMAAKdlGIYKSiuUV1KhvOJy5RaXq7C0QsXldhWX2VVS4VBJmV3F5XZV2B2yG4bsDsnucMjukByGoZti2yqmTZAp+eu9jGRmZio0NLTattDQUFVUVCg7O1utW7c+7ZjExEQ9++yz9R0NAIBGraisQkdOlujIyWJl5BYrK69UOYVlOlZQquz8UmUXVP46r7hcDuP8vldsRDPXLSPS6Y8ONgzjjNtPmTRpkiZOnFj167y8PIWHh9dfQAAATGB3GMrILdaBnCLtzyms/G92oQ6eqCwfJ4vKa/T1vKwWBfl6KdDXS/42T/l6WeXrba38r5dVNi+rvKwWeVgssnpUvir/X+rY0r+eRvn76r2MtGrVSpmZmdW2ZWVlydPTU8HBwWc8xmazyWaz1Xc0AAAaRIXdoQPHi7T7aL52HS3Q7qwC7T6ar33HClVmd/zmsQE2T4U19VXrpj4KDfBRsL+3QvxtCgmwKeTn/w/y9VKQr5dsnh5n/Yd+Y1bvZSQ+Pl6ffvpptW1Lly5VXFzcGeeLAADgzErK7dqekaetR069crUjM19lFWcuHV5Wi8Kb+6l9cBNFBFf+t11zv6oCEujj+p+VNS4jBQUF2rNnT9Wv09LStGHDBjVv3lzt2rXTpEmTdPjwYc2fP1+SNG7cOL322muaOHGi7rnnHq1evVqzZ8/WwoUL624UAACYwOEwtPdYgVLTTyr14Emlpp/Q7qwC2c8wgcPXy6pOof7q2NJfnUMD1DnUXx1bBKhNM19ZPZzvbEZdqnEZWb9+vS677LKqX5+a2zF69GjNmzdPGRkZSk9Pr3o/MjJSSUlJevjhh/X6668rLCxMr7zyCrf1AgCcTkm5XRsOntTatONat/+4NqSfVH5pxWn7hfh7q1tYkGLCAtU9LEjdwwLVrrmfPNy8dJyNxTg1m7QRy8vLU1BQkHJzcxUYGGh2HACAmygpt2v9/hP6cV+O1qYd14aDJ0+b4+HrZVXPtkHq3a6p+oQ3U+/wpgoNtDnl3I26dq6f3w1yNw0AAM7A4TC0LSNPK/dka+XubK3bf1ylv5rr0SLApgGRzdU/srliI5qpS2iAPK086u18UEYAAG4tt6hcy3Yf07fbj2rF7mzlFJZVez800KYLO4RoQFRz9Y8MVvtgP8561DHKCADArRiGoT1ZBfpmR5a+3ZGllAMnqk04beJt1QVRwbqoU4gGdQpRhxb+lI96RhkBALg8wzC08VCuvtySqa+2Ziotu7Da+51D/XVZ15a6vEtL9Y1oJi8uuzQoyggAwCU5HIbW7j9eVUAyckuq3vP29NDADsG6vGtLXdalpcKb+5mYFJQRAIDLOHUG5NONR/TZpiM6mve/h6428bbqsq4t9YeYVrq0S0v52/gIbCz4nQAAOL09WQX6OPWwPt10RAdyiqq2B/h46ururTQ4ppUu7BgiHy+riSlxNpQRAIBTOllUpk83HtEHKYe08VBu1XZfL6uu7BaqP/YK08WdQ2TzpIA0dpQRAIDTqLA7tGzXMX2QckjfbM+qWoDM6mHRpZ1b6No+bXRldEv5efPx5kz43QIANHqHThTp/9Yd1KL1B6vNA+nWOlA3xrbVtb3DFOLP096dFWUEANAoldsd+mZ7lhauTdfy3cd06uElzZt46/o+bXRj37bqFsYjQlwBZQQA0Khk5Zdo4ZqDWrDmgLLy/3cW5MKOwRrRv52u6hbKPBAXQxkBAJjOMAz9lH5S81fvV9LmDJXbK0+DhPh76+a4cN3aL1wRwU1MTon6QhkBAJimrMKhzzYd0dwf9mvz4f/dEdO3XVONHtheg2Nay9uT1VBdHWUEANDgcovKtWDtAb29an/VhFRvTw9d2ytMo+Lbq0fbIJMToiFRRgAADSY9p0hzfkjT/60/qKIyuySpZYBNowe214j+7dS8ibfJCWEGyggAoN5tz8jTjO/36vNNR3TqAbldWwXo7kFRGtarNRNS3RxlBABQb9bvP64Z3+/VtzuyqrZd3LmF/jQoShd2DJbFYjExHRoLyggAoE4ZhqEVu7P12nd7tDbtuCTJwyIN6dFaf760g7qHMR8E1VFGAAB14lQJefnrXfop/aQkyctq0Y192+reSzooMoRbc3FmlBEAwHkxDEPLfy4hqT+XEJunh24fEKE/XRylVkE+5gZEo0cZAQDU2g97svXi0p3VSsgdF0To3kui1DKAEoJzQxkBANTYT+kn9OJXO7Vqb44kSgjOD2UEAHDOdmbm68WlO5W87aikyjkhtw+I0H2XdlDLQEoIaocyAgD4XYdOFOmlpbv00YbDMozKu2Nu7NtW46/spLbN/MyOBydHGQEAnFVucblmfLdHc1ftV1mFQ5I0pEcrTbyqizq29Dc5HVwFZQQAcJqyCof+++MBvfrtbp0sKpckxUcFa9KQrurZtqm54eByKCMAgCqGYeiLLZl6/osdSj9eJEnq1NJfk4Z01WVdWrJiKuoFZQQAIEnaeiRXz366rWrV1BYBNk28qrNujm0rT6uHyengyigjAODmsgtK9c+lu/TeunQZRuVtuvde0kH3XhylJjY+JlD/+FMGAG6q3O7Q26v2699f71Z+aYUkaVivMD05uKvaNPU1OR3cCWUEANzQqr3ZevqTrdqTVSBJimkTqKev6a7+kc1NTgZ3RBkBADdyNK9Ez32+XUs2HpEkBTfx1uN/6KKbYsNl9WByKsxBGQEAN3Dqksy/knepsMwuD4t0xwUReuSqLgry8zI7HtwcZQQAXFzKgROa8tFm7cjMlyT1Dm+q6dfFKKZNkMnJgEqUEQBwUbnF5frHVzu0YE3lXTLN/Lz05OCuujk2XB5ckkEjQhkBABdjGIaSNmfqmU+36lh+qSTp5ti2mjwkWs2aeJucDjgdZQQAXMihE0V6+pOt+nZHliQpKqSJnru+h+I7BJucDDg7yggAuACHw9A7aw7o71/sUGGZXd5WD/350g6677IOsnlazY4H/CbKCAA4uX3HCvTk4s1au79yGfd+7Zsp8YaePFUXToMyAgBOqsLu0OyVaXopeZdKKxzy87bqiT901cgLIpigCqdCGQEAJ7T7aL4efX+jNh7KlSQN6hSiv13fQ+HN/UxOBtQcZQQAnIjdYWj2yn16cekulVU4FOjjqaeu6aabY9vKYuFsCJwTZQQAnMT+7EI9+v5GrT9wQpJ0WZcWev7GngoN9DE5GXB+KCMA0MidulMmMWmHisvt8rd56i/XROuWuHDOhsAlUEYAoBHLyC3WY+9v0so92ZKk+Khg/ePmnmrbjLkhcB2UEQBopD7flKHJH21WbnG5fLw8NGlwNHfKwCVRRgCgkckvKdczS7Zp8U+HJEk92wbpX8N7q0ML1g2Ba6KMAEAjsn7/cU1YtEGHThTLwyLdd2lHjb+yk7ysHmZHA+oNZQQAGoEKu0OvfLtHr327Ww5DatvMV/8a3lv92jc3OxpQ7ygjAGCyIyeLNf69VK3bX3nL7g192+jZP3ZXgI+XycmAhkEZAQATfbklU08s3qTc4nL52zz13PUxurZ3G7NjAQ2KMgIAJigpt+u5z7frvz8ekCT1ahukV0b0UURwE5OTAQ2PMgIADWzvsQLdv+An7cjMlyTde3GUHknoIm9PJqnCPVFGAKABfbLhsCZ9uFlFZXaF+Hvrn7f01iWdW5gdCzAVZQQAGkBJuV3TPtumd9ekS5IuiGquV0b0UcsAnisDUEYAoJ7tzy7UfQt+0raMPFks0oOXddT4KzvLykqqgCSpVhcoZ8yYocjISPn4+Cg2NlYrVqz4zf0XLFigXr16yc/PT61bt9add96pnJycWgUGAGfyxeYMXfPqSm3LyFPzJt56+87+mpjQhSIC/EKNy8iiRYs0YcIETZkyRampqRo0aJAGDx6s9PT0M+6/cuVKjRo1SnfddZe2bt2q999/X+vWrdPdd9993uEBoLGqsDv0t6Tt+vOCn1RQWqF+7Zsp6aFBupj5IcBpLIZhGDU5YMCAAerbt69mzpxZtS06OlrXXXedEhMTT9v/xRdf1MyZM7V3796qba+++qpeeOEFHTx48Jy+Z15enoKCgpSbm6vAwMCaxAWABncsv1QPLvxJP+47Lkn608VRevzqLvJkSXe4mXP9/K7R34yysjKlpKQoISGh2vaEhAStWrXqjMcMHDhQhw4dUlJSkgzD0NGjR/XBBx9o6NChZ/0+paWlysvLq/YCAGeQcuCErnl1hX7cd1xNvK2aeXtfTR4STREBfkON/nZkZ2fLbrcrNDS02vbQ0FBlZmae8ZiBAwdqwYIFGj58uLy9vdWqVSs1bdpUr7766lm/T2JiooKCgqpe4eHhNYkJAA3OMAzNX71ft85araN5perY0l+fPHCRBvdobXY0oNGrVVW3WKpPvDIM47Rtp2zbtk0PPfSQnn76aaWkpOjLL79UWlqaxo0bd9avP2nSJOXm5la9zvVyDgCYoaTcrkff36SnP9mqcruhoT1a6+P7L1THlv5mRwOcQo1u7Q0JCZHVaj3tLEhWVtZpZ0tOSUxM1IUXXqjHHntMktSzZ081adJEgwYN0vTp09W69en/arDZbLLZbDWJBgCmyMgt1rj/pmjjoVxZPSx68g9ddfegyLP+Aw3A6Wp0ZsTb21uxsbFKTk6utj05OVkDBw484zFFRUXy8Kj+baxWq6TKMyoA4KzW7z+uYa/+oI2HctXUz0vzx/bXPRdHUUSAGqrxomcTJ07UyJEjFRcXp/j4eM2aNUvp6elVl10mTZqkw4cPa/78+ZKkYcOG6Z577tHMmTN19dVXKyMjQxMmTFD//v0VFhZWt6MBgAby7pp0TV2yReV2Q11bBejNUXEKb+5ndizAKdW4jAwfPlw5OTmaNm2aMjIyFBMTo6SkJEVEREiSMjIyqq05MmbMGOXn5+u1117TI488oqZNm+ryyy/X3//+97obBQA0kHK7Q88s2aoFPy/rPrRHa/3j5p7y82ZBa6C2arzOiBlYZwRAY3CisEx/XpCiH/cdl8UiPZrQRfdd2oHLMsBZnOvnN1UeAM7Bnqx83fX2eh3IKZK/zVP/vrW3rog+88R9ADVDGQGA3/H9ziw9+G6q8ksrFN7cV7NH91Pn0ACzYwEugzICAGdhGIbm/LBfz32+TQ5D6t++uWbe0VfB/iw9ANQlyggAnEG53aGnP9mqhWsrJ6reEtdW06/rIW9PlnUH6hplBAB+Ja+kXPcv+EkrdmfLYpGmDInWXRexkBlQXygjAPALh04Uaey8ddp1tEC+Xla9MqKPrurGRFWgPlFGAOBnGw+e1F1vr1d2QalaBtg0Z0w/xbQJMjsW4PIoIwAg6cstGZqwaINKyh3q2ipAc8b0U1hTX7NjAW6BMgLA7c1emabpn2+TYUiXdWmhV2/rK38bPx6BhsLfNgBuy+Ew9FzSds1emSZJGnlBhKYO6yZPK3fMAA2JMgLALZWU2/XI/23U55szJEmTBnfVn3jiLmAKyggAt3OyqEx/mp+itfuPy8tq0Ys399K1vduYHQtwW5QRAG7l0IkijZm7TnuyChRg89R/RsVqYIcQs2MBbo0yAsBtbM/I0+g5a5WVX6rWQT6ae2c/dW3Fk8ABs1FGALiFNftydPf89covqVCX0ADNG9tPrYO4dRdoDCgjAFzeV1sz9eDCVJVVONSvfTO9Nbqfgny9zI4F4GeUEQAu7b216Zr80WY5DOnK6FC9dlsf+XhZzY4F4BcoIwBckmEYmvH9Xv3jq52SKp+6+7fre7CGCNAIUUYAuByHw9BfP9+muT/slyTdf1kHPZrQhTVEgEaKMgLApVTYHXpi8WYt/umQJOnpa7pp7EWRJqcC8FsoIwBcRmmFXQ8tTNVXW4/K6mHRP27qqRv6tjU7FoDfQRkB4BIKSyt0739TtHJPtrytHnrttj5K6N7K7FgAzgFlBIDTO1lUpjvnrVNq+kn5eVv11qg4DezIqqqAs6CMAHBqWfklGjV7rXZk5qupn5fm3dlfvcObmh0LQA1QRgA4rYzcYt3+5hrtyy5UywCb/nvXAHVpFWB2LAA1RBkB4JQOHi/SiDd/1KETxWrT1Ffv3jNAEcFNzI4FoBYoIwCczt5jBbr9zTXKzCtR+2A/LbjnArVpynNmAGdFGQHgVHZk5umOt9You6BMnVr6a8HdA9Qy0MfsWADOA2UEgNPYfChXI+es0cmickW3DtQ7d/VXsL/N7FgAzhNlBIBTSE0/oVFz1iq/pEK9wptq/p39FeTHk3cBV0AZAdDopRw4oTFz1iq/tEL92jfTnDH9FOBDEQFcBWUEQKO2fv9xjZ6zVoVldvWPbK65Y/qpiY0fXYAr4W80gEZrzb4c3TlvnYrK7IqPCtbsMXHy8+bHFuBq+FsNoFFavTdHY+etU3G5XRd1DNGbo+Lk6201OxaAekAZAdDorNqTrbFvr1NJuUMXd26hWSNj5eNFEQFclYfZAQDgl1bvzakqIpd1oYgA7oAzIwAajTX7Ki/NlJQ7dGmXFnpjZKxsnhQRwNVxZgRAo7Bu/3Hd+fMckYs7t9Abd1BEAHdBGQFgupQDxzVmzloVldk1qFMIl2YAN0MZAWCqn9JPaPScdSoss2tgh2DNGhlHEQHcDGUEgGk2HTqp0bPXqqC0QhdENdfs0f24fRdwQ5QRAKbYnpGnkbMrl3jv37655oyhiADuijICoMHtycrXHW+tUW5xufq0a6o5d/ZjZVXAjVFGADSo/dmFuu3NNcopLFNMm0DNu7O//HnWDODWKCMAGsyhE0W6/a01ysovVZfQAP137AAF+fL0XcDdUUYANIjM3BLd9uYaHT5ZrKgWTfTO3QPUrIm32bEANAKUEQD1LqegVHfMXqP040Vq19xP7959gVoE2MyOBaCRoIwAqFd5JeUaNWet9mQVqHWQjxbcPUCtgnzMjgWgEaGMAKg3xWV23TVvnbYeyVNwE2+9c/cAhTf3MzsWgEaGMgKgXpRVOHTvOylat/+EAnw8Nf+u/urQwt/sWAAaIcoIgDpXYXdowqJULd91TL5eVs27s5+6hwWZHQtAI0UZAVCnHA5Dkz7crKTNmfK2emjWqFjFRjQ3OxaARowyAqDOGIah55K26/2UQ7J6WPTKiD4a1KmF2bEANHKUEQB1Zsb3ezV7ZZok6e839tQfYlqZnAiAM6CMAKgTC9em6x9f7ZQkPTU0WjfFtjU5EQBnQRkBcN6+3JKhKR9tliTdd2kH3T0oyuREAJxJrcrIjBkzFBkZKR8fH8XGxmrFihW/uX9paammTJmiiIgI2Ww2dejQQXPmzKlVYACNy6o92Xpo4QY5DGlE/3A9dnUXsyMBcDI1flTmokWLNGHCBM2YMUMXXnih/vOf/2jw4MHatm2b2rVrd8ZjbrnlFh09elSzZ89Wx44dlZWVpYqKivMOD8Bcmw/l6p7561Vmd+gP3Vtp+nU9ZLFYzI4FwMlYDMMwanLAgAED1LdvX82cObNqW3R0tK677jolJiaetv+XX36pW2+9Vfv27VPz5rW7vS8vL09BQUHKzc1VYGBgrb4GgLqVll2om2auUk5hmeKjgjX3zn7y8bKaHQtAI3Kun981ukxTVlamlJQUJSQkVNuekJCgVatWnfGYJUuWKC4uTi+88ILatGmjzp0769FHH1VxcfFZv09paany8vKqvQA0Hln5JRo1Z41yCssU0yZQs0bFUkQA1FqNLtNkZ2fLbrcrNDS02vbQ0FBlZmae8Zh9+/Zp5cqV8vHx0UcffaTs7Gzdd999On78+FnnjSQmJurZZ5+tSTQADaSgtEJj563TwePFatfcT3PH9FeAj5fZsQA4sVpNYP31NWHDMM56ndjhcMhisWjBggXq37+/hgwZopdeeknz5s0769mRSZMmKTc3t+p18ODB2sQEUMfKKhz68zsp2nK48sF388f2V4sAm9mxADi5Gp0ZCQkJkdVqPe0sSFZW1mlnS05p3bq12rRpo6Cg/z2XIjo6WoZh6NChQ+rUqdNpx9hsNtls/IADGhOHw9DjH2zUit3Z8vO2au6d/dQ+pInZsQC4gBqdGfH29lZsbKySk5OrbU9OTtbAgQPPeMyFF16oI0eOqKCgoGrbrl275OHhobZtWRQJcBbPf7lDH284Ik8Pi2bc3lc92zY1OxIAF1HjyzQTJ07UW2+9pTlz5mj79u16+OGHlZ6ernHjxkmqvMQyatSoqv1vu+02BQcH684779S2bdu0fPlyPfbYYxo7dqx8fX3rbiQA6s3slWmatXyfJOmFm3rq0i4tTU4EwJXUeJ2R4cOHKycnR9OmTVNGRoZiYmKUlJSkiIgISVJGRobS09Or9vf391dycrIefPBBxcXFKTg4WLfccoumT59ed6MAUG+SNmdo+ufbJElP/KGrbujLGU0AdavG64yYgXVGAHOs339ct721RmUVDo2Kj9Czf+zOomYAzlm9rDMCwH3sPVagu+evV1mFQ1dGh2rqMIoIgPpBGQFwmmP5pRozd61OFpWrV3hTvTqij6weFBEA9YMyAqCaorIK3f32/xY1mz06Tr7erK4KoP5QRgBUsTsMPbQwVRsP5aqZn5fm3dlPIf6s+QOgflFGAFSZ9ulWfb09SzZPD701Ok5RLfzNjgTADVBGAEiS5v6QprdXH5AkvTy8t2IjaveUbQCoKcoIACVvO6ppn1WuJTJpcFcN7tHa5EQA3AllBHBzmw/l6qGFqTIMaUT/dvrTxVFmRwLgZigjgBs7crJYd729TsXldg3qFKJp17KWCICGRxkB3FRBaYXGzlunrPxSdQ711+u395WXlR8JABoeP3kAN1Rhd+iBd3/Sjsx8hfjbNGdMPwX6eJkdC4CboowAbmj659v1/c5j8vGqvIW3bTM/syMBcGOUEcDNvPPjAc1btV+S9NItvdU7vKmpeQCAMgK4kZW7szV1yVZJ0iNXddYQbuEF0AhQRgA3sfdYge5bkCK7w9B1vcP0wOUdzY4EAJIoI4BbOFFYprvmrVNeSYX6tmuq52/syS28ABoNygjg4soqHPrzghTtzylSm6a++s/IOPl48RReAI0HZQRwYYZhaOqSrfpx33E18bZq9pg4tQjgKbwAGhfKCODC5q8+oIVr02WxSK+M6KOurQLNjgQAp6GMAC7qhz3ZVQ+/e/IPXXVFdKjJiQDgzCgjgAs6kFOo+xb8JLvD0PV92vDwOwCNGmUEcDH5JeW6++31yi0uV6/wpkq8oQd3zgBo1CgjgAuxOwxNeG+DdmcVKDTQplkjY7lzBkCjRxkBXMiLS3fqmx1Zsnl6aNbIOIUG+pgdCQB+F2UEcBGfbDismd/vlSS9cFNP9eKZMwCcBGUEcAFbDufqicWbJEnjLumga3u3MTkRAJw7ygjg5HIKSnXvf1NUUu7QJZ1b6LGru5gdCQBqhDICOLFyu0MPvJuqwyeL1T7YT6/c2kdWD+6cAeBcKCOAE/tb0nat3pejJt5WvTkqTkF+XmZHAoAao4wATuqDlEOa+8N+SdJLw3urU2iAuYEAoJYoI4AT2njwpCZ/tFmSNP6KTrq6eyuTEwFA7VFGACeTXVCqce+kqKzCoSujQzX+ik5mRwKA80IZAZxIhd2h+xf8pIzcEkW1aKJ/De8lDyasAnBylBHAiTz/xQ6tSTsuf5unZo2MU4APE1YBOD/KCOAkPt14RG+tTJMkvXhzL3Vs6W9yIgCoG5QRwAnszMzX4x9UrrD650s76A8xTFgF4DooI0Ajl1tcrnv/u17F5XZd1DFEjyawwioA10IZARoxh8PQxEUbtD+nSG2a+uqVEaywCsD1UEaARuy17/bomx1Zsnl66D8jY9W8ibfZkQCgzlFGgEZq+a5j+tfXuyRJ06+LUUybIJMTAUD9oIwAjdDhk8Ua/16qDEMa0b+dbo4LNzsSANQbygjQyJRW2HXfgp90oqhcPdoEaeqwbmZHAoB6RRkBGpnpn23XxoMnFeTrpRm395WPl9XsSABQrygjQCPyUeoh/ffHA7JYpJdv7a3w5n5mRwKAekcZARqJHZl5mvRh5ZN4H7ysoy7r0tLkRADQMCgjQCOQX1Ku+975SSXlDg3qFKLxV3Y2OxIANBjKCGAywzD05OLN2pddqNZBPvr3rSxsBsC9UEYAk/33xwP6fHOGPD0sev32vixsBsDtUEYAE208eFJ//WybJGnSkGj1bdfM5EQA0PAoI4BJcovKdf+7P6ncbujq7qEae2F7syMBgCkoI4AJDMPQI+9v1KETxWrX3E8v3NRLFgvzRAC4J8oIYIK3VqTp6+1H5W310Izb+yrI18vsSABgGsoI0MDW7z+u57/cIUl6elg3HoAHwO1RRoAGdKKwTA8uTJXdYeja3mG6fUA7syMBgOkoI0ADcTgq54lk5JYoKqSJnru+B/NEAECUEaDBzF6Zpm93ZMnb00Ov3dZX/jZPsyMBQKNQqzIyY8YMRUZGysfHR7GxsVqxYsU5HffDDz/I09NTvXv3rs23BZzWT+kn9Pef54lMHdZN3cICTU4EAI1HjcvIokWLNGHCBE2ZMkWpqakaNGiQBg8erPT09N88Ljc3V6NGjdIVV1xR67CAM8otKteD76aqwmFoaM/Wuq0/80QA4JdqXEZeeukl3XXXXbr77rsVHR2tl19+WeHh4Zo5c+ZvHnfvvffqtttuU3x8fK3DAs7GMAw9vnijDp+sXE8k8QbmiQDAr9WojJSVlSklJUUJCQnVtickJGjVqlVnPW7u3Lnau3evpk6dek7fp7S0VHl5edVegDN6e9V+fbX1qLysFr12Wx8F+rCeCAD8Wo3KSHZ2tux2u0JDQ6ttDw0NVWZm5hmP2b17t5588kktWLBAnp7nNmEvMTFRQUFBVa/w8PCaxAQahS2Hc/W3pMp5IpOHRKtn26bmBgKARqpWE1h/fZrZMIwznnq22+267bbb9Oyzz6pz587n/PUnTZqk3NzcqtfBgwdrExMwTUFphR549yeV2R26qluoxgxsb3YkAGi0anRvYUhIiKxW62lnQbKysk47WyJJ+fn5Wr9+vVJTU/XAAw9IkhwOhwzDkKenp5YuXarLL7/8tONsNptsNltNogGNytOfbNH+nCKFBfnoHzf1ZJ4IAPyGGp0Z8fb2VmxsrJKTk6ttT05O1sCBA0/bPzAwUJs3b9aGDRuqXuPGjVOXLl20YcMGDRgw4PzSA43QR6mH9OFPh+VhkV6+tY+a+nmbHQkAGrUar7o0ceJEjRw5UnFxcYqPj9esWbOUnp6ucePGSaq8xHL48GHNnz9fHh4eiomJqXZ8y5Yt5ePjc9p2wBWkZRfqqY+2SJIeuqKT+kc2NzkRADR+NS4jw4cPV05OjqZNm6aMjAzFxMQoKSlJERERkqSMjIzfXXMEcEVlFQ49tDBVhWV29Y9srgcv72R2JABwChbDMAyzQ/yevLw8BQUFKTc3V4GBrFyJxum5z7fpzRVpaurnpS/GD1LrIF+zIwGAqc7185tn0wB14LudWXpzRZok6YUbe1JEAKAGKCPAecrKL9Gj/7dRkjQ6PkIJ3VuZnAgAnAtlBDgPDoehR9/fpJzCMnVtFaBJQ6LNjgQATocyApyHuav2a/muY7J5eujVEX3k42U1OxIAOB3KCFBL247k6e9fVC73/tQ13dQpNMDkRADgnCgjQC0Ul9n10HupKrM7dGV0qO4Y0M7sSADgtCgjQC08l7RNe7IK1CLApr/f2IPl3gHgPFBGgBpK3nZU7/xYubDfS7f0UrA/z1ECgPNBGQFq4GheiR7/oPI23nsGRWpQpxYmJwIA50cZAc5R5W28G3WiqFzdWgfq0au7mB0JAFwCZQQ4R/NW7deK3dny8fLQKyN6y+bJbbwAUBcoI8A52JmZr+e/rLyNd8rQburYktt4AaCuUEaA31FaYdf491JVVuHQZV1acBsvANQxygjwO/65dJd2ZOYruIm3XripF7fxAkAdo4wAv2HVnmy9uWKfJOn5G3uqRQC38QJAXaOMAGeRW1SuR97fKMOQRvRvp6u6hZodCQBcEmUEOAPDMDTl483KyC1RZEgT/eUansYLAPWFMgKcwScbjuizTRmyelj08vDe8vP2NDsSALgsygjwK0dOFusvn2yRJI2/opN6hTc1NxAAuDjKCPALDoehxz7YqPySCvUOb6r7Lu1gdiQAcHmUEeAX3l69Xz/syZGPl4deuqWXPK38FQGA+sZPWuBne7Ly9fwXP6+yOiRaUS38TU4EAO6BMgJIKrc79PCijSqtcOjizi10xwURZkcCALdBGQEkvfrtHm0+nKsgXy/946aerLIKAA2IMgK3t+HgSb3+3R5J0vTrYhQa6GNyIgBwL5QRuLXiMrsmLtogu8PQH3uFaVivMLMjAYDboYzArb3w1Q7tyy5UaKBN067tbnYcAHBLlBG4rVV7szX3h/2SpL/f2FNN/bzNDQQAbooyArdUUFqhx97fJKnyIXiXdmlpciIAcF+UEbil5z7fpsMni9W2ma+mDOUheABgJsoI3M73O7O0cO1BSdI/buolfxsPwQMAM1FG4FZyi8r1xOLKyzN3Xthe8R2CTU4EAKCMwK088+lWHc0rVVRIEz1+dVez4wAARBmBG/lyS6Y+Sj0sD4v04i295OttNTsSAECUEbiJ44VleurjzZKkP13cQX3bNTM5EQDgFMoI3MLTn2xRdkGZOof66+GrOpkdBwDwC5QRuLwvNmfos00ZsnpY9OLNvWTz5PIMADQmlBG4tOOFZfrLJ1skSeMuiVLPtk3NDQQAOA1lBC7tl5dnHrqCyzMA0BhRRuCyuDwDAM6BMgKXlFNQqqc+5vIMADgDyghc0tQlW5VTWKYuoQFcngGARo4yApfz5RYuzwCAM6GMwKWcKCzTUx9vlVR5eaZH2yCTEwEAfg9lBC7lr59tU3ZBqTq25O4ZAHAWlBG4jG93HNWHPz975h839eTyDAA4CcoIXEJeSbkmf1h598xdF0WqD8+eAQCnQRmBS0hM2q7MvBK1D/bTxKu6mB0HAFADlBE4vZW7s7Vw7UFJ0t9v7Clfby7PAIAzoYzAqRWWVujJDzdJkkbFR2hAVLDJiQAANUUZgVP7x1c7dehEsdo09dXjf+hqdhwAQC1QRuC0Ug4c19ur90uSnr+xh/xtnuYGAgDUCmUETqm0wq4nFm+WYUg3x7bVoE4tzI4EAKglygic0uvf7tGerAKF+Nv01NBuZscBAJwHygiczvaMPM34fq8kadq13RXk52VyIgDA+aCMwKnYHYaeXLxJFQ5DCd1CNTimldmRAADnqVZlZMaMGYqMjJSPj49iY2O1YsWKs+774Ycf6qqrrlKLFi0UGBio+Ph4ffXVV7UODPc294c0bTyUqwAfT/31uhhZLBazIwEAzlONy8iiRYs0YcIETZkyRampqRo0aJAGDx6s9PT0M+6/fPlyXXXVVUpKSlJKSoouu+wyDRs2TKmpqecdHu4lPadILy7dKUmaMiRaoYE+JicCANQFi2EYRk0OGDBggPr27auZM2dWbYuOjtZ1112nxMTEc/oa3bt31/Dhw/X000+f0/55eXkKCgpSbm6uAgMDaxIXLsIwDI2cvVYr92QrPipY794zgLMiANDInevnd43OjJSVlSklJUUJCQnVtickJGjVqlXn9DUcDofy8/PVvHnzs+5TWlqqvLy8ai+4t/dTDmnlnmzZPD2UeEMPiggAuJAalZHs7GzZ7XaFhoZW2x4aGqrMzMxz+hr//Oc/VVhYqFtuueWs+yQmJiooKKjqFR4eXpOYcDHH8kv13OfbJUkTr+qs9iFNTE4EAKhLtZrA+ut/lRqGcU7/Ul24cKGeeeYZLVq0SC1btjzrfpMmTVJubm7V6+DBg7WJCRfx18+2Kbe4XN3DAnXXRZFmxwEA1LEarZ8dEhIiq9V62lmQrKys086W/NqiRYt011136f3339eVV175m/vabDbZbLaaRIOL+m5nlpZsPCIPi/T8DT3laeVudABwNTX6ye7t7a3Y2FglJydX256cnKyBAwee9biFCxdqzJgxevfddzV06NDaJYXbKSyt0FMfbZEkjb0wUj3aBpmcCABQH2r8ZLGJEydq5MiRiouLU3x8vGbNmqX09HSNGzdOUuUllsOHD2v+/PmSKovIqFGj9O9//1sXXHBB1VkVX19fBQXx4YKzeyl5lw6frHwi78SEzmbHAQDUkxqXkeHDhysnJ0fTpk1TRkaGYmJilJSUpIiICElSRkZGtTVH/vOf/6iiokL333+/7r///qrto0eP1rx5885/BHBJGw+e1Nwf0iRJ06+PkZ83T+QFAFdV43VGzMA6I+6l3O7QH1/7Qdsz8vTHXmF6ZUQfsyMBAGqhXtYZARrC7JVp2p6RpyBfLz09jCfyAoCro4ygUUnPKdLLX++SJE0ZGq0Qf+6qAgBXRxlBo2EYhp76ZItKyh2KjwrWzbFtzY4EAGgAlBE0Gp9uytDyXcfkbfXQc9fzRF4AcBeUETQKucXlmvbpNknS/Zd1VFQLf5MTAQAaCmUEjcILX+5QdkGpolo00bhLo8yOAwBoQJQRmC7lwAktWFO5Ns3fru8hm6fV5EQAgIZEGYGpyu0OTf5wsyTppti2uiAq2OREAICGRhmBqd5akaadR/PVzM9Lk4dEmx0HAGACyghMk55TpH9/c2pNkW5q3sTb5EQAADNQRmAKwzD0l5/XFLkgqrlu7NvG7EgAAJNQRmCKL7ZkalnVmiI9WFMEANwYZQQNLr+kXM9+ulWSNO6SKHVgTREAcGuUETS4fyXv1tG8UkUE++m+yzqaHQcAYDLKCBrUlsO5mrcqTZI07doY+XixpggAuDvKCBqMw2HoqY+3yGFIQ3u21iWdW5gdCQDQCFBG0GAWrkvXhoMn5W/z1NPXdDM7DgCgkaCMoEEcyy/V37/YIUl6JKGzQgN9TE4EAGgsKCNoEIlJ25VXUqHuYYEaeUGE2XEAAI0IZQT1bvXeHH2YelgWi/Tc9T3kaeWPHQDgf/hUQL0qq3DoL59skSTd1r+deoc3NTcQAKDRoYygXs35IU17sgoU3MRbj1/d1ew4AIBGiDKCenPkZLFe+Wa3JOnJwV0V5OdlciIAQGNEGUG9mf75NhWV2RUX0Uw39m1rdhwAQCNFGUG9WL7rmJI2Z8rqYdFfr4uRhwcPwgMAnBllBHWutMKuqUsqH4Q3Or69olsHmpwIANCYUUZQ52Yt26e07EK1DLDp4as6mR0HANDIUUZQpw4eL9Jr3+2RJE0ZGq0AHyatAgB+G2UEderZT7eqtMKh+Khg/bFXmNlxAABOgDKCOvPN9qP6enuWPD0smnZtd1ksTFoFAPw+ygjqREm5Xc9+uk2SdNdFkeoUGmByIgCAs6CMoE7MWr5P6ceLFBpo04NXMGkVAHDuKCM4bwePF+n1qkmr3eRv8zQ5EQDAmVBGcN6mf75NpRUODYhsrmE9W5sdBwDgZCgjOC/Ldh3TV1uPyuph0bRrY5i0CgCoMcoIaq20wq5nf7HSapdWTFoFANQcZQS1NntlmvZlFyrE36YJrLQKAKglyghqJSO3WK9+UzlpdfKQrgpkpVUAQC1RRlArz32+XcXldsVFNNP1fdqYHQcA4MQoI6ix1Xtz9NmmDHlYpGdZaRUAcJ4oI6iRCrtDz35aOWn1tgHt1D0syOREAABnRxlBjby7Nl07MvPV1M9Lj1zVxew4AAAXQBnBOTteWKZ/Lt0lSXokoYuaNfE2OREAwBVQRnDOXly6U7nF5YpuHajb+rczOw4AwEVQRnBOthzO1cK16ZKkZ4Z1k9WDSasAgLpBGcHvMgxDzyzZKsOQ/tgrTAOigs2OBABwIZQR/K5PNhzR+gMn5Otl1aQhXc2OAwBwMZQR/KaC0gr9LWm7JOmByzuqdZCvyYkAAK6GMoLfNOO7PcrKL1VEsJ/uuijS7DgAABdEGcFZpecU6a0VaZKkKUOi5eNlNTkRAMAVUUZwVs8lbVOZ3aFBnUJ0VbdQs+MAAFwUZQRn9MOebH219aisHhb95ZpuPH8GAFBvKCM4TYXdoWmfbpMkjbwgQp1DA0xOBABwZZQRnGbh2nTtPJqvZn5emnBlJ7PjAABcHGUE1ZwsKtM/kyufPzMxoYua+vH8GQBA/aKMoJp/Je/SyaJydW0VoBH9ws2OAwBwA7UqIzNmzFBkZKR8fHwUGxurFStW/Ob+y5YtU2xsrHx8fBQVFaU33nijVmFRv3Zm5uudNZXPn3n6mm7ytNJVAQD1r8afNosWLdKECRM0ZcoUpaamatCgQRo8eLDS09PPuH9aWpqGDBmiQYMGKTU1VZMnT9ZDDz2kxYsXn3d41A3DMLR6b44e+2Cj7A5Df+jeSgM7hpgdCwDgJiyGYRg1OWDAgAHq27evZs6cWbUtOjpa1113nRITE0/b/4knntCSJUu0ffv2qm3jxo3Txo0btXr16nP6nnl5eQoKClJubq4CAwNrEhe/ISu/RB+kHNL/rTuo/TlFkiRfL6u+mnCx2gX7mZwOAODszvXz27MmX7SsrEwpKSl68sknq21PSEjQqlWrznjM6tWrlZCQUG3b1VdfrdmzZ6u8vFxeXl6nHVNaWqrS0tJqg6kPi1MOacuR3Hr52o3dwePF+m5nluyOyi7axNuqP/ZuozsvbE8RAQA0qBqVkezsbNntdoWGVl+NMzQ0VJmZmWc8JjMz84z7V1RUKDs7W61btz7tmMTERD377LM1iVYry3Yd05KNR+r9+zRmsRHNNLxfuIb2aK0mthr9cQAAoE7U6tPn16txGobxmyt0nmn/M20/ZdKkSZo4cWLVr/Py8hQeXvd3dlzVLVThzd3zKbS+XlYldG/FgmYAANPVqIyEhITIarWedhYkKyvrtLMfp7Rq1eqM+3t6eio4OPiMx9hsNtlstppEq5VhvcI0rFdYvX8fAABwdjW6m8bb21uxsbFKTk6utj05OVkDBw484zHx8fGn7b906VLFxcWdcb4IAABwLzW+tXfixIl66623NGfOHG3fvl0PP/yw0tPTNW7cOEmVl1hGjRpVtf+4ceN04MABTZw4Udu3b9ecOXM0e/ZsPfroo3U3CgAA4LRqPGdk+PDhysnJ0bRp05SRkaGYmBglJSUpIiJCkpSRkVFtzZHIyEglJSXp4Ycf1uuvv66wsDC98soruvHGG+tuFAAAwGnVeJ0RM7DOCAAAzudcP79Z7xsAAJiKMgIAAExFGQEAAKaijAAAAFNRRgAAgKkoIwAAwFSUEQAAYCrKCAAAMBVlBAAAmKrGy8Gb4dQisXl5eSYnAQAA5+rU5/bvLfbuFGUkPz9fkhQeHm5yEgAAUFP5+fkKCgo66/tO8Wwah8OhI0eOKCAgQBaLpc6+bl5ensLDw3Xw4EGXfeaNq4+R8Tk/Vx+jq49Pcv0xMr7aMwxD+fn5CgsLk4fH2WeGOMWZEQ8PD7Vt27bevn5gYKBL/gH7JVcfI+Nzfq4+Rlcfn+T6Y2R8tfNbZ0ROYQIrAAAwFWUEAACYyq3LiM1m09SpU2Wz2cyOUm9cfYyMz/m5+hhdfXyS64+R8dU/p5jACgAAXJdbnxkBAADmo4wAAABTUUYAAICpKCMAAMBUbl1GZsyYocjISPn4+Cg2NlYrVqwwO1KtLF++XMOGDVNYWJgsFos+/vjjau8bhqFnnnlGYWFh8vX11aWXXqqtW7eaE7YWEhMT1a9fPwUEBKhly5a67rrrtHPnzmr7OPsYZ86cqZ49e1YtOhQfH68vvvii6n1nH9+vJSYmymKxaMKECVXbnHmMzzzzjCwWS7VXq1atqt535rH90uHDh3XHHXcoODhYfn5+6t27t1JSUqred+Zxtm/f/rTfQ4vFovvvv1+Sc49NkioqKvTUU08pMjJSvr6+ioqK0rRp0+RwOKr2MXWMhpt67733DC8vL+PNN980tm3bZowfP95o0qSJceDAAbOj1VhSUpIxZcoUY/HixYYk46OPPqr2/vPPP28EBAQYixcvNjZv3mwMHz7caN26tZGXl2dO4Bq6+uqrjblz5xpbtmwxNmzYYAwdOtRo166dUVBQULWPs49xyZIlxueff27s3LnT2LlzpzF58mTDy8vL2LJli2EYzj++X1q7dq3Rvn17o2fPnsb48eOrtjvzGKdOnWp0797dyMjIqHplZWVVve/MYzvl+PHjRkREhDFmzBhjzZo1RlpamvH1118be/bsqdrHmceZlZVV7fcvOTnZkGR89913hmE499gMwzCmT59uBAcHG5999pmRlpZmvP/++4a/v7/x8ssvV+1j5hjdtoz079/fGDduXLVtXbt2NZ588kmTEtWNX5cRh8NhtGrVynj++eertpWUlBhBQUHGG2+8YULC85eVlWVIMpYtW2YYhmuO0TAMo1mzZsZbb73lUuPLz883OnXqZCQnJxuXXHJJVRlx9jFOnTrV6NWr1xnfc/axnfLEE08YF1100Vnfd5VxnjJ+/HijQ4cOhsPhcImxDR061Bg7dmy1bTfccINxxx13GIZh/u+fW16mKSsrU0pKihISEqptT0hI0KpVq0xKVT/S0tKUmZlZbaw2m02XXHKJ0441NzdXktS8eXNJrjdGu92u9957T4WFhYqPj3ep8d1///0aOnSorrzyymrbXWGMu3fvVlhYmCIjI3Xrrbdq3759klxjbJK0ZMkSxcXF6eabb1bLli3Vp08fvfnmm1Xvu8o4pcrPiHfeeUdjx46VxWJxibFddNFF+uabb7Rr1y5J0saNG7Vy5UoNGTJEkvm/f07xoLy6lp2dLbvdrtDQ0GrbQ0NDlZmZaVKq+nFqPGca64EDB8yIdF4Mw9DEiRN10UUXKSYmRpLrjHHz5s2Kj49XSUmJ/P399dFHH6lbt25VPwicfXzvvfeefvrpJ61bt+6095z993DAgAGaP3++OnfurKNHj2r69OkaOHCgtm7d6vRjO2Xfvn2aOXOmJk6cqMmTJ2vt2rV66KGHZLPZNGrUKJcZpyR9/PHHOnnypMaMGSPJ+f98StITTzyh3Nxcde3aVVarVXa7Xc8995xGjBghyfwxumUZOcVisVT7tWEYp21zFa4y1gceeECbNm3SypUrT3vP2cfYpUsXbdiwQSdPntTixYs1evRoLVu2rOp9Zx7fwYMHNX78eC1dulQ+Pj5n3c9Zxzh48OCq/+/Ro4fi4+PVoUMHvf3227rgggskOe/YTnE4HIqLi9Pf/vY3SVKfPn20detWzZw5U6NGjaraz9nHKUmzZ8/W4MGDFRYWVm27M49t0aJFeuedd/Tuu++qe/fu2rBhgyZMmKCwsDCNHj26aj+zxuiWl2lCQkJktVpPOwuSlZV1Wit0dqdm9LvCWB988EEtWbJE3333ndq2bVu13VXG6O3trY4dOyouLk6JiYnq1auX/v3vf7vE+FJSUpSVlaXY2Fh5enrK09NTy5Yt0yuvvCJPT8+qcTjzGH+pSZMm6tGjh3bv3u0Sv3+S1Lp1a3Xr1q3atujoaKWnp0tynb+HBw4c0Ndff6277767apsrjO2xxx7Tk08+qVtvvVU9evTQyJEj9fDDDysxMVGS+WN0yzLi7e2t2NhYJScnV9uenJysgQMHmpSqfkRGRqpVq1bVxlpWVqZly5Y5zVgNw9ADDzygDz/8UN9++60iIyOrve8KYzwTwzBUWlrqEuO74oortHnzZm3YsKHqFRcXp9tvv10bNmxQVFSU04/xl0pLS7V9+3a1bt3aJX7/JOnCCy887Zb6Xbt2KSIiQpLr/D2cO3euWrZsqaFDh1Ztc4WxFRUVycOj+ke+1WqturXX9DHW+xTZRurUrb2zZ882tm3bZkyYMMFo0qSJsX//frOj1Vh+fr6RmppqpKamGpKMl156yUhNTa26Tfn55583goKCjA8//NDYvHmzMWLECKe6Je3Pf/6zERQUZHz//ffVbr0rKiqq2sfZxzhp0iRj+fLlRlpamrFp0yZj8uTJhoeHh7F06VLDMJx/fGfyy7tpDMO5x/jII48Y33//vbFv3z7jxx9/NK655hojICCg6ueJM4/tlLVr1xqenp7Gc889Z+zevdtYsGCB4efnZ7zzzjtV+zj7OO12u9GuXTvjiSeeOO09Zx/b6NGjjTZt2lTd2vvhhx8aISEhxuOPP161j5ljdNsyYhiG8frrrxsRERGGt7e30bdv36pbRZ3Nd999Z0g67TV69GjDMCpv2Zo6darRqlUrw2azGRdffLGxefNmc0PXwJnGJsmYO3du1T7OPsaxY8dW/Vls0aKFccUVV1QVEcNw/vGdya/LiDOP8dR6DF5eXkZYWJhxww03GFu3bq1635nH9kuffvqpERMTY9hsNqNr167GrFmzqr3v7OP86quvDEnGzp07T3vP2ceWl5dnjB8/3mjXrp3h4+NjREVFGVOmTDFKS0ur9jFzjBbDMIz6P/8CAABwZm45ZwQAADQelBEAAGAqyggAADAVZQQAAJiKMgIAAExFGQEAAKaijAAAAFNRRgAAgKkoIwAAwFSUEQAAYCrKCAAAMBVlBAAAmOr/AZHQH6PycUoTAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# adapt from fortran code\n",
    "nlev = 80\n",
    "stretch_fac = 1.0\n",
    "top_height = 22000.0\n",
    "z_flat = 16000.0\n",
    "\n",
    "vct_b = np.zeros(nlev+1)\n",
    "for k in range(1,nlev+2):\n",
    "    print('')\n",
    "    print(k)\n",
    "    x1 = (nlev + 1 - k) / nlev\n",
    "    print(x1)\n",
    "    z_height = top_height * x1 * (stretch_fac * x1 + 1.0 - stretch_fac)\n",
    "    print(z_height)\n",
    "    if (z_height >= z_flat):\n",
    "        continue\n",
    "    else:\n",
    "        vct_b[k-1] = (z_flat - z_height) / z_flat\n",
    "        print(vct_b[k-1])\n",
    "\n",
    "plt.plot(np.arange(nlev+1), vct_b)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vctb_from_vcta(vct_a, flat_height):\n",
    "    n_levels = len(vct_a)\n",
    "    # get flat_level\n",
    "    flat_levels = np.sum(vct_a > flat_height)\n",
    "    nf_levels = 0\n",
    "    vct_b = np.zeros(len(vct_a))\n",
    "    for i,a in enumerate(vct_a):\n",
    "        if a > flat_height:\n",
    "            vct_b[i] = 0.0\n",
    "        else:\n",
    "            vct_b[i] = (flat_height - a) / flat_height\n",
    "\n",
    "    return vct_b\n",
    "\n",
    "        \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x149ca1ee4340>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8RklEQVR4nO3deVxXdaL/8feX7YsgoIJsCgiuuCek4VJNC445zfSrW7aplc3NmRbN6pY5dypvDc29M03Tom1WU1k6jeW1ckqazCU1lXA3VwRkEUHhyyJf4Ps9vz9QbgSWGHC+y+v5eHwf6vmeA+/PQ4G353zO51gMwzAEAABgEh+zAwAAAO9GGQEAAKaijAAAAFNRRgAAgKkoIwAAwFSUEQAAYCrKCAAAMBVlBAAAmMrP7ADnwul0qrCwUCEhIbJYLGbHAQAA58AwDFVWVio2NlY+Pmc//+EWZaSwsFBxcXFmxwAAAOchPz9fvXv3Puv7blFGQkJCJDUOJjQ01OQ0AADgXNhsNsXFxTX9HD8btygjZy7NhIaGUkYAAHAzPzbFggmsAADAVJQRAABgKsoIAAAwFWUEAACYijICAABMRRkBAACmoowAAABTUUYAAICpKCMAAMBUbS4ja9eu1dVXX63Y2FhZLBYtX778R49Zs2aNUlJSFBgYqKSkJL300kvnkxUAAHigNpeR6upqjRgxQi+88MI57Z+Tk6OrrrpKEyZMUHZ2th599FHdd999WrZsWZvDAgAAz9PmZ9NMmjRJkyZNOuf9X3rpJcXHx+vZZ5+VJCUnJ2vr1q3605/+pOuuu66tnx4AAHiYDn9Q3saNG5Went5s28SJE7Vo0SLV19fL39+/xTF2u112u73pzzabraNjAgDg8RxOQ2VVdpVU2lVSWatjNrtKbI2/n5qWoEHR5jyMtsPLSHFxsaKioppti4qKUkNDg0pLSxUTE9PimIyMDD3xxBMdHQ0AAI/gcBoqq24sFsdsp0vG6bJx/PSvx2y1Kq2yy2m0/jFGJ/bw3DIitXx0sGEYrW4/Y+7cuZozZ07Tn202m+Li4jouIAAALsgwDNlqG3TMVqviitrTRaOxXBTbalVy+vfHq+xynK1lfI+PRQrvalVUqFWRIYGKCrWqZ0ig+keGdPBozq7Dy0h0dLSKi4ubbSspKZGfn5/Cw8NbPcZqtcpqtXZ0NAAATGMYhk5U16moolZFFbUqrjh1+tfGPx+zNf56qt5xTh/PxyJFdLUqMtSqqJBARTaVjUBFhlgbfw21Kjw4QH6+rrWyR4eXkbS0NH300UfNtq1atUqpqamtzhcBAMAT1NY7VFh+SgXlp1Rw8pQKy0+psKK28dfyxuJhb3Ce08cK6+Kv6NBARYUFKjr0TLEIbNx2+s+uWDLOVZvLSFVVlQ4ePNj055ycHG3btk09evRQfHy85s6dq4KCAr311luSpJkzZ+qFF17QnDlz9Otf/1obN27UokWL9N5777XfKAAA6GT1DqcKy08p/8Qp5Z+sUd6JGuWfqNHRk40F5Hil/cc/iKSeIVbFhDUWi5iwQMV066Lo0EBFh50pG4HqEuDbwaMxV5vLyNatW/Wzn/2s6c9n5nZMnz5db775poqKipSXl9f0fmJiolauXKn7779fL774omJjY/Xcc89xWy8AwOXV1jt0pKxaR0prlHeiWkfKapRXVqMjZdUqLD911smgZ3Tx91Wv7l3Uq1uXpl9jwgIV263x91GhgQrwc8+zGe3JYpyZTerCbDabwsLCVFFRodBQc2b6AgA8k9NpqKD8lA4er9Lh49XKKa1STmljASkoP/WDx1r9fNS7exfF9QhSXPcgxfXoorjuQerVvYt6dw9S9yD/s96s4Q3O9ed3p9xNAwCA2RxOQ7ll1dp/rFL7iqt06HiVDpZU6XBplWrrzz53IzTQT30igpUQHqw+4UGK7xGkhPBgJYQHqWdXq3x8vLdstBfKCADAoxiGoeOVdu0psmlvUeXp8lGpg8erVHeWCaMBvj7qExGkvj27KqlnsPqEByupZ7ASI7p6/dmNzkAZAQC4LafT0OHSau0urNCuggrtLarU3iKbyqrrWt0/0N9HA6JC1D8yRP2juqpfz67qG9lVcd27uO2dKJ6AMgIAcAuGYSintFrbj5Zr51GbdhVUaHdhharrWq7D4WOREiOClRwTqkHRIRoYHaqBUSHq3b0Ll1VcEGUEAOCSTlbXKTv/pLbllWvb0Qptzy9Xxan6FvsF+vtocEyohvYK0+CYUCXHhGpAVIjH3w7rSSgjAADTnTnrsTX3pLYeOaGtuSd1+Hh1i/2sfj4aEhuq4b27aVivMA3rHaakiGAusbg5yggAoNM5nYb2l1Rq06EyfZ1zQptzTrQ6zyOpZ7BGxnXTBXHdNDKuuwbFhMif4uFxKCMAgA5nGIYOHa/S+gOl2ni4TJtzTuhkTfNLLlY/H43o3U0pfborJb67UhK6q3twgEmJ0ZkoIwCADlFSWauvDpZq/YEyfXWwVMW22mbvBwX4KrVPD41J7KGLknpoWK9urEbqpSgjAIB20eBw6pu8cn25r0Sr9x3X3iJbs/etfj4andhDaX3DdVFSuIb1CuOSCyRRRgAAP0FZlV2r9x3X6n0lWrf/uGy1DU3vWSzS0NgwjesXoQn9I5SS0F2B/tzhgpYoIwCANsktq1bmnmNatfuYtuaeaPawuO5B/rpkQE9dOjBSFw/oqR7M+cA5oIwAAH6QYRjad6xSK3cU6bPdx7TvWGWz9wfHhOqK5EhdOihSI3p3ky+LiqGNKCMAgBbOFJBPdhTpk51Fzdb88PWxaExiD6UPjtIVg6PUu3uQiUnhCSgjAIAmOaXVWp5doI92FDYrIAF+PrpkQE9NGhqtywZFqlsQl1/QfigjAODlSqvs+nh7oT7cVqjt+eVN288UkMnDYnR5cqRCAv3NCwmPRhkBAC9kb3Do8z0l+kdWvtYeKJXj9CxUH4s0oX9P/WpkrK4cHEUBQaegjACAF9lbZNPft+ZreXZBsxVQR/QO069G9tIvRsQoMiTQxITwRpQRAPBwVfYGLc8u0NIt+dpZUNG0PSrUqn9L6a1rR/VW355dTUwIb0cZAQAPtbfIpnc25Wp5doGq6xySJH9fi65IjtINqXGa0D+Cp93CJVBGAMCD2Bsc+ufOYr2zKVdbc082bU/qGaybR8fr/13QS+FdrSYmBFqijACAByitsmvxpjy9vSlXpVV2SZKfj0XpQ6J065gEpfUNl8XCYmRwTZQRAHBje4tseuOrHC3fVqi6BqckKTo0UDePideNF8YpMpTJqHB9lBEAcDOGYWjdgVK9vPaQvjpY1rR9RFw3zRifqElDo3kaLtwKZQQA3ITDaeifu4q08MtD2l1ok9S4NPvPh0brjnGJGhXfjUsxcEuUEQBwcbX1Dn3wTYFeXntIuWU1kqQu/r66aXS87hjfh2fDwO1RRgDARdXWO7Rkc54WfHlIJZWNk1K7BfnrtrF9ND2tj7oH83wYeAbKCAC4mNZKSGxYoO6ckKQbR8cpKIBv3fAs/IsGABdxthJy92X9dH1KnAL8mJQKz0QZAQCTNTic+uCbAj37+X4VVtRKooTAu1BGAMAkhmHos93F+p/P9unQ8WpJUkxYoO6hhMDLUEYAwAQbDpXqj5/u0/b8ckmNE1Pv+Vk/3XpRggL9fc0NB3QyyggAdKJDx6uUsXKvPt9bIkkKCvDVneMTdefFSQoN9Dc5HWAOyggAdILymjo9+/kBvbMpVw1OQ74+Ft0yJl73XtZfPUN4cB28G2UEADpQXYNTb2/K1XP/OqCKU/WSpMsHRWruVcnqF9nV5HSAa6CMAEAHWX+gVI+t2NU0OXVQdIh+N3mwxvePMDkZ4FooIwDQzgrLT+mpT/bqk51FkqSIrgF6IH2gbkiNk68Pz44Bvo8yAgDtpK7BqUXrc/Tcvw7oVL1DPhZpWlof3X/lAIV1YXIqcDaUEQBoBxsPlWne8p06fPqSTGpCd83/1VANjg01ORng+igjAPATVNTU6w8r92rp1nxJjZdk5k5K1rWjesli4ZIMcC4oIwBwHgzD0Mc7ivTER3tUWtX4HJlbxsTrP34+iEsyQBtRRgCgjQrKT+k/l+/SF982LlzWL7KrMq4dpgv79DA5GeCeKCMAcI4Mw9CSLfl68uM9qq5zyN/Xort/1k+/ubSvrH4s4Q6cL8oIAJyDwvJTenjZDq07UCqpcYLq09cNU7/IEJOTAe6PMgIAP8AwDP0j66jmf7RHlfYGWf189NDEgbp9XCJrhgDthDICAGdRYqvV3A926l+n54ZcEN9Nf7p+hPr2ZBl3oD1RRgCgFZ/tLtYjy3boZE29Anx9NCd9gH49IYmzIUAHoIwAwHfU1DXovz7eq/c250mShsSG6i9TRmpAFHNDgI5CGQGA03YVVOi+Jdk6fLxaFov07xcn6YErByrAz8fsaIBHo4wA8HpOp6FX1x3Wn1btU73DUFSoVX+5YaTG9uPpukBnoIwA8GplVXbd//ftWrv/uCRp4pAoPX3tcHUPDjA5GeA9KCMAvNaWIyd077vZKrbVKtDfR49dPUQ3XhjHM2WATkYZAeB1nE5DL69tvCzjcBrq2zNYC25J0cBoJqkCZqCMAPAqJ6rr9MDft2n1vsbLMteMjNVT/2+Ygq18OwTMwlcfAK+x42i5Zr6dpcKKWln9fPTEL4doCpdlANNRRgB4hfe35mve8l2qa3AqMSJYC24ZpeSYULNjARBlBICHq2tw6r8+3qO3N+VKkq5IjtQzU0YqNNDf5GQAzjivlXwWLFigxMREBQYGKiUlRevWrfvB/RcvXqwRI0YoKChIMTExuv3221VWVnZegQHgXJVU1urmVzc1FZH7rxigV6amUkQAF9PmMrJ06VLNnj1b8+bNU3Z2tiZMmKBJkyYpLy+v1f3Xr1+vadOmacaMGdq9e7fef/99bdmyRXfeeedPDg8AZ/NN3kn94rn12pp7UiGBflo0PVWzrugvH54tA7icNpeRZ555RjNmzNCdd96p5ORkPfvss4qLi9PChQtb3X/Tpk3q06eP7rvvPiUmJmr8+PG66667tHXr1p8cHgBa87/bCnTjK5tUUmlX/8iuWnHPeF2eHGV2LABn0aYyUldXp6ysLKWnpzfbnp6erg0bNrR6zNixY3X06FGtXLlShmHo2LFj+sc//qHJkyef9fPY7XbZbLZmLwD4MU6noT+v2qdZS7aprsGpKwdH6cO7xykxItjsaAB+QJvKSGlpqRwOh6Kimv8PIyoqSsXFxa0eM3bsWC1evFhTpkxRQECAoqOj1a1bNz3//PNn/TwZGRkKCwtresXFxbUlJgAvdKrOoXve+0bPf3FQkjTzkr56+dYUdWX9EMDlndcE1u/fk28Yxlnv09+zZ4/uu+8+/f73v1dWVpY+/fRT5eTkaObMmWf9+HPnzlVFRUXTKz8//3xiAvASxRW1uuHljVq5s1j+vhb96foRemTSIOaHAG6iTf9liIiIkK+vb4uzICUlJS3OlpyRkZGhcePG6aGHHpIkDR8+XMHBwZowYYKefPJJxcTEtDjGarXKarW2JRoAL7W7sEJ3vLlFx2x29QgO0MtTU3Rhnx5mxwLQBm06MxIQEKCUlBRlZmY2256ZmamxY8e2ekxNTY18fJp/Gl9fX0mNZ1QA4Hyt3X9cN7y0Ucdsdg2I6qr/vXscRQRwQ22+mDpnzhxNnTpVqampSktL0yuvvKK8vLymyy5z585VQUGB3nrrLUnS1VdfrV//+tdauHChJk6cqKKiIs2ePVujR49WbGxs+44GgNd4f2u+5n6wUw1OQ2P7huulqSmsHwK4qTaXkSlTpqisrEzz589XUVGRhg4dqpUrVyohIUGSVFRU1GzNkdtuu02VlZV64YUX9MADD6hbt2667LLL9Mc//rH9RgHAaxiGoee/OKhnMvdLanzQ3X//2wgF+J3XFDgALsBiuMG1EpvNprCwMFVUVCg0lGdJAN6q3uHUfy7fpSVbGie1//bSvnpo4kAedAe4qHP9+c09bwDcwqk6h367OEur9x2Xj0V64ldDNfWiBLNjAWgHlBEALq/iVL1mvLlFW3NPKtDfR8/fNEpXDmZFVcBTUEYAuLSSylpNf32L9hbZFBrop9dvu1Cp3DEDeBTKCACXlX+iRrcu+lq5ZTWK6GrV2zNGKzmGeWOAp6GMAHBJ+49Vauqir3XMZldcjy56Z8YYJYTzjBnAE1FGALic7fnlmv7GZpXX1GtAVFe9PWOMokIDzY4FoINQRgC4lKzcE5r++hZV2Rs0Mq6b3rz9QnULCjA7FoAORBkB4DI2HS7THW9uUU2dQxcl9dCi6RcqmKfuAh6Pr3IALuGrg6Wa8bctqq13any/CL06LVVdAnzNjgWgE1BGAJjuy30luuvtLNkbnLp0YE+9dGuKAv0pIoC3oIwAMNXne47pt4u/UZ3DqSsHR+mFmy+Q1Y8iAngTyggA03y+55hmvpOlBqehq4ZF6683XiB/Xx54B3gbyggAU6zeV6LfLv5GDU5DV4+I1V9uGCE/igjglfjKB9Dp1h04rrvezlKdw6mrhkVTRAAvx1c/gE614VCp7vzbVtU1NM4R+euNF1BEAC/HdwAAnWZzzgnNeHOr7A1OXTYoUi/czBwRAJQRAJ0kK/ekbn9js07VOzShf4QW3DKKu2YASKKMAOgEewptuu2Nzaquc2hs33C9Oi2VdUQANKGMAOhQOaXVmvb616qsbVBqQne9Np0iAqA5ygiADlNUcUq3vva1SqvqlBwTqkW3XaigAFYUANAcZQRAhzhRXaepizaroPyUEiOC9dYdoxXWxd/sWABcEGUEQLurrK3XbW9s1sGSKsWEBertGaPVM8RqdiwALooyAqBd1dY79Ou3tmrH0Qr1CA7Q2zPGqHf3ILNjAXBhlBEA7cbhNDR7yTZtOnxCXa1++tvto9UvsqvZsQC4OMoIgHZhGIae+Gi3Pt1drABfH706LVXDeoeZHQuAG6CMAGgXL605rLc25kqSnpkyQml9w01OBMBdUEYA/GQfZh/VHz/9VpL0n78YrF8MjzU5EQB3QhkB8JOsO3BcD72/Q5L06wmJmjE+0eREANwNZQTAedtdWKHfvPONGpyGrh4Rq7mTks2OBMANUUYAnJfC8lO6/Y0tqrI3KC0pXH+6frh8fCxmxwLghigjANqsyt6gO97copJKuwZGhejlaSk8gRfAeaOMAGgTh9PQfe9l69viSkV0ter12y9UaCDLvAM4f5QRAG3y5Cd79MW3JbL6+ei16anq1a2L2ZEAuDnKCIBz9vbGI3rjqyOSpL9MGamRcd1MzQPAM1BGAJyTNfuP6/GP9kiSHpo4UFcNizE5EQBPQRkB8KP2H6vUPYu/kcNp6LpRvfXbS/uaHQmAB6GMAPhBJ6vrNONvW1Rpb9CYxB7KuHaYLBZu4QXQfigjAM6qweHU3e9+o/wTpxTfI0gv3ZqiAD++bQBoX3xXAXBWT63cqw2HyhQU4KtXp6Wqe3CA2ZEAeCDKCIBW/X1rftOdM8/cMFIDo0PMDQTAY1FGALTwTd5J/e7DXZKk2Vf018+HRpucCIAno4wAaOaYrVYz385SncOpiUOidN9l/c2OBMDDUUYANKmtd+jf385qeubMn28YycPvAHQ4ygiAJo+v2K3t+eXqFuSvV6elqqvVz+xIALwAZQSAJOnvW/K1ZEu+LBbp+ZsuUHx4kNmRAHgJyggA7Sqo0O/+t3HC6gNXDtCE/j1NTgTAm1BGAC9XXlOnme9kqa7BqcsHReq3l/YzOxIAL0MZAbyY02no/qXbdPRk4wqrzzBhFYAJKCOAF3v+i4Nave+4rH4+WnjrKIUF+ZsdCYAXoowAXurLfSV69l/7JUlPXjNUQ2LDTE4EwFtRRgAvVFh+SrOXbpNhSDeNjtf1qXFmRwLgxSgjgJepdzh133vZKq+p17BeYXrs6sFmRwLg5SgjgJf5S+Z+bc09qRCrn164+QIF+vuaHQmAl6OMAF5kzf7jWvDlIUnS09cNV0J4sMmJAIAyAniNY7ZazVm6TZJ060Xxmjw8xtxAAHAaZQTwAg6noVlLslVWXafkmFD9bjLzRAC4jvMqIwsWLFBiYqICAwOVkpKidevW/eD+drtd8+bNU0JCgqxWq/r27avXX3/9vAIDaLvn/nVAmw6fUHCAr15knggAF9PmR3IuXbpUs2fP1oIFCzRu3Di9/PLLmjRpkvbs2aP4+PhWj7nhhht07NgxLVq0SP369VNJSYkaGhp+cngAP27DoVI998UBSdIfrh2mpJ5dTU4EAM1ZDMMw2nLAmDFjNGrUKC1cuLBpW3Jysq655hplZGS02P/TTz/VjTfeqMOHD6tHjx7nFdJmsyksLEwVFRUKDQ09r48BeKOT1XX6+V/X6pjNrimpcfrjvw03OxIAL3KuP7/bdJmmrq5OWVlZSk9Pb7Y9PT1dGzZsaPWYFStWKDU1Vf/93/+tXr16acCAAXrwwQd16tSps34eu90um83W7AWgbQzD0CMf7NAxm11JPYP12C+ZJwLANbXpMk1paakcDoeioqKabY+KilJxcXGrxxw+fFjr169XYGCgPvzwQ5WWluq3v/2tTpw4cdZ5IxkZGXriiSfaEg3A9yzZkq/Pdh+Tv69Fz914gYIC2nxVFgA6xXlNYLVYmj/V0zCMFtvOcDqdslgsWrx4sUaPHq2rrrpKzzzzjN58882znh2ZO3euKioqml75+fnnExPwWoeOV2n+R3skSQ9NHKihvXjuDADX1ab/KkVERMjX17fFWZCSkpIWZ0vOiImJUa9evRQW9n/fDJOTk2UYho4ePar+/fu3OMZqtcpqtbYlGoDT6hqcmrUkW6fqHRrXL1x3jk8yOxIA/KA2nRkJCAhQSkqKMjMzm23PzMzU2LFjWz1m3LhxKiwsVFVVVdO2/fv3y8fHR7179z6PyAB+yJ8z92lXgU3dgvz15+tHysen9bOWAOAq2nyZZs6cOXrttdf0+uuva+/evbr//vuVl5enmTNnSmq8xDJt2rSm/W+++WaFh4fr9ttv1549e7R27Vo99NBDuuOOO9SlS5f2GwkAbThYqlfWHpYkPX3tcEWHBZqcCAB+XJtntE2ZMkVlZWWaP3++ioqKNHToUK1cuVIJCQmSpKKiIuXl5TXt37VrV2VmZuree+9VamqqwsPDdcMNN+jJJ59sv1EAUHlNneb8fbsMQ7ppdJx+PjTa7EgAcE7avM6IGVhnBPhx976XrY+2FyopIlgf3zeeu2cAmK5D1hkB4Jo+3lGoj7YXytfHor9MGUkRAeBWKCOAmyux1ep3y3dJku6+tK9GxHUzNxAAtBFlBHBjjaus7lR5Tb2GxIbqnsta3ioPAK6OMgK4saVb8vXFtyUK8PXRMzeMVIAfX9IA3A/fuQA3lX+iRv/1ceMqqw9OHKCB0SEmJwKA80MZAdyQ02nowfe3q7rOodF9emgGq6wCcGOUEcANvf5Vjr7OOaGgAF/96foR8mWVVQBujDICuJnDx6v0P5/tkyTNm5ys+PAgkxMBwE9DGQHciNNp6OFlO2RvcGpC/wjdPDre7EgA8JNRRgA38vamXG05clLBAb7KuHaYLBYuzwBwf5QRwE3kn6jRHz/9VpL0yKRB6t2dyzMAPANlBHADjYub7VBNnUNjEnvoljEJZkcCgHZDGQHcwNIt+frqYJkC/X30x+uGy4e7ZwB4EMoI4OKKKk7pqU/2SpIeTB+oPhHBJicCgPZFGQFcmGEYevSDnaq0N+iC+G66fVyi2ZEAoN1RRgAXtnxbgVbvO64AXx/9z78NZ3EzAB6JMgK4qBPVdfqvjxsvz8y6or/6RfLsGQCeiTICuKinPtmrE9V1GhQdon+/mGfPAPBclBHABW04WKpl3xyVxSL94dph8vflSxWA5+I7HOBiausdevTDnZKkqRclaFR8d5MTAUDHoowALuaFLw7qSFmNokKtemjiQLPjAECHo4wALmRfcaVeWnNIkvTEL4cqJNDf5EQA0PEoI4CLcDoNzf1ghxqchq4cHKWfD402OxIAdArKCOAiFm/O0zd55QoO8NX8Xw0xOw4AdBrKCOACSmy1+u9/Nj6R96GJAxUT1sXkRADQeSgjgAt48pO9qrQ3aETvME1N62N2HADoVJQRwGQbDpZqxfZCWSzSk9cMY8l3AF6HMgKYqK7Bqf/8312SpFvHJGhY7zCTEwFA56OMACZ6bf1hHTperYiuAXownTVFAHgnyghgkoLyU3r+XwclSXMnJSssiDVFAHgnyghgkvkf7dapeodG9+mha0f1MjsOAJiGMgKYYPW+En22+5h8fSz6r2uGymJh0ioA70UZATpZbb1Dj6/YLUm6Y1wfDYwOMTkRAJiLMgJ0spfWHFLu6QfhzbpigNlxAMB0lBGgE+WfqNHCLxsfhPefvxisrlY/kxMBgPkoI0An+sPKvbI3OJWWFK7Jw2LMjgMALoEyAnSSDYdK9c9dxfKxSI/9cjCTVgHgNMoI0AkaHE49sWKPJOnWixI0KDrU5EQA4DooI0AnWPx1nvYdq1S3IH/NuZJJqwDwXZQRoIOdrK7TM5n7JUkPpA9Ut6AAkxMBgGuhjAAd7M+Z+1Rxql6DokN08+h4s+MAgMuhjAAdaE+hTe9+nSdJevyXQ+Trw6RVAPg+ygjQQQzD0OMf7ZbTkCYPj9FFSeFmRwIAl0QZATrIyp3F2pxzQoH+Pnr0qmSz4wCAy6KMAB2gtt6hpz/dK0m66+K+6tWti8mJAMB1UUaADvC3DUeUf+KUokKtuuuSJLPjAIBLo4wA7aysyq4XvjgoSXpo4iAFBfD8GQD4IZQRoJ399V8HVGlv0NBeobr2gl5mxwEAl0cZAdrRwZJKLT59K++8qwbLh1t5AeBHUUaAdvSHld/K4TR05eAopfXlVl4AOBeUEaCdrD9Qqi++LZGfj0VzJw0yOw4AuA3KCNAOHE5DT37yf0/lTerZ1eREAOA+KCNAO/hHVr6+La5UaKCfZl3e3+w4AOBWKCPAT1Rtb9CfVjU+lfe+y/urezBP5QWAtqCMAD/RovU5Ol5pV3yPIE1L62N2HABwO5QR4Ccoq7Lr5TWHJEkPTRyoAD++pACgrc7rO+eCBQuUmJiowMBApaSkaN26ded03FdffSU/Pz+NHDnyfD4t4HKe/+KgquscGtYrTJOHxZgdBwDcUpvLyNKlSzV79mzNmzdP2dnZmjBhgiZNmqS8vLwfPK6iokLTpk3T5Zdfft5hAVeSW1atxV/nSpLmThrEAmcAcJ7aXEaeeeYZzZgxQ3feeaeSk5P17LPPKi4uTgsXLvzB4+666y7dfPPNSktLO++wgCv506r9qncYunhAT43tF2F2HABwW20qI3V1dcrKylJ6enqz7enp6dqwYcNZj3vjjTd06NAhPfbYY+f0eex2u2w2W7MX4Ep2Hq3QR9sLZbFIj/ycBc4A4KdoUxkpLS2Vw+FQVFRUs+1RUVEqLi5u9ZgDBw7okUce0eLFi+Xnd25PL83IyFBYWFjTKy4uri0xgQ73x0+/lSRdM7KXBseGmpwGANzbeU1gtViaXxs3DKPFNklyOBy6+eab9cQTT2jAgAHn/PHnzp2rioqKpld+fv75xAQ6xNr9x7X+YKkCfH0058pz/3cNAGjduZ2qOC0iIkK+vr4tzoKUlJS0OFsiSZWVldq6dauys7N1zz33SJKcTqcMw5Cfn59WrVqlyy67rMVxVqtVVqu1LdGATuF0Gnr6n41nRaamJSiuR5DJiQDA/bXpzEhAQIBSUlKUmZnZbHtmZqbGjh3bYv/Q0FDt3LlT27Zta3rNnDlTAwcO1LZt2zRmzJiflh7oZCu2F2pPkU0hVj/d/bN+ZscBAI/QpjMjkjRnzhxNnTpVqampSktL0yuvvKK8vDzNnDlTUuMlloKCAr311lvy8fHR0KFDmx0fGRmpwMDAFtsBV1fX4NSfM/dJkmZe2lc9WPYdANpFm8vIlClTVFZWpvnz56uoqEhDhw7VypUrlZCQIEkqKir60TVHAHf09635yj9xShFdrbp9XB+z4wCAx7AYhmGYHeLH2Gw2hYWFqaKiQqGh3LmAzldb79Al/7Nax2x2PX71YN02LtHsSADg8s715zcP0gDOwdsbc3XMZlevbl1005h4s+MAgEehjAA/orK2Xgu+PChJmnV5f1n9fE1OBACehTIC/IjX1x/RyZp6JUUE69pRvcyOAwAehzIC/IDymjq9tu6wJOn+KwfIz5cvGQBob3xnBX7AS2sOq9LeoOSYUE0eFmN2HADwSJQR4CxKbLV6c0OOJOnB9AHy8Wn5yAMAwE9HGQHO4sXVB1Vb79QF8d102aBIs+MAgMeijACtOHqyRu9ubly876H0ga0+CBIA0D4oI0ArXvjioOodhsb2DdfYfhFmxwEAj0YZAb4n/0SN/pF1VJI058oBJqcBAM9HGQG+54UvDqrBaWhC/wil9ulhdhwA8HiUEeA78spqtOybxrMis6/ob3IaAPAOlBHgO15YfaDprEhKAmdFAKAzUEaA03LLqrXsmwJJ0uwrmCsCAJ2FMgKc9sIXB+VwGrp4QE+lJHQ3Ow4AeA3KCCDpSGm1Psg+c1aEuSIA0JkoI4CkF1Y3nhW5ZEBPjYrnrAgAdCbKCLzekdJqfchZEQAwDWUEXu/503NFLh3YUxdwVgQAOh1lBF4tt6xay7dxBw0AmIkyAq+28MtDTXNFRsZ1MzsOAHglygi8VkH5qabVVu+7vJ/JaQDAe1FG4LVeWXNI9Q5DaUnhrLYKACaijMArlVTW6r0t+ZKkey/jrAgAmIkyAq/02roc1TU4NSq+m9L6hpsdBwC8GmUEXudEdZ3e2ZQrSbr3sv6yWCwmJwIA70YZgdd546sc1dQ5NCQ2VJcO7Gl2HADwepQReJWKU/V686sjkhrninBWBADMRxmBV3l74xFV2hs0IKqr0gdHmx0HACDKCLxItb1Bi9bnSJLu/lk/+fhwVgQAXAFlBF7j3a/zdLKmXokRwfrF8Fiz4wAATqOMwCvYGxx6dd1hSdJvLukrX86KAIDLoIzAKyzPLlBJpV3RoYG65oJeZscBAHwHZQQez+E09PKaxrMid05IVIAf/+wBwJXwXRkeb9XuYh0urVZYF3/dNDre7DgAgO+hjMCjGYahhWsOSZKmpyUo2OpnciIAwPdRRuDRNhwq046jFQr099H0sX3MjgMAaAVlBB7tpdNnRaakxim8q9XkNACA1lBG4LF2Hq3QugOl8vWx6M4JSWbHAQCcBWUEHuvMWZFfjohVXI8gk9MAAM6GMgKPlFNarZW7iiRJd13CWREAcGWUEXikV9YekmFIlw2K1KDoULPjAAB+AGUEHqfEVqtlWQWSpN9c2tfkNACAH0MZgcd5c8MR1TmcSknorgv79DA7DgDgR1BG4FGq7Q16Z1OuJOnfL2auCAC4A8oIPMr7W/Nlq21QYkSwrkiOMjsOAOAcUEbgMRocTi36KkeSdMf4RPn6WExOBAA4F5QReIzPdh9T/olT6h7kr38b1dvsOACAc0QZgUcwDEOvrDssSZqa1kddAnxNTgQAOFeUEXiErbkntT2/XAF+PpqWlmB2HABAG1BG4BFeWdt4VuS6Ub0UwQPxAMCtUEbg9g4fr9Lne49JkmaM53ZeAHA3lBG4vUXrc2QY0hXJkeoX2dXsOACANqKMwK2VVdn1j6yjkqQ7J3BWBADcEWUEbu2dTXmyNzg1vHeYxiSy9DsAuCPKCNxWbb1Db286Ikn69YQkWSwscgYA7ui8ysiCBQuUmJiowMBApaSkaN26dWfd94MPPtCVV16pnj17KjQ0VGlpafrss8/OOzBwxkfbC1VaVafYsEBNGhptdhwAwHlqcxlZunSpZs+erXnz5ik7O1sTJkzQpEmTlJeX1+r+a9eu1ZVXXqmVK1cqKytLP/vZz3T11VcrOzv7J4eH9zIMQ69/dUSSNH1sH/n5cpIPANyVxTAMoy0HjBkzRqNGjdLChQubtiUnJ+uaa65RRkbGOX2MIUOGaMqUKfr9739/TvvbbDaFhYWpoqJCoaGhbYkLD7XxUJluenWTuvj7atPcyxUW5G92JADA95zrz+82/Xeyrq5OWVlZSk9Pb7Y9PT1dGzZsOKeP4XQ6VVlZqR49zj7Z0G63y2azNXsB3/X66QfiXZfSiyICAG6uTWWktLRUDodDUVHNH80eFRWl4uLic/oYf/7zn1VdXa0bbrjhrPtkZGQoLCys6RUXF9eWmPBweWU1TYuc3TY20eQ0AICf6rwutH//rgXDMM7pTob33ntPjz/+uJYuXarIyMiz7jd37lxVVFQ0vfLz888nJjzUmxuOyDCkSwb0ZJEzAPAAfm3ZOSIiQr6+vi3OgpSUlLQ4W/J9S5cu1YwZM/T+++/riiuu+MF9rVarrFaeL4KWKmvr9fetjeX0jvGcFQEAT9CmMyMBAQFKSUlRZmZms+2ZmZkaO3bsWY977733dNttt+ndd9/V5MmTzy8pIOn9rUdVZW9Qv8iuurh/hNlxAADtoE1nRiRpzpw5mjp1qlJTU5WWlqZXXnlFeXl5mjlzpqTGSywFBQV66623JDUWkWnTpumvf/2rLrrooqazKl26dFFYWFg7DgWezuE09LeNRyRJt4/rwyJnAOAh2lxGpkyZorKyMs2fP19FRUUaOnSoVq5cqYSEBElSUVFRszVHXn75ZTU0NOjuu+/W3Xff3bR9+vTpevPNN3/6COA1vvi2RLllNQrr4q9rL+htdhwAQDtp8zojZmCdEUjSTa9s0sbDZZp5SV89MmmQ2XEAAD+iQ9YZAcyyt8imjYfL5Otj0bS0BLPjAADaEWUEbuHN00u//3xotGK7dTE3DACgXVFG4PLKa+q0fFuBJOn2sX3MDQMAaHeUEbi8v2/Nl73BqSGxoUpJ6G52HABAO6OMwKU5nIbe3pQrSZqexu28AOCJKCNwaV/uK1H+iVPqFuSvX46MNTsOAKADUEbg0v62sfGsyJTUOAX6+5qcBgDQESgjcFmHj1dp7f7jslikWy/idl4A8FSUEbist06fFbl8UKTiegSZnAYA0FEoI3BJVfYGLcs6KkmaltbH3DAAgA5FGYFL+jC7QJX2BiVFBGt8P57OCwCejDICl2MYht7acESSNDUtQT4+3M4LAJ6MMgKXs/FwmQ6UVCkowFfXpfB0XgDwdJQRuJy3NjROXL1uVG+FBvqbnAYA0NEoI3ApheWntGpPsSTxdF4A8BKUEbiUJZvz5DSki5J6qH9UiNlxAACdgDICl1HvcGrJlnxJ0tSL+pgbBgDQaSgjcBmf7zmmkkq7IrpadeXgKLPjAAA6CWUELmPx13mSpBsvjFOAH/80AcBb8B0fLiGntFrrD5bKYpFuHB1ndhwAQCeijMAlLN7UeDvvZQMj1bs7z6EBAG9CGYHpausd+sc3jc+hueWieJPTAAA6G2UEpvtkR5HKa+rVq1sXXTIg0uw4AIBORhmB6RZ/3XiJ5uYx8fLlOTQA4HUoIzDVnkKbvskrl5+PRTekMnEVALwRZQSmOnNWZOLQaPUMsZqcBgBgBsoITFNlb9Dy7AJJ0q1jeA4NAHgryghMszy7QNV1DvXtGayLknqYHQcAYBLKCExhGEbTiqu3jEmQxcLEVQDwVpQRmGL70QrtLbLJ6uej60b1NjsOAMBElBGYYsnmxrMiVw2LUViQv8lpAABmooyg01XZG7Rie6GkxofiAQC8G2UEne6j7YWqqXMoqWewRicycRUAvB1lBJ3uzCWamy6MZ+IqAIAygs61u7BC249WyN/XomtH9TI7DgDABVBG0KmWbM6XJKUPiVZ4V1ZcBQBQRtCJTtU5tHxb44qrN10Yb3IaAICroIyg03yys0iVtQ2K69FFY/uGmx0HAOAiKCPoNGcmrt54Ybx8fJi4CgBoRBlBpzhwrFJbc0/K18ei61NYcRUA8H8oI+gU752euHrZoEhFhgaanAYA4EooI+hwtfUOfZB9VJJ002hWXAUANEcZQYf7bHexymvqFRMWqEsGRJodBwDgYigj6HDvb208K3J9apx8mbgKAPgeygg6VP6JGn11qFSSmLgKAGgVZQQdatk3R2UY0rh+4YrrEWR2HACAC6KMoMM4nUbTJZobUpm4CgBoHWUEHWbj4TIVlJ9SSKCfJg6JNjsOAMBFUUbQYf6+tXFtkV+NjFWgv6/JaQAArooygg5RUVOvf+4qlsQlGgDAD6OMoEOs2FGouganBkWHaFivMLPjAABcGGUEHeL905dork+Nk8XC2iIAgLOjjKDd7S2yacfRCvn7WnTNyFiz4wAAXBxlBO3uzO28VyRHKbyr1eQ0AABXRxlBu6prcGr5tgJJTFwFAJwbygja1RffHtOJ6jpFhVo1oX+E2XEAAG7gvMrIggULlJiYqMDAQKWkpGjdunU/uP+aNWuUkpKiwMBAJSUl6aWXXjqvsHB9S7c0Tly9blRv+fnSdQEAP67NPy2WLl2q2bNna968ecrOztaECRM0adIk5eXltbp/Tk6OrrrqKk2YMEHZ2dl69NFHdd9992nZsmU/OTxcx+acE5q66Gut3ndcUuNdNAAAnAuLYRhGWw4YM2aMRo0apYULFzZtS05O1jXXXKOMjIwW+z/88MNasWKF9u7d27Rt5syZ2r59uzZu3HhOn9NmsyksLEwVFRUKDQ1tS1x0IMMw9NXBMj33xQFtzjkhSfL1sei3l/bVA+kDTU4HADDbuf789mvLB62rq1NWVpYeeeSRZtvT09O1YcOGVo/ZuHGj0tPTm22bOHGiFi1apPr6evn7+7c4xm63y263NxtMR1iWdVS7Cis65GN7g2355crOK5ck+ftadH1qnH5zSV+ezgsAaJM2lZHS0lI5HA5FRUU12x4VFaXi4uJWjykuLm51/4aGBpWWliomJqbFMRkZGXriiSfaEu28rNl/XCu2F3b45/FkVj8f3TQ6XnddkqSYsC5mxwEAuKE2lZEzvr+ipmEYP7jKZmv7t7b9jLlz52rOnDlNf7bZbIqLa/85CFcOjlJcD36Anq/QQH9dO6q3eoawlggA4Py1qYxERETI19e3xVmQkpKSFmc/zoiOjm51fz8/P4WHh7d6jNVqldXa8T/grh4Rq6tHsEIoAABmatPdNAEBAUpJSVFmZmaz7ZmZmRo7dmyrx6SlpbXYf9WqVUpNTW11vggAAPAubb61d86cOXrttdf0+uuva+/evbr//vuVl5enmTNnSmq8xDJt2rSm/WfOnKnc3FzNmTNHe/fu1euvv65FixbpwQcfbL9RAAAAt9XmOSNTpkxRWVmZ5s+fr6KiIg0dOlQrV65UQkKCJKmoqKjZmiOJiYlauXKl7r//fr344ouKjY3Vc889p+uuu679RgEAANxWm9cZMQPrjAAA4H7O9ec363UDAABTUUYAAICpKCMAAMBUlBEAAGAqyggAADAVZQQAAJiKMgIAAExFGQEAAKaijAAAAFO1eTl4M5xZJNZms5mcBAAAnKszP7d/bLF3tygjlZWVkqS4uDiTkwAAgLaqrKxUWFjYWd93i2fTOJ1OFRYWKiQkRBaLpd0+rs1mU1xcnPLz8z32mTeePkbG5/48fYyePj7J88fI+M6fYRiqrKxUbGysfHzOPjPELc6M+Pj4qHfv3h328UNDQz3yH9h3efoYGZ/78/Qxevr4JM8fI+M7Pz90RuQMJrACAABTUUYAAICpvLqMWK1WPfbYY7JarWZH6TCePkbG5/48fYyePj7J88fI+DqeW0xgBQAAnsurz4wAAADzUUYAAICpKCMAAMBUlBEAAGAqry4jCxYsUGJiogIDA5WSkqJ169aZHem8rF27VldffbViY2NlsVi0fPnyZu8bhqHHH39csbGx6tKliy699FLt3r3bnLDnISMjQxdeeKFCQkIUGRmpa665Rvv27Wu2j7uPceHChRo+fHjTokNpaWn65z//2fS+u4/v+zIyMmSxWDR79uymbe48xscff1wWi6XZKzo6uul9dx7bdxUUFOjWW29VeHi4goKCNHLkSGVlZTW9787j7NOnT4u/Q4vForvvvluSe49NkhoaGvS73/1OiYmJ6tKli5KSkjR//nw5nc6mfUwdo+GllixZYvj7+xuvvvqqsWfPHmPWrFlGcHCwkZuba3a0Nlu5cqUxb948Y9myZYYk48MPP2z2/tNPP22EhIQYy5YtM3bu3GlMmTLFiImJMWw2mzmB22jixInGG2+8YezatcvYtm2bMXnyZCM+Pt6oqqpq2sfdx7hixQrjk08+Mfbt22fs27fPePTRRw1/f39j165dhmG4//i+a/PmzUafPn2M4cOHG7NmzWra7s5jfOyxx4whQ4YYRUVFTa+SkpKm9915bGecOHHCSEhIMG677Tbj66+/NnJycozPP//cOHjwYNM+7jzOkpKSZn9/mZmZhiRj9erVhmG499gMwzCefPJJIzw83Pj444+NnJwc4/333ze6du1qPPvss037mDlGry0jo0ePNmbOnNls26BBg4xHHnnEpETt4/tlxOl0GtHR0cbTTz/dtK22ttYICwszXnrpJRMS/nQlJSWGJGPNmjWGYXjmGA3DMLp372689tprHjW+yspKo3///kZmZqZxySWXNJURdx/jY489ZowYMaLV99x9bGc8/PDDxvjx48/6vqeM84xZs2YZffv2NZxOp0eMbfLkycYdd9zRbNu1115r3HrrrYZhmP/355WXaerq6pSVlaX09PRm29PT07VhwwaTUnWMnJwcFRcXNxur1WrVJZdc4rZjraiokCT16NFDkueN0eFwaMmSJaqurlZaWppHje/uu+/W5MmTdcUVVzTb7gljPHDggGJjY5WYmKgbb7xRhw8fluQZY5OkFStWKDU1Vddff70iIyN1wQUX6NVXX21631PGKTX+jHjnnXd0xx13yGKxeMTYxo8fr3/961/av3+/JGn79u1av369rrrqKknm//25xYPy2ltpaakcDoeioqKabY+KilJxcbFJqTrGmfG0Ntbc3FwzIv0khmFozpw5Gj9+vIYOHSrJc8a4c+dOpaWlqba2Vl27dtWHH36owYMHN30jcPfxLVmyRN988422bNnS4j13/zscM2aM3nrrLQ0YMEDHjh3Tk08+qbFjx2r37t1uP7YzDh8+rIULF2rOnDl69NFHtXnzZt13332yWq2aNm2ax4xTkpYvX67y8nLddtttktz/36ckPfzww6qoqNCgQYPk6+srh8Ohp556SjfddJMk88folWXkDIvF0uzPhmG02OYpPGWs99xzj3bs2KH169e3eM/dxzhw4EBt27ZN5eXlWrZsmaZPn641a9Y0ve/O48vPz9esWbO0atUqBQYGnnU/dx3jpEmTmn4/bNgwpaWlqW/fvvrb3/6miy66SJL7ju0Mp9Op1NRU/eEPf5AkXXDBBdq9e7cWLlyoadOmNe3n7uOUpEWLFmnSpEmKjY1ttt2dx7Z06VK98847evfddzVkyBBt27ZNs2fPVmxsrKZPn960n1lj9MrLNBEREfL19W1xFqSkpKRFK3R3Z2b0e8JY7733Xq1YsUKrV69W7969m7Z7yhgDAgLUr18/paamKiMjQyNGjNBf//pXjxhfVlaWSkpKlJKSIj8/P/n5+WnNmjV67rnn5Ofn1zQOdx7jdwUHB2vYsGE6cOCAR/z9SVJMTIwGDx7cbFtycrLy8vIkec7XYW5urj7//HPdeeedTds8YWwPPfSQHnnkEd14440aNmyYpk6dqvvvv18ZGRmSzB+jV5aRgIAApaSkKDMzs9n2zMxMjR071qRUHSMxMVHR0dHNxlpXV6c1a9a4zVgNw9A999yjDz74QF988YUSExObve8JY2yNYRiy2+0eMb7LL79cO3fu1LZt25peqampuuWWW7Rt2zYlJSW5/Ri/y263a+/evYqJifGIvz9JGjduXItb6vfv36+EhARJnvN1+MYbbygyMlKTJ09u2uYJY6upqZGPT/Mf+b6+vk239po+xg6fIuuiztzau2jRImPPnj3G7NmzjeDgYOPIkSNmR2uzyspKIzs728jOzjYkGc8884yRnZ3ddJvy008/bYSFhRkffPCBsXPnTuOmm25yq1vSfvOb3xhhYWHGl19+2ezWu5qamqZ93H2Mc+fONdauXWvk5OQYO3bsMB599FHDx8fHWLVqlWEY7j++1nz3bhrDcO8xPvDAA8aXX35pHD582Ni0aZPxi1/8wggJCWn6fuLOYztj8+bNhp+fn/HUU08ZBw4cMBYvXmwEBQUZ77zzTtM+7j5Oh8NhxMfHGw8//HCL99x9bNOnTzd69erVdGvvBx98YERERBj/8R//0bSPmWP02jJiGIbx4osvGgkJCUZAQIAxatSopltF3c3q1asNSS1e06dPNwyj8Zatxx57zIiOjjasVqtx8cUXGzt37jQ3dBu0NjZJxhtvvNG0j7uP8Y477mj6t9izZ0/j8ssvbyoihuH+42vN98uIO4/xzHoM/v7+RmxsrHHttdcau3fvbnrfncf2XR999JExdOhQw2q1GoMGDTJeeeWVZu+7+zg/++wzQ5Kxb9++Fu+5+9hsNpsxa9YsIz4+3ggMDDSSkpKMefPmGXa7vWkfM8doMQzD6PjzLwAAAK3zyjkjAADAdVBGAACAqSgjAADAVJQRAABgKsoIAAAwFWUEAACYijICAABMRRkBAACmoowAAABTUUYAAICpKCMAAMBUlBEAAGCq/w8wFGTn1YkvUgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.arange(81), vctb_from_vcta(vct['vct_a'].values, flat_height=12000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 1   22000.000   0.00000\n",
      " 2   21187.520   0.00000\n",
      " 3   20396.250   0.00000\n",
      " 4   19625.890   0.00000\n",
      " 5   18876.130   0.00000\n",
      " 6   18146.680   0.00000\n",
      " 7   17437.220   0.00000\n",
      " 8   16747.450   0.00000\n",
      " 9   16077.080   0.00000\n",
      "10   15425.800   0.00000\n",
      "11   14793.300   0.00000\n",
      "12   14179.290   0.00000\n",
      "13   13583.470   0.00000\n",
      "14   13005.530   0.00000\n",
      "15   12445.190   0.00000\n",
      "16   11902.130   0.00816\n",
      "17   11376.060   0.05200\n",
      "18   10866.690   0.09444\n",
      "19   10373.710   0.13552\n",
      "20    9896.840   0.17526\n",
      "21    9435.770   0.21369\n",
      "22    8990.210   0.25082\n",
      "23    8559.870   0.28668\n",
      "24    8144.450   0.32130\n",
      "25    7743.660   0.35470\n",
      "26    7357.210   0.38690\n",
      "27    6984.800   0.41793\n",
      "28    6626.150   0.44782\n",
      "29    6280.970   0.47659\n",
      "30    5948.960   0.50425\n",
      "31    5629.830   0.53085\n",
      "32    5323.310   0.55639\n",
      "33    5029.090   0.58091\n",
      "34    4746.910   0.60442\n",
      "35    4476.460   0.62696\n",
      "36    4217.460   0.64855\n",
      "37    3969.640   0.66920\n",
      "38    3732.700   0.68894\n",
      "39    3506.370   0.70780\n",
      "40    3290.360   0.72580\n",
      "41    3084.400   0.74297\n",
      "42    2888.200   0.75932\n",
      "43    2701.480   0.77488\n",
      "44    2523.960   0.78967\n",
      "45    2355.380   0.80372\n",
      "46    2195.450   0.81705\n",
      "47    2043.900   0.82968\n",
      "48    1900.450   0.84163\n",
      "49    1764.840   0.85293\n",
      "50    1636.780   0.86360\n",
      "51    1516.020   0.87366\n",
      "52    1402.270   0.88314\n",
      "53    1295.280   0.89206\n",
      "54    1194.770   0.90044\n",
      "55    1100.480   0.90829\n",
      "56    1012.140   0.91565\n",
      "57     929.490   0.92254\n",
      "58     852.280   0.92898\n",
      "59     780.230   0.93498\n",
      "60     713.090   0.94058\n",
      "61     650.610   0.94578\n",
      "62     592.520   0.95062\n",
      "63     538.570   0.95512\n",
      "64     488.520   0.95929\n",
      "65     442.110   0.96316\n",
      "66     399.090   0.96674\n",
      "67     359.210   0.97007\n",
      "68     322.240   0.97315\n",
      "69     287.930   0.97601\n",
      "70     256.040   0.97866\n",
      "71     226.330   0.98114\n",
      "72     198.570   0.98345\n",
      "73     172.540   0.98562\n",
      "74     148.000   0.98767\n",
      "75     124.730   0.98961\n",
      "76     102.510   0.99146\n",
      "77      81.130   0.99324\n",
      "78      60.390   0.99497\n",
      "79      40.070   0.99666\n",
      "80      20.000   0.99833\n",
      "81       0.000   1.00000\n"
     ]
    }
   ],
   "source": [
    "vct_a = vct['vct_a'].values\n",
    "vct_b = vctb_from_vcta(vct_a, 12000)\n",
    "for i in range(81):\n",
    "    print(f'{(i+1):2}  {vct_a[i]:10.3f}   {vct_b[i]:6.5f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.6 ('vert_coord')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "db3f4cc0794993cb75ebd5fb7aef00324c0169b209758b1df38822e455387f35"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}