frost-cracking_2D.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://csdms.colorado.edu/wiki/ESPIn2020\"><img style=\"float: center; width: 75%\" src=\"../espin/media/ESPIn.png\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Frost-Cracking 2D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will set up a model to estimate frost-cracking intensities. To this end, we assume that heat is only transported through conduction and bedrock temperatures can therefore be assessed solving the diffusion equation with (transient) bedrock surface temperature as the boundary condition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initializations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we import the necessary modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from landlab import RasterModelGrid\n",
    "from landlab import imshow_grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to set up the computational grid. We will use a raster. `Lx` and `Ly` are the depth and height of the bedrock domain, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "rows = 10\n",
    "cols = 30\n",
    "dx = 100 # in m\n",
    "rg = RasterModelGrid((rows, cols), xy_spacing=dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assign fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dem = rg.add_zeros('topographic__elevation', at='node')\n",
    "Ts_ini = rg.add_zeros('temperature_surface_ini', at='node')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will simply assume a vertical wall and assign a slope along the columns of `dem`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAADxCAYAAAA3MOvfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de7hkVX3m8e/bzf2mXJSn6SaCpr2AMSiEaIiKEqQ1xsY4mmYmSpAJxICCcYwgmWjGIZJEMdEJzLSKoqMgKgwkoohENBdu3cgIDRJaWrChQ9uQQURt6D7v/LH3oYvDOVVrn95Vp6p4P89Tz6m9a69dq+rA+fW6/NaSbSIiIkrMm+sKRETE6EjQiIiIYgkaERFRLEEjIiKKJWhERESxBI2IiCiWoBERMQcknSdpvaRbpnntv0iypL06zp0uabWk2yUd1XH+YEk31699VJL6We8EjYiIufFpYMnUk5L2BY4E7u44dwCwDDiwLnOOpPn1y+cCJwCL68cT7tmmBI2IiDlg+9vAA9O89BHgj4HOzOulwIW2N9peA6wGDpW0ANjN9jWuMrU/Axzdz3pv08+bR0SMkyVLlnjDhg1F165cufIK243+1S/pdcA9tv/vlF6mhcC1Hcdr63OP1s+nnu+bBI2IiEIbNmxgxYoVRddKeq6kzouX217e5fqdgDOAV0338jTn3OV83yRoREQ0Uvw3eYPtQxrc+FnA/sBkK2MRcKOkQ6laEPt2XLsIuLc+v2ia832TMY2IiAbsiaJH8/v6ZttPt72f7f2oAsKLbP8bcBmwTNL2kvanGvC+3vY64CFJL65nTb0FuLS1DzuNBI2IiGJu8OhO0gXANcBzJK2VdPyM72qvAi4CbgW+Bpxke3P98tuAT1ANjn8f+OpsPlkpZWn0iIgyhxxysK+//treFwLz52+3smH31EjImEZERCNP7n9oJ2hERDSSoBEREYWe7F36CRoREcUMNJ8ZNU4SNCIiCtlpaSRoREQ0kqARERHFEjQiIqKI0z011xWIiBgtGQiPiIhCaWlEREShsnWlxlmCRkREI+meioiIQumeioiIBhI0IiKiiGe1wdI4SdCIiGgkLY2IiCiWoBEREQWyYGGCRkREA8nTSNCIiGjA3jzXVZhTCRoREY2kpREREUWyyu28ua5ARMRoceGjO0nnSVov6ZaOc38l6XuSvivpEklP7XjtdEmrJd0u6aiO8wdLurl+7aOS1NYnnU6CRkREI+0EDeDTwJIp564Enm/7BcC/AqcDSDoAWAYcWJc5R9L8usy5wAnA4vox9Z6tStCIiGjAdtGj4D7fBh6Ycu7rtjfVh9cCi+rnS4ELbW+0vQZYDRwqaQGwm+1rXL3pZ4CjW/qo08qYRkREMQMDmz31VuAL9fOFVEFk0tr63KP186nn+yZBIyKigQYD4XtJWtFxvNz28pKCks4ANgGfmzw1XVW6nO+bBI2IiEaK/yZvsH1I07tLOhZ4LXCEt0SotcC+HZctAu6tzy+a5nzfZEwjIqKR1gbCn0DSEuA9wOts/7TjpcuAZZK2l7Q/1YD39bbXAQ9JenE9a+otwKWzevNCaWlERBQqHeQuIekC4HCqbqy1wPuoZkttD1xZz5y91vYf2F4l6SLgVqpuq5O8JTX9bVQzsXYEvlo/+kZP9kSViIhSL3zhgf7Wt77Q+0LgKU/5pZWz6Z4admlpREQ0kLWnIiKiUFa5TdCIiGgkQSMiIgplj/CIiGggLY2IiCjitDTmugIREaNldIOGpB2oss1fCuwD/Ay4BfiK7VUl90jQiIgoZDdae2qoSHo/8FvA1cB1wHpgB+DZwFl1QHmX7e92u0+CRkREsZGecnuD7ffP8NrZkp4O/EKvmyRoREQ0MppBw/ZXery+nqr10VWCRkREA6M+EC7pEOAM4BlUMUCA690Ce0rQiIgoZkZ5ILz2OeDdwM3M4sMkaERENDDqLQ3gR7Yvm23hBI2IiEZGc0yjw/skfQK4Ctg4edL2xSWFEzQiIhoZ+aBxHPBcYFu2dE8ZSNCIiGjXWGSE/7LtX5pt4Wz3GhHRSP+2ex2QayUdMNvCaWlERBTzOGzC9OvAsZLWUI1pZMptRET/DHUrosSSrSmcoBERUWjE157axfZPbN/V65pu98mYRkREE1Xk6P0YPpdK+rCkl0naefKkpGdKOl7SFRS0QtLSiIh4ErB9hKTXACcCh0naHdgE3A58BTjW9r/1uk9aGhERDXjCRY9eJJ0nab2kWzrO7SHpSkl31D9373jtdEmrJd0u6aiO8wdLurl+7aOSNGPd7ctt/yfb+9l+iu09bf+a7TNLAgYkaERENGDsskeBT/PE7qDTgKtsL6bK2D4NoJ4iuww4sC5zjqT5dZlzgROAxfVjqwa6e0nQiIgoVZqiURAzbH8beGDK6aXA+fXz84GjO85faHuj7TXAauBQSQuA3Wxf4ypSfaajTF9kTCMioon+DnLvbXtd9TZeV2+MBLAQuLbjurX1uUfr51PP902CRkREAw1ixl6SVnQcL7e9fJZvO904hbuc736zqmtrbzpigO27Syoy0kFjj9139aJ9ntbjqhnHhBpdUmTm8aeSl1t7n3ZvU3JRO9d0Gb9r+F4lhqk+bdWlnfdq7/dZcJvW6tPbypUrN9ju9Qejt/KoscH2IQ3vfp+kBXUrYwFbdtJbC+zbcd0i4N76/KJpzs9I0tuB9wH38fgFC8c/I3zRPk/jsgv+vOs1ml/wH/f83kM7mrf19ympy7ySurRV33nze1+j3v+JlFwzb952BfcZXH3Krimpz7Yt3KPk+yu5T1vXdP9M7b7XIK/RjEltTfQ5ue8y4FjgrPrnpR3nPy/pbGAfqgHv621vlvSQpBcD1wFvAT7W4z1OAZ5j+/7ZVHCkg0ZExEC1uBahpAuAw6m6sdZS/ev/LOAiSccDdwNvBLC9StJFwK1UuRUnecsiWG+jmom1I/DV+tHND4EHZ1vvBI2IiCZaamnYPmaGl46Y4fozgTOnOb8CeH6Dt74TuFrSV3j8JkxnlxRO0IiIKGSGdYWQRu6uH9vVj0YSNCIiig3tulLFbP8ZgKRdq8PuCxROleS+iIgGRne9woqk50v6DnALsErSSkkHlpZPSyMiopSBgnWlhtxy4I9sfxNA0uHAx4FfKymcoBER0YBHfxOmnScDBoDtqzuXSu8lQSMioomRjxncKem/Ap+tj38XWFNaeKSDxsSmCX5+/0+7XjNvm5aS5VpIEmwtua/oM5Uk9/W+z7xBJkeWfK6COpclyw0mSbC9RLlBJiwO8vtrK3FvgH/KhnnAosxbgT8DLqZKt/82cFxp4ZEOGhERgzbqMcP2vwPvmG35BI2IiFKmaIOlYSTpr22fKunvmKaTzfbrSu6ToBERUWzI59N2NzmG8aGtuUmCRkREEyMaM2yvrJ8eZPtvOl+TdArwrZL79C25T9JzJN3U8fixpFMlvV/SPR3nX9NRZto9cCMihkG1jEhr273OlWOnOfd7pYX71tKwfTtwEDy24cc9wCVUo/Qfsf24JtKUPXD3Ab4h6dkdKzlGRMytFle5HTRJxwD/Edhf0mUdL+0KFC+TPqjuqSOA79u+q8tmMo/tgQuskbQaOBS4ZkB1jIjoachbEd38C7AO2Av4cMf5h4Dvlt5kUEFjGXBBx/HJkt4CrADeVU8Bm2kP3MeRdAJwAsA+e+/ZtwpHRExrRGdP2b4LuAt4ydbcp+9BQ9J2wOuA0+tT5wIfoGrkfYAq4r2Vwr1u6z12lwMc+Kxf8MYeyX0liWfzth1MAmBRXdpKNCxIlGvtvQZ4n4HWeWC7NQ7bjoVt7dY4yMTG3vdpy+g2NCr1Ln8fA55HtTT6fOBh27uVlB/EKrevBm60fR+A7ftsb7Y9QbVI1qH1dTPtgRsRMRwmN9QY5WVu4X8AxwB3UO3295/pvUXsYwYRNI6ho2uq3ix90uuplueFag/cZZK2l7Q/9R64A6hfREQ5Fz6GmO3VwPz6H/CfAl5RWrav3VOSdgKOBE7sOP2Xkg6i+lp/MPlajz1wIyKGwNBPpy3x03rY4CZJf0k1OD4cq9za/imw55Rzb+5y/bR74EZEDItRXUakw5upeplOBt5JNSzwhtLCyQiPiCg1Al1PBV4EXG77x1Sr3TaS7V4jIpoY/YHw1wH/Kumzkn5TDdeVT9CIiCjU5uQpSe+UtErSLZIukLSDpD0kXSnpjvrn7h3Xt7LMku3jgF8EvkiVIf59SZ8oLZ+gERHRRAtRQ9JCqj0tDrH9fKpciWXAacBVthcDV9XHU5dZWgKco61ITrH9KPBV4EJgJdWKHEVGekxjYtMEG9c/3PWaoiSukh3jSpLleiQJtlWXot30WtqxsGznw7YSANt5r4HuftjjvYp+VyW7Ebb1O2/pvYq+45Z2WSzb/XBwf8pa7HnaBthR0qPATlR5aacDh9evnw9cDbyHFpdZkrSEKgC9or7/J4A3Nal0RESUsFtZRsT2PZI+BNwN/Az4uu2vS9rb9rr6mnWSnl4XKVpmqdDvUbUwTqyDUCPpnoqIaKDB0uh7SVrR8Thh8h71WMVSYH+qVb13lvS7Xd62aJmlwvovA74DvLSuy46Sdi0tn5ZGREQT5f1TG2wfMsNrvwGssf0jAEkXA78G3CdpQd3KWACsr69vbZklSb9PtejrHsCz6nv9T6rVyHtKSyMiooGWZk/dDbxY0k6q9os4AriNajmlyU2SjgUurZ+3uczSScBhwI+rz+M7gKd3LdEhLY2IiFKTc2639jb2dZK+BNxItWzSd6hW794FuEjS8VSB5Y319W0us7TR9iOTexvVeRrFHypBIyKiiZZmT9l+H/C+Kac3MkM3UYvLLH1L0nupZm4dCfwh8HelhdM9FRHRgCdc9BhipwE/Am6mWjD2cuBPSgunpRERUchjsMptx15GH59N+X4vjf4Dqv1nNwObbB8iaQ/gC8B+VEujv6ne7hVJpwPH19e/w/YV3e4/8egEP1v3k651KEpCKti5ryjJrcc187Yr2KWspURDbTO4JLj2dj5s53O1tgNgCwmSRcl9g9z5cJDJmi3sfFhan5LfZytGeMFCSTfTpfa2X1Byn0G0NF5he0PH8WSa/FmSTquP3zMlTX4f4BuSnp09NSJiqIxuS+O1bdxkLsY0llKlx1P/PLrj/IW2N9peA0ymyUdEDI0GyX1DxfZdk4/61OL6+XrggdL79DtoGPi6pJUd2ZCPS5Nny/zghcAPO8puTZp8RERfjPpAeJ3c9yXgf9WnFgH/p7R8v7unDrN9b71+ypWSvtfl2qI0+Tr4nACw9267tVPLiIgSBibmuhJb7SSqXpzroEru61jjqqcZWxqSLpe039bUzPa99c/1wCV1Re+r0+OZTZq87eW2D7F9yFN2Lt7WNiKiHaO/CdNG249MHjRN7uvWPfVpqq6lM1SyNvEUknaeXARL0s7Aq4BbGEyafEREX4x+zHhCct8XaZDcN2P3lO2LJH0F+FNghaTP0tEws312j3vvDVxSp6pvA3ze9tck3UD/0+QjIvpg+CNCgdOoUhs6k/uKd+7rNabxKPAwsD2wKw1682zfCfzyNOfvp/9p8hER7RuDmNEruU/Sl22/YabyMwaNenens6m6jV5k+6dbWdfWbdq8mR899FDXa7aZ3zuhbpuCHcbmt3BNyT1K6jJv296fqSiJqyQpryS5b0A7H5bep7UdEltIJCxKHB3groZtJVkOdkfHIUruA5gY/ZHwHp7Z7cVuLY0zgDfaXtVufSIiRpMBj33M6D4o3m1M46Xt1yUiYoS1tDT6KMuChRERDTzJY0aCRkREI+MfNZ7R7cUEjYiIYmMwfaq3u7u9mKAREVHK4M1jHzRmNxAeERFPNP4Nje4SNCIimhj/qNE1MWakg8amzZtZ/+CDXa8pSZbbdpveX0NRYl6PRMLWkggLEha3LbimjYTF0mtKPntRwl3B7odtJRsW7RLYa7fGoiTLks/UUqJca8mRBQl32xTsVNnWjoQD3LlvGPfKaKJeC/BndWY4kuYBO3QkcL+nW/m52IQpImJ0TRQ+htdVwE4dxzsB35g8sP31boUTNCIiClUZ4e1twiTpqZK+JOl7km6T9BJJe0i6UtId9c/dO64/XdJqSbdLOmqWH2MH2z957DNVz3fqcv3jJGhERJSyYaLwUeZvgK/Zfi7VAq+3Ua1Ce5XtxVStgtMAJB0ALAMOBJYA50jq3Qf4RA9LetHkgaSDgZ+VFu5b0JC0r6Rv1tFzlaRT6vPvl3SPpJvqx2s6yrQRRSMi+qatPcIl7Qa8DPhkfd9HbP8/YClwfn3Z+cDR9fOlwIW2N9peA6ym2tiuqVOBL0r6R0n/CHwBOLm0cD8HwjcB77J9Y70Z00pJV9avfcT2hzovnhJF9wG+IenZ2VMjIoZKe+PgzwR+BHxK0i8DK4FTgL1trwOwva5jK9aFwLUd5dfW5xqxfYOk5wLPoZop9T3bj5aW71vQqD/05Ad/SNJtdP+Aj0VRYI2kySh6Tb/qGBHRWPnsqb0kreg4Xm57ecfxNsCLgLfbvk7S31B3Rc1guqlmxZWR9Erb/yDpt6e8tFgSti8uuc9AptzWe42/kGoj88OAkyW9BVhB1Rr5d1qKohERfWOKB7mBDbYP6fL6WmCt7evq4y9RBY37JC2oWxkLgPUd1+/bUX4RcG955Xk58A/Ab03zmoHhCBqSdgG+DJxq+8eSzgU+QFXJDwAfBt5KYRSVdAJwAsBuO+7IvQ880PX9SzZhKslpKLnPdj3yPYo2hGpp06iS+5TkV5TksLSVf1J0n5Y+V1vfT686F+VFlGzUVHKf7draPKmtPJeWNroqyfcoeK+2NAga3e9j/5ukH0p6ju3bqXY0vbV+HAucVf+8tC5yGfB5SWdTdeEvBq5v8H7vq38etzX17mvQkLQtVcD43GTTx/Z9Ha9/HPj7+rAoitbNu+UAC5761NHOsomI0WLazsF4O/A5SdsBdwLHUU1QukjS8VSLB74RwPYqSRdRBZVNwEmzGfOVtD3wBmA/OmKA7f9WUr5vQUOSqGYF3Gb77I7zCyYHeYDXA7fUz7cqikZE9F/ZzKjiu9k3AdN1YR0xw/VnAmdu5dteCjxINfC+sWnhfrY0DgPeDNws6ab63HuBYyQdRBWzfwCcCO1F0YiIvmqpe2oOLbK9ZLaF+zl76p+Yfpzi8i5l2oiiERH9M/Ixg3+R9Eu2b55N4ZFesDAiYpA8wgsWSrqZKuRtAxwn6U6q7ikBtv2CkvskaERENDDCmzC9to2bJGhERJQyIzumYfuuyef12lO/TvWJ/tn2jaX3yYKFERHFytadGuYuLEl/SrWm1Z7AXlTLmPxJafmRbmk8smkTd953X9drtt9225736ZWUB2VJbtv1SBjbrqAuJYmGbSUjtnafkgTAASYttpUk2MZGVq1tvNXWZyq4T8mGRq0lJLaV/FiwkVVrhnuvjBLHAC+0/XMASWcBNwL/vaTwSAeNiIhBG+ZWRKEfADsAP6+Ptwe+X1o4QSMiopSB0R0In7QRWFWvOm7gSOCfJH0UwPY7uhVO0IiIaGAMWhqX1I9JVzcpnKAREVHIjH7QsH1+76tmlqAREdHAqI+DS1oMfBA4gGpsAwDbzywpnym3ERGlCqfbDnlr5FPAuVRr/L0C+Azw2dLCCRoREQ2MQdDY0fZVgGzfZfv9wCtLC6d7KiKiUJUQPtQBocTPJc0D7pB0MnAP8PQeZR4z0kHjJw8/zD9fd13Xa56y554977PLLrv0vGbHHXbofc1223V9vSTRcMftt+95zfYliYYtJSwWXVOQVFZS50EmJLZ1Ta+EzpKEu0HuxDjI3SNbS1psadfHtgx5K6LEqcBOwDuodk99BdUOgUWGrntK0hJJt0taLanbJusREQM3yt1TkuYDb7L9E9trbR9n+w22ry29x1AFjfoD/S3waqqR/WMkHTC3tYqI2MKFj2FUb2x3cL2z6qwMW/fUocBq23cCSLoQWEq1m19ExJwa5lZEA98BLpX0ReDhyZO2Ly4pPGxBYyHww47jtcCvdl4g6QTgBIABLlEWEQGMxUD4HsD9PH7GlIGRDBrTNZke9xuyvRxYDrC9NPK/vYgYLaPe0rB93NaUH7agsRbYt+N4EXDvHNUlIuIJRj1oTC5MOMWDwArbl/YqP1QD4cANwGJJ+0vaDlgGXDbHdYqIALasPdXW7ClJ8yV9R9Lf18d7SLpS0h31z907rj29nlV6u6SjtuJj7AAcBNxRP15A1WV1vKS/7lV4qIKG7U3AycAVwG3ARbZXzW2tIiK2mCh8FDqF6m/dpNOAq2wvBq6qj6lnkS4DDgSWAOfUs01n4xeBV9r+mO2PAb8BPA94PfCqXoWHrXsK25cDl5dc+whsuBvumnJ6L2DDY0f339/7RiXX9M/j6zsaUuf+G7X6wvDX+RlbfYcWZ09JWgT8JnAm8Ef16aXA4fXz86mWLX9Pff5C2xuBNZJWU802vWYWb70Q2JmqS4r6+T62N0va2Kvw0AWNJmw/beo5SStsHzIX9ZmNUasvpM6DMGr1hdGsc1MGJiZaW+f2r4E/BnbtOLe37XUAttdJmlzeYyHQmYC3tj43G38J3CTpaqrJRy8D/lzSzsA3ehUe6aARETFoDdoZe0la0XG8vJ79iaTXAuttr5R0eMG9es4sLWX7k5Iup2qpCHiv7ckJR+/uVT5BIyKigQbdUxu6tLwOA14n6TVUA9O7SfrfwH2SFtStjAXA+vr6tmeW/grw0vr55ib3GqqB8JYsn+sKNDRq9YXUeRBGrb4wmnVurI3ZU7ZPt73I9n5UA9z/YPt3qWaLTi4eeCwwOQX2MmCZpO0l7Q8sBq6fTf0lnUU1AH9r/XiHpA+Wlh+7lsZk829UjFp9IXUehFGrL4xmnZuy3e+M8LOAiyQdD9wNvLF+31WSLqL6I78JOKleR2o2XgMcZHsCQNL5VEuLnF5SeOyCRkREP7UdMmxfTTVLCtv3A0fMcN2ZVDOt2vBU4IH6+VOaFByb7qlhXlJd0g8k3SzppsmBsQEl8ZTW7zxJ6yXd0nGucf0kHVx/ztWSPro1K2nOss7vl3RP/T3fVPcXD0WdJe0r6ZuSbpO0StIp9fmh/Z671Hlov+dBmJiYKHoMsQ8C35H06bqVsRL489LCYxE0NBpLqr/C9kEdA2ODSOIp9en6vTrNpn7nUi0mubh+TL1nv+sM8JH6ez6ozvkZljpvAt5l+3nAi4GT6noN8/c8U51heL/nvhvl/TQAbF9A9fu8GPgy8BLbF5aWH4ugQceS6rYfASaXVB9mS6mSd6h/Ht1x/kLbG22vASaTePrG9rfZ0lSdVf3qmR672b7G1f8xn+koM6g6z2TO62x7ne0b6+cPUWUBL2SIv+cudZ7JnNe53ya3ey15DLmXUCURvrx+XmxcgsZ0S6rPNvGlHwx8XdJKVUu7w5QkHrbs0Tssn6Vp/RbWz6eeH7STJX237r6a7OoZqjpL2g94IXAdI/I9T6kzjMD33BeFrYxhbmlIOgf4A+Bm4BbgREl/W1p+XIJGa4kvfXKY7RdRdZ+dJOllXa4d9s8yU/2God7nAs+iWoxtHfDh+vzQ1FnSLlRdAqfa/nG3S6c5Nyx1HvrvuZ9Geee+2suBo2x/yvanqGZTHV5aeFyCxlAvqT6ZbWl7PXAJVXfTfXWzHfU3iWe2mtZvbf186vmBsX2f7c31VMKPs6VbbyjqLGlbqj++n/OWXdKG+nuers7D/j33k4HNExNFjyF2O/ALHcf7At8tLTwuQWNol1SXtLOkXSefU60ieQsDSOLZSo3qV3etPCTpxfXMmLd0lBmIyT++tddTfc9DUef6/p8EbrN9dsdLQ/s9z1TnYf6eB2HUu6eAPYHbJF2tav2pW4GnSbpMUs+/m2ORp2F7k6TJJdXnA+d5eJZU3xu4pJ5huA3wedtfk3QD/U/iKSLpAqrm6V6S1gLvY3ZJRm+jmtW0I/DV+jHIOh8u6SCqfxD+ADhxiOp8GPBm4GZJN9Xn3stwf88z1fmYIf6e+27IA0KJP92awhqDLyAiYiCesddePmNp2cTME887b6WHdNVfSXtTrT8FVWtwfbfrO41L91RExEC0vAnTwEl6E1WX9xuBNwHXSfoPpeXHonsqImJQxqB35gzgVyZbF5KeRrWPxpdKCidoREQUmpw9NeLmTemOup8GvU4JGhERDYxBS+Nrkq4ALqiPf4cGExMSNCIiSg3/dNqebL9b0m8Dv06VeLnc9iWl5TMQHmOhXpF1jaQ96uPd6+NnzHXdYnyMw9pTkv7C9sW2/8j2O21fIukvSssnaMRYsP1DquUtzqpPnUX1L6i75q5WMY7GILnvyGnOvbq0cLqnYpx8BFgp6VSqpvfb57g+MYaGPCDMSNLbgD8Enimpc9mQXYF/Lr1PgkaMDduPSno38DXgVfUy+RGtsT3Ks6c+TzXg/UHqfVtqD9ku3WYgQSPGzqupVl59PnDlHNclxtAwj1d0Y/tB4EHgmK25T8Y0YmzU6yEdSbUr2TunLKwX0YoxGNPYKgkaMRbq1VPPpdrz4W7gr4APzW2tYtyMw+yprZWgEePi94G7bU92SZ0DPFfSy+ewTjGG0tKIGAO2l9v+nY7jzbYPtv2tuaxXjJl6ILyNTZjq3KJvSrpN0ipJp9Tn95B0paQ76p+7d5Q5XdJqSbdLOqqPn3RGCRoREYVa7p7aBLzL9vOoxuFOknQA1cymq2wvBq6qj6lfWwYcCCwBzpE0v/1P2V2CRkREA211T9leZ/vG+vlDwG3AQmApcH592fnA0fXzpcCFtjfaXgOsZstWuwOTKbcREaX6NMgtaT/ghcB1wN71FrnYXifp6fVlC4FrO4qtrc8NVIJGREQh0ygjfC9JKzqOl9tePvUiSbsAX6aa+ffjemvo6Uz3wsBH3BM0IiIaaNDS2NBru1dJ21IFjM/Zvrg+fZ+kBXUrYwEwuffFWmDfjuKLgHvLa96OjGlERBSa3ISppdlTAj4J3Gb77I6XLgOOrZ8fC1zacX6ZpO0l7Q8sptq2daDS0oiIKNXumMZhwJuBmyXdVJ97L9UKzRdJOh64m2ovb2yvknQRcCvVzFvSrpIAAAH6SURBVKuTbG9uqzKlEjQiIhpoK2jY/iemH6cAOGKGMmcCZ7ZSgVlK0IiIKGTAo7vKbSsSNCIiGhjndaVKJGhERJQa88UISyRoREQUMrBp88DHnodKgkZERAPjvIJtiQSNiIhCTvdUgkZERBMTmT0VERElJpdGfzJL0IiIKGVnIHyuKxARMSoMbE5LIyIiSmVMIyIiimT2VIJGREQjaWlERESRyf00nswSNCIiCtnm0cyeioiIUumeioiIIrbTPTXXFYiIGBUZ00jQiIgol5ZGgkZERCmTMY0EjYiIQrbZ+Mgjc12NOZWgERFRyDabnuQtjXlzXYGIiFGyefPmokcvkpZIul3SakmnDaDqrUhLIyKikCcmeHTjxq2+j6T5wN8CRwJrgRskXWb71q2+eZ8laEREFLLNo+2MaRwKrLZ9J4CkC4GlQIJGRMS4mJiY4KcPPdTGrRYCP+w4Xgv8ahs37rcEjYiIQo/AFXfDXoWX7yBpRcfxctvL6+ea5vqRWHM9QSMiopDtJS3dai2wb8fxIuDelu7dV5k9FRExeDcAiyXtL2k7YBlw2RzXqUhaGhERA2Z7k6STgSuA+cB5tlfNcbWKyE/yrQsjIqJcuqciIqJYgkZERBRL0IiIiGIJGhERUSxBIyIiiiVoREREsQSNiIgolqARERHF/j/hm8pgOEhMOwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "slope = np.deg2rad(10) # degrees\n",
    "for i in range(rows):\n",
    "    dem[i*cols:(i+1)*cols] = i * dx\n",
    "\n",
    "dz = np.arange(0, np.sin(slope)*(dx*cols), np.sin(slope)*(dx*cols)/cols)\n",
    "for i in range(rows):\n",
    "    dem[i*cols:(i+1)*cols] += dz\n",
    "    \n",
    "imshow_grid(rg, 'topographic__elevation', colorbar_label='topographic__elevation (m)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial surface temperature `Ts_init` needs to be extrapolated. Assumptions:\n",
    " \n",
    " - `Ts`is sinusoidal\n",
    " - lapse rate `lr` = 6 °C/km.\n",
    " - Provide the temperature mean `Tm` and the amplitude `Ta` and the representative elevation `T_ele`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fec017b2be0>]"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO29e5Bk2V3f+fnlu97V1V39nJmemZaQGCTNIJrhIa8sWYidEbJkYm0sLQsKYD0IIwxeIoyACOyNDSIMtjFri0UrQAYCWYQWEAgzehnsEDikGfU8JI0YaTTTMz3d0696P7LynWf/uHmzsrLyce89v3NvT9X9RnR0VVZm/u7vnN/5nd/r/I4YY0iRIkWKFAcfmaQfIEWKFClSxINU4adIkSLFIUGq8FOkSJHikCBV+ClSpEhxSJAq/BQpUqQ4JMgl/QCjcOzYMXPnnXcm/RgpUqRI8bLBY489tmyMWRz0t1ta4d95551cuHAh6cdIkSJFipcNROTSsL+lIZ0UKVKkOCRIFX6KFClSHBKkCj9FihQpDglShZ8iRYoUhwSpwk+RIkWKQ4LACl9EPiwiN0XkqZ7X/pWIvCQiT3b+vW3IZx8Qka+LyLMi8n6NB0+RIkWKFOEQxsL/XeCBAa//e2PMfZ1/D/f/UUSywG8ADwL3AO8WkXuiPGyKFClSpIiOwArfGPM5YDUCjfuBZ40xF40xdeAPgXdG+B4VXNuo8JFHLlGuNWOla4zhz790lcdfXIuVLsCzN7f52Bcv02i1Y6XbaLX52IXLPHtzK1a6AE9eXufPnnyJuNt/l2tNPvLIJa5tVGKlC/A/nl3mL5++ETvdpa0af/CFS2xUGrHT/tRT13nk4krsdC+tlPnooy9Sa7Zip20DjYNX7xORHwYuAD9rjOnXaGeAyz2/XwG+Y9iXichDwEMAd9xxh8Lj7cVPf/RJHn1hlWeub/F/vvM16t8/DJ966jo/9dEnyGeFC7/4VuYm87HQbbbavOfDj/LSeoWNSoN/8sa7Y6EL8Lv/4wV++eGnOT1X4q9/7u+RzUgsdDd2GvzA//t56s02xVyGB15zKha6AL/yqa/x+5+/xP13LvCx935XbHSfW9rmB3/7EQD+/H1/h9feNhcb7Z/74y/zV1+7yeMvrvFrP3BfbHT/5hvLvPcPHgPgkV94CydmS7HQbbcN/+T3L/DMjW2ubVT5P976TbHQ1YBt0vY3gXPAfcA14N8NeM+gVT7U7DLGfMgYc94Yc35xceDp4Mh4cWWHR1/wnJSPP/ESrXZ81t8fP34FgEbL8Om/vR4b3UefX+Wl9cqeZ4gLPr2rG1UeeT4+K+xTX71Gvel5M3/02Eux0W21DR9/3KP36AurXF7diY32nz6xy+efPBHfPK9s1/irr90E4C++fI1qIz6L948e27Uj/8uXr8VG96mrGzxzYxuAP34s3jVlCyuFb4y5YYxpGWPawG/hhW/6cQW4vef324CrNnSjwlf2P/Gmc2xWm3zt+mYsdNttwxdfWOMHzt/GwlSBR5+PEhmLhkeeXyUj8ONvvJuvXd9iYycet3t9p87Xrm/x42+8GxFi5/noVIEfOH8bX3xhlXZMG/vT1zbZqjX5iTed6z5HXHjk+VXuvW2O7z53NNax/uILnkP/E286R63Z5kuX12Oj/ejzq3zfa09xx8Ikj8ZoUPjj+96/e46X1itdg+rlACuFLyK9vvL3A08NeNsXgVeKyF0iUgDeBXzChm5UPHZpldlSjnd/uxcqeuLFeITzuaVtNioNvv3OBV5/xzxPxBjHf/zFNV59cpa/+yrPW3rySjw8+2P75lcf51UnZmIba5/2t509wvmzC2xUGjy/Uo6Frp+f+V/vv4OZYi62eW622jx5eZ3zdy7wbWeP8PS1TSr1eCztx19co5DL8CPffScAT8Sk8K9vVLm6UeX8nUf4trNHYpWvxy6tcfvCBG977UmAWNezLcKUZX4U+DzwKhG5IiI/BvyqiHxFRL4MvBn45533nhaRhwGMMU3gfcCngaeBjxljvqrMRyB848Y2rz41y+0LE0wVsjx7czsWuj6dbz41y6tOzvDCyk435BAH7VefnOHVJ2cB+MaNeBKo3+gkal99coZXnZyJbaxrzRaXVsq8ujPWQKzzPFPMcduRCb4pRp6vrFWoN9vdeW4beH45nk3uGze2OLc4zfHZEidmi7HxvCtf3jzf3KrFljT21tQs33RiBpH45EsDYap03m2MOWWMyRtjbjPG/I4x5oeMMa81xrzOGPMOY8y1znuvGmPe1vPZh40x32SMOWeM+WUXjATB88tl7j42hYhw7vg0zy3FM1EXO4vvrmNTvOL4NK224VIMVudOvcm1jSp3L06xMFXgyGSe55biUQQXl8osTBWYnyzwisVpXlqvsFN3Xxl1aWWHtoG7j01x7vg0EN+CvLhU5q7FjnwtTsUoXx6duxenOHd8CoBnY5Ttu495NM8tTsc61uDx/IpFb57jGG9v7e5w97EpSvkstx2ZOJgK/+WOjUqDlXKdOzvC+YqYhfPkbImpYo5zi/Epod1F4dE8tzjNc3Hx3KsIOor3YgybTa8imC7mODlbio/npe09ym95u876Tj0Guh2ej01z59EpMjFZnfVmm8urO9y9uMvzcze3YymFfX65zFQhy/GZYqwb+5W1Heqt9h6eU4V/C+KFHisb4PaFSa5vVmMJrVxc3u7SPXvU+//ymvsKjot9PJ89OsWLMVWOXFwqdxfF2aOTALFUrfjW7i7Pk7GMdaXe4upGtbu5+obF5VX3Cb2Ly2XmJ/McmSpQymc5NTfBlRjG+sXVsudNdeb5zmNTbNWasYRWnlva5u7FaUSE245MkBFi4dlfU915PjrF5dWd2M97RMXhUfidEMqdHYV7Zn4CY+DGZtU57curla7Smy3lmC7muLrunu6VjqK7Y8GjfWa+xM2tKk3HB7CqjRbL27Uu3dNzEwCxVDNcWauwMFVgpuSdczg9PxHLWPu8JcWzTxfg9HwpFrqX1/p59urg46D9Ug/P+WyG4zMlXoplTfXxPF+iXG+xWY33IGdUHBqF7yv2kx2hPD0fz4JstNqslGvdQyEiwun5EldjWBQ3NqrMlHJMFb3zdafmJ2gbuLFVc0u3M9Y+z/OTeSbyWa5tuF+QNzaqew7gnJ4vcX2z6vzMRT/Pp+e9/+M4cbuf54nYxhp2eT7VWVPXHCteYwzXN/fPc1xjnc0Ix6aLHbodnhM4WR0Fh0jh15jIZ5ktecrPX5CuFe/NrRrG7G40AKfm4lmQ1zernOxTBOCe5+sbezdXEeFUTAvS47nY/f3U3ASttmHJ8SbXz/PCVIFiLpPIPHvyVXF+/uB6Z5M7PhPvJrdVa7JTb3FyrmeeY9rkrm9WWZwudk+Nn5qLZ5PTwiFS+FVOzBYR8SYqduXXp3jjsPCvb9b2bDS+y+2c5839PJ+Zn4jF5b6xWd3D85mYPLl+nj1PbsI53Uq9xUal0cdziUbLsLzt3pM7Nl2gkPPUyLGpIvmsOJ/nfs8CvHm+ul5xHku/sVnlRALypYVDpvB3J6qUz7IwVeCqY6vAd/WP91idp+dKrJTrzo+h97v6p7qbXDw8n9jj1ZS45nhR1JotlrfrfTzHY3Ve36gyW8oxUcju0o6B50Gbq291upbt633ylckIJ+fce3LXBhhRp+ZK1JptVstuq6Kubez1IBdniuQykoZ0bjXc2Kzta650bLrASgxWEOwVTl/5u7TAWm3D0nZtD93pYo7JQta55Xd9o8ZkIctMcbc33/EZb5NzGWa4uenxtWesO+EG5yGdPs/Co11kyflY7w0lwa58ued5r3yBN95xjDX08ezPs+v1vLE3fJbNCEenC8551sKhUPh+kqd/QR6bLrK87dYiuL5ZpZDNsDBV6L52dMpX+O5oL2/XaLXNHisbfJ7db3InZ0vd8JlHt0CrbVhzWJc+yLOYn8iTzUgsPO83KIosb7mVr/5ksU8XiMWY2S9fBfdjPSCkc2zaW18ux7tca7JVaw5ZU+7PW2jgUCj8jUqDerPN8ZnintdjUX4bVY735A4Ajs24X5CDPAuIaUEOUH5HfSXk0OW+McDCz2SEhakCK44X5I2+xCl4PFcaLacnjPurz4CuceFynmvNFqvl+kCenY/1VpX5yTyl/G74bFe+4l9THs+phX/L4GbH3To+0AJzO1HL2/WuxbVL1/2C9BedT8tHHAtypVzvbmo+/DFwOd7+Yt8/3m43dmMMK9uDeHZvda6U65TyGaZ7wmelfJaZUs6p1blW9g5XDRrr1Z2607MeKwPW1OK0+zCWb6wMWs+phX8LwU/kHJ3aq/yOzRQo11tOOwuuluv76U67D+n4PC8MoO3awl8t11nou+Clq/wcWvg+z/MDaLsc661ak2bbsDA5ZJ4dWp3eWBf2ve56nnfla/9YGwNrDttwD+J5diJHPitOPchxa+rlcNr2UCj8tc5EHRm2IB0ujLWdOkf6BKSUzzJdzMWyIPtpH5susFquOzuI1Gy12ag0BtB1b+GvlevMlnLks3vF2rXyWxs61vHw3E/Xo+02jOXnYpJbU3s3GhHh6JRbj334PBeoNdtsx3xtahQcCoW/ujN4Z150LJzGGM8aGbIgnVr4O3XyWdlTKePRLdI2OEue+pZdv1czF0PydHWn0Y3l9uLoVMGpBTbMgzzaDd25nedB8nV0yu0mtzLE2j0aQ/5gtdxgYWrAPDvOT3X1yGQ/z+49di0cCoW/NsTVd70gK40WtWZ7nxXk0XZvjRyZLOxJFnt03S7IruXXpwgyGeGo4+TparnGkQF3BR+bKVJttNlxFLob5k35Y+0yoTfUoJgpOA1vrA0Lb3QLEtzQbncqvfpDSeAZM05DOtt1JvLZPWctIJ4iDC0cCoW/Wm4wVcjuyerDbmZ/1VGMdVicEzout+vY7kDPwg8zuFkYXZ4TiSs3RvPsiPYwnos5r5WH69DdIIPi2HSRNYfJ09VyHRHPc+unC+7GeqvapNU2Q3l2aUQN86biKMLQQpgbrz4sIjdF5Kme1/6NiHxNRL4sIh8Xkfkhn32hczPWkyJyQePBw2BQHB28Gm2AdUcJJr+SYZBwHpksOKMLHZ4H0PUFdr3iVuEPGu8jU3nWHbbOXRui/Hyr39k8d72a/Rv7kamCM54brTZb1eZAJXRk0kueuuriuFquMzeRJ9eXL5kt5chmxNlYDwvPgjfPzuVr0Bx3ZM7letZCGAv/d4EH+l77LPAaY8zrgGeAnx/x+TcbY+4zxpwP94j2WBuyM08WshSyGWdC4gvn0en9tOc6wukqrrwyxMJ3vckNi2d7tAvOLgQxxngW2ICx9kN5zua53KCQ3Vsa2aU9kY9hoxnBs6PxXt0ZXB0kIsxN5J0bFANlezLPTr1FrekodLczOHfgWr40EeaKw88Bq32vfaZzZy3AF4DbFJ9NDcMsPxHxFK8zC39wJQN4yq/ebFNtuHG514Yo/NmOwnd1ScVuvmTwJueK7k69Rb3ZHqiE5iZ8C8xRorpj+fXnSwDmJt1Z+L4HOZhnb55dlUcOqw4Cx5vcCIU/1xkHl7LdX24MMJHvGI4HzMIfhx8FPjnkbwb4jIg8JiIPjfoSEXlIRC6IyIWlpSWVBxsWewNfOJOxRsBNaKXVNqwPKI0EryR0Ip91xvNKuc5MMdftoNgLXxG48GpGhZL8sXalCFaGGBTg8bzhbKy9mPGgMMN8V/m5k+1hPLvc2FdHGlGdeXbovQ6SL99wdDXWmlBR+CLyi0AT+MiQt7zBGPN64EHgJ0XkjcO+yxjzIWPMeWPM+cXFRY3HY63cGL4gXVr4O3UyArOlAQvSYWhlo9LAGAZaI+Ce56GW32SeZttQdlAtMypZPOc6VzPKoHAYV/Yt/KODwgwxhO4Ghe182knE8F2GVmrNFtu15kD5Arc8a8Ja4YvIe4C3Az9ohphuxpirnf9vAh8H7relGxTdiRpgBYFnCbmL7XpWUCYzyNV3tyBXu5bfEAtswp0SGlYdBF4YC9yEVrqKYEAMP9+Jr7sMM4zyIDcqDSddQldHJIvnHcqXMWbMxl5wFsNfK9cp5jJM9pVGQq986fPsf+cg+QK3RpQmrBS+iDwA/BzwDmPMwBuERWRKRGb8n4HvBZ4a9F4X8CdhdLzR1SGk+r7a/1267lxuP247yqtx5fau7zQG1sKD203On8OhYQaHicRhFVHgxZWN8coJtbHu50sm9tOeKeURcWPtlustGi0zfJ4dJ6oHnS8Bt4nqYSeLfcxNuDMcNRGmLPOjwOeBV4nIFRH5MeADwAzw2U7J5Qc77z0tIg93PnoC+BsR+RLwKPAXxphPqXIxAv5EDVoU4HZn3qg0BiYvYdcqc5FU2+wI3qjNxpXy26w29tVm+zjiMKm2WfEU6lDaU242OWMMm9XmCJ7d5Wo2Kg2v0mxAviSb8aplXOQPfPkaNc9b1aaTMwCbleFj7TJX48vOqHl2lavRxP46siEwxrx7wMu/M+S9V4G3dX6+CNwb6ekUME4RzE8WqDRaVButfQezNGj3d6vs0nXofm5Wve8clDsAt5vcZqXRrQQaRBcc8dxZ5DOlwSI978gCK9dbtNqG2YkhdHt4PntUl/ZmtTF0jqHjvbrYXH35GjPPm0POCNjSHjbW00V3ZwD88wwj19RBsvBfrvAVwTAh8TeCTUcLY9iiKOUzFHIZJ5afv8kNo+3qDIBv7Q5dFBPurN3Nqmft9jdO8+GV37qzdofx3C0JdeTVDJNr6JSEOtlcxys/cBNaGbXJiUhnk3M4z0M39oLTMwBaOPAKf6s2ekH6YQZXoZVxwukizBDE2nVxBmBnjLU767ByZLMyfKOB3eSpOt2A1m7cyg8cWvgBjSh3m9xwnl2dqxnnNfs8uypH1cKBV/i+NTJU+TlakF1rd4QF5iq0Ms7adXUGYNyi8M8AuFK8QcZa26sZa+06VASjPEjwk/NuNhoYZeF3cjWOZHt2yFoGhxt7QD3iqhhCC4dA4fvW7uidWdsa6Vq7Iy2wgpM2xUGsXdC3tMeFksBXvElYuwUnZwACW7vOvJrRys+thR9v6M4YMzJHBJ2SUEcbzVQhu693UJeuw9CdJg6+wq82vKPPAyoZYNdK0S6bG+fqe3/LOWlutVltDLVEwJ0SGmf5+bSdKb9Rrv6EG09uHM85h2cAxln4c47OAPgyO0zGXMlXud6ibQLIl6MYfjD5ShV+ohiX2PKFdquqO1FbY7L6Hu28Ol0YrwhmupuctoU/2tr1aOec1KSPc/V351l5Yx9j7fq0tce6a+2OkS9joKx8ifpmZXTIcNr5WCclX6PnGPTXlDYOvsIfM1EHUjjHuPrOlF8AC3+mlO8m0lVpj7HAZpx5cqOtXf9v2nS71m4gY0Z/nkfNcT6bYSKf1TcoAsmXN9YucjVJjLU2DoXCH7UYfeHULssMKpzbNQfCOdbCd2ONBInhu1B+48pBfbrgxqsZZe16tPU3uXHloD5dcGHMjFZ+Hm39eQ4mX3labUOloZyrGWvhu/GatXHgFf5WdXRsF5IXTu2r94K4+uDOq4nb2h1XDtr7THFbuz7tJHJEzja5pHgOtMk5nOcRY13IZSjmMqmFnzTGKT/oCKe2Bda18ONVQkHKQbvCWdNfFIGs3apueWTQUBK48WrGW7t5dwbFLaj8PNr57rxo0oVxYSyH8zxiLfu0Xd0wpoWDr/DHKD/wrCR31m68SihIOahP282iGL+5NlpG9dBX0FAS6F/5F9zaTSJB7rc4SEL56VeghbHwNWm324atAJvcrIN51saBVvhBKhnAzc68WW122ycMw6wD4Qzi6vu0nSi/ANYu6G5yQSz8UudWomSsXTdjDaN5nk3Qwp91YVAESJC74Llcb44tB/WfKw3pJIhKo0WzbQLG8PUtsCAbDSgrvwCuvkfbTVx5HF0nm1wAaxdczfN4a3e2lKfebKv2WRnXsRLc5Gp8I2oUXY+2mxj+qMNPHl0XBoXvQQYJ3aUWfmIYdxzax2wp132vGu2ALiDoLsggcU5wGNIJoAggfgvfp52UhQ/a8zxetkv5DLmMqI51kMNP/nM5SRYnMdYBQkk+7dTCTxBbgRWBq3h20PDGwRDOcYefPLoueB4fw/dpa85z8JChGyU0ztoVEfV5Du5N5ak22jQUe+IHyxG58JqDhUlThZ8wgsazZ4o5as029aaicIawRjSTasFj+C42uWQtsHGenPaCDHL4CWCm6CZvMW6sQX+TC+NNgb73Om6spwpZMuLGm0rCcNRGmBuvPiwiN0XkqZ7XFkTksyLyjc7/R4Z89gER+bqIPCsi79d48CDYjWePVwSgbxWME5DJQpassssdhmcn5aBJWGABykE92o6s3UQs/PFjDd6lIEl5U+BgkxvDs4g44Dl4jsi/EOdWRRgL/3eBB/peez/wl8aYVwJ/2fl9D0QkC/wG8CBwD/BuEbkn0tOGRGAL30WYIUA5qEvhHFUO6v99p95Su4bOLwcNYmWDvvIbR9ej7cbaDTLWoK/8gvHsSr5u7Xl24TUHneftWzisE1jhG2M+B6z2vfxO4Pc6P/8e8A8GfPR+4FljzEVjTB34w87nnCOocPobgqaQbFUbTBeDuNy6C3K71qSYG10O6tP1369FF3Z7Ew3DdCGHiK7y2641mS7Gr/z8hT2OZxe14du15li6Hm1d5ded5zHj7SJcGZxnN/M8VRx9BaoLnrVhG8M/YYy5BtD5//iA95wBLvf8fqXz2kCIyEMickFELiwtLVk9XPDYm641Umu2aLTGW7sebV2rM4zyAz2egyqCTEaYLujWpQfnOc92vanWLniX59GKwEUL7u1ak6kAPM86MChg/Dwny7Pymqo3KWQzFHPj5vnWb6AWR9JWBrw2dMUZYz5kjDlvjDm/uLhoRbhca5LNCMWA1q6WkJRrXr31VGH8pejah3LKAReF9inMckcRTBXit8CC8jxbymGMt4B16HbmeQztaWX58mg3mQ481rp0YTzP2gZFo+UVVQTnWVu+gqzlW7+Bmq3CvyEipwA6/98c8J4rwO09v98GXLWkGwjlWpOpQhaRQXvOLma7yk9LEQRbFB5tbQusFZgu6Fv4QTcbba8mGF1dnoNuctmMMFXIKiuhYPM8U8qrdmQt15qIeAUH4+iCphEVRr60FX7QsT74Fv4ngPd0fn4P8GcD3vNF4JUicpeIFIB3dT7nHNu11i0d3vBo6yq/cq05NsTg0wVN5edZu0nE0sv14CEd0FNCSc2zMabDczAPsm1Qu9pxu9ZiqpAba0QdpDUVJmQIOLnvQQthyjI/CnweeJWIXBGRHwP+NfBWEfkG8NbO74jIaRF5GMAY0wTeB3waeBr4mDHmq7psDEZQV9+fzANhjdTDWrvaPAdTQpqLwrPAgtEF1E5Vh51nLbo79RbGBPemYLeAwRZBwxv5bIZSPqMfJg061oqXoAQPk+rKlwuM56IDY8y7h/zpLQPeexV4W8/vDwMPh346S5TrwbL6uY5w+gvYFmHCG15ZptcueJzVFJT27QuTY9+XtAX23FJZha5PO4zyU7PwO8m8cRVRHm29TS7sRgOK8xzQoPBo59WU33Yog2L3EpTJADH/cSjXmsxNFgLQdXP/gCYO9EnbrWowVww8RbVd03F7w4Q3pjsut1a74KDJPH/RapVlhlFC06Wc2uYaJpk37YDnIAoIYLqUV5OvMJurb/Bo8hx0Tc0Uc4oJ8uR43g4YJi3mvI6sWvPsAgda4XtJ2+AKX8/C93b4IMpAXwkFSzAVc15jLS2egybzwN9cE9hoOu8pK27sQa3d6WJWcayDhzd2eVbc5AKuqSnFNRVunrOdzyjOc2Ce9ebZBQ6+wg+4IDWFczuEhe8LkgbtMMk8EVHnOUgyDzyea822yinfMNbuVFcRaFp+AeWroDnWwcMbmvLl0Q6+yWkqv1DzrMxzUnrEBQ60wg/qioE3UVpX/oWxRjRDK2GSedDJHyQQ3phStMDC5Et8RaAb0gmuCLS9mkDhDQdhrKBrSjNMGjYv1vsZG+waUWFCw6nCjx3eRIVxuXXdz2IuM7aZl0/X/4wtwiwKn7aaBRYimdddkArx3TDVQZmMMFnQDK2E47msVA9frocxKHS9mqSs3TDzPKW4pioNryNqKJ6V8hYucGAVfq3ZptU2CYV0QlgEnQSThpCEcXvBd7m14pwhwhuqm1zw8JlPW2tBhvUgtZLz4cJYvnzpWdqhxloxlJTPytj2Bj5d7zOaayr4PKdJ2wQQVvnpVumEsfw8QdKgHSaU5L9PrZIhZEUU6CzIbngjQPmtT1uzFDU4z/48KyihanDZ9jxNUaHbaLWpNduJhDfCGBSa8hW0Qd4u7TRpmwjCLArvfZoJpjCJLRchneAxVk2vJvAmV3LAcwLVWKGqdBR5DlMRpZmcD2tQTBe95LzGrVehjCjVsfb7YsUvXy5wcBV+BGu30tDpDx8msdV1PxWsTl84ZwK0ZfZpqym/EImtbvJUhedkwlhhk3maCeMwFVE+bY2xDh8y1N3Yg9KdzOt5zVF4PhD98F9uCKsIuslThVhn0PYGoKsIwiS2QNvlDtbewKcL2jzHG2aIUhEFejwHHWuftg7d4PX/Hl29MFaYNZXpNKtLyqsp1/XaOmjj4Cr8enjlB3rWSFAByWaEibyOcEZJ2mp1UgzDs2blyHatFbi9gUdbJ2kbJV/S+zkbhKmI8mhnVYsCgpff6h10CxMm9WmrKPwQFVH++9rGq+65FXFgFX6U6g3Qi7EGOervY7qUnBLSqBwJ097ApwtK3lRIa1dLEUSp3uj9nA3CJDB92ppFAWHXlB7P8Xs1UcNYt2ot/oFV+FFcMdASznDWiFaFUJhkHni9TsCe57Bj7bd10AtvhB3rZJJ5vZ+zox28vQF0LtdOKLzR+zlb2mF41k9UB40U6LZ10Eaq8DvQ2pnbbROqPtujrRXSCZnMU1qQYa0gzcqRMMk88BR0tWHf1iFK+Az0wlihwhtKbR1C86zY4iBMyBB2w5X2dMNt7NptHbRxYBX+brle2OP+dhO10wiX2IJOFUUCyTytTS5sMg8ULe0I8Wywt8BCGxTKyflwBoWefPnfFwRaXrMxJnQYS+zl/qQAACAASURBVNNrnixkyWSCGVHarSy0Ya3wReRVIvJkz79NEfmZvve8SUQ2et7zS7Z0x2G72mQinyUXoL0B7JYy2gpJ2EUBHeHUKJsLqfy0XO6wyTz/vSqWXzUaz7ZtHcIm8zTbOkQJY2m0dfBzLuF7JtnxHLa9gUdbL6QTlq7/uVsRwTkZAmPM14H7AEQkC7wEfHzAW//aGPN2W3pBEd3yize84dHWS9qGpQv2bR3CJvP892ol1c4cmQhOV+lQTpR51uQ57Dz7lSM2F4Js15qB2xuAXl/6sAly773JhAy1e/FrQzuk8xbgOWPMJeXvDQ3vPtskwhsRFIFSUm27Gi6x5Y+PbauBKF6NXuVI8F7lPl3QTFTHWznS7LQ3CGfh69TDh7V2i7lsp62DrdecYMgwwpkHODwK/13AR4f87btE5Esi8kkR+ZZhXyAiD4nIBRG5sLS0FPlBwgunTuVI2BO+oGv5RXM/7RZkVGs3qfCG/zkbhE3mgU6YIYry05znMPz6tO15jmZQaNy5ENWguFVDOmoKX0QKwDuA/2/Anx8Hzhpj7gX+I/Cnw77HGPMhY8x5Y8z5xcXFyM8TVvmJiIqlHeZ6Qx9alSNBLz/p0lUSzihejYYi8NsbzARsbAV6VRRemWDwZB7otHXwcw8zCWxy5Vq4sQadCiHfoAjDs+YmF4ZnzbYOLqBp4T8IPG6MudH/B2PMpjFmu/Pzw0BeRI4p0t6HsPFs0KmWieLqd/MHlgeRwtb/a1WOJOXVREnm+TJhG8YKmyz2aVuPdTXaWIPOPCfBc7RCiE640vLi+LA8a7Z1cAFNhf9uhoRzROSkdIrDReT+Dt0VRdr7ENbVB50wQ9TwBugkEsPQzSpVjoRtbwC7VTo2lSNRNhq15HyIxmm7tO2T89EqovTCWGHXlEZbh2jzrGPhR9Ejmk0JtWFdpQMgIpPAW4Ef73ntvQDGmA8C/xD4CRFpAhXgXcZxdyEvaRtBOBMoy9RIJPrtDSIJp0KVThgF5NP12zpMBDwrsZ+uHz6LoPysvalkFEHU8BnoJG3PzJdCfWaqmGPTuiggQphU0asJq0du5WsOVRS+MWYHONr32gd7fv4A8AENWkGxXWuEUgSgJZzh2huATilXFEXgv1+jSieKNwUez9EVvn+4Lv7kfJRNLrnwhqK1GzJpO13McW2jak0XwldE9X42CroVUQkkql3hQJ60bbbaVBvhrV2NniNbnUURtL0B6AhnlFAS6ByAimoFQfw8ayXno3iQ00X75Hy0cx56Z0ySDJOGKzm2l6/diqjwG3vaSydG+O56lKStSvVGWM9CoXIkSrmeT9va8gt5yA10XO4o1i7oJeej8mwz3pFChh352rLgOUp7A+ict1DwIMO0NwCdRLVfERWFZ5uxdomDqfCjKgIVlzt8YmtXOKMrgijJPJ+2fZzThmd7Cz8JqzNaGKtTsmeRMwnb3gB0KkeqjXboiijQuRAkimehkaiOqkdu5XttU4XfA42eI1HCGxoud9QYvlbSNkq+xP9sVNiEsZJI5mnxHKa9QS9tnbGOlpy3uRDEZk1pGBSR1lSq8ONDVOGcLtlfCBIpsaWYtE2qciR8Ms9+QUZJ5nnvt2vrYJPMA3uew84x2HtyNtYuaPAcbo412jrYGI63apXOgVT4YS+n8OFPrM1hjSju565wRheSrYjWyExJI6STTDw7SnsDsA/p2CTzvM/bWZ1h+QU9Cz+ZvEW49ga9tHVCOuENCo22Di5wIBV+18IPeQxc47aasO0NfOgJZ/gEpk3liE0yD+zDWGGTeT5tK+UXNZmnkpwPP9Zgf8bEJmTY+/koiBLSAfsijLDXpHbpKpXBusDBVviJLMhW6I3Gp51MeMNuk4uazNNo6/DyC2/YJ+e9ooDwBoU1zyH7//fSBct5jlAF5tNONIyl0PJcGwdS4dsvSEuXO6Jw2lojUZJ5theCRM2X+G0dbMc6iuVnm5yPalB0czXVeEOGYK/8tiOcavbpAlalmZE3dstwZeR59i9TUrjUSBsHUuHbZNchuoXvtzeYjhRvtHe5oywKW56jhs982smEN+wqR/yFHJZnjSZ5Ybs37tK2NCh8njvKLAxdsLtkJzGea01yGaEYokeUR1fn/gEXOJAKv1xrko00UXYWflTPwv+MdXgjwkZj69VEaW/QS9uO56jhDbsFGZVnjeS8zTwnETK0la/uqflIPNt5kL4RFebUvEfXPjTsCgdW4U8VshYTFc0C26pG8yz8zySS2NKy8CMnEpPkOaKFb8Wz/TxHNShsKkeitDfw6NqdMYly0KxL2/IUeVJryiUOpMKP0ucENITTzsK3Cm/Uw9cqe3QtebbxaiwXZNRknu2CjGrtgl1yPmpFFNhvcpErorrJ+eh0IZnNNUr9P9za1xweSIUfOZ5t2XPERhFoJNWiJvP8z0ejG13hW/Mc4RISny5YhHTq0foW+bSjKqGoFVEeXbvKkaiba8byzgVb+bJp6xClTQqkFn7sKNebkZKItj1Hotbt+p8p11sWwhk9seV/Phrd6DzbtnVIMowVJZnn0Y6enI9aEeXRteU5mtfs004mZGiZnLc48wD2dy64wIFU+FvVaBMFdsJpm7RttU3ktg6JJ21jPmwWtb0B6CRtoyTzwC45bytfYMtztHsLbDw5G55V5jmCfGkk511BReGLyAsi8hUReVJELgz4u4jIfxCRZ0XkyyLyeg26wxB1osBOOG2sEVvhjJrM8y8EsXa5Y66iiNrewPuMfdI2qkFhE9KxDW+AnVcTdU3ZJOdtDQrvO6LnD6KMtU/7VgzpRONmMN5sjFke8rcHgVd2/n0H8Jud/50gqYnSsMDKtSaLM8VQn7VJ5omItdUZJZkHHs9+W4dcNpztEbW9gU8X7JK2Ua1dDfmK2mag9zui0D45G+56w17a0TdXu5Ah2G1yUcJnYH9y3hXiCum8E/h94+ELwLyInHJFzGqiLGKsGtZIFCGpNFqRk3lg6XJHTOb5dL3vCD/ettVBYOPqR0vmgf1Yg52Fb1MtkwjPCl5NFNrGGMp1u3mOutF87MJlfunPnor02XHQUvgG+IyIPCYiDw34+xngcs/vVzqv7YOIPCQiF0TkwtLSUqSHOXt0ituOTEb6rF1IJ1p7A58uRLNGbJJ5YOdy2ybzwJbn8LSzGWEib1c5Ep3nbOTkfNT2Bj5dsEvaWnnN1q074rXwa802rbax4Dm64fjF51f57N/eiPTZcdAK6bzBGHNVRI4DnxWRrxljPtfz90H+/kCJN8Z8CPgQwPnz5yOVrPz5T/2dKB8D7ITTNpQE0Y6gR73esJe2XZwzejLP/44odMGSZ4t5jhze6CTna802pXy4cUs6aRvVoLC5QzjqqXmwy4vZbDTgjfdmxF46OxaexTioWPjGmKud/28CHwfu73vLFeD2nt9vA65q0NaGbYw1ejLPF854wxsebbtEtU2C3P+OsLAJn3m0s4mFN8CW53iT8622odJIKIxVi3ZqHuySthpryi5BHk2ux8Fa4YvIlIjM+D8D3wv0B6A+Afxwp1rnO4ENY8w1W9ouYCOcW1aufjLhDbDrGx61/h/srE5fWc+EbOblY7qUi9y1MmozL7DrHmlTESUikbtHli0S5IDVnQvbtRYzpYhz3JWv8PNsu6ZsFP6ORV5sHDS+9QTw8c4OnAP+szHmUyLyXgBjzAeBh4G3Ac8CO8CPKNB1gl7hDFs5ohLSSSq8EXM5qEc3elzZ1sKPWjlijOnwHL1KByJuclWvIioboSIKoleObFv0iILeeW4xNxllTUWfY4jmNdvzbJcLPDNfiPTZcbBW+MaYi8C9A17/YM/PBvhJW1pxwFY45yajTZRN5Yht0tauHl4jvBFhQSq43Nc2qqE/Z9PewKcLETc5S8svqtWpEd4Ar5R2bjKctW5jUNi0dditiIqen/LvXAgbjtqJeGteEBzIk7Y2sLkQxKYc1KZyRCVpG7lyJLkwVtT2Bj7tKElbjWQeREvO21REebSjVY6o8Rxxnu14jrbJ2dT/+3Tbhkgn58u1JpO3ctL2IMHvwRNV8UZNYILvBsafYOqtHAkDm17lYJ/AjNreAKIrApv+/2CfnI9qcUL0MIOtQWE9zxZrKmpOzt6rsasQstnkRiFV+H2wibHahDcAZiKWr9kk83y6vd8TFDa9ysGrHMlGrByxXRRRw1i2oSRba9dW+VnJl2XeIuoGa7Omop4x0TCier8nKPx+WjbzPAqpwu9D1BirdzLP1v2MpoTKtSYTebtknv89YelCdLdXJHp3Ug1rN0rliC3PtsovifCGPc/Rk/M2YVKInpzfNaLiTc7b5g7GIVX4fYiq/GzbG/i0o5bNRWkH3aUbVTgtrSDwXe4o4Q27wylR2zrYLkib5LyttRs5vGHR0sGnC+HDWLbtDXzaUce6lM+ErtTrpet/T1i6YLemRiFV+H2IKpy2lTI+7ShWkE07aJ8uhD+gYpvMA5ukmr21C+EXpG0yzyY5b9PeAKK3dUgqaWvb3sCnnUyCPFpyPlX4MSOq+2mb2PI/m0x4Izmep0vRFqRtMi+qEtJYkDbJeavwRsTkvE17A4ietE3SoLD3prx52gp5wM6mX1IQpAq/D0mGN6IrArvqIH9Bhb3a0TaZ59NOJrwRrYpCJ4wV3sK3bW/g0Y0q263I7Q0gelsHrbFOojooalsH2yqwcUgVfh984UzCGomiCHzaSYQ3bJN5EL2tg017A58uvHzCWLbtDSB6fspWvvw7F6LQBfuxjpKc367p5MWi8pyGdGJCVOHUsvArjRatdrgYq+0JzCSFM0qnzt1knl14A6JZ+DYVUT7t0AZFVWesIUJoJeJl8b2Ikpy3uS+5ly5ES85rbK5R5zmtw48RUcIMWklbiJbosasOihbe0PJqwtLVSOZFraKwTZz6tKMm85JIzttWgUG0enjbfkneZ6N6r3bzHDU53/XkLMd7GFKFPwDRhFMnaet9VxSXO/qiyGUzlPKZSAvSpr0B7Fq7YSpHthSsoOmIh81sxxo6BkXoZJ4mz+G6R2qc/Izi1WwpbnKhaVcV5jlCQYLGPI9CqvAHIEqYQSuk0/tdQdBtb5CIy23X3gCiVY5oJLaiJzDtwxtRkvMqFVER2zrYJjA92uEVvtY5D4g4zyo8hx9rm4qocUgV/gDYhHTslFD4BenHJjUssCjhDVu6URakRu4galsHm+6NPqIk5zUqopIKb0C05HxSRpRGRZRHO1qkwKYiahxShT8AUQ5AlWt2vcohWhWFRmzX/3wUnm2PgEdZkBo8R23rYNveAKIl51UqohIKGfq0I29yFrc/RTljolERBdFOzrtsnAapwh+IqGVzGq4+hLN2tU7mRYmxavAcpR5eq9+ITRjLli6ES87btjeAaJUjxhglnsMn523bG3h0w5+c11pT0Y2oW1jhi8jtIvLfRORpEfmqiPz0gPe8SUQ2ROTJzr9fsqXrEtFCOvbhjW7XyhAJPa0kT5TKEa1kHoSrHLFtb9BLO1oYK36vRmOeo1SO1JptmpYVUeC3OAjX1kFjTdl4kBreaxIhw1HQ+OYm8LPGmMc7d9s+JiKfNcb8bd/7/toY83YFes7R23MkaCxNNbwRQvFqHdSYKuYoL4e3dk/MlKzoRimP7IY3rEsFo5VHJqWENJJ5YRPGaiHD0m5yvpQPtk40xjpajqhzX7KCfEXxam7pkI4x5pox5vHOz1vA08AZ2+9NElEqR2x7lUM04dSyRqK43FoHciCZMFZYT04rmRclOb9dbaok88ImjLUMiqjzbEs3SnJeq71BpDWlYDiOgmoMX0TuBL4VeGTAn79LRL4kIp8UkW8Z8R0PicgFEbmwtLSk+XiB0e0tEyK0orEzRxFOtfBGhCoKrWQeRAtv2C7IsFUUmsk8CMuzfXgDwocZNA4UQlSe7RV+lOS8ptcctq2DRkXUKKgpfBGZBv4Y+BljzGbfnx8Hzhpj7gX+I/Cnw77HGPMhY8x5Y8z5xcVFrccLhajVMnrCGX+CaaqYY6feoh2wckSjV7lPF8JbfrbtDXzaSY01xG/t+rTD0bWv/+/9fNjkvMYmFzY5r1n5BuHaOrwsqnREJI+n7D9ijPmT/r8bYzaNMdudnx8G8iJyTIO2C0QRTo0j9xA+zKCZtIXg+QON9gawW3IXNmmbxFhrhpJ6vy8QbYWKKJ92qOogdZ7DKF6deQ7r1Whv7EFpa1VEjYJGlY4AvwM8bYz5tSHvOdl5HyJyf4fuii1tV4iaSNToYR1FOG3bG/h0Ifgm57/PNrHVbesQUglpjnXQyhGtXuVRw1i2Y+3TDre5duZZ4RAShOdZbZ5DFULoVYFBcJ79iiiXFr7GN78B+CHgKyLyZOe1XwDuADDGfBD4h8BPiEgTqADvMmGv3YkRXeEMKCRayTyPdnir07a9gUc33ILU7NsdxavRsnabISpHuh0r1ZLz4cIMthVRHu1wiURtCz/UPFftCyF82uHkq0FGoJS3M6LCnjHRCiWNgvU3G2P+BhipbYwxHwA+YEsrLoRdkFrJPP87kkjmheVZs2/3VDFcMzFNhe9/XyCFr8RzKe8l58M0MdMKb4Rt3KaZwOz9vnHQNKKmizlubFYDv98fa/uKqHz3+4LSBXe98CE9aTsQYV0xLSvI+47wSVuNMq7wPOu4vRC+QkirVjnqPNvSjpKc1wxvhGnr0FVCFu0NfLre9wUca0UjKkplkg7dcBa+VkXUKKQKfwD8Az1BhUTjcgofoUM6isk8CCOcnnWqsdlESZ7qeDUhXW7FXuVheDbGqHs1QcOV27WGdXsDCH/ngtbhOogWxtL0mgPrEUXDcRhShT8AYXuOaO7MYasotKyR8MKpaOEXs6GTalqbKwR3uTV7lYexOrUqony6EG6eNfgNe+eCrtccrq2D1uYa9uS8Js/DkCr8AQjbc2TX7dVTBEGFU6Nvt0/X/76gdHs/Z4PpUj50GEsrvOF/X1C6Wr3Kw3hy2hsNhONZSwFNF/MhckQ6FVEQ/uS8toUf9ACn68tPIFX4QxFlQWq53I1WGOHUS2xBiES1Ks/BXW7tZB6ECTPo9SoPk5zXHmsIN88aBoVPO6kqMAg5zwobTdiT86mFnyA8JRRO+SURWtmqNlSsoFI+Q0YixBstk3ned4RQfsrJPAjHs5b1FSY5rxkyDHuKXJfn4POsXQUG4XjWoBu2rUN3npU22EFIFf4QhBFOjV7lvXQhWFzZb2+gkdgSkVBejUavch9h2jqoWrshczWa4Y1wY61Xrhe2PNIrCtCpGgkrX6BlRIVPzmttcmHaOuzOc1qlEzuihHQ0TkKGEU7NZJ5HO4wFppPM8+lCsOSWdgms953BLW3VsU4gmRfWg9Rs5pUUz2GT89obexjDsZjTMaKGIVX4QxA2xqqZzINgyk87yRPmCLr2ovC+c/yC1Ezm5bIZirngbR00e5VHCW+8/MNYwds66FaBBee51mzRaOm1Nwizplw3ToNU4Q9FqJ1ZMZkXxuXWTGz5tMPkLTTj2RCMZ632Bj7C1MNr9irfTc6PH28XOaIk5jlMcl6rvYFHN/ia2pUvvXkOE8bSCM+OQqrwhyBM0lZzZw7jcmsf1AhTRbGlHN6AcDxrLYywG7t/XN6aboguoZrzHCY532obdhRaYPsIlZzvhAw1jaggtLsnyEtK8xxiTWn1DhqFVOEPQRjh1Lj5qUs3inAqLsigfVa0wxv+dwahC8phrFDhDb0EJoTc2BWszjDJec2KKNhNzgdp66BqRIVIzmu3NwgXxkpDOokhTM+RrVpDJWELvcIZxPLz2htoWbth3U9tC38rhBLS9GqCHIzR7lUe5lCOZkWUTztUyDCh5LyeERUiZOiA5yQqooYhVfhDEEY4t6tNZhRdQAhm+fnKQmuzCZNg2qrq9Gf36UIyPAetHKk0WjTbRm2eu/2aAtDeUpQvCF6Q4EK+IPg8a9EN09Zhq+oZUWrzXAx+cl57ngchVfhDEKaBmrZwFnPBhHNTW/mVgoc3NHkOE8PfrDYo5DIUc/G63K6UXxDrT3OsfdrB6PrKT0++IOiaaqhvckG8Zhfz3AzY1kF7ngdB64rDB0Tk6yLyrIi8f8DfRUT+Q+fvXxaR12vQdYnuggzgcm86sMDCLMhZRWskSOVItdGi3mqr0oVgYaytapNZxUURdqyT2uSSsPB3DQqteQ7e1sHFJhfOwted53EyZoxR3+QGQeOKwyzwG8CDwD3Au0Xknr63PQi8svPvIeA3bem6RpgDUFvVhqoSCi6cTfJZnfp/CF45om0Fhakc0XZ7g461r/y0Nrmw4Q1d+QrW1mGry7NeUQAE3eSU5zlgEUZS81xrtmm0zMvCwr8feNYYc9EYUwf+EHhn33veCfy+8fAFYF5ETinQdoZd4Ry9MOrNNrVmOzGXe6aUVyld8+nCeOHUtoLCVI54POuOdZC2Duq5gxDJeRc8h/NqdJVfEkZUcE9O14gKajhudr31W1/hnwEu9/x+pfNa2PcAICIPicgFEbmwtLSk8HjREFQ4tRcFeBdGJxHbDep+dpWfUk26TzvufIlHN9j9xfrKL1xyXn2sAyaLIf4wlhsjKtidC+6MqKBe8y0e0mHwfbb95lKQ93gvGvMhY8x5Y8z5xcVF64eLiqDCqb0oIJzLrW35QVI8B1VCDVXlF35Bxp+cTzKenc0Ik0qnTsN7kNqhu8O1pgZBQ+FfAW7v+f024GqE99xSCNrTxp8ozQMT06V88PCGovLzhW1cPbwTr6aUC1STrr0gfR7GXSiuHcbyaY8b60arTaXRUh/rRstQbYzf5LROu/p0IYh86Su/mVI+oHzphs9mE1xTg6Ch8L8IvFJE7hKRAvAu4BN97/kE8MOdap3vBDaMMdcUaDtDVzjHCImLiZot5brfO5q2rvKbnfB42KyMU376C3K2lO8mzMbR1h5rgI3KeCUkotfDB2B2Ijd2rLcdjTXsxo2HQVu+SvkshVyGzQBjDcpraiI3ll+ftqYR1R3rBNbUIFgrfGNME3gf8GngaeBjxpivish7ReS9nbc9DFwEngV+C/intnRdoyuc4xaFYmtkH7MTeTYr4w9r6Cu/YMK5qVwOCh7PW2PottreZd5ONrkAym+6mCOT0bF2wRu/jSQ218Abu36ZoLexx+9NzZby1JvtQF5NMvKlz/MgqHy7MeZhPKXe+9oHe342wE9q0IoTs6V8YGtEVfmV8tRbXuKqlB8eP91Udj/97xpnaXfDWMqu77hFoXnvwC7d4Juc5hyDpwzGKfxNRx4kjPdqNh0cBAri1WgfKPTo7ireUWtKe5Mr5jIUssl4NYOQnrQdgSDC6cQamego3hG02w6s3VI+SzGXCeR+ThWyZDWt3QBejfZBM49u8E1OXfmVcmO9Gu1aeAjn1WiXCQYJ3TmZZ9+YCaB4NedZRAKFk/wxSZunJYhg7qcLa9cT9FHWX7nexBh9F3B2IpjL7cLVr7faVBvDj6C7yh1A0PBGMmMN2hZ+giGdAF6NyzDWKNrttmG77maTG89zg+liTtWIGoRU4Y/A3EQ+0KKYyGfJK15LFsQCc+UCzpZysVtB0Gtpx8vzbiIxSALTQTx7rFfjQvm9PLwaTWs3SKJ6u2tE6c7zTCA94r6PDqQKfyQ8CyyZRQGj3U9XWf1AVqdiO+gu3QBWp6vEVlBPzsUm5+dqhtN1k8CE0WNtjH7IEIJ7NZOFrOrdrnMBwqTO1lQpFyiMlSr8hOFZuwGqNxwsChhn7bqp2/WszgSs3VBejYtEYvwLMtgm59CrGTHW/kUlSXk1zsZ6hOJ1tqYCVKDF0RoZUoU/En68cZRwancyhICKwEHFCnTCWLesV+NmQc6NsTq9Tob6C3IuyCZXa3qVHkq9XXppJ+NBBvBqam5yB5CMhT9OvnzaqYWfMOYm8jTbhsqI2t3tmn6SJ0h5pIvqDQhemZSEhe+iXA/GezXVRrtz+YkbT25UeaSrlrnjymBdepAwXvFqj3UQr8ZpyDBABVpq4SeMXeEcrXhdCOe48kinIZ3qOK/GTSUDjFcEhWxmZB11JNpjvBp3Yz0+Ue1irMEvg01gcw24sbvZ5IJ6NdrGTJBcTWrhJ44glSMblUbXNdfEODfQF1wXC9LrszJYOGvNFnXlToYQ1KtpqOdLYHyupqv8lGukg4YZXCiCcfXwu7XwbkJ3Y70aB/Xo4+rhXXqQMHyed0OGqcJPFEEmamOnwdxEQZ/2mBjrRqVBIZthQtvaHXMGwH99blKX5yBezUalwbyDzdWvHBnm1ezyHH94Y6PSUB9rGG/hd3lWlu1AFn6loT7WMD50t9nlOd5wZbXRpt5qM+9Aj/QjVfgjMO6wRqVz1d+8E+EcbY1sVOrMTer17e7SHePVbOx4r7tUvMOw4VARjPJqNip1AI4oK94gXs3GTt3NWI/xanyZ15btIEbU+o7LjX34WK/v1JnIZ/VDhmO8mvWOfLnQI/1IFf4IjKui2HBkEcD4E4nOFsWYBbnuSBF4tEeXR7pTBKM3uXVHm1wQr2a90nAz1mO8Gp9nfWt39CZXrrdoto0z+RpVHrm+426sIX75GoRU4Y/AuFJBlxM1zv1MXjgdhRlGWPjrlTrzLsIb4za5HYeb3Aie223jLow1xqtZ3/GO+mueIPfpwqix7li7iciXm3xcUPly4b32I1X4IzCT4ER5CaYR1m7FUe5g7Cbnzv0Mssm58qZgxCZXaSDippPhKK9mq+od9XcTwx/j1VTqTsZ6XHmk0zU1pjxyw5kRNdqr8UOGaQw/YRRyXlJ0eEjHoTXSUX7DhbOeiIXvKoHp0x62KJqtNlvVpjNXH4Zvchs7dWZLeSeNrUZZnd3YbgJWpyvl59MeOtYVh17zmPLI9Urd2VqGZDzIfqQKfwzmRsTS3Vr4ow99rTty9f1Eop+c3Ud3x7vn1EnZXCk3dKz9jcBVMg+GJ+ddxdFhdCdF16EkSIjn5PfpFgAAE/xJREFUEYf7dnl2p3hHjbcLnsc16HOZF+uHlcIXkX8jIl8TkS+LyMdFZH7I+14Qka+IyJMicsGGZtyYnRiuhFxaI747vT5A8daaLXbqLScCUsxlmSxku0LYD9/V164OAk/gNyoN2u39Xs1uKElfEcx3x7o+8O+uksXg8TxojsGtIpgfIV/e626sXZ+2773so+uwYsX/zkE8G2O8MKkjpTs/MWKed9yUWA+CrYX/WeA1xpjXAc8APz/ivW82xtxnjDlvSTNWzE8WWBuxIPNZYbKgP1ELU95iWy3vXxiuauF9HJkssDaALrhVfgtTRVptM/Ae4XWHoaT5yQIisDpinpMZa+91F7marnwN2eRclcB6tIuslsd4zQ5kbGFy+JqqNtrUm+5q4RemCiPG2k2J9SBYKXxjzGc6d9oCfAG4zf6Rbi0cnRqt/OYmCk4maqTCd1zGdXR6lHC6VATe9w6i7ZLnbEaYn8gPnWdXtfDgyddWrUl9QFzZVS087MrXIJ5d1sJ7tEeMdaVBKa/fPgNgYbrD8wD5cl0LvzBGj8RRkgm6MfwfBT455G8G+IyIPCYiD436EhF5SEQuiMiFpaUlxceLhiNThYFKF7yd2ZWA+Id8Bgun25hfUhb+ka4FVttPt+IupAOj59llPPtIR/EOCie5tHYnC15cedDm6rIWHnbHelBBgstQkm/hrwyYZ9e18CPly2GCvB9jFb6I/FcReWrAv3f2vOcXgSbwkSFf8wZjzOuBB4GfFJE3DqNnjPmQMea8Meb84uJiSHb0cXSqwNpOfUhc2U2ZoE8XBlv4LmvhwbNGBi0KcFcLD3B0qggw0N13vSAXJgcvSJe18LBraQ9TQi5q4cG7a3VhssDq9iC6bssEj04VqLfalOv7CxJcKj9fbgcZM65r4Rcmh3vNrkqsB2FsqYUx5ntG/V1E3gO8HXiLGVJDaIy52vn/poh8HLgf+Fz4x40fRyYLtI3navrWmI/1nQYn50pO6M5O5MnIMIWfrPvpapM70gnpjFqQsw4V74urO/ted1kL79OFITw7qoXvpT3Qg3Ss/I70KN7+awxdHX4Cr8x6ppQbkhdzu8ktTBXYqDRottr7bvLa2KnzLadnndDth22VzgPAzwHvMMbsXy3ee6ZEZMb/Gfhe4CkbunHi6PTw5NZqud61xLWRzQjzQ6xOf5G6i6UXKNdbVPtKQhsOa+F9ujDY2l3bqTNbcnfJ8zCvZrVr7brleZB8rZXdhQx92qPkKwmvJimefa/SpWwbw77qN2MMqw5zRP2w9RU/AMwAn+2UXH4QQEROi8jDnfecAP5GRL4EPAr8hTHmU5Z0Y8ORIZl9Ywwr5RpHp4vOaA8TzpXtumetOKiFh12e+8vI/Gc55ojnyUKOUj4z0Opc2a5zbMbdWB/peDX9TurKtpdPcEX7yIgww0q57mysocPzgMqklU6YxxXPo7wa5zxPDvZq/Hn2DTx1ukN43qm3qDbaTmW7F1YawxjziiGvXwXe1vn5InCvDZ0kMaxaZrPSpNEyHHMkIDA8rry0XWNxuuisjMuvllkp1/aErJa2OsrP5SY3gmfXdJttw1at2T2gA7DsK3xH8+xblIOs3eWtGq88PuOELsDCZL6r6PbQ3XY7z8PWVLPVZm3HrcJfmCpwY7O67/Xl7RqzpRzFnJta+GEloa7Huh/pSdsxGCacy+UYlN8IC9+VJeLR9Xha60uernQtfHe0h1UzrGzX3G6u/jxv9y9It15NPpthbkBJqDGG5XLdMc9FNqtNGq29JaHL23XyWXFy0xbsWrv987y6U8cYt/I1bE0tO/YshuqRjny5XM+9SBX+GAybqJUYJurIkKTaStmxtTukHn7X7Y1/k1vejmlB7uO5vufvrmj3H/ra7tTmu93YB588XdmucXTKnQc5U8yRz8q+sV7ecru5wq58DQrduR3rYfLVMRynUgv/lkAp77UaGOaKHXU4UQtTedZ29rcaWN5ylyyG4XFl1+ENGFw5Um+22ag0nI71sBjr8naN+cm8k9LILu3J/Ra+a88CenjuV7zbNY7NuJtjERl41mOl7N6gODJZoNZs7+tR5dqg8EN3Q+fZ4Xj3IlX4ATDI6lyJRfl5rQZ6e/nEkSyenyyQEfbFd/1kcX8pnSYWpgpdq9pHN1nscFH4G2g/7ZVyzenmCt48L+8b63i8KWA/7XLd6ebq017eFz5zv6aGzrNjC7+UzzJdzO3j2Z9nlx5kL1KFHwCLM8VuwtLHcgyu/vFO5v7G1m6SyYu5uk0WZzPCsekiNzb38uw6WQxwYrbEdq1JubbbTyeOxNaiP9Z9Cb3lLbeWH8CJ2SI398mXe+V3YtZLyN/sm+flLbchQ5/2za39Yw3uqoMAjs/un+dGq83aTsM5z8dni/t5dpws7keq8APg5GyJaxuVPa8tb9c4Mpnfd4hClW6nQub6xq6QxJXVPzlX4nqf8nOdLAZvrIE9tONQfqV8liOT+X08LzvOl4DH82q5vufcQxwhnUFjHUey2KfdK9fgjXUh667cGHrWVA/PfpjFpTcFw3h2b1D0IlX4AXByrrTP2r2xWeP4jJtTtl26nQXZa4341tii47rdEwOE8+aWZ+G7pgt7Nznf+l2cdjveg3he2qy5H+u5/Zb2za0aIm49yKlijplibg/PmxUvWRwHz0vbtT0VQv5Yu/QgT46UL/cKf5/XvFmLrQYfUoUfCKfmvDDDVs/NRNc2Kpyad6+AAK5v7AqJ72mcctTSwcfJ2f0Wfhw8D/Jqrq17P5+Yi9er2ao22Ko1Yxlr2Gt1XluvcHym6DRZDJ7i7R3rq135mnBK9+RsCWPYEyq9ulFxPtZzE3mKucxentc9nk+7Xs9zJW5sVvcUYVzdqHDaMc+9SBV+AJwYYGlf26g6XxSFXIajU4W9iqAjqM4X5FyJjUqjG2bYqTdZ32nEogigT/ltVDg2XXQe5zzVWZC7dDtjPe9+rKGfZ/fyBfs3dl/5ud/Yvc17H8+Ox1pEONW3sce2pmZLNNume56l1Tbc2HTPcy9ShR8AviD4glFttFgt12PZmb0ww27+4KX1CgtTBSYcXLrSTxd2Le2rHSv7jGPhnChkmS3l9ijel9Yrzq0v8Hhe3q53e9O/1FF+Z2Ly5G70WZ2ux9qn3TvWVzfimed+ntttw7X1amxr6kbfJucbV67pwq7huLxdo9EynE4V/q2F/rjf9ZgsP6BjjfSEdNbdu72w39KOK5QEndDKRr+1Gx/PfiWFH0pybfnNlnJM5LPdsTbGxBLeAM/SvrlVo9UJM1xbr5DrVGk5pdsnXyvlOvVWOz756tvkTs2VyDhqzNdLF3qNqE4oKQ3p3Fo4MVdEZNfiuxqj8js1X+KltZ3uycC4XH3for6y1uG5G+eMg/ZEl64xprPJxUMXdnm+tlEhI7vlsa4gIpyeL3FlzWs4u77ToNpox2JQnJ6foNU2XQV4db3CidmSs66kPhamCpTymT1jDfEYUafnJ7i+UaXZSRjHZUTtrilvnuMKJfUiVfgBUMxlOTM/wfPLZQBeWPYm7OzRSee07zo2zWa1yUrZu4Tl0soOty+4F5DbFybJZYSLS9sAXFrZIZcRZ/3/e3HXsSmeXy7T7sQ7y/UWty/EMdZTAFxc8ub50soOp+cnnJbe+rh7cXqXbqcv/+1H3M/zLs/bXdpxyJeIcOfRqT3yBXD7kXjmudEyXO5sNpdWd2KhuzhdZKaY4+LyrnwBsYy3j1ThB8Tdi9M81xHO55a2KeUznI5hZz63uKuErm1WqTRavOL4tHO6+WyGOxYmu5vcc0vbnD066bxqBLyxrjRa3Niq8txNb8zj4PnM/ATFXIbnl3fnOQ66AHcvTnFpZYdW28TK87lFj8bzy2WM8WjHxfO5xek98iXijYN7uh6N55e32ag0WNqqxcKziHD34lR3Y39uaZsTs0VmSvH0wodU4QfGuc5EGWO4uLTNXcemncf8PLqeIF5c2u4qAv8119grnOXY6J7rsbSf69A/F4MiyGSEu455PLfbhoux8jxNvdXmytoOzy1tk89KLF7N8ZkiU4UsF5fKLG/X2aw2Y5Wvy2sV6s02zy2Vue3IhJPLy/fRPeavqXLXw4iP5+kuzeeWtmOj6yNV+AFxbnGanXqLl9YrPHMjPivodMfq/MbNbZ65sdV9ljhwbnGa51fKbNeaXFopcy4uy69D55kbWzxzY4uJfDYWbwo8np+5ucWVtQqVRiu+sT7ubWjP3PDm+ezRqVi8KRHh3PHp7lhDvPLVahsuLm/zzPUtXhET3SNTBRamCnt5js2rmeLqRpWNSoNnb7zMFL6I/CsRealz29WTIvK2Ie97QES+LiLPisj7bWgmhXtvmwfgk1+5zkvrFe67fT4WutmMcO9t81y4tMYTL65zZn7C+SlIH996xzz1ZpuPPvIijZaJjefjM0VOzZV47NIaj7+4xr23z8XiTQHcd/s8l1crfPqr17u/x4F7Ts2RzwoXLq3yxIvrsdEF+Nbb53ny8jqPPr8K7Mq6a/g8/vevL/HMzS3uu/1ILHR92o9dWuPxS+vMT+Y5G4M35dH1ePzYFy+zVWvGOs+gY+H/e2PMfZ1/D/f/UUSywG8ADwL3AO8WkXsU6MaKe07PMlPM8Suf+hoA58/GJ5z337XAly6v8+mvXuf8nfHR/fY7FwD41U97PH9bTDyLCPfftcAnn7rOl69scP7sQix0wRtr8HieKeZ41Ul3N071YqKQ5XW3zfPhv3melXI9Zvk6yk69xf/z35/lVSdmnN2V3I+zRyc5PlPk333m6xhDrLJ9/10LPLdU5k+euMK33XEkNoPi9WfnyWWku6bi5BniCencDzxrjLlojKkDfwi8Mwa6qshmhAdec5Jm23DHwiSvu20uNtpve+0pAJptw99/3enY6B6dLvJddx+l0TJ897mjsTZ5evvrTndrw7/vdadio/uaM3OcPTpJo2V48LUnnZcn9uL7XnuKRsvj+a33nIiN7hu/6RiThSyNluHtMY61iPD2153udH8tdjfbOPDga04CeDzfGx/Pk4Ucb371cRotw2vPzHH2qPvcVC802tK9T0R+GLgA/KwxZq3v72eAyz2/XwG+Y9iXichDwEMAd9xxh8Lj6eFfPPBqZifyvP11p5w2eOrHPadn+ZX/5bWslhu85ZuPx0YX4Je//zX8/ucv8Z7vvjNWum959XF+/kFvvL/51GxsdLMZ4df/8X38+Zeu8U/ffC42ugA/+J13sLZT51tOzzrv3NiLmVKe3/jB1/PIxVV+7H+6Kza6AP/sLa8gm4G3fPOJWHIWPs4eneLX//F9XFrZ4R33nomNLsC//Pv3cGZ+gnfdf3usdAGk/6qvfW8Q+a/AyQF/+kXgC8AyYID/CzhljPnRvs//I+B/Nsb8753ffwi43xjzU+Me7vz58+bChQtB+EiRIkWKFICIPGaMOT/ob2MtfGPM9wQk8lvAfxnwpytA71Z2G3A1yHemSJEiRQo92Fbp9Aa/vh94asDbvgi8UkTuEpEC8C7gEzZ0U6RIkSJFeNjG8H9VRO7DC+m8APw4gIicBn7bGPM2Y0xTRN4HfBrIAh82xnzVkm6KFClSpAgJK4VvjPmhIa9fBd7W8/vDwL6SzRQpUqRIER/Sk7YpUqRIcUiQKvwUKVKkOCRIFX6KFClSHBKkCj9FihQpDgnGHrxKEiKyBFyK+PFjeIfCDhNSng8+Dhu/kPIcFmeNMYuD/nBLK3wbiMiFYafNDipSng8+Dhu/kPKsiTSkkyJFihSHBKnCT5EiRYpDgoOs8D+U9AMkgJTng4/Dxi+kPKvhwMbwU6RIkSLFXhxkCz9FihQpUvQgVfgpUqRIcUhw4BT+QbgwPQxE5HYR+W8i8rSIfFVEfjrpZ4oLIpIVkSdEZNA9DAcOIjIvIn8kIl/rzPd3Jf1MriEi/7wj10+JyEdFpJT0M2lDRD4sIjdF5Kme1xZE5LMi8o3O/yqX3x4ohX9QLkwPiSbe1ZLfDHwn8JOHgGcfPw08nfRDxIj/G/iUMebVwL0ccN5F5Azwz4DzxpjX4LVXf1eyT+UEvws80Pfa+4G/NMa8EvjLzu/WOFAKnwNyYXoYGGOuGWMe7/y8hacE4r2kMwGIyG3A9wG/nfSzxAERmQXeCPwOgDGmboxZT/apYkEOmBCRHDDJAbwtzxjzOWC17+V3Ar/X+fn3gH+gQeugKfxBF6YfeOXnQ0TuBL4VeCTZJ4kFvw78C6Cd9IPEhLuBJeA/dcJYvy0iU0k/lEsYY14C/i3wInAN2DDGfCbZp4oNJ4wx18Az6oDjGl960BS+DHjtUNSdisg08MfAzxhjNpN+HpcQkbcDN40xjyX9LDEiB7we+E1jzLcCZZTc/FsVnbj1O4G7gNPAlIj8b8k+1csbB03hH8oL00Ukj6fsP2KM+ZOknycGvAF4h4i8gBe2+3si8gfJPpJzXAGuGGN87+2P8DaAg4zvAZ43xiwZYxrAnwDfnfAzxYUb/p3hnf9vanzpQVP4h+7CdBERvLju08aYX0v6eeKAMebnjTG3GWPuxJvjvzLGHGjLzxhzHbgsIq/qvPQW4G8TfKQ48CLwnSIy2ZHzt3DAE9U9+ATwns7P7wH+TONLbS8xv6VwSC9MfwPwQ8BXROTJzmu/0LlHOMXBwk8BH+kYMxeBH0n4eZzCGPOIiPwR8DheNdoTHMA2CyLyUeBNwDERuQL8S+BfAx8TkR/D2/j+kQqttLVCihQpUhwOHLSQTooUKVKkGIJU4adIkSLFIUGq8FOkSJHikCBV+ClSpEhxSJAq/BQpUqQ4JEgVfooUKVIcEqQKP0WKFCkOCf5/IqRQGmb7NFUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Tm = 5\n",
    "Ta = 10\n",
    "lr = 6\n",
    "T_ele = 700\n",
    "\n",
    "one_day = 60*60*24\n",
    "one_year = 365 * one_day\n",
    "timestep = 0.5 * 60*60 # 0.5 hours\n",
    "tend = 10 * one_year\n",
    "period = one_year\n",
    "t = np.arange(0, tend, timestep)\n",
    "\n",
    "Ts = Tm + Ta * np.sin(2*np.pi * t/period)\n",
    "plt.plot(t/one_year,Ts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.])"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dz = dem - T_ele\n",
    "Ts_ini[:] = Tm - dz * lr/(1000/dx)\n",
    "imshow_grid(rg, 'temperature_surface_ini', colorbar_label='Ts_ini (°C)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Frost-cracking model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we set up the model. This is based on the lab [Simulation of ground temperature profile evolution](https://csdms.rc.colorado.edu/user/jhirschberg/notebooks/learning-with-1d-models/ground_temperature_profile.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "class HeatConductionSimulator(object):\n",
    "\n",
    "    def __init__(self,\n",
    "                 surface_temp,\n",
    "                 time,\n",
    "                 thermal_diffusivity=1.0e-6,\n",
    "                 profile_depth=10.0,\n",
    "                 geotherm=0.025,\n",
    "                 num_nodes=101,\n",
    "                 diffusion_number=0.2,\n",
    "                ):\n",
    "        \"\"\"Initialize the HeatConductionSimulator.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        surface_temp: float\n",
    "            Time series of temperature, degrees C\n",
    "        time: float\n",
    "            Time array, seconds\n",
    "        thermal_diffusivity: float\n",
    "            Thermal diffusivity, in m2/s (default 10^-6)\n",
    "        profile_depth: float\n",
    "            Depth of profile, m (default 1 m)\n",
    "        geotherm: float\n",
    "            Background geothermal gradient (default 0.025 deg C / m)\n",
    "        num_nodes: int\n",
    "            Number of nodes (default 101)\n",
    "        diffusion_number: float\n",
    "            Dimensionless diffusion number = D dt / dz^2 (default 0.2)\n",
    "        \"\"\"\n",
    "        self.surface_temp = surface_temp\n",
    "        self.time = time\n",
    "        self.profile_depth = profile_depth\n",
    "        self.alpha = diffusion_number  # must be <= 1/2 for numerical stability\n",
    "        self.thermal_diffusivity = thermal_diffusivity\n",
    "        \n",
    "        self.dz = profile_depth / (num_nodes - 1)\n",
    "        self.timestep_duration_should = self.alpha * self.dz * self.dz / thermal_diffusivity\n",
    "        self.timestep_duration_is = time[1]-time[0]\n",
    "        assert self.timestep_duration_should > self.timestep_duration_is\n",
    "        \n",
    "        self.depth = np.linspace(0, profile_depth, num_nodes)\n",
    "        self.temp = surface_temp[0] + geotherm * self.depth\n",
    "        self.current_time = 0.0\n",
    "        self.current_timestep = 0\n",
    "        self.max_act_layer_depth = 0.0    \n",
    "        \n",
    "    def run_one_step(self):\n",
    "        \"\"\"Advance for one time step\"\"\"\n",
    "        self.temp[0] = self.surface_temp[self.current_timestep]\n",
    "        q = -self.thermal_diffusivity * np.diff(self.temp)/self.dz    \n",
    "        self.temp[1:-1] = self.temp[1:-1] - self.timestep_duration_is*np.diff(q)/self.dz  \n",
    "        self.current_time += self.timestep_duration_is\n",
    "        self.current_timestep += 1\n",
    "        \n",
    "        #self.temp[0] = (self.surface_temp + \n",
    "        #                self.amplitude * np.sin(2.0 * np.pi * self.current_time / self.period))\n",
    "        #self.temp[1:-1] += self.alpha * (self.temp[2:] - 2 * self.temp[1:-1] + self.temp[:-2])\n",
    "        #self.current_time += self.timestep_duration\n",
    "\n",
    "    def run_all_steps(self, n):\n",
    "        \"\"\"Advance for all time steps\"\"\"\n",
    "        for i in range(len(surface_temp)):\n",
    "            self.run_one_step()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-126-3ac43a57b667>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      7\u001b[0m         \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mhcs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtemp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0mhcs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdepth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m         \u001b[0mdisplay\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclear_output\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mwait\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m         \u001b[0mdisplay\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdisplay\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgcf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/IPython/core/display.py\u001b[0m in \u001b[0;36mdisplay\u001b[0;34m(include, exclude, metadata, transient, display_id, *objs, **kwargs)\u001b[0m\n\u001b[1;32m    311\u001b[0m             \u001b[0mpublish_display_data\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmetadata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmetadata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    312\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 313\u001b[0;31m             \u001b[0mformat_dict\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmd_dict\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minclude\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minclude\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexclude\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mexclude\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    314\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mformat_dict\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    315\u001b[0m                 \u001b[0;31m# nothing to display (e.g. _ipython_display_ took over)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36mformat\u001b[0;34m(self, obj, include, exclude)\u001b[0m\n\u001b[1;32m    178\u001b[0m             \u001b[0mmd\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    179\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 180\u001b[0;31m                 \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mformatter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    181\u001b[0m             \u001b[0;32mexcept\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    182\u001b[0m                 \u001b[0;31m# FIXME: log the exception\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<decorator-gen-9>\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36mcatch_format_error\u001b[0;34m(method, self, *args, **kwargs)\u001b[0m\n\u001b[1;32m    222\u001b[0m     \u001b[0;34m\"\"\"show traceback on failed format call\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    223\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 224\u001b[0;31m         \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    225\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    226\u001b[0m         \u001b[0;31m# don't warn on NotImplementedErrors\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/IPython/core/formatters.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    339\u001b[0m                 \u001b[0;32mpass\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    340\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 341\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mprinter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    342\u001b[0m             \u001b[0;31m# Finally look for special method names\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    343\u001b[0m             \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_real_method\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/IPython/core/pylabtools.py\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(fig)\u001b[0m\n\u001b[1;32m    246\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    247\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0;34m'png'\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mformats\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 248\u001b[0;31m         \u001b[0mpng_formatter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfor_type\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mFigure\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mfig\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mprint_figure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfig\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'png'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    249\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0;34m'retina'\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mformats\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;34m'png2x'\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mformats\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    250\u001b[0m         \u001b[0mpng_formatter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfor_type\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mFigure\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mfig\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mretina_figure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfig\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/IPython/core/pylabtools.py\u001b[0m in \u001b[0;36mprint_figure\u001b[0;34m(fig, fmt, bbox_inches, **kwargs)\u001b[0m\n\u001b[1;32m    130\u001b[0m         \u001b[0mFigureCanvasBase\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfig\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    131\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 132\u001b[0;31m     \u001b[0mfig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcanvas\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_figure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbytes_io\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkw\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    133\u001b[0m     \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbytes_io\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgetvalue\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    134\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mfmt\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'svg'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/matplotlib/backend_bases.py\u001b[0m in \u001b[0;36mprint_figure\u001b[0;34m(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, **kwargs)\u001b[0m\n\u001b[1;32m   2077\u001b[0m                             print_method, dpi=dpi, orientation=orientation),\n\u001b[1;32m   2078\u001b[0m                         draw_disabled=True)\n\u001b[0;32m-> 2079\u001b[0;31m                     \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdraw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrenderer\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2080\u001b[0m                     \u001b[0mbbox_artists\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"bbox_extra_artists\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2081\u001b[0m                     bbox_inches = self.figure.get_tightbbox(renderer,\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/matplotlib/artist.py\u001b[0m in \u001b[0;36mdraw_wrapper\u001b[0;34m(artist, renderer, *args, **kwargs)\u001b[0m\n\u001b[1;32m     36\u001b[0m                 \u001b[0mrenderer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstart_filter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     37\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 38\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mdraw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0martist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrenderer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     39\u001b[0m         \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     40\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0martist\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_agg_filter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/matplotlib/figure.py\u001b[0m in \u001b[0;36mdraw\u001b[0;34m(self, renderer)\u001b[0m\n\u001b[1;32m   1734\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpatch\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdraw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrenderer\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1735\u001b[0m             mimage._draw_list_compositing_images(\n\u001b[0;32m-> 1736\u001b[0;31m                 renderer, self, artists, self.suppressComposite)\n\u001b[0m\u001b[1;32m   1737\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1738\u001b[0m             \u001b[0mrenderer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclose_group\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'figure'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/matplotlib/image.py\u001b[0m in \u001b[0;36m_draw_list_compositing_images\u001b[0;34m(renderer, parent, artists, suppress_composite)\u001b[0m\n\u001b[1;32m    135\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mnot_composite\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mhas_images\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    136\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0ma\u001b[0m \u001b[0;32min\u001b[0m \u001b[0martists\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 137\u001b[0;31m             \u001b[0ma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdraw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrenderer\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    138\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    139\u001b[0m         \u001b[0;31m# Composite any adjacent images together\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/matplotlib/artist.py\u001b[0m in \u001b[0;36mdraw_wrapper\u001b[0;34m(artist, renderer, *args, **kwargs)\u001b[0m\n\u001b[1;32m     36\u001b[0m                 \u001b[0mrenderer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstart_filter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     37\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 38\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mdraw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0martist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrenderer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     39\u001b[0m         \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     40\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0martist\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_agg_filter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/matplotlib/axes/_base.py\u001b[0m in \u001b[0;36mdraw\u001b[0;34m(self, renderer, inframe)\u001b[0m\n\u001b[1;32m   2628\u001b[0m             \u001b[0mrenderer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstop_rasterizing\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2629\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2630\u001b[0;31m         \u001b[0mmimage\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_draw_list_compositing_images\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrenderer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0martists\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2631\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2632\u001b[0m         \u001b[0mrenderer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclose_group\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'axes'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/matplotlib/image.py\u001b[0m in \u001b[0;36m_draw_list_compositing_images\u001b[0;34m(renderer, parent, artists, suppress_composite)\u001b[0m\n\u001b[1;32m    135\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mnot_composite\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mhas_images\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    136\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0ma\u001b[0m \u001b[0;32min\u001b[0m \u001b[0martists\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 137\u001b[0;31m             \u001b[0ma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdraw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrenderer\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    138\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    139\u001b[0m         \u001b[0;31m# Composite any adjacent images together\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/matplotlib/artist.py\u001b[0m in \u001b[0;36mdraw_wrapper\u001b[0;34m(artist, renderer, *args, **kwargs)\u001b[0m\n\u001b[1;32m     36\u001b[0m                 \u001b[0mrenderer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstart_filter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     37\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 38\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mdraw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0martist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrenderer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     39\u001b[0m         \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     40\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0martist\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_agg_filter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/matplotlib/axis.py\u001b[0m in \u001b[0;36mdraw\u001b[0;34m(self, renderer, *args, **kwargs)\u001b[0m\n\u001b[1;32m   1230\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1231\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mtick\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mticks_to_draw\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1232\u001b[0;31m             \u001b[0mtick\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdraw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrenderer\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1233\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1234\u001b[0m         \u001b[0;31m# scale up the axis label box to also find the neighbors, not\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/matplotlib/artist.py\u001b[0m in \u001b[0;36mdraw_wrapper\u001b[0;34m(artist, renderer, *args, **kwargs)\u001b[0m\n\u001b[1;32m     36\u001b[0m                 \u001b[0mrenderer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstart_filter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     37\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 38\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mdraw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0martist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrenderer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     39\u001b[0m         \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     40\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0martist\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_agg_filter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/matplotlib/axis.py\u001b[0m in \u001b[0;36mdraw\u001b[0;34m(self, renderer)\u001b[0m\n\u001b[1;32m    295\u001b[0m         for artist in [self.gridline, self.tick1line, self.tick2line,\n\u001b[1;32m    296\u001b[0m                        self.label1, self.label2]:\n\u001b[0;32m--> 297\u001b[0;31m             \u001b[0martist\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdraw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrenderer\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    298\u001b[0m         \u001b[0mrenderer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclose_group\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    299\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstale\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/matplotlib/artist.py\u001b[0m in \u001b[0;36mdraw_wrapper\u001b[0;34m(artist, renderer, *args, **kwargs)\u001b[0m\n\u001b[1;32m     36\u001b[0m                 \u001b[0mrenderer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstart_filter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     37\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 38\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mdraw\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0martist\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrenderer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     39\u001b[0m         \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     40\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0martist\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_agg_filter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/matplotlib/lines.py\u001b[0m in \u001b[0;36mdraw\u001b[0;34m(self, renderer)\u001b[0m\n\u001b[1;32m    864\u001b[0m                 \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    865\u001b[0m                     \u001b[0;31m# Don't scale for pixels, and don't stroke them\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 866\u001b[0;31m                     \u001b[0mmarker_trans\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmarker_trans\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscale\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mw\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    867\u001b[0m                 renderer.draw_markers(gc, marker_path, marker_trans,\n\u001b[1;32m    868\u001b[0m                                       \u001b[0msubsampled\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maffine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfrozen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.6/site-packages/matplotlib/transforms.py\u001b[0m in \u001b[0;36mscale\u001b[0;34m(self, sx, sy)\u001b[0m\n\u001b[1;32m   2000\u001b[0m             [[sx, 0.0, 0.0], [0.0, sy, 0.0], [0.0, 0.0, 1.0]], float)\n\u001b[1;32m   2001\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_mtx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscale_mtx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_mtx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2002\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minvalidate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2003\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2004\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeZRe6V3Y+e9z733u+u61q7Rvrdbaq7ttY2ODHcycHBiSMJkQFsMwBkPInGFsjPcN29iYsJnNNhwgMcMyJ5mEMEkGw8TGS7s3t7qlltRSq1u7an3Xuy/P/FFlY3Db7q36VUvP55z3VKvrqu7vLUnf96mn7vuWUEqhaZqmXZ+McQ+gaZqmbRwdeU3TtOuYjrymadp1TEde0zTtOqYjr2madh2zxj3A15qcnFTbt28f9xiapmkvKg888MCyUmrqqd53TUV++/bt3H///eMeQ9M07UVFCHHuG71Pb9domqZdx3TkNU3TrmM68pqmadcxHXlN07TrmI68pmnadWzDIy+EeJ0Q4pQQ4owQ4uc3+nyapmna39nQyAshTOA3ge8G9gP/QgixfyPPqWmapv2djb5O/iXAGaXUWQAhxJ8A3ws8+nyeZHGY8Dv//SxvmvwiiVR8Uq3yL482EDtu5+xKgCmOsmPmCFmoeEw+wZbtPR5TNf6/0WEO9kfsj79MkN/D0doedj9h4C02WJ3fg2d0cQqBDCboTpxn6PW56licTV/K5MoFOv1j7Kud53IoKRZctg82M5yYBVVhGTHBVJu+t8g5t0/UmuLqcBK6A3amZ5mtzlP0Y/rxDPPhDJWokzl1vHqGEVREbsrFiYxVwyNOoEhsZsNzNEfnaSchg3web1ijFgcUs5sxrYhKDsimalwKMhYVCCryRBIML9PuXWa2yojiGsUoYI5NZH6dyi4RjZC47XHVrlg2ciwD0iRBDYZsHS3TSBJECJmapVO0QPiULZuqHhMGsOI49J2cokopsiEyFEwMl/DDmHqmyMtp7LyGrwLySY/CTxk4OWGjwcBMSIoRBinmqMAdjmgOYxpFhShqlHmDltUiD0yiQJE0FAPXISIjMkcYVYIxjLFGJRNxjB8XOLlFaczhKRdhC6KOQxwo+lKQS0VoRahshBWPkCOLehTjRTkN5VJWLQxh49g2Ycth5AlGdkXpGIRE5EUfp0ixBgonzGkmJUEpEaVLZQZ4jkvkSwaBRepV5JYiJqKoQqw8wRplyMigHpcEhYFdSirZwLIsMtdi1AyInILMyMnMgkKFiCTEjnPsyMSPK5wU6sKmMgMqywTXY9hwiKQitVIKWVCUISoZ4WQVViiwY2jkBl5lIUyLwvZQjkfoWoSeQWFnZEZGVUaoIsZOCsxIYccm9dzArgwcLHLXo7QdEsciqtlkVkZGjDJyqjzEjHNkUuIkEi8zcHOBKy1K6ZBISeH4xL5JIgsKI0KRU+URIs5x8hI7NrEyi3phYJsGlWWRuh65dEkck8QtyY0EpRLIE0S2dj6ZCpzcxssNHGEgbJPE9YilTSFNMleQmDFVFWJWBSLNEHGJk4OTOziZwBYCaVsktkNsOVSOJHErUpGhjATKDCPNIM3xMgOZS2Rm4SMwbUEqJZHtU0mLTJakdkkuM0pRUIqM0szBKFGqACGpiQzfypmqQzx3mB0zh3nd9tc9n2kEQGzk68kLIf4Z8Dql1I+v//qHgLuUUv/qa455A/AGgK1bt95+7tw3vKb/Gzp5dcDrfvVvuWfiffzWtIX72QG3nXsFp7fsYdbsUZuYZ2H+UU5vynmo2od7peSm8IvMjE5irtSR/Z3U1BxWfQLRLCmmB1yYFJx2AvJY4UeLtAfnmOgN8Ic+VjSJyyzSayDaOeFkzqWmxYJlUxQJZrpKo79Aqxfjj0yctIGrJjDtJmYdoo5isWWybJtkKqMsBvjDFeqDkGCocBIHv2hhGw0M3yFvQq9lsuxJIqOgqEYYWQ9/OMQfFfhDAy/38VQD0/KhZhA1TLo1ycAWpEZCoUbIsI8/SvDCCi828fMAmxqmY1MGJqOGzaprEdslmZmgihEyHuCOcvxQ4SY2XuEh8TBdSV6XDANJ3zXX/lKLCJWPsJMYN8xwogovcXArD4mL6VqkvmRYswkdk1QWZCKGfISMI5y4wI4rvMTGrQJsJIZtkdZsQk8SuhaZzElFBHmETBPsqMBO1n6PU3nYwsSwJZknGfk28fp5UiNFZCFWlmAnBTKucFOJW7hIYWFJi9STxJ5N5FhkVkFmZagixEhTZFIgkxI3M/FyF4mFaZnkvk3kOmSOSWJVpFaKKmOMLEamJVZS4KYWbi6RSmJZJqUnSVyHWFrk9lqYqzJBFClmmmElFU5m4OUSWVqYpoVwJIkniWxJJqGwMnKRQp5gZQVmWmCnAje3sAsLyzARtqSwbSLHJnUgt0pykUCZYGQZZlph5yVOamGXEhsTwzLJHYfEkWTSJJUVhZmiygSV58i8WLtfmYlbSGRlYFkWhSvJpE0qTXKpyKycqkqhSjGzAjMpkLmBm1vI0sIyDZRtkds2qVy/T2ZBYaRQZoiiwExLZK6ws7UHQykMhGVSSEkibXJbUMiKwiioqgSKHDMrsHKFzARObiGViWkaVLZNakty26KwFKlVUlgFucgozYKSHCUySpWDyKnIURSo9beVKKjW35aiohQlpajIREUuFJlQZAYkQpAYglyIZ9yx797+3Xzk2z/yjH8fgBDiAaXUHU/5vg2O/PcD3/UPIv8SpdTPPNXxd9xxh3q2z3j96U89yD899/3c+1iH6fif0G5UFJsdntz3KP9VvZR9584xX/w1zfM5Tu92HGuKoNUhn+lycWvI/fYmZDdlLjpOe+VJmksNrHgzvjWHrNcoZ4ZcmSk45TXJ45xaeImJ3pPUu4pg0MIpZ7GdDkZLEE1nnG9LFi0Llca44RXa/QVqPYE3quGVk0i7jdEQRBMll9s2q5ZBUUYY6TKN/gqNXoE7cvDzNrbRxqhJ8jastCwWPYuMlDLv4o261PsR/kjgJQF+1cKUNUQDBi2LxZpNJEtSlSCyVfzhgNowxw1NgqyGQxPTsSmaJv2GxaorSa2CTAwxoz5eGOEPC9zYws9r2KqO6VtkdZN+3WbgmiRWTq5GmHEfJ4zxQ4UXG/h5gEWAdCyymmSwfnxqpmQixkgGOFGMF1V4kYGXuUjhY9mSIpAM6g6hY5DKgkTEGNkQGcW4UYUbK7zUw1EeppQo32LkO4w8i9SuSIwI8hCZxthRgRuvPUg5lYtt2uBaxL7D0JWkriIzUgoVIbIIO8px4gI3sfAKF4nEcCSFLxm5LrFjkMp87T7kIWacYicVTiLwMhu7srEsC+VJYtdm5EoyuyIzM6oyxMhSZJxjJxVuZuLmDraQCNsk9xxGrk1qCzJZkIsEo4gxkxQrLnHS9fgXFrZtUdiS2HNIbIvUVmRmiqoiRJZhpTl2WuGkFk5u4QiJkCa5axM5DqljkFkFhUhQZYyZZVhxhZMrnEziFBLTNBGOSWK7xK4ktxSZzClVgshTjDzHiguczMQtLGRlYVoWlSOJnbXoF3ZJZhaoMoYsw8oLrEzhpiZOYSGFufbVhWMTuza5JcjtitTMyM2MwsipyKhERsXaf0P21QCXIl8PcLke35JCVBSiWovveoBTIYifZYBlpXAUOAqkEthKYClj/Wau3TCxsLAMGykcbMvFkwGeXSNwGzT9Nk2/Q7s2QSdo065NEPhtPLeFY7lYxrPfWBln5F8KvEcp9V3rv34rgFLqQ091/HOJ/GMLQ/7i393Jzf/ttQw2N0j3CT57kyA+5fOy0afgiR208iPYkw3S7T0emLVYjl2mVh9j6+LjeCuTePkO3NoU5aaQJ+cE5ywH4h7t1SeZWg4JBgF+uQXLa1FOFyxMGVzyLPJyhD1aoLO8Sq1vEMQtXGMKs+4QTQouNyV9uyQvhzijRZrdkGBgUEsbOHQwA4e4Y7DQlAwcRapGmNEKtf6I2hD82CUoW1i2R9m26DZtlj2TzEzIyj7eqEfQz/BDg1pWw6aO6VmkDclSw2FkK1IjQmR9nFFIMCzwYosg95FGgOGZRA2bbs0hsitSESLyAXIU4YUVQWjgFj6O4WG6FlHdoRvYJOsxVfkIO4xw4wIvMvAzF3t9tZ8GNj3fIfEUqZFSVSOsOMIOc5wE/NTBKV2ktMgDm9BziDyLRGYURJCGyDjFiQr8xMYr14KoXJvEcxn6a3FLzRhVRFhpjIwLnFjgZzZ2aSOlSe67RK5N5BpksqAghiLEijJkUuGm5lqghY1hm6S+Q+g5ZFKRyLU4fzWcSYWXWjiFjTQscCSJZxPZNplTkBkZZRljZSlWUqyFOTOQpYNlmVSOTexKUscilevbD2WClWZrK95U4GfrK3/TpHTXv1qwBYlVUhoJqky/eryTmXi5iVQWhlwPv+OQOQaZUVCaCarIMPMcM85xMgM3N7GRCGvt+MRx11fgFYWRQPmVgJfYucDJLGwsjK8cLx0yR5DLkszKKYycwlwLcLkeZPXVIK+tjNdWw/l6iAsKUZIbilxUpEKRGJAYgkgIqmcQYqkUTrUWYFsJ5NcEWCoDU5lfyS9S2FjCxrZcXMvHc2rU3CZ1r0Ur6DBRm6C9fqutB9i1PEzDfFZdeqF8s8hv9J78fcAeIcQO4BLwPwM/sBEn2jtTx5SHWZnoUB5O+T+bu/juL/829Uf34wc/AHsTvnCTYvVqyb7+59j/Vz4yO0i9fZBo5zT33VajPxBs7d3LxKUu049tY5vchpxosbJlN8e21MjilE7/FJNXj1G70Gb63Fbm6w3yWYuz09u43NwFyTKTK0/SWlnG77eo9efYUzeIp0vOTTRZbTeIaj2aq5coVlfwBjH1wTRe1mBrnrEwKbnqNskCkyECs+phliEysjAyBysC3y2QrkkiLJQh1/5BiBJRKUQlMISFME0K2yCTBrmVUagUJ4tx4mx9xemsbdW4klHTYrluEzolKSOsqIs3jAkGiiBx8VQTy7FJG5LlusPAVaRGSFUMcfsj/GGBH1kERR3bDFA1i2HDpudLEpmTqh5WMsQdpfhhhZ+4uGUL6Ujyhk2v7jByITVTyrKLjEOcbk4zgiD1kDSxXEnSsOkHLrFTkVgJKl9FxjHucsFEZBJkDpZYu09J4NCfslmVFbGMEdkAM0pweyWN2Fx/IHIRfo3Y8xg1Jct2SWolqHyIjPvYoxInMWhn9loQnYC05jKsS0bTimUzQ5VDzGSEDCvcUKwdW9mYtiT3awxcl7xZsWSllGWEmUbIZIATKdyuTauUWJak8ANGboe0LognCgYqwsgTzCRae1DpWTQKm4n1B5XYaxE5Fnm7YmDEDMoYM02xkhAnMfBHNnVhIiybzHWInDqZo0iDgshMyY2MggQlBpRkKFIUGYpsfWX8dzHORUlhKDJRkRqKVEBsCOJnGGO7Urjrq2G7MtaDbOJUJoH6uwhLYWObHr4M8N06Na9FO5ikHUwy2Zhioj5BqzZF4E3g2cFzWgHfCDZ0JQ8ghPgfgF8FTOD3lVIf+EbHPpeVPMBDD3+GM7/+RT73j1fY9vBRJh+cxJs7zBN3XuXzxQ5ecf5PCZ7YStM5TLUr5p6tFkkvYvuVB5m4ElAr9mNOe1zemfGYayNHV9h0+Umayx71YjuyVae3ueRMwyYt+9S6F5lcjKkPmgRiDjEpWZg1ueybpGUPv3eZiZWM2sCnpqaw6h7daYurDcnIiDDiJeqrfRo9QSNpYlttqgmT5ZbDqgexGCBHKzS6McHQoJG2sO06edNkpe3S8xSJMcQKewT9iGAoqKc1HFFHNSTdtkPPN4mMCJH3cIdD/EFFLXLxVB3Ls4kaNqs1h9gpSNUIGfdwBzlBKKhlPrbhQuDQb7gMPIPESimKAXYc4o0KgtDCLVwcy6GouQwCh9ATxGYKxWBtrz0s8SMLP/ewpaQIXAa+Q+gLMjOmLEZYSYwTlnixQZB7SFNSejah7xJ6BqnMKKsRZhxhxyVOLKilDrawwbMJfY/YE8R2Ra5GGHmEFRW4SUWQSOzSwbQtct9l6DlrwVtf/ZtpiowKvMTCy9dW55VrE/s2kWuRWRm5ShB5gkyyta9AEoldSaRtkroeoWevf8wUVYYYaYYV5zipubYqFxa4ksTxiB2LzCnISaGM1vbTkwIntXALG2kYlJ5NbLtk9tp2VSpTcjOnNDJKEipSqvU4V2So9TAXIqcQBblRkYmS1FDEBkSGIDSe3sV0slJ4Xxfjr2xLrMcYiWU4OKaLK318u0Hdb9HyJ+jUJplsTDNZn6TdmMJ3dYw32ti2a56p5xp5gH/3o2/mL14xzSs/+zfMdF7Ff74rZ+9jn2H62CytziEu3LHAPcUMhy/9BTOPb6Hm7CffE/O5LT7m0hLbrzxE+/IWavYuqm0Vx7YIljPF9NLDTF/JqQ+249WmSbZmHJt2CPOE9vJpphZi6sMZfHsTahaemLVYsCrM8ArTi5dprLrU803IoEk0q3iyYzMSCWZ4hYmlLo2uTb2YQfp1ommDi22HoZWg4iVaK6vUe1CP27iyTTlhsdCx6XqKRPVx+8s0VxNqI5da2cHyXUYTNosNm1CmlPkqfr9HvV9SCwM8mpiBzbDjsFSTRFZMla/i94cEg4p65OKqtWNGbYeVwCa2U4q8izMaEAwKaqGDV9SQrkfclqzUXGInJ1UDZDzA62fURyZeEWCbHmXNol/3GfmKyIzXt44ivJGiFjt4ysN0HMKGQz9wiO2cTIww4wHuMMePBEHq4hoOle8wDFxG/tqDTlmNsKIQJyrxIxO/8LENi6LmMvSdtfhbKVUxxIoi3KjEjyV+4WCZkiLwGPgOiQuJ/MpKfi38fmLi5x7Ssii8tfOmLkRmAmWImax989dPLNzMwTYtysAhdF1iV6x9JWMlFEayHucEtf62Il276uIrYRYluVGRGIpEqPUwP73VsqkUXsVanCsDWxnYam0LRyKROLimiydrBE6DZjBBpzbNdH2WmfYsk40ZGrUZAqeJNOVz+jeovfBuqMj/9z/7I/5s8QR3/m1K77tKVi4vsufLAbW9O/jLwyXbTn6J2TMtmq2DnD884qQS7H7yPqYuT1H39rJyU87xuoW3+jhbLqzSGm7Dbk1yZbviyUBg984xe3mVdn8Sz95EuEVwtmMTVl2aSxeYWFK0k2msoEl3k+RS0yQuV6mvXKWzXNGMJ7DdJuGMzaWWJDQGyOEirZWY5tAjUBPQdFiadFkJSrJylaDbpdYtacY1XLNJ0XJY7LgMnIy86OIPetR7FfW4hmvUKFouKy2XoVeQVn2cQZdgUFIPXTxVwwxsuq2AfgCxGWJEPbxBTG1kUct9pHTJGi4rdZfIzSiKHjIa4Q4L6pFHUHkYnsOg4TIITGIrgbyHHCV4oaKWuHjKRQUOw5rL0DdIZITKhjijGC9c2wZycWA92KFvkMgYlY+QYYwTKeqxjYsLjmQYuES+RWInVMUIK0mxw5wglriFg2XbxIFL5NtEMqdUESILseMCNzbwMwdpWlS+y8izSdYvqRP5CDPJceKKIJY4SAxpE/kesWuS2AWpzMjNhFzEXxPpFEVKKdZCnYuCzChIjYrIUIQGDA3jWwbaqxR+BU4lcNbDbCkLW1lIbGxjbe/Yt+vU/Q6d2hRTjVlmmrNMN2dpNWapeR0c00E8i28oateHce7Jv+Be9T/9MP/5N9+HGzg8kJu86tgq4q5N/Ecv465Pf4lJ+9s59toBx5YvcPN9j/Oa8m6iQ7fx2b11dl74PHNHA15hH2C4bwv33raX9uIJNp9/kukHb2LTXIPHdm7nTGMzUwuPMntpleaZXRyeFjy+rcnVTT6ZcxZ1+SKtoaJzaRLTKDnT7DDoVNjZJWS2ghW51PqSll8SBT6Z55PLcG1vtFCYpcAUisowKCuBqAqsAsxq7ZLC3LOIHEFmpFjZ2raJH7u4VQ2j4dBrOgx8RaIGeL0VGv2S5qiGa7YoWjZXOg5DNycruvjdPo1uSS2p4xk18rbDQstl6JZk1RJur0+7X9CIAlyjhaq7rDQ9Rn5FJIbIqE+wnFEbWtRUA9O1GU55XKqZxDJBpYs4w4j2EOqxj2O2KOou3Y7Hkl+QMEKGl3BXcyZHEq/wsd02Ud2nP2Fx1U4pqwF22MVdKZiOHbzKRXgtBnWX0YTJkkxR+QA77OF2K6YSB0c5CL+19gDSEvTsjMyMyFilIkIRUxFROQlFPSM3cjJRkBolsaEIzbVIF98inO56pL3KwFEmTuUwVUnmhYNr+AR2nYbXYaI2w2xrM/OdeWbam2k35gmcht7C0Dbcdfk37OXuFLG9yKbewzQat/PvZwxe9oX/wlTjNfy3WzP2P/J5brrycsJbX8KnayX7T97Pdzx0CLHrDu57pYG/cIrtJ1Nezc309uzi5C0Cb+k0W86fZ9d9O9myxeHEzBEeb1xh+tITTC232R5NU9tmc7G9myvWBfIri0wMCuoXJ9miBFfrbQYTCWbRw+p1CUKL+tDBtyWh6ZI7DpmtKKoUQ0lUYWOWFVahEMpEYGIYFqUDQ8ugUjlmOsIbFTiJjaPqULNYbtr03Yo8762t8PvQiJrYjk/YcVhuSCJrhBF1afUS6gNJoNrQcFlsrcU7rVaQ/QETg5J66OGZHvmEw0LdZ+SklOUidi9kaqioJT6OWSNrB1yt2cR2QlmuYA9imiNoJC6O2SZt+PSmHSI7pSgHyHABvwvTsYNt1Nf26Wccem5JzgArvoLdLZmKJV7lINwOg5pHb9LgihORGjGlWKX8SqytmNJPKKdSMiMnMQois2JowsAwKL9BrA2lqFfgVwK3MrArl04hmS1sHMPFk3Xqbmtta6M5z3xnM3OdzUw052l4k3prQ7vmXZeR/77/5Y387n0/QaMfsrB/xJ1HP89M89X85Y6Uu//6OO0tr+EvXy247ejf8J3DV3Llrjv5Yppw8MQ9vOzcXawc2ckDt8Pus/cx88gOXjo7yyN79/JYsMz2J44z8cQubknaHN8xx+VtNoKzTC4rZs7NU+40ON+cZzVPcdJV5Cig07PoBxarTpPE65MNYvy8ws4UlhIow6I0KkRVYVYmCItcCgpRYRQxTlzgZS6W5RHXJYOaIDaG+N0B9X5JI21jBwHLky4rdUFaLhP0urRWBY20hVkL6E76rNRL0nwZf7VLsweNpIX0ffoTHis1QWT0cfpd2t2CVhRgy4B00uNy02YoI4z4EsFiTGNgE6gm1BxW532u+hVZ1cUeDqktV7RCD1s2yRo+yzMOkR1T5V38QUzjqkE99ZFui2EzoDtjMnJiMjGiEFeo1Ggt2m5INrEW7NgoiYyKwfrqWn2DYEulqJfgVwZuZdMsJdOli28F1J0OE/UZNrW2sW1yO5undtJpbqbuNPU2h3Zduy4jDyD2b4YryzzaCLjj6l4+97qKO7/wOZpbX81/nQp5xWfOI3e9ik8fyLjt6IO8St3N8dvu5snuOfY8HPOd/s0cP3g3i1MX2fn4Ixy4dy+LN03w5E0+4cUzzF0dsT+d5YntE1zdCqVxiZkVg5kL06RbbJabs/SSi1iLKzQHFu26R9R0Sf2AzIsoygQzl8hMYTkCAwuMEkwopCIzFErFOEmOnUks4VHUYNWzyFWEM+xRG0KQNbG8gNUJh15QUWSr1LsDWj2boKpRTfhcabmEXowYLtLqZbRGHp5VI5rxudqwCM0h1qjLRL+iEQa4jsNwNmChbhGLPlZ0hfbVkkbo4dgTRJM+V2prXxFY4UW8KznTkYMjGuT1gNVJm9BNSawBBUMKQpQV0vUiirmE1MgIzYKBqVg1n3pLxFKKZgnBerA7pWSu9PCtGg2nzWRjjrn2NnZM7WB+cgcTza140tfB1rR/4LqN/Bv+9Tt5+0f/N5w0xZ+co3PmIdqtV/GFxiVec1+LR75jD50Tn+O1972MR152J8vnH2LPZ7eSHtrCvXeUHDh+D/vvvYNLRzZx8kCDXaceZvrRQzg3uxzbehNmdYzZqybb5RzhzgmWZxLceJHJfp3Zhklvtsmw3qPR61OFBfW4wmjZpNKmFANUVWAoA2UqSqEwigIrNzCEQ+xJYltQZRHeKCWIa0jHp9dyGXgVZdqlvZLQHtZwnQbdGY+luiJPr9JYHtHpefhmi3AmYKFlEqkVguVV2qsmjWoC1fS52vHoOxHm6ArN1Yx26GM7dYYzAU82BInq4Q26THQNmnkNo+axvLnG5SAnMgcU9KjUgKo5JGuHpEb61XB3TYPM+PrYGkrRKqFWmvjKY3PucZOq0/an2NTcxs7pveyc28vsxB4abksHW9OeB9dt5AG2ZnXOOiWjTTB3SvKFb1cc+esRZ145yezDn2fGfAV/c2fC7Z8/jbvpFr748pBDR4/yKvdWHjj8Evrnj7HtwTlquzucvPkA2dlTbDq5nZt3Nzi7ZR+qOsv8imS32+HU7AwrUyOcbInGss1U3WTRaxAGIbV4hExcvExRCIvSsqhM1l4DQ4FRpFg5SOWCIxjZgrKK8MIQL5ZI0yduSbqeoMqWaKxENEcBttdgZdJjxS/WLrVcjpkYNZFBjcWJgJ6fQniJiW5GJ2xgeQHLEwEDL0WlF2leSmiNfFy7SX+uzuV6xcjuUdCnrAaMJodcngpJjIShlbNiQc/8+mf+NUpFozQIKpe5wmW3Cmi6k8w0trBzei+7Z/cyP30TTX8SQ+gfYaBpL6TrOvI/8bZf4Dd+/mcZdWpMW3uYPv0g+eFbqD/6JVrtb+MhcYrv+NJu7vm2A+z48hd55fDlfObO27jpoS9x50N38eXbDnFeHGXbaYObbm5xfMcu7JMnmT57mNHBgIW5OWrhZSYWGky1HS62pohXz1EbZjRim6t+QOoISpHh5QJTCCphIMoSs1p7saXMBpVnOGmFXXqUgWQQGJTFgNYgo560MFo1ljoOI6tHbbFHp2cTWE260wGLtRIRXqazmNGO24h2nQuTLiOrh7+6TGfZpMEk4VTApbbJyOqRs4yyuqzUepw0Y4ZWzqqpWLW+PuD1UtEsDepVwN6iRsuaZHN7JzfPHWLfliPMdnbjSm8Mf7qapj0d13XkAXJf4I1yRAtmLvs8sCXnVaO7+Nstq7z8+C4+c3fOHcZ26j0AAB6pSURBVJ+9QPfWuzm1/Agvv2cXj95xJ/HZhzhy7z5O3HKEM8YJdj8Gu/d1uLBrK+apx9l8dhejvR0WZ7q4568yubCJ5c11+i2fRtjFHrp4TUlhO+SyolI5orAxzApDSDBNYhuqssCJUpzUwbBtVmqSVI0IBiFB5GH5Pstth9AYEKws01l18J0Wy9M1VoIEa7BAZxGa1SSj2QbLLcXIWiJnhXCyy+XpPpEZsyoLrlp/f//bqRQThUGj8thZ+NwqJphrb2ff7EH2b7uF+cl9+HYwxj89TdOeq+s+8j/7rl/mo296A77RpGZvZd/xY9z/ksO87PM2n7sr5uWfE9zzyu0c+cIp5KGbOcMxDtx3Mw/ddRj7xIPsP3Yb99+6j6vRg8yfCejdOsXCXJfgyWW2rM7waHuWcOUste4s9Tmbbr1GYvbxY4ElDFJpgYowSpPKEqiqwMpLTGwS1yIzIupJjJ/WUR2HQU2g8j71Xk6tbBPO1FipZZj9JdpLgjptujN1FuspMYtUzUUW26sMrIhlWbJo/v2rT/xKMVmYdMomu80ptrX3cNu2u9i/4y6mG1v09ommXeeu+8gDEPhkqwl+zWNidTuD849y6faD3H7faR56xU5u+9x5Hn/pLjY9dJTpLbdyYutpjtw7z8O3H8Y48TBHTh3iyzv2Y594nO1n9/LIju2s9E+z6WqLdiugO+HROL9Kuz9Lv+mTegOqNMVObSJhUJkmhVlRKYWZF1iFg3IMQhvMbIg3MrGcgNWmTab61FciGlmTbMJnoVGRGsvkrQErrWVG1lGu2gULX7O14leKmdxiLm9wwJhm1+TN3L7zpezfdhctf0p/A1PTbmA3ROTf9N5f5f3veiNbBjlGp8b2S7M8IU6ydMt+9n3hcR5+2RZu+cISj7z8EDd96QTNw/u4FH6ZIydv4Yt79xA8eJY9nT08sbVJ88QSW6bnuDzTZOLUKtODec402iRmlyBUVBM+mSxQowqJojJMRFlhWCaZqTDzDDs3UcHaqz7K1RH1yIGWS68OkTFgOLfEE5sfZcFOuST/7ok8llLM5iZzZYc7vR3cue3l3L3/tcy1duqQa5r2lG6IyAM4hcWoinAaLi4ttp0vOOWcwT+wi/33nuD4y3Zx8AsLHLt7B/u/dJZLtx7h6pnjHFk4yOm9LgdOr+DfsokrUyfYfKHDlb0zdOuPM7mcYTdqRMEy9SjGL1xyKamMDFFWmKJEYJE7glLluFGFFDUGdZOciFEnZbnzBF15L487kK6/UmCngLmixivVHAfmbueVN7+O3ZtvRRr6GZaapj19N0zkf+6Dv8EH3/kzzEYZV6brNGuz7H9iwIWJY6zeeoBDX1jm/tsD7n5A8MUjTe46mvCl23dx270XqN8+z1LjEXZcaPPQ3Cwzjywxlc3Tm7CYOh8T0CT0gYUSW5iktoEqQQmBUmCVUEmD3CoYzUZcnD/BeSfhilzbcpFKsTmzeYma466tr+If3fr9zLV2jPcTpmnadeGGiTyAyEuGxYDUaHBuwmTaa7D97E4uFQ9z4q6bueOBlPtv6nL7+SmObz3P7Y9v5cQegwMnUo7u3UH72BUmN82y3Flg+krBwvwEsdmlPqiz7AfkIkGmdVLDorIqCqFIfTh28AIrzgOcdNZeNrZWKrbnDe6q3cx3Hfo+7rzpH+FY7rg/PZqmXYduqMi/9Rd/iw+/41/RGCYMnJQLbsnVwwE7F44wcf9ZvnygxuHlGZZ5hHb9IOHl4zS37qfHcTanB1nonGf+yhzHp9vMPjbE3l4ndK/QChUXJz1UFWOXFYO25P7WAzzmj1he/wbp1kzyWmsP//SO1/OSva+75n+cmKZp14cbKvIAZV5Q5auYLRe/v0hhwKNTM8y3tnPLyRFnp09Tmz9I5/iTnLplOweP9nj44DxHTvR5eP9WZh+5ijvfYSjPM9lv0a/bdIYpzPscPfwYF917OOuYyEqxL2/yfbOv4l98208x1dwy7ruuadoN6Lr7oSFPx4fe9jOU9TrmyjlaAxtExWojYnnLHNuHE/hPPMHxQwGHLna4aJ6k1TxAvvAIzB+gunwcb/IgK6OTbKn28ODWVRK+wINBQSEE2zKL75j9Tn78NW+n4bU3/L5omqbdUD805OnIy5RGtyAyTXzZonJsNseKqUeWuTS1yOVbtnHkjMGCeRq1eQ/B8Ys8enAzB06FPLptmoMXQk7vm+GS+kvuDXIcBS9T23jjq9/OwW0vG/fd0zRN+6obMvLv+vDH+eBbf4q05lNUgvNzDqmV41bzTPdMvKOLXOisEG/dxU3HM07PJmwO5wk5SeDv5r49f80XgxGGglebe3jn//hbTNTnxn23NE3Tvs4NGXmArCpoDB3KMicxRtjdBcpScd71qQ612dyfYv7oAmemQ7zmdpqPXeHBI5Nctv4DZxyTu8o53v99v8tce+e474qmado3dMNG/j0f/jgf+vmfomtmmLliatGlTg0wKIsRQ/sip3b41N1tzJ+Ouf/AkHu9h5HK4B27f5J//vKfHvdd0DRN+5Zu2MgDxGWK7/jYlcStoLu5Qc+vQHg00kl2L5eIwWXuP1jypdrDzBSSj3/fnzM/uWfco2uapj0tN/RLEL7vl36Pqihp9hRYFcteiLX4JM0nnkRcfpyr1jmOHgl4sHY/ndLkD//5/6MDr2nai8oNvZIHSIsE5UIlc0SpmO4GmLUWuSfxFJzhP1Eh+MT3/BmTjU3jHlfTNO0ZuaFX8gDv/ejvQ1ERNl2MUmDWa5zdpLhsX+LEzDm+HMDrgrvZNr1v3KNqmqY9Yzd85AGyNMJOSrzCI5MGcrjE/GWHWJwkKBVv/ye/Oe4RNU3TnhUdeeC9v/JHVKrE6UdQFdR7gtGmCU57EbvzBrbljHtETdO0Z0VHfl2apuQWpHmMgUkvKFi1TOb97eMeTdM07VnTkV/33l/5A4qqIqlJMkeQmkMAOv7smCfTNE179nTkv0aajiCtKH0XQ619arIsHvNUmqZpz96GRV4I8UtCiJNCiIeFEP9BCNHaqHM9X97/q59CUeKEBV5mA7AyWBjzVJqmac/eRq7k/wo4qJQ6DDwGvHUDz/W8CUcRuQHNbklQVfTLxXGPpGma9qxtWOSVUv+vUqpY/+U9wOaNOtfz6YMf+7dUChKrZHNmsGoNxz2Spmnas/ZC7cn/GPBfnuodQog3CCHuF0Lcv7S09AKN881F4QByRTvvcM4uubp6YdwjaZqmPSvPKfJCiE8LIY49xe17v+aYtwMF8Kmn+hhKqY8rpe5QSt0xNTX1XMZ53nzwNz5FWRU0im2UQvDLf/yWcY+kaZr2rDyn165RSr3mm71fCPEjwD8GvlNdSz9n8GmIohGTi21ajYor5plxj6NpmvasbOTVNa8D3gJ8j1Iq2qjzbJQPfeyPKYqc3XGDU27E6SceHfdImqZpz9hG7sl/DKgDfyWEeEgI8TsbeK4NEYUjJtObSAzBr/0nvWWjadqLz0ZeXbNbKbVFKXXL+u0nN+pcG+VDv/nH1BckWzLFBe/cuMfRNE17xvQzXr+FsD9ge7KVs47i5z7yg+MeR9M07RnRkf8WPvw7f8rs0laCquJy/bFxj6NpmvaM6Mg/DcOlFQ5EHY55EW9774+NexxN07SnTUf+afil3/5z5nv7Abg6eXrM02iapj19OvJPU7ZYcCjyOeav8ra3vn7c42iapj0tOvJP0y/+2u+xZXiI2DBYntcvc6Bp2ouDjvwzUFwoOBhJjvlXeevP/vC4x9E0TfuWdOSfgY987A/ZGh5iaBp0d14d9ziapmnfko78M+Sf77A/NjnuX+JtP/sD4x5H0zTtm9KRf4be/eu/wtbwED3TYHXH6rjH0TRN+6Z05J8F71yHfbHBCf8Cb/vZHxr3OJqmad+Qjvyz8L5f/zW2RYdZtQx6O5bHPY6mado3pCP/LPln2+xLDI7753nrz79+3ONomqY9JR35Z+l9H/t1toXrq/n5hXGPo2ma9pR05J+D+pNT3JQYPOqf5+1vef24x9E0Tfs6OvLPwbt//d+wfXSIVctgdYu+bl7TtGuPjvxz9COv+Nfsi02Oexd4x8/96LjH0TRN+3t05J+jQy95CdvDQ3Qtg+Utl8c9jqZp2t+jI/88+KX/49+ureaDi7ztrT8+7nE0TdO+Skf+ebJjdHjtWbDz58c9iqZp2lfpyD9PPvKmP2J/bHHMv8Q736ZX85qmXRt05J9H20aH6ZsGi3Pnxj2KpmkaoCP/vPrIm/6QA5HFseAy73rnG8Y9jqZpmo78823b6DAD02Bh5uy4R9E0TdORf759+M1/yMFI8rB/hfe8+43jHkfTtBucjvwG2DY6zMg0uDz12LhH0TTtBqcjvwF+8c1/wKHI5hH/Ku99x/867nE0TbuB6chvkO2DI4xMg4uzZ8Y9iqZpNzAd+Q3ywbf8Pocjh4eDRd7zLr2a1zRtPHTkN9DWwWFiIbg0c3rco2iadoPa8MgLId4khFBCiMmNPte15kNv+X0ORx5H/SXe9R79LFhN0154Gxp5IcQW4LXADfuCLtuGh0mE4PKkXs1rmvbC2+iV/K8APweoDT7PNesDb/k9jkQ+R4MV3vW+nxz3OJqm3WA2LPJCiO8BLimljn6L494ghLhfCHH/0tLSRo0zVlsHh8kEXJp4dNyjaJp2g3lOkRdCfFoIcewpbt8LvB1417f6GEqpjyul7lBK3TE1NfVcxrlmfeDnP8mRKOCov8q7369f00bTtBeO9Vx+s1LqNU/1/4UQh4AdwFEhBMBm4EEhxEuUUjfkD0Pd2j/Ew/49XOjo1bymaS+cDdmuUUo9opSaVkptV0ptBy4Ct92ogQf4hbd+kiNRjaN+j3e/X183r2naC0NfJ/8C2rJ6kEqgV/Oapr1gXpDIr6/ol1+Ic13LfuEdn+RIWOehoM+736+vm9c0bePplfwLbHP3AAo4p1fzmqa9AHTkX2C/8I5PckvY4Ggw4J0f0FfaaJq2sXTkx2DL6gEEcKH1yLhH0TTtOqcjPwbve+cnuCVscjQY8s4P6L15TdM2jo78mGxe3Y+p4Hz72LhH0TTtOqYjPybve+cnOBK1OOqPeMeH9Gpe07SNoSM/RltXb8ZScL6p9+Y1TdsYOvJj9O53foIjUYejfsjbf1Gv5jVNe/7pyI/Z1uV92ArONx4e9yiapl2HdOTH7N3v/jhHwgmO+hFv+7BezWua9vzSkb8GzC/tw1OK841v+tL7mqZpz5iO/DXgve/9XQ5HUzzsxbztwz867nE0TbuO6MhfI+av7sVXinMNfaWNpmnPHx35a8R73v9xDoUzPOynvPXDPzbucTRNu07oyF9DNl/dRa2sOKf35jVNe57oyF9D3v0Ln+BQNMsjfsbP/9Lrxz2OpmnXAR35a8ympb3Uy4pzNX3dvKZpz52O/DXmPe/9bQ5Fcxzzc97y0dePexxN017kdOSvQdOXt9EsK56s6b15TdOeGx35a9D7P/h7HIzmedQreMtHf3jc42ia9iKmI3+Nmr60jVZZ8YTem9c07TnQkb9Gve9Dn+BAuJkTXsmbf/mHxj2OpmkvUjry17DOhS20i4onA/0sWE3Tnh0d+WvYBz/ySQ5EWznplbzp3/zguMfRNO1FSEf+Gjd9cYaJouKJ4BHe/tM/Ne5xNE17kdGRv8a998N/wM3RNh5zK5JdvXGPo2nai4yO/ItA5/I0k0XFE8ExvZrXNO0Z0ZF/EfjAh/6Am6MdnHYV8e7uuMfRNO1FREf+RcK/3GC6qDjrH9ereU3TnjYd+ReJj37oj9kX7eRxVxHtXh33OJqmvUjoyL+INE/5zOQVZwO9mtc07enZ0MgLIX5GCHFKCHFcCPGRjTzXjeCDH/tT9kW7OOtAuGdl3ONomvYisGGRF0K8Gvhe4LBS6gDw0Y06142kcdJnLlc87j+qV/Oapn1LG7mSfyPwi0qpFEAptbiB57phfPC3/oSbot086UC4Z2nc42iado3byMjvBV4hhPiSEOIzQog7n+ogIcQbhBD3CyHuX1rS0Xo6nEctNuWKM/5J3vwT/2zc42iadg17TpEXQnxaCHHsKW7fC1hAG7gbeDPwZ0II8Q8/hlLq40qpO5RSd0xNTT2XcW4YH/2d/4uboj2ccyDbP+5pNE27lj2nyCulXqOUOvgUt/8IXAT+vVpzL1ABk8/H0BqYJ3w2Z4rT/kne/JN6Na9p2lPbyO2a/xv4DgAhxF7ABpY38Hw3lF/57U+xJ7qJC7YgvVmNexxN065RGxn53wd2CiGOAX8C/IhSStfoeSQeNdiSKU77p/jf3/gvxz2OpmnXoA2LvFIqU0r94Pr2zW1Kqb/ZqHPdqH7td/+cPdH/3969B0d13mcc//6EQGYXgwFzsTECbMAgQPIkimNn0lzaTJu603Hd0pomTeomDWmcyTRtsQlgMLLBxDZppxcTBzQTJp6mbuPETVonTsbNNJ6O47byDBJ3EJiLwDYYzE0rLpJ+/WMPUxnrstKy+549+3xmNJzVOdrzzFnp2Xdfds+ZQ9sIo2tuJnQcEYkhfeK1xI3YAdMuwO6Re1miuXkRuYJKvsQ9+a3nmNkxh6MjjAs1naHjiEjMqOQTIL1nAtMvwO5Uq0bzIvIuKvkEWPvUBm7J1PDGcON8zcXQcUQkRlTyCZHeO56bL8Cu1D6W378odBwRiQmVfEKsfWoDN7fP463hFZya2x46jojEhEo+QVKt47jlgrE7tZ+HvnJv6DgiEgMq+QS5PJo/VlnBydkazYuISj5xRraOZdZ5Y2fqdVYsuy90HBEJTCWfMNnR/Hzerqzg5I06hb9IuVPJJ1DVvuuYfb6CnamDPLz0vtBxRCQglXwCrX1qAzPaazlRWcGxqW+GjiMiAankE2r9Xz7Drecr2D7yMMsf/NPQcUQkEJV8gs04V8s7lRW8M/Vw6CgiEohKPsGe/KtnmNsxjG3pNpYv02hepByp5BNu+rkFnBpWwckph0JHEZEAVPIJ98SSZ6jpqGRb6ggrl38+dBwRKTKVfBmYfq6O08MqOHaDRvMi5UYlXwYeX7KZ+R3D2Zo+yqqVi0PHEZEiUsmXiWlnazk7rII3J+0LHUVEikglXya+/sBm5meGszX1JqtX6p02IuVCJV9Gpp+t49ywCo5M2h86iogUiUq+jKx78NvUZqrYmn6L1ZqbFykLKvkyM+3MAtorKjgyeU/oKCJSBCr5MvPY0m9Tl6miJXWchx/+Yug4IlJgKvkyVH3mNjrMODJhV+goIlJgKvky9NjSRuoyKZrTJ2ho0Ny8SJKp5MtU9ZlaLhocul6jeZEkU8mXqbVfi0bzqZM0PPqF0HFEpEBU8mWs+nQtnQaHxu0MHUVECqRgJW9mt5nZq2a2xcyazOz2Qu1LhmbNskbqMqNoTp1ilUbzIolUyJH8E0CDu98GrIpuS8xUvzOfLoO2cTtCRxGRAihkyTswOloeAxwt4L5kiB5d0Uhd+7VsSZ/WaF4kgSoLeN9fBX5qZuvJPpl8qLeNzGwxsBigurq6gHGkL9WnFtCcfoVD47aHjiIiV5m5+9B/2OwlYHIvq1YAvwb8wt2/b2Z/ACx290/0d3/19fXe1NQ05DwydPdtuJPm1Fl+q+1O1jy0KXQcERkEM3vN3et7W5fXSL6/0jaz7wB/Ht38HtCYz76ksKadmEdL6lXaxm4LHUVErqJCzskfBT4aLf8qsLeA+5I8NaxspK59DM3pszy0RuebF0mKQs7JfwH4WzOrBM4TzbtLfFWfrKEl9UsOj9NoXiQpCjaSd/f/cvf3u3udu3/Q3V8r1L7k6mhYuYm6zHU0p87x0DqN5kWSQJ94lXeZerKGSodDY7aGjiIiV4FKXt6lYeVG6jJjaU61s+LrGs2LlDqVvLzHlBPzGO5waHRL6CgikieVvLzHI6ue5rbMOFpSGVY8/rnQcUQkDyp56dWU47OpcufgtZqbFyllKnnpVcPqRuoyE2hJdbBMo3mRkqWSlz5NeWsWI901Ny9SwlTy0qfVj2yitn0iLakLLNdoXqQkqeSlXze9OZNRXd0cGN0cOoqIDIFKXvr18JpNLMhMZmvqIl978r7QcURkkFTyMqAbj89mVFc3B0dpbl6k1KjkZUCrG75JbeYGtqUusVSjeZGSopKXnEx662ZGd3VzcJTm5kVKiUpecvLIoxuZn7mR7alOlq7/49BxRCRHKnnJ2aQj0xnT1c0Bzc2LlAyVvOTskXWbWNB+EztGdvLgNz4TOo6I5EAlL4My9vBUxnZ283pa57QRKQUqeRmUx55oZF7HVHaN7GLJN/4odBwRGYBKXgZt4uHJjO/s5sAojeZF4k4lL4PW8Phm5mamsfuabpb89adCxxGRfqjkZUjGHZ3I+M5u9qe3seLL94eOIyJ9UMnLkKxdt5l57TPYe41zfuap0HFEpA8qeRmykW+MZkJnN/tTGs2LxJVKXoZs/brvMjczg9ZrnI6ZJ0PHEZFeqOQlL2P3pJl4qZv96R0azYvEkEpe8rLm7/+ZOZmb2VflZGadCB1HRK6gkpe8jdmVZvKlbvalNJoXiRuVvOTtsQ3PMiczk9eroH3W26HjiEgPKnm5Kqp2DOeGS920pnbywBcXho4jIhGVvFwV659+jlszszlYBRdrQqcRkctU8nLVVO5MMeWi05rapdG8SEzkVfJm9vtmtt3Mus2s/op1y8ys1cx2m9lv5BdTSsHffPMfmZ25lUMjjIs1HjqOiJD/SH4b8LvAyz2/aWY1wCJgHvBJYIOZDctzX1IChu26hpsuOntSu/mLL306dByRspdXybv7Tnff3cuqu4Fn3f2Cu78OtAK357MvKQ2XR/NtI4yuuR2h44iUvULNyU8BDve43RZ97z3MbLGZNZlZ0/HjxwsUR4ppxA5jQWYEFd3DQ0cRKXuVA21gZi8Bk3tZtcLdf9jXj/XyvV4nad19I7ARoL6+XhO5CfDkt54LHUFEIgOWvLt/Ygj32wZM7XH7JuDoEO5HRETyUKjpmh8Bi8ysysxmALOA/ynQvkREpA/5voXyHjNrA+4EXjCznwK4+3bgX4AdwIvAl929K9+wIiIyOANO1/TH3Z8Hnu9j3VpgbT73LyIi+dEnXkVEEkwlLyKSYCp5EZEEU8mLiCSYucfn80dmdhw4mMddXA/E+aoVypcf5cuP8uUnzvmmufuE3lbEquTzZWZN7l4/8JZhKF9+lC8/ypefuOfri6ZrREQSTCUvIpJgSSv5jaEDDED58qN8+VG+/MQ9X68SNScvIiLvlrSRvIiI9KCSFxFJsJIreTM7YGZbzWyLmTX1st7M7O+ii4i3mNn7YpZvjpn90swumNmSYmbLMd+no+PWYmavmFldDDPeHeXbEl1V7MNxytdjuw+YWZeZLYxTPjP7mJmdjtZvMbNVccrXI+MWM9tuZr+IUz4ze6DHsdsWPcbjiplxUNy9pL6AA8D1/ay/C/gJ2atT3QH8d8zyTQQ+QPYMnUtiePw+BIyNln+z2Mcvx4yj+P//T6oFdsUpX7TNMODnwI+BhXHKB3wM+PdiP66DyHcd2dOUV0e3J8Yp3xXb/jbw81DHMpevkhvJ5+Bu4Due9SpwnZndEDrUZe5+zN3/F7gUOktv3P0Vd38nuvkq2at6xYq7n/PoLwxI08elJQP7CvB94FjoICXoU8AP3P0QZP9mAufpzx8C/xQ6RH9KseQd+JmZvWZmi3tZn/NFxAtkoHyhDSbf58m+Kiq2ATNGF6zZBbwAfK6o6QbIZ2ZTgHuAp4uc67JcHuM7zazZzH5iZvOKGY6B880GxprZf0bbfDZm+QAwsxTwSbJP5vEV+qXEEF5K3Rj9OxFoBj5yxfoXgA/3uP0fwPvjkq/HdqsJM12Ta76PAzuB8XHNGG3zEeClOOUDvgfcES1vpvjTNQPlGw2MipbvAvbGLN8/kH0VmSZ7vpi9wOy45Oux3b3AvxXz2A3lq+RG8u5+NPr3GNmrUt1+xSZBLyKeQ76gcslnZrVAI3C3u58obsLBHUN3fxm4xcyuL1K8XPLVA8+a2QFgIbDBzH4nLvnc/Yy7n4uWfwwMj9nxawNedPd2d38beBko2hsABvH7t4iYT9VAiU3XmFnazK69vAz8OrDtis1+BHw2epfNHcBpd38jRvmCySWfmVUDPwA+4+57YppxpplZtPw+YARQlCejXPK5+wx3n+7u04HngPvd/V/jks/MJvc4freT7YHYHD/gh8CvmFllNCXyQbKvKuOSDzMbA3w0yhpreV3jNYBJwPPR72cl8F13f9HM/gzA3Z8m+26Gu4BWIAP8SZzymdlkoInsS+ZuM/sqUOPuZ+KQD1gFjCc7+gTo9OKeeS+XjL9H9on8EtAB3OvR6+eY5Aspl3wLgS+ZWSfZ47coTsfP3Xea2YtAC9ANNLp7sQZLuT6+9wA/c/f2IuUaMp3WQEQkwUpqukZERAZHJS8ikmAqeRGRBFPJi4gkmEpeRCTBVPIiIgmmkhcRSbD/Awk6zwmeOS4DAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython import display\n",
    "hcs = HeatConductionSimulator(surface_temp = Ts, time = t)\n",
    "i = 0\n",
    "while i<len(Ts):\n",
    "    hcs.run_one_step()\n",
    "    if i%(8700*2/12) == 0:\n",
    "        plt.plot(hcs.temp, -hcs.depth)\n",
    "        display.clear_output(wait=True)\n",
    "        display.display(plt.gcf())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "213"
      ]
     },
     "execution_count": 128,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1450.0"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}