#!/usr/bin/env python # coding: utf-8 # # Transmittance Spectra for Planewave at Normal Incidence of a Binary Grating # The mode-decomposition feature can also be applied to planewaves in homogeneous media with scalar permittivity/permeability (i.e., no anisotropy). This will be demonstrated in this example to compute the diffraction spectrum of a binary phase [grating](https://en.wikipedia.org/wiki/Diffraction_grating). To compute the diffraction spectrum for a finite-length structure, see [Tutorials/Near to Far Field Spectra/Diffraction Spectrum of a Finite Binary Grating](https://meep.readthedocs.io/en/latest/Python_Tutorials/Near_to_Far_Field_Spectra/#diffraction-spectrum-of-a-finite-binary-grating). # # In this example, we'll define a grating that is periodic in the $y$ direction with periodicity gp and has a rectangular profile of height gh and duty cycle gdc. The grating parameters are gh=0.5 μm, gdc=0.5, and gp=10 μm. There is a semi-infinite substrate of thickness dsub adjacent to the grating. The substrate and grating are glass with a refractive index of 1.5. The surrounding is air/vacuum. Perfectly matched layers (PML) of thickness dpml are used in the $\pm x$ boundaries. # # We'll do our analysis first using an ideal quartz that ignores the dispersion of the material. Then we'll include the dispersion to detect any differences. # # A pulsed planewave with Ez polarization spanning wavelengths of 0.4 to 0.6 μm is normally incident on the grating from the glass substrate. The eigenmode monitor is placed in the air region. We will use mode decomposition to compute the transmittance — the ratio of the power in the +x direction of the diffracted mode relative to that of the incident planewave — for the first ten diffraction orders. # # Two simulations are required: (1) an empty cell of homogeneous glass to obtain the incident power of the source, and (2) the grating structure to obtain the diffraction orders. At the end of the simulation, the wavelength, angle, and transmittance for each diffraction order are computed. # # Initially, we'll import our standard libraries: # In: import meep as mp import math import numpy as np import matplotlib.pyplot as plt # We first need to simulate the empty, homogenous glass (fuzed quartz) with a constant refractive index. # In: fused_quartz = mp.Medium(index=1.5) resolution = 50 # pixels/μm dpml = 1.0 # PML thickness dsub = 3.0 # substrate thickness dpad = 3.0 # padding between grating and PML gp = 10.0 # grating period gh = 0.5 # grating height gdc = 0.5 # grating duty cycle sx = dpml+dsub+gh+dpad+dpml sy = gp cell_size = mp.Vector3(sx,sy,0) pml_layers = [mp.PML(thickness=dpml,direction=mp.X)] wvl_min = 0.4 # min wavelength wvl_max = 0.6 # max wavelength fmin = 1/wvl_max # min frequency fmax = 1/wvl_min # max frequency fcen = 0.5*(fmin+fmax) # center frequency df = fmax-fmin # frequency width src_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub) sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, center=src_pt, size=mp.Vector3(y=sy))] k_point = mp.Vector3(0,0,0) symmetries=[mp.Mirror(mp.Y)] sim = mp.Simulation(resolution=resolution, cell_size=cell_size, boundary_layers=pml_layers, k_point=k_point, default_material=fused_quartz, sources=sources, symmetries=symmetries) nfreq = 21 mon_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad) flux_mon = sim.add_flux(fcen, df, nfreq, mp.FluxRegion(center=mon_pt, size=mp.Vector3(y=sy))) f = plt.figure(dpi=120) sim.plot2D(ax=f.gca()) plt.show() # Now, we'll run the simulation and record the fields. # In: sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-9)) # In: input_flux = mp.get_fluxes(flux_mon) # Next, we'll simulate the actual grating. # In: sim.reset_meep() geometry = [mp.Block(material=fused_quartz, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub))), mp.Block(material=fused_quartz, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh))] sim = mp.Simulation(resolution=resolution, cell_size=cell_size, boundary_layers=pml_layers, geometry=geometry, k_point=k_point, sources=sources, symmetries=symmetries) mode_mon = sim.add_flux(fcen, df, nfreq, mp.FluxRegion(center=mon_pt, size=mp.Vector3(y=sy))) f2 = plt.figure(dpi=120) sim.plot2D(ax=f2.gca()) plt.show() # In: sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-9)) # Now we can compute the diffraction orders as a function of wavelength: # In: freqs = mp.get_eigenmode_freqs(mode_mon) nmode = 10 res = sim.get_eigenmode_coefficients(mode_mon, range(1,nmode+1), eig_parity=mp.ODD_Z+mp.EVEN_Y) coeffs = res.alpha kdom = res.kdom mode_wvl = [] mode_angle = [] mode_tran = [] for nm in range(nmode): for nf in range(nfreq): mode_wvl.append(1/freqs[nf]) mode_angle.append(math.degrees(math.acos(kdom[nm*nfreq+nf].x/freqs[nf]))) tran = abs(coeffs[nm,nf,0])**2/input_flux[nf] mode_tran.append(0.5*tran if nm != 0 else tran) tran_max = round(max(mode_tran),1) # Note the use of the keyword parameter argument eig_parity=mp.ODD_Z+mp.EVEN_Y in the call to get_eigenmode_coefficients. This is important for specifying **non-degenerate** modes in MPB since the k_point is (0,0,0). ODD_Z is for modes with Ez polarization. EVEN_Y is necessary since each diffraction order which is based on a given kx consists of *two* modes: one going in the +y direction and the other in the -y direction. EVEN_Y forces MPB to compute only the +ky + -ky (cosine) mode. As a result, the total transmittance must be halved in this case to obtain the transmittance for the individual +ky or -ky mode. For ODD_Y, MPB will compute the +ky - -ky (sine) mode but this will have zero power because the source is even. If the $y$ parity is left out, MPB will return a random superposition of the cosine and sine modes. Alternatively, in this example an input planewave with Hz instead of Ez polarization can be used which requires eig_parity=mp.EVEN_Z+mp.ODD_Y as well as an odd mirror symmetry plane in *y*. Finally, note the use of add_flux instead of add_mode_monitor when using symmetries. # # The diffraction spectrum is then plotted and shown in the figure below. # In: plt.figure(dpi=200) plt.pcolormesh(np.reshape(mode_wvl,(nmode,nfreq)), np.reshape(mode_angle,(nmode,nfreq)), np.reshape(mode_tran,(nmode,nfreq)), cmap='Blues', shading='flat', vmin=0, vmax=tran_max) plt.axis([min(mode_wvl), max(mode_wvl), min(mode_angle), max(mode_angle)]) plt.xlabel("wavelength (μm)") plt.ylabel("diffraction angle (degrees)") plt.xticks([t for t in np.linspace(wvl_min,wvl_max,3)]) plt.yticks([t for t in range(0,35,5)]) plt.title("transmittance of diffraction orders") cbar = plt.colorbar() cbar.set_ticks([t for t in np.arange(0,tran_max+0.1,0.1)]) cbar.set_ticklabels(["{:.1f}".format(t) for t in np.arange(0,tran_max+0.1,0.1)]) # Now let's try using a dispersive medium instead of an ideal, constant refractive index: # In: from meep.materials import fused_quartz resolution = 50 # pixels/μm dpml = 1.0 # PML thickness dsub = 3.0 # substrate thickness dpad = 3.0 # padding between grating and PML gp = 10.0 # grating period gh = 0.5 # grating height gdc = 0.5 # grating duty cycle sx = dpml+dsub+gh+dpad+dpml sy = gp cell_size = mp.Vector3(sx,sy,0) pml_layers = [mp.PML(thickness=dpml,direction=mp.X)] wvl_min = 0.4 # min wavelength wvl_max = 0.6 # max wavelength fmin = 1/wvl_max # min frequency fmax = 1/wvl_min # max frequency fcen = 0.5*(fmin+fmax) # center frequency df = fmax-fmin # frequency width src_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub) sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez, center=src_pt, size=mp.Vector3(y=sy))] k_point = mp.Vector3(0,0,0) symmetries=[mp.Mirror(mp.Y)] sim = mp.Simulation(resolution=resolution, cell_size=cell_size, boundary_layers=pml_layers, k_point=k_point, default_material=fused_quartz, sources=sources, symmetries=symmetries) nfreq = 21 mon_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad) flux_mon = sim.add_flux(fcen, df, nfreq, mp.FluxRegion(center=mon_pt, size=mp.Vector3(y=sy))) sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-9)) input_flux = mp.get_fluxes(flux_mon) sim.reset_meep() geometry = [mp.Block(material=fused_quartz, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub))), mp.Block(material=fused_quartz, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh))] sim = mp.Simulation(resolution=resolution, cell_size=cell_size, boundary_layers=pml_layers, geometry=geometry, k_point=k_point, sources=sources, symmetries=symmetries) mode_mon = sim.add_flux(fcen, df, nfreq, mp.FluxRegion(center=mon_pt, size=mp.Vector3(y=sy))) sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mon_pt, 1e-9)) freqs = mp.get_eigenmode_freqs(mode_mon) nmode = 10 res = sim.get_eigenmode_coefficients(mode_mon, range(1,nmode+1), eig_parity=mp.ODD_Z+mp.EVEN_Y) coeffs = res.alpha kdom = res.kdom mode_wvl = [] mode_angle = [] mode_tran = [] for nm in range(nmode): for nf in range(nfreq): mode_wvl.append(1/freqs[nf]) mode_angle.append(math.degrees(math.acos(kdom[nm*nfreq+nf].x/freqs[nf]))) tran = abs(coeffs[nm,nf,0])**2/input_flux[nf] mode_tran.append(0.5*tran if nm != 0 else tran) tran_max = round(max(mode_tran),1) plt.figure(dpi=200) plt.pcolormesh(np.reshape(mode_wvl,(nmode,nfreq)), np.reshape(mode_angle,(nmode,nfreq)), np.reshape(mode_tran,(nmode,nfreq)), cmap='Blues', shading='flat', vmin=0, vmax=tran_max) plt.axis([min(mode_wvl), max(mode_wvl), min(mode_angle), max(mode_angle)]) plt.xlabel("wavelength (μm)") plt.ylabel("diffraction angle (degrees)") plt.xticks([t for t in np.linspace(wvl_min,wvl_max,3)]) plt.yticks([t for t in range(0,35,5)]) plt.title("transmittance of diffraction orders") cbar = plt.colorbar() cbar.set_ticks([t for t in np.arange(0,tran_max+0.1,0.1)]) cbar.set_ticklabels(["{:.1f}".format(t) for t in np.arange(0,tran_max+0.1,0.1)]) # Since fuzed quartz isn't highly dispersive in this range, we don't see much of a difference between the two plots.