Implementation of gyrotropic susceptibilities (#863)

* Implement gyrotropic susceptibility class.

* Add Python and Scheme support for gyrotropic media.

* Initialize bias vector in python susceptibility struct.

* Remove "bias" from gyrotropic_susceptibility; the information is already in gyro_tensor.

* Add gyrotropic media to docs

* Minor copyedit

* Fix logic in py_susceptibility_to_susceptibility

* In add_susceptibilities, always pass a 3-vector as gyrotropic bias

* First try at gyrotropic media tutorial

* Fix errors in gyrotropy formulas in notes

* Re-implemention of gyrotropy using LLG equation

* Tweak handling of sigma tensor in gyrotropic case

* Update Python doc.

* Drop 2pi factor from gyrotropic sigma

* Doc updates

* Fix last change to update_P

* Fix printf typo

* Remove spurious 2pi factor in alpha (which is not a rate)

* Minor code tweak

* Use a central-difference scheme for the LLG dynamics, which seems slighly more stable...

* Try implementing the full nonlinear LLG equation

* Add implicit static polarization to gyrotropy implementation

* Put static P back in in subtract_P

* Add gyrotropy example

* Fix; use LOOP_OVER_VOL instead of LOOP_OVER_VOL_OWNED to ensure updating of off-diagonal components

* Clamp the magnitude of the LLG polarization vector.

* Revert inadvertent unrelated change to meep.i

* Minor code cleanup

* Flag "needs_W_notowned" for gyrotropic media

* Update gyrotropic P components explicitly; don't use LOOP_OVER_VOL_OWNED

* Enable needs_P on all components for gyrotropic media

* Fix gyrotropy scheme to track 9 polarization components per unit cell.

* Revert unrelated last change to meep.i

* Avoiding need for allocation of P_tmp in gyrotropy_data.

* Implement num_cinternal_notowned stuff for gyrotropic media

* Update documentation for gyrotropic media, and relax some minor restrictions.

* Add virtual keywords to gyrotropic_susceptibility methods

* Merge latest changes from master

* Remove gyrotropic-dispersion.py (incomplete attempt)

* Complete merge

* Update Materials.md to discuss both Lorentzian and LLG gyrotropic models

* Introduce a new gyrotropy_model enum type, to allow for the LLG model.

* More plumbing to provide support for Landau-Lifshitz-Gilbert type gyrotropy model

* Fix typo in susceptibility update equation

* Fix typos in Faraday rotation formula in docs

* Merge from master

* Reimplement linearized-LLG updating equations

* Fix typo

* Minor code clarification

* Fix Faraday rotation example

* For Landau-Lifshitz-Gilbert model, ignore the magnitude of the bias vector.

* Fix minor hiccup in docs.

* Support dumping and undumping of gyrotropic susceptibilities

* Doc updates and minor tweaks accompanying last merge

* Translate Faraday rotation tutorial from Python to Scheme

* Fix typo in last change

* Fix missing 2pi factor in gyrotropic LLG susceptibility's sigma parameter

* Minor fixes for gyrotropy documentation

* Add Faraday rotation unit test

* Use absolute tolerance (in degrees) for Faraday rotation unit test

* Add faraday rotation test to python/Makefile.am

doc/docs/Materials.md

@@ -167,6 +167,51 @@ $$ \mathbf{E} \; (\textrm{SALT}) = \frac{2 |\theta|}{\hbar \sqrt{\gamma_\perp \g
 
 For a two level gain medium, $\gamma_\parallel = \gamma_{12} + \gamma_{21}$. For more details on applying SALT to atomic media with an arbitrary number of levels, see [Optics Express, Vol. 23, pp. 6455-77, 2015](https://www.osapublishing.org/oe/abstract.cfm?uri=oe-23-5-6455).
 
+Gyrotropic Media
+----------------
+
+(**Experimental feature**) Meep supports gyrotropic media, which break optical reciprocity and give rise to magneto-optical phenomena such as the [Faraday effect](https://en.wikipedia.org/wiki/Faraday_effect). Such materials are used in devices like [Faraday rotators](https://en.wikipedia.org/wiki/Faraday_rotator).
+
+In a gyrotropic medium, the polarization vector undergoes precession around a preferred direction. In the frequency domain, this corresponds to the presence of skew-symmetric off-diagonal components in the ε tensor (for a gyroelectric medium) or the μ tensor (for a gyromagnetic medium). Two different gyrotropy models are supported:
+
+### Gyrotropic Drude-Lorentz Model
+
+The first gyrotropy model is a [Drude-Lorentz](Materials.md#material-dispersion) model with an additional precession, which is intended to describe gyroelectric materials. The polarization equation is
+
+$$\frac{d^2\mathbf{P}_n}{dt^2} + \gamma_n \frac{d\mathbf{P}_n}{dt} - \frac{d\mathbf{P}_n}{dt} \times \mathbf{b}_n + \omega_n^2 \mathbf{P}_n = \sigma_n(\mathbf{x}) \omega_n^2 \mathbf{E}$$
+
+(Optionally, the polarization may be of Drude form, in which case the $\omega_n^2 \mathbf{P}_n$ term on the left is omitted.) The third term on the left side, which breaks time-reversal symmetry, is responsible for the gyrotropy; it typically describes the deflection of electrons flowing within the material by a static external magnetic field. In the $\gamma_n = \omega_n = 0$ limit, the equation of motion reduces to a precession around the "bias vector" $\mathbf{b}_n$:
+
+$$\frac{d\mathbf{P}_n}{dt} = \mathbf{P}_n \times \mathbf{b}_n$$
+
+Hence, the magnitude of the bias vector is the angular frequency of the gyrotropic precession induced by the external field.
+
+### Gyrotropic Saturated Dipole (Linearized Landau-Lifshitz-Gilbert) Model
+
+The second gyrotropy model is a linearized [Landau-Lifshitz-Gilbert equation](https://en.wikipedia.org/wiki/Landau%E2%80%93Lifshitz%E2%80%93Gilbert_equation), suitable for modeling gyromagnetic materials such as ferrites. Its polarization equation of motion is
+
+$$\frac{d\mathbf{P}_n}{dt} = \mathbf{b}_n \times \left( - \sigma_n \mathbf{E} + \omega_n \mathbf{P}_n + \alpha_n \frac{d\mathbf{P}_n}{dt} \right) - \gamma_n \mathbf{P}_n$$
+
+Note: although the above equation is written in terms of electric susceptibilities, this model is typically used for magnetic susceptibilities. Meep places no restriction on the field type that either gyrotropy model can be applied to. As usual, electric and magnetic susceptibilities can be swapped by substituting ε with μ, **E** with **H**, etc.
+
+The Landau-Lifshitz-Gilbert equation describes the precessional motion of a saturated point magnetic dipole in a magnetic field. In the above equation, the variable $\mathbf{P}_n$ represents the linearized deviation of the polarization from its static equilibrium value (assumed to be much larger and aligned parallel to $\mathbf{b}_n$). Note that this equation of motion is completely different from the [Drude-Lorentz equation](Materials.md#material-dispersion), though the constants σ$_n$, ω$_n$, and γ$_n$ play analogous roles (σ$_n$ couples the polarization to the driving field, ω$_n$ is the angular frequency of precession, and γ$_n$ is a damping factor).
+
+In this model, $\mathbf{b}_n$ is taken to be a unit vector (i.e., its magnitude is ignored).
+
+### Frequency Domain Susceptibility Tensors
+
+Suppose $\mathbf{b} = b \hat{z}$, and let all fields have harmonic time-dependence $\exp(-i\omega t)$. Then $\mathbf{P}_n$ is related to the applied field $\mathbf{E}$ by
+
+$$\mathbf{P}_n =  \begin{bmatrix}\chi_\perp & -i\eta & 0 \\ i\eta & \chi_\perp & 0 \\ 0 & 0 & \chi_\parallel \end{bmatrix} \mathbf{E}$$
+
+For the [gyrotropic Lorentzian model](Materials.md#gyrotropic-drude-lorentz-model), the components of the susceptibility tensor are
+
+$$\chi_\perp = \frac{\omega_n^2 \Delta_n \sigma_n}{\Delta_n^2 - \omega^2 b^2},\;\;\; \chi_\parallel = \frac{\omega_n^2 \sigma_n}{\Delta_n}, \;\;\; \eta = \frac{\omega_n^2 \omega b \sigma_n}{\Delta_n^2 - \omega^2 b^2}, \;\;\;\Delta_n \equiv \omega_n^2 - \omega^2 - i\omega\gamma_n$$
+
+And for the [gyrotropic saturated dipole (linearized Landau-Lifshitz-Gilbert) model](Materials.md#gyrotropic-saturated-dipole-linearized-landau-lifshitz-gilbert-model),
+
+$$\chi_\perp = \frac{\sigma_n (\omega_n - i \omega \alpha_n)}{(\omega_n - i \omega \alpha_n)^2 - (\omega + i \gamma_n)^2}, \;\;\; \chi_\parallel = 0, \;\;\;  \eta = \frac{\sigma_n (\omega + i \gamma)}{(\omega_n - i \omega \alpha_n)^2 - (\omega + i \gamma_n)^2}$$
+
 Materials Library
 -----------------
 

doc/docs/Python_Tutorials/Gyrotropic_Media.md

@@ -0,0 +1,124 @@
+---
+# Gyrotropic Media
+---
+
+In this example, we will perform simulations with gyrotropic media. See [Materials](../Materials.md#gyrotropic-media) for more information on how gyrotropy is supported.
+
+[TOC]
+
+### Faraday Rotation
+
+Consider a uniform gyroelectric medium with bias vector $\mathbf{b} = b \hat{z}$. In the frequency domain, the *x* and *y* components of the dielectric tensor have the form
+
+$$\epsilon = \begin{bmatrix}\epsilon_\perp & -i\eta \\ i\eta & \epsilon_\perp \end{bmatrix}$$
+
+The skew-symmetric off-diagonal components give rise to [Faraday rotation](https://en.wikipedia.org/wiki/Faraday_effect): when a plane wave linearly polarized along *x* is launched along the gyrotropy axis *z*, the polarization vector will precess around the gyrotropy axis as the wave propagates. This is the principle behind [Faraday rotators](https://en.wikipedia.org/wiki/Faraday_rotator), devices that act as one-way valves for light.
+
+A plane wave undergoing Faraday rotation can be described by the complex ansatz
+
+$$\begin{bmatrix}E_x \\ E_y\end{bmatrix} = E_0 \begin{bmatrix}\cos(\kappa_c z) \\ \sin(\kappa_c z)\end{bmatrix} e^{i(kz-\omega t)}$$
+
+where $\kappa_c$ is the Faraday rotation (in radians) per unit of propagation distance. Substituting this into the frequency domain Maxwell's equations, with the above dielectric tensor, yields
+
+$$|\kappa_c| = \omega \sqrt{\frac{\mu}{2} \, \left(\epsilon_\perp - \sqrt{\epsilon_\perp^2 - \eta^2}\right)}$$
+
+We model this phenomenon in the simulation script [faraday-rotation.py](https://github.com/NanoComp/meep/blob/master/python/examples/faraday-rotation.py). First, we define a gyroelectric material:
+
+```python
+import meep as mp
+
+## Parameters for a gyrotropic Lorentzian medium
+epsn = 1.5    # background permittivity
+f0   = 1.0    # natural frequency
+gamm = 1e-6   # damping rate
+sn   = 0.1    # sigma parameter
+b0   = 0.15   # magnitude of bias vector
+
+susc = [mp.GyrotropicLorentzianSusceptibility(frequency=f0, gamma=gamma, sigma=sigma,
+                                              bias=mp.Vector3(0, 0, b0))]
+mat = mp.Medium(epsilon=epsn, mu=1, E_susceptibilities=susc)
+```
+
+The `GyrotropicLorentzianSusceptibility` object has a `bias` argument that takes a `Vector3` specifying the gyrotropy vector. In this case, the vector points along *z*, and its magnitude (which specifies the precession frequency) is determined by the variable `b0`. The other arguments play the same role as in an ordinary (non-gyrotropic) [Lorentzian susceptibility](Material_Dispersion.md).
+
+Next, we set up and run the Meep simulation.
+
+```python
+tmax = 100
+L = 20.0
+cell = mp.Vector3(0, 0, L)
+fsrc, src_z = 0.8, -8.5
+pml_layers = [mp.PML(thickness=1.0, direction=mp.Z)]
+
+sources = [mp.Source(mp.ContinuousSource(frequency=fsrc),
+                     component=mp.Ex, center=mp.Vector3(0, 0, src_z))]
+
+sim = mp.Simulation(cell_size=cell, geometry=[], sources=sources,
+                    boundary_layers=pml_layers,
+                    default_material=mat, resolution=50)
+sim.run(until=tmax)
+```
+
+The simulation cell is one pixel wide in the *x* and *y* directions, with periodic boundary conditions. [PMLs](../Perfectly_Matched_Layer.md) are placed in the *z* direction. A `ContinuousSource` emits a wave whose electric field is initially polarized along *x*. We then plot the *x* and *y* components of the electric field versus *z*:
+
+```python
+import numpy as np
+import matplotlib.pyplot as plt
+
+ex_data = sim.get_efield_x().real
+ey_data = sim.get_efield_y().real
+
+z = np.linspace(-L/2, L/2, len(ex_data))
+plt.figure(1)
+plt.plot(z, ex_data, label='Ex')
+plt.plot(z, ey_data, label='Ey')
+plt.xlim(-L/2, L/2)
+plt.xlabel('z')
+plt.legend()
+plt.show()
+```
+
+<center>
+![](../images/Faraday-rotation.png)
+</center>
+
+We see that the wave indeed rotates in the *x*-*y* plane as it travels.
+
+Moreover, we can compare the Faraday rotation rate in these simulation results to theoretical predictions. In the [gyrotropic Lorentzian model](../Materials.md#gyrotropic-media), the ε tensor components are given by
+
+$$\epsilon_\perp = \epsilon_\infty + \frac{\omega_n^2 \Delta_n}{\Delta_n^2 - \omega^2 b^2}\,\sigma_n(\mathbf{x}),\;\;\; \eta = \frac{\omega_n^2 \omega b}{\Delta_n^2 - \omega^2 b^2}\,\sigma_n(\mathbf{x}), \;\;\;\Delta_n \equiv \omega_n^2 - \omega^2 - i\omega\gamma_n$$
+
+From these expressions, we can calculate the rotation rate $\kappa_c$ at the operating frequency, and hence find the $\mathbf{E}_x$ and $\mathbf{E}_y$ field envelopes for the complex ansatz given at the top of this section.
+
+```python
+dfsq = (f0**2 - 1j*fsrc*gamma - fsrc**2)
+eperp = epsn + sn * f0**2 * dfsq / (dfsq**2 - (fsrc*b0)**2)
+eta = sn * f0**2 * fsrc * b0 / (dfsq**2 - (fsrc*b0)**2)
+
+k_gyro = 2*np.pi*fsrc * np.sqrt(0.5*(eperp - np.sqrt(eperp**2 - eta**2)))
+Ex_theory = 0.37 * np.cos(k_gyro * (z - src_z)).real
+Ey_theory = 0.37 * np.sin(k_gyro * (z - src_z)).real
+
+plt.figure(2)
+plt.subplot(2,1,1)
+plt.plot(z, ex_data, label='Ex (MEEP)')
+plt.plot(z, Ex_theory, 'k--')
+plt.plot(z, -Ex_theory, 'k--', label='Ex envelope (theory)')
+plt.xlim(-L/2, L/2); plt.xlabel('z')
+plt.legend(loc='lower right')
+
+plt.subplot(2,1,2)
+plt.plot(z, ey_data, label='Ey (MEEP)')
+plt.plot(z, Ey_theory, 'k--')
+plt.plot(z, -Ey_theory, 'k--', label='Ey envelope (theory)')
+plt.xlim(-L/2, L/2); plt.xlabel('z')
+plt.legend(loc='lower right')
+plt.tight_layout()
+plt.show()
+```
+
+As shown in the figure below, the results are in excellent agreement:
+
+<center>
+![](../images/Faraday-rotation-comparison.png)
+</center>

doc/docs/Python_User_Interface.md

@@ -400,6 +400,38 @@ The noise has root-mean square amplitude σ $\times$ `noise_amp`.
 
 This is a somewhat unusual polarizable medium, a Lorentzian susceptibility with a random noise term added into the damped-oscillator equation at each point. This can be used to directly model thermal radiation in both the [far field](http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.93.213905) and the [near field](http://math.mit.edu/~stevenj/papers/RodriguezIl11.pdf). Note, however that it is more efficient to [compute far-field thermal radiation using Kirchhoff's law](http://www.simpetus.com/projects.html#meep_thermal_radiation) of radiation, which states that emissivity equals absorptivity. Near-field thermal radiation can usually be computed more efficiently using frequency-domain methods, e.g. via [SCUFF-EM](https://github.com/HomerReid/scuff-em), as described e.g. [here](http://doi.org/10.1103/PhysRevB.92.134202) or [here](http://doi.org/10.1103/PhysRevB.88.054305).
 
+### GyrotropicLorentzianSusceptibility or GyrotropicDrudeSusceptibility
+
+(**Experimental feature**) Specifies a single dispersive [gyrotropic susceptibility](Materials.md#gyrotropic-media) of [Lorentzian (damped harmonic oscillator) or Drude form](Materials.md#gyrotropic-drude-lorentz-model). Its parameters are `sigma`, `frequency`, and `gamma`, which have the [usual meanings](#susceptibility), and an additional 3-vector `bias`:
+
+**`bias` [`Vector3`]**
+—
+The gyrotropy vector.  Its direction determines the orientation of the gyrotropic response, and the magnitude is the precession frequency $|\mathbf{b}_n|/2\pi$.
+
+### GyrotropicSaturatedSusceptibility
+
+(**Experimental feature**) Specifies a single dispersive [gyrotropic susceptibility](Materials.md#gyrotropic-media) governed by a [linearized Landau-Lifshitz-Gilbert equation](Materials.md#gyrotropic-saturated-dipole-linearized-landau-lifshitz-gilbert-model). This class takes parameters `sigma`, `frequency`, and `gamma`, whose meanings are different from the Lorentzian and Drude case. It also takes a 3-vector `bias` parameter and an `alpha` parameter:
+
+**`sigma` [`number`]**
+—
+The coupling factor $\sigma_n / 2\pi$ between the polarization and the driving field. In magnetic ferrites, this is the Larmor precession frequency at the saturation field.
+
+**`frequency` [`number`]**
+—
+The Larmor precession frequency, $f_n = \omega_n / 2\pi$.
+
+**`gamma` [`number`]**
+—
+The loss rate $\gamma_n / 2\pi$ in the off-diagonal response.
+
+**`alpha` [`number`]**
+—
+The loss factor $\alpha_n$ in the diagonal response. Note that this parameter is dimensionless and contains no 2π factor.
+
+**`bias` [`Vector3`]**
+—
+Vector specifying the orientation of the gyrotropic response. Unlike the similarly-named `bias` parameter for the [gyrotropic Lorentzian/Drude susceptibilities](#gyrotropiclorentziansusceptibility-or-gyrotropicdrudesusceptibility), the magnitude is ignored; instead, the relevant precession frequencies are determined by the `sigma` and `frequency` parameters.
+
 ### Vector3
 
 Properties:

doc/docs/Scheme_Tutorials/Gyrotropic_Media.md

@@ -0,0 +1,90 @@
+---
+# Gyrotropic Media
+---
+
+In this example, we will perform simulations with gyrotropic media. See [Materials](../Materials.md#gyrotropic-media) for more information on how gyrotropy is supported.
+
+[TOC]
+
+### Faraday Rotation
+
+Consider a uniform gyroelectric medium with bias vector $\mathbf{b} = b \hat{z}$. In the frequency domain, the *x* and *y* components of the dielectric tensor have the form
+
+$$\epsilon = \begin{bmatrix}\epsilon_\perp & -i\eta \\ i\eta & \epsilon_\perp \end{bmatrix}$$
+
+The skew-symmetric off-diagonal components give rise to [Faraday rotation](https://en.wikipedia.org/wiki/Faraday_effect): when a plane wave linearly polarized along *x* is launched along the gyrotropy axis *z*, the polarization vector will precess around the gyrotropy axis as the wave propagates. This is the principle behind [Faraday rotators](https://en.wikipedia.org/wiki/Faraday_rotator), devices that act as one-way valves for light.
+
+A plane wave undergoing Faraday rotation can be described by the complex ansatz
+
+$$\begin{bmatrix}E_x \\ E_y\end{bmatrix} = E_0 \begin{bmatrix}\cos(\kappa_c z) \\ \sin(\kappa_c z)\end{bmatrix} e^{i(kz-\omega t)}$$
+
+where $\kappa_c$ is the Faraday rotation (in radians) per unit of propagation distance. Substituting this into the frequency domain Maxwell's equations, with the above dielectric tensor, yields
+
+$$|\kappa_c| = \omega \sqrt{\frac{\mu}{2} \, \left(\epsilon_\perp - \sqrt{\epsilon_\perp^2 - \eta^2}\right)}$$
+
+We model this phenomenon in the simulation script [faraday-rotation.ctl](https://github.com/NanoComp/meep/blob/master/scheme/examples/faraday-rotation.ctl). First, we define a gyroelectric material:
+
+```scm
+(define-param epsn 1.5)    ; background permittivity
+(define-param f0 1.0)      ; natural frequency
+(define-param g0 1e-6)     ; damping rate
+(define-param sn 0.1)      ; sigma parameter
+(define-param b0 0.15)     ; magnitude of bias vector
+
+(set! default-material
+      (make dielectric
+            (epsilon epsn)
+            (E-susceptibilities
+             (make gyrotropic-lorentzian-susceptibility
+                   (frequency f0)
+                   (sigma sn)
+                   (gamma g0)
+                   (bias (vector3 0 0 b0))))))
+```
+
+The `gyrotropic-lorentzian-susceptibility` object has a `bias` argument that takes a `vector3` specifying the gyrotropy vector. In this case, the vector points along *z*, and its magnitude (which specifies the precession frequency) is determined by the variable `b0`. The other arguments play the same role as in an ordinary (non-gyrotropic) [Lorentzian susceptibility](Material_Dispersion.md).
+
+Next, we set up and run the Meep simulation.
+
+```scm
+(define-param tmax 100)
+(define-param L 20.0)
+(define-param fsrc 0.8)
+(define-param src-z -8.5)
+(set-param! resolution 50)
+
+(set! geometry-lattice (make lattice (size 0 0 L)))
+
+(set! pml-layers (list (make pml (thickness 1.0) (direction Z))))
+
+(set! sources (list
+               (make source
+                     (src (make continuous-src (frequency fsrc)))
+                     (component Ex)
+                     (center (vector3 0 0 src-z)))))
+
+(run-until tmax
+           (to-appended "efields"
+                        (at-end output-efield-x)
+                        (at-end output-efield-y)))
+```
+
+The simulation cell is one pixel wide in the *x* and *y* directions, with periodic boundary conditions. [PMLs](../Perfectly_Matched_Layer.md) are placed in the *z* direction. A `ContinuousSource` emits a wave whose electric field is initially polarized along *x*.
+
+After running the simulation, the `ex` and `ey` datasets in `faraday-rotation-efields.h5` contain the values of $\mathbf{E}_x$ and $\mathbf{E}_y$. These are plotted against *z* in the figure below:
+
+<center>
+![](../images/Faraday-rotation.png)
+</center>
+
+We see that the wave indeed rotates in the *x*-*y* plane as it travels.
+
+Moreover, we can compare the Faraday rotation rate in these simulation results to theoretical predictions. In the [gyrotropic Lorentzian model](../Materials.md#gyrotropic-media), the ε tensor components are given by
+
+$$\epsilon_\perp = \epsilon_\infty + \frac{\omega_n^2 \Delta_n}{\Delta_n^2 - \omega^2 b^2}\,\sigma_n(\mathbf{x}),\;\;\; \eta = \frac{\omega_n^2 \omega b}{\Delta_n^2 - \omega^2 b^2}\,\sigma_n(\mathbf{x}), \;\;\;\Delta_n \equiv \omega_n^2 - \omega^2 - i\omega\gamma_n$$
+
+From these expressions, we can calculate the rotation rate $\kappa_c$ at the operating frequency, and hence find the $\mathbf{E}_x$ and $\mathbf{E}_y$ field envelopes for the complex ansatz given at the top of this section. As shown in the figure below, the results are in excellent agreement:
+
+<center>
+![](../images/Faraday-rotation-comparison.png)
+</center>