plot pb_soil vs alpha

alpha_vs_pbsoil.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "4806a516-bd2f-478f-85a0-ddea17f71b11",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAHFCAYAAAAe+pb9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJ9ElEQVR4nO3deVhV5f7//9dmHgRkFhQR59kUzDktlVLzZJ8GG45DaWVWptY5abPV79h0yjqlZZOnNPPbUcvKShqcUnMIpzTTVEBlEEFmGdfvD2InggoIbNjr+biufcm+973Wet/7tni5RothGIYAAABMyMHWBQAAANgKQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQimcv3118vd3V2nT58+b5/bb79dzs7OSk5OvuTtHT16VBaLRYsWLar2sidOnNDTTz+tnTt3Vvjs6aeflsViueT66lJBQYGmTJmikJAQOTo66rLLLrvkdS5atEgWi0VHjx6t9rJlc/Hyyy9fch11obKxffzxx5o3b16FvvU5lokTJ8pisVhfjo6OatGihW6++Wbt3bu3zrZb2fcxceJEtWrVqs62CXMiCMFUJk2apDNnzujjjz+u9POMjAytXLlS1157rYKDgy95eyEhIdq8ebNGjRpV7WVPnDihOXPmVBqEJk+erM2bN19yfXVpwYIFevvtt/XYY49p48aN+uijj2xdUoM2atQobd68WSEhIda28wWh+ubu7q7Nmzdr8+bNWrdunZ577jn98ssv6t+/v44fP14n26zs+wDqgpOtCwDq04gRIxQaGqr3339fU6dOrfD50qVLlZeXp0mTJl3SdoqLi1VUVCRXV1f17dv3ktZVmRYtWqhFixa1vt7atHfvXrm7u+v++++3dSmNQmBgoAIDA21dRqUcHBzK/T0eOHCgWrZsqaFDh+qrr77S3XffXevbbMjfB+wLe4RgKo6OjpowYYJ27NihPXv2VPj8gw8+UEhIiEaMGKGTJ09q6tSp6ty5s5o0aaKgoCBdddVV2rBhQ7llyg5TvPjii3ruuecUEREhV1dX/fjjj5UeGjt06JDuuOMOtWvXTh4eHmrevLlGjx5drp61a9eqd+/ekqQ77rjDelji6aefllT5obGSkhK9+OKL6tixo1xdXRUUFKTx48fr2LFj5foNGTJEXbt21bZt2zRo0CB5eHiodevWev7551VSUnLR7/DMmTOaPXu2IiIi5OLioubNm+u+++4rd7jRYrHo3XffVV5enrX2Cx0ejImJ0XXXXacWLVrIzc1Nbdu21T333KPU1NSL1lM2ng0bNqhv375yd3dX8+bN9cQTT6i4uLjSZV555RVFRESoSZMm6tevn7Zs2VLu8+3bt+uWW25Rq1at5O7urlatWunWW29VXFzcRevp3bt3hT2A3bp1k8Vi0bZt26xtK1askMVisc77uYeChgwZoq+++kpxcXHlDk1VdyyVKdtWTEyM7rjjDvn5+cnT01OjR4/W4cOHL7q8JPn4+EiSnJ2drW25ubl6+OGHFRERITc3N/n5+SkqKkpLly4tt+yqVavUr18/eXh4yMvLS8OHD6+wh/NSDoMC1UEQgunceeedslgsev/998u179u3T1u3btWECRPk6OiotLQ0SdJTTz2lr776Sh988IFat26tIUOGaO3atRXW+/rrr+uHH37Qyy+/rK+//lodO3asdPsnTpyQv7+/nn/+eX3zzTd688035eTkpD59+ujAgQOSpF69eumDDz6QJD3++OPWwxKTJ08+77juvfdePfLIIxo+fLhWrVqlZ599Vt9884369+9fIVAkJSXp9ttv19///netWrVKI0aM0OzZs7V48eILfneGYWjMmDF6+eWXNW7cOH311VeaOXOm/vvf/+qqq65Sfn6+JGnz5s0aOXJkuUMqFzo8+Mcff6hfv35asGCB1qxZoyeffFI///yzBg4cqMLCwgvWVDaeW265Rbfffrs+//xz3XjjjXruuef04IMPVuj75ptvKiYmRvPmzdOSJUuUk5OjkSNHKiMjw9rn6NGj6tChg+bNm6dvv/1WL7zwghITE9W7d++LhrNhw4Zp/fr11rqTk5Ote8diYmKs/b777jsFBwerW7dula5n/vz5GjBggJo1a2b9Ds8NC1UZy4VMmjRJDg4O1kNwW7du1ZAhQyo9h66oqEhFRUU6c+aM9u7dq3/84x/y9fUtN68zZ87UggULNG3aNH3zzTf66KOPdNNNN+nUqVPWPh9//LGuu+46eXt7a+nSpXrvvfeUnp6uIUOGaOPGjVWqG6hVBmBCgwcPNgICAoyCggJr20MPPWRIMn7//fdKlykqKjIKCwuNoUOHGtdff721/ciRI4Yko02bNuXWd/ZnH3zwwXlrKSoqMgoKCox27doZM2bMsLZv27btvMs+9dRTxtn/+e7fv9+QZEydOrVcv59//tmQZDz66KPlxi7J+Pnnn8v17dy5s3H11Veft07DMIxvvvnGkGS8+OKL5dqXLVtmSDIWLlxobZswYYLh6el5wfVVpqSkxCgsLDTi4uIMScbnn39u/eyDDz4wJBlHjhypMJ6z+xmGYdx1112Gg4ODERcXZxjGX3PRrVs3o6ioyNpv69athiRj6dKl562pqKjIyM7ONjw9PY3XXnvtgvV/9913hiRj/fr1hmEYxuLFiw0vLy9j6tSpxpVXXmnt165dO+O222674NhGjRplhIeHV9jGpYzl7G2d/ffYMAzjp59+MiQZzz33nLVtwoQJhqQKr5CQEGPjxo3llu/atasxZsyY8263uLjYCA0NNbp162YUFxdb27OysoygoCCjf//+F/w+JkyYUOn3AVwK9gjBlCZNmqTU1FStWrVKUum/dhcvXqxBgwapXbt21n5vvfWWevXqJTc3Nzk5OcnZ2Vnff/+99u/fX2Gdf/vb38odJjifoqIi/etf/1Lnzp3l4uIiJycnubi46ODBg5Wutyp+/PFHSaVX1Zzt8ssvV6dOnfT999+Xa2/WrJkuv/zycm3du3e/6KGfH374odLt3HTTTfL09KywnapKSUnRlClTFBYWZv2ew8PDJalK34mXl5f+9re/lWu77bbbVFJSovXr15drHzVqlBwdHa3vu3fvLknlxp6dna1HHnlEbdu2lZOTk5ycnNSkSRPl5ORctJ4BAwbIzc1N3333naTSw35DhgzRNddco02bNik3N1cJCQk6ePCghg0bdtGxXUhVxnIht99+e7n3/fv3V3h4uPXvUxl3d3dt27ZN27Zt088//6wVK1aoffv2GjlyZLm9VJdffrm+/vprzZo1S2vXrlVeXl659Rw4cEAnTpzQuHHj5ODw16+fJk2a6IYbbtCWLVuUm5tbtcEDtYQgBFO68cYb5ePjYz38tHr1aiUnJ5c7SfqVV17Rvffeqz59+mj58uXasmWLtm3bpmuuuabC/+AlVfnqlpkzZ+qJJ57QmDFj9MUXX+jnn3/Wtm3b1KNHj0rXWxVlhx4qqyE0NLTcoQlJ8vf3r9DP1dX1ots/deqUnJycKpzEarFY1KxZswrbqYqSkhJFR0drxYoV+uc//6nvv/9eW7dutZ7rUpXvpLIr/Jo1a2at+Wznjt3V1bXCdm677Ta98cYbmjx5sr799ltt3bpV27ZtU2Bg4EXrcXNz04ABA6xB6Pvvv9fw4cM1ZMgQFRcXa8OGDdZDZJcahKoylgsp+47ObTv3O3NwcFBUVJSioqJ0+eWX6/rrr9fq1avl5OSkmTNnWvu9/vrreuSRR/TZZ5/pyiuvlJ+fn8aMGaODBw9Kuvjf05KSEqWnp1epdqC2cNUYTMnd3V233nqr3nnnHSUmJur999+Xl5eXbrrpJmufxYsXa8iQIVqwYEG5ZbOysipdZ1Xv67N48WKNHz9e//rXv8q1p6amqmnTptUbyJ/KfiEmJiZWuJrsxIkTCggIqNF6K9tOUVGRTp48WS4MGYahpKQk6wne1bF3717t2rVLixYt0oQJE6zthw4dqvI6KrvnU1JSkrXm6sjIyNCXX36pp556SrNmzbK25+fnW88bu5ihQ4fqySef1NatW3Xs2DENHz5cXl5e6t27t2JiYnTixAm1b99eYWFh1aqttpV9R+e2tW3b9qLLenh4qE2bNtq1a5e1zdPTU3PmzNGcOXOUnJxs3Ts0evRo/fbbb+X+np7rxIkTcnBwkK+v7yWMCKg+9gjBtCZNmqTi4mK99NJLWr16tW655RZ5eHhYP7dYLNZ/YZfZvXv3Jd+/p7L1fvXVVxXux1Kdf91fddVVklThZOdt27Zp//79Gjp06KWUbFW2nnO3s3z5cuXk5NRoO2UB8tzv5O23367yOrKysqyHOct8/PHHcnBw0BVXXFHtegzDqFDPu+++e96r0M41bNgwFRUV6YknnlCLFi2sJ84PGzZM3333nX744Ycq7Q2qyl66S7FkyZJy7zdt2qS4uDgNGTLkostmZ2fr0KFDCgoKqvTz4OBgTZw4UbfeeqsOHDig3NxcdejQQc2bN9fHH38swzCsfXNycrR8+XLrlWRAfWKPEEwrKipK3bt317x582QYRoV7B1177bV69tln9dRTT2nw4ME6cOCAnnnmGUVERKioqKjG27322mu1aNEidezYUd27d9eOHTv00ksvVdiT06ZNG7m7u2vJkiXq1KmTmjRpotDQUIWGhlZYZ4cOHXT33XfrP//5jxwcHDRixAgdPXpUTzzxhMLCwjRjxowa13u24cOH6+qrr9YjjzyizMxMDRgwQLt379ZTTz2lnj17aty4cdVeZ8eOHdWmTRvNmjVLhmHIz89PX3zxRbkrrC7G399f9957r+Lj49W+fXutXr1a77zzju699161bNmyWvV4e3vriiuu0EsvvaSAgAC1atVK69at03vvvVflPXaRkZHy9fXVmjVrdMcdd1jbhw0bpmeffdb688V069ZNK1as0IIFCxQZGWk9RFVbtm/frsmTJ+umm25SQkKCHnvsMTVv3rzCPbZKSkqshypLSkp0/Phxvf7660pPT7fe0kGS+vTpo2uvvVbdu3eXr6+v9u/fr48++qhcwHnxxRd1++2369prr9U999yj/Px8vfTSSzp9+rSef/75WhsbUGU2PVUbsLHXXnvNkGR07ty5wmf5+fnGww8/bDRv3txwc3MzevXqZXz22WcVrlwpu4LnpZdeqrCOyq4aS09PNyZNmmQEBQUZHh4exsCBA40NGzYYgwcPNgYPHlxu+aVLlxodO3Y0nJ2dDUnGU089ZRhGxavGDKP0ipwXXnjBaN++veHs7GwEBAQYf//7342EhIRy/QYPHmx06dKlQq1VvSInLy/PeOSRR4zw8HDD2dnZCAkJMe69914jPT29wvqqetXYvn37jOHDhxteXl6Gr6+vcdNNNxnx8fHlxmwY579qrEuXLsbatWuNqKgow9XV1QgJCTEeffRRo7Cw0NrvQvN07naOHTtm3HDDDYavr6/h5eVlXHPNNcbevXuN8PBwY8KECVUa0/XXX29IMpYsWWJtKygoMDw9PQ0HB4cK31dlY0tLSzNuvPFGo2nTpobFYrHOeXXGUpmyba1Zs8YYN26c0bRpU8Pd3d0YOXKkcfDgwXJ9K7tqLCgoyBg8eLCxcuXKcn1nzZplREVFGb6+voarq6vRunVrY8aMGUZqamq5fp999pnRp08fw83NzfD09DSGDh1q/PTTTxf9PrhqDHXBYhhn7Z8EgEZmyJAhSk1NrdPnXtmbRYsW6Y477tC2bdtqdQ8T0BhxjhAAADAtghAAADAtDo0BAADTsukeofXr12v06NEKDQ2VxWLRZ599dsH+K1as0PDhwxUYGChvb2/169dP3377bf0UCwAA7I5Ng1BOTo569OihN954o0r9169fr+HDh2v16tXasWOHrrzySo0ePVqxsbF1XCkAALBHDebQmMVi0cqVKzVmzJhqLdelSxeNHTtWTz75ZN0UBgAA7FajvqFiSUmJsrKy5Ofnd94++fn5ys/PL7dMWlqa/P39q/xIBAAAYFuGYSgrK0uhoaHlHtp7qRp1EPr3v/+tnJwc3XzzzeftM3fuXM2ZM6ceqwIAAHUlISGhwp34L0WjPTS2dOlSTZ48WZ9//vkFb1V/7h6hjIwMtWzZUgkJCfL29r7UsgEAQD3IzMxUWFiYTp8+LR8fn1pbb6PcI7Rs2TJNmjRJn3766UWf1+Pq6lrh4YlS6fOECEIAADQutX1aS6O7oeLSpUs1ceJEffzxxxo1apStywEAAI2YTfcIZWdn69ChQ9b3R44c0c6dO+Xn56eWLVtq9uzZOn78uD788ENJpSFo/Pjxeu2119S3b18lJSVJktzd3Wt1NxkAADAHm+4R2r59u3r27KmePXtKkmbOnKmePXtaL4VPTExUfHy8tf/bb7+toqIi3XfffQoJCbG+HnzwQZvUDwAAGrcGc7J0fcnMzJSPj48yMjI4RwgAgEairn5/N7pzhAAAAGoLQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJgWQQgAAJiWTYPQ+vXrNXr0aIWGhspiseizzz676DLr1q1TZGSk3Nzc1Lp1a7311lt1XygAALBLNg1COTk56tGjh954440q9T9y5IhGjhypQYMGKTY2Vo8++qimTZum5cuX13GlAADAHjnZcuMjRozQiBEjqtz/rbfeUsuWLTVv3jxJUqdOnbR9+3a9/PLLuuGGG+qoSgAAYK8a1TlCmzdvVnR0dLm2q6++Wtu3b1dhYaGNqgIAAI2VTfcIVVdSUpKCg4PLtQUHB6uoqEipqakKCQmpsEx+fr7y8/Ot7zMzM+u8TgAA0Dg0qj1CkmSxWMq9Nwyj0vYyc+fOlY+Pj/UVFhZW5zUCAIDGoVEFoWbNmikpKalcW0pKipycnOTv71/pMrNnz1ZGRob1lZCQUB+lAgCARqBRHRrr16+fvvjii3Jta9asUVRUlJydnStdxtXVVa6urvVRHgAAaGRsukcoOztbO3fu1M6dOyWVXh6/c+dOxcfHSyrdmzN+/Hhr/ylTpiguLk4zZ87U/v379f777+u9997Tww8/bIvyAQBAI2fTPULbt2/XlVdeaX0/c+ZMSdKECRO0aNEiJSYmWkORJEVERGj16tWaMWOG3nzzTYWGhur111/n0nkAAFAjFqPsbGOTyMzMlI+PjzIyMuTt7W3rcgAAQBXU1e/vRnWyNAAAQG0iCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANMiCAEAANOyeRCaP3++IiIi5ObmpsjISG3YsOGC/ZcsWaIePXrIw8NDISEhuuOOO3Tq1Kl6qhYAANgTmwahZcuWafr06XrssccUGxurQYMGacSIEYqPj6+0/8aNGzV+/HhNmjRJv/76qz799FNt27ZNkydPrufKAQCAPbBpEHrllVc0adIkTZ48WZ06ddK8efMUFhamBQsWVNp/y5YtatWqlaZNm6aIiAgNHDhQ99xzj7Zv317PlQMAAHtgsyBUUFCgHTt2KDo6ulx7dHS0Nm3aVOky/fv317Fjx7R69WoZhqHk5GT973//06hRo867nfz8fGVmZpZ7AQAASDYMQqmpqSouLlZwcHC59uDgYCUlJVW6TP/+/bVkyRKNHTtWLi4uatasmZo2bar//Oc/593O3Llz5ePjY32FhYXV6jgAAEDjZfOTpS0WS7n3hmFUaCuzb98+TZs2TU8++aR27Nihb775RkeOHNGUKVPOu/7Zs2crIyPD+kpISKjV+gEAQOPlZKsNBwQEyNHRscLen5SUlAp7icrMnTtXAwYM0D/+8Q9JUvfu3eXp6alBgwbpueeeU0hISIVlXF1d5erqWvsDAAAAjZ7N9gi5uLgoMjJSMTEx5dpjYmLUv3//SpfJzc2Vg0P5kh0dHSWV7kkCAACoDpseGps5c6beffddvf/++9q/f79mzJih+Ph466Gu2bNna/z48db+o0eP1ooVK7RgwQIdPnxYP/30k6ZNm6bLL79coaGhthoGAABopGx2aEySxo4dq1OnTumZZ55RYmKiunbtqtWrVys8PFySlJiYWO6eQhMnTlRWVpbeeOMNPfTQQ2ratKmuuuoqvfDCC7YaAgAAaMQshsmOKWVmZsrHx0cZGRny9va2dTkAAKAK6ur3t82vGgMAALAVghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtghAAADAtp5oumJOTo3Xr1ik+Pl4FBQXlPps2bdolFwYAAFDXahSEYmNjNXLkSOXm5ionJ0d+fn5KTU2Vh4eHgoKCCEIAAKBRqNGhsRkzZmj06NFKS0uTu7u7tmzZori4OEVGRurll1+u7RoBAADqRI2C0M6dO/XQQw/J0dFRjo6Oys/PV1hYmF588UU9+uijtV0jAABAnahREHJ2dpbFYpEkBQcHKz4+XpLk4+Nj/RkAAKChq9E5Qj179tT27dvVvn17XXnllXryySeVmpqqjz76SN26davtGgEAAOpEjfYI/etf/1JISIgk6dlnn5W/v7/uvfdepaSkaOHChbVaIAAAQF2xGIZh2LqI+pSZmSkfHx9lZGTI29vb1uUAAIAqqKvf39xQEQAAmFaNglBycrLGjRun0NBQOTk5Wa8eK3sBAFAfDMNQcYmpDmygltXoZOmJEycqPj5eTzzxhEJCQqxXkAEAUBMFRSXKyCtURl6BTucWlr7yCnU6t+DP9tLX6dzSPzPPanv6b130977hth4CGqkaBaGNGzdqw4YNuuyyy2q5HABAY1ZSYigjr1DpuQVKzy0NMmk5Bda2v0JOgdJzysJNgXIKimu8zYy8wlocAcymRkEoLCxMtXWO9fz58/XSSy8pMTFRXbp00bx58zRo0KDz9s/Pz9czzzyjxYsXKykpSS1atNBjjz2mO++8s1bqAQCUMgxDmWeKlJ5ToFM5BUrPKVBa7jl/5pSFntL3GXmFqumRKotF8nF3lo+7s5q6O8vHw+Wvn92d1dTDWd7lPi/92dfDpXYHDlOpURCaN2+eZs2apbffflutWrWq8caXLVum6dOna/78+RowYIDefvttjRgxQvv27VPLli0rXebmm29WcnKy3nvvPbVt21YpKSkqKiqqcQ0AYBaGYSgzr0ipOfk6lV2gtJx8pWaX7rFJyylQana+9eey4FNUw1TTxNVJTT2c5efpoqYeLvL1KA0vTT1c1NSjNLw09fjz/Z9hxsvNSQ4OnGqB+lXly+d9fX3LnQuUk5OjoqIieXh4yNnZuVzftLS0Km28T58+6tWrlxYsWGBt69Spk8aMGaO5c+dW6P/NN9/olltu0eHDh+Xn51elbZyLy+cB2JOi4hKl5RToZHZpqEnNytepnL9+Ts0p0KnsfGvIKSyufrDxdHGUr6eL/D1d5OvpIj+P0nDj38RFvh4u8vMsDTSlocdZTd1d5OLERcmoXXX1+7vKe4TmzZtXaxuVpIKCAu3YsUOzZs0q1x4dHa1NmzZVusyqVasUFRWlF198UR999JE8PT31t7/9Tc8++6zc3d0rXSY/P1/5+fnW95mZmbU3CACoA2V7bk5mn1FKZr5OZufrZNZZr7Pep+UWqLpnKni5Osm/SWlw8W/iKn9Plz/fl/7sd87LzZmrgWG/qhyEJkyYUKsbTk1NVXFxsYKDg8u1BwcHKykpqdJlDh8+rI0bN8rNzU0rV65Uamqqpk6dqrS0NL3//vuVLjN37lzNmTOnVmsHgJowjNITiZMz85WceUbJmWeUkpWvlMwzSs7MV0rWn++z8lVQVFLl9TpYJD9PVwU0cVGgl6sC/gw3AV5//Rng6Sr/JqWBx9WJYAOUqdE5QpJUXFyslStXav/+/bJYLOrUqZOuu+46OTlVb5XnXnpvGMZ5L8cvKSmRxWLRkiVL5OPjI0l65ZVXdOONN+rNN9+sdK/Q7NmzNXPmTOv7zMxMhYWFVatGALiYgqISJWeeUVLmGSVllIacpIzS98l/Bp3kzDPKr0bA8XZzUqCXqwK9XBXk5Wb9ObCJq/XngCau8vN0kSPn1gA1UqMgtHfvXl133XVKSkpShw4dJEm///67AgMDtWrVqio9eDUgIECOjo4V9v6kpKRU2EtUJiQkRM2bN7eGIKn0nCLDMHTs2DG1a9euwjKurq5ydXWtzvAAoJyykHPidJ4SM87oREaekjLOKDHjzJ9/5ik1u6DK6/P1cFawt5uCvN0U5OWqIC/X0vderta2QC9XDkkB9aBGQWjy5Mnq0qWLtm/fLl9fX0lSenq6Jk6cqLvvvlubN2++6DpcXFwUGRmpmJgYXX/99db2mJgYXXfddZUuM2DAAH366afKzs5WkyZNJJUGMAcHB7Vo0aImQwFgcoZh6HRuoY6fztPx03k6YX2dsb4/mZ1fpfNwXJwc1MzbTc283RTs46Zm3qUBp5nPn23ebgQcoIGp0UNX3d3dtX37dnXp0qVc+969e9W7d2/l5eVVaT3Lli3TuHHj9NZbb6lfv35auHCh3nnnHf36668KDw/X7Nmzdfz4cX344YeSpOzsbHXq1El9+/bVnDlzlJqaqsmTJ2vw4MF65513qrRNrhoDzMUwDJ3KKVBCWq6OpefpWHqejp8u/fl4emn4ya3CzfxcnBwU4uOmEB83hfq4K6Spm0J83BXi81fQ8fN04U77QB2x+VVjZ+vQoYOSk5MrBKGUlBS1bdu2yusZO3asTp06pWeeeUaJiYnq2rWrVq9erfDw0lulJyYmKj4+3tq/SZMmiomJ0QMPPKCoqCj5+/vr5ptv1nPPPVeTYQCwEzn5RUpIz1X8qVzFp+UqIS1XCel51vCTV3jxoBPo5arQpu5q0dRdoU3dFNrUXSE+7mretDT0+BNyALtUoz1Cq1ev1j//+U89/fTT6tu3ryRpy5YteuaZZ/T8889r4MCB1r4Nba8Le4SAxscwDKVmFyg+LUdHU3MVl5ar+FM5ivsz9Fzs/ByLRWrm7aYWvqXBpoWvR+nPvqU/h/i4cbgKaODq6vd3jYKQg8NfN8oq+xdS2WrOfm+xWFRcXPPnx9QFghDQMJUdwjqSmqMjqTk6mpqjuFO5OpKao7hTORd9FlVTD2e19PNQmK+Hwvw8FObnrpZ+Hmrh66HQpm5cMg40cg3q0NiPP/5YawUAMJfcgiIdPpmjw6k5OnIyR4dTs0vDz8kcZeWf/3E5FosU6uOucH8Phft7qKWf559/eqilv4e83ZzPuywAnE+NgtDgwYNruw4AdqTsUNahlGz9cTLb+ufhkzk6fvr8F1OUhZ2IAE+1CvBQK39PRQR4KtzfU2F+7uzVAVDrqhyEdu/eXeWVdu/evUbFAGhcDMNQcma+fk/O0sGUbB1KydLB5GwdTMlWRl7heZdr6uGsNoFN1DrAUxGBnmod4KnWgU3U0s+Dc3UA1KsqB6HLLrtMFotFFzulqCGeFwTg0p3OLdBvSVn6PTmr9M8/f848U/nhLItFaunnoTaBTdQm0LP0z6AmahPYRH6eLvVcPQBUrspB6MiRI3VZB4AGorC4RH+czNb+xEz9llgaen5LylRyZn6l/R0sUit/T7ULbqL2wV5qG9RE7YK81DrQk707ABq8Kgehsnv7nG3fvn2Kj49XQcFfl65aLJZK+wJoeLLOFGp/YpZ+PZGhfScytS8xUweTs1VQXPnzsFr4uqtjMy+1D/ZSh2ZeBB4AjV6NTpY+fPiwrr/+eu3Zs6fc4bKyS+c5NAY0POk5BdpzPEN7T2To1xOZ+vV4ho6eyq20bxNXJ3Vs5qVOId7qGOKljs281T64iby4MguAnalREHrwwQcVERGh7777Tq1bt9bPP/+stLQ0PfTQQ3r55Zdru0YA1ZR5plB7jmVo17HT2nMsQ3uOZ+hYeuVXazVv6q5OId7qHOqtziGlrxa+7nLgaeYATKBGQWjz5s364YcfFBgYKAcHBzk6OmrgwIGaO3eupk2bptjY2NquE8B55BcVa39ilnYlnNbOhNPalXBah1NzKu3byt9DXZr7qFtzH3UN9VGXUG/5cuIyABOrURAqLi62Pv09ICBAJ06cUIcOHRQeHq4DBw7UaoEAyjtxOk+/xKcrNv60folP16/HMys9pyfMz13dWzRV9+Y+6tbCR11CfeTjzqEtADhbjYJQ165dtXv3brVu3Vp9+vTRiy++KBcXFy1cuFCtW7eu7RoB0yoqLtG+xExtP5quHfHp+iUuXYkZZyr08/Vw1mVhTdWj7NWiKZeoA0AV1CgIPf7448rJKd31/txzz+naa6/VoEGD5O/vr2XLltVqgYCZ5BUUKzY+XT8fSdP2uDTFxp9W7jnP2HJ0sKhTiJd6tfRVz5ZN1TPMV+H+HjwZHQBqoEYPXa1MWlqafH19G/z/jHnoKhqSnPwibY9L15bDp/Tz4VPaczxDhcXl/5P0cnNSZLivosJ9FRnupx5hPvJwqdG/YQCg0WpQD12tjJ+fX22tCrBbZwqLtf1oujb9karNh09pz7EMFZWUDz7NvN10eYSfekf4qXcrX7UP8uIKLgCoI/yzEqhDxSWGdh87rY0HU/XTH6n6Je50hRObW/i6q29rf/WJ8FPf1v5q4eve4PesAoC9IAgBtexYeq7W/56qDQdP6qdDqRWexRXs7aoBbQLUr42/+rb2V5ifh40qBQAQhIBLdKawWFuPpGntgZNa93uK/jhZ/h4+Xm5O6t/GXwPbBmhA2wBFBHiyxwcAGgiCEFADiRl5+n5/itYeSNFPh04pr/CvK7scLFKvlr66on2gBrULULfmPnJydLBhtQCA8yEIAVVgGIb2Hs9UzP5kfb8/Wb+eyCz3eZCXq4Z0CNSVHYLUv20ANy4EgEaCIAScR2FxiX4+nKY1+5IUsy+53I0MLRapZ1hTDe0UrCs7BKlTiBeHuwCgESIIAWc5U1isjQdT9fXeJH23P1kZeYXWzzxcHHVFu0AN7RSkKzsGKaCJqw0rBQDUBoIQTO9MYbE2HEzVl7tP6Lt9yco5607O/p4uGt45WNFdgtW/TYDcnB1tWCkAoLYRhGBKhcUl2ngoVV/sPKGYfcnKyv/rEvcQHzdd3aWZRnRtpqhWfnLkZoYAYLcIQjANwzC0Iy5dn+08rtV7kpSWU2D9rJm3m0Z1D9Go7iG6rEVT7uQMACZBEILdO5qao5Wxx7Uy9rji03Kt7f6eLhrVPUR/6xGqXi19CT8AYEIEIdil7Pwird6dqE93JGjb0XRru4eLo67p0kzX9WyuAW38ub8PAJgcQQh2o+zQ1yfbErR6T6Jy/zzp2cEiDWwXqP/r2VzRXYJ5cjsAwIrfCGj00nMKtCL2uD7ZGq+DKdnW9tYBnroxqoX+r2cLNfNxs2GFAICGiiCERskwDO06lqGPNsfpi90nVFBU+kR3d2dHXds9RGN7hyky3JebHAIALogghEblTGGxVu06oY82x2nP8Qxre+cQb93Wp6X+dlmovN14vAUAoGoIQmgUkjPP6KPNcfp4a7z1sncXJwdd2z1E4/qG67Kwpuz9AQBUG0EIDdre4xl6d8Nhfbk7UUUlhiQp1MdN4/q10tjeYfLzdLFxhQCAxowghAbHMAyt+/2k3tlwWD8dOmVtv7yVnyYOaKXozsFc9g4AqBUEITQYxSWGVu9J1Py1f2h/YqYkydHBolHdQnTXoNbq1sLHxhUCAOwNQQg2V1BUopWxx7Rg7R86eqr0zs8eLo66pXdL3TmwlVr4eti4QgCAvSIIwWYKikr06Y4Ezf/xDx0/nSdJaurhrDv6R2hC/3A19eD8HwBA3SIIod4VFpfo/20vH4ACvVx1zxWtdevlLeXpyl9LAED94DcO6k1xiaFVu47r1ZiD1oefBnm56t4hbXTr5S3l5uxo4woBAGZDEEKdMwxD3+9P0UvfHtCB5CxJUkATF00d0la39SEAAQBshyCEOrUz4bT+tXq/th5JkyR5uznpnsFtdMeAVjz8FABgc/wmQp1ISMvVC9/8pi93J0qSXJ0cdOfACE0Z3EY+7jwCAwDQMBCEUKty8ou0YO0fWrjhsAqKSmSxSDf0aqGHotsrxMfd1uUBAFAOQQi1wjAMfbbzuJ7/+jclZ+ZLkvq38dfjozqrc6i3jasDAKByBCFcst+Ts/T4Z3ut5wGF+bnr8VGdFd05mAehAgAaNIIQaiwnv0ivf39Q7208oqISQ+7Ojrr/qraaNDCCK8EAAI0CQQg18uNvKXps5R6dyDgjSYruHKwnR3fmcRgAgEaFIIRqOZWdrzlf7NOqXSckSS183fXMdV10VcdgG1cGAED1EYRQJYZhaNWuE3p61a9Kzy2Ug0WaPKi1pg9rx/2AAACNFr/BcFGp2fl6fOVeffNrkiSpYzMvvXhjd3Vv0dS2hQEAcIkcbF3A/PnzFRERITc3N0VGRmrDhg1VWu6nn36Sk5OTLrvssrot0OS+2Zuoq19dr29+TZKTg0Uzh7fXFw8MJAQBAOyCTYPQsmXLNH36dD322GOKjY3VoEGDNGLECMXHx19wuYyMDI0fP15Dhw6tp0rNJzu/SA9/uktTFv+iUzkF6tjMS5/dN0DThraTs6PN8zMAALXCYhiGYauN9+nTR7169dKCBQusbZ06ddKYMWM0d+7c8y53yy23qF27dnJ0dNRnn32mnTt3VnmbmZmZ8vHxUUZGhry9udFfZXYlnNaDn8Tq6KlcOVikewa30fRh7eTqxCXxAADbqKvf3zb7p31BQYF27Nih6Ojocu3R0dHatGnTeZf74IMP9Mcff+ipp56q6xJNp6TE0IK1f+iGBZt09FSuQn3ctPSuvnrkmo6EIACAXbLZydKpqakqLi5WcHD5y66Dg4OVlJRU6TIHDx7UrFmztGHDBjk5Va30/Px85efnW99nZmbWvGg7lpZToBnLdmrd7yclSSO7NdPc67vLx4MHpAIA7JfNrxo79xEMhmFU+liG4uJi3XbbbZozZ47at29f5fXPnTtXc+bMueQ67dkv8em6f8kvOpFxRq5ODprzty4a2zuMx2MAAOyezc4RKigokIeHhz799FNdf/311vYHH3xQO3fu1Lp168r1P336tHx9feXo+NchmpKSEhmGIUdHR61Zs0ZXXXVVhe1UtkcoLCyMc4RUGjoXbTqqf63er8JiQ638PTT/9kgekgoAaHDq6hwhm+0RcnFxUWRkpGJiYsoFoZiYGF133XUV+nt7e2vPnj3l2ubPn68ffvhB//vf/xQREVHpdlxdXeXq6lq7xduBM4XFemzlXi3/5Zik0kNhL9zQXV5uHAoDAJiHTQ+NzZw5U+PGjVNUVJT69eunhQsXKj4+XlOmTJEkzZ49W8ePH9eHH34oBwcHde3atdzyQUFBcnNzq9COC0vKOKN7Fu/QroTTcrBIj47spEkDIzgUBgAwHZsGobFjx+rUqVN65plnlJiYqK5du2r16tUKDw+XJCUmJl70nkKonh1x6ZqyeIdOZuXLx91Zb97WSwPbBdi6LAAAbMKm9xGyBTPfR+jL3Sc08//tUkFRiToEe2nh+EiF+3vauiwAAC7K7s4RQv0xDEML1v2hF785IEka1ilIr93SU56uTD8AwNz4TWjnCotL9MRne/XJtgRJ0h0DWunxUZ3l6MD5QAAAEITsWF5Bse77+Bf98FuKHCzSk9d21sQBlV9dBwCAGRGE7FRGbqEm/Xebtsely83ZQf+5tZeGdw6++IIAAJgIQcgOpWSe0fj3t+q3pCx5uznp/Ym9FdXKz9ZlAQDQ4BCE7ExCWq5ue3eLEtLyFOTlqv/eebk6hZjr6jgAAKqKIGRH4k7l6NaFW3Qi44zC/T300Z191NLfw9ZlAQDQYBGE7MSR1NIQlJR5Rm0CPbX0rr4K8nazdVkAADRoBCE7cCglW7e9s0UpWflqF9RES+7qoyAvQhAAABdDEGrkjqTm6NZ3tuhkVr46BHtpyV19FNCEh8wCAFAVBKFG7PjpPN3+Zwjq2MxLH9/VV36eLrYuCwCARsPB1gWgZlKyzuj2d0pPjG4d4KmPJvUhBAEAUE0EoUYoPadA497dqqOnctW8qbsWT+6jQC8OhwEAUF0EoUYmr6BYdyzapgPJWQryctXHd/VRaFN3W5cFAECjRBBqRIqKS/TA0ljtTDitph7OWjy5j8L9PW1dFgAAjRZBqJEwDENPf/GrvtufLBcnB707Pkrtg71sXRYAAI0aQaiRWLDuDy3eEi+LRXpt7GU8OwwAgFpAEGoEPt95XC9+c0CS9MSozhrRLcTGFQEAYB8IQg3c7mOn9c//7ZYkTRoYoTsHRti4IgAA7AdBqAFLyTyjuz/cofyiEl3VMUiPjuxk65IAALArBKEG6kxhse7+aIeSMs+obVATvXbLZXJ0sNi6LAAA7ApBqAEyDEOPrdyrnQmn5ePurHfHR8nLzdnWZQEAYHcIQg3Q4i1xWv7LMTk6WPTmbb3UKoB7BQEAUBcIQg3M7mOn9eyX+yVJs0d01MB2ATauCAAA+0UQakAycgs1dckvKiguUXTnYE3iCjEAAOoUQaiBMAxDD/9vl46l5ynMz10v3dRDFgsnRwMAUJcIQg3EexuPKGZfslwcHTT/tkj5uHNyNAAAdY0g1ADsOZah57/+TZL0xOjO6tbCx8YVAQBgDgQhG8srKNb0ZbEqKjE0omsz/b1PS1uXBACAaRCEbOz5r/frj5M5CvJy1b+u78Z5QQAA1COCkA2t+/2k/rs5TpL08k095OvpYuOKAAAwF4KQjaTnFOgfn+6SJE3s30pXtA+0cUUAAJgPQchGHv98r1Ky8tUm0FOPXNPR1uUAAGBKBCEb+PbXJH21O1GODhbNG9tT7i6Oti4JAABTIgjVs8wzhXry872SpLuvaM2l8gAA2BBBqJ698PVvSs7MVyt/Dz04tJ2tywEAwNQIQvVo65E0Lfk5XpL0r//rJjdnDokBAGBLBKF6cqawWLNW7JYk3dI7TP3b8FR5AABsjSBUT+av/UOHT+Yo0MtVs0d0snU5AABABKF6kZCWq7fW/SFJenp0F/l48EBVAAAaAoJQPXjuq30qKCpR/zb+Gtmtma3LAQAAfyII1bENB0/q21+T5ehg0VOju/AsMQAAGhCCUB0qLC7RnC/2SZLG9Q1Xh2ZeNq4IAACcjSBUhz7cHKdDKdny83TRjGHtbV0OAAA4B0GojqRm52tezO+SpH9c3YETpAEAaIAIQnXkP98fVFZ+kbo299bNUWG2LgcAAFSCIFQH4k/l6uOtpXeQfnREJzk6cII0AAANEUGoDrz63e8qLDY0qF2A+rflDtIAADRUBKFatj8xU5/tPC5JeuSajjauBgAAXAhBqJa99O0BGYZ0bfcQdW3uY+tyAADABRCEatHWI2n64bcUOTlY9FB0B1uXAwAALsLmQWj+/PmKiIiQm5ubIiMjtWHDhvP2XbFihYYPH67AwEB5e3urX79++vbbb+ux2vMzDEMvfPObJGls7zBFBHjauCIAAHAxNg1Cy5Yt0/Tp0/XYY48pNjZWgwYN0ogRIxQfH19p//Xr12v48OFavXq1duzYoSuvvFKjR49WbGxsPVde0YaDqdoRly43ZwdNG9rO1uUAAIAqsBiGYdhq43369FGvXr20YMECa1unTp00ZswYzZ07t0rr6NKli8aOHasnn3yySv0zMzPl4+OjjIwMeXt716juytz81mZtPZqmOwdE6MnRnWttvQAAoO5+f9tsj1BBQYF27Nih6Ojocu3R0dHatGlTldZRUlKirKws+fn5nbdPfn6+MjMzy71q28+HT2nr0TS5ODro7ita1/r6AQBA3bBZEEpNTVVxcbGCg4PLtQcHByspKalK6/j3v/+tnJwc3XzzzeftM3fuXPn4+FhfYWG1f5fnN348JEm6MaqFmvm41fr6AQBA3bD5ydIWS/m7LhuGUaGtMkuXLtXTTz+tZcuWKSgo6Lz9Zs+erYyMDOsrISHhkms+266E09pwMFWODhbdO7hNra4bAADULSdbbTggIECOjo4V9v6kpKRU2Et0rmXLlmnSpEn69NNPNWzYsAv2dXV1laur6yXXez5le4OuuyxUYX4edbYdAABQ+2y2R8jFxUWRkZGKiYkp1x4TE6P+/fufd7mlS5dq4sSJ+vjjjzVq1Ki6LvOCfkvKVMy+ZFks0tQhbW1aCwAAqD6b7RGSpJkzZ2rcuHGKiopSv379tHDhQsXHx2vKlCmSSg9rHT9+XB9++KGk0hA0fvx4vfbaa+rbt691b5K7u7t8fOr/Ls7zf/xDkjSya4jaBjWp9+0DAIBLY9MgNHbsWJ06dUrPPPOMEhMT1bVrV61evVrh4eGSpMTExHL3FHr77bdVVFSk++67T/fdd5+1fcKECVq0aFG91n7idJ6+2pMoSZp6JecGAQDQGNn0PkK2UFv3IXj+69/01ro/1K+1v5be3bcWKwQAAOeyu/sINWZ5BcVaurV0T9XEAa1sWwwAAKgxglANrIw9roy8QoX5uWtYpwtf4QYAABouglA1GYahRZuOSJIm9GslR4eL3/MIAAA0TAShavrp0Cn9npwtDxdH3RRV+3epBgAA9YcgVE0f/FS6N+jGyBbycXe2cTUAAOBSEISq4Whqjn44kCJJmtC/lW2LAQAAl4wgVA0fbo6TYUhDOgSqTSA3UAQAoLEjCFVRflGxVsQekySN7xdu42oAAEBtIAhVUcy+ZJ3OLVSwt6uuaBdo63IAAEAtIAhV0bJtCZJKT5J2cuRrAwDAHvAbvQqOpedq46FUSdLNXDIPAIDdIAhVwf92HJNhSP1a+yvc39PW5QAAgFpCELqIkhJDn24vPUl6bG/2BgEAYE8IQhfx0x+pOn46T15uTrqmazNblwMAAGoRQegiyk6SHnNZc7k5O9q4GgAAUJsIQheQnlOgNb8mS+KwGAAA9oggdAFf7j6hguISdQrxVtfmPrYuBwAA1DKC0AV8sStRkvR/PZvbuBIAAFAXCELnkZiRp21xaZKkUd1DbFwNAACoCwSh8/hqd6IMQ4oK91VoU3dblwMAAOoAQeg8vtxdelhsdI9QG1cCAADqCkGoEglpudqZcFoOFmlEN+4dBACAvSIIVaJsb1Df1v4K8nKzcTUAAKCuEIQq8cWuE5Kka7tzWAwAAHtGEDrHHyeztS8xU04OFh6pAQCAnSMInePLP+8dNLBdgPw8XWxcDQAAqEsEobMYhqEvdnNYDAAAsyAIneVQSrYOpWTLxdFB0V2CbV0OAACoYwShs3z/W4okqV8bf3m7Odu4GgAAUNcIQmf5fn/pk+aHdQqycSUAAKA+EIT+lJ5ToB1x6ZKkKzsShAAAMAOC0J/W/X5SJYbUsZmXWvh62LocAABQDwhCf/ruz8NiQzksBgCAaRCEJBUWl2jd7yclSUM7cbUYAABmQRCStO1omrLOFMnf00U9WjS1dTkAAKCeEIQk/bC/9LL5KzsGydHBYuNqAABAfSEI6a/7Bw3lajEAAEzF9EHo8MlsHUnNkbOjRQPbBdi6HAAAUI9MH4S+//OwWN/W/vLibtIAAJgKQei30svmr+KwGAAApmPqIJRbUPTX3aQ7EIQAADAbUweh7UfTVVhsqHlTd4X7czdpAADMxtRBaNMfpyRJ/dv4y2LhsnkAAMzG1EFo8x+pkqR+bfxtXAkAALAF0wahjLxC7TmeIYkgBACAWZk2CO2IS1eJIbUO8FSIj7utywEAADZg2iD085HS84PYGwQAgHmZNghtO5ImSerfhrtJAwBgVqYNQr8nZ0uS+rb2s3ElAADAVkwbhCSpYzMv+TdxtXUZAADARmwehObPn6+IiAi5ubkpMjJSGzZsuGD/devWKTIyUm5ubmrdurXeeuutGm+b84MAADA3mwahZcuWafr06XrssccUGxurQYMGacSIEYqPj6+0/5EjRzRy5EgNGjRIsbGxevTRRzVt2jQtX768Rtvn/CAAAMzNYhiGYauN9+nTR7169dKCBQusbZ06ddKYMWM0d+7cCv0feeQRrVq1Svv377e2TZkyRbt27dLmzZurtM3MzEz5+PgofMb/067/b4x83HniPAAADV3Z7++MjAx5e3vX2npttkeooKBAO3bsUHR0dLn26Ohobdq0qdJlNm/eXKH/1Vdfre3bt6uwsLBa2+8c4k0IAgDA5JxsteHU1FQVFxcrODi4XHtwcLCSkpIqXSYpKanS/kVFRUpNTVVISEiFZfLz85Wfn299n5FRejfp7s1clJmZeanDAAAA9aDsd3ZtH8iyWRAqc+7DTg3DuOADUCvrX1l7mblz52rOnDkV2l8Yf5VeqG6xAADApk6dOiUfH59aW5/NglBAQIAcHR0r7P1JSUmpsNenTLNmzSrt7+TkJH//yq8Amz17tmbOnGl9f/r0aYWHhys+Pr5Wv8iGLjMzU2FhYUpISKjVY6sNHeNm3GbAuBm3GWRkZKhly5by86vd+//ZLAi5uLgoMjJSMTExuv76663tMTExuu666ypdpl+/fvriiy/Kta1Zs0ZRUVFydq78fB9XV1e5ula8V5CPj4+p/gKV8fb2ZtwmwrjNhXGbi1nH7eBQu6c32/Ty+ZkzZ+rdd9/V+++/r/3792vGjBmKj4/XlClTJJXuzRk/fry1/5QpUxQXF6eZM2dq//79ev/99/Xee+/p4YcfttUQAABAI2bTc4TGjh2rU6dO6ZlnnlFiYqK6du2q1atXKzw8XJKUmJhY7p5CERERWr16tWbMmKE333xToaGhev3113XDDTfYaggAAKARs/nJ0lOnTtXUqVMr/WzRokUV2gYPHqxffvmlxttzdXXVU089VenhMnvGuBm3GTBuxm0GjLt2x23TGyoCAADYks2fNQYAAGArBCEAAGBaBCEAAGBaBCEAAGBadhmE5s+fr4iICLm5uSkyMlIbNmy4YP9169YpMjJSbm5uat26td566616qrR2VWfca9eulcViqfD67bff6rHiS7d+/XqNHj1aoaGhslgs+uyzzy66jD3Md3XHbQ/zPXfuXPXu3VteXl4KCgrSmDFjdODAgYsu19jnuybjtof5XrBggbp37269aWC/fv309ddfX3CZxj7XUvXHbQ9zXZm5c+fKYrFo+vTpF+xXG3Nud0Fo2bJlmj59uh577DHFxsZq0KBBGjFiRLn7EZ3tyJEjGjlypAYNGqTY2Fg9+uijmjZtmpYvX17PlV+a6o67zIEDB5SYmGh9tWvXrp4qrh05OTnq0aOH3njjjSr1t5f5ru64yzTm+V63bp3uu+8+bdmyRTExMSoqKlJ0dLRycnLOu4w9zHdNxl2mMc93ixYt9Pzzz2v79u3avn27rrrqKl133XX69ddfK+1vD3MtVX/cZRrzXJ9r27ZtWrhwobp3737BfrU254adufzyy40pU6aUa+vYsaMxa9asSvv/85//NDp27Fiu7Z577jH69u1bZzXWheqO+8cffzQkGenp6fVQXf2QZKxcufKCfexlvs9WlXHb43ynpKQYkox169adt489zndVxm2P820YhuHr62u8++67lX5mj3Nd5kLjtre5zsrKMtq1a2fExMQYgwcPNh588MHz9q2tOberPUIFBQXasWOHoqOjy7VHR0dr06ZNlS6zefPmCv2vvvpqbd++XYWFhXVWa22qybjL9OzZUyEhIRo6dKh+/PHHuiyzQbCH+b4U9jTfGRkZknTBBzDa43xXZdxl7GW+i4uL9cknnygnJ0f9+vWrtI89znVVxl3GXub6vvvu06hRozRs2LCL9q2tOberIJSamqri4uIKT68PDg6u8NT6MklJSZX2LyoqUmpqap3VWptqMu6QkBAtXLhQy5cv14oVK9ShQwcNHTpU69evr4+SbcYe5rsm7G2+DcPQzJkzNXDgQHXt2vW8/extvqs6bnuZ7z179qhJkyZydXXVlClTtHLlSnXu3LnSvvY019UZt73MtSR98skn+uWXXzR37twq9a+tObf5IzbqgsViKffeMIwKbRfrX1l7Q1edcXfo0EEdOnSwvu/Xr58SEhL08ssv64orrqjTOm3NXua7Ouxtvu+//37t3r1bGzduvGhfe5rvqo7bXua7Q4cO2rlzp06fPq3ly5drwoQJWrdu3XlDgb3MdXXGbS9znZCQoAcffFBr1qyRm5tblZerjTm3qz1CAQEBcnR0rLAXJCUlpUJqLNOsWbNK+zs5Ocnf37/Oaq1NNRl3Zfr27auDBw/WdnkNij3Md21prPP9wAMPaNWqVfrxxx/VokWLC/a1p/muzrgr0xjn28XFRW3btlVUVJTmzp2rHj166LXXXqu0rz3NdXXGXZnGONc7duxQSkqKIiMj5eTkJCcnJ61bt06vv/66nJycVFxcXGGZ2ppzuwpCLi4uioyMVExMTLn2mJgY9e/fv9Jl+vXrV6H/mjVrFBUVJWdn5zqrtTbVZNyViY2NVUhISG2X16DYw3zXlsY234Zh6P7779eKFSv0ww8/KCIi4qLL2MN812TclWls810ZwzCUn59f6Wf2MNfnc6FxV6YxzvXQoUO1Z88e7dy50/qKiorS7bffrp07d8rR0bHCMrU259U6tboR+OSTTwxnZ2fjvffeM/bt22dMnz7d8PT0NI4ePWoYhmHMmjXLGDdunLX/4cOHDQ8PD2PGjBnGvn37jPfee89wdnY2/ve//9lqCDVS3XG/+uqrxsqVK43ff//d2Lt3rzFr1ixDkrF8+XJbDaFGsrKyjNjYWCM2NtaQZLzyyitGbGysERcXZxiG/c53dcdtD/N97733Gj4+PsbatWuNxMRE6ys3N9faxx7nuybjtof5nj17trF+/XrjyJEjxu7du41HH33UcHBwMNasWWMYhn3OtWFUf9z2MNfnc+5VY3U153YXhAzDMN58800jPDzccHFxMXr16lXuMtMJEyYYgwcPLtd/7dq1Rs+ePQ0XFxejVatWxoIFC+q54tpRnXG/8MILRps2bQw3NzfD19fXGDhwoPHVV1/ZoOpLU3bp6LmvCRMmGIZhv/Nd3XHbw3xXNl5JxgcffGDtY4/zXZNx28N833nnndb/nwUGBhpDhw61hgHDsM+5Nozqj9se5vp8zg1CdTXnFsP488wiAAAAk7Grc4QAAACqgyAEAABMiyAEAABMiyAEAABMiyAEAABMiyAEAABMiyAEAABMiyAEAJImTpyoMWPGWN8PGTJE06dPt1k9AOoHQQhAgzdkyBBZLBZZLBY5ODgoODhYN910k+Li4mptG6+99poWLVpUa+sD0DgQhAA0CnfddZcSExN1/Phxff7550pISNDf//73Wlu/j4+PmjZtWmvrA9A4EIQA2NyQIUN0//336/7771fTpk3l7++vxx9/XGc/AcjDw0PNmjVTSEiI+vbtq/vuu0+//PKL9fP09HTdfvvtCgwMlLu7u9q1a6cPPvjA+vmePXt01VVXyd3dXf7+/rr77ruVnZ1t/fzcQ2MAzIEgBKBB+O9//ysnJyf9/PPPev311/Xqq6/q3XffrbRvWlqaPv30U/Xp08fa9sQTT2jfvn36+uuvtX//fi1YsEABAQGSpNzcXF1zzTXy9fXVtm3b9Omnn+q7777T/fffXy9jA9BwOdm6AACQpLCwML366quyWCzq0KGD9uzZo1dffVV33XWXJGn+/Pl69913ZRiGcnNz1b59e3377bfW5ePj49WzZ09FRUVJklq1amX9bMmSJcrLy9OHH34oT09PSdIbb7yh0aNH64UXXlBwcHD9DRRAg8IeIQANQt++fWWxWKzv+/Xrp4MHD6q4uFiSdPvtt2vnzp3atWuXNm7cqLZt2yo6OlpZWVmSpHvvvVeffPKJLrvsMv3zn//Upk2brOvav3+/evToYQ1BkjRgwACVlJTowIED9TRCAA0RQQhAo+Dj46O2bduqbdu2GjBggN577z0dPHhQy5YtkySNGDFCcXFxmj59uk6cOKGhQ4fq4YcfliQZhlEuZJ3tfO0AzIEgBKBB2LJlS4X37dq1k6OjY6X9y9rz8vKsbYGBgZo4caIWL16sefPmaeHChZKkzp07a+fOncrJybH2/emnn+Tg4KD27dvX9lAANCIEIQANQkJCgmbOnKkDBw5o6dKl+s9//qMHH3zQ+nlubq6SkpKUlJSkXbt2aerUqXJzc1N0dLQk6cknn9Tnn3+uQ4cO6ddff9WXX36pTp06SSo9rObm5qYJEyZo7969+vHHH/XAAw9o3LhxnB8EmBwnSwNoEMaPH6+8vDxdfvnlcnR01AMPPKC7777b+vk777yjd955R5Lk6+ur7t27a/Xq1erQoYMkycXFRbNnz9bRo0fl7u6uQYMG6ZNPPpFUeun9t99+qwcffFC9e/eWh4eHbrjhBr3yyiv1P1AADYrFOPtGHQBgA0OGDNFll12mefPm2boUACbDoTEAAGBaBCEAAGBaHBoDAACmxR4hAABgWgQhAABgWgQhAABgWgQhAABgWgQhAABgWgQhAABgWgQhAABgWgQhAABgWgQhAABgWv8/HqIgxEhMQasAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "def alpha(pBsoil, sN, pVal, sRootZoneMax):\n",
    "    \"\"\"\n",
    "    Calculates the value of alpha for a given value of pBsoil, sN, pVal, and sRootZoneMax.\n",
    "\n",
    "    Args:\n",
    "        pBsoil: The bulk density of the soil.\n",
    "        sN: The unavailable water fraction.\n",
    "        pVal: The water content of the soil.\n",
    "        sRootZoneMax: The maximum root zone storage capacity.\n",
    "\n",
    "    Returns:\n",
    "        The value of alpha.\n",
    "    \"\"\"\n",
    "    pCmax = sRootZoneMax * (pBsoil + 1)\n",
    "    coeff1 = 1.0 - ((pBsoil + 1.0) * (sN) / pCmax)\n",
    "    exp = 1.0 / (pBsoil + 1.0)\n",
    "    ctPrev = pCmax * (1.0 - np.power(coeff1, exp))\n",
    "    uT1 = max(pVal - pCmax + ctPrev, 0.0)\n",
    "    dummy = min((ctPrev + pVal - uT1) / pCmax, 1.0)\n",
    "    coeff2 = (1.0 - dummy)\n",
    "    exp2 = (pBsoil + 1.0)\n",
    "    xn = (pCmax / (pBsoil + 1.0)) * (1.0 - np.power(coeff2, exp2))\n",
    "    uT2 = max(pVal - uT1 - (xn - sN), 0)\n",
    "    alpha = (uT1 + uT2) / pVal\n",
    "    return alpha\n",
    "\n",
    "# Set the unavailable water fraction (sN), the water content of the soil (pVal), and the maximum root zone storage capacity (sRootZoneMax)\n",
    "sN = 100\n",
    "pVal = 40.0\n",
    "sRootZoneMax = 150.0\n",
    "\n",
    "# Create a list of pBsoil values\n",
    "pBsoil_values = np.linspace(0.0, 3.0, 100)\n",
    "\n",
    "# Calculate the alpha values for each pBsoil value\n",
    "alpha_values = []\n",
    "for pBsoil in pBsoil_values:\n",
    "    alpha_value = alpha(pBsoil, sN, pVal, sRootZoneMax)\n",
    "    alpha_values.append(alpha_value)\n",
    "\n",
    "# Plot the results\n",
    "\n",
    "# Set axis limits to [0, 1]\n",
    "plt.xlim(0, 4)\n",
    "plt.ylim(0, 1.2)\n",
    "plt.plot(pBsoil_values, alpha_values)\n",
    "plt.xlabel(\"pBsoil\")\n",
    "plt.ylabel(\"alpha\")\n",
    "plt.title(\"Variation of alpha with pBsoil\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "c01e6bdf-ab76-479d-8f06-32386447beea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.9999749999999992,\n",
       " 0.9999749999999992,\n",
       " 0.9999749999999992,\n",
       " 0.9999749999999992,\n",
       " 0.9999749999999992,\n",
       " 0.9999749999999992,\n",
       " 0.9999749999999992,\n",
       " 0.999976499616956,\n",
       " 0.9999801045508931,\n",
       " 0.999983638162945]"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "alpha_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "8860d9b5-2eeb-4456-8d2b-f70caa46623c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.        , 0.33333333, 0.66666667, 1.        , 1.33333333,\n",
       "       1.66666667, 2.        , 2.33333333, 2.66666667, 3.        ])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pBsoil_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acd8f41a-a722-4c39-8d33-8083cc75da65",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.12"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}