minor tweaks and fixes to docs (#1958)

* minor tweaks and fixes to docs

* add title to plot indicating large errors from finite truncation and numerical dispersion

* mention pymeep-extras Conda package in Installation page

README.md

@@ -3,7 +3,7 @@
 [![CI](https://github.com/NanoComp/meep/actions/workflows/build-ci.yml/badge.svg)](https://github.com/NanoComp/meep/actions/workflows/build-ci.yml)
 [![Sanitizers](https://github.com/NanoComp/meep/actions/workflows/build-san.yml/badge.svg)](https://github.com/NanoComp/meep/actions/workflows/build-san.yml)
 [![Latest Docs](https://readthedocs.org/projects/meep/badge/?version=latest)](http://meep.readthedocs.io/en/latest/)
-![Python versions 3.6–3.9](https://img.shields.io/badge/python-3.6%2C%203.7%2C%203.8%2C%203.9-brightgreen.svg)
+![Python versions 3.7–3.10](https://img.shields.io/badge/python-3.7%2C%203.7%2C%203.8%2C%203.10-brightgreen.svg)
 [![codecov](https://codecov.io/gh/NanoComp/meep/branch/master/graph/badge.svg?token=k88ZuW3795)](https://codecov.io/gh/NanoComp/meep)
 
 **Meep** is a free and open-source software package for [electromagnetics](https://en.wikipedia.org/wiki/Electromagnetism) simulation via the [finite-difference time-domain](https://en.wikipedia.org/wiki/Finite-difference_time-domain_method) (FDTD) method spanning a broad range of applications.

doc/docs/Installation.md

@@ -28,13 +28,13 @@ bash miniconda.sh -b -p <desired_prefix>
 export PATH=<desired_prefix>/bin:$PATH
 ```
 
-Next, we create a Conda environment for PyMeep to isolate it from other Python libraries that may be installed.
+Next, we create a Conda environment for PyMeep (serial version) to isolate it from other Python libraries that may be installed.
 
 ```bash
-conda create -n mp -c conda-forge pymeep
+conda create -n mp -c conda-forge pymeep pymeep-extras
 ```
 
-This creates an environment called "mp" (you can name this anything you like) with PyMeep and all its dependencies. This will default to the version of Python in your Miniconda installation (Python 3 for us since we installed Miniconda3), but if you want to work with Python 2, just add `python=2` to the end of the command.
+This creates an environment called "mp" (you can name this anything you like) with PyMeep and all its dependencies. This will default to the version of Python in your Miniconda installation. The optional "pymeep-extras" package contains scipy, matplotlib, autograd, and other packages necessary if you want to try out the tutorial examples.
 
 Next, we need to activate the environment before we can start using it.
 

doc/docs/Python_Tutorials/Adjoint_Solver.md

@@ -2,7 +2,7 @@
 # Adjoint Solver
 ---
 
-Meep has a density-based adjoint solver for efficiently computing the gradient of an arbitrary function of the mode coefficients (S-parameters), DFT fields, and far fields (using the analytic near-to-far field transformation) with respect to $\varepsilon$ on a discrete spatial grid (a [`MaterialGrid`](../Python_User_Interface.md#materialgrid) class object) at multiple frequencies over a broad bandwidth. Regardless of the number of degrees of freedom for the grid points, just **two** separate timestepping runs are required. The first run is the "forward" calculation to compute the objective function and the DFT fields of the design region. The second run is the "adjoint" calculation to compute the gradient of the objective function with respect to the design variables. The adjoint run involves a special type of source distribution used to compute the DFT fields of the design region. The gradient is computed in post processing using the DFT fields from the forward and adjoint runs. This is all done automatically, and is described in:
+Meep contains an adjoint-solver module for efficiently computing the gradient of an arbitrary function of the mode coefficients (S-parameters), DFT fields, and far fields (using the analytic near-to-far field transformation) with respect to $\varepsilon$ on a discrete spatial grid (a [`MaterialGrid`](../Python_User_Interface.md#materialgrid) class object) at multiple frequencies over a broad bandwidth. Regardless of the number of degrees of freedom for the grid points, just **two** separate timestepping runs are required. The first run is the "forward" calculation to compute the objective function and the DFT fields of the design region. The second run is the "adjoint" calculation to compute the gradient of the objective function with respect to the design variables. The adjoint run involves a special type of source distribution used to compute the DFT fields of the design region. The gradient is computed in post processing using the DFT fields from the forward and adjoint runs. This is all done automatically, and is described in:
 
 - A. M. Hammond, A. Oskooi, M. Chen, Z. Lin, S. G. Johnson, and S. E. Ralph, “[High-performance hybrid time/frequency-domain topology optimization for large-scale photonics inverse design](https://doi.org/10.1364/OE.442074),” *Optics Express*, vol. 30, no. 3, pp. 4467–4491 (2022).
 

doc/docs/Python_Tutorials/Mode_Decomposition.md

@@ -15,7 +15,7 @@ This example involves computing the reflectance of the fundamental mode of a lin
 ![](../images/waveguide-taper.png)
 </center>
 
-The structure, which can be viewed as a [two-port network](https://en.wikipedia.org/wiki/Two-port_network), consists of a single-mode waveguide of width 1 μm (`w1`) at a wavelength of 6.67 μm and coupled to a second waveguide of width 2 μm (`w2`) via a linearly-sloped taper of variable length `Lt`. The material is silicon with $\varepsilon=12$. The taper geometry is defined using a single [`Prism`](../Python_User_Interface.md#prism) object with eight vertices. PML absorbing boundaries surround the entire cell. An eigenmode current source with $E_z$ polarization is used to launch the fundamental mode. The dispersion relation (or "band diagram") of the single-mode waveguide is shown in [Tutorial/Eigenmode Source](Eigenmode_Source.md). There is an eigenmode-expansion monitor placed at the midpoint of the first waveguide. This is a line monitor which extends beyond the waveguide in order to span the entire mode profile including its evanescent tails. The Fourier-transformed fields along this line monitor are used to compute the basis coefficients of the harmonic modes. These are computed separately via the eigenmode solver [MPB](https://mpb.readthedocs.io/en/latest/). This is described in [Mode Decomposition](../Mode_Decomposition.md) where it is also shown that the squared magnitude of the mode coefficient is equivalent to the power (Poynting flux) in the given eigenmode. The ratio of the complex mode coefficients can be used to compute the [S parameters](https://en.wikipedia.org/wiki/Scattering_parameters). In this example, we are computing $|S_{11}|^2$ which is the reflectance (shown in the line prefixed by "refl:,"). Another line monitor could have been placed in the second waveguide to compute the transmittance or $|S_{21}|^2$ into the various guided modes (since the second waveguide is multi mode). The scattered power into the radiative modes can then be computed as $1-|S_{11}|^2-|S_{21}|^2$. As usual, a normalization run is required involving a straight waveguide to compute the power in the source.
+The structure, which can be viewed as a [two-port network](https://en.wikipedia.org/wiki/Two-port_network), consists of a single-mode waveguide of width 1 μm (`w1`) at a wavelength of 6.67 μm and coupled to a second waveguide of width 2 μm (`w2`) via a linearly-sloped taper of variable length `Lt`. The material is silicon with $\varepsilon=12$. The taper geometry is defined using a single [`Prism`](../Python_User_Interface.md#prism) object with eight vertices. PML absorbing boundaries surround the entire cell. An eigenmode current source with $E_z$ (or $\mathcal{S}$) polarization is used to launch the fundamental mode. The dispersion relation (or "band diagram") of the single-mode waveguide is shown in [Tutorial/Eigenmode Source](Eigenmode_Source.md). There is an eigenmode-expansion monitor placed at the midpoint of the first waveguide. This is a line monitor which extends beyond the waveguide in order to span the entire mode profile including its evanescent tails. The Fourier-transformed fields along this line monitor are used to compute the basis coefficients of the harmonic modes. These are computed separately via the eigenmode solver [MPB](https://mpb.readthedocs.io/en/latest/). This is described in [Mode Decomposition](../Mode_Decomposition.md) where it is also shown that the squared magnitude of the mode coefficient is equivalent to the power (Poynting flux) in the given eigenmode. The ratio of the complex mode coefficients can be used to compute the [S parameters](https://en.wikipedia.org/wiki/Scattering_parameters). In this example, we are computing $|S_{11}|^2$ which is the reflectance (shown in the line prefixed by "refl:,"). Another line monitor could have been placed in the second waveguide to compute the transmittance or $|S_{21}|^2$ into the various guided modes (since the second waveguide is multi mode). The scattered power into the radiative modes can then be computed as $1-|S_{11}|^2-|S_{21}|^2$. As usual, a normalization run is required involving a straight waveguide to compute the power in the source.
 
 The structure has mirror symmetry in the $y$ direction which can be exploited to reduce the computation size by a factor of two. This requires using `add_flux` rather than `add_mode_monitor` (which is not optimized for symmetry) and specifying the keyword argument `eig_parity=mp.ODD_Z+mp.EVEN_Y` in the call to `get_eigenmode_coefficients`. Alternatively, the waveguide could have been oriented along an arbitrary oblique direction which would require specifying `direction=mp.NO_DIRECTION` and `kpoint_func` as the waveguide axis. For an example, see [Tutorials/Eigenmode Source/Index-Guided Modes in a Ridge Waveguide](Eigenmode_Source.md#index-guided-modes-in-a-ridge-waveguide).
 
@@ -147,9 +147,9 @@ Diffraction Spectrum 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](Near_to_Far_Field_Spectra.md#diffraction-spectrum-of-a-finite-binary-grating). The unit cell geometry of the grating is shown in the schematic below. The grating 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.
 
-### Transmittance Spectra for Planewave at Normal Incidence
+### Transmissive Diffraction Spectrum for Planewave at Normal Incidence
 
-A pulsed planewave with $E_z$ 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 &mdash; the ratio of the power in the $+x$ direction of the diffracted mode relative to that of the incident planewave &mdash; 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.
+A pulsed planewave with $E_z$ (or $\mathcal{S}$) 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 &mdash; the ratio of the power in the $+x$ direction of the diffracted mode relative to that of the incident planewave &mdash; 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.
 
 The simulation script is in [examples/binary_grating.py](https://github.com/NanoComp/meep/blob/master/python/examples/binary_grating.py). The notebook is [examples/binary_grating.ipynb](https://nbviewer.jupyter.org/github/NanoComp/meep/blob/master/python/examples/binary_grating.ipynb).
 
@@ -247,7 +247,7 @@ for nm in range(nmode):
     print("grating{}:, {:.5f}, {:.2f}, {:.8f}".format(nm,mode_wvl[-1],mode_angle[-1],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 $E_z$ polarization. `EVEN_Y` is necessary since each diffraction order which is based on a given $k_x$ 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 $+k_y + -k_y$ (cosine) mode. As a result, the total transmittance must be halved in this case to obtain the transmittance for the individual $+k_y$ or $-k_y$ mode. For `ODD_Y`, MPB will compute the $+k_y - -k_y$ (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 $H_z$ instead of $E_z$ 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.
+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 $E_z$ (or $\mathcal{S}$) polarization. `EVEN_Y` is necessary since each diffraction order which is based on a given $k_x$ 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 $+k_y + -k_y$ (cosine) mode. As a result, the total transmittance must be halved in this case to obtain the transmittance for the individual $+k_y$ or $-k_y$ mode. For `ODD_Y`, MPB will compute the $+k_y - -k_y$ (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 $H_z$ instead of $E_z$ 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.
 

doc/docs/Python_Tutorials/Near_to_Far_Field_Spectra.md

@@ -135,7 +135,7 @@ plt.show()
 Focusing Properties of a Metasurface Lens
 -----------------------------------------
 
-This example demonstrates how to compute the far-field profile at the focal length of a metasurface lens. The lens design, which is also part of the tutorial, is based on a supercell of binary-grating unit cells. For a review of the binary-grating geometry as well as a demonstration of computing its phasemap, see [Tutorial/Mode Decomposition/Phase Map of a Subwavelength Binary Grating](Mode_Decomposition.md#phase-map-of-a-subwavelength-binary-grating). The far-field calculation of the lens contains two separate components: (1) compute the phasemap of the unit cell as a function of a single geometric parameter, the duty cycle, while keeping its height and periodicity fixed (1.8 and 0.3 μm), and (2) form the supercell lens by tuning the local phase of each of a variable number of unit cells according to the quadratic formula for planar wavefront focusing. The design wavelength is 0.5 μm and the focal length is 0.2 mm. The input source is an $E_z$-polarized planewave at normal incidence.
+This example demonstrates how to compute the far-field profile at the focal length of a metasurface lens. The lens design, which is also part of the tutorial, is based on a supercell of binary-grating unit cells. For a review of the binary-grating geometry as well as a demonstration of computing its phasemap, see [Tutorial/Mode Decomposition/Phase Map of a Subwavelength Binary Grating](Mode_Decomposition.md#phase-map-of-a-subwavelength-binary-grating). The far-field calculation of the lens contains two separate components: (1) compute the phasemap of the unit cell as a function of a single geometric parameter, the duty cycle (also referred to as the filling fraction), while keeping its height and periodicity fixed (1.8 μm and 0.3 μm, respectively), and (2) form the supercell lens by tuning the local phase of each of a variable number of unit cells according to the quadratic formula for planar wavefront focusing. The operating wavelength is 0.5 μm and the focal length is 0.2 mm. The input source is an $E_z$-polarized planewave at normal incidence.
 
 The simulation script is in [examples/metasurface_lens.py](https://github.com/NanoComp/meep/blob/master/python/examples/metasurface_lens.py). The notebook is [examples/metasurface_lens.ipynb](https://nbviewer.jupyter.org/github/NanoComp/meep/blob/master/python/examples/metasurface_lens.ipynb).
 
@@ -146,6 +146,7 @@ import meep as mp
 import numpy as np
 import matplotlib.pyplot as plt
 
+
 resolution = 50         # pixels/μm
 
 dpml = 1.0              # PML thickness
@@ -791,24 +792,24 @@ ff = np.array(ff)
 if mp.am_master():
     plt.figure()
     plt.subplot(1,3,1)
-    plt.plot(np.real(Hz_mon),'bo-',label='DFT')
-    plt.plot(np.real(ff),'ro-',label='N2F')
+    plt.plot(x,np.real(Hz_mon),'bo-',label='DFT')
+    plt.plot(x,np.real(ff),'ro-',label='N2F')
     plt.legend()
-    plt.xlabel('array index')
+    plt.xlabel('$x$ (μm)')
     plt.ylabel('real(Hz)')
 
     plt.subplot(1,3,2)
-    plt.plot(np.imag(Hz_mon),'bo-',label='DFT')
-    plt.plot(np.imag(ff),'ro-',label='N2F')
+    plt.plot(x,np.imag(Hz_mon),'bo-',label='DFT')
+    plt.plot(x,np.imag(ff),'ro-',label='N2F')
     plt.legend()
-    plt.xlabel('array index')
+    plt.xlabel('$x$ (μm)')
     plt.ylabel('imag(Hz)')
 
     plt.subplot(1,3,3)
-    plt.plot(np.abs(Hz_mon),'bo-',label='DFT')
-    plt.plot(np.abs(ff),'ro-',label='N2F')
+    plt.plot(x,np.abs(Hz_mon),'bo-',label='DFT')
+    plt.plot(x,np.abs(ff),'ro-',label='N2F')
     plt.legend()
-    plt.xlabel('array index')
+    plt.xlabel('$x$ (μm)')
     plt.ylabel('|Hz|')
 
     plt.suptitle(f'comparison of near2far and actual DFT fields\n dpad={dpad}, d1={d1}, d2={d2}')
@@ -830,3 +831,7 @@ A closed near-field surface will still disagree with a brute-force far-field cal
 In this example, to demonstrate agreement between the far fields and DFT fields, there are two requirements: (1) the cell size in the $x$ direction via `dpad` needs to be sufficiently large in order to minimize the impact of the spurious radiation from the edge of the near-field surface and (2) the far-field region needs to be sufficiently close to the near-field surface to minimize discrepancies caused by numerical dispersion.
 
 ![center|Comparison of the far fields from the near-to-far field transformation and the DFT fields at the same location for a holey-waveguide cavity.](../images/farfields_vs_DFTfields_holeycavity.png)
+
+When these two conditions are not met as in the example below involving a small `dpad` and large `d2`, the error from the finite truncation and numerical dispersion can be large and therefore result in a significant mismatch between the far fields computed using the near-to-far field transformation versus the actual DFT fields at the same location.
+
+![center|Comparison of the far fields from the near-to-far field transformation dominated by errors and the DFT fields at the same location for a holey-waveguide cavity.](../images/farfields_vs_DFTfields_holeycavity_mismatch.png)