{
 "metadata": {
  "name": "",
  "signature": "sha256:3decae47d4f6b9719477ef0d353ab055e5a195aa36b854cf14deb90fe52caa6a"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Diffusion in 2D\n",
      "==\n",
      "\n",
      "$$\\frac{\\partial u}{\\partial t} = \\nabla\\cdot D\\nabla x$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "from pylab import *\n",
      "from numpy import *"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#dimensions of the square computational domain:\n",
      "maxx = 10.\n",
      "maxt = 1.\n",
      "\n",
      "#difffusivity\n",
      "D = 0.5"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#discretization parameters\n",
      "nx = 100   # number of unknown grid points in spatial direction\n",
      "SC = 1.0  # stability criterion SC = 2*D*dt/dx^2, for FTCS should be SC <= 1\n",
      "\n",
      "def diffusion_init(maxx, maxt, D, nx, SC, stride=1):\n",
      "    #choose time step according to CFL condition\n",
      "    dx = maxx/(nx+1)\n",
      "    dt = SC*dx**2/(4*D)\n",
      "    nt = int(maxt/dt)+1\n",
      "    \n",
      "    #define array for storing the solution\n",
      "    U = zeros((int(nt/stride), nx+2, nx+2))\n",
      "    \n",
      "    x = arange(nx+2)*dx\n",
      "    t = arange(nt)*dt\n",
      "    return U, dx, dt, x, t"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def diffusion_solve(U, dx, dt, nx, nt, D, stride=1):\n",
      "    Ui = U[0,:,:]\n",
      "    for it in range(0,nt-1):\n",
      "        Ui[1:-1, 1:-1] = Ui[1:-1, 1:-1] + D*dt/dx**2*(Ui[2:, 1:-1] - 2*Ui[1:-1, 1:-1] + Ui[0:-2, 1:-1])\\\n",
      "            + D*dt/dx**2*(Ui[1:-1,2:] - 2*Ui[1:-1, 1:-1] + Ui[1:-1, 0:-2])\n",
      "        if it%stride == 0:\n",
      "            U[(it+1)/stride, :, :] = Ui"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's try to propagate a square initial condition"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "U, dx, dt, x, t = diffusion_init(maxx, maxt, D, nx, SC=0.9)\n",
      "X, Y = meshgrid(x, x)\n",
      "U[0, :, :] = 0.\n",
      "U[0, (X<maxx*0.4) & (X>maxx*0.2) & (Y<maxx*0.4) & (Y>maxx*0.2) ] = 1.\n",
      "pcolormesh(U[0,:,:], rasterized=True, vmin=-1, vmax=1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "<matplotlib.collections.QuadMesh at 0x7fd398ed5128>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD3tJREFUeJzt3H+sX3V9x/Hna6040Y2uMSlIm9vO2QhEjWQi2zRetbrO\nmMJfCJmkKvMfdaJ/ONsaBiGBKMaoycIfUyHVDBamhEBCsl6Rb3GZ8UdEhpSusKwKml6c+HOJsch7\nf3wP+PXu9kLP99t7e7+f5yNpcs7nnM+573f6zet+es73NFWFJGm6/d5KFyBJOvEMe0lqgGEvSQ0w\n7CWpAYa9JDXAsJekBiwZ9kluSDKf5P6RsY8neTDJfUluTXLayLHdSR5KcjDJm09k4ZKkZ++ZVvY3\nAtsXjO0DzqmqVwCHgN0ASc4G3gac3c25Pon/cpCkk8CSYVxVXwV+smBsrqqe7Ha/Dmzsti8Abq6q\no1V1GHgYOG+y5UqS+hh35f0u4M5u+0XAoyPHHgXOHPP6kqQJ6B32ST4C/LqqblriNP8vBkk6Cazt\nMynJO4C3AG8cGf4BsGlkf2M3tnCuvwAkqYeqSt+5x72yT7Id+BBwQVX9auTQ7cDFSU5JsgV4CfCN\nxa5RVVP758orr1zxGuzP/lrsb5p7qxp/jbzkyj7JzcDrgBcmeQS4kuG3b04B5pIAfK2q3lNVB5Lc\nAhwAngDeU5OoUJI0tiXDvqouWWT4hiXOvxa4dtyiJEmT5ffgJ2x2dnalSzih7G91m+b+prm3Schy\n32lJ4t0dSTpOSajlfEArSVp9DHtJaoBhL0kNMOwlqQGGvSQ1wLCXpAYY9pLUAMNekhpg2EtSAwx7\nSWqAYS9JDTDsJakBhr0kNcCwl6QGGPaS1ADDXpIaYNhLUgMMe0lqgGEvSQ0w7CWpAYa9JDXAsJek\nBhj2ktQAw16SGmDYS1IDDHtJasCSYZ/khiTzSe4fGVufZC7JoST7kqwbObY7yUNJDiZ584ksXJL0\n7KWqjn0weS3wS+DzVfWybuw64H+q6rokHwb+qKp2JTkbuAl4FXAm8GVga1U9ueCadUXtOTHdSNIi\nruaalS5hbEmoqvSdv+TKvqq+CvxkwfAOYG+3vRe4sNu+ALi5qo5W1WHgYeC8voVJkianzz37DVU1\n323PAxu67RcBj46c9yjDFb4kaYWN9YC2hveAjn0faOljkqRlsrbHnPkkp1fVkSRnAI914z8ANo2c\nt7Eb+3/2X3XP09szszNsnp3pUYYkTa/BYMBgMJjY9ZZ8QAuQZDNwx4IHtD+uqo8l2QWsW/CA9jx+\n+4D2T2rBD/ABraTl5gPaZ1jZJ7kZeB3wwiSPAH8PfBS4JcllwGHgIoCqOpDkFuAA8ATwnoVBL0la\nGUuGfVVdcoxD245x/rXAteMWJUmaLN+glaQGGPaS1ADDXpIaYNhLUgMMe0lqgGEvSQ0w7CWpAYa9\nJDXAsJekBhj2ktQAw16SGmDYS1IDDHtJaoBhL0kNMOwlqQGGvSQ1wLCXpAYY9pLUAMNekhpg2EtS\nAwx7SWqAYS9JDTDsJakBhr0kNcCwl6QGGPaS1ADDXpIaYNhLUgMMe0lqQO+wT7I7yQNJ7k9yU5Ln\nJlmfZC7JoST7kqybZLGSpH56hX2SzcC7gXOr6mXAGuBiYBcwV1Vbgbu6fUnSCuu7sv85cBQ4Ncla\n4FTgh8AOYG93zl7gwrErlCSNrVfYV9XjwCeA7zMM+Z9W1Rywoarmu9PmgQ0TqVKSNJa1fSYleTHw\nAWAz8DPgX5K8ffScqqoktdj8/Vfd8/T2zOwMm2dn+pQhSVNrMBgwGAwmdr1ULZrHS09K3ga8qar+\nptu/FDgfeAPw+qo6kuQM4O6qeumCuXVF7Rm/ckl6lq7mmpUuYWxJqKr0nd/3nv1B4Pwkz0sSYBtw\nALgD2NmdsxO4rW9hkqTJ6XUbp6ruS/J54FvAk8C3gX8E/gC4JcllwGHgognVKUkaQ6+wB6iq64Dr\nFgw/znCVL0k6ifgGrSQ1wLCXpAYY9pLUAMNekhpg2EtSAwx7SWqAYS9JDTDsJakBhr0kNcCwl6QG\nGPaS1ADDXpIaYNhLUgMMe0lqgGEvSQ0w7CWpAYa9JDXAsJekBhj2ktQAw16SGmDYS1IDDHtJaoBh\nL0kNMOwlqQGGvSQ1wLCXpAYY9pLUAMNekhrQO+yTrEvyxSQPJjmQ5NVJ1ieZS3Ioyb4k6yZZrCSp\nn3FW9p8G7qyqs4CXAweBXcBcVW0F7ur2JUkrrFfYJzkNeG1V3QBQVU9U1c+AHcDe7rS9wIUTqVKS\nNJa+K/stwI+S3Jjk20k+k+T5wIaqmu/OmQc2TKRKSdJY1o4x71zgfVX1zSSfYsEtm6qqJLXY5P1X\n3fP09szsDJtnZ3qWIUnTaTAYMBgMJna9VC2ax0tPSk4HvlZVW7r91wC7gT8GXl9VR5KcAdxdVS9d\nMLeuqD3jVy5Jz9LVXLPSJYwtCVWVvvN73capqiPAI0m2dkPbgAeAO4Cd3dhO4La+hUmSJqfvbRyA\nvwX+KckpwH8B7wTWALckuQw4DFw0doWSpLH1Dvuqug941SKHtvUvR5J0IvgGrSQ1wLCXpAYY9pLU\nAMNekhpg2EtSAwx7SWqAYS9JDTDsJakBhr0kNcCwl6QGGPaS1ADDXpIaYNhLUgMMe0lqgGEvSQ0w\n7CWpAYa9JDXAsJekBhj2ktQAw16SGmDYS1IDDHtJaoBhL0kNMOwlqQGGvSQ1wLCXpAYY9pLUAMNe\nkhowVtgnWZPk3iR3dPvrk8wlOZRkX5J1kylTkjSOcVf2lwMHgOr2dwFzVbUVuKvblyStsN5hn2Qj\n8Bbgs0C64R3A3m57L3DhWNVJkiZi7RhzPwl8CPjDkbENVTXfbc8DG8a4/qq3JteudAknzG9qz0qX\nIOk49Ar7JG8FHquqe5PMLnZOVVWSWuzY/qvueXp7ZnaGzbMzfcqQpKk1GAwYDAYTu16qFs3jpScl\n1wKXAk8Av89wdX8r8CpgtqqOJDkDuLuqXrpgbl3RyKrQlb10criaa1a6hLEloaryzGcurtc9+6ra\nU1WbqmoLcDHwlaq6FLgd2NmdthO4rW9hkqTJmdT37J/658FHgTclOQS8oduXJK2wcR7QAlBV+4H9\n3fbjwLZxrylJmizfoJWkBhj2ktQAw16SGmDYS1IDDHtJaoBhL0kNMOwlqQGGvSQ1wLCXpAYY9pLU\nAMNekhpg2EtSAwx7SWqAYS9JDTDsJakBhr0kNcCwl6QGGPaS1ADDXpIaYNhLUgMMe0lqgGEvSQ0w\n7CWpAYa9JDXAsJekBhj2ktQAw16SGmDYS1IDeoV9kk1J7k7yQJLvJnl/N74+yVySQ0n2JVk32XIl\nSX30XdkfBT5YVecA5wPvTXIWsAuYq6qtwF3dviRphfUK+6o6UlXf6bZ/CTwInAnsAPZ2p+0FLpxE\nkZKk8awd9wJJNgOvBL4ObKiq+e7QPLBh3OuvZr+pPStdgiQBYz6gTfIC4EvA5VX1i9FjVVVAjXN9\nSdJk9F7ZJ3kOw6D/QlXd1g3PJzm9qo4kOQN4bLG5+6+65+ntmdkZNs/O9C1DkqbSYDBgMBhM7HoZ\nLsCPc1IShvfkf1xVHxwZv64b+1iSXcC6qtq1YG5d4e0NScvoaq5Z6RLGloSqSt/5fVf2fwG8HfiP\nJPd2Y7uBjwK3JLkMOAxc1LcwSdLk9Ar7qvo3jn2/f1v/ciRJJ4Jv0EpSAwx7SWqAYS9JDTDsJakB\nhr0kNcCwl6QGGPaS1ADDXpIaYNhLUgMMe0lqgGEvSQ0w7CWpAYa9JDXAsJekBhj2ktQAw16SGmDY\nS1IDDHtJaoBhL0kNMOwlqQGGvSQ1wLCXpAYY9pLUAMNekhpg2EtSAwx7SWqAYS9JDTDsJakBEw/7\nJNuTHEzyUJIPT/r6kqTjN9GwT7IG+AdgO3A2cEmSsyb5M052hwffW+kSTij7W92mub9p7m0SJr2y\nPw94uKoOV9VR4J+BCyb8M05q35vyD5z9rW7T3N809zYJkw77M4FHRvYf7cYkSSto0mFfE76eJGkC\nUjW5fE5yPnBVVW3v9ncDT1bVx0bO8ReCJPVQVek7d9Jhvxb4T+CNwA+BbwCXVNWDE/shkqTjtnaS\nF6uqJ5K8D/hXYA3wOYNeklbeRFf2kqST07K+QTttL1wl2ZTk7iQPJPlukvd34+uTzCU5lGRfknUr\nXWtfSdYkuTfJHd3+NPW2LskXkzyY5ECSV09Zf7u7z+b9SW5K8tzV3F+SG5LMJ7l/ZOyY/XT9P9Rl\nzptXpupn7xj9fbz7fN6X5NYkp40cO67+li3sp/SFq6PAB6vqHOB84L1dT7uAuaraCtzV7a9WlwMH\n+O03raapt08Dd1bVWcDLgYNMSX9JNgPvBs6tqpcxvK16Mau7vxsZ5seoRftJcjbwNoZZsx24PsnJ\n/t/DLNbfPuCcqnoFcAjYDf36W87mp+6Fq6o6UlXf6bZ/CTzI8L2CHcDe7rS9wIUrU+F4kmwE3gJ8\nFnjqWwDT0ttpwGur6gYYPm+qqp8xJf0BP2e4GDm1++LEqQy/NLFq+6uqrwI/WTB8rH4uAG6uqqNV\ndRh4mGEGnbQW66+q5qrqyW7368DGbvu4+1vOsJ/qF666ldQrGf6FbKiq+e7QPLBhhcoa1yeBDwFP\njoxNS29bgB8luTHJt5N8JsnzmZL+qupx4BPA9xmG/E+rao4p6W/Esfp5EcOMeco05M27gDu77ePu\nbznDfmqfBCd5AfAl4PKq+sXosRo+AV91vSd5K/BYVd3Lb1f1v2O19tZZC5wLXF9V5wL/y4JbGqu5\nvyQvBj4AbGYYDC9I8vbRc1Zzf4t5Fv2s2l6TfAT4dVXdtMRpS/a3nGH/A2DTyP4mfvc306qU5DkM\ng/4LVXVbNzyf5PTu+BnAYytV3xj+HNiR5L+Bm4E3JPkC09EbDD97j1bVN7v9LzIM/yNT0t+fAv9e\nVT+uqieAW4E/Y3r6e8qxPo8L82ZjN7bqJHkHw9upfz0yfNz9LWfYfwt4SZLNSU5h+HDh9mX8+ROX\nJMDngANV9amRQ7cDO7vtncBtC+ee7KpqT1VtqqotDB/sfaWqLmUKeoPh8xbgkSRbu6FtwAPAHUxB\nfwwfNp+f5Hnd53Qbwwft09LfU471ebwduDjJKUm2AC9h+JLnqpJkO8NbqRdU1a9GDh1/f1W1bH+A\nv2L4hu3DwO7l/NknqJ/XMLyf/R3g3u7PdmA98GWGT8/3AetWutYx+3wdcHu3PTW9Aa8Avgncx3Dl\ne9qU9fd3DH+B3c/w4eVzVnN/DP+F+UPg1wyf/71zqX6APV3WHAT+cqXr79Hfu4CHgO+N5Mv1ffvz\npSpJasDJ/r1TSdIEGPaS1ADDXpIaYNhLUgMMe0lqgGEvSQ0w7CWpAYa9JDXg/wCtJdzDkCpM8AAA\nAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7fd39a970588>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "diffusion_solve(U, dx, dt, nx, len(t), D)\n",
      "pcolormesh(U[:, 30, :], rasterized=True, vmin=-1, vmax=1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "<matplotlib.collections.QuadMesh at 0x7fd398dd44e0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGOxJREFUeJzt3X2sZHV9x/HPp3NZUnzaEJuFhbVLCiSuaQXTUltbmbSE\nrk0D+A8PqQ2txpho1fpHW1cT3W1Somm0NmnwjwqG2EJLfCBLoikr9QLGKNoCossWSFjLIixatWJj\nZe/w7R9z5u7v3p1zZ+bMmTPn4f1KJnvmzHn4/fbO+Z7f/J6OI0IAgHb7uWUnAACweAR7AOgAgj0A\ndADBHgA6gGAPAB1AsAeADtgy2NveZftLtr9t+1u235Wt32/7mO0Hstcbkn322X7M9hHbly86AwCA\nybxVP3vbZ0k6KyIetP1iSf8u6SpJV0t6LiI+umn7PZJulfRrks6R9EVJF0bECwtKPwBgCluW7CPi\nmYh4MFv+iaRHNAzikuQxu1wp6baIOBERRyU9LumS8pILAChi6jp727slXSzpq9mqd9p+yPZNtrdn\n63ZKOpbsdkwnbw4AgCWZKthnVTiflvTurIT/cUnnSbpI0tOSPrLF7szHAABLtjJpA9unSfqMpH+M\niDskKSKeTT7/hKQ7s7dPSdqV7H5utm7zMbkBAEABETGuCn2iSb1xLOkmSYcj4mPJ+rOTzd4o6eFs\n+aCka21vs32epAsk3Z+T4Na+PvjBDy49DeSN/JG/9r3mMalk/zpJb5L0TdsPZOveJ+k62xdpWEXz\nhKS3ZQH8sO3bJR2WtCbp7TFvCgEAc9sy2EfElzW+9P+FLfa5QdINc6YLAFAiRtAuQL/fX3YSFqbN\neZPIX9O1PX/z2HJQ1cJOalO7AwAzsq1YRAMtAKAdCPYA0AEEewDoAII9AHQAwR4AOoBgDwAdQLAH\ngA4g2ANABxDsAaADCPYA0AEEewDoAII9AHQAwR4AOoBgDwAdMPEZtIvyAb1fPQ3W3/e0tr68kq3f\npue3/Hy4/uTyaPt029M3HGP8+U5f32+QrPvZludIj5GfnpPHGG3Ty9k2XU6l6RxnJWe/NfW23G+Q\n86cfZPsNkv3T5bUN608e43lt2/Lz9BiTth19nrfNpM/z0p+Xj1TeMbZat9Xxxh13krxzoLi/0l8v\nOwlLR8keADpgaSX7suSViJelivTklebHbTOpxAmgGxof7IuaJmBiej1uLkCtUY0DAB3Q2ZJ9V6S/\nYNJSd9rwO0vjIYBm4iqvoXE9cKh2AjAPgn0J6tZIDACbEeyBBaF6DHXCtzExaQATADQVwR6lStsW\nKNkC9cHVWHNlNszm9czB8jFFAhaNYD+FLjbA5uWZmwTQTAT7ipV948g7XtGS4uh4lDSBdiHYV2Rc\nUO7iLwYAy0Gwb6BpbhKLKqGX0QDLPDpA9Qj2cyiz8XSZ3T6ZOgFoPyZCA4AOoBi3QMxnA6AuCPYN\nUbQxt7ehjv1kHTl97oFu2bIax/Yu21+y/W3b37L9rmz9mbYP2X7U9l22tyf77LP9mO0jti9fdAba\naEWD9Ved9bS2/gJQb5Pq7E9Iek9EvErSayW9w/YrJb1X0qGIuFDS3dl72d4j6RpJeyTtlXSjbdoF\nAGDJtgzEEfFMRDyYLf9E0iOSzpF0haRbss1ukXRVtnylpNsi4kREHJX0uKRLFpDu2qK0e1JT/i8G\nyW8poK2mLnXb3i3pYklfk7QjIo5nHx2XtCNb3inpWLLbMQ1vDkDrrCW3M6DupirK2H6xpM9IendE\nPGd7/bOICNuxxe5jP7tn/71y9tHu/it0fn/n1IkGqsTUEViW1dVVra6ulnKsicHe9mkaBvpPRcQd\n2erjts+KiGdsny3p2Wz9U5J2Jbufm607xaX7X7+ph0m9f+q3QV7PHAD11O/31e/3198fOHCg8LEm\n9caxpJskHY6IjyUfHZR0fbZ8vaQ7kvXX2t5m+zxJF0i6v3DqoJ4GW3a77A0GY18AkJpUsn+dpDdJ\n+qbtB7J1+yR9SNLttt8i6aikqyUpIg7bvl3SYQ2L6m+PiK2qeNASzJYJ1NuWwT4ivqz80v9lOfvc\nIOmGOdNVW3XoWTJNyX20zaA3W/Adlz8CONB89DWrWFolU4cbB4BuINjnmGV6gjLmpV/maNlJUyfQ\nsAs0H8EetdNj3h6gdAT7isxb+i/awybdb9b6ewDtwbw1ANABlOxrqIpn01IPD3QLwV6LC65tfKA4\n9elAM1GNAwAdQMkepcor+Y+rNlrZsI6vIrBIXGE1N6kXTm9t/MCswcqpf1p65gDdRTUOAHQAJfsG\nyivNj9tmXAkfQPcQCVAYde5Ac3CFzqGNXStRLcY4oCqdCvZVzDI5y4Rmy5z8LDXupkVJHWgXruiK\n5d1wJv1KmKaeftJ+s9Tfp+mk9Dk9bpKoK76ZU6C6pju4saGtCPY1VMUzZOlzD3RLK4N9XerC24JJ\n04Dma2WwryOqggAsEyNoAaADKNk33MrghbHr13qn3seL9syZxqj3DtU8QD0R7GuuaJdLAEgR7EtQ\nxWCtVF5pftw240r485ilsZbpFID64ApcIBplAdQFDbQA0AEEewDoAKpxEnUejDVNPf2k/cquvwfQ\nHAT7GZV5Q0jr9KuYIiE3Hem5szbXNbpQAq1CUQ8AOoCSfcWm+WVQRd/6SQOs6DYJtAtXccP1cu4L\ngwr+smk11CzVPqP9GG0LVIdgX7I6962nsRboLoJ9zRXthVOmWZ5axXTIQD0R7CtS5pQKeVU347ap\nojoHQP1N/C1v+2bbx20/nKzbb/uY7Qey1xuSz/bZfsz2EduXLyrhAIDpTVNx+0lJezetC0kfjYiL\ns9cXJMn2HknXSNqT7XOjbSqH0Rpr6q2/gCaZGIgj4j5JPxzzkcesu1LSbRFxIiKOSnpc0iVzpbDG\nehrUukEWAEbmKXW/0/ZDtm+yvT1bt1PSsWSbY5LOmeMcyPTWTr6q2A9AuxQN9h+XdJ6kiyQ9Lekj\nW2wbBc+BivTW1tZfANqpUF+NiHh2tGz7E5LuzN4+JWlXsum52bpT3LP/Xjm7D+zuv0Ln93cWScrC\nlFE9M8sx6hJoR/PkDHrj66TpWglUZ3V1Vaurq6Ucq1Cwt312RDydvX2jpFFPnYOSbrX9UQ2rby6Q\ndP+4Y1y6//WbgmE9gl2V8iY/q6JvPQOsloebJKbV7/fV7/fX3x84cKDwsSYGe9u3SbpU0sttPynp\ng5L6ti/SsIrmCUlvk6SIOGz7dkmHNYzeb48IqnEWaJnTJaRW1qdAoGM/UEcTr8yIuG7M6pu32P4G\nSTfMkygsX9F5b/KPN/0o3LpoSjqBaVAMK8Esc9zP2hYwruRedARtuh8ja4FuobIWADqAYA8AHUCw\nB4AOqHXNbdpAVuaskbMq49zzPru26AjYMurpq3hq1SIbcGloBWoe7FFvTRtgxeRl6DKCfcmYGK2b\nGF+AuuMbioUpu68+gOI6G+zrVgJPpy+YpX7eOdmICbF14zlOnpv+90A7cWkvUV0mP0uN0pQ3ERqA\nZiLYN1BeaX7cNpNK+AC6gWA/o6qrf8p86AjTJQDdxSW/QOP65+dNa1w3aTrTKp2mTGg2KW30nkHX\n8I1HqZpyMwC6hmDfENPU00/aj/p7oLsI9nNY5hQOADALJkIDgA4g2ANAB1CN03R5NUkl/mWnmfZg\ntA3TIgD1RLCvyKT++Wkf+DL71k9zvtQasRpoJYJ9CYoOtErnw8kzthfONDeD0TbJX3ianjmjNDHo\nCmgX6uwBoAMov6Gwpj28BOgygn3J6jZ1cp0x3z1QHYJ9UxRttE3346+9NPzywbJx+VesjnPYj5Om\ns+1z2xOI0QU00AJAB1Cy1+LmuFmpov6+gkFVGw+bNsrW6+tDCR3IR8keADqgXkWzDktHtG4YSDWu\n5F5wUFW6n9NzJ8sMpgLaqfGXdvrTnW6PzTD6OzWl2oVuoWiDxgd7LF7eIwoBNAfBviLjnj1bxYRn\nsxqliQnRgHYh2DdFwwdV1bkXD9AFXHUJ6vwBtBXBHgvDRGlAfUzsZ2/7ZtvHbT+crDvT9iHbj9q+\ny/b25LN9th+zfcT25YtKODJrOa8lWdFg/QWgPqYZVPVJSXs3rXuvpEMRcaGku7P3sr1H0jWS9mT7\n3Gi7kQO3ehqsv9IANm59um6SlcEL66/Cpgnqcwb9WdM5bf4BLMfEQBwR90n64abVV0i6JVu+RdJV\n2fKVkm6LiBMRcVTS45IuKSepHbKoEnrJvwJ6Wlt/dVFaDADqrui3dEdEHM+Wj0vakS3vlPTVZLtj\nks4peI5GKhr4xj5+UKqkSmbDuYlbQCvNXcUSESEpttpk3nMAAOZTtBx33PZZEfGM7bMlPZutf0rS\nrmS7c7N1p7hn/71ydh/Y3X+Fzu/vLJiU+iq1Drsm/ex5uhRQndXVVa2urpZyrKKX/0FJ10v6cPbv\nHcn6W21/VMPqmwsk3T/uAJfuf/2mYNi9et86j6Dt6oRodBFFnfT7ffX7/fX3Bw4cKHysiZe07dsk\nXSrp5baflPQBSR+SdLvtt0g6KulqSYqIw7Zvl3RYw+j99qyaBwCwRBODfURcl/PRZTnb3yDphnkS\n1UYLexxhXk1RiQXULj2iEGirjv5Yb7hpmgJG2xCbAYgnVQFAJ1Cyr6O1nOUyj5sa8y1IG4+Z7hho\nPoJ9xerYAwdA+7Uy2K818VGFy7wJpOemFA+0UiuDfSsVvWel+xHIK0WffdQJDbQA0AGU7JtuhkbX\nNmLGSWA6XClLlDvT5RKtp4lvBtAqXNIlS5/Q1BssKJpP05g72qbkv3Cap5Ve+tjBGt65AKwj2Jdg\noT1+yjx03rH4FgCtx2VekZkeQ1izvvhp2rs6GybQdI25dNOGOB5mXV+NGdcAdExjgn0dFQlsuSNo\nJ5XmF/nwkjHbpOmkNA80H5dx03W86yWA6TCoCgA6gPIfFqa34WfHtrmPx/QDQHEE+xn1yuwqU0Y9\n/KRtZv0Lj/YjrgKtQrCvu0XNZ08wBzqFYF9HVfReTM/RkW8B8+igy/j2N8Uiu1623FrBnzFF9wPq\nqKOXfw0VDeZ5vwJmiVPcEIDW49JeoN7aWrK8xITMKU17b8Ob6XvYMLIWWC6CPTAD6v3RVHxzpzBL\nqXTStqXMYT/NMUbbFKx29hQNuJTWgeYg2KOTGKCFriHYL9Msg6PKPt+kSdFm/GacHGw2/0hZAOUj\n2Jeg1FG10vgAX7TGJN0vLcy2fIAVJXdgI4J90zHrJYApEBKAEvGLAnVFsK+LZXZs6eDUCUDXcGk3\n0SJnvawAjblA9WoYCtBk9L0H6olgP4eZBls1bbqEDVMkJOtPrzwlAEpAsM9RRnfKlcELp65cy1nO\n22YwZt0sptkv/RZMqP5J8/R8wSQBqB7Bvunoejk3etCgC+YKCbaPSvqxhuXPExFxie0zJf2LpF+U\ndFTS1RHxoznTCQCYw8/NuX9I6kfExRFxSbbuvZIORcSFku7O3ndGbzBYfy3MWvKatM0CVZJXAKWY\nN9hLkje9v0LSLdnyLZKuKuEc7TdNAG/juQFUooyS/Rdtf8P2W7N1OyLieLZ8XNKOOc8BAJjTvM14\nr4uIp23/gqRDto+kH0ZE2I5xO96z/15Zw49291+h8/s750wKALTL6uqqVldXSznWXME+Ip7O/v2e\n7c9JukTScdtnRcQzts+W9Oy4fS/d//pN/dTnr0MY9aoofRbKZVrUFMcVdEBZSf6+ZXbTzHsQOL1q\n0Db9fl/9fn/9/YEDBwofq3A1ju0zbL8kW36RpMslPSzpoKTrs82ul3RH4dTlWFNv/VW1FQ3WX4VM\nUz8+SF6T9lvLeU3ab9w5Zk1nYq7/EwALN0/Jfoekz9keHeefIuIu29+QdLvttyjrejl3KmtumiA3\n9nGERWPjIufGGaUp2a+URykCWKrCwT4inpB00Zj1P5B02TyJAgCUi3GWdbHMZoZpHlcIoNG4tBeo\ncZOfTWFDnpgUDWgMgn3dVf3AcQCtxCXfFHlBP6/xdNRRiQAPQFz+9VRF9U/Ffe4BLBfBvmS9tQoi\n9TRdIUfbLDCQb8hrr72PGBxwmaAF+BZXbZqHlyzTKE0FbxJNfyxh0YF6jN5F3RHs51BqYFtmjEzP\nzTcCaCUu7arUsRQ/rQV1t2z6rwCgSQj2dTcKtFvNYTPO6C+btx9/eaBTynh4CQCg5ijf1dGkKp8y\nJkKruP89VTbAci012Kc9GOo8B33RQMV0CQDqgmocAOiA2lTjpANXlvUQjMIl+Fn2y3s4V9FfAZMa\naGc9xujH1ozfjGVV09C/HZhObYI9coxi6DJvBgAaj2qctprhkYIA2o/yXVVmCbxl9MaZ5RyTvgU0\nygKNR7Av2crghWUnoTJpXp9fYjomGbUHUb+PLqMaBwA6oJYl+9HMgwzEAYBy1DLYz2LjwKzF3xwq\nvwGtbfp3s1meVDXu8wXiZg3UR+ODfZ15XKxLg29en/tJyn54ybjG2pxG2bF5AlB7BPs6qiKgtnwO\n+0kPIaGxFl1DAy0AdEALy3TVqfPkbcuy8f+kvc+lBZqm1sG+6sbXhZr3vlDGw0uK/rW5pwGNV+tg\nX4bG182WGWgrnsMeQH1wyZegNyj4q6PorJdlPLxkUjpmnBZh/f+g4ffWWTS+IIFOoYEWADqgUyX7\nwYbsLmY2l4lPp5qlVD6NolMYj5vDPu/znFJ+259aNejW5YGW49vcFGXMZ0+tA9BZrQz2dXjqlaRy\nGlfLTH7ZPXNaWJoH2qqVwb71ym6gbQgaRIHiWhQKGm6ZwwgaPoRhhJsBkG8hvXFs77V9xPZjtv+y\njGMO1Ft/NdJApwbVteQ1zfoyTTp3aqDx6a/AQCvrr0VZU2/9BbRV6cHedk/S30vaK2mPpOtsv7Ls\n89TZPfdVeLJBzmtBVr+yuGPXwROrTy47CQt1dPU7y07CQrU9f/NYRMn+EkmPR8TRiDgh6Z8lXbmA\n89TWPV9ewEGLlvZL/pVAsG+277Q8GLY9f/NYRLA/R1J6xRzL1kGavfRdtFpllmqjrfZvSX0+0HWL\nqAiNaTdMJzdLZ0scdZfclgx8Gvf5qccYnLJt3gRqPxszI2NaL3x68vlPdcb68nN6ySlpWuklaXjp\nQN87/fs6/NKXa9tLf3ZKmseld6t0TppZM69r6eT53Mf/6UdtImnbSLr8/Zc8pyNnbz/lGM9n/19r\nG/ZbSZZ7hbZN879tzJL0s2Tbk/8X49Kfrsv7/zlNJ3SGfnrK9uOOmyr6/z1+W9oOUD5HTB2bpzug\n/VpJ+yNib/Z+n6QXIuLDyTblnhQAOiIiXGS/RQT7FUn/Kel3JX1X0v2SrouIR0o9EQBgaqVX40TE\nmu0/lfSvGg7Qv4lADwDLVXrJHgBQP5VOcbyIwVbLZHuX7S/Z/rbtb9l+V7b+TNuHbD9q+y7b25ed\n1nnY7tl+wPad2fvW5M/2dtuftv2I7cO2f70t+bO9L/tuPmz7VtunNzlvtm+2fdz2w8m63Pxk+X8s\nizmXLyfV08vJ399k382HbH/W9suSz2bKX2XBvqWDrU5Iek9EvErSayW9I8vTeyUdiogLJd2dvW+y\nd0s6rJM9rdqUv7+T9PmIeKWkX5F0RC3In+3dkt4q6TUR8csaVqleq2bn7ZMaxo/U2PzY3iPpGg1j\nzV5JN9qu+/M7xuXvLkmviohXS3pU0j6pWP6qzHzrBltFxDMR8WC2/BNJj2g4puAKSbdkm90i6arl\npHB+ts+V9PuSPiFp1AugFfnLSkm/HRE3S8P2poj4H7Ujfz/WsDByRtZp4gwNO0w0Nm8RcZ+kH25a\nnZefKyXdFhEnIuKopMc1jEG1NS5/EXEoIl7I3n5N0rnZ8sz5qzLYt3qwVVaSuljDP8iOiDiefXRc\n0o4lJasMfyvpzyW9kKxrS/7Ok/Q925+0/R+2/8H2i9SC/EXEDyR9RNJ/aRjkfxQRh9SCvG2Sl5+d\nGsaYkTbEmzdL+ny2PHP+qgz2rW0Jtv1iSZ+R9O6IeC79LIYt4I3Mu+0/kPRsRDygk6X6DZqcPw17\no71G0o0R8RpJ/6tN1RpNzZ/tX5L0Z5J2axgYXmz7Tek2Tc1bniny09i82n6/pOcj4tYtNtsyf1UG\n+6ck7Ure79LGO1Mj2T5Nw0D/qYi4I1t93PZZ2ednS3p2Wemb029KusL2E5Juk/Q7tj+l9uTvmKRj\nEfH17P2nNQz+z7Qgf78q6SsR8d8RsSbps5J+Q+3IWyrvu7g53pybrWsc23+sYVXqHyarZ85flcH+\nG5IusL3b9jYNGxcOVnj+0tm2pJskHY6IjyUfHZR0fbZ8vaQ7Nu/bBBHxvojYFRHnadi4928R8Udq\nT/6ekfSk7QuzVZdJ+rakO9X8/B2R9FrbP599Ty/TsJG9DXlL5X0XD0q61vY22+dJukDDAZ6NYnuv\nhtWoV0bE/yUfzZ6/iKjsJekNGo6ufVzSvirPvaD8/JaGddkPSnoge+2VdKakL2rYen6XpO3LTmsJ\neb1U0sFsuTX5k/RqSV+X9JCGpd+XtSV/kv5Cw5vXwxo2Xp7W5Lxp+Ovyu5Ke17D970+2yo+k92Wx\n5oik31t2+gvk782SHpP0nSS+3Fg0fwyqAoAOqHu/UwBACQj2ANABBHsA6ACCPQB0AMEeADqAYA8A\nHUCwB4AOINgDQAf8P2jRXmK1fuqYAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7fd398e44710>"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The time step limit seems to be to strict, let's try increasing the time step twice"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "U, dx, dt, x, t = diffusion_init(maxx, maxt, D, nx, SC=1.1)\n",
      "X, Y = meshgrid(x, x)\n",
      "U[0, :, :] = 0.\n",
      "U[0, (X<maxx*0.4) & (X>maxx*0.2) & (Y<maxx*0.4) & (Y>maxx*0.2) ] = 1."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "diffusion_solve(U, dx, dt, nx, len(t), D)\n",
      "pcolormesh(U[:, 30, :], rasterized=True, vmin=-1, vmax=1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "<matplotlib.collections.QuadMesh at 0x7fd398db1b38>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXm8HUWV+L/HXAgEWWQLBIJhDGGIogQYQJYkSAQEAVkU\nUARBGPiAgLixKe89BvgBjiAoOowKMsimbMKIsoyERVlkHTAgMBAgLGHfTIQk1O+P7q46N6/73Xvf\nu3277+3z/Xze551Xr6tPVd3u0+eeqjotzjkMwzCM3uYDRTfAMAzDyB8z9oZhGBXAjL1hGEYFMGNv\nGIZRAczYG4ZhVAAz9oZhGBVgSGMvIuNF5GYR+auIPCwiR8TlK4rIjSLymIjcICIrqDrHisjjIvKo\niGybdwcMwzCMxshQ6+xFZDVgNefcAyLyQeBe4HPA/sArzrnTReRo4EPOuWNEZDJwMfAvwBrATcAk\n59z7eXfEMAzDyGZIz94596Jz7oFYfgd4hMiI7wxcEB92AdEDAGAX4BLn3ALn3GzgCWCTHNptGIZh\ntEDTMXsRmQBMAe4Cxjrn5sb/mguMjeVxwBxVbQ7Rw8EwDMMokKaMfRzCuQI40jn3tv6fi+JAQ+Vc\nsHwMhmEYBVNrdICILEFk6C90zl0dF88VkdWccy+KyOrAS3H5c8B4VX3NuGzxc9oDwDAMYxg452Q4\n9RqtxhHgF8As59wP1b+uAfaL5f2Aq1X5XiKypIisDawD3J3R4J796evrK7wN1jfrn/Wv935GQiPP\nfgtgH+B/ReT+uOxY4FTg1yLyVWA28IXYgM8SkV8Ds4CFwKFupC00DMMwRsyQxt45dzvZ3v+MjDqn\nAKeMsF2GYRhGG7EdtDkwffr0opuQG73cN7D+dTu93r+RMOSmqtyUilh0xzAMo0VEBJfHBK1hGIbR\nG5ixNwzDqABm7A3DMCqAGXvDMIwKYMbeMAyjApixNwzDqABm7A3DMCqAGXvDMIwKYMbeMAyjAjRM\ncZwXJ3A8oySk0Omnf5CcVtaMXOZ6eerQVH0srE/tq1fmtmXVW+SOQ3MiJ1N1zLM3DMOoALUilZfZ\nM+hGL0lT9bGwPrWvXpnbllXve7yHUU+tSOVlvli68cbRVH0srE/tq1fmtmXVW0R9GMcwY9/xennq\n0FR9LKxP7atX5rZl1TPPfjC1IpWX+WLpxhtHU/WxsD61r16Z25ZVzzz7wdSKVF7mi6UbbxxN1cfC\n+tS+emVuW1Y98+wHUytSeZkvlm68cTRVHwvrU/vqlbltWfXMsx+MLb00DMOoALUilZfZM+hGL0lT\n9bGwPrWvXpnbVl/vRl/2PW7GqKdWpPLyXSz518tTh6bqY2F9al+9MrdNy26HT/uyEyyMM4hakcrL\ndrF0+42jqfpYWJ/aV6/MbdOybOF82fc4HqOeWpHKy3axdPuNo6n6WFif2levzG2L5F8B4HYTX2ae\n/WBqRSovz8XSuXp56tBUfSysT+2rV+a2Abgt9gFg2bVe9mVHcSZGPbUilZflYumVG0dT9bGwPrWv\nXpnbBiBTo/DNfWMm+7Kr2BWjnlqRystysfTKjaOp+lhYn9pXr5xte8rLbs8ofLM1v/dlW3EbRj21\nIpV350U2snp56tBUfSysT+2rV8a2ua+s7eUxa78GwHkc4MtmEbx8I6JWpPJuvMjKfONoqj4W1qf2\n1StL236gZPlkWHlz33KRYT+To3zZBGZj1GM7aA3DMCpArUjl3eJRdIuXpKn6WFif2levLG375oeC\nN+92Cssst+IGAPbiMl82l7EY9dSKVN4tF1m33Diaqo+F9al99crSNndMMPBjV3/ay2dxJAC3sZUv\nW4lXMeqpFam8Wy6ybrlxNFUfC+tT++oV2banlCwbBM/+DqZ4+UccDsD6POTL5jEGo55akcrLfJF1\n442jqfpYWJ/aV6/Itq09Ixj4N7cZ7eXP8xsv78w1ADzFBF82hvkY9dSKVF7mi6wbbxxN1cfC+tS+\nep1v29m+zB0cQjcTRz3s5QH6vPxnNgdgdZ73Ze8RHgxGRK1I5eW7yPKvl6cOTdXHwvrUvnqdbpvb\n4AhfJpODZ3+DislfyL5enswsAN5gBV+2pL2pahC1IpWX7SLr9htHU/WxsD61r16n2yZbBwP/9OSw\nquZr/NjL23CTl59lPADL8o4vW8QojHpqRSov20XW7TeOpupjYX1qX73Ote1ZANyXQuhmCnd4+XB+\n5OUH2MDLqzIXgPlqUnYUizDqsU1VhmEYFaBWpPLyeBSdq5enDk3Vx8L61L56nWqb+0oUjhm99pu+\n7Ep287LOZDmRJ7ycxOotTj80tSKVl+Ui65UbR1P1sbA+ta9enjp+puQk383DK070ZX0MeHlz/uzl\n51ndy0ms/j2W9GUWxhlMrUjlZbiQO10vTx2aqo+F9al99fLUcdDKg1MgJOkPAPblQi8/xPpe1jtk\nEyNvBn5oakUqL8OF3Ol6eerQVH0srE/tq5enDvftwSkQfszXfNlNbOPl8fEELsA7LOvlJHyjV+CY\n4R9MrUjlZbiQO10vTx2aqo+F9al99dqt40Ulp6VASNIfAGzAA16ey6pe1jtkEyNvBn5oakUqL8OF\n3Ol6eerQVH0srE/tq9duHat9Jj0Fwm5cCcCuXOXLnmCil1fgDS+nxefNsx+aWpHKy3Ahd7penjo0\nVR8L61P76rVHx/ledgcNnQIhSX8A9SkQ0kI3YJ59s9SKVF6GC7nT9fLUoan6WFif2levHTrcZvt7\nWdYdOgWCzl75Kit5Oc3Ag3n2zVJrdICInAfsCLzknFs/LusHDgRejg87zjn3+/h/xwIHAIuAI5xz\nNww6aUwZLuRO18tTh6bqY2F9al+9duiQLZtPgTCb8H7ZRqEbMM++WZrZQXs+sP1iZQ44wzk3Jf5J\nDP1kYE9gclznJyJiu3QNwzAKptboAOfcbSIyIeVfklK2C3CJc24BMFtEngA2Ae5MO3cZvJZO18tT\nh6bqY2F9al+94esIa+HdXs3nu0ly3UB2vhsL47RObQR1DxeRfYF7gG86594AxlFv2OcAa2SdoAwX\ncqfr5alDU/WxsD61r95wdbhDQry9lRQIWamKs4y5hXGaozbMej8FTozlfwN+AHw141iXVnhL/63c\nyrT4rwlddwMMt16eOjRVHwvrU/vqtarj0liWKeHWf3DFdb18Iid4OS0Fgk5V3ChOr8t70bOfOXMm\nM2fObMu5xLlUW1x/UBTGuTaZoM36n4gcA+CcOzX+3x+APufcXYvVcd9zxzFKTvFl3XYDDLdenjo0\nVR8L61P76rWsY3xkU9xdIXSz1erNp0BoxcDr8qEM/ImcnPm/bkJEcM6lhdAbUhumwtWdcy/Ef+4K\nfq3UNcDFInIGUfhmHeDurPOU4ULudL08dWiqPhbWp/bVa1WH+1Zki5L0B5CdAmFtZns5LXtlIwOv\ny3vRs28ntUYHiMglwDRgZRF5FugDpovIBkQhmqeAgwGcc7NE5NfALGAhcKgb4qtDGS7kTtfLU4em\n6mNhfWpfvWaOnbdMkOVj0S2fpD+A1lIgtGLgdbkZ+KGpNTrAObd3SvF5Qxx/CnBK1v81ZbiQO10v\nTx2aqo+F9al99Zo5dsy2KgXCtCgFQpL+AFpLgdCKgdfl5tkPTa1I5WW4kDtdL08dmqqPhfWpffWy\n/3+Fl93+IYy87qgHATjBr+FoLQWCefb5UCtSeRku5E7Xy1OHpupjYX1qX72s/7tpu3s5LQVCkv4A\nWkuBYJ59PtjuVsMwjApQK1J5GbyWTtfLU4em6mNhfWpfvaz/y+Yq382kwfluklw30Fq+Gwvj5EOt\nSOVluJA7XS9PHZqqj4X1qX316stC2MXtOXQKhCT9AbSWAsHCOPlQK1J5GS7kTtfLU4em6mNhfWpf\nPV3mjgye+OjxQ6dASNIfQGspEMyzz4dakcrLcCF3ul6eOjRVHwvrU/vq3ajKZMPmUyAk6Q+gtRQI\n5tnnQ61I5WW4kDtdL08dmqqPhfWpffU+PTEYePcZlQKBoVMgrKSyXrYjx4159iOjVqTyMlzIna6X\npw5N1cfC+tS+eu7rwcCPW+VJL5/FkV5OS4HQjuyV5tm3j1qRystwIXe6Xp46NFUfC+vTyOu5VSNZ\n1g+e/e1s7OWfcKiX01IgtDt7pXn2I6NWpPJuvAHKbAw0VR8L69PI68k2kZF/bWpYPbMXl3p5Z67x\ncloKhHZkrzTPvn3YpirDMIwKUCtSeTd6O2X2/DRVHwvr03Dr3e5l9+UoVr8uD/qyVvLdtDtVsYVx\nRkatSOXdcwO0r16eOjRVHwvr0/DquR229LJ8JArjJLluoLV8N+0w2hbGaR+1IpV3yw3QLcZAU/Wx\nsD4Nr55sGiZjn5w0DoAjOcuXtZICod1G2zz7kVErUnm33ADdYgw0VR8L61Mr9dQ6+t3CMsuN45DO\nofzEl2WlQMgrVbF59u2jVqTyct8A+dTLU4em6mNhfWq+njs6GPgxa73m5UvZC4Br2NmXtZICwTz7\nclErUnmZb4BuNAaaqo+F9Wlo+V5VJh9XKRCWG5wCIUl/AK2lQDDPvlzUilRethug242BpupjYX0a\nWt5o42Dg398+ePZTU1IgJOkPoLUUCObZl4takcrLdgN0uzHQVH0srE9Z8mkAuENUCoQVh06BkKQ/\ngM5nrzTPvn3UilRenhugc/Xy1KGp+lhYn9Jlt87RAMh6zadASNIfQLHZK82zHxm2g9YwDKMC1IpU\nXhZvp1c8P03Vx8L6lC7LlpFH//LmYankl7jIy2n5bhq9RhAsjNMN1IpUXpYboFeMgabqY2F90vJj\nXnb7tJ4CIW0NPXTeaFsYZ2TUilTenTfOyOrlqUNT9bGwPgXZ7T3JyzJxAQA3sLUva5QCocgcN+bZ\nt49akcq78cYpszHQVH0sqt6nnylZ1DLLx9caD8A3OMOXNUqBUJZUxebZj4xakcq75cbpFmOgqfpY\nVL1PB31IpUDYeegUCPeykZfLlr3SPPv2UStSebfcON1iDDRVH4uq98kdFQz8suNe9vJFfAloLQVC\nWYy2efYjo1ak8m65cbrFGGiqPhZV7NOLSpaPBs/+vjGTvXwyxwOtpUAo0mibZ98+akUqL/ON043G\nQFP1sahin1abEQz8ghnBs996hCkQymK0zbMfGbUilZf5xulGY6Cp+lhUp09X+DL31WDgxy/3uJfP\n4BteHk4KBPPsewPbQWsYhlEBakUqL5+XlH+9PHVoqj4WVemT22x3XyaTGue72Yh7ge5MVWxhnJFR\nK1J52W6cbjcGmqqPRVX6JJ8MBv6FDUM4Zj/+y8sjTYFgYZzeoFak8rLdON1uDDRVH4ve71MUW3d7\nhjj9ZO7z8vGc7OWRpkAoi9E2z35k1IpUXp4bp3P18tShqfpY9Hqf3CGRNy7j/+HLfs/2Xr6YL3p5\npCkQzLPvDWpFKi/LjdMrxkBT9bHoxT7drGSZEoVvHhk3wZd9h+97uZ0pEMpitM2zHxm1IpWX4ebs\ndL08dWiqPha92KetJ6oUCJ+JwjebcosvO5hzvfwAG3h5VeZ6eTgpEMyz7w1qRSovw83Z6Xp56tBU\nfSx6p08nedkdHuLzK68xB4DzOMCX5ZUCoSxG2zz7kVErUnkZbs5O18tTh6bqY9ErfXJrfdfL8s/B\ns7/7A9EO2NM4xpfllQLBPPvewDZVGYZhVIBakcrL4Il1ul6eOjRVH4te6ZNsEbz5d6YFz/Wz/A6A\nL3KxL8sr301ZPHQL44yMWpHKy3Bzdrpenjo0VR+L7u7TLC+7L4c4/YTRj3j5VI4F4Gam+7K88t1Y\nGKc3qBWpvAw3Z6fr5alDU/Wx6OY+ud1CSmL5cPDsb2FTL5/LwQBswAO+bC6rermdKRDKYrTNsx8Z\ntSKVl+Hm7HS9PHVoqj4W3dan81WZbBQM/NOTx3r5EH7q5R25DgjpDyC/FAjm2fcGtSKVl+Hm7HS9\nPHVoqj4W3dan/ZdNf43g+tzt5W/x716+O/by014jCO1NgVAWo22e/cioFam8DDdnp+vlqUNT9bHo\ntj7p1wiOHveml3/LLl6+jD29PDmO6yfpDyC/jVLm2fcGtUYHiMh5wI7AS8659eOyFYHLgA8Ds4Ev\nOOfeiP93LHAAsAg4wjl3Q9p5oRw3Z6fr5alDU/Wx6JY+vVqLZJkcPPuHV5zo5e+qTVXTudnLSQqE\nVkI3urwbjbZ59iOj1sQx5wM/ApUzFY4BbnTOnS4iR8d/HyMik4E9gcnAGsBNIjLJOfd+2onLcHN2\nul6eOjRVH4tu6dNK0yMj//6ng2e/BX/08oH83MtpKRBaCd3o8m4x2ubZt49aowOcc7eJyITFincG\npsXyBcBMIoO/C3CJc24BMFtEngA2Ae5MO3cZbs5O18tTh6bqY1HuPt3oZbdfZOTHrfikL/sxh3n5\nOnb0cloKhFYMvC7vRqNtnv3IGO4O2rHOuSS70lwgWTIwDpijjptD5OEbhmEYBVIb6Qmcc05E3FCH\npBXe0n8rt/ovBxO6zoMbbr08dWiqPhZl7pOb9mkvyzrR7XEHU3zZmXzdy5uq1Thp+W5aidPr8m7x\n0Ksexpk5cyYzZ85sy7lqw6w3V0RWc869KCKrAy/F5c8B49Vxa8Zlg5jWP5VPDZwS/3VL193Uw62X\npw5N1ceibH06ScmyYfB/Xt40irnvxaW+bA8u9/IswgartBQInTDEZTHaVQzjTJ8+nenTp/u/BwYG\nhn2u2jDrXQPsB5wW/75alV8sImcQhW/WAeWaLEYZbs5O18tTh6bqY1G2Pn1XfcF1u4fJ2HV5EIAB\n+nxZKykQOmGIzbPvDWqNDhCRS4gmY1cWkWeBE4BTgV+LyFeJl14COOdmicivgVnAQuBQ51xmiKcM\nN2en6+WpQ1P1sShbn9yRwcDL+AVe/iNTAfg5B/qyVlIgmGff2559O6k1OsA5t3fGv2ZkHH8KcEra\n/xanDDdnp+vlqUNT9bEoQ58eU7JeR//4WiHSeQRnASH9AbSWAsE8e/Psm6VWpPIy3JydrpenDk3V\nx6IMfZq08eDXCAJszO1ePpyzgZD+AFpLgWCevXn2zVIrUnkZbs5O18tTh6bqY1Fcn67wZe7AYOBX\nGPuCly9iHy9fzh5ASH8AraVAMM/ePPtmqRWpvAw3Z6fr5alDU/WxKKqe23h3X5YsqwS4b8mwquZk\njvfyVtwKhPQHkF/2Sl3ejUbbPPuRYa8lNAzDqAC1IpWXwRPrdL08dWiqPhZF1ZMN1GsEtwie6Pb8\nwcv7cqGXk1cJJrluIL9Uxbq8Wzz0Vs9hZFMrUnkZbs5O18tTh6bqY9H5epGRd3unv0bw+3zHy9ez\nnZeTfDeNXiMIFsaxMM7IMGPf4Xpm7Lujba3Wc4dERl7WTH+N4Dkc6uW0FAiNXiMI5tmbZz8yzNh3\nuJ4Z++5oWzP1ZilZ1ouM/NOTwmsED+JnXt6Vq1S9wSkQOpHjRpeX2Wibgc8HM/YdrmfGvjva1ky9\nyVPUOvodIs9ev0bwGE7z8q1s5eW0FAidSFWsy8totFvRZ7SOGfsO1zNj3x1ty66n1tF/JcTnx6z6\nGgBXspsvu5gvenl9HvJyWgoE8+zNs88bM/YdrmfGvjvallXPbaDW0U8Mnv2Dy60L1L9GcDuu93Kj\nFAhl9J7Ns+8tzNh3uJ4Z++5om5br0hOrFAj/mB48+63jVwnq1wjey0ZebpQCoczec1naZowMM/Yd\nrmfGvjvapuW69MSfDwZ+/JjHvXwWRwBwLTv7slZSIJTRey5D24z2YTtoDcMwKoB59h2uZ559d7RN\ny8kaeqhfR387G3v5bA4HYHP+7MueVS9ta7SOviyhEgvj9C4FG/ufKfkgJffX/W5VLnM9M/bd0TaA\nObGcrKEHmDN5ZS/vxwVe3pPLgJD+AFpLgWBhHAvj5E2hxt6tFQy8PNPv5W4xBmUzcJqqj0U76q0Z\nT8a6bYNnP5n7vHwCJ3o5eZVgo9cIQrm95zK0zciHQo29fCJ4TP3PhBuqW4xB2QycpupjMfx6d3o5\nWUc/Ztxrvuy3agL2l+zv5Y24FwjpD6C1FAjm2Rt5U6ixd19QsdBr+71cbmMwsnpm7MvdNrfFZl5O\n8tE/vNxEX3Ys/8/LaevoW8lFr8vNszfDnzeFGntZK3j2jxEM/6QSG4MyGzhN1ceilXqXKlk2VOvo\nt4yuyanc4ssO5Rwvp71KsJX0xLrcPHsjbwo19u9sGj70D66vQjoPRTdZWYxBtxg4TdXHopV6ey2r\n1tF/bvA6+nM4zJddxa5e1uvokxQIZc5eOVx9thO2NyjU2P/z6L95WecZkW+eD0C/iol2oxExY1+8\njux6Z3vZ/Wv60so7mALAv/MtX5a8RhDSXyVY5hw3w9VnSyh7g0KN/W/4vJdlgvKu/inODf5kvy/r\nHiNSnA5N1ceiUT236hFelklqaeWkwUsrk2WV0HhpZZmzVw5Xn3n2vYHtoDUMw6gAhXr2OkOg20p9\nlY4nyfqfHLwcsxm5qt6spupj0aiebKq+SX566HX0NzHDlyWvEYT0dfRl9tCHq8/COL1Bocb+YM71\n8thVnvay2z0O41wevhL3q1homY2IGfty6aivN8+X6WW/y4592ctXqgnYZB19soYeGq+jL7PRHq4+\n2wnbGxRq7K/mc17+A9t7OVmSqV/eLJf0e7l8RqQcOjRVH4s02e09xpfpZb8Pj5noZb2OfgeuA+pf\nI9goBUKZjfZw9dl6+d6gUGO/Fbd5+QjO8vKbm44GQGaqiTO1Dn/NkhmRsujQVH0stOxz3HwsXE96\n2e+nucnLB/MfXv4zmwOtpUAos9Eerj4z8L1Bocb+fjbw8lGc6eWPj3oQALeN8uyvUuvw7xk6ll9V\nA6ep+lho2ee4UdfThNGPePlsjvRy2jr6tNcIQvdlrxyuPjPwvUGhxn5dHvPyZezp5YvYBwAZ/w9f\nVhfSuedZL/erNLLVNHDhIajpp/wb0/LVEa6R5NrR19PtbOHlMznKy2nr6FtJgVBmoz1cfWbge4NC\njb3O9z2D//FyHwMAPD1uLV+m3/fptlOG//p+L1fT2Iex0FRzLILstgvXVnLtPD1urC87gPO8rNfR\nP6C+bQ4nBUKZjfZw9ZmB7w0KNfZj1WTXPep9ncl7PHfht77MfVIZ+NuC4Z91vV421w9Uy8D9+9/D\nShJN/zKr5KKvzGMxS8mi0m8k184U7vBlR3Oal/XSynUJu7qTVwkuzXxfVkXP3ugNbFOVYRhGBSjU\ns39FvYRZx++vZScATufbvqxuHf6ndRinuxKotV1H7MEvTif7NNx67dYxef30jVLJtfMrvuTL0nLR\nQ/o6+l5JaNaqPqO3KNTY65UNs5ng5WRJ5o8I+Uvq1uHrPDqfV4b/oRcB6Gc1X9brBm68+yJp9Es+\n+so3Fi96ue5aUNfIffE6+WQuCMIaemi8jr5XUhU3o8/oXQo19vNZ2svj4skwCMmmdudyX6bX4b88\nKUyY6SRWblq86uKWk3xZP99Vcv+QcvcYOCVLkDWd7NNw6w1fR/h83bTwYNfXgr5GduVKIMwFQVhD\nD43X0Xej0bZJV2NxCjX2euJLh3QmxDffTLb2ZXod/lS1PM5topdkRjf7Y7cMfhEKdLuBS/9/30uk\ns2o++sowFo+pB7j8i3rYb6In60OOm5M4HoBr1SsF03LRw8hTIHSLZ29Uj0KNvY6F6pvsecYBsDH3\n+DK9Dv8/OMTLK6zxgpfd1Niz/5OK49/RfDK1Mhu4urJfqLX1q2YsvYyP6f9qeZPJDVfHpE8qAz81\n9E9fC9ewi5fP4VAApqod28lrBKHxOvoyG20z8EazFGrs9UWoQzor8QoAf2NdX6bX4Z/K0V6+c0n1\nztB1452SegL3DmX4GdrwldnA9fOzUKYMeF+ITNTz1WSyWtXjoA60M59jIznl8103fL6PLDnBy9/g\nDC/vytVA41z0MPIUCGX07A0DCjb2+uLUIZ3k5huvdkHqdfj7cqGXD+B8L7/5kTinzkZZG7Bu93I/\nWyq5v+53q3JnjH0w1G+4kKxrQI4ljeSYfmnewLenne0+NnxmyWepP9/kMwfYnj94+Qj1JqqbmQ7A\nRP7Pl72qwoaN1tGXxWjbRKsxEgo19vri1DdZcvPpWKpemnkdO3j5eE728maj7gTATVIGfutgGObd\nHsrH/P0kLyeTuOUxcCn1/kd9QxHl2X+FVJKHQF29bbolpBU+m3nLhIdy8lnqz3fyqBCb/75aqqvD\nfhvwANA4PTF0Jo+MTbQaRVCosc+6kJOQjo6l6tQKm/NnL5+v1kufxwEArDxpji9zk5Xh3y0Yvisu\nDOW7l8LAZcgz4ti7MtQHuuCVDsirpOGPUQ+G5FwA/TeVN6R1hZqAHaM+s+Sz1J/vf/NZL+t3xc5Q\nmSyTcGASHoT8UiDYRKtRVmwHrWEYRgUoTRhHey2Jp6W9L51HR0+07cQ1Xj6R7wHwlw/8iy+TjVTm\nzIeVl//JwSt2yuPZp3jgLkxK/1xCXpe+MD9dh/f4VT0d/qnTkTJxXeRY7K5X26ynPrP4s3zyA2v7\nssP4sZd1QrO72dTLydxPo1z0UL4wjsa8eWMklDKMk5RnrcP/iJpou42pXk42zezLBb7snZXU24k+\npYzIQ8qIzI5DJS/kv0wx+/9PKTm04yZ3AwAzZFtf1vcbL3Ly50klOWZAPRiSc8Hihl/rXrupfmTJ\nwz52dfXZTFCfjfrM3lkpui70RKzef3EdO3pZr6NPYvWtxOl1ucXejV6glJ59Uq5vyKz4/RTu9/LV\n8YsnkhdFA2w5+k9edispI7KXMi6PxKs83lDe7vx8JjPry2YpOXirOpPltyTKfXPUwjAWA7XglfaN\nI5WB+CGg66EeGFqHzq+TtKlfpRDIzdgvrT6DsRmfjfrMpoyOslbqnEkXsq+XN+UuL+t19Emsfj7h\nwd+JdMAWezfKRKk9e12WFdLRa/GTnDo/50Bf9p9qyeKEieHtRO5JZVwOj4zLm4eFsuVryvAvbN9k\nZj+vqrLwbaX/r0rfMkHf1m4TAM6Uu31ZX9hTxk/DW/TqSI7RD4bkXACIMvBa90eTkFZGO9th7OOx\nfdMN/gwGtQ38AAAZIElEQVQA3FLqjVLqM/sN0ROsfiI27L/Q4T2dfiMJ37QSutHltsnJ6AVGZOxF\nZDbwFrAIWOCc20REVgQuAz4MzAa+4Jx7I61+KzeODuno2Ktei5+8eGJnrvVl/8YJXr6VaV4eu63K\nonlu7Nn/IBicO5Xh32xlZQxf0Yb/uvj3DqrsUiXvFeTYwNU9OP402MgCjHNf9vLNEu0p6AvbCRgI\nC5A4ailSGYgfAnX11AND66gL6cRt6t9ClWU++Ab3tV8lGKsbFzWGd8ZjuPxZysAvUFkq1WfzF8L8\ny9fi+LyOzescNxNUjpu0dfSthG50uW1yMnqBkXr2DpjunHtNlR0D3OicO11Ejo7/Piatcis3Tto6\nfKhfiz8pfvHErWzly3Tyq4PUbtIn3w1hk2W/HIU03InK0zwnGKJLleHfK21id1lV9rYyhlup8tvi\nY+9UZZupVAdqIvV5PQEbZ4wY2Dgc2hfeqMdFIUpVR3KMfjD0hewTDEjYmJY2iZvVzrQ+QRiDuv6r\nsbpUpa3YLB5b90woW/aEEFZ68d3wRqndRl/l5cM4B4DfqYeIjs3r8N5Ic9zoctvkZPQC7QjjLJ6c\nZWfwLvQFwEwyjP1wb5xGsfwp8SYaqH+BtN6A9dnRv/Py209EIY0xx4VnljtEGf6LgtE6+0uh/Iij\nYgN3pjJw31PG8N9U+d/iY9cNZXon7ApqJ+wJaun8QOyg6s1TA78McnqCYxj4U0o99cDQOk5UD5ew\n83Zw2xdvf1pf+49SY6XGZS81hu6/o/Ix/xHGe97zK3p564m/D22mz8vJnor62PxHvKyvhbT4fCfW\nwJuBN8pKOzz7m0RkEXCuc+5nwFjnXBJUnwuMzao83JtMx17rY/lRCkgdx5/OzV7WsfyzONLLW02M\nVqnMuzUYnBV+qRKsbasMvwq9nBSHOr6rX6Ci0jOc9mYwZv3LR+UvuR/4slXlm14+2C3v5RPlTS/3\nxXPNAyEaVTcpe1kITdeRHKMfDH1h3to/RBbXfW780NHt1IY/rU8A/fEYnKT6f4QaK3ecSlh2QzS2\n8+4M473V1LBS6By+5mUdn98m3ij1AFN8mY7NN0qBYMsijSozUmO/hXPuBRFZBbhRRB7V/3TOORFx\naRVv6Q9pij88/cOMn/5P/u9WPHt9Uyc3u47j30twZ3dSsfzTVDK1S9kbqDc4b1wattaPn/l46JPa\nkbtsvKLlJDWher572MtHSzBmifHUBn4PF7zScyUsJ02Lz/cpr/x0FY6ZQTqnPz+4Xt0Doy6WHx4u\nSZsuV+3Uhl/36ULV15PiB8KpapWP2ziMy/hZYQyTsd1qrzDeyWcA8B1O9/Ln4iRmALfF4Tkdm39J\nhfHameNGl5shN4pi5syZzJw5sy3nGpGxd869EP9+WUSuAjYB5orIas65F0VkdSA14/q0/ql1f9vt\nZBiGUc/06dOZPn26/3tgYCD74AYM29iLyBhglHPubRFZBtgWGACuAfYDTot/X511jnbESrUH90He\nBuq9vXXjSVsIniHAl1XmzGOJ4tQ3zQvr0Gfs9d9efvaMdby88azbvfz22CjWv5XarHSofCzUUzH5\nxKPfzm3gyy6XMLfQd70XGdguyH3x+zYGwkZhlc0/mhBJY7/kXOpbQF94d0f9xK3WHX/D0O28Xnn5\nep5B9/WqeAySMQHYeG4YKz2GM74Rja0e74PG/KeX9Wejk94ln6WeiE0+c2i8jt7WwBtVZiSe/Vjg\nKom+vteAi5xzN4jIPcCvReSrxEsvs07QKP7ZzA2p4/fJzZ61AWsDNXGrjcih/ASAQ8b81Jfd9ERI\nsLXHN37l5XuOCVkY95z7SwBu+GAwWgOu38uT1KTrP7tPAnC93OHL+kJzGAi2tW4dvV9CGZ5f/FR9\nV/o46VyZnEvV0w+MOh364RK3aUA9iKbGbQd4VPXpadXXZAz2fOeXvkyP1R6nhjFMxnb/iWG8k88A\n4GI17aw/s8eYBNR/vm9nJDRrJYyjMQNv9CrDNvbOuaeADVLKXyM7lNwSrXr2yc2uDcCqKoqUGAuo\nz5yZpMPVq3UOnvhDL1/+p328fNipIX592f5fAWDgnRD//84S/V6+20338tMyE4C+Z3wRA2sFuU/F\n0wfURGrfetHvH4W9RXWD+7+kk2j+kXowJOeC8BAZpDv+RHU7f6keUGurPu2p+nr6gmgMkjEBOOz8\nMFZ6DA/eIhpbPd7ncJiX9WfzgLrEks9ST8Q2MvBgsXfDABDnUudP81Uq4r7njmv7eUNOnXm+TH+1\n1ys3tMef7Ly92ZtIOJ5TvKwnc0+9vD/Ie3wdgGN2CQ+G3/12Gy9vIGF358p/j35fvUxo757K4P5I\nedqHT1Plt0S/9dej3ys5a6lTshzqM6rs10pO0wFweNymy1R7Pvf3IL+i2v+AC33dcZeor6f+9uu+\n7JjLw7gcs0e/l48mWup5MuEa2FoFpHS4TU+2J6+rzPp8zZgbWZyoHItuRkRwTm09b6VuLxn7tK/r\nY5RhmKcMg17RMZsJAGxHCF7fpPznvpfDOvQzVjnUy984Owo9XHlEMKm7bRlM8Yu3hyWN/4hXvEy4\n0hdx7W5B3kklNPu5SnS2X3yKC8KCGZSDzhzSWTP+rb4Q+HMtfr4Dle5rY907qXbOVu1cSi3TXG3L\ncJIrb4/GYLezQ//POEKN1cshTDOwSvTw1DnnryfEktI+Gwif5bwmctwYhsaMfcG5cdpN2td1bRh0\nrFcbkWQXpvbs91SpAM5a5V+9/I0Hg9H6zRE7AfD548OSziduX9PLEycrUxzPVT4SwtjsFJb6c+1Z\nQT4wzGXyy3jFoto0q7YUoQIa9SQLHXW9i7SBVzquVQ+XpE2PKAO/XphnhcnhJE/MCn3d7fjIyP/m\n5J18mR6rsz4RxjAZ28vZI5xW7YSdpZKw6c8sSZNhBt4wWqenPPs0luRdLy9Sz7aV1VuLXo+NiN55\nqzdm7ct/eVl7oF/6/RUA3P+Z4GtP2V/50ioW7l95qwz8fSE6xIbKKl+pUiAkWdlnqlPp0M1c0hmb\n8v/pStYPjN2U7vti3Ruepg5QDyK04VdzC/efH43BlN+H/l/0md29rL81/VecqVKvlLpfxeY/pAz8\nK6zs5VEsBOA9wntnDaMZzLPvUWOf5fktq5bpvasm88YR7ejUE39JHB/qPc0DXr7Yy7euEmWRnHpu\nSDD2bsgvxmidaz4unx828bJ0WBDEn5R3vYkKt1wZO9LBh0YFObK/mi2Mf09QZTrks5vScbfy+LeI\n2zw/5DNj6ZBeCLUqkndVm0fH5bceHDJrTn05jMt5q4QVNokXr2PzK6ksm/pdsaPVBGwy8W6evdEq\nZuzttYSGYRiVoCc9e01WGGcFXvdysqJDvwFLe/lpL68G2PHeaAXKcxuFY9c4XmUYU+vXffhGrYJ5\nToVH1lAx9JseV+Xx7/vUqZZW8luks1z8e74q21DJzyl5htL9XKx7DRVuQq3WqQvpqM1Yz50cjcEa\n94b+/26jsFpHh2ySyW/tzf+fSmimV9u8wYe8bGEcY7iYZ9+jxj7ra75embNQla8cGx0d2tEThq+q\nuPEOi0J8Y/aoCQB8/KrHgnK986CPweU/UGXKyN6nDOoa6pAkfL+cKtPhmCVIZ0H8W4d/9INBT9xq\nw79h8jBSDxy+qWS1EYyBweX/u2vYyzBh0WwvXzcqxKySN0fp8JgO1+hXUNZSJtstjGO0ihn7HjX2\nmqy11zp+nxgMvQ7/XeU96lcf6henbP1UtNnorbXCQ2K5s4LRUm8aDB6x2kj12iVBXjHkF2Om2giV\nePFPqFPpOL323DVJvYWqbKKSdb3papfta3FSyxX3VgeoDVZ6glm9upa3jozGYLlnQv9vXjvsvNWr\nau6Ps1aOVt+6kjX0UP+Z6Q1yCWbgjVYxY18BY6/RIR2NXv2RoNd6a49/U8Kk4/zYpE78o/K19VpI\nPbGZGPNfqTLlrj+iPGZtzB+Kf2sPXgWKGqKbs0DJ6ytZPxDWS76B6K8B+yhZv6ZGTTYnjXriU+G7\nhM5CeRdh4jbx4mfXTR8HXlcPVI2Fb4zhYsa+AsY+y7PXIZ0EHUMepUzgRBXL12zx8r0ALFQ2qKa8\n9TrblKxcUZb8rZA/jSWUs3qX2rGatEIvodRGWxtqTaJGPyTGpvwfYFO1K3ZBPFzLbasO0Er0CiP1\n7FwYfxOoqbI/rbJRatuSF47oOZRXM3YM6H0SCebZG61ixr7HNlWlUZ80LVgi7a0nIZ35aupz9YyX\nYqzvfW14ZpUo/rHWXSruEhJBwlVKTuLzysAvFzI28IiKkauIjg/fZIVutOFPI2syV4d0ZquHy3pJ\nO3XMXhv+h5UcXgJGLV5e/8ymISakw2L6ZeDJeL+gQjf6s9GhG/3QNc/eMIZPzxv7+sm80F0dL06M\nyJiMbIprqfwsc5V/PGVeHHsJ712B3ylZx+wTI6/eMjX/ziBrA/+QkhNjnbXqphH6waAnefVLBnRI\nZ35sn5feTBWq2Hyd4dcPhB2jX2PnhTPfPybMVus5kmfivEQ6idk8FbrRk7J67iQtHYZhGM3R88Y+\nK1vmQganSX4vxduHesOvPf5XxsTLDR9XUXQ1Aau9eG/k1aqbpdXGJu1d6ynJ2QymmTBOgv6A9QNj\ngpJ1GH5s0ib9ukO1XFQ99+oNf3z8K58I34L08lb9joEkZv9OnQevDXz4HDRm5A1j+PS8sW+GtJw6\n2uCspFIr6HXf4xbF76lVq1nqMo/p8sQLVi78W8pj1uEWvRQyibm/zfDQDwMdv9fn07H8t2LLv5z+\nVqK/BqyTUR5njFh2UTjzs6NCnEqvnX/dctwYRsexHbSGYRgVwDx7RdqkLdTnTB+r3Nm3R0VhiOVe\nVWGcrOB7EjBXnv9yahXM7H8EudE6eu2tN5qg1d68PpduZl14J2mTju1ob16v+/zY4PK3Vw+hGb3i\nSYdxkph8VrjGMIz2Y8aeEELImhjUG4L0UkD/Yg0Ve+dRJaeVKwP/lorTZxnl5ANqZNSz0PX0XIDW\noXUnbdIPorqnwXpDl2sDr18QMyplbPV4z7cwjmHkSmWNvV6ZkxidhXWTucF/frfudXehPDH8H3pX\nmU5tJJUx988ObeCVfcvaCZtW3ozhT47J+pawXEa5b1Na2xcv132NFzfVv1hk6DFcmGHgF1X3sjSM\n3KjUXdXMC6cT9JrumvLsU5dvamOowjEspeTESKoRn5++oTd1hU3WqptGm6parZe0aWl9ZWgDr/uk\n+zo6+ZWeeE578a2sl7eJW8NoD5Uy9o3I8iizPFAf6ln4TijUp0gxhnVliixvvdHSylbIOleq7rS2\nL16u+xqf/N0MQ542hubBG0bnsLuNYIiamS5M9TRTjN6g8hQvfkGG9U0zvq2src86Ri/vzHq4pLZJ\nt72Zvsa04pVnPVANw2gPZuwV2jjVy+nD5I/JMnpphjMjepS1wqZRfD7r/2mpj7POm/rwyIpyNehr\nloHXY5gcYyEaw+gcZuybIOsh4NFx7BaC5Fk5bhp57sOdoM1qjj5f6kRxM31aavC/G46bYRgdwzZV\nGYZhVADz7DPIDkekeKuthHEyvORGzvNwJ2oXZsgNdTdzcIMwTjNjaBhGZ6issU8zOIsyAtX62NSJ\nxKxR1OXJqZWxzJwkzShPo9HSyywa6tYn1sPSTF/9KYY2/FnzIvYwMIz2U1lj34iFjQw8TRilBp59\nM2kPhrtztpVzpc4XtLpAPyZrTJoZT8Mw8sOMfYtoD9QbrWZGMcVIZk3KptnTdiy9zIowNXzoNKOw\nlhzaeBWTYRidx+5GRZahainE0GhEM5Y0trQscpi0tLxz6A3GEQ36mh2mqcXtMQ/fMDqFGfsmaLR2\n3Kl/Z74cMra0C1oM43SCtIeLbucSTTxxkjFouCfBMIxCMGOPfnnJwkFlzciL1CjWDWiKd7ywiQna\nVnLjNCLr5SWNUifodtZt0Mrw+JMxaGXc6k9rDwPDyBMz9or6UEO6VWu0wiSTxLNfNKhokKzJa4JW\np05Ia4du59ItPGkarlxSx1hM3zA6h91tTZAVe17YyLOvP0lUJ8OzbyVFQjtopFu3s5n4/aIGE7Tm\nuRtGsdgOWsMwjApgnn0GWevC07z8haPCM3P0qPf1SQbJWRO0zYR02kkj3QuaaZxy1pMxyPLmbZ29\nYRSLGXtF/ctN9GRtugELE7R6GN9LP/nCul9Ac6GbvJZe6knX1DBOi41IxqCVzKEW2jGMzmHGPoNm\nPFRv7Ec1YbTiuHczG6k6HbNP9ez1wc3E7EcNTlvc6BuRYRidw+66FtGv1Eveq1rnoTbYpqpTCM+j\nfCRtqkt1nPVUUn1NxkC/a7aV1w8ahpEvZuyboFFooqkYdBKzH1w0qLzTpG6qIqVwCBamLKe0fPaG\nUR7M2GfQTDjiPe/ZZ6y9XDRYbuaFJUVO0CZtqvPsU/oBLObZR3+8pzz7ZsbQMIzOYMa+CbQxf7cu\njBPJ7zXz9tr4Rd3aiGYZ/k6T9qaqOmOf8ZJ0TTIGaeMDFqc3jKKxO7BF0jx7Haeuc1q1yxy/tHt+\nxr/LQtKmOmP/bsoBUNfXZAyyPHvDMIrFjH0TZMWeg2evJiIzJmjnp3j2eWW3bJVGYZz5yrNfOmOC\nNhmDem/eQjeGURZsB61hGEYFMM++RbTnmoQs5jEmHJDh2b/19+h3Vsy+LKTF7JO2Q7Znn4yBDuO8\na0svDaM05GLsRWR74IdEUd2fO+dOy0NPEehwRGLg5uscktq+qVj3W4v9hnLH7HU7tTxWx+9VX5Mx\n0A8+C90YRnloexhHREYBPwa2ByYDe4vIeu3WUwbmsTTzWJrXWcH/sBTMfBBYisjYxz8vEf28rX7K\njG7nS+pn5hxCv5YKP0n/kzGZV5dEuXuYPfPpopuQK9a/6pJHzH4T4Ann3Gzn3ALgUmCXHPQUznuM\n5j1G8yor+x+Wh5kPAMsDb4afOUQ/89VPmdHtnKN+Zs4h9Gv58JP0PxmTbt09+3SPGwvrX3XJI4yz\nBvCs+nsOsGkOekrDs4z38uurLM38ZRbw+ipL8KGXgkl/tIiGtYm6tv+dyMUn6muCHgPDMMpHHsbe\nNXvgKDnFy/30D5LTypqRi6y3Iv3AzZw+sHXH26bJS9+0l/oZiI19vxw9In1l/HyncTujBm4ftr4y\n9imtf2VsWzvqCTBq4HYWueMw6hHnmrbNzZ1QZDOg3zm3ffz3scD7epJWRNqr1DAMoyI452Q49fIw\n9jXgb8A2wPPA3cDezrlH2qrIMAzDaJpau0/onFsoIl8DridaevkLM/SGYRjF0nbP3jAMwygfHU2X\nICLbi8ijIvK4iJrd61JEZLyI3CwifxWRh0XkiLh8RRG5UUQeE5EbRGSFots6EkRklIjcLyLXxn/3\nTP9EZAURuVxEHhGRWSKyaa/0T0SOja/Nh0TkYhEZ3c19E5HzRGSuiDykyjL7E/f/8djmbFtMq5sn\no3/fj6/NB0XkShFZXv2vpf51zNj36GarBcBRzrmPApsBh8V9Oga40Tk3Cfif+O9u5khgFmGlVS/1\n7yzgOufcesDHiVaadn3/RGQCcBCwoXNufaKQ6l50d9/OJ7IfmtT+iMhkYE8iW7M98BMRKXsusLT+\n3QB81Dn3CeAx4FgYXv862fme22zlnHvROfdALL8DPEK0z2Bn4IL4sAuAzxXTwpEjImsCOwA/B5JV\nAD3Rv9hL2so5dx5E803OuTfpjf69ReSMjIkXTYwhWjDRtX1zzt0GvL5YcVZ/dgEucc4tcM7NBp4g\nskGlJa1/zrkbnXPvx3/eBawZyy33r5PGPm2z1Rod1J8rsSc1hegDGeucmxv/ay4wtqBmtYMzgW8D\n76uyXunf2sDLInK+iNwnIj8TkWXogf45514DfgA8Q2Tk33DO3UgP9G0xsvozjsjGJPSCvTkAuC6W\nW+5fJ419z84Ei8gHgSuAI51zdWlvXDQD3pV9F5HPAi855+4nePV1dHP/iFajbQj8xDm3IdH+4Lqw\nRrf2T0Q+AnwdmEBkGD4oIvvoY7q1b1k00Z+u7auIHA+855y7eIjDhuxfJ439c1C3p3489U+mrkRE\nliAy9Bc6566Oi+eKyGrx/1fHJxjoOjYHdhaRp4BLgE+JyIX0Tv/mAHOcc3+J/76cyPi/2AP92xj4\ns3PuVefcQuBK4JP0Rt80Wdfi4vZmzbis6xCRrxCFUr+kilvuXyeN/T3AOiIyQUSWJJpcuKaD+tuO\niAjwC2CWc+6H6l/XAPvF8n7A1YvX7Qacc8c558Y759Ymmtz7o3Puy/RO/14EnhWRSXHRDOCvwLV0\nf/8eBTYTkaXj63QG0SR7L/RNk3UtXgPsJSJLisjawDpEGzy7ijhd/LeBXZxz+m3QrffPOdexH+Az\nRLtrnwCO7aTunPqzJVEs+wHg/vhne2BF4Cai2fMbgBWKbmsb+joNuCaWe6Z/wCeAvwAPEnm/y/dK\n/4DvED28HiKavFyim/tG9O3yeeA9ovm//YfqD3BcbGseBbYruv3D6N8BwOPA08q+/GS4/bNNVYZh\nGBWg7OtODcMwjDZgxt4wDKMCmLE3DMOoAGbsDcMwKoAZe8MwjApgxt4wDKMCmLE3DMOoAGbsDcMw\nKsD/ByDUfz1s69mQAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7fd398e199b0>"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.sparse import dia_matrix, eye, kron\n",
      "from scipy.sparse.linalg import factorized\n",
      "def d2matrix(nelem):\n",
      "    elements = ones((3,nelem))\n",
      "    elements[1,:] *= -2\n",
      "    return dia_matrix((elements, [-1,0,1]), shape=(nelem,nelem)).tocsc()\n",
      "\n",
      "def diffusion_solve_implicit(U, dx, dt, nx, nt, D):\n",
      "    alpha = D*dt/dx**2\n",
      "    mat2D = kron(d2matrix(nx), eye(nx)) + kron(eye(nx), d2matrix(nx))\n",
      "    M = eye(nx*nx)-mat2D*alpha\n",
      "    LU = factorized(M.tocsc())\n",
      "    for it in range(0,nt-1):\n",
      "        U[it+1,1:-1, 1:-1] = LU(U[it,1:-1, 1:-1].ravel()).reshape(nx, nx)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "U, dx, dt, x, t = diffusion_init(maxx, maxt, D, nx, SC=10)\n",
      "X, Y = meshgrid(x, x)\n",
      "U[0, :, :] = 0.\n",
      "U[0, (X<maxx*0.4) & (X>maxx*0.2) & (Y<maxx*0.4) & (Y>maxx*0.2) ] = 1."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "diffusion_solve_implicit(U, dx, dt, nx, len(t), D)\n",
      "pcolormesh(U[:, 30, :], rasterized=True, vmin=-1, vmax=1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Again, Crank-Nicolson follows:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def diffusion_solve_CN(U, dx, dt, nx, nt, D):\n",
      "    alpha2 = D*dt/dx**2*0.5\n",
      "    mat2D = kron(d2matrix(nx), eye(nx)) + kron(eye(nx), d2matrix(nx))\n",
      "    M1 = eye(nx*nx)-mat2D*alpha2\n",
      "    M2 = eye(nx*nx)+mat2D*alpha2\n",
      "    LU = factorized(M1.tocsc())\n",
      "    for it in range(0,nt-1):\n",
      "        U[it+1,1:-1, 1:-1] = LU(M2.dot(U[it,1:-1, 1:-1].ravel())).reshape(nx, nx)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "diffusion_solve_CN(U, dx, dt, nx, len(t), D)\n",
      "pcolormesh(U[:, 30, :], rasterized=True, vmin=-1, vmax=1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Neumann BC"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def diffusion_solve_CN_neumann_zero(U, dx, dt, nx, nt, D):\n",
      "    alpha = D*dt/dx**2\n",
      "    d2x = d2matrix(nx+2)\n",
      "    d2x[0, 1] += 1\n",
      "    d2x[-1, -2] += 1\n",
      "    mat2D = kron(d2x, eye(nx+2)) + kron(eye(nx+2), d2x)\n",
      "    M1 = eye((nx+2)**2)-mat2D*alpha/2\n",
      "    M2 = eye((nx+2)**2)+mat2D*alpha/2\n",
      "    LU = factorized(M1.tocsc())\n",
      "    for it in range(0,nt-1):\n",
      "        # U[it, boundary] -= alpha*dx*(NBC[it]+NBC[it+1])\n",
      "        U[it+1,:, :] = LU(M2.dot(U[it,:, :].ravel())).reshape(nx+2, nx+2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "U, dx, dt, x, t = diffusion_init(maxx, maxt*50, D, nx, SC=50)\n",
      "X, Y = meshgrid(x, x)\n",
      "U[0, :, :] = 0.\n",
      "U[0, (X<maxx*0.4) & (X>maxx*0.2) & (Y<maxx*0.4) & (Y>maxx*0.2) ] = 1."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "diffusion_solve_CN_neumann_zero(U, dx, dt, nx, len(t), D)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pcolormesh(U[:, 30, :], rasterized=True, vmin=0, vmax=.1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# constants taken from http://en.wikipedia.org/wiki/File:Brusselator_space.gif\n",
      "#initial field values\n",
      "U0 = 1.\n",
      "V0 = 1.\n",
      "noise = 2.\n",
      "\n",
      "#diffusion coefficients\n",
      "DU = 0.2\n",
      "DV = 0.02\n",
      "\n",
      "#constant concentrations\n",
      "A = 1\n",
      "B = 3\n",
      "\n",
      "def RHS(U, V):\n",
      "    fV = B*U - U**2*V\n",
      "    fU = A - fV - U\n",
      "    return fU, fV\n",
      "\n",
      "maxt = 100.\n",
      "maxx = 100.\n",
      "\n",
      "U, dx, dt, x, t = diffusion_init(maxx, maxt, max((DU, DV)), nx, SC=.05)\n",
      "V, dx, dt, x, t = diffusion_init(maxx, maxt, max((DU, DV)), nx, SC=.05)\n",
      "X, Y = meshgrid(x, x)\n",
      "U[0, :, :] = U0 + noise*randn(nx+2, nx+2)\n",
      "V[0, :, :] = V0 + noise*randn(nx+2, nx+2)\n",
      "#for i in range(10):\n",
      "#    U[0, randint(0, nx), randint(0, nx)] = 10\n",
      "#    V[0, randint(0, nx), randint(0, nx)] = 10"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def brusselator(U, V, dx, dt, nx, nt, DU, DV):\n",
      "    alphaU = DU*dt/dx**2\n",
      "    alphaV = DV*dt/dx**2\n",
      "    d2x = d2matrix(nx+2)\n",
      "    d2x[0, 1] += 1\n",
      "    d2x[-1, -2] += 1\n",
      "    mat2D = kron(d2x, eye(nx+2)) + kron(eye(nx+2), d2x)\n",
      "    M1U = eye((nx+2)**2)-mat2D*alphaU/2\n",
      "    M2U = eye((nx+2)**2)+mat2D*alphaU/2\n",
      "    M1V = eye((nx+2)**2)-mat2D*alphaV/2\n",
      "    M2V = eye((nx+2)**2)+mat2D*alphaV/2\n",
      "    LUU = factorized(M1U.tocsc())\n",
      "    LUV = factorized(M1V.tocsc())\n",
      "    for it in range(10):\n",
      "        U[0] = M2V.dot(U[it,:, :].ravel()).reshape(nx+2, nx+2)\n",
      "    for it in range(0,nt-1):\n",
      "        # U[it, boundary] -= alpha*dx*(NBC[it]+NBC[it+1])\n",
      "        # diffusion operator\n",
      "        U[it+1,:, :] = LUU(M2U.dot(U[it,:, :].ravel())).reshape(nx+2, nx+2)\n",
      "        V[it+1,:, :] = LUV(M2V.dot(V[it,:, :].ravel())).reshape(nx+2, nx+2)\n",
      "        \n",
      "        fU, fV = RHS(U[it+1], V[it+1])\n",
      "        U[it+1] += fU*dt\n",
      "        V[it+1] += fV*dt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "brusselator(U, V, dx, dt, nx, len(t), DU, DV)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "figure(figsize=(12,12))\n",
      "subplot(221)\n",
      "pcolormesh(U[:, 30, :], rasterized=True, vmin=0, vmax=5)\n",
      "subplot(223)\n",
      "pcolormesh(V[:, 30, :], rasterized=True, vmin=0, vmax=5)\n",
      "subplot(222)\n",
      "pcolormesh(U[1400, :-2, :-2], rasterized=True, vmin=0, vmax=5)\n",
      "subplot(224)\n",
      "pcolormesh(V[1400, :-2, :-2], rasterized=True, vmin=0, vmax=5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}