fdtd1d.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to run an electromagnetic simulation, with initial code and ideas borrowed from\n",
    " * Understanding the Finite-Difference Time-Domain Method, John B. Schneider, http://www.eecs.wsu.edu/~schneidj/ufdtd, 2010.\n",
    " \n",
    "We're going to run a simple, bare-bones 1D FDTD simulation with a hard source.\n",
    "\n",
    "The impedance of free space (or vacuum) is 377.0.\n",
    "\n",
    "We are going to model 400 mm of space, and run the simulation for 500 time units.\n",
    "\n",
    "We use a time step that matches the space step with a Courant number of 1, which means that $c * dt / dx = 1$, where c is the speed of light.  Given that $c$ is 299,792,458 m/s and our space step ($dx$) is 1 mm, our time step ($dt$) is $dx / c$ or $3.33*10^{-12}$ s. (Don't worry about this if it doesn't make sense - read more in an FDTD book if you want to.)\n",
    "\n",
    "Our source is the electric field at the left edge of the grid.\n",
    "We use a Gaussian source that peaks at 30 time units and has a standard deviation of 7.\n",
    "\n",
    "Our output will be the electric field measured at 250 mm over the time of the simulation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "from IPython.display import HTML\n",
    "\n",
    "\n",
    "imp0 = 377.0 # property of free space (vacuum)\n",
    "SIZE = 400 # dimension of space to model\n",
    "sensorLocation = 250 # location of output sensor\n",
    "maxTime = 400 # simulation time\n",
    "\n",
    "sourcePeakTime = 30 # peak of the Gaussian source\n",
    "sourceSdv = 7 # standard deviation of the Gaussian source\n",
    "sourceSigma = 2 * sourceSdv**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first initialize lists for the electric and magnetic field components.\n",
    "\n",
    "We then start at time 0, and update the magnetic fields at time 0.5 based on the magnetic fields at time -0.5 and the electric fields at time 0.  And then we update the electric fields at time 1 based on the electric field at time 0 and the magnetic field at time 0.5. At the left edge of the grid, we set the source electric field.  And we write out the electric field at the sensor location.\n",
    "\n",
    "We then repeat this for each time step, until the set number of time steps have been reached."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.00010270256461965522\n",
      "0.00018751857457483846\n",
      "0.0003354626279025118\n",
      "0.0005880047627378895\n",
      "0.0010098442341126815\n",
      "0.001699279365552657\n",
      "0.0028016380348484843\n",
      "0.004525807866766719\n",
      "0.007163364470919217\n",
      "0.011108996538242306\n",
      "0.01687988414878991\n",
      "0.025130488985274417\n",
      "0.036658041953378004\n",
      "0.05239314106982558\n",
      "0.0733696513683829\n",
      "0.10066889977289778\n",
      "0.1353352832366127\n",
      "0.17826397958504792\n",
      "0.23006629899380904\n",
      "0.2909238070461482\n",
      "0.360447788597821\n",
      "0.43756473770026916\n",
      "0.5204501210207021\n",
      "0.6065306597126334\n",
      "0.6925693242051977\n",
      "0.7748374288832494\n",
      "0.8493658165683124\n",
      "0.9122540768285452\n",
      "0.9600054412854777\n",
      "0.9898478033787339\n",
      "1.0\n",
      "0.9898478033787338\n",
      "0.9600054412854776\n",
      "0.9122540768285451\n",
      "0.8493658165683124\n",
      "0.7748374288832494\n",
      "0.6925693242051977\n",
      "0.6065306597126334\n",
      "0.5204501210207021\n",
      "0.43756473770026916\n",
      "0.360447788597821\n",
      "0.2909238070461482\n",
      "0.23006629899380907\n",
      "0.17826397958504792\n",
      "0.1353352832366127\n",
      "0.10066889977289775\n",
      "0.0733696513683829\n",
      "0.05239314106982558\n",
      "0.036658041953378004\n",
      "0.025130488985274414\n",
      "0.01687988414878991\n",
      "0.011108996538242306\n",
      "0.007163364470919217\n",
      "0.004525807866766719\n",
      "0.0028016380348484843\n",
      "0.001699279365552657\n",
      "0.0010098442341126815\n",
      "0.0005880047627378895\n",
      "0.00033546262790251185\n",
      "0.00018751857457483843\n",
      "0.00010270256461965522\n",
      "5.5113137010524656e-05\n",
      "2.89778274286711e-05\n",
      "1.4928403298804454e-05\n",
      "7.535251360782923e-06\n",
      "3.726653172078671e-06\n",
      "1.805830798460929e-06\n",
      "8.573773029336787e-07\n",
      "3.98844636715367e-07\n",
      "1.8179099965707518e-07\n",
      "8.118538349144053e-08\n",
      "3.552386103252793e-08\n",
      "1.522997974471263e-08\n",
      "6.397574269770496e-09\n",
      "2.6331050999142502e-09\n",
      "1.0618371101937938e-09\n",
      "4.1955070187358226e-10\n",
      "1.624231211810119e-10\n",
      "6.160955882896819e-11\n",
      "2.289734845645553e-11\n",
      "8.337947094201161e-12\n",
      "2.9748832347325257e-12\n",
      "1.0399622154910407e-12\n",
      "3.562066673900191e-13\n",
      "1.195427865933476e-13\n",
      "3.930805485965858e-14\n",
      "1.2664165549094176e-14\n",
      "3.9976839079982126e-15\n",
      "1.2364517461413708e-15\n",
      "3.7469917746990634e-16\n",
      "1.1125643881294095e-16\n",
      "3.2367146737857757e-17\n",
      "9.226150318955576e-18\n",
      "2.5767571091549857e-18\n",
      "7.051204120964661e-19\n",
      "1.8905577075969974e-19\n",
      "4.966534275776451e-20\n",
      "1.2783617954222807e-20\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n",
      "2.1666711874356403e-33\n"
     ]
    }
   ],
   "source": [
    "ez = [0.0] * SIZE\n",
    "hy = [0.0] * SIZE\n",
    "\n",
    "# do time stepping\n",
    "for qTime in range(maxTime):\n",
    "\n",
    "    # update magnetic field\n",
    "    for mm in range(SIZE-1):\n",
    "        hy[mm] = hy[mm] + (ez[mm + 1] - ez[mm]) / imp0\n",
    "\n",
    "    # update electric field\n",
    "    for mm in range(SIZE):\n",
    "        ez[mm] = ez[mm] + (hy[mm] - hy[mm - 1]) * imp0\n",
    "\n",
    "    # hardwire a source node */\n",
    "    if qTime < sourceSigma:\n",
    "        ez[0] = math.exp(-(qTime - sourcePeakTime)**2 / sourceSigma)\n",
    "    else:\n",
    "        ez[0] = 0.0\n",
    "\n",
    "    print(ez[sensorLocation])\n",
    "#done with time stepping loop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of printing the output electric field values at the sensor, let's save them to a list and then plot them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "ez = [0.0] * SIZE\n",
    "hy = [0.0] * SIZE\n",
    "output = [0.0] * maxTime\n",
    "\n",
    "epsR = [1.0] * SIZE\n",
    "\n",
    "# do time stepping\n",
    "for qTime in range(maxTime):\n",
    "\n",
    "    # update magnetic field\n",
    "    for mm in range(SIZE-1):\n",
    "        hy[mm] = hy[mm] + (ez[mm + 1] - ez[mm]) / imp0\n",
    "\n",
    "    # update electric field\n",
    "    for mm in range(SIZE):\n",
    "        ez[mm] = ez[mm] + (hy[mm] - hy[mm - 1]) * imp0 / epsR[mm]\n",
    "\n",
    "    # hardwire a source node */\n",
    "    if qTime < sourceSigma:\n",
    "        ez[0] = math.exp(-(qTime - sourcePeakTime)**2 / sourceSigma)\n",
    "    else:\n",
    "        ez[0] = 0.0\n",
    "\n",
    "    output[qTime] = ez[sensorLocation]\n",
    "#done with time stepping loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGYZJREFUeJzt3X+QHOV95/H3Vz8RIBCgBWMkENiiAoUxEB2mDirxBd8ZOBeKK3AlLnd2rjhThYPvLrZjQ9nmfLiSiu3K+ZIqHExyHI7vYsDYvuhccnAqxomTChhhfgqMWRA2awisAP0GhKTv/dG9w2g1szuzOz2zM/N+VW1td8+z3V+1dvezz/NMd0dmIkkSwLxeFyBJmjsMBUlSjaEgSaoxFCRJNYaCJKnGUJAk1RgKkqQaQ0GSVGMoSJJqFvS6gHYtX748V61a1esyJKmv3H///Vsyc2S6dn0XCqtWrWLjxo29LkOS+kpE/KyVdg4fSZJqDAVJUo2hIEmqMRQkSTWGgiSpprJQiIhbIuLFiHi0yesREX8cEaMR8XBEnFNVLZKk1lTZU7gVuGiK1y8GVpcfVwF/UmEtkqQWVBYKmfl3wMtTNFkL/HkW7gGWRcTxVdUjqfNeegn27Ol1FeqkXs4pnAA8W7c+Vm47SERcFREbI2Lj+Ph4V4qTNL3f+A343Od6XYU6qZehEA22ZaOGmXlzZq7JzDUjI9NepS2pC556CjZuhFtvhX37el2NOqWXoTAGrKxbXwE816NaJLXpr/4K1q2DQw+FTZt6XY06pZehsB74QPkupPOAbZn5fA/rkdSGp5+GU0+F1ath8+ZeV6NOqeyGeBHxdeDdwPKIGAP+K7AQIDNvAjYAlwCjwG7gP1RVi6TO27wZzjsPnn22CAgNhspCITOvmOb1BH67quNLqtbmzXDKKTA2Zk9hkHhFs6QZ2bwZTj65+DAUBoehIKltr7wCmXDUUYbCoDEUJLVtbAxWroSI4vPYWK8rUqcYCpLaNj4OE5cMLVsGu3Z5ZfOgMBQktW18HJYvL5bnzYNjjilueaH+ZyhIatuWLW/2FKAICO9AMxgMBUltqx8+gmLZUBgMhoKkttUPH0ERClu29K4edY6hIKltk3sKDh8NDkNBUtu2bLGnMKgMBUltc05hcBkKkto2uafg8NHgMBQktW3r1uIWFxOOOqrYpv5nKEhqy2uvFU9aW7LkzW3LlsG2bb2rSZ1jKEhqy7ZtRQhE3QN1jzzSnsKgMBQktWXr1iIU6i1bZigMCkNBUlu2bSt6BvUMhcFhKEhqS6OewiGHFMNJr73Wm5rUOYaCpLY0CgVwXmFQGAqS2jIx0TyZQ0iDwVCQ1JatWw+eUwBDYVAYCpLaMtXwkdcq9D9DQVJbHD4abIaCpLY06yksWwavvNL9etRZhoKktmzbBkcccfD2I46AHTu6X486y1CQ1JYdO2Dp0oO3L11qKAwCQ0FSW3bsaNxTOPxw2Lmz+/WoswwFSW3Zvt2ewiAzFCS1pdnw0eGHGwqDwFCQ1Jap5hQcPup/lYZCRFwUEU9ExGhEXNvg9RMj4u6IeCAiHo6IS6qsR9Ls7NtX3PTusMMOfs3ho8FQWShExHzgRuBi4HTgiog4fVKzTwN3ZObZwDrgy1XVI2n2duwohonqH7AzwYnmwVBlT+FcYDQzn87MPcBtwNpJbRKYeB/DkcBzFdYjaZaaDR2BPYVBsaDCfZ8APFu3Pga8a1KbzwLfi4iPAIcB76mwHkmzNFUo2FMYDFX2FBp0MMlJ61cAt2bmCuAS4GsRcVBNEXFVRGyMiI3j4+MVlCqpFc3ejgr2FAZFlaEwBqysW1/BwcNDVwJ3AGTmPwKHAMsn7ygzb87MNZm5ZmRkpKJyJU2n2YVrUEw+79oF+/d3tyZ1VpWhcB+wOiJOjohFFBPJ6ye1+TlwIUBEnEYRCnYFpDlqquGj+fNhyRLYvbu7NamzKguFzNwLXAPcBTxO8S6jTRFxQ0RcWjb7GPChiHgI+DrwW5k5eYhJ0hwx1fAROIQ0CKqcaCYzNwAbJm27vm75MeD8KmuQ1DlT9RTAC9gGgVc0S2rZzp1Th4K3uuh/hoKklu3a1fhq5gn2FPqfoSCpZTt3Fr2BZuwp9D9DQVLLpgsFJ5r7n6EgqWUOHw0+Q0FSyxw+GnyGgqSW7do1/fCRPYX+ZihIatnOnVMPH9lT6H+GgqSWOdE8+AwFSS1z+GjwGQqSWubw0eAzFCS1rJXhI3sK/c1QkNSSN96Afftg8eLmbZxT6H+GgqSWTFy4Fo2eqVjykZz9z1CQ1JLpho7AnsIgMBQktWS6dx6BE82DwFCQ1JLp3nkETjQPAkNBUkta6SksXAjz5sGePd2pSZ1nKEhqSStzCuBkc78zFCS1pJXhIzAU+p2hIKklrQwfgaHQ7wwFSS1ptadw2GGGQj8zFCS1xDmF4WAoSGqJw0fDwVCQ1BInmoeDoSCpJfYUhoOhIKklzikMB0NBUkvaGT7atav6elQNQ0FSSxw+Gg6GgqSWOHw0HAwFSS3x4rXhUGkoRMRFEfFERIxGxLVN2vybiHgsIjZFxF9UWY+kmXP4aDgsqGrHETEfuBH4l8AYcF9ErM/Mx+rarAauA87PzFci4tiq6pE0Ow4fDYcqewrnAqOZ+XRm7gFuA9ZOavMh4MbMfAUgM1+ssB5Js+DFa8OhylA4AXi2bn2s3FbvVODUiPiHiLgnIi5qtKOIuCoiNkbExvHx8YrKldTMG2/Avn2wePH0bQ2F/lZlKESDbTlpfQGwGng3cAXwZxGx7KAvyrw5M9dk5pqRkZGOFyppahPzCdHop3oSQ6G/VRkKY8DKuvUVwHMN2vxlZr6RmZuBJyhCQtIcsmtXa0NHYCj0uypD4T5gdUScHBGLgHXA+klt/i/wLwAiYjnFcNLTFdYkaQZanWQGr2jud5WFQmbuBa4B7gIeB+7IzE0RcUNEXFo2uwt4KSIeA+4GfjczX6qqJkkz004oLFkCr74K+/dXW5OqUdlbUgEycwOwYdK26+uWE/ho+SFpjmpn+GjePDj0UNi9u/Ug0dzhFc2SptVOTwG8qrmfGQqSptVuKDjZ3L8MBUnTamf4CAyFfmYoSJqWPYXhYShImpahMDwMBUnTcvhoeBgKkqY1k56CF7D1J0NB0rTsKQwPQ0HStJxTGB6GgqRpefHa8DAUJE2r1QfsTLCn0L9aCoWI2BcRfxDx5t3UI+LH1ZUlaS7ZuROWLm29vaHQv1rtKWwq234vIo4ut7XwuA1Jg8A5heHRaijszcxPAH8K/DAifpmDn6ImaUDt2GFPYVi0euvsAMjMOyJiE/B14MTKqpI0p9hTGB6thsJ/nFgoH5RzAfDr1ZQkaS7JNBSGyZTDRxHxCYDMvD8iLp/YnpnbgV+quDZJc8Brr8HChbCgjUdyeUVz/5puTmFd3fJ1k167qMO1SJqD2u0lgD2FfjZdKEST5UbrkgbQjh2GwjCZLhSyyXKjdUkDqN1rFMArmvvZdKOE74yI7RS9giXlMuX6IZVWJmlOmMnw0aJFxec9e95cVn+YMhQyc363CpE0N81k+AjeHEI6+ujp22ru8N5HkqY0k+EjcF6hXxkKkqY0k+EjMBT6laEgaUqzHT5SfzEUJE3J4aPhYihImtJsho+8qrn/GAqSpuTw0XAxFCRNyeGj4WIoSJrSTHsKXtXcnyoNhYi4KCKeiIjRiLh2inaXRURGxJoq65HUPt+SOlwqC4WImA/cCFwMnA5cERGnN2i3FPhPwL1V1SJp5hw+Gi5V9hTOBUYz8+nM3APcBqxt0O5zwBeA1yqsRdIMOdE8XKoMhROAZ+vWx8ptNRFxNrAyM79TYR2SZsHho+FSZSg0et5C7XbbETEP+BLwsWl3FHFVRGyMiI3j4+MdLFHSdBw+Gi5VhsIYsLJufQXwXN36UuAM4AcR8QxwHrC+0WRzZt6cmWsyc83IyEiFJUuabDbDR1681n+qDIX7gNURcXJELKJ4tOf6iRczc1tmLs/MVZm5CrgHuDQzN1ZYk6Q27N8Pu3cXby9tlz2F/lRZKGTmXuAa4C7gceCOzNwUETdExKVVHVdS5+zeDYceCvNm8JvCUOhP0z15bVYycwOwYdK265u0fXeVtUhq30yHjsCL1/qVVzRLamqm7zwCewr9ylCQ1NSOHTN75xEYCv3KUJDU1Gx6CoceCq++WkxWq38YCpKa2rYNjjhiZl87b14xr7BjR2drUrUMBUlNbdsGy5bN/OuPPLLYh/qHoSCpqW3bil/sM7VsmaHQbwwFSU3NNhTsKfQfQ0FSU50Iha1bO1ePqmcoSGpq61Z7CsPGUJDUlMNHw8dQkNSU7z4aPoaCpKbsKQwfQ0FSU040Dx9DQVJTXqcwfAwFSU05fDR8DAVJDe3dWzxkZ6Y3xANDoR8ZCpIa2r69uBneTJ66NsFQ6D+GgqSGZjt0BIZCPzIUJDXUqVDw3Uf9xVCQ1FAnQmHp0mJeYt++ztSk6hkKkhrqRCjMm1dMVG/f3pmaVD1DQVJDnQgF8FqFfmMoSGpo69bZ3fdogpPN/cVQkNRQp3oKhkJ/MRQkNdTJUPAdSP3DUJDUkD2F4WQoSGrIUBhOhoKkhgyF4WQoSGqoU+8+WrbMOYV+YihIauill+CYY2a/n2OOKfal/mAoSGrIUBhOlYZCRFwUEU9ExGhEXNvg9Y9GxGMR8XBE/E1EnFRlPZJas29fcWuKTgwfGQr9pbJQiIj5wI3AxcDpwBURcfqkZg8AazLzTOBO4AtV1SOpda+8UkwQz58/+30ZCv2lyp7CucBoZj6dmXuA24C19Q0y8+7M3F2u3gOsqLAeSS3asqUzQ0dgKPSbKkPhBODZuvWxclszVwLfbfRCRFwVERsjYuP4+HgHS5TUSKfmEwCOPrroeWR2Zn+qVpWhEA22Nfy2iIh/B6wBvtjo9cy8OTPXZOaakZGRDpYoqZFOhsKiRbBkidcq9IsqQ2EMWFm3vgJ4bnKjiHgP8Cng0sx8vcJ6JLWok6EAsHy5Q0j9ospQuA9YHREnR8QiYB2wvr5BRJwNfIUiEF6ssBZJbeh0KDiv0D8qC4XM3AtcA9wFPA7ckZmbIuKGiLi0bPZF4HDgGxHxYESsb7I7SV1URShs2dK5/ak6C6rceWZuADZM2nZ93fJ7qjy+pJl54QU45ZTO7e/YY+FFxwL6glc0SzrICy/Accd1bn/HHVfsU3OfoSDpIIbC8DIUJB3EUBhehoKkA2QW4/+GwnAyFCQdYOtWWLy4uOCsUwyF/mEoSDpAp4eOwFDoJ4aCpAO88AK85S2d3efy5cX9j/bu7ex+1XmGgqQDPP9850NhwYLiAjavVZj7DAVJBxgbg5Urp2/XrhUrin1rbjMUJB3g2WeLX+CdtmJFsW/NbYaCpANU1VNYudKeQj8wFCQdoKqewsqV9hT6gaEg6QBjY9UNH9lTmPsMBUk1b7xR3OL6+OM7v297Cv3BUJBU8/OfF4GwoIKb6q9aBc880/n9qrMMBUk1o6OwenU1+z7hBHj5Zdi1q5r9qzMMBUk1o6Pw9rdXs+9584oH9zz1VDX7V2cYCpJqqgwFKPY9Olrd/jV7hoKkmiefNBSGnaEgqeanP61uTgHg1FPhJz+pbv+aPUNBElBMAI+NFb+4q/KOd8DDD1e3f82eoSAJgEcegdNOg4ULqzvGO94Bjz3mLbTnMkNBEgAPPQTvfGe1x1i6FN761mKYSnOToSAJgAcegLPOqv44Z50FP/5x9cfRzBgKkgD44Q/hgguqP84FF8Df/331x9HMGAqSGB+HX/yi+uEjgF/9Vfjbv63+OJoZQ0ESP/gBnH8+zJ9f/bHOPBP+6Z/gueeqP5baZyhI4pvfhLVru3Os+fPhfe+Db3+7O8dTewwFacjt2gXf/S68//3dO+bll8Ptt3fveGqdoSANuVtugQsvhJGR7h3zve8tboz34IPdO6ZaU2koRMRFEfFERIxGxLUNXl8cEbeXr98bEauqrEfSgbZvhy98AT75ye4ed/Fi+J3fgc98BjK7e2xNrbJQiIj5wI3AxcDpwBURcfqkZlcCr2Tm24EvAZ+vqh5JB9q/H66+uvir/V3v6v7xP/KR4qE7N93U/WOruSp7CucCo5n5dGbuAW4DJk9lrQW+Wi7fCVwYEVFhTdLQy4RHHy0mlsfG4I/+qDd1LF4M3/oW/P7vw8c+Bs8/35s6dKAKHrpXcwJQ/0TWMWDy3yO1Npm5NyK2AccAWzpdzM9+Bp9vox/STpe2qrbWYR2dbL9vX/H85cceKx63efXV8PGPw6JF7R2zk1avLq5u/vSni/suHXssnHQSHH00HHpo8WCeiOLzxPLExzD6zGfgLW+p9hhVhkKj/7bJ38KttCEirgKuAjjxxBNnVMwhh8AZZ7T3Ne1841XV1jqso1Pt582D5cuLZxqcdNLc+cU6MgJf+Qp8+cvw+OPF9QsvvQSvvlqE3v79b36eWB5WixdXf4wqQ2EMWFm3vgKYfLnKRJuxiFgAHAm8PHlHmXkzcDPAmjVrZvQtcdxx8OEPz+QrJXXD/PnFH27t/vGmzqpyTuE+YHVEnBwRi4B1wPpJbdYDHyyXLwO+nznMfwdIUm9V1lMo5wiuAe4C5gO3ZOamiLgB2JiZ64H/CXwtIkYpegjrqqpHkjS9KoePyMwNwIZJ266vW34NuLzKGiRJrfOKZklSjaEgSaoxFCRJNYaCJKnGUJAk1US/XRYQEePAz2b45cup4BYaHTBX64K5W5t1tce62jOIdZ2UmdPeIL3vQmE2ImJjZq7pdR2TzdW6YO7WZl3tsa72DHNdDh9JkmoMBUlSzbCFws29LqCJuVoXzN3arKs91tWeoa1rqOYUJElTG7aegiRpCkMTChFxUUQ8ERGjEXFtj2t5JiIeiYgHI2Jjue3oiPjriHiy/HxUF+q4JSJejIhH67Y1rCMKf1yev4cj4pwu1/XZiPhFec4ejIhL6l67rqzriYh4b4V1rYyIuyPi8YjYFBH/udze03M2RV09PWcRcUhE/CgiHirr+m/l9pMj4t7yfN1e3lqfiFhcro+Wr6+qoq5pars1IjbXnbOzyu3d/P6fHxEPRMR3yvXunq/MHPgPilt3PwWcAiwCHgJO72E9zwDLJ237AnBtuXwt8Pku1PErwDnAo9PVAVwCfJfiaXnnAfd2ua7PAh9v0Pb08v9zMXBy+f88v6K6jgfOKZeXAj8tj9/TczZFXT09Z+W/+/ByeSFwb3ke7gDWldtvAq4ulz8M3FQurwNur/B7rFlttwKXNWjfze//jwJ/AXynXO/q+RqWnsK5wGhmPp2Ze4DbgLU9rmmytcBXy+WvAr9e9QEz8+84+El3zepYC/x5Fu4BlkXE8V2sq5m1wG2Z+XpmbgZGKf6/q6jr+cz8cbm8A3ic4jnjPT1nU9TVTFfOWfnv3lmuLiw/Evg14M5y++TzNXEe7wQujKjmoaFT1NZMV/4vI2IF8K+BPyvXgy6fr2EJhROAZ+vWx5j6h6ZqCXwvIu6P4vnTAMdl5vNQ/JADx/aotmZ1zIVzeE3Zdb+lbnitJ3WVXfWzKf7CnDPnbFJd0ONzVg6FPAi8CPw1Ra9ka2bubXDsWl3l69uAY6qoq1FtmTlxzn6vPGdfioiJpyJ365z9D+ATwP5y/Ri6fL6GJRQapWcv33Z1fmaeA1wM/HZE/EoPa2lVr8/hnwBvA84Cngf+sNze9boi4nDgm8B/ycztUzVtsK2y2hrU1fNzlpn7MvMsime0nwucNsWxu3q+JtcWEWcA1wG/BPwz4Gjgk92qLSLeB7yYmffXb57iuJXUNCyhMAasrFtfATzXo1rIzOfKzy8C36b4YXlhojtafn6xR+U1q6On5zAzXyh/iPcDf8qbwx1drSsiFlL84v0/mfmtcnPPz1mjuubKOStr2Qr8gGI8fllETDz1sf7YtbrK14+k9WHETtR2UTkUl5n5OvC/6O45Ox+4NCKeoRji/jWKnkNXz9ewhMJ9wOpyFn8RxaTM+l4UEhGHRcTSiWXgXwGPlvV8sGz2QeAve1HfFHWsBz5QvgvjPGDbxJBJN0wav30/xTmbqGtd+U6Mk4HVwI8qqiEoniv+eGb+97qXenrOmtXV63MWESMRsaxcXgK8h2K+427gsrLZ5PM1cR4vA76f5Sxql2r7SV24B8XYff05q/T/MjOvy8wVmbmK4nfU9zPzN+n2+erUjPlc/6B498BPKcY0P9XDOk6heOfHQ8CmiVooxgL/Bniy/Hx0F2r5OsWwwhsUf3Vc2awOiq7qjeX5ewRY0+W6vlYe9+Hyh+H4uvafKut6Ari4wrouoOiePww8WH5c0utzNkVdPT1nwJnAA+XxHwWur/sZ+BHFBPc3gMXl9kPK9dHy9VMq/L9sVtv3y3P2KPC/efMdSl37/i+P927efPdRV8+XVzRLkmqGZfhIktQCQ0GSVGMoSJJqDAVJUo2hIEmqMRSkKUTEsoj4cLn81oi4c7qvkfqZb0mVplDeS+g7mXlGj0uRumLB9E2kofYHwNvKG6c9CZyWmWdExG9RXPE6HziD4r5Ci4B/D7wOXJKZL0fE2yguehoBdgMfysyfdP+fIbXG4SNpatcCT2Vx47TfnfTaGcC/pbg/zu8BuzPzbOAfgQ+UbW4GPpKZvwx8HPhyV6qWZsiegjRzd2fx/IIdEbEN+H/l9keAM8u7lv5z4Bt1t7lffPBupLnDUJBm7vW65f116/spfrbmUdwL/6xuFybNlMNH0tR2UDzism1ZPNNgc0RcDrXn/L6zk8VJnWYoSFPIzJeAf4iIR4EvzmAXvwlcGRETd8Wda4+BlQ7gW1IlSTX2FCRJNYaCJKnGUJAk1RgKkqQaQ0GSVGMoSJJqDAVJUo2hIEmq+f+h+WtAi4K83wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113aa3f60>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(maxTime), output, color='blue', linewidth=1)\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('Ez')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's put in a slab of glass that's 50 mm thick, located betweeen 150 and 200 mm in our grid.  The relative permittivity of glass (pyrex) is 4.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up glass slab\n",
    "for mm in range(150,200):\n",
    "    epsR[mm] = 4.7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and rerun the code, storing the output in output2 this time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "ez = [0.0] * SIZE\n",
    "hy = [0.0] * SIZE\n",
    "output2 = [0.0] * maxTime\n",
    "\n",
    "# do time stepping\n",
    "for qTime in range(maxTime):\n",
    "\n",
    "    # update magnetic field\n",
    "    for mm in range(SIZE-1):\n",
    "        hy[mm] = hy[mm] + (ez[mm + 1] - ez[mm]) / imp0\n",
    "\n",
    "    # update electric field\n",
    "    for mm in range(SIZE):\n",
    "        ez[mm] = ez[mm] + (hy[mm] - hy[mm - 1]) * imp0 / epsR[mm]\n",
    "\n",
    "    # hardwire a source node */\n",
    "    if qTime < sourceSigma:\n",
    "        ez[0] = math.exp(-(qTime - sourcePeakTime)**2 / sourceSigma)\n",
    "    else:\n",
    "        ez[0] = 0.0\n",
    "\n",
    "    output2[qTime] = ez[sensorLocation]\n",
    "#done with time stepping loop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's plot the two outputs, for free space in blue and with a glass slab in green "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3X1wXXd95/H3V4+WLFnyg+I4thI7wd5i0jxhAiUpzQZ2SbLbmB0CTaAk2U2T2bbp7g6lEIadwNJh2sJ02e1MWki7bChLCeGh1GXMhA4JpANNiPNInEfFD1jxkyzryXqyZH33j3Pu9dXV1b3nSufce670ec14fO+5R+d8fSzr49/3dx7M3REREQGoq3YBIiKSHgoFERHJUiiIiEiWQkFERLIUCiIikqVQEBGRLIWCiIhkKRRERCRLoSAiIlkN1S6gXOvWrfPNmzdXuwwRkZry1FNPnXD3rlLr1VwobN68mT179lS7DBGRmmJmB6Osp/aRiIhkKRRERCRLoSAiIlkKBRERyVIoiIhIVmKhYGZfMbPjZvbCPJ+bmf2FmfWY2fNmdkVStYiISDRJjhQeAK4r8vn1wNbw113AXyVYi4iIRJBYKLj7Y8DJIqvsBP7WA48DnWa2Ial6RCR+/f1w+nS1q5A4VXNOYSNwKOd9b7hsDjO7y8z2mNmevr6+ihQnIqW9//3wx39c7SokTtUMBSuwzAut6O73u/sOd9/R1VXyKm0RqYCeHnjqKXjgAThzptrVSFyqGQq9QHfO+03A4SrVIiJlevhh+K3fgtZW2Lu32tVIXKoZCruAW8OzkN4BDLn7kSrWIyJl2LcPtm2DrVuD17I0JHZDPDP7BnANsM7MeoFPA40A7v4lYDdwA9ADjAH/MalaRCR++/fDr/0a9PYGr2VpSCwU3P2WEp878PtJ7V9EkrV/P2zZAocOKRSWEl3RLCILkgmFLVsUCkuJQkFEyjYwAO6werVCYalRKIhI2Xp7obsbzILfe3urXZHERaEgImXr64PMJUOdnTA6qiublwqFgoiULTcU6upg7drglhdS+xQKIlK2Eydg3bqz79etC4JCap9CQUTKljtSgOC1QmFpUCiISNn6+maPFLq6gtGD1D6FgoiUTSOFpUuhICJlO3FidihoTmHpUCiISNnUPlq6FAoiUrb8s4/UPlo6FAoiUraBgeAWFxmdnTA4WL16JD4KBREpy8QEzMxAS8vZZZ2dMDRUvZokPgoFESnL0FAQApbzQF2NFJYOhYKIlGVwMAiBXB0dCoWlQqEgImUZGgpCIJdGCkuHQkFEylJopLBiRdBOmpioTk0SH4WCiJSlUCiAWkhLhUJBRMqSmWjOpxbS0qBQEJGyzDdSUCgsDQoFESnL4ODciWZQKCwVCgURKct87aOOjtq4gG1gfIDhyeFql5FaCgURKUux9tHAQOXrKdcdu+7gvD8/j18O/bLapaSSQkFEyjI0BKtWzV2+ahWMjFS+nnJMz0zz6IFHuWLDFTx28LFql5NKCgURKcvISOFQaG9Pfyg8feRpNrZv5MZ/dSNP9D5R7XJSSaEgImUZGQkCIF9bW/pD4ck3nuSq7qt4+8a38/PDP692OamkUBCRsgwPFw6F9nY4dary9ZRj38A+LlpzEZeeeykvHH8Bd692SamjUBCRssw3UqiF9tG+wX1cuPpCVjWvoqGugYGJGpgZrzCFgoiUZb45hba29I8U9g/sZ0vnFgC6V3XTO9xb5YrSJ9FQMLPrzOwVM+sxs3sKfH6+mT1qZs+Y2fNmdkOS9YjI4kxPBze9a22d+1naRwruzr6BYKQAsGnVJg4NHapyVemTWCiYWT1wH3A9sB24xcy2563234GH3P1y4GbgL5OqR0QW79SpYESQ+4CdjLSPFE6On6TO6ljdEjxHtHtVN4eGFQr5khwpXAn0uPs+dz8NPAjszFvHgcxAtAM4nGA9IrJI880nQPpHCoeGD9Hd0Z19393RrZFCAUmGwkYg94j3hstyfQb4bTPrBXYDf5BgPSKySPOdeQTpD4Vjp45xbtu52fcaKRSWZCgUGGCSf/7XLcAD7r4JuAH4mpnNqcnM7jKzPWa2p6+vL4FSRSSK+SaZIf3to2Ojx1i/cn32/Yb2DRwbPVbFitIpyVDoBbpz3m9ibnvoDuAhAHf/F2AFsC5/Q+5+v7vvcPcdXV1dCZUrIqUUax+1tcHYGMzMVLamqI6eOjprpLCudR0nxk5UsaJ0SjIUngS2mtkWM2simEjelbfOL4F3A5jZmwlCQUMBkZQqFgp1ddDSEgRDGh07NXukoFAoLLFQcPdp4G7gYeAlgrOM9prZZ83sxnC1PwTuNLPngG8At7suMRRJrWJzCpDuW10cGz2mkUIEDUlu3N13E0wg5y67N+f1i8BVSdYgIvEpNlKAs5PNGzZUrqaojp46yvq2syOFlY0rOTNzhvGpcVoaW6pYWbroimYRiezUqdIjhbRONudPNJsZ61rX0T/eX8Wq0kehICKRZS5em0+aT0s9Pnqcc1aeM2uZWkhzKRREJLLRUVi5cv7P03qnVHfn5PhJ1raunbVcoTCXQkFEIis1UkjrRPPw5DAtDS001TfNWq5QmEuhICKRjY6Wbh+lcaTQP94/Z5QAQSj0jeos+FwKBRGJ7NSp0u2jNI4U+sf6WdOyZs7y1StWMzgxWIWK0kuhICKR1Wr7qH+8n7Utc0cKnSs6FQp5FAoiElnNto/GCrePFApzKRREJLJS7aNaGyl0rOhgcFKhkEuhICKRRblOIbUjhXnaR0MTQ1WoKL0UCiISWZT2URpHCifHTxacaFb7aC6FgohEVtPtI80pRKJQEJFIpqbgzBlobp5/nbS2jwYnBlm9YvWc5QqFuRQKIhJJpnVkhZ6pGErrSGFwYpCOFR1zlnc0dzA0OYTu2H+WQkFEIinVOoL0jhSGJofoaJ4bCo31jTTXNzM6NVqFqtJJoSAikZSaZIb0TjQPTgzSuaKz4GdqIc2mUBCRSEqdjgpnn6eQtm7M0MRQwfYRKBTyKRREJJIo7aPGxuBZzadPV6amKKZnphmfHqetqXCirWpexfDkcIWrSi+FgohEEqV9BOl7+trQxBCrmldRZ4V/3LU3tysUcigURCSSKO0jSGEozDPJnNHe1M7IZAonQqpEoSAikURpH0EKQ2FiaN5JZghGCiOnFQoZCgURiaSc9tFois7wnO8ahQyNFGZTKIhIJDU7UijRPlrVvEojhRwKBRGJpGbnFEq1jzRSmEWhICKRRG0frVyZrlAYnBgsPtGss49mUSiISCQ13T4qNaeg9lGWQkFEIlmy7SOdfTSLQkFEIqnVi9dKto80pzCLQkFEIlmy7SONFGZRKIhIJLU6UhiaLN4+0r2PZlMoiEgktTqnoPZReRINBTO7zsxeMbMeM7tnnnU+aGYvmtleM/u7JOsRkYUrp32Upiuai902G9Q+yteQ1IbNrB64D/g3QC/wpJntcvcXc9bZCnwSuMrdB8zsnKTqEZHFWarto7amNsamxpjxmXnvpLqcJHkErgR63H2fu58GHgR25q1zJ3Cfuw8AuPvxBOsRkUVYqu2jOqujpaGF0dMpGt5UUZKhsBE4lPO+N1yWaxuwzcx+amaPm9l1hTZkZneZ2R4z29PX15dQuSIyn6kpOHMGmppKr5umK5onpidwd1Y0rCi6nlpIZyUZClZgWf5D+hqArcA1wC3A35jZnHGeu9/v7jvcfUdXV1fshYpIcZnWkRX6V50nTSOFzLOZrUTh7U261UVGkqHQC3TnvN8EHC6wzj+4+5S77wdeIQgJEUmRqK0jSFcolJpkzljVvEpnIIWSDIUnga1mtsXMmoCbgV1563wP+NcAZraOoJ20L8GaRGQBop55BCkLhRKTzBlqH52VWCi4+zRwN/Aw8BLwkLvvNbPPmtmN4WoPA/1m9iLwKPBH7t6fVE0isjBRzzwCaGmBiYlgDqLaSk0yZ+hahbMSOyUVwN13A7vzlt2b89qBj4a/RCSlymkf1dVBayuMjUF7e7J1lRK1faSRwlk6KVdEShodjd4+gvS0kIYmh+hsjtA+0kRzlkJBREoqZ6QA6bmqudTzmTPUPjpLoSAiJS0kFFIxUpgo/nzmDD2n+SyFgoiUVNPto6hnH2mkACgURCSCckcKabmquaz2kUYKgEJBRCKo2fbRZLT2kc4+OkuhICIl1Wz7qMTzmTN09tFZCgURKalWRwqR20fN7Zw6nYKCU0ChICIl1WooRG0ftTW1aaI5pFAQkZKWQ/tII4WAQkFESqrFkYK7MzQ5xKrmVSXX1UTzWZFCwczOmNmfWs5Nyc3s6eTKEpE0GRkp7z5GaQiFU6dPsaJhBY31jSXXbWtq00ghFHWksDdc94dmtiZcFuFxGyKyFNTibS6iXrgG0NLQwukzp5memU64qvSLGgrT7v5x4K+BfzaztzL3KWoiskTVYvso6m2zAcxMo4VQ1FtnG4C7P2Rme4FvAOcnVpWIpEq57aM0XNEc9bbZGZmb4kUdXSxVUUPhdzIvwgflXA28L5mSRCRN3Mt78hqkZ6RQzg94jRQCRdtHZvZxAHd/ysw+kFnu7sPAryRcm4ikwMQENDVBQxmP5EpDKES9RiFDZyAFSs0p3Jzz+pN5n10Xcy0ikkLlto4gJaEQ8RqFDF3AFigVCjbP60LvRWQJKneSGdIRCuVMNIMuYMsoFQo+z+tC70VkCRoZKT8UMhPNXsWfEkOT5U00tzW1qX1E6YnmS81smGBU0BK+Jny/ItHKRCQVTp0qv33U1AR1dXD6NDQ3J1NXKUMTQ5zfEf0kSY0UAkVDwd3rK1WIiKTTQtpHcLaFVK1QGJwss32kp68BuveRiJSwkPYRVP+q5nKvU9ApqQGFgogUtZD2EVR/srmc21yAHsmZoVAQkaIW2z6qlnLPPtJIIaBQEJGiFto+qvatLsq+zYUuXgMUCiJSwnJpH2mkEFAoiEhRtdg+mp6ZZmxqjLam6IVnboi33CkURKSoxZx9VK1QGJ4cpr2pnTqL/iNOF68FFAoiUlQtto/Kve8RBHMKah8lHApmdp2ZvWJmPWZ2T5H1bjIzN7MdSdYjIuWrxZHC4MRgWZPMoBviZSQWCmZWD9wHXA9sB24xs+0F1msH/gvwRFK1iMjC1eKcQrm3zQbd5iIjyZHClUCPu+9z99PAg8DOAuv9MfB5YCLBWkRkgRbTPqrWFc0LaR+1NbUxOjWKV/MufimQZChsBA7lvO8Nl2WZ2eVAt7t/P8E6RGQRlkv7qL6unub6ZsamxhKqqjYkGQqFnreQjWAzqwO+CPxhyQ2Z3WVme8xsT19fX4wlikgptdo+6mwu/1nLuoAt2VDoBbpz3m8CDue8bwcuBn5sZgeAdwC7Ck02u/v97r7D3Xd0dXUlWLKI5Fto+6iaVzQvZKQAuoANkg2FJ4GtZrbFzJoIHu25K/Ohuw+5+zp33+zum4HHgRvdfU+CNYlIGWZmYGwMWlvL/9pqjhQGxgdYvWJ12V+nC9gSDAV3nwbuBh4GXgIecve9ZvZZM7sxqf2KSHxGR4NAqFvAT4qqhsLEAKtbyg8FXcBW+slri+Luu4HdecvunWfda5KsRUTKNzKysNYRVD8U1rSsKfvrdAGbrmgWkSIWOskM1Q2Fk+MnF9Q+0gVsCgURKaJWQ2FgfGHtI13AplAQkSIWeo0CBHMRExNw5ky8NUUxMLGwiWbNKSgURKSIoSHoKP/MTgDMgkAZqcLPWI0UFk6hICLzWkwoQPC1Q0Px1RPF5PQk0zPTrGxcWfbXtjfrlFSFgojMa2gIOsu/MDirGqGQOR3VrNBNFYrTxWsKBREpYrEjhc7OyofCQs88gvDiNc0piIgUFkf7aHAwvnqiWOh8AmikAAoFESlicLD25hQWeuYR6IZ4oFAQkSJqcaJ5YHxhVzODLl4DhYKIFFGTobCYkYJOSVUoiMj8avLso0XOKah9JCIyj1ocKSzq7CPdEE+hICLzq8mzjxZ4h1TQnAIoFESkiFq8TmGhz1IAaK5vBoKropcrhYKIFDQ1BePjC78hHlRxTmGB7SMzo2NFB0OTFS46RRQKIlLQ8DCsWhXc2G6hqnmbi4XqaO5gcKLCPa8UUSiISEGLPfMIam+kANC5opOhCY0URERmWex8AlRnovnk+MkFTzQDah9VuwARSac4QqGtLZiXmJ6Op6ZSRk+PYma0NLYseBtqH4mIFBBHKNTVQXt7MD9RCX1jfZyz8pxFbUPtIxGRAuIIBajsaanHR4/T1dq1qG10NKt9JCIyx+Dg4ieaobKTzX2j8YwU1D4SEckT10ihoqEw1kfXykWOFFZ0qH0kIpIvzlCo1BlIx0ePc07r4kYKah+JiBRQkyOF0cWPFNQ+EhEpoBZD4fjY8UXPKeg6BRGRAmoxFPpG+2I5+0gjBRGRPHGdfdTZWdk5hTjaR5poFhHJ098Pa9cufjtr1wbbqoQ4Ll5T+0hEpIBaCwV3j+/itYkh3D2mympLoqFgZteZ2Stm1mNm9xT4/KNm9qKZPW9mPzKzC5KsR0SimZ6GkZF42keVCoXRqVHqrI6VTSsXtZ3G+kaaG5oZnRqNqbLaklgomFk9cB9wPbAduMXMtuet9gyww90vAb4NfD6pekQkupMng0Coi+EnRKVCIY5RQkZmtLAcJTlSuBLocfd97n4aeBDYmbuCuz/q7mPh28eBTQnWIyIRxdU6Ali3rjKhEMctLjKW87UKSYbCRuBQzvvecNl87gB+UOgDM7vLzPaY2Z6+vr4YSxSRQuIMhTVrYGAAZmbi2d584jjzKGM5TzYnGQqFHuJXcObGzH4b2AF8odDn7n6/u+9w9x1dXfH8pYvI/OIMhYaG4LkKSZ+W2je2+GsUMpbzaalJhkIv0J3zfhNwOH8lM3sP8CngRnefTLAeEYkozlCAyswrxNk+Ws4XsCUZCk8CW81si5k1ATcDu3JXMLPLgS8TBMLxBGsRkTLUYigcGz0WayiofRQzd58G7gYeBl4CHnL3vWb2WTO7MVztC0Ab8C0ze9bMds2zORGpoBMn4g2FSkw2vzHyBhvbi01bRrec20cNSW7c3XcDu/OW3Zvz+j1J7l9EFub4cdi2Lb7tnXNOsM0kvTH8Bue1nxfLtjpWdDAwMRDLtmqNrmgWkTmOHYP16+Pb3vr1wTaTdHjkMBtXxTNSWNuylpPjJ2PZVq1RKIjIHEePxh8KR4/Gt7187s7hkcOxjRTWta6jf7xCN2xKGYWCiMwR90jhnHOSHSn0j/fT0thCa2NrLNtb27qWE2MnYtlWrVEoiMgsMzPQ1xf8II9L0u2jwyOHY5tkhmCkoFAQESG4+ri1FVasiG+bSYdCnJPMoFAQEcmKu3UEyYfCwaGDXNAR302WMxPNy/H22QoFEZkliVBYty64zcXUVLzbzTgweIDNnZtj215jfSOtja3L8gI2hYKIzHLkCGzYEO826+uTnWw+OHQw1lCA5dtCUiiIyCy9vdDdXXq9cm3aFGw7CQcGD3BBZ7zP6FIoiIgAhw4FP8Dj1t0dbDsJcbePQKEgIgLU3khhYnqCk+Mn2dAWb89r/cr1HD2V4BV3KaVQEJFZam2k8PrJ19nSuYX6uvpYt7uhbQNHRo7Eus1aoFAQkVkOHUpupJBEKLza/yrb1sZ4977QhvYNHDmlUBCRZWxyMrh4Le5TUiG5kUJiodCmUBCRZe7gweB/9PXxdmIA2LIFDhyIf7uJjhTUPhKR5aynB7ZuTWbbGzbAyAgMD8e73VdPJjdS0ESziCxrPT3wpjcls20zuOgieP31+Lbp7rxw/AXe0vWW+DYaOrftXI6eOrrsbnWhUBCRrCRDAYJt9/TEt71Dw4doaWiha2VXfBsNtTS20NLYsuyeq6BQEJGsWguFZ48+yyXrL4lvg3k2d27m4ODBxLafRgoFEcl6+eV4n82cb9u2YB9xee7oc1y6/tL4NphnS+cW9g/uT2z7aaRQEBEgmAA+dizZkcIll8Bzz8W3vT1H9vDW894a3wbzbO7czP4BhYKILEPPPw8XX5zM6agZF18Mr7wCp08vflvuzs8O/Yx3dr9z8Rubh0YKIrJsPfccXJpcJwYInuh2wQXxtJBe7X+V1sZWNq1K4J4coS2rt3Bg8EBi208jhYKIAPD003DZZcnv5/LLg30t1mMHH+Pq869e/IaKuHD1hfScjHFmvAYoFEQEgMceg1//9eT3c/XVwb4Wa3fPbq5/0/WL31ARW9dspXe4l7GpsUT3kyYKBRHh8GE4eRLeEv81YHP8xm/AT36yuG1MTk/yyP5HeO9F742nqHk01jeybe02Xux7MdH9pIlCQUR45BF417ugrgI/EbZvD850OriI0/9/0PMDLl1/aSIXreW7ZP0lPH/s+cT3kxYKBRHhO9+B972vMvuqq4OdO4N9LtQDzz7A7ZfdHltNxVyy/hKePfpsRfaVBgoFkWVucDAYKezcWbl9fvCD8OCDC/va1/pf46eHfsoHtn8g3qLmcfX5V/OTg4vsd9UQhYLIMvflL8Nv/iZ0dlZun9deCydOwM9+Vv7XfvrHn+but91Ne3N7/IUVsOO8HRwcPMixU8cqsr9qSzQUzOw6M3vFzHrM7J4Cnzeb2TfDz58ws81J1iMis504AV/8InziE5Xdb0MDfPzj8KlPwcxM9K/73svf4/Hex/nYOz+WXHF5GuoauHbLtex+bXfF9llNiYWCmdUD9wHXA9uBW8xse95qdwAD7v4m4IvAnyVVj4jMdvo03HYbfPjD8Ku/Wvn933knTEzA5z4Xbf0fvv5D7vzHO3nwpgdZ2bQy2eLy3HH5Hdz35H3L4jbaSY4UrgR63H2fu58GHgTyu5Y7ga+Gr78NvNvMLMGaRJa9yUn40Y/gmmugqQn+5E+qU0d9fTDZ/PWvw0c+ElxRnT9qGJ8a55H9j3Dr39/K7d+7ne9+8LtcufHKitd6/dbrmTwzyQPPPlDxfVeaJZV8ZnYTcJ27/074/iPA29397px1XgjX6Q3fvx6uc2K+7e7YscP37NlTdj0/eqaHW78WdYzsRD4s5hBhXTfPbDrS/qPyyOt6kXeV2H/07ca1zbmfFNlu7kdWev/Z748I61LGumVtl/KOlbszNRX877y1FTZuhA3nAgX+Gxb150K5+y/kzBk41At9x2FqCpqbob55nKkVRzjdcIKOscs5d2gnF/b/Lk3egVnwwJ5KG25+kX/e8m5Wj7+N5ukuputGmbH8mzgVL8xKfF7K9/7z/2bHto0L+loze8rdd5Rar2FBW49YQ4Fl+d8VUdbBzO4C7gI4//zzF1TM5vVr+NDFHyrjKyzyN55hpb4Xzq5H1G/osyuVWn/WN1oZ6xZbtZwBW6F15//qAusWWLnoPx7Lf1tkb2WsO/swlv7zn9129HXj3i5AXcS/KzPo7DTOXR+EQun1I263jB90pbY5MAD9/TA+0sxKP4+O+vXU04R7MIqYmSH6f9hit50Pn3mZZ07tZnJmlOa6lTRYU/bPXyog4/gP+Lmrk59cTzIUeoHunPebgMPzrNNrZg1AB3Ayf0Pufj9wPwQjhYUUc9F5a/jCf3r/Qr5URCrlgmoXUEoHt3FLtYtIVJJzCk8CW81si5k1ATcDu/LW2QXcFr6+CXjEl8NMjohISiU2UnD3aTO7G3gYqAe+4u57zeyzwB533wX8H+BrZtZDMEK4Oal6RESktCTbR7j7bmB33rJ7c15PAJW5LFFERErSFc0iIpKlUBARkSyFgoiIZCkUREQkS6EgIiJZid3mIilm1gcs9JlN64B5b6FRRWmtC9Jbm+oqj+oqz1Ks6wJ3L/moupoLhcUwsz1R7v1RaWmtC9Jbm+oqj+oqz3KuS+0jERHJUiiIiEjWcguF+6tdwDzSWhektzbVVR7VVZ5lW9eymlMQEZHilttIQUREilg2oWBm15nZK2bWY2b3VLmWA2b2CzN71sz2hMvWmNk/mdlr4e+rK1DHV8zsePgEvMyygnVY4C/C4/e8mV1R4bo+Y2ZvhMfsWTO7IeezT4Z1vWJm702wrm4ze9TMXjKzvWb2X8PlVT1mReqq6jEzsxVm9nMzey6s63+Ey7eY2RPh8fpmeGt9zKw5fN8Tfr45ibpK1PaAme3POWaXhcsr+f1fb2bPmNn3w/eVPV7uvuR/Edy6+3XgQqAJeA7YXsV6DgDr8pZ9HrgnfH0P8GcVqONdwBXAC6XqAG4AfkDwOLB3AE9UuK7PAB8rsO728O+zGdgS/j3XJ1TXBuCK8HU78Gq4/6oesyJ1VfWYhX/utvB1I/BEeBweAm4Ol38J+N3w9e8BXwpf3wx8M8HvsflqewC4qcD6lfz+/yjwd8D3w/cVPV7LZaRwJdDj7vvc/TTwILCzyjXl2wl8NXz9VeB9Se/Q3R9j7pPu5qtjJ/C3Hngc6DSzDRWsaz47gQfdfdLd9wM9BH/fSdR1xN2fDl+PAC8BG6nyMStS13wqcszCP/ep8G1j+MuBa4Fvh8vzj1fmOH4beLdZMk9jLlLbfCryd2lmm4B/B/xN+N6o8PFaLqGwETiU876X4v9okubAD83sKQuePw2w3t2PQPCPHDinSrXNV0cajuHd4dD9KznttarUFQ7VLyf4H2ZqjlleXVDlYxa2Qp4FjgP/RDAqGXT36QL7ztYVfj4ErE2irkK1uXvmmH0uPGZfNLPm/NoK1B2n/wV8HJgJ36+lwsdruYRCofSs5mlXV7n7FcD1wO+b2buqWEtU1T6GfwVcBFwGHAH+PFxe8brMrA34DvDf3H242KoFliVWW4G6qn7M3P2Mu19G8Iz2K4E3F9l3RY9Xfm1mdjHwSeBXgLcBa4BPVKo2M/v3wHF3fyp3cZH9JlLTcgmFXqA75/0m4HCVasHdD4e/Hwf+nuAfy7HMcDT8/XiVypuvjqoeQ3c/Fv4jngH+mrPtjorWZWaNBD94v+7u3w0XV/2YFaorLccsrGUQ+DFBP77TzDJPfczdd7au8PMOorcR46jturAV5+4+CfxfKnvMrgJuNLMDBC3uawlGDhU9XsslFJ4Etoaz+E0EkzK7qlGIma00s/bMa+DfAi+E9dwWrnYb8A/VqK9IHbuAW8OzMN4BDGVaJpWQ178cw5SbAAACbklEQVT9DwTHLFPXzeGZGFuArcDPE6rBCJ4r/pK7/8+cj6p6zOarq9rHzMy6zKwzfN0CvIdgvuNR4KZwtfzjlTmONwGPeDiLWqHaXs4JdyPo3eces0T/Lt39k+6+yd03E/yMesTdP0ylj1dcM+Zp/0Vw9sCrBD3NT1WxjgsJzvx4DtibqYWgF/gj4LXw9zUVqOUbBG2FKYL/ddwxXx0EQ9X7wuP3C2BHhev6Wrjf58N/DBty1v9UWNcrwPUJ1nU1wfD8eeDZ8NcN1T5mReqq6jEDLgGeCff/AnBvzr+BnxNMcH8LaA6Xrwjf94SfX5jg3+V8tT0SHrMXgP/H2TOUKvb9H+7vGs6efVTR46UrmkVEJGu5tI9ERCQChYKIiGQpFEREJEuhICIiWQoFERHJUiiIFGFmnWb2e+Hr88zs26W+RqSW6ZRUkSLCewl9390vrnIpIhXRUHoVkWXtT4GLwhunvQa82d0vNrPbCa54rQcuJrivUBPwEWASuMHdT5rZRQQXPXUBY8Cd7v5y5f8YItGofSRS3D3A6x7cOO2P8j67GPgQwf1xPgeMufvlwL8At4br3A/8gbu/FfgY8JcVqVpkgTRSEFm4Rz14fsGImQ0B/xgu/wVwSXjX0ncC38q5zX3z3M2IpIdCQWThJnNez+S8nyH4t1VHcC/8yypdmMhCqX0kUtwIwSMuy+bBMw32m9kHIPuc30vjLE4kbgoFkSLcvR/4qZm9AHxhAZv4MHCHmWXuipu2x8CKzKJTUkVEJEsjBRERyVIoiIhIlkJBRESyFAoiIpKlUBARkSyFgoiIZCkUREQkS6EgIiJZ/x+HSgLJeQ42bgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x116ceff98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(maxTime), output, color='blue', linewidth=1)\n",
    "plt.plot(range(maxTime), output2, color='green', linewidth=1)\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('Ez')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that the second (green) curve has been reduced in height and has been delayed by passing through the glass slab.  The reduction in height is a property of going from one material to another - some energy is reflected by at the interface between the two materials. The delay is a property of the glass itself - waves travel more slowly in glass than in vaccum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's calculate the sum of the electric field that passed through free space and compare it to the sum of the electric field that passed through the glass slab.  The ratio is the transmission, and the part of the field that didn't pass through the slab was reflected, and is the reflectance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Transmission = 0.8640776048486403\n",
      "Reflectance = 0.13592239515135973\n"
     ]
    }
   ],
   "source": [
    "base = 0\n",
    "transmitted = 0\n",
    "for qTime in range(maxTime):\n",
    "    base += output[qTime]\n",
    "    transmitted += output2[qTime]\n",
    "transmission = transmitted/base\n",
    "print('Transmission =',transmission)\n",
    "print('Reflectance =', (1 - transmission))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if we run the simulation for one more time step?  Remember that the source is turned off after sourceSigma time steps, so this does not inject any more energy into the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "newMaxTime = 501 # simulation time\n",
    "\n",
    "ez = [0.0] * SIZE\n",
    "hy = [0.0] * SIZE\n",
    "output3 = [0.0] * newMaxTime\n",
    "\n",
    "# free space\n",
    "epsR = [1.0] * SIZE\n",
    "\n",
    "# do time stepping\n",
    "for qTime in range(newMaxTime):\n",
    "\n",
    "    # update magnetic field\n",
    "    for mm in range(SIZE-1):\n",
    "        hy[mm] = hy[mm] + (ez[mm + 1] - ez[mm]) / imp0\n",
    "\n",
    "    # update electric field\n",
    "    for mm in range(SIZE):\n",
    "        ez[mm] = ez[mm] + (hy[mm] - hy[mm - 1]) * imp0 / epsR[mm]\n",
    "\n",
    "    # hardwire a source node */\n",
    "    if qTime < sourceSigma:\n",
    "        ez[0] = math.exp(-(qTime - sourcePeakTime)**2 / sourceSigma)\n",
    "    else:\n",
    "        ez[0] = 0.0\n",
    "\n",
    "    output3[qTime] = ez[sensorLocation]\n",
    "#done with time stepping loop\n",
    "\n",
    "\n",
    "ez = [0.0] * SIZE\n",
    "hy = [0.0] * SIZE\n",
    "output4 = [0.0] * newMaxTime\n",
    "\n",
    "# add a glass slab\n",
    "for mm in range(150,200):\n",
    "    epsR[mm] = 4.7\n",
    "\n",
    "# do time stepping\n",
    "for qTime in range(newMaxTime):\n",
    "\n",
    "    # update magnetic field\n",
    "    for mm in range(SIZE-1):\n",
    "        hy[mm] = hy[mm] + (ez[mm + 1] - ez[mm]) / imp0\n",
    "\n",
    "    # update electric field\n",
    "    for mm in range(SIZE):\n",
    "        ez[mm] = ez[mm] + (hy[mm] - hy[mm - 1]) * imp0 / epsR[mm]\n",
    "\n",
    "    # hardwire a source node */\n",
    "    if qTime < sourceSigma:\n",
    "        ez[0] = math.exp(-(qTime - sourcePeakTime)**2 / sourceSigma)\n",
    "    else:\n",
    "        ez[0] = 0.0\n",
    "\n",
    "    output4[qTime] = ez[sensorLocation]\n",
    "#done with time stepping loop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's calculate transmission and reflectance again - they should be the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Transmission2 = 0.8640775949984995\n",
      "Reflectance2 = 0.13592240500150055\n"
     ]
    }
   ],
   "source": [
    "base2 = 0\n",
    "transmitted2 = 0\n",
    "for qTime in range(newMaxTime):\n",
    "    base2 += output3[qTime]\n",
    "    transmitted2 += output4[qTime]\n",
    "transmission2 = transmitted2/base2\n",
    "print('Transmission2 =',transmission2)\n",
    "print('Reflectance2 =', (1 - transmission2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Are they the same?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Transmission = 0.8640776048486403\n",
      "Reflectance = 0.13592239515135973\n"
     ]
    }
   ],
   "source": [
    "print('Transmission =',transmission)\n",
    "print('Reflectance =', (1 - transmission))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not quite.  How different are they?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.1399602025146226e-08\n"
     ]
    }
   ],
   "source": [
    "difference = (transmission2 - transmission) / transmission\n",
    "print(difference)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. plot an animation of the fields\n",
    "1. bring in C version?\n",
    "1. library?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Planning notes\n",
    "\n",
    "* 1-D wave equation\n",
    " * Code in python\n",
    " * Define problem (inputs, outputs)\n",
    "* Use Jupyter notebook on workbench.\n",
    "\n",
    "* Talk - Highlight Importance of scriptability \n",
    "* Talk - Compare GUI and script. \n",
    "\n",
    "* Code in C\n",
    "* Talk: Explain difficulties in how to compile and how to run\n",
    " * Compile and code in notebook\n",
    "   * Add a Library dependency in the code \n",
    "* Maybe a wave equation library?\n",
    "* Write output.\n",
    "* Read input and plot \n",
    "\n",
    "* Could also suggest they download code and compile and run themselves \n",
    "\n",
    "\n",
    "* also talk about what can fail\n",
    "\n",
    "* time basis\n",
    "\n",
    "* software collapse\n",
    "* containers\n",
    "* VMs\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ims = []   \n",
    "\n",
    "# set up a plot\n",
    "fig, ax = plt.subplots()\n",
    "ax.axvspan(150, 200, alpha=0.1, color='red') # glass slab\n",
    "ax.axvspan(250, 250, alpha=0.9, color='black') # sensor\n",
    "ax.set_xlabel('time')\n",
    "ax.set_ylabel('Ez')\n",
    "\n",
    "\n",
    "ez = [0.0] * SIZE\n",
    "hy = [0.0] * SIZE\n",
    "\n",
    "# do time stepping\n",
    "for qTime in range(maxTime):\n",
    "\n",
    "    # update magnetic field\n",
    "    for mm in range(SIZE-1):\n",
    "        hy[mm] = hy[mm] + (ez[mm + 1] - ez[mm]) / imp0\n",
    "\n",
    "    # update electric field\n",
    "    for mm in range(SIZE):\n",
    "        ez[mm] = ez[mm] + (hy[mm] - hy[mm - 1]) * imp0 / epsR[mm]\n",
    "    # capture a snapshot of the ez field at this timestep to the animation\n",
    "    ims.append((plt.plot(range(SIZE), ez, color='green', linewidth=1)))\n",
    "\n",
    "    # hardwire a source node */\n",
    "    if qTime < sourceSigma:\n",
    "        ez[0] = math.exp(-(qTime - sourcePeakTime)**2 / sourceSigma)\n",
    "    else:\n",
    "        ez[0] = 0.0\n",
    "\n",
    "#done with time stepping loop\n",
    "\n",
    "#build and display the animation\n",
    "im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=5000, blit=True)\n",
    "HTML(im_ani.to_jshtml())\n",
    "\n"
   ]
  },
  {
   "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
}