various minor fixes and improvements to tutorials and docstrings

doc/docs/Materials.md

@@ -65,7 +65,7 @@ For a wavelength-dependent, purely-real permittivity (i.e., with no loss) which
 
 $$\varepsilon(\lambda) = 1 + \sum_n \frac{B_n \lambda^2}{\lambda^2 - C_n}$$
 
-where $\lambda$ is the vacuum wavelength, each term containing two coefficients ($B_n$ and $C_n$) can be directly transferred to a Lorentzian polarization field using a simple substitution of variables: $\omega_n=1/\sqrt{C_n}$, $\gamma_n=0$, and $\sigma_n=B_n$. Several examples of importing Sellmeier coefficients from published fitting data including [germanium](https://github.com/NanoComp/meep/blob/master/python/materials.py#L884-L901) (Ge) and [gallium nitride](https://github.com/NanoComp/meep/blob/master/python/materials.py#L1162-L1188) (GaN) are provided in the [Materials Library](#materials-library).
+where $\lambda$ is the vacuum wavelength, each term containing two coefficients ($B_n$ and $C_n$) can be directly transferred to a Lorentzian polarization field using a simple substitution of variables: $\omega_n=1/\sqrt{C_n}$, $\gamma_n=0$, and $\sigma_n=B_n$. Several examples of importing Sellmeier coefficients from published fitting data including [germanium](https://github.com/NanoComp/meep/blob/master/python/materials.py#L1046-L1065) (Ge) and [gallium nitride](https://github.com/NanoComp/meep/blob/master/python/materials.py#L1460-L1500) (GaN) are provided in the [Materials Library](#materials-library).
 
 Numerical Stability
 -------------------

doc/docs/Python_Tutorials/Custom_Source.md

@@ -384,7 +384,7 @@ The previous examples involved emission in the *normal* direction. To investigat
 
 Finally, for cases which involve computing the emission from a collection of $N$ random dipoles in a *single* direction, one can exploit [reciprocity](https://en.wikipedia.org/wiki/Reciprocity_(electromagnetism)) to obtain the same result using a single simulation rather than Monte Carlo sampling of a large number of runs simply by swapping the sources and monitors and running the simulation backwards. This approach is analogous to computing [thermal radiation](https://en.wikipedia.org/wiki/Thermal_radiation) from a blackbody via [Kirchhoff's law for thermal radiation](https://en.wikipedia.org/wiki/Kirchhoff%27s_law_of_thermal_radiation) which states that emissivity and absorptivity are equivalent under thermal equilibrium. However, the approach demonstrated in this tutorial, which was originally developed for computing the extraction efficiency of LEDs in [Optics Express, Vol. 18, pp. 24522-35 (2010)](https://opg.optica.org/oe/fulltext.cfm?uri=oe-18-24-24522), does not require lossy materials.
 
-The reciprocal calculation involves two parts: (1) for the simulation (with the setup shown in the figure below), send an input planewave in the opposite direction of the output mode of the forward calculation and monitor the DFT fields at the same location as the sources in the forward calculation, and (2) in post processing, evaluate an inner product of the DFT fields using a correlation operator of the random currents. Since the currents are uncorrelated in space and consist only of electric currents (because, in practice, they are generated by quantum wells), the correlation operator is a diagonal matrix. Also, since these are 2d simulations with $E_z$ polarization and involve an isotropic dielectric medium, the inner product of the fields is therefore just a sum of $|E_z|^2$ from the DFT monitors (equation 10 of the reference). In 3d, we would need to perform two separate calculations for the $\mathcal{S}$ and $\mathcal{P}$ polarizations and sum the results. (For more general cases, including multiple output channels, spatially correlated emitters, anisotropy, and other complications, see Section 2.3 of [arXiv:2111.13046](https://arxiv.org/abs/2111.13046).)
+The reciprocal calculation involves two parts: (1) for the backward simulation (with the setup shown in the figure below), send an input planewave in the opposite direction of the output mode of the forward calculation and monitor the DFT fields at the same location as the sources in the forward calculation, and (2) in post processing, evaluate an inner product of the DFT fields using a correlation operator of the random currents. Since the currents are uncorrelated in space and consist only of electric currents (because, in practice, they are generated by quantum wells), the correlation operator is a diagonal matrix. Also, since these are 2d simulations with $E_z$ polarization and involve an isotropic dielectric medium, the inner product of the fields is therefore just a sum of $|E_z|^2$ from the DFT monitors (equation 10 of the reference). In 3d, we would need to perform two separate calculations for the $\mathcal{S}$ and $\mathcal{P}$ polarizations and sum the results. (For more general cases, including multiple output channels, spatially correlated emitters, anisotropy, and other complications, see Section 2.3 of [arXiv:2111.13046](https://arxiv.org/abs/2111.13046).)
 
 ![](../images/LED_layout_reciprocity.png#center)
 

doc/docs/Python_User_Interface.md

@@ -1413,7 +1413,8 @@ def get_eigenmode_coefficients(self,
 
 <div class="method_docstring" markdown="1">
 
-Given a flux object and list of band indices `bands` or `DiffractedPlanewave`, return a `namedtuple` with the
+Given a flux object and either a list of indices `bands` (starts at one)
+or a `DiffractedPlanewave` object, return a `namedtuple` with the
 following fields:
 
 + `alpha`: the complex eigenmode coefficients as a 3d NumPy array of size
@@ -7297,22 +7298,22 @@ def __init__(self,
 
 Construct an `Animate2D` object.
 
-+ **`sim=None`** — Optional Simulation object (this has no effect, and is included for backwards compatibility).
++ **`sim`** — Optional `Simulation` object (this has no effect, and is included for backwards compatibility).
 
-+ **`fields=None`** — Optional Field component to record at each time instant.
++ **`fields`** — Optional field component to record at each time instant.
 
-+ **`f=None`** — Optional `matplotlib` figure object that the routine will update
++ **`f`** — Optional `matplotlib` figure object that the routine will update
   on each call. If not supplied, then a new one will be created upon
   initialization.
 
-+ **`realtime=False`** — Whether or not to update a figure window in realtime as
++ **`realtime`** — Whether or not to update a figure window in realtime as
   the simulation progresses. Disabled by default.
 
-+ **`normalize=False`** — Records fields at each time step in memory in a NumPy
++ **`normalize`** — Records fields at each time step in memory in a NumPy
   array and then normalizes the result by dividing by the maximum field value at a
   single point in the cell over all the time snapshots.
 
-+ **`plot_modifiers=None`** — A list of functions that can modify the figure's
++ **`plot_modifiers`** — A list of functions that can modify the figure's
   `axis` object. Each function modifier accepts a single argument, an `axis`
   object, and must return that same axis object. The following modifier changes
   the `xlabel`:
@@ -7325,9 +7326,9 @@ Construct an `Animate2D` object.
   plot_modifiers = [mod1]
 ```
 
-+ **`update_epsilon=False`** — Redraw epsilon on each call. (Useful for topology optimization)
++ **`update_epsilon`** — Redraw epsilon on each call. (Useful for topology optimization)
 
-+ **`nb=False`** — For the animation work in a Jupyter notebook, set to True and use the cell magic:
++ **`nb`** — For the animation work in a Jupyter notebook, set to True and use the cell magic:
     `%matplotlib ipympl`
 + **`**customization_args`** — Customization keyword arguments passed to
   `plot2D()` (i.e. `labels`, `eps_parameters`, `boundary_parameters`, etc.)

doc/docs/Python_User_Interface.md.in

@@ -176,7 +176,7 @@ Given a structure, Meep can decompose the Fourier-transformed fields into a supe
 
 @@ Simulation.get_eigenmode_coefficients @@
 
-The flux object should be created using `add_mode_monitor`.  (You could also use `add_flux`, but with `add_flux` you need to be more careful about symmetries that bisect the flux plane: the `add_flux` object should only be used with `get_eigenmode_coefficients` for modes of the same symmetry, e.g. constrained via `eig_parity`.  On the other hand, the performance of `add_flux` planes benefits more from symmetry.) `eig_vol` is the volume passed to [MPB](https://mpb.readthedocs.io) for the eigenmode calculation (based on interpolating the discretized materials from the Yee grid); in most cases this will simply be the volume over which the frequency-domain fields are tabulated, which is the default (i.e. `flux.where`). `eig_parity` should be one of [`mp.NO_PARITY` (default), `mp.EVEN_Z`, `mp.ODD_Z`, `mp.EVEN_Y`, `mp.ODD_Y`]. It is the parity (= polarization in 2d) of the mode to calculate, assuming the structure has $z$ and/or $y$ mirror symmetry *in the source region*, just as for `EigenModeSource` above. If the structure has both $y$ and $z$ mirror symmetry, you can combine more than one of these, e.g. `EVEN_Z+ODD_Y`. Default is `NO_PARITY`, in which case MPB computes all of the bands which will still be even or odd if the structure has mirror symmetry, of course. This is especially useful in 2d simulations to restrict yourself to a desired polarization. `eig_resolution` is the spatial resolution to use in MPB for the eigenmode calculations. This defaults to twice the Meep `resolution` in which case the structure is linearly interpolated from the Meep pixels. `eig_tolerance` is the tolerance to use in the MPB eigensolver. MPB terminates when the eigenvalues stop changing to less than this fractional tolerance. Defaults to `1e-12`.  (Note that this is the tolerance for the frequency eigenvalue $\omega$; the tolerance for the mode profile is effectively the square root of this.) For examples, see [Tutorial/Mode Decomposition](Python_Tutorials/Mode_Decomposition.md).
+The flux object should be created using `add_mode_monitor`.  (You could also use `add_flux`, but with `add_flux` you need to be more careful about `symmetries` that bisect the flux plane: the `add_flux` object should only be used with `get_eigenmode_coefficients` for modes of the same symmetry, e.g. constrained via `eig_parity`.  On the other hand, the performance of `add_flux` planes benefits more from symmetry.) `eig_vol` is the volume passed to [MPB](https://mpb.readthedocs.io) for the eigenmode calculation (based on interpolating the discretized materials from the Yee grid); in most cases this will simply be the volume over which the frequency-domain fields are tabulated, which is the default (i.e. `flux.where`). `eig_parity` should be one of [`mp.NO_PARITY` (default), `mp.EVEN_Z`, `mp.ODD_Z`, `mp.EVEN_Y`, `mp.ODD_Y`]. It is the parity (= polarization in 2d) of the mode to calculate, assuming the structure has $z$ and/or $y$ mirror symmetry *in the source region*, just as for `EigenModeSource` above. If the structure has both $y$ and $z$ mirror symmetry, you can combine more than one of these, e.g. `EVEN_Z+ODD_Y`. Default is `NO_PARITY`, in which case MPB computes all of the bands which will still be even or odd if the structure has mirror symmetry, of course. This is especially useful in 2d simulations to restrict yourself to a desired polarization. `eig_resolution` is the spatial resolution to use in MPB for the eigenmode calculations. This defaults to twice the Meep `resolution` in which case the structure is linearly interpolated from the Meep pixels. `eig_tolerance` is the tolerance to use in the MPB eigensolver. MPB terminates when the eigenvalues stop changing to less than this fractional tolerance. Defaults to `1e-12`.  (Note that this is the tolerance for the frequency eigenvalue $\omega$; the tolerance for the mode profile is effectively the square root of this.) For examples, see [Tutorial/Mode Decomposition](Python_Tutorials/Mode_Decomposition.md).
 
 Technically, MPB computes $\omega_n(\mathbf{k})$ and then inverts it with Newton's method to find the wavevector $\mathbf{k}$ normal to `eig_vol` and mode for a given frequency; in rare cases (primarily waveguides with *nonmonotonic* dispersion relations, which doesn't usually happen in simple dielectric waveguides), MPB may need you to supply an initial "guess" for $\mathbf{k}$ in order for this Newton iteration to converge.  You can supply this initial guess with `kpoint_func`, which is a function `kpoint_func(f, n)` that supplies a rough initial guess for the $\mathbf{k}$ of band number $n$ at frequency $f=\omega/2\pi$. (By default, the $\mathbf{k}$ components in the plane of the `eig_vol` region are zero.  However, if this region spans the *entire* cell in some directions, and the cell has Bloch-periodic boundary conditions via the `k_point` parameter, then the mode's $\mathbf{k}$ components in those directions will match `k_point` so that the mode satisfies the Meep boundary conditions, regardless of `kpoint_func`.) If `direction` is set to `mp.NO_DIRECTION`, then `kpoint_func` is not only the initial guess and the search direction of the $\mathbf{k}$ vectors, but is also taken to be the direction of the waveguide, allowing you to [detect modes in oblique waveguides](Python_Tutorials/Eigenmode_Source.md#oblique-waveguides) (not perpendicular to the flux plane).
 

python/examples/eps_fit_lorentzian.py

@@ -110,6 +110,9 @@ def lorentzfit(
     idx_end = end_idx[0][-1] + 1
 
     # The fitting function is ε(f) where f is the frequency, rather than ε(λ).
+    # Note: an equally spaced grid of wavelengths results in the larger
+    #       wavelengths having a finer frequency grid than smaller ones.
+    #       This feature may impact the accuracy of the fit.
     freqs = 1000 / wl  # units of 1/μm
     freqs_reduced = freqs[idx_start:idx_end]
     wl_reduced = wl[idx_start:idx_end]

python/simulation.py

@@ -4142,8 +4142,9 @@ def get_eigenmode_coefficients(
         direction: int = mp.AUTOMATIC,
     ) -> NamedTuple:
         """
-        Given a flux object and list of band indices (integers) `bands` or a `DiffractedPlanewave` object,
-        return a `namedtuple` with the following fields:
+        Given a flux object and either a list of indices `bands` (starts at one)
+        or a `DiffractedPlanewave` object, return a `namedtuple` with the
+        following fields:
 
         + `alpha`: the complex eigenmode coefficients as a 3d NumPy array of size
           (`len(bands)`, `flux.nfreqs`, `2`). The last/third dimension refers to modes

python/visualization.py

@@ -1407,22 +1407,22 @@ def __init__(
         """
         Construct an `Animate2D` object.
 
-        + **`sim=None`** — Optional Simulation object (this has no effect, and is included for backwards compatibility).
+        + **`sim`** — Optional `Simulation` object (this has no effect, and is included for backwards compatibility).
 
-        + **`fields=None`** — Optional Field component to record at each time instant.
+        + **`fields`** — Optional field component to record at each time instant.
 
-        + **`f=None`** — Optional `matplotlib` figure object that the routine will update
+        + **`f`** — Optional `matplotlib` figure object that the routine will update
           on each call. If not supplied, then a new one will be created upon
           initialization.
 
-        + **`realtime=False`** — Whether or not to update a figure window in realtime as
+        + **`realtime`** — Whether or not to update a figure window in realtime as
           the simulation progresses. Disabled by default.
 
-        + **`normalize=False`** — Records fields at each time step in memory in a NumPy
+        + **`normalize`** — Records fields at each time step in memory in a NumPy
           array and then normalizes the result by dividing by the maximum field value at a
           single point in the cell over all the time snapshots.
 
-        + **`plot_modifiers=None`** — A list of functions that can modify the figure's
+        + **`plot_modifiers`** — A list of functions that can modify the figure's
           `axis` object. Each function modifier accepts a single argument, an `axis`
           object, and must return that same axis object. The following modifier changes
           the `xlabel`:
@@ -1435,9 +1435,9 @@ def mod1(ax):
           plot_modifiers = [mod1]
         ```
 
-        + **`update_epsilon=False`** — Redraw epsilon on each call. (Useful for topology optimization)
+        + **`update_epsilon`** — Redraw epsilon on each call. (Useful for topology optimization)
 
-        + **`nb=False`** — For the animation work in a Jupyter notebook, set to True and use the cell magic:
+        + **`nb`** — For the animation work in a Jupyter notebook, set to True and use the cell magic:
             `%matplotlib ipympl`
         + **`**customization_args`** — Customization keyword arguments passed to
           `plot2D()` (i.e. `labels`, `eps_parameters`, `boundary_parameters`, etc.)