TIDY3D
LEARNING CENTER

Lectures

Tutorials

Applications

Using FDTD to Compute a Transmission Spectrum

By Tyler Hughes, Zongfu Yu and Shanhui Fan

In this lecture, we show how to use FDTD to solve a basic EM problem involving the transmission of light though a slab of material.

- How to set up the simulation using periodic boundary conditions and a plane wave source.
- How to use a broadband pulse to compute the normalized transmission over a range of wavelengths.
- How the FDTD results agree with those computed with analytical methods.

Share On:

Visit public projects:

Download .ipynb View Presentation Slides

Additional information: This was updated in Apr 14, 2022

FDTD_Python_Tutorial_02

Tutorial 2: Slab Transmission¶

This tutorial will walk you through computing transmission spectrum of a dielectric slab using Tidy3D.

For more details, refer to the Tidy3d documentation.

Setup¶

First we import Tidy3D and the other packages needed. If it is not installed, please perform the installation and account authentication by following the detailed instructions here.

Introduction / Setup¶

In [1]:

                                        # standard python imports
import numpy as np
import matplotlib.pyplot as plt

import tidy3d as td
from tidy3d import web

                                      

First, let us define some basic parameters setting up the slab and other simulation properties.

In [2]:

                                        # Wavelength and frequency range
freq_range = (200e12, 400e12)
lambda_range = (td.constants.C_0/freq_range[1], td.constants.C_0/freq_range[0])
freq0 = np.sum(freq_range)/2
Nfreq = 301

# frequencies and wavelengths of monitor
monitor_freqs = np.linspace(freq_range[0], freq_range[1], Nfreq)
monitor_lambdas = td.constants.C_0 / monitor_freqs

# central frequency, frequency pulse width and total running time
lambda0 = td.C_0 / freq0
bandwidth = 0.18  # bandwidth of source in units of delta frequency.  0.38 for broadband

freqw  = bandwidth * (freq_range[1] - freq_range[0])
t_stop = 100 / freq0

# Thickness and refractive index of slab
t_slab = 0.5
n_slab = 3.5
mat_slab = td.Medium(permittivity=n_slab**2, name='silicon')

# Grid resolution (cells per um)
# dl = lambda_range[0] / 30 / n_slab
dl = 6e-3

# space between slabs and sources and PML
spacing = 1 * lambda_range[-1]

# simulation size
sim_size = Lx, Ly, Lz = (2.0, 2.0, 2*spacing + t_slab)

                                      

Analytical Method¶

First, let's use the analytical expresion for slab transmission to see what we should expect.

In [3]:

                                        def slab_analytical(d, n, wvl):
    """ computes transmision as a function of slab thickness (d), refractive index (n), and wavelength (wvl). """
    rho = (n-1)/(n+1)
    t = ((1+rho)*(1-rho)*np.exp(-2j*np.pi*n*d/wvl)) / (1 - rho**2*np.exp(-4j*np.pi*n*d/wvl))
    return np.abs(t)**2

transmission_analytical = slab_analytical(t_slab, n_slab, monitor_lambdas)

                                      

In [4]:

                                        plt.figure()
plt.plot(monitor_freqs / 1e12, transmission_analytical, 'k', label='analytical')
plt.xlabel('frequency (THz)')
plt.ylabel('Transmitted')
plt.legend()
plt.show()

                                      

Create Simulation¶

Now we set everything else up (structures, sources, monitors, simulation) to run the example.

First, we define the multilayer stack structure.

In [5]:

                                        slab = td.Structure(
    geometry=td.Box(
        center=(0, 0, -Lz/2 + spacing + t_slab/2),
        size=(td.inf, td.inf, t_slab),
    ),
    medium=mat_slab,
    name='slab',
)

                                      

We must now define the excitation conditions and field monitors. We will excite the slab using a normally incident (along z) planewave, polarized along the x direciton.

In [6]:

                                        # Here we define the planewave source, placed just in advance (towards negative z) of the slab
source = td.PlaneWave(
    source_time = td.GaussianPulse(
        freq0=freq0,
        fwidth=freqw
    ),
    size=(td.inf, td.inf, 0),
    center=(0, 0, -Lz/2+spacing/2),
    direction='+',
    pol_angle=0,
    name='planewave',
)

                                      

Here we define the field monitor, placed just past (towards positive z) of the stack.

In [7]:

                                        # We are interested in measuring the transmitted flux, so we set it to be an oversized plane.
monitor = td.FluxMonitor(
    center = (0, 0, Lz/2 - spacing/2),
    size = (td.inf, td.inf, 0),
    freqs = monitor_freqs,
    name='flux',
)

                                      

Now it is time to define the simulation object.

In [8]:

                                        sim = td.Simulation(
    center = (0, 0, 0),
    size = sim_size,
    grid_spec = td.GridSpec.uniform(dl=dl),
    structures = [slab],
    sources = [source],
    monitors = [monitor],
    run_time = t_stop,
    boundary_spec = td.BoundarySpec.pml(z=True),
    normalize_index = None,
)

                                      

Plot The Structure¶

Let's now plot the permittivity profile to confirm that the structure was defined correctly.

First we use the Simulation.plot() method to plot the materials only, which assigns a different color to each slab without knowledge of the material properties.

In [9]:

                                        sim.plot(x=0)
plt.show()

Next, we use Simulation.plot_eps() to vizualize the permittivity of the stack. However, because the stack contains dispersive materials, we need to specify the freq of interest as an argument to the plotting tool. Here we show the permittivity at the lowest and highest frequences in the range of interest. Note that in this case, the real part of the permittivity (being plotted) only changes slightly between the two frequencies on the dispersive material. However, for other materials with more dispersion, the effect can be much more prominent.

In [10]:

                                        # plot the permittivity 
freqs_plot = (freq_range[0], freq_range[1])
fig, ax = plt.subplots(1, tight_layout=True, figsize=(8, 4))
sim.plot_eps(x=0, freq=None, ax=ax)
plt.show()

                                      

In [11]:

                                        times = np.arange(0, sim.run_time, sim.dt)
amps = np.real(sim.sources[0].source_time.amp_time(times))

We can also take a look at the source to make sure it's defined correctly over our frequency range of interst.

In [12]:

                                        # Check probe and source
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(8, 4))
ax1 = sim.sources[0].source_time.plot(times=np.linspace(0, sim.run_time, 1001), ax=ax1)
ax1.set_xlim(0, sim.run_time)
ax1.legend(('source amplitude',))
ax2 = sim.sources[0].source_time.plot_spectrum(times=np.linspace(0, sim.run_time, 1001), val='abs', ax=ax2)
fill_max = 30e-16 # 10e-16
ymax = 50e-16 # 20e-16
ax2.fill_between(freq_range, [-0e-16, -0e-16], [fill_max, fill_max], alpha=0.4, color='g')
ax2.legend(('source spectrum', 'measurement'))
ax2.set_ylim(-1e-16, ymax)
plt.show()

                                      

Define Normalization Simulation¶

To compare the results with and without the slab, let's make a copy of our simulation and remove the slab. This will be our reference for normalization.

In [13]:

                                        sim0 = sim.copy(update={'structures':[]})

                                      

Run the simulations¶

We will submit both simulations to run.

In [14]:

                                        sim_data0 = web.run(sim0, task_name='lecture02_slab_narrowband_normalization', path=f'data/data0_narrowband.hdf5', verbose=True)
sim_data = web.run(sim, task_name='lecture02_slab_narrowband_transmission', path=f'data/data_narrowband.hdf5', verbose=True)

[10:22:13] Created task 'lecture02_slab_narrowband_normalization' with task_id  webapi.py:139
           'fdve-06a38fb7-83b1-4621-ada3-27e4dab4c490v1'.

[10:22:17] status = queued                                                      webapi.py:269

[10:22:33] status = preprocess                                                  webapi.py:263

[10:22:37] Maximum FlexCredit cost: 0.603. Use 'web.real_cost(task_id)' to get  webapi.py:286
           the billed FlexCredit cost after a simulation run.

           starting up solver                                                   webapi.py:290

           running solver                                                       webapi.py:300

[10:23:27] early shutoff detected, exiting.                                     webapi.py:313

           status = postprocess                                                 webapi.py:330

[10:23:29] status = success                                                     webapi.py:337

[10:23:31] loading SimulationData from data/data0_narrowband.hdf5               webapi.py:509

[10:23:31] Created task 'lecture02_slab_narrowband_transmission' with task_id   webapi.py:139
           'fdve-4074db4c-e7c9-4b2e-9c7e-f273d8b7dbdfv1'.

[10:23:32] status = queued                                                      webapi.py:269

[10:23:45] status = preprocess                                                  webapi.py:263

[10:23:53] Maximum FlexCredit cost: 0.603. Use 'web.real_cost(task_id)' to get  webapi.py:286
           the billed FlexCredit cost after a simulation run.

           starting up solver                                                   webapi.py:290

[10:23:54] running solver                                                       webapi.py:300

[10:25:19] early shutoff detected, exiting.                                     webapi.py:313

[10:25:20] status = postprocess                                                 webapi.py:330

[10:25:25] status = success                                                     webapi.py:337

[10:25:26] loading SimulationData from data/data_narrowband.hdf5                webapi.py:509

Postprocess and Plot¶

Once the simulation has completed, we can download the results and load them into the simulation object.

Now, we compute the transmitted flux and plot the transmission spectrum.

In [15]:

                                        # Retrieve the power flux through the monitor plane.
transmission0 = sim_data0['flux'].flux
transmission = sim_data['flux'].flux
transmission_normalized = transmission / transmission0

f, (ax1, ax2) = plt.subplots(2, 1, figsize=(4, 7), tight_layout=True)
transmission0.plot(ax=ax1, label='no slab')
transmission.plot(ax=ax1, label='with slab')
transmission_normalized.plot(ax=ax2)
ax1.legend()
ax1.set_title('raw transmission')
ax2.set_title('normalized transmission')
plt.show()

                                      

We notice that at the edges of the frequency spectrum, the results look a bit strange.

In [16]:

                                        plt.figure()
plt.plot(monitor_freqs/1e12, transmission_analytical, 'k', label='analytical')
plt.plot(monitor_freqs/1e12, transmission_normalized, 'b--', label='FDTD')
plt.xlabel('frequency (THz)')
plt.xlim([200, 400])
plt.ylabel('Transmitted')
plt.legend()
plt.show()

                                      

And also dont match the analytical.

Inreasing Source Bandwdith¶

Those results show some strange behavior at the edge frequencies. As we discuss in the lecture, it is likely due to the source spectrum not having enough power in those regions. So let's increase the source bandwidth and try again.

In [17]:

                                        bandwidth = 0.38

freqw  = bandwidth * (freq_range[1] - freq_range[0])

# make a copy of the original simulation, but change the frequency width of the source
bw_source = source.copy(update={
    'source_time':source.source_time.copy(update={'fwidth':freqw})
    }
)
sim_bw = sim.copy(update={'sources':[bw_source]})

# make the normalization simulation for this
sim_bw0 = sim_bw.copy(update={'structures':[]})

                                      

In [18]:

                                        # Check probe and source, looks now like the source spectrum fully covers the green square
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(8, 4))
ax1 = sim_bw.sources[0].source_time.plot(times=np.linspace(0, sim_bw.run_time, 1001), ax=ax1)
ax1.set_xlim(0, sim_bw.run_time)
ax1.legend(('source amplitude',))
ax2 = sim_bw.sources[0].source_time.plot_spectrum(times=np.linspace(0, sim_bw.run_time, 1001), val='abs', ax=ax2)
fill_max = 15e-16
ymax = 20e-16
ax2.fill_between(freq_range, [-0e-16, -0e-16], [fill_max, fill_max], alpha=0.4, color='g')
ax2.legend(('source spectrum', 'measurement'))
ax2.set_ylim(-1e-16, ymax)
plt.show()

                                      

In [19]:

                                        sim_bw_data0 = web.run(sim_bw0, task_name='lecture02_slab_broadband_normalizatio', path=f'data/data0_broadband.hdf5')
sim_bw_data = web.run(sim_bw, task_name='lecture02_slab_broadband_transmission', path=f'data/data_broadband.hdf5')

[10:25:27] Created task 'lecture02_slab_broadband_normalizatio' with task_id    webapi.py:139
           'fdve-c6b72d5c-8938-4033-9178-6aa8a1545355v1'.

[10:25:28] status = queued                                                      webapi.py:269

[10:25:46] Maximum FlexCredit cost: 0.607. Use 'web.real_cost(task_id)' to get  webapi.py:286
           the billed FlexCredit cost after a simulation run.

           starting up solver                                                   webapi.py:290

           running solver                                                       webapi.py:300

[10:26:14] early shutoff detected, exiting.                                     webapi.py:313

           status = postprocess                                                 webapi.py:330

[10:26:17] status = success                                                     webapi.py:337

[10:26:18] loading SimulationData from data/data0_broadband.hdf5                webapi.py:509

[10:26:18] Created task 'lecture02_slab_broadband_transmission' with task_id    webapi.py:139
           'fdve-097ef95a-0709-4c43-be93-4216b0ff3f4ev1'.

[10:26:19] status = queued                                                      webapi.py:269

[10:26:34] status = preprocess                                                  webapi.py:263

[10:26:42] Maximum FlexCredit cost: 0.607. Use 'web.real_cost(task_id)' to get  webapi.py:286
           the billed FlexCredit cost after a simulation run.

           starting up solver                                                   webapi.py:290

           running solver                                                       webapi.py:300

[10:28:08] early shutoff detected, exiting.                                     webapi.py:313

           status = postprocess                                                 webapi.py:330

[10:28:13] status = success                                                     webapi.py:337

[10:28:15] loading SimulationData from data/data_broadband.hdf5                 webapi.py:509

In [20]:

                                        # Retrieve the power flux through the monitor plane.
transmission0 = sim_bw_data0['flux'].flux
transmission = sim_bw_data['flux'].flux
transmission_normalized = transmission / transmission0

f, (ax1, ax2) = plt.subplots(2, 1, figsize=(4, 7), tight_layout=True)
transmission0.plot(ax=ax1, label='no slab')
transmission.plot(ax=ax1, label='with slab')
transmission_normalized.plot(ax=ax2)
ax1.legend()
ax1.set_title('raw transmission')
ax2.set_title('normalized transmission')
plt.show()

                                      

In [21]:

                                        plt.figure()
plt.plot(monitor_freqs/1e12, transmission_analytical, 'k', label='analytical')
plt.plot(monitor_freqs/1e12, transmission_normalized, 'b--', label='FDTD')
plt.xlabel('frequency (THz)')
plt.xlim([200, 400])
plt.ylabel('Transmitted')
plt.legend()
plt.show()

                                      

These results look much better.

In [ ]:

Download FDTD Computed Transmission Spectrum

TIDY3D LEARNING CENTER