streampower_diffusion_test1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "8ed4f412-f4f4-4da4-a84a-5b8d468e9c3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from landlab import HexModelGrid, RasterModelGrid\n",
    "from landlab.components import (\n",
    "    ChannelProfiler,\n",
    "    ChiFinder,\n",
    "    FlowAccumulator,\n",
    "    SteepnessFinder,\n",
    "    StreamPowerEroder,\n",
    "    LinearDiffuser\n",
    ")\n",
    "\n",
    "from landlab.plot import imshow_grid\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec4e7a38-1e02-42b1-8d4a-6c319a763cbf",
   "metadata": {},
   "source": [
    "# Simple example using Linear Diffuser & Stream Power components\n",
    "No uplift,\n",
    "not considering precipitation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "70df8dd7-5041-47fa-bf5a-7d5abd5744af",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGOCAYAAAC9oPjrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAABnwElEQVR4nO2dfXRV5Z3vvyfnNQnJgSSSwykBYxvrS/Cl4FCpLTi8WCqgpTPUYpXecu/Ci1JTQNRhWtMuTSqzisyFSi9eRxgYBte6LdbeaZXYViiLccQoLaDVWqME4RjBkDeS87rvH5ST/H472ZvDCckO+X5cey2fs9+e/exnb548v+/+/lyGYRgghBBCCHEQOYNdAUIIIYQQDQcohBBCCHEcHKAQQgghxHFwgEIIIYQQx8EBCiGEEEIcBwcohBBCCHEcHKAQQgghxHFwgEIIIYQQx+EZ7AoQQgghxExXVxdisVjWx/H5fAgEAv1Qo4GFAxRCCCHEYXR1daG8vByRSCTrY4VCITQ0NAy5QQoHKIQQQojDiMViiEQiaGw8gsLCwvM+TmtrK8rKxiEWi3GAQgghhJD+obCwMKsBylCGAxRCCCHEoRiGgWxy+g7lfMAcoBBCCCGOxfjrks3+QxN+ZkwIIYQQx8EZFEIIIcSxDN8ZFA5QCCGEEIcynDUoDPEQQgghxHFwBoUQQghxLAzxEEIIIcRxcIBCCCGEEIdBDQohhBBCiIPgDAohhBDiWBjiIYQQQojjGL4DFIZ4CCGEEOI4OINCCCGEOJThLJLlAIUQQghxLAzxEEIIIYQ4Bs6gEEIIIY5l+M6gcIBCCCGEOJThrEFhiIcQQgghjoMzKIQQQohjYYiHEEIIIY6DAxRCCCGEOAxqUAghhBBCHARnUAghhBDHwhAPIYQQQhyHASCV5f5DE4Z4CCGEEOI4OINCCCGEOJThLJLlAIUQQghxLMNXg8IQDyGEEEIcB2dQCCGEEMcyfGdQOEAhhBBCHMpw1qAwxEMIIYQQx8EZFEIIIcSxDN8QD2dQCCGEEMdi9MNy7uzZswdz585FOByGy+XCc889Z9rmrbfewrx58xAMBlFQUIDPf/7zOHLkSHp9NBrFsmXLUFJSgvz8fMybNw9Hjx7N9MI5QCGEEEKcylkNSjZLJnR0dODaa6/Fhg0bel3/l7/8BTfddBOuuOIKvPzyy/jDH/6A733vewgEAultqqqqsHPnTuzYsQN79+5Fe3s75syZg2QymVFdXMZQVtAQQgghFyGtra0IBoOIRN5CYWFBFsdpQyh0JVpaWlBYWJjRvi6XCzt37sTtt9+e/u2OO+6A1+vF1q1be92npaUFl1xyCbZu3Yqvf/3rAIBjx46hrKwMv/rVr3DLLbec8/k5g0IIIYQ4loEN8ViRSqXwH//xH7j88stxyy23YPTo0Zg8ebIIA9XX1yMej2PWrFnp38LhMCorK7Fv376MzscBCiGEEOJQ+ivE09raKpZoNJpxXZqamtDe3o4f/ehH+PKXv4xdu3bhq1/9KubPn4/du3cDACKRCHw+H0aNGiX2LS0tRSQSyeh8HKAQQgghFzllZWUIBoPppba2NuNjpFJnsirfdttt+O53v4vrrrsODz30EObMmYOf/vSnlvsahgGXy5XR+fiZMSGEEOJYUn9dstkfaGxsFBoUv9+f8ZFKSkrg8Xhw1VVXid+vvPJK7N27FwAQCoUQi8XQ3NwsZlGampowZcqUjM7HGRRCCCHEsfSPBqWwsFAs5zNA8fl8uOGGG/D222+L39955x2MHz8eADBx4kR4vV7U1dWl1x8/fhyHDh3KeIDCGRRCCCGEAADa29vx7rvvpssNDQ04cOAAioqKMG7cODzwwAP4+te/ji996Uu4+eab8cILL+CXv/wlXn75ZQBAMBjE4sWLsWLFChQXF6OoqAgrV67EhAkTMGPGjIzqws+MCSGEEIdx9jPjY8f+kPVnxuHwtef8mfHLL7+Mm2++2fT7okWLsHnzZgDAv/zLv6C2thZHjx7FZz/7WfzgBz/Abbfdlt62q6sLDzzwALZv347Ozk5Mnz4dTz75JMrKyjKqOwcohBBCiMPoHqAc6IcBynXn5YMy2FCDQgghhBDHQQ0KIYQQ4liGb7JADlAIIYQQh3I++XT0/kMVDlAIIYQQRzN0BxnZQA0KIYQQQhwHZ1AIIYQQx9I/TrJDEQ5QCCGEEIcynDUoDPEQQgghxHFwBoUQQghxLPzMmBBCCCGOY/gOUBjiIYQQQojj4AwKIYQQ4lCGs0iWAxRCCCHEsTDEQwghhBDiGDiDQgghhDgWGrURQgghxGFQg0IIIYQQB0INCiGEEEKIY+AMCoBUKoVjx46hoKAALpdrsKtDCCHEwRiGgba2NoTDYeTkXOi/84fvDAoHKACOHTuGsrKywa4GIYSQIURjYyPGjh17Qc9hGCkYxvkLXbPZd7DhAAVAQUEBACAMxrwIIYRYkwJwDN3/dpALAwcoQDqskwMOUAghhJwbAyMJYIiHEEIIIQ6DnxkTAMDRHv//rZtuEus8Sgjl9/lE2ed2i7LXI5vWq9ZrOuNxUX793XdF+XOf+Yzl8ZKp7jhjVB0rpTpojhr1f9LeLspjRo2y3F7zmqrrdZddJuuWTIqyrl8skZDrVVmjr92n2trv9VqW32psFOUrlf4orutjU1/dvqmUjPl6VH11e+q+ElD11ev19fako6tLlPe++aYo33TVVaKcHwiIsn6Z6X6p2yah7q2+Vt32mvciEVEuKymR54/FRFm3rUaf3637is1zqttWt+crf/qTKE+6/HJRLszNtTxeUtVf9503Vd+8LBTqc3vdL3VddVslbf6h0s+V7qe6bHevdVlf6+t/+Ysof+7TnxZlt3rn6rbT9fng449FefwllyAT3njvPVG+Xr3H1vzHf2R0PJI9HKAQQgghjoUhHkIIIYQ4DgPZ2dUP3QEKNaGEEEIIcRycQemDq9S37Tqemuf3i3Ku0qTo7XW8Vsfum1paRPn4J5+I8mdULLp4xAhR7hnf7VKxaX0urTz/6NQpUdaxWx0L7ohGRTnS3CzKnx0zRpR1nF/rCrp0WdVfo2PlAdX2+t6YdBXqfJ8Nh0VZt4+O9euybl8da9f11X3Bp/rKCJu+1fN4MXXuIyoOX1xYKMpji4tFeZy617puuq3srt1Ol6B1EbptPq36eaYaFDtNid9Gr6R1DR+q5/BdpZm5RLVv+ejRojwyP9+yvqfVs6R1Fldb+DPptjnV0SHKWitlJ5bU907fG922envdT7WWqvX0aVFuUu+dy9V7Q7ddwube56vnpkI/12p7u/bT79zBgiJZQgghhDgQalAIIYQQ4jiG7wCFGhRCCCGEAAD27NmDuXPnIhwOw+Vy4bnnnutz2yVLlsDlcmHdunXi92g0imXLlqGkpAT5+fmYN28ejh492vtBLOAMSh/o+GWuiqfm50u/A+9IGf9051n7P8Q+6RRlHcsvUbHt0mBQlMPjZaw7x98dD052SJ1AslP5iqhhaTAvT5Q/NVbqElweuUPrR62iPKaoSJS1l0XBCHn8WFTWT8eCdVmjvSy0BiU3V96Lzk4Z5z/RKuuvNTeBgDxePCbbT9cvrjUoKlaudRC6bNKYFEtvEk++XC/q8mGbKJ9ok+UC5csxQvmelBaNlHUJS21Tskteu+5LqZi8drdfXltOQN6reIu8F1qDoTUc2ndF/y2odQVaz+PLlWX3CKU58SptWLt1XxyldBHHlf6qQukoisPSU8hbKPtm9GOpy2hROo3xYdkenoLuvpBoUxqKU9LPSOuFtL5FY6ff0ZoUrefR/dQdkOtPfCj1PCHltxRW75GScbKsb34qav3cjR9XKsr6Xrd/LNur8cQJUdZ+UINHdrl4Mv0CqKOjA9deey3+23/7b/ja177W53bPPfcc/uu//gth9W8lAFRVVeGXv/wlduzYgeLiYqxYsQJz5sxBfX29yZvICg5QCCGEEMcysCGe2bNnY/bs2ZbbfPjhh7jvvvvw4osv4tZbbxXrWlpa8PTTT2Pr1q2YMWMGAGDbtm0oKyvDSy+9hFtuueWc68IQDyGEEELOiVQqhbvuugsPPPAArr76atP6+vp6xONxzJo1K/1bOBxGZWUl9u3bl9G5OINCCCGEOJb+mUFpVaFtv98Pv/o0+1x4/PHH4fF48J3vfKfX9ZFIBD6fD6NUiKy0tBQR9Zm+HRyg9EH4MzJ+6Q3KG+kbpTQoI2T81ZUjo+OxFpknI94qY+/aa6RI+Zzob/xzPyXXB0Z3l1NxFZtVOgEjJTu7+0+yG+RdKvUuHn1tHnltpY1y+6KxsmPq46WUriGu4v7Jdhk719/xu5WuQNfPo/Q/gYiMNYeV/8LI8bK+AaXD0LFurVPQOg2o9s1RugytT/IoXYQvKPtWjlf2jZ59J94s+5XWCWgNisnLokB5yIyT90prSlIJGc82VFnX1eWW5U51L8rapBYrWCF9WuySxWp9lNY9mHQR6npMHjmqfqM/kd4Y1156qSh/pPyLTN4hSnOSVya1ZZ58ee/HNsvrH/EZqcMIjO7WwOh+F4hI/ZF+jvR7AOodpdvOnav0RLof6+0Duq/I88VPyb5aUlAgyoV5sq/ml49Ux5dtpZ/DcapcUCHbLiegr0feq1BEnq+oRN6rwcIwstOgnN23THnqPPLII6iurs7oWPX19fjnf/5nvP766xlncjYMI+N9OEAhhBBCLnIaGxtR2OPji/OZPfn973+PpqYmjBs3Lv1bMpnEihUrsG7dOrz//vsIhUKIxWJobm4WsyhNTU2YMmVKRuejBoUQQghxLEY/LEBhYaFYzmeActddd+GPf/wjDhw4kF7C4TAeeOABvPjiiwCAiRMnwuv1oq6uLr3f8ePHcejQoYwHKJxB6YPgVfLTU2+BCvEE5eea+hO2eJuczuz8SE4d60/u2lWqdJ22XdtG66nr3FB3WMLtU1OtOuSjpuV1OMpfJD8L1lO3eqp2tPoEOhCWU7cFl6rP9VQ4Nd4hp2YTp62t7vU0vZ4m1/cix6emct8fKco6rDFivFxvJGV76folVQjIHOJRVvc6RKXLKgSUjOpPfXucX82Y6lBhrnoJ6U9H9b3V/TwvJO+l7is6XKjXm0JCavtRY0eKsm57XT+X26XK1iEe3Vf0DHO0WYaYutRzqtGf6urQq/6M2ztSlk3Plqpf0ccqPKqepbxwj7CDCk/5Vdg5cVqlCUhYhwl0CMWt+q3pE3L1XOWoexE9Jdu2s1R+ol1yXIZQ/KNl2+i28hfJ64u3Wb83RFvBHA7UIaLRI0eKsn6PDR4D+xVPe3s73n333XS5oaEBBw4cQFFREcaNG4dilS7D6/UiFArhs5/9LAAgGAxi8eLFWLFiBYqLi1FUVISVK1diwoQJ6a96zhUOUAghhBACAHjttddw8803p8vLly8HACxatAibN28+p2M88cQT8Hg8WLBgATo7OzF9+nRs3rw5Iw8UgAMUQgghxLH0l0j2XJk2bVpGCQbff/9902+BQADr16/H+vXrMzq3hgMUQgghxLEM31w8HKD0QUG5/ETNqzQhHo+MbyaTyrq+Q35GrIPfWmOi7eb/ourjU7FsT56yR+9R9udL/YwmmZSfTupYtneEEk+puvsvUbFk9bmgjhUHSuT2Hq+81mRCtV2ntQZFx77dXllfw5CajYTSuATDUnOiY/e5JfIz45wcuT6ZlHbkKWWFb/qMW3+e6VGfEefI9tCf4nWmzt07IEftq2307T7z0+t1P/P5pPW6y6V19vLaY7GPRdk7QrZd/qUjRTlvjOxL+hNyt1drSpQuIkc+J5rYaWlNr3ULMaWvylGaHp1W4VJlzT/yUyNFWT8LPvXZsdbM6M+QtZ2Bv7C7fdxu2U+9I06KcqJLfa5vp0ExaUxUv3XLurhcUrNiGOqz5qS0jtdWDfqdpz+/11oy03tJkT9WtZ3WCar3RqxYnr9kpNw/oDQzg4VhGBnNaPS2/1CFX/EQQgghxHFwBoUQQghxLAYyTfhn3n9owgEKIYQQ4lioQSGK/OB4UdbxVx0dSyaPyrKylU4qXYX/Ems/hFEdUieibad1XLFn0e1Wx1ax6lRKxtkhw+im7Q1Dajj8o2TdCkukbkDrZXx++d28zyfLhqHaKiB1ChqtM9A6hGRSelnERkmNi/Y30PX1epX+yFui6qti7QHZPvqvHV3fnBwZS9c6jkRC1l/7N/T0/nApa3ltZa/9c7RGRXu8mLwyVD/zeGRc3u2WbanbRmuzAiWy72kNRq7JY0HZlZvaTn+2KOsfi0ldRrxdasNOH5P28Mf+/JFcH5Xb5/ltdBQqbYGdL4vWWeQq3xmzj0v39j6fbBuPR+6bCsi21tosjdZa6bbOyZH6G7PeSLZ1jkd6Pem20G2n0e9MrT9yj1J+SMrHxWs6vuz7gWJ5bwNjVfqQYuv6kQsPByiEEEKIQxnOIlkOUAghhBCnYhimmcyM9x+i8CseQgghhDgOzqD0gd9fKspmnYPyElHxWJ0K3ae9QUqVLkTlWxl1tFWuj8nYeleTPH9+j7wTRqE8ltsdUGUda5axW+1vkEjIOL0vKGPNgTFSl+DVKe7d1h4yWkeQStn4oOTobqs1HLL+vqDUFQSUj4vOZ6OvX+sudF/Qugc9parrq/dPpWQsXGtQTBq3HmWTb4n27dA5nGysprWHiy5r3G6tS5C6Bb9f+aYEVe6ffNnP/f4xcr1XetbYtb1uO623Sqjn0khI/ZNur1CZ1B/lBGT7aU+gXPVc539K5akaofVM8vymvExe+ax4vSPT/+92K38hj6ybdhDVZd137N4D+kHR70CTlkzp8HTunDyVa6fws7JttAZEa8F0fb1e/d6T7x2twYkHZf31vdP+SGTg4QCFEEIIcSjDWYMyqCGe6upquFwusYRCofR6wzBQXV2NcDiM3NxcTJs2DYcPHxbHiEajWLZsGUpKSpCfn4958+bh6NGj+lSEEEIIGUIMugbl6quvxvHjx9PLwYMH0+vWrFmDtWvXYsOGDdi/fz9CoRBmzpyJtrbukENVVRV27tyJHTt2YO/evWhvb8ecOXOQTCZ7Ox0hhBAydDD6YRmiDHqIx+PxiFmTsxiGgXXr1mH16tWYP38+AGDLli0oLS3F9u3bsWTJErS0tODpp5/G1q1bMWPGDADAtm3bUFZWhpdeegm33HJLFvUqUL/Iu5xKyXiq9hDQ3hq5o2V8U+seYi0yVj7i06Pk8VQ81JPft9+C9ivQsVqzL4ddvhZ57b4Cub+Ow3sLVL4R5ati1sDI89tl39S+IeZYu9Sw+HwyN1FeqfTm8I3UGhmdi0dreKxzgpjRfwco7xHVl3Qs36QL6Tllq26V1pj4bTQoJo8VdTz76WG5v/bO6KmZOHN8lUNKaUy8SnOhdRYanRdJ+64kE3K9EVd9S+cuUu2VN17WLzekdArqOfePkrqJQO6nRNnjkcfTfdfrlT4wuj09nu7ja78jsz5HY30vzZ4y+jnTHjfWeZ1yVN/S78TAaHlve+roAMCXr/IyefR7xNrnRL/3kkn5jvXly/P5i5XWrkC9swcJhngGkT//+c8Ih8MoLy/HHXfcgffeew8A0NDQgEgkglmzZqW39fv9mDp1Kvbt2wcAqK+vRzweF9uEw2FUVlamt+mNaDSK1tZWsRBCCCGOYxjPoAzqAGXy5Mn413/9V7z44ot46qmnEIlEMGXKFJw8eRKRyJkMrqWl8mua0tLS9LpIJAKfz4dRo0b1uU1v1NbWIhgMppeysrJ+vjJCCCGEZMOghnhmz56d/v8JEybgxhtvxKc//Wls2bIFn//85wGYww2GYdimjLfb5uGHH8by5cvT5dbWVg5SCCGEOI9hbNQ26BqUnuTn52PChAn485//jNtvvx3AmVmSMWO6vRGamprSsyqhUAixWAzNzc1iFqWpqQlTpkzp8zx+vx9+v7WOQMd+9Tf0Or6pNSv+oNQ1+ANSZ6Pjgm6fnPHJHydj1b6gjN9qPLnd5zPnK9EaFOuy9iHRsWmtC/ArPwOPV8fGtc5B6SBc2n8BlpgHrVKzoWPTOp9LoOiUrK9HzsBpnYS5vtrXJLOJSK05sdcGqHJO9w9aQ+K10aDo9aa6qVw8RsJaL6P1NHZ9xfxcyedGt71d39TPpdagpOLWYvkclcvIozx8tNYrV+mtAiNl3/L5Rquy7Ftaq6bvvb4e3dd6tof5ubG+t3b09segLOu219ovWfeU8kHxKy8oXc4tlO9Ij2ekKGu9kl1OKx0g0O3j9cp74y9qVuvlvR0shvH4ZPA1KD2JRqN46623MGbMGJSXlyMUCqGuri69PhaLYffu3enBx8SJE+H1esU2x48fx6FDhywHKIQQQghxNoM6g7Jy5UrMnTsX48aNQ1NTEx599FG0trZi0aJFcLlcqKqqQk1NDSoqKlBRUYGamhrk5eVh4cKFAIBgMIjFixdjxYoVKC4uRlFREVauXIkJEyakv+ohhBBChizDeAplUAcoR48exTe+8Q2cOHECl1xyCT7/+c/jlVdewfjx4wEAq1atQmdnJ5YuXYrm5mZMnjwZu3btQkFB97TwE088AY/HgwULFqCzsxPTp0/H5s2b4baZyiaEEEKIcxnUAcqOHTss17tcLlRXV6O6urrPbQKBANavX4/169f3a93MXhtaGCHLenuvT8Yvtb+B1k3E45+Icp7yBPCOkPFWHVvvmd/GLjZr53eQo/K56FisrnugRNY9N7dc7T9SnV/HznWk0UaEYkL+haDj9tp3xe8Pi7L2pvB4tE+LnebELlKq/4Kxvh8mHxuv0nX08Lzx5Km8Qx5ZV78uK01KskN5W6hcNUmVI0rntkmldP4TrUnR99r6lWMWt+vnTB5ft5XWFRj5qv4lsu/G2mQeJM9HOseWRHvS6PObPX5y1XprLZmZvt8z5uc40+fGDpW3yORPpO+F0jsVyvdQXki+0/R7Sut3zL4nsi3Nfcnu4wmtVVOeNkHthzTS8njkwuMokSwhhBBCuhnORm0coBBCCCFOJVuztaE7PnHWVzyEEEIIIQBnUCyw1pjo+KnbXWC5Xn/Drz0D/H7pAeC+RHp56OMlkzJvRM/Yttm/QMeSdRxdayL0tcuyyfNlxEjL9R6P9I7I1kfEjK6/PL45t4722vCrsp1uIjPNiX1uIa1jkPX1+qR+yT+yWwfiK5Zx87w8le/EJ499QqV1aOuUviF+lR8lMUbeS5PPSCqq1kuNitaqm++99pixblutBzIMfe+0/krduxyVr0X5yLhNOa5U2a/rr/uSnYeOnR5M07cmx86wMlstF6C1UVozon1C5LUlR0g9TyAw1nJ7rVWzfy6ttWv6udP7m/2SSkRZv7cGDX7FQwghhBCnMYzHJxygEEIIIY5lGI9QqEEhhBBCiOPgDEofmOO7WoOic/HoHBvaD0Hn/pGxap/SGWgdh463plLS76GnR4COrZrj4na+I1qDItdqPwLtX6C9Huw1J9n6nljfGzv9kFkXkKnPic5JIv0W7PLH6PUak+ZnVPf9zR+rtE0J2TZXhaWeJd4qNSPNx0/Jk6m/tnQ+lXin1LDk5MgcUtonRWuzdFzfrDOw66t2yPqa8sckVd9xq1w8ecpXxU6T4tZ5p6xfqWY9kl3fs9aDWW2bOXa6O3nt2htKP/da06E9asx+Rfq9Za39sssdZOepozHn+snUs+bCYCDLCZR+q8nAwxkUQgghxKmcDfFks2TAnj17MHfuXITDYbhcLjz33HPpdfF4HA8++CAmTJiA/Px8hMNh3H333Th27Jg4RjQaxbJly1BSUoL8/HzMmzcPR48ezfjSOUAhhBBCCACgo6MD1157LTZs2GBad/r0abz++uv43ve+h9dffx0///nP8c4772DevHliu6qqKuzcuRM7duzA3r170d7ejjlz5iCZtM4urmGIhxBCCHEqA2zUNnv2bMyePbvXdcFgEHV1deK39evX42/+5m9w5MgRjBs3Di0tLXj66aexdevWdNLebdu2oaysDC+99BJuueWWc64LByjnjF18VueJyMzvQOsM7LwzdKy9Z84TOx8PszeDXQ4LfS0q54aK3WaaIyNzrI9nH8fXT6x1nN9sFS3vjdaQaG+QREL6QZi9RLrU+nZVlvsnurrP51K+Hr6RSv/jdlmWvR/Le9WlctHo7ePt8tpyS2Xd/UUnZTlwiSinUlJr5fXK/Cwuly7Lvma+99pzRj4Xui2jzartYyqXUJ5sjxzle2LWYcREOZGQ59N9w+3W/kmZ6cX0syi3Nf3S57bnhp2GQ3nKqLr3zA8G9JZLR7alvc+Jdf309evHVh9P10ffG309g4Xx1/+y2f9C0tLSApfLhZEjRwIA6uvrEY/HMWvWrPQ24XAYlZWV2LdvHwcohBBCCOmmVZk0+v1++P36443M6OrqwkMPPYSFCxeisPDMgDQSicDn82HUKCmKLi0tRSQS6e0wfUINCiGEEOJUjH5YAJSVlSEYDKaX2trarKoVj8dxxx13IJVK4cknn7S/DMM4B/djCWdQCCGEEKfST0ZtjY2N6VkOAFnNnsTjcSxYsAANDQ347W9/K44bCoUQi8XQ3NwsZlGampowZcqUjM7DAUqf2HmD6DwV1t/omyerdD6cTDuL1EH0rI9Z72JXNztNhy7L47vddrl2+luDYod1bDpb3xWzzkHqMOLxU6Ici30syomEXB9X+XASHVLnEVPeJT29TGKn5L6xT5SeJSo1FsnTshzwKj1RocohpTQtBZfKadvckaVye+V1ofuGjvtrrwlz37W7V9aeNfr42tfEVyCvt+PDFlFOdEiNyemPT4lyPP+0KHty5fE8nqAqS12Gvf+R1ot5evy/nZ+QtdYqUzL3htKeN1pf09/+SBrrHF1mPyTd1nZ5koYWhYWFYiBxvpwdnPz5z3/G7373OxQXyz48ceJEeL1e1NXVYcGCBQCA48eP49ChQ1izZk1G5+IAhRBCCHEoA+10397ejnfffTddbmhowIEDB1BUVIRwOIy/+7u/w+uvv47/9//+H5LJZFpXUlRUBJ/Ph2AwiMWLF2PFihUoLi5GUVERVq5ciQkTJqS/6jlXOEAhhBBCnMoAj1Bee+013Hzzzeny8uXLAQCLFi1CdXU1nn/+eQDAddddJ/b73e9+h2nTpgEAnnjiCXg8HixYsACdnZ2YPn06Nm/eDLdOb24DByiEEEKIQxnoGZRp06b1Yq3Q83j2BwwEAli/fj3Wr1+f2ckVHKCcI+b4qx4J2uWBMB1RlHT81oxdp3D18f+Za07sMPsJaP2M1mxovYy1HqeXM5575fpl/8zqo3OG6JwkWoehfWO8XuV7ki+9SHJLpK4k0dWtO+lskvue6pJ6lxGfKZI1t2kKd67KNeOT/dxfKOvu948RZa2hyDTOn3meKOvnSPuMePJlfXx5SvPhUboF5QMTKNDtqb0/tO5C5+jSHjvaU0fea338nJxz78uZ5zHKFGuNh9YbmfVC/ftesteead2g9l3JVP9ELjQcoBBCCCFOZaCnUBwEByiEEEKIUxlgq3snQaM2QgghhDgOzqCcM3ZeIf2r68iOC5v7xpyHSGouTHtnqMcZemTmB+HxjFBlqYNIpaT3hvZdiXm6890kTh+Rx1I+H/5RUgfgHaE0E/lK4+CW15JKSI2E11siylrjYc4tY51Lx94vyA47fyJ1/bmyfj7faLn9KGs/Ir8/pH6x1mPZeX2YNSZ2ebN6lu38jfobO42Hfg6s8woNfn0H26/p3DAM45yEqVb7D1U4QCGEEEKcCkM8hBBCCCHOgTMohBBCiENhiIc4AGfGP3vHLrZrvf1Qw86bQ3tT6Fi29rrQugit4TF7Y8hcPMlktw+KO6C8J/Kk5iPWIn01UjGpZwmMsvY16XkuoDefEzuNhXWuHI257axfrnZtpY/n8YwUZbtcQfr8OteQnSeOWSeSWftY9T3nabusNSb9n2snU/o7R9cAwRAPIYQQQohzcMwApba2Fi6XC1VVVenfDMNAdXU1wuEwcnNzMW3aNBw+fFjsF41GsWzZMpSUlCA/Px/z5s3D0aNHB7j2hBBCyAXgrFFbNssQxREDlP3792PTpk245pprxO9r1qzB2rVrsWHDBuzfvx+hUAgzZ85EW1tbepuqqirs3LkTO3bswN69e9He3o45c+YgmUzq0xBCCCFDimE8Phl8DUp7ezvuvPNOPPXUU3j00UfTvxuGgXXr1mH16tWYP38+AGDLli0oLS3F9u3bsWTJErS0tODpp5/G1q1b02mct23bhrKyMrz00ku45ZZbBuWahh8Ojd2eN3aeN/oX67xMWmNi5w1izusk/47weLp1I/4RxWJdwXjpoaJz68ROSU2KzpWjNSY+n/Y90foZ7eOhr83O+yLTtrb2PdH1074sWkOirx9IwQr73ELWOgv767+webT6lyGq6Ujj9Pr9FQNZWt33W00GnEGfQbn33ntx6623pgcYZ2loaEAkEsGsWbPSv/n9fkydOhX79u0DANTX1yMej4ttwuEwKisr09v0RjQaRWtrq1gIIYQQ4hwGdQZlx44deP3117F//37TukgkAgAoLS0Vv5eWluKDDz5Ib+Pz+TBq1CjTNmf3743a2lr84Ac/yLb6hBBCyIWFX/EMPI2Njbj//vuxbds2BAJ9W6X39tmh3ed1dts8/PDDaGlpSS+NjY2ZVZ4QQggZAAwYaS+U81qG8Ahl0GZQ6uvr0dTUhIkTJ6Z/SyaT2LNnDzZs2IC3334bwJlZkjFjur0Zmpqa0rMqoVAIsVgMzc3NYhalqakJU6ZM6fPcfr8ffr91zg1C+ibTnCRSl2HvDaI1KPIx7ZnbR+eGMS6Jq32lRiIZbVLrdd36PhfQW34Vax+UC58XSufe8auyrK8Z67Y3nT3L63Wed0l/MpTrTpzIoM2gTJ8+HQcPHsSBAwfSy6RJk3DnnXfiwIEDuOyyyxAKhVBXV5feJxaLYffu3enBx8SJE+H1esU2x48fx6FDhywHKIQQQsiQwOiHZYgyaDMoBQUFqKysFL/l5+ejuLg4/XtVVRVqampQUVGBiooK1NTUIC8vDwsXLgQABINBLF68GCtWrEBxcTGKioqwcuVKTJgwwSS6JYQQQoYcWVrdD+XvjAf9M2MrVq1ahc7OTixduhTNzc2YPHkydu3ahYKC7k8hn3jiCXg8HixYsACdnZ2YPn06Nm/eDLdbf/pJCCGEkKGCyxjKmYT6idbWVgSDMicJm4WcO7Kv6PwwhiFNA806Buv8M3r/VKrb6ySZbBfrYrFmUU4k2kT5dIt0Wc4tlBqWQEDm4vF6R4qy9hm58BoU3RbWbWMYWoOjNTZ2viL6uc/Up0RDXcbFQm/3uqWlBYWFhRfkfGf/XXp9108xIj/Xfoc+aO/oxOdm3XNB63qhcPQMCiGEEDKsydYOdgj/sc0BCiGEEOJQhvH4ZPCdZAkhhBBCNJxBIaSfMesw9PpMvTJ0PpucHv+vc8HIfRMJqUlx++Qjb+drMvh/w9jlKdLrrdvDPjdQpvUh5AIzjKdQOEAhhBBCHMowHp8M+p9HhBBCCCEmOINCCCGEOJVhPIXCAQrpB7J9AIZ6XD9TL41M91dbCw2K9iGRvh85OTIRpzswQpbdeX0e+0zZzjdkYO+dneYkc9+Sod73yEUPsxkTQgghhDgHzqAQQgghDsXIMhfPUHZF5wwKIYQQ4lTOalCyWTJgz549mDt3LsLhMFwuF5577jlVHQPV1dUIh8PIzc3FtGnTcPjwYbFNNBrFsmXLUFJSgvz8fMybNw9Hj8o0G+cCByikH3DZLMMNef0ul1wybx+rts0Ri8vlFYvbnSeWQKBMLH5/qVjc7nyxuFwesQz+vbXuay5XjlqybfuBxriAi925CAE6Ojpw7bXXYsOGDb2uX7NmDdauXYsNGzZg//79CIVCmDlzJtrauvN+VVVVYefOndixYwf27t2L9vZ2zJkzB8lkstdj9gVDPIQQQohDGeiPeGbPno3Zs2f3cSwD69atw+rVqzF//nwAwJYtW1BaWort27djyZIlaGlpwdNPP42tW7dixowZAIBt27ahrKwML730Em655ZZzrgtnUAghhBCn0k8hntbWVrFEo9GMq9LQ0IBIJIJZs2alf/P7/Zg6dSr27dsHAKivr0c8HhfbhMNhVFZWprc5VzhAIYQQQpxKP0X3ysrKEAwG00ttbW3GVYlEIgCA0tJS8XtpaWl6XSQSgc/nw6hRo/rc5lxhiIc4gEx9Q4Ya2eZ76bt9tG9JTo72QckVZbdb+aK45Xpzbh/9N0ymPikDzWCf3w7r+Xb9xUXmuYLO/9zk4qaxsRGFhYXpst/vP+9j6X5pGIZtXz2XbTScQSGEEEIcytnPjLNZAKCwsFAs5zNACYVCAGCaCWlqakrPqoRCIcRiMTQ3N/e5zbnCAQohhBDiVC7EB1znSXl5OUKhEOrq6tK/xWIx7N69G1OmTAEATJw4EV6vV2xz/PhxHDp0KL3NucIQDyGEEEIAAO3t7Xj33XfT5YaGBhw4cABFRUUYN24cqqqqUFNTg4qKClRUVKCmpgZ5eXlYuHAhACAYDGLx4sVYsWIFiouLUVRUhJUrV2LChAnpr3rOFQ5QLkqGe26ci5uecVytGcnJkY+01pxojYpZc2KnQTHVxmY9scLs8plS689d42Mf3+e9GooMtJPsa6+9hptvvjldXr58OQBg0aJF2Lx5M1atWoXOzk4sXboUzc3NmDx5Mnbt2oWCgoL0Pk888QQ8Hg8WLFiAzs5OTJ8+HZs3b4bbrXNnWeMyhrIPbj/R2tqKYDAofhvazeK0AUqm9eGLtCeGkVK/dJdTqbhYk0x2iHI83iLKPl+xKJuTBXKAcmGRz4LdACWT5IwcoFxYemvflpYWITztT87+u/TK/12HEXm59jv0QfvpTnz+76ouaF0vFNSgEEIIIcRxMMRDCCGEOJWBtpJ1EBygXDR0d8JMw1Pmqcuh26EvRsyeA2Kt2tajyj51tGxDNpmGCS52jxs77EI61ut1iC0zH4nh1tYXJ8xmTAghhBDiIDiDQgghhDiVbL1Mhu4ECgcohBBCiFMZziEeDlCGLDp2nezx//pTRY3dp4vZ6RDOJSeD9fbDXbeQCbofyHvvdksNiv6MuP81J8Qafb8San027c17dVGSgvnr80z3H6IMqgZl48aNuOaaa9K5AW688Ub8+te/Tq83DAPV1dUIh8PIzc3FtGnTcPjwYXGMaDSKZcuWoaSkBPn5+Zg3bx6OHj060JdCCCGEkH5kUAcoY8eOxY9+9CO89tpreO211/C3f/u3uO2229KDkDVr1mDt2rXYsGED9u/fj1AohJkzZ6KtrS19jKqqKuzcuRM7duzA3r170d7ejjlz5iCZTPZ1WkIIIWRocPYz42yWIco5D1AuxKzE3Llz8ZWvfAWXX345Lr/8cjz22GMYMWIEXnnlFRiGgXXr1mH16tWYP38+KisrsWXLFpw+fRrbt28HcMbF7+mnn8aPf/xjzJgxA9dffz22bduGgwcP4qWXXur3+hJCCCEDyTAen5z7AKWyshJbt269YBVJJpPYsWMHOjo6cOONN6KhoQGRSASzZs1Kb+P3+zF16lTs27cPAFBfX494PC62CYfDqKysTG9z4eiPFJPnvxhGSiypVCK9JJOdYkmlomqJqSVuuRhG0nLpDpKeWczbyLqaWrKX1ODECld6cbly1OIRi9udL5acHK9YXC63WlxiyRy7NKoD/aw4G/Oz4VILIcOXcx6g1NTU4N5778XXvvY1nDx5st8qcPDgQYwYMQJ+vx/33HMPdu7ciauuugqRSAQAUFpaKrYvLS1Nr4tEIvD5fBg1alSf2/RGNBpFa2urWAghhBDHMYynUM55gLJ06VL84Q9/QHNzM66++mo8//zz/VKBz372szhw4ABeeeUV/M//+T+xaNEivPnmm+n1ZhdN45y+ErHapra2FsFgML2UlZVldxGEEELIBWAYj08y+8y4vLwcv/3tb7FhwwZ87Wtfw5VXXgmPRx7i9ddfz6gCPp8Pn/nMZwAAkyZNwv79+/HP//zPePDBBwGcmSUZM2ZMevumpqb0rEooFEIsFkNzc7OYRWlqasKUKVP6POfDDz+cTiENnMkayUEKIYQQ4hwy9kH54IMP8LOf/QxFRUW47bbbTAOUbDEMA9FoFOXl5QiFQqirq8P1118PAIjFYti9ezcef/xxAMDEiRPh9XpRV1eHBQsWAACOHz+OQ4cOYc2aNX2ew+/3w+/3Z1ozUz0l2X5s3n9+FFrnkZPjV2W7e9a/3hjnp2Xoid2fAMMtVt/zenWuFulzkpPjtVyf/b2xxk5T1NsMqc0Rbcq6PWwOZ3e2LPNa2efekc+qy2V3fcOtrxMmCzxHnnrqKaxYsQIzZszAoUOHcMkll2R18n/4h3/A7NmzUVZWhra2NuzYsQMvv/wyXnjhBbhcLlRVVaGmpgYVFRWoqKhATU0N8vLysHDhQgBAMBjE4sWLsWLFChQXF6OoqAgrV67EhAkTMGPGjKzqRgghhAw6HKDY8+UvfxmvvvoqNmzYgLvvvrtfTv7RRx/hrrvuwvHjxxEMBnHNNdfghRdewMyZMwEAq1atQmdnJ5YuXYrm5mZMnjwZu3btQkFBQfoYTzzxBDweDxYsWIDOzk5Mnz4dmzdvhtutHTMJIYQQMlRwGec4hzlz5kw888wzGDt27IWu04DT2tqKYDAofss0LfrAh3jk+VOpePr/k8lOsc7tzhXloRfisT3DBT7+YNN33+uZ4gAAUqmYKNuHeOx08pm2rd1zoo5+wUM82fWN/g/xyPdEMhkV5ezu18X+HAwuvfWllpYWFBYWXpDznf136ff/8jhG5OXa79AH7ac78cVvP3hB63qhOOcZlLq6ugtZDwdiHSs2/8MQhxXmzm1dNm+vX0y6fole/7+3bTM/92DjtPpcaDL5R1pua59rp7/bMrMBhdkHx65fW/8Db9Zw6OvXx7f7AlDXL1NNSGbnMw9Y5B8XZv1Y9wDGfnBJLgoMZBni6beaDDhMFkgIIYQ4lGEsQRncXDyEEEIIIb3BGRRCCCHEqQzjKRQOUPpAx761rkOL25LJDlE2x4et/Soy3V4HFuPxlnPe1qxbyHQibaB1DBcbmV2vtVBT3gstssxcX9TfolL9HOks43Zar4RlWWu/3O6AKLtcWnSqRax257Ouv352zD40Pkj0/ZD1S6VOq+3tNDDkomcYD1AY4iGEEEKI4+AMCiGEEOJQDAPoJQl8RvsPVThAIYQQQpzKMA7xcIDSJ9ax8VRKalDi8VOW25s1J7LptXmaeb21gVPP2LnbnQdJpt4MA+07Ym1Cd/Fh175aI3TuXib9b7yWHXYaDjtNiWFI47lkskuVpWbD6x0pym53vk399PGjar2un9zerDmRviUezwi1Xj7Hdu1h1shYeeAMN78gcrHDAQohhBDiUIbxBAoHKIQQQohjGcYjFH7FQwghhBDHwRmUPjBrTmRsXPuexGIfW64359SQZe2HoP0TzOutNStDm+EeSz93jcrAa04yTeZn7TOitVw6F43WmCQS7Wq9LOfkBFRZPmf6+KmU1LTEYp9AIq8nkWgRZfNzLc/v95eKsk7caU4WqJ/jofvXb+Z1H+7Pfe8YhpFx0kq9fyYkEglUV1fj3/7t3xCJRDBmzBh861vfwj/+4z8iJycnfcwf/OAH2LRpE5qbmzF58mT85Cc/wdVXX33e9ewNzqAQQgghTiXVD0sGPP744/jpT3+KDRs24K233sKaNWvwT//0T1i/fn16mzVr1mDt2rXYsGED9u/fj1AohJkzZ6KtrS3Li5VcTH92E0IIIRcXA6xB+c///E/cdtttuPXWWwEAl156Kf793/8dr7322l8PZ2DdunVYvXo15s+fDwDYsmULSktLsX37dixZsuT866rgDAohhBBykdPa2iqWaDTa63Y33XQTfvOb3+Cdd94BAPzhD3/A3r178ZWvfAUA0NDQgEgkglmzZqX38fv9mDp1Kvbt29evdeYMSh+Y/Qmk/4GOhcc6T4lyjkeO/XwBGYs2+yNozYldbh4Zr02lumP7Zn0Kx6EXFz3v/eBqFOxy72h0rhztC2KnIdHPZbQrIsopv9SUJBKy78fjzaKstWKdrU3y/B7lc+KVZb8/rOpn/Z7Qmheztixfre/7WR5qvid2Wgjz5Qyt67tQ9NcESllZmfj9kUceQXV1tWn7Bx98EC0tLbjiiivgdruRTCbx2GOP4Rvf+AYAIBI588yVlsp/00pLS/HBBx+cf0V7gQMUQgghxKmkjDNLNvsDaGxsRGFhYfpnv9/f6+bPPvsstm3bhu3bt+Pqq6/GgQMHUFVVhXA4jEWLFqW3MyfeNPp90MwBCiGEEHKRU1hYKAYoffHAAw/goYcewh133AEAmDBhAj744APU1tZi0aJFCIVCAJD+wucsTU1NplmVbOHcPyGEEOJQzn5mnM2SCadPn05/TnwWt9uNVOrM50Dl5eUIhUKoq6tLr4/FYti9ezemTJmS/QX3gDMofaL9G7QmRfqipBJyvTdQIMp2mhKzn4LWpMjYt8blivf4f31bzz2XCxlqXOh7qXPp2OXW0X5BUojndutpZVl/8xSxfA61j0miQ2o+4v5Tcr3yLYl2Sb+iHLd8VlIxeT25hSFR9vlKRNnrDcIa6+vT7Wn2R9LPvZOe3UzDDtZ6JcOw1tkNWwxkJzXLcN+5c+fisccew7hx43D11VfjjTfewNq1a/Htb38bwJn7UlVVhZqaGlRUVKCiogI1NTXIy8vDwoULs6ioGQ5QCCGEEAIAWL9+Pb73ve9h6dKlaGpqQjgcxpIlS/D9738/vc2qVavQ2dmJpUuXpo3adu3ahYKCAosjZw4HKIQQQohTGWAflIKCAqxbtw7r1q3rcxuXy4Xq6upevwLqTzhAOU9MadATsnz6xAlRThXJzws9Hjk17PePEWWfb5Qomz83lFO/bnfPsp5W5tTpxYvdyyfTe20d0jHbUlqnhNAhFmCEKsu+abbCj6myDPHE2+X6To/87NhIyvpHT0rr/Byf/oxYh2Ll9enQq/mzaPnZsB36PaL3t352s32Os/tEPVNtg+mdqcKD5rbTxx+m762UAaMfvuIZilAkSwghhBDHwRkUQgghxKkMsEjWSXCAQgghhDgUA1lmMx7CIxQOUPpAazxcLhl79nik4Y0/KNXLHRGZtr39qLTY9gVlLB2jtG5E3hqfT9fHqmz36aZmmMZ2hyzdLxw7q3m7NAf2Lz6tMZGaD60jSCalxiMalZ/1plK95//oXp9QvyTVeqnl0hqU0xGZTdWtnhu3Xz5XeWPUVwc5ds+CtZ7Lypq+d3T7X0i9mJ2+yHp78/5aj2R9Leb9E6ps/Un1sLXCP4+MxKb9hyjUoBBCCCHEcQzqAKW2thY33HADCgoKMHr0aNx+++14++23xTaGYaC6uhrhcBi5ubmYNm0aDh8+LLaJRqNYtmwZSkpKkJ+fj3nz5uHo0aMDeSmEEEJIvzPQTrJOYlAHKLt378a9996LV155BXV1dUgkEpg1axY6OrozjK5ZswZr167Fhg0bsH//foRCIcycORNtbd1TuVVVVdi5cyd27NiBvXv3or29HXPmzEEymezttIQQQsjQ4GyywGyWIcqgalBeeOEFUX7mmWcwevRo1NfX40tf+hIMw8C6deuwevVqzJ8/HwCwZcsWlJaWYvv27ViyZAlaWlrw9NNPY+vWrZgxYwYAYNu2bSgrK8NLL72EW2655Txrp63otQZFpkX3eotF2T9Sakw6jrWKcqeKlXvypAeAttBOpXJF2Wxn7+mxTntLwJJhG9sdMvStHTBbzcu4vp0vh52GRR8vkZD9VmtGtO9JV/sxUU7myudA1z/eKZ+bHI/131C+kdKHJNosNTAjxo0U5UCxfG4DgbGibNaQyGfB7ZbvAfu/8TLVg1nvnxmZaUzMmhJr3xJ9782+JvodJdtK65ns9Tu6LeiTcrHjKA1KS8uZl1tRUREAoKGhAZFIBLNmzUpv4/f7MXXqVOzbtw8AUF9fj3g8LrYJh8OorKxMb6OJRqNobW0VCyGEEOI4jH5YhiiOGaAYhoHly5fjpptuQmVlJYAz6ZwBmFI4l5aWptdFIhH4fD6MGjWqz200tbW1CAaD6aWsrKy/L4cQQgjJGmpQHMB9992HP/7xj/j3f/9307reMoDaTZVabfPwww+jpaUlvTQ2Np5/xQkhhBDS7zjCB2XZsmV4/vnnsWfPHowd2x0TDoXOpDqPRCIYM6Y7V01TU1N6ViUUCiEWi6G5uVnMojQ1NWHKlCm9ns/v98Pv12nfJeYcGCpnh8rB4fHIHCPeXFn2FUr/h9ORdlFOxa1zmtjHa5MW6/RATY9LtacKY7lDB9lvksmuPrY7i9YdaA2Lte9JLHZSnU9qPuJx6fdz+rjSWo2Qz4GRVDmsjikfk4Dsy9rHROcoScWsc2S53XmyPp4CVZYaFa2z0H5I5mfFTtMj15r3v3DPnl0uHP3OMa+XfUvfe+1xY85TJN+5+vhmzYlsLLOeSr+j1eqLRZOSrdB1CItkB3UGxTAM3Hffffj5z3+O3/72tygvLxfry8vLEQqFUFdXl/4tFoth9+7d6cHHxIkT4fV6xTbHjx/HoUOH+hygEEIIIUOBs8mMs1mGKoM6g3Lvvfdi+/bt+MUvfoGCgoK0ZiQYDCI3NxculwtVVVWoqalBRUUFKioqUFNTg7y8PCxcuDC97eLFi7FixQoUFxejqKgIK1euxIQJE9Jf9RBCCCFkaDGoA5SNGzcCAKZNmyZ+f+aZZ/Ctb30LALBq1Sp0dnZi6dKlaG5uxuTJk7Fr1y4UFHRPzT7xxBPweDxYsGABOjs7MX36dGzevBlut7ZOJoQQQoYQwzjEM6gDlHNRF7tcLlRXV6O6urrPbQKBANavX4/169f3Y+20f4GOd+qyV5WtmzYVl/HXRIeM38Zym0TZMGR82OORPik9Y+tut9bX5KhtZWzYPg5+kcRy+yTTB3iw26O7vmYfkna1rZ13RUKt1zoEub6zUwrK4x1Sh5Dskvu3Hzklyp4RSsOhct+c+uNHopzjt86lo9fnqPXxNvlcaR2FWecgnx23W2rJrHNg9VYe7L7SE2v9kdYbJZMdoqw9cOLxU6os84/pfGXmsmxbw9DvLWudoBkntXX/ke2XOEP5Kx5HiGQJIYQQ0gtMFkgIIYQQ4hw4g0IIIYQ4lDNf4mQT4unHygwwHKCcM9ZeIuZv+K0np7QGpeXtE6KcOC1j+b6RMv7rK5C5eXrmAvJ6pauu9iMw+6RY58xwVhy9P+g7t82ZsrU/g32OkP6tj1VOFK0b0DqAZFLG8bUGQ+9vp9HoPCFz7XQclWVNx1+kL4o7X2q1tCYl2Sk1LwVXyBxXuZcoTYjK1aM1Lb5CrSmRPid2fV97b2T6nDuLzK5N59LRbaV9T7ra5DvLH5SeMz5fiShrzxmzb8pwey/1gZGlSHYIj1CG0tNFCCGEkGECZ1AIIYQQp5Jtwr+hO4HCAQohhBDiVPiZMekF65tqzs0jY906x4d3hNQG+ItkfLbjiIzlt/zpY1H2FcvtAyUyfhso7tYO+Aukf4HPN1rVTe5rlz/E+T4h2T6AWuMRV+v1Y6I9cvTx7K5fa070+a3rk0x2x/4TiVaxLhbT/jnyWInOTlFOqVw1yajUgOj4dfQTuX/7e1Jj4nLLqHEkIvu9zyPb0qvMFH1epVHJkxqVvLDy0ghofZU8v+k59Bar9VoHIc+n76V+2btc1t9wGoa1p5B9Lp9MnqX+fk71tVvn5kkoD5zkCJVvTGlWAK0HyrSt7NpOM0w0KxcRHKAQQgghDoUzKIQQQghxHMPYp41f8RBCCCHEeXAGpQ90vFXH8vW41OwhIGPZOT4Vey+QmhVvUJY7jyvfk6CMted/SsbiA7nh7m19Ms7udkvPFK2XMY9TrX1BzHF4a18Q+9ixHdZTlOYpzMw0NHY5SXJytP9CpjoBO98V2Zd0rD6ZlPlu4vFuvVI0GhHrYm2y30RPSZ1AtFlqSLTGJNGu8rF0SU1KrFke7+T7J0U5zyf7fSIp2zZcLvVQgVLpa6LJGyM1JIERl4iyzyfL+l7ovFPaB8UuL5XW/+i+YueHZO+bkk3uHjv/HOuyzruk+525bO2hk+iQx0uM1BoUvb/cXr+X7JvC7volmWvFnAFDPIQQQghxHMN5gMIQDyGEEOJQjH5YMuXDDz/EN7/5TRQXFyMvLw/XXXcd6uvru+tkGKiurkY4HEZubi6mTZuGw4cPn/9F9gEHKIQQQggBADQ3N+MLX/gCvF4vfv3rX+PNN9/Ej3/8Y4wcOTK9zZo1a7B27Vps2LAB+/fvRygUwsyZM9GmQszZwhBPH5hzlOh4rYzlJxLyxsRiMrdOZ5Nc3/Gh9K9o/7P0i8jxydi0xpWj/R5G9Pr/QG+5eOS+hiF1BqlUZrpvc9zdLodGZrls7HPTWJftNSjKC0RpPvT1GIa+N9Z+DFpXYNYZWHvqaM8dl6vbK0T7eHg8Upvk8rwvyu3vS9+S6MfyWqNNstzWLsunTstyiyprXxO/8j3x5EmfkzylpXIH5PaBYnl9Xm+RKgdFWWu/rNruTFneC/2c676gy/peme+dR623Pr/uC2a6+1bmz4V8zhMJqRGJx+U7KZlsU+tl30nEpN9SokPrl9Q7MyA1LWY9j7Vmxk5jY5dDyzD0e0gdzqGalIEO8Tz++OMoKyvDM888k/7t0ksvFcdbt24dVq9ejfnz5wMAtmzZgtLSUmzfvh1Lliw577pqOINCCCGEOJSzA5RsFgBobW0VSzSqjfPO8Pzzz2PSpEn4+7//e4wePRrXX389nnrqqfT6hoYGRCIRzJo1K/2b3+/H1KlTsW/fvn69dg5QCCGEkIucsrIyBIPB9FJbW9vrdu+99x42btyIiooKvPjii7jnnnvwne98B//6r/8KAIhEznw5WFpaKvYrLS1Nr+svGOIhhBBCHEp/hXgaGxtRWNgdUvX7td3EGVKpFCZNmoSamhoAwPXXX4/Dhw9j48aNuPvuu9PbmT/JNzL8TN4eDlD6IJnsVGUZb+3pRQGYc6B0Ncv1bSr233pYalRyPHIyK+9SGVv3j5Q6Eq93pNw/p2euHjuNSZdar2PVOjas0ZoKGTc3+67IuDug4/DWZzP7sOjrsfOqsLseidYTmTU8HlW2ztdijm3b5Ryx9sboGVtPpbQ+RulpRsh+O/JK5bHyadk2OtdO/jGpS8j7k/Q9OalEcSUF0rfkElXWfj/eEVIzonNUaZ8T7WOi742dBsXOC8Ssy5DaMO074/WOEmWtATLrwTJ95fatw7B7LnQ5mZTPvUkn1/mBXN+q9EnKQ8eVI/tlXGtQovo5tdagmN9D1nmPzH5CKo9Uhu8Zp9JfyYwLCwvFAKUvxowZg6uuukr8duWVV+JnP/sZACAUCgE4M5MyZsyY9DZNTU2mWZVsYYiHEEIIIQCAL3zhC3j77bfFb++88w7Gjx8PACgvL0coFEJdXV16fSwWw+7duzFlypR+rQtnUAghhBCnkmWIR7tF2/Hd734XU6ZMQU1NDRYsWIBXX30VmzZtwqZNmwCcmdGtqqpCTU0NKioqUFFRgZqaGuTl5WHhwoXnX89e4ACFEEIIcSgD/ZnxDTfcgJ07d+Lhhx/GD3/4Q5SXl2PdunW4884709usWrUKnZ2dWLp0KZqbmzF58mTs2rULBSqkmy0coPSBjpfG46dEuavrQ1k+IWPXp4+pnCjKXyLVKeOl3pDyswjI+KlblXXsvafOQutjtD+C9njR+UZ0XN7sc2LtLaH1MeY4vHU+EnPcXcfSrT1o7PwfzHmUZP2134PbresvdQ7aN0af37797HxV9P7d63VT6rb2+0OybqNl3F97XWiiJ1UeIHWt8YRs278oFX+OCvwHPpL93F8sNSc6R5W9HkpibmutD7LuG/q515qU2GnZXu4CWX/zvfKp9Zl5/phFhz3XWz/X+jnRz1UiIfVJnR8rb6YPTsnye/K58Oh8YgXqWk3/Llp7vJjvjdaWWWu37HNoZeo5M3yZM2cO5syZ0+d6l8uF6upqVFdXX9B6cIBCCCGEOJT+EskORThAIYQQQhzKcE4WyAEKIYQQ4lA4QCEmtN+Czyfjl1qH4C6VsXftEaDzVLS9f0qUuxpl2T9ant/sKSA9DaLRj3us61Tbyrh6KqnLKpYdl3H/HLeM5eZ4ZezW41X5VJQPimHIazFjl2PDOtau9UHas0a3ldbYuN1SR5BQGh7tdeF2a52EtY+MOR+LT5W1X4OdRiXHYlulXVLXpq8lkZC6g1Rc9wVZ1pqSkfny3h77RPmGKI1Kol35jpxW+qeY9vLQOgStJ7LLs2Sdr0U/GzrXjtYjxVplX/Lny75oZ1Rl7/VhnX+m5/52/c6cN0h7xKh3mk976ijNSqu81hGflnmR8sJSIJlbLPua1qbpe6frm0pZ5+zS16M1N/Y5tKzfO2Tw4QCFEEIIcSjDeQZlUI3a9uzZg7lz5yIcDsPlcuG5554T6w3DQHV1NcLhMHJzczFt2jQcPnxYbBONRrFs2TKUlJQgPz8f8+bNw9GjRwfwKgghhJALQ6oflqHKoA5QOjo6cO2112LDhg29rl+zZg3Wrl2LDRs2YP/+/QiFQpg5cybaethrV1VVYefOndixYwf27t2L9vZ2zJkzB8lkZp8nEkIIIcQ5DGqIZ/bs2Zg9e3av6wzDwLp167B69WrMnz8fALBlyxaUlpZi+/btWLJkCVpaWvD0009j69atmDFjBgBg27ZtKCsrw0svvYRbbrnlvOumdQbaX0GjPQVyfKfk8fKkNiD3EqkN+OSojHWnlOYk2SXL8ZjcvqdGxTYnRlTlwFAalESnjPt7crUni+w2gRJ5fJ+vRO6v8pOYvS3sNCg6ti5jzfG4zCkSbZUakqTS1CTaZft4VD4Yl4p1B3JlbDyRkH1D10drZNxunS9G3nuPR2p2tHeGuS+iT7QmRed+scsjpI9tLuvtZTmmNCfaJyUVUxoQrXlJZqsLsNOkaA2KdS6eRKe8t7FTUoOSKJEanqRP9hWN1raZNSv6lax9XLrrq+tu9geSz4HHI/N7mab+UypPkXpPRE/H1Hp5b/V7Qfdz/Vzo+tnpqXTf1W1pGNoHRT93OieYtUeOUzQpDPE4kIaGBkQiEcyaNSv9m9/vx9SpU7Fv3z4AQH19PeLxuNgmHA6jsrIyvU1vRKNRtLa2ioUQQghxGmcHKNksQxXHDlAif3Wk1NkRS0tL0+sikQh8Ph9GjRrV5za9UVtbi2AwmF7Kysr6ufaEEEIIyQbHDlDO0lsqe/tP+ay3efjhh9HS0pJeGhsb+6WuhBBCSH8ynGdQHPuZcSh0JodIJBLBmDFj0r83NTWlZ1VCoRBisRiam5vFLEpTU5Nl2me/3w+/39/n+jPoWLt1Xge7QZMdiaR1bF5rUDo/lrHyrhPdse94i4qTd1rrAHS2y+jHMo6uPVlMfgchnSBKDyqV5sXkbaF2N42bdWzY2v8h2qx0A6o9Th+VIT1/iYyVB0rk9WqdgT5fLHZSrtc+Kx5rLxLDkNoA7V1i7ntWunzrl5G5n9oN9nVZeWOofhvX/dhG56Crk+PWz5W1hsZeRwDL9VoXkUopH5R26ZMSbVbrO2Tfj3k/hkTeO79/tGX9zPdeXl+yh4eROW+QzBMUi0ltlq6L1mykEkrrpd45p6PqfMrDxkjotpXPYTQq69d1SpZzPLJ+br+81x6f1m6NFGWtYdEalVRK+8BoPZJ1zq7BYjhb3Tt2BqW8vByhUAh1dXXp32KxGHbv3p0efEycOBFer1dsc/z4cRw6dMhygEIIIYQQZzOoMyjt7e1499130+WGhgYcOHAARUVFGDduHKqqqlBTU4OKigpUVFSgpqYGeXl5WLhwIQAgGAxi8eLFWLFiBYqLi1FUVISVK1diwoQJ6a96CCGEkKHKcP6KZ1AHKK+99hpuvvnmdHn58uUAgEWLFmHz5s1YtWoVOjs7sXTpUjQ3N2Py5MnYtWsXCgq6p1WfeOIJeDweLFiwAJ2dnZg+fTo2b94Mt5uptAkhhAxtDGQ3yBi6w5NBHqBMmzbNsuFdLheqq6tRXV3d5zaBQADr16/H+vXr+7Vu2lPAHKu3jpXreKrOzWMov4eYit0nlRdJ50dSc6K9Szo/7I7nxk6q2G9ceaqoc+nBnNYVjA6NEGXfSOkv4PHJstZMmHIBqbi/foTsvDrMOgLZVlqD0nlMxrqjTVJHENDXN0rnEpLn6+o8Lo9vc2+8I6QGJhmU5/f5Zf29XpnjRHuZWD0zZq8HiV1eIJfqtzkea41KKiWvVefeSar1UM+BPp/LrfU21vljtI4imdR9S+surHNaJeJKx9GqNCgn5L2LfqLyxyjPHa2jcI/SOgqpYbHTZ/XMV6NzUMViTbLcJq8FSiqmz6W1aokOub69S/VjpUFJKt+UREx6Q3Wptmp9R2pkvAVSF+gtlGVfUJ7fWyCP7wuMlOuV1iuV0jm0rPuObp/BImUYZi1XhvsPVRyrQSGEEELI8MWxX/EQQgghw53h/BUPByiEEEKIQ6FIlpiIxaSfgc4LYeuTYorl6zwT1jlMOj5UuomIjLcayk8icupU+v/bVKxYe1NojUm+8oQJ+GTcv1j5IZiTwciy1pyY9Twq54fyHzDnntFtL/fX3hGBS6T/QUxrUk7JckHcOrFkIiY1Jqcj8t60Nci8SJpAscpJovwmUCy1A+acJLIs21dqPHJydNvJe6nR+Uq0ZiLHK8/tybGOCuvcO/rlaPIPUpoTfT79XCWTSgeh+pb2qDG3pbw+U19V+WXibXJ9TOmXTn+oPXJE0eQZ5BshPXN0niqtp9L11Roaua9sm6jJD0keW7e19j1J2mhQEu1Kg6KOn1C5e+Jaz6P8lrxBpW1TOcD8I6U2zOeXLuNer2xLjc79k0i0qrJ6xzpEgzKc4QCFEEIIcSicQSGEEEKI88jWrn4ID1D4FQ8hhBBCHAdnUPrg9Om/iLKO7efkyHiojpXnqNi6O1fF9gMy/qtj9+8el14b2l8iV+lE3m/q1jFoPYs+dldcxlYvKZSx25DKDp1ScXmdTyWVVLHlqMwknYjJ2LX2gDGUV4bW6+hYufbuMHlHqPolTsv6t5yWse9RrTJWrmPnun5dJ2Ssuu1P0s9B5y7KLZU+K94Rsi95vDIXj9bY6Jwq8Xh32UhJzUWO0u/ofCS67TTugPLtyJM6gBy/tQFiVPWtRErnTVJ/zSk5k753PX0/zpQ7VFlrUOS91r4wbre170iiS5bjbbIvtDUrncKbsr456jn3F8n3RDIun5W4R+VxUp47Zk+gnD63TXTJusZOyedOa7G05iOh9DbxZrm/1qBETypPGHV83XcSHbJ+2o/I5VP+SdpTRvVNT7HSxCS0543UnOjcQLEupUkx5Ray1qYNFMPZB4UDFEIIIcShDGcNCkM8hBBCCHEcnEEhhBBCHAqN2oiJtvc/EWWXV0426W/0tX+Ezr3jHSFj//njR4ry5e1l8vwnpffGm0ePivKJVvkNv8/Tff7RQalpOKlycrSdlHHvS9T2nwmFRFnnqtE6BK3ZSHTJ2LL2P0glVC4dFWvWmhO3T5VVLDpHrdc5QnRs/VSH1BHETlrXV/tD6HwsLS3yeKFPSZ1Djuo7WmcRj56Sx0/IvpeKqfP3iPXrttT90JPbbLlee9qkYirubvN2S6rpY61/MvmguLV/kMqto64n2i77qqmvqXut20PnxPIWyOdGt4e+11p/1dopdQxHTkh90KeKZB6lEaeknivWIvdPdh2TZdX++j3S89nQddeaCb1vV5Psp9qHJNEu21Y/J61Ku6W1XCPUc+HJl+887SnT2iaP3/i6bMtLLxsjyl51PFNfhuwr2tNGl2Mtsj5aI6PfS4PFcA7xcIBCCCGEOJThLJKlBoUQQgghjoMDFEIIIcShnA3xZLOcL7W1tXC5XKiqqhL1qa6uRjgcRm5uLqZNm4bDhw/3w5WaYYinD069KfOjaJ2DW2lQPPmqrNbr3Dx6/YgKGbv2XSLzt0zKk54FnZ0yftrWIzbeMy8PALzZ2CjK2kNF5+LJHyXPrT1cdKxax3J1nF37KWidgy67tAZFaV60v4JHlZNRebzYx7I+WpOj10c/kWWtIen6SMbOP1Z6oKLmkfL4p6yvX+tAksqLwxQ7P9lDgxJXGpSAdT81+ZqY2lrlqrGJw+uXX2cs1seWvaOP39nUbrle97Xkaa0bUBoUde+094e3QPZ9jdZxFORKX5Oe/kMAcOno0ao++vpk39F9QesgoDQ7bl/3/dG6Np3XSOtzdD8/1SrbWj8XTS3SJ0Tn+Dqm3jNFjdJPyaPqp6/V75V98aRqy7FtxaIca5Xn7ziqfEyUfkjrlZLKDyl+SucWsu5Lg8VgaVD279+PTZs24ZprrhG/r1mzBmvXrsXmzZtx+eWX49FHH8XMmTPx9ttvo6CgoI+jnR+cQSGEEEJImvb2dtx555146qmnMKqHcadhGFi3bh1Wr16N+fPno7KyElu2bMHp06exffv2fq8HByiEEEKIQ+mvEE9ra6tYotFon+e89957ceutt2LGjBni94aGBkQiEcyaNSv9m9/vx9SpU7Fv375+v3YOUAghhBCHkuqHBQDKysoQDAbTS21tba/n27FjB15//fVe10ciZ9KYlJaWit9LS0vT6/oTalD64KM3ZWO7c+RYzqfipz6tSSlQ3+yr2L8eGurcPTr/i69Yxs5zOqR2INVj+0RSxnpPq5HytZdeKsrFI6TPidbLaL8E7aeQ7FCx7hMy1t3RpfwGVP3iquxRbe1xy2vVsWu/ym2jfVq0n8MJFWvXfg7+D6WmBEqH0KFy8XykYvFjT8jYubtBPWY6l5H2helU/g2q3Ha6hwZFxZe9qq10P/X6VK4dpS/yBK01GdonJKn6qe5r+l7r/Ts/bLMsm/ROSjfQqXL/6Pro5zY/IK/PWySfK5Mvitb/qPOdbJc6jk51/Z3H5PWYvEeU505cXW+Oqn/PfDVeda+05kNrUnwlUj9TrO59wGedp6lBaUT+6513RDmYJ7Vrl+YpbyhVH61X0t5Obe3qufxAak46j8q2NfkXtct7pfOZnVaaGp2jTHv8DHUaGxtR2CPvmt9vftYbGxtx//33Y9euXQgEAqb1ZzH5FxmG6bf+gAMUQgghxKH0l0i2sLBQDFB6o76+Hk1NTZg4cWL6t2QyiT179mDDhg14++23AZyZSRkzpttIr6mpyTSr0h8wxEMIIYQ4lIH8zHj69Ok4ePAgDhw4kF4mTZqEO++8EwcOHMBll12GUCiEurq69D6xWAy7d+/GlClT+v3aOYNCCCGEEBQUFKCyslL8lp+fj+Li4vTvVVVVqKmpQUVFBSoqKlBTU4O8vDwsXLiw3+vDAUofvH1M5sjQsWyvRzadT5UDWqOi1ps0LWq9juelVGw9ocrRHvHTqIqldqhYqz52u4qbJxtlTgzPh7KuWjOiY7f6fDrWbKdB0W1jp0HRvi56fZvSmGjNyPFmma9Gt59urzaVj+XDT2TunNKRI0W5QOf+0boMm/aIq9h5R4/7pTUXuh/ptrNbrz1xdD/XdelS91a3tV7/kfLW8Kq20X1D950O1Vd1ffT+WsORp/pKXpO8Xt139L3/ROmXdF/SHkRa96D/mtXXo58lff6e+qw8da/0O8d0LZDottH6Gn3uMaNkXqH9774rylpD4j+qcgWpa29W914/h8dU+bTqC/qdqNfrtrTruzqPlFMs4p1mdb9q1Sp0dnZi6dKlaG5uxuTJk7Fr165+90ABOEAhhBBCHMtgJwt8+eWXRdnlcqG6uhrV1dVZHfdc4ACFEEIIcSpZDlDgkJmg84EiWUIIIYQ4Ds6g9MGhI0dEOUfHglXs3q1j/aqsY/k6vqv9K/T5NDqu2HOE3a40IK0q1qvj5Nq7wu57dq2R0LFru9ivrrvWDei2tWt7ky+KKuvzHfn4Y1EeHQyKsr4XGq3/eU8ZFGltgK6/SZegjqc1KjrW3lMjo6/NzkNG91NdN63nsWsLu76m9Tnac0aj+5bWA+my7ju6PfT16edQ6zb0eq2H0vX78ORJUf5zxNo/SWuG9PXoZ0XT8/7pd4z2vPFn+M7RddOajFalvTqh3iPvq+dKa1L08XVb6udyZH6+KNvdC7u21H1DX58+nn7uBgunaVAGkotmBuXJJ59EeXk5AoEAJk6ciN///veDXSVCCCEkKwxk+anxYF9AFlwUA5Rnn30WVVVVWL16Nd544w188YtfxOzZs3FEzYIQQgghZGhwUQxQ1q5di8WLF+O///f/jiuvvBLr1q1DWVkZNm7cONhVI4QQQs6bgTRqcxpDXoMSi8VQX1+Phx56SPw+a9asjLMrjkX3iO2lizhE9Avl8TLcefev9s39xcGDB/v1eEOZncOsr7391luDXYVB4/9+8EG/Hu/Nw4f79XjZMq7H/6cAHB2g8w5nDcqQH6CcOHECyWQyo+yK0WhUpJpuVWIuQgghhAwuF0WIB8gsu2Jtba1IO11WVjYQVSSEEEIygiGeIUxJSQncbrdptsQqu+LDDz+M5cuXp8stLS0YN24cnPFRGSGEECdz9t+KgfjHnyGeIYzP58PEiRNRV1eHr371q+nf6+rqcNttt/W6j9/vh7+HV8XZEM/wipYTQgjJhra2NgSVjxLpP4b8AAUAli9fjrvuuguTJk3CjTfeiE2bNuHIkSO45557zmn/cDiMxsZGGIaBcePGobGxEYWFhRe41hcXra2tKCsrY9udJ2y/84dtd/6w7c4PwzDQ1taGcDg8IOcazFw8g8lFMUD5+te/jpMnT+KHP/whjh8/jsrKSvzqV7/C+PHjz2n/nJwcjB07Nj2TUlhYyIf1PGHbZQfb7/xh250/bLvMGaiZE4Z4LgKWLl2KpUuXDnY1CCGEkH5jOM+gXDRf8RBCCCHk4uGimUHpD/x+Px555BEhoCXnBtsuO9h+5w/b7vxh2zmf4RzicRlDef6HEEIIuQhpbW1FMBjEw3PnmrJuZ0JXPI7aX/4SLS0tQ05nxBAPIYQQQhwHQzyEEEKIQzGyDPEM5SAJByiEEEKIQ+FXPIQQQgghDoIDlL/y5JNPory8HIFAABMnTsTvf//7wa7SoFNbW4sbbrgBBQUFGD16NG6//Xa8/fbbYhvDMFBdXY1wOIzc3FxMmzYNh1Wa9Gg0imXLlqGkpAT5+fmYN28ejh4dqGTlzqC2thYulwtVVVXp39h21nz44Yf45je/ieLiYuTl5eG6665DfX19ej3br3cSiQT+8R//EeXl5cjNzcVll12GH/7wh0ilurONse2GDme/4slmGapwgALg2WefRVVVFVavXo033ngDX/ziFzF79mwcOXJksKs2qOzevRv33nsvXnnlFdTV1SGRSGDWrFno6OhIb7NmzRqsXbsWGzZswP79+xEKhTBz5ky0tbWlt6mqqsLOnTuxY8cO7N27F+3t7ZgzZw6SyeRgXNaAs3//fmzatAnXXHON+J1t1zfNzc34whe+AK/Xi1//+td488038eMf/xgjR45Mb8P2653HH38cP/3pT7Fhwwa89dZbWLNmDf7pn/4J69evT2/Dths6DOdsxjCI8Td/8zfGPffcI3674oorjIceemiQauRMmpqaDADG7t27DcMwjFQqZYRCIeNHP/pRepuuri4jGAwaP/3pTw3DMIxTp04ZXq/X2LFjR3qbDz/80MjJyTFeeOGFgb2AQaCtrc2oqKgw6urqjKlTpxr333+/YRhsOzsefPBB46abbupzPduvb2699Vbj29/+tvht/vz5xje/+U3DMNh2Q4WWlhYDgLFy9mxj9bx5572snD3bAGC0tLQM9iVlzLCfQYnFYqivr8esWbPE77NmzcK+ffsGqVbOpKWlBQBQVFQEAGhoaEAkEhFt5/f7MXXq1HTb1dfXIx6Pi23C4TAqKyuHRfvee++9uPXWWzFjxgzxO9vOmueffx6TJk3C3//932P06NG4/vrr8dRTT6XXs/365qabbsJvfvMbvPPOOwCAP/zhD9i7dy++8pWvAGDbDTWGc4hn2H/Fc+LECSSTSZSWlorfS0tLEYlEBqlWzsMwDCxfvhw33XQTKisrASDdPr213QcffJDexufzYdSoUaZtLvb23bFjB15//XXs37/ftI5tZ817772HjRs3Yvny5fiHf/gHvPrqq/jOd74Dv9+Pu+++m+1nwYMPPoiWlhZcccUVcLvdSCaTeOyxx/CNb3wDAPveUGM4O8kO+wHKWVwulygbhmH6bThz33334Y9//CP27t1rWnc+bXext29jYyPuv/9+7Nq1C4FAoM/t2Ha9k0qlMGnSJNTU1AAArr/+ehw+fBgbN27E3Xffnd6O7Wfm2WefxbZt27B9+3ZcffXVOHDgAKqqqhAOh7Fo0aL0dmy7oYHBz4yHLyUlJXC73aa/Cpqamkx/YQxXli1bhueffx6/+93vMHbs2PTvoVAIACzbLhQKIRaLobm5uc9tLkbq6+vR1NSEiRMnwuPxwOPxYPfu3fhf/+t/wePxpK+dbdc7Y8aMwVVXXSV+u/LKK9PCdfa9vnnggQfw0EMP4Y477sCECRNw11134bvf/S5qa2sBsO3I0GHYD1B8Ph8mTpyIuro68XtdXR2mTJkySLVyBoZh4L777sPPf/5z/Pa3v0V5eblYX15ejlAoJNouFoth9+7d6babOHEivF6v2Ob48eM4dOjQRd2+06dPx8GDB3HgwIH0MmnSJNx55504cOAALrvsMradBV/4whdMn7S/8847GD9+PAD2PStOnz6NnBz5ane73enPjNl2QwtqUIY5y5cvx1133YVJkybhxhtvxKZNm3DkyBHcc889g121QeXee+/F9u3b8Ytf/AIFBQXpv7iCwSByc3PTvh41NTWoqKhARUUFampqkJeXh4ULF6a3Xbx4MVasWIHi4mIUFRVh5cqVmDBhgkk4ejFRUFCQ1uqcJT8/H8XFxenf2XZ9893vfhdTpkxBTU0NFixYgFdffRWbNm3Cpk2bAIB9z4K5c+fisccew7hx43D11VfjjTfewNq1a/Htb38bANtuqDGcNSj8zPiv/OQnPzHGjx9v+Hw+43Of+1z6U9rhDIBel2eeeSa9TSqVMh555BEjFAoZfr/f+NKXvmQcPHhQHKezs9O47777jKKiIiM3N9eYM2eOceTIkQG+msGn52fGhsG2s+OXv/ylUVlZafj9fuOKK64wNm3aJNaz/XqntbXVuP/++41x48YZgUDAuOyyy4zVq1cb0Wg0vQ3bzvmc/cz43hkzjOWzZ5/3cu+MGRl9ZlxTU2NMmjTJGDFihHHJJZcYt912m/GnP/1JbHO2/4wZM8YIBALG1KlTjUOHDvV7G7gMYygPrwghhJCLj9bWVgSDQSydPh1+z/kHO6KJBJ78zW/Q0tKCwsJC2+2//OUv44477sANN9yARCKB1atX4+DBg3jzzTeRn58P4IwZ4GOPPYbNmzfj8ssvx6OPPoo9e/bg7bffRkFBwXnXVcMBCiGEEOIwzg5Q7vnbv816gPLT3/72nAcomo8//hijR4/G7t278aUvfQmGYSAcDqOqqgoPPvjgmXNEoygtLcXjjz+OJUuWnHddNcNeJEsIIYRc7LS2toolGo2e037nY9DZX3CAQgghhDgUI8sveM4GScrKyhAMBtPL2c/O7c6diUFnf5v48SseQgghxKH011c8jY2NIsTj9/tt9+1vg85M4QCFEEIIucgpLCzMSINy1qBzz549fRp0jhkzJv37hTDxY4iHEEIIcSipVCrrJROMfjDo7C84g0IIIYQ4FGOAc/H0h0Fnf8EBCiGEEOJQBtpJduPGjQCAadOmid+feeYZfOtb3wIArFq1Cp2dnVi6dCmam5sxefJk7Nq1q189UAD6oBBCepBMJvHFL34RY8aMwc9+9rP07y0tLaisrMSiRYvw6KOPDmINCRkenPVB+dZNN8GXhQ9KLJHA5r17z9sHZTChBoUQksbtdmPLli144YUX8G//9m/p35ctW4aioiJ8//vfH8TaETL8GGgNipNgiIcQIqioqEBtbS2WLVuGm2++Gfv378eOHTvw6quvwufzDXb1CBlWDOdkgRygEEJMLFu2DDt37sTdd9+NgwcP4vvf/z6uu+66wa4WIWQYwQEKIcSEy+XCxo0bceWVV2LChAl46KGHBrtKhAxLhvMMCjUohJBe+Zd/+Rfk5eWhoaEBR48eHezqEDIsGc4aFA5QCCEm/vM//xNPPPEEfvGLX+DGG2/E4sWLs/JiIISQTOEAhRAi6OzsxKJFi7BkyRLMmDED/+f//B/s378f//t//+/Brhohw46UYSCZxcIQDyHkouGhhx5CKpXC448/DgAYN24cfvzjH+OBBx7A+++/P7iVI2SYwRAPIYQA2L17N37yk59g8+bNyM/PT//+P/7H/8CUKVMY6iGEDBj8iocQkmbq1KlIJBK9rnvxxRcHuDaEkBSy+xJn6M6fcIBCCCGEOJZUKoWUy5XV/kMVDlAIIYQQh5JMpZCTxQAlOYQHKNSgEEIIIcRxcAaFEEIIcSjD2UmWAxRCCCHEoTDEQwghhBDiIDiDQgghhDgUfsVDCCGEEMeRSqWQHKYDFIZ4CCGEEOI4OINCCCGEOJRkKgXXMBXJcoBCCCGEOJThPEBhiIcQQgghjoMzKIQQQohD4Vc8hBBCCHEciSwHGNnuP5hwgEIIIYQ4lGQyCVcWdvXUoBBCCCGE9COcQSGEEEIcSrYzIEN5BoUDFEIIIcShJJNJgCEeQgghhBBnwBkUQgghxKHEo1Gkcs5/LoEzKIQQQgjpd+LxOOKx2Pkv8fh5nffJJ59EeXk5AoEAJk6ciN///vf9fGX2cIBCCCGEkDTPPvssqqqqsHr1arzxxhv44he/iNmzZ+PIkSMDWg+XYWShviGEEEJIv9Pa2opgMIhyrxc52TjJGgYa4nG0tLSgsLDwnPaZPHkyPve5z2Hjxo3p36688krcfvvtqK2tPe+6ZAo1KIQQQohDicfjWYU6MlWgxGIx1NfX46GHHhK/z5o1C/v27cuiJpnDAQohhBDiULKVuJ7dv7W1Vfzu9/vh9/tN2584cQLJZBKlpaXi99LSUkQikSxrkxkcoBBCCCEOw+fzIRQK4Vg/DApGjBiBsrIy8dsjjzyC6urqPvdxqbCSYRim3y40HKAQQgghDiMQCKChoQGxWCzrY/U2uOht9gQASkpK4Ha7TbMlTU1NplmVCw0HKIQQQogDCQQCCAQCA3pOn8+HiRMnoq6uDl/96lfTv9fV1eG2224b0LpwgEIIIYSQNMuXL8ddd92FSZMm4cYbb8SmTZtw5MgR3HPPPQNaDw5QCCGEEJLm61//Ok6ePIkf/vCHOH78OCorK/GrX/0K48ePH9B60AeFEEIIIY6DTrKEEEIIcRwcoBBCCCHEcXCAQgghhBDHwQEKIYQQQhwHByiEEEIIcRwcoBBCCCHEcXCAQgghhBDHwQEKIYQQQhwHByiEEEIIcRwcoBBCCCHEcXCAQgghhBDHwQEKIYQQQhzH/wekSIldlDNRDQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# PARAMS\n",
    "K = 0.0005  # stream power coefficient, bigger = streams erode more quickly\n",
    "U = 0.0004  # uplift rate in meters per year\n",
    "\n",
    "dt = 1000  # time step in years\n",
    "dx = 10  # space step in meters\n",
    "\n",
    "nr = 60  # number of model rows\n",
    "nc = 100  # number of model columns\n",
    "\n",
    "# Original K_sp value is 1e-5\n",
    "K_sp = 1.0e-5  # units vary depending on m_sp and n_sp\n",
    "m_sp = 0.5  # exponent on drainage area in stream power equation\n",
    "n_sp = 1.0  # exponent on slope in stream power equation\n",
    "\n",
    "D = 0.001  # Hillslope diffusivity and m2 per years\n",
    "\n",
    "\n",
    "#create grid\n",
    "grid = RasterModelGrid((nr, nc), xy_spacing=dx)\n",
    "grid.set_closed_boundaries_at_grid_edges(True, False, True, False)\n",
    "# add a topographic__elevation field with noise\n",
    "z = grid.add_zeros(\"topographic__elevation\", at=\"node\")\n",
    "z[grid.core_nodes] += 100.0 + np.random.randn(grid.core_nodes.size)\n",
    "\n",
    "#imshow_grid(grid, \"topographic__elevation\")\n",
    "\n",
    "# initializing flow accumulator\n",
    "fr = FlowAccumulator(grid, depression_finder=\"DepressionFinderAndRouter\", routing=\"D4\")\n",
    "\n",
    "# initializing stream power incision\n",
    "spr = StreamPowerEroder(\n",
    "    grid, K_sp=K_sp, m_sp=m_sp, n_sp=n_sp, threshold_sp=0.0\n",
    ")  \n",
    "\n",
    "# initializing linear diffuser\n",
    "dd = LinearDiffuser(\n",
    "    grid,\n",
    "    linear_diffusivity=D,\n",
    ")\n",
    "\n",
    "# Run this model for multiple iterations\n",
    "for i in range(200):\n",
    "    # Flow Accumulator\n",
    "    fr.run_one_step()\n",
    "    # Stream Power\n",
    "    spr.run_one_step(dt)\n",
    "    # linear Diffuser\n",
    "    dd.run_one_step(dt)\n",
    "\n",
    "    #z[grid.core_nodes] += U * dt\n",
    "\n",
    "# plot the final topography\n",
    "imshow_grid(grid, \"topographic__elevation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b961304-3d98-40ae-bb82-ac882a8ef30b",
   "metadata": {},
   "source": [
    "Vary diffusivity spatially"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "320fc1d4-76fc-4a43-ae75-fddd0bc23fde",
   "metadata": {},
   "outputs": [],
   "source": [
    "# PARAMS\n",
    "K = 0.0005  # stream power coefficient, bigger = streams erode more quickly\n",
    "U = 0.0004  # uplift rate in meters per year\n",
    "\n",
    "dt = 1000  # time step in years\n",
    "dx = 10  # space step in meters\n",
    "\n",
    "nr = 60  # number of model rows\n",
    "nc = 100  # number of model columns\n",
    "\n",
    "# Original K_sp value is 1e-5\n",
    "K_sp = 1.0e-5  # units vary depending on m_sp and n_sp\n",
    "m_sp = 0.5  # exponent on drainage area in stream power equation\n",
    "n_sp = 1.0  # exponent on slope in stream power equation\n",
    "\n",
    "D = 0.001  # Hillslope diffusivity and m2 per years\n",
    "\n",
    "\n",
    "# rate of EQs\n",
    "rate = np.array([0, 0, 0.01, 0.01, 0, 0, 0.01, 0.01, 0, 0, 0.01, 0.01, 0])\n",
    "\n",
    "\n",
    "#create grid\n",
    "grid = RasterModelGrid((nr, nc), xy_spacing=dx)\n",
    "grid.set_closed_boundaries_at_grid_edges(True, False, True, False)\n",
    "# add a topographic__elevation field with noise\n",
    "z = grid.add_zeros(\"topographic__elevation\", at=\"node\")\n",
    "z[grid.core_nodes] += 100.0 + np.random.randn(grid.core_nodes.size)\n",
    "\n",
    "#imshow_grid(grid, \"topographic__elevation\")\n",
    "\n",
    "# initializing flow accumulator\n",
    "fr = FlowAccumulator(grid, depression_finder=\"DepressionFinderAndRouter\", routing=\"D4\")\n",
    "\n",
    "# initializing stream power incision\n",
    "spr = StreamPowerEroder(\n",
    "    grid, K_sp=K_sp, m_sp=m_sp, n_sp=n_sp, threshold_sp=0.0\n",
    ")  \n",
    "\n",
    "# initializing linear diffuser\n",
    "dd = LinearDiffuser(\n",
    "    grid,\n",
    "    linear_diffusivity=D,\n",
    ")\n",
    "\n",
    "# Run this model for multiple iterations\n",
    "for i in range(200):\n",
    "    # Flow Accumulator\n",
    "    fr.run_one_step()\n",
    "    # Stream Power\n",
    "    spr.run_one_step(dt)\n",
    "    # linear Diffuser\n",
    "    dd.run_one_step(dt)\n",
    "\n",
    "    #z[grid.core_nodes] += U * dt\n",
    "\n",
    "# plot the final topography\n",
    "imshow_grid(grid, \"topographic__elevation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "159ffba0-bf1f-402f-9b46-b6f920d16d1a",
   "metadata": {},
   "source": [
    "# Linear Diffuser example with varying k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "4cf1916a-6c64-44e0-9760-76a53ffb4e3b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAG2CAYAAABxpo8aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAnl0lEQVR4nO3df3BU9bnH8c8SyAY0iQNINilJum21ZQzFkthIEH9QiMYOI0JruHQk9LaOGcBOTLGS0kpkHFKZ6tzbxkRQsdKil2kVyh1TNQ4SsOBtjEQRaMut0QRNzCXVTUghMcm5f8SsrEnIht3N2d3v+9U5M+bknLPPbhcenud7vt/jsCzLEgAAMMY4uwMAAABji+QPAIBhSP4AABiG5A8AgGFI/gAAGIbkDwCAYUj+AAAYhuQPAIBhSP4AABiG5A8AgGFI/gAA2GT//v1atGiRUlJS5HA4tHv37hHPqampUWZmpuLi4vSlL31Jjz766Khfl+QPAIBNOjs7NWvWLJWXl/t1fENDg26++WbNmzdPhw8f1k9/+lP96Ec/0rPPPjuq13XwYB8AAOzncDi0a9cuLV68eNhj7r33Xu3Zs0fHjx/37issLNSbb76pQ4cO+f1a4wMJNBL09fXpgw8+UHx8vBwOh93hAABGybIsdXR0KCUlRePGhaZhffbsWXV3dwflWpZlDco3TqdTTqcz4GsfOnRIubm5PvtuvPFGPfHEE/rkk080YcIEv64T9cn/gw8+UGpqqt1hAAAC1NTUpOnTpwf9umfPnpXb7VZLS0tQrnfxxRfr9OnTPvs2bNig0tLSgK/d0tKipKQkn31JSUnq6enRqVOnlJyc7Nd1oj75x8fH2x0CACAIQvX3eXd3t1paWtTU1KSEhISArtXe3q7U1NRB1wpG1T/g812FgdH70XS3oz750+oHgOgQ6r/PExLilZAQ6D8wrE+vlRDwPySG4nK5BnUoWltbNX78eE2ZMsXv60R98gcAwB+WZSnQe+BDfQ/9nDlz9N///d8++1566SVlZWX5Pd4vMdUPAIBPWUHa/Hf69GnV19ervr5eUv9Uvvr6ejU2NkqSSkpKtGLFCu/xhYWFeu+991RcXKzjx49r27ZteuKJJ7R27dpRvtUo5/F4gvX/JhsbGxubjZvH4wlpnvj4439afX09AW0ff/zPUcX6yiuvDPleCwoKLMuyrIKCAuu6667zOWffvn3WN77xDSs2Ntb64he/aFVWVo76PUf9PP/29nYlJibaHQYAIEAejyck4+gDeeLjj9uCcsPfJZdMCVmswcKYPwAAUkSM+QcLY/4AABiGyh8AAEmfDbkHeo3wR/IHAECSScmftj8AAIah8gcAQGbd8EfyBwBAEm1/AAAQtaj8R5BmdwAAECUa7Q5gROZU/iR/AADEmD8AAAYyp/JnzB8AAMNQ+QMAIMmkyp/kDwCAzBrzp+0PAIBhqPwBAJBE2x8AAOOYk/xp+wMAYBgqfwAAZNYNfyR/AAAk0fYHAABRi8ofAABJ/VV7XxCuEf5I/gAAiDF/AAAMxJg/AACIUlT+AABIMqnyJ/kDACCzxvxp+wMAYBgqfwAAJNH2BwDAOOYkf9r+AAAYhsofAACZdcMfyR8AAEm0/QEAQNSi8gcAQLT9AQAwUJ8Cf6pfoOePDZI/AACSGPMHAABRi8ofAAAx5g8AgIFo+wMAgChF5Q8AgCSTKn+SPwAAMmvMn7Y/AACGofIHAMArMir3QJH8AQCQZNIKf7T9AQAwDJU/AAAy64Y/kj8AAJKY6gcAgHHMSf6M+QMAYBgqfwAAxJg/AAAGou0PAACiFJU/AACSTFrkh+QPAIDMGvO3te1fVlamq666SvHx8Zo2bZoWL16sv/3tbz7HWJal0tJSpaSkaOLEibr++ut19OhRmyIGACDy2Zr8a2pqtHr1ar322muqrq5WT0+PcnNz1dnZ6T1m8+bNevjhh1VeXq7a2lq5XC4tXLhQHR0dNkYOAIg+VpC28OewwqhH8X//93+aNm2aampqdO2118qyLKWkpKioqEj33nuvJKmrq0tJSUl68MEHdeedd454zfb2diUmJl5wTGkXfCYA4FyNAZ7v8XiUkJAQlFjONZAnGhv3KyHh4gCvdVppadeGLNZgCau7/T0ejyRp8uTJkqSGhga1tLQoNzfXe4zT6dR1112ngwcPDnmNrq4utbe3+2wAAOAzYZP8LctScXGxrrnmGmVkZEiSWlpaJElJSUk+xyYlJXl/93llZWVKTEz0bqmpqaENHAAQFSyrLyhbJAib5L9mzRq99dZbeuaZZwb9zuFw+PxsWdagfQNKSkrk8Xi8W1NTU0jiBQBEG3PG/MNiqt9dd92lPXv2aP/+/Zo+fbp3v8vlktTfAUhOTvbub21tHdQNGOB0OuV0OkMbMAAg6jDVb4xYlqU1a9boueee0969e+V2u31+73a75XK5VF1d7d3X3d2tmpoa5eTkjHW4AACEREVFhdxut+Li4pSZmakDBw6c9/gdO3Zo1qxZmjRpkpKTk/X9739fbW1tfr+ercl/9erV+t3vfqenn35a8fHxamlpUUtLi86cOSOpv91fVFSkTZs2adeuXXr77be1cuVKTZo0ScuXL7czdABA1LGn7b9z504VFRVp/fr1Onz4sObNm6e8vDw1Ng49P+LVV1/VihUr9IMf/EBHjx7V73//e9XW1uqHP/yh369p61S/4cbtn3zySa1cuVJSf3fg/vvv15YtW/TRRx8pOztbjzzyiPemwJEw1Q8AwkO4T/V7992XlZBwUYDX6tQXv7hgVLFmZ2dr9uzZqqys9O6bMWOGFi9erLKyskHH//KXv1RlZaX+8Y9/ePf9+te/1ubNm/2+z832tv9Q20Dil/r/gVBaWqrm5madPXtWNTU1fid+AADs8Pkp511dXUMe193drbq6Op8p7ZKUm5s77JT2nJwcnTx5UlVVVbIsSx9++KH+8Ic/6Nvf/rbf8YXN3f4AANhpuIJ0tJskpaam+kw7H6qCl6RTp06pt7d3VFPac3JytGPHDuXn5ys2NlYul0uXXHKJfv3rX/v9Xkn+AABICuaYf1NTk8+085KSkvO+8mimtB87dkw/+tGPdN9996murk4vvPCCGhoaVFhY6Pc7DYupfgAARJOEhAS/xvynTp2qmJiYQVX++aa0l5WVae7cubrnnnskSV//+td10UUXad68eXrggQd8psYPh8ofAABJdtztHxsbq8zMTJ8p7ZJUXV097JT2f/3rXxo3zjd9x8TE9L8DP+/hp/IHAECSFIzleUd/fnFxsW6//XZlZWVpzpw52rp1qxobG71t/JKSEr3//vvavn27JGnRokW64447VFlZqRtvvFHNzc0qKirSN7/5TaWkpPj1miR/AABslJ+fr7a2Nm3cuFHNzc3KyMhQVVWV0tPTJUnNzc0+c/5Xrlypjo4OlZeX68c//rEuueQSzZ8/Xw8++KDfrxlWj/QNBeb5A0B4CPd5/g0NVYqPD2yef0dHp9zum8P+kb5U/gAASArOg3kio54m+QMAoM8e6RvoNSIBd/sDAGAYKn8AACTR9gcAwDjmJH/a/gAAGIbKHwAAmXXDH8kfAABJtP0BAEDUovIHAED9D8UJdNHbSFk0l+QPAICk/pZ9oGP2kZH8afsDAGAYKn8AACSZdMMfyR8AADHmDwCAeSyrfwv0GhGAMX8AAAxD8gcAwDC0/QEAkFlj/lT+AAAYhsofAADJpJl+JH8AACTa/gAAIIpR+QMAINH2BwDAOCzyAwAAohWVPwAAMqrwJ/kDACDJqOxP2x8AAMOQ/AEAMAxtfwAAZNYiPyR/AAAko+b50/YHAMAwVP4AAEhG3e1P8gcAQEblftr+AACYhsofAADJqNKf5A8AgD692T/Q3B+USEKPtj8AAIah8gcAQKLtDwCAcQxa5IfkDwCAJOvT/wV6jUjAmD8AAIah8gcAQKLtDwCAcQy64Y+2PwAAhqHyBwBARhX+JH8AACQZlf1p+wMAYBgqfwAAZFThT/IHAECSUdmftj8AAIah8gcAQGKRHwAATGNZlqwA2/aBnj9WSP4AAEhGVf6M+QMAYBhbk//+/fu1aNEipaSkyOFwaPfu3T6/X7lypRwOh8929dVX2xMsACCqDbT9A90iga3Jv7OzU7NmzVJ5efmwx9x0001qbm72blVVVWMYIQDAGFaQtghg65h/Xl6e8vLyznuM0+mUy+Uao4gAAIh+YT/mv2/fPk2bNk2XX3657rjjDrW2ttodEgAgGg0s8hPoFgHC+m7/vLw8ffe731V6eroaGhr085//XPPnz1ddXZ2cTueQ53R1damrq8v7c3t7+1iFCwCIYAYt8BfeyT8/P9/73xkZGcrKylJ6erqef/55LVmyZMhzysrKdP/9949ViAAARJywb/ufKzk5Wenp6Tpx4sSwx5SUlMjj8Xi3pqamMYwQABCxLAWh7W/3m/BPWFf+n9fW1qampiYlJycPe4zT6Rx2SAAAgGGxyM/YOH36tOrr61VfXy9JamhoUH19vRobG3X69GmtXbtWhw4d0rvvvqt9+/Zp0aJFmjp1qm699VY7wwYARCFLQZjnf4HZv6KiQm63W3FxccrMzNSBAwfOe3xXV5fWr1+v9PR0OZ1OffnLX9a2bdv8fj1bK//XX39dN9xwg/fn4uJiSVJBQYEqKyt15MgRbd++XR9//LGSk5N1ww03aOfOnYqPj7crZAAAgmrnzp0qKipSRUWF5s6dqy1btigvL0/Hjh1TWlrakOfcdttt+vDDD/XEE0/oK1/5ilpbW9XT0+P3azqsSFmO6AK1t7crMTHxgs8f+mMHAIxWY4DnezweJSQkBCWWcw3kiTde3qL4iyYGdK2OzjOaveDOUcWanZ2t2bNnq7Ky0rtvxowZWrx4scrKygYd/8ILL2jZsmV65513NHny5AuKM6Ju+AMAIGSCsbTvp/V0e3u7z3buFPRzdXd3q66uTrm5uT77c3NzdfDgwSHP2bNnj7KysrR582Z94Qtf0OWXX661a9fqzJkzfr/ViLrhDwCASJCamurz84YNG1RaWjrouFOnTqm3t1dJSUk++5OSktTS0jLktd955x29+uqriouL065du3Tq1CmtWrVK//znP/0e9yf5AwAgBfVu/6amJp+2/0iz0BwOh+9lLGvQvgF9fX1yOBzasWOHd1j74Ycf1ne+8x098sgjmjhx5KELkj8AAFJQl/hLSEjwa8x/6tSpiomJGVTlt7a2DuoGDEhOTtYXvvAFn/vZZsyYIcuydPLkSV122WUjvi5j/gAA2CQ2NlaZmZmqrq722V9dXa2cnJwhz5k7d64++OADnT592rvv73//u8aNG6fp06f79bokfwAAZN9zfYqLi/X4449r27ZtOn78uO6++241NjaqsLBQUv/KtStWrPAev3z5ck2ZMkXf//73dezYMe3fv1/33HOP/v3f/92vlr9E2x8AgH42PdknPz9fbW1t2rhxo5qbm5WRkaGqqiqlp6dLkpqbm9XY+NlEyYsvvljV1dW66667lJWVpSlTpui2227TAw884PdrMs9/BMzzB4DgCPd5/nV/qtDFAc7zP915Rpl5q0IWa7BQ+QMAIB7pCwCAeQzK/iR/AAAknuoHAACiF5U/AADSZ+vzB3iNSEDyBwBAMmrMn7Y/AACGofIHAEBGFf4kfwAAJBmV/Wn7AwBgGCp/AAAko+b5k/wBAJBZU/1o+wMAYBgqfwAAJNr+AACYxqS2P8kfAIABkZG7A8aYPwAAhqHyBwBAMmqRH5I/AAAya8yftj8AAIah8gcAQGKqHwAApqHtDwAAohaVPwAAktT36RboNSIAyR8AAMmoqX5+t/1PnjwZyjgAAMAY8Tv5Z2Rk6Le//W0oYwEAwDYDhX+gWyTwO/lv2rRJq1ev1tKlS9XW1hbKmAAAGHsGZX+/k/+qVav05ptv6qOPPtIVV1yhPXv2hDIuAADGlEG5f3Q3/Lndbu3du1fl5eVaunSpZsyYofHjfS/xxhtvBDVAAAAQXKO+2/+9997Ts88+q8mTJ+uWW24ZlPwBAIhIBt3tP6rM/dhjj+nHP/6xFixYoLfffluXXnppqOICAGBskfwHu+mmm/SXv/xF5eXlWrFiRShjAgAAIeR38u/t7dVbb72l6dOnhzIeAABsYfX1b4FeIxL4nfyrq6tDGQcAAPayFIS2f1AiCTke7AMAgGG4VR8AABl1vx/JHwAASUZlf9r+AAAYhsofAADJqMqf5A8AgD7N/YFO9YuM3E/yBwBAklGVP2P+AAAYhsofAAAZVfiT/AEAkGRU9qftDwCAYaj8AQCQZFmWrAAr90DPHyskfwAAJKnv0y3Qa0QA2v4AABiGyh8AAMmoG/5I/gAAyKjcT9sfAADTUPkDACBJfVb/Fug1IgDJHwAAMdUPAADzWJ9ugV4jAtg65r9//34tWrRIKSkpcjgc2r17t8/vLctSaWmpUlJSNHHiRF1//fU6evSoPcECABAlbE3+nZ2dmjVrlsrLy4f8/ebNm/Xwww+rvLxctbW1crlcWrhwoTo6OsY4UgBA1Bu43T/QLQLY2vbPy8tTXl7ekL+zLEv/8R//ofXr12vJkiWSpKeeekpJSUl6+umndeedd45lqACAaNdnyTLkhr+wnerX0NCglpYW5ebmevc5nU5dd911Onjw4LDndXV1qb293WcDAACfCdvk39LSIklKSkry2Z+UlOT93VDKysqUmJjo3VJTU0MaJwAgSlhB2iJA2Cb/AQ6Hw+dny7IG7TtXSUmJPB6Pd2tqagp1iACAKGDJ8k73u+AtQrJ/2E71c7lckvo7AMnJyd79ra2tg7oB53I6nXI6nSGPDwCASBW2lb/b7ZbL5VJ1dbV3X3d3t2pqapSTk2NjZACAqNQXpC0C2Jr8T58+rfr6etXX10vqv8mvvr5ejY2NcjgcKioq0qZNm7Rr1y69/fbbWrlypSZNmqTly5fbGTYAIAoF3PIPYIXAiooKud1uxcXFKTMzUwcOHPDrvD//+c8aP368rrzyylG9nq1t/9dff1033HCD9+fi4mJJUkFBgX7zm9/oJz/5ic6cOaNVq1bpo48+UnZ2tl566SXFx8fbFTIAAEG1c+dOFRUVqaKiQnPnztWWLVuUl5enY8eOKS0tbdjzPB6PVqxYoW9961v68MMPR/WaDitSFiK+QO3t7UpMTLzg84f/2AEAo9EY4Pkej0cJCQlBieVcA3mi+oGf6aK4uICu1Xn2rBb+7IFRxZqdna3Zs2ersrLSu2/GjBlavHixysrKhj1v2bJluuyyyxQTE6Pdu3d7u+j+CNsxfwAAxlQQp/p9fr2Zrq6uIV+yu7tbdXV1PmvaSFJubu5517R58skn9Y9//EMbNmy4oLdK8gcAQMEd809NTfVZc2a4Cv7UqVPq7e0d1Zo2J06c0Lp167Rjxw6NH39ho/dhO9UPAIBI1dTU5NP2H2kKur9r2vT29mr58uW6//77dfnll19wfCR/AACk/nX5g7S2f0JCgl9j/lOnTlVMTMygKn+4NW06Ojr0+uuv6/Dhw1qzZk3/S/b1ybIsjR8/Xi+99JLmz58/4uuS/AEAUHAeyjfa82NjY5WZmanq6mrdeuut3v3V1dW65ZZbBh2fkJCgI0eO+OyrqKjQ3r179Yc//EFut9uv1yX5AwBgo+LiYt1+++3KysrSnDlztHXrVjU2NqqwsFBS/7L177//vrZv365x48YpIyPD5/xp06YpLi5u0P7zIfkDACAFte0/Gvn5+Wpra9PGjRvV3NysjIwMVVVVKT09XZLU3NysxsZAJ0r6Yp7/CJjnDwDBEe7z/F/42bqgzPO/6YFfhCzWYGGqHwAAhqHtDwCAFJwH80TIg31I/gAAaOBu/8BGwiNlIJ22PwAAhqHyBwBA6i/bA73bP0JKf5I/AACSz4N5ArpGBCD5AwAg+TyYJ5BrRALG/AEAMAyVPwAAMqvyJ/kDACCjpvnT9gcAwDRU/gAAiLY/AADGMSn50/YHAMAwVP4AAMioNX5I/gAASLT9AQBAFKPyBwBAZlX+JH8AAETyBwDAOCbd8MeYPwAAhqHyBwBAkoLQ9hdtfwAAIodJY/60/QEAMAyVPwAAMuuGP5I/AACi7Q8AAKIYlT8AADKr8if5AwAgs5I/bX8AAAxD5Q8AgKS+T7dArxEJSP4AAMistj/JHwAAmZX8GfMHAMAwVP4AAMisyp/kDwCAzFrel7Y/AACGofIHAEC0/QEAMI6lwJN3ZKR+2v4AABiHyh8AAEl9lqW+ACv/QM8fKyR/AADE3f4AACCKUfkDACDu9gcAwDgkfwAATBOE5K8ISf6M+QMAYBgqfwAAxFQ/AACMY9KYP21/AAAMQ+UPAIDMWuSH5A8AgGj7AwCAKEblDwCAuNsfAADj0PYPE6WlpXI4HD6by+WyOywAACJa2Ff+V1xxhV5++WXvzzExMTZGAwCIViZV/mGf/MePH0+1DwAIOZOSf1i3/SXpxIkTSklJkdvt1rJly/TOO++c9/iuri61t7f7bAAAjKQvSFskCOvkn52dre3bt+vFF1/UY489ppaWFuXk5KitrW3Yc8rKypSYmOjdUlNTxzBiAADCX1gn/7y8PC1dulQzZ87UggUL9Pzzz0uSnnrqqWHPKSkpkcfj8W5NTU1jFS4AIIINtP0D3SJB2I/5n+uiiy7SzJkzdeLEiWGPcTqdcjqdYxgVACAaMOYfprq6unT8+HElJyfbHQoAABErrJP/2rVrVVNTo4aGBv3P//yPvvOd76i9vV0FBQV2hwYAiDIDK/wFul2IiooKud1uxcXFKTMzUwcOHBj22Oeee04LFy7UpZdeqoSEBM2ZM0cvvvjiqF4vrJP/yZMn9W//9m/66le/qiVLlig2Nlavvfaa0tPT7Q4NABBl7Brz37lzp4qKirR+/XodPnxY8+bNU15enhobG4c8fv/+/Vq4cKGqqqpUV1enG264QYsWLdLhw4f9fk2HFSkDFBeovb1diYmJF3x+WhBjAQCTDZ3K/OfxeJSQkBCUWM41kCceXLZME2NjA7rWme5u3ftf/zWqWLOzszV79mxVVlZ6982YMUOLFy9WWVmZX9e44oorlJ+fr/vuu8+v48O68gcAYMwEo+r/tJ7+/HozXV1dQ75kd3e36urqlJub67M/NzdXBw8e9Cvsvr4+dXR0aPLkyX6/VZI/AAAK7ph/amqqz5ozw1Xwp06dUm9vr5KSknz2JyUlqaWlxa+4H3roIXV2duq2227z+71G1FQ/AAAiQVNTk0/bf6Qp6A6Hw+dny7IG7RvKM888o9LSUv3xj3/UtGnT/I6P5A8AgCRLgc/THzg7ISHBrzH/qVOnKiYmZlCV39raOqgb8Hk7d+7UD37wA/3+97/XggULRhUnbX8AAGTP3f6xsbHKzMxUdXW1z/7q6mrl5OQMe94zzzyjlStX6umnn9a3v/3tUb9XKn8AAKSA5umfe43RKi4u1u23366srCzNmTNHW7duVWNjowoLCyX1L1v//vvva/v27ZL6E/+KFSv0n//5n7r66qu9XYOJEyf6PbuN5A8AgI3y8/PV1tamjRs3qrm5WRkZGaqqqvKuadPc3Owz53/Lli3q6enR6tWrtXr1au/+goIC/eY3v/HrNZnnPwLm+QNAcIT7PP+NS5cqbsKEgK519pNPdN+zz4Ys1mCh8gcAQPa1/e3ADX8AABiGyh8AAJn1SF+SPwAAou0PAACiGJU/AACi7Q8AgHFo+wMAgKhF5Q8AgGj7AwBgHCsIbX+SPwAAEcSkyp8xfwAADEPlDwCAzLrbn+QPAIBo+wMAgChG5Q8AgGj7AwBgHJOSP21/AAAMQ+UPAIDMuuGP5A8AgGj7AwCAKEblDwCAzKr8Sf6IWO8F8Ics3eEIYiRAZOHPztAsy5LV1xfwNSIByR8AAJlV+TPmDwCAYaj8AQBQf8s+0Mqdtj8AABGEtj8AAIhaVP4AAEjq6+tTX4B3+wd6/lgh+QMAILOW96XtDwCAYaj8AQCQWTf8kfwBAJBZY/60/QEAMAyVPwAAou0PAIBxSP4AABjGpDF/kj8iVjQ/WhQIJf7sgOQPAID6W/a9tP0BADCHSW1/pvoBAGAYKn8AACT1KfC2fWTU/SR/AAAkfdr2D/BmSNr+AAAgLFH5AwAgqbevT+MCrPx7I6TyJ/kDACCzVvij7Q8AgGGo/AEAEG1/AACMY9Ld/iR/AADUn7h7DUn+jPkDAGAYKn8AANQ/Xu9gzB+S1Gh3AACAMWFS8qftDwCAYaj8AQCQWXf7R0TlX1FRIbfbrbi4OGVmZurAgQN2hwQAiDI9fX3q6e0NbCP5B8fOnTtVVFSk9evX6/Dhw5o3b57y8vLU2MhoPAAAF8JhWeG9EHF2drZmz56tyspK774ZM2Zo8eLFKisrG/H89vZ2JSYmhjJEAMAY8Hg8SkhICPp1B/LEVampGj8usJq4p69PtU1NIYs1WMJ6zL+7u1t1dXVat26dz/7c3FwdPHhwyHO6urrU1dXl/bm9vT2kMQIAokMw7tTnbv8gOHXqlHp7e5WUlOSzPykpSS0tLUOeU1ZWpsTERO+Wmpo6FqECABAxwjr5D/j8vEvLsoadi1lSUiKPx+PdmpqaxiJEAECE6+3tDcoWCcK67T916lTFxMQMqvJbW1sHdQMGOJ1OOZ3OsQgPABBFPunqUl+AY/60/YMgNjZWmZmZqq6u9tlfXV2tnJwcm6ICAESjTz75RJ90dwe2ffKJ3W/DL2Gd/CWpuLhYjz/+uLZt26bjx4/r7rvvVmNjowoLC+0ODQCAoBjtejY1NTXKzMxUXFycvvSlL+nRRx8d1euFddtfkvLz89XW1qaNGzequblZGRkZqqqqUnp6ut2hAQCiSE9XV+Ar/F3A7PmB9WwqKio0d+5cbdmyRXl5eTp27JjS0tIGHd/Q0KCbb75Zd9xxh373u9/pz3/+s1atWqVLL71US5cu9es1w36ef6CY5w8A0SHU8/ynK/B2eJ+kkxpdrKNdz+bee+/Vnj17dPz4ce++wsJCvfnmmzp06JBfrxn2bf9ARfm/bQDAGKH++7wvSJvU/w+Kc7dz158518B6Nrm5uT77z7eezaFDhwYdf+ONN+r111/3+56DsG/7B6qjo8PuEAAAQdDR0RGSTm5sbKxcLpc+GGb9mNG6+OKLB60xs2HDBpWWlg469kLWs2lpaRny+J6eHp06dUrJyckjxhj1yT8lJUVNTU2Kj48ftDZAe3u7UlNT1dTUFNbLMNqNz8k/fE7+4XPyD5/TZyzLUkdHh1JSUkJy/bi4ODU0NKi7uzso1xtqLZqRpqCPZj2b4Y4fav9woj75jxs3TtOnTz/vMQkJCcb/4fIHn5N/+Jz8w+fkHz6nfqG+dysuLk5xcXEhfY2hXMh6Ni6Xa8jjx48frylTpvj1ulE/5g8AQLi6kPVs5syZM+j4l156SVlZWZowYYJfr0vyBwDARiOtZ1NSUqIVK1Z4jy8sLNR7772n4uJiHT9+XNu2bdMTTzyhtWvX+v2aUd/2Px+n06kNGzawHPAI+Jz8w+fkHz4n//A5mWOk9Wyam5vV2NjoPd7tdquqqkp33323HnnkEaWkpOhXv/qV33P8JQPm+QMAAF+0/QEAMAzJHwAAw5D8AQAwDMkfAADDGJ38R/sIRdOUlpbK4XD4bC6Xy+6wbLd//34tWrRIKSkpcjgc2r17t8/vLctSaWmpUlJSNHHiRF1//fU6evSoPcHaaKTPaeXKlYO+X1dffbU9wdqkrKxMV111leLj4zVt2jQtXrxYf/vb33yO4fuEUDA2+Q88QnH9+vU6fPiw5s2bp7y8PJ/pFJCuuOIKNTc3e7cjR47YHZLtOjs7NWvWLJWXlw/5+82bN+vhhx9WeXm5amtr5XK5tHDhQuOeMzHS5yRJN910k8/3q6qqagwjtF9NTY1Wr16t1157TdXV1erp6VFubq46Ozu9x/B9QkhYhvrmN79pFRYW+uz72te+Zq1bt86miMLPhg0brFmzZtkdRliTZO3atcv7c19fn+Vyuaxf/OIX3n1nz561EhMTrUcffdSGCMPD5z8ny7KsgoIC65ZbbrElnnDV2tpqSbJqamosy+L7hNAxsvK/kEcomurEiRNKSUmR2+3WsmXL9M4779gdUlhraGhQS0uLz3fL6XTquuuu47s1hH379mnatGm6/PLLdccdd6i1tdXukGzl8XgkSZMnT5bE9wmhY2Tyv5BHKJooOztb27dv14svvqjHHntMLS0tysnJUVtbm92hha2B7w/frZHl5eVpx44d2rt3rx566CHV1tZq/vz5wz73PNpZlqXi4mJdc801ysjIkMT3CaFj9PK+o32Eomny8vK8/z1z5kzNmTNHX/7yl/XUU0+puLjYxsjCH9+tkeXn53v/OyMjQ1lZWUpPT9fzzz+vJUuW2BiZPdasWaO33npLr7766qDf8X1CsBlZ+V/IIxQhXXTRRZo5c6ZOnDhhdyhha2A2BN+t0UtOTlZ6erqR36+77rpLe/bs0SuvvOLzCHK+TwgVI5P/hTxCEVJXV5eOHz+u5ORku0MJW263Wy6Xy+e71d3drZqaGr5bI2hra1NTU5NR3y/LsrRmzRo999xz2rt3r9xut8/v+T4hVIxt+xcXF+v2229XVlaW5syZo61bt/o8QhHS2rVrtWjRIqWlpam1tVUPPPCA2tvbVVBQYHdotjp9+rT+93//1/tzQ0OD6uvrNXnyZKWlpamoqEibNm3SZZddpssuu0ybNm3SpEmTtHz5chujHnvn+5wmT56s0tJSLV26VMnJyXr33Xf105/+VFOnTtWtt95qY9Rja/Xq1Xr66af1xz/+UfHx8d4KPzExURMnTpTD4eD7hNCwda6BzR555BErPT3dio2NtWbPnu2dXoN++fn5VnJysjVhwgQrJSXFWrJkiXX06FG7w7LdK6+8YkkatBUUFFiW1T89a8OGDZbL5bKcTqd17bXXWkeOHLE3aBuc73P617/+ZeXm5lqXXnqpNWHCBCstLc0qKCiwGhsb7Q57TA31+UiynnzySe8xfJ8QCjzSFwAAwxg55g8AgMlI/gAAGIbkDwCAYUj+AAAYhuQPAIBhSP4AABiG5A8AgGFI/gAAGIbkD0SY3t5e5eTkaOnSpT77PR6PUlNT9bOf/cymyABEClb4AyLQiRMndOWVV2rr1q363ve+J0lasWKF3nzzTdXW1io2NtbmCAGEM5I/EKF+9atfqbS0VG+//bZqa2v13e9+V3/5y1905ZVX2h0agDBH8gcilGVZmj9/vmJiYnTkyBHdddddtPwB+IXkD0Swv/71r5oxY4ZmzpypN954Q+PHG/uUbgCjwA1/QATbtm2bJk2apIaGBp08edLucABECCp/IEIdOnRI1157rf70pz9p8+bN6u3t1csvvyyHw2F3aADCHJU/EIHOnDmjgoIC3XnnnVqwYIEef/xx1dbWasuWLXaHBiACkPyBCLRu3Tr19fXpwQcflCSlpaXpoYce0j333KN3333X3uAAhD3a/kCEqamp0be+9S3t27dP11xzjc/vbrzxRvX09ND+B3BeJH8AAAxD2x8AAMOQ/AEAMAzJHwAAw5D8AQAwDMkfAADDkPwBADAMyR8AAMOQ/AEAMAzJHwAAw5D8AQAwDMkfAADDkPwBADDM/wMYJQqZr06e5gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from landlab import RasterModelGrid\n",
    "import numpy as np\n",
    "\n",
    "mg2 = RasterModelGrid((25, 25))\n",
    "z2 = mg2.add_zeros(\"topographic__elevation\", at=\"node\")\n",
    "z2.reshape((25, 25))[2, 8] = 1.\n",
    "z2.reshape((25, 25))[2, 22] = 1.\n",
    "mg2.set_closed_boundaries_at_grid_edges(True, True, True, True)\n",
    "imshow_grid(mg2,\"topographic__elevation\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "434d395a-26ca-4e83-a7d0-aca181e493fb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.          0.08333333  0.16666667  0.25        0.33333333  0.41666667\n",
      "  0.5         0.58333333  0.66666667  0.75        0.83333333  0.91666667\n",
      "  1.          1.08333333  1.16666667  1.25        1.33333333  1.41666667\n",
      "  1.5         1.58333333  1.66666667  1.75        1.83333333  1.91666667\n",
      "  2.          0.          0.08333333  0.16666667  0.25        0.33333333\n",
      "  0.41666667  0.5         0.58333333  0.66666667  0.75        0.83333333\n",
      "  0.91666667  1.          1.08333333  1.16666667  1.25        1.33333333\n",
      "  1.41666667  1.5         1.58333333  1.66666667  1.75        1.83333333\n",
      "  1.91666667  2.          0.          0.08333333  0.16666667  0.25\n",
      "  0.33333333  0.41666667  0.5         0.58333333  0.66666667  0.75\n",
      "  0.83333333  0.91666667  1.          1.08333333  1.16666667  1.25\n",
      "  1.33333333  1.41666667  1.5         1.58333333  1.66666667  1.75\n",
      "  1.83333333  1.91666667  2.          0.          0.08333333  0.16666667\n",
      "  0.25        0.33333333  0.41666667  0.5         0.58333333  0.66666667\n",
      "  0.75        0.83333333  0.91666667  1.          1.08333333  1.16666667\n",
      "  1.25        1.33333333  1.41666667  1.5         1.58333333  1.66666667\n",
      "  1.75        1.83333333  1.91666667  2.          0.          0.08333333\n",
      "  0.16666667  0.25        0.33333333  0.41666667  0.5         0.58333333\n",
      "  0.66666667  0.75        0.83333333  0.91666667  1.          1.08333333\n",
      "  1.16666667  1.25        1.33333333  1.41666667  1.5         1.58333333\n",
      "  1.66666667  1.75        1.83333333  1.91666667  2.          0.\n",
      "  0.08333333  0.16666667  0.25        0.33333333  0.41666667  0.5\n",
      "  0.58333333  0.66666667  0.75        0.83333333  0.91666667  1.\n",
      "  1.08333333  1.16666667  1.25        1.33333333  1.41666667  1.5\n",
      "  1.58333333  1.66666667  1.75        1.83333333  1.91666667  2.\n",
      "  0.          0.08333333  0.16666667  0.25        0.33333333  0.41666667\n",
      "  0.5         0.58333333  0.66666667  0.75        0.83333333  0.91666667\n",
      "  1.          1.08333333  1.16666667  1.25        1.33333333  1.41666667\n",
      "  1.5         1.58333333  1.66666667  1.75        1.83333333  1.91666667\n",
      "  2.          0.          0.08333333  0.16666667  0.25        0.33333333\n",
      "  0.41666667  0.5         0.58333333  0.66666667  0.75        0.83333333\n",
      "  0.91666667  1.          1.08333333  1.16666667  1.25        1.33333333\n",
      "  1.41666667  1.5         1.58333333  1.66666667  1.75        1.83333333\n",
      "  1.91666667  2.          0.          0.08333333  0.16666667  0.25\n",
      "  0.33333333  0.41666667  0.5         0.58333333  0.66666667  0.75\n",
      "  0.83333333  0.91666667  1.          1.08333333  1.16666667  1.25\n",
      "  1.33333333  1.41666667  1.5         1.58333333  1.66666667  1.75\n",
      "  1.83333333  1.91666667  2.          0.          0.08333333  0.16666667\n",
      "  0.25        0.33333333  0.41666667  0.5         0.58333333  0.66666667\n",
      "  0.75        0.83333333  0.91666667  1.          1.08333333  1.16666667\n",
      "  1.25        1.33333333  1.41666667  1.5         1.58333333  1.66666667\n",
      "  1.75        1.83333333  1.91666667  2.          0.          0.08333333\n",
      "  0.16666667  0.25        0.33333333  0.41666667  0.5         0.58333333\n",
      "  0.66666667  0.75        0.83333333  0.91666667  1.          1.08333333\n",
      "  1.16666667  1.25        1.33333333  1.41666667  1.5         1.58333333\n",
      "  1.66666667  1.75        1.83333333  1.91666667  2.          0.\n",
      "  0.08333333  0.16666667  0.25        0.33333333  0.41666667  0.5\n",
      "  0.58333333  0.66666667  0.75        0.83333333  0.91666667  1.\n",
      "  1.08333333  1.16666667  1.25        1.33333333  1.41666667  1.5\n",
      "  1.58333333  1.66666667  1.75        1.83333333  1.91666667  2.\n",
      "  0.          0.08333333  0.16666667  0.25        0.33333333  0.41666667\n",
      "  0.5         0.58333333  0.66666667  0.75        0.83333333  0.91666667\n",
      "  1.          1.08333333  1.16666667  1.25        1.33333333  1.41666667\n",
      "  1.5         1.58333333  1.66666667  1.75        1.83333333  1.91666667\n",
      "  2.          0.          0.08333333  0.16666667  0.25        0.33333333\n",
      "  0.41666667  0.5         0.58333333  0.66666667  0.75        0.83333333\n",
      "  0.91666667  1.          1.08333333  1.16666667  1.25        1.33333333\n",
      "  1.41666667  1.5         1.58333333  1.66666667  1.75        1.83333333\n",
      "  1.91666667  2.          0.          0.08333333  0.16666667  0.25\n",
      "  0.33333333  0.41666667  0.5         0.58333333  0.66666667  0.75\n",
      "  0.83333333  0.91666667  1.          1.08333333  1.16666667  1.25\n",
      "  1.33333333  1.41666667  1.5         1.58333333  1.66666667  1.75\n",
      "  1.83333333  1.91666667  2.          0.          0.08333333  0.16666667\n",
      "  0.25        0.33333333  0.41666667  0.5         0.58333333  0.66666667\n",
      "  0.75        0.83333333  0.91666667  1.          1.08333333  1.16666667\n",
      "  1.25        1.33333333  1.41666667  1.5         1.58333333  1.66666667\n",
      "  1.75        1.83333333  1.91666667  2.          0.          0.08333333\n",
      "  0.16666667  0.25        0.33333333  0.41666667  0.5         0.58333333\n",
      "  0.66666667  0.75        0.83333333  0.91666667  1.          1.08333333\n",
      "  1.16666667  1.25        1.33333333  1.41666667  1.5         1.58333333\n",
      "  1.66666667  1.75        1.83333333  1.91666667  2.          0.\n",
      "  0.08333333  0.16666667  0.25        0.33333333  0.41666667  0.5\n",
      "  0.58333333  0.66666667  0.75        0.83333333  0.91666667  1.\n",
      "  1.08333333  1.16666667  1.25        1.33333333  1.41666667  1.5\n",
      "  1.58333333  1.66666667  1.75        1.83333333  1.91666667  2.\n",
      "  0.          0.08333333  0.16666667  0.25        0.33333333  0.41666667\n",
      "  0.5         0.58333333  0.66666667  0.75        0.83333333  0.91666667\n",
      "  1.          1.08333333  1.16666667  1.25        1.33333333  1.41666667\n",
      "  1.5         1.58333333  1.66666667  1.75        1.83333333  1.91666667\n",
      "  2.          0.          0.08333333  0.16666667  0.25        0.33333333\n",
      "  0.41666667  0.5         0.58333333  0.66666667  0.75        0.83333333\n",
      "  0.91666667  1.          1.08333333  1.16666667  1.25        1.33333333\n",
      "  1.41666667  1.5         1.58333333  1.66666667  1.75        1.83333333\n",
      "  1.91666667  2.          0.          0.08333333  0.16666667  0.25\n",
      "  0.33333333  0.41666667  0.5         0.58333333  0.66666667  0.75\n",
      "  0.83333333  0.91666667  1.          1.08333333  1.16666667  1.25\n",
      "  1.33333333  1.41666667  1.5         1.58333333  1.66666667  1.75\n",
      "  1.83333333  1.91666667  2.          0.          0.08333333  0.16666667\n",
      "  0.25        0.33333333  0.41666667  0.5         0.58333333  0.66666667\n",
      "  0.75        0.83333333  0.91666667  1.          1.08333333  1.16666667\n",
      "  1.25        1.33333333  1.41666667  1.5         1.58333333  1.66666667\n",
      "  1.75        1.83333333  1.91666667  2.          0.          0.08333333\n",
      "  0.16666667  0.25        0.33333333  0.41666667  0.5         0.58333333\n",
      "  0.66666667  0.75        0.83333333  0.91666667  1.          1.08333333\n",
      "  1.16666667  1.25        1.33333333  1.41666667  1.5         1.58333333\n",
      "  1.66666667  1.75        1.83333333  1.91666667  2.          0.\n",
      "  0.08333333  0.16666667  0.25        0.33333333  0.41666667  0.5\n",
      "  0.58333333  0.66666667  0.75        0.83333333  0.91666667  1.\n",
      "  1.08333333  1.16666667  1.25        1.33333333  1.41666667  1.5\n",
      "  1.58333333  1.66666667  1.75        1.83333333  1.91666667  2.\n",
      "  0.          0.08333333  0.16666667  0.25        0.33333333  0.41666667\n",
      "  0.5         0.58333333  0.66666667  0.75        0.83333333  0.91666667\n",
      "  1.          1.08333333  1.16666667  1.25        1.33333333  1.41666667\n",
      "  1.5         1.58333333  1.66666667  1.75        1.83333333  1.91666667\n",
      "  2.        ]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAggAAAGwCAYAAADMjZ3mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy88F64QAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAv4klEQVR4nO3df3RU9Z3/8deQHzP8SKIhkElKEgdKD2hcLYmFRAPoajB2XanpLtauYqusOfgrpK4KuMeU7iGVIoeD4YcgVNlW5JzFKmebXYlVIi1xBQxqkfXg12AiZr5p+NqEHzqB5H7/QEZn7gCTewfuMDwfnnsOc+fzufPJODDvvN+f+/m4DMMwBAAA8A2DnB4AAACIPwQIAADAhAABAACYECAAAAATAgQAAGBCgAAAAEwIEAAAgEmy0wM42/r7+/XZZ58pLS1NLpfL6eEAAAbIMAwdOnRIubm5GjTo7Pxe++WXX6q3tzcm10pNTZXH44nJtZyU8AHCZ599pry8PKeHAQCwqb29XaNGjYr5db/88kv5fD75/f6YXM/r9aq1tfW8DxISPkBIS0tzeggAgBg4W/+e9/b2yu/3q729Xenp6bau1dPTo7y8PPX29hIgxDvKCgCQGM72v+fp6WlKT7cbhCTO7gUJHyAAABANwzBkd3uiRNreiAABAABJJ377t/sFnzgBArc5AgAAEzIIAABIIoMQigABAAAxByEcJQYAAGBCBgEAAEmUGEIRIAAAIIkAIRQlBgAAYEIGAQAAMUkxHAECAACSKDGEosQAAABMyCCcQb7TAwCABNHm9ADOiAzCNxEgAAAg5iCEI0AAAEASGYRQzEEAAAAmZBAAAJBEBiEUAQIAAGIOQjhKDAAAwIQMAgAAkigxhCJAAABAEgFCKEoMAADAhAwCAABikmI4AgQAACRRYghFiQEAAJiQQQAAQNKJ3/77Y3CNxECAAACAmIMQjgABAABJzEEIxRwEAABgQgYBAABJZBBCESAAACDmIISjxAAAAEzIIAAAIIkSQygCBAAAJBEghKLEAAAATMggAAAgJimGI0AAAEASJYZQlBgAAIAJGQQAAESJIRwBAgAAkk7s5Gh3N0e7/eMHAQIAAJKYgxCKOQgAAMCEDAIAAGIOQjgCBAAAJFFiCEWJAQAAmJBBAABAEhmEUAQIAACIOQjhKDEAAAATAgQAAIIMm4c1K1askM/nk8fjUVFRkbZt23ba9k1NTSoqKpLH49Ho0aO1atWqkOefe+45uVwu0/Hll19GPSYCBAAAJH29kqLdY2A2btyo6upqzZ8/Xy0tLSorK1NFRYXa2toitm9tbdVNN92ksrIytbS0aN68eXrwwQe1adOmkHbp6enq6OgIOTweT9TjYg4CAAAOWrJkie6++27dc889kqSlS5fq1Vdf1cqVK1VXV2dqv2rVKuXn52vp0qWSpPHjx2vnzp1avHixKisrg+1cLpe8Xq/lcZFBAABAX09StHtIUk9PT8gRCAQivmZvb6927dql8vLykPPl5eXavn17xD7Nzc2m9tOmTdPOnTt17Nix4LnDhw+roKBAo0aN0t/93d+ppaVlQO8HAQIAAJLszz/4eh5CXl6eMjIygkekTIAkdXV1qa+vT9nZ2SHns7Oz5ff7I/bx+/0R2x8/flxdXV2SpHHjxum5557T5s2btWHDBnk8Hl199dXat29f1O8GJQYAACTFch2E9vZ2paenB8+63e7T9nK5XKFXMQzTuTO1/+b5SZMmadKkScHnr776ak2YMEFPP/20li1bFsXPQYAAAEDMpaenhwQIp5KVlaWkpCRTtqCzs9OUJTjJ6/VGbJ+cnKzhw4dH7DNo0CBdddVVA8ogUGIAAECxnYMQrdTUVBUVFamxsTHkfGNjo0pLSyP2KSkpMbXfsmWLiouLlZKScsqfbffu3crJyYl6bAQIAABIiuUchIGoqanRs88+q3Xr1mnv3r2aM2eO2traVFVVJUmaO3eu7rzzzmD7qqoqffLJJ6qpqdHevXu1bt06rV27Vg8//HCwzc9//nO9+uqr+vjjj7V7927dfffd2r17d/Ca0aDEAACAg2bMmKGDBw9qwYIF6ujoUGFhoRoaGlRQUCBJ6ujoCFkTwefzqaGhQXPmzNHy5cuVm5urZcuWhdzi+Ne//lX//M//LL/fr4yMDH33u9/Vm2++qe9973tRj8tlJNLC0RH09PQoIyPDcv/8GI4FAC5kkZf9iV53d3dUdf2BOvk90d7+J6WnD7N5rcPKy7v6rI31XCKDAACA2KwpnKNzEOrq6nTVVVcpLS1NI0eO1PTp0/Xhhx+GtDEMQ7W1tcrNzdXgwYM1depU7dmzx6ERAwBwYXA0QGhqatJ9992nt956S42NjTp+/LjKy8t15MiRYJtFixZpyZIlqq+v144dO+T1enXDDTfo0KFDDo4cAJB4nJmkGK/iag7CX/7yF40cOVJNTU2aPHmyDMNQbm6uqqur9eijj0qSAoGAsrOz9eSTT+ree+894zWZgwAA8SHe5yC0tb0ZkzkI+fmTE2IOQlzd5tjd3S1JyszMlHRixyq/3x+y5rTb7daUKVNOuUZ1IBAwrYENAAAGJm4CBMMwVFNTo2uuuUaFhYWSFFwpaiBrVNfV1YWsf52Xl3d2Bw4ASAiG0R+TI1HETYBw//3367333tOGDRtMzw1kjeq5c+equ7s7eLS3t5+V8QIAEg1zEL4pLm5zfOCBB7R582a9+eabGjVqVPD8yX2s/X5/yPKQp1uj2u12n3FTDAAAwnGbYyhHMwiGYej+++/XSy+9pNdff10+ny/keZ/PJ6/XG7LmdG9vr5qamk65RjUAALDP0QzCfffdpxdeeEGvvPKK0tLSgvMKMjIyNHjwYLlcLlVXV2vhwoUaO3asxo4dq4ULF2rIkCG6/fbbnRw6ACDhxG6750TgaICwcuVKSdLUqVNDzv/617/WXXfdJUl65JFH9MUXX2j27Nn6/PPPNXHiRG3ZskVpaWnneLQAgMRmSLI7yTBxAoS4WgfhbGAdBACID/G+DsL+/a8pPX2ozWsd0SWXXJ8Q6yDExSRFAACcxiTFUAQIAABIYg5CqLhZBwEAAMQPMggAAEgigxCKAAEAAElSLJZKZqllAACQwMggAAAgiRJDKAIEAAAkESCEIkAAAECKyXbNbPcMAAASGhkEAAAkUWIIRYAAAIAkAoRQlBgAAIAJGQQAAMQkxXAECAAASKLEEIoSAwAAMCGDAACAJMMwZBj2MgB2+8cTAgQAACSdKA/YnUOQOAECJQYAAGBCBgEAAElMUgxFgAAAgJiDEI4AAQAASTKME4fdayQI5iAAAAATAgQAAGBCiQEAADEHIRwZBAAAYEIGAQAAibscwxAgAAAgSgzhKDEAAAATMggAAEiUGMIQIAAAILFQUhhKDAAAwIQMAgAAIoEQjgABAACJCCEMJQYAAGBCgAAAAEwoMQAAIBZKCkeAAACAxDoIYSgxAAAAEzIIAABI3MUQhgwCAAD6Oj6we1ixYsUK+Xw+eTweFRUVadu2badt39TUpKKiInk8Ho0ePVqrVq06ZdsXX3xRLpdL06dPH9CYCBAAAHDQxo0bVV1drfnz56ulpUVlZWWqqKhQW1tbxPatra266aabVFZWppaWFs2bN08PPvigNm3aZGr7ySef6OGHH1ZZWdmAx0WAAACA5FgKYcmSJbr77rt1zz33aPz48Vq6dKny8vK0cuXKiO1XrVql/Px8LV26VOPHj9c999yjn/70p1q8eHFIu76+Pv34xz/Wz3/+c40ePXrA4yJAAABAX93EYDc++OpaPT09IUcgEIj4mr29vdq1a5fKy8tDzpeXl2v79u0R+zQ3N5vaT5s2TTt37tSxY8eC5xYsWKARI0bo7rvvtvR+ECAAABBjeXl5ysjICB51dXUR23V1damvr0/Z2dkh57Ozs+X3+yP28fv9EdsfP35cXV1dkqQ//elPWrt2rdasWWP5Z+AuBgAApJjexdDe3q709PTgabfbfdpuLpcr7DKG6dyZ2p88f+jQIf3TP/2T1qxZo6ysrAEN/5sIEAAAkGK6UFJ6enpIgHAqWVlZSkpKMmULOjs7TVmCk7xeb8T2ycnJGj58uPbs2aP9+/fr5ptvDj7f398vSUpOTtaHH36oMWPGnHFslBgAAJBkxOi/gUhNTVVRUZEaGxtDzjc2Nqq0tDRin5KSElP7LVu2qLi4WCkpKRo3bpzef/997d69O3j8/d//va699lrt3r1beXl5UY2NDAIAAA6qqanRHXfcoeLiYpWUlGj16tVqa2tTVVWVJGnu3Lk6cOCA1q9fL0mqqqpSfX29ampqNGvWLDU3N2vt2rXasGGDJMnj8aiwsDDkNS666CJJMp0/HQIEAAAkx/ZimDFjhg4ePKgFCxaoo6NDhYWFamhoUEFBgSSpo6MjZE0En8+nhoYGzZkzR8uXL1dubq6WLVumyspKm4MP5TISaeupCHp6epSRkWG5f34MxwIAF7LIy/5Er7u7O6q6/kCd/J7489u/VtqwIbaudejwURV+7ydnbaznEnMQAACACSUGAADEXk3hCBAAAJCIEMJQYgAAACZkEAAAEAmEcAQIAABIRAhhKDEAAAATMggAAEiOLZQUrwgQAADQiR0R7a4dmEhrDxIgAAAgkUEIwxwEAABgQgYBAABRYgjnaAbhzTff1M0336zc3Fy5XC69/PLLIc/fddddcrlcIcekSZOcGSwAILEZMToShKMBwpEjR3TFFVeovr7+lG1uvPFGdXR0BI+GhoZzOEIAAC5MjpYYKioqVFFRcdo2brdbXq/3HI0IAHDBYqGkEHE/SXHr1q0aOXKkvvOd72jWrFnq7Ow8bftAIKCenp6QAwCAMzkZH9g9EkVcBwgVFRX67W9/q9dff11PPfWUduzYoeuuu06BQOCUferq6pSRkRE88vLyzuGIAQBIDHF9F8OMGTOCfy4sLFRxcbEKCgr0+9//XrfeemvEPnPnzlVNTU3wcU9PD0ECAODMDMWgxBCTkcSFuA4QwuXk5KigoED79u07ZRu32y23230ORwUASAgslBTivAoQDh48qPb2duXk5Dg9FABAgjEUg3UQEihCcDRAOHz4sD766KPg49bWVu3evVuZmZnKzMxUbW2tKisrlZOTo/3792vevHnKysrSD37wAwdHDQBA4nM0QNi5c6euvfba4OOTcwdmzpyplStX6v3339f69ev117/+VTk5Obr22mu1ceNGpaWlOTVkAECiosQQwtEAYerUqadN57z66qvncDQAgAtaDJZaTqT7HOP6NkcAAOCM82qSIgAAZw0lhhAECAAASCy1HIYSAwAAMCGDAACASCCEI0AAAEAiQghDiQEAAJiQQQAAQCQQwhEgAAAgESGEIUAAAEBiHYQwzEEAAAAmZBAAAJBkxGAvBtt7OcQRAgQAACTmIIShxAAAAEzIIAAAIBII4QgQAACQiBDCUGIAAAAmZBAAAJBYByEMAQIAAOI2x3CUGAAAgAkZBAAAJEoMYQgQAAAQJYZwBAgAAJyUON/vtjEHAQAAmJBBAABAYqGkMAQIAACIOQjhKDEAAAATMggAAEjc5hiGAAEAAFFiCEeJAQAAh61YsUI+n08ej0dFRUXatm3bads3NTWpqKhIHo9Ho0eP1qpVq0Kef+mll1RcXKyLLrpIQ4cO1ZVXXql///d/H9CYyCAAwHkoJTXV6SEMXG+v0yM4vf6vDrvXGKCNGzequrpaK1as0NVXX61nnnlGFRUV+uCDD5Sfn29q39raqptuukmzZs3Sb37zG/3pT3/S7NmzNWLECFVWVkqSMjMzNX/+fI0bN06pqan6z//8T/3kJz/RyJEjNW3atKjG5TISKR8SQU9PjzIyMiz3N/+vAQDnnY8Bwv+xGSB0d3crPT09RqP52snvie2/XaxhQwbbutbho1+o9McPD2isEydO1IQJE7Ry5crgufHjx2v69Omqq6sztX/00Ue1efNm7d27N3iuqqpK7777rpqbm0/5OhMmTND3v/99/eIXv4hqXJQYAACIsZ6enpAjEAhEbNfb26tdu3apvLw85Hx5ebm2b98esU9zc7Op/bRp07Rz504dO3bM1N4wDP3hD3/Qhx9+qMmTJ0f9M0QdIHz66adRXxQAgPPNyXWS7B6SlJeXp4yMjOARKRMgSV1dXerr61N2dnbI+ezsbPn9/oh9/H5/xPbHjx9XV1dX8Fx3d7eGDRum1NRUff/739fTTz+tG264Ier3I+o5CIWFhXr66ad1xx13RH1xAADOGzFcSbG9vT2kxOB2u0/bzeVyhV3GMJ07U/vw82lpadq9e7cOHz6sP/zhD6qpqdHo0aM1derUqH6UqAOEhQsX6r777tPLL7+s1atXa/jw4dF2BQAg7sVypeX09PSo5iBkZWUpKSnJlC3o7Ow0ZQlO8nq9EdsnJyeHfDcPGjRI3/72tyVJV155pfbu3au6urqoA4SoSwyzZ8/Wu+++q88//1yXXXaZNm/eHG1XAAAQQWpqqoqKitTY2BhyvrGxUaWlpRH7lJSUmNpv2bJFxcXFSklJOeVrGYZxyrkQkQzoNkefz6fXX39d9fX1qqys1Pjx45WcHHqJd955ZyCXBAAgPji0WVNNTY3uuOMOFRcXq6SkRKtXr1ZbW5uqqqokSXPnztWBAwe0fv16SSfuWKivr1dNTY1mzZql5uZmrV27Vhs2bAhes66uTsXFxRozZox6e3vV0NCg9evXh9wpcSYDXgfhk08+0aZNm5SZmalbbrnFFCAAAHBecihAmDFjhg4ePKgFCxaoo6NDhYWFamhoUEFBgSSpo6NDbW1twfY+n08NDQ2aM2eOli9frtzcXC1btiy4BoIkHTlyRLNnz9ann36qwYMHa9y4cfrNb36jGTNmRD2uAa2DsGbNGv3sZz/T9ddfr2eeeUYjRoyI+oWcwjoIABIR6yDEzsnviT8+92RM1kG45q5Hz9pYz6Wof/2/8cYb9fbbb6u+vl533nnn2RwTAADnnNF/4rB7jUQRdYDQ19en9957T6NGjTqb4wEAwBmGYlBiiMlI4kLUAUL4jEkAAJC4mGEIAIAcm6MYtwgQAMAmJyYMppxhZb6z9rp2ftaDB2M3kLOBCCEEmzUBAAATMggAAEhkEMIQIAAAoK/iA7u3OSZOfECAAACAJDIIYZiDAAAATMggAAAgEgjhCBAAAJCIEMJQYgAAACZkEAAAkGQYhgawwfEpr5EoCBAAAJCk/q8Ou9dIEJQYAACACRkEAAAkJimGIUAAAEDEB+EoMQAAABMyCAAge9sY29l62erreoYMsfyayYOs/26YlJRkuW/cb/fcb5w47F4jQRAgAAAgbnMMR4AAAIAkGV8ddq+RIJiDAAAATBwNEN58803dfPPNys3Nlcvl0ssvvxzyvGEYqq2tVW5urgYPHqypU6dqz549zgwWAJDYTt7GYPdIEI4GCEeOHNEVV1yh+vr6iM8vWrRIS5YsUX19vXbs2CGv16sbbrhBhw4dOscjBQAkvH5Dhs2DSYoxUlFRoYqKiojPGYahpUuXav78+br11lslSc8//7yys7P1wgsv6N57743YLxAIKBAIBB/39PTEfuAAACS4uJ2D0NraKr/fr/Ly8uA5t9utKVOmaPv27afsV1dXp4yMjOCRl5d3LoYLADjfGTE6EkTcBgh+v1+SlJ2dHXI+Ozs7+Fwkc+fOVXd3d/Bob28/q+MEACQGQ0bwVkfLRwJFCHF/m6PL5Qp5bBiG6dw3ud1uuW0sWgIAAOI4g+D1eiXJlC3o7Ow0ZRUAALCtP0ZHgojbAMHn88nr9aqxsTF4rre3V01NTSotLXVwZACARGS7vBCDlRjjiaMlhsOHD+ujjz4KPm5tbdXu3buVmZmp/Px8VVdXa+HChRo7dqzGjh2rhQsXasiQIbr99tsdHDUAAInP0QBh586duvbaa4OPa2pqJEkzZ87Uc889p0ceeURffPGFZs+erc8//1wTJ07Uli1blJaW5tSQAQCJis2aQjgaIEydOvW06RiXy6Xa2lrV1taeu0EBAC5M7MUQIu7vYgCAaNnZsnmIjcykna2Xh3o8lvqlJFv/5zvVRt9BNraKjnfs5hgqcf9PAwAAy8ggAAAgMQchDAECAACKzWaMCVRhoMQAAADMyCAAACBRYghDgAAAgLiLIRwlBgAAYEIGAQAAKTabLSXQZk0ECAAA6ORdDHZLDDEaTBygxAAAAEzIIAAAIJ349d/uXQgJlEIgQAAAQGKzpjAECAAAiNscwzEHAQAAmJBBgOOsbtGb4nbHeCRn37FAwFq/3t4YjyS+Wd162c52z8OGDbPe18Z2z2mDB1t7TYvbREv2tntOTkqy3DfekUEIRYAAAIBYBiEcJQYAAGBCgAAAgL4uMdg9rFixYoV8Pp88Ho+Kioq0bdu207ZvampSUVGRPB6PRo8erVWrVoU8v2bNGpWVleniiy/WxRdfrOuvv15vv/32gMZEgAAAgJwLEDZu3Kjq6mrNnz9fLS0tKisrU0VFhdra2iK2b21t1U033aSysjK1tLRo3rx5evDBB7Vp06Zgm61bt+pHP/qR3njjDTU3Nys/P1/l5eU6cOBA1ONyGYk0oyKCnp4eZWRkWO6fH8OxIDImKUbRj0mKUWGSYnScmqS46o03LPeVpO7ubqWnp9u6RiQnvyde+dnPNNTmvytHAgHd8tRTAxrrxIkTNWHCBK1cuTJ4bvz48Zo+fbrq6upM7R999FFt3rxZe/fuDZ6rqqrSu+++q+bm5oiv0dfXp4svvlj19fW68847oxoXGQQAAPT1Okl2D+lE0PHNI3CKXw56e3u1a9culZeXh5wvLy/X9u3bI/Zpbm42tZ82bZp27typY8eORexz9OhRHTt2TJmZmad9D76JAAEAAMW2xJCXl6eMjIzgESkTIEldXV3q6+tTdnZ2yPns7Gz5/f6Iffx+f8T2x48fV1dXV8Q+jz32mL71rW/p+uuvj/r94DZHAABirL29PaTE4D5D6cLlcoU8NgzDdO5M7SOdl6RFixZpw4YN2rp1qzwDKE0RIAAAoNgulJSenh7VHISsrCwlJSWZsgWdnZ2mLMFJXq83Yvvk5GQNHz485PzixYu1cOFCvfbaa/qbv/mbgfwolBgAAJCcuYshNTVVRUVFamxsDDnf2Nio0tLSiH1KSkpM7bds2aLi4mKlpKQEz/3qV7/SL37xC/33f/+3iouLBzQuiQABAABJsZ2kOBA1NTV69tlntW7dOu3du1dz5sxRW1ubqqqqJElz584NufOgqqpKn3zyiWpqarR3716tW7dOa9eu1cMPPxxss2jRIj3++ONat26dLrnkEvn9fvn9fh0+fDjqcVFiAADAQTNmzNDBgwe1YMECdXR0qLCwUA0NDSooKJAkdXR0hKyJ4PP51NDQoDlz5mj58uXKzc3VsmXLVFlZGWyzYsUK9fb26oc//GHIaz3xxBOqra2NalwECAAASFIM5iDIYv/Zs2dr9uzZEZ977rnnTOemTJmid95555TX279/v6VxfBMBAoLsLDJjdWEbyfoCNW4b47UjYGPRosAp7lE+k6OHDll+TauLMznJ6mfRzqJoGUOHWu57kY1FltItLpSUbmNxpiE2/u6wm+OZr5EomIMAAABMyCAAACDrkwzDr5EoCBAAABAlhnCUGAAAgAkZBAAARAYhHAECAAAiQAhHiQEAAJiQQQAAQFL/V4fdayQKAgQAAESJIRwBAgAAIkAIxxwEAABgQgYBAACRQQhHgAAAgFhqORwlBgAAYEIGIcHY2bI5IyvLct+RF19sue/w9HRL/QanpFh+zX4bacBDX3xhue9fjxyx1C/Jxha7Xx49arnvMRtbW9v5LFrdAtzOls3DbWxZnmVjm+nhFn/Wi2z8rEPdbst9U5IT92uDEkOoxP0/DQDAABiy/wWfOOEBJQYAABABGQQAAHSi9Gin/HjyGomCAAEAAHEXQzhKDAAAwIQMAgAA4i6GcAQIAACIACEcAQIAAJIUgwBBCRQgMAcBAACYkEEAAEDc5hiOAAEAADEHIRwlBgAAYEIGAQAAsVBSOAIEAABEiSEcAUKCGWJjy9qc4cMt9/12To7lvvkWt5ke4vFYfs0vbWxj3NXTY7lvx+efW+qXamOL3f83yHolMXDsmOW+draoHjZkiKV+6Rb7Sfa2e/ba2O55pMW+dsabZuPvTvJg69us4/xCgAAAgLiLIRwBAgAAosQQjrsYAACASVwHCLW1tXK5XCGH1+t1elgAgAR0MoNg90gUcV9iuOyyy/Taa68FH9uZ+AQAwKlQYggV9wFCcnLygLIGgUBAgUAg+LjHxoxzAMCFo/+rw+41EkVclxgkad++fcrNzZXP59Ntt92mjz/++LTt6+rqlJGRETzy8vLO0UgBAEgccR0gTJw4UevXr9err76qNWvWyO/3q7S0VAcPHjxln7lz56q7uzt4tLe3n8MRAwDOV8xBCBXXJYaKiorgny+//HKVlJRozJgxev7551VTUxOxj9vtltvtPldDBAAkCOYghIrrDEK4oUOH6vLLL9e+ffucHgoAAAntvAoQAoGA9u7dqxwby/oCABDJyZUU7R6JIq4DhIcfflhNTU1qbW3V//zP/+iHP/yhenp6NHPmTKeHBgBIMMxBCBXXcxA+/fRT/ehHP1JXV5dGjBihSZMm6a233lJBQYHTQwMAIKHFdYDw4osvOj0EAMCFIhYZADIIOJtSUlMt9x1hcetkyd6Wzd/1+Sz3HZNn7XVTMwdbfs3eri8s9237v3+x3Nfqts126pq9x49b7qujRy13ddv4HA+1uB1xhkPbPVvdslmScjMzLfVLz0m3/JqpWdbfp+TBifu1wW6OoeJ6DgIAAHBG4oaCAAAMgCH76xgkTv6AAAEAAEkslBSOAAEAADEHIRxzEAAAgAkZBAAARIkhHBkEAADk7FLLK1askM/nk8fjUVFRkbZt23ba9k1NTSoqKpLH49Ho0aO1atWqkOf37NmjyspKXXLJJXK5XFq6dOmAx0SAAACAgzZu3Kjq6mrNnz9fLS0tKisrU0VFhdra2iK2b21t1U033aSysjK1tLRo3rx5evDBB7Vp06Zgm6NHj2r06NH65S9/Ka/Xa2lcBAgAAMi5vRiWLFmiu+++W/fcc4/Gjx+vpUuXKi8vTytXrozYftWqVcrPz9fSpUs1fvx43XPPPfrpT3+qxYsXB9tcddVV+tWvfqXbbrtNbrfb0vtBgAAAgGJbYujp6Qk5AoFAxNfs7e3Vrl27VF5eHnK+vLxc27dvj9inubnZ1H7atGnauXOnjh07FoN34gQCBAAAYiwvL08ZGRnBo66uLmK7rq4u9fX1KTs7O+R8dna2/H5/xD5+vz9i++PHj6urqys2P4C4iwEAAEmxvYuhvb1d6elf75dxpjS/y+UyXSf83JnaRzpvBwECAACK7UJJ6enpIQHCqWRlZSkpKcmULejs7DRlCU7yer0R2ycnJ2v48OEWR25GiQEAAIekpqaqqKhIjY2NIecbGxtVWloasU9JSYmp/ZYtW1RcXKyUlJSYjY0MQhwaYmPb2RE2tp0dfYpoNRrjLrvEct8RV+db6peaYW1LYEk62nHIct+kHdb/AgYsTiA69OWXll+zx8aWzX39/Zb7Jg2y/vvHYItbRafb2O754mHDLPfNiuI3xVNJz7XWd+glF1l+TU/WUMt9k4fE7gso3ji1UFJNTY3uuOMOFRcXq6SkRKtXr1ZbW5uqqqokSXPnztWBAwe0fv16SVJVVZXq6+tVU1OjWbNmqbm5WWvXrtWGDRuC1+zt7dUHH3wQ/POBAwe0e/duDRs2TN/+9rejGhcBAgAAOvHlbrfEYCVAmDFjhg4ePKgFCxaoo6NDhYWFamhoUEFBgSSpo6MjZE0En8+nhoYGzZkzR8uXL1dubq6WLVumysrKYJvPPvtM3/3ud4OPFy9erMWLF2vKlCnaunVrVOMiQAAAQM4utTx79mzNnj074nPPPfec6dyUKVP0zjvvnPJ6l1xyie2fhTkIAADAhAwCAABiu+dwBAgAAIjdHMNRYgAAACZkEAAAECWGcAQIAACIACEcJQYAAGBCBgEAADFJMRwBAgAAosQQjhIDAAAwIYMAAIDIIIQjQIhDHhs70qV5rO9wONLGTpDpl46w3Hf4mEst9XO7R1p+zUPpeyz37f3rF5b7jvjsc0v90gcPtvyaw2x8JqzuPilJg2zs5ui2uGWtx8ZWt0Pdbst97fz/8WRb21lx8Ejru096RtjYzXFwgu/maGMH05PXSBQECAAAiAxCOOYgAAAAEzIIAADoRHnAbgaAEgMAAAmGEkMoSgwAAMCEDAIAAJL6+/vVb/MuBrv94wkBAgAAYqnlcJQYAACACRkEAADEJMVwBAgAAIg5COEoMQAAABMyCAAAiBJDOAIEAABEgBCOAAEAADEHIRwBQoKxs8VuarL1j0NqmvWtcj2eHEv93O5sy6957Fi35b6pF31mue/g1FRrr2nj/01yUpLlvnY+T0k2+iZb7GvrM2xjq+ikYdb7Jg+19plIHmatnySlDLX+9zU5+WLLfXF+IUAAAEAnygN9lBiCCBAAABAlhnDc5ggAAEzIIAAAIKlf9ksEiZM/IEAAAEDSVyUGl8v2NRIFJQYAAGBCBgEAAEl9/f0aZDOD0JdAGQQCBAAAxEqK4SgxAAAAEzIIAACIEkM4AgQAAMRdDOEIEAAA0Ikv9z4ChCDmIAAAABMyCAAA6MT8ARdzEIIIEM6gzYnXbG+33HeHjb6rt2613FdLrXcFgHhAgBCKEgMAADAhgwAAgLiLIRwBAgAAko7H4Ms9FteIF+dFiWHFihXy+XzyeDwqKirStm3bnB4SAAAJLe4DhI0bN6q6ulrz589XS0uLysrKVFFRobY2J6YPAgASVV9fX0yOROEyjPjeWWLixImaMGGCVq5cGTw3fvx4TZ8+XXV1dab2gUBAgUAg+Linp0d5eXnnZKwAgLOnu7tb6enpMb9uT0+PMjIyNOFb31LSIHu/N/f19+udAwfO2ljPpbjOIPT29mrXrl0qLy8POV9eXq7t27dH7FNXV6eMjIzgQXAAAMDAxXWA0NXVpb6+PmVnZ4ecz87Olt/vj9hn7ty56u7uDh7tNtYFAABcOCgxhDov7mIIX7jCMIxTLmbhdrvldrvPxbAAAAnkWCCg/hiUGBJFXAcIWVlZSkpKMmULOjs7TVkFAADsOHbsmO11EPrie1rfgMR1iSE1NVVFRUVqbGwMOd/Y2KjS0lKHRgUAQGwN9Hb+pqYmFRUVyePxaPTo0Vq1apWpzaZNm3TppZfK7Xbr0ksv1e9+97uBDcqIcy+++KKRkpJirF271vjggw+M6upqY+jQocb+/fuj6t/d3W1I4uDg4OA4z4/u7u6z8j1z8nvCl5JijElNtXX4UlIGPNaT33Nr1qwxPvjgA+Ohhx4yhg4danzyyScR23/88cfGkCFDjIceesj44IMPjDVr1hgpKSnGf/zHfwTbbN++3UhKSjIWLlxo7N2711i4cKGRnJxsvPXWW1GPK+4DBMMwjOXLlxsFBQVGamqqMWHCBKOpqSnqvgQIHBwcHIlxnO0AYZRk5Ns8RlkY6/e+9z2jqqoq5Ny4ceOMxx57LGL7Rx55xBg3blzIuXvvvdeYNGlS8PE//uM/GjfeeGNIm2nTphm33XZb1OOK6xLDSbNnz9b+/fsVCAS0a9cuTZ48Oeq+RgLVgwDgQna2/z3vj9EhnVhb4ZvHN9fn+SYrt/M3Nzeb2k+bNk07d+7UsWPHTtvmVNeMJK4nKcbCoUOHnB4CACAGDh06pIyMjJhfNzU1VV6vV5+d4vb5gRo2bJhpDZ4nnnhCtbW1prZWbuf3+/0R2x8/flxdXV3Kyck5ZZtTXTOShA8QcnNz1d7errS0NNOtkSdXWWxvbz/vV7w6m3ifosP7FB3ep+jwPn3NMAwdOnRIubm5Z+X6Ho9Hra2t6u3tjcn1jAi34p/p9vuB3M5/qvbh5wd6zXAJHyAMGjRIo0aNOm2b9PT0C/4vYDR4n6LD+xQd3qfo8D6dcDYyB9/k8Xjk8XjO6mtEYuV2fq/XG7F9cnKyhg8ffto2A1ki4LyYgwAAQCKycjt/SUmJqf2WLVtUXFyslJSU07YZyBIBCZ9BAAAgntXU1OiOO+5QcXGxSkpKtHr1arW1tamqqkrSiS0EDhw4oPXr10uSqqqqVF9fr5qaGs2aNUvNzc1au3atNmzYELzmQw89pMmTJ+vJJ5/ULbfcoldeeUWvvfaa/vjHP0Y9rgs6QHC73XriiSdYmvkMeJ+iw/sUHd6n6PA+XThmzJihgwcPasGCBero6FBhYaEaGhpUUFAgSero6FBbW1uwvc/nU0NDg+bMmaPly5crNzdXy5YtU2VlZbBNaWmpXnzxRT3++OP613/9V40ZM0YbN27UxIkTox5X3G/3DAAAzj3mIAAAABMCBAAAYEKAAAAATAgQAACAyQUdIAx0e80LTW1trVwuV8jh9XqdHpbj3nzzTd18883Kzc2Vy+XSyy+/HPK8YRiqra1Vbm6uBg8erKlTp2rPnj3ODNZBZ3qf7rrrLtPna9KkSc4M1iF1dXW66qqrlJaWppEjR2r69On68MMPQ9rweYJTLtgAYePGjaqurtb8+fPV0tKisrIyVVRUhNxKAumyyy5TR0dH8Hj//fedHpLjjhw5oiuuuEL19fURn1+0aJGWLFmi+vp67dixQ16vVzfccMMFty/Imd4nSbrxxhtDPl8NDQ3ncITOa2pq0n333ae33npLjY2NOn78uMrLy3XkyJFgGz5PcEzU+z4mmIFur3kheuKJJ4wrrrjC6WHENUnG7373u+Dj/v5+w+v1Gr/85S+D57788ksjIyPDWLVqlQMjjA/h75NhGMbMmTONW265xZHxxKvOzk5DUnBLez5PcNIFmUGwsr3mhWrfvn3Kzc2Vz+fTbbfdpo8//tjpIcW11tZW+f3+kM+W2+3WlClT+GxFsHXrVo0cOVLf+c53NGvWLHV2djo9JEd1d3dLkjIzMyXxeYKzLsgAwcr2mheiiRMnav369Xr11Ve1Zs0a+f1+lZaW6uDBg04PLW6d/Pzw2TqziooK/fa3v9Xrr7+up556Sjt27NB1112nQCDg9NAcYRiGampqdM0116iwsFASnyc464JeatnuVpiJrqKiIvjnyy+/XCUlJRozZoyef/551dTUODiy+Mdn68xmzJgR/HNhYaGKi4tVUFCg3//+97r11lsdHJkz7r//fr333nsR18rn8wQnXJAZBCvba0IaOnSoLr/8cu3bt8/pocStk3d58NkauJycHBUUFFyQn68HHnhAmzdv1htvvBGyPT2fJzjpggwQrGyvCSkQCGjv3r3Kyclxeihxy+fzyev1hny2ent71dTUxGfrDA4ePKj29vYL6vNlGIbuv/9+vfTSS3r99dfl8/lCnufzBCddsCWGM22vCenhhx/WzTffrPz8fHV2durf/u3f1NPTo5kzZzo9NEcdPnxYH330UfBxa2urdu/erczMTOXn56u6uloLFy7U2LFjNXbsWC1cuFBDhgzR7bff7uCoz73TvU+ZmZmqra1VZWWlcnJytH//fs2bN09ZWVn6wQ9+4OCoz6377rtPL7zwgl555RWlpaUFMwUZGRkaPHiwXC4Xnyc4x9F7KBy2fPlyo6CgwEhNTTUmTJgQvLUIJ8yYMcPIyckxUlJSjNzcXOPWW2819uzZ4/SwHPfGG28YkkzHzJkzDcM4cWvaE088YXi9XsPtdhuTJ0823n//fWcH7YDTvU9Hjx41ysvLjREjRhgpKSlGfn6+MXPmTKOtrc3pYZ9Tkd4fScavf/3rYBs+T3AK2z0DAACTC3IOAgAAOD0CBAAAYEKAAAAATAgQAACACQECAAAwIUAAAAAmBAgAAMCEAAEAAJgQIAAAABMCBOA809fXp9LSUlVWVoac7+7uVl5enh5//HGHRgYgkbDUMnAe2rdvn6688kqtXr1aP/7xjyVJd955p959913t2LFDqampDo8QwPmOAAE4Ty1btky1tbX685//rB07dugf/uEf9Pbbb+vKK690emgAEgABAnCeMgxD1113nZKSkvT+++/rgQceoLwAIGYIEIDz2P/+7/9q/Pjxuvzyy/XOO+8oOTnZ6SEBSBBMUgTOY+vWrdOQIUPU2tqqTz/91OnhAEggZBCA81Rzc7MmT56s//qv/9KiRYvU19en1157TS6Xy+mhAUgAZBCA89AXX3yhmTNn6t5779X111+vZ599Vjt27NAzzzzj9NAAJAgCBOA89Nhjj6m/v19PPvmkJCk/P19PPfWU/uVf/kX79+93dnAAEgIlBuA809TUpL/927/V1q1bdc0114Q8N23aNB0/fpxSAwDbCBAAAIAJJQYAAGBCgAAAAEwIEAAAgAkBAgAAMCFAAAAAJgQIAADAhAABAACYECAAAAATAgQAAGBCgAAAAEwIEAAAgMn/Bw/Wh/I4dARoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#print(mg2.node_x)\n",
    "#print(mg2.node_x.mean())\n",
    "kd = mg2.node_x/mg2.node_x.mean()\n",
    "\n",
    "print(kd)\n",
    "\n",
    "ld2 = LinearDiffuser(mg2, linear_diffusivity=kd)\n",
    "for i in range(10):\n",
    "    ld2.run_one_step(0.1)\n",
    "z2[mg2.core_nodes].sum() == 2.\n",
    "\n",
    "z2.reshape((25, 25))[2, 8] > z2.reshape((25, 25))[2, 22]\n",
    "imshow_grid(mg2, \"topographic__elevation\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c5eecbd-9f5e-4a36-9012-8c9e536405a0",
   "metadata": {},
   "source": [
    "# Include Precipitation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f847cf9-7ddd-4a64-85c3-22c7f81a1946",
   "metadata": {},
   "outputs": [],
   "source": [
    "# precipitation increase w/ elevation by factor of 3\n",
    "# P0 = 80 (reference precipitation rate) \n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "CSDMS",
   "language": "python",
   "name": "csdms-2023"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}