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 5 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);
std::complex<double> get_E0(int n) { return E0[n]; };
oskooi marked this conversation as resolved.
Show resolved Hide resolved

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]; // electric field-strength 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
207 changes: 207 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,211 @@ 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");
}
}

static int gaussianbeam_has_tm() {
return abs(global_gaussianbeam_object->get_E0(2)) > 0;
}

static int gaussianbeam_has_te() {
return abs(global_gaussianbeam_object->get_E0(0)) > 0 || abs(global_gaussianbeam_object->get_E0(1)) > 0;
}

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");
}
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,
gaussianbeam_has_tm(), gaussianbeam_has_te());
oskooi marked this conversation as resolved.
Show resolved Hide resolved
global_gaussianbeam_component = cH[np2];
add_volume_source_check(cE[np1], src, where, gaussianbeam_ampfunc, -1.0, cE[np1], d,
gaussianbeam_has_tm(), gaussianbeam_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,
gaussianbeam_has_tm(), gaussianbeam_has_te());
global_gaussianbeam_component = cE[np2];
add_volume_source_check(cH[np1], src, where, gaussianbeam_ampfunc, +1.0, cH[np1], d,
gaussianbeam_has_tm(), gaussianbeam_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
(file scuff-em/src/libs/libIncField/GaussianBeam.cc), which is GPL v2+. */

oskooi marked this conversation as resolved.
Show resolved Hide resolved
void gaussianbeam::get_fields(std::complex<double> *EH, const vec &x) {
double n = sqrt(eps * mu);
double k = 2 * pi * freq * n;
double impedance = sqrt(mu/eps);
double z0 = k * w0 * w0 / 2; // Rayleigh range
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*impedance);
EH[4] = H[1] / (Eorig*impedance);
EH[5] = H[2] / (Eorig*impedance);
}

} // namespace meep