Python interface for Gaussian beam source (#1310)

* Python interface for Gaussian beam source

* create gaussianbeam object in Python using functools.partial

* rename gaussianbeam constructor argument EO to E0

* pass E0 from Python as NumPy array

* add documentation, tutorial example, and test (Python)

* Update Python_User_Interface.md

* Update Python_User_Interface.md

Co-authored-by: Steven G. Johnson <[email protected]>

doc/docs/FAQ.md

@@ -142,15 +142,13 @@ For an example of the second approach, see [Tutorial/Eigenmode Source/Planewaves
 
 ### How do I create a focused beam with a Gaussian envelope?
 
-A focused beam with a Gaussian envelope can be created using the amplitude function (`amp_func`) of the [`Source`](Python_User_Interface.md#source) object. Examples are provided for [Python](https://github.com/NanoComp/meep/blob/master/python/examples/gaussian-beam.py) and [Scheme](https://github.com/NanoComp/meep/blob/master/scheme/examples/gaussian-beam.ctl). Four snapshots of the resulting field profile generated using this script for different values of the beam width (`sigma`) and rotation angle (`tilt_angle`) are shown in the following image:
+A focused beam with a Gaussian envelope can be created using the [`GaussianBeamSource`](Python_User_Interface.md#gaussianbeamsource) as demonstrated in this [script](https://github.com/NanoComp/meep/blob/master/python/examples/gaussian-beam.py). The following are four snapshots of the resulting field profile generated using a line source (red line) in 2d for different values of the beam waist radius/width $w_0$ (`beam_w0`), propagation direction $\vec{k}$ (`beam_kdir`), and focus (`beam_x0`):
 
 <center>
 ![](images/gaussian_beam.png)
 </center>
 
-Beams in a homogeneous material do not have a fixed width in Maxwell's equations; they always spread out during propagation. The [numerical aperture (NA)](https://en.wikipedia.org/wiki/Gaussian_beam#Beam_divergence) of a Gaussian beam of width w (2*`sigma` from the example script) and vacuum wavelength $\lambda$ in a medium of index $n$ is $n\sin(\lambda/(\pi nw))$.
-
-Note: in this example, the beam waist is at the source position (i.e., top center of the cell). If you want the beam waist to be at a position other than the position of the source, you need to adjust the *phase* of the beam accordingly. If you assume you have a Gaussian beam profile with zero phase at some plane $y=y_0$, then you can work out the beam profile (including phase) at any other plane $y=y_1$ by taking the Fourier transform and looking at the propagation of each planewave component, and then inverse Fourier transforming. In this way, you can work out the desired source profile at any plane $y=y_1$ to get a Gaussian beam waist at $y=y_0$.
+Note: beams in a homogeneous material do *not* have a fixed width in Maxwell's equations; they always spread out during propagation (i.e., diffraction). The [numerical aperture (NA)](https://en.wikipedia.org/wiki/Gaussian_beam#Beam_divergence) of a Gaussian beam with waist radius $w$ and vacuum wavelength $\lambda$ within a medium of index $n$ is $n\sin(\lambda/(\pi nw))$.
 
 ### How do I create a circularly polarized planewave source in cylindrical coordinates?
 

doc/docs/Python_User_Interface.md

@@ -5832,7 +5832,7 @@ it does support anisotropic ε. Any nonlinearities, magnetic responses μ,
 conductivities σ, or dispersive polarizations in your materials will be *ignored* when
 computing the eigenmode source. PML will also be ignored.
 
-The `src_time` object (`Source.src`), which specifies the time dependence of the
+The `SourceTime` object (`Source.src`), which specifies the time dependence of the
 source, can be one of `ContinuousSource`, `GaussianSource` or `CustomSource`.
 
 </div>
@@ -5961,6 +5961,59 @@ Returns the total power of the fields from the eigenmode source at frequency `fr
 
 </div>
 
+
+---
+<a id="GaussianBeamSource"></a>
+
+### GaussianBeamSource
+
+```python
+class GaussianBeamSource(Source):
+```
+
+<div class="class_docstring" markdown="1">
+
+This is a subclass of `Source` and has **all of the properties** of `Source` above. However, the `component` parameter of the `Source` object is ignored. The [Gaussian beam](https://en.wikipedia.org/wiki/Gaussian_beam) is a transverse electromagnetic mode for which the source region must be a *line* (in 2d) or *plane* (in 3d). For a beam polarized in the $x$ direction with propagation along $+z$, the electric field is defined by $\mathbf{E}(r,z)=E_0\hat{x}\frac{w_0}{w(z)}\exp\left(\frac{-r^2}{w(z)^2}\right)\exp\left(-i\left(kz + k\frac{r^2}{2R(z)}\right)\right)$ where $r$ is the radial distance from the center axis of the beam, $z$ is the axial distance from the beam's focus (or "waist"), $k=2\pi n/\lambda$ is the wavenumber (for a free-space wavelength $\lambda$ and refractive index $n$ of the homogeneous, lossless medium in which the beam propagates), $E_0$ is the electric field amplitude at the origin, $w(z)$ is the radius at which the field amplitude decays by $1/e$ of its axial values, $w_0$ is the beam waist radius, and $R(z)$ is the radius of curvature of the beam's wavefront at $z$. The only independent parameters that need to be specified are $w_0$, $E_0$, $k$, and the location of the beam focus (i.e., the origin: $r=z=0$).
+
+(In 3d, we use a ["complex point-source" method](https://doi.org/10.1364/JOSAA.16.001381) to define a source that generates an exact Gaussian-beam solution.  In 2d, we currently use the simple approximation of taking a cross-section of the 3d beam.  In both cases, the beam is most accurate near the source's center frequency.)
+
+The `SourceTime` object (`Source.src`), which specifies the time dependence of the source, should normally be a narrow-band `ContinuousSource` or `GaussianSource`.  (For a `CustomSource`, the beam frequency is determined by the source's `center_frequency` parameter.)
+
+</div>
+
+<a id="GaussianBeamSource.__init__"></a>
+
+<div class="class_members" markdown="1">
+
+```python
+def __init__(self,
+             src,
+             center=None,
+             volume=None,
+             component=mp.ALL_COMPONENTS,
+             beam_x0=Vector3(),
+             beam_kdir=Vector3(),
+             beam_w0=None,
+             beam_E0=Vector3(),
+             **kwargs):
+```
+
+<div class="method_docstring" markdown="1">
+
+Construct a `GaussianBeamSource`.
+
++ **`beam_x0` [`Vector3`]** — The location of the beam focus *relative* to the center of the source. This does *not* need to lie within the source region (i.e., the beam focus can be anywhere, inside or outside the cell, independent of position of the source).
+
++ **`beam_kdir` [`Vector3`]** — The propagation direction of the beam. The length is *ignored*. The wavelength of the beam is determined by the center frequency of the `Source.src` object and the refractive index of the homogeneous medium by its value at the center of the source region.
+
++ **`beam_w0` [`number`]** — The beam waist radius.
+
++ **`beam_E0` [`Vector3`]** — The polarization vector of the beam. Elements can be complex valued (i.e., for circular polarization). The polarization vector must be *parallel* to the source region in order to generate a transverse mode.
+
+</div>
+
+</div>
+
 ---
 <a id="ContinuousSource"></a>
 

python/Makefile.am

@@ -52,6 +52,7 @@ TESTS =                                   \
     $(TEST_DIR)/field_functions.py        \
     $(TEST_DIR)/force.py                  \
     $(TEST_DIR)/fragment_stats.py         \
+    $(TEST_DIR)/gaussianbeam.py           \
     $(TEST_DIR)/geom.py                   \
     $(TEST_DIR)/get_point.py              \
     $(TEST_DIR)/holey_wvg_bands.py        \

python/examples/gaussian-beam.py

@@ -1,48 +1,42 @@
-## generate a titled Gaussian beam profile by defining the amplitude function of the source
+## launch a Gaussian beam
 
 import meep as mp
 import math
-import cmath
-import numpy as np
+import matplotlib
+matplotlib.use('agg')
 import matplotlib.pyplot as plt
 
-resolution = 40 # pixels/μm
+s = 14
+resolution = 50
+dpml = 2
 
-cell_size = mp.Vector3(20,10,0)
+cell_size = mp.Vector3(s,s)
 
-pml_layers = [mp.PML(thickness=1.0,direction=mp.Y)]
+boundary_layers = [mp.PML(thickness=dpml)]
 
-fcen = 1.0 # center frequency of CW source (wavelength is 1 μm)
+beam_x0 = mp.Vector3(0,3.0)    # beam focus (relative to source center)
+rot_angle = 0  # CCW rotation angle about z axis (0: +y axis)
+beam_kdir = mp.Vector3(0,1,0).rotate(mp.Vector3(0,0,1),math.radians(rot_angle))  # beam propagation direction
+beam_w0 = 0.8  # beam waist radius
+beam_E0 = mp.Vector3(0,0,1)
+fcen = 1
+sources = [mp.GaussianBeamSource(src=mp.ContinuousSource(fcen),
+                                 center=mp.Vector3(0,-0.5*s+dpml+1.0),
+                                 size=mp.Vector3(s),
+                                 beam_x0=beam_x0,
+                                 beam_kdir=beam_kdir,
+                                 beam_w0=beam_w0,
+                                 beam_E0=beam_E0)]
 
-tilt_angle = math.radians(-10) # angle of tilted beam
-k = mp.Vector3(y=1).rotate(mp.Vector3(z=1),tilt_angle).scale(fcen)
+sim = mp.Simulation(resolution=resolution,
+                    cell_size=cell_size,
+                    boundary_layers=boundary_layers,
+                    sources=sources)
 
-sigma = 1.5 # beam width
+sim.run(until=20)
 
-def gaussian_beam(sigma, k, x0):
-    def _gaussian_beam(x):
-        return cmath.exp(1j*2*math.pi*k.dot(x-x0)-(x-x0).dot(x-x0)/(2*sigma**2))
-    return _gaussian_beam
+sim.plot2D(fields=mp.Ez,
+           output_plane=mp.Volume(center=mp.Vector3(),
+                                  size=mp.Vector3(s-2*dpml,s-2*dpml)))
 
-src_pt = mp.Vector3(y=4)
-sources = [mp.Source(src=mp.ContinuousSource(fcen, fwidth=0.2*fcen),
-                     component=mp.Ez,
-                     center=src_pt,
-                     size=mp.Vector3(20),
-                     amp_func=gaussian_beam(sigma,k,src_pt))]
-
-sim = mp.Simulation(cell_size=cell_size,
-                    sources=sources,
-                    k_point=k,
-                    boundary_layers=pml_layers,
-                    resolution=resolution)
-
-non_pml_vol = mp.Volume(center=mp.Vector3(), size=mp.Vector3(20,8,0))
-sim.run(until=50)
-
-ez_data = sim.get_array(vol=non_pml_vol, component=mp.Ez)
-
-plt.figure()
-plt.imshow(np.flipud(np.transpose(np.real(ez_data))), interpolation='spline36', cmap='RdBu')
-plt.axis('off')
-plt.show()
+plt.savefig('Ez_angle{}.png'.format(rot_angle),bbox_inches='tight',pad_inches=0)

python/meep.i

@@ -1028,6 +1028,15 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
     $1 = (std::complex<double> *)array_data($input);
 }
 
+// typemaps for gaussianbeam
+
+%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") std::complex<double> E0[3] {
+    $1 = is_array($input);
+}
+
+%typemap(in) std::complex<double> E0[3] {
+     $1 = (std::complex<double> *)array_data($input);
+}
 
 //--------------------------------------------------
 // typemaps needed for get_eigenmode_coefficients
@@ -1648,6 +1657,7 @@ PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where
         Source,
         SourceTime,
         check_positive,
+        GaussianBeamSource,
     )
     from .visualization import (
         plot2D,

python/simulation.py

@@ -18,7 +18,7 @@
 
 import meep as mp
 from meep.geom import Vector3, init_do_averaging
-from meep.source import EigenModeSource, check_positive
+from meep.source import EigenModeSource, GaussianBeamSource, check_positive
 import meep.visualization as vis
 from meep.verbosity_mgr import Verbosity
 
@@ -1997,14 +1997,27 @@ def get_epsilon_point(self, pt, frequency=0, omega=0):
         Given a frequency `frequency` and a `Vector3` `pt`, returns the average eigenvalue
         of the permittivity tensor at that location and frequency. If `frequency` is
         non-zero, the result is complex valued; otherwise it is the real,
-        frequency-independent part of ε (the $\\omega\\to\\infty$ limit).
+        frequency-independent part of $\\varepsilon$ (the $\\omega\\to\\infty$ limit).
         """
         v3 = py_v3_to_vec(self.dimensions, pt, self.is_cylindrical)
         if omega != 0:
             frequency = omega
             warnings.warn("get_epsilon_point: omega has been deprecated; use frequency instead", RuntimeWarning)
         return self.fields.get_eps(v3,frequency)
 
+    def get_mu_point(self, pt, frequency=0, omega=0):
+        """
+        Given a frequency `frequency` and a `Vector3` `pt`, returns the average eigenvalue
+        of the permeability tensor at that location and frequency. If `frequency` is
+        non-zero, the result is complex valued; otherwise it is the real,
+        frequency-independent part of $\\mu$ (the $\\omega\\to\\infty$ limit).
+        """
+        v3 = py_v3_to_vec(self.dimensions, pt, self.is_cylindrical)
+        if omega != 0:
+            frequency = omega
+            warnings.warn("get_mu_point: omega has been deprecated; use frequency instead", RuntimeWarning)
+        return self.fields.get_mu(v3,frequency)
+
     def get_filename_prefix(self):
         """
         Return the current prefix string that is prepended, by default, to all file names.
@@ -2234,6 +2247,20 @@ def add_source(self, src):
                 add_eig_src()
             else:
                 add_eig_src(src.amp_func)
+        elif isinstance (src, GaussianBeamSource):
+            gaussianbeam_args = [
+                py_v3_to_vec(self.dimensions, src.beam_x0, is_cylindrical=self.is_cylindrical),
+                py_v3_to_vec(self.dimensions, src.beam_kdir, is_cylindrical=self.is_cylindrical),
+                src.beam_w0,
+                src.src.swigobj.frequency().real,
+                self.fields.get_eps(py_v3_to_vec(self.dimensions, src.center, self.is_cylindrical)).real,
+                self.fields.get_mu(py_v3_to_vec(self.dimensions, src.center, self.is_cylindrical)).real,
+                np.array([src.beam_E0.x, src.beam_E0.y, src.beam_E0.z], dtype=np.complex128)
+            ]
+            gaussianbeam = mp.gaussianbeam(*gaussianbeam_args)
+            add_vol_src_args = [src.src.swigobj, where, gaussianbeam]
+            add_vol_src = functools.partial(self.fields.add_volume_source, *add_vol_src_args)
+            add_vol_src()
         else:
             add_vol_src_args = [src.component, src.src.swigobj, where]
             add_vol_src = functools.partial(self.fields.add_volume_source, *add_vol_src_args)

python/source.py

@@ -526,3 +526,57 @@ def eig_power(self,freq):
         if callable(getattr(self.src, "fourier_transform", None)):
            amp *= self.src.fourier_transform(freq)
         return abs(amp)**2
+
+
+class GaussianBeamSource(Source):
+    """
+    This is a subclass of `Source` and has **all of the properties** of `Source` above.
+
+    The [Gaussian beam](https://en.wikipedia.org/wiki/Gaussian_beam) is a transverse electromagnetic mode for which the source region *must* be a line in 2d or plane in 3d. For a beam polarized in the $x$ direction with propagation along $+z$, the electric field is defined by $\\mathbf{E}(r,z)=E_0\\hat{x}\\frac{w_0}{w(z)}\\exp\\left(\\frac{-r^2}{w(z)^2}\\right)\\exp\\left(-i\\left(kz + k\\frac{r^2}{2R(z)}\\right)\\right)$ where $r$ is the radial distance from the center axis of the beam, $z$ is the axial distance from the beam's focus (or "waist"), $k=2\pi n/\lambda$ is the wavenumber (for a free-space wavelength $\lambda$ and refractive index $n$ of the homogeneous, lossless medium in which the beam propagates), $E_0$ is the electric field amplitude at the origin, $w(z)$ is the radius at which the field amplitude decays by $1/e$ of its axial values, $w_0$ is the beam waist radius, and $R(z)$ is the radius of curvature of the beam's wavefront at $z$. In Meep, the only parameters that need to be specified are $w_0$, $E_0$, $k$, and the location of the beam focus (i.e., the origin $r=z=0$ from the previous formula). The `component` parameter of the `Source` object is ignored.
+
+    The `src_time` object (`Source.src`), which specifies the time dependence of the source, can be one of `ContinuousSource` or `GaussianSource`.
+    """
+
+    def __init__(self,
+                 src,
+                 center=None,
+                 volume=None,
+                 component=mp.ALL_COMPONENTS,
+                 beam_x0=Vector3(),
+                 beam_kdir=Vector3(),
+                 beam_w0=None,
+                 beam_E0=Vector3(),
+                 **kwargs):
+        """
+        Construct an `GaussianBeamSource`.
+
+        + **`beam_x0` [`Vector3`]** — The location of the beam focus. This does *not* need to lie within the source region (i.e., the beam can be focused at any arbitrary point, inside or outside the cell, independent of where the source is positioned).
+
+        + **`beam_kdir` [`Vector3`]** — The propagation direction of the beam. The length is *ignored*. (The wavelength of the beam is determined by the center frequency of the `src_time` object and the refractive index of the homogeneous medium by its value at the center of the source region.)
+
+        + **`beam_w0` [`number`]** — The beam waist radius.
+
+        + **`beam_E0` [`Vector3`]** — The polarization vector of the beam. Elements can be complex valued (i.e., for circular polarization). The polarization vector must be *parallel* to the source region (i.e., in order to generate a transverse mode).
+        """
+
+        super(GaussianBeamSource, self).__init__(src, component, center, volume, **kwargs)
+        self._beam_x0 = beam_x0
+        self._beam_kdir = beam_kdir
+        self._beam_w0 = beam_w0
+        self._beam_E0 = beam_E0
+
+    @property
+    def beam_x0(self):
+        return self._beam_x0
+
+    @property
+    def beam_kdir(self):
+        return self._beam_kdir
+
+    @property
+    def beam_w0(self):
+        return self._beam_w0
+
+    @property
+    def beam_E0(self):
+        return self._beam_E0

python/tests/gaussianbeam.py

@@ -0,0 +1,62 @@
+from __future__ import division
+
+import unittest
+import meep as mp
+import math
+import numpy as np
+
+class TestGaussianBeamSource(unittest.TestCase):
+
+    def gaussian_beam(self, rot_angle):
+
+        s = 14
+        resolution = 25
+        dpml = 2
+
+        cell_size = mp.Vector3(s,s)
+        boundary_layers = [mp.PML(thickness=dpml)]
+
+        beam_x0 = mp.Vector3(0,7.0)    # beam focus (relative to source center)
+        beam_kdir = mp.Vector3(0,1,0).rotate(mp.Vector3(0,0,1),math.radians(rot_angle))  # beam propagation direction
+        beam_w0 = 0.8  # beam waist radius
+        beam_E0 = mp.Vector3(0,0,1)
+        beam_x0 = beam_x0.rotate(mp.Vector3(0,0,1),math.radians(rot_angle))
+        fcen = 1
+        src_x = 0
+        src_y = -0.5*s+dpml+1.0
+        sources = [mp.GaussianBeamSource(src=mp.GaussianSource(fcen,fwidth=0.2*fcen),
+                                         center=mp.Vector3(src_x,src_y),
+                                         size=mp.Vector3(s),
+                                         beam_x0=beam_x0,
+                                         beam_kdir=beam_kdir,
+                                         beam_w0=beam_w0,
+                                         beam_E0=beam_E0)]
+
+        sim = mp.Simulation(resolution=resolution,
+                            cell_size=cell_size,
+                            boundary_layers=boundary_layers,
+                            sources=sources)
+
+        cell_dft_fields = sim.add_dft_fields([mp.Ez], fcen, 0, 1, center=mp.Vector3(), size=mp.Vector3(s-2*dpml,s-2*dpml))
+
+        sim.run(until_after_sources=50)
+
+        Ez_cell = sim.get_dft_array(cell_dft_fields, mp.Ez, 0)
+        [x,y,z,w] = sim.get_array_metadata(dft_cell=cell_dft_fields)
+
+        tol = 0.05
+        idx_x = np.nonzero((np.squeeze(x) > (src_x+beam_x0.x-tol)) & (np.squeeze(x) < (src_x+beam_x0.x+tol)))
+        idx_y = np.nonzero((np.squeeze(y) > (src_y+beam_x0.y-tol)) & (np.squeeze(y) < (src_y+beam_x0.y+tol)))
+
+        Ez_beam_x0 = Ez_cell[np.squeeze(idx_x)[0],np.squeeze(idx_y)[0]]
+
+        frac = np.abs(Ez_beam_x0)**2 / np.amax(np.abs(Ez_cell)**2)
+        print("ratio of the Gaussian beam energy at the focus over the maximum beam energy for the entire cell: {}".format(frac))
+
+        self.assertGreater(frac, 0.99)
+
+    def test_gaussian_beam(self):
+        self.gaussian_beam(-40)
+
+if __name__ == '__main__':
+  unittest.main()

src/meep.hpp

@@ -1488,7 +1488,7 @@ class gaussianbeam {
 
 public:
   gaussianbeam(const vec &x0, const vec &kdir, double w0, double freq,
-               double eps, double mu, std::complex<double> EO[3]);
+               double eps, double mu, std::complex<double> E0[3]);
   void get_fields(std::complex<double> *EH, const vec &x) const;
   std::complex<double> get_E0(int n) const { return E0[n]; };