From 8f4e972c9919778a155c4e9bdad0c7dd86dc865e Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Tue, 28 Jul 2020 16:59:57 -0700 Subject: [PATCH 1/6] complex source point Gaussian beam --- src/meep.hpp | 3 ++ src/sources.cpp | 98 +++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 101 insertions(+) diff --git a/src/meep.hpp b/src/meep.hpp index 146f38cab..d8a17d5aa 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1682,6 +1682,9 @@ class fields { void add_volume_source(component c, const src_time &src, const volume &, std::complex amp = 1.0); void require_component(component c); + void gaussianbeam_3d(std::complex *EH, const vec &x0, const vec &x, + const vec &kdir, double w0, double freq, double eps, + double mu, const vec &E0); // mpb.cpp diff --git a/src/sources.cpp b/src/sources.cpp index be07742b2..a561b82ef 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -436,4 +436,102 @@ void fields::add_volume_source(component c, const src_time &src, const volume &w require_component(c); } +/* Given the field-strength vector E0 for a Gaussian beam centered at x0 with beam waist w0 and propagation direction kdir, + compute the E/H fields EH[6] (6 components) at x for a frequency freq in the homogeneous 3d medium eps and mu. + + 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+. */ + +void fields::gaussianbeam_3d(std::complex *EH, const vec &x0, const vec &x, + const vec &kdir, double w0, double freq, double eps, + double mu, const vec &E0) { + if (x0.dim == Dcyl) abort("wrong dimensionality in gaussianbeam_3d"); + + double n = sqrt(eps * mu); + double k = 2 * pi * freq * n; + double impedance = 1/sqrt(eps); + double zR = k * w0 * w0 / 2; // Rayleigh range + double kzR = k*zR; + + vec xrel = x0 - x; + + 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 zc = complex(zhatdotxrel,-zR); + complex Rsq = rho*rho + zc*zc; + complex R = sqrt(Rsq); + complex kR = k*R, kR2 = kR*kR, kR3 = kR2*kR; + + complex f,g,fmgbRsq; + if (abs(kR) > 1e-4) { + complex coskR, sinkR; + if (fabs(imag(kR)) > 30.0) { + UseRescaledFG = true; + complex ExpI = exp(complex(0,1.0)*real(kR)); + complex ExpPlus = exp(imag(kR)-kzR); // should imag(kR) + kzR ? + complex ExpMinus = exp(-(imag(kR)+kzR)); + coskR = 0.5*(ExpI*ExpMinus + conj(ExpI)*ExpPlus); + sinkR = -0.5*complex(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 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 i2fk = 0.5*complex(0,1)*f*k; + + double rnorm = abs(E0); + vec xhat = E0 / 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 gb_Ex = g + fmgbRsq * xhatdotxrel * xhatdotxrel + i2fk * zc; + complex gb_Ey = fmgbRsq * xhatdotxrel * yhatdotxrel; + complex gb_Ez = fmgbRsq * xhatdotxrel * zc - i2fk * xhatdotxrel; + complex gb_Hx = EH[1]; + complex gb_Hy = g + fmgbRsq * yhatdotxrel * yhatdotxrel + i2fk * zc; + complex gb_Hz = fmgbRsq * yhatdotxrel * zc - i2fk * yhatdotxrel; + + vec E_real = (vec(xhatdotxrel,0,0)*real(gb_Ex) + vec(0,yhatdotxrel,0)*real(gb_Ey) + vec(0,0,zhatdotxrel)*real(gb_Ez))*rnorm; + vec H_real = (vec(xhatdotxrel,0,0)*real(gb_Hx) + vec(0,yhatdotxrel,0)*real(gb_Hy) + vec(0,0,zhatdotxrel)*real(gb_Hz))*rnorm; + + vec E_imag = (vec(xhatdotxrel,0,0)*imag(gb_Ex) + vec(0,yhatdotxrel,0)*imag(gb_Ey) + vec(0,0,zhatdotxrel)*imag(gb_Ez))*rnorm; + vec H_imag = (vec(xhatdotxrel,0,0)*imag(gb_Hx) + vec(0,yhatdotxrel,0)*imag(gb_Hy) + vec(0,0,zhatdotxrel)*imag(gb_Hz))*rnorm; + + double Eorig; + if (UseRescaledFG) + Eorig = 3/(2*kzR*kzR*kzR) * (kzR*(kzR-1) + 0.5*(1.0-exp(-2.0*kzR))); + else + Eorig = 3/(2*kzR*kzR*kzR) * (exp(kzR)*kzR*(kzR-1) + sinh(kzR)); + + E_real = E_real / Eorig; + E_imag = E_imag / Eorig; + EH[0] = complex(E_real.x(),E_imag.x()); + EH[1] = complex(E_real.y(),E_imag.y()); + EH[2] = complex(E_real.z(),E_imag.z()); + H_real = H_real / (Eorig*impedance); + H_imag = H_imag / (Eorig*impedance); + EH[3] = complex(H_real.x(),H_imag.x()); + EH[4] = complex(H_real.y(),H_imag.y()); + EH[5] = complex(H_real.z(),H_imag.z()); +} + } // namespace meep From 1440cbb75cad7fb7f3ab86013eb4cedcffd02f26 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Thu, 30 Jul 2020 12:37:43 -0700 Subject: [PATCH 2/6] add gaussianbeam class and add_volume_source function --- src/meep.hpp | 20 +++++- src/sources.cpp | 174 +++++++++++++++++++++++++++++++++++------------- 2 files changed, 143 insertions(+), 51 deletions(-) diff --git a/src/meep.hpp b/src/meep.hpp index d8a17d5aa..6010b5c38 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1484,6 +1484,22 @@ class chunkloop_field_components { std::vector > 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 EO[3]); + void get_fields(std::complex *EH, const vec &x); + +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 E0[3]; // electric field-strength vector +}; + /***************************************************************/ /* prototype for optional user-supplied function to provide an */ /* initial estimate of the wavevector of mode #mode at */ @@ -1681,10 +1697,8 @@ class fields { std::complex A(const vec &), std::complex amp = 1.0); void add_volume_source(component c, const src_time &src, const volume &, std::complex amp = 1.0); + void add_volume_source(const src_time &src, direction d, const volume &, gaussianbeam beam); void require_component(component c); - void gaussianbeam_3d(std::complex *EH, const vec &x0, const vec &x, - const vec &kdir, double w0, double freq, double eps, - double mu, const vec &E0); // mpb.cpp diff --git a/src/sources.cpp b/src/sources.cpp index a561b82ef..b9ce8028d 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -436,27 +436,77 @@ void fields::add_volume_source(component c, const src_time &src, const volume &w require_component(c); } -/* Given the field-strength vector E0 for a Gaussian beam centered at x0 with beam waist w0 and propagation direction kdir, - compute the E/H fields EH[6] (6 components) at x for a frequency freq in the homogeneous 3d medium eps and mu. +static gaussianbeam *global_gaussianbeam_object = 0; +static component global_gaussianbeam_component; + +static complex gaussianbeam_ampfunc(const vec &p) { + std::complex 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, direction d, const volume &where, gaussianbeam beam) { + global_gaussianbeam_object = &beam; + 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(cE[np2], src, where, gaussianbeam_ampfunc, +1.0); + global_gaussianbeam_component = cH[np2]; + add_volume_source(cE[np1], src, where, gaussianbeam_ampfunc, -1.0); + // Nx = +Ey, Ny = -Ex (for d==Z) + global_gaussianbeam_component = cE[np1]; + add_volume_source(cH[np2], src, where, gaussianbeam_ampfunc, -1.0); + global_gaussianbeam_component = cE[np2]; + add_volume_source(cH[np1], src, where, gaussianbeam_ampfunc, +1.0); +} + +gaussianbeam::gaussianbeam(const vec &x0_, const vec &kdir_, double w0_, double freq_, + double eps_, double mu_, std::complex 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+. */ -void fields::gaussianbeam_3d(std::complex *EH, const vec &x0, const vec &x, - const vec &kdir, double w0, double freq, double eps, - double mu, const vec &E0) { - if (x0.dim == Dcyl) abort("wrong dimensionality in gaussianbeam_3d"); - +void gaussianbeam::get_fields(std::complex *EH, const vec &x) { double n = sqrt(eps * mu); double k = 2 * pi * freq * n; - double impedance = 1/sqrt(eps); - double zR = k * w0 * w0 / 2; // Rayleigh range - double kzR = k*zR; - - vec xrel = x0 - x; + double impedance = sqrt(mu/eps); + double z0 = k * w0 * w0 / 2; // Rayleigh range + double kz0 = k*z0; + 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())); @@ -464,19 +514,18 @@ void fields::gaussianbeam_3d(std::complex *EH, const vec &x0, const vec bool UseRescaledFG = false; - complex zc = complex(zhatdotxrel,-zR); + complex zc = complex(zhatdotxrel,-z0); complex Rsq = rho*rho + zc*zc; complex R = sqrt(Rsq); complex kR = k*R, kR2 = kR*kR, kR3 = kR2*kR; - complex f,g,fmgbRsq; if (abs(kR) > 1e-4) { complex coskR, sinkR; if (fabs(imag(kR)) > 30.0) { UseRescaledFG = true; complex ExpI = exp(complex(0,1.0)*real(kR)); - complex ExpPlus = exp(imag(kR)-kzR); // should imag(kR) + kzR ? - complex ExpMinus = exp(-(imag(kR)+kzR)); + complex ExpPlus = exp(imag(kR)-kz0); + complex ExpMinus = exp(-(imag(kR)+kz0)); coskR = 0.5*(ExpI*ExpMinus + conj(ExpI)*ExpPlus); sinkR = -0.5*complex(0,1.0)*(ExpI*ExpMinus - conj(ExpI)*ExpPlus); } else { @@ -495,43 +544,72 @@ void fields::gaussianbeam_3d(std::complex *EH, const vec &x0, const vec } complex i2fk = 0.5*complex(0,1)*f*k; - double rnorm = abs(E0); - vec xhat = E0 / 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 gb_Ex = g + fmgbRsq * xhatdotxrel * xhatdotxrel + i2fk * zc; - complex gb_Ey = fmgbRsq * xhatdotxrel * yhatdotxrel; - complex gb_Ez = fmgbRsq * xhatdotxrel * zc - i2fk * xhatdotxrel; - complex gb_Hx = EH[1]; - complex gb_Hy = g + fmgbRsq * yhatdotxrel * yhatdotxrel + i2fk * zc; - complex gb_Hz = fmgbRsq * yhatdotxrel * zc - i2fk * yhatdotxrel; + complex E[3], H[3]; + for (int j=0; j<3; ++j) { + E[j] = complex(0,0); + H[j] = complex(0,0); + } - vec E_real = (vec(xhatdotxrel,0,0)*real(gb_Ex) + vec(0,yhatdotxrel,0)*real(gb_Ey) + vec(0,0,zhatdotxrel)*real(gb_Ez))*rnorm; - vec H_real = (vec(xhatdotxrel,0,0)*real(gb_Hx) + vec(0,yhatdotxrel,0)*real(gb_Hy) + vec(0,0,zhatdotxrel)*real(gb_Hz))*rnorm; + 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 gb_Ex = g + fmgbRsq * xhatdotxrel * xhatdotxrel + i2fk * zc; + complex gb_Ey = fmgbRsq * xhatdotxrel * yhatdotxrel; + complex gb_Ez = fmgbRsq * xhatdotxrel * zc - i2fk * xhatdotxrel; + complex gb_Hx = EH[1]; + complex gb_Hy = g + fmgbRsq * yhatdotxrel * yhatdotxrel + i2fk * zc; + complex 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()); + } - vec E_imag = (vec(xhatdotxrel,0,0)*imag(gb_Ex) + vec(0,yhatdotxrel,0)*imag(gb_Ey) + vec(0,0,zhatdotxrel)*imag(gb_Ez))*rnorm; - vec H_imag = (vec(xhatdotxrel,0,0)*imag(gb_Hx) + vec(0,yhatdotxrel,0)*imag(gb_Hy) + vec(0,0,zhatdotxrel)*imag(gb_Hz))*rnorm; + 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 gb_Ex = g + fmgbRsq * xhatdotxrel * xhatdotxrel + i2fk * zc; + complex gb_Ey = fmgbRsq * xhatdotxrel * yhatdotxrel; + complex gb_Ez = fmgbRsq * xhatdotxrel * zc - i2fk * xhatdotxrel; + complex gb_Hx = EH[1]; + complex gb_Hy = g + fmgbRsq * yhatdotxrel * yhatdotxrel + i2fk * zc; + complex gb_Hz = fmgbRsq * yhatdotxrel * zc - i2fk * yhatdotxrel; + + E[0] += complex(0,1.0) * inorm * (gb_Ex * xhat.x() + gb_Ey * yhat.x() + gb_Ez * zhat.x()); + E[1] += complex(0,1.0) * inorm * (gb_Ex * xhat.y() + gb_Ey * yhat.y() + gb_Ez * zhat.y()); + E[2] += complex(0,1.0) * inorm * (gb_Ex * xhat.z() + gb_Ey * yhat.z() + gb_Ez * zhat.z()); + H[0] += complex(0,1.0) * inorm * (gb_Hx * xhat.x() + gb_Hy * yhat.x() + gb_Hz * zhat.x()); + H[1] += complex(0,1.0) * inorm * (gb_Hx * xhat.y() + gb_Hy * yhat.y() + gb_Hz * zhat.y()); + H[2] += complex(0,1.0) * inorm * (gb_Hx * xhat.z() + gb_Hy * yhat.z() + gb_Hz * zhat.z()); + } double Eorig; if (UseRescaledFG) - Eorig = 3/(2*kzR*kzR*kzR) * (kzR*(kzR-1) + 0.5*(1.0-exp(-2.0*kzR))); + Eorig = 3/(2*kz0*kz0*kz0) * (kz0*(kz0-1) + 0.5*(1.0-exp(-2.0*kz0))); else - Eorig = 3/(2*kzR*kzR*kzR) * (exp(kzR)*kzR*(kzR-1) + sinh(kzR)); - - E_real = E_real / Eorig; - E_imag = E_imag / Eorig; - EH[0] = complex(E_real.x(),E_imag.x()); - EH[1] = complex(E_real.y(),E_imag.y()); - EH[2] = complex(E_real.z(),E_imag.z()); - H_real = H_real / (Eorig*impedance); - H_imag = H_imag / (Eorig*impedance); - EH[3] = complex(H_real.x(),H_imag.x()); - EH[4] = complex(H_real.y(),H_imag.y()); - EH[5] = complex(H_real.z(),H_imag.z()); + Eorig = 3/(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 From 6ca7c90c402e4811daa5eac41ddd51361193abb6 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Fri, 31 Jul 2020 12:43:01 -0700 Subject: [PATCH 3/6] automatically set direction of source volume to its normal vector and remove direction argument --- src/meep.hpp | 2 +- src/sources.cpp | 9 +++++---- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/meep.hpp b/src/meep.hpp index 6010b5c38..6b4c9b7aa 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1697,7 +1697,7 @@ class fields { std::complex A(const vec &), std::complex amp = 1.0); void add_volume_source(component c, const src_time &src, const volume &, std::complex amp = 1.0); - void add_volume_source(const src_time &src, direction d, const volume &, gaussianbeam beam); + void add_volume_source(const src_time &src, const volume &, gaussianbeam beam); void require_component(component c); // mpb.cpp diff --git a/src/sources.cpp b/src/sources.cpp index b9ce8028d..ee5ef74f6 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -439,7 +439,7 @@ void fields::add_volume_source(component c, const src_time &src, const volume &w static gaussianbeam *global_gaussianbeam_object = 0; static component global_gaussianbeam_component; -static complex gaussianbeam_ampfunc(const vec &p) { +static std::complex gaussianbeam_ampfunc(const vec &p) { std::complex EH[6]; global_gaussianbeam_object->get_fields(EH, p); switch (global_gaussianbeam_component) { @@ -453,8 +453,9 @@ static complex gaussianbeam_ampfunc(const vec &p) { } } -void fields::add_volume_source(const src_time &src, direction d, const volume &where, gaussianbeam beam) { +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) { @@ -600,9 +601,9 @@ void gaussianbeam::get_fields(std::complex *EH, const vec &x) { double Eorig; if (UseRescaledFG) - Eorig = 3/(2*kz0*kz0*kz0) * (kz0*(kz0-1) + 0.5*(1.0-exp(-2.0*kz0))); + Eorig = 3.0/(2*kz0*kz0*kz0) * (kz0*(kz0-1) + 0.5*(1.0-exp(-2.0*kz0))); else - Eorig = 3/(2*kz0*kz0*kz0) * (exp(kz0)*kz0*(kz0-1) + sinh(kz0)); + Eorig = 3.0/(2*kz0*kz0*kz0) * (exp(kz0)*kz0*(kz0-1) + sinh(kz0)); EH[0] = E[0] / Eorig; EH[1] = E[1] / Eorig; From 89cb5eaa0d7e75aa1129bd78f3988c378f08f85b Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Sat, 1 Aug 2020 12:55:09 -0700 Subject: [PATCH 4/6] move add_volume_source_check to private method of fields class and add two separate arguments for te/tm --- src/meep.hpp | 5 +++++ src/mpb.cpp | 28 ++++++++++++++++------------ src/sources.cpp | 20 ++++++++++++++++---- 3 files changed, 37 insertions(+), 16 deletions(-) diff --git a/src/meep.hpp b/src/meep.hpp index 6b4c9b7aa..5886635a5 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1490,6 +1490,7 @@ class gaussianbeam { gaussianbeam(const vec &x0, const vec &kdir, double w0, double freq, double eps, double mu, std::complex EO[3]); void get_fields(std::complex *EH, const vec &x); + std::complex get_E0(int n) { return E0[n]; }; private: vec x0; // beam center @@ -2009,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); + // mpb.cpp + void add_volume_source_check(component c, const src_time &src, const volume &where, + std::complex A(const vec &), std::complex amp, + component c0, direction d, int has_tm, int has_te); public: // monitor.cpp diff --git a/src/mpb.cpp b/src/mpb.cpp index 0d7c615ad..9e0596c4d 100644 --- a/src/mpb.cpp +++ b/src/mpb.cpp @@ -732,18 +732,18 @@ vec get_k(void *vedata) { /* 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; +void fields::add_volume_source_check(component c, const src_time &src, const volume &where, + cdouble A(const vec &), cdouble 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 (f->gv.dim == D2) // parity checks + if (gv.dim == D2) // parity checks { - if ((parity & EVEN_Z_PARITY) && is_tm(c)) return; - if ((parity & ODD_Z_PARITY) && !is_tm(c)) return; + if (has_te && is_tm(c)) return; + if (has_tm && !is_tm(c)) return; }; - f->add_volume_source(c, src, where, A, amp); + add_volume_source(c, src, where, A, amp); } /***************************************************************/ @@ -803,14 +803,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); diff --git a/src/sources.cpp b/src/sources.cpp index ee5ef74f6..455ef675a 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -453,6 +453,14 @@ static std::complex gaussianbeam_ampfunc(const vec &p) { } } +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); @@ -470,14 +478,18 @@ void fields::add_volume_source(const src_time &src, const volume &where, gaussia int np2 = (n + 2) % 3; // Kx = -Hy, Ky = Hx (for d==Z) global_gaussianbeam_component = cH[np1]; - add_volume_source(cE[np2], src, where, gaussianbeam_ampfunc, +1.0); + add_volume_source_check(cE[np2], src, where, gaussianbeam_ampfunc, +1.0, cE[np2], d, + gaussianbeam_has_tm(), gaussianbeam_has_te()); global_gaussianbeam_component = cH[np2]; - add_volume_source(cE[np1], src, where, gaussianbeam_ampfunc, -1.0); + 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(cH[np2], src, where, gaussianbeam_ampfunc, -1.0); + 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(cH[np1], src, where, gaussianbeam_ampfunc, +1.0); + 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_, From 37d05c10b325b1d4f0ab0f925ad856b8c615fb34 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Sat, 1 Aug 2020 18:41:20 -0700 Subject: [PATCH 5/6] move fields::add_volume_source_check definition from mpb.cpp to sources.cpp --- src/meep.hpp | 2 +- src/mpb.cpp | 18 ------------------ src/sources.cpp | 18 ++++++++++++++++++ 3 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/meep.hpp b/src/meep.hpp index 5886635a5..6c769b881 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -2010,7 +2010,7 @@ class fields { void step_source(field_type ft, bool including_integrated = false); void update_pols(field_type ft); void calc_sources(double tim); - // mpb.cpp + // sources.cpp void add_volume_source_check(component c, const src_time &src, const volume &where, std::complex A(const vec &), std::complex amp, component c0, direction d, int has_tm, int has_te); diff --git a/src/mpb.cpp b/src/mpb.cpp index 9e0596c4d..37a642490 100644 --- a/src/mpb.cpp +++ b/src/mpb.cpp @@ -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 fields::add_volume_source_check(component c, const src_time &src, const volume &where, - cdouble A(const vec &), cdouble 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); -} - /***************************************************************/ /* call get_eigenmode() to solve for the specified eigenmode, */ /* then call add_volume_source() to add current sources whose */ diff --git a/src/sources.cpp b/src/sources.cpp index 455ef675a..ee23ffc9d 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -436,6 +436,24 @@ 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 A(const vec &), complex 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; From 5fef7ae284ba6953e3c8ed4abe9fb5032e506a2b Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Mon, 3 Aug 2020 17:43:02 -0700 Subject: [PATCH 6/6] replace has_tm() and has_te() functions with local variables and add citation information to get_fields --- src/meep.hpp | 6 ++-- src/sources.cpp | 81 ++++++++++++++++++++++++++++++++++++------------- 2 files changed, 63 insertions(+), 24 deletions(-) diff --git a/src/meep.hpp b/src/meep.hpp index 6c769b881..a9391750f 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1489,8 +1489,8 @@ class gaussianbeam { public: gaussianbeam(const vec &x0, const vec &kdir, double w0, double freq, double eps, double mu, std::complex EO[3]); - void get_fields(std::complex *EH, const vec &x); - std::complex get_E0(int n) { return E0[n]; }; + void get_fields(std::complex *EH, const vec &x) const; + std::complex get_E0(int n) const { return E0[n]; }; private: vec x0; // beam center @@ -1498,7 +1498,7 @@ class gaussianbeam { double w0; // beam waist radius double freq; // beam frequency double eps, mu; // permittivity/permeability of homogeneous medium - std::complex E0[3]; // electric field-strength vector + std::complex E0[3]; // polarization vector }; /***************************************************************/ diff --git a/src/sources.cpp b/src/sources.cpp index ee23ffc9d..583fd4d08 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -471,14 +471,6 @@ static std::complex gaussianbeam_ampfunc(const vec &p) { } } -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); @@ -492,22 +484,24 @@ void fields::add_volume_source(const src_time &src, const volume &where, gaussia 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, - gaussianbeam_has_tm(), gaussianbeam_has_te()); + has_tm, has_te); 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()); + 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, - gaussianbeam_has_tm(), gaussianbeam_has_te()); + has_tm, 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()); + has_tm, has_te); } gaussianbeam::gaussianbeam(const vec &x0_, const vec &kdir_, double w0_, double freq_, @@ -520,20 +514,65 @@ gaussianbeam::gaussianbeam(const vec &x0_, const vec &kdir_, double w0_, double freq = freq_; eps = eps_; mu = mu_; - for (int j=0; j<3; ++j) + 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+. */ + 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 +*/ -void gaussianbeam::get_fields(std::complex *EH, const vec &x) { +/* 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 *EH, const vec &x) const { 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 ZR = sqrt(mu/eps); + double z0 = k * w0 * w0 / 2; double kz0 = k*z0; vec xrel = x - x0; @@ -638,9 +677,9 @@ void gaussianbeam::get_fields(std::complex *EH, const vec &x) { 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); + EH[3] = H[0] / (Eorig*ZR); + EH[4] = H[1] / (Eorig*ZR); + EH[5] = H[2] / (Eorig*ZR); } } // namespace meep