Computational Electromagnetics with MEEP Part 2: 1D MEEP
Lesson 7: Angular-Dependent Reflection
Everything up to this point has been at normal incidence, where the plane wave is hitting the interface head-on. But often we are interested in what happens at an angle - when we tilt the plane wave. FDTD simulators in general use what are called "Bloch-Periodic boundary conditoins". I don't assume you know what those are and won't be diving into them here, but suffice to say it's the easiest way to tilt planewaves :). It also lets us do very cool things with periodic structures (which we might get into later).
Great, but what does this mean?
In practice, this means a couple of things. First, we specify the angle of a plane wave not in the Source() object as you might expect, but by passing our desired wave vector \(\overrightarrow{k}\) directly in to the simulation. In case you've forgotten, \(\overrightarrow{k}\) is just the direction in which the wave is propagating multiplied by \(|\overrightarrow{k}|=\frac{2\pi}{\lambda}\), (the spatial frequency). At normal incidence, with the wave propagating along the z-direction in free space, \(\overrightarrow{k}\) is just \(\hat{z}*\frac{2\pi}{\lambda_0}\), where \(\hat{z}\) is the unit vector in the z-direction.
Second, it means the fields will be stored in memory in phasor notation - they will become complex. This doesn't matter except when we want to visualize the fields - as EE's we know what to do here - just take the real part. Internally, MEEP will make sure to deal with the complex numbers appropriately so that it correctly calculates things like power.
Enough jabber, let's get started
We are going to use the code from the previous lesson, but changing the material back from Si to have a refractive index of 2 (and the code will be below as usual). Let's also change the number of frequencies so that we only grab one frequency from our Fourier-Transformed fields. This is both so we can sanity-check our simulation, and also because it turns out doing multi-frequency multi-wavelength simulations has some unexpected subtleties in FDTD, which we will deal with once we make sure everything is working at a single frequency.
numberFrequencies = 1
Let's also going to drop the resolution back down to 32 so that our calculations are faster initially, because we don't need the high resolution from the previous lesson any more.
resolution = 32
Since we have to deal with angles now, we need to be more careful about our coordinates, and actually pay attention to things like polarization. So far, we have been working at normal incidence, and our electric field is has been polarized in the "x" direction:
(note: don't take the length on the drawing of the E and k vectors literally, they have different units than length, it's just meant to signify direction). If we want our planewave to come in at an angle, call it \(\theta\), things will now look something like this:
If we zoom in on our k-vector, we see that it has two components: \(k_x\) and \(k_z\). We can see visually that \(k_z\) is just equal to \(|\overrightarrow{k}|cos\theta\) and \(k_x\) is equal to \(|\overrightarrow{k}|sin\theta\). This makes sense because \(\overrightarrow{k}\) is completely along z when \(\theta=0\).
Similarly, the electric field will now have x and z components, instead of just an x-component. We should make sure to verify this in the simulation! Now, all we have to do is figure out how to enter our desired value of \(\overrightarrow{k}\) into MEEP. Remember from a previous lesson that MEEP uses normalized units, where a length we enter into MEEP \(x_{MEEP}\) is equal to a real length \(x\) divided by our chosen characteristic length \(a\).
\(x_{MEEP} = \frac{x}{a} \)
To figure out what to plug in to MEEP, we can start with the definition of (phase) velocity of an EM wave in free space, which we know is equal to the speed of light \(c\):
\(c = \frac{\omega}{|\overrightarrow{k}|} = \frac{2 \pi f}{|\overrightarrow{k}|}\)
If we plug in MEEP's normalized frequency \(f_{MEEP} \) for \(f\) we can rearrange the equation in terms of \(|\overrightarrow{k}|\):
\(|\overrightarrow{k}|=\frac{2\pi}{a}*f_{MEEP}\)
Now, since \(\overrightarrow{k}\) has units of inverse length, and MEEP normalizes everything to make it unitless, the all we need to do to get rid of the units of \(\overrightarrow{k}\) is multiply by our characteristic length \(a\), and we are left with what MEEP stores internally as the wavevector's length, \(2\pi * f_{MEEP}\). However, the designers chose to take care of the factor of \(2\pi\) themselves, so we only need to enter \(f_{MEEP}\):
\(|\overrightarrow{k}_{MEEP}|=f_{MEEP}\)
That's pretty darn simple! Maybe that's why they did it. 🤔
Suppose your frequency \(f_{MEEP}\) is 2. What is the magnitude of the k-vector you need to enter into MEEP?
Again, easier in practice
If all this talk of k-vectors has got you confused, fortunately this is another case where things are much easier in practice than in theory. We need only to add a couple lines of code, which define the k-Vector and scale it so it has the correct length:
theta = np.radians(30) kVector = mp.Vector3(np.sin(theta), 0, np.cos(theta)).scale(frequency)
Here I've entered a value of \(\theta\) of 30 degrees, and added a variable called 'kVector', which is just the \(\overrightarrow{k}\) we have been talking about before. We also need to add an argument in declaring our 'simulation' variable, changing this code:
simulation = mp.Simulation(
cell_size=cellSize,
sources=sources,
resolution=resolution,
boundary_layers=pmlLayers)
To this:
simulation = mp.Simulation(
cell_size=cellSize,
sources=sources,
resolution=resolution,
boundary_layers=pmlLayers,
k_point=kVector)
Make sure to also add it to the second simulation (because remember we are running two simulations), by changing this code:
simulation = mp.Simulation(
cell_size=cellSize,
sources=sources,
resolution=resolution,
boundary_layers=pmlLayers,
geometry=geometry)
To this:
simulation = mp.Simulation(
cell_size=cellSize,
sources=sources,
resolution=resolution,
boundary_layers=pmlLayers,
geometry=geometry,
k_point=kVector)
I would also take the real part of the fieldEx variable before you pass it into fieldData - the code should work as-written but python might complain. Now, let's run it, and see what happens:
You might notice a couple interesting things - first, the pulses are actually smaller than they were at normal incidence. Their amplitude is not as large. Does this make sense?. Well, yes, because our source has the same amplitude, but now the field it generates is split - part along the x direction (which is what we are seeing) and part along the z direction. You might also notice that the crests of the pulse appear to be moving faster than the pulse itself. This has to do with the difference between phase velocity, or how fast the 'peaks' of the pulse move, and group velocity, how fast the pulse itself moves. If you don't look along the direction the pulses are propagating, but along a different direction (as we did, since we are taking a slice along the z-axis, and our pulse is propatating at an angle of 30 degrees). It's hard to notice at this shallow angle - I recommend you try it at steeper angles - but you will also find the pulses at an angle are moving slower. All these phenomena can be explained in terms of phase and group velocity along the z direction:
\(v_{phase, z} = \frac{\omega}{k_z} = c / cos\theta \)
\(v_{group, z} = \frac{d\omega}{dk_z} = c * cos\theta \)
As the angle gets steeper, the group velocity falls and the phase velocity increases. In other words, the pulse slows down, but the crests within the pulse move faster. This is weird, and it takes some getting used to. I recommend you play around with the simulation parameters and see what happens.
You might also notice that the tilted pulse is shorter than the non-tilted one, and the crests are spaced further apart. This is a consequence of our viewing angle - because we are 'projecting' our pulse onto the z-axis, the pulse looks smaller. And since we are tilting a plane wave, and slicing through that tilted wave, the crests are further apart. You can observe these things more dramatically in free space: I ran a couple of simulations, one at 40 degrees incidence, and one at normal incidence, and ran them both for 20 MEEP units.
Normal Incidence
40 degrees
Give me the data
This is all fine and good, but what about some numbers? At an angle of incidence of 30 degrees, we can use Fresnel's equations along with Snell's Law to predict what the reflection / transmission coefficients should be (since this is polarized in the plane of incidence, it is p-polarized light):
\( R_p = |\frac{n_1 * cos\theta_2 - n_2 * cos\theta_1}{n_1 * cos\theta_2 + n_2 * cos\theta_1}|^2 \) (Fresnel's equation)
\( n_1 sin\theta_1 = n_2 sin\theta_2 \) (Snell's Law)
Full Code - Angular Reflection
Grab the code here as well.
import meep as mp
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
resolution = 64
frequency = 2.0
frequencyWidth = 1
numberFrequencies = 1
pmlThickness = 2.0
animationTimestepDuration = 0.05
materialThickness = 2.5
length = 2 * materialThickness + 2 * pmlThickness
powerDecayTarget = 1e-9
transmissionMonitorLocation = mp.Vector3(0, 0, materialThickness/4)
theta = np.radians(30)
kVector = mp.Vector3(np.sin(theta), 0, np.cos(theta)).scale(frequency)
cellSize = mp.Vector3(0, 0, length)
sources = [mp.Source(
mp.GaussianSource(frequency=frequency,fwidth=frequencyWidth),
component=mp.Ex,
center=mp.Vector3(0, 0, -materialThickness/2),
)]
geometry = [mp.Block(
mp.Vector3(mp.inf, mp.inf, materialThickness + pmlThickness),
center=mp.Vector3(0, 0, materialThickness/2 + pmlThickness/2),
material=mp.Medium(index=2))]
pmlLayers = [mp.PML(1.0)]
simulation = mp.Simulation(
cell_size=cellSize,
sources=sources,
resolution=resolution,
boundary_layers=pmlLayers,
k_point=kVector)
incidentRegion = mp.FluxRegion(center=mp.Vector3(0, 0, -materialThickness/4))
incidentFluxMonitor = simulation.add_flux(frequency, frequencyWidth, numberFrequencies, incidentRegion)
simulation.run(until_after_sources=mp.stop_when_fields_decayed(20, mp.Ex,
transmissionMonitorLocation, powerDecayTarget))
fieldEx = np.real(simulation.get_array(center=mp.Vector3(0, 0, 0), size=cellSize, component=mp.Ex))
fieldData = np.zeros(len(fieldEx))
incidentFluxToSubtract = simulation.get_flux_data(incidentFluxMonitor)
simulation = mp.Simulation(
cell_size=cellSize,
sources=sources,
resolution=resolution,
boundary_layers=pmlLayers,
geometry=geometry,
k_point=kVector)
transmissionRegion = mp.FluxRegion(center=transmissionMonitorLocation)
transmissionFluxMonitor = simulation.add_flux(frequency, frequencyWidth, numberFrequencies, transmissionRegion)
reflectionRegion = incidentRegion
reflectionFluxMonitor = simulation.add_flux(frequency, frequencyWidth, numberFrequencies, reflectionRegion)
simulation.load_minus_flux_data(reflectionFluxMonitor, incidentFluxToSubtract)
def updateField(sim):
global fieldData
fieldEx = np.real(sim.get_array(center=mp.Vector3(0, 0, 0), size=cellSize, component=mp.Ex))
fieldData = np.vstack((fieldData, fieldEx))
simulation.run(mp.at_every(animationTimestepDuration, updateField),
until_after_sources=mp.stop_when_fields_decayed(20, mp.Ex,
transmissionMonitorLocation, powerDecayTarget))
incidentFlux = np.array(mp.get_fluxes(incidentFluxMonitor))
transmittedFlux = np.array(mp.get_fluxes(transmissionFluxMonitor))
reflectedFlux = np.array(mp.get_fluxes(reflectionFluxMonitor))
R = -reflectedFlux / incidentFlux
T = transmittedFlux / incidentFlux
print(T)
print(R)
print(R + T)
#frequencies = np.array(mp.get_flux_freqs(reflectionFluxMonitor))
#reflectionSpectraFigure = plt.figure()
#reflectionSpectraAxes = plt.axes(xlim=(frequency-frequencyWidth/2, frequency+frequencyWidth/2),ylim=(0, 1))
#reflectionLine, = reflectionSpectraAxes.plot(frequencies, R, lw=2)
#reflectionSpectraAxes.set_xlabel('frequency (a / \u03BB0)')
#reflectionSpectraAxes.set_ylabel('R')
#plt.show()
fig = plt.figure()
ax = plt.axes(xlim=(-length/2,length/2),ylim=(-1,1))
line, = ax.plot([], [], lw=2)
xData = np.linspace(-length/2, length/2, fieldData.shape[1])
def init():
line.set_data([],[])
return line,
def animate(i):
line.set_data(xData, fieldData[i])
return line,
numberAnimationTimesteps = fieldData.shape[0]
fieldAnimation = animation.FuncAnimation(fig, animate, init_func=init,
frames=numberAnimationTimesteps, interval=20, blit=True)
#fieldAnimation.save('basic_animation.mp4', fps=30, extra_args=['-vcodec', 'libx264'])
plt.show()
Full Code - Angular Propagation
Grab the code here as well.
import meep as mp
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
resolution = 32
frequency = 2.0
frequencyWidth = 1
numberFrequencies = 1
pmlThickness = 2.0
animationTimestepDuration = 0.05
materialThickness = 5
length = 2 * materialThickness + 2 * pmlThickness
powerDecayTarget = 1e-9
transmissionMonitorLocation = mp.Vector3(0, 0, materialThickness/4)
theta = np.radians(0)
kVector = mp.Vector3(np.sin(theta), 0, np.cos(theta)).scale(frequency)
endTime=20
cellSize = mp.Vector3(0, 0, length)
sources = [mp.Source(
mp.GaussianSource(frequency=frequency, fwidth=frequencyWidth),
component=mp.Ex,
center=mp.Vector3(0, 0, -materialThickness/2),
)]
pmlLayers = [mp.PML(1.0)]
simulation = mp.Simulation(
cell_size=cellSize,
sources=sources,
resolution=resolution,
boundary_layers=pmlLayers,
k_point=kVector)
simulation.run(until=0)
fieldEx = np.real(simulation.get_array(center=mp.Vector3(0, 0, 0), size=cellSize, component=mp.Ex))
fieldData = np.zeros(len(fieldEx))
def updateField(sim):
global fieldData
fieldEx = np.real(sim.get_array(center=mp.Vector3(0, 0, 0), size=cellSize, component=mp.Ex))
fieldData = np.vstack((fieldData, fieldEx))
simulation.run(mp.at_every(animationTimestepDuration, updateField),until=endTime)
fig = plt.figure()
ax = plt.axes(xlim=(-length/2,length/2),ylim=(-1,1))
line, = ax.plot([], [], lw=2)
xData = np.linspace(-length/2, length/2, fieldData.shape[1])
def init():
line.set_data([],[])
return line,
def animate(i):
line.set_data(xData, fieldData[i])
return line,
numberAnimationTimesteps = fieldData.shape[0]
fieldAnimation = animation.FuncAnimation(fig, animate, init_func=init,
frames=numberAnimationTimesteps, interval=20, blit=True)
!fieldAnimation.save('basic_animation.mp4', fps=30, extra_args=['-vcodec', 'libx264'])
plt.show()
If you found this content helpful, it would mean a lot to me if you would support me on Patreon. Help keep this content ad-free, get access to my Discord server, exclusive content, and receive my personal thanks for as little as $2. :)
Become a Patron!