Archer climate class number cruncher 1

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Global Warming: The Science of Climate Change\n",
      "=============================================\n",
      "\n",
      "Number Cruncher #1: A time-dependent climate mode\n",
      "-------------------------------------------------\n",
      "\n",
      "This iPython Notebook presents the code to solve the first programming assignment in the [**Coursera** MOOC](https://class.coursera.org/globalwarming-001/class/index) taught by **David Archer (U Chicago)** in Oct. 2013.\n",
      "\n",
      "I am using numpy for the data structures and matplotlib for plotting."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "np.set_printoptions(precision=2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The goal of the exercise is to calculate the convergence to equilibrium for a simplified dynamic Earth climate model whereby the planet is covered by an ocean of depth D and starts out at an initial temperature T. In each time step, we calculate incoming and outgoing heat flux, total heat content and global temperature.\n",
      "\n",
      "First set up the basic variables. If we want to change parameters, it is sufficient to edit the next cell and then re-run those that follow."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "d_ocean = 100                    # ocean depth in meters\n",
      "L = 1350                         # solar constant on Earth surface, W/m^2\n",
      "sigma = 5.67e-8                  # Stefan-Boltzmann constant, W/(m^2*K^4)\n",
      "epsilon = 1                      # emissivity\n",
      "alpha = 0.3                      # albedo of the Earth, unitless\n",
      "c_water = 4.1855e3               # heat capacity of liquid water, J/(kg*K)  \n",
      "tstart = 0.                      # time in years: tstart, tstop, tstep\n",
      "tstep = 1.\n",
      "nstep = 2000                     # number of time steps after initialization\n",
      "tstop = tstep * nstep + tstart \n",
      "T_ini = 0.                       #initial temperature, K"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Next, we initialize some dependent variables and the arrays to hold the data of interest."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# intitialize\n",
      "c_surface = c_water * 1000 * d_ocean                   # heat capacity of the surface layer: J/(K*m^2) \n",
      "L_solar_in = L * (1-alpha) / 4.                        # solar heat flux in, W/m^2 - constant! \n",
      "tstep_sec = tstep * 365.25 * 24 * 3600                 # duration of a time step in sec\n",
      "hc_ini = T_ini * c_surface                             # initial heat content, or a proxy\n",
      "\n",
      "# set up arrays\n",
      "times_years = np.arange(tstart, tstop+tstep, tstep)    # time steps in years\n",
      "T = T_ini * np.ones(nstep+1)                           # temperature steps\n",
      "L_out = T_ini**4 * epsilon * sigma * np.ones(nstep+1)  # outgoing heat flux steps\n",
      "hc = hc_ini * np.ones(nstep+1)                         # heat content steps"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Some useful functions, as well as the (constant) solar heat that flows into the system in each time step."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def d_heatcontent(d_T):\n",
      "    ''' heat content change per square meter for given temperature change in K'''\n",
      "    return d_T * c_surface\n",
      "\n",
      "def d_T(d_hc):\n",
      "    '''temperature change for given heat content change'''\n",
      "    return d_hc / c_surface\n",
      "\n",
      "def d_heat_out(T, time):\n",
      "    ''' heat loss per unit area over time in sec for body at temperature T'''\n",
      "    return time * sigma * epsilon * T**4 \n",
      "\n",
      "d_heat_solar_in = tstep_sec * L_solar_in"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we can calculate for each time step:   \n",
      "1. the change in heat content   \n",
      "2. the new total heat content (from the change in heat content)   \n",
      "3. the new temperature (from the change in heat content)  \n",
      "4. the new outgoing heat flux  "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i in xrange(1, nstep+1):\n",
      "    dhc = d_heat_solar_in - d_heat_out(T[i-1], tstep_sec)\n",
      "    hc[i] = hc[i-1] + dhc\n",
      "    T[i] = T[i-1] + d_T(dhc)\n",
      "    L_out[i] = d_heat_out(T[i], tstep_sec) / tstep_sec"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that the incoming heat flux in each time step is constant:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "L_solar_in  "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "236.24999999999997"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we can plot whichever quantity we're interested in."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = plt.figure(1, figsize=(10, 6))\n",
      "plt.plot(times_years[:100], L_out[:100])\n",
      "plt.title('Outgoing heat flux over time')\n",
      "plt.xlabel('Time step in years')\n",
      "plt.ylabel('Outgoing heat flux in $W/m^2$')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "<matplotlib.text.Text at 0x10588f9d0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAnAAAAGHCAYAAAAwbG+fAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlYlXX+//HXQcgNEnLBBRXMBVEUknAhDTM1LZe0GJfM\nspqr/OVMU03WTBlao2bjlFZOfR1Ts3LJJnMXJ8Mxl9wwHVExlUJQXEANNxA/vz9OnsSVAxzuc+D5\nuK5zcZab+37D3YWvPqvNGHNKkp8AAADgEWzGGGN1EQAAACg8L6sLAAAAgHMIcAAAAB6GAAcAAOBh\nCHAAAAAehgAHAADgYQhwAFxizZo1Cg0Ndcm5Y2NjNW3aNJec+2bOnj2rXr16yd/fX3FxcZo5c6Y6\nduxoSS2lyc/PT6mpqVaXAeBXBDigDJkxY4bCw8NVtWpV1alTR8OHD9fJkycL/f3BwcFatWpVidTS\nsWNH7d69u0TOdSWbzSabzVbi501MTFT9+vVveMz8+fN15MgRZWVlad68eSqLKzFdKyD/8ssvCg4O\ntqYgAFchwAFlxMSJE/Xyyy9r4sSJOnXqlDZs2KCffvpJXbt2VV5eXqHOYbPZymQgKUk//fSTmjZt\nKi+vsvHn8+LFi1e954pwDKBklY2/QEA5d+rUKcXHx+v9999Xt27dVKFCBTVs2FDz5s1TamqqPv30\nU0nSY489ptdee83xfZe3OA0ZMkQ///yzevXqJT8/P/3973+XJH3yySdq2LChatSooTfffFPBwcH6\n5ptvJEnnz5/Xc889p3r16qlevXr605/+pNzc3KvOLdlb9yZOnKjWrVvL399fAwYM0Pnz5x2fT5gw\nQXXr1lVQUJD+9a9/ycvLS/v377/uz5yamqq77rpLt956q7p3767jx487PtuwYYM6dOiggIAARURE\naPXq1Y7Ppk+frrCwMN166626/fbb9X//93+SpNOnT6tHjx7KyMiQn5+fbr31Vh0+fLjANV9//XW9\n8cYbmjt3rvz8/PTxxx8XCDupqany8vIqEIoub8165pln9NBDDzk+GzlypO69995r/nzGGMfvOzAw\nUEOHDtWpU6ckST169NAHH3xQ4PjWrVtrwYIFkqTdu3era9euql69ukJDQ/XFF184jnvsscf0zDPP\nqGfPnvL19VViYmKB8/z1r3/VmjVr9Oyzz8rPz09/+MMfJKnA/Xjsscc0fPhw9ezZU35+furYsaMO\nHz6sP/7xjwoICFDz5s21bds2xzkzMjLUv39/1apVS40aNdJ77713zZ8ZgBMMAI+3bNky4+3tbfLz\n86/6bOjQoWbgwIHGGGMee+wx89prrzk++/bbb01QUJDjdXBwsPnmm28cr3fu3Gl8fX3N2rVrTW5u\nrnnxxReNj4+P45jXXnvNtG/f3hw9etQcPXrUdOjQwXH+a527bdu25tChQyYrK8s0b97cfPjhh476\na9eubZKTk82ZM2fM4MGDjZeXl9m3b981f967777b3H777Wbv3r3m7NmzJjY21rz88svGGGMOHjxo\nqlevbpYtW2aMMWblypWmevXq5tixY8YYY5YsWWL2799vjDFm9erVpkqVKmbr1q3GGGMSExML1Hwt\n8fHxZsiQIY7X06dPN3fddZcxxpgDBw4Ym81W4D7ExsaaadOmGWOMOXPmjGnatKmZMWOG+e9//2tq\n1Khh0tPTr3mdadOmmcaNG5sDBw6YnJwc069fP8d1P/nkExMTE1PgPvn7+5vc3FyTk5NjgoKCzIwZ\nM0x+fr5JSkoyNWrUMMnJycYY+38P1apVM+vWrTPGGHPu3Lmrrn15zZfYbDbH/Rg6dKipUaOG2bp1\nqzl37py55557TMOGDc2sWbPMxYsXzauvvmo6d+5sjDEmPz/f3HHHHeaNN94weXl5Zv/+/aZRo0Zm\nxYoVN/w9A7gxWuCAMuDYsWOqUaPGNbv1ateuXaB1yjjRRTp//nz17t1bHTp0kI+Pj8aMGVOgxenz\nzz/XqFGjVKNGDdWoUUOvv/66Zs2add3z/eEPf1Dt2rUVEBCgXr16OVpp5s2bp2HDhql58+aqXLmy\nRo8efcM6bTabhg0bpsaNG6tSpUqKi4tznOvTTz9Vz549dd9990mS7r33XkVFRWnJkiWSpJ49eyok\nJESS1KlTJ3Xr1k1r1qwp9O/GGFPkbubKlStr1qxZ+tOf/qQhQ4bo/fffV926da957GeffaYXXnhB\nwcHBqlq1qsaNG6c5c+bo4sWL6tu3r7Zt26a0tDTHsf3795ePj48WL16skJAQDR06VF5eXoqIiFC/\nfv0KtML17dtX7du3lyRVrFjxuj/n9dhsNvXr10+RkZGqWLGiHnzwQVWtWlWPPPKIbDab4uLilJSU\nJEnatGmTjh07pldffVXe3t4KCQnRk08+qTlz5hTpdwjAjgAHlAE1atTQsWPHrjme6dChQ6pRo0aR\nznvo0CEFBQU5XleuXFnVq1d3vM7IyFDDhg0drxs0aKCMjIzrnq927doFznX69GnHdS7vbr38moU9\nV05OjiT7GLUvvvhCAQEBjsfatWsd3aHLli1Tu3btVL16dQUEBGjp0qUFAq6rRUdHq1GjRpKkhx9+\n+LrHHTp06Krf7YULF5SZmSk/Pz/df//9mj17tiRpzpw5Gjx4sCT7z//9998X+Pk///xzZWZmSrKH\nr5tN1Lh03I3UqlXL8bxSpUoFXl95PzIyMgrUM27cOB05cuSmNQC4PgIcUAa0b99eFStW1Jdfflng\n/ZycHC1fvlxdunSRJFWtWlVnzpxxfH7lGK8r/9GuU6eODh486Hh99uzZAmGnbt26BZaW+Pnnn6/b\nonQjderUcbQmSSrw3FkNGjTQkCFDlJ2d7Xj88ssveumll3T+/Hn1799fL730ko4cOaLs7Gz17NnT\n0dpUmMH7NzqmatWqknTD3/EHH3yg3Nxc1a1bVxMmTLjuua71u/X29lZgYKAkaeDAgZo9e7bWr1+v\nc+fOqXPnzo6f/+67777q579yzFxRf0Zn1a9fXyEhIQXqOXXqlBYvXlxi1wDKIwIcUAZUq1ZNr7/+\nukaMGKEVK1YoLy9PqampiouLU/369TVkyBBJUkREhJYuXars7GwdPnxY7777boHzBAYGat++fY7X\nDz30kBYtWqT169crNzdX8fHxBbrWBg4cqDfffFPHjh3TsWPHNGbMGMe1CuPSueLi4jR9+nTt3r1b\nZ86c0RtvvFHo773SI488okWLFikhIUH5+fk6d+6cEhMTlZ6ertzcXOXm5jq6m5ctW6aEhIQCP//x\n48cdkwWcua4k1axZU/Xq1dOsWbOUn5+vjz/+uMDvMyUlRa+99po+++wzffLJJ5owYYJ++OGHa55r\n4MCBeuedd5SamqqcnBz95S9/0YABAxzd5D179tRPP/2k119/XQMGDHB83wMPPKCUlBR9+umnysvL\nU15enjZt2uRY0qUw3b9X/nfgzO/gStHR0fLz89OECRN09uxZ5efn63//+582b95c6HMAuBoBDigj\n/vznP2vs2LF68cUXVa1aNbVr104NGzbUN998Ix8fH0n2maatW7dWcHCw7rvvPg0YMKBAa8srr7yi\nN998UwEBAfrHP/6hsLAwvffeexowYIDq1q0rPz8/1apVyzFu6tVXX1VUVJRatWqlVq1aKSoqSq++\n+qrjfDdqybl8Lbf77rtPf/jDH9S5c2c1bdr0puOzrjz35ecKCgrS119/rbFjx6pWrVpq0KCBJk6c\nKGOM/Pz8NHnyZMXFxem2227T7Nmz1adPH8d5QkNDNXDgQDVq1Ei33XbbVa1nV17rWq+nTp2qt99+\nWzVq1FBycrJiYmIkSfn5+RoyZIhefvllhYeHq3Hjxho7dqyGDBlyzWVehg0bpiFDhqhTp05q1KiR\nqlSpUmD25i233KJ+/frpm2++0aBBgxzv+/r6KiEhQXPmzFG9evVUp04dvfLKK47ZwYVZQ++Pf/yj\n5s+fr9tuu03PPfec07+DS+9JUoUKFbR48WJt27ZNjRo1Us2aNfX73//+hiEZwM3ZTFFH4wIod3Jy\nchQQEKAff/yxwPiskrZr1y6Fh4crNze3zKy3BgAlyZK/jGlpaercubNatGihli1bavLkyZKk+Ph4\nBQUFKTIyUpGRkVq2bJnje8aNG6cmTZooNDS0QJcHANdatGiRzpw5o9OnT+vFF19Uq1atXBLevvrq\nK50/f17Z2dkaOXKkevfuTXgDgOuwpAXu8OHDOnz4sCIiIpSTk6M2bdpowYIFmjdvnvz8/PT8888X\nOD45OVmDBg3Spk2blJ6ernvvvVcpKSn8cQdKwVNPPaX58+fLGKM777xTU6ZMUZMmTUr8Oj169ND6\n9etVoUIFxcbGasqUKY4B+wCAgrytuGjt2rUdSwD4+vqqefPmSk9Pl3TtwbFff/21Bg4cKB8fHwUH\nB6tx48bauHGj2rVrV6p1A+XR1KlTNXXqVJdf5/IWdwDAjVnehJWamqqkpCRHGHvvvffUunVrPfHE\nEzpx4oQk+1pTl68LFRQU5Ah8AAAA5Y0lLXCX5OTk6KGHHtKkSZPk6+urZ555RqNGjZIkvfbaa3rh\nhRccewhe6VqzqNiAGQAAeJKijmSzLMDl5eWpf//+euSRR9S3b19JBVf2fvLJJ9WrVy9JUr169Qos\n7Hnw4EHVq1fvmudlUq1nio+PV3x8vNVloIi4f56N++e5uHeerTgNT5Z0oRpj9MQTTygsLKzAGkOH\nDh1yPP/qq68UHh4uSerdu7fmzJmj3NxcHThwQHv37lV0dHSp1w0AAOAOLGmBW7t2rT799FO1atVK\nkZGRkqSxY8dq9uzZ2rZtm2w2m0JCQvTRRx9JksLCwhQXF6ewsDB5e3trypQpdJcCAIByq0wt5Guz\n2ehC9VCJiYmKjY21ugwUEffPs3H/PBf3zrMVJ7cQ4AAAACxQnNxi+TIiAAAAcA4BDgAAwMMQ4AAA\nADwMAQ4AAMDDEOAAAAA8DAEOAADAwxDgAAAAPAwBDgAAwMMQ4AAAADwMAQ4AAMDDEOAAAAA8jLfV\nBQBwP8YU/nH58ZeeX/71es+v9/mVdThTc3GxlTIAV/DxkQICSvacBDigFJ0/L504IZ08af+akyOd\nPl3wceaMdO7c1Y/z56Xc3N8eeXm/fb1w4epHfv5vj4sXf/t65fNLD2N++3opyNhshXtcfuyl55d/\nvd7z631+ueu9X9xjXXkOALhcTIz01Vcle06bMWXn/zltNpvK0I8DD2CMPYj9/LOUlmZ/ZGZKR49K\nR47Yvx49Kh0/bj/u4kXJ31+qVs3+8PW1P6pW/e1RpYpUubJUqZL9UbmyVLGi/XHLLb89fHx+e3h7\n//bV21uqUOHqh5fXb18vf26z/fbepeeXhzEAgGsUJ7cQ4IBCOHlS2rWr4GPfPntg8/KSGjSQ6te3\nPwIDpVq1pJo1f/tavbo9uFWqRDACANgR4H5FgENJOH1a2rJFWr9e2rBB2rjRHuBCQ6XmzX97NG5s\nD27VqlldMQDAExHgfkWAQ1Hk5UnffSctXCitXi3t2SOFh0vt2tkfbdtKDRvaW9oAACgpBLhfEeBQ\nWCdOSMuX20Pb8uXS7bdLvXtLXbpId9xh7+oEAMCVCHC/IsDhRoyR1qyRPvhAWrZM6tTJHtoeeECq\nW9fq6gAA5Q0B7lcEOFxLTo40a5Y0ZYp9eY3hw6VHH2XsGgDAWsXJLawDhzIrK0v629+k6dOl2Fjp\n3Xele+5hFigAwPMxLBtlTl6e9N579lmjZ89KP/wg/fvf9vFthDcAQFlACxzKlGXLpOefl4KCpFWr\npJYtra4IAICSR4BDmXDggH1s2/790sSJ0v3309oGACi76EKFx1uyxL5eW+fO0o4d9lmlhDcAQFlG\nCxw8Vn6+NHq0fZLCv/9t3ywYAIDygAAHj3TsmDR4sJSbK23ebN9/FACA8oIuVHicTZukNm2kiAhp\n5UrCGwCg/KEFDh7l22+luDjp//5PevBBq6sBAMAa7MQAj7F+vX3rqy++sC/MCwCAJytObqELFR4h\nKUnq00f65BPCGwAABDi4veRkqWdP6cMPpR49rK4GAADrEeDg1n78UerWTXr7balfP6urAQDAPRDg\n4LbS0qR775VGjZIeecTqagAAcB9MYoBbysuzL8zbv780cqTV1QAAUPKKk1sIcHBLf/2rtG2btHgx\n22IBAMqm4uQW1oGD20lMtG+PtW0b4Q0AgGthDBzcSlaW9Oij0scfS7VqWV0NAADuiS5UuA1j7Lss\n1Ksnvfuu1dUAAOBadKGiTJg+XUpJkWbNsroSAADcGy1wcAspKfZZp6tXS2FhVlcDAIDrsZUWPFpe\nnjRokDR6NOENAIDCIMDBctOmSdWqSc88Y3UlAAB4BrpQYamcHKlJE2nJEumOO6yuBgCA0kMXKjzW\nxInSPfcQ3gAAcAYtcLDM4cNSixbS5s1SSIjV1QAAULrYSutXBDjPMny4VKmS9I9/WF0JAACljwD3\nKwKc59izx75syJ49UvXqVlcDAEDpYwwcPM5f/iL9+c+ENwAAioIWOJS69evtW2alpEiVK1tdDQAA\n1qAFDh7DGHvL2xtvEN4AACgqAhxK1cKF0qlT0pAhVlcCAIDnIsCh1BgjjRkjvfmmVKGC1dUAAOC5\nCHAoNRs2SCdOSA88YHUlAAB4NgIcSs3770v/7/9JXvxXBwBAsTALFaXi8GGpeXNp/34pIMDqagAA\nsB6zUOH2pk61Lx1CeAMAoPgsCXBpaWnq3LmzWrRooZYtW2ry5MmSpKysLHXt2lVNmzZVt27ddOLE\nCcf3jBs3Tk2aNFFoaKgSEhKsKBtFlJcnffihvfsUAAAUnyVdqIcPH9bhw4cVERGhnJwctWnTRgsW\nLND06dNVo0YNvfTSS3rrrbeUnZ2t8ePHKzk5WYMGDdKmTZuUnp6ue++9VykpKfK6YjAVXajuad48\n6YMPpNWrra4EAAD34XFdqLVr11ZERIQkydfXV82bN1d6eroWLlyooUOHSpKGDh2qBQsWSJK+/vpr\nDRw4UD4+PgoODlbjxo21ceNGK0pHEbz/vvTss1ZXAQBA2WH5GLjU1FQlJSWpbdu2yszMVGBgoCQp\nMDBQmZmZkqSMjAwFBQU5vicoKEjp6emW1Avn/PCDtG+f1Lev1ZUAAFB2eFt58ZycHPXv31+TJk2S\nn59fgc9sNptsNtt1v/d6n8XHxzuex8bGKjY2tiRKRRF98IH09NOSj4/VlQAAYK3ExEQlJiaWyLks\nC3B5eXnq37+/hgwZor6/Ns8EBgbq8OHDql27tg4dOqRatWpJkurVq6e0tDTH9x48eFD16tW75nkv\nD3CwVna29MUX0u7dVlcCAID1rmxYGj16dJHPZUkXqjFGTzzxhMLCwvTcc8853u/du7dmzpwpSZo5\nc6Yj2PXu3Vtz5sxRbm6uDhw4oL179yo6OtqK0uGE6dOl+++Xfu0VBwAAJcSSWajfffedOnXqpFat\nWjm6QseNG6fo6GjFxcXp559/VnBwsObNmyd/f39J0tixY/Xxxx/L29tbkyZNUvfu3a86L7NQ3cfF\ni1KTJtJnn0nt2lldDQAA7qc4uYWdGOASK1dKI0dKW7ZINxjKCABAueVxy4ig7Js7Vxo8mPAGAIAr\n0AKHEpeXJ9WuLW3dKjVsaHU1AAC4J1rg4Fa++UZq2pTwBgCAqxDgUOLmzbNvXA8AAFyDLlSUqNxc\nqU4dads2qX59q6sBAMB90YUKt/Gf/0ihoYQ3AABciQCHEkX3KQAArkcXKkpMbq599umOHdJ1djoD\nAAC/ogsVbmHlSqlFC8IbAACuRoBDiaH7FACA0kEXKkrE+fP27tOdO6W6da2uBgAA90cXKiyXkCC1\nakV4AwCgNBDgUCLoPgUAoPTQhYpiO3fOvnjvrl32blQAAHBzdKHCUitWSBERhDcAAEoLAQ7FRvcp\nAACliy5UFMuFC1LNmlJysr0bFQAAFA5dqLDMpk1SgwaENwAAShMBDsWSkCB17251FQAAlC8EOBTL\nihVSt25WVwEAQPnCGDgU2YkT9u7TI0ekSpWsrgYAAM/CGDhY4ptvpJgYwhsAAKWNAIciS0ig+xQA\nACsQ4FAkxtjHvzGBAQCA0keAQ5Hs3WtfA655c6srAQCg/CHAoUgutb7ZbFZXAgBA+UOAQ5Ew/g0A\nAOuwjAiclpsr1aghHTggVa9udTUAAHgmly8jkpOTI0nKy8tTfn5+kS6EsmPdOik0lPAGAIBVbhrg\nJkyYoDFjxuj555/XyZMn9fTTT5dGXXBjbJ8FAIC1vG92QNu2bdW2bVv5+Pho7ty5unjxYmnUBTe2\nYoX07rtWVwEAQPl10xa4qlWrasaMGapQoYIGDRqkTp06lUZdcFNHj0r79knt2lldCQAA5VexJjFs\n375d4eHhsrnJWhJMYnC9zz+X5s2TFiywuhIAADxbqe6F+sknn+i5557TjBkzVLVqVc2ePbtIF4Zn\nYvkQAACsV6R14EaNGqVatWrp7bff1t69e0u6JrgpY5jAAACAO7jpJIYr1ahRQ7fccot69uypnj17\nuqImuKn//U+qXFm6/XarKwEAoHxzOsAtX75cEyZMUPXq1RUdHa3OnTsrOjraFbXBzfznP1LXrlZX\nAQAAnO5CjY2NVWJiombNmqX27dtr8+bNrqgLbui77yQmIQMAYD2nZ6F+9dVXCgoK0p133umqmoqM\nWaiuY4xUu7a0aZPUoIHV1QAA4PmKk1uc7kJdvXq1JGnMmDGqVKmS7r77bj377LNFujg8x48/ShUr\nEt4AAHAHTge4/v37y2az6a677tLZs2e1c+dOV9QFN/Pdd9Jdd1ldBQAAkArRhbpnzx55eXmpSZMm\npVVTkdGF6jpPPCG1aSMNH251JQAAlA3FyS03DXAXLlxQYmKiI8jdeeedioqKKtLFXI0A5zrNmklf\nfCG1amV1JQAAlA0uDXBX2rhxo7Zs2aKLFy+qWbNmio2Nlbe30z2xLkGAc40jR6SmTaXjx6UKFayu\nBgCAsqFUA9zl9uzZo8TEROXm5qpevXrq3r27qlatWtTTFRsBzjW++kqaOlVautTqSgAAKDtcPgs1\nOTlZYWFhV73frFkzNWvWTJKUkZGhxYsX63e/+12RCoH7YgIDAADupVAL+T755JM6fvz4DY+pW7cu\n4a2MWruWAAcAgDspVBfqvHnzVL9+fR0/flwxMTEKCAgojdqcRhdqyTtzRqpVSzp61L4PKgAAKBku\n70KNi4tzPF+3bp2ysrJ01113yd/fv0gXhefYuNE+85TwBgCA+yhUF+rs2bMdzyMiIhQZGam5c+fq\nnXfe0alTp1xWHKzH+DcAANxPobpQfX19VaVKFVWsWFG+vr7y9/dXQECA/P391aRJE40ePbo0ar0p\nulBL3n332Rfv7d3b6koAAChbXL6MyNy5c9WtWzctXbpU1atX13333Veki7kaAa5k5edLt90m7dsn\n1ahhdTUAAJQtLg9wZ86cUZUqVSTZlwtZsGCBGjZsqPvvv79IF3UVAlzJ2rZNGjhQ2rXL6koAACh7\nipNbCjUG7vHHH9fnn3+uzz77TN9++62qVaumDRs26O6779bixYuLdGG4P8a/AQDgngo1C3XLli2S\n5Bj3FhAQoKCgIA0fPlyVmZ5YZn33ndSjh9VVAACAKxWqC3X79u1q5QG7mNOFWnKMkYKCpP/+V7r9\ndqurAQCg7LFsL1R3Q4ArOampUvv2UkaGZLNZXQ0AAGWPy8fAofy5NP6N8AYAgPuxJMANGzZMgYGB\nCg8Pd7wXHx+voKAgRUZGKjIyUsuWLXN8Nm7cODVp0kShoaFKSEiwouRyhwkMAAC4L6e6UM+dO6cv\nv/xSqampunDhgv0ENptGjRrl1EXXrFkjX19fPfroo9qxY4ckafTo0fLz89Pzzz9f4Njk5GQNGjRI\nmzZtUnp6uu69916lpKTIy+vq7EkXaslp2VKaOVNq08bqSgAAKJtKrQu1T58+WrhwoXx8fOTr6ytf\nX19VrVrV6Yt27NhRAQEBV71/rR/i66+/1sCBA+Xj46Pg4GA1btxYGzdudPqaKLyTJ+1j4Fq3troS\nAABwLYVaRuSS9PR0rVixwlW16L333tMnn3yiqKgoTZw4Uf7+/srIyFC7du0cxwQFBSk9Pd1lNUDa\nvFmKjJS8nfqvAwAAlBan/onu0KGDy5YUeeaZZxxdsa+99ppeeOEFTZs27ZrH2m4wsj4+Pt7xPDY2\nVrGxsSVZZrmwaZN0551WVwEAQNmSmJioxMTEEjmXUwFuzZo1mj59ukJCQlSxYkVJ9jC1ffv2YhdS\nq1Ytx/Mnn3xSvXr1kiTVq1dPaWlpjs8OHjyoevXqXfc8lwc4FM3GjVJcnNVVAABQtlzZsDR69Ogi\nn8upAHf5zNCSdujQIdWpU0eS9NVXXzlmqPbu3VuDBg3S888/r/T0dO3du1fR0dEuqwP2Fri337a6\nCgAAcD1OBbjg4OASuejAgQO1evVqHTt2TPXr19fo0aOVmJiobdu2yWazKSQkRB999JEkKSwsTHFx\ncQoLC5O3t7emTJlywy5UFM+hQ9KZM1KjRlZXAgAArqdQy4jExMRo7dq18vX1vSo82Ww2nTp1ymUF\nOoNlRIpv4UJpyhRp+XKrKwEAoGwrTm4pVAvc2rVrJUk5OTlFugg8BxMYAABwf2ylhQI2biTAAQDg\n7tjMHg7GSDVqSP/7n/TrfBIAAOAibGaPErF/v1SlCuENAAB351SAS05Ovuq9klqQDtaj+xQAAM/g\nVICLi4vTW2+9JWOMzpw5oxEjRujll192VW0oZZs2SSyxBwCA+3MqwH3//fdKS0tT+/btFR0drTp1\n6mjdunWuqg2ljBmoAAB4BqcCnLe3typXrqyzZ8/q3LlzatSokby8GEZXFly4ICUlSW3aWF0JAAC4\nGafSV3R0tCpVqqTNmzdrzZo1+vzzz/Xwww+7qjaUouRkKShI8ve3uhIAAHAzTm2lNW3aNEVFRUmS\n6tSpo4ULF2rWrFkuKQyli+5TAAA8h1MBbsmSJVqyZInjNXuSlh3MQAUAwHM41YVatWpV+fr6ytfX\nVxUqVNAkZOSTAAAgAElEQVTSpUuVmprqotJQmpiBCgCA5yjWTgznz59Xt27dtHr16pKsqcjYiaFo\nzp6VqleXsrKkSpWsrgYAgPLBsp0YTp8+rfT09OKcAm5g2zYpNJTwBgCAp3BqDFx4eLjj+cWLF3Xk\nyBGNGjWqxItC6aL7FAAAz+JUgFu0aNFv3+jtrcDAQPn4+JR4UShdmzZJsbFWVwEAAAqrWGPg3A1j\n4IqmWTNp/nzpsgZWAADgYsXJLYUKcL6+vtddMsRms+nUqVNFunhJI8A578QJqX59KTtb8naqPRYA\nABRHcXJLof7JzsnJKdLJ4f62bJEiIghvAAB4kkLNQh0yZIgkadKkSS4tBqVv40YmMAAA4GkKFeC2\nbNmijIwMTZs2TVlZWVc94Lm2bGEDewAAPE2hOs6efvppdenSRfv371ebK/61t9ls2r9/v0uKg+tt\n3Sr97W9WVwEAAJzh1CzUp59+Wh9++KEr6ykWJjE4JztbatBAOnlS8irWks4AAMBZpbYTgzuHNzhv\n2zb7BAbCGwAAnoV/usuxrVulyEirqwAAAM4iwJVjW7dKd9xhdRUAAMBZTgW4kSNHFuo9eAYCHAAA\nnsmpAJeQkHDVe0uXLi2xYlB6cnKkn36Smje3uhIAAOCsQi0j8s9//lNTpkzRvn37FH7Zhpm//PKL\nYmJiXFYcXOeHH6SWLSUfH6srAQAAzirUMiInT55Udna2Xn75Zb311luOKa9+fn6qXr26y4ssLJYR\nKbz33pN27pSYWAwAgDVcvpn95bKzs7V3716dO3fO8V6nTp2KdPGSRoArvMcfl9q3l37/e6srAQCg\nfHL5ZvaXTJ06VZMnT9bBgwcVERGhDRs2qH379lq1alWRLg7rbN0q/b//Z3UVAACgKJyaxDBp0iRt\n3LhRDRs21LfffqukpCRVq1bNVbXBRc6dk/butY+BAwAAnsepAFepUiVVrlxZknTu3DmFhoZqz549\nLikMrrNjh9S0qVSpktWVAACAonCqC7V+/frKzs5W37591bVrVwUEBCg4ONhFpcFVWP8NAADP5vQk\nhksSExN16tQp3XfffbrllltKuq4iYRJD4Tz9tL379Nlnra4EAIDyq9Q2s7948aJmzZqlMWPGKDY2\nVhEREdq2bVuRLgzr0AIHAIBnc6oF7umnn5aXl5dWrVql3bt3KysrS926ddPmzZtdWWOh0QJ3c3l5\nkr+/dOSIVLWq1dUAAFB+ldoyIt9//72SkpIUGRkpSbrtttuUl5dXpAvDGrt2SQ0bEt4AAPBkTnWh\n3nLLLcrPz3e8Pnr0qLy8nDoFLLZ1q/Rr/gYAAB7KqfQ1YsQIPfjggzpy5Ij+8pe/KCYmRq+88oqr\naoMLMP4NAADP5/Qs1F27dmnVqlUyxqhLly5q3ry5q2pzGmPgbu6uu6Q33pA6d7a6EgAAyrdS2wv1\n3Llz+vLLL5WamqoLFy44Lj5q1KgiXbykEeBuLD/fPoEhLc3+FQAAWKfUJjH06dNH/v7+atOmjSqx\njL/H2btXqlWL8AYAgKdzKsClp6drxYoVrqoFLsb4NwAAyganJjF06NBB27dvd1UtcDECHAAAZUOh\nWuDCw8MlSfn5+Zo+fbpCQkJUsWJFSfb+W0KdZ0hKkl56yeoqAABAcRVqEkNqauoNP3eXDe2ZxHB9\nxki33Sbt2WMfBwcAAKzl8kkM7hLQUHSpqfbdFwhvAAB4PrZRKCcY/wYAQNlBgCsnkpIIcAAAlBUE\nuHKCPVABACg7nFoHbuLEiQUG3NlsNlWrVk1t2rRRRESESwpEyUhKIsABAFBWOLWV1qBBg7R582b1\n6tVLxhgtWbJE4eHh+umnn/TQQw9p5MiRrqz1ppiFem2HDknh4dLRo5LNZnU1AABAKsW9UDt27Khl\ny5bJ19dXkpSTk6OePXtq+fLlatOmjXbt2lWkIkoKAe7ali6V3nlHWrnS6koAAMAlxcktTo2BO3r0\nqG655RbHax8fH2VmZqpKlSrsjerGmIEKAEDZ4tQYuMGDB6tt27bq27evjDFatGiRBg0apNOnTyss\nLMxVNaKYkpKkhx+2ugoAAFBSnOpClaRNmzZp3bp1kqSYmBhFRUW5pLCioAv12kJCpBUrpKZNra4E\nAABcUmpdqOfOnVNKSopycnJ04sQJLVmyRGPGjHH6osOGDVNgYKBjj1VJysrKUteuXdW0aVN169ZN\nJ06ccHw2btw4NWnSRKGhoUpISHD6euVZdrZ0/LjUuLHVlQAAgJLiVIDr06ePFi5cKB8fH/n6+srX\n11dVq1Z1+qKPP/64li9fXuC98ePHq2vXrkpJSVGXLl00fvx4SVJycrLmzp2r5ORkLV++XMOHD9fF\nixedvmZ5lZQktWolebHiHwAAZYZTY+DS09O1YsWKYl+0Y8eOSk1NLfDewoULtXr1aknS0KFDFRsb\nq/Hjx+vrr7/WwIED5ePjo+DgYDVu3FgbN25Uu3btil1HecAODAAAlD1Otct06NBB27dvd0khmZmZ\nCgwMlCQFBgYqMzNTkpSRkaGgoCDHcUFBQUpPT3dJDWUROzAAAFD2ONUCt2bNGk2fPl0hISGqWLGi\nJPsAvJIOdTabTbYbrDh7o8/i4+Mdz2NjYxUbG1uClXmepCTppZesrgIAACQmJioxMbFEzuVUgFu2\nbFmJXPRaAgMDdfjwYdWuXVuHDh1SrVq1JEn16tVTWlqa47iDBw+qXr161z3P5QGuvDt9WkpNlVjh\nBQAA613ZsDR69Ogin8upLtTg4OBrPkpC7969NXPmTEnSzJkz1bdvX8f7c+bMUW5urg4cOKC9e/cq\nOjq6RK5Z1m3fbg9vPj5WVwIAAEpSoVrgYmJitHbtWvn6+l7VfWmz2XTq1CmnLjpw4ECtXr1ax44d\nU/369TVmzBi9/PLLiouL07Rp0xQcHKx58+ZJksLCwhQXF6ewsDB5e3trypQpN+xCxW8Y/wYAQNnk\n9EK+7oyFfAt68kmpTRvpmWesrgQAAFypOLnFqTFwkvTDDz/ov//9r2w2mzp27KjWrVsX6cJwva1b\n7SEOAACULU6NgZs0aZIGDx6so0ePKjMzU4888ogmT57sqtpQDLm50u7d9kV8AQBA2eJUF2p4eLg2\nbNjg2H3h9OnTateunXbs2OGyAp1BF+pvkpKkRx6Rdu60uhIAAHAtpbYXqiR5XbYnkxf7M7mtpCQm\nMAAAUFY5NQbu8ccfV9u2bdWvXz8ZY7RgwQINGzbMVbWhGNhCCwCAssvpWahbtmzR2rVrJdn3NI10\no2YeulB/ExMjvfmm1Lmz1ZUAAIBrKU5ucSrATZw4scDFbDabqlWrpjZt2igiIqJIBZQkApxdfr7k\n7y+lpdm/AgAA91NqY+C2bNmiDz/8UBkZGUpPT9dHH32kZcuW6amnntJbb71VpAJQ8vbulWrWJLwB\nAFBWOdUC17FjRy1btky+vr6SpJycHPXs2VPLly9XmzZttGvXLpcVWhi0wNnNni3Nny99+aXVlQAA\ngOsptRa4o0eP6pZbbnG89vHxUWZmpqpUqaJKlSoVqQCUvK1bmcAAAEBZ5tQs1MGDB6tt27bq27ev\njDFatGiRBg0apNOnTyssLMxVNcJJW7dKL7xgdRUAAMBVnJ6FumnTJq1du1Y2m00xMTGKiopyVW1O\nowtVMka67Tb7LgyBgVZXAwAArqfUZqGOHDnyqskK13rPKgQ46ccfpXvukX7+2epKAADAjZTaGLiE\nhISr3lu6dGmRLgzX2LxZcqNGUQAA4AKFGgP3z3/+U1OmTNG+ffsUHh7ueP+XX35RTEyMy4qD8whw\nAACUfYXqQj158qSys7P18ssv66233nI09/n5+al69eouL7Kw6EKVYmOlv/xF6tbN6koAAMCNlNoY\nOHdX3gPcxYv2xXsPHJDcKFcDAIBrKE5ucWoZkdGjR1/z4qNGjSrSxVGyUlKkGjUIbwAAlHVOBbiq\nVavKZrNJks6ePavFixez/psbYfwbAADlQ7G6UM+fP69u3bpp9erVJVlTkZX3LtTnnpPq1pVeesnq\nSgAAwM2U2jIiVzp9+rTS09OLcwqUIFrgAAAoH5zqQr18CZGLFy/qyJEjjH9zExcuSNu2sQcqAADl\ngVMBbt68ebpw4YIkyd/fX3Xq1JG3t1OngIvs3m3vPvX3t7oSAADgaoVKX3l5efrrX/+qjz/+WA0a\nNJAk/fzzz3r88cc1duxY+fj4uLRI3BzdpwAAlB+FGgP35z//WVlZWTpw4IC2bt2qrVu3av/+/Tpx\n4oRefPFFV9eIQiDAAQBQfhRqFmrjxo2VkpIiL6+CeS8/P1/NmjXTjz/+6LICnVGeZ6G2aydNmCB1\n6mR1JQAAoDBcPgvVy8vrqvAmSRUqVLjm+yhdeXnSjh1SZKTVlQAAgNJQqPTVvHlzzZw586r3Z82a\npdDQ0BIvCs7ZuVNq2FDy87O6EgAAUBoK1YV68OBB9evXT5UrV1abNm0kSVu2bNGZM2f01VdfKSgo\nyOWFFkZ57UL917+k//5X+uQTqysBAACFVSqb2RtjtGrVKu3cuVM2m01hYWHq0qVLkS7qKuU1wD39\ntBQWJv3hD1ZXAgAACqtUApwnKK8BLipKmjxZ6tDB6koAAEBhEeB+VR4D3PnzUkCAdOyYVKWK1dUA\nAIDCsmwvVFhvxw6pcWPCGwAA5QkBzsOxgC8AAOUPAc7DEeAAACh/CHAejgAHAED5wyQGD3b2rFS9\nupSVJVWqZHU1AADAGUxiKKe2bZNCQwlvAACUNwQ4D7ZhgxQdbXUVAACgtBHgPNi6dVJMjNVVAACA\n0kaA81DGEOAAACivCHAe6uefpfx8KSTE6koAAEBpI8B5qHXr7Huf2mxWVwIAAEobAc5DrV3L5vUA\nAJRXBDgPxfg3AADKLxby9UA5OVJgoH0B34oVra4GAAAUBQv5ljMbN0oREYQ3AADKKwKcB7o0gQEA\nAJRPBDgPxPg3AADKN8bAeZiLF+0b2O/ebR8HBwAAPBNj4MqRXbvsAY7wBgBA+UWA8zCMfwMAAAQ4\nD0OAAwAABDgPwwQGAABAgPMgR49KmZlSWJjVlQAAACsR4DzI+vVS27ZShQpWVwIAAKxEgPMgjH8D\nAAASAc6jMP4NAABIbriQb3BwsG699VZVqFBBPj4+2rhxo7KysvS73/1OP/30k4KDgzVv3jz5+/tf\n9b1leSHf3FzpttukjAzp1lutrgYAABRXmVrI12azKTExUUlJSdq4caMkafz48eratatSUlLUpUsX\njR8/3uIqS9+2bVLjxoQ3AADghgFO0lVpdOHChRo6dKgkaejQoVqwYIEVZVlq7VrGvwEAADu3C3A2\nm0333nuvoqKiNHXqVElSZmamAn/dOyowMFCZmZlWlmgJJjAAAIBLvK0u4Epr165VnTp1dPToUXXt\n2lWhoaEFPrfZbLLZbNf9/vj4eMfz2NhYxcbGuqjS0mOM9N130oQJVlcCAACKKjExUYmJiSVyLreb\nxHC50aNHy9fXV1OnTlViYqJq166tQ4cOqXPnztq9e/dVx5fVSQz/+5/Up4+0b5/VlQAAgJJSZiYx\nnDlzRr/88osk6fTp00pISFB4eLh69+6tmTNnSpJmzpypvn37WllmqUtIkLp2tboKAADgLtyqCzUz\nM1MPPvigJOnChQsaPHiwunXrpqioKMXFxWnatGmOZUTKk5UrpaeesroKAADgLty6C9VZZbEL9dw5\nqWZNKS1NusbSdwAAwEOVmS5UXG3dOqllS8IbAAD4DQHOzTH+DQAAXIkA5+ZWrpS6dbO6CgAA4E4Y\nA+fGjh61b5917Jjk42N1NQAAoCQxBq6M+uYb6e67CW8AAKAgApwbo/sUAABcCwHOTRnDBAYAAHBt\nBDg3tWeP5OUlNW1qdSUAAMDdEODc1KXWN5vN6koAAIC7IcC5Kca/AQCA62EZETeUm2vfPmv/fql6\ndaurAQAArsAyImXMhg1SkyaENwAAcG0EODeUkED3KQAAuD4CnBtauZLlQwAAwPUxBs7NZGVJwcH2\nbbQqVrS6GgAA4CqMgStDVq2S7rqL8AYAAK6PAOdmvvpKuv9+q6sAAADujC5UN3LmjFS3rpSSItWq\nZXU1AADAlehCLSMWL5batiW8AQCAGyPAuZHZs6UBA6yuAgAAuDu6UN3EyZNSgwbSTz9J/v5WVwMA\nAFyNLtQy4KuvpM6dCW8AAODmCHBuYs4caeBAq6sAAACegC5UN3D0qH3v0/R0qWpVq6sBAAClgS5U\nDzd/vtSzJ+ENAAAUDgHODTD7FAAAOIMuVIsdPCi1bi1lZLB9FgAA5QldqB5s7lypb1/CGwAAKDwC\nnMWYfQoAAJxFgLPQ3r1SWpoUG2t1JQAAwJMQ4Cw0d6708MOSt7fVlQAAAE9CgLOIMcw+BQAARUOA\ns8iaNVJurtS+vdWVAAAAT0OAs8i4cdKf/yx5cQcAAICTWAfOAtu2SfffL+3fz/IhAACUV6wD52HG\nj5f+9CfCGwAAKBpa4ErZjz/ax73t3y/5+VldDQAAsAotcB5kwgTpmWcIbwAAoOhogStFGRlSixZS\nSopUs6bV1QAAACvRAuch3nlHevRRwhsAACgeWuBKSXa21LixlJQkNWhgdTUAAMBqtMB5gPffl3r3\nJrwBAIDiowWuFJw+LTVqJK1eLYWGWl0NAABwB7TAubmpU6W77iK8AQCAkkELnIulpkp33il9+63U\nsqXV1QAAAHdBC5ybys+3zzp96SXCGwAAKDkEOBf6xz8km016/nmrKwEAAGUJXagusn271KWLtGmT\nFBxsdTUAAMDd0IXqZs6flx55RHr7bcIbAAAoebTAucBLL9k3rf/yS3sXKgAAwJWKk1u8S7iWcm/1\naunTT6UffiC8AQAA16ALtQQdPy499ph93Tf2OwUAAK5CgCshP/8sdewoDRok3X+/1dUAAICyjABX\nAnbskGJipKeekv72N6urAQAAZR1j4Ipp9Wrp4YelyZOlAQOsrgYAAJQHtMAVw/z59vA2ezbhDQAA\nlB5a4IogN1eaNEl6910pIUGKiLC6IgAAUJ54TAvc8uXLFRoaqiZNmuitt96ypIbcXPsM06ZNpZUr\npe++I7wBAIDS5xEBLj8/X88++6yWL1+u5ORkzZ49W7t27Sq1618e3ObPlz7/3N7yFhJSaiWUeYmJ\niVaXgGLg/nk27p/n4t6VXx4R4DZu3KjGjRsrODhYPj4+GjBggL7++muXXvOXX6QVK6RXX5WaNfst\nuK1YIXXo4NJLl0v8EfJs3D/Pxv3zXNy78ssjxsClp6erfv36jtdBQUH6/vvvi31eY6ScHCkzUzpy\nRDp4UFq/XlqzRtq9W2rTxr6225w5Utu2xb4cAABAifCIAGdzYk+qBx6wB7NLLl6U8vIKPnJzpVOn\n7MHNZpMCA+2POnWkO++0T1CIipIqVnTBDwMAAFBMHrGZ/YYNGxQfH6/ly5dLksaNGycvLy+NHDmy\nwHHOBD0AAACrFTWGeUSAu3Dhgpo1a6ZvvvlGdevWVXR0tGbPnq3mzZtbXRoAAECp84guVG9vb73/\n/vvq3r278vPz9cQTTxDeAABAueURs1AlqUePHtqzZ49+/PFHvfLKKwU+c4c14lA4aWlp6ty5s1q0\naKGWLVtq8uTJkqSsrCx17dpVTZs2Vbdu3XTixAmLK8WN5OfnKzIyUr169ZLE/fMkJ06c0EMPPaTm\nzZsrLCxM33//PffPQ4wbN04tWrRQeHi4Bg0apPPnz3Pv3NiwYcMUGBio8PBwx3s3ul/jxo1TkyZN\nFBoaqoSEhJue32MC3PVYvUYcnOPj46N33nlHO3fu1IYNG/TBBx9o165dGj9+vLp27aqUlBR16dJF\n48ePt7pU3MCkSZMUFhbmGHfK/fMcf/zjH9WzZ0/t2rVL27dvV2hoKPfPA6Smpmrq1KnaunWrduzY\nofz8fM2ZM4d758Yef/xxx9j9S653v5KTkzV37lwlJydr+fLlGj58uC5evHjjCxgPt27dOtO9e3fH\n63Hjxplx48ZZWBGc0adPH7Ny5UrTrFkzc/jwYWOMMYcOHTLNmjWzuDJcT1pamunSpYtZtWqVeeCB\nB4wxhvvnIU6cOGFCQkKuep/75/6OHz9umjZtarKyskxeXp554IEHTEJCAvfOzR04cMC0bNnS8fp6\n92vs2LFm/PjxjuO6d+9u1q9ff8Nze3wL3LXWiEtPT7ewIhRWamqqkpKS1LZtW2VmZiowMFCSFBgY\nqMzMTIurw/X86U9/0ttvvy0vr9/+fHD/PMOBAwdUs2ZNPf7447rjjjv01FNP6fTp09w/D3Dbbbfp\nhRdeUIMGDVS3bl35+/ura9eu3DsPc737lZGRoaCgIMdxhckyHh/gWDrEM+Xk5Kh///6aNGmS/Pz8\nCnxms9m4r25q8eLFqlWrliIjI6879Z37574uXLigrVu3avjw4dq6dauqVq16VZcb98897du3T+++\n+65SU1OVkZGhnJwcffrppwWO4d55lpvdr5vdS48PcPXq1VNaWprjdVpaWoEUC/eTl5en/v37a8iQ\nIerbt68k+/+JHD58WJJ06NAh1apVy8oScR3r1q3TwoULFRISooEDB2rVqlUaMmQI989DBAUFKSgo\nSHfeeack6aGHHtLWrVtVu3Zt7p+b27x5szp06KDq1avL29tb/fr10/r167l3HuZ6fyuvzDIHDx5U\nvXr1bngujw9wUVFR2rt3r1JTU5Wbm6u5c+eqd+/eVpeF6zDG6IknnlBYWJiee+45x/u9e/fWzJkz\nJUkzZ850BDu4l7FjxyotLU0HDhzQnDlzdM8992jWrFncPw9Ru3Zt1a9fXykpKZKk//znP2rRooV6\n9erF/XNzoaGh2rBhg86ePStjjP7zn/8oLCyMe+dhrve3snfv3pozZ45yc3N14MAB7d27V9HR0Tc+\nWUkP2LPC0qVLTdOmTc3tt99uxo4da3U5uIE1a9YYm81mWrdubSIiIkxERIRZtmyZOX78uOnSpYtp\n0qSJ6dq1q8nOzra6VNxEYmKi6dWrlzHGcP88yLZt20xUVJRp1aqVefDBB82JEye4fx7irbfeMmFh\nYaZly5bm0UcfNbm5udw7NzZgwABTp04d4+PjY4KCgszHH398w/v1t7/9zdx+++2mWbNmZvny5Tc9\nv0fsxAAAAIDfeHwXKgAAQHlDgAMAAPAwBDgAAAAPQ4ADAADwMAQ4AKXm+PHjioyMVGRkpOrUqaOg\noCBFRkbKz89Pzz77bKnU8MMPP2jZsmUldr77779fp06dKrHzAUBhMAsVgCVGjx4tPz8/Pf/886V6\n3RkzZmjLli167733SvW6rnThwgV5e3tbXQaAUkQLHADLXPr/x8TERPXq1UuSFB8fr6FDh6pTp04K\nDg7Wv//9b7344otq1aqVevTooQsXLkiStmzZotjYWEVFRem+++5zrG5+uS+++ELh4eGKiIhQbGys\n8vLyNGrUKM2dO1eRkZH64osvdPr0aQ0bNkxt27bVHXfcoYULF0qyB70+ffqoc+fOatq0qcaMGXPN\nnyE4OFhZWVlKTU1V8+bN9fvf/14tW7ZU9+7dde7cuQLH/vLLL2rUqJHjZzh16pQaNWqk/Px87du3\nTz169FBUVJQ6deqkPXv2SJIWLVqkdu3a6Y477lDXrl115MgRx+9pyJAhuuuuuzR06FDt3LlT0dHR\nioyMVOvWrfXjjz8W9/YAcGcuWr8OAG4oPj7e/P3vfzfGGPPtt9+aBx54wBhjzOuvv246duxoLly4\nYH744QdTuXJlx6KWDz74oFmwYIHJzc017du3N8eOHTPGGDNnzhwzbNiwq64RHh5uMjIyjDHGnDx5\n0hhjzIwZM8yIESMcx7zyyivm008/NcYYk52dbZo2bWpOnz5tpk+fburUqWOysrLM2bNnTcuWLc3m\nzZuvukZwcLA5fvy4OXDggPH29jY//PCDMcaYuLg4x3kv9/jjj5sFCxYYY4z56KOPzIsvvmiMMeae\ne+4xe/fuNcYYs2HDBnPPPfc4arpk6tSp5oUXXnD8nqKiosy5c+eMMcaMGDHCfPbZZ8YYY/Ly8szZ\ns2ev/8sH4PFocwfgVmw2m3r06KEKFSqoZcuWunjxorp37y5JCg8PV2pqqlJSUrRz507de++9kqT8\n/HzVrVv3qnPFxMRo6NChiouLU79+/STZW/3MZSNHEhIStGjRIv3973+XJJ0/f14///yzbDabunXr\npoCAAElSv3799N1336lNmzbXrT0kJEStWrWSJLVp00apqalXHfPkk09qwoQJ6tOnj2bMmKF//etf\nysnJ0bp16/Twww87jsvNzZVk3985Li5Ohw8fVm5urho1auT4PfXu3VsVK1aUJLVv315/+9vfdPDg\nQfXr10+NGzcuxG8bgKciwAFwO7fccoskycvLSz4+Po73vby8dOHCBRlj1KJFC61bt+6G5/nnP/+p\njRs3asmSJWrTpo22bNlyzeP+/e9/q0mTJgXe+/777wu8NsbIy+vGo04uhSlJqlChgs6ePXvVMR06\ndFBqaqoSExOVn5+vsLAwnTp1SgEBAUpKSrrq+BEjRujFF1/UAw88oNWrVys+Pt7xWZUqVRzPBw4c\nqHbt2mnx4sXq2bOnPvroI3Xu3PmG9QLwXIyBA+BWTCHmVTVr1kxHjx7Vhg0bJEl5eXlKTk6+6rh9\n+/YpOjpao0ePVs2aNXXw4EHdeuut+uWXXxzHdO/eXZMnT3a8vhSijDFauXKlsrOzdfbsWX399deK\niYkp7o8nSXr00Uc1ePBgDRs2TJJ06623KiQkRPPnz3dce/v27ZLs4+QutS7OmDHDcY4rf08HDhxQ\nSEiIRowYoT59+mjHjh0lUisA90SAA2AZm83m+Hqt55cfc/lrHx8fzZ8/XyNHjlRERIQiIyO1fv36\nq87/0ksvqVWrVgoPD1dMTIxatWqlzp07Kzk52TGJ4bXXXlNeXp5atWqlli1b6vXXX3dcJzo6Wv37\n974lLzYAAADLSURBVFfr1q310EMP6Y477rjuz3C9Wq9l0KBBys7O1sCBAx3vffbZZ5o2bZoiIiLU\nsmVLx2SK+Ph4Pfzww4qKilLNmjWv+3uaN2+eWrZsqcjISO3cuVOPPvroNa8NoGxgGREAuAZXLjcy\nf/58LVq0SDNnzizxcwMoHxgDBwDXcGULV0kZMWKEVqxYoaVLl5b4uQGUH7TAAQAAeBjGwAEAAHgY\nAhwAAICHIcABAAB4GAIcAACAhyHAAQAAeBgCHAAAgIchwAEAAHiY/w8XcMsW8jwrfwAAAABJRU5E\nrkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x105880910>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig = plt.figure(2, figsize=(10, 6))\n",
      "plt.plot(times_years[:100], T[:100])\n",
      "plt.title('Temperature of the planet over time')\n",
      "plt.xlabel('Time step in years')\n",
      "plt.ylabel('Temperature in K')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "<matplotlib.text.Text at 0x1058b7e90>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAGHCAYAAADr18W9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VGXexvF7UkAgoUMICRiQGkIkgFFAJJTQpAmCgEIU\nXVd9RV30ta2sYdcF3NUFdYVlBQEbRVTaC6FJXEERpIgaMQhEIIUiNRJIe94/ZhkJaRAyOVO+n+vK\nlZkzZ875zRwuuHnasRljzkgKFAAAAFyGzRhjrC4CAAAABflYXQAAAAAKI6QBAAC4IEIaAACACyKk\nAQAAuCBCGgAAgAsipAFwOzNnzlRQUJCqV6+ukydPlrr/vHnz1LVr1wqoTAoLC9OGDRsq5FzuJDAw\nUCkpKVaXAbgVQhpQzgICAhQYGKjAwED5+PioatWqjucLFiywurwyCQsL06effmp1GZKknJwcPfnk\nk9qwYYPOnDmjWrVqFXg9JSVFPj4+ys/Pt6Q+m80mm83m9PP4+Pho//79Tj9PWcTExGjOnDkFtp09\ne1ZhYWHWFAS4KT+rCwA8TWZmpuNxkyZNNGfOHPXo0cPCikqWl5cnX1/fEvex2Wy6liUVr+QcVyoj\nI0Pnz59X69atS9zPG5aAdIXPmJ+fLx+fgv/fr4iQCngDWtKACpKfn6+pU6eqWbNmqlu3ru666y5H\nV93F1p958+apcePGqlOnjv71r39p27ZtioyMVK1atTR+/HjHsebNm6cuXbpo/Pjxqlmzplq3bl2g\npev06dO6//771bBhQ4WGhmrixImOlqWL750wYYLq1q2rSZMmaf/+/erRo4fq1q2revXq6Z577tHp\n06clSWPGjNHBgwc1cOBABQYG6pVXXlFiYqIaNWpU4PNd2toWHx+vO++8U2PGjFGNGjU0f/78Emu6\n3IULF/TEE08oJCREISEh+sMf/qDs7GwlJyc7wlnNmjXVq1evQu+97bbbHK9Xr15dW7ZscYSG//3f\n/1Xt2rXVtGlTJSQkXNH3dbmLn23kyJGqXr26OnTooN27dxe579atW9WpUyfVqlVLDRs21Pjx45WT\nk+N43cfHR7NmzVKLFi1Uq1YtPfroowXe//bbbys8PFy1a9dW3759dfDgwQKf8cYbb1RgYKA+/PDD\nQuc2xuill15SWFiYgoKCFBcXpzNnzkiS+vXrpzfffLPA/jfeeKOWLl0qSdqzZ49iY2NVp04dtWrV\nqsDx7733Xj388MPq37+/AgIClJiYWOA4f/zjH/X555/r0UcfVWBgoB577DHHZ73Y8nfvvffqkUce\nUf/+/RUYGKiuXbsqIyNDjz/+uGrVqqXWrVtr165djmOmpaVp2LBhql+/vpo2bao33nijyO8b8DgG\ngNOEhYWZDRs2GGOMmT59uunUqZNJTU012dnZ5ve//70ZNWqUMcaYAwcOGJvNZh5++GFz4cIFs3bt\nWlOpUiUzZMgQc+zYMZOammrq169vPvvsM2OMMXPnzjV+fn5m+vTpJjc31yxatMjUqFHDnDx50hhj\nzJAhQ8xDDz1kzp07Z44ePWqio6PNrFmzCrz3n//8p8nLyzNZWVnmp59+MuvXrzfZ2dnm2LFj5rbb\nbjNPPPFEkZ/DGGM2btxoQkNDi/2sL774ovH39zfLli0zxhiTlZVVYk2XmzhxounUqZM5duyYOXbs\nmOncubOZOHGiMcaYlJQUY7PZTF5eXpHvLer1uXPnGn9/fzN79myTn59vZs6caRo2bOh4/Wpqu/jZ\nPvroI5Obm2teeeUV06RJE5Obm1voe9i+fbv56quvTF5enklJSTGtW7c206dPdxzLZrOZgQMHmtOn\nT5uDBw+aevXqmYSEBGOMMUuXLjXNmjUze/bsMXl5eeall14ynTt3LvDeffv2FVmjMcbMmTPHNGvW\nzBw4cMBkZmaaoUOHmjFjxhhjjHnnnXdMly5dHPt+//33pmbNmiY7O9tkZmaa0NBQM2/ePJOXl2d2\n7txp6tata5KSkowxxsTFxZkaNWqYL774whhjzPnz5wudOyYmxsyZM6fAtkvrjYuLM3Xr1jU7duww\n58+fNz169DDXX3+9effdd01+fr554YUXTPfu3Y0xxuTl5Zn27dubv/zlLyYnJ8fs37/fNG3a1KxZ\ns6bYzw54CkIa4ESX/oPdunXrAkEnLS3N+Pv7m7y8PEdIS0tLc7xep04ds3jxYsfzYcOGOf6Bnzt3\nboGQYYwx0dHR5t133zUZGRmmcuXKJisry/HaBx984PhHb+7cuaZx48Yl1v3JJ5+YqKioIj+HMVcW\n0rp16+Z4rbSaLnfDDTeY1atXO56vWbPGhIWFGWN+C7TFhbSiXp87d65p1qyZ4/mvv/5qbDabOXLk\nyFXX9uKLL5pOnTo5nufn55vg4GCzadOmQt/D5aZNm2buuOMOx3ObzWY2b97seD5ixAjz8ssvG2OM\n6du3b4Ggk5eXZ6pWrWoOHjzoeG9JIa1Hjx5m5syZjuc//vij48/bmTNnTLVq1RzHev755839999v\njDFm4cKFpmvXrgWO9eCDD5pJkyYZY+wBKy4urtjzGmMPabNnzy6w7dJ67733XvPggw86XnvjjTdM\neHi44/nu3btNzZo1jTHGbNmypdCf18mTJ5v77ruvxBoAT8CYNKCCpKSk6I477igwfsfPz09Hjhxx\nPA8KCnI8rlKlSqHnv/76q+N5SEhIgeNff/31SktL08GDB5WTk6Pg4GDHa/n5+WrcuLHj+eVdlUeO\nHNHjjz+uTZs26ezZs8rPz1ft2rWv4dNKoaGhjsc///xzqTVdKi0tTddff73jeePGjZWWlnZN9TRo\n0MDxuGrVqpLs4wePHz9+VbVJBT+bzWZTaGhokfUlJydrwoQJ2r59u86dO6fc3Fx17NixxLoujmn8\n+eef9fjjj+vJJ58ssH9qamqh61eU9PT0Qt9hbm6ujhw5ouDgYN1+++1asGCBnn76aS1cuFCzZ892\nnPerr74qMCEjNzdXY8eOLfB5S1PauLT69es7Hl933XUFnlepUqXA95CWllagnry8PEeXL+DJCGlA\nBWncuLHmzp2rTp06FXqtLEsTpKamFnj+888/a/DgwWrUqJEqV66sX375pdCA7osu/wf0+eefl6+v\nr7777jvVrFlTS5cuLTAG7vL9q1WrpnPnzjme5+Xl6dixY8We40pqulTDhg2VkpLiGH928OBBNWzY\nsNT3FVVraa62Nkk6dOiQ43F+fr4OHz5cZH0PP/ywOnTooEWLFqlatWqaPn26Pvrooys6R+PGjTVx\n4kSNGjXqyj7IZS5+hxcdPHhQfn5+juA/atQoTZo0SV27dtX58+fVvXt3x3m7deumtWvXlum8UvlO\nHGjUqJGaNGmi5OTkcjsm4C6YOABUkIceekjPP/+8Y/D3sWPHtHz58qs6hrlkNt/Ro0f1+uuvKycn\nRx9++KH27Nmj/v37q0GDBurdu7cmTJjgaBXbt2+f/vOf/xR73MzMTFWrVk3Vq1dXamqq/v73vxd4\nPSgoSPv27XM8b9Gihc6fP69Vq1YpJydHL730ki5cuFDs8YODg6+qplGjRumll17S8ePHdfz4cf35\nz3/WmDFjrug7qlevnnx8fArUW5KrrU2Stm/frk8++US5ubmaPn26rrvuOt1yyy2F9svMzFRgYKCq\nVq2qPXv2aObMmSXWYuxDUCTZ/7xMnjxZSUlJkuyTGy4dwH/5NbncqFGjNG3aNKWkpCgzM1PPP/+8\nRo4c6Qii/fv3188//6wXX3xRI0eOdLxvwIABSk5O1nvvvaecnBzl5ORo27Zt2rNnj6PG0pRW25Uc\n46Lo6GgFBgbqb3/7m7KyspSXl6fvvvtOX3/99RUfA3BXhDSggjz++OMaNGiQevfurerVq6tTp07a\nunWr4/UraX24dJ+bb75Ze/fuVb169TRx4kR99NFHji6hd955R9nZ2Y6ZgcOHD1dGRobjGJef68UX\nX9SOHTtUo0YNDRw4UMOGDSuwz3PPPaeXXnpJtWrV0j/+8Q/VqFFDM2bM0AMPPKDQ0FAFBAQU6IIr\n6hwl1XS5F154QR07dlRkZKQiIyPVsWNHvfDCC1f0XVWtWlV//OMf1aVLF9WuXVtfffVVkfVc+vxq\narPZbBo8eLAWLVqk2rVr6/3339fHH39c5BIjr7zyij744ANVr15dDz74oEaOHFngvEXVdHHbkCFD\n9Mwzz2jkyJGqUaOG2rZtqzVr1jj2jY+PV1xcnGrVqqUlS5YUOve4ceM0ZswY3XbbbWratKmqVq1a\nYFZkpUqVNHToUG3YsEGjR492bA8ICNDatWu1cOFChYSEKDg4WM8995yys7ML1Vicxx9/XEuWLFHt\n2rX1xBNPFPkdXv49FHd9fH19tXLlSu3atUtNmzZVvXr19OCDDzpmqgKezGau5r80AFzCvHnzNGfO\nHH3++edWl+J1Jk2apJ9++knvvvuu1aUA8HBOa0k7f/68br75ZrVr107h4eF67rnnJEknTpxQbGys\nWrRood69e+vUqVOO90yZMkXNmzdXq1atrmk8BAA4C/+vBVBRnBbSrrvuOm3cuFG7du3S7t27tXHj\nRm3atElTp05VbGyskpOT1bNnT02dOlWSlJSUpEWLFikpKUkJCQl65JFHLLutC+DqKurWQyiM7x5A\nRamQ7s5z586pW7dumjdvnoYNG6bPPvtMQUFBysjIUExMjPbs2aMpU6bIx8dHzzzzjCSpb9++io+P\nL3IwLgAAgKdz6sSB/Px8tWvXTkFBQerevbvatGmjI0eOOKaABwUFOdaISktLK7D2TmhoaKElBgAA\nALyFU9dJ8/Hx0a5du3T69Gn16dNHGzduLPB6ad0GRb1GNwMAAHAnZe20rJDFbGvUqKHbb79d27dv\nd3RzNmjQQOnp6Y5VpkNCQgosEHn48OFCK6pfxMBd9xQfH6/4+Hiry0AZcf3cG9fPfXHt3Nu1NC45\nrbvz+PHjjpmbWVlZWrdunaKiojRo0CDNnz9fkjR//nwNGTJEkjRo0CAtXLhQ2dnZOnDggPbu3avo\n6GhnlQcAAODSnNaSlp6erri4OOXn5ys/P19jxoxRz549FRUVpREjRmjOnDkKCwvT4sWLJUnh4eEa\nMWKEwsPD5efnpxkzZtC1CQAAvJbbLWZrs9no7nRTiYmJiomJsboMlBHXz71x/dwX1869XUtuIaQB\nAAA4ybXkFu7dCQAA4IIIaQAAAC6IkAYAAOCCCGkAAAAuiJAGAADggghpAAAALoiQBgAA4IIIaQAA\nAC6IkAYAAOCCCGkAAAAuiJAGAADggghpAAAALoiQBgAA4IIIaQAAAC6IkAYAAOCCCGkAAAAuiJAG\nAADggghpAAAALoiQBgAA4IIIaQAAAC6IkAYAAOCCCGkAAAAuiJAGAADggghpAAAALoiQBgAA4III\naQAAAC6IkAYAAOCCCGkAAAAuiJAGAADggghpAAAALoiQBgAA4IIIaQAAAC6IkAYAAOCCCGkAAAAu\niJAGAADggghpAAAALoiQBgAA4IIIaQAAAC6IkAYAAOCCCGkAAAAuiJAGAADggghpAAAALoiQBgAA\n4IIIaQAAAC7Iz+oCAE9jjJSXZ//JzZXy8+2PL/6++Dg/377vxceXbrv4c/nzy38unu/SxyX9Lm7b\n5fVf6fbi9r1WzjouADhL9erSjTeW7zGdFtIOHTqksWPH6ujRo7LZbHrwwQf12GOPKT4+XrNnz1a9\nevUkSZMnT1a/fv0kSVOmTNHbb78tX19fvf766+rdu7ezyoMXyc2VTp2STp6UTpyQTp+Wzp61/2Rm\n/vb73DkpK+u331lZ0vnzUna2dOGC/efi45ycgj+5ufbfFwOYj4/k5yf5+tp/fHwK/778x2b77ffl\nj0v6kQo/Lul3cdsudTXbi9v3WjnruADgDJGR0ptvlu8xbcY45/+sGRkZysjIULt27ZSZmakOHTpo\n6dKlWrx4sQIDAzVhwoQC+yclJWn06NHatm2bUlNT1atXLyUnJ8vHp2CPrM1mk5NKhps5d05KSZH2\n75fS0qT0dCkj47ffR47YQ9mvv0o1aki1a0u1atkfBwbafwICfvtdtapUpcpvP1WrSpUrF/ypVOm3\nH3//3378/Oy/L4YyAgYAQLq23OK0lrQGDRqoQYMGkqSAgAC1bt1aqampklRkscuWLdOoUaPk7++v\nsLAwNWvWTFu3btUtt9zirBLhBoyRDh6Udu+2//zwgz2UHThgbxkLC5OaNJFCQ6UGDaSICCk21v64\nfn2pTh17KPNh9CUAwM1UyJi0lJQU7dy5U7fccos2b96sN954Q++88446duyoV199VTVr1lRaWlqB\nQBYaGuoIdfAeGRnSf/4jff65tGuX9O23UrVq9mbkyEipVy/phhukpk2l4GDCFwDAczk9pGVmZurO\nO+/Ua6+9poCAAD388MP605/+JEmaOHGinnzySc2ZM6fI99roM/J4R49K69ZJn31m/zl2TLr1Vum2\n26ShQ6W2baW6da2uEgCAiufUkJaTk6Nhw4bpnnvu0ZAhQyRJ9evXd7z+wAMPaODAgZKkkJAQHTp0\nyPHa4cOHFRISUuRx4+PjHY9jYmIUExNT/sXDaY4elT7+WFq8WNqxQ+rRQ+rWTXrkEXso8/W1ukIA\nAMomMTFRiYmJ5XIsp00cMMYoLi5OderU0bRp0xzb09PTFRwcLEmaNm2atm3bpg8++MAxcWDr1q2O\niQM//fRTodY0Jg64p7NnpQ8+kBYtsgez/v2l4cOlvn3tg/QBAPBELjlxYPPmzXrvvfcUGRmpqKgo\nSfblNhYsWKBdu3bJZrOpSZMmmjVrliQpPDxcI0aMUHh4uPz8/DRjxgy6Oz3A4cPS669Lc+ZI3btL\n48cTzAAAuBJOa0lzFlrS3MPOndKrr0qrVklxcdJjj9lnYQIA4E2uJbcwNw7l6scf7V2ZgwbZV17e\nv1+aNo2ABgDA1SKkoVxkZkrPPit16WJfJmPfPul//1eqWdPqygAAcE+ENFwTY+yzNMPDpdRU+7pm\nEybYV+QHAABlxw3WUWYpKdIDD9iX1Hj/falrV6srAgDAc9CShjJZu1a6+Wapd2/7khoENAAAyhct\nabgq+fnSlCnSm2/auzm7dbO6IgAAPBMhDVfs9Gn7chpHj0rbtknF3BACAACUA7o7cUW++0666SYp\nNFRKTCSgAQDgbLSkoVTbtkkDBkh//7s0dqzV1QAA4B0IaSjR9u32gDZ7tjRwoNXVAADgPejuRLF2\n7rTfPeDf/yagAQBQ0QhpKNI330j9+kkzZ0qDB1tdDQAA3oeQhkK+/Vbq21f65z+loUOtrgYAAO9E\nSEMBe/ZIffpI06dLd95pdTUAAHgvmzHGWF3E1bDZbHKzkt3GmTP2ZTaeeUYaN87qagAAcH/XklsI\naZBkv1H6sGFSUJB9HBoAALh215JbWIIDkuxroKWmSgsWWF0JAACQCGmQ9Omn0rRp0tatUuXKVlcD\nAAAkJg54vcOHpbvvlt57T2rUyOpqAADARYQ0L3bhgn0G5+OPSz17Wl0NAAC4FBMHvNj//I+UliZ9\n/LFks1ldDQAAnoeJA7hqa9dK//d/9jsLENAAAHA9dHd6oXPnpIcesi+1UaOG1dUAAICi0N3phZ55\nRjp4kOU2AABwNro7ccW++UaaO9d+f04AAOC66O70Inl50u9+J02ebL+zAAAAcF2ENC/y5ptSlSrc\nlxMAAHfAmDQvceiQFBUlbdoktWpldTUAAHiHa8kttKR5AWPsa6I99hgBDQAAd8HEAS/wySfS3r3S\nhx9aXQkAALhSdHd6uJwcqXVr6V//knr1sroaAAC8C92dKNbcuVJYGAENAAB3Q0uaB8vKkpo3t9+b\nMzra6moAAPA+tKShSDNmSDfdREADAMAd0ZLmoc6csbeibdggRURYXQ0AAN6JljQU8o9/SH36ENAA\nAHBXtKR5oOPHpZYtpW3bpKZNra4GAADvdS25hZDmgZ58Ujp/3n4bKAAAYB1CGhwOH5YiI6Xvv5eC\ng62uBgAA70ZIg8Pvfy/VrCm9/LLVlQAAAEIaJEkHDtiX3EhOlmrXtroaAADA7E5IkqZPl+6/n4AG\nAIAnoCXNQ5w8Kd1wg7R7txQaanU1AABAoiUNkt56S7r9dgIaAACegpY0D5CdbV8PbcUKKSrK6moA\nAMBFtKR5ucWLpRYtCGgAAHgSQpqbM0Z69VX7ArYAAMBzOC2kHTp0SN27d1ebNm0UERGh119/XZJ0\n4sQJxcbGqkWLFurdu7dOnTrleM+UKVPUvHlztWrVSmvXrnVWaR4lMdF+d4F+/ayuBAAAlCenjUnL\nyMhQRkaG2rVrp8zMTHXo0EFLly7V3LlzVbduXT399NN6+eWXdfLkSU2dOlVJSUkaPXq0tm3bptTU\nVPXq1UvJycny8SmYIxmTVtCAAdLgwdLvfmd1JQAA4HIuOSatQYMGateunSQpICBArVu3VmpqqpYv\nX664uDhJUlxcnJYuXSpJWrZsmUaNGiV/f3+FhYWpWbNm2rp1q7PK8wg//GC/ifo991hdCQAAKG8V\nMiYtJSVFO3fu1M0336wjR44oKChIkhQUFKQjR45IktLS0hR6yfoRoaGhSk1NrYjy3Na0adLDD0tV\nqlhdCQAAKG9+zj5BZmamhg0bptdee02BgYEFXrPZbLLZbMW+t6TXvN3Ro9KHH0o//mh1JQAAwBmc\nGtJycnI0bNgwjRkzRkOGDJFkbz3LyMhQgwYNlJ6ervr160uSQkJCdOjQIcd7Dx8+rJCQkCKPGx8f\n73gcExOjmJgYp30GVzVzpjR8uPTfrw8AALiAxMREJSYmlsuxnDZxwBijuLg41alTR9OmTXNsf/rp\np1WnTh0988wzmjp1qk6dOlVg4sDWrVsdEwd++umnQq1pTByQcnOl66+X1qyRIiKsrgYAABTnWnKL\n01rSNm/erPfee0+RkZGK+u8qq1OmTNGzzz6rESNGaM6cOQoLC9PixYslSeHh4RoxYoTCw8Pl5+en\nGTNm0N1ZjFWr7CGNgAYAgOfitlBuaOBAaehQ6b77rK4EAACU5FpyCyHNzRw+LEVGSocOSdWqWV0N\nAAAoiUuukwbnmDdPuusuAhoAAJ6OljQ3kp8v3XCD9NFHUvv2VlcDAABKQ0ual1i/XqpVi4AGAIA3\nIKS5kbfe4h6dAAB4C7o73cTRo1LLllJKilSjhtXVAACAK0F3pxd45x1pyBACGgAA3sLp9+7EtTNG\nmj1bevttqysBAAAVhZY0N/D555Kvr9Spk9WVAACAilKmkOaNY8Ks9NZb0gMPSNwlCwAA71FsSJs/\nf36R23NycjR69GinFYSCTp6UVqyQxoyxuhIAAFCRig1p06dP16xZswpsy8zM1O23364qVao4vTDY\nLVok9ekj1a1rdSUAAKAiFRvSNmzYoNmzZ+u1116TJB07dkzdu3dX+/bt9TYj2CvM++9L99xjdRUA\nAKCilbhO2unTp9W/f3917dpVS5cu1UMPPaQnnniiIusrxJvWSfv5Z6ljRyk1VapUyepqAADA1bqW\n3FJsSPvoo49ks9l09uxZTZgwQT169NDIkSMdJxw6dGjZK74G3hTSpkyRDh6UZs60uhIAAFAWTglp\n9957r2z/nU5ojHE8vmju3LllOuG18paQZozUtq30r39Jt95qdTUAAKAsnBLSXJW3hLRvvpEGD5b2\n75d8WM0OAAC3xG2hPND770ujRxPQAADwVrSkuaD8fOn666WEBKlNG6urAQAAZUVLmof5z3+kOnUI\naAAAeLMrusH65s2blZKSotzcXEn2VDh27FinFubNPvjA3tUJAAC8V6kh7Z577tH+/fvVrl07+fr6\nOrYT0pzjwgXpo4+kXbusrgQAAFip1JC2fft2JSUlFVqCA86xerV96Y1GjayuBAAAWKnUMWkRERFK\nT0+viFog+6zOu++2ugoAAGC1Umd3xsTEaNeuXYqOjlblypXtb7LZtHz58gop8HKePLvz9GmpcWMp\nJUWqVcvqagAAwLW6ltxSandnfHx8mQ6Mq/fxx1KPHgQ0AADAOmkupVcv6aGHpDvvtLoSAABQHpxy\nW6guXbpo8+bNCggIKDRpwGaz6cyZM2U64bXy1JB29KjUooWUni5VqWJ1NQAAoDxw704P8O9/Sxs3\nSgsWWF0JAAAoL9xxwAMsWUI3JwAA+A0taS7gl1+kpk3tXZ1Vq1pdDQAAKC+0pLm5Zcuk3r0JaAAA\n4DdXFNJSUlK0fv16SdK5c+csmzTgqejqBAAAlys1pP373//W8OHD9fvf/16SdPjwYd1xxx1OL8xb\nnDwpbd4s9e9vdSUAAMCVlBrS3nzzTW3atEnVq1eXJLVo0UJHjx51emHeYsUK+wK2gYFWVwIAAFxJ\nqSGtcuXKjttBSVJubi43Wy9HdHUCAICilBrSunXrpr/+9a86d+6c1q1bp+HDh2vgwIEVUZvHO3NG\nSkyUBgywuhIAAOBqSl2CIz8/X7Nnz9batWslSX369NEDDzxgWWuaJy3B8cEH9sVrV6ywuhIAAOAM\nTrvjQG5uriIiIrRnz54yF1fePCmkDR0qDR4sxcVZXQkAAHAGp62T5ufnp5YtW+rnn38u08FRvMxM\nacMGadAgqysBAACuyK+0HU6cOKE2bdooOjpa1apVk2RPhcuXL3d6cZ5s1Sqpc2epVi2rKwEAAK6o\n1JD2l7/8pSLq8DrM6gQAACXh3p0WOHdOCg6W9u2T6ta1uhoAAOAs15JbSm1JCwgIcMzkzM7OVk5O\njgICArg11DVISJCiowloAACgeKWGtMzMTMfj/Px8LV++XFu2bHFqUZ7uk08k7qwFAABKUqbuznbt\n2mnXrl3OqKdU7t7dmZMjNWgg7d4thYRYXQ0AAHAmp3Z3fvTRR47H+fn52r59u6pUqVKmk0HatElq\n2pSABgAASlbqbaFWrFihlStXauXKlVq7dq0CAwO1bNmyUg88btw4BQUFqW3bto5t8fHxCg0NVVRU\nlKKiorR69WrHa1OmTFHz5s3VqlUrx90NPNGyZfYFbAEAAEpSanfnpk2bdOuttxbYtnnzZnXp0qXE\nA3/++eekeg4KAAAdMUlEQVQKCAjQ2LFj9e2330qSJk2apMDAQE2YMKHAvklJSRo9erS2bdum1NRU\n9erVS8nJyfLxKZwh3bm70xipSRNp5UopIsLqagAAgLM57Y4DkvTYY48V2vboo4+WeuCuXbuqVhEr\ntRZV6LJlyzRq1Cj5+/srLCxMzZo109atW0s9h7vZvVvy9ZXatLG6EgAA4OqKHZP25Zdf6osvvtDR\no0f1j3/8wxGuzp49q/z8/DKf8I033tA777yjjh076tVXX1XNmjWVlpamW265xbFPaGioUlNTy3wO\nV3Wxq9Oie9MDAAA3UmxIy87O1tmzZ5WXl6ezZ886tlevXl1Lliwp08kefvhh/elPf5IkTZw4UU8+\n+aTmzJlT5L62EpJMfHy843FMTIxiYmLKVE9FW7pUmjbN6ioAAICzJCYmKjExsVyOVeqYtJSUFIWF\nhZXp4CkpKRo4cKBjTFpxr02dOlWS9Oyzz0qS+vbtq0mTJunmm28uXLCbjkk7eFBq317KyJD8Sp1T\nCwAAPIFTl+CoWrWqnnrqKSUlJSkrK8txwk8//fSqT5aenq7g4GBJ0ieffOKY+Tlo0CCNHj1aEyZM\nUGpqqvbu3avo6OirPr4rW75cGjCAgAYAAK5MqZHh7rvv1l133aWVK1dq1qxZmjdvnurVq1fqgUeN\nGqXPPvtMx48fV6NGjTRp0iQlJiZq165dstlsatKkiWbNmiVJCg8P14gRIxQeHi4/Pz/NmDGjxO5O\nd7RsmfTII1ZXAQAA3EWp3Z3t27fXjh07FBkZqd27d0uSOnbsqK+//rpCCrycO3Z3njolNW4spadL\n1apZXQ0AAKgoTu3urFSpkiSpQYMGWrlypRo2bKiTJ0+W6WTeatUqqVs3AhoAALhypYa0F154QadO\nndKrr76q8ePH68yZM5rGFMWrsmyZNGSI1VUAAAB3UmJ3Z15enl577bVCdwiwkrt1d164IAUFST/+\naP8NAAC8h9PuOODr66sFCxaU6cCw27jRfocBAhoAALgapXZ33nrrrXr00Ud11113qVq1ajLGyGaz\nqX379hVRn9vjhuoAAKAsSp3dGRMTU+RyGBs3bnRaUSVxp+7O/HypUSN7a1qLFlZXAwAAKtq15JZS\nQ5qrcaeQtn27NHq0fTwaAADwPk4bkyZJGRkZuv/++9W3b19JUlJSUrH320RBK1dKAwdaXQUAAHBH\npYa0e++9V71791ZaWpokqXnz5izBcYVWrLDfCgoAAOBqlRrSjh8/rrvuuku+vr6SJH9/f/lxA8pS\npaVJ+/dLXbpYXQkAAHBHpYa0gIAA/fLLL47nW7ZsUY0aNZxalCdYtUrq00fy97e6EgAA4I5KbRJ7\n9dVXNXDgQO3fv1+dO3fWsWPHtGTJkoqoza2tXCkNH251FQAAwF1d0ezO3Nxc/fjjjzLGqGXLlvK3\nsHnIHWZ3nj9vX7x2/36pTh2rqwEAAFZx6g3Ws7KyNGPGDG3atEk2m01du3bVww8/rOuuu65MJ/QG\niYlSZCQBDQAAlF2pLWnDhw9X9erVdc8998gYow8++ECnT5/Whx9+WFE1FuAOLWn/8z9S48bSM89Y\nXQkAALCSUxezDQ8PV1JSUqnbKoqrhzRjpLAwafVqKTzc6moAAICVnLqYbfv27fXll186nm/ZskUd\nOnQo08m8wXffSb6+UuvWVlcCAADcWaktaa1atVJycrIaNWokm82mgwcPqmXLlvLz85PNZtPu3bsr\nqlZJrt+SNmWKlJ4uvf661ZUAAACrOXXiQEJCQpkO7K1WrpRefNHqKgAAgLu7oiU4Tp48qUOHDik3\nN9exrX379k4trDiu3JJ2/LjUrJl05IhUubLV1QAAAKs5tSVt4sSJmjdvnpo2bSofn9+GsG3cuLFM\nJ/Rkq1dLPXsS0AAAwLUrNaQtWrRI+/btU6VKlSqiHrfGDdUBAEB5KXV2Z5s2bXTy5MmKqMWtZWdL\n69ZJ/ftbXQkAAPAEpbakPf/884qKilJERIQq/7cfz2azafny5U4vzp1s2iS1aGG/HRQAAMC1KjWk\njR07Vs8++6wiIiIcY9JsNpvTC3M3K1fS1QkAAMpPqbM7b7rpJm3btq2i6imVq87ubNVK+uADyaJJ\nrwAAwAU59bZQEyZMUOXKlTVo0CBHd6fEEhyX2rdPuvVWKTVV8il1lB8AAPAWTl2CY8eOHbLZbNqy\nZUuB7SzB8ZtVq+wTBghoAACgvFzRYrauxBVb0vr1kx54QBo2zOpKAACAK3HqDdYzMjJ0//33q2/f\nvpKkpKQkzZkzp0wn80TnzkmbN0u9elldCQAA8CSlhrR7771XvXv3VlpamiSpefPmmjZtmtMLcxcb\nN0odOkg1alhdCQAA8CTFhrSL9+k8fvy47rrrLvn6+kqS/P395edX6lA2r3FxPBoAAEB5KjakRUdH\nS5ICAgJ0/Phxx/YtW7aoBs1GkiRjCGkAAMA5im0SuzjI7dVXX9XgwYO1f/9+de7cWceOHdOSJUsq\nrEBX9sMPUn6+FB5udSUAAMDTFDu7MzQ0VBMmTJAxRsYYXbhwQcYYVa5cWb6+vpowYUJF1yrJtWZ3\nvvKKfY20mTOtrgQAALgip6yTlpeXp7Nnzxbafu7cuTKdyBOtWiVZlFUBAICHK7YlLSoqSjt37qzo\nekrlKi1pZ85IoaFSerpUrZrV1QAAAFfk1HXSULT166XOnQloAADAOYoNaevXr6/IOtwOszoBAIAz\ncVuoMjBGCgmR/vMfqVkzS0sBAAAujO7OCrZrlxQQQEADAADOQ0grA7o6AQCAsxHSyoCQBgAAnI0x\naVfpl1+kpk2lo0elypUtKwMAALgBxqRVoLVrpW7dCGgAAMC5CGlXafVqqV8/q6sAAACeju7Oq5Cf\nLwUHS199JYWFWVICAABwIy7Z3Tlu3DgFBQWpbdu2jm0nTpxQbGysWrRood69e+vUqVOO16ZMmaLm\nzZurVatWWrt2rbPKuiY7dki1axPQAACA8zktpN13331KSEgosG3q1KmKjY1VcnKyevbsqalTp0qS\nkpKStGjRIiUlJSkhIUGPPPKI8vPznVVamdHVCQAAKorTQlrXrl1Vq1atAtuWL1+uuLg4SVJcXJyW\nLl0qSVq2bJlGjRolf39/hYWFqVmzZtq6dauzSiszQhoAAKgoFTpx4MiRIwoKCpIkBQUF6ciRI5Kk\ntLQ0hYaGOvYLDQ1VampqRZZWqhMnpO++k7p2tboSAADgDfysOrHNZpPNZivx9eLEx8c7HsfExCgm\nJqYcKyva2rXSbbdJ113n9FMBAAA3lZiYqMTExHI5VoWGtKCgIGVkZKhBgwZKT09X/fr1JUkhISE6\ndOiQY7/Dhw8rJCSk2ONcGtIqCl2dAACgNJc3Hk2aNKnMx6rQ7s5BgwZp/vz5kqT58+dryJAhju0L\nFy5Udna2Dhw4oL179yo6OroiSytRfr6UkEBIAwAAFcdpLWmjRo3SZ599puPHj6tRo0b685//rGef\nfVYjRozQnDlzFBYWpsWLF0uSwsPDNWLECIWHh8vPz08zZswosbuzou3cKdWsab8dFAAAQEVgMdsr\n8NJL0vHj0vTpFXpaAADg5lxyMVtPwng0AABQ0WhJK8WJE/Y7DBw9ysxOAABwdWhJc6J16+xroxHQ\nAABARSKklYKuTgAAYAW6O0uQny81bCht3izdcEOFnBIAAHgQujudZNcuqXp1AhoAAKh4hLQS0NUJ\nAACsQkgrQUKC1Lev1VUAAABvxJi0Ypw6JTVqZF96o0oVp58OAAB4IMakOcH69dKttxLQAACANQhp\nxaCrEwAAWInuziIYY+/q3LBBatnSqacCAAAejO7Ocvbdd1KlSlKLFlZXAgAAvBUhrQgXuzptNqsr\nAQAA3oqQVgTGowEAAKsxJu0ymZlScLCUni4FBDjtNAAAwAswJq0cffqpFB1NQAMAANYipF2Grk4A\nAOAKCGmXMIb7dQIAANdASLvE3r1STo7Upo3VlQAAAG9HSLvE6tUsvQEAAFwDIe0SCQl0dQIAANfA\nEhz/lZUl1a8vHTok1axZ7ocHAABeiCU4ysFnn0nt2hHQAACAayCk/RddnQAAwJUQ0v6L9dEAAIAr\nIaRJOnBAOnnS3t0JAADgCghpktaskfr0kXz4NgAAgIsglojxaAAAwPV4/RIc2dlSvXrSvn1S3brl\ndlgAAACW4LgWmzdLrVoR0AAAgGvx+pDGrE4AAOCKCGmENAAA4IK8ekxaWpoUESEdPSr5+ZXLIQEA\nABwYk1ZGa9ZIsbEENAAA4Hq8OqTR1QkAAFyV13Z35uZK9etL330nNWxYDoUBAABchu7OMti2TWrU\niIAGAABck9eGtIQE+62gAAAAXJFXhzTGowEAAFfllWPSjh+XmjaVjh2TKlcup8IAAAAuw5i0q7Ru\nnRQTQ0ADAACuyytDGl2dAADA1Xldd2d+vn1G5xdf2Ls8AQAAnIXuzqvwzTdS9eoENAAA4Nq8LqTR\n1QkAANwBIQ0AAMAFWTImLSwsTNWrV5evr6/8/f21detWnThxQnfddZd+/vlnhYWFafHixapZs2bh\ngq+hb/fMGSkkRMrIkKpVu9ZPAQAAUDK3G5Nms9mUmJionTt3auvWrZKkqVOnKjY2VsnJyerZs6em\nTp1a7ufdsEHq1ImABgAAXJ9l3Z2Xp8rly5crLi5OkhQXF6elS5eW+znp6gQAAO7Cspa0Xr16qWPH\njnrrrbckSUeOHFFQUJAkKSgoSEeOHCnXcxpDSAMAAO7Dz4qTbt68WcHBwTp27JhiY2PVqlWrAq/b\nbDbZbLZi3x8fH+94HBMTo5iYmFLP+eOP9qDWunVZqwYAAChZYmKiEhMTy+VYli9mO2nSJAUEBOit\nt95SYmKiGjRooPT0dHXv3l179uwptH9ZB+BNny4lJUn//nd5VA0AAFA6t5o4cO7cOZ09e1aS9Ouv\nv2rt2rVq27atBg0apPnz50uS5s+fryFDhpTreenqBAAA7qTCW9IOHDigO+64Q5KUm5uru+++W889\n95xOnDihESNG6ODBg+W+BEdWllS/vnT4sFSjRrl8DAAAgFJdS0ua5d2dV6ssHzYhQZo8WfrPf5xU\nFAAAQBHcqrvTCgkJUp8+VlcBAABw5bwmpDEeDQAAuBOPD2kpKdKJE1JUlNWVAAAAXDmPD2lr1ti7\nOn08/pMCAABP4vHRha5OAADgjjx6dmd2tn3pjb17pXr1nFwYAADAZZjdWYwvv5SaNyegAQAA9+PR\nIY2uTgAA4K4IaQAAAC7IY8ekZWRIrVtLx45Jfn4VUBgAAMBlGJNWhLVrpV69CGgAAMA9eWxI41ZQ\nAADAnXlkd2denhQUJO3cKTVqVEGFAQAAXIbuzsts3y41aEBAAwAA7ssjQxqzOgEAgLsjpAEAALgg\njxuTduKEFBYmHT0qXXddxdUFAABwOcakXWL9eum22whoAADAvXlcSFuzhqU3AACA+/Oo7k5jpNBQ\nKTHRfmN1AAAAK9Hd+V/ffWfv5mzWzOpKAAAAro1HhbSLszptNqsrAQAAuDYeGdIAAADcnceMScvM\nlIKDpfR0KSDAgsIAAAAuw5g0SRs3SjfdREADAACewWNCGl2dAADAkxDSAAAAXJBHhLSffpKysqS2\nba2uBAAAoHx4REhj6Q0AAOBpPCqkAQAAeAq3X4LjwgWpXj0pJUWqXdu6ugAAAC7n1UtwfP651KYN\nAQ0AAHgWtw9pdHUCAABPREgDAABwQW4d0g4fljIypI4dra4EAACgfLl1SFuzRoqNlXx9ra4EAACg\nfLl1SKOrEwAAeCq3XYIjN9e+9EZSkhQcbHVVAAAAhXnlEhxffSVdfz0BDQAAeCa3DWl0dQIAAE/m\ntiFtzRpCGgAA8FxuOSbt6FGjZs2kY8ekSpWsrggAAKBoXjcmbd06qXt3AhoAAPBcbhnSGI8GAAA8\nnVt2d9avb7Rli9SkidXVAAAAFM/rujtr1iSgAQAAz+ZSIS0hIUGtWrVS8+bN9fLLLxe7H12dAADA\n07lMSMvLy9Ojjz6qhIQEJSUlacGCBfrhhx+K3JeQ5p4SExOtLgHXgOvn3rh+7otr571cJqRt3bpV\nzZo1U1hYmPz9/TVy5EgtW7asyH27davg4lAu+IvGvXH93BvXz31x7byXy4S01NRUNWrUyPE8NDRU\nqampRe5btWpFVQUAAGANlwlpNpvN6hIAAABchssswbFlyxbFx8crISFBkjRlyhT5+PjomWeeKbAf\nYQ4AALiTskYtlwlpubm5atmypTZs2KCGDRsqOjpaCxYsUOvWra0uDQAAoML5WV3ARX5+fvrnP/+p\nPn36KC8vT/fffz8BDQAAeC2XGZMmSf369dOPP/6on376Sc8991yB1650DTVY79ChQ+revbvatGmj\niIgIvf7665KkEydOKDY2Vi1atFDv3r116tQpiytFSfLy8hQVFaWBAwdK4vq5k1OnTunOO+9U69at\nFR4erq+++orr5yamTJmiNm3aqG3btho9erQuXLjAtXNh48aNU1BQkNq2bevYVtL1mjJlipo3b65W\nrVpp7dq1pR7fpUJaca5mDTVYz9/fX9OmTdP333+vLVu26M0339QPP/ygqVOnKjY2VsnJyerZs6em\nTp1qdakowWuvvabw8HDHOFCun/t4/PHH1b9/f/3www/avXu3WrVqxfVzAykpKXrrrbe0Y8cOffvt\nt8rLy9PChQu5di7svvvuc4ylv6i465WUlKRFixYpKSlJCQkJeuSRR5Sfn1/yCYwb+OKLL0yfPn0c\nz6dMmWKmTJliYUW4GoMHDzbr1q0zLVu2NBkZGcYYY9LT003Lli0trgzFOXTokOnZs6f59NNPzYAB\nA4wxhuvnJk6dOmWaNGlSaDvXz/X98ssvpkWLFubEiRMmJyfHDBgwwKxdu5Zr5+IOHDhgIiIiHM+L\nu16TJ082U6dOdezXp08f8+WXX5Z4bLdoSbuaNdTgWlJSUrRz507dfPPNOnLkiIKCgiRJQUFBOnLk\niMXVoTh/+MMf9Pe//10+Pr/9FcH1cw8HDhxQvXr1dN9996l9+/b63e9+p19//ZXr5wZq166tJ598\nUo0bN1bDhg1Vs2ZNxcbGcu3cTHHXKy0tTaGhoY79riTLuEVIY9kN95SZmalhw4bptddeU2BgYIHX\nbDYb19VFrVy5UvXr11dUVFSx08a5fq4rNzdXO3bs0COPPKIdO3aoWrVqhbrHuH6uad++fZo+fbpS\nUlKUlpamzMxMvffeewX24dq5l9KuV2nX0i1CWkhIiA4dOuR4fujQoQJpFK4nJydHw4YN05gxYzRk\nyBBJ9v9RZGRkSJLS09NVv359K0tEMb744gstX75cTZo00ahRo/Tpp59qzJgxXD83ERoaqtDQUN10\n002SpDvvvFM7duxQgwYNuH4u7uuvv1bnzp1Vp04d+fn5aejQofryyy+5dm6muL8rL88yhw8fVkhI\nSInHcouQ1rFjR+3du1cpKSnKzs7WokWLNGjQIKvLQjGMMbr//vsVHh6uJ554wrF90KBBmj9/viRp\n/vz5jvAG1zJ58mQdOnRIBw4c0MKFC9WjRw+9++67XD830aBBAzVq1EjJycmSpPXr16tNmzYaOHAg\n18/FtWrVSlu2bFFWVpaMMVq/fr3Cw8O5dm6muL8rBw0apIULFyo7O1sHDhzQ3r17FR0dXfLBynsA\nnbOsWrXKtGjRwtxwww1m8uTJVpeDEnz++efGZrOZG2+80bRr1860a9fOrF692vzyyy+mZ8+epnnz\n5iY2NtacPHnS6lJRisTERDNw4EBjjOH6uZFdu3aZjh07msjISHPHHXeYU6dOcf3cxMsvv2zCw8NN\nRESEGTt2rMnOzubaubCRI0ea4OBg4+/vb0JDQ83bb79d4vX661//am644QbTsmVLk5CQUOrxXeaO\nAwAAAPiNW3R3AgAAeBtCGgAAgAsipAEAALggQhoAAIALIqQBKFe//PKLoqKiFBUVpeDgYIWGhioq\nKkqBgYF69NFHK6SGb775RqtXry63491+++06c+ZMuR0PAK4EszsBOM2kSZMUGBioCRMmVOh5582b\np+3bt+uNN96o0PM6U25urvz8/KwuA0AFoiUNgFNd/H9gYmKiBg4cKEmKj49XXFycbrvtNoWFhenj\njz/WU089pcjISPXr10+5ubmSpO3btysmJkYdO3ZU3759Hat4X+rDDz9U27Zt1a5dO8XExCgnJ0d/\n+tOftGjRIkVFRenDDz/Ur7/+qnHjxunmm29W+/bttXz5ckn2MDd48GB1795dLVq00J///OciP0NY\nWJhOnDihlJQUtW7dWg8++KAiIiLUp08fnT9/vsC+Z8+eVdOmTR2f4cyZM2ratKny8vK0b98+9evX\nTx07dtRtt92mH3/8UZK0YsUK3XLLLWrfvr1iY2N19OhRx/c0ZswY3XrrrYqLi9P333+v6OhoRUVF\n6cYbb9RPP/10rZcHgCtz0vpuAGDi4+PNK6+8YowxZuPGjWbAgAHGGGNefPFF07VrV5Obm2u++eYb\nU6VKFcfCjnfccYdZunSpyc7ONp06dTLHjx83xhizcOFCM27cuELnaNu2rUlLSzPGGHP69GljjDHz\n5s0z48ePd+zz3HPPmffee88YY8zJkydNixYtzK+//mrmzp1rgoODzYkTJ0xWVpaJiIgwX3/9daFz\nhIWFmV9++cUcOHDA+Pn5mW+++cYYY8yIESMcx73UfffdZ5YuXWqMMWbWrFnmqaeeMsYY06NHD7N3\n715jjDFbtmwxPXr0cNR00VtvvWWefPJJx/fUsWNHc/78eWOMMePHjzfvv/++McaYnJwck5WVVfyX\nD8Dt0XYOoMLZbDb169dPvr6+ioiIUH5+vvr06SNJatu2rVJSUpScnKzvv/9evXr1kiTl5eWpYcOG\nhY7VpUsXxcXFacSIERo6dKgke+uduWQkx9q1a7VixQq98sorkqQLFy7o4MGDstls6t27t2rVqiVJ\nGjp0qDZt2qQOHToUW3uTJk0UGRkpSerQoYNSUlIK7fPAAw/ob3/7mwYPHqx58+Zp9uzZyszM1Bdf\nfKHhw4c79svOzpZkvx/xiBEjlJGRoezsbDVt2tTxPQ0aNEiVK1eWJHXq1El//etfdfjwYQ0dOlTN\nmjW7gm8bgLsipAGwRKVKlSRJPj4+8vf3d2z38fFRbm6ujDFq06aNvvjiixKPM3PmTG3dulX/93//\npw4dOmj79u1F7vfxxx+refPmBbZ99dVXBZ4bY+TjU/IokIuBSZJ8fX2VlZVVaJ/OnTsrJSVFiYmJ\nysvLU3h4uM6cOaNatWpp586dhfYfP368nnrqKQ0YMECfffaZ4uPjHa9VrVrV8XjUqFG65ZZbtHLl\nSvXv31+zZs1S9+7dS6wXgPtiTBqACmeuYL5Sy5YtdezYMW3ZskWSlJOTo6SkpEL77du3T9HR0Zo0\naZLq1aunw4cPq3r16jp79qxjnz59+uj11193PL8YlIwxWrdunU6ePKmsrCwtW7ZMXbp0udaPJ0ka\nO3as7r77bo0bN06SVL16dTVp0kRLlixxnHv37t2S7OPWLrYSzps3z3GMy7+nAwcOqEmTJho/frwG\nDx6sb7/9tlxqBeCaCGkAnMpmszl+F/X40n0ufe7v768lS5bomWeeUbt27RQVFaUvv/yy0PGffvpp\nRUZGqm3bturSpYsiIyPVvXt3JSUlOSYOTJw4UTk5OYqMjFRERIRefPFFx3mio6M1bNgw3Xjjjbrz\nzjvVvn37Yj9DcbUWZfTo0Tp58qRGjRrl2Pb+++9rzpw5ateunSIiIhwTGOLj4zV8+HB17NhR9erV\nK/Z7Wrx4sSIiIhQVFaXvv/9eY8eOLfLcADwDS3AA8FrOXKpjyZIlWrFihebPn1/uxwbgHRiTBsBr\nXd5SVV7Gjx+vNWvWaNWqVeV+bADeg5Y0AAAAF8SYNAAAABdESAMAAHBBhDQAAAAXREgDAABwQYQ0\nAAAAF0RIAwAAcEGENAAAABf0/4eDOyHyusihAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x105880c50>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Visibly, the temperature converges towards the \"bare rock\" value of $255 K$ while the outoing heat flux approaches the incoming value of $236.25 W/m^2$.\n",
      "\n",
      "We can also inspect some of the variables as needed:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "L_out[:10]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "array([  0.00e+00,   5.71e-03,   9.13e-02,   4.62e-01,   1.46e+00,\n",
        "         3.54e+00,   7.28e+00,   1.33e+01,   2.21e+01,   3.42e+01])"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i in range(10):\n",
      "    print \"%2.2f\\t%2.2e\\t%.2f\" % (T[i], hc[i], L_out[i])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.00\t0.00e+00\t0.00\n",
        "17.81\t7.46e+09\t0.01\n",
        "35.62\t1.49e+10\t0.09\n",
        "53.43\t2.24e+10\t0.46\n",
        "71.21\t2.98e+10\t1.46\n",
        "88.91\t3.72e+10\t3.54\n",
        "106.46\t4.46e+10\t7.28\n",
        "123.72\t5.18e+10\t13.28\n",
        "140.53\t5.88e+10\t22.11\n",
        "156.68\t6.56e+10\t34.17\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "L_out[1000]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "236.24999999999986"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "L_solar_in"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "236.24999999999997"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}