Python interface for Gaussian beam source

python/meep.i

@@ -1029,6 +1029,17 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
 }
 
 
+// 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] {
+    std::complex<double> *pyE0  = (std::complex<double> *)array_data($input);
+    $1[0] = pyE0[0]; $1[1] = pyE0[1]; $1[2] = pyE0[2];
+}
+
 //--------------------------------------------------
 // typemaps needed for get_eigenmode_coefficients
 //--------------------------------------------------
@@ -1648,6 +1659,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,
+                [src.beam_E0.x, src.beam_E0.y, src.beam_E0.z]
+            ]
+            gaussianbeam = functools.partial(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,39 @@ 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):
+
+    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):
+
+        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