diff --git a/doc/docs/FAQ.md b/doc/docs/FAQ.md
index 173d063bb..fe52c91cb 100644
--- a/doc/docs/FAQ.md
+++ b/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`):
![](images/gaussian_beam.png)
-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?
diff --git a/doc/docs/Python_User_Interface.md b/doc/docs/Python_User_Interface.md
index be4e2f6d4..d48dac143 100644
--- a/doc/docs/Python_User_Interface.md
+++ b/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`.
@@ -5961,6 +5961,59 @@ Returns the total power of the fields from the eigenmode source at frequency `fr
+
+---
+
+
+### GaussianBeamSource
+
+```python
+class GaussianBeamSource(Source):
+```
+
+
+
+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.)
+
+
+
+
+
+
+
+```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):
+```
+
+
+
+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.
+
+
+
+
+
---
diff --git a/doc/docs/images/gaussian_beam.png b/doc/docs/images/gaussian_beam.png
index b170c6bc0..4672df00e 100644
Binary files a/doc/docs/images/gaussian_beam.png and b/doc/docs/images/gaussian_beam.png differ
diff --git a/python/Makefile.am b/python/Makefile.am
index 2b213809d..f213d3c73 100644
--- a/python/Makefile.am
+++ b/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 \
diff --git a/python/examples/gaussian-beam.py b/python/examples/gaussian-beam.py
index 6a68c0f9e..bc2747718 100644
--- a/python/examples/gaussian-beam.py
+++ b/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)
diff --git a/python/meep.i b/python/meep.i
index 8ab010cd9..0b1d2fab2 100644
--- a/python/meep.i
+++ b/python/meep.i
@@ -1028,6 +1028,15 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
$1 = (std::complex *)array_data($input);
}
+// typemaps for gaussianbeam
+
+%typecheck(SWIG_TYPECHECK_POINTER, fragment="NumPy_Fragments") std::complex E0[3] {
+ $1 = is_array($input);
+}
+
+%typemap(in) std::complex E0[3] {
+ $1 = (std::complex *)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,
diff --git a/python/simulation.py b/python/simulation.py
index 6c3f80859..41838eb82 100644
--- a/python/simulation.py
+++ b/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,7 +1997,7 @@ 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:
@@ -2005,6 +2005,19 @@ def get_epsilon_point(self, pt, frequency=0, omega=0):
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)
diff --git a/python/source.py b/python/source.py
index ea5c34fc8..ae3cbea10 100644
--- a/python/source.py
+++ b/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
diff --git a/python/tests/gaussianbeam.py b/python/tests/gaussianbeam.py
new file mode 100644
index 000000000..9595f6e03
--- /dev/null
+++ b/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()
diff --git a/src/meep.hpp b/src/meep.hpp
index 74b0464a4..f8fd74efa 100644
--- a/src/meep.hpp
+++ b/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 EO[3]);
+ double eps, double mu, std::complex E0[3]);
void get_fields(std::complex *EH, const vec &x) const;
std::complex get_E0(int n) const { return E0[n]; };