From 7bcdee80fe99f1cf0e10e37b01f3fec2a038d2ce Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Tue, 4 Aug 2020 18:45:31 -0700 Subject: [PATCH] complex source point Gaussian beam (#1303) * complex source point Gaussian beam * add gaussianbeam class and add_volume_source function * automatically set direction of source volume to its normal vector and remove direction argument * move add_volume_source_check to private method of fields class and add two separate arguments for te/tm * move fields::add_volume_source_check definition from mpb.cpp to sources.cpp * replace has_tm() and has_te() functions with local variables and add citation information to get_fields --- src/meep.hpp | 22 +++++ src/mpb.cpp | 30 ++---- src/sources.cpp | 246 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 276 insertions(+), 22 deletions(-) diff --git a/src/meep.hpp b/src/meep.hpp index a8859d539..74b0464a4 100644 --- a/src/meep.hpp +++ b/src/meep.hpp @@ -1484,6 +1484,23 @@ 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) const; + std::complex 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 E0[3]; // polarization vector +}; + /***************************************************************/ /* prototype for optional user-supplied function to provide an */ /* initial estimate of the wavevector of mode #mode at */ @@ -1681,6 +1698,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, const volume &, gaussianbeam beam); void require_component(component c); // mpb.cpp @@ -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 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 fd00bb59b..d4e1857a4 100644 --- a/src/mpb.cpp +++ b/src/mpb.cpp @@ -730,24 +730,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 */ @@ -805,14 +787,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 be07742b2..583fd4d08 100644 --- a/src/sources.cpp +++ b/src/sources.cpp @@ -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 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; + +static std::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, 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 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 +*/ + +/* 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 ZR = sqrt(mu/eps); + double z0 = k * w0 * w0 / 2; + 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())); + double zhatdotxrel = zhat & xrel; + + bool UseRescaledFG = false; + + 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)-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 { + 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; + + complex E[3], H[3]; + for (int j=0; j<3; ++j) { + E[j] = complex(0,0); + H[j] = complex(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 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()); + } + + 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.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