complex source point Gaussian beam

[oskooi]

Fixes #711.

This PR adds a new function gaussianbeam_3d to the fields class which computes the fields for a Gaussian beam using the complex source point method. The new function is based on @HomerReid's implementation from SCUFF-EM in GaussianBeam::GetFields. The main difference is that the user-specified field-strength vector E0 is real (because it is based on Meep's vec object which only supports real-valued elements) whereas Homer supported the more general case of complex numbers. We could support a complex field-strength vector but this does not seem to be necessary.

The next thing to do is to add the function add_gaussianbeam_source which calls gaussianbeam_3d for every point in a given user-specified region where and uses these fields to compute the equivalent currents.

A couple of issues regarding the API for add_gaussianbeam_source:

• for a 3d simulation involving a 2d source object, should the center of the Gaussian beam be an input parameter (with arbitrary values) or just be set to the geometric center of the user-specified where object? Under what scenario would the beam center not lie in the 2d plane?
• is there a requirement that the source object be positioned within a lossless medium or should the imaginary part of ε be ignored along with a warning message?
• should we support a 3d source object?
• how can quantify how much more accurate generating a Gaussian beam is using this complex source point method versus just an amplitude function? this could form the basis for a test that we can add to the make check suite.

[stevengj]

I would define two things:

1. a gaussianbeam class that encapsulates the parameters (freq, E0, etc) of the beam, and a method to compute the fields E at a point x.

2. a new method fields::add_volume_source(const src_time &src, const volume &v, gaussianbeam beam); that then calls add_volume_source several times on the same volume v with different transverse components c and an amplitude function A that calls the beam to compute the corresponding field component of the beam.

[stevengj]

Note that the only purpose of this is to create an equivalent current source, which must be a surface (a line in 2d or a plane in 3d).

At the C++ level, the user will have to specify kdir, eps, mu to the gaussianbeam constructor. (In the higher-level Python API, we might infer defaults for these from the source volume.)

[stevengj]

Note that because our source volume is finite, the user will need to ensure that the gaussian beam (evaluated in the source volume) is significantly narrower than the source volume (so that the currents/fields decay to a negligible value at the edge of the source).

(We could even print out a warning if this is not the case, e.g. if the currents at the edge of the source volume are > 1% of the peak current.)

[oskooi]

A new gaussianbeam class and fields::add_volume_source(const src_time &src, direction d, const volume &v, gaussianbeam beam) function have been added. While everything compiles, it has not been tested yet. Perhaps that can wait until we add the Python API since it will simplify things a bit (i.e., computing the normal to the source volume object for the direction d argument, etc).

[oskooi]

The following C++ example tries to launch a Gaussian beam for three test cases involving different values for the beam waist radius and propagation direction. As shown in the accompanying figures, the beam profile is not correct although it does appear that the mode is unidirectional.

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <meep.hpp>

using namespace meep;

double one(const vec &) { return 1.0; }

int main(int argc, char **argv) {
  meep::initialize mpi(argc, argv);
  int sxy = 10;
  int resolution = 50;
  int dpml = 1;

  grid_volume gv = voltwo(sxy, sxy, resolution);
  gv.center_origin();
  structure the_structure(gv, one, pml(dpml));

  fields f(&the_structure);

  double fcen = 1.0;
  continuous_src_time src(fcen);
  src.is_integrated = true;

  vec x0(0,-0.5*sxy+dpml);         // beam center                                                                                                              
  double theta = -35 * pi/180;
  vec kdir(sin(theta),cos(theta)); // beam propagation direction                                                                                               
  double w0 = 1.0;                 // beam waist radius                                                                                                        
  std::complex<double> EO[3] = {1.0,1.0,1.0}; // electric field-strength vector                                                                                
  gaussianbeam beam(x0, kdir, w0, fcen, 1.0, 1.0, EO);

  volume vol(vec(-0.5*sxy,-0.5*sxy+dpml), vec(0.5*sxy,-0.5*sxy+dpml));

  f.add_volume_source(src, vol, beam);

  while (f.time() < 50.0)
    f.step();

  f.output_hdf5(Ez, gv.surroundings());

  return 0;
}

1. w0 = 0.5, theta = 0

2. w0 = 1.0, theta = 25 * pi/180

3. w0 = 1.0, theta = -35 * pi/180

[stevengj]

To use add_volume_source_check, you will need to add it as a private method of fields so that you can call it from other methods.

You will need to change the parity argument to two arguments, say int has_tm and int has_te. Then when you call it from mpb.cpp, pass parity & EVEN_Z_PARITY and parity & ODD_Z_PARITY for those two parameters; in contrast, when you call it from the gaussian-beam code, pass E0[2] != 0 for the has_tm parameter and E0[0] != 0 || E0[1] != 0 for the has_te parameter.

[oskooi]

Working now: unidirectional Gaussian beam. Tested for both polarizations in 2d.

[oskooi]

An example of a Gaussian beam source (green dotted line) and focus (black dot) which are at different locations. (The previous examples had the source and focus at the same location.)

[smartalecH]

Have you checked 3D yet?

[oskooi]

Yes, this is working in 3d as well.

(I also tried removing the field scaling that Homer had included in his GaussianBeam::GetFields to account for potential numerical stability issues related to large complex values for the argument to the sin function and found that it did not affect the results. We may therefore want to remove this from our implementation.)

[smartalecH]

Are you planning to add swig/python support in this PR too? Seems pretty straightforward given all the current machinery.

[oskooi]

I was planning to create a separate PR for the Python API which would also include a test, documentation, and tutorial example.

[oskooi]

Results for 3d Gaussian beams at two different angles: 0° and 30° (rotation about x axis, 0° along +z axis).

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <meep.hpp>

using namespace meep;

double one(const vec &) { return 1.0; }

int main(int argc, char **argv) {
  meep::initialize mpi(argc, argv);
  int s = 14;
  int resolution = 30;
  int dpml = 2;

  grid_volume gv = vol3d(s, s, s, resolution);
  gv.center_origin();
  structure the_structure(gv, one, pml(dpml));

  fields f(&the_structure);

  double fcen = 1.0;
  continuous_src_time src(fcen);

  vec x0(0,0,0);             // beam center
  double theta = 30 * pi/180; // rotation angle around +x axis (0: +z axis)
  vec kdir(0,sin(theta),cos(theta)); // beam propagation direction
  double w0 = 1.0;                   // beam waist radius
  std::complex<double> EO[3] = {1.0,0,0}; // polarization vector
  gaussianbeam beam(x0, kdir, w0, fcen, 1.0, 1.0, EO);

  volume vol(vec(-0.5*s,-0.5*s,0), vec(0.5*s,0.5*s,0));

  f.add_volume_source(src, vol, beam);

  while (f.time() < 10.0)
    f.step();

  volume output_vol(vec(-0.5*s+dpml,-0.5*s+dpml,-0.5*s+dpml),
                    vec(0.5*s-dpml,0.5*s-dpml,0.5*s-dpml));
  f.output_hdf5(Ex, output_vol);

  return 0;
}

$ h5topng -vZc bluered -o beam3d-theta30-xz.png -d ex.r -y 151 ex-000010.00.h5
$ h5topng -vZc bluered -o beam3d-theta30-xy.png -d ex.r -z 151 ex-000010.00.h5
$ h5topng -vZc bluered -o beam3d-theta30-xy.png -d ex.r -z 151 ex-000010.00.h5

[stevengj]

Maybe all of these parameters should be computed once in the gaussianbeam constructor, rather than each time gaussianbeam::get_fields is called?

It probably makes little difference in performance, however, so we might as well leave it as-is for simplicity.

[stevengj]

Would be nice to add a 2d C++ test. Maybe creating a focus like your example above (#1303 (comment)) but at a nonzero angle. Since the code seems to be working now, the test can just check to make sure that the intensity at the focal spot doesn't change in the future.

Or maybe the test and documentation can wait until a future PR adds a Python API?

[stevengj]

This is not quite right for 2d, because we are just taking a slice of a circular-profile Gaussian beam.

In principle, we need two width parameters (and axes thereof), so that we can set the z width to infinity for 2d calculations (and more generally this would allow an elliptical focus). The complex-point-source method doesn't support multiple width parameters in 3d. Instead, the right thing to do seems to be to apply the complex-point-source method to a line of dipoles, i.e. a "2d dipole" source. There are two ways to do this:

1. Just adapt the green2d function from near2far.cpp, extending it to allow a complex source point. However, this would require Hankel functions supporting complex arguments, which we don't (yet) link to. (GSL doesn't support this AFAIK.) There are free-software Bessel functions by Amos (in Fortran) that implement this, which are wrapped by libraries like openspecfun, but it would add another dependency.

2. We could brute-force it — since we are only interested in localized Gaussian beams, we could just integrate the source in the z direction, either using a simple trapezoidal rule and truncating when the field decays enough or using a fancier scheme like Gauss–Hermite quadrature (since the beam is Gaussian with a known width).

[stevengj]

In Python, we'll want a GaussianBeam object that takes the same parameters as gaussianbeam. Then you would pass this for the amp_func parameter of the source object. Then, when Python calls the add_volume_source method, it will check if amp_func is a GaussianBeam, and if so it will instantiate a meep::gaussianbeam object with the same parameters and pass it.

E0 can be a mp.Vector3 (which can be complex). When creating the meep::gaussianbeam object, you will need to convert it to a list [E0.x, E0.y, E0.z] to pass to C++, and you will need something like the following typemap:

%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];
}

[stevengj]

At some point in the future, we could include a function to compute the beam power, and maybe normalize the beam to unit power. Homer provided a routine for this. However, that function won't work in 2d (#1308), where we'd either have to re-derive it or do a numerical (2d) quadrature to get the power. It's probably not worth the trouble, since it will be more reliable to do what you normally do in Meep — perform a separate normalization run.

[oskooi]

Then, when Python calls the add_volume_source method, it will check if amp_func is a GaussianBeam, and if so it will instantiate a meep::gaussianbeam object with the same parameters and pass it.

How would this be set up in python/meep.i?

[oskooi]

An alternative to passing a GaussianBeam object as the amp_func parameter of the Source object (which is somewhat awkward) is to specify the source using a new GaussianBeamSource class object (a subclass of Source) as e.g.:

source = [mp.GaussianBeamSource(src=mp.GaussianSource(fcen,df),
                                center=mp.Vector3(),
                                size=mp.Vector3(sx,sy,0),
                                x0=mp.Vector3(bx,by,bz),
                                kdir=mp.Vector3(kx,ky,kz),
                                w0=w0,
                                E0=mp.Vector3(cx,cy,cz))]

The remaining parameters eps, mu, and freq of the GaussianBeam object can be automatically inferred from center and src.

This approach was described in #711:comment and is more consistent with the API for EigenModeSource.

[stevengj]

That's fine.