add phase map demonstration to binary grating tutorial example (#639)

* add phase map demonstration to binary grating tutorial example

* tweaks to text

* Update Mode_Decomposition.md

doc/docs/Python_Tutorials/Mode_Decomposition.md

@@ -306,6 +306,140 @@ In the limit where the grating periodicity is much larger than the wavelength an
 
 To convert the diffraction efficiency into transmittance in the *x* direction (in order to be able to compare the scalar-theory results with those from Meep), the diffraction efficiency must be multiplied by the Fresnel transmittance from air to glass and by the cosine of the diffraction angle. We compare the analytic and simulated results at a wavelength of 0.5 μm for diffraction orders 1, 3, 5, and 7. The analytic results are 0.3886, 0.0427, 0.0151, and 0.0074. The Meep results are 0.3891, 0.04287, 0.0152, and 0.0076. This corresponds to relative errors of approximately 1.3%, 0.4%, 0.8%, and 2.1% which indicates good agreement.
 
+We can also use the complex mode coefficients to compute the phase of the diffraction orders. This can be used to generate a phase map of the binary grating as a function of its geometric parameters. Phase maps are important for the design of subwavelength phase shifters such as those used in metalenses. In this demonstration, which is adapted from the previous example, we compute the transmittance spectra and phase map of the zeroth diffraction order (at 0°) for an E<sub>z</sub>-polarized planewave pulse spanning wavelengths of 0.4 to 0.6 μm which is normally incident on a binary grating with a periodicity of 0.35 μm and height of 0.6 μm. The duty cycle of the grating is varied from 0.1 to 0.9 in separate runs.
+
+Note that it is generally only the relative phase (the phase difference) between different structures that is useful. The overall mode coefficient α is multiplied by a complex number given by the source amplitude, as well as an arbitrary (but deterministic) phase choice by our mode solver (MPB)—but as long as you keep the current source fixed as you vary the parameters of the structure, the relative phases are meaningful. More generally, for metasurface design it is preferable to use the overall complex mode coefficient α (both phase and amplitude), as we describe in [this 2018 paper](https://www.osapublishing.org/oe/abstract.cfm?uri=oe-26-26-33732): you can use α to optimize for a desired wavefront, as in equation (3) in section 3.1, or more generally for some arbitrary function of the far field as in section 3.2.
+
+The simulation script is in [examples/binary_grating_phasemap.py](https://github.com/stevengj/meep/blob/master/python/examples/binary_grating_phasemap.py).
+
+```py
+import meep as mp
+import numpy as np
+import matplotlib.pyplot as plt
+import numpy.matlib
+import argparse
+
+resolution = 60         # pixels/μm
+
+dpml = 1.0              # PML thickness
+dsub = 3.0              # substrate thickness
+dpad = 3.0              # padding between grating and PML
+
+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
+nfreq = 21              # number of frequency bins
+
+k_point = mp.Vector3(0,0,0)
+
+glass = mp.Medium(index=1.5)
+
+def grating(gp,gh,gdc,oddz):
+  sx = dpml+dsub+gh+dpad+dpml
+  sy = gp
+
+  cell_size = mp.Vector3(sx,sy,0)
+  pml_layers = [mp.PML(thickness=dpml,direction=mp.X)]
+
+  src_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub,0,0)
+  sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez if oddz else mp.Hz, center=src_pt, size=mp.Vector3(0,sy,0))]
+
+  symmetries=[mp.Mirror(mp.Y, phase=+1 if oddz else -1)]
+
+  sim = mp.Simulation(resolution=resolution,
+                      cell_size=cell_size,
+                      boundary_layers=pml_layers,
+                      k_point=k_point,
+                      default_material=glass,
+                      sources=sources,
+                      symmetries=symmetries)
+
+  mon_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad,0,0)
+  flux_mon = sim.add_flux(fcen, df, nfreq, mp.FluxRegion(center=mon_pt, size=mp.Vector3(0,sy,0)))
+
+  sim.run(until_after_sources=100)
+
+  input_flux = mp.get_fluxes(flux_mon)
+
+  sim.reset_meep()
+
+  geometry = [mp.Block(material=glass, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub),0,0)),
+              mp.Block(material=glass, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh,0,0))]
+
+  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(0,sy,0)))
+
+  sim.run(until_after_sources=300)
+
+  freqs = mp.get_eigenmode_freqs(mode_mon)
+  res = sim.get_eigenmode_coefficients(mode_mon, [1], eig_parity=mp.ODD_Z+mp.EVEN_Y if oddz else mp.EVEN_Z+mp.ODD_Y)
+  coeffs = res.alpha
+
+  mode_wvl = [1/freqs[nf] for nf in range(nfreq)]
+  mode_tran = [abs(coeffs[0,nf,0])**2/input_flux[nf] for nf in range(nfreq)]
+  mode_phase = [np.angle(coeffs[0,nf,0]) for nf in range(nfreq)]
+
+  return mode_wvl, mode_tran, mode_phase
+
+if __name__ == '__main__':
+  parser = argparse.ArgumentParser()
+  parser.add_argument('-gp', type=float, default=0.35, help='grating periodicity (default: 0.35 μm)')
+  parser.add_argument('-gh', type=float, default=0.6, help='grating height (default: 0.6 μm)')
+  parser.add_argument('-oddz', action='store_true', default=False, help='oddz? (default: False)')
+  args = parser.parse_args()
+
+  gdc = np.arange(0.1,1.0,0.1)
+  mode_tran = np.empty((gdc.size,nfreq))
+  mode_phase = np.empty((gdc.size,nfreq))
+  for n in range(gdc.size):
+    mode_wvl, mode_tran[n,:], mode_phase[n,:] = grating(args.gp,args.gh,gdc[n],args.oddz)
+
+  plt.figure(dpi=150)
+
+  plt.subplot(1,2,1)
+  plt.pcolormesh(mode_wvl, gdc, mode_tran, cmap='hot_r', shading='gouraud', vmin=0, vmax=mode_tran.max())
+  plt.axis([wvl_min, wvl_max, gdc[0], gdc[-1]])
+  plt.xlabel("wavelength (μm)")
+  plt.xticks([t for t in np.arange(wvl_min,wvl_max+0.1,0.1)])
+  plt.ylabel("grating duty cycle")
+  plt.yticks([t for t in np.arange(gdc[0],gdc[-1]+0.1,0.1)])
+  plt.title("transmittance")
+  cbar = plt.colorbar()
+  cbar.set_ticks([t for t in np.arange(0,1.2,0.2)])
+  cbar.set_ticklabels(["{:.1f}".format(t) for t in np.arange(0,1.2,0.2)])
+
+  plt.subplot(1,2,2)
+  plt.pcolormesh(mode_wvl, gdc, mode_phase, cmap='RdBu', shading='gouraud', vmin=mode_phase.min(), vmax=mode_phase.max())
+  plt.axis([wvl_min, wvl_max, gdc[0], gdc[-1]])
+  plt.xlabel("wavelength (μm)")
+  plt.xticks([t for t in np.arange(wvl_min,wvl_max+0.1,0.1)])
+  plt.ylabel("grating duty cycle")
+  plt.yticks([t for t in np.arange(gdc[0],gdc[-1]+0.1,0.1)])
+  plt.title("phase (radians)")
+  cbar = plt.colorbar()
+  cbar.set_ticks([t for t in range(-3,4)])
+  cbar.set_ticklabels(["{:.1f}".format(t) for t in range(-3,4)])
+
+  plt.tight_layout()
+  plt.show()
+```
+
+Note that the phase of the zeroth diffraction order is simply the angle of its complex mode coefficient. The script is run from the shell terminal using: `python binary_grating_phasemap.py -gp 0.35 -gh 0.6 -oddz`. The figure shown below is produced of the transmittance spectra (left) and phase map (right). The transmittance is nearly unity over most of the parameter space mainly because of the subwavlength dimensions of the grating. The phase variation spans the full range of -2π to +2π but is weak as a result of the the relatively low index of the glass grating. Higher-index materials such as [titanium dioxide](https://en.wikipedia.org/wiki/Titanium_dioxide#Thin_films) (TiO<sub>2</sub>) generally provide more control over the phase.
+
+<center>
+![](../images/grating_phasemap.png)
+</center>
+
 ### Reflectance and Transmittance Spectra for Planewave at Oblique Incidence
 
 As an additional demonstration of the mode-decomposition feature, the reflectance and transmittance of all diffracted orders for any grating with no material absorption and a planewave source incident at any arbitrary angle and wavelength must necessarily sum to unity. Also, the total reflectance and transmittance must be equivalent to values computed using the Poynting flux. This demonstration is somewhat similar to the [single-mode waveguide example](#reflectance-of-a-waveguide-taper).

python/examples/binary_grating_phasemap.py

@@ -0,0 +1,121 @@
+# -*- coding: utf-8 -*-
+
+import meep as mp
+import numpy as np
+import matplotlib.pyplot as plt
+import numpy.matlib
+import argparse
+
+resolution = 60         # pixels/μm
+
+dpml = 1.0              # PML thickness
+dsub = 3.0              # substrate thickness
+dpad = 3.0              # padding between grating and PML
+
+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
+nfreq = 21              # number of frequency bins
+
+k_point = mp.Vector3(0,0,0)
+
+glass = mp.Medium(index=1.5)
+
+def grating(gp,gh,gdc,oddz):
+  sx = dpml+dsub+gh+dpad+dpml
+  sy = gp
+
+  cell_size = mp.Vector3(sx,sy,0)
+  pml_layers = [mp.PML(thickness=dpml,direction=mp.X)]
+
+  src_pt = mp.Vector3(-0.5*sx+dpml+0.5*dsub,0,0)
+  sources = [mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ez if oddz else mp.Hz, center=src_pt, size=mp.Vector3(0,sy,0))]
+
+  symmetries=[mp.Mirror(mp.Y, phase=+1 if oddz else -1)]
+
+  sim = mp.Simulation(resolution=resolution,
+                      cell_size=cell_size,
+                      boundary_layers=pml_layers,
+                      k_point=k_point,
+                      default_material=glass,
+                      sources=sources,
+                      symmetries=symmetries)
+
+  mon_pt = mp.Vector3(0.5*sx-dpml-0.5*dpad,0,0)
+  flux_mon = sim.add_flux(fcen, df, nfreq, mp.FluxRegion(center=mon_pt, size=mp.Vector3(0,sy,0)))
+
+  sim.run(until_after_sources=100)
+
+  input_flux = mp.get_fluxes(flux_mon)
+
+  sim.reset_meep()
+
+  geometry = [mp.Block(material=glass, size=mp.Vector3(dpml+dsub,mp.inf,mp.inf), center=mp.Vector3(-0.5*sx+0.5*(dpml+dsub),0,0)),
+              mp.Block(material=glass, size=mp.Vector3(gh,gdc*gp,mp.inf), center=mp.Vector3(-0.5*sx+dpml+dsub+0.5*gh,0,0))]
+
+  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(0,sy,0)))
+
+  sim.run(until_after_sources=300)
+
+  freqs = mp.get_eigenmode_freqs(mode_mon)
+  res = sim.get_eigenmode_coefficients(mode_mon, [1], eig_parity=mp.ODD_Z+mp.EVEN_Y if oddz else mp.EVEN_Z+mp.ODD_Y)
+  coeffs = res.alpha
+
+  mode_wvl = [1/freqs[nf] for nf in range(nfreq)]
+  mode_tran = [abs(coeffs[0,nf,0])**2/input_flux[nf] for nf in range(nfreq)]
+  mode_phase = [np.angle(coeffs[0,nf,0]) for nf in range(nfreq)]
+
+  return mode_wvl, mode_tran, mode_phase
+
+if __name__ == '__main__':
+  parser = argparse.ArgumentParser()
+  parser.add_argument('-gp', type=float, default=0.35, help='grating periodicity (default: 0.35 μm)')
+  parser.add_argument('-gh', type=float, default=0.6, help='grating height (default: 0.6 μm)')
+  parser.add_argument('-oddz', action='store_true', default=False, help='oddz? (default: False)')
+  args = parser.parse_args()
+
+  gdc = np.arange(0.1,1.0,0.1)
+  mode_tran = np.empty((gdc.size,nfreq))
+  mode_phase = np.empty((gdc.size,nfreq))
+  for n in range(gdc.size):
+    mode_wvl, mode_tran[n,:], mode_phase[n,:] = grating(args.gp,args.gh,gdc[n],args.oddz)
+
+  plt.figure(dpi=150)
+
+  plt.subplot(1,2,1)
+  plt.pcolormesh(mode_wvl, gdc, mode_tran, cmap='hot_r', shading='gouraud', vmin=0, vmax=mode_tran.max())
+  plt.axis([wvl_min, wvl_max, gdc[0], gdc[-1]])
+  plt.xlabel("wavelength (μm)")
+  plt.xticks([t for t in np.arange(wvl_min,wvl_max+0.1,0.1)])
+  plt.ylabel("grating duty cycle")
+  plt.yticks([t for t in np.arange(gdc[0],gdc[-1]+0.1,0.1)])
+  plt.title("transmittance")
+  cbar = plt.colorbar()
+  cbar.set_ticks([t for t in np.arange(0,1.2,0.2)])
+  cbar.set_ticklabels(["{:.1f}".format(t) for t in np.arange(0,1.2,0.2)])
+
+  plt.subplot(1,2,2)
+  plt.pcolormesh(mode_wvl, gdc, mode_phase, cmap='RdBu', shading='gouraud', vmin=mode_phase.min(), vmax=mode_phase.max())
+  plt.axis([wvl_min, wvl_max, gdc[0], gdc[-1]])
+  plt.xlabel("wavelength (μm)")
+  plt.xticks([t for t in np.arange(wvl_min,wvl_max+0.1,0.1)])
+  plt.ylabel("grating duty cycle")
+  plt.yticks([t for t in np.arange(gdc[0],gdc[-1]+0.1,0.1)])
+  plt.title("phase (radians)")
+  cbar = plt.colorbar()
+  cbar.set_ticks([t for t in range(-3,4)])
+  cbar.set_ticklabels(["{:.1f}".format(t) for t in range(-3,4)])
+
+  plt.tight_layout()
+  plt.show()