diff --git a/doc/docs/Materials.md b/doc/docs/Materials.md
index c71885402..5ea4a16d2 100644
--- a/doc/docs/Materials.md
+++ b/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
-----------------
diff --git a/doc/docs/Python_Tutorials/Gyrotropic_Media.md b/doc/docs/Python_Tutorials/Gyrotropic_Media.md
new file mode 100644
index 000000000..a10250f05
--- /dev/null
+++ b/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()
+```
+
+
+
+
+
+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:
+
+
+
+
diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md
index 88b0050e6..d8a7c86fe 100644
--- a/doc/docs/Python_User_Interface.md
+++ b/doc/docs/Python_User_Interface.md
@@ -392,6 +392,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:
diff --git a/doc/docs/Scheme_Tutorials/Gyrotropic_Media.md b/doc/docs/Scheme_Tutorials/Gyrotropic_Media.md
new file mode 100644
index 000000000..0164be888
--- /dev/null
+++ b/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:
+
+
+
+
+
+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:
+
+
+
+
diff --git a/doc/docs/Scheme_User_Interface.md b/doc/docs/Scheme_User_Interface.md
index c422c2945..871327c48 100644
--- a/doc/docs/Scheme_User_Interface.md
+++ b/doc/docs/Scheme_User_Interface.md
@@ -339,6 +339,38 @@ Specifies a single dispersive susceptibility of Lorentzian (damped harmonic osci
—
The noise has root-mean square amplitude σ $\times$ `noise-amp`.
+### gyrotropic-lorentzian-susceptibility or gyrotropic-drude-susceptibility
+
+(**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$.
+
+### gyrotropic-saturated-susceptibility
+
+(**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.
+
### geometric-object
This class, and its descendants, are used to specify the solid geometric objects that form the dielectric structure being simulated. The base class is:
diff --git a/doc/docs/images/Faraday-rotation-comparison.png b/doc/docs/images/Faraday-rotation-comparison.png
new file mode 100644
index 000000000..9cf49a63f
Binary files /dev/null and b/doc/docs/images/Faraday-rotation-comparison.png differ
diff --git a/doc/docs/images/Faraday-rotation.png b/doc/docs/images/Faraday-rotation.png
new file mode 100644
index 000000000..cfc9b210c
Binary files /dev/null and b/doc/docs/images/Faraday-rotation.png differ
diff --git a/mkdocs.yml b/mkdocs.yml
index b3bbe07dc..fc67d59c7 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -56,6 +56,7 @@ pages:
- 'Tutorial/Near-to-Far Field Spectra': 'Python_Tutorials/Near_to_Far_Field_Spectra.md'
- 'Tutorial/Local Density of States': 'Python_Tutorials/Local_Density_of_States.md'
- 'Tutorial/Optical Forces': 'Python_Tutorials/Optical_Forces.md'
+ - 'Tutorial/Gyrotropic Media': 'Python_Tutorials/Gyrotropic_Media.md'
- 'Tutorial/Multilevel Atomic Susceptibility': 'Python_Tutorials/Multilevel_Atomic_Susceptibility.md'
- 'Tutorial/Frequency-Domain Solver': 'Python_Tutorials/Frequency_Domain_Solver.md'
- 'Tutorial/Eigenmode Source': 'Python_Tutorials/Eigenmode_Source.md'
@@ -72,6 +73,7 @@ pages:
- 'Tutorial/Near-to-Far Field Spectra': 'Scheme_Tutorials/Near_to_Far_Field_Spectra.md'
- 'Tutorial/Local Density of States': 'Scheme_Tutorials/Local_Density_of_States.md'
- 'Tutorial/Optical Forces': 'Scheme_Tutorials/Optical_Forces.md'
+ - 'Tutorial/Gyrotropic Media': 'Scheme_Tutorials/Gyrotropic_Media.md'
- 'Tutorial/Multilevel Atomic Susceptibility': 'Scheme_Tutorials/Multilevel_Atomic_Susceptibility.md'
- 'Tutorial/Frequency-Domain Solver': 'Scheme_Tutorials/Frequency_Domain_Solver.md'
- 'Tutorial/Eigenmode Source': 'Scheme_Tutorials/Eigenmode_Source.md'
diff --git a/python/Makefile.am b/python/Makefile.am
index fd1197209..9e9ad2999 100644
--- a/python/Makefile.am
+++ b/python/Makefile.am
@@ -43,6 +43,7 @@ TESTS = \
$(TEST_DIR)/cyl_ellipsoid.py \
$(TEST_DIR)/dft_energy.py \
$(TEST_DIR)/dft_fields.py \
+ $(TEST_DIR)/faraday_rotation.py \
$(TEST_DIR)/field_functions.py \
$(TEST_DIR)/force.py \
$(TEST_DIR)/fragment_stats.py \
diff --git a/python/examples/faraday-rotation.py b/python/examples/faraday-rotation.py
new file mode 100644
index 000000000..449c18b10
--- /dev/null
+++ b/python/examples/faraday-rotation.py
@@ -0,0 +1,69 @@
+# From the Meep tutorial: plotting Faraday rotation of a linearly polarized plane wave
+
+import meep as mp
+
+## Parameters for a gyrotropic Lorentzian medium
+epsn = 1.5 # background permittivity
+f0 = 1.0 # natural frequency
+gamma = 1e-6 # damping rate
+sn = 0.1 # sigma parameter
+b0 = 0.15 # magnitude of bias vector
+
+susc = [mp.GyrotropicLorentzianSusceptibility(frequency=f0, gamma=gamma, sigma=sn,
+ bias=mp.Vector3(0, 0, b0))]
+mat = mp.Medium(epsilon=epsn, mu=1, E_susceptibilities=susc)
+
+## Set up and run the Meep simulation:
+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)
+
+## Plot results:
+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()
+
+## Comparison with analytic result:
+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()
diff --git a/python/geom.py b/python/geom.py
index a3dc12715..e3dc6ea61 100644
--- a/python/geom.py
+++ b/python/geom.py
@@ -280,6 +280,26 @@ def __init__(self, noise_amp=0.0, **kwargs):
super(NoisyDrudeSusceptibility, self).__init__(**kwargs)
self.noise_amp = noise_amp
+class GyrotropicLorentzianSusceptibility(LorentzianSusceptibility):
+
+ def __init__(self, bias=Vector3(), **kwargs):
+ super(GyrotropicLorentzianSusceptibility, self).__init__(**kwargs)
+ self.bias = bias
+
+class GyrotropicDrudeSusceptibility(DrudeSusceptibility):
+
+ def __init__(self, bias=Vector3(), **kwargs):
+ super(GyrotropicDrudeSusceptibility, self).__init__(**kwargs)
+ self.bias = bias
+
+class GyrotropicSaturatedSusceptibility(Susceptibility):
+
+ def __init__(self, bias=Vector3(), frequency=0.0, gamma=0.0, alpha=0.0, **kwargs):
+ super(GyrotropicSaturatedSusceptibility, self).__init__(**kwargs)
+ self.frequency = frequency
+ self.gamma = gamma
+ self.bias = bias
+ self.alpha = alpha
class MultilevelAtom(Susceptibility):
diff --git a/python/meep.i b/python/meep.i
index 12a020f52..af274db9d 100644
--- a/python/meep.i
+++ b/python/meep.i
@@ -1371,6 +1371,9 @@ PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where
Ellipsoid,
FreqRange,
GeometricObject,
+ GyrotropicDrudeSusceptibility,
+ GyrotropicLorentzianSusceptibility,
+ GyrotropicSaturatedSusceptibility,
Lattice,
LorentzianSusceptibility,
Matrix,
diff --git a/python/tests/faraday_rotation.py b/python/tests/faraday_rotation.py
new file mode 100644
index 000000000..140868b5b
--- /dev/null
+++ b/python/tests/faraday_rotation.py
@@ -0,0 +1,101 @@
+from __future__ import division
+
+import unittest
+import numpy as np
+import meep as mp
+
+## Farady rotation rate for gyrotropic Lorentzian medium
+def kgyro_lorentzian(freq, epsn, f0, gamma, sigma, b0):
+ dfsq = (f0**2 - 1j*freq*gamma - freq**2)
+ eperp = epsn + sigma * f0**2 * dfsq / (dfsq**2 - (freq*b0)**2)
+ eta = sigma * f0**2 * freq * b0 / (dfsq**2 - (freq*b0)**2)
+ return 2*np.pi*freq * np.sqrt(0.5*(eperp - np.sqrt(eperp**2 - eta**2)))
+
+## Farady rotation rate for gyrotropic Drude medium
+def kgyro_drude(freq, epsn, f0, gamma, sigma, b0):
+ dfsq = - 1j*freq*gamma - freq**2
+ eperp = epsn + sigma * f0**2 * dfsq / (dfsq**2 - (freq*b0)**2)
+ eta = sigma * f0**2 * freq * b0 / (dfsq**2 - (freq*b0)**2)
+ return 2*np.pi*freq * np.sqrt(0.5*(eperp - np.sqrt(eperp**2 - eta**2)))
+
+## Farady rotation rate for Landau-Lifshitz-Gilbert medium
+def kgyro_llg(freq, epsn, f0, gamma, sigma, alpha):
+ df1 = f0 - 1j*freq*alpha
+ df2 = freq + 1j*gamma
+ eperp = epsn + sigma * df1/(df1**2 - df2**2)
+ eta = sigma * df2 / (df1**2 - df2**2)
+ return 2*np.pi*freq * np.sqrt(0.5*(eperp - np.sqrt(eperp**2 - eta**2)))
+
+class TestFaradayRotation(unittest.TestCase):
+ ## Simulate a linearly polarized plane wave traveling along the gyrotropy axis.
+ ## Extract Faraday rotation angle by comparing the Ex and Ey amplitudes, and
+ ## compare to a theoretical result corresponding to rotation rate KPRED.
+ ## The default acceptable tolerance TOL is 1.5 degrees.
+ def check_rotation(self, mat, L, fsrc, zsrc, resolution, tmax, zout, kpred, tol=1.5):
+ cell = mp.Vector3(0, 0, L)
+ 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, zsrc))]
+
+ self.sim = mp.Simulation(cell_size=cell, geometry=[], sources=sources,
+ boundary_layers=pml_layers,
+ default_material=mat, resolution=resolution)
+
+ record_vol = mp.Volume(center=mp.Vector3(0, 0, zout))
+ record_Ex, record_Ey, record_t = [], [], []
+
+ def record_ex_ey(sim):
+ record_Ex.append(sim.get_array(vol=record_vol, component=mp.Ex))
+ record_Ey.append(sim.get_array(vol=record_vol, component=mp.Ey))
+ record_t.append(sim.meep_time())
+
+ self.sim.run(mp.after_time(0.5*tmax, mp.at_every(1e-6, record_ex_ey)), until=tmax)
+
+ ex_rel = np.amax(abs(np.fft.fft(record_Ex)))
+ ey_rel = np.amax(abs(np.fft.fft(record_Ey)))
+ result = np.arctan2(ey_rel, ex_rel) * 180/np.pi
+
+ Ex_theory = np.abs(np.cos(kpred * (zout - zsrc)).real)
+ Ey_theory = np.abs(np.sin(kpred * (zout - zsrc)).real)
+ expected = np.arctan2(Ey_theory, Ex_theory) * 180/np.pi
+
+ print("Rotation angle (in degrees): {}, expected {}\n".format(result, expected))
+ np.testing.assert_allclose(expected, result, atol=tol)
+
+ def test_faraday_rotation(self):
+ L, zsrc, zout = 12.0, -4.5, 4.0
+ freq, tmax = 0.8, 100.0
+ resolution = 24
+
+ ## Test gyrotropic Lorentzian medium
+ epsn, f0, gamma, sn, b0 = 1.5, 1.0, 1e-3, 0.1, 0.15
+ susc = [mp.GyrotropicLorentzianSusceptibility(frequency=f0, gamma=gamma, sigma=sn,
+ bias=mp.Vector3(0, 0, b0))]
+ mat = mp.Medium(epsilon=epsn, mu=1, E_susceptibilities=susc)
+ k = kgyro_lorentzian(freq, epsn, f0, gamma, sn, b0)
+ print('=' * 24)
+ print("Testing Faraday rotation for gyrotropic Lorentzian model...")
+ self.check_rotation(mat, L, freq, zsrc, resolution, tmax, zout, k)
+
+ ## Test gyrotropic Drude medium
+ susc = [mp.GyrotropicDrudeSusceptibility(frequency=f0, gamma=gamma, sigma=sn,
+ bias=mp.Vector3(0, 0, b0))]
+ mat = mp.Medium(epsilon=epsn, mu=1, E_susceptibilities=susc)
+ k = kgyro_drude(freq, epsn, f0, gamma, sn, b0)
+ print('=' * 24)
+ print("Testing Faraday rotation for gyrotropic Drude model...")
+ self.check_rotation(mat, L, freq, zsrc, resolution, tmax, zout, k)
+
+ ## Test Landau-Lifshitz-Gilbert medium
+ alpha = 1e-5
+ susc = [mp.GyrotropicSaturatedSusceptibility(frequency=f0, gamma=gamma, sigma=sn,
+ alpha=alpha,
+ bias=mp.Vector3(0, 0, 1.0))]
+ mat = mp.Medium(epsilon=epsn, mu=1, E_susceptibilities=susc)
+ k = kgyro_llg(freq, epsn, f0, gamma, sn, alpha)
+ print('=' * 24)
+ print("Testing Faraday rotation for Landau-Lifshitz-Gilbert model...")
+ self.check_rotation(mat, L, freq, zsrc, resolution, tmax, zout, k)
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/python/typemap_utils.cpp b/python/typemap_utils.cpp
index dc0b2e253..6e3d7029b 100644
--- a/python/typemap_utils.cpp
+++ b/python/typemap_utils.cpp
@@ -329,7 +329,10 @@ static int py_susceptibility_to_susceptibility(PyObject *po, susceptibility_stru
s->frequency = 0;
s->gamma = 0;
+ s->alpha = 0;
s->noise_amp = 0;
+ s->bias.x = s->bias.y = s->bias.z = 0;
+ s->saturated_gyrotropy = false;
s->transitions.resize(0);
s->initial_populations.resize(0);
@@ -345,6 +348,15 @@ static int py_susceptibility_to_susceptibility(PyObject *po, susceptibility_stru
if (!get_attr_dbl(po, &s->noise_amp, "noise_amp")) { return 0; }
}
+ if (PyObject_HasAttrString(po, "bias")) {
+ if (!get_attr_v3(po, &s->bias, "bias")) return 0;
+ }
+
+ if (PyObject_HasAttrString(po, "alpha")) {
+ s->saturated_gyrotropy = true;
+ if (!get_attr_dbl(po, &s->alpha, "alpha")) { return 0; }
+ }
+
if (PyObject_HasAttrString(po, "transitions")) {
// MultilevelAtom
PyObject *py_trans = PyObject_GetAttrString(po, "transitions");
@@ -459,7 +471,28 @@ static PyObject *susceptibility_to_py_obj(susceptibility_struct *s) {
PyObject *res;
PyObject *args = PyTuple_New(0);
- if (s->noise_amp == 0) {
+ if (s->saturated_gyrotropy || s->bias.x || s->bias.y || s->bias.z) {
+ if (s->saturated_gyrotropy) {
+ PyObject *py_gyrotropic_class = PyObject_GetAttrString(geom_mod, "GyrotropicSaturatedSusceptibility");
+ res = PyObject_Call(py_gyrotropic_class, args, NULL);
+ Py_DECREF(py_gyrotropic_class);
+ PyObject *py_alpha = PyFloat_FromDouble(s->alpha);
+ PyObject_SetAttrString(res, "alpha", py_alpha);
+ Py_DECREF(py_alpha);
+ } else if (s->drude) {
+ PyObject *py_gyrotropic_drude_class = PyObject_GetAttrString(geom_mod, "GyrotropicDrudeSusceptibility");
+ res = PyObject_Call(py_gyrotropic_drude_class, args, NULL);
+ Py_DECREF(py_gyrotropic_drude_class);
+ } else {
+ PyObject *py_gyrotropic_lorentz_class =
+ PyObject_GetAttrString(geom_mod, "GyrotropicLorentzianSusceptibility");
+ res = PyObject_Call(py_gyrotropic_lorentz_class, args, NULL);
+ Py_DECREF(py_gyrotropic_lorentz_class);
+ }
+ PyObject *py_bias = vec2py(vector3_to_vec(s->bias));
+ PyObject_SetAttrString(res, "bias", py_bias);
+ Py_DECREF(py_bias);
+ } else if (s->noise_amp == 0) {
if (s->drude) {
PyObject *py_drude_class = PyObject_GetAttrString(geom_mod, "DrudeSusceptibility");
res = PyObject_Call(py_drude_class, args, NULL);
diff --git a/scheme/examples/faraday-rotation.ctl b/scheme/examples/faraday-rotation.ctl
new file mode 100644
index 000000000..bd287eda3
--- /dev/null
+++ b/scheme/examples/faraday-rotation.ctl
@@ -0,0 +1,40 @@
+;; From the Meep tutorial: plotting Faraday rotation of a linearly polarized plane wave
+
+;; Parameters for a gyrotropic Lorentzian medium
+(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))))))
+
+;; Set up and run the Meep simulation:
+(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)))
diff --git a/scheme/meep.scm.in b/scheme/meep.scm.in
index a4ea946bd..7c312c83a 100644
--- a/scheme/meep.scm.in
+++ b/scheme/meep.scm.in
@@ -57,6 +57,16 @@
(define-property noise-amp no-default 'number))
(define-class noisy-drude-susceptibility drude-susceptibility
(define-property noise-amp no-default 'number))
+(define-class gyrotropic-lorentzian-susceptibility lorentzian-susceptibility
+ (define-property bias no-default 'vector3))
+(define-class gyrotropic-drude-susceptibility drude-susceptibility
+ (define-property bias no-default 'vector3))
+
+(define-class gyrotropic-saturated-susceptibility susceptibility
+ (define-property frequency no-default 'number)
+ (define-property gamma no-default 'number)
+ (define-property bias no-default 'vector3)
+ (define-property alpha no-default 'number))
(define polarizability lorentzian-susceptibility) ; backwards compat
(define omega frequency) ; backwards compat
diff --git a/scheme/structure.cpp b/scheme/structure.cpp
index 43c8f49a2..d71626bcf 100644
--- a/scheme/structure.cpp
+++ b/scheme/structure.cpp
@@ -1295,9 +1295,14 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s)
master_printf("noisy lorentzian susceptibility: frequency=%g, gamma=%g, amp = %g\n",
d->frequency, d->gamma, nd->noise_amp);
sus = new meep::noisy_lorentzian_susceptibility(nd->noise_amp, d->frequency, d->gamma);
+ } else if (d->which_subclass == lorentzian_susceptibility::GYROTROPIC_LORENTZIAN_SUSCEPTIBILITY) {
+ gyrotropic_lorentzian_susceptibility *gd = d->subclass.gyrotropic_lorentzian_susceptibility_data;
+ master_printf("gyrotropic lorentzian susceptibility: bias=(%g,%g,%g), frequency=%g, gamma=%g\n",
+ gd->bias.x, gd->bias.y, gd->bias.z, d->frequency, d->gamma);
+ sus = new meep::gyrotropic_susceptibility(vector3_to_vec(gd->bias), d->frequency, d->gamma,
+ 0.0, meep::GYROTROPIC_LORENTZIAN);
} else { // just a Lorentzian
- master_printf("lorentzian susceptibility: frequency=%g, gamma=%g\n", d->frequency,
- d->gamma);
+ master_printf("lorentzian susceptibility: frequency=%g, gamma=%g\n", d->frequency, d->gamma);
sus = new meep::lorentzian_susceptibility(d->frequency, d->gamma);
}
break;
@@ -1308,14 +1313,27 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s)
noisy_drude_susceptibility *nd = d->subclass.noisy_drude_susceptibility_data;
master_printf("noisy drude susceptibility: frequency=%g, gamma=%g, amp = %g\n",
d->frequency, d->gamma, nd->noise_amp);
- sus = new meep::noisy_lorentzian_susceptibility(nd->noise_amp, d->frequency, d->gamma,
- true);
+ sus = new meep::noisy_lorentzian_susceptibility(nd->noise_amp, d->frequency, d->gamma, true);
+ } else if (d->which_subclass == drude_susceptibility::GYROTROPIC_DRUDE_SUSCEPTIBILITY) {
+ gyrotropic_drude_susceptibility *gd = d->subclass.gyrotropic_drude_susceptibility_data;
+ master_printf("gyrotropic drude susceptibility: bias=(%g,%g,%g), frequency=%g, gamma=%g\n",
+ gd->bias.x, gd->bias.y, gd->bias.z, d->frequency, d->gamma);
+ sus = new meep::gyrotropic_susceptibility(vector3_to_vec(gd->bias), d->frequency, d->gamma,
+ 0.0, meep::GYROTROPIC_DRUDE);
} else { // just a Drude
master_printf("drude susceptibility: frequency=%g, gamma=%g\n", d->frequency, d->gamma);
sus = new meep::lorentzian_susceptibility(d->frequency, d->gamma, true);
}
break;
}
+ case susceptibility::GYROTROPIC_SATURATED_SUSCEPTIBILITY: {
+ gyrotropic_saturated_susceptibility *d = p->user_s.subclass.gyrotropic_saturated_susceptibility_data;
+ master_printf("gyrotropic Landau-Lifshitz-Gilbert-type susceptibility: bias=(%g,%g,%g), frequency=%g, gamma=%g, alpha=%g\n",
+ d->bias.x, d->bias.y, d->bias.z, d->frequency, d->gamma, d->alpha);
+ sus = new meep::gyrotropic_susceptibility(vector3_to_vec(d->bias), d->frequency, d->gamma,
+ d->alpha, meep::GYROTROPIC_SATURATED);
+ break;
+ }
case susceptibility::MULTILEVEL_ATOM: {
multilevel_atom *d = p->user_s.subclass.multilevel_atom_data;
sus = make_multilevel_sus(d);
diff --git a/src/material_data.hpp b/src/material_data.hpp
index 0610d7fdb..20c6f9823 100644
--- a/src/material_data.hpp
+++ b/src/material_data.hpp
@@ -59,10 +59,13 @@ struct transition {
typedef struct susceptibility_struct {
vector3 sigma_offdiag;
vector3 sigma_diag;
+ vector3 bias;
double frequency;
double gamma;
+ double alpha;
double noise_amp;
bool drude;
+ bool saturated_gyrotropy;
bool is_file;
std::vector transitions;
std::vector initial_populations;
diff --git a/src/meep.hpp b/src/meep.hpp
index 7689334e2..e387006cb 100644
--- a/src/meep.hpp
+++ b/src/meep.hpp
@@ -272,6 +272,45 @@ class noisy_lorentzian_susceptibility : public lorentzian_susceptibility {
double noise_amp;
};
+typedef enum { GYROTROPIC_LORENTZIAN, GYROTROPIC_DRUDE, GYROTROPIC_SATURATED } gyrotropy_model;
+
+/* gyrotropic susceptibility */
+class gyrotropic_susceptibility : public susceptibility {
+public:
+ gyrotropic_susceptibility(const vec &bias, double omega_0, double gamma, double alpha=0.0,
+ gyrotropy_model model=GYROTROPIC_LORENTZIAN);
+ virtual susceptibility *clone() const { return new gyrotropic_susceptibility(*this); }
+
+ virtual void *new_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], const grid_volume &gv) const;
+ virtual void init_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], double dt,
+ const grid_volume &gv, void *data) const;
+ virtual void *copy_internal_data(void *data) const;
+
+ virtual bool needs_P(component c, int cmp, realnum *W[NUM_FIELD_COMPONENTS][2]) const;
+ virtual void update_P(realnum *W[NUM_FIELD_COMPONENTS][2],
+ realnum *W_prev[NUM_FIELD_COMPONENTS][2], double dt,
+ const grid_volume &gv, void *P_internal_data) const;
+ virtual void subtract_P(field_type ft, realnum *f_minus_p[NUM_FIELD_COMPONENTS][2],
+ void *P_internal_data) const;
+
+ virtual int num_cinternal_notowned_needed(component c, void *P_internal_data) const;
+ virtual realnum *cinternal_notowned_ptr(int inotowned, component c, int cmp,
+ int n, void *P_internal_data) const;
+
+ virtual void dump_params(h5file *h5f, size_t *start);
+ virtual int get_num_params() { return 8; }
+ virtual bool needs_W_notowned(component c, realnum *W[NUM_FIELD_COMPONENTS][2]) const {
+ (void)c;
+ (void)W;
+ return true;
+ }
+
+protected:
+ double gyro_tensor[3][3];
+ double omega_0, gamma, alpha;
+ gyrotropy_model model;
+};
+
class multilevel_susceptibility : public susceptibility {
public:
multilevel_susceptibility() : L(0), T(0), Gamma(0), N0(0), alpha(0), omega(0), gamma(0) {}
diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
index 3c2ad82dd..a85a1250d 100644
--- a/src/meepgeom.cpp
+++ b/src/meepgeom.cpp
@@ -38,9 +38,11 @@ void check_offdiag(medium_struct *m) {
bool susceptibility_equal(const susceptibility &s1, const susceptibility &s2) {
return (vector3_equal(s1.sigma_diag, s2.sigma_diag) &&
- vector3_equal(s1.sigma_offdiag, s2.sigma_offdiag) && s1.frequency == s2.frequency &&
- s1.gamma == s2.gamma && s1.noise_amp == s2.noise_amp && s1.drude == s2.drude &&
- s1.is_file == s2.is_file);
+ vector3_equal(s1.sigma_offdiag, s2.sigma_offdiag) &&
+ vector3_equal(s1.bias, s2.bias) && s1.frequency == s2.frequency &&
+ s1.gamma == s2.gamma && s1.alpha == s2.alpha &&
+ s1.noise_amp == s2.noise_amp && s1.drude == s2.drude &&
+ s1.saturated_gyrotropy == s2.saturated_gyrotropy && s1.is_file == s2.is_file);
}
bool susceptibility_list_equal(const susceptibility_list &s1, const susceptibility_list &s2) {
@@ -1178,10 +1180,12 @@ double geom_epsilon::conductivity(meep::component c, const meep::vec &r) {
(must be updated manually, re-copying from ctl-io.cpp), if we
add new susceptibility subclasses) */
static bool susceptibility_equiv(const susceptibility *o0, const susceptibility *o) {
+ if (!vector3_equal(o0->bias, o->bias)) return 0;
if (o0->frequency != o->frequency) return 0;
if (o0->gamma != o->gamma) return 0;
if (o0->noise_amp != o->noise_amp) return 0;
if (o0->drude != o->drude) return 0;
+ if (o0->saturated_gyrotropy != o->saturated_gyrotropy) return 0;
if (o0->is_file != o->is_file) return 0;
if (o0->transitions != o->transitions) return 0;
@@ -1363,6 +1367,8 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s)
susceptibility *ss = &(p->user_s);
if (ss->is_file) meep::abort("unknown susceptibility");
bool noisy = (ss->noise_amp != 0.0);
+ bool gyrotropic = (ss->saturated_gyrotropy || ss->bias.x != 0.0
+ || ss->bias.y != 0.0 || ss->bias.z != 0.0);
meep::susceptibility *sus;
if (ss->transitions.size() != 0 || ss->initial_populations.size() != 0) {
@@ -1373,12 +1379,24 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s)
if (noisy) {
sus = new meep::noisy_lorentzian_susceptibility(ss->noise_amp, ss->frequency, ss->gamma,
ss->drude);
+ } else if (gyrotropic) {
+ meep::gyrotropy_model model
+ = ss->saturated_gyrotropy ? meep::GYROTROPIC_SATURATED
+ : ss->drude ? meep::GYROTROPIC_DRUDE : meep::GYROTROPIC_LORENTZIAN;
+ sus = new meep::gyrotropic_susceptibility(meep::vec(ss->bias.x, ss->bias.y, ss->bias.z),
+ ss->frequency, ss->gamma, ss->alpha, model);
} else {
sus = new meep::lorentzian_susceptibility(ss->frequency, ss->gamma, ss->drude);
}
- master_printf("%s%s susceptibility: frequency=%g, gamma=%g", noisy ? "noisy " : "",
- ss->drude ? "drude" : "lorentzian", ss->frequency, ss->gamma);
- if (noisy) master_printf(" amp=%g ", ss->noise_amp);
+ master_printf("%s%s susceptibility: frequency=%g, gamma=%g",
+ noisy ? "noisy " : gyrotropic ? "gyrotropic " : "",
+ ss->saturated_gyrotropy ? "Landau-Lifshitz-Gilbert-type"
+ : ss->drude ? "drude" : "lorentzian", ss->frequency, ss->gamma);
+ if (noisy) master_printf(", amp=%g ", ss->noise_amp);
+ if (gyrotropic) {
+ if (ss->saturated_gyrotropy) master_printf(", alpha=%g", ss->alpha);
+ master_printf(", bias=(%g,%g,%g)", ss->bias.x, ss->bias.y, ss->bias.z);
+ }
master_printf("\n");
}
diff --git a/src/structure_dump.cpp b/src/structure_dump.cpp
index b7cd91eff..7b9ca14e0 100644
--- a/src/structure_dump.cpp
+++ b/src/structure_dump.cpp
@@ -362,6 +362,32 @@ susceptibility *make_sus_list_from_params(h5file *file, int rank, size_t dims[3]
res = sus;
}
if (sus->next) sus = sus->next;
+ } else if (num_params == 8) {
+ // This is a gyrotropic_susceptibility and the next 8 values in the dataset
+ // are id, bias.x, bias.y, bias.z, omega_0, gamma, alpha, and model.
+ size_t gyro_susc_dims[3] = {8, 0, 0};
+ realnum gyro_susc_params[8];
+ file->read_chunk(rank, &start, gyro_susc_dims, gyro_susc_params);
+ start += gyro_susc_dims[0];
+
+ int id = (int)gyro_susc_params[0];
+ const vec bias(gyro_susc_params[1], gyro_susc_params[2], gyro_susc_params[3]);
+ const double omega_0 = gyro_susc_params[4];
+ const double gamma = gyro_susc_params[5];
+ const double alpha = gyro_susc_params[6];
+ const gyrotropy_model model = (gyrotropy_model)gyro_susc_params[7];
+ if (sus) {
+ sus->next =
+ new gyrotropic_susceptibility(bias, omega_0, gamma, alpha, model);
+ sus->next->ntot = ntot;
+ sus->next->set_id(id);
+ } else {
+ sus = new gyrotropic_susceptibility(bias, omega_0, gamma, alpha, model);
+ sus->ntot = ntot;
+ sus->set_id(id);
+ res = sus;
+ }
+ if (sus->next) sus = sus->next;
} else {
abort("Invalid number of susceptibility parameters in structure::load");
}
diff --git a/src/susceptibility.cpp b/src/susceptibility.cpp
index d45034352..f63332095 100644
--- a/src/susceptibility.cpp
+++ b/src/susceptibility.cpp
@@ -321,4 +321,266 @@ void noisy_lorentzian_susceptibility::dump_params(h5file *h5f, size_t *start) {
h5f->write_chunk(1, start, params_dims, params_data);
*start += num_params;
}
+
+
+gyrotropic_susceptibility::gyrotropic_susceptibility(const vec &bias, double omega_0, double gamma,
+ double alpha, gyrotropy_model model)
+ : omega_0(omega_0), gamma(gamma), alpha(alpha), model(model) {
+ // Precalculate g_{ij} = sum_k epsilon_{ijk} b_k, used in update_P.
+ // Ignore |b| for Landau-Lifshitz-Gilbert gyrotropy model.
+ const vec b = (model == GYROTROPIC_SATURATED) ? bias/abs(bias) : bias;
+ memset(gyro_tensor, 0, 9 * sizeof(double));
+ gyro_tensor[X][Y] = b.z(); gyro_tensor[Y][X] = -b.z();
+ gyro_tensor[Y][Z] = b.x(); gyro_tensor[Z][Y] = -b.x();
+ gyro_tensor[Z][X] = b.y(); gyro_tensor[X][Z] = -b.y();
+}
+
+/* To implement gyrotropic susceptibilities, we track three
+ polarization components (e.g. Px, Py, Pz) on EACH of the Yee cell's
+ three driving field positions (e.g., Ex, Ey, and Ez), i.e. 9
+ numbers per cell. This takes 3x the memory and runtime compared to
+ Lorentzian susceptibility. The advantage is that during update_P,
+ we can directly access the value of P at each update point without
+ averaging. */
+
+typedef struct {
+ size_t sz_data;
+ size_t ntot;
+ realnum *P[NUM_FIELD_COMPONENTS][2][3];
+ realnum *P_prev[NUM_FIELD_COMPONENTS][2][3];
+ realnum data[1];
+} gyrotropy_data;
+
+void *gyrotropic_susceptibility::new_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2],
+ const grid_volume &gv) const {
+ int num = 0;
+ FOR_COMPONENTS(c) DOCMP2 {
+ if (needs_P(c, cmp, W)) num += 6 * gv.ntot();
+ }
+ size_t sz = sizeof(gyrotropy_data) + sizeof(realnum) * (num - 1);
+ gyrotropy_data *d = (gyrotropy_data *)malloc(sz);
+ d->sz_data = sz;
+ return (void *)d;
+}
+
+void gyrotropic_susceptibility::init_internal_data(realnum *W[NUM_FIELD_COMPONENTS][2], double dt,
+ const grid_volume &gv, void *data) const {
+ (void)dt; // unused
+ gyrotropy_data *d = (gyrotropy_data *)data;
+ size_t sz_data = d->sz_data;
+ memset(d, 0, sz_data);
+ d->sz_data = sz_data;
+ d->ntot = gv.ntot();
+ realnum *p = d->data;
+ FOR_COMPONENTS(c) DOCMP2 {
+ if (needs_P(c, cmp, W)) {
+ for (int dd = X; dd < R; dd++) {
+ d->P[c][cmp][dd] = p; p += d->ntot;
+ d->P_prev[c][cmp][dd] = p; p += d->ntot;
+ }
+ }
+ }
+}
+
+void *gyrotropic_susceptibility::copy_internal_data(void *data) const {
+ gyrotropy_data *d = (gyrotropy_data *)data;
+ if (!d) return 0;
+ gyrotropy_data *dnew = (gyrotropy_data *)malloc(d->sz_data);
+ memcpy(dnew, d, d->sz_data);
+ realnum *p = dnew->data;
+ FOR_COMPONENTS(c) DOCMP2 {
+ if (d->P[c][cmp][0]) {
+ for (int dd = X; dd < R; dd++) {
+ dnew->P[c][cmp][dd] = p; p += d->ntot;
+ dnew->P_prev[c][cmp][dd] = p; p += d->ntot;
+ }
+ }
+ }
+ return (void *)dnew;
+}
+
+bool gyrotropic_susceptibility::needs_P(component c, int cmp, realnum *W[NUM_FIELD_COMPONENTS][2]) const {
+ if (!is_electric(c) && !is_magnetic(c)) return false;
+ direction d0 = component_direction(c);
+ return (d0 == X || d0 == Y || d0 == Z) && sigma[c][d0] && W[c][cmp];
+}
+
+// Similar to the OFFDIAG macro, but without averaging sigma.
+#define OFFDIAGW(g, sx, s) (0.25 * (g[i] + g[i - sx] + g[i + s] + g[i + s - sx]))
+
+void gyrotropic_susceptibility::update_P(realnum *W[NUM_FIELD_COMPONENTS][2],
+ realnum *W_prev[NUM_FIELD_COMPONENTS][2], double dt,
+ const grid_volume &gv, void *P_internal_data) const {
+ gyrotropy_data *d = (gyrotropy_data *)P_internal_data;
+ const double omega2pidt = 2 * pi * omega_0 * dt;
+ const double g2pidt = 2 * pi * gamma * dt;
+ (void)W_prev; // unused;
+
+ switch (model) {
+ case GYROTROPIC_LORENTZIAN:
+ case GYROTROPIC_DRUDE:
+ {
+ const double omega0dtsqr = omega2pidt * omega2pidt;
+ const double gamma1 = (1 - g2pidt / 2);
+ const double diag = 2 - (model == GYROTROPIC_DRUDE ? 0 : omega0dtsqr);
+ const double pt = pi*dt;
+
+ // Precalculate 3x3 matrix inverse, exploiting skew symmetry
+ const double gd = (1 + g2pidt / 2);
+ const double gx = pt * gyro_tensor[Y][Z];
+ const double gy = pt * gyro_tensor[Z][X];
+ const double gz = pt * gyro_tensor[X][Y];
+ const double invdet = 1.0 / gd / (gd*gd + gx*gx + gy*gy + gz*gz);
+ const double inv[3][3]
+ = {{ invdet*(gd*gd+gx*gx), invdet*(gx*gy+gd*gz), invdet*(gx*gz-gd*gy) },
+ { invdet*(gy*gx-gd*gz), invdet*(gd*gd+gy*gy), invdet*(gy*gz+gd*gx) },
+ { invdet*(gz*gx+gd*gy), invdet*(gz*gy-gd*gx), invdet*(gd*gd+gz*gz) }};
+
+ FOR_COMPONENTS(c) DOCMP2 {
+ if (d->P[c][cmp][0]) {
+ const direction d0 = component_direction(c);
+ const realnum *w0 = W[c][cmp], *s = sigma[c][d0];
+
+ if (!w0 || !s || (d0 != X && d0 != Y && d0 != Z))
+ abort("gyrotropic media require 3D Cartesian fields\n");
+
+ const direction d1 = cycle_direction(gv.dim, d0, 1);
+ const direction d2 = cycle_direction(gv.dim, d0, 2);
+ const realnum *w1 = W[direction_component(c, d1)][cmp];
+ const realnum *w2 = W[direction_component(c, d2)][cmp];
+ realnum *p0 = d->P[c][cmp][d0], *pp0 = d->P_prev[c][cmp][d0];
+ realnum *p1 = d->P[c][cmp][d1], *pp1 = d->P_prev[c][cmp][d1];
+ realnum *p2 = d->P[c][cmp][d2], *pp2 = d->P_prev[c][cmp][d2];
+ const ptrdiff_t is = gv.stride(d0) * (is_magnetic(c) ? -1 : +1);
+ const ptrdiff_t is1 = gv.stride(d1) * (is_magnetic(c) ? -1 : +1);
+ const ptrdiff_t is2 = gv.stride(d2) * (is_magnetic(c) ? -1 : +1);
+ realnum r0, r1, r2;
+
+ if (!pp1 || !pp2)
+ abort("gyrotropic media require 3D Cartesian fields\n");
+ if (sigma[c][d1] || sigma[c][d2])
+ abort("gyrotropic media do not support anisotropic sigma\n");
+
+ LOOP_OVER_VOL_OWNED(gv, c, i) {
+ r0 = diag*p0[i] - gamma1*pp0[i] + omega0dtsqr*s[i]*w0[i]
+ - pt*gyro_tensor[d0][d1]*pp1[i] - pt*gyro_tensor[d0][d2]*pp2[i];
+ r1 = diag*p1[i] - gamma1*pp1[i] + (w1 ? omega0dtsqr*s[i]*OFFDIAGW(w1,is1,is) : 0)
+ - pt*gyro_tensor[d1][d0]*pp0[i] - pt*gyro_tensor[d1][d2]*pp2[i];
+ r2 = diag*p2[i] - gamma1*pp2[i] + (w2 ? omega0dtsqr*s[i]*OFFDIAGW(w2,is2,is) : 0)
+ - pt*gyro_tensor[d2][d1]*pp1[i] - pt*gyro_tensor[d2][d0]*pp0[i];
+
+ pp0[i] = p0[i]; pp1[i] = p1[i]; pp2[i] = p2[i];
+ p0[i] = inv[d0][d0] * r0 + inv[d0][d1] * r1 + inv[d0][d2] * r2;
+ p1[i] = inv[d1][d0] * r0 + inv[d1][d1] * r1 + inv[d1][d2] * r2;
+ p2[i] = inv[d2][d0] * r0 + inv[d2][d1] * r1 + inv[d2][d2] * r2;
+ }
+ }
+ }
+ }
+ break;
+
+ case GYROTROPIC_SATURATED:
+ {
+ const double dt2pi = 2*pi*dt;
+
+ // Precalculate 3x3 matrix inverse, exploiting skew symmetry
+ const double gd = 0.5;
+ const double gx = -0.5 * alpha * gyro_tensor[Y][Z];
+ const double gy = -0.5 * alpha * gyro_tensor[Z][X];
+ const double gz = -0.5 * alpha * gyro_tensor[X][Y];
+ const double invdet = 1.0 / gd / (gd*gd + gx*gx + gy*gy + gz*gz);
+ const double inv[3][3]
+ = {{ invdet*(gd*gd+gx*gx), invdet*(gx*gy+gd*gz), invdet*(gx*gz-gd*gy) },
+ { invdet*(gy*gx-gd*gz), invdet*(gd*gd+gy*gy), invdet*(gy*gz+gd*gx) },
+ { invdet*(gz*gx+gd*gy), invdet*(gz*gy-gd*gx), invdet*(gd*gd+gz*gz) }};
+
+ FOR_COMPONENTS(c) DOCMP2 {
+ if (d->P[c][cmp][0]) {
+ const direction d0 = component_direction(c);
+ const realnum *w0 = W[c][cmp], *s = sigma[c][d0];
+
+ if (!w0 || !s || (d0 != X && d0 != Y && d0 != Z))
+ abort("gyrotropic media require 3D Cartesian fields\n");
+
+ const direction d1 = cycle_direction(gv.dim, d0, 1);
+ const direction d2 = cycle_direction(gv.dim, d0, 2);
+ const realnum *w1 = W[direction_component(c, d1)][cmp];
+ const realnum *w2 = W[direction_component(c, d2)][cmp];
+ realnum *p0 = d->P[c][cmp][d0], *pp0 = d->P_prev[c][cmp][d0];
+ realnum *p1 = d->P[c][cmp][d1], *pp1 = d->P_prev[c][cmp][d1];
+ realnum *p2 = d->P[c][cmp][d2], *pp2 = d->P_prev[c][cmp][d2];
+ const ptrdiff_t is = gv.stride(d0) * (is_magnetic(c) ? -1 : +1);
+ const ptrdiff_t is1 = gv.stride(d1) * (is_magnetic(c) ? -1 : +1);
+ const ptrdiff_t is2 = gv.stride(d2) * (is_magnetic(c) ? -1 : +1);
+ realnum r0, r1, r2, q0, q1, q2;
+
+ if (!pp1 || !pp2)
+ abort("gyrotropic media require 3D Cartesian fields\n");
+ if (sigma[c][d1] || sigma[c][d2])
+ abort("gyrotropic media do not support anisotropic sigma\n");
+
+ LOOP_OVER_VOL_OWNED(gv, c, i) {
+ q0 = -omega2pidt*p0[i] + 0.5*alpha*pp0[i] + dt2pi*s[i]*w0[i];
+ q1 = -omega2pidt*p1[i] + 0.5*alpha*pp1[i] + dt2pi*s[i]*(w1 ? OFFDIAGW(w1,is1,is) : 0);
+ q2 = -omega2pidt*p2[i] + 0.5*alpha*pp2[i] + dt2pi*s[i]*(w2 ? OFFDIAGW(w2,is2,is) : 0);
+
+ r0 = 0.5*pp0[i] - g2pidt*p0[i] + gyro_tensor[d0][d1]*q1 + gyro_tensor[d0][d2]*q2;
+ r1 = 0.5*pp1[i] - g2pidt*p1[i] + gyro_tensor[d1][d2]*q2 + gyro_tensor[d1][d0]*q0;
+ r2 = 0.5*pp2[i] - g2pidt*p2[i] + gyro_tensor[d2][d0]*q0 + gyro_tensor[d2][d1]*q1;
+
+ pp0[i] = p0[i]; pp1[i] = p1[i]; pp2[i] = p2[i];
+ p0[i] = inv[d0][d0] * r0 + inv[d0][d1] * r1 + inv[d0][d2] * r2;
+ p1[i] = inv[d1][d0] * r0 + inv[d1][d1] * r1 + inv[d1][d2] * r2;
+ p2[i] = inv[d2][d0] * r0 + inv[d2][d1] * r1 + inv[d2][d2] * r2;
+ }
+ }
+ }
+ }
+ break;
+ }
+}
+
+void gyrotropic_susceptibility::subtract_P(field_type ft,
+ realnum *f_minus_p[NUM_FIELD_COMPONENTS][2],
+ void *P_internal_data) const {
+ gyrotropy_data *d = (gyrotropy_data *)P_internal_data;
+ field_type ft2 = ft == E_stuff ? D_stuff : B_stuff; // for sources etc.
+ size_t ntot = d->ntot;
+ FOR_FT_COMPONENTS(ft, ec) DOCMP2 {
+ if (d->P[ec][cmp][0]) {
+ component dc = field_type_component(ft2, ec);
+ if (f_minus_p[dc][cmp]) {
+ realnum *p = d->P[ec][cmp][component_direction(ec)];
+ realnum *fmp = f_minus_p[dc][cmp];
+ for (size_t i = 0; i < ntot; ++i)
+ fmp[i] -= p[i];
+ }
+ }
+ }
+}
+
+int gyrotropic_susceptibility::num_cinternal_notowned_needed(component c,
+ void *P_internal_data) const {
+ gyrotropy_data *d = (gyrotropy_data *)P_internal_data;
+ return d->P[c][0][0] ? 3 : 0;
+}
+
+realnum *gyrotropic_susceptibility::cinternal_notowned_ptr(int inotowned, component c, int cmp,
+ int n, void *P_internal_data) const {
+ gyrotropy_data *d = (gyrotropy_data *)P_internal_data;
+ if (!d || !d->P[c][cmp][inotowned]) return NULL;
+ return d->P[c][cmp][inotowned] + n;
+}
+
+void gyrotropic_susceptibility::dump_params(h5file *h5f, size_t *start) {
+ size_t num_params = 9;
+ size_t params_dims[1] = {num_params};
+ double bias[] = { gyro_tensor[Y][Z], gyro_tensor[Z][X], gyro_tensor[X][Y] };
+ double params_data[] = {
+ 8, (double)get_id(), bias[X], bias[Y], bias[Z],
+ omega_0, gamma, alpha, (double)model};
+ h5f->write_chunk(1, start, params_dims, params_data);
+ *start += num_params;
+}
+
} // namespace meep