MIET2039_CL07-Modelling-1D-2D-Conduction-SOLN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4 style=\"color:#7171C6;\">Exercise 7.1 steady no heat generation</h4>\n",
    "\n",
    "<img src=\"img/conduction-numerical-steady-fixedT.png\" width=\"300px\"/>\n",
    "\n",
    "A plane wall of length 2m has a temperature of 60$^o$C at the left boundary and  20$^o$C at the right boundary. The conductivity of the material is k = 28W/mK.\n",
    "\n",
    "The solution is simple, since Fourier's conduction equation for a constant conductivity in this case is a linear gradient. Applying this to the nodes the solution temperature distribution becomes: <br>\n",
    "[60 50 40 30 20]\n",
    "\n",
    "We will solve the problem numerically using the finite difference discretization discussed in class.\n",
    "\n",
    "Left boundary condition\n",
    "$$T_0 = T_{\\textrm{left}}$$\n",
    "Right boundary condition\n",
    "$$T_4 = T_{\\textrm{right}}$$\n",
    "Internal nodes\n",
    "$$\\frac{T_{i-1} - 2T_i+T_{i+1}}{\\Delta x^2}=0 $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 60.  35.  10.  15.  20.]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7b70e48>]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAElRJREFUeJzt3X+snXV9wPH3h1s1VjZSN9proCJiRLZsQ4x1Cz96pK2i\nWyzZFjZdUFzUZXHBZMtCQUlvoktETRbn9ofN1FS0megyZUND21uvdUt0jIISLB1mkSn2XqOsDGR/\nkPrZH8+5cHo9997z+zw/3q/khHNPn3PPNycPH56+z/McIjORJNXDWdNegCRpdBzqklQjDnVJqhGH\nuiTViENdkmrEoS5JNdLTUI+IcyLi8xFxPCIejIjXRMSmiDgYESci4u6IOGfci5Ukra3XI/WPAl/O\nzEuA3wAeAvYAhzPzYuAIcPN4lihJ6lWsd/FRRPwicF9mXrTi8YeA7Zm5FBGzwEJmvmJ8S5UkraeX\nI/ULgR9HxKci4lhE7IuIjcCWzFwCyMxFYPM4FypJWl8vQ30DcBnwd5l5GfBTivSy8hDf7xuQpCnb\n0MM2PwC+n5n/0f75HymG+lJEbOnILz/q9uSIcNhL0gAyM/p9zrpH6u3E8v2IeHn7oR3Ag8CdwA3t\nx94GfGm13/GiFyUnTyaZ3oa97d27d+prqMvN99L3s8y3QfVypA5wI/DZiHgO8F/A24EZ4I6I+GPg\nEeC61Z78rnfBW94Chw7BzMzAa5UkraOnUxoz81uZ+erMvDQzfzczH8/MxzJzZ2ZenJmvy8xTqz3/\n1luLf77//aNZtCSpu4lcUTozAwcOwL59MD8/iVesr1arNe0l1Ibv5Wj5fpbDuuepD/0CEbn8GvPz\ncP31cOwYzM6O9WUlqdIighzHB6WjtGPHs3399OlJvrIkNcPEv9DLvi5J4zPR/LJscREuuwxuv704\nepcknakS+WXZ7Gwx0K+/vhjwkqTRmNr3qdvXJWn0pvo/ybCvS9JoTaWpd7KvS9LPq1RT72Rfl6TR\nmfpQB/u6JI1KKYY62NclaRSm3tQ72dclqVDZpt7Jvi5JwynVUAf7uiQNo3RDHezrkjSoUjX1TvZ1\nSU1Wi6beyb4uSf0r7VAH+7ok9avUQx3s65LUj9I29U72dUlNU7um3sm+Lkm9qcRQB/u6JPWiMkMd\n7OuStJ5KNPVO9nVJTVDrpt7Jvi5Jq6vcUAf7uiStppJDHezrktRN5Zp6J/u6pLpqTFPvZF+XpDNV\neqiDfV2SOlV+qIN9XZKWVbqpd7KvS6qTRjb1TvZ1SarRUAf7uiTVaqiDfV1Ss/XU1CPie8DjwM+A\npzNzW0RsAj4HXAB8D7guMx/v8tyJNPVO9nVJVTfupv4zoJWZr8zMbe3H9gCHM/Ni4Ahwc78vPi72\ndUlN1etQjy7b7gb2t+/vB64d1aJGwb4uqYl6HeoJHIqIeyLiHe3HtmTmEkBmLgKbx7HAYdjXJTXN\nhh63uzwzT0bEucDBiDhBMeg7rRrO5+bmnrnfarVotVp9LnMwMzNw4EDR16+80r4uqbwWFhZYWFgY\n+vf0ffFRROwFngTeQdHZlyJiFvhqZl7SZfuJf1C60vx80dePHSt6uySV3dg+KI2IjRFxdvv+C4DX\nAQ8AdwI3tDd7G/Clfl98Uuzrkppi3SP1iLgQ+CeKvLIB+GxmfjAiXgjcAWwFHqE4pfFUl+dP/Ugd\nimG+axdcdRV01CBJKqVBj9Rr890vvfD8dUlV0fjvfumF569LqrtGDXWwr0uqt8YNdfD8dUn11aim\n3sm+LqnMbOp9sq9LqqPGDnWwr0uqn0YPdbCvS6qXxjb1TvZ1SWVjUx+CfV1SXTjU2+zrkurAod7B\nvi6p6mzqK9jXJZWBTX1E7OuSqsyh3oV9XVJVOdRXYV+XVEU29TXY1yVNi019DOzrkqrGob4O+7qk\nKnGo98C+LqkqbOo9sq9LmiSb+pjZ1yVVgUO9D/Z1SWXnUO+TfV1SmdnUB2BflzRuNvUJsq9LKiuH\n+oDs65LKyKE+BPu6pLKxqQ/Jvi5pHGzqU2Jfl1QmDvURsK9LKguH+ojY1yWVgU19hOzrkkbFpl4C\n9nVJ0+ZQHzH7uqRpcqiPgX1d0rT0PNQj4qyIOBYRd7Z/3hQRByPiRETcHRHnjG+Z1TIzAwcOwL59\nMD8/7dVIapJ+jtTfA3yn4+c9wOHMvBg4Atw8yoVVnX1d0jT0NNQj4nzgjcDfdzy8G9jfvr8fuHa0\nS6s++7qkSev1SP2vgb8EOs9N3JKZSwCZuQhsHvHaasG+LmmS1h3qEfHbwFJm3g+sdc5kM05G75N9\nXdIkbehhm8uBN0XEG4HnA78QEbcDixGxJTOXImIW+NFqv2Bubu6Z+61Wi1arNdSiq6azrx87Vvws\nSZ0WFhZYWFgY+vf0dUVpRGwH/iIz3xQRHwJ+kpm3RcRNwKbM3NPlOY25onQ9c3Nw9CgcOlQcwUvS\naqZxRekHgV0RcQLY0f5Za7CvSxo3v/tlwvx+GEm98LtfKsLz1yWNk0N9Cjx/XdK4ONSnxL4uaRxs\n6lNkX5e0Gpt6BdnXJY2aQ33K7OuSRsmhXgL2dUmjYlMvCfu6pE429Yqzr0saBYd6idjXJQ3LoV4y\n9nVJw7Cpl5B9XZJNvUbs65IG5VAvKfu6pEE41EvMvi6pXzb1krOvS81kU68p+7qkfjjUK8C+LqlX\nDvWKsK9L6oVNvULs61Jz2NQbwL4uaT0O9Yqxr0tai0O9guzrklZjU68o+7pUbzb1hrGvS+rGoV5h\n9nVJKznUK86+LqmTTb0G7OtS/djUG8y+LmmZQ70m7OuSwKFeK/Z1STb1mrGvS/VgUxdgX5eazqFe\nQ/Z1qbkc6jVlX5eayaZeY/Z1qbrG1tQj4nkR8c2IuC8iHoiIve3HN0XEwYg4ERF3R8Q5gyxc42Nf\nl5qnpyP1iNiYmU9FxAzwb8CNwO8BP8nMD0XETcCmzNzT5bkeqU/Z3BwcPQqHDsHMzLRXI6kXYz37\nJTOfat99HrABSGA3sL/9+H7g2n5fXJNhX5eao6ehHhFnRcR9wCJwKDPvAbZk5hJAZi4Cm8e3TA1j\nZgYOHIB9+2B+ftqrkTROvR6p/ywzXwmcD2yLiF+lOFo/Y7NRL06jY1+XmmFDPxtn5v9GxAJwDbAU\nEVsycykiZoEfrfa8ubm5Z+63Wi1ardZAi9VwOs9ft69L5bKwsMDCwsLQv2fdD0oj4peBpzPz8Yh4\nPnA38EFgO/BYZt7mB6XVcfo07NoFV11VfIAqqZwG/aC0l6H+axQfhJ7Vvn0uM/8qIl4I3AFsBR4B\nrsvMU12e71AvGc9fl8pvbEN9WA71cpqfL/r6sWNFb5dULn6hl/ri98NI9eRQbzDPX5fqx/zScPZ1\nqZzMLxqI569L9eJQl31dqhGHugD7ulQXNnU9w74ulYdNXUOzr0vV51DXGezrUrU51PVz7OtSddnU\n1ZV9XZoum7pGyr4uVZNDXauyr0vV41DXmuzrUrXY1LUu+7o0eTZ1jY19XaoOh7p6Yl+XqsGhrp7Z\n16Xys6mrL/Z1aTJs6poI+7pUbg519c2+LpWXQ10Dsa9L5WRT18Ds69L42NQ1cfZ1qXwc6hqKfV0q\nF4e6hmZfl8rDpq6RsK9Lo2VT11TZ16VycKhrZOzr0vQ51DVS9nVpumzqGjn7ujQ8m7pKw74uTY9D\nXWNhX5emw6GusbGvS5NnU9dY2delwdjUVUr2dWmy1h3qEXF+RByJiAcj4oGIuLH9+KaIOBgRJyLi\n7og4Z/zLVRXZ16XJWTe/RMQsMJuZ90fE2cC9wG7g7cBPMvNDEXETsCkz93R5vvlFnD4Nu3bBVVfB\n3Ny0VyOV36D5pe+mHhFfBP62fduemUvtwb+Qma/osr1DXYB9XerHRJp6RLwEuBT4BrAlM5cAMnMR\n2Nzvi6tZ7OvS+p54Au66a/Dnb+h1w3Z6+QLwnsx8MiJWHn6vejg+1/H37VarRavV6m+Vqo3Ovn7o\nEMzMTHtF0nQ9/TTccw/s27fAkSMLLC7CeecN/vt6yi8RsQH4F+ArmfnR9mPHgVZHfvlqZl7S5bnm\nF53Bvq4my4SHHioOag4fhq99DV76Uti5s/j34oorYOPGMTf1iPg08OPM/POOx24DHsvM2/ygVP2y\nr6tJTp4sBvjybcOGYoDv2gVXXw3nnvvzzxnbUI+Iy4GjwAMUiSWBW4B/B+4AtgKPANdl5qkuz3eo\nq6v5+aKvHztW9HapLp54ojgCXx7iP/whvPa1xRDfuRMuughinXE9sbNf+n4Bh7rWMDcHR4/a11Vt\ny118Oancdx9s2/ZsUrnssv73b4e6Ksm+rirqtYsPw6GuyrKvqwoG6eLDcKir0uzrKpvVuvjy0Xgv\nXXwYDnVV3twcfP3rcPCgfV2Tt1YX37kTXvWqye6XDnVV3nJf374d9u6d9mpUd2t18Z074corh+/i\nw3CoqxaW+/pnPlN0SmmUVuviO3cWn+eMuosPw6Gu2pifh7e+Fe69176u4Uy7iw/Doa5asa9rEGXr\n4sNwqKtW7OvqRdm7+DAc6qod+7q6qVIXH4ZDXbVkX9daXXznTnjZy8rbxYfhUFdt2debpU5dfBgO\nddWWfb3elrv44cPFIF/Zxa+4Al7wgmmvcvIc6qo1+3q9nDxZpLXlo/HOLn711bDZ/zmmQ131Z1+v\nrieeePYrlpvUxYfhUFcj2NerYbmLLyeVpnbxYTjU1Qj29XLq1sUvvPDZpNLULj4Mh7oaw75eDnbx\n8XKoq1Hs65O3sos/+uiZ/99Nu/hoOdTVOPb18erWxV/96meHuF18vBzqahz7+mh5vni5ONTVSPb1\n4djFy8uhrsayr/fO88Wrw6GuRrOvd+f54tXlUFej2dcLdvH6cKir8Zra1+3i9eRQl2hGX7eLN4ND\nXWqrW1+3izeTQ11qq3pft4sLHOrSGarW1+3iWsmhLq1Q5r5uF9d6HOpSF2Xp63Zx9cuhLnUxrb6+\nsosfPVp8v/jyEL/ySti4cXLrUfU41KVVTKqvr9XFd+yAc88d32urfhzq0hrG0deXu/jy0XhnF9+1\nCy66yC6uwY1tqEfEJ4DfAZYy89fbj20CPgdcAHwPuC4zH1/l+Q51lcKwfb2zix8+/Oz3i9vFNQ6D\nDvWzetjmU8DrVzy2BzicmRcDR4Cb+31hDWZhYWHaS6isW28tWvcHPlD8vN57mQnHj8PHPga7dxf5\n5N3vLo7Q3/teWFqCI0fglluKDz2bPtDdN8th3aGemf8K/M+Kh3cD+9v39wPXjnhdWoX/4gxuZgYO\nHICPf7wYxt3ey5Mni/Z+ww2wdStccw1861vw5jfDww8XR+cf/jC8/vV+0LmS+2Y5bBjweZszcwkg\nMxcjwksjVAmzs3D77XD99cWgXquLv+99dnFVz6BDfSWjuSpjxw545zvhIx8pjtqXu/gnP2kXV/X1\ndPZLRFwA/HPHB6XHgVZmLkXELPDVzLxklec68CVpAIN8UNrrkXq0b8vuBG4AbgPeBnxplIuSJA2m\nl1MaDwAt4JeAJWAv8EXg88BW4BGKUxpPjXWlkqR1jf3iI0nS5PRynnpPIuKaiHgoIv4zIm5aZZu/\niYiHI+L+iLh0VK9dN+u9lxGxPSJORcSx9u1901hnFUTEJyJiKSK+vcY27pc9Wu/9dN/sT0ScHxFH\nIuLBiHggIm5cZbve99HMHPpG8R+H71JcYfoc4H7gFSu2eQNwV/v+a4BvjOK163br8b3cDtw57bVW\n4QZcAVwKfHuVP3e/HO376b7Z3/s5C1zavn82cGLY2TmqI/VtwMOZ+UhmPg38A8UFSp12A58GyMxv\nAudExJYRvX6d9PJewpkfXGsV2f3iuU7ul33o4f0E982eZeZiZt7fvv8kcBw4b8Vmfe2joxrq5wHf\n7/j5B10WtnKbR7tso97eS4Dfav9V7K6I+JXJLK2W3C9Hz31zABHxEoq/BX1zxR/1tY+O6uIjTda9\nwIsz86mIeAPF2Ugvn/KaJHDfHEhEnA18AXhP+4h9YKM6Un8UeHHHz+e3H1u5zdZ1tlEP72VmPpmZ\nT7XvfwV4TkS8cHJLrBX3yxFy3+xfRGygGOi3Z2a3a3762kdHNdTvAV4WERdExHOBP6S4QKnTncBb\nASLiN4FT2f7+GJ1h3feys6dFxDaKU1Mfm+wyK2XlxXOd3C/7t+r76b45kE8C38nMj67y533toyPJ\nL5l5OiL+DDhI8R+KT2Tm8Yj4k+KPc19mfjki3hgR3wV+Crx9FK9dN728l8DvR8SfAk8D/wf8wfRW\nXG6dF89FxH9TXDz3XNwvB7Le+4n7Zl8i4nLgj4AHIuI+iu/RuoXi7LeB9lEvPpKkGhnZxUeSpOlz\nqEtSjTjUJalGHOqSVCMOdUmqEYe6JNWIQ12SasShLkk18v9X6ykMK2SkXwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x681128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# ******************************************\n",
    "# 1. PRELOAD MODULES\n",
    "# ******************************************\n",
    "import numpy as np                \n",
    "from matplotlib import pyplot as plt   \n",
    "%matplotlib inline  \n",
    "# if you're using sagemathcloud the las line isn't needed\n",
    "\n",
    "\n",
    "# ******************************************\n",
    "# 2. Set up the geometry\n",
    "# ******************************************\n",
    "# length of material\n",
    "L = 2\n",
    "\n",
    "# number of nodes (0 to 4 inclusive)\n",
    "nx = 5\n",
    "\n",
    "# distance between nodes\n",
    "dx = L/(nx-1)\n",
    "\n",
    "\n",
    "# ******************************************\n",
    "# 3. Setup the nodes\n",
    "# ******************************************\n",
    "# create a 1D array of nodes where the nodes have a value of 10-degrees\n",
    "# Later test the effect of different initial temperature value instead of 10 \n",
    "\n",
    "T = np.ones(nx)*10\n",
    "\n",
    "# set the boundary condictions T_left = 60, T_right = 20\n",
    "T[0] = 60\n",
    "T[nx-1] = 20\n",
    "\n",
    "\n",
    "# ******************************************\n",
    "# 4. Solve the problem\n",
    "# ******************************************\n",
    "# write out the discretized equation using array slicing and array arithmetic \n",
    "# this is just like vector addition\n",
    "\n",
    "T[1:-1] = (T[2:] + T[0:-2])/2\n",
    "\n",
    "\n",
    "\n",
    "# ******************************************\n",
    "# 5. Post process data \n",
    "# ******************************************\n",
    "# print the array T, and then plot it using plt.plot(....)\n",
    "print(T)\n",
    "plt.plot(np.linspace(0,L,nx), T)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have only performed one iteration which shows that the solution hasn't reached a linear profile. We can say that the solution hasn't converged. We need to iterate continuously until convergence is achieved. \n",
    "\n",
    "We add a continuous loop to repeat the iterations and then introduce a convergence criteria based on residuals. We can calculate the residuals based on the change in values between the new and previous iteration values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 ) resid: 0.9375\n",
      "20 ) resid: 0.029296875\n",
      "[ 60.          49.99267578  39.99267578  29.99267578  20.        ]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7c12f98>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFaZJREFUeJzt3X+Mm/ddwPH3J9za9UdyDT+aVCsUBqIrbKOtWACt7Jzk\nkiZDohNMAw3YRtISDQqjRGztBMoVTdM62nQrYoRyMBW0ibHBaNMfSS5NvGl/rBTaNKG0BSQKK1rC\nxmjSrghty4c/Hrt1Ul/O9tn32I/fL8maz/c89leW9827X3+fJDITSVI1LCt7AJKk/nFSl6QKcVKX\npApxUpekCnFSl6QKcVKXpArpaFKPiMmI+HREPBERj0fEj0XEyojYFxFPRcTeiJgc9GAlSWfWaal/\nFLg/My8DfgR4ErgR2J+ZlwIHgJsGM0RJUqdioYuPImIF8Ghmfv9pjz8JTGXmsYhYDdQz8zWDG6ok\naSGdlPr3AV+NiI9HxCMRcWdEnAusysxjAJl5FLhwkAOVJC2sk0l9ArgS+MPMvBL4OsXSy+mJ7983\nIEklm+jgmGeAL2Xm3zd+/muKSf1YRKxqWX75r3YnR4STvST1IDOj23MWLPXGEsuXIuIHGw+tBx4H\n7gHe1XjsncDd8z1HvZ68+tXJL/1S8rWvJZneer3t2LGj9DFU5eZ76fs5zLdedbr75TeAT0TEIYrd\nLx8EbgE2RMRTFBP9h+Y7eWoKDh+GCy6A174Wdu/uebySpDPoZPmFzHwMeEObX013+kLnnQd33AE/\n+7OwZQt8+tPw0Y/CypWdPoMkaSFLfkVps9onJ+F1r4N7713qEYy2Wq1W9hAqw/eyv3w/h8OC+9QX\n/QIROd9rfO5zRbVfdRV85CNWuyQ1RQQ5iC9KB6lZ7StWWO2S1A+llnqrer2o9p/8Satdkkay1FvV\naqdW+333lT0iSRo9Q1Pqrep12LrVtXZJ42vkS71VrQaPPQbLl1vtktSNoSz1VgcPFtXuWrukcVKp\nUm+1dm2x1m61S9LChr7UWzWr/U1vgttvt9olVVdlS71Vs9rPP99ql6R2RqrUWzWrfWqqqPYLLuj7\nS0hSacai1Fs1q/3cc4tqv//+skckSeUb2VJvZbVLqpqxK/VWVrskFSpR6q0OHCiqvVaz2iWNrrEu\n9Vbr1sGRI1a7pPFUuVJvZbVLGlWWehvr1hVr7eecU1T7Aw+UPSJJGqyOSj0ingaOAyeBb2TmmojY\nAVwH/FfjsPdn5p4255ZW6q2a1b52LezcabVLGm6DLvWTQC0zr8jMNS2P78zMKxu3l03ow6RZ7a98\npdUuqbo6ndRjnmO7/lOkTMuXw8c+BnfdBb/6q8W/tPTss2WPSpL6p9NJPYG5iHg4Iq5refz6iDgU\nEbMRMTmA8Q1Es9rPPttql1Qtna6pX5SZX46I7wLmgOuBp4CvZmZGxAeAizJza5tzh2JNfT4PPgjX\nXltM9Lfd5lq7pOHQ65r6RCcHZeaXG//7lYj4LLAmM7/QcsifALvnO39mZubF+7VajVqt1u04B2b9\n+qLa3/veotrvvBM2by57VJLGTb1ep16vL/p5Fiz1iDgXWJaZz0fEecA+4GbgcGYebRxzA/CGzHx7\nm/OHutRbPfhgsUNm/fpih8zkyCwoSaqaQe5+WQV8ISIeBb4I7M7MfcCHI+JwRBwCpoAbun3xYbN+\nfXE16llnFdW+Z6j380jSy1X6itLFsNollckrSvusWe2veIXVLml0WOod2L+/2CFjtUtaKpb6AE1P\nn1rte/eWPSJJas9S71Kz2qeni33tVrukQbDUl0iz2icmrHZJw8dSX4RmtW/YALfearVL6h9LvQTT\n08XVqMuWWe2ShoOl3idzc0W1b9xotUtaPEu9ZBs2FGvtzWrft6/sEUkaR5b6ALRW+223wYoVZY9I\n0qix1IeI1S6pLJb6gFntknphqQ+pZrVHWO2SBs9SX0L79sF118HVVxc7ZKx2SfOx1EfAxo1FtYPV\nLmkwLPWSWO2SzsRSHzEbNxZXo2YW1T43V/aIJFWBpT4EmtW+aRP8/u9b7ZIs9ZHWrPaTJ612SYvT\nUalHxNPAceAk8I3MXBMRK4FPAZcATwNvy8zjbc611Luwd29R7Zs3W+3SOBt0qZ8Eapl5RWauaTx2\nI7A/My8FDgA3dfviermrry52yJw8Ca9/ffHX+0pSpzot9X8DfjQz/7vlsSeBqcw8FhGrgXpmvqbN\nuZZ6j6x2aXwNutQTmIuIhyPi2sZjqzLzGEBmHgUu7PbFdWbNav/Wt6x2SZ2Z6PC4N2bmlyPiu4B9\nEfEUxUTfat4cn5mZefF+rVajVqt1OczxNTkJs7NFtW/ZYrVLVVWv16nX64t+nq63NEbEDuB54FqK\ndfbm8svBzLyszfEuv/TJ8eOwfXtR7LOzxb+8JKmaBrb8EhHnRsT5jfvnARuBI8A9wLsah70TuLvb\nF1d3mtW+a1dR7du2wYkTZY9K0jDpZE19FfCFiHgU+CKwOzP3AbcAGxpLMeuBDw1umGq1aZNr7ZLa\n84rSEbdnD/zKr7jWLlWNV5SOKatdUitLvUKsdqk6LHVZ7ZIs9aqy2qXRZqnrFFa7NJ4s9TFgtUuj\nx1LXvKx2aXxY6mPGapdGg6WujljtUrVZ6mPMapeGl6WurlntUvVY6gKsdmnYWOpaFKtdqgZLXS9j\ntUvls9TVN1a7NLosdZ2R1S6Vw1LXQFjt0mix1NUxq11aOpa6Bs5ql4Zfx5N6RCyLiEcj4p7Gzzsi\n4pmIeKRx2zS4YWpYTE7C7Czs2gVbtsC2bXDiRNmjktTUTam/B3j8tMd2ZuaVjduePo5LQ85ql4ZT\nR5N6RFwMvBmYPf1XfR+RRobVLg2fTkv9duC3gdO/8bw+Ig5FxGxETPZ3aBoVVrs0PCYWOiAifgo4\nlpmHIqLW8quPAb+XmRkRHwB2AlvbPcfMzMyL92u1GrVard1hGmHNat+zp6h2d8hI3anX69Tr9UU/\nz4JbGiPig8AvAt8EzgGWA3+Tme9oOeYSYHdmvr7N+W5pHDPHj8P27UWxz87C9HTZI5JGT69bGrva\npx4RU8D2zPzpiFidmUcbj98AvCEz397mHCf1MeW+dql3ZexT/3BEHI6IQ8AUcMMinksV5Fq7tPS8\nolRLwmqXuuMVpRpqVru0NCx1LTmrXVqYpa6RYbVLg2Opq1RWu9Sepa6RZLVL/WWpa2hY7dJLLHWN\nPKtdWjxLXUPJate4s9RVKVa71BtLXUPPatc4stRVWVa71DlLXSPFate4sNQ1Fqx26cwsdY0sq11V\nZqlr7Fjt0stZ6qoEq11VY6lrrFntUsFSV+VY7aoCS11qsNo1zjou9YhYBvw98Exm/nRErAQ+BVwC\nPA28LTOPtznPUldprHaNqqUo9fcA/9Ty843A/sy8FDgA3NTti0uDZrVr3HQ0qUfExcCbgdmWh68B\n7mrcvwt4S3+HJvXH5CTMzsKuXbBlC2zbBidOlD0qaTA6LfXbgd8GWtdRVmXmMYDMPApc2OexSX11\nerXPzZU9Iqn/JhY6ICJ+CjiWmYcionaGQ+ddOJ+ZmXnxfq1Wo1Y709NIg9Os9r17YetW19o1POr1\nOvV6fdHPs+AXpRHxQeAXgW8C5wDLgc8CPwrUMvNYRKwGDmbmZW3O94tSDaXjx2H79mKdfXYWpqfL\nHpH0kl6/KO1qn3pETAHbG7tfPgz8d2beEhHvA1Zm5o1tznFS11Bzh4yGURn71D8EbIiIp4D1jZ+l\nkeMOGVWJV5RKLax2DQuvKJX6wGrXqLPUpXlY7SqTpS71mdWuUWSpSx2w2rXULHVpgKx2jQpLXeqS\n1a6lYKlLS8Rq1zCz1KVFsNo1KJa6VAKrXcPGUpf6xGpXP1nqUsmsdg0DS10aAKtdi2WpS0PEaldZ\nLHVpwKx29cJSl4aU1a6lZKlLS8hqV6csdWkEWO0aNEtdKonVrjOx1KURY7VrEBac1CPi7Ih4KCIe\njYgjEbGj8fiOiHgmIh5p3DYNfrhStUxOwuws7NoFW7bAtm1w4kTZo9IoW3BSz8z/A9Zm5hXA5cDm\niFjT+PXOzLyycdszyIFKVWa1q186Wn7JzBcad88GJoDmInnX6z2S2rPa1Q8dTeoRsSwiHgWOAnOZ\n+XDjV9dHxKGImI2IyYGNUhojVrsWo6vdLxGxAvgs8OvAV4CvZmZGxAeAizJza5tzcseOHS/+XKvV\nqNVqix23NBbcITM+6vU69Xr9xZ9vvvnmnna/dL2lMSJ+F/h6Zu5seewSYHdmvr7N8W5plBbh+HHY\nvr0o9tlZmJ4ue0RaCgPb0hgR39lcWomIc4ANwJMRsbrlsJ8B/rHbF5e0MNfa1Y1O1tQvAg5GxCHg\nIWBvZt4PfDgiDjcenwJuGOA4pbHnWrs64RWl0ghyrb36vKJUGiNWu+ZjqUsjzmqvJktdGlNWu1pZ\n6lKFWO3VYalLstplqUtVZbWPNktd0ims9vFkqUtjwGofPZa6pHlZ7ePDUpfGjNU+Gix1SR2x2qvN\nUpfGmNU+vCx1SV2z2qvHUpcEWO3DxlKXtChWezVY6pJexmovn6UuqW+s9tFlqUs6I6u9HJa6pIGw\n2kfLgqUeEWcDnwfOAiaAz2TmzRGxEvgUcAnwNPC2zDze5nxLXaoIq33pDKzUM/P/gLWZeQVwObA5\nItYANwL7M/NS4ABwU7cvLmm0WO3Dr6Pll8x8oXH3bIpaT+Aa4K7G43cBb+n76CQNnclJmJ2FXbtg\nyxbYtg1OnCh7VGrqaFKPiGUR8ShwFJjLzIeBVZl5DCAzjwIXDm6YkoaN1T6cJjo5KDNPAldExArg\nsxHxwxS1fsph850/MzPz4v1arUatVut6oJKGT7Pa9+wpqt219t7V63Xq9fqin6frLY0R8bvAC8C1\nQC0zj0XEauBgZl7W5ni/KJXGwPHjsH17UeyzszA9XfaIRtvAviiNiO+MiMnG/XOADcATwD3AuxqH\nvRO4u9sXl1QdzWr/4z92rb1MnaypXwQcjIhDwEPA3sy8H7gF2BARTwHrgQ8NbpiSRsXVVxdr7SdP\nFmvtc3Nlj2i8eEWppIHZuxeuu8619l54RamkodNa7a97ndW+FCx1SUti376i2jdtsto7YalLGmob\nN8Lhw1b7oFnqkpZcs9qvvhpuvdVqb8dSlzQymtUORbXv21fueKrEUpdUKqu9PUtd0kjauLHYIQNW\nez9Y6pKGRrPaN26E224b72q31CWNvGa1L1tmtffKUpc0lObm4Npri4n+1luLv1tmnFjqkiplw4ZT\nq33v3rJHNBosdUlDb26uWGvfsGF8qt1Sl1RZGzYU+9qt9oVZ6pJGyv79xVr79HSxQ6aq1W6pSxoL\n09NFtX/bt1nt7VjqkkZWlavdUpc0dqanix0yExNFte/ZU/aIymepS6qEZrWvXw87d45+tQ/yH56+\nOCIORMTjEXEkIn698fiOiHgmIh5p3Db1MnBJ6odmtZ911nhX+4KlHhGrgdWZeSgizgf+AbgG+Dng\nuczcucD5lrqkJfXgg7B162hX+8BKPTOPZuahxv3ngSeAVzVft9sXlKRBW7/+1Gp/4IGyR7R0uvqi\nNCK+F7gceKjx0PURcSgiZiNiBP8slFRVy5fDH/0RfPzj8O53F+X+7LNlj2rwOp7UG0svnwHe0yj2\njwGvzszLgaPAGZdhJKkM41btE50cFBETFBP6X2Tm3QCZ+ZWWQ/4E2D3f+TMzMy/er9Vq1Gq1HoYq\nSb1pVvtb31rskFm3rtjXfsEFZY/sJfV6nXq9vujn6WhLY0T8OfDVzPytlsdWZ+bRxv0bgDdk5tvb\nnOsXpZKGxnPPwXvfC/feC3feCZs3lz2i9nr9orST3S9vBD4PHAGycXs/8HaK9fWTwNPAtsw81uZ8\nJ3VJQ+fBB4tqX7u22CEzTNUOA5zUF8tJXdKweu45eN/7YPfu4at2J3VJ6tGBA8XumGGqdv/uF0nq\n0bp1xd/8+MpXFjtk7r+/7BH1zlKXpBbNaq/V4Pbby6t2S12S+mDdumJf+7nnjma1W+qSNI8yq91S\nl6Q+G8Vqt9QlqQMHDxbVPjW1NNVuqUvSAK1dW+yQGfZqt9QlqUvNan/Tm4pqX7my/69hqUvSEmlW\n+/nnF9V+331lj+gllrokLcKgqt1Sl6QSDFu1W+qS1Cf9rHZLXZJK1qz25cvLq3ZLXZIGoF4vqv2q\nq+AjH+m+2i11SRoitRo89hisWFFU+733Ls3rWuqSNGC9VLulLklDqlYr1tqXototdUlaQp/7HGzZ\nsnC1D6zUI+LiiDgQEY9HxJGI+I3G4ysjYl9EPBUReyNistsXl6RxMzV1arXv3t3f5+9k+eWbwG9l\n5g8DPwH8WkS8BrgR2J+ZlwIHgJv6OzS1U6/Xyx5CZfhe9pfvZ+fOOw/+4A/gE5+A3/xNeMc74H/+\npz/PveCknplHM/NQ4/7zwBPAxcA1wF2Nw+4C3tKfIelM/D9O//he9pfvZ/ea1X7BBfDa1/an2rv6\nojQivhe4HPgisCozj0Ex8QMXLn44kjRezjsP7rgDPvnJ/lR7x5N6RJwPfAZ4T6PYT//2029DJalH\nzWpfubJYa+9VR7tfImICuBd4IDM/2njsCaCWmcciYjVwMDMva3Ouk70k9aCX3S8THR73Z8A/NSf0\nhnuAdwG3AO8E7u7XoCRJvVmw1CPijcDngSMUSywJvB/4O+CvgO8G/h14W2Y+O9DRSpLOaOAXH0mS\nlk7f/pqAiNgUEU9GxD9HxPvmOeaOiPiXiDgUEZf367WrZqH3MiKmIuLZiHikcfudMsY5CiLiTyPi\nWEQcPsMxfi47tND76WezO/Nd3NnmuM4/o5m56BvFHw7/ClwCvAI4BLzmtGM2A/c17v8Y8MV+vHbV\nbh2+l1PAPWWPdRRuwFUU23APz/N7P5f9fT/9bHb3fq4GLm/cPx94arFzZ79KfQ3wL5n575n5DeAv\nKS5OanUN8OcAmfkQMBkRq/r0+lXSyXsJ4BfQHcjMLwBn2vXr57ILHbyf4GezY9n+4s5XnXZYV5/R\nfk3qrwK+1PLzM20Gdvox/9nmGHX2XgL8ROM/xe6LiB9amqFVkp/L/vOz2YOWizsfOu1XXX1GO93S\nqOHyD8D3ZOYLEbEZ+FvgB0sekwR+NnvS5uLOnvWr1P8T+J6Wny9uPHb6Md+9wDHq4L3MzOcz84XG\n/QeAV0TEty/dECvFz2Uf+dnsXuPizs8Af5GZ7a736eoz2q9J/WHgByLikog4C/h5iouTWt0DvAMg\nIn4ceDYbf3eMTrHge9m6nhYRayi2pn5taYc5UoL513n9XHZv3vfTz2ZP2l3c2aqrz2hfll8y81sR\ncT2wj+IPij/NzCciYlvx67wzM++PiDdHxL8CXwd+uR+vXTWdvJfAWyPi3cA3gP8Ffq68EQ+3iPgk\nUAO+IyL+A9gBnIWfy54s9H7iZ7MrjYs7fwE4EhGP8tLFnZfQ42fUi48kqUL8N0olqUKc1CWpQpzU\nJalCnNQlqUKc1CWpQpzUJalCnNQlqUKc1CWpQv4flLXNPqM82ggAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x78dc898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# ******************************************\n",
    "# 6. Settings for iterations\n",
    "# ******************************************\n",
    "# since the array gets updated after iterating we should store the old values to allow \n",
    "# a calculation to see the changes after each iteration\n",
    "# length of material\n",
    "L = 2\n",
    "# number of nodes (0 to 4 inclusive)\n",
    "nx = 5\n",
    "# distance between nodes\n",
    "dx = L/(nx-1)\n",
    "\n",
    "T = np.ones(nx)*10   \n",
    "T[0] = 60\n",
    "T[nx-1] = 20\n",
    "T_old = np.array(T)\n",
    "\n",
    "# set a convergence criteria value for the residual. After doing some runs you can\n",
    "# return here and see what effect the criteria has on the solution\n",
    "criteria = 0.01     \n",
    "\n",
    "# in case we never reach convergence, set a maximum number of iterations to solve for\n",
    "maxIteration = 5000   \n",
    "\n",
    "# set initial iteration counter to 0, and an initial residual counter larger than the criteria\n",
    "# initial iteration counter\n",
    "iteration =  0 \n",
    "\n",
    "# set any number larger than criteria to the iteration\n",
    "residual = 1  \n",
    "\n",
    "\n",
    "# ******************************************\n",
    "# 7. Solver\n",
    "# ******************************************\n",
    "# Use a while loop to let the iterations run until the residual is smaller than the criteria set\n",
    "# place the discretized formula defined in Step 4, into the while loop\n",
    "# store the array into T_old\n",
    "# increment the iteratrion number\n",
    "# check if the iteration has reached its maximum iterations and break the iterator if needed\n",
    "# e.g. if iteration > maxIteration:\n",
    "#        break\n",
    "# The residual equation can be written as np.sqrt(sum((T - T_old)**2))\n",
    "\n",
    "while residual > criteria:         # continuously loop while residual > criteria\n",
    "    T[1:-1] = (T[2:] + T[0:-2])/2\n",
    "    \n",
    "    residual = np.sqrt(sum(np.square(T - T_old)))   # calculate the residual\n",
    "    T_old = np.array(T)             # store the Temperature values into Told\n",
    "    iteration +=1                   # increment the iteration counter\n",
    "    if np.mod(iteration, 10) == 0:  # print every 10th iteration\n",
    "        print (iteration, ') '  'resid:', residual)  # print the results\n",
    "    \n",
    "    if iteration > maxIteration:   # ensure the while loop doesn't go on forever\n",
    "        break        \n",
    "    \n",
    "    \n",
    "            \n",
    "# ******************************************\n",
    "# 8. Process the results\n",
    "# ******************************************\n",
    "# plot the results\n",
    "print(T)\n",
    "plt.plot(np.linspace(0,L,nx), T)\n",
    "\n",
    "\n",
    "\n",
    "# ******************************************\n",
    "# 9. Conclusion\n",
    "# ******************************************\n",
    "# i) how confident are you about the solution reached is accurate?\n",
    "# ii) what effect does changing the residual criteria have on the solution?\n",
    "# iii) what happens if you increase the number of nodes?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "<h4 style=\"color:#7171C6;\">Exercise 7.2 steady state - heat flux and convection boundary conditions</h4>\n",
    "\n",
    "<img src=\"img/conduction-numerical-steady-heatFlux-conv-BC.png\" width=\"300px\"/>\n",
    "The problem is slightly more complicated whenre the boundary conditions are now a  heat flux at the left boundary, a convection heat transfer at the right boundary. The equations setup only changes for these boundaries and are given as: \n",
    "\n",
    "Left boundary condition\n",
    "$$T_0 = \\frac{q_{\\textrm{left}}\\Delta x}{k} + T_{1}$$\n",
    "Right boundary condition\n",
    "$$T_4 = \\frac{(hx/k)T_{\\infty} + T_3}{1 + hx/k} $$\n",
    "\n",
    "Note that the discretization for the internal nodes remain the same:\n",
    "\n",
    "Internal nodes\n",
    "$$\\frac{T_{i-1} - 2T_i+T_{i+1}}{\\Delta x^2}=0 $$\n",
    "\n",
    "* <span style=\"color:#00aa00;\"> The free stream temperature (T_inf = 20 degrees)</span>\n",
    "* <span style=\"color:#00aa00;\">The heat transfer coefficient is 15 W/m2K</span>\n",
    "* <span style=\"color:#00aa00;\">The heat flux at the left boundary is 200 W/m2</span>\n",
    "* <span style=\"color:#00aa00;\">The thermal conductivity of the material is 28 W/mK </span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "277 9.74525530378e-05\n",
      "[ 47.62059398  44.049117    40.47759463  36.90598434  33.33432867]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x5efe208>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEACAYAAABMEua6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHjhJREFUeJzt3X10lOW57/HvFQLZwdb3QlUaC/uIEhEEEXkRmAq0gt3o\namEfl6IHXS3YYkFFepTWJu4u2VhQBPSo3YiravGotSpW0Io62qMFQUCQxJZuXVTtFqsF2QICIdf5\n4x5pTBMyk8zkmXnm91krq5PJMzPXmjW9vPnN/WLujoiIxENJ1AWIiEj2qKmLiMSImrqISIyoqYuI\nxIiauohIjKipi4jESNpN3cxKzGy9mS1L/X66mf0+dd+rZjYgd2WKiEg6MhmpTwc2N/j9ZqDK3fsB\nVcDcbBYmIiKZS6upm1k3YCywuMHd9cARqdtHAu9ltzQREclUaZrXzQdm8vcmDnA18IyZ3QIYMCTL\ntYmISIZaHKmb2XnANnffQGjen/keMN3dKwgNfkluShQRkXRZS3u/mNlsYCJQB5QDXwQeA77p7kc1\nuO5jdz+iicdrcxkRkVZwd2v5qs9rcaTu7rPcvcLdewAXAs+7+yXAX8xsBICZjQT+eIjn0E+Wfqqq\nqiKvIS4/ei/1fubzT2ulm6k3ZTKwwMw6AJ+mfhcRkQhl1NTd/UXgxdTtlwHNTRcRySPtsqK0Df+S\nkEYSiUTUJcSG3svs0vuZH1r8orTNL2DmQ4c6t98Op5+e05cSEYkNM8Nz8UVpNlx6KXzjG3DllbB9\ne3u8oohIcWqXpj55MtTUQH099OoFixeH2yIikl3tEr80fI1162DqVDhwAO64A848M6cvLyJSkPI6\nfmmof394+WX4/vdh3Lgwiv/ww/auQkQkniLZT72kBCZNgtpa6NwZKivhzjvD6F1ERFqv3eOXpmzc\nCD/4AezcGSKZIdoaTESKXGvjl7xo6hDmsj/4IMycCaNHw803Q9euOS1NRCRvFUym3hwzuOgiePNN\n6NIFeveGBQugri7qykRECkfejNQbq60Nkcy2bXD77TBiRA6KExHJUwUfvzTFHX79a7jmGhg6FObO\nhRNOyHKBIiJ5qODjl6aYwbe/HRYu9egBffuGxr5vX9SViYjkp7weqTf2pz/B9Onwn/8JixaFL1RF\nROIolvFLU9zhN78Jzb1/f7j1VqioyNrTi4jkhVjGL00xg3/5F9i8Gfr0CY39ppvg00+jrkxEJHoF\n19Q/U14OP/kJrF0bfnr3hqeeiroqEZFoFVz80pynn4Zp0+Dkk+G22+Cf/znnLykikjNFE78059xz\nYdOmMPXxrLPCKH737qirEhFpX2k3dTMrMbN1ZraswX0/MLNaM9tkZnNyU2L6ysrguutg/Xr4wx/C\nRmGPPabj9ESkeKQdv5jZ1cAZwOHuPs7MvgZcD4x19zozO9bd/2ET3faKX5ry/PNhVWq3brBwYYhm\nREQKQU7jFzPrBowFFje4+wpgjrvXATTV0KN2zjmwYUM4Sm/o0DCK/+STqKsSEcmddOOX+cBMoOGQ\nuycw3MxWmdkLZjYg69VlQceOYZuBTZvgvffCcXoPPaRIRkTiqcWmbmbnAdvcfQPQ8J8CpcBR7j4I\n+CHwcG5KzI7jjoP774elS2H2bBg5Msx1FxGJk9I0rhkKjDOzsUA58EUzuw94B/g1gLuvMbN6MzvG\n3T9q/ATV1dUHbycSCRKJRBZKb51hw+C11+CuuyCRgEsugepqOPzwyEoSESGZTJJMJtv8PBnNUzez\nEcCM1BelU4Dj3b3KzHoCz7r7iU08JrIvSlvywQcwaxYsXx4O5Zg4MaxYFRGJWhTz1JcAPcxsE7AU\nuLQNzxWJLl1g8eIw7XHBgjCK37Ah6qpERFovNitK2+rAAbjnHrjhBpgwAX76UzjqqKirEpFiVfQr\nStuqQweYPDns3V5fH2bJ3HNPuC0iUig0Um/GunUwdWpo6nfcAQPycsKmiMSVRupZ1r8/vPwyfO97\nYavfyZPhw7xbXiUi8nlq6odQUgKTJoVDsDt3DnvJ3HlnyN9FRPKR4pcMbNwY9pLZuTNEMkOGRF2R\niMRV0RxnFzV3ePBBmDkznJF6883QtWvUVYlI3ChTbydmcNFF8OabYZ57795hjntdXdSViYhopN5m\ntbUhktm2DW6/HUaMiLoiEYkDxS8RcodHHw27QZ59NsydCyecEHVVIlLIFL9EyAzGjw+j9h49oG/f\n0Nj37Yu6MhEpNhqp58CWLTB9Orz1FixaFL5QFRHJhOKXPOMOv/lNaO79+8Ott0JFRdRViUihUPyS\nZ8zCStTNm6FPn9DYb7oJPv006spEJM7U1HOsvBx+8hNYuzb89O4NTz0VdVUiEleKX9rZ00/DtGlw\nyilw223hi1URkcYUvxSIc88Nh2APGQIDB0JVFezeHXVVIhIXauoRKCuD666D9evDytTKynD6kv5B\nIyJtpfglDzz/fFiV2q0bLFwIJ58cdUUiEjXFLwXsnHPC2ajf+AYMHRpG8Z98EnVVIlKI0m7qZlZi\nZuvMbFmj+2eYWb2ZHZ398opHx45hm4FNm+C998Jxeg89pEhGRDKTyUh9OlDT8A4z6waMBrZms6hi\ndtxxcP/9sHQpzJ4NI0eGue4iIulIq6mnmvdYYHGjP80HZma7KIFhw+C11+Bb34JEIozid+6MuioR\nyXfpjtQ/a94HwwAzOx94x9035aIwgdJSuPLKMFL/+OMwt/3++xXJiEjzSlu6wMzOA7a5+wYzS6Tu\nKweuJ0QvBy9t7jmqq6sP3k4kEiQSidZVW6S6dIF77oHVq2HqVLj77rB3++mnR12ZiGRLMpkkmUy2\n+XlanNJoZrOBiUAdUA58EVgBDAN2E5p5N+A9YKC7f9Do8ZrSmEUHDoQGf8MNMGEC/PSncNRRUVcl\nItmWsymN7j7L3SvcvQdwIfC8u09w9y+7ew937w68C/Rr3NAl+zp0gMmToaYmNPhevUKTr6+PujIR\nyQfZmqfuHCJ+kew75hi4805YvhwWL4bBg8OGYSJS3LSiNAbq6+G+++D668N2v7Nnw7HHRl2ViLSF\nVpQWsZISmDQpHKfXuXPYS+bOO0M8IyLFRSP1GNq4Mewl89//DXfcEaIZESksOs5OPscdHnwQZs4M\nZ6TefDN07Rp1VSKSLsUv8jlmcNFFYWvfLl3CiUsLFkBdXdSViUguaaReJGprQySzbVtYuDRiRNQV\nicihKH6RFrnDo4+GfWTOPhvmzYPjj4+6KhFpiuIXaZEZjB8fRu09ekCfPjB3LuzbF3VlIpItGqkX\nsS1bYPp0eOstWLQofKEqIvlB8Yu0ijs8+SRcdRX07w+33goVFVFXJSKKX6RVzGDcuLC972mnQb9+\ncNNN8OmnUVcmIq2hpi4AlJdDVVXYP2bt2jAF8qmnoq5KRDKl+EWa9PTTMG1aOJjjttvCF6si0n4U\nv0hWnXtuOAR78GA488wwit+9O+qqRKQlaurSrLKysPPjhg1hZeqpp8Ljj+s4PZF8pvhF0vbcc2FV\nakUFLFwIPXtGXZFIfCl+kZwbORJefz3MZx8yJIziP/kk6qpEpCE1dclIx44wY0bI2999Nxyn9/DD\nimRE8oXiF2mT3/0OrrwyHK+3aFHI3UWk7RS/SCSGDYPXXoNvfQsSiTCK37kz6qpEilfaTd3MSsxs\nvZktS/3+MzOrNbMNZvaomR2euzIln5WWhtH65s2wY0eIZB54QJGMSBTSjl/M7GrgDOBwdx9nZqOA\n59293szmAO7u1zfxOMUvRWb1apg6NaxSvf126Ns36opECk9O4xcz6waMBRZ/dp+7r3T3+tSvq4Bu\nmb64xNNZZ4XGfskl8PWvh2mQ27dHXZVIcUg3fpkPzASaG3JfDqzISkUSCx06wOTJUFMTjtDr1QuW\nLIH6+pYfKyKtV9rSBWZ2HrDN3TeYWQKwRn//EbDf3Zc29xzV1dUHbycSCRKJRCvLlUJzzDFw553w\n3e+GSObuu+GOO2DAgKgrE8kvyWSSZDLZ5udpMVM3s9nARKAOKAe+CPza3S81s0nAd4Fz3H1vM49X\npi5AGKXfd19YtDRuXNji99hjo65KJD/lLFN391nuXuHuPYALCV+OXmpm5xIimXHNNXSRhkpKYNKk\ncJzeP/0TVFbCXXfBgQNRVyYSH22Zp74I+ALwrJmtM7P/k6WaJOaOPBIWLICVK2Hp0rAL5O9/H3VV\nIvGgFaUSKXd48EGYOTPsKXPzzdC1a9RViURPK0qlIJnBRReFSOZLXwonLi1YEGbMiEjmNFKXvFJb\nG+a1b9sWFi6NGBF1RSLRaO1IXU1d8o47PPooXHMNnH02zJsHxx8fdVUi7Uvxi8SGGYwfH0bt3btD\nnz4wdy7s2xd1ZSL5TyN1yXtbtsD06fD222F731Gjoq5IJPcUv0isucOTT8JVV0H//nDrreFYPZG4\nUvwisWYWVqFu3gynnQb9+oUVqXu17E3kc9TUpaCUl0NVFaxdC2vWhCmQy5dHXZVI/lD8IgVtxYqQ\nt59yCtx2G/ToEXVFItmh+EWK0pgx4RDswYNh4MAwit+zJ+qqRKKjpi4Fr6ws7Py4fj28+WbYKOzx\nx3WcnhQnxS8SO889F1alnnhi2HKgZ8+oKxLJnOIXkZSRI+H118N89iFDwih+166oqxJpH2rqEksd\nO8KMGSFvf/fdcJzeww8rkpH4U/wiReF3v4MrrwwnLS1aFHJ3kXym+EXkEIYNg9degwsuCDs/zpgB\nO3dGXZVI9qmpS9EoLQ1foG7eDDt2hEjmgQcUyUi8KH6RorVqFUydCp07h73b+/aNuiKRv1P8IpKh\nQYPg1Vdh4kT4+tfDKH779qirEmmbtJu6mZWkDphelvr9KDP7rZn9wcyeMbMjclemSG506ABTpkBN\nDezfHyKZJUugvj7qykRaJ5OR+nSgpsHv1wEr3f1k4Hng+mwWJtKejjkG7roLnnoKfv7zsO3A2rVR\nVyWSubSaupl1A8YCixvcfT7wi9TtXwAXZLc0kfZ3xhnwyitwxRXwzW+GUfyHH0ZdlUj60h2pzwdm\nAg2/8ezq7tsA3P19oEuWaxOJREkJXHZZ2EemrCzMab/rLjhwIOrKRFpW2tIFZnYesM3dN5hZ4hCX\nNjvFpbq6+uDtRCJBInGopxHJD0ceCQsXwne+ExYu/cd/hFkygwdHXZnEUTKZJJlMtvl5WpzSaGaz\ngYlAHVAOfBF4DBgAJNx9m5l9GXjB3Xs18XhNaZSC5w5Ll8IPfxhmysyZA127Rl2VxFnOpjS6+yx3\nr3D3HsCFwPPufgnwJDApddn/Ap7I9MVFCoUZXHwx1NaGL1V79w6j+Lq6qCsT+by2zFOfA4w2sz8A\nI1O/i8Ta4YfDvHnw4ovwxBPhEOyXXoq6KpG/04pSkVZyh0cfhWuuCXvLzJ0Lxx8fdVUSF1pRKtLO\nzGD8+BDJfPWr0KdPGMXv2xd1ZVLMNFIXyZItW8Ih2G+/Hbb3HTUq6oqkkLV2pK6mLpJF7vDkk3DV\nVWEh0y23QEVF1FVJIVL8IpIHzGDcuLC9b+/e4YvU2bNh796oK5NioaYukgPl5VBVBWvWhJ0ge/eG\nFSuirkqKgeIXkXawYkXI23v1gvnzoUePqCuSfKf4RSSPjRkTDsEeNAgGDoTqatizJ+qqJI7U1EXa\nSVkZXH89rFsX9m+vrITHH9dxepJdil9EIvLcc+G0pRNPhAULoGfPqCuSfKL4RaTAjBwJGzaE+exD\nhoRR/K5dUVclhU5NXSRCnTrBjBmwcSO88074IvXhhxXJSOspfhHJIy+9FPZu/9KXwqrUysqoK5Ko\nKH4RiYHhw8MXqRdcACNGhFH8zp1RVyWFRE1dJM+UloYvUN94A7ZvD5HMAw8okpH0KH4RyXOrVsHU\nqdC5czhOr2/fqCuS9qD4RSSmBg0KWw1MnAijR4dR/I4dUVcl+UpNXaQAdOgAU6aEvdv374dTToEl\nS6C+PurKJN8ofhEpQGvXhlkyECKZAQOirUeyT/GLSBEZMABeeSWM3r/5zfC/H30UdVWSD1ps6mZW\nZmarzWy9mW0ys6rU/X3N7Pep+181M40VRNpRSQlcdlmIZMrKwiyZu+6CAweirkyilFb8Ymad3X23\nmXUAXgamA/8G3OLuvzWzMcAP3f1rTTxW8YtIO9i4MUQyu3aFSGbw4KgrkrbIafzi7rtTN8uAUqA+\n9XNE6v4jgfcyfXERyZ4+feDFF+Gaa8KB2JdfDh98EHVV0t7SaupmVmJm64H3gWfdfQ1wNTDPzP4M\n/Ay4Pndlikg6zODii0Mkc/TRcOqpsHAh1NVFXZm0l4xmv5jZ4cBjwDRgMvCCuz9uZuOBKe4+uonH\neFVV1cHfE4kEiUSirXWLSBpqasK89r/+NUQyw4dHXZE0J5lMkkwmD/5+4403tip+yXhKo5ndAOwG\nfuzuRzW4/2N3P6KJ65Wpi0TIHX71q7CPzPDh8LOfwfHHR12VtCRnmbqZHWtmR6RulwOjgVrgL2Y2\nInX/SOCPmb64iOSeGUyYECKZioqQvc+bFxYxSfy0OFI3s9OAXxD+A1ACPOTuN5nZUGAB0AH4FPi+\nu69v4vEaqYvkkT/+MRyCvXVr2N535MioK5KmtHakrhWlIkXIHZYtg6uugjPPDCP3ioqoq5KGtKJU\nRNJmBuef//cDsPv1g9mzYe/eqCuTtlJTFyli5eVQXQ1r1sDq1dC7N6xYEXVV0haKX0TkoOXLQ95e\nWQnz50OPHlFXVLwUv4hIm40dG05cOuuskLVXV8OePVFXJZlQUxeRzykrg1mzYP162Lw5jNqfeELH\n6RUKxS8ickgrV4ZVqV/9athy4KSToq6oOCh+EZGcGDUKXn89zGcfPDiM4nftiroqaY6auoi0qFMn\nuPbasL3v1q1h7/ZHHlEkk48Uv4hIxl56Kezd3qVLiGQqK6OuKH4Uv4hIuxk+HNatg3HjYMSIMIrf\nuTPqqgTU1EWklUpLYdq0MAXyo49CJPPLXyqSiZriFxHJilWrYOpUOOywsHd7nz5RV1TYFL+ISKQG\nDYJXXw0nL40eHUbxO3ZEXVXxUVMXkazp0AGmTAkbhe3bFyKZe++F+vqoKyseil9EJGfWrg2zZMxC\nJHPGGVFXVDgUv4hI3hkwAF55BSZPhvPOgyuuCF+qSu6oqYtITpWUwGWXheP0OnUKc9rvvhsOHIi6\nsnhS/CIi7er110Mks3t3iGQGD466ovyk+EVECkLfvmFF6jXXwPjxcPnl8MEHUVcVHy02dTMrM7PV\nZrbezDaZWVWDv/3AzGpT98/JbakiEhdmYepjbS0cfTScemrYbqCuLurKCl9a8YuZdXb33WbWAXgZ\nmAZ0BmYBY929zsyOdfcPm3is4hcROaSamrC971//GiKZ4cOjrih6OY1f3H136mYZUAo48D1gjrvX\npa75h4YuIpKOysqwb/sNN8DEieHnL3+JuqrClFZTN7MSM1sPvA886+5rgJ7AcDNbZWYvmNmAXBYq\nIvFmBhMmhEimoiJsMzBvHuzfH3VlhaU0nYvcvR7oZ2aHA4+Z2ampxx7l7oPM7EzgYaDJY2qrq6sP\n3k4kEiQSiTaWLSJxddhhMHs2TJoUDsFesgQWLQqHdMRZMpkkmUy2+XkyntJoZjcAu4GRwM3u/mLq\n/j8BZ7n7R42uV6YuIq3iDsuWwVVXhYOwb7kFvvKVqKtqHznL1M3sWDM7InW7HBgN1AKPA+ek7u8J\ndGzc0EVE2sIMzj8/fJFaWQn9+oVR/N69UVeWv9LJ1I8DXjCzDcBq4Bl3Xw7cC/Qws03AUuDS3JUp\nIsWsvByqq8MukKtXQ+/esGJF1FXlJ60oFZGCs3x5yNtPPRXmz4fu3aOuKPu0olREisbYseHEpYED\nQ9ZeXQ179kRdVX5QUxeRglRWBrNmhbNSN28OmfsTT+g4PcUvIhILK1eGVandu8OCBXDSSVFX1DaK\nX0SkqI0aFXaAPOecsPPjrFmwa1fUVbU/NXURiY1OneDaa2HjRti6NRyn98gjxRXJKH4Rkdh66aWw\nd3uXLmEXyMrKqCtKn+IXEZFGhg8PX6Sefz6MGBFG8Tt3Rl1Vbqmpi0islZaGL1DfeAP+9rcQyfzy\nl/GNZBS/iEhRWbUKpk4NG4fdfnvYDTIfKX4REUnDoEFhu4GLL4bRo2HaNNixI+qqskdNXUSKTocO\nMGVK2Chs374Qydx7L9TXR11Z2yl+EZGit3ZtmCVjFiKZM86IuiLFLyIirTZgALzyCkyeDOedB1dc\nAR8V6EbiauoiIkBJCVx2WThOr1OnMKf97rvhwIGoK8uM4hcRkSa8/nqIZPbsCZHMoEHt+/qKX0RE\nsqhv37Ai9eqr4dvfhssvhw8+iLqqlqmpi4g0wyxMfaythaOPDodyLFoEdXVRV9Y8xS8iImmqqQmr\nUz/8MEQyw4bl7rVyefB0mZmtNrP1ZrbJzKoa/X2GmdWb2dGZvriISCGprAz7tv/4x2EEP3Ei/Nd/\nRV3V57XY1N19L/A1d+8HnA6MMbOBAGbWDRgNbM1plSIiecIMJkwIkUxFBZx2GtxyC+zfH3VlQVqZ\nurvvTt0sA0qBz/KU+cDMHNQlIpLXDjsMZs8O89tXrgxfrD73XNRVpdnUzazEzNYD7wPPuvsaMxsH\nvOPum3JaoYhIHuvZE5Yvh3//d/jOd+Bf/xXeeSe6etIdqden4pduwEAzOw2YBTTM1zMO9EVE4sAs\n7NleUxNy9379wih+794Iasl0ZoqZ3UCIX64EdhOaeTfgPWCgu3/Q6Hqvqvp7708kEiQSibZVLSKS\nx956K8xvr6kJJy6NGdPyY5LJJMlk8uDvN954Y6tmv7TY1M3sWGC/u39sZuXAM8Acd1/e4Jq3gf7u\nvr2Jx2tKo4gUpeXLYfr0ML99/nzo3j39x+ZyRelxwAtmtgFYDTzTsKGnOIpfREQ+Z+zYcOLSwIFw\n5plQXR22HcglLT4SEWkHf/4zzJgRtvm97TYYNy5k8c1p7UhdTV1EpB2tXBlWpXbvDgsWwEknNX2d\nNvQSESkAo0aFHSBHjoTBg+FHP4Jdu7L3/GrqIiLtrFOnEMVs3Ahbt4bj9B55BLIRaih+ERGJ2Esv\nhb3bu3QJu0D26qX4RUSkYA0fDuvWhQVMw4fDtde2/rnU1EVE8kBpafgC9Y03YPs/rPhJn+IXEZE8\npPhFRETU1EVE4kRNXUQkRtTURURiRE1dRCRG1NRFRGJETV1EJEbU1EVEYkRNXUQkRtTURURiRE1d\nRCRG1NRFRGKkxaZuZmVmttrM1pvZJjOrSt3/MzOrNbMNZvaomR2e+3JFRORQWmzq7r4X+Jq79wNO\nB8aY2UDgt8Cp7n46sAW4PqeVCgDJZDLqEmJD72V26f3MD2nFL+6+O3WzDCgNd/lKd69P3b8K6JaD\n+qQR/R8ne/ReZpfez/yQVlM3sxIzWw+8Dzzr7msaXXI5sCLbxYmISGbSHanXp+KXbsBZZlb52d/M\n7EfAfndfmqMaRUQkTRmffGRmNwC73P1WM5sEfBc4J5W9N3W9jj0SEWmF1px81GJTN7NjCSPxj82s\nHHgGmAPUA7cAw939o1bUKyIiWVaaxjXHAb8wsxJCXPOQuy83sy1AJ+BZMwNY5e7fz12pIiLSkpwf\nPC0iIu0naytKzexcM3vTzP5oZv+7mWsWmtmW1IKl07P12nHT0ntpZiPMbIeZrUv9/DiKOguBmd1j\nZtvMbOMhrtHnMk0tvZ/6bGbGzLqZ2fNmtjm1uHNaM9el/xl19zb/EP7j8CfgRKAjsAE4pdE1Y4Cn\nUrfPIsQ1WXn9OP2k+V6OAJZFXWsh/ABnExbNbWzm7/pcZvf91Gczs/fzy8DpqdtfAP7Q1t6ZrZH6\nQGCLu2919/3A/wXOb3TN+cB9AO6+GjjCzLpm6fXjJJ33EiDjb8WLkbv/P2D7IS7R5zIDabyfoM9m\n2tz9fXffkLr9CVALnNDosow+o9lq6icA7zT4/d0mCmt8zXtNXCPpvZcAg1P/FHuq4boByZg+l9mn\nz2YrmNlXCf8KWt3oTxl9RtOZ/SL55zWgwt13m9kY4HGgZ8Q1iYA+m61iZl8AfgVMT43YWy1bI/X3\ngIoGv3dL3df4mq+0cI2k8V66+yee2o/H3VcAHc3s6PYrMVb0ucwifTYzZ2alhIZ+v7s/0cQlGX1G\ns9XU1wD/w8xONLNOwIXAskbXLAMuBTCzQcAOd9+WpdePkxbfy4Z5WmrHTHP3v7VvmQXFaD7n1ecy\nc82+n/pstsoSoMbdFzTz94w+o1mJX9z9gJldSdiOtwS4x91rzWxK+LP/3MOCpbFm9idgF3BZNl47\nbtJ5L4HxZvY9YD+wB/if0VWc38xsKZAAjjGzPwNVhEVz+ly2QkvvJ/psZsTMhgIXA5tSmyY6MIsw\n+61Vn1EtPhIRiREdZyciEiNq6iIiMaKmLiISI2rqIiIxoqYuIhIjauoiIjGipi4iEiNq6iIiMfL/\nAQauY5kIGlPpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x5eba278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# ******************************************\n",
    "# 1. PRELOAD MODULES (if you haven't done so from exercise 7.1)\n",
    "# ******************************************\n",
    "import numpy as np                \n",
    "from matplotlib import pyplot as plt   \n",
    "%matplotlib inline\n",
    "\n",
    "# ******************************************\n",
    "# 2. Set up the geometry (same geometry as exercise 7.1)\n",
    "# ******************************************\n",
    "# length of material\n",
    "L = 2\n",
    "# number of nodes (0 to 4 inclusive)\n",
    "nx = 5\n",
    "# distance between nodes\n",
    "dx = L/(nx-1)\n",
    "\n",
    "# ******************************************\n",
    "# 3. Define the variables\n",
    "# ******************************************\n",
    "# define the freestream temperature, heat flux, material conductivity, heat transfer coefficient\n",
    "# free stream temperature for convection [K]\n",
    "T_inf = 20\n",
    "\n",
    "# heat flux at left boundary [W/m2]\n",
    "qLeft = 200\n",
    "\n",
    "# material conductivity [W/mK]\n",
    "kCond = 28\n",
    "\n",
    "# heat transfer coefficient [W/m2K]\n",
    "hConv = 15 \n",
    "\n",
    "\n",
    "# ******************************************\n",
    "# 4. Settings for iterations\n",
    "# ******************************************\n",
    "# Set initial node values with a temperatur of 50-deg, later test the effect of different numbers\n",
    "T = np.ones(nx)*50\n",
    "\n",
    "# Store node values in a temporary array, called T_old\n",
    "T_old = np.array(T)\n",
    "\n",
    "# set a convergence criteria\n",
    "# return here and see what effect the criteria has on the solution\n",
    "criteria = 0.0001     \n",
    "\n",
    "# set maximum number of iterations\n",
    "maxIteration = 1000   \n",
    "\n",
    "# set initial iteration counter to 0\n",
    "iteration =  0    \n",
    "\n",
    "# set residual to 1\n",
    "residual = 1      \n",
    "\n",
    "\n",
    "# ******************************************\n",
    "# 5. Solve the problem (similar to exercise 7.1)\n",
    "# ******************************************\n",
    "# use a while loop, write the boundary condition equations, then the internal node equation\n",
    "# check residual, and convergence criteria, \n",
    "# store the array into T_old\n",
    "# increment the iteratrion number\n",
    "# check if the iteration has reached its maximum iterations and break the iterator if needed\n",
    "# e.g. if iteration > maxIteration:\n",
    "#        break\n",
    "# The residual equation can be written as np.sqrt(sum((T - T_old)**2))\n",
    "\n",
    "while residual > criteria:         # continuously loop while residual > criteria\n",
    "    T[0] = qLeft*dx/kCond + T[1]   # left heat flux boundary\n",
    "    T[nx-1] = (hConv*dx/kCond * T_inf + T[nx-2])/(1+hConv*dx/kCond) # right convection boundary\n",
    "    T[1:-1] = (T[2:] + T[0:-2])/2  # internal nodes\n",
    "    \n",
    "    residual = np.sqrt(sum(np.square(T - T_old)))   # calculate the residual\n",
    "    T_old = np.array(T)             # store the Temperature values into Told\n",
    "    iteration +=1\n",
    "    #print (iteration, residual)  # might be too many to print\n",
    "    if iteration > maxIteration:\n",
    "        break                      # stop the iteration if the maximum has been reached\n",
    "\n",
    "print(iteration, residual)\n",
    "print(T)\n",
    "plt.plot(np.linspace(0,L,nx), T)\n",
    "\n",
    "\n",
    "# ******************************************\n",
    "# 6. Conclusion\n",
    "# ******************************************\n",
    "# i) what are the left boundary and right boundary temperatures?\n",
    "# ii) what is the effect of a higher heat flux?\n",
    "# iii) what is the effect of a higher heat transfer coefficient?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<h4 style=\"color:#7171C6;\">ADDITIONAL STUDY EXERCISE: 2D and 3D models</h4>\n",
    "Solve the steady state solution for a 2D rectangular block with the following temperature boundary conditions:\n",
    "\n",
    "<img src=\"img/2D-heat-conduction.png\" width=\"200px\"/>\n",
    "\n",
    "The 2D-heat equation is:\n",
    "\n",
    "$$ \\frac{\\partial ^2 T}{\\partial x^2} + \\frac{\\partial ^2 T}{\\partial y^2} = 0$$\n",
    "\n",
    "Use the same discretization in the 1-D case, but now repeat in the other axis to complete the 2D form.\n",
    "\n",
    "$$ \\frac{T_{i+1,j} - 2 T_{i,j} + T_{i-1,j}}{\\Delta x^2} + \\frac{T_{i,j+1}-2 T_{i,j} + T_{i,j-1}}{\\Delta y^2} = 0$$\n",
    "\n",
    "We set the cell spacing $\\Delta x = \\Delta y$ which simplifies the equation\n",
    "\n",
    "$$ T_{i+1,j} - 2 T_{i,j} + T_{i-1,j} + T_{i,j+1}-2 T_{i,j} + T_{i,j-1} = 0$$\n",
    "\n",
    "\n",
    "Recall this represents an interior node and is given by $T_{i,j}$\n",
    "$$T_{i,j} =\\frac{ T_{i+1,j} + T_{i-1,j} + T_{i,j+1} + T_{i,j-1}}{4}$$\n",
    "\n",
    "The interpretation here is that the temperature of each interior node is the average of the temperatures of the four neighboring nodes. \n",
    "\n",
    "(What do you think it would be for 3D problem? Can you write out the discretized 3D equation??)\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAM0AAAEZCAYAAAA0ZKHkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGHJJREFUeJzt3X34VGWdx/H3ByWVUBIfqESwSFNzvcyKKFPwoVJsZeuy\nzCzTdjfXfGBrN2yxLiGzy9ptS5LWHswNS9MoH0qtNEVFi0wlSbGkTISUQkVCXEP57h/n/ukwzMxv\nzsycefjN53VdczEz55773OfHfOfc5z7nfG9FBGZWv2GdboBZr3HQmOXkoDHLyUFjlpODxiwnB41Z\nTg6aGiQ9KOmQTrej1SSdJeni9HxXSWslqUb5v0rarV3t63ZdFzSS/ihpffqPfEzSDyXt0ul2tYOk\n0yUtkbRO0nJJl0l6TUGrC4CIeDgitot0wk7STZI+tEnBiG0j4o8FtaPndF3QkP1nHhkR2wEvA/4M\nfLmzTcpH0hYNfGYOcBpwKrA9sAdwJXBka1tnTYuIrnoADwKHlLw+Ari/5PV2wDyyYHoQOLNk2VnA\nxSWvxwMbgWHp9U3Ap4GFwFrgx8DokvIfAP4I/AWYWdoW4A3A7cATwEqyQN6y5LMbgY8AvwN+D5wP\n/FfZtl0FTK+wza8CngVeV+PvUmu7PwjcCvwn8Hha/+Ely3cDFgBPAj9JbZ9X/jcCPpPasT79feaU\nbNsrW9COE9J7a9O/x3b6+9bQd7TTDagVNMAI4H+Bi0qWzwOuSMvGA78FTiwJmnllQfNcWdA8AEwA\ntkqvP5uW7Q38FTgAGA58AfhbSVv2ByYCAsYB9wKnlwXNT4BRqe43ACtKlu8ArAN2rLDNJwEPDvJ3\nqbXdHwSeAT6U2vcvwMqSz96evsjDgQPTl3Zejb/Rh8rW/VxJ0DTUjlT+SeBV6fUYYK9Of9+GUtCs\nTb9UfwNWAK9Jy4al/5RXl5T/MHBjjqCZWbL8ZODa9PxTwCUly0akdR1SpZ3Tge+XBc3ksjL3Aoem\n56cAP6pS10zg9hp/k8G2+4PA70qWbZPaszOwa/o7blOy/Ds5g2Yj8Mom2zEi/Z++E9i609+zZh7d\neEwDMC0iRpP9Yp8G3CJpZ2BHYEtgeUnZh4A8AwWPljxfD4xMz18OPDywICLWA48NvJa0exqUeETS\nGuCc1J5SK8pezwPen56/H7i4SpseIzt+q6ae7X5+uyLi6fR0JNl2PVHy3sBnG9FwO9Lf8xiyH6pH\n0t/y1Q22o6O6NWgEEJkryH4J3wKsJutzjy8pO57sGAPgKbJftAG1vojlHiH7Vc4aII0g61IN+B9g\nKTAhIl4CnDnQzhLll4x/G5gmaV9gT7ID+0p+BoyVtH+V5auBDVTf7loeAbaXtE3Je+NqlK912Xsz\n7SAiro+ItwEvJevWfb2ez3Wbbg2a50maBrwEuC8iNgKXAedIGilpPPBRXvgFXwwclM49jAI+kWNV\n84F3SHqzpOFkAwalQbEtsDYi1kvak+wXs6aIWAn8KrXv+xHxTJVyy4CvAJdKmixpuKStJB0jaUba\n7u/V2O5abVie2jA71fsW4O/LipVu5yqyrlilujYClzfSDkk7Szoq/RhtIDu+e26wz3Wjbg2aH6bz\nNE8CZwPHR8T9adnpZN2qPwC3AN+OiIsAIuIGsqC6B7gD+GFZvVV/RSPiPrLjjkuBP5F1mUq7W/8O\nHCdpLfBV4Lt11v0tYB+yrlpVETGdbMRtLtkI3TLgH0q24TSqbHe1KkueHwdMStv0qdSmamXPA96d\nzpF9qcLyqn//QdoxDPgY2V5pNXAQdfzwdCOlA7ZiVyINI/u1WxERR1VYPodsaPkp4ISIWFx4o9pE\n0oFkw+C7dbot/U7ShcA7gFURsW96b3uyH9rxZKcb3hMRT9aqp117munAfZUWSDqC7Dhhd7Kh1wva\n1KbCpW7edHq07z4EXQS8vey9TwA3RMSrgRuB/xisksKDRtJYYCrwjSpFppG6LhGxCBglaUzR7Spa\nOu55gux8xHkdbo4BEbGQ7P+k1DRe6K5+i6xLXNOWLW5XJV8EPk520q+SXSgZ6iXr8+5CdkDas9Ix\n2MhBC1qn7RwRqwAi4tF0aqOmQvc0ko4k6z8uJhuhqXolrVmXGPQgv+g9zQHAUZKmkp0d3lbSvIg4\nvqTMSkrOjwBjqTDuL8lpc3pcRDT1ozlaivK+VQ2rIuKl9ZSTNCYiVkl6Kdk1dTUVGjQRMZPsEhEk\nTQb+rSxgAK4mG+q9TNIkYM3A7rLcwngdb7n6V/nOvjTqL7Ngp1ltWNEQXt/74SMz/5uJLOIEXd50\ndU8An6uz7BnZsWQl5T2eq8kuJP0c2WVAVw1WdzuOaTYj6SSyE/5fi4hrJU2VtIxsyPnEap/7CW/P\nAmbp59vQyttgdTvWM4TXd+YM5v3r8UwcsajY9dRJ0iXAFGAHScvJrlU8F/heuofoIeA9g9XTtqCJ\niJuBm9Pzr5YtO7Vd7bD+FRHvq7LosDz1dOsVAV2g4pUkXp91pnvWGyZ4fV2oG+57957GLCcHjVlO\nDhqznBw0Zjl5IMB6igcCzHqQg8YsJweNWU4OGrOcPBBgPaUbbun1nsYsJweNWU4OGrOcHDRmOXkg\nwHrKLlvXWfD/Nn9L0nTgn9LLr0fEnEba4D2N9YU0DeM/Aq8H9iPL293QnXgOGusXewGLIuKZiHiO\nLA/1uxqpqOi8Z1tJWiTp7jQB61kVykyWtEbSXenxySLbZH3rN8CBkrZPMxdMZdPUYXUrOoXTM5IO\nTtNTbAHcJum6iPhlWdFbKiVGN6vXrc/Bwo3Vl0fE/ZI+B1xPNs3H3TQ41UfhAwFpBizIZjXbksoZ\nDJ150+qyXZWksUey6TTY5y7fvEyaEuQiAEnnsGk65Lq1IwH6MEl3k00rd31E3FGh2JskLZZ0jaS9\ni26T9SdJO6V/x5HN/XlJI/W0Y0+zEXitpO2AKyXtnSZQGnAnMC514Y4gm2Jvj6LbZX3p+5JGk83E\n9pGIWNtIJe1MFrhW0k3A4ZTMVRMR60qeXyfpK5JGR8Tj5XXcPOuWLJ0qt5Hl7erNNET94ffALJ45\nZz1XDP9NpxsDQEQc1Ip6ih492zHNfUmaKPWtwP1lZcaUPJ9INjvbZgEDMHnWQSn/8FtxwHS7CcAs\ntjpzBu+ctU+nG9NSRe9pXgZ8K00fOAy4LOVufj6XM3C0pJPJdplPk02bbVbZ6DrLVRgIaJWih5yX\nAJtN812ayzki5pJNzmrWE3xFgFlODhqznBw0Zjn51gDrLV2QJMB7GrOcHDRmOTlozHJy0Jjl5IEA\n6y1Vbg1oJ+9pzHJy0Jjl5KAxy8lBY5aTBwKstzQxECDpo2S5zzYCS4ATI+Jveevxnsb6gqSXA6cB\n+0fEvmQ7jPc2Upf3NNZPtgBeLGkjMAL4UyOVeE9jfSEi/gR8geyezpXAmoi4oZG6vKexIWHBw9mj\nGkkvAaYB44EngfmS3hcRudM4OWist+xQ+e0pO8CU/V54PfvnmxU5DPjDQNIWST8A3kwDuc86nss5\nlZsj6YGUMHC/SmXMmrQcmCRpa0kCDgWWNlJRoUETEc8AB0fEa8mmNzgipWl6XkoQOCEidgdOAi4o\nsk3Wn1L+8PlkOZx/TZYK+WuN1NUNuZynAfNS2UWSRkkaExGrim6b9ZeImA3MbraebsjlvAubJqJe\nmd4z60rdkMu5bk5L20sKSkvbT7cGpGTTA7mcS61k08l1xqb3NuO0tL1k6Kal7XguZ+Bq4PhUZhLZ\nSScfz1jX6ngu5/R6qqRlwFPAiQW3yawpHc/lnF6fWmQ7zFrJVwRYb+mngQCzocJBY5aTg8YsJweN\nWU4eCLDeUu/0gQXynsYsJweNWU4OGrOcHDRmOXkgwHpLE9MHStoDuIzsRkiR3V/yqYiYk6ceB431\njYj4HfBayG6OBFYAV+Stx90z61eHAb+PiBqJnypz0Fi/Oga4tJEPuntmPeXZKlc533wr3Lywvjok\nDQeOAj7RSBscNDYkTD4weww4+9yaxY8A7oyIvzSyLnfPrB8dS4NdM3DQWJ+RNIJsEOAHjdbh7pn1\nlZS8cqdm6ig0aCSNJcueOYZs9qmvl59IkjQZuAr4Q3rrBxHxmSLbZb1rzXYj6yy5rrA2FL2neRb4\nWEQsljQSuFPSTyOiPI3TLRFxVMFtMWuJohOgPxoRi9PzdWRZ2iulnFWR7TBrpbYNBEjajWzmgEUV\nFr8pTbNxjaS929Ums0a0ZSAgdc3mA9PTHqfUncC4iFifpt24EtijUj3O5dxLCsrl3AUUUT7zRYtX\nIG0J/Ai4LiLOq6P8g8DrBmasKnk/PhUzOXvvc2Dp5wtqrbXWDEY+tZrzR5zCCbqciGiqGy4plsb4\nusrupYeaXl817eiefRO4r1rASBpT8nwiWSA/XqmsWTcoesj5AOA4YEmaoyaAmWSThUZEfA04WtLJ\nwAbgabIL6cy6VtG5nG8jm7u9Vpm5wNwi22HWSr6MxiwnX0ZjPeUxdqyz5EOFtcF7GrOcHDRmOTlo\nzHJy0Jjl5IEA6ylreEnDn02TJn8D2IfsVpUPRUSlayFrctBYPzkPuDYi3p0u7xrRSCUOGusLkrYD\nDoyIEwAi4llgbSN1+ZjG+sUrgNWSLpJ0l6SvSdqmkYocNNYvtgT2B+ZGxP7Aepz3zPrBanao+P7S\nBX/m/gV/rvXRFcDDEfGr9Ho+cEYjbXDQ2JCw15Sd2WvKC+k3r5p93ybLI2KVpIcl7ZESoR8KbFqo\nTg4a6yenA99JaWn/AJzYSCUOGusbEfFr4A3N1uOBALOcvKexnlL/rQHF8Z7GLCcHjVlOhQaNpLGS\nbpR0r6Qlkk6vUm6OpAdSwsD9imyTWbM6nss5JQicEBG7S3ojcAEwqeB2mTWs6Gw0jwKPpufrJA3k\nci5NgD6NbGYBImKRpFGSxkTEqiLbZr2pmVsDWmXQ7pmk0yRt3+yKauRy3gUonWF3JZWTpJt1hXr2\nNGOAOyTdRZYt8yeRM5ftILmc6+Zczr1k6OZyHnRPExGfBHYHLgROAB6Q9FlJdX1j080+84GLI+Kq\nCkVWAruWvB6b3tvM5FkHwU6zgLfigOl2E4BZbHXmDN45a59ON6al6ho9S3uWgeOTZ4HtgfmS6slE\nXjOXM3A1cDyApEnAGh/PWDcbtHsmaTrZl3o12f3VH4+IDZKGAQ8AM2p8dtBczhFxraSpkpYBT9Hg\nRXTWHx6rcmtAO9VzTDMaeFdEbJKyMCI2SnpHrQ/Wk8s5lTu1jnaYdYVBgyYizqqxbGlrm2PW/XwZ\njVlODhqznHxrgPWUajkC6iXpj8CTZMkCN0TExLx1OGis32wEpkTEE41W4O6Z9RvR5PfeQWP9JoDr\nJd0h6Z8bqcDdM+s3B0TEI5J2IguepRGxME8FDhrrKdVyBDyx4B6eWLBk0M9HxCPp379IugKYCDho\nrP9sP2Vftp+y7/OvH5x96WZlJI0AhqV7u14MvA2YnXddDhrrJ2OAKyQF2Xf/OxHx07yVOGisb0TE\ng2Q3QjbFo2dmOXlPYz2lJ3IEmNmmHDRmOTlozHJy0Jjl5IEA6ynN3hrQCkXncr5Q0ipJ91RZPlnS\nmjTb7l2SPllke8xaoeg9zUXAl0lpZ6u4JSKOKrgdZi1T6J4mXT062M0+KrINZq3WDQMBb0pTbFwj\nae9ON8ZsMJ0eCLgTGBcR69OUG1cCe1Qr7FzOvaSYXM6Pre/z6QMjYl1ErE/PrwOGSxpdrbxzOfeS\nPs/l3CRR5bhF0piS5xMBRcTjbWiTWcMK7Z5JugSYAuwgaTlwFvAiUh5n4GhJJwMbgKeBY4psj1kr\nFD0T2vsGWT4XmFtkG8xardMDAWa5rFvd3K0BabaLXwErGj0/2A1DzmbtNB24r5kKHDTWNySNBaaS\nzbPUMAeN9ZMvAh8nSxjYMB/T2NDw85vhFzdXXSzpSGBVRCyWNIUmLt9y0FhveaTKV3a3Q7PHgC+d\nXV7iAOAoSVOBbYBtJc2LiOPzNsHdM+sLETEzIsZFxCuB9wI3NhIw4KAxy83dM+s7EXEzUP0AaBDe\n05jl5D2N9ZZVnW6A9zRmuTlozHJy0Jjl5KAxy8kDAdZbVne6Ad7TmOXmoDHLyUFjllNHczmnMnMk\nPZASBjY9H6JZ0TqayzklCJwQEbtLeiNwATCp4DZZL3u00w3ofC7naaSAiohFwKjSXGhm3ajTxzS7\nAA+XvF6Z3jPrWj11nsa5nHtJMbmcu0Gn9zQrgV1LXo9N71XkXM69ZOjmcm7HnqZqLmfgauAU4DJJ\nk4A1EdEFF39b12pwIEDSVsAtZGmRtwTmR8TsRurqaC7niLhW0lRJy4CngBOLbI/1r4h4RtLBaVqX\nLYDbJF0XEb/MW1dHczmnMqcW2QazAQPTugBbkX33G8p/1uljGrO2kTRM0t1knbzrI+KORurpqdEz\ns6pWLICVC2oWiYiNwGslbQdcKWnviMid19lBY72l2q0BW0+BCVNeeP3L6sf4EbFW0k3A4TSQDN3d\nM+sLknaUNCo934bsvMX9jdTlPY31i5cB30rz0wwDLouIaxupyEFjfSEilgD7t6Iud8/McvKexnrL\nUL81wGwoctCY5eSgMcvJQWOWkwcCrLd4IMCs9zhozHJy0Jjl5KAxy8kDAdZb+mEgQNLhku6X9DtJ\nZ1RYPlnSGkl3pccni26TWTOKTqwxDDgfOBT4E3CHpKsiovw+hlsi4qgi22LWKkXvaSYCD0TEQxGx\nAfguWSractVSPJl1naKDpjzt7Aoqp519U5o14BpJexfcJrOmdMNAwJ3AuJSP6gjgSmCPDrfJulWt\ndPo1SBpLlmx/DLAR+HpEzGmkrqKDZiUwruT1ZmlnI2JdyfPrJH1F0uiIeLy8Mudy7iVdl8v5WeBj\nEbFY0kjgTkk/rXB8Paiiu2d3AK+SNF7Si4D3kqWifV7p1BqSJgKqFDDgXM69pbtyOUfEoxGxOD1f\nByylwRkqis6w+ZykU4GfkgXohRGxVNJJpNS0wNGSTgY2AE8DxxTZJjNJuwH7AYsa+XzhxzQR8WPg\n1WXvfbXk+VxgbtHtsKFibZX3bwUWDvrp1DWbD0wvPTTIoxsGAsxa4MD0GHDuZiUkbUkWMBdHxFWN\nrsnXnlk/+SZwX0Sc10wlDhrrC5IOAI4DDpF0d7pk6/BG6nL3zPpCRNwGbNGKuhw01mOqzi7ZNu6e\nmeXkoDHLyUFjlpODxiwnDwRYj1nV6QZ4T2OWl4PGLCcHjVlODhqznDwQYD3GVwSY9RwHjVlODhqz\nnBw0Zjl1PJdzKjNH0gMpYeB+RbfJetnKOh+bk3ShpFWS7mmmBYUGTUku57cDrwGOlbRnWZkjgAkR\nsTtwEnBBkW2q3++9vqHnIrLvYlOKHnJ+PpczgKSBXM6lCdqmkWU+JCIWSRolaUxEbHaR0aeXfZaz\nLzoHzptRcLOBJbPg79qwnqG8vsPg/BGnFLuOHCJioaTxzdZTdNBUyuU8cZAyK9N7m1+Z9z0IxOpL\nRra4mZv7/KxnmDHrC4WvZ6iv7xqmFr6Oduupk5sLZsIvAWY2lK4ql7uBb87eUPh6hvr6Zoy7HH4G\nJxS+pvbpeC7n9HrXQcoAcHBLmza4G7y+pp2xHNi9ZdU9BGfU270q7B6CooPm+VzOwCNkuZyPLStz\nNXAKcJmkScCaSsczEeE5bPpcROzWgmpEk/MhFTp6FhHPAQO5nO8FvjuQy1nSh1OZa4EHJS0Dvgp8\npMg2Wf+SdAlwO7CHpOWSTmyonohobcvMhriuuyKg3SdD2zmRbj0n11q8bTXX1+pJgiWNlXSjpHsl\nLZF0epVyvX0yOyK65kEWxMuA8cBwYDGwZ1mZI4Br0vM3Ar8oeH2TgatbtH1vIZvi4Z4qy1u2bXWu\nr2Xblup7KbBfej4S+G2R/3+denTbnqaeiW03ORkKjCqdGKqA9UGLJtKNiIXUngCvldtWz/qghZME\nR30TJ7V0Gzuh24Kmnoltq50MLWp90L6JdFu5bfUqZNtqTJzUiW1sqZ46udkhQ3ki3UK2rRUTJ3Wz\nbtvTtPRkaCvWFxHrImJ9en4dMFzS6AbXV097WrVtgypi2+qYOKmt21iEbguaQSe2Ta+PB6h1MrRV\n68szkW6dap1ca+W2Dbq+ArYNBp84qYhtbKuu6p5FHRPbRsS1kqamk6FPAQ2doKp3fbRwIt10cm0K\nsIOk5cBZwIuK2LZ61keLJwkumThpiaS7gQBmko1OFrKNneCTm2Y5dVv3zKzrOWjMcnLQmOXkoDHL\nyUFjlpODxiwnB41ZTg4as5wcNB0k6fWSfi3pRZJeLOk3BV9FbS3gKwI6TNKngW3S4+GI+FyHm2SD\ncNB0mKThZBeOPg28Ofwf0vXcPeu8HcluDd4W2LrDbbE6eE/TYZKuAi4FXgG8PCJO63CTbBBddWtA\nv5H0AeBvEfHdNMPCbZKmRMSCDjfNavCexiwnH9OY5eSgMcvJQWOWk4PGLCcHjVlODhqznBw0Zjk5\naMxy+n8tP+XP4TtiFwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x629c710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "501 0.246965660902\n",
      "[[  8.          10.          10.         ...,  10.          10.           6.        ]\n",
      " [  8.           8.97988234   9.35515231 ...,   8.75571012   7.98251442\n",
      "    6.        ]\n",
      " [  8.           8.56454461   8.91986385 ...,   7.93035755   7.17449777\n",
      "    6.        ]\n",
      " ..., \n",
      " [  8.           7.0655085    6.44524157 ...,   5.45573527   5.67546166\n",
      "    6.        ]\n",
      " [  8.           6.48625159   5.87961569 ...,   5.2801735    5.48888367\n",
      "    6.        ]\n",
      " [  8.           5.           5.         ...,   5.           5.           6.        ]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANAAAAEZCAYAAADmNKu8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXucJUWV57+noVUQ6OrqpnuUFtAGRIZ1ERXwSfkcHg64\nrrP4RNFV1pHHRx0fy+gHmdEd3dUdZYRFHWUF34MPcARRV0sGHRCRVp7KU6CFhq7u6hZotaXP/pGR\n1VF5IzIjX/dm3orv53M/dW9mZGRk3fzdc+JExElRVSKRSDUWjboBkUifiQKKRGoQBRSJ1CAKKBKp\nQRRQJFKDKKBIpAZRQDmIyO0i8vxRt6NpROR0ETnfvH+ciGwWEckp/zsR2XtY7esTnROQiNwhIg+Z\nL3VGRL4lInuMul3DQEROEZFrReQBEblTRL4iIn/e0ukUQFXvUtXd1AwIisgPReQN8wqq7qqqd7TU\njl7TOQGRfLFHq+puwGOA+4B/Gm2TyiEiO1Q45kzgZOAkYCmwH/BN4OhmWxdpFFXt1Au4HXi+9flI\n4Cbr827AeSTCuh34W2vf6cD51ue9gG3AIvP5h8DfAZcDm4HvAJNW+dcCdwD3A6fZbQGeDvwE2Ais\nJRH1jtax24C/Bn4N3Ap8AvhI5touBE51XPM+wJ+Ap+b8X/Ku+3XAvwH/C9hgzn+EtX9vYBrYBFxq\n2n5e9n8EfMC04yHz/znTurYnNNCO15ttm83fV476fqt9v466AXkCAnYG/i9wrrX/POAbZt9ewK+A\nEywBnZcR0MMZAd0MrAYeaT7/D7PvAOB3wLOAxcBHgT9abTkYOAQQYE/geuCUjIAuBZaYup8O3G3t\nXwY8ACx3XPOJwO0F/5e8634d8AfgDaZ9/w1Yax37E3NTLwaeY27g83L+R2/InPthS0CV2mHKbwL2\nMZ9XAk8a9f02rgLabH7B/gjcDfy52bfIfEFPtMq/GfhBCQGdZu1/C3Cxef8+4IvWvp3NuZ7vaeep\nwNesz9uAwzNlrgdeYN6/FfhXT12nAT/J+Z8UXffrgF9b+3Yy7VkBPM78H3ey9n+hpIC2AU+o2Y6d\nzXf6n4BHjfo+a+rVxT4QwLGqOknyS34ycJmIrACWAzsCd1plfwOUCTLca71/CNjFvH8scFe6Q1Uf\nAmbSzyKyrwlo3CMis8AHTXts7s58Pg94jXn/GuB8T5tmSPp7PkKue+66VHWLebsLyXVttLalx1ah\ncjvM//M4kh+te8z/8okV29EZuiogAdCEb5D8Qj4bWE/io+9lld2LpE8C8CDJL11K3k2Z5R6SX+uk\nASI7k7hdKf8HuBFYraoTwN+m7bTITm3/PHCsiDwZ2J8kKODi/wGrRORgz/71wFb8153HPcBSEdnJ\n2rZnTvm86fl12oGqfk9VXwz8GYnr9+mQ47pMVwU0h4gcC0wAN6jqNuArwAdFZBcR2Qt4G9t/2dcA\nzzVjG0uA95Q41QXAS0TkmSKymCTYYAtkV2Czqj4kIvuT/JLmoqprgZ+Z9n1NVf/gKXcLcDbwJRE5\nXEQWi8gjReQ4EXmXue5/ybnuvDbcadpwhqn32cBfZorZ17mOxF1z1bUN+GqVdojIChE5xvwwbSXp\nDz5cdFzX6aqAvmXGgTYBfw8cr6o3mX2nkLhetwGXAZ9X1XMBVPX7JAL7JXAV8K1Mvd5fV1W9gaSf\n8iXgtyRule2S/Q3wahHZDHwS+HJg3Z8DDiRx57yo6qkkkbuzSCJ9twAvta7hZDzX7avSev9q4DBz\nTe8zbfKV/TjwV2YM7mOO/d7/f0E7FgFvJ7FW64HnEvAj1HXEdPbaPYnIIpJfwbtV9RjH/jNJwtUP\nAq9X1TWtN2pIiMhzSELre4+6LZHmGZYFOhW4wbVDRI4k6VfsSxLOPWdIbWod4wqeyhj4+hE3rQtI\nRFYBRwH/7ClyLMa9UdUrgSUisrLtdrWN6SdtJBnv+PiImxNpiR2HcI5/BN5JMsDoYg+s8DGJj7wH\nSWe2t5g+2y6FBSO9plULJCJHA+tMn0YYDPtGIr2mbQv0LOAYETmKZFR6VxE5T1WPt8qsxRp/AVbh\nGFcQkZg+qOeoaq0f0EkR3Rhe/DfDCNwMJQoHICKHA+/IRuGMuN6qqkeLyGHAx1T1MMfxum0Gpief\nwTU8hVtZzSwTbGSCWZYyy8Rc2XnvN09kq5pjy3r/Pj7+fjj1/aGXl3B/jfvj0++HN5U8X1l2t77r\nKteXYaflswPbJnZLtk0wa14bWcosX5Q31haQiOiHA8u+m/qCDWEYfaABROREkokGn1LVi0XkKBG5\nhSSMfYL3uA2w9+QdbDQCuZXVc/sm2DgnpAmSL3GWibkvFAbFlN4AuUJyUUcoo8Ru94Oy/fPu5X9E\ni8STspTBcuPE0ASkqj8CfmTefzKz76SgSmZgYp9EJrPMsppbmWWCGZbNszpZ0n3pF+wSUrCI+iqe\nPEoKKSse+0cqFU9qfcadkVigukwwy+O5ndt5vLfMLEuZYHbAIkHyhReK6NCp+RXWEc5McRH2mQor\nV8Sy4iIAHDw1uO1+KRSRy/Kk2OJpixWt1VyN3gloYsMWJibTL8r9C7eRiTmXLhVRUn42XESHTVVr\nYFURPLni+ULPnxXWUz3nCxCRjW19FiJdnQuXS/oLtzTjJixjZi6kkJTbOK/8wPvQL7/I+sxYr65S\npo2e6w113ZK/4+++QQ8tkGwAJhmwLMDA56qWqBRlRLO++mmCyK5O8pG2OdTdi3jplwXakPyZ2JCs\n09oewN44FzZdZt3ReZaoNiG/5uszr7Ype76GLGabfZ4mEZFTTdaja0XkFE+ZKRG5RkSuE5EfFtXZ\nOwuUYluRLKmIZplgKbNzYe+Q4weCCS53Ju/GqyqUkONCLYyrTt+xM7gtUUFfaFR9n6qTJE16sDcC\nTyNZlHmJiPyrqt5mlVlCspzkxaq6VkQK/+O9FRD4f/mcrprDlZurp64bB+HCqWOJ8o4t+qrTY13l\nfCIqSbb/07ExoCcBV6aLGkXkMuBlwEesMq8iWfi4FkBVC7+tfrlwEOx2uFy51s6f928elhsXeg7f\n/i4HQJrhOuA5IrLUrIo9ivlTyCDJxTdpkkteJSKvLaq0lxYoDSS4cAUUfFYozw10UkY8ZcRSpmyI\nG1fkuq0PrKdkSLvLqOpNIvJh4Hsky8mvYXBJ+Y4k6cueDzwa+HcR+Xez5N5JLwUEJpAwuV0oS+eE\ns3R4jXDd+FUtQJ1z5onB57q5RNSQK9cmvj7Qz4CrC441S8/PBRCRDzJ/GQ0kS/jXq+rvgd8bN+8/\nkiyvd9I/Fy5DGn1zUduNyxv/KSOetl24kPpD29tTV+5pJMuZ05cLEdnd/N2TJD/dFzNFLgSeLSI7\nGDfvUJJMTF56a4GK8LlnqRtXmpCQdcg2F/eXaMfuge0Idd1C3bnx4GsiMkmSFeivVXVzZmLzTSJy\nKUlSmoeBT5lkM17GUkBZ8ZTu6xSRFUZZ8ZQRTNGxPkGVcd3yRGT1g7asn8idC9d1VPW5jm3Zic0f\nYX5kLpfeu3BZQqeQ+Ny+0ksbXGJyied+69UkRfVWsYw9deNGwVhZINeg6TJmmKnbM7ZvqPWe967P\nEC6YkJu26DLSc2WtkssaZa1OT1y5rmWbGRsLFDKdpNWBvbLimaH8JNTQY3wWKUTwkVKMjYAaxxWB\n81mfUPE0PWu7qL4QEfn22XU6/he1Z26MCWPhwvkjbg0ED6rc7NkbN7SOOlN1fDOs7yc/ctcT162r\njIWAhkKo9SkjnjqzFcpMDs2KKK//03FBrfBlF8yyqdVmzNF2XrhHisiVZnr4tSJyuqPM4SIyKyI/\nN6/3ljnHSKfSVxFPU4OqefW43Lps2youdygdpRxzWrVAqvoHEXmeeSTIDsCPReQSVf1ppuhlrqTz\nncR34xWJp6q1KTP3rWhqTp4713HL01Vad+HMk8kgedrcjrgfA9KtVDdppznU/aojnipz53w3esj8\nNltEPtG4to/RxNImaV1A5tEmV5M82PcsVb3KUewZIrKGJCPpO4umT7SF1z1Zn/mbR6h4mlgXFHrz\n500SLbI8HZtguji0LePQB4LkqWaq+hSSlL2HisgBmSJXA3uq6kEkD5jyPQaxO/isT4h4mpxU6qvL\ntX1eWDqnvkgphplYcbNZY34E1rOCVPUB6/0lInK2iEyq6oZsHe//PEmG7Z3h8BfCQTV6TY3OjXNR\ndtAyZMZCE/Pe7PPluXIBfaLZzUnmV3ttVRrUmZm+jrumixYY9J9WBWTWlG9V1U3mIbcvAj6UKbNS\nVdeZ94eQ5OseEA/A+19DspBuGegkjcXfcgcF89w3n/UJFU/ZeXF2eZeYmnDHAoRTNKl0lqUsmzqQ\nfaZWzW277oyLCk7cT9q2QI8BPmf6QYuAr5hc2HNTyIGXi8hbSKaYbyF5FPro8QUQitbQNDk/Lo+8\neW95c9xsERUNsmaPiYGEAdoOY19LskQ2u/2T1vuzSDKhdIPQFL4hImhLPNn6yorIRV/WCnmW8g9w\nW3GRJlhQMxFc6a1qEzpTu2hiaSjOFFTmb+hsgypWKOIkTiY1OEPY2f5P+rfsJE2feKpMLC0zebRq\nVC1G44JZ0AJqbaJp3iBrelzdGdlVROSbbe3a3wCtWPyOsaBcuGDK3txFSxvq1B1SV9kZ2Fn65MbV\nGNQVkbeRZCfdBlwLnKCqf7T2/w3wapLZMotJkjEuV1VvyHGsLVAreRDKBAHKzMzWrcUvH2USnhRZ\noWy59G9a1gRZ+japVEQeC5wMHKyqTyYxHq+wy6jqR1T1Kap6MPDfgek88cCYC8hFbkaeokV0NkUu\nUKh4isQRWrZo9nUo493/2QF4tIjsCOwM/Dan7CuBLxVVuOAElMU7iBp6I1VcFlBKOKHH5omoygra\nAPqyMlVVfwt8FLiTZM7lrKp+31XWDPofAXytqN7YB2qK0PlxUF04rnpk8fxtZSd/dmyyaCGecaDp\ndTB9n/8wEZkAjgX2IplqeoGIvEpVs8kVAf4SuLzIfYMFICBXRp7SfaNs/2cms92m6Jc8VzzrPNtz\nctG4RJRtT9HyhTFgamXySjnjuoEiLwRuS6eJicjXgWcymJ0Ukr5RofsGY+DCNREomNchrhIlC+0n\necWzDr947P2eMtl6Q64hZI1SiWUcrU/Orc+dwGEi8igREeAFONL2mmcEHU6S5reQ3guoVep0qIPE\nUyQcFxVE5OsLZSmynj1OuGhWQV9A8lSGX5jNnxKRE0XkzVbRlwKXquqWkHrH3oVrnCrum1c8VVlH\noykG+9YPqoiqngGckdmcTe37OeBzoXUuCAG53Itcl6OMOIpmZzupI55sHZaQsv2hojlvVftEo5yV\n3TGhj4WAGvG/fbOwQ924YOtTJJ5sKGlFYAOs8+UFFWCsgwnDJvaBhkqeeO5jUDz2dl+MtkCQZfst\n6Q/BeA+oNsaCEZBrYmPlQcC86JR9wwaN9+SJw1XWRUZEvvPmzSIvs4DQ0LfpPG0wFi5cY+StLM3u\ntwkaxXdZilDhZI9xuXUNBxYqYudFaIXQBXVDYmwsUKvjEFUeUVJofaqIp8Sx9vmbWLoQXTonYyOg\nEIb6AOJ5ZK1PHfHk1REQ3Wt6SfkCZ+S5sU25M0XkZhFZIyIHtdmmLJX9+Lxf5Npz3UIHWJsQYkl6\nPJjaBq0KSFX/ADzPJFY8CDjSpK6aQ0SOBFar6r4kD1g+p+r5GlkBWTRruUyuA6Cc9ckKp2AKj7M+\nq2yRG5cysOYnp2xkHl3IjX0scJ4pe6WILLFzxXWGol/eoF/miqHouTINBwkaGg9KEywOhY4NpLbe\nBxKRRSJyDXAv8D1Hbuw9gLusz2vNtlZwBhtCU1mVJlQYdesr6cqFWpgYOChkGBZoG/AUEdkN+KaI\nHFA1eXxIat/RBQoM3v5PxU6/85giS2SVCZmZ0AIz09exYfp6HsXvh37uYTLy3NgkFudx1udVZtsA\nbaX2HaAoq07IvlZxicg3PmRRddJohcedLJs6kGVTB8492HlcU/u2HYVbbtZXpMtkXwTclCl2EXC8\nKXMYyVLbVvs/tceM8p58MEfRJXSri9cbJgNfQ2LkubHN56NE5BbgQeCEuietJJC6i9CghPvW1Gzs\nPFdu9G7cQmDkubHN55PabMfCIcCNy5K6Z+nfBbI2qCniXDgfrUWgipZuu/BZmhqh7SrJFOMyiAHG\ndipPkRtXaiZ26dF3WwihIeYqwqpA6LV0dTB1WeDLgYjsZ2bF/Nz83SQip2TK7CYiF5lZMdeKyOvz\nmjO2AuoXdcaL8mY6xECFjar+2so8+lSSPvc3MsXeClxvHjn6POCjJhGjkyggmzJuWzbFVRB1b+go\niAZ5IXCrqt6V2a7Arub9rsCMqv7JV8lY94FcOeFq4xNZ5QmkQxJFjMRlOQ537rdPABeJyG+BXSh4\nYuJYC6iI9ldUFvV/mpqJ0OJiup4EDqZ/mbxCEJHFwDHAexy7/wK4RlWfLyKrge+JyJPth2HbLGgB\nFVKpIx3aV4lUwjNIOjWVvFLOcOUb3c6RwNWq6vqGTwD+AUBVbxWR24H9gZ+5Khr7PlBjK1U7NbGy\nTGafKFwHeU9e+A1J/wgRWQnsR84TV8deQN2lAzd2Axl4+pZYRER2JhHI161tdnbSDwDPFJFfAt8D\n3pXm03YRXbgsnbI0ZSjRD/LNNgjt73iObz2hSAOY9Wm7Z7bZT42/h6QfFEQUUGu07UZ1IwvP0OnY\nNKPowtVZTJeOAZUOYY/IfWvquUSRORaEgBbC06Ijo2FBCKhxvOHtpixLmWyleXVUJGbeCSb2gYoY\n+s10X+Z93vKEbD9o/AdUNWYmXeiUSVHl2xbpClFAvaCuiDow5jSmLBgBlcrW4xsLqjRGFHrzRkvT\nRxaMgIZLFMNCodUggoisIsk6uhLYBnxaVc/MlEmfiJzON/q6qn6gzXYNMPKoU4jgKuQ7GENmJ3cK\nLBn0jODatB2F+xPwdlVdIyK7AFeLyHdVNZva6jJVdaRJjBSzQGckdIS2k8vfq6przPsHgBtxp+1t\nK7dux4id+XFjaH0gEdmb5AkNVzp2P8Mkcfi2iBwwrDZFCuhqYpEOMZSBVOO+XQCc6ljZdzWwp6o+\nZB518k2SNRgDhOTG7h9lAg4h/aBuuHRt5cYOX981Hn0gTEaTC4DzVfXC7H5bUKp6iYicLSKTrjUY\nQ8uNPTaMLvAQc2M3x2eBG1T1466dZtVf+v4QQPIWMNWh1eeo9o34MK1GaDuM/Szg1cC15hlBCpwG\n7IXJjQ28XETeAmwlsbu5WVAikS7Rdm7sHwM7FJQ5CzirzXa0Slxj0xvMk0L+GTiQZFzyDap6pbV/\ngsRjWk3yY/6GomdZxdnYkV4R7obf49r4ceBiVf0r0zffObP/NJKUVi8TkSeS/LC/MO8scSpPY7jG\neMpkzxkSC9RimickPkdVzwVQ1T+p6uZMsQOAH5j9vwL2FpHcFPxRQCG4JpGObPpPS6LrbTKVYB4P\nrBeRc01y+U+Zh77Z/AJ4GcwFtPYkeWKilyigsSDOcAhgR5JnVZ1lkss/xGBm0g8BS0Xk5yRJ5q8B\nHi6qNBJpnFmWMsHGFup194HWTG9izXTWI5vH3cBdqppmGL0AeLddQFV/B7wh/WyyknqTKkIUUGRM\nOGhqCQdNLZn7fN4Zd8/br6rrROQuEdlPVX8NvID5D7tOo3QPqepWEXkT8CNfTuyUKKChEd2sDnAK\n8AWTXP424AT7eb3Ak0ie6bsNuB54Y1GFUUA+4gj92KGqvwCentlsZyW9AnhimTqjgFLGPwo1FnQt\nx1+MwkUiNYgCikRqEAUUidQgCigSqUEMIkR6Ran8fkMgWqBIpAZRQJFIDaKAIpEaxD5QpFd0La9F\ntECRSA1aFZCIrBKRH4jI9SJyrYic4il3pojcbJIrHtRmmyKRJhl5bmyTTHG1qu4rIocC5wCHtdyu\nSKQRupAb+1iSJzhgMqQssXPFRSJdptACicjJwOdVtdbywpzc2HsAd1mf15ptcQFNZICuBRFCXLiV\nwFVmnfhngUtVVcucpCA3djDjmRt7PNk0vYb7p69oPDd21ygUkKq+V0TeB7wYOAH4hIh8FfiMqt5a\ndHxRbmwSi/M46/Mqs22AmBu7PyyZOoiJqb1ibmxI1rsC95rXn4ClwAUi8j8DDs/NjQ1cBBwPICKH\nAbOqGt23SOOIyBIR+RcRudFEhg91lCkVEQ7pA51KcoOvJ0mL+k6TdGERcDPwrpxjC3Njq+rFInKU\niNwCPEhi5SIRJzMsq3N4bmbSKhHhkD7QJPAyVf2NvVFVt4nIS/IODMmNbcqdFNCOdllOXNY9xliZ\nSV8PSWZSIJsHa15E2FislXkeUaELp6qnZ8Vj7bsxsP2RyKgJyUzqiwh7iVN5fORmRI70kJDMpJUq\njQyFlcShrfr4xoHWTd/Euulf5R1amJmUEhHhlCigyFiwcmp/Vk7tP/c5GzYPyUxKEhF+K/CV0Ihw\nFNBY0L2ZT23kxW6A3MykVSLCUUCRBUNRZlJTplREOAYRQlju2DYwHOGyAkWWocoTtFt66rbrGiOF\nRAu00JDFo25BLWJq33Gj5zdkpB5RQJFIDRaUgCbi/O3tpH24OGBci9gHGmtaCjiMkJiZNJKhzE0e\nUrZ7Y0LjTBTQUOnZzR3du0KigCKRGsQ+ECQd6plRN6IIn/vWM6tWk64lFYkWqBXKdt7Hr7O/UFgw\nAio1udE3raXSdJeFZSEWGgtGQN0nzwrVtVBRxG3Rdm7sz4jIOhH5pWf/4SIya5bY/lxE3ttmeypR\nK4eFi7yb2SWU6N51mbaDCOcC/4RJ1ODhMlUdkxSJTaw6XQHcZ70vOl/e5wbpyGztBRVEUNXLobDz\nIW22AZhL7tcYrY+PrKC+5alxfONWd3zpQh/oGSaJ3bdF5IChn333UlmK55PeaH2Zkd2XdvaIUY8D\nXQ3sqaoPmaR23wT28xXuV25s2xXLunZNuHrdDgzMTF/HhunrY27sNrETzavqJSJytohMquoGV/kq\nubEnmC3nNy+EBIs+Fy20n+M53p7tvmzqQJZNHdh4buzZzfX6QCJyB7AJ2AZsVdVDMvufSNJ3Pxg4\nTVX/d159wxCQ4Onn2FkfReQQQHziGTmNC6vpNFcVLFLal6sRINhpee+WiGwDpnIe1zMDnAy8NKSy\nVgUkIl8EpoBlInIncDrwCEwWFODlIvIWYCuwBTiuzfZk2Wn5LFvW5/yi7Q7cP7TmlKBMroVuu3oj\nQMjp+6vqepIMprlpq1NaFZCqvqpg/1nAWW2dfxkzdZORV8BnWVzbe5BssSPh6wZR4Hsi8jDwKVX9\ndJ3KRh1E6BYhbpqvjCwG3ZrZaAcSmqJK9p/x5+F/u5xtl18eUvRZqnqPiOxOIqQbzXBLJaKAqpK6\nd6Vmcte1QjWEMiYhbK/L/aSXJK+UD7kfXaWq95i/94vIN4BDgMoC6sI4UCsU5T+Y2K1E57d1LzBE\nGL4y2e0N9X/GcDGdiOxsHjeKiDya5KmL1+UdUlRntEA+akXdbKsS6sblWaIGXbQ6Pwb97w+tBL4h\nIkpy739BVb9rp/c1T4j/GbArsM08YO4A37N9x1JApcd+YIiL6poUSgsWZoyn8ajq7SRPis9u/6T1\nfh3zn9CQy1gKqIhKAoNiq+QMJMCgFWor+uaZ/5bX/1me+dt17m996mQpxq4P1MhTAeybyfVLvTyz\nz/mr3XZkrKj+ls9fZw7hGDF2AirCGVwYyc1Q5wZ3HRsw+7qqe9YX6zQCxkpAjS9bsAm9+ea5S/aN\n7rrBq4ioZMTO5b6NYYRtVIyVgIpwuXcDc7mWe97nbYNAgflEFCokX7mKa3+KLEug5Sk1JFCXmcDX\nkBgbAQ0973W2H2TjtUJ5+IS0MmcfDIrHUy5P4MszZfIs1BhH6aowNgIqopZ716jLE7JMu0g0gXWN\nyeyDLjMWAqpqfQZcD/vXNeipdDn7c61QE4lCKvapqv4YxECCk7EQUBGNu3fZsZOgm9IloqpCChCP\nLWCX8FsSxEJ7hEzvB1KrDooOHLe7bh+kswdMy07psWc0eAdWbcrM2A5M7+tz3fLGtLL9n4AB1pEs\npuvYauGxt0DLmgrJFFmZSlYoZQX5FqmGtRqi9VmIjL2AXGTdjHm/pKH9oLxfaG9fCIqDAyscLx+O\nunyuW571KUuchTBHrwU0sWFL/n6HP15pqk/ojWbfpLkigvpTbRzHV426hYamYwh7gF4LqAnmReLs\nX1aXaLJz31xWqLSIqszADhCPz/rk/Rj4+j8BPyBDCx6sD3wNiZHmxjZlzhSRm01yxYGp5lVxfaHZ\nsaCgLz3U+oSU81qIonGfgv154rEpmmURKU3bFuhc4C98O00yxdWqui9wInBOaMWSk/yqipvm7Qel\n+GZg51kqV12yuMDVyg6k5lgoV13Z8xUFN2o+rXuo03g6yKhzYx+LSTyvqlcCS8yKwMbJWptGrY+r\nfJ6IIEBIOfiOzRNPqPUpct+Mm9vDfHCtMOo+0B7AXdbntWZbZVxTdopC2d5+UIpv4NRlhfJElCek\nIjEVlasinjLWp0IAoZXZ8TX7QCKyyDxKZyBVqoi8SkR+YV6Xi8h/KGpOrwZS7dzYU4fC4UdWq2eC\njcyy1Lz3DMSWXeLtG3zNJmdMb0RX3VUskuvGDgkatDw+NDN9HXdNX91chc1xKnADsJtj323Ac1V1\nk4gcAXwaOCyvslELaC3z15+vMtuc2LmxIcmQZ+MSQ3bbUmbZ6Jm54MxUmooh+9eV1ipURJAvpBCK\nhJO2Ie+z6xi7XA1BLZs6kH2mVs19bio3dh1EZBVwFPBB4O3Z/ap6hfXxCgK8oWG4cN7c2MBFwPEA\nInIYMJvmys6lockFdj/I6cZVGR/Jc+dcN+sy/O6dr1xd8ZS1PnOu3vz+T/o/S/+P2/82sKy+Hf4R\neCeDv70u/itwSVGhkebGVtWLReQoEbkFeBA4oeq5QqxPyLy5SlYIwiwRbL/RXTm3qwxUhsww8Imn\nydkJo+Y303DntHe3iBwNrFPVNSIyRU7ONxF5Hsm9+Oyi0440N7Ypc1Ib57b7OTCYJ7uwH5ROLs32\nhaqKCNxPt+3rAAAO5UlEQVRCgvIJ7H2d/qIVtGUGV7OTTLuC73+18xTsP7X98+VnZEs8CzhGRI4i\n6UnvKiLnqerxdiEReTLwKeCInCc4zDHqKNxIsKNDLjfOGaINGcGH/Ll0yz317F7y5Wpb0blCxJPX\n9+l5+FpVT1PVPVX1CcArgB84xLMn8DXgtap6a0i9YyGg7JhO3oyD+e83OrfPke0L5Y2R5IkoVEhl\n8dWT3VZGPL7jMozLAKqInCgibzYf30cSpjpbRK4RkZ8WHT/qKFznmNhtltnNE/5nB/n6QzDozsF8\nlw7mu3HZmzbkyRBl9lcZG8qxPi76uIBOVX8E/Mi8t7OSvgl4U5m6xk5AeYEDux9kh7PTMqX6QuAX\nEbj7Rekx4BZLFasUajnKisdRR577NjQhDTHjTgi9d+HyljT4XDSfGze3LdsX8rly9vusO1U0C2F5\n5hVK3nGu82T7TWXF47A+4+K+NUH/BJQ7ibTcF+sKJsz9rSKi7PvszZs3jpMVhu/lIlQ4NcTjsj6N\nrbfqMf0TUCC+SJtrXlz6pRcK0CUiV0DBZY2y0bOigVEf2eOKhJNtm6utdt1z9bjF47I+fewHNUWv\nBZS3pMGFy41zTXj0WiFwz1LwuXE+IRXNSMh7ufDVW9SeGuJZyKKxGZsggh0AmP9+/oAqbA8muMpl\nAwrOqJwrsJA3cFo0kJpSNKAaMms6JKydG+bOd9vysN231vKUx6w8wyOvj+Mql3Xlgi2Ry6XL+1x1\nQNVFXp1F7clatRzx+KzPQrdEY2OBUtxz4oqtUBrWDrFEYB52a+eSs8d9fBbI3pbd3gS++grHh+ZH\n2qJ4whkLC+QKZbsjRPMti2+hXZElgkx0zr4BXRapjBUqQ4j1ybM4kCueid1mg/s9Cy36ltJfCzTD\nwM3gsj5Zy+Iia4Vc+2xLBLitEbgtErRvhfLqcQ6shlsdKLckvtXnNHWsD9RfAeWQN7Mguy/PlQMG\n3Llk23aXDggTErjFlKXudJ7sOW0cg6JlhePatlCtD4yBgGQD6GR+maL+jU9EMBidS7YNWiPwCAnc\nYkpxTQ8qQ9EYUoBoIGx8Jw6cDiKq/UjTKiKqF5sPqWDSpd3m8+zkTslfc6Nn/26c+7zUuT+dJ2db\nrewx2f3z3m8enHzqnJDa1pOmcyZ9+kLSoQOjoeLxuW9flDeiqrUuXESUJwXerzdK7fOF0HsL1AQ+\nSwR4Xbrk/aBbl5K1SsD8caQ8fAIrkZM6d+JnidkEZaxOq32flLKLD1tmLCwQVLdCrjL2ylWXNfId\n7/3ssEw2TisVSMhgp2/yZxnRJNv97lqReBqzQMsD79f1w7FArQvIpAf6GEnI/DOq+uHM/sOBC0lS\nCgF8XVU/4Khnu4CgtBtnvw8REeS7dNuP8wvJtw2KRVWFolnS+WHo8qKBcKsTBVSlcpFFwK+BFwC/\nBa4CXqGqN1llDgfeoarHFNQVJCCoLyL7fYg1ytaTLRuyvWmKBjnrjOlUcdXGVUBt94EOAW5W1d8A\niMiXSdL53pQpV/tCXdE4Vzg7JFxtl0sHW9O+UVJ2Yt5NZNdj15XF7jPlkSeysqP/xUIqjqINpW/T\nMiLySOAykqxQOwIXqOoZmTLHAH8PbAO2Am9T1R/n1du2gLKpe+8mEVWWZ4jIGpKkiu9U1RuctTkG\nT50Dqhu2zFmhEBEl2zaSFyBwCSndlyempMz27D8+sqtoQylXtn2xZNvTuMWtuCJVVf8gIs9T1YdE\nZAfgxyJyiaraeQ++r6oXAZi0vl8FnpRXbxeicFcDe5oLOxL4JrBflYpsK1RGREBwpM2e/lMkpmy9\nLkLEFUromExZkdRpm2swe1So6kPm7SNJ7n317AfYhcQS5dK2gNYCe1qfB1L3quoD1vtLRORsEZlU\n1YHVPu//OklGL2DqyTA1ZXZYVihERDAojDJCSj/niSlbbxbb8jXNMAXiY+30LaydDsoMNTRMn/xq\nYDVwlqpe5SjzUuAfSOa/H11YZ8tBhB2AX5EEEe4Bfgq8UlVvtMqsTNP5isghwFdVdW9HXarnM99d\ny85AcIS1YXtQAXIGQTO/kr6Qtaus6xd2xjNFYNi/xlXEUebBzL7rzHK2vKOZIIIE3q/qDyKIyG4k\nns5Jvu6CiDwbOF1VX5R3mrYzkz4sIicB32V7GPtGETkRk94XeLmIvIWk07YFOC74BBsYFJEhzxKB\nYxDUeg84LVKyf+mAJXJZHtdN6LJSw6Kxp5UX1BsqqMroVs+OuUxVxVWobhaRHwJHkDypwVXmchF5\ngs8bSunXQGrWAqUEWiKYb40g3wrlWaVk/+Ds7lGHr0NoQsR51+MSUWMWiD8Gln7EvPOJyHJgq3l0\nyU7ApcCH1BocEZHVaUZSETkYuFBVH5et2aYLQYT6ZC1Rpk8E860RkGuRfJ9Dom11wteu41y0ZcXK\nRNCy/xubbB7yjvAY4HOmH7QI+Ip5uIHtDf1nETmeRKVbgP9SVGn/LBD4ZyC73LlM2TIWybfNVcb3\nzKGkvHsdkr98mLUa5YrQkP+TLaJRW6C26KcFco0HwfaccVlrBAMWCRIx2atZZyd3cv4Kh/R58qJt\noZG2qiHtYS4pcLUxOzTQLsWPjxom/RQQ+EUE85MvpmKy+7qBYoLBG6XMzZ0nrKZp8zy2dc0GVJJt\ng+Nr0FlXrlH6KyAYHJXOs0pQWUzgtk5QvX/jqqfpWQB1yUYkU7JjWK7xsi4FTdqk3wLKUiSoIjGZ\nY7IJG12CAr+ooDlhlaXJCJvvB8MVTMlao4UiovESUBbX0EcqKldkf9J9jK8nOoE/sb1vfMpHUzdc\nHQHlCSfd74pM5iVtGXc3rv8CKkrvm72RGxKVy1LZ5IrLR0nROc+b87SKPPKsKfiDKS4RtWt97mup\n3mr0V0ChebHzyvncuJRlnuN9orKoEj9dam7+oiQpzvOVzBNu43NR52FPjSpwT3OfPTtm9FNArpul\nzEyVPIsDxcJqGSnqyzXcriKxDwgsXf1rWaH5kbiNc/vHnf4JKHvTV7mZ8o7xWZ2UAOvTODWsSxNk\nBTbnnnqsUt6g8rjRLwHZN1L2Jm7iJhuFOLpA3phaFqvv5xKSu/omgwhxILU+9k1e152DYpcuZaEK\nzMaad5hapgm2wKTfZWtrJngX6KeAUvIskqtMSkhkri45Sy06SWbK0zwKrkXYLiKbsnMA+0j/BJR+\n0Rsyn+1tRVRx96pYnzKuUc9xiSgJa493f6h/AoJEAFkhpXTNW3BNcO0yPktkW6H0hyFjmWwRpeND\ns2Meyu6XgGaYLxiXFYLmo1ZN9H0c8+86TVkX1Hw3sm8ioonJJCKXhLT3brBhMYhQn1RIOe7b1oYs\n0eKisHYVQvplXSDrgtrWNGuFrHIyCXvvew93pIsYR7huqW36JaBUNA4XziWY+zbVO92KJc0JsZAh\nu56LQ62gb42VLZ4ZBt25DTAxOctSZrmjbmM7TOsCKsqNbcqcCRwJPAi8XlXXeCu0RGTf3D6xVDX4\nK3PqHAsc17ZiSU55879ebIf8XZZoA7BvUmQpW5jYZ7YTFkhEVgHnkXy124BPq+qZnrJPB34CHKeq\nX8+rt1UBmfXnn8DKjS0iF2ZyYx8JrFbVfUXkUOAc4DBnhTPME096gxeJpIyIVpq/3waeVuK4uvys\nA+dblxHVShyYMiuWsF1UtpVKf+DMvr33uYNreErt9m6n8mTSPwFvV9U1IrILcLWIfNe+F2Hunv0Q\nSdKRQrqQG/tYkl8GVPVKEVli54qbxwbYenMinHRn013KtL4fAbnpWBqmi+fL/m9tQdliW5mKaga4\n2eo3boClk1uYOHT0T7FT1XuBe837B0TkRpLU09k87ScDFwBPD6m3C7mxs2XWmm0D2th6BXw/0OrU\n5cEhnKNv5/OVnxOW+W5WbgJug5WXwh6T8PhD76jQwvYQkb2Bg4ArM9sfC7xUVZ9nknwW0qsgwrWb\n4NohnetBhrvypM/ny9Zjf0fv+jzsfcYdDZ2pPsZ9uwA41U4rbfgY8G67eGGFqtrai6Qv8x3r83uA\nd2fKnEPSWUs/3wSsdNSl8dXvVwP30x0lznev4/gdge+QiMdV/23mdTvwOxKX75i8NrVtga4C9hGR\nvUhyY78CeGWmzEXAW4GviMhhwKyr/zOMHF+RbuPKmV6SzwI3qOrHPfU/IX0vIucC30ofd+Jj5Lmx\nTXbIo0TkFhLP4oQ22xRZmIjIs4BXA9eKyDUkVuo0YC+Yy0xqo0H19iUzaSTSRRaNugFZROQIEblJ\nRH4tIu/2lDlTRG4WkTUiclCb5xORw0VkVkR+bl7vrXGuz4jIOhH5ZU6ZJq8t93xNXpupb5WI/EBE\nrheRa0XkFE+5xq5x5LQZRKjQSVwE3EJiVhcDa4D9M2WOBL5t3h8KXNHy+Q4HLmro+p5NEj79pWd/\nY9cWeL7Grs3U92fAQeb9LiTPhmrt++vCq2sWaG7gVVW3AunAq828gVdgiYg4B80bOh808BBkAFW9\nHHITWTd5bSHng4auzZzvXjXTsDQJEaeDlTaNXuOo6ZqAXAOv2S/AN/Da1vnAPARZRL4tIgdUPFeV\n9tS5tlBauTbfYCWjucbW6NVA6oho7CHIHaSVaysYrBwrumaBCh9KbD4/rqBMY+dT1QfUPL1ZVS8B\nFotIW6t3mry2Qtq4NhHZkUQ856vqhY4iQ73GtumagOYGXkXkESQDr9mBrIuA4wHyBl6bOp/tn5v5\nUaI5z8wMQPD3O5q8tsLztXBtUDBYSTvXODI65cLpkAdeQ85HnYcgZxCRLwJTwDIRuRM4HXhEG9cW\ncj4avDZzvsLByqavcdTEgdRIpAZdc+EikV4RBRSJ1CAKKBKpQRRQJFKDKKBIpAZRQJFIDaKAIpEa\nRAFFIjWIAhohIvI0EfmFiDxCRB4tIte1PNs70jBxJsKIEZG/A3Yyr7vUkfo40l2igEaMiCwmmdS6\nBXimxi+kV0QXbvQsJ1n+vCvwqBG3JVKSaIFGjIhcCHwJeDzwWFU9ecRNipSgU8sZFhoi8lrgj6r6\nZfNUgB+LyJSqTo+4aZFAogWKRGoQ+0CRSA2igCKRGkQBRSI1iAKKRGoQBRSJ1CAKKBKpQRRQJFKD\nKKBIpAb/H3bhqz16A5CKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x95a7940>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# To solve you need to setup a 2D mesh. Use numpy's 'np.meshgrid'\n",
    "# example:\n",
    "#=========\n",
    "# x = np.linspace(0,2,nx)\n",
    "# y = np.linspace(0,4,ny)\n",
    "# T = np.ones((ny,nx)) ##create a 1xn vector of 1's\n",
    "# X,Y = np.meshgrid(x,y)\n",
    "# to plot surfaces use plt.contourf\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "###variable declarations\n",
    "xLength = 2;   yLength = 4;\n",
    "nx = 50\n",
    "ny = 100\n",
    "dx = xLength/(nx-1)\n",
    "dy = yLength/(ny-1)\n",
    "# in our discretization formulation we set dx = dy, \n",
    "# therefore our in this section we must set our mesh to dx = dy as well\n",
    "\n",
    "x = np.linspace(0,2,nx)\n",
    "y = np.linspace(0,4,ny)\n",
    "\n",
    "T = np.ones((ny,nx)) ##create a 1xn vector of 1's\n",
    "T_old = T.copy()   # store current values of Temperature\n",
    "\n",
    "X,Y = np.meshgrid(x,y)\n",
    "\n",
    "# ******************************************\n",
    "# Assign boundary conditions\n",
    "# ******************************************\n",
    "T[0,:]=10\n",
    "T[-1,:]=5\n",
    "T[:,0]=8\n",
    "T[:,-1]=6\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "cp = plt.contourf(X, Y, T, 50)\n",
    "plt.colorbar(cp)\n",
    "plt.title('Boundary Conditions')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.axes().set_aspect('equal')\n",
    "plt.show()\n",
    "\n",
    "# ******************************************\n",
    "# ITERATION SETTINGS\n",
    "# ******************************************\n",
    "criteria = 0.00001     # convergence criteria\n",
    "maxIteration = 500    # in case we never reach convergence\n",
    "iteration = 1        # initial iteration counter\n",
    "residual = 1         # set any number larger than criteria to the iteration\n",
    "\n",
    "\n",
    "# ******************************************\n",
    "# SOLVE THE PROBLEM\n",
    "# ******************************************\n",
    "while residual > criteria:         # continuously loop while residual > criteria\n",
    "    \n",
    "    T[1:-1, 1:-1] = (T[1:-1,2:] + T[1:-1,0:-2] +T[2:,1:-1] +T[0:-2,1:-1] )*0.25\n",
    "        \n",
    "    residual = np.sqrt(sum(sum(np.square((T - T_old)))))   # calculate the residual\n",
    "    T_old = T.copy()             # store the Temperature values into Told\n",
    "    iteration +=1\n",
    "    \n",
    "    #print (iteration, residual)\n",
    "    if iteration > maxIteration:\n",
    "        break                      # stop the iteration if the maximum has been reached\n",
    "\n",
    "print(iteration, residual)\n",
    "print(T)\n",
    "\n",
    "\n",
    "cp = plt.contourf(X, Y, T,50)\n",
    "plt.colorbar(cp)\n",
    "plt.title('Boundary Conditions')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.axes().set_aspect('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# np.sum(np.sqrt(sum((T - T_old)/T)**2)) \n",
    "np.sqrt(sum(sum(np.square((T - T_old)))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4 style=\"color:#7171C6;\">NOTES: How to setup other types of boundary conditions</h4>\n",
    "\n",
    "\n",
    "<em>Treatment of boundary nodes - energy balance method</em>\n",
    "* treat all heat transfer into the volume element from all surfaces except for specified heat flux, since its direction is already specified.\n",
    "* examples for on convection boundaries were given in class, which are useful for  investigating fins\n",
    "\n",
    "for a boundary node on a plane surface with convection we get:\n",
    "$$T_{i,j} = \\frac{1}{\\left(2 + \\phi \\right)}\\left[  \\left( T_{i-1,j} + \\frac{T_{i,j+1} + T_{i,j-1}}{2}\\right ) +\n",
    "\\phi T_\\infty \\right]$$\n",
    "\n",
    "where\n",
    "$$\\phi = \\frac{h\\Delta x}{k}$$"
   ]
  }
 ],
 "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.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}