Reformulate Simulation.get_sources() and add 2D GaussianBeam #2333
Conversation
Codecov Report
@@ Coverage Diff @@
## master #2333 +/- ##
==========================================
+ Coverage 71.99% 72.55% +0.55%
==========================================
Files 17 17
Lines 5032 5152 +120
==========================================
+ Hits 3623 3738 +115
- Misses 1409 1414 +5
|
There was a problem hiding this comment.
Thanks for putting this together.
It would be useful to add a unit test for this new feature, perhaps similar to the existing test for GaussianBeamSource in https://github.com/NanoComp/meep/blob/master/python/tests/test_gaussianbeam.py.
|
|
||
| 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.) | ||
| (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.) To use the true solution for a 2d Gaussian Beam, use the `GaussianBeam2DSource` class instead. |
There was a problem hiding this comment.
Given the new GuassianBeam2DSource class, it might be appropriate to rename GaussianBeamSource to GaussianBeam3DSource to clearly distinguish the two?
There was a problem hiding this comment.
This makes sense to me. My only worry is will people struggle to realize the name change when they next update meep and their old scripts fail to recognize GaussianBeamSource?
There was a problem hiding this comment.
they next update meep and their old scripts fail to recognize GaussianBeamSource
I agree. We should probably avoid any breaking changes for something like this... especially since this is preliminary. A future PR will add support for multiple frequencies etc.
There was a problem hiding this comment.
In other words, we shouldn't need to break the API at all. In the future, we can refactor and call GaussianBeamSource in either 2D or 3D and it just "does what it's supposed to do".
Right now that's a challenge due to the complex Hankel functions etc which aren't (readily) compatible with the current 3D codebase in C++.
|
The most backward-compatible option would be something like:
|
|
Performed the following changes:
|
First, I reformulated get_sources. Before, simulation defined logic for each source type separately with
if isinstance(source, GaussianBeam),elif...and then provided the specific logic to add the source parameters to the backend.Now, it simply loops through the list of sources belonging to the simulation and calls add_source on each, passing in self to give the source access to the simulation. Source.py now contains the base logic for adding sources, but
EigenModeSource,IndexedSource,GaussianBeamSource, and (new)GaussianBeam2DSourceoverload this logic.In addition, a new class
GaussianBeamSource2Dis created. This allows a rigorous solution to Maxwell's equations in 2-dimensions via the complex point source method (similar to 3D). All the logic for finding the solution is provided in the python alone, finding the field at a slice. This is done very similar togreen2din near2far.cpp by using Hankel functions of both kinds to generate forward and backward propagating beams when the point source is complex in location. Using this slice, equivalent sources are calculated (also in python).Once the equivalent sources are made (python meep objects),$|E|=1$ at the location of the point source.
GaussianBeam2dSource.add_source()simply callsadd_source()on these equivalent source objects. The fields are normalized to haveThis is all done in python. Because I wrote a green2d python function and one for equivalent source (for single frequency), we can run our code once in python to generate the sources at the beginning of the run and then pass this information in to the backend before the actual meep runs start. Thus, we do not need to write our code in c++ the way 3D was done. Here is a demonstration of it working, under different conditions.
We can change our plane and provide rotations of k, which is done by rotating the imaginary component of the location of the complex point source between x and y.
We can also change our polarization. In order for this to work, an E source and H source are placed in the same location oriented 90 degrees from each other, similar to what's done for 3D.