Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

complex source point Gaussian beam #1303

Merged
merged 6 commits into from
Aug 5, 2020
Merged

complex source point Gaussian beam #1303

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
8f4e972
complex source point Gaussian beam
oskooi Jul 28, 2020
1440cbb
add gaussianbeam class and add_volume_source function
oskooi Jul 30, 2020
6ca7c90
automatically set direction of source volume to its normal vector and…
oskooi Jul 31, 2020
89cb5ea
move add_volume_source_check to private method of fields class and ad…
oskooi Aug 1, 2020
37d05c1
move fields::add_volume_source_check definition from mpb.cpp to sourc…
oskooi Aug 2, 2020
5fef7ae
replace has_tm() and has_te() functions with local variables and add …
oskooi Aug 4, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 22 additions & 0 deletions src/meep.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1484,6 +1484,23 @@ class chunkloop_field_components {
std::vector<std::complex<double> > values; // updated by update_values(idx)
};

class gaussianbeam {

public:
gaussianbeam(const vec &x0, const vec &kdir, double w0, double freq,
double eps, double mu, std::complex<double> EO[3]);
void get_fields(std::complex<double> *EH, const vec &x) const;
std::complex<double> get_E0(int n) const { return E0[n]; };

private:
vec x0; // beam center
vec kdir; // beam propagation direction
double w0; // beam waist radius
double freq; // beam frequency
double eps, mu; // permittivity/permeability of homogeneous medium
std::complex<double> E0[3]; // polarization vector
};

/***************************************************************/
/* prototype for optional user-supplied function to provide an */
/* initial estimate of the wavevector of mode #mode at */
Expand Down Expand Up @@ -1681,6 +1698,7 @@ class fields {
std::complex<double> A(const vec &), std::complex<double> amp = 1.0);
void add_volume_source(component c, const src_time &src, const volume &,
std::complex<double> amp = 1.0);
void add_volume_source(const src_time &src, const volume &, gaussianbeam beam);
void require_component(component c);

// mpb.cpp
Expand Down Expand Up @@ -1992,6 +2010,10 @@ class fields {
void step_source(field_type ft, bool including_integrated = false);
void update_pols(field_type ft);
void calc_sources(double tim);
// sources.cpp
void add_volume_source_check(component c, const src_time &src, const volume &where,
std::complex<double> A(const vec &), std::complex<double> amp,
component c0, direction d, int has_tm, int has_te);

public:
// monitor.cpp
Expand Down
30 changes: 8 additions & 22 deletions src/mpb.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -728,24 +728,6 @@ vec get_k(void *vedata) {
return vec(edata->Gk[0], edata->Gk[1], edata->Gk[2]);
}

/***************************************************************/
/* helper routine for add_eigenmode_source that calls */
/* add_volume_source only if certain conditions are met */
/***************************************************************/
void add_volume_source_check(component c, const src_time &src, const volume &where,
cdouble A(const vec &), cdouble amp, fields *f, component c0,
direction d, int parity) {
if (!f->gv.has_field(c)) return;
if (c0 != Centered && c0 != c) return;
if (component_direction(c) == d) return;
if (f->gv.dim == D2) // parity checks
{
if ((parity & EVEN_Z_PARITY) && is_tm(c)) return;
if ((parity & ODD_Z_PARITY) && !is_tm(c)) return;
};
f->add_volume_source(c, src, where, A, amp);
}

/***************************************************************/
/* call get_eigenmode() to solve for the specified eigenmode, */
/* then call add_volume_source() to add current sources whose */
Expand Down Expand Up @@ -803,14 +785,18 @@ void fields::add_eigenmode_source(component c0, const src_time &src, direction d
int np2 = (n + 2) % 3;
// Kx = -Hy, Ky = Hx (for d==Z)
global_eigenmode_component = cH[np1];
add_volume_source_check(cE[np2], *src_mpb, where, meep_mpb_A, +1.0 * amp, this, c0, d, parity);
add_volume_source_check(cE[np2], *src_mpb, where, meep_mpb_A, +1.0 * amp, c0, d,
parity & ODD_Z_PARITY, parity & EVEN_Z_PARITY);
global_eigenmode_component = cH[np2];
add_volume_source_check(cE[np1], *src_mpb, where, meep_mpb_A, -1.0 * amp, this, c0, d, parity);
add_volume_source_check(cE[np1], *src_mpb, where, meep_mpb_A, -1.0 * amp, c0, d,
parity & ODD_Z_PARITY, parity & EVEN_Z_PARITY);
// Nx = +Ey, Ny = -Ex (for d==Z)
global_eigenmode_component = cE[np1];
add_volume_source_check(cH[np2], *src_mpb, where, meep_mpb_A, -1.0 * amp, this, c0, d, parity);
add_volume_source_check(cH[np2], *src_mpb, where, meep_mpb_A, -1.0 * amp, c0, d,
parity & ODD_Z_PARITY, parity & EVEN_Z_PARITY);
global_eigenmode_component = cE[np2];
add_volume_source_check(cH[np1], *src_mpb, where, meep_mpb_A, +1.0 * amp, this, c0, d, parity);
add_volume_source_check(cH[np1], *src_mpb, where, meep_mpb_A, +1.0 * amp, c0, d,
parity & ODD_Z_PARITY, parity & EVEN_Z_PARITY);

delete src_mpb;
destroy_eigenmode_data((void *)global_eigenmode_data);
Expand Down
246 changes: 246 additions & 0 deletions src/sources.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -436,4 +436,250 @@ void fields::add_volume_source(component c, const src_time &src, const volume &w
require_component(c);
}

/***************************************************************/
/* helper routine for add_eigenmode_source that calls */
/* add_volume_source only if certain conditions are met */
/***************************************************************/
void fields::add_volume_source_check(component c, const src_time &src, const volume &where,
complex<double> A(const vec &), complex<double> amp, component c0,
direction d, int has_tm, int has_te) {
if (!gv.has_field(c)) return;
if (c0 != Centered && c0 != c) return;
if (component_direction(c) == d) return;
if (gv.dim == D2) // parity checks
{
if (has_te && is_tm(c)) return;
if (has_tm && !is_tm(c)) return;
};
add_volume_source(c, src, where, A, amp);
}

static gaussianbeam *global_gaussianbeam_object = 0;
static component global_gaussianbeam_component;

static std::complex<double> gaussianbeam_ampfunc(const vec &p) {
std::complex<double> EH[6];
global_gaussianbeam_object->get_fields(EH, p);
switch (global_gaussianbeam_component) {
case Ex: return EH[0];
case Ey: return EH[1];
case Ez: return EH[2];
case Hx: return EH[3];
case Hy: return EH[4];
case Hz: return EH[5];
default: abort("invalid component in gaussianbeam_ampfunc");
}
}

void fields::add_volume_source(const src_time &src, const volume &where, gaussianbeam beam) {
global_gaussianbeam_object = &beam;
direction d = normal_direction(where);
component cE[3] = {Ex, Ey, Ez}, cH[3] = {Hx, Hy, Hz};
int n = (d == X ? 0 : (d == Y ? 1 : 2));
if (d == NO_DIRECTION) {
n = where.in_direction(X) == 0
? 0
: where.in_direction(Y) == 0 ? 1 : where.in_direction(Z) == 0 ? 2 : -1;
if (n == -1)
abort(
"can't determine source direction for non-empty source volume with NO_DIRECTION source");
}
bool has_tm = abs(beam.get_E0(2)) > 0;
bool has_te = abs(beam.get_E0(0)) > 0 || abs(beam.get_E0(1)) > 0;
int np1 = (n + 1) % 3;
int np2 = (n + 2) % 3;
// Kx = -Hy, Ky = Hx (for d==Z)
global_gaussianbeam_component = cH[np1];
add_volume_source_check(cE[np2], src, where, gaussianbeam_ampfunc, +1.0, cE[np2], d,
has_tm, has_te);
global_gaussianbeam_component = cH[np2];
add_volume_source_check(cE[np1], src, where, gaussianbeam_ampfunc, -1.0, cE[np1], d,
has_tm, has_te);
// Nx = +Ey, Ny = -Ex (for d==Z)
global_gaussianbeam_component = cE[np1];
add_volume_source_check(cH[np2], src, where, gaussianbeam_ampfunc, -1.0, cH[np2], d,
has_tm, has_te);
global_gaussianbeam_component = cE[np2];
add_volume_source_check(cH[np1], src, where, gaussianbeam_ampfunc, +1.0, cH[np1], d,
has_tm, has_te);
}

gaussianbeam::gaussianbeam(const vec &x0_, const vec &kdir_, double w0_, double freq_,
double eps_, double mu_, std::complex<double> E0_[3]) {
if (x0.dim == Dcyl) abort("wrong dimensionality in gaussianbeam");

x0 = x0_;
kdir = kdir_;
w0 = w0_;
freq = freq_;
eps = eps_;
mu = mu_;
for (int j = 0; j < 3; ++j)
E0[j] = E0_[j];
}

/* Compute the E/H fields EH[6] (6 components) at point x in 3d.

Adapted from code by M. T. Homer Reid in his SCUFF-EM package:
https://github.com/HomerReid/scuff-em/blob/master/libs/libIncField/GaussianBeam.cc
*/

oskooi marked this conversation as resolved.
Show resolved Hide resolved
/* Copyright (C) 2005-2011 M. T. Homer Reid
*
* This file is part of SCUFF-EM.
*
* SCUFF-EM is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* SCUFF-EM is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/

/*
* GaussianBeam.cc -- routine for computing the electric and magnetic
* fields of a focused gaussian beam
*
* homer reid -- 4/2011
*
* new approach:
* johannes feist -- 11/2011
*
* note: this calculation follows CJR Sheppard & S Saghafi, J Opt Soc Am A 16, 1381,
* approximating a Gaussian beam as the field of the sum of an
* x-polarized electric and y-polarized magnetic dipole, but located at
* the _complex_ source point X = (0,0,i z0).
* This produces a field that looks like a Gaussian beam in paraxial approximation
* and exactly solves the field-free Maxwell equations everywhere in (real) space.
*
* note that the core of the calculation is carried out in a coordinate
* system in which the beam propagates in the +z direction and is
* polarized such that the E-field is (mostly) in the +x direction. after
* carrying out the computation described in the memo to obtain the field
* components in this coordinate system, we then rotate the field components
* as appropriate for the user's specified propagation and polarization
* vectors.
*/

void gaussianbeam::get_fields(std::complex<double> *EH, const vec &x) const {
double n = sqrt(eps * mu);
double k = 2 * pi * freq * n;
double ZR = sqrt(mu/eps);
double z0 = k * w0 * w0 / 2;
double kz0 = k*z0;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

vec xrel = x - x0;

vec zhat = kdir / abs(kdir);
double rho = abs(vec(zhat.y()*xrel.z()-zhat.z()*xrel.y(),
zhat.z()*xrel.x()-zhat.x()*xrel.z(),
zhat.x()*xrel.y()-zhat.y()*xrel.x()));
double zhatdotxrel = zhat & xrel;

bool UseRescaledFG = false;

complex<double> zc = complex<double>(zhatdotxrel,-z0);
complex<double> Rsq = rho*rho + zc*zc;
complex<double> R = sqrt(Rsq);
complex<double> kR = k*R, kR2 = kR*kR, kR3 = kR2*kR;
complex<double> f,g,fmgbRsq;
if (abs(kR) > 1e-4) {
complex<double> coskR, sinkR;
if (fabs(imag(kR)) > 30.0) {
UseRescaledFG = true;
complex<double> ExpI = exp(complex<double>(0,1.0)*real(kR));
complex<double> ExpPlus = exp(imag(kR)-kz0);
complex<double> ExpMinus = exp(-(imag(kR)+kz0));
coskR = 0.5*(ExpI*ExpMinus + conj(ExpI)*ExpPlus);
sinkR = -0.5*complex<double>(0,1.0)*(ExpI*ExpMinus - conj(ExpI)*ExpPlus);
} else {
coskR = cos(kR);
sinkR = sin(kR);
}
f = -3.0 * (coskR/kR2 - sinkR/kR3);
g = 1.5 * (sinkR/kR + coskR/kR2 - sinkR/kR3);
fmgbRsq = (f-g)/Rsq;
} else {
complex<double> kR4 = kR2*kR2;
// Taylor series expansion for small R
f = kR4/280.0 - kR2/10.0 + 1.0;
g = 3.0*kR4/280.0 - kR2/5.0 + 1.0;
fmgbRsq = (kR4/5040.0 - kR2/140.0 + 0.1)*(k*k);
}
complex<double> i2fk = 0.5*complex<double>(0,1)*f*k;

complex<double> E[3], H[3];
for (int j=0; j<3; ++j) {
E[j] = complex<double>(0,0);
H[j] = complex<double>(0,0);
}

double rnorm = sqrt(real(E0[0])*real(E0[0])+real(E0[1])*real(E0[1])+real(E0[2])*real(E0[2]));
if (rnorm > 1e-13) {
vec xhat = vec(real(E0[0]),real(E0[1]),real(E0[2])) / rnorm;
vec yhat = vec(zhat.y()*xhat.z()-zhat.z()*xhat.y(),
zhat.z()*xhat.x()-zhat.x()*xhat.z(),
zhat.x()*xhat.y()-zhat.y()*xhat.x());
double xhatdotxrel = xhat & xrel;
double yhatdotxrel = yhat & xrel;

complex<double> gb_Ex = g + fmgbRsq * xhatdotxrel * xhatdotxrel + i2fk * zc;
complex<double> gb_Ey = fmgbRsq * xhatdotxrel * yhatdotxrel;
complex<double> gb_Ez = fmgbRsq * xhatdotxrel * zc - i2fk * xhatdotxrel;
complex<double> gb_Hx = EH[1];
complex<double> gb_Hy = g + fmgbRsq * yhatdotxrel * yhatdotxrel + i2fk * zc;
complex<double> gb_Hz = fmgbRsq * yhatdotxrel * zc - i2fk * yhatdotxrel;

E[0] += rnorm * (gb_Ex * xhat.x() + gb_Ey * yhat.x() + gb_Ez * zhat.x());
E[1] += rnorm * (gb_Ex * xhat.y() + gb_Ey * yhat.y() + gb_Ez * zhat.y());
E[2] += rnorm * (gb_Ex * xhat.z() + gb_Ey * yhat.z() + gb_Ez * zhat.z());
H[0] += rnorm * (gb_Hx * xhat.x() + gb_Hy * yhat.x() + gb_Hz * zhat.x());
H[1] += rnorm * (gb_Hx * xhat.y() + gb_Hy * yhat.y() + gb_Hz * zhat.y());
H[2] += rnorm * (gb_Hx * xhat.z() + gb_Hy * yhat.z() + gb_Hz * zhat.z());
}

double inorm = sqrt(imag(E0[0])*imag(E0[0])+imag(E0[1])*imag(E0[1])+imag(E0[2])*imag(E0[2]));
if (inorm > 1e-13) {
vec xhat = vec(imag(E0[0]),imag(E0[1]),imag(E0[2])) / inorm;
vec yhat = vec(zhat.y()*xhat.z()-zhat.z()*xhat.y(),
zhat.z()*xhat.x()-zhat.x()*xhat.z(),
zhat.x()*xhat.y()-zhat.y()*xhat.x());
double xhatdotxrel = xhat & xrel;
double yhatdotxrel = yhat & xrel;

complex<double> gb_Ex = g + fmgbRsq * xhatdotxrel * xhatdotxrel + i2fk * zc;
complex<double> gb_Ey = fmgbRsq * xhatdotxrel * yhatdotxrel;
complex<double> gb_Ez = fmgbRsq * xhatdotxrel * zc - i2fk * xhatdotxrel;
complex<double> gb_Hx = EH[1];
complex<double> gb_Hy = g + fmgbRsq * yhatdotxrel * yhatdotxrel + i2fk * zc;
complex<double> gb_Hz = fmgbRsq * yhatdotxrel * zc - i2fk * yhatdotxrel;

E[0] += complex<double>(0,1.0) * inorm * (gb_Ex * xhat.x() + gb_Ey * yhat.x() + gb_Ez * zhat.x());
E[1] += complex<double>(0,1.0) * inorm * (gb_Ex * xhat.y() + gb_Ey * yhat.y() + gb_Ez * zhat.y());
E[2] += complex<double>(0,1.0) * inorm * (gb_Ex * xhat.z() + gb_Ey * yhat.z() + gb_Ez * zhat.z());
H[0] += complex<double>(0,1.0) * inorm * (gb_Hx * xhat.x() + gb_Hy * yhat.x() + gb_Hz * zhat.x());
H[1] += complex<double>(0,1.0) * inorm * (gb_Hx * xhat.y() + gb_Hy * yhat.y() + gb_Hz * zhat.y());
H[2] += complex<double>(0,1.0) * inorm * (gb_Hx * xhat.z() + gb_Hy * yhat.z() + gb_Hz * zhat.z());
}

double Eorig;
if (UseRescaledFG)
Eorig = 3.0/(2*kz0*kz0*kz0) * (kz0*(kz0-1) + 0.5*(1.0-exp(-2.0*kz0)));
else
Eorig = 3.0/(2*kz0*kz0*kz0) * (exp(kz0)*kz0*(kz0-1) + sinh(kz0));

EH[0] = E[0] / Eorig;
EH[1] = E[1] / Eorig;
EH[2] = E[2] / Eorig;
EH[3] = H[0] / (Eorig*ZR);
EH[4] = H[1] / (Eorig*ZR);
EH[5] = H[2] / (Eorig*ZR);
}

} // namespace meep