PLaSK library
Loading...
Searching...
No Matches
expansion3d.cpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
14#include <boost/algorithm/clamp.hpp>
15using boost::algorithm::clamp;
16
17#include "../meshadapter.hpp"
18#include "expansion3d.hpp"
19#include "solver3d.hpp"
20
21#define SOLVER static_cast<FourierSolver3D*>(solver)
22
23namespace plask { namespace optical { namespace modal {
24
25constexpr const char* GradientFunctions::NAME;
26constexpr const char* GradientFunctions::UNIT;
27
28ExpansionPW3D::ExpansionPW3D(FourierSolver3D* solver)
29 : Expansion(solver), initialized(false), symmetry_long(E_UNSPECIFIED), symmetry_tran(E_UNSPECIFIED) {}
30
31void ExpansionPW3D::init() {
32 auto geometry = SOLVER->getGeometry();
33
34 RegularAxis long_mesh, tran_mesh;
35
36 periodic_long = geometry->isPeriodic(Geometry3D::DIRECTION_LONG);
37 periodic_tran = geometry->isPeriodic(Geometry3D::DIRECTION_TRAN);
38
39 back = geometry->getChild()->getBoundingBox().lower[0];
40 front = geometry->getChild()->getBoundingBox().upper[0];
41 left = geometry->getChild()->getBoundingBox().lower[1];
42 right = geometry->getChild()->getBoundingBox().upper[1];
43
44 size_t refl = SOLVER->refine_long, reft = SOLVER->refine_tran, nMl, nMt;
45 if (refl == 0) refl = 1;
46 if (reft == 0) reft = 1;
47
48 if (symmetry_long != E_UNSPECIFIED && !geometry->isSymmetric(Geometry3D::DIRECTION_LONG))
49 throw BadInput(solver->getId(), "longitudinal symmetry not allowed for asymmetric structure");
50 if (symmetry_tran != E_UNSPECIFIED && !geometry->isSymmetric(Geometry3D::DIRECTION_TRAN))
51 throw BadInput(solver->getId(), "transverse symmetry not allowed for asymmetric structure");
52
53 if (geometry->isSymmetric(Geometry3D::DIRECTION_LONG)) {
54 if (front <= 0.) {
55 back = -back;
56 front = -front;
57 std::swap(back, front);
58 }
59 if (back != 0.)
60 throw BadInput(SOLVER->getId(), "longitudinally symmetric geometry must have one of its sides at symmetry axis");
61 if (!symmetric_long()) back = -front;
62 }
63 if (geometry->isSymmetric(Geometry3D::DIRECTION_TRAN)) {
64 if (right <= 0.) {
65 left = -left;
66 right = -right;
67 std::swap(left, right);
68 }
69 if (left != 0.)
70 throw BadInput(SOLVER->getId(), "transversely symmetric geometry must have one of its sides at symmetry axis");
71 if (!symmetric_tran()) left = -right;
72 }
73
74 if (!periodic_long) {
75 if (SOLVER->getLongSize() == 0)
76 throw BadInput(solver->getId(),
77 "Flat structure in longitudinal direction (size_long = 0) allowed only for periodic geometry");
78 // Add PMLs
79 if (!symmetric_long()) back -= SOLVER->pml_long.size + SOLVER->pml_long.dist;
80 front += SOLVER->pml_long.size + SOLVER->pml_long.dist;
81 }
82 if (!periodic_tran) {
83 if (SOLVER->getTranSize() == 0)
84 throw BadInput(solver->getId(),
85 "Flat structure in transverse direction (size_tran = 0) allowed only for periodic geometry");
86 // Add PMLs
87 if (!symmetric_tran()) left -= SOLVER->pml_tran.size + SOLVER->pml_tran.dist;
88 right += SOLVER->pml_tran.size + SOLVER->pml_tran.dist;
89 }
90
91 double Ll, Lt;
92
93 eNl = 2 * SOLVER->getLongSize() + 1;
94 eNt = 2 * SOLVER->getTranSize() + 1;
95
96 if (!symmetric_long()) {
97 Ll = front - back;
98 Nl = 2 * SOLVER->getLongSize() + 1;
99 nNl = 4 * SOLVER->getLongSize() + 1;
100 nMl = refl * nNl;
101 double dx = 0.5 * Ll * double(refl - 1) / double(nMl);
102 long_mesh = RegularAxis(back - dx, front - dx - Ll / double(nMl), nMl);
103 } else {
104 Ll = 2 * front;
105 Nl = SOLVER->getLongSize() + 1;
106 nNl = 2 * SOLVER->getLongSize() + 1;
107 nMl = refl * nNl;
108 if (SOLVER->dct2()) {
109 double dx = 0.25 * Ll / double(nMl);
110 long_mesh = RegularAxis(dx, front - dx, nMl);
111 } else {
112 size_t nNa = 4 * SOLVER->getLongSize() + 1;
113 double dx = 0.5 * Ll * double(refl - 1) / double(refl * nNa);
114 long_mesh = RegularAxis(-dx, front + dx, nMl);
115 }
116 } // N = 3 nN = 5 refine = 5 M = 25
117 if (!symmetric_tran()) { // . . 0 . . . . 1 . . . . 2 . . . . 3 . . . . 4 . .
118 Lt = right - left; // ^ ^ ^ ^ ^
119 Nt = 2 * SOLVER->getTranSize() + 1; // |0 1 2 3 4|5 6 7 8 9|0 1 2 3 4|5 6 7 8 9|0 1 2 3 4|
120 nNt = 4 * SOLVER->getTranSize() + 1; // N = 3 nN = 5 refine = 4 M = 20
121 nMt = reft * nNt; // . . 0 . . . 1 . . . 2 . . . 3 . . . 4 . . . 0
122 double dx = 0.5 * Lt * double(reft - 1) / double(nMt); // ^ ^ ^ ^
123 tran_mesh = RegularAxis(left - dx, right - dx - Lt / double(nMt), nMt); // |0 1 2 3|4 5 6 7|8 9 0 1|2 3 4 5|6 7 8 9|
124 } else {
125 Lt = 2 * right; // N = 3 nN = 5 refine = 4 M = 20
126 Nt = SOLVER->getTranSize() + 1; // # . 0 . # . 1 . # . 2 . # . 3 . # . 4 . # . 4 .
127 nNt = 2 * SOLVER->getTranSize() + 1; // ^ ^ ^ ^
128 nMt = reft * nNt; // |0 1 2 3|4 5 6 7|8 9 0 1|2 3 4 5|6 7 8 9|
129 if (SOLVER->dct2()) {
130 double dx = 0.25 * Lt / double(nMt);
131 tran_mesh = RegularAxis(dx, right - dx, nMt);
132 } else {
133 size_t nNa = 4 * SOLVER->getTranSize() + 1;
134 double dx = 0.5 * Lt * double(reft - 1) / double(reft * nNa);
135 tran_mesh = RegularAxis(-dx, right + dx, nMt);
136 }
137 }
138
139 SOLVER->writelog(LOG_DETAIL, "Creating expansion{3} with {0}x{1} plane-waves (matrix size: {2})", Nl, Nt, matrixSize(),
140 (!symmetric_long() && !symmetric_tran()) ? ""
141 : (symmetric_long() && symmetric_tran()) ? " symmetric in longitudinal and transverse directions"
142 : (!symmetric_long() && symmetric_tran()) ? " symmetric in transverse direction"
143 : " symmetric in longitudinal direction");
144
145 if (symmetric_long())
146 SOLVER->writelog(LOG_DETAIL, "Longitudinal symmetry is {0}", (symmetry_long == E_TRAN) ? "Etran" : "Elong");
147 if (symmetric_tran()) SOLVER->writelog(LOG_DETAIL, "Transverse symmetry is {0}", (symmetry_tran == E_TRAN) ? "Etran" : "Elong");
148
149 auto dct_symmetry = SOLVER->dct2() ? FFT::SYMMETRY_EVEN_2 : FFT::SYMMETRY_EVEN_1;
150 auto dct_asymmetry = SOLVER->dct2() ? FFT::SYMMETRY_ODD_2 : FFT::SYMMETRY_ODD_1;
151
152 matFFT = FFT::Forward2D(6, nNl, nNt, symmetric_long() ? dct_symmetry : FFT::SYMMETRY_NONE,
153 symmetric_tran() ? dct_symmetry : FFT::SYMMETRY_NONE);
154
155 if (SOLVER->expansion_rule == FourierSolver3D::RULE_COMBINED) {
156 cos2FFT = FFT::Forward2D(2, nNl, nNt, symmetric_long() ? dct_symmetry : FFT::SYMMETRY_NONE,
157 symmetric_tran() ? dct_symmetry : FFT::SYMMETRY_NONE);
158 cssnFFT = FFT::Forward2D(2, nNl, nNt, symmetric_long() ? dct_asymmetry : FFT::SYMMETRY_NONE,
159 symmetric_tran() ? dct_asymmetry : FFT::SYMMETRY_NONE);
160 }
161
162 // Compute permeability coefficients
163 if (!periodic_long || !periodic_tran) {
164 SOLVER->writelog(LOG_DETAIL, "Adding side PMLs (total structure dimensions: {0}um x {1}um)", Ll, Lt);
165 }
166 if (periodic_long) {
167 mag_long.reset(nNl, Tensor2<dcomplex>(0.));
168 mag_long[0].c00 = 1.;
169 mag_long[0].c11 = 1.; // constant 1
170 } else {
171 mag_long.reset(nNl, Tensor2<dcomplex>(0.));
172 double pb = back + SOLVER->pml_long.size, pf = front - SOLVER->pml_long.size;
173 if (symmetric_long())
174 pib = 0;
175 else
176 pib = std::lower_bound(long_mesh.begin(), long_mesh.end(), pb) - long_mesh.begin();
177 pif = std::lower_bound(long_mesh.begin(), long_mesh.end(), pf) - long_mesh.begin();
178 for (size_t i = 0; i != nNl; ++i) {
179 for (size_t j = refl * i, end = refl * (i + 1); j != end; ++j) {
180 dcomplex s = 1.;
181 if (j < pib) {
182 double h = (pb - long_mesh[j]) / SOLVER->pml_long.size;
183 s = 1. + (SOLVER->pml_long.factor - 1.) * pow(h, SOLVER->pml_long.order);
184 } else if (j > pif) {
185 double h = (long_mesh[j] - pf) / SOLVER->pml_long.size;
186 s = 1. + (SOLVER->pml_long.factor - 1.) * pow(h, SOLVER->pml_long.order);
187 }
188 mag_long[i] += Tensor2<dcomplex>(s, 1. / s);
189 }
190 mag_long[i] /= double(refl);
191 }
192 // Compute FFT
193 FFT::Forward1D(2, nNl, symmetric_long() ? dct_symmetry : FFT::SYMMETRY_NONE)
194 .execute(reinterpret_cast<dcomplex*>(mag_long.data()));
195 // Smooth coefficients
196 if (SOLVER->smooth) {
197 double bb4 = PI / Ll;
198 bb4 *= bb4; // (2π/L)² / 4
199 for (std::size_t i = 0; i != nNl; ++i) {
200 int k = int(i);
201 if (!symmetric_long() && k > int(nNl / 2)) k -= int(nNl);
202 mag_long[i] *= exp(-SOLVER->smooth * bb4 * k * k);
203 }
204 }
205 }
206 if (periodic_tran) {
207 mag_tran.reset(nNt, Tensor2<dcomplex>(0.));
208 mag_tran[0].c00 = 1.;
209 mag_tran[0].c11 = 1.; // constant 1
210 } else {
211 mag_tran.reset(nNt, Tensor2<dcomplex>(0.));
212 double pl = left + SOLVER->pml_tran.size, pr = right - SOLVER->pml_tran.size;
213 if (symmetric_tran())
214 pil = 0;
215 else
216 pil = std::lower_bound(tran_mesh.begin(), tran_mesh.end(), pl) - tran_mesh.begin();
217 pir = std::lower_bound(tran_mesh.begin(), tran_mesh.end(), pr) - tran_mesh.begin();
218 for (size_t i = 0; i != nNt; ++i) {
219 for (size_t j = reft * i, end = reft * (i + 1); j != end; ++j) {
220 dcomplex s = 1.;
221 if (j < pil) {
222 double h = (pl - tran_mesh[j]) / SOLVER->pml_tran.size;
223 s = 1. + (SOLVER->pml_tran.factor - 1.) * pow(h, SOLVER->pml_tran.order);
224 } else if (j > pir) {
225 double h = (tran_mesh[j] - pr) / SOLVER->pml_tran.size;
226 s = 1. + (SOLVER->pml_tran.factor - 1.) * pow(h, SOLVER->pml_tran.order);
227 }
228 mag_tran[i] += Tensor2<dcomplex>(s, 1. / s);
229 }
230 mag_tran[i] /= double(reft);
231 }
232 // Compute FFT
233 FFT::Forward1D(2, nNt, symmetric_tran() ? dct_symmetry : FFT::SYMMETRY_NONE)
234 .execute(reinterpret_cast<dcomplex*>(mag_tran.data()));
235 // Smooth coefficients
236 if (SOLVER->smooth) {
237 double bb4 = PI / Lt;
238 bb4 *= bb4; // (2π/L)² / 4
239 for (std::size_t i = 0; i != nNt; ++i) {
240 int k = int(i);
241 if (!symmetric_tran() && k > int(nNt / 2)) k -= int(nNt);
242 mag_tran[i] *= exp(-SOLVER->smooth * bb4 * k * k);
243 }
244 }
245 }
246
247 // Allocate memory for expansion coefficients
248 size_t nlayers = solver->lcount;
249 coeffs.resize(nlayers);
250 gradients.resize(nlayers);
251 coeffs_ezz.resize(nlayers);
252 coeffs_dexx.assign(nlayers, cmatrix());
253 coeffs_deyy.assign(nlayers, cmatrix());
254 diagonals.assign(nlayers, false);
255
256 mesh = plask::make_shared<RectangularMesh<3>>(plask::make_shared<RegularAxis>(long_mesh),
257 plask::make_shared<RegularAxis>(tran_mesh), solver->verts,
258 RectangularMesh<3>::ORDER_102);
259
260 initialized = true;
261}
262
263void ExpansionPW3D::reset() {
264 coeffs.clear();
265 coeffs_ezz.clear();
266 coeffs_dexx.clear();
267 coeffs_deyy.clear();
268 gradients.clear();
269 initialized = false;
270 k0 = klong = ktran = lam0 = NAN;
271 mesh.reset();
272 temporary.reset();
273}
274
275void ExpansionPW3D::beforeLayersIntegrals(dcomplex lam, dcomplex glam) { SOLVER->prepareExpansionIntegrals(this, mesh, lam, glam); }
276
277template <typename T1, typename T2>
278inline static Tensor3<decltype(T1() * T2())> commutator(const Tensor3<T1>& A, const Tensor3<T2>& B) {
279 return Tensor3<decltype(T1() * T2())>(A.c00 * B.c00 + A.c01 * B.c01, A.c01 * B.c01 + A.c11 * B.c11, A.c22 * B.c22,
280 0.5 * ((A.c00 + A.c11) * B.c01 + A.c01 * (B.c00 + B.c11)));
281}
282
283inline void add_vertex(int l, int t, ExpansionPW3D::Gradient& val, long double& W, const ExpansionPW3D::Gradient::Vertex& vertex) {
284 if (l == vertex.l && t == vertex.t) return;
285 int dl = l - vertex.l, dt = t - vertex.t;
286 long double w = 1. / (dl * dl + dt * dt);
287 w *= w;
288 val += vertex.val * w;
289 W += w;
290}
291
292inline double cf(const ExpansionPW3D::Coeff& c) { return c.c00.real() + c.c11.real(); }
293
294void ExpansionPW3D::layerIntegrals(size_t layer, double lam, double glam) {
295 auto geometry = SOLVER->getGeometry();
296
297 auto long_mesh = static_pointer_cast<RegularAxis>(mesh->lon()), tran_mesh = static_pointer_cast<RegularAxis>(mesh->tran());
298
299 const double Lt = right - left, Ll = front - back;
300 const size_t refl = (SOLVER->refine_long) ? SOLVER->refine_long : 1, reft = (SOLVER->refine_tran) ? SOLVER->refine_tran : 1;
301 const size_t nN = nNl * nNt, NN = Nl * Nt, nNl1 = nNl - 1, nNt1 = nNt - 1;
302 double normlim = min(Ll / double(nNl), Lt / double(nNt)) * 1e-9;
303
304#if defined(OPENMP_FOUND) // && !defined(NDEBUG)
305 SOLVER->writelog(LOG_DETAIL, "Getting refractive indices for layer {}/{} (sampled at {}x{} points) in thread {}", layer,
306 solver->lcount, refl * nNl, reft * nNt, omp_get_thread_num());
307#else
308 SOLVER->writelog(LOG_DETAIL, "Getting refractive indices for layer {}/{} (sampled at {}x{} points)", layer, solver->lcount,
309 refl * nNl, reft * nNt);
310#endif
311
312 if (isnan(lam)) throw BadInput(SOLVER->getId(), "no wavelength given: specify 'lam' or 'lam0'");
313
314 double matv;
315 for (size_t i = 0; i != solver->stack.size(); ++i) {
316 if (solver->stack[i] == layer) {
317 matv = solver->verts->at(i);
318 break;
319 }
320 }
321
322 if (epsilon_connected && solver->lcomputed[layer]) {
323 SOLVER->writelog(LOG_DEBUG, "Layer {:d} takes some materials parameters from inEpsilon", layer);
324 if (isnan(glam)) glam = lam;
325 }
326 if (gain_connected && solver->lgained[layer]) {
327 SOLVER->writelog(LOG_DEBUG, "Layer {:d} has gain", layer);
328 if (isnan(glam)) glam = lam;
329 }
330
331 DataVector<Coeff>& coeffs = this->coeffs[layer];
332
333 // Make space for the result
334 coeffs.reset(nN);
335 std::fill_n(reinterpret_cast<char*>(coeffs.data()), nN * sizeof(Coeff), 0);
336
337 // Normal cos² ans cos·sin for proper expansion rule
338 if (SOLVER->expansion_rule == FourierSolver3D::RULE_COMBINED) {
339 gradients[layer].reset(nN, Gradient(NAN, 0.));
340 }
341 size_t normnans = nN;
342
343 // Average material parameters
344 DataVector<Tensor3<dcomplex>> cell(refl * reft);
345 double nfact = 1. / double(cell.size());
346
347 double pb = back + SOLVER->pml_long.size, pf = front - SOLVER->pml_long.size;
348 double pl = left + SOLVER->pml_tran.size, pr = right - SOLVER->pml_tran.size;
349
350 bool anisotropic = SOLVER->expansion_rule != FourierSolver3D::RULE_INVERSE && !(symmetric_long() || symmetric_tran()),
351 nondiagonal = false;
352
353 // We store data for computing gradients
354 // TODO: it is possible to progressively store only neighbour cells and thus reduce the memory usage,
355 // but it is too hard for now
356 std::unique_ptr<double[]> vals_lock;
357 std::unique_ptr<TempMatrix> vals_temp_matrix;
358 double* vals;
359 size_t nMl = refl * nNl, nMt = reft * nNt;
360 if (SOLVER->expansion_rule == FourierSolver3D::RULE_COMBINED) {
361 if (nMl * nMt > 2 * matrixSize() * matrixSize()) {
362 vals_lock.reset(new double[nMl * nMt]);
363 vals = vals_lock.get();
364 } else {
365 size_t N = matrixSize();
366 vals_temp_matrix.reset(new TempMatrix(&temporary, N, N));
367 vals = reinterpret_cast<double*>(vals_temp_matrix->data());
368 }
369 }
370
371 for (size_t it = 0; it != nNt; ++it) {
372 size_t tbegin = reft * it;
373 size_t tend = tbegin + reft;
374 double tran0 = 0.5 * (tran_mesh->at(tbegin) + tran_mesh->at(tend - 1));
375
376 size_t cto = nNl * it;
377
378 for (size_t il = 0; il != nNl; ++il) {
379 size_t lbegin = refl * il;
380 size_t lend = lbegin + refl;
381 double long0 = 0.5 * (long_mesh->at(lbegin) + long_mesh->at(lend - 1));
382
383 // Store epsilons for a single cell and compute surface normal (only for old rules)
384 Vec<2> norm(0., 0.);
385 for (size_t t = tbegin, j = 0; t != tend; ++t) {
386 for (size_t l = lbegin; l != lend; ++l, ++j) {
387 std::set<std::string> roles;
388 if (epsilon_connected && solver->lcomputed[layer] || gain_connected && solver->lgained[layer])
389 roles = geometry->getRolesAt(vec(long_mesh->at(l), tran_mesh->at(t), matv));
390 bool computed = epsilon_connected && solver->lcomputed[layer] && roles.find("inEpsilon") != roles.end();
391 if (computed) {
392 double W = 0.;
393 Tensor3<dcomplex> val(0.);
394 for (size_t k = 0, v = mesh->index(l, t, 0); k != mesh->vert()->size(); ++v, ++k) {
395 if (solver->stack[k] == layer) {
396 if (isnan(epsilons[v]))
397 throw BadInput(solver->getId(), "complex permittivity tensor got from inEpsilon is NaN at {}",
398 mesh->at(v));
399 double w = (k == 0 || k == mesh->vert()->size() - 1)
400 ? 1e-6
401 : solver->vbounds->at(k) - solver->vbounds->at(k - 1);
402 val += w * epsilons[v];
403 W += w;
404 }
405 }
406 cell[j] = val / W;
407 } else {
408 double T = 0., W = 0., C = 0.;
409 for (size_t k = 0, v = mesh->index(l, t, 0); k != mesh->vert()->size(); ++v, ++k) {
410 if (solver->stack[k] == layer) {
411 double w = (k == 0 || k == mesh->vert()->size() - 1)
412 ? 1e-6
413 : solver->vbounds->at(k) - solver->vbounds->at(k - 1);
414 T += w * temperature[v];
415 C += w * carriers[v];
416 W += w;
417 }
418 }
419 T /= W;
420 C /= W;
421 {
422 OmpLockGuard lock; // this must be declared before `material` to guard its destruction
423 auto material = geometry->getMaterial(vec(long_mesh->at(l), tran_mesh->at(t), matv));
424 lock = material->lock();
425 cell[j] = material->Eps(lam, T, C);
426 if (isnan(cell[j]))
427 throw BadInput(solver->getId(),
428 "complex permittivity tensor (Eps) for {} is NaN at lam={}nm, T={}K n={}/cm3",
429 material->name(), lam, T, C);
430 }
431 if (gain_connected && solver->lgained[layer]) {
432 if (roles.find("QW") != roles.end() || roles.find("QD") != roles.end() ||
433 roles.find("gain") != roles.end()) {
434 Tensor2<double> g = 0.;
435 W = 0.;
436 for (size_t k = 0, v = mesh->index(l, t, 0); k != mesh->vert()->size(); ++v, ++k) {
437 if (solver->stack[k] == layer) {
438 double w = (k == 0 || k == mesh->vert()->size() - 1)
439 ? 1e-6
440 : solver->vbounds->at(k) - solver->vbounds->at(k - 1);
441 g += w * gain[v];
442 W += w;
443 }
444 }
445 Tensor2<double> ni = glam * g / W * (0.25e-7 / PI);
446 double n00 = sqrt(cell[j].c00).real(), n11 = sqrt(cell[j].c11).real(),
447 n22 = sqrt(cell[j].c22).real();
448 cell[j].c00 = dcomplex(n00 * n00 - ni.c00 * ni.c00, 2 * n00 * ni.c00);
449 cell[j].c11 = dcomplex(n11 * n11 - ni.c00 * ni.c00, 2 * n11 * ni.c00);
450 cell[j].c22 = dcomplex(n22 * n22 - ni.c11 * ni.c11, 2 * n22 * ni.c11);
451 cell[j].c01.imag(0.);
452 cell[j].c02.imag(0.);
453 cell[j].c10.imag(0.);
454 cell[j].c12.imag(0.);
455 cell[j].c20.imag(0.);
456 cell[j].c21.imag(0.);
457 }
458 }
459 }
460
461 // Add PMLs
462 if (!periodic_long) {
463 dcomplex s = 1.;
464 if (l < pib) {
465 double h = (pb - long_mesh->at(l)) / SOLVER->pml_long.size;
466 s = 1. + (SOLVER->pml_long.factor - 1.) * pow(h, SOLVER->pml_long.order);
467 } else if (l > pif) {
468 double h = (long_mesh->at(l) - pf) / SOLVER->pml_long.size;
469 s = 1. + (SOLVER->pml_long.factor - 1.) * pow(h, SOLVER->pml_long.order);
470 }
471 cell[j].c00 *= 1. / s;
472 cell[j].c11 *= s;
473 cell[j].c22 *= s;
474 }
475 if (!periodic_tran) {
476 dcomplex s = 1.;
477 if (t < pil) {
478 double h = (pl - tran_mesh->at(t)) / SOLVER->pml_tran.size;
479 s = 1. + (SOLVER->pml_tran.factor - 1.) * pow(h, SOLVER->pml_tran.order);
480 } else if (t > pir) {
481 double h = (tran_mesh->at(t) - pr) / SOLVER->pml_tran.size;
482 s = 1. + (SOLVER->pml_tran.factor - 1.) * pow(h, SOLVER->pml_tran.order);
483 }
484 cell[j].c00 *= s;
485 cell[j].c11 *= 1. / s;
486 cell[j].c22 *= s;
487 }
488
489 if (SOLVER->expansion_rule == FourierSolver3D::RULE_COMBINED)
490 vals[nMl * t + l] = abs(real(cell[j].c00)) + abs(real(cell[j].c11));
491 else
492 norm += (real(cell[j].c00) + real(cell[j].c11)) * vec(long_mesh->at(l) - long0, tran_mesh->at(t) - tran0);
493 }
494 }
495
496 Tensor3<dcomplex> eps(0.);
497 Coeff& coeff = coeffs[cto + il];
498
499 // Compute avg(eps)
500 for (size_t t = tbegin, j = 0; t != tend; ++t) {
501 for (size_t l = lbegin; l != lend; ++l, ++j) {
502 eps += cell[j];
503 }
504 }
505 eps *= nfact;
506 if (SOLVER->expansion_rule == FourierSolver3D::RULE_OLD && eps.c22 != 0.) eps.c22 = 1. / eps.c22;
507
508 double a = abs(norm);
509 if (a >= normlim && SOLVER->expansion_rule != FourierSolver3D::RULE_COMBINED) {
510 norm /= a;
511
512 // Compute avg(eps**(-1))
513 Tensor3<dcomplex> ieps(0.);
514 for (size_t t = tbegin, j = 0; t != tend; ++t) {
515 for (size_t l = lbegin; l != lend; ++l, ++j) {
516 ieps += cell[j].inv();
517 }
518 }
519 ieps *= nfact;
520 // Average permittivity tensor according to:
521 // [ S. G. Johnson and J. D. Joannopoulos, Opt. Express, vol. 8, pp. 173-190 (2001) ]
522 Tensor3<double> P(norm.c0 * norm.c0, norm.c1 * norm.c1, 0., norm.c0 * norm.c1);
523 Tensor3<double> P1(1. - P.c00, 1. - P.c11, 1., -P.c01);
524 eps = commutator(P, ieps.inv()) + commutator(P1, eps);
525 }
526
527 if (!is_zero(eps.c00 - eps.c11)) anisotropic = true;
528 if (!is_zero(eps.c01)) {
529 if (!is_zero(eps.c01 - conj(eps.c10)))
530 throw BadInput(solver->getId(), "complex permittivity tensor (Eps) must be Hermitian for this solver");
531 nondiagonal = true;
532 }
533
534 coeff = eps;
535 }
536 }
537
538 // Compute surface normals using larger cell sizes
539 if (SOLVER->expansion_rule == FourierSolver3D::RULE_COMBINED) {
540 // First increase tolerance to avoid artifacts
541 double vmax = 0., vmin = std::numeric_limits<double>::max();
542 double* v = vals;
543 for (size_t i = 0, nM = nMl * nMt; i < nM; ++i, ++v) {
544 if (*v > vmax) vmax = *v;
545 if (*v < vmin) vmin = *v;
546 }
547 if (vmax != vmin) normlim *= 0.1 * (vmax - vmin);
548
549 bool nowrapl = !periodic_long || symmetric_long(), nowrapt = !periodic_tran || symmetric_tran();
550
551 for (size_t it = 0; it != nNt; ++it) {
552 std::ptrdiff_t tbegin, tend;
553 if (it == 0 && nowrapt) {
554 tbegin = 0;
555 tend = tbegin + reft;
556 } else if (it == nNt - 1 && nowrapt) {
557 tbegin = reft * it - reft / 2;
558 tend = nMt;
559 } else {
560 tbegin = reft * it - reft / 2;
561 tend = tbegin + 2 * reft;
562 }
563 double tran0 = tran_mesh->first() + 0.5 * (tbegin + tend - 1) * tran_mesh->step();
564
565 for (size_t il = 0; il != nNl; ++il) {
566 std::ptrdiff_t lbegin, lend;
567 if (il == 0 && nowrapl) {
568 lbegin = 0;
569 lend = lbegin + refl + refl / 2;
570 } else if (il == nNl - 1 && nowrapl) {
571 lbegin = refl * il - refl / 2;
572 lend = nMl;
573 } else {
574 lbegin = refl * il - refl / 2;
575 lend = lbegin + 2 * refl;
576 }
577 double long0 = long_mesh->first() + 0.5 * (lbegin + lend - 1) * long_mesh->step();
578
579 Vec<2> norm(0., 0.);
580 for (std::ptrdiff_t jt = tbegin; jt != tend; ++jt) {
581 size_t t = (jt < 0) ? jt + nMt : (jt >= nMt) ? jt - nMt : jt;
582 for (std::ptrdiff_t jl = lbegin; jl != lend; ++jl) {
583 size_t l = (jl < 0) ? jl + nMl : (jl >= nMl) ? jl - nMl : jl;
584 norm += vals[nMl * t + l] * vec(long_mesh->at(l) - long0, tran_mesh->at(t) - tran0);
585 }
586 }
587
588 double a = abs(norm);
589 if (a >= normlim) {
590 gradients[layer][nNl * it + il] = norm / a;
591 --normnans;
592 }
593 }
594 }
595 vals_temp_matrix.reset(); // probably not really necessary
596 }
597
598 // int nn = 0;
599 // for (size_t il = 0; il < nNl; ++il) {
600 // for (size_t it = 0; it < nNl; ++it) {
601 // double val = gradients[layer][nNl * it + il].cs.real();
602 // if (isnan(val)) nn++;
603 // std::cerr << str(val, "{:4.1f} ");
604 // }
605 // std::cerr << "\n";
606 // }
607 // std::cerr << ":" << normnans << "/" << nn << "\n";
608
609 // Check if the layer is uniform
610 if (periodic_tran && periodic_long && !nondiagonal) {
611 diagonals[layer] = true;
612 for (size_t i = 1; i != nN; ++i) {
613 if (coeffs[i].c00 != coeffs[i].c22 || coeffs[i].c11 != coeffs[i].c22 || coeffs[i].differs(coeffs[0])) {
614 diagonals[layer] = false;
615 break;
616 }
617 }
618 } else {
619 diagonals[layer] = false;
620 }
621
622 if (diagonals[layer]) {
623 SOLVER->writelog(LOG_DETAIL, "Layer {0} is uniform", layer);
624 std::fill_n(reinterpret_cast<dcomplex*>(coeffs.data() + 1), 6 * (coeffs.size() - 1), 0.);
625 } else {
626 // Perform FFT
627 if (SOLVER->expansion_rule == FourierSolver3D::RULE_INVERSE) {
628 matFFT.execute(reinterpret_cast<dcomplex*>(coeffs.data()), 1);
629 matFFT.execute(reinterpret_cast<dcomplex*>(coeffs.data()) + 2, 1);
630 if (anisotropic) {
631 matFFT.execute(reinterpret_cast<dcomplex*>(coeffs.data()) + 4, 1);
632 } else {
633 for (size_t i = 0; i != nN; ++i) coeffs[i].ic11 = coeffs[i].ic00;
634 }
635 } else if (SOLVER->expansion_rule == FourierSolver3D::RULE_COMBINED) {
636 if (anisotropic) {
637 matFFT.execute(reinterpret_cast<dcomplex*>(coeffs.data()), 5);
638 } else {
639 matFFT.execute(reinterpret_cast<dcomplex*>(coeffs.data()), 3);
640 for (size_t i = 0; i != nN; ++i) {
641 coeffs[i].c11 = coeffs[i].c00;
642 coeffs[i].ic11 = coeffs[i].ic00;
643 }
644 }
645 } else {
646 matFFT.execute(reinterpret_cast<dcomplex*>(coeffs.data()), 2);
647 if (anisotropic) {
648 matFFT.execute(reinterpret_cast<dcomplex*>(coeffs.data()) + 3, 1);
649 } else {
650 for (size_t i = 0; i != nN; ++i) coeffs[i].c11 = coeffs[i].c00;
651 }
652 }
653 if (nondiagonal) matFFT.execute(reinterpret_cast<dcomplex*>(coeffs.data()) + 5, 1);
654
655 // Smooth coefficients
656 if (SOLVER->smooth) {
657 double bb4l = PI / ((front - back) * (symmetric_long() ? 2 : 1));
658 bb4l *= bb4l; // (2π/Ll)² / 4
659 double bb4t = PI / ((right - left) * (symmetric_tran() ? 2 : 1));
660 bb4t *= bb4t; // (2π/Lt)² / 4
661 for (std::size_t it = 0; it != nNt; ++it) {
662 int kt = int(it);
663 if (!symmetric_tran() && kt > int(nNt / 2)) kt -= int(nNt);
664 for (std::size_t il = 0; il != nNl; ++il) {
665 int kl = int(il);
666 if (!symmetric_long() && kl > int(nNl / 2)) kl -= int(nNl);
667 coeffs[nNl * it + il] *= exp(-SOLVER->smooth * (bb4l * kl * kl + bb4t * kt * kt));
668 }
669 }
670 }
671
672 if (SOLVER->expansion_rule != FourierSolver3D::RULE_OLD) {
673 TempMatrix tempMatrix = getTempMatrix();
674 cmatrix work1(NN, NN, tempMatrix.data());
675 cmatrix work2(NN, NN, tempMatrix.data() + NN * NN);
676
677 int ordl = int(SOLVER->getLongSize()), ordt = int(SOLVER->getTranSize());
678
679 char symx = char(symmetric_long() ? 2 * int(symmetry_long) - 3 : 0),
680 symy = char(symmetric_tran() ? 2 * int(symmetry_tran) - 3 : 0);
681 // +1: Ex+, Ey-, Hx-, Hy+
682 // 0: no symmetry
683 // -1: Ex-, Ey+, Hx+, Hy-
684
685 makeToeplitzMatrix(work1, ordl, ordt, layer, 0, -symx, symy);
686 coeffs_ezz[layer].reset(NN, NN);
687 make_unit_matrix(coeffs_ezz[layer]);
688 invmult(work1, coeffs_ezz[layer]);
689
690 if (SOLVER->expansion_rule == FourierSolver3D::RULE_INVERSE) {
691 makeToeplitzMatrix(work1, ordl, ordt, layer, 2, symx, symy);
692 coeffs_dexx[layer].reset(NN, NN);
693 make_unit_matrix(coeffs_dexx[layer]);
694 invmult(work1, coeffs_dexx[layer]);
695 if (anisotropic) {
696 makeToeplitzMatrix(work1, ordl, ordt, layer, 4, -symx, -symy);
697 coeffs_deyy[layer].reset(NN, NN);
698 make_unit_matrix(coeffs_deyy[layer]);
699 invmult(work1, coeffs_deyy[layer]);
700 } else {
701 coeffs_deyy[layer] = coeffs_dexx[layer];
702 }
703
704 } else if (SOLVER->expansion_rule == FourierSolver3D::RULE_COMBINED) {
705 // Fill gaps in cos² and cos·sin matrices
706 if (nN == normnans)
707 throw ComputationError(SOLVER->getId(), "cannot compute normals - consider changing expansion size");
708 std::vector<Gradient::Vertex> vertices;
709 vertices.reserve(nN - normnans);
710 size_t i = 0;
711 for (int t = 0; t < nNt; ++t) {
712 for (int l = 0; l < nNl; ++l, ++i) {
713 if (!gradients[layer][i].isnan()) vertices.emplace_back(l, t, gradients[layer][i]);
714 }
715 }
716 i = 0;
717 for (int t = 0; t < nNt; ++t) {
718 for (int l = 0; l < nNl; ++l, ++i) {
719 if (gradients[layer][i].isnan()) {
720 Gradient val(0., 0.);
721 long double W = 0.;
722 for (const Gradient::Vertex& v : vertices) {
723 add_vertex(l, t, val, W, v);
724 if (symmetric_long()) {
725 add_vertex(l, t, val, W, Gradient::Vertex(-v.l, v.t, v.val));
726 if (periodic_long) add_vertex(l, t, val, W, Gradient::Vertex(2 * nNl - v.l, v.t, v.val));
727 } else if (periodic_long) {
728 add_vertex(l, t, val, W, Gradient::Vertex(v.l + nNl, v.t, v.val));
729 add_vertex(l, t, val, W, Gradient::Vertex(v.l - nNl, v.t, v.val));
730 }
731 if (symmetric_tran()) {
732 add_vertex(l, t, val, W, Gradient::Vertex(v.l, -v.t, v.val));
733 if (periodic_tran) add_vertex(l, t, val, W, Gradient::Vertex(v.l, 2 * nNt - v.t, v.val));
734 } else if (periodic_tran) {
735 add_vertex(l, t, val, W, Gradient::Vertex(v.l, v.t + nNt, v.val));
736 add_vertex(l, t, val, W, Gradient::Vertex(v.l, v.t - nNt, v.val));
737 }
738 }
739 gradients[layer][i] = val / W;
740 }
741 }
742 }
743 for (auto& grad : gradients[layer]) {
744 double c2 = grad.c2.real();
745 if (c2 > 1.) {
746 grad.c2 = c2 = 1.;
747 } // just for safety
748 double cs = sqrt(c2 - c2 * c2);
749 if (grad.cs.real() >= 0)
750 grad.cs = cs;
751 else
752 grad.cs = -cs;
753 }
754
755 // for (size_t il = 0; il < nNl; ++il) {
756 // for (size_t it = 0; it < nNt; ++it) {
757 // auto val = gradients[layer][nNl * it + il];
758 // std::cerr << str(val.c2, "{:.6f}{:+.6f}j/") << str(val.cs, "{:.6f}{:+.6f}j ");
759 // }
760 // }
761 // std::cerr << "\n\n";
762
763 cos2FFT.execute(reinterpret_cast<dcomplex*>(gradients[layer].data()), 1);
764 cssnFFT.execute(reinterpret_cast<dcomplex*>(gradients[layer].data()) + 1, 1);
765
766 // makeToeplitzMatrix(work1, work2, gradients[layer], ordl, ordt, symx, symy);
767 // for (size_t r = 0; r != work1.rows(); ++r) {
768 // for (size_t c = 0; c != work1.cols(); ++c)
769 // std::cerr << str(work1(r,c), "{:7.4f}{:+7.4f}j ");
770 // std::cerr << "\n";
771 // }
772 // std::cerr << "\n";
773 // for (size_t r = 0; r != work2.rows(); ++r) {
774 // for (size_t c = 0; c != work2.cols(); ++c)
775 // std::cerr << str(work2(r,c), "{:7.4f}{:+7.4f}j ");
776 // std::cerr << "\n";
777 // }
778
779 // Smooth coefficients
780 if (SOLVER->grad_smooth) {
781 double bb4l = PI / ((front - back) * (symmetric_long() ? 2 : 1));
782 bb4l *= bb4l; // (2π/Ll)² / 4
783 double bb4t = PI / ((right - left) * (symmetric_tran() ? 2 : 1));
784 bb4t *= bb4t; // (2π/Lt)² / 4
785 for (std::size_t it = 0; it != nNt; ++it) {
786 int kt = int(it);
787 if (!symmetric_tran() && kt > int(nNt / 2)) kt -= int(nNt);
788 for (std::size_t il = 0; il != nNl; ++il) {
789 int kl = int(il);
790 if (!symmetric_long() && kl > int(nNl / 2)) kl -= int(nNl);
791 gradients[layer][nNl * it + il] *= exp(-SOLVER->grad_smooth * (bb4l * kl * kl + bb4t * kt * kt));
792 }
793 }
794 }
795
796 makeToeplitzMatrix(work1, ordl, ordt, layer, 2, symx, symy);
797 // for (size_t r = 0; r < NN; ++r) {
798 // for (size_t c = 0; c < NN; ++c) {
799 // std::cerr << str(work1(r,c), "{:9.6f}{:+9.6f}j") << " ";
800 // }
801 // std::cerr << "\n";
802 // }
803 // std::cerr << "\n";
804 coeffs_dexx[layer].reset(NN, NN);
805 make_unit_matrix(coeffs_dexx[layer]);
806 invmult(work1, coeffs_dexx[layer]);
807 addToeplitzMatrix(coeffs_dexx[layer], ordl, ordt, layer, 1, symx, symy, -1.);
808
809 if (anisotropic || !(symmetric_long() && symmetric_tran())) {
810 makeToeplitzMatrix(work1, ordl, ordt, layer, 4, -symx, -symy);
811 coeffs_deyy[layer].reset(NN, NN);
812 make_unit_matrix(coeffs_deyy[layer]);
813 invmult(work1, coeffs_deyy[layer]);
814 addToeplitzMatrix(coeffs_deyy[layer], ordl, ordt, layer, 3, -symx, -symy, -1.);
815 } else {
816 coeffs_deyy[layer] = coeffs_dexx[layer];
817 }
818 } else {
819 coeffs_dexx[layer].reset();
820 coeffs_deyy[layer].reset();
821 }
822 }
823 }
824}
825
826LazyData<Tensor3<dcomplex>> ExpansionPW3D::getMaterialEps(size_t lay,
827 const shared_ptr<const typename LevelsAdapter::Level>& level,
828 InterpolationMethod interp) {
829 assert(dynamic_pointer_cast<const MeshD<3>>(level->mesh()));
830 auto dest_mesh = static_pointer_cast<const MeshD<3>>(level->mesh());
831
832 if (interp == INTERPOLATION_DEFAULT || interp == INTERPOLATION_FOURIER) {
833 return LazyData<Tensor3<dcomplex>>(dest_mesh->size(), [this, lay, dest_mesh](size_t i) -> Tensor3<dcomplex> {
834 Tensor3<dcomplex> eps(0.);
835 const int nt = symmetric_tran() ? int(nNt) - 1 : int(nNt / 2), nl = symmetric_long() ? int(nNl) - 1 : int(nNl / 2);
836 double Lt = right - left;
837 if (symmetric_tran()) Lt *= 2;
838 double Ll = front - back;
839 if (symmetric_long()) Ll *= 2;
840 for (int kt = -nt; kt <= nt; ++kt) {
841 size_t t = (kt >= 0) ? kt : (symmetric_tran()) ? -kt : kt + nNt;
842 const double phast = kt * (dest_mesh->at(i).c1 - left) / Lt;
843 for (int kl = -nl; kl <= nl; ++kl) {
844 size_t l = (kl >= 0) ? kl : (symmetric_long()) ? -kl : kl + nNl;
845 const double phasl = kl * (dest_mesh->at(i).c0 - back) / Ll;
846 if (SOLVER->expansion_rule != FourierSolver3D::RULE_INVERSE)
847 eps += Tensor3<dcomplex>(coeffs[lay][nNl * t + l]) * exp(2 * PI * I * (phasl + phast));
848 else
849 eps += coeffs[lay][nNl * t + l].toInverseTensor() * exp(2 * PI * I * (phasl + phast));
850 }
851 }
852 if (SOLVER->expansion_rule == FourierSolver3D::RULE_OLD) {
853 eps.c22 = 1. / eps.c22;
854 } else if (SOLVER->expansion_rule == FourierSolver3D::RULE_INVERSE) {
855 eps.c00 = 1. / eps.c00;
856 eps.c11 = 1. / eps.c11;
857 }
858 return eps;
859 });
860 } else {
861 size_t nl = symmetric_long() ? nNl : nNl + 1, nt = symmetric_tran() ? nNt : nNt + 1;
862 DataVector<Tensor3<dcomplex>> params(nl * nt);
863 for (size_t t = 0; t != nNt; ++t) {
864 size_t op = nl * t, oc = nNl * t;
865 if (SOLVER->expansion_rule != FourierSolver3D::RULE_INVERSE)
866 for (size_t l = 0; l != nNl; ++l) params[op + l] = coeffs[lay][oc + l];
867 else
868 for (size_t l = 0; l != nNl; ++l) params[op + l] = coeffs[lay][oc + l].toInverseTensor();
869 }
870 auto dct_symmetry = SOLVER->dct2() ? FFT::SYMMETRY_EVEN_2 : FFT::SYMMETRY_EVEN_1;
871 FFT::Backward2D(9, nNl, nNt, symmetric_long() ? dct_symmetry : FFT::SYMMETRY_NONE,
872 symmetric_tran() ? dct_symmetry : FFT::SYMMETRY_NONE, nl)
873 .execute(reinterpret_cast<dcomplex*>(params.data()));
874 shared_ptr<RegularAxis> lcmesh = plask::make_shared<RegularAxis>(), tcmesh = plask::make_shared<RegularAxis>();
875 if (symmetric_long()) {
876 if (SOLVER->dct2()) {
877 double dx = 0.5 * (front - back) / double(nl);
878 lcmesh->reset(back + dx, front - dx, nl);
879 } else {
880 lcmesh->reset(0., front, nl);
881 }
882 } else {
883 lcmesh->reset(back, front, nl);
884 for (size_t t = 0, end = nl * nt, dist = nl - 1; t != end; t += nl) params[dist + t] = params[t];
885 }
886 if (symmetric_tran()) {
887 if (SOLVER->dct2()) {
888 double dy = 0.5 * right / double(nt);
889 tcmesh->reset(dy, right - dy, nt);
890 } else {
891 tcmesh->reset(0., right, nt);
892 }
893 } else {
894 tcmesh->reset(left, right, nt);
895 for (size_t l = 0, last = nl * (nt - 1); l != nl; ++l) params[last + l] = params[l];
896 }
897 for (Tensor3<dcomplex>& eps : params) {
898 if (SOLVER->expansion_rule == FourierSolver3D::RULE_OLD) {
899 eps.c22 = 1. / eps.c22;
900 } else if (SOLVER->expansion_rule == FourierSolver3D::RULE_INVERSE) {
901 eps.c00 = 1. / eps.c00;
902 eps.c11 = 1. / eps.c11;
903 }
904 }
905 auto src_mesh = plask::make_shared<RectangularMesh<3>>(
906 lcmesh, tcmesh, plask::make_shared<RegularAxis>(level->vpos(), level->vpos(), 1), RectangularMesh<3>::ORDER_210);
907 return interpolate(
908 src_mesh, params, dest_mesh, interp,
909 InterpolationFlags(SOLVER->getGeometry(),
910 symmetric_long() ? InterpolationFlags::Symmetry::POSITIVE : InterpolationFlags::Symmetry::NO,
911 symmetric_tran() ? InterpolationFlags::Symmetry::POSITIVE : InterpolationFlags::Symmetry::NO,
912 InterpolationFlags::Symmetry::POSITIVE));
913 }
914}
915
916void ExpansionPW3D::getMatrices(size_t lay, cmatrix& RE, cmatrix& RH) {
917 assert(initialized);
918 if (isnan(k0)) throw BadInput(SOLVER->getId(), "wavelength or k0 not set");
919 if (isinf(k0.real())) throw BadInput(SOLVER->getId(), "wavelength must not be 0");
920
921 bool diagonal = diagonals[lay];
922
923 int ordl = int(SOLVER->getLongSize()), ordt = int(SOLVER->getTranSize());
924
925 char symx = char(symmetric_long() ? 2 * int(symmetry_long) - 3 : 0),
926 symy = char(symmetric_tran() ? 2 * int(symmetry_tran) - 3 : 0);
927 // +1: Ex+, Ey-, Hx-, Hy+
928 // 0: no symmetry
929 // -1: Ex-, Ey+, Hx+, Hy-
930
931 assert(!(symx && klong != 0.));
932 assert(!(symy && ktran != 0.));
933
934 assert(!isnan(k0.real()) && !isnan(k0.imag()));
935
936 double Gx = 2. * PI / (front - back) * (symx ? 0.5 : 1.), Gy = 2. * PI / (right - left) * (symy ? 0.5 : 1.);
937
938 dcomplex ik0 = 1. / k0;
939
940 const size_t N = Nl * Nt;
941 assert((symx ? ordl + 1 : 2 * ordl + 1) == Nl);
942 assert((symy ? ordt + 1 : 2 * ordt + 1) == Nt);
943
944 if (SOLVER->expansion_rule == FourierSolver3D::RULE_COMBINED && !coeffs_dexx[lay].empty() && !diagonal) {
945 TempMatrix tempMatrix = getTempMatrix();
946 cmatrix workxx(N, N, tempMatrix.data());
947 cmatrix workxy(N, N, tempMatrix.data() + N * N);
948 cmatrix workyx(N, N, tempMatrix.data() + 2 * N * N);
949 cmatrix workyy(N, N, tempMatrix.data() + 3 * N * N);
950 cmatrix workc2(N, N, RE.data());
951 cmatrix workcs(N, N, RE.data() + N * N);
952
953 makeToeplitzMatrix(workc2, workcs, gradients[lay], ordl, ordt, symx, symy);
954
955 mult_matrix_by_matrix(coeffs_dexx[lay], workc2, workxx);
956 mult_matrix_by_matrix(coeffs_dexx[lay], workcs, workxy);
957 if (!(symmetric_long() || symmetric_tran())) {
958 if (coeffs_deyy[lay].data() != coeffs_dexx[lay].data()) {
959 mult_matrix_by_matrix(coeffs_deyy[lay], workc2, workyy);
960 mult_matrix_by_matrix(coeffs_deyy[lay], workcs, workyx);
961 } else {
962 workyy = workxx;
963 workyx = workxy;
964 }
965 } else {
966 makeToeplitzMatrix(workc2, workcs, gradients[lay], ordl, ordt, -symx, -symy);
967 mult_matrix_by_matrix(coeffs_deyy[lay], workc2, workyy);
968 mult_matrix_by_matrix(coeffs_deyy[lay], workcs, workyx);
969 }
970
971 zero_matrix(RE);
972 zero_matrix(RH);
973
974 for (int iy = (symy ? 0 : -ordt); iy <= ordt; ++iy) {
975 dcomplex gy = iy * Gy - ktran;
976 for (int ix = (symx ? 0 : -ordl); ix <= ordl; ++ix) {
977 dcomplex gx = ix * Gx - klong;
978 size_t iex = iEx(ix, iy), iey = iEy(ix, iy);
979 size_t ihx = iHx(ix, iy), ihy = iHy(ix, iy);
980 size_t i = idx(ix, iy);
981
982 for (int jy = -ordt; jy <= ordt; ++jy) {
983 dcomplex py = jy * Gy - ktran;
984 int ijy = iy - jy;
985 if (symy && ijy < 0) ijy = -ijy;
986 for (int jx = -ordl; jx <= ordl; ++jx) {
987 bool toeplitz = true;
988 dcomplex px = jx * Gx - klong;
989 int ijx = ix - jx;
990 if (symx && ijx < 0) ijx = -ijx;
991 size_t jex = iEx(jx, jy), jey = iEy(jx, jy);
992 size_t jhx = iHx(jx, jy), jhy = iHy(jx, jy);
993 double fx = 1., fy = 1.;
994 if (jx < 0 && symx) {
995 fx *= symx;
996 fy *= -symx;
997 toeplitz = false;
998 }
999 if (jy < 0 && symy) {
1000 fx *= symy;
1001 fy *= -symy;
1002 toeplitz = false;
1003 }
1004 RH(iex, jhy) += fx * k0 * muyy(lay, ijx, ijy);
1005 RH(iey, jhx) += fy * k0 * muxx(lay, ijx, ijy);
1006 dcomplex imu = imuzz(lay, ijx, ijy) * ik0;
1007 RE(ihy, jex) += fx * (-gy * py * imu + k0 * epsxx(lay, ijx, ijy));
1008 RE(ihx, jex) += fx * (gx * py * imu + k0 * epsyx(lay, ijx, ijy));
1009 RE(ihy, jey) += fy * (gy * px * imu + k0 * epsxy(lay, ijx, ijy));
1010 RE(ihx, jey) += fy * (-gx * px * imu + k0 * epsyy(lay, ijx, ijy));
1011 if (toeplitz) {
1012 size_t j = idx(jx, jy);
1013 dcomplex iepszz = coeffs_ezz[lay](i, j) * ik0;
1014 RH(iex, jhy) += fx * (-gx * px * iepszz);
1015 RH(iey, jhy) += fx * (-gy * px * iepszz);
1016 RH(iex, jhx) += fy * (-gx * py * iepszz);
1017 RH(iey, jhx) += fy * (-gy * py * iepszz);
1018 RE(ihy, jex) += fx * (k0 * workxx(i, j));
1019 RE(ihx, jex) += fx * (k0 * conj(workyx(i, j))); // Should this conjugate be here?
1020 RE(ihy, jey) += fy * (k0 * workxy(i, j));
1021 RE(ihx, jey) += fy * (k0 * (coeffs_deyy[lay](i, j) - workyy(i, j))); // ([[1]]-[[c²]]) [[Δε]]
1022 }
1023 }
1024 }
1025
1026 // Ugly hack to avoid singularity
1027 if (RE(iex, iex) == 0.) RE(iex, iex) = 1e-32;
1028 if (RE(iey, iey) == 0.) RE(iey, iey) = 1e-32;
1029 if (RH(ihx, ihx) == 0.) RH(ihx, ihx) = 1e-32;
1030 if (RH(ihy, ihy) == 0.) RH(ihy, ihy) = 1e-32;
1031 }
1032 }
1033
1034 } else if (SOLVER->expansion_rule == FourierSolver3D::RULE_INVERSE && !coeffs_dexx[lay].empty() && !diagonal) {
1035 zero_matrix(RE);
1036 zero_matrix(RH);
1037
1038 for (int iy = (symy ? 0 : -ordt); iy <= ordt; ++iy) {
1039 dcomplex gy = iy * Gy - ktran;
1040 for (int ix = (symx ? 0 : -ordl); ix <= ordl; ++ix) {
1041 dcomplex gx = ix * Gx - klong;
1042 size_t iex = iEx(ix, iy), iey = iEy(ix, iy);
1043 size_t ihx = iHx(ix, iy), ihy = iHy(ix, iy);
1044 size_t i = idx(ix, iy);
1045
1046 for (int jy = -ordt; jy <= ordt; ++jy) {
1047 dcomplex py = jy * Gy - ktran;
1048 int ijy = iy - jy;
1049 if (symy && ijy < 0) ijy = -ijy;
1050 for (int jx = -ordl; jx <= ordl; ++jx) {
1051 bool toeplitz = true;
1052 dcomplex px = jx * Gx - klong;
1053 int ijx = ix - jx;
1054 if (symx && ijx < 0) ijx = -ijx;
1055 size_t jex = iEx(jx, jy), jey = iEy(jx, jy);
1056 size_t jhx = iHx(jx, jy), jhy = iHy(jx, jy);
1057 double fx = 1., fy = 1.;
1058 if (symx && jx < 0) {
1059 fx *= symx;
1060 fy *= -symx;
1061 toeplitz = false;
1062 }
1063 if (symy && jy < 0) {
1064 fx *= symy;
1065 fy *= -symy;
1066 toeplitz = false;
1067 }
1068 RH(iex, jhy) += fx * k0 * muyy(lay, ijx, ijy);
1069 RH(iey, jhx) += fy * k0 * muxx(lay, ijx, ijy);
1070 dcomplex imu = imuzz(lay, ijx, ijy) * ik0;
1071 RE(ihy, jex) += fx * (-gy * py * imu);
1072 RE(ihx, jex) += fx * (gx * py * imu + k0 * epsyx(lay, ijx, ijy));
1073 RE(ihy, jey) += fy * (gy * px * imu + k0 * epsxy(lay, ijx, ijy));
1074 RE(ihx, jey) += fy * (-gx * px * imu);
1075 if (toeplitz) {
1076 size_t j = idx(jx, jy);
1077 dcomplex iepszz = coeffs_ezz[lay](i, j) * ik0;
1078 RH(iex, jhy) += fx * (-gx * px * iepszz);
1079 RH(iey, jhy) += fx * (-gy * px * iepszz);
1080 RH(iex, jhx) += fy * (-gx * py * iepszz);
1081 RH(iey, jhx) += fy * (-gy * py * iepszz);
1082 RE(ihy, jex) += fx * (k0 * coeffs_dexx[lay](i, j));
1083 RE(ihx, jey) += fy * (k0 * coeffs_deyy[lay](i, j));
1084 }
1085 }
1086 }
1087
1088 // Ugly hack to avoid singularity
1089 if (RE(iex, iex) == 0.) RE(iex, iex) = 1e-32;
1090 if (RE(iey, iey) == 0.) RE(iey, iey) = 1e-32;
1091 if (RH(ihx, ihx) == 0.) RH(ihx, ihx) = 1e-32;
1092 if (RH(ihy, ihy) == 0.) RH(ihy, ihy) = 1e-32;
1093 }
1094 }
1095
1096 } else {
1097 zero_matrix(RE);
1098 zero_matrix(RH);
1099
1100 for (int iy = (symy ? 0 : -ordt); iy <= ordt; ++iy) {
1101 dcomplex gy = iy * Gy - ktran;
1102 for (int ix = (symx ? 0 : -ordl); ix <= ordl; ++ix) {
1103 dcomplex gx = ix * Gx - klong;
1104 size_t iex = iEx(ix, iy), iey = iEy(ix, iy);
1105 size_t ihx = iHx(ix, iy), ihy = iHy(ix, iy);
1106
1107 for (int jy = -ordt; jy <= ordt; ++jy) {
1108 dcomplex py = jy * Gy - ktran;
1109 int ijy = iy - jy;
1110 if (symy && ijy < 0) ijy = -ijy;
1111 for (int jx = -ordl; jx <= ordl; ++jx) {
1112 dcomplex px = jx * Gx - klong;
1113 int ijx = ix - jx;
1114 if (symx && ijx < 0) ijx = -ijx;
1115 size_t jex = iEx(jx, jy), jey = iEy(jx, jy);
1116 size_t jhx = iHx(jx, jy), jhy = iHy(jx, jy);
1117 double fx = 1., fy = 1.;
1118 if (symx && jx < 0) {
1119 fx *= symx;
1120 fy *= -symx;
1121 }
1122 if (symy && jy < 0) {
1123 fx *= symy;
1124 fy *= -symy;
1125 }
1126 dcomplex ieps = ((SOLVER->expansion_rule == FourierSolver3D::RULE_OLD) ? iepszz(lay, ijx, ijy)
1127 : diagonal ? ((ix == jx && iy == jy) ? 1. / coeffs[lay][0].c22 : 0.)
1128 : iepszz(lay, ix, jx, iy, jy)) *
1129 ik0;
1130 RH(iex, jhy) += fx * (-gx * px * ieps + k0 * muyy(lay, ijx, ijy));
1131 RH(iey, jhy) += fx * (-gy * px * ieps);
1132 RH(iex, jhx) += fy * (-gx * py * ieps);
1133 RH(iey, jhx) += fy * (-gy * py * ieps + k0 * muxx(lay, ijx, ijy));
1134 dcomplex imu = imuzz(lay, ijx, ijy) * ik0;
1135 RE(ihy, jex) += fx * (-gy * py * imu + k0 * epsxx(lay, ijx, ijy));
1136 RE(ihx, jex) += fx * (gx * py * imu + k0 * epsyx(lay, ijx, ijy));
1137 RE(ihy, jey) += fy * (gy * px * imu + k0 * epsxy(lay, ijx, ijy));
1138 RE(ihx, jey) += fy * (-gx * px * imu + k0 * epsyy(lay, ijx, ijy));
1139 }
1140 }
1141
1142 // Ugly hack to avoid singularity
1143 if (RE(iex, iex) == 0.) RE(iex, iex) = 1e-32;
1144 if (RE(iey, iey) == 0.) RE(iey, iey) = 1e-32;
1145 if (RH(ihx, ihx) == 0.) RH(ihx, ihx) = 1e-32;
1146 if (RH(ihy, ihy) == 0.) RH(ihy, ihy) = 1e-32;
1147 }
1148 }
1149 }
1150
1151 assert(!RE.isnan());
1152 assert(!RH.isnan());
1153}
1154
1155void ExpansionPW3D::prepareField() {
1156 if (field_interpolation == INTERPOLATION_DEFAULT) field_interpolation = INTERPOLATION_FOURIER;
1157 if (symmetric_long() || symmetric_tran()) {
1158 Component syml = (which_field == FIELD_E || !symmetry_long) ? symmetry_long : Component(3 - symmetry_long),
1159 symt = (which_field == FIELD_E || !symmetry_tran) ? symmetry_tran : Component(3 - symmetry_tran);
1160 size_t nl = (syml == E_UNSPECIFIED) ? Nl + 1 : Nl;
1161 size_t nt = (symt == E_UNSPECIFIED) ? Nt + 1 : Nt;
1162 if (field_interpolation != INTERPOLATION_FOURIER) {
1163 int df = SOLVER->dct2() ? 0 : 4;
1164 FFT::Symmetry x1, xz2, yz1, y2;
1165 if (symmetric_long()) {
1166 x1 = FFT::Symmetry(3 - syml + df);
1167 yz1 = FFT::Symmetry(syml + df);
1168 } else {
1169 x1 = yz1 = FFT::SYMMETRY_NONE;
1170 }
1171 if (symmetric_tran()) {
1172 xz2 = FFT::Symmetry(3 - symt + df);
1173 y2 = FFT::Symmetry(symt + df);
1174 } else {
1175 xz2 = y2 = FFT::SYMMETRY_NONE;
1176 }
1177 fft_x = FFT::Backward2D(3, Nl, Nt, x1, xz2, nl);
1178 fft_y = FFT::Backward2D(3, Nl, Nt, yz1, y2, nl);
1179 fft_z = FFT::Backward2D(3, Nl, Nt, yz1, xz2, nl);
1180 }
1181 field.reset(nl * nt);
1182 } else {
1183 if (field_interpolation != INTERPOLATION_FOURIER)
1184 fft_z = FFT::Backward2D(3, Nl, Nt, FFT::SYMMETRY_NONE, FFT::SYMMETRY_NONE, Nl + 1);
1185 field.reset((Nl + 1) * (Nt + 1));
1186 }
1187}
1188
1189void ExpansionPW3D::cleanupField() {
1190 field.reset();
1191 fft_x = FFT::Backward2D();
1192 fft_y = FFT::Backward2D();
1193 fft_z = FFT::Backward2D();
1194}
1195
1196// TODO fields must be carefully verified
1197
1198LazyData<Vec<3, dcomplex>> ExpansionPW3D::getField(size_t l,
1199 const shared_ptr<const typename LevelsAdapter::Level>& level,
1200 const cvector& E,
1201 const cvector& H) {
1202 Component syml = (which_field == FIELD_E) ? symmetry_long : Component((3 - symmetry_long) % 3);
1203 Component symt = (which_field == FIELD_E) ? symmetry_tran : Component((3 - symmetry_tran) % 3);
1204
1205 bool diagonal = diagonals[l];
1206
1207 size_t nl = (syml == E_UNSPECIFIED) ? Nl + 1 : Nl, nt = (symt == E_UNSPECIFIED) ? Nt + 1 : Nt;
1208
1209 const dcomplex kx = klong, ky = ktran;
1210
1211 int ordl = int(SOLVER->getLongSize()), ordt = int(SOLVER->getTranSize());
1212
1213 double bl = 2 * PI / (front - back) * (symmetric_long() ? 0.5 : 1.0),
1214 bt = 2 * PI / (right - left) * (symmetric_tran() ? 0.5 : 1.0);
1215
1216 assert(dynamic_pointer_cast<const MeshD<3>>(level->mesh()));
1217 auto dest_mesh = static_pointer_cast<const MeshD<3>>(level->mesh());
1218 double vpos = level->vpos();
1219
1220 int dxl = 0, dyl = 0, dxt = 0, dyt = 0;
1221 if (field_interpolation != INTERPOLATION_FOURIER) {
1222 if (symmetric_long()) {
1223 if (syml == E_TRAN)
1224 dxl = 1;
1225 else
1226 dyl = 1;
1227 for (size_t t = 0, end = nl * nt; t != end; t += nl) field[nl - 1 + t] = Vec<3, dcomplex>(0., 0., 0.);
1228 }
1229 if (symmetric_tran()) {
1230 if (symt == E_TRAN)
1231 dxt = 1;
1232 else
1233 dyt = 1;
1234 for (size_t l = 0, off = nl * (nt - 1); l != Nl; ++l) field[off + l] = Vec<3, dcomplex>(0., 0., 0.);
1235 }
1236 }
1237
1238 if (which_field == FIELD_E) {
1239 for (int it = symmetric_tran() ? 0 : -ordt; it <= ordt; ++it) {
1240 for (int il = symmetric_long() ? 0 : -ordl; il <= ordl; ++il) {
1241 // How expensive is checking conditions in each loop?
1242 // Fuck it! The code is much more clear this way.
1243 size_t iex = nl * (((it < 0) ? Nt + it : it) - dxt) + ((il < 0) ? Nl + il : il) - dxl;
1244 size_t iey = nl * (((it < 0) ? Nt + it : it) - dyt) + ((il < 0) ? Nl + il : il) - dyl;
1245 size_t iez = nl * (((it < 0) ? Nt + it : it) - dxt) + ((il < 0) ? Nl + il : il) - dyl;
1246 if (!(it == 0 && dxt) && !(il == 0 && dxl)) field[iex].lon() = E[iEx(il, it)];
1247 if (!(it == 0 && dyt) && !(il == 0 && dyl)) field[iey].tran() = E[iEy(il, it)];
1248 if (!(it == 0 && dxt) && !(il == 0 && dyl)) {
1249 field[iez].vert() = 0.;
1250 for (int jt = -ordt; jt <= ordt; ++jt)
1251 for (int jl = -ordl; jl <= ordl; ++jl) {
1252 double fhx =
1253 ((jl < 0 && symmetry_long == E_LONG) ? -1. : 1) * ((jt < 0 && symmetry_tran == E_LONG) ? -1. : 1);
1254 double fhy =
1255 ((jl < 0 && symmetry_long == E_TRAN) ? -1. : 1) * ((jt < 0 && symmetry_tran == E_TRAN) ? -1. : 1);
1256 field[iez].vert() +=
1257 ((SOLVER->expansion_rule == FourierSolver3D::RULE_OLD) ? iepszz(l, il - jl, it - jt)
1258 : diagonal ? ((il == jl && it == jt) ? 1. / coeffs[l][0].c22 : 0.)
1259 : iepszz(l, il, jl, it, jt)) *
1260 ((bl * double(jl) - kx) * fhy * H[iHy(jl, jt)] + (bt * double(jt) - ky) * fhx * H[iHx(jl, jt)]);
1261 }
1262 field[iez].vert() /= k0;
1263 }
1264 }
1265 }
1266 } else { // field == FIELD_H
1267 for (int it = symmetric_tran() ? 0 : -ordt; it <= ordt; ++it) {
1268 for (int il = symmetric_long() ? 0 : -ordl; il <= ordl; ++il) {
1269 size_t ihx = nl * (((it < 0) ? Nt + it : it) - dxt) + ((il < 0) ? Nl + il : il) - dxl;
1270 size_t ihy = nl * (((it < 0) ? Nt + it : it) - dyt) + ((il < 0) ? Nl + il : il) - dyl;
1271 size_t ihz = nl * (((it < 0) ? Nt + it : it) - dxt) + ((il < 0) ? Nl + il : il) - dyl;
1272 if (!(it == 0 && dxt) && !(il == 0 && dxl)) field[ihx].lon() = -H[iHx(il, it)];
1273 if (!(it == 0 && dyt) && !(il == 0 && dyl)) field[ihy].tran() = H[iHy(il, it)];
1274 if (!(it == 0 && dxt) && !(il == 0 && dyl)) {
1275 field[ihz].vert() = 0.;
1276 for (int jt = -ordt; jt <= ordt; ++jt)
1277 for (int jl = -ordl; jl <= ordl; ++jl) {
1278 double fex =
1279 ((jl < 0 && symmetry_long == E_TRAN) ? -1. : 1) * ((jt < 0 && symmetry_tran == E_TRAN) ? -1. : 1);
1280 double fey =
1281 ((jl < 0 && symmetry_long == E_LONG) ? -1. : 1) * ((jt < 0 && symmetry_tran == E_LONG) ? -1. : 1);
1282 field[ihz].vert() += imuzz(l, il - jl, it - jt) * (-(bl * double(jl) - kx) * fey * E[iEy(jl, jt)] +
1283 (bt * double(jt) - ky) * fex * E[iEx(jl, jt)]);
1284 }
1285 field[ihz].vert() /= k0;
1286 }
1287 }
1288 }
1289 }
1290
1291 if (field_interpolation == INTERPOLATION_FOURIER) {
1292 const double lo0 = symmetric_long() ? -front : back, hi0 = front, lo1 = symmetric_tran() ? -right : left, hi1 = right;
1293 DataVector<Vec<3, dcomplex>> result(dest_mesh->size());
1294 double Ll = (symmetric_long() ? 2. : 1.) * (front - back), Lt = (symmetric_tran() ? 2. : 1.) * (right - left);
1295 dcomplex bl = 2. * PI * I / Ll, bt = 2. * PI * I / Lt;
1296 dcomplex ikx = I * kx, iky = I * ky;
1297 result.reset(dest_mesh->size(), Vec<3, dcomplex>(0., 0., 0.));
1298 for (int it = -ordt; it <= ordt; ++it) {
1299 double ftx = 1., fty = 1.;
1300 size_t iit;
1301 if (it < 0) {
1302 if (symmetric_tran()) {
1303 if (symt == E_LONG)
1304 fty = -1.;
1305 else
1306 ftx = -1.;
1307 iit = nl * (-it);
1308 } else {
1309 iit = nl * (Nt + it);
1310 }
1311 } else {
1312 iit = nl * it;
1313 }
1314 dcomplex gt = bt * double(it) - iky;
1315 for (int il = -ordl; il <= ordl; ++il) {
1316 double flx = 1., fly = 1.;
1317 size_t iil;
1318 if (il < 0) {
1319 if (symmetric_long()) {
1320 if (syml == E_LONG)
1321 fly = -1.;
1322 else
1323 flx = -1.;
1324 iil = -il;
1325 } else {
1326 iil = Nl + il;
1327 }
1328 } else {
1329 iil = il;
1330 }
1331 Vec<3, dcomplex> coeff = field[iit + iil];
1332 coeff.c0 *= ftx * flx;
1333 coeff.c1 *= fty * fly;
1334 coeff.c2 *= ftx * fly;
1335 dcomplex gl = bl * double(il) - ikx;
1336 for (size_t ip = 0; ip != dest_mesh->size(); ++ip) {
1337 auto p = dest_mesh->at(ip);
1338 if (!periodic_long) p.c0 = clamp(p.c0, lo0, hi0);
1339 if (!periodic_tran) p.c1 = clamp(p.c1, lo1, hi1);
1340 result[ip] += coeff * exp(gl * (p.c0 - back) + gt * (p.c1 - left));
1341 }
1342 }
1343 }
1344 return result;
1345 } else {
1346 if (symmetric_long() || symmetric_tran()) {
1347 fft_x.execute(&(field.data()->lon()), 1);
1348 fft_y.execute(&(field.data()->tran()), 1);
1349 fft_z.execute(&(field.data()->vert()), 1);
1350 double dx, dy;
1351 if (symmetric_tran()) {
1352 dy = 0.5 * (right - left) / double(nt);
1353 } else {
1354 for (size_t l = 0, off = nl * Nt; l != Nl; ++l) field[off + l] = field[l];
1355 dy = 0.;
1356 }
1357 if (symmetric_long()) {
1358 dx = 0.5 * (front - back) / double(nl);
1359 } else {
1360 for (size_t t = 0, end = nl * nt; t != end; t += nl) field[Nl + t] = field[t];
1361 dx = 0.;
1362 }
1363 auto src_mesh = plask::make_shared<RectangularMesh<3>>(plask::make_shared<RegularAxis>(back + dx, front - dx, nl),
1364 plask::make_shared<RegularAxis>(left + dy, right - dy, nt),
1365 plask::make_shared<RegularAxis>(vpos, vpos, 1),
1366 RectangularMesh<3>::ORDER_210);
1367 LazyData<Vec<3, dcomplex>> interpolated = interpolate(
1368 src_mesh, field, dest_mesh, field_interpolation,
1369 InterpolationFlags(SOLVER->getGeometry(),
1370 symmetric_long() ? InterpolationFlags::Symmetry::POSITIVE : InterpolationFlags::Symmetry::NO,
1371 symmetric_tran() ? InterpolationFlags::Symmetry::POSITIVE : InterpolationFlags::Symmetry::NO,
1372 InterpolationFlags::Symmetry::NO),
1373 false);
1374
1375 return LazyData<Vec<3, dcomplex>>(interpolated.size(),
1376 [interpolated, dest_mesh, syml, symt, kx, ky, this](size_t i) -> Vec<3, dcomplex> {
1377 Vec<3, dcomplex> result = interpolated[i];
1378 if (symmetric_long()) {
1379 double Ll = 2. * front;
1380 if (syml == E_TRAN) {
1381 double x = std::fmod(dest_mesh->at(i)[0], Ll);
1382 if ((-front <= x && x < 0) || x > front) {
1383 result.lon() = -result.lon();
1384 result.vert() = -result.vert();
1385 }
1386 } else {
1387 double x = std::fmod(dest_mesh->at(i)[0], Ll);
1388 if ((-front <= x && x < 0) || x > front) {
1389 result.tran() = -result.tran();
1390 }
1391 }
1392 } else {
1393 dcomplex ikx = I * kx;
1394 result[i] *= exp(-ikx * dest_mesh->at(i).c0);
1395 }
1396 if (symmetric_tran()) {
1397 double Lt = 2. * right;
1398 if (symt == E_TRAN) {
1399 double y = std::fmod(dest_mesh->at(i)[1], Lt);
1400 if ((-right <= y && y < 0) || y > right) {
1401 result.lon() = -result.lon();
1402 result.vert() = -result.vert();
1403 }
1404 } else {
1405 double y = std::fmod(dest_mesh->at(i)[1], Lt);
1406 if ((-right <= y && y < 0) || y > right) {
1407 result.tran() = -result.tran();
1408 }
1409 }
1410 } else {
1411 dcomplex iky = I * ky;
1412 result *= exp(-iky * dest_mesh->at(i).c1);
1413 }
1414 return result;
1415 });
1416 } else {
1417 fft_z.execute(reinterpret_cast<dcomplex*>(field.data()));
1418 for (size_t l = 0, off = nl * Nt; l != Nl; ++l) field[off + l] = field[l];
1419 for (size_t t = 0, end = nl * nt; t != end; t += nl) field[Nl + t] = field[t];
1420 auto src_mesh = plask::make_shared<RectangularMesh<3>>(
1421 plask::make_shared<RegularAxis>(back, front, nl), plask::make_shared<RegularAxis>(left, right, nt),
1422 plask::make_shared<RegularAxis>(vpos, vpos, 1), RectangularMesh<3>::ORDER_210);
1423 LazyData<Vec<3, dcomplex>> interpolated =
1424 interpolate(src_mesh, field, dest_mesh, field_interpolation,
1425 InterpolationFlags(SOLVER->getGeometry(), InterpolationFlags::Symmetry::NO,
1426 InterpolationFlags::Symmetry::NO, InterpolationFlags::Symmetry::NO),
1427 false);
1428 dcomplex ikx = I * kx, iky = I * ky;
1429 return LazyData<Vec<3, dcomplex>>(interpolated.size(), [interpolated, dest_mesh, ikx, iky](size_t i) {
1430 return interpolated[i] * exp(-ikx * dest_mesh->at(i).c0 - iky * dest_mesh->at(i).c1);
1431 });
1432 }
1433 }
1434}
1435
1436double ExpansionPW3D::integratePoyntingVert(const cvector& E, const cvector& H) {
1437 double P = 0.;
1438
1439 int ordl = int(SOLVER->getLongSize()), ordt = int(SOLVER->getTranSize());
1440
1441 for (int iy = -ordt; iy <= ordt; ++iy) {
1442 for (int ix = -ordl; ix <= ordl; ++ix) {
1443 P += real(E[iEx(ix, iy)] * conj(H[iHy(ix, iy)]) + E[iEy(ix, iy)] * conj(H[iHx(ix, iy)]));
1444 }
1445 }
1446
1447 double dlong = symmetric_long() ? 2 * front : front - back, dtran = symmetric_tran() ? 2 * right : right - left;
1448 return P * dlong * dtran * 1e-12; // µm² -> m²
1449}
1450
1451void ExpansionPW3D::getDiagonalEigenvectors(cmatrix& Te, cmatrix& Te1, const cmatrix& RE, const cdiagonal&) {
1452 size_t nr = Te.rows(), nc = Te.cols();
1453 std::fill_n(Te.data(), nr * nc, 0.);
1454 std::fill_n(Te1.data(), nr * nc, 0.);
1455
1456 // Ensure that for the same gamma E*H [2x2] is diagonal
1457 assert(nc % 2 == 0);
1458 size_t n = nc / 2;
1459 for (std::size_t i = 0; i < n; i++) {
1460 // Compute Te1 = sqrt(RE)
1461 // https://en.wikipedia.org/wiki/Square_root_of_a_2_by_2_matrix
1462 // but after this normalize columns to 1
1463 dcomplex a = RE(2 * i, 2 * i), b = RE(2 * i, 2 * i + 1), c = RE(2 * i + 1, 2 * i), d = RE(2 * i + 1, 2 * i + 1);
1464 if (is_zero(b) && is_zero(c)) {
1465 Te1(2 * i, 2 * i) = Te1(2 * i + 1, 2 * i + 1) = Te(2 * i, 2 * i) = Te(2 * i + 1, 2 * i + 1) = 1.;
1466 Te1(2 * i, 2 * i + 1) = Te1(2 * i + 1, 2 * i) = Te(2 * i, 2 * i + 1) = Te(2 * i + 1, 2 * i) = 0.;
1467 } else {
1468 dcomplex s = sqrt(a * d - b * c);
1469 a += s;
1470 d += s;
1471 // Normalize
1472 s = 1. / sqrt(a * a + b * b);
1473 a *= s;
1474 b *= s;
1475 s = 1. / sqrt(c * c + d * d);
1476 c *= s;
1477 d *= s;
1478 Te1(2 * i, 2 * i) = a;
1479 Te1(2 * i, 2 * i + 1) = b;
1480 Te1(2 * i + 1, 2 * i) = c;
1481 Te1(2 * i + 1, 2 * i + 1) = d;
1482 // Invert Te1
1483 s = 1. / (a * d - b * c);
1484 Te(2 * i, 2 * i) = s * d;
1485 Te(2 * i, 2 * i + 1) = -s * b;
1486 Te(2 * i + 1, 2 * i) = -s * c;
1487 Te(2 * i + 1, 2 * i + 1) = s * a;
1488 }
1489 }
1490}
1491
1492double ExpansionPW3D::integrateField(WhichField field,
1493 size_t layer,
1494 const cmatrix& TE,
1495 const cmatrix& TH,
1496 const std::function<std::pair<dcomplex, dcomplex>(size_t, size_t)>& vertical) {
1497 Component syml = (field == FIELD_E) ? symmetry_long : Component((3 - symmetry_long) % 3);
1498 Component symt = (field == FIELD_E) ? symmetry_tran : Component((3 - symmetry_tran) % 3);
1499
1500 bool diagonal = diagonals[layer];
1501
1502 size_t nl = (syml == E_UNSPECIFIED) ? Nl + 1 : Nl, nt = (symt == E_UNSPECIFIED) ? Nt + 1 : Nt;
1503 const dcomplex kx = klong, ky = ktran;
1504
1505 int ordl = int(SOLVER->getLongSize()), ordt = int(SOLVER->getTranSize());
1506
1507 double bl = 2 * PI / (front - back) * (symmetric_long() ? 0.5 : 1.0),
1508 bt = 2 * PI / (right - left) * (symmetric_tran() ? 0.5 : 1.0);
1509
1510 assert(TE.rows() == matrixSize());
1511 assert(TH.rows() == matrixSize());
1512
1513 size_t M = TE.cols();
1514 assert(TH.cols() == M);
1515
1516 size_t NN = Nl * Nt;
1517 TempMatrix temp = getTempMatrix();
1518 cmatrix Fz(NN, M, temp.data());
1519
1520 if (field == FIELD_E) {
1521 PLASK_OMP_PARALLEL_FOR
1522 for (openmp_size_t m = 0; m < M; m++) {
1523 for (int it = symmetric_tran() ? 0 : -ordt; it <= ordt; ++it) {
1524 for (int il = symmetric_long() ? 0 : -ordl; il <= ordl; ++il) {
1525 dcomplex vert = 0.;
1526 for (int jt = -ordt; jt <= ordt; ++jt) {
1527 for (int jl = -ordl; jl <= ordl; ++jl) {
1528 double fhx =
1529 ((jl < 0 && symmetry_long == E_LONG) ? -1. : 1) * ((jt < 0 && symmetry_tran == E_LONG) ? -1. : 1);
1530 double fhy =
1531 ((jl < 0 && symmetry_long == E_TRAN) ? -1. : 1) * ((jt < 0 && symmetry_tran == E_TRAN) ? -1. : 1);
1532 vert += ((SOLVER->expansion_rule == FourierSolver3D::RULE_OLD) ? iepszz(layer, il - jl, it - jt)
1533 : diagonal ? ((il == jl && it == jt) ? 1. / coeffs[layer][0].c22 : 0.)
1534 : iepszz(layer, il, jl, it, jt)) *
1535 ((bl * double(jl) - kx) * fhy * TH(iHy(jl, jt), m) +
1536 (bt * double(jt) - ky) * fhx * TH(iHx(jl, jt), m));
1537 }
1538 }
1539 Fz(idx(il, it), m) = vert / k0;
1540 }
1541 }
1542 }
1543 } else { // field == FIELD_H
1544 PLASK_OMP_PARALLEL_FOR
1545 for (openmp_size_t m = 0; m < M; m++) {
1546 for (int it = symmetric_tran() ? 0 : -ordt; it <= ordt; ++it) {
1547 for (int il = symmetric_long() ? 0 : -ordl; il <= ordl; ++il) {
1548 dcomplex vert = 0.;
1549 for (int jt = -ordt; jt <= ordt; ++jt) {
1550 for (int jl = -ordl; jl <= ordl; ++jl) {
1551 double fex =
1552 ((jl < 0 && symmetry_long == E_TRAN) ? -1. : 1) * ((jt < 0 && symmetry_tran == E_TRAN) ? -1. : 1);
1553 double fey =
1554 ((jl < 0 && symmetry_long == E_LONG) ? -1. : 1) * ((jt < 0 && symmetry_tran == E_LONG) ? -1. : 1);
1555 vert += imuzz(layer, il - jl, it - jt) * (-(bl * double(jl) - kx) * fey * TE(iEy(jl, jt), m) +
1556 (bt * double(jt) - ky) * fex * TE(iEx(jl, jt), m));
1557 }
1558 }
1559 Fz(idx(il, it), m) = vert / k0;
1560 }
1561 }
1562 }
1563 }
1564
1565 double result = 0.;
1566
1567 if (field == FIELD_E) {
1568 PLASK_OMP_PARALLEL_FOR
1569 for (openmp_size_t m1 = 0; m1 < M; ++m1) {
1570 for (openmp_size_t m2 = m1; m2 < M; ++m2) {
1571 dcomplex resxy = 0., resz = 0.;
1572 for (int t = -ordt; t <= ordt; ++t) {
1573 for (int l = -ordl; l <= ordl; ++l) {
1574 size_t ix = iEx(l, t), iy = iEy(l, t), i = idx(l, t);
1575 resxy += TE(ix, m1) * conj(TE(ix, m2)) + TE(iy, m1) * conj(TE(iy, m2));
1576 resz += Fz(i, m1) * conj(Fz(i, m2));
1577 }
1578 }
1579 if (!(is_zero(resxy) && is_zero(resz))) {
1580 auto vert = vertical(m1, m2);
1581 double res = real(resxy * vert.first + resz * vert.second);
1582 if (m2 != m1) res *= 2;
1583#pragma omp atomic
1584 result += res;
1585 }
1586 }
1587 }
1588
1589 } else { // field == FIELD_H
1590 PLASK_OMP_PARALLEL_FOR
1591 for (openmp_size_t m1 = 0; m1 < M; ++m1) {
1592 for (openmp_size_t m2 = m1; m2 < M; ++m2) {
1593 dcomplex resxy = 0., resz = 0.;
1594 for (int t = -ordt; t <= ordt; ++t) {
1595 for (int l = -ordl; l <= ordl; ++l) {
1596 size_t ix = iHx(l, t), iy = iHy(l, t), i = idx(l, t);
1597 resxy += TH(ix, m1) * conj(TH(ix, m2)) + TH(iy, m1) * conj(TH(iy, m2));
1598 resz += Fz(i, m1) * conj(Fz(i, m2));
1599 }
1600 }
1601 if (!(is_zero(resxy) && is_zero(resz))) {
1602 auto vert = vertical(m1, m2);
1603 double res = real(resxy * vert.second + resz * vert.first);
1604 if (m2 != m1) res *= 2;
1605#pragma omp atomic
1606 result += res;
1607 }
1608 }
1609 }
1610 }
1611
1612 const double area = (front - back) * (symmetric_long() ? 2. : 1.) * (right - left) * (symmetric_tran() ? 2. : 1.);
1613 return 0.5 * result * area;
1614}
1615
1616LazyData<double> ExpansionPW3D::getGradients(GradientFunctions::EnumType what,
1617 const shared_ptr<const typename LevelsAdapter::Level>& level,
1618 InterpolationMethod interp) {
1619 double z = level->vpos();
1620 const size_t lay = SOLVER->stack[solver->getLayerFor(z)];
1621 const int which = int(what);
1622
1623 assert(dynamic_pointer_cast<const MeshD<3>>(level->mesh()));
1624 auto dest_mesh = static_pointer_cast<const MeshD<3>>(level->mesh());
1625
1626 if (diagonals[lay] || SOLVER->expansion_rule != FourierSolver3D::RULE_COMBINED)
1627 return LazyData<double>(dest_mesh->size(), [](size_t i) { return 0.0; });
1628
1629 if (interp == INTERPOLATION_DEFAULT || interp == INTERPOLATION_FOURIER) {
1630 return LazyData<double>(dest_mesh->size(), [this, lay, which, dest_mesh](size_t i) -> double {
1631 double res = 0.;
1632 const int nt = symmetric_tran() ? int(nNt) - 1 : int(nNt / 2), nl = symmetric_long() ? int(nNl) - 1 : int(nNl / 2);
1633 double Lt = right - left;
1634 if (symmetric_tran()) Lt *= 2;
1635 double Ll = front - back;
1636 if (symmetric_long()) Ll *= 2;
1637 // double f = 1.;
1638 for (int kt = -nt; kt <= nt; ++kt) {
1639 size_t t = (kt >= 0) ? kt : (symmetric_tran()) ? -kt : kt + nNt;
1640 const double phast = kt * (dest_mesh->at(i).c1 - left) / Lt;
1641 for (int kl = -nl; kl <= nl; ++kl) {
1642 size_t l = (kl >= 0) ? kl : (symmetric_long()) ? -kl : kl + nNl;
1643 res += (*(reinterpret_cast<dcomplex*>(gradients[lay].data() + nNl * t + l) + which) *
1644 exp(2 * PI * I * (kl * (dest_mesh->at(i).c0 - back) / Ll + phast)))
1645 .real();
1646 }
1647 }
1648 if (which && (dest_mesh->at(i).c0 < 0) != (dest_mesh->at(i).c1 < 0)) res = -res;
1649 return res;
1650 });
1651 } else {
1652 size_t nl = symmetric_long() ? nNl : nNl + 1, nt = symmetric_tran() ? nNt : nNt + 1;
1653 DataVector<double> grads;
1654 {
1655 DataVector<dcomplex> work(nl * nt);
1656 for (size_t t = 0; t != nNt; ++t) {
1657 size_t op = nl * t, oc = nNl * t;
1658 for (size_t l = 0; l != nNl; ++l) {
1659 work[op + l] = *(reinterpret_cast<dcomplex*>((gradients[lay].data() + oc + l)) + which);
1660 }
1661 }
1662 auto dct_symmetry = SOLVER->dct2() ? (which ? FFT::SYMMETRY_ODD_2 : FFT::SYMMETRY_EVEN_2)
1663 : (which ? FFT::SYMMETRY_ODD_1 : FFT::SYMMETRY_EVEN_1);
1664 FFT::Backward2D(1, nNl, nNt, symmetric_long() ? dct_symmetry : FFT::SYMMETRY_NONE,
1665 symmetric_tran() ? dct_symmetry : FFT::SYMMETRY_NONE, nl)
1666 .execute(reinterpret_cast<dcomplex*>(work.data()));
1667 grads.reset(nl * nt);
1668 auto dst = grads.begin();
1669 for (const dcomplex& val : work) *(dst++) = val.real();
1670 }
1671 shared_ptr<RegularAxis> lcmesh = plask::make_shared<RegularAxis>(), tcmesh = plask::make_shared<RegularAxis>();
1672 if (symmetric_long()) {
1673 if (SOLVER->dct2()) {
1674 double dx = 0.5 * (front - back) / double(nl);
1675 lcmesh->reset(back + dx, front - dx, nl);
1676 } else {
1677 lcmesh->reset(0., front, nl);
1678 }
1679 } else {
1680 lcmesh->reset(back, front, nl);
1681 for (size_t t = 0, end = nl * nt, dist = nl - 1; t != end; t += nl) grads[dist + t] = grads[t];
1682 }
1683 if (symmetric_tran()) {
1684 if (SOLVER->dct2()) {
1685 double dy = 0.5 * right / double(nt);
1686 tcmesh->reset(dy, right - dy, nt);
1687 } else {
1688 tcmesh->reset(0., right, nt);
1689 }
1690 } else {
1691 tcmesh->reset(left, right, nt);
1692 for (size_t l = 0, last = nl * (nt - 1); l != nl; ++l) grads[last + l] = grads[l];
1693 }
1694 auto src_mesh = plask::make_shared<RectangularMesh<3>>(
1695 lcmesh, tcmesh, plask::make_shared<RegularAxis>(level->vpos(), level->vpos(), 1), RectangularMesh<3>::ORDER_210);
1696
1697 auto sym =
1698 (what == GradientFunctions::COS2) ? InterpolationFlags::Symmetry::POSITIVE : InterpolationFlags::Symmetry::NEGATIVE;
1699
1700 return interpolate(
1701 src_mesh, grads, dest_mesh, interp,
1702 InterpolationFlags(SOLVER->getGeometry(), symmetric_long() ? sym : InterpolationFlags::Symmetry::NO,
1703 symmetric_tran() ? sym : InterpolationFlags::Symmetry::NO, InterpolationFlags::Symmetry::POSITIVE));
1704 }
1705}
1706
1707}}} // namespace plask::optical::modal