Created
March 26, 2025 11:53
-
-
Save grinsted/d5dcfdc729f3db804f78673527d5939d to your computer and use it in GitHub Desktop.
shallow ice approximation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"# Glacier SIA" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 29, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"%matplotlib inline\n", | |
"import numpy as np\n", | |
"import matplotlib.pyplot as plt\n", | |
"from IPython import display \n", | |
"from time import sleep\n", | |
"plt.style.use('https://www.glaciology.net/aslak.mplstyle')\n", | |
"plt.rcParams['figure.figsize']=[7,3]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 30, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"#physical constants\n", | |
"sec_per_year = 365.25 * 24 * 60 * 60\n", | |
"rhoi = 917\n", | |
"n = 3\n", | |
"R = 8.31446261815324 # J/mol/K\n", | |
"def A_fun(T):\n", | |
" A1 = 3.985e-13 * np.exp(-60e3 / (R * T))\n", | |
" A2 = 1.916e3 * np.exp(-139e3 / (R * T))\n", | |
" return np.maximum(A1, A2)\n", | |
"A = A_fun(273.15-10) #assume constant temperature -5C\n", | |
"g = 9.82\n" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 31, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"Text(0.5, 0, 'bdot (m/yr)')" | |
] | |
}, | |
"execution_count": 31, | |
"metadata": {}, | |
"output_type": "execute_result" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAvAAAAE1CAYAAACBa9TgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3IUlEQVR4nO3deVyU9f7+8QtlKSxzyV1pJEU7goIZGmn7plmnjtXJr5WWhueUlWmZndR20zJbzPY0cyH3TDPNyjX3CBVRRBEZXABRVPZl7t8fHvnJSdmEuWfmfj0fDx9H5nPDXHPejl3ecy9ehmEYAgAAAOAWapkdAAAAAEDFUeABAAAAN0KBBwAAANwIBR4AAABwIxR4AAAAwI1Q4AEAAAA3QoEHAAAA3Ii32QFqUmZmpvbt2yc/Pz+zowAAAADnlZ+fryuvvFL16tUrd1uP3gO/b98+Rf+x2ewYcJKCggIVFBSYHQNOxMythXlbC/O2HqvPPDk5Wfv27avQth69B97Pz0/xMZv12ONPyNvbo18qdPpfrpL4xMVCmLm1MG9rYd7Ww8wrzqP3wEtS6pGD+vmHuWbHAAAAAKqFxxd4SZrz7ec6mXnc7BgAAADABbNEgc/JztJ333xqdgwAAADgglmiwEvSr0sXKmlvvNkxAAAAgAtimQJvGIamfjJBhmGYHQUAAACoMssUeEnateNPbVi9wuwYAAAAQJVZqsBL0vQvPlR+Xp7ZMQAAAIAqsVyBz0hP1aLZ08yOAQAAAFSJUwv8a6+9Ji8vL8XGxpY8lpCQoIiICAUFBSk8PFxxcXEXvFaeRXO+VXrq4ep5UQAAAIATOa3AR0dHa+PGjQoICCj1+ODBgxUZGak9e/ZoxIgRGjhw4AWvlaewIF/Tv/jwwl8UAAAA4GRehhMuy5Kfn68bb7xRs2bN0k033aQlS5YoODhYaWlpCgoK0tGjR+Xt7S3DMNSsWTNt3LhR/v7+VVqz2WwlzxsbG6txo59TQda5b+I0avwnatehU02/fDhJQUGBJMnX19fkJHAWZm4tzNtamLf1WH3mO3fulK+vr4KDg8vd1il74MeMGaOHH35YrVu3LvW43W5X8+bN5e3tLUny8vJSQECAkpOTq7x2RlRUlIYMGaLsrOzz5pr77edcVhIAAABuxbumn2DDhg3asmWLxo0bd851Ly+vUl+fXairuiZJffv2VUhIiMaPGab8U8fO+dzxO7dp1/Y/FBZ+XfkvBG7Dz8/P7AhwMmZuLczbWpi39Vh15pX55KHG98CvXr1au3fvVuvWrWWz2ZSSkqI77rhDP/30k1q1aqWUlBQVFRVJOl3C7Xa7AgICqrxWWd9N/UQOh6P6XjAAAABQg2q8wI8cOVKHDh1SUlKSkpKS1LJlSy1fvlw9e/ZU48aNFRYWphkzZkiS5s+fL5vNJpvNVuW1ytq/N14b1/5aba8XAAAAqEk1fghNeT7//HMNGDBAY8eOVd26dTVt2rQLXqus2d98qq7db1Lt2qb/3wEAAACUyemNNSkpqdTX7dq104YNG865bVXXKutwSrJW/7xEN/e8t1p+HgAAAFBTLHcn1vOZO/1LFRTkmx0DAAAAKBMF/r8y0lO1Ysl8s2MAAAAAZaLAn2Vh1FTl5pz/uvEAAACA2SjwZzmZeVw/LogyOwYAAABwXhT4/7F47nSdOplpdgwAAADgnCjw/yM3J1uLZlf9kpQAAABATaLAn8NP38/R0bQjZscAAAAA/oICfw6FBfmaMvlds2MAAAAAf0GBP4+t61dry/pVZscAAAAASqHAl2HKx+8qLzfH7BgAAABACQp8GTLSUzV3+pdmxwAAAABKUODL8eP8WTqQmGB2DAAAAEASBb5cDkexvvhgrBwOh9lRAAAAAAp8RSTs2qHffvre7BgAAAAABb6iZn41SSeOHzM7BgAAACyOAl9B2VmnNP2LD8yOAQAAAIujwFfCml+WKvbPLWbHAAAAgIVR4Cvpy4/eVm5OttkxAAAAYFEU+Eo6nJKsd14ZroL8PLOjAAAAwIIo8FWwM2arJr4xUkWFhWZHAQAAgMVQ4KsoetM6TRo/Ro7iYrOjAAAAwEIo8Bdgw+oV+uz9N7nJEwAAAJyGAn+BVi1frG8+fU+GYZgdBQAAABZAga8Gy76fragpn5gdAwAAABZAga8m3383VQujppodAwAAAB6OAl+NoqZM1m/LFpkdAwAAAB6MAl/Nvnh/rGK2bDA7BgAAADwUBb6aORzFmvjGi9q/d7fZUQAAAOCBKPA1IC83R+NeHqqjaUfMjgIAAAAPQ4GvIcePHdXbLz+j7KxTZkcBAACAB6HA1yB7UqImvPq8CgsKzI4CAAAAD0GBr2E7t/2hT997nbu1AgAAoFpQ4J1g3W/LFDVlMndrBQAAwAWjwDvJotnT9OHY/ygnO8vsKAAAAHBjFHgnWr9qhUY++YiS9sabHQUAAABuyikF/vbbb1fHjh0VGhqqHj16KCYmpmTNZrOpffv2Cg0NVWhoqGbPnl2ylpCQoIiICAUFBSk8PFxxcXEVWnNlRw7Z9fIzj2nFkvkcUgMAAIBKc0qBnzNnjrZv366YmBgNHz5cjz/+eKn1efPmKSYmRjExMfrnP/9Z8vjgwYMVGRmpPXv2aMSIERo4cGCF1lxdYWGBvvzwbQ6pAQAAQKV5O+NJ6tWrV/L7EydOqFat8v/dkJaWpujoaP3888+SpD59+mjIkCFKSkqSv7//eddsNlupn+PKe7nXr1qhffG79PRLb+iKwCCz47i9Ai7XaTnM3FqYt7Uwb+ux+swLCgrk6+tboW2ddgz8o48+qlatWmnUqFGaNm1aqbV+/fopJCREgwYNUnp6uiTJbrerefPm8vY+/W8MLy8vBQQEKDk5ucy1M6KiojRkyBBlZ2U76RVWTerhFI196RllZ500OwoAAADcgFP2wEvSt99+K0maNm2aXnjhBS1dulSStGbNGgUEBKiwsFCjRo1S//79S9a8vLxK/Yyz96aXtSZJffv2VUhIiMaPGab8U8eq/fVUp5zsU1r2/Rz938CnzI7iEfz8/MyOACdj5tbCvK2FeVuPVWde0b3vkglXoenfv79WrlypjIwMSVJAQIAkycfHR0OHDtXatWslSa1atVJKSoqKiooknS7odrtdAQEBZa65q6ULZ+l4xlGzYwAAAMDF1XiBP3nypA4dOlTy9cKFC9WwYUM1aNBA2dnZyszMLFmLiopSWFiYJKlx48YKCwvTjBkzJEnz58+XzWaTzWYrc81dFeTna/7Mr82OAQAAABdX44fQnDhxQn369FFubq5q1aqlRo0aacmSJfLy8lJqaqr69Omj4uJiGYahwMDAkkNtJOnzzz/XgAEDNHbsWNWtW7fUsfNlrbmrX5cuUO/7+6lp85ZmRwEAAICL8jJc+TItFyg2NtYtjoE/W/eb79QzL71pdgy3lJ+fL8m6x85ZETO3FuZtLczbeqw+89jYWElScHBwudtyJ1YXs+63ZUrat8fsGAAAAHBRFHgX9N3UT8yOAAAAABdFgXdB0ZvWadeOP82OAQAAABdEgXdRs77+2KXvIgsAAABzUOBdVPzObYretM7sGAAAAHAxFHgXFjVlshwOh9kxAAAA4EIo8C4sef9e/b5yudkxAAAA4EIo8C5uyuR3tXHNL2bHAAAAgIugwLu47FMnNfGNkfr4nVeUk51ldhwAAACYjALvJtas+FEvDO6ruO3RZkcBAACAiSjwbiQ99bBee36wZnz5oQoLCsyOAwAAABNQ4N2MYRj6Yc50/efp/krev9fsOAAAAHAyCrybOpCYoDHPDZI9aZ/ZUQAAAOBEFHg3lpOdpbdfflbHM46aHQUAAABOQoF3c0fTjmjc6KHKy80xOwoAAACcgALvAfYn7NYHY19WcXGR2VEAAABQwyjwHiJ641pNnfyeDMMwOwoAAABqEAXeg/y8eK4Wz5thdgwAAADUIAq8h5nxxYfasPoXs2MAAACghlDgPdDH48dod2yM2TEAAABQAyjwHqiwsEBvvfS0Nv++yuwoAAAAqGYUeA+Vn5erCa8+rwWzpnBiKwAAgAehwHu476Z+oo/eHqWC/DyzowAAAKAaUOAt4PeVy/XK8EgdO5pudhQAAABcIAq8ReyLj9NLQx7Vvvg4s6MAAADgAlDgLeR4RrrGDHtCm9etNDsKAAAAqogCbzGFBfn6cOzLStgVa3YUAAAAVAEF3oIKCwv0zivDdTTtiNlRAAAAUEkUeIs6cTxD74wZprzcXLOjAAAAoBIo8BaWtG+PPn5njBwOh9lRAAAAUEEUeIvbvG6l5nz7udkxAAAAUEEUeGjBzK+17rdlZscAAABABVDgIUn6dMLrXJkGAADADVDgIen0lWnefZUr0wAAALg6CjxKZB7L0MtPD1D0pnVmRwEAAMB5OKXA33777erYsaNCQ0PVo0cPxcTElKwlJCQoIiJCQUFBCg8PV1xc3AWvoeqOHzuqcaOG6rP33lBOdpbZcQAAAPA/nFLg58yZo+3btysmJkbDhw/X448/XrI2ePBgRUZGas+ePRoxYoQGDhx4wWu4cL8tW6TnIx9S7J9bzI4CAACAs3gZhmE48wmnTZumSZMmaevWrUpLS1NQUJCOHj0qb29vGYahZs2aaePGjfL396/Sms1mK3mu2NhYjRv9nAqyjjvzJXqc2+6+X//s/2/5XXSR2VHKVFBQIEny9fU1OQmchZlbC/O2FuZtPVaf+c6dO+Xr66vg4OByt/V2Qh5J0qOPPqqVK1dKkpYtO33JQrvdrubNm8vb+3QMLy8vBQQEKDk5WXXq1KnS2pkCHxUVpc8//1z1fQrk46wX6aFWLJ6nHX9s0lMvvi7blUFmxwEAALA0pxX4b7/9VtLpPfAvvPCCli5dKul0+T7b2R8IVHVNkvr27auQkBCNHzNM+aeOXfgLsLgjh+ya8MpwvTXpGzVu2tzsOGXy8/MzOwKcjJlbC/O2FuZtPVadeWU+eXD6VWj69++vlStXKiMjQ61atVJKSoqKiooknS7hdrtdAQEBVV5DzTmReUzjRg3l5FYAAAAT1XiBP3nypA4dOlTy9cKFC9WwYUM1aNBAjRs3VlhYmGbMmCFJmj9/vmw2m2w2W5XXULNSDiRq4usvlvzjCQAAAM5V4yex2u129enTR7m5uapVq5YaNWqkCRMmKDQ0VJIUHx+vAQMGKCMjQ3Xr1tW0adPUoUOHC1o7IzY2lkNoasitd92nJ579z18OZTJTfn6+JOt+9GZFzNxamLe1MG/rsfrMY2NjJalCJ7E6/So0zkSBr1kPRz6rex54xOwYJaz+xrciZm4tzNtamLf1WH3mlSnw3IkVVTbzy4+0ed1Ks2MAAABYCgUeVWYYhj4aN0r74rkLLgAAgLM47TKS8EwF+fkaP+Y5/WvYaGVnndLxjHQdy0jXsaPpJb/39fFVWNfuurpbd7XrEFpy/X4AAABUHk0KFyzzWIbGjRpa5jYH7UlaMm+G/OtcotBrrlXnrj0UFh6hS+vWc0pGAAAAT0GBh1PlZGdp/aoVWr9qhbxq1VKPW3qq/7+eo8gDAABUEMfAwzSGw6E1K37UcwMf0O8rl//lbroAAAD4Kwo8THcy87g+HPuyxo8Zpoz0VLPjAAAAuDQKPFxG9Ma1GjboQS3/Ya4cDofZcQAAAFwSBR4uJTcnW19PGq9Xhj0he9I+s+MAAAC4HAo8XFL8zm16YfD/6ZtP3lN21imz4wAAALgMCjxclsNRrKULo/TsY//Qbz99z2E1AAAAosDDDZzMPK7PJr6pl58ZoIRdsWbHAQAAMFWlrgO/f/9+ff311/rtt9+UkpKiiy++WJ06dVKfPn3Up08f7rCJGrUvPk4vPzNAN95+tzp16VZy2UnDMGQYhgoLC+Xl5aXrbrxdfhddZHJaAACAmlHhxv2vf/1LW7du1QMPPKDx48eradOmysvL065du7Rs2TKNGzdOn376qbp161aTeQGt+nmxVv28+Lzrq39erBffeF/+dS5xYioAAADnqHCB7927tz777LO/PB4SEqIHH3xQGRkZ2rePq4bAfLt2/Kk3Rjyp/7z9EXd4BQAAHqfCx8D37t27zPWGDRsqPDz8ggMB1WHfnji9OjxSxzOOmh0FAACgWlX6oPWcnBzNmjVLiYmJKioqKnn8nXfeqdZgwIWyJyVqzLBBGj3+EzVu2tzsOAAAANWi0lehue+++7Rw4UJ5e3urTp06Jb8AV5R6KEWvDHtCh+xJZkcBAACoFpXeA5+SkqKdO3fWRBagRmSkp2rMsCc06u2PZWvTzuw4AAAAF6TSe+BDQkJ0+PDhmsgC1JiTmcf1yvBIff/dN8rLzTU7DgAAQJVVeg/86NGj1bVrV4WGhuqis661PWfOnGoNBlS33Jxszfr6Y/24IEr39X1Mt931D/n4+podCwAAoFIqXeD79++ve+65R507d1bt2rVrIhNQo04cz9A3n0zQ4rnT1efhQbrx9ru5CRkAAHAblW4tBQUF+vjjj2siC+BUGemp+uL9t7Ro9jQ98Eikut90h2rxj1IAAODiKn0M/HXXXacdO3bURBbAFKmHUvTx+DEa9sSDWvfbMjmKi82OBAAAcF6VLvAbN25Uly5d1LFjR4WHh5f8AtzdIfsBffT2KIo8AABwaZU+hOaDDz6ogRiA6zhT5OfN+FL3P/yEIm64jUNrAACAy6h0gb/hhhtqIgfgcs4U+bnTv1DPex/SDbfdpYv9uWkZAAAwV4UPoXnyySdlt9vPuWYYhubPn69Zs2ZVWzDAVRxOSdaUj9/R4Id66utJ45VyYL/ZkQAAgIVVeA98z549ddddd6l+/fqKiIhQ06ZNlZubq927d2vt2rXq2bOnXn/99ZrMCpgqLzdHy3+Yq+U/zFVw6DW6454H1CXietWuzSUoAQCA83gZhmFU5hvWrVunVatWKSUlRf7+/urYsaN69+6tyy+/vKYyVllsbKzGjxmm/FPHzI4CD9W4aQs9+fwY/a3T1WZHsaT8/HxJkp+fn8lJ4AzM21qYt/VYfeaxsbGSpODg4HK3rfSuw+7du6t79+6VTwV4oLQjB/XaC//SXX366aHH/i1fX2v+pQMAAJyn0peRBFCaYRhaMm+GXnrqEe3fu9vsOAAAwMNR4IFqYk9K1H+eHqAFs6aouLjI7DgAAMBDUeCBalRcVKTvpn6iV4ZF6sjBc1+1CQAA4EJUusBPnz5d2dnZpR5bsmRJtQUCPMGeuO0aNfRxSjwAAKh2lS7wkZGRuv7665Wamlry2JgxY867fV5enu69914FBQUpNDRUd955p5KSkkrWbTab2rdvr9DQUIWGhmr27NklawkJCYqIiFBQUJDCw8MVFxdXoTXAFZzMPK63XhqizGNHzY4CAAA8SKULfPv27TVs2DD16NFDu3efPmGvvCtRRkZGKj4+XjExMerdu7ciIyNLrc+bN08xMTGKiYnRP//5z5LHBw8erMjISO3Zs0cjRozQwIEDK7QGuIrUwwc19j/PKCc7y+woAADAQ1T6MpJeXl7q16+fmjdvrl69emnq1Kny8vI67/YXXXSRevXqVfJ1t27d9MEHH5T7PGlpaYqOjtbPP/8sSerTp4+GDBmipKQk+fv7n3fNZrOV+jmVvMw9UO2S9u3RO2OG6/nXJsjHx9fsOB6loKDA7AhwIuZtLczbeqw+84KCAvn6VqwnVHoP/JlCfNNNN+mHH37QoEGDtH9/xW8t/9FHH+nuu+8u9Vi/fv0UEhKiQYMGKT09XZJkt9vVvHlzeXuf/jeGl5eXAgIClJycXObaGVFRURoyZIiys0ofrw+YIW77H/p84htyOBxmRwEAAG6u0nvgx44dW/L74OBgrVq1Sl9//XWFvzchIUGfffZZyWNr1qxRQECACgsLNWrUKPXv319Lly6VpL/s2T97b3pZa5LUt29fhYSEcCdWuIxNa39T/YaN9NiTz5f5qRUqz6p37bMq5m0tzNt6rDrziu59l6pQ4Hv27Fnq6xYtWpR5EusZEyZM0IIFC/TLL7/I39+/5PGAgABJko+Pj4YOHaqgoCBJUqtWrZSSkqKioiJ5e3vLMAzZ7XYFBATI39//vGuAK1v2/WzVb3C57uv7mNlRAACAm3LKdeAnTpyoqKgorVixQvXq1St5PDs7W5mZmSVfR0VFKSwsTJLUuHFjhYWFacaMGZKk+fPny2azyWazlbkGuLqoKZO1eO4MFRVxsycAAFB5XkYNn+WZkpKiVq1aKTAwUJdeeqmk0x+NbNq0SYmJierTp4+Ki4tlGIYCAwP14YcflhTx+Ph4DRgwQBkZGapbt66mTZumDh06lLt2RmxsLIfQwGU1btpC9/V9TDfcdpe8fXzMjuOW8vPzJVn341arYd7Wwrytx+ozj42NlXT6EPXy1HiBNxMFHu6gUZNmuq/vY7rx9rsp8pVk9b/srYZ5Wwvzth6rz7wyBd4ph9AAOL/01MP64oOxembAffp58TwVWvwyWgAAoGwUeMBFHE07oq8+GqcXn3xYh1OSy/8GAABgSRR4wMWkHEjUS0MeVfSmdWZHAQAALogCD7ignOwsjR/9nBbM/Jq7CQMAgFIo8ICLMgxD333zqd57fYRyc7ijMAAAOI0CD7i4zetW6uVnHjvncfFFhYU6lHJA0ZvWKXn/XhPSAQAAZ6v0nVgBON+Z4+J79+mnzOPHlHrIriOH7EpPPSKHo1iSdMmll+m1iV+ole1Kk9MCAICaRIEH3EROdpbmfPv5edezTp3QmyOH6I33v1LjZi2cmAwAADgTh9AAHuR4RrreGPmUMo8dNTsKAACoIRR4wMOkHkrRWy89raxTJ82OAgAAagAFHvBABxITNH7UUOXl5pa53dG0Izp5ItM5oQAAQLWgwAMeKj5uu957fYSKCgtLPZ55PEPLvp+t0UMH6sl+vfXeayNUVFRkUkoAAFBZnMQKeLBtWzdo0vgxGjhkhLasX631q35W7LatMhyOkm127YjWjC8/1IB/DzcxKQAAqCgKPODhNqxeoQ2rV5S5zdIFUWrTroO633ynk1IBAICq4hAaAJKkzya+oQOJCWbHAAAA5aDAA5AkFeTna8Krzyvr5AmzowAAgDJQ4AGUSD18UB+NGy1HcbHZUQAAwHlQ4AGUErNlveZO//K86yczj2vdb8t05KDdiakAAMAZnMQK4C/mz/xKgUHtdU3EjSoqKlLCrh3atnWDtm3dqMSEXTIMQ/51LtGQF19Xl2uvNzsuAACWQoEHcE4fj39FIWHh2vHnZuXmZP9lPSc7S++MGaZ/9BuoBx+JVK3atU1ICQCA9XAIDYBzys3J1ubfV56zvJ9twcyv9faooZz8CgCAk1DgAVywbVs36MWnHlFiwm6zowAA4PEo8ACqRfqRQxo9dKBWLvtBjrPu9FpRRUVFNZAKAADPwzHwAKpNYUG+Pn3vdX3z6XsKaN1GVwQGyXZlWwW0bquA1m100cUXyzAMZR7L0IHEPUrad/rXgX17dORwih554ln1+kdfs18GAAAujQIPoNrl5mQrfuc2xe/cVvKYl5eXmjRrqdycbJ3IPHbO7/vm0/eUlnpIj0YO5aRYAADOg0NoADiFYRg6csh+3vJ+xtIFUZr45kgV5Oc5KRkAAO6FAg/A5Wxet1KvvfBvncw8bnYUAABcDgUegEtK2LVDo559nDu+AgDwPyjwAFzWkUN2vfzMgFLH0gMAYHUUeAAu7dTJE3r1+cGaO/0LFRUWmh0HAADTUeABuLzioiLN/fYLvfhkPyXsijU7DgAApqLAA3Ab9qREjXr2MX3z6XvKy801Ow4AAKbgOvAA3IphGFq6IEpb1q/WY0+9oJCw8JK1/Lw8HUzeL/uBfbInJepo2mFdFdJZ10TcoAaXNzYxNQAA1cfLMAzD7BA1JTY2VuPHDFP+qbKvOw3AfXW59npJkv1AotIOH9T5/kpr076Dwq+7SeHX3ajmrWxOTIjqkp+fL0ny8/MzOQmcgXlbj9VnHht7+hDR4ODgcrdlDzwAt7Z1w5oKbbd3907t3b1Ts77+WC0CWqtrj5t142291bRFqxpOCABA9arxY+Dz8vJ07733KigoSKGhobrzzjuVlJRUsp6QkKCIiAgFBQUpPDxccXFxF7wGAGU5mLxfC2Z+rWcG3KdXhkVq9YolHFMPAHAbTjmJNTIyUvHx8YqJiVHv3r0VGRlZsjZ48GBFRkZqz549GjFihAYOHHjBawBQUbt2RGvyO69q8EN36osP3lLCrtjzHoYDAIArcPox8Fu3btVDDz2kvXv3Ki0tTUFBQTp69Ki8vb1lGIaaNWumjRs3yt/fv0prNput5LliY2M1bvRzKsjiduwAKu6KwLa6/5FIdepyrby8vMyOg/8qKCiQJPn6+pqcBM7AvK3H6jPfuXOnfH19K3QMvNMvI/nRRx/p7rvvliTZ7XY1b95c3t6nD8X38vJSQECAkpOTq7x2RlRUlIYMGaLsrGwnv0IA7u5AYoLee+0Fvfnik4qPjTE7DgAApTj1JNaxY8cqISFBn332Wclj/7t36+wPBKq6Jkl9+/ZVSEgIV6EBUGV74rbrzZFPKfSaCPV9/Em1btPe7EiQda9QYVXM23qsOvPKfPLgtAI/YcIELViwQL/88ov8/f0lSa1atVJKSoqKiopKDoWx2+0KCAiQv79/ldYAoLrFbFmvmC3rFXHjbep570Nq1KSZ6tZrUPIpIAAAzuSU//pMnDhRUVFR+uWXX1SvXr2Sxxs3bqywsDDNmDFDAwYM0Pz582Wz2UqOY6/qGgDUhPWrVmj9qhUlX19a9zJdVr+h6tVvqMvqN1Drtu11a6/75F/nEhNTAgA8XY2fxJqSkqJWrVopMDBQl156qaTTH41s2rRJkhQfH68BAwYoIyNDdevW1bRp09ShQ4cLWjuDGzkBcDb/Opfozr8/qF739VXdevXNjuNRrH6TF6th3tZj9ZlX5kZO3IkVAGqAr5+fbu31D/W+v58ub9y00t9fVFSkHdGblJuTrU5drlWdSy6tgZTuxer/cbca5m09Vp85d2IFAJMV5Odr6cIoLV88VzfcepduvONuBbRuU+bhNY7iYsVtj9bvq37WprW/KevUCUlS7dq1dVVIZ3W59npdfW0PNWnW0lkvAwDggtgDDwBO1LBRE7W8IlCtrghUS1ugAmxXqqioSOtX/awNa37VieMZ5f6MllcE6upuPXR1tx5qe1Wwate2xr4Yq++dsxrmbT1Wnzl74AHARWWkpyojPVXbtm6o8s9IOZColAOJWjR7mupccqk6Xt1NYeHXKfSaa1WvfsNqTAsAcEUUeABwY9lZp7Rh9QptWH366jiBba9SWPh1uibiBrVu2547yQKAB6LAA4AHSUzYpcSEXZo/8ys1btpC3a6/Wd2uv1VXBv2NMg8AHoICDwAeKu3IQf0wZ7p+mDNdjZo0U9cet6hbj1vU9qpgyjwAuDEKPABYQHrqYS2ZN0NL5s1Qm/Yd1KffIHXu2p0iDwBuqJbZAQAAzrV3906NH/2cRj75sDat/U0Oh8PsSACASqDAA4BF7d8br/deH6EXBj+k31cul6O42OxIAIAKoMADgMXZkxL14diX9dygB7R47gxlHjtqdiQAQBk4Bh4AIEk6nJKs6V98oJlfTVJoeIRuvL23ru7aQz6+vmZH8ygb1/6qtCOHFBx6jWyBbVWrdm2zIwFwMxR4AEApDkexojeuVfTGtbrk0svU/eY7dP1td8kWGCRvHx+z47ktR3GxZk2ZrB/mfFvyWJ1L66pDx6sVHNpFHUK7qOUVgZxYDKBcFHgAwHllnTqhZYvmaNmiOfKqVUuXN2qiJs1aqmmLlmrSrKWaNG+pJs1a6GL/S+Tr6ye/iy6Sr6+fvH18KKJnOXUyUx++9bK2R28q9Xj2qZPa/PtKbf59pSSpXoOG6n7Tnbql131qEWAzISkAd0CBBwBUiOFwKD31sNJTDys2ZkuZ23p5ecnXz0/1GzTS4GEvq0OnLk5K6XqS9u3RhFdfUNqRg+Vum3ksQ0vmz9SS+TN1VUhn3XrXfera42b5+vo5ISkAd8FJrACAamcYhvLz8nTkkF2vv/BvfffNpyouLjI7ltP9vnK5Rj37WIXK+//atSNak8aN1uCHeuqbT97TgcQE5WRnqaiwUIZh1EBaAO6CPfAAgBplGIYWzPxaO//coqdfelONmzY3O1KNKyoqUtSUyVo8d/oF/6zsUye1dGGUli6MKvW4j6+ffHx85OPrp8vq1dfQl99WyytaX/DzAXB97IEHADhFfNx2jfjX/2n96hVmR6kxJzOPa8GsKRryyD3VUt7LUliQr5zsLJ04nqHk/Xv1yrBB2rs7tkafE4BrYA88AMBpcrKz9MGbL2n7H5s04N/DddHFF5sdqVrs37tbPy2crd9XLldhYYEpGU6dPKHXXvi3Xnhtgjp27mpKBgDOQYEHADjdbz99r7htf2jg0y+qU5duZsepkuLiIm1et0o/ff+ddsfGmB1HkpSfl6txo4bqmZfeVLcet5gdB0AN4RAaAIApjhyy662Xhuj9N1/SsaNpZsepMMMwtGX9Kj0f2VfvvznSZcr7GUWFhXr/zZf0y9KFZkcBUEPYAw8AMNWG1Sv05+bf9WD/wep57z9Vu7br/qcpfuc2zfjyI8Xv3GZ2lDIZDoe+eP8tZZ08oXsfGmB2HADVzHX/lgQAWEZebo6+/ex9rf55iQY9M1LtOnSqlp+bn5enXTuite2PjUo9fFABtisVENhWAbY2atayVYX/sXAwOUlRUyaX3HDJXcz6+mMdOZSiG267S23bB3MnXcBDeBkefDHZ2NhYjR8zTPmnjpkdBQBQCbfffb8eHTxUvn4XlTyWn58vSfLzO/9NjRwOhw4k7tG2rRu1/Y9N2r0zRkWFhefc1sfHVy1tgbqidVtd3qSpvGt7y9vHR7W9vVX7v7/3ru2tPbu267effpDDUVy9L9LJfP381K5DJ3Xo1EXBodcoMOgqeXu77n68iswbnsXqM4+NPX0VqeDg4HK3dd13LgDAsn5ePE/xO7frudFvq3nLK8rd3jAMrf1lqWZ+NUnHjx2t0HMUFhZof8Ju7U/YfaFx3UJBfr52RG/WjujNkqSLLvbX3zpdrWuvv0XXRNwo/zqXmJwQQEVR4AEALulA4h6NfPIRDX7uZV130x3n3S71cIq+/OBtbY/e5MR07i8vN0fRG9cqeuNa+fj4Kqzrdep+0x3q3LV7qU8+ALgeCjwAwGXl5ebow7EvK257tB56/Cn5+v7/j9aLi4v044IozZn2mQr++9E7qqawsECb163U5nUrddHF/uoScYNu7XWf/taxs9nRAJwDl5EEALi8FUvm6/XnB+vIQbskKTFht/7z9ADN+OJDyns1y8vN0bpff9KrwyM1/YsPz3sOAQDzsAceAOAWDiQmaPTQxxXe/Sat+3WZ259U6g4Wz52uuG1/6Nn/vKWmLVpV6WcUFxcpYVes/tz8u9JTD+uqkDB1CO2iZi0C5OXlVc2JAWugwAMA3EZebo7WrPjR7BiWsm9PnEb8u5+eePYl9bilZ4W+52TmccVsWa/ozb9r2x8blX3qZMnaut+WSZIaXN5YHTpdrQ6hp6+K07hp8xrJD3giCjwAAChTXm6OJo0brW1/bNTAISN0sX+dkjXDMJR6KEV743dqb/xOxe/cpsQ9u1TeVaqPHU3T2l9/0tpff5IktWnfQX36DVLnrt3ZMw+UgwIPAAAqZM2KH7Vn53b1eXigDqcka198nPbuiSu1h72q9u7eqfGjn1Prtu3Vp98ghXTuqlq1OFUPOBdu5AQAAFxOK1sb3fvQAF130+2linxxcZEy0tOUdvig0o4cVKMmzdUhtAtl3wNwIydu5AQAANyYPWmvJo0bpe+/m6K2V3VU+pFDSjtyUEfTjqi4uPQJzE2atdBNd/5dN95+txpc3sikxIDzUOABAIDLsiclyp6UWOY2qYcP6rupn2j2tM/UOfw63dLzXoV1vU61a1Nz4Jn4kw0AADyC4XDoj41r9cfGtarf4HJ1uuZatQ8OVfsOndSs5RWcHAuPQYEHAAAe5/ixo1q1fLFWLV8sSbr0snpq97eOatehk9oHh6rtVSEcNw+35ZQ/uc8884xsNpu8vLxKDtA/w2azqX379goNDVVoaKhmz55dspaQkKCIiAgFBQUpPDxccXFxFVoDAAA426kTmdq6YY1mfjVJo4cO1HOP36+fF89Tfl6e2dGASnPKHvj7779fI0aMUPfu3c+5Pm/evHOecTt48GBFRkZqwIABmjdvngYOHKgNGzaUu3Y2D77IDgAAqKLDB5P11Ufj9N03n+qWnvfqtt7367L6DZzy3IWFBcrLzVV+Xq4cDocaNWlWqcN7DMNQYWGBfH0v7GotxcVFOpp6RIcPJuvIIbsMw1CHTl3U8opAUz6dKCgocPpzupKCggL5+vpWaFunXkbSZrNpyZIlpcr6uR6TpLS0NAUFBeno0aPy9vaWYRhq1qyZNm7cKH9///Ou2Ww2SVJUVJQ+//xz1fcpkI+s/QcCAACUzdvbRxE33aGe9/5TLa8IvOCfZxiGUg4kasv61dq2Zb1OZB5TXm6O8vJyVVxUVGrbRk2aq3O37rq62/UK+lvIOU++LSoq0p64bfpz0+/6c/Pvys/P1fAx78rWpl2FM2VnndSP82fpoD1JR1KSlXrk4F+ySFLdevXVoVOX079Cu+jyxk3/so3D4VBuTpayTp2Uv/8luvSyehXOcT5nCnxFS6yn2blzp3x9fd3nMpL9+vWTw+FQ165d9fbbb6tRo0ay2+1q3ry5vL1PR/Ty8lJAQICSk5NVp06d866dKfB9+/ZVSEgI14EHAADlKioq1JoVS7RmxRIFtG6jq7v1UJdrb9CV7f5W4b3RDodDe3fHavPvK7Vp3UqlHkqp0Pelpx7S8kVztHzRHF1y6WXq3K27rrn2BrW9KkS7dkRr64Y1itmyXtlZp0p935sjn9Jzo95W567nPsLhbPv37tZ7r72otCMHy932ZOZxbVi9QhtWr5AkNWsRoEZNmik765Sysk4q69RJ5WSdKnWUQ8NGTdS6TTu1bttera88/b8NLm9cpROHrXod+Mr8w8X0Ar9mzRoFBASosLBQo0aNUv/+/bV06VJJ+svQz/6DUtYaAABAVSXv36vk/Xu1MGqqLqvfUJ27Xqeru12vjp27qnbt2so8nqHjGek6nnH09P8eO6qM9FTtiN6s48eOXtBzZ506oTUrftSaFT+Wu21+Xq7Gjxmmx596QXfc88B5t1u57Ad9NWm8Cgvyq5Tp8MFkHT6YXOY2GempykhP1dYNa0oeq1uvvv5v4BDdfOffq/S8OD/TC3xAQIAkycfHR0OHDlVQUJAkqVWrVkpJSVFRUVHJYTJ2u10BAQHy9/c/7xoAAEB1OXE8QyuX/aCVy35QrVq15XAUl/9NTmQ4HPp60nilHTmkfoOeLvVpQUFBvqZOnqBfly40JdvJzOP67L03lLQ3Xo/+a1jJkRO4cKZePyk7O1uZmZklX0dFRSksLEyS1LhxY4WFhWnGjBmSpPnz58tms8lms5W5BgAAUBNcrbyfbfHc6Xr/zZEqyD99VZ301MMaM3SQaeX9bMsWzdFbI5/SyROZZkfxGE45ifWpp57SokWLdOTIEV1++eW65JJLtHfvXiUmJqpPnz4qLi6WYRgKDAzUhx9+WFLE4+PjNWDAAGVkZKhu3bqaNm2aOnToUO7aGbGxsRwDDwAALKPtVSHqdd9D+nrSO8o6dcLsOKU0atJML7w64bwn3ubnnz7Ex6rHwJ+51HpFTmJ16lVonI0CDwAA4Dr8LrpI/37+FUXccNtf1ijwFS/w3IIMAAAATpGfl6cP3nxJUVMmcwGSC0CBBwAAgFMtjJqqJfNmmh3DbVHgAQAA4HQzvvpI26M3mR3DLVHgAQAA4HSGw6EP3vyP0g6Xf3MplEaBBwAAgCmyTp3QhNdeUH5entlR3AoFHgAAAKZJ2rdHn018g5NaK4ECDwAAAFP9vnK5ln0/2+wYboMCDwAAANNFTZ2snTFbzY7hFrzNDgAAAAAYDoc+fmeMxk2ersZNm1f4+w7Zk7Ry+WKtX/Wz6tarr6CrOqrtVcEK+ltHNWrSTF5eXjWY2hwUeAAAALiErJMnNH70UN3S6x8KsF2pgNZtVLde/b9sl5OdpfWrV2jV8sXaE7e95PH01MPaFx+nn74//XW9Bg3V9qoQBYd20W139ZG3j4+TXknN8jI8+IyB2NhYjR8zTPmnjpkdBQAAAFVwWb0GCmjdRq1sV6p5qysUv3ObNq37TQX5+ZX6OVcEBumpEa/KdmVQDSW9MLGxsZKk4ODgcrdlDzwAAABc1onMY9rx52bt+HPzBf2cA4l79NJTj6jPw0/o3ocGyNvbfWswJ7ECAADAEoqLizVn2md6+ZkBSt6/1+w4VUaBBwAAgKXsT9itF598WAtmTVFxcZHZcSqNAg8AAADLKS4q0ndTP9GIf/XTsu9n69TJTLMjVRgFHgAAAJZlT9qnKZPf1eCHemri6y8qetM6l98r775H7wMAAADVpKiwUBvX/qqNa39V/YaNdP2tvXTj7XerRYDN7Gh/QYEHAAAAznI8I12LZk/TotnT1KZ9B11/Sy9F3HSH6l5Wz+xokijwAAAAwHnt3b1Te3fv1LTPJiqsa3fdcOtd6ty1u3x8fU3LRIEHAAAAylFcXKyt61dr6/rVqnNpXXXtfrM6h1+n4LBr5F/nEqdm8egCn5+fr8yTWSrIzjU7CgAAADzEiaxcLZwbpYVzo1SrVi21btNO7YJD1T44VM1bXiEvLy9J0skTmTpkT1LKgUQdTE7SQXuSTp3IVKOmzdS4aQs1adZCjZs2V5NmLZV5MkutAwMr9PxehmEYNfkCzZSZmanJkyfr73//u9lR4CRLly5Vr169zI4BJ2Lm1sK8rYV5W4+VZ56fn68rr7xS9erVK3dbjy7wknTPPffohx9+MDsGnIR5Ww8ztxbmbS3M23qYecV4/HXg+/bta3YEOBHzth5mbi3M21qYt/Uw84rx+D3wAAAAgCfx+D3wAAAAgCehwAMAAABuhAIPAAAAuBGPLfAJCQmKiIhQUFCQwsPDFRcXZ3YkVIHNZlP79u0VGhqq0NBQzZ49u2StrBlXdQ3O9cwzz8hms8nLy0uxsbGl1mpivszefGXNnPe7Z8nLy9O9996roKAghYaG6s4771RSUlLJOu9xz1LevHl/VzPDQ910003G1KlTDcMwjLlz5xrdunUzNxCq5IorrjB27NhxzrWyZlzVNTjX6tWrDbvdfs4518R8mb35ypo573fPkpuba/z444+Gw+EwDMMwJk2aZNx2220l67zHPUt58+b9Xb08ssCnpqYal112mVFYWGgYhmE4HA6jSZMmxv79+80Nhko73xu+rBlXdQ3m+d8518R8mb1rqUyBZ+aeYcuWLcaVV15pGAbvcSs4e96Gwfu7unnkITR2u13NmzeXt7e3JMnLy0sBAQFKTk42ORmqol+/fgoJCdGgQYOUnp4uqewZV3UNrqMm5svs3QPvd8/10Ucf6e6775bEe9wKzp73Gby/q49HFnjp9CDPZnC5e7e0Zs0abdu2TdHR0WrYsKH69+9fslbWjKu6BtdRE/Nl9q6N97vnGjt2rBISEvTWW2+VPMZ73HOda968v6uXt9kBakKrVq2UkpKioqIieXt7yzAM2e12BQQEmB0NlXRmZj4+Pho6dKiCgoIklT1jf3//Kq3BddTEfJm96+P97pkmTJigBQsW6JdffpG/v78k3uOe7Fzzlnh/VzeP3APfuHFjhYWFacaMGZKk+fPny2azyWazmRsMlZKdna3MzMySr6OiohQWFiap7BlXdQ2uoybmy+xdG+93zzRx4kRFRUVpxYoVqlevXsnjvMc90/nmzfu7BtTsIfbm2b17t9GtWzejbdu2xtVXX23ExsaaHQmVtG/fPiM0NNQICQkxgoODjXvuuafUySllzbiqa3CuJ5980mjRooVRu3Zto0mTJqVOeKqJ+TJ7851v5rzfPY/dbjckGYGBgUanTp2MTp06GeHh4SXrvMc9S1nz5v1d/bwMw0IHDAEAAABuziMPoQEAAAA8FQUeAAAAcCMUeAAAAMCNUOABAAAAN0KBBwBUyLBhw/Tdd99V+vuGDx+uqKioGkgEANbEVWgAAOU6ePCgbr31VsXFxf3l7oflSUtLU48ePbRr1y7VqsV+IwC4UPxNCgAWtXv3brVs2VKJiYmSpHfffVe9evU65+3Ip0yZovvvv7+kvL/66qvq27evevfurTZt2ujBBx/Un3/+qZtvvlmBgYEaNmxYyfc2btxYrVu31q+//uqcFwYAHo4CDwAW1b59e7377rt68MEHtWrVKn3yySf69ttvz7mHfdWqVYqIiCj12NatWzVz5kzFx8crPj5eI0eO1E8//aQdO3ZoxowZ2rNnT8m2ERERFHgAqCYUeACwsL59+6pz58664447NH36dF1++eXn3C4lJUVNmzYt9dgdd9yhyy67TLVr11bHjh112223yc/PT3Xq1FG7du1K9uxLUtOmTZWSklKjrwUArIICDwAWVlRUpNjYWDVo0EAHDx4873b+/v7Kzc0t9dhFF11U8vvatWv/5euioqKSr/Py8nTxxRdXY3IAsC4KPABY2MiRI9WuXTutWbNGw4cP1969e8+5XceOHbV79+4qP8+uXbvUqVOnKn8/AOD/o8ADgEUtWbJEy5Yt0+TJk9W2bVtNmDBBDzzwgPLy8v6y7f3336+ffvqpSs9jGIZ+/fVX/f3vf7/QyAAAcRlJAEAFOBwOXXPNNVq0aJFatmxZqe9dtmyZZs6cqenTp9dQOgCwFvbAAwDKVatWLX3++edKSkqq9PeeOHFC48ePr/5QAGBR7IEHAAAA3Ah74AEAAAA3QoEHAAAA3AgFHgAAAHAjFHgAAADAjVDgAQAAADdCgQcAAADcCAUeAAAAcCMUeAAAAMCN/D+9JQc0SiIGdwAAAABJRU5ErkJggg==", | |
"text/plain": [ | |
"<Figure size 700x300 with 1 Axes>" | |
] | |
}, | |
"metadata": {}, | |
"output_type": "display_data" | |
}, | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAvAAAAE1CAYAAACBa9TgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAn20lEQVR4nO3df3BU9b3/8dcJIcHFphFMDJHsXWMNfC8JbODb8EOwMrVRVKh+09pLcYRbeoO13F7ndqTXGa8z996KttdxOno71/SXDVIjFKTypZSq36rYFlQupDZEJKKZ3ZBfJBrQGLJJ9nz/WLJJJD/2bPbX2X0+ZhzJfs6e/ex7PkxevPez5ximaZoCAAAAYAtp8Z4AAAAAgNAR4AEAAAAbIcADAAAANkKABwAAAGyEAA8AAADYCAEeAAAAsBECPAAAAGAj6fGeQDR1dXXp1KlTyszMjPdUAAAAgDH19vbq6quvVnZ29oTHJnWAP3XqlDwej6655pqIn9vn80mSMjIyIn7uZEXNrKNm1lEz66hZeKibddTMuo4P+/TO+2ky0qbEeyq2YfoHJMl2NTvt9WhNubRo0aIJj03qAJ+ZmalrrrlGxcXFET93b29v8DUQGmpmHTWzjppZR83CQ92so2bWeZp9OtnKP3isMOI9gXBlSKH+1WAPPAAAAGAjBHgAAADARgjwAAAAgI0Q4AEAAAAbIcADAAAANkKABwAAAGyEAA8AAADYCAEeAAAAsJGYBvh/+7d/k2EYqqurCz7W0NCgZcuWqaioSGVlZaqvr5/0GAAAAJCsYhbgjx49qsOHD8vpdI54fNOmTaqsrNTJkye1ZcsWbdy4cdJjAAAAQLIyTNM0o/0ivb29uv766/XMM89o5cqV2rdvn4qLi9Xe3q6ioiJ1dHQoPT1dpmlq1qxZOnz4sBwOR1hjLpcr+Lp1dXXy+XyaN29exN+Tz+eTJGVkcHvjUFEz66iZddTMOmoWHupmHTWzrql1QP/vDUe8p4EYOO2p01dvlIqLiyc8NiYd+AcffFB33nmnrrrqqhGPe71e5efnKz09XZJkGIacTqc8Hk/YY4Nqamq0efNmHThwIBZvEQAAAIiJ9Gi/wKFDh/Tmm2/qkUceGXXcMIwRPw//QCDcMUlau3atSkpKJEmZmZnWJx6iaJ47WVEz66iZddTMOmoWHupmHTUL3dSpvnhPAQko6gH+1Vdf1YkTJ4Ld96amJt1444362c9+pkWLFqmpqUn9/f3BrTBer1dOp1MOhyOsMQAAACCZRX0Lzb/8y7+oublZjY2Namxs1OzZs/X73/9eq1atUm5urkpLS7V9+3ZJ0u7du+VyueRyucIeAwAAAJJZ1DvwE6mqqtKGDRu0detWZWVlqbq6etJjAAAAQLKKeYBvbGwc8fOcOXN06NChUY8NdwwAAABIVtyJFQAAALARAjwAAABgIwR4AAAAwEYI8AAAAICNEOABAAAAGyHAAwAAADZCgAcAAABshAAPAAAA2AgBHgAAALARAjwAAABgIwR4AAAAwEYI8AAAAICNEOABAAAAGyHAAwAAADZCgAcAAABshAAPAAAA2AgBHgAAALARAjwAAABgIwR4AAAAwEYI8AAAAICNEOABAAAAGyHAAwAAADZCgAcAAABshAAPAAAA2AgBHgAAALARAjwAAABgIwR4AAAAwEYI8AAAAICNEOABAAAAGyHAAwAAADZCgAcAAABshAAPAAAA2AgBHgAAALARAjwAAABgIwR4AAAAwEZiEuDLy8s1f/58ud1urVixQrW1tcExl8uluXPnyu12y+12a8eOHcGxhoYGLVu2TEVFRSorK1N9fX1IYwAAAECySo/Fi+zcuVPZ2dmSpN/85jf6xje+oaNHjwbHd+3apeLi4ouet2nTJlVWVmrDhg3atWuXNm7cqEOHDk04BgAAACSrmAT4wfAuSWfPnlVa2sSN//b2dh09elQvvPCCJKmiokKbN29WY2OjHA7HmGMul2vEeXw+n3p7eyP2XoafF9ZQM+uomXXUzDpqFh7qZh01s66vb0BSRryngQQTkwAvSXfddZdefvllSdKBAwdGjK1bt05+v1+LFy/Www8/rJycHHm9XuXn5ys9PTBFwzDkdDrl8Xg0ffr0MccGA3xNTY2qqqpUXl6uefPmxeptAgAAAJZddWW/Qo3mMQvw27ZtkyRVV1frvvvu0/79+yVJBw8elNPpVF9fnx544AGtX78+OGYYxohzmKYZ/PN4Y5K0du1alZSUSJIyMzMj+2aGiea5kxU1s46aWUfNrKNm4aFu1lGz0E2dyqcWyS4tTVqxUOrrDj2Wx/wqNOvXr9fLL7+szs5OSZLT6ZQkTZ06Vffee69ee+01SVJBQYGamprU398vKRDQvV6vnE7nuGMAAACAHTimSWuul+ZcZe15UQ/w586dU3Nzc/DnPXv2aObMmZoxY4a6u7vV1dUVHKupqVFpaakkKTc3V6Wlpdq+fbskaffu3XK5XHK5XOOOAQAAAInuipnS/7lByp1p/blR30Jz9uxZVVRUqKenR2lpacrJydG+fftkGIba2tpUUVGhgYEBmaapwsLC4FYbSaqqqtKGDRu0detWZWVlqbq6OqQxAAAAIFHNvUq6tlSaMiW85xvmpzePJ5G6ujpJGvUSlZM1eGUb9vGFjppZR82so2bWUbPwUDfrqJl1nmafDvyJq9AkC8MIBPf/VRj483BWcmvMvsQKAAAApKppmdKXlkqzciZ/LgI8AAAAEEWXXyaVL5MudUTmfAR4AAAAIEo+55Su+99Sepj73UdDgAcAAAAizJC0eIFUcs3F+90niwAPAAAARFBmhvTFJdLsK6JzfgI8AAAAECEzPhvY7551afRegwAPAAAARMBVs6XrPy9NjXLCJsADAAAAYZqZLRXkSc5ZgburRnq/+2gI8AAAAECIMqYG9rYX5AX+c1wS+zkQ4AEAAIBxDO+y586Q0tLiOx8CPAAAADBMsMs+Syq4Ij5d9vEQ4AEAAJDyEq3LPh4CPAAAAFLOiC57nuSYFu8ZhY4ADwAAgJRwefZQYE/0Lvt4CPAAAABISnbuso+HAA8AAICkMeOzfl2Z69dVs9Nt3WUfDwEeAAAAtpU5VboyT3LmSbPzpClGX+DxzOSNucn7zgAAAJCUxtvL3tsbt2nFDAEeAAAACe3TXfZk2cseLgI8AAAAEs7llwU67Ha/Ykw0EOABAAAQd5lTA931ArrsEyLAAwAAIC7osoeHAA8AAICYyMy4cF12uuyTQoAHAABA1OQMdtlnSTkzpDQj3jOyPwI8AAAAIibYZZ8lFVwhXUKXPeII8AAAAJgUuuyxRYAHAACAJXTZ44sADwAAgAnlXDZ091O67PFFgAcAAMBFMjOGLvE4O0+6JDPeM8IgAjwAAAAk0WW3CwI8AABAiqLLbk8EeAAAgBSSMyMQ2J150uV02W2JAA8AAJDE6LInHwI8AABAgjIkTZliXvhT6GZ8li57MiPAAwAAJKjcmabuvMWnzEza5hiSFosXKS8v1/z58+V2u7VixQrV1tYGxxoaGrRs2TIVFRWprKxM9fX1kx4DAAAAklVMAvzOnTv11ltvqba2Vt/97nf1jW98Izi2adMmVVZW6uTJk9qyZYs2btw46TEAAAAgWRmmaZqxfMHq6mo98cQTOnLkiNrb21VUVKSOjg6lp6fLNE3NmjVLhw8flsPhCGvM5XIFX6uurk4+n0/z5s2L+Pvw+XySpIyMjIifO1lRM+uomXXUzDpqFh7qZh01s46aWWfXmh0/flwZGRkqLi6e8NiY7YG/66679PLLL0uSDhw4IEnyer3Kz89XenpgGoZhyOl0yuPxaPr06WGNDQb4mpoaVVVVqby8PCoBHgAAIOJ6zyvtw3bJ75ckTUlL18CM3DhPCokmZgF+27ZtkgId+Pvuu0/79++XFAjfww3/QCDcMUlau3atSkpKJCmqX/zgSyXWUTPrqJl11Mw6ahYe6mYdNRvGNKUP2qVWr9TilT7skDSUaabMyFXf8lXULAx2q5mVTwxifhWa9evX6+6771ZnZ6cKCgrU1NSk/v7+4FYYr9crp9Mph8MR1hgAAEBC6+2RWpukFo/U1iT5euM9I9hM1L/Eeu7cOTU3Nwd/3rNnj2bOnKkZM2YoNzdXpaWl2r59uyRp9+7dcrlccrlcYY8BAAAkFNOUOtuk40ekl/ZIe5+W3nhZ8p4ivCMsUe/Anz17VhUVFerp6VFaWppycnK0b9++4BaYqqoqbdiwQVu3blVWVpaqq6uDzw13DAAAIK7O90htF7bF0GVHhMX8KjSxVFdXJ0khfZvXqt7ewF9Eu+2viidqZh01s46aWUfNwkPdrEvqmg3uZW/xBvazf3gmIqf1swfeMruuMyu5lTuxAgAAhIMuO+KEAA8AABCKKHXZAasI8AAAAGOhy44ERIAHAAAYRJcdNkCABwAAqe18TyCst9Jlhz0Q4AEAQGoZ0WX3XLj7KWAfBHgAAJD86LIjiRDgAQBA8jH9Umf7UGiny44kQoAHAADJgS47UgQBHgAA2BNddqQoAjwAALCP859IrU2BL5+2nabLjpREgAcAAIlreJe9xSt10WUHLAX4999/Xz//+c/1hz/8QU1NTbrkkku0YMECVVRUqKKiQunp/HsAAABM0vAue+tpqY8uOzBcyIn77rvv1pEjR/TVr35VP/jBD5SXl6fz58/r7bff1oEDB/TII4/ov//7v7VkyZJozhcAACQbuuyAJSEH+FtvvVVPPvnkRY+XlJTojjvuUGdnp06dOhXRyQEAgCR1/pOhwN5Glx2wwlKAH8/MmTM1c+bMSU8IAAAkIdMvo7Ndae1N0plmqasz3jMCbMvypvVPPvlEzzzzjN577z319/cHH//hD38Y0YkBAACb+1SXPYMuOxARlgP87bffrvT0dC1atEiZmZnRmBMAALCjwb3sLZ5AcKfLDkSF5QDf1NSk48ePR2MuAADAbkZ02ZukPl+8ZwQkPcsBvqSkRC0tLZo1a1Y05gMAABIZXXYg7iwH+H/913/V4sWL5Xa7NW3atODjO3fujOjEAABAgui50GVvpcsOJALLAX79+vVas2aNFi5cqClTpkRjTgAAIJ5Mv9TRNhTa6bIDCcVygPf5fPqv//qvaMwFAADEC112wDYsB/hrr71Wf/3rX1VSUhKN+QAAgFjw+6VOuuyAHVkO8IcPH9YvfvELzZkzZ8Qe+DfeeCOiEwMAABHW84nUeuHLp22n6bIDNmU5wP/oRz+KwjQAAEDE0WUHkpLlAP+FL3whGvMAAACR0NM9bC87XXYgGYUc4O+55x7df//9KigouGjMNE0999xz6u3t1de//vWIThAAAIxjeJe9xSudpcsOJLuQA/yqVat0yy236LLLLtOyZcuUl5ennp4enThxQq+99ppWrVqlf//3f4/mXAEAgDTUZW/xSu102YFUE3KAX716tVavXq0//vGPeuWVV/T222/L4XDo+uuv16OPPqrLL788mvMEACB10WUHMIzlPfDLly/X8uXLozEXAAAwiC47gDFYDvAAACAK/H6pszUQ2Fu90tkP4j0jAAmKAA8AQLzQZQcQBgI8AACxQpcdQASkWX3C008/re7u7hGP7du3L2ITAgAgqfR0S++dkP78gvR8tfTKPumdvxDeAYTNcoCvrKzUddddp7a2tuBjDz744JjHnz9/XrfddpuKiorkdrt10003qbGxMTjucrk0d+5cud1uud1u7dixIzjW0NCgZcuWqaioSGVlZaqvrw9pDACAuPH7pTPN0luvSy/skvb9Svqfg9LpRqm/L96zA5AELAf4uXPn6p//+Z+1YsUKnThxQlLgRk7jqays1DvvvKPa2lrdeuutqqysHDG+a9cu1dbWqra2Vl/72teCj2/atEmVlZU6efKktmzZoo0bN4Y0BgBATNFlBxBDlvfAG4ahdevWKT8/XzfffLOeeuopGYYx5vHTpk3TzTffHPx5yZIl+tGPfjTh67S3t+vo0aN64YUXJEkVFRXavHmzGhsb5XA4xhxzuVwjzuPz+dTb22v1bU7I5+OLRlZRM+uomXXUzDpqFga/XwNtpzWlo0X+jlalffRhvGeEJOU3Tf6OWmTXevl8PmVkZIR0rOUAP9htX7lypfbu3avbb79dHR0dIT//8ccf1+rVq0c8tm7dOvn9fi1evFgPP/ywcnJy5PV6lZ+fr/T0wBQNw5DT6ZTH49H06dPHHBsM8DU1NaqqqlJ5ebnmzZtn9W0CADBST7fS2k8rre200jpaZLAdBkCcWA7wW7duDf65uLhYr7zyin7+85+H/NyGhgY9+eSTwccOHjwop9Opvr4+PfDAA1q/fr32798vSRd19odv1RlvTJLWrl2rkpISSVJmZmZI8wtHNM+drKiZddTMOmpmHTX7FL9f6mgNXC2GK8YgTtIMQxkZGfz9DIPdahZq910KI8CvWrVqxM9XXnnluF9iHfToo4/queee00svvSSHwxF83Ol0SpKmTp2qe++9V0VFRZKkgoICNTU1qb+/X+np6TJNU16vV06nUw6HY8wxAADC9snHQ4G97TRfOgWQkCx/iTUcjz32mGpqavTiiy8qOzs7+Hh3d7e6urqCP9fU1Ki0tFSSlJubq9LSUm3fvl2StHv3brlcLrlcrnHHAAAImd8vtV+4Yszvfy399hnpf17jijEAElrUb+TU1NSk7373uyosLNTKlSslBT7SeP3119XW1qaKigoNDAzINE0VFhZq27ZtwedWVVVpw4YN2rp1q7KyslRdXR3SGAAAY6LLDsDmoh7gZ8+ePeZlJgsLC3Xs2LExnztnzhwdOnTI8hgAAEHBveyewB1Qz3HFGAD2FvUADwBAzNFlB5DECPAAAPujyw4ghRDgAQD2NNhlb/FK7XTZAaQOAjwAwB7osgOAJAI8ACCR0WUHgIsQ4AEAicPvlzpaAoG9lS47AIyGAA8AiK9gl90TuKkSXXYAGBcBHgAQW/6BwF52uuwAEBYCPAAg+j75ONBhb/XSZQeASSLAAwAizz8gnWkdupkSXXYAiBgCPAAgMuiyA0BMEOABAOGhyw4AcUGABwCErvujocBOlx0A4oIADwAY2/Aue4tH+qgr3jMCgJRHgAcAjESXHQASGgEeAFKdf0BGZ5vUeeHa7HTZASChEeABIBUNdtlbPMpob5Yx0B/vGQEAQkSAB4BU4B+QzrRcCO0ju+xG/GYFAAgDAR4AktWwLrvamyW67ACQFAjwAJAsxumyAwCSBwEeAOyMLjsApBwCPADYycCA1NES6LC30mUHgFREgAeAREeXHQAwDAEeABLNiC67R/robLxnBABIIAR4AEgE3eeGtsXQZQcAjIMADwDxMDB4xRjPhb3sdNkBAKEhwANArNBlBwBEAAEeAKKFLjsAIAoI8AAQSXTZAQBRRoAHgMkY3mVv8Uof02UHAEQXAR4ArPr4XKDDTpcdABAHBHgAmMjAgHSm+cLNlOiyAwDiiwAPAKOhyw4ASFAEeACQ6LIDAGyDAA8gddFlBwDYUFq0X+D8+fO67bbbVFRUJLfbrZtuukmNjY3B8YaGBi1btkxFRUUqKytTfX39pMcAYFQD/YGwXvtn6Xc7pN89Kx37k9TiIbwDAGwj6gFekiorK/XOO++otrZWt956qyorK4NjmzZtUmVlpU6ePKktW7Zo48aNkx4DgKCPz0nv1kmv/U56vjrw/4Y6tsgAAGzLME3TjOULHjlyRH/3d3+nd999V+3t7SoqKlJHR4fS09NlmqZmzZqlw4cPy+FwhDXmcrmCr1VXVyefz6d58+ZF/H34fD5JUkZGRsTPnayomXXUzDpfT4+mfNCmjA/aldZ+Wmnd5+I9JQAIW/9lOepZfAO/Byyw6+/O48ePKyMjQ8XFxRMeG/M98I8//rhWr14tSfJ6vcrPz1d6emAahmHI6XTK4/Fo+vTpYY0NBviamhpVVVWpvLw8KgEeQALpPhcI622nldHRKsM/EO8ZAQAQNTEN8Fu3blVDQ4OefPLJ4GOGYYw4ZvgHAuGOSdLatWtVUlIiScrMzJzcxMcRzXMnK2pmHTX7lIH+wN1PWy58AZXtMACSVJphKCMjg98DYbBbzax8YhCzAP/oo4/queee00svvSSHwyFJKigoUFNTk/r7+4NbYbxer5xOpxwOR1hjAJLUx+ekVk8gtJ9pDlz2EQCAFBSTL7E+9thjqqmp0Ysvvqjs7Ozg47m5uSotLdX27dslSbt375bL5ZLL5Qp7DECSGLxizLHhV4z5c+AxwjsAIIVF/UusTU1NKigoUGFhoT7zmc9ICnyk8frrr0uS3nnnHW3YsEGdnZ3KyspSdXV1cM96uGOD6urqJCmkLwNY1dvbG3wvCA01sy7lavbx2aFtMXTZAUD+GbnqW74qdX4PRIBdf3daya0xvwpNLBHgEws1sy7pazbQH7iB0uDNlD7mijEAMBwB3jq7/u60klu5EyuA2KLLDgDApBDgAURX8IoxHrrsAABEAAEeQOR9dHZoW0x7s8R12QEAiBgCPIDJG76XvcUrcfdTAACihgAPIDx02QEAiAsCPIDQ0GUHACAhEOABjG2wy97iCXwRlS47AABxR4AHMIQuOwAACY8AD6S6j7qGAjtddgAAEh4BHkg1/f2BGygNXpe9+6N4zwgAAFhAgAdSwUddw+5+SpcdAAA7I8ADyWhwL3uLR2ptYi87AABJhAAPJIuPzkqtHvayAwCQ5AjwgF1xxRgAAFISAR6wE/ayAwCQ8gjwQCLr71daR4vU2coVYwAAgCQCPJB4znUF9rK3epVxpkWG3x/vGQEAgARCgAfirb9faj8d6LB/qstuxHFaAAAgMRHggXg49+FQYD/Tyl52AAAQMgI8EAv9fSOvy/4Je9kBAEB4CPBAtAzby06XHQAARAoBHoiU/r4Le9mbApd6pMsOAACigAAPTMZgl73FK3W0SFwxBgAARBkBHrBicC/74N1P6bIDAIAYI8ADExm8YkyLV+pgLzsAAIgvAjzwaXTZAQBAAiPAAxJddgAAYBsEeKSmwSvGtHi5LjsAALAVAjxSx4guO1eMAQAA9kSAR/Ia0WX3Sp98HO8ZAQAATBoBHsnl3IdDgZ0uOwAASEIEeNhbf5/UdjoQ2OmyAwCAFECAh/2c/WAosHe00mUHAAAphQCPxEeXHQAAIIgAj8RElx0AAGBUBHgkBrrsAAAAIUmLxYt85zvfkcvlkmEYqqurGzHmcrk0d+5cud1uud1u7dixIzjW0NCgZcuWqaioSGVlZaqvrw9pDDZx9gPpnb9Ir+6Tnq+W/vyC9N7bhHcAAIBxxKQD/5WvfEVbtmzR8uXLRx3ftWuXiouLL3p806ZNqqys1IYNG7Rr1y5t3LhRhw4dmnBsOJ/Pp97e3si+oQvnhTW+T7qV3tmmgQ/alNZ+WkZPd7ynBABAQvObJpnDIrvWy+fzKSMjI6RjYxLgr7vuOsvPaW9v19GjR/XCCy9IkioqKrR582Y1NjbK4XCMOeZyuSRJNTU1qqqqUnl5uebNmxex9wJrjLMfKK39tNLam5XxQbsMk73sAAAAk5EQe+DXrVsnv9+vxYsX6+GHH1ZOTo68Xq/y8/OVnh6YomEYcjqd8ng8mj59+phjgwF+7dq1KikpkSRlZmZGbe7RPLct9fmk9map1RO4oRJddgAAwpZmGMrIyCBvhMFuNQu1+y4lQIA/ePCgnE6n+vr69MADD2j9+vXav3+/pEAwH840zeCfxxtDjA1eMablwhVj6LIDAABETdwDvNPplCRNnTpV9957r4qKiiRJBQUFampqUn9/v9LT02Waprxer5xOpxwOx5hjiIE+n9R+OhDYW+myAwAAxFJMrkIzlu7ubnV1dQV/rqmpUWlpqSQpNzdXpaWl2r59uyRp9+7dcrlccrlc444hSs5+IJ2olV75v9Lz26Q/vyi9f4LwDgAAEGMx6cB/+9vf1vPPP6/W1lbdcMMNuvTSS/Xuu++qra1NFRUVGhgYkGmaKiws1LZt24LPq6qq0oYNG7R161ZlZWWpuro6pDFEAF12AACAhGSYSbx5fPCa86NdonKyBi9NabcvSIzr7AdSi+fC3U/b2MsOAECc+Wfkqm/5quTKG1Fm14xmJbfGfQ884oguOwAAgO0Q4FMNXXYAAABbI8Anuz6f1HY6ENjpsgMAANgeAT4Z0WUHAABIWgT4ZECXHQAAIGUQ4O2qq3MosNNlBwAASBkEeLugyw4AAAAR4BMbXXYAAAB8CgE+kQS77B6ptYkuOwAAAC5CgI83uuwAAACwgAAfa30+qa3pQminyw4AAABrCPCxMNhlb/FKna2SacZ7RgAAALApAnw0DHbZWy5sjTn/SbxnBAAAgCRBgI+UYJfdI3W20WUHAABAVBDgw9XnU9qZZqmjlS47AAAAYoYAb0VXR3BbTEZnmwy67AAAAIgxAvx4+nyBK8UMXpd9WJfdiOO0AAAAkLoI8J82rMvOXnYAAAAkGgJ88IoxF3fZAQAAgESTmgGeLjsAAABsKqkDfG9vrzweT6DL3tkWuGJMR6vUez7eUwMAAJiQ/2yP+i87royMjHhPxTZ8Pp8k2a5mDQ0NcjqdIR1rmGbytp+7urp06tQpZWZmRuX8+/fv18033xyVcycramYdNbOOmllHzcJD3ayjZtZRM+vsWLPe3l5dffXVys7OnvDYpA7w0bZmzRrt3bs33tOwFWpmHTWzjppZR83CQ92so2bWUTPrkr1mafGegJ2tXbs23lOwHWpmHTWzjppZR83CQ92so2bWUTPrkr1mdOABAAAAG6EDDwAAANgIAR4AAACwEQI8AAAAYCMEeIvKy8s1f/58ud1urVixQrW1taMe53K5NHfuXLndbrndbu3YsSO2E00godasoaFBy5YtU1FRkcrKylRfXx/biSaI8+fP67bbblNRUZHcbrduuukmNTY2jnos6yzASs1YZ0O+853vyOVyyTAM1dXVjXkc62xIqDVjnY0Uaj1Sda2FWh/W1Ugpva5MWPLhhx8G/7xnzx6ztLR01OP+5m/+xvzrX/8ao1kltlBrtnLlSvOpp54yTdM0f/3rX5tLliyJwewST09Pj/nb3/7W9Pv9pmma5hNPPGF+6UtfGvVY1lmAlZqxzoa8+uqrptfrnXAdsc6GhFoz1tlIodYjVddaqPVhXY2UyuuKAD8Jv/zlL81FixaNOpaMiyUSxqpZW1ub+dnPftbs6+szTdM0/X6/ecUVV5jvv/9+jGeYeN58803z6quvHnWMdTa6sWrGOhsdAd668WrCOhvJSj1Sca2FWh/W1Uipvq7YQhOGu+66SwUFBXrggQdUXV095nHr1q1TSUmJvvnNb+rMmTMxnGHimahmXq9X+fn5Sk9PlyQZhiGn0ymPxxPrqSacxx9/XKtXrx5znHV2sbFqxjoLH+ssdKyzkazWI9XWWqj1YV2NlOrrigA/zIoVK3T55ZeP+p/X6w0et23bNnm9Xn3/+9/XfffdN+q5Dh48qL/85S86evSoZs6cqfXr18fqbcRUJGtmGMaIn80kvUVBqDWTpK1bt6qhoUEPPfTQqOdinVmvGevs4ppNhHVmvWapss6k0OoWaj1SZa19Wqj1SaV1FYqUXldx7P4nhWnTppkdHR3jHtPc3GxeeumlMZpR4hutZm1tbWZWVhYfDQ7zn//5n+aiRYtGfIdgPKyziWvGOhudlY+XWWcBE22hYZ0NCbceqbLWQq0P62qkVF9XdOAtOHfunJqbm4M/79mzRzNnztSMGTNGHNfd3a2urq7gzzU1NSotLY3VNBNKqDXLzc1VaWmptm/fLknavXu3XC6XXC5XLKebMB577DHV1NToxRdfVHZ29qjHsM5GCqVmrDPrWGfWsc5GCrUeqbrWQq0P62qklF9X8f4XhJ14PB7z85//vFlcXGzOnz/f/OIXv2geO3YsOL5q1SrzzTffNE+dOmW63W6zpKTELC4uNtesWZOy/0IOtWamaZonTpwwlyxZYl5zzTXmokWLzLq6ujjNOr68Xq8pySwsLDQXLFhgLliwwCwrKwuOs84uFmrNTJN1Ntw999xjXnnlleaUKVPMK664YsQXf1lnowulZqbJOvu08erBWgutPhMdl4pSeV0ZppniG6gAAAAAG2ELDQAAAGAjBHgAAADARgjwAAAAgI0Q4AEAAAAbIcADgI0ZhqGPP/44Ysc2NjbqJz/5ybjHNDc3q6ysTH6/P+R5Dvftb39bO3bssPSc1tZWLV68WP39/WG9JgAkEwI8ACAolAD//e9/X5s3b1ZamvVfIaZp6sCBA7r55ptDfk5/f7/y8vK0ePHi4DWfASCVEeABwOYeffRRXXvttSoqKlJNTU3w8eeee05z587V0qVL9R//8R8jnnPgwAEtXLhQ8+fP1xe+8AXV19dLku6++27V19fL7XZrzZo1F73W+fPntWPHDn3lK18JPmYYhh5++GGVlZWpsLBQL730ku6//36VlpZq3rx5On78ePDYI0eOaM6cOZo6dary8vLk9XqDY/fff7++973vSZJcLpceeughrVy5Mnjb869//ev66U9/GoGKAYC9EeABwOYMw9Cf/vQnHThwQP/4j/8or9er9vZ2/cM//IOef/55HTp0SJmZmcHj29vbdeedd6q6ulpvvfWWKisrdccdd0iSnnzySf3t3/6tamtrtXfv3ote680339TnPvc5ORyOEY9nZWXpjTfe0A9+8AN9+ctf1vLly3Xs2DGtX79eDz30UPC4PXv26LbbbtO0adO0ceNGVVVVSZJ6e3v11FNP6Vvf+lbwWI/Hoz/84Q/61a9+JUlatGiRjh07pu7u7sgVDwBsiAAPADb3zW9+U5JUWFio5cuX67XXXtPhw4e1cOFCzZkzR5JUWVkZPP7111+X2+1WSUmJJGndunVqampSS0vLhK/V1NSkvLy8ix7/2te+JklauHCh0tLSdMstt0gKhO733nsveNzevXuDnf177rlHv/zlL+Xz+fTss89q8eLFI26D/vd///cyDCP489SpU5WdnR3SPAEgmRHgASDJGIah8W6ybZrmiGA8/HkTcTgc6unpuejxadOmSZKmTJkyots/ZcqU4BdPT548qezs7OA/AK688kqtWLFCu3bt0o9//GNt3rx5xDkvvfTSi17n/PnzuuSSSyacJwAkMwI8ANjcL37xC0mBL6D+8Y9/1PLly7V06VIdO3ZMJ0+elCT97Gc/Cx6/dOlS1dbW6u2335YkPfvss5o9e7by8vKUlZWls2fPjvlaCxYs0IkTJ8Ka5+D2meH+6Z/+Sd/73vd07tw53XDDDeM+v62tTenp6crPzw/r9QEgWaTHewIAgMnJzMzUtddeqzNnzuiJJ55QQUGBJOknP/mJVq9erZkzZ4740mlOTo6efvpprVu3TgMDA8rOztbOnTslSfPnz9ecOXNUXFyswsLCi/bBu1wu5eTk6Pjx45o3b56lef7mN7/R008/PeKxJUuWKDs7W5WVlRN+AnDgwAHdfvvtIX1SAADJzDDH+5wVAIBP2blzp1599VX9+Mc/Dvk5LS0tuvHGG/XWW2+NeNzr9aqsrEwnT57UZz7zmXHPsWLFCv30pz/V3Llzw5o3ACQLttAAACy54447NHfuXEs3cpo1a9ZF4f3BBx/U0qVL9cgjj0wY3tva2vStb32L8A4AogMPAAAA2AodeAAAAMBGCPAAAACAjRDgAQAAABshwAMAAAA2QoAHAAAAbIQADwAAANgIAR4AAACwEQI8AAAAYCP/H+/eJ452tU+cAAAAAElFTkSuQmCC", | |
"text/plain": [ | |
"<Figure size 700x300 with 1 Axes>" | |
] | |
}, | |
"metadata": {}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"#Inspired by Aletch glacier\n", | |
"vol, area, length = 15e9, 81.7e6, 23e3\n", | |
"width = area/length\n", | |
"minelev,maxelev = 1650, 4160\n", | |
"\n", | |
"x = np.linspace(0,length*1.2,100) #x-location of the grid nodes\n", | |
"dx = x[1]-x[0]\n", | |
"\n", | |
"zbed = maxelev - x*(maxelev-minelev)/length - np.sin((x/length)**.8*np.pi)*900 + np.cumsum(np.random.randn(len(x)))*4000/len(x)\n", | |
"def bdot(z, z_ela=3000, P=0.8):\n", | |
" return np.minimum((z-z_ela)*0.002,P)\n", | |
"\n", | |
"\n", | |
"plt.figure(dpi=100)\n", | |
"plt.fill_between(x,zbed,np.min(zbed),fc='#554433')\n", | |
"plt.ylabel('z (m)')\n", | |
"plt.xlabel('x (m)')\n", | |
"\n", | |
"plt.figure(dpi=100)\n", | |
"plt.fill_betweenx(zbed,bdot(zbed),0*zbed,where=bdot(zbed)>=0,fc='#99aaff')\n", | |
"plt.fill_betweenx(zbed,bdot(zbed),0*zbed,where=bdot(zbed)<=0,fc='#ffaa99')\n", | |
"plt.ylabel('z (m)')\n", | |
"plt.xlabel('bdot (m/yr)')\n" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [] | |
}, | |
{ | |
"cell_type": "markdown", | |
"metadata": {}, | |
"source": [ | |
"# Run a number of timesteps of the model. \n", | |
"... and visualize the results while running the model." | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 32, | |
"metadata": {}, | |
"outputs": [ | |
{ | |
"data": { | |
"image/png": "iVBORw0KGgoAAAANSUhEUgAABGcAAAHiCAYAAABFkqGkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAABcSAAAXEgFnn9JSAACBSUlEQVR4nOzdd3xUVfrH8W8qCYQAIbTQEZBOkI7SZBEQLLCigBXWVRFULFjoLIjsKmsX7CAlCAIqLfTeO4Tea0glvc/c3x/5MUtMgJlkkpkkn/frNS/Ivc8995mZM5fMw7nnuBiGYQgAAAAAAAAO4eroBAAAAAAAAIozijMAAAAAAAAORHEGAAAAAADAgSjOAAAAAAAAOBDFGQAAAAAAAAeiOAMAAAAAAOBAFGcAAAAAAAAciOIMAAAAAACAA1GcAQAAAAAAcCCKMwAAAAAAAA5EcQYAAAAAAMCBKM4AAAAAAAA4EMUZAAAAAAAAB6I4AwAAAAAA4EAUZwAAAAAAAByI4gwAAAAAAIADUZwBANzWhAkT5OLikqtjZ86cKRcXF124cMG+STmh9evXa8iQIWrQoIFKlSqlqlWr6rHHHtO+ffuyxSYkJGjEiBEKCAiQl5eXAgMDNX/+/BzbtTbWljZzUqtWLU2YMMHqeEn64osv5OLioiZNmuS4/8KFC3JxcbHqcfLkSZvOfScbN2687Xl27txpVRu2vp6HDx/W4MGDVbt2bXl5ecnHx0f33Xef/vOf/yg6Ovqu58vN628PefmMHjx4UL1791aNGjXk7e0tPz8/tW/fXnPmzLG6jd27d6tHjx4qXbq0fHx81LVrV23bti1b3M3rUGRkpM152sIez8kR76U9X5/t27drwoQJiomJyXti+cja98qWazMAOJq7oxMAAKCwmz59uqKiovTGG2+oUaNGioiI0LRp09SuXTutWrVKDz74oCW2X79+2rNnj6ZOnar69etr3rx5GjhwoMxmswYNGpSlXWtjbWnTXn766Se5uLjo6NGj2rVrl9q2bZtlv5+fn3bs2GH5OSEhQd27d9fjjz+u9957L0ts/fr17Z7flClT1LVr1yzbbldI+itbXs/vv/9er776qu69916NHDlSjRo1Unp6uvbu3asZM2Zox44dWrJkid2el7OIiYlR9erVNXDgQFWtWlWJiYmaO3eunn32WV24cEFjxoy54/F79uxRp06d1KZNG82ePVuGYeg///mPunXrpg0bNqh9+/YF9Ez+J6/PqSjYvn27Jk6cqBdeeEFly5Z1dDq3Ze17Zcu1GQAczgAA4DbGjx9v5Pafip9//tmQZJw/f96+SeWzxMREm48JCwvLti0+Pt6oVKmS0a1bN8u25cuXG5KMefPmZYnt3r27ERAQYGRkZNgca0ubt1OzZk1j/Pjxd427ac+ePYYk49133zU8PT2Nf/7zn3c9ZuvWrYYk48svv7T6PH/VuXNn4/nnn79jzIYNGwxJxsKFC3N1Dltez+3btxtubm5Gz549jZSUlGxtpaamGn/88cddz2nr628v+fEZbdu2rVG9evW7xvXo0cOoVKlSls9bXFyc4e/vb3To0CFL7M3rUEREhN3ytIW1z8kwHPNe2vP1+fjjjwvldfumv75X1l6bAcAZcFsTADixm8PVDx8+rP79+6tMmTLy8/PTW2+9pYyMDJ08eVI9e/ZU6dKlVatWLf3nP//J1sbWrVvVrVs3lS5dWiVLllSHDh20fPnybHHLly9XYGCgSpQoodq1a+uTTz65bV6nT5/WoEGDVLFiRZUoUUINGzbU119/bfPz27Jli1xcXBQUFJRt3y+//CIXFxft2bPH5vOeOXNGgwcPVr169VSyZElVrVpVjzzyiI4cOZIl7ubru3//fj3xxBMqV66c7rnnHpufR8WKFbNt8/HxUaNGjXT58mXLtiVLlsjHx0f9+/fPEjt48GBdu3ZNu3btsjnWljbt5ccff5Sbm5vefPNN9enTR/Pnz1dSUtIdj9m/f78k6b777rN7PvZky+s5ZcoUubi46LvvvlOJEiWyteXp6alHH33ULnn9/vvvcnFx0bp167Ltmz59uuU6cZO1n3t78/f3l7v73Qdmb9u2TV26dFHJkiUt20qXLq1OnTpp+/btCg0NvePxJ06cUJ06ddS2bVuFh4fb5VqZ1+dkrYiICL300kuqXr26SpQooQoVKuj+++/X2rVrLTEvvPCCatWqle3YO91qevnyZfXr10++vr4qU6aMnnnmGUVERFh93gkTJmjkyJGSpNq1a1tuB9y4caMk66+/N3M8cOBAnvLJjb++V9ZemwHAGVCcAYBC4Mknn1Tz5s21aNEi/fOf/9Snn36qN998U48//rh69+6tJUuW6MEHH9R7772nxYsXW47btGmTHnzwQcXGxurHH39UUFCQSpcurUceeUS//vqrJW7dunV67LHHVLp0ac2fP18ff/yxFixYoJ9//jlbLseOHVPr1q0VEhKiadOmadmyZerdu7def/11TZw40abn1bFjR7Vo0SLHX/C/+uortW7dWq1bt7b5vNeuXVP58uU1depUBQcH6+uvv5a7u7vatm2b4/wm/fr1U926dbVw4ULNmDHDst3FxUVdunSx6TndFBsbq/3796tx48aWbSEhIWrYsGG2L3rNmjWz7Lc11pY27SE5OVlBQUHq1auXKleurMGDBys+Pl4LFy6843EHDhyQq6urmjdvnuP+vLzWORk2bJjc3d3l6+urHj16aOvWrVYdZ+3raTKZtH79erVs2VLVq1e3W96306dPH1WsWDHHz+TMmTN13333WXK09nOfE1vfB7PZrIyMDEVEROibb77RqlWrst22lpO0tLQcC1o3t/21kHqrTZs2qUOHDmrWrJk2bNiQ5Qt4bq+V9nhO1nr22Wf1+++/a9y4cVq9erV++OEH/e1vf1NUVFSe2u3bt6/q1q2r3377TRMmTNDvv/+uHj16KD093arzvvjii3rttdckSYsXL9aOHTu0Y8cO3Xfffbm67uc1n5vu1Cdz817ldG0GAKfg6KE7AIDbuzlcfdq0aVm2BwYGGpKMxYsXW7alp6cbFSpUMPr162fZ1q5dO6NixYpGfHy8ZVtGRobRpEkTo1q1aobZbDYMI3MoeEBAgJGcnGyJi4uLM/z8/LLd1tSjRw+jWrVqRmxsbJbtw4cPN7y8vIzo6GjDMKy/ZeJm3IEDByzbdu/ebUgyZs2aZfN5c5KRkWGkpaUZ9erVM958803L9puv77hx43I8zs3NzXjwwQfvmP/tPP3004a7u7uxd+9ey7Z69eoZPXr0yBZ77do1Q5IxZcoUm2NtafN2bLkV45dffjEkGYsWLTIMI/O1rVy5stGxY8c7HhcYGGg0aNDgtvv/+lqbzWYjPT09y6NTp07Gc889l237rfbv32+88cYbxpIlS4zNmzcbP/30k9GwYUPDzc3NCA4Ovuvzs/b1vH79uiHJGDBgwF3bvBtrX/+33nrL8Pb2NmJiYizbjh07lu12MWs/9zl9Rm3t8y+//LIhyZBkeHp6Gt98841VxwUGBhr169c3TCaTZVt6erpRp06dbLeV3XrbzuzZsw1PT0/j9ddfz3JsXq+V9nhOhmHde+nj42OMGDHijjHPP/+8UbNmzWzbc7rV9Oa2W69thmEYc+fONSQZc+bMsfq8t7utyZbrrz3zMYw798ncvFc5XZsBwBkwcgYACoE+ffpk+blhw4ZycXFRr169LNvc3d1Vt25dXbx4UZKUmJioXbt26YknnpCPj48lzs3NTc8++6yuXLmikydPKjExUXv27FG/fv3k5eVlibv5P+23SklJ0bp169S3b1+VLFlSGRkZlsfDDz+slJSU266Ic2tsRkaGDMOQJA0cOFAVK1bMMnrmyy+/VIUKFfTUU0/l6rwZGRmaMmWKGjVqJE9PT7m7u8vT01OnT5/W8ePHs+X297///bY553Qbyd2MHTtWc+fO1aeffqqWLVtm2Xen1a/+us/aWFvazKsff/xR/v7+lj55sz9t2bJFp0+fzvGYtLQ0HT169I63NP31td60aZM8PDyyPDZv3qxffvkl2/ZbVxtq0aKFPvvsMz3++OPq2LGjBg8erO3bt6tKlSp69913rXqOBfl62mLIkCFKTk7OMvrl559/VokSJSwTFVv7ub8dW/v8qFGjtGfPHi1fvlxDhgzR8OHD73hL5E2vvfaaTp06peHDh+vq1au6fPmyXnnlFcv1y9U1+6+oH374oV544QVNnTpVn3/+eY4xublW2us5WatNmzaaOXOmJk+erJ07d1pGkuTV008/neXnJ598Uu7u7tqwYUOezpvb67698rlTn7T1vbrTtRkAHI3iDAAUAn5+fll+9vT0VMmSJbMUU25uT0lJkSTduHFDhmGoSpUq2doLCAiQJEVFRenGjRsym82qXLlytri/bouKilJGRoa+/PLLbF+QH374YUnKcTnXCxcuZIvftGmTpMzbGF5++WXNmzdPMTExioiI0IIFC/Tiiy9abnGw9bxvvfWWxo4dq8cff1xLly7Vrl27tGfPHjVv3lzJycnZ8svpNcqtiRMnavLkyfrwww81fPjwLPvKly+f460LN5dbvvV9tjbWljbz6syZM9q8ebOefvppeXp6WrYPHjxYUuYKTjkJCQlRenq6TfPNtGzZUnv27MnyuO+++9SnT59s22/259spW7as+vTpo8OHD+f4/t/K2tfT399fJUuW1Pnz561+TnnVuHFjtW7d2nJrk8lk0pw5c/TYY49Z8rL2c28vNWrUUKtWrfTwww9r+vTpeumll/TBBx9kmVskJ0OGDNHUqVM1e/ZsVatWTTVq1NCxY8f0zjvvSJKqVq2a7Zg5c+aoatWqGjBgwG3bzc210l7PyVq//vqrnn/+ef3www9q3769/Pz89Nxzz+n69et5avev12t3d/cs/Tm3583tdT+/8rmVLe/Vna7NAOAMWEobAIqocuXKydXVNceJNa9duyYp8wtmuXLl5OLikuMvxH/dVq5cOcv/wA8bNizH89auXTvbtoCAgCwT+0rSvffea/n70KFDNXXqVP30009KSUlRRkaGXnnllVyfd86cOXruuec0ZcqULDGRkZE5Lg9rr9EQEydO1IQJEzRhwgSNGjUq2/6mTZsqKChIGRkZWeY0uTm/xq1LPVsba0ubefXTTz/JMAy98MILWbY3bNhQbdu21axZszR58mS5ubll2X/gwAFJtk0GXLp0abVq1SrbtvLly2fbbo2bI7Xu9l5b+3q6ubmpW7duWrlypa5cuaJq1arZnFNuDB48WK+++qqOHz+uc+fOKTQ01FIck6z/3OeXNm3aaMaMGTp37pwqVKhwx9j33ntPI0aM0OnTp1W6dGnVrFlTL7/8skqVKpXjqIbg4GA99dRT6tixo9atW6eaNWvm19PIwpbnZA1/f3999tln+uyzz3Tp0iX9+eefev/99xUeHq7g4GBJkpeXl1JTU7Mdm1MR5Kbr169nKWplZGQoKipK5cuXt/q8OcntdT+/8rmT271Xd7s2A4AzYOQMABRRpUqVUtu2bbV48eIsowXMZrPmzJmjatWqqX79+ipVqpTatGmjxYsXZ/mf5Pj4eC1dujRLmyVLllTXrl114MABNWvWTK1atcr2uPmL9608PT2zxZUuXdqyv0qVKurfv7+++eYbzZgxQ4888ohq1KiR6/O6uLhkm2x0+fLlunr1au5f0LuYNGmSJkyYoDFjxmj8+PE5xvTt21cJCQlatGhRlu2zZs1SQECA2rZta3OsLW3mhclk0qxZs9SiRQsFBgZm2z948GCFhoZq5cqV2fbdXKmpRYsWdsnFVjdu3NCyZcsUGBiYbQTFX9nyen7wwQcyDEP//Oc/lZaWlq2t9PT0bJ+hvBo4cKC8vLw0c+ZMzZw5U1WrVtVDDz1k2W/t5z6/bNiwQa6urqpTp45V8SVKlFCTJk1Us2ZNXbp0Sb/++qv++c9/ytvbO1tszZo1tWXLFpUoUUIdO3a87W109mbrc7JFjRo1NHz4cHXv3t3yOZGkWrVqKTw8XGFhYZZtaWlpWrVq1W3bmjt3bpafFyxYoIyMjBwn073deW9eN2/tO7m97tsjH1vl9F5Zc20GAGfAyBkAKMI++ugjde/eXV27dtU777wjT09PffPNNwoJCVFQUJBlFMGkSZPUs2dPde/eXW+//bZMJpP+/e9/q1SpUpbbOW76/PPP9cADD6hjx44aOnSoatWqpfj4eJ05c0ZLly7V+vXrc5XrG2+8Yfnim9OKNLact0+fPpo5c6YaNGigZs2aad++ffr4449tHt3g7u6uzp0733UOjmnTpmncuHHq2bOnevfunW3+hXbt2kmSevXqpe7du2vo0KGKi4tT3bp1FRQUpODgYM2ZMyfLiBNrY21pMy9Wrlypa9euqUuXLvr999+z7b9Z9Pjxxx+zzftx4MAB1a5dO8dRSzdZ+1rfzaBBgyy3Ovj7++v06dOaNm2awsLCNHPmzCyxmzZtUrdu3TRu3DiNGzdOkm2vZ/v27TV9+nS9+uqratmypYYOHarGjRsrPT1dBw4c0HfffacmTZpkm7spL8qWLau+fftq5syZiomJ0TvvvJNt7hVrP/c5sfZ9eOmll+Tr66s2bdqoUqVKioyM1MKFC/Xrr79q5MiRWUYt5PQ6h4SEaNGiRWrVqpVKlCihQ4cOaerUqapXr54mTZp02/NWqVJFmzZtUo8ePdSpUyetWbPGbqPDbHlOuRUbG6uuXbtq0KBBatCggUqXLq09e/YoODhY/fr1s8Q99dRTGjdunAYMGKCRI0cqJSVFX3zxhUwm023bXrx4sdzd3dW9e3cdPXpUY8eOVfPmzfXkk09afd6mTZtKyrzePv/88/Lw8NC9996bq+u+PfKRcu6T1r5X1l6bAcApOHI2YgDAnd26Ssmtnn/+eaNUqVLZ4jt37mw0btw4y7YtW7YYDz74oFGqVCnD29vbaNeunbF06dJsx/75559Gs2bNDE9PT6NGjRrG1KlTc1wZxDAM4/z588aQIUOMqlWrGh4eHkaFChWMDh06GJMnT7bEWLta061q1aplNGzY8Lb7rTmvYRjGjRs3jH/84x9GxYoVjZIlSxoPPPCAsWXLFqNz585G586dLXG3e31vkpQl/nY6d+5sWTEkp8et4uPjjddff92oXLmy4enpaTRr1swICgrKsV1rY21pMyfWrDDz+OOP3/E53ny4u7sb169ftxxnMpmMUqVKGX//+9/v2L41r3Xnzp2N559//o4xH330kREYGGiUKVPGcHNzMypUqGD07dvX2L17d7bYDRs2GJKyPXdbX8+DBw8azz//vFGjRg3D09PTKFWqlNGiRQtj3LhxRnh4+B3zNQzbVssyDMNYvXq15fU+depUjjHWfO5z+oxa2+d/+ukno2PHjoa/v7/h7u5ulC1b1ujcubMxe/bsbLE5vc4nT540OnXqZPj5+Rmenp5G3bp1jTFjxhgJCQnZjs/pcxoTE2Pcf//9hp+fn7Fnzx67XCtteU63c7f3MiUlxXjllVeMZs2aGb6+voa3t7dx7733GuPHjzcSExOzxK5YscIIDAw0vL29jTp16hhfffXVHVdr2rdvn/HII48YPj4+RunSpY2BAwcaYWFhNp/3gw8+MAICAgxXV1dDkrFhwwbDMKy//to7n5z6pLXvlS3XZgBwNBfD+P+bsAEAcKDDhw+refPm+vrrr/Xqq686Op1ipVatWnrhhRc0YcIER6dSLPH6Fx28l9KECRM0ceJERURE5Ov8RgBQ1HBbEwDAoc6ePauLFy9q1KhRqlKlSrbJZgEAAICijgmBAQAONWnSJHXv3l0JCQlauHChSpYs6eiUAAAAgAJFcQYA4FAzZ86UyWRSSEiI7r//fkenAwDIgwkTJsgwDG5pAgAbMecMAAAAAACAAzFyxo6+/fZbR6cAAAAAAAAcIC81AUbO2FGTRo3UtfsjKl3KzdGpAPkqPT1dHh4ejk4DyFf0cxQX9HUUB/RzFAf0c8c7cOCAVq5cmatjWa3JjsqX8Va79k+od+/WKlva0dkA+ScyMpJ7yVHk0c9RXNDXURzQz1Ec0M8db9SoUbk+ltua7Gz1r6O1bnuaklIcnQkAAAAAACgMKM7YWXJCtK6f/lPBW6T0DEdnAwAAAAAAnB3FmXyw6Y+pSktN1Jrtktns6GwAAAAAAIAzoziTT45v/Y+uhBnavM/RmQAAAAAAAGdGcSafHN+3XKXdL+vUBWlviKOzAQAAAAAAzoriTD7avPhdebgb2n9cOn7O0dkAAAAAAABnRHEmH0VdPyMjbqskaet+6eI1BycEAAAAAACcDsWZfLZ2wQcq5ZUuw5DW7ZTCox2dEQAAAAAAcCYUZ/JZRlqKQo/9lPl3kxS8VYpLcHBSAAAAAADAaVCcKQC7134vX6/MITMpqdKqbSyxDQAAAAAAMlGcKSD7Vo2Tq4shSboRJx096+CEAAAAAACAU3D64swPP/wgFxcX+fj4ZNuXkJCgESNGKCAgQF5eXgoMDNT8+fNzbMfaWFvatMXFkzvlbRy3/Lz/aOYoGgAAAAAAULy5OzqBO7l69areeecdBQQEKDY2Ntv+fv36ac+ePZo6darq16+vefPmaeDAgTKbzRo0aFCuYm1p01brfn1HnZ9dptQ0V6WmS3uPSg/cl6cmAQAAAABAIefUxZlXXnlFnTp1kp+fn3777bcs+1asWKE1a9ZYiieS1LVrV128eFEjR47UU089JTc3N5tibWkzNxJiw5Ucukyu5R+VJB0/JzW6R/Irk+smAQAAAABAIee0tzXNmTNHmzZt0jfffJPj/iVLlsjHx0f9+/fPsn3w4MG6du2adu3aZXOsLW3m1obFk1XaO1mSZBjS9oN5bhIAAAAAABRiTlmcCQ8P14gRIzR16lRVq1Ytx5iQkBA1bNhQ7u5ZB/80a9bMst/WWFvavMlkMik9PV3p6ekyrHhuhmHWmV2fWX6+Fi5duGrFgQAAAAAAoEhyytuaXn31Vd17770aOnTobWOioqJUp06dbNv9/Pws+22NtaXNmyZNmqSJEydKkprUr63KPn63zfmmIzsXqfZ9/1BSRkVJ0rYDJnl73JCbU5bKgOySkpIUGRnp6DSAfEU/R3FBX0dxQD9HcUA/L9ycrjizaNEiLV26VAcOHJCLi8sdY++0/6/7rI21pU1JGjt2rEaPHi1J6tap3W2P/avjWz9VzXZTJLkoMdlNodH+Cmxg9eGAQ0VGRsrf39/RaQD5in6O4oK+juKAfo7igH5euDnVWI2EhAQNGzZMr732mgICAhQTE6OYmBilpaVJkmJiYpSYmChJKl++fI4jWaKjoyX9b7SLLbG2tHmTm5ubPDw85OHhoTuXkrI6fXiNfD2vW34+cFxKSrGhAQAAAAAAUCQ4VXEmMjJSYWFhmjZtmsqVK2d5BAUFKTExUeXKldPTTz8tSWratKmOHz+ujIyMLG0cOXJEktSkSRPLNmtjbWnTHvaumiQXl8yZatIzpN1H7No8AAAAAAAoBJyqOFO5cmVt2LAh26NHjx7y8vLShg0bNHnyZElS3759lZCQoEWLFmVpY9asWQoICFDbtm0t26yNtaVNe7h0erd83C5afj51QYqItuspAAAAAACAk3OqOWe8vLzUpUuXbNtnzpwpNze3LPt69eql7t27a+jQoYqLi1PdunUVFBSk4OBgzZkzR25ubjbH2tKmvWxfOlaBvX+RyZx5U9T2g9JjD9r9NAAAAAAAwEk5VXHGVosXL9bo0aM1btw4RUdHq0GDBgoKCtKAAQNyHWtLm/YQdvm4vMzHlahGmT9HSWcuSXVr5MvpAAAAAACAk3ExDMNwdBJFRef2LVXJx5ZpgTOV86+udk8uVnpG5rGlvKWnekruhbp0hqKMmeBRHNDPUVzQ11Ec0M9RHNDPHW/UqFGaMmVKro51qjlniqsbkZflmrzP8nNisnTwpAMTAgAAAAAABYbijJPY+NtoeXqYLT8fOiklJDkwIQAAAAAAUCAozjiJxPgomWO2WH42maRdhx2YEAAAAAAAKBAUZ5zI+kXj5V3CZPn57GXp4jUHJgQAAAAAAPIdxRknkpaSoKTrq7JsW79Liol3UEIAAAAAACDfUZxxMhsXT1Ypr3TLz+kZ0qptUlr6HQ4CAAAAAACFFsUZJ2MypSn6/KIs22LjM0fQsOg5AAAAAABFD8UZJ7Rl6Scq7Z2SZdulUGlPiIMSAgAAAAAA+YbijJO6dvSXbNsOnsicJBgAAAAAABQdFGec1K6138nXOyHb9k17pKiYgs8HAAAAAADkD4ozTuzsnm+ybcswSau3SSmpDkgIAAAAAADYHcUZJ3Zo+wKVKRGRbXt8krR2h2Q2OyApAAAAAABgVxRnnNy2P96Ru1v2ZZquRUg7DzkgIQAAAAAAYFcUZ5zc9UtHpbjNOe4LOSOdvFCw+QAAAAAAAPuiOFMIrJn/QbaltW/auk8KjyrghAAAAAAAgN1QnCkETKY0ndr+saTstzeZzNLq7VJScsHnBQAAAAAA8o7iTCFxdM8fKu12Psd9SSnSyq1SSloBJwUAAAAAAPKM4kwhsm7+GyrhmfMSTVEx0opNUioFGgAAAAAAChWKM4VI3I1QxV9efNv9kTHSis0UaAAAAAAAKEwozhQym/+YKl/v2Nvuj7iRWaBJSy/ApAAAAAAAQK5RnCmE9q0cJVfX7JMD3xRxQ1pOgQYAAAAAgEKB4kwhdPHULnmm7b9jTEQ0I2gAAAAAACgMKM4UUqvmvqlSXneuvIRHSyu3UKABAAAAAMCZUZwppNJTk3Tl8PS7xoVFZRZo0jMKICkAAAAAAGAzijOF2L6Nv8jXM/SucRRoAAAAAABwXhRnCrlNv70pD/fbTw580/VIKZgCDQAAAAAATofiTCEXdf2MMqJWWxUbGikFb6VAAwAAAACAM6E4UwSsXTBWpb2TrIoNjZCWb5JS0/I5KQAAAAAAYBWKM0WAYZi1Y8mr8i5h3ZCY8Ghp6UYpKSVf0wIAAAAAAFagOFNEXL8UokPBr8vL02xVfHSstHSDlGDdgBsAAAAAAJBPKM4UIZdO79bxjSPl6WFdgSY2QfpjvRQTn8+JAQAAAACA26I4U8ScDdmk8zsnWrWCkyQlJmeOoImKyd+8AAAAAABAzijOFEHH9y3XtYMfy83VugJNcmrmHDRhUfmaFgAAAAAAyAHFmSLq0PYFijr5jVytLNCkpUsrNktXw/I5MQAAAAAAkAXFmSJs74aflXDhF7m4WFegSc+QgrdKF67lc2IAAAAAAMCC4kwRtz34S6VdXyTJugKNySyt2S6dvpi/eQEAAAAAgEwUZ4qBTX9MlWJWWR1vGNLG3RRoAAAAAAAoCBRniom1C8bIPWmL1fGGpM37pMgb+ZcTAAAAAACgOFOsBM95U14Z+62ON5kyb3FKScvHpAAAAAAAKOYozhQzy356SaXdzlkdH58krd+ZeasTAAAAAACwP4ozxdDv3w1SGa8oq+OvhEl7QvIxIQAAAAAAijGKM8WQYcrQip+eUmnvZKuPOXhCunA1H5MCAAAAAKCYojhTTCUnxGjz/OdVskSG1cds2C3FxOdjUgAAAAAAFEMUZ4qxqLBzOrzmLXl6mK2KT8+QVm/L/BMAAAAAANgHxZli7sKJ7bq8799ydbVuxt+YeGnj7nxOCgAAAACAYoTiDHRk5yIlXZ4nyboCzfmrmXPQAAAAAACAvKM4A0nS1uWfyj1xs9Xxe0IyV3ECAAAAAAB5Q3EGFsFz35aPyymrYg1DWr9Tik/M56QAAAAAACjiKM4giz++e0ZlSoRbFZuSJq3ZIWWY8jkpAAAAAACKMIozyMIwzFr2w1Py9bZuSEzkDWnr/nxOCgAAAACAIoziDLJJTY7XhrnPqJRXmlXxpy4wQTAAAAAAALlFcQY5uhF5WWt/fkS+Ja5bFb/7iHTuSj4nBQAAAABAEURxBreVGB+lxV/2kWKC5ep692W2N+yWwqMLIDEAAAAAAIoQpyvOHDx4UL1791aNGjXk7e0tPz8/tW/fXnPmzMkSt3HjRrm4uOT42LlzZ7Z2ExISNGLECAUEBMjLy0uBgYGaP39+ruOKk7ULxujCjlEq5ZV+xziTSVq1TUpIKqDEAAAAAAAoAtwdncBfxcTEqHr16ho4cKCqVq2qxMREzZ07V88++6wuXLigMWPGZImfMmWKunbtmmVbkyZNsrXbr18/7dmzR1OnTlX9+vU1b948DRw4UGazWYMGDbI5rrg5dWiNrl04qIee+1lxqZVvG5ecIgVvlR7tKnl6FGCCAAAAAAAUUi6GYdz9fhUn0K5dO127dk2XLl2SlDlypmvXrlq4cKGeeOKJOx67YsUK9e7d21Jouemhhx7S0aNHdenSJbm5uVkddzud27dUJR+XPD5T5/e3JyfL1a+HzObbP9fqlaUeD0iuRf/lKJYiIyPl7+/v6DSAfEU/R3FBX0dxQD9HcUA/d7xRo0ZpypQpuTrW6W5ruh1/f3+5u+duoM+SJUvk4+Oj/v37Z9k+ePBgXbt2Tbt27bIprriz5jany9el7QcKMCkAAAAAAAoppy3OmM1mZWRkKCIiQt98841WrVql9957L1vcsGHD5O7uLl9fX/Xo0UNbt27NFhMSEqKGDRtmK+40a9bMst+WuFuZTCalp6crPT1dhWIIkp2cOrRG62Y9esfVnI6dlY6cLsCkAAAAAAAohJxuzpmbXn31VX377beSJE9PT33xxRd6+eWXLfvLlCmjN954Q126dFH58uV15swZffzxx+rSpYuWL1+uHj16WGKjoqJUp06dbOfw8/Oz7Lcl7laTJk3SxIkTJUlN6tdWZR+/3D7lQichNkKLv+yjvw34RPLtkmPMjoOGZIpTFf87TyaMwiUpKUmRkZGOTgPIV/RzFBf0dRQH9HMUB/Tzws1pizOjRo3Siy++qPDwcC1dulTDhw9XYmKi3nnnHUlSixYt1KJFC0t8x44d1bdvXzVt2lTvvvtuluKMJLm43H7yk1v3WRt309ixYzV69GhJUrdO7ax7ckXM2vnvqO8bWxWf7JXDXhftPV5Gj3SV/MsWdGbIL9zPiuKAfo7igr6O4oB+juKAfl64Oe1tTTVq1FCrVq308MMPa/r06XrppZf0wQcfKCIi4rbHlC1bVn369NHhw4eVnJxs2V6+fPkcR71ER0dL+t/IGGvjbuXm5iYPDw95eHioOM99ezVk5m33pWdIq7ZKicm3DQEAAAAAoNhy2uLMX7Vp00YZGRk6d+7cHeNuLj516yiXpk2b6vjx48rIyMgSe+TIEUn/W3rb2jhkt3vdD/L1jr/t/sTkzAJNGnc3AQAAAACQRaEpzmzYsEGurq45zglz040bN7Rs2TIFBgbKy+t/t9j07dtXCQkJWrRoUZb4WbNmKSAgQG3btrUpDjk7vfPLO+6PjJGWbZKSUwomHwAAAAAACgOnm3PmpZdekq+vr9q0aaNKlSopMjJSCxcu1K+//qqRI0eqQoUKkqRBgwZZbn3y9/fX6dOnNW3aNIWFhWnmzJlZ2uzVq5e6d++uoUOHKi4uTnXr1lVQUJCCg4M1Z84cubm52RSHnB3ZtVj12r6suJTyt42JvCH9sV7q1Ukq41OAyQEAAAAA4KScrjjTvn17/fzzz5o1a5ZiYmLk4+Oj5s2ba/bs2XrmmWcscc2aNdOvv/6qGTNmKCEhQX5+fnrggQc0e/ZstW7dOlu7ixcv1ujRozVu3DhFR0erQYMGCgoK0oABA3IVh5wdWv+R6nT4WMYdZuCJS5T+XC/17ChVKFeAyQEAAAAA4IRcjJuTtCDPOrdvqUo+xXla4Ez9hv+huLSqd43zcJe6t5eqVS6ApGBXzASP4oB+juKCvo7igH6O4oB+7nijRo3SlClTcnVsoZlzBoXH7hUT5epy95pfeoYUvFU6fbEAkgIAAAAAwElRnIHdXTm3XyVdzloVazakDbulQyfzOSkAAAAAAJwUxRnkiy2/j5G7m/V3zO06LO04JHGTHQAAAACguKE4g3wRdf2MPNIO2XTMkVPS+l2SyZxPSQEAAAAA4IQoziDfbPhtlDw9bKu0nL0sBW+R0tLzKSkAAAAAAJwMxRnkm4TYcBlx220+7mq4tGyjlJxq/5wAAAAAAHA2FGeQr9b/NlZeniabj4uMkZZukBKS7J8TAAAAAADOhOIM8lVqcrxSwlbl6tiYeOnPDZl/AgAAAABQVFGcQb7bsHiySnnlbhKZhKTMETSRMfbNCQAAAAAAZ0FxBvnOZEpT1LmFuT4+OTVzDprrkfbLCQAAAAAAZ0FxBgVi67L/qrR3Sq6PT0uXlm+WLoXaMSkAAAAAAJwAxRkUmMtHfsrT8SaTtHpb5nLbAAAAAAAUFRRnUGD2rv9J7omb5Opq5LoNsyGt3ykdO2vHxAAAAAAAcCCKMyhQwXPf1rlt78nHOzXXbRiStu6XDp6wX14AAAAAADgKxRkUuDNH1mv5jO4qaRyVi3I/imb3EWnXYTsmBgAAAACAA1CcgUOkpybpz++fV/jRafIukZHrdg6dlNbvksxmOyYHAAAAAEABojgDhzq0bb42/PKofD2u5LqNM5eklVsyV3QCAAAAAKCwoTgDh0uIDdfirx9X8pXZ8vTI3RCYq+HSnxukhCQ7JwcAAAAAQD6jOAOnsW3F59q75BmV8bqRq+OjY6U/1mf+CQAAAABAYUFxBk4l/NopLfqiu0qk7srV8YnJmQWaq2F2TgwAAAAAgHxCcQZOafmsYfJxPZ2rY9MzMuegOXXRzkkBAAAAAJAPKM7Aaf3x/bMq4xWZq2PNhrRxt7T/mJ2TAgAAAADAzijOwGkZpgwt/+FJlfZOznUbe49Km/ey1DYAAAAAwHlRnIFTS0mK0+b5z6tkiYxct3HivLRqW+btTgAAAAAAOBuKM3B6UWHndGTtSHm4G7lu4/J1aclaKSrGfnkBAAAAAGAPFGdQKJw/vkWhh6fJ1TX3BZqYeOn3dVJI7uYZBgAAAAAgX1CcQaFxcOt8pYQuylMbJrO0/aAUvFVKSbVPXgAAAAAA5AXFGRQqm/+YqhJpu/PczqVQ6bfV0rVwOyQFAAAAAEAeUJxBobN85qsq7X4xz+0kpUjLN0l7QljNCQAAAADgOBRnUCj9+e1A+XrF5LkdQ9KB49LSjVJ8Yp6bAwAAAADAZhRnUCiZTGlaPWugfLztM3FMWJS0aI107rJdmgMAAAAAwGoUZ1BoJcRGaPtv/7BbgSYtXVq7U9q0R0pJs0uTAAAAAADcFcUZFGrhV09ozU+PqoxXpN3aPHlB+nWldPSMZM79yt0AAAAAAFiF4gwKvcT4KC3+8mGVNB+xW5upadK2A9LiNdK1CLs1CwAAAABANhRnUCQYhll//jBYaWG/ydXVfsNdomOlZRultTukhCS7NQsAAAAAgAXFGRQpm/+Yqqv7JsvL075rY5+7Ii0IlvYfkzJMdm0aAAAAAFDMUZxBkXN0zx/a++dglfZOsWu7GSZp71Fp4Srp/BW7Ng0AAAAAKMYozqBIun7pqFZ+/7DKlAize9vxidKaHdLyTdJ1+81DDAAAAAAopijOoMhKSYrToi97yytjf760fzVc+nOD9NvqzJWd0tLz5TQAAAAAgCKO4gyKvGU/vaSUa/PkZseJgm8VHZu5stOcpdLmvVLkjXw5DQAAAACgiKI4g2Jh67L/Kvbcj3JR/hRopMw5aU6clxavlZaszfx7Rka+nQ4AAAAAUERQnEGxsXP1DLknby+Qc0XcyBxFM2dZ5qiaG3EFcloAAAAAQCFEcQbFysrZb8jX43KBnS8tPXM+moWrpKUbpbOXJbN9V/kGAAAAABRyFGdQ7Cz94Wn5eicW+HlDI6R1O6W5y6TdRzJXfQIAAAAAgOIMip301CRtnP+CvEs4ZkKY5FTp4Alp/gopeKt0KVQy8m8qHAAAAACAk6M4g2IpOuy8jm8cJXc3x1VFDGUWZoK3SvNXZhZsklMdlg4AAAAAwEEozqDYOnNkvWLO/iTl4wpO1opPzLzVae6yzImE4xIcnREAAAAAoKBQnEGxtnP1dHmk7HB0GhZmc+YS3L8GZ85PEx3r6IwAAAAAAPmN4gyKvZW/vC5fjyuOTiMLw8hc2em31Zm3PYVFOTojAAAAAEB+oTgDSFr6wyCHrOBkjUuh0h/rM5fivhLm6GwAAAAA5CcfHx+Fh4dLkrp06aL58+dLkqZMmaIRI0ZIkmbOnKmePXta3eaWLVvUsmVLu+f6V7Vq1dLOnTslSa+88or++9//5vs5iwqKM4AyV3DaNH+ww1ZwskZohLRis7RkrXT+Kis8AQAAAEVRQkKCKlasmG37qFGj9Nlnn+WqzY4dO2rfvn15zMw2M2bM0FtvvWXXNs+dO6fWrVurXLly8vPz0+OPP67Q0FDL/sjISPXp00elSpVSvXr1tHr16izH//vf/1bFihVVrlw5vf322zJu+VK1d+9eBQYGqmTJkurYsaPOnz9v19zvhuIM8P+iws7pxKbRDl3ByRoRN6Q126XFazJH1QAAAADA7aSnpzs6BbupUKGCFixYoOjoaF2/fl0NGjSwjCaSpFdffVVVqlRRRESEPvvsMw0YMECRkZGSpBUrVujrr7/Wjh07dPLkSa1bt07fffedJCk1NVV9+/bV66+/rujoaHXt2lWDBg0q0OdGcQa4xenD6xR54ht5epgdncpdRcVmzkezbKMUHu3obAAAAIDi5dKlS+rdu7fKly+vhg0bKjg42LLv7Nmz6tChg0qXLq1+/frpySef1NSpUyVJEyZM0CuvvGKJ3bhxoxo0aGD52cXFRdevX892vr8eZzKZ9Pzzz6t06dJq06aNzp49K0m6cOGCvLy89NVXXykgIEAvvfRSlnPc3H+rW29H6tKli8aPH6/AwED5+PjonXfe0ZkzZ9S6dWuVLVtW77zzjlWvzwsvvJDlOT/33HPq37+/SpcurXbt2unixYuW2CNHjqhTp04qV66cWrZsqb179+bYZunSpVW7dm25uLhIklxdXS0jXBISEvTHH3/oX//6l0qWLKnevXurVatWWrJkiSRp9uzZGjZsmO655x5VrFhR7777rn755RfLe+Dj46MhQ4bIy8tLo0eP1tGjR3X69GnL6/PJJ5+ofv368vX11aeffqrdu3erUaNG8vPz06effmrVa3InTlecOXjwoHr37q0aNWrI29tbfn5+at++vebMmZMtNiEhQSNGjFBAQIC8vLwUGBhouR8vt7G2tImiae+Gn7Vr0SD5ehWOWXivRUi/r5PW7JBi4x2dDQAAAFD0mc1mPfLII3r44YcVFhamn376Sc8884ylqDJo0CB169ZNUVFReu655ywFAnvasGGDunTpoqioKPXq1StL4SYtLU3Hjx/XuXPnNH36dJvbXrJkiVasWKGQkBDNmDFDQ4cO1ZIlSxQSEqIffvhBhw4dsrnNxYsX6/XXX9eNGzdUv359/etf/5IkxcfHq1evXnrzzTcVGRmpsWPHqm/fvkpJSbltW2XLlpW3t7emTZumt99+W5J0+vRplSlTRlWqVLHENW/eXEePHpUkHTt2TE2bNrVqX4kSJXTvvfda9kvSypUrtWfPHq1du1bvv/++Pv74Y23btk0bNmzQqFGjFBERYfNrciunK87ExMSoevXqmjJlilasWKFffvlFtWrV0rPPPqvJkydnie3Xr59mzZql8ePHa+XKlWrdurUGDhyoefPmZWvX2lhb2kTRFXX9jBZ/0UOmyD/k4e7ctznddP6KtGCVtGWflHT76xgAAACAPNq9e7fS09M1bNgwubu7q3379urSpYtWrlypixcvKiQkRGPHjpWnp6cef/xxtW3b1u451K5dW4MHD5anp6dGjRqlU6dOWUajGIahCRMmyMvLK9soGWu8+OKLCggIUK1atdSyZUs99NBDqlatmqpVq6a2bdvq8OHDNrf50EMPqWPHjnJ3d9eAAQMsBZ7ly5erWbNm6tu3r9zc3PT444+rUqVK2rFjx23biomJ0Y0bN/TRRx+pdu3akjIHWvj6+maJ8/X1VUJCQo7777Tvr/sl6Y033lCZMmXUpk0bVa5cWU8++aTKlSun5s2bq0aNGjpx4oTNr8mt3PN0dD7o0qWLunTpkmVbnz59dP78eX333XcaM2aMpMz7xdasWaN58+Zp4MCBkqSuXbvq4sWLGjlypJ566im5ubnZFGtLmygeNiyepMo1lqhD388Vl1zG0enclWFIx89Jpy9KTetLze+VPD0cnRUAAABQtFy6dEmnT59W2bJlLdsyMjLUsmVLhYaGqmLFivL09LTsq169ut1zuLXNEiVKyN/fX6GhoapcubI8PT1VoUKFXLd964TE3t7eWdry9vZWYqLtK93e2mbJkiUthY9Lly5p3bp1WV7L9PT0LBP95sTX11fPPfecmjVrpqtXr8rHx0dxcXFZYuLi4uTj4yNJ2fbfad9f9/81f3u9JrdyupEzt+Pv7y939//VkpYsWSIfHx/1798/S9zgwYN17do17dq1y+ZYW9pE8XH9UogWf95NLnFrnH6y4JsyTNKB49L8FdLhU1KG8y5CBQAAABQ6VatWVdOmTRUTE2N5JCQk6IMPPlCVKlUUHh6utLQ0S/zly5ctfy9VqpSSkpIsP+c0v4w1rly5Yvl7WlqaIiMjVblyZUmyzMmSk1KlSik9PV0Z//8lwWQy5fmWnLyoWrWqevfuneW1TExMtGpCXpPJpNDQUCUlJalevXqKjY3N8noeOnRIjRs3liQ1atRIR44csWpfamqqTp48adlfEJy2OGM2m5WRkaGIiAh98803WrVqld577z3L/pCQEDVs2DBLwUaSmjVrZtlva6wtbd5kMpmUnp6u9PR0GYbBowg/Vge9r5DVr6i0V7zDc7H2kZxqaMdBQ3OXGdp7NPNnR+fEgwcPHjx48ODBg0dhf7Rp00bp6en69ttvlZqaqtTUVG3evFkXL15UjRo11KhRI3344YdKS0vTH3/8od27d1uObd68udatW6fQ0FCFhYXpiy++kCTL/tv9/a+Pc+fOaebMmUpLS9OUKVNUt25d1axZ87bH3dzm7++vKlWqaO7cuUpPT9eUKVOUmpp623Pe7efb7btT3K37e/furb179+qPP/5QRkaGkpKStHLlSsXExGQ7ZtOmTdq3b58yMjIUHR2tt99+W23atFHJkiVVqlQpPfrooxo7dqySkpK0YsUK7d27V48//rgMw9DTTz+t6dOn6+zZswoLC9PHH3+sZ599VoZhqHPnzkpISNDPP/+slJQUTZkyRY0bN1bdunVtek3ywulua7rp1Vdf1bfffitJ8vT01BdffKGXX37Zsj8qKkp16tTJdpyfn59lv62xtrR506RJkzRx4kRJ0r11qquEUcq6J4hCKebAZh052E5tur8uV5+WMuv2FWlnc+qk9Ieroar+qapeKVVeJXJ/8UhJScnVvatAYUI/R3FBX0dxQD9Hfvjggw/09ddf67333pNhGKpfv75GjBihypUra/jw4fr3v/+tjz/+WC1btlSHDh109uxZrV+/Xq6urrrvvvtUr149+fv7q2fPnlqxYoXWr19vaXvr1q3y8/PTjRs3dPToUa1fv17nz5/XjRs3tH79eh0/flzNmzdXUFCQXn31VVWvXl1vvvmm1q9fr+vXr8tsNmdp7+DBg0pKSrJse/XVV/Xuu+/qtdde05NPPil/f3/t3btXSUlJWc4pZX4PPn78uOXniIgInTx5Mkv7N6WkpFjaCQ0Nlbu7e7bcc8pn7Nixmjhxop599lm5ubmpSZMmeuedd7LcViRJu3bt0rfffqvw8HB5eXmpWbNmeuONNyztDBw4UP/+97/l5+cnf39/vffee5b5cby9vdW9e3e1bNlSGRkZ6tWrl+rWrWs59oMPPtDkyZP1yiuvqH79+nrvvfcs+259XpKUlJSk/fv3y2zOXOU3ISFBBw8ezG1XkiS5GHkt7+STS5cuKTw8XOHh4Vq6dKm+++47/fvf/7Ys21W/fn3dc889WrlyZZbjQkNDFRAQoI8++kjvv/++TbG2tHmTyWSyvCHdOrZVJZ/C82UdeVOnUSc17vYfJac6bY3ztlxdpLo1pOYNpHK+d4//q8jISPn7+9s/McCJ0M9RXNDXURzQz+FogwcP1r333pvt+6Q90c8db/To0ZoyZUqujnXab5U1atRQjRo1JEkPP/ywpMxK1vPPP68KFSqofPnyOY5kiY6OlvS/0S6SrI61pc2b3NzcLJMEu7i43PHePhQt549vkYvLe7q3yzSlZxSu992QdPpS5qNmgBTYQKpU3vrj6esoDujnKC7o6ygO6OdwBvndD+nnhZvTzjnzV23atFFGRobOnTsnSWratKmOHz9umcToppuT+DRp0sSyzdpYW9oEJOncsc2KOPaVXFyccgCaVS5ek/5YLy3dIF2PdHQ2AAAAAFD8FJrizIYNG+Tq6mqZE6Zv375KSEjQokWLssTNmjVLAQEBWdaRtzbWljaBm/ZunCUjOtjRaeRZaKT05wZp/S4pMdnR2QAAAABFx8yZM/P1liYUfk53W9NLL70kX19ftWnTRpUqVVJkZKQWLlyoX3/9VSNHjrSsJd6rVy91795dQ4cOVVxcnOrWraugoCAFBwdrzpw5lluNbIm1pU3gVmsXjtVjL9VUoho5OpU8O3NJunBVatFQalZfotsDAAAAQP5yuuJM+/bt9fPPP2vWrFmKiYmRj4+PmjdvrtmzZ+uZZ57JErt48WKNHj1a48aNU3R0tBo0aKCgoCANGDAgW7vWxtrSJnCrP79/QX1fW6G4lMI/CVeGSdoTIp04L7VvLtWq6uiMAAAAAKDoctrVmgqjzu1bslpTMedV0le9/rlc8cnejk7FrqpWkjoE/m9lJ2aCR3FAP0dxQV9HcUA/R3FAP3e8UaNG5Xq1pkIz5wxQGKQkxWnrwiHyLmFydCp2dTVM+m21tP2glJbu6GwAAAAAoGihOAPYWcS10zqxaZTc3YrWoDTDkEJOS/NXShdCSzg6HQAAAAAoMijOAPng9OF1ij71baFeYvt2UlKl/SdKa+lGKTbe0dkAAAAAQOFHcQbIJ7vX/SDFrnV0GvkmNCLzVqcDxyWz2dHZAAAAAEDhRXEGyEdr5n8gH5dTd4zxcDfk45UmX68Y+XpcVmnvZEmFY8SNyZy5qtPitVJ4tKOzAQAAAIDCyemW0gaKmj9+eE6PDPlBLi5uSk0MV2JcqOKiryg67Iwirp1SSlJctmPKV66r5vcPUrmA9krM8FeGyblXAYuOlf5YJzWqK7VpKnlwZQEAAAAAq/EVCshnhilDf37/gk3HRF0/o/WL/iVJ8ihRUs07PKVq93aX2aOOklKd82NrSDp6RrpwVXrgPqlmgKMzAgAAAIDCwTm/5QGwSE9N0t4NP2vvhp8lSXUaddK9Lf8uz3KtlZji6eDssktMllZtk+pUkzq0kEp6OTojAAAAAHBuFGeAQubcsc06d2yz3D291OPpz5ReoqXMZue77encFelKmNSysdT4HsmVGa4AAAAAIEd8XQIKqYy0FC3/+RWd2fy2fL0THZ1OjtLSpR0HM1d1unzd0dkAAAAAgHOiOAMUcueObdaSL7rKNX6d3N2cc5WnmHhp5ZbMR0y8o7MBAAAAAOdCcQYoAgzDrNVB7+nQiiEq43XD0enc1uXr0m+rMkfTpKY5OhsAAAAAcA4UZ4AiJPTiES36orvSw36Tp4fZ0enkyGxIR05Lv66Ujp3N/BkAAAAAijOKM0ARtOmPqdqxsL98Pa85OpXbSkmTtu6XFq+RroU7OhsAAAAAcByKM0ARdSP8ohZ/9ajizs5QSa90R6dzW9Gx0rJN0podUkKSo7MBAAAAgIJHcQYo4nav+0Erv+0mz5SdTjthsCSdvyItCJYOnpDMznlHFgAAAADkC4ozQDGQnpqkFb8M14Flz8nXM9TR6dxWhknafSRz6e0rYY7OBgAAAAAKBsUZoBgJu3xci796RBFHP5GPd6qj07mtmHhpxWZpLbc6AQAAACgG3O3V0Pnz57VixQpt27ZNV69eVXJysvz9/dWoUSM9+OCD6t69uzw8POx1OgB5cGjbfB3e+Zv+9sS/5FH+b0rLcM467bkr0qVQ6b5GUrP6kqtzpgkAAAAAeZLnrzobN25Uz549Va9ePb322mvasmWLEhIS5OHhofPnz2vGjBnq06ePqlWrpnHjxikuLs4eeQPII8OUoTW/jtKOBX9XabdzcpFzzkdz661OV1nVCQAAAEARlKeRM3379tXy5cvVs2dPBQUFqUuXLqpQoUKWGJPJpMOHD2vJkiWaM2eOvv32W82dO1d/+9vf8pQ4APu4EXlZS6Y/qXsDe6hu80dlyFBmncb4/78b+v8NMv7/7x4lqypZ9xRonjHx0vJN0j3Vpfvvk7w8C/T0AAAAAJBv8lScKV26tE6cOKE6dercNsbNzU0tWrRQixYtNGHCBM2ePVtXr17Ny2kB5IOTB1fp5MFVVsf3fuEbpXq2yceMcnb2shQaIXVpLVWrXOCnBwAAAAC7y1Nx5pdffrEp3tXVVc8//3xeTgnASSyf+ap6DPpYJp8uklwK9NxJKdKKLVLjulLbZpK7W4GeHgAAAADsiuk1AeTaqnkjZY5a6rD5ao6ekRavkSJuOOT0AAAAAGAXdlut6aajR4/q4sWLSklJybavX79+9j4dAAdbv+hfeqBPgkpWHSizUbAjaKTMuWj+WJe5olNgQ8m14FMAAAAAgDyxW3Hm7NmzeuKJJ3T48GFJNycO/R8XFxeZTCZ7nQ6AE9m67L9q87d4lav7kkzmgq+OmA1p71Hp8nWpaxvJ16fAUwAAAACAXLNbceall17S9evX9emnn6phw4by9GQpFaA42b32e7VITVClJm8pw+SY4SthUdKiNVK75lLD289TDgAAAABOxW7Fmd27d+v777/XgAED7NUkgELmwJYgNUpJUK3WY5WW4ZgprdIzpC37pEuhUufWLLkNAAAAwPnZ7dtThQoVVKZMGXs1B6CQOrZnqU5teV9enmaH5nHxmvT7OikmzqFpAAAAAMBd2a04M3ToUH3//ff2ag5AIXbmyHodWf2avEs4dp6puATp9/WZc9EAAAAAgLOy221NI0eO1Ntvv62WLVuqV69e8vPzy7LfxcVFb775pr1OB8DJXTy1SynJL+j+x/6t+LQqMuSYeWjS0qXgrZnz0DSt55AUAAAAAOCO7Fac2bVrl2bNmqXo6GgdOHAg236KM0DxE3b5uBZ/9ajuadJZgV3fV2xqBYfkYRjSjoPSjVjpgfskV8dMhwMAAAAAObJbcWb48OHy9/fXTz/9xGpNALI4G7JJZ0M2qcF9D6vxA28qNqWcQ/I4cV6KTZC6t5e8SjgkBQAAAADIxm7FmaNHj2r+/Pl69NFH7dUkgCLmxP4VOrF/hZq176+6bV5VXHLpAs8hNCJzouAeD0jlfAv89AAAAACQjd0G99eoUUOGYdirOQBF2OEdC7X48666ceoLlfZOKvDzxyVmFmguhRb4qQEAAAAgG7sVZ95//3198sknSklJsVeTAIq4fRt/0ZLPOynhwg8q7V2w1470DGnVVunwqQI9LQAAAABkY7fbmvbv36+rV6/qnnvuUdeuXXNcrenzzz+31+kAFCE7V8+QVs9Qh56vqdK9Tyk+2atAzmtI2nlIiryROVGwp0eBnBYAAAAAsrBbcearr76y/H3evHnZ9lOcAXA324O/lIK/LPAizZlLUliU1LWNVNm/QE4JAAAAABZ2u63JbDbf8WEymex1KgBF3PbgL7Xk8weUdGlWgd3uFJ8oLd0o7Q2RzOYCOSUAAAAASLJjcQYA7K2gizSGIe0/Lv25QYpLyPfTAQAAAICkPBZnEhMTC/Q4AMVTQRdpwqOlRWukE+fz/VQAAAAAkLfiTO3atfXpp58qLi7Oqvg9e/bo0Ucf1X//+9+8nBZAMXWzSBNx9BP5lgiTi4uRb+dKz5A275VWb5dS0vLtNAAAAACQt+LMJ598ov/+97+qUqWKBgwYoG+//VZ79uzRxYsXFRYWpuPHj2vZsmUaM2aMmjZtqnbt2qlUqVIaMmSIvfIHUAwd2jZfi7/srUPLnpVXxn55eebfJDEXrkq/rZKuhOXbKQAAAAAUc3larem5555T//79NXPmTM2YMUMLFiyQi4tLlhjDMOTt7a0nnnhCM2fOVMuWLfOUMADcFH71hJb99JLcPb3UvsdwVajTR3HJPnY/T1KKtGKz1KSe1LqJ5GG3de4AAAAAwA5LaXt7e2vo0KEaOnSorl69qu3bt+vatWtKTk6Wv7+/GjRooLZt28rDw8Me+QJANhlpKdqy9BNJn6hes+5q8sDLSjLXlMnsctdjbRFyWjp/RWrbTKpbw65NAwAAACjG7Pr/v1WrVlX//v3t2SQA2OT04TU6fXiNfMtVUbcBnynedI9d209Mltbvko6flTq0kMqXtWvzAAAAAIohltIGUCTF3QjVkulPKfbsdHl5muzefmiktHittHW/lMqEwQAAAADygOIMgCJtz7oftTXoCfmWuG73tg1DOnZW+nVl5p9G/i0eBQAAAKAIozgDoMiLibqsxV/2Udr1hfJwt38FJSUtcwTNkrXS9Ui7Nw8AAACgiKM4A6DY2Pznv7V/6XPy9Y7Nl/YjY6Q/N2TOSZOYnC+nAAAAAFAEUZwBUKyEXT6uxZ93k2v8Orm55s99SGcuSfNXSNsPZi7DDQAAAAB3YrfizKVLl5SQkJDjvvT0dF26dMlepwKAPFsd9J5ObnxDpb3zZ4iLyZy59Pb8FdKuw1JKar6cBgAAAEARYLfiTK1atdSgQQMdOnQo2779+/erdu3a9joVANjFhRPb9edX3eTjejrfzpFhkg6dlIJWSHtDpLT0fDsVAAAAgELKrrc1ubi4qFOnTlqzZk2u21i/fr2GDBmiBg0aqFSpUqpataoee+wx7du3L0vcxo0b5eLikuNj586d2dpNSEjQiBEjFBAQIC8vLwUGBmr+/Pm5jgNQNJhMafp9xkD5el7N1/OkZ0j7j0tByzP/TM/I19MBAAAAKETsWpyZPXu2unbtqj59+mjWrFm5amP69Om6cOGC3njjDa1YsUKff/65wsPD1a5dO61fvz5b/JQpU7Rjx44sjyZNmmSL69evn2bNmqXx48dr5cqVat26tQYOHKh58+blKg5A0fLnjP7y9YrJ9/OkpmeOoAlanjmiJoMiDQAAAFDsuduzsZIlS2rx4sUaNmyYhgwZosuXL2vMmDE2tfH111+rYsWKWbb17NlTdevW1ZQpU/Tggw9m2VevXj21a9fujm2uWLFCa9as0bx58zRw4EBJUteuXXXx4kWNHDlSTz31lNzc3KyOA1D0ZGSkaeVPT6rnP/5QfLJ3vp8vJS1zLpojp6Q2TaV6NSUXl3w/LQAAAAAnZPfVmlxdXTV9+nRNnDhR48aN00svvaQMG/5r+K+FGUny8fFRo0aNdPny5VzltGTJEvn4+Kh///5Ztg8ePFjXrl3Trl27bIoDUDQlJ0RrY9AzKuVVcBPDJKVIG/dIv6+TrkcW2GkBAAAAOJF8W0p7zJgx+umnnzRz5ky9+OKLeWorNjZW+/fvV+PGjbPtGzZsmNzd3eXr66sePXpo69at2WJCQkLUsGFDubtnHSjUrFkzy35b4m5lMpmUnp6u9PR05c+ivAAK0o3wi9q77FV5eZoL9LwRN6Q/N0jrdkoJSQV6agAAAAAOZtfbmv7qhRdeUOXKlbONRLHVsGHDlJiYqNGjR1u2lSlTRm+88Ya6dOmi8uXL68yZM/r444/VpUsXLV++XD169LDERkVFqU6dOtna9fPzs+y3Je5WkyZN0sSJEyVJTerXVmUfvzw8UwDO4Oq5A/La+J7qdfy3Mkz5VsPO0dnL0oWrhurVSFb9Gkly505Kh0lKSlJkJMOZUPTR11Ec0M9RHNDPCze7FWfM5pz/l7lnz546cuSILl68mKt2x44dq7lz5+rLL79Uy5YtLdtbtGihFi1aWH7u2LGj+vbtq6ZNm+rdd9/NUpyRMleSup1b91kbd2t+N4tG3Trdee4bAIXH2ZAN8i33qSo0fktmc8FOBmMyu+jEhZK6HFZSbZpKdWswH40jREZGyt/f39FpAPmOvo7igH6O4oB+XrgVyH8J16pVS507d7b5uIkTJ2ry5Mn68MMPNXz48LvGly1bVn369NHhw4eVnJxs2V6+fPkcR71ER0dL+t/IGGvjbuXm5iYPDw95eHiI705A0XJgS5ASL82Ri4NuWkxMljbslv5YL4VHOyQFAAAAAAWgYMfr22DixImaMGGCJkyYoFGjRll9nGFkfom6dZRL06ZNdfz48WwTEx85ckSSLEtvWxsHoPjYtuJzKXaNQ3MIj86cMHj9LuajAQAAAIoipyzOTJo0SRMmTNCYMWM0fvx4q4+7ceOGli1bpsDAQHl5eVm29+3bVwkJCVq0aFGW+FmzZikgIEBt27a1KQ5A8bLm11Hyyjjg6DR05pK0IFjae1SyYRE8AAAAAE4uXycEzo1p06Zp3Lhx6tmzp3r37q2dO3dm2d+uXea8LoMGDVKNGjXUqlUr+fv76/Tp05o2bZrCwsI0c+bMLMf06tVL3bt319ChQxUXF6e6desqKChIwcHBmjNnjtzc3GyKA1D8LPvpn+o7dKHiTbUdmkeGSdp/TDp5XsxHAwAAABQRTlecWbp0qSQpODhYwcHB2fbfvG2pWbNm+vXXXzVjxgwlJCTIz89PDzzwgGbPnq3WrVtnO27x4sUaPXq0xo0bp+joaDVo0EBBQUEaMGBAruIAFD9LpvdXqy7Pq0bzwYpL9nFoLjfnozl6RmofKFUq79B0AAAAAOSBi3Gz2oE869y+pSr58F/YQHHQsvOzqhn4D4cXaW66p7rUtpnkU9LRmRQdrHiA4oK+juKAfo7igH7ueKNGjdKUKVNydaxTzjkDAM5u36bZWvx5F0UdnyZf73hHp6Ozl6Vfg6W9IVI689EAAAAAhQrFGQDIgwNbgrT4866KOPqJfL3jHJqLySTtP545afC5Kw5NBQAAAIANKM4AgB0c2jZfiz9/UOFH/i1f71iH5pKYLK3dIa3YLMU4flAPAAAAgLugOAMAdnR4x0It/rybUq7OlbubY6f0uhIm/bZa2n2EpbcBAAAAZ0ZxBgDywdbln+rM1vdU0suxVRGzWTp4QlqwSjp/1aGpAAAAALgNijMAkE/OHFmvTXP/Ll+vGEenooQkac12aeUWKS7B0dkAAAAAuBXFGQDIR7FRV7Xk654qaT7i6FQkSZevSwtXSXuPShkmR2cDAAAAQKI4AwD5zjBl6M8fBjvFPDSSZDJL+49Jf6yXDMenAwAAABR7FGcAoIA4yzw0N0XFZI6kAQAAAOBYFGcAoADldh4aN9f8GeJy5FS+NAsAAADABhRnAKCAWeahMUKybHd3M+TjnSrfEuEqaT4ic/RyRR2fppDg57VqRmt5Zey3ey5XwzNH0AAAAABwHHdHJwAAxZFhytCf37+g5h2elCSFXjio8Gt3Hsay7KeX9Pgr85Rgrm/XXI6clrq0tmuTAAAAAGzAyBkAcKBD2xfo0PYFdy3M3PT7jEHy9bxq1xzOXJKSUuzaJAAAAAAbUJwBgELmj+n95esVZbf2zGbp2Bm7NQcAAADARhRnAKCQMZnStOy7v8vXO95ubR47K2WY7NYcAAAAABtQnAGAQigtJUGrZz6p0t7JdmkvJU06fdEuTQEAAACwEcUZACikEmIjtDHoGZXySrdLeyGnJSN/VuwGAAAAcAcUZwCgELsRflG7fv+nvEvYdk+Sp4dZvp6hctH/qjE34qQrYfbOEAAAAMDdUJwBgELu+qUQHV79hkp4mO8Y5+FuyNfjsmLPTteqGQ9o8VePyDNtb5aYw9YtGgUAAADAjtwdnQAAIO8untwpL++xqtNhstIzXCzb3d0MlXQPU+jpYO3d8JPSU5OyHLd85lD9/fVgxab4S5KuhknRsZJfmQJNHwAAACjWKM4AQBFx8uAqeZUqp4Dmb6uke5TCz6/Vtg0/KDkh5o7HrZr1nLo9/4cSUzwkSUdOSZ1bF0DCAAAAACRRnAGAIuXQtvkK2blYJlOa1cckxIbryNr3VL/zNGWYXHTmktSmqeTtlY+JAgAAALBgzhkAKGJsKczcdO7YZsVfnJt5vFk6etbeWQEAAAC4HYozAABJ0rYVn6mkcVSSdPyslGHbAlAAAAAAconiDADAYulP/5Cvd6ySU6UzlxydDQAAAFA8UJwBAFgYpgxtmPuCvEtk6AjLagMAAAAFguIMACCLG5GXdWrrBMUlGLpy3dHZAAAAAEUfxRkAQDYnDwQrLfxPHT5lWH2M2WRSZDjVHAAAAMBWLKUNAMjRhsWT1LdSE92Iu0flfHOOiY4M18G9O3Rwzw4dObBbSYkJevofw/VI/2fl4uJSsAkDAAAAhRTFGQDAbS398XkFVJqnzvfXkCSlpaXqxJEDOrh3pw7t3aHLF7KvuT3n+y906vgRvfrOeJUs5VPQKQMAAACFDsUZAMBtZaSl6Nfpbynu2uMKObBbRw/vU1pq6l2P2711g65ePK+3x3+sajVrF0CmAAAAQOHFnDMAgDuKDLug2d99pgN7tltVmLnp6uULGvXa89qxaW0+ZgcAAAAUfhRnAAD5JiU5SZ9Ofl+zv/tcJlOGo9MBAAAAnBLFGQBAvlu6cLYmvzdMMTeiHJ0KAAAA4HQozgAACsTRQ/v0/qvP6tSxI45OBQAAAHAqFGcAAAUmOjJc49/+p1b9udDutzmlpabYtT0AAACgoLBaEwCgQJkyMvTjl//WL99+puo166jmPfVUs0591apTTzXvqa9SPqXveLzZbFZ46FVdOHtSF86d1sWzp3Th7CndiIrUkOEj9dAjTxTQMwEAAADsg+IMAMAh0tNSde70cZ07fTzLdv+KlVXz/ws1terUU1k/f125eE4Xzp7SxbOndPH8GaUkJ+XY5g9fTFX49Wsa9I/hcnVlcCgAAAAKB4ozAACnEhl+XZHh17Vv55ZcHf/ngl8UGR6qV0dOkKdnCTtnBwAAANgf/60IAChytm9co8nvDVNCXKyjUwEAAADuiuIMAKBIOhFyUGNGDFFY6BVHpwIAAADcEcUZAECRde3yRY15fYjOnDzq6FQAAACA26I4AwAo0mJjojXh7Ze0d/smR6cCAAAA5IjiDACgyEtLTdXHE0cq+I8Fjk4FAAAAyIbiDACgWDDMZv301X/0y7efyWw2OzodAAAAwILiDACgWFn22xxNfOdlXbt8wdGpAAAAAJIozgAAiqHjRw5o5MuDtCToZ2VkZDg6HQAAABRzFGcAAMVSenqagn76WqOGP6dzp447Oh0AAAAUYxRnAADF2oWzpzTqtRc05/svlJaa4uh0AAAAUAxRnAEAFHtms0l/LvhF77w0UEcP7bX5+IyMDBmGkQ+ZAQAAoDhwd3QCAAA4i+vXLmviO6/ob737qle/ZyT5Z9mfkZGh61cv68rFs7p88Zwun8/88/rVS/KrUEmtO3RWm/u7qkHj5nJ1c3PMkwAAAEChQ3EGAIC/WLt8ifZs36z+z/5TCXGxunThrK5cPKdrVy4qIz09x2Mirl/TisVBWrE4SKXLlFWr9p3U5v6uanpfG3l6lijgZwAAAIDChOIMAAA5iL0RpR++mJqrY+NjY7Qh+E9tCP5TXt4lFdi6g9rc30X3tX1AJUv52DlTAAAAFHZON+fM+vXrNWTIEDVo0EClSpVS1apV9dhjj2nfvn3ZYhMSEjRixAgFBATIy8tLgYGBmj9/fo7tWhtrS5sAANxNSnKSdm5eqy8+GqN/PPE3TR0zQru2rL/tCBwAAAAUP043cmb69OmKiorSG2+8oUaNGikiIkLTpk1Tu3bttGrVKj344IOW2H79+mnPnj2aOnWq6tevr3nz5mngwIEym80aNGhQlnatjbWlTQAAbGHKyND+XVu1f9dW+ZYtp07dHlbXno+qeq17HJ0aAAAAHMjFcLLlJcLDw1WxYsUs2xISElS3bl01adJEa9eulSStWLFCvXv3thRPbnrooYd09OhRXbp0SW7/PxmjtbG2tJmTzu1bqpKPi11eBwBA8VGvQRN17fmoOnR5iNueUOAiIyPl7+9/90CgEKOfozignzveqFGjNGXKlFwd63S3Nf21MCNJPj4+atSokS5fvmzZtmTJEvn4+Kh///5ZYgcPHqxr165p165dNsfa0iYAAPZy+kSIvvtsil56qoe++s94HTu0j6W5AQAAihGnK87kJDY2Vvv371fjxo0t20JCQtSwYUO5u2e9M6tZs2aW/bbG2tLmTSaTSenp6UpPTxe/RgMA8iItNVWb1yzXhHde1ojBf9fmNctlNpkcnRYAAADymdPNOZOTYcOGKTExUaNHj7Zsi4qKUp06dbLF+vn5WfbbGmtLmzdNmjRJEydOlCQ1qV9blX38rH5eAADcTujVS/rqP+O1OOhnPfLk82rWsp1cXLh1FvaXlJSkyMhIR6cB5Cv6OYoD+nnh5vTFmbFjx2ru3Ln68ssv1bJlyyz77vRL6l/3WRtrS5s387tZNOrWqd1tjwUAIDeuXb6gb6dNVL2GTTVwyDA1CWzl6JRQxDBHAYoD+jmKA/p54ebUtzVNnDhRkydP1ocffqjhw4dn2Ve+fPkcR7JER0dL+t9oF1tibWnzJjc3N3l4eMjDw0P8fyYAIL+cPn5E/xr5iia/N0xnTx5zdDoAAACwI6ctzkycOFETJkzQhAkTNGrUqGz7mzZtquPHjysjIyPL9iNHjkiSmjRpYnOsLW0CAOAIh/fv0gfDn9MnE0fqysXzjk4HAAAAduCUxZlJkyZpwoQJGjNmjMaPH59jTN++fZWQkKBFixZl2T5r1iwFBASobdu2Nsfa0iYAAI60e+sGvf3SU/rm44k6sHubzp85oRtRkTKZMu5+MAAAAJyK0805M23aNI0bN049e/ZU7969tXPnziz727XLnNelV69e6t69u4YOHaq4uDjVrVtXQUFBCg4O1pw5c+Tm5mY5xtpYW9oEAMDRDLNZG1cv1cbVSy3bXFxcVLpMWZUtV15ly5VXmXJ+Kuvnr7LlyqtK1eq6r+0DcuXfMwAAAKfidMWZpUszf8EMDg5WcHBwtv2G8b8FqxcvXqzRo0dr3Lhxio6OVoMGDRQUFKQBAwZkO87aWFvaBADA2RiGobiYG4qLuaFL589k2185oLoee+p5dfrbw/Lw9HRAhgAAAPgrF+PWagfypHP7lqrkw7TAAADn5+dfUY/0f0bdevWVl7e3o9OBA7G6B4oD+jmKA/q5440aNUpTpkzJ1bFOOecMAADIX9GR4Zo1/b8a9kwfLZr7gxLi4/L1fBFhodkm3AcAAEAmp7utCQAAFJz4uFj9OnOG/lwwWw898oR6/32QypYrb5e2r1+7ou0bV2v7xtW6dP6MSvmUVmDrDmrZrqNatLlfpXxK2+U8AAAAhR3FGQAAoOSkRP3x6yytWDJfXXs8osBW7VWt1j2qWDlArq7WD7SNDL+u7ZvWaMfGNTp76liWfYkJ8dq2YZW2bVglNzc3NWjaQi3bdVSr9p1VOaCavZ8SAABAoUFxBgAAWKSnpWr10t+0eulvkqQSXl6qWr22qtWqoxq17lG1mnVUvdY98q9YWS4umfOsxURHasfmtdq+cY1OHj1k1XlMJpOOHtyrowf36pcZn6pqjdr/X6jppPoNm7KiFAAAKFYozgAAgNtKTUnRudPHde708SzbvbxLqnrNOnL38NCJo4dkmM15Os/VS+d19dJ5/bngF/mULqPmrdqpRZv7FdiqvXzLlstT2wAAAM6O4gwAALBZSnKSTp8IyZe2E+JjLbc/ubi46J57GymwdQfd1+Z+1anfyKbbrAAAAAoDijMAAMBpGYahMyeO6syJo/pt9vfyLVtOzVv+b1SNj28ZR6cIAACQZxRnAABAoREXc0Nb1q3UlnUr5erqpiaBrdSuUze17tBFZcr5OTo9AACAXKE4AwAACiWz2aTD+3fp8P5d+v6LqWrcrKXadeqmNvd3UVk/f0enBwAAYDWKMwAAoNAzzGaFHNyjkIN79OOX/1aDJi3UvlM3tXngQfn5V3B0egAAAHdEcQYAABQphmHo+JH9On5kv37+5hPd27i5Hu47QG0eeJDJhAEAgFOiOAMAAIoswzB0IuSgToQcVLWaddRv0BB16Nxdrm5ujk4NAADAgv8+AgAAxcKVi+f0xUdj9OY/+mvjqqXKyMhwdEoAAACSKM4AAIBiJvTqJX3zyUSNGPx3rV2+WBnp6Y5OCQAAFHMUZwAAQLEUfv2qvvtsil57/nEF/7FAaWmpjk4JAAAUUxRnAABAsRYVEaafvvqPhj/7qBbP+0lREWGOTgkAABQzFGcAAAAkxURHaf7P3+jVp/vow/eHa+v6YKWlpjg6LQAAUAywWhMAAMAtDMPQoX07dWjfTnmXLKX7uz6kLg89onoNm8rFxcXR6cEBDu/bqVkzPlWdeg3VpEUrNW7eSv4VKzs6LQBAEUJxBgAA4DaSkxK1dvkSrV2+RFWq1VCXhx5R5+695edf0dGpoQAYhqGlC2dr7o9fyTCbdfnCWW1as0ySVDmgupq0aKUmga3VqHlLlS1X3sHZAgAKM4ozAAAAVgi9cklBP32t+TOnq3GzlqpRp64qB1RX5YBqqlSlmipUDpC7O79aFRUpycmaPu1f2rFpTY77r1+7rOvXLmvt8iWSpOq16qhxYGs1bdFaga06yMPTsyDTBQAUcvwGAQAAYAPDbFbIwT0KObgny3YXV1f5V6xsKdbc/LNMOT95lvCSp2cJeZYooRIlvORZIvPvbm78KuaMrl+7ok8mvKNL589YfczlC+d0+cI5Bf/+q0r7llHn7n3U7eG+qlqjVv4lCgAoMviNAAAAwA4Ms1kR168p4vo1HdFuq45xc3PLLNyU8FK1GrX0yttjValKtXzOFHdycM8Off7RaCXGx+W6jfi4WC1bNFfLFs1Vw6b36W+9+6ptxwfl6VnCjpkCAIoSVmsCAABwEJPJpOSkRMXeiNLRQ/v07itPa+v6YEenVSwZhqHf58/UR6Nfz1Nh5q+OH9mvL6eO1SsDHtbM6dN05eI5u7UNACg6GDkDAADgJJKTEvXFR2N0eN8uDRk+Ul7eJR2dUrGQkpykbz6eqJ1b1uXbORLiY7VicZBWLA5SgyaB6vbw42oS2FqeJUrIw8NTHp6e3OYGAMUY/wIAAAA4mY2rl+rk0UN6Y/QU1anXwNHpFGnXr17WxxPe1uULBTei5UTIQZ0IOZhtu6urmzw8PeThUUIeHh7y8CwhD09PVawcoOHvTVRp37IFliMAoGBxWxMAAIATCr16SaNff0HLfpsjs9ns6HSKnMSEeC1bNFcfDH+uQAszd2I2m5SakqKE+FjdiI5U+PWrunrpvA7s3qbxb72k6MhwR6cIAMgnFGcAAACclCkjQ798+5mmjhmhmBtRjk6nSLhy8bx++GKqXhn4sH6Z8akSE+IdnZJVrlw8p7Ej/qFrVy46OhUAQD6gOAMAAODkDu7ZrpEvD9KhvTsdnUqhZDaZtHfHZk1+b5jeerG/Vi/9TakpyY5Oy2YRYaEa9+aLOn/mhKNTAQDYGcUZAACAQiD2RpQ+/GC4Zn/3udJSUxydTqFw89alNwb303/GvaXD+3c5OqU8i4u5oQlvv6xjh/Y5OhUAgB1RnAEAAChEli6crbf/+ZQO7N7m6FSc1rXLF7LcuhQWetXRKdlVclKiPvzgNe3dvsnRqQAA7ITiDAAAQCETFnpVH41+Q59MHKnI8OuOTsdpREeGa8a0SXrzxScL7a1L1kpPT9MnE9/VxtXLHJ0KAMAOWEobAACgkNq9dYMO7d2p/s+9pIf7DpS7e/H81S4pMUG/z5+lFUvmKS011dHpFBiz2aRvPp6ghLgY9XniGUenAwDIg+L5LzgAAEARkZqSrDnffa5Nq5fqxdc/UMOmLRydUoFJT0vTqqULtXjuT0qIj3V0Og7zy7efKT4uVgMGvyoXFxdHpwMAyAWKMwAAAEXA5QvnNP6tf6pz9z569qU35Fu2nEPyiAy/rsP7durCudMKqFpDNerUU43adeVT2tdu5zCbzdq2YZXmz5yuiOvX7NZuYbYk6GfFx8Vo4JBhKu1b1tHpAABsRHEGAACgCNm0Zpn27tysQUOGqdvDfeXqmr9TDKYkJ+nooX06vG+XDu/bqauXL+QYV75CJdWsU081atdTzTp1VbNOPVWpVkNubrb9Onp4307N/eFLnT9z0g7ZFy1rly/RuhW/q2ademoc2EqNm7dSo2b3qWQpH0enBgC4C4ozAAAARUxifJy+//wj7di0Vq9/MEll/fzt1rbZZNL5Myd1aN9OHd63UyePHZYpI+Oux0VFhCkqIkz7d221bPPw8FTVGrVUqrSv3N3c5ebhkfmne+bD/f8fbm7ucvfw0KVzZ3TkwG67PZeiyDAMXTh7ShfOntLyRfPk4uqqOnUbqHFgKzUJbKUGTQLl5V3S0WkCAP6C4gwAAEARFXJwj94d+rRe/+BDNQlslae2MtLT9ceCWVqxOEjxcfaZ3yU9PU0Xzp6yS1vImWE26+ypYzp76pj+XPCL3NzcdM+9jdW8VTvd3+UhBVSv5egUAQCiOAMAAFCkxURHadJ7r6r/M/9Uv0FD5OrmZnMbp44d1refTtblC+fyIUMUJJPJpFPHDuvUscNa+Mt3qnVPfXXo8pDu79pDFSpVcXR6AFBsUZwBAAAo4gyzWQt++VbHQw7otfcnqWy58lYdl5SYoKCfvtbqpb/JMIx8zhKOcPMWqHk/fqX6jZqpQ5eH1KHz3+x6KxwA4O7yd4Y4AAAAOI0j+3fr3VcG6eihvXeN3bt9k9568Umt+nMhhZli4tSxw5r5zSd6eeDD+tfIoVq7YokS7HQLGwDgzijOAAAAFCMx0VH617uvatHcH2Q2m7PtvxEVqf/+6z39Z/zbio4Md0CGcDTDbFbIwT367tMP9cqgh7V2+WIKdACQzyjOAAAAFDOG2axfZ87QlA9eU3xsjCTJbDZr7YolevMfT2jnlnWOTRBOIy01Vd99NkWfTnpfCfFxjk4HAIos5pwBAAAopg7v36WL77+q5155U2uXL9HxI/sdnRKc1M4t63Tm5FG9/sFkNWgSmK/nMgxDF8+e0oE921WihJeatGitajXryNWV/1cGUHRRnAEAACjGYmOi9eXUsY5OA4VAZPh1TXj7ZfV/7p/qO2Bwrlb+up2kxAQd3r9LB3dv14Hd23QjOjLL/tJlyqpxs5Zq0qK1GjdvpYDqNeXi4mK38wOAo1GcAQAAAGAVs9mkX2fO0JH9u/Xa+5NUvkKlXLVjGIauXDynA7u36cDubToRclAmk+m28fGxMdq5ZZ3llrtyfv5qHNhKjZu3UuPAVqpUpSrFGgCFGsUZAAAAADY5dni/Rr48SK++M06tOnS+a7zZZNKVS+d15uRRnT4eokN7dygy/Hquz38jOlJb1wdr6/pgSVKN2nXV7+l/qN0DD9p1RA8AFBSKMwAAAABslhAfq/+Mf1s9H3tSz7z0hjw9S0jKHBUTERaqMyeP6uzJYzpz4qjOnT6u1JTkfMvl0vkz+mzyB6pao7b+/vQ/1KFzd4o0AAoVijMAAAAAci34jwU6fuSAWnfoorOnjunsqWOKi7nhkFyuXjqvLz4ao4Wzv1O/QUP0wIM95ebGVx4Azo8rFQAAAIA8uXjutC6eO+3oNCxCr1zS1/+ZoN9mf6++g4aocWBbq481m82KiY6Sd8mS8i5ZKh+zBID/oTgDAAAAoEgKC72qGdMmyc+/op545kV17t5HHp6eSkyIV/j1qwoPvaaw//8z4vpVhV+/pvDr15SenqYSXt7q0KW7uvV6XPUaNmXCYQD5ysUwDMPRSRQVndu3VCUfLtoAAACAMypdpqzMZrMS4+NsOq5azTp6sNdj6vS33vItUzZ/kgPyKDIyUv7+/o5Oo1gbNWqUpkyZkqtjXe2cCwAAAAA4pfjYGJsLM5J05eI5/TLjU708oKc+nfyBDu/bKbPZnA8ZAiiuuK0JAAAAAKxgysjQjk1rtGPTGlWoVEVdezyqTt17q2LlAEenBqCQc7qRM/Hx8Xr33Xf10EMPqUKFCnJxcdGECROyxW3cuFEuLi45Pnbu3JktPiEhQSNGjFBAQIC8vLwUGBio+fPn5zoOAAAAQPEVERaqBb98q+HPPqphzzyiLz4ao1V/LtSFs6dkNpkcnR6AQsbpRs5ERUXpu+++U/PmzfX444/rhx9+uGP8lClT1LVr1yzbmjRpki2uX79+2rNnj6ZOnar69etr3rx5GjhwoMxmswYNGmRzHAAAAABImYWaiLBQbV0fLEnyLllK9Rs21b2Nm+vexs1Vr2ETeXmXdHCWAJyZ0xVnatasqRs3bsjFxUWRkZF3Lc7Uq1dP7dq1u2PMihUrtGbNGkuhRZK6du2qixcvauTIkXrqqafk5uZmdRwAAAAA3E5yUqIO7dupQ/syR/S7urqp1j311fmh3ura41EKNQCycbrbmm7emmRPS5YskY+Pj/r3759l++DBg3Xt2jXt2rXLpjgAAAAAsJbZbNK508f189efaOig3pr341eKjoxwdFoAnIjTFWdsNWzYMLm7u8vX11c9evTQ1q1bs8WEhISoYcOGcnfPOlCoWbNmlv22xN3KZDIpPT1d6enpYk1yAAAAAHeSmBCv3+fP1LBnH9FX/xmvi+dOOzqlYs0wDEVFhOnoob2KDL/u6HRQjDndbU3WKlOmjN544w116dJF5cuX15kzZ/Txxx+rS5cuWr58uXr06GGJjYqKUp06dbK14efnZ9lvS9ytJk2apIkTJ0qSmtSvrco+fnl/cgAAAACKNFNGhjavWa7Na5arQZMW6tb772rUvKXd7yKwRUZGui6cOanEhHilpiT//yNFKSlJSk1J+d+21BSlpaSoWq06atayverUbyhX19xP/5CclKijh/bq4pmT6tV3oEr6lLbjs/rfOcKuXVFY6BWFX7+qsGtXFB56VeHXrygtNdUSV7FKVTVo0kINmrZQ/YbN8iWX/JKUlKTIyEhHp4FcKrTFmRYtWqhFixaWnzt27Ki+ffuqadOmevfdd7MUZyTd8SJ36z5r424aO3asRo8eLUnq1unOc98AAAAAwF+dCDmgEyEHVL3WPerz96f1wIM95eHpWSDnTklO1sG927V76wbt37VVSYkJVh97IuSA1i5bJN+y5XRf2wfUukNnNbuvnUp4ed312LDQK9q3Y4v27dyiY4f3yfT/K1ydOLJfH3z4uSpWqZrr53RTRFiofvhiqs6dOq7YmGirjgkPvarw0KvavGaZXFxddU/9hmraoo2a3tdW9zZqdsf3xTAMpaakKCE+VgnxcUpNSVaVqjXkW7Zcnp+LNSIjI+Xv718g54L9FdriTE7Kli2rPn36aMaMGUpOTpa3t7ckqXz58jmOeomOzvyA3hwZY23crdzc3CyTBDuuxg0AAACgsLt84aymT/uXZk6fpuYt26lV+05q0fZ+lfYta9fzJMTFau/OLdq9db0O7dul9LTUux90B3ExN7Rx1VJtXLVUniVKqFnLdmrdvrNatutoKUyYTSadPhGifTu3aO+Ozbpy8VyObV29fEGjXx+s9yb9V3UbZF+F11oH9+zQFx+NUUJ8bK7bMMxmnTlxVGdOHNWSoJ/lWaKEGja9TwHVaioxMV4J8XFKjI9Twv8/EhPilJGenq2d8hUqqXbde1W7XgPVrttAteveKz//ig4dJQXnU6SKM1JmtVLKOsqladOmCgoKUkZGRpb5ZI4cOSLpf0tvWxsHAAAAAPklOSlRO7es084t6+Ti6qp7GzVTy3Yd1bJdJ1WtUStXX+qjI8O1e9tG7dm2UUcP7ZPZbMqHzKW01FTt3b5Je7dvsuReoVIVHdy7Q/GxMVa1ERsTrQnvvKzX35+sNg90ten8ZrNZi+f+qIWzv7N8N7SXtNRUHdq7Q4f27rDpuKiIMEVFhGnvjs2Wbb5ly6n2Pf8r2LRs31GeniXsmi8KlyJVnLlx44aWLVumwMBAed0ylK5v3776/vvvtWjRIj311FOW7bNmzVJAQIDatm1rUxwAAAAAFATDbNaJkIM6EXJQc3/4UpUCqqllu45q1a6TGjRtIVcXF8XG3NCNqAjdiI7UjagIRUdFKCY6UtFRkYqJirTsc2TutkpLTdW0f72rZ18eod79BllVkEqIi9WXU8fqwJ7tuci2YMXF3Miy3HrdBo31zviP5edf0cGZwVGcsjizcuVKJSYmKj4+XpJ07Ngx/fbbb5Kkhx9+WCVLltSgQYNUo0YNtWrVSv7+/jp9+rSmTZumsLAwzZw5M0t7vXr1Uvfu3TV06FDFxcWpbt26CgoKUnBwsObMmWO5LcnaOAAAAABwhLBrV7RicZBWLA6Sh2cJZWSkyzCbHZ1WvjAMQ7/M+FThoVf1wtC35XqH72PnTp/QtH+9q4jr1wowQ/s5c+Ko3h/2rN4Z/7HqN2rm6HTgAC6Gvcd62UGtWrV08eLFHPedP39etWrV0tSpU/Xrr7/q/PnzSkhIkJ+fnx544AF98MEHat26dbbjEhISNHr0aC1YsEDR0dFq0KCBPvjgAw0YMCBXcTnp3L6lKvlw3yAAAAAA2NN97TpqxKgP5eVdMtu+9St/149f/kfp6WkOyMy+3D089OJr7+nBXo/bfCwTAjveqFGjNGXKlFwd65TFmcKK4gwAAAAA5I/a9Rro/UmfqVz5zAJEWmqKfvzqP9oQ/KeDM7O/no89qedeeSvLXKh3Q3HG8fJSnHG1cy4AAAAAANjd+dMnNOq153Xp/BmFh17V2BH/KJKFGUkK/mOBPnx/mOJibjg6FRQQp5xzBgAAAACAv4qKCNPYEf+Qq5ubEuPjHJ1Ovjp6aJ8+GP6cRk74RLXq3uvodJDPGDkDAAAAACg0kpMSi3xh5qaIsFCNGTFE2zeudnQqyGcUZwAAAAAAcFJpqan67MNRmvfjVzKbTI5OB/mE4gwAAAAAAE7u9/kzNeO/k8WaPkUTxRkAAAAAAAqBjauXKviPXx2dBvIBxRkAAAAAAAqJWdM/1bFD+xydBuyM4gwAAAAAAIWE2WzSfye/r8jw645OBXZEcQYAAAAAgEIkLuaGpk18V2mpKY5OBXZCcQYAAAAAgELm7Klj+v7zqUwQXERQnAEAAAAAoBDatGYZEwQXERRnAAAAAAAopJgguGigOAMAAAAAQCF1c4Lg6KgIR6eCPKA4AwAAAABAIRYXc0Pf/fdfTBBciFGcAQAAAACgkLt07jQTBBdi7o5OAAAAAAAA5N2mNct0T/2G6vn4U3Zr8+K509qw6k/t2bZRFSsHqF7DpqrfsKnqNWyqMuX87Hae4o7iDAAAAAAARcSsGf9VjTr11KjZfbluIyEuVlvXB2vD6qU6f/qEZXtEWKiO3jL5cKUqVVWvUTPVb9BE9Rs1U4069eTuTpkhN3jVAAAAAAAoIkwmk/499k01anafqteqqxq171GN2nUVUK2m3D08bnuc2WTSoX07tWHVUu3dsUkZ6el3PVdY6FWFhV7V1nUrJUmeJUronvqN9NhTz+u+tg/Y7TkVBxRnAAAAAAAoQpKTErVv5xbt27nFss3NzU0B1WtmFmxqZRZsqte+RxnpGdq0Zpk2rVmuG3lc8SktNVXHjxzQ8SMH1KXHI3ph6NsqWconr0+nWKA4AwAAAABAEWcymXT5wjldvnBO2wvgfBtXLdWR/bv18ptjFNi6fQGcsXBjtSYAAAAAAGB3URFhmjLqNX376YdKSkxwdDpOjeIMAAAAAADIN+tWLNE7Lw3Q4f27HJ2K06I4AwAAAAAA8lVk+HVNfm+Yvv/8IyUnJTo6HadDcQYAAAAAABSINcsW6Z2XByrk4F5Hp+JUKM4AAAAAAIACE3H9mv418hV9MXWsjh3aJ7PZ7OiUHI7VmgAAAAAAQIHbum6ltq5bqUpVqqpz9z7q/FAfVahUxdFpOQTFGQAAAAAA4DBhoVe14JdvtXD2d2oS2FpdejyiNvd3VQkvL0enVmAozgAAAAAAAIczDENHDuzWkQO75V2ylO7v+pC6PPSo6jVsIhcXF0enl68ozgAAAAAAAKeSnJSotcuXaO3yJapavZY6de+tjt16yb9iZUenli8ozgAAAAAAAKd19fIFBf30teb//I0aN2+lTt17q+0DXeVdspSjU7MbijMAAAAAAMDpGYahkIN7FHJwj3744iO1eeBBder2sJrd10aubm6OTi9PKM4AAAAAAIBCJS011bLaUzk/fz3QrZc6d++tGrXrOjq1XKE4AwAAAAAACq0b0ZFaunC2li6crYDqNRXYuoMCW3VQo+b3ydOzhKPTswrFGQAAAAAAUCRcu3xR1y5f1IrFQfIsUUKNm7VUYJvMYk2VajXueKzZbNb1q5d0/sxJnT9zUhfOnNT5syfl6uqmqtVrqmqN2pmP6rVUtUYtla9QyW6rSFGcAQAAAAAARU5aaqoO7NmuA3u2S5IqBVRTi9YdFNi6g+5t3Fzh16/+rwhz5qQunD2l1JTkHNuKvRGlY4f3Z9lWwstbVavXUkCNWqpWvVaecnUxDMPIUwuQJJlMJgUGBqp3795ydXV1dDpAvjGbzdq2bZvuv/9++jqKLPo5igv6OooD+jmKA/q5c7h48aLmzp2bq2MpzthJenq6PD09lZaWJg8PD0enA+Qb+jqKA/o5igv6OooD+jmKA/p54UdJDQAAAAAAwIEoztiJq6urxo8fzxAyFHn0dRQH9HMUF/R1FAf0cxQH9PPCj9uaAAAAAAAAHIiyGgAAAAAAgANRnAEAAAAAAHAgijN2kJCQoBEjRiggIEBeXl4KDAzU/PnzHZ0WYLFx40a5uLjk+Ni5c2eWWFv6s7WxfEZgb/Hx8Xr33Xf10EMPqUKFCnJxcdGECRNyjHV0n6b/Iy+s7eu2XOcl+jqcy/r16zVkyBA1aNBApUqVUtWqVfXYY49p37592WK5pqOwsrafcz0vxgzkWffu3Y2yZcsaM2bMMNavX2+8+OKLhiRj7ty5jk4NMAzDMDZs2GBIMqZMmWLs2LEjyyM+Pj5LrC392dpYPiOwt/PnzxtlypQxOnXqZOlP48ePzzHW0X2a/o+8sLav23KdNwz6OpzLE088YXTt2tX45ptvjI0bNxoLFy402rVrZ7i7uxvr1q3LEss1HYWVtf2c63nxRXEmj5YvX25IMubNm5dle/fu3Y2AgAAjIyPDQZkB/3PzIr9w4cI7xtnSn62N5TOC/GA2mw2z2WwYhmFERETc9guro/s0/R95ZW1ft/Y6bxj0dTifsLCwbNvi4+ONSpUqGd26dbNs45qOwszafs71vPjitqY8WrJkiXx8fNS/f/8s2wcPHqxr165p165dDsoMsJ0t/dnaWD4jyA83h/fejaP7NP0feWVtX7cFfR3OpmLFitm2+fj4qFGjRrp8+bJlG9d0FGbW9nNb0M+LFoozeRQSEqKGDRvK3d09y/ZmzZpZ9gPOYtiwYXJ3d5evr6969OihrVu3ZtlvS3+2NpbPCBzJ0X2a/o+CdrfrvERfR+EQGxur/fv3q3HjxpZtXNNR1OTUz2/iel78UJzJo6ioKPn5+WXbfnNbVFRUQacEZFOmTBm98cYb+vbbb7VhwwZ9/vnnunz5srp06aJVq1ZZ4mzpz9bG8hmBIzm6T9P/UVCsvc5L9HUUDsOGDVNiYqJGjx5t2cY1HUVNTv2c63nx5X73ENzNnYYb23soMpAbLVq0UIsWLSw/d+zYUX379lXTpk317rvvqkePHpZ9tvRna2P5jMCRHN2n6f8oCLZc5yX6Opzb2LFjNXfuXH355Zdq2bJlln1c01FU3K6fcz0vvhg5k0fly5fPsXoYHR0tSTlWHQFnULZsWfXp00eHDx9WcnKyJNv6s7WxfEbgSI7u0/R/OFJO13mJvg7nNnHiRE2ePFkffvihhg8fnmUf13QUFXfq5znhel48UJzJo6ZNm+r48ePKyMjIsv3IkSOSpCZNmjgiLcAqhmFI+l8F3Jb+bG0snxE4kqP7NP0fjvbX67xEX4fzmjhxoiZMmKAJEyZo1KhR2fZzTUdRcLd+fjtcz4sBB68WVeitWLHCkGTMnz8/y/aePXuy1BicWnR0tFG1alUjMDDQss2W/mxtLJ8R5Lc7LS/s6D5N/4c93amv5ySn67xh0NfhnP71r38ZkowxY8bcNoZrOgo7a/p5TrieFw/MOZNHvXr1Uvfu3TV06FDFxcWpbt26CgoKUnBwsObMmSM3NzdHpwho0KBBqlGjhlq1aiV/f3+dPn1a06ZNU1hYmGbOnGmJs6U/WxvLZwT5ZeXKlUpMTFR8fLwk6dixY/rtt98kSQ8//LBKlizp8D5N/4c9WNPXrb3OS/R1OJ9p06Zp3Lhx6tmzp3r37q2dO3dm2d+uXTtJjv89hX6OvLC2n3M9L8YcXR0qCuLj443XX3/dqFy5suHp6Wk0a9bMCAoKcnRagMVHH31kBAYGGmXKlDHc3NyMChUqGH379jV2796dLdaW/mxtLJ8R5IeaNWsaknJ8nD9/3hLn6D5N/0deWdPXbbnOGwZ9Hc6lc+fOt+3jf/26wjUdhZW1/ZzrefHlYhj/f/MaAAAAAAAAChwTAgMAAAAAADgQxRkAAAAAAAAHojgDAAAAAADgQBRnAAAAAAAAHIjiDAAAAAAAgANRnAEAAAAAAHAgijMAAAAAAAAORHEGAAAgl3755RdVqFBB8fHxdmvz2Wef1eOPP2639gAAgPNzMQzDcHQSAAAAhU1SUpLq16+vESNG6J133rFbu2fPnlWDBg20atUqPfjgg3ZrFwAAOC+KMwAAALkwffp0vfXWWwoNDVXZsmXt2vYjjzyi1NRUrV692q7tAgAA58RtTQAAoNhJSUlRixYtVLduXcXGxlq2X79+XZUrV1aXLl1kMpnu2Mb06dP1yCOPZCvMuLi4aPjw4fr555917733ytvbW61atdLOnTtlGIY+/vhj1a5dWz4+PnrwwQd15syZbG0/++yzWrt2rc6ePWuX5wsAAJwbxRkAAFDseHl5acGCBQoPD9eQIUMkSWazWU8//bQMw1BQUJDc3Nxue/yVK1d05MgRde3aNcf9y5Yt0w8//KCpU6cqKChI8fHx6t27t95++21t27ZNX331lb777jsdO3ZMf//73/XXgcxdunSRYRhasWKF/Z40AABwWu6OTgAAAMAR6tWrpx9++EFPPfWUPv/8c0VHR2vjxo0KDg5WlSpV7njs9u3bJUn33Xdfjvtv3pJU6v/au3tQavs4gOM/pRPdiySLVfISm8wGyqaIUiIZDJwMpOw2MtkYFKuQnUWdJAM2E2JRJyUv5e0e9Oi5e+635+6cLrfz+Wznqt/v/K/123V1ffkSEW9P03R2dsb29nYcHBxEUVFRRERcXV3F+Ph4HB8fR2Nj4/t8ZWVlVFVVxe7uboyNjeXidgGAD0ycAQAKVk9PT+zs7MTk5GQ8Pz/H9PR0tLW1/XLu8vIyIt4iyve0tra+h5mIiLq6uoiI6OjoeA8z/75+enr6TZz5Z/fFxcX/uyEA4K/ktSYAoKANDQ3F4+NjFBcXRzqd/q2Z+/v7iHh7Pep7ysvLv/mdSqV+ev3h4eE/O0pKSt7/BwD43MQZAKBg3d7eRn9/f9TU1ERpaWkMDw//1lxFRUVERGSz2bydLZvNvv8PAPC5iTMAQMEaGRmJs7OzWFtbi6Wlpdjc3Iz5+flfztXW1kZE5O1rSk9PT3F+fh719fV52Q8AfCziDABQkBYXF2NlZSUWFhaioaEhurq6YnR0NKampmJvb++nsy0tLVFaWhqZTCYvZzs8PIy7u7sffg0KAPhcxBkAoOAcHR1FOp2OgYGBGBwcfL8+OzsbTU1N0dvbG9fX1z+cT6VS0d3dHRsbG3k53/r6elRUVER7e3te9gMAH0vR6+vra9KHAAD42+zv70dzc3NkMploaWnJ2d7n5+eorq6Ovr6+mJmZydleAODjEmcAAP5Qb29v3N7extbWVs52Li8vx8TERJycnERZWVnO9gIAH5fXmgAA/tDc3Fw0NzfHzc1Nzna+vLzE6uqqMAMABcSTMwAAAAAJ8uQMAAAAQILEGQAAAIAEiTMAAAAACRJnAAAAABIkzgAAAAAkSJwBAAAASJA4AwAAAJAgcQYAAAAgQeIMAAAAQILEGQAAAIAEiTMAAAAACRJnAAAAABIkzgAAAAAk6CuMlWJqBykyUAAAAABJRU5ErkJggg==", | |
"text/plain": [ | |
"<Figure size 1050x450 with 1 Axes>" | |
] | |
}, | |
"metadata": {}, | |
"output_type": "display_data" | |
} | |
], | |
"source": [ | |
"h = zbed\n", | |
"\n", | |
"years = np.arange(0,2001,25)\n", | |
"\n", | |
"for i in range(1,len(years)):\n", | |
" dT = 5 * (np.tanh((years[i]-1200)/150)*.5+.5)\n", | |
" dTdz = 5/1000 #https://rmets.onlinelibrary.wiley.com/doi/full/10.1002/joc.3444\n", | |
" z_ela = 2000 + dT/dTdz\n", | |
" remaining_seconds = (years[i] - years[i-1]) * sec_per_year\n", | |
" Nsubsteps=0\n", | |
" while remaining_seconds>0:\n", | |
" #for plotting later.\n", | |
" thickness = h - zbed\n", | |
"\n", | |
" dhdx = np.diff(h)/dx\n", | |
" Hm = (thickness[:-1]+thickness[1:])/2\n", | |
" Dx = 2 * A * np.abs(dhdx)**(n-1) * (rhoi*g)**n * (Hm**(n+2)) / (n+2); #ralf greve for constant A\n", | |
" Qx = -Dx*dhdx\n", | |
"\n", | |
" dhdt = bdot(h,z_ela=z_ela)/sec_per_year - np.diff(Qx, append=np.inf, prepend=-np.inf)/dx\n", | |
"\n", | |
" maxDx = np.max(Dx)\n", | |
" if maxDx>0:\n", | |
" maxdt = 0.25 * dx**2 / np.max(Dx)\n", | |
" else:\n", | |
" maxdt = 10\n", | |
" dt = np.min((remaining_seconds,maxdt))\n", | |
"\n", | |
" h = h + dhdt*dt\n", | |
" \n", | |
" h = np.maximum(h,zbed)\n", | |
" \n", | |
" remaining_seconds=remaining_seconds-dt\n", | |
" Nsubsteps+=1\n", | |
" #remaining_seconds=remaining_seconds*0\n", | |
" \n", | |
" vol = np.sum(thickness)*dx*width/1e9\n", | |
" plt.figure(dpi=150)\n", | |
" plt.fill_between(x,zbed,np.min(zbed),fc='#554433')\n", | |
" plt.fill_between(x,zbed,h,fc='#99aaff')\n", | |
" plt.axhline(z_ela,color='k',lw=0.5,alpha=0.2)\n", | |
" plt.text(x[-1],z_ela,f'equilibrium line {z_ela:.0f}m',ha='right',fontsize='x-small')\n", | |
" plt.ylabel('z (m)')\n", | |
" plt.xlabel('x (m)')\n", | |
" plt.title(f'model-year: {years[i]} | $\\Delta T$:{dT:+.1f}C | vol: {vol:.1f}km3 | substeps:{Nsubsteps:.0f}')\n", | |
" plt.axis('tight')\n", | |
"\n", | |
" #old-school way to make an animation that is updating as it is being calculated\n", | |
" plt.show()\n", | |
" sleep(1/100)\n", | |
" display.clear_output(wait=True) \n", | |
" " | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [ | |
"\n", | |
"\n" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": null, | |
"metadata": {}, | |
"outputs": [], | |
"source": [] | |
} | |
], | |
"metadata": { | |
"kernelspec": { | |
"display_name": "base", | |
"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.11.9" | |
}, | |
"vscode": { | |
"interpreter": { | |
"hash": "a7ec67b8da850908bf022555edd6af3178010e99e99669fcaf3330f57f4b214a" | |
} | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 4 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment