Python interface for Gaussian beam source

[oskooi]

Initial attempt at a Python interface for the Gaussian beam source based on adding a new class GaussianBeamSource as described in #701:comment. Documentation, test, and tutorial will eventually be added to this PR once the API is working properly.

It compiles but a simple test in 2d aborts with the error:

Using MPI version 3.1, 1 processes
Traceback (most recent call last):
  File "gaussianbeam.py", line 27, in <module>
    beam_E0=beam_E0)]
  File "/home/install/meep/python/meep/source.py", line 545, in __init__
    self.beam_x0 = beam_x0
AttributeError: can't set attribute

The error is occuring during the call to the GaussianBeamSource constructor. It's unclear what is causing this error since the beam_x0 parameter is just a Vector3 object.

import meep as mp
import math
import matplotlib.pyplot as plt

s = 14
resolution = 50
dpml = 2

cell_size = mp.Vector3(s,s)

boundary_layers = [mp.PML(thickness=dpml)]

beam_x0 = mp.Vector3()  # beam center
rot_angle = math.radians(0)  # CW rotation angle about z axis (0: +y axis)
beam_kdir = mp.Vector3(0,1).rotate(mp.Vector3(0,0,1),rot_angle)  # beam propagation direction
beam_w0 = 0.8  # beam waist radius
beam_E0 = mp.Vector3(1,0,0) # beam polarization vector
fcen = 1
sources = [mp.GaussianBeamSource(src=mp.ContinuousSource(fcen),
                                 center=mp.Vector3(),
                                 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)

sim.run(until=20)

sim.plot2D(fields=mp.Hz,
           output_plane=mp.Volume(center=mp.Vector3(),
                                  size=mp.Vector3(s-2*dpml,s-2*dpml)))
plt.savefig('Hz.png')

[oskooi]

Getting closer: the problem now is creating a C++ gaussianbeam object to pass to void fields::add_volume_source(const src_time &src, const volume &, gaussianbeam beam) from Python.

Currently, the line involving gaussianbeam() in python/simulation.py is producing a type error:

        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()]

Using MPI version 3.1, 1 processes
-----------
Initializing structure...
time for choose_chunkdivision = 0.000320795 s
Working in 2D dimensions.
Computational cell is 14 x 14 x 0 with resolution 50
time for set_epsilon = 0.509921 s
-----------
Traceback (most recent call last):
  File "gaussianbeam.py", line 39, in <module>
    sim.run(mp.at_end(mp.output_efield_z),until=20)
  File "/home/install/meep/python/meep/simulation.py", line 3501, in run
    self.init_sim()
  File "/home/install/meep/python/meep/simulation.py", line 1864, in init_sim
    self.add_source(s)
  File "/home/install/meep/python/meep/simulation.py", line 2261, in add_source
    add_vol_src_args = [src.src.swigobj, where, gaussianbeam()]
  File "/home/install/meep/python/meep/__init__.py", line 3650, in __init__
    this = _meep.new_gaussianbeam(x0, kdir, w0, freq, eps, mu, EO)
TypeError: in method 'new_gaussianbeam', argument 7 of type 'std::complex< double > [3]'

From the error message, the problem is the last argument std::complex<double> EO[3] of the gaussianbeam constructor:

  gaussianbeam(const vec &x0, const vec &kdir, double w0, double freq,
               double eps, double mu, std::complex<double> EO[3]);

Python is unable to convert the list of three complex values ([src.beam_E0.x, src.beam_E0.y, src.beam_E0.z] from above) into std::complex<double> EO[3] using the SWIG typemaps.

Any thoughts?

[oskooi]

Renaming the gaussianbeam constructor argument from EO to E0 so that it is recognized by the SWIG typemap produces a segmentation fault when the constructor is called:

Using MPI version 3.1, 1 processes
-----------
Initializing structure...
time for choose_chunkdivision = 0.00032058 s
Working in 2D dimensions.
Computational cell is 14 x 14 x 0 with resolution 50
time for set_epsilon = 0.507612 s
-----------
*** Process received signal ***
Signal: Segmentation fault (11)
Signal code: Address not mapped (1)
Failing at address: 0x3
[ 0] /lib/x86_64-linux-gnu/libpthread.so.0(+0x128a0)[0x7f3cd0a168a0]
[ 1] /home/install/meep/python/meep/_meep.so(+0xfd300)[0x7f3ccf16b300]
[ 2] python3.5(PyCFunction_Call+0x4f)[0x4ea10f]

I also tried passing numpy.array([src.beam_E0.x, src.beam_E0.y, src.beam_E0.z],dtype=numpy.complex128) instead of the list [src.beam_E0.x, src.beam_E0.y, src.beam_E0.z] but that still produces the same error.

Perhaps the SWIG typemap for E0 is somehow still not set up correctly?

[oskooi]

Passing the E0 argument from Python as a NumPy array and updating the SWIG typemaps fixed the problem.

[smartalecH]

How is the beam waist defined here? i.e. is there a 2 in the gaussian exponential?

Would be nice to grab a cross-section of the fields, and plot an analytical gaussian with the same beam waist on top.

[oskooi]

For definitions of all the parameters as well as the formula for the beam itself, see the wikipedia article. There is also a helpful figure that shows the beam waist w0.

[oskooi]

Documentation, tutorial example (part of the revised FAQ), and test (Python) have been added. Unfortunately, because the doc/generate_py_api.py was not working as intended, the Python_User_Interface.md markdown page had to be manually edited rather than building this markdown page directly from the docstrings in python/source.py via Python_User_Interface.md.in. We will need to correct this in a future PR to be consistent with the rest of the user manual.

The test verifies that the energy at the beam focus (for an oblique source) is more than 99% of the maximum energy over the entire cell.

Here are some additional demonstrations included in the tutorial example.