We are going to run an electromagnetic simulation, with initial code and ideas borrowed from
We're going to run a simple, bare-bones 1D FDTD simulation with a hard source.
The impedance of free space (or vacuum) is 377.0.
We are going to model 400 mm of space, and run the simulation for 500 time units.
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.)
Our source is the electric field at the left edge of the grid. We use a Gaussian source that peaks at 30 time units and has a standard deviation of 7.
Our output will be the electric field measured at 250 mm over the time of the simulation.
import math
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
imp0 = 377.0 # property of free space (vacuum)
SIZE = 400 # dimension of space to model
sensorLocation = 250 # location of output sensor
maxTime = 700 # simulation time
sourcePeakTime = 30 # peak of the Gaussian source
sourceSdv = 7 # standard deviation of the Gaussian source
sourceSigma = 2 * sourceSdv**2
We first initialize lists for the electric and magnetic field components.
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.
We then repeat this for each time step, until the set number of time steps have been reached.
ez = [0.0] * SIZE
hy = [0.0] * SIZE
# do time stepping
for n in range(maxTime):
# update magnetic field
for i in range(SIZE-1):
hy[i] = hy[i] + (ez[i + 1] - ez[i]) / imp0
# update electric field
for i in range(SIZE):
ez[i] = ez[i] + (hy[i] - hy[i - 1]) * imp0
# hardwire a source node */
if n < sourceSigma:
ez[0] = math.exp(-(n - sourcePeakTime)**2 / sourceSigma)
else:
## PEC Boundaries
ez[0] = 0.0
ez[(SIZE - 1)] = 0.0
print(ez[sensorLocation])
#done with time stepping loop
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00010270256461965522 0.00018751857457483846 0.0003354626279025118 0.0005880047627378895 0.0010098442341126815 0.001699279365552657 0.0028016380348484843 0.004525807866766719 0.007163364470919217 0.011108996538242306 0.01687988414878991 0.025130488985274417 0.036658041953378004 0.05239314106982558 0.0733696513683829 0.10066889977289778 0.1353352832366127 0.17826397958504792 0.23006629899380904 0.2909238070461482 0.360447788597821 0.43756473770026916 0.5204501210207021 0.6065306597126334 0.6925693242051977 0.7748374288832494 0.8493658165683124 0.9122540768285452 0.9600054412854777 0.9898478033787339 1.0 0.9898478033787338 0.9600054412854776 0.9122540768285451 0.8493658165683124 0.7748374288832494 0.6925693242051977 0.6065306597126334 0.5204501210207021 0.43756473770026916 0.360447788597821 0.2909238070461482 0.23006629899380907 0.17826397958504792 0.1353352832366127 0.10066889977289775 0.0733696513683829 0.05239314106982558 0.036658041953378004 0.025130488985274414 0.01687988414878991 0.011108996538242306 0.007163364470919217 0.004525807866766719 0.0028016380348484843 0.001699279365552657 0.0010098442341126815 0.0005880047627378895 0.00033546262790251185 0.00018751857457483843 0.00010270256461965522 5.5113137010524656e-05 2.89778274286711e-05 1.4928403298804454e-05 7.535251360782923e-06 3.726653172078671e-06 1.805830798460929e-06 8.573773029336787e-07 3.98844636715367e-07 1.8179099965707518e-07 8.118538349144053e-08 3.552386103252793e-08 1.522997974471263e-08 6.397574269770496e-09 2.6331050999142502e-09 1.0618371101937938e-09 4.1955070187358226e-10 1.624231211810119e-10 6.160955882896819e-11 2.289734845645553e-11 8.337947094201161e-12 2.9748832347325257e-12 1.0399622154910407e-12 3.562066673900191e-13 1.195427865933476e-13 3.930805485965858e-14 1.2664165549094176e-14 3.9976839079982126e-15 1.2364517461413708e-15 3.7469917746990634e-16 1.1125643881294095e-16 3.2367146737857757e-17 9.226150318955576e-18 2.5767571091549857e-18 7.051204120964661e-19 1.8905577075969974e-19 4.966534275776451e-20 1.2783617954222807e-20 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 2.1666711874356403e-33 -0.00010270256461974819 -0.00018751857457457714 -0.00033546262790305303 -0.0005880047627374668 -0.0010098442341133016 -0.0016992793655521334 -0.0028016380348489666 -0.004525807866765932 -0.007163364470920133 -0.011108996538241368 -0.016879884148790963 -0.025130488985273244 -0.03665804195337945 -0.052393141069823894 -0.07336965136838487 -0.10066889977289599 -0.13533528323661484 -0.17826397958504586 -0.23006629899381092 -0.2909238070461463 -0.36044778859782295 -0.4375647377002675 -0.5204501210207039 -0.6065306597126316 -0.692569324205199 -0.7748374288832484 -0.8493658165683129 -0.9122540768285452 -0.9600054412854772 -0.9898478033787339 -0.9999999999999998 -0.989847803378734 -0.9600054412854773 -0.9122540768285456 -0.8493658165683119 -0.77483742888325 -0.692569324205196 -0.6065306597126342 -0.5204501210207013 -0.43756473770027016 -0.36044778859782023 -0.29092380704614873 -0.23006629899380893 -0.17826397958504803 -0.13533528323661234 -0.10066889977289822 -0.07336965136838249 -0.05239314106982596 -0.036658041953378004 -0.025130488985274486 -0.0168798841487898 -0.011108996538242379 -0.007163364470919203 -0.004525807866767026 -0.0028016380348483486 -0.0016992793655528666 -0.0010098442341126093 -0.0005880047627381064 -0.0003354626279023046 -0.00018751857457503318 -0.00010270256461948155 -5.511313701084151e-05 -2.8977827428482807e-05 -1.492840329913143e-05 -7.535251360550328e-06 -3.7266531723314633e-06 -1.805830798319928e-06 -8.573773028800223e-07 -3.9884463696009075e-07 -1.8179099984669746e-07 -8.118538329569771e-08 -3.552386127878456e-08 -1.5229979920243373e-08 -6.397574135376856e-09 -2.633105522910678e-09 -1.0618367782550657e-09 -4.195512984041569e-10 -1.6242255781521704e-10 -6.161021889596403e-11 -2.2897057200183176e-11 -8.33821057152514e-12 -2.9745580095154963e-12 -1.0404084189817011e-12 -3.559995218267198e-13 -1.1962410037142686e-13 -3.914164190663356e-14 -1.3090657970339918e-14 -3.832442778776947e-15 -1.2997946275071567e-15 -1.0281557467809818e-16 -1.1072212629837988e-16 -2.242179767451839e-16 -2.740596874716629e-16 1.982826975774268e-16 -2.4547265402142366e-16 7.696226635620134e-17 1.0870526729182805e-17 -5.270100755689163e-17 1.8683545272834908e-17 -1.627095409240092e-17 6.403536026251142e-19 -3.0407217891192303e-19 -7.924391321512499e-22 3.0821181777443354e-22 -1.9868599516498162e-23 8.527779037800367e-24 -3.284305282347064e-27 -1.714783985311744e-27 -3.4238219586344085e-30 -2.290893668848617e-31 1.8840409680968135e-31 -8.387424907806412e-32 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33 6.885199551184368e-33
Instead of printing the output electric field values at the sensor, let's save them to a list and then plot them
ez = [0.0] * SIZE
hy = [0.0] * SIZE
output = [0.0] * maxTime
epsR = [1.0] * SIZE
# do time stepping
for n in range(maxTime):
# update magnetic field
for i in range(SIZE-1):
hy[i] = hy[i] + (ez[i + 1] - ez[i]) / imp0
# update electric field
for i in range(SIZE):
ez[i] = ez[i] + (hy[i] - hy[i - 1]) * imp0 / epsR[i]
# hardwire a source node */
if n < sourceSigma:
ez[0] = math.exp(-(n - sourcePeakTime)**2 / sourceSigma)
else:
## PEC Boundaries
ez[0] = 0.0
ez[(SIZE-1)] = 0.0
output[n] = ez[sensorLocation]
#done with time stepping loop
plt.plot(range(maxTime), output, color='blue', linewidth=1)
plt.xlabel('time')
plt.ylabel('Ez')
plt.show()
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
# set up glass slab
for i in range(150,200):
epsR[i] = 4.7
and rerun the code, storing the output in output2 this time
ez = [0.0] * SIZE
hy = [0.0] * SIZE
output2 = [0.0] * maxTime
# do time stepping
for n in range(maxTime):
# update magnetic field
for i in range(SIZE-1):
hy[i] = hy[i] + (ez[i + 1] - ez[i]) / imp0
# update electric field
for i in range(SIZE):
ez[i] = ez[i] + (hy[i] - hy[i - 1]) * imp0 / epsR[i]
# hardwire a source node */
if n < sourceSigma:
ez[0] = math.exp(-(n - sourcePeakTime)**2 / sourceSigma)
else:
## PEC Boundaries
ez[0] = 0.0
ez[(SIZE-1)] = 0.0
output2[n] = ez[sensorLocation]
#done with time stepping loop
Now let's plot the two outputs, for free space in blue and with a glass slab in green
plt.plot(range(maxTime), output, color='blue', linewidth=1)
plt.plot(range(maxTime), output2, color='green', linewidth=1)
plt.xlabel('time')
plt.ylabel('Ez')
plt.show()
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.
Let's see what's happening by visualizing the Ez field over the simulation, first with no slab
ims = []
# set up a plot
fig, ax = plt.subplots()
ax.axvspan(250, 250, alpha=0.9, color='black') # sensor
ax.set_xlabel('space (mm)')
ax.set_ylabel('Ez')
ez = [0.0] * SIZE
hy = [0.0] * SIZE
# free space
epsR = [1.0] * SIZE
# do time stepping
for n in range(maxTime):
# update magnetic field
for i in range(SIZE-1):
hy[i] = hy[i] + (ez[i + 1] - ez[i]) / imp0
# update electric field
for i in range(SIZE):
ez[i] = ez[i] + (hy[i] - hy[i - 1]) * imp0 / epsR[i]
# capture a snapshot of the ez field at this timestep to the animation
ims.append((plt.plot(range(SIZE), ez, color='blue', linewidth=1)))
# hardwire a source node */
if n < sourceSigma:
ez[0] = math.exp(-(n - sourcePeakTime)**2 / sourceSigma)
else:
## PEC Boundaries
ez[0] = 0.0
ez[(SIZE-1)] = 0.0
#done with time stepping loop
#build and display the animation
im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=5000, blit=True)
HTML(im_ani.to_jshtml())
# # # Set up formatting for the movie files
# Writer = animation.writers['ffmpeg']
# writer = Writer(fps=50, metadata=dict(artist='Me'), bitrate=2000)
# im_ani.save('fdtd_PEC.mp4', writer=writer,dpi=600)
and then with a glass slab
ims = []
# set up a plot
fig, ax = plt.subplots()
ax.axvspan(150, 200, alpha=0.1, color='red') # glass slab
ax.axvspan(250, 250, alpha=0.9, color='black') # sensor
ax.set_xlabel('Space (mm)')
ax.set_ylabel('Ez')
ez = [0.0] * SIZE
hy = [0.0] * SIZE
# free space
epsR = [1.0] * SIZE
# add a glass slab
for i in range(150,200):
epsR[i] = 4.7
# do time stepping
for n in range(maxTime):
# update magnetic field
for i in range(SIZE-1):
hy[i] = hy[i] + (ez[i + 1] - ez[i]) / imp0
# update electric field
for i in range(SIZE):
ez[i] = ez[i] + (hy[i] - hy[i - 1]) * imp0 / epsR[i]
# capture a snapshot of the ez field at this timestep to the animation
ims.append((plt.plot(range(SIZE), ez, color='green', linewidth=1)))
# hardwire a source node */
if n < sourceSigma:
ez[0] = math.exp(-(n - sourcePeakTime)**2 / sourceSigma)
else:
## PEC Boundaries
ez[0] = 0.0
ez[(SIZE-1)] = 0.0
#done with time stepping loop
#build and display the animation
im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=5000, blit=True)
HTML(im_ani.to_jshtml())
# writer = Writer(fps=50, metadata=dict(artist='Me'), bitrate=2000)
# im_ani.save('reflection_PEC.mp4', writer=writer,dpi=600)
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
# Absorbing Boundary Conditions and additive source in a homogeneous medium
ims = []
# set up a plot
fig, ax = plt.subplots()
# ax.axvspan(150, 200, alpha=0.1, color='red') # glass slab
# ax.axvspan(250, 250, alpha=0.9, color='black') # sensor
ax.set_xlabel('Space (mm)')
ax.set_ylabel('Ez')
ez = [0.0] * SIZE
hy = [0.0] * SIZE
# free space
epsR = [1.0] * SIZE
# # add a glass slab
# for i in range(150,200):
# epsR[i] = 4.7
# do time stepping
for n in range(maxTime):
# additive Gausssian source node */
ez[50] += math.exp(-(n - sourcePeakTime)**2 / sourceSigma);
# update magnetic field
for i in range(SIZE-1):
hy[SIZE - 1] = hy[SIZE - 2] # ABC
hy[i] = hy[i] + (ez[i + 1] - ez[i]) / imp0
# update electric field
for i in range(SIZE):
if i ==0:
ez[0] = ez[1] # ABC
else:
ez[i] = ez[i] + (hy[i] - hy[i - 1]) * imp0 / epsR[i]
# capture a snapshot of the ez field at this timestep to the animation
ims.append((plt.plot(range(SIZE), ez, color='blue', linewidth=1.5)))
#done with time stepping loop
#build and display the animation
im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=5000, blit=True)
HTML(im_ani.to_jshtml())
# writer = Writer(fps=50, metadata=dict(artist='Me'), bitrate=2000)
# im_ani.save('reflection_PEC.mp4', writer=writer,dpi=600)
# Absorbing Boundary Conditions and additive source in an inhomogeneous medium
ims = []
# set up a plot
fig, ax = plt.subplots()
ax.axvspan(150, 200, alpha=0.1, color='red') # glass slab
# ax.axvspan(250, 250, alpha=0.9, color='black') # sensor
ax.set_xlabel('Space (mm)')
ax.set_ylabel('Ez')
ez = [0.0] * SIZE
hy = [0.0] * SIZE
# free space
epsR = [1.0] * SIZE
# add a glass slab
for i in range(150,200):
epsR[i] = 4.7
# do time stepping
for n in range(maxTime):
# additive Gausssian source node */
ez[50] += math.exp(-(n - sourcePeakTime)**2 / sourceSigma);
# update magnetic field
for i in range(SIZE-1):
hy[SIZE - 1] = hy[SIZE - 2] # ABC
hy[i] = hy[i] + (ez[i + 1] - ez[i]) / imp0
# update electric field
for i in range(SIZE):
if i ==0:
ez[0] = ez[1] # ABC
else:
ez[i] = ez[i] + (hy[i] - hy[i - 1]) * imp0 / epsR[i]
# capture a snapshot of the ez field at this timestep to the animation
ims.append((plt.plot(range(SIZE), ez, color='blue', linewidth=1.5)))
#done with time stepping loop
#build and display the animation
im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=5000, blit=True)
HTML(im_ani.to_jshtml())
# writer = Writer(fps=50, metadata=dict(artist='Me'), bitrate=2000)
# im_ani.save('reflection_PEC.mp4', writer=writer,dpi=600)
base = 0
transmitted = 0
for n in range(maxTime):
base += output[n]
transmitted += output2[n]
transmission = transmitted/base
print('Transmission =',transmission)
print('Reflectance =', (1 - transmission))
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.
newMaxTime = 501 # simulation time
ez = [0.0] * SIZE
hy = [0.0] * SIZE
output3 = [0.0] * newMaxTime
# free space
epsR = [1.0] * SIZE
# do time stepping
for n in range(newMaxTime):
# update magnetic field
for i in range(SIZE-1):
hy[i] = hy[i] + (ez[i + 1] - ez[i]) / imp0
# update electric field
for i in range(SIZE):
ez[i] = ez[i] + (hy[i] - hy[i - 1]) * imp0 / epsR[i]
# hardwire a source node */
if n < sourceSigma:
ez[0] = math.exp(-(n - sourcePeakTime)**2 / sourceSigma)
else:
ez[0] = 0.0
ez[(SIZE-1)] = 0.0
output3[n] = ez[sensorLocation]
#done with time stepping loop
ez = [0.0] * SIZE
hy = [0.0] * SIZE
output4 = [0.0] * newMaxTime
# add a glass slab
for i in range(150,200):
epsR[i] = 4.7
# do time stepping
for n in range(newMaxTime):
# update magnetic field
for i in range(SIZE-1):
hy[i] = hy[i] + (ez[i + 1] - ez[i]) / imp0
# update electric field
for i in range(SIZE):
ez[i] = ez[i] + (hy[i] - hy[i - 1]) * imp0 / epsR[i]
# hardwire a source node */
if n < sourceSigma:
ez[0] = math.exp(-(n - sourcePeakTime)**2 / sourceSigma)
else:
ez[0] = 0.0
ez[(SIZE-1)] = 0.0
output4[n] = ez[sensorLocation]
#done with time stepping loop
Let's calculate transmission and reflectance again - they should be the same.
base2 = 0
transmitted2 = 0
for n in range(newMaxTime):
base2 += output3[n]
transmitted2 += output4[n]
transmission2 = transmitted2/base2
print('Transmission2 =',transmission2)
print('Reflectance2 =', (1 - transmission2))
Are they the same?
print('Transmission =',transmission)
print('Reflectance =', (1 - transmission))
Not quite. How different are they?
difference = (transmission2 - transmission) / transmission
print(difference)