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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "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 see what's happening by visualizing the Ez field over the simulation, first with no slab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ims = []   \n",
    "\n",
    "# set up a plot\n",
    "fig, ax = plt.subplots()\n",
    "ax.axvspan(250, 250, alpha=0.9, color='black') # sensor\n",
    "ax.set_xlabel('space (mm)')\n",
    "ax.set_ylabel('Ez')\n",
    "\n",
    "\n",
    "ez = [0.0] * SIZE\n",
    "hy = [0.0] * SIZE\n",
    "\n",
    "# free space\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",
    "    # capture a snapshot of the ez field at this timestep to the animation\n",
    "    ims.append((plt.plot(range(SIZE), ez, color='blue', 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",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and then with a glass slab"
   ]
  },
  {
   "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('Space (mm)')\n",
    "ax.set_ylabel('Ez')\n",
    "\n",
    "\n",
    "ez = [0.0] * SIZE\n",
    "hy = [0.0] * SIZE\n",
    "\n",
    "# free space\n",
    "epsR = [1.0] * SIZE\n",
    "\n",
    "# add a glass slab\n",
    "for mm in range(150,200):\n",
    "    epsR[mm] = 4.7\n",
    "\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": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at this numerically, not graphically.  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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "difference = (transmission2 - transmission) / transmission\n",
    "print(difference)"
   ]
  },
  {
   "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
}