dedalus_timevarying_1d_adv_diff.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solve 1D advection diffusion problem\n",
    "\n",
    "$$\n",
    "\\frac{\\partial c}{\\partial t} + w_s \\frac{\\partial c}{\\partial z} - \\frac{\\partial}{\\partial z}\\left( F\\right) =0\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "F=K_v \\frac{\\partial c}{\\partial z}\n",
    "$$\n",
    "\n",
    "with boundary conditions\n",
    "\n",
    "$$\n",
    "F = f(t),\\ z =0\n",
    "$$\n",
    "\n",
    "$$\n",
    "F = 0,\\ z=H\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gotta find my FFTW_PATH\n",
    "!export FFTW_PATH=/path/to/your/fftw_prefix "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[sudo] password for suntans: \n"
     ]
    }
   ],
   "source": [
    "!sudo pip install dedalus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "import dedalus.public as de\n",
    "from dedalus.tools  import post\n",
    "\n",
    "from dedalus.extras.plot_tools import quad_mesh, pad_limits\n",
    "\n",
    "\n",
    "import logging\n",
    "import matplotlib.pyplot as plt\n",
    "logger = logging.getLogger(__name__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Constants\n",
    "\n",
    "H = 5. # Depth\n",
    "\n",
    "nz = 64\n",
    "\n",
    "K_0 = 1e-1 # Diffusivity - had to boost this to actually get the sediment into the water column\n",
    "ws = 1e-3 # settling velocity\n",
    "\n",
    "E = -1e-3 # constant erosion rate [c/s]\n",
    "E = -1e-2 # constant erosion rate [c/s]\n",
    "\n",
    "uniform_diffusion = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CReate the domain\n",
    "\n",
    "z = de.Chebyshev('z',nz,interval=[0.,H])\n",
    "domain = de.Domain([z],grid_dtype='float')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the problem\n",
    "problem = de.IVP(domain, variables=['c','F'])\n",
    "problem.parameters['ws'] = ws\n",
    "problem.parameters['E'] = E\n",
    "\n",
    "if uniform_diffusion:\n",
    "    problem.parameters['K'] = K_0\n",
    "else:\n",
    "    # Add vertically varying diffusion\n",
    "    ncc = domain.new_field(name = 'K')\n",
    "    ncc['g'] = (H-z)*K_0\n",
    "    problem.parameters['K'] = ncc\n",
    "\n",
    "problem.add_equation('dt(c) + ws * dz(c) - dz(F) = 0')\n",
    "problem.add_equation('F - K * dz(c) = 0')\n",
    "\n",
    "# Add the time varying boundary conditions\n",
    "problem.add_bc('left(F) = E*abs(sin(t))')\n",
    "problem.add_bc('right(F) = 0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2019-11-06 12:14:21,146 pencil 0/1 INFO :: Building pencil matrix 1/1 (~100%) Elapsed: 0s, Remaining: 0s, Rate: 4.8e+01/s\n",
      "2019-11-06 12:14:21,269 __main__ 0/1 INFO :: Time step 20\n",
      "2019-11-06 12:14:21,378 __main__ 0/1 INFO :: Time step 40\n",
      "2019-11-06 12:14:21,461 __main__ 0/1 INFO :: Time step 60\n",
      "2019-11-06 12:14:21,557 __main__ 0/1 INFO :: Time step 80\n",
      "2019-11-06 12:14:21,647 __main__ 0/1 INFO :: Time step 100\n",
      "2019-11-06 12:14:21,750 __main__ 0/1 INFO :: Time step 120\n",
      "2019-11-06 12:14:21,846 __main__ 0/1 INFO :: Time step 140\n",
      "2019-11-06 12:14:21,960 __main__ 0/1 INFO :: Time step 160\n",
      "2019-11-06 12:14:22,082 __main__ 0/1 INFO :: Time step 180\n",
      "2019-11-06 12:14:22,234 __main__ 0/1 INFO :: Time step 200\n",
      "2019-11-06 12:14:22,397 __main__ 0/1 INFO :: Time step 220\n",
      "2019-11-06 12:14:22,528 __main__ 0/1 INFO :: Time step 240\n",
      "2019-11-06 12:14:22,659 __main__ 0/1 INFO :: Time step 260\n",
      "2019-11-06 12:14:22,799 __main__ 0/1 INFO :: Time step 280\n",
      "2019-11-06 12:14:22,944 __main__ 0/1 INFO :: Time step 300\n",
      "2019-11-06 12:14:23,071 __main__ 0/1 INFO :: Time step 320\n",
      "2019-11-06 12:14:23,194 __main__ 0/1 INFO :: Time step 340\n",
      "2019-11-06 12:14:23,289 __main__ 0/1 INFO :: Time step 360\n",
      "2019-11-06 12:14:23,392 __main__ 0/1 INFO :: Time step 380\n",
      "2019-11-06 12:14:23,482 __main__ 0/1 INFO :: Time step 400\n",
      "2019-11-06 12:14:23,566 __main__ 0/1 INFO :: Time step 420\n",
      "2019-11-06 12:14:23,682 __main__ 0/1 INFO :: Time step 440\n",
      "2019-11-06 12:14:23,833 __main__ 0/1 INFO :: Time step 460\n",
      "2019-11-06 12:14:23,968 __main__ 0/1 INFO :: Time step 480\n",
      "2019-11-06 12:14:24,115 __main__ 0/1 INFO :: Time step 500\n",
      "2019-11-06 12:14:24,223 __main__ 0/1 INFO :: Time step 520\n",
      "2019-11-06 12:14:24,346 __main__ 0/1 INFO :: Time step 540\n",
      "2019-11-06 12:14:24,496 __main__ 0/1 INFO :: Time step 560\n",
      "2019-11-06 12:14:24,637 __main__ 0/1 INFO :: Time step 580\n",
      "2019-11-06 12:14:24,767 __main__ 0/1 INFO :: Time step 600\n",
      "2019-11-06 12:14:24,896 __main__ 0/1 INFO :: Time step 620\n",
      "2019-11-06 12:14:24,999 __main__ 0/1 INFO :: Time step 640\n",
      "2019-11-06 12:14:25,080 __main__ 0/1 INFO :: Time step 660\n",
      "2019-11-06 12:14:25,186 __main__ 0/1 INFO :: Time step 680\n",
      "2019-11-06 12:14:25,274 __main__ 0/1 INFO :: Time step 700\n",
      "2019-11-06 12:14:25,375 __main__ 0/1 INFO :: Time step 720\n",
      "2019-11-06 12:14:25,514 __main__ 0/1 INFO :: Time step 740\n",
      "2019-11-06 12:14:25,667 __main__ 0/1 INFO :: Time step 760\n",
      "2019-11-06 12:14:25,806 __main__ 0/1 INFO :: Time step 780\n",
      "2019-11-06 12:14:25,920 __main__ 0/1 INFO :: Time step 800\n",
      "2019-11-06 12:14:26,040 __main__ 0/1 INFO :: Time step 820\n",
      "2019-11-06 12:14:26,116 __main__ 0/1 INFO :: Time step 840\n",
      "2019-11-06 12:14:26,231 __main__ 0/1 INFO :: Time step 860\n",
      "2019-11-06 12:14:26,342 __main__ 0/1 INFO :: Time step 880\n",
      "2019-11-06 12:14:26,436 __main__ 0/1 INFO :: Time step 900\n",
      "2019-11-06 12:14:26,541 __main__ 0/1 INFO :: Time step 920\n",
      "2019-11-06 12:14:26,634 __main__ 0/1 INFO :: Time step 940\n",
      "2019-11-06 12:14:26,740 __main__ 0/1 INFO :: Time step 960\n",
      "2019-11-06 12:14:26,826 __main__ 0/1 INFO :: Time step 980\n",
      "2019-11-06 12:14:26,922 __main__ 0/1 INFO :: Time step 1000\n",
      "2019-11-06 12:14:27,010 __main__ 0/1 INFO :: Time step 1020\n",
      "2019-11-06 12:14:27,091 __main__ 0/1 INFO :: Time step 1040\n",
      "2019-11-06 12:14:27,205 __main__ 0/1 INFO :: Time step 1060\n",
      "2019-11-06 12:14:27,307 __main__ 0/1 INFO :: Time step 1080\n",
      "2019-11-06 12:14:27,397 __main__ 0/1 INFO :: Time step 1100\n",
      "2019-11-06 12:14:27,480 __main__ 0/1 INFO :: Time step 1120\n",
      "2019-11-06 12:14:27,567 __main__ 0/1 INFO :: Time step 1140\n",
      "2019-11-06 12:14:27,680 __main__ 0/1 INFO :: Time step 1160\n",
      "2019-11-06 12:14:27,770 __main__ 0/1 INFO :: Time step 1180\n",
      "2019-11-06 12:14:27,871 __main__ 0/1 INFO :: Time step 1200\n",
      "2019-11-06 12:14:27,961 __main__ 0/1 INFO :: Time step 1220\n",
      "2019-11-06 12:14:28,046 __main__ 0/1 INFO :: Time step 1240\n",
      "2019-11-06 12:14:28,148 __main__ 0/1 INFO :: Time step 1260\n",
      "2019-11-06 12:14:28,240 __main__ 0/1 INFO :: Time step 1280\n",
      "2019-11-06 12:14:28,328 __main__ 0/1 INFO :: Time step 1300\n",
      "2019-11-06 12:14:28,432 __main__ 0/1 INFO :: Time step 1320\n",
      "2019-11-06 12:14:28,520 __main__ 0/1 INFO :: Time step 1340\n",
      "2019-11-06 12:14:28,621 __main__ 0/1 INFO :: Time step 1360\n",
      "2019-11-06 12:14:28,711 __main__ 0/1 INFO :: Time step 1380\n",
      "2019-11-06 12:14:28,802 __main__ 0/1 INFO :: Time step 1400\n",
      "2019-11-06 12:14:28,895 __main__ 0/1 INFO :: Time step 1420\n",
      "2019-11-06 12:14:28,988 __main__ 0/1 INFO :: Time step 1440\n",
      "2019-11-06 12:14:29,078 __main__ 0/1 INFO :: Time step 1460\n",
      "2019-11-06 12:14:29,191 __main__ 0/1 INFO :: Time step 1480\n",
      "2019-11-06 12:14:29,342 __main__ 0/1 INFO :: Time step 1500\n",
      "2019-11-06 12:14:29,495 __main__ 0/1 INFO :: Time step 1520\n",
      "2019-11-06 12:14:29,637 __main__ 0/1 INFO :: Time step 1540\n",
      "2019-11-06 12:14:29,768 __main__ 0/1 INFO :: Time step 1560\n",
      "2019-11-06 12:14:29,900 __main__ 0/1 INFO :: Time step 1580\n",
      "2019-11-06 12:14:29,998 __main__ 0/1 INFO :: Time step 1600\n",
      "2019-11-06 12:14:30,096 __main__ 0/1 INFO :: Time step 1620\n",
      "2019-11-06 12:14:30,216 __main__ 0/1 INFO :: Time step 1640\n",
      "2019-11-06 12:14:30,379 __main__ 0/1 INFO :: Time step 1660\n",
      "2019-11-06 12:14:30,534 __main__ 0/1 INFO :: Time step 1680\n",
      "2019-11-06 12:14:30,658 __main__ 0/1 INFO :: Time step 1700\n",
      "2019-11-06 12:14:30,812 __main__ 0/1 INFO :: Time step 1720\n",
      "2019-11-06 12:14:30,936 __main__ 0/1 INFO :: Time step 1740\n",
      "2019-11-06 12:14:31,078 __main__ 0/1 INFO :: Time step 1760\n",
      "2019-11-06 12:14:31,182 __main__ 0/1 INFO :: Time step 1780\n",
      "2019-11-06 12:14:31,277 __main__ 0/1 INFO :: Time step 1800\n",
      "2019-11-06 12:14:31,386 __main__ 0/1 INFO :: Time step 1820\n",
      "2019-11-06 12:14:31,479 __main__ 0/1 INFO :: Time step 1840\n",
      "2019-11-06 12:14:31,581 __main__ 0/1 INFO :: Time step 1860\n",
      "2019-11-06 12:14:31,684 __main__ 0/1 INFO :: Time step 1880\n",
      "2019-11-06 12:14:31,716 solvers 0/1 INFO :: Simulation stop time reached.\n",
      "2019-11-06 12:14:31,723 __main__ 0/1 INFO :: Total time: 10.572451\n",
      "2019-11-06 12:14:31,724 __main__ 0/1 INFO :: Iterations: 1885\n",
      "2019-11-06 12:14:31,728 __main__ 0/1 INFO :: Average timestep: 0.010000\n"
     ]
    }
   ],
   "source": [
    "solver = problem.build_solver(de.timesteppers.RK443)\n",
    "\n",
    "solver.stop_sim_time = 6*np.pi\n",
    "solver.stop_wall_time = 1000.\n",
    "solver.stop_iteration = np.inf\n",
    "dt = 0.01\n",
    "\n",
    "# Store data for final plot\n",
    "# Initial conditions\n",
    "z = domain.grid(0)\n",
    "c = solver.state['c']\n",
    "\n",
    "c.set_scales(1, keep_data=True)\n",
    "c_list = [np.copy(c['g'])]\n",
    "t_list = [solver.sim_time]\n",
    "\n",
    "start_time = time.time()\n",
    "while solver.ok:\n",
    "    solver.step(dt)\n",
    "    if solver.iteration % 20 == 0:\n",
    "        logger.info(\"Time step {}\".format(solver.iteration))\n",
    "        c.set_scales(1, keep_data=True)\n",
    "        c_list.append(np.copy(c['g']))\n",
    "        t_list.append(solver.sim_time)\n",
    "end_time = time.time()\n",
    "\n",
    "# Print statistics\n",
    "logger.info('Total time: %f' %(end_time-start_time))\n",
    "logger.info('Iterations: %i' %solver.iteration)\n",
    "logger.info('Average timestep: %f' %(solver.sim_time/solver.iteration))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_array = np.array(c_list)\n",
    "t_array = np.array(t_list)\n",
    "zmesh, tmesh = quad_mesh(x=z, y=t_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEMCAYAAAAyO4a7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2dfZAlV3mfn/fe+d4PwSJUhhKWY6QCQ2RBkGPLMsaK5IC9BhtI7IqDKWI5qsgYIztYUDZFIUeYRZFjCxxZbPhywIUprKRIUBxFFtgQeUESYFDicmSIkGxivFoJmNXu7Hy++eN2T789031v90xPz+2d31M1tX1Pn+5zbt/Zd879vR/H3B0hhBDdpLfbExBCCLF1ZMSFEKLDyIgLIUSHkREXQogOIyMuhBAdRkZcCCE6zERbA5nZZ4AzyctVd7+yrbGFEOJspTUjDvx3d39ri+MJIcRZj7WV7GNmtwP3ArPAfe5+RysDCyHEWUybK/F3uPu9ZtYHPmVmJ939U7GDmV0DXDN4MfECmzmnxekJIWphttszAMBPnzjh7k/dzj2eYbN+hrVKfU+wdKe7v2Q74zVJayvx3KBmR4AFd7+hrE9v7lyfeNbLWpyVEKIO1uvv9hQAWPr8ez7n7pdu5x5PtWl/JU+r1PfdPLzt8ZqklegUM3u2mV0dmi4CvtzG2EIIMQoD+lbtZ9xoS06ZB37MzJ4OHAT+GvhwS2MLIcRQDJjqVbTQqzs6ldq0YsTd/f8BL29jLCHEzjEuEkrTDFbiY7jMrkCbjk0hhBhPxlQqqYKMuBBiz6OVuBBCdJjUsdnY/cyuAl4BHAe8KBLPzH4SeDvwenf/eGh/FfB8Bur7V9z93cPGkhEXQgissZW4mc0BtwHPdfdFM7vdzK5097tDn78HPMogyCNeez7wBuD57u5mdp+ZfcLd/6psPBXAEkJ0Duv1138auR8waVbppwKXAQ+7+2Ly+h7gcOzg7g+5+ycLrn0x8DnPEniOAT8ybDCtxIUQex5r1rF5HnAyvJ5P2nbkWhlxIcRQztawwo3UkFPONbP7w+uj7n40vD4OHAivDyZtVTgOXLjh2qGJkTLiQog9T03H5okRaffHgAvMbDqRVC4HbjWzQ8CKu88PufZO4HVmZomkchnwrmGTkREXQux5mgwxdPfTZnYt8E4zexT4krvfbWY3AY8DR8zMgF8DLgB+ysyW3f1Od/8bM7sZ+C0zWwXeM8ypCTLiQogSxk1G2cn5mNVIu6+Au98F3LWh7fpw7MCNyc/Gaz8EfKjqWDLiQgiBMjaFEGcBe2n1nRsHGXEhhOgs1mCyT9vIiAshBFqJCyE6yrhIKLs5j0GyTzetuIy4EGLPU2tTiDFDRlwIseeRY1MIMfaMi2wSGac5SU4RQoiOYgY9GXEhhOgqhnVUT5ERF0LsecygPzU+0k4dZMSFOMsZB915HOYwFEMrcSGE6Cxm9GTEhRCiu1ivm7tVyogLcZYwLpJFG/NoegwztBIXQoguI01cCCE6ipkpOkUIITqLgal2ihBip9hLeneb44QR6fXl2BRCiG6iOHEhhOguJiMuhBDdRnKKEKIRznb9e1zeX8TM6E/KiAshRDcxMK3EhRCiuyhjUwghuoqpnnglzGwW+CzwP9z9DW2OLcQ4MC568E7Oo8l7txaXLjmlMjcCX2h5TCGEGI4hx+YozOxngHuA7wb2tzWuEEKMwjqcsdnKrM3sOcB3uft/GtHvGjO738zu95UzbUxNCCHWMzar/Iwbba3EXw6cMbM3AT8ATJnZde7+27GTux8FjgL05s71luYmxI6wm/p3V2K8x8VHoBDDEbj729JjM5sB9m804EIIsXuYdvapgpm9EvhBBivxf+buH25zfCGEKGKws4+M+Ejc/Xbg9jbHFEKIkZjRm+pm2kw3Zy3ELnM26t1N3buRe/TbryfepJxiZlcBrwCOA+7uN2w4PwPcDHwNuAg44u4PJuf+HbDMIPBkDnidu6+VjSUjLoQQ1twfDjObA24Dnuvui2Z2u5ld6e53h27XAY+4+01mdjHwXuCFZva9wJXufklyry8ClzEIzy6kmyKQEEI0iGFYv1fppwKXAQ+7+2Ly+h7g8IY+h4FjAO7+AHCJmR0EHgP2m9mEmU0ADjw0bDCtxIUQwqBXXU4518zuD6+PJuHRKecBJ8Pr+aSNUX3c/ctmdhT4KLAG/DHw6LDJyIgLIQS14sRPuPulQ84fBw6E1weTtpF9zOxlwBXu/hIAM7sd+JfArWWDyYgLUUDbjstxd1Zu915b1Zt7rRXAMnqTjZnDY8AFZjadSCqXA7ea2SFgxd3ngTsYyC6fTjTxL7r7vJk9A/h6uNffAjPDBpMRF0IIs8YyNt39tJldC7zTzB4FvuTud5vZTcDjwBHgFuBmM3szcCFwdXL5B4DLzOxtwCpwDvDuYePJiAshRMNp9+5+F3DXhrbrw/EC8NqC604Br6ozloy4EEKgjE0hOsNuJOqMc0GqrujcO+o3MNVOEUKI7qK0eyGE6DZaiQshREcxM3qt12tpBhlxIYRAm0IIMTZ0PVHnbHdW1plba5+ldvYRQoguo+gUIYToLNZTdIoQQnQarcSFEKKrmO3qbk3bQUZcdJa2/tONw5Zl27m+jrOyjWzKKn13ZVUsIy6EEF3FQHKKEEJ0lAb32GwbGXEhhDCDiandnsWWkBEXQux5THHiQjTPOJdvbeJ+bWVV7mY2ZR3DOOp+O+rINuTYFEKI7mIy4kII0WUkpwghRFexnhybQgjRWRRiKEQ9dnS/xN0s5dqCs7JuVmWTDsMqkkMT4428R+MGV8k+QgjRXRSdIoQQXUYFsIQQottIThFCiI5iPUzRKUJkjGO25bg5K5vOpOyKg7KKU7L17E1DK/FhmFkP+K/AZ4Ep4JnAz7r7QhvjCyHEMAxTiGEFjrn7jQBm9jHgFcDvtzi+EEIUo+iU4bj7GpAa8AngfOD/tDG2EEKMRrVTKmFmLwZ+Cfi4u99fcP4a4BoAJve1OTUhxF7GDJuY3O1ZbIlWjbi73wncaWb/0cx+3t1v3XD+KHAUoDd3rrc5N7E1Ol/WtY4zsIJm2razcicdlE1kTRY6Niu876LnuONx3NZNx2Yrszaz55jZ4dD0EPCdbYwthBCjsYERr/IzZrS1El8Erjaz5wOTwHcBv9jS2EIIMRJv0ECb2VUMgjeOA+7uN2w4PwPcDHwNuAg44u4PJue+D/hhYA24AvgX7v7XZWO15dj8CoM3JIQQ44fR2CrbzOaA24Dnuvuimd1uZle6+92h23XAI+5+k5ldDLwXeKGZHQR+xd1fmdzrw8Djw8Ybv+8GQgjROjbYLLnKz2guAx5298Xk9T3A4Q19DgPHANz9AeCSxID/KPCEmf2ymb0F+AfufmrYYMrYFELseRzwfmPm8DzgZHg9n7RV6XMB8L3AzwGrwCfN7DF3/2TZYDLiYihNRgS0lTLfdsRJG1Em2xljq1EmW41q2WpkSZ1IlubT7q2OnHKumcUQ6aNJZF3KceBAeH0waaNCn3ngC+6+PJiWHQNeBMiICyHEUKob8RPufumQ88eAC8xsOpFULgduNbNDwIq7zwN3MJBdPp1o4l9093kz+yTw6nCvCxiULClFRlwIIai1Eh+Ku582s2uBd5rZo8CX3P1uM7uJgZPyCHALcLOZvRm4ELg6ufYvzeyDSd9l4G+BPxg2noy4EELQbIihu98F3LWh7fpwvAC8tuTaW4vay5ARF0IIGMtEnirIiO9hdjNlvu7YTTorq6TGN1NLe7NRaNrBl+tTw0FZxylZ5bo67bX6btHpWhtTASwhhOg0TcopbSIjLoQQmHb2EUKIztJg2n3byIgLIUSDIYZtIyO+R2jCEbRjWZPbcNql7KSzcqvZlLWceiPeX12HYp2syZ1ySkLx+6pyj6I+vV6luiVbxnvdNIfdnLUQQjRJvbT7sUJGXAghoGqFwrFDRlwIITqsiVeetZl93Myes6Htd83ss81PSwgh2sWtV+ln3Kgzo+8D7jGzf502uPu1wCONz0oIIdpmD+yx+QCDgi0fMrOXAa9294cZ1FMXY0Iru743EQmyg3W8t1q7u4kojVyfLUZmREZFmexmZEmd+5VFllhotwJNusp1TeAYa3RTE6/zZ8Xd/S+Af8hgu6EvmNnVOzMtIYRoE2fNq/2MG3VW4k82s+939z8DftXMPg78HvCMnZmaEEK0x/iZ52rUWYn/NvDt6YvEmD8PuLHpSQkhRJs4sObVfsaNyitxd/+9grZTyIgLIc4CfAylkiooTrxjjIPjcifT5NtwVraVOt5G6nvbTskyh2J0ShY5I0uvC+29AsemVRi7CdKVeBeRERdCCIdVGXEhhOguklOEEKKjOLC225PYIjLiQggBdHQhLiM+ruyUA3Mnsy3rOCub2Gh43B2UW3WqNuGUrOVcrOGULLtHdEoW9ck5MEudnJvn1FbGJsixKYQQncUdVju6FJcRF0IIJKcIIURnGcSJd9OKy4gLIQTdrZ0iIy6EEMixKbZI3SiU7UafbHW3ddh+GnmMNmkk/by/tUiVXs1IjyZS37d6jzQ6o25kycgIkZL9JIsiUepEk2ykaP51IlxK++5AdEpH1RQZcSGEcHdFpwzDzJ7JoNrh54Hzgcfc/dfbGFsIIaogOWU4h4A/cPePAZjZX5jZHe7+uZbGF0KIUhzJKUNx9/s2NPWAU22MLYQQVVjraHxK65q4mb0cuNPd/7Lg3DXANQBM7mt5ZjtPGzW9t+q4rOvsy9pG1/Ee6dQrmduo66qk+TfpoNxO6nuho7EBp2SlOt3rTsLCriPnUTfVvqh/7r6564r7FN53RHmA7aKVeAXM7ArgCuC6ovPufhQ4CtCbO7ejj1QI0TWaTvYxs6uAVwDHGWwyf8OG8zPAzcDXgIuAI+7+YDh/HvAF4O3u/jvDxmrNiJvZYeCFwOuBp5nZBe5+rK3xhRCiDHdYbmhXCDObA24Dnuvui2Z2u5ld6e53h27XAY+4+01mdjHwXgb2ETPrAW8D7q8yXp2NkreMmb0A+AjwfcAngY8Bz2pjbCGEGM0gxLDKTwUuAx5298Xk9T3A4Q19DgPHANz9AeASMzuYnHsj8B7gG1UGa8ux+TlgfxtjCSFEXWrKKeeaWVwlH02k4JTzgJPh9XzSxqg+yYL3tLt/1syurTIZJfvsMI07KEdkLkaa2ah3eJZlqRNxizWxyxystcbeooOyTtZkFSdbYfZjA07J1seu4IgclVFa1rdfcjzqfFHfbeGwWn1rnxPufumQ88eBA+H1waStSp9fAL5uZm8CLgaebGan3P39ZYPJiAsh9jwNOzaPAReY2XQiqVwO3Gpmh4AVd58H7mAgu3w60cS/mLSvB32Y2bOB+4cZcJARF0IIHFhuKGXT3U8nUsg7zexR4EvufreZ3QQ8DhwBbgFuNrM3AxcCV8d7mNnPAt8NPMXMvuLuf1Q2noy4EEI4rDaYd+/udwF3bWi7PhwvAK8dcv37gPdVGUtGXAix53Fcm0Lsdba6AfFWMy/rbEoc2+tmWNYpAVtnc+G62Y/r868wdhMOypEZiBWyFfP3a3m8EmfkyLEL5lfFEVnHATlR6bpe6fXD2rdDQ2HirSMjLoTY82h7NiGE6DINa+JtIiMuhNjzNBmd0jYy4kKIPY/klD3GdrMw65aLHVUmNu/U236GZdPOyqL3XaX8bNFxHYfcpj4jMhfrOCUjVRyU6ZhbdUpu532nTsAqGY9Ffeo4Isv65Nps9P0mCq4rm2cjuLOmlbgQQnQTR9EpQgjRaSSnCCFERxnUE69eAWuckBEXQux5JKfsAbaabTkq47Hu3pap47Lu2KPKsFYpAVvooKzhHK2y/2Sd/R7byprcrlMy9il7r/2JGo7BCs6+ImdkHUdk7rowz6kwzyoOz6pjlI7dVilaJKcIIURncSrv2jN2yIgLIYQyNoUQors4MuJCCNFZ3GFpRdEpZx11Mi+rZGH2RjkGS7ItR2Ze1si2zLVNTI28ruje5WMMd1ZWcSL2+73C9lF7TUbH4CjnZ9NOybL3NVXgrByVoVjWp44jEoqzIquMV2fOtRyUVvyZ9IN/slcw/yp9m8BxrcSFEKKzSBMXQojuIk1cCCE6jGslLoQQ3UZGvMPU3udyRGnVXqnzcHO2Zd3sx5R+FadkUSnXydHX9ScmQvvWHJSjHIr9iZrlVPu9TW3lDtH1w+KszxrlW+NxmQOwSubi1ER/6PlRWZF1MxebcFBOpvtcjnA+wmgHZOn5krK0aXuVvk2w5s6iolOEEKK7aCUuhBAdRZq4EEJ0HNVO6Ri1qwAW6NWjtO+N16V96m6FVqeSYJnmnY4dNc1e0JFz7UFPrZMs0+sXt4+qQJgbb8T9yuY8KskmHtdJshkc9za15caoUeVv1NyqzLmKrl6kQU/2iz/XqDtP9Hub2mPfyYLzG/tE+uuaeGgrm2f0ZSTtk7kkp9HjbRUl+wghRIdR2r0QQnSYQbKPjLgQQnQTl5wihBCdRWn3HWJkRcASJ2GRo7E3MZmdL0m+KXI0VtkKrbTCYME2a7FvmeMyrfJXViUwXwWQrL2gf90t0rJEney+o+YJxY7EKkk2aWLNxvZCJ2FJAkmR87PMKVmnfaqCc3Gr1fzK+qTOyJyDsjTBpyDhpsShWMdBGf2QPcLvSc45yqb2/BjZ+X7Tjk2HlY4a8d7oLtvHzL7NzN5jZve1MZ4QQtQhXYlX+Rk32lqJ/wDwMeB5LY0nhBCVcXdFpwzD3f/QzH6ojbGEEGIrNLnKNrOrgFcAxwF39xs2nJ8Bbga+BlwEHHH3B83se4DrgC8AzwLudff/MGyssdLEzewa4BoAJvc1d98ROnhOfw6JOLHIVJH+XaZ956/brG1XKXpVqrEXJL3kdOSoucb2pOBUqSZecl3RzjejCkzFvnF+ZXpxPJ7OtW/WtnPX9Yu16NG75IzWwQv17JIEn6gjlyXUTPZ3vhBU1IyLkmRi22R/8xgb31d6WKRVQzVte/088XzxM7ACXT3qvfG9NquIN5t2b2ZzwG3Ac9190cxuN7Mr3f3u0O064BF3v8nMLgbeC7wQeBpwi7vfa2aTwHEz+8/ufqJsvFY08aq4+1F3v9TdL7WJmd2ejhBiD+FrXumnApcBD7v7YvL6HuDwhj6HgWMA7v4AcImZHXT3/+Lu94Z+K8DysMHGaiUuhBC7gTusVV+Jn2tm94fXR939aHh9HnAyvJ5P2qjQZz60/QLwG+7+rWGTacWIm9mLgJ8BnmZmbwZ+090X2hhbCCFG43j1Algn3P3SIeePAwfC64NJW+U+ZvbTwD53v3HUZNpybP4p8KdtjJVSVuCqSM/uTxZr2GVx4Kn+PTE1W3h+ZJx4SVz3KE059umXaOITkyWaeNI/9o0FpnJjjNg0oY6ODJnOXRbLHa+bnexvui5eW6alV9koYVSceC5GumBzhKgjl+ngRXHW8doqsdWRTBMffV15MazNbbnrCnTw2F62/0JeHy/rs/ke8XZlmv762GSG1dZWsg7ecCSJw2pz0SnHgAvMbDqRVC4HbjWzQ8CKu88DdzCQXT6daOJfTNoxs58D9rv7jcm5RXd/sGwwySlCiD2P09zfBXc/bWbXAu80s0eBL7n73WZ2E/A4cAS4Bbg5USYuBK4GMLMfB34T+IKZ/QTwFOB1gIy4EEIMo4acUuVedwF3bWi7PhwvAK8tuO5jwDl1xpIRF0KIeo7NsUJGXAghqBw+OHbsCSNe5HQsc2b2p7L49P505rjs5/rMbr4u9M3vHL/ZYThqR/dN7QXOyLwzs3i8nJOzwLFZxUlYlDhTtOvNxvaipJzpkvFmp4qdnEX3KHM+lhVEKkosKSsaNcpZmXdsFifOlLWvJ86UpKmMTJYpSbIpc1DGJJr1z62gDUY7GssSa8p214nNlkoUQXDOOSjXShyX6fHqSvH5hrdSc4fVVaXdCyFEZ9FKXAghOoyM+BhQZZPjVEbJySYh3jvKIjEOvKjPxFQWR15ae6SBjYajXJLeO0ol8fxUiTQxl7RPjahNsvE4kn5dL+tbJpekx1E2ifHgk7l7hPaCuOyyuiGR0k17bfP5OhLJ1ESxfFNF0iiaU5WS2Omti+qKDBuvqEZ4pZjr1dXN86wSe1fWJx1nrUxOWS1pHxzbaoW+DeDucmwKIUSXaTLEsE1kxIUQguaTQNtCRlwIsefxZtPuW6XzRrzKHpUxnDDVwidm96+3Re17Yia05/pk+neqQZeF9sV620V1SEbtRbnxfkX693TUl8PxgZnJwvZRYX6RWrVHSsIKi2qgzIT3NxGOZ2I98X6xLl2412SFQspFtbLLNPGoLwdpfl1Tzp3vjdaoy+qF1GG97klZXZQynXs1q15qS0tJW4VwvRgKmBy7jX7QOd284Dg/t6Btry4VztmXzgz+XQ5tK6HvWvO1U+TYFEKIzuKsSRMXQohuMiiAJSMuhBDdRHLK7lFFB5+YyfbrTGO8J4P2PbkvKxo2OTOXHU9PFB6nGnXZ3pZWoimn7UX7VgJMTBXHfkdtO4333h+07/1hbnMlWnmqV5el11ehX6ANl8V4R+091byj3j0X3l/UxHMadYyhT1TlGFtd5dtvUTp7WYz3REFsdWyfLOmbS2Ff2xxnXT65Eq25QJcu1ZGD1mxBM7agL5P0L425LtOzk/ltSxNP5uErQdtePJNNYyk79ni8vLTpOsL786Y1cVQASwghOou7s6baKUII0V20Em+ZVEYp22atrBrh1NxAOokSytRsdn56NpMpJqeDHBHllESmKEuTL59zKgkUVygskk2gWDo5MDMRzhcfz+TC9QbHZSnpkTIvfXptlDymy0IFJzZLJ1FCKZcuNo8HmSwSp1YlmqDoHkVSyaA9u24qllHwgfRgRWF7kJM0SkPf0pT+XPzpaDklDQvMjb26VHjMSnGfLFyvuO9alFbitobpcShHYP3ibQ/XQihglGpSuSRKKF4moSwubGpfW86ere/wStnrSGFjRGeNuBBCNIa7jLgQQnQVR0a8dYrklBidErMtUwkFYOrAIQCmZ6fX22b2ZddNzWaPZDrIGEWRI2UZj6sjtLWi6oIAs1NBFgnyzZPmsnmckxznIlJiJEtJdcBUTomV7lZLpImy+adRGGUSSi7iJGgTU8mgUzG7M7fpwuYxoDjLMc45Fiwq+6IdP6H1KJNYoTBmVa4sZseLQbJYTr7yl0gahdEkbIjq6A0+L++N/i9XWM0vSiUrQYI4cyo7DpLFamxPpY4K0R3WiyE4g/8XNjlZ2DeSy6wskEuiVLJ25vT68cqprO/Kmez5ry0N3vfqcrPVCktxZy3KTR2is0ZcCCEaQ3KKEEJ0GxnxFigqdjURN3EIST1lyTwz+2aSf7OviDNzU+F81j49XRz1kUZexK/+UYKIxyvhOP06H+WUGHlyIIx3ToGEAnAgkVz2x0iWyeLjKBsURaUsB21iOXy9LlODUnlmlGwC+ciY6YnNckqc20SVZJlEmlgJc1v17LoorVjJRgnpng452SRIE7a0UNy+POhfWqypRJroTWSfm/cTaaJXHN0RySXqJHNdW8jkkShHxPactFIgp+QKSMXx4pzCnC35v2UTJXJK+KzyUSZBIjk9eKbLpxZCW3Z+OSenZPNbW0rmXEX2aQBp4kII0WVcK/HW6RVtsxZW5fnVd5ZKP7t/cN3cgeDYDCvdJ+3PVuWjyrpGlkIM61KoSxxX5enKPTofi5yWAE8O8eoH42o9TbufKnZslm0hlhJDq+Pq+8zK8C3ZIPsGElfc8Xg2V1I2rMrTby5rwTm3nK2GKdsqrMAxODmRfW4TwUkYnYhxNZ+Lr05W1+uOSsCWw+o7rMRZDCvZZLWbSxEv+Q+f+7YYfjdtdvAtsWxVG9PL12K8dLLqzq24T89nfU+dXD9eXSh2GK6mK/HwO2qx7HFwXPZnQnDAXHJdybeHuLJfDavouOpeOjmY0/J8Nrfl0wuhb1iJL2TPYHV58wo85mVMzDRtujwXL98lOmvEhRCiKdydtRK5adyRERdCCHe8YJPoLtApI55LsU9iWItS6gGm92WOzdkD2VfEVEbZfzD7Wv6UIKEc2pe1x9T2oup/USpZWMp+AU5b8S9Deo943yihPClIKE8OEso5sf/0oH22xKE4mauQmI2dyijLYc62GlLLS3eI3yyd5CWUcBzn5KG63uITg3+Xsq/UOedd2eaGQSLxIKNkkw6ySbwsyjNBOukl0kluHsncANZOfjM7PrNZTonOuzKKJBQI8epBuohGI+cYLHBcrp38xnrb6qlszkvz2TxT6QLy8kYac52bZ9xlaV8256kDc0P7RklmlIQCsPjNk8k8Q9t8JqWtLGRzWz4TU+w3e9f7UdLcNzp2vS7SxIUQoqsoTrwdYkZmGr5VHkqY9Z3dn63i0hX4056UreAPhb5xVV5UjzsSV99VSO9XZfV9aLbYyTk3me5XGfeDHO7MBFhNluLxS0JZ/ahYy2uit3klHh2Y+8Lqux/C8npnMudb6jzMORTLVuJx9d3fvNrKOTDD70MuazKOk8vCTFatC9ncVr71WHaLU8FheDpzGK47NmNt60AMJSxcfQNraYbxSliJl9TYzjkuk28HS9/K5pOubiG/wl06ma3KVxY2Z2fGsLx+cIzH1XWvwGnfC0WvVkNmZnRKRsdlnN+Zbw7aF76R9V38VliJh9V3PE6Jc54I2dTNl42VERdCiM4y2J5N9cR3hFzIVlgR9JMd6uPu9LGMbEzamQur6/OSlfh5IcTwqUEfj+F8sS5I1IzT0Ly46h0VVghZwtDBmNQTjuPq+1BYrcfVbuFKPO6AE+YU66GsFPx+ruSyeoq19KIQwtzqezlbgcXVd28xW42loXu+GLTekjoV8TO22QPrx96fSM6H1fly8W4ypWGDpwa68mpYfa/F46CJr+TC4JL5h//k+RVi9vsTvzt41MdXCna4Cc8grvyj/r30zYH+vfDYt9bbFr/5RDjOrls+FUL+wrfEVF/uhW+Wk8HPkgs3DCGG6XuM69OYkJPTvr8RV9/Zszt9YtDnTFyJz4d7nA5hhQVfDSfDc55cDqGQk5u/MWwLRafsHEXOTMgcmsOeNDwAAArCSURBVLHQVc6Il8R+P/XA4D9WNNxPCXJKjLmeLvhqCbC4OjBsayXGOkovuezMxGDH4lXnzBTLJtFQ7g/HqfNwMpR8srUoTcQ5heJgifGLseE56aW3ObMU8n8sZhInZplskjPcC5slgRhnHetZ5+pVB8PXi3/AE8dm7iuvRSdtLFgVYtDDPFbTeQTDHQ36UjCOZdLE+hjhdyPWvI7GfXJflEvCH5O0LTozg5SzHIxjarzPBCN+5rEgrQQn4dKpYBALpL7oGPSD4f9S+B2Mf6jS4lNRblnJSSjZMyoy3AALJxaStuz8/HI2t4XoKA02PP21m+1njecshAzpgjjybeGKEx+JmV0FvAI4Dri739DW2EIIMQwHhRgOw8zmgNuA57r7opndbmZXuvvdde4THUgTiZwSd9zJlZENq9qn7N8cThhXveeUlHWdDKutovKsi/FraMlKfLagxknMtow1UKLDMHdcELoXJYNcxmOkHxzB6ebOBTvdAPSDnDJZ4MyEbPVvQUKxpWw1FlffOckiySrMhejFVU/8xhBWgjk5Iumf25A3bhIcV+JhTqunCr4RBNnkTIlMkc94TFek2ZwnZrLfqbgqn1wu3ox4fYUbd70JWZhRvokSSTqnOLfoJIwyRQzRWwsrcUs+w5jlGOWI+E0ihiOm3yri+ZhtuTifHZ8pmVO6Av9GuO+3wip6ISy/o5yS1tmJq/PZsIKfrRlUMJKGo1NGLVrNbAa4GfgacBFwxN0fTM69Cng+AyXrK+7+7mFjtbUSvwx42N3T7373AIeBkUa8dDf7VBOPKeclBaRiBEjaHrXvshT2qKbEYlGp1rywXLzjeWQqVwgq0bNjdEcw4tMx5npis/GEzEDFCIyyMJOc4p3IKT3L3ndOLIoGvSQ6ZT1tPejMvXAcU8DjcRptkZMP4u7tIdbfY2GjWEM7Md5elqK/UlLPOsRcp/OIkR5R143HKyHuebUgzjrKDr2SSI/cH6qVzbu3x3nGOOsY9ZHGgZcZyTI5JRdxMrU5bT7KEUUx2ZAZ75VcPHicczbe0hNZnzi/VDp5fCkb74ngoDlT4kxcWkuku+h3CguptZI5b53mjHjFRet1wCPufpOZXQy8F3ihmZ0PvAF4vru7md1nZp9w978qG68tI34ecDK8nk/aRlJUuRCy2ikT4Rd0IqwuYt2TWIEwrQJYVvmvzJkZ3Tu9kjojKf2S/SMnCzZVKEvUye3xuJStwlLjHUPncsYsbjqwGjTQNMQwlwxUkuAT96CMmyYkq93cpgTRiJdU0VurkSzjsfLfiP9UcVUeE3zWQrhervpfcpyrqBdrd5QY0lRfjuF3vZIEmLWSr+SFK/Gwwl0tNZSDzzka6KXgwIztMXEmOrX7icHrh9/zqJnHcL21+MdpfZ6hrku4Lq78o0E/E9pTgx0N96nVuBIv8dEkU42GO7chSNORJM06NqssWg8DvzoY2h8ws0vM7CDwYuBznpXlPAb8CLDrRvw4cCC8Ppi05TCza4BrkpdPLP/5+x9bhhPp+egaejz595GSAT+19bmWcS5hLrvMOM0Fxms+mksx4zQXGDWf9O9AWKtwKhznr7xgu5PxhcfuXP7z959bsfuMmd0fXh9196PhdZVFa1mf2gvetoz4MeACM5tO/jpdDty6sVPyINYfhpnd7+6XtjTHoWgu5YzTfDSXYsZpLjB+83H3lzR4uyqL1rI+x4ELN7R/edhgDQdbFuPup4FrgXea2Y3Al+o6NYUQoiOsL1qT15cDd5jZoUQyAbiDgexCool/0d3ngTuBF1imdV4G/NGwwVoLMXT3u4C72hpPCCF2A3c/bWbpovVRkkWrmd3EQAk+AtwC3Gxmb2aw8r46ufZvzOxm4LfMbBV4zzCnJox/ss/R0V1aQ3MpZ5zmo7kUM05zgfGbT6MULVrd/fpwvAC8tuTaDwEfqjqWeVkVJCGEEGNPK5q4EEKInWEs5JTtZDc1PI9nAjcCnwfOBx5z91/f0Oc1wL8C0mDe97r7B5ueSzLWZ8I4q+5+5YbzrTyXZKzvYBDn+tdJ00EGWt9rQp/XsEPPxsy+jcFnc4m7f0/SdoiBvvh/Gbz/X3X3vyu4tlYG3Bbn8lvAaeAJ4BLgOnf/esG1XwW+mrz8mrv/8+3MZch83gr8UOj2tuQr/sZr23g2dwD7QrfvBp7u7mc2XPtVGn42ewJ339UfYI5BCM108vp24MoNfd4EXJ8cXwx8eofm8j3Aj4fXfwG8YEOf1wDf0dKzeeuI8608l+T+TwGuCq9vAH6grWcD/BPgpcD9oe024CeT45cCHyy47nzgz8mkw/uAi3ZgLjeG4zcC79rKZ9rgfEaO0+Kz+alw/J3Au9t6NnvhZxzklLLspshhBmE7uPsDwCUhVKcx3P0+d/9YaOqRTzFI+QUze4OZvSVZDe4UF5vZG83srWa28ZlAS88luf9j7v7HAEno1KXu/j8Luu7Is3H3PySfBAHh/VP8ewPlGXCNzsXd3xxe9hisyIv4QTO73sz+jZl9/3bmMWw+AGb2a8ln8cYkFXwjbT2bj4SXvwi8q+Tyxp/NXmAc5JTtZDfNs0OY2cuBO939Lzec+lPgDnd/1Mx+FPgocOWmGzTDO9z9XjPrA58ys5PuHpNRW38uCT8NfLigvc1nA/n3Pw882cwm3H2lpE/ar1LJh61gZk8C/jHwypIub0o+0zng82b2Y+4+NJlji3wU+Kq7nzKzn2dgOK/e0KftZ3MQ+HZ3/18lXdp6NmcV47AS3052045gZlcAVwC/tPGcuz/k7o8mLz8BvCgxso3j7vcm/64Cn07mFGn1uQT+KfCRjY1tPpuE+P4PAt/YYMA39kn77cgzMrNzGGQi/6y7P17UJ3ympxlIGZfvxFzc/X+7e/ot8hPAPyro1vbvz9XA+8pOtvVszjbGwYhvJ7upcRLZ4sXA64FvM7PL4lzM7O1mln6DuQh4KDGyTc/j2WYWV04XAV/erecS5nUF8Gfuvpy8bv3ZBNbfP8nvTTKPnpl9e9JeOwNuK5jZucC/B37F3R8ys1dunIuZXWlmMb37QuArTc8lGevfhpcXkaRu78azScdl8P/qjti2G8/mbGMs4sTN7IcZOEQeBZbd/YY0u8ndj5jZLIMojL9l8OH+hu9MdMoLGEgCaXGbfQz+Yz4nzOX1wN8HHmLgTLzF3T+zA3N5ejL25xmskCaBX2YQjdHqc9kwrw8Dr3P3E8nr+Dnt2LMxsxcBrwZeAvwu8JvALPAO4GHgmQy+jv+dmT2PgZPz4uTaVwGXMojAeNC3H4FRNJd7GMiT6Qr8pLu/NM4l+UP7VuBzwNMZRGC8fTtzGTKftzAIGjjO4LN4i7s/uBvPxt0XzOwngPPd/XdC3x1/NnuBsTDiQgghtsY4yClCCCG2iIy4EEJ0GBlxIYToMDLiQgjRYWTEhRCiw8iICyFEh5ERF53AzP4kqZIohAjIiAshRIeRERdjj5m9HXge8KZkRV5UrVCIPYkyNkUnMLM/AT7g7h/Y5akIMVZoJS6EEB1GRlwIITqMjLgQQnQYGXHRFU4Cc2Z20YZa2ULsaWTERVd4H4ONOn4f+G+7PBchxgZFpwghRIfRSlwIITqMjLgQQnQYGXEhhOgwMuJCCNFhZMSFEKLDyIgLIUSHkREXQogOIyMuhBAdRkZcCCE6zP8H2auy6mdg0IkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.pcolormesh(tmesh.T, zmesh.T, c_array.T, cmap='RdBu_r')\n",
    "plt.axis(pad_limits(tmesh, zmesh))\n",
    "plt.colorbar()\n",
    "plt.ylabel('z')\n",
    "plt.xlabel('t')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((63, 64), (64,), 63)"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c_array.shape, z.shape, len(t_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f55e9120a58>]"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD3CAYAAAA0Vx7KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3zV5d3/8dd1RvbJ3jtkAIGwtwLKUkSljrpr6+3obbV3e9f2V+1ttfWuo629W+24W+vd1lW1ddWBgIhIVFYIeyZkELL3OFlnXL8/viGAIgQyzknyeT4eeUDIyTmf8xXfXLm+n+u6lNYaIYQQ3svk6QKEEEKcmQS1EEJ4OQlqIYTwchLUQgjh5SSohRDCy1kG40kjIyN1amrqYDy1EEKMSNu3b6/TWked7muDEtSpqank5eUNxlMLIcSIpJQq/bKvydSHEEJ4OQlqIYTwchLUQgjh5SSohRDCy0lQCyGEl5OgFkIIL9en9jyl1Gags+dTl9Z68eCVJIQQ4mR97aNerbX+yWAWAvDCNxfSGZXOVf/+NEfya5i0KBH/IJ/BflkhhPBqfZ36yFFK/VAp9ROl1IrTPUApdZdSKk8plVdbW3texWTk1RBycBcNFW3krSqhy+48r+cRQoiRpK8j6p9rrbcqpczARqVUq9Z648kP0Fo/AzwDMGPGjAE5jUAONRBCiD6OqLXWW3t+dQG5wMWDWZRSqud1B/NVhBBieDhrUCulximlbj/pjzKBwsErCZSpJ6jdktRCCNGXqY8W4HKlVDwQDJQBLw9mUaaeoHZLUAshxNmDWmtdAVw1BLX0vCCYzD1B7ZKgFkIIr1rwoo18xmTpCWqn24PVCCGEd/CqoD7OYjXKcklQCyGENwa1xmwxA+B0SFALIYRXBbUGlAaLj1GWs1uCWgghvCqo6ZmjPhHULg8WI4QQ3sGrgvr4iNrqazSjOLokqIUQwquC+viI2upnzFFLUAshhJcFtcbIaovVhDIpujtkUyYhhPCqoEYB2tjrw8fPLEEthBB4WVCfvA7RN8BClwS1EEJ4WVArerfM8w2w0tUuQS2EEN4V1BhdHwB+gRY67Q6P1iOEEN7Au4Ja9TZ+4BtopbNNgloIIbwuqI9PVPvbfOiQoBZCCO8KahSonjlq/yAr3R1O2ZhJCDHqeVVQnzz1ERBsnD7e3tLtuYKEEMILeFdQc2JEHRDiC0B7swS1EGJ086qgdpvonaMODDFG1PbmLs8VJIQQXsCrgtqYozZ+GxhqjKjtTRLUQojRzauC2n1SUAfYfDCZFW2NnZ4tSgghPMyrglorMPWcPK5MiqAwX1obZEQthBjdvCqoTx5RAwSF+dHWICNqIcTo5lVBrU0nuj4AgiP9aKnr8GBFQgjheV4V1G4FppPWt4RE+WNv7pYjuYQQo5pXBfUXR9T+ADTLqFoIMYp5VVC7TaeOqENjAgBorpGgFkKMXt4V1J+7mRgabQR1Y5XdQxUJIYTneVVQa9OJ9jwAH38LgSE+NFW1e7AqIYTwLK8K6s9PfQCExgbSIEEthBjFvC6ozZ9r8IiID6Sh0o5269N/kxBCjHBeF9SmzwVyRGIQzi4XLfVyQ1EIMTr1OaiVUv5Kqd1KqScHqxh9mqmPiPggAOqPyQ1FIcTodC4j6p8BOwarEACXGcyuU0fU4QmBKAW1Za2D+dJCCOG1+hTUSqmvAZ8CxYNZjDaB+XMjaquPmbC4QGqPSlALIUanswa1UiobGK+1fuMsj7tLKZWnlMqrra09r2JcZoXlNKvFo1ODqSltQWu5oSiEGH36MqK+CuhUSt0PXAjMUkp99/MP0lo/o7WeobWeERUVdV7FaDNYXHwhkKOTbXS0OmiVnfSEEKOQ5WwP0Fo/evz3Sik/IEhr/ZvBKMZtVoBGOxwoH5/eP49JCwaguriF4Aj/wXhpIYTwWufS9XENsACYo5S6cTCKcZmNX3X3qQfaRiYGYfE1U1nYPBgvK4QQXu2sI+rjtNavA68PYi3o40Hd1QVBQb1/bjKbiE0LpvJI02C+vBBCeCXvWvBiUUBPUH9OXHoI9cfa6OpwDnVZQgjhUd4V1GYjqN2dX7xpGJ8ZitZQWSCjaiHE6OJVQa0tRjnuji8uF49ND8FsNVF2oGGoyxJCCI/yqqDunfo4zYjaYjUTnxkqQS2EGHW8KqhPjKhP3y+dNC6cxqp22hqln1oIMXp4V1Afn6O2n34DpqTscACO7pNRtRBi9PCuoLYa/XkO++n39YhICMQW7kfx7rqhLEsIITzKy4LaaOvubms57deVUqROjqTsQAOOrtNsCiKEECOQlwV1z4j6S4IaIG1yJC6HW24qCiFGDa8KaswWXOrMQR2fGYpvgIWinee3Q58QQgw3XhXUZmWl3RecLV8e1GazibQpURTtrMXZLdMfQoiRz7uCWlux+4Gz9cybL2XNjMHR6aJ0b/0QVSaEEJ7jdUHd7gfuljOf5pIwNgz/YB8Ob6seosqEEMJzvCqoTfjQ7qvQZwlqk0mROSOakj11dLU7hqg6IYTwDO8Kau1Lmx/o1razPnbs7FjcTs3hrTKqFkKMbN4V1PjRGgCm5rMHdXRKMFHJNvbllstZikKIEc3LgtqfFn8wt7Sj3e6zPj77wnjqy+1Ul3x5l4gQQgx3XhXUmgBa/RVKa9xnaNE7LmtWDBZfM/tyK4agOiGE8AyvCmqX8qXL35jGcDacfeWhj5+FrJkxFG6rpqOt+6yPF0KI4cirgrpL+eEI6Anqur5tvDRpUSJOh5u9H5cPZmlCCOExXhXUncoPl78xN+2s7dsS8Yj4IFImRrBnwzFZqSiEGJG8Kqi7lB/qHIMaYOrSZDpaHRzcXDVYpQkhhMd4XVD7Wt04LOcW1PFZoUSn2Ni57ihut7TqCSFGFq8K6g7lT6jbTZPNhLOq7wtZlFJMXZZCc00HBVtlVC2EGFm8KqjbTYGEu11UB2u6K87t5mD61Cgik4LY+m4xLtfZe7CFEGK48Kqg7lABRLhc1AZDd/m5BbUyKWZfOYaWuk4OflY5SBUKIcTQ86qg7lT+hLrc1IWAu7YO3X1uvdEpEyOISQsmb1UJTod0gAghRgavCmqUIshtpSZEgdY4Ks9tZKyUYs7KMbQ1drFng/RVCyFGBu8KaiDQ7U9VmAKgu7T0nL8/cVw4yRMi2PZeMfbmroEuTwghhpzXBbWvtlEZbvy+u6TkvJ5j/nWZuBxuNr91ZOAKE0IID/G6oHaqELS/ojvAh67i4vN6jtCYACYvTuLgpiqqis58rJcQQng7rwvqNmUj1g2N0f50F5ec9/PMuCyVgBAfcl89jJZFMEKIYeysQa2UMiml3lNKPaSU+plS6mWllP9gFdRqCiHW6aAiQtF95PynLnz8LMy7OoOa0lb2fHxsACsUQoih1dcR9Sat9SNa6weBAODqwSqoyRxGUlcnB8M7cNbW4qw//5PGs2bFkDwhnE1vHqG5tn0AqxRCiKFz1qDWWru11j8DUEpZgETg0Ocfp5S6SymVp5TKqz2HfTo+r9EUTpLTyeFI49DazoMHz/u5lFJcfMs4TCbF+ucPyhSIEGJY6vMctVLqEuBd4F2tdd7nv661fkZrPUNrPSMqKuq8C2oyh5PscFAabbTodfUjqAGCwvy44KuZVBQ0sXej9FYLIYafPge11nqN1vpSIE0p9a3BKqjJFEaS00lbgMIRFUrngf4FNcD4eXEkTwjnszcKZQpECDHs9OVmYrZSasVJf1QMjBmsghrMkSQ6nFgw0ZgcSueePf1+TqUUF908DrPFxNpn9+FyyKZNQojhoy8j6i7gdqXUg0qpnwLjgScHq6A2ZcPqE0SqOYDCJAvdpaX9uqF4nC3cj0W3jqemtJVP3ygcgEqFEGJoWM72AK31EQaxy+MLlILQZDK1mW1RrcwDOnbswLZkSb+fesyUKCYvSmLX+jISMkNJnxbd/3qFEGKQed2CFwBCU8js7GRraB1YrbTn7xiwp557dTrRqcGsf/6AzFcLIYYFLw3qZDJb63FYFO6xaXTk5w/YU5stJi65YwLKpFjz53045EBcIYSX886gDk9jvN3Yo6NmbBQde/fiam0dsKcPjvRnyTeyqS1r5cO/7Zf+aiGEV/POoI7MIsblIsYnlPw0DU4n9s2bB/QlUidFMu/qDI7k17LlnaIBfW4hhBhI3hnUUeMAmOQbyYfBxzAFBWHfmDvgLzNlSRLZF8Sx/f1SDm6W47uEEN7JO4M6OB58bEx2mynrqMAyazptn3yC1gM7RaGUYsFNY0kYG8ZHLxykoqBpQJ9fCCEGgncGtVIQmcnk4/PUOfE4KyvpKigY8Jcym01cetdEgiP9ef+Pe2iosA/4awghRH94Z1ADxOaQXXkYf4sfn6Z1g8lE65q1g/JSfoFWLr93Eiaz4u2ndtBc2zEoryOEEOfDe4M6fio+nU3MCJ/Aho6dBMycScuqVQM+/XFcSFQAV35nCk6nm3/9ZgdtjZ2D8jpCCHGuvDeoE6YBMM8nitKWUlyL5tJdXNzv3fTOJCIhiCv/Ywqddgf/+s1O2lu6B+21hBCir7w3qKOzwezLvE7jJPGdE3zBYqFl1arBfdmUYC6/dzJtDZ28/fROOu2OQX09IYQ4G+8NarMVYnNIqz5EbGAsH7fsIHDeXJrfew/tGtzVhPEZoVx29yQaq+y8/dROOtpkZC2E8BzvDWqAlHmoY9tZlLCATys+xX/lCpwVlbRt3DjoL52UHc7yb+bQUGnnzV/toK2xa9BfUwghTse7gzptIbgdXOIbS5eri62ZCkt0NI0vvzwkL5+aE8kV3zamQd781XbpBhFCeIR3B3XKXDBZmFJXRnRANGuOrSP0uuuw535C99GjQ1JCQlYYK/9zKl0dTt54cjv1FW1D8rpCCHGcdwe1TyAkzsRUksuylGV8Uv4J1q8sB5OJxldeHbIyYlKDueo+owvlzV/lU13SMmSvLYQQ3h3UYEx/VO7k0ti5ONwO1nfsxLZkCU2vvYarbehGtxHxQVz9/en4+lt469c7KNp5/ietCyHEufD+oB53GWg3k2qLyQjN4JWDrxBx5x24W1pofPHFIS0lJMqfq38wnfDYAN7/0x7y15YO2gIcIYQ4zvuDOnYShKagDrzD9WOv50DDAQpjNEEXXUT9X/82pKNqgMAQX666bxoZ06LZ9MYRPnrhIC6nHJYrhBg83h/USkH2lVC0gcvj5xNgCeDVQ68Sec89uJubaXzxpSEvyeJjZtntE5hxWSoHPqvk7ad20tkmC2OEEIPD+4MaIPsr4HYQVJzL5WMuZ3XJaroyEwlauJCGv/51yEfVAMqkmH3lGJbclk1VcTOv/TyPxirZeU8IMfCGR1AnTIeQJNj9KtePu54uVxcvH3qZyHvvwdXcTP0zf/ZYaWNnx/KV/5xGd6eTfz6eR0FetcdqEUKMTMMjqJWCqV+DI+vJ0lYuTrqYF/a/gHNsKiErV1L/17/SVVzssfLi0kO47kcziUgIYu2z+9j48iFcDpm3FkIMjOER1ADTvgbKDNv/xt2T76a1u5WXDrxE9Pfvw+TjQ/Vjj3u0AyMozI+v3DeVKUuS2PNxOW88uZ2WOlnJKITov+ET1MHxMHY57HiR8SHpXJx0Mc/vf56OED8iv30v9txc2tav92iJZrOJC67NZPm/59BU08E/HttG8S7ptxZC9M/wCWqA6bdBex0ceLt3VP3igRcJv/lmfDMzqH7scdwdnh/FjpkSxXU/mklwpD+r/ncPn75eKFMhQojzNryCOn0RhI+Bz55mfPg4Ficv5q97/0p1dz0xP/4xjvJyan75pKerBI4vjpnGxIUJ7PzgKP98Io+6Y7JPiBDi3A2voDaZYP73oXIXHF7N92d8H7d282TekwTOmkX4179O49//TtvHH3u6UgAsVjMLbxzLinsm0d7azT8f38b21SW43bKaUQjRd8MrqAEmXQ9hqbDhCRKDErgj5w7WlKxhU8Umor73n/hmZVHxXw/ibGjwdKW9UnMiufGhWaRNjmTzW0W8+WQ+TTXtni5LCDFMDL+gNltgwQ+gciccXsNtE28jyZbEY1sew2UxEf/LX+Jubqbyxw951T4c/kE+XHLnRJb+WzaNVXZe/dlW9m4s96oahRDeafgFNZw0qn4cX2XlgVkPUNJSwt/2/Q2/sVlE3fc92j78kKZXXvF0padQSpE1K5YbfjybuIxQPv77Id757S5p4xNCnNFZg1opla6Uelkp9QOl1FNKqYeGorAzMlvhoh8Zo+qdLzI/cT7LUpbxh11/YF/9PsJvvZXABfOpevQx7Fu3erraLwgK8+WKb09m4Y1ZVB1p5uWfbmH76hJcLukMEUJ8UV9G1OHAK1rrX2qtvwPcoJSaPsh1nd2k6yDlAvjgIbDX89Dch4jwi+D+jffT4eok4Ve/wic5mfL/+A7dZWWervYLlFJMXJjITT+ZTfLECDa/VcQ/Ht1G5ZFmT5cmhPAyZw1qrfU2rfW/Pvc9nt99SCm47EnoaoUPf0KIbwiPz3+c0pZSfrHtF5htNpL+9w+gNWV33+2RjZv6IijMj+XfzOGyb02iu9PJG7/czkcvHaTTLrvxCSEM5zRHrZS6ClijtT54mq/dpZTKU0rl1dYO0Wq8mGyYczfkPw9l25gZO5Pbc27n9YLX+aD0A3xSUkh46im6S0opv+8+tMs1NHWdh7RJkdz40GymLEniwKeV/P0nmzm8tUpuNgoh+h7USqmLgYuB/zzd17XWz2itZ2itZ0RFRQ1UfWe38H6wxcM73wFHJ9+a8i0mRkzk4c8eprSllMA5s4l98L+wf7yRqv/+b68OPh8/Cxdcm8lXH5iBLcKfD/6yn3/9ege1Za2eLk0I4UF9Cmql1ArgEuA7QKxSau6gVnUufIPgiqegZh988GOsJiu/WPgLzMrMvR/eS3NXM2E33EDEHbfT9Mqr1DzxhFeHNUBUko1r/t90Ft6YRX25nX88to31LxzA3tzl6dKEEB7Ql66P6cCrwBzgI+BfwNhBruvcZC2DOffA1mfgwLsk2ZJ46uKnKG8r574N9+FwO4i67z7Cbv0aDc89T+3//I/Xh7XJZNxsvPmROUxenMShTVW89NBmtq8uwenw3ikcIcTA68vNxO1a6yCt9UU9HzO11n8bgtrOzZKHIW4K/OseaD7GtJhp/GTeT9hStYVHNz8KQMwDDxB6w/XU//lZ6n73ew8X3Dd+gVYuvDaTGx+eTeK4MDa/VcTff7KFgrxqr//HRggxMIbngpfTsfjCtX8BtxNevwNcTq5Mv5I7c+7k9YLXeX7/8yiliH3oIUKuvpq63/+euj894+mq+yw0JoDL7p7Eyu9OwcfPwtpn9/Hmk/lUF7d4ujQhxCAbOUENEJEOl/8ajm6C1feD1tw79V6WpizlV3m/4s2CN1EmE3H//QjBV1xB7a9/Tc2TT6Ldw2ehSeK4cK77r5lcfMs4mmraee3nebz/xz3Ul3tn+6EQov8sni5gwE26Dqp2w2e/hbAUTPO+zePzH8fusPPwZw9jUiZWZqwk/vHHMAUGUP/s/+GoqibusUcx+fh4uvo+MZkU2RfGkzE9ml3ry9jxwVGKdtWSNTOGmZenERod4OkShRADaOQFNcCSR6CpDNY+CCFJ+E74Ck9d/BTfXv9tfvzpjzEpE1ekX0Hsww9jjU+g9n/+B2dtLYm/fRpzcLCnq+8zH38LM1ekkbMwkfy1pez56BgFeTWMnxfHjMtSsYX7ebpEIcQAGFlTH8eZTHDVnyBpNrxxFxzdgp/Fj6cXPc2s2Fk8+OmDvFf0HkopIu+6k/ifP0H79u2U3nwLjqoqT1d/zvyCrMy7OoNbfjaXiQsSOLipkpce2kzuPw7T3tLt6fKEEP00MoMawOoHN7wMIYnw8g1Qewh/iz9PL3qa6THT+dEnP+LdoncBCFm5kuRn/oSjooKS62+gY98+Dxd/fgJDfFlwQxY3PzKHrFkx7PnoGC88+Bmb3iyUwBZiGBu5QQ0QGAG3vAYmC/xtBVTvJ8AawO8W/Y7pMdN5IPcBnt3zLFprAufNI+WlF0EpSm+8icZX/zFs29+CI/xZdOt4bvrJHNImRZK/9ijP/9dnbHz1MK0NnZ4uTwhxjkZ2UINxxuJtq06EdeVuAqwB/HHJH7ks7TKeyn+KRzY/gtPtxG/cONLefIOAmTOpevhhKu+/H3f78D2JJTQmgGV3TOSmh2eTNTOGfR+X8+KDm/jw+QM0Vnl+Xy0hRN+M/KAGiMyEb7wH1gB47gooz8fH7MPj8x/nzpw7ee3wa9y7/l7sDjuWsDCSnvkTkd++l+a336Hk+uvpKiry9Dvol7DYQBbdOt6Yw16YQOG2av7+0y2sfmav7CMixDAwOoIajB7r21aBXzA8vxLKtmJSJv5j2n/w0NyH2FyxmW+s/gY17TUos5moe+4h6dk/46yrp+Tar9L83nuefgf9Zgv3Y/71WXzt0XlMuySFsv31/OPRbbz7u11UFjZ5ujwhxJcYPUENEJYCt70PgZHw/Ffg0PsAfDXrq/x20W852nKUG9+9ke3V2wEIuuAC0t58wzgw977vU3H/A7hahv9KwIBgH+Z+JZ1bH5vH7JVjqC5p4Y0n83nt53kU5FXjlpNmhPAqoyuowegCue19iMqCl2+ET58GrZmfOJ/nlz+Pn8WP29fczl/2/gW3dmONjSXlheeJ+OY3aX7nHYouv4K2jz/29LsYEL4BVmYsT+XWx+Yx//osOtscrH12Hy88uIn8NaVyeIEQXkINRmfDjBkzdF5e3jl/3/V/2gTAq98cgl1Uu9vhrbth/1sw9RZY8Wuw+NDW3cbDnz3M2tK1LEhcwKMXPEqoXygAHXv2UvmjB+gqKCTkqquIeeD+YbVA5mzcbk3p3np2fVhG+aFGLD4mxs2NY9LFiYTFBnq6PCFGNKXUdq31jNN+bdQGNYDbDRseh42/gJQL4foXICAcrTUvH3yZX+b9kij/KJ5c+CSToiYZ39LdTd0f/kD9n5/FEhFB3H8/QtDChUNT7xCqO9bKrvXHOLy1CrdTkzIxgsmLk0gcF4ZSytPlCTHiSFCfze5/GtujBsfDDS9BzAQA9tbt5b4N91HTUcN3p32XW8bfgtlkBo6Prn9EV0EBIStXEv3/foAlImJo6x4C7S3d7N1Yzt6Pj9HR6iA8PpCcixLJmhWDj9/I3IFACE+QoO6Lsm3w6s3Q2QyXPAozbgelaO5q5sef/piPyj5iStQUHrngEdJC0oCTRtfP/h8mf3+i7r2HsJtuQlmtQ1//IHM53BTkVbNrfRl1ZW1Yfc1kzYphwoIEopJsni5PiGFPgrqv2mrhrX+HwnUw/gq48rfgH4bWmneL3uWJrU/Q5erinin3cGv2rb2j666iIqofexz7J5/gk55OzI8eIOiCCzzzHgaZ1prq4hb25ZZTkFeDy+EmOjWYiQviyZgRg9XH7OkShRiWJKjPhdsNm38P634Ctji45llIngNAbXstP9v8M9aXrWdS5CQeueAR0kPTASPA2j7aQPUTT+A4epSgJYuJ+eEP8UlK8tx7GWSddgeHNlexL7ecxqp2fPwtjJ0Ty4T58UTEB3m6PCGGFQnq81G+HV77N2O71IsegPnfA5MZrTWrS1bz2JbHsDvsfGvKt/j6hK9jNRnTHe7ubhr+9hx1f/wjOJ2E/9ttRN55J6bAkds1obWmsrCJvRsrOLKjBrdTE5cRwoT5CaRPjcIio2whzkqC+nx1tsB734M9/4SEGbDydxA9HoD6jnoe2/IYa0vXkh6Szg9n/ZC58SfqdlRXU/OrX9Hy9juYIyKIvOtOQm+4AZOvr6fezZDoaO3mwKZK9udW0FzbgY+/hcwZ0YybG0dMWrB0jAjxJSSo+0Nr2PMarP6hEdwLvg8Xfg8sxmkwHx39iF/m/ZKy1jIuTrqYH8z4AUnBJ6Y7OnbtouY3v6F902YsMTFE3n03oVdfhRomp8mcL+3WHDvcyMFNlRTl1+J0uAmNCWDc3FjGzo4lKEwONRDiZBLUA8FeZ5zDuOefEDXeuNGYNBOAblc3L+x/gWd2P4PD7eDW7Fu5c9KdBFpPTHfYN2+h9qmn6NixA2tiIpH33kPIFVegzCN/WqC7w0lhfg0HN1VSWdgMCpLGhTFubhxpU6LkBqQQSFAPrMNr4N3vQUs5zP53WPQg+Bo3zmraa3gq/ynePvI2Uf5RfHf6d7l8zOWYlLFSX2uNPTeX2t88Ref+/fiMGUPUt+/FdsklKNPoWM3fXNvOwc1VHNpcRWt9Jz5+ZjKmG1MjsekhMjUiRi0J6oHW1Qrrfgrb/gwhybDsEcj+CvSEzJ7aPTyx9Ql21+0mOyKbe6bcw/yE+b0hpLWm9YMPqH36aboLj+CbmUnEHbcTfNllI7IH+3S0W1NR0MTBTZUU7qjF2eUiOMqfzBnRZM6Mka4RMepIUA+W0k2w6vtQvRdSLoBLn4C4nqXm2s17Re/x+52/p7ytnElRk7h3yr3MiZtzIrBdLlpWraL+mWfoKijEEhdHxDe+Tui1147oLpHP6+50UrSjlsNbqzh2sBGtITw+kMyZMWTOiCYkSk5VFyOfBPVgcrsg/zlY/zNob4Bpt8KiH0NQFAAOl4O3jrzFM7ufocpexfSY6dwz5R5mxs7sfQrtdtO2cSMNz/4f7Xl5mEJCCLvxBsJvuQVLZKSn3plHtLd0cyS/hoK8amM+G4hODSZzRjQZ02MIChvZXTNi9JKgHgodTfDxL2Drn4yTZBb8wJjD7ukO6XZ183rB6/x595+p7ahldtxs7p1yL1Oip5z6NDt3Uv9/f6F13TqU1UrI1VcRcdtt+KSkeOJdeVRrQyeFeUZo1x5tBQXxGaFkzowhfWoU/raR3TkjRhcJ6qFUVwBrfgQFa43zGhc9CNlXQc/Nwk5nJ/88/E+e3fMsDZ0NzIufxzcmfOOUKRGArqJiGv76V5rfegvtdGJbsoSwm28mYPasUXnDram6nYK8agq2VdNY1Y4yKZLGhTFmahRpk6MICJbQFsObBLUnFKyDtQ9C7QGImQgX/xeMXd57w7Hd0c4rh17h+X3PU99Zz9iwsXx9wte5NO3S3lWOAI6aGhpfeJGmf/wDV3MzPq7yD+wAABruSURBVBnphN14IyErV2IOGn033LTW1JfbKdhWTWF+DS21Hb0j7TFTokibEklwhL+nyxTinElQe4rbBXvfgA2PQUMRJEw3RthjLu4N7C5XF6uKVvHcvuc40nyE6IBobh5/M9dmXUuwz4lDCdydnbSsep/Gv/+dzr17MQUEELzySsJvugnfzExPvUOPOh7aRTtqKNpZR315GwBRyTbGTI0ifWqUHHgghg0Jak9zOWHXy/Dxz6G5zOgQWfRjSDnxPt3azafln/LcvufYUrWFAEsAV2dezS3Zt5AQlHDK03Xs3k3jS3+n5f330d3dBMycSdjNN2FbvHjUtPedTlNNO0U7aynaUUt1sXG2ZVhsQE9oRxOZFDQqp43E8CBB7S2cXbD9Och9EtqqIX0xzL8PUub1jrABDtQf4Ln9z7GmeA0azeLkxVw39jpmxs7sXTwD4GxspOm112h6+RUcFRVYoqMJueZqQq++ekTv2tcXbY1dFO+q5ciOWioKmtBujS3cj7TJkaTmRBKfGYrZOjoWGYnhoV9BrZSKBX4GTNZazzzjg3tIUJ9Fd7uxWOaz34K9FhJnwoX/CVnLe286AlTZq3jpwEu8UfAGLd0tJNuSuTbrWq5Mv5II/xOnyWiXi7aNG2l8+WXsn3wKbjcBs2cTes3V2JYuxeQ/uudsO9q6KdldR9GOWsoONuJyuLH6mknKDic1J4KUiZFyM1J4XH+D+lqgC3j4y57k8ySo+8jRATtfMk5CbyqFyLFwwXcg56u9bX1gdIp8UPoBrx1+jfyafCwmC4uTF/PVrK9+YZTtqKqi+a1/0fTGGziOHsUUFETwihWEXnM1fjk5o/5Hf0e3i/KDjZTsqaNkTz32pi7A6NVOzYkgdVIkkYkyRSKGXr+nPpRSFwFPnimolVJ3AXcBJCcnTy8tLT3nQkddUB/nchqnoX/ya2OVY3ACzL3XWDzje2pnx5GmI7x2+DXePvL2mUfZWtORl0fTa6/TsmYNurMTn4x0Qq++hpCVV47I8x3PldaaumNtlPaEdnVJC2gIDPU1QjsnkoRxYbJplBgSQxLUJ5MR9XnS2jgG7JPfQOkn4B8G074OM2+H0ORTHtrl6mJtydpTRtkXJV7E5WMuZ37ifHzMJ0bkrrY2Wlatovn1N+jYtQssFoIWLCDkissJuvhiTH6y5SgYqyJL99ZTuqeOo/sbcHS5MFtNJGSFkpwdQVJ2OGGxATLaFoNCgno4KtsKnz4Fh1YZn4+9DGbdBWkLTrnxCFDUVMRrBa+xqmgV9Z312HxsLEtZxooxK5geM/2UqZGuI0doev0NWt59F2dNDabAQGzLlhFyxeUEzJ49KrZd7QuXw01FYRMlu43QbqpuByAo3Jfk8eEkT4ggcVwYvgGjt8tGDCwJ6uGs6Sjk/cXoFuloMPbCnnUnTLr+C9MiTreTLZVbeK/oPdYdXUeHs4PYwFiWpy1nRdoKxoaP7X2sdrlo37qV5nfepXXtWtxtbViiogi+7DKCr7wCv+xsGTmepKWug6P7Gyjb38Cxgw10d7pQJkVMajDJE8JJyg4nOiUYk0mumTg//b2ZuBC4FbgU+F/gV1rrjjN9jwT1IHB0GItntv4JKneBbwhMvRlm3gER6V94eLujnQ1lG3iv+D0+K/8Mp3aSEZrBijErWJG2griguN7Hujs7advwMc3vvEPbxo3gcOAzZgwhV1xO8OWXj/pWv89zudxUF7dQtr+Bo/vqqTnaChp8Ay0kjTNCOzk7XE6xEedE+qhHEq2NaZGtzxg3IN1OSJ1v3HgcfwVYv9iK19DZwNqStbxb9C67ancBkBOZw+LkxSxOXkxqSGrvY11NTbSsWUvLO+/Q3vPf0G/SJIIvvZTgSy/BGh8/JG9zOOlo6+bYgUaO7q/n6P4G2pu7AQiNCSBxXBiJY8NIGBuGX6BMk4gvJ0E9UrVWQf4LsOMFo73PNwRyroWpt0D81C/MZQOUtZaxpmQN60rXsa9+HwAZoRksTl7MkpQljA0b2zvl4Sgvp3nVKlrfX03n/v0A+E+ZQvDyS7FdeinWmJihe6/DhNaahgo7R/c3cOxgIxWFTTi7XKAgMjGIxHHhJI4LIy49BB8/i6fLFV5Egnqkc7uNLpEdL8L+f4Gz09gIauotkHMdBJ6+Fa+yrZL1ZetZV7qO/Jp83NpNQlACS5KXsDhlMZOjJvfeiOwuLaXl/dW0rF5N18GDAPhPn07w8uXYli3FGh09ZG93OHG53NQUt3DsUCPHDjZSVdyM26kxmRQxY4JJGBtG0rgwYtJCMFtkpeRoJkE9mnQ0wd7XjVF2xQ4w+xgdI5NvMJasW06/Aq+hs4ENZRtYV7qOzZWbcbgdRPpHsihpEQuTFjIzdib+FmNapauomJbV79P6/mq6CgpAKQJmzMC2/FJsS5ZIaJ+Bo9tFZWET5T3BfXx+2+JjIi4jlMSxYcRnhRKVbMNsluAeTSSoR6uqvcbKx12vGB0jfqGQvdKYHkm5AEynb8Vr624jtzyXdaXryC3PpcPZga/ZlxmxM5ifMJ8FCQtICjZuMHYVFNCyeg0t779Pd1ERKIX/tGkEX7IM29KlWOPiTvsawtBpd1BR0NQ74m6stANg8TUTlx5CfGYoCZmhRKcGy4h7hJOgHu2c3VD0Eex5DQ6+Bw47BMXCxKth4rWQMO2089lgLKzZXr2d3GO5fFL+CSUtJQCkBqdyYcKFzE+cz4yYGVhNVroLC2lZs5bWtWvpOnwYAL/JkwhetgzbsmXSPdIH9uYuKgqaej8aKnqC22oiZkwICVmhxGeGEpMWjMUqPe8jiQS1OKG7HQ6vNqZHCtaCqxvC0mDiNcZIO3r8Gb/9aMtRcstzyS3PJa8qjy5XF/4Wf2bHzmZ+4nwuTLiQ+KB4uoqLaV37Aa1r1vTeiPTNHt8T2pfgOyZtKN7tsNfR2k1FYRMVh5soL2gy9tzWYLaYiEkLNkbcWaHEjAmRpe7DnAS1OL2OJjj4rjHSLv4YtBsis2DcChh3OcRPO2U3vy98u7ODbVXbyD1mBHd5WzlgjLZnx81mdtxsZsXOwr+2tTe0O3YZ7YG+mRnYli7FtmQJvuPHy+KaPuq0O4w57gIjvOvKWtEaTGZFVLKN+IxQ4jJCiEsPxS9I2gGHEwlqcXZtNUbHyMF3oeQToz/bFmfciBy3wujV/pIbkWC0pRW3FPPJsU/YUrWFvKo82p3tKBTjI8YzO242c+LmkOOOx7E+l9Z164w+bbcba0KCEdrLluI/ZQrqDP84iFN1dTipLGyisrCJioJmakpbcLuM/6fD4gKJywgxwjs9BFuEn/yD6MUkqMW56WiEw2uN0C5cB452o0c7a5kR2hlLwNd2xqdwuB3srdvL5srNbKncwq7aXTjdTqwmK1OipzAnbg6z/caTuKsS+7oPsX+2CRwOzFGR2BYvxrZ0KYGzZo3qE2vOh7PbRU1pCxWFzVQWNlN1pInuThcAQWG+xKWHEJcRSlxGKBHxgShZ8u41JKjF+XN0QNEGI7QPvQ/t9UbLX9oCyFhqhHZE+pfejDyu3dFOfk0+Wyq3sLlyMwcbjF7sIGsQ02OmMytoIlOLNCGbDtCe+wm6owNTcDC2iy8yQvuCC0b9AQjnw+3WNFS0UVHQTOWRJioLmrD3rJz08bcQOybECO/0EKJTg7H6yjy3p0hQi4HhdkHZFjjwrnEjsr7A+PPQFCOwM5YYAe579tPRGzsb2Vq1lc2Vm8mryuvtJvG3+DM9ZCKLqiIYu6cJ3027cbe0oPz9CZo/H9vSJQQtXIg5OPjMLyBOS2tNa32nMVVS2EzlkebelkCTSRGZFERsekhPgIcSFObr4YpHDwlqMTgaiuHIh1D4IRR9bLT9mazGob3Hgzs6+6yjbYC6jjryq/PZXr2d7dXbOdx4GI3GT1tY3pTChYUWEnYcw1TfDFYrgXPmYFuyBNviRVgiI4fgzY5cnXYHVUXNVB0xgrumpAWnww2ALdyP2J4Rd2x6CBEJQbJD4CCRoBaDz9kNZZuh4AMjuGuMfUSwxUH6IuNmZOqFENq3XurmrmZ21OzoDe799ftxu52MrTBxaWkokw90EljTaiywmT6N4KVLjQU2smlUv7lcburK2nqCu4nKI829G01Z/czEpgUTO8YI7pi0EHz9Zc+SgSBBLYZeS4UR2IUfQPFG4wYlQFiqEdppC4xfg/u2crHd0c7O2p29wb23dg8xVV3MPuTmggIL8dUO44Hj0olYfjkhyy7BN016tQdC73TJkROj7oaKNrQGFETEB/YGd+yYEEKi/KW75DxIUAvPcruhZj+U5EJxrrGBVGez8bWIDGOknTrf+LD1bUc+h8vBwYaD7Kzdyc6anZQfyCNjVx2zDrnJrDQe05YUjumieaRceT0RE6dLeAyg7g4n1SUtvVMmVUXNvd0l/jarEdw94R2dbMMii3HOSoJaeBe3C6r2GP3aJblQ+hl0tRhfi8yC5DmQNNv4iMjo0xy31poqexU7a3dy8MAnuD/eRFp+FWPLNCagLsJKzcw0fJdeTMacS8kIy8Bikh/ZB4rbrWmstJ+Y6y5qprnGOF/k+GKc4x0msekhBIbITcrPk6AW3s3lhKpdPcH9qdFZ0tlkfM0/vCe0ZxkBHj/1tIcjnE67o519hz+lcs3b+GzcTuLhRsxuqA2G7eOt1M3JJGLGPCZFTyYnMoeYQNlfeyC1t3RTXWxMlVQVNVNT0orLadykDI706x11x2WEEB4vNyklqMXw4nYbrX9lW4yPo1tOtAKaLBA3+cSIO2kWBPftBqKzsZGj779Bw+r38Ms/hNnppiEItoxVbBmrqB8bw6SYKeRE5pATmUN2RDYB1oBBfKOji8vppvZo6ykdJu0tp79JGZsWgs8ou0kpQS2GP3s9HNt2IrzLtxsHJICxE2D8VGMXwPipxh4lX3JYwnGutjbaPtpA85r3aduYi+p20GHzIX+8L+vS7RxIVmA2kx6azsSIiUyMnMiEyAlkhWZhNctqyYFwupuU9RXGplNKQURiEHEZob37l4z06RIJajHyOLuheg+UbTMOSKjIh7oCoOfvc2jyidBOmGaMwv1CTvtUbrudttxcWtasoW3Dx+iODtyhNiqnJ7N1nJm14RXUO42pGB+TD+PCxzEhcgI5kTlMiJxAanBq70k4on+6O5xU9UyXVBY2U13cjLP7xHTJycEdGhMwom4QS1CL0aGzxTihvSLfCO/yfOMsyeMiMo3wjptkHFUWO+kLI293Rwdtubm0rl5D24YNuNvbMYeGYlo4h2Mzk8lPdLK7aR/76/fT4TRulgVZg8iOyGZC5AQmRhgj7/jA+BEVIp7icrmpO9pm9HMXGn3dHa1GK6ZfkPXE4QpZYUQkDu95bglqMXrZ66FyB5T3jLordkJrxYmv2+IhNufUj7A0MJlwd3Zi//RTWlavoW39etx2O6aQEGyLFxO4dAk1E2LZ23yQffX72Fu3l0ONh3C6nQCE+YaRHZnNhIgT4R0dIEeU9ZfWmuaaDiqO7xhY2ExLrfEPpo+/pXd/7vjMUCKTbMMquCWohTiZvd6YNqk66aP2EGijDxifIIiZYIR2zESImYA7OA379r20rllN64frcbe1YbLZsC1ahO2SSwi8YB5Oi6KgsYC9dXvZV7+PffX7ONJ0BFfP80b7R58I78iJTIiYQJhfmAcvxMjQ1thJ+eEmKg43Ul7Q1NsW6ONvIS4jhITMMBLHhRGZGOTVuwVKUAtxNo5OqD14anhX7z3R3w0QkgTR43GHZmKv9qV1VyWtn+Xjbm3FFBhI0EUXYVu2jKAF83t3+utwdnCo4dAp4V3SXILumUtPCEowpk16wnt8xHiCfWTDqf6wN3VRXtDYE95NNFW3A8ZUSdK4MBLHh5M0PhxbuJ+HKz2VBLUQ50NrY4675oCxsrLmoPH7ukPGEWaAdinsHcm0VgTTeqgNl70L5etD0IL52C5ZTtDCBZhtp+7d3dbdxoGGA+yr28fe+r3sq9vHsbZjvV9PCU4hOyK7d8pkfPh4aRPsB3tTF8cONlB2oJGyAw29LYGhMQEkjQ8nKTucxLFhHt/iVYJaiIHkckJDkRHetQd7QvwAuvYI7TVmWsv8aTnmh6vTDCZF4NhobPOmYlt6KZasmRAQ/oWnbOpsYn/9/t5R9966vVS3VwNgUibGhIxhYuTE3k4TaRM8P1prGirslB1ooOxAAxWHm3A63JitJpLGhZE6KZLUnEgCQ4e+FVCCWoih4OyCusNQdxhdfYiOHfm05h+h9bAdR5sZ0PhHOLCNMRE0OQWfrGxU1FijGyUi3djX23xikUddR13vqHtvnTHybuwyNreymqxGm2DEBHKicpgYOVHaBM+Dy+Gm4kgTJbvrKNldR0ud0ZsfnWIjdVIkGdOjCYsNHJJaJKiF8CDtctKVn0vr++/Q9mkenaW1AFhtbmxx7QQldBIQ2Y2yWoyOk4gMI7hP/tUWhwbK28p7p0v21u1lf/1+2p3GHKzNx8bEiInkROUwKXISOVE5hPt9cfQuTu/4aLtkTx3Fu+qoLmkBDVHJNrJmxZAxPWZQD1KQoBbCizgqKmj7+GNa139E+5Yt6O5uTP6+BI6LJijJRGBkE9bu4hMrLwGsgRAxBsJP+ghLwxWWQrG7kz31+9hdt5s9tXsobCrs7TRJCEpgUtQkJkdNZkrUFLLCs7CaZMqkL+xNXRTkVVOwrZqa0lZQkDg2jJyLEkmdFDngrX8S1EJ4KXd7O/ZNm2j96CPsG3Nx1tQA4JuVReDMyQRlx+EfA6bWUqgvhMZiaCwFt+PEk1j8jJF4eBqEj6E9JIEDVjN7XG3stpexq24PNe3G8/qZ/ZgQOYHJUZOZGj2VqdFTCfE9/YpNcUJjlZ3D26o5+FklbY1d2ML9mHhRAhMXJODjNzB7kkhQCzEMaK3pOlyA/ZNc2jbm0p6fDw4HytcX/2lTCZw7j8C5c/DLykS1Vxk3NBuKT/218XMjcWWC4ESqwhLYGRjMLgvsdDZzsKMaZ8+oOyM0g6nRU5kWM43p0dOJC+rbYQ6jkdvlpnh3HXs+Okb54Sb8bVZmrkgje348ZnP/7g/0O6iVUkuAq4EaQGutf3qmx0tQC9F/rjY77XnbaN+0GfumTXQdPgyAKSiIgOnTCZg1i4BZM/EbPx5l6RnVud3QWmm0FTYUQ2PJqR92Y2TdqRR7fX3YEWBje6CNXRZNG8aeGom+4cyMnMTMxAuZmbSQ2MDYIX/vw0FVcTOb3jhCRUET0Sk2lt0xgZCo82+j7FdQK6UCgN3ABK11l1LqdeAPWusPv+x7JKiFGHjOujrsm7fQvm0b7Vu30l1cDIAKCMB/8iQCps/Af+oU/CdPxhz0JSfBd7UZId5YYkyhNJVC01FcjSUU2MvZbnaxzc+PPD9fms1GX3Gq28QcSxhzbWnMjp5CYFi6selVaBL4hfbpYIeRSmvNkfxaNrx0EJNZcdV90867S6S/Qb0Y+JHWenHP598DErXW3/uy75GgFmLwOWpq6Ni+nfa87bRv307XoUPGIh2l8M3MJOl//4A1IaHvT6g1tDdAUwnuxhIKanaxuWEfW9rLyXPb6VBg0ZoF7R08VVNnfI+PzQjsr73V52PURqKm6nbeeHI7wZH+XPP/zu/YtzMFdV9mwaOB1pM+b+n5s8+/yF3AXQDJycnnXCRAdrwsnRWir6zR0ViXLyd4+XLA2GO7Y9cuOnbspGPPbixRUef2hEoZuwkGRmBKmM5YrmEs8HWMMyp31uwgt2QddDXDRfOg6Sg0lUFzGfiP7j1LQmMCmHtVBjWlLbicbizWgV3l6FUjaiGEGK3ONKLuy23KTUCKUup4p/cFwHsDVZwQQogzO+vUh9a6XSl1N/C0UqoW2H2mG4lCCCEGVp86tbXWHwAfDHItQgghTkN2cBFCCC8nQS2EEF5OgloIIbycBLUQQng5CWohhPByg7J7Xk8bX+l5fnskUDeA5QxXch0Mch0Mch0MI/k6pGitT7ucdFCCuj+UUnlftjpnNJHrYJDrYJDrYBit10GmPoQQwstJUAshhJfzxqB+xtMFeAm5Dga5Dga5DoZReR28bo5aCCHEqbxxRC2EEOIkEtRCCOHlBuac8z442wG5Sik/4EmgHMgEntBaH+752i3AVMAFHNFa/2mo6h5o53sdlFIzge8CO4CxwFat9Z+HtPgB1p+/Ez1fj8a4Ho9rrX83ZIUPsH7+vzEHWAq4gYuB27TWZUNY/oDp53X4H8CBMfgMAL6ttXYPYfmDS2s96B8YF64Q8O35/HVg8ececz/w/3p+nwPk9vw+EdjJifn0bUDmUNTtZdfhSmBWz++tQCMQ6en35Ilr0fO5Cfgz8C/gXk+/Hw/9nQgGXj/pcWOAQE+/Jw9ch9nArpMetwu4wNPvaSA/hmrqYy5QqrXu6vn8U2DF5x6zAuM0GbTWe4DJSqlg4BJgu+75L9DzmOWDX/KgOO/roLV+W2u99aTHOTFGEMNVf/5OAPwQeBbjH6zhrD/X4TKgTSn1PaXUQ8A0rbV9iOoeaP25DvVAkFLKopSyABooHpqyh8ZQTX305YDcL3tMnw7XHSb6cx1aTvqze4HHtNbNg1HkEDnva6GUmg60a6239Jw+NJz15+9ECsZo8g6MacGPlFL1WuuPBq/cQXPe10FrXaiUegb4J8YU0DqgdhBrHXJDNaKuAWwnfR7c82d9eUxfvne46M91AEApdRPGj7e/Hqwih0h/rsVKwF8pdT/Gj8BLlVK3DWKtg6k/16EF2KG1dmhjPnYTsHAQax1M530dlFJXAhdrra/SWl8DpAF3DmaxQ22ogvq0B+QqpcJP+lH2PYwff1BK5WDMObUAa4DpSinV87i5wPtDVPdA6891QCl1B8YI4mdKqRylVNYQ1z+QzvtaaK2/q7V+Qmv9BLAH+EBr/dehfgMDpD9/Jz4CUk96rhTgMMNTf65DElB10nNVAn5DU/bQGLIFL0qppcC1GD+SOLTWP1VK/QJo0Fo/oZTyx7ijWwlkYPxof3LXxwyMH+8O6+Hd9XFe10EptRJ4HqPLASAC4872hiF/EwOkP38ner7/3zCmgcqBP2ith+U/4P38f+NbGGHtAPyB7+th2u3Qj/83AoE/YezY6cIYUf/7MJ6v/wJZmSiEEF5OFrwIIYSXk6AWQggvJ0EthBBeToJaCCG8nAS1EEJ4OQlqIYTwchLUQgjh5f4/OHGP1VV77eEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(c_array[0,:],z)\n",
    "plt.plot(c_array[40,:],z)\n",
    "plt.plot(c_array[45,:],z)\n",
    "plt.plot(c_array[50,:],z)\n",
    "plt.plot(c_array[-1,:],z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}